Like many other additive manufacturing (AM) processes, fused deposition modeling (FDM) process is driven by a moving heat source, and temperature history plays an important role in determining the mechanical properties and geometry of the final parts. Thermal simulation of FDM is challenging due to geometric complexity of manufacturing process and inherent computational complexity which requires numerical solution at every time increment of the process. We describe a new approach to thermal simulation of the FDM process, formulated as an explicit finite difference method that is applied directly on as-manufactured model described by a typical manufacturing process plan. The thermal model accounts for most relevant thermal effects including heat convection and radiation to the environment, heat conduction with build platform and between adjacent roads (and adjacent layers). We show that the proposed simulation method achieves linear time complexity both theoretically and numerically. This implies that the simulation not only scales to handle three-dimensional (3D) printed components of arbitrary complexity but also can achieve real-time performance. The approach is fully implemented, validated against known analytic solutions, and is tested on realistic complex shapes.

Introduction

Motivation.

Additive manufacturing (AM) can fabricate parts with complicated shapes and internal structures without a significant increase in cost or turnaround time. It has been used to build prototypes and even functional parts in many industries. As one of the most widely used and rapidly growing AM technologies, fused deposition modeling (FDM) has the advantages of flexibility, low manufacture time, and cost-effectiveness. Like many other additive manufacturing processes (e.g., selective laser sintering (SLS)), FDM process is driven by a moving heat source, and temperature history plays an important role in determining the mechanical properties and geometry of the printed parts. The formation of bonds in FDM process is induced by the thermal energy of semi-molten material [1], and the repeated heating and cooling of the materials can aggregate nonuniform thermal gradient and cause stress accumulated that consequently result in part distortion [2]. Thermal analysis not only is important to understand different manufacturing failures but also provides basis for the design of manufacturing process plan.

To be manufactured using FDM, a designed solid model is converted into a process plan that specifies the printer head tool path (e.g., G-code) together with the printing process specifications, such as build direction, nozzle diameter, and infill percentage [3]. Perhaps, a most common example of this scenario is a process plan generated by process planning software such as Slice3r from stereolithography format describing the solid's boundary representation. In the course of manufacturing, as the printer head moves, a semi-molten filament is extruded through a heated nozzle. For each layer, the nozzle moves following a piecewise linear path horizontally. The material deposited along each line segment is commonly referred as a road. A schematic of the FDM extrusion process is shown in Fig. 1.

Thermal simulation of FDM is challenging at least for two reasons. The first one is the inherent geometric complexity of the manufacturing process that discretizes the design model and leads to numerous artifacts, such as “stair-stepping” on the surface, voids between roads, and underfill [4]. The second challenge is the high computational complexity, stemming from the fact that the transient thermal simulation requires numerical solution at every time increment of the process. Common simplified approaches to thermal simulation in additive manufacturing rely on voxelization of the designed geometry followed by transient finite element analysis. However, since the design geometry differs significantly from the manufactured shape [3], this approach is neither accurate nor efficient. For example, heat conductance along the road is usually greater than that between the adjacent roads, which suggests that accurate thermal simulation must account for detailed road geometry and all relevant thermal effects.

Approach and Contribution.

We describe a new approach to thermal simulation of the FDM process, formulated as an explicit finite difference method that is applied directly on as-manufactured model described by a typical manufacturing process plan. In our model, the manufacturing primitive-deposited road is discretized into road elements according to the tool path plan. The thermal model accounts for most relevant thermal effects including heat convection and radiation to the environment, heat conduction with build platform, between adjacent roads as well as between adjacent layers. Heat transfer between adjacent road elements is determined by their contacts. The simulation achieves linear time complexity by exploiting spatial and temporal locality, based on known empirical data, which allows updating at most a constant number of road elements at each time-step. This implies that our simulation not only scales to handle three-dimensional (3D) printed components of arbitrary complexity but also can achieve real-time performance. The output of the simulation is the complete temperature history of each road element over the course of manufacturing process.

Outline.

The remainder of this paper is organized as follows: Related work is briefly discussed in Sec. 2. The main contributions of the paper are contained in Sec. 3. The first principles model and dimensional reduction are introduced in Secs. 3.1 and 3.2. Sections 3.3 and 3.4 describe the preprocessing of the manufacturing process plan into an efficient data structure, including road discretization and contact detection. Numerical scheme is given in Sec. 3.5. A basic stability and convergence analysis are discussed in Sec. 3.6. Section 3.7 describes the concepts of active body and lazy update as means to utilize spatial and temporal locality to attain real time performance. Section 4 is devoted to experimental validation of the proposed method by comparing it with basic analytic solution in Sec. 4.1 and studying influence of active body in Sec. 4.2. Section 4.3 shows a realistic numerical example. Section 5 is reserved for some concluding remarks and discussion of future work.

Related Work

Broadly speaking, the previous research on thermal simulation of FDM can be divided into two categories that we briefly discuss below: simplifying the problem such that it can be solved analytically or approximating transient heat transfer on voxelized design models.

