Experimental and numerical investigations on heat transfer in fused filament fabrication 3D-printed specimens

A good understanding of the heat transfer in fused filament fabrication is crucial for an accurate stress prediction and subsequently for repetitive, high quality printing. This work focuses on two challenges that have been presented when it comes to the accuracy and efficiency in simulating the heat transfer in the fused filament fabrication process. With the prospect of choosing correct thermal boundary conditions expressing the natural convection between printed material and its environment, values for the convective heat transfer coefficient and ambient temperature were calibrated through numerical data fitting of experimental thermal measurements. Furthermore, modeling simplifications were proposed for an efficient numerical discretization of infill structures. Samples were printed with varying infill characteristics, such as varying air void size, infill densities and infill patterns. Thermal measurements were performed to investigate the role of these parameters on the heat transfer and based on these observations, possible modeling simplifications were studied in the numerical simulations.


Introduction
Fused filament fabrication (FFF), or fused deposition modeling, is one of the most widely applied additive manufacturing methods. It is a three-dimensional (3D) printing method in which a thermoplastic material is extruded through a nozzle to construct a layer-by-layer structure [1]. Once lauded for its potential in prototype manufacturing, it has now found its way into a plethora of fields and applications, such as biomedical, electrical and aerospace engineering [2]. Advantages of FFF printing are the relatively low costs, wide applicability and availability, large variety of suitable materials and the ease of use [3,4]. Some of these advantages have caused a surge of the use of FFF in four-dimensional (4D) printing, where the fourth dimension represents the change of shape over time of the smart materials [5,6]. By combining FFF printing, and the use of shape memory polymers (SMP), structures can be printed that can maintain a temporary shape and that can return to their original shape after being exposed to an external stimulus, such as heat [7]. The thermal sensitive nature of printed SMP parts enables potential applications such as fasteners in active assembly/disassembly, smart actuators and deployable structures for aerospace applications [3,8,9].
However, deposition of the filament at high temperatures, followed by quick cooling to enforce solidification results in significant thermal gradients and subsequently, residual stresses [1]. The presence of residual stresses can lead to warping during printing which can result in a failed print [10]. This is a common problem for Acrylonitrile Butadiene Styrene (ABS). Residual stresses can also cause distortion of parts and a loss of strength [11]. In case of printed parts with strict geometrical tolerances or structural requirements, this effect can lead to a loss of functionality.
During FFF printing the material is extruded at a temperature which is typically higher than the glass transition temperature. It is either deposited on the building platform or cooled existing layers, causing the material to be (partially) stretched as it bonds, to cool down and to finally solidify [12]. The induced pre-strain will be released as soon as the material is heated above the glass transition temperature and a 'new' permanent shape emerges. An example of such programming can be seen in figure 1 where two rectangular specimens are printed flat, but with a different orientation of the printed filament. After reheating the samples above the glass transition temperature, they either bend upwards, or twist upwards depending on the printing orientation. Such approach has been used in 4D printing to design structures with self-folding, self-bending, self-twisting and shape-shifting mechanisms [12][13][14][15].
Whether the objective is to repetitively print parts of high quality, or to print SMP parts with a 4D effect, understanding the heat transfer is crucial for an accurate stress and bond strength prediction [16]. Several efforts have been made to predict thermal gradients and the development of the residual stresses in printed parts by performing thermo-mechanical finite element simulations. Such simulations often make use of sequential element activation to represent the deposition sequence that takes place during FFF printing. Zhang and Chou [17] used such a model to perform a parametric study to predict part distortions in ABS. The analysis was thermo-mechanically coupled and the influence of various process parameters on the residual stresses was studied. Zhou et al. [18] presented a numerical model in which the heat transfer in ABS material was analyzed during the FDM process. In their model temperaturedependent material properties were used for the specific heat capacity and thermal conductivity. Cattenone et al. [1] performed thermo-mechanical simulations in which the FFF process was simulated to predict part distortions in ABS. Special attention was paid to the constitutive modeling of the polymer, and the influence of numerical considerations on the predicted residual stresses, such as time step size and meshing strategies. Yin et al. [16] performed similar heat transfer simulations but with a different objective. They used such models to predict the inter-facial bonding strength between printed filaments.
This work focuses on two challenges that have been presented when it comes to the accuracy and efficiency in simulating the heat transfer in the FFF process. The first challenge that arises is the correct choice of thermal boundary conditions, particularly the convective heat transfer coefficient. The exact magnitude of this parameter is not always explicitly stated, its empirical determination can be quite cumbersome and often the focus is on forced convection. Pereira et al. [19] focused on forced convection in their investigation of the effect of surface roughness on the convective heat transfer on the surfaces of FFF printed ABS cylindrical specimens. Zhou et al. [20] also calculated a value of the heat transfer coefficient based on forced convection over a rectangular body model. Costa et al. [21] did focus on natural convection between printed filament and the ambient air, but the determined values for the heat transfer coefficient still varied between a fairly large range (5-60 W/m 2 ·K) and this range was not experimentally validated. Lepoivre et al. [22] determined the heat transfer coefficient by an empirical correlation for an external free convection flow. The magnitude of the parameters used in said correlation were not listed or further elaborated.
In this work the heat transfer coefficient is determined experimentally. Experimental thermal measurements are numerically simulated and a value for the heat transfer coefficient that describes natural convection is determined through data fitting.
The second challenge is related to the computational effort which is required for an accurate heat transfer simulation. In previous work, the printed material is often modeled as a continuum in which the finite element size is determined by the filament cross-sectional dimensions. This is a computationally expensive exercise as the cross-section of a filament is much smaller than the global dimensions of printed geometries. As a result, the question arises how the modeling can be simplified to speed up simulations. The first question that must be answered is whether the assumption of discretizing the material as a continuum without accounting for the inherent air voids is an accurate one. If this is the case, methods to simplify the discretization of the characteristic mesostructure can be investigated and the computational effort

