Multiscale Eulerian CFD of Chemical Processes: A Review

: This review covers the scope of multiscale computational fluid dynamics (CFD), laying the framework for studying hydrodynamics with and without chemical reactions in single and multiple phases regarded as continuum fluids. The molecular, coarse-grained particle, and mesoscale dynamics at the individual scale are excluded in this review. Scoping single-scale Eulerian CFD approaches, the necessity of multiscale CFD is highlighted. First, the Eulerian CFD theory, including the governing and turbulence equations, is described for single and multiple phases. The Reynolds-averaged Navier–Stokes (RANS)-based turbulence model such as the standard k -  equation is briefly presented, which is commonly used for industrial flow conditions. Following the general CFD theories based on the first-principle laws, a multiscale CFD strategy interacting between micro- and macroscale domains is introduced. Next, the applications of single-scale CFD are presented for chemical and biological processes such as gas distributors, combustors, gas storage tanks, bioreactors, fuel cells, random- and structured-packing columns, gas-liquid bubble columns, and gas-solid and gas-liquid-solid fluidized beds. Several multiscale simulations coupled with Eulerian CFD are reported, focusing on the coupling strategy between two scales. Finally, challenges to multiscale CFD simulations are discussed. The need for experimental validation of CFD results is also presented to lay the groundwork for digital twins supported by CFD. This review culminates in conclusions and perspectives of multiscale CFD.

Multiscale modeling has been used in the form of sequential and concurrent couplings [4,15]. In a sequential coupling framework, often referred to as parameter passing between scales, the macroscale model uses a small number of parameters calculated in the microscale model. The concurrent coupling strategy computes several scale models simultaneously [1,4,16]. Concurrent coupling is preferred over sequential coupling when the missing information is a function of many variables. Concurrent coupling is a powerful approach but it is also computationally prohibitive [1].
Delgado-Buscalioni et al. (2005) presented a multiscale hybrid scheme coupling MD in a subdomain with computational fluid dynamics (CFD) in a continuum fluid domain [1]. Nanofluids typically employed as heat transfer fluids include multiphysics such as drag and lift forces between liquid and nanoparticles and Brownian, thermophoretic, van der Waals, and electrostatic doublelayer forces [18]. Challenge remains in nanofluid flow simulations using a CFD model that includes multiphysics and a multiscale CFD model. Tong et al. (2019) reviewed multiscale methods divided into a domain decomposition scheme and a hierarchical scheme in fluid flows with heat transfer for coupling MD with particle-based mesoscale methods such as the lattice Boltzmann method (LBM) and CFD [15]. LBM using the particle distribution function is suitable for simulating fluid flows involving interfacial dynamics and complex boundaries [15,19]. The fundamental idea of LBM is that the macroscopic dynamics of a fluid are the result of the collective behavior of many microscopic particles in the system [19]. Thus, by integrating the particle momentum space, the particle distribution function (Boltzmann transport equation) at the microscopic scale is linked to the continuous fluid dynamic properties such as density, velocity, and energy at the macroscopic scale [3,20]. LBM is an alternative approach in CFD, which has been applied to fluid flows in porous media, microfluidics, and particulate flows [19,20]. However, because CFD based on LBM focuses on the microscopic domain with complex boundary conditions, macroscopic Eulerian CFD is preferred for investigating hydrodynamics on the process level [2,15].
Lagrangian CFD approaches such as the discrete element method (DEM) and the discrete phase model (DPM) solves the Newton's second law of motion for each particle to identify the trajectories of the particles. The Lagrangian model is normally limited to a relatively small number of particles (less than 10 5 particles) because of the computational expense [17,21]. This review focuses on the Eulerian CFD for the continuum phase in chemical processes, excluding the more sophisticated methods such as the LBM and Lagrangian approaches.
Eulerian CFD has become a powerful and important tool for simulating and predicting fluid behaviors in chemical and biological processes [22]. Eulerian CFD provides the average quantities at the equipment-scale that are of practical values to engineers [23]. CFD is used for investigating the hydrodynamics of a process following geometrical and/or operational modifications [22,24]. CFD has become crucial for understanding physical phenomena in two-(2D) and three-dimensional (3D) geometries and for scaling up, optimizing, and designing chemical and biological processes [25]. However, CFD models must be validated by comparing the simulation results with the experimental data to provide meaningful information [26]. Modeling, meshing, the physical properties, and selection of a suitable turbulence model play an important role in CFD analysis. The Eulerian CFD model is represented by the mass, momentum, and energy conservation laws described by partial differential equations (PDEs) in 2D or 3D space. As the numerical method to convert PDEs into a set of algebraic equations (AEs), the finite element method (FEM) [27,28] has been traditionally used to investigate the mechanical and structural properties of materials, whereas the finite volume method (FVM) [29], which ensures conservative fluxes within a finite volume, has been often used for fluid dynamics such as Eulerian CFD.
Fluid dynamics are governed by the Navier-Stokes (NS) equation representing the conservation of momentum. Turbulence is involved in the NS equation to consider random fluctuations of fluid motion [3]. There are several turbulence equations, such as the direct numerical simulation (DNS), large eddy simulation (LES) [30], and Reynolds-averaged Navier-Stokes (RANS)-based k- models. The uncertainty of RANS-based k- models, which originates from information loss in the Reynoldsaveraging process, was discussed by Xiao and Cinnella (2019) [31].
Multiphase physics is omnipresent in both environmental and engineering flows [32]. As the precision and accuracy of manufacturing techniques progress, micro-and nanotechnologies become crucial for providing engineering solutions to problems across different industrial sectors [3]. The multiscale CFD approach has been applied to various disciplines for more accurate solutions than can be obtained using single-scale modeling. Ngo (2018) presented multiphase, multiphysics, and multiscale CFD simulations for gas-solid fluidized beds, steam methane reforming reactors, and impregnation dies for carbon fiber production [29,33]. Da Rosa and Braatz (2018) [34] proposed a multiscale CFD model for a continuous flow tubular crystallizer. The micromixing model, energy balance, and population balance equation (PBE) were coupled with the CFD model. Haghighat et al. (2018) [35] combined a 3D CFD model for smoke flow and a 1D fire event model as a far-field boundary condition to predict hydrodynamics in a road tunnel. Hohne et al. (2019) [36] combined a generalized two-phase flow boiling model with CFD to predict the breakup, coalescence, condensation, and evaporation mechanisms in a heated pipe. Uribe et al. (2019) [37] compared three different CFD models (heterogeneous micropores model, pseudo-homogenous catalyst particle model, and single-scale reactor model) in a trickle bed reactor (TBR) to show its multiphysics and multiscale nature.
The behavior of the macroscale flow is affected by microscale physical processes, which also leads to multiscale CFD modeling [3]. Multiscale behavior with complex physical phenomena, which are highly interrelated [17], appears in process engineering. For example, Figure 1 illustrates a carbon fiber (CF) production process [29,33] from a multiscale point of view. Figure 1a shows CF tape with a width of 100 mm as the final product. Fifteen tows enter the impregnation die and are impregnated with a polymer resin (see Figure 1b). The CF tape is produced from the die after it is grooved, spread, pressed, and cooled [29]. One tow having a width of 6 mm and a thickness of 130 m is shown in Figure 1c, magnified from Figure 1b. A microscale representative element volume (REV) with 128 randomly distributed CF filaments (7 m in diameter) is depicted in Figure 1d. A macroscale Eulerian CFD can be used for the impregnation die ( Figure 1b) to investigate the hydrodynamics of the process, whereas a microscale Eulerian CFD can be applied to the creeping flow of the resin inside the tow (Figure 1d). In the concurrent coupling framework, the macro-and microscale CFD models are solved simultaneously or the two models exchange information at each time and space. If the tow domain is assumed to be a uniform porous medium, the permeability of the resin through the tow is needed in the macroscale porous CFD model. In the sequential coupling framework, the permeability computed from the microscale porous CFD model is used in the macroscale CFD model [29].
When the different models at different scales are coupled together, sequentially or simultaneously, some errors often appear at the interface where the two models meet [4]. One of the challenges is to determine how the multiscale models can be coupled smoothly. Section 3.5 shows an example of the coupling strategy. Many reviews in the Eulerian CFD simulation have been recently reported. Ferreira et al. (2015) reviewed gas-liquid CFD simulations without electrochemical reactions in proton exchange membrane (PEM) fuel cells [38]. Karpinska and Bridgeman (2016) [39] reported CFD studies on activated sludge systems in a wastewater treatment plant (WWTP). Pan et al. (2016) presented gasliquid-solid CFD models for fluidized-bed reactors [23]. Sharma and Kalamkar (2016) [40] used gas-phase CFD models to optimize geometrical and flow parameters that led to designing a roughened duct for the best thermohydraulic performance. Yin and Yan (2016) [41] reviewed gas-phase CFD studies on oxy-fuel combustors in a power plant. Uebel et al. (2016) reported the CFD-based multiobjective optimization of a syngas conversion reactor [42]. The interplay between electrostatics and hydrodynamics in gas-solid fluidized bed was reviewed by Fotovat et al. (2017) [43]. Pires et al. (2017) [44] addressed the recent progress in CFD modeling of photo-bioreactors for microalgae production. Bourgeois et al. (2018) reviewed CFD models in gas-filling processes, investigating the temperature evolution of the gas and the tank wall during H2 filling [26]. Malekjani and Jafari (2018) presented recent advances in CFD simulations with heat and mass transfers of food-drying processes [45]. Pinto et al. (2018) [28] applied a CFD approach to physical vapor deposition (PVD) coating processes. Ge et al. (2019) reviewed the general features of multiscale structures in particle-fluid systems [17].   [5,18] reported recent advances in modeling and simulation of nanofluids which are a mixture of a common liquid and solid particles less than 100 nm in size, focusing on 3D Eulerian CFD models for thermal systems [5]. Drikakis et al. (2019) presented the application of CFD to the energy field [3], integrated together with molecular dynamics and LBMs. Lu et al. (2019) [46] and Wang (2020) [47] reviewed Eulerian CFD approaches for dense gas-solid flows, providing a concise introduction to multiscale methods and highlighting the effects of mesoscale structures (gas bubbles and particle clusters) on the gas-solid Eulerian CFD model, focusing on the energy minimization multi-scale (EMMS) method.
Nevertheless, few researchers have reviewed Eulerian CFD studies applied to micro-and macroscale hydrodynamic problems coupled with each other. This review covers the broad scope of CFD for chemical and biological processes. In particular, it focuses on i) multiscale CFD studies in the past decade, ii) continuous phase rather than discrete phase, which can be formulated by Eulerian and volume of fluid (VOF) models, iii) incompressible fluids rather than compressible fluids, and iv) applications to chemical and biological processes. In this review, multiscale CFD encompasses a physical domain described by a continuum, excluding molecular, coarse-grained particles, and mesoscale particle dynamics at the individual scale. Eulerian CFD models in single and multi-phases are first presented, followed by applications of the Eulerian CFD to chemical and biological processes, and finally the challenges and perspectives of multiscale CFD are discussed.