Several analytical models have been developed to predict the thermal history of a road [5]. Thomas and Rodriguez [6] presented a simplified two-dimensional thermal model that treated the roads as rectangle in shape. The analytic solution for the temperature averaged over width of the road could be obtained. The model neglected the effects of conduction to the build surface and any contact resistance between filaments, the former of which is the predominant heat transfer mechanism in the system [5]. Observe that the two-dimensional analysis showed temperature gradients that rapidly become negligible along the width and height of the filament led Bellehumeur et al. [7,8] to propose a lumped capacity model which assumed uniform temperature distribution of cross-sectional area of filament, semi-infinite filament length, and constant heat transfer and convection coefficients. The cooling process is thus simplified into a one-dimensional (1D) transient heat transfer model. A single deposition road is modeled as a simple block (1D sweep with ellipse cross section). The governing equation is simplified to a 1D transient heat transfer equation which can be solved analytically. This approach has the advantage of simplicity but does not extend to simulation of realistic process plans with complex tool paths. Costa et al. [9] proposed an analytical solution for the transient heat transfer during filament deposition, taking into account contacts between filaments and assuming a relatively simple deposition sequence. Both radial and axial heat conduction are ignored in their model. In a more recent paper [10], the same authors examined the contribution of various thermal phenomena during the manufacturing process to overall heat transfer. The weakness of all analytic approaches is that simplified closed-form solutions are limited to simple geometry (e.g., only one road or simple deposition strategy) and cannot be applied to realistic parts and manufacturing processes.

Approaches based on voxelization of the design geometry rely on voxel-based numerical methods (e.g., finite element method) to solve a transient heat transfer problem in each time-step. Zhang and Chou [11] developed a finite element analysis model using “element activations” to simulate coupled mechanical and thermal phenomena in FDM as well as resulting stresses and distortions. Ji and Zhou [12] proposed three-dimensional transient thermal finite element model taking into a account temperature-dependent thermal conduction and heat capacity. Yardimci and coworkers [13,14] proposed a 1D thermal model in which each road is modeled as a 1D array of blocks. The thermal interactions with the environment and between roads are considered by including sink terms in the enthalpy form energy equation. Finite volume method was utilized for spatial and temporal discretization. More recently, a number of researchers adopted voxelization of the design geometry as a basis for simulating moving heat source problems common in additive manufacturing [1517]. As we explained earlier, such approaches cannot account for significant difference between design and as-manufactured models and tend to be computationally prohibitive.

Formulation

The First Principles Model.

In FDM, the thermal-plastic filament is heated past its glass transition temperature and is then deposited by an extrusion head which follows the tool path and the part is built from bottom up, one layer at a time. After the semi-molten material has been deposited onto the build platform, it gradually cools down to envelope temperature (the temperature of the air around the part) and bonds to the neighboring material diffusively. The bonding between adjacent roads in the same layer and in adjacent layers is driven by the thermal energy absorbed by the building material during heating and extrusion. The cooling process is transient and physically complex in its nature [7]. Thus in most cases, it is impossible to obtain the analytical solution.

The standard approach for the formulation of transient thermal problem is based on the conservation of energy, which states that the energy variation in a volume over a period of time is equal to the difference between heat produced and the outflowing heat. The conservation equation is written as
$ΔU=∫τ∫VσdVdt−∫τ∫S=∂Vq·ndt$
(1)
where U is the specific internal energy, V is the volume, τ represents an arbitrary time interval, σ is heat production density rate, $q$ is heat flux density, $S=∂V$ stands for the boundary (surface) of an arbitrary volume V, and $n$ is the unit vector points to the outer normal. Because of the arbitrary time and space setting, by applying the divergence theorem, one obtains
$∂U∂t=σ−∇·q$
(2)
Moreover, $q$ and U can be expressed as functions of the temperature T. According to Fourier's law, heat flux is proportional to temperature gradient
$q=−λ∇T$
(3)
where λ is thermal conductivity ($W/(m·K)$), $∇T$ is the temperature gradient, and the pointwise constitutive equation that links specific internal energy U to temperature T is
$U=U0+ρcT$
(4)
where ρ is the density of the material (kg/m3) and c is its specific heat capacity ($J/(kg·K)$). Assuming heat capacity is not temperature dependent and substituting Eqs. (3) and (4) in Eq. (2), one obtains the partial differential equation (PDE) as
$ρc∂T∂t=∇·λ∇T+σ$
(5)
In FDM, it is reasonable to assume that the thermal conductivity λ is not spatially dependent within the road, and σ = 0 because there is no inner heat source. Then, Eq. (5) simplifies to
$ρc∂T∂t=λ∇2T$
(6)

where $∇2$ is the Laplace operator which is defined as $∇2=(∂2/∂x2)+(∂2/∂y2)+(∂2/∂z2)$ in 3D Cartesian coordinates. Equation (6) is solved as an initial boundary value problem meaning that the equation must be accompanied by (i) a domain, here the evolving printing part, (ii) initial conditions specifying the temperature at every point in the domain at the starting time for the problem, and (iii) a set of conditions specifying either the temperature or its derivatives (heat flux) along the boundary.

Simplification and Dimensional Reduction.

The above PDE (Eq. (6)) expresses a general statement of conservation of energy on arbitrary time and space setting which usually does not admit analytic solutions and is challenging to solve numerically. However, the FDM process is essentially a one-dimensional serial rasterization of three-dimensional space by a moving heat source. This observation allows dramatic simplification of the formulated model.

Since the cross-sectional area of road is quite small, the temperature difference on cross sections of individual roads can be neglected. As a result, the road “collapses” to its axis (see Fig. 2). And the thermal problem for a continuous deposition of material on single road is simplified to a one-dimensional problem
$∂∂tT(t,x)=K∂2∂x2T(t,x)+f(T(t,x))$
(7)

where $K=(λ/ρc)$. Each cross section of road is now represented by a point on the road's axis. Hence, the original 3D model is reduced to 1D. The diffusion term on RHS acts as the heat conduction along road axis and term f(T(t,x)) accounts for the heat transferred on the perimeter of the cross section, which includes convection, radiation to the ambient environment, heat conduction to build platform as well as heat conductance with neighboring roads (see Fig. 3).

