Focalization of electromagnetic power at the interface between two composites materials for induction welding

This paper presents a fast-numerical approach to optimize the localization of induced power in highly anisotropic composite materials with a large-scale factor. A one-dimensional hybridmulti-physical tool allows us to find the best ply sequence and the best strategy to constrain the induced currents to circulate in an interest area defined respecting specific temperature objectives. A three-dimensional multi-physical modeling is then used on a few selected candidates to refine results.


Introduction
Electromagnetic induction welding of carbon fiber composites with thermoplastic matrix is a very promising technology for aircraft and space industry over the coming years. Compared with existing processes, it requires less energy, and the cycle times are shortened which allows us to overcome some issues associated with classical assembly techniques [1]. Figure 1 shows a classic configuration of a structural part of the aircraft (fuselage). The skin (composite 2) is reinforced by omega (composite 1). The skin is the bottom plate; omega is the upper piece in shape of omega which strengthen the structure. To obtain a good assembly, it is necessary to heat the interface up to the melting point of the thermoplastic matrix (350°C). To avoid material deconsolidation, the temperature of the remaining volume is kept below the target temperature. To achieve this objective, the induced electromagnetic power must be mainly located nearby the interface area.
Composite carbon fiber material is a piling up of unidirectional plies layered with different orientations (Fig. 2). In each ply, the electrical conductivity ratio between the conductivity along fiber direction and transverse direction is higher than s // /s T > 10 4 . No flow of induced current loops in the same ply during the induction phenomenon, if not the induced power would be extremely low. Due to electrical percolation phenomena, a current path between plies is then established and could be different according to the ply sequence [2].
Induction heating ensures fast and contactless injection of electromagnetic power into composite materials and paves the way to new processes in technological breakthrough either for assembly (welding) or for the evaluation of material health (non-destructive testing) [2]. However, due to the overly complex nature of these materials, the absorption zone of the electromagnetic fields is still too diffuse. To have a better control of electromagnetic power location, we have three degrees of freedom: the orientation of the plies, the addition of an insert material at the interface, and/or the insulation between plies (i.e. to have a degree of freedom on the conductivity of the plies). To investigate these possibilities, modeling tools are required. A 3D multi-physical model is extremely relevant but is very time-consuming because many cases must be investigated. In this paper, an original approach based on a coupling between 1D and 3D models is presented. The first one helps investigate different families of possibilities to optimize the ply sequence and the location of the power to realize the desired functionality. A three-dimensional magneto-thermal modeling [3] is then used to optimize pre-selected configurations.

1D model description
This study involves the electromagnetic and thermal modeling of the system. Two models are thus released: the electromagnetic model that goes through the recall of Maxwell's equations and the behavior laws of materials and the thermal model that is defined by the heat diffusion equation.
For 1D modeling, the welding system in Figure 1 is reduced to a 1D system as shown in Figure 2. The layer colored in blue is the interface between the two composites. It represents the insert layer in the case of an insert addition.

Electromagnetic model
The magnetic field at a layer (fold or ply) k of the configuration is governed by equation (1) [4].
with m k the magnetic permeability, s k the electrical conductivity of the layer k and f the frequency. Since the source field H 0 is assumed to be tangential and the plies infinitely long, (1) admits a solution in the form in (2).
where g k = (1 + j)/d k with the skin depth d k ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1=ðpfs k m k Þ p . The constants A k and B k of a layer k are calculated according to electromagnetic transition conditions between layers k + 1 and k − 1. For example, between the layer k (1 < k < n) and the layer k + 1, we have the relations in (3).
with E k = ∂H k /(s k ∂z k ) the electric field at a layer k. At the borders of the system (layers 1 et n), we have the conditions in (4).
The resolution of (2) gives the value of the magnetic field H at any point z of the system. Starting from these values, we can deduce all the useful quantities such as the induced power density q and the induced current density J defined by equation (5).

Thermal model
The first principle of thermodynamics translates the energy conservation equation (Eq. (6)) which defines the thermal behavior of the material as a function of time tt [5].
with r the specific mass, Cp the specific heat of the material, l its thermal conductivity and T the temperature at the given point. In this study, the temperature T ð tt z Þ varies along a single axis z and the time tt, thus (6) becomes (7).
∂lðzÞ ∂z : ∂T tt z ∂z þ lðzÞ: This equation is solved by the finite difference method which considers all the non-linearities of physical properties. This resolution passes through the discretization of the layer thickness in Dz and the time in Dtt. Which allows us to write the relations in (8).
To find solutions to this problem, it is necessary to set boundary conditions. In our case, the conditions are heat convection type defined in (9) and explained in [5]. h is the thermal convection coefficient (Wm À2 K À1 ) and T a the ambient temperature.

3D formulation
The phenomena related to electromagnetic induction welding are multi-physical and multi-scale. Moreover, the important scale factor between the macroscopic scale and the microscopic scale of composite materials could lead to important calculation cost. The 3D computational codes in A-' formulation (Eq. (10)) developed in [2] are used for the electromagnetic simulation of our system.
where [s] tensor is obtained with the help of homogenization techniques. C and G are the discrete counterparts of curl and grad operators,  Equation (6) is used to solve the thermal problem in the composite plate with the Fourier boundary condition defined in (9). The heating problem is solved using a 3-D anisotropic finite-element method with different thermal conductivities along the fiber axis and transverse fiber axis.
To allow a reasonable comparison between 1D and 3D model, a U-shaped inductor, following Figure 3, is used to obtain a tangential field under the active part of the inductor.

