Implementation of viscosity and density models for improved numerical analysis of melt flow dynamics in the nozzle during extrusion-based additive manufacturing

Fused Filament Fabrication (FFF) is an Additive Manufacturing (AM) process that builds up a part via layer by layer deposition of polymeric material. The purpose of this study is to implement viscosity and density models for improving the assessment of melt flow behavior inside the nozzle during deposition. Numerical simulations are carried out for different combinations of important process parameters like extrusion velocity Ve, extrusion temperature Te, and filament material (Acrylonitrile Butadiene Styrene (ABS) and Polylactic Acid (PLA)). Cross-Williams–Landel–Ferry (Cross-WLF) viscosity and Pressure–Volume–Temperature (PVT) density models are incorporated to get realistic results. Distribution of printing parameters like pressure, temperature, velocity and viscosity inside the nozzle are observed at steady state and their relationship with the print quality is discussed. Effect of the PVT model on polymer deposition is illustrated by comparing it with deposition considering a constant density. Velocity profiles are obtained for the different cases considered and locations where the flow is fully developed, along the axial distance of the nozzle, are determined and termed as stable zones. A direct correlation between the position of the developed melt flow profile and printing quality is established and the best combination of printing parameters is proposed for ABS and PLA. Extended stable zones are obtained for the polymer melt in the nozzle at Ve = 60 mm/s, Te = 220 °C for ABS and Ve = 30 mm/s and Te = 195 °C for PLA and hence, these can be considered as the optimum values of the printing parameters.


Introduction
Fused Filament Fabrication (FFF) is by far the most common Additive Manufacturing (AM) technology in commercial use. In this advanced manufacturing process, an extruder is used to deposit the material layer by layer on a print bed. The motion of the extruder is controlled by motors that allow it to move in X, Y and Z directions. The material to be printed is fed to the extruder in the form of filament or beads which is melted in the liquefier. The movement of the extruder is controlled by G-code generated by slicing software.

Strength and quality of 3D printed parts
Material extrusion-based Additive Manufacturing (AM) technology is not widely used for industrial applications 1 3 because parts produced by this process lack in strength and quality in terms of dimensional accuracy. Many research works have been done in recent years to identify the factors affecting the strength and quality and most of them relate to the rheology and microstructure of the material used and the process parameters.

Influence of inter-strand bonding on strength
McIlroy and Graham [1] established a relationship between the nucleation in the polymer melt at a molecular level and its cooling time by introducing the concept of flow-enhanced crystallization in the melt dynamics. They showed that at certain low temperatures, the bulk of the deposited polymer cools down while the backbone is still in stretched condition. This phenomenon reduces crystallization time and has a significant impact on the inter-strand bonding, thus affecting the mechanical strength of the print. Ravoori et al. [2] performed thermal conductivity measurements of the extruded strand and used high-speed imaging techniques to observe microstructure of the build path. They observed that thermal conductivity increases with increasing infill percentage, layer height and decreasing extruder speed. An increase in thermal conductivity results in higher bonding of individual strands and therefore higher mechanical strength. Printing at high temperatures could especially affect the bonding between deposited filaments and hence the strength. Flow simulations by Xia et al. [3], considering the Cross-WLF viscosity model, concluded that at high temperatures, the reheated zone depth is higher and it lasts for a longer time, resulting in lower deformation of the printed part. Morales et al. [4] highlighted the effect of interlayer cooling on the strength and quality of parts produced by FFF process. Decreasing the temperature of the deposited layers results in weaker inter layer bonding. But if the layers are not allowed to cool down, the form and shape of the part gets affected and hence the print quality gets deteriorated. Pollard et al. [5] followed an optical measurement methodology to study the effect of temperature fluctuations on inter layer bonding of the extruded material. They observed that during step changes and retraction motions, the bonds weaken due to sudden federate and temperature fluctuations, causing the deposited structure to be mechanically weak at these places. Wang et al. [6] investigated the effect of infill rate, extrusion velocity and layer height on tensile strength. Their study showed that the deposited layer height is the most dominant factor affecting the interlayer bonding and hence the tensile strength of the material. Bellehumeur et al. [7] directly correlated the extrusion temperature with the bonding between extruded filaments, and the mechanical strength, by conducting sintering experiments and observing neck formation. It was majorly influenced by the extrusion temperature and not the envelope temperature that is respsonsible for cooling down the material after extrusion. Moreover, the change in neck radius was insignificant when the extrudate cooled down below 200 °C. This is way higher than the glass transition temperature of ABS material used in this study.