Discretization of Tool Path.

Conceptually, the proposed approach to simulation consists of two stages: preprocessing and execution. In the preprocessing stage, the manufacturing process plan is used to create a contact graph—a spatial data structure representing adjacency between discretized road elements. In the execution stage, the contact graph is used to perform transient thermal simulation, updating thermal history for all elements that are affected in each time-step.

In preprocessing, the roads in the tool path are discretized into road segments, which we call elements. The road elements will be used to simulate the actual deposition process, with road segments being deposited sequentially at distinct time steps. Note that the tool head can also move without extruding any material.

All sequential collinear road segments are merged into one road segment. Two user-defined parameters are used to discretize the road segments: maximum road element length Δxmax and minimum road element length Δxmin. A long road segment with length $L>Δxmax$ is divided into $⌊L/Δxmax⌋+1$ road elements. Since very short road element may cause numerical problems, a road element with length $L<Δxmin$ is ignored. This happens rarely because typically Δxmin is very small relative to machine's accuracy (e.g., Δxmin = 0.5 mm) and such segments should not appear in the process plan. As a result, it is guaranteed that the length of each road element lies in a range from Δxmin to Δxmax.

Discretization of a single layer of part from Ref. [18] are shown in Fig. 4, where tool path is represented by solid line and element centers are denoted by small circles. Note that the boundary road segments tend to be shorter, as they control the accuracy of the part's geometry set by the manufacturing processing software.

Initially, the elements are stored in an array sequentially according to the process plan. Each element record also includes its index, position in the global Cartesian coordinate, length, deposition time, and temperature history. Subsequently, as we describe next, this array is transformed into a contact graph data structure that adds explicit adjacency information for each road element.

Contact Graph Construction.

The contact information between adjacent elements is represented by a weighted contact graph, $G=(V,E)$, comprising a set $V$ of vertices (road elements) together with a set $E$ of edges corresponding to contacts between adjacent road elements. Contact of road elements i and j is represented by an edge $ei,j∈E$. The weight of edge ei,j is the contact area Ai,j. It is reasonable to assume that the position of an element and its contact information with neighboring elements remain unchanged during simulation. Hence, the contact graph is a static data structure that can be built once for a given process plan. It should also be clear that the size of the contact graph is linear in the number of road elements ($O(|V|)$), since every road element can contact at most constant number of neighbors.

We also note that the road elements are already ordered by the process plan, and any given element can be in contact only with elements in adjacent roads within the same layer, or with elements above or below its layer. The contact queries can be further localized using the usual binning method, without additional preprocessing, requiring $O(|V|)$ time to construct a complete contact graph. Each element i in the constructed graph stores the list of its contacting elements {j} and the contact areas Ai,j computed as described below.

In order to accurately estimate the contact area, consider the microphotograph of the cross section of an FDM part from Ref. [8] reproduced in Fig. 5. It suggests that the contacting surfaces of the road can be reasonably approximated by the four planar regions: two vertical sides of the road where they contact other roads within the same layer, and two horizontal sides that can contact the roads in the layers immediately above and below the road. The precise manner in which these sides are connected are not important, as long as we can estimate the contact and noncontact areas of the road. For concreteness, we approximate the road's cross section by a rectangle with four corners cut (Fig. 6). The shape of the cut corner is assumed to be an isosceles triangle. More generally, the size of cross-sectional shape depends on specific process parameters (e.g., print head velocity and rate of deposition), though we will not deal with this dependence in the present paper.

The contact parameters w and h are easily determined from elementary geometry of the cross section. H and W are preset manufacture parameters obtained from the process plan: H is the layer height and W is the road width. Another manufacture parameter, called extrusion factor e (0 < e ≤ 1), determines how full the rectangle W × H is filled by the road. So, the material conservation is expressed as
$WH−12(W−w)(H−h)=eWH$
(8)
where h and w represent the interface neck length between adjacent roads within the same layer and between adjacent layers, respectively (see Fig. 6). As shown in Fig. 5, h = 2y and the isosceles triangle implies
$W−w=H−h$
(9)
Combining Eqs. (8) and (9), we can obtain
$w=W−2(1−e)WH$
(10)

$h=H−2(1−e)WH$
(11)

With these parameters, the actual contact areas are determined by considering the contacts made by the sweeps of these cross sections along the road's trajectories. We must distinguish contacts between the roads within the same layer (Figs. 7 and 8) and contacts between the roads in adjacent layers (Fig. 9). In the first case, the contact area is approximated by a product of neck length h and either contact length $L′$ when the roads are parallel or road width W when the road end contacts another road.

For contacting road elements in adjacent layers, the contact area is determined by the overlap of the top view (build direction) projections of the road elements, as illustrated in Fig. 9. Computing area of two overlapping rectangles is a classroom problem in computational geometry (e.g., see Ref. [19]).

Numerical Scheme.

