# Comparisons of ray and finite element simulations of ultrasonic wave fields in smoothly inhomogeneous austenitic welds of thick-walled components

Nicolas Leymarie, Alexandre Imperiale, Thibaud Fortuna, Edouard Demaldent

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Nicolas Leymarie, Alexandre Imperiale, Thibaud Fortuna, Edouard Demaldent. Comparisons of ray and finite element simulations of ultrasonic wave fields in smoothly inhomogeneous austenitic welds of thick-walled components. NDT & E International, 2024, 147, 103177 (15 p.). 10.1016/j.ndteint.2024.103177. cea-04674499

HAL Id: cea-04674499

https://cea.hal.science/cea-04674499v1

Submitted on 21 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research Paper

# Comparisons of ray and finite element simulations of ultrasonic wave fields in smoothly inhomogeneous austenitic welds of thick-walled components

Nicolas Leymarie *, Alexandre Imperiale, Thibaud Fortuna, Edouard Demaldent

Université Paris-Saclay, CEA, List, F-91190, Palaiseau, France

## ARTICLE INFO

**Keywords:**

Ultrasonic modelling

Anisotropy

Weldment

Transient finite element

Dynamic ray tracing

## ABSTRACT

This work addresses the crucial challenges of simulating ultrasonic inspections in welded parts. We focus on the simulation of radiated ultrasonic fields in welded materials, incorporating spatially varying anisotropy and damping properties. Using macroscale descriptions to capture variations in weld properties as a function of grain orientation, two simulation methods are evaluated: a ray method and a hybrid finite element solution within the CIVA software. Simple simulations validate both models, but complexities emerge in real-world scenarios, notably caustic phenomena and beam splitting due to specific grain orientations. The proposed hybrid FE model can therefore provide reliable reference solutions, even for 3D configurations, allowing the biases arising from radius approximations to be quantified.

## 1. Introduction and motivation

Combining mechanical strength and corrosion resistance, austenitic steels are widely used, in particular for the assembly of components in the nuclear industry. The ability to detect fault-related weaknesses near the welding areas is a major issue for the safety of critical thick-walled components in the primary and auxiliary circuits of nuclear power plants, and remains a major industrial challenge in order to meet the new safety and durability requirements of these structures. Ultrasonic Testing (UT) is a common Non Destructive Testing (NDT) technique to inspect their structural integrity, but conventional UT may suffer from performance limitations into the weld due to the polycrystalline grain texture inducing attenuation, deviation or splitting of the ultrasonic beam. UT simulation can help understanding the influence and variability of these phenomena and therefore improving the inspection. To meet this need, it is necessary to offer versatile simulation tools dedicated to parametric studies, and to propose numerical models adapted to the applications of representative industrial cases.

Whereas simulation of ultrasonic wave propagation in elastic solids can be approached through various methods, available options become significantly constrained when it comes to proposing practical NDT solutions applied to complex industrial geometries. Two main families of models can therefore be distinguished. The first concerns the set of semi-analytical models based on a High Frequency (HF) approximation using a modal decomposition of the solution of the wave equation in elastodynamics, which will be referred to hereafter as the “ray-based approaches” [1–11]. The second corresponds to the large group

of Finite Element (FE) or finite difference numerical methods [12–16]. Even if these methods are now widespread and widely used, obtaining dedicated solutions adapted to the propagation of bulk waves in 3D welded structures requires the implementation of well-adapted strategies to reduce the memory footprint and the computational cost.

In addition to the particular geometries associated with the weld bevel, the modelling of the propagation of ultrasonic waves in the polycrystalline medium associated with the welded zone requires particular attention. The welding material properties are then generally modelled either by a heterogeneous description at the grain scale (microscale approach) or by a simplified description considering equivalent effective material properties at the wavelength scale ( macroscale approach). Inspection simulations based on a microscale approach are essentially performed with the help of FE solution [17–20]. Although well suited to the analysis of interactions between waves and grain microstructure on small Representative Elementary Volumes (REVs) [21,22], this microscale approach applied to the simulation of weld inspection has several limitations. First, a deterministic definition of the weld is impossible since we can never have perfect knowledge of the characteristics of each grain. The simulated microstructure can only be consistent with observed statistics, primarily concerning grain size distributions and crystallographic orientation distribution functions. Furthermore, simulating wave propagation demands an extensive computational capacity, coupled with a time frame that renders 3D parametric investigations challenging, if not impossible, without simplifying the problem (such as assuming geometric consistency along the weld deposit axis or imposing dimensional restrictions [20,23,24]).

<footnote>* Corresponding author.
E-mail address: nicolas.leymarie@cea.fr (N. Leymarie).</footnote>
Consequently, the macroscale approach is preferred in the majority of inspection simulations conducted on industrial welded structures of realistic dimensions [25–27]. The underlying assumption is that we consider grain scattering regimes at frequency ranges where the concept of coherent wavefronts remains relevant. This choice is justified when the wavelengths involved are slightly larger than the characteristic dimensions of the grains. In practice, this is the case for many UT applications in welded parts, where the wavelengths involved are several millimetres, while the grain size is less than one millimetre. In the context of austenitic welds, the dendritic solidification processes of austenite lead to columnar grain growth, which, upon solidification, forms an elongated grain microstructure that is highly textured along this growth axis [28]. The weld therefore has the characteristics of an anisotropic viscoelastic medium with a spatially variable anisotropic orientation according to the grain elongation axis also called fibre axis. In particular, a continuous representation of fibre axis variations leads to a smoothly inhomogeneous medium description.

Both FE and ray-based approaches are of well-suited to the macroscale approach. In particular, the ray-based solution implemented in the CIVA software [29] enables fast 3D simulations that are useful for variability and performance analysis, and has made it possible to improve the description and comprehension of various phenomena observed in weld inspections [30]. However, such paraxial ray tracing techniques suffer from certain limitations, particularly when it comes to the quantitative modelling of surface wave phenomena and of potentially observed intricate caustic occurrences [1,2,31]. As these artefacts can be frequently observed in welding structures, and although ray approaches significantly reduce computation times compared to FE methods, especially in 3D scenarios, it remains essential to propose simulation tools capable of validating or rejecting the ray approximations of the radiated field within the weld.

The initial motivation for this work was to provide a 3D FE reference solution to identify and quantify these potential artefacts. Based on the transient spectral FE code studied in [32–34], we propose an hybrid procedure using a ray calculation of the source field radiated by the probe at the surface of the component, while the rest of the calculation is performed by FEs on deformed structured quad meshes including viscoelastic anisotropic constitutive laws. This paper presents and compares results of ultrasonic field simulations performed with the two models available in CIVA. Our modelling considerations concerning the smoothly inhomogeneous weld description and the two different modelling tools are introduced in the first part of this paper. In the second part, a benchmarking analysis of the two simulation approaches is proposed. The analysis focuses on different cases of increasing complexity, to quantitatively validate the effects of anisotropy and attenuation and to identify situations where the limits of the ray models can be observed on unrealistic academic cases and on realistic weld descriptions.

## 2. Modelling approaches

The simulation of UT for industrial welding structures requires a comprehensive knowledge of the main geometric and mechanical characteristics associated with the weld, while the influence of uncertainties in the weld description should be taken into account by parametric studies. Fast, effective solutions are therefore needed for modelling welds and UT, and the most common approach consists in simulating the behaviour of the coherent wave used in ultrasonic inspection and imaging techniques. This work focuses on a macroscale (or wavelength scale) smoothly inhomogeneous description of the weld material, enabling fast ray-based simulations and optimized FE solutions. In this section, we first explain the practical application of this macroscale approach commonly used by the NDT community for simulation. We then introduce the two simulation models used in this study, the CIVA paraxial ray-tracing model and a new transient finite element solution, and describe their adaptations to the proposed smoothly inhomogeneous weld. Hereafter, geometries are limited to V-shaped welds for the sake of simplicity.

<caption>**Table 1**</caption>

Stiffness coefficients (Voigt form) used for effective material properties of austenitic weld including damping properties inspired from [44] with grain orientation along the first axis.

<table><thead><tr><th>Effective stiffnesses (GPa)</th><th>C<sub>11</sub></th><th>C<sub>22</sub> = C<sub>33</sub></th><th>C<sub>12</sub> = C<sub>13</sub></th><th>C<sub>23</sub></th><th>C<sub>44</sub></th><th>C<sub>55</sub> = C<sub>66</sub></th></tr></thead><tbody><tr><th>Real part</th><td>220</td><td>245</td><td>135</td><td>110</td><td>80</td><td>110</td></tr><tr><th>Imaginary part</th><td>1.6</td><td>12</td><td>2</td><td>2.5</td><td>2.7</td><td>4.3</td></tr></tbody></table>

