Non-Newtonian modeling of contact pressure in fused filament fabrication

The nozzle pressure was monitored in a fused filament fabrication process for the printing of high impact polystyrene. The contact pressure, defined as the pressure applied by the newly deposited layer onto the previous layer, is experimentally calculated as the difference between the pressure during printing and open discharge at the same volumetric flow rates. An analytical method for estimating the contact pressure, assuming one-dimensional steady isothermal flow, is derived for the Newtonian, power-law, and Cross model dependence of shear rates. A design of experiments was performed to characterize the contact pressure as a function of the road width, road height, and print speed. Statistical analysis of the results suggests that the contribution of the pressure driven flow is about twice that of the drag flow in determining contact pressure, which together describe about 60% of the variation in the observed contact pressure behavior. Modeling of the elastic and normal stresses at the nozzle orifice explains an additional 30% of the observed behavior, indicating that careful rheological modeling is required to successfully predict contact pressure. © 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1122/8.0000052


I. INTRODUCTION
Fused filament fabrication (FFF, also referred to as material extrusion additive manufacturing and fused deposition modeling or FDM ® ) is a material extrusion process in which the filament feedstock is driven through a hot end to programmatically extrude roads and layers to form complex three-dimensional structures. A salient feature of most FFF processes is the use of the driven solid filament to convey compressive stresses that yield melt pressures in the hot end for the extrusion of the filament. Pressure losses through the hot end and nozzle are related to the geometry of the flow channels, processing temperatures and flow rates, and the material rheology [1,2].
Recently, there has been significant interest in understanding the states and rheology of melts in FFF processes. Serdeczny et al. experimentally characterized and analytically modeled the feeding force of the polymer melt flow through the hot end in FFF, finding that it followed expected isothermal melt rheology until melt temperatures dropped at higher melt flow rates [3]. Build rate limits in FFF including nonisothermal effects have been formally investigated by Hart and co-workers with findings that limits on the traction force exerted on the filament and conduction heat transfer to the filament core are dominating constraints [4]. Osswald et al. modeled the heat transfer of FFF nozzles, finding that the thermal contact resistance between the filament and the hot end, as well as the low thermal conductivity of the filament, also limits the energy transfer from the heater to the filament [5]. As such, the transient heat transfer in hot ends and nozzles is an active area of the current research [3,5,6].
In FFF processes, the cohesion between layers is largely affected by the contact stresses between the incoming road and underlying substrate [7]. In conventional polymer processing, bonding and lamination are achieved by direct pressurization or indirect tensioning. For example, multilayer films are formed between rollers in which the contact pressure at the roller nips causes the polymer to integrally weld [8]. Schach and Creton found that cohesion in laminated layers is closely related to the ratio between the contact time and the reptation time and the Deborah number [9]. Similarly, in the fiber placement process, which is a popular additive manufacturing method for thermoplastic composites, the layers are fusion bonded by a compaction roller [10]. Also, in selective laser sintering of powders, compaction is an important process step [11,12]. However, the only compaction that occurs in the FFF process is caused by the pressure required to force the incoming material into the space between the nozzle and the previous layer. Recently, rheological and heat transfer effects in FFF have been studied to investigate the correlation between the Nusselt and Graetz numbers and estimate the nozzle pressure by monitoring the power used to drive the hot end [13]. While nozzle pressure monitoring is important for viscosity estimation and quality assurance, there is increasing interest in the use of nozzle pressure to estimate the contact pressure between the deposited extrudate and substrate.
The melt pressure formed between the nozzle and the underlying substrate and adjacent roads influences the intimate contact surface area as well as the structural boundary conditions for fusion bonding. The structural boundary conditions, given the impinging flow field, will result in a contact pressure. This contact pressure is not directly controlled but can be influenced based on the process variables such as the polymer viscosity, flow rates, and local geometry between the layers and roads. The contact pressure has been investigated in the course of developing the interlayer contact model. Coogan and Kazmer first defined the contact pressure as the difference between the pressure during printing and the pressure of open discharge at the same volumetric flow rates [14]. In that study, a method was developed to measure the pressure inside the nozzle using a load pin delivering the resultant force to a calibrated, cantilever load cell, enabling real-time pressure measurements [15]. The results showed that higher pressure differences were observed for narrower gaps (smaller layer heights) and higher flow rates [14], which are consistent with engineering physics and intuition.
In real FFF processes, the melt flow is inherently three dimensional and involves complex, nonisothermal melt rheology as well as free and moving boundaries. Mackay discussed some of the non-Newtonian and viscoelastic effects and heat transfer limitations that challenge extrusion operations coming to bear in FFF [16]. For example, it is known that the extruded melt in FFF nozzles swells up to 30% in diameter [7], which is not easily regulated. If the land segment of the nozzle is long enough, the swell might be suppressed to a certain degree due to sufficient relaxation times within the nozzle, thereby reducing the memory effect. Longer lands with reduced swell are constrained in actual application by the induction of the additional pressure loss. Accordingly, compressibility and viscoelastic effects pertaining to die swell will affect the contact pressure. The normal stress development in the nozzle is empirically correlated with shear stress, which is proportional to the pressure loss in the nozzle [17]. Moreover, the internal flow driven by the filament feeding has not been fully identified yet and there are still ongoing studies [18]. Therefore, the nozzle pressure measurement is a versatile method for experimentally analyzing the flow within the nozzle in FFF and subsequent extrusion.
The flow beyond the nozzle orifice greatly affects the surface and weld properties of the extruded materials [19,20]. With respect to flow and pressure in the deposition zone, Agassant et al. recently modeled the spreading of the flow according to viscous non-Newtonian behavior through an approximate shear thinning power-law model [21]. Still, the contact pressure is quite difficult to predict. The flow directly below the nozzle orifice is the most complicated regime in the entire flow. A free boundary tends to first be radially formed with the flow direction changing in the direction of the plate movement. Then, the flow between the nozzle bottom and the substrate follows, which has been treated as a fully developed Poiseuille and Couette flow of a Newtonian fluid with the nonzero pressure gradient [14]. Leaving this regime, the pressure is released while solidification further progresses.
This work experimentally and theoretically investigates the flow between the nozzle and the substrate to estimate the contact pressure. The flow between the nozzle bottom and the substrate has been analytically solved based on the Cross viscosity model. A constant pressure gradient is assumed considering the short flow length. Moreover, the flow is assumed to be fully developed since the high viscosity renders a rapid transition to a steady flow profile. In this work, the pressures during open discharge and extrudate deposition (i.e., printing) are characterized according to a design of experiments (DOE) by the response surface method. An approximate model for the flow and pressure directly below the nozzle is proposed to constitute the contact pressure model. Comparing the experimental measurements with analytical estimates, the statistical validity of the contact pressure model is assessed relative to contributions due to pressure flows, drag flows, elastic stresses, and normal stresses.

