k-Wave User Forum » Topic: Using intensity as heat source term (nonlinear acoustics)

k-Wave User Forum » Topic: Using intensity as heat source term (nonlinear acoustics) http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:28:31 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Using intensity as heat source term (nonlinear acoustics)" http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics#post-7232 Wed, 12 Feb 2020 16:19:48 +0000 Bradley Treeby 7232@http://www.k-wave.org/forum/ Hi DLam, I'd suggest recording a time varying pressure signal for an integer number of periods at the end of the simulation. Then spectrally decompose the signal, and sum the absorption coefficient times intensity (calculated using p^2) at each harmonic. Brad. DLam on "Using intensity as heat source term (nonlinear acoustics)" http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics#post-7200 Tue, 28 Jan 2020 17:26:28 +0000 DLam 7200@http://www.k-wave.org/forum/ Hi Brad, I'm replying to this thread since the work here is similar to what I wanted to ask. I have been using the C++ version of the simulator, where there is no I_avg output flag available. Currently, I am using the approximation that I_avg = (P_max^2)/(2*density*soundspeed) to calculate the scalar intensity array. Can I still use Q = -div() in this case? I suppose the vector field will have to be assumed somehow, e.g. to be along the transducer beam axis. Bradley Treeby on "Using intensity as heat source term (nonlinear acoustics)" http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics#post-4096 Wed, 20 Nov 2013 10:24:38 +0000 Bradley Treeby 4096@http://www.k-wave.org/forum/ Hi Anthony, Thanks for your code - I see what you mean. You could try using the <code>gradientSpect</code> function in k-Wave to see if that makes a difference to your calculation of the divergence. This computes the derivative in the same way as within the simulation functions. In regards to the accuracy of the intensity calculation, I've added some code to shift the velocity to the non-staggered spatial and temporal grid before multiplying by the pressure, and this seems to improve things. This will go into the V1.1 release. Brad. Anthony on "Using intensity as heat source term (nonlinear acoustics)" http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics#post-4092 Fri, 15 Nov 2013 15:45:57 +0000 Anthony 4092@http://www.k-wave.org/forum/ This small piece of code illustrates what I have just written (you just have to specify the correct path in the 3 first lines). Anthony ********************************* CODE ****************************** clear all close all clc infilename = 'C:\Users\Grisey\Desktop\checkSource\input.h5'; outfilename = 'C:\Users\Grisey\Desktop\checkSource\output.h5'; pathToKwave = '"C:\Users\Grisey\k-wave\k-wave-toolbox-version-1.0-cpp-windows-binaries\Windows Binaries\kspaceFirstOrder3D-OMP.exe"'; %(don't forget the ") %% Paramètres myPMLsize = 20; medium.density = 1000; medium.sound_speed = 1500; source_mag = 1e5; %% Grille de calcul f = 1e5; N = 128; h = 1e-3; kgrid = makeGrid(N, h, N, h, N, h); %% Definition du transducteur R = 10e-3; transducerMask = zeros(N,N,N); xgrid = h*((1:N)-N/2); ygrid = h*((1:N)-N/2); zgrid = h*((1:N)-N/2); [X,Y,Z] = meshgrid(xgrid,ygrid,zgrid); Rmap = (X.^2+Y.^2+Z.^2).^0.5; transducerMask(Rmap<=R)=1; %% Simulation % create the time array [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed); periode = 1/f; recStartTime = numel(kgrid.t_array)-2*round(periode/dt)+1; source.p_mask = logical(transducerMask); source.p = source_mag*sin(2*pi*f*kgrid.t_array); source.p_mode = 'dirichlet'; % assign the input options input_args = {'PlotPML', true, 'PMLSize', myPMLsize,'Smooth',[false true true]}; % run the simulation toto = ones(N, N, N); sensor.mask = logical(toto); clear toto; kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:},'SaveToDisk',infilename); status = system([pathToKwave ' -i ' infilename ' -o ' outfilename ' --p_rms --I_avg -s' num2str(recStartTime)]); %% Lecture du résultat p_rms = reshape(h5read(outfilename,'/p_rms'),N,N,N); Ix_avg = reshape(h5read(outfilename,'/Ix_avg'),N,N,N); Iy_avg = reshape(h5read(outfilename,'/Iy_avg'),N,N,N); Iz_avg = reshape(h5read(outfilename,'/Iz_avg'),N,N,N); H = -divergence(Ix_avg,Iy_avg,Iz_avg)/h; % Calcul de Ir Ianalytical = source_mag^2*R^2./(2*medium.sound_speed*medium.density.*Rmap.^2); theta = acos(Z./Rmap); phi = atan(Y./X)+pi*(X<0); Ix_analytical = Ianalytical.*cos(phi).*sin(theta); Iy_analytical = Ianalytical.*sin(phi).*sin(theta); Iz_analytical = Ianalytical.*cos(theta); figure(1) subplot(1,3,1) imagesc(squeeze(H(round(N/2),:,:)).*(1-squeeze(transducerMask(round(N/2),:,:)))) caxis([-0.1e5 2e4]) title('-div(k-wave intensity) (order 1 uncentered divergence)') subplot(1,3,2) Hanalytical = -divergence(Ix_analytical,Iy_analytical,Iz_analytical)/h; imagesc(squeeze(Hanalytical(round(N/2),:,:)).*(1-squeeze(transducerMask(round(N/2),:,:)))) caxis([-0.1e5 2e4]) title('-div(analytical intensity) (order 1 uncentered divergence)') subplot(1,3,3) divx = (Ix_avg(3:end,2:end-1,2:end-1)-Ix_avg(1:end-2,2:end-1,2:end-1))/2/h; divy = (Iy_avg(2:end-1,3:end,2:end-1)-Iy_avg(2:end-1,1:end-2,2:end-1))/2/h; divz = (Iz_avg(2:end-1,2:end-1,3:end)-Iz_avg(2:end-1,2:end-1,1:end-2))/2/h; HnumericalOrder2 = -divx-divy-divz; imagesc(squeeze(HnumericalOrder2 (round(N/2),:,:)).*(1-squeeze(transducerMask(round(N/2),2:end-1,2:end-1)))) caxis([-0.1e5 2e4]) title('-div(k-wave intensity) (order 2 centered divergence)') Anthony on "Using intensity as heat source term (nonlinear acoustics)" http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics#post-4091 Fri, 15 Nov 2013 13:54:18 +0000 Anthony 4091@http://www.k-wave.org/forum/ Hi Brad, According to some tests I have made, I think it comes from the numerical calculation of the divergence. In the simple case of a radiating sphere, with a linear medium and no attenuation, -div(pu)=0 theoretically. If I compute it numerically with an uncentered scheme of order 1, I get strong positive and negative zones along the axis of the cartesian grid. If I use a 2 order centered scheme, it becomes much better. However, it does not seem sufficient to solve the problem in the case of a focused field... Is it possible for you to try to do the calculation of the divergence in the spectral domain in k-wave in order to see if it solves the problem? Best regards, Anthony Bradley Treeby on "Using intensity as heat source term (nonlinear acoustics)" http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics#post-4047 Thu, 31 Oct 2013 20:12:43 +0000 Bradley Treeby 4047@http://www.k-wave.org/forum/ Hi Anthony, We are still looking into the heat deposition issue. We've had a few blackboard sessions, but I think more questions have been raised than answered! Certainly you wouldn't expect to see a heat sink, but exactly how the heat deposition should be computed depends partly on where the energy loss terms appear in the conservation equations. If we get any closer, I'll let you know. Regarding the variation in the expected power, this might be due to approximations in the way intensity is calculated in V1.0 (see <a href="http://www.k-wave.org/forum/topic/defining-incident-acoustic-plane-wave-intensity#post-3723">here</a>). We've changed the way this is calculated ready for V1.1. Brad. Anthony on "Using intensity as heat source term (nonlinear acoustics)" http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics#post-4044 Thu, 31 Oct 2013 17:35:17 +0000 Anthony 4044@http://www.k-wave.org/forum/ Hi Brad, I must admit that I couldn't solve the mystery (still using 2*alpha*I in the meantime). Did you have some new insight about this? I wrote a piece of code that shows the problem, if someone has an idea... By the way, the acoustic power I get by integrating intensity over a plane is smaller than the value I should theoretically get : I guess it comes from the way I compute the amplitude of the velocity source : u0 = sqrt( 2 * acousticPower / ( rho * soundSpeed * transducerSurface ) ) Is there another way to compute it (without the plane wave approximation, maybe...)? Best regards, Anthony ******************************** SAMPLE CODE ************************************ clear all close all clc %% Files work_dir = uigetdir('Set working directory (where the files will be saved) '); kspaceFirstOrder3D_OMP_path = uigetdir('Where is kspaceFirstOrder3D-OMP.exe?'); acousticSimulationInputFileName = 'input.h5'; acousticSimulationOutputFileName = 'output.h5'; %% Parameters % Source acousticPower = 80; % [W] source_freq = 2.5e5; %[Hz] % Physical parameters rho = 1000; soundSpeed = 1500; %% Geometry R = 40e-3; % [m] Rap = 30e-3; % [m] h = 0.5e-3; % [m] N = ceil(62e-3/h)+20; % acoustic simulation grid size %% k-wave parameters % Domain geometry myPMLsize = 10; % Source transducerSurface = 2*pi*R*(R-sqrt(R^2-Rap^2)); source_mag = sqrt(2 * acousticPower / (rho * soundSpeed * transducerSurface)); %velocity source kgrid = makeGrid(N, h, N, h, N, h); % define a curved transducer element xCenter = round(N/2); yCenter = round(N/2); zCenter = round(R/h)+2; transducerMask = zeros(N,N,N); xGrid = (1:N)'*ones(1,N); yGrid = ones(N,1)*(1:N); for zmap=1:(zCenter-(sqrt(R^2-Rap^2)/h)); rmap2D = round(((xGrid-xCenter).^2+(yGrid-yCenter).^2+(zmap-zCenter).^2).^0.5); transducerMask(:,:,zmap) = 0+(rmap2D==round((R)/h)); end clear xGrid yGrid rmap2D % voxelPlot(transducerMask) % define the properties of the propagation medium medium.density = 1000; % [kg/m3] medium.sound_speed = 1500; % [m/s] % create the time array [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed); periode = 1/source_freq; recStartTime = numel(kgrid.t_array)-round(periode/dt)+1; source.u_mask = logical(transducerMask); % define a time varying sinusoidal source [I,J,K] = ind2sub([N,N,N],find(source.u_mask)); I = xCenter - I; J = yCenter - J; K = zCenter - K; distCenter = (I.^2+J.^2+K.^2).^0.5; timeSine = sin(2*pi*source_freq*kgrid.t_array); source_ux0 = abs(I./distCenter)*source_mag; source_uy0 = abs(J./distCenter)*source_mag; source_uz0 = abs(K./distCenter)*source_mag; source.ux = source_ux0*timeSine; source.uy = source_uy0*timeSine; source.uz = source_uz0*timeSine; % assign the input options input_args = {'PlotPML', true, 'PMLSize', myPMLsize,'Smooth',[false true true]}; % run the simulation sensor.mask = logical(ones(N, N, N)); kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:},'SaveToDisk',[work_dir acousticSimulationInputFileName]); clear kgrid medium source sensor status = system(['"' kspaceFirstOrder3D_OMP_path '\kspaceFirstOrder3D-OMP.exe" -i ' [work_dir acousticSimulationInputFileName] ' -o ' [work_dir acousticSimulationOutputFileName] ' --I_avg -s' num2str(recStartTime)]); %% Load the results Ix_avg = h5read([work_dir acousticSimulationOutputFileName],'/Ix_avg'); Ix_avg = reshape(Ix_avg,N,N,N); Iy_avg = h5read([work_dir acousticSimulationOutputFileName],'/Iy_avg'); Iy_avg = reshape(Iy_avg,N,N,N); Iz_avg = h5read([work_dir acousticSimulationOutputFileName],'/Iz_avg'); Iz_avg = reshape(Iz_avg,N,N,N); %% Compute heat source term heatSource = -divergence(Ix_avg,Iy_avg,Iz_avg)/h; figure(1) imagesc(squeeze(heatSource(:,round(N/2),:))') xlabel('X') ylabel('Z') title('Heat source term : -div(I_{avg})') % Acoustic power Pac_from_I = max(squeeze(sum(sum(Iz_avg(myPMLsize:end-myPMLsize,myPMLsize:end-myPMLsize,myPMLsize:end-myPMLsize),2),1))*h^2); disp(['Acoustic power : ' num2str(Pac_from_I) 'W']) Anthony on "Using intensity as heat source term (nonlinear acoustics)" http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics#post-3829 Fri, 26 Jul 2013 09:27:48 +0000 Anthony 3829@http://www.k-wave.org/forum/ Bradley Treeby on "Using intensity as heat source term (nonlinear acoustics)" http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics#post-3825 Thu, 25 Jul 2013 18:46:45 +0000 Bradley Treeby 3825@http://www.k-wave.org/forum/ Hi Anthony, Interesting question! Leave it with us for a short time, I'd like to check a few things before I give you an answer. Brad. Anthony on "Using intensity as heat source term (nonlinear acoustics)" http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics#post-3818 Thu, 25 Jul 2013 14:19:15 +0000 Anthony 3818@http://www.k-wave.org/forum/ I realize my question is not clear. I understand that the formula simply comes from the Green formula, but... when I apply it to a focused field, I get a strong heat source right before the focal point and a strong (but weaker due to absorption) heat sink right after the focal point, what results in a strange thermal behaviour. (<a href="http://www.casimages.com/img.php?i=130725031106339995.jpg" title="example">see example</a>) From your experience, is this normal? I must probably miss some point... Regards, Anthony on "Using intensity as heat source term (nonlinear acoustics)" http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics#post-3802 Tue, 23 Jul 2013 08:21:37 +0000 Anthony 3802@http://www.k-wave.org/forum/ Hi Bradley, Just to be sure, when you write Q = - div(), the source term become negative at some places, is it normal? If yes, how do you physically explain it? Anthony PS : I am sorry, I suppose it is probably in Pierce's book but I don't have it and there is nothing in the book of Hamilton & Blackstock... Do you know another ref about this particular topic? Anthony on "Using intensity as heat source term (nonlinear acoustics)" http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics#post-3645 Thu, 16 May 2013 10:23:39 +0000 Anthony 3645@http://www.k-wave.org/forum/ Thanks a lot for your reply Bradley! Anthony Bradley Treeby on "Using intensity as heat source term (nonlinear acoustics)" http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics#post-3639 Wed, 15 May 2013 14:33:20 +0000 Bradley Treeby 3639@http://www.k-wave.org/forum/ Hi Anthony, The formula you mention is only valid for a single-frequency plane wave. In the more general case, you can compute the volume rate of heat deposition due to acoustic absorption from the acoustic intensity using Q = - div(), where I is the acoustic intensity vector and <> denotes the time-average over an acoustic period. You can use k-Wave to return the average acoustic intensity by setting <code>sensor.record = {'I_avg'}</code> and then setting <code>sensor.record_start_index</code> to the appropriate value to force the average to be over an acoustic period. In the linear case with a single frequency, this should give something close to Q = alpha * p0^2 / rho0 * c0, where p0 is the acoustic amplitude. Brad. Anthony on "Using intensity as heat source term (nonlinear acoustics)" http://www.k-wave.org/forum/topic/using-intensity-as-heat-source-term-nonlinear-acoustics#post-3638 Wed, 15 May 2013 09:54:35 +0000 Anthony 3638@http://www.k-wave.org/forum/ Hi everyone, I have a simple question about using the Intensity as a source term in the bioheat equation. In the linear case, there is no problem using the 2*alpha*I classical formula but, in the case of nonlinear acoustics, as alpha varies with frequency, I am a bit confused... Do I have to apply a temporal Fourier transform to the intensity field in order to separate the different frequencies and multiply each of them by the corresponding alpha ? By the way, thanks a lot for k-wave ! It is fast and user-friendly... awesome :-) Anthony