The Fast Computation of Eddy Current Distribution and Probe Response in Homogenized Composite Material Based on Semi-Analytical Approach

. Due to the excessive use of composites in the industrial ﬁeld, many numerical modeling approaches dedicated to the characterization of such complex material by means of Non Destructive Testing Techniques were developed. In this paper, we present a numerical model dedicated to simulate the inspection of unidirectional Carbon Fiber Reinforced Polymer using Eddy Current technique for detecting ﬁber disorientation. A semi-analytical model based on a modal approach is developed for the fast computation of quasi-static ﬁeld induced by an arbitrary 3D Eddy Current probe in the material. Because of the high anisotropy and strong heterogeneity of such material, a prior phase of homogenization is assumed and the material is then considered as homogeneously anisotropic. The modal approach consists in resolving Maxwell’s equations in the Fourier domain. Therefore, the electromagnetic ﬁeld is expressed as a sum of eigen-modes. To take into account the wave propagation through the multilayered structure and boundary conditions at each separating interface, a stable and recursive scattering matrix algorithm has been implemented. The impedance of the probe is computed analytically using Auld’s formula in orders to identify the main orientation of the ﬁbers in the inspected zone. For numerical validation, simulated data provided by the model are compared to ﬁnite element data.


Introduction
The use of unidirectional composite materials as Carbon Fiber Reinforced Polymer (CFRP) in manufacturing aeronautical structure has been significantly increased during the last decades thanks to their lightness and stiffness characteristics. Therefore, a wide variety of NDT techniques have been used, including ultrasonic testing, infrared thermography testing, radiographic testing, shearography testing, and Eddy Current Testing [1], for the characterization and the detection of defects such fiber disorientation in intermediate layers, delamination, out-of-plan waviness or fiber breakage.
Despite of the low electrical conductivity along carbon fibers with respect to metallic structure, eddy currents remain sensitive to the abrupt variation of the conductivity and sufficiently enough for detecting eventual damages. In a favorable manner, this method has an ability of fast and contact-free inspection of defects in conductive materials and also the global characterization of the sample as the estimation of the bulk conductivity using a multifrequency analysis and the identification of fiber orientation [2]. Therefore, many research works show interest Send offprint requests to: a E-mail: houssem.chebbi@cea.fr in developing inspection procedure using new specific sensors and signal processing [3,4]. Moreover, CFRP is a nonhomogeneous and highly anisotropic conductive material due to the random positioning of the carbon fibers which it challenging for modeling. To overcome this problem, homogenization method where proposed in [5] to determine an equivalent conductivity tensor. Besides, new modeling approaches dedicated to the analysis of the response of such complex material to a EC probe excitation have been developed: Finite Element (FE) models were proposed as in [6,7] for the characterization and flaw detection in homogenized composite materials but still need costly computing resources. Other models were introduced in [9, 10] based on integral equations and Green's Dyads. In addition, analytical models were presented in [11] to deal with the response of such material to an incident plane wave.
In the present paper, we present a Semi-Analytical (S-A) model based on a modal approach dedicated to the fast computation of quasi static fields and EC density distribution induced by an arbitrary EC probe. Contrarily to the S-A model in [20] where the second order potential formalism where used, our S-A method consists in resolving Maxwell's equations in the Fourier domain in order to obtain an algebraic form. Hence, the Electromagnetic (EM) field is expanded in each layer of the laminate as sum of eigen-modes. Applying boundary conditions (BC) leads us to estimate the unknown coefficients of the modal expansion. In practice, we only have access to the response of the probe to the presence of the plate or the eventual defect, thus we compute the variation of its impedance using Auld's formula [12] involving the interaction between the incident field from the coil and its interaction with the anisotropic material.
To deal with the stratified media, we made use of a S-matrix formalism, used generally in optics for numerical modeling of multilayered diffraction gratings [13,14], in which the multilayered structure is considered as a black box taking into account the propagation of the fields through the layers. Face to the physical complexity of the problem due to the high anisotropy and non-homogeneity of the material, a prior phase of homogenization is assumed and the material is considered as anisotropically homogeneous layer by layer while the rotation of the fibers is considered in the conductivity tensor.
This paper is organized as follows, first we present the formalism used for the S-A modeling approach and we illustrate the keystone of the modal decomposition of the EM field in both isotropic and anisotropic media. Then we express the BC and the S-matrix formalism used for stratified media. Next, the computation of eddy current distribution and probe response are explained. Finally, the comparison between the numerical results provided by the model and Finite Element data is presented to demonstrate the effectiveness of the method.  As described in Fig.1, a CFRP structure is made of of several unidirectional plies with different orientation of around 125 µm of thickness. Each ply is composed of carbon fibers embedded in an electrically non conductive polymer matrix. Since the fibers are in reality not perfectly aligned as displayed in Fig.2, the contact between them give rise to a transversal and cross-ply conductivities. Depending on the volume fraction and the orientation of the fibers in the ply, the electrical conductivity along the fiber varies between 5.10 3 and 5.10 4 S/m, and between 1 and 5.10 2 S/m in the transverse and cross-ply direction [5]. The conductivity effects are therefore introduced by (σ l , σ t , σ n ). The rotation of the fibers is referred by the angle θ introduced in the rotation matrix R: Thus, the conductivity tensor in each ply is obtained as follow:σ = R(θ).

