Optimal computation of integrals in the Half-Space Matching method for modal simulation of SHM/NDE in 3D elastic plate

Simulating structural health monitoring (SHM) or nondestructive evaluation (NDE) methods based on elastic guided waves (GW) is very helpful to handle their complexity (co-existence of several GW modes, frequency dependence of wavespeed) and to further design optimal methods of inspection offering high sensitivity to the sought flaws. The half-space matching (HSM) method has been established for the development of a model that hybridizes local finite element (FE) computations for GW scattering by a flaw, with a modal semi-analytical model for GW radiation and propagation in flawless plate-like structures. Highly oscillatory Integral formulae appear in the HSM method that radiate the scattered field away from the FE zone as the superimposition of modal contributions, which computation can be time-consuming. The present work is concerned with their optimal computation. Integral of this form can be efficiently computed under the far-field approximation but this classical technique fails at predicting accurately wavefields at relatively short distances (small number of wavelengths). The method developed herein relies on the complexification of the integrals to be computed and on specific deformation of integration paths in the complex plane, as detailed in the paper. It allows the evaluation of the integrals without approximation other than that of numerical quadratures, ensuring high accuracy while offering high computing performances. It indifferently applies in the far-field and in the near-field. The method of computation is validated by comparing its predictions with a reference solution of GW scattering. Its computational performances are also demonstrated, compared to those of the standard computation of the HSM integral formulae.


Introduction
Structural health monitoring (SHM) of composite materials is expected to become an alternative to some nondestructive evaluation methods (NDE) requiring periodical tests.For NDE or SHM of thin structures made of metallic or composite materials, elastic guided waves (GW) are of high interest as they propagate over long distances.Therefore, large parts that are monitored or under test can be scanned at the speed of elastic waves rather than at that of a robot or a human.The multimodal and dispersive nature of GW make them challenging to handle (to master, to interpret) without prior deep understanding of their complex radiation, propagation and scattering, all the more when anisotropic structures are concerned (e.g., composite plates).Needs for accurate and computationally efficient tools for simulating NDE/SHM based on GW led to the development of many modelling approaches [1].
2 Among the various problems to be solved for a full simulation of SHM/NDT experiments, that of GW scattering by an arbitrary defect in an unbounded domain is specifically considered here.For example, its solution is of great interest to understand which guided mode is best adapted to interact with a given kind of defects, to ease the optimal design of NDT/SHM methods that would consist of selecting this mode among all possible modes to ensure higher sensitivity to these defects.
A quantitative simulation of the scattering of guided waves by realistic defects can only rely locally onto a pure numerical scheme, all the more when the overall structure is a composite plate.For such a complex case, it is impossible to find an approximate solution.However, these numerical schemes are challenged when dealing with large structures.A natural idea to overcome these difficulties is to develop a hybrid method of calculation consisting in treating the scattering by (for example) a finite element (FE) method locally in a small zone surrounding the defect and the propagation from a source to this zone and from this zone to a receiver by another method of calculation of higher numerical performance optimally developed to deal with GW propagation in the undamaged plate.Different models have been presented in the literature with various strategies of hybridization (see [1] where some are reviewed).In a long term collaboration of some of the present authors, a specific hybrid approach is developed with the special constraint of obtaining results (scattering coefficients, signal received etc.) decomposed onto the modal basis of the uniform (flawless) guides involved.The case of unidirectional waveguides has first been treated resulting in a modal hybrid model using the semi-analytical finite element (SAFE) method for propagation in unidirectional guide portions of arbitrary section and FE computations in "boxes" to treat the scattering by one or several features such as defects, junctions, etc. [2][3][4], special (transparent) boundary conditions being developed [5] so that waves can escape FE boxes without artificial reflections.A simulation tool based on this modelling approach is available in CIVA1 .The model has been extended to deal with the scattering of obliquely incident GW in a composite on an infinite straight feature [6] or with pipeworks including elbows [7].One of the most interesting properties of these specific transparent boundary conditions is that they also allows us to express the scattering coefficients of the FE boxes (reflection from or transmission through) as modal coefficients.Therefore, the various scattering phenomena are quantified over the modal basis of the guides involved in the simulated configuration.
The natural and expected extension of this modelling approach is that consisting in hybridizing different methods to solve the case of scattering by a localized feature (inner defect, delamination, surface corrosion etc.) in a multi-directional waveguide such like a plate-like structure.For this, a modal semi-analytical solution for radiation by an arbitrary source in an isotropic [8] or in an anisotropic [9] plate of finite size can be used to propagate GW outside localized small regions where scattering by various features arise.In these regions, a FE code can be used to compute the scattering without any other approximation than that of numerical quadrature.This is made possible only on the condition of having a means of defining transparent boundary conditions to the finite element calculation zone.
Such boundary conditions adapted to solve the problem of local scattering within an unbounded anisotropic plate have been developed in the so-called half-space matching (HSM) method [10,11].The method consists of coupling a finite elements (FE) computation on a small domain surrounding the defect with analytic representations in four overlapping half-plates.The HSM method uses integral representation formulae, expressed as modal expansions of Fourier integrals in each half plate.The system of equations ensures the compatibility of representations formulae in the overlaps.It results in the definition of transparent boundary conditions expressed on the GW modal basis.The HSM method has been first successfully implemented but it appeared that numerical evaluation of these integrals required long computation times to ensure the expected accuracy.
Speeding up their evaluation is crucial for the further use of HSM method in hybrid NDE/SHM simulations.In Sec. 2 of this paper, after a brief recall about the HSM, the method for speeding-up the integral computations is described.The method relies on a deformation in the complex plane of a contour integration resulting from the complexification of the integral to be computed.In Sec. 3, the method is used to treat a problem for which a reference solution is available [12] (only in the isotropic case).Comparing results obtained by the HSM calculated with this method with the reference solution makes it possible to validate the very principle of its calculations.A discussion of the performances obtained is also made.We conclude the paper with a summary of the main results obtained and some prospects.

