Issue 
Eur. Phys. J. Appl. Phys.
Volume 89, Number 1, January 2020



Article Number  10901  
Number of page(s)  11  
Section  Physics of Energy Transfer, Conversion and Storage  
DOI  https://doi.org/10.1051/epjap/2020190220  
Published online  10 April 2020 
https://doi.org/10.1051/epjap/2020190220
Regular Article
The fast computation of eddy current distribution and probe response in homogenized composite material based on semianalytical approach
CEA LIST, Laboratoire de Simulation et de Modélisation en Électromagnétisme,
GifsurYvette 91191, France
^{*} email: houssem.chebbi@cea.fr
Received:
23
July
2019
Received in final form:
19
December
2019
Accepted:
4
February
2020
Published online: 10 April 2020
Due to the excessive use of composites in the industrial field, 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 fiber disorientation. A semianalytical model based on a modal approach is developed for the fast computation of quasistatic field 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 field is expressed as a sum of eigenmodes. 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 fibers in the inspected zone. For numerical validation, simulated data provided by the model are compared to finite element data.
© H. Chebbi and D. Prémel, EDP Sciences, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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, outofplan 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 contactfree 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 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 [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 SemiAnalytical (SA) 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 SA model in [11] where the second order potential formalism where used, our SA 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 eigenmodes. 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 Smatrix 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 nonhomogeneity 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 SA 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 Smatrix 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 nonconductive 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 crossply 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 crossply 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 MaxwellFaraday law for induction and Ampere’s law equations. We introduce the permeability and the permittivity of vacuum respectively μ_{0} and ε_{0}: (3) (4)
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)
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)
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_{1} and v_{2}. 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 “airplate” 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 semiinfinite 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 Smatrix algorithm
The Smatrix 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)
In a recursive way, all the intermediate Smatrices are concatenated taking into account the relationship between the coefficient describing the wave propagation through the layers.
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)
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_{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 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.
Numerical parameters of the model.
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 crosssection 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 crosssection. 
Parameters of the first configuration.
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 semianalytical 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 ℜ(E_{x}) and ℑ(E_{x}) at several depths. 
Fig. 9 Slice view of ℜ(E_{y}) and ℑ(E_{y}) at several depths. 
Fig. 10 Slice view of ℜ(E_{z}) and ℑ(E_{z}) at several depths. 
Fig. 11 Slice view of ℜ(H_{x}) and ℑ(H_{x}) at several depths. 
Fig. 12 Slice view of ℜ(H_{y}) and ℑ(H_{y}) at several depths. 
Fig. 13 Slice view of ℜ(H_{z}) and ℑ(H_{z}) 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. 
Quadratic error on the module of electric field.
Quadratic error on the module of magnetic field.
Computation time comparison between the SA model and the FEM.
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 antisymmetry 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. 
Parameters of the second configuration.
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. 
Quadratic error on the EC density on the top surface.
Real part of the total impedance at different frequencies with relative quadratic error.
Imaginary part of the total impedance at different frequencies with relative quadratic error.
5 Conclusion
In the present paper, we adopt a semianalytical 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 semianalytical 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 SklodowskaCurie 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
 S. Gholizadeh, Proc. Struct. Integr. 1, 50 (2016) [Google Scholar]
 X. Li, W. Yin, Z. Liu, P.J. Withers, A.J. Peyton, in 17th World Conference on Nondestructive Testing, 2008 [Google Scholar]
 K. Kiyoshi, H. Hiroshi, K. Gouki, J. Press. Vessel Technol. 135, 1 (2013) [Google Scholar]
 B. Salski, W. Gwarek, P. Korpas, IEEE Trans. Microw. Theory Tech. 62, 1535 (2013) [Google Scholar]
 G. Wasselynck, D. Trichet, J. Fouladgar, IEEE Trans. Magn. 49, 1825 (2013) [Google Scholar]
 W. Yin, P.J. Withers, U. Sharma, A.J. Peyton, IEEE Trans. Instrum. Measur. 58, 738 (2009) [Google Scholar]
 J. Cheng, J. Qiu, H. Ji, E. Wang, T. Takagi, T. Uchimoto, Composites Part B 110, 141 (2017) [Google Scholar]
 T.M. Roberts, H.A. Sabbagh, L.D. Sabbagh, Math. Phys. 29, 2675 (1988) [Google Scholar]
 T.M. Roberts, IEEE Trans. Magn. 26, 3064 (1990) [Google Scholar]
 I. Bardi, R. Remski, D. Perry, Z. Cendes, IEEE Trans. Magn. 38, 641 (2002) [Google Scholar]
 F. Caire, D. Prémel, G. Granet, Eur. Phys. J. Appl. Phys. 64, 24511 (2013) [CrossRef] [EDP Sciences] [Google Scholar]
 B.A. Auld, J.C. Moulder, J. Nondestruct. Eval. 18, 3 (1999) [Google Scholar]
 J. Chandezon, G. Raoult, D. Maystre, J. Opt. 11, 235 (1980) [Google Scholar]
 N.P.K.Cotter, T.W. Preist, J.R. Sambles, J. Opt. Soc. Am. A 12, 1097 (1995) [Google Scholar]
 T.P. Theodoulidis, IEEE Trans. Magn. 41, 2447 (2005) [Google Scholar]
 G. Granet, J.P. Plumey, J. Chandezon, Pure Appl. Opt. 4, 1 (1995) [Google Scholar]
 L. Li, J. Opt. Soc. Am. A 13, 1024 (1996) [Google Scholar]
 H.M. Wen, Compos. Sci. Technol. 61, 1163 (2001) [Google Scholar]
 J. Cheng, H. Ji, J. Qiu, T. Takagi, T. Uchimoto, N. Hu, NDT&E Int. 68, 1 (2014) [Google Scholar]
Cite this article as: Houssem Chebbi, Denis Prémel, The Fast Computation of Eddy Current Distribution and Probe Response in Homogenized Composite Material Based on SemiAnalytical Approach, Eur. Phys. J. Appl. Phys. 89, 10901 (2020)
All Tables
Real part of the total impedance at different frequencies with relative quadratic error.
Imaginary part of the total impedance at different frequencies with relative quadratic error.
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 crosssection. 

In the text 
Fig. 8 Slice view of ℜ(E_{x}) and ℑ(E_{x}) at several depths. 

In the text 
Fig. 9 Slice view of ℜ(E_{y}) and ℑ(E_{y}) at several depths. 

In the text 
Fig. 10 Slice view of ℜ(E_{z}) and ℑ(E_{z}) at several depths. 

In the text 
Fig. 11 Slice view of ℜ(H_{x}) and ℑ(H_{x}) at several depths. 

In the text 
Fig. 12 Slice view of ℜ(H_{y}) and ℑ(H_{y}) at several depths. 

In the text 
Fig. 13 Slice view of ℜ(H_{z}) and ℑ(H_{z}) 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.