k-Wave

1. unifi
   

   Member
   Joined: Nov '13
   Posts: 1

   Hi, I'm trying to model an interface between two mediums (air and some solid) and construct a matched impedance surface.
   If the characteristic impedance is Z=rho*c, I try to increase the density of the second medium while reducing by the same factor the sound speed, so that the values of Z in both mediums are identical.
   I would expect that if the two impedances match, there would be no reflection. Nevertheless I see a lot of reflection. If the rescaling factor R is big enough (c/R, rho*R) almost all the wave is reflected.

   I'm testing it with a sine signal

   I guess the answer is really simple but still I can't figure it.
   Can you please explain me why?

   Thanks

   Posted 10 years ago #

2. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi unifi,

   Using the simple 1D example below, I can see the same behaviour that you describe when the rescaling factor R is very large. In this case, the spurious reflection is due to frequency content of the incident wave not being supported on the grid with the much lower sound speed. You can calculate the maximum supported frequency in each part of the domain using f_max = c_0/(2*dx). What happens if you change your source frequency (or grid spacing) so that it is supported in both the air and solid media?

   Hope that helps,

   Brad.

   --

   % create the computational grid
   Nx = 512;       % 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
   R    = 10;
   c1   = 1500;
   rho1 = 1000;
   c2   = c1/R;
   rho2 = rho1*R;
   
   % create medium
   medium.sound_speed         = c1*ones(Nx, 1);
   medium.sound_speed(1:Nx/2) = c2;
   medium.density             = rho1*ones(Nx, 1);
   medium.density(1:Nx/2)     = rho2;
   
   % create initial pressure distribution using a smoothly shaped sinusoid
   x_pos = 320;    % [grid points]
   width = 50;     % [grid points]
   height = 1;     % [pa]
   in = (0:pi/(width/2):2*pi).';
   source.p0 = [zeros(x_pos, 1); ((height/2)*sin(in-pi/2)+(height/2)); zeros(Nx - x_pos  - width - 1, 1)];
   
   % run the simulation
   sensor_data = kspaceFirstOrder1D(kgrid, medium, source, []);

   Posted 10 years ago #