Eulerian CFD Models
The main CFD equations are derived from the first-principles laws for mass, momentum, and energy conservations in an infinitesimal volume of space [3], which are called the governing equation. The Reynolds averaged Navier-Stokes (RANS) fluid governing equations are often used in Eulerian CFD studies. The CFD models are described for single and multiple phases with the gas, liquid, and solid phases. The standard k- turbulence equation is presented as an example of the RANS-based turbulence equations.

Single-Phase Eulerian CFD Model
A gas or liquid phase as a continuum is modeled by a continuity equation, a Navier-Stokes (NS) momentum equation, and an energy equation [24,48]. Table 1 shows the Eulerian CFD model for an incompressible fluid. The density () is constant, and the velocity vector ( ⃗) in 2D or 3D is solved from the NS equation. The stress tensor ( ̿ ) is defined with the molecular () and turbulent ( ) viscosities for the turbulent flow. The turbulent kinetic energy (k) is described in Section 2.5, and ⃗ is the acceleration of gravity. The temperature (T) is obtained from the energy equation, Equation (T1-3), including the convective energy (the first term in the right-hand side), the diffusive energy with conduction and viscous dissipation (the second term in the right-had side), and the heat sources (Sh) such as reaction heat and radiation. The cp and  are the heat capacity and effective conductivity, respectively. Table 1. Single-phase Eulerian computational fluid dynamics (CFD) model [24,48].

Continuity equation:
⃗ • ⃗ = 0 (T1-1) Momentum equation: Energy equation: The equation of state (EOS) is needed to connect the state variables of pressure (P), volume (V), and temperature (T). As a result that Equation (T1-2) represents three equations for four unknowns (ux, uy, uz, and P) with Equation (T1-1) as a constraint on the velocity, the semi-implicit method for a pressure-linked equation (SIMPLE) or the pressure implicit method with splitting the operators (PISO) is often used to link the velocity and pressure [3,49].

