Mesoscale Modeling of Extrusion and Solidification During Material Extrusion Additive Manufacturing

Jeffrey B. Allen ^1,2,*, Guillermo A. Riveros ¹, Ivan P. Beckman ¹, Elton L. Freeman ¹

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

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

1. Introduction

Of the various forms of additive manufacturing (AM), extrusion-based 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 fiber-reinforced 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 thermo-mechanical 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 two-dimensional 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 front-tracking, or level-set 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 front-tracking 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 non-viscoelastic 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 non-Newtonian 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 multi-phase approach. The free surface, liquid/vapor interface will be simulated using a front-tracking, level-set 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 thermos-physical sample properties, was not performed in this work for purposes of scope limitations. Post-processing 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:

1. The geometry is assumed two dimensional with coordinates (x, y).
2. The surface tension coefficient (γ) is assumed constant, and therefore the Marangoni effect is neglected.
3. 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].
4. 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₁(T)).
5. A modified heat capacity is used within the energy equation ($C_{p}^{eq}$) which accounts for the latent heat of fusion.
6. 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}T-h(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 front-tracking method (in this case the level-set method) to track the solid-liquid-gas interface. As opposed to other front-tracking methods, the level-set 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 re-initialization, 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}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 )\\ 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₁(T), and c₁ and c₂ are constants (with c₁ and c₂ 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₁/c₂. If c₁ is sufficiently large (with c₂ 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₁(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{(T-T_{\mathrm{m}})^{2}}{(T_{1}-T_{s})^{2}}\right)}}{\sqrt{\pi(T_{1}-T_{s})^{2}}} \tag{9} \]

Where T₁ = 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 non-Newtonian, 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{n-1}{2}} \tag{11} \]

Where μ₀ 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 μ₀. For intermediate shear rates (i.e., $\dot{\gamma}$ > 1/λ), the fluid behaves as a Power-law fluid (i.e. μ = K (∂u/∂y) ^n-1), 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 stand-off distance (h_s) and nozzle diameter (D_N) are specified in Table 1. As indicated, the nozzle translates along the x-direction 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 outward-directed, 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 Thermo-physical 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 (λ, μ₀, 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₀ = 283 K subject to a fixed temperature of 253 K at its left-most 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 thermo-physical 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 stand-off 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 stand-off 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 multi-phase approach. The free surface, liquid/vapor interface was simulated using a front-tracking, level-set 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:

1. FDM Solidification can be effectively modeled using a front-tracking method (i.e. the level set method) in concert with a temperature dependent liquid fraction function (porosity function) and a modified heat capacity. 
2. 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.
3. 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 write-up, respectively.

Funding

The simulations described and the resulting data presented herein, unless otherwise noted, were funded under PE 0603734, Project T15 "Additive Manufacturing-Enhanced 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 non-financial interests to disclose.