Effects of extrudate swell and nozzle geometry on fiber orientation in Fused Filament Fabrication nozzle flow

Accepted Manuscript Title: Effects of extrudate swell and nozzle geometry on fiber orientation in Fused Deposition Modeling nozzle flow Author: Blake P. Heller Douglas E. Smith David A. Jack PII: DOI: Reference:

S2214-8604(16)30116-6 http://dx.doi.org/doi:10.1016/j.addma.2016.06.005 ADDMA 103

To appear in: Received date: Revised date: Accepted date:

31-12-2015 4-5-2016 1-6-2016

Please cite this article as: {http://dx.doi.org/ This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Effects of Extrudate Swell and Nozzle Geometry on Fiber Orientation in Fused Deposition Modeling Nozzle Flow

Blake P. Heller, Douglas E. Smith1, David A. Jack

Department of Mechanical Engineering, Baylor University, One Bear Place #97356, Waco, TX, 767987356 [email protected], [email protected], [email protected]

1

Corresponding Author

Abstract The viability of Fused Deposition Modeling (FDM) as an Additive Manufacturing (AM) method is improved by including discrete carbon fibers in the polymer feedstock. The resulting carbon fiber polymer composite enhances mechanical and thermal properties, where improvements are directly related to the orientation of the carbon fibers within the extruded polymer composite. Fiber orientation in the polymer melt is affected by the flow field which is defined by the melt flow geometry including flow in the nozzle convergence and extrudate swell. This paper presents a computational approach for modeling the polymer melt flow and the associated effects of nozzle geometry and extrudate swell on fiber orientation in the extruded polymer. The finite element method is used to evaluate Stokes flow in the axisymmetric flow field of an FDM nozzle. Fiber orientation tensors are calculated throughout the fluid domain using the Orthotropic Fitted Closure and the Folgar-Tucker isotropic fiber interaction model. The radial and axial modulus of the extruded polymer is then calculated from fiber orientation tensors within the composite at the nozzle exit. Fiber orientation along the extrudate axis is shown to increase due to elongation flow within the nozzle convergence zone and decrease in the extrudate swell such that the modulus of elasticity of the extruded polymer is moderately affected by extrudate swell and substantially affected by the nozzle geometry.

Keywords: Fused Deposition Modeling, Fiber Orientation, Additive Manufacturing, Polymer Composite, Elastic Properties

Nomenclature v – velocity vector 𝛔 – stress tensor p – fluid pressure μ – Newtonian fluid viscosity n – normal vector rne – radius of the nozzle exit S(α) – extrudate swell free surface α – radial parameters defining the free surface 𝜎n – normal stress 𝜎t – tangential stress p – unit vector along fiber axis Ω – vorticity tensor Γ – rate of deformation tensor λ – coefficient related to fiber geometry AR – cylindrical fiber aspect ratio AR* – equivalent aspect ratio l – length of fiber d – diameter of fiber ψ(p) – probability density function Dr – diffusion function Aij, 𝔸ijkl – fiber orientation tensors CI – fiber interaction constant 𝛾̇ – scalar magnitude of the rate of deformation tensor VF – volume fraction of fibers ̃ 𝑂𝑅𝑇 𝔸 𝑚𝑚 – orthotropic fitted closure approximation of the fourth order orientation tensor 𝜆1 , 𝜆2 – eigenvalues of the second order orientation tensor 1 15 𝐶𝑚 -𝐶𝑚 – fitted coefficients from different flow types 〈𝐶𝑖𝑗𝑘𝑙 〉 – equivalent stiffness matrix C11, C12, C22, C23, and C66 – elasticity tensor coefficients for unidirectional composites 𝑑𝑣 𝑑𝑣 𝑣 𝑑𝑣 𝑑𝑣 𝑣𝑟 , 𝑣𝑧 , 𝑟 , 𝑟 , 𝑟 , 𝑧 , 𝑧 – velocity gradients 𝑑𝑟

𝑑𝑧

𝑟

𝑑𝑟

𝑑𝑧

𝑣𝑧𝑖𝑛𝑙𝑒𝑡 , 𝑣𝑟𝑖𝑛𝑙𝑒𝑡 , 𝑎𝑛𝑑 𝑣𝜃𝑖𝑛𝑙𝑒𝑡 – inlet velocity components 𝑣𝑚𝑎𝑥 – maximum velocity 𝑟𝑖 – radial location at inlet 𝑟𝑚𝑎𝑥 – max radius at inlet ρ – fluid density CZS – start of convergence zone CZE – end of convergence zone NE – nozzle exit STE – straight tube end 𝐴33 – average fiber alignment 𝐸𝑓 , 𝐸𝑚 – modulus of elasticity of the fiber and matrix 𝑣𝑓 , 𝑣𝑚 – Poisson’s ratio of the fiber and matrix 𝐺𝑟𝑧 - Shear modulus of the fiber filled composite

1 Introduction Fused Deposition Modeling (FDM) is transitioning from a small scale rapid prototyping technology to large scale additive manufacturing of industrially viable parts and tooling. FDM is used to fabricate cost-effective, intricate three dimensional parts that can be created with little wasted material and reduced energy consumption. Recently the addition of discrete carbon fibers to ABS filaments has been shown to greatly improve the mechanical properties of the printed part [1]. The carbon fiber filled ABS provides increased strength and toughness and a decreased CTE which allows for better dimensional stability, decreased warping, and reduced delamination in the printed part [1]. The improvement of mechanical properties is directionally dependent on the orientation of the discrete fibers suspended in the polymer melt making it is necessary to accurately predict the fiber orientation state in the extruded polymer. Little research has been done on the orientation of discrete fibers in the FDM printing process. Nixon, et al. [2] recently predicted fiber orientation in three FDM nozzle configurations (straight, converging, and diverging) within Moldflow (Moldflow Corporation, Framingham, MA) using the AdvaniTucker [3] orientation tensor approach. Unfortunately, their approach did not consider the extruded polymer beyond the nozzle exit and they neglected the evaluation of extrudate swell and its effects on the resulting fiber orientation state. Fiber orientation modeling in injection and compression molding has been an area of substantial research in the previous decades and is well understood. Software packages such as Moldflow (Moldflow Corporation, Framingham, MA) and Moldex3D (Core Tech Systems Co., Ltd., Chupei City, Taiwan) simulate these molding processes and predict fiber orientation tensor values using methods first proposed by Advani and Tucker [3]. The effect of fiber orientation on melt flow rheology is commonly not considered for the thin cavity, shear dominated flows found in injection molding [4]. However, full coupling between melt flow and fiber orientation in three-dimensional extrusion die problems is important as seen in the simulation of rods and tubular channels (see e.g., Lee and Lee [5], VerWeyst and Tucker [6], Baloch [7] and Ausias, et al. [8]). Of particular interest here is the work by Lee and Lee [5] who simulated die swell of a fiber filled polymer from a square die using a fully coupled simulation of a Newtonian fluid. They showed that the swell ratio decreased considerably when suspended fibers are present, and that the fiber orientation state at the die inlet was important. Decreases in die swell due to the presence of suspended fibers has been seen experimentally in various studies as well (see e.g., Uematsu, et al., [9], Wagner, et al., [10], Stabik [11], and Becraft and Metzner [12]), further supporting the need for a fully coupled simulation approach in future work. The swell of an extrudate at the die exit has received much attention in previous decades. Georgiou, et al. [13] and Reddy and Tanner [14] studied extrudate swell for Newtonian fluids exiting a die and have shown that a swell of 13% is expected from the original exit diameter. More recently, Mitsoulis, et al., [15] used a similar numerical approach with a Newtonian fluid model to show that gravity, inertia, slip at the wall, surface tension, and pressure dependence of viscosity all tend to decrease swell while melt compressibility increases it after a small initial reduction. Mitsoulis and Hatzikiriakos [16], and Kanvisi, et al., [17] have simulated die swell using viscoelastic material models to include the elastic effect in polymer melts beyond the die

exit. These works predict a larger swell ratio as compared to the Newtonian fluid predictions, as expected. However, simulations of the die swell with a viscoelastic fluid coupled with fiber orientation have yet to be solved. The effect of fiber orientation on directionally dependent mechanical properties has been studied in great depth for short fiber reinforced composites. The majority of the research focuses on unidirectional short fiber models [18,1 9], but there exist methods that can predict elastic mechanical properties for non-aligned discontinuous fibers using the associated unidirectional properties [3]. Advani and Tucker proposed an orientation homogenization method which uses associated unidirectional stiffness tensors to define the stiffness tensor of randomly oriented fibers using the homogenization method which is implemented by Camacho [20]. This method was later fully derived and extended using the spherical harmonic expansion approach by Jack and Smith [21]. This paper presents a computational method for predicting fiber orientation in an axisymmetric FDM nozzle including the evaluation of Newtonian extrudate swell at the nozzle exit. The effects of the FDM nozzle convergence zone and the extrudate swell at the nozzle exit on the fiber orientation state and, therefore, the resulting material stiffness of the extruded polymer, are the focus of the present study. Here we assume that the fiber length and volume fraction of the fibers are such that nozzle clogging can be neglected. The computationally efficient fiber orientation tensor approach from Advani and Tucker [3] is used with Folgar-Tucker [22] isotropic rotary diffusion for modeling fiber orientation. The Orthotropic Fitted Closure [23, 24] is employed to address the orientation tensor closure problem. Fiber orientation is calculated through an FDM nozzle with geometry similar to that of most commercially available desktop FDM printers, such as the Makerbot 2X (Makerbot Industries, LLC, Brooklyn, NY). A parametric study of the nozzle geometry is then presented to more fully understand the effect of nozzle geometry on fiber orientation. In addition, fiber orientation-dependent mechanical properties are calculated at the nozzle exit and just outside the nozzle exit in the extrudate swell to quantify the effect of the extrudate swell. The Tandon and Weng [19] approach for calculating uniaxial compliance tensor components is used along with the orientation homogenization method by Jack and Smith [21] to predict the axial and radial moduli of the extruded composite. This work is a first attempt to model fluid flow, fiber orientation, and predict elastic properties in the FDM nozzle flow of a polymer composite. Assumptions made to provide a reasonable initial approach to modeling this phenomena include the steady flow of a Newtonian fluid, uncoupled fiber orientation calculation, and an axisymmetric fluid domain which does not include the contact of the polymer bead with the moving platform as is typical in the FDM process. As a result, details associated with bead deposition such as changes in speed are not considered here. The effect of suspended fibers on flow rheology and extrudate swell, and fluid elasticity will also be the subject of future model development.

2 Methodology 2.1 Polymer Melt Flow and Extrudate Swell Models

The fluid in the axisymmetric FDM nozzle geometry shown in Figure 1 is modeled as an incompressible Stokes flow. Stokes flow neglects the inertial term in the Navier-Stokes equation and is used for flows with low Reynolds number and high viscosity such as those found in polymer melt flows. The continuity and momentum equations for Stokes flow are given respectively as ∇∙𝐯 = 0

(1)

and ∇∙𝝈=0 In the above, 𝐯 is the fluid velocity vector and 𝝈 is the stress tensor given as

(2)

(3) 𝝈 = −𝑝𝑰 + 𝜇[(∇𝐯) + (∇𝐯)T ] where 𝑝 is the pressure, 𝑰 is the identity tensor, and 𝜇 is the Newtonian viscosity of the fluid.

The shape of the extrudate swell that occurs immediately after the nozzle exit is dependent on the stress state within and on the free surface of the fluid. The free surface of the extrudate swell is defined by a zero surface stress (𝐧 ∙ 𝝈 ∙ 𝐧 = 0) and a zero fluid penetration condition (𝐧 ∙ 𝐯 = 0) as shown in Figure 1. The extrudate swell free surface shape can be computed by setting either the zero surface stress (see e.g., [14]) or zero fluid penetration condition (see e.g., [13]) using a minimization process to reduce the complimentary condition to zero. In the current study, a zero fluid penetration condition is defined on the free boundary and the free surface stress is minimized by adjusting the shape of the flow domain as shown in Figure 1. The normal and tangential stress components are calculated on element edges along the free surface using the finite element solution to Eqs. (3) and (4). The objective function and constraints for the minimization problem are 𝑚𝑖𝑛 𝛼 𝑓(𝛼) = ∫𝑆 (𝜎𝑛 (𝑆(𝛼)) + 𝜎𝑡 (𝑆(𝛼))) 𝑑𝑧 = 0 (4) subject to

1.2 ∗ 𝑟𝑛𝑒 ≥ 𝛼1 ≥ 𝛼2 ≥ 𝛼3 ≥ … ≥ 𝛼𝑛 ≥ 𝑟𝑛𝑒

In Eq. (4), 𝑟𝑛𝑒 is defined as the radius at the nozzle exit which occurs immediately before the extrudate swell, and 𝑆(𝛼) is the free surface defined by the 𝛼 values. Here, the 𝛼 values serve as design variables in the minimization and directly define the expansion of the polymer extrudate. The normal stress 𝜎𝑛 and tangential stress 𝜎𝑡 components on the free surface are computed from 𝜕𝑣𝑟 2 𝜕𝑣𝑟 𝜕𝑣𝑧 ) 𝑛𝑟 + 𝜇 ( + )𝑛 𝑛 𝜕𝑟 𝜕𝑧 𝜕𝑟 𝑧 𝑟