Influence of process parameters on strength and print quality
Duty, Ajinjeru et al. [8] elaborately discussed the print quality in terms of distortions and warpage caused by different levels of printing parameters. They developed an analytical approach, based on the power law viscosity model that could be followed to define the optimal printing conditions. Nuñez et al. [9], through their experimental study established that density was the most influential parameter controlling the part's dimensional precision. Hence, a proper density model that can accurately account for the compressible behavior of the polymer should be considered while modeling material extrusion-based additive manufacturing process. Zhu et al. [10] concluded that the Cross-WLF is the best viscosity model that can be used for plastics because it not only accounts for a wider range of shear rates, but also considers the effect of pressure on viscosity. Their finite element analysis of the warpage of a considered part along with experimental results verified this assumption. Parts produced with material extrusion-based Additive Manufacturing (AM) technologies often have poor quality at positions where the extruder changes printing direction. Study by Ertay et al. [11] showed that this was caused due to irregular deposition patterns and fluctuation of printing parameters like temperature at these positions. They proposed a methodology to improve the quality of the printed parts by synchronizing the extrusion velocity with the tangential printing velocity and adaptively controlling the extrusion temperature. Verma et al. [12] developed a Computational Fluid Dynamic (CFD) model to simulate the melt flow of PLA after extrusion. Higher extrusion velocities resulted in uneven and distorted strand deposition because of the less time available for interlayer cooling. Agassant et al. [13] analytically and numerically modeled the extrusion process and studied the effect of the process parameters and material rheology on the shape and spread of the extrudate. Material spread in the lateral direction decreased when the print head velocity increased. When the extrusion velocity is as low as half the print head velocity, the extrudate loses contact with the print head bottom surface and spreads in the upward direction instead. These results show the importance of material extrusion velocity for achieving dimensionally accurate prints without significant material spread and distortions. In a similar work, Serdeczny et al. [14] observed that the cross-section of the extruded material changes from circular to almost a flat-rectangular one when the velocity ratio (ratio between the print 1 3 head velocity (V) and velocity of material inside the nozzle (U) or the stand-off distance is reduced. However, they did not consider the non-Newtonian effects of the polymer used. The importance of printing head velocity was further elucidated by Prajapati et al. [15] in their analytical and experimental work. Higher print head velocity is required to reduce the temperature loss of the filament in the stand-off region, before the material gets deposited on to the print bed. To achieve high production rate in the extrusion process, Luo, Wang et al. [16] obtained the maximum value of feed rate that could be used for extrusion at a given temperature for ABS material. They also correlated the extrusion speed with the print speed, i.e. the speed at which the print head moves, and inferred that jamming occurs when the extrusion speed is higher than the print speed. The work of Comminal et al. [17] reveals a trajectory planning methodology to control the print quality at the corner when the nozzle takes a 90° turn during material extrusion. Their CFD model considers an incompressible fluid and introduces a non-dimensional acceleration blending factor to control the X and Y velocities of the print head, hence controlling the print quality. Devlin et al. [18] conducted a wide range of experiments and came to the conclusion that increasing the print head speed and the extrusion temperature increases the strength of the FFF printed parts, although print head speed has lesser effect as compared to extrusion temperature. Kousiatza and Karalekas [19] used optical sensors to measure the temperature profiles and strain fields at different layers of the deposited material during fabrication in real time. The strain has significant values during the solidification of the material and remains constant thereafter. They also recorded the variation of temperature with time and layer number.

