Numerical Study on the Temperature-Dependent Viscosity E ﬀ ect on the Strand Shape in Extrusion-Based Additive Manufacturing

: A numerical model that incorporates temperature-dependent non-Newtonian viscosity was developed to simulate the extrusion process in extrusion-based additive manufacturing. Agreement with the experimental data was achieved by simulating a polylactic acid melt ﬂow as a non-isothermal power law ﬂuid using experimentally ﬁtted parameters for polylactic acid. The model was used to investigate the temperature e ﬀ ect on the ﬂow behavior, the cross-sectional area, and the uniformity of the extruded strand. OpenFOAM, an open source simulation tool based on the ﬁnite volume method, was used to perform the simulations. A computational module for solving the equations of non-isothermal multiphase ﬂows was also developed to simulate the extrusion process under a small gap condition where the gap between the nozzle and the substrate surface is smaller than the nozzle diameter. Comparison of the strand shapes obtained from our model with isothermal Newtonian simulation, and experimental data conﬁrms that our model improves the agreement with the experimental data. The result shows that the cross-sectional area of the extruded strand is sensitive to the temperature-dependent viscosity, especially in the small gap condition which has recently increased in popularity. Our numerical investigation was able to show nozzle temperature e ﬀ ects on the strand shape and surface topography which previously had been investigated and observed empirically only. discrepancy validation, B.B. and J.P.; formal analysis, B.B.; investigation, B.B., M.S., and J.P.; resources, L.M.; data curation, B.B.; writing—original draft preparation, B.B.; writing—review and editing, M.S., L.M., M.L. and J.P.; visualization, B.B.; supervision, M.L. and J.P.; project administration, M.L. and J.P.; funding acquisition, M.L. and J.P.


Introduction
Extrusion-based additive manufacturing, also called fused deposition modeling (FDM) or fused filament fabrication (FFF), is a type of additive manufacturing (AM) process which fabricates a 3D part by building up extruded material layer by layer. A common choice of raw material for FDM is thermoplastic polymers such as polylactic acid which is used in this study. A typical FDM process consists of the following steps: (1) A thermoplastic filament is fed into the printer head, which consists of a liquefier and a nozzle. (2) The filament is heated in the liquefier to just above its melting temperature to reach a semi-molten state. (3) The molten filament is extruded through the nozzle. (4) The extruded polymer is deposited on a substrate or on a previously extruded strand. (5) The deposited strand cools from different process conditions obtained from our simulations are compared with the experimental observations and the isothermal Newtonian simulation results. In Sections 3.3 and 3.4, the effects of the nozzle temperature on the strand shape and the surface roughness are investigated. Section 4 summarizes the findings and describes the future works.

Governing Equations
The extrusion and deposition of a polymer melt flow is simulated as a free surface flow with temperature and shear rate dependent viscosity. The following governing equations can describe the flow behaviors in both liquid and gas phases: where ρ, u, p and g are density, fluid velocity vector, pressure, and gravitational acceleration, respectively. The gas phase is assumed as air and the liquid phase as a PLA melt. More specifically, we chose the properties of a PLA filament used in the work by Phan et al. [24] (see the next section). The stress tensor, τ, is described by a power-law model: where K and n are power-law model consistency and index, respectively. The velocity gradient tensor is defined as ∆ ≡ ∇u + ∇u T . We used the Arrhenius relation for the consistency of a shear thinning fluid, K, to describe the temperature effect on the polymer viscosity according to a study on PLA rheology [24]: where T is temperature in Kelvin and K 0 is the value of K at a reference temperature, T 0 . The variation of K with T is governed by the parameter, f, which is obtained by fitting experimental data. The temperature distribution in the PLA flow is solved by the energy equation: ∂(ρc P T) ∂t + ∇·(uρc P T) = k∇ 2 T where c p and k are specific heat capacity and thermal conductivity, respectively. We assume that there is no heat generation or viscous dissipation since it was found that its contribution to a 1-K temperature rise is less than 5% (personal communication with Dr. M. Mackay, the author of Ref [24,26]). Our simulation uses the volume of fluid (VOF) approach to simulate the two-phase system. In VOF, as an Eulerian approach, a color function at each position x and time t is defined as: φ(x, t) = 1 for x ∈ D l at time t 0 for x ∈ D g at time t (6) where D l and D g are domains of liquid and gas, respectively [27]. If the color function is averaged over the volume of a cell, V cell , the liquid fraction of the cell is calculated as: ρ(x, t) = ρ l φ(x, t) + ρ g (1 − φ(x, t)) (8) By replacing the density in the continuity Equation (1) with the density defined in Equation (8), another governing equation for the liquid fraction of the cell is derived:

Numerical Model Solver
We performed the CFD simulation utilizing OpenFOAM, an open source CFD software based on a finite volume method. Various modules, based on object-oriented programming, for solving different equations for various CFD problems are available online. For velocity-pressure coupling (solving Equations (1) and (2)), the PIMPLE algorithm was adapted. The PIMPLE algorithm is a combination of PISO (Pressure Implicit with Splitting of Operators) and SIMPLE (Semi-Implicit Method for Pressure-Linked Equations). The PISO algorithm works well in solving transient problems while the SIMPLE algorithm works for steady-state problems. Thus, the PIMPLE algorithm can solve steady-state problems as well as transient problems and can do so using larger time-steps than the PISO algorithm. Further information regarding the PIMPLE algorithm can be found in [28].
In OpenFOAM, users can modify existing modules or develop their own modules for their own problems. In our study, we developed a new solver for non-isothermal free surface flows by improving an existing solver called "InterFoam" which is used for isothermal, incompressible, immiscible, two-phase fluids based on VOF. The new solver accounts for the non-isothermal nature of the extrusion process by including the Arrhenius model for the power law model consistency (Equation (4)) and the energy equation (Equation (5)). It is aimed to simulate transient free surface problems where transport properties vary with temperature. To accomplish this, the energy equation was coupled with the continuity (Equation (1)) and the momentum equations (Equation (2)). Similar to "InterFoam", both phases are considered incompressible. It was assumed that the dependence of liquid viscosity on temperature follows Equation (4). The validation of our new solver is shown in Appendix A.

