I was simulating a structure suspended in air with varying density (o.5, 1, 1.5, 2 gm/cm3) and speed (1.2 power of density) . I was changing value of CFL to make simulations stable when value density value was at 1.5 and 2 as it started to break, without changing it. However, I noticed that I could also change the reference sound speed to do the same thing, to make the simulation more stable. Interestingly, it decreases computational time drastically.

1. Changing Ref sound speed: 1500 m/s, computing time: 9 mins.

2. Changing CFL to 0.04, computing time: 1hr 1 minutes (when CFL is changed, reference speed changes automatically to 577 m/s.)

I calculated transmission loss in both cases. In case 1, I am getting 80 dB with formula used in the code. In case 2, I am getting the same value 0f 80 dB.

Does this mean changing ref speed is more efficient than CFL?

Any insight on this would be much appreciated. Thank you.

Here is my code.

clearvars;

close all;

% create the computational grid

Nx = 3032 ;

Ny = 1496;

dx = 6.6e-06;

dy = dx;

kgrid = kWaveGrid(Nx, dx, Ny, dy);

%load image

img1 = loadImage('1.bmp'); %https://i.imgur.com/hmBLu0q.png

%resize

img1 = imbinarize(resize(img1, [Nx, Ny]));

img_indices1 = find(img1 ==1);

% define the medium properties

%Medium: Air

medium.density = 10* ones(Nx, Ny); % [kg/m^3]

medium.sound_speed = 330*ones(Nx, Ny); % [m/s]

medium.alpha_coeff = 1.64 * ones(Nx,Ny);

medium.alpha_power = 2;

medium.BonA = 20* ones(Nx,Ny);

% create the time array

cfl = 0.05;

t_end = 0.5e-4; % [s]

kgrid.makeTime(max(medium.sound_speed(:)), cfl);

% define source mask for a linear transducer with an odd number of elements

x_offset = 50; % [grid points]

sensor_start = round(Ny *0.1);

sensor_end = round(Ny*0.9);

source.p_mask = makeLine(Nx, Ny, [1 sensor_start], [1 sensor_end]);

% define the properties of the tone burst used to drive the transducer

sampling_freq = 1/kgrid.dt; % [Hz]

tone_burst_freq = 0.5e6; % [Hz]

tone_burst_cycles = 3;

source_strength = 1; %[Pa]

%Appending the source

source.p = source_strength * toneBurst(sampling_freq, tone_burst_freq, tone_burst_cycles);

%Sensor

x_offset1 = round(Nx*0.4); % [grid points]

x_offset2 = round(Nx*0.6);

width = 10; % [grid points]

sensor.mask = zeros(Nx, Ny);

sensor.mask(x_offset1, Ny/2 - width/2 + 1:Ny/2 + width/2) = 1;

sensor.mask(x_offset2, Ny/2 - width/2 + 1:Ny/2 + width/2) = 1;

% define the frequency response of the sensor elements

center_freq = 0.5e6; % [Hz]

bandwidth = 80; % [%]

sensor.frequency_response = [center_freq, bandwidth];

% assign the input options

input_args = { 'PlotPML', false,'PlotLayout', false...

'DataCast', 'gpuArray-single', 'PMLInside',false,...

'RecordMovie', false, 'MovieName', 'example_movie',...

'MovieProfile', 'MPEG-4', 'MovieArgs', {'FrameRate', 10} };

counter = 0;

for i = 0.5:0.5:2

clear sensor_data;

%Defining properties of structure as scaling function

medium.density(img_indices1) = i*100;

a = medium.density(img_indices1);

medium.sound_speed(img_indices1)= a .^ 1.2;

b = mean(medium.sound_speed(img_indices1));

medium.alpha_coeff(img_indices1) = 10;

medium.BonA(img_indices1) = 400;

%reference.sound_speed = 1500;

%run simulations

[sensor_data] = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% compute the amplitude spectra of the recorded time series

[~, sensor_data1] = spect(sensor_data(1,:), 1/kgrid.dt);

[f, sensor_data2] = spect(sensor_data(2,:), 1/kgrid.dt);

% applying gaussian filter to the data

sigma = 3; % pick sigma value for the gaussian

gaussFilter = gausswin(6*sigma + 1)';

gaussFilter = gaussFilter / sum(gaussFilter); % normalize

filteredY1 = conv(sensor_data1, gaussFilter, 'same');

filteredY2 = conv(sensor_data2, gaussFilter, 'same');

max_Y1 = max(filteredY1);

max_Y2 = max(filteredY2);

loss1 = -10* log((max_Y2)/ max_Y1) ;

loss = round(loss1);

disp(loss);

counter = counter +1;

%graphs

subplot(2,2,counter);

plot(f*1e-6, filteredY1, '-k',...

f*1e-6, filteredY2, '-r', 'linewidth', 1);

xlim([0 4]);

xlabel('Frequency (MHz)', 'FontWeight', 'bold');

ylabel('Amplitude (a.u)','FontWeight', 'bold');

legend('Sensor-1', 'Sensor-2', 'Fontsize', 7,'FontWeight', 'bold','FontSize', 9);

title(['Loss = ', num2str(loss), ' dB'], 'FontSize', 9, 'FontWeight', 'bold');

end