A Fully Coupled Simulation of Planar Deposition Flow and Fiber Orientation in Polymer Composites Additive Manufacturing

Numerical studies for polymer composites deposition additive manufacturing have provided significant insight promoting the rapid development of the technology. However, little of existing literature addresses the complex yet important polymer composite melt flow–fiber orientation coupling during deposition. This paper explores the effect of flow–fiber interaction for polymer deposition of 13 wt.% Carbon Fiber filled Acrylonitrile Butadiene Styrene (CF/ABS) composites through a finite-element-based numerical approach. The molten composite flow in the extrusion die plus a strand of the deposited bead contacting the deposition substrate is modelled using a 2D isothermal and incompressible Newtonian planar flow model, where the material deposition rate is ~110 mm/s simulating a large scale additive manufacturing process. The Folgar–Tucker model associated with the Advani–Tucker orientation tensor approach is adopted for the evaluation of the fiber orientation state, where the orthotropic fitted closure is applied. By comparing the computed results between the uncoupled and fully coupled solutions, it is found that the flow-orientation effects are mostly seen in the nozzle convergence zone and the extrusion-deposition transition zone of the flow domain. Further, the fully coupled fiber orientation solution is highly sensitive to the choice of the fiber–fiber interaction coefficient CI, e.g., assigning CI as 0.01 and 0.001 results in a 23% partial relative difference in the predicted elastic modulus along deposition direction. In addition, Structural properties of deposited CF/ABS beads based on our predicted fiber orientation results show favorable agreements with related experimental studies.


Introduction
Due to superior structural properties and thermal dimensional stability, short fiber reinforced polymer composites have experienced widespread applications in Polymer Deposition Additive Manufacturing (PDAM), otherwise widely known as Fused Filament Fabrication (FFF) for desktop 3D printing. PDAMs open up new possibilities for rapid prototyping of spare parts for engineering applications [1] and continues to see increasing implementations in building large dimension parts and tool with advanced large area screw-extruder-based PDAM systems [2]. While significant progress has been made in the development of the PDAM technology, only recently has simulation research helped to improve the understanding of the PDAM process. Finite-volume-based numerical approaches have been used to simulate the melt flow dynamics in the FFF process of ABS polymers [3], where the residual stresses within the printed part [4], the shape of deposited beads [5], processing parameters (e.g., extrusion temperature and spacing of the filaments) [5], and viscoelastic stress fields of the extrudate [6] were investigated. This research showed that viscoelastic effects in the polymer melt flow played an important role in defining the shape of the deposited bead and inter-layer bonding. Urbanic and DiCecco [7] investigated the stair-case effect resulting from the layered deposition process on the surface roughness of an FFF fabricated part and suggested that adaptive slicing and surface offsets should be included in surface roughness prediction tools. Daniyan et al. [8] studied the thermal stress and displacement fields of the fused deposition process using Abaqus (Dassault Systems, Waltham, MA, USA) and showed that satisfactory material strength under different loading conditions could be obtained with the help of their numerical analyses. Nasirov and Fidan [9] investigated the effects of infill-voids on the anisotropic structural properties using the micromechanical representative volume element approach and their predicted tensile moduli for a line infill pattern showed good agreement with experimental data. Osswald et al. [10] numerically evaluated the melt extrusion process of ABS materials in FFF and showed that the maximum melting rate is controlled by filament insertion force. Nienhaus [11] studies effect of the nozzle geometry in FFF processing of ABS and suggested that the filament feed velocity, extrusion temperature, and filament insertion force were closely inter-related. Zhou et al. [12] optimized FFF processing parameters to reduce the warpage of printed ABS parts through a voxel-based finite element simulation and showed that the print chamber temperature has a significant effect on thermal deformation. D'Amico and Peterson presented an adaptive finite element method to evaluate heat transfer in polymer deposition additive manufacturing [13,14]. They found that large volume deposition processes retain more heat than desktop size printers which improved the interlayer diffusion, but promoted slumping or sagging. Mcllroy and Olmsted numerically evaluated the polymer molecules disentanglement by employing a modified RoliePoly model and found that the polymer melt significantly disentangled when the deposition process involved corner flow geometries [15]. Comminal and co-workers developed a CFD-based algorithm that simulated the shape of deposited bead [16,17], voids between the interlayer beads [18], and material deposition in corners [19] for the FFF polymer deposition process. Balani et al. [20] studied the melt extrusion stability of PLA polymer in an FFF nozzle flow and provided a theoretical model that facilitates selections of printing parameters for manufactured parts with improved qualities. McIlroy and Graham [21] stated that the flow-enhanced crystallization of the semi-crystalline polymer (i.e., PCL) significantly affected the interlayer bead bonding strength in FFF processing. Phan et al. [22] investigated the flow dynamics in the nozzle region with the Phan-Thien-Tanner viscoelastic model through the finite element suite ANSYS-Polyflow (Ansys, Inc., Canonsburg, PA, USA) and showed that increased extrusion rate of FFF was limited by nozzle pressure. Yang and Zhang [23] also investigated the thermal history an FFF printed part with different infill patterns through the finite element approach, where the honeycomb infill was shown to result a minimum temperature gradient. Shadvar et al. [24] studied the die swell of ABS polymer during the FFF extrusion process through the finite element suite ANSYS-Polyflow and showed that die swell played a key role in improving the dimensional accuracy of FFF parts. Additional earlier efforts beyond the recent three years for the modelling of polymer deposition processes can be found in related review articles such as [25,26]. From the above review, it can be seen that numerical modelling of PDAM and its related processes has seen significant study, but the results have not provided adequate insight to sufficiently improve the qualities of PDAM-fabricated parts. Unfortunately, a majority of the prior literature focuses on the virgin polymer materials, which may not be directly applied to the PDAM processing of fiber reinforced composites. Additional research is needed to better understand PDAM processing of polymer composites which have seen increasingly application in aerospace, automotive and marine applications, especially for the recent-emerging large scale PDAM systems (e.g., see Figure 1).  PDAM processing of short-fiber polymer composites typically yields an anisotropic material microstructure which affects the bead's thermal expansion coefficient [2], elastic moduli [27], and thermal conductivity [28]. Bead elastic modulus, for example, is often much higher in the direction of printing as compared to transverse directions. Brenken et al. [29] studied the thermal history of material deposition additive manufacturing for the 50 wt.% CF/PPS, where the anisotropic thermal conductivity of the short-fiber polymer composite was modelled. Compton et al. [30] studied the macro thermal history of a large area additively manufactured carbon fiber filled ABS part through the finite element method assuming an anisotropic thermal conductivity. Hoskins et al. [31] simulated the residual thermal stress in a printed cuboidal part of CF/ABS, where the coefficient of thermal expansion of the deposited beads was defined through a non-homogenized approach based on locally measured fiber orientation states. To explain the factors that produce anisotropic material properties, in-depth analyses for the polymer composite deposition flow are needed. Nixon et al. [32] are likely the first to investigate the fiber orientation state within a FFF extrusion nozzle, where Moldflow (Autodesk, Inc., San Rafael, CA, USA) was applied to examine the effects of three different nozzle geometries on the resulting fiber orientation at the nozzle exit. A team lead by Smith presented numerical approaches for evaluating the nozzle extrusion flow [33][34][35][36] and polymer deposition flow [37,38] in the context of PDAM for short-fiber polymer composite materials. Their work has provided valuable insight into the extrudate swell [34,35,38] as well as the processproperty map. Elastic constants [33][34][35][36][37][38] and thermal expansion coefficients [37,38] of a deposited CF/ABS were evaluated and favorable agreement with reported data from related experimental work [26] was reported. Their work largely employed a weakly coupled formulation to simulate fiber reinforced polymer flow, where the presence of fibers is neglected in the computation of the melt flow kinematics. Recently, efforts have been made to evaluate the effect of the fiber orientation in the polymer melt flow within the PDAM process. Mezi et al. [39] studied the die swell of a fully coupled Newtonian fiber suspension flow for the extrusion FFF process. Yang et al. [40], Bertevas et al. [41], and Ouyang et al. [42,43] employed the Smoothed Particle Hydrodynamics (SPH) approach that effectively simulates the flow-fiber orientation coupling behavior in the bead deposition process of fiber reinforced composites. Additionally, Wang and Smith recently developed a finite-element-based algorithm that captures the mutually dependent effect between the polymer flow rheology and fiber reinforcement orientation in the PDAM nozzle-extrudate flow [44].
Nevertheless, Wang and Smith [44] only depict the fully coupled flow/orientation interactions axisymmetric material extrusion, where the important deposition feature is not included due in part to the high degree of divergence issue encountered in their simulation. This paper extends the finite-element-based fully coupled flow/orientation method to the polymer extrusion deposition process of concentrated fiber-reinforced polymer composites melt flow, which provides more insight into the material properties of PDAM processing of short-fiber polymer composites typically yields an anisotropic material microstructure which affects the bead's thermal expansion coefficient [2], elastic moduli [27], and thermal conductivity [28]. Bead elastic modulus, for example, is often much higher in the direction of printing as compared to transverse directions. Brenken et al. [29] studied the thermal history of material deposition additive manufacturing for the 50 wt.% CF/PPS, where the anisotropic thermal conductivity of the short-fiber polymer composite was modelled. Compton et al. [30] studied the macro thermal history of a large area additively manufactured carbon fiber filled ABS part through the finite element method assuming an anisotropic thermal conductivity. Hoskins et al. [31] simulated the residual thermal stress in a printed cuboidal part of CF/ABS, where the coefficient of thermal expansion of the deposited beads was defined through a non-homogenized approach based on locally measured fiber orientation states. To explain the factors that produce anisotropic material properties, in-depth analyses for the polymer composite deposition flow are needed. Nixon et al. [32] are likely the first to investigate the fiber orientation state within a FFF extrusion nozzle, where Moldflow (Autodesk, Inc., San Rafael, CA, USA) was applied to examine the effects of three different nozzle geometries on the resulting fiber orientation at the nozzle exit. A team lead by Smith presented numerical approaches for evaluating the nozzle extrusion flow [33][34][35][36] and polymer deposition flow [37,38] in the context of PDAM for short-fiber polymer composite materials. Their work has provided valuable insight into the extrudate swell [34,35,38] as well as the process-property map. Elastic constants [33][34][35][36][37][38] and thermal expansion coefficients [37,38] of a deposited CF/ABS were evaluated and favorable agreement with reported data from related experimental work [26] was reported. Their work largely employed a weakly coupled formulation to simulate fiber reinforced polymer flow, where the presence of fibers is neglected in the computation of the melt flow kinematics. Recently, efforts have been made to evaluate the effect of the fiber orientation in the polymer melt flow within the PDAM process. Mezi et al. [39] studied the die swell of a fully coupled Newtonian fiber suspension flow for the extrusion FFF process. Yang et al. [40], Bertevas et al. [41], and Ouyang et al. [42,43] employed the Smoothed Particle Hydrodynamics (SPH) approach that effectively simulates the flow-fiber orientation coupling behavior in the bead deposition process of fiber reinforced composites. Additionally, Wang and Smith recently developed a finite-element-based algorithm that captures the mutually dependent effect between the polymer flow rheology and fiber reinforcement orientation in the PDAM nozzle-extrudate flow [44].
Nevertheless, Wang and Smith [44] only depict the fully coupled flow/orientation interactions axisymmetric material extrusion, where the important deposition feature is not included due in part to the high degree of divergence issue encountered in their simulation. This paper extends the finite-element-based fully coupled flow/orientation method to the polymer extrusion deposition process of concentrated fiber-reinforced polymer composites melt flow, which provides more insight into the material properties of beads formed with large area additive manufacturing polymer deposition. A two-dimensional (2D) planar Newtonian creeping flow is assumed to analyze extrudate swell of the freely deposited bead where the flow domain is computed using a one-dimensional remeshing technique [45]. Fiber orientation is modeled with Advani-Tucker fiber orientation tensors [46] using the orthotropic fitted closure approximation [47] and isotropic rotary diffusion [48]. The nonlinear fiber orientation equations are solved through the Newton-Raphson iterative method as in [44]. Unlike prior transient simulations (e.g., [16][17][18][19][41][42][43]), the finite element nodal solutions provide quasi-steady solutions for flow velocity and fiber orientation field of the deposition process, which is then used to compute the structural properties of a solidified composite bead. Our results from the fully coupled simulation are compared to a weakly coupled result to expose the importance of incorporating the effect of fiber orientation on the polymer melt velocity solution.

