Analysis of the power loss along propagation in lossless medium

k-Wave User Forum » General Acoustics » Acoustic Simulations

Analysis of the power loss along propagation in lossless medium

(9 posts) (2 voices)

Started 8 years ago by Anthony
Latest reply from Bradley Treeby

Anthony

Member
Joined: Apr '13
Posts: 96

Hi,

I am working on the simulation of HIFU treatments and I wanted to check the power which is dissipated by the numerical scheme itself (in the case of nonlinear propagation) in order to correct the deposited energy in the absorbing case.

Let z denote the main propagation of my field (represented here for information: https://www.dropbox.com/s/089etiubtkhkkwt/XZ-field.png?dl=0).
I computed the integral of I_z on each xy-plane and plotted the derivative of the power with respect to z (see https://www.dropbox.com/s/lgynn4yaa3esg6o/deposited_power.png?dl=0).

My question is: before the focus, the energy conservation is approximately OK but downstream to it, there is some significant power loss: could you tell me why?

A few precisions:
- of course, I removed the power which leaves the domain by the sides (using I_x and I_y).
- my source frequency is 3 MHz and my grid step is 40 microns, which is OK to properly propagate several harmonics.
--> At first, I thought that the energy loss was due to the generation of harmonics which are not supported by the grid but, even at the end of my domain, the high harmonics are not carrying that much energy (a few percents). More info about this in PS.
- I set my source as a Dirichlet boundary condition. I first simulated my acoustic field with the linear equations up to 6 mm upstream to the focus and then used this as my source term for the nonlinear simulation because I couldn't simulate the whole field in the nonlinear case due to memory limitations. By the way, why do I get these power oscillations? (my medium is homogeneous).
- I didn't forget to unstagger the velocity fields in time using timeShift before computing the intensity.

Don't hesitate to ask me for more details or figures if I am not clear. I could also send you the acoustic field itself (18 GB) or only the intensity values I computed (~100MB along each direction).

Best regards and sorry to nitpick ;-)
Anthony

PS: in order to assess the harmonic content of the field, I computed the module of the time-fft of my pressure field. Then for each voxel, I divided the module of each harmonics by the sum of the modules (taking the square gives approximately the same results). Finally, I took the mean of this proportion over each xy-plane and plotted this curve: https://www.dropbox.com/s/9ee9vfy4vd6w8lo/harmonic_content.png?dl=0
Little detail: in each xy-plane, I only considered the pixels where the rms pressure was greater than -6dB of the max in this plane i.e. only the pixels "inside the beam".
--> As you can see, the high harmonics are not very important... (well, I admit that my metric is not wonderful, if you want me to measure this using another formula, feel free to ask).

Posted 8 years ago #

Developer
Joined: Mar '10
Posts: 1,640

Hi Anthony,

I'm not sure if I follow exactly what you are trying to do? If you take a simple measure of the kinetic and potential energy, does that work out? I've written a simple script to demonstrate what I mean. This calculates the kinetic and potential energy, and the Lagrangian density. Once the source has been added, the energy is constant (there is no dissipation by the numerical scheme until the waves reach the PML), and the Lagrangian density is zero as expected. I would expect the same in nonlinear case.

Hope that helps,

Brad.

clear all;

% create the computational grid
Nx = 128;
Ny = 128;
dx = 0.1e-3;
dy = 0.1e-3;
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500;
medium.density = 1000;

% create time array
kgrid.t_array = makeTime(kgrid, medium.sound_speed);

% create initial pressure distribution using makeDisc
source.p_mask = zeros(Nx, Ny);
source.p_mask(Nx/2, 50:60) = 1;

% define source
source_mag = 10e6;
source.p = source_mag*toneBurst(1/kgrid.dt, 3e6, 3);

% zero pad and filter
source.p = filterTimeSeries(kgrid, medium, [source.p, zeros(1, length(source.p))]);

% define sensor
sensor.mask = ones(Nx, Ny);
sensor.record = {'p', 'u_non_staggered'};

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,...
    'DisplayMask', 'off', 'PlotScale', [-0.5, 0.5]*source_mag);

% time shift
ux = timeShift(sensor_data.ux_non_staggered, kgrid.dt);
uy = timeShift(sensor_data.uy_non_staggered, kgrid.dt);

% compute velocity magnitude
u = sqrt(ux.^2 + uy.^2);

% compute Lagrangian density
kinetic = sum(u.^2, 1) * 0.5 * medium.density;
potential = sum(sensor_data.p.^2, 1) / (2*medium.density*medium.sound_speed.^2);
Lagrangian_density = kinetic - potential;

% plot change in energy over time
figure;
subplot(2, 1, 1);
plot(kgrid.t_array, kinetic, 'b-');
hold on;
plot(kgrid.t_array, potential, 'r-');
legend('kinetic', 'potential')
title('Energy Components');

