An efficient multi-grid finite element fictitious boundary method for particulate flows with thermal convection

International Journal of Heat and Mass Transfer 126 (2018) 452–465

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

An efficient multi-grid finite element fictitious boundary method for particulate flows with thermal convection Khuram Walayat a,b, Zekun Wang a,b, Kamran Usman c, Moubin Liu a,b,⇑ a

BIC-ESAT, College of Engineering, Peking University, Beijing 100871, China StateKey Laboratory for Turbulence and Complex Systems, Department of Mechanics and Engineering Science, Peking University, Beijing 100871, China c Department of Mathematics, Faculty of Basic and Applied Sciences, Air University, Islamabad 44000, Pakistan b

a r t i c l e

i n f o

Article history: Received 16 October 2017 Received in revised form 30 March 2018 Accepted 1 May 2018

Keywords: Direct numerical simulation Fictitious boundary method Multigrid Particulate flows Thermal convection Sedimentation

a b s t r a c t This paper presents a direct numerical simulation (DNS) technique, the Finite Element Fictitious Boundary Method (FEM-FBM) for the simulation of fluid-solid two-phase flows with heat transfer. The heat transfer equation is introduced to study thermal convection in fluid-solid two-phase flows. The Boussinesq approximation is considered for the coupling of momentum and temperature flow fields. Multi-grid finite element solver is used to compute flow equations of mass, momentum, and energy on a fixed Eulerian mesh which is independent of time and the solid particles are allowed to move freely in the whole computational domain. Fictitious boundary method (FBM) is used to treat the particles inside the fluid, FBM takes account of the thermal and momentum interaction between the fluid and the particles. The accuracy and stability of presented method are validated by comparing our test cases results with the results reported in the available literature. Numerical tests are performed to show that this method is potentially powerful and provides an efficient approach to simulate complex thermal convective particulate flows with a large number of particles. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction Particulate flows or motion of solid particles in fluids have a wide range of industrial applications, such as fluidized suspensions, lubricated transport, sedimentation, hydraulic fracturing of reservoirs, slurry flow, paper pulp, food products etc. These types of flows are common in many natural processes such as sand or dust particles in air blown by the wind, ocean current interaction with offshore structures, lava flow and sedimentation in estuary etc. Particulate flows in biological processes have been a subject of great importance with research contributions coming from the field of biology, chemistry, physics, engineering, and mathematics. The sedimentation of suspended particles has a great importance in the chemical, petroleum, paper pulp, wastewater, food, pharmaceutical, ceramic and other industries as a way of separating particles from the fluid as well as separating particles of different types settling with different speeds. Here we are particularly interested in the gravitational settling of particles with convection heat transfer.

⇑ Corresponding author at: BIC-ESAT, College of Engineering, Peking University, Beijing 100871, China. E-mail address: [email protected] (M. Liu). https://doi.org/10.1016/j.ijheatmasstransfer.2018.05.007 0017-9310/Ó 2018 Elsevier Ltd. All rights reserved.

Particulate flows are quite complex and hard to simulate numerically because frequent generation and deformation of the computational grid are required in many cases when the particle boundaries are complex and moving with time. The problem becomes more complex in the case with a large number of particles due to fluid-particle interaction as well as due to particle-particle and particle-wall collisions. Rapid advancements in computational power make the direct numerical simulation (DNS), an important and practical tool to study particulate flow mechanism. It treats the fluid and solid objects separately. The DNS approach is based on Navier–Stokes equation for the fluid and equations of rigid body motion for particles. A variety of DNS numerical schemes have been proposed over the past decade to simulate fluid-particle flow problems. These methods are broadly classified into two types, one is based on Lagrangian approach while other is a Eulerian approach. In Lagrangian approach, the mesh moves and follows the moving boundaries of the particles in the fluid. Since the motion of the mesh can be defined arbitrarily within the fluid, therefore this approach is usually called Arbitrary Lagrangian Eulerian (ALE). Hu et al. [14,15], Maury [12,13] and Feng et al. [16] have used the ALE method to study particulate flows. ALE method normally requires generating a new mesh at every time step in the case of moving particles, so it is computationally expensive especially for the simulations of problems with large number of

453

K. Walayat et al. / International Journal of Heat and Mass Transfer 126 (2018) 452–465

particles. Whereas the Eulerian approach is more efficient than the Lagrangian one but the resulting accuracy is often not as clear as in case of Lagrangian. Therefore, the overall objective is to deal with the moving boundaries in the fluid, improve the accuracy of the numerical approximation and reduce the computational cost. Eulerian methods do not require re-meshing, a fixed Cartesian mesh is required which covers the whole computational domain comprising of both particles and fluid. Peskin [17] introduced immersed boundary method (IBM), based on a Eulerian approach to study fluid-solid interaction problems. Similar to IBM, Glowinski et al. [7,18,19] developed a finite element fictitious domain method to simulate fluid-particle flow problems. Turek et al. [1,2,3,10] presented a multigrid finite element fictitious boundary method (FBM) for the simulation of particulate flows. In recent past, many hybrid methods have been developed to simulate sedimentation of suspended particles in fluid, For example, lattice Boltzmann method (LBM) is coupled with other numerical methods, such as IBM [27], Finite Element Method (FEM) [28] and Discrete Element Method (DEM) [29]. Sun et al. [31] and Hager et al. [30] used another hybrid method CFD-DEM, in which computational fluid dynamics (CFD) is coupled with DEM for the study of sediment transport. Although, heat transfer in particulate flows is involved in many industrial applications, yet only a few results are found in the literature. Among the reported studies, Kim et al. [35] proposed an Immersed Boundary Method (IBM) in context of a Finite Volume (FV) fixed grid scheme for the solution of heat transfer in complex two-dimensional geometries without considering the solid-phase motion. Pacheco et al. [25] also presented FV-IB method with a non-staggered grid technique for heat transfer and fluid flow around fixed particles. Demirdzˇic´ et al. [24] used finite volume nonorthogonal boundary-fitted grids to study the thermal convection around fixed solid-phase. Gan et al. [20,21] used ALE scheme for the simulations of two-dimensional particle sedimentation problem with heat transfer between particles and the surrounding fluid. Yu et al. [22] employed fictitious domain method to solve particulate flow problem and extended it to study heat convection at the interface between the fluid and solid particles. Feng et al. [11,34,37] developed finite volume IBM method for heat transfer between the fluid and moving solid particles in particle-laden flows. The numerical study of thermal convection in particulate flows has been attempted by some other researchers as well [36,38,39]. In the present work, we apply the multi-grid FEM fictitious boundary method [1,2,3,10,33] to simulate particle sedimentation problems and extend this method to study heat transfer in solidliquid two-phase flows. The considerable advantage of multi-grid FEM fictitious boundary method which makes it an efficient scheme is that it is based on a fixed FEM background grid which is independent of flow features, hence re-meshing is not required with the edge of high convergence rate of the multi-grid solver. By applying boundary conditions at the interface between fluid and particles which become an additional constraint to the governing equations, so the fluid domain can be extended into the whole domain which covers both fluid and particles. This paper is organized as: in the second section, we present the mathematical modeling and governing equations of the problem and discussed fluid-particle interaction. The computational scheme is described in the third section. The fourth section consists of numerical experiments to check the validation of our method and test its efficiency and potential to simulate real particulate flows. The conclusion of the presented study is given in the final section. 2. Mathematical modeling We consider the flow of N number of solid particles of mass M i ði ¼ 1; 2; . . . ; NÞ in an incompressible Newtonian fluid, as shown