Numerical and experimental studies on melt dynamics inside the nozzle
Although a lot of research work has been done to study the rheological properties of polymers used in the FFF process, few experimental works observe the rheology of the polymer inside the nozzle. For this purpose, Coogan and Kazmer [20] designed an in-line rheometer and used it to plot the shear rate dependency of viscosity for Polycarbonate material. They also used the Cross-WLF viscosity model to extrapolate the viscosity to a wider range of shear rates. Anderegg, Bryant et al. [21] carried out in-situ temperature and pressure measurements and observed significant variations in temperature as compared to theoretically obtained values from power law viscosity model, thus proving its inefficiency to be a robust viscosity model for this process. Temperature fluctuations increased with increasing the flow rate, showing the flow rate dependency of the temperature. The viscoelastic behavior of polymers used in the Fused Filament Fabrication process was discussed by Mendes et al. [22]. By experimentally obtaining the dimensionless Reynolds number (Re) at different sections of the flow in a microfluidic channel, they concluded that at low values of Re, the polymer showcased a somewhat Newtonian behavior but as the Re was increased, the effects of non-Newtonian behavior, and as a result the shear thinning, were profound. To experimentally observe the velocity profile distribution inside the nozzle, Peng et al. [23] pigmented the filament to be extruded and then observed the deformation of the pigment using an optical microscope. They discovered that the shear rate is small near the axial position and this gives rise to a blunt velocity profile. They measured the temperature distribution as well and concluded that the extrusion process is non-isothermal as opposed to the characteristics of a power law fluid. Based on this framework, Pigeonneau, Xu et al. [24] used the Carreau-Yasuda's viscosity law to model the temperature distribution inside the liquefier, even though they believed that Cross-WLF viscosity model would have been a better candidate. The results of this investigation revealed that formation of distinct temperature profiles was evident from the contour plots and this helped them develop a metric to find the axial distances at which the flow was completely developed i.e. when the polymer melt reached its glass transition temperature. FEA studies, assuming a constant viscosity, conducted by Sa'ude et al. [25] to observe the melt flow behavior of various polymers in the nozzle showed that the melt velocity was higher at the nozzle centre as compared to the walls, giving rise to a velocity gradient. The velocity increased when the melt entered the constricted region of the nozzle.

Purpose of this study
The extensive review of state-of-the-art research in the field of extrusion-based Additive Manufacturing (AM) has revealed that process parameters control the strength, quality, productivity and the overall efficiency of this process. While most of the work done in the past is experimental and analytical in nature, the literature lacks numerical analysis that can be used to predict process outcome. The relevant numerical studies that have been done are based on various assumptions and models that ultimately reduce the robustness and reliability of the results obtained. The most prevalent assumption is that the polymers used as extrudate are Newtonian fluids and have constant or temperature independent density. Moreover, the numerical studies which have been done in the past for simulating the extrudate deposition, assume the extrusion velocity to be constant throughout the whole nozzle outlet section. This assumption ignores the presence of a velocity profile distribution across the nozzle outlet.
The purpose of this study is to conduct numerical simulations to model the melt flow dynamics inside the nozzle and avoid the abovementioned assumptions by considering the Cross-WLF model to account for the non-Newtonian property of the polymers and the Pressure-Volume-Temperature (PVT) model for the density. Implementation of these models results in a velocity profile distribution at the nozzle outlet. A single strand deposition is simulated, taking into account this velocity profile, and its shape is compared with cases where these models are not implemented. ANSYS® Fluent is used to carry out 36 individual simulations to obtain data related to different permutations process parameters. In the end, a unique concept of a stable zone, existing inside the nozzle, is presented and used to recommend the optimal printing conditions for ABS and PLA. Figure 1a shows the internal cross-sectional dimensions of the nozzle. A filament wire of 1.75 mm diameter is considered for this study and hence the diameter of the inlet is taken to be 1.8 mm, considering the geometrical allowances. It is assumed that the walls of the nozzle domain considered in this study are uniformly heated by the heater block. The outlet nozzle diameter (D n ) is 0.4 mm. To reduce computational time, an axisymmetric profile is modeled using ANSYS® SpaceClaim with the axis of the geometry representing the axis of the real FFF extruder. The axisymmetric model allows the use of structured mesh which gives faster convergence and better results compared to a 3D unstructured mesh. To maximize the accuracy of the simulation results, a structured mesh is achieved using the Multi Zone Quad/Tri Method. The 2D representation of mesh with details is given in Fig. 1b. This ensures that the results do not deviate in the small cross-section near the outlet. Element size is set to be 0.02 mm which produces a total of 10 elements in the outlet section and a total of 38,586 nodes and 37,668 elements in the axisymmetric geometry of nozzle. Following boundary conditions are imposed in all the studies.  Also, a constant temperature boundary condition (equal to the extrusion temperature) is imposed to simulate the temperature-controlled heating effect of the heater block. Moreover, a coupled scheme solver is implemented on the model because it provides more accurate results as compared to the other solvers offered by ANSYS® [27] Figure 1c shows the schematic setup of extruder used for single strand deposition tests, highlighting the various printing variables namely average extrusion velocity (V e ), Table Velocity (V t ), nozzle stand-off distance (h) and the nozzle diameter (D n ).

