Application of Computational Fluid Dynamics (CFD) Simulation for the Effective Design of Food 3D Printing (A Review)

The progress of food 3D printing (3DP) applications demands a full understanding of the printing behavior of food materials. Computational fluid dynamics (CFD) simulation can help determine the optimum processing conditions for food 3DP such as layer height, deposit thickness, volume flow rate, and nozzle shape and diameter under varied material properties. This paper mainly discusses the application of CFD simulation for three core processes associated with 3DP: (1) flow fields in the nozzle during the extrusion process; (2) die swelling of materials at the die (the exit part of the nozzle); and (3) the residual stress of printed products. The major achievements of CFD simulation in food 3DP with varied food materials are discussed in detail. In addition, the problems and potential solutions that modelers encountered when utilizing CFD in food 3DP were explored.


Introduction
Three-dimensional printing (3DP) is a rapid prototyping technology that combines material science with computing and control technologies [1]. The application of 3DP to food began primarily with digitally controlled extrusion operations that could accurately build up the material deposition layer by layer from a pre-designed file such as stereolithography (STL) or computer-aided design (CAD) [2,3]. Several studies have emerged with strong efforts in applying 3DP technology to various designs of food products. In the early stage, the material used in food 3DP was limited to chocolate products [4][5][6], but recently, various materials have been applied successfully and their characteristics in the 3DP have been reported (i.e., potato starch with pea protein, wheat, corn starch, lemon juice gel, alginate with pea protein, cereal food, oat protein, meat product, and surimi [7][8][9][10][11]). 3DP applications were accompanied by a large number of trial-and-error attempts and conservative designs, which resulted in a waste of resources and time during the experiments. This is in part due to a lack of accurate instruments for predicting the results of the process [9]. Recently, patents associated with food 3DP technology have presented computational methods, mainly computational fluid dynamic (CFD) simulation, to predict and quantify the performance of 3DP [12][13][14].
CFD is one of the computer simulation methods used to characterize fluid flow behavior under certain geometries with boundary conditions [12,14]. The criterion for its performance is how well numerical simulation outcomes align with experimental results conducted under specific conditions, and how well simulations can predict extremely complicated processes that cannot be analyzed in the real world [9]. CFD has gained recognition as an emerging science around the world, particularly since the late 2000s when

Overview of Studies on the CFD Analysis of 3DP
Over the last 20 years, a total of 44 research articles on CFD research and its applications in 3D food printing and extrusion analysis have been reported when searched using specific keywords (CFD, food, numerical simulation, 3D printing, and extrusion process) on Google Scholar. The annual distribution of findings reflects the growing importance of the topic, with more than half of the publications occurring in the last four years (Figure 1). The majority of them focused on the development of novel food items. Several researchers have recently begun exploring issues related to 3D food printing processes such as the die swell effect at the 3D printer nozzle, residual stress component in the extruded products, and the shape of printed products during and after the 3DP process. However, such information is disseminated in several publications with differing technical foci. The related articles published during the last 10 years were analyzed (Table 1).  To predict the residual stress, total deformation, and printing precision requirement.
Transient [9] 3D printing of cereal grains. ANSYS POLYFLOW To investigate the use of cereal grains as alternative printing material. Transient [28] 3D printing of lemon juice gel ANSYS POLYFLOW Evaluation of the fluid characteristics of lemon juice gel in the 3D printer flow channel.
Steady [1] Comparative study of a syringe and screw-based 3D food printer

COMSOL Multiphysics
Developing CFD models to investigate and compare the flow field characteristic of two 3D extrusion-based printing units (syringe and Screw based printing).
Transient [21] Soy white flakes based dough in a single screw extruder

ANSYS POLYFLOW
To characterize the rheological properties of soy white flakes-based dough in a single screw extruder and analyze the flow behavior in the high shear rate zone extruder.
Transient [29] Extrusion process of soy white flakes Design Expert To analyze the extrusion process of Jatropha seeds by a single screw extruder.
Steady [30] 3D deposition of porous scaffolds for cartilage tissue engineering To establish the flow of biomaterials through the nozzle. [22] Mixing in food extrusion

ANSYS POLYFLOW and ANSYS Fluent
To characterize the dispersive mixing of oil in the plasticized starch matrix within a twin-screw extruder.

