Sound Thesis
Sound Thesis
University of London
by
January 2008
Friedrich Nietzsche
2
Abstract
The aim of the work presented in this thesis is to provide tools to extend modelling
capacities and improve quality and reliability of bulk and guided wave propagation
models using commercially available finite element (FE) packages.
During the development process of NDT inspection techniques, the knowledge of the
interaction of waves with defects is key to the achievement of robust and efficient
techniques as well as identifying potential weaknesses. The reflection of ultrasound
from cracks and notches of simple geometry and orientation is already well
understood, but there are few results for more complex cases. A discrete approach is
needed to model how the waves interact with discontinuities, including structural
features, cracks, corrosion or other forms of defects. FE methods have been used to
model a wide range of bulk and guided waves problems and have successfully
provided important information about wave interaction with discontinuities. In these
studies, defects were strongly simplified. One reason for this is that initial work is
bound to focus on the simplest cases, but many modellers are ready to go on to more
complex problems. The reason that so little of that is happening is that, despite rapid
growth in computer power, many of the more complex realistic problems are still
beyond the capacity of the models. The more complex problems require much larger
models than the simplified ones, and so have remained out of reach.
This can be changed by using innovative techniques and improving the quality and
reliability of modelling by taking the right decisions during the modelling process.
Perfectly matched layers (PML) and absorbing layers using increasing damping
(ALID) enabling a reduction in the model geometric size are implemented in
commercially available FE packages. Analytical models are developed in order to
facilitate the achievement of high computational efficiency. Demonstrator cases
highlight the gains achieved by the use of these techniques.
As the choice of mesh density is crucial in defining the resources necessary to solve a
model, a study of the influence of meshing parameters for various element types and
numerical schemes on the propagation velocity is performed. This provides
information helping modellers to reach the right modelling compromises thanks to an
improved understanding of the consequences of the decisions made.
3
The accuracy of defect modelling is investigated for a range of situations and
modelling strategy. The weight of the choice of the right strategies is demonstrated.
The outcome of the use of the techniques and information presented in this thesis is a
significant improvement in FE modelling of waves in elastic media.
4
Acknowledgements
Firstly, I would like to thank my supervisors Professor Mike Lowe and Professor Peter
Cawley for their guidance, help and patience during the course of this project. I would
also like to express my gratitude for giving me the opportunity to be part of the NDT
laboratory at Imperial College.
My thanks extend to all my colleagues in the NDT lab. It was a pleasure to be part of
the team.
I would like to thank my family and my wife’s family for their support.
This work was funded by the EPSRC and Rolls Royce plc.
5
Table of contents
Chapter 1
Introduction
1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2. Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3. Outline of the thesis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
Chapter 2
Theoretical background
1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2. Theory of wave propagation in elastic media . . . . . . . . . . . . . . . . . . . . . . . 24
2.1 Bulk wave propagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.1.1 Bulk wave propagation in infinite isotropic elastic media. . . . . . 24
2.1.2 Bulk propagation in a semi-infinite isotropic elastic medium . . . 26
2.2 Guided wave propagation in a plate . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3. Finite elements modelling of wave propagation . . . . . . . . . . . . . . . . . . . . . 34
3.1 Explicit method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.2 Implicit method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2.1 ABAQUS/Standard procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.2.2 COMSOL Multiphysics procedure . . . . . . . . . . . . . . . . . . . . . . . 40
4. Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2. Review of non-investigated techniques. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.1 Infinite element methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.2 Non reflecting boundary condition . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3. Absorbing layer theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.1 Concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.2 Perfectly matched layer (PML) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.3 Absorbing layer using increasing damping (ALID) . . . . . . . . . . . . . . 50
4. Efficient layer parameters’ definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.1 Analytical model for bulk waves. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.1.1 General definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.1.2 Validation procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.1.3 PML analytical model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.1.4 ALID analytical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.2 Analytical model for 2D guided wave cases . . . . . . . . . . . . . . . . . . . . 66
4.2.1 Consideration for guided wave PML implementation . . . . . . . . 66
4.2.2 Validation procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.2.3 PML analytical model for guided wave cases . . . . . . . . . . . . . . 68
4.2.4 ALID analytical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5. Demonstrators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.1 Computational efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.1.1 Bulk wave model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.1.2 Guided wave model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6
Table of contents
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
2. Explicit solving . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
2.2 Linear quadrilateral elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
2.2.1 Square elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
2.2.2 Rectangle elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
2.2.3 Rhombus elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
2.2.4 Parallelogram elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
2.2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
2.3 Linear triangular elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
2.3.1 Equilateral triangle elements. . . . . . . . . . . . . . . . . . . . . . . . . . . 106
2.3.2 Isosceles triangle elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
2.3.3 Scalene triangle elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
2.3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
2.4 Modified quadratic triangular elements. . . . . . . . . . . . . . . . . . . . . . . 114
2.4.1 Equilateral triangle elements. . . . . . . . . . . . . . . . . . . . . . . . . . . 114
2.4.2 Isosceles triangle elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
2.4.3 Scalene triangle elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
2.4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
3. Implicit solving . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
3.2 Linear quadrilateral elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
3.2.1 Square elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
3.2.2 Rectangle elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
3.2.3 Rhombus elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
3.2.4 Parallelogram elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
3.2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
3.3 Quadratic quadrilateral elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
3.3.1 Square elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
3.3.2 Rectangle elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
3.3.3 Rhombus elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
3.3.4 Parallelogram elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
3.3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
3.4 Linear triangular elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
3.4.1 Equilateral triangle elements. . . . . . . . . . . . . . . . . . . . . . . . . . . 133
3.4.2 Isosceles triangle elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
3.4.3 Scalene triangle elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
3.4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
3.5 Quadratic triangular elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
3.5.1 Equilateral triangle elements. . . . . . . . . . . . . . . . . . . . . . . . . . . 136
3.5.2 Isosceles triangle elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
3.5.3 Scalene triangle elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
7
Table of contents
Chapter 5
Accurate modelling of defects using Finite Elements
1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
2. Model definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
3. Reflection from a straight edge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
4. Reflection from a straight crack at an angle . . . . . . . . . . . . . . . . . . . . . . . 156
4.1 Crack of unit length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
4.2 Crack of length 0.25 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
4.3 Crack of length 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
5. Reflection from circular defects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
5.1 Hole of unit diameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
5.2 Hole of diameter 0.25 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
5.3 Hole of diameter 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
5.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178
6. Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
Chapter 6
Local mesh refinement
1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
2. Fictitious domain technique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
2.1 Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
2.2 Presentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
2.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183
3. Abrupt mesh density variation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183
3.1 1D wave propagation models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183
3.1.1 Model definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183
3.1.2 L wave 1D model using theoretical material properties . . . . . . 185
3.1.3 L wave 1D model with matched acoustic impedance. . . . . . . . 186
3.1.4 L and S wave 1D model with matched acoustic impedance. . . 187
3.1.5 L and S wave 1D model with varying acoustic impedance . . . 190
3.1.6 L wave 1D model with different mesh ratio . . . . . . . . . . . . . . . 193
3.2 2D wave propagation models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195
3.2.1 Model definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195
4. Gradual mesh density variation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199
4.1 1D wave propagation model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199
4.2 2D wave propagation model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201
5. Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202
8
Table of contents
Chapter 7
Conclusions
1. Review of thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 204
2. Summary of findings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205
2.1 Absorbing layers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205
2.2 Influence of mesh parameters on the elastic bulk wave velocities . . 207
2.3 Accurate modelling of complex defects using Finite Elements . . . . 208
2.4 Local mesh refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 209
3. Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 210
3.1 Absorbing layers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 210
3.2 Influence of mesh parameters on the elastic bulk wave velocities . . 210
3.3 Accurate modelling of complex defects using Finite Elements . . . . 211
3.4 Local mesh refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211
References
9
List of figures
Figure 1.1 a) 2D plane strain model of a plate including a defect, b) Time signal at
the monitoring point.............................................................................. 19
Figure 3.1 a) 2D plane strain model of a plate including a defect, b) Time signal at
the monitoring point.............................................................................. 42
Figure 3.2 Illustration of use of infinite elements .................................................. 44
Figure 3.3 ABAQUS benchmark model: a) Model geometry, b) vertical displace-
ment at point A, Extended model (reference): c) Model geometry, d) ver-
tical displacement at point A ................................................................ 45
Figure 3.4 Absorbing layer concept for 2D models: a) infinite medium, b) semi in-
finite medium, c) plate .......................................................................... 46
Figure 3.5 Variation of αx(x) and αy(y) in a 2D model ........................................ 48
Figure 3.6 Spatial spread of the reflection and transmission for a single layer (no
mode conversion shown for simplicity)................................................ 54
Figure 3.7 Illustration of extreme angles defining the range of angles to consider
when dimensioning an absorbing layer ................................................ 54
Figure 3.8 FE model used to validate the analytical models a) normal incidence
model, b) angled incidence model ........................................................ 56
Figure 3.9 a) Reflection coefficient against αx b) Reflection coefficient against the
number of elements per wavelength ..................................................... 57
Figure 3.10 Reflection coefficient for a given PML obtained with bulk wave analyti-
cal and FE models................................................................................. 59
Figure 3.11 Reflection coefficient for a given ALID obtained with bulk wave analyt-
ical and FE models................................................................................ 64
Figure 3.12 FE model used to validate the guided wave analytical models ............ 67
Figure 3.13 Reflection coefficient for a given PML obtained with guided wave ana-
lytical and FE models ........................................................................... 70
Figure 3.14 Definition of the multi layered system ................................................. 71
Figure 3.15 Reflection coefficient for a given ALID obtained with guided wave ana-
lytical and FE models ........................................................................... 76
Figure 3.16 a) bulk wave demonstrator, FE model: b) without absorbing layer, c) with
ALID, d) with PML .............................................................................. 77
10
List of figures
Figure 3.17 Absolute displacement field for the bulk demonstrator with ALID at
time: a)5msec b)10msec c)15msec d)20msec. Colour scale extends from
0 (blue) to 0.1% (red) of the maximum absolute displacement. Grey indi-
cates out of scale (0.1% to 100%). White dashed line indicates the bound-
ary between area of study and ALID .................................................... 78
Figure 3.18 a) guided wave demonstrator, FE model: b) without absorbing layer, c)
with ALID, d) with PML ...................................................................... 79
Figure 3.19 Absolute displacement field for the guided demonstrator with ALID at
time: a)150msec b)300msec c)450msec d)600msec. Colour scale is var-
ied and extends from 0 (blue) to 2% or 10% (red) of the maximum abso-
lute displacement as indicated on the figure. Grey indicates out of scale
(2% or 10% to 100%). White dashed line indicates the boundary between
area of study and ALID ........................................................................ 80
Figure 3.20 Input preprocessing............................................................................... 81
Figure 3.21 Model geometry for time reconstruction case ...................................... 81
Figure 3.22 Normal displacement monitored 700mm away from the defect. a) Classi-
cal time domain analysis with ABAQUS, b) Frequency domain analysis
with ABAQUS, c) Frequency domain analysis with COMSOL .......... 82
Figure 3.23 a) dispersion curve data used for input definition, b) input definition . 82
Figure 3.24 Representation of model used for guided wave scattering validation .. 83
Figure 3.25 Example of a typical spatial FFT curve ................................................ 83
Figure 3.26 Reflection coefficient against notch width ........................................... 84
Figure 3.27 Energy reflection coefficient for A0 incident on a 2mm square notch in
an 8mm thick aluminium plate from 140kHz to 500kHz ..................... 86
11
List of figures
Figure 4.12 a) Shape of the different rhombic elements used in the mesh; Variation of
the longitudinal (b and d) and shear (c and e) velocity error against the
angle of incidence for various shearing angle g plotted in a polar (b and
c) and linear (d and e) fashion. The coloured circles indicate the error pre-
diction along the element side and diagonal....................................... 102
Figure 4.13 a) Shape of the different parallelogramatic elements used in the mesh;
Variation of the longitudinal (b and d) and shear (c and e) velocity error
against the angle of incidence for various shearing angle g plotted in a po-
lar (b and c) and linear (d and e) fashion. The coloured circles indicate the
error prediction along the element side and diagonal ......................... 104
Figure 4.14 Schematic defining L0, L90, L30 and Lq in a mesh of equilateral-trian-
gular elements ..................................................................................... 105
Figure 4.15 Variation of the longitudinal (a and c) and shear (b and d) velocity error
against the angle of incidence for various mesh densities plotted in a lin-
ear (a and b) and polar (c and d) fashion ............................................ 106
Figure 4.16 Velocity error against mesh density for shear and longitudinal waves at 0
and 30 degrees .................................................................................... 107
Figure 4.17 Velocity errors against CFLX for various mesh densities at a) 0 and b) 30
degrees ................................................................................................ 108
Figure 4.18 a) Shape of the different isosceles-triangular elements used in the mesh;
Variation of the longitudinal (b and d) and shear (c and e) velocity error
against the angle of incidence for various values of f plotted in a polar (b
and c) and linear (d and e) fashion. The coloured circles indicate the error
prediction along the element side and diagonal.................................. 110
Figure 4.19 a) Shape of the different scalene-triangular elements used in the mesh;
Variation of the longitudinal (b and d) and shear (c and e) velocity error
against the angle of incidence for various values of g plotted in a polar (b
and c) and linear (d and e) fashion. The coloured circles indicate the error
prediction along the element side and diagonal.................................. 112
Figure 4.20 Variation of the longitudinal (a and c) and shear (b and d) velocity error
against the angle of incidence for various mesh densities plotted in a polar
(a and b) and linear (c and d) fashion ................................................. 114
Figure 4.21 Schematic defining L0, L90, L30 and Lq in a mesh of quadratic equilat-
eral-triangular elements ...................................................................... 114
Figure 4.22 Velocity error against mesh density for shear and longitudinal waves at 0
and 30 degrees .................................................................................... 115
Figure 4.23 Velocity errors against CFLX for various mesh densities at a) 0 and b) 30
degrees ................................................................................................ 117
Figure 4.24 a) Shape of the different quadratic isosceles-triangular elements used in
the mesh; Variation of the longitudinal (b and d) and shear (c and e) ve-
locity error against the angle of incidence for various values of f plotted
in a polar (b and c) and linear (d and e) fashion. The coloured circles in-
dicate the error prediction along the element side and diagonal......... 118
Figure 4.25 a) Shape of the different quadratic scalene-triangular elements used in the
mesh; Variation of the longitudinal (b and d) and shear (c and e) velocity
error against the angle of incidence for various values of g plotted in a po-
lar (b and c) and linear (d and e) fashion. The coloured circles indicate the
error prediction along the element side and diagonal ......................... 120
Figure 4.26 Variation of the longitudinal (a and c) and shear (b and d) velocity error
against the angle of incidence for various values of mesh density plotted
in linear (a and b) and polar (c and d) plots ........................................ 123
12
List of figures
Figure 4.27 a) Shape of the different rectangular elements used in the mesh; Variation
of the longitudinal (b and d) and shear (c and e) velocity error against the
angle of incidence for various values of R plotted in a polar (b and c) and
linear (d and e) fashion ....................................................................... 124
Figure 4.28 a) Shape of the different rhombic elements used in the mesh; Variation of
the longitudinal (b and d) and shear (c and e) velocity error against the
angle of incidence for various values of g plotted in a polar (b and c) and
linear (d and e) fashion ....................................................................... 125
Figure 4.29 a) Shape of the different parallelogramatic elements used in the mesh;
Variation of the longitudinal (b and d) and shear (c and e) velocity error
against the angle of incidence for various values of g plotted in a polar (b
and c) and linear (d and e) fashion...................................................... 126
Figure 4.30 Variation of the longitudinal (a and c) and shear (b and d) velocity error
against the angle of incidence for various values of mesh density plotted
in polar (a and b) and linear (c and d) plots ........................................ 127
Figure 4.31 a) Shape of the different rectangular elements used in the mesh; Variation
of the longitudinal (b and d) and shear (c and e) velocity error against the
angle of incidence for various values of R plotted in a polar (b and c) and
linear (d and e) fashion ....................................................................... 129
Figure 4.32 a) Shape of the different rhombic elements used in the mesh; Variation of
the longitudinal (b and d) and shear (c and e) velocity error against the
angle of incidence for various values of g plotted in a polar (b and c) and
linear (d and e) fashion ....................................................................... 130
Figure 4.33 a) Shape of the different parallelogramatic elements used in the mesh;
Variation of the longitudinal (b and d) and shear (c and e) velocity error
against the angle of incidence for various values of g plotted in a polar (b
and c) and linear (d and e) fashion...................................................... 131
Figure 4.34 Variation of the longitudinal (a and c) and shear (b and d) velocity error
against the angle of incidence for various mesh densities plotted in a lin-
ear (a and b) and polar (c and d) fashion ............................................ 132
Figure 4.35 a) Shape of the different quadratic isosceles-triangular elements used in
the mesh; Variation of the longitudinal (b and d) and shear (c and e) ve-
locity error against the angle of incidence for various value of f plotted in
a polar (b and c) and linear (d and e) fashion ..................................... 133
Figure 4.36 a) Shape of the different scalene-triangular elements used in the mesh;
Variation of the longitudinal (b and d) and shear (c and e) velocity error
against the angle of incidence for various values of g plotted in a polar (b
and c) and linear (d and e) fashion...................................................... 134
Figure 4.37 Variation of the longitudinal (a and c) and shear (b and d) velocity error
against the angle of incidence for various mesh density plotted in a linear
(a and b) and polar (c and d) fashion .................................................. 135
Figure 4.38 a) Shape of the different quadratic isosceles-triangular elements used in
the mesh; Variation of the longitudinal (b and d) and shear (c and e) ve-
locity error against the angle of incidence for various values of f plotted
in a polar (b and c) and linear (d and e) fashion ................................. 136
Figure 4.39 a) Shape of the different scalene-triangular elements used in the mesh;
Variation of the longitudinal (b and d) and shear (c and e) velocity error
against the angle of incidence for various values of g plotted in a polar (b
and c) and linear (d and e) fashion...................................................... 137
13
List of figures
Figure 4.40 Variation of the longitudinal (a and c) and shear (b and d) velocity error
against the angle of incidence for various mesh densities plotted in a lin-
ear (a and b) and polar (c and d) fashion ............................................ 138
Figure 4.41 a) Shape of the different quadratic isosceles-triangular elements used in
the mesh; Variation of the longitudinal (b and d) and shear (c and e) ve-
locity error against the angle of incidence for various values of f plotted
in a polar (b and c) and linear (d and e) fashion ................................. 139
Figure 4.42 a) Shape of the different scalene-triangular elements used in the mesh;
Variation of the longitudinal (b and d) and shear (c and e) velocity error
against the angle of incidence for various values of g plotted in a polar (b
and c) and linear (d and e) fashion...................................................... 140
Figure 5.1 (a) Longitudinal and (b) shear wave excitation for a square element mesh
and (c) longitudinal and (d) shear excitation for a triangular element mesh
............................................................................................................ 147
Figure 5.2 Straight edge model: a) with edge, b) without edge ........................... 149
Figure 5.3 a) square mesh at 0 degrees aligned with the edge, b) square mesh at 45
degrees, c) triangular mesh ................................................................. 149
Figure 5.4 Implicit models for straight edge: Monitored absolute displacement for a
longitudinal wave excitation using CPE4 and CPE4R meshes at 0 de-
grees, CPE4 and CPE4R meshes at 45 degrees and CPE3, CPE6 and
CPE6M triangular elements. Thin red line is reference for N=30 for each
case...................................................................................................... 151
Figure 5.5 Implicit models for a straight edge: Monitored absolute displacement for
a shear wave excitation using CPE4 and CPE4R meshes at 0 degrees,
CPE4 and CPE4R meshes at 45 degrees and CPE3, CPE6 and CPE6M
triangular elements. Thin red line is reference for N=30 for each case
............................................................................................................ 152
Figure 5.6 Explicit models for a straight edge: Monitored absolute displacement for
a longitudinal wave excitation using CPE4R meshes at 0 degrees, CPE4R
meshes at 45 degrees and CPE3 and CPE6M triangular elements. Thin
red line is reference for N=30 for each case ....................................... 153
Figure 5.7 Explicit models for a straight edge: Monitored absolute displacement for
a shear wave excitation using CPE4R meshes at 0 degrees, CPE4R mesh-
es at 45 degrees and CPE3 and CPE6M triangular elements. Thin red line
is reference for N=30 for each case .................................................... 154
Figure 5.8 Straight crack model: a) with crack, b) without crack ........................ 156
Figure 5.9 Definition of unit long cracks with triangular and square meshes. Blue
line shows modelled crack and red line theoretical crack (which is the
same line with triangular element meshes but not with regular square el-
ement meshes)..................................................................................... 157
Figure 5.10 Implicit models for a crack of unit length: Monitored absolute displace-
ment for a longitudinal wave excitation using mesh made of CPE3,
CPE6, CPE6M, CPE4 and CPE4R elements. Thin red line is reference for
N=30 for each case ............................................................................. 158
Figure 5.11 Implicit models for a crack of unit length: Monitored absolute displace-
ment for a shear wave excitation using mesh made of CPE3, CPE6,
CPE6M, CPE4 and CPE4R elements. Thin red line is reference for N=30
for each case........................................................................................ 158
14
List of figures
Figure 5.12 Explicit models for a crack of unit length: Monitored absolute displace-
ment for a shear and longitudinal wave excitation using mesh made of
CPE3, CPE6M and CPE4R elements. Thin red line is reference for N=30
for each case........................................................................................ 159
Figure 5.13 0.25 unit long crack definition with triangular and square meshes. Blue
line shows modelled crack and red line theoretical crack (which is the
same line with triangular element meshes but not with regular square el-
ement meshes)..................................................................................... 161
Figure 5.14 Implicit models for a crack of length 0.25: Monitored absolute displace-
ment for a longitudinal wave excitation using mesh made of CPE3,
CPE6, CPE6M, CPE4 and CPE4R elements. Thin red line is reference for
N=30 for each case ............................................................................. 162
Figure 5.15 Implicit models for a crack of length 0.25: Monitored absolute displace-
ment for a shear wave excitation using mesh made of CPE3, CPE6,
CPE6M, CPE4 and CPE4R elements. Thin red line is reference for N=30
for each case........................................................................................ 162
Figure 5.16 Explicit models for a crack of length 0.25: Monitored absolute displace-
ment for a shear and longitudinal wave excitation using mesh made of
CPE3, CPE6M and CPE4R elements. Thin red line is reference for N=30
for each case........................................................................................ 163
Figure 5.17 4 unit long crack definition with triangular and square meshes. Blue line
shows modelled crack and red line theoretical crack (which is the same
line with triangular element meshes but not with regular square element
meshes) ............................................................................................... 165
Figure 5.18 Implicit models for a crack of length 4: Monitored absolute displacement
for a longitudinal wave excitation using mesh made of CPE3, CPE6,
CPE6M, CPE4 and CPE4R elements. Thin red line is reference for N=30
for each case........................................................................................ 165
Figure 5.19 Implicit models for a crack of length 4: Monitored absolute displacement
for a shear wave excitation using mesh made of CPE3, CPE6, CPE6M,
CPE4 and CPE4R elements. Thin red line is reference for N=30 for each
case...................................................................................................... 166
Figure 5.20 Explicit models for a crack of length 4: Monitored absolute displacement
for a shear and longitudinal wave excitation using mesh made of CPE3,
CPE6M and CPE4R elements. Thin red line is reference for N=30 for
each case ............................................................................................. 166
Figure 5.21 Circular defect model: a) with circular defect, b) without circular defect
............................................................................................................ 168
Figure 5.22 Unit diameter hole definition with triangular and square meshes ...... 169
Figure 5.23 Implicit models for a hole of unit diameter: Monitored absolute displace-
ment for a longitudinal wave excitation using mesh made of CPE3,
CPE6, CPE6M, CPE4 and CPE4R elements. Thin red line is reference for
N=30 for each case ............................................................................. 170
Figure 5.24 Implicit models for a hole of unit diameter: Monitored absolute displace-
ment for a shear wave excitation using mesh made of CPE3, CPE6,
CPE6M, CPE4 and CPE4R elements. Thin red line is reference for N=30
for each case........................................................................................ 170
Figure 5.25 Explicit models for a hole of unit diameter: Monitored absolute displace-
ment for a shear and longitudinal wave excitation using mesh made of
CPE3, CPE6M and CPE4R elements. Thin red line is reference for N=30
for each case........................................................................................ 171
15
List of figures
Figure 5.26 0.25 unit diameter hole definition with triangular and square meshes 173
Figure 5.27 Implicit models for a hole of diameter 0.25: Monitored absolute displace-
ment for a longitudinal wave excitation using mesh made of CPE3,
CPE6, CPE6M, CPE4 and CPE4R elements. Thin red line is reference for
N=30 for each case ............................................................................. 173
Figure 5.28 Implicit models for a hole of diameter 0.25: Monitored absolute displace-
ment for a shear wave excitation using mesh made of CPE3, CPE6,
CPE6M, CPE4 and CPE4R elements. Thin red line is reference for N=30
for each case........................................................................................ 174
Figure 5.29 Explicit models for a hole of diameter 0.25: Monitored absolute displace-
ment for a shear and longitudinal wave excitation using mesh made of
CPE3, CPE6M and CPE4R elements. Thin red line is reference for N=30
for each case........................................................................................ 174
Figure 5.30 4 units diameter hole definition with triangular and square meshes... 175
Figure 5.31 Implicit models for a hole of diameter 4: Monitored absolute displace-
ment for a longitudinal wave excitation using mesh made of CPE3,
CPE6, CPE6M, CPE4 and CPE4R elements. Thin red line is reference for
N=30 for each casE............................................................................. 176
Figure 5.32 Implicit models for a hole of diameter 4: Monitored absolute displace-
ment for a shear wave excitation using mesh made of CPE3, CPE6,
CPE6M, CPE4 and CPE4R elements. Thin red line is reference for N=30
for each case........................................................................................ 176
Figure 5.33 Explicit models for a hole of diameter 4: Monitored absolute displace-
ment for a shear and longitudinal wave excitation using mesh made of
CPE3, CPE6M and CPE4R elements. Thin red line is reference for N=30
for each case........................................................................................ 177
16
List of Tables
Table 6.1 Table of reflection coefficients due to the tie between two meshes in %
........................................................................................................... 197
Table 6.2 Table of reflection coefficients due to the impedance difference between
two meshes in % ................................................................................. 197
17
Chapter 1
Introduction
1.1. Introduction
During the development process of ultrasonic NDT techniques, the knowledge of the
interaction of waves with defects is key to the achievement of robust and efficient
procedures as well as identifying potential weaknesses. Any problems such as poor
sensitivity to certain shapes and orientation needs to be identified as early as possible.
The reflection of ultrasound from cracks and notches of simple geometry and
orientation is already well understood, but there are few results for more complex
cases, for example from the multiple-crack geometry of stress corrosion cracking, or
from realistic profiles of corrosion.
Against this background, experience has shown that, as NDT is developed for
increasingly challenging applications, it has become increasingly important to start
with careful modelling studies in order to be confident that the procedures which are
developed are optimized and robust. Indeed in the particular topic of guided waves it
is universally accepted that this is the only viable approach. Whereas various wave
models, such as ray models for bulk wave problems or DISPERSE [1] for guided wave
problems can predict the properties of waves in continuous uniform structures, a
discrete approach is often needed to model how the waves interact with discontinuities,
including structural features, cracks, corrosion or other forms of defects. Finite
Element (FE) methods have been used to model a wide range of guided waves
problems [2, 3, 4, 5, 6, 7] and have successfully provided important information about
wave interaction with discontinuities [8, 9, 10, 11, 12, 13]. In these studies, defects
were simplified and straight cracks, notches or cylindrical holes were used to represent
real defects such as fatigue cracks or corrosion patches. One reason for this is that
initial work is bound to focus on the simplest cases, but many modellers are ready to
go on to more complex problems. The reason that so little of this is happening is that,
despite rapid growth in computer power, many of the more complex realistic problems
are still beyond the capacity of the models. The more complex problems require much
18
Chapter 1
Introduction
larger models than the simplified ones. In a similar way, modelling of bulk wave
propagation and scattering has also received a huge amount of interest with many
techniques being used but many cases remain out reach.
1.2. Objectives
One aspect is that related to the existence of unwanted reflections from the boundaries
of a FE model. This has been a limiting factor for FE wave modelling. In time domain
solving, it leads to a large increase in the model geometric size (and therefore a large
increase in the number of degrees of freedom to be solved) as it is generally desirable
19
Chapter 1
Introduction
to separate in time the interaction of the waves with defects from unwanted boundary
reflections. This point is illustrated in Figure 1.1 where it can be seen that the modelled
area is much larger than the area of interest in order to fulfil this. Moreover, in
frequency domain solving, removal of unwanted reflections is a requirement in order
to correctly represent wave propagation in the system [16].
b)
Displacement
time
Time signal
of interest
Figure 1.1. a) 2D plane strain model of a plate including a defect, b) Time signal at the monitoring
point
Being able to represent total radiation outside the area of study in numerical wave
propagation models at a reasonable cost would enable a large increase in
computational efficiency. This topic has attracted a vast amount of interest in the last
25 years [17]. Well known techniques include infinite elements [18], boundary integral
methods [19], non reflecting boundary conditions (NRBC) [20] and absorbing layer
techniques [5, 16, 21, 22, 23, 24].
Another key aspect in the field of NDT is the precise knowledge of wave propagation
velocities. The time of flight of a wave packet has been widely used to determine the
position of defects in structures in many ultrasound industrial applications. Recently,
imaging algorithms [25] and super resolution of defects [26, 27] have been
investigated for NDT applications. Research in these techniques is very much active as
their use holds promises of improved applications not only in NDT but also in the
medical field. In the process of moving from the genesis of a numerical model to its
conclusion and the obtention of results enabling a better understanding of a physical
phenomenon, a modeller takes a series of design decisions. The choice of the solving
20
Chapter 1
Introduction
technique, mesh density or time step influences the successful outcome of the exercise
but also the level of accuracy with which the phenomenon is represented. Discreti-
zation of the space and time leads to a discrepancy between theoretical physical values
and their actual values measured in the model. One of the parameters suffering from
this is the wave propagation velocity in elastic media. It would be beneficial to provide
modellers with information enabling the improvement of not only the quality of their
modelling (adequate decision making) but also the quality of the analysis made based
on the numerical results (quantitatively taking into account the numerical deviation).
There is also an interest in understanding the reasons for mesh scattering where a
change in mesh size causes some numerical reflection.
Achieving better representation of the interaction of waves with complex defects has
led to the development of specialist ways of resolving this interaction. A method
generating interest recently is the fictitious domain technique [28, 29, 30].
Investigation of this technique shows that its implementation in a commercially
available FE packages is not possible. Instead, in this thesis, the path of local mesh
refinement in the context of a regular node grid, is pursued as this could potentially
lead to a large increase in model efficiency and improved representation of the defects.
Overall, the knowledge gained from this work should enable modellers to improve the
quality of their modelling not only by increasing the computational efficiency of their
21
Chapter 1
Introduction
In Chapter 3, the use of two different absorbing layer techniques for elastodynamic
problems is investigated. Implementation of perfectly matched layers [24] (PML) in
the frequency domain in COMSOL Multiphysics and absorbing layer using damping
[31, 32] (ALID) for time and frequency domain in ABAQUS and COMSOL
Multiphysics is devised. Analytical models are developed to enable a quick and
accurate determination of adequate layer parameters. Typical example FE models of
NDT applications are used to demonstrate the added value of using absorbing layers.
The modelling options available to modellers are discussed.
The work presented in Chapter 5 aims to provide the foundations to more complex
defect modelling. The cases investigated are the reflection of bulk waves from straight
edges, straight cracks at an angle and circular holes. For each case, the use of a regular
square mesh approach where the defects are defined by removing or disconnecting
elements is compared with the use of automatically generated meshes made of
22
Chapter 1
Introduction
The work is Chapter 6 focuses on local refinement of a regular square mesh using
conventional tools available in the commercially available FE package ABAQUS. The
method investigated in this part consists of creating 2 separate areas having different
mesh densities and tying them together. The effect of the change in element size in a
regular mesh of elements is investigated. The knowledge of the propagation velocities
gained in Chapter 4 is used. This enables the influence of the way meshes of different
size are tied together on reflections at the interface between the 2 meshes to be
determined. A range of 1D and 2D FE models are created and studied in order to better
understand and quantify the influence of these sources of error. In particular, elastic
properties in parts of the models are adjusted in order to minimize the level of
reflection. Parametric studies are performed in order to provide modellers with a
database of information enabling them to determine the level of reflection which exist
in a range of situations.
23
Chapter 2
Theoretical background
2.1. Introduction
This chapter presents the theoretical background to the work presented in this thesis.
The theory of wave propagation in elastic media is developed and is followed by a
presentation of the FE techniques used in this work.
The theory of wave propagation in elastic media is presented in this section. Bulk
waves and guided waves in a 2D plate are considered. The aim of this section is not
only to present the cases which are of interest to us but also to define the notations used
in following chapters.
The propagation of acoustic waves in unbounded isotropic media has been the subject
of many studies and its theory is widely available [33, 34]. The general equation of
motion relates the particle displacement u in a material of density ρ to the stress field
tensor σ by:
2
∂ ui ∂σ ij
ρ ---------
- = (2·1)
∂t
2 ∂ xj
It is well known that Hooke’s law links the stress σ and strain ε in an isotropic material:
σ ij = λδ ij ε kk + 2με ij (2·2)
1 ∂u i ∂u j
with ε ij = --- ⎛⎝ ------- + -------⎞⎠ (2·3)
2 ∂x j ∂x i
24
Chapter 2
Theoretical background
The Navier governing equation is derived from equations 2·1, 2·2 and 2·3:
2 2
∂ ui ∂u ∂ u
- = ( λ + μ ) ∂ -------j + μ ----------i
ρ --------- (2·4)
∂t
2 ∂ x i ∂x j ∂x j
2
2
ρu·· = ( λ + μ )∇ ( ∇ ⋅ u ) + μ∇ u (2·5)
u = ∇φ + ∇ × ψ (2·6)
Using equation 2·6 to substitute the displacement in the Navier equation 2·4, the
following relation is obtained:
2 ⎛ ∂ 2 φ⎞ 2 ∂ ψ
2
∇ ( λ + 2μ )∇ φ – ρ ⎜ --------2 ⎟ + ∇ × μ∇ ψ – ρ ---------
2
= 0 (2·7)
⎝ ∂t ⎠ ∂t
This equation is fulfilled if both terms are zero which can be interpreted as the
Helmholtz differential equations:
2 2
∂--------
φ ( λ + 2μ )∇ φ 2 2
= ------------------------------- = c L ∇ φ (2·8)
∂t
2 ρ
2 2
μ∂ ψ
2 2∂ ψ
and ∇ ψ = --- --------- = c S --------- (2·9)
ρ ∂t 2
∂t
2
cL and cS are the velocities of longitudinal (dilatational) and shear (rotational) waves
in the infinite isotropic medium and are equal to:
(--------------------
λ + 2μ )- and c = μ
cL = S --- (2·10)
ρ ρ
25
Chapter 2
Theoretical background
ω ω
k L = ----- and k S = ----- (2·11)
cL cS
cL 2πc L 2π cS 2πc S 2π
λ L = ----- = ------------ = ------ and λ S = ----- = ------------ = ------ (2·12)
f ω kL f ω kS
The solution to the Helmholtz differential equations (2·8) and (2·9) for longitudinal or
shear harmonic waves propagating in a given direction is:
i ( ± k Lx x ± k Ly y – ωt ) i ( ± k Sx x ± kSy y – ωt )
φ = Φe and ψ = Ψe (2·13)
with Φ and Ψ the amplitude of the wave and kLx, kLy, kSx and kSy, the projection of the
wavenumbers of the longitudinal and shear waves on the x and y axes.
Let us consider a semi-infinite medium as shown in Figure 2.1. The arrows represent
the longitudinal and shear waves that can exist in the material. In this work, positive
waves are defined as waves propagating in the positive x direction and negative waves
are defined as waves propagating in the negative x direction.
neg
ativ ave y
e lo lw
ngi ud ina
tud git x
n
nega
ina
lw e lo ave
tive ave itiv ear w
shea pos e sh
r wav
e po sitiv
solid
free edge
vacuum
26
Chapter 2
Theoretical background
i ( k Ly y + k Lx x – ωt )
φp = Φp e (2·14)
i ( k Sy y + k Sx x – ωt )
ψp = Ψp e (2·15)
i ( – k Ly y + k Lx x – ωt )
φn = Φn e (2·16)
i ( – kSy y + k Sx x – ωt )
ψn = Ψn e (2·17)
It is also considered that these waves are linked by the Descartes-Snell’s law. This
implies that kLx=kSx=kx. The previous expressions can therefore be simplified, omitting
i ( k x x – ωt )
the common factor e .
ik Ly y
φp = Φp e (2·18)
ik Sy y
ψp = Ψp e (2·19)
– i k Ly y
φn = Φn e (2·20)
– i k Sy y
ψn = Ψn e (2·21)
Using equations 2·2, 2·3 and 2·6, the displacements and stresses are derived:
ik Ly y – i k Ly y ik Sy y – i k Sy y
u x = ik Lx Φ p e + ik Lx Φ n e + ik Sy Ψ p e – i k Sy Ψ n e (2·22)
ik Ly y – i k Ly y ik Sy y – i k Sy y
u y = ik Ly Φ p e – i k Ly Φ n e – i k Sx Ψ p e – i k Sx Ψ n e (2·23)
2 2 ik Ly y 2 2 – i k Ly y
σ yy = – ( ( λ + 2μ )k Ly + λk Lx )Φ p e – ( ( λ + 2μ )k Ly + λk Lx )Φ n e
ikSy y – i k Sy y
(2·24)
+ 2μk Sx k Sy Ψ p e – 2 μk Sx k Sy Ψ n e
ik Ly y – i k Ly y
σ xy = – 2 μk Lx k Ly Φ p e + 2μk Lx k Ly Φ m e
2 2 ik Sy y 2 2 – i k Sy y (2·25)
– 2μ ( k Sy – k Sx ) Ψ p e – 2μ ( k Sy – k Sx ) Ψ n e
27
Chapter 2
Theoretical background
These expressions will be needed in Chapter 3 to model a multi-layered system and can
be linked to the potentials in the following way:
ux Φp
uy Ψp
= M (2·26)
σ yy Φn
σ xy Ψn
with M=
ik Ly y ik Sy y – i k Ly y – i k Sy y
ik Lx e ik Sy e ik Lx e – i k Sy e
ik Ly y ik Sy y – i k Ly y – i k Sy y
ik Ly e – i k Sx e – i k Ly e – i k Sx e
2 2 ik Ly y ik Sy y 2 2 – i k Ly y – i k Sy y
– ( λk + 2μk Lx )e 2μk Sx k Sy e – ( λk + 2μk Lx )e – 2 μk Sx k Sy e
ik Ly y 2 2 ikSy y – i k Ly y 2 2 – i k Sy y
– 2 μk Lx k Ly e – 2μ ( k Sy – k Sx ) e 2μk Lx k Ly e – 2μ ( k Sy – k Sx ) e
(2·27)
This expression enables to link all displacements, stresses and wave potentials and can
be used to describe a multi-layered system straightforwardly.
y λx
ψn
φn
h
x
h φp
ψp
28
Chapter 2
Theoretical background
Guided waves in plates propagate in the same way as bulk waves in the infinite
medium except the fact that the internal longitudinal and shear waves are reflected at
the free surfaces. The theory of guided waves is treated in detail in several textbooks
[33, 35]. Guided waves result from this interaction of both wave types with the top and
bottom surfaces of the plate where Snell’s law is satisfied. This implies that both wave
types have the same wave number kx in the x direction. As represented in Figure 2.2,
in general there can be up to two shear waves and two longitudinal waves. Following
the notation used previously for bulk waves, the field F can be expressed as:
i ( k x x – ωt )
F = φ + ψ = ( φ y ( y ) + ψ y ( y ) )e (2·28)
ik Ly y – i k Ly y
where φy ( y ) = Φp e + Φm e (2·29)
ik Sy y – i k Sy y
ψy ( y ) = Ψp e –Ψm e (2·30)
It is essential to note that, at this stage, Φp, Φn, Ψp, Ψn and kx are unknown.
The symmetry of the system imposes that the amplitude of waves going up is either
equal to or opposite to the amplitude of waves going down:
This leads to the existence of two types of guided modes: symmetric and anti
symmetric modes as represented in Figure 2.3.
a) b)
Figure 2.3. Typical deformation caused by symmetric (a) and anti-symmetric (b) modes.
29
Chapter 2
Theoretical background
i ( k x x – ωt )
When the amplitudes are equal, symmetric modes exist. The common factor e
is omitted and the following expressions are obtained:
φ = Φ s cos ( k Ly y ) (2·32)
ψ = Ψ s sin ( k Sy y ) (2·33)
Equations 2·2, 2·3 and 2·6 are still valid in this case. The displacements and stresses
are derived:
2 2
σ xx = – ( k x ( λ + 2μ ) + k Ly λ ) Φ s cos ( k Ly y ) – 2μik x k Sy Ψ s cos ( k Sy y ) (2·36)
2 2
σ xy = – 2 μik x k Ly Φ s sin ( k Ly y ) + μ ( k Sy – k x )Ψ s sin ( k Sy y ) (2·38)
2 2
σ yy = – ( ( λ + 2μ )k Ly + λk x ) Φ s cos ( k Ly y ) + 2μik x k Sy Ψ s cos ( k Sy y ) (2·40)
When the amplitudes are opposite, anti symmetric modes exist. The common factor
i ( k x x – ωt )
e is omitted and the following expressions are obtained:
φ = Φ a sin ( k Ly y ) (2·42)
ψ = Ψ a cos ( k Sy y ) (2·43)
30
Chapter 2
Theoretical background
2 2
σ xx = – ( k x ( λ + 2μ ) + k Ly λ ) Φ a ( sin ( k Ly y ) + 2μik x k Sy Ψ a sin ( k Sy y ) ) (2·46)
2 2
σ xy = 2μik x k Ly Φ a cos ( k Ly y ) + μ ( k Sy – k x )Ψ a cos ( k Sy y ) (2·48)
2 2
σ yy = – ( ( λ + 2μ )k Ly + λk x ) Φ s sin ( k Ly y ) – 2μik x k Sy Ψ s sin ( k Sy y ) (2·50)
Φa, Φs, Ψa, Ψs and kx are unknown. In order to calculate them, the traction-free
boundary conditions on the top and bottom surfaces (where y = ± h ) are used:
σ yy = 0 and σ xy = 0 (2·52)
I s cos ( k Ly h ) G s sin ( k Ly h )
Ψ s = – ----------------------------- and Ψ s = – ------------------------------ (2·54)
J s cos ( k Sy h ) H s sin ( k Sy h )
This gives the dispersion equation that enables the calculation of kx:
I s cos ( k Ly h ) G s sin ( k Ly h )
----------------------------
- = -----------------------------
- (2·55)
J s cos ( k Sy h ) H s sin ( k Sy h )
2 2
– ( ( λ + 2μ )k Ly + λk x ) cos ( k Ly h ) – 2μi k x k Ly sin ( k Ly h )
----------------------------------------------------------------------------- = --------------------------------------------------
- (2·56)
2μik x k Sy cos ( k Sy h ) 2 2
μ ( k Sy – k x ) sin ( k Sy h )
31
Chapter 2
Theoretical background
2
tan ( k Sy h ) 4μk x k Ly k Sy
or ------------------------ = ---------------------------------------------------------------------------
- (2·57)
tan ( k Ly h ) 2
– ( ( λ + 2μ )k Ly + λk x ) ( k Sy – k x )
2 2 2
2
tan ( k Sy h ) 4k x k Ly k Sy
which can be further simplified to: ------------------------ = – -------------------------2 (2·58)
tan ( k Ly h ) 2 2
( k Sy – k x )
For anti symmetric modes, in a similar way, the dispersion equation is:
2 2 2
tan ( k Sy h ) ( k Sy – k x )
-----------------------
- = – ------------------------- (2·59)
tan ( k Ly h ) 2
4k x k Ly k Sy
Equation 2·58 and 2·59 were first published by Lamb [36]. These are solved
numerically as there are no analytical solutions to these equations. The solutions give
values of kx which can be either real (propagating mode), imaginary (evanescent mode)
or complex (propagating evanescent mode). Figure 2.4 illustrates the deformation of a
plate caused by these various types of modes.
a) b) c)
y
x
Propagating waves are those which have received the largest amount of research in the
field of NDT as they can be used for long range inspection of structures such as rails,
plates or pipes. For a given set of properties (plate thickness, bulk wave numbers and
frequency), there is a finite number of real and imaginary wave numbers and an infinite
number of complex wave numbers. It implies the existence of a finite number of
propagating and evanescent waves and an infinite number of propagating evanescent
modes. DISPERSE [1, 37], a software tool developed at Imperial College, provides a
useful tool to find these solutions. Figure 2.5 shows an example of dispersion curves
of propagating modes in a typical engineering structure (3mm thick steel plate).
32
Chapter 2
Theoretical background
0.6
0.5
0.4
k (1/mm)
0.3
0.2
0.1
0
0 0.5 1.0 1.5 2.0 2.5
Frequency (MHz)
Figure 2.5. Wave number against frequency for a 3mm thick steel plate
These can be presented in many forms but a common form is to represent the phase
velocity of the modes against the frequency-thickness product, as shown in Figure 2.6.
In this figure, it can be seen that the phase velocities of guided waves in a plate vary
with the frequency thickness product. Guided waves in a plate are therefore described
as dispersive.
12
10
8
Vph (m/ms)
0
0 1.5 3.0 4.5 6.0 7.5
Frequency.thickness (MHz.mm)
Figure 2.6. Phase velocity against frequency.thickness for a 3mm thick steel plate
In practical NDT, wave packets (often tone bursts) are usually used. Since guided
waves are dispersive, the wave packet will not retain its original shape as it travels. In
this case, it is convenient to define the velocity of the packet using the group velocity,
found from the wave number curves.
33
Chapter 2
Theoretical background
The two previous figures illustrate the existence of several guided modes whose
existence depends on the frequency and whose velocity varies against the frequency.
In order to identify the modes, each mode is named by a letter and a number.
Symmetric modes use ‘S’ and anti-symmetric modes ‘A’. The number is a counter
variable starting at 0 used to distinguish modes in their family. S0 is the fundamental
symmetric mode, S1, the first symmetric mode to appear as the frequency increases,
and so on.
Another aspect of a guided mode is its mode shape. As can be seen from the
displacement and stress equations derived previously in this section, the amplitudes of
these values vary through the thickness of the plate. Knowledge of the mode shapes is
essential to NDT and these are commonly plotted as shown in Figure 2.7.
Bottom
of plate
DISPERSE In plane displacement DISPERSE out of plane displacement
Figure 2.7. Example of A0 mode shapes for a free plate case at different frequencies, shown for a 3mm
thick steel plate
In this work, wave propagation problems are solved using both implicit and explicit
solvers. Time domain explicit models are run with ABAQUS/Explicit [14]. For
implicit models, ABAQUS/Standard [14] and COMSOL Multiphysics [15] are used.
It is important to emphasize that ABAQUS/Explicit and ABAQUS/Standard are two
separate programs which are part of the ABAQUS package. This section presents the
details of the numerical implementation in these cases.
34
Chapter 2
Theoretical background
u·· , u· , u are the acceleration, velocity and displacement respectively. [M] is the
diagonal lumped mass matrix whose values are determined by the density of the
material used. [K] is the static stiffness matrix whose values are defined by the Young’s
modulus and Poisson’s ratio. [C] is the viscous damping matrix which is determined
by the Rayleigh damping. Stiffness or mass proportional damping can be introduced:
[ C ] = CM [ M ] + CK [ K ] (2·61)
where CM and CK are the mass and stiffness proportional damping coefficient. [Fa] is
the external force.
Wave propagation occurs when the initial equilibrium is disturbed by the application
of forces or displacement constraints on nodes. Commonly, these are applied in the
form of a tone burst. The central difference operator links the displacement, velocity
and acceleration in the following way:
(i + 1) (i)
( i + 0.5 ) ( i – 0.5 ) Δt + Δt ( i )
u· = u· + ----------------------------------- u·· (2·62)
2
(i + 1) (i) ( i + 1 ) · ( i + 0.5 )
u = u + Δt u (2·63)
with Δt, the time increment and i, the time increment number.
The procedure is explicit because the process is advanced using known values from the
previous time step.
35
Chapter 2
Theoretical background
The computational efficiency comes from two aspects of the explicit scheme. Firstly,
unlike in the implicit scheme, there is no need to assemble and invert the global mass
matrix. Secondly, the diagonal mass matrix is used when calculating the acceleration
in the following way:
where [F] is the applied load vector and [I] is the internal force vector. [M]-1 is easily
calculated because [M] is diagonal.
This implementation is conditionally stable and the time step Δt has to be smaller than
the critical time step Δtcr which in an undamped system depends on the highest
frequency in the smallest element [14]:
2
Δt ≤ Δ t cr = ------------ (2·65)
ω max
ΔL
Δt ≤ Δt cr = ------- (2·66)
cL
where ΔL is the smallest element size and cL is the velocity of the dilatational wave.
Figure 2.8 illustrates the value of ΔL for a linear square element, a linear equilateral-
triangular element and a quadratic equilateral-triangular element. It can be noted that
this distance corresponds to the shortest distance between 2 nodes for a square element
but this is not the case for triangular elements. This will be verified in Section 4.2.3.1.
36
Chapter 2
Theoretical background
a) b) c)
ΔL
ΔL ΔL
Figure 2.8. Illustration of ΔL for a a) linear square element, b) linear triangle element and c) quadratic
triangle element.
The condition in equation 2·66 is referred to as the CFL condition after R. Courant, K.
Friedrichs, and H. Lewy [39]. The actual time step used in a model can be expressed
in terms of its ratio by the critical time step. This ratio is called the Courant (or CFL)
number:
Δt
CFL = --------- (2·67)
Δt cr
ΔL
Δt = CFLΔt cr = CFL ------- (2·68)
cL
Given that the space usually needs to be discretized in such a way that the shortest
wavelength in the model is finely discretized (usually, at least 7 nodes per shortest
wavelength are used), a large number of increments is required to solve a model. The
number of increment in an explicit scheme is larger than the one required with a time
domain implicit scheme but, since global mass and stiffness matrices do not need to be
formed nor inverted, each increment is computationally inexpensive. This makes the
explicit scheme particularly attractive for wave propagation.
Although it is possible to use implicit schemes to solve time domain wave propagation
problems, the central difference explicit scheme has proved to be far more efficient
than the implicit one. Therefore, in this work, time domain implicit schemes are not
investigated and the implicit method is only used to solve wave problems in the
37
Chapter 2
Theoretical background
frequency domain (i.e. subjected to continuous harmonic excitation). The steady state
response of the model forced by an harmonic excitation (i.e. at a single frequency) is
calculated. In this case, it is necessary for waves to radiate out of the area of study in
order to model wave propagation. The particular aspect of radiation outside the area of
study is the topic of Chapter 3.
This procedure provides frequency domain results which can be used to reconstruct
time domain results. The reconstruction procedure is presented in Section 3.5.2.
where σ is the stress, t is the surface traction and δε is the strain variation compatible
with the displacement variation δu . The : operation means that corresponding
conjugate components of the stress and strain rate matrices are multiplied as pairs and
the products summed.
N NM ·· M NM · M N N
δu { M u + CM M u + I –P } = 0 (2·70)
MN N M
with the mass matrix M = ∫ ρN ⋅ N dV (2·71)
V
N N
the internal load vector I = ∫β : σ dV (2·72)
V
N N
and the external load vector P = ∫N ⋅ t dS (2·73)
St
38
Chapter 2
Theoretical background
N
NN is the shape function of the specified elements and β defines the strain variation
from the variation of the kinematic variables. The change in internal force vector is:
N N N
ΔI = ∫ [ Δβ : σ + β : Δσ ] dV (2·74)
V
el ·
Δσ = D : ( Δε + C K Δε ) (2·75)
el
where D is the elasticity tensor. The strain and strain rate become:
M M · M M
Δε = β Δu and Δε = β Δu· (2·76)
NM N el M
K = ∫β : D : β dV (2·77)
V
N NM ·· M NM NM M N N
δu { M u + ( CM M + CK K )u· + K u – P } = 0 (2·78)
M M M
Δu = ( ℜ ( u ) + iℑ ( u ) ) exp ( iωt ) (2·79)
M M M
Δu· = ( – ω ℑ ( u ) + iωℜ ( u ) ) exp ( iωt ) (2·80)
M 2 M M
Δu·· = – ω ( ℜ ( u ) + ℑ ( u ) ) exp ( iωt ) (2·81)
M M M
ΔP = ( ℜ ( P ) + iℑ ( P ) ) exp ( iωt ) (2·82)
39
Chapter 2
Theoretical background
NM 2 NM NM NM
K –ω M –ω ( CM M ) ℜ ( uM )
+ CK K M
= ℜ(P ) (2·83)
NM NM NM 2 NM M M
–ω ( CM M + CK K ) –( K – ω M ) ℑ(u ) ( – ℑ )P
Solving this equation provides the complex response of the system. Provided that
radiation out of the area of study occurs, this response describes the wave propagation
phenomenon occurring in the model.
∇ ⋅ ( – c ∇u ) + au = 0 in Ω (2·84)
n ⋅ ( c∇u ) = g in ∂Ω (2·85)
hu = r in ∂Ω (2·86)
Ω is the computational domain, ∂Ω is the domain boundary and n is the outward unit
normal vector on ∂Ω . u is a dependent variable unknown on the computational domain
and can be a scalar or a vector, c, a, g, h and r are constants or matrices depending on
the form of u.
u = u1 x + u2 y (2·87)
40
Chapter 2
Theoretical background
a 11 a 12
a = (2·89)
a 21 a 22
∂u 1 ∂u 1 ∂u 2 ∂u 2
- ∂ ⎛ c 1111 + c 1112 + c 1211 + c 1212 ⎞ (2·90)
∂x⎝ ∂x ∂y ∂x ∂y ⎠
∂u 1 ∂u 1 ∂u 2 ∂u 2
- ∂ ⎛ c 1121 + c 1122 + c 1221 + c 1222 ⎞ + a 11 u 1 + a 12 u 2 = 0 (2·91)
∂y⎝ ∂x ∂y ∂x ∂y ⎠
∂u 1 ∂u 1 ∂u 2 ∂u 2
- ∂ ⎛ c 2111 + c 2112 + c 2211 + c 2212 ⎞ (2·92)
∂x⎝ ∂x ∂y ∂x ∂y ⎠
∂u 1 ∂u 1 ∂u 2 ∂u 2
- ∂ ⎛ c 2121 + c 2122 + c 2221 + c 2222 ⎞ + a 21 u 1 + a 22 u 2 = 0 (2·93)
∂ y⎝ ∂x ∂y ∂x ∂y ⎠
Since the equations of dynamic equilibrium in the frequency domain for a 2D isotropic
elastic medium are:
∂ ⎛ ( λ + 2μ ) ∂u 1 + λ ∂u 2⎞ + ∂ ⎛ μ ∂u 1 + μ ∂u 2⎞ + ρω 2 u = 0 (2·94)
∂x⎝ ∂x ∂y ⎠ ∂y⎝ ∂y ∂x ⎠ 1
∂ ⎛ ∂u 1 ∂u 2⎞ ∂ ⎛ ∂u 1 ∂u 2⎞ 2
μ + μ + λ + ( λ + 2μ ) + ρω u 2 = 0 (2·95)
⎝
∂x ∂y ∂x ⎠ ⎝
∂y ∂x ∂y ⎠
41
Chapter 2
Theoretical background
It is therefore possible to solve Equations 2·94 and 2·95 using Equation 2·84 by
defining c and a as
λ + 2μ 0 00000 000λ000
0 00μ00 00μ00 000
c = (2·96)
00000 00μ00 00μ00 0
000λ000 000 0 λ + 2μ
2
a = ρω 0 (2·97)
2
0 ρω
The model is excited by either applying a force on a boundary with Equation 2·85 or
prescribing a displacement using Equation 2·86. The model is solved using one of
COMSOL built-in solvers and solutions can be output to be post processed in
MATLAB. As with ABAQUS/Standard, provided that radiation out of the area of
study exists, the response obtained represents the wave propagation phenomenon
occurring in the system.
2.4. Conclusions
In this chapter, the theoretical background necessary to the work presented in this
thesis was presented. The theory of wave propagation in elastic media was developed
and was followed by a presentation of the FE techniques used in this work.
42
Chapter 3
Modelling waves in unbounded elastic
media using absorbing layers
3.1. Introduction
Unwanted reflections from the boundaries of the system have been a limiting factor for
FE modelling of waves. In time domain solving, this leads to a large increase in the
model geometric size (and therefore a large increase in the number of degrees of
freedom to be solved) as it is generally desirable to separate in time the interaction of
the waves with defects from unwanted boundary reflections. This issue is illustrated in
Figure 3.1 where it can be seen that the plate had to be extended to achieve separation
of the signal of interest and the unwanted reflections. Moreover, in frequency domain
solving, removal of unwanted reflections is a requirement in order to correctly
represent wave propagation in the system.
b)
Displacement
time
Time signal
of interest
Figure 3.1. a) 2D plane strain model of a plate including a defect, b) Time signal at the monitoring
point
43
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
(NRBC) [20] and absorbing layer techniques [5, 21, 22, 23, 24]. These various
techniques have proved successful in their context but, as engineers involved in
developing ultrasound NDT applications, our goal is to find computationally efficient
techniques for removing unwanted reflections that can be used with commercial FE
packages. Commercial packages were chosen from the outset of the study not only
because they offer memory efficient and robust solvers but also because they remove
the necessity of developing and maintaining specialist codes. It is key that the
technique that achieves the removal of unwanted boundary reflections is easily
implementable in the commercial FE packages. It needs to be able to remove the
reflections with a high degree of accuracy (typically, in this work, 60dB ≡ 99.99% of
the incident wave amplitude removed) as ultrasound NDT relies on quantitatively
evaluating relatively low amplitude signals. A significant reduction in model
geometric size compared with the classical technique (i.e. increasing model geometric
size) needs to be achieved in order to justify the use of the technique.
In this work, I have focused on two different absorbing layer techniques for
elastodynamic problems. I implemented “Perfectly Matched Layers” [24] (PML) in
the frequency domain in COMSOL [15] and “Absorbing Layer Using Damping” [31,
32] (ALID) for time and frequency domain in ABAQUS [14] and COMSOL.
Optimum definition of the layer parameters is essential to improve modelling
capabilities, but it becomes counter efficient to invest a large amount of time to
determine the best parameters to use for these layers. This issue is raised in several
papers [40,41]. In order to resolve this, we have developed analytical models to enable
a quick and accurate determination of adequate layer parameters.
One of the challenges of this chapter is to stimulate interest from fellow researchers,
engineers and modellers. To convince them of the advantage of the techniques, two
typical examples demonstrate the added value of using absorbing layers.
44
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
Infinite elements are a special type of element with modified properties which can be
used in conjunction with standard finite elements and which simulate an infinite space.
A general presentation of the development of infinite elements can be found in [18]. A
single row of infinite elements is positioned outside the boundaries of the area of study
as shown in Figure 3.2.
Infinite elements have been proved to work relatively well for static cases as well as in
certain domains of wave propagation: electromagnetism, acoustics and elastic bulk
waves with incidence at angles close to the normal. It should be noted that elastic wave
problems are more complex than acoustic or electromagnetic ones because in 2D two
wave types exist in this case. In the general case of wave propagation, the simulation
of an infinite expanse of material can be defined using the Sommerfeld radiation
condition [42] which defines the condition for total radiation of a wave from a source.
Satisfying exactly the radiation condition in the infinite element results in perfect
absorption and no reflection from the boundary. For solid media, ABAQUS provides
infinite elements based on [43]. A benchmark problem for the use of infinite elements
is available in the ABAQUS user manual [14]. This example is similar to the one
analysed in [44]. The geometry of the model is presented in Figure 3.3.a. The problem
45
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
is an infinite half-space subjected to a vertical pulse line load in the form of a 10MHz
raised-cosine function 1-sin(ωt/3)*cos(ωt) with a spatially constant amplitude of
1GPa. The vertical displacement is monitored at point A 2mm below the edge of the
load as indicated on Figure 3.3.a. The resulting displacement is plotted on Figure 3.3.b.
A comparative model is created. This model has the same properties as the benchmark
but does not use infinite elements. Instead, the boundaries of the area of study where
the infinite elements were placed are moved away from point A so that the monitored
signal at point A does not include any interaction of the waves with these boundaries.
This represents an exact reference and the difference between the two models is the
reflection from the infinite element boundary. This model is shown in Figure 3.3.c. and
the displacement monitored at point A is plotted in Figure 3.3.d.
U2 (m)
A Shear incident
0
Area of study -1
0 1 2 3 4
time (sec) x 10 -6
c) Excitation d) x 10 -6
1
Constraint
U2 (m)
A
0
-1
0 1 2 3 4
time (sec) x 10 -6
Figure 3.3. ABAQUS benchmark model: a) Model geometry, b) vertical displacement at point A,
Extended model (reference): c) Model geometry, d) vertical displacement at point A
There is a clear discrepancy between the 2 monitored signals. This indicates that a
noticeable reflection occurs as the wave reaches the infinite elements. This example
confirms that infinite elements are not suitable for high accuracy removal of unwanted
reflection of bulk waves. Studies [45, 46] have proved that similar conclusions can be
drawn for the use of these elements with guided waves. These findings justify why this
technique is not investigated further in this chapter.
46
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
Non reflecting boundary conditions (NRBC) are special boundary conditions used in
FE or finite difference (FD) methods to model wave propagation in unbounded media.
The dimensions of the model are the same as the area of study; only the boundary
conditions are changed. The modified boundary conditions generally use extra
variables in order to approximate the infinite expanse of the medium. Many different
methods have been developed for a wide range of fields. Early work in the ‘70s on
NRBC was dominated by the use of the radiation condition given by Sommerfeld [42]
but this method had limits as demonstrated in [47]. Other methods were developed
such as the Engquist-Majda [48], Bayliss-Turkel [49], and are still widely used today.
Exact non-local methods based on the Dirichlet-to-Neumann map [50, 51] or the
difference potential method were developed from the late ‘80s. In the last 10 years,
high order local NRBC have been developed for various cases [52]. These techniques
were proved to work correctly but, as they require modification of the standard solving
procedure, the development of specialist FE codes is required. To the best of the
author’s knowledge, they cannot be implemented in commercially available FE
packages and are therefore not investigated here.
3.3.1 Concept
Absorbing layers are finite regions “attached” at the extremities of a model - see Figure
3.4. Their objective is to approximate the case of an unbounded problem by absorbing
waves entering them. Small reflections from the absorbing region exist but these can
be made acceptably small by correctly defining the layer’s parameters.
a) b) c)
Figure 3.4. Absorbing layer concept for 2D models: a) infinite medium, b) semi infinite medium, c)
plate
47
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
The PML technique was created in 1994 by Berenger for electromagnetism [24] and
has been extended to other fields such as acoustics [53, 54, 55] and seismological and
other elastic waves [56, 57, 58, 59]. As its name indicates, a PML matches perfectly
the impedance of the area of study. This means that, in theory, a wave enters a PML
without reflection. Once inside it, the wave decays exponentially. A PML can therefore
be used to achieve total radiation of a wave out of the area of study.
– iωt ikx
u ( x, y, t ) = u 0 e e (3·1)
For an elastic isotropic medium in 2D plane strain, the equations of equilibrium are:
∂⎛ ∂u x ∂u y⎞ ∂ ⎛ ∂u x ∂u y⎞ 2
( λ + 2μ ) + λ + μ + μ + ρω u x = 0 (3·2)
∂x ⎝ ∂x ∂y ⎠ ⎝
∂y ∂y ∂x ⎠
∂ ⎛ μ ∂u x ∂u y ∂u x ∂u y
+ μ ⎞ + ∂ ⎛λ + ( λ + 2μ ) ⎞ + ρω u y = 0
2
(3·3)
∂x⎝ ∂y ∂ x ⎠ ∂ y⎝ ∂ x ∂y ⎠
In the PML, we want the wave to decay exponentially and therefore to have the
following form.
– iωt ikx – kα x x
u ( x, y, t ) = u 0 e e (3·4)
The modification from equation 3·1 to equation 3·4 can be achieved by performing a
change of variable in the PML.
x → x ( 1 + iα x ) (3·5)
48
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
∂ 1 ∂ ∂
→ ⎛ -----------------⎞ = Sx (3·6)
∂x ⎝ 1 + iα x ∂ x ⎠ ∂x
After this change of variable for x and y and some algebra, the equations of equilibrium
3·2 and 3·3 become:
∂ ⎛ ∂u x S-----x ∂u y⎞ ∂ ⎛ ∂u x S-----y ∂u y⎞ ρω 2
⎜μ + μ ⎟ + ⎜λ + ( λ + 2μ ) ⎟ + ---------- u y = 0 (3·8)
∂ x⎝ ∂ y Sy ∂ x ⎠ ∂ y⎝ ∂ x Sx ∂ y ⎠ Sx Sy
In the area of study, αx and αy are equal to zero, hence Sx and Sy are equal to 1 and
equation 3·7 and 3·8 become equation 3·2 and 3·3. In the PML, αx and αy are
respectively x and y dependent and are defined as illustrated in Figure 3:
p p
α x ( x ) = A x ⋅ x and α y ( y ) = A y ⋅ y (3·9)
αx(x)
x
Layer
αy(y) thickness
Area of study
The same approach can be readily extended to 3 dimensional cases by following the
same reasoning.
49
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
In order to avoid numerical issues due to the discontinuity of the attenuation parameter
and its first and second derivatives at the interface between the area of the study and
the PML, the attenuation parameters are taken to vary at least quadratically with the
distance which implies that p in equation 3·9 has to be superior or equal to 2.
At the end of the PML, the decaying wave is reflected and generally partially mode
converted. The reflected waves propagate back toward the area of study while
continuing to decay and then re-enter the area of study. The PML has to be defined so
that these reflections are negligible.
Correct definition of the PML parameters (thickness of the layer, α and p) is essential
to achieve an efficient model. These points are looked at in detail later in this chapter.
Solving of the model is done in the frequency domain but time domain results can be
obtained easily by solving the model over a range of frequency and performing an
inverse Fourier transform.
The ALID is an absorbing layer which is made of a material with the same properties
as those of the area of study apart from having a gradually increasing damping. The
general concept was mentioned in 1980 by Israeli and Orszag in [22] and was recently
revived by Liu et al [21] and Castaings et al [16].
– iωt ikx
By convention, we consider harmonic waves of the form: u ( x, y, t ) = u 0 e e . We
have u· → – i ωu . The equation of equilibrium in the frequency domain is:
2
– [ M ] ω u – [ C ] iωu + [ K ]u = f (3·11)
50
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
[ C ] = CM [ M ] + CK [ K ] (3·12)
where CM and CK are the mass and stiffness proportional damping coefficients.
Equation 3·11 becomes:
2
– [ M ] ω u – C M [ M ] + C K [ K ] iωu + [ K ]u = f (3·13)
CM 2
– [ M ] ⎛ 1 + i -------⎞ ω u + [ K ] ( 1 – iω C K )u = f (3·14)
⎝ ω⎠
In an ALID with a boundary perpendicular to the x axis, the value of CM and CK are
gradually increased in the x direction. We set the following formulation:
p p
C M ( x ) = C Mmax ⋅ X ( x ) and C K ( x ) = C Kmax ⋅ X ( x ) (3·15)
CMmax and CKmax are positive real numbers and X(x) varies from 0 at the interface
between the ALID and the area of study to 1 at the end of the ALID following a power
law whose order is defined by p. In a time domain model, CMmax and CKmax are
constant and therefore do not vary with the frequency. This leads to a variation of the
wave attenuation against the frequency. In a frequency domain model, it can be made
constant over the range of frequencies by adjusting CMmax and CKmax accordingly with
frequency.
The effect of Rayleigh damping can be better understood by using a formulation using
complex density and complex stiffness modulus in the ALID:
CM ( x )
ρ ALID = ρ ⎛ 1 + i ---------------⎞ or E ALID = E ( 1 – iω C K ( x ) ) (3·16)
⎝ ω ⎠
51
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
The wavenumber is proportional to the square root of the density ρ and inversely
proportional to the square root of the Young’s modulus, E:
k ∝ ω --ρ- (3·17)
E
From this, we see that, in the ALID, the wave number is complex and the solution to
wave equation becomes:
k’ and k’’ are the real and imaginary parts of the wave number and are both positive
numbers for a wave propagating in the positive direction. k’’ induces the decay of the
wave in the layer. The real part of the wave number in the ALID varies with the
damping and does not match perfectly the real part of the wave number in the area of
study. This denotes a change in acoustic impedance which will cause reflections at the
interface of the ALID as well as inside it. To simulate unboundedness with high
accuracy, these reflections need to be negligible. Increasing the damping gradually
enables us to achieve this. The technique is not as neat as the PML technique since
there is a change of impedance between the area of study and the ALID, and inside the
ALID itself, but it has the advantage of being easily implementable in both time and
frequency domain solving in most commercially available FE packages.
One important point to note is that the introduction of damping decreases the value of
the stable time increment when solving the model with the central difference explicit
scheme [14]. The value of the damping at the end of the ALID is usually very large
compared to values commonly used in structures. A high value of CM causes a
relatively small decrease in the stable increment whereas one of CK usually has a very
strong effect leading to a great loss in computational efficiency (e.g. time increment
divided by a thousand or more). Therefore, it is preferable to avoid using CK to define
ALID with an explicit scheme. For this work, we have only used CM for time and
frequency domain studies:
p
C M ( x ) = C Mmax ⋅ X ( x ) and C K ( x ) = 0 (3·19)
52
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
It is easy to see that correct definition of the layer parameters (length of the layer L,
variation of the attenuation parameter CM and the power law p) is essential to achieve
an efficient model.
Running FE test cases to determine layer parameters has proved to be highly time
consuming and inefficient. It would be beneficial for modellers to be able to use tools
to define the layer parameters quickly, reliably and efficiently. Analytical models have
therefore been created to evaluate the reflection coefficient of absorbing layers for bulk
and guided waves cases. In the following parts, the explanation of what these models
53
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
represent and why layers parameters can be defined using them are given. These
models guarantee the fulfilment of the acceptability criterion and greatly speed up the
correct definition of layer parameters.
All the following models are defined in the frequency domain, as harmonic waves are
considered, and they can be extended readily to work in the time domain. For a time
domain model, the modeller needs to know the frequency content of its signal and
validate the absorbing layer over the range of interest. The parts of the frequency
spectrum whose amplitude is smaller than the acceptability criterion can be ignored.
The rest of the frequency spectrum defines the frequency range of interest. It is
advantageous to take into account the actual frequency spectrum as it relaxes the
acceptability criterion for the low amplitude parts of the frequency spectrum.
Analytical models presented in this section have been developed to evaluate the
reflection coefficient of harmonic longitudinal and shear plane waves from absorbing
layers over a range of frequencies and angles. Although all models in this part are 2
dimensional, they are also valid for 3 dimensional cases. The intention is that these
models can be used to validate the acceptability of a given absorbing layer for both
longitudinal and shear waves over any desired range of frequencies and angles of
incidence. As one considers the geometry of an area of study, it is easy to see that shear
and longitudinal waves coming from a range of sources may be incident at the
absorbing layers. As will be seen later, knowledge of the range of angles at which these
waves are incident at the absorbing layer is one consideration needed to define the
absorbing layers. The range of angles of incidence is determined by the geometry of
the model and the position of wave sources such as excitation points or scatterers in the
model, as illustrated in Figure 3.6.
Frequency domain technique using monochromatic plane waves rather than transient
waves generated by point sources are used in the model as this greatly simplifies the
approach and provides a conservative value of the reflection from the layer. In the real
model, waves decay (beam spread) before reaching the absorbing layer and reflections
54
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
θ
θ
Area of study
Absorbing layer
Wave source
θ
θ
Figure 3.6. Illustration of extreme angles defining the range of angles to consider when dimensionning
an absorbing layer.
spread spatially (unless the incidence is purely normal) as the reflection phenomenon
will occur at different locations in the layer - see Figure 3.7. In the time domain, the
reflection will also be spread in time as a consequence of the previous point.
Reflected waves
{ { Transmitted waves
Incident wave
Layer number 1 2 3
Figure 3.7. Spatial spread of the reflection and transmission for a single layer (no mode conversion
shown for simplicity)
In order to confirm the validity of the analytical models, their results are compared with
results from FE models. The challenge is to reproduce exactly the plane wave
interaction with the absorbing layers. As, exceptionally in this case, the reflected
amplitude from an absorbing layer is the one of interest, it is necessary to use extremely
55
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
low (-120dB) reflection ALID referred to below as “ultra” ALID in order to correctly
monitor the reflection from the absorbing layer studied.
One model shown in Figure 3.8.a is used for normal incidence. This is a 2D plane strain
model where the displacements on the top and bottom surfaces are constrained to zero
in the vertical direction to have plane longitudinal wave propagation or in the
horizontal direction to have plane shear wave propagation. The excitation is applied
uniformly over the thickness of the model on one side of the area of study where an
“ultra” ALID is placed. On the other side is the layer studied, either an ALID or a PML.
The complex displacement is monitored along a line in the area of study and the
amplitude of the incident and reflected signal are determined by performing a spatial
FFT [61].
For other angles of incidence, another 2D plane strain model is defined as illustrated
in Figure 3.8.b. The area of study is triangular. On one edge is the layer studied. On the
other two, “ultra” ALID are placed. The excitation is uniformly applied on one edge.
The displacements are monitored along 3 lines aligned with the predicted propagation
direction of the incident and reflected waves. Spatial FFTs are performed to obtain
reflection coefficients. Models are non-dimensional and the material is defined so that
a)
b)
Area of study
PML or ALID
“Ultra” ALID
Excitation
Spatial FFT
Constraint
Figure 3.8. FE model used to validate the analytical models a) normal incidence model, b) angled
incidence model
56
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
the longitudinal and shear wavelengths are 2 and 1 units respectively: Young’s
modulus is 8/3, Poisson’s ratio is 1/3 and density is 1.
Shear and longitudinal wave cases are treated separately. The lengths L of the layers
are L=3 for longitudinal wave incidence and L=1.5 for shear wave incidence for PML
and ALID. The attenuation parameters as defined previously are Ax=1 and p=3 for
PML and CMmax=10π and p=3.
Performance of the two techniques should not be compared with each other using these
examples, as the layers are not optimized.
The PML analytical model for bulk waves is based on the knowledge of wave decay
in the PML [60] and the reflection/mode conversion occurring at the end of the PML.
The PML is defined in the x direction and αx is defined as in equation 3·9. Waves in
the PML have the following form:
p
u ( x, y, t ) = u 0 exp ( – iωt ) exp ( ik ( x cos θ + y sin θ ) ) exp ( – A x k cos θx ) (3·20)
where θ indicates the angle of incidence of the wave as shown in Figure 3.8. Note that
the decay rates of longitudinal and shear waves are unequal as their wave numbers k
differ. If the amplitude u0 of a wave entering the PML is 1, then the amplitude of the
wave after travelling the length of the PML is:
p+1
L
Amplitude = exp ⎛ – A x k cos ( θ ) -----------------⎞ (3·21)
⎝ ( p + 1 )⎠
At the free end of the PML, reflection and mode conversion occurs. Calculation of the
amplitude of each reflected wave is done using the stress free boundary conditions at
the end of the PML:
σ xx = 0 and σ xy = 0 (3·22)
57
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
The reflected waves exponentially decay again over the length of the PML following
equation 3·21. The theoretical reflection coefficient for bulk waves incident at a PML
over a range of angles for a given case can therefore be quickly calculated using this
method.
As mentioned before, in theory, no reflection should occur at the interface of, or inside,
the PML. In practice, numerical experiments have shown that this does occur and is
evident in extreme cases when a limited number of elements cannot accurately
represent a strong decay of the wave (spatial under-sampling) induced by a large value
of Αx. Figure 3.9.a presents the reflection obtained at normal incidence (θ=0) for a case
where a PML with a quadratic variation of α is used. Αx is varied from 0 to 50. For low
values of Αx the FE reflection coefficient at normal incidence matches the prediction
of the theoretical model but as Αx increases a non-negligible numerical reflection
occurs. One can see that the numerical reflection grows with Αx and is negligible for
low values of Αx where the theoretical reflection is dominant.
a) 0 b) 0
Reflection coefficient (dB)
Reflection Coefficient (dB)
-20 -20
-60 -60
-80
-80
-100
-100 Theoretical prediction
of the reflection coeffcient -120
-120
0 10 20 30 40 50 0 5 10 15 20 25 30 35 40
Ax Number of element per wavelength
Figure 3.9. a) Reflection coefficient against Ax b) Reflection coefficient against the number of
elements per wavelength.
It is also shown that the numerical reflection also depends on the number of elements
per wavelength. A series of models are run while only varying the number of elements
per wavelength. Results are shown in Figure 3.9.b. As the number of elements per
wavelength increases, the numerical reflection decreases. Further studies have shown
that the numerical reflection can also be decreased by increasing the polynomial order
of the shape function of the element or reducing the decay rate per wavelength. In most
cases, the element size and order are given by other parameters of the model and
altering them to reduce numerical reflection would decrease computational efficiency.
58
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
Starting by taking a layer whose length Lor is equal to the shortest wavelength, the
reflection coefficient is:
p+1
L
RC dB = 20 ⋅ log 10 ⎛ exp ⎛ – 2A x k -----------------⎞ ⎞ (3·23)
⎝ ⎝ ( p + 1 )⎠ ⎠
Numerical reflection will be minimum for the lowest value of Αx that satisfies the
acceptability criterion RCdB. Experience has shown that for this PML length, the
calculation can be based on the theoretical value and the numerical reflection can be
confidently neglected. Αx is taken equal to:
RC dB⎞
⎛ ⎛⎝ ------------
- ⎞
( p + 1 ) ⎜ 20 ⎠ ⎟
A x = – 0.5 ---------------------------
- log 10 (3·24)
p+1 ⎜ ⎟
k max L or ⎝ ⎠
Keeping this value of Ax, the modeller can define the length of the PML so that it meets
the acceptability criterion over the range of interest by using the theoretical reflection
coefficients obtained with the analytical model. Clearly, Lor influences the final length
of the PML as it limits the value of Ax. Modellers can reduce the value of Lor to less
than the advised value but should be cautious that numerical reflection does not pollute
their results.
A FE model using the parameters defined in Section 3.4.1.2 is used to verify the
validity of the analytical model. The comparison of these 2 models is shown in Figure
3.10.
Looking at Figure 3.10, one can see that as expected there is no mode conversion at
normal incidence (θ=0). The influence of the angle of incidence is highlighted as the
59
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
0 0
−60 −60
−80 −80
0 15 30 45 60 75 90 0 15 30 45 60 75 90
Longitudinal wave incidence (degree) Shear wave incidence (degree)
Analytical results: Longitudinal wave reflection FE results: L-wave
Shear wave reflection T-wave
Evanescent wave reflection
Figure 3.10. Reflection coefficient for a given PML obtained with bulk wave analytical and FE models
non converted reflection of the incident mode smoothly tends towards unity as the
angle of incidence tends towards 90 degrees for both wave types. The general shape of
these curves for waves at incidence above 60 degrees indicates that PMLs generally
cope well with high angle incidence. This confirms that the range of angle of incidence
plays a critical role in determining the layer’s parameters. One can also note the
presence of critical angles indicating that evanescent waves propagating along the back
surface of the PML and decaying perpendicularly to the free surface (likely to be a
Rayleigh wave) are generated. No reflection of this wave type has been noticed in the
area of study. Overall, there is an excellent agreement between the FE model results
and the prediction obtained with the analytical model. This confirms that the analytical
model can be used to define PML parameters.
This analytical model enables the quick evaluation of the acceptability of a given PML
for a particular model. Experience has shown that a small value of p, the variable
defining the power dependence, is usually preferable. Combined with previous
comments concerning this variable, it is recommended to use: 2 ≤ p ≤ 3 . No great gain
in model efficiency can be achieved by optimizing p. It is advised to only concentrate
on minimising L for the recommended Ax. In the cases where time domain results are
required, it is advantageous to take into account the shape of the frequency spectrum.
60
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
The ALID analytical model for bulk waves is based on a multi-layered approach. As
space is discretized in a FE model, the ALID is considered as a series of sub layers
(referred to as layers in this Section) with increasing damping. The variation of CM is
implemented by changing its value in each of these layers. Using layers one element
thick is the thinnest possible and clearly minimizes the change in damping between any
2 consecutive layers. It is of course not a problem to use layers several elements thick,
but best results were obtained with the thinnest layer possible. The ALID can be
analytically modelled using the continuity of displacements and stresses at the
interface of each of these layers. Such an analytical model will enable a modeller to
easily optimize the layer parameters.
Let us consider a multi-layered system. We define the displacement and stresses at the
front and back with the subscript f (front) or b (back). At the front of layer n, the
following equations are valid:
u x [ n ]f
u y [ n ]f u [ n ]f
= (3·25)
σ xx [ n ]f σ [ n ]f
σ xy [ n ]f
Φp[ n ]
Ψp[ n ] P[n ]
= (3·26)
Φn[ n ] N[ n ]
Ψn[ n ]
M [ n ]f = M P [ n ]f M N [ n ]f (3·27)
u [ n ]f P[ n ]
and = M P [ n ]f M N [ n ]f (3·28)
σ [ n ]f N[ n ]
61
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
At the interface between layer n and n-1, there is continuity of displacements and
stresses:
u [ n – 1 ]b u [ n ]f
= (3·29)
σ [ n – 1 ]b σ [ n ]f
Several techniques exist to deal with wave propagation through multi-layered systems
[37]. The technique used in this section is based on the one used by Pialucha [62]. It
does not use a transfer matrix but a global matrix algorithm; amplitudes of the waves
are calculated at each interface in the system. The strength of the technique is that
instabilities are removed by placing the origin of all waves at their entry point in the
layer.
ik Lx x ik Sx x
ik Lx – ik Sy – i k Lx e – i k Sy e
ik Lx x ik Sx x
ik Ly ik Sx ik Ly e – i k Sx e
M [ n ]f = '
2 2 2 2 ik x ik Sx x
– ( λk + 2μk Ly ) 2μk Sx k Sy – ( λk + 2μk Ly )e Lx – 2 μk Sx k Sy e
2 2 ik Lx x 2 2 ikSx x
– 2 μk Lx k Ly – 2μ ( k Sx – k Sy ) 2μk Lx k Ly e – 2μ ( k Sx – k Sy ) e
(3·30)
ik Lx x ik Sx x
ik Lx e – ik Sy e – i k Lx – i k Sy
ik Lx x ikSx x
ik Ly e ik Sx e ik Ly – i k Sx
M [ n ]b =
2 2 ik Lx x ik Sx x 2 2
– ( λk + 2μk Ly )e 2μk Sx k Sy e – ( λk + 2μk Ly ) – 2 μk Sx k Sy
ik Lx x 2 2 ik Sx x 2 2
– 2 μk Lx k Ly e – 2μ ( k Sx – k Sy ) e 2μk Lx k Ly – 2μ ( k Sx – k Sy )
(3·31)
u [ n ]b P[ n ]
= M [ n ]b (3·32)
σ [ n ]b N[ n ]
62
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
u [ n + 1 ]f P[ n + 1 ]
= M [ n + 1 ]f (3·33)
σ [ n + 1 ]f N[ n + 1 ]
As there is continuity of stresses and displacements - see equation 3·29, the following
equation is defined:
P[ n] P[ n + 1 ]
M [ n ]b = M [ n + 1 ]f (3·34)
N[ n ] N[ n + 1 ]
P[n ] P[ n + 1]
hence M [ n ]b – M [ n + 1 ]f = 0 (3·35)
N[ n ] N[ n + 1 ]
P[1 ] P[ 2]
M [ 1 ]b – M [ 2 ]f = 0 (3·36)
N[ 1 ] N[ 2 ]
P[2 ] P[ 3]
M [ 2 ]b – M [ 3 ]f = 0 (3·37)
N[ 2 ] N[ 3 ]
and so on...
P[1] is known as it is the model input. Equation 3·32 can be rewritten as:
P[ 1] P[2 ]
M P [ 1 ]b M N [ 1 ]b – M [ 2 ]f = 0 (3·38)
N[ 1 ] N[ 2 ]
P[2 ]
M N [ 1 ]b N [ 1 ] – M [ 2 ]f = – M P [ 1 ]b P [ 1 ] (3·39)
N[ 2 ]
63
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
σ [ n ]b = 0 (3·40)
P[ n ]
hence L [ n ]b = 0 (3·41)
N[ n ]
with
2 2 ik Lx x ik Sx x 2 2
– ( λk + 2μk Ly )e 2μk Sx k Sy e – ( λk + 2μk Ly ) – 2 μk Sx k Sy
L [ n ]b =
ik Lx x 2 2 ikSx x 2 2
– 2 μk Lx k Ly e – 2μ ( k Sx – k Sy ) e 2μk Lx k Ly – 2μ ( k Sx – k Sy )
(3·42)
N[ 1 ]
M N [ 1 ]b – M [ 2 ]f 0 … 0 P[2 ]
N[ 2 ] – M P [ 1 ]f P [ 1 ]
0 M [ 2 ]b – M [ 3 ]f … 0
… 0
… … … … 0 =
…
P[n – 1 ]
0 0 0 M [ n – 1 ]b – M [ n ]f 0
N[ n – 1 ] 0
0 0 0 0 L [ n ]b
P[n ]
N[ n ]
(3·43)
By solving this system, [N[1]] is obtained. It gives the amplitude of the wave potentials
reflected into the first layer (the area of study in the FE model). The amplitude of the
wave displacements is therefore derived from this data using the expressions defined
in Section 2.2.1. Finally, reflection coefficients are deduced from this.
64
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
Results for a sample case are compared with the FE models and are presented in Figure
3.11.
0 0
−20 −20
Reflection Coefficient (dB)
−60 −60
−80 −80
0 15 30 45 60 75 90 0 15 30 45 60 75 90
Longitudinal wave incidence (degree) Shear wave incidence (degree)
Analytical results: Longitudinal wave reflection FE results: L-wave
Shear wave reflection T-wave
Evanescent wave reflection
Figure 3.11. Reflection coefficient for a given ALID obtained with bulk wave analytical and FE
models
There is an excellent agreement between the FE and analytical results. As for the PML
case, one can note the importance of the range of angle of incidence in validating the
acceptability criterion. Although the ALID and PML used in this section should not be
quantitatively compared, the general shape of the curves in Figure 3.10 and 3.11 at high
angles of incidence indicates that the efficiency of absorbing layers for both non mode
converted reflected waves deteriorates faster for ALID than PML at these angles.
Therefore ALID can be considered as more sensitive to high angles of incidence than
PML. The presence of a critical angle for incident shear waves is also noticeable and
should be considered when validating the layer.
The excellent agreement between FE and analytical model confirms that the analytical
model can be used in the process of quickly determining the ALID parameters. When
defining an ALID, we recommended to use a value of p equal to 3 and find the right
combination of thickness and CMmax to minimize the thickness of the ALID.
For time domain models, the frequency spectrum of a signal is not constant. Taking
this into account relaxes the constraints to define an acceptable layer and therefore
leads to thinner layers.
65
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
The aim of the analytical model for guided wave cases is to evaluate the reflection
coefficient of guided waves from a given absorbing layer over a range of frequencies.
This part presents the case of guided waves in a 2D plate. The same approach could
readily be taken for wave guiding phenomena in pipes or cylinders.
i ( k x x – ωt )
u ( x, y, t ) = u x ( y )e (3·44)
p
x
– A x k x ----p-
i ( k x x – ωt ) L
u ( x, y, t ) = u x ( y )e e (3·45)
It is important to mention that PMLs have not been created for evanescent waves
(imaginary or complex wavenumber) and do not work well with them [63, 64]. At best,
the natural decay of these waves carry on as they propagate in the PML. Let us take the
example of a wave with an imaginary wave number kx=i.k with k a real number. In the
PML, the wave is described by:
p
x
– iA x k ----p-
– iωt – k x x L
u ( x, y, t ) = u x ( y )e e e (3·46)
The decay is the same as in the elastic area so the wave continues its normal decay in
the PML but now there is also a propagating part added to it. This part varies in the
PML and leads to a change in impedance. This can also lead to spatial under sampling
causing numerical reflection. This confirms that PML is not suited to evanescent
modes. A large amount of research has been done regarding this issue in the domain of
electromagnetism [65, 66, 67]. Evanescent waves by essence decay and usually do so
66
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
rapidly so they are usually very localised in practical cases. A pragmatic approach is
therefore to make sure that PML are placed so that only evanescent waves of negligible
amplitude (which have sufficiently decayed) reach them or that their potential
reflection from the PML will not influence the area of interest in the model.
Another type of guided wave which has proved to be problematic with PML is
backwardly propagating modes which in essence are waves whose group and phase
velocities are of opposing signs. A backward travelling wave incident at a PML will
have a negative wave number kx=-k with k a real positive number. In the PML, the
wave is hence described by:
p
x
A x k ----p-
i ( – k x – ωt ) L
u ( x, y, t ) = u x ( y )e e (3·47)
It is easy to see that this expression will lead to the wave growing exponentially inside
the PML, instead of decaying, which clearly is a problem.
Since these waves only occur over a very limited range of frequencies, which usually
is of no practical interest in NDE applications, one way of avoiding the problem is
simply to avoid these frequencies. This can be achieved easily in the PML implemen-
tation presented here as the models are solved in the frequency domain. A study
performed at Imperial College [68] showed that it is possible to fix the problem by
performing a modal decomposition using the bi-orthogonality relations of the guided
modes [69] and channelling them to separate PMLs defined such that all wave types
are absorbed correctly. This is an attractive approach in cases where access to the
coding is possible, but is out of reach for the standard packages used in this study. If
the frequency range of study includes the possibility of backwardly propagating waves,
one pragmatic alternative is to use ALID instead of PML as ALID does not have any
difficulty absorbing backwardly propagating waves.
In order to confirm the validity of the analytical models, their results are compared with
results from FE models. 2D plane strain models shown in Figure 3.12 are used. As for
bulk waves, the use of “ultra” ALID is necessary. The excitation is applied over the
67
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
thickness of the model on one side of the area of study where an “ultra” ALID is placed.
On the other side is the layer studied, either an ALID or a PML. The displacement is
monitored along a line in the area of study and the amplitude of the incident and
reflected signal are determined by performing a 2D FFT.
a)
Figure 3.12. FE model used to validate the guided wave analytical models.
The plate is 8mm thick and made of aluminium (ρ=2780kg/m3, E=70GPa, ν=0.33). It
is excited at frequencies ranging from 114kHz to 186kHz. At these frequencies, the
two fundamental modes, A0 and S0, are studied. For these cases, the PMLs are defined
with L= 13.1mm, Ax=2.66e6 and p=3 for S0 and L= 6.55mm, Ax=2.66e6 and p=3 for
A0. The ALID are defined with L= 13.1mm, CMmax=2.66e6 and p=3 for S0 and L=
76.5mm, CMmax=1e6 and p=3 for A0.
Performance of the two techniques should not be compared using these examples as
the layers are not optimized.
The analytical model used to determine the reflection coefficient of a PML is based on
the knowledge of the decay of the mode in the layer and the reflection of modes at the
end of the layer.
Propagating modes are treated using a model similar to the PML bulk wave model. One
difference is that any mode conversion that may occur at the end of the PML is not
taken into account in the analytical model and total reflection is assumed. Each
propagating mode is considered independently. The reflection coefficient for each
propagating mode is equal to:
p+1
L
RC dB = 20 ⋅ log 10 exp ⎛ – 2A x k -----------------⎞ (3·48)
⎝ ( p + 1 )⎠
As for the case of bulk waves, the issue of numerical reflection needs to be taken into
account. In a given model, the wave which suffers the most from numerical reflection
68
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
is the one with the highest wave number (and hence the shortest wave length) as it is
the most strongly attenuated (i.e.: the decay occurs on fewer elements than the other
waves).
As a starting point, a layer whose length Lor is equal to the shortest wavelength is taken.
Numerical reflection will be the lowest for the lowest value of Ax that satisfies the
acceptability criterion RCdB. Experience has shown that for this PML length, the
calculation can be based on the theoretical value and the numerical reflection can be
confidently neglected. Αx is taken equal to:
RC dB⎞
⎛ ⎛⎝ ------------
- ⎞
(p + 1) 20 ⎠ ⎟
- log ⎜ 10
A x = – 0.5 --------------------------- (3·49)
p+1 ⎜ ⎟
k max L or ⎝ ⎠
Keeping this value of Ax, all the other waves will be less strongly attenuated. The wave
with the longest wavelength (and so the shortest wave number) will have the slowest
decay in the layer. In order to correctly define the PML over the range of wave
numbers, the length of the layer for the previously defined Ax is increased to achieve
the theoretical reflection criterion using the following formula:
1 -
RC dB⎞ p + 1 -----------
⎛ ⎛ ⎛⎝ ------------
- ⎞⎞
⎜ ( p + 1 ) ⎜ 20 ⎠ ⎟ ⎟
L PML = ⎜ – 0.5 ----------------- log ⎜ 10 (3·50)
k min A x ⎟⎟
⎝ ⎝ ⎠⎠
Using equation 3·49 and 3·50, the length of the PML is:
1 -
-----------
k max p + 1
L PML = L or ⎛ ----------⎞ (3·51)
⎝ k min ⎠
Validating that the reflection coefficient is acceptable for each propagating mode
ensures that the PML is correctly defined. As mentioned before, this formula is only
valid for propagating modes. Evanescent and backwardly propagating modes need to
be dealt with separately.
Clearly, Lor directly influences the final length of the PML but this can be reduced to
less than the advised value if required although caution is necessary in this case. The
69
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
formula also highlights that the final length of the PMLs depends on the ratio of the
extreme wave numbers. High ratios typically occur at frequencies just above cut-off
and lead to potentially long PMLs.
As shown previously, evanescent waves are not correctly dealt with by PMLs and it is
necessary to make sure they do not have a significant impact on the displacement field
of the area of study. The evanescent wave with the smallest imaginary part k″ min
(which has the lowest decay rate) is considered. Its amplitude is conservatively taken
equal to 1 at any feature which can generate evanescent waves (e.g. defect, source). It
is also conservatively considered that an evanescent wave incident at a PML would be
totally reflected. In order to have an acceptable displacement field between the PML
and the feature, the wave needs to have decayed sufficiently (i.e. its amplitude is equal
to RCdB) when it reaches the PML.
The distance between the feature and the PML needs to be:
RC dB⎞
⎛ ⎛⎝ ------------
- ⎞
1 20 ⎠ ⎟
L = – ------------ log ⎜ 10 (3·52)
k″ min ⎜ ⎟
⎝ ⎠
In cases where the displacement field between the PML and the feature is not
monitored, one can use:
RC dB⎞
⎛ ⎛⎝ ------------
- ⎞
1 ⎜ 20 ⎠ ⎟
L = – ---------------
- log 10 (3·53)
2k″ min ⎜ ⎟
⎝ ⎠
In these cases, the wave amplitude is equal to RCdB as it reaches the feature after being
reflected and is therefore acceptable.
Results from the analytical model are compared with FE results obtained for the first
fundamental propagating modes A0 and S0 with the parameters defined in Section
3.4.2.2. Figure 3.13 shows that there is an excellent agreement between both sets of
results and validates the use of the analytical approach to dimension PMLs for guided
waves.
70
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
0 0
−40 −40
S0
A0
−60 −60
−80 −80
120 130 140 150 160 170 180 120 130 140 150 160 170 180
Frequency (kHz) Frequency (kHz)
Analytical: S0 A0 FE: S0 A0
Figure 3.13. Reflection coefficient for a given PML obtained with guided wave analytical and FE
models
The following analytical model determines the reflection coefficient of ALID for
guided waves in plates modelled in 2D plane strain. It is based on a multi-layered
approach and makes use of an impedance matrix.
An infinity of modes can exist in a given guided wave system. A finite number of
modes are selected. This selection is based on the characteristics of the wave number
of the mode in the elastic part of the model. All real (propagating) and imaginary
(evanescent) modes and a finite number of complex modes (inhomegeneous) are
included in the model. As explained before, each layer in the ALID has the same elastic
properties as the area of study save for having a gradually increasing damping. As for
bulk waves, the damping causes the properties of the waves to vary as they propagate
through the ALID. The actual properties of modes existing in each layer inside an
ALID need to be calculated by solving the Rayleigh equation using the properties of
each layer. All modes have a complex wave number in the ALID. As for bulk waves,
the real part of the wave number varies as the damping is increased, leading to a change
of acoustic impedance inside the layer and therefore causing reflection. Contrary to the
PML, due to the physical properties of the ALID materials, all modes including
backwardly propagating modes decay in the ALID [31].
71
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
In guided wave cases as in bulk wave cases, ALID can be seen as a series of finite
layers with varying properties. In 2D, the symmetric and anti-symmetric modes can be
treated independently as ALID are symmetric with respect to the mid plane of the plate.
The model defined in this section is a development of a model by Pagneux [70]. The
idea is to define a global impedance matrix of the ALID that can be used to determine
reflection coefficients. The matrix is defined starting from the end of the ALID and
propagated through each layer and each interface between these layers until the area of
study. Let us consider a multi layered system as represented in Figure 3.14 where layer
1 represents the area of study.
y
x xb xf xb xf xb xf xb xf xb x
f xb Front or back face
One starts by noting that at any point in the plate and ALID the total displacements and
stresses (subscript tot) are equal to the sum of the displacement and stresses of the
modes travelling in the positive (po) and the negative (ne) x-direction. Combining this
with the symmetry properties of the Lamb wave mode shape, we have:
u x ( tot ) ( x, y ) u x ( po ) ( x, y ) u x ( ne ) ( x, y )
u y ( tot ) ( x, y ) u y ( po ) ( x, y ) u y ( ne ) ( x, y )
= + (3·54)
σ xx ( tot ) ( x, y ) σ xx ( po ) ( x, y ) σ xx ( ne ) ( x, y )
σ xy ( tot ) ( x, y ) σ xy ( po ) ( x, y ) σ xy ( ne ) ( x, y )
n
u x ( po ) ( x, y ) ux ( y )
n n
u y ( po ) ( x, y ) uy ( y )
with
σ xx ( po ) ( x, y )
= ∑ An ( x ) n
(3·55)
1 σ xx ( y )
σ xy ( po ) ( x, y ) n
σ xy ( y )
72
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
n
u x ( ne ) ( x, y ) –ux ( y )
n n
u y ( ne ) ( x, y ) uy ( y )
and
σ xx ( ne ) ( x, y )
= ∑ Bn ( x ) n
(3·56)
1 σ xx ( y )
σ xy ( ne ) ( x, y ) n
– σ xy ( y )
n n n n
u x ( tot ) ( x, y ) ux ( y ) u y ( tot ) ( x, y ) uy ( y )
σ xx ( tot ) ( x, y )
= ∑ an ( x ) n
and
σ xy ( tot ) ( x, y )
= ∑ bn ( x ) n
(3·57)
1 σ xx ( y ) 1 σ xy ( y )
T T
One defines a ( x ) = a 1 ( x ) … a n ( x ) and b ( x ) = b 1 ( x ) … b n ( x ) and the
impedance matrix as:
a ( x ) = Z ( x )b ( x ) (3·58)
Inside a layer, the impedance matrix can be propagated between the back (b) and the
front (f) of a layer d. One has:
d d d d
A n ( x f ) = A n ( x b ) exp ( ik n ( x f – x b ) ) and B n ( x f ) = B n ( x b ) exp ( ik n ( x f – x b ) ) (3·59)
therefore:
d d
a ( xf ) a ( xb )
= C iS (3·60)
d
b ( xf ) iS C b d ( x )
b
d –1
Z ( x f ) = iS + CZ d ( x b ) C + iSZ d ( x b ) (3·61)
73
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
This formulation shows instabilities [70] and the following equation is a better
alternative:
d –1 –1 d –1 –1 –1
Z ( x f ) = – iT +S Z ( xb ) –i T S (3·62)
At the boundary between the back face of layer d and the front face of layer d+1, one
has:
n n, d n n, d + 1
d ux ( y ) d+1 ux (y)
∑ an ( xb )
n, d
= ∑ an ( xf )
n, d + 1
(3·63)
1 σ xx ( y ) 1 σ xx (y)
n n, d n n, d + 1
d – σ xy ( y ) d+1 – σ xy (y)
∑ bn ( xb )
n, d
= ∑ bn ( xf )
n, d + 1
(3·64)
1 uy ( y ) 1 uy (y)
T
These expressions are multiplied by m, d m, d and
– σ xy ( y ) uy ( y )
T
m, d + 1 m, d + 1 respectively and integrated over the thickness of the plate:
ux ( y ) σ xx (y)
n h n, d m, d n h n, d + 1 m, d
d ux ( y ) – σ xy ( y ) d+1 ux (y) – σ xy ( y )
∑ an ( xb ) ∫ n, d
⋅
m, d
dy = ∑ an ( xf ) ∫ n, d + 1
⋅
m, d
dy
1 – h σ xx ( y ) uy ( y ) 1 – h σ xx (y) uy ( y )
(3·65)
n h m, d + 1 n, d n h m, d + 1 n, d + 1
d ux (y) – σ xy ( y ) d+1 ux (y) – σ xy (y)
∑ bn ( xb ) ∫ m, d + 1
⋅
n, d
dy = ∑ bn ( xf ) ∫ m, d + 1
⋅
n, d + 1
dy
1 – h σ xx (y) uy ( y ) 1 – h σ xx (y) uy (y)
(3·66)
One has:
d d, d d+1 d, d + 1
a ( x b )J n, m = a ( x f )J n, m (3·67)
d d + 1, d d+1 d + 1, d + 1
b ( x b )J m, n = b ( x f )J m, n (3·68)
74
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
–1
d d+1 d, d + 1 d, d
a ( xb ) = a ( x f ) J n, m J n, m (3·69)
–1 –1 –1
d = bd + 1 ( x ) d + 1, d + 1 d, d + 1 (3·70)
b ( xb ) f J m, n J m, n
d –1 d+1 –1 –1
Z ( xb ) = ad ( x ) bd ( x ) = Z ( x f ) J d, d + 1 J d, d d + 1, d + 1
J m, n
d + 1, d
J m, n
b b n, m n, m
(3·71)
d, d d, d
J n, m = J n, m δ n, m (3·72)
The final missing part to describe an ALID is the determination of the behaviour at the
end of the last layer. At this point, there is a stress free condition:
end end
σ xx ( tot ) ( x b ) = 0 and σ xy ( tot ) ( x b ) = 0 (3·73)
n n, end n n, end
end σ xx (y) end σ xx (y)
∑ An ( xb )
n, end
+∑ Bn ( xb )
n, end
= 0 (3·75)
1 σ xy ( y ) 1 – σ xy ( y )
n n
0 n, end
+ ∑ b n ( x b ) σ xx (y) = 0
end end
∑ an ( xb ) n, end
σ xy ( y )
(3·76)
1 1 0
75
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
T
One multiplies equation 3·76 by: u m, end ( y ) u m, end ( y ) and integrates through the
x y
thickness
n h m, end n h m, end
0 ux (y) n, end (y)
end end σ xx ( y ) ⋅ ux
∑ an ( xb ) ∫ n, end
σ xy (y)
⋅
m, end
dy + ∑ bn ( xb ) ∫ m, end
dy = 0
1 –h uy (y) 1 –h 0 uy (y)
(3·77)
end end
a ( x b )G n, m + b ( x b )H n, m = 0 (3·78)
end –1
Z ( x b ) = – H n, m G n, m (3·79)
Zend(xb) is used in equation 3·62 to obtain Zend(xf) which is used to obtain Zend-1(xb)
with equation 3·71 and so on until Z1(xb)=ZALID, which represents the total impedance
of the ALID.
1 T
The following expressions are defined: A ( x b ) = A 1 ( x b ) … A n ( x b ) and
1 T
B ( x b ) = B 1 ( x b ) … B n ( x b ) . At the entrance of the ALID, A1(xb) is proportional
to the amplitude of the modes incident to the ALID and B1(xb) to the reflected modes:
–1
1 1 1
B ( x b ) = I – Z ALID I + Z ALID A ( x b ) = RC A ( x b ) (3·80)
In order not to over-complicate the model, each mode is studied separately and it is
conservatively assumed that each mode has a unit amplitude as it reaches the ALID.
The diagonal of matrix [RC] gives the ratio of the reflected to incident amplitude for
each mode (propagating and evanescent). This value is used to evaluate the compliance
of all modes with the acceptability criterion. As all modes are considered, the ALID
can be placed very close to any feature (e.g. defect, source) in the model.
Results obtained with this analytical model are compared with FE models with the
parameters defined in Section 3.4.2.2. Results for propagating modes are presented in
Figure 3.15. There is an excellent agreement between both techniques so the analytical
models can be confidently used to determine parameters of the ALID.
76
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
0 0
−40 −40
A0
−60 −60 S0
−80 −80
120 130 140 150 160 170 180 120 130 140 150 160 170 180
Frequency (kHz) Frequency (kHz)
Analytical: S0 A0 FE: S0 A0
Figure 3.15. Reflection coefficient for a given ALID obtained with guided wave analytical and FE
models
3.5. Demonstrators
These demonstrators show the possibilities and gains that can be achieved using
absorbing layers.
This model is a realistic simplified practical bulk wave case represented in Figure
3.16.a. The material used is steel (ρ=7800kg/m3, E=200Gpa, ν=0.33). The model is
excited by a point force on the top free surface acting at 45 degrees to the normal (in
order to extend the generality of the model). It is applied in the form of a 2MHz 5
cycles tone burst (Blackman Harris 3 terms). A simplified horizontal 4.5mm crack is
located 30mm under and to the right of the excitation point. The aim of this is to
monitor the interaction of the waves with the free surface and the defect by monitoring
the displacement at a chosen point. In this example, the monitoring point is located at
the excitation point.
77
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
a) b)
100mm
200mm
c) d)
33.5mm
58.5mm
Area of study
Absorbing layer
66.5mm 13.5mm
Figure 3.16. a) bulk wave demonstrator, FE model: b) without absorbing layer, c) with ALID, d) with
PML.
Without using absorbing layers, the model needs to be 100mm by 200mm (Figure
3.16.b) in order to achieve separation of the first reflections from the defect and the
unwanted reflections at the monitoring point.
For models using absorbing layers, the acceptability criterion is set at -60dB meaning
that the maximum reflection from the absorbing layers will be smaller than 0.1% of the
maximum displacement of the excitation point. Using the models described previously
in this chapter to define the required PML or ALID and taking into account the
frequency spectrum of the excitation, the model size can be reduced to 66.5mm by
56mm (Figure 3.16.c) using ALID and 33.5mm by 13.5mm (Figure 3.16.d) using
PML. The thickness of the ALID is 16mm, CMmax is 11e6 and p is 3. The PML is 3mm
thick, Amax is equal to 2.3e9 and p is 3. The thickness of the PML is smaller than the
ALID one because of the perfectly matched impedance. It should also be noted that the
area of study is smaller for the PML case than for the ALID case because ALID is more
sensitive than PML to waves incident at high angles. The geometric size of the 2D
model is divided by 5.35 with ALID and 44 with PML compared with the classical
78
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
technique. This reasoning can be extended to a 3D model of the same case where the
model size would be divided by 15 and 655 respectively.
An added benefit of models using absorbing layers is that the complete area of study
is clean of any unwanted reflection. This can help a better understanding of the
physical phenomena occurring in the model. Figure 3.17 shows the displacement field
obtained for the ALID case described above. The absolute value of the displacement is
represented and, in order to highlight low amplitude reflection, the colour scale only
spans up to 0.1% of the maximum absolute displacement in the model. These figures
show that there is no visible reflection coming back from the ALID into the area of
study. This confirms the correct dimensioning obtained using the analytical model.
a) c)
b) d)
Figure 3.17. Absolute displacement field for the bulk demonstrator with ALID at time: a)5μsec
b)10μsec c)15μsec d)20μsec. Colour scale extends from 0 (blue) to 0.1% (red) of the maximum
absolute displacement. Grey indicates out of scale (0.1% to 100%). White dashed line indicates the
boundary between area of study and ALID
79
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
The guided wave model represented in Figure 3.18.a is a simplified practical case
which has been used extensively in guided wave research. An 8mm thick aluminium
(ρ=2780kg/m3, E=70Gpa, ν=0.33) plate is excited by a 150kHz 12 cycles tone burst
(Blackman Harris 3 terms) with a point force at 45 degrees in order to strongly excite
both fundamental modes. Over the frequency range excited, the fundamental
symmetric and anti-symmetric modes can propagate in the medium. The plate contains
a 2mm by 2mm slit located 700mm from the excitation. If the first reflections from the
defect were studied without using absorbing layers, the model would need to be
2700mm long (Figure 3.18.b) to avoid the presence of unwanted reflections reaching
the monitoring point before the reflections from the slit.
Using the analytical models devised above and taking into account the frequency
spectrum, an optimized ALID is 80mm thick. In this case, CMmax is 2.1e6 and CM
varies following a cubic law. The left hand side ALID starts at the excitation point and
the right hand side ALID 3mm behind the slit. The model size is divided by 3.12
(Figure 3.18.c).
Using the analytical model, an optimized PML is 18.5mm thick. Amax is 8.36e5 and α
varies following a cubic law. As only the field between the slit and the excitation point
is of interest, it is calculated (taking advantage of the shape of the frequency spectrum)
that the PML needs to be located 29mm away from these features. The model size is
hence divided by 3.38 with PML (Figure 3.18.d).
2mm
a) 150kHz 700mm Area of study
Absorbing layer
8mm
b) 2700mm
865mm
c)
797mm
d)
Figure 3.18. a) guided wave demonstrator, FE model: b) without absorbing layer, c) with ALID, d)
with PML.
Figure 3.19 shows the displacement field for the ALID case described above. A scale
factor of 10 is applied in the vertical direction in order to improve the visibility. The
80
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
colour scale is adjusted to highlight the low amplitude part of the absolute
displacement field. Figure 3.19.a shows the 2 modes propagating in the plate. On
Figure 3.19.b one can see the reflection of these 2 modes from the defect. As mode
conversion occurs at the defect, 4 wave packets exist. On Figure 3.19.c the last wave
packet is about to enter the ALID and the absence of any other visible wave packets
reflected from the ALID can be noted, although some numerical noise in the form of
low amplitude long wavelength waves can be observed. Other models were run and
confirmed that this noise is not related to the ALID. Finally, Figure 3.19.d shows that
all modes have been correctly absorbed indicating that the ALID worked correctly.
Figure 3.19. Absolute displacement field for the guided demonstrator with ALID at time: a)150μsec
b)300μsec c)450μsec d)600μsec. Colour scale is varied and extends from 0 (blue) to 2% or 10% (red)
of the maximum absolute displacement as indicated on the figure. Grey indicates out of scale (2% or
10% to 100%). White dashed line indicates the boundary between area of study and ALID.
As mentioned in the introduction, our implementation of the PML can only be used in
the frequency domain. Wave propagation time domain results can be obtained by using
absorbing layers (either ALID or PML) and running the model in the frequency
domain. Results are then processed with an inverse Fourier transform. The aim of this
case is to describe the method used to do this, using an example of modelling of guided
waves. This case was developed with Ludovic Moreau and formed part of a publication
[31].
81
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
In this case, the first anti-symmetric mode A0 is generated in an 8mm thick aluminium
plate with a 2mm by 2mm rectangular notch as shown in Figure 3.20. The excitation
time signal desired is a 150kHz 10 cycle tone burst applied in the normal direction to
the plate. Frequency domain input data is defined by FFT of the time excitation shown
in Figure 3.21. Most of the energy is concentrated between 120kHz and 180kHz. Only
frequencies located in this range will be used to solve the model.
Excitation
Amplitude
φ(f)
Magnitude
A(f)
Phase
A0 mode generation is achieved by applying the input data as a force in the normal
direction over a region on the top and bottom surfaces of the plate, in a spatially
windowed sine wave whose wavelength matches the wavelength (obtained using
DISPERSE for example) of the desired mode. This is shown in Figure 3.22. For this
case, ALIDs are used to absorb incident waves but PMLs could have been used
similarly. The complex normal displacements for each frequency are monitored
700mm away from the defect and an inverse FFT is applied to this complex spectrum.
The real part of the inverse FFT result gives the time traces. In order to validate the
method, they are compared with time traces obtained with a classical time domain
model.
82
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
a) 45 b) 135mm
λS0(f)
λ(mm)
λmax
A(f)
λA0(f) φ(f)
0
120 Freq. (kHz)180 λA0(f)
Figure 3.22. a) dispersion curve data used for input definition, b) input definition
Results are presented in Figure 3.23. At the monitoring point, results from the different
models agree with high accuracy in terms of propagation velocities and relative
amplitudes. Differences on both of these parameters are smaller than 1%. This
confirms that it is possible to obtain wave propagation time domain results by using
absorbing layers and solving the model in the frequency domain.
Reflected A0 wave
0 Time (ms) 0.6
b)
Amplitude
c)
Amplitude
83
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
Using absorbing layers in models designed to study wave scattering can greatly
simplify them in addition to reducing their size compared to classical models.
The model used in this section is similar to the one used in [72]. The plate used is 3mm
thick and made of steel as defined in the published paper [72]. A notch 0.5mm deep is
positioned 1m away from the excitation. The width of the notch is varied from 0.5mm
to 5mm in separate analyses. The COMSOL FE model is represented in Figure 3.24. It
is excited by an harmonic normal force whose amplitude is constant through the
thickness of the plate. The frequency of the harmonic excitation is varied from 400kHz
to 500kHz. ALID are placed on each side of the area of study to absorb any incoming
wave.
Figure 3.24. Representation of model used for guided wave scattering validation
−20
Reflected peak Incident peak
Amplitude (dB)
−40
A0ref A0inc
−60
−80
−100
−120
−60 Wavenumber K (1/m) 60
84
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
A0 ref
RC = ------------- (3·81)
A0 inc
The notch width varies from 0.5 to 5mm and the wavelength from 5.17mm to 6.16mm.
The notch width is therefore varied from 8.1% to 96.7% of the wavelength. In Figure
3.26, the reflection is plotted against the notch width. Results from the COMSOL
models are compared with the published experiments and published finite element
models (FINEL [73]).
0.2
0.18
0.16
0.14
Reflectioncoefficient RC
0.12
0.1
0.08
0.06
0.04
0.02
0
0 10 20 30 40 50 60 70 80 90 100
Notch width (% of input wavelength)
Published experimental results [72] COMSOL FE results Published FINEL FE results [72]
It can be seen that results obtained from time domain modelling with FINEL and
frequency domain modelling with COMSOL agree very well. Incidentally, the
agreement of the COMSOL results with experimental ones is similar to the agreement
of FINEL results with experimental ones.
85
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
The aim of this case is to show that multi-mode reflection coefficients from a scatterer
can be easily obtained using frequency domain analysis and absorbing regions. This
case was developed with Ludovic Moreau and formed part of a publication [31].
The model set up is similar to the previous case although the dimensions are changed
and the excitation is defined so that only the fundamental anti symmetric mode is
excited as was done in section 3.5.2. The plate is an 8mm thick aluminium plate with
a 2mm by 2mm notch. The excitation frequency is varied from 140kHz to 500kHz by
steps of 15kHz. ALID are placed on each side of the area of study. The normal
displacement is monitored on the top surface of the plate over a region whose length is
equal to 3 times the longest wavelength in the model. It is used to perform a spatial FFT
[61] which gives the relative amplitude of the normal displacements for each mode. As
the mode shapes of each mode is different, it is necessary to use power-normalized
mode shapes in order to obtain the reflection coefficients in terms of energy.
For a particular mode, the value of the energy Ei carried by each mode [33]:
Ai 2
E i = ⎛ ------⎞ (3·82)
⎝ M i⎠
with Ai, the amplitude of the displacement in a given direction at a through thickness
position given by the spatial FFT and Mi, the value of the displacement in the same
through thickness position and direction of the power normalized mode shape of the
mode of interest. In this way, it is possible to verify that the summation of the energies
carried by the different modes reflected and transmitted from a defect is equal to the
energy of the incident mode. Reflection coefficients are calculated by dividing the
value of the energy carried by each mode by the value of the energy carried by the
incident mode (A0 in this case).
Results obtained with ABAQUS and COMSOL (both using ALID) are shown in
Figure 3.27. The evolution of the reflection coefficient from an A0 incident mode for
frequencies ranging from 140kHz to 500kHz is presented for each possible mode. They
agree well, showing the same patterns, the maximum difference over the range of
86
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
values being 3.2%. Similar results were obtained for transmission coefficients
calculation.
60
Reflection coefficient(%)
COMSOL results
0
140 Frequency (kHz) 500
60
ABAQUS results
Reflection coefficient(%)
0
140 Frequency (kHz) 500
A0 S0 A1 S1 S2 S-2
Figure 3.27. Energy reflection coefficient for A0 incident on a 2mm square notch in an 8mm thick
aluminium plate from 140kHz to 500kHz
In this case, as several modes having a large range of velocities are present, the use of
absorbing regions greatly simplifies the obtention of the reflection curves.
87
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
3.6. Discussion
The previous demonstrators have highlighted the gains obtained by using ALID or
PML and have also shown that PML generally achieves smaller model size than ALID.
This can be explained by the fact that with PML the layer matches the impedance of
the area of study perfectly whereas for the ALID the strength of the decay is limited by
the impedance change it causes. As a consequence, for a given case, ALID is generally
thicker than PML. It is also more sensitive to waves incident at high angles leading to
a necessary increase of the area of study in order to minimize the model size. In short,
PML is superior to ALID from a theoretical point of view for all wave types except
evanescent and backwardly ones. As the aim is to improve modelling capabilities, this
gives an advantage to the PML technique but it has to be balanced by the fact that, to
our knowledge at the time of writing, COMSOL is the only mainstream FE package
that allows implementation of PML. The implementation is limited to frequency
domain solving with implicit solvers whereas ALID can be implemented to be used
with any FE package which allows Rayleigh damping. The actual performance
comparison depends on the type of result required.
A user wishing to obtain frequency domain results such as the reflection coefficient
from a defect has a choice of directly using either PML or ALID. The two best
combinations investigated are ABAQUS with ALID and COMSOL with PML.
ABAQUS models were solved using the “Direct-solution steady-state dynamic
analysis” procedure of ABAQUS/Standard which uses a direct implicit solver.
COMSOL offers a range of direct and iterative implicit solvers. Direct solvers have
proved to be robust and efficient but limited in terms of the number of degrees of
freedom that can be solved. Iterative solvers showed improvement in terms of speed
and memory but convergence proved difficult. Overall, in our experience, using
COMSOL 3.2 and ABAQUS 6.6, COMSOL did not match the performance of
ABAQUS for a given number of degrees of freedom. Given these two combinations it
is therefore not possible to clearly recommend one of them, as the best solution
depends on the specific details of the model considered.
For time domain results, the most straightforward approach is to use ABAQUS/
Explicit with ALID. The explicit solver of ABAQUS uses the central difference
88
Chapter 3
Modelling waves in unbounded elastic media using absorbing layers
algorithm. It is particularly attractive for large models as it is faster and more memory
efficient than an implicit solver. This technique is recommended for most cases where
time domain results are required. An alternative is to solve the model in the frequency
domain and reconstruct the time domain results by performing an inverse fourier
transform.
3.7. Conclusion
89
Chapter 4
On the influence of mesh parameters on
elastic bulk wave velocities
4.1. Introduction
In the process of moving from the genesis of a numerical model to its conclusion and
the obtention of results enabling a better understanding of a physical phenomenon, a
modeller takes a series of design decisions. In the context of wave propagation
modelling, the choice of the solving technique, mesh density or time step influences
the successful outcome of the exercise but also the level of accuracy with which the
phenomenon is represented. Approximation by discretization of the space and time
leads to a discrepancy between theoretical physical values and their actual values in the
model. One of the parameters suffering from this is the propagation velocity of bulk
waves in elastic media. This study looks into how this is influenced by design decisions
and aims to provide modellers with information enabling them to improve not only the
quality of their modelling (adequate decision making) but also the quality of the
analysis made based on the numerical results (quantitatively taking into account the
numerical deviation). There is also an interest in understanding the reasons for mesh
scattering where a change in mesh size causes some numerical reflection. This is an
important point as a better understanding of this could lead to vast increase of
modelling capabilities by enabling local mesh refinement to be implemented.
In this chapter, both explicit and implicit solvers used with a range of element types are
studied. Regular meshes made of square or equilateral triangles are investigated first.
These meshes are then distorted in order to understand the influence of element
deformation on wave velocities. The general approach is to run a series of FE models.
The results of these models enables the identification of simple functions to describe
the errors as generally as possible and enable users to establish guidelines to suit their
modelling requirements.
90
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
Details from all the studies performed are presented in this chapter. This is done in
complete details in order to provide a useful general reference to users for future use.
Being aware that providing such a level of details for all cases studied does not help
with the flow of the chapter, it is suggested for users interested in a particular case to
investigate findings for this particular case and then to proceed to general discussion
and conclusions. Users interested in the principles will find them described in the
following sections: 4.2.1 (explicit solving), 4.3.1 (implicit solving), 4.2.2.1 (square
elements) and 4.2.3.1 (equilateral triangle elements).
4.2.1 Introduction
FE models are run for each case and post processing leads to the obtention of a precise
value of actual propagation velocity in the model. For simplicity, all models are non-
dimensional: density ρ=1, frequency f=1, Young’s modulus E=8/3, Poisson’s ratio
η=1/3 so that we have a longitudinal velocity cL=2 and shear velocity cS=1. Models are
made with a regular mesh of elements and represent an area covering at least 20
wavelengths away from the excitation point (different model sizes were taken for
longitudinal and shear waves) as shown in Figure 4.1.
Displacement
monitoring area
ALID
Excitation
point
Area of study
Models are 2D plane strain. Both shear and longitudinal wave generations are
considered, as models are excited by nodal forces being applied on four elements as
91
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
shown on Figure 4.2. This generates an almost pure wave - shear or longitudinal - and
as the size of the source is only a fraction of a wavelength, it is close to a point source.
The point located at the centre of the four elements is referred to as the excitation point.
a) b) c) d)
Figure 4.2. a) Longitudinal and b) shear wave excitation for a square element mesh and c) longitudinal
and d) shear excitation for a triangular elements mesh
The time domain excitation signal is a sine function of unit frequency. Displacements
are monitored over a disc area centred on the excitation point of 18 wavelengths radius.
Spatial FFTs are performed on the displacements over a range of radii at different
angles in a similar way to the example presented in Section 3.5.3.1. On these radii, only
the displacements from 3 to 18 wavelength from the centre are used. The signals used
for the spatial FFTs are Hanning windowed and padded in order to improve the
velocity precision. As some coarse meshes are used, cubic interpolation is performed
in order to increase the ability to look correctly at the influence of the direction of
propagation.
ALID are placed outside the area of study in order to absorb incoming waves.
This study looks into the velocity error for models made with square elements. The side
of the square element is L=L0=L90 as shown in Figure 4.3.
92
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
L90
L45
Lθ
θ
L0
Figure 4.3. Schematic defining L0, L90, L45 and Lθ in a mesh of square elements.
The mesh density N is expressed in term of nodes per wavelength and is varied from 6
to 30 along the side of the element:
λ λ
N = N0 = ------ = N90 = --------- (4·1)
L0 L90
In this chapter, the mesh density is defined relative to the considered wave. As the
longitudinal wavelength is twice the shear one, for a given mesh density, the element
size for longitudinal waves is twice that of shear ones.
In the first series of models run, the Courant number CFL as defined in Section 2.3.1
is varied between 1 and 0.05. Figure 4.4.a and 4.4.b represent the variation of the error
on the longitudinal velocity EcL and shear velocity EcS along the element side (at 0
degrees) against CFL for various mesh densities N. Some shear wave models did not
complete when run at the stability limit (CFL=1). The maximum CFL used for shear
waves is 0.98.
On Figure 4.4.a, it can be seen that the error on the longitudinal wave EcL is nil when
the CFL is equal to 1. This is mentioned in [38] and can be explained by the fact that,
in this case at 0 and 90 degrees, there is a match between the time stepping frequency
and the highest frequency of a single element:
ΔL 2
Δt = ------- = ------ (4·2)
cL ωn
2c L
hence ω n = -------- = ω max (4·3)
ΔL
93
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
5
a) N=6
3
8
2 9
10
11
12
1
15
20
30
0
b) 0 0.2 0.4 0.6 0.8 1
CFL
5
N=6
4
Shear velocity error EcS (%)
3
8
2 9
10
11
12
1
15
20
30
0
0 0.2 0.4 0.6 0.8 1
CFL
Figure 4.4. a) Longitudinal and b) shear velocity errors against CFL for various mesh densities at 0
degrees.
On Figure 4.4.b, a similar phenomenon is observed although the error does not tend
towards zero as CFL tends towards 1. This can be explained by the fact that the highest
frequency of elements is determined by the (faster) longitudinal wave and not the
(slower) shear wave (see equation 2·67). For CFL=1, the time stepping frequency does
not match the highest frequency for a shear wave. This explains why the error is not
zero for shear waves.
94
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
Let us consider a hypothetical critical time step for a shear wave at 0 degrees, ΔtS0,
which would correspond to the maximum frequency for this wave, ωmax,S0, such as:
L0 2 cL
Δt s0 = ------ = ------------------- = ----- Δt cr (4·4)
cS ω max, L0 cS
Δts0 is larger than Δtcr, so the model could not be run as it would be unstable. We define
MS0 as the ratio between ΔtS0 and Δtcr:
Δt cr
M S0 = ---------- (4·5)
Δt S0
L0 2
Δt L0 = ------ = ------------------- = Δt cr (4·7)
cL ω max, L0
Δt cr
and M L0 = ---------- = 1 (4·8)
Δt L0
We also define:
Plotting the curves of Figure 4.4.a and 4.4.b against CFLX leads to a good agreement
as shown on Figure 4.5.
This confirms that the amplitude of the velocity errors on both waves are subject to the
same mechanism. A general analytical function of CFLX fitting these curves is:
2
E c ( N, CFL X ) = E c ( N ) ( CFLX → 0 ) ⋅ ( 1 – CFL X ) (4·10)
where E c ( N ) ( CFLX → 0 ) is the value of the velocity error when CFLX tends to zero. This
value is dependent on the mesh density and can rightly be considered as the worst
95
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
5 N=6
3
8
2 9
10
11
12
1
15
20
30
0
0 0.2 0.4 0.6 0.8 1
CFLX
Shear wave Longitudinal wave c=K(1-CFLX2) curve shape
Figure 4.5. Velocity errors against CFLX for various mesh densities at 0 degrees.
possible error for a given mesh density. In order to establish the dependency, the
velocity error for both the longitudinal and shear waves at 0 degrees are plotted against
the mesh density N in Figure 4.6. For these cases, CFL is taken equal to 0.025.
Following equation 4·10, for this value of CFL, the velocity error is less than 0.1%
away from E c ( N ) ( CFLX → 0 ) .
Figure 4.6 highlights the strong influence of the mesh density on the velocity error. The
shape of the curve is hyperbolic and shows how the error varies from a value of about
5% for N=6 to a value of less than 0.2% for N=30.
4
Velocity error (%)
0 Shear wave
5 10 15 20 25 30 Longitudinal wave
Mesh density N Fitted curve
Figure 4.6. Velocity error against mesh density for shear and longitudinal waves at 0 degree with a
CFL of 0.025.
96
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
As mentioned previously, the mesh density value has a great impact on the
computational model size and therefore the amount of memory required to solve a
model. For example, in 3D, a model with N=30 would require 125 times more memory
than a model with N=6. It would also take much longer to solve as the number of time
steps would be multiplied by 5, and each time step would take longer to solve. This
illustrates the importance of carefully choosing N to match the accuracy required.
In order to facilitate the choice of N, the error variation can be fitted to the following
analytical expression:
180 180
E c ( N ) ( CFLX → 0 ) = --------
2
- = -----------
2
(4·11)
N ⎛ --λ-⎞
⎝ L⎠
This expression gives the value of the velocity error within 0.15% in a conservative
manner based on the results in this work.
So far we have only considered the error along the element side (at 0 degrees). In a
mesh of square elements such as the one in Figure 4.3, the mesh density is not constant
for all angles of propagation of a wave. In particular, at element level, L0=L90 but
L45= 2 L0. In order to study this, the longitudinal and shear velocities were monitored
from 0 to 90 degrees for various mesh densities N defined at zero degrees. As seen
above, the error varies quite strongly for CFL close to 1 whereas it is almost constant
at low CFL. For this series of models, a low CFL value of 0.05 is taken in order to avoid
a variation of the error due to this parameter.
Figure 4.7 shows the variation of the error against the angle for both wave types and a
mesh density varying from 6 to 30. The expected symmetry exists at 45 degrees. The
curves for both wave types and all mesh densities follow a similar pattern. The
minimum error is located at 0 degrees and increases as the angle increases towards 45
degrees where it reaches its maximum. There is a good match between shear and
longitudinal errors which indicates that the same mechanism determines the error
variation against the angle.
97
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
90
c) d)
Figure 4.7. Variation of the longitudinal (a and c) and shear (b and d) velocity error against the angle of
incidence for various values of mesh density plotted in polar (a and b) and linear (c and d) plots.
On Figure 4.8, it is very interesting to note that the actual error at 0, 45 and 90 degrees
matches the value given closely by the analytical expression of equation 4·11 when it
is calculated using L as L0, L45 and L90.
12
10
Error on velocity (%)
4
Shear wave at 0 degree
2
Longitudinal wave at 0 degree
Shear wave at 45 degrees
0
0 5 10 15 20 25 30 Longitudinal wave at 45 degrees
Mesh density N (λ/L) Fitted curve
Figure 4.8. Velocity error against mesh density for shear and longitudinal waves at 0 and 45 degree.
The agreement is relatively good for all values of N. The maximum difference between
the predicted values and the actual value is 0.6% for N=4.24 (mesh of density 6 at 45
degrees). The analytical expression used for the prediction gives a conservative
estimate for all cases apart from extremely low mesh densities (N<6). This should not
98
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
be a cause for concern. Such low mesh densities are only rarely used as the wave
phenomenon is likely to be incorrectly represented.
The agreement between predicted and measured values does not exist at other angles
than 0, 45 and 90 degrees if we were to use Lθ to calculate the mesh density N but it
can be precisely predicted using another expression. The variation of the error between
these values calculated with the analytical expression can be fitted to the following
expression:
π 2
E ( θ ) = E ( 0 ) + ( E ( 45 ) – E ( 0 ) ) ⎛ sin ⎛ 2θ ---------⎞ ⎞ (4·12)
⎝ ⎝ 180⎠ ⎠
with E(0) and E(45) the errors at 0 and 45 degrees calculated using equation 4·11 and
θ in degrees. Therefore, we are able to predict the error for any direction of the mesh
by using the fitted analytical expressions.
As discussed previously, the value of CFLX influences the value of the error. At 45
degrees, the value of ΔtL45 and ΔtS45 is higher than ΔtL0 and ΔtS0 and is the highest in
the model. This restricts the range of CFLX as illustrated on Figure 4.9 and implies that
the maximum error due to the CFL occurs at 45 degrees.
12
N=6
10
Error on longitunal velocity
at 45 degree (%)
8
7
6
8
9
4 10
11
12
2 15
20
30
0
0 0.2 0.4 0.6 0.8 1
CFLX
Shear wave Longitudinal wave c=K(1-CFLX2) curve shape
Figure 4.9. Velocity errors against CFLX for various mesh densities at 45 degrees.
Figure 4.10 presents the results of Figure 4.5 over a surface that shows the variation of
the velocity error against the scaled Courant number CFLX and mesh density N. This
99
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
surface can be used to predict the error on the velocity in the principal directions of the
mesh (0, 45, 90 degrees) for any mesh density and stability factor.
2
( 1 – CFLX )-
E c = 180 -------------------------------
2
(4·13)
N
5
4
Velocity error (%)
3
2
1
0
1
0.8 10
0.6 15
0.4 20
Scaled Courant number CFLX 0.2 25 Mesh density N
0 30
Figure 4.10. Velocity error against the scaled Courant number CFLX and mesh density N
For a given model, as shear waves have shorter wavelength than longitudinal waves,
the mesh density is lower for shear waves than longitudinal waves. Also, as shown in
this section, CFLX will be lower than for longitudinal waves. As the mesh density is
lowest at 45 degrees, the overall maximum error will be indicated by equation 4·13 for
a shear wave travelling at 45 degrees with the lowest CFLX.
The models studied in this part are similar to those used in the previous part apart from
the fact that the elements are “stretched” in one direction to become rectangular. The
element size in the vertical direction y is kept the same for all models so that the mesh
density N in this direction is 15. The element size in the horizontal direction x is varied
so that the mesh density N varies from 7.5 to 30. The ratio between the element side
R=Lx/Ly is hence varied from 0.5 to 2 as illustrated on Figure 4.11.a. The error on the
100
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
velocity is measured over a range of angles from 0 to 90 degrees and plotted on Figure
4.11.b,c,d,e. Small circles are plotted on the linear plots. These represent the amplitude
of the error predicted using the expression given in equation 4·13.
90 90
b) c)
60 60
4 R= 4
0 0 0 0
5 5
Longitudinal velocity error (%)
d) e)
Shear velocity error (%)
4 R 4 R
3 3
2 2
1 1
0 0
0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90
θ (deg) θ (deg)
Figure 4.11. a) Shape of the different rectangular elements used in the mesh; Variation of the
longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence for various R
plotted in a polar (b and c) and linear (d and e) fashion. The coloured circles indicate the error
prediction along the element side and diagonal.
It can be noted that, as seen in the case of square elements, longitudinal and shear
waves react in a similar way to the element distortion. At 90 degrees, the mesh density
is constant for all models. At this angle, the error is almost the same for all models and
matches the prediction of equation 4.13 very closely, within 0.08%. As the angle θ is
decreased, it can be seen that the error increases smoothly - in a similar way to that
observed for square elements - until it reaches its maximum. As would be expected, the
maxima are reached at different angles depending on R, but it is interesting to note that
the angle at which the maximum is reached matches the angle of the element diagonal
only for R equals 1. In other cases, an angle difference exists and increases non linearly
101
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
as R moves away from 1. Interestingly, despite the drift of the position of the maximum
error, the actual value of the maximum very closely matches the predicted error given
using the diagonal of the element to evaluate the mesh density. The maximum
predicted error is within 0.2% of the actual error. Following the maximum value, the
error decreases smoothly to reach a minimum at 0 degrees. As expected, at this angle,
the error varies as the mesh density differs. The error prediction at this angle agrees
well with the actual result. It is therefore clear that the analytical predictions match the
results obtained at 0 and 90 degrees extremely well.
The study of square elements showed that the prediction on the diagonal of a square
element could be used. For the other cases, it is also noticed that the analytical
prediction of the velocity error on the diagonal closely matches the actual maximum,
despite the fact that the diagonal angle does not match the angle at which the maximum
occurs. Over the cases studied, the difference between predicted errors and actual
maximum errors is 0.2%. It is therefore justified to apply the model devised for square
elements to rectangular elements.
In this case, the model used for the square elements study is sheared so that all elements
are deformed from a square to a rhombus. The models studied have a mesh density of
15 along the sides of the elements. A low CFL of 0.05 is used. The shearing angle γ is
varied from 0 to 45 degrees by 15 degrees increments. The error on the velocity is
monitored from -90 to +90 degrees, as in this case there is no symmetry in the
horizontal direction. Results for longitudinal and shear waves are plotted on Figure 4·7.
Small circles are plotted on the linear plots. These represent the amplitude of the error
predicted using the expression given in equation 4·13.
All cases for both longitudinal and shear waves exhibit two maxima and two minima.
Contrary to the case with a mesh of rectangular elements, the minima occur at the same
angle as the diagonals. This corresponds with the minima and maxima for Lθ. As
shown on Figure 4.12.b and c, the values of the maxima and minima are well predicted
in a conservative manner. It can be noted that the difference between the prediction and
the measured error increases for maxima on the longest diagonal and decreases for
maxima on the shortest diagonal. This is more pronounced for shear waves than
102
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
γ γ
γ
b) 3 90 c) 3 90
60 60
θ (deg) θ (deg)
2 γ=0 deg 2
γ=15 deg
Longitudinal velocity error (%)
30 γ=30 deg 30
γ=45 deg
0 0 0 0
-30 -30
-60 -60
-90 -90
3 3
Longitudinal velocity error (%)
d) e)
Shear velocity error (%)
2 γ 2 γ
1 1
0 0
-80 -60 -40 -20 0 20 40 60 80 -80 -60 -40 -20 0 20 40 60 80
θ (deg) θ (deg)
Figure 4.12. a) Shape of the different rhombic elements used in the mesh; Variation of the longitudinal
(b and d) and shear (c and e) velocity error against the angle of incidence for various shearing angle γ
plotted in a polar (b and c) and linear (d and e) fashion. The coloured circles indicate the error
prediction along the element side and diagonal.
103
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
The mesh in this case is sheared and stretched. This means that the elements are
transformed from a square to a parallelogram shape. The mesh density is kept constant
in both the horizontal and the vertical directions. Contrary to the rhombus case, the
element sides are not of equal length.
The models studied have a mesh density of 15 at 0 and 90 degrees. A low CFL of 0.05
is used. The shearing angle γ is varied from 0 to 45 degrees by 15 degree increments.
The error on the velocity is monitored from -90 to +90 degrees as in this case there is
no symmetry in the horizontal direction. Results for longitudinal and shear waves are
plotted on Figure 4.13. Small circles are plotted on the linear plots. These represent the
amplitude of the error predicted using the expression given in equation 4·13.
The general evolution of the error agrees with the rhombus case. The error is the same
for all cases at 0 degrees and agrees well with the prediction. Contrary to the rhombus
case, but as for the rectangle case, there is a shift between the actual position of the
maximum and the diagonal where it would be expected. This confirms that the shift is
likely to be due to the difference in length of the two pairs of opposite edges. Despite
this shift, the actual value of the error and the predicted value on the diagonal agree
closely. As seen before, the maximum value for the shear wave error is increasingly
overestimated as the shearing angle increases. This is interpreted as a stiffening of the
element along its longest diagonal as it is deformed. The longitudinal wave does not
seem to suffer from this effect as strongly as the shear wave.
104
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
γ γ
γ
4 90 c) 4 90
b)
60 60
γ=0 deg
3 θ (deg) 3 θ (deg)
γ=15 deg
γ=30 deg
Longitudinal velocity error (%)
2 30 γ=45 deg 2 30
0 0 0 0
-30 -30
-60 -60
-90 4 -90
4
Longitudinal velocity error (%)
d) e)
Shear velocity error (%)
3 3
2 γ 2 γ
1 1
Figure 4.13. a) Shape of the different parallelogramatic elements used in the mesh; Variation of the
longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence for various
shearing angle γ plotted in a polar (b and c) and linear (d and e) fashion. The coloured circles indicate
the error prediction along the element side and diagonal.
4.2.2.5 Conclusion
The study of square element meshes has shown that the error along the element sides
and diagonals can be predicted by the expression:
2
( 1 – CFLX )-
E c = 180 -------------------------------
2
(4·14)
N
where CFLX is the scaled Courant number and N is the mesh density.
105
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
Studies of deformed meshes made of rectangles, rhombi and parallelograms show that
this expression also enables the obtention of a correct conservative approximation of
the velocity error along the element sides and diagonals. This study looks at element
deformation where square elements are stretched with element side ratios up to 2 and
sheared up to 45 degrees. The main consequence of this point is that it indicates that as
square elements are deformed in such a way, they do not lose their precision
dramatically. Therefore, it is acceptable to apply such deformation to square elements
in practical models.
In this part, the influence of CFL and N on wave propagation in a regular mesh made
of linear equilateral triangular elements is investigated. A typical mesh is shown in
Figure 4.14.
L0
L90
L30
Lθ
Figure 4.14. Schematic defining L0, L90, L30 and Lθ in a mesh of equilateral-triangular elements.
The critical time step for this case is determined by the value L0 and not the element
side length L30 or L90 as in the square element case, as explained in Section 2.3.1. The
Courant number CFL is varied from 1 to 0.05 and the mesh density N ranges from 6 to
30 along the element side.
106
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
A first series of models is run. The error on both shear and longitudinal velocities is
monitored over a range of angles from 0 to 90 degrees for a range of mesh densities
using a low CFL of 0.05. Results are presented in Figure 4.15.
c)
N N
3 3
2 2
1 1
0 0
0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90
θ (deg) θ (deg)
Figure 4.15. Variation of the longitudinal (a and c) and shear (b and d) velocity error against the angle
of incidence for various mesh densities plotted in a linear (a and b) and polar (c and d) fashion.
As expected there are symmetries at 0, 30, 60 and 90 degrees in both cases. The
longitudinal error is maximum at 30 and 90 degrees in the direction aligned with the
mesh where Lθ is maximum, and the error is minimum at 0 and 60 degrees where Lθ
is minimum. Between these points, the longitudinal error varies smoothly in a similar
fashion to the phenomena observed for quadrilateral elements. On the other hand, the
shear wave velocity error reaches its maximum at 0 and 60 degrees where Lθ reaches
its minimum and is unexpectedly close to zero at 30 and 90 degrees in the direction
aligned with the mesh. Between these extremes, the variation is smooth but, unlike in
previous cases, this variation cannot be described by a squared sine relation.
107
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
It is interesting to note that the value of maxima (at 30 and 90 degrees) and minima (at
0 and 60 degrees) of the error for the longitudinal wave and the maxima (at 0 and 60
degrees) of the error for the shear wave can be predicted accurately using the equation
defined previously for quadrilateral elements:
180- 180
E c = --------
2
= ----------2- (4·15)
N ⎛ --λ-⎞
⎝ L⎠
with L the element size defined as in Figure 4.14. This is clearly shown in Figure 4.16.
This expression gives a general conservative value of the maximum error for a regular
mesh of linear equilateral triangular elements. Given that the minimum shear error is
zero at 30 and 90 degrees, all minima and maxima can be precisely predicted.
4
Velocity error (%)
1
Shear wave at 0 degree
Longitudinal wave at 0 degree
0
5 10 15 20 25 30 35 Longitudinal wave at 30 degrees
Mesh density N Fitted curve
Figure 4.16. Velocity error against mesh density for shear and longitudinal waves at 0 and 30 degrees.
For longitudinal waves, it can be noted that the variation of the longitudinal wave error
is precisely described by the following fitted expression:
π 2
E ( θ ) = E ( 0 ) + ( E ( 30 ) – E ( 0 ) ) ⎛ sin ⎛ 3θ ---------⎞ ⎞ (4·16)
⎝ ⎝ 180⎠ ⎠
with E(0) and E(30) the error at 0 and 30 degrees and θ the angle in degrees.
Figure 4.17 shows the influence of CFLX on the value of the error on the longitudinal
velocity at 0 and 90 degrees and the shear velocity at 0 degrees. The error tends towards
zero as CFLX tends to 1. This confirms that the stability limit is not defined by the
smallest gap between nodes but by the transit time of a dilatational wave through the
108
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
smallest element in the model as described in Section 2.3.1. Once again, it can be seen
that the errors can be fitted to the analytical expression:
2
E c ( N, CFL X ) = E c ( N ) ( CFLX → 0 ) ⋅ ( 1 – CFL X ) (4·17)
where E c ( N ) ( CFLX → 0 ) is value of the velocity error when CFLX tends to zero.
a) 5
4
N=6
at 0 degree (%)
Velocity error
3
7
2 8
9
10
1 11
12
15
20 30
0
0 0.2 0.4 0.6 0.8 1
CFLX
b) 5
N=6
4
7
at 30 degree (%)
Velocity error
3
8
2 9
10
11
12
1
15
20
30
0
0 0.2 0.4 0.6 0.8 1
CFLX
Shear wave Longitudinal wave c=K(1-CFLX2) curve shape
Figure 4.17. Velocity errors against CFLX for various mesh densities at a) 0 and b) 30 degrees.
It is therefore possible to predict the maximum longitudinal error along the element
side, the minimum longitudinal error along the altitude of the element and the
maximum shear error along the altitudes of the triangular element for any mesh density
N and relative Courant number CFLX by using the following expression:
2
( 1 – CFLX )
E c = 180 -------------------------------
2
- (4·18)
N
109
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
This section looks into the influence of having a mesh of isosceles triangular elements
rather than equilateral ones. It is essential to be able to evaluate the impact of
deforming elements from the equilateral form as automatic meshing algorithms used
by standard FE package will commonly contain non-equilateral triangular elements. In
this study, we used a mesh similar to the one used previously but the mesh density is
kept the same in the y (vertical) direction whereas it is varied in the x (horizontal)
direction. This leads to a mesh constituted of same-size isosceles triangles in a similar
way to the equilateral case considered previously. The 7 isosceles triangles shapes used
for this study are shown in Figure 4.18.a. For each triangle, one angle is varied from
30 to 90 degrees while the two other angles are equal. It can be seen that the mesh
density will be constant in the vertical direction and is taken as 10 elements per
wavelength whereas the density in the horizontal direction is varied from 5.35 to 20.
The mesh is stretched in the horizontal direction and the mesh density is kept the same
in the vertical direction. The equal angles are varied from 75 degrees to 45 degrees
leading to a variation of the other angle φ from 30 to 90 degrees. Results of the
variation of the velocity error against the angle are shown on Figure 4.18.
Looking at the longitudinal wave error, it can be seen that the constant mesh density at
90 degrees leads to the error being between 1.35% and 1.70% and relatively close to
the estimate of 1.8% predicted by equation 4·18. Despite the constant mesh density in
this direction, one can see that the error drops systematically as φ is increased at this
angle. At 0 degrees, a much greater difference exists between cases. This is justified
by the fact that the mesh density varies strongly in this direction. Once again, the
prediction agrees closely with the actual values. Between these values, a cycle similar
to that for equilateral linear elements is observed. The maximum and minimum occurs
close to - but not exactly at - the location of the local maximum and minimum of Lθ.
These occur when Lθ is aligned with the mesh and when Lθ is perpendicular to the
opposite face of an element. Although the angular position of the maximum is slightly
incorrect, the predicted value of the error proves to be close to the actual minimum and
maximum. It is therefore possible for the modeller to predict the longitudinal velocity
error quite precisely using equation 4·18.
110
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
φ φ φ φ φ φ φ
Longitudinal velocity error (%)
90 6 90
b) c)
60 60
d) e) 6
6
0 -4
0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90
θ (deg) θ (deg)
Figure 4.18. a) Shape of the different isosceles-triangular elements used in the mesh; Variation of the
longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence for various
values of φ plotted in a polar (b and c) and linear (d and e) fashion. The coloured circles indicate the
error prediction along the element side and diagonal.
For shear waves, a maximum (local for cases 5 to 7) occurs at 0 degrees. The value of
this maximum matches the prediction closely as represented on the Figure 4.18. For all
cases, another maximum (local for cases 1 to 3) occurs at an angle close to the direction
of the median of the upper triangle. The value of this maximum is also closely
predicted. Between these maximum values, the error varies smoothly toward minimum
values located at angles where Lθ is aligned with the mesh. For case 4 where φ is 60
degrees the triangle is equilateral, this minimum occurs at θ equals 30 degrees and is
close to zero. As φ is increased, the value of the minimum increases despite the fact
that Lθ is reduced. When φ is decreased from 60 degrees, the minimum error becomes
negative, meaning that the actual velocity is higher than the theoretical one. This
phenomenon occurs despite the fact that Lθ is increased and indicates a stiffening
effect along the longest edge of the element. The same effect occurs at 90 degrees when
minima (local for cases 1 to 3) occur. At this angle, a decrease of φ leads to an increase
111
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
These results strongly highlight that using two triangular linear elements to create a
“composite quadrilateral element” leads to results which will not match those of a real
quadrilateral linear element and should therefore be avoided.
This part looks into the influence of having a mesh of scalene triangular elements. The
mesh used for the equilateral case is sheared in a similar way to the parallelogram study
for quadrilateral elements. The mesh density is kept the same in the y (vertical) and x
(horizontal) direction at 10 and 11.55 elements per wavelength respectively while the
element is sheared in the y direction. The shearing angle varies from 0 to 45 degrees.
The 4 scalene triangular shapes used for this study are shown in Figure 4.19.a. Case 1
is the equilateral case while case 3 uses right (but not isosceles) triangles. Results of
the variation of the velocity error against the angle are shown on Figure 4.19.
For longitudinal waves, the variation of the error for cases 1, 2 and 3 is similar to what
was experienced in previous studies. For case 4, an unexpected minimum occurs at 22
degrees and an unexpected maximum occurs at -63 degrees. These could be explained
by the cumulation of the effects of mesh density variation and a stiffening and
softening of the element along the longest edge and the altitude from the longest edge
respectively. The effect explains why the overall maximum for case 4 is strongly
underestimated and the overall minimum slightly overestimated. A similar milder
effect is occurring for case 3 where the difference between the predicted and actual
maximum is 0.75%.
For shear waves, the variation of the error highlights the same issue observed with
isosceles triangles with large angle opposite long edges as the error becomes negative.
112
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
Case 4 is the most critical case studied for triangle elements so far as the ratio of the
longest to the shortest edge is 1.62 and the long edge is opposite to an angle of 107
degrees. The stiffening effect is the most pronounced as the error reaches a maximum
negative value of 16.4% which is more than 10 times the maximum predicted error for
this case. This emphasizes that using highly deformed triangular linear elements for
wave propagation should be avoided.
φ
φ φ
φ
90 2 90
b) c)
60 60
2 γ=0 deg
θ (deg) 0 θ (deg)
γ=15 deg
γ=30 deg
Longitudinal velocity error (%)
30 γ=45 deg -2 30
1
Shear velocity error (%)
-4
0 0 -6 0
-30 -30
-60 -60
-90 -90
5 2
Longitudinal velocity error (%)
d) e)
0
Shear velocity error (%)
4 -2
-4
3 -6 γ
-8
2 γ -10
-12
1 -14
-16
0 -18
-80 -60 -40 -20 0 20 40 60 80 -80 -60 -40 -20 0 20 40 60 80
θ (deg) θ (deg)
Figure 4.19. a) Shape of the different scalene-triangular elements used in the mesh; Variation of the
longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence for various
values of γ plotted in a polar (b and c) and linear (d and e) fashion. The coloured circles indicate the
error prediction along the element side and diagonal.
113
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
4.2.3.4 Conclusion
The study of equilateral triangular element meshes has shown that the error along the
element sides and altitudes could be accurately predicted. The influence of the Courant
number on the error is confirmed to be similar to the one observed with quadrilateral
elements.
Studies of deformed meshes made of isosceles and scalene triangles shows that it is
possible to use predicted errors to conservatively evaluate the maximum error for
longitudinal waves. The situation with shear waves is rather different as it is observed
that an apparent stiffening of the deformed elements leads to a large increase in the
propagation velocity leading to a large velocity error that cannot be predicted in the
same way as was done in previous cases. This effect is a noticeable concern regarding
the use of linear triangular elements with automatic meshing algorithms.
In this part, we investigate the influence of CFL and N on wave propagation in a regular
mesh made of quadratic equilateral-triangular elements as shown in Figure 4.20. The
mesh density N and CFL is expressed in terms of the node grid and not the actual
element size as it will be shown that this is what actually matters with this element type.
As for linear elements, the critical time step is determined by the value L0 and not the
element side length L30 or L90. The Courant number CFL is varied from 0.9 to 0.05
and the mesh density N ranges from 6 to 30 along the element side.
114
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
L0
L90
L30
θ
Lθ
Figure 4.20. Schematic defining L0, L90, L30 and Lθ in a mesh of quadratic equilateral-triangular
elements.
A first series of models is run. The error on both shear and longitudinal velocity is
monitored over a range of angle from 0 to 90 degrees for a range of mesh densities
using a low CFL of 0.05. Results are presented in Figure 4.21.
90 90
a) b)
Longitudinal velocity error (%)
60 N= 60
4 6 4
Shear velocity error (%)
θ (deg) 7 θ (deg)
8
30 9 30
2 10 2
11
12
15
0 20 0
0 30 0
5 d) 5
c)
Longitudinal velocity error (%)
N N
4 4
Shear velocity error (%)
3 3
2 2
1 1
0 0
0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90
θ (deg) θ (deg)
Figure 4.21. Variation of the longitudinal (a and c) and shear (b and d) velocity error against the angle
of incidence for various mesh densities plotted in a polar (a and b) and linear (c and d) fashion.
The results for longitudinal waves indicate that for any mesh density higher than 7 the
error remains close to constant at any angle. This is an interesting feature as it means
that the wave propagates in all directions at a constant velocity as one would expect in
an isotropic material. For a mesh density of 6 nodes per wavelength, the number of
elements per wavelength is 3. The inconsistency of the results in this particular case
115
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
indicates that this is not sufficient to correctly represent longitudinal wave propagation.
As the mesh density is increased, the strength of these inconsistencies decreases
indicating that the wave propagation is correctly represented. This point raises the
issue of the minimum number of elements necessary to represent the wave propagation
correctly and will be discussed further.
For shear waves, the pattern observed for all mesh densities is similar to the one
observed for linear quadrilateral elements as the maxima occur at 0 and 60 degrees and
the minima at 30 and 90 degrees with a squared sine variation between them.
The improvement in the consistency of the velocity for both longitudinal and shear
waves against the angle of propagation is an interesting feature. The fact that the mesh
density needs to be at least equal to 7 to correctly represent the wave propagation is of
limited importance. It indicates that the use of quadratic elements increases the
minimum computational needs, but this only affects a limited range of models. The
influence of the mesh density on the error is represented in Figure 4.22.
3
Velocity error (%)
0
5 10 15 20 25 30 35
Mesh density N
Shear wave at 0 degree
Longitudinal wave at 30 degrees
2
Fitted curve 180/N
Longitudinal wave at 0 degree
2
Fitted curve 180/(Ncos(π/6))
Shear wave at 30 degrees
2
Fitted curve 180/(N/cos(π/6))
Figure 4.22. Velocity error against mesh density for shear and longitudinal waves at 0 and 30 degrees.
116
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
As the longitudinal velocity is constant against the angle, the error is plotted against
the mesh density defined by the node spacing on the edge of the elements. Plotting the
error for both longitudinal and shear waves at 0 and 30 degrees, it can be seen that the
errors for shear waves at 0 degrees and for longitudinal waves at 30 degrees agree well
with the prediction given for linear elements:
180- 180
E c = -------- = -------------2- (4·19)
N
2 λ⎞
⎛ -----
⎝ L θ⎠
The plot also indicates the trend followed by the longitudinal error at 0 degrees as:
180 180
E c = --------------------------2- = ----------------------------------2 (4·20)
π λ⎞ 2 π
N cos ⎛ ---⎞ ⎛ ----- cos ⎛ ---⎞
2
⎝ 6⎠ ⎝ L θ⎠ ⎝ 6⎠
π 2 π 2
180 cos ⎛ ---⎞ 180 cos ⎛ ---⎞
⎝ 6⎠ ⎝ 6⎠
E c = ----------------------------- = ----------------------------- (4·21)
N
2 λ ⎞2
⎛ -----
⎝ L θ⎠
This validates that the maximum errors for quadratic and linear elements vary in the
same way and with the same amplitude as the mesh density varies but that the
minimum error follows different rules.
The next point to look into is the influence of the Courant number CFL on the error. It
is varied from 0.9 to 0.05. The scaled Courant number CFLX varies from 0.78 to 0.043
for a longitudinal wave at 30 degrees and from 0.45 to 0.025 for a shear wave at 0
degrees. The errors for longitudinal waves at 90 degrees and shear waves at 0 degrees
are plotted against CFLX over the range of mesh density studied. Results are shown on
Figure 4.23. The variation of the error against CFLX is similar to the one experienced
for linear elements.
117
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
2
( 1 – CFLX )-
E c = 180 -------------------------------
2
(4·22)
N
a) 5
N=6
4
Longitudinal velocity error
7
at 30 degree (%)
3
8
2 9
10
11
12
1
15
20
30
0
0 0.2 0.4 0.6 0.8 1
CFLX
b) 5
4 N=6
Shear velocity error
at 0 degree (%)
3
7
8
2
9
10
11
1 12
15
20
0 30
0 0.2 0.4 0.6 0.8 1
CFLX
Shear wave Longitudinal wave c=K(1-CFLX2) curve shape
Figure 4.23. Velocity errors against CFLX for various mesh densities at a) 0 and b) 30 degrees.
118
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
This study is similar to that done for linear isosceles triangular elements. The same
node grid is used but is filled with quadratic rather than linear elements. This means
that the mesh density is the same and the element size is doubled. The mesh density is
kept at 10 nodes per wavelength in the y (vertical) direction whereas it is varied in the
x (horizontal) direction from 5.35 to 20. The mesh is hence constituted of same-size
isosceles triangles. The 7 isosceles triangle shapes used for this study are shown in
Figure 4.24.a. For each triangle, one angle is varied from 30 to 90 degrees while the
two other angles are kept equal. The equal angles are varied from 75 degrees to 45
degrees leading to a variation of the other angle φ from 30 to 90 degrees. Results of the
variation of the velocity error against the angle are shown on Figure 4.24.
φ φ φ φ φ φ φ
8 90
Error on longitudinal velocity (%)
b) 8 c)
60 60
Error on shear velocity (%)
φ=
6 θ (deg) 30 6 θ (deg)
40
4 30 50 4 30
60
70
2 80 2
90
0 0 0 0
9 9
e)
Error on longitudinal velocity (%)
d) 8 8
Error on shear velocity (%)
7 7
6 φ 6 φ
5 5
4 4
3 3
2 2
1 1
0 0
0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90
θ (deg) θ (deg)
Figure 4.24. a) Shape of the different quadratic isosceles-triangular elements used in the mesh;
Variation of the longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence
for various values of φ plotted in a polar (b and c) and linear (d and e) fashion. The coloured circles
indicate the error prediction along the element side and diagonal.
119
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
For longitudinal waves, as seen previously, the error for the equilateral case is constant.
For other cases, the maximum and minimum are located at 0 and 90 degrees and a
smooth variation occurs between these 2 points. As the angle φ is reduced, the
maximum error is located at 0 degrees and the minimum at 90 degrees. For φ greater
than 60 degrees, the maximum is at 90 degrees and the minimum at 0 degrees. The
predictions given by equation 4·18 agree well with the actual minima and maxima and
are conservative for the maxima.
For shear waves, the variation follows a similar pattern to the one observed for bulk
waves although the undulation occurring for the equilateral case is visible on the other
cases but does not seem to influence the overall results. The prediction of the maxima
is accurate for all cases apart for case 1 for which the prediction is 0.9% below the
actual value. It is also noted that the prediction is not conservative for case 2 being
0.1% below the correct value. For all other cases, the estimate is conservative and
within 0.2%. It is therefore correct to use the estimates as outlined previously for
equilateral triangles.
It is very interesting to note the significant improvement observed for shear waves
compared to linear elements. The shear error does not deteriorate dramatically as the
element is deformed. Quadratic elements can be used without any fear of a small
deformed element giving a larger than expected error for shear waves. This allows the
modeller to use automatic meshing algorithms more confidently.
This part looks into the influence of having a mesh of quadratic scalene triangular
elements. The node mesh is the same as used previously for the linear case. The mesh
density is kept the same in the y (vertical) and x (horizontal) direction at 10 and 11.55
elements per wavelength respectively while the element is sheared in the y direction.
The shearing angle varies from 0 to 45 degrees. The 4 scalene triangular shapes used
for this study are shown in Figure 4.25.a. Case 1 is the equilateral case while case 3
uses right (but not isosceles) triangles. Results of the variation of the velocity error
against the angle are shown on Figure 4.25.
120
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
φ
φ φ
φ
b) 4 90 c) 4 90
60 60
γ=0 deg
3 θ (deg) 3 θ (deg)
γ=15 deg
γ=30 deg
Error on longitudinal velocity (%)
2 30 γ=45 deg 2 30
0 0 0 0
-30 -30
-60 -60
-90 -90
5 5
e)
Error on longitudinal velocity (%)
d)
Error on shear velocity (%)
4 4
3 3
γ γ
2 2
1 1
0 0
-80 -60 -40 -20 0 20 40 60 80 -80 -60 -40 -20 0 20 40 60 80
θ (deg) θ (deg)
Figure 4.25. a) Shape of the different quadratic scalene-triangular elements used in the mesh; Variation
of the longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence for various
values of γ plotted in a polar (b and c) and linear (d and e) fashion. The coloured circles indicate the
error prediction along the element side and diagonal.
For longitudinal waves and all cases apart from case 1, the error varies smoothly and
describes a full sine like cycle. Maxima are reached between 0 and 90 degrees and a
minimum between -90 and 0 degrees. The prediction for the maxima is calculated
along the longest edge of the element. The estimation is relatively accurate although,
as γ is increased, the margin grows larger. This is in line with the relative error that was
observed at 90 degrees for isosceles triangles.
121
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
For shear waves, a similar pattern is observed although the undulation occurring for
case 1 (equilateral triangle), exists for all cases and modifies slightly the overall shape
of the curves. The predictions for the maximum error are calculated along the longest
element side and agree well with the actual results in a conservative way. Again, using
modified quadratic elements instead of linear ones dramatically improves the precision
of the shear wave propagation velocity.
4.2.4.4 Conclusion
The study of equilateral triangular meshes made of modified quadratic elements has
shown that the longitudinal velocity error is very close to constant for N superior to 7.
For shear waves, a velocity error variation exists but is limited. This variation is closer
to the one observed with square elements than with equilateral triangular elements. The
maximum error for both wave types is identified and an expression can be used to
predict its amplitude. The influence of the Courant number on the error is confirmed
to be similar to the one observed previously in this chapter.
Studies of deformed meshes made of isosceles and scalene triangles show that the error
predictions derived from equation 4.26 are quite accurate although not always
conservative in cases where deformation is high. The most important point of this
study is that, contrary to linear triangular elements, no dramatic degradation of the
error is noted for any wave type and deformation level. Therefore, it can be said that
this element type can be used deformed without concern. It can be noted that the
modified quadratic elements react to deformation in a similar fashion as linear
quadrilateral elements. In the same way, the velocity error for a mesh of linear square
elements and modified quadratic equilateral triangle elements is relatively close. Based
on these two facts, it can be said that modified quadratic triangular elements behave in
a similar way to linear quadrilateral elements.
122
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
4.3.1 Introduction
In this section, the influence of the mesh density and angle of incidence are
investigated for a range of element types.
The geometry and properties of the models are the same as in the explicit cases. The
excitation is applied in the same way apart from the fact that it is applied as a real
harmonic signal of unit frequency. Models are solved in the frequency domain by the
implicit direct-solution steady state scheme of ABAQUS/Standard [14]. As the model
is solved in the frequency domain, no time increment is used and no Courant number
applies to these models. The monitored displacements are processed in the same way
in order to measure the propagation velocity.
All studies start with a mesh made of square or equilateral elements. Following this,
the meshes are deformed in order to understand how the deformation impacts the
velocity error. The mesh density is varied from 6 to 30 and is defined in the same way
as in Section 2.3.1.
The velocity errors for both shear and longitudinal velocities monitored over a range
of angles from 0 to 90 degrees for a range of mesh densities are plotted on Figure 4.26.
It can be noted that the general variation of the error is the same as the one observed
with reduced integration linear square elements in Section 4.2.2.1. The amplitude of
the error for the longitudinal wave superimposes almost perfectly with the results from
123
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
Section 4.2.2.1. For shear waves, the agreement is similar at 0 and 90 degrees but the
increase in error with reduced integration elements in the explicit case is stronger than
the full integration ones in the implicit case. This discrepancy is likely to be due to the
change in integration technique rather than the change in solving process.
90
Longitudinal velocity error (%)
N= b) 90
a)
60 6
c) d)
Figure 4.26. Variation of the longitudinal (a and c) and shear (b and d) velocity error against the angle
of incidence for various values of mesh density plotted in linear (a and b) and polar (c and d) plots.
The meshes studied in this part are “stretched” in one direction to create rectangular
elements. The element size in the vertical direction y is kept the same for all models so
that the mesh density N in this direction is 15. The element size in the horizontal
direction x is varied so that the mesh density N varies from 7.5 to 30. The ratio between
the element side R=Lx/Ly is hence varied from 0.5 to 2 as illustrated on Figure 4.27.a.
The error on the velocity is measured over a range of angles from 0 to 90 degrees and
plotted on Figure 4.27.b,c,d,e.
As seen with square elements, the longitudinal wave case matches the explicit case
almost perfectly and the shear wave case matches it at 0 and 90 degrees, but the
increase in error between these points is smaller than in the explicit case.
124
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
0 0 0 0
5 5
Longitudinal velocity error (%)
d) e)
3 3
2 2
1 1
0 0
0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90
θ (deg) θ (deg)
Figure 4.27. a) Shape of the different rectangular elements used in the mesh; Variation of the
longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence for various
values of R plotted in a polar (b and c) and linear (d and e) fashion.
Meshes in this case are sheared so elements are deformed from a square to a rhombus.
The models studied have a mesh density of 15 along the element sides. The shearing
angle γ is varied from 0 to 60 degrees by 15 degree increments. The error on the
velocity is monitored from -90 to +90 degrees and is plotted on Figure 4.28.
In these cases, the amplitude of the error is similar to the explicit case for low
deformation. It diverges quite strongly (in particular for shear waves) as the mesh is
deformed but stays smaller than the explicit case error.
125
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
γ γ
γ
b) 3 90 c)
90
60 60
θ (deg) θ (deg)
2 γ=0 deg 1
γ=15 deg
Longitudinal velocity error (%)
30 γ=30 deg 30
γ=45 deg
0 0 0 0
-30 -30
-60 -60
-90 -90
3 3
Longitudinal velocity error (%)
d) e)
Shear velocity error (%)
2 2
γ
γ
1 1
0 0
-80 -60 -40 -20 0 20 40 60 80 -80 -60 -40 -20 0 20 40 60 80
θ (deg) θ (deg)
Figure 4.28. a) Shape of the different rhombic elements used in the mesh; Variation of the longitudinal
(b and d) and shear (c and e) velocity error against the angle of incidence for various values of γ plotted
in a polar (b and c) and linear (d and e) fashion.
The mesh in this case is sheared and stretched with a mesh density of 15 at 0 and 90
degrees. The shearing angle γ is varied from 0 to 45 degrees by 15 degree increments
as shown on Figure 4.29.a. Results are plotted on Figure 4.29.
The findings for the parallelogram case are in line with those of the rhombus cases.
126
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
γ γ
γ
4 90 c)
90
b)
60 60
γ=0 deg
3 θ (deg) 1 θ (deg)
γ=15 deg
γ=30 deg
Longitudinal velocity error (%)
2 30 γ=45 deg 30
0 0 0 0
-30 -30
-60 -60
-90 4 -90
4
Longitudinal velocity error (%)
d) e)
Shear velocity error (%)
3 3
2 γ 2
γ
1 1
Figure 4.29. a) Shape of the different parallelogramatic elements used in the mesh; Variation of the
longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence for various
values of γ plotted in a polar (b and c) and linear (d and e) fashion.
4.3.2.5 Conclusion
The study of the implicit cases with fully integrated elements shows that results for
these cases are similar to the ones obtained for the reduced integration cases solved
with the explicit solver. Two differences are noticed: the amplitude of the error tends
to be lower for fully integrated elements than reduced integration elements with shear
waves for all element shapes; deformation of elements leads to a reduction in error
compared to the explicit cases. Although not optimum, the general rule to determine
127
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
the maximum velocity error established for explicit quadrilateral elements can be
adapted to implicit cases as:
180-
E cmax = ----------
2
(4·23)
N min
It is important to mention that, in this case, as in Section 4.2.4.1, the mesh density is
determined in terms of node spacing not element side size.
The velocity errors for both shear and longitudinal velocities monitored over a range
of angles from 0 to 90 degrees for a range of mesh densities are plotted on Figure 4.26.
90
Longitudinal velocity error (%)
N= b) 90
a)
60 6
Shear velocity error (%)
60
-1 7 -1
θ (deg) θ (deg)
8
9
30 10 30
-0.5 -0.5
11
12
15
0 20 0
0 0
30
0 0
Longitudinal velocity error (%)
c) d)
-0.2 -0.2
-0.4 -0.4
-0.6 -0.6
-0.8 -0.8
N N
-1.0 -1.0
-1.2 -1.2
0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90
θ (deg) θ (deg)
Figure 4.30. Variation of the longitudinal (a and c) and shear (b and d) velocity error against the angle
of incidence for various values of mesh density plotted in polar (a and b) and linear (c and d) plots.
128
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
Looking at the results for quadratic square elements, it can be seen that the absolute
value of the error is much lower than for linear elements as it is smaller than 0.5% for
a mesh density superior to 8. Unlike all cases seen previously, the velocity is higher
than the theoretical value which explains why a negative error is plotted. As seen in
previous cases, there is a variation in the error but here the smallest error occurs at 45
degrees where the mesh density is at its lowest. This indicates that the approach to
predict the velocity error that was used for linear and modified quadratic elements
cannot be applied to this case.
The meshes studied in this part are “stretched” in one direction to create rectangular
elements. The element size in the vertical direction y is kept the same for all models so
that the mesh density N in this direction is 15. The element size in the horizontal
direction x is varied so that the mesh density N varies from 7.5 to 30. The ratio between
the element side R=Lx/Ly is hence varied from 0.5 to 2 as illustrated on Figure 4.31.a.
The error on the velocity is measured over a range of angles from 0 to 90 degrees and
plotted on Figure 4.31.b,c,d,e.
Looking at the results for both longitudinal and shear waves, it can be seen that the
amplitude of the reflection is impacted by the deformation of the elements. The general
shape of the curves is similar to what was seen for linear quadrilateral elements
although the amplitude of the error stays in line with quadratic quadrilateral results. It
can be noted that the error is positive along the shortest edge of case 5. This indicates
a softening of the element in this direction and the error is higher than what would have
been expected at the mesh density in this direction.
129
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
0.1 0 0.2 0
0.1 0.1
Longitudinal velocity error (%)
d) e)
-0.2 -0.2
-0.4 R -0.4 R
-0.6 -0.6
0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90
θ (deg) θ (deg)
Figure 4.31. a) Shape of the different rectangular elements used in the mesh; Variation of the
longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence for various
values of R plotted in a polar (b and c) and linear (d and e) fashion.
Meshes in this case are sheared so that all elements are deformed from a square to a
rhombus. The models studied have a mesh density of 15 along the element sides. The
shearing angle γ is varied from 0 to 60 degrees by 15 degree increments as seen in
Figure 4.32.a. The error on the velocity is monitored from -90 to +90 degrees and is
plotted on Figure 4.32.
The deformation of the element from a square to a rhombus does not lead to a dramatic
deterioration of the error and is in line with amplitudes experienced before with
quadratic elements.
130
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
γ γ
γ
b) -0.4 90 c) -0.4
90
60 60
θ (deg) θ (deg)
γ=0 deg
-0.2 γ=15 deg -0.2
Longitudinal velocity error (%)
30 γ=30 deg 30
γ=45 deg
0.1 0 0.1 0
-30 -30
-60 -60
-90 -90
0.1 0.1
Longitudinal velocity error (%)
d) e)
Shear velocity error (%)
0 0
γ γ
-0.2 -0.2
-0.4 -0.4
-80 -60 -40 -20 0 20 40 60 80 -80 -60 -40 -20 0 20 40 60 80
θ (deg) θ (deg)
Figure 4.32. a) Shape of the different rhombic elements used in the mesh; Variation of the longitudinal
(b and d) and shear (c and e) velocity error against the angle of incidence for various values of γ plotted
in a polar (b and c) and linear (d and e) fashion.
The mesh in this case is sheared and stretched with a mesh density of 15 at 0 and 90
degrees. The shearing angle γ is varied from 0 to 45 degrees by 15 degree increments
as shown on Figure 4.33.a. Results are plotted on Figure 4.33.
Results are relatively similar to the rhombus results except that the velocity error has a
higher amplitude. This can be explained by the fact that the mesh density in this case
is more strongly modified than in the rhombus one.
131
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
γ γ
γ
b) -0.8 90 c) -0.8 90
60 60
γ=0 deg
θ (deg) θ (deg)
γ=15 deg
γ=30 deg
Longitudinal velocity error (%)
0.2 0 0.2 0
-30 -30
-60 -60
d) e)
Shear velocity error (%)
0 0
-0.2 γ -0.2 γ
-0.4 -0.4
-0.6 -0.6
Figure 4.33. a) Shape of the different parallelogramatic elements used in the mesh; Variation of the
longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence for various
values of γ plotted in a polar (b and c) and linear (d and e) fashion.
4.3.3.5 Conclusion
132
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
The error on both shear and longitudinal velocity is monitored over a range of angles
from 0 to 90 degrees for a range of mesh densities. Results are presented in Figure 4.34.
90
Longitudinal velocity error (%)
b) 90
a) N=
60 60
c)
3 3
2 2
1 1
0 0
0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90
θ (deg) θ (deg)
Figure 4.34. Variation of the longitudinal (a and c) and shear (b and d) velocity error against the angle
of incidence for various mesh densities plotted in a linear (a and b) and polar (c and d) fashion.
Results for this case superimpose almost perfectly with those of the explicit case with
linear equilateral triangle meshes. It is interesting to note that the implicit results match
the explicit ones obtained with a low CFLX. Findings from the explicit case can
therefore be applied to the implicit case, in which case, the CFLX in the former should
be replaced by zero in the latter.
The meshes studied in this part are “stretched” in one direction to have isosceles
triangle elements. The element size in the vertical direction y is kept the same for all
models so that the mesh density N in this direction is 10. The element size in the
horizontal direction x is varied so that the mesh density N varies from 5.35 to 20.
Results plotted on Figure 4.35 show an excellent agreement with the explicit case.
133
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
φ φ φ φ φ φ φ
Longitudinal velocity error (%)
90 6 90
b) c)
60 60
d) e) 6
6
0 -4
0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90
θ (deg) θ (deg)
Figure 4.35. a) Shape of the different quadratic isosceles-triangular elements used in the mesh;
Variation of the longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence
for various value of φ plotted in a polar (b and c) and linear (d and e) fashion.
In this case, the mesh density is kept the same in the y (vertical) and x (horizontal)
direction at 10 and 11.55 elements per wavelength respectively while the element is
sheared in the y direction. The shearing angle varies from 0 to 45 degrees. Results of
the variation of the velocity error against the angle for a range of mesh density are
shown on Figure 4.36 and agree extremely well with results obtained in the explicit
case.
134
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
φ
φ φ
φ
90 2 90
b) c)
60 60
2 γ=0 deg
θ (deg) 0 θ (deg)
γ=15 deg
γ=30 deg
Longitudinal velocity error (%)
30 γ=45 deg -2 30
1
0 0 -6 0
-30 -30
-60 -60
-90 -90
3 2
Longitudinal velocity error (%)
d) e)
0
Shear velocity error (%)
-2
γ -4
2
-6 γ
-8
-10
1 -12
-14
-16
0 -18
-80 -60 -40 -20 0 20 40 60 80 -80 -60 -40 -20 0 20 40 60 80
θ (deg) θ (deg)
Figure 4.36. a) Shape of the different scalene-triangular elements used in the mesh; Variation of the
longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence for various
values of γ plotted in a polar (b and c) and linear (d and e) fashion.
4.3.4.4 Conclusion
Results obtained with all shapes of linear triangular elements in the implicit case match
almost perfectly those obtained with the explicit case with low CFLX. This highlights
that the solving scheme has a limited influence compared to the element type used.
135
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
The error on both shear and longitudinal velocities is monitored over a range of angles
from 0 to 90 degrees for a range of mesh densities. Results are presented in Figure 4.37.
90
Longitudinal velocity error (%)
N= b) 90
a)
60 6
Figure 4.37. Variation of the longitudinal (a and c) and shear (b and d) velocity error against the angle
of incidence for various mesh density plotted in a linear (a and b) and polar (c and d) fashion.
The error on velocity measured in this case is much lower than that observed with
linear triangular elements. The absolute value of the error on velocity is lower than
0.4% for a mesh density greater than 8. As for the case of quadratic quadrilateral
elements, the measured velocity is higher than the theoretical one. The error is the same
for both wave types at 0 and 60 degrees. The velocity variation is not the same between
these two points as it increases for shear waves but decreases for longitudinal waves.
As for the quadratic quadrilateral elements, the prediction of the error devised for
linear and modified quadratic elements cannot be applied. The fact that the error is so
low means that the error is practically negligible.
136
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
The meshes studied in this part are “stretched” in one direction to create isosceles
triangle elements. The element size in the vertical direction y is kept the same for all
models so that the mesh density N in this direction is 10. The element size in the
horizontal direction x is varied so that the mesh density N varies from 5.35 to 20.
Results are plotted on Figure 4.38.
φ φ φ φ φ φ
Longitudinal velocity error (%)
90 90
b) c)
60 60
30 30
-0.4 -0.4
0 0 0 0
0 0
d)
Error on shear velocity (%)
-0.1 e) -0.1
-0.2 -0.2
-0.3 -0.3
φ
-0.4 -0.4
-0.5 -0.5
-0.6 -0.6 φ
-0.7 -0.7
-0.8 -0.8
-0.9 -0.9
0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90
θ (deg) θ (deg)
Figure 4.38. a) Shape of the different quadratic isosceles-triangular elements used in the mesh;
Variation of the longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence
for various values of φ plotted in a polar (b and c) and linear (d and e) fashion.
Deformation of the elements leads to a modification of the error curves. The maximum
error is in line with what is expected and no loss of precision seems to occur in this
case.
137
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
In this case, the mesh density is kept the same in the y (vertical) and x (horizontal)
direction at 10 and 11.55 elements per wavelength respectively while the element is
sheared in the y direction. The shearing angle varies from 0 to 45 degrees. Results of
the variation of the velocity error against the angle for a range of mesh density are
shown on Figure 4.39.
φ
φ φ
φ
90 90
b) c)
60 60
-0.6 γ=0 deg
θ (deg) 1 θ (deg)
γ=15 deg
γ=30 deg
Longitudinal velocity error (%)
30 γ=45 deg 30
0
Shear velocity error (%)
0
0.2 0 -0.2 0
-30 -30
-60 -60
0.2 0.5
d)
Error on shear velocity (%)
0.1 e) 0
0 -0.5
-0.1 -1.0
-0.2 -1.5
γ
-0.3 -2.0
-0.4 γ -2.5
-0.5 -3.0
-0.6 -3.5
-0.7 -4.0
-80 -60 -40 -20 0 20 40 60 80 -80 -60 -40 -20 0 20 40 60 80
θ (deg) θ (deg)
Figure 4.39. a) Shape of the different scalene-triangular elements used in the mesh; Variation of the
longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence for various
values of γ plotted in a polar (b and c) and linear (d and e) fashion.
138
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
Results for both wave types up to case 3 show a relatively good behaviour but results
for case 4 clearly show a great increase in the velocity error along the longest triangle
edge. The absolute error reaches a value of 3.5% which is much higher than expected.
4.3.5.4 Conclusion
Quadratic triangular elements (CPE6 [14]) have been shown to provide good results
for equilateral and isosceles triangle shapes. Velocity error in these cases is very low
and is close to the amplitudes observed with quadratic quadrilateral elements. The
shearing of the element did not prove problematic up to the case of the right angle
triangle. For larger shearing, the velocity error dramatically increased along the longest
edge of the triangles. This phenomenon is similar to what was observed for linear
triangular elements. Element deformation is acceptable but should be limited.
The error on both shear and longitudinal velocity is monitored over a range of angles
from 0 to 90 degrees for a range of mesh densities. Results are presented in Figure 4.40.
90 90
a) b)
Longitudinal velocity error (%)
60 N= 60
4 7 4
Shear velocity error (%)
θ (deg) 8 θ (deg)
9
30 10 30
2 11 2
12
15
20
0 30 0
0 0
5 d) 5
c)
N
Longitudinal velocity error (%)
4 4
Shear velocity error (%)
N
3 3
2 2
1 1
0 0
0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90
θ (deg) θ (deg)
Figure 4.40. Variation of the longitudinal (a and c) and shear (b and d) velocity error against the angle
of incidence for various mesh densities plotted in a linear (a and b) and polar (c and d) fashion.
139
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
It is interesting to note that the results agree very well with the ones obtained with the
same element type solved with an explicit scheme. The amplitude of the error for the
modified quadratic triangle elements is much higher than the quadratic triangle or
quadrilateral elements and is generally in line with the amplitude observed for linear
elements. Velocity error prediction as defined in Section 4.2.4 is valid in this case.
The meshes studied in this part are “stretched” in one direction to have isosceles
triangle elements. The element size in the vertical direction y is kept the same for all
models so that the mesh density N in this direction is 10. The element size in the
horizontal direction x is varied so that the mesh density N varies from 5.35 to 20.
Results are plotted on Figure 4.42. As for the equilateral case, an excellent agreement
exists between the explicit and implicit cases.
φ φ φ φ φ φ φ
8 90
Error on longitudinal velocity (%)
b) 8 c)
60 60
Error on shear velocity (%)
φ=
6 θ (deg) 30 6 θ (deg)
40
4 30 50 4 30
60
70
2 80 2
90
0 0 0 0
9 9
e)
Error on longitudinal velocity (%)
d) 8 8
Error on shear velocity (%)
7 7
6 6
φ φ
5 5
4 4
3 3
2 2
1 1
0 0
0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90
θ (deg) θ (deg)
Figure 4.41. a) Shape of the different quadratic isosceles-triangular elements used in the mesh;
Variation of the longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence
for various values of φ plotted in a polar (b and c) and linear (d and e) fashion.
140
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
In this case, the mesh density is kept the same in the y (vertical) and x (horizontal)
direction at 10 and 11.55 elements per wavelength respectively while the element is
sheared in the y direction. The shearing angle varies from 0 to 45 degrees. Results are
plotted on Figure 4.42. Comparing implicit and explicit results shows an almost perfect
superimposition of results.
φ
φ φ
φ
b) 4 90 c) 4 90
60 60
γ=0 deg
3 θ (deg) 3 θ (deg)
γ=15 deg
γ=30 deg
Error on longitudinal velocity (%)
2 30 γ=45 deg 2 30
Error on shear velocity (%)
1 1
0 0 0 0
-30 -30
-60 -60
-90 -90
4 4
e)
Error on longitudinal velocity (%)
d)
Error on shear velocity (%)
3 3
γ
2 2 γ
1 1
0 0
-80 -60 -40 -20 0 20 40 60 80 -80 -60 -40 -20 0 20 40 60 80
θ (deg) θ (deg)
Figure 4.42. a) Shape of the different scalene-triangular elements used in the mesh; Variation of the
longitudinal (b and d) and shear (c and e) velocity error against the angle of incidence for various
values of γ plotted in a polar (b and c) and linear (d and e) fashion.
141
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
4.3.6.4 Conclusion
As for the linear triangular elements, results for modified quadratic triangular elements
show an excellent agreement between explicit and implicit results. It is interesting to
note that this element type does not suffer from a great loss of precision due to the
deformation of the element. Modified triangular elements behave in a way closer to
linear quadrilateral elements than linear or quadratic triangular ones.
4.4. Conclusions
The study of reduced integration linear quadrilateral elements solved with an explicit
solver showed the influence of the scaled Courant number CFLX and mesh density N
on the velocity. Patterns were identified and permitted to define a fitted analytical
expression to evaluate with high accuracy the amplitude of the velocity error Ec at any
angle, any scaled Courant number and for a mesh density higher than 6:
2
( 1 – CFLX )-
E c = 180 -------------------------------
2
(4·24)
N
Following this, the study of deformed elements highlighted that the maximum error in
a mesh could be predicted using the fitted analytical expression defined for square
elements (see Equation 4·24). This indicates that this element type is not the subject of
dramatic loss of performance when deformed. Therefore using deformed quadrilateral
elements in a model is not an issue on its own but other factors should also be
considered. As quadrilateral elements are deformed, the maximum error will
consequently increase with the longest diagonal. Therefore, for a given maximum
acceptable error, the overall node density of a model will increase with the level of
element deformation. In other words, the number of degrees of freedom necessary to
cover the geometry of a model with deformed elements will be greater than with square
ones. On top of this, the reduction in the length of the shortest diagonal leads to a
reduction in the stable increment meaning that, for the same number of degrees of
freedom, a model with deformed elements can not be solved as quickly as one made of
square elements.
142
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
With linear equilateral triangular elements, it was demonstrated that the error could be
correctly predicted using rules established for quadrilateral elements (see Equation
4·24) along the altitudes and element sides for longitudinal waves and along the
element sides for shear waves. The influence of the scaled Courant number CFLX and
mesh density N was confirmed at these positions and wave types. An error of zero was
measured for shear waves along the altitudes of the triangular elements. Deformed
mesh cases showed that this element type is sensitive to high levels of deformations as
the predicted error was more than tripled in one case (16.4% instead of 4.8%). Using
highly deformed linear triangular elements is therefore not advised. The author of this
work recommends that linear triangular element internal angles should all have values
higher than 50 degrees and lower than 70 degrees in order not to experience significant
loss of precision.
For implicit solving, full integration linear square element results for longitudinal
waves showed very good agreement with the ones obtained with reduced integration
elements and the explicit solver. For shear waves, the plot had the same shape but the
error was marginally smaller (factor of 0.75). In a similar way to the explicit case, the
following fitted analytical expression can be used to define the velocity error:
180-
E c = --------
2
(4·25)
N
143
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
elements cope well with being deformed and the use of deformed ones in a model is
not a problem.
For a given mesh density, quadratic quadrilateral elements proved to strongly reduce
the velocity error to less than 0.5% for a node density higher than 8. Therefore, the
velocity error with this type of element is not a cause for concern compared to other
cases. Deformation of elements did not prove to be a problem. The two weaknesses of
this element type are that in ABAQUS it can only be used with an implicit solver and
that, as it is a quadratic element, geometric resolution (defined by the smallest
geometric size of an element) is halved compared to linear elements for a given node
density.
Linear and modified quadratic elements results for the implicit case matched almost
exactly results obtained with the explicit solver and therefore the findings are the same.
Quadratic equilateral triangular elements were shown to lead to low velocity error as
were quadratic square elements, but proved to be sensitive to high levels of
deformation in the same way as linear triangular elements. Therefore the author
recommends limiting the internal angle range to 50 and 70 degrees.
Overall, the influence of the solver is ruled out, as element type and shapes are the
factor determining the velocity error. Unlike quadrilateral elements, all triangular
elements apart from modified quadratic ones were shown to be sensitive to element
deformation. This justifies why automatic meshing algorithms seek to minimize the
deformation and why modellers need to exercise caution when using linear and
quadratic triangular elements.
The angle of propagation in a regular mesh was shown to play an important part in the
amplitude of the velocity error as it leads to a non-negligible variation.
The mesh density was shown to be the driving factor influencing the velocity error for
all element types. The influence is not so noticeable for quadratic elements as the
velocity is highly accurate for any mesh density. For linear and modified quadratic
elements, its influence is much stronger.
144
Chapter 4
On the influence of mesh parameters on elastic bulk wave velocities
145
Chapter 5
Accurate modelling of defects using
Finite Elements
5.1. Introduction
In the same way as a stress analysis engineer studying the stresses in a mechanical
component, an NDT engineer developing an ultrasonic inspection technique has to
make assumptions and simplifications in order to be able to numerically model the
interaction of a wave with a structural defect. In the past, the vast majority of wave
propagation work focused on strongly simplified defects. One reason for this is that
initial work is bound to focus on the simplest cases, but many modellers are ready to
go on to more complex problems. As computer power has developed exponentially,
new opportunities have arisen and have highlighted the fact that only a limited amount
of information is available regarding how to accurately model the interaction of elastic
waves with complex defects using FE.
One specific issue related to modelling wave propagation using FE is that a change in
element shape or size leads to a change of acoustic impedance that generates unwanted
reflections. As these reflections pollute the signal of interest, most NDT modellers
favour the use of regular meshes usually constituted of linear square elements. In this
context, it is well known that the correct representation of the wave propagation
phenomenon requires the mesh density to be at least 7 nodes per shortest wavelength.
When using a regular square mesh, defects and features are represented by adding or
removing elements. This can be thought of as a “LEGO” approach where the elements
are LEGO bricks. Clearly, the down side to this is that the geometric definition of
features and defects is dependent on the mesh size and orientation. As a consequence,
theoretical geometric edges not aligned with the mesh are approximated and
sometimes represented in a jagged fashion. This leads to discrepancies between the
theoretical and modelled geometries.
146
Chapter 5
Accurate modelling of defects using Finite Elements
It is generally accepted - although not thoroughly verified - that refinement of the mesh
leads to a convergence of the FE results towards the exact results. In some cases, mesh
refinement is applied in order to correctly model defects. The direct consequence of
this is an increase in the number of degrees of freedom solved and hence the
computational requirements of the model. Keeping the classic regular square mesh
approach, this raises the question of how fine a mesh needs to be refined - if at all - for
a particular defect or feature to be realistically represented from a wave propagation
perspective.
Another approach which has been less popular in the past is to use automatic meshing
algorithms using triangular elements. These are now widely available in commercial
FE packages. Aside from the lack of availability in the past, one reason why these
algorithms have not been commonly used in NDT wave modelling is that, in most
cases, elements are deformed when automatically generated. As deformation of
elements leads to a change of acoustic impedance and consequently reflections, the
mesh structure is key to appropriate modelling of wave propagation. On the positive
side, when using an automatic meshing algorithm, the geometry is represented more
closely than with the regular mesh approach presented above as the edges of the
elements follow the edges of the geometry. It would be beneficial to perform a
meticulous study to verify whether this approach is well adapted for the modelling of
the interaction of elastic waves with complex defects.
In this chapter, the two meshing approaches described above are tested for a range of
defects, element types and solving methods with both elastic wave types. The general
approach is to refine the mesh size for all these cases until an acceptable reference is
obtained in order to evaluate the evolution of the error as the mesh is refined. The
defects modelled included straight edges, circular defects of various sizes and cracks
of various sizes. These cases are chosen as they are commonly modelled in NDT and
should provide the foundations to more complex defect modelling.
147
Chapter 5
Accurate modelling of defects using Finite Elements
As in the previous chapter, models are non-dimensional and the elastic material is
chosen so that the longitudinal and shear wavelengths are 2 and 1 units respectively:
Young’s modulus is 8/3, Poisson’s ratio is 1/3 and density is 1.The cases investigated
here are the reflection of bulk waves from:
• straight edges
• straight cracks at an angle of 80 degrees and of length 0.25, 1 and 4 units
• circular holes of diameter 0.25, 1 and 4 units
Models are 2D plane strain. Both shear and longitudinal wave generations are
considered, as models are excited by body forces being applied on four elements as
shown on Figure 5.1. Body forces rather than point forces are used in this Chapter as
they offer better mode purity. This generates an almost pure wave - shear or
longitudinal - and as the size of the source is only a fraction of a wavelength is a close
representation of a point source. The point located at the centre of the four elements is
referred to as the excitation point. As the body force load is dependent on the area of
the elements excited, the load is normalised so that the amplitude of the incident wave
is the same in all cases studied in this chapter.
a) b) c) d)
Figure 5.1. (a) Longitudinal and (b) shear wave excitation for a square element mesh and (c)
longitudinal and (d) shear excitation for a triangular element mesh.
Models are run with ABAQUS 6.6-1 in both the time domain (explicit) and the
frequency domain (implicit). For the latter, the excitation is applied as a real harmonic
excitation whose frequency is 1; the complex displacement field is obtained directly
and used later to plot the absolute displacement field. For the explicit case, in order to
obtain an absolute displacement field similar to the frequency domain one, it is decided
to run two separate models: one with cosine and one with sine excitation of unit
frequency . The complex displacement field is obtained by adding the result of the
148
Chapter 5
Accurate modelling of defects using Finite Elements
cosine case to the result of the sine case multiplied by i. The first half of a Hanning
window is used over the first 5 cycles in order to ramp up the amplitude of the sine and
cosine to a constant value in order to avoid exciting a broad range of frequencies. The
model is run sufficiently long so that the steady state is reached at all locations where
the displacements are monitored.
Implicit models are run using linear square elements with full (CPE4) and reduced
(CPE4R) integration, linear triangular elements (CPE3), quadratic triangular elements
(CPE6) and modified quadratic triangular elements (CPE6M). Explicit models are
only run using CPE4R, CPE3 and CPE6M as they are the only elements available for
this type of solver.
Appropriately dimensioned ALID are used for all models. In the case of square
meshes, the width of each sub layer in the ALID is one element. For triangular mesh
cases, as the model is generated using the interactive user interface of ABAQUS CAE
6.6-1, a fixed number of sub layers of width 0.25 are defined. Free meshing of these
layers is then performed.
In order to evaluate the influence of the mesh density on the accuracy of the results, the
following numbers of nodes per shortest wavelength are investigated: 6, 10, 15, 20, 25
and 30. This parameter is identified as the mesh density and is assigned the letter N. In
the case of regular square meshes, N is exact over the whole model. For triangular
element meshes generated by the free meshing algorithm of ABAQUS CAE 6.6-1, the
node spacing is specified on the internal and external boundaries of the model. A
variation in node spacing over the model exists but is limited overall. This leads to a
limited variation in the mesh density. There are two cases of interest where significant
local deviations can exist between the prescribed and actual node spacing: small linear
edges and small circular holes. When a linear edge is smaller than the node spacing, it
is defined by a single element. The node spacing on the edge is directly dependent on
the edge length and locally over-rules the prescribed node spacing condition. All
circular holes receive a special treatment that guarantees that they are not incorrectly
represented. The default setting of the free meshing algorithm of ABAQUS CAE 6.6.1
[14] means that a circular hole is defined by at least 8 nodes. This condition is used in
this work and locally over-rules the prescribed node spacing condition for small holes.
149
Chapter 5
Accurate modelling of defects using Finite Elements
The first case investigated is that of the reflection from a straight edge. The geometry
of the model is presented in Figure 5.2.a.
a) b)
Area of study
ALID
Excitation point
Monitoring line
Monitoring direction
24 Shear waves
6
Monitoring direction
y Longitudinal waves
x
a) b) c)
Figure 5.3. a) square mesh at 0 degrees aligned with the edge, b) square mesh at 45 degrees, c)
triangular free mesh
As shown on Figure 5.2.a, the excitation point is located 6 units away from the edge (6
shear wavelengths or 3 longitudinal ones). As the edge is jagged when the mesh is at
45 degrees, an error on the exact position of the edge exists. In this study, the outside
limit is used as the reference. The complex displacement field is monitored up to 24
units away from the excitation point along a line coincident with the excitation point
and parallel to the edge (see Figure 5.2). For a longitudinal wave excitation, the
displacements are monitored in the x direction and for a shear wave excitation, the
150
Chapter 5
Accurate modelling of defects using Finite Elements
Results obtained for the implicit models for longitudinal and shear wave excitation are
presented in Figure 5.4 and 5.5 respectively. Results for the explicit models for
longitudinal and shear waves are presented in Figure 5.6 and 5.7. This represents a very
large amount of data. In order to simplify the identification of cases, each plot is
associated with a letter and a number (e.g. D15). The number is the mesh density of the
case (row) and the letter identify the excitation, mesh and element types as indicated
on the column of the figures. A letter on its own (e.g. D) indicates the whole column.
The first point to be noted is that the results for all cases with a density of 30 agree very
well with each other and that, in most cases, the difference between the 25 and 30 cases
is negligible. These plots can therefore be confidently taken as references. In order to
highlight the discrepancies between these cases (i.e. the reference) and all other cases
where the density is lower than 30, the reference curves are plotted in red for each case
modelled.
151
Chapter 5
Accurate modelling of defects using Finite Elements
20
G) CPE6M 0deg
0
20 -20
F) CPE6 0deg
0
20 -20
E) CPE3 0deg
0
20 -20
D) CPE4R 45deg
0
20 -20
C) CPE4 45deg
0
20 -20
B) CPE4R 0deg
0
20 -20
A) CPE4 0deg
0
-20
1
0
1
0
1
0
1
0
1
0
1
|U| N=6 |U| N=10 |U| N=15 |U| N=20 |U| N=25 |U| N=30
Figure 5.4. Implicit models for straight edge: Monitored absolute displacement for a longitudinal wave
excitation using CPE4 and CPE4R meshes at 0 degrees, CPE4 and CPE4R meshes at 45 degrees and
CPE3, CPE6 and CPE6M triangular elements. Thin red line is reference for N=30 for each case.
152
Chapter 5
Accurate modelling of defects using Finite Elements
20
N) CPE6M 0deg
0
20 -20
M) CPE6 0deg
0
20 -20
L) CPE3 0deg
0
20 -20
K) CPE4R 45deg
0
20 -20
J) CPE4 45deg
0
20 -20
I) CPE4R 0deg
0
20 -20
H) CPE4 0deg
0
-20
5
0
5
0
5
0
5
0
5
0
5
|V| N=6 |V| N=10 |V| N=15 |V| N=20 |V| N=25 |V| N=30
Figure 5.5. Implicit models for a straight edge: Monitored absolute displacement for a shear wave
excitation using CPE4 and CPE4R meshes at 0 degrees, CPE4 and CPE4R meshes at 45 degrees and
CPE3, CPE6 and CPE6M triangular elements. Thin red line is reference for N=30 for each case.
153
Chapter 5
Accurate modelling of defects using Finite Elements
|U| N=6
0
1
|U| N=10
0
1
|U| N=15
0
1
|U| N=20
0
1
|U| N=25
0
1
|U| N=30
0
-20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.6. Explicit models for a straight edge: Monitored absolute displacement for a longitudinal
wave excitation using CPE4R meshes at 0 degrees, CPE4R meshes at 45 degrees and CPE3 and
CPE6M triangular elements. Thin red line is reference for N=30 for each case.
A and B show the results obtained with square element meshes at 0 degrees. The
difference between using full (A) or reduced (B) integration elements is negligible for
all mesh densities and they are therefore assessed jointly. It can be seen that there is a
good convergence towards the reference case (30) with an error of less than 3% for a
density of 20. A mesh of density 6 gives inaccurate results but increasing this value to
10 leads to a relatively correct representation of the phenomena occurring despite a
large amplitude error.
C and D show the results obtained with square element meshes at 45 degrees. Results
for fully integrated elements (C) are very similar to the results obtained with a mesh at
0 degrees (A) whereas reduced integration element ones (D) differ noticeably. In this
case, one can see that the amplitude of side lobes oscillates on each side of the
154
Chapter 5
Accurate modelling of defects using Finite Elements
|V| N=6
0
5
|V| N=10
0
5
|V| N=15
0
5
|V| N=20
0
5
|V| N=25
0
5
|V| N=30
0
-20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.7. Explicit models for a straight edge: Monitored absolute displacement for a shear wave
excitation using CPE4R meshes at 0 degrees, CPE4R meshes at 45 degrees and CPE3 and CPE6M
triangular elements. Thin red line is reference for N=30 for each case.
reference ones as the mesh is refined. This indicates that, despite a fine mesh (D25),
the stair-cased edge interferes with the correct representation of the reflection.
E, F and G show the results obtained with automatically generated triangular element
meshes. A convergence towards the reference occurs in the same way as for square
elements results (A). Quadratic elements (F) offer the best results with a good
agreement achieved for a mesh density of 15. Linear triangular elements (E) cause an
oscillation to occur on the side lobes and a mesh density of 25 (E25) is required to have
acceptable results. A similar, although milder, phenomenon occurs with modified
quadratic elements (G).
155
Chapter 5
Accurate modelling of defects using Finite Elements
Results for shear wave excitation are presented in Figure 5.5 and show a similar
behaviour to longitudinal waves except from the lack of oscillation (observed in E and
G but absent in L and M).
Explicit results presented in Figure 5.6 and 5.7 agree very well with the results from
the implicit solver with their respective element type. In other words, the results of the
CPE4R element type are the same when solved with an implicit (B) or an explicit
method (Q). Consequently the findings obtained for the implicit solver are also valid
for the explicit one.
Overall, based on the plots for square element meshes at 0 degrees, it can be said that
this case provides superior results to square element meshes at 45 degrees and
triangular element meshes for the modelling of the reflection from a straight edge. The
use of square element meshes at an angle gives acceptable results but the oscillation of
the side lobes is a concern and could be explained by diffraction occurring due to the
jagged nature of the edge. Quadratic triangular elements give results close to square
elements with a mesh at 0 degrees. Linear triangular elements need a finer mesh than
the other triangular element types in order to provide accurate results and the
noticeable undulation on the side lobes is a concern. Modified quadratic element
meshes seem to offer a compromise between linear and quadratic triangular elements.
In this part, the longitudinal and shear wave reflection from a crack at an angle are
studied. It should be noted that the angle and length of the modelled cracks have not
particular significance and are just sensibly chosen examples. The geometry of the
models is presented in Figure 5.8.a. The mid-point of the crack is located 10 units away
from the excitation point and the crack is at an angle of 10 degrees relative to the
vertical. Three crack lengths are used: 0.25, 1, 4. This permits a wide range of
reflection regimes for both wave types.
156
Chapter 5
Accurate modelling of defects using Finite Elements
a) b)
Area of study
ALID
Excitation point
10
Monitoring line
Monitoring direction
24 24 Shear waves
o Monitoring direction
y 10 Longitudinal waves
x
Two meshing strategies are used. On one side, a regular square mesh is used and nodes
are disconnected in order to create a crack following the edges of the elements in the
mesh. It can be seen as creating a free edge following the lines of the regular square
grid. This is achieved by using a MATLAB routine created by a summer placement
student, Sumeet Kale, under the guidance of the author. It defines a crack as close as
possible to the exact crack within the particular regular square mesh. On the other side,
a triangular element mesh generated by the free meshing algorithm from ABAQUS
CAE 6.6-1 is used. In this situation, the geometry of the crack is determined first and
then free meshing of the geometry is performed. The particular shape and mesh of each
crack generated with both methods is shown in each section where the particular
straight crack geometry is studied.
157
Chapter 5
Accurate modelling of defects using Finite Elements
The first case studied is the case of a crack of unit length. The actual cracks used for
these models are represented in Figure 5.9. When using a triangular mesh, the crack
geometry is exactly defined whereas the crack obtained with the square mesh differs
from the exact crack. This is highlighted in Figure 5.9 where the theoretical crack as
defined in Figure 5.8 is shown in red and the actual modelled crack in blue. Both lines
superimpose for triangular element meshes but do not for the regular square element
meshes. In both cases, as the mesh is refined, the number of elements defining the
crack increases. This should increase the accuracy of the reflection.
CPE3 N=6 CPE3 N=10 CPE3 N=15 CPE3 N=20 CPE3 N=25 CPE3 N=30
CPE6/6M N=6 CPE6/6M N=10 CPE6/6M N=15 CPE6/6M N=20 CPE6/6M N=25 CPE6/6M N=30
CPE4/4R N=6 CPE4/4R N=10 CPE4/4R N=15 CPE4/4R N=20 CPE4/4R N=25 CPE4/4R N=30
Figure 5.9. Definition of unit long cracks with triangular and square meshes. Blue line shows modelled
crack and red line theoretical crack (which is the same line with triangular element meshes but not with
regular square element meshes).
Results obtained for the implicit models are presented in Figure 5.10 and 5.11. Results
for the explicit models for longitudinal and shear waves are presented in Figure 5.12.
In order to simplify the identification of cases, each plot is associated with a letter and
a number as in the previous sections.
158
Chapter 5
Accurate modelling of defects using Finite Elements
|U| N=6
2
0
4
|U| N=10
2
0
4
|U| N=15
2
0
4
|U| N=20
2
0
4
|U| N=25
2
0
4
|U| N=30
0
-20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.10. Implicit models for a crack of unit length: Monitored absolute displacement for a
longitudinal wave excitation using mesh made of CPE3, CPE6, CPE6M, CPE4 and CPE4R elements.
Thin red line is reference for N=30 for each case.
0
30
|V| N=10
0
30
|V| N=15
0
30
|V| N=20
0
30
|V| N=25
0
30
|V| N=30
0
-20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.11. Implicit models for a crack of unit length: Monitored absolute displacement for a shear
wave excitation using mesh made of CPE3, CPE6, CPE6M, CPE4 and CPE4R elements. Thin red line
is reference for N=30 for each case.
159
Chapter 5
Accurate modelling of defects using Finite Elements
|V| N=6
|U| N=60 0
4 30
|V| N=10
|U| N=10
0 0
4 30
|V| N=15
|U| N=15
0 0
4 30
|V| N=20
|U| N=20
0 0
4 30
|V| N=25
|U| N=25
0 0
4 30
|V| N=30
|U| N=30
0 0
-20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.12. Explicit models for a crack of unit length: Monitored absolute displacement for a shear
and longitudinal wave excitation using mesh made of CPE3, CPE6M and CPE4R elements. Thin red
line is reference for N=30 for each case.
Looking at the implicit results for a longitudinal wave excitation presented in Figure
5.10, it can be seen that, for all element types, there is a good convergence towards the
models with 30 elements per wavelength. In this case, the displacement fields are
almost perfectly matching for all element types. For all triangular elements, a mesh
density of 20 gives results very close to the reference. CPE6 elements (B) provide
higher quality than CPE3 (A) and CPE6M (C), with acceptable results being obtained
for N=10. Square elements (D-E) - in particular the reduced type (E) - have a slower
convergence. For shear wave excitation, the results in Figure 5.11 show a similar
pattern, although the rate of convergence for all element types is slower. A mesh
density of 25 is necessary for the error on the displacement amplitude to be negligible.
This can be explained by the fact that, for a given mesh density, the shear wave velocity
error is higher than the longitudinal one but also that the shear wavelength is shorter
than the longitudinal one making it more sensitive to numerical “defects” (changes in
element size in triangular meshes, steps in square meshes).
160
Chapter 5
Accurate modelling of defects using Finite Elements
For all mesh densities, explicit results shown in Figure 5.12 (K-P) superimpose almost
perfectly with the results obtained for the implicit ones with the same element type.
These are presented separately for completeness.
Overall, representing an average size straight crack at an angle of 80 degrees does not
seem to be challenging for either meshing strategy. The physics are modelled quite
accurately for any mesh density over 10. The main source of error is likely to be the
actual velocity error and its consequences on the reflection angles. Indeed, as angle of
reflection of a wave depends on the wave propagation velocity, an error on the
propagation velocity leads to an error on the angle of reflection. In the regular mesh
approach, the actual geometric difference between the theoretical crack and the actual
crack accounts for another source of error which can explain the slower convergence.
Once again, it is interesting to note that the use of free meshing with triangular
elements of various sizes does not harm the accuracy of the results.
The second crack case studied is the case of a crack of length 0.25. The actual cracks
used for these models are represented in Figure 5.13. This defect is more challenging
than the previous one as the number of nodes defining the defect is very limited in all
cases with only 9 nodes for the highest density. As in the previous case, when using a
triangular mesh, the crack geometry is exactly defined whereas the crack obtained with
the square mesh differs from the exact crack. As the defect is smaller than in the
previous case, the definition of the crack using a square mesh differs more strongly in
terms of shape but also in terms of length from the theoretical crack. When a square
mesh of density 6 is used, the algorithm does not generate any crack and therefore the
results for this case are of no interest as no reflection occurs.
161
Chapter 5
Accurate modelling of defects using Finite Elements
CPE3 N=6 CPE3 N=10 CPE3 N=15 CPE3 N=20 CPE3 N=25 CPE3 N=30
CPE6/6M N=6 CPE6/6M N=10 CPE6/6M N=15 CPE6/6M N=20 CPE6/6M N=25 CPE6/6M N=30
CPE4/4R N=6 CPE4/4R N=10 CPE4/4R N=15 CPE4/4R N=20 CPE4/4R N=25 CPE4/4R N=30
Figure 5.13. 0.25 unit long crack definition with triangular and square meshes. Blue line shows
modelled crack and red line theoretical crack (which is the same line with triangular element meshes
but not with regular square element meshes).
Results obtained for the implicit models are presented in Figure 5.14 and 5.15. Results
for the explicit models for longitudinal and shear waves are presented in Figure 5.16.
In order to simplify the identification of cases, each plot is associated with a letter and
a number as in the previous section.
Looking at the results for the longitudinal cases presented in Figure 5.14, it is
interesting to see that there is a convergence towards the results obtained with the finest
mesh for both triangular and square meshes, but this convergence is much slower than
the one observed for the crack of unit length. For quadratic triangular meshes of density
20, the error on the amplitude relative to the reference varies between 6 and 8%, but
for a density of 25 this error is almost negligible. The most likely explanation for this
is that independently of the mesh density the improvement in the quality of results is
caused by the increase in the number of nodes defining the crack. Using the free
meshing algorithm, the number of nodes on the crack is 5, 5, 9 and 9 for densities of
15, 20, 25 and 30. It can therefore be said that it is necessary to have at least 7 nodes
over the length of a crack to have an error smaller than 5% relative to the case with a
mesh of density 30. As mentioned previously, the agreement is very good between
162
Chapter 5
Accurate modelling of defects using Finite Elements
|U| N=6
0
|U| N=10
1
0
|U| N=15
0
|U| N=20
0
|U| N=25
0
|U| N=30
0
-20 0 20 -20 0 20 -20 0 20-20 0 20 -20 0 20
Figure 5.14. Implicit models for a crack of length 0.25: Monitored absolute displacement for a
longitudinal wave excitation using mesh made of CPE3, CPE6, CPE6M, CPE4 and CPE4R elements.
Thin red line is reference for N=30 for each case.
0
|U| N=10
0
|U| N=15
0
|U| N=20
0
|U| N=25
0
|U| N=30
0
-20 0 20 -20 0 20 -20 0 20-20 0 20 -20 0 20
Figure 5.15. Implicit models for a crack of length 0.25: Monitored absolute displacement for a shear
wave excitation using mesh made of CPE3, CPE6, CPE6M, CPE4 and CPE4R elements. Thin red line
is reference for N=30 for each case.
163
Chapter 5
Accurate modelling of defects using Finite Elements
|U| N=6
|V| N=6
2
0 0
|U| N=10
|V| N=10
1
2
0 0
|U| N=15
|V| N=15
1
2
0 0
|U| N=20
|V| N=20
1
2
0 0
|U| N=25
|V| N=25
1
2
0 0
|U| N=30
|V| N=30
1
2
0 0
-20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.16. Explicit models for a crack of length 0.25: Monitored absolute displacement for a shear
and longitudinal wave excitation using mesh made of CPE3, CPE6M and CPE4R elements. Thin red
line is reference for N=30 for each case.
results for N=25 and N=30 for quadratic (B) and modified (C) triangular element types
indicating that the physical phenomena are very accurately represented. The fact that
the convergence is not as good with linear triangular elements (A) and that there is an
error of 13% between the linear and quadratic element cases for N=30 means that with
linear triangular elements the mesh needs to be refined further in order to model the
physical phenomena very accurately. Therefore, one can say that quadratic triangular
elements offer superior results to linear ones when modelling small defects.
For square elements (D-E), cases for density 6 (no crack) and 10 (short vertical crack)
are ruled out. Convergence towards the reference occurs for both square element types.
Unexpectedly, the reference for the reduced integration elements (E) matches more
closely our best reference (for triangular quadratic elements) than the full integration
ones. No strong convergence has been reached for N=30; it is therefore difficult to
draw definite conclusions on the quality of the representation using square elements. It
can nevertheless be noted that, despite the coarse geometry of the crack and the
inaccuracy on the exact length of the crack, the reflection is relatively well modelled.
164
Chapter 5
Accurate modelling of defects using Finite Elements
The findings from shear wave excitation for all element types are in line with the
longitudinal ones.
For all mesh densities, explicit results in Figure 5.16 superimpose closely with the
results obtained from the implicit ones with the same element type and are presented
separately for completeness.
Overall, it is clear that not only the mesh density but also the number of nodes used to
define a short straight crack at an angle is crucial to accurately represent its reflections
from both longitudinal and shear waves. Quadratic triangular meshes offer a quick
convergence towards accurate results. The use of square meshes leading to only
roughly represented cracks give acceptable results. These are not as accurate as the
ones obtained with quadratic triangular elements. Linear triangular meshes would need
to be refined over 30 nodes per shortest wavelength to obtain the same result quality as
with quadratic triangular elements. The use of quadratic triangular elements is
therefore preferable to linear triangular and square ones in this case.
The last crack case studied is the case of a crack of length 4. The actual cracks used for
these models are represented in Figure 5.17. This defect is not as challenging as the two
previous ones in terms of correctly defining the geometry of the crack. With a square
mesh approach even with a mesh density of 6, the crack is well represented geometri-
cally.
Results obtained for the implicit models for longitudinal and shear wave excitation are
presented in Figure 5.18 and 5.19 respectively. Results for the explicit models for
longitudinal and shear waves are presented in Figure 5.20.
For both longitudinal and shear wave excitation, all element types show a quick
convergence towards their reference. The references are very closely matched for all
element types. For triangular elements, a mesh density of 10 gives acceptable results
and a mesh density of 15 limits the error on the amplitude to 2%. The error is negligible
for N superior or equal to 20. For square elements, the criterion is slightly more
stringent, with results for a mesh density of 20 being acceptable.
165
Chapter 5
Accurate modelling of defects using Finite Elements
CPE3 N=6 CPE3 N=10 CPE3 N=15 CPE3 N=20 CPE3 N=25 CPE3 N=30
CPE6/6M N=6 CPE6/6M N=10 CPE6/6M N=15 CPE6/6M N=20 CPE6/6M N=25CPE6/6M N=30
CPE4/4R N=6 CPE4/4R N=10 CPE4/4R N=15 CPE4/4R N=20 CPE4/4R N=25 CPE4/4R N=30
Figure 5.17. 4 unit long crack definition with triangular and square meshes. Blue line shows modelled
crack and red line theoretical crack (which is the same line with triangular element meshes but not with
regular square element meshes).
0
10
|U| N=10
0
10
|U| N=15
0
|U| N=20
10
0
10
|U| N=25
0
10
|U| N=30
0
-20 0 20 -20 0 20 -20 0 20 -20 0 20-20 0 20
Figure 5.18. Implicit models for a crack of length 4: Monitored absolute displacement for a
longitudinal wave excitation using mesh made of CPE3, CPE6, CPE6M, CPE4 and CPE4R elements.
Thin red line is reference for N=30 for each case.
166
Chapter 5
Accurate modelling of defects using Finite Elements
|U| N=6
50
|U| N=10 0
50
0
|U| N=15
50
0
|U| N=20
50
0
|U| N=25
50
0
|U| N=30
50
0
-20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.19. Implicit models for a crack of length 4: Monitored absolute displacement for a shear wave
excitation using mesh made of CPE3, CPE6, CPE6M, CPE4 and CPE4R elements. Thin red line is
reference for N=30 for each case.
50
0 0
10
|U| N=10
|U| N=10
50
0 0
10
|U| N=15
|U| N=15
50
0 0
10
|U| N=20
|U| N=20
50
0 0
10
|U| N=25
|U| N=25
50
0 0
10
|U| N=30
|U| N=30
50
0 0
-20 0 20 -20 0 20-20 0 20 -20 0 20 -20 0 20-20 0 20
Figure 5.20. Explicit models for a crack of length 4: Monitored absolute displacement for a shear and
longitudinal wave excitation using mesh made of CPE3, CPE6M and CPE4R elements. Thin red line is
reference for N=30 for each case.
167
Chapter 5
Accurate modelling of defects using Finite Elements
Explicit results for all mesh densities superimpose closely with the results obtained
from the implicit ones with the same element type and are presented separately on
Figure 5.20 for completeness.
In this case, the relaxation of constraints regarding the number of nodes defining the
crack leads to a faster convergence of results for all techniques and element types.
5.4.4 Conclusion
The influence of the mesh density on the accuracy of results is confirmed by the present
study but most significantly the importance of the number of elements defining a defect
is demonstrated. The fact that it is much more difficult to obtain acceptable results for
a short crack than for a long one is easily understandable. It is interesting to note that
a large improvement in result quality was achieved when a limited refinement of the
triangular mesh caused a large increase in the number of nodes defining the smallest
crack (from 5 to 9). It shows that the number of nodes defining a crack has a strong
impact on the accuracy of the reflection. It is also worth noting that the regular square
mesh approach is more adapted to the long crack case than the shorter ones. However,
even with long cracks, a square mesh needs to be finer than a triangular one in order to
achieve the same accuracy. Local mesh refinement would be advantageous in models
where short defects are present if there are no issues with the mesh density transition.
This last point will be examined in the next chapter.
The reflection of longitudinal and shear waves from circular defects is studied. In each
case, a hole is located 10 units away from the excitation point. Three hole diameters
are used: 0.25, 1 and 4 units. This permits to have a wide range of reflection regimes
for both wave types. The general geometry of the models is shown in Figure 5.21.a.
Two meshing strategies are used. On one side, a regular square mesh is used and
elements are removed in order to create a hole. This is achieved by using a MATLAB
routine created by a summer placement student, Sumeet Kale, under the guidance of
the author. It defines a hole as close as possible to the exact hole within the particular
regular square mesh. The criterion used to determine whether to generate an element
or not is the distance between the centre point of a potential square element in the node
168
Chapter 5
Accurate modelling of defects using Finite Elements
a) b)
Area of study
ALID
Excitation point
10
Monitoring line
Monitoring direction
24 24 Shear waves
Monitoring direction
y Longitudinal waves
x
Figure 5.21. Circular defect model: a) with circular defect, b) without circular defect.
grid to the centre of the exact hole. If this distance is lower than the radius, no element
is created. This leads to a hole whose edge is jagged. On the other side, a triangular
mesh generated by the free meshing algorithm from ABAQUS CAE 6.6-1 is used. In
this situation, the geometry of the hole is determined first and then free meshing of the
geometry is performed. When using a triangular mesh, nodes are placed on the edge of
the exact hole and are then used to create elements. The number of elements on the hole
depends on the diameter of the hole. It is never lower than 8 as the curvature control is
set to the default setting of 0.1. It is important to note that with linear triangular
elements, the hole will be made of straight segments whose end points are located on
the exact circle whereas, with quadratic triangular elements, the hole is represented
almost exactly with curved segments.
The displacement field is monitored on a line perpendicular to the line joining the hole
centre point to the excitation point and over a length of 24 units on each side of the
excitation point as shown on Figure 5.21. For a longitudinal wave excitation, the
displacements are monitored in the x direction and for a shear wave excitation, the
displacements are monitored in the y direction. In theory, the incident displacement
field is zero in these directions and only the reflection from the edge is present. In
practice, each case was run with (see Figure 5.21.a) and without (see Figure 5.21.b) the
crack being modelled in order to achieve high quality reflection fields as this resolved
some slight numerical issues (low amplitude noise, oscillation close to the excitation
point). The complex displacement field of the case without the edge is substracted from
the one with the edge. The absolute displacement in the direction specified above is
then plotted.
169
Chapter 5
Accurate modelling of defects using Finite Elements
The first case studied is the case of a hole of unit diameter. The actual holes used for
these models are represented in Figure 5.22. The hole definition is very good with
triangular meshes whereas, with square elements, a noticeable difference exists for
coarse meshes. For a regular square mesh of density 15, two models are run: one using
a diameter of 1 and one with a slightly larger diameter of 1.01. The consequence of this
change in diameter is that an extra 8 elements are not generated in the second case.
These elements are shown in Figure 5.22 and are crossed in order to identify them.
CPE3 N=6 CPE3 N=10 CPE3 N=15 CPE3 N=20 CPE3 N=25 CPE3 N=30
CPE6/6M N=6 CPE6/6M N=10 CPE6/6M N=15 CPE6/6M N=20 CPE6/6M N=25 CPE6/6M N=30
CPE4/4R N=6 CPE4/4R N=10 CPE4/4R N=15 CPE4/4R N=20 CPE4/4R N=25 CPE4/4R N=30
Figure 5.22. Unit diameter hole definition with triangular and square meshes
Results obtained for the implicit models for longitudinal and shear wave excitation are
presented in Figures 5.23 and 5.24 respectively. Results for the explicit models for
longitudinal and shear waves are presented in Figure 5.25.
Results for the triangular meshes indicate a good convergence of the solution for both
longitudinal and shear wave excitation towards a common reference. Acceptable
displacement fields are obtained with 15 nodes per wavelength. A negligible error is
achieved with a mesh density of 20, although some unwanted low amplitude
oscillations are still present for linear and modified quadratic elements.
Looking at the results given by square meshes using the standard algorithm (blue and
red lines), one can see that the general physical phenomena are correctly represented
for N superior to 10 and 15 with longitudinal and shear wave excitation respectively.
170
Chapter 5
Accurate modelling of defects using Finite Elements
|U| N=6 0
3
|U| N=10
0
3
|U| N=15
0
3
|U| N=20
0
3
|U| N=25
0
3
|U| N=30
0
-20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.23. Implicit models for a hole of unit diameter: Monitored absolute displacement for a
longitudinal wave excitation using mesh made of CPE3, CPE6, CPE6M, CPE4 and CPE4R elements.
Thin red line is reference for N=30 for each case.
0
15
|V| N=10
0
15
|V| N=15
0
15
|V| N=20
0
15
|V| N=25
0
15
|V| N=30
0
-20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.24. Implicit models for a hole of unit diameter: Monitored absolute displacement for a shear
wave excitation using mesh made of CPE3, CPE6, CPE6M, CPE4 and CPE4R elements. Thin red line
is reference for N=30 for each case.
171
Chapter 5
Accurate modelling of defects using Finite Elements
|V| N=6
|U| N=6 0 0
3 15
|U| N=10
|V| N=10
0 0
3 15
|U| N=15
|V| N=15
0 0
3 15
|U| N=20
|V| N=20
0 0
3 15
|U| N=25
|V| N=25
0 0
3 15
|U| N=30
|V| N=30
0 0
-20 0 20-20 0 20-20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.25. Explicit models for a hole of unit diameter: Monitored absolute displacement for a shear
and longitudinal wave excitation using mesh made of CPE3, CPE6M and CPE4R elements. Thin red
line is reference for N=30 for each case.
The reference matches the one for triangular elements indicating that the displacement
field is accurate at this density. There is an issue with the amplitude of the signal for N
up to 25. The convergence is not constant as the amplitude error is higher for a density
of 25 and 15 than for one of 20. This is much stronger with shear wave cases than with
longitudinal ones. Two possible sources of error combine to explain this. The first is
that the diameter at 0 and 90 degrees is smaller than the exact ones for densities of 15
(diameter of 0.88) and 25 (diameter of 0.96). This should give lower displacement
amplitude which is not the case here. Against this, with a shear wave excitation, a
density of 20 and an exact diameter at 0 and 90 degrees, the amplitude is lower than
the reference. This shows that the error on the diameter is not the major contributor to
the amplitude error. The low slope of the front face of the hole can not be resolved
accurately by the square mesh and a small variation on the density or diameter leads to
very different shape. At a density of 15, the front face of the hole is made of a flat area
of 8 elements out of the 14 on the diameter (57%) and in the same way, for a density
of 25, a flat area of 10 out of 24 elements (42%) whereas for a density of 20, it is only
6 out of 20 (30%).
172
Chapter 5
Accurate modelling of defects using Finite Elements
In order to look at the influence of this, a hole of diameter 1.01 is generated. The small
change in the exact diameter means that the shape of the front face changes as the
crossed elements in Figure 5.22 are omitted. The flat central part of 8 elements changes
into 2 central elements bordered by 2 flat parts of 3 elements each. The results for this
modified hole are plotted with a green line on Figures 5.23 and 5.24. The change
causes the reflected amplitude at the excitation point to drop by 40% for the shear wave
despite the change in the actual diameter at 0 and 90 degrees from 0.88 to 1.12 units.
It is likely that this dramatic change can be explained by the fact that the front face of
the hole does not behave as a flat reflector in the new configuration. This may also be
the consequence of a change in the reflection interference pattern of the various steps
of the edge of the hole. The shear wave case is far more sensitive to this than the
longitudinal case. This shows that the way the jagged edge of the hole is defined
strongly influences the reflection from a hole in particular in the case of shear wave
excitation.
The second case studied is the case of a hole of diameter 0.25. The actual holes used
for these models are represented in Figure 5.26. The effect of the curvature control is
clearly seen for linear and quadratic triangular meshes up to a density of 15 and 20
respectively. It is worth mentioning that using the default setting, the number of
elements is at least 8 on the edge of the hole. This means that at least 8 and 16 nodes
are used for linear and quadratic elements respectively. The definition of holes using a
square mesh is much approximated. Given the findings of the previous part, this should
be reflected in the results.
Results obtained for the implicit models for longitudinal and shear wave excitation are
presented in Figures 5.27 and 5.28 respectively. Results for the explicit models for
longitudinal and shear waves are presented in Figure 5.29.
The results obtained with quadratic and modified quadratic triangular meshes for both
type of excitation show a good convergence with a negligible difference between N=25
and N=30. Convergence is slower with linear elements and errors still exist between
N=25 and N=30. An error of 4% exists between the reference for quadratic and linear
173
Chapter 5
Accurate modelling of defects using Finite Elements
CPE3 N=6 CPE3 N=10 CPE3 N=15 CPE3 N=20 CPE3 N=25 CPE3 N=30
CPE6/6M N=6 CPE6/6M N=10 CPE6/6M N=15 CPE6/6M N=20 CPE6/6M N=25 CPE6/6M N=30
CPE4/4R N=6 CPE4/4R N=10 CPE4/4R N=15 CPE4/4R N=20 CPE4/4R N=25 CPE4/4R N=30
Figure 5.26. 0.25 unit diameter hole definition with triangular and square meshes
0
|U| N=10
0
|U| N=15
0
|U| N=20
0
|U| N=25
0
|U| N=30
0
-20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.27. Implicit models for a hole of diameter 0.25: Monitored absolute displacement for a
longitudinal wave excitation using mesh made of CPE3, CPE6, CPE6M, CPE4 and CPE4R elements.
Thin red line is reference for N=30 for each case.
elements. This indicates that an even finer mesh is required for the linear elements in
order to represent the reflection very accurately. For square meshes, the reflection
pattern is correctly represented for all mesh densities but, as seen previously, a
variation on the amplitude is caused by the actual diameter of the hole and the shape
174
Chapter 5
Accurate modelling of defects using Finite Elements
|V| N=6 0
10
|V| N=10
0
10
|V| N=15
0
10
|V| N=20
0
10
|V| N=25
0
10
|V| N=30
0
-20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.28. Implicit models for a hole of diameter 0.25: Monitored absolute displacement for a shear
wave excitation using mesh made of CPE3, CPE6, CPE6M, CPE4 and CPE4R elements. Thin red line
is reference for N=30 for each case.
|V| N=6
0 0
10
|U| N=10
|V| N=10
0 0
10
|U| N=15
|V| N=15
0 0
10
|U| N=20
|V| N=20
0 0
10
|U| N=25
|V| N=25
0 0
10
|U| N=30
|V| N=30
0 0
-20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.29. Explicit models for a hole of diameter 0.25: Monitored absolute displacement for a shear
and longitudinal wave excitation using mesh made of CPE3, CPE6M and CPE4R elements. Thin red
line is reference for N=30 for each case.
175
Chapter 5
Accurate modelling of defects using Finite Elements
of the jagged edge. This variation is far from being as strong as in the previous case.
This is likely to be due to the fact that the curvature of the hole is higher relative to the
element size and that the behaviour of the hole is closer to a point scatterer.
The last case studied is the case of a hole of diameter 4. The actual holes used for these
models are represented in Figure 5.30. The curvature control has no effect even for the
coarsest of mesh. The definition of the holes using the square mesh approach provides
a jagged edge which seems to be geometrically close to the exact hole.
CPE3 N=6 CPE3 N=10 CPE3 N=15 CPE3 N=20 CPE3 N=25 CPE3 N=30
CPE6/6M N=6 CPE6/6M N=10 CPE6/6M N=15 CPE6/6M N=20 CPE6/6M N=25 CPE6/6M N=30
CPE4/4R N=6 CPE4/4R N=10 CPE4/4R N=15 CPE4/4R N=20 CPE4/4R N=25 CPE4/4R N=30
Figure 5.30. 4 units diameter hole definition with triangular and square meshes
Results obtained for the implicit models for longitudinal and shear wave excitation are
presented in Figures 5.31 and 5.32 respectively. Results for the explicit models for
longitudinal and shear waves are presented in Figure 5.33.
Results for the triangular meshes show a good convergence for both longitudinal and
shear wave excitation towards a common reference. The error is negligible with a mesh
density of 20, although some unwanted low amplitude oscillation is still present for
linear and modified quadratic elements. For N=15, acceptable results are obtained.
With square elements, longitudinal wave excitation results converge towards the
reference of the triangular elements, but results for shear wave excitation are only
acceptable for N=30. For mesh densities lower than this, the expected single peak at
176
Chapter 5
Accurate modelling of defects using Finite Elements
|U| N=6 0
5
|U| N=10
0
5
|U| N=15
0
5
|U| N=20
0
5
|U| N=25
0
5
|U| N=30
0
-20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.31. Implicit models for a hole of diameter 4: Monitored absolute displacement for a
longitudinal wave excitation using mesh made of CPE3, CPE6, CPE6M, CPE4 and CPE4R elements.
Thin red line is reference for N=30 for each case.
0
30
|V| N=10
0
30
|V| N=15
0
30
|V| N=20
0
30
|V| N=25
0
30
|V| N=30
0
-20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.32. Implicit models for a hole of diameter 4: Monitored absolute displacement for a shear
wave excitation using mesh made of CPE3, CPE6, CPE6M, CPE4 and CPE4R elements. Thin red line
is reference for N=30 for each case.
177
Chapter 5
Accurate modelling of defects using Finite Elements
|U| N=6
|V| N=6
0 0
5 30
|U| N=10
|V| N=10
0 0
5 30
|U| N=15
|V| N=15
0 0
5 30
|U| N=20
|V| N=20
0 0
5 30
|U| N=25
|V| N=25
0 0
5 30
|U| N=30
|V| N=30
0 0
-20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20
Figure 5.33. Explicit models for a hole of diameter 4: Monitored absolute displacement for a shear and
longitudinal wave excitation using mesh made of CPE3, CPE6M and CPE4R elements. Thin red line is
reference for N=30 for each case.
the excitation point is replaced by twin peaks that gradually disappear as the mesh is
refined. This confirms that shear waves are particularly sensitive to the jagged edge
that occurs when modelling a hole with a regular square mesh.
5.5.4 Conclusion
The study of the modelling of the reflection of waves from circular defects showed that
the use of triangular element meshes gives accurate results for small, average and large
holes. The strength of the method was demonstrated when modelling a small hole of
diameter 0.25 (a quarter of the shortest wavelength). The curvature control provided
by the free meshing algorithm set at the default setting (0.1) refined the mesh around
the location of the hole allowing a precise representation of the feature despite the
overall coarseness of the mesh. In this case, excellent results were obtained with
quadratic and modified quadratic elements for a mesh density of 20. A finer mesh is
required to achieve similar quality using linear elements. For other cases, all triangular
elements performed in a similar way providing high quality results for a mesh density
of 20.
178
Chapter 5
Accurate modelling of defects using Finite Elements
It was demonstrated that the use of a regular square mesh to model circular defects is
not appropriate unless a very fine mesh of density superior to 30 is used. Although the
method performed well with longitudinal excitation, shear wave excitation cases were
shown to be particularly sensitive to the way the edge of the hole is defined.
Performance was in fact better for the small hole case than for the average and large
hole cases.
5.6. Conclusions
This study has highlighted the challenges of modelling a range of defects using finite
elements. The use of a regular mesh was shown to perform very well to model a
straight edge when aligned with the mesh, but its inefficiency was demonstrated when
modelling features whose geometries do not follow the grid exactly. The issue of
correctly dimensioning a feature in a fixed grid (e.g exact crack length) length was
highlighted but can be easily resolved by adjusting the mesh density to fit the size of
the defect. It was demonstrated that this is not the main issue influencing the quality of
results. The stair-cased definition of the edges of defects proved to be the crucial
parameter influencing the quality of results. Even with large defects, relatively fine
meshes (N=25) performed poorly in cases dominated by shear waves. Therefore, the
regular mesh approach can only be confidently used to model defects not aligned with
the mesh when the mesh density is extremely fine (N>=30).
Local mesh refinement where only a small area around the defect is extremely finely
meshed would be advantageous in this situation. The use of two meshes of different
densities tied together, which would enable such a refinement to be achieved
economically, is investigated in the next chapter.
Triangular meshes generated by the free meshing algorithm have proved to perform
well in most cases. The fact that the element size is not exactly constant did not cause
a deterioration of results. Dimensioning of the defect was better than with the regular
square mesh approach. As the geometry was defined using the interactive tools of
ABAQUS/CAE before meshing, straight defects were defined exactly. Circular
defects were well defined in all cases; this was helped by the curvature control which
in fact provided a local mesh refinement in the case where a small hole was modelled.
179
Chapter 5
Accurate modelling of defects using Finite Elements
Overall, the author recommends using meshes of regular linear square elements for
models where features can be aligned with the node grid. Automatically generated
meshes made of modified quadratic triangular elements should be used in cases where
this can not be done or where features include curved edges.
The number of nodes that defines a small feature plays an important role in the quality
of the modelling of the interaction of a wave with the feature. In this situation, the
author would recommend using at least 8 nodes to define a feature in order to be
confident that the physical phenomenon is correctly represented.
With regards to the mesh density, it is difficult to provide a general rule. Nevertheless,
it can be said that a mesh density of 10 generally provides good qualitative results and
that a mesh density of 30 generally provides excellent qualitative and quantitative
results. It is therefore down to the modeller to decide what mesh density to use for a
particular model based on what is required from the model.
180
Chapter 6
Local mesh refinement
6.1. Introduction
As seen in the previous chapter, when using a mesh made of regular square elements,
it is challenging to accurately model defects whose geometry is not aligned with the
node grid. One technique targeted at correctly modelling the interaction of waves with
complex defects is known as the fictitious domain method. This technique is presented
in Section 6.2. Research is still active in this topic and has successfully shown to
improve modelling of complex defects but, to the best of our knowledge, the imple-
mentation of the technique cannot be achieved in commercially available FE packages
and requires the use of specialist FE codes.
The work in this chapter solely focuses on local mesh refinement of a regular square
mesh using conventional tools available in the commercial FE package ABAQUS. The
advantage of opting for a commercially available solution is that it removes the need
to create and maintain a specialist code and strongly facilitates the transfer of the
capability to industry. The method investigated in this part consists in creating two
separate areas having a different mesh density and tieing them together.
181
Chapter 6
Local mesh refinement
There are two main issues which have not allowed this to be implemented straightfor-
wardly in FE models in the past. The first one is the fact that changing the element size
in a regular mesh of elements changes the acoustic impedance of the mesh. This fact is
confirmed by the findings of chapter 4 as it was demonstrated that wave velocities
depend on the type of element and the density of the mesh. The other issue is related
to the way meshes of different size are tied together. At the interface between meshes
of different densities, the displacements of both meshes on each side of this boundary
are matched closely. As this is not done perfectly, it causes reflections.
In this chapter, a range of FE models are created and studied in order to better
understand and quantify the influence of these sources of error.
6.2.1 Review
A recent review by Dr Elizabeth Skelton [28] offers a good overview of the technique.
The development of the fictitious domain method started in the early 1960s but only
gained momentum in the 1990s through the work performed at the INRIA research
centre in France. Early developments and a portion of the research that followed
concentrated on the implementation of the technique for elliptic equations [75, 76, 77].
Since our focus is on wave propagation problems, it is noted that developments of the
technique for the Helmholtz problem were performed in several papers [78, 79, 80].
Time domain methods were also developed for time dependent acoustic [29] and
elastodynamic [30] problems.
6.2.2 Presentation
182
Chapter 6
Local mesh refinement
equation in order to solve the problem. The technique proved to work well for a variety
of cases and is still the subject of active research.
6.2.3 Conclusion
Despite the success of the technique, it requires the modification of the way the solver
computes the results. This unfortunately makes this technique not adapted to an imple-
mentation into a commercial FE package. Although the technique seems promising, as
the implementation in commercial FE packages is not possible, it is not investigated
further in this work.
In this section, the reflection from an abrupt change in mesh density is investigated.
As in the previous chapter, models are non-dimensional and the elastic material is
defined so that the longitudinal and shear wavelengths are 2 and 1 units respectively:
Young’s modulus is 8/3, Poisson’s ratio is 1/3 and density is 1. The general geometry
of the models is presented in Figure 6.1. Models are 2D plane strain, 100 units long and
1 unit thick but are constrained (i.e. displacement equal to zero) in order to have 1D
plane wave propagation. For 1D longitudinal wave propagation, the vertical
displacement is set to zero on the top and bottom surfaces of the model. In the case of
shear wave propagation, the horizontal displacement is constrained.
A regular square mesh is used over the whole model but it is split halfway (50 units
from excitation) into two parts. The density of the mesh is different between the left
and right hand sides. Both parts are tied together using the TIE function of ABAQUS
6.6-1 [14]. This function applies a rule that links the displacement of the nodes of the
boundaries of the two parts.
The excitation is applied using a uniform force of unit amplitude over the whole
thickness at the left hand extremity of the geometry. An ALID is used on the right hand
183
Chapter 6
Local mesh refinement
extremity of the model to absorb incoming waves. Figure 6.1 shows a schematic of the
model geometry.
Left part of model: Element size=0.1 Right part of model: Element size=0.2
{
{{
ALID
Constraint u=0 or v=0
Excitation
Models are run with ABAQUS 6.6-1 in both the time domain (explicit) and the
frequency domain (implicit). For the latter, the excitation is applied as a real harmonic
excitation whose frequency is 1; the complex displacement field is obtained directly
and used later to plot the absolute displacement field. For the explicit case, in order to
obtain an absolute displacement field similar to the frequency domain one, two
separate models are run: one with cosine and one with sine excitation of unit frequency.
The complex displacement field is obtained by adding the result of the cosine case to
the result of the sine case multiplied by i. The first half of a Hanning window is used
over the first 5 cycles in order to ramp up the amplitude of the sine and cosine to a
constant value in order to avoid exciting a broad range of frequencies. The model is run
sufficiently long so that the steady state is reached at all locations where the
displacements are monitored.
For longitudinal waves, the horizontal displacement is monitored at the half thickness
line from 10 to 40 units away from the excitation. For shear waves, the vertical
displacement is monitored. A spatial FFT [61] is applied to the complex displacement
field in a similar way to the example presented in Section 3.5.3.1. As shown in Figure
3.25, the positive peak gives the amplitude of the incident wave generated by the
excitation and the negative one the amplitude of the wave travelling in the opposite
direction caused by the interface between the two parts of the model. The ratio of these
two amplitudes gives the reflection coefficient of the waves from the interface. As
explained in Section 3.5.3.1, the propagation velocity is given by the x-position of the
184
Chapter 6
Local mesh refinement
peak on the spatial FFT graph. In this case, padding is applied to the FFT to have a
precise value of the velocity.
Implicit models are run using linear square elements with full (CPE4) and reduced
(CPE4R) integration. Explicit models are only run using reduced integration linear
square elements (CPE4R) as they are the only linear square elements available for this
type of solver. The tie function is used with its default settings for implicit
(TYPE=SURFACE TO SURFACE) and explicit (TYPE=NODE TO SURFACE)
solvers with no adjustment on the position of the nodes as the two tied surfaces
superimpose perfectly. Element based surfaces are used with explicit and implicit
models. In both cases, following the user manual instructions, the surface with coarsest
mesh is set as the master surface. Details about the TIE function can be found in [14].
In the first model, only longitudinal waves are generated. The element size is 0.1 on
the left part and 0.2 on the right part. This corresponds to a mesh density of 20 and 10
elements per longitudinal wavelength. The theoretical material properties given above
are used.
The results obtained from the model give a reflection of 1.972% (-34.1dB) of the
incident wave in all cases. The fact that the same reflection is obtained for all solvers
and element types means that the integration type (full or reduced) and solver type
(explicit or implicit) do not influence the amplitude of the reflection at the interface
between the two meshes. This also confirms that implicit and explicit implementation
of the tie are correct. The following models are therefore run just with the implicit
solver and reduced integration elements. The propagation velocities measured on each
side of the interface are 1.992 and 1.966. This agrees very well with results obtained in
Chapter 4, as the predicted velocities at 0 degrees are 1.991 and 1.964 for a regular
mesh of 20 and 10 elements per longitudinal wavelength respectively.
185
Chapter 6
Local mesh refinement
As presented in [34], the reflection at 0 degrees based on the impedance is given by:
Z 20 – Z 10 ρ 20 c 20 – ρ 10 c 10
RC = ---------------------- = ------------------------------------- (6·1)
Z 20 + Z 10 ρ 20 c 20 + ρ 10 c 10
where Z is the impedance, ρ the density and c the velocity. The subscript indicates the
density of the mesh.
Using the velocities measured in the model, the theoretically predicted reflection due
to the change in impedance is 0.657%. As the amplitude of the reflection is higher than
this value, the remaining part of the reflection may be due to the way the interface is
implemented.
The second model is similar to the first model but in order to separate the contribution
of the velocity error and the tie of the mesh, the acoustic impedances of both meshes
at 0 degrees are matched by adjusting the longitudinal velocities in both parts of the FE
model to be equal to the theoretical one (cLth=2) The velocity in the original FE model
cLor is defined by the theoretical velocity cLth multiplied by a ratio RL:
1 1 E th ( 1 – ν th )
c Lor = ------ c Lth = ------ ------- ---------------------------------------------- (6·2)
RL R L ρ th ( 1 + ν th ) ( 1 – 2ν th )
c Lth
with R L = --------- (6·3)
c Lor
The original velocity in the FE model cLor necessary to define RL is either measured
by spatial FFT in the model or derived by using the estimate given in Chapter 5.
As the error between cLor and cLth is small, the adjusted velocity cLadj is defined as:
E th ( 1 – ν th ) E adj ( 1 – ν adj )
c Ladj = R L c Lth = R L ------- ---------------------------------------------- = ---------- --------------------------------------------------
- (6·4)
ρ th ( 1 + ν th ) ( 1 – 2ν th ) ρ adj ( 1 + ν adj ) ( 1 – 2ν adj )
186
Chapter 6
Local mesh refinement
The density and Poisson’s ratio are kept the same but the Young’s moduli are changed:
2
ν adj = ν th , ρ adj = ρ th and E adj = R L E th (6·5)
The adjusted Young’s moduli for mesh densities of 20 and 10 are 2.765 and 2.688. The
adjusted velocities measured with the spatial FFT are 2.000 and 2.000 for densities of
10 and 20 elements. Therefore the theoretical reflection due to the change of
impedance should be 0%. The reflection coefficient given by the spatial FFT is 1.240%
(-38.134dB). As expected, matching the impedance has reduced the reflection
coefficient. Interestingly, the difference in the reflection coefficient of the first and
second models is 0.732% which is close to the predicted contribution of the change of
impedance in the first model (0.657%). It indicates that the total reflection from the
interface is equal to the sum of the change of impedance and the error due to the tie of
the 2 meshes.
In this model, both shear and longitudinal waves are modelled in separate runs. The
element size is 0.05 and 0.1 on the left and right part of the model respectively. The
mesh density is therefore 20 and 10 elements per shear wavelength and 40 and 20
elements per longitudinal wavelength.
Theoretical material properties are used to start with. The measured reflection
coefficients in this case are 0.470% (-46.56dB) and 1.971% (-34.108dB) for
longitudinal and shear waves. The measured velocities are: cLleft=1.998, cSleft=0.996,
cLright=1.9918 and cSright=0.983. The reflection caused by the change of impedance is
predicted to be 0.155% and 0.657% of the incident signal. The reflection coefficient
for shear waves in this case is the same as the longitudinal one in the first 1D model.
The number of elements per considered wavelength is the same in both cases. This
leads to the error on velocities being the same and therefore the change of impedance
being the same. It indicates that both shear and longitudinal waves interact in the same
way with the interface at equivalent number of elements per wavelength.
187
Chapter 6
Local mesh refinement
The material properties are adjusted in order to match the longitudinal wave acoustic
impedance on both sides of the model. The velocities then are: cLleft=2.000,
cSleft=0.997, cLright=2.000 and cSright=0.983. The reflection caused by the change of
impedance is predicted to be 0% and 0.484% of the incident signal. The measured
reflection coefficient for the longitudinal wave is 0.309% (-50.19dB) and for the shear
wave is 1.795% (-34.92dB). In both cases, the drop in reflection coefficient matches
the change in the contribution of the change of impedance. This confirms the findings
of the previous model.
Following this, the material properties are adjusted in order to match the shear wave
acoustic impedance in both sides of the model. The velocities then are: cLleft=2.006,
cSleft=1.000, cLright=2.026 and cSright=1.001. In this case, the velocities in the right
hand side of the model are higher than in the left one. This leads to a negative reflection
contribution due to the impedance change which indicates a phase shift of 180 degrees.
It is predicted to be -0.496% and -0.050% of the incident wave. The reflection
coefficient for the shear wave is 1.240% and for the longitudinal wave is -0.197%. The
negative sign indicates a phase shift. The results for the shear wave are in line with the
previous cases. For the longitudinal wave, it is very interesting to note that the
reflection coefficient (-0.197%) is also very close to be the sum of the contribution due
to the change of impedance (-0.496%) and the interface (0.309%). This indicates that
it should be possible to completely remove the reflection by matching the reflection
due to the interface with a reflection due to an impedance change of equal amplitude
and opposite phase.
Finally, the material properties are adjusted in order to match both the shear and
longitudinal wave acoustic impedance in both sides of the model.
188
Chapter 6
Local mesh refinement
It is necessary to adjust the Young’s modulus and the Poisson’s ratio to do this.
1 1 E th ( 1 – ν th )
c Lor = ------ c Lth = ------ ------- ---------------------------------------------- (6·6)
RL R L ρ th ( 1 + ν th ) ( 1 – 2ν th )
c Lth
with R L = --------- (6·7)
c Lor
1 1 E th 1
c Sor = ------ c Sth = ------ ------- ------------------------ (6·8)
RS R S ρ th 2 ( 1 + ν th )
c Sth
with R S = --------- (6·9)
c Sor
As the error between the FE and theoretical velocities is small, the adjusted velocities
cLadj and cSadj are defined as:
E th ( 1 – ν th ) E adj ( 1 – ν adj )
c Ladj = R L c Lth = R L ------- ---------------------------------------------- = ---------- --------------------------------------------------
- (6·10)
ρ th ( 1 + ν th ) ( 1 – 2ν th ) ρ adj ( 1 + ν adj ) ( 1 – 2ν adj )
E th 1 E adj 1
c Sadj = R S c Sth = R S ------- ------------------------ = ---------- --------------------------- (6·11)
ρ th 2 ( 1 + ν th ) ρ adj 2 ( 1 + ν adj )
The density is kept equal to 1 but the Young’s moduli and Poisson’s ratios are changed.
Replacing cLadj and cSadj in 6·12 with 6·10 and 6·11 and doing some algebra, we have:
2 2
2R L – R S 2 2 ( 1 + ν adj )
ν adj = --------------------- , ρ adj = ρ th and E adj = R S --------------------------- (6·13)
2
4R L – R S
2 ρ adj
189
Chapter 6
Local mesh refinement
of the incident signal with longitudinal and shear waves respectively. The measured
reflection coefficient for the longitudinal wave is 0.309% (-50.19dB) and for the shear
wave is 1.240% (-38.13dB). In both cases, the drop in reflection coefficient matches
the change in the contribution of the change of impedance compared to the non
adjusted case. In this situation, the reflection is only due to the interface between the 2
parts.
This model is similar to the previous one except that in this case, the Young’s modulus
and Poisson’s ratio are kept at their theoretical values in the left part of the model
where the element size is 0.05 whereas the ones for the right part where the element
size is 0.1 are varied parametrically.
To start with, the Poisson’s ratio is kept constant and the Young’s modulus is varied
from 2.65 to 2.8 for the longitudinal wave and 2.65 to 2.95 for the shear wave. The
reflection coefficient is measured as the Young’s modulus is varied. The results are
plotted in Figure 6.2. The reflection coefficient due to the impedance mismatch
between the 2 parts of the model is calculated using Equation 6·1 and is plotted against
the Young’s modulus of the right part of the model. The difference between this and
the total reflection is also plotted.
The first point of interest is that with both wave types the variation of the reflection
from the FE model is linear against the Young’s modulus. The impedance part of the
reflection also varies linearly and with a slope close to the one of the total reflection.
This is confirmed by the fact that the difference between these 2 curves, although not
constant, varies linearly with a low slope. As mentioned in the previous part, it can be
considered that the total reflection is made of the sum of the contribution due to the
impedance mismatch and the one from the tie. In the previous plot, the difference
between the FE and impedance related reflection is plotted as the reflection due to the
tie. The assumption is made that the calculated and actual reflection due to the
impedance mismatch are equal. As it is likely that a small discrepancy between these
2 exists, the plot for the tie reflection therefore includes this source of error. It is
190
Chapter 6
Local mesh refinement
a) 0.8
0.4
-0.4
-0.8
-1.2
2.65 El 2.7 2.75 2.8
Young’s modulus
b)
2
S wave reflection (%)
-1
-2
2.65 El 2.7 2.75 2.8 2.85 2.9 2.95
Young’s modulus
FE Impedance FE-Impedance (Tie?)
Figure 6.2. Reflection coefficient for longitudinal and shear waves against Young’s modulus.
possible that this discrepancy is the reason for the reflection due to the tie not being
constant but given the results from the previous model, this is unlikely. It can be said
with a fair degree of confidence that the variation in the reflection due to the tie is
mainly due to the variation of the Young’s modulus rather than a change in the
difference between the calculated and actual reflection due to the impedance
mismatch.
The variation in Young’s modulus leads to a variation in the overall reflection which
is largely driven by the change in the impedance. It can be seen that the reflection due
to the interface is higher for the shear wave case than the longitudinal one. Using
results from 6.3.1.2, it is clear that this is not related to the wave type but to the mesh
density and indicates that the interface reflection varies with the mesh density and that
191
Chapter 6
Local mesh refinement
it is dependent on the number of elements per wavelength of the considered wave. This
point is investigated in detail further in this chapter.
Following this, the Poisson’s ratio is varied from 0.32 to 0.335 and the Young’s
modulus from 2.6 to 2.9. The reflection coefficient is measured as the Young’s
modulus and Poisson’s ratio are varied. The results are plotted in Figure 6.3. The
reflection coefficient due to the impedance mismatch between the 2 parts of the model
is calculated using equation 6·1 and is plotted against the Young’s modulus in the right
part of the model. The difference between this and the total reflection is also plotted in
Figure 6.3.
a) 2
L wave reflection (%)
-1
-2
0.32
Po
0.325
iss
on
0.33
’s
r
ati
0.335 2.9
o
-1
0.32
Po
0.325
iss
on
0.33
’s
r
ati
Figure 6.3. Reflection coefficient for longitudinal and shear waves against Young’s modulus and
Poisson’s ratio.
The variation of the reflection coefficient against the Young’s modulus and Poisson’s
ratio is close to linear for both parameters and wave types. The reflection due to the
192
Chapter 6
Local mesh refinement
impedance changes varies in a similar fashion. This means that, as can be seen on the
figure, the reflection due to the interface is relatively stable over the whole range of
parameters. As seen in the previous case, the amplitude of the reflection due to the
interface is not the same for both wave types and depends on the mesh density.
Given the shape of the total reflection surfaces, it should be possible to find a point
where there is no reflection for both wave types. Making the assumption that the
surfaces are flat ones, the equations for the reflection coefficient for longitudinal and
shear wave are found to be:
The total reflection is therefore zero for both waves when the elastic properties in the
right part of the model are: Young’s modulus=2.8372 and Poisson’s ratio=0.3208.
A FE model is run for these parameters. The reflection is 0.006% (-84.52dB) and
0.005% (-85.75dB) for longitudinal and shear waves respectively. These values are
negligible and confirm that it is possible to remove the reflection from both wave types
at a given angle by adjusting the Poisson’s ratio and Young’s modulus.
As seen previously, the reflection due to the interface is close to constant when the
acoustic impedance in one side of the model is changed but it varies noticeably with
the mesh density. In this model, only longitudinal waves are studied as it was
demonstrated that the reflection from the interface is independent of the wave type.
The theoretical material properties as defined in Section 6.3.1 are used. The mesh
density of both parts of the model is varied from 6 to 40 elements per longitudinal
wavelength. The reflection coefficient is measured over this range of values and is
plotted in Figure 6.4.a. The reflection due to the change in impedance calculated by
using the measured velocities in each part of the model with equation 6·1 is plotted in
Figure 6.4.b. The difference between these 2 reflections is considered to be the
reflection due to the tie. This is plotted on a linear and log scale in Figure 6.4.c and d.
193
Chapter 6
Local mesh refinement
a)
3
2
1
0
-1
-2
-3
5
10
15 35 40
20 30
25
30 15 20 25
Nright 35 10 Nleft
40 5
c)
6
(Total - impedance)
Tie reflection (%)
4
2
0
-2
-4
-6
5
10 40
15 35
20 25 30
25 20
Nright 30 35 10 15 Nleft
40 5
d)
-20
-30
(Total - impedance)
Tie reflection dB)
-40
-50
-60
-70
-80
-90
-100 15 10 5
5 20
10 15 20 25 35 30 25
30 35 40 Nleft
Nright 40
Figure 6.4. a) Total reflection, b) Reflection due to the impedance change, c) and d) Reflection due to
the tie (linear scale and log scale).
194
Chapter 6
Local mesh refinement
The total reflection varies smoothly as the mesh densities are varied. As expected,
central to all graphics in Figure 6.4 is a line of symmetry of zero reflection. Figure 6.4.d
highlights that the reflection due to the tie is largely dependent on the mesh densities
of both parts of the model. The reflection coefficient due to the tie and the impedance
for a range of combinations of meshes in the situation where material properties are the
same are presented in Table 6.1 and 6.2. It is interesting to note that the reflection
caused by the tie contributes about twice as much as the one due to the impedance.
Table 6.1.Table of reflection coefficients due to the tie between two meshes in %
N 15 20 25 30 35 40
10 0.995% 1.325% 1.478% 1.561% 1.608% 1.639%
15 0.329% 0.482% 0.565% 0.612% 0.643%
20 0.152% 0.235% 0.281% 0.312%
25 0.082% 0.128% 0.159%
30 0.045% 0.076%
35 0.030%
Table 6.2.Table of reflection coefficients due to the impedance difference between two
meshes in %
N 15 20 25 30 35 40
10 0.481% 0.646% 0.720% 0.759% 0.786% 0.802%
15 0.165% 0.238% 0.278% 0.304% 0.321%
20 0.073% 0.113% 0.1393% 0.156%
25 0.039% 0.066% 0.082%
30 0.026% 0.043%
35 0.016%
As seen in the previous part, the reflection from the interface between 2 meshes can be
made equal to zero at normal incidence. In this part, the effect of removing the
reflection at one angle is investigated over the full angle range of 2D models.
The 2D model is similar to the models used for the study of Chapter 4. At the centre of
the model is a 70 unit square made of a fine mesh of density 20 elements per shortest
wavelength. This square is bordered by a 9 unit thick layer made of a coarser mesh of
195
Chapter 6
Local mesh refinement
density 10 elements per shortest wavelength which is bordered by a 6 unit thick ALID
of the same mesh density. This is shown in Figure 6.5.
ALID
Coarse mesh
Fine mesh
Excitation
point
Both shear and longitudinal wave generations are considered. Models are excited by
body forces being applied on four elements as in Chapter 5. This generates an almost
pure wave - shear or longitudinal - and, as the size of the source is only a fraction of a
wavelength, is close to a point source. The point located at the centre of the four
elements is referred to as the excitation point. This point is located at the centre of the
model. A 10 cycle Hanning windowed tone burst of centre frequency 1Hz is used.
The model is run for two cases. In one case, it is run using the reference material
properties defined in Section 6.3.1.1 in both the coarse and fine mesh parts. In the other
case, the material properties are the theoretical ones in the fine mesh part and the
adjusted ones in the coarse mesh part in order to have zero reflection for both wave
types at 0 degrees as defined in Section 6.3.1.5. The displacement fields in the top right
corner at t=30 and t=60 for longitudinal and shear waves respectively are plotted in
Figure 6.6.
Unlike in the 1D model, the waves decay as they propagate, due to beam spreading.
The colour scale is adjusted so that the colour range covers the range of amplitude from
196
Chapter 6
Local mesh refinement
a) ALID
Coarse mesh
Fine mesh
c)
b)
e)
Figure 6.6. Absolute displacement field in the top right corner of the 2D models. Longitudinal wave
excitation with a) theoretical and b) adjusted material properties. c) Definition of wave packet
positions. d), e), f) same with shear wave excitation
197
Chapter 6
Local mesh refinement
0 to 1% of the desired incident wave when it hits the interface between coarse and fine
meshes. All data above 1% is coloured in grey. The position and propagation direction
of the different wave packets are schematically represented in order to simplify under-
standing.
Looking at the results for longitudinal wave excitation, it can be seen that there is a
noticeable longitudinal wave reflection occurring at 0 degrees in the case using the
theoretical material properties everywhere. The amplitude of this wave at 0 degrees is
about 0.5%. It confirms that in the case of an abrupt change in mesh density the implicit
and explicit cases give similar results. The use of adjusted properties significantly
reduces the amplitude of the reflection for low angles of incidence but the amplitude
of this reflection seems to increase with the angle of incidence. The shear wave
reflection is reduced by the change in properties but has a similar amplitude variation.
Shear wave results carry similar findings. The shear wave reflection is not seen at 0
degrees but can be clearly identified at higher angles of incidence.
It can be concluded that the use of adjusted material properties based on the 1D model
does not provide a significant gain in 2D compared to the non adjusted case. The
explanation for this is two fold. First, as the angle of incidence varies, the visible
acoustic impedance of the meshes changes because the mesh density is altered.
Therefore the exact opposition of the reflections due to the tie and the impedance
change does not exist anymore. Secondly, as mode conversion exists for angles
different to zero, the reflection characteristics at each angle changes. Unique
combinations of waves exist for each angle of incidence and means that there is a
strong increase in the complexity of the problem.
Minimising the reflection due to the change of impedance over the full range of angle
of incidence is an achievable target. Significantly reducing the reflection due to the tie
below its original value is unlikely to be easily achievable if at all possible.
198
Chapter 6
Local mesh refinement
As discussed previously, in a model where 2 meshes using the same material properties
but having different mesh densities are tied together, a reflection occurs at the
interface; this reflection is caused by the change of impedance and the error in the tying
of the meshes. This latter contributes twice as much as the former. It would be
beneficial if the reflection due to the tie could be reduced significantly. A solution
which is investigated in this part is to have a gradual change in mesh density rather than
an abrupt one.
The 1D model is similar to the 1D models used in the previous part in all aspects apart
from the way the mesh density change is achieved. In this case, the change in mesh
density is achieved over a number of layers so that the mesh density is varied by 1 unit
per layer. The layers are made of square elements whose dimension vary as the mesh
density is varied. This is shown in Figure 6.7.
Left part of model: Element size=0.1 Right part of model: Element size=0.2
{
Constraint
{
{{
Gradual interface ALID
Excitation
The effect of the impedance change is removed by matching the acoustic impedance
for both wave types at 0 degrees in every part of the model as described in Section 6·13.
The impedances are adjusted to 2 for longitudinal waves and 1 for shear waves in all
parts of the model. It is run in the time domain with the explicit solver. The excitation
consists of a 10 cycles tone burst of unit centre frequency. The thickness is 1 unit. The
left and right parts are 75 and 100 units long respectively. The central part where the
gradual change of mesh density occurs is 0.62 units long.
199
Chapter 6
Local mesh refinement
The absolute displacement fields obtained at times 64 and 128 for longitudinal and
shear waves are plotted in Figure 6.8. At these times, the incident wave packets have
travelled through the interface. As previously, the colour scale is adjusted so that the
colour range covers the amplitude from 0 to 1% of the incident wave. This is done to
highlight low amplitude signals. All data above 1% is coloured in grey.
{ {
interface
{
b)
Reflection Gradual Incident wave packet
{
interface
{
X scale: 10 units =1.18cm Y scale: 1 units =1.4cm
Figure 6.8. Absolute displacement field for a) longitudinal and b) shear wave excitation models with a
gradual change of mesh density at t=34 (longitudinal) and t=68 (shear).
In both cases, the incident wave packet goes through the interface with almost no
change to its amplitude but it can be seen that reflection from the interface occurs in
both cases. Their amplitude is 0.52% and 0.32% of the longitudinal and shear incident
wave packet. In the longitudinal wave case, it can be noted that the reflection is more
spatially spread than the incident signal. This does not appear to be the same in the
shear wave case. After the interface, the longitudinal wave packet is followed by a trail
of waves similar to the reflection one. The striking feature in both cases is the presence
of unexpectedly high displacements localized around the interface. For the fields
plotted above, the maximum displacement is 10.5% and 0.65% of the incident
longitudinal and shear wave packets. Given that boundary conditions are set such that
the model should behave as a 1D one, the fact that these displacements are not constant
in the vertical direction indicates numerical problems caused by the gradual interface.
This shows that the use of a gradual change of mesh density does not significantly
reduce the reflection from the interface. It also causes numerical issues in the model.
200
Chapter 6
Local mesh refinement
The 2D model is similar to the 2D models used in Section 6.3.2 apart from the way the
mesh density change is treated. In this case, the mesh density change occurs over a
range of layers as in the 1D study of this section. Figure 6.9 shows the way this is
achieved in the 2D model. The mesh density is varied by 1 unit for each layer to
gradually change the mesh density from 10 to 20. The Young’s modulus and Poisson’s
ratio are adjusted to minimize the reflection due to the impedance change at 0 degrees.
Coarse mesh
{
Transition zone
{
Fine mesh
{
Figure 6.9. Gradual mesh density change for the 2D model.
The displacement fields in the top right corner of the model at t=35 and t=70 for
longitudinal and shear wave respectively are plotted in Figure 6.10. The colour scale
is adjusted so that the colour range covers the amplitude from 0 to 1% of the incident
wave when it hits the gradual interface. All data above 1% is coloured in grey.
a) Inc b) Inc
ide ide
nt nt
w wa
a
ve
ve
p
p
ac
ac
ke
ke
t
t
Figure 6.10. Absolute displacement field in the top right corner of the 2D model with a) theoretical and
b) adjusted material properties.
201
Chapter 6
Local mesh refinement
As for the 1D model, it can be seen that gradual mesh density change as implemented
in this case leads to a deterioration of the displacement field. This confirms that gradual
change of mesh density with ties is not a viable option to achieve local mesh
refinement.
6.5. Conclusions
The study of 1D models with abrupt change in density shows that the reflection is due
to the change in impedance and the tie of the meshes. It is shown that the reflection due
to the impedance change can be predicted using results from Chapter 5. This part of the
reflection can be adjusted so that the overall reflection is zero for both wave types in
the 1D model. In such case, it has the same amplitude but opposite phase as the part
due to the tie of the meshes. 2D models show that this adjustment of the impedance
leads to a significant reduction in the reflection at 0 and 90 degrees of incidence but
that it is still noticeable at other angles of propagation different to the ones simulated
in the 1D model (0 and 90 degrees). It is unlikely that it is possible to remove the
reflection for all angles using this method. The contribution of the impedance change
and tie at 0 degrees for a range of combination of mesh density in cases where
theoretical material properties are used are presented in Tables 6.1 and 6.2. It gives
modellers an estimate of the amplitude of reflection to expect if the technique is used
and shows that the amplitude, although not negligible, is limited.
As the impedance reflection can be made equal to zero at a given angle and as the tie
is the main contributor to the reflection, the idea of reducing the reflection due to the
tie by using a gradual change in mesh density is investigated. The gradual change in
mesh density is achieved by using a series of layers of varying mesh densities between
the 2 meshes. Results from the 1D and 2D models show that this technique does not
reduce the reflection significantly and causes localized numerical problems close to the
interface although there is room for more detailed parametric studies.
Overall, local mesh refinement of a square mesh as implemented in this chapter with
an abrupt change in mesh density has shown to work correctly but generates reflection
at the interface. This can be an acceptable configuration in some cases. In other cases,
202
Chapter 6
Local mesh refinement
the use of an automatically generated modified quadratic triangular mesh (which does
not use ties) is likely to be the best option. This is discussed in Chapter 7.
203
Chapter 7
Conclusions
The work in this thesis investigates ways of improving elastic wave modelling in
isotropic solids using commercially available FE packages. The focus is on developing
new techniques as well as gaining a better understanding of the existing techniques.
Knowledge of the interaction of waves with defects and structural features is crucial in
order to develop satisfactory ultrasonic NDT applications. As applications are
developed for increasingly challenging situations, the use of numerical modelling has
proved in many cases to be the most effective and cost efficient way of assessing
potential problems early in the development process. Despite the exponential increase
in computational capacities in recent years, there are many cases which remain out of
reach.
The choice of commercially available FE packages is motivated by the fact that they
enable new techniques to be quickly transferred to industry where these products are
already or could easily be available and where a specific unsupported research code
would be difficult to maintain.
Chapter 2 presents the theoretical basis of the bulk and guided wave propagation
phenomena modelled in this work. Details of the numerical scheme used in this work
are provided.
The implementation of the radiation of waves outside the area of study is studied in
Chapter 3. Following the assessment of the available techniques, 2 types of absorbing
layers are selected: Perfectly Matched Layers (PML) and Absorbing Layers Using
Increasing Damping (ALID). The implementation of both methods is discussed. It is
shown that the correct definition of layers parameters is crucial to the achievement of
numerical efficiency. Analytical models are developed in order to facilitate the
determination of these parameters for bulk wave and 2D plate guided wave problems.
Analytical model results are validated against FE validation cases results.
204
Chapter 7
Conclusions
Demonstrator cases for bulk and guided waves are developed to demonstrate gains
achieved by the use of the techniques. The procedure of time domain reconstruction of
frequency domain results is explained and validated. The ease of use of the frequency
domain method to obtain reflection coefficient curves is presented and results are
compared with published results.
The use of local mesh refinement techniques in regular square meshes is investigated
in Chapter 6. One method investigated consists in tying together 2 meshes of different
mesh densities. 1D and 2D models are used to evaluate the technique. Results obtained
in Chapter 4 are used to improve the performance of the technique. The second method
investigated achieves the mesh density change gradually through the use of a series of
tied layers. The performance of this method is measured using 1D and 2D models.
The study of PML and ALID in wave propagation models has shown that a significant
reduction in a model’s geometric size can be achieved through their use. Analytical
205
Chapter 7
Conclusions
models provide a fast and efficient way of defining the layer’s parameters and therefore
offer a significant improvement over trial and error methods. Demonstrator cases have
highlighted the gains obtained by using ALID or PML and have also shown that PML
achieves smaller model sizes than ALID. This can be explained by the fact that with
PML the layer matches the impedance of the area of study perfectly whereas for the
ALID the strength of the decay is limited by the impedance change it causes. As a
consequence, for a given case, ALID is generally thicker than PML. It is also more
sensitive to waves incident at high angles, leading to a necessary increase of the area
of study in order to minimize the model size. As the aim is to stretch modelling
capabilities, this gives an advantage to the PML technique but it has to be balanced by
the fact that, to our knowledge at the time of writing, COMSOL is the only mainstream
FE package that allows implementation of PML. The implementation is limited to
frequency domain solving with implicit solvers whereas ALID can be implemented to
be used with any FE package which allows Rayleigh damping. The actual performance
comparison depends on the type of result required.
A user wishing to obtain frequency domain results such as the reflection coefficient
from a defect has a choice of directly using either PML or ALID. The two best
combinations investigated are ABAQUS with ALID and COMSOL with PML.
ABAQUS models were solved using the “Direct-solution steady-state dynamic
analysis” procedure of ABAQUS/Standard which uses a direct implicit solver.
COMSOL offers a range of direct and iterative implicit solvers. Direct solvers have
proved to be robust and efficient but limited in term of the number of degrees of
freedom that can be solved. Iterative solvers showed improvement in terms of speed
and memory but convergence proved difficult. Overall, in our experience, using
COMSOL 3.2 and ABAQUS 6.6, COMSOL did not match the performance of
ABAQUS for a given number of degrees of freedom. Given these two combinations it
is therefore not possible to clearly recommend one of them, as the best solution
depends on the specific details of the model considered.
For time domain results, the most straightforward approach is to use ABAQUS/
Explicit with ALID. The explicit solver of ABAQUS uses the central difference
algorithm. It is particularly attractive for large models as it is faster and more memory
efficient than an implicit solver. This technique is recommended for most cases where
time domain results are required. One alternative is to solve the model in the frequency
206
Chapter 7
Conclusions
domain and reconstruct the time domain results by performing an inverse Fourier
transform.
The study of reduced integration linear quadrilateral elements solved with an explicit
solver showed the influence of the scaled Courant number CFLX and mesh density N
on the velocity. Patterns were identified. These permitted rules to be defined to
evaluate with high accuracy the amplitude of the velocity error at any angle, any scaled
Courant number and for a mesh density higher than 6. Following this, the study of
deformed elements highlighted that the maximum error in a mesh could be predicted
using the rules defined for square elements. This indicates that this element type can
cope well with being deformed.
With linear equilateral triangular elements, it was demonstrated that the error could be
correctly predicted using rules established for quadrilateral elements along the
altitudes and element sides for longitudinal waves and along the element sides for shear
waves. The influence of the scaled Courant number CFLX and mesh density N was
confirmed at these positions and wave types. A zero error was measured for shear
waves along the altitudes of the triangular elements. Deformed mesh cases showed that
this element type is sensitive to high levels of deformation as the predicted error was
more than tripled in one case (16.4% instead of 4.8%).
207
Chapter 7
Conclusions
quadratic quadrilateral elements proved to strongly reduce the velocity error. Rules
established in previous cases could not be applied to this case but given that the
velocity error is less than 0.5% for a mesh density higher than 8, it can be understood
that the velocity error with this type of element is not a cause for concern. Deformation
of elements did not prove to be a problem.
Linear and modified quadratic elements results for the implicit case matched almost
exactly results obtained with the explicit solver. Quadratic equilateral triangular
elements were shown to lead to low velocity error as were quadratic square elements,
but proved to be sensitive to high levels of deformation in the same way as linear
triangular elements.
Overall, the influence of the solver is ruled out as element type is the factor
determining the velocity error. Unlike quadrilateral elements, all triangular elements
apart from modified quadratic ones were shown to be sensitive to element deformation.
This justifies why automatic meshing algorithms seek to minimize the deformation and
why modellers need to exercise caution when using triangular elements.
The angle of propagation in a regular mesh was shown to play an important part in the
amplitude of the velocity error as it leads to a non-negligible variation.
The mesh density was shown to be the driving factor influencing the velocity error for
all element types. The influence is not so remarkable for quadratic elements as the
velocity is highly accurate for any mesh density. For linear and modified quadratic
elements, its influence is much stronger.
This study has highlighted the challenges of modelling a range of defects. The use of
a regular mesh was shown to perform very well to model a straight edge when aligned
with the mesh, but its limitations were demonstrated when modelling features whose
geometry does not follow the grid exactly. The issue of correctly dimensioning a
feature in a fixed grid was highlighted and can be easily resolved by adjusting the mesh
density to fit the size of the defect. It was demonstrated that this is not the main issue
influencing the quality of results. The stair-cased definition of the edges of defects
proved to be the crucial parameter influencing the quality of results. Even with large
defects, relatively fine meshes (N=25) performed poorly in cases dominated by shear
208
Chapter 7
Conclusions
waves. Therefore, the regular mesh approach can only be used with extremely fine
meshes (N>=30).
Triangular meshes generated by the free meshing algorithm have proved to perform
well in most cases. The fact that the element size is not exactly constant did not cause
a deterioration of results. Dimensioning of the defect was better than with the regular
square mesh approach. Straight defects were defined exactly. Circular defects were
well defined in all cases; this was helped by the curvature control which in fact
provided a local mesh refinement in the case where a small hole was modelled.
Quadratic and modified quadratic elements proved to perform better in situations
where only a limited number of nodes defined a defect.
The study of 1D models with an abrupt change in density shows that the reflection is
due to the change in impedance and the tie of the meshes. It is shown that the reflection
due to the impedance change can be predicted using results from Chapter 4. This part
of the reflection can be adjusted so that the overall reflection is zero for both wave
types in the 1D model. In such case, it has the same amplitude but opposite phase as
the part due to the tie of the meshes. 2D models show that this adjustment of the
impedance leads to a significant reduction in the reflection at 0 and 90 degrees of
incidence but that it is still noticeable at other angles of propagation different to the
ones simulated in the 1D model (0 and 90 degrees). It is unlikely that it is possible to
remove the reflection for all angles using this method. The contribution of the
impedance change and tie at 0 degrees for a range of combination of mesh densities in
cases where theoretical material properties are used are presented for reference in
Chapter 6. It gives modellers an estimate of the amplitude of reflection to expect if the
technique is used and shows that the amplitude, although not negligible, is limited.
As the impedance reflection can be made equal to zero at a given angle and as the tie
is the main contributor to the reflection, the idea of reducing the reflection due to the
tie by using a gradual change in mesh density is investigated. The gradual change in
mesh density is achieved by using a series of layers of varying mesh densities between
the two meshes. Results from the 1D and 2D models show that this technique does not
reduce the reflection significantly and causes numerical problems localized close to the
interface.
209
Chapter 7
Conclusions
Overall, local mesh refinement of a square mesh as implemented in this chapter with
an abrupt change in mesh density has shown to work correctly but generates reflection
at the interface. This can be an acceptable configuration in some cases. In other cases,
the use of an automatically generated modified quadratic triangular mesh (without ties)
is likely to be the best option.
It is mentioned that mass and stiffness proportional damping can be used in ALID. As
stiffness proportional damping can lead to a dramatic reduction in the time increment
in an explicit scheme, only mass proportional damping is used. For the implicit
scheme, this is not the case and the simultaneous use of both damping types may lead
to increased performance and could be investigated.
The bulk wave analytical models are valid for all cases (up to 3D) but the guided wave
model is only valid for 2D plates. It would be beneficial to develop similar models for
other guiding structures. In particular, in the field of Structural Health Monitoring
(SHM), the development of a model for 3D plates would be of particular interest. In a
similar way, a pipe analytical model would enable the creation of absorbing layers for
this type of system.
Models using PML with viscoelastic materials will need to be investigated as PML
proved to perform poorly when evanescent waves are present.
The work presented in this thesis is limited to a range of 2D plane strain elements. The
study of other element types would increase the scope of the current study. It would be
of particular interest to perform 3D elements studies.
As there are subtle differences in the way packages are implemented, small differences
may exist between packages. Reproducing some of the models presented in this thesis
with other FE packages should show that the results presented here can be used for
other packages. It should be noted that despite this, it is expected that the principles
developed here apply to all programs.
210
Chapter 7
Conclusions
Another aspect that would require additional work is the investigation of irregular
triangular element meshes such as those generated by automatic algorithms. Regarding
this, it would be beneficial to evaluate whether the irregularity of the mesh causes some
noise. If this occurs, the level of noise should be quantified.
The local mesh refinement technique of regular quadrilateral meshes presented in this
work did not lead to the achievement of extremely low amplitude reflection at the
interface and it seems unlikely that it would be possible to achieve it. The use of
automatically generated meshes made of triangular elements seems to be the most
promising technique for the achievement of local mesh refinement. Indeed, based on
results obtained in Chapter 5, no numerical problem was noticed when mesh
refinement in a triangular mesh was applied by the free meshing algorithm to define
small circular defects. As it is expected that numerical problems will occur in extreme
cases of mesh refinement, it would be interesting to determine the limits of this
particular refinement technique.
211
References
214
12. O. Diligent, T. Grahn, A. Bostrom, P. Cawley and M. J. S. Lowe, The low-
frequency reflection and scattering of the S0 Lamb mode from a circular through-
thickness hole in a plate: finite element, analytical and experimental studies,
J. Acoust. Soc. Am., vol. 112, pp. 2589-2601 (2001).
13. M. J. S. Lowe and O. Diligent, Low-frequency reflection characteristics of the
S0 Lamb wave froma rectangular notch in a plate, J. Acoust. Soc. Am., vol. 111,
pp. 64-74 (2002).
14. ABAQUS v6.6 Analysis User's Manual. (2006), www.simulia.com.
15. COMSOL User's Guide version 3.2. (2006), www.femlab.com.
16. M.Drozdz, M.J.S. Lowe, E. Skelton, R.V. Craster, Modeling bulk and guided
wave propagation in unbounded elastic media using absorbing layers in
commercial FE packages, Review of Progress in Quantitative NDE (2007), Vol.
26, pp. 87 - 94.
17. D. Givoli, High-order local non-reflecting boundary conditions: a review, Wave
Motion, vol. 39, pp. 319-326 (2004).
18. P. Bettess, Infinite Elements, Penshaw Press, 1992.
19. A.J. Burton and G.F. Miller, The application of integral equation methods to
the numerical solution of some exterior boundary-value problems, Proc. R. Soc.
Lond. A, vol. 323, pp. 201-210 (1971).
20. D. Givoli and J.B. Keller, Non-reflecting boundary conditions for elastic waves,
Wave Motion, 12, 261-279 (1990).
21. G.R. Liu and Jerry S.S. Quek, A Non-Reflecting Boundary for Analysing Wave
Propagation Using the Finite Element Method, Finite Elements in Analysis and
Design, vol. 39, pp. 403-417 (2003).
22. M. Israeli and S.A. Orszag, Approximation of radiation boundary conditions,
Journal of Computational Physics, vol. 41, pp. 115-135 (1981).
23. J. Sochacki, R. Kubichek, J. George, W.R. Fletcher and S. Smithson,
Absorbing boundary condition and surface waves, Geophysics, vol. 52, pp. 60-71
(1987).
24. J.P. Berenger, A Perfectly Matched Layer for the Absorption of Electromagnetic
Waves, Journal of Computational Physics, vol. 114, pp. 185-200 (1994).
25. J. Davies and P. Cawley, The Application of Synthetically Focused Imaging
Techniques for High Resolution Guided Wave Pipe Inspection, Review of
Progress in Quantitative NDE, vol. 27, pp. 681-688 (2007).
215
26. F. Simonetti, Multiple scattering: The key to unravel the subwavelength world
from the far-field pattern of a scattered wave, Phys. Rev. E, vol. 73, 036619
(2006).
27. M. Fleming, M. J. S. Lowe, F. Simonetti and P. Cawley, Super Resolution
Imaging: Performance Studies, Review of Progress in Quantitative NDE, vol. 26,
pp. 820 (2006)
28. E.A. Skelton, An Overview of the Fictitious Domain Method, internal report,
Imperial College of Science Technology and Medicine, London, 2007.
29. F. Collino, P. Joly and F. Millot, Fictitious domain method for unsteady
problems: application to electromagnetic scattering, Journal of Computational
Physics, vol. 138, pp. 907–938 (1997).
30. E. Becache, P. Joly and C.Tsogka, Fictitious domains, mixed finite elements and
perfectly matched layers for 2D elastic wave propagation, I.N.R.I.A. report 3889,
2000.
31. M. Drozdz, L. Moreau, M. Castaings, M.J.S. Lowe and P Cawley, Efficient
numerical modelling of absorbing regions for boundaries of guided waves
problems, in Review of Progress in Quantitative NDE, vol. 25, pp. 126 (2005).
32. M. Castaings and C. Bacon, Finite Element modeling of torsional wave modes
along pipes with absorbing materials, J. Acoust. Soc. Am., vol. 119, pp. 3741-
3751 (2006).
33. B. A. Auld, Acoustic fields and waves in solids. (Krieger Publishing Company,
Florida, 1990), 2nd edition.
34. J.L. Rose, Ultrasonic Waves in Solid Media, Cambridge University Press ed.,
1999.
35. J.D. Achenbach, Wave propagstion in elastic solids. (North-Holland Publishing
Compagny, Amsterdam, 1984).
36. H. Lamb, On waves in an elastic plate, Conference of the Royal Society, 1917,
London.
37. M.J.S. Lowe, Matrix techniques for modelling ultrasonic waves in multilayered
media, IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency control,
vol. 42, pp. 525-542 (1995).
38. K.J. Bathe, Finite element procedures, Ed. Prentice Hall, 1996.
39. R. Courant, K. Friedrichs and H. Lewy, Über die partiellen
Differenzengleichungen der mathematischen Physik, Mathematische Annalen,
vol. 100, no. 1, pages 32–74, 1928
216
40. U. Basu and A. K. Chopra, Perfectly matched layers for time-harmonic
elastodynamics of unbounded domains: theory and finite-element
implementation, Computer Methods in Applied Mechanics and Engineering, vol
192, pp. 1337-1375 (2003).
41. F. Collino and C. Tsogka, Application of the PML Absorbing Layer model to the
Linear Elastodynamic Problem in Anisotropic Heterogeneous Media, Technical
Report 3471, INRIA, 1998.
42. A. Sommerfeld, Partial differencial equations in physics, 1949, New York:
Academic Press.
43. J. Lysmer and R.L. Kuhlemeyer, Finite Dynamic Model for Infinite Media,
ASCE Journal of the Engineering Mechanics 1969. 95(August): p. 859-877.
44. M. Cohen and P. C. Jennings, Silent Boundary Methods for Transient Analysis,
Computational Methods for Transient Analysis, Ed. T. Belytschko and T. R. J.
Hughes, Elsevier, 1983.
45. T. Dehoux, Finite Elements Modelling of Dissipative Materials, internal report,
Imperial College London, 2004.
46. G. Liu and J.S. Quek, A Non-Reflecting Boundary for Analysing Wave
Propagation Using the Finite Element Method, Finite Elements in Analysis and
Design, vol. 39, pp. 403-417 (2003).
47. D. Givoli, Non-reflective boundary conditions: A review, J. Comput. Phys., vol
94, pp. 1-29 (1991).
48. B. Engquist and A. Madja, Radiation boundary conditions for acoustic and
elastic calculation, Commun. Pure Appl. Math., vol. 32, pp. 313-357 (1979).
49. A. Bayliss and E. Turkel, Radiation boundary conditions for wave-like
equations, Commun. Pure Appl. Math., vol. 33, pp. 707-725 (1980).
50. D. Givoli and J.B. Keller, Non-reflecting boundary conditions for elastic waves,
Wave Motion, vol. 12, pp. 261-279 (1990).
51. J.B. Keller and D. Givoli, Exact non-reflecting boundary conditions, J. Comput.
Phys., vol. 82, pp. 172-192 (1989).
52. D. Givoli, High-order local non-reflecting boundary conditions: a review, Wave
Motion, vol. 39, pp. 319-326 (2004).
53. Q. Qi and T. L. Geers, Evaluation of the perfectly matched layer for
computational acoustics, J. Comput. Phys., vol. 139, pp. 166–183 (1994).
217
Chapter 4
References
54. E. Becache, A.-S. Bonnet Ben Dhia, and G. Legendre, Perfectly matched
layers for the convected Helmholtz equation, SIAM J. Numer. Anal. vol. 42, pp.
409–433 (2004).
55. I. Singer and E. Turkel, A perfectly matched layer for the Helmholtz equation in
a semi infinite strip, J. Comput. Phys., vol. 201, pp. 439–465 (2004).
56. F. D. Hastings, J. B. Schneider and S. L. Broschat, Application of the perfectly
matched layer (PML) absorbing boundary conditions to elastic wave
propagation, J. Acoust. Soc. Am., vol. 100, pp. 3061–3069 (1996).
57. F. Collino and C. Tsogka, Application of the perfectly matched absorbing layer
to the linear elastodynamic problem in anisotropic heterogeneous media,
Geophysics, vol. 66, pp. 294–307 (2001).
58. D. Komatitsch and J. Tromp, A perfectly matched layer absorbing boundary
condition for the second-order seismic wave equation, Geophysical Journal
International, vol. 154, pp. 146–153 (2003).
59. U. Basu and A. K. Chopra, Perfectly matched layers for transient
elastodynamics of unbounded domains, Int. J. Numer. Meth. Eng., vol. 59, pp.
1039–1074 (2004).
60. E.A. Skelton, A perfectly matched layer (PML) for elastic wave propagation,
internal report, Imperial College of Science Technology and Medicine, London,
2007.
61. D. Alleyne and P. Cawley, A two-dimensional Fourier transform method for the
measurement of propagating multimode signals, J Acoust Soc Am, vol. 89, pp.
1159-1168 (1991).
62. T.P. Pialucha, The reflection coefficient from interface layers in NDT of
adhesive joints, PhD thesis, Imperial College of Science Technology and
Medicine, London, 1992.
63. J. De Moerloose and M.A. Stuchly, Behavior of Berenger's ABC for evanescent
waves, Microwave and Guided Wave Letters, IEEE [see also IEEE Microwave
and Wireless Components Letters] , vol.5, no.10, pp.344-346, Oct 1995.
64. J. P. Berenger, Evanescent waves in PML's: origin of the numerical reflection in
wave-structure interaction problems, Antennas and Propagation, IEEE
Transactions on , vol.47, no.10, pp.1497-1503, Oct 1999.
218
Chapter 4
References
219
Chapter 4
References
77. L. Tomas, Une estimation d’erreur pour une formulation par domaine fictif avec
multiplicateurs de volume du probleme de Dirichlet, C. R. Acad. Sci. Paris Ser.
I325, pages 793–796, 1997.
78. A. Yu, O. Kuznetsov, O. D. Trufanov, Two-stage fictitious components method
for solving the wave Helmholtz equation, Sov. J. Numer. Anal. Math. Modelling,
5, 1988.
79. E. Heikkola, Y. A. Kuznetsov, P. Neittaanmaki, J. Toivanen, Fictitious domain
methods for the numerical solution of two dimensional scattering problems,
Journal of Computational Physics, vol. 145, pp 89–109, 1998.
80. A. Bespalov, Application of fictitious domain method to the Helmholtz equation
in unbounded domain, I.N.R.I.A. Report 1797, 1992.
220