Immersed boundary method for multiphase transport phenomena

2.3.1 Momentum transfer boundary conditions

For momentum transfer at the boundary surface, the no-slip and no-penetration boundary conditions are usually enforced on solid surfaces. Those two boundary conditions indicate there are no relative velocities at the surface (Day 1990), mathematically u|Γ=0 and v|Γ=0, where u|Γ and v|Γ are relative velocity in the tangential and perpendicular directions at the interface Γ.

Another extensively used boundary condition for momentum transfer is the slip boundary condition (Kempe et al. 2015). In the cases of rarefied-gas dynamics and non-Newtonian fluids in which the Newtonian law of viscosity cannot hold, modeling the solid surface as slip boundary would be more physically accurate (Lauga et al. 2007). Investigation of flows over such slip surface is quite important as non-Newtonian fluid is frequently encountered in chemical and process engineering (Chhabra and Richardson 1999). Another example is in the numerical study of bubbly flows without any contaminants in a clean system, in which it is usually more appropriate to assume that the surface boundary condition of bubbles is free-slip (Kempe et al. 2015). For slip boundary condition, the tangential slip velocity is determined by the surface shear τw and the wall-normal stress σ _w as u|Γ=f(τw,σw) (Rao and Rajagopal 1999). The wall-normal velocity v|Γ is not affected by the slip boundary condition and is still zero for the solid impermeable surface.

2.3.2 Energy transfer boundary conditions

For the energy transfer at the interface, additional thermal boundary conditions usually need to be specified. The most extensively used thermal boundary conditions in simulations of multiphase flows are the isothermal and isoflux boundary conditions. The former refers to a Dirichlet boundary condition of T|Γ=T0, and the latter refers to a Neumann boundary condition of ∂T∂n‾|Γ=Φ. When Φ=0, the heat flux across the interface is zero. This is the so-called adiabatic boundary condition.

For more realistic consideration, the heat conduction in the solid side can also be taken into account. In the previous work of Iaccarino et al. (2002), it was shown that heat transfer in heated ribbed passages was significantly modified when the heat conduction through the passage walls was included. Therefore, to include the effect of heat conduction, the thermal boundary conditions should be enforced on both sides of the interface. This kind of heat transfer problem is called as the conjugate heat transfer (CHT) (Perelman 1961). The boundary conditions for CHT are Tfluid=Tsolid and λfluid∂Tfluid∂n‾=λsolid∂Tsolid∂n‾, which indicate that the continuity of temperature field and balance of heat flux on both sides are associated at the interface. Essentially, the boundary condition for CHT problems is a combination of the Dirichlet type and the Neumann type of boundary conditions. Due to the constraint of temperature and temperature gradient at the boundary, these two boundary conditions must be decoupled first in numerical implementation, which will be discussed in Section 4.

2.3.3 Mass transfer boundary conditions

For multiphase flows with mass transfer, the boundary conditions for two typical types of mass transfer problems are discussed. The first type is the mass transfer among a membrane without heterogeneous surface reaction or permeation problem (De Jong et al. 2012a; De Jong et al. 2012b; Miyauchi et al. 2017; Takeuchi et al. 2019; Yao and Mori 2017), and the second type is the mass transfer associated with heterogeneous surface reaction (Lu et al. 2019a; Qu and You 2019; Yousefzadeh and Battiato 2019; Zhu et al. 2019).

For the permeation problem, the mass transfer process is governed by the convection and diffusion process across the membrane interface. At the permeable interface, the solute flux in the normal direction is simply related to the concentration difference across the interface as [c]Γ=cΓ+−cΓ− (Gong et al. 2014). Thus the solute flux is expressed as J⋅n‾|Γ=−P[c]Γ, where P is the permeable coefficient, which is determined by the properties of the membrane. The diffusive fluxes J across the interface on both sides of the membrane are identical due to mass conservation. For more details of calculating the flux of solute, it is recommended to refer to the Kedem–Katchalsky relations (Kedem and Katchalsky 1958).

For the second mass transfer problem we consider here, the existence of surface heterogeneous reaction would strongly couple the momentum, energy, and mass transfer at the interface. Heterogeneous reaction introduces Stefan flow, which makes the problem more complicated. In fully resolved simulations of interphase mass transfer with the heterogeneous surface reaction, the convection and diffusion of the gas phase are balanced by the reaction rate at the interface as (Kee et al. 2005),

where V _k , Y _k , ωk, and W _k are the diffusion velocity, mass fraction, reaction rate, and molar mass of species k, and u _s is the velocity of the Stefan flow. The velocity of Stefan flow could be obtained by considering the mass balance equation

The sum of diffusion fluxes of all species at the surface is zero, thus giving

where K _g is the total number of species. At the surface, the heat balance contains the convective and diffuse flux of the gas phase, the reaction heat, the radiation, and the heat conduction across the interface as

where h _k is the enthalpy of species k, σ and ε are Stefan–Boltzmann constant and emissivity, respectively. Considering the computational efficiency, some of these parameters or terms in the equations may be simplified or ignored. These will be further discussed in Section 4.

Another noteworthy boundary condition related to the mass transfer problem at the body surface is that of a reactive flow over an embedded body without heterogeneous reactions. Under this condition, there is no surface reaction and mass transfer, thus only no-penetration condition for chemical species should be specified (McGurn et al. 2013), which is ∂Y∂n‾|Γ=0 at the boundary surface.

3.1.1 The original immersed boundary method

In the pioneering work of Peskin (1972), the embedded body surface is represented by a series of Lagrangian points. At each Lagrangian point, the flow exerts an elastic force f to the body, which is calculated under an elastic forcing model like Hooke’s law. This elastic force f is then spread backward to the surrounding Euler grid nodes by a Dirac delta function. Thus the interphase coupling conditions are represented by an additional forcing term F added to the right-hand side of governing equations in this artificially created area near the embedded body. By introducing the additional force F, the no-slip boundary condition at the boundary surface could be satisfied. We term this kind of IB method that spreads a Lagrangian force calculated from forcing models to surrounding Euler grids to exert boundary conditions as the original IB method. This method can be expressed in the incompressible Navier–Stokes momentum equation as follows,

where f(s,t) is the force at the s-th Lagrangian point X‾(s,t) at the interphase surface, δ is the Dirac delta function, μ is the dynamic viscosity, Γ is the immersed boundary, x‾ is the Eulerian coordinate. The formulation of the Dirac delta function is crucial. It would significantly alter the thickness of the artificial boundary as well as the accuracy of the original IB method. In the original study of Peskin (1972), the projection method was performed to achieve the continuity. The effect of the force term was only considered in the advection–diffusion step, which is shown as the following equation

where u∗ is the intermediate velocity and u ⁿ is the velocity at the time step n. No additional manipulations were required in the solution of Poisson solver to account for the presence of immersed objects.

One advantage of the original IB method is the consistency in handling the interior and exterior flow fields, i.e., the additional forcing F is added to the grid points on both sides of the boundary. This feature is essential for the case when the flows inside and outside the boundary interface are both of interest. The other advantage is the capability of treating deformable interfaces due to the elastic properties of the boundary specified in the original IB method. These merits make it especially suitable for simulations of the FSI problems. Similar ideas are applied to some other methods for simulations of deformable particles too, such as the front-tracking (FT) method (Unverdi and Tryggvason 1992) and the immersed-finite-element method (Saadat et al. 2018; Zhang et al. 2004).

The construction of the Dirac delta function that acts as a kernel in transferring information between Lagrangian and Eulerian fields is quite important and must follow several rules so as to give a smooth and accurate distribution of the additional force. Firstly, the discrete delta function should satisfy the zeroth and first discrete moment conditions so that the conservation of momentum and angular momentum could be achieved in interpolating the Lagrangian force to Eulerian grid nodes. Secondly, the discrete delta function should be continuous and vanish in far fields for computational efficiency. Based on these rules, Peskin (2002) introduced a four-point one-dimensional discrete delta function as ϕ(x/h)/h, where h is the grid width and ϕ(r) reads

