Excerpt from the Proceedings of the COMSOL Conference 2010 India

Two Dimensional FEM Simulation of Ultrasonic Wave Propagation in Isotropic Solid Media using COMSOL
Bikash Ghose1*, Krishnan Balasubramaniam2*, C V Krishnamurthy3, A Subhananda Rao1
1 2

High Energy Materials Research Laboratory, Sutarwadi, Pune, India PIN 411021 Center for Non-Destructive Evaluation, Dept of Mechanical Engg, IIT, Madras, Chennai - 36 3 Department of Physics, IIT Madras, Chennai 36
* Author for correspondence: [Link]@[Link], balas@[Link] Abstract: Ultrasonic wave propagation in solid can be modeled using Finite Element Method (FEM) that helps in understanding of the interaction of wave with material. The FEM uses various parameters which need to be optimized properly to obtain the solution closer to the exact one. A two dimensional FEM model is designed in COMSOL and has been simulated for ultrasonic wave propagation in an isotropic solid media to optimize the FEM parameters for getting closer solution. The change in the shape of incident wave, amplitude and frequency with respect to the change in length of element and time steps has been studied. It was observed that free meshing provided in COMSOL with triangular element is capable of simulating the wave propagation for the ratio of wavelength to the element length of more than 8. The time steps required for the simulation to obtain proper result should be better than that is merely satisfied by CFL criteria. Keywords: Ultrasonic, CFL Criteria, Element length, time steps, FEM COMSOL. The ultrasonic wave characteristics like shape of wave, frequency content, and wave amplitude of the onward and reflected propagated wave in material with respect to important simulation parameters like time step, element length has been studied for getting proper solution.

2. Two Dimensional FEM Model of ultrasonic wave propagation

A two-dimensional finite element model for propagation of ultrasonic wave in an isotropic solid media has been studied with COMSOL. 2.1 Material and ultrasonic wave properties The material chosen was steel with following properties Youngs Modulus (E) = 2 1011 Pa Poissons ratio () = 0.33 Density () = 7850 kg/m3 Ultrasonic velocity (CL) (longitudinal wave) = 5850 m/s Wavelength of the longitudinal ultrasonic wave (L) = CL/20103 s-1 = 0.2925 m 2.2 Incident ultrasonic wave A 20 kHz frequency signal has been applied to the isotropic material. The signal is chosen as the 3 cycles of cosine function and operated with a hanning window. The applied disturbance as the input signal is applied to a line of length 0.04 m which means the transducer has been modeled as the disturbances on a line. The initial disturbances is given as the displacements in the ve y direction only for modeling the normal incidence of ultrasonic wave. The input signal is provided from a data file and the linear interpolation is considered for the value of displacement at the undefined time in the data file. The time domain representation of input displacement pulse used in the model is shown in the figure 1. The frequency content of the input signal was evaluated as the Fourier transform of the time domain signal and is as shown in figure 2.

1. Introduction
Ultrasonic Testing (UT) is one of the important Non-Destructive Evaluation (NDE) technique widely used for characterisation of materials as well as detection and characterisation of flaws present in the material used in various industries. Understanding of ultrasonic wave propagation and its capability for flaw detection is an important aspect for proper evaluation of wave characteristics for characterisation of material or detection of flaws. The ultrasonic wave propagation can be simulated using Finite Element Method (FEM). The simulation parameters used in FEM plays very important role for the correctness of the obtained data. This paper presents the Finite Element Method (FEM) simulation of ultrasonic wave propagation in isotropic solid media using

8 6 4 Displacem ent (m ) 2 0 -2 -4 -6 -8

x 10

-7

-10 0

0.5 Time (s)

1.5 -4 x 10

Figure 1: Input displacement pulse as source excitation

x 10 5
-7

Single-Sided Amplitude Spectrum of y(t)(Hanning windowed signal)