The formulation simplifies even further when it is applied to the discretization of roads into elements as discussed earlier in Sec. 3.3. Consider the Biot number of the road element, a dimensionless parameter which is defined as the ratio of the conductive heat resistance within the object to the convective heat transfer resistance across the object's boundary with a uniform bath of different temperatures [20]. In general, problems involving small Biot numbers (much smaller than one) are thermally simple, due to uniform temperature fields inside the body. Considering the thermal properties of the common material used in FDM, short element length, and small road cross section, it is easy to check that Biot number of the road element is less than 0.1, which implies that the heat conduction inside the element is much faster than the heat convection away from its surface, hence temperature gradient is negligible in the element. So, it is reasonable to assume the temperature within the element is completely uniform in space, although this spatially uniform temperature value changes over time. This leads to a kind of lumped capacitance model. Each element (lump) is treated as a capacitive reservoir which transfers heat in the system. This simplification reduces the state space of the system to a finite dimension. And the original PDE of a continuous time–space model of the thermal system is changed into an ordinary differential equation with a finite number of parameters. Writing the energy conservation equation for a specific road element during a specific time interval Δt, we can get
$ΔU=ρcVΔT=−∫Δt∫S=∂Vq·ndt$
(12)
where V is the volume of the road element and $S=∂V$ represents the boundary (includes the surface and two cross sections) of the element. Equation (12) is the integral form of the energy conservation. In order to calculate temperature change $ΔT$ in this time interval, we need to integrate the heat flux $q$ on the boundary of road element. This can be accomplished by calculating different thermal effects on different portions of the boundary.
Figure 3 shows different thermal effects on the boundary of a road element. Hence, the integration on RHS of Eq. (12) could be written as
$∫Δt∫S=∂Vq·ndt=Qcond+Qconv+Qradi+Qplat+Qroads$
(13)

where Qcond represents the heat conduction along the road, Qconv and Qradi denote convection and radiation from free surface to ambient environment, respectively, Qplat serves as heat conductance to build platform, and Qroads covers heat conductance to contact road elements within the same layer and in adjacent layers.

We approximate the temperature along the road by a piecewise constant function, which assumes that temperature is spatially uniform within each element due to its small Biot number. The expressions for each thermal effect included in Eq. (13) are described as follows:

First, consider the heat conduction along the road. Second spatial derivative is approximated by central difference, which is expressed as
$Qicond=∫Δtλ(Ai,i−1Ti(t)−Ti−1(t)0.5(Li+Li−1)+Ai,i+1Ti(t)−Ti+1(t)0.5(Li+Li+1))·dt$
(14)

where Ti(t) is the temperature of road element i at time t, Ai,j is the cross-sectional area between road elements i and j, Li is the length of road element i, and Δt is the current time-step.

In addition to heat conduction along the road, there is energy dissipated from free surface to ambient environment through heat convection, which plays a major role during the cooling process. And heat convection on the free surface can be expressed as
$Qiconv=∫ΔthconvAifree(Ti(t)−T∞)·dt$
(15)

where hconv is the convective heat transfer coefficient ($W/(m2·K)$). In Ref. [21], Rodriguez et al. adopted the Churchill correlation for natural convection to estimate $hconv=67W/(m2·K)$ for an elliptical road deposited at $270$ °C toward an environment at $70$ °C. $Aifree$ is the free surface area of road element i and $T∞$ denotes the ambient temperature.

Thermal radiation has usually been neglected in the previous work. The energy carried by the emission of electromagnetic waves is given by Stefan–Boltzmann law
$Qiradi=∫ΔtεσAifree(Ti(t)4−T∞4)·dt$
(16)

where ε is a dimensionless parameter which is called emissivity coefficient. In our simulation, ε equals to 0.96 [10]; σ is the Stefan–Boltzmann constant ($σ=5.670×10−8W/(m2·K4)$). Kelvin scale should be used when calculating radiation.

Furthermore, as a deposited road in the first layer contacts the cooler build platform, heat transfer by conduction develops between them. A temperature difference arises at the interface due to the thermal contact resistance. Since the mass of build platform is much higher than that of the road, temperature of build platform can be considered as constant. The conduction heat transfer with the build platform can be considered in the form of convection [7], with the thermal contact conductance hc. In this case, the expression for heat conduction between road and the build platform becomes
$Qiplat=∫ΔthcAiplat(Ti(t)−Tplat)·dt$
(17)

where $Aiplat$ is the contact area of the road element i with build platform, Tplat is the temperature of build platform, and hc is the thermal contact conductance ($W/(m2·K)$) between road and build platform.

The contact between two adjacent roads with different temperatures can also be formulated as thermal contact conductance. If $Si$ is a set of road elements contacting element i, heat transfer to contacting road elements could be written as
$Qiroads=∑j∈Si∫ΔthrAi,jc(Ti(t)−Tj(t))·dt$
(18)

where $Ai,jc$ is the contact area between road elements i and j, and hr is the thermal contact conductance coefficient ($W/(m2·K)$) between two contacting roads.

Applying the forward Euler method in time to approximate energy transferred in Eqs. (14)(18), the heat conservation equation (Eq. (12)) of a road element is transformed into the following finite difference approximation:
$ρcVi(Tin+1−Tin)=−Δtn(λ(Ai,i−1Tin−Ti−1n0.5(Li+Li−1)+Ai,i+1Tin−Ti+1n0.5(Li+Li+1))+hconvAifree(Tin−T∞)+εσAifree((Tin)4−T∞4)+hcAiplat(Tin−Tplat)+∑j∈SihrAi,jc(Tin−Tjn)),$
(19)

where Vi represents the volume of road element i (see Fig. 10), $Tin$ denotes the temperature of road element i at time tn, and $Δtn$ is the nth time-step, $Δtn=tn+1−tn$. It can be seen that the above numerical scheme is essentially forward-time central-space (FTCS) method applied on a nonuniform grid.

