k-Wave User Forum » Topic: Simulation of Ultrasonic Testing (Borehole in Solid)

k-Wave User Forum » Topic: Simulation of Ultrasonic Testing (Borehole in Solid) http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 22:33:30 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6735 Tue, 29 Jan 2019 22:36:57 +0000 Bradley Treeby 6735@http://www.k-wave.org/forum/ Hi Rola, No, keep the velocity the same. The sound speeds are much closer in magnitude than the density (factor of 5 vs 1000). Brad. Rola on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6709 Tue, 15 Jan 2019 19:18:22 +0000 Rola 6709@http://www.k-wave.org/forum/ Hi Brad and Ben, If I inflate the air density artificially by a factor of 100, do I need to inflate the inflate the air velocity too? Best bencox on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6687 Sun, 09 Dec 2018 23:28:27 +0000 bencox 6687@http://www.k-wave.org/forum/ HI fabwo, You're putting in a toneburst at 40 MHz, so it's not surprising it's changing when you filter at 100 kHz or so. Perhaps use the <code>...'PlotSpectrums', true)</code> option in <code>filterTimeSeries</code> to see the spectrum of your signal before and after filtering. If you are only interested in propagation in one direction, then could you use the 1D code? It would be much more computationally efficient. Best wishes, Ben fabwo on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6574 Fri, 31 Aug 2018 18:01:45 +0000 fabwo 6574@http://www.k-wave.org/forum/ Hi Ben, I adjusted my grid dimensions in x-direction (wave propagation is only in x-direction), shortened the time steps and tried different tone burst frequencies. Unfortunately I didn't get better results. I also have some questions about the k-wave filter function. My grid dimension without the PML are Nx=1004, Ny=16, Nz=16 with dx=2.5e-5, dy=1e-3, dz=1e-3 My time steps are created using the function kgrid.makeTime with cfl=0.3, t_end=2e-5. This leads to 8002 time steps in the kgrid.t_array with 2.5e-9 seconds apart, but I have also tried shorter ones (e.g. dt=1e-9). I modelled time varying stress at each of the source positions only in xx using the toneBurst function. toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst-cycles) with tone_burst_freq=4e7 (4Mhz), tone_burst-cycles=3 The filterTimeSeries function uses by default 3 PPW of the largest grid point spacing to compute the cutoff frequency. Is there a way to consider only the grid spacing in x-direction (because it is the only direction the waves are propagating)? At the moment my filter cutoff frequency is 114.3333 kHz (3PPW). But somehow this filter changes even a tone burst with a frequency of 100 kHz. The amplitude of the tone burst is reduces by a factor of 100. How is this possible? When I simulated with a tone burst of 100 kHz the pressure in sensor_data still blows up. Do you know where my mistakes are? <a href="https://ibb.co/iQfJee" rel="nofollow">https://ibb.co/iQfJee</a> <a href="https://ibb.co/cOZdee" rel="nofollow">https://ibb.co/cOZdee</a> <a href="https://ibb.co/cRx9kK" rel="nofollow">https://ibb.co/cRx9kK</a> <a href="https://ibb.co/nD9dee" rel="nofollow">https://ibb.co/nD9dee</a> Best wishes, Fabian % ========================================================================= % 3D ULTRASONIC SIMULATION - BOREHOLE IN A FRP COMPONENT % CIRCULAR TRANSDUCER % ========================================================================= % ========================================================================= % DEFINITION OF THE K-WAVE GRID % ========================================================================= % PML size (perfectly matched layer) on the outside of the k-wave grid % Grid points == GP PML_X_SIZE = 10; % [GP] PML_Y_SIZE = 8; % [GP] PML_Z_SIZE = 8; % [GP] PML_SIZE = [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE]; % Computational grid: number of GP without PML % The number of GP (computational grid + PML) in each dimension should be % given by a power of 2 or have a small primefactor because of % computational advantages for later FFT use. Nx = 1024 - 2*PML_X_SIZE; % [GP] Ny = 32 - 2*PML_Y_SIZE; % [GP] Nz = 32 - 2*PML_Z_SIZE; % [GP] % Distance between GP dx = 2.5e-5; % [m] dy = 1e-3; % [m] dz = 1e-3; % [m] % Creation of the k-space grid kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); %% ======================================================================== % TEST SET-UP & MEDIUM PARAMETER & TIME ARRAY % ========================================================================= % Test set-up frp_thickness = 5.5e-3; % [m] FRP plate thickness hole_diameter = 4e-3; % [m] borehole diameter hole_depth = 2e-3; % [m] borehole depth wedge_diameter = 6e-3; % [m] PMMA wedge diameter wedge_length = 11e-3; % [m] PMMA wedge length % Transformation [m] into [GP] frp_thickness_GP = single(frp_thickness / dx); % [GP] plate thickness hole_radius_GP = single((hole_diameter/2) / dy); % [GP] borehole radius hole_depth_GP = single(hole_depth / dx); % [GP] borehole depth wedge_radius_GP = single((wedge_diameter/2) / dy); % [GP] wedge radius wedge_length_GP = single(wedge_length / dx); % [GP] wedge lenth %% Medium properties c_gel = 1500; % [m/s] speed of sound gel rho_gel = 1002; % [kg/m^3] density gel c_air = 343; % [m/s] speed of sound air % The simulation is unstable for the correct density of air (1.2 [kg/m^3]) % because of the very large impedance contrast. To fix this problem the air % density is artificially inflated by factor 100. The reflection % coefficient is still close to 1. rho_air = 120; % [kg/m^3] density air cp_FRP = 3000; % [m/s] speed of sound FRP Compression (longitudinal) cs_FRP = 3000; % [m/s] speed of sound FRP Shear (transversal) rho_FRP = 1500; % [kg/m^3] density FRP cp_PMMA = 2700; % [m/s] speed of sound PMMA Compression (longitudinal) cs_PMMA = 1300; % [m/s] speed of sound PMMA Shear (transversal) rho_PMMA = 1180; % [kg/m^3] density PMMA %% Medium model: % Creating matrices with the dimension of the k-space grid for the density % and the sound speed of longitudinal and transversal waves. % Initializing medium model matrices with air medium.sound_speed_compression = c_air * ones(Nx, Ny, Nz); % [m/s] medium.sound_speed_shear = c_air * ones(Nx, Ny, Nz); % [m/s] medium.density = rho_air * ones(Nx, Ny, Nz); % [kg/m^3] % FRP plate for m = (Nx - frp_thickness_GP) : Nx medium.sound_speed_compression(m,:,:) = cp_FRP; % [m/s] medium.sound_speed_shear(m,:,:) = cs_FRP; % [m/s] medium.density(m,:,:) = rho_FRP; % [kg/m^3] end % Defects in the FRP plate (1 borehole) disc2 = makeDisc(Ny, Nz, Ny/2, Nz/2, hole_radius_GP); % position, radius flaw1 = repmat(disc2, [1, 1, Nx]); % 3D disc extrusion flaw1 = permute (flaw1, [3,1,2]); % rearranging dimensions for n = 1 : (Nx - hole_depth_GP - 1) % borehole model using 0 and 1 flaw1(n,:,:) = 0; end % Filling the flaw model with air properties in the medium matrices medium.sound_speed_compression(flaw1==1) = c_air; % [m/s] medium.sound_speed_shear(flaw1==1) = c_air; % [m/s] medium.density(flaw1==1) = rho_air; % [kg/m^3] % Gel layer between FRP plate and PMMA wedge gel_layer1 = (Nx - frp_thickness_GP - 1); % [GP] gel layer position medium.sound_speed_compression(gel_layer1,:,:) = c_gel; % [m/s] medium.sound_speed_shear(gel_layer1,:,:) = c_air; % [m/s] medium.density(gel_layer1,:,:) = rho_gel; % [kg/m^3] % PMMA wedge disc1 = makeDisc(Ny, Nz, Ny/2, Nz/2, wedge_radius_GP); % position, radius wedge = repmat(disc1, [1, 1, Nx]); % 3D disc extrusion wedge = permute (wedge, [3,1,2]); % rearranging dimensions for h = gel_layer1 : Nx % PMMA model using 0 and 1 wedge(h,:,:) = 0; end for k = 1 : (gel_layer1 - wedge_length_GP) wedge(k,:,:) = 0; end % Filling the wedge model with PMMA properties in the medium matrices medium.sound_speed_compression(wedge==1) = cp_PMMA; % [m/s] medium.sound_speed_shear(wedge==1) = cs_PMMA; % [m/s] medium.density(wedge==1) = rho_PMMA; % [kg/m^3] % Gel layer above PMMA wedge (transducer position) gel_layer2 = ((gel_layer1 - wedge_length_GP) - 1); % [GP] gel layer position medium.sound_speed_compression(gel_layer2,:,:) = c_gel; % [m/s] medium.sound_speed_shear(gel_layer2,:,:) = c_air; % [m/s] medium.density(gel_layer2,:,:) = rho_gel; % [kg/m^3] %% Time array cfl = 0.3; % default = 0.3 t_end = 2e-5; % [s] kgrid.makeTime(max(medium.sound_speed_compression,medium.sound_speed_shear), cfl, t_end); % ========================================================================= % INPUT SIGNAL % ========================================================================= % Tone burst tone_burst_freq = 4e7; % [Hz] tone_burst_cycles = 3; % Time varying stress at each of the source positions source_prefilter = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles); source.syy = 0; source.szz = 0; source.sxy = 0; source.sxz = 0; source.syz = 0; %% Filter the source to remove high frequencies not supported by the grid source.sxx = filterTimeSeries(kgrid, medium, source_prefilter); source.syy = filterTimeSeries(kgrid, medium, source.syy); source.szz = filterTimeSeries(kgrid, medium, source.szz); source.sxy = filterTimeSeries(kgrid, medium, source.sxy); source.sxz = filterTimeSeries(kgrid, medium, source.sxz); source.syz = filterTimeSeries(kgrid, medium, source.syz); % ========================================================================= % TRANSDUCER DEFINITION [SOURCE AND SENSOR] % ========================================================================= % Transducer model using 0 and 1 % Source source.s_mask = zeros(Nx, Ny, Nz); source.s_mask(1,:,:) = disc1; % same dimensions as the wedge % Sensor sensor.mask = source.s_mask; %% ======================================================================== % RUN THE SIMULATION % ========================================================================= % Assign additional input options % place PML outside of the computational grid input_args = { ... 'PMLInside', false, 'DisplayMask', sensor.mask, 'PlotPML', false, ... 'PMLSize', PML_SIZE}; % Run the simulation sensor_data = pstdElastic3D(kgrid, medium, source, sensor, input_args{:}); save('Workspace_SmallGrid.mat'); bencox on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6556 Thu, 16 Aug 2018 20:45:50 +0000 bencox 6556@http://www.k-wave.org/forum/ Hi Fabian, It's not an instability problem, in that using too high a source frequency won't lead to the simulation blowing up. However, if you use a source with a frequency higher than the grid can support, then most of the energy will just not propagate anywhere. It's lost. Look at the amplitude of your source before and after you filter it to remove the frequencies not supported on the grid. It is almost set to zero. So I suspect you're just seeing an artifact rather than a real signal. Perhaps try a lower frequency source and see what happens. Best wishes, Ben fabwo on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6552 Thu, 16 Aug 2018 18:05:25 +0000 fabwo 6552@http://www.k-wave.org/forum/ Hi Ben, thanks for your response. I will adjust my grid to prevent instability after a certain time. However I still do not understand why the peak at 0.3 microseconds in the source signal occours at 0.4 microsenconds in the sensor signal. In my oppinion the signal returns way too fast (see calculation in my second post). Will this problem dissapear with the grid changes as well? Best wishes, Fabian bencox on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6549 Thu, 16 Aug 2018 14:14:56 +0000 bencox 6549@http://www.k-wave.org/forum/ Hi Fabwo, Sorry for the delay in getting back to you. If you look at the toneburst signal you use as a source signal, 0.3 microseconds is where the peak occurs. However, the more serious problem is that the frequency you use is too high to be supported by a grid of that spacing. Best wishes, Ben fabwo on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6524 Sat, 14 Jul 2018 12:52:06 +0000 fabwo 6524@http://www.k-wave.org/forum/ Thanks Brad, inflating the density of air worked perfectly. Unfortunately I encountered a different problem which I can not solve by myself. I modelled a PMMA wedge of a simple circular transducer on top of a solid block with a borehole drilled from underneath. The source/sensor (same grid points) is placed on top of the wedge in an additional gel layer, another gel layer exists between the wedge and the solid block. The source signal is a tone burst with a frequency of 5 Mhz and 3 cycles. When I run my pstdElastic3D simulation the sensor detects an echo-signal which has a plausible shape but is reflected way too fast. I expected an echo-signal after ~14 µs (calculation below) , not after 0.3µs. Furthermore I experience very high negative pressure values after a certain simulation time. Therefore I have two questions: 1) Why does the echo-signal return so quickly to the transducer? 2) Why do I get very high negative pressure values after ~ simulation time? Graph of source and sensor signal: <a href="https://imgur.com/a/0EdSQU4" rel="nofollow">https://imgur.com/a/0EdSQU4</a> Fabian Calculation of the echo-signal time difference: t = 2 * (length_wedge / soundspeed_wedge + length_solid / soundspeed_solid) The length of the different materials is given in gridpoints... t = 2 * (110 * dx / soundspeed_wedge + 101 * dx / soundspeed_solid) With soundspeed_wedge = 2700 m/s, soundspeed_solid = 3000 m/s, dx = 0,0001 m/s t = 1.481e-5 = 14.81 µs Code (with the longer simulation time the simulation takes some hours): % ========================================================================= % 3D ULTRASONIC SIMULATION - BOREHOLE IN A FRP COMPONENT % CIRCULAR TRANSDUCER % ========================================================================= %% ======================================================================== % DEFINITION OF THE K-WAVE GRID % ========================================================================= % PML size (perfectly matched layer) on the outside of the k-wave grid PML_X_SIZE = 10; % [gridpoints] PML_Y_SIZE = 10; % [gridpoints] PML_Z_SIZE = 10; % [gridpoints] PML_SIZE = [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE]; % Computational grid: number of gridpoints without PML % The number of gridpoint (computational grid + PML) in each dimension % should be given by a power of 2 or have a small primefactor because of % computational advantages for FFT use. Nx = 256 - 2*PML_X_SIZE; % [gridpoints] Ny = 64 - 2*PML_Y_SIZE; % [gridpoints] Nz = 64 - 2*PML_Z_SIZE; % [gridpoints] % Distance between gridpoints dx = 1e-4; % [m] dy = 1e-3; % [m] dz = 1e-3; % [m] % Creation of the k-space grid kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); %% ======================================================================== % MEDIUM PARAMETER & TIME ARRAY % ========================================================================= % Medium properties c_gel = 1500; % [m/s] speed of sound gel rho_gel = 1002; % [kg/m^3] density gel c_air = 343; % [m/s] speed of sound air % The simulation is unstable for the correct density of air (1.2 [kg/m^3]) % because of the very large impedance contrast. To fix this problem the air % density is artificially inflated by factor 100. The reflection coefficient % is still close to 1. rho_air = 120; % [kg/m^3] density air cp_FRP = 3000; % [m/s] speed of sound FRP Compression (longitudinal) cs_FRP = 3000; % [m/s] speed of sound FRP Shear (transversal) rho_FRP = 1500; % [kg/m^3] density FRP cp_PMMA = 2700; % [m/s] speed of sound PMMA Compression (longitudinal) cs_PMMA = 1300; % [m/s] speed of sound PMMA Shear (transversal) rho_PMMA = 1180; % [kg/m^3] density PMMA % Medium model % Creating matrices with the dimension of the k-space grid for the density % and the sound speed of longitudinal and transversal waves medium.sound_speed_compression = cp_FRP * ones(Nx, Ny, Nz); % [m/s] medium.sound_speed_shear = cs_FRP * ones(Nx, Ny, Nz); % [m/s] medium.density = rho_FRP * ones(Nx, Ny, Nz); % [kg/m^3] % Air layer (on top of the plate) air_layer = 110; % [gridpoints] air layer thickness for m = 1 : air_layer medium.sound_speed_compression(m,:,:) = c_air; % [m/s] medium.sound_speed_shear(m,:,:) = c_air; % [m/s] medium.density(m,:,:) = rho_air; % [kg/m^3] end % PMMA wedge disc1 = makeDisc(Ny, Nz, Ny/2, Nz/2, 3); % disc position, radius wedge = repmat(disc1, [1, 1, Nx]); % 3D disc extrusion wedge = permute (wedge, [3,1,2]); % rearranging dimensions % length of wedge +1 because of later added gel layer on top of the wedge wedge_length = 110+1; % [gridpoints] borehole depth for n = (wedge_length+1) : Nx % PMMA model using 0 and 1 wedge(n,:,:) = 0; end % Filling the medium model with PMMA properties medium.sound_speed_compression(wedge==1) = cp_PMMA; % [m/s] medium.sound_speed_shear(wedge==1) = cs_PMMA; % [m/s] medium.density(wedge==1) = rho_PMMA; % [kg/m^3] % Defects in the FRP component (borehole) disc2 = makeDisc(Ny, Nz, Ny/2, Nz/2, 5); % disc position, radius flaw1 = repmat(disc2, [1, 1, Nx]); % 3D disc extrusion flaw1 = permute (flaw1, [3,1,2]); % rearranging dimensions hole_depth = 25; % [gridpoints] borehole depth for n = 1 : (Nx-hole_depth-1) % borehole model using 0 and 1 flaw1(n,:,:) = 0; end % Filling the medium model with air properties medium.sound_speed_compression(flaw1==1) = c_air; % [m/s] medium.sound_speed_shear(flaw1==1) = c_air; % [m/s] medium.density(flaw1==1) = rho_air; % [kg/m^3] % Two gel layers (on both sides of the wedge) % Filling the medium model with gel properties % Gel layer above the wedge medium.sound_speed_compression(1,:,:) = c_gel; % [m/s] medium.sound_speed_shear(1,:,:) = c_air; % [m/s] medium.density(1,:,:) = rho_gel; % [kg/m^3] % Gel layer under the wedge medium.sound_speed_compression((wedge_length+1),:,:) = c_gel; % [m/s] medium.sound_speed_shear((wedge_length+1),:,:) = c_air; % [m/s] medium.density((wedge_length+1),:,:) = rho_gel; % [kg/m^3] %% Time array dt = 1e-9; % [s] simulation timestep t_end = 2e-5; % [s] time end kgrid.t_array = 0:dt:t_end; % create time array % ========================================================================= % INPUT SIGNAL % ========================================================================= % Tone burst tone_burst_freq = 5e6; % [Hz] tone_burst_cycles = 3; % Time varying stress at each of the source positions source.sxx = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles); source.syy = 0; source.szz = 0; source.sxy = 0; source.sxz = 0; source.syz = 0; % Filter the source to remove high frequencies not supported by the grid source.sxx = filterTimeSeries(kgrid, medium, source.sxx); source.syy = filterTimeSeries(kgrid, medium, source.syy); source.szz = filterTimeSeries(kgrid, medium, source.szz); source.sxy = filterTimeSeries(kgrid, medium, source.sxy); source.sxz = filterTimeSeries(kgrid, medium, source.sxz); source.syz = filterTimeSeries(kgrid, medium, source.syz); % ========================================================================= % TRANSDUCER DEFINITION [SOURCE AND SENSOR] % ========================================================================= % Transducer model using 0 and 1 % Source source.s_mask = zeros(Nx, Ny, Nz); source.s_mask(1,:,:) = disc1; % same dimensions as the wedge % Sensor sensor.mask = source.s_mask; %% ======================================================================== % RUN THE SIMULATION % ========================================================================= % Assign additional input options % place PML outside of the computational grid input_args = { ... 'PMLInside', false, 'DisplayMask', sensor.mask, 'PlotPML', false, ... 'PMLSize', PML_SIZE}; % Run the simulation sensor_data = pstdElastic3D(kgrid, medium, source, sensor, input_args{:}); Bradley Treeby on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6497 Sun, 24 Jun 2018 19:19:33 +0000 Bradley Treeby 6497@http://www.k-wave.org/forum/ Hi Fabian, 1) The simulation is unstable because of the very large impedance contrast to air. You can solve this by reducing the size of the time step. Alternatively, I normally artificially inflate the density of the air by a factor of 100 or so. The reflection coefficient is still close to one, but the simulation should then run without any problems. 2) If you have a single transducer, you can just sum or average the data, i.e., output = mean(sensor_data, 1); 3) Depends on the situation - you can always try it with and without the layer to see if makes any difference. 4) Yes that's correct. Brad. fabwo on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6473 Fri, 25 May 2018 15:25:55 +0000 fabwo 6473@http://www.k-wave.org/forum/ Hi Brad and Ben, my goal is to simulate UT (TOFD) of CFRP Components. Because I am new to K-Wave I want to start with simulating the CFRP as homogenous material. I used pstdElastic3D to simulate propagation in a solid plate. The modelled defect in the CFRP plate is a simple borehole from beneath. Ultimately I want to model a phased array. For now the transducer is circular and directly centered above the borehole. I used a sinusoidal source. In this code I cut of the time steps after 61 because of lower computation time. 1) My sensor_data always oscillates to very hight values after a few time steps. I have no idea what the reason is. I tried to play with the PML and the absoption coefficients but did not receive better results. Can you help me with this? 2) I created the circular transducer with the makeDisc command. How can i put all sensor_data together to create an A-Scan at the transducer position? 3) I modelled a gel layer (sound speed and density) where the transducer is placed. Is this necessary? 4) When i tried to model a phased array I received a notice telling me that the kWaveTransducer class can not be used in pstdElastic. Is this correct? Thanks, Fabian 1 % ========================================================================= 2 % 3D ULTRASONIC SIMULATION BOREHOLE IN A CFRP COMPONENT 3 % CIRCULAR TRANSDUCER 4 % ========================================================================= 5 6 7 clearvars; 8 9 10 % ========================================================================= 11 % DEFINITION OF THE K-WAVE GRID 12 % ========================================================================= 13 14 % PML size (perfectly matched layer) 15 PML_X_SIZE = 20; % [gridpoints] 16 PML_Y_SIZE = 20; % [gridpoints] 17 PML_Z_SIZE = 20; % [gridpoints] 18 PML_SIZE = [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE]; 19 20 % Number of gridpoints without PML 21 Nx = 64 - 2*PML_X_SIZE; % [gridpoints] 22 Ny = 256 - 2*PML_Y_SIZE; % [gridpoints] 23 Nz = 256 - 2*PML_Z_SIZE; % [gridpoints] 24 25 % Distance between gridpoints 26 dx = 1e-4; % [m] 27 dy = 1e-4; % [m] 28 dz = 1e-4; % [m] 29 30 % K-pace grid 31 kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); 32 33 34 % ========================================================================= 35 % MEDIUM PARAMETER & TIME ARRAY 36 % ========================================================================= 37 38 % Medium properties 39 c_gel = 1500; % [m/s] speed of sound gel 40 rho_gel = 1200; % [kg/m^3] density gel 41 c_air = 343; % [m/s] speed of sound air 42 rho_air = 1.2; % [kg/m^3] density air 43 cp_cfrp = 3000; % [m/s] speed of sound CFRP Compression (longitudinal) 44 cs_cfrp = 3000; % [m/s] speed of sound CFRP Shear (transversal) 45 rho_cfrp = 1500; % [kg/m^3] density CFRP 46 47 medium.sound_speed_compression = cp_cfrp * ones(Nx, Ny, Nz); % [m/s] 48 medium.sound_speed_shear = cs_cfrp * ones(Nx, Ny, Nz); % [m/s] 49 medium.density = rho_cfrp * ones(Nx, Ny, Nz); % [kg/m^3] 50 51 % % Absorption properties 52 % medium.alpha_coeff_compression = 1; % [?] 53 % medium.alpha_coeff_shear = 1; % [?] 54 55 % Defects in the CFRP component (borehole) 56 disc1 = makeDisc(Ny, Nz, Ny/2, Nz/2, 5); % Disc position, radius 57 flaw1 = repmat(disc1, [1, 1, Nx]); % 3D disc extrusion 58 flaw1 = permute (flaw1, [3,1,2]); 59 hole_depth = 25; % [gridpoints] borehole depth 60 for n = 1 : (Nx-hole_depth-1) 61 flaw1(n,:,:) = 0; 62 end 63 medium.sound_speed_compression(flaw1==1) = c_air; % [m/s] 64 medium.sound_speed_shear(flaw1==1) = c_air; % [m/s] 65 medium.density(flaw1==1) = rho_air; % [kg/m^3] 66 67 % gel layer 68 medium.sound_speed_compression(1,:,:) = c_gel; % [m/s] 69 medium.sound_speed_shear(1,:,:) = c_gel; % [m/s] 70 medium.density(1,:,:) = rho_gel; % [kg/m^3] 71 72 % Time array 73 cfl = 0.1; % default = 0.3 74 t_end = 0.2e-6; % [s] 75 kgrid.makeTime(max(medium.sound_speed_compression,medium.sound_speed_shear), cfl, t_end); 76 77 78 % ========================================================================= 79 % INPUT SIGNAL 80 % ========================================================================= 81 82 % time varying sinusoidal source 83 source_freq = 1e6; % [Hz] 84 source_mag = 0.5; % [Pa] 85 source.sxx = source_mag * sin(2 * pi * source_freq * kgrid.t_array); 86 source.syy = 0; 87 source.szz = 0; 88 source.sxy = 0; 89 source.sxz = 0; 90 source.syz = 0; 91 92 % filter the source to remove high frequencies not supported by the grid 93 source.sxx = filterTimeSeries(kgrid, medium, source.sxx); 94 source.syy = filterTimeSeries(kgrid, medium, source.syy); 95 source.szz = filterTimeSeries(kgrid, medium, source.szz); 96 source.sxy = filterTimeSeries(kgrid, medium, source.sxy); 97 source.sxz = filterTimeSeries(kgrid, medium, source.sxz); 98 source.syz = filterTimeSeries(kgrid, medium, source.syz); 99 100 101 % ========================================================================= 102 % TRANSDUCER DEFINITION [SOURCE AND SENSOR] 103 % ========================================================================= 104 105 % Transducer element (source) 106 % source.s_mask: time varying stress source 107 % source.u_mask: time varying verlocity source 108 109 source.s_mask = zeros(Nx, Ny, Nz); 110 source.s_mask(1, Ny/2, Nz/2) = 1; 111 112 % Transducer element (sensor) 113 sensor.mask = source.s_mask; 114 115 116 117 % ========================================================================= 118 % RUN THE SIMULATION 119 % ========================================================================= 120 121 % assign the input options 122 input_args = { ... 123 'PMLInside', false, 'DisplayMask', sensor.mask, 'PlotPML', false, ... 124 'PMLSize', PML_SIZE}; 125 126 % run the simulation 127 sensor_data = pstdElastic3D(kgrid, medium, source, sensor, input_args{:}); 128 129 % size_sensor_data = size(sensor_data); 130 % unmasked_sensor_data = zeros(Ny, Nz, size_sensor_data(2)); 131 132 % for i = 1: size_sensor_data(2) 133 % unmasked_sensor_data(sensor.mask(1,:,:) ~= 0) = sensor_data(:,i); 134 % end