T4F3: temperature for fused filament fabrication

Temperature fields and their variations in printed parts are the basis for understanding the physical process of fused filament fabrication (FFF). However, reliable temperature data are still rather limited to date. This article presents a three-dimensional transient-state model to simulate the temporal and spatial temperature variations in FFF printed parts. Model variables range from geometry dimensions and (dynamic) material properties to process parameters, covering all important physical phenomena, including conduction anisotropy and radiant heat transfer. The validation of the model is performed against six sets of experimental temperature data obtained with different geometries, machines, materials, processes, temperature measuring methods, etc. Insights in the thermal process are also reported. For example, the heat penetration depth in printing with poly(lactic acid) is limited to 3 mm, and the Biot number intimately characterises the reheating peaks in temporal profiles. This model shows the potential to become a standardised tool to study the thermal characteristics of FFF printed parts. It is made openly available on website https://iiw.kuleuven.be/onderzoek/aml/technologyoffer.


Introduction
Fused filament fabrication (FFF) is a popular thermally driven additive manufacturing (AM) technique based on material extrusion. In a typical process, a thermoplastic filament experiences transitions from the glass state to polymer melt in the hot-end. Then, the material is pushed through the nozzle under pressure by the filament itself. The pressure transfer is enabled by the presence of a cooler polymer recirculation region located nearby the neck of the heat breaker, sealing the gap between the filament and the hot-end barrel [1][2][3]. Afterwards, the material is selectively deposited onto the build plate, and the extrudate solidifies back into the vitreous state as printed parts with cooling. The solidification process primarily determines the bond formation, bond quality and the ultimate integrity of the parts. As a measure of thermal energy, temperature characterises the non-isothermal solidification process [4][5][6]. As such, insights into the spatial and temporal temperature in the printed parts are essential to understand this AM technology, to tailor the process for optimised material performance.
On the one hand, online temperature monitoring is possible with both contacting (e.g. thermocouple [6,7], fibre Bragg grating sensor [8]) and non-contacting techniques (e.g. infrared thermography [9,10]). However, each technique has its limitations, and they can provide only partial information with limited accuracy [11]. On the other hand, numerical modelling provides a fast, economical yet powerful means to probe the temperature variations in FFF parts. The fundamental strategy is to view the additive deposition process as a sequential union of voxel elements for the geometry modelling and then describe pertinent physical phenomena based on the changing geometry [12][13][14]. As reviewed in Das et al. [15], 28 models were developed to simulate the FFF thermal process between 2015 and 2020. Most of the models for printed parts rely on the element activation function from commercial finite element software to account for the voxel element deposition. Alternatively, the finite difference method (FDM) with adjusted boundaries has also been developed [16,17]. For most of these models, a single voxel element usually represents the whole cross-section of a deposited strand (track/bead/road), and it is assumed uniform in temperature. Such an assumption is supported by the typical Biot number Bi < 0.1, indicative of a negligible temperature gradient within the strand when the layer thickness is small. Nevertheless, a multi-scaled crosssection and/or bulk geometry model was also discussed [12,16,18] and can provide information in multi-resolution at the expense of calculation time and data storage.
Apart from the geometry modelling, the physical process considered reveals more about the nature of this AM technique. Despite some inspiring results from the above models, not all critical physical phenomena were considered. For instance, the air temperature above the plate (near environment temperature, associated with the convective heat transfer) was never distinguished from the far environment temperature (room/chamber temperature, associated with the radiant heat transfer) [19]. Besides, all currently available models await sufficient experimental validations before being applied to simulate temperature development in printing different 3D geometries with different materials or process parameters. To the best of the authors' knowledge, a comprehensive and sufficiently validated model for temperature distributions and variations in FFF printed parts is still missing.
This article bridges such a gap, offering a validated and open-access numerical model for temperature fields and variations in FFF printed parts. Details of the heat transfer inside the hot-end [2,20] and during the typical 90°-turn deformation of the extrudate [21] are not included. The validation was performed concerning spatial and temporal temperature profiles in FFF with different geometries, machines, materials, processes, temperature measuring methods, etc. Because of the similarities in the fabrication process and heat transfers after the fused deposition, this model can also be applied in a few other thermal energydriven AM techniques, such as the pellet-based hot-melt extrusion (HME) [22], big area additive manufacturing (BAAM) [23], molten soda-lime glass [24,25] and molten sugar [26] AM.

