Eur. Phys. J. Appl. Phys.
Volume 89, Number 1, January 2020
|Number of page(s)||11|
|Section||Physics of Energy Transfer, Conversion and Storage|
|Published online||10 April 2020|
The fast computation of eddy current distribution and probe response in homogenized composite material based on semi-analytical approach
CEA LIST, Laboratoire de Simulation et de Modélisation en Électromagnétisme,
Gif-sur-Yvette 91191, France
* e-mail: firstname.lastname@example.org
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 semi-analytical model based on a modal approach is developed for the fast computation of quasi-static 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 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 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.
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 , 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 . 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  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  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  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  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.
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 × 103 and 5 × 104 S/m, and between 1 and 5 × 102 S/m in the transverse and cross-ply direction . 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)
Structure of multilayered CFRP.
Eddy current path through contact points in transversal and cross ply directions.
We suppose that the field components have eiωt time dependency and assuming the absence of any external electric source and external charge density, we have Js = 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: (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.
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 Ex, Ey and magnetic field denoted by Hx, Hy can be expressed in function of two potential Ez and Hz so called the TE/TM decomposition: (6) (7)
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)
The Mu and Mv 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 = (2Mu + 1)(2Mv + 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)
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 v1 and v2. So Ex and Ey 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 Ht components: (27)
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 ep. 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 .
General diagram: multilayered structure.
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 (Js = 0), the BC at a separating interface are defined by the continuity of the tangential components of the EM field:
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 .
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 Sg matrix is such as: (36)
The global Sg matrix is obtained through classical recursion formulas between intermediate Sp matrices at each separating interface. As displayed in Figure 4, the intermediate scattering matrix Sj 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 S-matrices 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 ep is the thickness inside the material.
Illustration of the field propagation in the pth layer.
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)
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 . It has been proved that the variation ΔZ of the coil impedance is given by the following formula: (42)
where I0 is the driven current intensity circulating in the coil. Ei and Hi 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) , 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 SF 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)
where Zair is the impedance of the coil in air.
Closed surface for impedance calculation.
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.
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 , 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.
Config1: flawless unidirectional CFRP sample.
cylindrical pancake coil with a rectangular cross-section.
Parameters of the first configuration.
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 (Mu, Mv) and the complexity of the geometry (will be discussed in future work).
Slice view of ℜ(Ex) and ℑ(Ex) at several depths.
Slice view of ℜ(Ey) and ℑ(Ey) at several depths.
Slice view of ℜ(Ez) and ℑ(Ez) at several depths.
Slice view of ℜ(Hx) and ℑ(Hx) at several depths.
Slice view of ℜ(Hy) and ℑ(Hy) at several depths.
Slice view of ℜ(Hz) and ℑ(Hz) at several depths.
Top view of J at the surface z = −0.5 mm.
Top view of J at the surface z = −1.5 mm.
Top view of J at the surface z = −2.5 mm.
Top view of J at the surface z = −3.5 mm.
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 S-A model and the FEM.
As presented in Figure 19 (see ), 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.
Config2: a planar anisotropic slab excited by a rectangular 3D EC probe.
Parameters of the second configuration.
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.
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.
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.
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.
All authors have contributed equally to the paper. H.C. and D.P. are responsible for analytic and numerical calculations.
- S. Gholizadeh, Proc. Struct. Integr. 1, 50 (2016) [CrossRef] [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) [CrossRef] [Google Scholar]
- J. Cheng, J. Qiu, H. Ji, E. Wang, T. Takagi, T. Uchimoto, Composites Part B 110, 141 (2017) [CrossRef] [Google Scholar]
- T.M. Roberts, H.A. Sabbagh, L.D. Sabbagh, Math. Phys. 29, 2675 (1988) [CrossRef] [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) [CrossRef] [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) [CrossRef] [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) [CrossRef] [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 Semi-Analytical Approach, Eur. Phys. J. Appl. Phys. 89, 10901 (2020)
Eddy current path through contact points in transversal and cross ply directions.
|In the text|
Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.