Coupled thermo-mechanical numerical model to minimize risk in large-format additive manufacturing of thermoplastic composite designs

The collapse of deposited thermoplastic composite material under self-weight presents a risk in large-format extrusion-based additive manufacturing. Two critical processing parameters, extrusion temperature and deposition rate, govern whether a deposited layer is stable and bonds properly with the previously deposited layer. Currently, the critical parameters are determined via a trial-and-error approach. This research work uses a simplified physics-based numerical simulation to determine a suitable combination of the parameters that will avoid the collapse of the deposited layer under self-weight. The suitability of the processing parameters is determined based on the maximum plastic viscous strains computed using a sequentially coupled thermo-mechanical numerical model. This computational tool can efficiently check if a combination of temperature and extrusion rate causes layer collapse due to self-weight, and hence minimize the manufacturing risk of large-format 3D-printed parts.


Introduction
Large-format extrusion-based additive manufacturing has become increasingly widespread in research and is headed towards industrial use. Srivastava and Rathee [1] reviewed the recent trends in additive manufacturing and highlighted that the size of the parts that can be manufactured using large-format 3D equipment has increased consistently over the years. BAAM (big area additive manufacturing), one of the first large-format 3D printers, was completed in January 2015. Love and Duty [2] reported that BAAM increased the maximum size of parts that can be manufactured using additive manufacturing by over 10 times, increased the manufacturing speed by over 1000 times, and reduced the cost by over 100 times compared to the next largest 3D printer available at the time. Polymer extrusion-based large-format additive manufacturing process has since been used in a variety of applications. Kunc et al. [3] used large-format additive manufacturing for producing composite tooling. Love [4] presented the use of large-format additive manufacturing for producing automotive bodies. Post et al. [5] used polymer extrusion-based large-format additive manufacturing to manufacture a boat hull mold. Bhandari et al. [6] and Roschli et al. [7] demonstrated the use of large-format 3D-printed forms for manufacturing precast concrete formwork. Bhandari et al. [8] manufactured a highway culvert diffuser using a large-format polymer extrusion-based additive manufacturing process. Davids et al. [9] used a large-format additively manufactured mold to manufacture fiber-reinforced polymer girders for highway bridges. With large-format polymer extrusion-based additive manufacturing technology moving towards widespread industrial applications, new techniques need to be developed to ensure the production of consistently high-quality additively manufactured parts. Process qualification techniques can aid additive manufacturing of consistent parts that meet acceptance criteria. ISO/ASTM 52903-2:2020 [10] describes a method for defining requirements and assuring component integrity for plastic parts created using material extrusion-based additive manufacturing processes.
While scaling up enables extrusion-based additive manufacturing to produce larger parts at faster speed and reduced costs, it also introduces a new set of challenges to extrusionbased additive manufacturing. Shah et al. [11] highlighted the challenges of designing large stiff frames, controlling the temperature of the surrounding environment, maintaining an effective melt pool of the printing material for extrusion, and supplying mass of material for print. Roschli et al. [12] discussed the design problems encountered in large-format additive manufacturing including challenges during the preparation of a print, bead width constraints, layer time limitations, and the challenge of maximizing throughput.
One of the major challenges in large-format additive manufacturing is the selection of a suitable layer time. Roschli et al. [12] discussed the concepts of minimum layer time and maximum layer time in polymer extrusion-based large-format additive manufacturing. If a layer does not cool down sufficiently, it does not become rigid enough to support another layer above it. The time required for the layer to cool down enough to be able to support the layer above it is called minimum layer time. If the layer cools down too much, the bonding with the layer above it is hindered and the layers do not adhere properly. The time for which the bonding between the layers remains at acceptable levels is called maximum layer time. A suitable layer time is bounded by the minimum layer time and the maximum layer time. Compton et al. [13] created a 1D thermal model to predict temperature evolution and found that adopting a suitable layer time minimizes the likelihood of build failure. Borish et al. [14] used thermal cameras to monitor the temperatures of the deposited beads and ensure that sufficient cooling has occurred to avoid build failure. Fathizadan et al. [15] and Wang et al. [16], in two different studies, used real-time thermal analysis to control the layer time and ensure successful manufacturing. Bhandari [17] discussed the toolpath generation process in extrusion-based additive manufacturing. Once the 3D model of the part is sliced and a toolpath is generated, two critical process parameters can be changed that affect the minimum and maximum layer time. The two critical process parameters are extrusion rate and extrusion temperature. Extrusion temperature is the temperature at which the molten thermoplastic material is extruded from the nozzle of the AM equipment. The extrusion rate is the rate at which the molten thermoplastic material is deposited during the additive manufacturing process.
The concept of layer time is based on the change in mechanical properties of polymers with temperature. The mechanical properties of the polymers used in extrusion-based additive manufacturing change with temperature. Tan et al. [18] reviewed different polymers used in thermoplastic polymer extrusion-based additive manufacturing and broadly divided them into two categories, namely, semi-crystalline and amorphous polymers. Hiemenz and Lodge [19] and Ishinabe [20] highlighted the several transitions that thermoplastic polymers undergo when cooled down from high extrusion temperatures to low room temperatures. Two notable transitions that affect the extrusion-based additive manufacturing process are melting and glass transition. Glass transition is marked by a sharp change in the stiffness of the polymer for both semi-crystalline and amorphous polymers. McIlroy and Olmsted [21] demonstrated that inter-penetration depth and re-entanglement are arrested by the glass transition marking the lower limit of the temperature at which bonds can develop between deposited beads. Melting transition is prominent in semi-crystalline polymers where the crystalline regions in the polymer melt. Since there are no crystalline regions in amorphous polymers, melting transition is not observed in amorphous polymers. Stansbury [22] pointed out that amorphous polymers are more amenable to additive manufacturing since amorphous polymers do not undergo significant dimensional contraction associated with polymer crystallization as in semi-crystalline polymers. In amorphous polymers, a steady decrease in storage modulus and viscosity is noted with increasing temperature.
Dynamic mechanical analysis (DMA) is carried out to characterize the viscoelastic behavior of a polymer. Sinusoidal stress is applied to the polymer and the mechanical response of the polymer to the stress is measured. The elastic response of the polymer is characterized by the storage modulus (G′), which accounts for the potential energy stored by the polymer. The viscous response of the polymer is characterized by the loss modulus (G″), which accounts for the energy lost by the polymer. The elastic response is in phase with the applied sinusoidal stress and the viscous response is out of phase with the applied sinusoidal stress. Tan δ is the tangent of the phase angle between the applied stress and the viscoelastic response of the polymer and is equal to the ratio of loss modulus to the storage modulus. Menard [23] pointed out that tan δ is an indicator of how efficiently the material loses energy to molecular rearrangements and internal friction. Based on tan δ at different temperatures, Duty et al. [24] identified the four stages that a polymer undergoes during the cooldown after extrusion during additive manufacturing, namely, viscous liquid, viscoelastic liquid, viscoelastic solid, and elastic solid. Figure 1 shows the different stages in temperature-dependent viscoelastic behavior of a typical polymer during AM process. The polymer behaves as a viscous liquid at high extrusion temperatures. As it cools down, the polymer begins to solidify and the polymer exhibits more solid-like behavior with decreasing temperature.
1. Viscous liquid: at near extrusion temperature, the fluidlike nature of the polymer is dominant and is characterized by high values for tan δ. The loss modulus is very high compared to the storage modulus for viscous liquid. 2. Viscoelastic liquid: As the polymer cools down, the elastic nature of the polymer becomes significant and is categorized as a viscoelastic liquid. Maxwell model [25] is used to calculate the plastic viscous strain in the material. 3. Viscoelastic solid: As the polymer cools down even further, the elastic behavior of the polymer is more prominent than the viscous behavior. Kelvin-Voigt model [25] is used to calculate the viscoelastic strains in the material. 4. Elastic solid: When the polymer cools down below its glass transition temperature, the viscous behavior is insignificant for short-term strain and elastic behavior dominates the overall short-term deformation of the polymer.
Duty et al. [24] chose a value of tan δ greater than 10 to mark the viscous liquid behavior region, between 1 and 10 for the viscoelastic liquid region, between 0.1 and 1 for the viscoelastic solid behavior region, and less than 0.1 for the solid behavior region for polymers used in extrusion-based additive manufacturing.
Byard et al. [26] raised concerns about environmental issues regarding the disposal of massive waste and possible toxic byproducts with the increase in the use of polymer extrusion-based additive manufacturing and presented recycling as a potential solution. Park et al. [27] highlighted bio-based polycarbonate to address the environmental issues. Mulhaput [28] discussed several renewable bio-based thermoplastic composites that have been developed to reduce the consumption of petrochemical-based polymers and reduce waste. Khoramnejadian et al. [29] presented bio-based plastics as a way to reduce municipal solid waste. Tümer et al. [30] and Wickramasinghe et al. [31] highlighted that PLA (polylactic acid) thermoplastic composites are the most used bio-based thermoplastic composites used in additive manufacturing. Mazzanti et al. [32] emphasized that PLA is a semi-crystalline thermoplastic polymer and its degree of crystallization depends on the cooling rate during its processing as well as the reinforcement in its composites. Choo et al. [33] created a heat retention model for large-area additive manufacturing and showed that in large-format additive manufacturing, the cooling rate is slower compared to smallscale additive manufacturing due to the smaller surface area to volume ratio. Quelho de Macedo et al. [34] highlighted the fact that the crystallization process for PLA is dependent on cooling rate and slower cooling leads to a greater crystallinity in the additively manufactured PLA part. During the large-format additive manufacturing process, the PLA composites cool down slower and undergo significant crystallization resulting in a thermal history-dependent mechanical behavior.
Different methods have been used to monitor the thermal history of the polymers and the effects of corresponding changes in mechanical properties during the additive manufacturing process. Yin et al. [35] inserted thermocouples into 3D-printing specimens to measure temperature at the interface between layers to predict interlayer diffusion in the specimens. Borish et al. [14] used thermal imaging to monitor the single-layer time alteration in large-format additive manufacturing. Fathizadan et al. [15] used a set of thermal images to generate temperature data on the surface and used a control model to predict the best time for printing the next layers. Wang et al. [16] used temperature data from thermal images to enable real-time control of largeformat 3D printing. Oleff et al. [36] carried out a state-of-art review of process monitoring for the extrusion-based additive manufacturing process that highlighted the importance of time-temperature data in the manufacturing process. Numerical models have also been used to determine the temperature changes of additively manufactured parts during the printing process. The thermal models solve the heat equation, shown in Eq. 1, in combination with the changing geometry of the part during additive manufacturing to generate the thermal history of the additively manufactured part. where is the density (kg/m 3 ), c p is the specific heat capacity (K/kg/K) , T is the temperature (K), k x , k y , k z is the conductivity in the x, y, and z directions, respectively (J/s/m/K), and Q is the volumetric heat flow (J/m 3 /s). The temperature-time history is used in combination with the temperature-dependent material properties to determine the deformations during the thermoplastic polymer extrusion-based additive manufacturing process. Sun et al. [37] described a mathematical model and used the model to numerically calculate the temperature of 3D-printed parts and to model the bond formation between ABS (acrylonitrile butadiene styrene) filaments. Wang et al. [38] used a simplified mathematical model to calculate the temperature and the resulting warp deformation in 3D-printed ABS parts. Zhang et al. [39] used an adaptable, boundary-adjusting finite-difference method to numerically generate the thermal history of a 3D-printed PLA (polylactic acid) part. Stockman et al. [40] used the finite-difference method to create a thermal model tailored for additive manufacturing. Ji and Zhou [41] described the use of a finite-element model with temperaturedependent material properties. Finite-element method-based models have been used recently to simulate the thermal history of polymer extrusion-based 3D-printed materials. Some of the models are layer-by-layer simulations that use simplified assumptions about the thermal behavior of materials. D'Amico and Peterson [42] presented an adaptable FEA simulation of material extrusion additive manufacturing in 3D. Zhou et al. [43] presented a finite element-based temperature analysis considering the change in geometry during the additive manufacturing process. The more recent finite-element modeling approach uses the concept of element activation, whereby an activation variable is multiplied with the constitutive equations to simulate material deposition. El Moumen et al. [44] modeled the temperature and residual stress fields during the 3D printing of polymers using the concept of activated elements. The activation variable can vary from zero to one depending on whether the region defined by finite-element mesh is empty, partially filled, or completely filled with the deposited material at a given time. Zhou et al. [45] used voxelization-based modeling to activate elements and perform thermo-mechanical finite-element analysis. Brenken et al. [46] used progressive element activation where elements can be partially activated based on the temporal and spatial position history of the deposition head from the physical AM process. Das et al. [47] discussed highly accurate but computationally expensive models that couple momentum, heat, and mass transfer. Mishra et al. [48] described an implementation of viscosity and density models for improved numerical analysis of melt flow dynamics in the nozzle during extrusion-based additive manufacturing. Colon Quintana et al. [49] discussed a roadmap for characterizing thermo-mechanical properties that should be considered for high-fidelity numerical simulations of the large-format additive manufacturing process. Bertevas et al. [50] presented a smoothed particle hydrodynamic (SPH) model to predict the orientation of fiber orientation in a polymer extrusion-based additive manufacturing process. The computationally expensive model accounted for the flow of the extruded thermoplastic material during the additive manufacturing process and the orientation of fiber. Afriasiabi et al. [51] used an SPH model to simulate a laser powder bed fusion additive manufacturing process. Bhandari and Lopez-Anido [52] presented a simplified and computationally inexpensive discrete-event simulation-based thermal model to predict the thermal history of a 3D-printed part and compared it with the results from the finite-element model. The discrete-event simulation results were found to be comparable in accuracy with the finite-element model results. The discrete-event simulation was faster by two orders of magnitude in generating the thermal history compared to the comparable implicit finite-element model.
A computationally efficient numerical model enables the generation of plastic viscous deformation of the extruded material layers due to the weight of the layers above it, even for large parts. Based on the plastic viscous deformation calculated by the numerical model, a range of critical process parameters (extrusion temperature and deposition rate) that avoid layer collapse can be determined. Such a numerical model can be used for expediting process qualification by aiding additive manufacturing of dimensionally consistent parts and preventing layer collapse in both amorphous and semicrystalline thermoplastic composite material systems. The objectives of the research work presented in this manuscript are to: 1. Present a simplified, computationally efficient physicsbased numerical model that can reduce risk in large-format additive manufacturing by predicting whether layer collapse will occur for a given set of printing parameters. 2. Apply the numerical model to generate a design envelope of acceptable values that prevent layer collapse for critical printing parameters (extrusion temperature and deposition speed) for a given print job. 3. Demonstrate the use of the simplified physics-based simulation to prevent excessive plastic viscous deformations in a large-format additively manufactured biobased polymer composite part.

