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