Hi,
I am getting unexpected heating at focus while using kWaveDiffusion. I have also checked the answer with bioheatExact but it also gave me the same heating. I think it's not correct because I only opened source for 1 sec and it gave me a temperature around 75 deg C.
I don't think this is a legitimate answer and I am lost about how I should correct it.
Code:
clear all
clc
clearvars;
% medium properties
c0 = 1585;
% source parameters
source_f0 = 0.5e6; % source frequency [Hz]
source_roc = 63.2e-3; % bowl radius of curvature [m]
source_diameter = 64e-3; % bowl aperture diameter [m]
source_amp = 0.5e6; % source pressure [Pa]
% grid parameters
axial_size = 130e-3; % total grid size in the axial dimension [m]
lateral_size = 80e-3; % total grid size in the lateral dimension [m]
% computational parameters
ppw = 3; % number of points per wavelength
record_periods = 3; % number of periods to record
cfl = 0.3; % CFL number
source_x_offset = 4; % grid points to offset the source
% =========================================================================
% GRID
% =========================================================================
% calculate the grid spacing based on the PPW and F0
dx = c0 / (ppw * source_f0); % [m]
% compute the size of the grid
Nx = roundEven(lateral_size / dx);
Ny = roundEven(lateral_size / dx);
Nz = roundEven(axial_size / dx);
% create the computational grid
kgrid = kWaveGrid(Nx, dx, Ny, dx, Nz, dx);
% compute points per temporal period
PPP = round(ppw / cfl);
% compute corresponding time spacing
dt = 1 / (PPP * source_f0);
% create the time array using an integer number of points per period
t_end = sqrt( kgrid.x_size.^2 + kgrid.y_size.^2 + kgrid.z_size.^2 )...
/ c0; % total compute time [s] (this must be long enough to reach steady state)
Nt = round(t_end / dt);
kgrid.setTime(Nt, dt);
% calculate the actual CFL and PPW
disp(['PPW = ' num2str(c0 / (dx * source_f0))]);
disp(['CFL = ' num2str(c0 * dt / dx)]);
% =========================================================================
% SOURCE
% =========================================================================
% create time varying source
source_sig = createCWSignals(kgrid.t_array, source_f0, source_amp, 0);
% set bowl position and orientation
bowl_pos = [0, 0, kgrid.z_vec(1) + source_x_offset * kgrid.dx];
focus_pos = [0, 0, 0];
% create empty kWaveArray
karray = kWaveArray('BLITolerance', 0.01, 'UpsamplingRate', 10);
% add bowl shaped element
karray.addBowlElement(bowl_pos, source_roc, source_diameter, focus_pos);
% assign binary mask
source.p_mask = karray.getArrayBinaryMask(kgrid);
% assign source signals
source.p = karray.getDistributedSourceSignal(kgrid, source_sig);
% =========================================================================
% MEDIUM
% =========================================================================
% define the properties of the propagation medium
medium.sound_speed = 1481*ones(Nx,Ny,Nz); % [m/s]
medium.density = 998*ones(Nx,Ny,Nz); % [kg/m^3]
medium.alpha_coeff = 0.025*ones(Nx,Ny,Nz); % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;
for ii = (Nx/2)-20:(Nx/2)+20
for jj = (Ny/2)-20:(Ny/2)+20
for kk = (Nz/2)-20:(Nz/2)+20
% rr = ceil(sqrt((Nx/2 - ii)^2 + (Ny/2 - jj)^2) - 20);
rr = (sqrt((Nx/2 - ii)^2 + (Ny/2 - jj)^2) + ((Nz/2 - kk)^2) - 20);
if(rr<0)
%
medium.sound_speed(ii, jj, kk) = 1585; % [m/s]
medium.density(ii, jj, kk) = 1040; % [kg/m^3]
medium.alpha_coeff(ii, jj, kk) = 0.555;
end
end
end
end
% --------------------
% SENSOR
% --------------------
% set sensor mask to record central plane, not including the source point
sensor.mask = ones(Nx, Ny, Nz);
% sensor.mask(source_x_offset + 2:end, :, Nz/2 + 1) = 1;
% record the pressure
sensor.record = {'p'};
% average only the final few periods when the field is in steady state
sensor.record_start_index = kgrid.Nt - record_periods * PPP + 1;
% =========================================================================
% SIMULATION
% =========================================================================
% set input options
input_args = {...
'PMLSize', 'auto', ...
'PMLInside', false, ...
'PlotPML', false, ...
'DisplayMask', 'off', 'PlotScale', [-1, 1] * source_amp};
% run code
% sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ...
% input_args{:}, ...
% 'DataCast', 'single', ...
% 'PlotScale', [-1, 1] * source_amp);
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
%, ...
%'Dim', 2, 'Window', 'Rectangular', 'FFTPadding', 1);
% =========================================================================
% % CALCULATE HEATING
% =========================================================================
%
% extract amplitude from the sensor data
p = extractAmpPhase(sensor_data.p, 1/kgrid.dt, source_f0);
alpha_np = db2neper(medium.alpha_coeff, medium.alpha_power) * ...
(2 * pi * source_f0).^medium.alpha_power;
% % reshape data
p = reshape(p, Nx, Ny, Nz);
Q = alpha_np .* p.^2 ./ (medium.density .* medium.sound_speed);
% =========================================================================
% THERMAL SIMULATION
% =========================================================================
% % set the background temperature and heating term
source.Q = Q;
source.T0 = 37*ones(Nx,Ny,Nz);
% % define medium properties related to diffusion
medium.density = 1040; % [kg/m^3]
medium.thermal_conductivity = 0.4683; % [W/(m.K)]
medium.specific_heat = 3210; % [J/(kg.K)]
% % create kWaveDiffusion object
kdiff = kWaveDiffusion(kgrid, medium, source, []);
% % set source on time and off time
on_time = 1; % [s]
% off_time = 20; % [s]
% % set time step size
dt = 0.1;
% % take time steps
kdiff.takeTimeStep(round(on_time / dt), dt);
% % store the current temperature field
T1 = kdiff.T;
% % % turn off heat source and take time steps
% kdiff.Q = 0;
% kdiff.takeTimeStep(round(off_time / dt), dt);
% %
% % % store the current temperature field
% T2 = kdiff.T;
% =========================================================================
% =========================================================================
% THERMAL SIMULATION Validation by Pennes Bioheat Transfer equation
% =========================================================================
D = medium.thermal_conductivity / (medium.density * medium.specific_heat);
S = Q/(medium.density * medium.specific_heat);
T_exact = bioheatExact(source.T0, S, [D, 0, 0], kgrid.dx, on_time);
% =========================================================================
% % plot the acoustic pressure
subplot(3, 3, 1);
imagesc(kgrid.z_vec * 1e3, kgrid.x_vec * 1e3, squeeze(p(:, (kgrid.Ny)/2,:)*1e-6));
h = colorbar;
xlabel(h, '[MPa]');
ylabel('x-position [mm]');
xlabel('z-position [mm]');
axis image;
title('Acoustic Pressure Amplitude in Axial');
% % set colormap and enlarge figure window
colormap(hot(256));
scaleFig(1.5, 1);
% plot the acoustic pressure Transverse direction
subplot(3, 3, 2);
imagesc(kgrid.x_vec * 1e3, kgrid.y_vec * 1e3, squeeze(p(:, :,(kgrid.Nz)/2)*1e-6));
h = colorbar;
xlabel(h, '[MPa]');
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('Acoustic Pressure Amplitude in Radial');
% % set colormap and enlarge figure window
colormap(hot(256));
scaleFig(1.5, 1);
% % plot the volume rate of heat deposition
subplot(3, 3, 3);
imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, squeeze(Q(:, :, (kgrid.Nz)/2)*1e-9));
h = colorbar;
xlabel(h, '[kW/cm^3]');
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('Volume Rate Of Heat Deposition');
% % set colormap and enlarge figure window
colormap(hot(256));
scaleFig(1.5, 1);
% % plot the temperature after heating
subplot(3, 3, 4);
imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, squeeze(T1(:, (kgrid.Ny)/2, :)));
h = colorbar;
xlabel(h, '[degC]');
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('Temperature After Heating');
% % set colormap and enlarge figure window
colormap(hot(256));
scaleFig(1.5, 1);
% % % %verification of temperature by Pennes Bioheat Transfer equation
subplot(3, 3, 5);
imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, squeeze(T_exact(:, (kgrid.Ny)/2, :)));
h = colorbar;
xlabel(h, '[degC]');
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('Temperature After Heating By Pennes Bioheat Eq.');
% % set colormap and enlarge figure window
colormap(hot(256));
scaleFig(1.5, 1);
% %
% % plot thermal dose
subplot(3, 3, 6);
imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, squeeze(kdiff.cem43((kgrid.Nx)/2, :, :)), [0, 1000]);
h = colorbar;
xlabel(h, '[CEM43]');
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('Thermal Dose');
colormap(hot(256));
scaleFig(1.5, 1);
% % plot lesion map
subplot(3, 3, 7);
imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, squeeze(kdiff.lesion_map(:, (kgrid.Ny)/2, :)), [0, 1]);
colorbar;
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('Ablated Tissue');
% % set colormap and enlarge figure window
colormap(hot(256));
scaleFig(1.5, 1);
%