System Description
Figure 1 depicts the simulated system which consists of a cylindrical nozzle, a substrate plate, and the surrounding atmosphere. The nozzle has an inner diameter of d = 0.4mm and outer diameter of 2d. We chose the cylindrical nozzle for the purpose of comparing with the published experimental data by Serdeczny et al. [15].The gap distance between the nozzle exit and the substrate plate is g; we are particularly interested in the small gap condition g/d = 0.4, which is a condition where the isothermal Newtonian simulation showed a considerable discrepancy from the experimental data [15]. A viscous polymer melt inside the nozzle with an average flow velocity (the z-component of u is averaged over a cross-section of the nozzle) of <u> = 0.02 m/s is extruded through the nozzle. In this simulation, instead of moving the nozzle horizontally, the nozzle is stationary while the substrate plate moves horizontally at three different velocities, (U = 0.02, 0.015 and 0.01 m/s). The no-slip boundary condition is imposed on the nozzle wall as well as the substrate surface. The liquid-air interface has the boundary condition of equal stress and velocity on each phase. The air velocity at the top of the surrounding atmosphere domain is set to 0. The distance between the top of the surrounding atmosphere and the substrate plate was chosen as 0.76 mm, which does not affect the result of the simulated strand shape while maintaining a reasonable computational load. As shown in Figure 2, a hexahedral grid was used for meshing the 3D simulation domain which resulted in 242,000 cells. Velocity profiles simulated in finer meshes did not give noticeable differences.  Temperatures inside the nozzle and in the surrounding air were initially set to 463-503 and 298 K, respectively. Nozzle temperatures were kept at the inlet temperature by imposing a no heat flux boundary condition on the nozzle wall. Heat flux from the extruded liquid surface was continuously transferred to the air phase by conduction because the bulk air phase is stationary. The polymer melt flow is assumed to have the properties of a PLA melt, as reported in the literature [23,24]. The density, heat capacity, and thermal conductivity of PLA were set to 1300 kg/m 3 , 1641 J/kg⋅K and 0.2 W/m⋅K, respectively. The polymer viscosity data was fitted to Equation (3) as well as Equation (4) for the shear rate and temperature dependencies. The fitting parameters were found to be n = 0.8, K0 = 5.87 × 10 8 Pa⋅s 0.8 , T0 = 298 K, and f = 1.11 × 10 4 ( Figure 3). Based on an experimental investigation of PLA melts, the temperature dependence of K can be described by Equation (4) in the temperature range right above the glass transition temperature of PLA (423K) [25]. Therefore, we extrapolated the fit in Figure 3 near the glass transition temperature. It should be noted that, for the  Temperatures inside the nozzle and in the surrounding air were initially set to 463-503 and 298 K, respectively. Nozzle temperatures were kept at the inlet temperature by imposing a no heat flux boundary condition on the nozzle wall. Heat flux from the extruded liquid surface was continuously transferred to the air phase by conduction because the bulk air phase is stationary. The polymer melt flow is assumed to have the properties of a PLA melt, as reported in the literature [23,24]. The density, heat capacity, and thermal conductivity of PLA were set to 1300 kg/m 3 , 1641 J/kg⋅K and 0.2 W/m⋅K, respectively. The polymer viscosity data was fitted to Equation (3) as well as Equation (4) for the shear rate and temperature dependencies. The fitting parameters were found to be n = 0.8, K0 = 5.87 × 10 8 Pa⋅s 0.8 , T0 = 298 K, and f = 1.11 × 10 4 ( Figure 3). Based on an experimental investigation of PLA melts, the temperature dependence of K can be described by Equation (4) in the temperature range right above the glass transition temperature of PLA (423K) [25]. Therefore, we extrapolated the fit in Figure 3 near the glass transition temperature. It should be noted that, for the Temperatures inside the nozzle and in the surrounding air were initially set to 463-503 and 298 K, respectively. Nozzle temperatures were kept at the inlet temperature by imposing a no heat flux boundary condition on the nozzle wall. Heat flux from the extruded liquid surface was continuously transferred to the air phase by conduction because the bulk air phase is stationary. The polymer melt flow is assumed to have the properties of a PLA melt, as reported in the literature [23,24]. The density, heat capacity, and thermal conductivity of PLA were set to 1300 kg/m 3 , 1641 J/kg·K and 0.2 W/m·K, respectively. The polymer viscosity data was fitted to Equation (3) as well as Equation (4) for the shear rate and temperature dependencies. The fitting parameters were found to be n = 0.8, K 0 = 5.87 × 10 8 Pa·s 0.8 , T 0 = 298 K, and f = 1.11 × 10 4 ( Figure 3). Based on an experimental investigation of PLA melts, the temperature dependence of K can be described by Equation (4) in the temperature range right above the glass transition temperature of PLA (423 K) [25]. Therefore, we extrapolated the fit in Figure 3 near the glass transition temperature. It should be noted that, for the simulation, the minimum and maximum viscosities of the liquid phase were set to 0.001 and 10 6 Pa·s to avoid infinite viscosity, particularly when the melt inside the nozzle starts to move. For the simulation of a fluid of which viscosity is so high that it behaves like a solid, its maximum viscosity must be chosen carefully. If the maximum viscosity is set too high, a much smaller time step is required. In contrast, a very low maximum viscosity, such as 5000 Pa⋅s set by Xia et al. [23], may result in a seemingly more liquid-like strand shape than the actual one. In order to handle the high viscosity, we utilized the non-constant time step option of OpenFOAM. At the beginning of the simulation, the range of the time step size was set as 10 −8 −10 −7 s. After approximately 10 time steps, the time step size reached 10 −7 s. At the moment when the liquid surface in the nozzle started to move, the maximum time step size was changed to 10 −6 s. This method enabled us to reduce the computational load while maintaining reasonable and convergent numerical results. Under these simulation conditions, the polymer melt part maintained its liquid state since the temperature of the liquid part did not fall to the melting temperature of PLA (318 K). As such, solidification can be neglected in the simulation. In future simulations, however, such as a multilayer deposition simulation, we plan to include solidification as in the work by Liu et al. [22] which did not consider the temperature-dependent viscosity in the liquid phase.