Hypotheses
The heat transfer in FFF printed parts is taken as an initial and boundary value problem on a changing domain. Four hypotheses on the physical processes are adopted to model the temperature fields and their variations: H1. The FFF process is assumed a sequential deposition of ideal cuboid voxel elements. H2. The printed part after the fused deposition is stationary with reference to the build plate. H3. The volumetric heat generation by the polymer melt crystallisation is negligible. H4. The only significant radiant heat transfer is between the printed parts and the far environment.
H1 assumes FFF printed parts to be compact and porosity-free when the target infill density is 100%. Accordingly, the cross-section of the printed part is idealised as a union of tightly connected rectangles (Fig. 1, top-right). The reality contradicts such idealisation-the real cross-section contains almost unavoidable voids since the building elements, typically of elliptical cross-section, cannot compactly fill in the whole space. On the one hand, those trapped voids will not significantly contribute to the local thermal convection nor thermal radiation between different polymer strands [27]. On the other hand, they can cause anisotropy in thermal conduction due to the uneven bond length, based on which, the (extrinsic) thermal contact resistance can be defined [28,29]. Ultimately, the voids affect temperature variations in a macroscopic sense. Hence, deviations from the ideal geometry in H1 can be compensated in the modelling by introducing uncertainties in the heat transfer coefficients and anisotropy in the thermal conductivity ( x , y , z ). As an example, the conductive heat flow Q z in the model (Fig. 1) can be balanced by defining a reduced thermal conductivity according to the ratio between the real and assumed bond length, i.e. z ≡ Δy * ∕Δy . For demonstration purpose, this article only uses homogeneous thermal conductivity ( x = y = z = ), but the model developed supports anisotropy in the simulation. The authors intend to discuss their impacts and relationship with thermal contact resistance in a separate article.
H2 hypothesises FFF printed parts to be perfectly stationary, with no change in the volume/shape/velocity after the deposition. In particular, the phenomena of bond neck Fig. 1 Top: cross-section of a printed double-wall and its geometry idealisation. The heat transfer is equivalent when the geometrical factors are compensated through the heat flux by adjusted process coefficients and material properties. Bottom: a demonstration of the heat transfer mechanisms of element J20 (deposition sequence: ⋯ , J19, P19, J20, P20, J21, ⋯ ) growth are neglected considering the limited time duration and volume change. In many cases, the neck growth was not well observed in a dimensional scale comparable to the layer thickness in experiments, e.g. in Fig. 1, the sharp edges of the void kept their shape. Thus, the temperature modelling (energy conservation) decouples from the mass and momentum conservations. In the meantime, the work done by the pressure/viscous/gravitational forces is neglected.
H3 is irrelevant for amorphous polymers; for semi-crystalline polymers, it holds in most cases. For example, neither poly(lactic acid) (PLA) nor poly(ether-ketone-ketone) (PEKK TI-60/40) can considerably crystallise at a cooling rate above 10 °C/min [30,31]. In contrast, the typical (instantaneous or average) cooling rate in FFF is at least one order of magnitude higher [9,16]. In addition, PLA has a material melt-crystallisation half time t 1∕2 of about 150 s at the melt-crystallisation temperature T mc = 100 °C [32], and the t 1∕2 of PEKK is 12 s in an isotherm at 305 °C [33]. But the duration for those polymers to remain hot enough (e.g. ≥ T mc ) in FFF is usually considerably lower than t 1∕2 . However, one needs to treat H3 with discretion when a heated chamber or a high layer thickness (e.g. ≥ 4 mm, as in BAAM) is used. In either case, the overall cooling rates can be significantly reduced.
H4 simplifies the radiant heat transfer in FFF. The thermal radiation between different locations in the printed part ( Fig. 1) is negligible due to their low significance [27]. Meanwhile, the radiant heat transfer between the hot-end and printed part is also ignored for now (despite that Thézé et al. [34] suggested otherwise), since it can only contribute to a reheating in the surface temperature by no more than 5 °C [34,35]. When there are no additional external/internal heat sources for pre/post-heating [36][37][38][39], the only mechanism retained is the radiant heat transfer between the printed part and the far environment.

Heat transfer equation
Under hypotheses H1, H2 and H3, the governing partial differential equation (PDE) for the temperature field T = T(r, t) in FFF printed parts takes the form where the mass density and specific heat capacity c p can be a pressure/temperature-dependent form such that = (p, T) and c p = c p (p, T) ; t denotes the time, the thermal conductivity, r = (x,y,z) the spatial location, ∇T = T x , T y , T z the temperature gradient. The changing domain Ω(t) is specified by the three-dimensional geometry and the building strategy in FFF. Hypothesis H4 specifies the boundary condition (BC) to the combined type (convection + radiation). Across the boundary, the directional heat flux is

Geometry modelling
The basic geometry of interest is a cuboid of length L ( x direction), width W ( y direction) and height H ( z direction, Fig. 2a). Under hypothesis H1, it can be built by depositing ideal voxel elements in a zigzag pattern (a rectilinear infill pattern at 0 • and an infill density of 100%). The building element has dimensions Δx × Δy × Δz , where Δy represents the strand width (assumed proportional to the nozzle (1) (2) q = q conv + q rad . diameter such Δy = , where is the shape factor [3]); and Δz equals the layer thickness ℏ , Δz = ℏ . The product vΔyΔz provides an estimation of the volume flow rate Q v in FFF, where v is the nozzle travelling speed relative to the build plate. The choice of the element dimension Δx is arbitrary, but it determines the time step between pseudo-depositions at adjacent locations such that Δt = Δx∕v . Accordingly, element numbers in each direction are N L = L∕Δx , N W = W∕Δy and N H = H∕Δz . The total element number adds up to N = N L N W N H . The inter-layer time Δ , indicating the time required to print one layer, can be determined by Δ = N L N W Δt.
Voxel elements in FFF are chronologically deposited. Accordingly, each has a unique time index n indicating its birth-moment t 0 ( t 0 ≡ nΔt , n = 1, 2, ..., N ). The birthmoment of t he last element is denoted by t end ≡ NΔt (= N H Δ ) , which is also the time required to finish the printing. In the meanwhile, every element also has a set of unique spatial indexes (i, j, k) 1 describing its position r = (iΔx, jΔy, kΔz) , (i = 1, 2, … , N L ; j = 1, 2, … , N W ; k = 1, 2, … , N H ) . Given the specific zigzag travelling path, a bijection exists between indexes (i, j, k) and n [16]. These voxel elements are symbolically denoted by E ijk , and the FFF building process is essentially a sequential union ⋃ E ijk . For each element E ijk , its temperature is assumed uniform within its volume (typically, the Biot number Bi < 0.1) and thus denoted by T nΔt ijk (or simply T n ijk ) 2 at the discrete time nΔt . The time index n here should be higher than or at least equal to the birth-moment index uniquely determined by its spatial indexes (if not, it means that the element has not been deposited yet). n can go beyond N , indicating that FFF is finished, but the heat transfer continues. For the sake of discussion, six directions are defined with respect to the axis orientations: front-positive x ; back-negative x ; left-positive y ; right-negative y ; up-positive z ; down-negative z , as shown in Fig. 2b.
Two special geometries are introduced: double-wall ( N W = 2 ) and single-wall ( N W = 1, Fig. 2a). Their temperature variations can be effectively monitored, as demonstrated in [9,10], respectively. For both geometries, the inter-layer reheating is an important heat transfer characteristic; for the double-wall, the intra-layer reheating can also be well observed. Temperature information in these geometries is representative of printing other geometries as well. When depositing the first perimeter on a new layer in an arbitrary geometry and before any intra-layer or inter-layer reheating, the conductive heat transfer takes place only in the downwards direction in the plane perpendicular to the velocity vector (as anywhere in the single-wall or the right track of the double-wall, e.g. J20 in Fig. 1). When depositing at any locations other than the first deposited perimeter, the conductive heat transfer typically occurs in two directions (downwards and left/right), as represented by the second deposited (left) track of the double-wall (e.g. P20 in Fig. 1). Although reheatings upon thermal contact may bring up the temperature above the material glass transition temperature T g , its contribution to bond formation can be negligible [4]. Therefore, the temperature information in the doublewall can shed light on the bond formation in printing bulky geometries, as long as the intra-layer and inter-layer times are comparable for different locations of interest.

