Tensile and Compression Strength Prediction and Validation in 3D-Printed Short-Fiber-Reinforced Polymers

In the current study, a methodology is validated for predicting the internal spatially varying strength properties in a single 3D-printed bead composed of 13%, by weight, carbon-fiber-filled acrylonitrile butadiene styrene. The presented method allows for the characterization of the spatially varying microstructural behavior yielding a local anisotropic stiffness and strength that can be integrated in a finite element framework for a bulk estimate of the effective stiffness and strength. The modeling framework is presented with a focus on composite structures made from large area additive manufacturing (LAAM). LAAM is an extrusion-based process yielding components on the order of meters, with a typical raster size of 10 mm. The presented modeling methods are applicable to other short-fiber-reinforced polymer processing methods as well. The results provided indicate the modeling framework yields results for the effective strength and stiffness that align with experimental characterization to within ∼1% and ∼10% for the longitudinal compressive and tensile strength, respectively, and to within ∼3% and ∼50% for the longitudinal compressive and tensile stiffness, respectively.


Introduction
The development of better manufacturing techniques promotes the development of new engineering materials, and vice versa.This concept is well embodied in the area of 3D-printed short-fiber-reinforced polymers (SFRPs).These materials offer the attractive combination of good mechanical properties compared to virgin polymers and a relatively easy processing.In addition, they help enable large area additive manufacturing (LAAM), which often uses carbon fiber to provide crucial thermal and mechanical property advantages (see, e.g., [1]).However, 3D-printed SFRPs introduce added difficulty in engineering design since their anisotropic properties are more complicated to model.Additive manufacturing allows one to customize tool paths to tailor that anisotropy to one's advantage (for example, one can orient the way a part is built so that it has greater stiffness and strength in a load-bearing direction), but good methods for being able to predict these anisotropic properties must be established so that one can generate designs with high confidence.
The anisotropic properties of 3D-printed SFRPs are a function of the internal, spatially varying fiber orientation state.Thus, predicting the properties of SFRP parts oftentimes involves the problem of predicting their internal fiber orientation state first, then predicting the properties as functions of that orientation state.Research on predicting the spatially varying fiber orientation state in a fluid flow has its roots going back about a hundred years ago to George B. Jeffery [2].The works of Charles Tucker III and his students in the past few decades have contributed greatly to the current understanding of this topic as well.Tucker and his student Suresh Advani brought orientation tensors to light as a general but very compact way to describe orientation states, and these have since become mainstream, allowing for the prediction of material properties as functions of those orientation tensors [3].Folgar and Tucker [4] introduced a fiber interaction parameter to allow for the modeling of dense suspensions, and in the past decade, multiple models have been developed to advance the initial work of Folgar and Tucker to implement fiber orientation kinematics for industrial flow processes (see, e.g., [5][6][7]), each of which uses the orientation tensor approach to characterize the fiber alignment state.This topic is discussed in detail in Section 4 while introducing the mathematical framework.Modern modeling techniques now give the ability to predict the final properties of SFRP parts based on the processing conditions used to build said parts.We wish to bring this capability to light in the additive manufacturing industry so that SFRP parts can be evaluated prior to printing, saving time and money and improving engineers' confidence level in their designs.
The next two sections of this paper discuss experimental specimen preparation and testing, respectively.A miniature LAAM system was constructed for printing specimens.Section 4 then covers some of the mathematical models used to make predictions.Finally, Section 5 compares the experimental findings with the predicted structural behavior.The novelty of the present work is the experimental validation of a methodology, also demonstrated in [8], which fully integrates fiber orientation kinetics with structural property predictions for additively manufactured composites made from the LAAM process.This work is extended to compare against experimental observations and highlights the need to properly account for both the proper fiber orientation kinetics slowness parameter and the void content within the deposited material, the latter of which is able to be neglected in other industrial processes that study fiber orientation.