Another popular ϕ(r) function is formulated as (Peskin 1977),

However, the rules for constructing discrete delta function are not unique. In the framework of locally refined and staggered grid system, Roma et al. (1999) modified these rules into weaker forms and proposed a discrete delta function with three-point width. Stokie (1997) considered the second discrete moment condition and proposed a new discrete delta function with better accuracy but more computational cost. Based on these rules, a lot of efforts have been made in comparatively studying the performance of specific delta functions and/or constructing a proper delta function with good performance (Bao et al. 2016; Delouei et al. 2015; Kang and Hassan 2011; Lee and Kim 2012; Yang et al. 2009).

Besides the original model developed by Peskin, many other models have been developed to calculate the surface Lagrangian forcing, such as the PID (proportional–integral–derivative) controller forcing model (Goldstein et al. 1993), the penalty IB method (Kim and Peskin 2007; Kim and Peskin 2009), and the volume penalty method (Angot et al. 1999; Kadoch et al. 2012). According to the classification of Mittal and Iaccarino (2005), the original IB method belongs to the continuous forcing method, where the force terms are firstly calculated according to some constitutive equations and then interpolated into the governing equation of the whole flow field before discretization. As a result, it is independent of the spatial discretization. On the contrary, the force term in the discrete forcing method is introduced only in the boundary region after the discretization of the governing equations, which will be discussed in the following sections.

3.1.2 The direct-forcing method

The direct-forcing method was originally derived by Mohd-Yusof (1997), but was first termed by Fadlun et al. (2000). It is also called as the discrete forcing method according to Mittal and Iaccarino (2005). In the direct-forcing method, the presence of the immersed object is represented through introducing additional forcing terms to the governing equations after numerical discretization. This forcing is direct in the sense that it is directly solved from the discrete momentum equations with the desired value of velocity imposed at the IB (Kajishima et al. 2001), rather than from constitutive models as applied in the original IB method. Therefore, the discretization scheme has a great influence on the accuracy, stability, and the conservation performance of the direct-forcing IB method.

The basic principle of the direct-forcing method is described as follows. For temporal advancement of the incompressible Navier–Stokes equation, one has

where u ⁿ is the velocity at the nth time step and uⁿ⁺¹ for the (n+1)-th time step, F^n+1/2 is the boundary force, RHS^n+1/2 contains convective, viscous, and pressure gradient terms at the right-hand side of the momentum equation at the intermediate time step between n and n+1 time step. RHS^n+1/2 could be represented by the prediction velocity u∗ as u∗=un+RHSn+1/2Δt. At the interphase boundary, to yield uⁿ⁺¹ = u _d (u _d represents the desired velocity at the boundary), the solution can be used backward to give

Under this formulation, the forcing term F^n+1/2 that directly compensates the error between the calculated velocity and the desired velocity on the IB is derived from the discrete-time Navier–Stokes equations. The fractional-step method is often used to enforce continuity. When calculating the intermediate velocity in the advection–diffusion step, the force term is involved (Uhlmann 2005), which is similar to the treatment in the original IB method. The k-th Runge–Kutta step is shown as

where α _k , γk, and ζk are coefficients of Runge–Kutta scheme and ϕ is the pseudo-pressure (Rai and Moin 1991). The force term f is zero in the original fractional-step method without IB.

The idea of calculating the direct-forcing term is easily understood, while the key issue is how to spread this forcing term to the Eulerian flow field due to the employment of non-body-conformal grids in the IB methods. Many schemes have been proposed to overcome this issue (Chiu et al. 2010; Fadlun et al. 2000; Kajishima et al. 2001; Kim et al. 2001; Luo et al. 2007a; Mohd-Yusof 1997; Silva et al. 2003; Zhang and Zheng 2007). However, it was pointed out by Uhlmann (2005) that the interpolation method applied by Fadlun et al. (2000) could lead to the spurious force oscillations due to insufficient smoothing (this feature is more related to the ghost-cell method and will be further discussed in the next section), which makes the method not suitable for particulate flows with moving boundaries. One important modification to the direct-forcing method is made by Uhlmann (2005) by combining the direct-forcing method with the original IB method. The basic idea is similar to the original IB method by representing the boundary by a series of Lagrangian points, each of which exerts a direct forcing to the surrounding flow fields both inside and outside the boundary. Similarly, a regularized delta function is used as the kernel for smoothly transferring quantities between Lagrangian and Eulerian grid points. Employment of this regularized delta function could effectively damp the pressure oscillations when simulating moving particles (Uhlmann 2005; Yang et al. 2009).

However, the velocity on the Lagrangian points at the interphase boundary may not adequately satisfy the no-slip boundary condition due to the mutual influence of the direct forcing at the neighboring points through the interpolation of the Dirac delta function in the direct-forcing method of Uhlmann (2005). After the additional force added, one desires to get un+1=ud but always gets un+1≠ud instead. Although un+1≈ud could be achieved after a long period, the no-slip boundary condition may be still unsatisfactory. To remedy this problem, Luo et al. (2007b) and Wang et al. (2008) proposed the multi-direct-forcing method to reduce the difference between uⁿ⁺¹ and u _d by iteratively updating the velocity un+Δt1 at the sub-time step n+Δt1 based on the difference between un+Δt1 and u _d . This iterative method can largely improve the numerical accuracy in imposing the no-slip and no-penetration interphase coupling conditions (Luo et al. 2007b). Wang et al. (2008) found out that three additional direct-forcing iterations could reduce the l₂-norm of the velocity at the surface Lagrangian points by nearly 90%. In the study of Kempe and Fröhlich (2012), they found that for the flow past a single cylinder, three additional forcing loops could reduce nearly 95% velocity error at the cylinder surface, with only 8% computational cost increased. Tang et al. (2014) tested the performance of different forcing loop numbers in simulations of flow past an array of spheres and found that two additional iterations could reduce nearly 95% of the error of the nondimensional drag force.

Both the direct-forcing method of Uhlmann (2005) and the original IB method of Peskin (1972) apply the Dirac delta function to spread the Lagrangian forcing. One significant disadvantage of using the Dirac delta function is that the effective particle diameter is increased, resulting in inauthentic boundary capturing (Breugem 2012). Using the Dirac delta function would also overestimate the drag force experienced by the fluid phase. To alleviate this issue, some studies proposed to retract the Lagrangian points slightly inward to correct the boundary size, and a better drag force and a more accurate flow field could be obtained (Breugem 2012; Luo et al. 2019; Tang et al. 2014).

Besides the direct-forcing method of Uhlmann (2005), some other popular direct-forcing methods are also briefly introduced here. Nakayama and Yamamoto (2005) developed a smoothed-profile method to model the particle phase as Eulerian fields of particle concentration and velocity, which are all zero except in the immediate vicinity of the IB. Based on these Eulerian fields, additional forcing term is added to the governing equation to represent the presence of the IB (Jafari et al. 2011; Mino et al. 2017). Tenneti et al. (2011) developed a particle-resolved uncontaminated fluid reconcilable immersed boundary method (PUReIBM), which successfully alleviates the smear of boundary by restricting the additional direct-forcing term only to the grid points inside the particles, while maintaining the solution of the Navier–Stokes equation uncontaminated in the fluid phase.

3.1.3 The ghost-cell method

The idea of the ghost-cell IB method originates from the studies of Mohd-Yusof (1997) and Fadlun et al. (2000). They pointed out that the explicit forcing terms in governing equations could be replaced equivalently by reconstructing the velocity field in the vicinity of the boundary region in the time advancement of the governing equations. Later on, Tseng and Ferziger (2003) further used the ghost cells inside the boundary to implicitly represent the direct forcing terms and named the method as the ghost-cell IB method.