A. Methodology
In this study, the printing patterns and conditions were determined according to a DOE. First, a base layer of 0.2 mm was deposited to provide a level substrate relative to the traversal of the print head. Then, there was an open discharge at the same volumetric flow rate to be used during the print (this is referred to as the purge step). Next, the preset pattern was printed on top of the base layer, all while measuring the melt pressures inside the nozzle. The contact pressure was experimentally obtained from the pressure difference, which will be further described in Sec. II G. Non-Newtonian flow analysis was performed to estimate the pressure around the nozzle orifice using similar conditions. To construct a model that can closely predict the contact pressure, statistical corrections are proposed, and the model coefficients were determined by regression analysis using the analytically determined pressure and the observed contact pressure. The workflow performed in the contact pressure modeling is presented in Fig. 1.

B. Variables
The flow field in FFF is greatly affected by the geometric setup and process conditions. Consider the cross section of the flow during deposition as shown in controlled melt temperature, T melt , and one of the most important process parameters, the print speed, S. When the layer height, H, is set, the flow rate, Q, determines the road width, W, as (approximating as a rectilinear flow when W ) H) The nominal shear rate is defined as In this work, S, H, and W are taken as the independent experimental variables. The printing pressure, P print , and the equilibrium open discharge pressure, P open , were measured during the normal printing process and purging (open discharge) at the same volumetric flow rate wherein Q ¼ HWS.
The selection of S, H, and W was proscribed according to a full factorial DOE with the three factors examined at three levels as later detailed.