Gas-Liquid Eulerian CFD Model
Gas-liquid phases are modeled by the Eulerian approach assuming that the two phases flow as non-interpenetrating or interpenetrating continua. The Eulerian model assuming noninterpenetrating continua is often called the volume of fluid (VOF) method, which is a surfacetracking technique for immiscible fluids (hereafter VOF-CFD). The Eulerian model assuming interpenetrating continua is the Eulerian multiphase approach solving the continuity and momentum equations at each phase (hereafter EM-CFD).
The VOF-CFD model in an incompressible gas-liquid phase [33,50] is presented in Table 2. The VOF-CFD model relies on the assumption that the two phases do not interpenetrate [33], as mentioned earlier. The gas and liquid phases are considered as the primary and secondary phases, respectively. The VOF-CFD typically solves a single set of continuity, momentum, and energy (E) equations weighted by the phase fraction () [51]. The volume fraction of the secondary phase (L) is solved by Equation (T2-2), which is a continuity equation for the liquid phase where the interface between the phases is captured.
The surface tension force ( ) can be defined as a continuum surface stress (CSS) [50] tensor ( ̿ ) or a continuum surface force (CSF) [33]. The CSS method represents the surface tension force in a conservative manner and does not require an explicit calculation of the surface curvature. For turbulent flows, the turbulence equation is added to the VOF-CFD model (see Section 2.5).   Table 3 lists the EM-CFD model for the gas-liquid phase. The governing equation includes the continuity equation for the total mass conservation, the momentum equation with drag and non-drag forces, and the energy equation for each phase. The gas phase is considered as the secondary phase (or dispersed phase). The interfacial forces ( ⃗ ) between the gas and liquid phases in Equation (T3-3) can be formulated by the linear sum of the drag force ( ⃗ ), the lift force ( ⃗ ) emerging from the interaction between the gas (or bubble) and vorticity in the liquid [52], the wall lubrication force ( ⃗ ) pushing bubbles away from the wall, and the turbulent dispersion force ( ⃗ ) accounting for the turbulent momentum transfer between turbulent eddies and bubbles [53]. ⃗ is proportional to the difference between the gas and liquid velocities in Equation (T3-4). The gas-liquid momentum exchange coefficient (KGL) contains the drag coefficient (CD) [53]. QGL in Eqs. (T3-8) and (T3-9) is the heat transfer between the gas and liquid phases.
Energy equation: The mass, momentum, and heat transfers between the two phases play a key role in multiphase CFD simulations. The interfacial transfer terms require correct assessment using analytical models and empirical correlations.

Gas-Solid Eulerian CFD Model
The gas-solid CFD model is often expressed by Eulerian CFD with the kinetic theory of granular flow (KTGF) [55,56]. The KTGF-CFD model in Table 4 is the most practical choice for industrial-scale simulation of particle-fluid systems [17]. Equations (T4-1)-(T4-12) are the continuity and momentum conservation equations for the gas and solid phases. Equations (T4-17)-(T4-18) are the energy equations for each phase [57]. Here, is the solid volume fraction, and are the velocities of the gas and solid phases, respectively, and are the densities of the gas and solid, respectively, is the pressure shared by all phases,  is the stress tensor associated with the two phases, and ⃗ , is the momentum exchange between the phases, which is defined as the interaction force per unit volume of the bed [55,58].
Assuming incompressible flow (constant ), no chemical reaction, and no phase change, the first term of Equation (T4-1) accounts for the rate of total mass accumulation per unit volume, and the second term is the net rate of convection mass flux. and are expressed in Eqs. (T4-4) and (T4-5), respectively, for a Newtonian fluid with constant viscosity (). Table 4. Kinetic theory of granular flow (KTGF)-CFD model for gas-solid phases [55,58].

Continuity equation (gas phase):
( ) Momentum equation (gas phase): where the solid-gas drag [59] the gas-phase stress tensor is the solid-phase stress tensor is the solid shear viscosity is Energy equation: The Gidaspow drag model (KGS) [59] in Equation (T4-3), which is a combination of the Ergun and Wen-Yu drag models, is applicable to both dilute and dense solids flows [21]. The drag model was validated by a number of studies for bubbling fluidized bed with acceptable accuracy [55]. The granular temperature ( ) in Equation (T4-6) for the solid pressure (pS) is defined with the energy generation by the solid stress in Equation (T4-14), the dissipation rate ( ) of collisional energy in Equation (T4-15), and the random fluctuation kinetic energy transfer ( ) from the solid to gas phases in Equation (T4-16) [55,56].
The mesoscale structure of clustered particles has a significant effect on the flow and transport behavior of particle-fluid systems [17]. The drag coefficient (CD in Equation (T4-3)) measured in gassolid fluidization may be orders of magnitude lower than that of uniform suspension because of the presence of bubbles and particle clusters [17]. To elucidate the complex behavior of clustered particles, energy minimization multiscale (EMMS) methods have been proposed, where the energy consumption for the suspension and transport of particles is minimized under stability constraints [64].
Using the operating conditions (superficial gas velocity and solid circulation flux) and the physical properties of the gas (density and viscosity) and the solid (diameter and density) as input parameters, the drag coefficient of the EMMS method is calculated in sequential or concurrent strategy coupling with the KTGF-CFD model in the macroscopic scale [17]. The table-looking method as a sequential coupling strategy can be used for the EMMS method in a computationally efficient way (refer to Section 3.5 for details) [58].

Three-Phase Eulerian CFD Model
In the three-phase Eulerian CFD model, three phases are treated as fully interpenetrating continuous phases. The solid pressure (pS) and viscosity (S) are calculated by the KTGF model. All the phases have the same conservation equations, where additional closures (or constitutive equations) are required for interface interactions such as gas-liquid ( ⃗ , ), solid-liquid ( ⃗ , ), and gas-solid ( ⃗ , ) drags [23], and gas-liquid (QGL), solid-liquid (QLS), and gas-solid (QGS) heat exchanges [57]. The drag forces ( ⃗ ) caused by the slip velocity between phases are predominant in the interphase momentum exchange forces, therefore Table 5 shows only the drag forces in the momentum equation. Pan et al. (2016) [23] summarized the drag coefficient ( , ) between the gas and liquid phases such as in the Grace, Tomyama, Schiller-Naumann, and Ishii-Zuber models. The gas-liquid drag force influences the gas holdup. Hamidipour et al. (2012) [65] used the Gidaspow drag model [59] to describe the solid-liquid interaction force ( ⃗ , ), which has been applied to gas-solid systems. The gas-solid drag force ( ⃗ , ) similar to that between the gas and liquid phases is used in Equation (T5-3). In addition to the drag models, the granular temperature ( ) from the KTGF is calculated (see Equation (T4-13)Equation (T4-16)) for the solid pressure (pS) and viscosity (S) [65]. Table 5. EM-KTGF-CFD model for gas-liquid-solid phases [65].

