Hi!
I am currently working on my master thesis and want to simulate the propagation of a pressure distribution in a structure of alterning Al and Pb foils (a few um thick, circular, 2cm radius, the whole thing is in vacuum). I am having problems with my result, as it is diverging: the mean pressure on a thin detector on the other side of the stack is going towards 10^293Pa. I assumed it accured because my geometry is not smooth enougth, so to check I just built a block of Pb in vacuum and put a circular pressure source inside of it, and run the simulation in 2D. As I understood from one of the last posts (point 4) it shouldn't be a problem if the source originates from within the high impedance material. Why does my code still diverge? is there a solution to this problem? i tried to change the material properties, but I have to change them so much, that the simulation is not rapresentative anymore.
Here is my code:
%% % create masks
z_lim = 0.03;
x_lim = 0.01;%in m
Nx = 200; % number of grid points in the x direction
Nz = 600; % number of grid points in the z direction
dz = 0.00005;
dx = 0.0001;
[Z,X] = meshgrid((0:dz:z_lim-dz),(-1*x_lim:dx:x_lim-dx));
geom = zeros(size(Z));
det = zeros(size(Z));
lz = 0.02; % [m], z length
lx = 0.01; % [m], x length
z0 = 0.015; % [m], z center
x0 = 0; % [m], x center
ind = ((Z >= z0-lz/2) & (Z < z0+lz/2) & (X >= x0-lx/2) & (X < x0+lx/2));
[num, ~] = materials(4); %add the material of interest
geom(ind) = num;
geom = permute(geom,[2,1]);
% Create the computational grid for the 2D slice
kgrid = kWaveGrid(Nz, dz, Nx, dx);
[densitymask,speedmask,alphamask,ymask] = get_masks(geom);
%%
% Define the properties of the upper layer of the propagation medium
medium.sound_speed = speedmask;
medium.density = densitymask;
% Create time array
t_end = 3e-6; % [s]
kgrid.makeTime(medium.sound_speed, [], t_end);
%%
%create SOURCE
% define a single source point
disc_magnitude = 300; % [Pa]
disc_x_pos = Nx/2; % [grid points]
disc_z_pos = Nz/4;
disc_radius = 0.4; % [grid points]
disc = disc_magnitude * makeDisc(Nz, Nx, disc_z_pos, disc_x_pos, disc_radius);
source.p0 = zeros(Nz, Nx);
source.p0 = disc ;
% %%
% %SENSOR
det = makeLine(Nz, Nx, [round(0.023/dz), 90],[round(0.023/dz), 110]);
%det = permute(det,[2,1]);
sensor.mask=det;
%%
% Run the simulation with PlotLayout and PlotScale
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,...
'PlotLayout', true);
%%
% Assuming your pressure data is stored in a variable named 'pressure_data'
figure
% Calculate the mean pressure across all pixels at each time step
mean_pressure = mean(sensor_data, 1);
% Plot pressure versus time
plot(kgrid.t_array, mean_pressure, 'LineWidth', 2);
xlabel('Time');
ylabel('Mean Pressure on detector');
title('Mean Pressure vs Time');
grid on;
the materials are defined as follows:
function [num, mat, rho, c, gamma, alpha_coeff,y] = materials(arg)
% attenuation data: Appl. Sci. 2020, 10(7), 2230; https://doi.org/10.3390/app10072230
if arg == 0 | strcmpi(arg,'vacuum') | strcmpi(arg,'vac')
mat = 'vacuum';
num = 0;
rho = 0.0012 ;% 1.2e-3; % [g/cm^3] density
c = 0.4; %0.4; % [mm/us] speed of sound
beta = 3400e-6; % [1/K] vol. thermal expansion
c_p = 1012; % [J/kg/K] spec. heat capacity
gamma = beta*(c*1000)^2/c_p; % [Pa/(J/m^3)] Grüneisenparameter
alpha_coeff = 1e-6; %0 [dB/cm/MHz^y] acoustic attenuation coef.
y = 1; % attenuation exponent
elseif arg == 4 | strcmpi(arg,'lead') | strcmpi(arg,'Pb')
mat = 'Pb';
num = 4;
A = 207; % [u] eff. Massezahl
Z = 82; % [u] eff. Ordnungszahl
rho = 11.29;%11.29; % [g/cm^3] density
c = 1.96; % 1.96 [mm/us] speed of sound
beta = 29e-6; % 87e-6; % [1/K] vol. thermal expansion
c_p = 129; % [J/kg/K] spec. heat capacity
gamma = beta*(c*1000)^2/c_p; % [Pa/ (J/m^3)] Grüneisenparameter
I = 823; % [eV] mean excitation energy
alpha_coeff = 4.33; %alpha = 4.33; % [dB/cm/MHz] acoustic attenuation coef.
y = 1; % attenuation exponent
end
and get masks is a function:
function [densitymask,speedmask,alphamask,ymask] = get_masks(geom)
% create mask arrays
densitymask = zeros(size(geom));
speedmask = zeros(size(geom));
alphamask = zeros(size(geom));
ymask = zeros(size(geom));
% initialise counter
n = 0;
while ~all(densitymask(:))
% assign values of materials to masks
[~, ~, rho, c, ~, alpha_coeff,y] = materials(n);
densitymask(geom == n) = rho*1e3; % [kg/m^3]
speedmask(geom == n) = c*1e3; % [m/s]
alphamask(geom == n) = alpha_coeff; % [dB/cm/MHz^y]
ymask(geom == n) = y; %
% increase counter
n = n+1;
end
end
Thank you in advance for any help.
Best regards,
Gianna