Prediction of warping in thermoplastic AFP-manufactured laminates through simulation and experimentation

Abstract The thermoplastic automated fiber placement (T-AFP) process is a non-autoclave method for in situ consolidation of thermoplastic composite material on a piecewise constructed laminate. High thermal gradients and nonlinear material behavior, especially due to crystallization, make predictions of process-induced stress and warping difficult. This article describes a method for simulating parts manufactured by T-AFP using a detailed material model to capture the dynamic nature of the process. The material model is flexible and can be altered to describe different semi-crystalline matrices, in this study focusing on low-melt polyaryletherketone. Two laminate panels are simulated within this work and assess the impact of a heated tooling on overall part warping. Panel warping is validated by performing 3D-scans of T-AFP-manufactured laminates produced using the same parameters as the simulation. The results show a good match between numeric and experimental warping, especially for heated tools, thus, providing a useful method for predicting laminate warping and reducing the demand on manufacturing experimentation. Graphical Abstract


Introduction
The use of fiber reinforced plastics (FRP) was steadily increased over the past decades. Advantages in strength-to-weight ratio and the possibility to align the load carrying fibers along the load paths separate FRP-structures from traditional metal constructions in lightweight designs. The primarily used thermoset matrix entails relevant disadvantages. The required time to process the matrix is very long and can take multiple hours or even days. Additionally, most thermosets are highly toxic and carcinogenic, requiring special personal protective equipment (PPE). These properties have slowed widespread introduction of FRP into lightweight structures.
Thermoplastic matrix materials can alleviate these challenges. They are processed by heating the prepreg material above the matrix melting point. After that, they can be formed and welded depending on the requirements and finally cooled down to solidify the matrix again. Unlike thermoset materials, which keep their form after curing, this process can be repeated multiple times, so parts can be remoulded and recycled [1]. Additionally, thermoplastic materials are mostly nontoxic and can be stored for a longer duration. As a disadvantage, high-performance thermoplastics usually require higher processing temperatures than thermosets and often include semi-crystalline properties, which induce additional stress and warping during cooldown.
An innovative non-autoclave method for processing thermoplastic composites is the automated fiber placement (T-AFP) with in situ consolidation (ISC).
In T-AFP a prepreg tape of fiber reinforced thermoplastic is heated above the matrix melting point fusing it onto the preceding layer with a compaction roller ( Figure A1). In the case of ISC no further consolidation steps are necessary, thus, reducing both manufacturing time and initial investment cost.
However, the high heating and cooldown rates and the inhomogeneous temperature distribution inherent to the T-AFP process induce high stress into the part. This is particular true for semi-crystalline matrix materials, with crystallization kinetics inducing high warping of the part, and hence, counteracting the advantageous effect on matrix strength properties. Qureshi et al. [2] reporting increased lap shear strength and double cantilever beam test performance with higher tool temperature. The exact temperature of the tooling depends on the material but needs to stay below the melting point. The part can cool down and solidify, but the temperature differences, and thus, the cool down speed are lower. Crystallization processes have more time to complete and the temperature distribution along the part is less extreme. The tool heating needs to be very homogeneous and has to incorporate heat bleeding at the tool edges.
Optimal part quality depends on the combination of all processing parameters, with Dreher et al. [3] describing the process of optimizing these parameters using DoE methods for an experimental setup. However, the cost and effort of conducting such an extensive experimental campaign is high, making a simulative alternative an attractive option. Multiple such approaches have been conducted in the past.
Stokes-Griffin et al. [4] conducted two-dimensional simulations of the thermal distribution during the process and compared their simulation results to temperature measurements. The model was expanded for determining the required laser heat flux for the desired temperature distribution [5]. These results were evaluated by determining the resulting bonding strength of the part. Another analysis was conducted by Tierney and Gillespie [6] with a focus on interlaminar shear strength. The bonding strength of two layers was simulated and compared to experimental results. Khan et al. [7] also created a thermal model and used it to investigate thermal degradation of the bonding zone. This expanded the work of Khan [8]. The thermal history of a part in two dimensions and its influence on the aging effects of the material was conducted by Regnier et al. [9]. While this work focuses on continuous welding of composites, the challenges of that process can be compared to T-AFP. The twodimensional modelling of the thermal history of a welded panel shows similarities to the thermal history of a flat T-AFP panel. Lemarchand et al. [10] finally investigated the residual stress of the T-AFP process in small scale two-dimensional simulations. The degree of crystallinity and its influence on residual stress is explained. Khan et al. [11] also investigated the influence of the tool temperature and the consolidation pressure on the final part quality in an experimental environment (i.e. without additional numerical support), while Schaefer et al. [12] used a design of experiments approach to optimize T-AFP process parameters with a fast 2dimensional simulation framework and Hosseini et al. [13] looked at the influence of non-uniform degrees of crystallinity during tape winding of pipes and predicted the resulting final part crystallinity. This work focused on the crystallization process though neglected its influence on the final part. In the work of Celik et al. [14] the three-dimensional temperature history of T-AFP parts, specifically in the contact zone between individual tapes, was numerically analysed and Zaami et al. [15] detailed the influence of different laser and manufacturing parameters on the process temperature of tape winding.
While the T-AFP process has been investigated in great detail at small scale and for individual process aspects (degree of crystallization and consolidation quality between individual tapes and layers) a feasible method for simulating large scale threedimensional warping of T-AFP panels has not been introduced. This work aims to close this gap by presenting a method for using a detailed material model in a simulative setup created to replicate the highly dynamic nature of the process at minimum effort. As an example, two panels manufactured on heated and cold tools are investigated.