Focusing of electromagnetic field
In this section, we study the impact of plies orientation, plies insulation, and addition of insert on the welding process. 1D model allows to pre-optimize the configuration. 3D model allows us to refine the results.
First, it is considered that the composites 1 and 2 consist of 15 unidirectional plies successively oriented at 0 and 90°([0/90/0/90/0/90/0/90/0/90/0/90/0/90/0] or [[0/90] 7 /0]). The idea here is to compare 1D and 3D models. A comparison between 1D and 3D models shows the same order of magnitude and physical behavior shown in Figure 4. We notice that all the power density is in the 0°p lies. This power distribution prevents any welding between the two composites. To obtain a good welding, a ply sequence where the power is more localized nearby the interface is aimed.

Plies orientation
Our first investigation of the electromagnetic power focusing consists modifying the ply-sequence of the composite material. Because the 1D approach is extremely fast, all possible configurations could be investigated with an automatic process. As a result, a candidate that seems serious to us for a focus of power at the interface of the two composites is obtained. The ply sequence is the same for both composites and is ½0 2 =45 2 = À45 2 =90 2 s . 0 2 means two adjacent layers oriented 0°from X-axis and subscript "s" means the next 8 plies has symmetric orientation from previous ones. Figure 5 shows its electromagnetic response with our one-dimensional tool.
Here, the power is brought to the interface between the two composites, but a large part of the power is also located at the outer ends of the composites.
With the 3D tool, the same configuration is simulated. The observation remains the same as illustrated in Figure 6. There is some difference between the power distributions with the 1D tool and those with the 3D model. Nevertheless, with the results of 1D we know what to expect from 3D.
This investigation track has its limits because the power located on the plies on the outer ends of the two composites is especially important. To overcome this problem, the modification of the local properties of the composites is investigated in the following section.

Insulation of the plies
The idea here is to insert thin layers of resin (with a thickness of 10 mm) between the plies located at the outer    ends of the two composites. Figure 7 shows the electromagnetic power distribution obtained in this case. We can notice that the resin layers have forced the induced currents to flow differently. In this way there was a clearer focus at the interface of the two composites. Figure 8 shows the thermal response of this configuration. The temperature is well centered at the interface, but the welding width is very wide.

Addition of an insert
Our last investigation is to add a highly conductive layer with target properties between the two composites to be welded. The insert must be able to concentrate the heating power, thus permitting the welding of the two composites. The composites are like those in Figure 4. To find the right properties of the insert an optimization simulation was done with our 1D tool. The optimization variables, in this case, are the magnetic field H0, the electrical conductivity s, the magnetic permeability m of the insert, and the convection coefficients h1 and h2 respectively of the inductor side and the external side of the composite 2. Playing on the convection coefficients means that a cooling system can be added to each side of the plate composed to the two composites. A coefficient of 10 W/m 2 /K means that the corresponding side does not need a cooling system, so it can be left in the open air. NSGA-II, Non-dominated Sorting Genetic Algorithm II, is used for better robustness [6] to meet two contradictory objectives. On the first hand, it is necessary to have a temperature of 400°C ± 10 in the insert layer (f1 = 1000* (T − 400) 2 ) and, on the other hand, to have a temperature that is as low as possible in the plies adjacent to the insert (f2 = exp(3* (T/100 − 2)).
The realization of this optimization is done using MATLAB. After 5201 number of function calls and 103 generations, a Pareto front is generated where the different optimal combinations of the optimization variables are highlighted. In these series of combinations, the one that best meets our objectives is chosen (Tab. 1). By exploiting the data in [7], an insert with these target properties can be obtained with cobalt alloys.
With these optimized parameters of the insert, we simulate the configuration shown in Figure 2. All the electromagnetic power is located at the insert layer as shown in Figure 9. This is justified by the superior properties of this layer. Figure 10 shows the temperature map according to thickness and time. At a heating time of 10 seconds, the temperature in the insert layer is 392°C. This temperature is the highest in the studied system and does not exceed the material degradation limits (<400°C). This is in the image of the power induced in this layer which    represents 81% of the total power induced as shown in Figure 9.
To refine the results found, a simulation with the 3D tool developed was done. The current of the inductor is set to have the same magnetic field profile in the center of the plate as in the one-dimensional case (Fig. 3). The thermal response of this simulation is shown in Figure 11.
We can see that the curve obtained follows the same logic as that in 1D. The slight difference between the curves is mainly due to the shape of the inductor, the phenomena of edges, and the hypothesis made in 1D. Indeed, the magnetic field on the external phase of the composite 1 is assumed constant and tangential to it, which is not the case in 3D simulation. Nevertheless, a good concordance of the results is obtained.

Conclusion
In this article, a 1D hybrid method (analytical for the electromagnetic part and digital for the thermal part) for solving electromagnetic induction heating problems was proposed. The different results allowed us to see the robustness and the speed (1.86s by iteration) of this method coupled to an optimization algorithm. This tool is very effective for predetermining the pair (s, m r ) of the insert, in particular for a very large number of plies. Nevertheless, the use of a 3D numerical simulation is necessary to refine the results and consider the inter-plies phenomena.
Author contribution statement