2.5 Application Mode For modeling of the two dimensional geometry, the case of plane strain was considered because of the dimension of the material considered to be large enough in the third dimension as compared to the x and y directions. Time dependent analysis, Lagrange quadratic type of element and time dependent solver was used for the solution. The duration of time span for the solution was chosen such that back-wall reflected signal comes back to the front wall or from where excitation started. All the simulations were carried out in the 64-bit XP environment.

4.5

3.5

3 |Y(f)|

3. Simulation Results for Effect of Maximum Length of Element and Time Steps
3.1 Effect of Maximum Length of Element The length of element in the mesh used for solution of any FEM model plays a crucial role for correctness of the obtained solution. The length of the element needs to be smaller for evaluation of proper solution whereas as the length of element decreases the cost of computation increases. So it is always necessary to evaluate the optimum element length for correct solution along with the lesser computational efforts. In the current study, the maximum element size (xmax) was varied from L /2 to L /16 and the model was simulated for the ultrasonic wave propagation. L is the wavelength of the longitudinal ultrasonic wave propagated in the direction of incidence. Following table shows the L different values of for which the
xmax

2.5

1.5

0.5

10 Frequency (Hz)

12

14

16

18

20 x 10
4

Figure 2: Input signal in frequency domain

2.3 Model Geometry As the study was mainly focused on the bulk longitudinal wave propagation, a simple rectangular geometry has been chosen and the dimension of the geometry was chosen such that no side wall reflection reaches to the point of observation within the specified time. The initial disturbance was applied on a line of length 0.04 m located at the centre of the top horizontal line. For study of point sources as the model for the transducer, nine numbers of equally spaced points of initial disturbances are located within the line length of 0.04 m which was earlier considered as the line source. As all the observations were made for the longitudinal wave only the signal was observed on the perpendicular line that passes through the middle of the source. 2.4 Meshing The mesh was generated automatically with the triangular elements. Free meshing provided in COMSOL was used for the generation of mesh. The elements were chosen as uniform through out the sub-domain. The mesh has been refined at the region below the excitation source. Minimum mesh quality that has been ensured for each of the model that has been solved is 0.7 with an element area ratio of about 0.1.

simulation was carried out along with the corresponding maximum element size for each case. The time steps chosen initially for simulation for different xmax is 2.510-6 s. The reason behind Ratio of Wavelength to xmax (L /xmax ) 2 3 4 5 8 9 12 16 Maximum element size (xmax) (m) 0.1462 0.0975 0.0731 0.0585 0.0366 0.0325 0.0244 0.0182

Displacement (m)

the chosen time step is as follows. As per CFL criteria, the critical time steps to be used for simulation on FEM model for a time dependent solver is x , where Cph is the
C ph

4.00E-07 3.00E-07 2.00E-07 1.00E-07 0.00E+00 0 -1.00E-07 -2.00E-07 -3.00E-07 -4.00E-07 Time (s) 0.0001 0.0002 0.0003 0.0004

Ratio=2 Ratio=3 Ratio=4 Ratio=5 Ratio=8 Ratio = 9 Ratio = 12 Ratio = 16 0.0005

ultrasonic phase velocity. As only longitudinal ultrasonic wave was considered for this study, the Cph value was taken as the longitudinal ultrasonic velocity in the steel material which is about 5850 m/s. In this regard, for different xmax values different t should be evaluated which is to be used for solution. The critical time t would be smallest for the smallest value of xmax. If the smallest xmax was chosen for evaluating t then the t must be sufficiently smaller for other cases. Hence t was evaluated for the least value of xmax . The least value of xmax is 0.0182 m and correspondingly the critical value of t was calculated to be 0.0182m x x tcritical = max = max = = 3.1 10 6 s C ph CL 5850m / s
The value of t used in the simulation was 2.510-6 s which is less than the critical time steps for the case of L /xmax =16. The ultrasonic signals were plotted at different perpendicular distances from the middle of line source for different values of L at t =
xmax