Methods and materials
For the simulations, the commercial FEM-tool Ansys was used in the 2020R1 version. At first, the Ansys Composite PrePost-tool was used to create a three-dimensional lagrangian mesh of the laminate. The mesh has to contain at least one element per UD-layer in z-direction and at least three elements per tape width in the other spatial directions (compare Figure A2). In addition to that, a separate simplified mesh for the tool was designed, which only needed to contain one element in z-direction.
To limit the number of nodes, first-order elements with eight nodes per element are used. The analysis is split into a thermal and a structural part with the element type adjusted to the Ansys element types SOLID278 for thermal and SOLID185 for structural analysis.
The thermal simulation creates the thermal history of the part during the process. As such, it has to replicate the conditions during manufacturing.
The T-AFP process uses a heat source, e.g., a laser, to heat up a tape and weld it to the laminate. The simulation has to incorporate a moving heat source to simulate this effect. As a simplification, the heat source is presumed to be controlled to heat the tape up to the defined process temperature. In a given time step, the nodes occupying the current heat source position are forced on a heat-up temperature curve from tool temperature to process temperature. The duration of the time step is defined by the element size and the processing speed of the T-AFP machine. In the subsequent time step, the affected nodes are released from the forced temperature boundary condition, which is instead applied to the next nodes along the tape. If this moving boundary condition reaches the end of a tape, a short pause with no heating is inserted to simulate the time needed by the robot to move to the start of the next tape. The necessary geometric and temporal information are automatically calculated during the process by a pre-programmed script based on the initial geometric and process information. The script differentiates between standard heating time steps, tape switching steps with no heating, layer switching steps for cases when the laminate is examined after each layer, and the final long cooldown step before demoulding. Figure A3 shows the inner structure of the script, which is repeated for each time step.
In addition to the heat source, the laminate is also subject to thermal influences of the surrounding air and the tool. Also, the elements belonging to tapes which are not yet placed must not influence the process. This is achieved by deactivating the relevant elements, which is represented in Ansys by setting all thermal and mechanical properties as close to zero as possible. Heat conduction between active and inactive elements is eliminated. When a tape is due to be processed, its affected elements are reactivated and the material properties are reset. The current surface of the laminate is defined as all element surfaces with neighboring deactivated elements or without neighboring elements. On these element surfaces, a convective boundary condition using standard room temperature as reference is applied to simulate the influence of the surrounding air. The tool uses literature values for its thermal properties (e.g., of steel or aluminium). The interface to the laminate uses a constant thermal resistance, while the opposite surface is heated or cooled by a constant temperature boundary condition set to the tool temperature. During the final cooldown, this boundary condition is swapped to a convective boundary condition at room temperature. Figure A4 demonstrates the applied boundary conditions. The convection coefficients and heat transfer values were calibrated with thermocouple experiments. Between the panel and the tool a thermal contact with an experimentally-determined resistance of 750 W/m 2 K is defined. The internal heat transfer calculation within the panel uses Equation (1), which includes the density q and specific heat c p as well as the directional thermal conductivity K x , K y , and K z of the material for the transient temperature T calculation.
This part of the simulation calculates the thermal history of the part, which is transferred to a structural simulation afterward as an input parameter. The structural simulation uses the same script structure displayed in Figure A3 and the same mesh geometry to calculate the necessary time steps and the currently active and inactive elements. It also applies a boundary condition to fix the laminate on the tooling during the process and cooldown as well as insert a pressure boundary condition for the compaction roller. As a simplification, the compaction roller uses the same position calculation as the heat source, transferred one layer in positive z-direction. After the cooldown, an additional time step is inserted where the laminate is separated from the tool and can transform freely. Only the corners at the origin of the panel coordinate system and at two additional corners are fixed to inhibit free body motion. The thermal history influences the local degree of crystallinity of the material model and induces local stress, which transfer into warping afterwards. The detailed calculations are shown in Equation (2) through (12) below.
The geometry of the laminate used as an example was based on a standard laminate size used at the local T-AFP machine. The geometric parameters are shown in Table A1.
Within this study, carbon fiber reinforced tapes with low-melt polyaryletherketone matrix (CF/LM-PAEK) are used as the reference material owing to its popularity in contemporary T-AFP. Process parameters for the simulations were chosen based on the best process parameter set developed in experimental investigations within the context of Schiel et al. [16]. The tape properties and process parameters are listed in Table A2. For this investigation, a similar setup is used to process a laminate both on a heated and an unheated tool. The heat transfer coefficients have been determined as artificial values by thermocouple measurements.
The semi-crystalline nature of the matrix material has to be considered in the numeric material model. A material model for CF/PEEK was developed by Gordnian [17] and subsequently modified by Teltschik et al. [18] for CF/LM-PAEK. This material model iteratively calculates the relevant material parameters dependant on the current degree of crystallinity X and temperature T. The calculation for the degree of crystallinity of the thermoelastic CHILE material model is calculated by Gordnian and Poursartip [19] on the basis of Godovsky and Slonimsky [20], Mantell and Springer [21], and Maffezzoli et al. [22] with the rate of crystallization defined as where K(T) defines the temperature dependence of the crystallization rate and f(X) the crystallinity dependence. K(T) is defined as while f(X) is determined experimentally and is provided as a lookup table. E m /R and E g /R are model parameters for the temperature dependence of the crystallization rate and are also provided as lookup tables. T M and T g are the melting and glass transition temperature, respectively. All additional parts of the equation are fitting parameters. Teltschik et al. [18] describe how to use different experiment, especially differential scanning caliometry (DSC), to derive these parameters and lookup tables. The new degree of crystallinity X i after step i with the duration Dt i is defined based on the previous step n -1 as According to Godovsky and Slonimsky [20], the crystallization starts after reaching a defined induction time, which is defined in Equation (4). t m and c are additional fitting parameters.
ð t 1 0 dt t m T M À T ð Þ Àc ¼ 1 After each time step of the simulation, the temperature is compared with the relevant crystallization temperature range and the degree of crystallization is updated when necessary. In this way, the model incorporates both isothermal and non-isothermal crystallization. Subsequently further crystallinity-dependant material parameters are calculated. These include the matrix modulus E m and Poisson's ratio m according to Chapman et al. [23].
These parameters are calculated depending on the current temperature and its relation to T g and T M as well as the rubber-elastic transition temperature T 1 . the matrix moduli and corresponding Poisson's ration values for glass-like (G), rubberelastic (R), and liquid (Liq) states E m G , E m R , and E m Liq as well as m G , m R , and m Liq , respectively, are further refined by the temperature and the fitting values A 1 through A 6 .
The stress-strain-relationship is given as a differential equation by Gordnian [17] as a thermo-viscoelastic model based on the integral equation of Schapery [24] for creep of thermo-rheologically complex materials.
Here, r X and X refer to the stress and strain of the material dependent on the degree of crystallinity. The formula integrates terms for the relaxed modulus E e , the relaxed magnitude of stress due to temperature change DT in a completely constrained body b e ¼ aE e with the coefficient of thermal expansion a, the relaxed stress due to crystallinity change DX in a constrained body c e ¼ a CS E e with the free strain due to crystallinity changes a CS and the state variable q i defining the individual stresses of the elements for a generalized Maxwell model as a mechanical analogue for the viscoelastic material behavior.
For the final warping calculation, the volumetric shrinkage dVCS and linear shrinkage dLCS due to the degree of crystallization are calculated according to Chapman et al. [23] with the theoretical matrix densities of q amo for a pure amorphous matrix and q cry for a pure crystalline matrix resulting in a real matrix density q m . For the shrinkage, the densities of the current (i) and previous (i -1) time step are necessary.
In addition to shrinkage, the coefficient of thermal expansion is also included in the material model. The resulting mechanical matrix properties are used in the micromechanical equations according to Bogetti and Gillespie [25] combined with the temperature independent fiber material values and fiber volume content to calculate the combined unidirectional layer properties. The final laminate properties are the result of the classic laminate theory.
The material parameters for the LM-PAEK material relevant for the simulation are listed in the Appendix of this paper. The heat conduction coefficient perpendicular to the tape surface has been modified according to thermocouple measurements to reflect the reduced heat transfer due to incomplete consolidation between tapes.
The warping obtained by the simulations of the laminate on the heated and unheated tool after demoulding is compared to the measured warping of real panels manufactured with the same set of process parameters. The panels are measured after demoulding using an optical noncontact GOM Atos 5 3D scanning device, which uses two stereo cameras to take pictures of a projected pattern on the target to calculate the coordinates of each reflected point in relation to the sensor position. The location of the sensor itself is determined by special reference points adhered to the panel whose relative position to each other is recorded by the machine. The sensor cameras can use different lenses to either capture more details or measure a larger volume at once. For this analysis, lenses with medium focal lengths are used to capture objects of up to 500 mm with one scan while still capturing fine details. Multiple scans from different directions and angles are carried out until the whole part is captured.