Experimental set-up
All specimens used in this work were printed with a Prusa i3 MK3 printer. The material used was Poly Lactid Acid (PLA) and the thermal properties as provided by the manufacturer Fillamentum are listed in table 1. A JADE CW infrared camera from Cedip Infrared Systems was used for the thermal measurements. The experimental set-up used in this paper is shown in figure  2. The measuring procedure was as follows. First the printing bed was heated up to a nominal temperature of 60 • C. A printed specimen was then placed on the heated printing bed for a sufficiently long time to allow for time to reach a steady-state. The heating of the specimens was recorded with the thermal   camera and the temperature profiles were recorded with a sampling frequency of 1 Hz. The temperature evolution in time was recorded for five points on the top surface of the specimens (figure 2). All measurements were performed three times and they were spaced at least 4 hours apart to ensure full cooling prior to reheating. Even though the temperature of the printing bed was set to 60 • C, the temperature which was measured on the surface of the bed with the thermal camera was equal to 56 • C. The temperature of the ambient air in the room was 25 • C.

Geometries and infill structure
To facilitate the study of the two objectives in this paper different specimens printed with varying global and infill geometries were used during the thermal measurements. A reference specimen was printed which was used in the fitting process of the convective thermal boundary condition. All of the other printed specimens had a different infill geometry (air content, infill pattern) to investigate the influence of different infills on heat flow through the specimen. There are two different ways to vary the air content in the FFF printed specimens: 1. By variation of the extrusion factor 2. By variation of the infill density Both approaches are shown schematically in figure 3. Varying the extrusion factor in a densely packed setting of filaments influences the size of the air voids between those filaments. The infill density controls the gap size between the printed filaments.
The Prusa slicer recommends a default value for the extrusion factor. This will be referred to as the default extrusion factor of 1.0 (EF=1.0). One sample was printed with the default extrusion factor, and the second sample was printed with an extrusion factor which was 10% higher than the default value. The dimensions of these samples were 30x30x2 mm and both specimens were printed with an infill density of 100%. These specimens were heated for five minutes. For the numerical discretization of the geometries with varying extrusion factor, a model as shown in figure 4 was used. The diamond shaped air voids were assumed based on a CT-scan of one of the printed samples ( figure  4). The air void size is governed by parameter a and it was used to express the area of the air voids A air as a fraction of the total area A tot : where w and h are the width and the height of a printed filament respectively, as can be seen in figure 4. The width of the filament equaled 0.45 mm and the layer height equaled 0.2 mm for all samples used in this work. By calculating the real volume fractions of air and filament in the printed samples, eq. 1 was solved for a. The real volume fractions were calculated as follows (eq. 2): 1. The mass (m) and total volume (V tot ) of each printed sample were measured. The total volume was calculated by measuring the length, width and thickness of the printed specimens with a vernier caliper. 2. From the mass, the (real) extruded volume V extr was calculated. The density of PLA was used (table 1). 3. The volume fractions of the filament v f r; f and air v f r; a were then calculated.
The samples for which the infill density and pattern were varied were printed with an extrusion factor of 1.0 and an infill density ≤100%. These samples were printed in the shape of a block with the dimensions of 30x30x20 mm. Two infill patterns were chosen; a rectilinear and a gyroid infill pattern (figure 5). The total steady-state heating time for these specimens was 40 minutes.
An overview of all the FFF printed samples with their dimensions, infill characteristics and heating time is given in table 2. The additional parameters required to describe the geometry of the cross-section of the printed specimens with varying air void size are listed in table 3.       Generally, conduction is well described by utilizing conductivity parameters found either in literature or provided by filament manufacturers. Filament cooling due to radiative heat exchange between the filaments is negligible [21]. Radiation with the environment can have an influence on filament temperature when the value of the heat transfer coefficient, the parameter which expresses convection, is relatively low (5 W/m 2 ·K) [21]. However, in most practical applications this value is much higher, which means that overall filament cooling becomes convection controlled [21]. Thus, in this work heat transfer by radiation was neglected and the focus was on determining the heat transfer coefficient and the temperature of the ambient air that describes the convection with the environment. The temperature field T (x, t) is described by the heat equation: where c p [J/kg·K] is the specific heat capacity, ρ [kg/m 3 ] is the material density, K 0 [W/m·K] is the conductivity of the material, and q [W/m 3 ] is the internal heat source. The thermal properties of PLA shown in table 1 were used here. The initial temperature of the specimens equaled the room temperature T a , thus the initial condition is expressed as: where Ω represents the domain of the printed specimen. Distinction is made between the boundary conditions at the interface between the heated printing bed and the sample Γ b , and at the free surfaces of the sample Γ f . The temperature of the heated printing bed is applied as a Dirichlet boundary condition: 5) At the free surfaces, it is assumed that the heat exchange between the sample and the environment is governed by convection. The Neumann boundary conditions are expressed as: where h [W/m 2 ·K] is the heat transfer coefficient, and T c is the temperature of the ambient air. These two parameters were calibrated by fitting the numerical simulations to the experimental data acquired in the thermal measurements.
Due to the nature of the thermal measurements, it was not entirely clear how the temperature was distributed in the vicinity of the outer surface of the sample. Since the printed specimens were relatively small, the heated printing bed might have influenced the temperature of the air surrounding the samples. This effect was taken into account during the calibration of the heat transfer coefficient. The process was structured in such a way that three situations were taken into account for the temperature of the ambient air: 1. Case 1 It was assumed that the printing bed did not have a significant influence on the temperature of the ambient air. The temperature of the ambient air equaled room temperature T a .

