Hi Im trying to simulate open baffle speaker box. I based my code on examples with Slit Diffraction (for replacing boundary conditions) and from Transducer Field Patterns Example.
I'm simulating the membrane by putting two sources - one in front of baffle, one on the other side with opposite phase - In comments of the code. The code below is for one side of membrane only.
The problem is, that the source placed inside of the box, creates much bigger pressure field, than the source outside. I'm not sure, if this is caused by standing wave resonance of the open box, or by some computational mistake I missed. But the difference is massive, the sources have completely same initial conditions, on the link below is the image of recorded RMS pressure in Pascals for both cases.
Graph link:
https://ibb.co/G382xJR
THE CODE:
clearvars;
% =========================================================================
% SIMULATION
% =========================================================================
% create the computational grid
Nx = 200; % number of grid points in the x (row) direction
Ny = Nx; % number of grid points in the y (column) direction
dx = 0.025; % grid point spacing in the x direction - vždy 2,5cm [m]
dy = 0.025; % grid point spacing in the y direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy);
t_end=0.5; % set end of calculation [s]
%define properties of test signal
source_freq = 60; % [Hz]
source_mag = 1; % [Pa]
% define the properties of the propagation medium
medium.sound_speed = 344; % [m/s]
%kgrid.makeTime(medium.sound_speed);
%kgrid.t_array = makeTime(kgrid, medium.sound_speed, 0.3, 0.050);
% create the speaker
slit_mask = zeros(Nx, Ny);
% front baffle of speaker
slit_mask(Nx/2,41:58) = 1;
% left baffle of speaker
slit_mask(Nx/2+1:Nx/2+52,41) = 1;
% right baffle of speaker
slit_mask(Nx/2+1:Nx/2+52,58) = 1;
% define a transducer element
source.p_mask = zeros(Nx, Ny);
source.p_mask(Nx/2+1, 43:55) = 1;
% source.p_mask(Nx/2-1, 43:55) = 1; FOR BOTH SIDES OF SOURCE
% assign the slit to the properties of the propagation medium
barrier_scale = 15;
c0=344;
rho0=1.29;
medium.sound_speed = c0 * ones(Nx, Ny);
medium.density = rho0 * ones(Nx, Ny);
medium.sound_speed(slit_mask == 1) = barrier_scale * c0;
medium.density(slit_mask == 1) = 310 * rho0;
% assign the reference sound speed to the background medium
medium.sound_speed_ref = c0;
% create the time array
%kgrid.makeTime(medium.sound_speed);
kgrid.t_array = makeTime(kgrid, medium.sound_speed, 0.3, 0.050);
%define a time varying sinusoidal source
signal_A=[source_mag * sin(2 * pi * source_freq * kgrid.t_array)];
% =========================================================================
% signal for Both sides of source
% signal_A=[source_mag * sin(2 * pi * source_freq * kgrid.t_array); source_mag2 * sin(2 * pi * source_freq * kgrid.t_array+pi);source_mag * sin(2 * pi * source_freq * kgrid.t_array); source_mag2 * sin(2 * pi * source_freq * kgrid.t_array+pi);source_mag * sin(2 * pi * source_freq * kgrid.t_array); source_mag2 * sin(2 * pi * source_freq * kgrid.t_array+pi);source_mag * sin(2 * pi * source_freq * kgrid.t_array); source_mag2 * sin(2 * pi * source_freq * kgrid.t_array+pi);source_mag * sin(2 * pi * source_freq * kgrid.t_array); source_mag2 * sin(2 * pi * source_freq * kgrid.t_array+pi);source_mag * sin(2 * pi * source_freq * kgrid.t_array); source_mag2 * sin(2 * pi * source_freq * kgrid.t_array+pi);source_mag * sin(2 * pi * source_freq * kgrid.t_array); source_mag2 * sin(2 * pi * source_freq * kgrid.t_array+pi);source_mag * sin(2 * pi * source_freq * kgrid.t_array); source_mag2 * sin(2 * pi * source_freq * kgrid.t_array+pi);source_mag * sin(2 * pi * source_freq * kgrid.t_array); source_mag2 * sin(2 * pi * source_freq * kgrid.t_array+pi);source_mag * sin(2 * pi * source_freq * kgrid.t_array); source_mag2 * sin(2 * pi * source_freq * kgrid.t_array+pi);source_mag * sin(2 * pi * source_freq * kgrid.t_array); source_mag2 * sin(2 * pi * source_freq * kgrid.t_array+pi);source_mag * sin(2 * pi * source_freq * kgrid.t_array); source_mag2 * sin(2 * pi * source_freq * kgrid.t_array+pi);source_mag * sin(2 * pi * source_freq * kgrid.t_array); source_mag2 * sin(2 * pi * source_freq * kgrid.t_array+pi)];
% =========================================================================
source.p=signal_A;
% create a display mask to display the transducer
display_mask = source.p_mask;
% create a sensor mask covering the entire computational domain using the
% opposing corners of a rectangle
sensor.mask = [1, 1, Nx, Ny].';
% set the record mode capture the final wave-field and the statistics at
% each sensor point
sensor.record = {'p_final', 'p_max', 'p_rms'};
% assign the input options
input_args = {'DisplayMask',slit_mask, 'PMLInside', false, 'PlotPML', false};
% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
% =========================================================================
% VISUALISATION
% =========================================================================
% add the source mask onto the recorded wave-field
sensor_data.p_final(source.p_mask ~= 0) = 1;
sensor_data.p_max(source.p_mask ~= 0) = 1;
sensor_data.p_rms(source.p_mask ~= 0) = 1;
% plot the final wave-field
figure;
subplot(1, 3, 1);
imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.p_final, [-1 1]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('Final Wave Field');
% plot the maximum recorded pressure
subplot(1, 3, 2);
imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.p_max, [-1 1]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('Maximum Pressure');
% plot the rms recorded pressure
subplot(1, 3, 3);
imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.p_rms, [0 0.5]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('RMS Pressure');
scaleFig(2, 1);