Rheology Effects on Predicted Fiber Orientation and Elastic Properties in Large Scale Polymer Composite Additive Manufacturing

Short fiber-reinforced polymers have recently been introduced to large-scale additive manufacturing to improve the mechanical performances of printed-parts. As the short fiber polymer composite is extruded and deposited on a moving platform, velocity gradients within the melt orientate the suspended fibers, and the final orientation directly affects material properties in the solidified extrudate. This paper numerically evaluates melt rheology effects on predicted fiber orientation and elastic properties of printed-composites in three steps. First, the steady-state isothermal axisymmetric nozzle melt flow is computed, which includes the prediction of die swell just outside the nozzle exit. Simulations are performed with ANSYS-Polyflow, where we consider the effect of various rheology models on the computed outcomes. Here, we include Newtonian, generalized Newtonian, and viscoelastic rheology models to represent the melt flow. Fiber orientation is computed using Advani–Tucker fiber orientation tensors. Finally, elastic properties in the extrudate are evaluated based from predicted fiber orientation distributions. Calculations show that the Phan–Thien–Tanner (PTT) model yields the lowest fiber principal alignment among considered rheology models. Furthermore, the cross section averaged elastic properties indicate a strong transversely isotropic behavior in these composites, where generalized Newtonian models yield higher principal Young’s modulus, while the viscoelastic fluid models result in higher shear moduli.


Introduction
Recently, extrusion-based Additive Manufacturing (AM) (otherwise known as fused filament fabrication or fused deposition modeling) has moved rapidly from small scale rapid prototyping to the manufacture of large-scale parts and tooling such as the Big Area Additive Manufacturing (BAAM) system developed by Oak Ridge National Laboratories (ORNL) (Oak Ridge, TN, USA) [1].Typical thermoplastic polymer extrusion-based AM is a process where polymer feedstock materials are melted and deposited on a heated platform, layer-by-layer, to form three-dimensional (3D) objects [2].To achieve a relatively high dimensional accuracy and superior mechanical performances in large-scale parts, carbon fiber filled polymers are employed.Duty et al. show that adding short carbon fibers into the neat Acrylonitrile Butadiene Styrene (ABS) polymer yields a composite with improved elastic properties, especially along the printing direction, and less distortion in the printed part following the bead deposition process of the BAAM system [3,4].
A key factor in the polymer composite deposition AM process is the flow-induced fiber orientation within the printed composite (cf. Figure 1) since material properties of solidified parts depend on the fiber alignment within the printed bead [5].Therefore, the prediction of fiber orientation during the polymer extrusion process, and the subsequent evaluation of mechanical properties in the short fiber polymer composite extrudate is of great importance.
Evans et al. [6] and Lipscomb et al. [7] considered a fully coupled approach where the motion of suspended fibers depend on the flow field, and the fiber orientation influences flow kinematics, typically through the suspension melt viscosity.A computationally more efficient approach that is often employed is a one-way weakly coupled formation that ignores the effect of the fiber suspension on viscosity.This approach has been effective in applications having shear dominant narrow gap flows such as injection and compression modelling simulations [8][9][10], and is the approach we use in this study.Fiber orientation studies in polymer composite AM applications have recently become of interest.Nixon et al. [11] simulated fiber orientation in three Fused Deposition Modeling (FDM) nozzle geometries (convergent, straight and divergent) using Moldflow (Moldflow Corporation, Framingham, MA, USA) and the Folgar-Tucker Isotropic Rotary Diffusion (IRD) model [12].Their work, which ignored die swell, showed that a converging geometry yielded the highest principal fiber alignment and the divergent geometry resulted in the lowest.Additionally, at the exit of the straight and the converging nozzle, a higher alignment was predicted near the center than at the edge, unlike the experimental result reported by Kunc [13].Heller et al. [14] computed the fiber orientation tensors in a conventional small scale FDM nozzle and extruded filament.In their work, die swell was computed by minimizing the integrated normal stress on the free surface using COMSOL Multiphysics (Comsol, Inc., Burlington, MA, USA).Their approach modeled the molten polymer as an isothermal Newtonian fluid in a creeping flow, and assumed an axisymmetric velocity field.Orientation tensors (c.f.e.g., Advani and Tucker [15]) were computed along streamlines within the flow domain from velocity and velocity gradient information.Their results showed that fiber alignment reached its peak at the outer edge of the nozzle, and then decreased towards the center of the flow.
Extrudate swell occurs in many extrusion-based polymer processing applications and is known to be highly influenced by the non-Newtonian behavior of the melt.Crochet et al. [16] theoretically analyzed the die swell of an upper-convected Maxwell fluid based on the mixed finite element method for fluids with implicit constitutive equations.Luo and Tanner [17] applied the Streamline Finite Element Method (SFEM) to the die swell problem, which avoided the numerical instability in high Weissenberg number problems.Luo and Mitsoulis [18] extended the SFEM by adding in a particle-tracking scheme along the streamlines with a Picard iterative scheme.Béraudo et al. [19] applied a finite-element-based method to investigate the extrudate swell of Linear Low-Density Fiber orientation studies in polymer composite AM applications have recently become of interest.Nixon et al. [11] simulated fiber orientation in three Fused Deposition Modeling (FDM) nozzle geometries (convergent, straight and divergent) using Moldflow (Moldflow Corporation, Framingham, MA, USA) and the Folgar-Tucker Isotropic Rotary Diffusion (IRD) model [12].Their work, which ignored die swell, showed that a converging geometry yielded the highest principal fiber alignment and the divergent geometry resulted in the lowest.Additionally, at the exit of the straight and the converging nozzle, a higher alignment was predicted near the center than at the edge, unlike the experimental result reported by Kunc [13].Heller et al. [14] computed the fiber orientation tensors in a conventional small scale FDM nozzle and extruded filament.In their work, die swell was computed by minimizing the integrated normal stress on the free surface using COMSOL Multiphysics (Comsol, Inc., Burlington, MA, USA).Their approach modeled the molten polymer as an isothermal Newtonian fluid in a creeping flow, and assumed an axisymmetric velocity field.Orientation tensors (cf.e.g., Advani and Tucker [15]) were computed along streamlines within the flow domain from velocity and velocity gradient information.Their results showed that fiber alignment reached its peak at the outer edge of the nozzle, and then decreased towards the center of the flow.
Extrudate swell occurs in many extrusion-based polymer processing applications and is known to be highly influenced by the non-Newtonian behavior of the melt.Crochet et al. [16] theoretically analyzed the die swell of an upper-convected Maxwell fluid based on the mixed finite element method for fluids with implicit constitutive equations.Luo and Tanner [17] applied the Streamline Finite Element Method (SFEM) to the die swell problem, which avoided the numerical instability in high Weissenberg number problems.Luo and Mitsoulis [18] extended the SFEM by adding in a particle-tracking scheme along the streamlines with a Picard iterative scheme.Béraudo et al. [19] applied a finite-element-based method to investigate the extrudate swell of Linear Low-Density Polyethylene (LLDPE) and Low-Density Polyethylene (LDPE) melts using a multi-mode Phan-Thien-Tanner (PTT) model [20].Their approach provided an accurate die swell prediction for die geometries of a 2D slit die and a 2D axisymmetric capillary die in low and intermediate shear conditions.Ganvir et al. [21] applied an Arbitrary Lagrangian Eulerian (ALE) algorithm to calculate the extrudate free surface, which enabled the die swell simulations to be performed in both steady state and transient problems.Alternatively, Limtrakarn et al. [22] employed the Simplified Viscoelastic (SV) model implemented in ANSYS-Polyflow (ANSYS, Inc., Canonsburg, PA, USA) to predict die swell of a 3D circular die flow of LDPE.A good agreement between the numerical results and the experimental data was achieved.Clemeur et al. [23] found that the SV model was a cost-effective approach for evaluating the flow-viscoelasticity as compared to conventional viscoelastic fluid flow models including the Oldroyd-B and the Phan-Thien-Tanner models.This paper presents a numerical approach to study the effect of assumed polymer melt rheology on predicting fiber orientation and elastic properties of short fiber polymer composites extruded in large-scale AM.The weakly coupled formation is used to compute the fiber orientation within the polymer melt flow where we first obtain the flow kinematic in an isothermal axisymmetric large-scale AM nozzle.Our flow model is created in two dimensions and solved with the finite element suite ANSYS-Polyflow (version 17.1, ANSYS, Inc.) [24], and includes melt flow within the nozzle in addition to a short section of post-nozzle extrudate, which enables the prediction of die swell at the nozzle exit.We consider a Newtonian fluid model, a Power law model, a Carreau-Yasuda model, a multi-mode Phan-Thien-Tanner (PTT) model, and a Simplified Viscoelastic (SV) model in separate flow simulations.Secondly, the fiber orientation along streamlines within the flow domain is computed from the velocity field computed within the melt flow domain.The Advani-Tucker fiber orientation tensor evaluation equation [15] and the Folgar-Tucker Isotropic Rotary Diffusion (IRD) model [12] are employed to solve the fiber orientation problem.In addition, the Orthotropic Closure (ORT) [5] is used to address the closure problem encountered in the fiber orientation computation.Finally, elastic properties are computed from fiber orientation predictions for each melt rheology using the Tandon-Wang approach with fiber orientation averaging [25,26].This paper uses computational methods alone to gain useful insight into the effect of assumed rheology on properties of an extruded composite bead.It does not include the difficult, if not impossible, experimental procedures that would be required to validate these results.

