k-Wave User Forum » Forum: Ultrasound Simulations - Recent Posts
http://www.k-wave.org/forum/forum/ultrasound-simulations
Support for the k-Wave MATLAB toolboxen-USSat, 07 Sep 2024 23:05:48 +0000http://bbpress.org/?v=1.0.2<![CDATA[Search]]>q
http://www.k-wave.org/forum/search.php
bencox on "Unwanted pressure oscillation"
http://www.k-wave.org/forum/topic/unwanted-pressure-oscillation#post-9134
Wed, 28 Aug 2024 18:35:27 +0000bencox9134@http://www.k-wave.org/forum/<p>Hi a_floquet,</p>
<p>What kind of input function are you using? Have you done a convergence test, ie. reducing <code>dx</code> and <code>dt</code> until the output converges? If that doesn't help please could you post a version of your simulation highlighting the issue in the absence of a scatterer.</p>
<p>Best wishes<br />
Ben
</p>bencox on "Alpha Power Coefficient for HiFU propagating through a Heterogenous Medium"
http://www.k-wave.org/forum/topic/alpha-power-coefficient-for-hifu-propagating-through-a-heterogenous-medium#post-9132
Wed, 28 Aug 2024 18:23:53 +0000bencox9132@http://www.k-wave.org/forum/<p>Hi scmmw,</p>
<p>Yes, `medium.alpha_power' must be a scalar. It's a restriction of the method that k-Wave uses to compute the absorption term. As you're doing a HIFU simulation it's presumably at a single driving frequency, so you just have to ensure that the value of the absorption coefficient alpha = alpha_coeff * f^alpha_power is the correct value for that frequency, eg. set alpha_power to 2 and then adjust alpha_coeff until alpha is the value you need.</p>
<p>Hope that helps.</p>
<p>Best wishes<br />
Ben
</p>ArteryLee on "Simulation Unstable"
http://www.k-wave.org/forum/topic/simulation-unstable#post-9128
Wed, 21 Aug 2024 16:47:26 +0000ArteryLee9128@http://www.k-wave.org/forum/<p>Hi Bradley,</p>
<p>Thanks very much for your reply. I have tried to remove the line <code>medium.alpha_mode = 'no_dispersion'</code>. However, it still produces all NaN matrix in the sensor in this case as long as I assign alpha coefficient to the collimator mask by <code>medium.alpha_coeff(collimator_mask) = alpha_collimator</code>. I mean, once I comment out this line, it can produce non-NaN results with the C++/CUDA code.</p>
<p>BTW, I just modified my code according to your advice by setting the power law exponent to 2. In this case, it can also produce non-NaN results with the C++/CUDA code. But I am not sure how to get the correct power law absorption coefficients based on my previous numbers now. Would you mind explaining more about this?</p>
<p>Thank you!</p>
<p>Best,<br />
Artery
</p>Bradley Treeby on "Simulation Unstable"
http://www.k-wave.org/forum/topic/simulation-unstable#post-9126
Sun, 18 Aug 2024 10:36:50 +0000Bradley Treeby9126@http://www.k-wave.org/forum/<p>Hi Artery, the reason is that the no dispersion input is not supported in the C++ code. Unfortunately the docs are not super clear on this. However, if you’re doing a CW simulation, you can also get no dispersion but setting the power law exponent to 2, and then adjusting the pre factor so you get the same level of absorption at your driving frequency.
</p>ArteryLee on "Simulation Unstable"
http://www.k-wave.org/forum/topic/simulation-unstable#post-9123
Wed, 14 Aug 2024 19:39:54 +0000ArteryLee9123@http://www.k-wave.org/forum/<p>Hi Bradley,</p>
<p>I encountered a situation similar to the one Toni described above. In my simulation, I added a collimator into the medium. With <code>kspaceFirstOrder3DC</code> or <code>kspaceFirstOrder3DG</code>, I got all NaN in the sensor grid as long as I added the collimator. While with common <code>kspaceFirstOrder3D</code>, whether casting data to <code>single</code> or <code>gpuArray-single</code>, the results were normal. Meanwhile, without adding the collimator, <code>kspaceFirstOrder3DC</code> and <code>kspaceFirstOrder3DG</code> work fine with pure water medium.</p>
<p>I'm not sure how this would happen as I have used the function <code>checkStability</code> to make sure my dt (3.3333e-08) is smaller than the dt_stability_limit (5.0456e-08).</p>
<p>Thank you very much for taking a look into this!</p>
<p>Best,<br />
Artery</p>
<p>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br />
% =========================================================================<br />
% DEFINE LITERALS<br />
% =========================================================================<br />
% medium parameters<br />
c0 = 1500; % sound speed [m/s]<br />
rho0 = 1000; % density [kg/m^3]<br />
alpha0 = 3.0227e-05; % water power law absorption coefficient [dB/(MHz^y cm)]</p>
<p>% collimator parameters<br />
c_collimator = 2410;<br />
rho_collimator = 1160;<br />
alpha_collimator= 9.54;</p>
<p>% source parameters<br />
source_f0 = 5e5; % source frequency [Hz]<br />
source_roc = 59.3e-3; % bowl radius of curvature [m]<br />
source_diameter = 25.4e-3; % bowl aperture diameter [m]<br />
source_amp = 1.4*2e4; % source pressure [Pa]<br />
focus_depth = 38.1e-3; % [m]</p>
<p>% grid parameters<br />
axial_size = 70e-3; % total grid size in the axial dimension [m]<br />
lateral_size = 70e-3; % total grid size in the lateral dimension [m]</p>
<p>% computational parameters<br />
ppw = 12; % number of points per wavelength<br />
t_end = 200e-6; % total compute time [s]<br />
record_periods = 1; % number of periods to record<br />
cfl = 0.2; % CFL number<br />
source_x_offset = 20; % grid points to offset the source<br />
bli_tolerance = 0.01; % tolerance for truncation of the off-grid source points<br />
upsampling_rate = 10; % density of integration points relative to grid</p>
<p>% =========================================================================<br />
% DEFINE KGRID<br />
% =========================================================================</p>
<p>% calculate the grid spacing based on the PPW and F0<br />
dx = c0 / (ppw * source_f0); % [m]</p>
<p>% compute the size of the grid<br />
Nx = roundEven(axial_size / dx) + source_x_offset;<br />
Ny = roundEven(lateral_size / dx);<br />
Nz = Ny;</p>
<p>% create the computational grid<br />
kgrid = kWaveGrid(Nx, dx, Ny, dx, Nz, dx);</p>
<p>% compute points per temporal period<br />
PPP = round(ppw / cfl);</p>
<p>% compute corresponding time spacing<br />
dt = 1 / (PPP * source_f0);</p>
<p>% create the time array using an integer number of points per period<br />
Nt = round(t_end / dt);<br />
kgrid.setTime(Nt, dt);</p>
<p>% =========================================================================<br />
% DEFINE KWAVEARRAY FOR TRANDUCER<br />
% =========================================================================</p>
<p>% set bowl position and orientation<br />
bowl_pos = [kgrid.x_vec(1) + source_x_offset * kgrid.dx, 0, 0];<br />
focus_pos = [bowl_pos(1) + focus_depth, 0, 0];</p>
<p>% create empty kWaveArray<br />
karray = kWaveArray('BLITolerance', bli_tolerance, 'UpsamplingRate', upsampling_rate, 'SinglePrecision', true);</p>
<p>% add bowl shaped element<br />
karray.addBowlElement(bowl_pos, source_roc, source_diameter, focus_pos);</p>
<p>% =========================================================================<br />
% DEFINE SOURCE SIGNAL<br />
% =========================================================================</p>
<p>% create time varying source<br />
source_sig = createCWSignals(kgrid.t_array, source_f0, source_amp, 0);</p>
<p>% assign binary mask<br />
source.p_mask = karray.getArrayBinaryMask(kgrid);</p>
<p>% assign source signals<br />
source.p = karray.getDistributedSourceSignal(kgrid, source_sig);</p>
<p>% =========================================================================<br />
% DEFINE BINARY COLLIMATOR MASK<br />
% =========================================================================</p>
<p>[OUTPUTgrid] = VOXELISE(166,166,152,'Collimator_14mm.STL','xyz');</p>
<p>OUTPUTgrid = imrotate3(double(OUTPUTgrid), 90, [1 0 0]);<br />
OUTPUTgrid = logical(OUTPUTgrid);</p>
<p>collimator_mask = zeros(Nx, Ny, Nz, 'logical');</p>
<p>collimator_mask(ceil(Nx/2)-ceil(size(OUTPUTgrid,1)/2)-59-14:ceil(Nx/2)+ceil(size(OUTPUTgrid,1)/2)-1-59-14, ...<br />
ceil(Ny/2)-ceil(size(OUTPUTgrid,2)/2)+1:ceil(Ny/2)+ceil(size(OUTPUTgrid,2)/2)-1+1, ...<br />
ceil(Nz/2)-ceil(size(OUTPUTgrid,3)/2)+2:ceil(Nz/2)+ceil(size(OUTPUTgrid,3)/2)-1+2) = OUTPUTgrid;</p>
<p>% =========================================================================<br />
% DEFINE MEDIUM<br />
% =========================================================================</p>
<p>% assign medium properties<br />
medium.sound_speed = c0.*ones(Nx, Ny, Nz);<br />
medium.sound_speed(collimator_mask) = c_collimator;<br />
medium.density = rho0.*ones(Nx, Ny, Nz);<br />
medium.density(collimator_mask) = rho_collimator;<br />
medium.alpha_coeff = alpha0.*ones(Nx, Ny, Nz);<br />
medium.alpha_coeff(collimator_mask) = alpha_collimator;</p>
<p>medium.alpha_power = 1.001;<br />
medium.alpha_mode = 'no_dispersion';</p>
<p>% =========================================================================<br />
% DEFINE SENSOR<br />
% =========================================================================</p>
<p>% set sensor mask to record central plane, not including the source point<br />
sensor.mask = zeros(Nx, Ny, Nz);<br />
sensor.mask(source_x_offset + 2:end, :, Nz/2 + 1) = 1; % mid-plane/central plane</p>
<p>% record the pressure<br />
sensor.record = {'p', 'p_rms'};</p>
<p>% record only the final few periods when the field is in steady state<br />
sensor.record_start_index = kgrid.Nt - record_periods * PPP + 1;</p>
<p>% =========================================================================<br />
% SIMULATION<br />
% =========================================================================</p>
<p>% set input options<br />
input_args = {...<br />
'PMLSize', 'auto', ...<br />
'PMLInside', false, ...<br />
'PlotPML', false, ...<br />
'DisplayMask', 'off'};</p>
<p>% run code<br />
sensor_data = kspaceFirstOrder3DG(kgrid, medium, source, sensor, input_args{:});<br />
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
</p>healsea on "Embeddings with larger sound speed and density"
http://www.k-wave.org/forum/topic/embeddings-with-larger-sound-speed-and-density#post-9121
Thu, 01 Aug 2024 01:27:27 +0000healsea9121@http://www.k-wave.org/forum/<p>Hello,</p>
<p>I am using k-Wave to simulate some embeddings within tissue and have been referencing your example, example_us_bmode_phased_array.m. I encountered an issue when adjusting the wave speed and density of the embedding.</p>
<p>When the embedding's wave speed and density are half that of the tissue, the simulation works perfectly. However, when I set the wave speed and density to twice that of the tissue, the sensor data becomes unstable, increasing exponentially and eventually resulting in NaN values. Consequently, the B-mode image is not visible.</p>
<p>I have attached my code defining the properties for your reference. Thank you very much for your software and assistance!</p>
<p>Best wishes<br />
Haohui Zhang</p>
<p>'% define properties<br />
sound_speed_map = c0 * ones(Nx, Ny, Nz) ;<br />
density_map = rho0 * ones(Nx, Ny, Nz) ;<br />
sound_speed_iron = c0*2;<br />
density_iron = rho0*2;</p>
<p>% define a cylinder for Iron oxide<br />
radius = 4e-3;<br />
x_pos = 28e-3;<br />
y_pos = dy * Ny/2;<br />
x_loc = round(x_pos / dx);<br />
height = 4e-3;<br />
x_height = round(height / dx);<br />
neighbor = 12e-3;</p>
<p>scattering_plane1 = makeDisc(Ny, Nz, round(y_pos / dx), Nz/2, round(radius / dx)); % Center array<br />
scattering_plane2 = makeDisc(Ny, Nz, round((y_pos-neighbor)/ dx), Nz/2, round(radius / dx)); % Left array<br />
scattering_plane3 = makeDisc(Ny, Nz, round((y_pos+neighbor)/ dx), Nz/2, round(radius / dx)); % Right array<br />
scattering_plane = scattering_plane1 + scattering_plane2 + scattering_plane3;</p>
<p>scattering_region = repmat(scattering_plane,[1,1,Nx]);<br />
scattering_region = permute(scattering_region,[3,1,2]);<br />
tmp = scattering_region;<br />
tmp([1:x_loc,x_loc+x_height:end],:,:) = 0;<br />
% assign region<br />
sound_speed_map(tmp == 1) = sound_speed_iron;<br />
density_map(tmp == 1) = density_iron;</p>
<p>% assign to the medium inputs<br />
medium.sound_speed = sound_speed_map;<br />
medium.density = density_map;'
</p>bencox on "Cylindrical transducer array"
http://www.k-wave.org/forum/topic/cylindrical-transducer-array#post-9120
Wed, 31 Jul 2024 10:52:22 +0000bencox9120@http://www.k-wave.org/forum/<p>Hi p_acoustics,</p>
<p>Interesting question. There's not an obvious way to do it. </p>
<p>One way to make a source that is directional is to put a monopole source (<code>source.p</code>) and a dipole source (<code>source.ux</code>, ...) at the same source point, ie. <code>source.p_mask = source.u_mask</code>. If the scaling is right (pressure and velocity related by the acoustic impedance), the waves propagating in one direction will destructively interfere and in another constructively interfere. Also note that to have the sources at exactly the same point it will be necessary to turn off the staggered grid, setting the <code>UseSG</code> option to false. Also note that the pressure source is advanced in time compared to the velocity source by half-a-timestep, so <code>source.p</code> may have to be shifted back by that much.</p>
<p>I'd be interested to know if it works!</p>
<p>Best wishes<br />
Ben
</p>Bradley Treeby on "Viscoelastic simulation on kWave"
http://www.k-wave.org/forum/topic/viscoelastic-simulation-on-kwave#post-9118
Tue, 30 Jul 2024 08:14:31 +0000Bradley Treeby9118@http://www.k-wave.org/forum/<p>For viscoelastic simulations using <code>pstdElastic2D</code> / <code>pstdElastic3D</code> it is not possible to set the power law. The Kelvin-Voigt model is used to model the loss, which gives an power law dependence of 2 for small attenuation.
</p>masud407 on "Viscoelastic simulation on kWave"
http://www.k-wave.org/forum/topic/viscoelastic-simulation-on-kwave#post-9113
Tue, 16 Jul 2024 16:33:14 +0000masud4079113@http://www.k-wave.org/forum/<p>Hello,</p>
<p>If I want to simulate the viscoelastic medium, I just need to define α0 ≡ medium.alpha_coeff and y ≡ medium.alpha_power which ranges from 0 to 3. Is there any specific things I should do?</p>
<p>Thanks
</p>scmmw on "Alpha Power Coefficient for HiFU propagating through a Heterogenous Medium"
http://www.k-wave.org/forum/topic/alpha-power-coefficient-for-hifu-propagating-through-a-heterogenous-medium#post-9109
Thu, 11 Jul 2024 15:51:43 +0000scmmw9109@http://www.k-wave.org/forum/<p>Hi, I am trying to model HiFU through a heterogeneous medium but I am struggling with how best to approach assigning the alpha power coefficient as I am propagating the ultrasound through three different regions: water, fat and tissue. </p>
<p>Each region has a different alpha power value based on experimental data values. However, kWave only allows for a scalar value of alpha_power rather than a matrix. I discovered on the forum the function fitPowerLawParamsMulti, but this still requires picking a reference value y_ref.</p>
<p>Currently, the three values for alpha power are the following,</p>
<p> % define different absorption properties for different regions<br />
medium.alpha_power = 1 * ones(Nx, Ny); % [dB/(MHz^y cm)] - Water<br />
medium.alpha_power(Region_2{j}) = 0.4; % [dB/(MHz^y cm)] Region 2 - Fat<br />
medium.alpha_power(Region_3{j}) = 1.41; % [dB/(MHz^y cm)] Region 3 - Tissue</p>
<p>When I use the function, it plots a graph, and based on the point where the lines overlap for each three regions, I changed the y_ref to be a value in between the points of overlap for Region 2 and 3, which is roughly 0.55. Is this the correct approach to define a y_ref, or am I misunderstanding how best to define y_ref? Any advice on how to go about this would be much appreciated!
</p>Tomaubier on "Derivation of the transducer surface pressure from the total radiated power"
http://www.k-wave.org/forum/topic/derivation-of-the-transducer-surface-pressure-from-the-total-radiated-power#post-9105
Tue, 25 Jun 2024 09:08:16 +0000Tomaubier9105@http://www.k-wave.org/forum/<p>Hello,<br />
I'm currently performing simulations of a spherical mono element transducer radiating in a multilayered axisymmetric domain. My goal is to evaluate temperature elevations produced by this transducer in the presence of those layers. The input of these simulations is the total acoustic power emitted by the transducer measured in water using a radiation force balance.<br />
To be able to run the simulation I have to go from this acoustic power <code>tx_emitted_ac_power</code> to pressure magnitudes <code>tx_p_mag</code> that can be applied to the transducer surface.<br />
Considering the plane wave assumption, a pressure magnitude can be obtained such that<br />
<code>tx_p_mag = sqrt(2) * sqrt((tx_emitted_ac_power * rho_tx_coupling_medium * c_tx_coupling_medium) / tx_surface_area)</code>.</p>
<p>Regarding the transducer surface, my first instinct was to evaluate it as a section of a sphere such that<br />
<code>tx_surface_area = 2pi * r (r- sqrt(r^2 - a^2))</code> where <code>r</code> corresponds to the radius of curvature and <code>a</code> to half of the transducer aperture.</p>
<p>After several validation with Rayleigh-based simulation models in homogeneous domains, I realized that the pressures obtained with kWave were consistently under-estimated.</p>
<p>I later realized that evaluating <code>tx_surface_area</code> as the planar projection of the transducer aperture tx_surface_area = pi * (tx_aperture / 2)^2 led to an agreement of kWave with the reference Rayleigh simulations.</p>
<p>Does anyone know what is going on here?</p>
<p>Best,</p>
<p>Tom
</p>cyl on "Transcranial ultrasound attenuation"
http://www.k-wave.org/forum/topic/transcranial-ultrasound-attenuation#post-9102
Tue, 04 Jun 2024 13:53:29 +0000cyl9102@http://www.k-wave.org/forum/<p>%%<br />
% 定义网格<br />
Nx = 300; % number of grid points in the x direction<br />
Ny = 300; %number of grid points in the y direction<br />
Nz = 300; % number of grid points in the z direction<br />
dx =0.1e-3; % grid point spacing in the x direction [m]<br />
dy =0.2e-3; % grid point spacing in the y direction [m]<br />
dz =0.1e-3; % grid point spacing in the z direction [m]<br />
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);<br />
%%<br />
% define the properties of the propagation medium<br />
medium.sound_speed = 1540; % [m/s]<br />
medium.density = 1000; % [kg/m^3]<br />
medium.sound_speed = medium.sound_speed * ones(Nx, Ny, Nz); % [m/s]<br />
medium.sound_speed(:,75:80,:) = 3500 ;%填颅骨的宽度和声速<br />
medium.density = medium.density * ones(Nx, Ny, Nz); % [kg/m^3]<br />
medium.density(:,75:80,:) = 2000 ;%填颅骨的宽度和密度<br />
medium.alpha_coeff=0.05*ones(Nx,Ny,Nz); %[dB/(MHz^y cm)]<br />
medium.alpha_coeff(:,75:80,:) = 2.7; %颅骨衰减系数<br />
%%<br />
% tempo e step de simulação<br />
kgrid.makeTime(medium.sound_speed);<br />
%%<br />
% target ping signal<br />
source_strength = 1e3; % [Pa]<br />
tone_burst_freq = 1e6; % [Hz]<br />
tone_burst_cycles = 1;<br />
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);<br />
%%<br />
%Setting up ultrasonic transducers<br />
Transducer.number_elements=1; %阵元个数<br />
Transducer.element_width = 5e-3; %阵元宽度<br />
Transducer.element_length = 12 ;%阵元高度<br />
Transducer.element_spaceing = 0 ;%阵元间距<br />
Transducer.radius = inf;%换能器曲率半径<br />
%计算换能器总宽度<br />
transducer_width = Transducer.number_elements*Transducer.element_width;<br />
transducer_x_pos = Nx/23;%grid points<br />
transducer_y_pos = 25;%grid points<br />
transducer_z_pos = 170;%grid points<br />
Transducer.positoin = round([transducer_x_pos,transducer_y_pos,transducer_z_pos]);<br />
% 用于波束合成延迟叠加的参数<br />
Transducer.sound_speed = 1500; % 声速 [m/s]<br />
Transducer.focus_distance = 20e-3; % 焦距 [m]<br />
Transducer.elevation_focus_distance = 19e-3; % elevation plane 仰角平面焦距 [m]<br />
Transducer.steering_angle = 0; % steering angle 偏转角 [degrees]<br />
% 设置变迹 apodization<br />
Transducer.transmit_apodization = 'Rectangular';<br />
Transducer.receive_apodization = 'Rectangular';<br />
% 设置当前激活的换能器阵元 define the transducer elements that are currently active<br />
Transducer.active_elements = 1;<br />
%设置用于驱动换能器的输入信号<br />
Transducer.input_signal = input_signal;<br />
%创建换能器<br />
transducer = Transducer(kgrid, Transducer);<br />
%%<br />
%创建具有检测传感器<br />
sensor.mask = zeros(Nx,Ny,Nz);<br />
sensor.mask(Nx/2,275,170)=1;<br />
%run the simulation<br />
[sensor_data] = kspaceFirstOrder3D(kgrid , medium , Transducer, sensor,input_args{:});<br />
%计算输入信号和每个传感器位置记录信号的幅度<br />
[f_input, as_input] = spect([input_signal,zeros(1,2*length(input_signal))],1/kgrid.dt);<br />
[~,as_1]=spect(sensor_data(1,:),1/kgrid.dt);<br />
[sensor_data] = kspaceFirstOrder3D(kgrid, medium, transducer, sensor, input_args{:});
</p>cyl on "Transcranial ultrasound attenuation"
http://www.k-wave.org/forum/topic/transcranial-ultrasound-attenuation#post-9100
Tue, 04 Jun 2024 13:50:02 +0000cyl9100@http://www.k-wave.org/forum/<p>Hello, I am currently researching the use of a single element transducer to generate ultrasound penetration through the skull, and then testing the sound field beneath the skull to obtain the attenuation coefficient of the skull. Here is my code. Currently, I have encountered an issue: during simulation, it will display‘transducer=Transducer (kgrid, Transducer)；there is an error Cannot use a value of type kWaveGrid as an index’ , which I have researched a lot of information but still cannot solve. I hope to receive your help. Thank you. Additionally, if there are any issues with my code, I hope you can point them out. Thank you very much again
</p>cyl on "Defining An Ultrasound Transducer Example Question"
http://www.k-wave.org/forum/topic/defining-an-ultrasound-transducer-example-question#post-9099
Tue, 04 Jun 2024 13:31:14 +0000cyl9099@http://www.k-wave.org/forum/<p>Hello, I would like to know how your problem was resolved
</p>karanmanoj on "Underwater acoustic simulation in shallow waters with bathymetry"
http://www.k-wave.org/forum/topic/underwater-acoustic-simulation-in-shallow-waters-with-bathymetry#post-9095
Tue, 21 May 2024 18:01:44 +0000karanmanoj9095@http://www.k-wave.org/forum/<p>Hi, I am trying to use the k-Wave toolbox for my master's thesis on Underwater Acoustics. I am interested in simulating the wave propagation of acoustic waves, up to about 50 - 100 kHz in shallow waters. The bathymetry (seafloor) of shallow waters acts as a frequency-dependent absorber. I do not know how to model this condition in the k-Wave toolbox. I am not interested in modelling the wave propagation into the subsurface, I just want to apply some kind of absorption/impedance condition on a surface (with varying depth) in a 3D domain. Can somebody please help me on this? Could not find anything in the documentation or in this forum.</p>
<p>Thanks,<br />
Karan
</p>g.c. on "A Question about The Reflection at The Hetrogeneous Media Boundary"
http://www.k-wave.org/forum/topic/a-question-about-the-reflection-at-the-hetrogeneous-media-boundary#post-9089
Tue, 07 May 2024 15:09:01 +0000g.c.9089@http://www.k-wave.org/forum/<p>Hi!<br />
I am currently working on my master thesis and want to simulate the propagation of a pressure distribution in a structure of alterning Al and Pb foils (a few um thick, circular, 2cm radius, the whole thing is in vacuum). I am having problems with my result, as it is diverging: the mean pressure on a thin detector on the other side of the stack is going towards 10^293Pa. I assumed it accured because my geometry is not smooth enougth, so to check I just built a block of Pb in vacuum and put a circular pressure source inside of it, and run the simulation in 2D. As I understood from one of the last posts (point 4) it shouldn't be a problem if the source originates from within the high impedance material. Why does my code still diverge? is there a solution to this problem? i tried to change the material properties, but I have to change them so much, that the simulation is not rapresentative anymore. </p>
<p>Here is my code:<br />
<pre><code>%% % create masks
z_lim = 0.03;
x_lim = 0.01;%in m
Nx = 200; % number of grid points in the x direction
Nz = 600; % number of grid points in the z direction
dz = 0.00005;
dx = 0.0001;
[Z,X] = meshgrid((0:dz:z_lim-dz),(-1*x_lim:dx:x_lim-dx));
geom = zeros(size(Z));
det = zeros(size(Z));
lz = 0.02; % [m], z length
lx = 0.01; % [m], x length
z0 = 0.015; % [m], z center
x0 = 0; % [m], x center
ind = ((Z >= z0-lz/2) & (Z < z0+lz/2) & (X >= x0-lx/2) & (X < x0+lx/2));
[num, ~] = materials(4); %add the material of interest
geom(ind) = num;
geom = permute(geom,[2,1]);
% Create the computational grid for the 2D slice
kgrid = kWaveGrid(Nz, dz, Nx, dx);
[densitymask,speedmask,alphamask,ymask] = get_masks(geom);
%%
% Define the properties of the upper layer of the propagation medium
medium.sound_speed = speedmask;
medium.density = densitymask;
% Create time array
t_end = 3e-6; % [s]
kgrid.makeTime(medium.sound_speed, [], t_end);
%%
%create SOURCE
% define a single source point
disc_magnitude = 300; % [Pa]
disc_x_pos = Nx/2; % [grid points]
disc_z_pos = Nz/4;
disc_radius = 0.4; % [grid points]
disc = disc_magnitude * makeDisc(Nz, Nx, disc_z_pos, disc_x_pos, disc_radius);
source.p0 = zeros(Nz, Nx);
source.p0 = disc ;
% %%
% %SENSOR
det = makeLine(Nz, Nx, [round(0.023/dz), 90],[round(0.023/dz), 110]);
%det = permute(det,[2,1]);
sensor.mask=det;
%%
% Run the simulation with PlotLayout and PlotScale
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,...
'PlotLayout', true);
%%
% Assuming your pressure data is stored in a variable named 'pressure_data'
figure
% Calculate the mean pressure across all pixels at each time step
mean_pressure = mean(sensor_data, 1);
% Plot pressure versus time
plot(kgrid.t_array, mean_pressure, 'LineWidth', 2);
xlabel('Time');
ylabel('Mean Pressure on detector');
title('Mean Pressure vs Time');
grid on;</code></pre>
<p>the materials are defined as follows: </p>
<pre><code>function [num, mat, rho, c, gamma, alpha_coeff,y] = materials(arg)
% attenuation data: Appl. Sci. 2020, 10(7), 2230; <a href="https://doi.org/10.3390/app10072230" rel="nofollow">https://doi.org/10.3390/app10072230</a>
if arg == 0 | strcmpi(arg,'vacuum') | strcmpi(arg,'vac')
mat = 'vacuum';
num = 0;
rho = 0.0012 ;% 1.2e-3; % [g/cm^3] density
c = 0.4; %0.4; % [mm/us] speed of sound
beta = 3400e-6; % [1/K] vol. thermal expansion
c_p = 1012; % [J/kg/K] spec. heat capacity
gamma = beta*(c*1000)^2/c_p; % [Pa/(J/m^3)] Grüneisenparameter
alpha_coeff = 1e-6; %0 [dB/cm/MHz^y] acoustic attenuation coef.
y = 1; % attenuation exponent
elseif arg == 4 | strcmpi(arg,'lead') | strcmpi(arg,'Pb')
mat = 'Pb';
num = 4;
A = 207; % [u] eff. Massezahl
Z = 82; % [u] eff. Ordnungszahl
rho = 11.29;%11.29; % [g/cm^3] density
c = 1.96; % 1.96 [mm/us] speed of sound
beta = 29e-6; % 87e-6; % [1/K] vol. thermal expansion
c_p = 129; % [J/kg/K] spec. heat capacity
gamma = beta*(c*1000)^2/c_p; % [Pa/ (J/m^3)] Grüneisenparameter
I = 823; % [eV] mean excitation energy
alpha_coeff = 4.33; %alpha = 4.33; % [dB/cm/MHz] acoustic attenuation coef.
y = 1; % attenuation exponent
end</code></pre>
<p>and get masks is a function: </p>
<pre><code>function [densitymask,speedmask,alphamask,ymask] = get_masks(geom)
% create mask arrays
densitymask = zeros(size(geom));
speedmask = zeros(size(geom));
alphamask = zeros(size(geom));
ymask = zeros(size(geom));
% initialise counter
n = 0;
while ~all(densitymask(:))
% assign values of materials to masks
[~, ~, rho, c, ~, alpha_coeff,y] = materials(n);
densitymask(geom == n) = rho*1e3; % [kg/m^3]
speedmask(geom == n) = c*1e3; % [m/s]
alphamask(geom == n) = alpha_coeff; % [dB/cm/MHz^y]
ymask(geom == n) = y; %
% increase counter
n = n+1;
end
end</code></pre>
<p>Thank you in advance for any help.<br />
Best regards,<br />
Gianna
</p>a_floquet on "Unwanted pressure oscillation"
http://www.k-wave.org/forum/topic/unwanted-pressure-oscillation#post-9076
Mon, 15 Apr 2024 09:01:05 +0000a_floquet9076@http://www.k-wave.org/forum/<p>Hello,</p>
<p>I am using k-wave to simulate ultrasound echography. Among other things, I wish to simulate a PSF.<br />
I work in 2D, so I don't use the transducer class, I instead "manually" place my elements and delay to create my probe. However, did test in 3D using the transducer class and my problem remained the same.</p>
<p>My problem is that they are unwanted pressure oscillation in the received signal. My sensor is, of course, at the same position as my emiiting elements. These oscillations happens even when propagating a wave in an perfectly homogeneous medium, without any scatterer at all.</p>
<p>If my understanding is correct, this is due to the finite length FFT used to compute spatial gradient, as described in section 2.8 "Smoothing and the Band-Limited Interpolant" of the user manual. Whatever option I use for the 'smooth' optionnal input parameter, this problem persists. </p>
<p>My question is then : is there any other way to mitigate this issue ? I assume that with enough smoothing I could get rid of the oscillation, but my results would not be accurate anymore ?</p>
<p>Thanks in advance
</p>p_acoustics on "Cylindrical transducer array"
http://www.k-wave.org/forum/topic/cylindrical-transducer-array#post-9075
Fri, 12 Apr 2024 01:03:16 +0000p_acoustics9075@http://www.k-wave.org/forum/<p>I saw that maybe I did not write it clearly:</p>
<p>If I have an array of rectangular sources arranged in a way (with some angle to create a convex shape) that they shape an imaginary cylinder together, and the purpose is to create sound waves outside of this cylinder, like a large cylindreical source made of mosaic of transducers.</p>
<p>How could we stop the acoustic sources from emitting the sound waves inward. Because, if the sound propagates inward (toward the inside of the cylinder), then the sound waves will eventually reach the other side of the cylinder and pass it and they will interfere with the outward propagation. I am interested in both far-field and near-field results.</p>
<p>Also, how could we create a baffle between the square sources.</p>
<p>many thanks
</p>p_acoustics on "Cylindrical transducer array"
http://www.k-wave.org/forum/topic/cylindrical-transducer-array#post-9074
Fri, 12 Apr 2024 01:03:11 +0000p_acoustics9074@http://www.k-wave.org/forum/<br />p_acoustics on "Cylindrical transducer array"
http://www.k-wave.org/forum/topic/cylindrical-transducer-array#post-9073
Thu, 11 Apr 2024 18:56:34 +0000p_acoustics9073@http://www.k-wave.org/forum/<p>Hi,</p>
<p>Thanks for the great software. k-wave is not only a great software, the research that comes with it and the papers published using k-wave really helps the users to understand the physics of acoustics.</p>
<p>I have a question, if you have a cylinder of sources (multiple sources are arranged in a way that they shape an imaginary cylinder), how could we stop the acoustic sources to emit sound from their face in-ward and not out-ward. Because if the sound propagates in-ward (toward the inside of the cylinder), then it will eventually interfere with the outward propagation. I am interested in both far field and nearfield results.</p>
<p>I looked into the research on bowl transducer (multiple sources are arranged to shape a bowl) but I did not see anything related to acoustic emission from the back of the bowl).</p>
<p><a href="http://bug.medphys.ucl.ac.uk/papers/2016-Martin-IEEETUFFC.pdf" rel="nofollow">http://bug.medphys.ucl.ac.uk/papers/2016-Martin-IEEETUFFC.pdf</a></p>
<p>I tried putting hard material (high speed of sound and high density) inside the cylinder but it did not work. I also tried to increase the absorbsion inside the cylinder and was not successful. There is no PML that we could put inside the cylinder either. I also tried to go with the half of the cylinder, but the interference from inside of the cylinder eventually will happen.</p>
<p>I appreciate your help.</p>
<p>many thanks </p>
<p>p.s: I use k-wave for python.
</p>