Methodology
Our numerical simulation approach for fiber reinforced polymer composite flow in PDAM, including the identifications of the flow fields of the molten material, die swell, and fiber orientation state within the flow suspension, are described in this section.

Governing Equations
For polymer melt flows in die extrusion and deposition, an incompressible, isothermal, and highly viscous creeping Newtonian fluid flow is often assumed (e.g., [33,37,38]), such that the equations of mass and momentum conservation within the flow are written as and where v is the velocity vector, ρ is the density of the continuum, and f is the body force vector [49]. In the above, σ is the Cauchy stress tensor, which may be expressed as where P is the pressure, I is the identity tensor, and τ is part of the stress tensor associated with the viscosity [49]. In fiber suspensions, τ may be written as a function of the fiber orientation state within the flow, such that where η is the viscosity of the Newtonian fluid, and N p is the particle number that characterizes the intrinsic anisotropic effect of fibers on the flow rheology [50]. In Equation (4), D is the second-order rate of deformation tensor and A is the fourth-order fiber orientation tensor which is described in detail below. We note that a non-Newtonian flow model can better represent the rheology of polymer composite melts (e.g., shear thinning). However, incorporating a shear-dependent viscosity models (e.g., power law model, Carreau law model, etc.) increases the nonlinearity of the fully coupled problem considered in this study, leading to additional convergence issues beyond those discussed below. Therefore, the flow/fiber-orientation coupling is modeled as a highly viscous Newtonian flow in this work as in prior PDAM research (see e.g., [33,37,38]), leaving the use of a non-Newtonian fluid for future work. Jeffery's work [51] which describes the motion of a single rigid ellipsoidal particle in a pure shear flow forms the basis for most fiber orientation studies. More recently, Folgar and Tucker [48] extended Jeffery's theory to analyze the interaction between fibers in a non-dilute fiber suspension using a fiber orientation distribution function. Further, Advani and Tucker proposed the fiber orientation tensor approach to quantify the fiber alignment state for concentrated suspension systems [46], which requires fewer independent variables than that of the Folgar-Tucker model [48]. The Advani-Tucker orientation tensor evolution equation with Folgar-Tucker orientation diffusion is written as where the second-and fourth-order fiber orientation tensors are defined respectively as Here, p is the unit vector defining the orientation of a single rigid fiber [46]. The angle bracket "< >" in Equation (6) denotes orientation average over all directions, weighted by the probability distribution function of the orientation [50]. The fourth order orientation tensor is approximated through the widely applied orthotropic fitted closure [47]. We note that the effect of applying different closure approximations in our simulation would be an interesting future work. In addition, Favaloro and Tucker [52] summarized the fiber orientation kinetic equations including the more recently developed RSC and ARD models. Nevertheless, applying these orientation models brings in other parameters in the simulation, which lies beyond the major scope of this study and may be addressed in separate in-depth work.
The tensors W and D are the vorticity tensor and rate-of-deformation tensor (as introduced above) of the suspension flow, respectively, which can be written as where ∇v is the velocity gradient and the superscript T refers to the tensor transpose [49]. In addition, the parameter λ in Equation (5) is a function of fiber geometry. For an ellipsoidal fiber, λ can be evaluated as where a r (g r ) is the hydrodynamic aspect ratio of the ellipsoidal fiber, and g r is the geometrical aspect ratio of the fiber (i.e., the ratio of the fiber length to diameter) [53]. Note that a r = g r for ellipsoidal fibers which is not the case for fibers of other shapes (cf. Zhang et al. [53]). The last term appearing in Equation (5) is referred as the isotropic rotary diffusion that was first proposed by Folgar and Tucker [48], which includes the empirically obtained fiber interaction coefficient C I , to provides a means for incorporating the effect of fiber-fiber interaction. Bay [54] proposed C I is a function of volume fraction v f and hydrodynamic aspect ratio as The dimensionless particle number N p in Equation (4) is of great importance as it defines the amount of anisotropy of the melt viscosity. Several models have been proposed to address material anisotropy for fiber suspensions for different fiber concentrations (see e.g., [55][56][57]). Herein, we adopt the rheological constitutive equation presented by Dinh and Armstrong [56] due to its simple implementation and relatively high reliability [44]. The Dinh-Armstrong model describes the effect of fibers on viscosity in terms of fiber orientation parameters including fiber aspect ratio and fiber volume fraction. For narrowgap shear flows, Dinh and Armstrong defined N p as where v f is the fiber volume fraction, and H and S are material coefficients related to the fiber suspension system [44]. For suspensions with a fully aligned orientation state (which is often approached in PDAM nozzle flow, e.g., [32][33][34][35]), H is given by (cf. [58]) and S = 0, as suggested by the Dinh-Armstrong model assuming that the particle's thickness can be ignored [58]. The model assumes fibers behave as slender-bodies, which allows for hydrodynamic interactions between fibers, at least in an average sense. In this approach, simulations can be performed beyond the dilute regime, while Jeffery's solution for flow around a single ellipsoid are limited to the dilute case [58]. This is pivotal in the analyses the short-fiber polymer composite systems used in PDAM, where most filled polymers are highly concentrated fiber suspensions.