(a) At 0.1 m
1.00E-07 8.00E-08 6.00E-08 Displacements (m) 4.00E-08 2.00E-08 0.00E+00 -2.00E-08 -4.00E-08 -6.00E-08 -8.00E-08 Time (s) 0 0.0001 0.0002 0.0003 0.0004 0.0005 Ratio = 2 Ratio = 3 Ratio = 4 Ratio= 5 Ratio = 8 Ratio = 9 Ratio = 12 Ratio = 16

(b) At 0.3 m
4.00E-08 3.00E-08 2.00E-08 Displacement 1.00E-08 0.00E+00 0 -1.00E-08 -2.00E-08 -3.00E-08 -4.00E-08 Time (t) 0.00005 0.0001 0.00015 0.0002 0.00025 0.0003 0.00035 0.0004 0.00045 Ratio = 2 Ratio = 3 Ratio = 4 Ratio = 5 Ratio = 8 Ratio = 9 Ratio = 12 Ratio = 16

410-4 s. A cross sectional plot (line profile for the plot of displacements Vs distances or position) was also plotted at t = 410-4 s. The time t = 410-4 s is chosen for the above said plot as by that time the signal just reaches to the back-wall. The above plots were used to evaluate the minimum value of L required for the
xmax

(c) At 1.0 m
2.50E-08 2.00E-08 1.50E-08 Displacement (m) 1.00E-08 5.00E-09 0.00E+00 -5.00E-09 -1.00E-08 -1.50E-08 -2.00E-08 -2.50E-08 Distance from Source center (m) 0 0.5 1 1.5 2 2.5 3 Ratio = 2 Ratio = 3 ratio = 4 ratio = 5 Ratio = 8 Ratio = 9 Ratio = 12

solution to converge irrespective of time steps if the solution remains consistent for the various element lengths. The signal at a various perpendicular distances from the middle of the line source was plotted for different L /xmax as shown below. The following plots (figure 3) show the forward propagating signal at a time of 410-4 s plotted along the perpendicular line passing through the centre of the line source.
Figure 3: Onward propagating wave signal at various distances for different wavelength to element length ratio

Figure 4: Line profile (Displacement Vs Distance from source) at 0.0004s for different wavelength to element length ratio

The plot indicates that the solution converges for L /xmax = 8 or more for constant time steps which is sufficient for L /xmax = 16 as per CFL criteria. In each of the above cases, the time step was taken as t = 2.5 10-6 s. The t taken is sufficient as per CFL criteria i.e., t (x/CL). xmax for condition of L /xmax = 16 is 0.0182 which accounts for the t as x/CL = 0.0182m/5850 m/s =3.110-6 s. The plot indicates that oscillations still persists after complete passing out of the signal which may be because of the lesser x. This plot also indicates the converging of the solution for the ratio (L /xmax) 8. It has been observed that whether the time steps are taken from solver or exactly what has been given does not make any difference to the propagating signal. 3.2 Effect of time steps All of the above solutions have been obtained for the input signal of 20 kHz. The sampling of the signal was done at 1 MHz and the time steps used was 2.510-6 s for all of the above solution. As per CFL criteria the time steps should be less than x/CL. As per the criteria, for example, for the maximum element length of L /10 the time steps should be (L /10CL) that means if
10 x L t max = 10 C L CL 1 10 f In the above cases L varied from 2 to 16 xmax t xmax =

media. The amplitude was expected to change consistently but at least not the shape of the input signal. In this regard, it may be thought that probably time steps used so far was not sufficient for a correct solution although the CFL criteria was satisfied well for all of the above cases. In this regard, the time steps used were further reduced for three values of L to obtain the
xmax

solution and checked for the convergence of solution as per the expected result. Case I: L =8 For the case of