### 2.1. Macroscale weld material description

When considering wave phenomena in the Rayleigh scattering regime for polycrystalline materials, it is common to describe the properties of the weld as those of a homogenized material, with variable properties that are spatially correlated with the local orientation of grain elongation. Although each austenite grain is characterized by a face-centred cubic (FCC) structure, the preferred texture associated with the grain elongation axis imparts orthotropic or transversely isotropic macroscopic properties to the material [18]. Therefore, the equivalent material can be described using such anisotropic properties, with the axis of symmetry normal to the (quasi-)isotropic plane aligned with the grain elongation axis [35,36]. By means of this approach, weld inspection simulations can be performed on representative geometries of a weld component. In assigning to the effective medium the properties of an anisotropic and attenuating material, note that this strategy only accounts for the behaviour of the coherent wave. In this case, the structural noise is not directly modelled but can be added a posteriori using a set of fictitious scatterers into the weld [37].

#### Effective material properties in the columnar grain frame

Let us consider a REV of the weld in which the columnar grains have similar morphologies and orientations. Then, the analysis of the coherent wavefronts in various directions shows wave velocity and attenuation variations representative of an homogeneous quasi-hexagonal viscoelastic material where the elongation axis of the grains corresponds to the normal direction of the quasi-isotropic plane. The viscoelastic character or damping properties helps to account for the energy loss of the coherent wavefront due to scattering at each grain interface. It is often transcribed in the form of a complex-valued stiffness tensor at a fixed frequency, and different rheological viscosity models are available with the help of parametric relaxation functions [34,38,39]. As an example, the Maxwell’s model provides a quasi-constant attenuation around the frequency of interest whereas the Zener one is almost linear.

Such effective anisotropic and viscosity properties can be obtained in different ways. Either by means of experimental measurements on samples taken from representative blocks of the welded material [35, 36] or directly deduced from a model of the microstructure (statistics on grain geometries and crystallographic textures) by using effective wavenumbers estimations performed by scattering models [40–43]. In what follows, the weld material is expressed in a reference coordinate system denoted $(e_1^R, e_2^R, e_3^R)$ with columnar grain orientation along the first axis. Its properties are given by the complex stiffness tensor in Table 1 [44] estimated at the central frequency of 2 MHz where the plane of isotropic symmetry is perpendicular to the elongation axis. Note that in the following we use $c_{ijkl}$ to define the components of the fourth order stiffness tensor and $C_{ij}$ in upper case to express the same values in the contracted Voigt form.

#### Grain orientation variation into the weld

The next step involves the ability to describe the variations in columnar grains throughout the entire weld. These descriptions are generally defined in 2D, the third dimension being deduced using an extrusion according to the geometry of the component (axial or revolution extrusion). The medium is considered homogeneous for most industrial imaging purposes, but more advanced approaches involve inhomogeneous descriptions as illustrated in Fig. 1. A first one roughly
<figure><img src="image_3.png" /><figcaption>Fig. 1. Main types of representations of grain orientation variations in the weld from the most simplistic homogeneous and piecewise homogeneous model to orientation mapping descriptions (closed-form expression or MINA grid data).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Example of a Gaussian blur filter applied to raw MINA data to regularize the description inside the welded part (standard deviation 2 mm).</figcaption></figure>

describes the weld as a set of homogeneous domains with different crystallographic orientations (macrographic image analysis, see e.g. [45]). Other approaches offer smooth, idealized variations based on observations and defined using closed-form expressions such as Ogilvy's laws [46]. Intermediate approaches provide grid representations in which large pixels fix the dominant orientation locally (macrographic image analysis or MINA-like model [47]). Excessively noisy (especially in the case of macrographic analysis) or too coarse orientation data are not uncommon when employing such grid descriptions, as illustrated in Fig. 1. This motivates the application of a regularization procedure, such as a spatial Gaussian blur whose aperture size is adjustable according to the observed orientation variations, illustrated in Fig. 2. In what follows, smoothly inhomogeneous variation in the orientation of columnar grains is assumed, such as the Ogilvy or regularized grid descriptions.

Defining a coordinate system ($e_1, e_2, e_3$) at each position $M$ in the weld material, where $e_1(M)$ corresponds to the local orientation of the grains, the local stiffness values are given by

$$c_{ijkl}(M) = J_{im}(M)J_{jn}(M)J_{kp}(M)J_{lq}(M)c_{mnpq}^R \quad (1)$$

where $J_{im}(M) = e_i(M) \cdot e_m^R$ and $c_{mnpq}^R$ corresponds to the effective properties of Table 1.

## 2.2. The CIVA ray-based model for smoothly inhomogeneous welds: RAY BEAM

The aim is not to describe in detail the ray model calculation procedures used in CIVA [48], but rather to outline them, as they

<figure><img src="image_5.png" /><figcaption>Fig. 3. Ray calculation principles and application in the case of welding specimen.</figcaption></figure>

are largely inspired by techniques found in the geophysical literature and summarized exhaustively in these reference works [1,2]. The ray-based model implemented in CIVA was initially based on the so-called “pencil method” [4,6,49]. This technique actually uses the same principle as the dynamic ray tracing model in a wave front orthonormal coordinates system early developed in the geophysical literature. In order to determine the ultrasonic field observed at an observation point $M$ radiated by the active surface of the transducer, we calculate a set of ray paths linking $M$ and the active surface, discretizing it as a set of source points (see Fig. 3). The calculation of these ray paths can be summarized as a succession of propagation in an homogeneous or smoothly inhomogeneous medium and refraction (or reflection) at the interface between two distinct media. In an homogeneous fluid or elastic isotropic media, the calculation principle for ray tracing is well known, and is expressed in terms of wave propagation velocities with a straight line trajectory. For the interface crossing (discontinuity in material properties), Snell's laws are used to determine the new properties of the refracted or reflected wave according to the new mode considered (with or without conversion). Solving the system is much more complex in an inhomogeneous anisotropic material properties. Concerning refraction, the calculation is more onerous in an anisotropic medium but there is no additional complexity due to the inhomogeneity of the medium under consideration. In the following, we focus on the particularities of propagation management in smoothly inhomogeneous parts.

A ray is usually defined by both its position vector and its slowness vector, where $x_i$ and $p_i$ denote the position and slowness components respectively. Using the general elastodynamic equation (see Cerveny [2] §3.6), the ray trajectory for a fixed wave mode $m$ (with polarization components denoted $g_i$) can be estimated by solving two coupled ordinary differential equations (ODEs), called axial ray system, with respect to the travel time $T$ such as

$$\begin{cases} \frac{dx_i}{dT} = a_{ijkl} p_l g_j g_k = \mathcal{V}_i, \\ \frac{dp_i}{dT} = -\frac{1}{2} a_{ijkl,i} p_k p_n g_j g_l = \eta_i. \end{cases} \quad (2)$$

where $a_{ijkl} = \rho^{-1}c_{ijkl}$ with $\rho$ the constant density of the material considered. $\mathcal{V}_i$ stands for the components of the energy velocity vector. In the second ODE, the perturbation of the slowness components along the ray path is related to the gradient of the elasticity tensor and can be interpreted as a continuous refraction along the trajectory of the wavefront. This procedure is equivalent to the one proposed in [50,51]. Note that in the case of an homogeneous medium, the components of $\eta_i$ are null and the ray tracing solution corresponds to a straight line trajectory in the energy direction estimated for the given initial value of the wavefront direction ($p_i$ components remaining constant). Hence, while trajectory calculations for piecewise homogeneous descriptions are estimated using straight ray paths, smoothly inhomogeneous descriptions lead to trajectories of curved ray paths (see Fig. 3).
In order to estimate the amplitude variation due to geometrical spreading of the wavefront, we have to use a second system, called the paraxial ray system, describing the evolution of the ray tube (or pencil) with respect to the axial ray defining by ray tube infinitesimal parameters ($\delta x_l$ and $\delta p_l$, i.e. 6 parameters in 3D with $i = 1, 2, 3$). In the case of the pencil method used in CIVA, the paraxial ray quantities are managed in a wavefront orthonormal coordinate system (see Cerveny [2] §4.2.2 and §4.14.2) reducing the infinitesimal ray tube properties to projected pencil quantities in the wavefront plane denoted here ($\delta x_l, \delta p_l$) with $I = 1, 2$ in 3D, i.e. 4 parameters in all. The evolution of paraxial quantities from a source point $S$ to an observation point $M$ can be expressed in the form of a $2 \times 2$ block matrix, called the propagator matrix $\Pi$, such as