The thermal simulation is executed in a stepwise manner by activating the road elements in the contact graph that are affected by heat exchange during the FDM process. Before the process begins, all the road elements are inactive. The thermal analysis initiates with the current temperature distribution as the initial condition. In each time-step, the printer head either moves without extruding any material or extrudes amount of material associated with some road element in the graph. In the latter case, this road element is also activated, as well as all elements that are either affected by the newly deposited element or remain active from the previous time-step. We will refer to the collection of all active elements as “active body.” In each time-step, the above numerical scheme (Eq. (19)) is used to update temperature of elements in the active body. The update assumes that material extrusion takes place instantly at the start of each time-step, which does not require solving a system of equations and can be performed simultaneously (and in parallel) for all elements. The initial (final) cross section of the first (last) element on each road is treated as free boundary.

Elements that are not in the active body are thermally inactive at that time-step. An efficient but less accurate temperature update approach called “lazy update” is applied to update temperature of inactive elements. Detailed definition of both “active body” and “lazy update” are described in Sec. 3.7.

Stability and Convergence.

The following analysis shows that the proposed numerical scheme is conditional stable. From Eq. (19), we can see that the most strict stability condition is enforced by using FTCS method to approximate the diffusion term. FTCS can yield unstable solutions that oscillate and grow if $Δt$ is too large. The stability condition [22] for uniform grid is given by
$Δt<(Δx)22K$
(20)

where $K=(λ/ρc)$. Since in our simulation the element length lies in the $[Δxmin,Δxmax]$ range, we show that the numerical scheme is stable for the more restrictive case of uniform grid with $Δxmin=0.5 mm$ and acrylonitrile butadiene styrene (ABS) P400 material that is widely used in FDM. The thermal properties of ABS P400 were reported by Bellehumeur et al. [7], which are listed in Table 1.

Substituting the thermal properties and Δxmin into Eq. (20), stability condition becomes $Δt<1.5424 s$. We now consider two separate cases. During the manufacturing process, when the nozzle travels without extruding any material, travel time $ttravel=(distance/vtravel)$, where $vtravel$ is the speed of printer head during traveling without extruding any material. If $ttravel>0.1 s$, it is discretized such that time-step $Δt<0.1 s$. And when there is an element extruded in this time-step, time-step $Δt=(Δx/vinfill)$ where $vinfill$ is the speed of printer's head while extruding amount of material. Substituting this into Eq. (20), we conclude that the stability condition requires $vinfill>(2K/Δx)≥0.3242mm/s$. Since the deposition speed of FDM is usually much greater than 1 mm/s, it is safe to claim that the proposed numerical scheme is stable during simulating the manufacturing process. The thermal simulation may also continue after the 3D printing process completes; in this case, we simply choose the constant time-step of Δt = 0.1 s, which obviously lies in the stability region. Thus, in both cases, the proposed numerical scheme is stable. Therefore, the Lax equivalence theorem [23] implies that the proposed numerical procedure converges.

Active Body and Lazy Update.

Assuming that a single road element is deposited at each time-step, the complete thermal simulation of an FDM process with n elements will require $O(n2)$ temperature updates. This approach is not likely to be practical for realistic simulation scenarios with hundreds of thousands or million time steps.

We overcome this difficulty by observing that the majority of the road elements are not affected by deposition of a new element, because they either already cooled down or are far away from the newly deposited element. For example, empirical data suggest that for ABS P400 material, it takes about $p=8 s$ to cool down (close) to the envelop temperature [7] and the thermal influence of the newly deposited element does not propagate through more than q contacts with neighboring road elements (due to thermal contact resistance on interface). Thus, an element that is deposited long time ago or far away from the current deposition area has reached thermal equilibrium with its surroundings. So, updating temperature of all the deposited elements in every time-step is neither efficient nor necessary. A natural question arises then: What is the smallest set of elements that should be updated at every time-step in order to maintain reasonable accuracy of simulation? We answer the question by introducing the concept of active body.

Active body ($At⊂V$) of road elements at time t is the union of two sets $Atempt∪Aspatt$. The temporal active body $Atempt$ includes all elements deposited in the last p seconds immediately preceding time t, i.e., during the time interval $[t−p,t]$. The parameter p is a predefined constant based on physical observation of the cooling profile of a road. The elements in the temporal active body $Atempt$ are in contact with their neighbors, which may or may not be in $Atempt$. As these elements get reheated, they conduct the heat to their neighbors and so on. Formally, we define the contact q-neighborhood $Nqj$ of element j in the contact graph $V$ to be the set of elements ${k | dist(k,j)≤q}$, where dist(k, j) is the length of the shortest path between the elements k and j in the contact graph. However, as the elements in $Atempt$ cool off, their effect on neighboring elements diminishes. Hence, we assume that the reheating effects are significant only in the contact q-neighborhoods of elements in the set $Acoret⊂Atempt$, which is defined to be the set of elements deposited in the last s time-steps prior to time t. In other words, we define the spatial active body at time t as
$Aspatt=∪j∈AcoretNqj$
(21)

Figure 11 shows an example of a printing path, as well as the core, spatial, and temporal active bodies at a time instant t (element t is the newly deposited element which is in the third layer).

Informally, since the core elements are deposited most recently, they have high temperature and thus interact actively with their neighbors. The elements of the active body that are not in the core are about to leave the active body and have relatively low temperatures, which means their thermal interaction with the neighboring elements outside of the active body is relatively insignificant. These interactions are still captured during the update, but we assume that the thermal energy does not propagate any further.

