Advances in modeling transport phenomena in material-extrusion additive manufacturing: Coupling momentum, heat, and mass transfer

Material-extrusion (MatEx) additive manufacturing involves layer-by-layer assembly of extruded material onto a printer bed and has found applications in rapid prototyping. Both material and machining limitations lead to poor mechanical properties of printed parts. Such problems may be addressed via an improved understanding of the complex transport processes and multiphysics associated with the MatEx technique. Thereby, this review paper describes the current (last 5 years) state of the art modeling approaches based on momentum, heat and mass transfer that are employed in an effort to achieve this understanding. We describe how specific details regarding polymer chain orientation, viscoelastic behavior, and crystallization are often neglected and demonstrate that there is a key need to couple the transport phenomena. Such a combined modeling approach can expand MatEx applicability to broader application space, thus we present prospective avenues to provide more comprehensive modeling and therefore new insights into enhancing MatEx performance.


Introduction
Additive manufacturing (AM) is a processing technique that employs computer-aided design (CAD) files to fabricate three-dimensional objects using metals, ceramics, polymers, and composite materials [1][2][3][4]. However, filament-based material extrusion (MatEx) of thermoplastic polymers is one of the most common and economical AM techniques due to the easy operability, low material and energy requirements [5][6][7]. MatEx is also a relatively uncomplicated forming method that can be used to manufacture three-dimensional parts with intricate geometric features [8]. MatEx is being widely employed in different industries mainly because it can be used to print functional prototypes like brackets, wire harnesses, ducts, fixtures, and jigs rather than just lab-scale prototypes by carefully selecting the material and processing conditions [9]. Although MatEx has progressed significantly, its complexity demands further research to gain an in-depth understanding of the process and consequently ensure the repeatable production of consistent parts [10].
The typical MatEx-based AM process has been schematically illustrated in Fig. 1 (reprinted with permission from Ref. [11]). In this method, a thermoplastic polymer feed (either in the form of solid pellets or filament) is used to print parts of desired shapes. Briefly, the polymer is melted and extruded through a heated liquefier block (heated section in Fig. 1) via a pinch-roller mechanism (that maintains the feed flow rate). The polymer melt experiences shear rates ~ 100 s −1 inside the nozzle that can have an appreciable effect on the polymer chain orientation and stretch (denoted by Region 1 in Fig. 1) [7,11]. After deposition on the build plate the extrudate cools down and the polymer chains relax (denoted by Region 2 in Fig. 1). Subsequently, the solid filament pushes more material through the nozzle which gets laid on the initially deposited layer where the material undergoes solidification Arit Das and Claire McIlroy are co-first authors and contributed equally.
and weld together to fabricate the final product. The print head and the print bed then moves in the x-y direction (or the z-direction) to continue the print in a layer by layer manner till the entire part has been built. The interactions among the process controls, material properties, process parameters, and printed part properties (mechanical properties, surface roughness, surface finish, etc.) have been reviewed in great detail in previously reported literature [12][13][14].
In summary, the flow of a viscoelastic material through the nozzle, typically at elevated shear rates, coupled with rapid cooling and solidification (and crystallization in the case of some thermoplastics), affects part geometry, porosity, and mechanical strength [12,13]. Ultimately, the part micro-structure determines the buildup of residual stresses and inter-layer (weld) strength. Thus, MatEx involves the simultaneous and co-dependent transfer of momentum, heat, and mass, which are intimately linked to the end-part properties. An overview of these transport phenomena typical to polymer processing, in general, is provided by Bird [15]. A systematic way to analyze MatEx is via multi-scale multi-physics models that can simulate the entire print, accounting for the four following fundamental processes: (a) flow behavior of the polymer melt [15,16], (b) heat transfer [17], (c) viscoelastic properties (which affect anisotropy and inter-diffusion) [18], and (d) solidification and crystallization [19,20].
These processes occur across various length scales; from the size of a polymer molecule (nm) to the scale of the part itself (mm to m).
The modeling approaches to be reviewed here provide powerful tools to explore these complex phenomena occurring across disparate length scales and aim to predict material-specific optimal printing parameters, as well as the mechanical performance of the final printed part. Although there is a substantial amount of current literature that investigate the effect of process variables on the print part properties through various experimental measurements [21][22][23][24][25][26]; such approaches typically focus on empirical data fits rather than a fundamental understanding and model analysis of the key process physics. Therefore, the ultimate goal of such models is to provide tools for quality control during MatEx that respond to the physical phenomena to help guide machine and material design (e.g. changes to heat capacity and thermal conductivity) for improved printing performance. However, deriving equations that link processing parameters directly to part properties is non-trivial.
In general, transport phenomena are described by three governing equations [15]: conservation of mass, conservation of momentum, and the heat equation. That is for an evolving velocity field u(x, t) , stress field (x, t), and temperature field T(x, t) , which vary in both time t and space x , the governing equations are given by where ρ is the density, Du Dt is the material derivative of the flow velocity, p is the pressure, C P is the specific heat capacity, is the thermal diffusivity of the material, and q V is the volumetric heat source term. Analytic solutions to these equations prove challenging even for simple fluids. For more complex materials, such as the polymer melts used in MatEx, the velocity, stress, and temperature fields must also be coupled to an appropriate constitutive equation [27]. Whilst simple relations can be employed to capture sheardependent viscosity, these do not capture the complex (viscoelastic) phenomenon arising from the microstructure of the material. Also, semi-crystalline polymers require further consideration, since the crystallization kinetics [28] affect Fig. 1 Schematic of the MatEx based 3D printing process. The solid polymer filament is fed into the heated section (at a temperature of T N ) at a constant rate (U 0 ) using the pinch roller mechanism. The polymer melt is deposited on the build plate (print bed) at a constant extrusion speed (U N ) which is usually kept equal to the print bed speed (U L ). The flow profile inside the nozzle is highlighted by the dotted lines. Region 1 highlights the melt flow region in the nozzle and subsequent deposition on the print bed while Region II represents the rapid cooling of the polymer extrudate and relaxation of the polymer chains. For a crystalline material, the crystallization process is initiated in Region II (image reprinted from Ref. [11] with permission from Elsevier) 1 3 both the rheology and temperature (via, shear heating and latent heat). Moreover, flow can have interesting implications on the crystallization process [29] and hence affect the ultimate properties of the printed products. Figure 2 provides an overview of the multi-physics associated with MatEx and further illustrates how improper control of one or more of these processes can lead to failure of the MatEx printing process. A full computational model of the MatEx process would require evolving Eqs. (1)-(3) forward in time across the spatial domain (which is also time-dependent due to the layer-by layer process), together with full non-linear coupling to a constitutive equation that is molecularly-aware and accounts for the multi-phase nature of semi-crystalline polymers. Despite the increase in computational power, this is an unwieldy task. Instead, a fundamental understanding of MatEx may be gained by simplifying the problem in various ways (for example, considering steady-state behavior, reducing the spatial dimension, and decoupling the processes) In this paper, we review recent advances (over the last 5 years) in modeling these complex transport phenomena for applications in MatEx based AM. Table 1 highlights the current state of the art in modeling approaches for MatEx with reference to the fundamental processes modeled as well as the scale at which the simulations were performed.
We discuss models that consider material flow to predict the shape of the deposited layers in Sect. 2. Heat-transfer models that predict the thermal profiles over multiple length scales (from inside the nozzle to the final part) are discussed in Sect. 3. In Sect. 4, we review how these models can be coupled with classical polymer constitutive analyses to provide valuable insights into the micro-structural evolution and ultimate strength of a printed part after solidification. Classical polymer crystallization models can also be incorporated where appropriate to reveal details of the crystal content and morphology, which depend on both thermal and flow history. Finally, we address the challenges associated with these approaches, potential methodologies to address these challenges, and recommendations on developing comprehensive models to accurately represent the intricate transport phenomena occurring during the MatEx process. In traditional polymer processing techniques (e.g. injection molding, extrusion etc.) where both shear and elongational flow of the polymer melt are of critical importance [15,16], geometry of the extrudate is not usually significant since well-defined molds typically dictate the final geometry of the manufactured parts. In contrast, polymer flow through the nozzle and subsequent deposition determines the geometry of extruded roads in MatEx, and thus the fully assembled shape of the part itself [30,31]. The road geometry depends on both the print resolution and print speed [32] and is ultimately a result of the polymer's viscoelastic properties. Agassant et al. [33] highlight the relationship between print head velocity and pressure at the nozzle outlet. This relationship leads to an optimal operation window ( Fig. 3a, Table 1 Summary of the modeling approaches in material extrusion based additive manufacturing along with the polymer materials that were used to validate the models The checkmarks indicate the fundamental models that were analyzed in the corresponding studies. The terms within parenthesis in the Materials column indicate the scale (road/part) of the simulation performed. It is worth mentioning that the term "road" has been used to denote both the extrudate exiting the nozzle and the material that is deposited on the printer bed