Turbulence Equations
The dynamics of fluid flow are described by the NS momentum equations [31]. As the Reynolds number increases, the flow reaches a state of motion characterized by strong and unsteady random fluctuations of the velocity and pressure fields, which is referred to as the turbulent regime. Turbulent flow can be modeled by the direct numerical simulation (DNS) of the NS equations, which is computationally prohibitive but accurate, RANS based on empirical models, and large eddy simulation (LES) as a compromise between DNS and RANS [31].
The instantaneous velocity and pressure are decomposed into the sum of the mean and fluctuation components in the RANS-based equation. Substituting the decomposition into the NS equations and taking the ensemble-average leads to the RANS equations. RANS modeling includes uncertainties due to information loss in the Reynolds-averaging process [31]. Here, the Reynolds stress tensor ( ̿ ) including the turbulent viscosity (t) is solved with constitutive equations describing the turbulent kinetic energy (k) and its dissipation rate (). The RANS equation for the turbulent momentum equation is expressed as where the Reynolds stress tensor ( ̿ ) is defined as The standard k- turbulence model [50,66,67] is expressed for incompressible fluid as follows: where Gk and Gb represent the generation of the turbulent kinetic energy (k) caused by the mean velocity gradients and the buoyancy, respectively.
The k- model involves coefficients C, C1, C2, k, and . The nature of these coefficients leads to ambiguity regarding their values. However, these coefficients have been calibrated empirically to reproduce the results of a few flows [31]. Xiao and Cinnella (2019) presented the uncertainty of these parameters [31]. The RANS-based turbulence model [26] includes the standard k- [67], modified k-, realizable k- for complex flows [22,24,55,68], renormalization group (RNG) k-, and shear stress transport (SST) k- equations [53]. The RNG k-model is more accurate and reliable for a wider range of flows than the standard k-model, which is mostly used in turbulence modeling of duct channel flows with heat transfer [40]. The realizable k- turbulence model satisfies the mathematical constraints on the Reynolds stresses, which is consistent with the physics of turbulence flow. The realizable k- model has exhibited substantial improvements over the other k- type turbulence models when the flow features include strong streamline curvature, vortex, and rotation [68]. The SST k- model is widely used for laminar and turbulent mixed flows [53,69]. For turbulent two-phase flows in microscale CFD, LES may offer a good compromise between RANS-based closures and DNS due to suitable subgrid closures [32], as mentioned earlier. However, wide application of LES is still limited because of the large mesh requirement and high computational cost [39]. The DNS and LES may be impractical as a general-purpose design tool for most industrial flow conditions [39].

Single-Phase CFD Simulations
Gas-phase CFD with various turbulence equations has been used to find the optimal geometries of equipment such as a vortex tube producing a temperature separation effect [81], a syngas conversion reactor [42], a duct with rib-roughened walls [40], an industrial-scale gas distributor with injectors [24], and a sleeve-type steam methane reformer for H2 production [22]. For oxy-fuel combustion in power plants, a gas-phase Eulerian CFD model was incorporated with combustion reactions (pyrolysis, gasification, char reactions, and gas-phase reactions), turbulence, convective and radiative heat transfers, and convective mass transfer [41].
Using a liquid-phase 3D CFD model with heat transfer, the effect of various hollow fiber geometries on hydrodynamics through the flow channels was investigated to enhance the heat transfer coefficient (HTC) between the feed and the permeate [82]. Liquid-phase CFD under a laminar flow condition was used to evaluate the HTC of nanoparticle-added ionic liquids that are nonflammable and non-volatile at ambient conditions and recyclable [83]. The effect of the Reynolds number and nanoparticle concentration on the HTC was investigated using the Eulerian CFD model with the thermos-physical properties of ionic nanofluids such as thermal conductivity, heat capacity, density, and viscosity [83].

Gas Distributors
CFD is considered a pivotal technology in the design and performance evaluation of chemical processes [84]. Appropriate performance indexes should be obtained from CFD results for both the process design and optimization. Pham et al. (2018) [24] compared several gas distributors used in an industrial-scale CO2 amine absorber using a gas-phase steady-state 3D CFD with the realizable k- turbulence equation. Three performance indexes, including the pressure drop, dead-area ratio, and coefficient of gas distribution, were used to describe the quality of gas distribution [24]. Figure 2 shows a Schoepentoeter gas distributor as a gas injection device, the velocity vector around the Schoepentoeter, and the radial velocity contour above the chimney tray. The velocity contour obtained from the CFD simulation (Figure 2c) is able to quantify the uniformity of the gas distribution.

Gas-Filling Storage Processes
The maximum temperature is regulated throughout H2 refueling to preserve the integrity of the tanks [26]. CFD simulations of H2 tanks with experimental validation were performed in the range of 50  P  800 bar [26]. CFD made it possible to investigate the effect of feed temperature, density, pressure, and flow speed on the gas and tank wall temperatures in space and time [26]. A gas-phase CFD model with internal and external heat transfers, which was validated with experimental data, was used to identify the temperature profile of the H2 tank.

