Solving KZK numerically (Directional speaker)

Why this matters

The simulation is how you predict what the speaker should sound like before you measure it, and later how you explain what you measured. A solver you have not validated is just a plot generator, so a good part of this step is proving the code right.

The physics

Split each step $h$ in $z$ into pieces that are each easy to solve exactly (Strang splitting: half nonlinear, full linear, half nonlinear):

Diffraction and absorption. In the frequency domain each harmonic and transverse wavenumber $k_x$ evolves independently, so the linear step is just a phase and a decay, $\hat{P} \to \hat{P}\, e^{\,(i k_x^2/2k_n - \alpha_n)\,h}$ . In air $\alpha_n$ comes from ISO 9613-1 for every harmonic separately, so the molecular-relaxation absorption is right at every frequency.
Nonlinearity. $\partial p/\partial z = \frac{\beta}{\rho_0 c_0^3}\, p\,\partial p/\partial \tau$ has the implicit exact solution $p(z{+}h,\tau) = p(z,\ \tau + b\,h\,p(z{+}h,\tau))$ : every point of the waveform slides in $\tau$ in proportion to its pressure.

AM mode. The source is $p_0\,[1 + m\cos(2\pi f_m \tau)]\sin(2\pi f_c\tau)$ with $f_c = 60$ kHz and $f_m = 1$ kHz. The time window is one period of the audio, 1 ms, so “harmonic $n$ ” means $n \times 1$ kHz: the carrier is harmonic 60, and harmonics 1, 2, 3 are the audible tone and its distortion. The code low-passes the on-axis waveform below 20 kHz to get what a listener hears, and reports the total harmonic distortion $\mathrm{THD} = \sqrt{A_2^2 + A_3^2 + \dots}/A_1$ .

Four analytic solutions check the code:

A linear Gaussian beam, whose on-axis amplitude in 2D falls as $(1+z^2/z_R^2)^{-1/4}e^{-\alpha z}$ .
The Fubini solution for a lossless plane wave before shock formation, $B_n = 2J_n(n\sigma)/(n\sigma)$ .
The Fay solution for an old-age sawtooth, $B_n = (2/\Gamma)/\sinh\!\bigl(n(1+\sigma)/\Gamma\bigr)$ .
An AM plane wave in air at short range, from second-order theory. For equal absorption $\alpha$ of the carrier and sidebands it reduces to

p_\text{audio}(z,\tau) = \frac{\beta p_0^2}{4\rho_0 c_0^3}\,\frac{1-e^{-2\alpha z}}{2\alpha}\,\frac{\partial E^2}{\partial \tau}.

In a plane wave the audio follows the first derivative of $E^2$ . Berktay’s second derivative only appears in the far field of a beam.

AM simulation in air at 3 m: the modulated carrier, its spectrum with the audible lines in red, the audible waveform against a pure 1 kHz tone, and the audible spectrum

Audible level of the 1, 2 and 3 kHz lines and the total harmonic distortion against distance on the beam axis

Hints

Hint 1Run the validation before anything else

python kzk_sim.py --validate should print ALL CHECKS PASSED. If you change the solver, run it again. The validation table belongs in your Methods section: it is the evidence that the later plots mean something.

Hint 2Watch it build up with distance

python kzk_gui.py --am opens the viewer in AM mode. Set the carrier, the tone and how far to go, then drag the distance slider to watch the audible tone appear, grow and distort, step by step. The other parameters (carrier amplitude, modulation index, source size, air temperature and humidity) are fixed in the viewer; change them with kzk_sim.py --am options if you need to.

Hint 3Turn the nonlinearity off

Run python kzk_sim.py --am 3 --linear. Where did the audible lines go? That one comparison shows the paper’s main idea: the air, not the speaker, makes the sound.

Hint 4Is that distortion real?

Plain AM with $m = 1$ is expected to distort a lot. Second-order plane-wave theory gives $A_2/A_1 = m/2$ , and Berktay’s far field gives $A_2/A_1 = m$ . Where does your distance sit compared with the carrier’s Rayleigh distance $z_R$ , and does the THD land between the two?

Hint 5Try square-root AM

source_waveform() builds the envelope $1 + m\cos(2\pi f_m\tau)$ . Change it to $\sqrt{1 + m\cos(2\pi f_m\tau)}$ and rerun. How much does the THD drop? Compare with what you measure from the department speaker.

Hint 62D is not the department’s speaker

The code models a strip source, so its beam spreads as $z^{-1/2}$ while a real panel spreads as $z^{-1}$ . Shapes and trends are meaningful; absolute levels need the axisymmetric version of the diffraction step.

What past teams did

2026–27 team

We wrote kzk_sim.py following Lee & Hamilton: FFTs in $x$ and $\tau$ for the linear step, a Newton-solved implicit shift for the nonlinear step, and an absorbing sponge at the transverse edges. Then we added the AM mode (--am), with ISO 9613-1 absorption for every harmonic and figures of the audible waveform, its spectrum and its THD. Two numerical lessons from that:

The nonlinear step originally used linear interpolation, which smooths the waveform a little every step. That is harmless with ~1000 samples per period (the pure tone) but not with the ~64 samples per carrier cycle that the AM window allows, so AM mode uses cubic interpolation.
The audible tone spreads much faster than the 60 kHz beam, so the transverse box has to be sized for the audio, not the carrier. Otherwise the audio reflects off the box edges and corrupts the result after a metre or two.

kzk_gui.py has both modes too: it imports the solver pieces from kzk_sim.py, so the two programs always use the same validated physics.