L = 8, the solution is xmax obtained for the following time steps (i) As per CFL criteria x 0.0366m t = = 6.26 10 6 s CL 5850m / s (ii) t = 5 106 s (iii) t = 2 106 s (iv) t = 2.5 106 s (v) t = 1.0 106 s (vi) t = 0.5 106 s (vii) t = 0.2 106 s

xmax

The following plots (figure 5) show the signal at a perpendicular distance of 1.0 m from the middle of the line source for different time steps. For all the following cases the ratio L is taken
xmax

as 8.
Figure 5: Time domain signal at different time steps for L = 8 at 1.0 m
xmax

whereas in all the cases time step was taken 1 as t = 1 = = 2.5 10 6 s which 3 1
20 f 20 20 10 s

was considerably lesser than required as per the criteria. As already been observed, irrespective of the better time steps the solution converges only for L 8 i.e. there is no more change in the
xmax

solution beyond L = 8.
xmax

Unfortunately, although there is no change in the signal for any further reduction of element length but the signal so far obtained is not as per expected one. The ultrasonic wave form should not change after propagated through an isotropic

(a) Original applied signal at source

(b) t = 6.26 106 s (Critical time steps as per CFL Criteria)

(e) t = 2.0 10 s

(c)

t = 5 106 s

(f) t = 1.0 106 s

(d) t = 2.5 10 6 s

(g) t

= 0.5 106 s

(h) t = 0.2 10 6 s

Case II: (i)

(b) t = 10 106 s (Critical time steps as per CFL Criteria)

xmax
As

=5

per CFL criteria x 0.0244m t = = 10 10 6 s C L 5850m / s

(ii) t = 5 106 s (iii) t = 2.5 10 6 s (iv) t = 2 106 s (v) t = 1 106 s (vi) t = 0.5 10 6 s (vii) t = 0.2 106 s The following plots (figure 6) show the signal at a perpendicular distance of 1.0 m from the middle of the line source for different time steps. For all the following cases the ratio L is taken
xmax

(c)

t = 5 106 s

as 5.
Figure 6: Time domain signal at different time steps for

xmax

L = 5 at 1.0 m

(d) t = 2.5 10 6 s

(a) Original applied signal at source

(h) t = 0.2 10 6 s

Case III: (i)

xmax

= 12

As

per CFL criteria x 0.0244m t = = 4.17 10 6 s 5 10 6 s CL 5850m / s

(e) t = 2.0 10 s

(ii) t = 2.5 10 6 s (iii) t = 1 106 s (iv) t = 0.5 10 6 s The plots in figure 7 show the signal at a perpendicular distance of 1.0 m from the middle of the line source for different time steps. For all the following cases the ratio L is taken as 12. xmax
Figure 7: Time domain signal at different time steps for

xmax

L = 12 at 1.0 m

(f) t = 1.0 106 s

(a) t = 5 106 s

(g) t = 0.5 106 s

(b) t = 2.5 10 6 s

(c) t = 1.0 106 s

(b) L = 8 t = 0.5 10 6 s xmax

(d) t = 0.5 106 s To compare the plots for convergence of the solution in terms of time step of 0.510-6 s, the total signal is shown below for three different L values.

(c) L = 5 t = 0.5 10 6 s xmax

The signal again indicates that the optimum value of the ratio L = 8 . More ripples in the signal are seen for higher xmax in spite of very less t value. 4. Effect of line and point source excitations In all of the earlier cases the source was modeled as the line source which gets excitation as a whole. It was also possible to model the source as per Huygens principle i.e., considering the source as many secondary sources of excitations. Few of the simulations were carried out by considering point sources as the detector. In that case, within total 0.04 m length of source, nine numbers of equidistant source excitation points were placed. The simulation was done for line as L well as point sources for =8
xmax
xmax

xmax

Figure 8: Comparisons between converged solution at

xmax

L = 8, 12 and 5 at t = 0.5 106 s

(a)

xmax

6 = 12 t = 0.5 10 s

