Introduction

Acoustic wave simulation is fundamental to many applications, from medical ultrasound imaging to non-destructive testing. Understanding how pressure waves propagate, reflect, and transmit through layered media is essential for designing transducers, interpreting imaging data, and validating theoretical models.

This post explores 1D acoustic wave propagation in a heterogeneous medium—a classic problem that demonstrates key phenomena such as impedance mismatch, partial reflection, and wave speed variation. The simulation is implemented using two complementary tools:

• k-Wave: An open-source MATLAB toolbox using the k-space pseudospectral method
• COMSOL Multiphysics: A commercial finite element software suite

Theoretical Background

The Acoustic Wave Equation

The propagation of acoustic waves in a fluid medium is governed by the linearized equations of fluid mechanics. For a compressible, inviscid fluid, these can be written as coupled first-order partial differential equations:

$$\frac{\partial p}{\partial t} = -\rho_0 c_0^2 \nabla \cdot \mathbf{u}$$

$$\frac{\partial \mathbf{u}}{\partial t} = -\frac{1}{\rho_0} \nabla p$$

where:

• $p$ is the acoustic pressure [Pa]
• $\mathbf{u}$ is the particle velocity [m/s]
• $\rho_0$ is the ambient density [kg/m³]
• $c_0$ is the sound speed [m/s]

In one dimension, these simplify to:

$$\frac{\partial p}{\partial t} = -\rho_0 c_0^2 \frac{\partial u}{\partial x}$$

$$\frac{\partial u}{\partial t} = -\frac{1}{\rho_0} \frac{\partial p}{\partial x}$$

Initial Value Problem (IVP)

An Initial Value Problem specifies the state of the system at $t = 0$ and calculates its evolution over time. In acoustics, this corresponds to scenarios like photoacoustic imaging, where an initial pressure distribution $p_0(x)$ is generated instantaneously (e.g., by laser absorption) and then propagates through the medium.

The initial conditions are: $$p(x, t=0) = p_0(x)$$ $$u(x, t=0) = 0$$

Heterogeneous Medium

In a heterogeneous medium, the material properties $c_0(x)$ and $\rho_0(x)$ vary spatially. This creates:

• Reflection at interfaces due to acoustic impedance mismatch ($Z = \rho_0 c_0$)
• Transmission of the remaining wave energy through the interface
• Speed variation as waves travel faster or slower depending on local properties

The reflection coefficient at a normal interface between two media is:

$$R = \frac{Z_2 - Z_1}{Z_2 + Z_1} = \frac{\rho_2 c_2 - \rho_1 c_1}{\rho_2 c_2 + \rho_1 c_1}$$

Problem Setup

The simulation uses a 1D domain with the following configuration:

Heterogeneous Medium Properties

The medium consists of three regions with different acoustic properties:

Initial Pressure Distribution

A smooth sinusoidal pulse is placed in region 2, starting at x ≈ 1.2 mm with a width of 5 mm (centered at x ≈ 3.6 mm):

$$p_0(x) = \frac{1}{2}\left[\sin\left(\frac{2\pi(x - x_{start})}{w} - \frac{\pi}{2}\right) + 1\right]$$

for $x_{start} \leq x \leq x_{start} + w$, and $p_0 = 0$ elsewhere.

Sensor Locations

Two sensors are placed symmetrically about the domain center:

• Sensor 1: x = -10 mm
• Sensor 2: x = +10 mm

k-Wave Implementation

k-Wave is an open-source MATLAB toolbox for time-domain acoustic and ultrasound simulations. It uses a k-space pseudospectral method that is highly efficient for wave propagation problems.

Simulation Setup

The k-Wave simulation requires defining four main structures:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

% Create computational grid
Nx = 512;           % number of grid points
dx = 0.05e-3;       % grid spacing [m]
kgrid = kWaveGrid(Nx, dx);

% Define heterogeneous medium
% Note: k-Wave grid is centered at x=0, so x ranges from -12.8 to +12.75 mm
medium.sound_speed = 1500 * ones(Nx, 1);    % default [m/s]
medium.sound_speed(1:round(Nx/3)) = 2000;   % x < -4.3 mm (region 1)
medium.density = 1000 * ones(Nx, 1);        % default [kg/m³]
medium.density(round(4*Nx/5):end) = 1500;   % x > +7.7 mm (region 3)

% Create initial pressure (sinusoidal pulse)
x_pos = 280;    % start at grid point 280 → x ≈ 1.2 mm
width = 100;    % 100 grid points → 5 mm width
in = (0:pi/(width/2):2*pi).';
source.p0 = [zeros(x_pos, 1); 
             (sin(in - pi/2) + 1)/2; 
             zeros(Nx - x_pos - width - 1, 1)];

% Define sensor positions
sensor.mask = [-10e-3, 10e-3];  % at x = -10 mm and +10 mm

Simulation Layout

The following figure shows the simulation setup: initial pressure, sensor locations, sound speed profile, and density profile:

Running the Simulation

1
2
3
4
5
6