C. Printer setup
A FFF printer that allows the measurement of temperature and pressure in the nozzle was employed. Figure 3(a) illustrates the overall experimental setup for extrusion control, measurement, and data acquisition; refer to [14] for its design and specifications. The nozzle geometry was also characterized by x-ray computed tomography [15]. The geometry and dimensions required for the analysis are presented in Fig. 3(b). Several thermocouples were installed in the system to measure temperatures of the melt, build plate, and ambient air. In addition, the load cell temperature was measured so that temperature corrections could be applied to the pressure sensor. To verify the feed rate of the filament, an additional encoder was mounted upstream of the feeder gear as indicated. Each run in the DOE had a print time and open discharge time of 10 s. Data were sampled at 250 kHz with down sampling to 1 kHz so that each experimental run provides approximately 20 000 pressure observations from which the steady state printing and open discharge pressure can be reliably ascertained.

D. Filament
An amorphous thermoplastic filament made of a high impact polystyrene (HIPS) was chosen. The filament is denoted as HIPS Natural from eSun (Shenzhen, China). Its diameter was reported to be 2.85 mm with a tolerance of ±0.05 mm, which is sufficiently controlled to avoid significant errors in the flow rate and shear rate. The glass transition temperature was measured as 97.9°C via differential scanning calorimetry (DSC, Q2000, TA Instruments, New Castle, DE) with a ramp rate of 20°C/min; a heat-cool-heat cycle was used with the glass transition temperature in the first heat cycle measured as 104.5°C. The melt viscosity of the filament was measured by a capillary rheometer (Rosand RH10, Malvern Panalytical, Westborough, MA).
The non-Newtonian viscosity of the molten filament can be modeled by the Cross model without the truncation term [22], where η 0 is the zero-shear viscosity, n is the Cross law index, λ is the time constant, and j _ γj is the shear rate that will be defined later in this work. The fitted constants for the filament are presented in Table I [15].

E. Design of experiments
The experiments in this work were performed under steady state conditions with a build plate temperature of 80°C and a nozzle temperature of 250°C. The 80°C plate temperature was set according to standard printing practices for this material to balance the solidification rate with bed adhesion. Three different conditions for each independent variable (W, H, S) were selected. To avoid excessive pressures in the nozzle that could damage the instrumentation at melt pressures greater than 10 MPa, the DOE excluded those cases with volumetric flow rate, Q, above 0.01 ml/s from the full factorial design that would have had 3 3 = 27 runs. The performed runs are listed in Table II together with the apparent shear rates at the nozzle inlet/exit and the nominal shear rate at the gap between nozzle and plate. It is noted that the bulk melt temperature at flow rates below 0.01 ml/s were previously found with an intrusive thermocouple to vary less than 1°C from the set nozzle temperature [23]. Accordingly, the isothermal flow is subsequently assumed in the analysis.

F. Print patterns
The printed road geometries for the DOE are shown in Fig. 4. Because of three different speeds, S, in the DOE with the fixed run time of 10 s, there are three different travel distances of the roads. This work only considers printing of roads on the second layer because the first layer of roads was deposited to ensure a level substrate relative to the extrusion crosshead (x, y) travel. Each road will have different cross sections due to the different flow rates and print speeds. The circles in Fig. 4 represent the location and relative volume of the purge performed at the flow rate and the duration for each road to define the equilibrium open discharge pressure, P open . To prevent possible interference, these purges were manually removed during the continuous printing process.

G. Pressure measurement
The goal of this work is to measure and predict the pressure at the nozzle exit. The pressures observed in the FFF  system during the printing process are influenced by many factors. Apart from the flow resistance outside the nozzle, the flow and thermal conditions inside the nozzle can affect the measured pressure. Still, there may be systematic biases that are difficult to identify. For example, the sensor is placed upstream to the converging channel as can be seen in Fig. 3; the nozzle geometry as characterized by x-ray computed tomography is available [15]. The pressure loss in the nozzle section upstream of the orifice is certainly not zero. Accordingly, this work assumes that these losses and biases are the same during printing and open discharge at the same volumetric flow rates. Therefore, a reasonable pressure at the nozzle exit can be obtained by subtracting the open discharge pressure, P open , from the pressure during printing, P print , to cancel errors such as the entrance losses upstream of the nozzle orifice. Thus, the base assumption is that the flow downstream of the nozzle orifice does not disturb the flow inside the nozzle. As long as the flow rate and temperatures are consistent, the estimated contact pressures should be valid. This will be further described in Sec. III F.

