k-Wave User Forum » Topic: A Question about The Reflection at The Hetrogeneous Media Boundary

k-Wave User Forum » Topic: 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 Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 00:12:03 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php 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 +0000 g.c. 9089@http://www.k-wave.org/forum/ Hi! 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. Here is my code: <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> the materials are defined as follows: <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> and get masks is a function: <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> Thank you in advance for any help. Best regards, Gianna Bradley Treeby 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-5690 Thu, 29 Sep 2016 13:40:09 +0000 Bradley Treeby 5690@http://www.k-wave.org/forum/ Hi peet2521, We don't have the data presented in exactly that form, but there are some plots in <a href="http://www.k-wave.org/papers/2014-Robertson-IEEEIUS.pdf">this paper</a> that might be useful. Brad. peet2521 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-5684 Mon, 19 Sep 2016 21:41:43 +0000 peet2521 5684@http://www.k-wave.org/forum/ Brad, Are there any case studies or something on the grid resolution we should use to get a correct transmission/reflection for a given impedance mismatch ratio? cge 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-4643 Wed, 16 Jul 2014 06:52:27 +0000 cge 4643@http://www.k-wave.org/forum/ Hi Brad Thank you very much. Now it works and I have got the correct result Bradley Treeby 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-4641 Tue, 15 Jul 2014 11:21:35 +0000 Bradley Treeby 4641@http://www.k-wave.org/forum/ Hi cge, There are a few problems with your simulation: (1) Your source is defined within the perfectly matched layer (PML). The PML is a thin layer that surrounds the domain and absorbs the outgoing waves to simulate free-field conditions. The default PML thickness for 2D simulations is 20 grid points. You can either move your source to be outside the PML, or set the PML to be outside the domain you specify by using the optional input <code>'PMLInside', false</code>. (2) You are using a initial pressure source at the very edge of the grid. By default, the <code>source.p0</code> is smoothed inside the simulation functions. This means that parts of the source will be wrapped to the other side of the domain. To see what I mean, try running the following code: <pre><code>% define grid kgrid = makeGrid(20, 1); % define source and smooth p0 = [1; zeros(19, 1)]; p0_smooth = smooth(kgrid, p0, true); % plot figure; stem(p0, 'bx'); hold on; stem(p0_smooth, 'r');</code></pre> (3) The medium interfaces are defined within the PML (see discussion above). (4) As Anthony has already mentioned, the third layer has a huge impedance mismatch changing from 3.2 MRayls to 0.4 kRayl (a factor of nearly 10,000). The problem we have found with such large changes is that the transmission coefficient when waves travel from the low to high impedance material is vastly over-estimated (except using high numbers of points per wavelength), resulting in very large amplitude waves. However, in your case, this shouldn't be a problem as the source originates from within the high-impedance material. The CFL number also affects the accuracy of the simulation as mentioned above. A revised script is below. Brad. <pre><code>clear all; % ========================================================================= % SIMULATION % ========================================================================= % create the computational grid, with the PML outside the grid PML_size = 20; % [grid points] Nx = 128 - 2*PML_size; % [grid points] dx = 0.05e-3; % [m] kgrid = makeGrid(Nx, dx); % define position of source and interfaces source_pos = 15; % [grid points] sensor_pos = 30; % [grid points] interface_1_pos = 35; % [grid points] interface_2_pos = 60; % [grid points] % define material properties c1 = 1540; % [m/s] rho1 = 1000; % [kg/m^3] c2 = 2690; % [m/s] rho2 = 1190; % [kg/m^3] c3 = 340; % [m/s] rho3 = 1.2; % [kg/m^3] % define sound speed map medium.sound_speed = c1*ones(Nx, 1); medium.sound_speed(interface_1_pos:interface_2_pos - 1) = c2; medium.sound_speed(interface_2_pos:end) = c3; % define density medium.density = rho1*ones(Nx, 1); medium.density(interface_1_pos:interface_2_pos - 1) = rho2; medium.density(interface_2_pos:end) = rho3; % set the simulation time to capture the reflections t_end = Nx*dx/min(medium.sound_speed); % define the time array CFL = 0.05; [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, CFL, t_end); % define source, offset from the edge to prevent smoothing from wrapping % the source source_mag = 1e6; source.p0 = zeros(Nx,1); source.p0(source_pos) = source_mag; % create a binary sensor mask sensor_mask = zeros(Nx,1); sensor_mask(sensor_pos) = 1; sensor.mask = sensor_mask; % run the simulation sensor.record = {'p'}; sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor,... 'PMLInside', false, 'PlotLayout', true, ... 'PlotScale', [-1, 1]*source_mag, 'PlotFreq', 10); % ========================================================================= % VISUALISATION % ========================================================================= % plot the recorded time signals figure; plot(kgrid.t_array, sensor_data.p); xlabel('time step'); ylabel('pressure at the sensor position');</code></pre> cge 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-4638 Mon, 14 Jul 2014 05:41:26 +0000 cge 4638@http://www.k-wave.org/forum/ Hi Hypermoi Thank you very much for your suggestion. I will take a try Best Regards Anthony 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-4635 Fri, 11 Jul 2014 11:57:13 +0000 Anthony 4635@http://www.k-wave.org/forum/ Hi cge, I am not part of the development team but I would say that k-wave is accurate only for weakly heterogeneous media. What if you try other values for the second medium sound speed, density and attenuation? Is there a threshold on the difference of impedance under which this strange peak disappears? Intuitively, it seems natural that a spectral method is not comfortable with discontinuities. I have not investigated in details the references mentioned at the bottom of the 21st page of the user's manual but maybe you can try to play with the CFL number. You could also smooth your sound speed and density distributions (see p63 of the user's manual). Maybe you could also set a source mask at the place where the interface should be and use Dirichlet boundary conditions (see p34), in order to have a reflecting interface. Best regards, Anthony cge 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-4622 Tue, 08 Jul 2014 06:57:34 +0000 cge 4622@http://www.k-wave.org/forum/ Hi everyone Recently I am helping another guy in our group to simulate the reflection at the hetrogeneous media boundary. I started from a simple 1D simulation, the code is provided below. clear all; % ========================================================================= % SIMULATION % ========================================================================= % create the computational grid Nx = 64; % number of grid points in the x (row) direction dx = 0.05e-3; % grid point spacing in the x direction [m] kgrid = makeGrid(Nx, dx); % define the properties of the propagation medium medium.sound_speed = 1540*ones(Nx, 1); % [m/s] medium.sound_speed(11:50) = 2690; medium.sound_speed(51:64) = 340; % [m/s] medium.density = 1000*ones(Nx, 1); % [kg/m^3] medium.density(11:50) = 1190; % [kg/m^3] medium.density(51:64) = 1.204; medium.alpha_power=2; medium.alpha_coeff=0.002*ones(Nx,1); medium.alpha_coeff(11:50)=1.19; medium.alpha_coeff(51:64)=1.64; % create initial pressure distribution using a smoothly shaped sinusoid source_mask=zeros(Nx,1); source_mask(1)=1; source.p0 =10^6*source_mask; % create a Cartesian sensor mask sensor_mask=zeros(Nx,1) sensor_mask(9)=1; sensor.mask = sensor_mask; % [mm] % set the simulation time to capture the reflections t_end = dx*10/1540+dx*40/2690+dx*14/340; % define the time array [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, [], t_end); % run the simulation sensor.record={'p'}; sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor, 'PlotLayout', true); % ========================================================================= % VISUALISATION % ========================================================================= % plot the recorded time signals figure; plot(kgrid.t_array, abs(sensor_data.p)); xlabel('time step'); ylabel('pressure at the sensor position'); As for the result, a interesting happens, the excitation has a amplitude only of 10^6 Pa, however, the result I got in this end is quite huge, almost 10^60~70 pa, almost infinity and impossible. I checked the code but can not find the problem. Can anyone helo? Best Regards