Theory
The present section is divided into two subsections.In the first, we give a schematic summarized description of the HSM method for elastic wave scattering.The full description is presented in chapter 3 of Tjandrawidjaja [13] .The oscillatory integral formulae to be evaluated are shown.In the second subsection, the method to compute efficiently these integrals is presented.

The half-space matching method (HSM)
The HSM method is a domain decomposition-like method using finite elements, with semianalytical representations to solve unbounded wave scattering problems.The general theory of the HSM method has been already detailed in the literature [10,11].The aim of this part is to highlight the key concepts of the HSM and to extract integral formulae to be computed optimally.
To solve the wave scattering problem defined in an unbounded domain with the HSM method, the first step is to reduce the FE domain to a parallelepipedic box containing the defect.The unbounded properties of the initial problem are kept by introducing some unknowns defined on four infinite strips (Fig. 1a).The displacement d j and the normal stress s j of the scattered field are the unknowns considered on the strips Σ j , j := 1..4.Therefore, the HSM method is a multi-unknown formulation using The FE unknown in the box, the displacement and the normal stress unknowns on each strip.The principle of the HSM method is to consider several representations of the displacement field: the FE representation in the box and semi-analytical representations in infinite (half) plates delimited by the strips.Some conditions must be satisfied to ensure that the different representations coincide where ❏ The compatibility conditions By construction, there are four infinite quarters of plates where two representations coexist (Fig. 2a demonstrates the potential pairs of unknowns that can be employed for reconstructing the field in one overlapping domain).It is enough to impose some compatibility relations on the boundaries to ensure the two representations coincide in each quarter of the plate.These conditions involve relationships between the unknowns d j and s j based on integral representation integrals.