Case 2
The temperature of the ambient air was higher than the room temperature due to the heated printing bed. The same elevated temperature was assumed for all of the free surfaces of the heated specimen.

Case 3
In this case, it was still assumed that the free surfaces on the sides of the specimen, Γ f ; s were exposed to air of which the temperature was higher than room temperature. However, since the top surface of the specimen, Γ f ; t was more distant from the printing bed, a different temperature was assumed for this surface.

Numerical set-up
All finite element (FE) analyses simulating the heat transfer in the printed specimens were set up in Ansys (Mechanical APDL 19.2). The analysis type was a transient thermal analysis. The initial and boundary conditions are listed in equations 4-6. The FE samples were discretized with 8-node thermal finite elements (SOLID70). The default element size coincided with the width of the extruded filament and layer height of the printed specimens. This is the case unless it is explicitly stated that a different mesh discretization was used. Convergence analyses were performed prior to all simulations presented Fig. 6: Summary of experimental and numerical methodology used in this paper in this paper to ensure that the chosen mesh sizes were appropriate. As the air content in the printed samples was included in some simulations, distinction was made between the material properties of the air and of the PLA. Both materials were discretized with the same element type. All analyses that were performed in this work were assumed to be physically linear. This assumption is deemed valid as the maximum temperature that the specimens were heated up to is lower than the glass transition temperature of PLA. Thus, the constant thermal properties as given in table 1 were used in the numerical simulations.
In the simulations, the measured value of 56 • C was used for the temperature of the heated printing bed, unless it is explicitly stated otherwise. The overall experimental and numerical methodology is summarized in figure 6. On the left side of the chart it is shown that the influence of the infill density (and pattern) is measured experimentally in specimens S1 -S5. Specimens S6 and S7 are used for the measurement of the influence of the air voids. Prior to numerical validation of these two factors, the convective thermal parameters h and T c are numerically calibrated by using the experimental measurements on specimen S1, and numerically validated by comparing numerical and experimental results on specimens S2 and S3.

