k-Wave User Forum » Topic: Unexpected heating at the focus

k-Wave User Forum » Topic: Unexpected heating at the focus http://www.k-wave.org/forum/topic/unexpected-heating-at-the-focus Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 01:38:14 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Unexpected heating at the focus" http://www.k-wave.org/forum/topic/unexpected-heating-at-the-focus#post-8286 Sat, 07 Aug 2021 08:06:30 +0000 Bradley Treeby 8286@http://www.k-wave.org/forum/ Why do you think it's incorrect? The amount of heating will depend on the volume rate of heat deposition, which depends on the acoustic intensity / pressure. It will also depend on the level of perfusion you set. You could try tuning these parameters. Deepak Sonker on "Unexpected heating at the focus" http://www.k-wave.org/forum/topic/unexpected-heating-at-the-focus#post-8284 Fri, 06 Aug 2021 21:13:13 +0000 Deepak Sonker 8284@http://www.k-wave.org/forum/ 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: <pre><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); %</code></pre>