Materials and material property characterization
HiFill ® PLA 1816 3DP from Techmer PM was used as the feedstock material for additive manufacturing. The natural fiber-reinforced polylactic acid (PLA) thermoplastic composite is bio-based and consists of 20% wood flour filler by weight.
The thermal and mechanical properties required as input by the numerical model were characterized experimentally. The specific heat capacity of the polymer composite was characterized using differential scanning calorimetry (DSC). The shear modulus and viscosity of the material were characterized using dynamic mechanical analysis (DMA) in shear mode using a parallel plate rheometer. The thermal conductivity of the material was characterized using a heat flow meter (HFM).

Differential scanning calorimetry
DSC was used to determine the specific heat capacity of the material at different temperatures using ASTM E1269-11 [53]. DSC was carried out using a TA Q2000 calorimeter. Approximately 5 mg samples were prepared from the feedstock material. A heating cycle followed by a cooling cycle from 25 °C to 220 °C was carried out at a rate of 5 °C/min. Three replicates were tested. The specific heat capacities obtained from the cooling cycle were adopted.

Dynamic mechanical analysis
DMA was carried out to determine the storage modulus and viscosity as defined in ASTM D4021-21 [54] of the PLA/wood thermoplastic composite at different temperatures. DMA of the feedstock material was performed using a Discovery HR 30 from TA instruments. A 25 mm diameter parallel plate set was used. A specimen of 25 mm diameter and 1 mm thickness was prepared and placed in between the parallel plates. The specimen was heated from 75 °C to 220 °C at a rate of 5 °C/min with cyclical strain rates of 10 Hz, 1.0 Hz, and 0.1 Hz.