Varying infill density
Specimens S1, S2 and S3 were heated to capture the influence of the infill density on the heat flow. Figure 7 shows the measured average temperature and envelope for each specimen as well as a comparison of the measured average temperature between the specimens with different infill density. The envelope data consists of three measurement repetitions on five data points for each measurement (15 data points total). It can be seen that the infill density affects the steady-state temperature and the heating rate. The significant increase in   Figure 8 shows the results of the thermal measurements on the specimens S6 and S7 with varying air void size. Both the average values and the envelopes are plotted. Additionally, a comparison is made between the results of samples S6 and S7 to see what the influence is of the air void content. It can be seen that there is no significant difference between the temperatures measured on the samples printed with different extrusion factors. Varying the extrusion factor and thus the air void size does not seem to influence the heat transfer in the printed samples.

Varying infill pattern
To see whether the type of infill pattern influences the heat transfer, a comparison is made between samples printed with a rectilinear infill pattern and samples with a gyroid infill pattern. The comparison is done at an infill density of 50% (S2 and S4) and 25% (S3 and S5). The experimental results of the thermal measurements are shown in figure 9. At both infill densities, the differences between the measured temperatures of the rectilinear and gyroid specimen are negligible. Both the heating rate and steady-state temperature are similar for the samples with a rectilinear and gyroid infill pattern.  The three cases as described in section 2.3 were used in the calibration of the convective parameters, h and T c . The numerical results presented below were fitted to the the experimental results from the thermal measurements on sample S1. The experimental thermal measurements of samples S2 and S3 were used for validation of the determined parameters.

Case 1
The first situation considered in the fitting process was the one in which the ambient temperature was assumed to be equal to the room temperature of 25 • C. This value was also used as the initial temperature of the FE sample in the numerical simulations. The heat transfer coefficient in the simulations was varied between 10 and 30 W/m 2 ·K. The results are plotted in figure 10 for h=10, 20, 30 W/m 2 ·K. It is obvious that none of the simulations show agreement with the measurements. In each of the simulations, the converged temperature is significantly lower than the steady-state temperature found in the measurement (within the time frame of 40 minutes). Increasing the heat transfer coefficient above 30 W/m 2 ·K leads to an even lower steady-state temperature and a faster heating rate. For h=10 W/m 2 ·K, the steady-state is not reached within the considered time span, so using a smaller value for the heat transfer coefficient will delay this even further.

Case 2
The assumption of the printed specimen being subjected to room temperature convection in the vicinity of its free surfaces seems to be incorrect. The simulations were therefore repeated for the same values of h as in figure 10,    Fig. 11: Numerical simulations compared to experiments (case 2) but in combination with higher values for the ambient temperature. In figure  11 the results are plotted for the simulations in which the ambient temperature was varied between 25-32 • C, with heat transfer coeffients of 20 W/m 2 ·K and 30 W/m 2 ·K. The same convective parameters were applied at all free surfaces. It can be seen that the correct steady-state temperature is reached for h=20 W/m 2 ·K and T c =32 • C. However, the temperature at the start of the simulation increases much faster than in the experimental measurements. This problem is not solved by choosing a different value for h than those displayed in figure 11.

Case 3
In an attempt to tackle the overestimated heating rate at the start of the simulations, the ambient air temperature at the top surface T c; t , was varied between 25 • C and 32 • C. The temperature at the side surfaces of the specimens T c; s , was chosen to be higher than room temperature, again under the assumption that the printing bed heated up the air above it too. This temperature was varied between 32 • C and the bed temperature. Heat transfer coefficients of 20, 25 and 30 W/m 2 ·K were prescribed at all free surfaces. In figure 12 the results are plotted for the case in which the temperature at the side surfaces is equal to the bed temperature and the temperature at the top surface equals 27 • C. It can be seen that these values for the air temperature combined with a heat transfer coefficient of 25 W/m 2 ·K provide a good match with the experimental measurements. Typically, the exact temperature on the surface of the printing bed is not measured but assumed to be equal to the value prescribed in the printing settings. To account for this, the data fitting process was repeated for T b = T c; s =60 • C. These results are shown in figure 12. It can be seen that h=30 W/m 2 ·K combined with T c; s =60 • C and T c; t = 27 • C also provides a good fit with the experimental results. In the remainder of this paper, the following convective parameters are used: h=25 W/m 2 ·K, T c; s =56 • C, T c; t =27 • C.

Validation h & T c
The accuracy of the calibrated convective parameters was validated by simulating the thermal measurements performed on the specimens printed with 50% and 25% infill density. The results of the numerical simulations are plotted against the experiments in figure 13 for both of the specimens. Good agreement is found between the simulations and experiments.