% Calculate time array
t_end = 2.5 * kgrid.x_size / max(medium.sound_speed(:));
kgrid.makeTime(medium.sound_speed, [], t_end);

% Run simulation
sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor);

Wave Propagation Animation

The following animation shows the acoustic wave propagating through the heterogeneous medium. Notice the reflections at the interfaces and the different propagation speeds in each region:

Key observations:

1. The initial pulse splits into left-going and right-going waves
2. At x ≈ -4.3 mm (interface 1), partial reflection occurs due to impedance mismatch
3. At x ≈ 7.7 mm (interface 2), another reflection occurs
4. The wave travels faster in region 1 (c = 2000 m/s) than in regions 2 and 3 (c = 1500 m/s)

Sensor Results

The pressure recorded at the two sensor positions over time:

COMSOL Implementation

COMSOL Multiphysics uses the Finite Element Method (FEM) to solve the acoustic wave equation. While more computationally expensive, it offers a graphical interface and supports complex geometries.

Model Setup in COMSOL

Create New Model

1. Open COMSOL Multiphysics 6.4
2. Select Model Wizard → 1D geometry
3. Add physics: Acoustics → Pressure Acoustics, Transient
4. Select study: Time Dependent

Define Parameters

Go to Global Definitions → Parameters and add:

Build Geometry

1. Right-click Geometry 1 → Interval
2. Create three intervals corresponding to the three material regions:
   • Region 1: [-12.8 mm, -4.3 mm]
   • Region 2: [-4.3 mm, 7.7 mm]
   • Region 3: [7.7 mm, 12.8 mm]
3. Click Build All

Define Materials

Create three materials with properties matching each region:

1. Right-click Materials → Blank Material
2. For each material, set:
   • Sound speed: c1 or c2 as appropriate
   • Density: rho1 or rho3 as appropriate
3. Assign each material to its corresponding domain

Configure Physics

1. Initial Values:

   • Click on Initial Values 1
   • Set Pressure: Enter the initial pulse expression

   1 2	p0_amp * (sin(2*pi*(x-x_pulse_start)/pulse_width - pi/2) + 1)/2 * (x >= x_pulse_start) * (x <= x_pulse_end)

2. Boundary Conditions:

   • Right-click Pressure Acoustics, Transient → Plane Wave Radiation
   • Apply to both end points to simulate absorbing boundaries (equivalent to PML)

Create Mesh

1. Right-click Mesh 1 → Edge
2. Set Maximum element size to 0.05 mm (matching k-Wave grid)
3. Click Build All

Configure Study

1. Click Study 1 → Step 1: Time Dependent
2. Set Times: range(0, 0.0125[μs], 32[μs])
3. Under Solver Configurations, adjust relative tolerance to 1e-5 for accuracy

Run and Post-Process

1. Click Compute
2. Create 1D Plot Group for sensor data:
   • Add Cut Point 1D datasets at x = -10 mm and x = +10 mm
   • Create Point Graph to plot pressure vs. time

COMSOL Results

The COMSOL animation shows the same wave propagation physics:

The COMSOL simulation produces comparable results to k-Wave:

Results Comparison

Sensor Data Overlay

The following figure directly compares the k-Wave and COMSOL results at both sensor positions:

Both simulations show excellent agreement, with the wave arrivals and reflection patterns matching closely.

Efficiency Comparison

One of the most significant differences between k-Wave and COMSOL is computational efficiency.

Benchmark Results

Speedup Analysis

k-Wave is approximately 155× faster than COMSOL for this 1D acoustic simulation.

This dramatic speedup is due to several factors:

1. Algorithm Efficiency: k-Wave uses a k-space pseudospectral method with FFT-based spatial derivatives, which is inherently efficient for wave propagation problems. COMSOL uses FEM, which requires assembling and solving large matrix systems.
2. Overhead: COMSOL has significant overhead for model management, GUI updates, and general-purpose infrastructure. k-Wave is purpose-built for acoustic simulations.
3. Time Stepping: k-Wave uses an explicit k-space corrected scheme, while COMSOL uses an implicit Generalized-α method. The implicit scheme is more stable but computationally heavier per step.
4. Dimensionality: The speedup advantage of k-Wave increases with problem dimensionality (2D, 3D) due to the superior scaling of spectral methods.

When to Use Each Tool

Conclusion

Both k-Wave and COMSOL successfully simulate 1D acoustic wave propagation in heterogeneous media, producing results that agree well with each other. However, they have different strengths:

• k-Wave excels in computational efficiency, making it ideal for large-scale simulations, parameter sweeps, and inverse problems in medical ultrasound and photoacoustic imaging.
• COMSOL provides greater flexibility for complex geometries, multi-physics coupling, and scenarios where a GUI-based workflow is preferred.

References

1. Treeby, B. E., & Cox, B. T. (2010). k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields. Journal of Biomedical Optics, 15(2), 021314.
2. COMSOL Multiphysics Reference Manual, version 6.4.
3. Kinsler, L. E., Frey, A. R., Coppens, A. B., & Sanders, J. V. (2000). Fundamentals of Acoustics. John Wiley & Sons.