Physical equations
We suppose that the field components have e iωt time dependency and assuming the absence of any external electric source and external charge density, we have J s = 0 and ρ = 0. Thus, we focus on Maxwell-Faraday law for induction and Ampere's law equations. We introduce the permeability and the permittivity of vacuum respectively µ 0 and ε 0 : Therefore we define: The permittivity tensor is denoted byε r refers to the uniform electric polarization of the material (in our case we chooseε r =Ī but it is not required in general). The relative permeability is a constant µ r = 1.

Modal decomposition in isotropic media
In isotropic media whereσ is reduced to a constant (σ = 0 in air), so that from Eq (3) and Eq(4), the tangential components of the Electric field donated by E x , E y and Magnetic field denoted by H x , H y can be expressed in function of two potential E z and H z so called the TE/TM decomposition: At this stage, we suppose the variable separation f (x, y, z) = f (x, y)e −iγz , thus, one can obtain: Where k 2 c = k 2 µ r ε r and Z 0 = µ0 0 . The two potentials E z and H z (referred by φ) verify Helmholtz's equation to be resolve:

Modal Decomposition in anisotropic media
In the general case whereε takes the Eq(3) and Eq(4), we obtain: Contrary to the isotropic case, the TE/TM decomposition cannot be adopted since the tangential components of the fields are no longer decoupled. From Maxwell's equations, we obtain: where L EG and L GE is differential operators defined by: Applying the operator ∂ z on Eq (12) and Eq(13) we obtain two decoupled differential systems verified by the transversal components of the electric field in one hand and the magnetic field in the other hand. we obtain: Finally we can recompute the E z and H z from the last two equalities of Eq(10) and Eq(11).

Fourier Transform and Truncating
To obtain an algebraic form, we aim to resolve Maxwell's equation in the Fourier domain. Thus, we apply 2D discrete Fourier Transform (FT) along Ox and Oy donated by: The M u and M v denote the number of modes used to represent a wave function in the Fourier domain. Numerically, each component is a matrix of dimension L, with In the Fourier domain, the solution of Helmholtz's equation can be expressed as sum of Fourier basic functions: These expansions are known as Rayleigh expansions. From Eq(8) we obtain the modal expansion: Adopting the Einstein convention we obtain: In the above compact relations, the superscripts ± correspond to wave that propagate or decay in the positive (+) or negative (-) Oz direction respectively. The unknown vector coefficients A ±T E et A ±T M are determined by Boundary Conditions. On the other hand, the eigenvalues system (Eq (16)) provides two eigenvalues ± √ λ 1 et ± √ λ 2 corresponding to the direction of the propagation of the wave, and two associated eigenvectors v 1 and v 2 . So E x and E y are a linear combination of eigenvectors. In a vector writing: The eigenvalues are used to calculate the wave attenuation factors φ ± i (γ). From Eq.13 we reconstruct the H t components: Finally the modal decomposition in the anisotropic media is defined by: Where the modes are:

