k-Wave
A MATLAB toolbox for the time-domain
simulation of acoustic wave fields
- Getting Started
- Examples
- Initial Value Problems
- Example: Homogenous Propagation Medium
- Example: Using A Binary Sensor Mask
- Example: Defining A Sensor Mask By Opposing Corners
- Example: Loading External Image Maps
- Example: Heterogeneous Propagation Medium
- Example: Saving Movie Files
- Example: Recording The Particle Velocity
- Example: Defining A Gaussian Sensor Frequency Response
- Example: Comparison Of Modelling Functions
- Example: Setting An Initial Pressure Gradient
- Example: Simulations In One Dimension
- Example: Simulations In Three Dimensions
- Example: Photoacoustic Waveforms in 1D, 2D and 3D
- Time Varying Source Problems
- Example: Monopole Point Source In A Homogeneous Propagation Medium
- Example: Dipole Point Source In A Homogeneous Propagation Medium
- Example: Simulating Transducer Field Patterns
- Example: Steering A Linear Array
- Example: Snell's Law And Critical Angle Reflection
- Example: The Doppler Effect
- Example: Diffraction Through A Slit
- Example: Simulations In Three-Dimensions
- Sensor Directivity
- Example: Focussed Detector in 2D
- Example: Focussed Detector in 3D
- Example: Modelling Sensor Directivity in 2D
- Example: Modelling Sensor Directivity in 3D
- Example: Sensor Element Directivity in 2D
- Example: Focussed 2D Array with Directional Elements
- Photoacoustic Image Reconstruction
- Example: 2D FFT Reconstruction For A Line Sensor
- Example: 3D FFT Reconstruction For A Planar Sensor
- Example: 2D Time Reversal For A Line Sensor
- Example: 2D Time Reversal For A Circular Sensor
- Example: 3D Time Reversal For A Planar Sensor
- Example: 3D Time Reversal For A Spherical Sensor
- Example: Image Reconstruction With Directional Sensors
- Example: Image Reconstruction With Bandlimited Sensors
- Example: Iterative Image Improvement Using Time Reversal
- Example: Attenuation Compensation Using Time Reversal
- Example: Attenuation Compensation Using Time Variant Filtering
- Example: Automatic Sound Speed Selection
- Diagnostic Ultrasound Simulation
- Example: Defining An Ultrasound Transducer
- Example: Simulating Ultrasound Beam Patterns
- Example: Using An Ultrasound Transducer As A Sensor
- Example: Simulating B-mode Ultrasound Images
- Example: Simulating B-mode Images Using A Phased Array
- Numerical Analysis
- Example: Controlling The Absorbing Boundary Layer
- Example: Source Smoothing
- Example: Filtering A Delta Function Input Signal
- Example: Modelling Power Law Absorption
- Example: Modelling Nonlinear Wave Propagation
- Example: Optimising k-Wave Performance
- Using The C++ Code
- Elastic Wave Propagation
- Example: Explosive Source In A Layered Medium
- Example: Plane Wave Absorption
- Example: Shear Waves And Critical Angle Reflection
- Example: Simulations In Three Dimensions
- Functions - By Category
- Functions - Alphabetical List
- Release Notes
- License
k-Wave Toolbox |
Explosive Source In A Layered Medium Example
On this page… |
---|
Overview
This example provides a simple demonstration of using k-Wave for the simulation and detection of compressional and shear waves in elastic and viscoelastic media within a two-dimensional heterogeneous medium. It builds on the Homogenous Propagation Medium and Heterogeneous Propagation Medium examples.
Defining the medium and source properties
In addition to the simulation functions for modelling compressional waves in fluid media (kspaceFirstOrder1D
, kspaceFirstOrder2D
, and kspaceFirstOrder3D
), k-Wave also includes functions for simulating the propagation of compressional and shear waves in isotropic elastic and viscoelastic media. These functions are based on the pseudospectral time domain (PSTD) method, and are called pstdElastic2D
, pstdElastic3D
. These are used in a very similar fashion to the fluid codes, and require four input structures which define the properties of the computational grid, the material properties of the medium, the properties and locations of any acoustic sources, and the properties and locations of the sensor points used to record the evolution of the wavefield over time. The kgrid
and sensor
inputs are defined in an identical fashion to the fluid code, while the medium
and source
input structures have slightly different field names to reflect the nature of solid materials.
In an isotropic elastic medium, the material properties can be characterised by the shear and compressional sound speeds, and the mass density. These are assigned to the medium
structure as the fields sound_speed_compression
, sound_speed_shear
, and density
. In a homogeneous medium, the material parameters are set as scalar values in SI units. In a heterogeneous medium, these are instead given as matrices with the same dimensions as the computational grid. In this example, a heterogeneous medium is created with two material layers.
% define the properties of the upper layer of the propagation medium medium.sound_speed_compression = 1500*ones(Nx, Ny); % [m/s] medium.sound_speed_shear = zeros(Nx, Ny); % [m/s] medium.density = 1000*ones(Nx, Ny); % [kg/m^3] % define the properties of the lower layer of the propagation medium medium.sound_speed_compression(Nx/2:end, :) = 2000; % [m/s] medium.sound_speed_shear(Nx/2:end, :) = 800; % [m/s] medium.density(Nx/2:end, :) = 1200; % [kg/m^3]
For viscoelastic media, the absorption coefficient in units of dB/(MHz^2 cm) can also be defined for both compressional and shear waves. The elastic codes are based on the classical Kelvin-Voigt absorption model, which gives absorption proportional to frequency squared in the low frequency limit.
% define the absorption properties medium.alpha_coeff_compression = 0.1; % [dB/(MHz^2 cm)] medium.alpha_coeff_shear = 0.5; % [dB/(MHz^2 cm)]
The elastic codes support three types of sources: (1) an initial pressure distribution, (2) time varying velocity sources, and (3) time varying stress sources. These are used in an analogous fashion to the fluid codes. In this example, an explosive source is assigned to source.p0
. Within the simulation function, this is then assigned to the normal components of the stress.
% create initial pressure distribution using makeDisc disc_magnitude = 5; % [Pa] disc_x_pos = 30; % [grid points] disc_y_pos = 64; % [grid points] disc_radius = 5; % [grid points] source.p0 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);
Running the simulation
The elastic codes are based on the pseudospectral time domain method. This does not use a k-space corrected scheme for time integration like the fluid codes, which means smaller time steps are typically required for stable and accurate simulations. In this example, the Courant-Friedrichs-Lewy (CFL) number is set to 0.1.
% create the time array cfl = 0.1; % Courant-Friedrichs-Lewy number t_end = 8e-6; % [s] kgrid.t_array = makeTime(kgrid, max(medium.sound_speed_compression(:)), cfl, t_end);
The simulation is then started by calling pstdElastic2D
with the four input structures analogous to kspaceFirstOrder2D
. Optional inputs can similarly be defined as 'string', value
pairs after the main inputs. By default, as the simulation runs, a visualisation of the propagating wave field and a status bar are displayed with frame updates every ten time steps. Both the normal and shear components of the stress are displayed. The plot scales for these can be defined individually using the 'PlotScale'
input. This is defined as [sii_min, sii_max, sij_min, sij_max]
, where sii
and sij
denote the plot scales for the normal and shear stress.
% define input arguments input_args = {'PlotScale', [-0.75, 0.75, -0.15, 0.15], 'PlotPML', false, 'DisplayMask', display_mask}; % run the simulation sensor_data = pstdElastic2D(kgrid, medium, source, sensor, input_args{:});
A visualisation of the running simulation and the recorded sensor data are given below. Both compressional and shear waves in the lower layer of the medium are clearly visible.
Elastic Wave Propagation | Plane Wave Absorption |
© 2009-2014 Bradley Treeby and Ben Cox.