During the simulation, the active body is updated in every time-step. If no new elements are deposited, both $Atempt$ and $Aspatt$ shrink in size as elements leave the active body. If a new road element is deposited at time t, it is added to the temporal active body $Atempt$ and its contact neighborhood is added to the spatial active body $Aspatt$. It should be clear that each update of active body takes O(1) time. Furthermore, we claim that, for any given discretization and fixed parameter values (p, q, and s), the number of elements in the active body is bounded by a constant at any given time. It is fairly obvious that the size of temporal active body $|Atempt|$ is bounded by a constant determined by p. To see that it is also true for spatial active body $Aspatt$, observe that the maximum number of elements in contact with a given road element is bounded. This stems from the requirement on minimum size (and hence finite possible contact area) of an individual element. It follows that the maximum degree k of a node in the contact graph is also bounded. A loose upper bound on the spatial active body is then obtained immediately as $|Aspatt|≤skq$. However, the bound is fairly conservative. For example, in a typical simulation, with p = 8 s, q = 3, s = 150, the size of active body $|At|$ is about 1000 elements.

In every time-step of the simulation, the temperatures of all elements that are active at time t (i.e., elements in active body $At$) are updated according to numerical scheme of Eq. (19). Since the cost of such an update is bounded by a constant, the whole simulation takes O(n) time, where n is the total number of time steps (which is also approximately equal to the number of deposited road elements).

Elements that are inactive at time t (elements that are not in active body $At$) have relatively low temperatures, and the net heat transfer of an inactive element with other elements is low. We ignore the heat exchange of inactive element with other elements. Instead, their temperatures are updated in a lazy fashion using Newton's law of cooling [20], which has an analytic solution. In other words, temperatures of inactive elements are updated on demand, when these elements are queried (e.g., for visualization purposes or when an element becomes active again). This method of updating temperatures of inactive elements does not affect the running time of simulation.

The size of active body is determined by the parameters p, q, and s. Their values are selected based on physical tests and desired tradeoff between the accuracy and the efficiency of simulation. A larger active body will require more temperature updates and hence longer running time, but it may improve accuracy of the computed results.

Validation and Tests

Simple Accuracy Test.