.1 Air voids
In the simulations of the samples with an infill density of 100% and an extrusion factor of 1.0 presented in section 4.1, the FE meshes were made completely dense, ignoring the air voids that were present between the printed filaments. At the same time, the experiments performed on the samples with 100% infill density and a varying extrusion factor (section 3.2), implied that role of the air voids was negligible. In this section, we investigate whether a simplified material discretization as used in aforementioned simulations is accurate, or if it is necessary to model the air voids as well.
The simulations presented in this section were performed with a FE mesh which didn't include air voids, and with FE meshes which included air voids according to the extrusion factors that were used in the printed specimens (table 3). A sample of the detailed mesh including air voids is shown in figure  4. It consists of 13 hexahedral and prismatic elements for the cross-section of a single printed filament with surrounding air compared to one brick element when air voids are not included in the discretization.
The numerical results are shown in figure 14. First, the experimental results for the specimens printed with different extrusion factors are compared with the simulations obtained with the meshes including air voids. It can be seen that there is good agreement between the experimental and numerical results. The difference between the steady-state temperatures found in the simulations and experiments is less than 2%. A comparison is also made between the numerical results for the meshes in which air voids of varying size are included with a mesh in which the material is discretized as a continuum without air voids ( figure 15). No significant difference is found between the measured temperatures for the different meshes. This coincides with the observation found   in the experimental measurements of the air voids size being of insignificant influence on the heat flow. Lastly, the meshes without air voids are coarsened to a mesh where the element dimensions are twice as large as the filament cross-sectional dimensions (mesh 2), and to a mesh where the element dimensions are five times as large as the filament dimensions (mesh 3). As can be seen in figure 15, no significant difference is found between the results obtained with these three meshes.

Infill density and pattern
From section 3.1 we know that varying the infill density has a significant influence on the heating rate and steady-state temperature (figure 7), whereas the   9). This means that it might be possible to simulate the accurate heat transfer within complex infill geometries by using much coarser and simplified FE meshes. To test this hypothesis, numerical models were set-up which simulated the heat transfer in the samples printed with a 25% rectilinear infill.
In the models the geometry was discretized with the correct infill density and type, but with meshes of varying coarseness. In the default mesh the element dimensions equaled the dimensions of the filament cross-section. Based on the default mesh, two other, coarser FE meshes were used. In mesh 2, the element dimensions were twice as large as the filament cross-sectional dimensions. In the coarsest mesh (mesh 3) the element dimensions were five times as large as those of the filament cross-section. Since the empty spaces between the printed material were modeled with air elements, the air elements were also coarsened accordingly. The three different meshes are shown in figure 16.
The results are plotted in figure 17. It can be seen that all of the aforementioned mesh discretizations yield a similar temperature profile. By using the simplified and coarse rectilinear mesh (mesh 3), the computational time was reduced significantly. The number of elements and the CPU time required to solve the numerical models with the various meshes are listed in table 4. When compared to the experimental results of the thermal measurements performed on the printed samples with rectilinear infill, there is also a good agreement.
The next modeling simplification was to discretize a complex infill pattern, with both a simplified pattern and a coarser mesh. This has been done by comparing the experimental results of the specimen with a gyroid infill pattern with the coarsest rectilinear finite element model, both with a 25% infill density. These results are shown in figure 18. Again, there is good agreement between the experiments and the simulations. Both results in figures 17 and 18 support the experimental findings, namely that exact infill pattern does not influence the temperature profiles, as long as the correct infill density is taken into account.

Conclusions & Outlook
The overall objective of this work was to provide more clarity on two aspects which affect the accuracy and efficiency of heat transfer simulations of FFF printing. On one hand, a closer look has been taken at the prescription of thermal convective boundary conditions. A value for the heat transfer coefficient and ambient temperature, expressing the natural convection between FFF  On the other hand, simplifications for modeling air-filled and complex infill structures have been investigated. It has been found that accurate heat transfer can be simulated in such structures, as long as the infill density is respected. Discretization of the exact infill geometry does not lead to a significant increase in accuracy of the heat transfer simulations when compared to the experimental data. In densely packed geometries, the printed material can be modeled as a continuum. Discretization of the air voids between the printed filament is not necessary for an accurate heat transfer prediction.
The work performed in this paper was limited to the use of PLA and a Prusa i3 MK3 printer. For practical applications it would be of interest to expand the research to other materials such as ABS and to industrial 3D printers. It  must also be noted that the influence of the surface roughness of the printed samples on the convective heat transfer was not taken into account. Lastly, the modeling simplifications proposed in this work have only been applied in thermal simulations. The next step is to model the extrusion process in FFF printing and to investigate the influence of the process parameters on the heat transfer and residual stresses.