# The fast computation of eddy current distribution and probe response in homogenized composite material based on semi-analytical approach

Houssem Chebbi, Denis Prémel

► To cite this version:

Houssem Chebbi, Denis Prémel. The fast computation of eddy current distribution and probe response in homogenized composite material based on semi-analytical approach. European Physical Journal: Applied Physics, 2020, 89 (1), pp.10901. 10.1051/epjap/2020190220. hal-02903884

HAL Id: hal-02903884

https://hal.science/hal-02903884

Submitted on 21 Jul 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The fast computation of eddy current distribution and probe response in homogenized composite material based on semi-analytical approach

Houssem Chebbi* and Denis Prémel

CEA LIST, Laboratoire de Simulation et de Modélisation en Électromagnétisme, Gif-sur-Yvette 91191, France

Received: 23 July 2019 / Received in final form: 19 December 2019 / Accepted: 4 February 2020

**Abstract.** 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.

## 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 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

<footnote>* e-mail: houssem.chebbi@cea.fr</footnote>
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.

# 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 ($\sigma_l, \sigma_t, \sigma_n$). The rotation of the fibers is referred by the angle $\theta$ introduced in the rotation matrix $R$:

$$R(\theta) = \begin{bmatrix} \cos(\theta) & -\sin(\theta) & 0 \\ \sin(\theta) & \cos(\theta) & 0 \\ 0 & 0 & 1 \end{bmatrix}. \quad (1)$$

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

$$\bar{\sigma} = R(\theta) \begin{bmatrix} \sigma_l & 0 & 0 \\ 0 & \sigma_t & 0 \\ 0 & 0 & \sigma_n \end{bmatrix} R^{-1}(\theta). \quad (2)$$

## 2.2 Physical equations

We suppose that the field components have $e^{i\omega t}$ time dependency and assuming the absence of any external electric source and external charge density, we have $J_s = 0$ and $\rho = 0$. Thus, we focus on Maxwell-Faraday law for induction and Ampere's law equations. We introduce the

<figure><img src="image_2.png" /><figcaption>Fig. 1. Structure of multilayered CFRP.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Eddy current path through contact points in transversal and cross ply directions.</figcaption></figure>

permeability and the permittivity of vacuum respectively $\mu_0$ and $\epsilon_0$:

$$\nabla \times \mathbf{E} = -i\omega\mu_0\mu_r\mathbf{H} \quad (3)$$

$$\nabla \times \mathbf{H} = i\omega\epsilon_0\bar{\epsilon}_r\mathbf{E} + \bar{\sigma}\mathbf{E}. \quad (4)$$

Therefore we define:

$$\bar{\bar{\epsilon}} = \bar{\bar{\epsilon}}_r - \frac{i\bar{\bar{\sigma}}}{\omega\epsilon_0} = \begin{bmatrix} \epsilon_{xx} & \epsilon_{xy} & 0 \\ \epsilon_{xx} & \epsilon_{xy} & 0 \\ 0 & 0 & \epsilon_{zz} \end{bmatrix}. \quad (5)$$

The permittivity tensor is denoted by $\bar{\bar{\epsilon}}_r$ refers to the uniform electric polarization of the material (in our case we choose $\bar{\bar{\epsilon}}_r = \bar{\bar{I}}$ but it is not required in general). The relative permeability is a constant $\mu_r = 1$.

# 3 Methodology

## 3.1 Modal decomposition in isotropic media