Specimen Preparation 2.1. Baylor's Large Area Additive Manufacturing System
For fabricating test specimens, a miniature LAAM system was constructed at Baylor University.This system has been utilized in a variety of studies by various researchers (see e.g., [9,10]).The LAAM system has a print area of approximately 1.2 m by 1.2 m and the extruder has approximately 15 cm of vertical travel space.The bed can also be manually unfastened from the perforated tubing it is attached to, lowered, and refastened for added vertical space.The materials for the main structure of the LAAM system include cold-rolled steel railings, a 6.4 mm aluminum plate for the print bed, and 76 mm steel tubing for the frame.The LAAM system also includes a heated print bed composed of multiple zones.
The extruder on the Baylor LAAM system is a Strangpresse Extruder Model 19, which has a pellet-fed hopper at the top which was a custom-made part fabricated at Baylor University.The Strangpresse extruder shown Figure 1 is about 1 m tall with the hopper.The Strangpresse human-machine interface (HMI) allows an operator to adjust the temperatures in three zones along the length of the extruder as well as the extruder screw revolutions per minute (RPM) to control the deposition rate.The gantry system used for the motion control of the extruder is a Magnum II system from Precision Plasma LLC.The Magnum II system was made from an assembly kit for computer numerical control (CNC) plasma cutters, but its ability to move in three dimensions makes it suitable for custom LAAM systems as well, if the plasma cutter head is replaced with a suitable extruder.The gantry system is controlled using the Mach software Mach3 by Newfangled Solutions (Livermore Falls, ME, USA) managed on a local computer.All the experimental specimens for this study were fabricated using G-codes manually written by the author and ported to Mach3.The individual pellets, shown in Figure 2a, are fed into the hopper and then the filament is extruded, as shown in Figure 2b, and then deposited on to the print bed as shown in Figure 2c.During each of these various steps, the polymer is subjected to shearing and elongation, thus increasing the potential for fiber damage.The beads in the present study were composed of 13%, by weight, carbon-fiber-reinforced acrylonitrile butadiene styrene (CF-ABS).

Reduction of Additively Manufactured Bead to a Test Specimen
For the fiber length distribution and aspect ratio characterization, single beads of acrylonitrile butadiene styrene (ABS) pellets with 13% carbon fiber, by weight, were printed at 250 • C on a 95 • C bed. Single beads were selected so as to keep the present focus on the individual deposited bead.Future studies proposed by the authors will focus on multiple deposited beads to study issues presented due to the presence of an interface.Specimens were gathered at the three different points in the printing process indicated in Figure 2, the preprocessed pellet, predeposition extrudate, and the deposited bead.Obtaining the fiber aspect ratio at each stage enables one to identify the most damaging part of the printing process and quantify its effect on the material property predictions.It was found that the volume fraction of fibers was consistent across all three stages of the process and over multiple samples processed over various batches.Thus, it was reasonably assumed that the fibers were well dispersed in the end product, but this did not suggest that there would not exist a small number of local variations in the fiber packing density as demonstrated in the companion study [11].The present paper is focused on local variations in orientation and for simplicity, assumes that the volume fraction of fibers is homogeneous.
For the strength specimens, the 13% CF-ABS was dried for 4 h at 82 • C prior to printing and stored in an environmental chamber with a −40 • C dew point.For printing the specimens, the three extruder temperature zones were set to 200 • C, 205 • C, and 210 • C, and the extruder RPM was set to 2250 with a nozzle-to-bed distance of 3 mm.
After printing, the upper surface of the beads were milled to remove the roughness of the surface of the bead.After machining, the specimens were sectioned with a Buehler IsoMet TM Low Speed Precision Cutter (Lake Bluff, IL, USA).The tensile specimens were 165 mm in length with 25 mm long tabs and a nominal thickness of 3 mm.The dimensions of the compression specimens were 81 mm in length with 38 mm long tabs and a nominal thickness of 3 mm.Specimen dimensions were of a single deposited raster and the length dimension was along the raster deposition dimension.After the specimens were cut, they were then cleaned, and the tabs were adhered using a Hysol 3039A adhesive and cured following the manufacturer's recommended curing cycle.Tensile specimens were then heat-treated at 110 • C for 10-15 min to reduce the warpage induced by milling the top surface and adhesive bonding.For the compression samples there was no measurable warpage, so they were not subjected to any heat treatment.After assembly of the tabs, samples were milled on their sides to ensure a constant, rectangular cross section and dimensional consistency between specimens.Figure 3 shows the prepared tensile specimens.

Experimental Characterization of Fiber-Reinforced Additive Manufactured Composite
The key experimental results of this paper are discussed in this section.Each subsection both describes an experimental setup and presents the results from that experimental setup.The fiber length distribution and aspect ratio study is discussed first.Next, the tensile stiffness and strength results are given.Finally, the compression stiffness and strength results are presented.