Filament material and melt flow modeling
Polylactic Acid (PLA) and Acrylonitrile Butadiene Styrene (ABS) are common FFF materials that show non-Newtonian and shear thinning characteristics. Considering the variation of viscosity is important in extrusion processes because of the strong shear flow feature of the polymers. To improve the numerical predictions, a viscosity model has to be implemented that can account for this variation. Out of the many viscosity models available for polymeric materials, like the power law and Carreau-Yasuda model, the Cross-WLF model is chosen to be implemented for this study because it accounts for a wider range of shear rates and considers the shear rate as well as temperature dependency of viscosity [10]. The Cross-WLF model is represented by Eq. (1) [26].
-shear rate (1/s). τ * -critical stress level at the transition to shear thinning, determined by curve fitting.
n-power law index in the high shear rate regime, determined by curve fitting.
The shear rate ( ̇ ) is defined as: [28] where D is the tensor representing rate-of-deformation and is given by where u i and u j are the velocities in the i and j directions. The zero-shear viscosity is given below: where, T-temperature (K), T * -reference temperature (K), [29].
The Cross-WLF parameters for the selected grades of ABS and PLA are given in Table 1.
The viscosity is determined by the Cross-WLF [26] relation which takes temperature and shear rate as inputs. Since FLUENT does not have a built-in function for this model of viscosity, it is included in the solver through a User-Defined Function.
In this study, apart from the usual continuity and momentum equations, the energy balance equation is also activated to account for the temperature dependency of viscosity. The specific heat capacity (C p ) and thermal conductivity (k) are considered constant throughout and density (ρ) is modeled using the PVT model [30]. This model expresses the specific volume of the polymer melt as a piece-wise function of Pressure and Temperature. The coefficients of this model are obtained by fitting the experimental data using a non-linear regression.

Selection of extrusion parameters
The major purpose of this study is to correlate the quality of melt deposition with the important process parameters of FFF process. To obtain a good quality of print, a correct combination of these parameters must be identified for the nozzle. Following are the parameters considered for this study: • the nozzle outlet diameter D n ; • the average extrusion velocity V e of the extruded filament at the nozzle outlet; • the extrusion temperature T e .
The nozzle diameter D n is considered to be 0.4 mm and the ratio of the stand-off distance and the nozzle diameter is taken as 0.6 for these simulations.
The experiments are designed keeping in mind the entire range of extrusion velocities and extrusion temperatures viable for ABS and PLA materials. Three average extrusion velocities (30 mm/s, 60 mm/s and 90 mm/s), representing a lower limit, an upper limit and a mid-value are considered so that the entire range can be studied and possible conclusion about the best average extrusion velocity can be proposed. For each average extrusion velocity, three extrusion temperatures are studied, as shown in Table 2, to identify the best combination of these 2 parameters for each material. The results obtained are used to identify response stable zones (mm) of melt, as illustrated in later sections. The resulting plan of experiments considers 18 combinations: 3 levels of V e , 3 levels of T e and two materials. Moreover, simulations are done for 2 cases considering: 1. Cross-WLF model for viscosity and constant density 2. Cross-WLF model for viscosity and PVT model for density [30] This is done because the state of the art research studies [13,14,17] consider constant viscosity and density for their simulations and only a few of them consider a viscosity model but constant density [12,17,24]. Study of the material rheology reveals that viscosity varies in the order of 10 7 to 10 2 Pa-s (from viscosity contour in Fig. 2), and density varies around 10% inside the nozzle. Following this observation, it is necessary to employ best-fit models for these quantities in the simulations. Hence, the abovementioned cases are considered to study the extent of influence these models have on the melt flow. Later in this study, a comparison between the results of these two cases is presented. Also, the velocity and temperature profiles obtained at the nozzle outlet are used to simulate a single strand deposition

