k-Wave User Forum » Topic: 3D simulations of Lamb wave in plate

k-Wave User Forum » Topic: 3D simulations of Lamb wave in plate http://www.k-wave.org/forum/topic/3d-simulations-of-lamb-wave-in-plate Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 00:45:02 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "3D simulations of Lamb wave in plate" http://www.k-wave.org/forum/topic/3d-simulations-of-lamb-wave-in-plate#post-6897 Sat, 22 Jun 2019 20:32:14 +0000 Bradley Treeby 6897@http://www.k-wave.org/forum/ Sounds like a classic instability. Have you tried making the time step (CFL) smaller? Note, your simulations will run much faster if you choose grid sizes with small prime factors (59 is a prime number, so a particularly bad choice). Brad. pakulam on "3D simulations of Lamb wave in plate" http://www.k-wave.org/forum/topic/3d-simulations-of-lamb-wave-in-plate#post-6858 Wed, 22 May 2019 07:39:15 +0000 pakulam 6858@http://www.k-wave.org/forum/ I have the problem with 3D simulation of Lamb wave in plate wich is embended in artificial fluid (which should simulate air). My simulations were based on existing in help code for simulations of shearwave and Snells law. I have performed simulations in 2D and the results are very good and in agreement with theoretical (analitical) model. Then I tried to perform such anaysis in 3D and in fact simulations is working but the results are bad (I receive NAN in the values of UX, uy, uz, p). Below please find well working 2D code and the 3D code I have developed based on 2D code. I would appreciate for any sugestions Kind Regards Michal 1. 2D code - GOOD clear all; clc close all % ========================================================================= % SIMULATION PARAMETERS % ========================================================================= scale = 1; % create the computational grid PML_size = 20; % size of the PML in grid points Nx = 59*scale - 2*PML_size; % number of grid points in the x direction Ny = 300*scale - 2*PML_size; % number of grid points in the y direction % Nz = 300*scale - 2*PML_size; % number of grid points in the y direction dx = 0.5e-3/scale; % grid point spacing in the x direction [m] dy = 0.5e-3/scale; % grid point spacing in the y direction [m] % dz=dx; kgrid = makeGrid(Nx, dx, Ny, dy); % kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz); % define the medium properties cp1 = 600; % compressional wave speed [m/s] cs1 = 0; % shear wave speed [m/s] rho1 = 200; % density [kg/m^3] alpha0_p1 = 5.0; % compressional wave absorption [dB/(MHz^2 cm)] alpha0_s1 = 5.0; % shear wave absorption [dB/(MHz^2 cm)] cp2 = 2656; % compressional wave speed [m/s] cs2 = 1078; % shear wave speed [m/s] rho2 = 1140; % density [kg/m^3] alpha0_p2 = 0.5; % compressional wave absorption [dB/(MHz^2 cm)] alpha0_s2 = 0.5; % shear wave absorption [dB/(MHz^2 cm)] % create the time array cfl = 0.1; t_end = 200e-6; kgrid.t_array= makeTime(kgrid, cp1, cfl, t_end); fs=1/kgrid.dt; % define position of heterogeneous slab % slab = ones(Nx, Ny, Nz); % slab (1,:,:)=0 ; % slab (19,:,:)=0 ; slab = ones(Nx, Ny); slab (1,:)=0 ; slab (19,:)=0 ; % define the source properties source_freq = 0.05e6; % [Hz] source_strength = 2e6; % [Pa] source_cycles = 2; % number of tone burst cycles %3D % source_mask=zeros(Nx,Ny,Nz); % source_mask(2,30,6)=1; % source_mask(18,30,6)=1; %2D source_mask=zeros(Nx,Ny); source_mask(2,6)=1; source_mask(18,6)=1; % assign the source source.s_mask = source_mask; source.sxx(1,:) = source_strength*toneBurst(1/kgrid.dt, source_freq, source_cycles); source.sxx(2,:) = -source_strength*toneBurst(1/kgrid.dt, source_freq, source_cycles); % define the sensor mask sensor.mask = zeros(Nx, Ny); for k=171:1:221 sensor.mask(2,k:k+10) = 1; %2D end sensor.record = {'p','u'}; % set the input arguments input_args = {'PMLSize', PML_size, 'PMLAlpha', 2, 'PlotPML', false, ... 'PMLInside', false, 'PlotScale', [-1, 1]*source_strength/2, ... 'DisplayMask', 'off', 'DataCast', 'single'}; % ========================================================================= % ELASTIC SIMULATION % ========================================================================= % define the medium properties medium.sound_speed_compression = cp1*ones(Nx, Ny); medium.sound_speed_compression(slab == 1) = cp2; medium.sound_speed_shear = cs1*ones(Nx, Ny); medium.sound_speed_shear(slab == 1) = cs2; medium.density = rho1*ones(Nx, Ny); medium.density(slab == 1) = rho2; medium.alpha_coeff_compression = alpha0_p1*ones(Nx, Ny); medium.alpha_coeff_compression(slab == 1) = alpha0_p2; medium.alpha_coeff_shear = alpha0_s1*ones(Nx, Ny); medium.alpha_coeff_shear(slab == 1) = alpha0_s2; % run the elastic simulation sensor_data_elastic = pstdElastic2D(kgrid, medium, source, sensor, input_args{:}); p=sensor_data_elastic.p; subplot(2,2,1),imagesc(p),title('Preassure') ux=sensor_data_elastic.ux; subplot(2,2,2),imagesc(ux),title('Ux') uy=sensor_data_elastic.uy; subplot(2,2,3),imagesc(uy),title('Uy') save('SIM_Results2D.mat','p','ux','uy','fs','source_freq') 2. 3D code - BAD clear all; clc close all % ========================================================================= % SIMULATION PARAMETERS % ========================================================================= scale = 1; % create the computational grid PML_size = 20; % size of the PML in grid points Nx = 59*scale - 2*PML_size; % number of grid points in the x direction Ny = 100*scale - 2*PML_size; % number of grid points in the y direction Nz = 300*scale - 2*PML_size; % number of grid points in the y direction dx = 0.5e-3/scale; % grid point spacing in the x direction [m] dy = 0.5e-3/scale; % grid point spacing in the y direction [m] dz=dx; kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz); % define the medium properties cp1 = 600; % compressional wave speed [m/s] cs1 = 0; % shear wave speed [m/s] rho1 = 200; % density [kg/m^3] alpha0_p1 = 5.0; % compressional wave absorption [dB/(MHz^2 cm)] alpha0_s1 = 5.0; % shear wave absorption [dB/(MHz^2 cm)] cp2 = 2656; % compressional wave speed [m/s] cs2 = 1078; % shear wave speed [m/s] rho2 = 1140; % density [kg/m^3] alpha0_p2 = 0.5; % compressional wave absorption [dB/(MHz^2 cm)] alpha0_s2 = 0.5; % shear wave absorption [dB/(MHz^2 cm)] % create the time array cfl = 0.1; t_end = 100e-6; kgrid.t_array= makeTime(kgrid, cp1, cfl, t_end); fs=1/kgrid.dt; % define position of heterogeneous slab slab = ones(Nx, Ny, Nz); slab (1,:,:)=0 ; slab (19,:,:)=0 ; % define the source properties source_freq = 0.05e6; % [Hz] source_strength = 2e6; % [Pa] source_cycles = 2; % number of tone burst cycles source_mask=zeros(Nx,Ny,Nz); source_mask(2,30,6)=1; source_mask(18,30,6)=1; % assign the source source.s_mask = source_mask; source.sxx(1,:) = source_strength*toneBurst(1/kgrid.dt, source_freq, source_cycles); source.sxx(2,:) = -source_strength*toneBurst(1/kgrid.dt, source_freq, source_cycles); % define the sensor mask sensor.mask = zeros(Nx, Ny, Nz); for k=171:1:221 sensor.mask(2,30,k:k+10) = 1; end sensor.record = {'p','u'}; % set the input arguments input_args = {'PMLSize', PML_size, 'PMLAlpha', 2, 'PlotPML', false, ... 'PMLInside', false, 'PlotScale', [-1, 1]*source_strength/2, ... 'DisplayMask', 'off', 'DataCast', 'single'}; % ========================================================================= % ELASTIC SIMULATION % ========================================================================= % define the medium properties clear medium medium.sound_speed_compression = cp1*ones(Nx, Ny, Nz); medium.sound_speed_compression(slab == 1) = cp2; medium.sound_speed_shear = cs1*ones(Nx, Ny, Nz); medium.sound_speed_shear(slab == 1) = cs2; medium.density = rho1*ones(Nx, Ny, Nz); medium.density(slab == 1) = rho2; medium.alpha_coeff_compression = alpha0_p1*ones(Nx, Ny, Nz); medium.alpha_coeff_compression(slab == 1) = alpha0_p2; medium.alpha_coeff_shear = alpha0_s1*ones(Nx, Ny, Nz); medium.alpha_coeff_shear(slab == 1) = alpha0_s2; % run the elastic simulation sensor_data_elastic = pstdElastic3D(kgrid, medium, source, sensor, input_args{:}); p=sensor_data_elastic.p; subplot(2,2,1),imagesc(p),title('Preassure') ux=sensor_data_elastic.ux; subplot(2,2,2),imagesc(ux),title('Ux') uy=sensor_data_elastic.uy; subplot(2,2,3),imagesc(uy),title('Uy') save('SIM_Results3D.mat','p','ux','uy','fs','source_freq')