subplot(2, 1, 2);
plot(kgrid.t_array, Lagrangian_density, 'k-');
title('Lagrangian Density');

Posted 8 years ago #

Anthony

Member
Joined: Apr '13
Posts: 96

Hi Brad,

Thanks for your very convincing answer. However, I still have a problem of energy conservation. I slightly modified your code in order to show you what I mean. There is probably something wrong about my calculation but, even with alpha = 0, I calculate a power loss of 2.6 W along the beam path, which represents 30 percents of the power loss in the case of a classical absorbing material (liver) in the same situation... As you can see, I take into account the power leaving the domain through the sides. By the way, you are right, nonlinearity does not change anything.

Best regards,
Anthony

clear all
close all
clc

%% PARAMETERS TO CHANGE
medium.alpha_coeff =  0; % liver: 0.37
medium.alpha_power =  0; % liver: 1.266

% create the computational grid
Nx = 128-40;
Ny = 128-40;
dx = 0.1e-3;
dy = 0.1e-3;
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500;
medium.density = 1000;
medium.BonA = 0;

% create time array
kgrid.t_array = makeTime(kgrid, medium.sound_speed);

% create initial pressure distribution using makeDisc
R = Ny*0.45;
Rap = 28/38*R;
source.p_mask = zeros(Nx, Ny);
xc = round(Nx/2);
yc = round(Ny/2);
xgrid = (1:Nx)-xc;
ygrid = ((1:Ny)-yc)';
rgrid = sqrt(bsxfun(@plus,xgrid.^2,ygrid.^2));
source.p_mask(logical(bsxfun(@times,(abs(rgrid-R)<1/2), bsxfun(@times,ygrid<0,abs(xgrid)<Rap)))) = 1;
source.p_mode = 'dirichlet';

figure(1)
clf
imagesc(source.p_mask)
axis equal

% define source
source_mag = 0.1e6;
source_f0 = 3e6;
source.p = source_mag*cos(2*pi*kgrid.t_array*source_f0);

% zero pad and filter
source.p = filterTimeSeries(kgrid, medium, [source.p, zeros(1, length(source.p))]);

% define sensor
sensor.mask = ones(Nx, Ny);
sensor.record = {'p', 'u_non_staggered'};

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,...
    'DisplayMask', 'off', 'PlotScale', [-0.5, 0.5]*source_mag,'PMLInside',false,'PlotPML',false);

% time shift
ux = timeShift(sensor_data.ux_non_staggered, kgrid.dt);
uy = timeShift(sensor_data.uy_non_staggered, kgrid.dt);

% compute velocity magnitude
u = sqrt(ux.^2 + uy.^2);

% Intensity
NpointsToKeep = 50; % in order to have an integer number of points)
Ix = reshape(sensor_data.p.*ux,[Nx,Ny,size(ux,2)]);
Iy = reshape(sensor_data.p.*uy,[Nx,Ny,size(uy,2)]);
Ix = mean(Ix(:,:,end-NpointsToKeep+1:end),3);
Iy = mean(Iy(:,:,end-NpointsToKeep+1:end),3);

% Display
figure(2)
clf
sidePowerLoss = ( cumsum(-Iy(1,:),2)'+cumsum(Iy(end,:),2)' ) * dx;
compensatedPower =  squeeze(sum(Ix,2)) * dy + sidePowerLoss;
plot(compensatedPower, '-', 'LineWidth',2)
line([0 numel(compensatedPower)],[max(compensatedPower) max(compensatedPower)],'Color',[0 0 0],'LineStyle','--')
ylim([0 1.1*max(compensatedPower)])
ylabel('Integrated intensity [W]')
xlabel('Main propagation axis [pixels]')
set(gca,'FontSize',18)
title(['Total power loss: ' num2str(max(compensatedPower)-compensatedPower(end),2) 'W'])

Posted 8 years ago #

Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,640
Hi Anthony,

Sorry for not coming back to you sooner.

I haven't run your code, however, if you want to calculate time average intensity, then it's important to make sure that you are using an integer number of time points per period (remember the formula for time average intensity is defined as the integral over a period). Otherwise, your time average intensity will be incorrect.

I normally use some code like this:
```
% compute points per temporal period
PPP = round(points_per_wavelength / CFL);

% compute corresponding time spacing
dt = 1/(PPP*source_freq);    

% create the time array using an integer number of points per period
t_end = 30e-6;
kgrid.t_array = 0:dt:t_end;  

% only record the final few periods when the field is in steady state
RECORD_PERIODS = 3;
sensor.record_start_index = kgrid.Nt - RECORD_PERIODS*PPP + 1;
```
Let me know if that helps.

Brad.
Posted 8 years ago #
Anthony

Member
Joined: Apr '13
Posts: 96

Hi Brad,

I usually do that also in my real codes. Here I simply adjusted the number of points to keep: as there were 16.66666667 points per wavelength with the default CFL, I simply kept 50 points to have 3 entire periods. I admit this is very unelegant ;-)