Fiber Length Distribution and Aspect Ratio Characterization
The process of characterizing the fiber length distribution consisted of five major steps.The first involved obtaining specimens.Specimens were gathered from the three points in the printing process shown in Figure 2.That is, a pellet, a 25.4 mm extrudate specimen, and a 25.4 mm printed specimen were gathered, with the printed specimen being of primary importance in this study.The second step was to isolate the fibers from each sample without damaging them.The third step was to capture micrographs of the fibers, stitch them together, and measure the fibers using a custom MATLAB code.The fourth and final step was to plot the fiber length distribution, get the number-average and weight-average fiber aspect ratios, and predict other material properties of interest.
Extracting the individual fibers from the matrix was accomplished by using a thermal digestion of each specimen.The complete test setup is shown in Figure 4, and the test setup where pellets were placed in a TA Instruments Q50 Thermogravimetric Analyzer (TGA) is shown in Figure 4a (New Castle, DE, USA).The burn-off procedure consisted of (1) ramping up the temperature to 600 • C, (2) holding it isothermal for 1 h in nitrogen at 600 • C, and (3) holding it isothermal for 10 min in air at 600 • C. The inert, nitrogen atmosphere prevented the carbon fibers from degrading.However, to get rid of remaining residue after step 2, step 3 was added so that the fibers could completely separate from each other.Adding step 3 is not ideal because the fibers may degrade slightly, skewing the measured fiber length distribution.Therefore, a correction factor was applied to the measurement data to mitigate these effects, as discussed in [12].To decrease the effect of cut fibers on skewing the fiber length distribution of the extrudate and printed specimens, 25.4 mm specimens were used since they were relatively large compared to the nominal length of a fiber, 100∼500 µm.These specimens were too large to fit in the TGA machine though, so a Ney Vulcan 3-1750 box furnace (Scotia, NY, USA) was used as shown in Figure 4b.The specimens were enclosed in a Petri dish with an aluminum lid.Two pipes attached to the aluminum lid purged the enclosure with nitrogen.A similar burn-off procedure as for the pellet was repeated for the extrudate and printed specimens.
The third step in characterizing the fiber length distribution was to capture images of the fibers under an MZ7 microscope, (Scienscope, New Haven, CT, USA) shown in Figure 5a, and stitch them together.After burn-off testing, the burn-off specimens were dropped into beakers of acetone and water and dispersed using a Branson Digital Sonifier 450 (Branson, MO, USA).The solutions were poured into Petri dishes before the dispersed fibers had time to settle completely.After pouring the dispersed fibrous solution into a Petri dish, the Petri dish was examined under a microscope and several micrographs were captured and stitched together.Figure 5b shows a micrograph, stitched in Adobe Photoshop CC 2017 and the resulting image was used for the fiber length measuring.After stitching, a custom code written in the MATLAB environment was used to measure the fibers.The code allows one to measure a fiber by performing a mouse-click on both ends of the fiber.First, however, a microscopic ruler must be measured to find the proper pixel-to-millimeter scaling factor.Once a fiber is measured, the length of the fiber is automatically scaled and logged.The code allows a user to measure several fibers at once, close out of the measurement session, and return to it at a later time without losing any data.
The fourth step was to plot the probability distribution function (PDF) of the fiber length and calculate other quantities of interest.Figure 6 shows the PDF of the fiber length, with the PDF being cast in units of µm −1 for the three manufacturing stages.From the data in Figure 6, the number-and weight-average fiber lengths were calculated, respectively, by In the above equations, N i is the number of fibers in the ith histogram bin, and L i is the average fiber length within the ith histogram bin.Similarly, the number-and weight-average aspect ratios were approximated, respectively, as where d n and d w are taken from [12].The results for the fiber diameters, lengths, and aspect ratios for both the weight and number averages are given in Table 1 for this study.Note, these data were corrected for fiber degradation due to thermal heating in the presence of oxygen to remove the polymer matrix as detailed in [12].

Stiffness and Strength Characterization
This section discusses the stiffness and strength testing.The results for the tensile properties are given first, followed by the compression results.

Tensile Testing
A Test Resources 810 Series fatigue tester with a F2500-B load cell was used for the tensile testing.Figure 7a shows a mounted specimen with an extensometer for measuring strain.The specimens were tested at a speed of 0.059 in/min with a data collection rate of 20 Hz and at a temperature of approximately 21 • C (70 • F).Three width and three thickness measurements were taken for each of the five specimens used in the tensile study using micrometers in or near the gage region prior to testing.The width and thickness measurements were averaged to obtain the cross-sectional area of each specimen, from which the stress was calculated.A typical tensile stress-strain plot is shown in Figure 7b.The effective longitudinal tensile stiffness of the tensile specimens was defined as E e f f LT , where the LT subscript indicates the longitudinal tensile component and the e f f superscript indicates the bulk stiffness of the locally varying stiffness.The effective longitudinal tensile stiffness was defined over the strain range 0.003-0.005.In addition, the effective longitudinal yield strength, defined as the term σ e f f LT , was determined by using an offset as defined in ASTM D638 [13].In the present study, a 0.2% offset was selected due to the initial nonlinear loading region and to align with the literature for which the models were based upon.That is, a line parallel to the elastic region was extended from (0, 0.002) up to the stress-strain curve, and the stress at the intersection was taken as σ e f f LT .A toe region sometimes appeared at the beginning of the stress-strain curves, which can be seen in Figure 7b.Therefore, the data were horizontally shifted to where the E e f f LT lines would intersect the origin, prior to applying the 0.2% offset method.A summary of the results of four tests is given in Table 2, where x, s n−1 , and C V are, respectively, the average, standard deviation, and coefficient of variation for the five specimens studied.