$$ \left\{ \begin{aligned} \delta x_I(M) &= \Pi_{IJ}^{xx}(M/S) \delta x_J(S) + \Pi_{IJ}^{xp}(M/S) \delta p_J(S), \\ \delta p_I(M) &= \Pi_{IJ}^{px}(M/S) \delta x_J(S) + \Pi_{IJ}^{pp}(M/S) \delta p_J(S). \end{aligned} \right. \qquad (3) $$

Contrary to homogeneous media where the solution is estimated analytically with respect to the rectilinear trajectory of the axial ray (see e.g. [4]), the calculation of the propagator matrices along a ray path in inhomogeneous material (see Cerveny [2], §4.3.1) is done by solving another ODEs along the curved ray path using a adaptive step size procedure [52] according to local material properties (stiffness values and gradients). In practice, the schemes used to solve these ODEs depend on the application. To determine the set of rays linking our extended source (transducer) to the zone of interest (defect or ultrasonic field zone), we use the fast Cash–Karp method combining Runge–Kutta schemes of order 4–5. On the other hand, to determine the evolution of the paraxial quantities and easily increment the KMAH index (following procedure given in [53]), we retain a low-order Euler-type scheme, essentially to ensure the stability of the calculations. Assuming the conservation of energy flux across a ray tube section, the geometrical spreading can also be estimated, even in inhomogeneous anisotropic media. Then, the geometrical spreading $\mathcal{L}(M/S)$ of the ray tube is evaluated thanks to the matrix $\Pi^{\text{xp}}(M/S)$ with $\mathcal{L}(M/S) = |\det \Pi^{\text{xp}}(M/S)|^{\frac{1}{2}}$.

In the case of a viscoelastic medium, considered later, this approach remains valid and leads us to manipulate complex quantities implying a specific ray tracing procedure as presented in the papers of Vavryčuk [54,55]. This procedure can be approximated by a simpler approach neglecting dispersion effects for materials presenting sufficiently high quality factors. For a given real phase direction, only the real part of the energy velocity vector is then considered to define the associated real ray direction. However, the slowness deduced from the Christoffel equation using a complex stiffness tensor remains a complex-valued solution. The associated attenuation is then related to the imaginary part of a complex travel time. In the case of an inhomogeneous medium, the complex time of flight $\tau(M/S)$ is estimated by the integration along the ray path $\ell(M/S)$ of the real travel time quantity such as

$$ \tau(M/S) = \int_{\ell(M/S)} 1 + i \frac{\mathfrak{F}(p_l)\mathcal{V}_l}{\Re(p_l)\mathcal{V}_l} dT. \qquad (4) $$

Finally, at each interface between two distinct materials, in addition to estimating the evolution of the axial and paraxial quantities applying the Fermat's principle, it is also necessary to calculate the transmission or reflection coefficients of the mode conversion. The product of all these coefficients calculated along the path $\ell(M/S)$ is cumulated in a global factor denoted $A_{\text{inter}}(M/S)$.

In practice for NDT applications, the probe is not considered as a point source but as an extended surface. The emitting surface of the ultrasonic probe is discretized using a backward ray strategy starting from the observation point (as illustrated Fig. 3), imposing as initial conditions a phase direction $p_{l0}(M)$ and an aperture cone $\delta p_{l0}(M)$ for each paraxial ray. Once the set of trajectories $\ell(M/S_n)$ linking the probe surface to the observation point $M$ has been determined, the

<figure><img src="image_6.png" /><figcaption>Fig. 4. Typical 2D description used to model ultrasonic testing of a welded joint component.</figcaption></figure>

displacement field in the positive frequency domain, with $\omega$ standing for the circular frequency, is given by

$$ u_l(M, \omega) = K V_0(\omega) \sum_{S_n \in I_p} \frac{A_{\text{inter}}(M/S_n)}{\mathcal{L}(M/S_n)} g_l(M/S_n) \\ \times e^{i[\omega\tau(M/S_n) - \frac{\pi}{2}\sigma(M/S_n)]} \Delta S_n, \qquad (5) $$

$\Delta S_n$ being the estimated effective surface area of the probe determined using the paraxial ray parameters and $K$ a common constant determined by the material properties of the source and the observation point $M$. In order to introduce potential phase shift due to caustics, the KMAH index $\sigma$ counts, as mentioned in [1,2], the caustics along the ray path $\ell(M/S)$ incrementing by different integer values with respect to local properties of the matrix $\Pi^{\text{xp}}$ checked for each step of the iterative procedure concerning propagation in inhomogeneous media (see [2] §4.14.13). The limitations of these high-frequency approaches are numerous: caustics and localized over-intensities, geometric shadowing and incorrect representation of frequency-dependent diffraction phenomena, surface phenomena at the interfaces between two different media. If most of the time these approximations have a negligible effect on the ultrasonic response of a localized defect, it is nevertheless difficult to affirm that the approximations of the ray model are without influence, especially when dealing with complex geometrical and material properties such as those encountered for welding parts.

## 2.3. Proposed ray-FE hybrid solution: FE BEAM

An alternative model has been developed to overcome the limitations of the previous RAY BEAM model. It combines the efficiency of a ray model to simulate the field radiated by a sensor in its surrounding environment with the robustness of a FE model to describe the propagation phenomena within the welded part.

### Model problem

The addressed problem is now defined as illustrated Fig. 4 with: a coupling medium denoted $\Omega_f$, here considered homogeneous and fluid; the two domains on either side of the weld denoted $\Omega_l$ and $\Omega_r$, made of the same base metal considered as a homogeneous isotropic steel; the welded part, denoted $\Omega_w$, with anisotropic and potentially smoothly inhomogeneous material properties. The complete domain of propagation in the solid is denoted by $\Omega_s = \Omega_l \cup \Omega_w \cup \Omega_r$. The continuity of the normal velocity and traction forces is assumed at the fluid/solid interface denoted $\Gamma_{fs}$ while a free surface condition is assumed at the backwall $\Gamma_b$. The probe surface included in $\Omega_f$ is denoted $\Gamma_p = \bigcup_{n=1}^N \Gamma_{p,n}$, $N$ being the number of elements of a phased-array probe. The surface source is considered to be spatially uniform for a given element $n$, with a common time dependence $V_p(t)$ corresponding to the particle velocity normal to the element probe surface $\Gamma_p$.
Thus, formally, the model problem can be summarized using scalar potentials $\phi$ in acoustics and displacement vector $u = (u_1, ..., u_d)^T$ (with $d=2$ in 2D, or $d=3$ in 3D) in solids such as

$$ 
\left\{ 
\begin{array}{l} 
\partial_{tt}^2 \phi - c_f \partial_{tt}^2 \phi = 0 \text{ in the fluid part } \Omega_f, c_f \\ 
\text{being the pressure wave velocity,} \\ 
\rho \partial_{tt}^2 u_i - \partial_j \sigma_{ij} = 0 \text{ with } \sigma_{ij} = c_{ijkl} \partial_l u_k \text{ in the solid part } \Omega_s, \\ 
c_f \partial_i \phi n_i = V_p \text{ on } \Gamma_p, \\ 
\partial_l \phi n_i = \partial_l u_i n_i \text{ and } \sigma_{ij} n_j = \partial_i \phi n_i \text{ on } \Gamma_{fs}, \\ 
\sigma_{ij} n_j = 0 \text{ on } \Gamma_b, \text{ the backwall surface of the part.} 
\end{array} 
\right. 
\quad (6) 
$$

### FE problem

For a computational domain $\Omega$, we denote by $\mathcal{T}_h(\Omega)$ a tessellation of $\Omega$ made of non-overlapping quadrangles (in 2D) or hexahedra (in 3D). The associated finite element approximation space reads

$$ V_h(\Omega) = \{ v_h \in C^0(\bar{\Omega}) \mid \forall K \in \mathcal{T}_h(\Omega), v_h|_K \in Q^p(K) \}, \quad (7) $$

where $C^0(\bar{\Omega})$ is the space of continuous function over the complete domain (boundary included), and $Q^p(K)$ is the space of polynomial of order $p$ in a mesh element $K$. FE formulations are based on the weak form of (6) – also referred to as the principle of virtual works – built upon discrete approximation spaces. Namely, for any time $t \ge 0$ we seek the discrete scalar potential $\phi_h(t) \in V_h(\Omega_f)$ and the discrete displacement vector $u_h \in [V_h(\Omega_s)]^d$ such that for any test functions $\psi_h \in V_h(\Omega_f)$ and $v_h \in [V_h(\Omega_s)]^d$

$$ 
\left\{ 
\begin{array}{l} 
\displaystyle \int_{\Omega_f} \left( \partial_{tt}^2 \phi_h \psi_h + c_f \partial_i \phi_h \partial_i \psi_h \right) d\Omega = \int_{\Gamma_p} V_p \psi_h d\Gamma \\ 
\displaystyle \quad + \int_{\Gamma_{fs}} c_f (\partial_i \phi_h n_i) \psi_h d\Gamma, \\ 
\displaystyle \int_{\Omega_s} \left( \rho \partial_{tt}^2 u_{h,i} v_{h,i} + c_{ijkl} \partial_l u_{h,k} \partial_j v_{h,j} \right) d\Omega = \int_{\Gamma_{fs}} (\sigma_{h,ij} n_j) v_{h,i} d\Gamma, \\ 
\displaystyle \quad \partial_i \phi_h n_i = \partial_i u_{h,i} n_i \text{ and } \sigma_{h,ij} n_j = \partial_i \phi n_i \text{ on } \Gamma_{fs}. 
\end{array} 
\right. 
\quad (8) 
$$

Following the previous work of [56], the fluid-solid transmission conditions are taken into account using the mortar element method [57–59]. This amounts to introducing two auxiliary scalar unknowns $\lambda_{h,f}$ and $\lambda_{h,s}$ both in $V_h(\Gamma_{fs})$ and such that

$$ \partial_t \lambda_{h,f} = \partial_t \phi_{h,n_i}, \quad \partial_t \lambda_{h,s} = \sigma_{h,ij} n_j n_i. \quad (9) $$

The complete weak form of the discrete problem reads, for any tests functions $\mu_h$ and $v_h$ in $V_h(\Gamma_{fs})$

$$ 
\left\{ 
\begin{array}{l} 
\displaystyle \int_{\Omega_f} \left( \partial_{tt}^2 \phi_h \psi_h + c_f \partial_i \phi_h \partial_i \psi_h \right) d\Omega - \int_{\Gamma_{fs}} c_f \partial_t \lambda_{h,f} \psi_h d\Gamma = \int_{\Gamma_p} V_p \psi_h d\Gamma, \\ 
\displaystyle \int_{\Omega_s} \left( \rho \partial_{tt}^2 u_{h,i} v_{h,i} + c_{ijkl} \partial_l u_{h,k} \partial_j v_{h,j} \right) d\Omega - \int_{\Gamma_{fs}} \partial_i \lambda_{h,s} v_{h,i} n_i d\Gamma = 0, \\ 
\displaystyle \int_{\Gamma_{fs}} (\lambda_{h,f} - u_{h,i} n_i) \mu_h d\Gamma = 0, \\ 
\displaystyle \int_{\Gamma_{fs}} (\lambda_{h,s} - \phi_h) v_h d\Gamma = 0. 
\end{array} 
\right. 
\quad (10) 
$$

In practice, the FE approximation spaces $V_h$ are built using Spectral Element Method (SEM). The SEM is a specific high-order mass-lumped FE method that provides after time discretization a fully explicit time- transient solver as the mass matrix is diagonal by construction. This technique has been initially used in computational fluid [60] dynamics then extended to 3D elastodynamics for seismic wave applications [14]. Thanks to its compact memory footprint, the SEM can manage 3D calculation configurations in domains covering dozens of wavelengths in all three dimensions.

### Hybridization with a ray model

The distance between the emitting boundary $\Gamma_p$ and the solid parts is usually large, which implies from a FE computation perspective a great number of degrees of freedom to be added for the sole purpose of propagating the acoustic field within an homogeneous fluid medium. To overcome this computational bottleneck we resort to an hybrid ray/FE solution based on the principle of separation of the incident field generated by the source and the total field, which encompasses every contributions of the solid part. The incident field is obtained using a ray model described in the previous section, and is denoted by $\phi_{\omega}^{\text{inc}}$. Note that the subscript $\omega$ is used to refer to the fact that ray-based approach are in essence built upon a high-frequency asymptotic approach. The incident field satisfies the following problem

$$ \partial_{tt}^2 \phi_{\omega}^{\text{inc}} - c_f \partial_{tt}^2 \phi_{\omega}^{\text{inc}} = 0 \text{ in } \mathbb{R}^d, \quad c_f \partial_i \phi_{\omega}^{\text{inc}} n_i = V_p \text{ on } \Gamma_p. \quad (11) $$

For the total field computation in the solid parts we use the previous FE discrete weak form. To inject the incident field in this formulation, we follow the strategy proposed in [33], Appendix B. It consists in solving in the fluid domain the difference between the incident and total field, referred to as the diffracted field, while the boundary conditions at the fluid-solid coupling are modify using the incident field. Hence, the fluid-solid coupling interface also becomes the coupling interface, as depicted in Fig. 5. More precisely, we consider the following decomposition

$$ \phi_h^{\text{tot}} = \phi_h^{\text{diff}} + \phi_{\omega}^{\text{inc}}, \quad \partial_t \lambda_{h,f} = \partial_t \phi_h^{\text{diff}} n_i + \partial_i \phi_{\omega}^{\text{inc}} n_i, \quad (12) $$

so that the final hybrid discrete formulation corresponds to finding the discrete diffracted scalar potential $\phi_h^{\text{diff}} \in V_h(\Omega_f)$, the discrete total displacement vector $u_h^{\text{tot}} \in [V_h(\Omega_s)]^d$, and the auxiliary unknowns $\lambda_{h,f}$ and $\lambda_{h,s}$ in $V_h(\Gamma_{fs})$ such that

$$ 
\left\{ 
\begin{array}{l} 
\displaystyle \int_{\Omega_f} \left( \partial_{tt}^2 \phi_h^{\text{diff}} \psi_h + c_f \partial_i \phi_h^{\text{diff}} \partial_i \psi_h \right) d\Omega - \int_{\Gamma_{fs}} c_f \partial_i \lambda_{h,f} \psi_h d\Gamma \\ 
\displaystyle \quad = - \int_{\Gamma_{fs}} c_f \partial_i \phi_{\omega}^{\text{inc}} n_i \psi_h d\Gamma, \\ 
\displaystyle \int_{\Omega_s} \left( \rho \partial_{tt}^2 u_{h,i}^{\text{tot}} v_{h,i} + c_{ijkl} \partial_l u_{h,k}^{\text{tot}} \partial_j v_{h,j} \right) d\Omega - \int_{\Gamma_{fs}} \partial_i \lambda_{h,s} v_{h,i} n_i d\Gamma = 0, \\ 
\displaystyle \int_{\Gamma_{fs}} (\lambda_{h,f} - u_{h,i}^{\text{tot}} n_i) \mu_h d\Gamma = 0, \\ 
\displaystyle \int_{\Gamma_{fs}} (\lambda_{h,s} - \phi_h^{\text{diff}}) v_h d\Gamma = \int_{\Gamma_{fs}} \phi_{\omega}^{\text{inc}} v_h d\Gamma. 
\end{array} 
\right. 
\quad (13) 
$$

Note that in the hybrid formulation (13), both the incident acoustic potential and its normal space derivative appear as right-hand sides, which is typical of transmission conditions. Finally, the fluid domain is bounded using Perfectly Matched Layers (PMLs) [61], see Fig. 5, entailing a significant gain of performances, especially when the distance between $\Gamma_p$ and $\Omega_s$ is large.

### Implementation strategy within the CIVA platform

Emphasis has been placed on two critical aspects essential for NDT studies within the CIVA platform: enabling 2D and 3D calculations on standard workstations and efficiently handling parametric variations in configurations of interest. In order to meet these objectives with the SEM model, the geometry associated with the UT configuration is decomposed into parametric subdomains composed of macro-elements (MEs). Each ME is meshed and defined as a deformation of a reference cube using a predefined hexahedral mesh and is considered as a physical domain or a perfectly matched absorbing layers for which specific characteristics are defined according to the properties of the material (fluid, solid). Two neighbour MEs communicate with each other using the mortar element method defined at their interface. Building the final mesh through a structured refinement of each ME allows unassembled
<figure><img src="image_7.png" /><figcaption>Fig. 5. Typical domain decomposition of a UT scene defined in CIVA for weld sample and identification of MEs with their specific features according to material properties.</figcaption></figure>

FE operations that improve the performance of standard FE procedures in terms of memory footprint and computational load with an optimal data structure for multi-threaded computations [32–34].

Three main stages allow an automated calculation procedure without requiring expertise on the meshing procedure and other constraints of the FE model. Initially, the process entails constructing the mesh based on the UT scene properties, which include geometry, materials, and central frequency. In the specific case of welded parts, a dedicated graphical user interface (GUI) is proposed in CIVA to parametrically describe a selection of weld geometries decomposed into MEs. Additionally, determining the optimal time step for FE computation, ensuring stability conditions, is crucial at this stage. Subsequently, the calculation of the incident field at the coupling points is conducted through the ray model. Finally, the FE calculation is performed across the entire sample, with the total field (displacement or velocity) recorded at specified observation points.

The calculation of the radiated field by the ray method is simple and we can benefit from all the possibilities offered by CIVA on the management of the various forms of probes encountered with phased array transducers. Apart from the near-field effects close to the emitting source, the solution of the asymptotic ray model provides an accurate approximation of the radiated field without including issues related to ray artifacts previously described. In addition, the smoothly inhomogeneous effective description of the weld material suits the high-order structured (easy) automatic meshing of the weld bevel. The hybrid solution thus implemented avoids certain approximations and singularities of the ray model by providing a solution serving as a reference while allowing the execution of 3D calculations on component sizes of several tens of wavelengths with optimal performance.

### 3. Benchmarking of models

For a given defect, the simulated echo amplitude error is proportional to the square of the incident field intensity error (factor 2 in dB) in a pulse echo acquisition mode. The accuracy of the simulated field therefore has a direct impact on the simulated prediction. This section

focuses on a comparative analysis of the ultrasonic field radiated by a transducer between the ray-based model (RAY BEAM) and the proposed hybrid FE solution (FE BEAM). The complexity of the reference cases studied is gradually increased: first on simple and homogeneous blocks to validate the incorporation of anisotropy and attenuation phenomena, then on heterogeneous blocks before studying 2D and 3D welded assemblies configurations.

#### 3.1. Viscoelastic and anisotropic effects on a simple geometry

##### Homogeneous parts

The impact of anisotropy and attenuation on the coherent wave is first studied on a homogeneous parallelepiped block under normal incidence (L0). This approach eliminates any supercritical phenomena and reduces the influence of shear waves, which facilitates the modal analysis of the main longitudinal wave while preserving energy deviation phenomena. Different orientations of the same hexagonal medium were obtained by rotating the axis of symmetry normal to the plane of transverse isotropy (here the axis 1 for the material given in Table 1). Applied to a uniform microstructure of elongated grains, this amounts to varying the orientation of the principal axis of grain elongation with respect to the principal wavefront, as illustrated in Fig. 6. The shape of the simulated probe is circular with a diameter of 12.7 mm and the centre frequency is 2 MHz with a 6 dB bandwidth of 75%. Shear modes (qS) and reflections/conversions at the bottom of the part are not taken into account within RAY BEAM so that only the contribution of the direct quasi-longitudinal (qL) wave is calculated, while simulated propagation time is limited to 23 µs within FE BEAM in order to avoid the reflection on the backwall of the part, which is 40 mm thick.

As the wavefront remains mainly parallel to the surface, the beam energy shifts to the right or left depending on the orientation of the axis considered. This phenomenon, as well as beam attenuation due to viscosity (imaginary part of the stiffness tensor), are clearly observed on the modulus of the displacement field simulated with RAY BEAM (Fig. 7) and FE BEAM (Fig. 8). The results are in good accordance, with increasing attenuation as a function of the orientation angle
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Illustration of the different scenarios proposed for an effective homogeneous anisotropic medium. The grain structure is shown to illustrate the correspondence between the direction of grain elongation and the **1** axis of hexagonal symmetry (red arrows).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Maximum displacement modulus in time for the different cases simulated with RAY BEAM: top without attenuation (elastic model), below with attenuation (viscoelastic model).</figcaption></figure>

of the grains with respect to the wavefront. Minor differences are observed, especially when the material is considered inviscid, and when grain orientations are such that the wave propagates away from the symmetry axis, as in the cases of the 30° and 60° orientations. What may appear to be artefacts are in fact traces of slower shear waves generated with FE that have not had time to propagate completely in the region of interest. This is confirmed by the ray paths illustrated in Fig. 8 (green line for qL and red lines for qS), which show the positions of the wavefronts at the final simulation time considered. This is no longer observed when material viscosity is considered, as the qS modes are much more attenuated than the qL one.

A quantitative assessment of beam deflection and energy attenuation with respect to variations in orientation in increments of 5 degrees is given Fig. 9. The two models are compared with theoretical values obtained with a plane wave (PW) model using the probe's focal axis as the main direction of the wavefront, always normal to the sample surface (LO). These results are obtained by solving the Christoffel equation in the viscoelastic case [38], which allows us to determine the complex wave numbers associated with the complex elasticity tensor defined in Table 1. We consider this to be an excellent agreement that validates the models. The slight deviations observed are mainly located in the presence of strong attenuations with grain orientation close to 90°, essentially on the RAY BEAM field whose attenuation is

underestimated. These differences are attributed to a shift to lower frequencies of the ray model whose propagator has a linear frequency dependence of its viscosity. FE BEAM is less sensitive to this frequency shift as the Maxwell rheological viscosity law selected provides a quasi-constant attenuation around the frequency of interest [34,39]. Note that the estimation of the attenuation was not done by a frequency spectral analysis procedure of the displacement field according to z, as it would require a tedious process to distinguish the effects of losses related to material viscosity from those related to the divergence of the beam. Instead, a ratio between the measurement of the maximum amplitude of the displacement modulus with and without taking into account the material viscosity was applied to remove the loss effects related to the beam divergence.

### Smoothly inhomogeneous parts

Similar blocks whose grain orientations varies with depth, similar to a continuously variable stratification according to the z direction, are now considered. This assumption has the merit of highlighting the effects of local variations while allowing easy interpretation of variations in energy direction by applying Snell's laws in the case of a stratified plane medium. Two descriptions are proposed Fig. 10. The first exhibits a linear variation in grain orientation from 90° to 0°. The second maintains a global grain orientation of 90° with a local Gaussian variation for a 60° grain orientation centred at a depth of 12.7 mm.
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Maximum displacement modulus in time for the different cases simulated with FE BEAM: top without attenuation (elastic model), below with attenuation (viscoelastic model). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Measurements at 30 mm depth on simulated results (RAY BEAM and FE BEAM) of beam energy deflection and apparent attenuation, comparison with the theoretical solution using a plane wave assumption at 2 MHz.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Global linear (top, red) and local Gaussian (bottom, green) grain orientations within the smoothly inhomogeneous blocks, L0 (left) and L60 (right ray path).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 11. Simulation results on smoothly inhomogeneous blocks: global linear (1st row) vs. local Gaussian (2nd row) variations, normal incidence L0 (1st column) vs. oblique incidence L60 (2nd column).</figcaption></figure>

Simulations were performed for a normal incidence (L0) and an oblique incidence where the beam is refracted at 60° of the normal (L60). The energy direction is deviated by about 60° in oblique incidence (12.5° in a coupling medium made of water) for a grain orientation of 90° with respect to the z-normal at z = 0 mm. Getting closer to the grain orientation of 60°, the energy directions become tangent to the surface corresponding to a set of “critical layers”, highlighting phenomena that penalize the ray approach. The norm of the projected slowness in the z-normal plane is conserved according to Snell’s laws. In the oblique case this norm is of the order of 0.175 µs mm⁻¹ and implies an over-critical incidence for a range of grain orientations between 16° and 62° with respect to the z axis. On Fig. 10, for an oblique incidence, one can anticipate with the help of the ray tracings that the linear description leads to a diffuse grazing wave phenomenon whereas the local Gaussian variation reveals a turning point with potential shadow zones (phenomena well described e.g. in [62] chapter 9 for acoustic applications).

Simulated radiated field are grouped in Fig. 11. While the RAY BEAM and FE BEAM results are similar regardless of the layer variation considered at normal incidence (L0), this is not necessarily the case at oblique incidence (L60). Indeed, for the linear grain orientation variation, the beam energy distribution patterns remain similar and the amplitudes are very close. However, for the Gaussian grain orientation variation, the amplitudes are less similar, particularly to the right of the region of interest, where a splitting of the beam is observed. This is confirmed by the quantitative comparisons of the amplitude profiles given Fig. 12. While the magnitudes of the field are quite consistent between the two models under a linear variation description, this consistency breaks down under the local Gaussian variation, leading to an observed over-intensity at 9 mm depth and a too strong decay at 13 mm of the ray model with respect to FE results. These two zones correspond to distinct areas in the ray model on each side of a fold caustic curve: a wavefront folding zone above and a penumbra on the other side. An analysis of ray tracing indeed reveals the high sensitivity to incidence and thus the orientation of the wavefront for this particular case. As shown in Fig. 13, a variation of a few tenths of a degree drastically changes the behaviour of the rays, with rays transmitted through the Gaussian layer, a fold caustic, and a folding of the wavefront. Due to the bounded nature of the beam, the field is nothing but a combination of these phenomena, which explains the observed beam splitting. If such deviations can be observed in the presence of variations in layered anisotropy orientation, there is a high probability that these RAY BEAM biases will be reproduced in weld descriptions and for certain wavefront orientations.

## 3.2. Application to representative welded part configurations

### 2D configuration

Configurations closer to industrial cases with a V-shaped weld bevel are now considered. The ability of the CIVA ray model to provide consistent results with the 2D FE calculation performed with the ATHENA code has already been demonstrated in previous studies [30] for a symmetrical Ogilvy’s description of the weld material. Hereafter, the same symmetrical description of the grain orientation is considered and an asymmetrical description is introduced Fig. 14. The latter is derived from the grain orientations shown in [63], which is characteristic of a welding process performed in a horizontal position. A linear phased array probe was chosen considering the same delay law for both grain descriptions. This law is estimated in an isotropic steel part in order to impose a refraction of the beam with a 45° L-wavefront (L45) and a focusing at the depth of the backwall (here 40 mm). Comparisons of the normalized (in space) maximum (in time) simulated displacement modulus are shown in Figs. 15 and 16 for symmetrical and asymmetrical weld descriptions, respectively. FE BEAM simulated propagation time is 19 µs and the rheological viscosity model used is Maxwell’s.

Fig. 15 shows an excellent agreement for both beam deflection and energy loss between the two models on the symmetrical weld. At a depth of 30 mm in the weld part, the maximum relative difference is less than 8% and its average is 3%. Moving on the asymmetric description, Fig. 16 shows an anomalous concentration of the energy simulated by RAY BEAM in comparison with FE BEAM with a maximum deviation of about 30% at -8 mm and an average deviation of 7% at a depth of 30 mm. These disparities correspond well to the location of strong gradients in the orientation properties of the anisotropy into the welded part. This observation of caustic area is similar to the phenomena observed in the previous test cases with Gaussian variation of the grain orientation along the z-axis.

It is worth noting that most of the discrepancies observed are mainly on the left-hand side of the asymmetrical weld. Hence, comparable results for the incident field estimation would be obtained if one attempts to determine the response of a defect located on the right side of the weld. By contrast, a significant bias is likely to occur with RAY BEAM if the crack is located on the left side. Assessing these biases of RAY BEAM with FE BEAM is a major asset in these particular testing configurations. It should be mentioned that CPU calculation times of the two modelling approaches are of the same order of magnitude in 2D, i.e. of the order of a second or a minute depending on the frequencies considered. These computational performances open the possibility of using FE BEAM to perform the full NDT study when the hypothesis of a 2D model is relevant.
<figure><img src="image_14.png" /><figcaption>Fig. 12. Maximum displacement modulus in time at the right side of the region of interest with respect to depth for the two simulations models in oblique incidence: at left the linear variation and at right the local Gaussian variation.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13. Analysis of ray tracing according to the incidence of the wavefront for the proposed Gaussian orientation variation highlighting folding phenomena when the grains are oriented at 60°.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. Symmetrical (left) and asymmetrical (right) Ogilvy's descriptions considered. Delay laws applied to the linear phased-array probe are illustrated by red bars above the probe. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

### 3D configuration

A welded component defined as a cylindrical extrusion of a 2D CAD profile, inspected with a TRL (Transmit–Receive-Longitudinal) probe, is

finally considered (see Fig. 17). The inclination of the probe elements and the curvature of the part induce 3D effects that must be taken into account by the simulation. The roof angle of 5° for the TRL probe
<figure><img src="image_17.png" /><figcaption>Fig. 15. Normalized maximum displacement modulus in time for RAY BEAM and FE BEAM simulations on the symmetrical weld description (left), extraction at depth z = 30 mm into the welded part (right).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Normalized maximum displacement modulus in time for RAY BEAM and FE BEAM simulations on the asymmetrical weld description (left), extraction at depth z = 30 mm into the welded part (right).</figcaption></figure>

is designed to naturally focus on the inspection plane at a depth of approximately 40 mm. Delay laws corresponding to an L45 deviation in a homogeneous specimen are applied to focus on the backwall in this plane approximately 5 mm to the right of the centre of the weld in the lateral position. An equivalent immersion configuration was assumed with density $ρ = 1.18 \text{ g cm}^{-3}$ and longitudinal wave velocity $v_L = 2.34 \text{ mm s}^{-1}$ in the coupling fluid (the same as for a Rexolite wedge). Both symmetrical and asymmetrical variations of Ogilvy's descriptions of the weld material are studied. The simulated ultrasonic field is extracted on a rectangular area at a depth of 30 mm inside the weld, with dimensions of 30 mm in the cylinder axis (lateral positions) and 40 mm in the extrusion axis perpendicular to the inspection plane (extension positions). Note that for this cylindrical 3D case, the lateral axis is now along y and the extension along x. Propagation simulation time is limited to 18 µs in FE BEAM to avoid contributions related to the reflection on the backwall, while only the direct wavefront contributions of the compression wave are considered in RAY BEAM.

Again, the description of the weld material has a strong influence on the deflection of the beam energy. As shown in Figs. 18 and 19, the lateral position of the maximum amplitude is located at –2 mm from the central position of the weld for the symmetrical case, whereas a significant concentration of energy is observed at –7.5 mm from this axis for the asymmetrical case. The beam deflection with out-of-plane divergence remains well represented for both models, but the quantitative estimates of the amplitude differ significantly. Considering

FE BEAM as a reference, the error is again mainly located on the left side of the weld (see Fig. 20). Concerning the symmetric configuration, the deviations of approximately 13% observed near the left bevel interface are mainly due to two reasons: mode conversions at the internal interface between the weld and the base metal are not taken into account in this simulation, and surface wave phenomena are not properly modelled by ray approaches. Concerning the asymmetric configuration, the difference of almost 40% can only be explained by the presence of caustics already observed in the previous test cases.

## 4. Conclusion

The aim of this work is to study and compare the fast ray-based solver available in CIVA to a new, more robust hybrid FE solution for simulating ultrasonic fields radiated into welds for the purposes of NDT. The weld material is modelled by an anisotropic viscoelastic medium with a spatially variable anisotropic orientation according to the grain elongation axis. This type of representation aims to reproduce an average description of the variations in weld properties in order to model the coherent wavefront behaviour commonly used for UT imaging techniques. The introduction of a non-zero imaginary part for the components of the stiffness tensor allows to account for the effects of energy loss of the coherent wavefront due to scattering at grain boundaries. Then, the two modelling solutions are introduced. Both solutions can deal with wave propagation and interactions in
<figure><img src="image_19.png" /><figcaption>Fig. 17. 3D weld configuration with a TRL probe and the simulated Region Of Interest (ROI).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 18. Normalized maximum displacement modulus in time at depth z = 30 mm for the 3D configurations.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 19. Normalized maximum displacement modulus in time at depth z = 30 mm w.r.t. lateral position x (top) and extension position y (bottom) for symmetric (left) and asymmetric (right) Ogilvy's descriptions.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Fig. 20.** Errors between RAY BEAM and FE BEAM in the weld area with symmetric (right) and asymmetric (left) Ogilvy's descriptions. The grey-scale is normalized w.r.t. the maximum error observed in the asymmetric (resp. symmetric) case on top (resp. bottom).</figcaption></figure>

these smoothly inhomogeneous anisotropic media, including damping behaviours to describe the coherent wavefront. The proposed hybrid FE solution combines a source field calculation using a ray solution and a transient FE calculation for estimating the wave propagation in the weld component.

We then carried out validations and comparisons of the two models, gradually increasing their complexity. In the first simple simulations involving a homogeneous anisotropic medium, we demonstrated the consistency of these models by comparing them with the solutions expected in a plane-wave framework. We have shown how variations in grain orientation affect both the energy deviation and the estimated attenuation with specific effects due to anisotropy. In these simple scenarios, we obtained strong agreement between the two models, both consistent with the expected plane-wave behaviour. However, in a subsequent phase, we sought to highlight the limitations of HF approximations, which can introduce unpredictable errors, particularly around caustics. To do this, we established domains keeping a simple geometric shape (planar half-space), but we introduced variations in the orientation properties of the grains only along a direction perpendicular to the plane surface. Depending on the wavefront direction, the orientation variations led to more or less significant deviations. The discrepancies between the two solutions were attributed to both property gradients and the orientation of the main wavefront. Notably, when encountering “critical layers” in oblique incidence cases, the ray method struggled to accurately handle the observed caustic phenomena, resulting in pronounced deviations.

Similar, yet more challenging to predict, phenomena arise in weld inspection configurations closer to real-world applications. Anisotropy variations within the weld, influenced by welding positions, have been observed to cause significant deviations in the ray model, especially in regions where both specific wavefront directions and abrupt changes in grain orientation coincide. These discrepancies can result in the erroneous simulation of ultrasound indications for defects located in

these regions of the weld. To address such issues and ensure the accuracy of fast ray solutions in both 2D and 3D cases, reference solutions established with the assistance of the hybrid FE model can be employed, providing interpretation and dispelling uncertainties regarding the reliability of the results.

The current mesh strategy for the FE scene results in a high-order structured mesh with excellent performance on straight or V-shaped bevels. Special attention is now given to expanding the range of weld geometries that can be simulated. Our primary focus is to strike a balance between maintaining high performance and increasing flexibility to handle a wider range of realistic, complex bevel shapes (such as backwall caps or complex weld roots). In addition, while the currently implemented hybrid solution addresses the biases associated with the HF approximation of the radiated field solution, there is no current equivalent for simulating the ultrasound echographic response. We are actively working on this task and will provide updates in a future paper.

### CRediT authorship contribution statement

**Nicolas Leymarie:** Writing – original draft, Software, Methodology, Formal analysis, Conceptualization.
**Alexandre Imperiale:** Writing – review & editing, Software, Methodology, Conceptualization.
**Thibaud Fortuna:** Writing – review & editing, Software, Methodology.
**Edouard Demaldent:** Writing – review & editing, Supervision, Methodology, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### Acknowledgements

This study was partly funded by the European ADVISE project, which received funding from the Euratom research and training programme 2014–2018 under grant agreement no. 755500.

### References

[1] Chapman C. Fundamentals of seismic wave propagation. Cambridge University Press; 2004.

[2] Cerveny V. Seismic ray theory. Cambridge University Press; 2005.

[3] Ogilvy JA. Computerized ultrasonic ray tracing in austenitic steel. NDT Int 1985;18(2):67-77. http://dx.doi.org/10.1016/0308-9126(85)90100-2.

[4] Gengembre N, Lhémery A. Calculation of wideband ultrasonic fields radiated by water-coupled transducers into heterogeneous and anisotropic media. AIP Conf Proc 2000;509(1):977-84. http://dx.doi.org/10.1063/1.1306150.

[5] Lhémery A, Calmon P, Lecoeur-Taibi I, Raillon R, Paradis L. Modeling tools for ultrasonic inspection of welds. NDT E Int 2000;33(7):499-513. http://dx.doi.org/10.1016/S0963-8695(00)00021-9.

[6] Gengembre N, Lhémery A. Pencil method in elastodynamics: application to ultrasonic field computation. Ultrasonics 2000;38(1):495-9. http://dx.doi.org/10.1016/S0041-624X(99)00068-2.

[7] Spies M. Modeling of transducer fields in inhomogeneous anisotropic materials using Gaussian beam superposition. NDT E Int 2000;33(3):155-62. http://dx.doi.org/10.1016/S0963-8695(99)00036-5.

[8] Spies M. Semi-analytical elastic wave-field modeling applied to arbitrarily oriented orthotropic media. J Acoust Soc Am 2001;110(1):68-79. http://dx.doi.org/10.1121/1.1380440.

[9] Schmerr LW. Fundamentals of ultrasonic nondestructive evaluation: A modeling approach. Springer series in measurement science and technology, Cham: Springer International Publishing; 2016, http://dx.doi.org/10.1007/978-3-319-30463-2.

[10] Connolly GD. Modelling of the propagation of ultrasound through austenitic steel welds (Ph.D. thesis), Imperial College London (University of London) London; 2009.
[11] Kolkoori SR, Rahman MU, Chinta PK, Ktrentzbruck M, Rethmeier M, Prager J. Ultrasonic field profile evaluation in acoustically inhomogeneous anisotropic materials using 2D ray tracing model: Numerical and experimental comparison. Ultrasonics 2013;53(2):396-411. http://dx.doi.org/10.1016/j.ultras.2012.07.006.

[12] Virieux J. P-SV wave propagation in heterogeneous media: Velocity-stress finite- difference method. Geophysics 1986;51(4):889-901. http://dx.doi.org/10.1190/ 1.1442147.

[13] Fellinger P, Marklein R, Langenberg KJ, Klaholz S. Numerical modeling of elastic wave propagation and scattering with EFIT — elastodynamic finite integration technique. In: Honor and memory of Julius miklowitz, Wave Motion In: Honor and memory of Julius miklowitz, 1995;21(1):47-66. http://dx.doi.org/10.1016/ 0165-2125(94)00040-C.

[14] Komatitsch D, Tromp J. Introduction to the spectral element method for three- dimensional seismic wave propagation. Geophys J Int 1999;139(3):806-22. http: //dx.doi.org/10.1046/j.1365-246x.1999.00967.x.

[15] Bossy E, Talmant M, Laugier P. Three-dimensional simulations of ultrasonic axial transmission velocity measurement on cortical bone models. J Acoust Soc Am 2004;115(5):2314-24. http://dx.doi.org/10.1121/1.1689960.

[16] Chassignole B, Duwig V, Ploix MA, Guy P, El Guerjouma R. Modelling the attenuation in the ATHENA finite elements code for the ultrasonic testing of austenitic stainless steel welds. Ultrasonics 2009;49(8):653-8. http://dx.doi.org/ 10.1016/j.ultras.2009.04.001.

[17] Van Pamel A, Huthwaite P, Brett CR, Lowe MJS. Numerical simulations of ultrasonic array imaging of highly scattering materials. NDT E Int 2016;81:9-19. http://dx.doi.org/10.1016/j.ndteint.2016.02.004.

[18] Lhuillier PE, Chassignole B, Oudaa M, Kerhervé SO, Rupin F, Fouquet T. Investigation of the ultrasonic attenuation in anisotropic weld materials with finite element modeling and grain-scale material description. Ultrasonics 2017;78:40-50. http://dx.doi.org/10.1016/j.ultras.2017.03.004.

[19] Kalkowski MK, Lowe MJS, Barth M, Rjelka M, Köhler B. How does grazing incidence ultrasonic microscopy work? A study based on grain-scale numerical simulations. Ultrasonics 2021;114:106387. http://dx.doi.org/10.1016/j.ultras. 2021.106387.

[20] Dunlap M, Connolly G, Dobson J. Modeling and simulation of cast austenitic stainless steel with onscale. In: 2020 IEEE international ultrasonics symposium. IUS, 2020, p. 1-4. http://dx.doi.org/10.1109/IUS46767.2020.9251616, ISSN: 1948-5727.

[21] Huang M, Sha G, Huthwaite P, Rokhlin SI, Lowe MJS. Maximizing the accuracy of finite element simulation of elastic wave propagation in polycrystals. J Acoust Soc Am 2020;148(4):1890-910. http://dx.doi.org/10.1121/10.0002102.

[22] Huang M, Huthwaite P, Rokhlin SI, Lowe MJS. Finite-element and semi-analytical study of elastic wave propagation in strongly scattering polycrystals. Proc R Soc A 2022;478(2258):20210850. http://dx.doi.org/10.1098/rspa.2021.0850.

[23] Carpentier C, Nageswaran C, Tse YY. Evaluation of a new approach for the inspection of austenitic dissimilar welds using ultrasonic phased array techniques. In: Proc 10th ECNDT conference,ECNDT, moscow, 2010). 2010. p. 442-4.

[24] Jacob RE, Prowant MS, Hughes MS, Diaz AA. Modeling and simulation of austenitic welds and coarse-grained specimens: Part II. Tech. Rep. PNNL-32702, Pacific Northwest National Laboratory (PNNL), Richland, WA (United States); 2022, http://dx.doi.org/10.2172/1988074.

[25] Silk MG. A computer model for ultrasonic propagation in complex or- thotropic structures. Ultrasonics 1981;19(5):208-12. http://dx.doi.org/10.1016/ 0041-624X(81)90004-4.

[26] Moysan J, Apfel A, Corneloup G, Chassignole B. Modelling the grain orientation of austenitic stainless steel multipass welds to improve ultrasonic assessment of structural integrity. Int J Press Vessels Pip 2003;80(2):77-85. http://dx.doi.org/ 10.1016/S0308-0161(03)00024-3.

[27] Chassignole B, El Guerjouma R, Ploix MA, Fouquet T. Ultrasonic and structural characterization of anisotropic austenitic stainless steel welds: Towards a higher reliability in ultrasonic non-destructive testing. NDT E Int 2010;43(4):273-82. http://dx.doi.org/10.1016/j.ndteint.2009.12.005.

[28] Ledbetter HM, Austin MW. Texture in stainless steel welds: an ultrasonic study. J Mater Sci 1985;20(5):1720-4. http://dx.doi.org/10.1007/BF00555276.

[29] CIVA. EXTENDE, Inc., URL https://www extende com/civa-in-a-few-words.

[30] Gardahaut A, Jezzine K, Cassereau D. Paraxial ray-tracing approach for the simulation of ultrasonic inspection of welds. AIP Conf Proc 2014;1581(1):529-36. http://dx.doi.org/10.1063/1.4864865.

[31] Kravtsov YA, Orlov YI, Edelev MG. Caustics, catastrophes and wave fields. 2nd ed. 1993 édition. Berlin ; New York: Springer-Verlag Berlin and Heidelberg GmbH & Co. K; 1998.

[32] Carrascal-Manzanares C, Imperiale A, Rougeron G, Bergeaud V, Lacassagne L. A fast implementation of a spectral finite elements method on CPU and GPU applied to ultrasound propagation. In: International conference on parallel computing. Parallel computing is everywhere, vol. 32, Bologna, Italy: IOS Press; 2017, p. 339-48. http://dx.doi.org/10.3233/978-1-61499-843-3-339.

[33] Imperiale A, Demaldent E. A macro-element strategy based upon spectral finite elements and mortar elements for transient wave propagation modeling. Application to ultrasonic testing of laminate composite materials. Internat J Numer Methods Engrg 2019;119(10):964-90. http://dx.doi.org/10.1002/nme.6080.

[34] Imperiale A, Leymarie N, Demaldent E. Numerical modeling of wave propagation in anisotropic viscoelastic laminated materials in transient regime: Application to modeling ultrasonic testing of composite structures. Internat J Numer Methods Engrg 2020;121(15):3300-38. http://dx.doi.org/10.1002/nme.6359.

[35] Chassignole B, El Guerjouma R, Ploix MA, Fouquet T. Ultrasonic and structural characterization of anisotropic austenitic stainless steel welds: Towards a higher reliability in ultrasonic non-destructive testing. NDT E Int 2010;43(4):273-82. http://dx.doi.org/10.1016/j.ndteint.2009.12.005.

[36] Ploix M-A, Guy P, Chassignole B, Moysan J, Corneloup G, Guerjouma RE. Measurement of ultrasonic scattering attenuation in austenitic stainless steel welds: Realistic input data for NDT numerical modeling. Ultrasonics 2014;54(7):1729-36. http://dx.doi.org/10.1016/j.ultras.2014.04.005.

[37] Dorval V, Jenson F, Corneloup G, Moysan J. Accounting for structural noise and attenuation in the modeling of the ultrasonic testing of polycrystalline materials. AIP Conf Proc 2010;1211(1):1309-16. http://dx.doi.org/10.1063/1.3362219.

[38] Hosten B, Deschamps M, Tittmann BR. Inhomogeneous wave generation and propagation in lossy anisotropic solids. Application to the characterization of viscoelastic composite materials. J Acoust Soc Am 1987;82(5):1763-70. http: //dx.doi.org/10.1121/1.395170.

[39] Carcione JM. Wave fields in real media, volume 38: Wave propagation in anisotropic, anelastic, porous and electromagnetic media. third ed.. Amsterdam: Elsevier Science; 2014.

[40] Stanke FE, Kino GS. A unified theory for elastic wave propagation in poly- crystalline materials. J Acoust Soc Am 1984;75(3):665-81. http://dx.doi.org/10. 1121/1.390577.

[41] Stanke FE. Spatial autocorrelation functions for calculations of effective propaga- tion constants in polycrystalline materials. J Acoust Soc Am 1986;80(5):1479-85. http://dx.doi.org/10.1121/1.394403.

[42] Sha G. Analytical attenuation and scattering models for polycrystals with uniform equiaxed grains. J Acoust Soc Am 2018;143(5):EL347-53. http://dx.doi.org/10. 1121/1.5037619.

[43] Huang M, Sha G, Huthwaite P, Rokhlin SI, Lowe MJS. Elastic wave velocity dispersion in polycrystals with elongated grains: Theoretical and numerical analysis. J Acoust Soc Am 2020;148(6):3645-62. http://dx.doi.org/10.1121/10. 0002916.

[44] Guy P, Mascaro B, Darmon M, Leymarie N, Clausse B. Accurate experimental determination of the complex elastic tensor of anisotropic materials of unknown orientation using ultrasonic transmitted bulk waves. Bruges, Belgium; 2019, URL https://hal.archives-ouvertes.fr/hal-02316999.

[45] Chassignole B. Influence de la structure métallurgique des soudures en acier inoxydable austénitique sur le contrôle non destructif par ultrasons (Ph.D. thesis), Lyon, INSA; 2000.

[46] Ogilvy JA. Ultrasonic beam profiles and beam propagation in an austenitic weld using a theoretical ray tracing model. Ultrasonics 1986;24(6):337-47. http://dx.doi.org/10.1016/0041-624X(86)90005-3.

[47] Moysan J, Corneloup G, Chassignole B, Gueudré C, Ploix MA. Modelling welded material for ultrasonic testing using MINA: Theory and applications. AIP Conf Proc 2012;1430(1):1219-26. http://dx.doi.org/10.1063/1.4716358.

[48] Calmon P, Mahaut S, Chatillon S, Raillon R. CIVA: An expertise platform for simulation and processing NDT data. In: Proceedings of ultrasonics international (UI'05) and world congress on ultrasonics (WCU). Ultrasonics In: Proceedings of ultrasonics international (UI'05) and world congress on ultrasonics (WCU), 2006;44:e975-9. http://dx.doi.org/10.1016/j.ultras.2006.05.218.

[49] Deschamps G. Ray techniques in electromagnetics. Proc IEEE 762;60(9):1022- 35. http://dx.doi.org/10.1109/PROC.1972.8850.

[50] Ogilvy JA. A layered media model for ray propagation in anisotropic inhom- geneous materials. Appl Math Model 1990;14(5):237-47. http://dx.doi.org/10. 1016/0307-904X(90)90014-V.

[51] Jezzine K, Gardahaut A, Leymarie N, Chatillon S. Evaluation of ray- based methods for the simulation of UT welds inspection. AIP Conf Proc 2013;1511(1):1073-80. http://dx.doi.org/10.1063/1.4789162.

[52] Gardahaut A, Jezzine K, Cassereau D, Leymarie N, Iakovleva E. Advanced simulation of ultrasonic inspection of welds using dynamic ray tracing. In: Proc. 13th Int. Symp. Nondestructive Characterization of Materials (NDCM-XIII). 2013.

[53] Klimeš L. Phase shift of the Green tensor due to caustics in anisotropic media. Stud Geophys Geod 2010;54(2):269-89. http://dx.doi.org/10.1007/s11200-010- 0014-x.

[54] Vavryčuk V. Real ray tracing in anisotropic viscoelastic media. Geophys J Int 2008;175(2):617-26. http://dx.doi.org/10.1111/j.1365-246X.2008.03898.x.

[55] Vavryčuk V, Svitek T, Lokajiček T. Anisotropic attenuation in rocks: theory, modelling and lab measurements. Geophys J Int 2017;208(3):1724-39. http: //dx.doi.org/10.1093/gji/ggw476.

[56] Imperiale A. An energy preserving time scheme based on the mortar element method for effective transmission conditions between fluid and solid domains in transient wave propagation problems. Internat J Numer Methods Engrg 2023. http://dx.doi.org/10.1002/nme.7246.

[57] Ben Belgacem F, Maday Y. The mortar element method for three dimensional finite elements. ESAIM: M2AM 1997;31(2):289-302.
[58] Casadei F, Gabellini E, Fotia G, Maggio F, Quarteroni A. A mortar spectral/finite element method for complex 2D and 3D elastodynamic problems. Comput Methods Appl Mech Eng 2002;191(45):5119-48. http://dx.doi.org/10.1016/S0045-7825(02)00294-3.

[59] Hauret P, Le Tallec P. A discontinuous stabilized mortar method for general 3D elastic problems. Comput Methods Appl Mech Eng 2007;196(49-52):4881-900. http://dx.doi.org/10.1016/j.cma.2007.06.014.

[60] Patera AT. A spectral element method for fluid dynamics: Laminar flow in a channel expansion. J Comput Phys 1984;54(3):468-88. http://dx.doi.org/10.1016/0021-9991(84)90128-1.

[61] Demaldent E, Imperiale S. Perfectly matched transmission problem with absorbing layers: Application to anisotropic acoustics in convex polygonal domains. Internat J Numer Methods Engrg 2013;96(11):689-711. http://dx.doi.org/10.1002/nme.4572.

[62] Pierce AD. Ray acoustics. In: Pierce AD, editor. Acoustics: an introduction to its physical principles and applications. Cham: Springer International Publishing; 2019. p. 427-86. http://dx.doi.org/10.1007/978-3-030-11214-1_8.

[63] Gueudré C, Mailhé J, Ploix M-A, Corneloup G, Chassignole B. Influence of the uncertainty of elastic constants on the modelling of ultrasound propagation through multi-pass austenitic welds. Impact on non-destructive testing. Int J Press Vessels Pip 2019;171:125-36. http://dx.doi.org/10.1016/j.ijpvp.2019.02.011.