Composite Material Characterization using Eddy Current by 3 D FEM Associated with Iterative Technique

In this paper, an iterative technique, employing the ⃗ formulation associated with the finite element method, based on Maxwell equations and the Biot and Savart law, is used for analyzing the density of eddy currents in composite carbon fiber reinforced polymer (CFRP) materials. For this purpose, a code has been developed for solving an electromagnetic 3D non-destructive evaluation problem. This latter permits the characterization of this CFRP and determinate of fibers orientation using the impedance variation which is implanted in polar diagram. Firstly, the obtained results are compared with those of the analytical model. This comparison reveals a high concordance which proves the validity of the proposed method. Secondly, three different applications are shown for illustrating the characterization of unidirectional, bidirectional and multidirectional piece using a rectangular coil plotted in normalized impedance variation diagram.


Introduction
Nowadays, carbon fiber composite materials (CFRP) occupy a privileged place in industrial fields, [1] and are more and more popular compared to isotropic materials, because of their mechanical properties, such as high mechanical tensile strength, light weight, corrosion resistance and ability to adapt to complex geometric shapes.This is what makes their use widespread in several industrial sectors such as: aerospace, air transport, maritime and rail, health, sports, wind turbines, aircraft fuselage and many others, [2][3][4][5][6].
Nevertheless, contrary to the above advantages of CFRP materials, these can be altered in many ways during manufacturing, assembly and life stages by a variety of different damages and / or structural failures.As a result, the wide applicability of these materials has created the need for the development of fast and reliable systems for non-destructive testing (NDT), both for quality control during a manufacturing process as well as for monitoring the state of their structures in use [1].
The electrical conductivity of carbon fiber composite materials depends on the volume fraction, the structure and the orientation of the fibers in the material.In laminated carbon fibers composite materials Figure 2 (a); each layer contains unidirectional fibers embedded in an insulating matrix.
The electrical conductivity of the CFRP in the transverse direction to the fibers is not zero, [7][8][9] because there are contacts between the fibers and this is due to the fact that the fibers are not arranged in a perfectly straight way and also some fibers are not completely covered by the insulating matrix [10,11].This conductivity varies between 10 and 100 S/m in the transverse direction to the fibers and between 5.10 3 and 5.10 4 S/m in the longitudinal direction of the fibers [12].For the same reasons, there is conductivity according to the thickness of a composite laminated material [10].
Several scientific works, in the literature, use the ⃗iterative method for modeling and simulating of problems of Eddy Current non Destructive Testing (EC-NDT) [13][14][15][16].In 2009, [15] exploited this iterative method associated with the finite difference method to characterize the composite carbon fiber materials (CFRP), and in 2012, [16] used the same method coupled with the finite volume method to characterize isotropic materials and CFRP composites possibly containing defects.
In this context, we have developed a Matlab calculation code based on ⃗ -iterative technique solved by 3D finite element method, to solve an EC-NDT problem that contains an absolute circular sensor and unidirectional plate.The obtained results are compared with the analytical model ones in the aim to validate the proposed method [17].Hence, a second rectangular shape sensor is used with a composite unidirectional and multi directional plate for characterizing the fibers orientations, based on the impedance variation.The mesh of the sensor / plate is made separately using Gmsh software [18].

Modelisation of a composite material in eddy current NDT domain
The problem of eddy currents NDT to be studied consists of an inductor (source) where the excitation ( ) is imposed and a composite material region representing the plate to be controlled in absence of air region which makes it possible to minimize the mesh.The system is shown in Figure 1.The homogenized conductivity tensor of ply where the fibers are aligned with the x-axis of reference system is expressed as follows: Where σ L , σ T and σ Z respectively represent the electrical conductivities of the composite material (CFRP), such that: σ L : is the conductivity along the longitudinal direction.σ T : is the conductivity according to the cross section of the fibers.σ Z : is the conductivity according to the thickness of the plate.θ : is the orientation angle of the fibers in relation to the chosen reference system.
In the reference linked to the fibers of each ply (L, T, Z), the Ohm law is written: Where and ⃗ When the fibers are oriented at some arbitrary angle θ, (as shown in Figure 2 (b)), with respect to the x-axis, the conductivity matrix is no longer diagonal and there is cross coupling of the components.To determinate the new conductivity tensor, we make use of the rotation matrix ̿ which relates the components in the principal and fibers linked coordinates system, [19]: By replacing equation (4) in equation (2), we get: We therefore deduce In this case, the homogenized conductivity tensor of a CFRP plate is given by the following expression:

Eddy current testing through Maxwell equations
Maxwell equations describe the physical model used for electromagnetic EC problems solved with the finite element method (FEM), that allows determining the response of eddy current sensor.The magnetic vector potential, electric and magnetic field or the pointing vectors are the most quantities widely used to solve the field equations.Quasi-stationary Maxwell equations are given hereafter: The electric vector potential formulation ( ⃗ ) is based on the conservation law of the electric current which is given by: By replacing equation ( 12) into equation ( 9), we obtain the following relation: Where ⃗ ⃗ is the intensity of the magnetic field.It is given by the following expression: Such as ⃗ ⃗ and ⃗ ⃗ are the intensities of the magnetic field in the source and in the piece respectively.The developed formulation is based on the ⃗ formula given by equation ( 13), and the Biot and Savart law which is expressed by the equation below [14,15]: Whose r is the distance between a point situated in the source region, and a point in the region that represents the plate.
is the volume in a source region.When replacing equation (15) into equation (13) we can obtain the following expression: ⃗ ⃗ can be calculated by Biot and Savart law as follows: represents the distance between two elements in the conductor and is the volume of each element in conductor.By replacing equation ( 17) into equation ( 16) we get the following expression: To solve this equation, we have used the ⃗ -iterative method which was already proposed by some authors to calculate eddy currents in NDT problems [13][14][15][16].
The advantage of this method is that it uses a reduced number of unknown variables and the air is not included in the mesh of the global system.
The flowchart of this method is illustrated in the Figure 3.
In the first step, the induced current density is initialized to zero in each meshing element of the piece.Then, in the second step, we proceed to solve equation (18) in an iterative process using the equation ⃗⃗⃗⃗⃗⃗⃗⃗ ⃗ to calculate the new density of induced currents.In the third step, a convergence test has been carried out on all the real and imaginary parts separately in order to avoid some problems of order of magnitude between both.
To control the convergence, a relaxation is used which is expressed by: Hence, the introduction of a factor α called "relaxation factor" such as " ", [16] aims to make it possible to ensure system convergence.The maximum convergence speed should be achieved with .The algorithm may usually diverge, so a damped algorithm choosing a parameter improves stability (convergence) but at the cost of higher number of iterations, leading to a longer computation time.The process stops when the convergence is reached.

Impedance variation calculation
There are several formulations in the literature for calculating the impedance variation where the inductor is fed by a current Is, the variation of its impedance due to the currents induced in the piece for a conductive material can be evaluated by the following expression [15,20,21]: Where ⃗ represents the opposite electric field produced by the induced currents in the piece and represents the source current density in the coil.Similarly, ⃗ represents the electric field in the source and expresses the density of induced currents in the piece.By using the Biot and Savart formula, we can calculate the electric field ⃗ .Finally; we get the following relation [9]: Where | | is the distance between a point marked by the vector belonging to the region and a point marked by the vector belonging to the region .

Validation of developed code
First of all, we validated the developed code by using the ( ⃗ -iterative) formulation.This method has the advantage that we only link the active zone compared to the other formulations ( formulation and the ⃗ formulation).
In this application, we consider the system described in Figure 4, which represents a coil excited by an alternating current placed in the vicinity of a carbon fiber composite plate having the characteristics described in the table below (Table 1).Figure 5 depicts the meshing of the plate and the coil respectively and Figure 6, whose objective is to compare the results obtained by the proposed method with analytical ones represented in ref [17].

Figure 4: Geometrical configuration of the validation device
The simulation of any electromagnetic system requires knowledge of all physical and geometrical characteristics in different regions.
The physical and geometrical parameters of the system to be studied are given in Table 1: The use of the ( ⃗ -iterative) formulation allows meshing the plate and the coil separately only once with a hexahydric mesh element.
The coil is divided into 2560 volumes and the plate is divided into 40000 volumes.Air is not included.Figure 5 illustrate the mesh of the plate and the coil respectively.In this section, we will calculate the impedance variations of the sensor as a function of different frequencies for a single fiber orientation θ = 0° as indicated in Figure 6.Through this figure one can remark a large agreement between the numerical results and those of analytic ones on all frequency intervals, the difference between them is denoted by a mean absolute percentage error of 2.9952%.
Figure 7 shows the distribution of eddy currents in a unidirectional carbon fiber plate for different fiber orientations.The induced currents circulate on a small section with a high density in the direction of the high conductivity σ L , and on a wider section with a lower density in the direction of the low conductivity σ T [8,15].