Finite Element Simulation
The fully coupled flow and fiber orientation problem has been solved in prior research (e.g., [41,44]). It is common to use the finite element method to compute the melt flow velocity and the fiber orientation tensor field is solved by integrating orientation evaluation equation (i.e., Equation (5)) along the flow streamlines (e.g., [33][34][35][41][42][43]). Alternatively, in this work, we use an approach as in Wang and Smith [44] (and also VerWeyst and Tucker [50]) who employed the finite element method directly to both the flow and fiber orientation governing equations. In this approach, we employ the standard Galerkin finite element approach to compute elemental nodal velocities d e , such that where the element stiffness matrix K e is and the element nodal force vector F e is where B e s = ∇ s N e , N e is the elemental interpolation function matrix, and ∇ s is symmetric gradient operator for a 2-D axisymmetric coordinate system (see. Reddy [59]). In the above, V is an integrated form of viscosity matrix that simplifies the weak form of the momentum equation to the concise format shown here (see e.g., [44]).
Similarly, the finite element form of the fiber orientation governing equation can be written as K e a d e a = F e a , where the element fiber orientation system matrix K e a is and d e a is the elemental nodal velocity and second-order orientation tensor solutions, respectively [44] [44]. Lastly, the element matrices and vectors appearing in Equation (15) are assembled to form the global finite element matrix equations in the usual manner. The fiber orientation finite element Equation (15) is a nonlinear system, which we solve using the Newtonian-Raphson (N-R) iterative method [60] where the tangent matrix of Equation (15) is computed using the forward finite difference method with a step size of 10 −5 for d e a , which yields trivial difference as compared to results computed by applying 10 −6 as the step size. In addition, we note that the steady state fiber orientation solution within the flow domain is computed with an implicit time marching method (i.e., first order backward finite difference method) as in [44,50] and index q appearing in Equation (17) refers to the N-R iteration. We note that applying finite difference approximation in computing the tangent matrix of Equation (15) reduces the convergence rate from the quadratic form of regular N-R method, while the accuracy of the solved nodal solutions is retained as indicated in [50]. A detail derivation process for the above is omitted for conciseness and can be found in the dissertation of the first author [61], including the elemental stiffness matrices and force vectors, and the tangent matrix of the finite element system, etc.
The Galerkin finite element solution of the second-order fiber orientation tensor (cf. Equations (15)-(17)) may exhibit spatial instabilities due to the lack of a diffusion term in the governing equation (cf. Equation (5)) [44]. Herein, we apply the Streamline Upwind Petrov Galerkin (SUPG) method [62] to stabilize the nodal orientation tensor solution of d e a in Equation (15). The SUPG method replaces the Galerkin weight function (i.e., w = N e a ) with where α is a scale factor controlling the magnitude of the streamline upwinding [63]. Here, we set α = 0.5 as in Wang and Smith for a similar problem [44] which provides just enough upwinding to eliminate spatial instabilities without degrading the orientation tensor solution.

Extrusion-Deposition Flow Modelling
This paper considers large area extrusion deposition additive manufacturing of short fiber polymer composites. The flow domain of interest includes the internal geometry of the extrusion die and a single layer of a planar deposited bead on a translating substrate as shown in Figure 2a. The extrusion-deposition process modeled here occurs in a short period of time justifying an isothermal assumption such that thermal effects are not included in the computation. Note that non-isothermal studies for similar flow geometries were performed using other simulation approaches (e.g., [42,43]). A 2D planar flow domain is employed where a unit width (i.e., the out-of-plane dimension) is assumed for the deposited bead (similar to that appearing in [38,[41][42][43]). The nozzle internal flow geometry is based on the Strangpresse large-scale additive manufacturing Model 19 single screw extruder nozzle (cf. Figure 2b). The height of the deposited bead is 3 mm (i.e., distance between nozzle end and deposition substrate), where the outer surface of the nozzle is assumed to yield a negligible effect on the flow surface of the planar deposited bead [38]. The ratio of the deposited bead length to its height is 10:1, which ensures a steady-state condition of the flow and fiber orientation is achieved within the modeled flow domain. It should be noted that the relative distance between the nozzle exit and the substrate is of great importance in form the bead height, however, this factor is out of the main scope of this research and are discussed in prior works (see e.g., [38,42]). between nozzle end and deposition substrate), where the outer surface of the nozzle is assumed to yield a negligible effect on the flow surface of the planar deposited bead [38]. The ratio of the deposited bead length to its height is 10:1, which ensures a steady-state condition of the flow and fiber orientation is achieved within the modeled flow domain. It should be noted that the relative distance between the nozzle exit and the substrate is of great importance in form the bead height, however, this factor is out of the main scope of this research and are discussed in prior works (see e.g., [38,42]). The planar deposition flow model appears in Figure 3, where the flow domain is separated into three subdomains including the internal nozzle flow domain Ω 1 , the deposition turning flow domain Ω 2 , and the deposited bead domain Ω 3 . The boundary conditions for the flow domain in Figure 3 are as follows:

Results and Discussions
In this study, we use the 13 wt.% CF/ABS as the polymer composite melt material model that has seen widespread applications in the large area polymer deposition additive manufacturing [2,27,36]. We assume the Newtonian viscosity of the melt is = 817 Pa•s, which is the shear viscosity of the ABS polymer at 230 ℃ and shear rate 100 , and the fluid density is = 1154 kg/m 3 , both of which are from Heller et al. [38]. In addition, we assume an isotropic fiber orientation field for the entire flow domain as a starting point for the fiber orientation calculation iterations when computing the fully coupled flow-orientation solutions, as in [44,50]. Moreover, the hyperbolic form of Equation (5) requires an initial condition of the second-order fiber orientation tensor at the flow inlet. Herein, we imposed the fiber orientation state of a fully developed flow, which assumes that the orientation state reaches steady state upstream of the flow inlet. The initial condi- In the above, f s is the tangential traction, f n is the normal traction, v s is the tangential velocity, and v n is the normal velocity. Note that the prescribed volumetric flow rate Q and the normal velocity of boundaries Γ 4 and Γ 5 are specified to simulate the large scale deposition flow rate of the Model 19 extruder, as in [38]. Further, v is the velocity vector at the free surface, and n is the free surface unit normal vector. A zero normal velocity boundary condition on the free surface is imposed by implementing the 1-D streamlinewise remeshing technique proposed by Tanner et al. [45], which assumes that the free surface boundary of the planar extrudate forms a streamline along the boundary of the flow domain. The direction of the polymer composite melt flow goes from a vertical to a horizontal flow in the deposition turning flow domain Ω 2 , making it helpful to divide the deposited bead into two subdomains Ω 2 and Ω 3 to simulate the free surface boundaries before and after the flow contacts the substrate, respectively, as shown in Figure 4a. The points P Ω 2 0 and P Ω 3 0 denotes the nodes at the intersection of the fixed-wall and free surface boundaries. To enforce the zero normal velocity on the free boundary nodes, we use the relation where ω and δ refer to the coordinates in the direction of flow and perpendicular to the direction of the flow, respectively. The subscript i indicates successive node numbers which start at 0 with the fixed node at the end of the no-slip wall [45], e.g., P Ω 2 0 , P Ω 2 1 , P Ω 2 2 appearing in Figure 4b. From above, it is seen that the ω coordinates of nodes in Ω 2 do not vary. Note, ω refers to different directions in Ω 2 and Ω 3 , due to the 90 degree turn of the deposition flow. The free surface is computed in an iterative manner where we convergence is assessed with the die swell ratio of the extrudate, which is written as where B donates the die swell ratio, δ (j) k refers to the normal-to-flow-direction coordinate of the k-th node of the free surface along the flow direction (i.e., typically the last node of the surface boundary which indicate a steady state) at the j-th iteration. The convergence of the flow surface is achieved when the error between the steady state die swell ratio values in two successive iterations B (j) − B (j−1) ≤ 10 −2 , from which the computed result shows good agreement with the free surface shape computed with a commercial finite element solver for the uncoupled flow problem as shown in the Appendix A. The δ coordinates of internal nodes are updated by a linear interpolation between those of the free surface and the fixed line of the subdomain (i.e., the intersection line between Ω 2 and Ω 3 is the fixed line for Ω 2 , and the contact surface between the deposition flow and the substrate serves as the fixed line for Ω 3 ). Additionally, the initial mesh of the flow domain is generated using ANSYS Polyflow module assuming a domain boundary without die swell. The mesh of the flow domain remains unchanged except that affected by the change in location of the free surface boundaries. The overall flow computation continues until both free surface shapes meet the convergence criterion described above.

Results and Discussions
In this study, we use the 13 wt.% CF/ABS as the polymer composite melt material model that has seen widespread applications in the large area polymer deposition additive manufacturing [2,27,36]. We assume the Newtonian viscosity of the melt is = 817 Pa•s, which is the shear viscosity of the ABS polymer at 230 ℃ and shear rate 100 , and the fluid density is = 1154 kg/m 3 , both of which are from Heller et al. [38]. In addition, we assume an isotropic fiber orientation field for the entire flow domain as a starting point for the fiber orientation calculation iterations when computing the fully coupled flow-orientation solutions, as in [44,50]. Moreover, the hyperbolic form of Equation (5) requires an initial condition of the second-order fiber orientation tensor at the flow inlet. Herein, we imposed the fiber orientation state of a fully developed flow, which assumes that the orientation state reaches steady state upstream of the flow inlet. The initial condi-

Results and Discussions
In this study, we use the 13 wt.% CF/ABS as the polymer composite melt material model that has seen widespread applications in the large area polymer deposition additive manufacturing [2,27,36]. We assume the Newtonian viscosity of the melt is η =817 Pa·s, which is the shear viscosity of the ABS polymer at 230 • C and shear rate 100 s −1 , and the fluid density is ρ = 1154 kg/m 3 , both of which are from Heller et al. [38]. In addition, we assume an isotropic fiber orientation field for the entire flow domain as a starting point for the fiber orientation calculation iterations when computing the fully coupled flow-orientation solutions, as in [44,50]. Moreover, the hyperbolic form of Equation (5) requires an initial condition of the second-order fiber orientation tensor at the flow inlet. Herein, we imposed the fiber orientation state of a fully developed flow, which assumes that the orientation state reaches steady state upstream of the flow inlet. The initial condition imposed at the flow inlet is fixed as A ini , such that which is computed using the uncoupled flow kinematics of a long tube flow (in the same diameter as the nozzle die inlet). Moreover, to reduce the singularity issue along the flow boundaries, fiber orientation along the no-slip wall is specified as fully aligned in the extrusion direction, and fiber-wall interaction is neglected as suggested by [64]. This fiber alignment assumption enhances the convergence behavior of the iterative process for the fiber orientation tensor field, while having little effect on the accuracy of velocity and fiber orientation solution [44]. We note that the boundary conditions specified above for the fiber orientation calculation are fixed throughout the iterative process. The fully developed velocity profile imposed at the flow inlet (cf. Figure 3) is computed in a separate analysis using the ANSYS Polyflow module based on the prescribed flow rate Q and the selected rheology model, and is fixed throughout the flow computation process. The flow domain in Figure 3 is meshed with 4-node quadrilateral elements using ANSYS. There are a total of 3550 elements with 3736 nodes in our flow model, which was found to provide a sufficient and efficient mesh quality through our mesh sensitivity study. Additionally, the fully coupled problem is solved in a decoupled fashion, where the flow and fiber orientation finite element matrices are computed independently with the counterpart fixed. The two sub-solutions are altered iteratively until both solutions converge, such that where " " refers to the Frobenius vector norm [60], and h is the iteration index of the overall coupling scheme. The convergence criterion appearing in Equation (22) indicates the absolute relative difference of two successive solutions are less than 1%, which yields sufficient accurate results as shown below. By comparing results between uncoupled and coupled solutions, mutually dependent effects between the flow and fiber orientation in the planar deposition flow can be clearly illustrated. Below, we consider the following computed outputs for comparison: flow kinematics throughout the flow fields, orientation tensor fields within the flow domain, and the mechanical properties at the flow end Γ 5 that also represents the mechanical performances of deposited beads. In addition, a parametric study is performed on C I to assess the sensitivity of the computed coupling effects with respect to the fiber-fiber interaction coefficient.