III. FLOW MODEL
A. Drag flow with pressure gradient Consider the melt flow during the FFF printing process per Fig. 2, which shows the cross section of the flow outside  the nozzle considering a control volume attached to the nozzle. The domain of the melt downstream the nozzle orifice can be divided into four subdomains: Ω 0 , Ω 1 , Ω 2 , and Ω 3 . In Ω 1 , the melt flow changes the direction and is subject to high shear, whereas the melt velocity is regarded negligible at Ω 0 due to the static free boundary at @Ω A . At @Ω 1 , the pressure, P 1 , which drives the whole flow together with the plate movement, is specified. This pressure is defined as the contact pressure, P contact , and the details will be described in Sec. III F. This work assumes that a unidirectional steady drag flow with a pressure gradient (Couette-Poiseuille flow) takes place in Ω 2 . At the inlet of Ω 2 , which is denoted as @Ω 2 , a pressure, P 2 , can be specified while the pressure at @Ω 3 is P atm . The pressure, P 2 , represents the pressure required for the flow through Ω 2 . This Couette-Poiseuille flow should dominate the stress field during printing including the contact pressure at the location directly below the nozzle on @Ω B . Subsections III B-III E provide the modeling method in Ω 2 for estimating P 2 on @Ω 2 . In Ω 3 , the pressure will be released and the melt solidifies forming the structure. Other details regarding other domains and boundaries of Fig. 2 are subsequently detailed.

B. Velocity and flow rate
The velocity, V(z), in Fig. 2, as a function of the distance, z, from the build plate can be obtained from the shear rate. It is expressed as Note that _ γ(z) above is a directional quantity. By force balance for a thin fluid element, Assuming a generalized Newtonian, isothermal fluid, the shear stress is Combining Eqs. (4)-(6) yields [14] @p @x ¼ @ @z η @V @z : However, for a non-Newtonian fluid, the above differential equation is not trivial to solve analytically. From Eq. (4), the velocity for the moving substrate case can be written as where the boundary condition is imposed. Integration of Eq. (5) gives where τ 0 ¼ τ(0) is the shear stress at the substrate. The pressure gradient is expressed as where The above expression gives Rewriting Eq. (8) using Eqs. (4) and (13) gives [24] To provide an analytical solution, we introduce an integral function The velocity is then described as The shear rate at the substrate, _ γ 0 , has to satisfy The flow rate per unit width is described as Integration by parts with Eqs. (9) and (17) gives This expression can be written again as To solve the model analytically, a second integral function is introduced As a result, Eq. (20) can be further rewritten as Replacing _ γ 0 in Eq. (16) gives Then, the pressure on @Ω 2 can be approximated by C. Verification in the Newtonian case When η ¼ μ, the integrals by Eqs. (15) and (21) are evaluated as From Eq. (6), is obtained. Then, by Eq. (22), q becomes Moreover, Eq. (16) gives Using Eq. (26), Eq. (28) reduces to As a result, are reached. Substitution of the above shear rates in Eq. (27) gives Describing Eq. (31) again in terms of p x gives where P 2N denotes the Newtonian approximation of the pressure on @Ω 2 . This is the well-known solution to the Navier-Stokes equation and is the same result as the literature [14].

D. Power-law and Cross model
For a fluid following the power-law model, and the stress is expressed as Substituting the above two equations in Eqs. (15) and (21)  results in Note that the velocity for the power-law model is conventionally obtained by directly integrating d(dV/dy) n /dy ¼ p x /K such as in [25], which is not possible for the Cross fluids. With Eqs. (33) and (34), the pressure gradient by Eq. (11) is obtained as where _ γ H and _ γ 0 must satisfy the following velocity and flow rate conditions: Substitution of Eq. (36) into Eq. (38) followed by the solution with Eq. (37) gives _ γ H and _ γ 0 . Then, back substitution into Eq. (36) provides the pressure gradient, which is used to estimate the pressure, P 2P , for the power-law viscosity model.
For a Cross fluid without the truncation term, the shear stress is where Integration in Eqs. (15) and (21) yields where 2 F 1 is the Gaussian or ordinary hypergeometric function and m ¼ 1 À n. Substitution of Eqs. (41) and (42) into Eq. (24) provides the estimated pressure on @Ω 2 , referred to as P 2C for the pressure predicted by the Cross model.

