k-Wave User Forum » Forum: Thermal simulations

k-Wave User Forum » Forum: Thermal simulations - Recent Topics http://www.k-wave.org/forum/forum/thermal-simulations-1 Support for the k-Wave MATLAB toolbox en-US Mon, 25 May 2026 04:57:46 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php s148275 on "Axisymmetric diffusion" http://www.k-wave.org/forum/topic/axisymmetric-diffusion#post-9178 Fri, 07 Feb 2025 17:16:26 +0000 s148275 9178@http://www.k-wave.org/forum/ Hello, Is there an axisymmetric variant of the kWaveDiffusion code available? Kind regards, Frank Crispi99 on "Thermo-acoustic simulation with annular transducer" http://www.k-wave.org/forum/topic/thermo-acoustic-simulation-with-annular-transducer#post-9177 Thu, 06 Feb 2025 11:24:57 +0000 Crispi99 9177@http://www.k-wave.org/forum/ Hi Brad, hi community, I'm simulating the thermo-acoustic field of an annular transducer, in a heterogeneous medium (a water tank and a small piece of fat in the focal region, in order to evaluate the thermal lesion). The source pressure is 175 kPa, the frequency is 1.5 MHz. The problem is that I obtain an extremely high (wrong) field of temperature, up to 1000 °C with a heating period of 5 s. I obtain a max intensity of 4691.1406 W/cm^2 and a max heat volume rate Q of 1212.9214 W/cm^3. Thus, maybe Q is excessive but I cannot explain why. This is the formula that I have considered: alpha_np = db2neper(medium.alpha_coeff, medium.alpha_power) .* (2 * pi * freq).^medium.alpha_power; Q = alpha_np .* amp.^2 ./ (medium.density .* medium.sound_speed); Briefly, these are some parameters that I'm considering: % water medium parameters c0 = 1500; % sound speed [m/s] rho0 = 1000; % density [kg/m^3] BonA = 5.2; % non linear parameter of the medium (B over A) alpha_0 = 0; % frequency-dependent absorption coefficient alpha_y = 2; % power-law exponent of water % tissue medium parameters (fat tissue) c0_t = 1450; % sound speed of tissue [m/s] rho0_t = 911; % density of tissue [kg/m^3] BonA_t = 10; % non linear parameter of tissue (B over A) alpha_0_t = 0.6; % frequency-dependent absorption coefficient of tissue [dB/cm*MHz] alpha_y_t = 1.0861; % power-law exponent of tissue axial_size = 100e-3; % total grid size in the axial dimension [m] lateral_size = 120e-3; % total grid size in the lateral dimension [m] x_tissue = 30e-3; % tissue dimension along x coordinate [m] y_tissue = 120e-3; % tissue dimension along y coordinate [m] z_tissue = 120e-3; % tissue dimension along z coordinate [m] distance_t = 70e-3; % tissue distance from the transducer external surface [m] % thermal parameters T0 = 20; % water initial temperature [°C] thermal_conductivity = 0.6045; % water thermal conductivity [W/(m.K)] specific_heat = 4178; % water specific heat [J/(kg.K)] T0_t = 20; % tissue initial temperature [°C] thermal_conductivity_t = 0.211; % tissue thermal conductivity [W/(m.K)] specific_heat_t = 2348.33; % tissue specific heat [J/(kg.K)] on_time = 5; % source on time [s] off_time = 20; % source off time [s] dt_thermal = 0.1; % time step size of thermal simulation I thank you in advance for your attention! vich on "C++ kWaveDiffusion" http://www.k-wave.org/forum/topic/c-kwavediffusion#post-7040 Tue, 10 Sep 2019 03:06:50 +0000 vich 7040@http://www.k-wave.org/forum/ Hi, I was wondering if the modified Pennes bioheat equation (<code>kWaveDiffusion</code>) could be optimized in theory with compiled binaries in the same fashion as <code>kspaceFirstOrder3D</code>? If so, would it be feasible to use the C++ code for <code>kspaceFirstOrder3D</code> (e.g. <code>KSpaceFirstOrder3DSolver.cpp</code>) as a starting point for a C++-based <code>kWaveDiffusion</code>, or is the problem much more complicated than that? Thank you! otyunis on "Converting perfusion values from literature for use in k-Wave thermal sims" http://www.k-wave.org/forum/topic/converting-perfusion-values-from-literature-for-use-in-k-wave-thermal-sims#post-9042 Tue, 20 Feb 2024 21:05:50 +0000 otyunis 9042@http://www.k-wave.org/forum/ Hi, I want to ensure that I am specifying perfusion values that map to realistic physiological ranges. In the literature, I see that perfusion rates are typically quoted as a volume of blood that flows through either a volume of tissue OR mass of tissue per some amount of time (e.g., 1 to 20 mL blood/100g tissue/min). In k-Wave, I see that I can specify a (possibly heterogeneous) perfusion distribution as either a single matrix of perfusion coefficients or combination of 3 matrices (blood density, heat capacity, and perfusion rate). I've chosen the latter option, but am a bit confused on how to specify the units properly. May I request an explanation of how to convert these literature values of perfusion to the 1/s units that k-Wave expects? Do I need to take into account the grid spacing? Thank you! k-Wave is AWESOME! Best regards, Omar Wave_Base on "EM wave SAR (specific absorbation rate) W/Kg unit conversion to [dB/(MHz^y cm)]" http://www.k-wave.org/forum/topic/em-wave-sar-specific-absorbation-rate-wkg-unit-conversion-to-dbmhzy-cm#post-8994 Fri, 29 Dec 2023 05:26:26 +0000 Wave_Base 8994@http://www.k-wave.org/forum/ Dear all, Greetings! I'm trying to simulate thermoacoustic response in k-wave by using EM waves as a source (already done with CST software), From CST I have SAR values for a meadium/ phantom, however i need to generate acoustic signal based on SAR using K-wave, I have gone through all the response in this forum but i still need clarity to integrate CST SAR values which i have obtained as W/kg or dB(W/kg) unit scale, and how this can be embed into K-wave-'kspaceFirstOrder3D' function, as i can see i need to play with (<code>'medium.BonA', ' medium.alpha_power', 'medium.alpha_coeff', 'medium.alpha_mode'</code>), variables for adaptation, however unites that I need to deal with is [dB/(MHz^y cm)] but here i have SAR values in W/kg or dB(W/kg), how this can be mearged or any other alternitive ways? Pls suggest your valuable comments. Mor on "HITU analysis" http://www.k-wave.org/forum/topic/hitu-analysis#post-8945 Tue, 14 Nov 2023 10:32:19 +0000 Mor 8945@http://www.k-wave.org/forum/ has anyone tried performing HITU analysis (high-intensity therapeutic ultrasound)? we are using 1 rectangular transducer (unfocused) and I'm trying to build its heatmap. I would appreciate pieces of advice. ellianarv on "Converting temperature-time data to CEM43?" http://www.k-wave.org/forum/topic/converting-temperature-time-data-to-cem43#post-8857 Fri, 07 Jul 2023 18:40:57 +0000 ellianarv 8857@http://www.k-wave.org/forum/ Hi – I am pretty inexperienced with MATLAB and k-wave, so apologies if this question has a simple answer. I have a few data sets of temperature over time recorded during various ablations, and I'm looking to convert that data to a CEM 43 value for each treatment. I see CEM43 mentioned in various k-wave instructions, but I can't seem to figure out how to convert my data. Any help would be appreciated! a.stark on "kdiff object creating T matrix of NaN values" http://www.k-wave.org/forum/topic/kdiff-object-creating-t-matrix-of-nan-values#post-8708 Sun, 26 Feb 2023 22:42:59 +0000 a.stark 8708@http://www.k-wave.org/forum/ Hello there! In my code, I set up a tone burst signal that reaches a steady state; however, whenever I move onto the thermal simulation using kwavediffusion as an object and taking time steps, it results in T being computed as a matrix of NaN values, despite there being heat deposition values. In addition to this issue, I have also used kWaveDiffusion as an object in a for loop to take time steps when the signal is on and time steps when the signal is off to simulate a pulsed wave, and was wondering if that was okay. Thanks alcebrui on "Volume rate of heat deposition units" http://www.k-wave.org/forum/topic/volume-rate-of-heat-deposition-units#post-6925 Tue, 25 Jun 2019 19:09:51 +0000 alcebrui 6925@http://www.k-wave.org/forum/ First of all, thanks a lot for this amazing tool, it is really powerful and intuitive! My question is regarding the units shown for the volume rate of heat deposition in the plot of the example "Heating by a focused ultrasound transducer". In subplot 232, the y label indicates kW/cm^3 and in the corresponding code the variable for the volume rate of heat deposition is multiplied by a factor 1e-7, i.e., Q * 1e-7. Looking into the documentation Q is supposed to be in S.I units (W/m^3), so I would expect to have Q * 1e-9 if the desired units are kW/cm^3. Am I missing something there? Thanks! Manuel Vielma on "Volume rate of heat deposition as divergence of intensity" http://www.k-wave.org/forum/topic/volume-rate-of-heat-deposition-as-divergence-of-intensity#post-8692 Wed, 11 Jan 2023 00:26:48 +0000 Manuel Vielma 8692@http://www.k-wave.org/forum/ In equation 6 of their recent (October 2022) paper [1], Kleparnik, Zemcik, Treeby and Jaros define the volume rate of heat deposition, Q, as Q = -div(I_avg) = - d(Ix_avg)/dx - d(Iy_avg)/dy - d(Iz_avg)/dz (1), (the second equality has been introduced by me for explicitness). This is in line with early expressions such as e.g. the one given by Nyborg in [2]. To my knowledge, eq. (1) is not used in any of the k-wave examples. When I use it to compare with the temperature output that other expressions, such as Q = p^2/(density*sound_speed) (2) -- plane wave approximation or Q = alpha_np * abs(I_avg) (3), lead to in example_diff_focused_ultrasound_heating, I get very different results -- the ones given by (1) seem particularly odd. More precisely, what I am doing is running this example for the different Qs shown below: ============================ <pre><code>% Heating By A Focused Ultrasound Transducer [...] % set the sensor mask to cover the entire grid sensor.mask = ones(Nx, Ny); sensor.record = {'p', 'p_max_all', 'I_avg'}; %Include I_avg in recorded data [...] % ========================================================================= % CALCULATE HEATING % ========================================================================= % convert the absorption coefficient to nepers/m alpha_np = db2neper(medium.alpha_coeff, medium.alpha_power) * ... (2 * pi * freq).^medium.alpha_power; % extract the pressure amplitude at each position p = extractAmpPhase(sensor_data.p, 1/kgrid.dt, freq); % reshape the data, and calculate the volume rate of heat deposition p = reshape(p, Nx, Ny); %Q = alpha_np .* p.^2 ./ (medium.density .* medium.sound_speed); % Plane wave approximation I_avg_x = reshape(sensor_data.Ix_avg, Nx,Ny); I_avg_y = reshape(sensor_data.Iy_avg, Nx,Ny); %Q = alpha_np .* sqrt(I_avg_x.^2 + I_avg_y.^2); Q = -divergence(I_avg_x, I_avg_y); [...]</code></pre> ==================================== MY QUESTION IS: how should eq. (1) be used in a thermal simulation in k-wave? Any idea why seemingly inconsistent answers are obtained when using it in the aforementioned example? This post mentions the gradient (?) of the intensity, but no further details are provided: <a href="http://www.k-wave.org/forum/topic/doubts-converting-pressure-into-heat-deposition" rel="nofollow">http://www.k-wave.org/forum/topic/doubts-converting-pressure-into-heat-deposition</a> Many thanks in advance for any clarification you may provide. Manuel ______________________________ [1] Kleparnik, P., Zemcik, P., Treeby, B. E., & Jaros, J. (2022). On-the-Fly Calculation of Time-Averaged Acoustic Intensity in Time-Domain Ultrasound Simulations Using a k-Space Pseudospectral Method. IEEE transactions on ultrasonics, ferroelectrics, and frequency control, 69(10), 2917–2929. <a href="https://doi.org/10.1109/TUFFC.2022.3199173" rel="nofollow">https://doi.org/10.1109/TUFFC.2022.3199173</a> [2] Nyborg (1981). Heat generation by ultrasound in a relaxing medium. The Journal of the Acoustical Society of America 70, 310 (1981); <a href="https://doi.org/10.1121/1.386778" rel="nofollow">https://doi.org/10.1121/1.386778</a> s148275 on "Ultrasound heating and Intensity" http://www.k-wave.org/forum/topic/ultrasound-heating-and-intensity#post-8617 Tue, 04 Oct 2022 14:26:47 +0000 s148275 8617@http://www.k-wave.org/forum/ I am a bit confused about the equation that needs to be used for the heating rate Q. From my simulation I get the outputs: <code>sensor_data.p</code>, <code>sensor_data.ux_non_staggered</code> and <code>sensor_data.uy_non_staggered</code>. The following code is used to get the average intensity: <pre><code>ux_sgt = fourierShift(sensor_data.ux_non_staggered, 1/2) uy_sgt = fourierShift(sensor_data.uy_non_staggered, 1/2) Ix = sensor_data.p .* ux_sgt Iy = sensor_data.p .* uy_sgt I = sqrt(Ix^2+Iy^2) I_avg = mean(I, ndims(I))</code></pre> If I want to find the heat deposition Q, what is the equation that I need to use? For now I am using: <pre><code>alpha_np = db2neper(alpha_coeff_p_db,alpha_power) ... * (2 * pi * f0).^alpha_power; Q = alpha_np .* I_avg;</code></pre> However, sometimes <code>Q=2*alpha_np*I_avg</code> is used. What is the case here and why? Thanks in advance! gabrielcathoud on "Doubts converting pressure into heat deposition" http://www.k-wave.org/forum/topic/doubts-converting-pressure-into-heat-deposition#post-8565 Wed, 13 Jul 2022 14:51:26 +0000 gabrielcathoud 8565@http://www.k-wave.org/forum/ I'm simulating the heating generated by the use of ultrasound in the brain. The ultrasound simulation is working fine, but I am having difficulties with the temperature simulation. How do I convert pressure aplied into deposited heat? I found a conversion here on the forum, but I don't understand why: source.Q = alpha_np.*sensor_data.p_max_all.^2./(medium.density .* medium.sound_speed); I appreciate any help. The whole temperatura simulation is: clear source; medium.thermal_conductivity = thermal_conductivity_water*ones(Nx,Ny); medium.thermal_conductivity(head_mask==1) = thermal_conductivity_bone; medium.specific_heat = specific_heat_water*ones(Nx,Ny); medium.specific_heat(head_mask==1) = specific_heat_bone; alpha_np = db2neper(medium.alpha_coeff, medium.alpha_power)... *(2 * pi * source_freq).^medium.alpha_power; dt = 0.001; source.Q = alpha_np.*sensor_data.p_max_all.^2./(medium.density .* medium.sound_speed); source.T0 = 37; input_args = {'PlotSim', true}; kdiff = kWaveDiffusion(kgrid, medium, source, [], input_args{:}); on_time = duration_of_stimulus; % [s] off_time = 60; % [s] kdiff.takeTimeStep(round(on_time / dt), dt); T1 = kdiff.T; kdiff.Q = 0; dt = 0.1; kdiff.takeTimeStep(round(off_time / dt), dt); T2 = kdiff.T; kate1233 on "The difference between 2 dimension and 3dimension" http://www.k-wave.org/forum/topic/the-difference-between-2-dimension-and-3dimension#post-8338 Tue, 09 Nov 2021 07:16:17 +0000 kate1233 8338@http://www.k-wave.org/forum/ Dear Bradley Treeby I am a student studying HIFU simulation .Rencently ,I use k-wave to simulate the two dimension and the three dimension of laminated tissue therapeutic process. Now ,I found the the temperature and heat dose with great differneces between 2 dimension and 3 dimension .I am confused about this result because thermal damage of three dimension is about 624.8mm^2 and two dimension is about 10.175mm^2（three periods with 1 second on time and 3 second off time ）at the same time the temperature of three is over 968°C only in one period and the temperature of two is over 73.68°C in one period.I also simulate in water between the two and three dimension and the difference of focus pressure is 15 times . 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> deepak11 on "Differences in focus pressure between k-Wave and analytical transducer model" http://www.k-wave.org/forum/topic/differences-in-focus-pressure-between-k-wave-and-analytical-transducer-model#post-8244 Tue, 20 Jul 2021 23:51:29 +0000 deepak11 8244@http://www.k-wave.org/forum/ Hi, First of all thank you for building this amazing toolbox. I have simulated the bowel transducer in 3D and then compared the focus pressure amplitude with the analytical model of a concave spherical transducer (from Theory of Focusing Radiators by H.T. O'Neil). I also checked this analytical model results from the Sonic Concepts transducer (H-276) and it was correct. I couldn't find the mistake in my code I think I have done everything correctly but the results shouldn't deviate that much (like the peak pressure at the focus must be 58 MHz but the simulated pressure is 46.8 MHz ). I have posted the code below: (`% define the grid parameters 3D Nx = 150; % [grid points] Ny = 150; % [grid points] Nz = 92; dx = 0.1e-3; % [m] dy = 0.1e-3; % [m] dz = 0.1e-3; % [m] % create the computational grid kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); % define the properties of the propagation medium water medium.sound_speed = 1500; % [m/s] medium.density = 1000; % [kg/m^3] medium.alpha_coeff = 0.025; % [dB/(MHz^y cm)] medium.alpha_power = 1.5; % define parameters grid_size = [150,150,92]; bowl_pos = [Nx/2, Ny/2, 2]; diameter = 15e-3; % [m] radius = 8.5e-3; % [m] focus_pos = [Nx/2, Ny/2, 87]; freq = 1.5e6; % [Hz] amp = 2.06e6; % [Pa] % create bowl source.p_mask = makeBowl(grid_size, bowl_pos, round(radius / dx), round(diameter / dx)+1, focus_pos, 'Plot', true); % calculate the time step using an integer number of points per period ppw = max(medium.sound_speed(:)) / (freq * dx); % points per wavelength cfl = 0.3; % cfl number ppp = ceil(ppw / cfl); % points per period T = 1 / freq; % period [s] dt = T / ppp; % time step [s] % calculate the number of time steps to reach steady state t_end = sqrt( kgrid.x_size.^2 + kgrid.y_size.^2 + kgrid.z_size.^2 ) / max(medium.sound_speed(:)); Nt = round(t_end / dt); % create the time array kgrid.setTime(Nt, dt); % define the input signal source.p = createCWSignals(kgrid.t_array, freq, amp, 0); % set the sensor mask to cover the entire grid sensor.mask = ones(Nx, Ny, Nz); sensor.record = {'p','p_rms','u','u_rms', 'I','I_avg'}; % record the last 3 cycles in steady state num_periods = 3; T_points = round(num_periods * T / kgrid.dt); sensor.record_start_index = Nt - T_points + 1; % set the input arguements input_args = {'PMLInside', false, 'PlotPML', false, 'DisplayMask', ... 'off', 'PlotScale', [-1, 1] * amp}; % run the acoustic simulation sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:}); % convert the absorption coefficient to nepers/m alpha_np = db2neper(medium.alpha_coeff, medium.alpha_power) * ... (2 * pi * freq).^medium.alpha_power; % extract the pressure amplitude at each position p = extractAmpPhase(sensor_data.p, 1/kgrid.dt, freq); % reshape the data, and calculate the volume rate of heat deposition p = reshape(p, Nx, Ny, Nz); Alexander on "Temperature Changes with Grid Size?" http://www.k-wave.org/forum/topic/temperature-changes-with-grid-size#post-8086 Thu, 18 Mar 2021 18:15:01 +0000 Alexander 8086@http://www.k-wave.org/forum/ I came across an issue running thermal sims with a ring transducer. It would appear that if I change only Nx/Ny that my temperature (T1 and T2) changes. It may be an issue with how extractAmpPhase(...) is working/ how Q is being calculated. But I just took the default example and changed it to use a ring array. Note: When I change Nx/Ny I also change dx/dy so the grid is the same size in [mm]. Also a pdf of the results and the code can be see at the link below. <a href="https://waynestateprod-my.sharepoint.com/:f:/g/personal/eu7209_wayne_edu/Et6Pvc10AstKmGzRdtM-KJsBcNYd4TW9l1P59UEWsQqLPA?e=nRqyJZ" rel="nofollow">https://waynestateprod-my.sharepoint.com/:f:/g/personal/eu7209_wayne_edu/Et6Pvc10AstKmGzRdtM-KJsBcNYd4TW9l1P59UEWsQqLPA?e=nRqyJZ</a> Andrea Addison on "Staircase Problem" http://www.k-wave.org/forum/topic/staircase-problem#post-7991 Sun, 27 Dec 2020 19:08:55 +0000 Andrea Addison 7991@http://www.k-wave.org/forum/ Hello Dr.Treeby, I want to use a plane piston transducer with incidence angles at different degrees. When I use an angled transducer mask, I observe some staircase issue which doesn't represent angled transducer accurately. Then I encountered this example (<a href="http://www.k-wave.org/documentation/example_tvsp_steering_linear_array.php#heading2)" rel="nofollow">http://www.k-wave.org/documentation/example_tvsp_steering_linear_array.php#heading2)</a>. I think this might be a good solution but I was wondering why you wouldn't prefer to use the steering method in your study: "The contribution of shear wave absorption to ultrasound heating in bones: Coupled elastic and thermal modeling". Do you think that the physical mask and electronic steering are different representations? Best regards, Andrea DLam on "Memory Issues and a Possible Workaround?" http://www.k-wave.org/forum/topic/memory-issues-and-a-possible-workaround#post-7718 Sat, 01 Aug 2020 00:15:48 +0000 DLam 7718@http://www.k-wave.org/forum/ Hi, I'm planning to use the kWaveDiffusion to simulate thermal diffusion. I will be switching between alternating between several source.Q arrays as I step through time, and I am using timesteps of 0.05 sec to simulate for a total time of 16 sec. However, currently the simulation runs out of memory -- I think this is because of the sensor_data field. If all I want at the end is the CEM43 lesion map and the final temperature map, should I just set sensor.mask to be FALSE everywhere, or is there something else I can do? Best regards, Daniel renardChien on "Perfusion Value in BioHeatExact" http://www.k-wave.org/forum/topic/perfusion-value-in-bioheatexact#post-7695 Thu, 09 Jul 2020 13:05:01 +0000 renardChien 7695@http://www.k-wave.org/forum/ I am confused about the perfusion value of the Bioheatexact() function. The function seems to accept a constant value, a 1x3 (x,y,z) array, and a 3D (x,y,z) matrix. For HIFU, the heat maps in a plane perpendicular to the HIFU beam are always symmetrical. I want to input a perfusion value and have the heat map shape change in the direction of fluid flow. Can you provide any additional information that will help me implement this feature? tahere on "heat generated by HIFU" http://www.k-wave.org/forum/topic/heat-generated-by-hifu#post-6222 Thu, 07 Dec 2017 23:17:39 +0000 tahere 6222@http://www.k-wave.org/forum/ hi dear professor bradley treeby I considered the heat source as a 3-d Gaussian function and I want to know the temperature in the xy-plane in front of HIFU.but I think the temperature after heating and cooling is unreasonable.may you tell me what is the problem? thanks in advance here is my code %% set the size of the perfectly matched layer (PML) PML_X_SIZE = 10; % [grid points] PML_Y_SIZE = 10; % [grid points] PML_Z_SIZE = 10; % [grid points] %% set total number of grid points not including the PML Nx = 120- 2*PML_X_SIZE; % [grid points] Ny = 120- 2*PML_Y_SIZE; % [grid points] Nz = 120- 2*PML_Z_SIZE; % [grid points] %% set desired grid size in the x-direction not including the PML x = .1; % [m] %% calculate the spacing between the grid points dx = x/Nx; % [m] dy = dx; % [m] dz = dx; % [m] %% create the k-space grid kgrid =kWaveGrid(Nx, dx, Ny, dy, Nz, dz); x1=linspace(0,.1,100); y1=x1; z1=x1; [X,Y,Z] = meshgrid(x1,y1,z1); Q=((2.34e4)*exp(-power(((X-.05)/.004284),2))).*((2.34e4)*exp(-power(((Y-.05)/.004284),2))).*((2.008e4)*exp(-power(((Z-.01153)/.02147),2))); clear medium source sensor; % set the background temperature and heating term source.Q = Q; source.T0 = 37; % define medium properties related to diffusion medium.density = 1020; % [kg/m^3] medium.thermal_conductivity = 0.5; % [W/(m.K)] medium.specific_heat = 3600; % [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; % ========================================================================= % VISUALISATION % ========================================================================= % plot the thermal dose and lesion map figure; % plot the volume rate of heat deposition subplot(3, 1, 1); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, Q(:,:,50) * 1e-7); h = colorbar; xlabel(h, '[kW/cm^3]'); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Volume Rate Of Heat Deposition'); % plot the temperature after heating subplot(3, 1, 2); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, T1(:,:,50)); h = colorbar; xlabel(h, '[degC]'); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Temperature After Heating'); % plot the temperature after cooling subplot(3, 1, 3); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, T2(:,:,50)); h = colorbar; xlabel(h, '[degC]'); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Temperature After Cooling');