The ghost-cell method owns the distinguishing feature of not requiring the explicit addition of forcing terms to the governing equations (Mittal and Iaccarino 2005). The values of the ghost cells are determined by interpolation of boundary conditions and mirror points. The mirror points, also called as image points in some literatures, are fluid points inside the boundary, whose values are interpolated from the surrounding Eulerian fluid grids. By directly solving the governing equations in the vicinity of the IB without any additional terms, the presence of the immersed object could be represented with the help of these reconstructed ghost cells. As for the more interior grid points (than the ghost cells) inside the boundary, their values could be generally left unsolved (Nagendra et al. 2014; Tseng and Ferziger 2003). In the time advancement of the fractional-step method, the detection of boundary interface and determination of values of ghost cells should be carried out before calculating the intermediate velocity. With the help of the ghost cells, the predicted velocity fields including those at the vicinity of the IB can be obtained. Therefore, the pressure Poisson equation could be solved straightforwardly, together with the update of the velocity fields.

The interpolation scheme is an essential component in the implementation of the ghost-cell method. Values at both the mirror points and the ghost cells are determined by interpolation. Many interpolation and extrapolation schemes are available in constructing the ghost-cell method, such as the linear interpolation (Gilmanov et al. 2003; Mittal et al. 2008), the bilinear/trilinear interpolation (Luo et al. 2008; Luo et al. 2016b; Pan 2010; Seo and Mittal 2011a), the high-order interpolation (Seo and Mittal 2011a; Xia et al. 2014), and even the wall function model (Bernardini et al. 2016). Utilizing a suitable interpolation/extrapolation scheme with higher accuracy and low computational cost is quite crucial in this method.

Pressure oscillation in moving/deformable boundaries is another frequently investigated issue for the ghost-cell methods. Just as commented by Uhlmann (2005), the interpolation procedures of Fadlun et al. (2000) relating values at fixed grid nodes and values at arbitrarily located boundaries can lead to spurious force oscillations due to insufficient smoothing. Some potential sources have been explored, such as implementation of the pressure boundary conditions at the IB (Mark and van Wachem 2008; Shinn et al. 2009), violation of geometric conservation law (Seo and Mittal 2011b), spatial discontinuity in pressure when a grid point inside the solid region becomes that of fluid and temporal discontinuity of velocity vice versa, namely “fresh cell” and “dead cell,” respectively (Lee et al. 2011; Seo and Mittal 2011b). Many researchers have paid attention to this problem of the ghost-cell method. Improvements include combining the ghost-cell method with a cut-cell IB method (Merlin et al. 2013; Seo and Mittal 2011b), introducing additional mass source/sink terms together with a larger time step in fully implicit numerical scheme (Lee and You 2013), correcting the fictitious velocity and pressure fields at the ghost cells to ensure no penetration and mass conservation through boundaries (Mark et al. 2013).

3.2.1 The cut-cell method

The cut-cell method is a successful application of the IB methods in the finite-volume framework. It was first developed by Udaykumar et al. (1997) and Ye et al. (1999). In this method, the cells near the boundary are cut and reshaped to fit the boundary curves, and the boundary conditions are imposed on the reshaped boundary interface. For the grid points inside the geometry, they could be just left unsolved if the interior flows are not of interest or solved based on specific problem of interest due to the convenience of the finite-volume solver (Ye et al. 1999). One of the most highlighted advantages of the cut-cell method is the strict conservation of mass and momentum in the vicinity of the embedded boundary due to the usage of the finite-volume solver (Schneiders et al. 2013). This makes the cut-cell method a good choice in maintaining a divergence-free velocity field (Muldoon and Acharya 2008) and damping the spurious unphysical pressure oscillation (Cecere and Giacomazzi 2014; Seo and Mittal 2011b) in incompressible flows.

However, in the generation of cut cells, exceptional cells of small size may exist, which could result in impractical small time step and numerical instability (Yang et al. 1999). Different techniques have been developed to handle this problem. Zeuuw (1993) used a special local time-stepping scheme to overcome the small time step problem. Ye et al. (1999) reshaped and merged cut cells to form a bigger cell. Muralidharan and Menon (2016) developed a novel cell clustering algorithm for polynomial reconstruction in the vicinity of small cells to achieve smoothness. Kirkpatrick and Armfield (2003) introduced a cell linking method to tag cells as master and slave cells to enable the same treatment for all cells.

When dealing with the time advancement, the boundary cells need to be reshaped into the trapezoidal boundary cells according to the intersection of the boundary with the Cartesian mesh, as shown in Figure 4. In most cases, the face centers of trapezoidal boundary cells are not located at the midpoint of neighboring cell centers. Therefore, to maintain the accuracy, some complicated polynomial interpolating functions may be employed instead of using linear approximation. Consequently, the convergence of the Poisson equation becomes more difficult as the stencil needed by the interpolating functions of the trapezoidal boundary cells cannot be included in the line-successive-overrelaxation (SOR) sweeps.

Figure 4:

The trapezoidal boundary cells and stencil points for interpolation to obtain the flux. Reprinted from Udaykumar (2001) with the permission of Elsevier.

In addition, there exist some inherent difficulties in employing the cut-cell method. The cut-cell generation is very challenging when dealing with complex three-dimensional problems since complicated cut-cell geometry may occur. To capture the boundary correctly with high efficiency, the cut-cell method is usually combined with the level-set method to identify the boundary interface (Cheny and Botella 2010; Schneiders et al. 2013). Another problem is that numerical oscillations will occur if the cut-cell method for fixed bodies is directly applied to solve moving boundary problems (Schneiders et al. 2013). To overcome this problem, Chung (2006) developed a fully implicit scheme for convective and diffusion terms integration along with a pressure-free projection method to handle moving boundary problems. Muralidharan and Menon (2018) combined the cell clustering technique with the cell mass redistribution for moving boundaries to ensure mass and momentum conservation near the boundary.

3.2.2 The immersed interface method

In the original IB method of Peskin (1972) and the direct-forcing IB method of Uhlmann (2005), the Dirac delta function would cause diffusive boundary problems, as discussed before. To alleviate this problem, Leveque and Li (1994) first proposed the IIM for elliptical equations to incorporate the forcing terms into jump conditions at the boundary interface and modify the numerical discretization algorithms, instead of spreading the influence of singular forces. In this method, the forcing terms that are a function of the jump conditions and their derivatives are derived from the Taylor expansions rather than approximating the Dirac delta function. This manipulation restricts direct modifications only to the boundary surface. This method was then extended to Stokes flows (LeVeque and Li 1997) and full Navier–Stokes flows (Lee and LeVeque 2003; Li and Lai 2001).

In the implementation of the IIM, the first step is to obtain the singular force acting on the interface. At the early stage, this singular force is calculated as that in the original IB method for elastic boundary or by PID-forcing model for rigid boundary (Goldstein et al. 1993). To avoid numerical instability due to the introduction of some constant parameters in these forcing models, Xu (2012, 2008, 2009, 2011) presented an explicit approach to decompose the singular force into tangential and normal components to enforce the prescribed motion of a rigid boundary in an incompressible viscous flow. Le et al. (2006), Tan et al. (2009), and Li et al. (2016) proposed an implicit method to calculate the singular force for the rigid boundary that shared the same basic idea with the direct-forcing scheme. In the work of Xu and Wang (2008), the singular force is calculated under the spring model or surface tension force model for different problems. Liu et al. (2019) modified the surface tension model to a strictly inextensible model, in which the surface tension was not calculated as an unknown variable but fully satisfied the inextensibility of the boundary surface.

