Using intensity as heat source term (nonlinear acoustics)

k-Wave User Forum » General Acoustics » Nonlinear Acoustics

Using intensity as heat source term (nonlinear acoustics)

(14 posts) (3 voices)

Started 11 years ago by Anthony
Latest reply from Bradley Treeby

Anthony

Member
Joined: Apr '13
Posts: 96

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

Posted 11 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,643

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(<I>), 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 sensor.record = {'I_avg'} and then setting sensor.record_start_index 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.

Posted 11 years ago #
Anthony

Member
Joined: Apr '13
Posts: 96

Thanks a lot for your reply Bradley!
Anthony

Posted 11 years ago #
Anthony

Member
Joined: Apr '13
Posts: 96

Hi Bradley,

Just to be sure, when you write Q = - div(<I>), 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?

Posted 11 years ago #
Anthony

Member
Joined: Apr '13
Posts: 96

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. (see example)

From your experience, is this normal? I must probably miss some point...
Regards,

Posted 11 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,643

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.

Posted 11 years ago #
Anthony

Member
Joined: Apr '13
Posts: 96

Posted 11 years ago #
Anthony

Member
Joined: Apr '13
Posts: 96

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'])

Posted 11 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,643

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 here). We've changed the way this is calculated ready for V1.1.

Brad.

Posted 11 years ago #
Anthony

Member
Joined: Apr '13
Posts: 96

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

Posted 11 years ago #
Anthony

Member
Joined: Apr '13
Posts: 96

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)')

Posted 11 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,643

Hi Anthony,

Thanks for your code - I see what you mean. You could try using the gradientSpect 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.

Posted 11 years ago #
DLam

Member
Joined: Oct '16
Posts: 19

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(<I>) in this case? I suppose the vector field will have to be assumed somehow, e.g. to be along the transducer beam axis.

Posted 4 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,643

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.

Posted 4 years ago #

RSS feed for this topic

Reply

You must log in to post.

k-Wave

A MATLAB toolbox for the time-domain simulation of acoustic wave fields

Using intensity as heat source term (nonlinear acoustics)

Reply

A MATLAB toolbox for the time-domain
simulation of acoustic wave fields