k-Wave User Forum » Topic: Simulation of a flat reflector

k-Wave User Forum » Topic: Simulation of a flat reflector http://www.k-wave.org/forum/topic/simulation-of-a-flat-reflector Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 00:42:23 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php rainbowwhale on "Simulation of a flat reflector" http://www.k-wave.org/forum/topic/simulation-of-a-flat-reflector#post-7374 Tue, 07 Apr 2020 03:24:36 +0000 rainbowwhale 7374@http://www.k-wave.org/forum/ CFL number should be calculated using highest speed of sound. By using highest speed of sound, current CFL number is about 4.13. I tried using dt as 0.05/(PPP*F0) which has CFL number of about 0.2 and it worked. Ryuji on "Simulation of a flat reflector" http://www.k-wave.org/forum/topic/simulation-of-a-flat-reflector#post-7363 Sat, 04 Apr 2020 02:07:08 +0000 Ryuji 7363@http://www.k-wave.org/forum/ Hi, thanks for providing the great library. I would like to simulate a small 40kHz-ultrasound transducer with a flat reflector in 3D. When I set the medium heterogeneous (air and acrylic), I cannot see any propagation even in the air. Is this because of the large change in impedance? I would appreciate it if you could tell me how to make my code work. Thanks advance! Best regards, Ryuji % ========================================================================= % sorce code % ========================================================================= clear all; % ========================================================================= % DEFINE LITERALS % ========================================================================= % select which code to run % 1: MATLAB CPU code % 2: MATLAB GPU code % 3: C++ code MODEL = 2; % scale parameter (changes the resolution of the simulation) SC = 1; % grid parameters PML_SIZE = 10; % size of the perfectly matched layer [grid points] Nx = 128*SC - 2*PML_SIZE; % number of grid points in the x direction Ny = 64*SC - 2*PML_SIZE; % number of grid points in the y direction Nz = 64*SC - 2*PML_SIZE; % number of grid points in the z direction % sampling parameters PPW = 10*SC; % points per wavelength (this defines the grid size) PPP = 10*SC; % points per period (this defines the CFL) T_END = 1e-3; % total compute time [s] (this must be long enough to reach steady state) RECORD_PERIODS = 2; % number of periods to record % medium parameters C0 = 346; % sound speed [m/s] RHO0 = 1.2; % density [kg/m^3] CR = 1.43e3; % sound speed in the reflector [m/s] RHOR = 1.18e3; % density of the reflector [kg/m^3] % source parameters F0 = 40e3; % source frequency [Hz] SOURCE_RADIUS = 5e-3; % piston radius [m] SOURCE_MAG = 0.5e6; % source pressure [Pa] % ========================================================================= % CREATE GRID % ========================================================================= % calculate the grid spacing based on the PPW and F0 dx = C0 / (PPW * F0); % [m] dt = 1 / (PPP*F0); % [s] SOURCE_RADIUS = SOURCE_RADIUS / dx; % calculate the CFL disp(['CFL = ' num2str(C0*dt/dx)]); % create the computational grid kgrid = makeGrid(Nx, dx, Ny, dx, Nz, dx); % create the time array kgrid.t_array = 0:dt:T_END; % assign medium properties % medium.sound_speed = C0; % medium.density = RHO0; medium.sound_speed = C0 * ones(Nx, Ny, Nz); medium.density = RHO0 * ones(Nx, Ny, Nz); % assign reflector properties medium.sound_speed(Nx, :, :) = CR; medium.density(Nx, :, :) = RHOR; % assign the reference sound speed to the background medium % medium.sound_speed_ref = C0; % ========================================================================= % CREATE SOURCE % ========================================================================= % create piston source mask piston = makeDisc(Ny, Nz, Ny/2, Nz/2, SOURCE_RADIUS); source.p_mask = zeros(Nx, Ny, Nz); source.p_mask(1, :, :) = piston; % create time varying source source.p = SOURCE_MAG * sin(2*pi*F0*kgrid.t_array); % filter source with up-ramp ramp_length = round((2*pi/F0)/dt); ramp = (-cos( (0:(ramp_length-1))*pi/ramp_length ) + 1)/2; source.p(1:ramp_length) = ramp .* source.p(1:ramp_length); % ========================================================================= % CREATE SENSOR MASK % ========================================================================= % set sensor mask to record central axis, not including the source point % sensor.mask = zeros(Nx, Ny, Nz); % sensor.mask(2:end, Ny/2, Nz/2) = 1; sensor.mask = [2, 1, Nz/2, Nx, Ny, Nz/2].'; % record the maximum pressure sensor.record = {'p', 'p_rms'}; % average only the final few periods when the field is in steady state sensor.record_start_index = kgrid.Nt - RECORD_PERIODS*PPP + 1; % ========================================================================= % RUN k-WAVE SIMULATION % ========================================================================= % run code switch MODEL case 1 % MATLAB CPU sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ... 'DataCast', 'single', 'PMLInside', false, ... 'DisplayMask', source.p_mask, 'PlotScale', [-1, 1]*SOURCE_MAG); case 2 % MATLAB GPU sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ... 'DataCast', 'gpuArray-single', 'PMLInside', false, ... 'DisplayMask', source.p_mask, 'PlotScale', [-1, 1]*SOURCE_MAG); case 3 % C++ sensor_data = kspaceFirstOrder3DC(kgrid, medium, source, sensor, 'PMLInside', false); otherwise end % ========================================================================= % ANALYTICAL SOLUTION % ========================================================================= % calculate the wavenumber k = 2*pi*F0./C0; % define radius and axis a = SOURCE_RADIUS*dx; x = (kgrid.x_vec(2:end, :) - min(kgrid.x_vec(:))); % calculate the analytical solution for a piston in an infinite baffle % for comparison (Eq 5-7.3 in Pierce) r_a = sqrt(x.^2 + a^2); p_ref = SOURCE_MAG*abs(-2i*exp(1i*k*(x+r_a)/2).*sin((k*r_a - k*x)/2)); % ========================================================================= % VISUALISATION % ========================================================================= % plot the pressure along the focal axis of the piston figure; plot(sensor_data.p(end,:), 'bx'); img = gather(sensor_data.p_rms); image(img, 'CDataMapping', 'scaled') daspect([1 1 1]) colormap default colorbar