General configuration
We consider the configuration in Fig.3 which consists of a multilayer stack with N layers, labeled by p=1,2,...,N of the wave number k cp and thickness e p . each layer is bounded by the p th and (p + 1) th interfaces in which the eigenmodes are denoted by Ψ ± p . At the p th interface, we have the outgoing waves corresponding to the coefficients (a +up p ,a −dn p ) and the incoming waves corresponding to the coefficients (a −up p ,a +dn p ). the superscripts -and + refers to transmission and reflection respectively. Moreover, the superscripts (up) and (dn) refer to the upper and downer coefficients regarding the p th interface. Assuming that the 0 th and the N th layers are in air, a − 0 results from the incident field due to the coil and a + N = 0.

Boundary Conditions
In this paragraph, we express the boundary conditions (BC) at the first separating interface " air-plate " to determine the unknown coefficients of the modal expansions.
In the absence of any surface current (J s = 0), the BC at a separating interface are defined by the continuity of the tangential components of the EM field: This leads us to the equality between the two modal expansions. We consider the relationship between outputs and inputs coefficients: In the case of semi-infinite case A + 1 = 0. A − 0 is computed using the relationship: where Γ − 0 is the part of the mode Ψ − 0 that corresponds to the magnetic field and the incident EM field provided by the coil in air E inc z , H inc z is calculated by means of analytical model already implemented into CIVA software [15].

S-matrix algorithm
The S-matrix formalism [16,17] is now well established for the study of modulated or planar stratified media. It is not affected by numerical instabilities linked to the number of the thickness of the layers. This is precisely because the boundary conditions are writing by dividing the waves according to the propagation medium. For a given structure as shown in Fig.3 a (4x4) matrix connects the incoming waves with the outgoing waves. Globally the sought global S g matrix is such as: The global S g matrix is obtained through classical recursion formulas between intermediate S p matrices at each separating interface. The intermediate scattering matrix where In a recursive way, all the intermediate S-matrices are concatenated taking into account the relationship between the coefficient describing the wave propagation through the layers. So we can write: where φ ± = e ±ikc p γep is th attenuation factor through the p th layer and e p is the thickness inside the material.

Eddy Current Density and skin effect
In metallic structure, Eddy currents are closed loops of induced current circulating in planes perpendicular to the Fig. 4. Illustration of the field propagation in the p th layer magnetic flux. However in CFRP, Eddy Currents have a different behaviour: they lengthen along the fibers due the high anisotropy and its decreasing depends on the nature of the plies and how they are stacked. Eddy currents remain travelling parallel to the coil's winding and its flow is limited to the area of the inducing magnetic field [18,19]. The depth that eddy currents penetrate into a material known as skin effect is affected by the frequency of the excitation current and the electrical conductivity and magnetic permeability of the specimen. To keep the δ reasonable with respect to the low values of σ zz , we work in the frequency range 0.8 -10 MHz. Although, δ depends also on the number and the sequence of plies. Thus, at fixed depth 'p' we consider:

Impedance of the probe
We are interested here in the calculation of the impedance of the probe in the presence of anisotropic flawless plate.
In particular, we focus on the variation of the impedance ∆Z. For that, the method we chose is based on the reci- procity theorem proved by Lorentz and adapted by Auld [12]. It has been proved that the variation ∆Z of the coil impedance is given by the following formula: (42) Where I 0 is the driven current intensity circulating in the coil. E i and H i are respectively the incident electric field and incident magnetic field producing by the coil in air. Numerically, they are computed by means of a module already implemented into CIVA software based on analytic approach (Dyades and Deeds) [15], while E and H are respectively the total electric and magnetic field resulting from the interaction with the anisotropic multilayered structure. The choice of the closed surface S F being arbitrary, we can extend R ∞ to the infinity as presented in Fig. 5. In practice, we reach to regions where the EM field is null. The only contribution left is the surface (z = 0) separating the plate and air. Thus, we obtain: The total impedance is then given by: Where Z air is the impedance of the coil in air.