Simulation of the Extrusion Process
The process where a PLA melt is extruded through the nozzle and deposited on the substrate surface is simulated in this procedure, as described in Section 2. The simulated time-evolution of the liquid fraction (α ≥ 0.5: red) is presented in Figure 4. Note that the color fraction 0 ≤ α < 0.5 is set to the atmosphere (blue). This color definition of each phase gives a clear surface contour to compare to the actual shape of the PLA melt flow [15]. At t = 0, a PLA melt in the nozzle with T = 473 K (for atmosphere, T = 298 K) starts flowing out of the nozzle. As time lapses, the PLA melt is extruded with <u> = 0.02 m/s and subsequently is being deposited along the substrate surface which is moving at U = 0.02 m/s. For the simulation of a fluid of which viscosity is so high that it behaves like a solid, its maximum viscosity must be chosen carefully. If the maximum viscosity is set too high, a much smaller time step is required. In contrast, a very low maximum viscosity, such as 5000 Pa·s set by Xia et al. [23], may result in a seemingly more liquid-like strand shape than the actual one. In order to handle the high viscosity, we utilized the non-constant time step option of OpenFOAM. At the beginning of the simulation, the range of the time step size was set as 10 −8 −10 −7 s. After approximately 10 time steps, the time step size reached 10 −7 s. At the moment when the liquid surface in the nozzle started to move, the maximum time step size was changed to 10 −6 s. This method enabled us to reduce the computational load while maintaining reasonable and convergent numerical results. Under these simulation conditions, the polymer melt part maintained its liquid state since the temperature of the liquid part did not fall to the melting temperature of PLA (318 K). As such, solidification can be neglected in the simulation. In future simulations, however, such as a multilayer deposition simulation, we plan to include solidification as in the work by Liu et al. [22] which did not consider the temperature-dependent viscosity in the liquid phase.

Simulation of the Extrusion Process
The process where a PLA melt is extruded through the nozzle and deposited on the substrate surface is simulated in this procedure, as described in Section 2. The simulated time-evolution of the liquid fraction (α ≥ 0.5: red) is presented in Figure 4. Note that the color fraction 0 ≤ α < 0.5 is set to the atmosphere (blue). This color definition of each phase gives a clear surface contour to compare to the actual shape of the PLA melt flow [15]. At t = 0, a PLA melt in the nozzle with T = 473 K (for atmosphere, T = 298 K) starts flowing out of the nozzle. As time lapses, the PLA melt is extruded with <u> = 0.02 m/s and subsequently is being deposited along the substrate surface which is moving at U = 0.02 m/s. Figure 5 presents the simulation results of the velocity magnitude distributions in the extruded strand as well as in the nozzle under the same condition as Figure 4. In Figure 5 (top), as the extruded liquid is deposited away from the nozzle, the velocity distribution in the deposition flow becomes almost uniform due to the free surface effect (the stress on the liquid surface becomes negligibly small). This transition of the flow pattern in the moving substrate direction indicates that a shear rate dependent property, such as the shear thinning viscosity, is dominant near the nozzle but becomes weaker in the deposition flow. The flow in the nozzle behaves like a pressure-driven flow in a cylinder which has the centerline velocity close to 2<u> of a Newtonian fluid and a maximum shear rate on the inner nozzle surface (Figure 5 bottom). The simulated velocity profile also shows good agreement with the corresponding analytical expression for a power-law fluid.      As soon as the hot melt in the nozzle (473 K) is extruded through the nozzle, it is exposed to the colder atmosphere (298 K), and its temperature distribution becomes non-uniform indicating that the temperature-dependent viscosity has also changed. weaker in the deposition flow. The flow in the nozzle behaves like a pressure-driven flow in a cylinder which has the centerline velocity close to 2<u> of a Newtonian fluid and a maximum shear rate on the inner nozzle surface (Figure 5 bottom). The simulated velocity profile also shows good agreement with the corresponding analytical expression for a power-law fluid. Figure 6 presents the simulated temperature distribution of an extruded strand under the same condition in Figures 4 and 5. As soon as the hot melt in the nozzle (473K) is extruded through the nozzle, it is exposed to the colder atmosphere (298K), and its temperature distribution becomes nonuniform indicating that the temperature-dependent viscosity has also changed.  Figure 7 shows the viscosity distribution in the extruded strand. As shown in Figures 6 and 7, the non-uniform temperature distribution directly affects the viscosity distribution which is expected to affect the flow behavior in the strand and the strand shape [8].