E. Procedures
To determine P 2 by Eq. (24) for any of the three viscosity models, the two wall shear rates and the pressure gradient, _ γ 0 , _ γ H and p x are obtained first by solving three nonlinear algebraic equations at the same time. The first condition originates from the force balance described by Eq. (11), which is described with the three unknown variables as Then, a second condition is obtained from the no slip condition on the bottom surface between the extrudate and substrate according to Eq. (17), which can be written as Then, the flow rate condition by Eq. (23) gives By simultaneously solving Eqs. (43)-(45), _ γ 0 , _ γ H , and p x can be determined; caution should be taken not to be trapped in a trivial solution of p x ¼ 0 and _ γ H ¼ _ γ 0 ¼ S/H. A beneficial aspect of this modeling approach is that the methodology can work for any generalized Newtonian fluid, regardless of the viscosity model.

F. Contact pressure and print pressure
The intimate contact between layers in the FFF process is governed by the pressure, P 1 , incoming to Ω 1 as shown in Fig. 5. We define ΔP 0 as the pressure drop between the sensor and the nozzle orifice during printing and ΔP as the pressure drop between the sensor and the nozzle orifice during open discharge. Thus, the observed pressure during printing is P print ¼ ΔP 0 þ P contact , while the observed pressure during open discharge is just P open ¼ ΔP. As proposed in [14], assuming ΔP 0 ¼ ΔP for the same flow rate, the contact pressure can be experimentally estimated by the difference between the printing and open purge pressure, i.e., Note that ΔP 0 does not equal ΔP when the flow in the nozzle-plate gap affects the upstream intranozzle flow. Otherwise, ΔP 0 ¼ ΔP as long as the rheological conditions are the same in both flows. For example, although the shear history or thermal degradation can vary the material properties, Eq. (46) is still valid since the changes in material properties occur in both P print and P open at the same volumetric flow rates. Consider a virtual boundary, @Ω 0 , shown in Fig. 5, where the magnitude of the velocity is negligible. Note that the flow through @Ω 0 still does not vanish. Ignoring the momentum, elastic effects, and the pressure gradient in the y direction in Ω 1 , the force balance in the x and z directions are P 0 % P 2 and P 0 % P 1 , respectively. As a result, P 2 can be approximated as the estimated contact pressure on @Ω 0 ,P contact . The experimentally observed P contact and modeled P 2 by Eq. (24) are next compared.

A. Pressure measurements
All the print traces in Fig. 4 were printed in a single printing process. The FFF process was conducted continuously to minimize any variability in the environment, material, or process that could affect results. The pressure history was obtained while executing the printing process. As previously described, the cyclic process of purging and printing is repeated for every run in the DOE. Figure 6 shows all the acquired print pressures sampled at 250 kHz during this experiment and down sampled to 1 kHz such that each reported pressure observation is the average of 250 acquired readings in a 1 ms interval. For each DOE run, the open purge and second layer of printing are conducted after the base printing. Figure 7 shows one cycle of data for the DOE's run number 9.
It is observed that the nozzle pressure is typically on the order of 1-2 MPa for this material at a variety of print conditions. During open purging of the melt as shown in Fig. 8, the nozzle pressure exhibits a relatively fast step response with consistent noise and little long-term variation. This response suggests that the flow readily achieves steady state during open discharge. However, the pressure response during road printing suggests a developing flow with free boundaries. In the leading edge, the pressure overshoots to overcome the initial flow resistance. For some cases, such as runs 4-6 and 12, these overshoots are extreme as can be seen in Figs. 6 and 7. Then, the pressure stabilizes and mildly decreases toward the falling edge. To assess the contact pressure, both the open purge and print pressures are averaged in each plateau (after the initial transient) to provide characteristic pressures P open and P print , respectively.
In Fig. 6, the first six runs are for S = 4 m/min, the next eight are for S = 2.5 m/min, and the next nine are for S = 1 m/min. In each group, the flow rate decreases along with the print velocity and thus the pressure does as well. It is observed that the open and print pressures are relatively similar in magnitude as shown in Fig. 6. For example, the open pressures are 98.4% and 94.2% of the print pressures for the run 2 and run 9, respectively.