3. camille
   

   Member
   Joined: Apr '15
   Posts: 1

   Hi,
   I’m trying to model a transducer which is composed of a ceramic layer and a matching layer and two others mediums (tissues +bone). I would like to change the matching layer impedance to see if there are less reflections if we adapt the matching layer for bone ( Z matching layer = sqrt(Z ceramic*Z bone) ).
   I also change the characteristics of the matching layer impedance (Z = rho*c). Indeed, I increase the sound speed and I decrease the density to keep the same impedance, so normally my signal shouldn’t change because reflections depend of the impedance but it changes…
   The signal from the transducer is 2 MHz. A sensor is located in the bone medium to capture the amplitude of the signal.

   The source code of the program is given below.

   clear all;

   % simulation settings
   DATA_CAST = 'single';

   % =========================================================================
   % DEFINE THE K-WAVE GRID
   % =========================================================================

   % set the size of the perfectly matched layer (PML)
   PML_X_SIZE = 20; % [grid points]
   PML_Y_SIZE = 10; % [grid points]
   PML_Z_SIZE = 10; % [grid points]

   % set total number of grid points not including the PML
   Nx = 128 - 2*PML_X_SIZE; % [grid points]
   Ny = 128 - 2*PML_Y_SIZE; % [grid points]
   Nz = 64 - 2*PML_Z_SIZE; % [grid points]

   % set desired grid size in the x-direction not including the PML
   x = 20e-3; % [m]

   % calculate the spacing between the grid points
   dx = x/Nx; % [m]
   dy = dx; % [m]
   dz = dx; % [m]

   % create the k-space grid
   kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);

   % =========================================================================
   % DEFINE THE MEDIUM PARAMETERS
   % =========================================================================

   % define the properties of the propagation medium
   v_ceramique = 3947; %% speed of sound in the ceramic [m/s]
   v_water = 1500; % speed of sound in water [m/s]
   v_bone = 2850; % speed of sound in bone [m/s]
   v_ML = 2600; % speed of sound in the matching layer [m/s]
   v_tissue = 1540; % speed of sound in tissue [m/s]
   den_ceramique = 7600; % density of ceramic [kg/m^3]
   den_water = 1000; % density of water [kg/m^3]
   den_bone = 2596; % density of steel [kg/m^3]
   den_ML = 5730; % density of matching layer [kg/m^3]
   den_tissue = 1052; % density of tissue [kg/m^3]
   medium.sound_speed = v_ceramique*ones(Nx,Ny,Nz); % [m/s]
   medium.sound_speed(2:3,:,:) = v_ML;
   medium.sound_speed(3:12,:,:) = v_tissue;
   medium.sound_speed(12:40,:,:) = v_bone; %[m/s]
   medium.sound_speed(40:60,:,:) = v_bone; %[m/s]
   medium.sound_speed(60:end,:,:) = v_water;
   medium.density = den_ceramique*ones(Nx,Ny,Nz); % [m/s]
   medium.density(2:3,:,:) = den_ML;
   medium.density(2:12,:,:) = den_tissue;
   medium.density(12:40,:,:) = den_bone; %[m/s]
   medium.density(40:60,:,:) = den_bone; %[m/s]
   medium.density(60:end,:,:) = den_water;
   % create the time array
   t_end = 40e-6; % [s]
   kgrid.t_array = makeTime(kgrid, 3947, [], t_end);

   % =========================================================================
   % DEFINE THE INPUT SIGNAL
   % =========================================================================

   % define properties of the input signal
   source_strength = 1e6; % [Pa]
   tone_burst_freq = 2e6; % [Hz]
   tone_burst_cycles = 4.5;

   % create the input signal using toneBurst
   input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);

   % scale the source magnitude by the source_strength divided by the
   % impedance (the source is assigned to the particle velocity)
   input_signal = (source_strength./(1500*1000)).*input_signal;

   % =========================================================================
   % DEFINE THE ULTRASOUND TRANSDUCER
   % =========================================================================

   % physical properties of the transducer
   transducer.number_elements = 64; % total number of transducer elements
   transducer.element_width =1.155; % width of each element [grid points]
   transducer.element_length = 20; % length of each element [grid points]
   transducer.element_spacing =0; % spacing (kerf width) between the elements [grid points]
   transducer.radius = inf; % radius of curvature of the transducer [m]

   % calculate the width of the transducer in grid points
   transducer_width = transducer.number_elements*transducer.element_width ...
   + (transducer.number_elements - 1)*transducer.element_spacing;

   % use this to position the transducer in the middle of the computational grid
   transducer.position = round([1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]);

   % properties used to derive the beamforming delays
   transducer.sound_speed = 1540; % sound speed [m/s]

   % apodization
   transducer.transmit_apodization = 'Rectangular';
   transducer.receive_apodization = 'Rectangular';

   % define the transducer elements that are currently active
   transducer.active_elements = zeros(transducer.number_elements, 1);
   transducer.active_elements(30:30) = 1;

   % append input signal used to drive the transducer
   transducer.input_signal = input_signal;

   % create the transducer using the defined settings
   transducer = makeTransducer(kgrid, transducer);

   % print out transducer properties
   transducer.properties;

   % =========================================================================
   % DEFINE SENSOR MASK
   % =========================================================================

   % create a Sensor

   sensor.mask = zeros(Nx, Ny, Nz);
   sensor.mask(40, 40, 20) = 1;
   % =========================================================================
   % RUN THE SIMULATION
   % =========================================================================

   % set the input settings
   input_args = {'PlotLayout', true,'DisplayMask', sensor.mask, ...
   'PMLInside', true, 'PlotPML', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
   'PMLInside', false, 'DataCast', DATA_CAST, 'PlotScale', [-source_strength/2, source_strength/2]};

   % run the simulation
   sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, sensor, input_args{:});

   % extract a single scan line from the sensor data using the current
   % beamforming settings
   scan_line = transducer.scan_line(sensor_data);

   %p_xyz = kspacePlaneRecon(sensor_data_rs, kgrid.dy, kgrid.dz, dt, medium.sound_speed,...
   % 'DataOrder', 'yzt', 'PosCond', true, 'Plot', true);

   % =========================================================================
   % VISUALISATION
   % =========================================================================

   figure;
   plot(kgrid.t_array*1e6, sensor_data, 'k-');
   xlabel('Time [\mus]');
   ylabel('Pressure [au]');
   title('Scan Line After Beamforming');

   Posted 9 years ago #

4. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Camille,

   In what way does your signal change that you don't expect?

   If the impedance remains identical, the plane wave reflection coefficient won't change. However, for non-normal incidence, the direction of the wave will change due to refraction (which depends specifically on the sound speed, not the impedance). If you have a non-plane wave, the shape of the field will thus change when you change the properties of the layer.

   Hope that helps,

   Brad.

   Posted 9 years ago #

RSS feed for this topic

Reply

You must log in to post.