Fixed Bed with Random Packing
Ahmadi and Sefidvash (2018) [72] presented a gas-phase steady-state 3D CFD with the standard k- turbulence equation to predict the pressure drop in a fixed bed with random packing. A supercritical steam flow with high Reynolds numbers was simulated in a column randomly packed with 800 fixed particles of 3 cm in diameter. After validation of the CFD model against experimental data, the CFD model was used to develop a CFD-based empirical equation for the pressure drop at high Reynolds numbers in random packing fixed beds, which is useful in industry [72]. Unlike porous media CFD simulations for fixed beds [75,85], this CFD study considered the random structure of 800 spherical particles using 15 million computational cells. Wang et al. (2020) [73] reported a gas-phase steady-state 3D CFD with the realizable k- turbulence equation to predict the pressure drop in a rotating fixed bed with random packing. N2 gas flow was simulated in a rotating column randomly packed by 470 fixed particles of 5 mm in diameter. Approximately 10 million computational cells were used to build a fixed bed randomly packed with the spherical particles. The CFD model was first validated with experimental data. Subsequently, the effects of rotating speed and gas velocity on the pressure drop were examined using the validated CFD model. A semi-empirical correlation based on the CFD simulation was proposed to predict the dry pressure drop [73], which is of practical use in industry.
Using VOF-CFD models, the effects of operating conditions (temperature and gas flow rate), material properties, and gas channel geometries on the hydrodynamics of PEM fuel cells (cathode side) were reported by Ferreira et al. (2015) [38]. Multiple length scales from a few micrometers to meters were also present in the PEM fuel cell. Coupling of the VOF model with electrochemical reactions and heat transfer is needed to better understand the hydrodynamics of PEM fuel cells [38]. Furthermore, the coupling of gas-liquid flows in both anode/cathode channels and porous electrodes is vital in current CFD models, which leads to multiscale simulations [86]. Zhang and Jiao (2018) [86] presented the necessity of a 3D CFD simulation at cell and stack levels to accurately simulate water and thermal management in PEM fuel cells. Li et al. (2018) [87] used a 3D VOF-CFD model with CSF to investigate the pressure drop, liquid holdup, and gas-liquid interfacial area in a microporous tube-in-tube microchannel reactor. The pressure drop obtained from the CFD was validated against that of an experiment. The effects of the microchannel width, micropore diameter, and liquid flow rate on the pressure drop and interfacial area were examined using the CFD model [87]. The CFD study addressed the potential for optimization of operating conditions and design of microchannel reactors.
CFD tools have been used as a design tool in the development of new aeration devices and optimization of their geometry in WWTPs [39]. The gas-liquid EM-CFD model was applied to activated sludge tanks of WWTP to investigate the gas holdup, flow velocity, and O2 concentration [39]. The fluid flow in such bioreactors depended on the tank geometry, physical properties (density, viscosity, and solubility), and operating conditions (inlet flow rate, inlet concentration, and temperature) [39]. A CFD model combining hydrodynamics, mass transfer, and biochemical reaction kinetics will be one of the major challenges in WWTP [39]. Collision, coalescence, breakage, and deformation of bubbles promote diversity in the shapes and sizes of bubbles in aeration tanks with a bubbly flow regime. The population balance equation (PBE) for bubbles may be helpful in accounting for the bubble-size distribution (BSD) and its dynamics [39].
Soman and Shastri (2015) [76] presented a gas-liquid 3D Eulerian CFD with radiation in a largescale photobioreactor for microalgae cultivation. The dimensions of a new photobioreactor (width, height, and clearances) were optimized using the CFD model validated against experimental data. The bubble size was set to 5 mm in the CFD. Assuming that the microalgae was distributed uniformly and was well diluted, the solid phase was ignored [76]. Kim et al. (2017) [75] investigated the effect of the center of gravity (CoG) position on the CO2 removal efficiency for an amine absorber with structured packing subject to pitching motion. A gasliquid transient Eulerian 3D CFD model was used. The sliding mesh method observing at the outside of the amine absorber was implemented for the pitching motion [75]. The structured packing was assumed as a porous medium, ignoring complex flow channels of the packing. However, porous resistance and liquid dispersion forces were added into the NS momentum equation in the porous medium zone to mimic key features of the packing [54]. Figure 3a shows the layout of the three CoGs used for comparison in the amine absorber (0.125 m in diameter and 4.325 m in height). The equipment motion is the least severe when the equipment is located as close as possible to the CoG. When the column was inclined to the left, the liquid holdup of CoG3 in the middle of the absorber (h = 2.1 m) was strongly skewed to the left wall side (see Figure  3b). The CFD study successfully identified the effect of the installed location on the hydrodynamics of the process under ship motion. The hydrodynamics of an air-water-oil separator under three angular ship motions (pitch, roll, and yaw) were also presented recently [77].

Gas-Liquid Flow in Random Packing
Kang et al. (2019) [74] presented an unsteady-state 3D VOF-CFD that included a surface tension force (i.e., CSF) to investigate gas-liquid flow behaviors in random packing with 80 Rasching rings (inner diameter 4 mm and outer diameter 6 mm). The VOF-CFD results were compared with experimental data such as the gas-liquid interfacial area, liquid holdup, and pressure drop. The CFD study represented well the hydrodynamics of the small-scale packing column [74]. However, a VOF-CFD approach that considers the detailed geometry of packing particles may be prohibitive for a large-scale industrial column, as mentioned earlier. The microscale VOF-CFD simulation may be useful for developing correlations that are required in a macroscale CFD simulation for industrial packing columns [74].

Gas-Liquid Bubble Columns
Bubble columns have been widely used in chemical, pharmaceutical, petrochemical, biochemical, and metallurgical processes for gas-liquid and gas-liquid-solid contact or chemical reactions due to good heat and mass transfer [25,53,78,88]. CFD has become an important approach in understanding the hydrodynamics of multiphase flow, bubble-size distribution, and gas-liquid interfacial area, thus facilitating the design, scaling-up, or optimization of bubble columns [25]. Many researchers developed CFD models with various constitutive equations or sub-models for bubble columns.  [78] presented a transient 3D EM-CFD without the population balance equation (PBE) to describe bubble breakage and coalescence. The drag model was modified considering bubble swarms, and the mass transfer of O2 from the gas to liquid phase was considered via the linear driving force model. The CFD model was validated against experiments (radial gas holdup and radial liquid velocity). In the CFD study, the use of a single bubble size was a gross simplification. The impact of bubble swarms and wakes on the drag force must be better understood. Bubble-induced turbulence [79] may improve the level of agreement between modeling and experimental results [78].
A wide range of bubble sizes and shapes exist in bubbly flows, which leads to different transport characteristics. Sharma et al. (2019) [79] proposed turbulence-induced bubble collision diffusion for two-group bubble simulations. The bubble collision force for the two groups of bubbles was added to the gas-phase momentum equation of each group [79]. The gas-liquid Eulerian CFD model without PBE was applied to a turbulent flow regime. Bubble interaction mechanisms such as bubble breakup and coalescence would be useful to better predict hydrodynamics in bubble columns. Tran et al. (2019) [53] presented a gas-liquid Eulerian CFD model coupled with a PBE for an airwater homogeneous bubble column under high pressures. The drag coefficient and breakage kernel were modified to consider the impact of bubble swarms and pressure on the hydrodynamics of the bubble column. The gas holdup and mean bubble size obtained from the modified CFD model were compared with experimental data. Figure 4 displays the gas holdup contours in the air-water bubble column at 1-35 bar. The gas holdup increases as the pressure increases. It is clearly shown at the bottom of Figure 4 that the radial gas holdup is higher at higher pressure.  [80] proposed three drag models used for a gas-liquid transient Eulerian CFD model without PBE in an air-water bubble column under elevated pressure. The drag models covering the flow regime from homogeneous to turbulent included a pressure correction term to represent the effect of pressure on the gas holdup [80]. Gas-organic liquid bubble columns at high pressures and temperatures are common in industrial applications [89], therefore attention has been paid to the CFD simulation of those bubble columns. However, prior to CFD simulations for highpressure and high-temperature bubble columns, reliable drag, breakage, and coalescence models have to be developed that consider microscopic bubble behaviors [90].
Through the stability condition at a mesoscale flow structure, Yang and Xiao (2017) [25] presented an approach based on the EMMS concept to calculate correction factors of the drag coefficient (CD) and the coalescence rate. A CFD-PBE model coupled with the EMMS concept was solved for air-water bubble columns at various gas velocities and compared to experimental data such as the axial bubble size and bubble size distribution. The CFD-PBE model considering the mesoscale bubble structure showed improved prediction of the experimental data over the original CFD-PBE [25].