t = 2.5 10 s . The following plots in figure 9 show the ultrasonic signal for both cases (i) for line source (ii) for point sources observed at different perpendicular distances from middle of source.

The line profile (displacement Vs distance from source) was also plotted at a time instant of 410-4 s.
Figure 9: Comparison between source with excitation on line and excitation on points
L
xmax =8

L
xmax

=8

t = 2.5 106 s (point source)

At 1.0 m

t = 2.5 106 s (line source)

At 1.0 m

At 0.5 m

At 0.5 m

Line profile at t = 410-4 s

Line profile at t = 410-4 s

5. Frequency content of the forward propagated signal

The frequency content of the onward propagating wave has been evaluated by Fourier transform of the onward propagating wave. The figure 10 shows the onward propagating wave in

time and frequency domain for three cases of 6 L = 12, 8 and 5 for t = 0.5 10 s
xmax

Frequency domain
14 12 10 8 6 4 2 x 10
-9

Single-Sided Amplitude Spectrum

(converged final solution for each case)

Figure 10: Time and Frequency domain signal for onward propagating wave
L
xmax = 12

t = 0.5 106 s

Time domain
8.00E-08

|Y(f)|

4 5 6 Frequency (Hz)

9 x 10
4

4.00E-08 Displacement (m)

L
xmax
2.00E-04 4.00E-04 6.00E-04

=5

t = 0.5 10 s

0.00E+00 0.00E+00 -4.00E-08

Time domain
8.00E-08

-8.00E-08
4.00E-08 Displacement (m)

-1.20E-07 Time (s)

Frequency domain
10 9 8 7 6
|Y(f)|

0.00E+00 0.00E+00 -4.00E-08

2.00E-04

4.00E-04

6.00E-04

x 10

-9

Single-Sided Amplitude Spectrum

-8.00E-08

-1.20E-07 Time (s)

5 4 3 2 1 0 1 2 3 4 5
Frequency (Hz)

Frequency domain

x 10
6 7 8 9 x 10
4

-9

Single-Sided Amplitude Spectrum

L
xmax

12 10 8 6 4 2

=8

t = 0.5 10 s
|Y(f)|

Time domain
8.00E-08

4.00E-08 Displacement (m)

0.00E+00 0.00E+00 -4.00E-08

4 5 6 Frequency (Hz)

9 x 10
4

2.00E-04

4.00E-04

6.00E-04

-8.00E-08

-1.20E-07 Time (s)

6. Frequency content of the back-wall reflected signal The figure 11 shows the back-wall reflected wave in time and frequency domain for three cases of L 12, 8 and 5 for t = 0.5 106 s
xmax =

(converged final solution for each case)

Figure 11: Time and Frequency domain signal for reflected wave

Frequency domain
x 10 8 7 6 5 |Y(f)| 4 3 2
-9

Single-Sided Amplitude Spectrum

L
xmax

= 12

t = 0.5 106 s

Time domain
4.00E-08

2.00E-08 Displacement (m)

1

0.00E+00 0.00E+00 -2.00E-08

4.00E-04

8.00E-04

1.20E-03

4 5 6 Frequency (Hz)

9 x 10
4

-4.00E-08

L
xmax
Time (s)

=5

t = 0.5 10 6 s

-6.00E-08

Time domain
4.00E-08

Frequency domain
Displacement (m)
7 6 5 4 3 2 1 x 10
-9

Single-Sided Amplitude Spectrum

2.00E-08

0.00E+00 0.00E+00 -2.00E-08

4.00E-04

8.00E-04

1.20E-03

|Y(f)|

-4.00E-08

-6.00E-08 Time (s)

Frequency domain
1 2 3 4 5 6 Frequency (Hz) 7 8 9 x 10
4

x 10 8

-9

Single-Sided Amplitude Spectrum

L
xmax

=8

t = 0.5 106 s
|Y(f)|

7 6 5 4 3 2

Time domain
4.00E-08

2.00E-08
1