❏ The transparent boundary conditions
In order to restrict the finite element computation within a limited volume surrounding the defect, it is necessary to establish transparent conditions on the artificial boundary.In the HSM approach, this is accomplished by defining Dirichlet-to-Neumann (DtN) operators, which connect the unknowns in the strips to the finite element unknowns located at the boundaries of the box, as shown in Fig. 2b.
The HSM approach relies heavily on integral representation formulae.Indeed the compatibility operators and the DtN operator used in HSM method depend on these formulae.The various strips serve as boundaries of elementary half-space domains.Figure 1b depicts the domain to which strip Σ 1 belongs.By using the integral representation formula, it becomes possible to reconstruct the wave field in the elementary domains based on the displacement and normal stress (d, s) on the strips.Integral formulae are obtained by combining a Fourier transform in the transverse direction (i.e., parallel to the strip) and a modal decomposition in the thickness.Tjandrawidjaja [13] presents the formulae in the following form: with A k (ξ; d; s), the modal contribution that can be computed using the bi-orthogonality of the where the subscript 'ˆ' indicates a 1D Fourier transform.ξ corresponds to the Fourier variable.The terms β ξ , Ûk (ξ; z) and Tk (ξ; z) are respectively a wave vector component, the normalized modal displacement and the normalized modal normal stress in the thickness, for a fixed value of ξ.In our developments, these terms are obtained using SAFE method [14].The normalization uses the Poynting vector computed by post-processing the SAFE results.Substituting the expression of A k (2) into Eq.( 1) gives: The observation point and the data on the strip position are denoted by x(x, y, z) and x'(x',y',z'), respectively.The integral kernels K DD and K DN are 3 by 3 matrice-valued functions, defined as: and These terms are repeatedly calculated in the HSM method.Specifically, the linear matrix term of the method necessitates the computation of these kernels for numerous combinations of degrees of freedom that belong to the strips.Computing K DD and K DN numerically can be computationally burdensome.To demonstrate this, let's take one term of the integral kernel K DN associated with a given mode (k fixed).We denoted this term by K, and it is expressed by: From ( 6), one can define the associated integrand as f : ξ → a(ξ, z)e i (β ξ x+ξy ) .
The following figures depict the integrand function f with respect to the variable ξ for various positions of x and y.As x or y increases, oscillations of the integrand become more pronounced.Specifically, when x is large, the oscillations tend to gather at two points that correspond to the wavenumber and its opposite value.Consequently, the accurate numerical evaluation of the integrand necessitates a very high discretization of the integration domain.This is the primary factor of high computational costs for building the HSM linear matrices.The complex number coefficients a(ξ, z) appearing in the expression of K ( 6) do not significantly affect the oscillatory behavior of the integrand, as shown by Fig. 4 displaying only the exponential part of the integrands, which present oscillations very similar to those shown on the previous figure for the whole integrands.Improving computational performances in the evaluation of this type of integrals can be studied on a kernel involving these oscillations but of simpler expression to ease the subsequent writing.
In what follows, we temporarily consider the kernel K ϵ , which corresponds to that obtained when applying the HSM method to the simple 2D isotropic acoustic case, for which a closed form solution exists involving Hankel functions: Here r and θ are the polar coordinates of a point (x, y).Moreover, the square-root is chosen such that the integrand is exponentially decaying.The acoustic isotropic kernel, represented by (7), behaves globally the same way as the 3D elastic kernels.To improve both the accuracy and the computational efficiency, we introduce a technique called contour deformation, which involves re-formulating the kernel in the complex plane.The subsequent section provides a detailed explanation of this method, which results in a less oscillatory integrand that is less challenging to compute, leading to the expected improvements.