After obtaining singular force acting on the interface, the next step of the IIM is to incorporate this singular force into the jump conditions and relevant derivatives. The relationship between the jump conditions together with their derivatives and the singular forces has been studied by a larger number of studies over recent decades (Lee and LeVeque 2003; LeVeque and Li 1997; Li and Lai 2001). Lai and Li (2001) derived the principal jump conditions across an immersed moving membrane of velocity, pressure, and relevant normal derivatives. Xu and Wang (2006) systematically derived necessary spatial and temporal jump conditions for simulating viscous flows subject to moving boundaries with the third-order and second-order accuracy near the boundaries. Xu and Pearson (2015) derived jump conditions for the triangular mesh system in three-dimensional simulations.

Finally, the discretization of governing equations in the vicinity of the immersed object, both inside and outside the interface, is modified implicitly with all these jump conditions and relevant derivatives. During the time advancement, the force term is not enforced only in the advection–diffusion step of the fractional-step scheme, compared with the original IB method. It is divided into tangential component F₁ and normal component F₂. The tangential component F₁ term is involved in the advection–diffusion step, while the normal component F₂ that indicates the proper jump condition is solved together with the Poisson equation (Lee and LeVeque 2003) as follows:

Advection–diffusion step:

Projection step:

To obtain high-order accuracy, more jump conditions of the variable and its high-order derivatives are favored. Zhong (2007) derived six versions of new IIM up to fourth-order accuracy by using only two jump conditions together with boarder grid stencils. Liang et al. (2008) and Jiang et al. (2012) combined the spectral scheme with the IIM to take advantage of the high accuracy of the spectral method. Hosseinverdi and Fasel (2017) used a compact difference scheme to derive the first- and second-order derivatives with seventh-order and sixth-order accuracy, respectively.

4.2.1 Natural convection

There are two major configurations in the investigation of the natural convection process by using the IB methods. One is the natural convection within an enclosure, and the other is the sedimentation process accompanied by natural convection.

Natural convection within an enclosure is substantial to many industrial and environmental applications such as heat exchangers, chemical reactors, cooling of electronic equipment, and refrigerator condensers. Due to the exclusion of external flow motion, the Rayleigh number is demonstrated to be a dominating factor in characterizing flow and thermal patterns and distribution of local Nusselt number. Kim et al. (2008) studied two-dimensional natural convection of a heated cylinder within a cold square enclosure by using the ghost-cell method. When the Rayleigh number was small at Ra = 10³ and 10⁴, the thermal field was fully symmetric, and four counterrotating vortices occurred, as shown in Figure 7. In the case with a stronger buoyancy effect, such as Ra = 10⁵ and 10⁶, the intense heat convection generated a thermal plume, which modified the flow field significantly, and only two recirculating eddy cores existed so that the dominant flow became upward. The thermal boundary layer at the cylinder surface also became thinner due to the strong convection. Similar phenomena were also observed in many other studies with different complex geometries (Liao and Lin 2012a; Park et al. 2013a; Park et al. 2013b).

Figure 7:

Isothermals (left) and streamlines (right) at four different Rayleigh numbers for the natural convection of a circular cylinder at the center within a square enclosure. Reprinted from Kim et al. (2008) with permission of Elsevier.

Recently, natural convection in a three-dimensional cubic enclosure immersed with cylinder has received much attention (Feldman 2018b; Lee et al. 2016; Seo et al. 2016; Spizzichino et al. 2019b). Comparing with the two-dimensional calculations, significant difference exists for the thermal and flow fields (Choi et al. 2015b; Pandey et al. 2019) as well as the surface-averaged Nusselt number (Choi et al. 2015a) in the three-dimensional simulations, even the overall computational parameters are similar. One typical and excellent three-dimensional study has been done by Seo et al. (2019a; 2019b). They demonstrated that the λ vortex structures (Jeong and Hussain 1995) in the natural convection of circular cylinder in a long rectangular enclosure are evidently three-dimensional and much more complicated than those in the two-dimensional cases, as shown in Figure 8.

Figure 8:

The lambda vortex structure (green isosurfaces) in the natural convection of circular cylinder in a long rectangular enclosure at Ra = 10⁶. Reprinted from Seo et al. (2019a) with the permission of Elsevier.

One interesting and notable feature of the natural convection within an enclosure is the unsteadiness. For an enclosure with simple symmetric bifurcated geometry, the flow field and thermal field are generally symmetrically bifurcated too. If the buoyancy effect is strong enough, even though the geometries of embedded object and enclosure are symmetric, the flow and thermal patterns could be asymmetric and unsteady (Cho et al. 2019; Cho et al. 2018; Choi et al. 2014). The asymmetric distribution of flow and thermal fields also leads to an asymmetric distribution of local Nusselt number at the enclosure surface due to the domination of convection (Park et al. 2013b). However, currently there is still no quantitative indicator for estimating the unsteadiness of flow and thermal fields, as this feature is influenced by geometries of enclosure and interior object, thermal boundary conditions, relative position, inclination angle, and aspect ratio of the embedded object, as well as primarily the Rayleigh number.

The other frequently studied problems are the sedimentation process of particles with natural convection. In these cases, the thermal flux between particles and fluid would modify the flow fields around particles under the influence of buoyancy force and then alter the motion of the particles. This combination effects would result in a huge change of particle settling speed compared with that without heat transfer or buoyancy effect. In the work of Wu et al. (2017), sedimentations of a hot, cold, and isothermal single particle in an infinite channel were studied by using the hybrid IB-LBM. Heat transfer of particle to fluid was proved to influence the buoyancy effect significantly, thus making the trailing vortices of hot settling particle disappear, as shown in Figure 9. It was reported that the hot particle owned the lowest settling velocity, which was only 50% of that in the isothermal case. While for cooled particle, the settling velocity is 40% more than that of the isothermal case. Xia et al. (2015) applied the high-order ghost-cell method to study the sedimentation of a hot spherical particle in the cool fluid. It was found out that for all the Reynolds number considered (Re = 50, 71, and 100), particle settling velocity considering buoyancy force (Grashof number Gr = 5000) was much smaller compared with those cases without buoyancy force (Gr = 0) at the same Reynolds number. The heat transfer and hydrodynamics interaction were strongest in the case of Re = 50 and Gr = 5000, which formed the highest Richardson number in their simulations.

Figure 9:

Dimensionless velocity fields of the fluid in the sedimentation of single hot particle (left) and cold particle (right) cases at Re = 40.5. The contours and the vectors represent the magnitude and the direction, respectively. Particles are represented by white spheres. Reprinted from Wu et al. (2017) with the permission of Elsevier.

Sedimentation of particle clusters has great importance in industrial applications. Compared with single-particle settling, sedimentation of clusters of particles becomes more complicated due to particle–particle interaction and hydrodynamics–thermal interaction. Feng and Michaelides (2009) studied the settling and cooling of 56 particles with the direct-forcing IB method. The particle–fluid density ratio was set as 1.5, Prandtl number as Pr = 1, Grashof number as Gr = 200. It was pointed out that there exists strong interaction between particle trajectory and heat transfer, and the thermal wake of leading particles could help the following particles keep higher temperatures.

Musong et al. (2015) studied the sedimentation of vertically and horizontally placed bispheres and triangular clusters at Pr = 0.72 and Gr = 100 with the direct-forcing method. They concluded that the orientation and distance between particles significantly affected the average Nusselt number and temperature patterns. For in-line placed two particles system, increasing the distance between particles would always increase the average heat transfer rate. For tangentially placed two particles system, the increase of distance between each particle would initially lead to an enhanced heat transfer rate. Further increasing would make the combined plume structure isolated, and the bispheres act like two single spheres. For multiple particles system, the average heat transfer rate per sphere particle always decreased with the increase of the number of spherical particles in a cluster.

Eshghinejadfard and Thévenin (2016) investigated the movements of 60 spherical particles under natural convection with a low solid–fluid density ratio of 1.1 at Pr = 0.7 and Gr = 1000 by using the direct-forcing IB-LBM. Particles were released at the bottom of the computational domain and firstly settled down and then quickly rose up to the top of the enclosure, as shown in Figure 10. Particle–particle and particle–fluid interactions generated strong vortical structures and made some individual particles move upward quicker than the single-particle case. Sometimes this strong vortex could even drive the particles to move back to the bottom of the enclosure.