Flow Kinematics
Results from the uncoupled flow kinematics simulated with N p = 0 computed using our custom melt flow simulation presented above show good agreement with those obtained from the commercial finite element suite ANSYS Polyflow. For conciseness and completeness, the data comparison is given in the Appendix A. In the following, we discuss the results of the extrusion/deposition flow and planar extrudate swell computed using both our uncoupled and our fully coupled computational methods.
Contours of computed melt flow velocity v x and v y appear in Figures 5 and 6, respectively. Upon comparing the uncoupled and fully coupled results, significant variations can be seen in both v x and v y as the melt flow passes through the Nozzle Convergence Zone (NZR) and Deposition Transition Zone (DTZ). Further, contours of the velocity differences v c x − v unc x and v c y − v unc y appear in Figures 7 and 8, respectively, where the NCZ and DTZ are separately specialized (note superscripts c = coupled and unc = uncoupled). It can be seen that the maximum difference in flow contours appears in the DTZ of v x (cf. Figure 7b), where the fully coupled velocity is about 14 mm/s higher than that of the uncoupled solution. This velocity difference is significant since it is 13.8% of the deposition rate of the flow (i.e., 101.6 mm/s, as imposed at the flow end). The fiber orientation coupling effects also induce a notable change in v y near the boundaries of the convergence zone, while the center of the melt flow slows over 10 mm/s (as indicated in Figure 8a). Moreover, the polymer melt die swell is also affected by the coupling. The shape of the free surfaces of deposited beads appearing in Figure 9 show that the maximum die swell ratios (cf. Equation (20)) of the bead front swell profile are 1.31 and 1.24 as computed with the uncoupled and coupled formulation, respectively. That is, the presence of fibers reduces the bead front swell by~6% as compared to that of virgin ABS polymer deposited bead. In addition, the free surface profile of the turning flow (i.e., Ω 3 in Figure 3) predicted by the fully coupled formulation increases slightly as compared to that of the uncoupled solution and the die swell ratios at the flow end solved by the uncoupled and fully coupled formulations are 0.95 and 0.96, respectively. The computed results shown above indicate that processing parameters such as the extrusion flow rate, deposition rate as well as the inter-beads spacing for virgin polymer and their filled composite counterpart should be different, and proper adjustment are needed to compensate for the differences yielded by the fully coupled interactions between the flow kinematics and the fiber orientations. shape of the free surfaces of deposited beads appearing in Figure 9 show that the maximum die swell ratios (cf. Equation (20)) of the bead front swell profile are 1.31 and 1.24 as computed with the uncoupled and coupled formulation, respectively. That is, the presence of fibers reduces the bead front swell by ~6% as compared to that of virgin ABS polymer deposited bead. In addition, the free surface profile of the turning flow (i.e., Ω in Figure 3) predicted by the fully coupled formulation increases slightly as compared to that of the uncoupled solution and the die swell ratios at the flow end solved by the uncoupled and fully coupled formulations are 0.95 and 0.96, respectively. The computed results shown above indicate that processing parameters such as the extrusion flow rate, deposition rate as well as the inter-beads spacing for virgin polymer and their filled composite counterpart should be different, and proper adjustment are needed to compensate for the differences yielded by the fully coupled interactions between the flow kinematics and the fiber orientations.  shape of the free surfaces of deposited beads appearing in Figure 9 show that the maximum die swell ratios (cf. Equation (20)) of the bead front swell profile are 1.31 and 1.24 as computed with the uncoupled and coupled formulation, respectively. That is, the presence of fibers reduces the bead front swell by ~6% as compared to that of virgin ABS polymer deposited bead. In addition, the free surface profile of the turning flow (i.e., Ω in Figure 3) predicted by the fully coupled formulation increases slightly as compared to that of the uncoupled solution and the die swell ratios at the flow end solved by the uncoupled and fully coupled formulations are 0.95 and 0.96, respectively. The computed results shown above indicate that processing parameters such as the extrusion flow rate, deposition rate as well as the inter-beads spacing for virgin polymer and their filled composite counterpart should be different, and proper adjustment are needed to compensate for the differences yielded by the fully coupled interactions between the flow kinematics and the fiber orientations.

Second-Order Fiber Orientation Tensor Fields
Our fully coupled algorithm is first tested with the same parameter inputs as in Heller et al. [38]. A fixed time step is adopted in time marching scheme in order to obtain the steady state fiber orientation state in the flow starting from an isotropic fiber alignment (i.e., fully random fiber orientation state). Herein, the time increment is fixed at 0.01 which

Second-Order Fiber Orientation Tensor Fields
Our fully coupled algorithm is first tested with the same parameter inputs as in Heller et al. [38]. A fixed time step is adopted in time marching scheme in order to obtain the steady state fiber orientation state in the flow starting from an isotropic fiber alignment (i.e., fully random fiber orientation state). Herein, the time increment is fixed at 0.01 which

Second-Order Fiber Orientation Tensor Fields
Our fully coupled algorithm is first tested with the same parameter inputs as in Heller et al. [38]. A fixed time step is adopted in time marching scheme in order to obtain the steady state fiber orientation state in the flow starting from an isotropic fiber alignment (i.e., fully random fiber orientation state). Herein, the time increment is fixed at 0.01 which