(5)

𝜕𝑣𝑧 𝜕𝑣𝑧 𝜕𝑣𝑟 𝜎𝑡 = −𝑝𝑛𝑧2 + 2𝜇 ( ) 𝑛𝑧2 + 𝜇 ( + )𝑛 𝑛 𝜕𝑧 𝜕𝑟 𝜕𝑧 𝑟 𝑧

(6)

𝜎𝑛 = −𝑝𝑛𝑟2 + 2𝜇 (

When the integral of the sum of the normal and tangential stress components along the free surface as defined in Eq. (4), reaches zero the correct extrudate swell shape is obtained. The final geometry from the minimization defines the entire flow field for a given FDM nozzle and prescribed processing conditions. 2.2 Fiber Orientation Models One goal of this work is to calculate the orientation of fibers within the polymer melt during the FDM polymer extusion process. Here, we calculate fiber orientation along streamlines throughout the entire flow domain of an FDM nozzle including the extrudate swell region. The evaluation of fiber orientation was first defined by Jeffery [25] who derived the periodic tumbling motion of a single ellipsoid in a viscous, incompressible simple shear flow. The motion of a unit vector 𝐩 along the axis of an ellipsoid may be evaluated from Jeffery’s equation as [3] 𝐷𝐩 = 𝜴 ⋅ 𝐩 + 𝜆(𝜞 ∙ 𝐩 − 𝜞: 𝐩𝐩𝐩) 𝐷𝑡

(7)

where 𝜴 is the vorticity tensor defined as 1 [(∇𝐯) − (∇𝐯)T ] 2 and 𝜞 is the rate of deformation tensor given as

(8)

1 [(∇𝐯) + (∇𝐯)T ] 2 In the above, 𝜆 is a coefficient related to fiber geometry as

(9)

𝜴=

𝜞=

λ=

𝐴𝑅 ∗ 2 − 1

(10) 𝐴𝑅 ∗ 2 − 1 where 𝐴𝑅 ∗ is the equivalent aspect ratio which is calculated for a cylindrical fiber as (see e.g., Zhang and Smith [27]) (11) 𝐴𝑅 ∗ = 0.000035𝐴𝑅 3 − 0.00467𝐴𝑅 2 + 0.764𝐴𝑅 + 0.404 here 𝐴𝑅 = 𝑙/𝑑 is the geometric aspect ratio of the cylindrical fiber with length 𝑙 and diameter 𝑑. Jeffery’s equation was extended by Folgar and Tucker [22] to represent multiple fibers in a fluid flow using the probability density function for fiber orientation 𝜓(𝐩). The rate of change in 𝜓(𝐩) is given as [22] 𝑑𝜓 (12) = −∇ ∙ 𝐩 ∙ (( 𝜴 ∙ 𝐩 + λ𝜞 ∙ 𝐩 − λ𝜞: 𝐩𝐩𝐩)𝜓(𝐩)) + 𝑫𝒓 𝑑𝑡 Unfortunately, this form is computationally prohibitive for flows such as those encountered in an FDM nozzle, and it is computationally more efficient to employ the orientation tensor approach first proposed by Advani and Tucker [3]. The fiber orientation tensor evolution equation may be written as [3]

𝐷𝑨 (13) = −(𝜴 ∙ 𝑨 − 𝑨 ∙ 𝜴) + 𝜆(𝜞 ∙ 𝑨 + 𝑨 ∙ 𝜞 − 2𝔸: 𝜞) + 𝑫𝒓 𝐷𝑡 where 𝑨 and 𝔸 are the second and fourth order orientation tensors, respectively, given as 𝐴𝒊𝒋 = ∮ 𝑝𝑖 𝑝𝑗 𝜓(𝒑) 𝑑𝒑

(14)

𝔸𝑖𝑗𝑘𝑙 = ∮ 𝑝𝑖 𝑝𝑗 𝑝𝑘 𝑝𝑙 𝜓(𝒑) 𝑑𝒑

(15)

𝑆

𝑆

Here 𝑆 is the surface of the unit sphere. Eqs. (12) and (13) include the fiber orientation diffusion term, 𝑫𝒓 , which accounts for fiber-fiber interactions. Several fiber orientation diffusion functions have been proposed including the IRD-RSC by Wang et al. [29], the ARD-RSC by Phelps and Tucker [26], and the IRD by Folgar-Tucker [22]. Our calculations employ the commonly used Folgar-Tucker isotropic rotary diffusion (IRD) function which is written as 𝑫𝒓 = 2𝐶𝐼 𝛾̇ (𝑰 − 3𝑨)

(16)

where 𝛾̇ is the scalar magnitude of the rate of deformation tensor (cf. Eq. 9) defined as 1

(17)

𝛾̇ = (𝜞: 𝜞)2

and 𝐶𝐼 is the empirically derived interaction coefficient. There exist several methods to predict 𝐶𝐼 and we will use the form presented by Bay [30] 𝑉𝐹∗𝑙 𝑑

𝐶𝐼 = 0.018408𝑒 −0.714807

(18)

where 𝑉𝐹 is the volume fraction of fibers in the polymer suspension and 𝑙/𝑑 = 𝐴𝑅 which is the cylindrical fiber aspect ratio for fibers. Observe in Eq. (13) the appearance of the fourth-order orientation tensor, 𝔸 which must be known to solve Eq. 13. Unfortunately, the equation of motion for 𝔸 contains the sixth-order orientation tensor (not shown), and every even ordered equation of motion contains the next higher even ordered orientation tensor. This necessitates the need for a closure whereby a higher ordered tensor is approximated as a function of lower ordered tensors. There exist many closure methods such as the Fast Exact Closure by Montgomery-Smith, et al. [28] which is a computationally efficient version of the exact closure by Montgomery-Smith, et al. [31], the Invariant Based Orthotropic Fitted Closure by Chung and Kwon [32], and the Orthotropic Fitted Closure (ORT) by Cintra and Tucker [23], which was improved by Wetzel [24]. The ORT is used in the current study due to its computational efficiency and stability. The ORT closure computes 𝔸 from the eigenvalues of the second order orientation tensor 𝑨 using 15 independent coefficients, as given by Verweyst [33], to obtain the approximate fourth order orientation tensor 2 1 2 3 4 5 2 6 2 7 2 8 ̃ 𝑂𝑅𝑇 𝔸 𝑚𝑚 = 𝐶𝑚 + 𝐶𝑚 𝜆1 + 𝐶𝑚 𝜆2 + 𝐶𝑚 𝜆1 𝜆2 + 𝐶𝑚 𝜆1 + 𝐶𝑚 𝜆1 + 𝐶𝑚 𝜆1 𝜆2 + 𝐶𝑚 𝜆1 𝜆2 9 3 10 3 11 2 2 12 3 13 14 4 15 4 + 𝐶𝑚 𝜆1 + 𝐶𝑚 𝜆2 + 𝐶𝑚 𝜆1 𝜆2 + 𝐶𝑚 𝜆1 𝜆2 + 𝐶𝑚 𝜆1 𝜆32 +𝐶𝑚 𝜆1 + 𝐶𝑚 𝜆2

(19)

1 15 where 𝜆1 and 𝜆2 are the eigenvalues of the second-order fiber orientation tensor, and 𝐶𝑚 -𝐶𝑚 are fitted coefficients from different flow types.

Our work quantifies fiber alignment using components of 𝑨 where orientation diffusion is modeled with Eq. (16) and the fourth order orientation tensor closure approximation is given by Eq. (19). It is worth noting that an important property of 𝑨 is that its trace equates to unity. It follows that a fiber orientation state for a uniformly random yields a tensor 𝑨 with components equal to zero except the main diagonal which has values of 𝑨𝑖𝑗 = 1/3, 𝑖 = 1,2,3 (no summation implied). Similary, an orientation tensor 𝑨 for fibers perfectly aligned in the 𝑥3 direction would yield all components having a value of zero except 𝐴33 = 1. 2.3 Calculation of Mechanical Properties from Fiber orientation The calculated fiber orientation state may be used to evaluate mechanical properties of a fiber filled polymer. Several models have been proposed to predict the elastic properties of a unidirectional composite as summarized by Tucker and Liang [34]. These approaches assume that the fiber and the polymer matrix are linearly elastic, the matrix is isotropic, the fibers are either isotropic or transversely isotropic, the fibers are axisymmetric and of the same size which is defined by the fiber’s geometric aspect ratio and volume fraction, and the fibers are well bonded to the matrix and remain that way during deformation. The initial work done by Eshelby [35] evaluated the elastic field around a single ellipsoidal inclusion. This was followed by Mori and Tanaka [36] who considered dilute suspensions, which was later used by Tandon and Weng [19] who obtained expressions for the elastic constants of unidirectional short fiber composites. The use of the orientation homogenization method allows for the mean stiffness of a composite with known fiber orientation. A representative volume of the composite is divided into a set of smaller aggregates with unidirectional fiber alignment where the stiffness of each aggregate is evaluated using a unidirectional stiffness model. The aggregates are homogenized using a constant mean strain (Voigt) or constant mean stress (Reuss) across all aggregates [37]. The stiffness of the representative volume is simply the combined stress of the aggregates. The orientation homogenization method introduced by Advani and Tucker [3] and derived in [21] is used in the current study to calculate the elastic constants of the short fiber composite and is written as 〈𝐶𝑖𝑗𝑘𝑙 〉 = 𝑏1 (𝔸𝑖𝑗𝑘𝑙 ) + 𝑏2 (𝐴𝑖𝑗 𝛿𝑘𝑙 + 𝐴𝑘𝑙 𝛿𝑖𝑗 ) + 𝑏3 (𝐴𝑖𝑘 𝛿𝑗𝑙 + 𝐴𝑖𝑙 𝛿𝑗𝑘 + 𝐴𝑗𝑙 𝛿𝑖𝑘 + 𝐴𝑗𝑘 𝛿𝑖𝑙 ) + 𝑏4 (𝛿𝑖𝑗 𝛿𝑘𝑙 ) + 𝑏5 (𝛿𝑖𝑘 𝛿𝑗𝑙 + 𝛿𝑖𝑙 𝛿𝑗𝑘 )

(23)

where 𝑏1 = 𝐶11 + 𝐶22 − 2𝐶12 − 4𝐶66

𝑏4 = 𝐶23

𝑏2 = 𝐶12 + 𝐶22

𝑏5 =

𝑏3 = 𝐶66 +

1 2(𝐶23 + 𝐶22 )

1 2(𝐶22 − 𝐶23 )

(24)

In the above, C11, C12, C22, C23, and C66 are the elasticity tensor coefficients of the associated unidirectional short fiber composite. In this study, these coefficients are calculated using the Tandon-Weng [19] approach.

3 Results 3.1 FDM Nozzle Description The cross section of a common FDM extrusion nozzle, such as that found in a Makerbot 2X (Makerbot Industries, LLC, Brooklyn, NY), is shown in Figure 2. The dimensions shown in the given cross section are used to create the polymer melt flow domain for our FDM nozzle fiber orientation simulations. An axisymmetric slice is taken from the cross section appearing in Figure 2 and is used to define the flow domain in the COMSOL Multiphysics finite element software (Comsol, Inc., Burlington, MA). The axisymmetric fluid domain built in COMSOL is shown in Figure 3 along with the associated boundary conditions. The flow enters the nozzle at z = 5.625 mm with a laminar flow having parabolic inlet velocity 𝑣𝑧𝑖𝑛𝑙𝑒𝑡 profile defined as 𝑣𝑧𝑖𝑛𝑙𝑒𝑡 = −𝑣𝑚𝑎𝑥 [1 − (

𝑟𝑖 𝑟𝑚𝑎𝑥

2

) ]

(25)

with 𝑣𝑟𝑖𝑛𝑙𝑒𝑡 = 𝑣𝜃𝑖𝑛𝑙𝑒𝑡 = 0. In Eq. (25), 𝑣𝑚𝑎𝑥 is the maximum velocity at the inlet of the fluid domain, 𝑟𝑖 is the radial location, and 𝑟𝑚𝑎𝑥 is the maximum radius of the fluid domain at the inlet. For this study we use 𝑣𝑚𝑎𝑥 = 4 𝑚𝑚/𝑠 and 𝑟𝑚𝑎𝑥 = 0.875 𝑚𝑚. This inlet velocity yields an average nozzle exit velocity of 65 𝑚𝑚/𝑠 which is within the range of extrusion speeds for most commercially available desktop FDM printers. The outlet at 𝑧 = 0 𝑚𝑚 is at the reference pressure 𝑃 = 0 𝑃𝑎 . The interior wall of the FDM nozzle appears as the right boundary of the fluid domain in Figure 3, and is prescribed as a noslip wall boundary with 𝑣𝑟 = 𝑣𝑧 = 0 for 𝑧 ≥ 0.875 𝑚𝑚. A slip wall condition is prescribed along the free surface defining the extrudate swell boundary for 𝑧 < 0.875 𝑚𝑚, and is given as 𝐧 ∙ 𝐯 = 0 where 𝐧 is the surface normal and 𝐯 is the velocity vector. The entire left boundary seen in Figure 3 of the model is the axis of symmetry which is assigned the axisymmetric boundary. An incompressible Newtonian fluid is assumed in our simulations having a density of 𝜌 = 1040 𝑘𝑔/𝑚3 and dynamic viscosity of 𝜇 = 350 𝑃𝑎 ∙ 𝑠, which is typical for an ABS plastic at 230 𝑑𝑒𝑔𝑟𝑒𝑒𝑠 𝐶 under a shear rate of 575 𝑠 −1. It is understood that ABS, is a shear-thinning non-Newtonian fluid, but assuming an isothermal Newtonian fluid simplifies the simulation process and allows for extrudate swell results to be compared with known data from Georgiou [13] and Reddy and Tanner [14]. The finite element model for our simulations is defined using the boundary conditions appearing in Figure 3 and described above. The axisymmetric fluid flow problem is evaluated in Comsol Multiphysics using LiveLink which employs Matlab (The Mathworks, Inc., Natick, MA). The objective function in Eq. (4) is minimized using the optimization approach described in Heller [38]. This approach yields a radial extrudate swell of approximately 13% for the FDM nozzle flow which agrees well with values obtained by Georgiou, et al. [13] and Reddy and Tanner [14].

3.2 Fiber Orientation Calculation The velocity and velocity gradients of the fluid within the entire flow domain are calculated along the ten streamlines appearing in Figure 4. The streamlines used here are positioned at equally spaced locations across the outlet of the fluid domain and are defined in COMSOL with a numerical integration between the flow exit and entrance. At every point along each streamline, values of the velocity and velocity gradients (𝑣𝑟 , 𝑣𝑧 ,

𝑑𝑣𝑟 𝑑𝑟

,

𝑑𝑣𝑟 𝑑𝑧

,

𝑣𝑟 𝑟

,

𝑑𝑣𝑧 𝑑𝑟

,

𝑑𝑣𝑧 𝑑𝑧

) are obtained for the

evaluation of fiber orientation. We assume a cylindrical fiber aspect ratio 𝐴𝑅 = 15 and a volume fraction VF = 15% in all fiber orientation calculations. The equivalent aspect ratio, 𝐴𝑅 ∗ is calculated using Eq. (11) proposed by Zhang [27] and returns a λ value of 0.9834. This yields an interaction coefficient 𝐶𝐼 = 0.0057 based on Eq. (18) as given by Bay [30]. In a related study we have investigated the pure shear steady state values of the second order fiber orientation tensor, 𝑨, for the given 𝐶𝐼 and 𝜆 values to provide a more accurate inlet condition and to shorten the required length of the upstream region (cf. Heller [38]). The pure shear steady state inlet condition is used for all simulations and the associated plots given here. The results of the fiber orientation calculations within the FDM nozzle flow including the extrudate swell region calculated as described above are shown in Figures 4 and 6. For this study the flow domain is extended 5R past the nozzle exit to capture the effects of the extrudate swell, where R is the radius of the tubular section between CZE and NE of the FDM nozzle in Figure 4. In Figure 4 (a), CZS represents the start of the convergence zone, CZE represents the end of the convergence zone, and NE represents the nozzle exit. Figure 4 (b) shows the A33 component for the 10 streamlines evaluated in the fluid domain. The A33 component represents the fiber alignment in the direction of the flow through the nozzle and across the width of the nozzle. Figure 4 (b) shows that the A33 component of the orientation tensor increases in the convergence zone which spans from CZS to CZE. Note that A33 reaches a value near unity, which represents a nearly uniaxial alignment in the z-direction. The maximum value of A33 at its peak is A33 = 0.972 along streamline 1. This increase in fiber orientation is due to elongation flow which is defined by an increase in

𝑑𝑣𝑧 𝑑𝑧

as seen in Figure 5 (d). The velocity gradient

𝑑𝑣𝑧 𝑑𝑧

begins to increase before the convergence zone which causes a similar increase in the A33 component slightly before the convergence zone. Immediately after the end of the convergence zone, in the straight portion of the nozzle from CZE to NE, there is a decrease in the A33 component which can be related to the velocity gradients shown in Figure 5 (a & d). The velocity gradients

𝑑𝑣𝑧 𝑑𝑧

and

𝑑𝑣𝑟 𝑑𝑟

return to zero in this region which results in a shear dominated flow. The

high shear causes the fiber orientation to begin to return to pure shear steady state orientation. After the nozzle exit NE there is a decrease in the A33 component of the fiber orientation which is due to the negative elongation component

𝑑𝑣𝑧 𝑑𝑧

and an expansion flow

𝑑𝑣𝑟 𝑑𝑟

shown in Figure 5 (a &

d), both of which decrease the fiber alignment in the z direction. Figure 4 (b) also shows that fibers along streamlines 9 and 10 undergo decreasing alignment over a short distance in the convergence zone of the nozzle. This is due to an increase in the expansion flow, as seen with increasing faster than the axial elongation velocity gradient

𝑑𝑣𝑧 𝑑𝑧

𝑑𝑣𝑟 𝑑𝑟

. The magnitude of the elongation

velocity gradient quickly surpasses that for radial expansion and the A33 component increases similar to streamlines 1-8 until CZE. Figure 4 (b) shows the A33 component of the orientation tensor throughout the fluid domain along all ten streamlines. It can be seen from Figure 4 (b) that the velocity gradient effects on the fiber orientation differ in magnitude depending on the spatial relation to the orientation changing flow. Because the magnitude of fluid velocity is greater at the center of the nozzle the elongation flow causes the fibers along streamline 1 to align more quickly in the axial direction than the fibers at streamline 10. This same effect causes the decrease of fiber alignment along the axis of flow at the nozzle exit to be greater than other streamlines further from the middle of the nozzle. Figure 5 also illustrates the velocity gradients along streamline 1. These velocity gradients are shown here since they directly influence changes in fiber orientation along streamline 1 throughout the nozzle. The velocities (𝑣𝑟 , 𝑣𝑧 ) as well as 𝑣 the velocity gradient ( 𝑟 ) along streamline 1 are not included in Figure 5 as these components 𝑟

have been found to have very little or no effect on the fiber orientation. Figure 6 shows values of A33 at steady state at a cross section of the flow beyond the nozzle exit for the ten streamlines. Values of A33 that appear in Figure 6 provide a good indication of fiber alignment in the extruded composite. Figure 6 also shows that the fiber orientation is somewhat highly aligned for all values of radial position r, all showing A33 greater than 0.55 across the flow. Note that there is a high degree of alignment near the wall of the extrusion nozzle at streamline 10 and high alignment at the center of the flow at streamline 1. There is also an obvious reduction in alignment between the wall and the center of the flow. The orientation of the fibers leaving the nozzle are highly aligned even considering the decrease in orientation right after the nozzle exit. This is due to the fact that the convergence zone has greatly increased the fiber alignment and that the reduction in alignment in the extrudate swell region cannot completely overcome the high fiber alignment. 3.3 Parametric Study of FDM Nozzle Geometry A parametric study was performed to evaluate the effect of nozzle geometry on fiber orientation. Four parameters were selected for the study which define the nozzle geometry as shown in Figure 7. Parameter A is the convergence zone length, B is the straight tube length, C is the amount of radial nozzle expansion, and D is the length of the nozzle expansion. In all cases, parameters are defined as multiples of the nozzle tube radius, R. The average A33 computed over a cross section of the free stream flow beyond the nozzle exit, denoted here as 𝐴33 is used as a singular metric to describe the effect of the nozzle geometry parameters on fiber alignment. Note that the 𝑥3 direction is the longitudinal flow direction along the main axis of the FDM nozzle, therefore, the A33 component is a measure of fiber alignment along the axial direction of the extrudate. 𝐴33 is calculated using Simpson’s 3/8 rule from the 10 values of A33 such as those appearing in Figure 6 calculated on streamlines defined above. The parametric study is performed in two steps: first, thirty-six possible individual combinations are run for parameters A and B, second, the minimum fiber alignment case for flows defined by values of A and B is then used for the computation of each of the thirty-six possible combinations of C and D. All parameters used in the study are listed in Table 1. From

the parametric study the maximum fiber alignment case and the minimum fiber alignment case are of the greatest interest and will be discussed further.

Figure 8 illustrates the effect of changing the nozzle geometry on the average fiber alignment by plotting values of 𝐴33 for the thirty-six combinations of the A and B geometry parameters described above. It is shown that the maximum fiber alignment state occurs at (A = 5R, B = 0.1R) and the minimum fiber alignment state is at (A = 5R, B = 10R). Figure 8 also shows that parameter A, the length of the convergence zone, has very little effect on the degree of fiber alignment. The fiber alignment, as measured by values of 𝐴̅33 , varies by only 0.1 over a range of 10R radii for the convergence zone length. Figure 8 also shows that parameter B, the length of the straight portion of the nozzle, has a much larger effect on fiber alignment than parameter A. As parameter B increases fiber alignment substantially decreases from 𝐴33 = 0.886 to 𝐴33 = 0.630 over B values from 0.1R to 10R. This decrease in fiber alignment is due to the pure shear flow in the straight portion of the nozzle which returns the fiber alignment to the pure shear steady state orientation yielding a lower fiber alignment than that found immediately after the convergence zone ends. Figure 9 shows the A33 component for all 10 streamlines in the fluid domain defined by (A = 5R, B = 0.1R, C = 1R, and D = 0R) and illustrates the amount of alignment as fluid travels through the length of the nozzle.

Note that A33 values along streamlines 1-5 are very similar to those in Figure 4 in that fiber alignment increases slightly before CZS and then continues to increase through the convergence zone to CZE. Fibers along streamlines 6-10, however, show a decrease in alignment from shortly before CZS to slightly after CZS, then increase quickly through the convergence zone. The larger decrease in fiber alignment and the inclusion of more streamlines showing this effect is due to the lower slope of the convergence zone wall than in the previous nozzle calculation. The reduced slope causes the magnitude of remains higher than

𝑑𝑣𝑧 𝑑𝑧

𝑑𝑣𝑟 𝑑𝑟

to be greater than seen in previous calculations, which

farther into the interior of the nozzle. The straight portion of the tube in

this nozzle design is only 0.1R, therefore the decrease in fiber alignment from CZE to NE is negligible. The expansion flow (as indicated by high values of indicated by high values of

𝑑𝑣𝑧 𝑑𝑧

𝑑𝑣𝑟 𝑑𝑟

), and contraction flow (as

), in the extrudate swell at the nozzle exit decreases fiber

alignment as seen above. However, due to the initial fiber alignment being very high, the overall fiber alignment remains high after the nozzle exit. Figure 10 shows values of 𝐴̅33 for thirty-six possible combinations of the A, B, C, and D geometry parameters. From this figure, the minimum fiber alignment with nozzle expansion case

(A = 5R, B = 10R, C = 1.15R, and D = 0.5R) is selected for further study to show the effects of an expanding nozzle near its exit. Figure 10 also shows that parameter C, the amount of radial expansion of the nozzle, decreases the fiber alignment in the nozzle. Parameter D, the length of the nozzle expansion, is shown to have little effect on the fiber alignment due to the small fiber alignment change over the span of expansion lengths. Parameter C shows a more substantial decrease in fiber alignment but it is relatively small in comparison to parameter B. Larger values of expansion may cause a greater decrease in the fiber alignment, however, a maximum nozzle radial expansion of C = 1.15R is used here to ensure that the flow does not separate from the wall inside the nozzle. Figure 11 shows the A33 component for all 10 calculated streamlines in the fluid domain defined by (A = 5R, B = 10R, C = 1.15R, and D = 0.5R) which quantify the degree of fiber alignment throughout the length of the nozzle and over the width of the nozzle. In Figure 11, STE defines the Straight Tube End of the nozzle. Due to the expansion zone the end of the straight tube is no longer the nozzle exit, NE.

Fiber orientation along streamlines 1-10 in Figure 11 show trends similar to those in Figure 4 (b) from inlet to NE due to similarities in nozzle geometry. From STE to NE a decrease in fiber alignment is seen for all streamlines which is due to expansion flow in the nozzle expansion zone. After NE there is a further decrease in fiber alignment due to radial expansion flow and axial contraction flow in the extrudate swell region. Results appearing in both Figures 9 and 11 show that changes in nozzle shape can have a significant effect on fiber alignment within the extruded polymer. These changes in geometry are expected to have a significant effect on the deposited FDM bead just beneath the nozzle exit. 3.4 Mechanical Properties from Fiber Orientation Data Mechanical properties of the fiber filled polymer extrudate are calculated from the second and fourth order orientation tensor as given in Eqs. (23) and (24). Of particular interest are the composite’s moduli of elasticity and Poisson’s ratio. We assume that fiber and polymer matrix in the FDM composite are isotropic having a carbon fiber modulus 𝐸𝑓 = 240 GPa and Poisson’s ratio 𝑣𝑓 = 0.2. The modulus of elasticity and the Poisson’s ratio of the matrix are 𝐸𝑚 = 2 GPa and 𝑣𝑚 = 0.4 respectively. In addition, we assume a volume fraction, 𝑉𝐹 = 15% and the aspect ratio, 𝐴𝑅 = 15 as stated earlier for the calculation of λ and CI. The mechanical properties of the fiber filled polymer in the FDM nozzle are evaluated at two locations in the fluid domain. The first location is inside the tubular section of the nozzle slightly before the nozzle exit at z = 0.9 mm and the second is in the steady state region of the extrudate swell after the nozzle exit at z = 0.05 mm. This allows for the evaluation of the effect that the extrudate swell has on the printed mechanical properties. Figures 12 (a) and (b) show the radial and axial modulus of elasticity (Err and Ezz, respectively) with and without the consideration of the extrudate swell phenomena. The modulus is evaluated from the center of the polymer extrudate outward along the radius of the axisymmetric fluid

domain. Figures 12 (c) and (d) show the Poisson’s ratio and shear modulus of the extrudate. The swell and non-swell values occur at different radial values due to the 13% swell increase. Figure 12 shows that there is a noticeable difference between the moduli of the swell and nonswell polymer composites, therefore the inclusion of extrudate swell in the model is necessary due to its effect on the mechanical properties. The shear modulus and Poisson’s ratio show a substantial change as well between the In Nozzle and In Swell results. The decrease in the fiber alignment due to extrudate swell causes the extruded polymer composite to be more compliant in the axial direction and to increase the radial modulus due to higher alignment in the transverse direction. The mechanical properties for the extruded polymer composite are given in Table 2. A decrease in the axial modulus of elasticity of 19.9% occurs when the extrudate swell is included in the model. An increase in the Poisson’s ratio of 27.3% is realized as well when the extrudate swell is included. The mechanical properties of the fiber filled polymer were then calculated for the maximum and minimum fiber alignment states found in the parametric study. This allows for the limits of the possible spectrum from the parametric study to be seen. In this study the mechanical properties are only calculated at the steady state region after the nozzle exit and extrudate swell. The effect of the swell has already been determined which makes the inner nozzle value unnecessary. Figure 13 shows the two extremes we modeled for degree of fiber alignment. These two extremes were used to determine how greatly the material stiffness could be affected by changing the nozzle geometry and thus the resulting fiber orientation. Figure 13 shows that there is a substantial increase in the axial modulus and the correlated decrease in the radial modulus for the extruded fiber filled polymer of the maximum 𝐴̅33 fiber alignment state. From this result it is understood that the nozzle geometry can be changed to fit a preferred fiber alignment state that allows for specific mechanical properties to be had in a completed part. The mechanical properties for the extruded composite for the maximum and minimum fiber alignment state are given in Table 3.

Table 3. shows that by changing the nozzle geometry an increase in the axial modulus of elasticity of 56.1% can be realized. This substantial increase shows that an FDM nozzle could be designed and modified to produce specific material properties of a printed composite material.

5 Conclusions It is shown in this study that the convergence zone and extrudate swell have a large effect on the fiber orientation state and mechanical properties for short fiber suspensions within the polymer melt during FDM processing. The effect of the extrudate swell free surface on the fiber orientation was seen to be substantial. A decrease in the axial modulus of elasticity of 19.9% and a related increase in the radial modulus of elasticity of 6.62% is seen for the given FDM nozzle geometry. Since the extrudate swell is the last source of major fiber orientation change in the polymer melt during extrusion, it has a significant effect on the final orientation of the fiber

suspension and the resulting mechanical properties of the FDM bead. Parametric studies were conducted to further explore and quantify the effect of nozzle geometry on the resulting fiber alignment state and the resulting mechanical properties. It was shown that a change in the nozzle convergence length (Parameter A) and the nozzle expansion length (Parameter D) had very little effect on the fiber alignment in the nozzle flow. The fiber alignment was seen to undergo more substantial decreases by changing the length of the straight tube portion (Parameter B) and the amount of nozzle radial expansion (Parameter C). This result allows for the augmentation of nozzle geometry to provide a preferred fiber alignment state in the extruded polymer composite. It was shown that the elastic moduli could be increased by 56.1% in the axial direction using the parameters given in Figure 7. This substantial change in the elastic moduli shows that the range of possible fiber orientation states and resulting mechanical properties is quite large. These results suggest that a large range of mechanical properties are possible by changing the FDM nozzle geometry. This study shows that the extrudate swell has a large effect on the fiber orientation and the resulting mechanical properties. It also shows that the FDM nozzle geometry can be modified to obtain a desired fiber alignment state and mechanical properties. Further development that includes viscoelastic fluids, coupled fluid flow with fiber orientation models, deposition onto a moving platform, and transient effects will help to make this modeling approach more useful in understanding FDM processing of fiber-filled polymer composites.

References [1] L.J. Love, V. Kunc, O. Rios, C.E. Duty, A.M. Elliot, B.K. Post, R.J. Smith, C.A. Blue, The Importance of Carbon Fiber to Polymer Additive Manufacturing, J. Mater. Res. 29 (17) (2014) 1893-1898. [2] J. Nixon, B. Dryer, D. Chiu, I. Lempert, D.I. Bigio, Three Parameter Analysis of Fiber Orientation In Fused Deposition Modeling Geometries, (2014) 985–995. [3] S.G. Advani, C.L. Tucker III, The Use of Tensors to Describe and Predict Fiber Orientation in Short Fiber Composites, J. Rheol. 31 (1987) 751-784. [4] C.L. Tucker and S.G. Advani, Processing Short-Fiber Systems, in Flow and Rheology in Polymer Composites Manufacturing, Ed. S.G. Advani, Elsevier Science (1994). [5] S.J. Lee and S.J. Lee, Numerical Stimulation of Extrudate Swell of Semi-concentrated Fiber Suspensions, in Recent Advances in Non Newtonian Flows, AMD Vol. 153 (ASME, 1992) 87. [6] B.E. VerWeyst, C.L. Tucker III, Fiber Suspensions in Complex Geometries: FlowOrientation Coupling, Can. J. Chem. Eng. 80 (6) (2002) 1093–1106. [7] A. Baloch, M.F. Webster, A Computer Simulation of Complex Flows of Fibre Suspensions, Computers and Fluids, 24 (2) (1995) 135-151. [8] G. Ausias, J.F. Agassant, M. Vincent, Flow and Fiber Orientation Calculations in Reinforced Thermoplastic Extruded Tubes, Intern. Polymer Processing 9 (1994) 51-59. [9] H. Uematsu, N. Horisawa, T. Horikida, S. Tanoue, Y. Iemoto, Effect of Carbon Fiber on the Capillary Extrusion Behaviors of High-Density Polyethylene, Polym. J. 45 (2013) 449-456. [10] A. Wagner, R. Yazici, D.M. Kalyon, Extrudate Swell Behavior of Glass Fiber Filled Polyamide 6, J. Polym. Composite. 17 (6) (1996) 840-849. [11] J. Stabik, Influence of Filler Particle Geometry on Die Swell, Intern. Polymer Processing 19 (2004) (4) 350-355. [12] M.L. Becraft and A.B. Metzner, The Rheology, Fiber Orientation, and Processing Behavior of Fiber-filled Fluids, J. Rheol. 36 (1) (1992) 143-174. [13] G.C. Georgiou, The Compressible Newtonian Extrudate Swell Problem, Int. J. Numer. Meth. Fl. 20 (1995) 255-261. [14] K.R. Reddy, R.I. Tanner, Finite Element Solution of Viscous Jet Flows with Surface Tension, Computers and Fluids, 6 (2) (1978) 83-91. [15] E. Mitsoulis, G.C. Georgiou, Z. Kountouriotis, A Study of Various Factors Affecting Newtonian Extrudate Swell, Computers and Fluids 57 (2012) 195-207. [16] E. Mitsoulis and S.G. Hatzikiriakos, Annular Extrudate Swell of Fluoropolymer Melt Intern. Polymer Processing 27 (5) (2012) 535-546.

[17] M. Kanvisi, S. Motahari, B. Kaffashi, Numerical Investigation and Experimental Observation of Extrudate Swell for Viscoelastic Polymer Melts, Intern. Polymer Processing 29 (2) (2014) 227-232. [18] J.C. Halpin, J.L. Kardos, The Halpin-Tsai Equations: A Review, Polym. Eng. Sci. 16 (5) (1976) 344–352. [19] G.P. Tandon, G.J. Weng, The Effect of Aspect Ratio of Inclusions on the Elastic Properties of Unidirectionally Aligned Composites, Polym. Composite. 5 (4) (1984) 327-333. [20] C.W. Camacho, C.L. Tucker III, S. Yalvac, R.L. McGee, Stiffness and Thermal Expansion Predictions for Hybird Short Fiber Composites, Polym. Composite. 11 (4) (1990) 229–239. [21] D.A. Jack, D.E. Smith, The Effect of Fiber Orientation Closure Approximations on Mechanical Property Predictions, Compos. Part A-Appl. S. 38 (3) (2006) 975-982. [22] F. Folgar, C.L. Tucker III, Orientatrion Behavior of Fibers in Concentrated Suspensions, J. Reinf. Plast. Comp. 3 (2) (1984) 98-119. [23] J.S. Cintra, C.L. Tucker III, Orthotropic Closure Approximations for Flow-Induced Fiber Orientation, J. Rheol. 39 (6) (1995) 1095–1122. [24] E.D. Wetzel, C.L. Tucker III, Area Tensors for Modeling Microstructure During Laminar Liquid-Liquid Mixing, Int. J. Multiphas. Flow. 25 (1999) 35-61. [25] G.B. Jeffery, The Motion of Ellipsoidal Particles Immersed in a Viscous Fluid, P. R. Soc. Lond. A-Conta. 102 (715) (1923) 161–179. [26] J. H. Phelps, C.L. Tucker III, An Anisotropic Rotary Diffusion Model for Fiber Orientation in Short and Long Fiber Thermoplastics, J. Non-Newton. Fluid. 156 (3) (2009) 165-176. [27] D. Zhang, Flow-Induced Micro- and Nano-Fiber Suspensions in Short-Fiber Reinforced Composite Materials Processing, PhD thesis, University of Missouri, 2013. [28] S. Montgomery-Smith, D.A. Jack, D.E. Smith, The Fast Exact Closure for Jeffery’s Equation with Diffusion, J. Non-Newton. Fluid. 166 (7-8) (2011) 343-353. [29] J. Wang, J.F. O’Gara, C.L. Tucker III, An Objective Model for Slow Orientation Kinetics in Concentrated Fiber Suspensions: Theory and Rheological Evidence, J. Rheol. 52 (5) (2008) 1179-1200. [30] R.S. Bay, Fiber Orientation in Injection Molded Composites: A Comparison of Theory and Experiment. PhD Thesis, University of Illinois at Urbana-Champaign, 1991. [31] S. Montgomery-Smith, W. He, D.A. Jack, D.E. Smith, Exact Tensor Closures for the Three Dimensional Jeffery’s Equation, J. Fluid. Mech. 680 (2011) 321-335.

[32] S.T. Chung, T.H. Kwon, Coupled Analysis of Injection Molding Filling and Fiber Orientation, Including In-Plane Velocity Gradient Effect, Polym. Composite. 17 (6) (1996) 859– 872. [33] B.E. VerWeyst, Numerical Predictions of Flow Induced Fiber Orientation in ThreeDimensional Geometries. PhD thesis, University of Illinois at Urbana Champaign, 1998. [34] C.L. Tucker III, E. Liang, Stiffness Predictions for Unidirectional Short-Fiber Composites: Review and Evaluation, Compos. Sci. Technol. 59 (5) (1999) 655–671. [35] J.D. Eshelby, The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems, P. R. Soc. Lond. A. Mat. 241 (1226) (1957) 376–396. [36] T. Mori, K. Tanaka, Average Stress in Matrix and Average Elastic Energy of Materials with Misfitting Inclusions, Acta Metall. Mater. 21 (5) (1973) 571–574. [37] R.W. Hill, The Elastic Behaviour of a Crystalline Aggregate, P. R. Soc. Lond. A. 65 (5) (1952) 349–354. [38] B.P. Heller, Effects of Nozzle Geometry and Extrudate Swell on Fiber Orientation in Fused Deposition Modeling Nozzle Flow, Master’s Thesis, Baylor University, 2015.

Figure 1. Extrudate swell boundary conditions and expansion value representation

Figure 2. Cross section of FDM extrusion nozzle, Dimensions are in millimeters (mm)

Figure 3. Boundary conditions for nozzle domain

Figure 4. (a) Nozzle configuration, dimensions, and streamlines, (b) Fiber alignment throughout the FDM nozzle

Figure 5. Velocity gradients along streamline 1

Figure 6. Steady state fiber alignment values for ten streamlines after nozzle exit

Figure 7. Geometry and parameter values for parametric study

Figure 8. Average fiber alignment values for parametric study of nozzle straight portion length and convergence zone length, (C = 1R, D = 0R).

Figure 9. Fiber alignment throughout the FDM nozzle which produces the maximum average fiber alignment, (A = 5R, B = 0.1R, C = 1R, D = 0R).

Figure 10. Average fiber alignment values for parametric study of nozzle expansion and nozzle expansion length.

Figure 11. Fiber alignment throughout the FDM nozzle which produces the minimum average fiber alignment with nozzle exit expansion, (A = 5R, B = 10R, C = 1.15R, D = 0.5R)

Figure 12. Radial and axial moduli of the printed fiber filled polymer with and without the effect of extrudate swell.

Figure 13. Radial and axial moduli of the printed fiber filled polymer using the maximum and minimum fiber alignment state

Table 1. Parameter Values Parameter A B C D

Value 5R, 7R, 9R, 11R, 13R, and 15R 0.1R, 2R, 4R, 6R, 8R, and 10R 1.025R, 1.05R, 1.075R, 1.1R, 1.125R, and 1.15R 0.5R, 1R, 1.5R, 2R, 2.5R, and 3R

Table 2. Elastic Properties in nozzle and in swell immediately outside the nozzle exit Property

In Nozzle

In Swell

% Difference

Ezz (GPa)

8.96

7.34

19.9

Err (GPa) νrz Grz (GPa)

3.21 0.15 1.25

3.43 0.20 1.39

6.6 27.3 10.5

Table 3. Elastic Properties in swell for maximum and minimum fiber alignment Property

Maximum A33

Minimum A33

% Difference

Ezz (GPa)

10.6

5.95

56.1

Err (GPa) νrz Grz (GPa)

3.08 0.12 1.13

3.55 0.26 1.40

14.4 72.2 21.5