Anthony

Posted 8 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,640

Ok, no problems.

Could you walk me through exactly what you are trying to calculate?

Posted 8 years ago #
Anthony

Member
Joined: Apr '13
Posts: 96

Hi Brad,

My goal is to calculate the heat deposition in the context of HIFU treatments.

I don't know if you remember an old topic I created about heat source term calculation. I would like to avoid the plane wave approximation, so I tried to compute -div(<I>), without success. Then I moved to a compromise:
- Let z denote the main propagation axis (in 3D).
- I want to compute how much power is lost by the acoustic wave within each xy plane. To do that, I compute \int Iz dx dy and differentiate it with respect to z to get the deposited power per meter (of course, I also compensate for the power leaving the domain through the sides).
- Then, I use the plane wave approximation (or anything else) to distribute this deposited power within each xy plane.

At least, I am sure that the energy deposited within each plane exactly corresponds to the energy that the acoustic wave has lost. Then, the exact repartition within each xy plane is approximated but it should not change much my simulated lesions.

My problem is:
when I test this method with lossless propagation, there should be no power power loss and therefore no heating. However, in the 2D example I reported, there seems to be some significant power loss along the propagation (Note: the main propagation axis is noted x instead of z in my 2D example).

I am not 100% sure it comes from k-wave: maybe the way I reason or calculate is wrong and yet it seems simple in the 2D example...

Thanks for taking the time to answer!
Anthony

Posted 8 years ago #

Anthony

Member
Joined: Apr '13
Posts: 96

Sorry, I made a stupid mistake in my 2D example: everything is perfect.

clear all
close all
clc

%% PARAMETERS TO CHANGE
medium.alpha_coeff =  0; % liver: 0.37
medium.alpha_power =  0; % liver: 1.266

% create the computational grid
Nx = 128-40;
Ny = 128-40;
dx = 0.1e-3;
dy = 0.1e-3;
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500;
medium.density = 1000;
medium.BonA = 0;

% create time array
kgrid.t_array = makeTime(kgrid, medium.sound_speed);

% create initial pressure distribution using makeDisc
R = Ny*0.45;
Rap = 28/38*R;
source.p_mask = zeros(Nx, Ny);
xc = round(Nx/2);
yc = round(Ny/2);
xgrid = (1:Nx)-xc;
ygrid = ((1:Ny)-yc)';
rgrid = sqrt(bsxfun(@plus,xgrid.^2,ygrid.^2));
source.p_mask(logical(bsxfun(@times,(abs(rgrid-R)<1/2), bsxfun(@times,ygrid<0,abs(xgrid)<Rap)))) = 1;
source.p_mode = 'dirichlet';

figure(1)
clf
imagesc(source.p_mask)
axis equal

% define source
source_mag = 0.1e6;
source_f0 = 3e6;
source.p = source_mag*cos(2*pi*kgrid.t_array*source_f0);

% zero pad and filter
source.p = filterTimeSeries(kgrid, medium, [source.p, zeros(1, length(source.p))]);

% define sensor
sensor.mask = ones(Nx, Ny);
sensor.record = {'p', 'u_non_staggered'};

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,...
    'DisplayMask', 'off', 'PlotScale', [-0.5, 0.5]*source_mag,'PMLInside',false,'PlotPML',false);

% time shift
ux = timeShift(sensor_data.ux_non_staggered, kgrid.dt);
uy = timeShift(sensor_data.uy_non_staggered, kgrid.dt);

% Intensity
NpointsToKeep = 50; % in order to have an integer number of points)
Ix = reshape(sensor_data.p.*ux,[Nx,Ny,size(ux,2)]);
Iy = reshape(sensor_data.p.*uy,[Nx,Ny,size(uy,2)]);
Ix = mean(Ix(:,:,end-NpointsToKeep+1:end),3);
Iy = mean(Iy(:,:,end-NpointsToKeep+1:end),3);

% Display
figure(2)
clf
sidePowerLoss = ( cumsum(-Iy(:,1),1)+cumsum(Iy(:,end),1) ) * dx;
compensatedPower =  squeeze(sum(Ix,2)) * dy + sidePowerLoss;
plot(compensatedPower, '-', 'LineWidth',2)
line([0 numel(compensatedPower)],[max(compensatedPower) max(compensatedPower)],'Color',[0 0 0],'LineStyle','--')
ylim([0 1.1*max(compensatedPower)])
ylabel('Integrated intensity [W]')
xlabel('Main propagation axis [pixels]')
set(gca,'FontSize',18)
title(['Total power loss: ' num2str(max(compensatedPower)-compensatedPower(end),2) 'W'])

Posted 8 years ago #

Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,640

Great! Nice to see the non-staggered velocity calculations working in practice.

Posted 8 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

Analysis of the power loss along propagation in lossless medium

Reply

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