Figure 10:

Temperature distribution of the center plane and particle position evolution at different times due to the natural convection in the enclosure; red and blue represent the higher and lower temperatures, respectively. Particles are represented by white spheres. Reprinted from Eshghinejadfard and Thévenin (2016) with permission of Elsevier.

In summary, multiphase flow with natural convection has been well studied using the IB methods. However, a systematic understanding of natural convection is still on the way. There is no criterion for determining whether the flow could eventually reach a steady state for the natural convection within an enclosure. Further investigations are required to explain the mechanism for the three-dimensional unsteadiness in natural convection within an enclosure. Also, the multiparticle sedimentation with heat transfer is far from clear. The distinct particle behaviors in the multiparticle system due to buoyancy force need further investigations (Eshghinejadfard and Thévenin 2016). Mechanism studies of the interaction between particle motion and interphase heat exchange also need more efforts.

4.2.2 Forced convection

Forced convection is usually encountered in gas turbines, combustors, fluidized beds, and other industrial processes. Here we mainly consider the following types of forced convection, namely the benchmark problems of uniform flow past single object with simple geometry, the heat transfer of flow past clusters, the investigations of heat transfer coefficient in a complex configuration, the turbulent heat transfer, and the heat transfer in fluidized bed systems.

Heat transfer of flow past single object with stationary or moving boundary has been extensively studied mainly as a benchmark problem for the validation of the thermal IB methods. A parametric summary of these studies is presented in Table 4. Most of these studies focus on low or moderate Reynolds numbers ranging from 10 to 200 with the direct-forcing method and the ghost-cell method, and all three types of thermal boundary conditions were considered. With the increase of the Reynolds number, the average Nusselt number is also increased due to a higher flow velocity and the resultant vortex shedding, thus leading to stronger interphase heat transfer. Considering the heat transfer rate among the three types of thermal boundary conditions, it is believed that the heat transfer coefficient is the highest for the Neumann-type thermal boundary condition (iso-heat-flux), while lowest for the Dirichlet (isothermal) boundary condition. The distributions of the local surface Nusselt number and the temperature of the heat transfer in the uniform flow past a single cylinder are shown in Figure 11. As can be seen that the surface temperature and Nusselt number in a Robin type case are smaller than those in a Neumann case (Luo et al. 2016b; Pan 2012). This is associated with the variation of surface temperature changing inversely with the local Nusselt number in a Neumann case (Zhang et al. 2008).

Figure 11:

The local distribution of the Nusselt number (A) and temperature (B) at the surface of the cylinder particle. Reprinted from Luo et al. (2016b) with permission of Elsevier.

Heat transfer in clusters is much more complicated compared to that of flow past a single object. The particle-averaged Nusselt number could be increased or decreased due to particle–particle interaction. Tavassoli et al. (2013) investigated a thermal flow of T∞=275 K past 54 randomly packed stationary sphere particles with a constant temperature of 320 K by using the direct-forcing IB method. The Reynolds number was chosen as Re = 10, 50, and 100, and the solid fraction as ϕ = 0.1, 0.3, and 0.5. The small changes in the temperature and velocity fields in the vicinity of particles were well captured, as shown in Figure 12. It was found that great discrepancy existed in the local Nusselt number for every single particle, but a good agreement of averaged Nusselt number with Gunn’s correlation (Gunn 1978) was obtained. Xia et al. (2014) studied the forced convection in a cluster of 20 randomly packed stationary sphere particles by using the high-order ghost-cell method. The solid fraction is ϕ = 0.16, and the Reynolds number is Re = 20, 50, and 100. They reported that the flow patterns and drag force coefficients were like those in flow past a single sphere particle with an equivalent diameter. Nevertheless, the particle-averaged Nusselt number of the cluster was much lower than that of an isolated sphere particle under the same equivalent Reynolds number, showing that particle clusters severely lowered the heat transfer of each particle.

Figure 12:

The temperature (A) and velocity (B) distributions of fluid phase in a randomly packed cluster domain (54 particles) at Re = 50 and ϕ = 0.1. Particles are represented by black spheres. Reprinted from Tavassoli et al. (2013) with permission of Elsevier.

Many studies focus on the correlations of local and averaged Nusselt number in particle clusters with forced convection heat transfer in order to establish more accurate multiphase heat transfer models. With the multi-direct-forcing method, Wang et al. (2009) derived a correlation for the local Nusselt number distribution at each column of the staggered cylindrical tube bank. The temperature at the cylinder surface was kept constant and the Reynolds number was set as Re = 20. The local Nusselt number distribution could be written as:

where Nu(θ,m) is the Nusselt number at the θ degree at the m-th cylinder’s surface. In the work of Tenneti et al. (2013), the heat transfer in a fixed cluster of sphere particles at constant temperature was studied by the PUReIBM. They concluded that for gas–solid systems with solid volume fraction ϕ < 0.4 and particle mean slip Reynolds number Re _m >10, the bulk fluid temperature would decay over a few particle diameters. However, in the case of ϕ > 0.4 and Re _m  < 10, particle heat transfer would significantly alter the bulk temperature. This result showed that the scale separation assumption in the modeling of the two-fluid method was not always valid. Later on, they studied the heat transfer in fixed assemblies of particles under the Dirichlet thermal boundary condition and proposed a new correlation of Nusselt number for volume fraction within 0.1–0.5 and particle slip Reynolds number within 1–100 (Sun et al. 2015). The newly proposed correlation was

where ε = 1−ϕ is the porosity. Based on this correlation, they modified the heat transfer model for the two-fluid model. He and Tafti (2017) investigated the heat transfer in the flow past a fixed random assembly of spheres and ellipsoid particles at low to moderate Reynolds numbers under Neumann thermal boundary condition by using the ghost-cell method. The Reynolds number was varied from 10 to 200. The Prandtl number and solid volume fraction were set as 0.2 and 0.74, respectively. The temperature variation at particle surfaces for both cases at Re = 200 was clearly shown in Figure 13. It could be seen that the temperature at the particle surface could vary significantly with the changing of the azimuthal angle, and the particle position could alter the surface temperature. The average Nusselt number of spherical particles was in good comparison with existing data. In general, the mean Nusselt number was larger for ellipsoid particles than sphere particles, and particle orientation had a significant influence on the modulation of the flow and thermal fields around ellipsoid particles. A new correlation for the ellipsoidal particles of aspect ratio 1:2.5 was given as follows:

Figure 13:

The dimensionless temperature variation at the surfaces of a fixed random assembly of spherical particles (A) and ellipsoidal particles (B) at solid fraction of ϕ = 0.2 and Re = 200. Reprinted from He and Tafti (2017) with permission of Elsevier.

This correlation was suitable when the Reynolds number and solid fraction were in the range of 10 ≤ Re ≤ 50 and 0<ϕ≤0.35, or of 50 < Re ≤ 200 and 0.1≤ϕ≤0.35.

Most of the aforementioned studies focus on heat transfer of particles under stationary status at a relatively low Reynolds number without turbulence. However, the actual particle motion and heat transfer processes in a fluidized bed and other chemical and industrial equipment are usually free and strongly interacted with turbulence. Hence, the study of fully resolved particles interaction with turbulence is imperative.