To test the accuracy of proposed simulation, we applied our simulation to a simplified problem—deposition of a single road—which has analytic solution [7] and compared the computed cooling profile with that given by analytic solution. A single deposition road is modeled as a one-dimensional elliptic sweep. The head moves at a constant speed of v along the x-axis when extruding. The origin of the reference coordinate system is set at the extrusion nozzle (see Fig. 1). Considering only heat conduction and convection, the governing partial differential equation reduces to
$ρcA∂T∂t=Aλ∂2T∂x2−hP(T−T∞)$
(22)
where A and P are the area and perimeter of the elliptic cross section, respectively. The boundary conditions are defined as
$T={T0x=0 and t≥0T∞x=∞ and t≥0$
(23)
where T0 is the deposition temperature and $T∞$ is the environment temperature. The origin is moving at a velocity $v (x=vt)$, and the time dependence term $(∂T/∂t)$ is transformed to
$∂T∂t=∂T∂x∂x∂t=∂T∂xv$
(24)
As shown in Ref. [7], Eq. (22) is then reduced to an ordinary differential equation with analytic solution
$T=T∞+(T0−T∞)e−mx$
(25)
with $m=((1+4αβ)−1)/2α$ and x = vt, where $α=(λ/ρcv)$ and $β=(hP/ρcAv)$.

We applied the proposed numerical simulation to this problem to compute the cooling profile of the road. As shown in Fig. 12, the simulation result using time-step Δt = 0.1 s is in good agreement with the analytic solution of Eq. (25).

Influence of Active Body.

In this section, the influence of active body on computation results and efficiency are studied using a more realistic part commonly used for testing FDM processes. The tool path for octopus model is shown in Fig. 13. After the path discretization, there are 59,595 road elements in the whole part.

We applied thermal simulation to the octopus part with and without active body. Parameters of active body are chosen as p = 8 s, q = 3, and s = 150. The running times of both simulations are plotted in Fig. 14. It can be easily seen that, without active body, the running time grows much faster ($O(n2)$) with regard to the number of time steps as opposed to linear time (O(n)) performance when using active body.

Next, the influence of active body on computation result is studied. In this part, we want to show that proper parameters (i.e., p, q, and s) of active body could be chosen based on physical observations, and using proper parameters to define active body, the simulation could obtain nearly the same result as that without active body. However, active body can reduce the number of temperature updates per time-step dramatically.

The value of p influences the initial cooling during which temperature drops under glass transitional temperature after deposition. Roughly speaking, q and s together control the shape of active body, among them q controls the depth of active body, whereas s controls the size of the horizontal cross section of active body. The choice of proper parameters depends on specific process plan (e.g., material properties, deposition temperature, and envelop temperature). We can choose p by observing the cooling profile of a single road deposition through either experiments or simulations. The value of p should be large enough for the road to cool from deposition temperature to transitional temperature. The value of q could be set as (or more than) the maximum penetration depth which could also be obtained by experiments or simulations. The choice of s depends on desired tradeoff between accuracy and efficiency. Based on our experiments, Acore should cover at least one-fourth to one-third of the hottest elements in Atemp to capture the reheating phenomena well.

The proposed simulation is applied to simulate the manufacturing of first ten layers of the octopus part, which has 24,340 road elements and 24,965 time steps in the discretized tool path. Based on our numerical experiments, active body parameters are chosen as $p=8 s$, q = 3, and s = 150. A road element in the third layer is picked to study its temperature evolution. The reference solution is given by simulation without active body. The comparison of the computed temperature evolutions is shown in Fig. 15.

It can be seen from Fig. 15 that the blue dashed line (simulation with active body) agrees well with the red line (simulation without active body). The mean absolute error is $0.3559 °C$ and the mean absolute percentage error is $0.48%$. Also, the average number of temperature updates per time-step with active body is 1149, whereas the average number goes up to 12,162 when active body is not used.

To confirm that the choice of parameters does not depend on geometry, we applied the same active body parameters (i.e., p = 8 s, q = 3, and s = 150) to simulate two other parts with the same material and manufacture parameters as the octopus part. Figure 16 shows the result of simulating the manufacturing of first ten layers of the mobius arm part (shown in Fig. 17). An element in the second layer is chosen to study its temperature evolution. It can be seen in Fig. 16 that the computed temperature evolution agrees well with the reference solution ($MAE=0.3823$ °C and $MAPE=0.49%$). Figure 18 shows the result of simulating the manufacturing of first five layers of the aerospace rocker part (shown in Fig. 19). An element in the second layer is picked to compare its temperature evolutions. From Fig. 18, we could see that active body also performed well with $MAE=0.4676$ °C and $MAPE=0.63%$.

Realistic Numerical Simulation.

The 3D model for a realistic part is shown in Fig. 17 and the tool path of the mobius arm part is shown in Fig. 20.

The tool paths for printing support material are shown in blue, whereas tool paths for printing the part's shape are shown in mixed colors. This part was discretized into 69,635 road elements. The temperature fields predicted by the thermal simulation at different times are shown in Fig. 21, where road elements are displayed as line segments with different colors according to their respective temperatures.

Complete simulation of 70,410 time steps in the printing process took 24.052 s to execute. The reference computer system configuration included 2.5 GHz quad-core Intel Core processor with 16 gigabytes of RAM, running Java 1.8 atop Apple macOS sierra 10.12.3. In comparison, it takes about 821.913 s to 3D print the part using an FDM printer whose nozzle moves at 6 cm/s.

Conclusion and Future Directions

Summary and Significance.

The main contribution of the paper is a new approach to thermal simulation of FDM process. The simulation is performed directly on as-manufactured model, derived from the manufacturing process plan, and includes detailed geometry of the printed part as well as the deposition order in which this geometry is created. The proposed simulation has three novel ingredients. First, our thermal model is formulated directly in terms of manufacturing process geometry (deposited materials) and covers all the important thermal effects developed during the manufacturing process, including heat convection and radiation to the environment, heat conduction with build platform and between adjacent roads (and adjacent layers). Second, the model is discretized and solved using explicit finite difference method which does not require solving system of linear equations at every time-step. Every manufacturing primitive—a deposited road in the case of FDM—is discretized into road elements which are used in our simulation. Finally, the discretized model is organized into efficient spatial-temporal data structure that allows efficient localization of computations in space and time, based on known physical evidence. The proposed method has been fully implemented and tested. As we showed both theoretically and experimentally, the simulation is numerically stable, converges, and achieves linear time complexity. This also implies that the simulation not only scales to handle 3D printed components of arbitrary complexity but also can achieve real-time performance.

Extensions and Promising Directions.

The proposed thermal simulation is not limited to FDM. Many of the results are applicable to other additive manufacturing methods which are driven by moving heat source. For example, SLS is driven by a moving laser beam. Unlike FDM, SLS does not appear to have the concept of road explicitly. However, the whole part can be seen as built by sweeping the melt pool which is controlled by laser power, scan speed, etc. Hence, the road could be defined as the sintered material by one scan of a laser beam. Another possible extension is to use the insights gained from simulation to design tool path plan to achieve desired temperature field(s). One challenging aspect of this research is that we are not aware of any other reference implementations or experimental measurements that could be used to validate the correctness of computed results. Finally, the proposed simulation has real time performance, which could be used to build a feedback control of manufacturing parameters.

Finally, the ability to efficiently simulate thermal history of 3D printing process opens up many possibilities for investigating mechanical properties of the manufactured part that are influenced by the thermal history. These include bond strength between adjacent roads [24,25] and warping [26].

Acknowledgment

The authors are also grateful to Paul Whitherell (NIST) and John Michopoulos (NRL) for many useful discussions. The responsibility for errors and omissions lies solely with the authors.

Funding Data

• National Institute of Standards and Technology (Grant No. 70 NANB14H232).

• National Science Foundation (Grant Nos. CMMI-1344205 and CMMI-1361862).

References

References
1.
Sun
,
Q.
,
Rizvi
,
G.
,
Bellehumeur
,
C.
, and
Gu
,
P.
,
2003
, “
Experimental Study of the Cooling Characteristics of Polymer Filaments in FDM and Impact on the Mesostructures and Properties of Prototypes
,”
14th Solid Freeform Fabrication Symposium
(SSF), Austin, TX, Aug. 4–6, pp.
170
178
.
2.
Zhang
,
Y.
, and
Chou
,
Y. K.
,
2006
, “
A Parametric Study of Part Distortions in FDM Using 3D FEA
,”
17th Solid Freeform Fabrication Symposium
(
SSF
), Austin, TX, Aug. 14–16, pp.
410
420
.
3.
Liu
,
X.
, and
Shapiro
,
V.
,
2016
, “
Homogenization of Material Properties in Additively Manufactured Structures
,”
Comput.-Aided Des.
,
78
, pp.
71
82
.
4.
Vasudevarao
,
B.
,
Natarajan
,
D. P.
,
Henderson
,
M.
, and
Razdan
,
A.
,
2000
, “
Sensitivity of RP Surface Finish to Process Parameter Variation
,”
Solid Freeform Fabrication Proceedings
(
SSF
), Austin, TX, Aug. 7–9, pp.
251
258
.
5.
Turner
,
B. N.
,
Strong
,
R.
, and
Gold
,
S. A.
,
2014
, “
A Review of Melt Extrusion Additive Manufacturing Processes—I: Process Design and Modeling
,”
Rapid Prototyping J.
,
20
(
3
), pp.
192
204
.
6.
Thomas
,
J.
, and
Rodriguez
,
J.
,
2000
,
Modeling the Fracture Strength Between Fused Deposition Extruded Roads
,
University of Texas at Austin
,
Austin, TX
.
7.
Bellehumeur
,
C.
,
Li
,
L.
,
Sun
,
Q.
, and
Gu
,
P.
,
2004
, “
Modeling of Bond Formation Between Polymer Filaments in the Fused Deposition Modeling Process
,”
J. Manuf. Process
,
6
(
2
), pp.
170
178
.
8.
Sun
,
Q.
,
Rizvi
,
G.
,
Bellehumeur
,
C.
, and
Gu
,
P.
,
2008
, “
Effect of Processing Conditions on the Bonding Quality of FDM Polymer Filaments
,”
Rapid Prototyping J.
,
14
(
2
), pp.
72
80
.
9.
Costa
,
S. F.
,
Duarte
,
F. M.
, and
Covas
,
J. A.
,
2008
, “
Towards Modelling of Free Form Extrusion: Analytical Solution of Transient Heat Transfer
,”
Int. J. Mater. Form.
,
1
(
S1
), pp.
703
706
.
10.
Costa
,
S. F.
,
Duarte
,
F. M.
, and
Covas
,
J. A.
,
2015
, “
Thermal Conditions Affecting Heat Transfer in FDM/FFE: A Contribution Towards the Numerical Modelling of the Process
,”
Virtual Phys. Prototyping
,
10
(
1
), pp.
35
46
.
11.
Zhang
,
Y.
, and
Chou
,
Y.
,
2006
, “
Three-Dimensional Finite Element Analysis Simulations of the Fused Deposition Modelling Process
,”
Proc. Inst. Mech. Eng., Part B
,
220
(
10
), pp.
1663
1671
.
12.
Ji
,
L. B.
, and
Zhou
,
T. R.
,
2010
, “
Finite Element Simulation of Temperature Field in Fused Deposition Modeling
,”
,
97–101
, pp.
2585
2588
.
13.
Yardimci, A.
,
M.
, and
Geri
,
S.
,
1996
, “
Conceptual Framework for the Thermal Process Modelling of Fused Deposition
,”
Rapid Prototyping J.
,
2
(
2
), pp.
26
31
.
14.
Yardimci
,
A. M.
,
Hattori
,
T.
,
Guceri
,
S.
, and
Danforth
,
S.
,
1997
, “
Thermal Analysis of Fused Deposition
,”
Solid Freeform Fabrication Conference
(
SFF
), Austin, TX, Aug. 11–13, pp.
689
698
.
15.
Zeng
,
K.
,
Pal
,
D.
, and
Stucker
,
B.
,
2012
, “
A Review of Thermal Analysis Methods in Laser Sintering and Selective Laser Melting
,”
Solid Freeform Fabrication Symposium
(
SFF
), Austin, TX, Aug. 6–8, pp.
796
814
.
16.
Krishnakumar
,
A.
,
Suresh
,
K.
, and
Chandrasekar
,
A.
,
2015
, “
Towards Assembly-Free Methods for Additive Manufacturing Simulation
,”
ASME
Paper No. DETC2015-46356.
17.
Romano
,
J.
,
,
L.
, and
,
M.
,
2015
, “
Thermal Modeling of Laser Based Additive Manufacturing Processes Within Common Materials
,”
Procedia Manuf.
,
1
, pp.
238
250
.
18.
McCarthy
,
B.
,
2015
, “
Characterization of Geometric Deviations in FDM
19.
Sutherland
,
I. E.
, and
Hodgman
,
G. W.
,
1974
, “
Reentrant Polygon Clipping
,”
Commun. ACM
,
17
(
1
), pp.
32
42
.
20.
Incropera
,
F. P.
, and
DeWitt
,
D. P.
,
1996
,
Fundamentals of Heat and Mass Transfer
,
Wiley
, New York.
21.
Rodriguez
,
J. F.
,
Thomas
,
J. P.
, and
Renaud
,
J. E.
,
2000
, “
Characterization of the Mesostructure of Fused-Deposition Acrylonitrile-Butadiene-Styrene Materials
,”
Rapid Prototyping J.
,
6
(
3
), pp.
175
186
.
22.
Recktenwald
,
G. W.
,
2004
, “
Finite-Difference Approximations to the Heat Equation
,”
Mech. Eng.
,
10
, pp.
1
27
.
23.
LeVeque
,
R. J.
,
1992
,
Numerical Methods for Conservation Laws
,
, Cham, Switzerland.
24.
Li
,
L.
,
Sun
,
Q.
,
Bellehumeur
,
C.
, and
Gu
,
P.
,
2001
, “
Composite Modeling and Analysis of FDM Prototypes for Design and Fabrication of Functionally Graded Parts
,”
Solid Freeform Fabrication Symposium
(
SFF
), Austin, TX, Aug. 6–8, pp.
187
194
.
25.
Rosenzweig
,
N.
, and
Narkis
,
M.
,
1981
, “
Sintering Rheology of Amorphous Polymers
,”
Polym. Eng. Sci.
,
21
(
17
), pp.
1167
1170
.
26.
Wang
,
T.-M.
,
Xi
,
J.-T.
, and
Jin
,
Y.
,
2007
, “
A Model Research for Prototype Warp Deformation in the FDM Process
,”