B. Observation of the contact pressure
The mean values in each segment are presented in Table III. Even after averaging, the open pressure is higher than the print pressure for some of the runs, resulting in negative P contact . Negative contact pressures are sometimes observed in cases with S = 4 m/min and S = 2.5 m/min. The negativity of P contact is especially common for smaller values of W. This result implies that the behavior is not due to an instability or a random phenomenon, but that the drag flow component has lowered the observed print pressures within the nozzle at the location of the load pin. Figure 8 plots P print , P open , and P contact as a function of H, W, and S. Both P print and P open increase with increasing H, W, and S. This behavior is expected as each of H, W, and S contributes to increasing Q that drives increasing pressures. However, P contact decreases with increasing H but increases with increasing W as shown in Fig. 8(b). Note that S does not significantly influence the contact pressure as can be seen in the bottom right subplot of Fig. 8(a). The strong dependence on H in Eq. (32) matches the results. The trends for P contact with H, W, and S also agree with the previous results [14].
The contact pressures shown in Fig. 8(a) and Table III indicate P contact is near zero or negative in nine runs of the DOE. Among these nine cases, eight of them occur when W = 0.35 mm. Only one negative contact pressure is observed for W = 0.5 mm with none for W = 0.65 mm. Thus, the negative value is speculated to be due to narrowness of the printing road relative to the nozzle diameter (W < D) for reasons later modeled and discussed.

C. Pressure estimation
The pressure estimation of P 2 by Eq. (24) is performed for the Newtonian, power-law, and Cross viscosity models. The Newtonian viscosity is specified as η( _ γ N ) by Eq. (3). The viscosities by the power-law and Cross models are plotted in Fig. 9 together with the apparent viscosities of the DOE runs at _ γ N ; the power-law consistency and the index at 250°C are provided in the caption of Fig. 9. With these viscosities, values for the modeled contact pressure,P contact , are obtained by Eq. (24).
The estimated pressures for all the cases are presented along with P contact in Fig. 10 and in supplementary material Fig. S1 [32]. Generally, the Newtonian assumption overesti-matesP contact , whereas the power-law and the Cross model approaches the values of P contact , indicating improved fidelity. It is observed that the Cross model provides slightly lower pressure than the power-law, since the power-law model overestimates the viscosity at low shear rates as shown in the viscosity plot of Fig. 9. While the accuracy of the power-law model seems sufficient for this material, we recommend the use of the Cross model in general application where the existence of a strong Newtonian plateau (e.g., for polycarbonate) is not known a priori.
Because the analytical model does not consider any effects of W, large errors are observed for the runs with narrow widths. Accounting for elastic effects and pressure gradients in the width direction would improve estimations. The observed tendency of the pressure to decrease with increasing thickness as shown in Fig. 8 is consistent with the estimated pressures in Table III. Specifically, for some runs with H = 0.4 mm, the observed and estimated contact pressures are quite close. However, for runs with H = 0.1 mm, discrepancies are evident. Modeling of elastic and three-dimensional effects are later investigated.

