Custom Transient Finite Element Method and Ray Tracing Hybridization Strategies for Ultrasonic Testing Modelling

The transient spectral element method (SEM) is a specific high order finite element method that is particularly accurate and fast with a low memory footprint, hence enabling 3D ultrasonic testing simulations on a standard PC. However, particular care in its settings and hybridization with a semi-analytical solution remains a major challenge when seeking for both practicality and performance. A natural hybrid strategy consists of coupling a field calculated by ray tracing in the defect-free object to a SEM subdomain wherein the perturbation is enclosed, then synthesizing the signal with a reciprocity argument. The difficulty is to define the suited coupling strategy when the flaw response interacts with the back wall. In a highly heterogeneous medium, it may even become necessary to widen the numerical domain to the entire thickness of the inspected object and, therefore, to ensure the performance of the SEM through a custom discretization strategy. Here we present different customized uses of the SEM for practical ultrasonic inspection applications and discuss how they complement the ray tracing solution.


Introduction
Modelling provides a valuable aid to understanding, evaluating and qualifying ultrasonic testing (UT) techniques. It relies on both the evolution of physical models and their interfacing with software to make representative simulations of various industrial applications available to non-destructive testing experts. The strategy developed at the CEA List through the CIVA software consists in sequencing the overall UT modelling as a simulation chain to treat the source (emitter), the propagation (specimen), the interaction (flaw) and then the response (receiver). Coupled with a calibration protocol and a reciprocity-based signal synthesis formalism [1], it facilitates the constitution of a library of sensors, geometries, materials and defects.
Semi-analytical approaches, in particular asymptotic ray tracing models for propagation and diffraction models for defect response, address most industrial needs and allow simulation of NDT studies on standard PCs. However, they are based on far-field and high-frequency assumptions, simplification of field components (modes) and interactions between defects, etc. Measuring the impact of taking into account or violating these assumptions becomes a challenge to qualify these models on new inspection configurations, or to complete them. This is why CEA List also proposes various forms of hybridization between the asymptotic ray-tracing model and the finite element method (FEM). Initiated with the introduction of 2D partner models [2], a 2D/3D FEM model with various hybridization strategies has since been developed. Section 2 introduces the main features of the method, section 3 presents the coupling for a flaw-restricted FEMbox while section 4 discusses its extension to the specimen.

The FEM model
The FEM solution aims at providing access to full 2D and partial (few shots or scans) 3D NDT studies on a standard PC. It is based on the transient spectral element method SEM [3] that combines a time-explicit scheme with a small memory footprint and a high degree polynomial representation of the simulated field. Furthermore, a domain decomposition method (DDM) by mortar elements [4] reinforces the stability of the calculation and facilitates the introduction of new variants of the FEM propagators: fluid, homogeneous or continuously variable solids, viscous, thin layers, absorbing layers... Particularly in an NDT context, the disconnection of a mortar element represents an idealized infinitely fine notch, modelled by a free boundary condition between two domains. We refer to [5] [6] [7] for more details on the model.

Hybridization
For each emitting element and inspection mode combining Longitudinal (L) or Transverse (T) waves, the ray field is calculated in the defect-free part at each time step of the FEM solver and at each degree of freedom (DoF) of a loading area. It is assimilated to the source of the FEM calculation. The FEM solution field is saved at each time step and at each DoF of a convolution area . After FEM calculation, the ray field for all receiving channels is calculated in the sound part at FEM time steps and convolution DoF. Then the reciprocity formula is applied to synthesize the signal ( ): where and correspond to the displacement fields with and without defect for the emitter and receiver, respectively ( and are the associated stress tensors whose projection along the local normal is calculated). This is repeated for each emission channel. A first distinction in the hybridization mode lies in the choice of the calculated quantity of the FEM: the total field = or diffracted field = − where is the emission ray field without defect in the FEM domain. In the total field formulation, the coupling loading and convolution zones are located on the contour of the FEM simulation box, excluding the absorbing layers for which = − is imposed. Using the diffracted field formulation, the coupling zones are limited to the disturbance of the healthy part, in practice the flaw surface. It allows reducing the time window that is imposed according to the timing of the ray energy measured at the coupling zones and a buffer defined empirically according to the size of the defect or the FEM box. Furthermore, note in (1) that imposing ⋅ = 0 on an idealized notch reduces the memory load. Hereafter we explain the choices made in the presented work to define the FEM scene according to the computational configuration as illustrated Fig.1.