Displacement (m)

0.00E+00 0.00E+00 -2.00E-08

4.00E-04

8.00E-04

1.20E-03

4 5 6 Frequency (Hz)

9 x 10
4

-4.00E-08

-6.00E-08 Time (s)

7. Discussion All results shown in figure 3 were obtained for the time steps of 2.510-6 s which is well below the value calculated using CFL criteria. The time domain signal shown in figure 3 (a) indicates that for all the cases of L /xmax , the signal follows the same time profile path (no change in time scale) till the major signal passes through the point of observation. For all the cases, after the major signal passes through, considerable amount of the residual displacement are still shown to be exists. At that region, the

curve for L /xmax = 2, 3 &4 does not follow to curves for L /xmax = 5, 8, 9, 12 &16. There is minor variation seen for the case of 5 with others. For the main signal, although there is no change in terms of time scale but it is seen that the amplitude is less in case of higher value of L /xmax . When the same signal is observed at a further distance it is seen from figure 3(b) that there is change in the shape of wave form that was propagated from the earlier point. From figure 3(a) to 3(c), it is observed that the solution for the major signal always follows for all L /xmax whereas at the region of residual displacement, the solution getting converged only from L /xma = 5 onwards. But in spite of very poor value of L /xmax = 2, 3 & 4, the solution evolved correctly in terms of time for all L /xmax . Again for all plots under figure 3, the signal is not of the shape as expected. For a nondispersive media like steel, the ultrasonic signal must not change its shape during propagation. This indicates that in spite of following the CFL criteria also the signal could not really converge for the amplitude details. Figure 4 shows that for L /xmax = 2, 3 & 4 there is considerable error in the residual signal (which remains after the actual signal pass through the point of observation). Figure 5 indicates that the time steps used earlier for obtaining the solution was not sufficient although it followed the CFL criteria. Because of change in time steps, there is a change in both frequency information as well as intensity information of the propagating signal. At the time step of 0.510-6 s only, the propagated signal resumes the shape of original signal. Further decrease in the time step below 0.510-6 s does not make any further improvement of the signal which indicates the possible convergence of the solution. Similarly the solution converges for the case of L /xmax = 8 only for time steps of 0.510-6 s. Even for L /xmax = 5 also the solution converges with the time step of 0.510-6 s. This clearly indicates that time step is a very important parameter for getting correct solution for wave propagation. For very poor value of L /xmax like 2 or 3 can at least give some information regarding the signal with a better value of t but any good value of x is of no use if a proper value of t is not used for obtaining the solution.

At least in three cases it is seen that the required time steps is about (1/100*f) which is independent of the value of x. The value of x has to be taken independently at about /8. From figure 9 it is seen that at least for the normal incidence of ultrasonic beam no difference in terms of shape of the signal is observed between the solutions obtained for two ways of excitation viz. excitation on a line and excitation on points to resemble the source. The only difference observed is in terms of the intensity. Intensity observed in case of point sources as the excitations is lesser than that observed in case of line source. The observation is obvious as power input in the case of line source will always be more than that from the point sources as the excitations. Figure 10 shows that there is no substantial change in the frequency domain signal for the onward propagating wave for L /xmax = 8, 12 & 5 in the final converged solution. Similarly no difference is observed in the frequency domain reflected signal. However a little difference in intensity is observed in both onward and reflected wave. Intensity is higher for L /xmax = 5 and least for L /xmax = 12.

8. Conclusions Ultrasonic wave propagation can well be modeled with COMSOL with excitation on a line segment as transducer or excitations on points on a line segment. The triangular element free meshing can be used for the simulation of ultrasonic wave propagation for a value of L /xmax = 8 or more. The time steps should be used near to 1/100*f. Smaller value of L /xmax such as 5 can also be used with expected smaller variation of amplitude information with use of proper time steps. No substantial difference in the frequency content of the onward as well as back-reflected converged solution is observed.