Fig. 3
a Phase map of the pressure across the printer nozzle as a function of print velocity showing the risks associated that can lead to print failure (image reprinted from Ref. [33] with permission from Elsevier); b cross-sectional shapes of the extruded filaments from CFD simulations and as observed experimentally using optical microscopy using PLA (image reprinted Ref. [32] with permission from Elsevier); c SEM image of a cross-section of a single layer of print using polyetherimide (PEI) (image reprinted with permission from Ref. [40]. Copyright 2019 American Chemical Society); d experimental 3D printed lattices and simulations of the lattice models as generated by the VOLCO model as a function of layer thickness. The filament spreading as a function of layer thickness was accurately captured by the VOLCO model (image reprinted from Ref. [34] under the terms of the Creative Commons CC-BY license) reprinted with permission from Ref. [33]) that eliminates the risk of buckling (as illustrated in Fig. 2), ensures continuous road deposition, and provides sufficient pressure to encourage inter-layer diffusion and subsequent adhesion. The shape of the roads also provides information regarding porosity [34], inter-road bonding [35], and thermal history [36]. The effect of the print parameters on the road geometry [32,33] and surface roughness [37] of a printed part has also been investigated. The print speed, extrusion flow rate, and, gap between the nozzle head and print bed were found to have a pronounced effect on the shape of the deposited roads as well as the mechanism of extrude spreading on the print bed (hydrodynamically controlled vs. surface tension controlled) [32,33]. On the other hand, layer thickness was found to have a significant effect in determining the surface roughness of a printed part [37]. Bakrani Balani et al. [38] studied the effect of the printing parameters, rheological and thermal characteristics on the geometry of deposited roads of polylactic acid (PLA) during MatEx. The authors further developed a two-phase flow model that was employed to model the geometry of the deposited roads [38]. The simulation results (suggesting deformation of the extrudate during deposition on the printer bed) were in close agreement with the experimental findings using an optical microscope, demonstrating the robustness of the modeling approach [38]. Simple analytic models based on conservation of mass have been developed to predict the height, width and crosssectional area of a single deposited road [30,31]. While this approach gives a reasonable agreement with CT-scanning images, it is unable to provide a full description of the ovality of the road, which impacts the weld contact area between adjacent roads (in both the x and z directions), or the part shape. Thus, a number of full computational fluid dynamics (CFD) models have recently been presented [31,39] and validated against optical microscopy experiments [32], as well as SEM images [40], as shown in Fig. 3b, c (reprinted with permission from Refs. [32,40]), respectively. These analyses illustrate a range of road shapes, from circular to elongated rectangle, depending of print parameters (Fig. 3b, reprinted with permission from Ref. [32]). A new computationally-efficient volume conserving (VOLCO) model has also been developed [34] to capture the geometry at roadroad interfaces and the porosity within a full 3D lattice. The flow model described was able to represent the printing of a full-size 3D lattice, rather than just a single filament and successfully described the experimental observations (Fig. 3d, reprinted with permission from Ref. [34]).
Other important geometrical aspects that have recently been considered via numeral modeling include surface instabilities along the road [38] and sagging of freely-supported roads [41]. While this progress has significantly enhanced our ability to simulate geometry, temperature-dependent viscoelasticity and crystallization must be coupled to these models to fully understand the complex flow dynamics. While power-law models have provided some insight into the spreading dynamics [33], a constitutive model that incorporates normal stress differences is necessary to capture extrudate swell effects [42]. Moreover, density changes during cooling, which are exaggerated during crystallization, are responsible for part shrinkage and warping (as highlighted in Fig. 2 as well).