Gas-Solid Phase CFD Simulations
Fluidized-bed combustors (FBCs) such as bubbling and circulating FBCs have been used for power production from coal, biomass, and organic wastes owing to excellent heat and mass transfer, low pollutant emission, and high combustion efficiency [91,92]. Fluidized-bed chemical looping combustion (FB-CLC) has shown promise as a solution for high-efficiency low-cost carbon capture from fossil-fueled power plants [93]. 2D and 3D Eulerian CFD simulations with chemical reactions of the FB-CLC were performed to understand hydrodynamics coupled with chemical reactions [93].
Fluidized beds can be formulated using the two-fluid CFD model where the solid phase is considered to be a continuum and fully interpenetrating with the gas phase. A Eulerian-Eulerian CFD framework involving a granular solid and a gas phase is often used for computational investigations of fluidized beds [91]. The Eulerian CFD model includes the governing equations (the continuity equation, NS momentum equations, and energy equation) and the constitutive equations for the solid-phase pressure and viscosity provided by the KTGF [93], as listed in Table 4. Zhou and Wang (2015) investigated mixing and segregation in circulating fluidized-bed (CFB) risers using the KTGF-CFD model with the EMMS drag model. The EMMS approach showed improved hydrodynamics over homogeneous drag forces compared to the experimental data [58]. Fotovat et al. (2017) [43] presented Eulerian and Lagrangian CFD models to investigate the electrostatic effects on the hydrodynamics of fluidized beds. The CFD models assumed constant particle charges without allowing for the dynamic nature of tribo-electrification. Thus, a multiscale approach where charge generation and dissipation mechanisms at the particle scale are coupled with the distribution of charged particles at the mesoscale was suggested for the future work to predict the effect of electrostatic charges on the overall reactor performance [43].
The computational particle fluid dynamics (CPFD) model as a kind of Eulerian-Lagrangian method has been applied to gas-solid fluidized beds [92]. The solid phase is modeled by the Liouville equation for a particle distribution function (f) in CPFD [94,95] instead of the KTGF equation in Eulerian CFD. The solid properties are calculated from the f of representative computational particles, which can increase the computational efficiency for dense solid phases [92]. Kim et al. (2020) [96] presented a CPFD coupled with the PBE, which was applied to a multiphase and multiscale crystallization. Wu et al. (2020) [92] developed a cluster-based drag model to consider the multiscale characteristics of particles in a CPFD simulation. Yang et al. (2001) [97] investigated numerically and experimentally the bubble formation in a gas-liquid-solid fluidized bed at elevated pressures. The 2D VOF-CFD successfully predicted the bubble size with respect to the pressure and solid fraction using the liquid density, viscosity, and surface tension adjusted with the pressure. The VOF-CFD was performed in a small area (8 cm width and 5 cm height) with two orifices at the bottom [97]. Hamidipour et al. (2012) [65] presented an unsteady-state 3D EM-KTGF-CFD simulation of gasliquid-solid fluidized-beds. The CFD results obtained from various turbulence equations were compared with the experimental data such as the axial solid velocity and gas holdup. This study indicated that the constitutive equations and numerical schemes involved in the CFD model should be carefully selected to accurately capture the flow patterns. Pan et al. (2016) reviewed CFD simulations of gas-liquid-solid flow in fluidized beds [23]. A three-phase fluidized bed is intrinsically multiscale in time and space, including the macroscale at the reactor level, the mesoscale with particle clusters and bubble swarms, and the microscale with particle and bubble size [23]. Hydrodynamic properties such as the bubble size, bubble-rising velocity, liquid velocity, and gas and solid holdups were examined according to operating conditions in threephase fluidized beds. The authors suggested additional work required for three-phase fixed-bed CFDs: i) new drag models for interphase momentum exchange, ii) CFD modeling including chemical reactions and heat and mass transfers in industrial-scale reactors, and iii) multiscale modeling on the balance between computation cost and accuracy [23].

Multiscale CFD Simulations
The critical task and major difficulty in multiscale modeling is how to exchange information at the interface of neighboring sub-domains [98]. The applications of the multiscale CFD simulation are introduced in terms of concurrent and sequential coupling strategies.

Concurrent Coupling Strategy
Chen et al. (2013) [98] presented a multiscale modeling framework combining a gas-phase Eulerian CFD and a D2Q5/D2Q9 LBM to predict electrochemical transport phenomena in the cathode of PEM fuel cells. At each time step, the density, velocity, and concentrations were interchanged between the CFD and LBM domains, including the gas channel, porous gas diffusion layer, and catalytic layer. The results obtained from LBM were used as the boundary conditions of the CFD domain. At each time step, the CFD and LBM simulations were repeated until the two results converged [98]. Crose et al. (2017) [13] developed a multiscale CFD framework in the plasma-enhanced chemical vapor deposition (PECVD) of thin-film solar cells. A macroscopic gas-phase transient 2D CFD model with plasma chemistry and transport phenomena was solved to calculate the pressure, velocity, temperature, and species concentration. The CFD model was concurrently coupled with a kinetic Monte Carlo scheme describing the microscopic thin-film growth. Indeed, decoupled CFD and surface interaction models were unable to capture the spatially nonuniform deposition of silicon films. The multiscale CFD was capable of predicting the codependent behavior of the macroscopic gas phase and the microscopic thin-film growth rate. Using three user-defined functions (UDFs), including gas-phase reaction rates, an electron density profile in the gas phase, and a hybrid kinetic Monte Carlo (KMC) algorithm on the wafer surface, the macro-and microscale models were coupled at each time step [13]. The multiscale computing was parallelized within a message-passing interface (MPI) architecture to reduce the computational time.
Pozzetti and Peters (2018) [99] proposed a multiscale DEM-CFD method for the simulation of three-phase flows using a dual-grid method for the bulk and fluid fine scales, which is different from the classical CFD-DEM approach. The classical CFD-DEM approach based on the volume averaging technique aims to couple CFD with a DEM solution at a single scale using cell values or locally interpolated values to evaluate the fluid-particle interactions. The concurrent coupling between CFD and DEM can lead to high computational cost. The dual-grid multiscale DEM-CFD method aimed to capture the most significant phenomena of the two different scales while reducing the computational burden as much as possible. However, the new method required a complex interpolation between the two scales [99].
A heterogeneity concept for mesoscale particle-cluster behaviors was introduced into Eulerian CFD computational cells to modify the interphase drag coefficient in the fluidized bed. Qi et al. (2007) [100] presented EMMS-based multiscale drag models, where flow variables such as the voidage () and the gas and solid velocities (ug and us, respectively) obtained from CFD were used as input values to the EMMS model. The drag force (FD) considering the heterogeneity of the mesoscale was simultaneously updated for each computational cell and flow time [100]. Li et al. (2018) [87] proposed a cluster structure-dependent (CSD) drag model by introducing convective and temporal accelerations of both gas and particles in the dilute and dense phases. Although multiscale drag models showed better accuracy than common single-scale drag models, they required greater computational time.
The concurrent coupling strategy of the gas-solid Eulerian CFD with the multiscale drag models is shown in Figure 5. In each computational cell, an EMMS or CSD model is solved with the input variables (, ug, and us) from the CFD. The drag force (FD) computed by the EMMS or CSD model is returned to each cell.

