Optical thermal model for LED heating in thermoset-automated fiber placement

Abstract Control of material temperature distribution and governing phenomena during automated fiber placement is an important factor. Numerical modeling of the radiative heat transfer for a newly presented LED-based heating unit is developed and analyzed in theory. An optical model allows taking into account the radiative energy output of every individual LED. By adjusting the electrical input to the multiple LED arrays on the heating unit, the irradiance distribution on the substrate’s surface can be controlled. To investigate the capability to adjust the surface temperature distribution resulting from this feature, thermal models for two and three dimensions are developed and employed for the calculated irradiance distributions. The resulting temperature distributions show that temperature gradients can be avoided or created, depending on the input to the heating unit. The results from the two models are compared and a method to select an appropriate model in general is proposed.


Introduction
The latest generation of long range passenger aircrafts, especially Airbus' A350 and Boeing's 787, feature more than 50% of advanced composite materials. In manufacturing of the large structural components such as fuselage panels and wing covers made of carbon fiber-reinforced plastics, automated fiber placement (AFP) of thermoset prepreg plays a key role [1]. In these processes, the controlled input of heat is one of the most important parameters [2,3]. Heating the material immediately ahead of the compaction roller to approximately 35 À 40 C is required to increase the material's tack to reliably adhere the newly placed tows to the substrate [2,4,5]. Historically, all mechanisms of heat transfer have been used for the heating, but today usually radiative heat transfer by means of infrared lamps is employed [2,3,5]. Nonetheless, even though fast reacting medium wave infrared lamps are used, their reaction times limit the achievable process control and lead to constructive measures to prevent overheating and possible damaging of the material in case of an emergency stop [4,6,7]. Furthermore, the distribution of irradiation on the surface, the irradiance, is inhomogeneous, not adjustable and overshooting the current placement path [4,8].
Therefore, a new radiative heating unit based on LEDs as energy sources has been proposed [9,10].
Similar to other radiative heating technologies, their electromagnetic radiation was shown to be well absorbed by the prepreg material and thus converted to heat [9,10]. Further, as LEDs are semiconductors, they offer negligible reaction times for automated layup [10]. Similar to laser diode arrays and related technologies, utilizing a multitude of LEDs arranged in individually controllable arrays allows to adjust the radiation output distribution [9][10][11][12]. This basic capability was demonstrated by first experiments in Orth [10], leading to the suggestion to employ numerical modeling for a detailed analysis in a second step. In this context, optical modeling allows calculating the irradiance on the material surface and subsequently, this irradiance distribution can be used in thermal models to predict the substrate's temperature distribution [2].
Optical models based on the geometrical view factor have been presented in the literature for infrared lamps [8,[13][14][15]. However, due to the infrared lamps' nature, these assumed a thermal radiator with a fixed output pattern and thus intensity distribution. The concept of utilizing a multitude of individually controllable LED arrays requires taking into account each LED individually and in order to control the irradiance in both directions of the substrate surface, a three-dimensional approach for the optical model is required.
The thermal aspect of numerical modeling in automated layup processes has been the subject of many studies, with the majority of models focusing on thermoplastic prepreg layup heated with lasers [2,5,13]. However, many of these models used a two-dimensional approach, neglecting the dimension along the length of the compaction roller as the irradiance was considered homogeneous in this direction [2,[16][17][18]. This was shown to be acceptable through three-dimensional optical modeling and experimental validation by Stokes-Griffin [16]. The new LED heating unit, however, was shown to create and control temperature gradients in this direction in experiments [10]. On the one hand, a threedimensional approach to the problem, therefore, seems natural, as neglecting it was considered a possible cause of deviation between experiment and the two-dimensional model for thermoset AFP by H€ ormann [13]. On the other hand, as two-dimensional models often are more efficient and have seen great use, it is of interest to investigate the differences between both for the given problem.
Therefore, following up on the experiments in Orth [10] it is the first aim of this study to develop a suitable numerical modeling approach for the new LED heating technology. For the optical component, different electrical power input distributions are analyzed and the resulting surface irradiance distributions compared. Subsequently, a two-and a three-dimensional thermal model are developed, their simulation results compared with each other and the findings discussed. Based on these findings, the second aim of this study then lies in developing a tool to help in choosing the required modeling approach in heat transfer during automated layup, but regardless of the applied heating unit.

Optical model
The purpose of the optical model is to calculate the conversion, transfer and absorption of the electrical input energy into each individual LED to the resulting heat input on the material's surface.

Modeling assumptions
The net radiative energy transfer rate _ Q 1!2 between two black body surfaces A 1 and A 2 is generally defined by [19] with F 1!2 as the view factor from surface A 1 to A 2 , r is the Stefan-Boltzmann constant and T are the surfaces' respective temperatures. Several assumptions are required for the application to the stated topic of a new LED-based heating unit. First, in contrast to infrared lamps, LEDs are non-thermal radiators, thus their electromagnetic energy output depends on the applied electrical input to an LED and the device's characteristics such as conversion efficiency, rather than surface temperature. Neglecting the influence of device temperature on conversion efficiency the radiative power output P LED of the OSRAM OSLON Square LD CQAR LEDs considered for this study can be assumed to be linearly dependent on the input current [20]. Second, comparing the radiation characteristics given by the LED's manufacturer to a perfect cosine function justifies the assumption of a Lambert radiator, even though it features a lens [20,21]. Third, the substrate is also assumed to be a grey radiator, as its absorptivity 2 remains almost constant up to an incident angle of over 70 and the change in absorptivity due to the temperature's influence is small [17,22,23]. The net heat input _ Q LED;1!2 from an LED to the substrate is, therefore, simplified to ignoring input from second-order reflections and radiative losses to the environment due to low reflectivity and low difference in temperature, respectively.

Model implementation
For this study, an LED heating unit in a placement situation at constant speed and on a flat surface is considered, as this allows focusing on the relevant effects. The optical model is symmetrical around the plane parallel to the x-z plane (Figure 1, left) and its geometry is based on the information given by Calawa [4], which was adapted to fit to a lab scale demonstrator currently under development. As the heating unit is aimed at an angle a away from the compaction roller, no significant shadowing effects are to be expected and its presence is, therefore, neglected. The utilized LED heating unit features 18 Â 12 LEDs in the lengthwise and crosswise directions, respectively, arranged in 10 individually controllable arrays as marked by the black rectangles in the cross-sectional view A-A (Figure 2, left) and denoted AR, AL, BR, BL, Z1-Z6. The LEDs are spaced equidistantly in both directions and the die area is centered in their contour, forming the discrete emitting elements for numerical simulation (Figure 2, right). All parameters are given in Table  1, including element sizes of emitters and substrate. The utilized nomenclature is detailed in Table 2. Discretizing the problem for both the emitting LEDs and the irradiated CFRP substrate allows an implementation in MATLABV R by calculating the radiative energy input from every LED into every substrate element.
The first step is the calculation of the view factor between the two surfaces, generally defined for an arbitrary configuration as As the LEDs emitting area A D ¼ dA 1 and the substrate element size dA 2 are small compared to the distance s between them, applying the differential view factor given by can be justified and, therefore, is shown in Figure 3 [19]. The radiative energy input from Equation (2)     (5) where u i;j ¼ ½0; 1 is the linear scaling factor for the input current. This current is applied individually for each of the arrays shown in Figure 2 and thus all LEDs within. P max ¼ 2:55 W is the radiative power output at the maximum input current of an LED. Finally, the heat flux _ q or absorbed irradiance E 1 for substrate element k, l is obtained by summation over all LEDs repeating the procedure for all substrate elements gives the desired irradiance distribution.

Optical simulation results
To evaluate the new heating unit's capability to adjust the surface irradiance distribution simulations for three different electrical input, power distributions were performed. In the first setup, all LED arrays were operated at the same level, thus making its output homogeneous and most similar to an infrared heating unit with reflectors [4]. The other two setups provide a heterogeneous electrical input power distribution, aiming to reduce and create irradiance gradients, respectively. For all three setups, the linear power scaling factors applied for each of the arrays shown in Figure 2 are given in Table 3, as well as the total electrical power input, which remained constant.
In the homogeneous setup, all arrays are operated at the same input level. As expected, this causes an inhomogeneous irradiance distribution over the width of the placement path w tp (red dotted lines in Figure 4) in the x-y plane in Figure 1. The irradiance contour lines are displayed at multiples of 5000 W=m2 from 0 to 50000 W=m2 and highlight a concentration in the center area directly underneath the heating unit (black dotted rectangle in Figure 4).
The first heterogeneous setup aims to compensate for this localized effect by providing a higher output power near the edges of the placement path and less in the center area. The selected power factors were determined iteratively and are not optimized through an algorithm. Nonetheless, the resulting irradiance distribution in Figure 5 shows no concentration inside the placement path and the innermost contour lines for E ¼ 2 Â 10 4 W=m 2 form an almost rectangular spot around the placement path.
The second heterogeneous setup follows the opposite approach. By only utilizing the inner arrays (Z1-Z6), the irradiance in the center area is even higher than for the homogeneous setup. Nonetheless, the edges of the placement path are exposed to a similar level of irradiation, as the high intensity in the center spreads up to the edges ( Figure 6). Although the resulting strong gradients are usually not desired in an actual placement situation, they are well suited to study the influence of irradiance gradients on the thermal model and the heat dissipation along the y-axis.
For better comparison of the irradiance along the y-axis for the spot of highest concentration at x ¼ 0:075 m are plotted for all setups in Figure 7. While the irradiance levels at the edges of the placement path are similar for all input setups, the homogeneous and heterogeneous 2 distributions cause a 60% and 124%, respectively, higher irradiance in the center than at the edges. The heterogeneous 1 setup provides a much more homogeneous distribution inside the placement path and only a 7% higher  irradiance in the center. Comparing the edges y ¼ 0:05 m and at the center/symmetry y ¼ 0:075 m over the x-axis in Figure 8 shows that the heterogeneous 1 setup also provides a homogeneous irradiance distribution along the x-axis in the direction of movement.

Thermal model
To study the effects of the different irradiance distributions on material temperatures during placement, both a two and three-dimensional thermal model are developed.

Governing equation
For this study, operation at a constant placement head movement speed is assumed, neglecting transient effects and thus resulting in a steady flow. Therefore, an Eulerian framework with a coordinate system fixed to the heating unit is used, avoiding constant updating of boundary conditions [2]. The heat release from exothermic reactions is neglected, as no curing is expected for temperatures below 65 C. Therefore, for movement in the x-direction and assuming the conductivities k to be constant over their relative direction, the governing equation for the three-dimensional case becomes For the two-dimensional case, the y component is omitted, thus neglecting diffusive heat transfer in this direction. Furthermore, only the irradiated substrate ahead of the nip-point is considered, neglecting heat flow from the substrate to the incoming tape and the compaction roller.

Boundary conditions
As the heating unit and the resulting irradiance distribution are symmetrical, a symmetry plane is introduced for the thermal model to reduce computational efforts. The locations of the boundary conditions and the symmetry plane are displayed in Figure 9.
The boundary conditions for x and z directions in regions S1 to S3 and S5 are well established for     two-dimensional models from the literature [15,16,24]. In region S4, no gradient is assumed, as it is sufficiently far away from the irradiated area. The symmetry is imposed by also setting the temperature gradient perpendicular to the symmetry plane to zero. All boundary condition and applied parameters are listed in Table 4.

Model implementation
Due to the straightforward geometry of the domain, a finite-difference approach with a uniform discretization may be selected, using forward and central differences for the advective and diffusive terms, respectively. As the initial governing equation contains no time derivative, a linear system of the form can be obtained, where K is a matrix containing the influence factors, the vector t holds the temperatures at every point and R the boundary conditions. Although this system can be solved directly, an iterative approach for the implementation in MATLAB V R is used in this study. In the case of a one-dimensional problem, K is tridiagonal, holding the influence factors of each element's two adjacent neighbors and of size m Â m, where m is the number of elements. However, for two-and, especially, three-dimensional considerations, a very sparse two-dimensional matrix is obtained and the non-zero values are not as clustered around the main diagonal anymore [25]. Therefore, in order to apply the simple structure from the one-dimensional approach to the three-dimensional problem, K and t are implemented as threedimensional matrices m Â n Â l representing the discretized elements in the spatial directions. Further, each matrix element in K holds a 3 Â 3 Â 3 matrixK for the local influence factors and thus making application of an iterative Jacobi method easy [25]. The same approach is chosen for the two-dimensional model with the respective matrix sizes.

Thermal simulation results
All parameters used are presented in Table 5, including the element size of the discretization. As material behavior and influences on thermal modeling are not the focus of this study and found to be of minor significance compared to processing parameters [13], thermophysical properties were obtained from the literature and assembled to represent a virtual, yet characteristic thermoset material with temperature independent values [2,8,13,26].

Three-dimensional model
Thermal simulations were performed for the different irradiance distributions, with the same discretization in the x-y plane as in the optical mode. In the z-direction, it was assumed that placement occurs on top of seven previously placed unidirectional plies of the same orientation that are in perfect contact, featuring two more plies than the experiments in Orth [10] to avoid the influence from the tooling. Each ply consists of five elements to a total of 35 elements. Figure 10 shows a contour plot of the calculated surface temperature in the original x-y plane resulting from the homogeneous input power setup. Inside the placement path, temperature gradients in both x and y directions are visible reflecting the inhomogeneous irradiance. Most notably, the area of highest temperature is located more towards lower x values or "downstream" than the irradiance concentration in Figure 4. While the desired temperature of more than 35 C is reached everywhere inside the placement path, a maximum temperature of over 55 C is present.
The temperature distribution for the heterogeneous 1 setup in Figure 11 shows similar surface temperatures outside the placement path as Figure 10, but significantly lower gradients on the inside. Overall, the maximum temperatures are much lower, but the 35 C isotherm extends well up to the border of the modeling domain. Therefore, the desired temperature is achieved over the entire width of the compaction Figure 9. Thermal modeling domain (laminate) with coordinate directions for calculation, locations of heating unit, and compaction roller (not modeled), boundary conditions, and symmetry plane. roller at x ¼ 0 m to achieve uniform adherence of the placed material to the substrate. Finally, the heterogeneous 2 input setup causes very high gradients inside the placement path in Figure 12 that also extend up to the location of the compaction roller. These exist even though the heat release path in the y-direction is covered by the three-dimensional model. Therefore, it is of interest to compare these results with the simulations of the two-dimensional model.

Two-dimensional model
The two-dimensional simulations were performed for all three input setups, featuring the same discretization in x-and z-direction as for the three-dimensional model. As the irradiance is not constant over the y-direction, which is omitted in the two-dimensional model, two simulations at the most relevant locations y ¼ 0:05 m and y ¼ 0:075 m, the edge of the placement path and the center symmetry line, respectively, were performed for each setup. Thus, the utilized irradiance distributions were exactly these shown in Figure 8. The calculated temperature distributions for all setups and locations are given in Figure  13 and may be compared with the results from the three-dimensional model in the previous section, indicating similar trends. Again, the heterogeneous 2 setup causes the highest gradient between center and edge, which is lower for the homogeneous setup and significantly reduced for the heterogeneous 1 input.

Comparison of models
In fact, when extracting the temperature distributions at these discrete locations from the threedimensional models and comparing these directly to the ones obtained from the two-dimensional model, they align extremely well. Therefore, instead of overlaying the graphs, Figure 14 shows the difference between the two models at the two locations over x for the heterogeneous 2 setup, where the gradients and thus the expected heat flow in the y-direction are the highest. The difference given is calculated as The positive values of the center curve, therefore, indicate a higher calculated temperature for the two-dimensional model and vice versa at the edge.     This result fits well to the expected behavior, where the gradients in the three-dimensional model cause heat flow from the center to the edge in the y-direction. The drop in differences at x ¼ 0 is most likely caused by the application of the boundary conditions, but of little significance for this study. Finally, the most notable observations in Figure  14 are the absolute values of the calculated temperature differences between À0.02 and 0.03 K. Their magnitude is about 0:1% of the absolute temperature, indicating that the heat flow in this direction is rather insignificant for the utilized parameters, despite the strong temperature gradients. Therefore, it may be concluded that for the considered model and application, the two-dimensional model captures the temperature distribution sufficiently well.

Interpretation of results
In order to better understand the reason for the small influence of the heat release path in the y-direction in this specific case, the concept of P eclet number is employed. An important factor associated with choosing the modeling dimensions in automated layup is the ratio of advective to diffusive heat transfer, expressed by the dimensionless P eclet number [22,27] It consists of the movement velocity v, a characteristic length l that is usually the discretization element size and the material's thermal diffusivity j, all of which are specified for the one spatial direction where the advection occurs. For a P eclet number greater unity, it was assumed that thermal modeling is dominated by the advective component [27,28], even leading some authors to neglect conductive heat transfer in the direction of movement for simplicity [22,29,30].
For the presented discretization and parameters, Pe ¼ 7 is obtained and indicates that the problem is dominated by the advective component, even for the very slow selected speed. Nonetheless, this value only describes the ratio in the x-direction. In the numerical implementation of the model, these parameters forming the P eclet number are incorporated in the influence matrix K, together with the values for the other spatial directions. The temperature of an element on the inside of the modeling domain without exposure to boundary conditions is only determined by its six adjacent neighbors. Their dimensionless influence factors are expressed by thê K matrix and the sum over all elements inK is always one for a steady flow. For the selected parameters the utilized 3 Â 3 Â 3K matrix with the notation of its components iŝ number. As calculation ofK can be performed for any material and discretization parameters without any actual numerical implementation and calculation, it is well suited to help in the selection of the modeling dimensions beforehand, thus extending the concept of P eclet number's previously presented application to determine the dominant factor in the problem.

Conclusion
In this study, numerical modeling for a new LEDbased heating unit for AFP was presented. The differences to existing radiative heat sources and their modeling were discussed, leading to an optical model that is able to account for the new LED heating's ability to actively control irradiance distributions. Three different input setups were presented achieving significantly differing irradiance distributions, aiming at both process optimization and investigation of numerical modeling conditions. These were then applied in both a twoand a three-dimensional thermal model. The thermal simulation results indicate that both absolute temperatures and temperature gradients inside the placement path may be controlled via the utilized input setups, explaining the experimental observations by Orth [10]. Nonetheless, to conclude the development and understanding of the new LED heating unit, experimental validation of the presented modeling approach is required and will be published soon. From the comparison of the two-and threedimensional thermal modeling approaches, it was concluded that for the applied set of parameters, both achieve similar results even in different irradiation situations, which was attributed to the low thermal diffusity in the planar directions. Based on these findings, a tool was developed to help in the selection of modeling dimensions for heat transfer problems in automated layup, regardless of the applied heating unit.

Disclosure statement
No potential conflict of interest was reported by the authors.

Notes
1. The quantity irradiance E describes the total incident energy, but for the remainder of this work the term irradiance instead of the heat flux is used to describe the absorbed portion, in order to better represent the radiative nature of the power input and as they are linked via a linear relationship due to the assumption of a fixed absorptivity 2 .