Fundamentals of the CFD Extrusion Process
The extrusion process is the core process for the food 3DP [3]. The extrusion in the food 3DP is a thermomechanical technology that was specifically designed for mass production [1,[9][10][11]. A correct combination of variables involving product formulation, process, and equipment settings are required to obtain the final product. There is a sizeable number of process variables involved in the extrusion process. Thus, it causes difficulty in controlling and analyzing the extrusion process theoretically. In extrusion processes involving either screws or pistons, materials are compressed and sheared with extremely high pressure, which causes thermal and mechanical stresses on the material [9,21]. The thermomechanical effects in the screw or piston section have a considerable influence on the printing process and, as a result, play an important role in the structural and mechanical characteristics of the extruded material. Strong interactions between mass, energy, and momentum transfer, coupled with complex physicochemical transformations governing thermochemical properties are present, making control of the extrusion process difficult. Therefore, understanding the mechanisms involved in the process is an essential step to apply the CFD simulation in the design of 3DP products or 3D printers for food.

Governing Equations
The mathematical equations of fluid motion have been developed for approximately two centuries. First, the Euler equations, which describe the motion of fluid based on the conservation of energy, mass, and momentum, were formulated in 1756-1757 [31]. The Navier-Stokes equations describe fluid motions based on the stress tensor with the Euler equations [32,33]. The Navier-Stokes equations are mathematical formulations of the conservation laws of fluid dynamics (CFD), which are referred to as the governing equations of fluid flow and heat transfer. When applied to a fluid continuum, these conservation laws attribute the rate of change of a desired fluid property to external forces and may be thought of as follows [27,[34][35][36]: i.
The law of conservation of mass (continuity), which states that the mass flows entering and exiting a fluid product must be equal (Equation (1)); ∂ρ ∂t where ρ and t are density (kg/m 3 ) and time (s); x is the Cartesian coordinates (m); u is the velocity component (m/s); and i is the Cartesian coordinate index. ii. The law of conservation of momentum (Newton's second law of motion), which states that the number of external forces acting on a fluid particle equals the rate of change in linear momentum (Equation (2)); where j is the Cartesian coordinate index; δ is the Kronecker delta; µ is the dynamic viscosity (kg/ms); and g is the acceleration due to gravity (m/s 2 ). iii. The law of conservation of energy (the first law of thermodynamics), which states that the rate of change of energy of a fluid particle is equal to the heat addition and the work done on the particle (Equation (3)); where C a is specific heat capacity (W/kgK); T is the temperature (K); λ is the thermal conductivity (W/mK); and s T is the thermal source W/m 3 [36].
By applying these conservation laws over discrete spatial volumes in a fluid domain, it is possible to achieve a systematic account of the changes in mass, momentum, and energy as the flow crosses volume boundaries. Although direct solving of Navier-Stokes equations is possible for laminar flow, computationally solving fluid motion in the Kolmogorov microscales in turbulent flow is not yet possible; hence, turbulence models need to be solved in addition to the Navier-Stokes equations in the case of turbulent flow regime, which is generally encountered in the modeling of a low viscosity food process because of the high flow rates and complex geometry involved [37]. As a result, choosing a suitable turbulence model is crucial and has a direct impact on the CFD results of the 3DP process. Several review studies are available on selecting the suitable turbulence models for various food processing, comparable to the 3D printing process, based on accuracy, computing time, and cost, thus this review will not go into detail [25,37].
Specifically in the extrusion process in 3DP, when developing a thermomechanical finite-element (FE) model for the stress and deformation analysis involved in the 3DP, it should be noted that the calculation of 3DP material deformations and residual stresses is dependent on the temperature field variation and the phase transfer according to the material's mechanical properties such as plastic strain and plasticity [38]. Thus, the initial mechanical properties and boundary conditions are required to develop the structuralmechanical problems. To simulate the heat effect during the 3DP process, the total heat flow introduced into the component must be estimated. This can be computed using the 3DP material deposition parameters, taking into account the thermal efficiency and radial dispersion of the material's heat flux density. As a result, as indicated in Equation (4), the temperature-field distribution may be characterized by the transient and nonlinear heat-balance equations of the Fourier-Kirchhoff heat conduction with the thermophysical characteristics, namely thermal conductivity (λ), specific heat capacity (c), and density (ρ) being dependent on the temperature and 3DP material sample [9,38].
If the convection is ignored for a moving heat source with a constant printing speed (v w ) and a torch power (Q T ), Equation (5) is expressed in a moving coordinate system: For boundary conditions used in simulations, the various temperature loading conditions such as uniform, time-dependent, or spatially varying temperature throughout the selected geometry are used as boundary conditions to calculate the heat losses due to conduction and convection heat transfer associated with the 3DP process. The heat source for 3DP on a localized distribution spot is usually assumed to be an ellipsoid-moving heat source with a normally distributed heat flux density in the form of a Gaussian function ( Figure 2). The localized heat distribution strongly affects the 3DP process as well as the 3DP products. Thus, temperature-dependent heat losses around the geometry must be analyzed. Radiation and convection heat transfer are the major modes for the localized heat distribution for all boundaries in the geometry. In this case, the heat transfer can be modeled by Equation (6), where α k is the heat transfer coefficient due to free convection; ε is the emission coefficient; C 0 is the radiation coefficient; T is the component temperature; and T 0 is the ambient temperature. Radiation and convection are significant only on the upper part of the geometry, but relatively less on the lower part. Thus, ε can be assumed to be 0. The local heating of the component causes an uneven thermal expansion, whereby the colder environment hinders the expansion of the warm areas. This leads to the component-stress formation in the material, which results in deformation [9]. The 3DP products usually deform as soon as the cooling process is completed [39]. With many food materials, a linear-elastic material behavior only applies in a limited initial range of the deformation or mechanical load. The deformation of the 3DP product caused by mechanical loading during 3DP is reversible in many cases. This means that no deformations remain after the material has been relieved or solidified [9]. However, as a result of high mechanical loads, the stress-strain behavior becomes non-linear. If the stress exceeds the yield stress, the material is plastically deformed. This plastic deformation reduces the stress and appears after the relief as a permanent deformation in the 3DP product. To determine the plastic strain, it is necessary to combine the flow condition, the flow law, and the solidification process to determine the load limit of the material used [38]. Model deformation during 3DP of the mixture of alginate and pea protein was induced when the stress values were above the yield stress. However, a quick adjustment in environmental temperature during printing to a level that leads to solidification of the material was reported to result in a substantial reduction in product deformation [9]. A detailed overview of 3DP product deformation is presented in Section 4.1.

Numerical Analysis
The selection of appropriate approaches to discretize the modeled fluid continuum is a critical factor for the application of CFD. The most prominent of these include finite difference (FD), finite elements (FE), and finite volumes (FV). Due to difficulties in the handling of complex geometry, FD techniques are of limited use in engineering flows. FE has been used to successfully model the piston-based 3DP of the mixture of alginate and pea protein [9] as well as the screw-based 3D extrusion process of potato paste [21]. However, the FE requires special care to ensure a conservative solution. Although the FE method must be carefully formulated to be conservative, it is far more stable than the FD method [40]. The FE method may demand more memory and may require more time to obtain solutions than the FV method [34]. Because of the complexities inherent in programming and implementing this technique, there are few commercial FE software available. Fortunately, such challenges may be avoided by employing FV approaches [41] that provide ease in understanding, programming, memory usage, solution speed, and versatility of finite volumes, particularly for large problems, high Reynolds number turbulent flows, and source term dominated flows. It is now the most commonly used numerical technique in CFD code development, and as a result, it is considered the most commonly used method to form design solutions for the 3DP process. In the FV method, the governing partial differential equations such as the Navier-Stokes equations, the mass and energy conservation equations, and the turbulence equations are developed in a conservative form, and then solved over discrete control volumes.

Generation of the Computational Grid (Meshing Process)
In CFD simulations for 3DP, the product geometry, often called a model, must be built before simulations. The small cells or elements are developed to fill the volume of the model. They constitute mesh, with each cell representing a discrete space that defines the local flow. The flow physics is then represented mathematically by equations that are applied to each cell of the mesh. As a result, each cell may store time-dependent flow properties such as pressure, velocity, shear rate, or temperature [42]. Generating a high-quality mesh is critical for obtaining reliable solutions and ensuring numerical stability.
In geometry, cells are arranged into a single grid that is usually one of three grid types based on the shape and complexity of the geometry: structured grids, unstructured grids, and hybrid grids ( Figure 3). Structured grids are identified by regular connectivity. When modeling a structure grid in a 3D model, a hexahedron cell is a possible choice element type. This grid is very space-efficient and has superior convergence and resolution than an unstructured grid [43]. Cell elements can be irregularly connected to form an unstructured grid. This allows for any possible elements that a solver might be able to use. The unstructured grids typically employ tetrahedral cells in 3D models [43]. A hybrid grid is made up of both structured and unstructured grid types. This can be helpful in the optimization of grids with complicated geometries. A structured grid can be part of the geometry with regular structure, while an unstructured grid can be part of the geometry with an irregular structure. Unstructured or hybrid grids are widely utilized in 3DP due to the relatively complicated geometry of the 3DP models (Table 2).  The evaluation of the generated grid's accuracy, sensitivity, and robustness is a critical stage in CFD simulation for the 3DP. The grid quality should be evaluated using a specific grid sensitivity that measures the model prediction's sensitivity to mesh size. Guo et al. [28] described the effect of mesh density on the simulation outcome of a 3DP model for a variety of cereal products. The number of meshes in their study affected the total mesh quality, and as a result, the simulation result for the local shear rate showed different values according to the element sizes ( Figure 4). To guarantee that the simulated outcomes are independent of the mesh density, the number of element sizes is varied until the simulated variables (shear rate, pressure, or velocity) remain unaffected. Next, the minimal mesh criteria among the applicable element sizes are chosen. Furthermore, several mesh quality metrics supplied by the meshing program may be used to enable grid quality measurement: the orthogonal (particularly at boundaries), relative grid spacing (15% to 20% stretching is considered a maximum value), and grid skewness are a few examples.

Boundary Conditions
All CFD problems are defined by initial and boundary conditions, which must be specified correctly. The initial values of heat and mass flow variables must be specified at all solution points in the computational domain. When modeling heat and mass transfer in the 3DP model, the distribution of all flow variables at the inlet boundary must be specified. In the case of 3DP with molten chocolate, the flow variables include inlet volume flow rate, pressure, turbulent intensity, and temperature [47]. The 3D printer's piston speed along with the dimensional parameters of the 3D printer chamber is used to define the inlet volume flow rate. Temperature and flow property distributions at the inlet are usually specified based on the experimental data [9]. Outlet boundary conditions can be stated and used in conjunction with the inlet boundary conditions. This type of boundary condition is specified mostly where outlet velocity is known. At the outlet of the flow domain, the force conditions are set (Equation (7)); where subscript n and s represent the normal and tangential components of the force and velocity, respectively. When the 3D printer's outlet position is far away from the chamber, the flow reaches a fully developed state with no change in the flow direction. As long as a fully developed state is maintained, the gradient of all variables could be equated to zero in the flow direction, except pressure [41]. In the chamber of the 3D printer, the continuity and the momentum conservation equations (Equations (8) and (9)) are usually used: where − ν is the velocity vector; ρ is the material density; p is the hydrostatic pressure; and = τ is the extra stress tensor. Mostly, the effect of gravity is neglected [10,21,28]. This seems to be a reasonable assumption in materials with high viscosity. At the wall boundary condition, models are usually established as no-slip between the material and the wall of the channel during the extrusion process wherein the normal component of velocity is fixed at zero, and the tangential component is set equal to the velocity of the wall (Equation (10) where V s and V n represent tangential and normal components, respectively. Although it may seem counterintuitive, the no-slip condition has been firmly established in both experiment and theory [48]. The wall is the most common boundary condition considered when solving heat and mass flow problems, specifically during the structural analysis of 3D printed food components. The bottom face of the printing bedplate on which molten materials are printed is supplied with a steady pre-heat temperature. Conduction heat transfer between the molten material and the bedplate occurs, and the temperature difference is important to control the cooling rate of the bottom layer. For the boundary conditions involving heat transfer, it is necessary to specify the rate of heat transfer through the walls. In general, the heat flux . q is defined by (Equation (11)): .
where T bed and T sample are the bedplate temperature and molten sample temperature, respectively. The overall heat transfer coefficient K depends on the thermal conductivity and thickness of the conductive solid material as well as the fluid material properties (density, heat capacity, etc.) and flow characteristics (geometry specification, turbulence, etc.). The radiant heat coefficient r s is determined by factors specific to the printing environment or the ambient temperature.
The symmetric boundary condition assumes that identical physical processes exist on both sides of the boundary. In the symmetric boundary conditions, all the variables have the same value and gradients at the same distance from the axis or plane of symmetry. It acts as a mirror that reflects all the flow distribution to the other side [49]. The conditions at the symmetric boundary allow no flow and no scalar flux across the boundary (Equations (12) and (13)); v n = o (mm/s) (12) f t = 0 (N) (13) where subscript t and n represent the tangential and normal components of the velocity and force, respectively. Here, n indicates the normal direction to the free surface. This is a great tool to reduce the computational effort and resource if the flow can be envisaged to be symmetrical about a plane or pair of planes. However, it must be noted that the geometrical symmetry does not guarantee the symmetry of the flow. Similarly, in cases where micro-structure of flow eddies are being captured such as "Large Eddy Simulation" or "DES-Detached Eddy Simulation", symmetry cannot be used due to the inherent nature of the eddies. As a result, this approach has not been used in many research publications on CFD analysis for 3DP processes.

CFD Applications in Food 3DP
CFD offers research tools for improving the process design and understanding of the fundamental physical basis of fluid dynamics. It provides benefits to food 3DP in many areas including evaluating the flow characteristics in the nozzle, die-swell analysis, residual stress analysis in the printed sample, material deformation after printing, and solidification of the printed sample, among others [9,21,28]. Vast improvements have been made in these fields in recent years. For food 3DP, the flow field analysis usually requires no mass transfer calculation, unlike CFD simulations for other unit operations used in the food industries such as spray drying of tomato pulp [17,25,44]. However, 3-dimensional flow field modeling is necessary to provide essential information regarding the 3D printing system design and the final product design. In addition, CFD can model several situations influenced under different circumstances faced during 3DP including the characteristics of the fluid at the inlets, inside the 3D printer chamber, nozzle domain, and outlets at the nozzle exit. Characterizing the flow field can also be utilized to validate a solution's degree of confidence before moving on to additional transport models or large-scale printing procedures. The following subsections are some of the properties that have been examined using CFD in 3D printing applications [9,21,28,44].

Thermal Mechanical Behavior in 3D Printed Food Material
The thermo-mechanical characteristics of food materials are studied in the food 3DP process to predict material behavior during or after printing in a wide range of factors describing their internal states (i.e., temperature, residual stress, and deformations) and structures (i.e., model shape or permeability). Changes in state parameters and structural characteristics of printed food material are induced by the energy exchange and mechanical interaction of a material with the environment. When a material is deformed, the strength describes the critical stresses. The residual stress is usually taken into account throughout the 3D printing process [9].
Residual stresses are the stresses maintained by a printed solid model after external or internal force loads have been removed and an equilibrium condition has been reached. These stresses are caused by misalignments in the form of components from different regions and phases. In addition, the residual stress is influenced by local variations of material properties such as elasticity and thermal properties. Residual stress can have several negative effects on the characteristics of a model part such as bending of the surface of the printed model or, in extreme cases, the entire deformation of the printed model [9,50]. The primary source of residual stresses in printed foods could be directly related to the material properties (viscosity, storage modulus, and loss modulus, etc.) and the parameters of the printing process such as nozzle size, printing speed, and layer height. Furthermore, the residual stress can be influenced by printing processes that incorporate temperature changes, which vary the material properties as well as cause phase transitions. Oyinloye and Yoon [9] investigated the mechanical behavior of printed alginate and pea protein paste at 45 • C on a 22 • C printing bed using numerical simulations. The finite element method (FEM) based numerical simulation was conducted to predict the residual stress and distortion of the printed material as the nozzle size was varied ( Figure 5). The stress of the samples was typically concentrated at the layers in the center of the sample toward the edges, and the stress value increased with increasing nozzle size. As the diameter of the nozzle increases, the extrusion flow rate increases and the printing time decreases. This limits the solidification and stability of each printed layer before printing another layer on top. Furthermore, it was observed that the increasing stress value became constant as soon as the build-up phase ended and the cooling phase began ( Figure 5). This was due to a decrease in the thermal gradient in the printed sample. In general, the major cause of stress is the presence of a high-temperature gradient during the build-up phase. The stress acts as an impediment to the expansion of molten material at the top layers due to the quick cooling process of the lower layers via thermal conduction with the printing bedplate. As a result, compressive stress is formed in these layers. In extreme temperatures, the yield stress of the material decreases, which causes stress formation that consequently produces plastic deformations. Additionally, shrinkage occurs during the material's cooling phase, affecting the material previously deposited under the layer. This results in a tensile residual stress in the top layer and compressive residual stress in the bottom layer [10].
In a temperature gradient environment, the generated stresses can be represented by Equation (14): In the simplified Equation (14), the residual stress σ (functions of temperature) at an arbitrary temperature T during cooling is given. It is a variant of Hooke's law σ = E ∈ , in which the strain replaces the linear shrinkage strain. The temperature interval over which residual stresses can build up is limited by the sample melting temperature [50]. Therefore, materials with a lower melting temperature reduce the residual stress [50]. The rate at which strain accumulates during cooling is then defined by thermal expansion (∝ CTE ). Finally, stiffness governs the conversion of these strains to stress (E). ∝ CTE , and especially E, are both dependent on temperature, with ∝ CTE typically increasing with increasing temperature, and E reducing to 0 at the melting temperature. E(T) and ∝(T) are temperature-dependent, and these parameters, along with T m , are material-specific properties. Hence, the only controllable parameter in the above equation is T, which may be either the bedplate preheat temperature or the temperature of the printing environment. In the CFD process, it may not be adequate to completely define the stress state since additional parameters affect the residual stress. These are entirely dependent on user input such as nozzle speed, mesh quality, and layer thickness [21,28]. However, it should be noted that the residual stress is also bounded by the material's melting temperature, which is a material property that cannot be user-defined.
A less severe but still detrimental consequence of the residual stress is the deformation of the printed parts. The bending moment caused by the stresses leads to components curling upward, regardless of the shape of the part itself. This is explained in Oyinloye and Yoon's [9] study for the mixture of alginate and pea protein paste, where the edges of the horizontally built cylindrical model show a substantial positive deformation in the Z direction ( Figure 6). The residual stress, which increases with increasing nozzle sizes, accounted for a greater total deformation in the sample printed with a 1.0 mm nozzle size than in the sample printed with a 0.20 mm nozzle size ( Figure 6). The high-stress value (7.1 × 10 −3 Pa) obtained from the 1.0 mm nozzle size caused the deformation of the printed material. Furthermore, the high flow rate associated with the bigger nozzle size limited the required time for each printed layer to solidify before another layer is printed on top of it. Therefore, the deformation of the product is usually associated with 3DP, regardless of whether the magnitude of deformation is big or small. It is possible to conclude that extending the downtime between the printing of successive layers allows the entire structure to reach thermal equilibrium. Keeping downtime to a minimum enables some heat to accumulate inside the printed product, which lowers the temperature gradients and consequently causes residual stress. This method can also be used in a printing process that does not involve temperature change.

Die Swell Effect in the 3D Food Printing Process
The paste in the 3D printer nozzle is under stress, with some of the deformation energy being stored elastically [51]. As the food paste leaves the printing nozzle, the free boundary of the die allows the paste to assume a plug-flow velocity profile (i.e., constant velocity in the radial direction [52]). Because the paste is no longer constrained by the walls of the 3D printer nozzle, the stresses acting on the paste in the nozzle are relaxed, and the elastically stored energy is released, resulting in radial expansion of the paste known as die swell [51]. The swelling of the extrudate emerging from a die is typical of viscoelastic liquids that contribute to the quality achieved during extrusion-based food 3DP [1]. Liang [53] analyzed the relationship between extrudate swell ratio and the entry stored elastic strain energy during die flow by regression analysis (Equation (15)); where β 1 and β 2 are constants characterizing the material elastic properties. The parameters β 1 and β 2 are calculated using Gaussian quadrature. β 1 represents the swell value at a shear strain rate of . γ = 1 s −1 . The dependence of the die swell on the shear rate, which is not negligible, is described by β 2 . Figure 7 shows the relationship between the extrudate swell and wall shear rate for capillary 3D printing experiments with various barrel diameter to die diameter ratios ( D B/D D ) for a viscoelastic liquid [54]. The result showed that the extrudate swell increased with the wall shear rate. This behavior is explained by an increase in recoverable elastic energy. Furthermore, as shown in Figure 7, the extrudate swell decreased as the actual barrel/die diameters increased from 20/4 to 30/6, and increased sharply for larger barrel/die diameters (35/7 and 40/8). This meant that the actual diameters of the barrel and die used affected the material's extrudate swell behavior. The dependence of the swell value on the 3D printer barrel diameter to die diameter ratio is intriguing. This indicates that the amount of elastic energy stored in the paste depends on the diameter of the 3D printer barrel up to a certain value, after which the conical zone of the converging flow streamlines into the die and is unaffected by the barrel wall. This fact demonstrates that the swelling of the extrudate is caused by the recoverable elastic energy stored in the paste, which is in turn caused by the large pressure drops at the capillary die's entrance.
According to Müllner et al. [55], the die swell ratio for a viscoelastic liquid is also affected by the temperature of the material and the length of the printing chamber up to the nozzle tip (residence time). However, this has not been evaluated in the food 3DP process, which could be an interesting area for further research. In general, understanding die swell during 3D printing aids in producing a high-quality printed sample. However, experimentally quantifying the die swell ratio during the food 3DP process is difficult due to several challenges including the small and closed 3DP environment and the small size of the swelling effect, which is sometimes not visible to the naked eye.
Yang et al. [1] demonstrated that CFD could provide information on the elastic nature of food paste. CFD was used to investigate the effect of inlet volume rate on the die swell ratio for lemon juice gel during 3DP. As shown in Figure 8, the extrusion at the exit of the die exhibited a certain level of swelling, with varying degrees of swell at different inlet volume rates. As the flow rate of the entrance volume increased, the extrusion swell ratio decreased first and then increased. This is most likely because the extrudate swell ratio is affected by the volume flow rate, shear rate, and pressure, all at the same time. According to the report, increasing the entrance volume flow rate increased the pressure in the flow channel, resulting in the fluid storing more elastic potential energy. In general, a precise and adequate quantitative description of flow is required for the 3DP process. Many CFD studies have demonstrated that die swell can contribute to unequal tensile stress in polymer samples, resulting in sample deformation during 3D printing [54,[56][57][58]. This is because when die swell occurs, the averaged velocities at different points around a profile are frequently not uniform, resulting in varying degrees of change to the local thickness. It should also be noted that as the velocities across the paste vary, the material size changes differentially across the profile wall thickness [54]. As a result, producing uniform averaged exit velocities at the nozzle tip is a desirable aspect in 3D extrusion, although it is difficult to achieve due to the limit of design specifications. Excessive change in extrudate shape may cause the material to stretch unacceptably, resulting in tensile stresses that can cause profile defects. Therefore, CFD simulations are used to precisely define the specific values of changes in extrudate size as a result of the die swell effect from the printing process parameter or the material properties. Thus, quantitative analyses can link the viscoelasticity of material being printed to the swelling of the extrudate, thereby improving the quality of 3DP. Unfortunately, this field has not yet been fully explored for food 3DP and could be a useful area for further research.

Assessment of Printing Quality with CFD
The quality of the final products can be influenced by various process parameters in food 3DP. The primary process parameters in 3DP are nozzle size/diameter, layer height, printing speed, extrusion rate, and material rheological properties (viscosity, storage, and loss moduli). To conduct CFD simulations for 3DP, the following variables must be set: preheat temperature, deposition thickness, dwell time, scanning speed, hatch spacing, layer thickness, density, number of heat sources, cool down condition, and working atmosphere. Previous research for extrusion-based food printing has focused on two critical factors in determining product quality: the printability of the materials during extrusion and the stability of the products after printing [9].
Food material printability is determined by their rheological (dynamic viscosity), thermal (i.e., melting point and glass transition temperature), and gel properties, which are important criteria for CFD processes [9,59,60]. Most commercial CFD software packages provide visible 3D/2D fluid characteristics values such as pressure, velocity, shear rate, and displacement with the input of fluid parameters (density and dynamic viscosity), assisting in a better understanding of the extruding processes. The Carreau model is popularly used to estimate the viscosity distribution of the material in the extruder from the experimental values. The model also describes the effect of shear rate on material viscosity (Equation (16)).
where η is the viscosity; η 0 is the zero shear viscosity; η ∞ is the infinite shear viscosity; λ is the relaxation time; γ is the shear rate; and n is the rheological index. Many publications on CFD analysis for extrusion processes have successfully used this model for the simulation and analysis of the flow properties in the 3DP extrusion process [28]. Liu et al. [44] used CFD to examine the effects of rheological properties on the printability of high and low gluten wheat flour, and found that pressure had an impact on the quality of 3D printed products. When the pressure applied to the material was too low during the extrusion process, the material could not be extruded out of the nozzle, and when the pressure was too high, a poorly textured model was printed. This was also justified by a reduction in viscosity as the printing pressure increased (Figure 9). As a result, the material extruded has limited time to expel the sample's inherent shear stress. Guo et al. [28] used the Bird-Carreau model to simulate and investigate the 3D printability performance of five cereal grains: buckwheat, brown rice, mung bean, job's tear seeds, and black rice. The result of the CFD simulation was validated with real printing experiments and rheological tests over a narrow range of shear rates. The viscosity of the material at eight different points in the 3D printer chamber is described in Figure 10. Because the model geometry causes the fluid flow to converge toward the centerline, the shear rate along the wall is expected to increase downward toward the nozzle tip. Figure 10 also shows that the viscosities of the materials were higher in the syringe, particularly at the corners and top center, where the shear rate value was expected to be lower. The viscosities gradually decreased toward the nozzle tip. Within the nozzle itself, viscosities were predicted to decrease gradually from the nozzle center to the wall. Furthermore, the CFD analysis revealed a similar trend to the rheological test in terms of descending order of viscosity for the five materials at the inlet position, which were buckwheat, mung bean, black rice, brown rice, and job's tear seeds. However, the corresponding values at the nozzle tip were mung bean, brown rice, black rice, job's tear seeds, and buckwheat. This was not the case with the rheological test results. This could be because the range of shear rates experienced by the grain gels during syringe flow (2.5 × 10 −4 to 726.6 s −1 ) was far beyond the range that the rheological test used in the study could cover (0.1 to 100 s −1 ). Therefore, the dynamic viscosity test results cannot fully and accurately describe the changes and differences in the viscosities of different materials during the extrusion. According to Yang et al. [1], a similar experiment was conducted with CFD simulation to investigate the distribution of velocity, shear rate, and pressure in the 3D printer chamber and nozzle at different process parameters such as material viscosity, relaxation time, inlet volume flow, and nozzle diameter for lemon juice gel. The above-mentioned study's simulation result was used to improve the 3DP of the gels and solve the problems from the die swelling phenomenon during extrusion (Figure 8). CFD analysis revealed that as the volume flow rate increased, so did the velocity and shear rate field. The pressure in the flow channel increased as the material viscosity and flow rate increased, while relaxation time and nozzle diameter decreased ( Figure 11). The greatest influence on the pressure field was discovered to be a change in nozzle diameter. Based on the aforementioned assertions, it is clear that numerical simulation could be used to better understand the distribution of flow properties (viscosity, shear rate, and pressure) of food materials during extrusion, particularly at the nozzle exit (outlet). Furthermore, numerical techniques can quickly and visually simulate and predict an object's stress distribution and fluid state without the need for repeated printing experiments.

Challenging Issues Confronting the Application of CFD in 3DP and Possible Solutions
In many ways, CFD has already been proven to be a valuable tool for the food processing industry. However, despite the overwhelming number of possibilities and benefits provided by the current CFD commercial codes, there are limitations. The calculated results represent a flow model that adheres to the physics and boundary conditions specified by the user. Proper physical modeling of the phenomena under investigation is thus a critical step in the preparation of a CFD simulation and will dominate the relevance and accuracy of the results obtained later. The effectiveness and practicability of a CFD simulation in food 3DP are dependent on several factors including specific food material properties, accurate equations of motion algorithms, powerful CFD packages, and hardware with high computing power. This necessitates a thorough understanding and justification of the models of all physical and chemical processes considered in the computations.
The materials used in food 3DP (fruit juice, milk, gels, meat paste, etc.) differ in many ways (rheological flow properties, thermodynamic properties, physical properties, etc.) from those to which CFD is conventionally applied (air, water, etc.). Conventional CFD has been able to predict only those 'mixture average' properties; however, the quality of food material is closely associated with its particle size, elastic properties, texture, and composition. The 3DP process can only be understood and managed effectively if physical-chemical, organic-chemical, and micro-biological models are available; unfortunately, traditional CFD programs do not support such processing. Furthermore, for 3DP simulation, numerous complexities must be considered such as non-Newtonian and viscoelastic materials, yield stress, frictional force along the wall, and simultaneous analysis of fluid flow and solid stresses (for example, fluid-structure interaction). Unfortunately, to make the numerical problems more approachable, several of these characteristics are considered to play a less crucial role in 3D CFD applications. Researchers, for example, make assumptions such as (i) no-slip conditions at the wall; (ii) that the thermal and physical properties during analysis are constant; (iii) that the material is at rest in the initial condition; and (iv) that the material observes a steady and laminar flow in the flow channel, among others [1,9,21,28]. However, these assumptions are often improper. A typical example is a 3DP research in which laminar flow is completely assumed for the material flow channel [1]. However, turbulent and laminar flow could still exist in the same problem. This has been reported in a simulation study of annular Couette flow in which turbulent and laminar flow regimes coexist [61]. Recent modeling advances have addressed this issue by developing a predictive laminar to turbulent flow transition model, which is now included in the ANSYS CFX ® 10.0 program upward [62]. However, no study using this model for the 3DP has been reported as of yet. Both the transient process and the steady-state process are used to analyze several flow regimes in 3DP [1,9,21,28]. The transient process arises as a result of the moving inlet boundary (i.e., the piston in the 3D printer chamber). In these cases, a steady-state flow regime does not exist and numerical difficulties are often encountered when attempting to solve the problem using steady governing equations [63]. Nonetheless, a steady solution can be 'forced', which means that some limiting constraints can be imposed to inhibit the flow regime's unsteady characteristics. These typically include one or more of the following: first-order convection schemes, dissipative turbulence models (standard ε-k model), and two-dimensional or symmetrical boundary conditions. However, because this driving action does not exist in nature, the credibility of solutions derived from such calculations has been doubted [63]. Before starting CFD modeling, it should be noted that the limitations of the available turbulence models and their associated wall functions must be considered. Models suitable for the study should be chosen based on previous experiences with similar applications in the literature.
Other challenges faced might be connected to the complexity of the geometry and the mesh quality, both of which have a significant role in deciding the correctness of the final result. The computational mesh offers spatial discretization of the governing equations and is a major vehicle for ensuring correct predictions of the fluid continuum, particularly in steady-state simulation. To obtain a high degree of precision, the mesh must be properly refined in areas of interest and where gradients exist in the flow field. This has been extensively discussed in Section 3.3. However, there are also situations when low mesh quality is reported after mesh refinement [64]. These issues can be unavoidable in some situations, but in many others, identifying and repairing problematic mesh areas can eliminate such difficulties. This means that cells with high aspect ratios or highly skewed cells must be eliminated from the simulation to obtain an accurate and efficient representation of convective fluxes. Some CFD programs provide methods for locating these cells. However, in many situations, the CFD modeler will rely on their prior knowledge with a CFD software program to validate the mesh quality. The use of CFD in food 3DP has not become widespread because, to date, several obstacles have stood in the way of its development in this area: a lack of basic knowledge of local kinetic processes at small scales (paste droplet or flow), the complexity of biological processes with the scarcity of models for those that are compatible with CFD and the nature of the food-processing industry and its relation to CFD, the difficulties linked to the cost of purchasing or leasing the software, the cost of high-powered computer hardware, and the scarcity, cost, and managerial difficulties of ensuring adequately trained and experienced personnel.
More efforts are needed to promote the use of mathematical models for 3DP in the food sector. Model prediction accuracy might be enhanced by using correct heat and mass transport coefficients. To obtain a more precise outcome, the CFD simulation in 3DP might be coupled with other techniques such as optimal shape design (OSD) and discrete element method (DEM). Sarghini [65] described the use of OSD in conjunction with CFD to analyze the extrusion process of paste dough. In OSD, the essential element concerning numerical simulations in fixed geometrical configurations is the introduction of a certain number of geometrical degrees of freedom as one of the unknowns, meaning that the geometry is not completely defined, but part of it is allowed to move dynamically to minimize or maximize the objective function. The applications of optimal shape design (OSD) are found in many places: structure mechanics, electromagnetism, and fluid mechanics, or a combination of these [65,66]. Furthermore, the DEM model technique enhances the modeling of flows involving an equal amount of fluid and powders that cannot be analyzed solely by the continuum flow based models used in CFD simulation [67]. DEM has recently been used to model a wide range of granular mixing applications with high qualitative accuracy [67]. This method can account for powder cohesion and has been combined with CFD to simulate the transport of powder materials through pneumatic pipes [68]. This will be useful for simulating mixture materials, particularly samples containing both fluid and solid phases such as meat products and starch paste. However, the computational cost of DEM is very expensive, and simulations often take many days to complete. This has made it unrealistic to extend this model to other modes of dense gas-solid flow exhibited by fine powders (particle size less than 100 mm). As a result, this technique may not be utilized in food 3DP processing immediately [69].
Additionally, validation of simulation with the experimental result is paramount for the future success of CFD modeling in the 3DP. As of now, most publications have validated the 3DP simulation model by comparing the images of the printed sample with the simulation results or the material properties measured [1,28]. However, a more accurate method such as the image process such as that used in the validation of a printed model of alginate and pea protein [9] may be introduced. This method could provide time-dependent changes occurring during the printing process such as the die swell effect discussed in Section 4.2 and the post-printing process such as time-dependent deformation in the printed sample.

Conclusions
This paper reviewed the application of CFD in food 3DP. The examined studies demonstrate the need to maintain a high degree of accuracy through careful decisions made during model development such as selecting an acceptable numerical model, refining the mesh, and selecting appropriate boundary conditions. The numerical simulation could predict the distributions of average velocity, local shear rate, and pressure of the materials during the extrusion process. The link between the rheological properties of the materials and their printing qualities was established, and it is clear that material properties play a significant role in achieving printability and application. To simulate the flow characteristics of food material using CFD, it is necessary to understand how the essential constituent of food material behaves during the printing process. The present state of CFD has been demonstrated to be of significant value in its use in 3DP, which extends from the study of residual stress to total deformation in the printed product and die swell effect at the nozzle tip, all of which contribute to a successful printing process. Additionally, future study areas in 3D CFD analysis such as the effect of die swell on material deformation is also emphasized. This field has the potential to improve the printing process with food materials.

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