Hi,
I am new user of the this software, what I am want to do is to use this toolbox to simulate very strong non-linear wave propagation in water or in a layered medium, but I didn't see any discussions in this forum and didn't find any example that shows how to simulate the strong nonlinear wave propagation, is this tool box applicable for the strong non-linear wave case? can anyone give me some examples or told me what should I take care of in the simulation?
regards
Alex
Hi Alex,
From my personal experience (I am a simple user, not a dev), I would say that you have two main points to check:
- your grid is fine enough to enable the propagation of the high harmonics
- the medium heterogeity: see some topics on this forum and paragraph V.B of Brad's paper in JASA 2012 (Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method).
I am not sure of what you call "strong" nonlinear case, I made a few simulations in water with B/A = 5 at quite high power (100W) at it seems to compare well with hydrophone measurements with 50 MPa peak positive pressure. However, I only have very old experimental data I did not acquire myself so it is not rigorous at all: don't consider this as a validation! See "Modelling nonlinear ultrasound propagation in absorbing media using the k-Wave toolbox: experimental validation" (2012) for a validation but the peak pressures seems quite low (a few MPa figure 2).
By the way, use the C code instead of the Matlab version (it's much faster).
Good luck,
Anthony
Hi Anthony,
I am trying to model a shock wave as it is used in shock wave lithotripsy. You mentioned 50 MPa, this sounds like you simulated something similar? Do you think you could send me the code you used for the simulation? It would help me a lot and I would really appreciate it! I have played with a lot of parameters (power, medium.BonA, grid spacing, wavelength etc.) but I have never managed to see the formation of a shock front. It does not matter if I define the medium.BonA value or not, the result is the same. Do you have a suggestion for the grid point spacing at a wavelength of 1.5e-3 m? Does it matter if I do the simulation in 2D or 3D?
Thanks in advance,
Timo
Hi Timo,
The nonlinearity depends on the value of the pressure. In principle, you could see strong nonlinearity in either 2D or 3D. However, if you are using a focused transducer, the focusing gain in 2D (cylindrical focusing) is much weaker than in 3D (spherical focusing), so you will need much higher transmit pressures. The grid point spacing should be at least half the wavelength at the highest frequency harmonic. To assess this, you should do a convergence analysis (see here, Fig. 3 for an example).
Hope that helps.
Brad.
Hi,
I am playing with kWave to see whether I can use it to estimate the effect of HIFU on tissues when used for ablation. I have a couple of questions regarding the transducer:
- I have a single transducer with the shape of a spherical cap with radius of curvature 80 mm, radius of the base 25 mm and height 4 mm. I used the function makeSphericalSection(radius, height, width) but the width can be at maximum equal to half of the radius of curvature. Is there another way to create it with the correct dimensions?
- These transducers are used at a central frequency of the order of MHz, however I see that the maximum supported frequency is 750 kHz, did I understand correctly?
Thanks in advance,
Silvia
The source frequency is of the order of MHz.
However in the log messages I see that the maximum supported frequency is 750kHz.
Hi again,
I found all these info in manual/articles. Sorry to have bothered!
Silvia
Hi Silvia,
Glad you got it sorted. Let us know if you get stuck.
Brad.
hello dear silvia
for my project i need to acquire the pressure due to a HIFU transducer in a heterogeneous soft tissue. i have a couple of questions. may you help me? thanks in advance
hi tahere,
can you help me ?are you a chinese,i am doing the same some like yours
hi jackYANG,
I have omes problems with the High intensity focused ultrasound simualtion, can you give some help to me?I'm a Chinese too.
hi everyone，
I'm doing simulation about a HIFU transducer, but I found that when I set the grids' size smaller, the pressure of focus doesn't converge. But when I set the medium linear(delete medium.BonA=7), the results show convergency. Can anyone help me to solve this problem?
Here are my codes:
Nx = 380; % number of grid points in the x direction
Ny = 240; % number of grid points in the y direction
Nz = 240; % number of grid points in the z direction
dx =0.3e-3; % grid point spacing in the x direction [m]
dy =0.3e-3; % grid point spacing in the y direction [m]
dz =0.3e-3; % grid point spacing in the z direction [m]
x_size=dx*Nx;
y_size=dy*Ny;
z_size=dz*Nz;
sound_speed=1500; % [m/s]媒介声速
density = 1000 ; % [kg/m^3]媒介密度
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
% define the properties of the propagation medium
medium.BonA=7; %水里的超声非线性参数
medium.sound_speed = sound_speed * ones(Nx, Ny, Nz); % [m/s]
medium.density = density * ones(Nx, Ny, Nz); % [kg/m^3]
kgrid.makeTime(medium.sound_speed,0.3,80e-6);
% define a centered circular sensor
% define a series of Cartesian points to collect the data
sensor.mask = zeros(Nx, Ny, Nz);
sensor.mask(:, :, Nz/2) = 1;
sensor.record = {'p_max'};
% % input arguments
input_args = {'PlotPML', true, ...
'DataCast', 'single', 'CartInterp', 'linear','PlotSim',true,'PlotLayout',true,'PMLInside',true,'PMLSize',[20,20,20]};
% define a time varying sinusoidal source
source_mag = 543.36e3; % [Pa]
% create empty array
karray = kWaveArray('SinglePrecision',true);
%position = [1, Ny/2, Nz/2];
position = [-x_size/2+20e-3,0,0];
% define arc properties
radius = 63.2e-3; % [m]
diameter = 64e-3; % [m]
focus_pos = [x_size/2, 0, 0];
karray.addBowlElement(position, radius, diameter, focus_pos);
source.p_mask = karray.getArrayBinaryMask(kgrid);
voxelPlot(double(source.p_mask));
source_freq = 1e6;
sig1 = source_mag * sin(2 * pi * source_freq * kgrid.t_array);
% smooth the source
sig1 = filterTimeSeries(kgrid, medium, sig1,'PlotSignals',true,'PlotSpectrums',true);
source_signal = zeros(1,length(sig1));
source_signal(1,1:length(sig1)) = sig1;
source.p = karray.getDistributedSourceSignal(kgrid, source_signal);
sensor_data = kspaceFirstOrder3DG(kgrid, medium, source, sensor, input_args{:});
sensor_data.p_max = reshape(sensor_data.p_max, [Nx, Ny]);
sensor_data.p_max = reshape(sensor_data.p_max, [Nx, Ny]);
figure;
imagesc(kgrid.x_vec, kgrid.y_vec,sensor_data.p_max);
xlabel('x [m]');
ylabel('y [m]');
%zlabel('z [m]');
title('Maximum Pressure');
c = colorbar;
c.Label.String = 'pressure()';
max = sensor_data.p_max(1,1);
maxi = [1,1];
[max, maxi];
for a=floor(0.3*Nx):Nx
for b=1:Ny
if sensor_data.p_max(a,b)>max
max= sensor_data.p_max(a,b);
maxi=[a,b];
end
end
end
focus_x = 0;
focusx_max = 0;
focusy_max = 0;
focus_y = 0;
focus_l = 0;
focus_r = 0;
focus_u = 0;
focus_d = 0;
[focus_x,focus_y];
for a=floor(0.3*Nx):(Nx - 1)
for b=1:(Ny - 1)
if sensor_data.p_max(a,b)<0.71*max&&sensor_data.p_max(a,b+1)>0.71*max
%sensor_data.p_max(a,b)=0;
focus_l = b+1;
end
if sensor_data.p_max(a,b)>0.71*max&&sensor_data.p_max(a,b+1)<0.71*max
%sensor_data.p_max(a,b)=0;
focus_r = b+1;
end
if sensor_data.p_max(a,b)<0.71*max&&sensor_data.p_max(a+1,b)>0.71*max
focus_u = a+1;
end
if sensor_data.p_max(a,b)>0.71*max&&sensor_data.p_max(a+1,b)<0.71*max
%sensor_data.p_max(a,b)=0;
focus_d = a+1;
end
end
focus_x = focus_r - focus_l;
focus_y = focus_d - focus_u;
if focus_x > focusx_max
focusx_max = focus_x;
end
if focus_y > focusy_max
focusy_max = focus_y;
end
end
%ONeil solution
% define transducer parameters
radius = 63.2e-3; % [m]
diameter = 64e-3; % [m]
frequency = 1e6; % [Hz]
sound_speed = 1500; % [m/s]
density = 1000; % [kg/m^3]
Z = density*sound_speed;
p = 543.36e3;
velocity = p/Z; % [m/s]
% define position vectors
axial_position = 0:5e-4:150e-3; % [m]
lateral_position = -15e-3:5e-4:15e-3; % [m]
% evaluate pressure
[p_axial, p_lateral] = focusedBowlONeil(radius, diameter, velocity, ...
frequency, sound_speed, density, axial_position, lateral_position);
% plot
figure;
subplot(4, 1, 1);
plot(axial_position .* 1e3, p_axial .* 1e-6, 'k-');
xlabel('Axial Position [mm]');
ylabel('Pressure [MPa]');
subplot(4, 1, 2);
plot((0:dx:(Nx-1)*dx) .* 1e3, sensor_data.p_max(:,Ny/2).* 1e-6, 'k-');
xlabel('Axial Position [mm]');
ylabel('Pressure [MPa]');
subplot(4, 1, 3);
plot(lateral_position .* 1e3, p_lateral .* 1e-6, 'k-');
xlabel('Lateral Position [mm]');
ylabel('Pressure [MPa]');
subplot(4, 1, 4);
plot((0:dy:(Ny-1)*dy).* sensor_data.p_max(Ny/2,:) .* 1e-6, 'k-');
xlabel('Lateral Position [mm]');
ylabel('Pressure [MPa]');
figure;
imagesc(kgrid.x_vec, kgrid.y_vec,sensor_data.p_max);
xlabel('x [m]');
ylabel('y [m]');
title('Maximum Pressure');
c = colorbar;
c.Label.String = 'pressure()';
You must log in to post.