The Effect of the Viscosity Model on the Relation between Strand Shape and U/<u>
Serdeczny et al. [15] studied the relation between the strand shape and the process conditions. They experimentally measured the width, W, and the thickness, H, of a strand extruded with various conditions of g/d and U/<u>. Their CFD simulation assumed the polymer melt to be an isothermal Newtonian fluid and was able to match the experimentally observed relation for a large gap of g/d = 1. However, it showed a discrepancy for a smaller gap, g/d = 0.4. We hypothesized that the discrepancy is due to the non-uniform viscosity and investigated the strand shape and U/<u> for g/d = 0.4 using our numerical model based on a non-isothermal power-law fluid. As demonstrated in Figure 8, W and H are measured at the cross-section of the liquid fraction of an extruded strand (at x = d).   Figure 7 shows the viscosity distribution in the extruded strand. As shown in Figures 6 and 7, the non-uniform temperature distribution directly affects the viscosity distribution which is expected to affect the flow behavior in the strand and the strand shape [8]. weaker in the deposition flow. The flow in the nozzle behaves like a pressure-driven flow in a cylinder which has the centerline velocity close to 2<u> of a Newtonian fluid and a maximum shear rate on the inner nozzle surface (Figure 5 bottom). The simulated velocity profile also shows good agreement with the corresponding analytical expression for a power-law fluid. Figure 6 presents the simulated temperature distribution of an extruded strand under the same condition in Figures 4 and 5. As soon as the hot melt in the nozzle (473K) is extruded through the nozzle, it is exposed to the colder atmosphere (298K), and its temperature distribution becomes nonuniform indicating that the temperature-dependent viscosity has also changed.  Figure 7 shows the viscosity distribution in the extruded strand. As shown in Figures 6 and 7, the non-uniform temperature distribution directly affects the viscosity distribution which is expected to affect the flow behavior in the strand and the strand shape [8].

The Effect of the Viscosity Model on the Relation between Strand Shape and U/<u>
Serdeczny et al. [15] studied the relation between the strand shape and the process conditions. They experimentally measured the width, W, and the thickness, H, of a strand extruded with various conditions of g/d and U/<u>. Their CFD simulation assumed the polymer melt to be an isothermal Newtonian fluid and was able to match the experimentally observed relation for a large gap of g/d = 1. However, it showed a discrepancy for a smaller gap, g/d = 0.4. We hypothesized that the discrepancy is due to the non-uniform viscosity and investigated the strand shape and U/<u> for g/d = 0.4 using our numerical model based on a non-isothermal power-law fluid. As demonstrated in Figure 8, W and H are measured at the cross-section of the liquid fraction of an extruded strand (at x = d).

The Effect of the Viscosity Model on the Relation between Strand Shape and U/<u>
Serdeczny et al. [15] studied the relation between the strand shape and the process conditions. They experimentally measured the width, W, and the thickness, H, of a strand extruded with various conditions of g/d and U/<u>. Their CFD simulation assumed the polymer melt to be an isothermal Newtonian fluid and was able to match the experimentally observed relation for a large gap of g/d = 1. However, it showed a discrepancy for a smaller gap, g/d = 0.4. We hypothesized that the discrepancy is due to the non-uniform viscosity and investigated the strand shape and U/<u> for g/d = 0.4 using our numerical model based on a non-isothermal power-law fluid. As demonstrated in Figure 8, W and H are measured at the cross-section of the liquid fraction of an extruded strand (at x = d). weaker in the deposition flow. The flow in the nozzle behaves like a pressure-driven flow in a cylinder which has the centerline velocity close to 2<u> of a Newtonian fluid and a maximum shear rate on the inner nozzle surface (Figure 5 bottom). The simulated velocity profile also shows good agreement with the corresponding analytical expression for a power-law fluid. Figure 6 presents the simulated temperature distribution of an extruded strand under the same condition in Figures 4 and 5. As soon as the hot melt in the nozzle (473K) is extruded through the nozzle, it is exposed to the colder atmosphere (298K), and its temperature distribution becomes nonuniform indicating that the temperature-dependent viscosity has also changed.  Figure 7 shows the viscosity distribution in the extruded strand. As shown in Figures 6 and 7, the non-uniform temperature distribution directly affects the viscosity distribution which is expected to affect the flow behavior in the strand and the strand shape [8].