Materials
In this study, we consider the material rheology of Acrylonitrile Butadiene Styrene (ABS) fabricated by the PolyOne Corporation (Avon Lake, OH, USA).The rheological properties, including the complex shear viscosity η, storage shear modulus G and loss shear modulus G , are measured using a HAAKE MARS 40 rheometer (Thermo Fisher Scientific, Waltham, MA, USA) at 210 • C. Once the experimental data is obtained, curve-fitting for the various rheology models is performed using ANSYS-Polymat (version 17.1, ANSYS, Inc.) [27].
The experimental data appears in Figure 2, which shows apparent shear shinning behavior as expected for ABS.The shear thinning behavior of polymer melts can be expressed through various Generalized Newtonian Fluid (GNF) models.Here, we consider the Power law model written as [28] η and Carreau-Yasuda model given as [28] η( . In the above, K is the consistency index, n is the power-law index, η ∞ is the infinite-shear-rate viscosity, η 0 is the zero-shear-rate viscosity, and κ is the natural time.In Equation (2), the constant a controls the transition from the Newtonian plateau to the Power-law region in the Carreau-Yasuda model.For measured data appearing in Figure 2, the fitted rheology model data for the Power law model is K = 16761 Pa•s n and n = 0.4503.For the Carreau-Yasuda model, we obtained a , η ∞ = 0, µ = 0.3333 s −1 , n = 0.000001455, and a = 0.2398.For comparison, we also consider the Newtonian fluid model, also appearing in Figure 2, having a viscosity µ = 3200 Pa•s, which corresponds to a shear rate of ~30 s −1 from our measured rheology.This shear rate was used to determine a Newtonian viscosity value based on a typical shear rate of 30~40 s −1 given by Duty et al. [29] for Big Area Additive Manufacturing (BAAM) systems.
We also consider a multi-mode PTT rheology model, which is a differential-type viscoelastic fluid model.The exponential form of the PTT model is expressed as [20] exp ελ η 1 tr(T 1 ) with and and Here, T 1 and η 1 are the stress tensor and the viscosity component associated with the viscoelasticity, D is the strain rate tensor, v is the velocity tensor, λ is the mode relaxation time, η is the mode viscosity, ξ controls the shear viscosity behavior, and ε controls the elongational behavior.The fitting results in Figure 3 are in good agreement with the experimental rheology data.The fitted parameters of the PTT model appear in Table 1.
Alternatively, ANSYS-Polyflow includes the Simplified Viscoelastic (SV) model that reduces computational expense when predicting die swell in viscoelastic flows.In the SV formulation, it is understood that extrudate swell in polymer extrusion is associated with the first normal stress difference in the fluid.Hence, the SV model extends the Generalized Newtonian Fluid (GNF) model, where the total stress tensor is given as [24] where the off-diagonal terms given as η( .
γ are the shear stress components.In this form, η( .γ) is expressed by a typical generalized Newtonian model, and .γ is the magnitude of the strain rate tensor D (cf.Equation ( 6)).In the above, Ψµ( .

χ)
. χ represents the first normal stress component, in which µ( .χ) is described in a similar fashion as is done for the shear strain rate.In Equation ( 7), .χ is the specialized viscoelastic variable, which is evaluated with the transport equation where θ .
γ is the relaxation time of the melt which controls the development of the extrudate swell diameter once the melt flow exits the nozzle.In addition, Ψ appearing in Equation ( 7) is an artificial weighting factor, which controls the swelling enhancement versus the input flow rate [24].To use Equation ( 7), a description of the shear viscosity is required.Since the Power law model exhibits an unbounded viscosity at a near-zero shear rate, it is not realistic, especially when free surface prediction is of primary interest.Hence, we employ the Carreau-Yasuda law to represent the behavior of shear viscosity.It has been shown that defining an independent law for the normal stress viscosity increases computational cost, but does not greatly enhance the accuracy of die swell calculation.Therefore, the same Carreau-Yasuda model form is used to describe the first normal stress viscosity term.
The SV model is an empirical construction defined in terms of θ and Ψ, each of which is typically defined to obtain a known flow domain property such as die swell.In our simulation, different sets of parameters θ and Ψ are attempted and values are selected to provide a predicted die swell profile that is in good agreement with results obtained using the PTT viscoelastic fluid model.In this study, θ and Ψ are thus defined as 0.26 and 0.47, respectively.

Flow Kinematics and Die Swell Evaluation
We use ANSYS-Polyflow [24] to evaluate the flow kinematics in the polymer melt flow domain based on conservation of momentum and conservation of mass where we have assumed an isothermal incompressible fluid represents the polymer melt as often appeared in extrusion die flow numerical studies [19,21,30].In the above, p is the pressure, T is the total stress tensor, f is the body force, ρ is the density of the fluid and a is acceleration.Note that non-isothermal effects such as the temperature gradients within the free extrudate or viscous heating in the melt flow may result in nonuniform melt rheology properties such as K in Equation ( 1), which is often addressed with an Arrhenius-type temperature dependency (e.g., the Williams-Landel-Ferry equation) [28].Consequently, we understand that our isothermal assumption may yield some inaccuracy in the predicted numerical data.However, it is expected, as suggested by others [19, 21,30], that any temperature-related variation in the rheology properties would not significantly alter the trends that appear in our result section.Moreover, for a Generalized Newtonian Fluid (GNF) model, the total stress tensor may be written as and for a viscoelastic fluid model where T 1 and η 1 are the non-viscous contributions, and T 2 and η 2 are related to the viscous effect of the flow and η is the total viscosity.The geometry of the flow domain in our study is based on the large-scale Additive Manufacturing Strangpresse Model-19 extruder nozzle appearing in Figure 4.In addition, we include a 1-inch section of free extrudate beyond the nozzle exit in the simulation to capture die swell.Due to the axisymmetry of the nozzle geometry and assumed flow, we are able to simplify the flow domain as a 2D axisymmetric model, which saves significant computational expense.Furthermore, our axisymmetric assumption ignores any swirling motion in the flow that may result from the extruder screw.Since creeping flow is assumed, the inertia contribution "ρa" in Equation ( 9) is ignored in our simulations.In addition, since only a short section of the free extrudate material is considered, gravitational effects are ignored, which results in the body force f in Equation (9) being zero as well.

Flow Kinematics and Die Swell Evaluation
We use ANSYS-Polyflow [24] to evaluate the flow kinematics in the polymer melt flow domain based on conservation of momentum and conservation of mass where we have assumed an isothermal incompressible fluid represents the polymer melt as often appeared in extrusion die flow numerical studies [19, 21,30].In the above, p is the pressure, is the total stress tensor, is the body force, ρ is the density of the fluid and is acceleration.Note that non-isothermal effects such as the temperature gradients within the free extrudate or viscous heating in the melt flow may result in nonuniform melt rheology properties such as K in Equation ( 1), which is often addressed with an Arrhenius-type temperature dependency (e.g., the Williams-Landel-Ferry equation) [28].Consequently, we understand that our isothermal assumption may yield some inaccuracy in the predicted numerical data.However, it is expected, as suggested by others [19, 21,30], that any temperature-related variation in the rheology properties would not significantly alter the trends that appear in our result section.
Moreover, for a Generalized Newtonian Fluid (GNF) model, the total stress tensor may be written as and for a viscoelastic fluid model with where and η are the non-viscous contributions, and and η are related to the viscous effect of the flow and η is the total viscosity.
The geometry of the flow domain in our study is based on the large-scale Additive Manufacturing Strangpresse Model-19 extruder nozzle appearing in Figure 4.In addition, we include a 1-inch section of free extrudate beyond the nozzle exit in the simulation to capture die swell.Due to the axisymmetry of the nozzle geometry and assumed flow, we are able to simplify the flow domain as a 2D axisymmetric model, which saves significant computational expense.Furthermore, our axisymmetric assumption ignores any swirling motion in the flow that may result from the extruder screw.Since creeping flow is assumed, the inertia contribution "ρ " in Equation ( 9) is ignored in our simulations.In addition, since only a short section of the free extrudate material is considered, gravitational effects are ignored, which results in the body force in Equation (9) being zero as well.In the above, F s is the tangential force, F n is the normal force, v s is the tangential velocity, v n is the normal velocity, v is the velocity vector at the free surface, and n is a unit vector normal to the free surface [24].The die swell of the free surface is predicted using the methods of spines in ANSYS-Polyflow, which is an efficient remeshing rule often applied to 2D free surface problem [31].The finite element domain is discretized into 704 nodes and 630 elements using 4-node quadrilateral elements as shown in Figure 5.The mesh size is reduced near the flow boundary as well as the nozzle exit to avoid potential singularity issues.Additionally, results obtained using a coarse mesh (448 nodes, 378 elements) and a fine mesh (960 nodes, 882 elements) were compared with those obtained with the model in Figure 5. Elastic moduli predictions appearing in the results section below were obtained using the model in Figure 5 and are within 1% absolute relative difference to the fine mesh model output.We therefore use the model in Figure 5 in the remainder of the paper to avoid the extra computational expense that would be required for a model having a finer mesh.
The boundary conditions of the flow domain appear in Figure 5, in which

•
: Flow domain inlet, where the prescribed volumetric flow rate Q is specified.In addition, a fully developed velocity profile is computed and imposed at the inlet by ANSYS-Polyflow based on Q and the selected rheology model.In the above, F is the tangential force, F is the normal force, v is the tangential velocity, v is the normal velocity, is the velocity vector at the free surface, and is a unit vector normal to the free surface [24].The die swell of the free surface is predicted using the methods of spines in ANSYS-Polyflow, which is an efficient remeshing rule often applied to 2D free surface problem [31].The finite element domain is discretized into 704 nodes and 630 elements using 4-node quadrilateral elements as shown in Figure 5.The mesh size is reduced near the flow boundary as well as the nozzle exit to avoid potential singularity issues.Additionally, results obtained using a coarse mesh (448 nodes, 378 elements) and a fine mesh (960 nodes, 882 elements) were compared with those obtained with the model in Figure 5. Elastic moduli predictions appearing in the results section below were obtained using the model in Figure 5 and are within 1% absolute relative difference to the fine mesh model output.We therefore use the model in Figure 5 in the remainder of the paper to avoid the extra computational expense that would be required for a model having a finer mesh.

Fiber Orientation Distribution Prediction
The direction of a single rigid fiber within a polymer matrix is commonly described by a unit vector (φ, ϕ), as shown in Figure 6, with coordinates [15] (φ, ϕ) =

Fiber Orientation Distribution Prediction
The direction of a single rigid fiber within a polymer matrix is commonly described by a unit vector p(ϕ, φ), as shown in Figure 6, with coordinates [15] The direction of a single rigid fiber within a polymer matrix is commonly described by a unit vector (φ, ϕ), as shown in Figure 6, with coordinates [15] (φ, ϕ) = sinφcosϕ sinφsinϕ cosφ .Advani and Tucker [15] considered the statistical behavior of the fibers using the computationally efficient orientation tensor approach where the orientation tensor evolution equation is written as which assumes isotropic rotary diffusion first proposed by Folgar and Tucker [12].Here, A and A are the second and fourth order orientation tensors, respectively, written as and where δ(ϕ, φ) is a probability distribution function and S is unit sphere surface.Note that, due to the normalization condition, the integral of δ(ϕ, φ) over the surface S equates to unity, making the trace of A equal to 1 (see e.g., [25,26]).It can also be shown that A is symmetric, yielding just five independent components in Equation (15).In additional, the vorticity tensor W appearing in Equation ( 15) is given as and the rate of deformation tensor D is evaluated by Equation (6).We evaluate the tensors W and D, from the velocity vector v computed along streamlines within the polymer melt flow field obtained from our ANSYS-Polyflow simulation result.The constant β in Equation ( 15) depends on the fiber aspect ratio as where α is the fiber aspect ratio.The interaction coefficient C I is used to capture the effect of fiber-fiber interaction, and .
γ represents the scalar magnitude of the rate of deformation tensor D. The last term in Equation ( 15) written as 2 C I .γ(I − 3A) results from the the Folgar-Tucker Isotropic Rotary Diffusion (IRD) model [12].Fu et al. [32] experimentally observed that molten short fiber polymer composite exhibited an asymmetric profile of fiber length distribution with a peak skewed toward small values of fiber length.In addition, we experimentally measured the fiber aspect ratio of the sample prepared by performing a burn-off test on the 13 wt % carbon fiber filled ABS (manufactured by Polyone, Avon Lake, OH, USA) and found that the values of the fiber aspect ratio are in a range of 10 to 60.Following Fu et al., we define β = 0.9802 (corresponding to α = 10 in Equation ( 19)).Bay and Tucker [9] defined an empirical formula for evaluating C I , which depends on the values of the fiber volume fraction and aspect ratio as where v f is the fiber volume fraction.Equation ( 20) has been effective in concentrated suspensions flows, which are typical for short fiber polymer composites in large-scale deposition processes.We also assume a fiber volume fraction of 13% following that used by ORNL in prior study during the development of the Big Area Additive Manufacturing (BAAM) [3].Therefore, we use C I = 0.0073, which is computed with Equation ( 20) for v f = 0.13 and α = 10.The fourth order fiber orientation tensor, A, is typically computed with a closure approximation in fiber orientation simulations.Prior studies have focused on the natural-type closure [33,34] and the orthotropic-type closure [35,36].In this paper, we employ the Orthotropic Closure (ORT) to compute for A from A as defined in [5].
Note that the second order orientation tensor A has seen widespread use for statistically describing the orientation of suspended fibers in narrow gap shear dominant applications [8][9][10]13].In this study, flow within the nozzle has a significant shear component; however, there is considerable extensional flows within the converging section of the nozzle, and also just outside the nozzle exit where die swell begins to form.Here, we have chosen to demonstrate the use of this common orientation tensor formulation in large-scale AM flows realizing that the limitations of the model in these flows is still to be determined.It is important to note that the orientation tensor approach does not track each individual fiber, but instead provides an indication of the degree of alignment through the nine components of tensor A. The second order orientation tensor A is the second moment of the orientation distribution function δ(ϕ, φ).Specifically, Figure 7 gives two important examples of A, in which a diagonal component of A having a value of one represents full alignment in the corresponding direction, and three diagonal components all equal to 1/3 represents the case of a uniformly random orientation.
In addition, the assumption of the initial fiber orientation state at the nozzle inlet directly influences the fiber orientation throughout the flow domain.We assume that the fiber orientation state prior to entering the nozzle has attained a fully developed steady state as the flow reaches the nozzle inlet.It is important to note that the initial condition of the fiber orientation has been found to have an influence on predicted fiber orientation in injection molding processes by Baird et al. [37].We note that the complexity in the melt flow during the extrusion process before the flow reaches the nozzle will indeed influence the fiber orientation state as the melt enters the nozzle.To better understand the effect of inlet fiber orientation on computed outputs, we performed other simulations using a uniformly random fiber orientation at the nozzle inlet.In this case, we found that using the alternate inlet fiber orientation condition had little effect on the trends in predicted extrudate fiber orientation and mechanical properties shown below.In addition, the steady state fiber orientation tensor is obtained by setting the left-hand side of Equation ( 15) to zero and then solving for components of A.
common orientation tensor formulation in large-scale AM flows realizing that the limitations of the model in these flows is still to be determined.It is important to note that the orientation tensor approach does not track each individual fiber, but instead provides an indication of the degree of alignment through the nine components of tensor .The second order orientation tensor is the second moment of the orientation distribution function δ(φ, ϕ).Specifically, Figure 7 gives two important examples of , in which a diagonal component of having a value of one represents full alignment in the corresponding direction, and three diagonal components all equal to 1/3 represents the case of a uniformly random orientation.
In addition, the assumption of the initial fiber orientation state at the nozzle inlet directly influences the fiber orientation throughout the flow domain.We assume that the fiber orientation state prior to entering the nozzle has attained a fully developed steady state as the flow reaches the nozzle inlet.It is important to note that the initial condition of the fiber orientation has been found to have an influence on predicted fiber orientation in injection molding processes by Baird et al. [37].We note that the complexity in the melt flow during the extrusion process before the flow reaches the nozzle will indeed influence the fiber orientation state as the melt enters the nozzle.To better understand the effect of inlet fiber orientation on computed outputs, we performed other simulations using a uniformly random fiber orientation at the nozzle inlet.In this case, we found that using the alternate inlet fiber orientation condition had little effect on the trends in predicted extrudate fiber orientation and mechanical properties shown below.In addition, the steady state fiber orientation tensor is obtained by setting the left-hand side of Equation ( 15) to zero and then solving for components of .

Printed Bead Elastic Properties
Advani and Tucker [15] first computed the volume-averaged stiffness tensor (C ) using the fiber orientation tensors A (A ) and (A ), where components of (C ) are given as

Printed Bead Elastic Properties
Advani and Tucker [15] first computed the volume-averaged stiffness tensor C (C ijkl ) using the fiber orientation tensors A (A ij ) and A (A ijkl ), where components of C (C ijkl ) are given as where the constants M i is computed from the five independent components of the related unidirectional composite stiffness tensor C ijkl and δ il is the Kronecker delta [38].We follow the work of Jack et al. [25] who used the Tandon-Wang analytical equations [39] to evaluate the unidirectional elastic moduli of the discontinuous fiber-reinforced polymer.The elastic properties of the fiber and the matrix used in our calculations below are given in Table 2, where we assume that the matrix and fiber are both isotropic materials.In addition, as in Section 3.2, we use a fiber aspect ratio of 10 and the fiber volume fraction is assumed to be 13% when calculating the unidirectional elastic moduli with the Tandon-Wang approach (see e.g., [26]).

Results and Discussions
The objective of this paper is to demonstrate how different representations of the melt rheology properties affects the predicted fiber orientation in the polymer melt and also the elastic properties of the resulting extrudate material.In this section, we first consider the melt flow die swell predictions computed with ANSYS-Polyflow.Then, fiber orientation tensors are obtained along streamlines within the polymer melt flow from the computed velocity field.Furthermore, the elastic properties within the solidified extrudate are evaluated from the fiber orientation tensors.We also evaluate the average elastic properties in the extrudate by numerical integration of the data points over the cross-section.

Predicted Die Swell
Ajinjeru and Duty showed that the typical wall shear rate appearing in Big Area Additive Manufacturing (BAAM) systems is between 30 and 40 s −1 , and reaches a peak value near 100 s −1 at the nozzle exit [29].In our simulations, the average wall shear rate using the PTT model with material constants from Table 1 and an inlet flow rate of Q = 100 mm 3 /s is calculated to be 36 s −1 , with a peak value of 87 s −1 at nozzle exit, which agrees well with the literature data [29].Here, the average wall shear rate .γ w is computed from the wall shear rate .γ w on Γ 4 (cf.Figure 5) as .
where L is the length of the nozzle exit tube defined by Γ 4 in the x 2 direction.
Die swell profiles for the nozzle flow problem defined in Section 2 using each of the rheology models defined above appear in Figure 8.The die swell just downstream of the nozzle exit is assessed using the apparent swell ratio B defined as where d is the steady state swell flow diameter evaluated along the free surface downstream of the die exit (length of Γ 6 appearing in Figure 5), and d 0 is the nozzle exit diameter.The computed data for B at the Γ 6 surface is given in Table 3.The apparent die swell ratio B = 1.133 computed using the Newtonian fluid model agrees with the swell ratio of 1.13 proposed by Reddy and Tanner [40].
The steady state die swell ratios calculated using the Power law and Carreau-Yasuda rheology model are nearly identical and significantly lower than the swell ratio of B = 1.199 computed using the PTT model.Furthermore, the die swell profile obtained using the SV model converges to that computed with the PTT model as the profiles reach steady state.Note that simulation time using the SV model was 75 s which is much less than the 328 s when using the PTT model with the mesh given in Figure 5. Therefore, for a larger size of flow domain or a finer mesh quality, the SV model is a good candidate to qualitatively solve the flow problem with less computational cost.

Computed Fiber Orientation Distribution
In this work, the fiber orientation is computed as described in Section 3.2 above.Our primary interest is in the direction of extrusion, i.e., in the direction of the positive x axis in Figure 5. Therefore, our primary focus is on fiber alignment in x , which is best represented by the A component of the second order orientation tensor .
The solution of A along various streamlines shown in Figure 9 is computed based on the flow kinematics solved using the PTT rheology model with Q = 100 mm 3 /s.Values of A computed for each of the rheology models showed similar trends as that appearing in Figure 9, so that these other results are omitted here for conciseness.ANSYS-CFD-Post [41] generates the 2D surface streamlines

Computed Fiber Orientation Distribution
In this work, the fiber orientation is computed as described in Section 3.2 above.Our primary interest is in the direction of extrusion, i.e., in the direction of the positive x 2 axis in Figure 5. Therefore, our primary focus is on fiber alignment in x 2 , which is best represented by the A 22 component of the second order orientation tensor A.
The solution of A 22 along various streamlines shown in Figure 9 is computed based on the flow kinematics solved using the PTT rheology model with Q = 100 mm 3 /s.Values of A 22 computed for each of the rheology models showed similar trends as that appearing in Figure 9, so that these other results are omitted here for conciseness.ANSYS-CFD-Post [41] generates the 2D surface streamlines based on the mesh quality.Due to the mesh defined in Figure 5, we consider flow velocity fields along eight streamlines computed from the finite element results achieved by ANSYS-Polyflow.Note that A 22 values near unity indicate fibers are highly aligned along the x 2 direction.It can be seen that the fiber orientation starts from steady state at flow inlet (appearing as FI in Figure 9) as assumed.The A 22 components then starts to separate as the flow approaches the nozzle convergent zone (appearing as CZS in Figure 9).Orientation states along all streamlines increase before the flow reaches the exit of the convergent zone (appearing as CZE in Figure 9).The peak value of the A 22 component occurs at the convergent zone exit.Then, the orientation tensors at inner region streamlines decrease while those located at outer region increase as the flow propagating to the nozzle exit (appearing as NE in Figure 9).Once the polymer melt passes the nozzle exit at NE, values of A 22 in the outer region increase immediately and those more central begin to increase.This change occurs due to the shear rate limitation vanishing at the outer boundary just after nozzle exit.The velocity along the outer boundary accelerates first, causing fibers nearby to orientate in the flow direction.In addition, the elongational flow near the center of the nozzle accelerates so that the extrudate attains a uniform speed at some point not far from the nozzle exit.The final state of fiber orientation is set once variation across the bead ceases and a plug flow develops.domain (i.e., across , in Figure 5) for all simulations considered here appear in Figure 10.Results indicate that fibers are highly aligned near the edge of the flow where shear rates are high.Fiber alignment then decreases just inside this outer band forming an intermediate slightly misaligned region.An increase in alignment then occurs within the core region towards the center of the extrudate.In addition, it can be seen that the PTT model yields the lowest alignment in x among the applied rheology models.Alternatively, the Power law and Carreau-Yasuda law result in a similar steady state fiber orientation, which are the highest among these results.Moreover, the orientation result obtained using the SV model shows a good agreement with that of the PTT model in the shear dominant region but varies at other locations within the flow.Finally, the Newtonian model yields an intermediate value of the fiber orientation, somewhat positioned between the GNF laws and the viscoelastic model results.Similar to the results appearing in Figure 9, we also compute the fiber orientation tensor along streamlines in the flow domain using other rheology models.Values of A 22 at the end of the flow domain (i.e., across Γ 6 , in Figure 5) for all simulations considered here appear in Figure 10.Results indicate that fibers are highly aligned near the edge of the flow where shear rates are high.
Fiber alignment then decreases just inside this outer band forming an intermediate slightly misaligned region.An increase in alignment then occurs within the core region towards the center of the extrudate.In addition, it can be seen that the PTT model yields the lowest alignment in x 2 among the applied rheology models.Alternatively, the Power law and Carreau-Yasuda law result in a similar steady state fiber orientation, which are the highest among these results.Moreover, the orientation result obtained using the SV model shows a good agreement with that of the PTT model in the shear dominant region but varies at other locations within the flow.Finally, the Newtonian model yields an intermediate value of the fiber orientation, somewhat positioned between the GNF laws and the viscoelastic model results.

Elastic Properties across the Exrudate
Finally, we calculate the elastic properties for the fiber-reinforced polymer using the volume-averaged stiffness tensor computed from the steady state fiber orientation tensor (cf.Jack et al. [25]) and the Tandon-Wang analytical equation (cf.Tandon et al. [39]).Computed elastic properties across the extrudate at boundary Γ 6 (cf. Figure 5) for each rheology model considered here appear in Figure 11.For conciseness, we omit the results obtained using the Carreau-Yasuda rheology model since these results are very similar to that obtained in our Power law model simulations.
Our results show that the elastic properties of the composite extrudate are enhanced by the fiber reinforcement in comparison to the properties of neat ABS (dash lines appearing in Figure 11), particularly the E 22 component, which is elastic modulus along the extrusion direction.In addition, the elastic modulus along the extrusion direction (E 22 ) is not quite uniform while the moduli at other directions (E 11 , E 33 , G 12 , G 23 , G 13 ) are in small variance across the extrudate.In comparison, the Power law yielded relatively higher estimation in the principal modulus (E 22 ) and the PTT model results in the lowest prediction.In addition, the Power law model, the PTT model as well as the SV model show a more sharp variation in the E 22 component than that yielded by the Newtonian model, which shows trends that are similar to results calculated for a small FDM nozzle in Heller et al. [14].
the elastic modulus along the extrusion direction (E ) is not quite uniform while the moduli at other directions (E , E , G , G , G ) are in small variance across the extrudate.In comparison, the Power law yielded relatively higher estimation in the principal modulus (E ) and the PTT model results in the lowest prediction.In addition, the Power law model, the PTT model as well as the SV model show a more sharp variation in the E component than that yielded by the Newtonian model, which shows trends that are similar to results calculated for a small FDM nozzle in Heller et al. [14].Furthermore, we average the elastic properties of the extrudate by numerically integrating the data appearing in Figure 11 as written in Equation ( 23) by the trapezoidal rule [42] and the results are given in Table 4: Furthermore, we average the elastic properties of the extrudate by numerically integrating the data appearing in Figure 11 as written in Equation ( 24) by the trapezoidal rule [42] and the results are given in Table 4: where X refers to any of the data points appearing in Figure 11, X avg is the area-averaged data, and r o is the outer radius of the extrudate.Calculated values of X avg corresponding to the elastic properties in Figure 11 appear in Table 4.Note that we assume that the centre of the extudate exhibits the same orientation behavior as the data in the streamline nearest to the centre-line, so as to achieve the fiber orientation over the entire cross-section of the extrudate (cf. Figure 11).In addition, due to the flow domain being built as a 2D axis-symmetric model, the integration is carried out in polar axis in order to include the effects over the entire cross-section area of the extrudate.The data computed in Table 4 indicate that the printed composite exhibit quasi transverse isotropic mechanical behavior.In detail, the Carreau-Yasuda law (appearing as Carreau-Y.in Table 4) yields the highest averaged E 22 value, while the PTT model results in the lowest.In contrast, the shear moduli yielded by the PTT model is higher than the GNF models.Finally, as we consider more non-viscous effects of the flow (see data from the Power law (as well as the Carreau-Yasuda law), the Newtonian model, the SV model to the PTT model), it can be seen that the principal modulus (E 22 ) decreases while other moduli increase.In addition, the transverse isotropic behavior of results predicted by the PTT model is the most obvious one, in which the E 22 .and G 13 are unique values and E 11 , E 33 as well as G 12 , G 23 show high agreement.Our material stiffness predictions agree remarkably well with test data appearing in the literature.Duty et al. [4] measured the Young's modulus of a 13 per cent carbon fiber reinforced ABS printed bead, reporting a mean value of 7.24 GPa and standard deviation of 0.59 GPa.Our predictions of the elastic modulus for the same material system shown as E 22 in Table 4 are in good agreement with these previously published experimental results.Note that all predicted values of E 22 in Table 4 are within one standard deviation of Duty's mean experimental value, regardless of rheology model used, except for the value of 6.57GPa obtained using the PTT model, which is at 1.1 standard deviation from the experimental mean.In addition, Duty found the stiffness of a BAAM printed tensile test sample to be highly anisotropic, which is also seen in our results.Furthermore, the previously reported test sample transvers moduli show some anisotropy in the transverse plane (i.e., Ey = 2.26 GPa and Ez = 2.56 GPa).Our results obtained using GNF fluids (i.e., power law and Carreau-Yasuda fluids) show similar trends.While we are not able to show experimental data specific to each rheology model, the overall favorable comparison with previously published experimental work supports our computational approach.

Conclusions
Polymer melt flow through a large-scale polymer deposition extrusion nozzle was simulated with the finite element method using ANSYS-Polyflow.These simulations included flow within the nozzle in addition to the free surface die swell flow just outside the nozzle exit.Several rheology models were compared in this work including a Newtonian model, a Power law model, a Carreau-Yasuda model, as well as a multi-mode PTT model and a Simplified Viscoelastic (SV) material model.Rheology data obtained experimentally using the HAAKE MARS 40 rheometer.
It was found that characterizing the melt flow by different rheology models yielded noticeable variation in predicted die swell, fiber orientation distribution and the ultimate elastic behavior of the extruded composites.The predicted die swell yielded by the PTT model was higher than those resulted by the Generalized Newtonian Fluid (GNF) models including the Newtonian model.The SV model yielded die swell results that agreed well with those from the PTT model by careful adjustments of the rheology model parameters.
Through the weakly couple formulation, the fiber orientation distribution within the extrudate was calculated from the melt flow velocity field.High fiber alignment in the direction of extrusion occurred near the high-shear flow edge region of the extrudate as well as the near-center region, which was due to the elongational effects of the free flow.Among the applied rheology models, the PTT model yielded the lowest principal fiber alignment while the Power law model resulted in the highest fiber orientation in polymer extrusion direction.
The elastic properties of a printed extrudate were evaluated based on the predicted fiber orientation distributions, in which the estimated elastic modulus along extrusion direction showed noticeable variance across the extrudate.The numerically-integrated averaged elastic moduli showed a good agreement with the published experimental work.The estimation indicated the composite extrudate exhibited quasi transverse isotropic behavior.In detail, the GNF models yielded higher Young's modulus along the principal direction while the PTT model resulted in a lower principal Young's modulus but higher values of shear moduli.This indicates that, by considering the non-viscous rheology effects, the elastic properties of extrudate through Additive Manufacturing (AM) systems reduced at a longitudinal direction but increased at shear directions.
In addition, the SV model yielded relatively similar data of fiber orientation distribution as well as elastic properties in comparison with the PTT model, especially in the shear dominant flow boundary, yet cost less computational time than the PTT model.In the future study of 3D deposition modelling of large-scale AM, the computationally cost-effective SV model is a reasonable alternative for conventional viscoelastic fluid models (e.g., PTT model).

Figure 1 .
Figure 1.Flow-induced fiber orientation occurs during the process of large-scale additive manufacturing.

Figure 1 .
Figure 1.Flow-induced fiber orientation occurs during the process of large-scale additive manufacturing.

J
. Compos.Sci.2018, 2, 10 5 of 17 die swell profile that is in good agreement with results obtained using the PTT viscoelastic fluid model.In this study, θ and Ψ are thus defined as 0.26 and 0.47, respectively.

Figure 3 .
Figure 3. Curve fitting the experimental rheology data using the Phan-Thien-Tanner model.

Figure 3 .
Figure 3. Curve fitting the experimental rheology data using the Phan-Thien-Tanner model.

Figure 3 .
Figure 3. Curve fitting the experimental rheology data using the Phan-Thien-Tanner model.

Figure 5 .
Figure 5. Mesh and boundary condition of the flow domain.

Figure 5 .
Figure 5. Mesh and boundary condition of the flow domain.

Figure 7 .
Figure 7. Fiber orientation tensors of (a) fully alignment along x direction and (b) uniformly random alignment.

Figure 7 .
Figure 7. Fiber orientation tensors of (a) fully alignment along x 2 direction and (b) uniformly random alignment.

Figure 8 .
Figure 8. Die swell profiles predicted by using different rheology models.

Figure 8 .
Figure 8. Die swell profiles predicted by using different rheology models.

Figure 9 .
Figure 9.A component of fiber orientation solution computed using the flow kinematics solved by the PTT model at Q equates 100 mm 3 /s.Figure 9.A 22 component of fiber orientation solution computed using the flow kinematics solved by the PTT model at Q equates 100 mm 3 /s.

Figure 9 .
Figure 9.A component of fiber orientation solution computed using the flow kinematics solved by the PTT model at Q equates 100 mm 3 /s.Figure 9.A 22 component of fiber orientation solution computed using the flow kinematics solved by the PTT model at Q equates 100 mm 3 /s.

Figure 9 .
Figure 9.A component of fiber orientation solution computed using the flow kinematics solved by the PTT model at Q equates 100 mm 3 /s.

Figure 10 .
Figure 10.A component at the flow domain exit solved using employed rheology models.

Figure 10 .
Figure 10.A 22 component at the flow domain exit solved using employed rheology models.

Figure 11 .
Figure 11.Elastic properties across the printed extrudate computed with flow kinematics obtained by (a) the Newtonian model; (b) the Power law model; (c) the PTT model; (d) the Simplified Viscoelastic (SV) model.Note that the legend shown in (b) works for all four sub-figures.

Figure 11 .
Figure 11.Elastic properties across the printed extrudate computed with flow kinematics obtained by (a) the Newtonian model; (b) the Power law model; (c) the PTT model; (d) the Simplified Viscoelastic (SV) model.Note that the legend shown in (b) works for all four sub-figures.

Table 2 .
Elastic properties of the fiber and the matrix of the composite.

Table 3 .
Apparent swell ratio values resulted by applied rheology models.

Table 3 .
Apparent swell ratio values resulted by applied rheology models.

Table 4 .
Averaged elastic properties of the printed extrudate.