The numerical scheme
The finite difference method (FDM, [40]) was adopted to discretise the governing Eq. (1) into where = ∕ c p denotes the thermal diffusivity, and T n ijk the element temperature in position (iΔx, jΔy, kΔz) at nΔt . The complete form of T n i+1⋅⋅ is T n i+1,j,k , where the dots ' ⋅ ' indicate the positions of trivial spatial indexes. The truncation error on the right-hand side of Eq. (3) is of O Δx 2 + Δy 2 + Δz 2 . The only unknown T n+1 ijk for the next moment can be expressed as a linear combination of the temperature of 6 neighbour elements and itself (Fig. 2b) at the current moment. Hence, Eq. (3) is also referred to as a 7-point scheme. This linear operation can be represented by a tri-diagonal block matrix [16] or more conveniently by a linear operator F such that A sufficient condition for a stable 7-point scheme for elements away from the boundaries requires Since Δt = Δx∕v , setting Δx = [(Δy) −2 + (Δz) −2 ] 1∕2 can deliver a most probable stable scheme. In hundreds of simulations already performed, the scheme stability was well observed. In case of instability, the inter-layer time similarity rule [3] can be applied to simultaneously scale up the part length and travelling speed so that the time step in Eq. (5) can guarantee the inequality. Practical tips on the scheme stability are also available in [41]. The initial condition is given by the temperature of every element at its deposition as where T i is the actual polymer temperature at the nozzle outlet, and it is assumed to be equal to the nominal nozzle temperature T n in the stable filament feed region [3]. It is noteworthy that T i can greatly deviate from T n (up to 30 °C [20,42]) at high volume flow rates when the heating power of the hot-end lags behind the polymer melt flow.
The main difficulty in applying Eq. (4) in the calculation is the appearance and disappearance of the boundaries, which originates from the additive nature of FFF. Consequently, the neighbours of E ijk may not have been deposited at a specific moment or might never exist (e.g. element J20 does not have the neighbour P20 before its intra-layer reheating, Fig. 1). This consequence requires boundary condition Eq. (2) to determine the temperature of those pseudo-neighbours.

Boundary condition treatment
For an element surface having an outward normal e x = [1, 0, 0] , the BC gives a constant heat flux across the boundary, equating that by thermal conduction, where h denotes the convective heat transfer coefficient; the whole spectrum hemi-spherical material emissivity; the Stefan-Boltzmann constant (5.670 × 10 -8 W/m 2 /K 4 ); T ∞ the far environment (room/chamber) temperature; and T a the air temperature above the build/heat plate. T a depends on the vertical height z and decays exponentially from the plate temperature T p towards T ∞ by a characteristic dimension l c such that Reference [19] provides more details for Eq. (7). For a parallel surface but negatively oriented (outward normal e −x = [−1, 0, 0] ), the BC gives In the context of the FDM, Eq. (6) is discretised by a forward scheme [40] into By rearranging, Equation (9) is already an iteration-ready form for the calculation. By defining the local radiant heat transfer coef- (9) is neatly rearranged into assuming that the air temperature equals the far environment temperature, i.e. T a (z) = T ∞ 3 . Recalling the direction convention in Fig. 2b, the unknown neighbour temperature is given by Similarly, Eq. (8) is discretised by a backward scheme [40], and ultimately rearranged into (counterpart of Eq. (10)) The truncation error in Eqs. (10) and (11) is of the order O(Δx) , which exceeds that of the PDE by Eq. (3). Hence, to realise better error control, the central difference schemes are applied and deliver for Eq. (6) and for Eq. (8). By rearrangement, the counterparts for Eqs. (10) and (11) are which are also linear combinations of the temperature of the boundary element ( T n ijk ), the far environment (T ∞ ) and the neighbours in the opposite direction ( T n i∓1⋅⋅ ). The truncation error now is of the order O Δx 2 , matching that of the PDE scheme. However, Eqs. (12) and (13) require the temperatures at two physical locations to calculate the temperature of boundary elements. The calculation is not always applicable when a neighbour element has not been deposited yet (e.g. when printing the first strand of a new layer), or when such neighbour element never exists (e.g. when printing a single-wall, Fig. 2a). In those cases, the lower order discretisation by Eqs. (10)(11) is the only choice.
The BC of the first layer downwards is of the first type (Dirichlet), and the neighbour temperature is assumed to be equal to the plate temperature T p : T n i,j,0 = T p . Alternatively, a thermal contact resistance between the parts and build plate may be defined. The schemes for BC at all other surface orientations resemble those for the ±x directions (Eqs. (10) and (11) for the O(Δx) scheme, Eqs. (12) and (13) for the O Δx 2 scheme).