The Effect of the Viscosity Model on the Relation between Strand Shape and U/<u>
Serdeczny et al. [15] studied the relation between the strand shape and the process conditions. They experimentally measured the width, W, and the thickness, H, of a strand extruded with various conditions of g/d and U/<u>. Their CFD simulation assumed the polymer melt to be an isothermal Newtonian fluid and was able to match the experimentally observed relation for a large gap of g/d = 1. However, it showed a discrepancy for a smaller gap, g/d = 0.4. We hypothesized that the discrepancy is due to the non-uniform viscosity and investigated the strand shape and U/<u> for g/d = 0.4 using our numerical model based on a non-isothermal power-law fluid. As demonstrated in Figure 8, W and H are measured at the cross-section of the liquid fraction of an extruded strand (at x = d).   Figure 9 shows the relations between H/d and U/<u> obtained from various simulation methods and the study by Serdeczny et al. [15]. The overall qualitative trend is that H/d decreases with increasing U/<u> (as U gets slower than <u>, the strand becomes thicker). The result from the isothermal Newtonian simulation (with the same approach as Comminal et al. [14] and Serdeczny et al. [15]) shows the largest quantitative discrepancy from the experimental data. However, the result from our non-isothermal simulations show improved agreement with the experimental data. The non-isothermal power law simulation results in the smallest discrepancy from the experiment. We also performed a non-isothermal Newtonian simulation (set n = 1 in Equation (3)) for comparison. The improvement by the non-isothermal Newtonian simulation is not as much as that of the non-isothermal power-law simulation. J. Manuf. Mater. Process. 2020, 4, x FOR PEER REVIEW 9 of 15 Figure 9 shows the relations between H/d and U/<u> obtained from various simulation methods and the study by Serdeczny et al. [15]. The overall qualitative trend is that H/d decreases with increasing U/<u> (as U gets slower than <u>, the strand becomes thicker). The result from the isothermal Newtonian simulation (with the same approach as Comminal et al. [14] and Serdeczny et al. [15]) shows the largest quantitative discrepancy from the experimental data. However, the result from our non-isothermal simulations show improved agreement with the experimental data. The non-isothermal power law simulation results in the smallest discrepancy from the experiment. We also performed a non-isothermal Newtonian simulation (set n = 1 in Equation (3)) for comparison. The improvement by the non-isothermal Newtonian simulation is not as much as that of the nonisothermal power-law simulation.    Figure 10 shows the relationship between W/d and U/<u> obtained from the same simulations and experiments as Figure 9. Again, our non-isothermal power-law simulation shows the most improved agreement with the experimental data and the non-isothermal Newtonian simulation results in limited improvement. The differences shown in Figures 9 and 10 according to different models indicate that both the effect of temperature and shear rate contributes to the strand shape. As shown in Figures 5-7, the temperature and the velocity profiles in the extruded strand are non-uniform, which affects the viscosity, resulting in the shape of the cross section. Figure 5 (top) showed that the shear rate in the strand diminishes along the distance from the nozzle. However, the formation of the cross-sectional shape happens mainly near the nozzle where the liquid fraction has higher fluidity and shear rate than the deposited part away from the nozzle. Figure 11 directly compares the cross-sections of an extruded strand obtained from various simulations and experiment. The cross-section obtained from the isothermal Newtonian simulation is wider than the experimental measurement and this work's simulation. This can be explained by considering the temperature and viscosity profiles seen in Figures 6 and 7. The surface of a strand has a lower temperature and higher viscosity than its inside. The exterior's higher viscosity has low fluidity which limits the spreading of the strand. Therefore, the width should be smaller than that of the isothermal Newtonian simulation where a constant viscosity of 1000 Pa·s is imposed. We believethat the quantitative agreement with the experimental data is improved because the temperature-dependent power law viscosity model is fitted to the experimental data measured from the PLA commonly used in the FDM, and the time step size is properly adjusted to handle very high viscosity. However, there still remains the quantitative discrepancy from the experimental data. Direct experimental measurement rather than using data from literature as well as including viscoelasticity is expected to be able to improve the quantitative agreement [29].    Figure 10 shows the relationship between W/d and U/<u> obtained from the same simulations and experiments as Figure 9. Again, our non-isothermal power-law simulation shows the most improved agreement with the experimental data and the non-isothermal Newtonian simulation results in limited improvement. The differences shown in Figures 9 and 10 according to different models indicate that both the effect of temperature and shear rate contributes to the strand shape. As shown in Figures 5-7, the temperature and the velocity profiles in the extruded strand are nonuniform, which affects the viscosity, resulting in the shape of the cross section. Figure 5 (top) showed that the shear rate in the strand diminishes along the distance from the nozzle. However, the formation of the cross-sectional shape happens mainly near the nozzle where the liquid fraction has higher fluidity and shear rate than the deposited part away from the nozzle. Figure 11 directly compares the cross-sections of an extruded strand obtained from various simulations and experiment. The cross-section obtained from the isothermal Newtonian simulation is wider than the experimental measurement and this work's simulation. This can be explained by considering the temperature and viscosity profiles seen in Figures 6 and 7. The surface of a strand has a lower temperature and higher viscosity than its inside. The exterior's higher viscosity has low fluidity which limits the spreading of the strand. Therefore, the width should be smaller than that of the isothermal Newtonian simulation where a constant viscosity of 1000 Pa⋅s is imposed. We believethat the quantitative agreement with the experimental data is improved because the temperature-dependent power law viscosity model is fitted to the experimental data measured from the PLA commonly used in the FDM, and the time step size is properly adjusted to handle very high viscosity. However, there still remains the quantitative discrepancy from the experimental data. Direct experimental measurement rather than using data from literature as well as including viscoelasticity is expected to be able to improve the quantitative agreement [29].

