© H. Chebbi and D. Prémel, EDP Sciences, 2020

1 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 multi-frequency analysis and the identification of fiber orientation [2]. Therefore, many research works show interest in developing inspection procedure using new specific sensors and signal processing [3,4]. Moreover, CFRP is a non-homogeneous 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 [8,9] based on integral equations and Green’s Dyads. In addition, analytical models were presented in [10] 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 [11] 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 demonstratethe effectiveness of the method.

2 Modeling framework

2.1 Material characteristics

As described in Figure 1, a CFRP structure is made 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 Figure 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³ and 5 × 10⁴ S/m, and between 1 and 5 × 10² 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: (1)

Thus, the conductivity tensor in each ply is obtained as follow: (2)

Fig. 1 Structure of multilayered CFRP.

Fig. 2 Eddy current path through contact points in transversal and cross ply directions.

2.2 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 μ₀ and ε₀: (3) (4)

Therefore we define: (5)

The permittivity tensor is denoted by refers to the uniform electric polarization of the material (in our case we choose but it is not required in general). The relative permeability is a constant μ_r = 1.

3 Methodology

3.1 Modal decomposition in isotropic media

In isotropic media where is reduced to a constant (σ = 0 in air), so that from equation (3) and equation (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: (6) (7)

At this stage, we suppose the variable separation f(x, y, z) = f(x, y)e^−iγz, thus, one can obtain: (8)

where and . The two potentials E_z and H_z (referred by ϕ) verify Helmholtz’s equation to be resolve: (9)

3.2 Modal decomposition in anisotropic media

In the general case where takes the equations (3) and (4), we obtain: (10) (11)

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: (12) (13)

where and is differential operators defined by: (14) (15)

Applying the operator ∂_z on equations (12) and (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: (16) (17)

Finally we can recompute the E_z and H_z from the last two equalities of equations (10) and (11).

3.3 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: (18)

where (19)

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 L = (2M_u + 1)(2M_v + 1). Besides, the partial derivatives () become ∂_x ≡−iα and ∂_y ≡−iβ. In the Fourier domain, the solution of Helmholtz’s equation can be expressed as sum of Fourier basic functions: (20) (21)

with (22)

These expansions are known as Rayleigh expansions. From equation (8) we obtain the modal expansion: (23)

adopting the Einstein convention we obtain: (24) (25)

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^±TE et A^±TM are determined by Boundary Conditions.

On the other hand, the eigenvalues system (Eq. (16)) provides two eigenvalues et corresponding to the direction of the propagation of the wave, and two associated eigenvectors v₁ and v₂. So E_x and E_y are a linear combination of eigenvectors. In a vector writing: (26)

The eigenvalues are used to calculate the wave attenuation factors . From equation (13) we reconstruct the H_t components: (27)

Finally the modal decomposition in the anisotropic media is defined by: (28)

where the modes are: (29) (30) (31)

3.4 General configuration

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

Fig. 3 General diagram: multilayered structure.

3.4.1 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: (34)

In the case of semi-infinite case . is computed using the relationship: (35)

where is the part of the mode that corresponds to the magnetic field and the incident EM field provided by the coil in air is calculated by means of analytical model already implemented into CIVA software [15].

3.4.2 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 Figure 3 a (4 × 4) matrix connects the incoming waves with the outgoing waves. Globally the sought global S_g matrix is such as: (36)

The global S_g matrix is obtained through classical recursion formulas between intermediate S_p matrices at each separating interface. As displayed in Figure 4, the intermediate scattering matrix S_j takes into account the modal decomposition of the field on both sides of the pth interface and the BC to calculate the unknown intermediate coefficients. Thus we have the output–input relationship: (37)

where (38)

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: (39)

where is th attenuation factor through the pth layer and e_p is the thickness inside the material.

Fig. 4 Illustration of the field propagation in the pth layer.

3.5 Eddy current density and skin effect

In metallic structure, Eddy currents are closed loops of induced current circulating in planes perpendicular to the 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: (40)

where (41)

3.6 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 reciprocity 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₀ is the driven current intensity circulating in the coil. Eⁱ and Hⁱ 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 Figure 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: (43)

The total impedance is then given by: (44)

where Z_air is the impedance of the coil in air.

Fig. 5 Closed surface for impedance calculation.

4 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, numerical parameters are in Table 1 and we use a finite element commercial software (COMSOL Multiphysics) to compare our simulated data with FE results.

4.1 First configuration

4.1.1 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 sequence [0°, −45°, 90°, 45°, 0°] as shown in Figure 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 = 39 000 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 Figure 7 and its characteristics are stored in Table 2. It is to highlight that the adopted numerical model is generic and any shape of probe can be used.

Fig. 6 Config1: flawless unidirectional CFRP sample.

Fig. 7 cylindrical pancake coil with a rectangular cross-section.

4.1.2 Results

We present in Figures 8–13, the slice views of the electric field and magnetic field at 4 different depths from the top surface z = 0 to z = −4 mm. We obtained a good agreement between the FE data and simulation results. In Figures 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 equation (45). To illustrate the accuracy of the model, the error is calculated at three different depths and stored in Tables 3 and 4. (45)

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) and the complexity of the geometry (will be discussed in future work).

Fig. 8 Slice view of ℜ(Ex) and ℑ(Ex) at several depths.

Fig. 9 Slice view of ℜ(Ey) and ℑ(Ey) at several depths.

Fig. 10 Slice view of ℜ(Ez) and ℑ(Ez) at several depths.

Fig. 11 Slice view of ℜ(Hx) and ℑ(Hx) at several depths.

Fig. 12 Slice view of ℜ(Hy) and ℑ(Hy) at several depths.

Fig. 13 Slice view of ℜ(Hz) and ℑ(Hz) at several depths.

Fig. 14 Top view of J at the surface z = −0.5 mm.

Fig. 15 Top view of J at the surface z = −1.5 mm.

Fig. 16 Top view of J at the surface z = −2.5 mm.

Fig. 17 Top view of J at the surface z = −3.5 mm.

Fig. 18 Top view of J at the surface z = −4.5 mm.

4.2 Second configuration

4.2.1 Description

As presented in Figure 19 (see [6]), we consider an EC probe made of 2 identical square coils, the parameters of the configuration are stored in Table 6. 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.

Fig. 19 Config2: a planar anisotropic slab excited by a rectangular 3D EC probe.

4.2.2 Results

Thanks to the symmetry of the anisotropy with respect to the Oyz plane, one can apply a quarter turnscan from 0° to 90° rotation with 20 steps and we plot the polar diagram of the impedance as shown in Figure 20. The simulated data are compared to FE data. Using equation (45), the error on the current induced on the top surface is calculated and stored in Table 7. 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 Tables 8 and 9.

Fig. 20 Polar diagram of ΔZ.

5 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.

Acknowledgements

This work is financially supported by the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No 722134. The research has been undertaken as a part of NDTonAir project.

Author contribution statement

All authors have contributed equally to the paper. H.C. and D.P. are responsible for analytic and numerical calculations.

References

All Tables

All Figures

	Fig. 1 Structure of multilayered CFRP.
In the text

	Fig. 2 Eddy current path through contact points in transversal and cross ply directions.
In the text

	Fig. 3 General diagram: multilayered structure.
In the text

	Fig. 4 Illustration of the field propagation in the pth layer.
In the text

	Fig. 5 Closed surface for impedance calculation.
In the text

	Fig. 6 Config1: flawless unidirectional CFRP sample.
In the text

	Fig. 7 cylindrical pancake coil with a rectangular cross-section.
In the text

	Fig. 8 Slice view of ℜ(Ex) and ℑ(Ex) at several depths.
In the text

	Fig. 9 Slice view of ℜ(Ey) and ℑ(Ey) at several depths.
In the text

	Fig. 10 Slice view of ℜ(Ez) and ℑ(Ez) at several depths.
In the text

	Fig. 11 Slice view of ℜ(Hx) and ℑ(Hx) at several depths.
In the text

	Fig. 12 Slice view of ℜ(Hy) and ℑ(Hy) at several depths.
In the text

	Fig. 13 Slice view of ℜ(Hz) and ℑ(Hz) at several depths.
In the text

	Fig. 14 Top view of J at the surface z = −0.5 mm.
In the text

	Fig. 15 Top view of J at the surface z = −1.5 mm.
In the text

	Fig. 16 Top view of J at the surface z = −2.5 mm.
In the text

	Fig. 17 Top view of J at the surface z = −3.5 mm.
In the text

	Fig. 18 Top view of J at the surface z = −4.5 mm.
In the text

	Fig. 19 Config2: a planar anisotropic slab excited by a rectangular 3D EC probe.
In the text

	Fig. 20 Polar diagram of ΔZ.
In the text