Pseudocode and software application
The pseudocode in Fig. 3 demonstrates the main algorithm structure of the T4F 3 model. The primary difficulty resides in step 7. One needs to judge whether an element within the boundary has been deposited and select which type of BC discretisation scheme for the calculation. Those details are further outlined in Appendix 1- Fig. 13. Executable software packages (Windows APP and Matlab APP) were also developed to allow more user-friendly applications of the model, which are openly available on https:// iiw. kuleu ven. be/ onder zoek/ aml/ techn ology offer. A screenshot of the main graphic user interface is provided in Appendix 2- Fig. 14.

Materials and methods
The T4F 3 model is first validated through a basic set of experiments (Set 1). This dataset verifies the capability of the model to predict the temporal and spatial temperature profiles of printed parts. Printing experiments were carried out and monitored via thermography. The experimental temperature data were then compared with simulations.
The geometry in Set 1 is a double-wall of dimensions (Fig. 4). It was printed on a Prusa i3 MK3, using an E3D V6 brass nozzle with outlet diameter = 400 μm ( = 1/2 W ) and the transparent Prusa PLA filament (diameter 1.75 mm). A thermographic infrared (IR) camera (Optris PI640), working in the spectral range 7.5-13 μm , was employed for the front 4 view realtime monitoring. The camera worked at a spatial resolution of 31.25 μm/pixel in a field of view (FOV) of 20.0 × 15.0 [mm]. The sampling frequency was 32 Hz. The accuracy of the IR camera is ± 2 °C or ± 2%, whichever is greater. The temperature at a given location was taken as the mean value of 6 × 3 pixels, based on an emissivity of 0.78 for PLA [9]. Table 1 summarises the FFF process parameters and pertaining materials properties. In particular, the part cooling fan was set to its maximum volumetric flow rate of 3.80 cubic Algorithm T4F 3 main structure Require: geometry, process, material properties, [i, j, k] ↔ n Ensure: N H ≥ 2, T n i,j,k known ∀ existing elements 1: for n = 2, 3, · · · , N, · · · (at each time step) do 2: if n < N (FFF in progress) then 3: for J = 1, 2, · · · , n (∀ existing elements) do 4: take T n i,j,k at the current moment (t = n∆t) 6: (update material properties if necessary) 7: (whether a neighbour element has been deposited? ) 8: new element deposition 10: else (FFF is finished)

11:
(update process parameters if necessary) 12: take T n i,j,k at the current moment (t = n∆t) 15: (update material properties if necessary) 16:

Fig. 3
The main structure of the T4F 3 algorithm Fig. 4 Schematic of the double-wall geometry in Set 1, the printing process and its optical and thermographic images feet per minute, and the convective heat transfer coefficient h was set at 60 W/m 2 /K as in [27].
In addition, five more datasets (Set 2-6) were collected to study the robustness of the model under different experimental conditions, including changes in the processes, materials, part geometries and equipment. These studies focus only on the temporal profiles at selected locations. The experimental details are summarised in Table 1. The features of each dataset are briefly introduced hereby.
The experimental work in Set 2 is constructed as Set 1 except for the nozzle travelling speed. The temporal temperature profiles at the geometry centre (front view) are captured and analysed. The roles of radiant and convective heat transfer are discussed. Set 3 makes use of an own-developed robot arm assisted FFF experimental setup [3]. A standard E3D V6 hot-end was mounted on a robot arm and controlled for movement. The printing was performed on a non-heated poly(methyl methacrylate) (PMMA) plate covered with paper tape in the open air. The geometry printed was a single-walled box. Set 4 and 5 include experimental data kindly shared by Lepoivre et al. [10], where a heated chamber was used. The heated chamber can effectively reduce the convective heat flux intensity and, ultimately, the overall cooling rate. Besides, the nozzle temperature in Set 5 exceeded 350 °C in printing PEKK. Set 6 includes experimental data kindly shared by Compton et al. [23], where the BAAM

Results
The model validation against Set 1 experimental data are discussed with respect to (1) temporal temperature profiles at given locations, (2) temperature fields at different times and, (3) temperature gradient at selected locations and times.