Thermal conductivity measurement
A Netzch 446 Lambda Series heat flow meter was used to measure the thermal conductivity of the PLA/wood thermoplastic composite. Three specimens with approximately 150 mm × 150 mm width and 16.0 mm height were tested. ASTM C518-21 [55] standard test method was followed for the measurement of thermal conductivity. The top plate was maintained at a temperature 5 °C higher than the mean temperature and the bottom plate was maintained at a temperature 5 °C lower than the mean temperature. Thermal conductivity was measured at seven different temperatures, 25 °C, 35 °C, 45 °C, 55 °C, 65 °C, 75 °C, and 85 °C. The mean of the thermal conductivities at the different temperatures was used as input for the numerical model.

Discrete-event simulation (DES) model
The DES model was used to generate thermal history and the layer deformation in the z-direction due to plastic viscous deformation of the 3D-printed material based on the toolpath and the print parameters. The DES model was previously used to generate thermal history for desktop scale 3D-printed parts and compared with finite-element model results.
The DES model is based on solving finite-difference equations for the heat equation mentioned in Eq. 1 using a forward Euler method. The model discretizes the 3D-printed part into hexahedral (brick) voxels and simulates the deposition of material as discrete events that activate the voxels. The model updates the temperature of the activated element by accounting for heat exchange between adjacent material and the environment in small time increments.
The DES thermal model presented by Bhandari and Lopez-Anido [52] accounts for the loss of thermal energy to the surrounding via convection. The DES model is extended in this research work to account for the heat losses due to thermal radiation using the following equation: where dQ r is the heat loss to the environment via radiation (J ), is the Stefan Boltzmann constant ( 5.670374 × 10 −8 W ⋅ m −2 ⋅ K −4 ), is the emissivity of the material ( adimensional ), A is the area of the surface radiating surface ( m 2 ), T is the absolute temperature at the centroid of the element ( K ), T env is the absolute temperature of the environment ( K ), and dt is the time increment ( s).
This research work extends the DES thermal model to evaluate the local plastic viscous deformation in each layer due to the self-weight and the weight of layers above it. The temperature of each element is calculated at every time step. Based on the temperature, tan δ is calculated for the element. Shear modulus, viscosity, and conductivity at a given temperature are interpolated from the experimental values.
A finite-element mechanical model is used to determine the plastic viscous strains in the elements. 8-noded hexahedral (brick elements) were used. The elements have the same geometry as the hexahedral voxels used in the DES model and are used in the finite-element mechanical model. The finite-element model was implemented in the Rust programming language. An isotropic material model with temperature-dependent Young's modulus and viscosity was adopted. A Poisson's ratio of 0.33 was used for the model. Only self-weight loads were considered in this analysis. The bottom nodes, that contact the bed, were assigned fixed boundary conditions. Since the viscoelastic liquid region of the material behavior affects the stability of the layer, the element was considered active if the tan δ was greater than 1. Once tan δ dropped below 1, the element was considered inactive, and all the nodes associated with the element were assigned a fixed boundary condition. On reheating, if the temperature increased so that tan δ increases above a value of 1, the element was considered active again, and the boundary conditions to the node were removed. Such activating and deactivating of the elements limits the degrees of freedom in the finite-element model and allows the simulation to complete in a reasonable amount of time.
The viscoelastic strain can be divided into two components, plastic viscous strain that increases and accumulates over time, and elastic strain that does not accumulate over time. Since the elastic strain is dependent on the Young's modulus, which is a function of temperature that is changing with time, elastic strain changes with time. However, the elastic strain does not accumulate over time with each time step. The total strain is the sum of the plastic viscous strain and the elastic strain as expressed in the following equation: where e is the elastic strain and p is the plastic viscous strain.
8-noded hexahedral elements were used for the finiteelement mechanical model. The 8-noded hexahedral element can be directly mapped to the elements from the DES thermal model.
The total load on an element increases with the number of layers above the element. The plastic viscous deformation increases with time. The deposited layer deforms in the z-direction under its own weight and the weight of the layers above it. Meanwhile, the layer also cools down, and the shear modulus and the viscosity of the material increase. Thus, the rate of viscoelastic deformation in the z-direction is reduced. If the deformation in the z-direction is too high, the new layer to be deposited does not have a solid base layer under it. As a result, the layer collapses and print failure occurs. Duty et al. [24] used a maximum limit of 10% strain in the layer to be used as a criterion for the collapse of the layer.
Convective heat transfer film coefficient is an important parameter for the DES thermal model. Convective heat transfer film coefficient characterizes the heat loss per unit area at an interface with fluid due to convection. For largeformat extrusion-based additive manufacturing, convective heat losses occur at the interface on the boundary of the 3D-printed part and the air surrounding it. Costa et al. [56] pointed out that for polymer extrusion-based processes, the convective heat losses control the cooling process and heat losses due to radiation are negligible, especially at the low extrusion temperatures. Khalifa [57] reviewed several empirical relations that have been proposed to estimate a convective heat transfer coefficient for isolated vertical and horizontal surfaces. Kunc et al. [13] and D'Amico and Peterson [42] adapted these empirical relations for simulating heat transfer in additive manufacturing processes. A convective heat transfer of 12.0 W/(m 2 K) was used as recommended by Brenken et al. [46,58] based on experimental observations.

