k-Wave User Forum » Topic: Failure using dirichlet source

k-Wave User Forum » Topic: Failure using dirichlet source http://www.k-wave.org/forum/topic/failure-using-dirichlet-source Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 00:59:55 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Failure using dirichlet source" http://www.k-wave.org/forum/topic/failure-using-dirichlet-source#post-5999 Mon, 19 Jun 2017 21:02:01 +0000 Bradley Treeby 5999@http://www.k-wave.org/forum/ Hi Bill, Full-wave models like k-Wave that implement heterogeneities directly in the material properties (rather than as explicit boundary conditions) do not deal with with such strong density contrasts. If you're trying to model the response of a PZT transducer in air, k-Wave or other similar models are probably not the best tool. You'd be better served using a FEM model. Brad. bwalker000 on "Failure using dirichlet source" http://www.k-wave.org/forum/topic/failure-using-dirichlet-source#post-5984 Wed, 14 Jun 2017 06:17:35 +0000 bwalker000 5984@http://www.k-wave.org/forum/ The following code is numerically unstable and blows up. I'm not sure why. The same code 'works' if I change to an additive source. That said, the pressures inside the additive source region are very inhomogeneous, which I do not understand. This instability doesn't make much sense to me and I'm having a hard time developing confidence in the model. Any feedback would be appreciated. Best, Bill Walker clear cfl = 0.1; % Courant-Friedrichs-Lewy number %% Set geometry N = 64; Nx = N; Ny = Nx; Nz = N; dx = 0.2e-3; dy = dx; dz = dx; % number of samples for the lossy boundary layers n_bnd = 6; %% x = ((-(Nx-1)/2):((Nx-1)/2))*dx; y = x; z = ((-(Nz-1)/2):((Nz-1)/2))*dx; z = z; [X,Y,Z] = meshgrid(x,y,z); R = ((X).^2+(Y).^2).^0.5; kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz); %% Set pulse characteristics and timing f0 = 0.5e6; N_cycles = 8; t_max = 20e-6; %% Define the object regions % PZT 5A Source r0(1) = 3.29e-3 / 2; r1(1) = 9e-3 / 2; z0(1) = -1e-3; z1(1) = 1e-3 + eps; material(1) = 1; %% Fill out the geometry matrix mat = zeros(size(X)); for j = 1:length(material), ind = find((r0(j) < R) & (r1(j) >= R) & (z0(j) < Z) & (z1(j) >= Z)); mat(ind) = material(j); end %% Set the material properties % define the properties of the surrounding medium medium.sound_speed_compression = 2*2^2*330*ones(size(mat)); % [m/s] ----- Made up medium.sound_speed_shear = 0*ones(size(mat)); % [m/s] ----- Made up medium.density = 1.225*ones(size(mat)); % [kg/m^3] ----- Made up % define the properties of the PZT 5A ind = find(mat == 1); medium.sound_speed_compression(ind) = 4350; % [m/s] Szabo book % poisson's ratio = 0.35 medium.sound_speed_shear(ind) = sqrt(1/(2*(1+0.35)))*4350; % [m/s] medium.density(ind) = 7750; % [kg/m^3] Szabo book figure(1) imagesc(squeeze(medium.density(:,Nx/2,:))); axis equal axis off %% Define the source % create the time array [kgrid.t_array, dt] = makeTime(kgrid, max(medium.sound_speed_compression(:)), cfl, t_max); clear kgrid.t_array kgrid.t_array = 0:dt:t_max; % define a ring source element source.u_mask = zeros(size(mat)); ind = find(mat == 1); source.u_mask(ind) = 1; tone_burst_freq = f0; tone_burst_cycles = N_cycles; source_signal = sin(2*pi*tone_burst_freq*(0:dt:(tone_burst_cycles/tone_burst_freq))); source_signal = source_signal.*hann(length(source_signal))'; fir_coeff = fir1(128,1.5e6*dt*2,'low'); source_signal = conv(fir_coeff,source_signal); ind = find(mat == 1); source.uz = zeros(length(ind),length(source_signal)); for j = 1:length(ind), source.uz(j,:) = source_signal; end source.u_mode = 'dirichlet'; %% Define the sensor points % define a series of Cartesian points to collect the data ind = find(abs(X - x(Nx/2)) < dx/2); sensor.mask = [kgrid.x(ind)'; kgrid.y(ind)'; kgrid.z(ind)']; %% Clean up some more memory clear ind mat R %% compute the pressure field % define the field parameters to record sensor.record = {'p','u'}; % input arguments input_args = {'DisplayMask', source.u_mask, 'DataCast', 'single', 'CartInterp', 'nearest'}; % run the simulation sensor_data = pstdElastic3D(kgrid, medium, source, sensor, input_args{:}); save(['simple_test']) close all