Temporal temperature profiles
The experimental and simulated temporal temperature profiles at different locations on the 20-th layer in the doublewall are shown in Fig. 5. Rapid cooling takes place at all three locations right after the deposition, followed by periodic intra-layer and inter-layer reheatings (denoted by ↑ and ⇑ , respectively). The intra-layer reheating manifests the thermal conduction between strands on the same layer. Given the dimensions and parameters, the intra-layer time depends on the x coordinates such = 2(L − x)∕v . At travelling speed v = 10 mm/s, the intra-layer reheating (denoted by the first ↑ from left) takes place 1.8 s after the voxel element deposition at the geometry centre ( r = (9, 0, 6) [mm], Fig. 5C). The inter-layer time Δ = 2L∕v = 3.6 s, which doubles the intra-layer time. Thus, intra-layer and inter-layer reheatings are equally spaced but differ in intensity because of different conduction conditions (e.g. bond length). On the left 5 side, e.g. at r = (2, 0, 6) [mm], is comparable to Δ ; the peaks due to intra-layer and inter-layer reheatings overlap with each other (Fig. 5L). On the right 5 side, e.g. at r = (16,0,6) [mm], ≪ Δ , the intra-layer reheating takes place right after the deposition. As a result, only a tiny peak/disturbance is observed (Fig. 5R); consequently, the first interlayer reheating peak overlaps with the secondary intra-layer reheating (Fig. 1, P20 → (J20&P19) → J19). At all locations, reheatings can penetrate through a few layers. However, the deeper they penetrate, the lower their intensities, as suggested by the greyscales of the ↑ and ⇑ arrows in Fig. 5.
From Fig. 5, the model can predict similar temperature patterns in thermography. However, considerable discrepancies are visible between the experimental monitoring and simulations. Particularly, the maximum temperature upon deposition is not correctly captured in the monitoring, possibly due to the thermal reflections in the IR signal transmission. At such a low volume flow rate of 4.3 cm 3 /h (Table 1), the extrudate temperature is expected to be effectively raised to the nominal nozzle temperature [42,46]. Thus, the majority of the discrepancy comes from the experimental error in this specific dataset. A critical discussion on the possible sources of inaccuracy (such as thermal reflections, sampling frequency) in IR thermography is given in [9]. It awaits dedicated calibrations before accurate temperature data can be extracted from IR monitoring. Furthermore, the intra-layer reheating peak at = 0.4 s is not correctly predicted by the model (Fig. 5R); but in the simulation for a next location, The diagrams below show the moment before the intra-layer reheating at each location e.g. r = (14.2, 0, 6) [mm] ( = 0.76 s), it appears. This observation indicates that the uncertainties in input variables may also contribute to the discrepancy between the experimental and simulation data.

Temperature fields
A few characteristic times were chosen to examine the capability of the model to predict the temperature fields T(x, y = 0, z, t) of the printed part. The experimental observations are shown on the left side of Fig. 6. Specifically, t end (= 144 s) denotes the total time required to print the double-wall in FFF, and 2t end is a typical time indicating the steady state after the printing.
As observed in Fig. 6-left, newly deposited layers remain relatively hot upon finishing ( t = t end ). The temperature gradient mainly exists in the vertical direction ( |∇T| ≈ | T∕ z| ) and directs downwards ( T∕ z > 0 ) for the top 8-12 layers. Recalling the layer thickness ℏ = 300 μm , the impact reaches roughly 3 mm-a characteristic heat penetration depth for PLA. On the contrary, the gradient always directs upwards ( T∕ z < 0 ) for the bottom 8-12 layers above the build plate. These gradients suggest that newly deposited layers have limited influence over the energy equilibrium of the bottom layers, when the nozzle prints the 24-th layers onwards (or the 16-th, as a bolder estimation. See white boxes in Fig. 6). In addition, the energy balance will be locally reached within their above and below 8-12 layers. Such a local equilibrium is almost independent of the printed part height and provides a physical proof of the active body concept used by Zhang and Shapiro [17] in their simulations. These observations suggest that a local thermal disturbance exerts only limited influence in space.
On the other hand, the temperature field upon finishing ( t = t end ) is not symmetric with respect to the plane x = L∕2 . Due to the sequential deposition, the corresponding thermal boundary is asymmetrical in spite of the symmetry in x coordinates. At the steady state (e.g. t = 2t end ), the temperature field from the monitoring is neither in perfect symmetry because of some subtle differences in the mesostructures at different locations [9]. These results show that the gradient component T∕ x exists in printed parts (but not obvious on a colour map of a wide range in thermography), and it sets a reminder that a simplification of the FFF thermal simulation from 3D to 2D or even 1D heat transfer problem may need further discretion.
The simulations in Fig. 6-right capture all characteristics aforementioned, as they can manifest the moving of the active body, the concept of heat penetration depth and the temperature gradients at different times. However, some discrepancies still exist. For example, the temperature fields at the steady state ( t = 2t end ) do not show equally apparent variations along the x axis, indicating the material properties and process parameters may require fine-tuning. Despite so, these results earn more credibility to the T4F 3 model.

More on the temperature gradient
It is clear from Fig. 6 that the temperature gradient component in the x direction is far less significant than that in the z direction, i.e. | T∕ x| ≪ | T∕ z| . Hence x = L∕2 is fixed, and more results on T = T(z, t) at different times are discussed. Here, t takes different values after the element deposition at r = (9, 0, 12) [mm] at t 0 = 141.30 s. Figure 7 presents corresponding results.  Fig. 7, the monitoring and simulations can capture similar evolutions of the temperature profiles. Since the deposition at t 0 , considerable temperature gradients exist in the top few layers along the z direction. The maximum gradient component | T∕ z| reached 920 °C/mm from the monitoring; while in the simulation, a maximum | T∕ z| = 568 °C/mm was predicted. Although the exact value differs, the magnitude matches. Such gradients quickly faded out as a consequence of fast cooling. Eventually, the temperature profile takes the shape of the near environment temperature T a (z).
As the temperature gradient directs the conductive heat flow, one can observe that most thermal energy flows downwards for the last 8-10 layers ( T∕ z > 0 in general). This layer number can also be used to define the heat penetration depth and agrees with the observations in Fig. 6. The heat penetration depth depends on the process, material and, more profoundly, the geometry. The thin double-wall in this dataset will cause energy loss to the environment faster than bulky geometries of smaller specific surface areas. On the bottom layers, temperature gradients had already reached a quasi-steady state at t 0 , suggesting that the penetration depth is limited. The simulations in Fig. 7b can capture the profile evolutions since the deposition. The heat penetration depth and the eventual shape of the profiles are well demonstrated, showing that the model can reveal temperature gradient information as well.