Results
Reiterating from the previous section, two separate configurations were analysed: a quasi-isotropic laminate manufactured on a room temperature and a 200 C tooling. Both variants were simulated and their real counterparts measured. Figure A5 shows the temperature and degree of crystallinity while placing three different tapes measured in z-direction of the hot tool panel. A similar depiction for the cold tool panel is shown in Figure A6. Both figures include a diagram for the final simulated degree of crystallinity of each panel. Figures A7 and A8 show the warping of both the simulated panels and their material counterparts. The warping scale the colors is optimized for each image and varies: for the heated tooling in Figure  A7 the warping of the simulated panel on the left side ranges from À6.7 mm to þ5.1 mm while the measured warping has minima and maxima of À5.4 mm and þ5 mm, respectively. Table A3 lists specific warping values for the points specified in Figure A9. On average a deviation of 0.48 mm is observed between simulated and measured warping with a standard deviation of 1.36 mm. The same principle applies to the results from the cold tooling in Figure A8, where the simulation predicts warping between À13.7 mm and 0.5 mm while the scan indicates À19 mm to þ2 mm. Similar to Table A3 the difference between simulated and measured warping for this panel is listed in Table A4 which shows an average deviation of 3.55 mm with a standard deviation of 3.06 mm.

Discussion
The results of the previous section show a good match between the simulated and experimental results. In particular, the results from the heated tooling simulation are close to the measured warping of the real panel. The distribution of the warping shown in Figure A7 is very similar for both results. Table A3 also indicates a good match of the results with little variation.
Panels manufactured on cold tools behave differently. While the maxima and minima of the deformation of simulation and experiment are still in the same order of magnitude, the difference between both are larger than those encountered with the heated tooling. The warping distribution shown in Figure A8 is also different, although it is similar enough for quantitative predictions. Table  A4 shows that the largest deviation between simulated and measured warping lies in the top left area. This is to be expected since the top left corner is the only one not used to inhibit free body motion. But the high general deviation between the results indicate a systematic error in the simulation results.
Since the faster cooldown due to the cold tool poses the main differences between both panels, the error originates probably from the crystallization formulation of the material model. Since the lookup tables for the material were created with a maximum cooldown rate of 500 C/min during the DSC measurements, this may not be sufficient for an encompassing material model. Since Qureshi et al. [2] already showed the value of heated tools, T-AFP manufacturing on unheated tools is already unfavorable. For this case, refined material models are necessary.

Conclusion
During this study, a method was developed to model the T-AFP process and simulate the processinduced warping including the effects of crystallization. The simulated warping has been verified with measurements of physical panels manufactured by the same parameters as those used in the simulations. Two panel configurations were chosen for this comparison, with the difference of the tool temperature as the only different parameter. While the panel manufactured on a heated tool showed matching warping patterns for numeric and experimental results, the laminate from the unheated tool experiment showed larger differences. These can be traced back to the fast cooldown rate of this case, which is not yet covered by the material model.
Despite this limitation, the simulation can prove valuable for deriving optimal process parameters for parts manufactured by T-AFP on heated tools. These numeric optimizations can also include nonstandard variants, e.g., changes in the tape placement sequence, without the high expenses and efforts of physical experiments.

Disclosure statement
No potential conflict of interest was reported by the authors.           Figure A8. Simulated (top) and measured (bottom) warping of the cold tooling panel. Table A3. simulated and measured warping at specific points for heated tooling.