Hi,
I’d like to say I’ve found the k-wave software really useful. Because of COVID I couldn’t access lab equipment to run an experiment and I’m using k-wave to simulate it.
It’s really good work. What I’m using now is the elastic solver in 2d to simulate surface motion from a laser ultrasound setup for NDT.
The code is the following:
close all
clear all
%Physical parameters (all in SI units)
%Air properties
factor = 2.0;
rho_air = 1.225/factor^2;
c_long_air = 343.0*factor;
c_shear_air = 0.0;
%Aluminum properties (pg 273 Scruby book)
rho_Al = 2700.0;
c_long_Al = 6400.0;
c_shear_Al = 3150.0;
%Laser properties
t_pulse_width = 10.0e-9; %(10 ns)
f_max = 0.1874/t_pulse_width
%Simulation domain
x_Al = 2.00e-3;
x_air = 0.56e-3;
y_size_length = 7.68e-3; %Factor of 3 instead of 2
%Generation laser pulse width
%Compute dx and Nx based on a desired x_size and f_max
points_per_wavelength = 3; %Can use min of 3, better to use higher number
max_possible_dx = min([c_long_air c_long_Al c_shear_Al])/f_max/points_per_wavelength
dx = 5.0e-6*factor
%Try to make these numbers a power of 2 (pg 30 of k-Wave manual)
Nx_Al = round(x_Al/dx)
Nx_air = round(x_air/dx)
Nx = Nx_Al+Nx_air
Ny = round(y_size_length/dx)
%Create computational grid
kgrid = kWaveGrid(Nx,dx, Ny,dx);
%Define the compressional sound speed [m/s]
medium.sound_speed_compression = c_long_Al*ones(Nx, Ny);
medium.sound_speed_compression((Nx_Al+1):end,:) = c_long_air;
%Define the shear sound speed [m/s]
medium.sound_speed_shear = c_shear_Al*ones(Nx, Ny);
medium.sound_speed_shear((Nx_Al+1):end,:) = c_shear_air;
%Define the mass density [kg/mˆ3]
medium.density = rho_Al*ones(Nx, Ny);
medium.density((Nx_Al+1):end,:) = rho_air;
%Define the absorption coefficients [dB/(MHzˆ2 cm)]
medium.alpha_coeff_compression = 0.0;
medium.alpha_coeff_shear = 0.0;
%Define circular through hole (along width)
depth = 1.0e-3;
diam = 0.5e-3;
y_position = round(0.5e-3/dx);
x_position = round((x_Al-depth)/dx);
radius = round(diam/2.0/dx);
plot_disc = 0;
disc = makeDisc(Nx,Ny, x_position,y_position, radius, plot_disc);
medium.sound_speed_compression(disc == 1) = c_long_air;
medium.sound_speed_shear(disc == 1) = c_shear_air;
medium.density(disc == 1) = rho_air;
figure; surf(medium.density)
%Create the time array
cfl = 0.008; %(0.02 with air props) % Courant-Friedrichs-Lewy number
t_end = 1.2e-6; % [s]
kgrid.makeTime(max(medium.sound_speed_compression(:)), cfl, t_end);
%Define an initial pressure distribution using makeBall
height_from_bottom = 1.9e-3;
from_side = 5.0e-3;
source.p0 = 2.0*makeDisc(Nx,Ny,round(height_from_bottom/dx),round(from_side/dx), 5); %Radius should be 10
%Define a 3D binary sensor mask for a single grid point
from_side = 2.0e-3;
sensor.mask = zeros(Nx,Ny);
sensor.mask(Nx_Al,round(from_side/dx)) = 1;
sensor.record = {'p','u'};
%Plot sound source and sensor locations
figure; surf(logical(source.p0)+sensor.mask)
%Run the simulation
input_args = {'PlotScale',[-0.75, 0.75, -0.15, 0.15], 'PlotPML',false, 'DisplayMask',zeros(Nx,Ny), 'DataCast','single', ...
'RecordMovie',true, 'MovieName','example_movie', 'PlotScale',[-0.1, 0.1], 'PlotFreq',50, 'MovieProfile','MPEG-4'};
sensor_data = pstdElastic2D(kgrid, medium, source,sensor, input_args{:});
%Plot A-scan data
figure; plot(kgrid.t_array, sensor_data.ux)I have a few questions regarding a more accurate simulation and computation speed.
• Right now I’m simulating an Aluminum sample with a through hole and a free surface (with air properties).
Instead of having an air layer, how could I enforce a free layer (reflecting?) boundary condition instead of a sample air interface. Are the BCs absorbing by default?
• I tried using the gpu option but it’s not working at the moment. Do I have to install something else to make this work?
• Is there an example (papers and/or code) of using a point, line (or half circle?) with a certain directivity function and time profile to more accurately simulate a laser pulse?
Thank you,
Alex