Fig. 2.1. Rigid particle and fluid.

in Fig. 2.1. The density of the fluid is qf and its viscosity is m. Xf ðtÞ and Xi ðtÞ denotes the domain occupied by the fluid and the ith particle at time t respectively. where XT is the total domain and is given as,

XT ¼ Xf ðtÞ [ Xi ðtÞ 8 i 2 ð1; 2; . . . ; NÞ

ð2:1Þ

XT as an entire computational domain is independent of t. As Xf and Xi are always depended on time t we denote Xf ðtÞ ¼ Xf and Xi ðtÞ ¼ Xi dropping t in the notations. Where @ Xi represents the boundary of the ith particle. 2.1. Fluid flow model The motion of an incompressible fluid with density qf is governed by the equations of continuity, momentum, and energy in the domain Xf as,

r  u ¼ 0 ðaÞ

 @u þ u  ru  r  r ¼ f ðbÞ @t
 @T qf c f þ u  rT ¼ kf r2 T 8 X 2 XT @t

qf




ð2:2Þ ðcÞ

r is the total stress tensor in the fluid phase defined as,

r ¼ pI þ lf ½ru þ ðruÞT :

ð2:3Þ

where u is the fluid velocity, p is the pressure, lf is the coefficient of viscosity, f is the source term and I is the identity tensor. T is the temperature, cf is specific heat and kf is the thermal conductivity of the fluid. 2.2. Particle motion model The rigid particles are free and allowed to move in the fluid domain. The particles have both translational and rotational motion under the action of gravity, forces due to fluid called hydrodynamic forces and collision forces due to interactions between particle-particle or particle-wall. The motion of solid particles is governed by the Newton-Euler equations [2,9], i.e. if U i and xi are the translational and angular velocities of the ith particle respectively, then the particle must satisfy the following equations,

mi

dU i ¼ ðDmi Þg þ F i þ F 0i ; dt

Ii

dxi þ xi  ðIi xi Þ ¼ si ; dt

ð2:4Þ

where Mi is the mass of the ith particle and if Mf is the mass of fluid occupying the same volume as Mi then DMi is given by the mass difference between the mass of the ith particle Mi and the mass of the fluid Mf ,

DM i ¼ M i  M f ;

ð2:5Þ

F i represents resultant hydrodynamic i.e. drag and lift forces acting on the ith particle, F 0i are the collision forces acting on the ith particle due to particle-particle and particle-wall collision, I i

454

K. Walayat et al. / International Journal of Heat and Mass Transfer 126 (2018) 452–465

is the moment of inertia tensor of the ith particle, si is the resultant torque acting about the center of mass of the ith particle and g is the gravitational acceleration. 2.3. Fluid-particle interactions In particulate flows, the particles and fluid interact with each other and exert forces on one another. In this interaction, the exchange of momentum and energy takes place between fluid and particles. 2.3.1. Hydrodynamic forces and torque The hydrodynamic drag and lift forces F i and the torque si about the center of mass of the ith particle can be obtained by [23],

Z

F i ¼ ð1Þ

@ Xi

Z ¼ ð1Þ

@ Xi

ðr  nÞdCi ;

si

ðX  X i Þ  ðr  nÞdCi ð6Þ

ð2:6Þ

where r is the total stress tensor in the fluid phase defined by the Eq. (2.3), X i is the position of the center of mass of the ith particle, @ Xi is the boundary of the ith particle and n is the outward drawn unit normal vector on the boundary @ Xi of the ith particle. Let S be the surface of rigid particles ns be the inward drawn unit vector nor! mal to the surface S of the particles and T is the tangential vector given by,

! T ¼ ðny ; nx Þ

ð2:7Þ

Then using Eqs. (2.6) and (2.3) the drag and lift forces are calculated as,

Z 

 @us ny  pnx ds; @ns s  Z  @u ¼ l s nx þ pny ds; @ns s

l

FD ¼

FL ð2:8Þ

The drag and lift coefficients are given by,

Cd ¼

2F D

qU 2 D

;

Cl ¼

2F L

qU 2 D

;

ð2:9Þ

where C d and C l are the coefficients of the drag and lift forces respectively, U is the characteristic velocity and D is the characteristic length. From Eqs. (2.8) and (2.9) it is clear that surface integral in Eq. (2.8) should be conducted for the calculation of C d and C l . 2.3.2. Momentum interaction Let X i be the position of the center of mass of the ith particle and hi be its angle, then the position X i and angle hi can be obtained by integrating the following kinematic equations [2,10],

dX i ¼ Ui; dt

dhi ¼ xi : dt

ð2:10Þ

By the application of no-slip boundary conditions at the interface @ Xi between the fluid and the ith particle, the velocity UðXÞ 8 X 2 Xi is given by,

UðXÞ ¼ U i þ xi  ðX  X i Þ;

ð2:11Þ

where U i is the translational velocity of the ith particle as shown in Fig. 2.2 and the second term of the sum in Eq. (2.11) is the tangential part of the angular velocity of the ith particle. 2.3.3. Energy interaction Similar to momentum interaction, the energy equation Eq. (2.2 (c)) describes the temperature field of the whole computational region [11]. We consider a uniform temperature on particles sur-

Fig. 2.2. fluid particle interface.

face and near to the particles surface the fluid and particles temperatures are equal, i.e. T f ¼ T p .

qp V p cp

dT p ¼ dt

I

@X

kf ðrT f  nÞds þ

Z

X

ð2:12Þ

qp dV;

this temperature equation is obtained by using energy balance for each particle. Where qp is energy sources within the particles, cp ; V p and qp are the specific heat, volume, and density of the particles respectively. 2.3.4. Boussinesq approximation To study the momentum and heat interactions between the fluid and particles, we employ the Boussinesq approximation. The combined set of momentum equations with heat transfer in terms of Boussinesq approximation are written as follows [1]

8 ru¼0 > < qf @u þ qf u  ru ¼ f  rp þ mr2 u þ qf jT @t > : @T þ u  rT ¼ mT r2 T @t

ðaÞ 8 X 2 XT ; ðbÞ 8 X 2 Xf ; ðcÞ 8 X 2 XT ; ð2:13Þ

where mT is the diffusion coefficient for temperature T and j is the unit vector along the y-direction. The Boussinesq approximation for the buoyancy forces is often used for coupling of momentum and energy equations and is valid for relatively small temperature differences between the fluid and particles. 2.4. Collision model To handle particle-particle and particle-wall interaction we employed a collision model presented by Singh et al. [32]. In this model, a new method of short-range forces of repulsion has been introduced which not only stop the particles from getting very close, but it can also deal with the overlapping of particles when numerical simulations bring these particles very close to each other due to inevitable numerical errors. 2.4.1. Particle-particle interaction For the particle-particle collisions, the force of repulsion can be obtained as,

F Pi;j

¼

8 > > < > > :

0;

for Di;j > Ri þ Rj þ q;

2 ep ðX i  X j ÞðRi þ Rj þ q  Di;j Þ ; for Ri þ Rj 6 Di;j 6 Ri þ Rj þ q; 1

e0p ðX i  X j ÞðRi þ Rj  Di;j Þ; 1

for Di;j 6 Ri þ Rj ; ð2:14Þ

455

K. Walayat et al. / International Journal of Heat and Mass Transfer 126 (2018) 452–465

where X i and X j are the coordinates of the centers of the ith and jth particle respectively and Di;j ¼ jX i  X j j is the distance between the centers of mass of the particles. Ri and Rj are the radius of the ith and jth particle respectively. q is called the range of the force of repulsion and usually q ¼ 0:5  2:5Dh; where Dh is the size of the mesh element. ep and e0p are small positive stiffness parameters for particle-particle collisions. The selection of the values of stiffness parameters ep and e0p are such that they do not cause a discontinuity or singularity. For a sufficiently viscous fluid, and q ’ Dh as well as qqi are of order 1 (qi is the density of the ith particle and qf is f

the fluid density), then the values of ep ’ ðDhÞ and e0p ’ Dh are used in the calculations 2

2.4.2. Particle-wall interaction For the particle-wall collisions, the corresponding force of repulsion is expressed by,

FW i

¼

8 0; > > <

ew ðX i  1

> > :

e0w ðX i 1



X 0i Þð2Ri X 0i Þð2Ri

þq

2 D0i Þ ;

 D0i Þ;

for

D0i > 2Ri þ q;

for

2Ri 6 D0i 6 2Ri þ q;

for

D0i 6 2Ri ;

 Multigrid starts with an initial guess on fine grid level i i.e., u0i and executes pre-smoothing to obtain more accuracy

ujþ1 ¼ Si ðuij Þ: i

ð2:15Þ X 0i

where is the coordinate of the center of mass of the nearest imaginary particle P0i imagined on the boundary wall C with respect to the ith particle, D0i ¼ jX i  X 0i j is the distance between the center of the imaginary particle P 0i and the mass center of the particle. ew and e0w are small positive stiffness parameters for particle-wall cole

Hackbush [4]. We can use multi-grid approach to solve the fluid particle interaction problems which are based on a number of grids, obtained by regularly refining the coarse mesh. For CFD problems, multi-grid is one of the fastest linear solver [5]. In multi-grid, the restriction is applied to the residual after smoothing on all mesh levels, and direct sparse linear solver [6] is used to get the solution on the coarsest grid if the number of degrees of freedom is sufficiently small. Then prolongation is applied followed by post-smoothing to obtain a better approximation. These steps continue until the iterations of the multi-grid cycle (V or F-cycle) is finished. With the help of some operators, here we explained how multi-grid works on a problem using different grid levels. The multi-grid procedure to solve the linear system Ai ui ¼ bi is presented in the following steps.

e0

lisions, usually, their values can be taken as ew ¼ 2p and e0w ¼ 2p in the calculations. Fig. 2.3 shows all the forces acting on a particle. 3. Numerical implementation A brief description of the multi-grid solver is presented in this section and the FBM strategy to deal with the coupled temperature and velocity fields for both fluid and particles is also discussed here. 3.1. Multi-grid solver Multi-grid methods were invented for PDEs like Poisson’s equation, however, these methods also work on a large number of problems. Contrary to other iterative methods, the convergence rate of the multi-grid method is independent of problem size. For an introduction to multi-grid, the readers are referred to the book of

j ¼ 0; . . . ; m  1

ð3:1Þ

where Si is the smoothing operator and it computes the first improved approximation to the system Ai ui ¼ bi .  The high frequency of the residual is sufficiently smoothed by pre-smoothing so that the remaining error shows high frequency on a coarser grid m ri1 ¼ Ri1 i ðbi  Ai ui Þ;

ð3:2Þ

Ri1 i

where is the restriction operator from a finer grid to coarser grid which gives coarser grid approximation.  Solve the system

Ai1 ui1 ¼ r i1

ð3:3Þ

ui1 .

on the coarser grid to get the correction  Prolongate the correction ui1 to the finer grid level and apply

i  umþ1 ¼ um i þ ai P i1 ui1 ; i

Pii1

where is the prolongation operator from the coarser grid to next finer grid and ai is the damping parameter.  Execute the post-smoothing steps to get the final solution umþ1þn . i By applying these steps recursively on different grid levels, faster reduction of error is achieved. Appropriate algorithms for restriction, prolongation, smoother and solver components have to be chosen to get accuracy and efficiency. Fig. 3.1 shows a schematic diagram of multi-grid V-cycle (MGV), multi-grid W-cycle (MGW) and multi-grid F-cycle (MGF). The basic idea for the construction of multigrid V-cycle (MGV) or multigrid F-cycle (MGF) algorithm is as follows. Algorithm 1. Multigrid V-Cycle i

function MGV ðb ; ui Þ replace the approximate solution ui of the system Ai ui ¼ bi with an improved solution if i ¼ 1 compute the exact solution u1 return u1 else i

 ui ¼ Sðb ; ui Þ; i

i

 r ¼ A :ui  b ; i

i

 d ¼ PðMGVðRðri Þ; 0ÞÞ; i

 ui ¼ xi  d ; i

 ui ¼ Sðb ; ui Þ; return ui end if Fig. 2.3. Forces acting on a particle.

ð3:4Þ

456

K. Walayat et al. / International Journal of Heat and Mass Transfer 126 (2018) 452–465

Fig. 3.1. MGV, MGW, and MGF cycles.

3.3. Calculation of hydrodynamic forces and torque Algorithm 2. Full Multigrid (FMG) In particulate flows, the calculation of the hydrodynamic forces and torque acting on the particles moving in the fluid is very important for the study of interactions between the fluid and the

i

function MGF ðb ; ui Þ return an accurate solution ui of the system Ai ui ¼ bi solve A1 u1 ¼ b1 exactly to get u1 for i ¼ 2 to k i

ui ¼ MGVðb ; Pðui ÞÞ end for

3.2. Multigrid FEM fictitious boundary method Several approaches have been presented using fictitious boundary method to deal with the particles in the fluid and to calculate the hydrodynamic forces acting on the particles. Glowinski et al. [7] described a semi-implicit approach to calculate the drag and the lift forces acting on the particles in the fluid and their movement in the fluid. Patankar et al. [8] also used an implicit scheme for the particle treatment. Wan et al. [3] proposed an explicit way to treat the fluid, moving particles inside the fluid and the explicit calculation of the drag and the lift forces acting on the particles. Here we present a brief description of multi-grid FEM fictitious boundary method, for more details reader is referred to the articles [1,2,3]. This multi-grid fictitious boundary method is based on FEM background grid which covers the whole computational domain XT . The grid can be chosen independent of solid particles arbitrary shape, size, and number. The multi-grid FBM starts with a coarse mesh which may already describe the geometrical details of particles Xi ði ¼ 1; 2; . . . ; NÞ and a boundary parametrization with respect to the boundary conditions of Eq. (2.11) which sufficiently depict all fine-scale structures of the particles. The internal solid objects are introduced in the corresponding components in all matrices and vectors in the solution process as unknown degrees of freedom. Then, the extra conditions arising due to the interior objects are incorporated implicitly in all the iterative solution steps. Hence, by using Eq. (2.11), the computations can be carried out on the whole domain XT for the fluid. The FBM has a considerable advantage that the computational domain does not require to be changed in time, and no re-meshing is required. More precisely, the mesh and the flow features can be handled independently of each other [1,3]. Hence using the FBM, the domain of definition of the fluid velocity u is extended according to Eq. (2.11) which can be seen as an additional constraint to the Navier-Stokes equations. So, the coupled system of Eqs. (2.13) can be written as

8 ru¼0 > > > > < q @u þ q u  ru ¼ f  rp þ mr2 u þ q jT f @t f f > @T þ u  rT ¼ mT r2 T > > @t > : uðXÞ ¼ U i þ xi  ðX  X i Þ

particles. Let XT ¼ Xf [ fXi gNi¼1 be the total computational domain including domain occupied by the particles and the fluid. Let n be the unit vector drawn normal to the boundary @ Xi of ith particle pointed outward to the flow region. To calculate the hydrodynamic forces F i acting on the surface of the ith particle and the torque si acting about the mass center of the ith particle, the surface integrals on the wall surface of the particle given by Eqs. (2.6) should be conducted. Wan et al. [3] proposed a volume integral approach rather than the surface integral approach for the calculation of hydrodynamic forces and torque acting on the solid bodies moving in the fluid. In [3], they replaced the surface integral in Eqs. (2.6) by a volume integral, which is computationally less expensive, by defining an auxiliary function ai ,

ai ¼



1 for

X 2 Xi ;

0

X 2 Xf ;

for

:

The gradient of ai is zero everywhere except on the wall surface of the ith particle, and at the wall surface, the gradient equals to the unit vector n normal to the wall surface of the ith particle as shown in Fig. 3.2.,

ðaÞ 8 X 2 XT ; ðbÞ 8 X 2 Xf ; ðcÞ 8 X 2 XT ðdÞ 8 X 2 Xi ;ði ¼ 1;2;...;NÞ ð3:5Þ

ð3:6Þ

Fig. 3.2. Cells where ai ¼ 1 and n ¼ rai on the particle’s boundary.

K. Walayat et al. / International Journal of Heat and Mass Transfer 126 (2018) 452–465

n ¼ rai :

ð3:7Þ

Then for the calculation of hydrodynamic forces and torque, the surface integrals in Eqs. (2.6) can be replaced by volume integral as

Z F i ¼ ð1Þ

XT

Z ¼ ð1Þ

XT

ðr  rai ÞdX;

Ti

ðX  X i Þ  ðr  rai ÞdX:

ð3:8Þ

The volume integral over each element covers the total domain

XT can be truly calculated with the help of a standard Gaussian

quadrature of sufficiently higher order. As the gradient of ai i.e.

457

4. Simulation results and discussions For the validation of our numerical method and to check its efficiency and robustness, we performed number of numerical tests. In the first case, we simulate natural thermal convection due to a fixed hot circular cylinder placed in a square cavity. Second, we perform numerical experiments for a single particle with fixed temperature on its boundary, settling in an infinite channel and then we investigate the motion of a catalyst particle with freely varying temperature on its boundary. Next, we simulate the sedimentation of a circular cluster of 107 hot and cold particles and finally, the case for the sedimentation of 10,000 small hot particles.

rai is non-zero on the wall surface of the ith particle, therefore the volume integrals in Eq. (3.8) required to be evaluated only in one layer of mesh cells around the ith particle, which is very effective treatment. 3.4. Calculation of drag and lift coefficients By using Eqs. (2.3) and (3.7), the surface integral in Eq. (2.8) is replaced by volume integral for the calculation of drag and lift forces acting on the moving particles in the fluid. Such that,

XT


   @u @ a @u @ a @a l þ p dX; @x @x @y @y @x

ð3:9Þ

XT


   @v @a @v @a @a p dX: l þ @x @x @y @y @y

ð3:10Þ

Z FD ¼  Z FL ¼ 

Therefore, through Eqs. (3.9), (3.10) and (2.9) the drag and lift coefficients C d and C l can be calculated by conducting the volume integral over the total domain XT rather than the surface integral conducted over the wall surface of the rigid bodies in Eq. (2.8). 3.5. Energy balance equation Similar to the equations of hydrodynamic forces and torque, by using Eq. (3.7) the surface integral in energy balance equation Eq. (2.12) can be replaced by volume integral as

qp V p cp

dT p ¼ dt

Z XT

kf ðrT f  rai ÞdX þ

Z X

qp dV

ð3:11Þ

4.1. A particle placed eccentrically in a square cavity (benchmark) The test case of natural thermal convection around a fixed cylinder has been studied by many researchers both experimentally and numerically [11,22,24,25,26]. This problem is used for the validation of current method against the results of different methods used for the same problem. In this numerical experiment, the fluid motion is induced by heat transfer. The fluid motion due to heat transfer has a significant impact on the motion of particles, and therefore these phenomena have great importance in fluidsolid heat transfer problems. The computational domain is a 2D square cavity of length 1 and it is filled with Newtonian fluid. Fig. 4.1 shows the schematic illustration of the problem, where a particle of diameter d ¼ 0:4 is represented by a circular cylinder and is placed eccentrically in the square enclosure with eccentric distance e ¼ 0:1 above the center of the enclosure. The particle is fixed and the dimensionless temperature T ¼ 1 on the surface of the particle and 0 for the fluid. A constant 0 temperature is applied on both side walls, whereas adiabatic boundary conditions are applied at the top and bottom walls. The values of governing dimensionless flow numbers such as Reynolds number, Prandtl number, and Grashoff number are taken as Re ¼ 0:1; Pr ¼ 10 and Gr ¼ 105 respectively. We used a coarse mesh with number of elements Nel ¼ 625 at mesh refinement lev el  1, this coarse mesh is refined by connecting midpoints of opposite edges of each element, so that, an element divides into four elements after every refinement. The simulations are performed using mesh refinement lev el  2; lev el  3 and lev el  4. The number of elements Nel ¼ 25; 00 and element size Dx ¼ 0:02

3.6. Governing dimensionless physical parameters The dimensionless numbers and physical parameters involved and used in our work to simulate particulate flows with heat transfer are defined as

Reynolds number Re ¼

Prandtl number Pr ¼

lf cf

Grashoff number Gr ¼

qr ¼

qf U ref Lref lf

kf

qf gbðT s  T 0 ÞL3ref lf

qp cp ; cr ¼ qf cf

where qr and cr are the density and heat capacity ratios.

ð3:12Þ

ð3:13Þ

ð3:14Þ

ð3:15Þ Fig. 4.1. Schematic diagram.

458

K. Walayat et al. / International Journal of Heat and Mass Transfer 126 (2018) 452–465

at lev el  2; Nel ¼ 10; 000 and Dx ¼ 0:01 when lev el  3 refined mesh is used, where Nel ¼ 40; 000 and Dx ¼ 0:005 at lev el  4: The thermal contour plots are shown in Fig. 4.2 which shows, that the results of lev el  3 and lev el  4 are very much consistent. Therefore, we use meshes having an element size Dx 6 0:01 in the simulation cases of the other sections. We computed the local Nusselt numbers along one of the side walls and plot the results together with the results obtained by Demirdzic et al. [24] and Pacheco et al. [25]. Fig. 4.3 shows that the simulation results obtained using present method agrees very well with the results presented by Demirdzic et al. [24] and Pacheco et al. [25]. This comparison provided us a confirmation that the energy balance equation is implemented accurately in the present method. Isothermal patterns are shown in Fig. 4.4. Fig. 4.5 shows the streamlines patterns and fluid velocity vector plots. It is seen that the temperature near the top wall is much higher than the temperature close to the bottom wall. Streamlines pattern and velocity vector plots show that the fluid flow is dominated by two circular regions on both sides of the cavity. 4.2. Settling of a particle in an infinite channel (with fixed temperature on its boundary) Fig. 4.2. Thermal contour plots using refinement lev el  2; lev el  3 and lev el  4.

Fig. 4.3. Local Nusselt numbers along a side wall.

The sedimentation of a particle in an infinite vertical channel with constant temperature was first simulated and analyzed in detail by Gan et al. [20] using ALE finite element method. Later, Yu et al. [22] and Feng et al. [11] reexamined the case by using Eulerian grids. They found that the particle settles along the center-line of the channel at low Grashoff numbers and it began to oscillate about the center-line of the channel while settling as the Grashoff number increases. The computational domain is a rectangular channel of width 4 and height 8. The domain is shifting while the particle is settling, which makes this problem as a flow simulation in an infinite channel. The radius of the particle is 0:5 and its initial position is ð1:5; 4Þ. The dimensionless temperature on the particle surface is unity and initially 0 for the whole fluid domain. The values of dimensionless parameters, in this case, are same as used in [11,22], i.e. Re ¼ 40:5; Pr ¼ 0:7; qr ¼ 1:00232 and cr ¼ 1: The diameter of the particle is considered as reference length and pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U ref ¼ paðqr  1Þg is the reference velocity, where a is the radius of the particle. We used a grid with number of elements Nel ¼ 1152

Fig. 4.4. Isotherms, Left. Feng et al. [11] Right: Present results.

459

K. Walayat et al. / International Journal of Heat and Mass Transfer 126 (2018) 452–465

Fig. 4.5. Streamlines and velocity vector plots.

Table 1 The Lateral amplitude of oscillations at different Grashoff numbers. Maximum lateral displacement of the oscillating particle from the equilibrium position Grashoff number

0 100 564 1000 2000 2500 4500

Fig. 4.6. Time history of the lateral position of the settling particle.

at grid refinement lev el  1; Nel ¼ 73; 728 at lev el  4 and Nel ¼ 294; 912 at lev el  5: Fig. 4.6 shows the lateral equilibrium positions of the particle according to the Grashoff number regimes established by Gan et al. [20]. The plots show the positions of the particle about the center-line of the channel at Grashoff number Gr ¼ 0; 100; 564; 1000; 2000; 2500 and 4500: We can see that at Gr ¼ 0 and Gr ¼ 100 the particle settles along the center-line of the channel and wake behind the particle is steady and symmetric. At Gr ¼ 564 the particle experiences very weak oscillations about the center-line and small Vortices begin to shed from the particle. The center-line remain the equilibrium position of the particle but the regular oscillations are of different amplitudes at Gr ¼ 1000; 2000; 2500 and 4500: In the simulations of Gan et al. [20], the particle migrates away from the center-line and settles steadily close to one of the walls at Gr ¼ 1000 and 2000: Yu et al. [22] found that this migration of the particle away from the center-line is due to random numerical disturbances. However, our simulation results show that the center-line remains the equilibrium position for an oscillating particle in all the cases, which demonstrate that our method is more stable and accurate. Fig. 4.6 shows the lateral X-position of a set-

Amplitude of oscillations from center-line of channel Present results

Yu et al. [22]

Feng et al. [11]

0 0 0.0013 0.3785 0.7807 0.9110 1.3509

0 0 0.0015 – 0.75 0.89 1.33

0 0 0.0015 – 0.73 0.90 1.35

tling particle with respect to time. Table 1 shows the maximum lateral displacement of an oscillating particle from the centerline of a channel at different Grashoff numbers. From the data, it is observed that the particle begins to oscillate slightly at Gr ¼ 564 and the amplitude of these oscillations increases with increase in the value of Grashoff number. We simulate cases in a channel whose width equals to sixteen times the diameter of particle where particle fluid density ratio varies as 1:001 6 qr 6 1:02 and plot the drag coefficient C d of the particle as a function of Re at terminal velocity V T when the Gr ¼ 0 and 100. Fig. 4.7 shows the comparison of our results with the results obtained by Feng et al. [34]. The comparison shows good agreement of the results. It is observed that more hydrodynamic drag force acts on a hot particle ðGr ¼ 100Þ, at low Reynolds number than the isothermal particle. However, this difference in drag coefficients significantly decreases as the Re increases and it becomes almost independent of the thermal convection at higher Re. 4.3. Motion of a particle in a vertical channel (with freely varying temperature on its boundary) In this case, we consider a particle settling under the action of gravity, with a constant heat source inside the particle. Heat is generated homogeneously inside the particle due to the internal heat source which increases particle’s temperature. The temperature of the surrounding fluid also began to increase eventually due to heat transfer between particle and fluid. This heat transfer effect also

460

K. Walayat et al. / International Journal of Heat and Mass Transfer 126 (2018) 452–465

Fig. 4.8. Time history of vertical component of the velocity of catalyst particle.

Fig. 4.7. Drag coefficient against the particle’s Reynolds number.

has a strong influence on the motion of the particle along with gravitational effects. This problem represents the motion of a catalyst particle when the temperature of the particle is allowed to change with time and space. We used a computational domain of width 4 and height 8, which represents a rectangular channel. The particle of radius 0:5 is initially placed at the center ð2; 4Þ of the channel. Initially, velocity and temperature of both fluid and particle are zero. The uniform heat source qp ¼ 1 is over the entire surface of the particle. We set the governing dimensionless physical parameters for this problem as, Re ¼ 40; Pr ¼ 0:7 and Gr ¼ 1000. The reference velocity is pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U ref ¼ paðqr  1Þg , where a is the radius of the particle. The grid used has number of elements Nel ¼ 1152 at grid refinement where Nel ¼ 4608; Nel ¼ 18; 432; Nel ¼ 73; 728; lev el  1; Nel ¼ 294; 912 and Nel ¼ 1; 179; 648 at refinement lev el  2; lev el  3; lev el  4; lev el  5 and lev el  6 respectively. We simulate the following three cases Fig. 4.9. Temporal evolution of the height of particle from the bottom wall.

1. qr ¼ 1:1; kr ¼ 15 and cr ¼ 1 2. qr ¼ 1:6; kr ¼ 5 and cr ¼ 1 3. qr ¼ 100; kr ¼ 15 and cr ¼ 1. We monitored the vertical velocity of the particle as the time progresses. The plots in Fig. 4.8 shows vertical velocity evolution of the particle with respect to time along with the results obtained by Yu et al. [22] and Feng et al. [11] for the respective cases. The graphs of the case (1) and (2) shows that the particle initially accelerates towards the bottom of the channel under the action of gravitational force. However, due to the internal heat source, the temperature of the particle increases and because of heat exchange the temperature of surrounding fluid also increases. The temperature gradient at the boundary layer induced a buoyancy force that balances the gravitational acceleration and then at a certain point it overcomes the force of gravity. Therefore, the particle reverses its direction of motion and begins to accelerate in the upward direction. The particle continues its upward motion until the resistance caused by the upper wall comes into play and slowed it down. The geometry restrictions weaken the natural convection currents which are responsible for the particle’s upward motion. As in case (3), the particle-fluid density ratio is much higher, so the particle motion is dominated by gravitational acceleration. The comparison of our results with available literature results shows a good agreement. Fig. 4.9 shows the spatial and temporal correlation of the particle in Case 2 using three different grid refinements. The ratio of the

area of a particle covered by mesh to the actual area of the particle is denoted by VEF and is equal to 98.46%, 99.16%, 99.53% and 99.84% for refinement lev el  3; lev el  4; lev el  5 and lev el  6 respectively. The plots in Fig. 4.9 show the height of the particle from the bottom wall against time. We noticed that the results of lev el  4 and level  5 are very consistent. To further study the effects of grid refinement on the particle motion, we conduct the simulations of Case 1 using five different refinements of the grid. The plots of vertical velocities of the catalyst particle against time are shown in Fig. 4.10. It is observed that the plot of lev el  4 results is very close to the reference results where the plots of lev el  5 and lev el  6 shows a perfect agreement with the reference results. Also, the difference between the results of lev el  5 and lev el  6 is very insignificant. We conclude that the mesh refinement having element size Dx ¼ 0:02 gives satisfactory results while the desired accuracy can be obtained by using refinement level having element size Dx 6 0:01. Fig. 4.11 shows the temperature contours and position of the particle at different time steps. 4.4. Heat conductivity of suspensions In context of the particulate flow simulation of many particles, not only the fluid-solid interaction effects are important but also

K. Walayat et al. / International Journal of Heat and Mass Transfer 126 (2018) 452–465

Fig. 4.10. Vertical velocity of the catalyst particle w.r.t time using different grid refinement levels.

461

the solid-solid interaction forces play a very important role in capturing the physical behavior of the system. Because the motion of a solid particle in the fluid strongly depends on the drag forces acting on them as well as the particle-particle and particle-wall collision forces. From the literature survey, it is found that there are numbers of promising schemes available for the simulation of fluidparticle interaction problems. However, most of the research works are limited to small scale because of the relatively crude descriptions of the particle-particle collisions and also due to the lack of physical validation. Numerous researchers [18,27,32,40,41] proposed the inter-particle and particle-wall treatments based on the repulsive force model. In all these collision models, a common drawback is that they require many user-defined parameters and a test, trial and empirical approach is required when using in any new configuration even with the same method [42]. Therefore, due to the above-mentioned reasons, the comparison at large-scale may not be possible. However, here we first present a comparison of our results with the results obtained by Feng et al. [34] for two interacting hot particles and

Fig. 4.11. Temperature fields of a catalyst particle at different time steps.

462

K. Walayat et al. / International Journal of Heat and Mass Transfer 126 (2018) 452–465

The values of stiffness parameters for particle-particle and particle-wall collisions in our simulations are selected as

ep ¼ 9:7  104 , e0p ¼ 3:125  102 , ew ¼ 4:8  104 and 2 0 ew ¼ 1:56  10 where the range of repulsive force q ¼ 3:125  102 . Fig. 4.12 shows the dimensionless settling veloc-

Fig. 4.12. Dimensionless settling velocity w.r.t the dimensionless time.

then the results of 172 settling hot particles inside a channel with detailed parameters. In the two-particle case, we used the same physical parameters in our simulations for fluid and particles provided by Feng et al. [34]. But for collision model, we defined our own parameters due to the insufficient or even no information about the collision.

ity of the two particles with respect to the dimensionless time. The present results agree very well with the reference results before and during the collision process but after the collision there is a significant difference in velocities due to the different repulsive forces used. This after collision difference may lead to a slightly different local physical behavior of the system but at large scale engineering problem the system shows much similar overall physical behavior. To study the effects of thermal conductivity on the large number of settling particles, we initially arranged 172 particles of radius a ¼ 0:5 in the form of a circular cluster of radius 7:5 inside a rectangular enclosure of width 20 and height 40. The center of the circular cluster in the channel is at point ð10; 30:5Þ. The fluid-solid density ratio qr ¼ 1:5; and cr ¼ 1, where the governing dimensionless numbers are chosen as Re ¼ 100; Gr ¼ 100 and Pr ¼ 5. The dimensionless temperature T ¼ 1 on the particle surface. Fig. 4.13 shows the results of present simulations where Fig. 4.14 shows the results obtained by Feng et al. [34]. It can be seen in comparison of the results that the main features of the system are qualitatively similar. As the time progresses the particles start to settle due to gravity. Under the influence of wall effects, the particles near the side walls move slower and swept upward while the particles closer to the center-line moves faster. The

Fig. 4.13. Present snapshots at t ¼ 0; 2:5; 5; 7:5; 10 and 20s.

Fig. 4.14. Feng et al. [34] snapshots at t ¼ 0; 2:5; 5; 7:5; 10 and 20s.

K. Walayat et al. / International Journal of Heat and Mass Transfer 126 (2018) 452–465

circular arrangement of particles gradually becomes ellipsoidal in shape and afterward changed to highly irregular shapes. Many complex vortices and buoyancy currents of different shapes and sizes are formed, which are responsible for mixing and boosting the particles. We can notice a significant change in the temperature of the fluid as the hot particles heat up the fluid. However, some regions still appeared blue in color indicating cold or lowtemperature regions. The buoyancy forces keep many of the particles floating for a longer period even after several leading particles settle down in the bottom of the channel. In the simulation of particulate flows, a huge part of CPU time is required for the calculation of hydrodynamic forces and fictitious boundary settings especially in case of large number of particles.

Table 2 Computational details of multigrid FEM-FBM and FEM-FBM. Computational details Number of particles

172 Number of particles

10 100 1000

Total CPU time (hours) Multigrid FEM-FBM

FEM-FBM

3.059

34.17

CPU time for 10-time steps (seconds) Multigrid FEM-FBM

FEM-FBM

520 618 680

1055 6158 56,385

463

Special time efficient techniques are used in the presented multigrid FEM-FBM method to simulate the particulate flows with large number of particles. A detailed discussion about the efficiency and robustness of the present multigrid FEM-FBM is presented by Wan et al. [24]. However, a brief review of the time efficient scheme is presented in the following points  Find particles inside the controlling area of each element on the coarsest mesh.  Use the information obtained from the coarser level on the next finer mesh level because every finer element is within the previous coarser element.  The midpoints of coarser elements become the vertices of the elements of next finer level mesh, assign the obtained information of the midpoints to the vertices of next finer level mesh elements.  Assign special values to a new array on finest mesh i.e. if a node is occupied by ith particle its value is i if not its 0. These hierarchical techniques help to reduce the CPU time for volume integral based force calculations. To show the efficiency and robustness of the present method in handling multiple particles the simulations of the same case with 172 particles are performed on equally refined meshes with and without using timeefficient techniques. For multigrid V-cycle mesh refinement lev el  5 and coarsest lev el  2 are used where without the time efficient multigrid solver the simulations are performed on lev el  5 refinement. The number of elements is Nel ¼ 204; 800 at

Fig. 4.15. Temperature distributions at time t ¼ 0; 1; 2; 3; 4; 5; 6; and 7 respectively while sedimentation of 10,000 particles.

464

K. Walayat et al. / International Journal of Heat and Mass Transfer 126 (2018) 452–465

Fig. 4.16. Temperature distributions at time t ¼ 8; 9; 10; 11; 12; 13; 14; and 15 respectively while sedimentation of 10,000 particles.

refinement lev el  5 and Nel ¼ 3200 at the coarsest lev el  2. The Jacobi relaxation is used in all cases and the computational details of the two solvers are presented in Table 2. Further, in Table 2 we present the comparison of CPU times consumed by both solvers in 10-time steps for different number of particles. The analysis and comparison of the data in Table 2 provided us a confirmation that the multigrid FEM-FBM is an efficient and robust technique. From the data of full simulations of 172 particles, it is observed that the time efficient scheme of multigrid FEM-FBM is ten times efficient then FEM-FBM. There is a significant difference between the CPU times consumed by both solvers with the same accuracy on an equally refined grid. We can see there is a small increment in the CPU time for multigrid FEM-FBM when the number of particles increased but for FEM-FBM the CPU time raises enormously as the number of particles increases.

4.5. Sedimentation of 10000 small particles The intent of this case is to show that this method is suitable for more complex problems and it can handle a very large number of particles very well. We consider 10,000 circular particles in a rectangular cavity of width 5 and height 8. The particles are placed in rows at the top of the cavity and there are 100 particles in each row. The radius of each particle is 0.025 and initially, at time t ¼ 0 the fluid and all the particles are at rest. The range of repulsive force q ¼ 0:01 where the stiffness parameter p ¼ 105 . A uniform fixed mesh is used with mesh element size Dh ¼ 0:00625 and number of elements NEL ¼ 1; 024; 000 at mesh refinement level

level  5. The temperature of the particles T ¼ 1, higher than the fluid. The flow parameters for this case are Re ¼ 100, Gr ¼ 10, Pr ¼ 1, qr ¼ 1:5, m ¼ 0:01 and cr ¼ 1. From Figs. 4.15 and 4.16, it can be observed that, as the particles move downwards they increase the temperature of surrounding fluid and cause a strong fluid motion. The development of Rayleigh-Taylor instability is clearly shown in these pictures. This instability develops into a fingering and text-book phenomenon, and many symmetry breaking and other bifurcation phenomena including drafting, kissing and tumbling take place at various scales in space and time. We can see that many complex vortices of different size have been formed, these vortices pull the particles down and mix them with each other. Some very strong eddies are formed, and we can see that these eddies boost the particles and push them almost back to the top of the channel. After some time, the overall temperature field within the channel turn into same and become stable, which illustrates quasi-steady state of heat transfer between fluid and particles. At the end, all the particles settle down at the bottom of the channel and the fluid again comes to rest. After the particles settled down, one can observe blue region at the top of the channel indicating that the fluid cools down again. The simulations are done on a computer with a 3.60 GHz processor and 12 GB (RAM) and the particles took around 360 simulation hours to settle down completely. 5. Conclusions In this paper, the multi-grid finite element fictitious boundary method for the simulation of fluid-solid two-phase flows with heat

K. Walayat et al. / International Journal of Heat and Mass Transfer 126 (2018) 452–465

transfer is presented. We used this method to resolve the motion of 2-D solid particles along with heat transfer effects in an incompressible Newtonian fluid flow. The Boussinesq approximation is employed to couple the momentum and energy equations. Fictitious boundary method is responsible to handle the velocity and energy interactions between fluid and particle. In comparison to other fictitious domain or Immersed boundary methods, the presented fictitious boundary method along with multi-grid solver is much more robust as the convergence rate of multi-grid solver does not depend on problem size. This method is very easy to implement and flexible for different boundary conditions. In particular for moving boundaries such as particles in fluid, the presented approach exhibits much more promising behavior in terms of efficiency, robustness, flexibility and more important accuracy. First, we simulate a benchmark case with a fixed particle placed eccentrically in a square enclosure and then perform test cases for a particle settling in a vertical channel and study the motion of a particle under the influence of both constant and time-dependent temperature on its surface. We validated our method by comparing our results of test cases with the results reported in the literature. The results obtained using the presented method shows superb agreement with those found in the literature. Further, we simulate more complex problems with large number of particles to test the efficiency of our method and inspect its potential to simulate real particulate flows. Results from the simulation problem of heat conductivity of suspensions and sedimentation of 10,000 particles show that multi-grid fictitious boundary method is computationally efficient, and it can handle a very large number of particles easily. We conclude that the presented method is computationally cheap and its implementation is simple and straightforward, as it treats the fluid part, the calculation of hydrodynamic forces, temperature and the motion of a particle in a sub-sequential way. Conflict of interest The authors declare that there are no conflict of interests. Acknowledgements This work has been supported by the Beijing Innovation Center for Engineering Science and Advanced Technology (BIC-ESAT), and the National Natural Science Foundation of China (Grant Nos. U1530110 and 51779003). References [1] S. Turek et al., The Fictitious Boundary Method for the implicit treatment of Dirichlet boundary conditions with applications to incompressible flow simulations, in: Challenges in Scientific Computing-CISC 2002, Springer, 2003, pp. 37–68. [2] D. Wan, S. Turek, Direct numerical simulation of particulate flow via multigrid FEM techniques and the fictitious boundary method, Int. J. Numer. Meth. Fluids 51 (5) (2006) 531–566. [3] D. Wan et al., An efficient multigrid FEM solution technique for incompressible flow with moving rigid bodies, in: Numerical Mathematics and Advanced Applications, Springer, 2004, pp. 844–853. [4] W. Hackbusch, Multi-grid methods and applications, Springer Science & Business Media, 2013. [5] P. Wesseling, C.W. Oosterlee, Geometric multigrid with applications to computational fluid dynamics, J. Comput. Appl. Math. 128 (1) (2001) 311–334. [6] T.A. Davis, I.S. Duff, A combined unifrontal/multifrontal method for unsymmetric sparse matrices, ACM Trans. Math. Software (TOMS) 25 (1) (1999) 1–20. [7] R. Glowinski et al., A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow, J. Comput. Phys. 169 (2) (2001) 363–426. [8] N.A. Patankar et al., A new formulation of the distributed Lagrange multiplier/fictitious domain method for particulate flows, Int. J. Multiph. Flow 26 (9) (2000) 1509–1524.

465

[9] T.I. Zohdi, An introduction to modeling and simulation of particulate flows, SIAM, 2007. [10] D. Wan, S. Turek, An efficient multigrid-FEM method for the simulation of solid–liquid two phase flows, J. Comput. Appl. Math. 203 (2) (2007) 561–580. [11] Z.-G. Feng, E.E. Michaelides, Heat transfer in particulate flows with direct numerical simulation (DNS), Int. J. Heat Mass Transf. 52 (3) (2009) 777–786. [12] B. Maury, Direct simulations of 2D fluid-particle flows in biperiodic domains, J. Comput. Phys. 156 (2) (1999) 325–351. [13] B. Maury, Characteristics ALE method for the unsteady 3D Navier-Stokes equations with a free surface, Int. J. Comput. Fluid Dyn. 6 (3) (1996) 175–188. [14] H.H. Hu et al., Direct simulation of fluid particle motions, Theor. Comput. Fluid Dyn. 3 (5) (1992) 285–306. [15] H.H. Hu, Direct simulation of flows of solid-liquid mixtures, Int. J. Multiph. Flow 22 (2) (1996) 335–352. [16] J. Feng et al., Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid Part 1. Sedimentation, J. Fluid Mech. 261 (1994) 95–134. [17] C.S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25 (3) (1977) 220–252. [18] R. Glowinski, Finite element methods for incompressible viscous flow, in: Handbook of Numerical Analysis, vol. 9, 2003, pp. 3–1176. [19] R. Glowinski et al., A distributed Lagrange multiplier/fictitious domain method for particulate flows, Int. J. Multiph. Flow 25 (5) (1999) 755–794. [20] H. Gan et al., Direct numerical simulation of the sedimentation of solid particles with thermal convection, J. Fluid Mech. 481 (2003) 385–411. [21] H. Gan et al., Simulation of the sedimentation of melting solid particles, Int. J. Multiph. Flow 29 (5) (2003) 751–769. [22] Z. Yu et al., A fictitious domain method for particulate flows with heat transfer, J. Comput. Phys. 217 (2) (2006) 424–452. [23] S. Kim, S.J. Karrila, Microhydrodynamics: Principles and Selected Applications, Courier Corporation, 2013. [24] I. Demirdzˇic´ et al., Fluid flow and heat transfer test problems for nonorthogonal grids: bench-mark solutions, Int. J. Numer. Meth. Fluids 15 (3) (1992) 329–354. [25] J. Pacheco et al., Numerical simulations of heat transfer and fluid flow problems using an immersed-boundary finite-volume method on nonstaggered grids, Numerical Heat Transf., Part B: Fundam. 48 (1) (2005) 1–24. [26] H. Ström, S. Sasic, DNS of dispersed multiphase flows with heat transfer and rarefaction effects, J. Comput. Multiph. Flows 6 (3) (2014) 193–206. [27] Z.-G. Feng, E.E. Michaelides, The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems, J. Comput. Phys. 195 (2) (2004) 602–628. [28] R.M. MacMECCAN et al., Simulating deformable particle suspensions using a coupled lattice-Boltzmann and finite-element method, J. Fluid Mech. 618 (2009) 13–39. [29] Y. Feng et al., Combined three-dimensional lattice Boltzmann method and discrete element method for modelling fluid-particle interactions with experimental assessment, Int. J. Numerical Methods Eng. 81 (2) (2010) 229– 245. [30] A. Hager et al., Parallel open source CFD-DEM for resolved particle-fluid interaction, J. Energy Power Eng. 7 (9) (2013) 1705. [31] R. Sun, H. Xiao, SediFoam: a general-purpose, open-source CFD–DEM solver for particle-laden flow with emphasis on sediment transport, Comput. Geosci. 89 (2016) 207–219. [32] P. Singh et al., Distributed Lagrange multiplier method for particulate flows with collisions, Int. J. Multiph. Flow 29 (3) (2003) 495–509. [33] K. Walayat, Z. Wang, K. Usman, M. Liu, A multigrid finite element fictitious boundary method for fluid-solid two-phase flows, in: The 8th International Conference on Computational Methods (ICCM2017), 2017, April. [34] Z.-G. Feng, E.E. Michaelides, Inclusion of heat transfer computations for particle laden flows, Phys. Fluids 20 (4) (2008) 040604. [35] J. Kim, H. Choi, An immersed boundary finite-volume method for simulations of the heat transfer in complex geometries, Korean Soc. Mech. Eng. Int. J. 18 (2004) 1026–1035. [36] T. Tsutsumi, S. Takeuchi, T. Kajishima, Heat transfer and particle behaviours in dispersed two-phase flow with different heat conductivities for liquid and solid, Flow Turbulence Combust. 92 (1–2) (2014) 103–119. [37] Z.G. Feng, S.G. Musong, Direct numerical simulation of heat and mass transfer of spheres in a fluidized bed, Powder Technol. 262 (2014) 62–70. [38] S. Haeri, J.S. Shrimpton, A new implicit fictitious domain method for the simulation of flow in complex geometries with heat transfer, J. Comput. Phys. 237 (2013) 21–45. [39] C. Dan, A. Wachs, Direct numerical simulation of particulate flow with heat transfer, Int. J. Heat Fluid Flow 31 (6) (2010) 1050–1057. [40] Z.G. Feng, E.E. Michaelides, Proteus: a direct forcing method in the simulations of particulate flows, J. Comput. Phys. 202 (1) (2005) 20–51. [41] X.D. Niu, C. Shu, Y.T. Chew, Y. Peng, A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating incompressible viscous flows, Phys. Lett. A 354 (3) (2006) 173–182. [42] H. Zhang, Y. Tan, S. Shu, X. Niu, F.X. Trias, D. Yang, H. Li, Y. Sheng, Numerical investigation on the role of discrete element method in combined LBM–IBM– DEM modeling, Comput. Fluids 94 (2014) 37–48.