Hi Brad,
Regarding my previous question about simulating plane waves of different frequencies but still having the same time-averaged incident intensities, I wanted to show you some simple code I have been working with that seems to be yielding results that differ from what we would expect. The code is shown below. As a simple test case, I simulated the passage of an acoustic plane wave through a uniform, inviscid medium. I defined two sensor points within the simulation grid, and I computed the mean intensity measured at these sensor points over approximately 2000 acoustic cycles. As you mentioned in your previous response, computing the average acoustic intensity over exactly one acoustic cycle would yield the correct average intensity I am looking for. However, since I did not discretize the time step to ensure sampling of the wave over exactly one cycle, I chose to average over a large number of cycles. This way, the mean intensity computed would have a relatively small error if I happened to sample only a fraction of a complete cycle.
As you mentioned in your previous response, for the speed of sound, mass density, and the source acoustic pressure amplitude I defined in the simulation (1510 m/s, 1020 kg/m^3, and 5 Pa, respectively), I should expect a mean measured acoustic intensity very nearly equal to
I = P^2/(2*c*rho) = ((5 Pa)^2)/(2*1510 m/s*1020 kg/m^3) = 8.12e-6 W/m^2.
However, when I run the simulation, the mean acoustic intensity recorded by both detector elements is 1.532e-6 W/m^2, more than an 80% error with respect to the theoretical mean intensity. I was hoping you could help show me where my mistake is, either in the code shown below or my understanding of your previous responses.
Thanks again for your help!
Sincerely,
Jonathan
*******************************************************************************
function planewave2D
dx = 1.5e-4; % grid point spacing in the x direction [meters]
dy = 1.5e-4; % grid point spacing in the y direction [meters]
Nx = 391; % number of pixels in x-direction
Ny = 333 + 40; % number of pixels in y-direction
kgrid = makeGrid(Nx, dx, Ny, dy); % make computational grid
middle = 196;
% define properties of the background medium
medium.sound_speed = 1510*ones(Nx,Ny);
medium.density = 1020*ones(Nx,Ny);
% define sensor
sensor.mask = zeros(Nx,Ny);
sensor.mask(middle,50) = 1;
sensor.mask(middle,353) = 1;
sensor.record = {'I'};
sensor.record_start_index = ceil(.051/(1510*dt)) + 100;
%define source
source.p_mask = zeros(Nx,Ny);
source.p_mask(:,20) = 1;
source.p = source_mag*sin(2*pi*source_freq.*kgrid.t_array);
source.p_mode = 'dirichlet';
frequency = 4; % transducer frequency [MHz]
source_freq = frequency*(1e6); % [Hz]
% define a time varying sinusoidal source
source_mag = 5;
tmax = (2000/(frequency*1e6)) + (3.4e-5);
[kgrid.t_array, dt] = makeTime(kgrid, 1510, 0.3, tmax); % time sampling intervals
source.p = filterTimeSeries(kgrid, medium, source.p, 'PlotSpectrums', false, 'PPW', 2);
medium.sound_speed = smooth(kgrid, medium.sound_speed, true);
medium.density = smooth(kgrid, medium.density, true);
input_args = {'DisplayMask', (source.p_mask + sensor.mask), 'DataCast', 'single', 'PMLSize', [0 20], 'PMLAlpha', [0 2], 'PMLInside', false};
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
results = sensor_data.Iy;
I0 = results(1,:);
meanI0 = mean(I0)
It = results(2,:);
meanIt = mean(It)
end