The Effect of the Nozzle Temperature on the Strand Shape
As shown in the previous section, the temperature is an important factor in accurately modeling the strand shape. We performed the simulation with different temperatures of the PLA in the nozzle to learn more of the temperature effect on the strand shape. Figure 12 shows the strand height and width obtained from the simulation with three different nozzle temperatures and a processing condition of g/d = 0.4 and U/<u> = 1. As the temperature increases, H/d decreases and W/d increases. The polymer flow inside the nozzle has the assigned temperature. As it is extruded into the air and onto the substrate surface, the flow is cooled to the ambient temperature (298 K in this simulation) and the viscosity significantly rises. Higher fluid temperatures in the nozzle result in a higher temperature and a lower viscosity of the deposited fluid which spreads more (larger W/d) resulting in a wider strand. Since the cross-sectional area of the strand must remain constant, H/d shows the opposite trend. Figure 11. Comparison of the cross-sections obtained from simulations and experiment (The images were adapted from Serdeczny et al. [15] by permission from Elsevier).

The Effect of the Nozzle Temperature on the Strand Shape
As shown in the previous section, the temperature is an important factor in accurately modeling the strand shape. We performed the simulation with different temperatures of the PLA in the nozzle to learn more of the temperature effect on the strand shape. Figure 12 shows the strand height and width obtained from the simulation with three different nozzle temperatures and a processing condition of g/d = 0.4 and U/<u> = 1. As the temperature increases, H/d decreases and W/d increases. The polymer flow inside the nozzle has the assigned temperature. As it is extruded into the air and onto the substrate surface, the flow is cooled to the ambient temperature (298 K in this simulation) and the viscosity significantly rises. Higher fluid temperatures in the nozzle result in a higher temperature and a lower viscosity of the deposited fluid which spreads more (larger W/d) resulting in a wider strand. Since the cross-sectional area of the strand must remain constant, H/d shows the opposite trend.

The Effect of the Nozzle Temperature on the Strand Uniformity
One of the problems of FDM is due to the surface roughness of the final product [30]. The process requires a finalizing stage so that resolution of the final product meets customers' expectations. One of the reasons for the surface roughness of the final product is the non-uniformity of layers. In other words, the thickness and width of the final product is not constant and changes along the layers. In previous studies, layer uniformity had been adjusted by introducing additives [31][32][33][34]. One study also considered the temperature effect on the strand uniformity: the surface roughness of the final product is reduced with increasing nozzle temperature [34]. We performed a simulation using our numerical model to investigate the temperature effect on surface uniformity. Figure 13 shows a simulation result of the strand width change along the strand axis (distance from the nozzle) with three different fluid temperatures in the nozzle. The results show the width becomes more uniform with higher temperatures. Standard deviations at each temperature are measured as 0.144, 0.0635, and 0.0319 mm for 463, 483, and 503 K, respectively. Our simulation can show the relation of surface uniformity and nozzle temperature similar to the experimental observation in [34]. This observation can be explained by the viscositytemperature relationship discussed in the previous sections. Higher temperatures in the nozzle lead

The Effect of the Nozzle Temperature on the Strand Uniformity
One of the problems of FDM is due to the surface roughness of the final product [30]. The process requires a finalizing stage so that resolution of the final product meets customers' expectations. One of the reasons for the surface roughness of the final product is the non-uniformity of layers. In other words, the thickness and width of the final product is not constant and changes along the layers. In previous studies, layer uniformity had been adjusted by introducing additives [31][32][33][34]. One study also considered the temperature effect on the strand uniformity: the surface roughness of the final product is reduced with increasing nozzle temperature [34]. We performed a simulation using our numerical model to investigate the temperature effect on surface uniformity. Figure 13 shows a simulation result of the strand width change along the strand axis (distance from the nozzle) with three different fluid temperatures in the nozzle.