Application of the numerical model for process qualification
A 711 mm (28 in.) high 762 mm (30 in.) wide tapering section was analyzed using the numerical model. Figure 3 shows the geometry of the 3D-printed part adopted to (3) = e + p implement the qualification process. The selected part is one of the four segments of a culvert outlet diffuser prototype. The part is a thin-walled structure designed to improve the flow of water through a culvert in a highway. An Ingersoll MasterPrint large-format additive manufacturing equipment was used to manufacture the part. A bead width of 12.7 mm and a layer height of 5.08 mm were chosen for this study. Three different layer times of 2 min 30 s, 3 min 10 s, and 4 min 45 s were used for the study. Two different extrusion temperatures of 200 °C and 207 °C were used for the experimental study. The part has a thickness of 25.4 mm and is 2 beads wide.
Mesh size of 4.23 mm × 4.23 mm × 1.69 mm was chosen so that a single bead (12.7 mm width and 5.07 mm height) is divided into a 3 × 3 × 3 mesh. An emissivity value of 0.92 was based on experiments carried out by Morgan et al. [59] for 3D-printed polymer composites.
The limitations of the assumptions in the numerical model are: The molten material was modeled as solid assuming viscoelastic behavior.  a build failure. Hence, the processing parameters that result in excessively large strains are not used for the process design. 3. Shear locking due to the use of 8-noded linear hexahedral elements with full integration results in higher bending stiffness in the model.
These limitations cause the numerical model to have reduced accuracy. For this reason, we need to demonstrate that the model has adequate accuracy to generate a pass/ fail acceptance criteria regarding the collapse of the layer for a set of printing parameters for the given geometry. The low-fidelity model has the advantage of lower computational demands. Figure 4 shows the different orthographic views of the meshed model and the boundary conditions at the end of the printing process. The nodes, elements, and boundary conditions were exported as an Abaqus ".inp" input file. The Abaqus input file was imported in Abaqus 2020 CAE to display the mesh and the boundary conditions. In the thermal model, the bottom nodes in contact with the print bed are assigned a boundary condition of temperature equal to 25 °C. The surfaces exposed to the environment during the additive manufacturing process lose thermal energy to the environment via convection and radiation. In the mechanical