In isotropic media where $\bar{\sigma}$ is reduced to a constant ($\sigma = 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:

$$ \left\{ \begin{array}{l} [\partial_z^2 + k^2 \mu_r \varepsilon] E_x = \partial_z \partial_x E_z - ik\mu_r Z_0 \partial_y H_z \\ [\partial_z^2 + k^2 \mu_r \varepsilon] E_y = \partial_z \partial_y E_z + ik\mu_r Z_0 \partial_x H_z \end{array} \right. \qquad (6) $$

$$ \left\{ \begin{array}{l} [\partial_z^2 + k^2 \mu_r \varepsilon] H_x = \frac{ik\varepsilon}{Z_0} \partial_y E_z + \partial_z \partial_x H_z \\ [\partial_z^2 + k^2 \mu_r \varepsilon] H_y = \partial_z \partial_y H_z - \frac{ik\varepsilon}{Z_0} \partial_x E_z. \end{array} \right. \qquad (7) $$

At this stage, we suppose the variable separation $f(x, y, z) = f(x, y)e^{-i\gamma z}$, thus, one can obtain:

$$ (k_c^2 - \gamma^2) \begin{bmatrix} \mathbf{E}_t \\ \mathbf{H}_t \end{bmatrix} = \begin{bmatrix} \pm\gamma \partial_y \\ \pm\gamma \partial_x \\ -\varepsilon \partial_y \\ -\varepsilon \partial_x \end{bmatrix} E_z + \begin{bmatrix} -\mu_r \\ \mu_r \partial_x \\ \pm\gamma \partial_x \\ \pm\gamma \partial_y \end{bmatrix} H_z \quad (8) $$

where $k_c^2 = k^2 \mu_r \varepsilon_r$ and $Z_0 = \sqrt{\frac{\mu_0}{\varepsilon_0}}$. The two potentials $E_z$ and $H_z$ (referred by $\phi$) verify Helmholtz's equation to be resolved:

$$ (\partial_x + \partial_y + k_c^2 - \gamma^2)\phi = 0. \qquad (9) $$

## 3.2 Modal decomposition in anisotropic media

In the general case where $\bar{\bar{\varepsilon}}$ takes the equations (3) and (4), we obtain:

$$ \left\{ \begin{array}{l} \partial_y E_z - \partial_z E_y = -ik\mu_r Z_0 H_x \\ \partial_z E_x - \partial_x E_z = -ik\mu_r Z_0 H_y \\ \partial_x E_y - \partial_y E_x = -ik\mu_r Z_0 H_z \end{array} \right. \qquad (10) $$

$$ \left\{ \begin{array}{l} \partial_y H_z - \partial_z H_y = \frac{ik}{Z_0}(\varepsilon^{xx} E_x + \varepsilon^{xy} E_y) \\ \partial_z H_x - \partial_x H_z = \frac{ik}{Z_0}(\varepsilon^{yy} E_y + \varepsilon^{yx} E_x). \end{array} \right. \qquad (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:

$$ \partial_z \begin{bmatrix} E_x \\ E_y \end{bmatrix} = \mathcal{L}_{EG} \begin{bmatrix} H_x \\ H_y \end{bmatrix} \qquad (12) $$

$$ \partial_z \begin{bmatrix} H_x \\ H_y \end{bmatrix} = \mathcal{L}_{GE} \begin{bmatrix} E_x \\ E_y \end{bmatrix} \qquad (13) $$

where $\mathcal{L}_{EG}$ and $\mathcal{L}_{GE}$ are differential operators defined by:

$$ \mathcal{L}_{EG} = \begin{bmatrix} -\frac{Z_0 \partial_x \partial_y}{ik\varepsilon^{zz}} & \frac{Z_0 \partial_x^2}{ik\varepsilon^{zz}} - ik\mu_r \\ -\frac{Z_0 \partial_y^2}{ik\varepsilon^{zz}} + ik\mu_r & \frac{Z_0 \partial_y \partial_x}{ik\varepsilon^{zz}} \end{bmatrix} \qquad (14) $$

$$ \mathcal{L}_{GE} = \begin{bmatrix} \frac{\partial_x \partial_y}{iZ_0 k \mu_r} + ik\varepsilon^{yx} & -\frac{\partial_x^2}{iZ_0 k \mu_r} + ik\varepsilon^{yy} \\ \frac{\partial_y^2}{iZ_0 k \mu_r} - ik\varepsilon^{xx} & -\frac{\partial_x \partial_y}{iZ_0 k \mu_r} - ik\varepsilon^{xy}. \end{bmatrix} \qquad (15) $$

Applying the operator $\partial_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:

$$ -\gamma^2 [\mathbf{E}_t] = \mathcal{L}_{EG} \mathcal{L}_{GE} [\mathbf{E}_t] \qquad (16) $$

$$ -\gamma^2 [\mathbf{H}_t] = \mathcal{L}_{GE} \mathcal{L}_{EG} [\mathbf{H}_t]. \qquad (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:

$$ F(\alpha, \beta) = \sum_{u=-M_u}^{u=M_u} \sum_{v=-M_v}^{v=M_v} f(x, y) e^{-i\alpha_u x} e^{-i\beta_v y} \qquad (18) $$

where

$$ \left\{ \begin{array}{l} \alpha_u = \frac{2\pi u}{dx}, \quad -M_u < u < M_u \\ \beta_v = \frac{2\pi v}{dy}, \quad -M_v < v < M_v. \end{array} \right. \qquad (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 ($\partial_i \equiv \frac{\partial}{\partial i}$) become $\partial_x \equiv -i\alpha$ and $\partial_y \equiv -i\beta$. In the Fourier domain, the solution of Helmholtz's equation can be expressed as sum of Fourier basic functions:

$$ \tilde{E}_z = \sum_{u,v} \Gamma_{uv}^{\pm TE} F^{TE}(\alpha, \beta) e^{-i(\gamma_{uv}^{\pm} z)} \qquad (20) $$

$$ \tilde{H}_z = \sum_{u,v} \Gamma_{uv}^{\pm TM} F^{TM}(\alpha, \beta) e^{-i(\gamma_{uv}^{\pm} z)} \qquad (21) $$

with

$$ \alpha_u^2 + \beta_v^2 + \gamma_{uv}^2 = k_c^2 = k^2 \mu_r \varepsilon_r. \qquad (22) $$

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

$$ \begin{bmatrix} \tilde{E}_x \\ \tilde{E}_y \\ \tilde{H}_x \\ \tilde{H}_y \end{bmatrix} = \Psi^{\pm TE} A^{\pm TE} \phi + \Psi^{\pm TM} A^{\pm TM} \phi \qquad (23) $$

adopting the Einstein convention we obtain:

$$ \Psi^{\pm TE} = \frac{1}{k_{c_p}^2 - \gamma_{uv}^2} \begin{bmatrix} \pm \gamma_{uv} \beta_v \\ \pm \gamma_{uv} \alpha_u \\ \varepsilon \beta_v \\ -\varepsilon \alpha_u \end{bmatrix} \qquad (24) $$
$$ \Psi^{\pm TM} = \frac{1}{k_c^2 - \gamma_{uv}^2} \begin{bmatrix} -\mu_r & & \\ & \mu_r \alpha_u & \\ & \pm \gamma_{uv} \alpha_u & \\ & & \pm \gamma_{uv} \beta_v \end{bmatrix}. \quad (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^{\pm TE}$ et $A^{\pm TM}$ are determined by Boundary Conditions.

On the other hand, the eigenvalues system (Eq. (16)) provides two eigenvalues $\pm\sqrt{\lambda_1}$ et $\pm\sqrt{\lambda_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:

$$ E_t = \begin{bmatrix} E_x \\ E_y \end{bmatrix} = V_E^- a_E^- + V_E^+ a_E^+. \quad (26) $$

The eigenvalues are used to calculate the wave attenuation factors $\phi_i^{\pm}(\gamma)$. From equation (13) we reconstruct the $H_t$ components:

$$ \begin{bmatrix} H_x \\ H_y \end{bmatrix} = \frac{i}{\gamma} \mathcal{L}_{GE} \begin{bmatrix} E_x \\ E_y \end{bmatrix} = V_{H}^{-} a_{H}^{-} + V_{H}^{+} a_{H}^{+}. \quad (27) $$

Finally the modal decomposition in the anisotropic media is defined by:

$$ \begin{bmatrix} \mathbf{E}_t \\ \mathbf{H}_t \end{bmatrix} = \begin{bmatrix} \mathbf{W}^{-} & \mathbf{W}^{+} \end{bmatrix} \begin{bmatrix} \mathbf{A}^{-} \\ \mathbf{A}^{+} \end{bmatrix} \quad (28) $$

where the modes are:

$$ \mathbf{W}^{\pm} = \begin{bmatrix} \mathbf{V}_{\mathbf{E}}^{\pm} \\ \mathbf{V}_{\mathbf{H}}^{\pm} \end{bmatrix} \quad (29) $$

$$ \mathbf{A}^{-} = \begin{bmatrix} \mathbf{a}_{\mathbf{E}}^{-} \\ \mathbf{a}_{\mathbf{H}}^{-} \end{bmatrix} \quad (30) $$

$$ \mathbf{A}^{+} = \begin{bmatrix} \mathbf{a}_{\mathbf{E}}^{+} \\ \mathbf{a}_{\mathbf{H}}^{+} \end{bmatrix} \quad (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, \dots, 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 $\Psi_p^{\pm}$. At the $p$th interface, we have the outgoing waves corresponding to the coefficients $(a_p^{+up}, a_p^{-dn})$ and the incoming waves corresponding to the coefficients $(a_p^{-up}, a_p^{+dn})$. 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 0th and the Nth layers are in air, $a_0^{-}$ results from the incident field due to the coil and $a_N^{+} = 0$.

<figure><img src="image_4.png" /><figcaption>Fig. 3. General diagram: multilayered structure.</figcaption></figure>

### 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:

$$ \mathbf{E}_0 \times \mathbf{z} = \mathbf{E}_1 \times \mathbf{z} \quad (32) $$

$$ \mathbf{H}_0 \times \mathbf{z} = \mathbf{H}_1 \times \mathbf{z}. \quad (33) $$

This leads us to the equality between the two modal expansions. We consider the relationship between outputs and inputs coefficients:

$$ \begin{bmatrix} \mathbf{A}_0^+ \\ \mathbf{A}_1^- \end{bmatrix} = [\mathbf{W}_0^+ - \mathbf{W}_1^-]^{-1} [\mathbf{W}_0^- \mathbf{W}_1^+] \begin{bmatrix} \mathbf{A}_0^- \\ \mathbf{A}_1^+ \end{bmatrix} \quad (34) $$

In the case of semi-infinite case $A_1^+ = 0$. $A_0^-$ is computed using the relationship:

$$ A_0^- = \left[ \frac{(\varepsilon\beta - \gamma\alpha)}{k_c^2 - \gamma^2} \right]^{-1} \begin{bmatrix} E_z^{\text{inc}} \\ H_z^{\text{inc}} \end{bmatrix} = [\Gamma_0^-]^{-1} \begin{bmatrix} E_z^{\text{inc}} \\ H_z^{\text{inc}} \end{bmatrix} \quad (35) $$

where $\Gamma_0^-$ is the part of the mode $\Psi_0^-$ that corresponds to the magnetic field and the incident EM field provided by the coil in air $E_z^{\text{inc}}, H_z^{\text{inc}}$ 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.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Illustration of the field propagation in the pth layer.</figcaption></figure>

Globally the sought global $S_g$ matrix is such as:

$$ \begin{bmatrix} a_0^+ \\ a_N^- \end{bmatrix} = S_g \begin{bmatrix} a_N^+ \\ a_0^- \end{bmatrix} \qquad (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:

$$ \begin{bmatrix} a_p^{+up} \\ a_p^{-dn} \end{bmatrix} = \begin{bmatrix} S_p^{11} & S_p^{12} \\ S_p^{21} & S_p^{22} \end{bmatrix} \begin{bmatrix} a_p^{+dn} \\ a_p^{-up} \end{bmatrix} = S_p \begin{bmatrix} a_p^{+dn} \\ a_p^{-up} \end{bmatrix} \quad (37) $$

where

$$ S_p = [\Psi_{p-1}^{+} - \Psi_p^{-}]^{-1} [\Psi_p^{+} - \Psi_{p-1}^{-}] \quad (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:

$$ \begin{bmatrix} a_p^{+dn} \\ a_p^{-up} \end{bmatrix} = \begin{bmatrix} \phi_p^+ & 0 \\ 0 & \phi_p^- \end{bmatrix} \begin{bmatrix} a_{p+1}^{+up} \\ a_p^{-dn} \end{bmatrix} \quad (39) $$

where $\phi^\pm = e^{\pm ik_{cp}} \gamma e_p$ is the th attenuation factor through the pth layer and $e_p$ is the thickness inside the material.

## 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 to 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 $\delta$ reasonable with respect to the low values of $\sigma_{zz}$, we work in the frequency range 0.8–10 MHz. Although, $\delta$ depends

<figure><img src="image_6.png" /><figcaption>Fig. 5. Closed surface for impedance calculation.</figcaption></figure>

also on the number and the sequence of plies. Thus, at fixed depth 'p' we consider:

$$ |\mathbf{J}_{\mathbf{p}}| = \sqrt{\mathbf{J}_{px}^2 + \mathbf{J}_{py}^2 + \mathbf{J}_{pz}^2} \quad (40) $$

where

$$ \left\{ \begin{array}{l} \mathbf{J}_{px} = \sigma_{xx}^p \mathbf{E}_x^p + \sigma_{xy}^p \mathbf{E}_y^p \\ \mathbf{J}_{py} = \sigma_{yx}^p \mathbf{E}_x^p + \sigma_{yy}^p \mathbf{E}_y^p \\ \mathbf{J}_{pz} = \sigma_{zz}^p \mathbf{E}_z^p \end{array} \right. \quad (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 $\Delta 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 $\Delta Z$ of the coil impedance is given by the following formula:

$$ \Delta Z = \frac{1}{I_0^2} \iint_{S_F} (\mathbf{E}^i \times \mathbf{H} - \mathbf{E} \times \mathbf{H}^i) \cdot \mathbf{n} \, dx \, dy \quad (42) $$

where $I_0$ is the driven current intensity circulating in the coil. $\mathbf{E}^i$ and $\mathbf{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 $\mathbf{E}$ and $\mathbf{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_\infty$ 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:

$$ \Delta Z = \frac{1}{I_0^2} \iint_{S_{z=0}} (\mathbf{E}^i \times \mathbf{H} - \mathbf{E} \times \mathbf{H}^i) \cdot \mathbf{z} \, dx \, dy. \quad (43) $$

The total impedance is then given by:

$$ Z = Z_{air} + \Delta Z \quad (44) $$

where $Z_{air}$ is the impedance of the coil in air.
<table><caption>**Table 1.** Numerical parameters of the model.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Number of modes along O<sub>x</sub></td><td>M<sub>u</sub></td><td>35</td></tr><tr><td>Number of modes along O<sub>y</sub></td><td>M<sub>v</sub></td><td>35</td></tr><tr><td>Spatial bound along x</td><td>x<sub>max</sub></td><td>32 mm</td></tr><tr><td>Spatial bound along y</td><td>y<sub>max</sub></td><td>32 mm</td></tr></tbody></table>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Config1: flawless unidirectional CFRP sample.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** cylindrical pancake coil with a rectangular cross-section.</figcaption></figure>

# 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

<table><caption>**Table 2.** Parameters of the first configuration.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Specimen</td><td></td><td></td></tr><tr><td>longitudinal Conductivity</td><td>σ<sub>l</sub></td><td>39<sup>-3</sup> MS/m</td></tr><tr><td>transversal Conductivity</td><td>σ<sub>t</sub></td><td>7.9<sup>-6</sup> MS/m</td></tr><tr><td>normal Conductivity</td><td>σ<sub>n</sub></td><td>7.9<sup>-6</sup> MS/m</td></tr><tr><td>Total thickness</td><td>e</td><td>1 mm</td></tr><tr><td>Probe</td><td></td><td></td></tr><tr><td>Frequency</td><td>f</td><td>1 MHz</td></tr><tr><td>Excitation current</td><td>I<sub>0</sub></td><td>1 A</td></tr><tr><td>Internal radius</td><td>r<sub>int</sub></td><td>4 mm</td></tr><tr><td>External radius</td><td>r<sub>ext</sub></td><td>6 mm</td></tr><tr><td>Height</td><td>H</td><td>1 mm</td></tr><tr><td>Number of turns</td><td>N</td><td>50</td></tr><tr><td>Liftoff</td><td>l<sub>0</sub></td><td>0.8 mm</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Slice view of $\Re(E_x)$ and $\Im(E_x)$ at several depths.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Slice view of $\Re(E_y)$ and $\Im(E_y)$ at several depths.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 10. Slice view of $\Re(E_z)$ and $\Im(E_z)$ at several depths.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Slice view of $\Re(H_x)$ and $\Im(H_x)$ at several depths.</figcaption></figure>

conductivities along the principle axis in each ply to be $\sigma_t = 39000$ S/m, $\sigma_l = \sigma_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.

### 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

<figure><img src="image_14.png" /><figcaption>Fig. 12. Slice view of $\Re(H_y)$ and $\Im(H_y)$ at several depths.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13. Slice view of $\Re(H_z)$ and $\Im(H_z)$ at several depths.</figcaption></figure>

<table><caption>Table 3. Quadratic error on the module of electric field.</caption><thead><tr><th>Depth (mm)</th><th>ξ (%) E<sub>x</sub></th><th>ξ (%) E<sub>y</sub></th><th>ξ (%) E<sub>z</sub></th></tr></thead><tbody><tr><td>0</td><td>1.2</td><td>1.04</td><td>4.2</td></tr><tr><td>0.4</td><td>2.34</td><td>2.03</td><td>6.76</td></tr><tr><td>0.8</td><td>2.5</td><td>2.36</td><td>5.8</td></tr></tbody></table>

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.

$$ \xi = \sqrt{\frac{\| F_{FE} - F_{model} \|_{L_2}}{\| F_{FE} \|_{L_2}}}. \quad (45) $$
<figure><img src="image_16.png" /><figcaption>Fig. 14. Top view of **J** at the surface *z* = -0.5 mm.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 15. Top view of **J** at the surface *z* = -1.5 mm.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Top view of **J** at the surface *z* = -2.5 mm.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 17. Top view of **J** at the surface *z* = -3.5 mm.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 18. Top view of **J** at the surface *z* = -4.5 mm.</figcaption></figure>

<table><caption>Table 4. Quadratic error on the module of magnetic field.</caption><thead><tr><th>Depth (mm)</th><th>ξ (%) H<sub>x</sub></th><th>ξ (%) H<sub>y</sub></th><th>ξ (%) H<sub>z</sub></th></tr></thead><tbody><tr><td>0</td><td>0.88</td><td>0.75</td><td>0.79</td></tr><tr><td>0.4</td><td>1.0</td><td>0.987</td><td>1.1</td></tr><tr><td>0.8</td><td>1.5</td><td>1.45</td><td>1.44</td></tr></tbody></table>

<table><caption>Table 5. Computation time comparison between the S-A model and the FEM.</caption><thead><tr><th>Computation time</th><th>FEM</th><th>S-A model</th></tr></thead><tbody><tr><td></td><td>≈180 mn</td><td>≈2 mn</td></tr></tbody></table>

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*<sub>u</sub>, *M*<sub>v</sub>) and the complexity of the geometry (will be discussed in future work).

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

### 4.2.2 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 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
<figure><img src="image_21.png" /><figcaption>**Fig. 19.** Config2: a planar anisotropic slab excited by a rectangular 3D EC probe.</figcaption></figure>

<table><caption>**Table 6.** Parameters of the second configuration.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Specimen</td><td></td><td></td></tr><tr><td>longitudinal Conductivity</td><td>σ<sub>l</sub></td><td>1 MS/m</td></tr><tr><td>transverse Conductivity</td><td>σ<sub>t</sub></td><td>0.01 MS/m</td></tr><tr><td>normal Conductivity</td><td>σ<sub>n</sub></td><td>0.01 MS/m</td></tr><tr><td>Total thickness</td><td>E</td><td>1 mm</td></tr><tr><td>orientation</td><td>θ</td><td>0°</td></tr><tr><td>Probe</td><td></td><td></td></tr><tr><td>Frequency</td><td>f</td><td>0.8 MHz</td></tr><tr><td>Excitation current</td><td>I<sub>0</sub></td><td>1 A</td></tr><tr><td>Length</td><td>L</td><td>10 mm</td></tr><tr><td>Width</td><td>W</td><td>5 mm</td></tr><tr><td>Height</td><td>e</td><td>1 mm</td></tr><tr><td>Number of turns</td><td>N</td><td>3</td></tr><tr><td>Liftoff</td><td>l<sub>0</sub></td><td>1 mm</td></tr><tr><td>Distance</td><td>d</td><td>10 mm</td></tr></tbody></table>

<table><caption>**Table 7.** Quadratic error on the EC density on the top surface.</caption><thead><tr><th>Depth (mm)</th><th>ξ (%) || J ||<sub>L<sub>2</sub></sub></th></tr></thead><tbody><tr><td>0</td><td>1.76</td></tr><tr><td>0.4</td><td>1.81</td></tr><tr><td>0.8</td><td>1.87</td></tr></tbody></table>

<table><caption>**Table 8.** Real part of the total impedance at different frequencies with relative quadratic error.</caption><thead><tr><th>Frequency (kHz)</th><th>Ξ(Z<sub>FEM</sub>)</th><th>Ξ(Z)</th><th>ξ (%)</th><th>Ξ(Z)</th></tr></thead><tbody><tr><td>800</td><td>0.00219</td><td>0.0022</td><td>0.4545</td><td></td></tr><tr><td>10<sup>3</sup></td><td>0.00342</td><td>0.0034</td><td>0.58</td><td></td></tr><tr><td>10<sup>4</sup></td><td>0.2235</td><td>0.2227</td><td>0.35</td><td></td></tr></tbody></table>

<table><caption>**Table 9.** Imaginary part of the total impedance at different frequencies with relative quadratic error.</caption><thead><tr><th>Frequency (kHz)</th><th>Ξ(Z<sub>FEM</sub>)</th><th>Ξ(Z)</th><th>ξ (%)</th><th>Ξ(Z)</th></tr></thead><tbody><tr><td>800</td><td>0.5608</td><td>0.562</td><td>0.213</td><td></td></tr><tr><td>10<sup>3</sup></td><td>0.533</td><td>0.529</td><td>0.75</td><td></td></tr><tr><td>10<sup>4</sup></td><td>5.2012</td><td>5.1687</td><td>0.624</td><td></td></tr></tbody></table>

<figure><img src="image_22.png" /><figcaption>**Fig. 20.** Polar diagram of ΔZ.</figcaption></figure>

part of the impedance of the coil for three different frequencies. The results are presented in Tables 8 and 9.

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

This work is financially supported by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 722134. The research has been undertaken as a part of NDToN Air 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

1. S. Gholizadeh, Proc. Struct. Integr. **1**, 50 (2016)

2. X. Li, W. Yin, Z. Liu, P.J. Withers, A.J. Peyton, in *17th World Conference on Nondestructive Testing*, 2008

3. K. Kiyoshi, H. Hiroshi, K. Gouki, J. Press. Vessel Technol. **135**, 1 (2013)

4. B. Salski, W. Gwarek, P. Korpas, IEEE Trans. Microw. Theory Tech. **62**, 1535 (2013)

5. G. Wasselynck, D. Trichet, J. Fouladgar, IEEE Trans. Magn. **49**, 1825 (2013)

6. W. Yin, P.J. Withers, U. Sharma, A.J. Peyton, IEEE Trans. Instrum. Measur. **58**, 738 (2009)

7. J. Cheng, J. Qiu, H. Ji, E. Wang, T. Takagi, T. Uchimoto, Composites Part B **110**, 141 (2017)

8. T.M. Roberts, H.A. Sabbagh, L.D. Sabbagh, Math. Phys. **29**, 2675 (1988)

9. T.M. Roberts, IEEE Trans. Magn. **26**, 3064 (1990)

10. I. Bardi, R. Remski, D. Perry, Z. Cendes, IEEE Trans. Magn. **38**, 641 (2002)

11. F. Caire, D. Prémel, G. Granet, Eur. Phys. J. Appl. Phys. **64**, 24511 (2013)

12. B.A. Auld, J.C. Moulder, J. Nondestruct. Eval. **18**, 3 (1999)

13. J. Chandezon, G. Raoult, D. Maystre, J. Opt. **11**, 235 (1980)

14. N.P.K. Cotter, T.W. Preist, J.R. Sambles, J. Opt. Soc. Am. A **12**, 1097 (1995)

15. T.P. Theodoulidis, IEEE Trans. Magn. **41**, 2447 (2005)

16. G. Granet, J.P. Plumey, J. Chandezon, Pure Appl. Opt. **4**, 1 (1995)

17. L. Li, J. Opt. Soc. Am. A **13**, 1024 (1996)

18. H.M. Wen, Compos. Sci. Technol. **61**, 1163 (2001)

19. J. Cheng, H. Ji, J. Qiu, T. Takagi, T. Uchimoto, N. Hu, NDT&E Int. **68**, 1 (2014)

**Open Access** 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.

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)