Mesoscale Modeling of Extrusion and Solidification During Material Extrusion Additive Manufacturing
Jeffrey B. Allen ^{1,2,*}, Guillermo A. Riveros ^{1}, Ivan P. Beckman ^{1}, Elton L. Freeman ^{1}

Information Technology Laboratory, U.S. Army Engineer Research and Development Center, Vicksburg, MS 39180, USA

Department of Mechanical and Aerospace Engineering, Utah State University, Logan, UT 84322, USA
* Correspondence: Jeffrey B. Allen
Academic Editor: Salman Pervaiz
Special Issue: Research and Development of Subtractive and Additive Manufacturing Technologies
Received: August 22, 2023  Accepted: October 19, 2023  Published: October 25, 2023
Recent Progress in Materials 2023, Volume 5, Issue 4, doi:10.21926/rpm.2304033
Recommended citation: Allen JB, Riveros GA, Beckman IP, Freeman EL. Mesoscale Modeling of Extrusion and Solidification During Material Extrusion Additive Manufacturing. Recent Progress in Materials 2023; 5(4): 033; doi:10.21926/rpm.2304033.
© 2023 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.
Abstract
In this work, we apply a multiphysics approach to fused deposition modeling to simulate extrusion and solidification. Restricting the work to a single line scan, we focus on the application of polylactic acid. In addition to heat, momentum and mass transfer, the solid/liquid/vapor interface is simulated using a fronttracking, levelset method. The results focus on the evolving temperature, viscosity, and volume fraction and are cast within a set of parametric studies, to include the printing and extrusion speed, as well as the extrusion temperature. Among other findings, it was observed that fused deposition modeling can be effectively modeled using a fronttracking method (i.e. the level set method) in concert with a temperature dependent porosity function. The use of the levelset method for discriminating the phase change interface in this context is relatively new and offers considerable advantages over existing methods.
Keywords
Fused deposition modeling; level set method; multiscale modeling
1. Introduction
Of the various forms of additive manufacturing (AM), extrusionbased approaches, such as fused deposition modeling (FDM), or fused filament fabrication (FFF), have become very popular [1]. The attraction to FDM is primarily due to its high reliability, low maintenance, large variety of applicable filament materials, and low initial investment cost [1]. FDM is known to work with a variety of different polymers such as acrylonitrile butadiene Styrene (ABS) and polylactic acid (PLA), as well as polymer matrix composites, polymer ceramic composites, nanocomposites, and fiberreinforced composites [2].
The primary liabilities associated with the FDM process typically relate to poor product quality and inconsistent reproduction efforts, which arise from various defects, including void formations, surface roughness, and inadequate bonding between layers. To some degree, these defects can be managed by proper control of the processing parameters (i.e. printing speed, extrusion speed, extrusion temperature, etc.), but in many instances an additional pre and/or post treatment is required (i.e. chemical, heat, laser, or ultrasound) [3].
While there have been a large number of experimental studies related to various aspects of FDM, including the potential for defect mitigation [4,5,6], complementary theoretical/numerical studies are comparatively few in number and offer considerable avenues for advancement in the field. Typically, the numerical studies incorporate a specific thermodynamic or thermomechanical type of analysis and are developed for the meso/macro scales using the finite element method (FEM). Studies of this type include for example, the inner flow dynamics of the nozzle [7,8], filament cooling subsequent to deposition [9], parts distortion and roughness [10,11,12], filament bonding [13,14,15], and postprocess mechanical property analysis [16,17,18,19,20,21]. Attempts to simulate the entire deposition process (including nozzle flow, deposition, cooling and solidification) are rare. Possibly the earliest effort was conducted by Bellini [22] wherein a twodimensional geometry was used to simulate heat transfer and fluid flow of ceramic filaments. Later more complicated models include the work of Du et al. [23] which adopted the Volume of Fluid (VOF) method for capturing the evolving fluid interface and modeled the polymer viscosity using shear rate and temperature dependence. The use of the fronttracking, or levelset method for interface simulation, conducted for example by Xia et al. [24] facilitated a more robust multiphysics approach, wherein the conservation equations (including mass and momentum) could be coupled and solved on a fixed grid. In particular these fronttracking methods utilize an implicit (rather than explicit) formulation of the multiphase interface that is represented by a separate time dependent partial differential equation for the solution of the level set field variable. Other recent simulations of this type may be found in [24,25,26,27]. These latter simulations typically assume Newtonian nonviscoelastic fluid flow conditions, and incorporate melt viscosities that vary in complexity (i.e. ranging from a constant value to functional dependence on temperature and shear rate). Additional effects, such as solidification (evolving simultaneously with the cooling process), are often neglected.
In this work we further simulate the FDM extrusion and deposition process at the mesoscale, while including solidification effects and nonNewtonian viscoelastic flow conditions. Restricting the work to a single line scan with a continuous extrusion flux of PLA polymer, the numerical simulations presented herein will be conducted with the assistance of the FEM based, multiphysics software, COMSOL [28], and incorporate a 2D, heat and mass transfer analysis through a multiphase approach. The free surface, liquid/vapor interface will be simulated using a fronttracking, levelset method, and the melt viscosity will assume both temperature and shear rate dependency. The inclusion of microscale phenomena, for example those induced by flow induced crystallization (crystallization kinetics), while important for the determination of various thermosphysical sample properties, was not performed in this work for purposes of scope limitations. Postprocessing outcomes will focus on the evolving temperature, viscosity, and volume fraction and combine parametric studies associated with the various process control parameters, to include printing and extrusion speed, as well as extrusion temperature. Initial benchmark applications focusing on an evolving melting front will serve to validate the heat transfer/phase change numerical framework.
2. Model Development & Numerical Considerations
2.1 Model Assumptions
As with all numerical/theoretical assessments attempting to replicate a physical process, there will inevitably exist some limiting assumptions. In this work, we assume the following:
 The geometry is assumed two dimensional with coordinates (x, y).
 The surface tension coefficient (γ) is assumed constant, and therefore the Marangoni effect is neglected.
 Heat losses due to vaporization and radiation effects are omitted because of the relatively low fusion temperatures (i.e., the vaporization temperature of PLA, ~644 K, is not reached), and because of the relatively large distance between the nozzle and substrate, respectively. This is in conformance with the assumptions from other works [29,30].
 Model solidification assumes that the fluid flow within the transition zone (between solid and fluid) is similar to the flow in a porous media (i.e. we define a temperature dependent porosity function f_{1}(T)).
 A modified heat capacity is used within the energy equation ($C_{p}^{eq}$) which accounts for the latent heat of fusion.
 The liquid phase is assumed incompressible (density is assumed constant) and laminar.