A. Validity of the contact pressure
It is suggested that the contact pressure and stress field are influenced by elastic effects. Given a moving substrate, the leading contact point of the melt may be pulled toward the nozzle exit and, in theory and practice, may even detach from the wall of the nozzle's bore. Even before detachment or slip, induced elastic strains can cause a reduction in the compressive normal stress and observed melt pressures.
Three flow regimes are illustrated in Fig. 11. In a stable normal operating mode, the printing bead is expected to be formed as shown in Fig. 11(a). The trailing contact point, C 2 , may also migrate but C 1 is the matter of concern here. In some processes such as slot coating, the upstream region is vacuumed to hold C 1 on the leading edge. By comparison, if the printing speed in FFF increases such that the crosssectional area of the extrudate is less than the cross-sectional area of the nozzle orifice, then C 1 can move toward the nozzle exit as shown in Fig. 11(b) and slip can occur in the inner nozzle wall. When the printing speed is further increased, the flow can be separated from the inner nozzle wall and C 1 can move into the nozzle wall as in Fig. 11(c). Most simulations of this process (e.g., [21,26,27]) present results similar to Figs. 11(a) and 11(b). Such detachment is not always but can be unstable. In Fig. 11(c), the contact point can vary due to small changes in the processing variables. Since no instability was observed in the experiments presented in this work, detachment is not believed to have occurred, although the possibility cannot be excluded.
Wall slip during printing will incur ΔP 0 , ΔP [24]. As aforementioned, the contact pressure measurement in this study relies on the assumption that ΔP 0 ¼ ΔP, the slip will cause P contact inaccurate. Further experimental investigations and modeling could elucidate the issue but, practically, print settings are often chosen to avoid this regime by selecting conditions so that the draw ratio (DR) is around 1 and wall slip does not occur. The draw ratio can be defined by DR ¼ S U , where U is the mean velocity in the nozzle exit section defined by U ¼ Q Let us further describe the flow phenomenon in the gap between the nozzle and substrate. Consider the melt just leaving the nozzle orifice. It will be accelerated to S along a  Table I. Power-law constants are K ¼ 6:04 Â 10 3 Pa s 2Àn and n = 0.405. streamline and the pressure will decrease. This streamline might influence the separation of the contact point, C 2 in Fig. 11(b). As in a plastication screw, a strong drag flow can result in an adverse pressure gradient [28]. This has been also found in previous studies on wire coating and spinning [29,30]. Acceleration and elongation could affect the pressure more strongly than the slip inside the nozzle since the material time in the gap is an order of magnitude larger than that inside the nozzle due to cooling and solidification. In other words, it is suggested that ΔP À ΔP 0 ( jP contact j would be legitimate even when P contact , 0, and the wall slip is feasible. The pressure energy could be reduced while increasing the kinetic and strain energies, which might result in P contact , 0.

B. Models for contact pressure
The predicted pressure on @Ω 2 can be compared to a first linear model (LM1) based on the pressure and drag flow terms of Eq. (32) as where the pressure flow term provides and the drag term provides Here, a 1 and a 2 are the coefficients that should be evaluated by the regression analysis. As is evident from Eq. (32), the ratio of the coefficients in a Newtonian flow should be which will be examined after the analysis of the results. The model here pursues a minimum of residual of the sum of squares between the modeled and observed contact pressures that is expressed as is assumed and the residual is of the form The model can be improved by considering the elastic and three-dimensional effects. So far,P contact ¼ P 2 ¼ P 1 ¼ P 0 has been assumed based on the force equilibrium in Ω 1 : Let us modify the relationship to deduce a more accurate model. The discrepancy between P 2 and P 0 can be estimated by the elastic stress, E 1 , due to volumetric strain in Ω 1 , and the pressure loss in the width direction, P corr , Moreover, the shear rate near @Ω 1 causes a normal stress, N 1 , with no normal stress on @Ω 0 due to negligible shear rate in the vicinity of the free boundary. The force balance in the z direction, also ignoring the momentum, gives Relating Eqs. (52) and (53) by elimination of P 0 gives Recall thatP contact is the estimated contact pressure on @Ω 1 and P 2 is analytically obtained on @Ω 2 by Eq. (24). As illustrated in Fig. 12, the melt can experience significant changes in shape and cross-sectional area at the nozzle exit. The volumetric strain can be expressed in terms of the area ratio of the cross sections before and after the discharge as ln(WH/πR 2 inner ). Because the axial strains would be different according to the directions, E 1 is estimated as The constant, c 0 , can provide a rough correction for elastic errors such as change in orientation from the vertical flow direction in the nozzle to a horizontal flow direction during deposition as well as the surface tensions on the free surfaces. The second and third terms of Eq. (55) represents the Hencky strains in the y-and z directions, respectively. A linear regression for a second linear model (LM2) considering E 1 will determine values of a i and c i that minimize the following residual of the sum of squares between the modeled and observed contact pressures, SS 2 , where i denotes the index for the DOE run and E 1, with a shear rate-dependent material function, ψ 1 , as The flow in Ω 1 is not unidirectional but develops significantly with the x-direction velocity. A suitable material function according to the White-Metzner formulation is of the form [17] where λ 1 is the constant relaxation time. Taking λ 1 ¼ c 3 , Eq. (57) is written again as Consider a third regression model (LM3) that can estimate the contact pressure incorporating the elastic stress by Eq. (55) and the normal stress by Eq. (59). Its least square form is given by where P 2,i is obtained by Eq. (24) with the corresponding W i , H i , and S i . The estimator of Eq. (60) can be further improved by introducing P corr that compensates for the transverse pressure losses due to shear stresses, which are expressed as Then, the regression for this fourth linear model (LM4) is accomplished by minimizing

C. Model estimations
The variance between the predicted and observed contact pressures suggests that the contact pressure cannot be accurately assessed simply by calculatingP contact . Thus, to predict P contact fromP contact , the correction terms provided in Eq. (54) are believed necessary. The objective of the statistical modeling in this section is to explore the topology of the system response and identify potential sources of variation that may be explained by rheological modeling. The intent is primarily to identify statistical significance of the error sources with causal (ab initio) modeling based on material characterization data.
The linear regression of the observed contact pressures according to Eq. (51) has been conducted with the results provided in Table IV. In the first case of linear model one (LM1) without elastic stresses, the ratio of the first two coefficients is This result is very close to the Newtonian theoretical value, −2, implying that an accurate model requires analysis of the elastic stress term. Thus, this term has been implemented in the predictive model described by Eqs. (60) and (62). The accuracy of the estimation has been significantly improved by including the elastic stress term as can be noticed by comparing the R 2 of LM1 and LM2 in Tables IV and V. It is acknowledged that the ratio a 1 a 2 is subject to variation by other factors such as wall slip. As previously discussed, slip can occur within the nozzle and would reduce the melt pressure on @Ω 2 . The wall shear stress estimated for the DOE runs ranges from 25 to 82 kPa, which is smaller that the reported critical slip stress for polystyrene, 120 kPa [31]. As such, slip is not an expected source of the predictive error. Figure 13 compares the estimated and observed pressures. For both the models with (LM2) and without (LM1), the elastic strain are very close to each other with only slight differences with respect to the values of a 1 and a 2 . These linear models correspond with the experimental observations quite well for the plot with respect to H. The trends are also well represented for varying W and S. The analytical estimations do not fully follow the experimental observations, but the results for H with the non-Newtonian viscosities are quite  example, the pressure data for DOE run 2 if a predicted pressure P 2 C ¼ 0:1572 MPa and a calculated contact pressure of P contact ¼ 0:0279 MPa. The shear work in the gap would generate a normal stress at the nozzle orifice resulting in a pressure increase at @Ω 2 . This effect has been accounted for by Eq. (59), giving a dramatic improvement of the correlation in LM3. In addition, it is observed that the model coefficient c 3 does not change much between LM3 and LM4, which suggests that the normal stress term bears substantial physical meaning related to the material relaxation time. Accordingly, the proposed model should be able to estimate the contact pressure within the tested range of W, H, and S. This approach thereby supports the direct comparison of the full analytical model including the Cross model, elastic stresses, normal stresses, and transverse pressure losses. It is observed that the addition of each modeled behavior improves the statistical fidelity; supplementary material Fig. S2 [32] compares all the observed contact pressure data with the estimates of models LM3 and LM4. The trade-off, of course, is that the added modeling complexity increases the characterization requirements as well as computational time for realtime control in ad hoc applications.

VI. CONCLUSIONS
Contact pressure is a critical determinant of the weld strength, with higher contact pressures (corresponding to higher print speeds, wider roads, and thinner layers) providing improved mechanical properties. As such, the modeling and in-line monitoring of the contact pressure is likely to be critical to future process and quality control paradigms. Experimental observations have indicated that the nozzle pressures in FFF are on the order of 1-5 MPa, which are relatively small compared to many conventional extrusion processes with higher flow rates. The comparison of the printing and open purging of melts suggests that the contact pressures are a fraction of the nozzle pressures, typically on the order of 1 MPa or less. Analysis of the flow suggests that the pressure driven flow (Poiseuille) tends to dominate the drag flow (Couette) with a ratio of −2. Statistically as well as experimentally, the obtained ratio of a1 to a2 − 2 suggests that the flow is not much deviant from no-slip flow even when the possibility of slip is accepted. The described solution methodology was implemented for the Newtonian, power-law, and Cross model dependencies of the shear rate with improved pressure prediction for the power-law and Cross models. Still, the analysis of the residual error suggests that elastic and normal stresses at the nozzle orifice are significant and need to be modeled in order to explain 90% or more of the observed variation in the contact pressure.
The results suggest that on-line characterization of the melt with instrumented processes remains a strong candidate for process and quality control. The suggested approach is to characterize the melt rheology for a given material feedstock as a function of temperature and flow rates for at least two nozzles with varying lengths and diameters. After fitting a suitable rheological model, the contact pressure can be reasonably estimated from the road dimensions and processing conditions. In the future, it is likely that FFF preprocessors could consider rheology models in the path planning to predict nozzle pressures and part properties with the goal to enable single part quality assurance.