Robustness study and discussions
Another five sets of temperature data were collected (Table 1) to further validate the model and probe its robustness. These data cover different aspects in heat transfers in FFF and similar technology (e.g. BAAM), including different machines, materials, geometries, processes, volume flow rates and temperature measuring methods. Simulations were performed, analysed and compared with experimental monitoring data. They validate the robustness and the applicability of the model to a variety of typical FFF situations and demonstrate the potential of the T4F 3 model as a learning and analysing tool for FFF. In each dataset, multiple simulations were performed, and key parameters were varied for comparisons, as summarised in Table 2.

Set 2 and an example of different processes and the influence of radiation and convection
The experimental and simulation work presented in Set 1 also apply to other FFF situations. Currently, at least 18 variables (part geometry dimensions L × W × H , nozzle temperature T n , plate temperature T p , room/chamber temperature T ∞ , air temperature above the build plate T a (z) , travelling speed v , strand width , layer thickness ℏ , material density (p, T) , specific heat capacity c p (p, T) , thermal conductivity ( x , y , z ), convective coefficient h , emissivity , and flat nozzle tip width [41]) can be varied in the model while experimentally monitored.
As an example of different processes, data Set 2 was collected, including temporal temperature profiles at different nozzle travelling speeds. The monitoring and simulation details coincide with Set 1 (Table 1). Critical parameters in each subset are summarised in Table 2. In addition to that, the voxel element dimension Δx was increased in the simulations to reduce the storage and calculation time; however, its influence on the result is trivial.
The temporal temperature profiles at the geometry centre in the front view are presented (Fig. 8a). The travelling speeds v were 5 and 20 mm/s in the monitoring (Exp.2-1 and Exp.2-2). The normalised temperature  and inter-layer reheatings share similar patterns for both speeds. However, the exact intra-layer and inter-layer times differ, resulting in different reheating moments. In the faster case, both inter-layer and intra-layer time are shorter, inducing higher neighbour temperature upon reheatings. Ultimately, the temperature gradients upon reheatings are smaller, and the reheating peak heights are lower. These thermal characteristics in the monitoring are well reflected in the simulations, suggesting the model can capture the essential heat transfer mechanisms in FFF with different process parameters. However, if the travelling speed further increased to 60 mm/s, the printing showed repeatable failures due to slow cooling. The simulation Sim.2-3 intuitively explains the situation. As an indicator, the temperature upon inter-layer reheating was ~ 160 °C, comparable to the melt temperature of PLA.
Likewise, some discrepancies exist between the monitoring and simulations in Fig. 8a. Apart from the aforementioned errors in the IR camera, the remaining subsection exams the roles of radiant and convective heat transfers in the simulations based on the experimental work at v = 5 mm/s. The significance of radiant heat transfer was investigated with Sim.2-4, which was identical to Sim.2-1 (details in Table 2) except that the material emissivity was 0. Figure 8b presents corresponding results. Visual inspection shows that the temperature profile is slightly higher than Table 2 A summary of non-trivial variables in the simulations *: "-" indicates trivial variables that are identical within the dataset. They can be found in Table 1 Dataset Simulations  Sim.2-1 when the radiant heat transfer is neglected. To quantify the difference in the two temporal temperature profiles, the ||•|| ∞ norm was applied within Δ seconds after the deposition: where • denotes the difference between any two temporal profiles. By calculation, ||T 2-4 − T 2-1 || ∞ = 4.1 °C. Such a difference can be neglected in FFF, suggesting the contribution of radiant heat transfer in the cooling is negligible in this dataset. This suggestion agrees with the conclusion in [27] that the radiant heat transfer from printed parts to the environment can be ignored when forced convection is anticipated in FFF. Indeed, apart from the forced convection, the low layer thickness ( ℏ = 300 µm) also characterises this dataset. The significance of radiant heat transfer will be reexamined with other datasets. The convection intensity can be subjected to considerable variations in FFF since (1) the part cooling fan only works in the nozzle vicinity and has directional influence; (2) the air temperature above the build plate also affects the temperature difference in Newton's law of cooling, thus the convection flux. In Sim.2-1, the convective heat transfer coefficient h was set to 60 W/m 2 /K when the fan was activated and set to its maximum speed. Regarding the uncertainty in this coefficient, consider Sim.2-5 as a limiting case with h = 8.5 W/m 2 /K for natural convection [23]. Sim.2-5 differs from Sim.2-1 only by h . Logically, a lower convective coefficient results in a slower cooling profile (Fig. 8b). At the same time, all other heat transfer characteristics (e.g. inter-layer and intra-layer reheating peaks locations and heights) are preserved. The ||•|| ∞ norm calculates that ||T 2-5 − T 2-1 || ∞ = 26 °C, meaning that the accuracy in h can exert a significant influence on the cooling in FFF, as also reported in [7,27]. Hence, an accurate description of the convection coefficient h under different fan working conditions and surface orientations is vital to the model accuracy.

