Linear-Time Thermal Simulation of As-Manufactured Fused Deposition Modeling Components

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.

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.

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.

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 [15–17]. As we explained earlier, such approaches cannot account for significant difference between design and as-manufactured models and tend to be computationally prohibitive.

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.

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.

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.

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).

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 Δx_max and minimum road element length Δx_min. 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 Δx_min is very small relative to machine's accuracy (e.g., Δx_min = 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 Δx_min to Δx_max.

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.

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 e_i,j is the contact area A_i,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 A_i,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.

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]).

where Q^cond represents the heat conduction along the road, Q^conv and Q^radi denote convection and radiation from free surface to ambient environment, respectively, Q^plat serves as heat conductance to build platform, and Q^roads 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:

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

where h_conv 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.

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.

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

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

where V_i represents the volume of road element i (see Fig. 10), $Tin$ denotes the temperature of road element i at time tⁿ, 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.

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 Δx_min 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.

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.

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.

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).

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, A_core should cover at least one-fourth to one-third of the hottest elements in A_temp 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%$⁠.

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.

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.

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].

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.

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

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