Contour deformation in the complex plane
The contour deformation method exploits the analytic properties of the integrand and allows to replace the real path integral with an integral path in the complex domain, while ensuring that the integral value remains unchanged.This is justified by the Cauchy integration theorem which says that the integral of a function on a closed path within which the function is analytic is zero [15].The goal is to find a new path where the integrand is less oscillating.Such approach is used in [16] for a phase function of the type φ(ξ, θ) = ξ p with p a positive integer.
In our case, a special care must be taken to study the analyticity properties of the phase function ξ → φ(ξ, θ) and to the definition of branch cuts.
Let us first present the method for the simple case of the kernel (7).We are looking for a parameterized path C = {ξ(η) ∈ C \ η ∈ R} such that: Roughly speaking, identity ( 8) is a consequence of Cauchy's integral theorem as soon as the integrand e iφ(ξ,θ)r is an analytic function of ξ in the domain S between the real axis and the path C, and if it decays fast enough at infinity.In particular, the square root function ξ(η) −→ β 2 0 − ξ 2 (η) has to be analytic in S. We define the complex square root as follows: For z = ηe iθ with η ≥ 0 and 0 so that the branch-cut is the positive real axis and that Using this definition of the function √ z, the function ξ(η) −→ β 2 0 − ξ 2 is analytic outside the branch cut iR ∪ [−β 0 , β 0 ] (represented by the red cross in Fig. 5.Note that by (10), the imaginary part of β 2 0 − ξ 2 is strictly positive for |ξ| > β 0 , which ensures the convergence of the integral (7).On the other hand, for |ξ| < β 0 , we want the complex square root to coincide with the usual real square root.We check that this corresponds to the bold black line on Fig. 5.This black line and the branch-cuts determine the permissible region where the new path C can be defined.
For the sake of simplicity of the method, the contour C is deduced from the real axis by a rotation of angle α, as illustrated in Fig 5: This angle α has to fulfill two conditions.First, due to the branch cut and to the position of the black line as explained above, we have to take 0 < α < π/2.On the other hand, the integrand as to decay at infinity on C.This is achieved if ℑm(φ(ηe −iα , θ)) > 0 for large η ∈ R, leading to the condition tan α | tan θ| < 1. (11) Graphically this condition means that the path C (in orange) has to between the real axis and the dashed line (Fig 5).Conversely, if we fix α, this condition defines the domain of the points (r, θ) where the deformation can be applied.Choosing α in accordance with the validity condition (11), the following illustrations depict the impact of complexification on the integrand.By deforming the contour, it becomes possible to reduce drastically the oscillations of the integrand, resulting in a more precise numerical evaluation that requires fewer integration points than previously needed.Figure 7 shows the kernel K ξ computed over a 2D zone, before and after complexification.Here, N f ou denotes the number of intervals used integral discretization.For N f ou = 200, the map exhibits artifacts that disappear after complexification.Along the vertical line x = 8, the relative error is 20% for N f ou = 200, whereas it reduces to less than 10 −5 % for N f ou = 100, after complexification.After complexification, some artifacts however appear (as seen on Fig. 7 right).This occurs here at points where Eq. ( 11) is not satisfied.When using the HSM method for a FE box inside a whole plate, field reconstruction is ensured by using all the four strips and this kind of artefacts disappears.On Fig. 8, the accuracy of the complex path deformation is illustrated for various numbers of integration points.We have plotted the relative error of this new method along the vertical line x = 8 with respect to the number of integration points N f ou .

Results
In the present section, one wants to validate the method of calculation of integrals appearing in the HSM method, based on the complexification as explained in the previous section Moreover, as our actual goal is to apply the HSM method to 3D elastic guides, results for such configuration are studied, being more representative, thus more demonstrative.One exact solution of the literature is used in what follows as a reference solution, allowing us to determine the accuracy of our numerical methods by comparing its predictions with exact results.The reference solution is presented in the first subsection.In the second, the various results (reference, computed by means of the proposed method) are compared for validation.Numerical performances relatively to the integral computation without complexification are discussed.