With the high-order ghost-cell method, Xia et al. (2017) studied the heat transfer of a group of free-moving particles subjected to incoming isotropic turbulence with a characteristic Reynolds number of Re = 271.2. The particle temperature was set as constant, and the solid volume fraction was set as 0.01 and 0.02 for two different cases. They found that the classical drag and Nusselt correlations for stationary particles, such as the Schiller–Naumann formula (Naumann and Schiller 1935) and Gunn’s correlation (Gunn 1978), could hardly predict those for moving particles. In all the cases considered, the particle-averaged Nusselt number was higher than that predicted by Gunn’s correlation (Gunn 1978). Compared to heavier particles, lighter particles were more pronouncedly affected by turbulent fluctuations so that the interphase heat exchange was much more significant, thus resulting in a higher average Nusselt number for lighter particles. Soon afterward, they carried out a fully resolved direct numerical simulation of spatially developing boundary layer laden with 1242 spherical particles at a characteristic momentum thickness Reynolds number Reδ2=800 by using the high-order ghost-cell IB method (Xia et al. 2018). The solid volume fraction was very low at ϕ=0.25%. It was pointed out that finite-size particles were strong thermal and kinetic energy sources to the turbulent boundary layer, thus leading to the increase in velocity fluctuations and wall-normal turbulent heat flux. However, since particle surface temperature was kept the same as the wall temperature, the existence of particles in the near-wall region damped temperature fluctuation. This could also be demonstrated in the isosurfaces of dimensionless temperature of 0.5, as shown in Figure 14. Compared with the no particle case in Figure 14A, the isosurfaces in Figure 14B are thicker and the temperature near the wall is more uniform due to the existence of particles so that temperature fluctuations in the near-wall regions are suppressed.

Figure 14:

Isosurfaces of the dimensionless temperature T^* = 0.5 in the turbulent thermal boundary layer for the no particle case (A) and multiphase case (B), colored by wall-normal height. Reprinted from Xia et al. (2018) with the permission of Elsevier.

Although the computational cost is very high, many studies directly focus on the heat transfer in fluidized bed systems aiming at revealing the complex multiscale heat transfer process by using the IB methods. Feng and Musong (2014) studied the heat transfer of 225 heat spheres in a fluidized bed using the direct-forcing method. Five cases with fluidization velocity varying from 0.035 m/s to 0.055 m/s were studied with a particle diameter of d = 0.00635 m. Characteristic hydrodynamic and thermal properties showed reasonable agreements with other literatures. It was found that the particle-averaged Nusselt number increased with the enhancement of the fluidization velocity and the bed height. This method was further extended to study the fluidization in a narrow channel with 1024 particles (Feng and Roig 2014). Fluidization velocity was set as 0.04 m/s and 0.045 m/s with the same particle diameter as that in the work of Feng and Musong (2014), and similar results were obtained. Deen et al. (2012) studied a pseudo-two-dimensional fluidized bed system with 1296 sphere particles with a diameter of dp=0.005 m by using the ghost-cell method. The minimum fluidization velocity was 1.17 cm/s, the terminal velocity was 39.6 cm/s, and the superficial liquid velocity was 8 cm/s. The flow field was highly uniform under the influence of particles, and significant bed expansion occurred due to the high superficial velocity, as shown in Figure 15. This high superficial velocity also led to the great spatial variation of fluid–particle heat transfer coefficients.

Figure 15:

The particle configuration in a subdomain of the fluidized bed (left) and the velocity profile of the fluidized agent in a partial enlarged view (right); blue spheres represent the fluidized particles while red and blue represent the higher and lower velocities, respectively. Reprinted from Deen et al. (2012) with the permission of Elsevier.

In summary, the IB methods, particularly the direct-forcing method and the ghost-cell method, have been extensively used in fully resolved simulations of multiphase flows with forced convection. The particles are generally assumed to be stationary in the study of particle clusters with heat transfer, and it has been demonstrated that the particle motion would significantly alter the heat transfer process as well as the Nusselt number (Xia et al. 2017). In addition, the distribution of local Nusselt number varies a lot in different locations at the particle surface, and the rotation of particles should be considered in the future work in deriving a general correlation for Nusselt number. The interaction of fully resolved particles and heat transfer in turbulent flow is also far from clear, and more investigations on particle–turbulence–heat transfer interactions should be carried out. Furthermore, most of the aforementioned studies are limited to small-scale system loaded with few particles. Fully resolved numerical simulations of multiphase heat transfer in a large-scale system laden with thousands of particles are on the way ahead.

4.2.3 Conjugate heat transfer

In the aforementioned studies, the effects of heat conduction within the solid were generally neglected. However, in realistic industrial applications, especially when the thermal conductivity of the solid phase is smaller than that of flow, considering conjugate heat transfer at the interphase surface would be more physically necessary.

Because of the coupling of thermal boundary conditions at both sides of the boundary surface in conjugate heat transfer, numerically modeling these coupled equations is the first step in the IB methods. Kang et al. (2009) transferred the CHT boundary condition to the Dirichlet boundary condition of temperature at the fluid side and Neumann boundary condition of heat flux at the solid side. The temperature and flux at each target boundary are interpolated from the opposite side, respectively. With a similar idea, Crocker et al. (2014) interpolated the heat flux through the one-side difference in each side of the interface under a cut-cell method. Under the framework of the ghost-cell method, Nagendra et al. (2014) and Das et al. (2016) applied the ghost cells at both fluid and solid regions to satisfy the thermal balance. Recently, a simplified numerical algorithm was proposed to solve the conjugate heat transfer problem. By assuming that the products of density and heat capacity ρCp are identical in both fluid and particle phases, the original governing equation of energy, i.e., ρCp(∂T∂t+u⋅∇T)=∇⋅(λ∇T), could be reduced to ∂T∂t+u⋅∇T=∇⋅(α∇T), where α=λ/(ρCp) is the thermal diffusivity. This simplified equation could be solved in the whole computational domain including the particle region. Consequently, the calculation of the heat flux at the interphase surface is the only remaining question. Ardekani et al. (2018a) used the volume of fluid method to calculate the velocity and thermal diffusivity at the interphase surface. In the work of Sato et al. (2016), the heat flux at the interfacial boundary was decomposed into the normal and tangential component, so that one more constraint of ∂T/∂τ‾|fluid=∂T/∂τ‾|solid was applied.

The first attempt to employ the IB method to deal with the CHT problem was made by Iaccarino and Moreau (2005). They combined the ghost-cell method with a RANS solver and local grid refinement technique to investigate the CHT problem of an automotive engine cooling fan system under natural and forced convective conditions. Their results agreed well with those obtained under the body-conformal grid system, which showed the capability of the IB method to study CHT problems. Kang et al. (2009) studied CHT problems of turbulent flows on a heated cylinder with the ghost-cell method. They found that the inclusion of solid conduction could lead to better results of local Nusselt numbers compared with experiments.

Takeuchi et al. (2013) extended the direct-forcing method of Kajishima et al. (2001) to include heat conduction to study the natural convection of a group of movable two-dimensional cylinders in a square enclosure with a cold bottom and a hot top considering the CHT at the cylinder surfaces. The Prandtl number, particle–fluid density ratio, and specific heat ratio are all set as unity. When the heat conductivity ratio is low, such as λ _s /λ _f  = 10⁻³ and 10⁰, the fluid and particles were found to circulate the domain center in one direction. However, at a high heat conductivity ratio, i.e., λ _s /λ _f  = 10¹ and 10², particles tended to oscillate around the domain center at a frequency relevant to the thermal diffusion time scale. They further investigated the effects of solid conductivity to heat transfer (Tsutsumi et al. 2014) and found that in the case of particles with lower solid conductivity, such as λ _s /λ _f  = 10⁻³, particles tended to accumulate near the wall region; while in the case of higher solid conductivity, such as λ _s /λ _f  = 10³, particles distributed more uniformly. Particles near the wall region would lower the heat transfer from hot wall to fluid, thus leading to a decreased Nusselt number. Later on, this method was extended to a three-dimensional case considering the effects of particle contacting heat transfer (Takeuchi et al. 2015). They found that when particles formed a chain contacting structure, a burst of heat flux occurred in the region across several particle diameters, as highlighted in Figure 16. This burst would cause a large temperature gradient near the particle chain and efficiently enhance fluid heat transfer and interphase heat exchange.