2.2 Model Equations
Subject to the above assumptions, the conservation equations for the energy, mass and momentum may be expressed as:
Energy:
$$\begin{array}{}\text{(1)}& \rho {C}_{p}(\frac{\mathrm{\partial}T}{\mathrm{\partial}t}+\mathrm{\nabla}\cdot (\mathit{u}T))=k{\mathrm{\nabla}}^{2}Th(T{T}_{0})\delta (\varphi )\end{array}$$
Mass:
$$\begin{array}{}\text{(2)}& \mathrm{\nabla}\cdot \mathit{u}=0\end{array}$$
Momentum:
$$\begin{array}{}\text{(3)}& \rho (\frac{\mathrm{\partial}\mathit{u}}{\mathrm{\partial}t}+\mathit{u}\cdot \mathrm{\nabla}\mathit{u})=\mathrm{\nabla}\cdot [pI+\mu (\mathrm{\nabla}\mathit{u}+(\mathrm{\nabla}\mathit{u}{)}^{T})]+\rho \mathit{g}+{\mathit{F}}_{s}++\gamma \mathit{n}\kappa \delta (\varphi )\end{array}$$
Where T is the temperature, u is the velocity, t is the time, ρ is the density, and k is the thermal conductivity. Additionally, p is the pressure, I is the identity matrix, K is the curvature, g is the acceleration due to gravity, γ is the surface tension coefficient, φ is the level set variable, n is the interface normal unit vector (where δ(φ) is used to track the interface and is defined in Equation 5), F_{s} is a momentum force due to solidification, and (∇u) ^{T} is the transposed vector of velocity gradients.
The conservation equations are complemented with a fronttracking method (in this case the levelset method) to track the solidliquidgas interface. As opposed to other fronttracking methods, the levelset method uses a fixed mesh and dependent variable φ defined over the entire computational domain. In this work, φ is prescribed a value of 0 within the liquid (and solid phases) and a value of 1 within the vapor (air) phase. The dynamics of the interface, and the associated advection of φ, is determined in part from solutions of velocity. Further, the phase transition (represented by the prescribed thickness of the interface) is governed by a smooth step function, allowing for the mitigation of numerical divergence.
The evolution of φ may be expressed as:
$$\begin{array}{}\text{(4)}& \frac{\mathrm{\partial}\varphi}{\mathrm{\partial}t}+\mathit{u}\cdot \mathrm{\nabla}\varphi ={\gamma}_{\mathrm{l}\mathrm{s}}\mathrm{\nabla}\cdot ({\epsilon}_{\mathrm{l}\mathrm{s}}\mathrm{\nabla}\varphi \varphi (1\varphi )\frac{\mathrm{\nabla}\varphi}{\mathrm{\nabla}\varphi })\end{array}$$
Where γ_{ls}, and ε_{ls} are constants representing the distance of transition, and speed of reinitialization, respectively. Additionally, the interface delta function (used to track the solid/liquid/gas interface), may be expressed as:
\[ \delta(\phi)=6\phi(1\phi)\nabla\phi \tag{5} \]
Solidification:
The solidification model used in this work is based on the model developed by Voller and Prakash [31], and is used to compute the force F_{s} shown in Equation 3. First, a temperature dependent liquid fraction function (f_{l}(T)) is defined as follows:
$$\begin{array}{}\text{(6)}& {f}_{1}(T)=\{\begin{array}{l}{\textstyle \phantom{\rule{1em}{0ex}}}{\textstyle \phantom{\rule{1em}{0ex}}}{\textstyle \phantom{\rule{1em}{0ex}}}{\textstyle \phantom{\rule{thickmathspace}{0ex}}}0,T\ge ({T}_{\mathfrak{m}}+\epsilon )\\ \frac{({T}_{\mathfrak{m}}+\epsilon T)}{2\epsilon},({T}_{\mathfrak{m}}+\epsilon )>T\ge ({T}_{\mathfrak{m}}\epsilon )\\ {\textstyle \phantom{\rule{1em}{0ex}}}{\textstyle \phantom{\rule{1em}{0ex}}}{\textstyle \phantom{\rule{1em}{0ex}}}{\textstyle \phantom{\rule{thickmathspace}{0ex}}}1,T<({T}_{\mathfrak{m}}\epsilon )\end{array}\end{array}$$
Where T_{m} is the average melting temperature, and 2ε represents the temperature of the transition zone; namely the temperature interval between a purely liquid and solid phase. Note that for pure materials ε = 0.
An intermediate function, A, is then defined as:
\[ A=\frac{c_{1}(1\psi)^{3}}{(\psi^{3}+c_{2})} \tag{7} \]
Where ψ = 1  f_{1}(T), and c_{1} and c_{2} are constants (with c_{1} and c_{2} prescribed with comparatively large and small values, respectively). F_{s} may then be written as:
$$\begin{array}{}\text{(8)}& {F}_{s}=A\mathit{u}\varphi \end{array}$$
Where φ ensures that F_{s} acts only on the liquid phase and not the gas phase.
Thus we see that when T is greater than T_{m} + ε (i.e. greater than the liquidus temperature), ψ = 1 and A = 0, representing a fully fluid phase. However, when T is less than T_{m}  ε, ψ = 0, and A = c_{1}/c_{2}. If c_{1} is sufficiently large (with c_{2} small), F_{s} dominates the convective and diffusive terms within the momentum equation (see Equation 3) and thus forces the velocity to become zero. For the intermediate case within the transition zone (i.e., (T_{m} + ε) > T ≥ (T_{m}  ε)), f_{1}(T) takes on a value between 0 and 1 and the flow behaves similar to that of a porous media.
Finally, to account for the latent heat of fusion of the polymer, a modified heat capacity equation is used:
\[ C_{p,\mathrm{PLA}}^{eq}=C_{p,s}+L_{f}\frac{e^{\left(\frac{(TT_{\mathrm{m}})^{2}}{(T_{1}T_{s})^{2}}\right)}}{\sqrt{\pi(T_{1}T_{s})^{2}}} \tag{9} \]
Where T_{1} = T_{m} + ε, T_{s} = T_{m}  ε, L_{f} is the latent heat of fusion, and C_{p, s} is the solid phase, isobaric heat capacity.
The total heat capacity, accounting for both the air (c_{p, a}) and the polymer can then be written as:
\[ C_{p}=C_{p,\mathrm{a}}+\bigl(C_{p,\mathrm{PLA}}^{eq}C_{p,\mathrm{a}}\bigr)\phi \tag{10} \]
Viscosity:
The polymer dynamic viscosity follows the nonNewtonian, Carreau model [32,33], allowing for temperature and shear rate dependence, and is expressed as:
\[ \mu_{\mathrm{PLA}}=\mu_{\infty}+(\mu_{0}\mu_{\infty})[1+(\lambda\dot{\gamma})^{2}]^{\frac{n1}{2}} \tag{11} \]
Where μ_{0} and μ_{∞} are the zero shear rate viscosity and infinite shear rate viscosity, respectively (note for this work μ_{∞} = 0), n is the power index, and λ is the relaxation time. In Equation 11, we observe that for low shear rates (i.e., $\dot{\gamma}$ ≪ 1/λ), the fluid behaves as a Newtonian fluid with viscosity μ_{0}. For intermediate shear rates (i.e., $\dot{\gamma}$ > 1/λ), the fluid behaves as a Powerlaw fluid (i.e. μ = K (∂u/∂y) ^{n1}), and at high shear rates (i.e., $\dot{\gamma}$ ≫ 1/λ), the fluid again behaves as a Newtonian flued with viscosity μ_{∞}.
The shear rate, $\dot{\gamma}$, is expressed in terms of the rate of deformation tensor (S):
$$\begin{array}{}\text{(12)}& \dot{\gamma}=\sqrt{2\mathit{S}:\mathit{S}}\end{array}$$
Where
$$\begin{array}{}\text{(13)}& \mathit{S}=\frac{1}{2}[\mathrm{\nabla}\mathit{u}+(\mathrm{\nabla}\mathit{u}{)}^{T}]\end{array}$$
Finally, the total dynamic viscosity (μ), density (ρ), and thermal conductivity (K) may be expressed as:
\[ \mu=\mu_{\mathrm{a}}+(\mu_{\mathrm{PLA}}\mu_{\mathrm{a}})\phi \tag{14} \]
\[ \rho=\rho_{\mathrm{a}}+(\rho_{\mathrm{PLA}}\rho_{\mathrm{a}})\phi \tag{15} \]
\[ \kappa=\kappa_{\mathrm{a}}+(\kappa_{\mathrm{PLA}}\kappa_{\mathrm{a}})\phi \tag{16} \]
2.3 Geometry Mesh and Initial/Boundary Conditions
The geometry and initial/boundary conditions are shown in Figure 1(a). With respect to the geometry, the magnitudes for the standoff distance (h_{s}) and nozzle diameter (D_{N}) are specified in Table 1. As indicated, the nozzle translates along the xdirection with a constant printing speed (V) made feasible by the use of COMSOL’s moving mesh facility [28]. The continuity boundary conditions along this moving interface assumes an outwarddirected, convective heat transfer, as well as a zero pressure condition. At the nozzle boundaries, the extrusion temperature (T_{e}) is maintained over the total nozzle transit time. The remaining boundaries are prescribed as shown in Figure 1(a). A portion of the prescribed mesh is shown in Figure 1(b). A mesh density of 12,834 triangular mesh elements with increased levels of refinement along the moving nozzle/substrate interface was deemed appropriate subsequent to a series of preliminary numerical stability and convergence tests involving different levels of mesh refinement.
Table 1 Thermophysical properties of PLA & Operational Parameters.
Figure 1 Computational domain; axisymmetric coordinate reference, boundary and initial conditions (a) and mesh (b).
2.4 Thermophysical Properties
The thermophysical material and rheological properties corresponding to PLA are shown in Table 1. As shown, several of the properties (λ, μ_{0}, n) relating to the viscosity are functions of temperature. Note that in this work, the thermal properties of PLA are constant for both liquid and solid phases. The operational parameters, corresponding to extrusion speed (U), extrusion temperature (T_{e}), printing speed (V), nozzle diameter (D_{N}), and standoff distance (h_{s}) are also indicated in Table 1 and include in certain cases (i.e. U, T_{e}, V) a range of values for purposes of parameterization studies.
3. Discussion and Results
3.1 Preliminary Model Validation
3.1.1 1D Heat Transfer with Phase Change
Preliminarily, and as a form of validation for the aforementioned equations applicable to heat transfer with phase change, a one dimensional benchmark application involving a melting front is considered (see Figure 2(a), inset). This particular problem was selected due to the availability of an analytical solution [39]. Specifically, we consider a column of ice/water of length x = 1 m, and initial temperature T_{0} = 283 K subject to a fixed temperature of 253 K at its leftmost boundary (x = 0 m). The problem was discretized using 50 linear finite elements over a total simulation time of 72.0E3 s.
Figure 2 Melting front example; Temperature profiles along 1D domain showing phase transition (a) and comparisons of current numerical and analytical solutions [39] at t = 72.0E3 s. (b).
As shown in Figure 2(a), the interface separating the frozen and liquid water is seen to transition from left to right. Figure 2(b) shows the good agreement between the current numerical model with the analytical results [39].
3.2 FDM Simulations
The following results pertain to the aforementioned 2D FDM simulations, for which, unless otherwise specified, the input thermophysical properties and operation properties are given in Table 1 (i.e., V = 0.10 m/s; U = 0.20 m/s, and T_{e} = 200°C). Figure 3 shows the time evolution of the PLA volume fraction over a time period of 0.12 seconds, which corresponds to a nozzle transit distance of approximately 80% of the total length of the underlying substrate. Here the surrounding vapor (air) phase is shown in red (i.e., φ = 1) while the extruded PLA liquid/solid phase is shown in blue (i.e., φ = 1). The interface region is clearly distinguished by the intermediate liquid/vapor color transition. As indicated (for the input properties specified), the transiting nozzle extrudes the liquid PLA in a nearly continuous bead with a height of approximately 0.75 h_{s}.
Figure 3 Time evolution of the PLA volume fraction (V = 0.10 m/s; U = 0.20 m/s; T_{e} = 476 K).
Figure 4 shows the effect of varying the extrusion speed from 0.10 m/s to 0.25 m/s, with all other input conditions as prescribed in Table 1. As indicated, the PLA volume fraction becomes a series of irregularly spaced, disjointed circular discs for U = 0.10 m/s. At U = 0.15 m/s, a continuous bead of height ~0.5 h_{s} eventually develops subsequent to the formation of several discontinuous circular bead formations. Not until the extrusion speed approaches 0.20 m/s (see Figure 4(c)) do we observe a nearly continuous uninterrupted PLA bead of height ~0.75 h_{s}. Finally, for U = 0.25 m/s, the PLA bead height becomes nearly equivalent to h_{s}.
Figure 4 Effect of variable extrusion speed on PLA volume fraction (V = 0.10 m/s; T_{e} = 473 K; t = 0.12 s.).
In light of the discontinuous bead formations shown in Figure 4(a), two additional simulations were conducted to ascertain the effect of standoff distance (h_{s}) for equivalent values of U and V (i.e. U = V = 0.10 m/s). As shown in Figure 5(a), for h_{s} = 1.0 D_{N} while the formation of circular disk formations persist, the degree of disjointedness is much less pronounced. For smaller values of h_{s} (i.e. h_{s} = 0.5 D_{N}), shown in Figure 5(b), the bead becomes approximately continuous in accordance with the results found in various previous research [26,40,41].
Figure 5 Effect of variable standoff distance, h_{s} on PLA volume fraction (U = V = 0.10 m/s; T_{e} = 473 K; t = 0.12 s.).
3.2.1 Solidification & Parametric Analysis
To quantify solidification effects, a series of additional, parametric simulations were conducted relative to the evolution of the temperature and apparent viscosity. Specifically, these parametric studies considered variations in extrusion speed (U) and printing speed (V). The apparent viscosity (η), defined as: η = τ/$\dot{\gamma}$, was used (i.e. in lieu of dynamic or kinematic viscosity) in order to include the effects of shear rate ($\dot{\gamma}$). As the extruded PLA fluid begins to solidify, in addition to the expected reduction in temperature, there is also a concomitant increase in the apparent viscosity.
The results pertaining to the variation in U are shown in Figure 6. In particular, four separate simulations were conducted with U ranging from 0.1 m/s to 0.25 m/s and with all other parameters as stated in Table 1. Precise measurements were considered at the location x = 3.7 mm, y = 0.35 mm. As indicated in Figure 6(a), in addition to the expected decrease in temperature that occurs subsequent to the initial rapid increase in temperature (corresponding to the initial PLA extrusion from the nozzle), there also appears a strong dependence on U. In general, it is observed that smaller values of U correspond to lower values of temperature. As indicated, the largest decrease in temperature (~35 K) corresponds to U = 0.10 m/s. This result is reasonable given that the increase of $U$ corresponds to an increase in mass flow rate (i.e., ṁ $\propto$ ρUD_{N}), thus a larger amount of energy is required for cooling purposes.
Figure 6 Effect of variable extrusion speed on temperature and apparent viscosity (V = 0.10 m/s; T_{e} = 473 K; t = 0.12 s.; measured at x = 3.7 mm, y = 0.35 mm).
The results pertaining to η are shown in Figure 6(b). Here η is observed to increase with decreasing values of U. As indicated, a maximum value of η = 5000 Pa·s is achieved for both U = 0.1 m/s and U = 0.15 m/s relatively early on (i.e. at t = ~0.4 seconds in the case of U = 0.1 m/s). This relatively rapid increase in apparent viscosity is a result of the incongruous, circular PLA beads observed in Figure 4 (a and b), which are more susceptible to the convective cooling effects than the continuous bead patterns observed for higher extrusion speeds and flow rates.
Figure 7 shows the results of variable printing speed (V). In this case, three separate simulations were conducted for V ranging from 0.10 m/s to 0.20 m/s (with all other parameters as specified in Table 1). As indicted in Figure 7(a), the local PLA temperature (measured at x = 3.7 mm, y = 0.35 mm) decreases for V = 0.20 m/s, while it remains approximately unvaried for V = 0.15 m/s and V = 0.10 m/s. The observed temperature decrease with increasing nozzle velocity is attributable to the increased time for cooling and solidification that a higher value of UN affords. The measurement location is taken within the PLA melt and is the same for all values of UN. The apparent viscosity (see Figure 7(b)) is shown to achieve a maximum value of 5000 Pa·s at 0.025 seconds for V = 0.20 m/s, while a gradual increase in apparent viscosity (reaching a maximum of approximately 2500 Pa·s at t = 0.12 seconds) is observed for V = 0.15 m/s and V = 0.10 m/s.
Figure 7 Effect of variable printing speed on temperature and apparent viscosity (U = 0.20 m/s; T_{e} = 473 K; t = 0.12 s.; measured at x = 3.7 mm, y = 0.35 mm).
Finally, Figure 8 shows the results of variable extrusion temperature (T_{e}), ranging from 463 K to 473 K (with all other parameters as nominal per Table 1). From Figure 8(a), the local PLA temperature (measured at x = 3.7 mm, y = 0.35 mm) marginally increases for increasing extrusion temperatures, with differences in magnitude nearly equivalent to the difference in extrusion temperature. In contrast, the apparent viscosity (see Figure 8(b)) is shown to increase with decreasing extrusion temperature. Achieving a maximum value of 3480 Pa·s at 0.025 seconds for T_{e} = 463 K.
Figure 8 Effect of variable extrusion temperature on temperature and apparent viscosity (V = 0.10 m/s; t = 0.12 s.; measured at x = 3.7 mm, y = 0.35 mm).
4. Conclusions
In this work we simulated the FDM mesoscale extrusion and deposition process with solidification effects. Restricting the work to a single line scan with a continuous extrusion speed of PLA polymer, the numerical simulations were conducted using COMSOL [28], and incorporated a 2D heat and mass transfer analysis through a multiphase approach. The free surface, liquid/vapor interface was simulated using a fronttracking, levelset method and the melt viscosity assumed both temperature and shear rate dependency. The results focused on the evolving temperature, viscosity, and volume fraction and were cast within a set of parametric studies associated with the printing and extrusion speeds, as well as the extrusion temperature. Some of the principal findings may be summarized as follows:
 FDM Solidification can be effectively modeled using a fronttracking method (i.e. the level set method) in concert with a temperature dependent liquid fraction function (porosity function) and a modified heat capacity.
 Within the sampled parameter space, the continuity of the bead melt was largely dependent on the extrusion speed, with discontinuities observed for U ≤ 1.5V. This was further accompanied by an increase in solidification, as quantified by the significant increase in apparent viscosity for U ≤ 1.5V.
 With respect to printing speed, it was observed that a significant increase in apparent viscosity occurred for V ≥ U. Like the extrusion speed, this also corresponded in a discontinuous melt pattern.
Author Contributions
Dr. Allen was responsible for the model development, Dr. Riveros oversaw the work and provided technical direction, and Drs. Freeman and Beckman assisted with the technical review and writeup, respectively.
Funding
The simulations described and the resulting data presented herein, unless otherwise noted, were funded under PE 0603734, Project T15 "Additive ManufacturingEnhanced AM Modeling and Simulation”, Task 02 managed and executed at the US Army Engineer Research and Development Center. Permission was granted by the Director, Information and Technology Laboratory to publish this information.
Competing Interests
The authors have no relevant financial or nonfinancial interests to disclose.
References
 Turner BN, Strong R, Gold SA. A review of melt extrusion additive manufacturing processes: I. Process design and modeling. Rapid Prototyp J. 2014; 20: 192204. [CrossRef]
 Mohan N, Senthil P, Vinodh S, Jayanth N. A review on composite materials and process parameters optimisation for the fused deposition modelling process. Virtual Phys Prototyp. 2017; 12: 4759. [CrossRef]
 Wickramasinghe S, Do T, Tran P. FDMbased 3D printing of polymer and associated composite: A review on mechanical properties, defects and treatments. Polymers. 2020; 12: 1529. [CrossRef]
 Sun Q, Rizvi GM, Bellehumeur CT, Gu P. Experimental study of the cooling characteristics of polymer filaments in FDM and impact on the mesostructures and properties of prototypes. Proceedings of the Fourteenth Solid Freeform Fabrication (SFF) Symposium; 2003 August 46; Austin, TX, US. Austin, TX, US: The University of Texas.
 Ziemian C, Sharma M, Ziemian S. Anisotropic mechanical properties of ABS parts fabricated by fused deposition modelling. Mech Eng. 2012; 23: 159180. [CrossRef]
 Nikzad M, Masood SH, Sbarski I. Thermomechanical properties of a highly filled polymeric composites for fused deposition modeling. Mater Des. 2011; 32: 34483456. [CrossRef]
 Ramanath HS, Chua CK, Leong KF, Shah KD. Melt flow behaviour of polyεcaprolactone in fused deposition modelling. J Mater Sci Mater Med. 2008; 19: 25412550. [CrossRef]
 Mostafa N, Syed HM, Igor S, Andrew G. A study of melt flow analysis of an ABSIron composite in fused deposition modelling process. Tsinghua Sci Technol. 2009; 14: 2937. [CrossRef]
 Ji LB, Zhou TR. Finite element simulation of temperature field in fused deposition modeling. Adv Mat Res. 2010; 97101: 25852588. [CrossRef]
 Zhang Y, Chou K. A parametric study of part distortions in fused deposition modelling using threedimensional finite element analysis. Proc Inst Mech Eng B J Eng Manuf. 2008; 222: 959968. [CrossRef]
 Ahn D, Kweon JH, Kwon S, Song J, Lee S. Representation of surface roughness in fused deposition modeling. J Mater Process Technol. 2009; 209: 55935600. [CrossRef]
 Boschetto A, Bottini L. Roughness prediction in coupled operations of fused deposition modeling and barrel finishing. J Mater Process Technol. 2015; 219: 181192. [CrossRef]
 Bellehumeur C, Li L, Sun Q, Gu P. Modeling of bond formation between polymer filaments in the fused deposition modeling process. J Manuf Process. 2004; 6: 170178. [CrossRef]
 Sun Q, Rizvi GM, Bellehumeur CT, Gu P. Effect of processing conditions on the bonding quality of FDM polymer filaments. Rapid Prototyp J. 2008; 14: 7280. [CrossRef]
 Costa SF, Duarte FM, Covas JA. Estimation of filament temperature and adhesion development in fused deposition techniques. J Mater Process Technol. 2017; 245: 167179. [CrossRef]
 Lieneke T, Denzer V, Adam GA, Zimmer D. Dimensional tolerances for additive manufacturing: Experimental investigation for fused deposition modeling. Procedia CIRP. 2016; 43: 286291. [CrossRef]
 Bellini A, Güçeri S. Mechanical characterization of parts fabricated using fused deposition modeling. Rapid Prototyp J. 2003; 9: 252264. [CrossRef]
 Ziemian C, Sharma M, Ziemian S. Innovation and intellectual property rights. Mechanical Engineering. London, UK: IntechOpen; 2012.
 Singamneni S, Roychoudhury A, Diegel O, Huang B. Modeling and evaluation of curved layer fused deposition. J Mater Process Technol. 2012; 212: 2735. [CrossRef]
 Domingo Espin M, Puigoriol Forcada JM, Garcia Granada AA, Llumà J, Borros S, Reyes G. Mechanical property characterization and simulation of fused deposition modeling polycarbonate parts. Mater Des. 2015; 83: 670677. [CrossRef]
 Shahrain M, Didier T, Lim GK, Qureshi AJ. Fast deviation simulation for ‘fused deposition modeling’ process. Procedia CIRP. 2016; 43: 327332. [CrossRef]
 Bellini A. Fused deposition of ceramics: A comprehensive experimental, analytical and computational study of material behavior, fabrication process and equipment design. Philadelphia, PA, US: Drexel University; 2002.
 Du J, Wei Z, Wang X, Wang J, Chen Z. An improved fused deposition modeling process for forming largesize thinwalled parts. J Mater Process Technol. 2016; 234: 332341. [CrossRef]
 Xia H, Lu J, Dabiri S, Tryggvason G. Fully resolved numerical simulations of fused deposition modeling. Part I: Fluid flow. Rapid Prototyp J. 2018; 24: 463476. [CrossRef]
 Heller BP, Smith DE, Jack DA. Effects of extrudate swell and nozzle geometry on fiber orientation in fused filament fabrication nozzle flow. Addit Manuf. 2016; 12: 252264. [CrossRef]
 Comminal R, Serdeczny MP, Pedersen DB, Spangenberg J. Numerical modeling of the strand deposition flow in extrusionbased additive manufacturing. Addit Manuf. 2018; 20: 6876. [CrossRef]
 Serdeczny MP, Comminal R, Pedersen DB, Spangenberg J. Experimental validation of a numerical model for the strand shape in material extrusion additive manufacturing. Addit Manuf. 2018; 24: 145153. [CrossRef]
 COMSOL Documentation. COMSOL Multiphysics Reference Manual, version 5.5. Burlington, MA, US: COMSOL Documentation; 2019.
 Seppala JE, Migler KD. Infrared thermography of welding zones produced by polymer extrusion additive manufacturing. Addit Manuf. 2016; 12: 7176. [CrossRef]
 Costa SF, Duarte FM, Covas JA. Thermal conditions affecting heat transfer in FDM/FFE: A contribution towards the numerical modelling of the process: This paper investigates convection, conduction and radiation phenomena in the filament deposition process. Virtual Phys Prototyp. 2015; 10: 3546. [CrossRef]
 Voller VR, Prakash C. A fixed grid numerical modelling methodology for convectiondiffusion mushy region phasechange problems. Int J Heat Mass Transf. 1987; 30: 17091719. [CrossRef]
 Masud A, Kwack J. A stabilized mixed finite element method for the incompressible shearrate dependent nonNewtonian fluids: Variational Multiscale framework and consistent linearization. Comput Methods Appl Mech Eng. 2011; 200: 577596. [CrossRef]
 Kim S. A study of nonNewtonian viscosity and yield stress of blood in a scanning capillarytube rheometer. Philadelphia, PA, US: Drexel University; 2002.
 Farah S, Anderson DG, Langer R. Physical and mechanical properties of PLA, and their functions in widespread applicationsA comprehensive review. Adv Drug Deliv Rev. 2016; 107: 367392. [CrossRef]
 Le Marec PE, Quantin JC, Ferry L, Bénézet JC, Guilbert S, Bergeret A. Modelling of PLA melt rheology and batch mixing energy balance. Eur Polym J. 2014; 60: 273285. [CrossRef]
 Khoo RZ, Ismail H, Chow WS. Thermal and morphological properties of poly (lactic acid)/nanocellulose nanocomposites. Procedia Chem. 2016; 19: 788794. [CrossRef]
 Balani SB, Chabert F, Nassiet V, Cantarel A. Influence of printing parameters on the stability of deposited beads in fused filament fabrication of poly (lactic) acid. Addit Manuf. 2019; 25: 112121. [CrossRef]
 Mehta R, Kumar V, Bhunia H, Upadhyay SN. Synthesis of poly (lactic acid): A review. J Macromol Sci Part C Polym Rev. 2005; 45: 325349. [CrossRef]
 Hu H, Argyropoulos SA. Mathematical modelling of solidification and melting: A review. Model Simul Mat Sci Eng. 1996; 4: 371. [CrossRef]
 Xia H, Lu J, Tryggvason G. A numerical study of the effect of viscoelastic stresses in fused filament fabrication. Comput Methods Appl Mech Eng. 2019; 346: 242259. [CrossRef]
 Mollah MT, Comminal R, Serdeczny MP, Šeta B, Spangenberg J. Computational analysis of yield stress buildup and stability of deposited layers in material extrusion additive manufacturing. Addit Manuf. 2023; 71: 103605. [CrossRef]