Reference solution
Lamb waves can be decomposed as sums of Hankel functions (see for example Achenbach and Xu [17]).This decomposition was used by Diligent et al [12] in the development of an exact solution for the scattering of the fundamental symmetric Lamb wave S 0 by a through thickness hole in an isotropic plate.The contribution of the symmetric Lamb mode in the exact solution is expressed by: u r (r, θ, z) = n=0,1.. m=0,2,..
u θ (r, θ, z) = n=0,1.. m=0,2,.. with V n (z) = s 1 cos(p n z) + s 2 cos(q n z), W n (z) = s 3 sin(p n z) + s 4 sin(q n z), s 1 = 2 cos(q n h), where 2h, k n , ω, c L and c T denote the thickness, the wavenumber, the angular frequency, the longitudinal and transversal bulk wave velocities, respectively.The reference solution can be expressed as the linear combination of contributions from all harmonic modes.By fixing the mode n = 0 and the harmonic m = 0, the solution satisfies the scattering problem.It can be used as a reference solution for the validation.With the time convention e −iωt , the function H 0 is the zero order Hankel function of first kind.Figure 9 shows the 3-components of the reference solution plotted in the cartesian domain.Figure 11 represents the reconstructed domain with the initial integral for a number of integration points N f ou = 30.The relative error obtained is 5.49%, while by using the the new formulation with a complexification angle α = π 60 , the relative error is 1.18% (Fig. 12).The impact of complexification on the relative error with respect to the discretization N f ou is depicted in Fig. 13.The figure demonstrates that the new formulation is more stable than the initial one.Specifically, the convergence of the complexification is already achieved for N f ou = 20, whereas for the classical integral evaluation, the relative error varies randomly.Figure 13: Average of the relative error on the domain reconstructed before ( α = 0 ) and after (α = π 60 ) the complexification.

Discussion
Hybrid modelling is promising for the development of simulation tools most helpful for designing optimal SHM/NDE methods.As far as elastic guided waves are concerned, simulation tools providing results in terms of modal amplitudes are desirable as practitioners generally interpret results as modal contributions.Developing such hybrid models requires both different submodels to treat different subproblems and means for them to exchange intermediate results.The HSM method offers a way of doing it with many advantages, the two main ones being that i) it is compatible with a modal description of scattered fields and ii) it allows to use numerical boxes of minimal size around the scatterers.The HSM method yields the numerical evaluation of integral formulae.Their oscillating integrands can lead to lengthy computation times to ensure the required accuracy obtained with sufficiently small variations of the integration variable, that is to say, sufficiently fine discretization.In the present paper, one way of reducing these computation times has been developed and proved to lead to accurate predictions, through the complexification of integral formulae.Another way of doing it but of restricted applicability is to use classical far-field approximation, which reduces the integral computation to the evaluation of the integrand around the path of stationary phase.This has been also implemented in the framework of the HSM method applied to elastodynamics, but not discussed herein and it constitutes the most efficient way of computing the integral formulae of the HSM method when it applies (far-field).
These means for evaluating HSM integrals pave the path for the development of a performant hybrid model for simulating a whole SHM/NDE experiment.Moreover, it is interesting to notice Note that the theoretical model of hybridization presented herein has been derived without any assumption about elastic symmetry of the part.However, for validating the mathematical method for computing the various integrals involved in the formulation with higher accuracy and in a shorter time, comparisons with published results of a semi-analytical solution have been made.Such solutions are not available in anisotropic cases so that presented results concern the case of an isotropic plate.In practice, our numerical developments are made to deal with anisotropic plates, as long as modal solutions can be computed using the SAFE method.
With the possibility to compute integral formulae of the HSM both accurately and efficiently thanks to the complexification procedure or to the far-field approximation (when it applies), as explained in the present paper, computations hybridizing FE to deal with the scattering by complex defects and the GW pencil method [8,9] will gives the possibility to optimally simulate complex NDT/SHM configurations involving GW, allowing the modal interpretation of complex results that could not be interpreted otherwise.
(a) u x (b) u y (c) u z

Figure 9 :
Figure 9: Real part of the reference solution S 0 plotted for ω = 2π MHz, h = 0.5 mm.

Figure 10 :
Figure 10: Configuration for the validation.

Figure 11 :
Figure 11: (a) the displacement u and (b) the error u ref − u x for N f ou = 30.

Figure 12 :
Figure 12: (a) the reconstructed u and (b) the error u ref − u x for N f ou = 30 and α = π 60 .

20th
Anglo-French Physical Acoustics Conference IOP Conf.Series: Journal of Physics: Conf.Series 2768 (2024) 012004 IOP Publishing doi:10.1088/1742-6596/2768/1/01200413 that the technique of contour deformation can be useful in other contexts of wave modelling where similar integral formulae with oscillating integrands are involved.