Hi Brad,
You probably recall that you were helping me some time ago to simulate monochromatic plane waves propagating through simple phantoms. I would like to set my time-averaged incident intensity to be 1 W/cm^2 (1e4 W/m^2) in my current simulation. I am including the code for your reference. When I set the incident intensity to be very low (1e-4 W/m^2), the simulation visually appears to run correctly. However, if I set the incident intensity to be the desired 1e4 W/m^2, the simulation does not seem to run correctly. In particular, the movie of the 2D simulation appears to show that the pressure field propagates the entire length of the phantom in one time step, and drastic changes in pressure can be observed. You only have to look at the first few time steps to see this. Could you please run the attached code and let me know what might be causing the apparent instability in the higher-intensity simulation?
Thanks for your help.
Sincerely,
Jon
function planewave(distance)
% distance defines lesion-to-detector distance (cm)
dx = 1e-4; % grid point spacing in the x direction [meters]
dy = 1e-4; % grid point spacing in the y direction [meters]
Nx = 587; % number of pixels in x-direction
Ny = 500; % number of pixels in y-direction
kgrid = makeGrid(Nx, dx, Ny, dy); % make computational grid
middle = 294;
x0 = middle;
frequency = 3.81;
% define properties of the lesion
y0 = round(500 - (distance*(1/100)*(1/dy)));
medium.sound_speed = 1550*makeDisc(Nx,Ny,x0,y0,6); % [m/s]
medium.density = 1040*makeDisc(Nx,Ny,x0,y0,6); % [kg/m^3]
medium.alpha_coeff = .570*(frequency^1.3)*makeDisc(Nx,Ny,x0,y0,6); % [dB/(MHz^y*cm)]
medium.alpha_power = 0; % value of power y
lesionmask = makeDisc(Nx,Ny,x0,y0,6);
medium.sound_speed(medium.sound_speed==0) = 1550;
medium.density(medium.density==0) = 1015;
medium.alpha_coeff(medium.alpha_coeff==0) = ((0.25*0.158*(frequency^1.7))...
+ (0.75*0.870*(frequency^1.5)));
medium.alpha_power = 0; % value of power y
% define sensor
sensor.mask = zeros(Nx,Ny);
sensor.mask(:,500) = 1;
sensor.record = {'p'};
%define source
source.p_mask = zeros(Nx,Ny);
source.p_mask(:,1) = 1;
source_intensity = 1e-4; % [W/m^2]
rho0 = 1015; % [kg/m^3]
c0 = 1550; % [m/s]
source_mag = sqrt(2*source_intensity*rho0*c0); % [Pa]
source_freq = frequency*(1e6); % [Hz]
% calculate the time length of one acoustic period
T = 1/source_freq;
% create the time array using a nice number of points
points_per_period = 10;
dt = T/points_per_period;
t_end = (10/source_freq) + (3.4e-5);
kgrid.t_array = 0:dt:t_end;
% define a time varying sinusoidal source
source.p = source_mag*sin(2*pi*source_freq.*kgrid.t_array);
% smooth startup using a cosine ramp
ramp_length = 40;
ramp = (-cos((0:(ramp_length-1))*pi/ramp_length ) + 1)/2;
source.p(1:length(ramp)) = source.p(1:length(ramp)) .* ramp;
source.p = filterTimeSeries(kgrid, medium, source.p, 'PlotSpectrums', false, 'PPW', 2);
num_cycles = 1;
sensor.record_start_index = kgrid.Nt - round(num_cycles * T / dt) + 1;
medium.sound_speed = smooth(kgrid, medium.sound_speed, true);
medium.density = smooth(kgrid, medium.density, true);
input_args = {'DisplayMask', (source.p_mask + sensor.mask + lesionmask), 'DataCast', 'single', 'PMLSize', [0, 20], 'PMLAlpha', [0, 2], 'PMLInside', false};
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
pressure = sensor_data.p;
pressure = pressure.^2;
pressure = transpose(pressure);
pressure = mean(pressure);
Iavg = pressure./(rho0*c0);
xaxis = -293*0.1:0.1:293*0.1;
plot(xaxis, Iavg)
xlabel('Lateral Position (mm)')
ylabel('Acoustic Intensity (W/m^2)')
hold on
end