The speed of sound in a bar: the square root of stiffness over density¶

The speed of a longitudinal wave along a thin solid bar, c = sqrt(E / rho), from its Young's modulus E and its density rho (a phase's, lusine.materials.mechanical_phases' stiffness and density).

It is the thin-rod speed, a bar much thinner than the wavelength, free to swell sideways; in bulk solid the longitudinal wave is faster, sqrt((K + 4/3 G) / rho), by a factor that depends on Poisson's ratio (about 1.15 for steel, whose Poisson ratio is 0.29). Source: Kinsler, Frey, Coppens & Sanders, Fundamentals of Acoustics, 4th ed., 3.2 and appendix A10. Checked against a reference in numpy, which has its own tests.

Source

Kinsler, Frey, Coppens & Sanders, Fundamentals of Acoustics, 4th ed., 3.2 and appendix A10.

Proof

reference, invariant.

speed_of_sound.trait

namespace lusine.properties;

float<m/s> speed_of_sound(float<Pa> stiffness, float<kg/m3> density)
{
  return sqrt(stiffness / density);
}

The speed of a longitudinal wave along a thin solid bar, c = √(E/ρ), from Young's modulus E and the density, both a phase's (lusine.materials.mechanical_phases' stiffness and density). It is the thin-rod speed, not the bulk longitudinal speed √((K + 4G/3)/ρ), which is faster. The library stores no speed of sound: this kernel derives it.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
from lusine.library_sources import parse_rows
EPS = float(np.finfo(np.float32).eps)   # Trait's float is binary32

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from lusine import trait_run from lusine.library_sources import parse_rows EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32

The kernel¶

The law as the library states it.

In [2]:

print(trait_run.source('properties/speed_of_sound.trait'))

print(trait_run.source('properties/speed_of_sound.trait'))

// lusine.properties.speed_of_sound: the speed of a longitudinal wave along a thin solid bar, c = sqrt(E / rho), from
// its Young's modulus E and its density rho (a phase's, lusine.materials.mechanical_phases' stiffness and density). It is
// the thin-rod speed, a bar much thinner than the wavelength, free to swell sideways; in bulk solid the longitudinal
// wave is faster, sqrt((K + 4/3 G) / rho), by a factor that depends on Poisson's ratio (about 1.15 for steel, whose
// Poisson ratio is 0.29).
// Source: Kinsler, Frey, Coppens & Sanders, Fundamentals of Acoustics, 4th ed., 3.2 and appendix A10.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.properties;

float<m/s> speed_of_sound(float<Pa> stiffness, float<kg/m3> density)
{
  return sqrt(stiffness / density);
}

The program¶

Each section runs the kernel over a lattice of items, one a case: the host sets every parameter, and main() binds the kernel over the items.

In [3]:

print(trait_run.law_program('properties/speed_of_sound.trait', 'speed_of_sound'))

print(trait_run.law_program('properties/speed_of_sound.trait', 'speed_of_sound'))

namespace probe;

lusine.version version = { { number = 1 } };

struct item {
  host float<Pa> stiffness;
  host float<kg/m3> density;
  float<m/s> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.properties.speed_of_sound(items.stiffness, items.density); }

Reference: every solid in the library, from its own E and ρ¶

Each phase with a stiffness and a density (each substance has one such phase: the solid's, or a packed sand bed's), read from the library's tables as committed, is held to √(E/ρ) of those two numbers, to binary32's rounding. Steel's 5047.5 m/s and copper's 3809 m/s are within 0.1% of the bar speeds Kinsler et al. tabulate (Fundamentals of Acoustics, 4th ed., appendix A10: steel 5050 m/s, copper 3810 m/s), which are taken from their own E and ρ.

In [4]:

mech = trait_run.source('materials/generated/mechanical.trait')
phases = [r for r in parse_rows(mech, 'mechanical_phases') if 'stiffness' in r and 'density' in r]
stiff = {r['material']: r['stiffness'] for r in phases}
dense = {r['material']: r['density'] for r in phases}
solids = sorted(m for m in stiff if m in dense)
E = np.array([stiff[m] for m in solids])
rho = np.array([dense[m] for m in solids])
c = trait_run.law('properties/speed_of_sound.trait', 'speed_of_sound', stiffness=E, density=rho)
want = np.sqrt(E / rho)
assert np.all(np.abs(c - want) <= 2 * EPS * want), np.abs(c / want - 1).max()
for m, e, r, v in zip(solids, E, rho, c):
    print(f'{m:14s} E {e:9.3g} Pa  rho {r:7.0f} kg/m3  c {v:7.1f} m/s')
kinsler = {'steel': 5050.0, 'copper': 3810.0}
for m, v in kinsler.items():
    got = float(c[solids.index(m)])
    assert abs(got / v - 1) < 1e-3, (m, got)
    print(f"{m}: {got:.1f} m/s against Kinsler et al.'s {v:.0f} m/s ({got / v - 1:+.2%})")
fig, ax = plt.subplots(figsize=(7, 3.2))
ax.bar(solids, c)
ax.set_ylabel('bar speed (m/s)')
plt.setp(ax.get_xticklabels(), rotation=30, ha='right')
plt.show()

mech = trait_run.source('materials/generated/mechanical.trait') phases = [r for r in parse_rows(mech, 'mechanical_phases') if 'stiffness' in r and 'density' in r] stiff = {r['material']: r['stiffness'] for r in phases} dense = {r['material']: r['density'] for r in phases} solids = sorted(m for m in stiff if m in dense) E = np.array([stiff[m] for m in solids]) rho = np.array([dense[m] for m in solids]) c = trait_run.law('properties/speed_of_sound.trait', 'speed_of_sound', stiffness=E, density=rho) want = np.sqrt(E / rho) assert np.all(np.abs(c - want) <= 2 * EPS * want), np.abs(c / want - 1).max() for m, e, r, v in zip(solids, E, rho, c): print(f'{m:14s} E {e:9.3g} Pa rho {r:7.0f} kg/m3 c {v:7.1f} m/s') kinsler = {'steel': 5050.0, 'copper': 3810.0} for m, v in kinsler.items(): got = float(c[solids.index(m)]) assert abs(got / v - 1) < 1e-3, (m, got) print(f"{m}: {got:.1f} m/s against Kinsler et al.'s {v:.0f} m/s ({got / v - 1:+.2%})") fig, ax = plt.subplots(figsize=(7, 3.2)) ax.bar(solids, c) ax.set_ylabel('bar speed (m/s)') plt.setp(ax.get_xticklabels(), rotation=30, ha='right') plt.show()

boxwood        E   1.3e+10 Pa  rho     950 kg/m3  c  3699.2 m/s
copper         E   1.3e+11 Pa  rho    8960 kg/m3  c  3809.1 m/s
iron           E  2.11e+11 Pa  rho    7870 kg/m3  c  5177.9 m/s
silicon        E   4.7e+10 Pa  rho    2330 kg/m3  c  4491.7 m/s
silicon_dioxide E   6.5e+07 Pa  rho    1826 kg/m3  c   188.7 m/s
steel          E     2e+11 Pa  rho    7850 kg/m3  c  5047.5 m/s
tungsten       E  4.11e+11 Pa  rho   19300 kg/m3  c  4614.7 m/s
steel: 5047.5 m/s against Kinsler et al.'s 5050 m/s (-0.05%)
copper: 3809.1 m/s against Kinsler et al.'s 3810 m/s (-0.02%)

Invariant: ρc² is the modulus again¶

Over 64 solids from 10⁸ to 10¹² Pa, ρc² is E to a few roundings; four times as stiff carries the wave exactly twice as fast and four times as dense exactly half as fast, bit for bit (a power of four scales E/ρ without rounding, and the square root is correctly rounded).

In [5]:

rng = np.random.default_rng(320)
E = np.float32(10 ** rng.uniform(8, 12, 64)).astype(float)
rho = np.float32(10 ** rng.uniform(2.5, 4.5, 64)).astype(float)
c = trait_run.law('properties/speed_of_sound.trait', 'speed_of_sound', stiffness=E, density=rho)
assert np.all(np.abs(rho * c * c / E - 1) <= 4 * EPS)
c4 = trait_run.law('properties/speed_of_sound.trait', 'speed_of_sound', stiffness=4 * E, density=rho)
c_4 = trait_run.law('properties/speed_of_sound.trait', 'speed_of_sound', stiffness=E, density=4 * rho)
assert np.array_equal(c4, 2 * c) and np.array_equal(c_4, c / 2)
print('rho c^2 / E - 1: worst', np.abs(rho * c * c / E - 1).max())

rng = np.random.default_rng(320) E = np.float32(10 ** rng.uniform(8, 12, 64)).astype(float) rho = np.float32(10 ** rng.uniform(2.5, 4.5, 64)).astype(float) c = trait_run.law('properties/speed_of_sound.trait', 'speed_of_sound', stiffness=E, density=rho) assert np.all(np.abs(rho * c * c / E - 1) <= 4 * EPS) c4 = trait_run.law('properties/speed_of_sound.trait', 'speed_of_sound', stiffness=4 * E, density=rho) c_4 = trait_run.law('properties/speed_of_sound.trait', 'speed_of_sound', stiffness=E, density=4 * rho) assert np.array_equal(c4, 2 * c) and np.array_equal(c_4, c / 2) print('rho c^2 / E - 1: worst', np.abs(rho * c * c / E - 1).max())

rho c^2 / E - 1: worst 1.062204951773893e-07