Confined Compression Testing
The same Test Resources 810 Series fatigue tester with the same F2500-B load cell was also used for compression testing.Figure 8a shows a specimen in a modified ASTM D695 compression test fixture (Boeing BSS 7260) [14], which was placed between two compression platens mounted on the testing machine.Three width measurements of the five specimens were taken for each compression specimen, two on either end and one in the middle.The compression specimens had tabs covering most of their length, so only one thickness measurement was taken in the middle of the specimens, and this was done with calipers since the micrometer anvils did not fit between the tabs.The average width and thickness measurements were used to calculate the specimen's cross-sectional areas from which the stress was derived.The length of each of the five specimens was also measured and used to calculate the approximate strain.Before initiating a test, the upper compression platen was lowered in contact with the specimen.Each test was conducted at 0.051 in/min A typical compression stress-strain curve is shown in Figure 8b.Five compression tests were performed.The toe region of each curve was trimmed, this time by determining the stiffness over the stress range 25-35 MPa and horizontally shifting the data.In addition, the effective longitudinal compression yield strength, σ e f f LC , was determined using a 0.2% strain offset.The results of the compression tests are given in Table 3.The ultimate compression strength is not reported because the upper platen would often come into contact with the top of the metal fixture before the ultimate failure of a specimen.In addition, it is noteworthy that both the experimental compression stiffness and strength surpass the experimental tensile stiffness and strength, respectively.Over the past few decades there has been considerable work in the area of fiber motion kinetics and the resulting structural performance of the processed composite.The following section pulls together a summary of the high points along with various components that were used in the present study.Consider a fluid flow domain with the velocity profile v(x) and velocity gradients L(x), where x denotes position.Jeffery's foundational work in expressing fiber motion for a dilute suspension in such a domain describes the time rate of change in orientation of a single, rigid, inertialess ellipsoid [2], In the above equation, the material derivative of the unit vector p, which points along the long axis of the fiber, is used.That is, ṗ = dp dt = ∂p ∂t + v • ∇p.In addition, the vorticity tensor is defined as W = 1  2 L − L T , and the rate of deformation tensor is similarly defined as D = 1  2 L + L T .ξ is the fiber geometric term, which allows Equation ( 5) to be used for other axisymmetric shapes besides ellipsoids as long as an equivalent ellipsoidal aspect ratio, r e , is properly defined.Expressing the fiber geometric term ξ in terms of r e gives (see, e.g., [15]) Jeffery's equation as expressed in Equation (5) may be used in dilute fiber solutions, but a fiber interaction model must be used for concentrated solutions.The concentrated regime may be defined as V f > (d/L), where V f is the fiber volume fraction and d and L are the fiber diameter and length.For the carbon fibers in the CF-ABS used in the present study d/L = 6.4 µm/278 µm = 2.3%.The CF-ABS used in the present study had V f = 8.11% (this corresponds to a 13% weight fraction) and 8.11% > 2.3%, thus the CF-ABS flow was considered concentrated and therefore a fiber interaction model was used.In addition, the orientation tensor approach popularized by Advani and Tucker [3] was used to describe the orientation state of a population of fibers and save orders of magnitude in computational time.Cast in terms of orientation tensors, there are multiple options for the proper characterization of the fiber interaction kinetics, such as the retarding principal rate model from Tseng et al. [16], the interaction coefficient approach of Folgar and Tucker [4], and the anisotropic rotary diffusion model of Phelps and Tucker [6].In the present study, we focused on the fiber interaction model of Wang et al.
where the second-and fourth-order orientation tensors, as discussed by Advani and Tucker [3], are defined as In the above, ψ(p) is the orientation probability density function.In addition, κ is the slowness factor as discussed by Wang et al. [5], C I is the fiber interaction constant introduced by Folgar and Tucker [4], γ is the scalar magnitude of D, L 4 = ∑ 3 i=1 λ i e i e i e i e i , M 4 = ∑ 3 i=1 e i e i e i e i , and λ i and e i are the ith eigenvalue and eigenvector of A. Based on the work of Bay and Tucker [17], we take Note, that a r is the geometric fiber aspect ratio (a r = L/D) and is related to the equivalent ellipsoidal aspect ratio of a cylindrical fiber as given by Zhang et al. [15].r e = 0.000035a 3 r − 0.00467a 2 r + 0.764a r + 0.404 (10) Furthermore, based upon earlier work by the present authors (see, e.g., [8]), a value for the slowness parameter of κ = 0.125 was selected.Finally, the eigenvalue-based orthotropic closure discussed by Wetzel [18] and VerWeyst [19] was used to estimate A 4 due to its accuracy and efficiency.Equation ( 7) and all the equations in this section which it depends on were custom-coded by the authors in MATLAB, such that all flow kinematics were performed in the MATLAB environment.This was conducted along several streamlines in the flow domain to capture a fine-resolution solution for the fiber orientation state across the flow domain.