Results and discussion
When only Cross-WLF model is applied for numerical modeling of nozzle dynamics in Case 1, the melt properties show a non-linear dependency on shear rate as well as temperature, while the density and heat capacity are considered constant. In Case 2, when both Cross-WLF and PVT models are applied, the PVT model accurately predicts the rate-dependent volumetric variation in the nozzle. These contours are plotted for PLA material for V e = 60 mm/s and T e = 195 °C. While comparing the contour profiles for Case 1 and Case 2, there is no clear difference for melt flow inside the nozzle.
Although the nozzle walls are maintained at T e , the filament core does not reach the melt temperature until it reaches the nozzle neck region, for any considered V e . This is apparent in Fig. 2 where a color contour of the temperature distribution is visualized. In any case, the simulated temperature profile of PLA shows that the polymer melt is not sufficiently heated before it reaches the nozzle neck region and this is confirmed with earlier reported results [24]. Furthermore, it essentially confirms the requirement of higher force for high feed rates or extrusion velocities because of limited temperature transport within the nozzle.
Velocity distribution is also obtained for extrusion of PLA material, considering integration of Cross-WLF and PVT models. Figure 3 shows gradual development of velocity profile at the nozzle neck. This leads to an increase in axial velocity as shown in the magnified contour (Fig. 3b). The velocity profile at the outlet is almost developed for the extrusion setting of V e = 60 mm/s and T e = 195 °C and can be attributed to the complete melting of polymer before going out of the nozzle.
However, Fig. 2 cannot justify the importance of using both Cross-WLF and PVT models for simulating the melt dynamics outside the nozzle because the melt deposited on table may need to consider the convex (almost parabolic shape) nature of velocity profile at the nozzle outlet. Thus, single strand deposition experiments are carried out to understand the combined effect of the viscosity and density models on the numerical predictions. The numerical simulation for the deposition process was carried out in ANSYS FLUENT with a two-phase VOF model. The outlet velocity profile from the nozzle simulation was imposed as the inlet profile of the deposition domain for cases II and III in Fig. 4. An implicit VOF scheme was adopted to ensure a robust solution and better convergence. Here, we assume that the outlet velocity profile (V e profile in model) is one of the key parameters controlling the geometry and quality of deposited material during FFF [23]. A high peak velocity at the outlet profile implies more divergence from average velocity, which gives rise to non-uniform melt deposition during 3D printing. The peak velocity can be lowered by decreasing the printing temperature or using a high viscosity material. The values of the printing parameters used in single strand 3D simulations are material: PLA, V e = 60 mm/s, T e = 195 °C, h = 0.24 mm. Stand-off distance (h), is the gap between the nozzle and the print table, V e is the average extrusion  Figure 4a shows the deposited strand and the positions of inlet and various cross-sectional planes. The cross-sections of the strand for the above said sliced positions with different combinations of model integration have been shown in Fig. 4b. A flattening is observed in the shape for all cases because of the combined effects of swelling and change of the flow direction i.e. the polymer melt leaves the nozzle vertically and gets deposited horizontally on the table. Thus, the approximation of an elliptic shape of the deposited strand is not accurate for the given process parameters. The cross-section of the strand is an elongated rectangular shape with rounded edges for most of the sections [14].
For Case III, velocity profile obtained from 2D nozzle simulations is imposed on the nozzle outlet along with Cross-WLF and PVT models. Strand deposited in this case accurately captures polymer melt velocity profile, viscosity and density variations and the swelling effect [15], giving the most accurate simulation results. To demonstrate the importance of these models and the velocity profile, Cases I is simulated without implementing the velocity profile and Case II is simulated without the PVT model.
The differences between these three cases can be seen in Fig. 5, where the strand height and width are plotted. Figure 5 shows the measured strand width and height as a function of the position along the strand length. The comparison of the dimensions of the strand along the deposited length shows a decreasing trend because of the shrinkage caused due to the increase in density in the strand as the melt temperature decreases. It is also evident from measurements that both strand width and height for Case III are always lower as compared to case I and II because rheological effects, volumetric shrinkage and outlet velocity distribution are simultaneously considered for case III. This improves its numerical prediction. The huge leap at the end of the plots is due to the fact that it represents the starting point of deposition and the temperature is the lowest in this region due to prolonged cooling. Shrinkage due to cooling and high density at this region (based on the PVT law) gives rise to a steep drop in dimensions. The printed strand dimensions obtained in this study are comparable with those obtained from empirical relationships proposed by Hebda et. al. [31], that take material swelling into consideration.