Application
In the literature there are several forms of sensors in the field of EC-NDT for the inspection of isotopic or anisotropic parts [6,8,15,16], the most important in the domain of composite materials is the rectangular sensor [15], Figure 8, since the orientation of the fibers in each ply constituting the CFRP plate is determined directly from the analysis of the normalized impedance variation of the sensor (see equation ( 22)) as a function of its angle of rotation.In this paper, three different applications are proposed to illustrate the characterization of a unidirectional, bidirectional and multidirectional carbon fiber plate and to show the effect of fiber orientation using the normalized impedance variation diagram in the polar coordinates.
The normalized impedance variation is calculated by the following relation: The different physical and geometrical characteristics of the system (plate/coil) are presented in Table 2.The studied system has been meshed by using a non-regular hexahydric mesh refined in the center of the plate, as shown in Figure 9.Among the major problems in the characterization of composite materials is the knowledge of fiber direction.In this part of this paper, we used the normalized impedance variation polar pattern to determine the orientation of fibers in a plate carbon fiber.This diagram is one of the most effective means for the characterization of composite materials, because of its simplicity in determining the position and fibers orientation in a plate with one or more plies with great readability.Firstly, we shall study the case of an anisotropic unidirectional plate.The second example is depicted for an anisotropic bidirectional plate and the third application is destined for an anisotropic multidirectional plate.

Characterization of a unidirectional anisotropic plate
In order to have the configurations depicted in Figure 10, each time the normalized impedance variation is plotted as a function of the rotation angle of the sensor, it can be seen that the shape of the normalized impedance variation takes the form of large lobes due to the higher conductivity in each direction (θ = 0 °, θ = 45 °, θ = 90 °, θ = 135 °), see Figure 11, which led us to know the exact direction fibers in each plate.

Characterization of a bidirectional anisotropic plate
Figure 12 shows different configurations that will be studied in this section, consider an anisotropic two-ply (bidirectional) CRFP plate, whose plies have several fiber directions (θ= 0°, θ = 45°), (θ = 0°, θ =90 °) and (θ = 0°, θ = -45°) respectively, according to the Figure 13.It can be seen that the shape of the variation of the normalized impedance as a function of the rotation angles of the sensor takes the form of two lobes, each of which is directed in a direction the maximum one (ply above) and the other minimum (ply below), which guarantees the presence of two plies in the carbon fiber plate, on the other hand the precise knowledge of the orientation of the fibers and the position of the plies in each plate.

Characterization of an anisotropic multidirectional plate
In this third application, we consider a four-ply (multidirectional) anisotropic CRFP plate, of two configurations of plies (see Figure 14), (θ = 0°, θ = 45°, θ = 90°, θ = -45°), (θ = 90°, θ = 45°, θ = 0°, θ = -45°) respectively.From Figure 15, it is found that the shape of the variation of the normalized impedance as a function of the rotation angles of the sensor takes the form of four lobes oriented in the directions of the fibers in the four plies.The amplitude of each lobe gives the position (depth) of the ply which corresponds to it.This ensures the presence of four plies in the plate.These last ones varied between four amplitudes in the polar diagram.

Conclusion
Throughout this article, a benchmark problem of a nondestructive evaluation has been solved.We have developed a calculation code based on the ⃗ -iterative technique associated with the finite element method.The Biot and Savart law is a means we have used to avoid meshing the air when moving the coil.Eddy currents were evaluated in the plate for showing the fiber orientation.The obtained simulation results were compared with the analytics ones.A great concordance has been noticed between these results and therefore justifies the utility of this method.In addition, the advantage to use the rectangular shape sensor, that at each position (rotation) gives us a new impedance variation value different from the previous one, so it clearly shows the direction and position of the fibers.The yielded results reveal that one can characterize the plate and determinate the number and position of the ply, hence the direction (orientation) of the fibers in the composite plate.

Figure 1 :
Figure 1: Different parts of the studied system

Figure 2 :
Figure 2 : a) Structure of a composite material, b) Principal and fibers-linked references

Figure 3 :
Figure 3: Flowchart of the proposed iterative method

Figure 5 :
Figure 5: 3D mesh of the studied device

Figure 6 :
Figure 6: Impedance variations of the sensor as a function of different frequencies

Figure 7 :
Figure 7: Eddy currents repartition in the plate for different orientations (a) θ=0°, (b) θ=45°, (c) θ=90°, (d) θ=135°) From the above figure, one can notice the effect of anisotropy.The induced currents circulate on a small section with a high density in the direction of the high conductivity σ L , and on a wider section with a lower density in the direction of the low conductivity σ T[8,15].

Figure 8 :
Figure 8: Geometrical configuration of the studied device

Table 2 :
Physical and geometrical characteristics of a studied device

Figure 9 :
Figure 9: 3D Mesh of the multidirectional plate

Table 1 :
Physical and geometrical characteristics of a considered device

Table 3 :
CPU time for each problem