The Effect of the Nozzle Temperature on the Strand Uniformity
One of the problems of FDM is due to the surface roughness of the final product [30]. The process requires a finalizing stage so that resolution of the final product meets customers' expectations. One of the reasons for the surface roughness of the final product is the non-uniformity of layers. In other words, the thickness and width of the final product is not constant and changes along the layers. In previous studies, layer uniformity had been adjusted by introducing additives [31][32][33][34]. One study also considered the temperature effect on the strand uniformity: the surface roughness of the final product is reduced with increasing nozzle temperature [34]. We performed a simulation using our numerical model to investigate the temperature effect on surface uniformity. Figure 13 shows a simulation result of the strand width change along the strand axis (distance from the nozzle) with three different fluid temperatures in the nozzle. The results show the width becomes more uniform with higher temperatures. Standard deviations at each temperature are measured as 0.144, 0.0635, and 0.0319 mm for 463, 483, and 503 K, respectively. Our simulation can show the relation of surface uniformity and nozzle temperature similar to the experimental observation in [34]. This observation can be explained by the viscositytemperature relationship discussed in the previous sections. Higher temperatures in the nozzle lead  The results show the width becomes more uniform with higher temperatures. Standard deviations at each temperature are measured as 0.144, 0.0635, and 0.0319 mm for 463, 483, and 503 K, respectively. Our simulation can show the relation of surface uniformity and nozzle temperature similar to the experimental observation in [34]. This observation can be explained by the viscosity-temperature relationship discussed in the previous sections. Higher temperatures in the nozzle lead to a more fluidic material being deposited on the substrate. This fluid will spread more and has a larger chance of obtaining a uniform surface before cooling to the point of solid-like viscosity. Although our current model can show the temperature-dependent uniformity, there are many rheological factors to be considered in the extrusion process in FDM for controlling the surface properties. For example, shark skin and die swelling are attributed to viscoelasticity and affect the surface uniformity of printed products [26]. Therefore, the addition of viscoelasticity to the numerical model will be the next step in further investigation.

Conclusions
We used OpenFOAM to simulate the extrusion process of FDM considering the temperature and shear-dependent viscosity which was fitted from experimental data. We compared our simulated strand shape extruded under a small gap condition and straight nozzle with published experimental data under the same conditions. The comparison showed that temperature-dependent viscosity is effective in improving the agreement with experimental observation: As the surface of the extruded strand becomes cooler, the higher viscosity surface hinders the strand from spreading across the substrate which is not observed in the isothermal Newtonian simulation. The parametric studies with various nozzle temperatures also showed that non-isothermal flow behavior is significant in strand shape and surface uniformity: The less viscous liquid portion with the higher nozzle temperature makes the extruded strand spread more across the substrate with a higher strand uniformity.
In our particular choice of geometry and process, due to the comparison with the available published data, the temperature and shear rate-dependent viscosity model was the major contribution to the improved prediction of extruded strand shape. The addition of viscoelasticity to the numerical model is ongoing for further improving its fit to experimental data and its predictive capabilities in other geometries and processes, such as welding of extruded strands. Funding: This work was funded by the Intellectual Systems Center at the Missouri University of Science and Technology.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
The ability of the new solver to calculate the temperature distribution in multiphase fluids was validated by comparing the results of a simple benchmark simulation using COMSOL Multiphysics, a popular commercial CFD software. For the benchmark simulation, a two-phase system with transient heat transfer was chosen. Detailed geometry of the two-phase system is depicted in Figure A1. The liquid phase is a hypothesized material with ρ = 800 kg/m 3 , η = 0.002 Pa·s, c p = 1800 J/kg·K and k = 13 W/m·K. The temperature of the entire system was initially set to T = 298 K. At t ≥ 0, all the wall temperatures are maintained at 298 K except the wall at the bottom of the liquid domain which was set to 500 K. The gas phase is set as air with 1 atm. Dimensionless temperature was defined as: Dimensionless temperature distribution in the liquid as a function of the distance from the bottom wall is plotted for both COMSOL and OpenFOAM at t = 1 min ( Figure A2). Figure A2 shows that the simulation results from both OpenFOAM and COMSOL are almost identical. Therefore, we deemed our model valid for describing the heat transfer in the open surface flow problems.