Heat transfer
Melt-based processing of polymers is associated with complex interactions between fluid mechanics and heat transfer in a non-Newtonian fluid system [17,43]. Similarly, in the MatEx process the filament must be molten to flow through the nozzle and hence heat transfer significantly impacts the print speed [44]. The thermal profile of the polymer from nozzle through to deposition (Fig. 4a, reprinted with permission from Ref. [45]) also controls inter-road adhesion (also known as welding) through polymer diffusion, significantly impacting final part properties [40,[46][47][48][49]. Thus, the ability to accurately predict heat transfer is crucial to fully understanding MatEx.
Transient one-dimensional heat-transfer models considering conduction and natural convection, which can be solved via finite-difference methods, have been implemented to calculate inter-road bond strength of acrylonitrile butadiene styrene (ABS) with reasonable accuracy [50]. A similar approach, with the addition of radiation, has also been developed for big-area AM (BAAM), and is able to predict the temperature profiles and thermal gradients in thin walls which was validated by infrared imaging of carbon fiber (CF) reinforced ABS prints [51]. A 2-D approach with order-of-magnitude analysis was adopted by Osswald et al. [52] for a shear-thinning quasi-Newtonian fluid, which was validated using a custom piston-driven capillary setup. However, these models are insufficient to properly represent the temperature gradient within the deposited roads [51,52]. Accurate prediction of this thermal gradient is key to determining the cooling rate after deposition, and hence the evolution of weld strength in the printed parts [40].
With the improvement in computational capabilities, 3-D heat-transfer models for MatEx, which can be solved using finite-element (FE) methods, have been investigated [36,53]. For example, a simulated thermal profile within the nozzle head is shown in Fig. 4b (reprinted with permission from Ref. [54]) [54]. FE models have also been developed to capture rapid cooling during the deposition process [55], however assumptions regarding interlayer bonding resulted in predicted trends which overestimate the heat transfer measured using infrared imaging techniques. The entire MatEx process (die extrusion, road deposition, and subsequent cooling and consolidation) has also been modeled recently; the thermal profile of an extruded road from such an analysis is illustrated in Fig. 4c (reprinted with permission from Ref. [40]). Zhang et al. [56] developed a space resolved, time-dependent 3-D model to capture the thermal behavior during and after MatEx. The integrity of the welds formed between the layers was found to strongly dependent on the printing speed while residual stresses inside a print could be controlled by manipulating the cooling rate [56]. The thermal history associated with the extrudate and the already deposited layers during printing was found to be dependent on the print temperature, print speed, layer thickness, and environmental temperature [56]. Zhou et al. [57] also developed a 3-D FE thermal model to describe temperature and stress distribution in a printed part of PLA. The model predicted the extruded filament to remain above the glass transition temperature (T g ) for ~ 0.1s (time scale during which inter-layer bonding occurs) which was validated experimentally. The simulation studies revealed that lower print temperature, reduced print speed, and lower layer thickness can mitigate the residual thermal stress and part distortion. Further examples of 3-D FE heat-transfer models can be found in Refs. [58,59]. The models discussed in this section are only valid for materials with isotropic properties and their application is severely limited for fiber reinforced composite materials with anisotropic thermal conductivities. Moreover, thermal modeling for the multi-material MatEx process is still in its nascent stage. Arguably more important is the validation of analytical and numerical models via comparison to experimental measurements. Luo et al. [60] recently developed a heat transfer model based on previous modeling efforts on polymer melt flow in tubes [61] and 3D printer nozzles [62]. The thermal model provides an upper limit of the polymer feed rates that can be successfully printed using MatEx without jamming [60]. The authors validated their simulation results by performing experiments on a modified extruder setup and found that their observations were in good agreement with model predictions, especially at a nozzle wall temperature range of 210-260 °C [60]. A non-isothermal welding model was developed by Coasey et al. [63] and experimentally verified by interlaminar fracture toughness measurements. The model describes the effect of processing parameters and rheological properties on the degree of healing of the printed parts [63]. The melting behavior of PLA inside a MatEx nozzle was computationally modeled and validated using experiments [64]. The authors coupled a generalized Newtonian fluid model with the Phan-Thein-Tanner viscoelastic model to model the melting behavior and subsequent melt flow in the nozzle under high shear rates [64]. Edwards and Mackay [65] also developed a post-extrusion heating model for both amorphous and semi-crystalline polymers to minimize the formation of "shark-skin" defects during MatEx.
Key works that account for radiative heat transfer effects have shown excellent agreement between the temperature profile in the nozzle prior to and after deposition [66], in the standoff region (between nozzle and bed) [67], and on the print bed [45].
Despite this progress, it is apparent that the rheological properties of the melt significantly affect the heat transfer by influencing the flow characteristics (e.g. viscous heating) [55]. Hence, to completely understand the complexities of the heat-transfer process, appropriate polymer constitutive models must be considered. Mackey et al. analyzed the problem in the nozzle (where melting and extrusion occurs) using a generalized Newtonian fluid framework to model the rheological behavior [68], reporting that heat-transfer affects the pressure drop across the nozzle and limits print speed. On the scale of a final printed part, Brenken et al. coupled a 3-D heat-conduction model with temperature-dependent viscoelasticity to simulate the transient thermal profile in a functional part as highlighted in Fig. 4d (reprinted with permission from Ref. [69]). Moreover, polymer crystallization and material shrinkage are also incorporated into their model that are crucial ingredients for predicting latent heat effects and part geometry.

Anisotropy and inter-diffusion
Viscoelasticity plays a crucial role in polymer processing [70] and hence it is of considerable importance in MatEx as well [68,71]. It has been shown that typical shear rates are sufficient to orient and stretch the polymer molecules [72]. Thus, as well as affecting heat transfer in the nozzle via shear heating [73], and influencing the geometry of a deposited road via extrudate swell [42], polymer deformation (anisotropy) may persist at solidification. For an amorphous polymer, Fig. 5a illustrates the typical degree of residual alignment across a single road locked in due to the glass transition, as predicted by the numerical model of McIlroy et al. [72].
Moreover, the key to ensuring part strength is successful welding of the polymer chains at road-road interfaces via inter-diffusion [74]. The poor interlayer bonding between the printed layers is highlighted in Fig. 2. As a rule of thumb, if there is sufficient time above the glass-transition temperature, T g , for the polymer to diffuse its radius of gyration (R g ), then bulk strength should be recovered in the weld region [75]. Figure 5b (reprinted with permission from Ref. [55]) demonstrates a typical cooling profile found in MatEx; for these print conditions (and thermal properties) there is approximately one second above T g available for interdiffusion [55]. Polymer alignment leads to anisotropic diffusion (in the direction of chain orientation) [76] and may hinder this welding process, although this effect is found to be small for polycarbonate (PC) [72] and for PLA as shown in Fig. 5c [72].
The role of polymer entanglements is also crucial to ensuring strength at road-road interfaces [74]. It is proposed that MatEx can lead to flow-induced disentanglement of the polymer network [77] and that only partial re-entanglement of the network is achieved during cooling [72]. McIlroy et al. [72] suggest that it is this partial entanglement, rather than anisotropic diffusion, that is responsible for weld weakness in amorphous parts. However, this is yet to be validated experimentally. Nevertheless, residual viscoelastic stresses affect part properties, in particular reducing the tensile properties of printed parts (as highlighted in Fig. 5d, reprinted with permission from Ref. [78]).
Since power-law models are not sufficient to capture this molecular-scale behaviour, recent works by McIlroy et al. have employed a molecularly-aware constitutive model [11,72,79,80], which requires knowledge of both the reptation and Rouse times of the polymer (obtained from linear rheology). In particular, McIlroy and Olmsted [80] prescribe the geometry of the deposited road, which enables the flow field to be solved directly from conservation of mass independently of the constitutive (Rolie-Poly) equation using finite-difference methods. The authors find that the melt experiences significant strain rates within the nozzle, which are able to substantially stretch and orient the polymers. Furthermore, the deposition process involving the 90° turn dominates the deformation and significantly disentangles the melt. Although this approach [80] provides a computationally efficient method to investigate a range of printing conditions, the interplay between viscoelasticity and the weld geometry is neglected. Since the Rolie-Poly model provides a simple one-mode constitutive equation for the stress tensor, it can be readily integrated into full CFD calculations to capture this behavior. Recent advances by Xia et al. [39] incorporate a reduced Rolie-Poly model into a full finite-element computation of the flow field to investigate the effect of these viscoelastic stresses on the deposition geometry. Fig. 5 a Residual alignment in a single deposited road of PLA printed at 120 mm/s and 190 °C, as predicted by the numerical model described in Ref. [72]; b prediction of the thermal profiles of the deposited layers using FE simulations; the different layers spend varying time periods above T g which determines the extent of their mobility (image reprinted from Ref. [55] with permission from Elsevier); c the polymer diffusion distance, χ, scaled by the polymer size, R g (radius of gyration), as a function of time. Time scaled by the time taken to reach the glass transition temperature, so that at t/t g = 1, T = T g . The effects of polymer alignment in the flow give rise to anisotropic diffusion, which hinders diffusion between filaments, i.e., in the transverse z-direction. See Ref. [72] for details of the model; d reduction in tensile properties of printed ABS and polyether-etherketone (PEEK) samples compared to injection molded parts (image reprinted from Ref. [78] under the Creative Commons Attribution License) 1 3

Polymer crystallization
The fundamental mechanisms of polymer crystallization are detailed in a review article elsewhere [81]. The use of semi-crystalline polymers in MatEx has sparked significant interest but is largely relegated to academic studies [82][83][84][85][86]. Briefly, the crystalline regions in these polymers act as physical crosslinks that restrict the movement of the polymer chains, thereby strengthening the material [19,20]. Thus, crystallinity is a major factor in determining tensile properties, and ultimately dictates the application of printed parts. For example, Fig. 6a (data adapted from Ref. [87]) highlights the variation in tensile strength of a common engineering thermoplastic, polyphenylene sulfide (PPS), with the percent crystallinity of the printed parts. The tensile strength increases by ~ 37% with increasing crystalline content due to variations in printing conditions. The most common approach to modeling polymer crystallization considers the formation of nuclei and subsequent growth into a stable crystalline framework. However, traditional isothermal models fail to capture the complex effects of the steep temperature gradients and cooling rates found in MatEx. As noted in the heat transfer section, Brenken et al. [69] tracked the thermal profile and crystallization process during MatEx using a dual (Velisaris) non-isothermal model which has been summarized in Fig. 4d (reprinted with permission from Ref. [69]). Further to this work, a finite-difference model has been developed to account for Fig. 6 a Effect of changes in percent crystallinity of printed part on the tensile strength of polyphenylene sulfide (data adapted from Ref. [87] with permission from Elsevier); b flow-induced crystallization plot (image reprinted from Ref. [11] with permission from Elsevier); c optical microscopy of annealed PLA printed cross-sections. Darker regions in the weld correspond to smaller sized spherulites due to flow-enhanced crystallization. The scale bar is 200 µm. (image reprinted with permission from Ref. [79]. Copyright 2019 American Chemical Society); d improvement in tensile and flexural properties post-annealing for 3D printed PLA and PLA/cellulose nanofiber (CNF) composites (data adapted from Ref. [101] with permission from Elsevier) multiple crystallization mechanisms, shear stress evolution, and more complicated thermal histories [88]. The model predicts both the thermal distribution and evolution of the crystal fraction throughout a printed part which are then applied to novel models for predicting the degree of healing and residual stress at each layer. This is the first healing model that accounts for crystallinity and the experimentally determined bond strength was found to be close to the theoretical predictions.
In the models discussed thus far, the polymers are assumed to crystallize quiescently, i.e. the effect of flow is neglected. However, it is well known that even slight alignment of the polymers at the onset of nucleation can strongly enhance the nucleation rate [89], and consequently lead to accelerated crystallization times. McIlroy et al. [11,79] have recently extended their finite-difference modelling approach (incorporating the Schneider rate equations) to show that there exists a MatEx-operation window in which crystallization times are reduced due to shear in the nozzle (Fig. 6b, reprinted with permission from Ref. [11]). Furthermore, it is shown that this flowenhanced crystallization is localized to a thin boundary layer near the road surface, leading to the development of smaller crystal structures in the weld region (Fig. 6c, reprinted with permission from Ref. [79]). These results are in agreement with in-situ Raman spectroscopy and optical microscopy measurements.
Shmueli et al. [90] coupled in-situ wide-angle X-ray scattering with infrared imaging during MatEx of PLA filaments to study the morphology development. The authors found that thermal profile as well as the printing direction (with respect to the X-ray beam) significantly affected the structure of the deposited roads [90]. Subsequent work performed by the same group [91] using polypropylene (PP) involved connecting a MatEx printer with a X-ray scattering beamline. The authors observed the formation of shish-kebab crystal structure post printing, which nucleated at the surface of the filaments and gradually advanced toward the center of the filament [91]. The experimental method provided both time and spacedependent data for the crystallization process and can be a valuable resource for capturing the detailed kinetics as well as the microstructure development during a print. Welding at the interface of two successive printed layers was found to occur due to polymer chain relaxation and diffusion [91]. In a more recent study, Nogales et al. [92] monitored the microstructure development and weld formation during MatEx of PP by performing in-situ X-ray scattering experiments in the synchrotron beamline. The crystallization mechanism during the printing of the layers was found to be similar to that encountered during quiescent crystallization of PP. Moreover, the printed layer had increased crystalline content in the bulk with respect to the interfacial regions indicating a gradient in crystallinity across a single layer [92]. The authors suggested that this can improve interlayer welding and hence improve the mechanical properties of MatEx printed parts since interlayer welding is directly related to the extent of polymer chain diffusion between the layers [92]. Whilst inter-diffusion must precede the crystallization to ensure sufficient inter-road adhesion, the development of smaller crystal structures may be considered beneficial to improving weld properties [79]. Moreover, McIlroy et al. [79,80] propose a simple comparison of time scales, which can be applied to any printing material, so that print conditions can be optimized to ensure (or eliminate) flow-enhanced crystallization.

Discussion and future directions
Numerous limitations exist in 3D printing technology, as discussed in Ref. [93]; multi-scale molecular modeling approaches are critical in addressing challenges in material diversity and print speed. A useful tool to evaluate the 'printability' of a range of materials, including high-temperature thermoplastics, fiber-reinforced thermoplastics, and low-viscosity thermosets has been presented by Duty et al. [94], and a range of approaches have been tested for improving part properties: • Since molecular weight can significantly affect the welding dynamics [33], tuning the molecular-weight distribution may be used as an avenue for optimization. In particular, Levenhagen et al. [95] show that a bimodal distribution can improve mechanical properties and isotropy of printed parts; the addition of a smallmolecular weight tail enhances the inter-diffusion at the welds. Current printing models consider only a single mode description of the molecular-weight distribution, however, new capabilities are emerging [96], offering exciting opportunities for a more detailed modeling approach. • The addition of particulate additives can also improve part properties. In particular, Brenken et al. [97] and Fallon et al. [98] review a range of different composite systems, comprising of discontinuous and continuous fiber-reinforced polymers. The materials are characterized by the fiber alignment, which is oriented parallel to the direction of the deposited road, thus can lead mechanical limitations in the transverse direction. Some progress has been made in modeling the material flow and orientation in these systems [99], but many open questions remain including wetting properties, bond formation, and thermo-mechanical solidification behavior.

3
• There is also precedent for the employment of postprinting annealing to remove residual viscoelastic stresses and increase the final crystal fraction, thereby improving thermal [100] and mechanical properties (as highlighted in Fig. 6d, data adapted from Ref. [101]) [102]. However, it has been shown that annealing may not result in the desired uniformity of spherulite sizes [79]. In particular, the model of McIlroy et al. [79] shows that a higher nucleation density may be 'templated' into the road by the flow during printing. Thus, on annealing, smaller flow-enhanced spherulites are revealed in the weld regions, as seen by optical microscopy (Fig. 6c, reprinted with permission from Ref. [79]). Nevertheless, it is proposed that these smaller spherulites can give rise to a more ductile fracture.
Although great progress has been in made understanding MatEx using the modeling approaches discussed, many challenges remain; foremost being that coupling viscoelastic flow, heat transfer, and crystallization kinetics are computationally demanding [39,40]. Table 1 provides a list of recent literature on modeling the MatEx process. As evident from the data, there are very few works that take into account all four of the fundamental processes that occur simultaneously during the printing process. Future research must address this concern even though it requires substantial computational power. Current efforts have investigated only a limited material catalog (PLA, ABS, and their fiber reinforced versions) for validating the model predictions; we expect research with other polymer materials (such as, high temperature engineering thermoplastics) to surge in the recent future. Another important observation from Table 1 is the amount of research on modeling the extrudate exiting the printer nozzle or on the deposition of roads on the printer bed which comfortably outnumber the research works that deal with simulating actual functional parts. It is expected that this trend will start to change in the near future with more industries and academic groups investigating large area AM [103][104][105][106] and bridge the gap between the current simulation scale (~ μm-mm) and part dimensions (~ m).
Quantitative predictions also rely on thorough materials characterization, including linear viscoelasticity, nucleation and crystal growth rates, as well as thermal properties, which remains non-trivial since the thermoplastic polymer properties for elevated temperature processing in MatEx are usually highly temperature dependent [107]. The Additive Manufacturing Benchmark (AM-Bench), focused primarily on MatEx and selective laser sintering, has been created as a platform for performing extensive in-situ and ex-situ to generate a huge repository of experimental data. The goal of the AM-Bench is to gain valuable insights into the current modeling endeavors by comparing the simulation results against the rigorously controlled benchmark data [108]. One of the current limitations of the simulation capabilities is the failure to describe spatial evolution of microstructure in a printed part but significant progress have been made in modeling the same for smaller build volumes. It is expected that incorporating the multiphysics described in this paper with the mechanics (e.g. residual strain/stress) of the printed part and optimizing the computational demands will lead to the development of models that can be of great industrial relevance [108].
Moreover, the profound effects of flow on polymer crystallization are still not fully understood [89]-in fact, even assumptions must be made on the nature of quiescent nucleation-and continuum-level models that are able to incorporate the full complexity of the crystallization process remain in their infancy. Finally, there is experimental evidence of a phenomenon known as co-crystallization [109][110][111], which results in improved mechanical properties of a part, but models that account for this mechanism during MatEx are yet to be developed.