Material characterization tests
The specific heat capacity as a function of temperature for the semi-crystalline PLA/wood material obtained with DSC is shown in Fig. 5. Negative specific heat values are observed at about 100 °C due to the exothermic peak of the material due to cold crystallization. The specific heat capacity peaks at the melting region of the material which is at around 150 °C. Gracia-Fernández et al. [60] and Shi et al. [61] discussed similar dual exothermic peaks that indicate the melt recrystallization during the processing of the material.
The thermal conductivity of the 3D-printed PLA/wood material was found to be 0.205 W/(m K) with a COV of 8.7%. The storage shear modulus of PLA/wood material during the cooling cycle at a cooling rate of 1 °C/min at an oscillation frequency of 0.1 Hz obtained with DMA is shown in Fig. 6. Figure 7 shows the tan δ for the PLA/wood material obtained from DMA. The tan δ values are used for determining the viscoelastic state of the material at a given temperature. The material is in the viscoelastic liquid region from 220 to 120 °C, in the viscoelastic solid region below 120 °C. This model uses Maxwell-type elements to determine the non-recoverable viscous deformation of the material in the viscoelastic region of the material behavior. The deformation in the viscoelastic solid and the elastic solid region of the material behavior is ignored. The elastic deformation is recoverable and has a negligible effect on the collapsing behavior of the 3D-printed material during molten material deposition. Figure 8 shows the viscosity of PLA/wood material during the cooling cycle at a cooling rate of 1 °C/min measured using DMA in shear mode. The viscosity increases with decreasing temperature, and also with decreasing shear rate. The shear rate of 0.1 Hz is adopted as zero-shear viscosity for this study. The temperature-dependent material properties for this study were chosen for the cooling cycle because the viscous deformation that can lead to the collapse of the  Table 1 summarizes the input material properties used for the coupled thermo-mechanical model. Specific heat capacity, viscosity, and shear modulus are used in the model as temperature-dependent properties in tabular form.
The maximum strain in the z-direction versus the number of elements in the numerical model is shown in Fig. 9. The maximum strain for the entire model in the z-direction ( z,max ) converges with the increasing number of elements. There is a − 0.0362% change in z,max when the number of elements is increased from 311,456 to 1,051,164 compared to − 0.3611% change in z,max when the number of elements is increased from 38,932 to 311,456 elements. Figure 10 shows the results from the simulation for maximum strain in the z-direction versus extrusion rate at different extrusion temperatures. The black cross marks show the simulation results for the temperature and extrusion rates at which the test parts were printed. The red line marks the acceptable unrecoverable plastic viscous strain of 10% as recommended by Duty et al. [24]. The numerical model predicts that the combination of extrusion rate and temperature below the red line will result in print failure due to the collapse of the layers. The model predicts that printing at 36.04 kg/h (143 s layer time) and 27.3 kg/h (190 s layer time) at 200 °C will result in a print failure while printing at 18.2 kg/h (285 s layer time) at 210 °C will prevent the collapse of layers. Figure 11 shows the histogram of plastic viscous strains in the z-direction for all elements at the different printing conditions that were carried out in the verification experiments. Figure 11a shows the histogram of plastic viscous strains for a part manufactured with an extrusion    Figure 11c shows the histogram of plastic viscous strains for a part manufactured with an extrusion temperature of 200 °C and a layer time of 143 s. It can be observed that increasing the layer time reduces the maximum plastic viscous strain in the z-direction. Figure 12a shows the temperature-time history of the first element that exceeds the 10% viscous strain for an extrusion temperature of 200 °C and a layer time of 143 s. Figure 12c shows the temperature-time history of the first element that undergoes the 10% viscous strain for an extrusion temperature of 200 °C and a layer time of 190 s. A small reheating peak is observed in both figures after the layer time has elapsed when a layer is deposited over the element. The rate of cooling also decreases as free surfaces where convection and radiation were previously occurring are now covered by the deposition of layers next to the surface. Figure 12b and d shows the viscous strain accumulation in the first element that exceeds the 10% viscous strain for extrusion temperatures/layer times of 200 °C/143 s and 200 °C/190 s, respectively. Four distinct regions can be observed in the curves. The first region is the viscous deformation of the layer under its self-weight. The material is hotter and has lower viscosity during the initial phase. The material cools down considerably and increases in viscosity and starts to deform at a slower rate with time. The second and third regions are marked by a sudden increase in the viscous deformation rate as a layer of material is deposited over the element in consideration. The deposited layers increase the temperature, thus reducing viscosity. Furthermore, the load over the element increases as well which increases the rate of deformation of the material. The fourth region begins when the material has cooled down to a point where the tan δ of the material becomes less than 1. The material behavior changes to viscoelastic solid and the strains are recoverable as the material cools down. The numerical model implemented does not account for the elastic strains and the accumulated viscous strains remain constant. Figure 13 shows the accumulated viscous strain in the additively manufactured part for different layers. The elements that are deposited later tend to have higher viscous strain in the z-direction compared to the elements that are deposited earlier. The farther away the elements are from the printing bed, the slower they cool down. The printing bed acts as a heat sink, which cools down the elements near it. Consequently, the viscosity of the layers that are deposited later remains higher for an extended period of time, resulting in higher viscous deformation. Figure 13a  in the additively manufactured part. The trend can be seen more clearly with the higher layer having greater viscous strain in the z-direction. The final few layers have a lower strain as they have fewer layers above them and consequently withstand a lower load.

