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