Sequential Coupling Strategy
Raynal and Royon-Lebeaud (2007) [101] developed a typical sequential multiscale CFD approach from small to large scales of structured-packing columns for achieving an optimum column design. A 2D VOF-CFD model was used to calculate the liquid holdup and velocity at a small scale of the structured-packing. A gas-phase 3D CFD model with a moving wall boundary condition was applied to an REV of the structured-packing to calculate the relationship between the pressure drop and the gas superficial velocity. Once the liquid holdup, gas-liquid interaction, and pressure drop were determined at the small scales, a gas-solid 3D porous CFD simulation was performed for the whole absorber to examine the influence of internals (e.g., gas and liquid distributors) on the hydrodynamics of the absorber. As a direct CFD simulation that completely considers the complex geometry of structured packing requires significant computational time, the three-step multiscale approach is useful for investigating the hydrodynamics of industrial-scale absorbers. Qi et al. (2017) [102] proposed a sequential multiscale CFD method to predict the liquid distribution in Mellapak 350Y structured packing. A microscale VOF-CFD model was solved for six REVs at different positions of the packing. Using a macroscale unit network model (UNM), the liquid distribution at the macroscale was evaluated from the distribution coefficient computed by the microscale VOF-CFD model. The macroscale UNM assisted in the design and optimization of structured packing [102]. Ngo et al. (2017) [29] presented the multiscale CFD model of an impregnation die for carbon fiber (CF) prepreg production. A liquid-phase Eulerian CFD was used at the microscale tow domain to obtain resin permeability (K). A liquid-phase steady-state Eulerian CFD with K was applied to the porous media tow domain at the macroscale. Figure 6 shows a sequential coupling strategy between the two Eulerian CFD models. The resin penetration velocity (up) at the two tow scales was interconnected, therefore several iterations were necessary to obtain a converged K value. Ngo et al. (2018) [33] proposed another multiscale simulation approach to predict the degree of impregnation in the CF prepreg production process (see Figure 1). The multiscale simulation included three scales: a 2D unsteady-state VOF-CFD for an REV of the microscopic tow domain (see Figure 1d), a liquid-phase Eulerian 3D CFD for the entire impregnation die (see Figure 1b), and a process-scale simulation combining data from the micro-and macroscale CFDs. Figure 7 shows a sequential coupling strategy between the three scales.  In gas-solid fluidization systems, several multiscale concepts for the Eulerian CFD model have been proposed: a filtered drag force (FDF) model using constitutive equations from highly resolved subgrid models [103], a mesoscale structure-based drag model, where the flow structure and stability conditions were considered to capture the clustering behavior of solid particles using the EMMS model [46,73,104,105], and a CSD drag model [87,100]. Lu et al. (2019) and Wang (2020) [46,47] reviewed several sequential coupling strategies between the EMMS drag and Eulerian CFD models to speed up the multiscale simulation using a lookup table [105] or a fitting function [104] of the drag force. Figure 8 shows the sequential coupling strategy of the multiscale Eulerian CFD with EMMS in fluidized beds. The flow regime with superficial velocity ( ) and solid flux ( ), solid material property with particles density ( ) and diameter ( ), and a broad range of flow conditions with the Reynolds number ( ) and voidage ( ) are used as the input values of the EMMS model for generating a fitting function of the heterogeneity index (HD) that is directly multiplied by the original Wen and Yu drag force [104].

Macro-scale steady
Eulerian CFD

Discussion of Multiscale CFD Simulations
In recent CFD studies, multiphase, multiphysics, and multiscale simulations with complex geometries were performed with the increase in computational capacity. Using CFD models validated with experimental data, efforts to replace the experiments are continuing to reduce the time and cost required during process development. However, current CFD simulations are still far from meeting practical demands for process development.
The main challenges limiting the use of CFD simulations are accuracy and computational cost [17]. Despite the successful spread and usage of CFD simulations, there are potential risks of its misuse because of poor model choice, bad assumptions, and wrong interpretation of the results [39]. Reliable theories and numerical methods are needed for accurate simulation. Particular emphasis is placed on the validation of CFD models via experiments and suitable analytic solutions. Extensive validation with experimental data is still required for closure models under a wide range of flow conditions [51]. Table 6 summarizes the recent progress of multiscale CFD simulations in the concurrent and sequential coupling frameworks. Macroscale CFD models have been coupled with LBM, KMC, DEM, EMMS, CSD, and PBE at smaller scales to more accurately calculate physical phenomena at the cost of computational time. In the sequential coupling approach, model parameters required in the macroscale CFD were evaluated at smaller scales or vice versa. Many current works on multiscale modeling focuses on concurrent coupling [4]. Several libraries have been presented to facilitate the concurrent coupling effort for independently developed solvers. Bernaschi et al. (2009) [106] presented a multi-physics/scale code (MUPHY) based on the combination of a microscopic MD with a macroscopic LBM, using MPI as the communication interface for parallelization. Tang et al. (2015) [107] reported a multiscale universal interface (MUI) which is capable of creating an easily customizable framework for solver-dependent data interpretation. Neumann et al. (2016) [108] developed a macro-micro-coupling tool (MaMiCo) for the modulation of molecular-continuum simulations in fluid dynamics. The computational domain was split into a continuum region and a molecular region. At the interface between the two regions, flow quantities such as density, velocity, and temperature were exchanged and matched. MaMiCo supports 2D and 3D simulations, interface definitions for arbitrary MD and CFD codes, MPI-parallel execution, and time-dependent simulations [109].
In addition to multiscale CFD simulations, the mesh quality in 2D or 3D geometries is a prerequisite for achieving accurate CFD simulations. An open-source package, OpenFOAM, based on FVM, has been widely used in CFD because one can customize solvers and easily extend the numerical libraries [110]. NASA defined digital twins (DT) as an integrated multiphysics, multiscale, probabilistic simulation of a system that uses the best available physical models and sensor updates to mirror the life of its flying twin [111]. Simulation tools such as CFD are a key aspect related to DT [112]. This section presents the current research on mesh independency, OpenFOAM, and DT.

