k-Wave User Forum » Topic: model shear wave scattering from hole

k-Wave User Forum » Topic: model shear wave scattering from hole http://www.k-wave.org/forum/topic/model-shear-wave-scattering-from-hole Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:29:15 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php ad000000 on "model shear wave scattering from hole" http://www.k-wave.org/forum/topic/model-shear-wave-scattering-from-hole#post-6401 Mon, 09 Apr 2018 19:34:55 +0000 ad000000 6401@http://www.k-wave.org/forum/ Thanks Brad! Bradley Treeby on "model shear wave scattering from hole" http://www.k-wave.org/forum/topic/model-shear-wave-scattering-from-hole#post-6393 Sun, 08 Apr 2018 18:27:36 +0000 Bradley Treeby 6393@http://www.k-wave.org/forum/ Hi ad000000, The dispersion is due to the low-order numerical schemes typically used in FEM. Switching to COMSOL is not likely to solve this issue. For a scatterer in an otherwise homogeneous domain, a combined FEM-BEM approach could work well. There are lots of open source FEM projects (e.g., <a href="https://fenicsproject.org">FEniCS</a>), but I haven't really used them enough to make a recommendation. Brad. ad000000 on "model shear wave scattering from hole" http://www.k-wave.org/forum/topic/model-shear-wave-scattering-from-hole#post-6373 Mon, 26 Mar 2018 18:14:23 +0000 ad000000 6373@http://www.k-wave.org/forum/ Thanks Brad. I really like K-wave, for the versatility. I'll try the convergence test. However I take your point about FEM being a better choice for elastic circles, large discontinuities in material properties etc. Can you recommend a FEM program for this? I know that PZflex has a core solver that is susceptible to considerable numerical dispersion over distances on the order of hundreds of wavelengths. I don't know if Comsol has same problem, but it's a bit pricey... Bradley Treeby on "model shear wave scattering from hole" http://www.k-wave.org/forum/topic/model-shear-wave-scattering-from-hole#post-6363 Wed, 21 Mar 2018 21:43:21 +0000 Bradley Treeby 6363@http://www.k-wave.org/forum/ Hi ad000000, Glad you managed to solve most of the issues. In general, modelling a circle in an elastic medium with a very large impedance mismatch is a challenging problem for the type of numerical model used in k-Wave (FEM would handle this much better). There are two issues. The first is staircasing, which effects elastic simulations more than fluid simulations (see <a href="http://bug.medphys.ucl.ac.uk/papers/2017-Robertson-JASA.pdf">this paper</a> for some details). The second is the large discontinuity in medium properties. Both of these can be solved by increasing the sampling. If you want to gain some confidence in your numerical results, I would suggest a convergence test, i.e., reduce dx and dt, and check you get the same answer. Brad. ad000000 on "model shear wave scattering from hole" http://www.k-wave.org/forum/topic/model-shear-wave-scattering-from-hole#post-6359 Wed, 21 Mar 2018 02:18:32 +0000 ad000000 6359@http://www.k-wave.org/forum/ ...continued: 1. Assigning uy=signal, ux=0 is NOT the way to make a shear wave. Need to use sxy. 2. I was able to set the hole material to very low density and pretty low velocity. Max attenuation I could set was pretty low, at about 20 dB/cm. This still gave substantial oscillations withing the hole, but I don't *think* they altered the scattering of the main wave much. Thus I think I got a reasonable approximation for scattering from a hole. I also tried using Dirichlet BCs and setting the hole as a source with u=0, but this only worked for a longitudinal incident wave, since I could not use the "u" source to make a shear wave, and could not use Dirichlet BC for "s" sources. Setting u=0 actually models an infinitely hard "hole", rather than air. Even so, the program still crashed. Apparently using the Dirichlet approach to make hard boundaries requires much care. ad000000 on "model shear wave scattering from hole" http://www.k-wave.org/forum/topic/model-shear-wave-scattering-from-hole#post-6358 Wed, 21 Mar 2018 02:04:15 +0000 ad000000 6358@http://www.k-wave.org/forum/ Follow-up: I solved these issues, mostly. ad000000 on "model shear wave scattering from hole" http://www.k-wave.org/forum/topic/model-shear-wave-scattering-from-hole#post-6354 Tue, 20 Mar 2018 05:23:45 +0000 ad000000 6354@http://www.k-wave.org/forum/ Hello, I'm trying to model scattering from a plane shear wave (signal=wavelet) from a hole. My material is steel, and I set the hole material velocity and density to be less than steel. I get several problems: 1. I can make a compression plane wave emanating from the left edge of the grid, by setting source.ux=0 and source.uy=signal. But if I reverse them, to make a shear wave, it gives nonsense results. 2. To simulate a hole (air inside), I set the density very low. However the "hole" seems to generate its own energy. I tried to kill that by setting attenuation high, but the program ground to a halt. Thanks for any thoughts! Code: clear source sensor medium kgrid % define variables freq=5e6; epw=12; cp=5850; cs=3230; rho=8030; SDHdia=1.5e-3; %==define simulation area (remember, horizontal is Y direction)================== xlen=20e-3; ylen=16e-3; % ==================== setup grid ========================================= lambda=cp/freq; dx=lambda/epw; dy=lambda/epw; Nx= ceil(xlen/dx); Ny= ceil(ylen/dy); kgrid = makeGrid(Nx, dx, Ny, dy); %====set up material regions========================================= medium.sound_speed_compression =ones(kgrid.Nx, kgrid.Ny)*cp; medium.sound_speed_shear =ones(kgrid.Nx, kgrid.Ny)*cs; medium.alpha_coeff_compression = ones(kgrid.Nx, kgrid.Ny)*.03/(freq/1e6); %to give dB/(MHz^2 cm) medium.alpha_coeff_shear = ones(kgrid.Nx, kgrid.Ny)*.03/(freq/1e6)*4; medium.density=ones(kgrid.Nx, kgrid.Ny)*rho; SDHloc(1)=0; SDHloc(2)=2e-3; SDHdia=SDHdia/dx; [SDHcenterGridx,SDHcenterGridy]=cartPt2Grid(kgrid,[SDHloc(1) SDHloc(2)]); disc = makeDisc(Nx, Ny, SDHcenterGridx, SDHcenterGridy, SDHdia/2); medium.sound_speed_compression(disc == 1) = cp*.9; medium.sound_speed_shear(disc == 1) = cs*.9; medium.density(disc==1)=60; medium.alpha_coeff_compression(disc==1) = 1/(freq/1e6); %to give dB/(MHz^2 cm) medium.alpha_coeff_shear(disc==1) = 4/(freq/1e6)*4; t_end=.3*2*ylen/cs; %if ReverseFire==1; t_end=t_end/2;end dt=.3*dx/(cp); %set time steps (C should be Cmax) kgrid.t_array = 0:dt:t_end; Nt=length(kgrid.t_array); PML=20*dx; %======================= setup source =================================================== signal = wavelet(freq,4,0,2/freq,kgrid.t_array); source.u_mask = zeros(kgrid.Nx,kgrid.Ny); % or source.p_mask source.u_mask(:,1)=1; source.ux=signal/(rho*cp); source.uy=0; %====setup sensor ===================================================== sensor.mask=zeros(kgrid.Nx, kgrid.Ny); %sensor.mask(1:5:kgrid.Nx,1:5:kgrid.Ny)=1; sensor.mask(1:2:kgrid.Nx,1:2:kgrid.Ny)=1; sensor.record = {'u_max_all','u'}; %{'p','p_max'}; %===== display_mask ========================================= display_mask=disc; %==================== RUN ================================================= figpos=get(0,'DefaultFigurePosition'); set(0,'DefaultFigurePosition',[340,275,980,735]); input_args = {'DisplayMask', display_mask, 'PMLInside', true, 'PlotPML', false,'PlotScale','auto','DataCast', 'single'};%,'RecordMovie',true}; sensor_data = pstdElastic2D(kgrid, medium, source, sensor, input_args{:}); %sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); set(0,'DefaultFigurePosition',figpos); kmodel=sensor_data; clear sensor_data; kmodel.freq=freq; kmodel.Nx=Nx;kmodel.Ny=Ny; kmodel.x_vec=kgrid.x_vec;kmodel.y_vec=kgrid.y_vec; kmodel.Nt=length(kgrid.t_array); kmodel.epw=epw; save kmodel.mat kmodel -v7.3 uall=sqrt(kmodel.ux_max_all.^2+kmodel.uy_max_all.^2); u=sign(kmodel.ux.*kmodel.uy).*sign(kmodel.ux).*sqrt(kmodel.ux.^2+kmodel.uy.^2); %uh=abs(hilbert(u'));uh=uh'; %nox=ceil(kmodel.Nx/5);noy=ceil(kmodel.Ny/5); nox=ceil(kmodel.Nx/2);noy=ceil(kmodel.Ny/2); dat=zeros(nox,noy);%may need to reverse this is Nx, Ny are in wrong order for t=1:size(u,2); for j=1:nox for k=1:noy; %dat(j,k,t)=uh((k-1)*nox+j,t); dat(j,k,t)=u((k-1)*nox+j,t); %daty(j,k,t)=kmodel.uy((k-1)*nox+j,t); end end end % Record the movie %maxdat=max(max(max(dat))); %ax=[0 maxdat/500]; nframes=50; step=floor((kmodel.Nt-1)/nframes); figure;imagesc(kmodel.y_vec,kmodel.x_vec,(dat(:,:,200)));axis equal;%caxis(ax) set(gca,'nextplot','replacechildren'); ax=caxis; for j= 1:nframes; t=1+(j-1)*step; imagesc(kmodel.y_vec,kmodel.x_vec,(dat(:,:,t)));caxis(ax) F(j) = getframe; end % Play the movie movie(F,1,1)