Modelling the Heating Process in the Transient and Steady State of an In Situ Tape-Laying Machine Head

: As the use of composite materials increases, the search for suitable automated processes gains relevance for guaranteeing production quality by ensuring the uniformity of the process, minimizing the amount of scrap generated, and reducing the time and energy consumption. Limitations on production by traditional means such as hand lay-up, vacuum bagging, and in-autoclave methods tend not to be as efﬁcient when the size and shape complexity of the part being produced increases, motivating the search for alternative processes such as automated tape laying (ATL). This work aims to describe the process of modelling and simulating a composite ATL with in situ consolidation by characterizing the machine elements and using the ﬁnite differences method in conjunction with energy balances in order to create a digital twin of the process for further control design. The modelling approach implemented is able to follow the process dynamics when changes are made to the heating element and to predict the composite material temperature response, making it suitable for use as a digital twin of a production process using an ATL machine.


Introduction and Related Works
The search for more efficient and automated manufacturing processes for composite materials has found that the automated tape-laying (ATL) process with in situ consolidation for thermoplastic laminates is a good alternative for reducing material scrap and increasing the manufacturability of complex geometries by out-of-autoclave processes [1][2][3][4].
This process requires a machine mainly made up of a heating source and a compaction mechanism in order to raise the composite temperature, guaranteeing the quality of the final part by ensuring that all the laid composite layers have been welded among themselves [5]. The main parameters that affect the quality of the final product are the temperature and the compaction roll pressure at the nip point [3,6]. The temperature is the most critical parameter that must be controlled in order to reduce defects such as voids or delamination.
The heating process is the most critical stage. Reaching the necessary composite temperature ensures a good bonding condition between the layers [7,8]. The exposure of the thermoplastic matrix to inadequate temperatures may cause the degradation of the material and residual thermal stresses [2,[9][10][11], making the final part defective and unable to meet the required performance.
Most modelling works focussing on the use of a heating head for the tape-laying process that can be found in the literature use lasers as heat sources [12][13][14] due to their controlled power delivery at a fixed wavelength and localized region. The position of the heat source relative to the composite tape is one of the most relevant parameters, as presented by [15], because the heating process depends on optical properties such as reflectivity for the head assembly, the wavelength absortibity of the heated material, and surface irregularities.
Other studies using similar heat sources study the temperature distribution by measuring the temperature just before the nip point using thermal cameras [16] and analysing the results for three different power values of the heat source. The authors of [4] propose a model that uses lasers to heat a set of thinner composite tapes in order to study the influence of the width of the tapes on the overall heating process when the machine has to follow a non straight path.
Another heat source used for this process is hot gas torches [14,17], where the heat transfer into the composite mainly depends on the convection coefficient computed using empirical correlations.
A study that is not related to the manufacturing of composite materials and that discusses the topic of heat sources was carried out by [18,19]. This work presents a strategy to model a radiation heat source using infrared lamps for a continuous paper drying process. This work suggests an energy balance between the process parts involved, taking into account their properties.
Independently of the heat source, all models share the same phenomenological problem-that is, understanding the temperature response of the composite material during the process. The majority of authors in this field have focussed on numerical and experimental studies of the problem, assuming constant process parameters such as heat power and feed velocity without considering the use of control strategies to improve the overall process by adjusting its parameters.
Another common characteristic of the models proposed in the literature is related to the mathematical approach used. The use of a three-dimensional heat exchange was proposed by [20]; this implies a high computational cost that increases when the model results are used in combination with resin cure models, as proposed by [21]. A simplified model was proposed by [22]; this model uses a hot air gun as a heating source, in which a 1.5 D thermal approach is used; this model does not take into account variations in the material properties with temperature. On the other hand, an analytical model was developed by [23]; this approach relates process parameters such as laser beam power, process speed, and mold temperature, regardless of the geometrical characteristics of the material being processed. This analytical approximation makes assumptions such as uniform heat flux exposition, no heat loss to the environment, and constant thermophysical properties.
This work focuses on the development of a digital twin model for an automated tape-laying (ATL) machine. The model takes into account the heat source, the composite material, the compaction roll, and the surroundings as an enclosure to study the temperature distribution along the material under different process conditions. The mathematical model is 1.5 dimensional, as both the material and the reflector thickness are modelled with only one element in that direction. The thermo-optical properties of the elements with its temperature dependence were also considered within the infrared range to obtain an accurate description of the process response.
In Section 2, the ATL machine head components are described. Section 3 presents a 1.5 dimensional mathematical model for the components described in the previous section, emphasising the radiation heat exchange model, which involves optical and temperature dependant properties as well as their spatial distribution. In this, a mathematical model for the material machine processes is also presented. Section 4 presents the thermal properties equations for the machine head elements as a function of temperature. Section 5 presents the composite material to be used in terms of its thermal and optical properties; it also presents the equation solver method to simulate the ATL process and describes a strategy that can be used to measure the ATL process variables needed to feed the mathematical model. Section 6 presents the model validation by comparing the obtained measurements of the ATL process variables with the results calculated using the mathematical model.

The ATL Machine Head
The main components of the machine head assembly, numbered from 1 to 5 in Figure 1, are: (1) a material feeder, consisting of an unwinding mechanism that keeps the material straight along its pathway through the heating section; (2) a heat source, consisting of an infrared lamp and (3) a back plate that acts as a reflector; (4) an optical temperature sensor that is used to measure the material temperature; (5) a compaction roll in which a fluid circulates at a controlled temperature and presses the heated material against (6) a mold. The goal of the machine head is to create a laminate by heating up a continuous material tape to a temperature above its melting point that is to be welded to a previous material layer or the mold (if it is the first layer of the process) near the nip point, as shown in Figure 1.

Infrared Heating Element
The infrared heater in the machine head assembly consists of a tungsten filament coil heated by an electric current that serves as the emitter. This emitter is surrounded by a quartz glass envelope. This quartz glass acts as an enclosure containing an inert gas to prevent filament oxidation at high temperatures as well as small amounts of halogen to inhibit the evaporation of tungsten via the halogen cycle.
To maintain the atmosphere inside the quartz envelope, the current is passed through a pair of Molybdenum foil pads, ensuring the gas chamber is sealed to the outside atmosphere, for an increased service life. Figure 2 shows the heating element assembly.

Tungsten Filament Inert Gas Quartz
Lamp leads Molybdenum foil

Compaction Roll
The compaction roll is a cylinder that rotates with the motion of the entire machine head assembly while applying a fixed amount of pressure on top of the new material layer and the previous laid one, or the mold surface for the case of the first laid material, ensuring a contact bonding between the two surfaces at the nip point. To avoid any damage being caused to the heated material by the cylinder, a 5 mm coating of silicone rubber from Silex [24] covers the 100 mm diameter aluminium cylinder, preventing it from sticking to the material matrix during the process. On the other hand, to avoid the overheating of the compaction roll silicon cover, a fluid is circulated through the metal housing of the cylinder to control its temperature, which, according to the silicon manufacturer, can reach a maximum of 300 • C.

Material Feeder
The material feeder consists of a cylinder with a tension mechanism that ensures the material roll placed around it will always be in tension while it is being fed towards the compaction roll.

1.5D Mathematical Model for the ATL Heating Process
As stated in the previous section, the ATL head is mainly made up of an infrared lamp acting as a heating element with an aluminum back plate acting as a radiation reflector and a consolidation roll that presses the heated material against the mold. The process involves two main phenomenons-namely, the heating process to achieve the adequate consolidation temperature in the incoming material and the consolidation process itself.
In this work, we focus on the heating process of the incoming material because this is the most critical to the process quality [8,11,25,26]. Both the reflector and the material will be modelled using a single element along its thickness. For the reflector, as its main function is to redirect the heating energy towards the material, it is not necessary to know the temperature distribution along its thickness; a compensation factor is implemented for heat losses due to convection at the surface facing the composite material and at the opposite surface. For the material, its thickness is significantly smaller than its length, meaning that the temperature gradient along its thickness is negligible. Hence, a 1.5D heat transfer model of the ATL head is used to develop a process digital twin, which would allow us to run simulations and design proper control strategies. Figure 1 presents the layout of the machine head assembly. It moves horizontally, parallel to the heated mold, and the material is fed vertically to the compaction roll, which changes its orientation to horizontal when the nip point is reached. This causes the material feed velocity to match the machine head assembly velocity. The other assembly elements, such as the infrared lamp, reflector, compaction roll, and optical temperature sensor, move horizontally along with the machine head.
Understanding the relative movements between the machine head components and the material is relevant, as this determines the selection of the appropriate correlation for estimating the convection coefficients. Figure 3 shows a schematic of the simplified model of the ATL presented in the previous section. The curved path was simplified to a straight line, where y convection corresponds to the region before reaching the roll, y roll corresponds to the curved path in contact with the roll with its equivalent length, and y mold corresponds to the region after the nip point.
The ATL model involves an open-cavity radiation problem where the heat emitted by the lamp and reflected by the reflector is used to increase the temperature of the entering material tape before it reaches the nip point. Ambient heat loss due to convection and radiation also takes place and is considered in the model. Figure 3 shows three different zones of heat exchange over the material: The first one, delimited by y convection , involves radiation from the lamp and convection to the surrounding air, for which the convection coefficient depends on the air temperature and air flow velocity and direction. The second one, delimited by y roll , involves heat exchange from radiation, convection, and conduction to the compaction roll. The third zone, delimited by y mold , involves only conduction to both the compaction roll and mold.
To develop the proposed model, all of the components are analyzed to obtain the governing equations and the corresponding finite differences schemes, following a similar approach to that presented in the literature by [18] to obtain a mathematical expression based on the literature, as follows.

Heater
The heater element is a halogen lamp composed of a tungsten filament enclosed by a cylindrical quartz envelope filled with neon gas. A simplification of the lamp geometry is presented in Figure 4, where the concentric arrangement of the components was proposed by [18]. Figure 4a presents a front view of the lamp model as the tungsten is modeled as a solid cylinder of length l w and diameter d coil as a simplification of its spiral construction (Figure 4c), while the quartz envelope is modeled as a hollow cylinder of length l q , thickness t q (Figure 4a), and external diameter d l (Figure 4b).

Tungsten Filament
The tungsten filament is considered as an homogeneous element with an internal heat generation given by its electrical resistance, which is a function of temperature T f , and the applied voltage. It exchanges heat with the neon by conduction and with the inner face of the quartz lamp by radiation. No convection between the filament and the neon is considered because of the low value for the Grashof number, as suggested by [18], indicating conduction to be the predominant phenomenon in this particular heat exchange process. The heat transfer model used for the tungsten filament includes the energy balance given in (1) [27,28], where: • The electric power of the resistance, can be expressed using (2): where U e stands for the applied voltage. The electrical resistance R(T f ) can be expressed as a function of the resistivity, filament length l f , and filament cross-sectional area a f (3).

Neon
The neon has a very low thermal capacitance compared to the tungsten and quartz; hence, its thermal capacitance is neglected. Accordingly, the neon transfers heat from the tungsten to the lamp internal surface by conduction. Taking advantage of the cylindrical configuration, the heat conducted by the neon can be approximated using the cylindrical solution for radial conduction (4), as suggested by [18,27]. where: • Q c,n : Conduction heat, W.

Quartz Glass Envelope
The quartz envelope is considered to have a constant temperature T q across the thickness because it has a small thickness. The energy balance for the envelope is given in (5), as suggested by [18]. where: • m q : Quartz cylinder mass, kg.
• cp q (T q ): T ∞ : Air Temperature inside the radiation cavity, K Some of the data required for the lamp simulation are expressed in Table 1 from measurements performed on a series of lamps.

Reflector
The lamp reflector is considered as a thin metal sheet with a constant temperature across thickness t r ; it is modelled as a 1.5D finite volume problem along its width W r . The reflector transfers heat with the ambient air at its right hand side surface (T ∞ , h 1 (T)), transfers heat with the internal air in the cavity (T ∞ , h 2r (T)), and reflects radiation heat from the lamp. An energy balance over a reflector cell of volume V r is given in (6), taking into account its surfaces S r for heat exchange. where: • ρ r : Reflector material density, kg m 3 . • T r : Reflector temperature, K.
• q conv,r (T r , T ∞ ): Convection heat W m 2 . A finite volume technique is used to discretize the lamp reflector geometry, as shown in Figure 5. The application of the energy balance to each reflector cell gives the cell equation (7) for the unknown cell temperature, with the coefficients given by (8) through (11). The material properties inside the volume cell are constant and are evaluated at the nodal temperature. with: where: • r: Reflector cell number 0,1,. . . ,n. • L r : Length of reflector cell, m. • t r : Reflector thickness, m. • T ∞ : Air temperature, K.
The lateral conditions for the reflector are defined as adiabatic (12a) (12b) due to W r t r , making the heat exchange in those surfaces negligible in comparison with the heat exchange of the remaining surfaces.

Material
The material fiber tape is placed at a velocity U over a mold at a constant temperature T mold while it is heated by the lamp. The material domain begins at the section where the material is fed (Z 1 ), exchanging heat by radiation and convection, followed by the section Z 2 , where heat transfer also occurs by conduction due to the contact between the material and the compaction roll. This ends with Z 3 , where the consolidation point its located.
The tape is considered as a thin sheet with a constant temperature across its thickness and supplied at a constant temperature T 0 . The energy balance for a cell volume V m of the material tape is given in (13): where: A finite volume technique, as suggested by [13], is used to represent the tape geometry (see Figure 6), and the application of the energy balance to each tape cell gives the cell equation (14) for the unknown cell temperature using the coefficients given by (15) through (21). where: •  (14), taking into account each material section Z i , as shown in Figure 6. This allows us to obtain the corresponding boundary conditions for the cell. After evaluating all the terms for all the material cells, a system of differential equations is built to compute the new material temperatures.
At section Z 1 , the material arrives at an initial temperature T 0 = T ∞ , where convection (C1, C2) and radiation (D3) are the main heat exchange phenomena. At this section, the convection's relation is defined using the inside enclosure air properties (T ∞ , h 2m ), and the radiation phenomena are added to the balance to complete the relation. The radiation term added comes from an energy balance further discussed in this work; at section Z 2 radiation (D3), convection (C1), and conduction (D1) are the heat exchange phenomena present; at section Z 3 the material is located between the compaction roll and the mold, which means that conduction (D1, D2) is the main heat exchange phenomena occurring between the compaction roll and the mold simultaneously. For this model, the temperature of the nip point is located at the material cell T n , and the previously laid material temperature is modelled assuming a constant value of T mold . The terms R mold and R roll are explained in the following section.
The boundary condition related to the y axis at the first material cell is represented using (22); the last material cell, in the same y axis, is represented using (23).

Compaction Roll and Mold
The compaction roll is used to press the incoming material, once heated, against the mold or the material previously laid during consolidation. It consists of a hollow aluminum cylinder in which water circulates at a constant temperature to avoid mechanical damage to the seals and bearings. The exterior of the cylinder is covered with a layer of solid high-temperature silicone rubber Silex GP60THT [24].
The model strategy used for simulating the heat exchange process between the compaction roll and the material uses a combination of a composite hollow cylinder with 1D radial conduction for the compaction roll internal fluid, the aluminium structure, and the rubber layer and the 1D conduction of a plane wall approximation for the material, as presented in [27], because the material is thinner than the rubber layer or even the aluminium structure of the cylinder, making the relation ln D roll /2 + t rubber + ∆x where D roll is the compaction roll diameter without the rubber layer, t rubber is the thickness of the rubber layer, and ∆x is the thickness of the material.
For the case presented at section Z 2 in Figure 6, the heat exchange phenomena are convection, radiation, and conduction. The energy balance for conduction to the compaction roll at this section can be obtained from Figure 7a, where the heat exchange by conduction is mathematically expressed as an array of resistors connected in series relations (24) and (25).
From Figure 7a, the term R roll is deduced, taking into account the convection coefficient of the internal fluid (giving relation (24)). For a constant internal temperature, the relation can be deduced (25).
At section Z 3 in Figure 6, the scheme for the energy balance involves conduction through the mold, as presented in Figure 7b. For the mold, a film of Polytetrafluoroethylene (PTFE) acting as a mold-releasing agent between the laid material and the moldit is considered; from Figure 7b, the term R mold gives the relation (26):

Radiation
Radiation as heat transfer has been well studied when gray surfaces participate inside an enclosure; it has been described in several text books [27][28][29][30] and by several researchers [18,[31][32][33]. For this work, the general enclosure is described in Figure 3, where the reflector, the lamp, and the material are involved. The dashed lines represent a hypothetical surface acting as the surroundings of the enclosure. In order to simplify calculations, the radiation flux problem is modeled as a two-dimensional system allowing one to calculate the heat flux per unit area; Equations (7) and (14) can be solved using the radiation heat input expressed in (W/m 2 ). It is worth mentioning that the view factors for two-dimensional problems have been well studied and equations have been developed to simplify the calculations [28][29][30]33,34], these equations are discussed further in this work.
This work also considers another enclosure, in this case a lamp, which is composed as described in Section 3.1 by a tungsten filament inside a quartz glass envelope. In this enclosure, the filament exchanges heat, mostly in the form of radiation with the inner surface of the quartz envelope. At the same time, part of this radiation is also exchanged with the external glass surface due to the glass being a transparent medium; this is fundamental when defining the radiosity balance among the surfaces involved in the general enclosure.
To calculate the radiation energy of a gray surface, first the radiation energy for an equivalent black body is calculated using the Planck's Law (27); then, a correcting factor that is specific for each surface nature and temperature is applied to transform the black body radiation into a gray body radiation. Those correcting factors are discussed later in this work. with: where: • λ: Wavelength, µm. For this work, the radiation energy is focused on the infrared region of the spectrum, narrowing the calculation region between λ 1 and λ 2 (0.4 µm and 20 µm, respectively); then, the total energy for a black body is obtained using (30) from Planck's Law (27).
As the energy calculation depends of the wavelength and the integration must take place along a finite band of wavelengths, the numerical solution is achieved evaluating the first 10 terms of the infinite series (31), where σ stands for the Stefan-Boltzmann constant, 5.670 · 10 −8 W m 2 ·K 4 . By dividing the studied wavelength interval into 20 equally distributed evaluation domains, it is possible to obtain a good approximation while improving the calculation time in comparison with numerical integration routines [28,31]. with:

Radiative Fluxes
To define the radiative balance, an identification of the surfaces involved is presented in Figure 8, as well as the radiosity, J i , and irradiation, G i , notations for each surface with their respective indices.

Material
Reflector Lamp The boundary condition labeled as 6 is defined as a black body representing a surface that encloses the enclosure. This allows us to determine the net radiating heat transfer (34) using relation (35) for each surface.

Black Body Surroundings
The first step is consider the radiosity for each surface, giving the set of relations (36) through (41): as suggested in [18,34]. In generic matrix form, for n involved surfaces, the relation is given by (42): where relations (37) and (38), which are related to the quartz envelope, contemplate the transmission τ λ,i between both the inner and outer surfaces due to the glass being a semitransparent medium. The term r λ,i refers to the reflectivity of the body and the term ε λ,i refers to the emissivity; all surfaces are considered diffuse gray bodies. Relation (41) refers to the surroundings considered as black bodies, acting as a heat sink as the temperature is known and kept constant over time. According to literature [18,28,29,35], it is defined r λ,6 = 0. The second step is to consider the irradiation balance on the different surfaces by using (43): where F j−i is the view factor from surface j to surface i and A i and A j are the areas of surfaces i and j, respectively. From the view factor reciprocity property (44), the expression (43) can be written as (45) or in its compact form (46).
In matrix form, for n involved surfaces the relation can be written as (47): In order to solve this equation system, Equations (42) and (47) are combined, leading to (48) and enabling the computation of the irradiation.
Then, the results obtained using (48) are used in (42) to calculate the radiosity J λ,i . This allows the computation of the net radiative flux, leaving each involved surface per unit area defined as (34), which are the boundary conditions needed to solve the lamp (1), (4), and (5); the reflector (7); and the material (14) temperatures.
As mentioned before, the lamp is modeled as a solid tungsten cylinder inside a quartz glass envelope (enclosure); thus, the view factor for infinitely long concentric cylinders can be expressed as (50), (51), and (52), whose indices refer to the notation shown in Figure 8.
The view factor for an infinitely long cylinder to an infinitely long rectangle can be obtained using the scheme presented in Figure 9a, with the one presented in Figure 9b being a particular case when b 1 = b/2.
The view factors between two planar surfaces, in the 2D case, can be expressed using Hottel's cross string method. Figure 10 shows the different components of this method. The values L 1 , L 2 , L 3 , L 4 , L 5 , and L 6 are lengths, with L 1 and L 2 being the surfaces involved in the view factor and the remaining values distances between the edges of the surfaces-in this case, the cells of the finite volume method. The view factor obtained using the Hottel's method is expressed as (54) for the case F L 1 −L 2 and (55) for the case F L 2 −L 1 .

Properties
The properties of the elements involved in the heat exchange process are temperature-dependant and can be classified as thermal or optical properties, directly affecting the heat flux distribution model.

Thermal Properties
In Equations (1), (4), and (5), the specific heat of the tungsten, cp w (T f ), as well as the quartz, cp q (T q ), are needed. These data were collected from the literature [18,29,37] and summarized in (56) and (57) for the tungsten with a temperature range of 298 K < T < 1900 K and 1900 K < T < 3680 K, respectively, and for the quartz (58) with a temperature range of 298 K < T < 2000 K [29,38,39].
Additionally, the tungsten resistivity and thermal conductivity of the neon, which are both temperature-dependant, are required. From [18], the data for the tungsten resistivity are given by (59) for a temperature range 300 K < T < 3655 K, while the conductivity for neon gas is given by (60) for a temperature range of 100 K < T < 2500 K.  (7) requires the specific heat of the reflector, cp r (T r ), whose expression is presented in (61) for a temperature range of 298 K < T < 933 K [37,40,41]. Its thermal conductivity, k r (T r ), is expressed in (62) for a temperature range of 250 K < T < 800 K [42].

Optical Properties
The optical properties needed by the system of Equation (42) for each surface type are wavelength-dependant. For general surfaces, the relation among its optical properties is expressed by (63); for the case of opaque surfaces which have no transmissivity τ = 0, the relation among its optical properties is expressed in (64) [27,41].
To estimate the optical properties of the opaque surfaces, the Maxwell electromagnetic wave theory gives the emissivity value (65), which, when used in (64), leads to the reflectivity: where n and k are given by Equations (66) and (67), respectively: with: • 0 : relative material i permitivity.
. For a filament the resistivity is given by (59); for a reflector the resistivity is given by (68) [42,43].
For the tungsten filament, the value for the relative permitivity is / 0 = 1.000068 and for the reflector / 0 = 1.00000065 [44].
The optical properties of the quartz glass envelope are gathered from specialized literature on the study of infrared dryers and from papers specializing in the study of material optical properties [34,45].
The emissivity values used in (41) is ε λ = 1 for the surroundings are considered as black bodies.

Convection Coefficients
For the convection parameters needed for solving (5), (7), and (14), empirical correlations are used to estimate h 1 , h 2r , h 2q , h 2m , and h 3 [27]. These correlations estimate an average convection coefficient for the air h i by calculating an average Nusselt number Nu. The properties for the air at atmospheric pressure are computed at T f ilm (69) according to the information presented in the literature [27].
To estimate h 1 , the set of correlations (70), (71), and (72) for perpendicular air flow were used, where the value U of the movement of the machine head is used to calculate the Reynolds number.
For the left-hand side of the reflector h 2r , the correlations for a vertical plane with free convection (73), (74), and (75) were used: where g = 9.807 m/s 2 . For the lamp h 2q , the correlations for a cylinder with free convection were used (76): (77) and (78). The position of the reflector relative to the lamp blocks the air flow due to the machine head assembly movement.
For the case of h 2m and h 3 , the same correlations are used-(79), (80), and (81)-which consider an external air flow of velocity U parallel to the material surface due to the material motion relative to the air, which is assumed be in a quiescent state.

Material Description
As a case of study for the ATL process, the employed material is a carbon fiber-reinforced Polyamide 6 tape from Toray with the commercial denomination Cetex TC910 ® . This has a 50 mm width, 0.16 mm thickness, and a fiber content of 60% [46].

Thermal Properties
Equation (13) requires the material heat capacity as a function of temperature, cp m (T m ). An experimental procedure was performed according to ASTM E 1269-9901 [47] for material samples, giving the results shown in Figure 11 for a temperature range of 300 K < T < 508 K. Experimental data Another property required in (13) is the material conductivity, k m (T m ). An in house procedure was used to determine this property as a function of temperature based on theoretical and practical works by [48][49][50][51]. The resulting data are presented in Figure 12 for a temperature range of 300 K < T < 503 K.

Optical Properties
In order to obtain the emissivity, ε m (T m ), and reflectivity, r m (T m ) of the composite material as a function of temperature, an in house test was designed based on the procedure indicated by pirometer manufacturers to determine the emissivity of the object being measured [52]. The results are presented in Figure 13.

Time Integration Scheme
The mathematical model involves a set of non-linear ordinary differential equations (ODE) for the temperature time derivative of the different components-namelym the tungsten filament (1), the quartz lamp (5), the reflector cell temperature (8), and the material cell temperature (14). To solve this ODE system, the initial condition is set to be the ambient temperature for all the components. The implicit solver Radau IIA is selected from the Scientific Python library [53] to avoid time-step limitations due to equation stiffness. The ODE equation system is structured in a array-like way as an explicit function of the given trial temperature suitable for the external solver. To evaluate the ODE system, the following steps are performed: Step 1.
For the given trial temperature.
Evaluate all material properties at the trial temperature.
Solve the radiation heat flux system for the trial temperature.
Step 4. Evaluate the ODE system array with the temperature derivatives.
The model solver consists on two parts, with the first one being the calculation of the net radiation flux using the current temperature as a reference to determine all the properties required. The second part consists of implementing the result of the net radiation flux as an input to predict the process component temperature at the next time step. Finally, the new temperatures are used to calculate a new net heat flux and then the cycle repeats until the defined stop time.
The relative error tolerance for the solver is set to 10 −4 and the absolute error tolerance to 10 −6 .

Convergence Analysis
A convergence analysis is performed to define a proper cell size for the material simulation in order to lower the required computational resource consumption using the lowest mesh size. Mesh sizes from 11 cells up to 239 cells were tested for the material. The mesh for the reflector was defined using the same cell size used for the material, regardless of the amount of cells, in order to preserve the stability during computation. Taking into account the temperature percentage variation using three selected points as a reference-namely, the first cell at which the material has its first interaction with the heat exchange process, the middle cell along the length of the material, and the last cell of the material-the mesh size is selected based on the lowest number of cells which gives a temperature variation under 0.5% with respect to the mesh size of 239 cells. The number of cells defined using odd values is intended to guarantee a central cell along the mesh at the same position, regardless of the mesh size.

Measures and Instrumentation
In order to measure the process parameters and the variables needed to feed the mathematical model, two types of data acquisition cards from National Instruments were used-namely, an NI 9234 for voltage measurements and an NI 9211 for temperature measurements. The NI 9234 data acquisition card allows high sampling rates and is adequate for measuring the 50 Hz frequency mains voltage that feeds the heating element.
To compute the electrical power delivered to the heating element, the applied voltage and the current consumption must be measured. A voltage transformer is used to step down the mains voltage sine wave to a safe range for the data acquisition card, as well as to provide galvanic isolation. For the current, a non-invasive current transformer is placed around one of the lines that powers the heating element; the output of this transformer is a voltage signal that is read by the same data acquisition card that reads the transformed mains voltage signal. The measurement strategy used for the voltage and current is presented in Figure 14. To measure the temperatures involved, a set of thermocuples is placed at the following locations: the compaction roll fluid inlet and outlet, and on the top surface of the mold. The composite material temperature is measured using a pyrometer PyroNFC K from CALEX facing the composite material at a point h above the center line of the compaction roll, as presented in Figure 15. The output of the pyrometer is a thermocouple-like signal that can be read by the same data acquisition card, NI 9211, as the other thermocouples. The last parameter to measure is the process speed; this value is obtained directly as a voltage from the analogue output of the servo motor driver OMROM R88D-KN08H-ECT and is read by the NI 9234 data acquisition card.

Model Validation
The experimental validation was performed in two parts. The first part consists of validating the temperature dynamics of the material at the measuring point, as shown in Figure 15, using a constant set of parameters-namely: process speed, compaction roll temperature, and mold temperature-and imposing changes in the voltage applied to the heating element. This voltage variation is conducted using a solid state relay that takes advantage of the phase angle strategy to reduce the delivered root mean square voltage and hence the electrical power delivered to the heating element.
The compaction roll temperature and the mold temperature are kept constant because any change to their respective reference generates a slow response in comparison with any change in reference for the voltage delivered to the heating element. The process speed is kept constant throughout the experiment because the driver controlling the servo motor only has signal outputs to monitor its value and any change will likely stop the process and change it manually. In order to see the influence of the speed upon the process, two sets of tests were carried out with different speed values, as presented in Table 2. To understand the machine parameters' influence on the operation of the ATL head, regarding the consolidation temperature of the material tape at the nip point, the proposed test cases presented in Table 3 will be simulated.

Mesh Convergence
The results obtained for the convergence analysis are shown in Figure 16, for a speed value of 5 mm/s in Figure 16a, and for a speed value of 15 mm/s in Figure 16b, both with a simulation time of 30 s and the same process conditions of initial temperature and electrical power delivered to the heating element, as well as compaction roll, mold, and ambient temperature .   11  17  23  29  35  41  47  53  59  65  71  77  83  89  95  101  107  113  119  125  131  137  143  149  155  161  167  173  179  185  191  197  203  209  215  221  227  233  First node Central node Last node (b) Figure 16. Temperature error as a function of the mesh size for the material. (a) speed value 5 mm/s. (b) speed value 15 mm/s. The data presented in Figure 16a show that, for a process speed of 5 mm/s, the material temperature value, at the selected locations, has a tendency to stabilize with a relative error below 0.5% for mesh sizes above 47 cells, and for the case of a process speed of 15 mm/s, the error below 0.5% is achieved with a mesh size above 35 cells. From this results, the selected mesh size is 51 cells, corresponding to the scenario presented in Figure 16a, with a cell size of ∆y = 2.87 mm.
The reflector cell size is fixed to have the same size as the material, although the reflector mesh size should have less influence.

Model Validation
In order to compare the dynamic response of the ATL machine with the predicted response of the proposed numerical model, different tests which parameters were defined according to the values presented in Table 2 as Set 1 and Set 2, were performed with multiple changes in the voltage supplied to the heating element. The results for the parameters Set 1 are presented in Figures 17 and 18, and for the Set 2 in Figures 19 and 20. This was carried out to simulate a possible control action and to analyze the transient response of the machine and the material. Measured Temperature Model Temperature (d) Figure 17. Results of Power consumption and temperature variation at the measuring point according to parameter set 1 from First, the predicted power consumption agrees with the experimental measured response, as shown in Figure 17a,c. This means that the proposed model is capable of predicting the overall energy flow of the system, which is a desirable result because it is based on energy conservation laws.
Comparing the measured temperature with the predicted temperature at the corresponding "machine measurement point", it is verified that the proposed model can predict the heating and cooling dynamics of the system, as shown in Figure 17b,d, according to the process speed value, compaction roll temperature, and mold temperature presented in Table 2 as Set 1.
The prediction of the material temperature was adjusted by defining a compensation factor for the convection coefficient of the heating element of 5.5, correlation (78), and a compensation factor of 2.5 for the convection coefficient of the material, correlation (81), in order to compensate the correlation approximations related to geometry and relative movement.
Additionally, the corresponding measurements for the compaction roll temperatures, mold temperature, and process speed for test set 1 are shown in Figure 18a The test and simulation results for the process parameters presented in Table 2 as Set 2 are presented in Figure 19, showing the predicted and measured power consumption for the heating element, Figure 19a,c, and predicted and measured temperature for the material at the measuring point, Figure 19b,d.
The model used to predict the composite material temperature at the machine measuring point shown in Figure 19b,d includes the same compensation factors in the convection coefficient for the heating element as well as in the material convection coefficient, showing an acceptable geometrical and relative movement compensation for those values.
The corresponding measurements for the compaction roll temperature, mold temperature, and process speed are shown in Figure 20.
As shown in Figures 17b,d and 19b,d, the model is capable of predicting the temperature of the material for values above 150 • C (423.15 K).
Comparing the measured temperatures at the same machine measuring point for both process speeds, it is noticed that when defining a machine speed of 15 mm/s, temperatures at this point are lower than the temperatures at the same point when tests are performed with a machine speed of 5 mm/s, showing the influence of the transport phenomena on the temperature evolution over time.
It can be seen that the compaction roll inlet and outlet temperatures are similar at both tests; this fact, added to an internal fluid flow through the compaction roll of 40 L/min produced by a temperature control unit Tool-Temp 137BP, indicates that the internal flow is sufficient to guarantee a constant temperature along the compaction roll, validating the proposed model assumption of constant temperature given in Section 3.4.   This result is expected and is related to the transport phenomena due to the feed speed and the heat conduction to the roll.

Sensitivity Analysis
The results obtained for the simulations described in Table 3 are presented in Figure 22 for the simulation in which the speed of the process is kept constant (Simulation 1) and in Figure 23 for the simulation in which the voltage applied to the heating element is kept constant (Simulation 2).  Table 3.
For the case presented in Figure 22, the temperature at the nip point also presents a delay with respect to the measuring point as expected. This is controlled by ambient losses and the heat transferred to the roll. Figure 23 shows that the process speed has a relevant influence on the nip point temperature as its value increases; this effect can be noticed for the simulated speed of 10 mm/s. After the nip point temperature starts to rise, it almost matches the temperature of the measuring point, and for the case of a speed of 15 mm/s the temperature experiences a start delay and rapidly reaches higher values with respect to the measuring point. This behaviour demonstrates the high sensitivity of the system to the feed speed, meaning that the transport phenomena dominates over the conduction and convection losses.  Table 3.

Conclusions
The energy balance model proposed for the heating element is able to replicate the dynamics of its electrical energy consumption, taking into account its components' properties and geometrical parameters. It is also capable of predicting the outgoing energy flux towards the material being processed, causing the temperature to rise along the material length.
The compaction roll can be modelled assuming isothermal conditions along its internal fluid section by ensuring a sufficient fluid flow, causing a temperature gradient along the compaction roll thickness. This assumption simplifies the heat transfer model for this element.
The proposed 1.5D heat transfer model for the automated tape-laying machine head can be used to simulate the temperature distribution along the processed material and its dynamics upon changes in the heating element input voltage and process speed, predicting the temperature along the heating zone domain.
The most relevant parameters that influence the machine behavior are the input voltage of the heating element, which increases the temperature of the processed composite material by directly heating its surface using radiation, and the process speed, which influences the material temperature through the transport effects.
The current machine configuration has an appropriate design and is capable of raising the material temperature to the required consolidation temperature value if a suitable control system is implemented.
With this model, which takes into account the thermal and optical properties of the material involved and their changes with temperature, a digital twin of the process can be formulated to develop control strategies using the machine head speed and the lamp power as the variables manipulated in order to control the nip point temperature. This model also allows one to estimate the composite temperature profile along the heating zone, acting as a supervisory function for a control system and allowing it to handle imposed constraints related to the material.
As the digital twin of the process is capable of estimating the composite temperature at the nip point, it lends itself to creating a temperature registry over time for further thermo-mechanical analysis for analyzing defects due to residual stress. The digital twin model can be used for selecting the best process elements and their capabilities in terms of the heat source delivering capabilities, the process speed, and the mold temperature as a function of the composite material being processed in order to meet quality requirements of the final product. Thus, the study of those effects on an already built machine can be avoided and large time and economic benefits can be obtained.