Description of the flaw
The first approach is to restrict the simulation of the FEM to a small sub-domain, as close as possible to the defect. This defect is described as a 2D geometry extruded over a limited distance and then placed freely in the specimen to simulate a 3D defect. A dedicated GUI allows editing of the defect or a defect network. Shown in Fig.2, this editor allows to introduce new vertical or horizontal edges into the FEM box (each cell corresponding to a FEM sub-domain), to declare certain edges of the grid as cracks (by deleting the associated mortar) and to move the internal nodes of the grid to get closer to the desired shape. An inclusion of air or any vacuum-like material is achieved by activating a closed loop of edges in the grid. Among other options, a color palette alerts the user to a strong local deformation of the grid that will result in an additional cost to the FEM calculation. The FEM mesh is obtained by refining the grid according to the rule of 2 edges of order 4 (i.e. an average of 8 DoFs) per wavelength, which is sufficient for the SEM.

Embedded defect
The FEM calculation is performed in diffracted field: the loading and convolution areas are restricted to the surface of the defect (Fig.1, (b)). Under this formalism, the inclusion of a non-vacuum medium would require extending the convolution to the volume of the inclusion without benefiting from the free boundary simplification, which could be less efficient than a full field coupling (Fig.1, (a)). Before adding the absorbing layers, the FEM-box is automatically adjusted to impose a margin of one wavelength between the defect and absorbing layer (gray box in (b.2), Fig. 2).

Border defect
Taking into account the interaction between a border defect and the specimen edge requires a more advanced approach. In the first instance, the same scheme is applied, this time including a free boundary condition at the bottom of the FEM box. The simulation takes into account the bounces of the ray field on the bottom of the part, as well as the possible bounces of the FEM surface waves along the crack. However, the incident surface waves on the bottom of the part are not simulated by the ray calculation, and compliance with the far-field criterion between the bottom and the defect is not ensured.
One option is to extend the loading and convolution area to the bottom of the FEM box so that these contributions are taken into account by the FEM calculation. From a coupling point of view, this requires removing the ray contributions after bouncing in the FEM box so as not to duplicate them. On the other hand, the signal now takes into account the defect and part of the workpiece backwall.
These two coupling models, illustrated Fig. 3, are called "partial" and "full" respectively in CIVA. The former is sufficient to simulate a specular inspection, while the latter is necessary when the surface waves on the bottom of the part are critical.

Hybrid scheme for part and defect
The hybrid approach restricted to the defect provides elements of understanding and analysis, and increased confidence on certain defects (e.g. 60° T inspection on a complex defect) but the complexity of the coupling protocol rapidly limits the method's extension capabilities (material inclusions, complex heterogeneous media…). This motivates the introduction of a more generic hybrid mode including the part and the defect, which is possible due to the performance of the SEM. Initially introduced at CEA List to deal with the UT of laminated composite materials, such a process is now generalized to 2.5D configurations: 2D extruded specimen, contact sensor and defect constrained by this extrusion orientation.
The FEM computational domain corresponds to the truncated workpiece beyond the area of influence of the ultrasound beam (Fig. 1, (c)). Less efficient than the hexahedral approach (by default of order 4), a prismatic triangle-based mesh (of order 2 in the 2D triangle plane mesh [8] and 4 along the extrusion) now facilitates the automatic CAD meshing process, see Fig. 4. The integration in a commercial release of CIVA is underway and R&D work to lift the restrictions on the orientation of the sensor and the defect will be carried out in the coming years.

Numerical comparison
As instructional examples, several simulations are performed for different incidence angles of an immersed circular probe (12.7 mm diameter, 2 MHz signal with a -6 dB bandwidth of 50%). The sample is a steel plate ( = 5.9 mm/µs, = 3.2 mm/µs, = 7.8 g/cm 3 ) of 30 mm thickness including a 5 mm high crack at the bottom. The different orientations of the transducer are chosen to favour the excitation of incident wave modes that are more or less appropriate to the simplifications carried by each hybrid models, as illustrated in Fig. 5. The "partial" model is sufficient to give a good approximation of the signal at L45°, whereas a discrepancy becomes visible at T60° due to surface wave conversions on the defect in the vicinity of this critical angle (from the point of view of the diffracted field), which are better taken into account by the "full" modelling approach. However, for the T30° incidence close to the longitudinal critical angle (connection zone between the shear wave and the head wave), the "partial" and "full" simulations suffer from the ray approximation in the volume of the part, which explains the significant difference with the reference FEM simulation using the whole part including the defect. As an indication, the calculation times obtained with the L45° inspection were approximately 2 sec. in partial, 5 sec. in full and 20 sec. with the whole part.