Figure 16:

Instantaneous flow field (left) and conductive heat flux (right) of a 2D natural convection problem with multiple particles in a square domain. The vectors represent fluid velocity and conductive heat flux, respectively, while in both figures the isocontours represent the temperature. Particles are represented by black circles. Reprinted from Takeuchi et al. (2015) with permission of Begell House.

Das et al. (2017) applied the ghost-cell IB method to simulate the heat transfer of a randomly packed structure of spherical particles in slender fixed-bed reactors. A CHT boundary condition was imposed on the surface of the spherical particles. The particle diameter was set as constant, while the diameter of the cylindrical domain was varied from 4 to 12 times of particle radius. Particle Reynolds number was varied from 1 to 500. It was found that increasing the solid conductivity from 1 to 100 would slightly enhance the overall heat transfer coefficient. This enhancement was further flattened at a high particle Reynolds number of Re = 500, which indicated that conjugate heat transfer had little effect on the overall wall-to-bed heat transfer. Further, this method was extended to simulate flow through a three-dimensional highly complex random solid foam structure under the CHT boundary condition (Das et al. 2018). The conductivity ratio was chosen as λ _s /λ _f  = 10, 100, and 1000. Although the flow field was little influenced by the conjugate heat transfer process, a higher temperature distortion could still be observed at λ _s /λ _f  = 1000, as shown in Figure 17.

Figure 17:

(A) The temperature at the solid surface of random foam structure at λ _s /λ _f  = 1000. (B) Temperature contours at mid-plane for λ _s /λ _f  = 10(left) and λ _s /λ _f  = 1000(right). The irregular hollow structures are crossed section of random foam. Reprinted from Das et al. (2018) with permission of Elsevier.

Ardekani et al. (2018a) combined the direct-forcing method and the volume of fluid method to study the conjugate heat transfer of the rigid spherical particle in a laminar Couette flow. The interphase heat transfer was studied by varying the particle Reynolds number from 1 to 16, the volume fraction of particles from 3 to 30%, and ratio of the thermal diffusivity of the particle to the fluid from 0.1 to 1. They found that when the particle volume fraction is lower than 10%, the increase of the particle Reynolds number could effectively enhance the interphase heat transfer. However, at higher volume fractions, the increase of Reynolds number cannot prompt the overall interphase heat transfer due to the larger enhancement of viscosity effect. It was also pointed out that the contribution of heat conduction in solid phase to the overall heat transfer was reduced as the particle Reynolds number increased. Using the same method, Ardekani et al. (2018b) studied the motion and heat transfer of spherical particles in a pipe flow at laminar and turbulent regions with Re = 370 and 5300, respectively. The solid volume fraction was varied from 0 to 40%. In the laminar regime, fully resolved particles greatly enhanced heat transfer (with a maximum of 330%) and wall drag compared with the case of no particle. However, in the turbulent regime, particles would relaminarize the turbulent pipe flow and severely damp the heat transfer. Only in the case of a volume fraction of 5%, the heat transfer coefficient was higher than that in the single-phase turbulent flow.

Apart from application in multiphase problems, CHT also gains popularity in characterizing the complex thermal boundary conditions of a wall or confined geometry. One noteworthy application is considering CHT in investigating flame–wall interaction (FWI) problems (Kedia et al. 2014; Michaels and Ghoniem 2016). Using the IB method to study CHT and FWI problems is quite straightforward. Consider the FWI problem, only the no-penetration boundary condition of species should be additionally satisfied at the immersed interface as ∂Y/∂n‾|Γ=0, and this Neumann boundary condition is easy to implement in the framework of the IB methods (Kedia et al. 2014).

To sum up, significant efforts have been made to apply the IB methods to simulations of the CHT problems. Numerical results demonstrate that including the conduction heat transfer within the solid could reproduce better collapse on experimental results. However, the multiphase heat transfer model considering the effect of CHT is rarely investigated. Due to the complexity in characterizing CHT problems, the influence of CHT on fluid convective heat transfer, interphase momentum, and energy transfer needs further investigations.

4.3.1 Permeation problems

The permeation problems are widely encountered in the studies of medical science, such as blood cells or medicine delivery. Gong et al. (2014) extended the idea of Huang et al. (2009) to the multiple moving interface problems governed by the convection–reaction–diffusion equation by using the original IB method. The effects of mechanical and geometry properties were analyzed. Figure 18 shows the oxygen transfer over the IB. Simulations indicated that the effect of the hematocrit value became a constant when the value was above 0.362. The stiffness influenced the oxygen transfer rate of the symmetrically distributed cells significantly, while it showed little effect on the asymmetric configuration. Wei et al. (2014) presented a direct-forcing IB method for simulations of sorption behaviors of porous media. A sink term was added in the species transport equation to represent the porous effect of the smaller-scale microporous material. It was found that the species could not be absorbed by a single layer of the porous structure, which could be used to determine the number of microscopic porous layers and provide valuable information for the design of the gas filter. Roghair et al. (2016) studied heat and mass transfer in a bubble swarm by combining the original IB method with the FT method under the framework of a finite-volume solver. The method was suitable to simulate the mass transfer of bubble swarms at a gas holdup range of 4–40%. Results showed that the mass transfer rate to the liquid phase just increased slightly with the increase of the gas holdup.

Figure 18:

Contours of the dimensionless concentration for oxygen transfer over the passively moving and deforming interface immersed in a steady flow field at different time instants. The yellow region represents the interior of the interface. Reprinted from Gong et al. (2014) with permission of Elsevier.

The key to the permeation problem is how to enforce the jump conditions of the concentration and pressure on the surface. Jayathilake et al. (2010) studied the effect of deformation and swelling of the membrane on the permeability properties by using the IIM to implement the jump conditions of solute concentration at the boundary surface. The jump conditions were calculated based on the Kedem–Katchalsky relations (Kedem and Katchalsky 1958). The accuracy of this method was validated with numerical errors within 1%. The simulations showed that the osmotic load influenced the capsule swelling and deformation observably. The capsule took a longer time to reach an equilibrium state when the osmotic load was smaller. Besides, the initial solute concentration ratio had a positive effect on the capsule swelling. This method was also applied by Jayathilake et al. (2011) to simulate the adhesive capsule, where the capsule–substrate adhesion was enforced using an adhesion potential (Seifert 1991). The solute absorbed by the vessel can be enhanced by increasing the diffusive permeability of the membrane. The capsule moving along the vessel increased the solute absorption by the vessel, and the adhesion strength influenced the solute transfer positively. Figure 19 shows the deformation and movement of the capsule. The solute concentration was much larger in a nonadhesive capsule, while the concentration inside an adhesive capsule was little different from that outside. The main reason was that adhesive capsules tended to move toward the wall, and more solute mass would be absorbed.

Figure 19:

The deformation and movement of nonadhesive ((A) and (B)) and adhesive capsules with different adhesive coefficients W _ad ((C) and (D)); the contours represent the solute concentration at the time of t = 88.36 s. The dashed/solid circles represent the capsules at different time instants. Reprinted from Jayathilake et al. (2011) with permission of Elsevier.

Musielak et al. (2009) studied the nutrient acquisition of diatom chains in different kinds of background flows and nutrient sources by using the original IB method. The uptake of nutrients was transmitted to the fluid domain by using a reaction term in the species concentration equation. The role of the mechanical properties of the chain between the diatoms was clarified, and some noteworthy points about the evolution of diatoms were pointed out. For the chain with four or five cells, the stiff chains took in more nutrients than solitary cells, and this effect was enhanced when the stiffness increased. Lambert et al. (2013) studied the active suspension of particles in a thin film by using the direct-forcing IB method. The nutrient uptake and swimming accumulation tendency of particles were investigated. A detailed analysis of the effect of the absorption rate, the swimming speed, and the volume fraction on nutrient uptake was performed. They clarified that although swimming took energy, it was still a winning strategy for microorganism to survive in the relatively dense suspensions.