Effect of Reynolds Number on outlet velocity profile
The radial velocity at outlet is a convex (almost parabolic shape) distribution. Area under the velocity distribution Reynolds number is the ratio of inertial forces to viscous forces and its value determines the dominant force, and hence the extent of laminar flow of the extruded material. Figure 6 can be referred to study the velocity front development at different values of Re. The dimensionless number, Re [24] is defined by Eq. (5) The difference between the maximum velocity (V max ) and average velocity (V avg ) is also calculated and used as a response variable in Fig. 6. The effect of V e is dominant on A vp as three distinct levels are seen in Fig. 6. In case of PLA, each line from left to right represents, increasing T e from 185 °C to 205 °C. It can be seen that the effect of T e is marginal on A vp , for a given level of V e . While a higher A vp represents a more developed melt profile, an increase in the difference between V max and V avg corresponds to higher fluctuation of V e and leads to a non-uniform melt deposition. Hence, both of these derived parameters are considered to assess the flow at the outlet and achieve the best quality of print.

Effect of melt viscosity on pressure
In Fig. 7, ΔP, the difference between inlet and extrusion pressure is plotted with respect to the average outlet viscosity (η). Pa-s. The extrusion pressure range obtained in this study is comparable to that obtained using mathematical formulation discussed by Phan et al. [32]. The pressure at inlet varies from 5 to 120 bar depending on the material and V e [33]. This corresponds to an inlet thrust force and torque rating consistent with most stepper motors used in extruders of FFF printers. As η increases, the pressure required to push the melt forward increases due to increasing viscous friction and constriction at nozzle throat, until P reaches the outlet pressure (which is 1 Pa). It is observed from Fig. 2 that a temperature gradient is formed, where the temperature is highest at the nozzle walls and lowest along the centre line due to imposed temperature constrain at the walls. Increasing V e , while keeping T e constant, leads to larger temperature gradient in the nozzle. At high V e , the temperature transport between the wall and the centerline is reduced. It is confirmed by earlier studies [3,5,24] that the print quality can be improved by achieving low temperature gradient at the outlet. If η is maintained between 150 and 250 Pa-s for PLA and 700-1000 Pa-s for ABS, the temperature gradient and ΔP are reduced and this corresponds to improved print quality.

Effect of melt temperature on Prandtl Number
For thermodynamic considerations, Pr the Prandtl number, [24] is defined by Eq. (6), where, C p is the heat capacity per unit volume of the fluid, is dynamic viscosity, and k is thermal conductivity. Thermal properties of PLA and ABS are also provided in Table 1. Pr is the ratio of momentum diffusivity and thermal diffusivity.
The high values of Pr observed in this study represent that heat convection dominates over conduction and as the value of Pr increases, the effect of heat convection prevails. It can also be concluded from the high Pr values that the thermal boundary layer is much thicker in case of ABS than in PLA. According to Eq. (6), Pr number is directly proportional to viscosity. In Fig. 8, Pr for V e = 60 mm/s is lower than 30 mm/s because the shear rate is higher for V e = 60 mm/s and PLA has a shear thinning behavior. For V e = 90 mm/s, the curve has an odd behavior. This is because, at such high extrusion speed, a higher temperature gradient is developed between the melt core and the walls. The existence of such temperature gradient was confirmed by the CFD analysis performed by Phan, Horner et al. [34]. Thus, the average outlet temperature decreases significantly compared to the set T e and this leads to a high viscosity and hence a high value of Pr compared to the case with V e = 60 mm/s. It can be concluded that the effect of temperature starts becoming dominant as V e exceeds 60 mm/s and is the major parameter affecting viscosity as V e reaches 90 mm/s (Fig. 8).