Mesh Independency
CFD models are solved numerically on meshes discretized in a 2D or 3D computational domain. The quality of the mesh plays an important role in the accuracy and stability of the numerical computation [22,113]. The discretization error vanishes only if the grid spacing tends toward zero [114]. A proper number of meshes should be determined as a compromise between computation efficiency and accuracy [68]. Thus, a mesh independence test should be performed in the verification step of the CFD calculation [115].
In a simple mesh independence test, three mesh systems-coarse, medium, and fine-were used [53,54,68,75]. The mesh number of the finer mesh system was double or triple that of the coarser mesh system. To ensure good mesh quality, the minimum orthogonal index was more than 0.01, and the average value of the orthogonal index was significantly higher than 0.01. The maximum skewness index was maintained below 0.95. In addition, the maximum aspect ratio was less than 100, which is much less for flow exhibiting strong gradients [22]. Hydrodynamic values such pressure, velocity, phase fraction, and temperature were compared with respect to the three mesh systems. The medium mesh system was selected, considering the accuracy of the CFD solution and the computation cost [22,54,68,116].
Based on the Richardson extrapolation theory, Roache (1994) [114] proposed a grid convergence index (GCI) to evaluate the numerical uncertainties of the CFD model. The GCI between two grid systems was evaluated. When the ratio of the two GCIs was close to unity, a proper mesh system was selected [115].
For gas-solid fluidized beds, Cloete et al. (2016) [117] investigated the grid independence behavior via the particle relaxation time for several common drag models. The practical rule for the hydrodynamic grid size (Δ ) expressed by a function of slip velocity ( ) was found as [117] Δ ≈ 0.0064 The required grid size for a reactive fluidized-bed (Δ ) was also proposed [117]. are the gas viscosity, particle diameter, solid-phase density, gas-phase density, bubble-to-emulsion mass transfer coefficient, and Archimedes number, respectively.

Open Source Code
In an open-source CFD such as OpenFOAM (open source field operation and manipulation), the user can access the source code and customize it to suit the user's needs [71,118]. Bhusare et al. (2017) [88] investigated the effect of internals on hydrodynamics in a bubble column using a gas-liquid Eulerian CFD model in OpenFOAM, which was compared with experiments and ANSYS Fluent (as a commercial CFD code). Kone et al. (2018) developed an OpenFOAM code to calculate a single-phase 3D CFD model in a PEM fuel cell [71]. Da Rosa and Braatz (2018) [34] presented a multiscale CFD model combined with a micromixing, energy balance and PBE in a tubular crystallizer using OpenFOAM. Farias et al. (2019) presented a liquid-solid Eulerian CFD model coupled with KTGF and PBE to simulate a lovastatin/methanol crystallizer using OpenFOAM [49]. The number of publications using OpenFOAM is increasing.

Digital Twins with CFD
A reliable CFD model is required to obtain meaningful results from the CFD [26], as mentioned previously. Wang and Xiao (2016) [119] proposed a data-driven CFD modeling of turbulent flows. The drag force (FD) for a porous medium domain, which was involved in a liquid-phase Eulerian CFD model, was inferred by assimilating experimental data (e.g., velocity and turbulence kinetic energy) and the CFD model prediction in an iterative ensemble Kalman method as an inversion approach [119]. If the experimental data are measured in real-time and provided to the CFD, a virtual machine realized by the CFD can be synchronized with the real machine (experimental apparatus). The main feature of DT is the real-time communication between the real and virtual machines, providing reliable CFD results. However, the real and virtual machines may be different from each other because errors can come from many sources, including a numerical discretization error or a modeling error. The computational cost increases owing to the iterative ensemble method [119]. Lu et al. (2020) [111] reviewed the recent development of DT for smart manufacturing in the context of Industry 4.0. DT, acting as a mirror of the real world, provides a smart tool for simulating, predicting, and optimizing physical manufacturing processes. Studies on standards, communication protocols, time-sensitive data processing, and the reliability of models need to be priorities for the next stage of research in DT [111].

Conclusions
This review explored exclusively Eulerian CFD simulations of chemical and biological processes at the multiscale in time and space. To shed light on the progress of CFD technologies in the past decade, CFD theories and their applications were reviewed in single and multiple phases. The following are the key conclusions from the review: 1) Chemical and biological processes have multiscale, multiphysics, and multiphase features in nature.
2) The Eulerian CFD model considering fluid as a continuum includes the volume of fluid (VOF) model, the Eulerian multiphase (EM) model, and the kinetic theory of granular flow (KTGF).
3) The Eulerian CFD has been used for designing, scaling up, and optimizing at the processlevel. Moreover, the Eulerian CFD is suitable for investigating the hydrodynamics (e.g., pressure, velocity, and volume fraction) of a process following geometrical and operational modifications. 4) RANS-based turbulence equations such as k- two-equation types are often used as a generalpurpose design tool. The direct numerical simulation (DNS) and large eddy simulation (LES) may be impractical for most industrial flow conditions because of computational time.
5) The Eulerian CFD has been applied to chemical and biological processes such as gas distributors, combustors, gas storage tanks, bioreactors, fuel cells, random and structured-packing columns, gas-liquid bubble columns, and gas-solid and gas-liquid-solid fluidized beds.
6) The critical task in multiscale modeling is to exchange information at the interface of neighboring sub-domains. The multiscale CFD strategy interacting between micro-and macroscale domains can be divided into concurrent and sequential couplings. Concurrent coupling is a powerful approach, but it is more computationally prohibitive than the sequential approach. 7) In the concurrent coupling framework, the Eulerian CFD at a macroscale is solved simultaneously with a lattice Boltzmann method (LBM), kinetic Monte Carlo (KMC), or discrete element model (DEM) at a micro-or mesoscale, updating information at each time step. In the sequential coupling framework, several parameters required for the Eulerian CFD model are calculated from a VOF or EMMS model at a micro-or mesoscale. When the model parameters are interconnected at the two scales, several iterations are necessary to obtain a converged parameter value.
8) Modeling, meshing, physical properties, and selection of a suitable turbulence model play an important role in CFD simulations. In particular, the proper number of meshes should be determined as a compromise between computation efficiency and accuracy. The mesh independence test is often performed to determine a proper number of meshes in the verification step of the CFD calculation. 9) Open-source packages such as OpenFOAM are widely used in CFD because one can customize solvers and easily extend the numerical libraries.
10) CFD models have to be validated by comparing the simulation results with experimental data or suitable analytic solutions to provide meaningful results. If the experimental data are measured in real-time and provided to the CFD, a virtual machine realized by the CFD can be synchronized with the real machine (experimental apparatus). The main feature of digital twins (DT) is the real-time communication between the real and virtual machines, providing reliable CFD results. DT acting as a mirror of the real world provides a smart tool for simulating, predicting, and optimizing physical processes.
11) The multiscale methods are attractive for industrial applications. However, substantial efforts in physical modeling and numerical implementation are still required before their widespread implementation. A powerful parallelized computational capacity is needed for achieving multiscale simulation and DT.