Direct Resistive Heating Simulation Tool for the Repair of Aerospace Structures through Composite Patches
Bonded composite patches are very appropriate for aircraft structural repair, showing very good properties when compared with traditional mechanical fastening of metal sheets. The curing process of the composite patch must be done “onsite” and a direct resistive heating method has been proposed. The heat produced by the electric current through the Joule effect when crossing the patch carbon fibre bundles has been modelled with a Finite Element Program code, developed for the electric current equation. The heat conduction equation has also been modelled in the program, as well as the kinetics of the curing reaction of the composite resin. The electric resistivity of the materials is function of the temperature, so a nonlinear coupled system has been developed. Therefore, a complete simulation tool able to study different configurations, current intensities, materials, etc. for the composite patches is presented. A study case has been run with the developed code and results have been compared with experimental values. Good agreement is found.
It is well known that small cracks tend to appear in different structural parts of aerospace structures mainly due to the cyclic loads that these components undergo. The measurement and control of these small cracks is an important part of the day to day maintenance activities of aircraft. When the length of the cracks reaches certain predetermined values, it is necessary to repair the structural part (wings, tail, etc.) which means that the whole aircraft must remain for some time on earth for this operation. It is clear that these repairs account for an important quantity of money, not only by the repairs themselves but also for the nonavailability of the aircraft for passenger or cargo services.
Traditional repair of these defects has been accomplished through mechanically fastened patches. However, adhesively bonded patches appear to be a good alternative, as they have a large contact area for the bonding and a very good value of specific stiffness. Composite patches also prevent crack initiation and propagation, have good fatigue performance, prevent corrosion, and conform more easily to complex geometries. The application of the composite patch on the aircraft involves either the precuring of the patch and the bonding on the substrate or the on-aircraft curing and bonding of the patching system (hot bonding). If the damaged component cannot be removed and placed in an autoclave, hot bonding seems to be almost mandatory. To provide compaction, vacuum bagging systems are used. For out of autoclave repairs, heat from portable equipment can be delivered through heat blankets, lamps, or hot air.
It has also been proposed  that electrical resistance curing could be used for carbon fibre/epoxy composites. The idea is to pass electrical current through a layer of patch carbon fibre bundles and make use of the heat generated by the Joule effect to reach the curing of the resin and adhesive present in the patch. This novel procedure seems to bring some advantages when compared with other methods: heat is distributed more evenly in the whole patch, and no specially tailored devices are needed, so a more uniform temperature distribution could be expected, which is very convenient for a uniform curing process.
Modelling of the curing process of composites through numerical methods has already been tackled by different authors. Yi, Hilton, and Ahmad  presented a complete approach through the Finite Element Method (FEM), considering the nonlinear heat transfer equation and the cure kinetics equation, with 2D examples. Joshi, Liu, and Lam  studied the curing process through a commercial heat transfer FEM code with some user developed code for the cure kinetics, employing nodal control volumes. Park, Goo, Min, and Yoon  presented a complete FEM formulation valid also for materials with different thermal conductivities in different directions, with special emphasis on the 3D nature of the process. The use of resistive heating elements for the curing process has been studied, along with the curing kinetics, by Zhu and Pitchumani . They solve analytically the 1D through thickness heat transfer equation with nondimensional numbers and make a complete study of processing windows for different parameters. Companion information about the experimental work is provided by Ramakrishnan, Zhu, and Pitchumani . Resistive heating curing experimental work in aerospace applications with a very simple heat transfer model is presented by Rider, Wang, and Cao  and the prediction and measurement of electrical conductivity for carbon unidirectional layers are reported by Athanasopoulos and Kostopoulos . The study of the cure kinetics for aerospace employed resins has been reported by Kim, Moon, and Howell  and Shin and Hahn . They provide, through Differential Scanning Calorimetry (DSC) experimental work, values for the cure kinetics equations. Finally, the study of the cure cycle influence on the mechanical properties and curing kinetics through a commercial FEM code and users subroutines is presented by Zhang, Xu, and Huang .
It is interesting to note that in all the simulation work previously reported only the nonlinear transient heat transfer equation and the cure kinetics equation have been used. Even for the work where resistive heating was employed [5–7], the Joule effect power was considered constant through the whole process. However, Athanasopoulos and Kostopoulos  report that the electrical resistivity depends on the temperature, although this fact is not taken into account in their modelling work because it is outside the scope of their investigations. It is also reported in the experimental work of Joseph and Viney  that there is a change in the voltage in their resistive heating setup (it goes from 6.0 to 6.4 V) when the current intensity is maintained constant with a 15 A value. This means clearly that the electrical resistivity of the composite material changes with the temperature and that it would be very interesting to take into account in the simulations the electric current equation also, coupled with the heat transfer and with the cure kinetics equations.
The role that computer modelling can play to determine particular values for imposed intensity electrical currents in particular geometries, expected curing times, foreseen maximum temperatures in the composite patch during the curing process, etc. is quite clear. Different alternatives can be studied in detail if a simulation tool is available and different materials behaviour could be compared for a specific needed repair. As a consequence, a numerical simulation tool based on the Finite Element Method has been developed to model the curing process of composites patches used in aerospace repair.
2. Theoretical Background
The equations that govern the curing process of the composite bonded patch are the electrical current equation (Ohm’s law in its most simple form), the conduction heat transfer equation, and the cure reaction kinetics equation. The three of them are detailed below with some comments for the particular specific case that is being studied.
Electric Current Equation. For simplicity, it has been considered that the material formed by the fibre bundle and the resin is homogeneous. In this context, electrical properties are the same in all parts of the patch itself, although different material properties are accepted in different directions (orthotropic). With this in mind, the equation that governs the electric current iswhere stands for the voltage drop [volts], is a current source [amp/m3], and is the electrical conductivity [Siemens/m] in the direction (inverse of the more common known value of electrical resistivity ). In the present process, is not present, and the equation is clearly stationary (no time dependent). It is possible however that the electrical conductivity is a function of the temperature of the material , which is in itself a source of nonlinearity for the equation. When solved through the Finite Element Method, this equation gives voltages in each node and current density [amp/m2] in each element and each direction. Once both are known, it is fairly easy to calculate the energy by unit volume and unit time [W/m3] that is produced by the Joule effect.
Heat Transfer Equation. The same simplifying consideration about homogeneity that has been accepted for the electric current equation applies for the heat transfer equation as follows.where stands for the temperature [°C], is a heat source [W/m3], is the specific heat [J/°Cm3], and is the thermal conductivity [W/m°C] in the direction. This equation is transient, that is, time dependent, as it is the curing process itself. The thermal conductivity of the material can be itself function of the temperature, which is a source of nonlinearity. Another source of nonlinearity is the boundary condition of radiation, where the heat flow interchanged by the material with the surrounding is function of the fourth power of the absolute temperature . When solved through the Finite Element Method, this equation gives temperatures in each node and heat fluxes [W/m2] in each element and each direction. These results will be provided for each time step. It is also interesting to note that the term will have two parts, one coming from the power generated by the direct resistive heating (Joule effect) and another part coming from the curing heat. Both terms must be calculated and used in the resolution of the heat transfer equation. Finally, it is also worth pointing out that the formal appearance of both electric current and heat transfer equation is the same ( and are equivalent as and are), which simplifies the computer programming of the method.
Curing Kinetics Equation. This equation has a very different aspect, as can be seen in (3), when compared to (1) and (2).where is the resin curing degree [no units], is the curing rate , is the absolute temperature [K], and , , , and are materials constants, found experimentally. Equation (3) is completed by (4), which shows the time integration scheme adopted and also completed by (5) which shows the generated power [W] as a function of the curing rate , the curing degree [no units], the fibre volume fraction [no units], the resin density [kg/m3], the total volume [m3], and the resin cure reaction heat [J/kg].
3. Numerical Implementation
The discretization and solution of the differential equations that govern the curing process of the composite bonded patch have been done through the Finite Element Method. The developed program has been named JOULE.
Finite Element Discretization. The resulting electric current equation is shown in (6), with the conductivity matrix expression provided in (7) and the nodal electric intensity vector in (8)Apart from values previously described and from the gradient operator and the interpolation functions , values of the electrical conductivity matrix and imposed surface electrical fluxes appear in (7) and (8).
The resulting heat transfer conduction equation is shown in (9), with the thermal inertia matrix expression provided in (10), the thermal stiffness matrix in (11), and the nodal flow vector in (12).Apart from values previously described, from the gradient operator and the interpolation functions , values like specific heat per unit volume , the thermal conductivity matrix , the convection coefficient , outside temperature , and imposed surface thermal fluxes appear in (10), (11), and (12). Radiation boundary condition is treated as a convection boundary condition with a coefficient which varies with the temperature. This implies a nonlinear heat transfer calculation, if present.
It is worth mentioning that the cure kinetics equation can be directly numerically integrated with the proposed scheme in (4) and that there is no need of employing the techniques shown for the electric current and heat transfer equations. In fact, at the beginning of each time step a space distribution of temperatures is known, as well as a space distribution of cure degrees and curing rates. The power generated in that time step is calculated through (5) and used to calculate the converged distribution of temperatures at the end of the time step, when updated values of the curing rates and cure degrees will be calculated for each node in the FEM mesh.
The power generated through the Joule effect can be expressed asThe dot product of the electric field vector and the current density vector expressed in (13) is calculated for each integration point for each element, which is a very convenient scheme, as the values can be used directly when solving the heat transfer equation.
The time step employed for the integration must be small enough to capture adequately the curing power. When the kinetics of the reaction imposes a very high curing rate, high values of power are calculated and, as a consequence, it is very convenient to decrease the size of the time step. It is not easy to foreseen before doing any simulation in what time step the curing rate will increase rapidly, so a posteriori check of the calculated heat produced by the program and the one that could be easily calculated by hand (total volume × density × curing heat) should be made. A second simulation with convenient nonequal time steps is then recommended; if error is bigger than certain limit (we have used a 1% error criteria in our simulations), time step should be decreased.
Solution Procedure. Figure 1 shows the flowchart to solve the set of coupled equations.
It is clear that there is a coupling between the electric current equation and the heat transfer equation as one of the inputs for this last equation (the power generated through the Joule effect) comes from the electric current equation. But it is also true that there is also a coupling between the electric current equation and the heat transfer equation when the electric properties of the materials are nonlinear, that is, when electric conductivity is known as a function of the temperature. This bidirectional coupling implies that the staggered solution scheme presented in Figure 1 makes two convergences checks, one for the thermal problem (with its nonlinearity) and another for the coupling with the electric problem. As a consequence, a significant number of equation resolutions must be made for each time step before a consistent set of temperatures is reached as a solution. The numerically intensive nature of the solution of the problem has been addressed by a very careful minimization of the disk writings, a dynamic RAM memory scheme for the management of data, and the equation solution and finally for a specific subroutine that takes advantage of the possible linear parts of the problem that do not need a convergence calculation if, for instance, material properties are not function of the temperature.
4. Case Study: Resistive Heating Curing Test of a Rectangular Composite Sample
It has not been possible to test a real composite aerospace patch at this stage but several direct resistance trials have been carried out in Tecnalia (see Figure 2). The samples included a fibreglass reinforced laminate to avoid electric current and heat leaks, a CFRP, composed of 5 carbon layers (HexForce® 43280 S 1070 TCT) (about 50% of filler content) and the Epocast 52-A/B epoxy and an insulation blanket to avoid heat leaks from the top. Doped (1% CNT) and nondoped laminates were built during the testing stage.
Concerning the measuring devices, two electrodes were connected at the edges of the CFRP and an electric current was applied. For temperature measurements, different thermocouples were placed. Cross section of the sample can be seen in Figure 3.
Regarding the reproduced test, a particular one where an electric current of 14.9 amperes was applied has been chosen. As shown in Figure 3, 3 thermocouples were placed on the CFRP: one on the bottom (A), the second in the middle (B), and the last one on the top (C). Finally, the specimen was inserted into a vacuum bag and another thermocouple (D) was placed on the top. The dimensions of the coupon are 250 x 110 x 1.615 mm.
5. Computer Models
Three different computer models have been generated for this case study, namely, the 2 DOF (Degree of Freedom) model, the 2D model, and the 3D model. Temperature results will be shown for the 2 DOF model and the 2D model, but only the 2D model will be studied in detail. The reason is that the 2 DOF model gives interesting preliminary results but it is not accurate enough to describe the real situation and the 3D model requires an important computer work with almost no improvement of the results.
2 DOF Model. This model takes into account the resistive heating, the exothermic curing effect, and the boundary conditions of convection and radiation to the surrounding air. All the components of the system (vacuum bag, Teflon, composite, and insulation blanket) are assumed to have a single temperature, which is the one calculated and shown in Figure 5. No FEM model has been used here, but only the heat transfer, electric current, and curing equations. Results are good enough for an initial approximation, as they give a good idea of the time moment of the maximum temperature in the system.
3D FEM Model. This model accounts for the whole coupon and a mesh of 15 000 eight node elements with a total of 17 776 nodes was created. The idea of this big model was to study whether the width of the coupon had an influence in the temperature levels reached in the middle of the coupon. Simulations showed that almost no difference of temperatures was present between the 2D and 3D models, so this 3D model was no longer studied.
2D FEM Model. The mesh shown in Figure 4 reflects the 2D FEM model, where the vacuum bags are represented by the red elements, the insulation blanket by the blue elements, the Teflon layers by the green elements, the CFRP by the yellow elements, and the fibreglass reinforced laminate by the pink elements.
Only one part of the employed mesh is shown, as the length of the case (250 mm) is much higher than the coupon thickness (1.615 mm) and the total thickness (5 mm). It is also necessary to mention that a 2D calculation has been done, as there was no variation in the width direction of the behaviour of the sample. In fact, a 250 x 5 mm rectangle has been drawn and divided into 7 500 four node elements with a total of 8 016 nodes. The width of the coupon has been included in the input file. It has been considered that the vacuum bag had a 0.1 mm thickness (1 layer of elements), the glass fibre 0.8 mm (2 layers of elements), the Teflon 0.17 mm (1 layer of elements), the CFRP and Epoxy 1.615 mm (4 layers of elements), and the thermal insulation blanket 2.045 mm (5 layers of elements).
It is interesting to note that electrical conductivity as a function of temperature was used for material number 4, the dry fibre plus epoxy resin. It is also worth noting that the thermal conductivity of the dry fibre plus resin is ten times lower in the vertical direction than in the on-plane direction.
For the cure kinetics calculation, (3) has been used with the following values: , , , . These data come from experiments made at Tecnalia. For the heat power generated during the cure process, (5) was used, with the following values: , , . Concerning the boundary conditions, a convection coefficient and an emissivity value for radiation of were used. Ambient temperature was fixed at °C according to measurements in the test room.
6. Results and Discussion
Figure 5 collects the experiment measured temperatures and results from simulations. It is important to say that the thermocouple B (the one in the middle of the composite, see Figure 3) had a malfunction during the test and, consequently, the recorded temperatures were wrong. For the test temperatures shown in Figure 5 for the point B, we have calculated the arithmetic mean between the temperatures recorded by thermocouples A and C. The difference between the temperatures in these two thermocouples was always lower than 2.5°C (at the same time step), which justifies the decision made.
Figure 5 also collects temperature values coming from the above explained 2D simulation. Companion results from the very simple preliminary two nodes model (2 DOF model previously explained) are also represented in Figure 5. The results obtained for the 3D simulation have not been included in the graph due to the fact that there are almost no differences when compared with the results of the 2D model.
The results of the 2 DOF model can be considered good for initial estimates of needed time of curing for a particular test geometry and electric current value. The difference between the temperatures shown by this preliminary model and the ones obtained with the 2D model are due to the fact that a single temperature is representing a set of different materials and their different thermal and electrical properties. The 2D FEM model predicts quite well the maximum temperature and the time when it happens for the three points (A, B, and C), although it overestimates somehow the temperature velocity gradient.
Regarding the variation with time of the voltage drop measured and predicted, results can be seen in Figure 6, where the change of electrical conductivity with temperature is the main responsible for the shape of the curve. The general trend of the voltage drop has been captured well but the real values are somehow different. This mismatch is attributed to the difficulty to have real electrical conductivity values properties as a function of the temperature.
Finally the expected values for the variation with time of the curing rate and the curing degree in the middle of the coupon are presented in Figure 7. Total cure (higher than 90%) is predicted for 1.5 hours of process, while the maximum curing rate is expected at time 1 000 s. These values are consistent with what was observed in the test.
A methodology to simulate the direct resistive heating of aerospace composites patches for onsite hot bonding trough the Finite Element Method has been presented. Tests results show good agreement with predicted values.
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
This work was developed under the European Seventh Framework Program, Theme 7 Transport, Project IAPETUS [Grant Agreement no. ACP7-GA-2008-234333].