Heat transfer analysis of automated fiber placement of thermoplastic composites using a hot gas torch

Abstract With more and more use of composites in engineering applications, the need for automated composites manufacturing is evident. The use of automated fiber placement (AFP) machine for the manufacturing of thermoplastic composites is under rapid development. In this technique, a moving heat source (hot gas torch, laser, or heat lamp) is melting the thermoplastic composite tape and consolidation occurs in situ. Due to the rapid heating and cooling of the material, there are many issues to be addressed. First is the development of the temperature distribution in different directions which gives rise to temperature gradients. Second is the quality of the bond between different layers, and third is the rate of material deposition to satisfy industrial demand. This paper addresses the first issue. The temperature distribution affects the variation in crystallinity, and residual stresses throughout the structure as it is being built. The result is the distortion of the composite laminate even during the process. In order to address this problem, first the temperature distribution due to a moving heat source needs to be determined. From the temperature distribution, the development and distribution of crystallinity, residual stresses and deformation of the structure can then be determined. As the first phase of the work, this paper investigates the temperature distribution due to a moving heat source for thermoplastic composites, without considering the material deposition. A finite difference (FD) code based on energy balance approach is developed to predict the temperature distribution during the process. Unidirectional composite strips are manufactured using AFP and fast-response K-type thermocouples (response time of 0.08 s, as compared to normal thermocouples with response time of 0.5 s) are used to determine the thermal profiles in various locations through the thickness of the composite laminate subjected to a moving heat source. It is shown that temperature variations measured experimentally during the heating pass, using thermocouples embedded into the composite substrate, underneath layers of the composite material, are consistent with the generated thermal profiles from the numerical model. The temperature distribution, in both the direction of the tape and through-thickness direction can be predicted numerically.