Second-Order Fiber Orientation Tensor Fields
Our fully coupled algorithm is first tested with the same parameter inputs as in Heller et al. [38]. A fixed time step is adopted in time marching scheme in order to obtain the steady state fiber orientation state in the flow starting from an isotropic fiber alignment (i.e., fully random fiber orientation state). Herein, the time increment is fixed at 0.01 which was found to converge well and to be computationally efficient. Convergence is assumed when the absolute relative difference (cf. Equation (22)) of orientation tensor solutions between two successive time steps is less than 5 × 10 −3 . We note that convergence issues were experienced when imposing a tighter convergence criterion while the fiber orientation tensor solution at the flow end with current convergence criterion show a good agreement with those explicitly given in Heller et al. [38] (cf. Figure 10, where "myfea" refers to results obtained by our finite-element-based algorithm). Therefore, it can be expected that the adopted convergence criterion for the fiber orientation solution yields reliable numerical results and thus is applied in the following studies. was found to converge well and to be computationally efficient. Convergence is assumed when the absolute relative difference (cf. Equation (22)) of orientation tensor solutions between two successive time steps is less than 5 × 10 . We note that convergence issues were experienced when imposing a tighter convergence criterion while the fiber orientation tensor solution at the flow end with current convergence criterion show a good agreement with those explicitly given in Heller et al. [38] (cf. Figure 10, where "myfea" refers to results obtained by our finite-element-based algorithm). Therefore, it can be expected that the adopted convergence criterion for the fiber orientation solution yields reliable numerical results and thus is applied in the following studies. For the material model of this study, we consider a CF/ABS composite having a weight fraction of 13% (i.e., the volume fraction of carbon fiber is roughly = 8.4% [35]), a constant fiber aspect ratio = 30, as defined in related literature (i.e., mean fiber length is 150 μm [65], and a mean fiber diameter of 5 μm [27]). Consequently, , , and are determined as 0.9950, 0.0055, and 13.92, respectively, using Equations (8)- (11). Computed contours of second order orientation tensor components , , and are given in Figures 11-13, respectively. Note that the fiber orientation tensor components and change significantly as the flow passed through the DTZ. To make the variations within a certain flow region more obvious, the color scale of the presented contours is adjusted. It is seen that the fibers are highly aligned along the direction of the flow downstream of the nozzle convergence zone. Specifically, the vortex in the nozzle convergence corner grows larger in the fully coupled simulation as compared to the uncoupled solution. In addition, there is a notable difference between the uncoupled and fully coupled solution near the wall boundaries of the NCZ in both the results of and (cf. Figures 11 and 12), which may be a result of the pre-defined wall-alignment fiber orientation boundary condition (cf. [64]). For the material model of this study, we consider a CF/ABS composite having a weight fraction of 13% (i.e., the volume fraction of carbon fiber is roughly v f = 8.4% [35]), a constant fiber aspect ratio g r =30, as defined in related literature (i.e., mean fiber length is 150 µm [65], and a mean fiber diameter of 5 µm [27]). Consequently, λ, C I , and N p are determined as 0.9950, 0.0055, and 13.92, respectively, using Equations (8)- (11). Computed contours of second order orientation tensor components A 22 , A 12 , and A 11 are given in Figures 11-13, respectively. Note that the fiber orientation tensor components A 22 and A 12 change significantly as the flow passed through the DTZ. To make the variations within a certain flow region more obvious, the color scale of the presented contours is adjusted. It is seen that the fibers are highly aligned along the direction of the flow downstream of the nozzle convergence zone. Specifically, the vortex in the nozzle convergence corner grows larger in the fully coupled simulation as compared to the uncoupled solution. In addition, there is a notable difference between the uncoupled and fully coupled solution near the wall boundaries of the NCZ in both the results of A 22 and A 12 (cf. Figures 11 and 12), which may be a result of the pre-defined wall-alignment fiber orientation boundary condition (cf. [64]).    The fiber orientation formation in the DTZ is of special interest. To this end, we plot the strain rate contours (i.e., velocity gradients) in Figures 14-17 to better explain the fiber alignment variations that occur during the transition between extrusion and planar deposition. It is clearly seen that the fiber alignment changes rapidly in the DTZ, which is a result of the extensional flow forming during the transition from extrusion to deposition (cf. Figures 14a and 17a, where significant variations appear in the DTZ). Figure 13 indicates that the fully coupled solution and the uncoupled solutions of A 11 exhibit a similar orientation pattern as the deposition flow reaching a quasi-steady-state, except that the upper region in the DTZ exhibits an increment of~0.1 in the fully coupled solution (cf. Figure 13c). This is expected to result from the non-uniform shear force that is applied to the melt as the extrudate contacts the substrate. This change is captured by the fully coupled simulation (cf. the strain rate variation along the y-direction, i.e., Figures 15 and 17, where large increments appear in the contacting region). Fibers near the free surface boundary of the deposition flow (i.e., free surface of Ω 3 ) recover to high alignment along flow direction, which is likely due to that the deposition flow evolves into a shear dominate flow as shown in Figure 15). Similar results are also reported in Ouyang et al. [43]. Specifically, variations between the uncoupled and coupled solutions can be seen near the DTZ (cf. Figure 13c), where the fully coupled solution exhibits a reduced fiber alignment in the front of the deposition flow and the values of A 11 increase notably in the middle of the DTZ as the flow is redirected towards front of the deposition. The noticeable differences between the uncoupled and fully coupled fiber orientation solutions within the extrusion-deposition flow domain imply the importance of the fully coupled flow/fiber-orientation simulation in interpreting the process-structure-properties mapping of PDAM applications of composite materials, as suggested in [66]. The fiber orientation formation in the DTZ is of special interest. To this end, we plot the strain rate contours (i.e., velocity gradients) in Figures 14-17 to better explain the fiber alignment variations that occur during the transition between extrusion and planar deposition. It is clearly seen that the fiber alignment changes rapidly in the DTZ, which is a result of the extensional flow forming during the transition from extrusion to deposition (cf. Figures 14a and 17a, where significant variations appear in the DTZ). Figure 13 indicates that the fully coupled solution and the uncoupled solutions of exhibit a similar orientation pattern as the deposition flow reaching a quasi-steady-state, except that the upper region in the DTZ exhibits an increment of ~0.1 in the fully coupled solution (cf. Figure 13c). This is expected to result from the non-uniform shear force that is applied to the melt as the extrudate contacts the substrate. This change is captured by the fully coupled simulation (cf. the strain rate variation along the y-direction, i.e., Figures 15 and 17, where large increments appear in the contacting region). Fibers near the free surface boundary of the deposition flow (i.e., free surface of Ω ) recover to high alignment along flow direction, which is likely due to that the deposition flow evolves into a shear dominate flow as shown in Figure 15). Similar results are also reported in Ouyang et al. [43]. Specifically, variations between the uncoupled and coupled solutions can be seen near the DTZ (cf. Figure 13c), where the fully coupled solution exhibits a reduced fiber alignment in the front of the deposition flow and the values of increase notably in the middle of the DTZ as the flow is redirected towards front of the deposition. The noticeable differences between the uncoupled and fully coupled fiber orientation solutions within the extrusion-deposition flow domain imply the importance of the fully coupled flow/fiber-orientation simulation in interpreting the process-structure-properties mapping of PDAM applications of composite materials, as suggested in [66].

Fiber-Wall Interaction in Nozzle Convergence Zone
Rosenberg and Denn [67] (also see VerWeyst and Tucker [50]) stated that fiber-wall interaction at the no-slip wall boundary of a flow suspension causes a singularity where the finite length fibers exhibit a periodic tumbling motion within the flow and the fiber orientation tensor yields (e.g., Equation (5)) a unique pattern at the wall boundary. To reduce the intensity of this singularity, we set the fiber orientation to align in the direction of the flow at the wall boundaries as suggested by Ranganathan and Advani [64]. Unfortunately, the flow contraction in the nozzle convergence zone exhibited numerical instabilities in the fiber orientation tensor solution, especially in the fully coupled simulation (cf. Figures 11b and 12b). Figure 18 shows a plot of fiber orientation from the uncoupled and fully coupled solutions near the nozzle convergence zone, where the length and direction of each line segment is determined by the largest eigenvalue of the second order orientation tensor and its associated eigenvector, respectively [50]. In the uncoupled solution, fibers near the NCZ boundary are mis-aligned, especially in the vicinity of the nozzle exit. While in the case of the fully coupled simulation, fiber exhibit a higher misalignment

Fiber-Wall Interaction in Nozzle Convergence Zone
Rosenberg and Denn [67] (also see VerWeyst and Tucker [50]) stated that fiber-wall interaction at the no-slip wall boundary of a flow suspension causes a singularity where the finite length fibers exhibit a periodic tumbling motion within the flow and the fiber orientation tensor yields (e.g., Equation (5)) a unique pattern at the wall boundary. To reduce the intensity of this singularity, we set the fiber orientation to align in the direction of the flow at the wall boundaries as suggested by Ranganathan and Advani [64]. Unfortunately, the flow contraction in the nozzle convergence zone exhibited numerical instabilities in the fiber orientation tensor solution, especially in the fully coupled simulation (cf. Figures 11b and 12b). Figure 18 shows a plot of fiber orientation from the uncoupled and fully coupled solutions near the nozzle convergence zone, where the length and direction of each line segment is determined by the largest eigenvalue of the second order orientation tensor and its associated eigenvector, respectively [50]. In the uncoupled solution, fibers near the NCZ boundary are mis-aligned, especially in the vicinity of the nozzle exit. While in the case of the fully coupled simulation, fiber exhibit a higher misalignment