4.3.2 Mass transfer with heterogeneous reactions

Researches involving heterogeneous reactions are very complicated since the heat and mass transfers are fully coupled, and the boundary conditions might be influenced by the reactions. The early studies of mass transfer with heterogeneous reactions were simplified by ignoring the thermal properties and complicated transport phenomena (e.g., Stefan flow) of the heterogeneous reactions. The coupling effects of these multiphysics were only investigated in the recent years.

In simplified studies, the Stefan flow phenomenon and energy equation, as well as the chemical reaction heat, are sometimes ignored in imposing boundary conditions. Yao and Fogelson (2012) formulated the direct-forcing IB method as an augmented system to simulate the chemical transport in the fluid field outside reactive particles. The forcing term and the value of species concentration were solved simultaneously in an augmented system of equations, which consisted of discrete chemical transport equations and discrete boundary conditions so that the boundary conditions were enforced explicitly. It was demonstrated to be more stable and accurate than the traditional direct-forcing IB method, which calculated the correcting forcing term using a prediction species concentration value of the next time step (the so-called predictor–corrector forcing method). The Stefan flow was neglected, and the energy equation was not included. Sulaiman et al. (2019) used the IB method coupling with the fictitious domain method to simulate the convective mass transfer around reactive particles in an incompressible flow. The ranges of Reynolds number and Schmidt number were from 10 to 100 and 0.5 to 10, respectively. The species diffusion and reactions inside the particle were considered. The chemical reaction was assumed to be first order and had no backward influence on the flow field. The convective heat transfer and interphase heat transfer, as well as the Stefan flow, were neglected. The effect of the reaction was treated as a conjugate mass transfer problem due to the continuity of species concentration and species normal flux. Figure 20 shows the isosurfaces of concentration within and outside three aligned reactive spheres, which are significantly dependent on the Damköhler number. Accordingly, a modified correlation of Sherwood number for a particle with an internal first-order reaction was proposed, which was based on an existing nonreactive Sherwood number correlation.

Figure 20:

Steady-state convection–diffusion process in the flow past three aligned reactive spheres at different Damköhler numbers ϕ²; isosurfaces represent the species concentration. Reprinted from Sulaiman et al. (2019) with permission of Elsevier.

In the studies including heat transfer effect, the coupling of energy equation and heterogeneous reaction at the immersed surface poses more difficulties in constructing the coupled boundary conditions. Deen and Kuipers (2013, 2014a, 2014b) extended the direct-forcing method to simulate coupled heat and mass transfer in reactive dense particles-laden flows. The Reynolds numbers of the exothermal chemical reactive particles were from 120 to 240 (Deen and Kuipers 2014b), and the mass transfer was simplified as a completely mass-transfer-limited condition. Several cases were presented under forced convection accompanied by a fast exothermic surface reaction. It was found that the mass transfer coefficients agreed well with the empirical equation (Deen and Kuipers 2013) and the correlation proposed by Gunn (Gunn 1978) for the cases of isolated particles and dense arrays, respectively. However, the mass transfer of a single particle in the dense array had a remarkable variation due to the preferred pathways of the fluid (Deen and Kuipers 2014b). Thereafter, a ghost-cell IB method was developed to simulate the mass transfer in the dense fluid–particle system (Lu et al. 2018a). The Reynolds numbers varied from 20 to 400, and the Damköhler numbers varied from 0.001 to 1000. The first-order surface reaction was considered instead of an infinitely fast reaction in the previous study. However, the Stefan flow was still neglected. As a result, the species boundary condition has the form of the Robin boundary condition. Figure 21 shows the concentration field of the dense particle system at different Damköhler numbers. With the increasing of Damköhler numbers, the reactant decayed quickly in the front part of the array, and full conversion of the reactant was achieved. The Sherwood number was influenced by the Damköhler number in the particle array. The temperature and species concentration profile of the dense particle array were compared with the results calculated according to the 1D heterogeneous model (Lu et al. 2018c), and good qualitative agreement was achieved. Nevertheless, the discrepancy caused by the heterogeneity of the array was amplified in the conditions of a high Reynolds number. The same method was also utilized to study the mass transfer in a catalyst cluster when diluted by inert particles (Lu et al. 2018b). The distribution of two kinds of particles led to a large variation of the mass transfer behavior. When the catalyst particle was surrounded by inert particles, the mass transfer became more active, especially at small cluster proximities.

Figure 21:

Concentration contours in the central slice of flow past the fixed particle array at three Damköhler numbers (0.1, 1, 10 from top to bottom). Reprinted from Lu et al. (2018a) with permission of Elsevier.

Nevertheless, the aforementioned studies do not consider the Stefan flow effect. Inclusion of Stefan flow and heat transfer caused by heterogeneous reactions would be much more physically realistic. McGurn et al. (2013) investigated the conjugate heat transfer and mass transfer processes associated with charring solids. The moving interface was described by a level-set method, and the ghost-cell method was used to enforce the boundary conditions. The effects of Stefan flow were superimposed through an explicit source term deposition into the Eulerian flow field. Lupo et al. (2019) presented a new direct-forcing IB method based on Breugem’s study (2012), which can fully capture the momentum, energy, and mass exchanges at the interface. Several validation cases of evaporating droplets were performed, including fixed rate and free evaporation of a static droplet, displacement of a droplet by Stefan flow, and evaporation of a hydrocarbon droplet in homogeneous isotropic turbulence. The Stefan flow velocity, temperature, and species distribution showed good agreements with analytical solutions and experimental data for both laminar and turbulent flows. Luo et al. (2016b) improved a ghost-cell IB method so that it could handle all the Dirichlet, Neumann, and Robin types of boundary conditions and further extended it to perform a fully resolved simulation of char particle combustion (Luo et al. 2018). In this study, the surface species transfer was a balance of diffusion, Stefan flow, and the heterogeneous reaction. The velocity of the interface was under a combined effect of the Stefan flow and the particle shrinking velocity. The heat transfer on the surface was coupled with a heterogeneous reaction. The char surface temperature was kept constant and varied from 1500 to 1900K in different cases, and the Reynolds numbers were in a range of 2.5–30.0. The method was successfully validated against the experiment (Makino et al. 2003), and the jump of the combustion rate due to the changing of the dominant reaction can be captured accurately. Figure 22 shows the flame shape around the char particle under different conditions. It can be seen that the CO flame was detached from the surface when the particle temperature was high enough. The gas-phase reaction in the boundary was always transport-limited, and the drag force coefficient and the Nusselt number became larger due to the chemical reactions. Furthermore, this method was extended by Zhang et al. (2020) to investigate the effects of heterogeneous reactions and homogeneous reactions on the drag force of a char particle. Simulation results indicated that both heterogeneous and homogeneous reactions tended to increase the pressure drop along the particle because of the accumulation of reactants at the rear region of the particle. Meanwhile, the high particle temperature led to high viscosity of gas phase in the particle boundary layer, which further increased the drag force.

Figure 22:

CO flame (shown as the contour of CO consumption rate) around the burning char particle (green cylinder in the figure) at different particle temperatures. Reprinted from Luo et al. (2018) with permission of John Wiley and Sons.

In summary, the study of multiphase flow with mass transfer and the heterogeneous reaction has been carried out by many researchers, and a lot of achievements have been obtained. However, how to model the boundary condition physically as much as possible remains an open question. So far, there is no boundary model that considers the effect of finite reaction rate, the Stefan flow, and heat and mass transfer within in the particle at the same time. Although sophisticated, developing a realistic boundary model for the heterogeneous surface reaction is in need. Moreover, improved correlations for multiphase transport phenomena can be obtained with the help of the IB methods, which will further improve the predictive capability of traditional coarse models for industrial applications.