Results of the process qualification for the demonstration part
The print failure of the selected part is illustrated in Fig. 14. Print failure was caused by the collapse of layers due to high plastic viscous strains for two parts printed with different process parameters: (a) 200 °C extrusion temperature, 143 s layer time, and (b) 200 °C extrusion temperature, 190 s layer time. Comparing the two parts, print failure occurs at a higher layer in the part printed at 190 s layer time, compared to the part printed with a layer time of 143 s. In both cases, the observed print failure was consistent with the numerical model prediction of process parameters in Fig. 10. Based on the results from the numerical model predictions of processing parameters (see Fig. 10), the demonstration part was printed with an extrusion rate of 18.2 kg/h and an extrusion temperature of 207 °C. In addition, the toolpath for the 3D-printed part was modified to create a continuous path connecting four segments of the culvert diffuser while maintaining the extrusion rate of 18.2 kg/h and layer time of 285 s. This set of process parameters met the acceptance criteria. Figure 15 shows the part that was successfully manufactured using the processing parameters that met the acceptance criteria.

Conclusions
The following conclusions were drawn from the research work: 1. A simplified, computationally efficient physics-based numerical model was developed to minimize risk in large-format additive manufacturing of thermoplastic composite designs. The model can predict if layer collapse occurs for a given set of printing parameters. 2. The numerical model was used to generate a design envelope of acceptable values of critical printing parameters (extrusion temperature and deposition speed) for a given print job to prevent layer collapse. 3. The simplified physics-based simulation was implemented to prevent excessive plastic viscous deformations in a large-format additively manufactured biobased polymer composite part.  The limitations discussed in Sect. 2.4 cause the numerical model to have low accuracy in predicting strains. The strains calculated using the model are useful for generating the pass/ fail criteria for a given set of printing parameters. However, the strains calculated by the model cannot be used for predicting dimensional accuracy or residual stresses in the additively manufactured part. An opportunity for future research is adapting the numerical model for heterogeneous computing to generate pass/fail criteria results in a relatively short amount of time. A potential direction for future research is implementing a higher fidelity mechanical model to generate accurate strains and predict the dimensional accuracy and residual stresses for additively manufactured parts.