Introduction
Over the past few years, the advent of automated fiber placement (AFP) has increasingly attracted the interests of industry due to faster rate of deposition, lower cost and repeatability of the process [1][2][3][4][5][6][7][8][9][10][11][12]. Thermoplastic composites offer many attractive characteristics such as high fracture toughness, better recyclability and good temperature resistance as compared to thermoset matrix composites [1,13,14]. However, the use of AFP machine for the automated manufacturing of thermoplastic composite structures with free edges has many technical challenges that hinder the fabrication of parts that are distortion free.
In AFP of thermoplastic composites as shown in Figure 1, there is a heat source (either laser, infrared, heat lamp or hot gas torch) which is used to heat up the incoming tape and a certain area of the substrate. The head assembly moves at a constant speed to provide heat while the compaction roller consolidates the incoming tape onto the substrate. As the head assembly completes one pass for depositing the first layer, it initiates the next pass to place the subsequent layer. Through a cyclic process, laying up is performed to achieve the required thickness [2,3,14,15]. Due to the rapid heating and cooling of the thermoplastic composite material, there are many issues to be addressed. First is the development of the temperature distribution in different directions which gives rise to temperature gradients. Second is the quality of the bond between different layers, and third is the rate of material deposition to satisfy industrial demand. This paper addresses the first issue. The temperature distribution affects the variation in crystallinity, and residual stresses throughout the structure as it is being built. In fact, residual stresses are developed as a result of temperature gradients because of different cooling rates in different layers along the thickness of the laminate. Consequently, the distortion of the composite laminate takes place even during the process. There are other reasons as well which gives rise to the development of residual stresses. These include accuracy of the tape material and the lay-up process, i.e., tape tension, tape feed, pressure control, and the nonuniformity of the quality of the bond between layers etc. The end result is that it may not even be possible to manufacture a flat unidirectional composite plate to make a test coupon. This is the inherent problem associated with thermoplastic laminates made by AFP discussed by several authors in literature [2,[16][17][18][19][20][21]. In order to address this problem, to make laminates with free edges with no or minimum distortion, one needs to understand the effect of different parameters on the process. As such, the first thing to look at would be to determine the temperature distribution history during the process. The focus of this paper, specifically, is to gain more understanding of the way temperature gradients develop during the AFP process of thermoplastic composites.
The governing heat transfer equation for orthotropic materials (e.g., composites) is presented as follows [14]: In this equation, q is the mass density (kg=m 3 ), C p is the specific heat (J=kg C), T is the temperature at a material point ( C), U ðJ=kg: SÞ is the heat generation term, and K x , K y , K z ðW=m CÞ are thermal conductivities in longitudinal, transversal and through-thickness directions, respectively.
In recent years, several studies have been devoted to develop heat transfer models for AFP process. Some of these have been focused on the development of temperature across the thickness of the laminate (z direction) at one particular location. The transient and steady state one-dimensional (1 D) problems were solved either analytically or numerically by various authors [22][23][24][25][26][27]. Although these studies can predict the overall temperature gradient through the thickness using a 1 D thermal model, due to several simplifications and approximations applied on the model, they cannot reflect the temperature distribution in the entire structure. As such, several studies focused on extending the 1 D model to a more sophisticated two-dimensional (2 D) thermal model by considering the temperature distribution along the lay-up direction (x-direction) and through the thickness (z-direction) as shown below.
Earlier studies by Beyeler et al. [28] and Nejhad et al. [29] in 1980s and 1990s addressed the 2 D heat transfer analysis of the tape laying process. In these studies, a mapping technique was used to transform the orthotropic domain into an isotropic domain to obtain a set of differential equations. Using finite difference (FD) techniques, the expressions are discretized and temperature profiles along both lengthwise and through-thickness directions were predicted under the steady state condition. Finite element (FE) techniques as an alternative numerical method were also implemented by Grove [30], Mantell and Springer [31] as well as Sonmez and Hahn [32] to investigate the principal thermal phenomenon taking place in both tape winding and tape placement processes. The model was later used for investigating the effect of process variables on the quality of thermoplastic composites to establish a processing window [20,33]. The previous studies were limited to steady-state condition by simplifying the transient problem. The identical FE approach was later extended to include the time dependency of the temperature (transient condition) in a 2D model. In 2011, Zhao et al. [6] constructed an ANSYS simulation model to perform a transient heat transfer analysis in the AFP process. A novel strategy based on the birth and death of elements technique was used to model the moving heat source and step-by-step growth of composite substrate (material deposition) in the lay-up process. The predicted thermal history was used in a thermoelastic FE model to compute thermal stresses induced during the process. The manufacturing technique they aimed at modeling, was processing composite rings using AFP in which structures do not have free edges and thus the problem of distortion does not arise. The identical birth and death of elements technique was later utilized by Li et al. [15] to obtain the temperature field during the processing of flat thermoplastic composite parts. The effect of heater temperature and roller speed on the bonding temperature was investigated. The numerical results were compared with analytical solution of 1 D heat conduction differential equation. Further investigations were performed by Han et al. [34] to study the effect of preheating on heat distribution inside the laminate. They stated that preheating, as an influential process parameter, can expand the melting area in both longitudinal and through thickness directions. They introduced boundary conditions; however, FE methodology implemented for the problem was not presented. Thermal conductivity values for the orthotropic domain were assumed to be constant and independent of the processing temperature. This might lead to erroneous results, since for the case of thermoplastic composites, the processing temperature is extremely high and there are dramatic changes in thermal properties at elevated temperatures which needs to be taken into account. Moreover, verification of their model was not performed and thus the accuracy of results was debatable. More recently in 2018, Kollmannsberger et al. [35] performed a 2 D transient simulation of laser-assisted AFP. In previous studies, perfect thermal contact between the plies of composite laminate was assumed; however, they took into account thermal contact resistance in the model. They compared temperature results for two cases of with/without the presence of thermal contact resistance and they highlighted that in the process modeling for an accurate analysis, it needs to be considered. The focus of their study was on investigating the effect of thermal contact resistance between layers.
In an attempt to provide a more precise prediction tool for composites produced by AFP, a small group of researchers, developed multidimensional thermal models to study spatial distribution of temperature in all x, y and z directions. In 1996, James and Black [36] developed a three-dimensional (3 D) thermal model to predict temperature profile in the filament winding process using an infrared energy source. The generated differential equations were solved numerically using a FD technique and results were compared with experimental temperature measurements using an infrared camera in conjunction with the winding apparatus. Unlike the previous model, Toso et al. [37] considered transient 3 D effects of the tape winding process. The numerical simulation of the process was constructed based on FE in ANSYS. The predicted temperature results during the lay-up process were compared with experimental measurements using infrared pyrometry to provide validation. The tape winding process was the focus of their study in which structures are made in the form of cylinders or spheres having no free edges. As such, the issue of distortion does not arise, and further studies are required to make structures in the form of strips with free edges where distortion plays an active role. Hassan et al. [38] considered the same analysis strategy as [37], yet for thermoset materials. They discretized the model into eight-node brick elements to obtain coupled ordinary differential equations. A 3 D FE code was written to obtain nodal temperatures during the lay-up process of a composite ring. In order to verify the FE code, carbon/epoxy composite rings were fabricated using a filament winding machine and K-type thermocouples were used to measure the temperature during the process. Although the 3 D transient analysis of the problem was introduced, the model was limited to manufacturing of thermoset composites using the filament winding technique. A few years later in 2016, Jeyakodi [17] performed a FE simulation of the in-situ AFP process for manufacturing flat thermoplastic composites. The coupled temperature displacement analysis was performed in ABAQUS to be able to apply the transient laser heat and roller pressure in simulation. The main focus of his study was on predicting the residual stresses developed during the in-situ process. To do so, thermal profiles induced in the laminate for different process configurations were determined and induced stresses in the laminate were calculated accordingly. His study provides a detailed numerical model of laser-assisted AFP process. Although the simulation tool could be used for studying different process parameters and identifying the interactions numerically, experiments have not been performed to provide validation. This necessitates more experimental research on in-situ processing of thermoplastic composites with free edges to gain more understanding of the process and to fill the gaps in existing literature. More recently, Tafreshi et al. [14] constructed a FE model based on ANSYS to determine the spatial temperature distribution in thermoplastic composites made by AFP using hot gas torch. In their study, a few different approaches were taken, including closed form (analytical) and FE, in order to obtain understanding of the heat transfer problem. The FE model was developed for isotropic case and then extended to the orthotropic case in which experiments were performed to provide validation of the results. Each heat transfer model associated with the process in literature has its own particular application, and thus, thermal models and corresponding process parameters change for different composite materials (thermoset or thermoplastic), specifications (composite rings or flat strips), and heating mechanism (laser, infrared, heat lamp or hot gas torch). Looking at existing literature, it is evident that there are number of works addressing processing of thermoplastic composites along with experimental validation of the results. Nevertheless, they are limited to either-filament/tape winding process or laser-assisted AFP. The focus of this paper is on processing thermoplastic composite parts with free edges using hot gas-assisted-AFP. There are differences between the previous manufacturing processes: both filament winding process and fiber placement process are considered additive in nature because the material (thermoset/thermoplastic) is deposited layer by layer. Nevertheless, there are some features that distinguish the fiber placement process. Firstly, filament winding is applicable to make structures in the form of cylinders or spheres (mainly convex surfaces) while fiber placement technique can be used for structures with both convex and concave surfaces. Secondly, filament winding can only be used to make structures having no free edges, where the problem of distortion does not arise. Fiber placement can be used to make structures with free edges, where the issue of distortion is important. Thirdly, in filament winding it is required to apply tension on the tows, while in fiber placement process, tows are placed on the surface of the tool using a consolidation roller [1,21,39]. Hot gas torch, on the other hand, has been implemented since 1986 as the primary heat source for processing thermoplastic composites using AFP due to its low cost and wide process window [40]. More recently, laser heating has garnered attention among manufacturers as it provides high energy density, faster processing rates and a better surface finish [40]. Even though this approach may be appropriate for certain applications and people tend to think that laser is better than hot gas torch (due to the fact that it is more focused, and heating can be more efficient), there are benefits from the hot gas torch that laser does not have. One thing is that laser cannot be used for glass fiber composites. This is because laser heats the carbon fibers and the heat from the carbon fiber heats the resin. Glass fibers do not absorb the laser energy and as such it cannot be heated by laser. Second, laser imposes safety hazard. Third, hot gas torch spreads out the heat and this can help to preheat the tape. This is more effective for heating for joints. And fourthly, laser needs to be precisely controlled for the end of the beam to be at the NIP point while the hot gas torch has less of a problem.
In summary, although the models proposed in literature provide insight into the transient heat transfer problem of the AFP process, they are either limited to the local region in the vicinity of the nip point and did not take into account the temperature distribution in the entire structure; or owing to several simplifications in the model and lack of experimental data, the predicted temperature history did not reflect the actual AFP process where the problem of distortion arises in manufacturing of composite parts with free edges. The focus of this paper is to develop a two-dimensional transient heat transfer model for moving heat source problem associated with AFP thermoplastic composites manufacturing with experimental validation of results to gain more understanding of the way temperature gradients develop during the process.
Heat transfer mechanism in AFP thermoplastic composites process As the first step in heat transfer analysis of AFP process using a hot gas torch, it is necessary to understand the heat transfer mechanism between the hot gas torch and the composite substrate. In fact, the heat transfer mechanism is quite complex not only because of the turbulent nature of the thermal fluid flow, but also due to the simultaneous heat transfer that occurs in both the thermal fluid and the composite substrate. The hot gas flow coming out of the torch nozzle, has very high temperature in the range of 800-1000 C for the case of AFP of thermoplastics. As such, noticeable temperature difference exists between the hot gas flow and the composite surface. This results in convective heat transfer between the aforementioned media as shown in Figure 2. A local heat flux equation is utilized to model the convective heat transfer as q " ¼ h HGT ðDTÞ where h HGT is the convective heat transfer coefficient between the hot gas torch and the composite surface, and for the conversion of a portion of the kinetic energy of the gas upon impact with the surface of the composite material. DT accounts for the temperature difference between the two media [41,42]. The localized heat flux diffuses through the composite substrate ( Figure 2) and the speed of diffusion depends upon the value of thermal diffusivity a; the larger the value of a, the faster heat will diffuse through the material. In fact, a temperature gradient exists in the laminate which results in an energy transfer from the upper surface subjected to the hot gas flow (high-temperature region) to the underneath layer (low-temperature region) [41].
In this paper, a moving heat source is considered in the transient heat transfer problem. Along with the heat transfer through the thermal fluid (heat convection) and through the composite substrate (heat diffusion), the unheated regions are exposed to the surrounding ambient which plays an important role in the cooling stage. The convection heat transfer follows the Newton cooling law as convective heat transfer coefficient between the composite substrate and the surrounding ambient, and DT is the temperature difference between the aforementioned media [41][42][43]. The composite substrate surface is subdivided into two separate zones as shown in Figure 2. Zone 1 represents the region of the substrate exposed to the surrounding ambient and cooled by ambient air, while zone 2 demonstrates the region of the composite substrate which is exposed to the hot gas flow. Therefore, at any moment during the process, the convective heat transfer coefficients can be expressed as follows: The convective heat transfer coefficient h HGT may vary along the length of the zone 2. As a preliminary analysis, it is taken to be a constant in this study.

Transient heat transfer model of the AFP process based on finite difference formulation in MATLAB
The 2D rectangular domain and the coordinate system considered in this paper are illustrated in Figure 3. There is a moving heat source which travels with a specified speed from the left edge to the right along the upper side of the rectangular domain. The heat energy provided by the heat source, is transferred to the composite substrate which results in nodal temperature increase. At the same time, the upper side is subjected to the surrounding ambient (room temperature), which accounts for the cooling stage. In this section, we aim at developing the FD formulation of the transient heat transfer problem in a 2D rectangular domain (the composite substrate-aluminum tool assembly) subjected to a moving heat source using the energy balance approach. First, the domain is subdivided into a sufficient number of volume elements and the energy balance is applied on each element. An appropriate formulation is derived for the interior nodes. The formulation can be used for all the interior nodes regardless of the boundary conditions, since the boundary conditions have no effect on the FD formulation of these nodes. Second, for nodes on the boundaries, depending on the type of the boundary condition, an appropriate formulation is derived for each case. The explicit methodology is introduced as a solving method for the transient heat transfer problem. The stability criteria are defined to confine the variation and oscillation of the nodal temperatures. Eventually, the formulations are coded to be able to generate the nodal temperatures after each time step and plot the temperature distribution.
Assumptions Prerequisite for a valid FD model is a feasible simplification of the problem, in order to be able to perform the analysis with short computational time with mathematical tools available. The FD model proposed in this paper is developed based on the following assumptions: a. In the actual AFP process, the material deposition is part of the model. In other words, it is necessary to address the placement of the incoming tape after passing through the AFP guided assembly on the substrate with the help of compaction roller. However, in this paper, it is assumed that there is no material deposition and the model is developed for the case of the moving heat source on a composite substrate being in contact with an aluminum tool. b. The two-dimensional heat transfer problem is considered in this paper. There are some reasons highlighting this assumption: the composite substrate is relatively long with respect to the width and one can neglect the temperature variation through the width. Second, the heat input is almost constant over the tape width and the material is uniformly heated in that direction. As such, one can reasonably ignore the temperature variation in width direction (y-direction) and the assumption is rated tolerable [29,31,32,44]. c. The process is regarded as a transient heat transfer problem which is necessary in order to investigate the temperature profiles in different regions of the composite substrate as time progresses. d. The aluminum tool is considered in the model as it cannot be ignored in the actual manufacturing process. e. The composite substrate is modeled as an orthotropic domain and the thermal conductivity depends on the direction (longitudinal and transversal). However, for the aluminum tool, only one value is assigned as it is isotropic. f. The thermal conductivity values for the composite domain are temperature dependent. Because in higher temperatures, the molecules move faster, and energy will be transported faster. As such, the value of thermal conductivity, as an indicator of how fast heat flows in a material, needs to be temperature dependent [41]. g. The hot gas flow exposing on the composite substrate surface is modeled as a convective boundary condition with a constant heat convection coefficient. h. A convective boundary condition is prescribed to the composite substrate surfaces (exposed to hot gas flow) and left, right and bottom boundaries (exposed to ambient air). i. The heat absorption and generation resulting from melting and crystallization respectively is assumed to be negligible.

Energy balance approach
Consider a transient 2D heat conduction problem in a rectangular domain. The length of the domain, the thickness of composite substrate and the thickness of tool are subdivided into M, N c and N m sections respectively. In the case of the transient condition, in addition to the discretization in space, one needs to consider the discretization in time as well. At time t ¼ 0, initial condition is applied, i.e., it is assumed that before starting the application of the heat source, the temperature is maintained constant throughout the entire domain and it is the same as room temperature. After one time step Dt, the nodal temperatures need to be updated based on the previous step. Solving for the unknown nodal temperatures repeatedly for each Dt, the final nodal temperatures at desirable time can be calculated. The selection of time step will be discussed in the following section. In this paper, the nodal temperatures are denoted using double subscript notation (m, n) along with the index i as T i m, n which means the temperature at a node (m, n) at time step i [42]. If one applies the energy balance for a volume element during a time interval Dt, the nodal temperatures at time step i þ 1, can be calculated as shown in Table 1.
The explicit method used in this paper to formulate the transient heat transfer problem is not unconditionally stable, and results may fluctuate and become unstable if some conditions are not met. Stability criterion indicates that the explicit solution is stable if the time step Dt is sufficiently small. The value of Dt must be maintained below a certain value throughout the solution to avoid instability in nodal temperatures. On the basis of thermodynamics, the stability criterion is satisfied if the coefficients of all T i m, n in the explicit representation of T iþ1 m, n are greater than or equal to zero for all nodes. Based on this definition, one can obtain the required step time to have a stable solution [42].

Input values directory
At the first stage, the input values are entered by the user. These include physical and thermal properties of the material, dimensions of domain, number of divisions in x and z directions and the attributes of the boundary conditions. For the case of this study, AS4/APC-2 carbon fiber/PEEK Table 1. Nodal temperatures for interior and boundary nodes based on energy balance approach.
Interior nodes T iþ1 m, n ¼ Dt Dx : mesh size in the x direction, Dz: mesh size in the z direction, h C : convective coefficient of zone 1, T 1 : surrounding temperature, h HGT : convective coefficient of zone 2, T HGT : HGT temperature.
composite material is selected for the composite substrate and aluminum is considered as the tool's material. The physical and thermal properties for AS4/APC-2 were reported in the literature [17] and those values are used in this paper as input. The corresponding properties of the aluminum tool were reported based on [6]. Experimental values and boundary condition attributes are also required to proceed the model. Temperature of the hot gas torch and the nitrogen gas flow rate is selected to be 875 C and 60 standard liter per minute respectively. The convective heat transfer coefficient for free convection (h c ) is selected to be 10 W=m 2 K [14]. The heating length, where the h HGT is applied, is selected to be 10 mm based on the experimental setup. Different values of the h HGT were reported by several authors in the literature [22,37,38,44]. In these studies, an initial value of h HGT was assumed to start the calculations. A few thermocouples were utilized to obtain the temperature data experimentally. A comparison between experimental and numerical results was done; a few iterations were carried out until a good match is found to obtain the h HGT : In this paper, following a procedure similar to what was presented in the literature, the maximum temperature corresponding to the position of TC3, embedded one layer beneath the top surface, obtained experimentally. It was then compared with the maximum temperature predicted numerically at that location. Following this procedure, the value of h HGT of 990 W=m 2 K was obtained. In order to ensure that the selection of this value can lead to numerical results with reasonable accuracy, the temperatures at two other through-thickness locations were also monitored using two other embedded thermocouples (TC1 and TC2) and results were compared with numerical predictions. A summary of all input parameters for the FD code is presented in the appendix.

Procedure for the analysis of heat transfer by a moving heat source
In the actual AFP process, the heat source moves continuously. However, in the FD simulation, the continuous movement of the heat source needs to be discretized into specified steps. Assume that the sample has a total length of L ¼ 20 inch (508 mm). The speed of travel of the torch (moving heat source) is v ¼ 1 inch/sec or 25.4 mm/sec. The total time for the complete process is L/v ¼ 20 s. The length of the heating region is 10 mm. Concerning the fact that the length of the sample is subdivided into 508 nodes (the distance between nodes is 1 mm), the heating region encompasses 10 successive nodes as shown schematically with red circles in Figure 4. The way the moving heat source is modeled in this paper is based on the fact that the heat source moves one node at a time. This means if at step 1 the center of the heating length is at x (t 1 ) ¼ 5 mm, at step 2, the center of the heating length would be at x (t 2 ) ¼ 6 mm (since the distance between nodes is 1 mm). As such, the time which takes for the moving heat source to move from step 1 to step 2 is equal to the distance between nodes (i.e., 1 mm) divided by the speed (i.e., 25.4 mm/s), so the time interval is 1/25.4 ¼ 0.04 s. In other words, the heat source moves every 0.04 s from one node to the adjacent node to simulate the speed of 25.4 mm/s of the heat source.
Based on the above analysis, the procedure for the analysis of the moving heat source for L ¼ 20 inch (508 mm), speed v ¼ 1 in/sec (25.4 mm/ sec) and heating length 10 mm is as follows: Step 1: During the first time increment (first stop), the center of the heating length is at x (t 1 ) ¼ 5 mm. Duration of stop ¼ 0.04 seconds. The time step at which the temperature of each node is re-calculated is equal to 0.0004 seconds which is determined by stability criterion for the problem. This means that at each stop, there are 100 time increments for the heat transfer analysis. Subsequently, the heat source moves one node (1 mm) and the calculation continues. The temperatures at different points in the structure at the end of this stop period are used as input to the next heating analysis in step 2 below.
Step 2: Heating length travels so that the new position of the center of the heating length is at x (t 2 ) ¼ 6 mm. Duration of stop ¼ 0.04 seconds. Using the temperature from the previous step (step 1) as initial conditions, the finite difference analysis is carried out. The temperature at the end of this step is used as input for the initial conditions for the analysis in the next step.
The same procedure is carried out until as far as necessary. This procedure is shown schematically for four first successive steps in Figure 4.

Mesh size study (mesh convergence)
In this work, a continuous process is converted into a discrete process, we do not know exactly which mesh size would be sufficient to obtain the numerical results from the heat transfer analysis. Using a sufficiently refined mesh is important to ensure that numerical results do not deviate from the actual results. Coarse mesh may yield erroneous results and needs to be refined. This will guarantee that numerical results are independent of the mesh size and results are valid. Table 2 presents different meshes used for the mesh size study (mesh convergence), including the number of divisions and total number of nodes for the heating length for each mesh size, the jump distance (the distance the heat source is moved each time), and the stop time (the time which the heat source rests at each stop).
In order to see how the mesh size would affect the numerical results, temperature distribution results are presented in Figure 5 for four values of number of divisions for the heating length, i.e., 10, 8, 6 and 4 (corresponding to the mesh sizes of 1, 1.25, 1.66, and 2.5 mm). Assume that we have a total length of 20 inches (508 mm) and the heat source (heating length of 10 mm) moves with the speed of 1 inch/sec (25.4 mm/s) from the left to the right of the sample. Temperature variation for the node located in the middle of the sample (x ¼ 10 inch), underneath 9 th layer is presented in Figure 5. The horizontal axis is limited from 8000 to 14,000 ms to be able to distinguish the difference between results from different numbers of divisions. Looking at Figure 5, it can be noticed that the difference between the temperatures from different mesh sizes varies with time. The difference is largest when the maximum temperature occurs (at around 11,000 ms) and decreases as the time moves on. As such, if one waits long enough for the structure to cool down, the mesh size has smaller effect.
A summary of numerical results for all mesh sizes considered in this study is provided in Figure  6. As can be seen from Figure 6 (Left), as the number of divisions for the heating length increases from 1 to 10 (corresponding to variation of mesh size from 10 to 1 mm), the maximum temperature results converge. Furthermore, in order to see how fast or how slow is the rise from room temperature to the maximum temperature at each layer, times to reach maximum temperature values are plotted for different mesh sizes in Figure 6 (Right). A similar trend can be observed, and results converge for the mesh size of 1 mm (number of divisions for the heating length ¼ 10). Overall, it can be noticed that as the mesh size decreases and reaches 1 mm, the difference between the results is insignificant and it can be concluded that the number of divisions for the heating length of 10 (corresponding to mesh size of 1 mm) is satisfactory to be used to generate numerical results in subsequent sections in this paper.

Experiments
AFP workcell and flat paddle tool at Concordia Centre for Composites (CONCOM) developed by Automated Dynamics Corporation (ADC) were used to perform the experimental studies in this paper. The AFP system employs a nitrogen hot gas torch to melt the incoming tape and a steel roller for compaction. Fully impregnated 1 = 4 -inch wide AS4/APC-2 carbon fiber/PEEK tapes supplied by CYTEC Engineered Materials were used for the layup trials. The unidirectional tape consists of a 68:32 weight percentage mixture of carbon fiber (AS4) and Polyetheretherketone (PEEK) matrix (APC-2). The nominal value of fiber volume fraction is 61% [12,45]. Figure 7 depicts the experimental trials setup used in this paper. Unidirectional composite strips of 20 inch long making up 20 layers were manufactured to investigate the temperature distribution in various locations through the thickness of the composite laminate when the moving heat source travels across the upper side of the laminate. To do so, 8 unidirectional composite strips of 20 inch long making up 8 layers were first manufactured, and the laminate was allowed to cool down to room temperature. Afterwards, on the upper side layer, a fast-response K-type thermocouple (TC; Omega CHAL-002; response time 0.08 s) was mounted on top of the eighth composite strip using a Kapton V R tape. Subsequently, the ninth layer was deposited. When the NIP point of the roller arrives on top of the TC, it deposits the composite layer and the TC became embedded into the composite substrate. The TC data were measured by Data Acquisition system (DAQ) at 100 Hz. It is important to have a TC with fast response time.
Otherwise, the temperature reading is not accurate. The aforementioned TC is denoted by TC1 in Figure 8. Then, 3 more layers were deposited on top of the previously laid composite substrate along with another TC (denoted by TC2 in Figure 8). The system was allowed to cool down to reach room temperature. Another 8 layers were deposited afterwards on top of TC2 to form 19 layers of composite laminate with two embedded TCs. Eventually, the third TC was mounted, and the last layer was deposited, and the entire system was allowed to cool down for 24 hrs. The entire layup consists of 20 layers of unidirectional composites with three embedded TCs of the same type to capture the temperature variation during the experimental trials. In order to compare the simulation results from moving heat source with experimental results, the material deposition was deactivated in AFP setup and the nozzle temperature was set at 875 C and was kept constant during the experimental trials. The distance between the torch beginning position and TC was 10 inches. The heat source moved from one end to the other end with speed of 1 in/s without material deposition and the temperature was recorded by TCs which were embedded into the composite substrate, underneath layers of composite material. After completing one pass, the setup was allowed to cool down to room temperature. The process was repeated five times to ensure readings are repeatable, and the data were consistent.

Results and discussion
In this section, for the case of a moving heat source travelling across the upper side of the composite substrate, the nodal temperature results were generated from the developed code based on FD formulation (node to node distance of 1 mm). For the selected nodes, experimental trials were performed to validate numerical results. The nodal temperature results are presented for different locations through the thickness of the composite substrate, all located on the mid-line at x ¼ 10 in (254 mm). Nodes are selected in such a way that each node is located at the interface of two adjacent layers. As such, each node will be representative of the corresponding interface temperature in the composite substrate.

Numerical and experimental through-thickness temperature distribution results
Temperature distribution results are shown in Figure 9. Family of curves with solid lines represent numerical results while the dashed lines describe experimental results. Looking at the graph, it can be noticed that the transient temperature profile can be divided into 3 different stages. As the heat source moves from left towards the embedded TCs, the temperature maintains room temperature (25 C) for the first few seconds (Stage 1). This is because of    the fact that the heat source is at the beginning of the pass and the upstream region is not heated up yet. When the heat source approaches the location of the thermocouple, the temperature starts to rise up and reaches a maximum temperature (Stage 2). Once the heat source (the center of the heating length) passes the spot, that spot begins to cool down (Stage 3). Although the temperature of the gas is about 875 C, the maximum temperature of TC3 is only 277 C. This is because TC3 is embedded below one layer of the composite, and this layer of composite shields TC3 from the high temperature. It was shown in [14] that the temperature of a thermocouple placed on top of the laminate reaches to around 800 C when the hot gas torch passes over it. The speed of the torch (moving heat source) is 1 inch/sec (25.4 mm/s); as such it would take 10 s for the heat source (center of the heating length) to reach the embedded TCs. Looking at the results for TC3, it can be noticed that the maximum temperature occurs at around 10.3 s from the numerical results and 10.5 s from the experimental results. The time difference between points "B" (the point on the graph corresponding to 10 s) and "C" (the point where the temperature reaches maximum) is denoted as delay time. Note that at 10 s, the center of the heating length passes over the position x ¼ 10 inch (254 mm). The remaining five nodes of the heat length would take 0.2 s (5 mm/ 25 mm/s) to pass away from the position of x ¼ 10 inch (254 mm). As such out of the delay of 0.3 s, the non-heating time is only 0.1 s. The heat energy provided by the heat source needs to be transferred and absorbed by the material and thus the maximum temperature occurs with a time delay. The experimental data indicate the effect of other influential parameters which result in a longer time delay (0.5 s) as compared to the numerical results. This is because of the fact that TC has a response time and cannot measure the temperature instantaneously. As such, due to the response time of TC, there is a slight difference between the time delay prediction based on numerical results and experiments. Another interesting point is the fact that temperature starts to rise up at around 6 s (for TC3), i.e., before the heat source reaches the embedded TC. This phenomenon is denoted as rise time which is the time difference between points "A" (where temperature deviates from the horizontal line) and "B" as shown in Figure 9. The rise time is attributed to the heat transfer through the material along the direction of travel of the moving heat source [14].
In Figure 9, the numerical method is able to capture the maximum temperature with good agreement with experimental results. However, there is significant difference in the rate at which the temperature rises. The numerical results show much sharper rise in temperature from the ambient as compared to experimental results. This can be due to the way how the h HGT is applied. In Figure 2, the top surface is divided into region 1 and region 2.
The boundary between these regions is sharp in the numerical method whereas in reality, it may not be that sharp. Also, h HGT was assumed to be constant over the heating length but it may not be. These can be the reasons for the significant difference in the rate of cooling between the numerical results and experimental results. This is the subject of a further investigation.
Looking at temperature variations through the thickness (Figure 9), it can be noticed that as one goes through the thickness, the maximum temperature becomes lower and lower. For instance, the maximum temperature recorded by TC2 was 124 C, i.e., 55% lower temperature than TC3. The corresponding nodal temperature predicted by numerical code was around 112 C (9.6% error). Furthermore, TC1 would reach the peak temperature of 93 C which is almost 66% lower temperature compared with TC3. Overall, the temperatures measured experimentally during the heating pass, using embedded TCs into the composite substrate, underneath layers of the composite material, show consistent trends with the generated temperatures from the numerical model.
Using the developed numerical model, temperature distribution for the structure was determined. Figure 10 represents the temperature distribution in the vicinity of the heat source along the lay-up direction at different through-thickness locations. Zeta (n ¼ x À vt) is a transformation on the x-axis indicating the distance between the heat source and the target point in the lay-up direction (heat source moves with the speed of v and its mid length has the same x position as the target point when n ¼ 0). The negative values (n < 0) indicate the upstream region which was heated up by the heat source while the positive values (n > 0) demonstrate the downstream region where the mid length of the heat source has not reached yet [10]. As it can be seen from Figure 10, the temperature gradient in downstream region is sharper than the one in upstream region. The reason is because the upstream region was affected by the heat source before the center of the heat arrives at x ¼ 10 inch, while the downstream region was not.
Looking at Figure 10, it can be noticed that as one goes through the thickness, the maximum temperature becomes lower and lower until after five layers, the temperature is not above glass transition temperature T g (143 C) anymore. For regions where the temperature is below T g , the material becomes stiff and this contributes to the development of residual stresses and distortion due to temperature gradients. This is also of practical importance for repass treatment (reconsolidation pass) which has been reported in the literature [12,27] to improve bonding quality, reduce void content and improve surface finish of in-situ consolidated thermoplastic composites. Repass refers to Figure 11. Temperature variations through the thickness (z direction) of the composite substrate at five different x locations at t ¼ 10 s. the application of heat and pressure using the AFP head to the laminate without addition of a new layer. It can be seen from Figure 10 that repass treatment can be effective in improving bonding interface between layers where the temperature is above glass transition temperature T g (i.e., up to 5 layers below the top surface). Another interesting point, which is obtained from the observation of Figure 10, is the fact that the maximum temperature occurs at later time as one goes through the thickness. This is due to the heat diffusion mechanism inside the composite domain. A localized heat flux resulting from the hot gas flow diffuses through the thickness of the composite substrate. As such, it takes time for the heat energy provided by the hot gas flow at the upper side layer to be diffused, transferred and absorbed by the material at underneath layers. This can be interpreted as the delay time in peak temperatures at underneath layers.
Further investigation of through-thickness temperature distribution is illustrated in Figure 11. When the mid length of the heat source is located at x(t n ) ¼ 10 inch (at t ¼ 10 s), the variations of temperature through the thickness (z direction) from 0 mm (top of the composite laminate) to 2.5 mm (bottom of the composite laminate) at five different x locations are presented in Figure 11. For x(t n ) ¼ 10 inch (254 mm), this location of the composite laminate experiences the maximum temperature because it is directly below the mid-length of the heat source. This temperature variation demonstrates a sharp temperature drop until 1 mm (8 layers below) through the thickness where the temperature reaches below 100 C followed by a gradual temperature reduction until it stabilizes. The curve for x ¼ 9.8 inch (249 mm), shows the temperature variation in the upstream region which was previously heated up by the heat source and it is exposed to the surrounding ambience in the cooling stage. As can be seen, the temperature curve for x ¼ 9.8 inch initially shows lower temperature than the curve at x(t n ) ¼ 10 inch. However, as z increases, there is a cross over ("A") and within the lower half of the thickness, the temperature at x ¼ 9.8 inch is higher than the temperature at x(t n ) ¼ 10 inch. A similar trend can be observed from the curve at x ¼ 9.6 inch (plotted at x(t n ) ¼ 10 inch). The temperature is initially lower than the temperature at both x(t n ) ¼ 10 inch and x ¼ 9.8 inch; however, as z increases, there is a cross over with the curve at x(t n ) ¼ 10 inch at z ¼ 0.75 mm ("B") and another cross over with the curve at x ¼ 9.8 inch at z ¼ 1 mm ("C"). The curve at x ¼ 10.2 inch (259 mm) represents the downstream region as elaborated in Figure 10. Looking at this temperature distribution, it can be noticed that the temperature at x ¼ 10.2 inch is lower than the temperature at x(t n ) ¼ 10 inch. The curve at x ¼ 10.4 inch is also plotted to confirm the trend for the curve at x ¼ 10.2 inch.
The reasons for these are as follows: Considering the surrounding of two positions, x ¼ 9.8 inch and x ¼ 10 inch. At time t ¼ 10 s, both the left side and right side of position x ¼ 9.8 inch are at high temperature because they were previously heated, particularly for positions of higher z values. However, at time t ¼ 10 s, the left side of x ¼ 10 inch is hot while its right side has not been heated previously. As such, while for positions at lower z values, the temperatures are high due to immediate exposure to heating, positions at higher z values have lower temperatures since there was no previous heating to allow the heat to soak in. This same reason applies to all other positions.

Conclusion
This paper presented results of an investigation on the temperature distribution due to a moving heat source in AFP of thermoplastic composites. A FD code was generated in MATLAB based on energy balance approach to predict the temperature gradient in both lay-up and through-thickness directions. Procedure for the analysis of the heat transfer by a moving heat source was presented. In order to see the effect of the duration of each stop of the heating length on temperature distribution results, the duration of each stop was varied with the corresponding change in the space increment along the x direction. The analysis was performed for different number of mesh sizes and it was shown that for the mesh size of 1 mm, results would converge.
Experiments were performed to obtain the thermal profiles in various locations through the thickness of the composite laminate. It was shown that numerical method was able to capture the maximum temperature with good agreement with experimental results. The general trend for the temperature distribution also shows agreement. However, there is significant difference in the rate of rise of the temperature. This is attributed to the way how the boundaries between region exposed to hot gas and region exposed to air temperature are interfaced. It was shown that as one goes through the thickness, the maximum temperature becomes lower and lower until after five layers where the temperature is below glass transition temperature for the setup and process parameters considered in this study. Furthermore, temperature gradients in the vicinity of the heat source were obtained. It was shown that the downstream region experienced sharper temperature gradient than the upstream region. Temperature variation through the thickness showed a sharp temperature drop in the first few layers where the temperature reached below 100 C followed by a gradual temperature reduction until it is stabilized.
The paper provides significant insight into the distribution of the temperature due to a moving heat source, along with reasons for its behavior. The results show significant temperature gradients which explain why there are so much distortion even during the manufacturing process. This study can open the way for understanding the development and distribution of crystallinity as a function of space, time and temperature in addition to determination of residual stresses and deformation of the structure. Finite difference size and mesh configuration Composite Analyzed Length in the x direction (L C ) 20 inch (508 mmÞ Analyzed Thickness in the z direction (W C ) 1 mm (equivalent to 8-layer composite) Number of nodes in the x direction (M C ) 508 (mesh size ¼ 1 mm) Number of nodes in the z direction (N C ) 3 nodes per layer Tool Analyzed Length in the x direction (L m ) 20 inch (508 mmÞ Analyzed Thickness in the z direction (W m ) 2 inch (50.8 mmÞ Number of nodes in the x direction (M m ) 508 Number of nodes in the z direction (N m ) 2 0