Fiber Micromechanics
The linear elastic nature of a material is often described using Hooke's law to relate the stress, σ ij , to the strain, ij , through the stiffness, C ijkl , or the compliance, S ijkl , as Fibers and matrix have different values for their respective stiffness and due to the length scales considered, an effective stiffness is considered for the composite, C ijkl .This effective stiffness, also termed the homogenized stiffness, due to its relationship to the associated underlying unidirectional stiffness tensor of the composite, Cijkl , and the fiber orientation distribution function, is expressed as The underlying unidirectional stiffness tensor of the associated composite Cijkl may be found in a variety of methods, the Halpin-Tsai and rule of mixtures often being the most popular due to their ease of implementation.In the present study, the authors used the Tandon and Wang [20] model with the closed-form solution suggested in Tucker and Liang [21] and with the mathematical form expressed in Zhang [22].The choice of the Tandon and Wang model for the underlying unidirectional stiffness tensor of the associated composite was based upon the results from [21] in comparison to other available unidirectional composite models.Thus, with knowledge of the underlying stiffness tensor Cijkl and the fiber orientation state from Equation ( 7), the effective stiffness tensor at a point could be computed.Advani and Tucker [3] provided a mathematical form of the effective stiffness tensor in terms of the orientation tensors as where the coefficients B i can be cast in terms of the underlying unidirectional stiffness tensor Cijkl as presented in [3].Observe that in Equation ( 13), the fourth-order orientation tensor A ijkl appears.The same orthotropic closure of VerWeyst [19] and Wetzel [18] used in Equation ( 7) is also used in Equation ( 13).The Tsai-Wu [23] model for the onset of failure was used and is expressed in terms of the second-order and fourth-order strength tensors, respectively, f i and F ij , as In Equation ( 14), a contracted notation is used where the index pairs [11, 22, 33, 23, 13, 12]  are replaced by [1, 2, 3, 4, 5, 6].For example, the stress tensor in index notation σ ij can be ex- pressed in contracted notation as σ m where [σ 11 , σ 22 , σ 33 , σ 23 , σ 13 , σ 12 ] → [σ 1 , σ 2 , σ 3 , σ 4 , σ 5 , σ 6 ].
In Equation ( 14), the summation convention continues to be used where {m, n} ∈ {1, 2, . . . ,6}.Alternatively, the failure envelope in terms of the contracted strain tensor is where g m and G mn are, respectively, the contracted second-and fourth-order strength tensors.One can obtain the relationship between the tensors f m , F mn , g m , and G mn as (see e.g., [23]) where the summations on the indices o and p are from one to six, and C mn is the contracted form of the stiffness tensor discussed in Equation ( 13).As discussed in [8], the strength tensor relationship for the failure envelope can be recast in terms of the underlying unidirectional strength tensors, fm , Fmn , ḡm , and Ḡmn , where the bar indicates the unidirectional composite property, using a homogenization similar to that of the stiffness tensor of Equation ( 12) as For a unidirectional-fiber-reinforced composite with all fibers aligned along the x 1 axis, the unidirectional strength tensors fm and Fmn may be expressed in terms of the tensile and compression strengths in the form where the first subscript of the term σ indicates the longitudinal (L) axis along the fiber or the transverse (T) axis, and the second subscript indicates a tensile (T) load or compression (C) load.For example, σTC is the experimental transverse strength found from compression loading.In the present study, we assumed the tensile and compression behavior for the transverse load direction of the underlying unidirectional composite was identical, σTT = σTC = σ m , where σ m is the experimentally obtained strength of the matrix.We also assumed the corresponding shear strength of the matrix was found as τ = σ m / √ 3. Using the form suggested in [24], the remaining nonzero terms for Fmn were expressed as F12 = F13 = − F11 /4 and F23 = − F22 .Then, using Equation ( 16), the unidirectional strength tensors ḡm and Ḡmn were obtained.As shown in [8], the homogenized fourth-order strength tensor was found using an identical form as Equation ( 13) in terms of the orientation tensors.The homogenized second-order strength tensor in terms of the second-order orientation tensor and second-order strength tensor of the underlying unidirectional SFRP was found as (see, e.g., [25]) In the present study, the composite was considered at the edge of failure when the equality of Equation ( 17) was satisfied.For simulation purposes, a loading state was imposed on the composite and the following expression was evaluated ϕ( ) = g m m + G mn m n (20) where is the strain state from the imposed load on the structure.The SFRP is not predicted to fail as long as ϕ( ) < 1.When ϕ( ) reaches one, the SFRP is predicted to fail.The longitudinal tensile strength σLT for the unidirectional composite was found using the modified rule of mixtures suggested by Van Hattum and Bernardo [25] as where the above expression is only valid when the fiber length, L, is less than the critical fiber length, L c .The fibers in the current study had a typical fiber length of 200-300 µm in the final processed composite structure, whereas the critical fiber length assumed in this study, L c 0.9 mm (taken from [25]), was much greater than the largest fiber length identified in Figure 6.The term σ m is the stress of the matrix at the failure strain of the fiber and is experimentally obtained once the strain at failure for the fiber is known.The longitudinal compression strength was cast in the form suggested by Hayashi and Koyama [26] as expressed by Bajracharya et al. [27] σLC where E f and E m are the moduli of, respectively, the fiber and the matrix, and , where * m is the strain at which the matrix yields.The parameter χ is expressed in terms of the weight-average length of the fibers, L w , and the critical fiber length, L c , where χ = L w /(2L w ).
It is well known that the presence of porosity in the composite reduces the stiffness and the strength.In the thesis of Nargis [10], the average void fraction identified in her specimen was 13.84%, a value that was used in this paper.Using the estimation of Zhang et al. [28] that the knockdown factor for the stiffness and the tensile strength is the same and is a function of the void volume fraction, the effective stiffness, E e f f , and strength, σ e f f , were expressed as where in the present study n = 2 was used.Equation ( 23

Deposition Flow Domain
The flow domain in the present study is depicted in Figure 9 with the dimensions selected based upon the nozzle shown in Figure 2. Solutions for the Navier-Stokes equations were performed using COMSOL Multiphysics.The interior walls of the nozzle were assumed to be no-slip boundary conditions along with the surface of the print bed, whereas the surfaces of the melt flow that are outside of the confined nozzle were considered to have slip boundary conditions.The geometry of the flow domain at the leading edge of the flow was constructed such that the surface was stress-free, and the normal component of the velocity was also zero using the approach from Heller et al. [29] to identify the surface using an iterative approach.Once the velocity field was obtained from the finite element simulation, the spatially varying fiber orientation state, A ij (x), was obtained using Equation (7) with C I = 1.9 × 10 −3 , κ ∈ {0.05, 0.2}, and λ cast as a function of a rw of the 3D-printed bead as given in Table 1.The fiber equations of motion were solved numerically using a custom program created in the MATLAB environment written by the authors.Then, using the material properties for the fiber and the matrix listed in Table 4, Equation (20) was evaluated across the flow domain using a custom script created by the authors that took the spatial location and returned the failure parameter at the specified location.
Fiber volume fraction [37] † † These references are for the equations used to calculate the property.‡ Fiber failure strain taken from [32] and using the stress-strain curve from [35] for ABS, the matrix stress was obtained.† † Refers to experimental results in the present paper.

Structural Simulation
The orientation state at the end of the flow domain, A ij (x 1 = 0, x 2 ), shown at the origin on the bottom left of Figure 9, was taken as the steady-state orientation and was used for all structural simulations.This orientation state was then projected along the direction of printing to form a specimen with an orientation that changed along x 2 but was constant along x 1 .The tensile and compression specimen model domains are shown in Figure 10.Each domain is of a single deposited bead.The bead subjected to tension or compression is allowed to freely slide in the vertical direction, except in the lower, right corner.That is, the displacement u 2 is zero in the lower, right corner and unconstrained everywhere else.The right edge of the bead is fixed horizontally, i.e., u 1 = 0, and equal to some static displacement value in the x 1 direction, δ, on the left edge.
For the applied displacement δ, the spatially varying strain, (x), was computed from the finite element results along with the surface stresses.The effective stiffness values for tension or compression, E e f f LT and E e f f LC , respectively, were obtained by taking the average of the stress divided by the applied strain on the right end of the loaded specimen.Then, using the spatially varying strain values along with Equation (20), the spatially varying value for ϕ( ) was obtained at every point in the solution domain.If the maximum value for ϕ was less than one, then the applied displacement δ was increased, whereas if ϕ at any point in the domain was more than one, the applied load was decreased.This process continued until the maximum spatial value for ϕ was one, thus suggesting that the member was at the onset of failure.At this point, σ e f f LT and σ e f f LC were calculated.A representative solution is shown in Figure 11a,b for a displacement resulting in a tensile and compression failure, respectively.
The finite element process described above was extended to study the stiffness and strength response over a range of values for κ.This first study fixed the fiber length and diameter to be that in the deposited bead, yielding a weight-average fiber aspect ratio of 38.8 as given in Table 1.The fibers were assumed to be randomly orientated initially for the flow domain of Figure 9, and once deposited, the structural domain shown in Figure 10 was analyzed for the stiffness and strength estimates for tension and compression.Using the fiber aspect ratio of 38.8 and a weight fraction of 13% (i.e., a volume fraction of 8.11%), the interaction coefficient from Equation ( 9) evaluated to C I = 1.9 × 10 −3 .The range for κ was taken to be 0.05 to 0.2, within the range suggested by Wang et al. [5].The remaining modeling inputs in this study are given in Table 4.
The results for the stiffness and strength parameters as a function of fiber slowness parameter κ are provided in Figure 12 for both the scenario when the porosity is neglected and the scenario when the porosity is accounted for in the structural simulations.Notice in Figure 12 that the tensile and compression stiffness values are graphically indistinguishable regardless of the porosity.This is by construction in the model, as the stiffness response is assumed to be equivalent in tension and compression.Conversely, the compression strength is measurably more than that of the tensile strength.Observe that the stiffness and strength increase as a function of the increasing value of the slowness parameter.This correlates to the increase in the orientation.As the slowness parameter, κ, increases from zero to one, the rate of alignment increases as well; thus, a higher value of the parameter κ correlates with a higher value of the alignment along the flow direction, and thus a higher value for the stiffness and strength for fibers that are stiffer and stronger than the surrounding matrix.It is also worth noting the significant reduction in the part performance as a result of porosity.This highlights the importance of accounting for porosity and the need to identify a means to reduce porosity generation in the deposited bead during processing.The last modeling study focused on the importance of including the reduction in the fiber length in the structural simulations.Taking the three length distributions presented in Table 1 for the case of the pellet, the extrudate, and the deposited bead, the weightaverage length is shown to decrease from 424 µm to 348 µm to 279 µm.Performing the flow simulation with the corresponding fiber aspect ratios, followed by the structural simulations led to the stiffness and strength predictions provided in Table 5. Observe that by properly including the reduction in the fiber length, there is a reduction of nearly 12% in the stiffness, a reduction of over 30% in the tensile strength, and a reduction of over 14% in the compression strength.

Comparison between Model and Experiment
Table 6 summarizes the experimental results of this study and Table 7 summarizes the model predictions for comparison.These results show that the model predictions are on the high side of the experimental results.However, there are several possible culprits for this.On the experimental side, the top surfaces of the specimens were machined down to 10-20% of the specimen height to eliminate waviness and to ensure uniform cross sections.While uniform cross sections are desirable, the skin of the specimens contained a high fiber alignment in the load-bearing direction.Thus, removing this skin from the tops of the specimens could have potentially caused a noticeable drop in the mechanical properties.
On the modeling side, Table 7 shows that a relatively large range of predictions can be found depending on the parameters fed into the models.The lower range for e f f LT and E e f f LC should be selected in the present study based on the recommendations found in recent fiber interaction model papers to decrease the value of κ. κ has a greater effect than C I , but there are other parameters that were estimated in the present study that could also influence the accuracy of the models.These include the initial orientation state and the constituent material properties in Table 4. Thus, it is important to zero in on accurate modeling parameters when using the proposed models in order to tighten the prediction window.In addition, the recent work by Nargis and Jack in [11] showed that the porosity was found to be significantly larger than anticipated, coupled with a broad spectrum of void sizes.Thus, accounting for porosity and selecting the lower range for the slowness factor κ, the predicted value for the compressive stiffness, E e f f LC , came out to be 4.2 GPa, whereas the measured value was 4.06 GPa; the predicted value for the tensile stiffness, E e f f LT , was 4.2 GPa, whereas the measured value was 2.71 GPa; the predicted value for the compressive strength σ e f f LC , was 56 MPa, whereas the measured value was 59 MPa; and the predicted value for the tensile strength, σ e f f LT , was 30 MPa, whereas the measured value was 27 MPa.Each of these values, with the exception of the tensile stiffness, is well within an acceptable range.This study highlights a limitation of the present modeling that should be addressed in a future study.The present model assumes that the tensile and compressive stiffness are identical, which is not the case based upon the results presented.Thus, a future expansion of the model to allow for discontinuous tensile and compressive behavior will be required prior to full-scale industrial implementation and will be the scope of a future study.

Conclusions
Having so many modeling parameters increases the difficulty to fully validate the proposed models.However, the results provided indicate the modeling framework yields results for the effective strength and stiffness that align with experimental characterization to within ∼1% and ∼10% for the longitudinal compressive and tensile strength, respectively, and to within ∼3% and ∼50% for the longitudinal compressive and tensile stiffness, respectively.In addition, a framework has been developed that allows for improvement by model substitution as more accurate modeling inputs are found.The modeling methodology goes from the properties of the individual fibers and surrounding matrix, through the fiber kinematic equations of motion through the final part performance.Results show the sensitivity of the final predicted response to subtle changes in the fiber interaction behavior as well as a reduction in the final part performance through the incorporation of the known porosity of the final processed part.
For future work on the experimental side, improvement in specimen preparation should be sought.For example, nonmachined, multibead tensile and compression specimens could be used to reduce the effect anomalous defects might have, or injection-molded specimens (which have less void content) could be used to reduce the effects of porosity.On the modeling side, in general, efforts to isolate and independently verify each underlying model and model input would be helpful.The fiber orientation could be physically measured as in Nargis' work [10], for example, and directly input into the stiffness and strength models to validate these models independently from a fiber orientation model.In addition, it is well understood that polymers are sensitive to changes in temperature, and performing tests across a temperature range would be of interest.The present model is linearly elastic, and extending the work to viscoelastic models would be of interest.A first pass at a 3D model for stiffness prediction which incorporated Nargis' results was performed in [12], but this area still needs more work.In addition, the effects of porosity need to be accounted for better, both in modeling the fiber orientation state and predicting the final mechanical properties as a function of the orientation state.The flow and orientation kinematics relationship was decoupled in the present study, but coupling them may make a difference in the predicted final orientation state.Another consideration for a future study may be a fiber aspect ratio distribution, rather than a single-value aspect ratio.

Figure 4 .
Figure 4. Burn-off testing setups.(a) TA Instruments Q50 TGA machine with a 13% CF-ABS pellet.(b) Ney Vulcan 3-1750 box furnace and 1-inch extrudate and printed specimens.These specimens were enclosed in a Petri dish with a custom-machined aluminum lid and inserted into the furnace.

Figure 5 .
Figure 5. (a) An MZ7 microscope, used to capture micrographs of the fibers in a Petri dish.(b) Fibers being measured in a stitched micrograph.Green graph paper was placed under the Petri dish to provide a helpful frame of reference.The measurements in blue were measured in a current session, whereas the green measurements were from a previous measurement session.
(1.3 mm/min, as called for by ASTM D695) with a data collection rate of 20 Hz and an environmental temperature of about 21-22 • C (70-72 • F).

Figure 8 .
Figure 8.(a) A compression specimen mounted in a modified ASTM D695 compression test fixture (a Boeing BSS 7260) and between two platens.(b) A typical compression stress-strain plot for a 13% CF-ABS specimen.
) was used to find the effective stiffness values E