Set 3 and the influence of temperature-dependent material properties
Set 3 includes experimental data collected in printing with a six-axis robot arm. The geometry is a single-walled square box, equivalent to a single-wall of quadruple length in heat transfers. Details in the experimental monitoring and simulations are in line with [3]. Figure 9 presents the experimental temporal temperature profile at the geometry centre (Exp.3) with different simulations. The shifted time 0 indicates the deposition moment at the location of interest. In Sim.3-1, the mass density and specific heat capacity c p took the temperature-dependent form. In contrast, the thermal conductivity was a constant as its dependence on T is weak (Appendix 3-Fig. 15). In general, Sim.3-1 and monitoring data agree well with each other regarding the initial rapid cooling, inter-layer reheatings, overall cooling, etc., despite that some discrepancies are observed.
Although both and c p strongly depend on temperature T , their product (and ultimately the thermal diffusivity ) only has a weak dependence on T , especially above T g (Appendix 3- Fig. 15). To investigate the influence of temperature-dependent material properties, Sim.3-2 with constant * = 1.226 g/cm 3 and c * p = 1.801 J/g/K was performed and presented in Fig. 9. These material properties were chosen based on an equivalent entropy density concept at a to-bedetermined temperature T * such that With the experimental data in [44], T * was found to be 58.11 °C. Consequentially, * and c * p were determined by interpolation with the 'pchip' method. However, the maximum temperature difference between the simulations with dynamic and static material properties is which is hardly discernible in Fig. 9. On the other hand, if simply taking � = 1.25 g/cm 3 and c � p = 1.314 J/g/K at * (T * )c * p (T * ) = the room temperature 25 °C in Sim.3-3 (the product ′ c ′ p is lower than * c * p by 26%), the profile does not significantly differ from the previous two. The maximum temperature differences are Hence, taking constant material properties in simulations in this article and many other publications may be justified.
Since the thermal diffusivity = ∕( c p ) is what ultimately used in the calculation, the uncertainties in or c p --in a dynamic/static form, or at the room/elevated temperature--can also be reflected in the uncertainty in , and vice versa. Data Set 4 sheds more light on this respect.

Set 4 and the uncertainty in thermal convection and conduction
The experimental data in Set 4 come from Lepoivre et al. [10], where a heated chamber was used. The chamber temperature was raised to reduce the convective heat flux. The temporal profile at the centre on the sixth layer in printing a single-wall is shown in Fig. 10a (Exp.4). The temperature profile decreases towards the elevated chamber temperature ( T ∞ = 95 °C). During the cooling, only one inter-layer reheating peak is noticeable, suggesting that the heat penetration lasts for only two (or at most three) layers at an intermediate layer thickness ℏ of 800 µm. The simulation performed by Lepoivre et al. [10] used a different numerical scheme, where they considered the thermal contact resistance but ignored the radiant heat transfer ( = 0 ). Following the exact process parameters and material properties in [10], the T4F 3 model overpredicted the cooling (Sim.4-0, Fig. 10a); however, the reheating characteristics in the profile (e.g. the reheating peak height and penetration depth) matched the monitoring data. In particular, the authors used h = 30 W/m 2 /K for natural convection. Considering the uncertainty in the convection intensity, two levels of lower h were chosen and used in Sim.4-1 ( h = 15 W/ m 2 /K) and Sim.4-2 ( h = 8.5 W/m 2 /K). Logically, the lower the value of h , the slower is the overall cooling (Fig. 10a). The profile in Sim.4-2 was substantially higher than the monitoring, but if the radiant heat transfer was considered (additionally taking = 0.91 ), Sim.4-3 provided a better result.
As analysed, the uncertainties in material properties ( , c p ) and mesostructure dimensions (hypothesis H1) can be directly transferred to thermal conductivity . Thus, an effective for proper calculations can considerably deviate from the intrinsic material property. To investigate the impact of such uncertainty in on the temperature profile, Fig. 10b presents a series of simulations (based on Sim.4-3, 0 = 0.2 W/m/K), ranging from 1% 0 to 10 0 . For all simulations, algorithm stability was observed. In the figure, is positively correlated to the cooling rate. The higher the , the higher is the reheating peak height, and the more pronounced is the heat penetration depth (or, the bigger is the reheating peak number). Furthermore, the magnified Fig. 10c indicates that when is higher, the reheating peak width (full width at half maximum) is smaller, and intuitively, the time to reach the peak summit due to reheating is shorter. These observations are consistent with the simulations in Fig. 9, where the product c p in Sim.3-3 is lower than that in Sim.3-2 by 26%, which is equivalent to a higher in Sim.3-3 (but with the same c p ). Similarly, the reheating peak is higher, and the time required to reach the peak summit is shorter in Sim. 3-3. These correlations between reheating peak characteristics and thermal conductivity can help to double check the magnitude of material properties used in simulations. One such example is given in data Set 5.

Set 5 and the uncertainty in the emissivity
Temperature data in Set 5 come from FFF with PEKK in [10], where the nozzle temperature T n was 356 °C. The experimental temporal profile at the centre on the sixth layer in printing a single-wall is displayed in Fig. 11a (Exp.5). The maximum temperature upon deposition was not properly captured, as in Set 1.
Following the exact material properties and process parameters presented in [10] ( = 0.5 W/m/K, h = 30 W/ m 2 /K, = 0), Sim.5-0 from the T4F 3 model only gave a rough prediction, with the absolute difference reaching up to 50 °C. Recalling the correlations in data Set 4 (higher --higher reheating peak height/number--lower peak width), all reheating characteristics in Sim.5-0 indicate that = 0.5 W/m/K is excessively higher than the effective material conductivity. (The original simulation in [10] also considered the thermal contact resistance.) Hence, = 0.25 W/m/K from [43] was adopted for PEKK in all further simulations. The subsequent Sim.5-1 delivered a better prediction than Sim.5-0, with slower cooling rates and a reheating peak of lower height but larger width. Furthermore, following the same choice of h = 8.5 W/m 2 /K for natural convection, Sim.5-2 provided an upper bound for the monitoring. If the radiant heat transfer is additionally considered, Sim.5-3 ( = 0.94) delivered a lower bound for the monitoring. These results suggest that an intermediate emissivity could deliver a better simulation.
However, accurately determining the emissivity (or any other material property/process parameter) is beyond the scope of this article. Figure 11b intuitively demonstrates how the emissivity affects the temperature profiles: the higher the emissivity, the faster is the cooling and the lower the temperature approaches the steady state. Unlike the case in Set 2 where the radiant heat transfer is insignificant, radiant heat transfer in Set 5 is still essential even at a low emissivity (e.g. = 0.25) when the nozzle temperature is high.
It is worthwhile to note that Sim.5-1 and Sim.5-3 give almost identical results: ||T 5-1 − T 5-3 || ∞ = 2.5 °C. This coincident appeared because the effective heat transfer coefficient h ef f ≡ h + h rad in Sim.5-3 is higher than h = 30 W/m 2 /K in Sim.5-1 at the deposition temperature; but, it decreases with the temperature during cooling. Such a coincident suggests that the radiant heat transfer somehow acts as if an additional source of thermal convection with a coefficient h = 21.5 W/m 2 /K is present in this dataset. Hence, its significance during the cooling is self-explanatory. More results on the significance of radiant heat transfer are presented in Set 6.