Fiber-Wall Interaction in Nozzle Convergence Zone
Rosenberg and Denn [67] (also see VerWeyst and Tucker [50]) stated that fiber-wall interaction at the no-slip wall boundary of a flow suspension causes a singularity where the finite length fibers exhibit a periodic tumbling motion within the flow and the fiber orientation tensor yields (e.g., Equation (5)) a unique pattern at the wall boundary. To reduce the intensity of this singularity, we set the fiber orientation to align in the direction of the flow at the wall boundaries as suggested by Ranganathan and Advani [64]. Unfortunately, the flow contraction in the nozzle convergence zone exhibited numerical instabilities in the fiber orientation tensor solution, especially in the fully coupled simulation (cf. Figures 11b and 12b). Figure 18 shows a plot of fiber orientation from the uncoupled and fully coupled solutions near the nozzle convergence zone, where the length and direction of each line segment is determined by the largest eigenvalue of the second order orientation tensor and its associated eigenvector, respectively [50]. In the uncoupled solution, fibers near the NCZ boundary are mis-aligned, especially in the vicinity of the nozzle exit. While in the case of the fully coupled simulation, fiber exhibit a higher misalignment where the flow wall boundary above the NCZ also shows such effect in addition to that of the NCZ. These results indicate that the flow-fiber coupling effect promotes the vortex generation in contraction region (i.e., corner of a converging flow as indicated by the quad block appearing in Figure 14). It is noted that a similar flow geometry (i.e., 4.5:1 contraction flow without nozzle convergence zone) was studied in Lipstomb et al. [68] and VerWeyst and Tucker [50] where the corner vortex exhibited a higher increment when the fully coupled effect was considered, which compares well to our result in Figure 18.
This implies that flow in the nozzle convergence zone may be a factor in reducing the corner vortex phenomenon of the fiber suspension, especially when highly concentrated composites are employed. We expect that the relaxed convergence criterion adopted for the fiber orientation iterative solution (i.e., error of orientation tensor solutions between two successive time steps reduces below 5 × 10 −3 ) likely results in the numerical instability in the computed solution near the nozzle convergence zone boundaries (see e.g., Figure 11b). In addition to enhance the convergence behavior of the simulation, we also check the converged solution at each time step, where the non-physical orientation tensor eigenvalues (i.e., diagonal terms of the tensor is above one) are adjusted to be within the physical bounds (i.e., within an interval from zero to one) and the associated eigenvectors are modified as well as in [50].
where the flow wall boundary above the NCZ also shows such effect in addition to that of the NCZ. These results indicate that the flow-fiber coupling effect promotes the vortex generation in contraction region (i.e., corner of a converging flow as indicated by the quad block appearing in Figure 14). It is noted that a similar flow geometry (i.e., 4.5:1 contraction flow without nozzle convergence zone) was studied in Lipstomb et al. [68] and Ver-Weyst and Tucker [50] where the corner vortex exhibited a higher increment when the fully coupled effect was considered, which compares well to our result in Figure 18. This implies that flow in the nozzle convergence zone may be a factor in reducing the corner vortex phenomenon of the fiber suspension, especially when highly concentrated composites are employed. We expect that the relaxed convergence criterion adopted for the fiber orientation iterative solution (i.e., error of orientation tensor solutions between two successive time steps reduces below 5 × 10 ) likely results in the numerical instability in the computed solution near the nozzle convergence zone boundaries (see e.g., Figure 11b). In addition to enhance the convergence behavior of the simulation, we also check the converged solution at each time step, where the non-physical orientation tensor eigenvalues (i.e., diagonal terms of the tensor is above one) are adjusted to be within the physical bounds (i.e., within an interval from zero to one) and the associated eigenvectors are modified as well as in [50].

Sensitivity on CI Parameter
The fiber-fiber interaction coefficient in Equation (5) has a significant influence on the fiber orientation prediction, especially in fully coupled flow-fiber simulations [50]. Equation (9) proposed by Bay [54], which we use to compute , empirically determined from injection molding studies, has not been fully validated for polymer deposition. To explore the sensitivity of the fully coupled simulation results to the fiber interaction coefficient, we consider various values of in the interval from 0.001 to 0.01, which is typical region for short fiber polymer composites studies [33][34][35][36][37][38]. Computed orientation tensor component values for 0.001 ≤ ≤ 0.01 at the nozzle exit and the flow end appear in Figures 19-24, respectively, in which the standard deviations among the numerical data are evaluated as well. In all cases, computed results show that fibers are highly aligned along the flow direction. Figure 21