Figure 9 .
Figure 9. Finite element modeling domain at the tip of the LAAM nozzle including the deposition region.

Figure 10 .
Figure 10.Specimen domains used for the structural analysis for the (a) tensile test and (b) compression test.

Figure 11 .
Figure 11.Structural simulation results showing the failure parameter ϕ(x) for the (a) tensile test and (b) compression test.

Figure 12 .
Figure 12.Effective stiffness and strength predictions of the deposited bead for a r = 38.8.(a) Effective stiffness and (b) effective strength.Table 5. Tensile and compression stiffness and strength predictions as a function of the weight-average aspect ratios from Table 1 along the flow domain for κ = 0.125.Parameter E e f f LT = E e f f LC (GPa) σ e f f LT (MPa) σ e f f LC (MPa) L w = 424 µm (pellet) 5.58 44.8 75.4 L w = 348 µm (extrudate) 5.34 39.3 71.2 L w = 279 µm (bead) 5.00 34.4 66.1

Table 1 .
Fiber length and aspect ratio data, corrected for fiber degradation.

Table 4 .
Properties of ABS matrix and carbon fiber used for modeling.

Table 6 .
Summary of experimental results, with the error taken from the standard deviation of experimental results.

Table 7 .
Tensile and compression stiffness and strength predictions of the deposited bead using the RSC model with a random initial orientation state (A ij = 1 3 δ ij ), a rw = 38.8, and κ ∈ {0.05, 0.20}.Equation (23) was used to generate the second row from the first row.