Stable zone identification
Pigeonneau et al. [24] characterized the melt profile development by introducing a parameter L Iso-Tg , which is the axial distance from the inlet where the melt temperature reaches its glass transition temperature. Their study showed that the value of L Iso-Tg is a function of the inlet velocity and a higher velocity results in a cold core area and plugged flow, which is undesirable. Thus, they correlated the process parameters with the melt profile development. In the present study, an alternative approach is followed to optimize printing parameters based on melt viscosity (η) and axial velocity (V) rather than melt temperature in nozzle. A combination of V e and T e Fig. 7 Effect of melt viscosity (η) on extrusion pressure (ΔP) for PLA and ABS for which the longest stable zone is obtained, is considered to be optimal. A stable zone is defined by the region where the melt is homogenized i.e. the slope of d 2 V/dX 2 (X being the axial position) is zero and steady state axial velocity is attained, as shown in Fig. 9a. When slope of first order differential of axial velocity (dV/dX) becomes constant and the slope of second order (d 2 V/dX 2 ) becomes zero, the melt is treated as homogenized and melt viscosity becomes practically constant at lowermost value. Figure 9b shows the variation of melt viscosity η (Pas) with stable zones of PLA, for different printing conditions. Ideally, the average η inside the nozzle and at its outlet should decrease with increasing V e as the considered polymers are non-Newtonian fluids and follow shear thinning behavior [13,35,36]. But, from Fig. 9b, it can be observed that as V e is increased, the stable zone is formed at higher η. This is because, with increasing V e , temperature transport from wall to core of polymer is delayed and the reduction of temperature dependent η slows down. Table 3 shows the variation of stable zone with printing parameters used in this study. The combination of printing parameters, for which the stable zone is the longest, is considered to be optimal and is underlined. Intermediate level of average extrusion velocity (V e = 60 mm/s) for ABS and lower level of average extrusion velocity (V e at 30 mm/s) for PLA show the longest stable zones. Studies referring to mechanical strength and print quality also confirm the abovementioned settings to be optimal for 3D printed parts [13,[37][38][39][40][41].

Conclusions
The study of melt dynamics inside the FFF nozzle is one of the most important areas of research for development of 3D printing technology. The Cross-WLF and PVT models are not yet integrated with the velocity profile at nozzle outlet to perform numerical studies [42]. This study reveals some valuable insights into the flow behavior of most common polymers such as ABS and PLA inside a typical FFF nozzle. Following are the conclusions resulting from this work: 1. Compared to the earlier studies assuming constant density and viscosity [4,13,17], the implementation of Cross-WLF and PVT models give significantly different results for melt deposition, which are comparable with the experimental results [37,38]. This is because the PVT model considers material shrinkage as well as thermal history of melt and the Cross-WLF model considers the dependency of viscosity on shear rate as well as temperature. 2. At the outlet of nozzle, the melt parameters like viscosity, temperature and velocity are not uniform; they vary across the section. A convex (nearly parabolic) profile of velocity is obtained at the nozzle outlet. The peak of this profile increases with increasing wall temperature and decreasing average extrusion speed. This fact is assumed to play a key role in print quality and cross-section of strand gets more elliptical and closer to earlier results [13,43]. 3. At high average extrusion velocities, the effect of the developed temperature gradient on the viscosity of the melt is more profound than the effect of shear thinning behavior of the non-Newtonian polymeric material. This can be explained by the reduced temperature transport between the walls and the nozzle core at high melt velocities and, as a result, a decrease of the average temperature at the outlet. 4. A unique concept of a stable zone existing inside the FFF nozzle is proposed in this study. Melt stability zone is determined where the steady state axial velocity is attained and slope of d 2 V/dX 2 is zero. A robust prediction of the optimal process parameters can be made to obtain the best print quality in FFF process. Based on simulation trials, V e = 60 mm/s, T e = 220 °C for ABS and V e at 30 mm/s, T e = 195 °C for PLA show elongated stable zone.
Funding Open access funding provided by Politecnico di Milano within the CRUI-CARE Agreement. Not applicable.
Availability of data and material Sources of data used for this paper have been included in references.
Code availability ANSYS FLUENT v19.2 Student's License was used for the numerical analysis. Custom codes were written to impose models of material properties-viscosity and density.

Conflict of interest
The authors have no conflicting or competing interests in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.