Numerical validation
We adopt two test cases: firstly we aim to calculate the EM field inside the anisotropic material and the eddy current distribution in each ply. Second, we aim to detect the fiber orientation via the computation of the probe response. For that reason, two configurations were adopted as explored in the next section. For the numerical validation, we use a finite element commercial software (COMSOL Multiphysics) to compare our simulated data with FE results.

Description
The considered unidirectional CFRP sample is provided by the electric department in Newcastle University as a partner of the project. It is composed of 5 layers with the  Fig.6. Based on homogenization model described in [5], electric department in Nantes university estimate the conductivities along the principle axis in each ply to be σ t = 39000 S/m, σ l = σ n = 7.8 S/m. For this configuration, we consider a cylindrical pancake 3D Eddy Current probe with a rectangular cross-section as shown in Fig.7 and its characteristics are stored in Tab.2. It is to highlight that the adopted numerical model is generic and any shape of probe can be used.

Results
We present in Fig.8-13, the slice views of the electric field and magnetic field at 4 different depths from the top surface z = 0 to z = −4mm. We obtaine a good agreement between the FE data and simulation results. In Fig.14-18, we display the 2D distribution of eddy current density in the middle of each layer. Obviously the main orientation of the fibers can be identified. As perspective opening, one can study the evolution of J by crossing a separating surface of two layers with different fiber orientation.  In order to quantify the error, the quadratic relative error is computed. we define the quantity as in Eq(45). To illustrate the accuracy of the model, the error is calculated  at three different depths and stored in Table 3 and Table  4.  Table 3. Quadratic error on the module of electric field.
The computation time for both cases is presented in Table 5. It is obvious that the adopted semi-analytical approach is well advantageous over a classic FEM. The calculation time depends on the number of modes (M u , M v )    Table 4. Quadratic error on the module of magnetic field.    and the complexity of the geometry (will be discussed in future work).
Computation time FEM S-A model 180 mn 2 mn Table 5. Computation time comparison between the S-A model and the FEM.

Description
As presented in [6], we consider an EC probe made of 2 identical square coils. The anisotropic slab is made of one layer with fibers aligned along the Ox axis (θ = 0 • ). The anti-symmetry of the coil regarding the vertical axis leads to a difference in the Joule loss in the material during a rotated scan. Thus, we are able to identify the main orientation of the fibers thanks to the variation of the impedance.

Results
Thanks to the symmetry of the anisotropy with respect to the Oyz plane, one can apply a quarter turn scan from 0 • to 90 • rotation with 20 steps and we plot the polar diagram of the impedance. The simulated data are compared to FE data. Using Eq(45), the error on the current induced on the top surface is calculated. To complete the validation part, we end with the error on the real part and imaginary part of the impedance of the coil for three different frequencies. the results are presented in Table 8 and Table 9.

Conclusion
In the present paper, we adopt a semi-analytical approach to simulate the response of an anisotropic, planar and multilayered structure to the excitation of 3D Eddy current probe. The proposed technique can be so helpful in detecting the fiber disorientation in intermediate plies. It represents the first part of the PhD work in which we aim to compute the field induced inside the material under test and to calculate the probe response via the variation of its impedance. The semi-analytical approach is based on a modal decomposition of the field in different layers. The numerical validation consist in using Finite Element solver to compare data with numerical results. The validity and the efficiency of the formalism was demonstrated as it shows a good agreement with the FE data. This work is the first step in a set of other future work in which we aim to generalize this knowledge to deal with anisotropic material presenting a complex geometry,like a flaw or a delamination or fiber breakage.