Sensitivity on CI Parameter
The fiber-fiber interaction coefficient C I in Equation (5) has a significant influence on the fiber orientation prediction, especially in fully coupled flow-fiber simulations [50]. Equation (9) proposed by Bay [54], which we use to compute C I , empirically determined from injection molding studies, has not been fully validated for polymer deposition. To explore the sensitivity of the fully coupled simulation results to the fiber interaction coefficient, we consider various values of C I in the interval from 0.001 to 0.01, which is typical region for short fiber polymer composites studies [33][34][35][36][37][38]. Computed orientation tensor component values for 0.001 ≤ C I ≤ 0.01 at the nozzle exit and the flow end appear in Figures 19-24, respectively, in which the standard deviations among the numerical data are evaluated as well. In all cases, computed results show that fibers are highly aligned along the flow direction. Figure 21 (i.e., A 22 refers to flow-direction alignment) and Figure 22 (i.e., A 11 refers to flow-direction alignment) show that C I = 0.001 yields the highest flowdirection alignment among all simulations, which implies that an increased fiber-fiber interaction reduces fiber rotation that aligns fibers the direction of flow. While the in-nozzle fiber orientation is symmetric with respect to the central axis of the nozzle, it is clearly seen that the orientation state across the bead is not, and instead changes as the melt comes into contact with the substrate. As a result, flow-direction alignment is higher in the lower region of the deposition flow (i.e., subdomain Ω 3 ) than along the upper surface as seen in . This pattern of fiber alignment is due in part to the shear force applied by the moving substrate. By reducing the value of C I , the flow-direction fiber alignment increases which can be as high as that near the substrate-contact boundary for C I = 0.001. In addition, the standard deviations appearing in the Figures 19-24 indicate that the employed value of C I mostly affects the flow-direction alignment, as evidenced by the deviation of~0.07 for A 22 at nozzle exit and~0.10 for A 11 at deposition flow end.
fiber interaction reduces fiber rotation that aligns fibers the direction of flow. While the in-nozzle fiber orientation is symmetric with respect to the central axis of the nozzle, it is clearly seen that the orientation state across the bead is not, and instead changes as the melt comes into contact with the substrate. As a result, flow-direction alignment is higher in the lower region of the deposition flow (i.e., subdomain Ω ) than along the upper surface as seen in Figures 22-24. This pattern of fiber alignment is due in part to the shear force applied by the moving substrate. By reducing the value of , the flow-direction fiber alignment increases which can be as high as that near the substrate-contact boundary for = 0.001. In addition, the standard deviations appearing in the Figures 19-24 indicate that the employed value of mostly affects the flow-direction alignment, as evidenced by the deviation of ~0.07 for at nozzle exit and ~0.10 for at deposition flow end.
(a) (b) Figure 19.  in-nozzle fiber orientation is symmetric with respect to the central axis of the nozzle, it is clearly seen that the orientation state across the bead is not, and instead changes as the melt comes into contact with the substrate. As a result, flow-direction alignment is higher in the lower region of the deposition flow (i.e., subdomain Ω ) than along the upper surface as seen in Figures 22-24. This pattern of fiber alignment is due in part to the shear force applied by the moving substrate. By reducing the value of , the flow-direction fiber alignment increases which can be as high as that near the substrate-contact boundary for = 0.001. In addition, the standard deviations appearing in the Figures 19-24 indicate that the employed value of mostly affects the flow-direction alignment, as evidenced by the deviation of ~0.07 for at nozzle exit and ~0.10 for at deposition flow end.
(a) (b) Figure 19.       Finally, the fiber orientation state in the polymer melt flow is a pivotal factor in determining the material properties of the deposited beads. To assess the effect of processing on mechanical properties, we apply the orientation-homogenization approach [46] to compute the material stiffness and thermal expansion coefficient of the composite bead using second order fiber orientation tensor values computed at the flow end (e.g., computed results appearing in Figures 22-24). The flow end is chosen here since this location represents the steady-state orientation pattern of a deposited bead. Properties of the constituent materials used in this study are given in Table 1. Furthermore, the computed properties at the nodal positions are numerically integrated along the y-axis over the end Finally, the fiber orientation state in the polymer melt flow is a pivotal factor in determining the material properties of the deposited beads. To assess the effect of processing on mechanical properties, we apply the orientation-homogenization approach [46] to compute the material stiffness and thermal expansion coefficient of the composite bead using second order fiber orientation tensor values computed at the flow end (e.g., computed results appearing in Figures 22-24). The flow end is chosen here since this location represents the steady-state orientation pattern of a deposited bead. Properties of the constituent materials used in this study are given in Table 1. Furthermore, the computed properties at the nodal positions are numerically integrated along the y-axis over the end of the bead in order to obtain the effective mechanical property of interest. More detail on material property prediction method used here can be found in [33]. Computed results are given in Tables 2-4, where the effective axial and shear moduli are computed based on local predicted material stiffness values. The Partial Relative Differences (PRD) among computed data are considered here, where the properties evaluated by employing C I = 0.001 is designated as the reference data. The PRD may be evaluated as [37]   The resulting PRDs are given in Tables 5-7, where the differences yielded by varying C I values clearly shown. It can be seen that the elastic constant E 11 is most influenced by the value of C I , since E 11 is the effective axial modulus along the bead which is also the major fiber alignment direction. Specifically, the PRD between E 11 computed with C I = 0.01 and that computed with C I = 0.001 is 23.16%. It is also seen that the impact of C I choice is larger for elastic constants than the thermal expansion coefficient, where the PRDs for moduli (cf. Tables 5 and 6) are~10-20%, and those in α are less than 10% (cf. Table 7). In addition, the transverse properties exhibit some skewness, especially for the computed thermal expansion coefficients. For example, in the case of C I = 0.01, α 22 is 13% higher than α 33 . Uneven thermal expansions transverse to the bead axis may be a critical factor affecting the inter-bead void formation and distortion during manufacturing. Additionally, we note that our computed values of effective modulus E 11 show good agreement with experimentally reported values on the same material system with slightly different processing parameters (i.e., measured tensile modulus parallel to the printing direction of a 13 wt.% CF/ABS bead is reported as 8.18 GPa [27]). In addition, the predicted y-axis thermal expansion coefficient (i.e., α 22 ) properties also agree well with reported data for the same material system in a related study (i.e., 7.41 × 10 −5 / • C for FEA results and 5.77 × 10 −5 / • C for DIC measurement, see Table 1. in [31]). Herein, we consider the above favorable agreements as a support to the computed results from our fully coupled flow-fiber simulation. Nevertheless, we note that the employed fiber orientation constitutive equations (i.e., Equations (5)-(9)) describe the fiber-fiber interaction through a phenomenological manner, which we believe provides valuable insight into the flow of fibers in the large area additive manufacturing polymer composite deposition process. More advanced physical-based approaches have yet to be developed. It is important to note that, beyond our simulation, the microstructural voids formation [69], intense fiberfiber interaction in highly concentrated fiber composites [65], and reduced fiber length distribution [36] can also significantly affect the material properties of polymer deposition additive manufactured composites.

Summary
This paper numerically characterizes the mutually dependent effects between polymer melt kinematics and suspended fiber orientation in polymer deposition additive manufacturing of short fiber reinforced composites. The numerical solution of the governing equations for flow and fiber orientation problems are computed iteratively using the finite element method. Calculated results indicate that the uncoupled and fully coupled flow field solutions exhibit significant differences near the nozzle convergence zone and the extrusion-deposition transition zone where intense flow/orientation mutual interaction occurs. The orientation tensor component A 22 exhibits notable variation in between the uncoupled and fully coupled solution along the nozzle convergence zone boundaries, where the fully coupled results experience significant numerical instability in the flow contraction region. The pattern of the A 11 component computed from both simulations are similar, where the largest difference is seen at the extrusion-deposition transition zone. The sensitivity study on C I indicates that fiber-fiber interaction is of great importance in determining the amount of fiber alignment along the flow domain, especially for the steady-state orientation at the flow end, which results in a 23% PRD error in E 11 computed using C I = 0.01 and C I = 0.001. These results indicate that an in-depth numerical and experimental combined study on C I would be of great value for further explaining the fully coupled flow/orientation features in polymer deposition additive manufacturing applications. In addition, the transverse thermal expansion coefficients exhibit notable degree of skewness which may lead to uneven distortion during manufacturing. Lastly, our predicted material properties of 13 wt.% CF/ABS show favorable agreement with experiment data reported in the literature on the same material model which support the simulation approach and associated computed results.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Flow Computation Verification
Flow kinematics of the extrusion-deposition flow presented in Figure 3 are separately solved using the ANSYS Polyflow and our formulated finite element algorithm (i.e., assigning N p = 0 with no iterative solutions between the flow and fiber orientation problems). Computed results appear in Figures A1-A4 are obtained using two simulation methods, where the ANSYS Polyflow and our finite element solutions are referred as "Ansys" and "myfea", respectively. Note, the locations indicated in legends of Figures A1-A3 can be referred to Figures 9 and A4 in above. The data presented here shows that the velocity profiles and predicted extrudate swell profile computed by our finite element simulation are in good agreement with those obtained with ANSYS Polyflow.
The computed results presented here also show that flow velocities near the free surface boundaries differ between the two simulation results. This is due in part to the onedimensional streamline free surface remeshing approach we adopted in our simulations which is not as effective in resolving the singularity issue that occurs at the intersection of the no-slip wall boundary and the free surface boundary as compare to the related algorithm in ANSYS Polyflow. However, we expect the singularity issue has minimal effects on free surface identification as most of the flow field results match ANSYS Polyflow well as shown through Figures A1-A3. In addition, Figure 10 shows that the fiber orientation results predicted with our uncoupled flow kinematics show a good agreement with prior literature.
which is not as effective in resolving the singularity issue that occurs at the intersection of the no-slip wall boundary and the free surface boundary as compare to the related algorithm in ANSYS Polyflow. However, we expect the singularity issue has minimal effects on free surface identification as most of the flow field results match ANSYS Polyflow well as shown through Figures A1-A3. In addition, Figure 10 shows that the fiber orientation results predicted with our uncoupled flow kinematics show a good agreement with prior literature.   which is not as effective in resolving the singularity issue that occurs at the intersection of the no-slip wall boundary and the free surface boundary as compare to the related algorithm in ANSYS Polyflow. However, we expect the singularity issue has minimal effects on free surface identification as most of the flow field results match ANSYS Polyflow well as shown through Figures A1-A3. In addition, Figure 10 shows that the fiber orientation results predicted with our uncoupled flow kinematics show a good agreement with prior literature.