Set 6 and the influence of resolution (or layer thickness)
Temperature data in Set 6 come from Compton et al. [23], where a BAAM system working on fused pellets was used in the printing. The BAAM system is characterised by the capability to build geometries of big dimensions (e.g. > 1 m 3 ) at a high material flow rate (e.g. > 10 4 cm 3 /h) and a high layer thickness (e.g. > 1000 µm. A higher layer thickness means a lower resolution). Although the BAAM system is not a fused filament technique, the processes after the material leaves the nozzle share a great similarity.
In printing a single-wall of 88 layers, the temperature profiles of a few characteristic layers were monitored and demonstrated as the noisy curves in Fig. 12 (Exp.6). It is interesting to note that the reheating peaks disappeared. This phenomenon appreciably resembles the simulations in Fig. 10b when the conductivity took only 25% or 1% of the true material properties ( = 25% 0 , the Biot number Bi = hℏ∕ = 0.14; = 1% 0 , Bi = 3.4). However, the disappearance here is due to the layer thickness ℏ reaching 4 mm: on the one hand, ℏ is comparable to the typical heat penetration depth in FFF, meaning that thermal energy cannot efficiently penetrate through such a distance; on the other, the temperature is no longer homogeneous across the layer thickness in the millimetre range. The Bi = 0.20 indicates that the energy transferred by conduction across adjacent layers may lose its significance to the combined effects of thermal diffusion within the strand, convection and radiation. As a result, one can only observe discontinuities in the cooling rates before and after the reheating but no longer the reheating peaks.
The Sim.6-0 took the exact material properties and process parameters as in [23]. Despite some discrepancies, all characteristics (e.g. disappearance of the reheating peak) in the temperature profiles were well simulated. Except for the first layer (which was heavily influenced by the build plate), Sim.6-0 overpredicted the cooling with a material emissivity = 0.87. However, the predictions by Sim.6-1 gave an upper bound for the experimental monitoring when ignoring the radiant heat transfer. Logically, an intermediate emissivity can give a better prediction. In addition, its significance is self-explanatory again, even at a low-moderate nozzle temperature of 200 °C.

Conclusions
This article presents a three-dimensional transient-state numerical model T4F 3 for temperature fields and their variations in fused filament fabrication (FFF) printed parts. It admits at least 18 input variables among geometry dimensions, material properties and process parameters, including dynamic material properties, material anisotropy and radiant heat transfer. The model validation is performed against six sets of experimental data, covering data obtained with different machines, geometries, materials, processes, temperature measuring methods, etc. Certain discrepancies between the experimental data and simulations are observed, but they can be reasonably ascribed to the errors in experimental monitoring and uncertainties in the input variables. All critical heat transfer characteristics in the monitoring are satisfactorily simulated. The T4F 3 model has been made openly available on the website https:// iiw. kuleu ven. be/ onder zoek/ aml/ techn ology offer. During the investigations, significant insights into the thermal process of FFF are also obtained: • A concept of heat penetration depth has been identified for FFF printing with poly(lactic acid), ranging from 2.4 to 3 mm (or 8-12 layers at a layer thickness of 300 μm ). Any local thermal disturbance (due to deposition, reheatings, unique thermal boundary, etc.) has little influence on locations 3 mm away. • The reheating peaks in temporal temperature profiles are closely related to the Biot number Bi. A lower Bi-due to lower layer thickness and/or higher thermal conductivity-correlates to a higher reheating peak, a lower peak width, a shorter time to reach the summit, and a higher heat penetration depth. • The radiant heat transfer between the printed parts and the far environment can be ignored in FFF only when printing at a medium/high-resolution (e.g. layer thickness ℏ ≤ 300 μm ) and subjected to forced convection. In all other cases, it cannot be ignored. • Temperature-dependent material properties may not outperform constants in the heat transfer. Constant material properties at a temperature close to the glass transition temperature are recommended for simulations.
So far, the T4F 3 model has been successfully applied to explain the role of a closed chamber (on the decaying air temperature T a (z) above the build plate [19]), the possibility of FFF in vacuum (on the dependence of thermal convection on the air pressure [47]), the inter-layer time similarity rule [3], etc. With minimal modifications, it can be used to quantify the effect of pre/post-heating with a hot plate attached to the nozzle [39] (on the role of a specialised form of T a (z) ), to study the hot-end designs [34,36] (on the role of radiant heat transfer from the hot-end), or to study the thermal process of continuous fibre-reinforced printing [48] (on the role of conduction anisotropy), etc. It is noteworthy that the model still suffers from heavy calculation and data storage.
The simulation process may lag behind the physical process of printing. Hence, it does not suffice for online feedback control with the current single-scale geometry modelling strategy. Despite so, it establishes a firm step towards future investigations into the non-isothermal bond quality, crystallinity development (for semi-crystalline materials only), thermal deformations, etc., where reliable temperature data are the prerequisite.