An improved immersed boundary method with direct forcing for the simulation of particle laden flows

An improved immersed boundary method with direct forcing for the simulation of particle laden flows

Journal of Computational Physics 231 (2012) 3663–3684 Contents lists available at SciVerse ScienceDirect Journal of Computational Physics journal ho...

1MB Sizes 1 Downloads 36 Views

Journal of Computational Physics 231 (2012) 3663–3684

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

An improved immersed boundary method with direct forcing for the simulation of particle laden flows Tobias Kempe ⇑, Jochen Fröhlich Institute of Fluid Mechanics, TU Dresden, George-Bähr Str. 3c, 01062 Dresden, Germany

a r t i c l e

i n f o

Article history: Received 27 April 2011 Received in revised form 1 November 2011 Accepted 13 January 2012 Available online 1 February 2012 Keywords: Multiphase flow Particles Immersed boundary method Euler–Lagrange coupling Collision modeling

a b s t r a c t An efficient approach for the simulation of finite-size particles with interface resolution was presented by Uhlmann [M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, J. Comput. Phys. 209 (2005) 448–476.]. The present paper proposes several enhancements of this method which considerably improve the results and extend the range of applicability. An important step is a simple low-cost iterative procedure for the Euler–Lagrange coupling yielding a substantially better imposition of boundary conditions at the interface, even for large time steps. Furthermore, it is known that the basic method is restricted to ratios of particle density and fluid density larger than some critical value above 1, hence excluding, for example, non-buoyant particles. This can be remedied by an efficient integration step for the artificial flow field inside the particles to extend the accessible density range down to 0.3. This paper also shows that the basic scheme is inconsistent when moving surfaces are allowed to approach closer than twice the step size. A remedy is developed based on excluding from the force computation all surface markers whose stencil overlaps with the stencil of a marker located on the surface of a collision partner. The resulting algorithm is throughly validated and is demonstrated to substantially improve upon the original method. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction A vast number of engineering applications involve flows in which solid particles are transported. These range from mechanical and chemical process engineering to energy conversion, biological flows, environmental flows, etc. The accurate and efficient numerical simulation of this type of flow hence is of substantial importance for fundamental research as well as for understanding and optimizing complex flow problems. The three-dimensional, fully coupled and interface-resolving simulation of flows with a huge number of particles of finite size to date is an area of active research and development. An efficient approach to model this situation is an immersed boundary method (IBM) as originally proposed by Peskin [26] which has since then been substantially advanced as documented in a large number of research papers and several recent reviews [27,13,20,28,47]. The basic idea of an IBM is to employ a numerically efficient Cartesian grid for the discretization of the fluid phase and to represent the immersed fluid–solid interface which does not coincide with any grid line by modifications of the governing equations. This can be formulated as additional source terms in the momentum equation and can be performed in the continuous setting, i.e. before discretization, or by altering the linear system resulting from discretizing in time and space which is termed discrete forcing [20]. The IBM approach is very appealing since it relieves the user from the task of grid generation

⇑ Corresponding author. E-mail address: [email protected] (T. Kempe). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2012.01.021

3664

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

which can be extremely difficult and tedious for complex configurations. Moreover, certain geometries such as closely packed spheres, for example, just can not be handled with body-fitted grids, except if the discretization is fully unstructured, while they pose no problem for an IBM [34]. For immobile fluid–solid boundaries many different variants of an IBM can be applied and generally work well as documented in the cited review papers. Iaccarino and Verzicco [13], for example, focus on methods applied to situations where the immersed boundary does not move with respect to the Eulerian grid. An exception is the work of Verzicco et al. [49], addressed in that review, where the flow in an IC engine with moving piston is simulated with, however, the velocity of the piston being constant and the time step chosen so that the piston surface coincides with a Cartesian grid line in each time step. Le et al. [16] also consider stationary rigid boundaries embedded in the Cartesian grid with local forces applied at the boundaries to impose the no-slip condition. These are computed implicitly by solving a small system of equations at each time step derived from a second order projection method and distributed to the nearby Cartesian grid points. Grid generation is an even more pressing issue for moving boundaries so that an IBM is particularly interesting in such cases. This was in fact the motivation for the earliest IBM proposed by Peskin [26] which is of the continuous-forcing type. For the same reason, several other variants of an IBM for moving boundaries have been proposed in the more recent literature. Tryggvason and co-workers [46,36,28,37] developed a method for multiphase flow computations based on a mobile unstructured two-dimensional grid located at the interface between the phases. A single set of equations is solved on a fixed Eulerian grid with the material properties, like density and viscosity, determined from the phase indicator. The latter is generated based on the position of the interface grid and is smeared out to avoid numerical problems. The forces between the phases are represented by the stress terms in the equations solved on the Eulerian grid, so that this is somewhat different compared to other methods introducing additional localized forces for the coupling of the phases. Interfacial terms resulting from surface tension are accounted for by adding the appropriate sources as delta functions at the immersed boundary. The method, for example, was applied to bubbles with their shape changing as a function of the surrounding flow [46]. Mittal, Udaykumar and co-workers [38,39,20,19] over several papers developed sharp-interface IBMs. In [39], a cut-cell method for stationary and moving boundaries was proposed as well as a remedy for the problem of freshly cleared cells. In [19], a ghost-cell method was developed and applied to moving bodies of very complex shape. Gilmanov and Sotiropoulos [9] also developed an IBM with sharp interface using appropriate reconstruction of the solution at points next to the boundary by means of the information from the boundary itself and points further away. The algorithm also is capable of treating moving boundaries. A similar method was proposed by Roman et al. [31] for curvilinear Eulerian grids, further extended in [30]. Very recently, Cheney and Botella [4] developed an energy conserving cut-cell method in two dimensions, also supplemented with the possibility to account for moving boundaries. The extension to three dimensions seems under way. While allowing for sharp interfaces, cut cell methods are a little cumbersome due to the necessity to handle the different stencils which can occur for geometrical reasons when the interface cuts a hexahedral cell and to treat the freshly cleared cells. One of the reasons why moving boundaries constitute a severe challenge to an IBM [20] is the following: If the local stencil near the boundary changes too abruptly when the boundary moves by a minute amount, substantial oscillations of the computed force exerted by the fluid onto the solid may be generated if the scheme is not sufficiently regularized, as documented for some methods in [42]. For situations where the trajectory of moving objects is prescribed this is uncritical, as with most of the applications cited above where IBMs with sharp interfaces were used. The same holds if freely moving objects are by, say, two or more orders of magnitude larger than the step size of the Cartesian grid. For particle-laden flows, in contrast, the least possible number of grid points per diameter is desired if large numbers of particles are to be discretized. Most of all, the trajectory of the particles is computed as a function of the forces acting upon them, so that abrupt, unrealistic fluctuations of the computed forces are intolerable. A variant to cope with moving bodies in an IBM without this problem is to solve the Navier–Stokes equations in a noninertial coordinate system attached to the moving body, at the price of introducing additional terms in the equations, as proposed in [15]. This indeed can be a solution for single objects but is not applicable to suspensions transporting a large number of moving particles without resorting to Chimera grids. If the price of interpolation between stationary and nonstationary grids is paid, however, it seems more advantageous for simple geometries, such as spheres in a viscous fluid, to employ local body-fitted Chimera grids for the moving bodies as, e.g., proposed in [23], while an IBM would only be applied for very complex particle shapes. If this kind of solution is not desired, it is the continuous forcing approach which allows to obtain a smooth behavior of the fluid forces onto moving particles in a relatively easy way, since regularization of the forces coupling fluid and solid is applied prior to discretization [20]. Based on various tests of existing methods supporting this point of view, Uhlmann [42] developed an IBM meeting these requirements. It is characterized by Lagrangian markers representing the phase boundary and smoothed variants of the Dirac delta function for the transfer of information between Cartesian and Lagrangian points. Furthermore, the algorithm treats the flow field as a constant-density field, hence avoiding performance problems of the Poisson solver due to high density ratios as for example observed by Francois et al. [8]. Instead, the ODE solved for the linear and angular velocities of the moving body accounts for the density of the particle in a global way by means of the resulting inertia and bouyancy forces. No empirical correlations are required for the forces acting on the particle since the interface is fully resolved. The results presented below were obtained in a project devoted to the development of realistic collision models for particle laden flow. The IBM of [42] was selected as a framework and physical collision modeling attempted in various ways. It

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

3665

soon turned out, however, that several discretization issues of the IBM dramatically overwhelmed all these attempts. Hence, before any further development could be continued these had first of all to be remedied. The resulting enhanced algorithm is presented in the present paper. The critical issues are highlighted and illustrated by prototypical situations. Appropriate remedies then are proposed and the simulation of the examples repeated to demonstrate the improvement. The new algorithm shows better performance and extends the parameter range were the method can be applied so that it is of general interest. The modifications described in the present paper are prerequisites for realistic collision models, a topic which will be addressed in a companion paper [14]. The structure of the remaining text is as follows. After recalling the basic algorithm a new forcing scheme is presented which substantially reduces the parasitic velocity at the phase boundary without noticeable increase of the computational cost. Subsequently, a variant of the method is proposed allowing to treat particles with density close to the fluid density and smaller, which was not possible so far. Finally, discretization issues related to the contact of particle surfaces are addressed and a simple but efficient enhancement of the scheme is discussed for this situation. Generally, there is no restriction for the shape of the interface, but here we focus on spherical geometries for ease of presentation. 2. Basic immersed boundary method In this section we describe the fluid solver and recall the method of Uhlmann [42] which is the starting point. The equations to be solved are the unsteady three-dimensional Navier–Stokes equations for Newtonian fluids of constant density. The momentum equation reads

@u þ r  ðuuÞ ¼ r  s þ f; @t

ð1Þ

where s is the hydrodynamic stress tensor divided by the fluid density qf, i.e.

s¼

p

qf

I þ mðru þ ðruÞT Þ:

ð2Þ

The continuity equation reads

r  u ¼ 0:

ð3Þ T

Nomenclature is as usual, with u = (u, v, w) designating the velocity vector in Cartesian components, i.e. along the Cartesian coordinates x, y, z, while p represents the pressure, f = (fx, fy, fz)T a volume force, I the identity matrix, m the kinematic viscosity of the fluid, qf the density of the fluid, and t the time. The spatial discretization of (1)–(3) is performed by a second-order finite-volume scheme on a staggered grid [12], as illustrated in Fig. 1a. The Cartesian grid is supposed to have the same constant step size h in all three directions. The fluid solver alone has been throughly validated and convergence of second order in space and time has been checked. Coupling of the fluid and the solid phase in an IBM is realized by inserting additional volume forces at or close to the interface. In the present approach, the fluid–solid interface is represented by discrete surface markers located on the phase boundary. Since they move with the particle they are termed Lagrangian points. Fig. 1 provides an illustration. The transfer of information between Eulerian and Lagrangian points and vice versa is performed by interpolation which is implemented

v p

u

y x

(a)

(b)

Fig. 1. Sketch of discretization and interface representation employed by the immersed boundary method. (a) Staggered Cartesian grid and embedded interface with marker points in the two-dimensional case. (b) Surface of a three-dimensional spherical particle represented by 1000 marker points using the method of Leopardi [17].

3666

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

via a weighted sum of regularized Dirac delta functions dh and is given by (7b) and (7d) below. The three-dimensional function dh is generated by the tensor product of three one-dimensional functions d1D h applied in each coordinate direction, 1D 1D dh ðx  XÞ ¼ d1D h ðx  XÞdh ðy  YÞdh ðz  ZÞ

ð4Þ

where

d1D h ðx  XÞ ¼

1 uðrÞ h

ð5Þ

is the regularized one-dimensional delta function. It is based on the continuous function u, which in the present implementation is defined by the three-point version of Roma et al. [29]

8  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > 1 > 5  3krk  3ð1  krkÞ2 þ 1 ; 0:5 6 krk 6 1:5 > >6 < qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uðrÞ ¼ 1 2 > krk 6 0:5 > 3 ð1 þ 3krk þ 1Þ; > > : 0; otherwise;

ð6Þ

with r = (x  X)/h and hence is 3h wide. It is shown in Fig. 2 and compared to the four point variant of [26]. The function (6) was selected here for its smaller support as this provides a good balance between numerical efficiency and smoothing properties. The present work is concerned with spherical particles. While Uhlmann [42] used an iterative method to equidistribute a given number NL of Lagrangian points on the surface of the sphere, we employ here the non-iterative method of Leopardi [17] illustrated by the distribution of NL = 1000 markers in Fig. 1b. For later reference the basic scheme is recalled here in equations (7a)–(7m). As in [42], capital letters designate values at the Lagrangian marker points while small letters identify values at the Eulerian grid points.

 k1  ~  uk1 p u  ck r  ðuuÞk1  fk r  ðuuÞk2 ¼ 2ak mr2 uk1  2ak r Dt q UC ðXl Þ ¼

Ny X Nz Nx X X i¼1

FðXl Þ ¼

ð7aÞ

~ ðxi;j;k Þdh ðxi;j;k  Xl Þh3 u

ð7bÞ

j¼1 k¼1

Ud ðXl Þ  UC ðXl Þ Dt

fðxi;j;k Þ ¼

Np X Nl X

ð7cÞ

FðXl Þdh ðxi;j;k  Xl ÞDV l

ð7dÞ

p¼1 l¼1

r2 u 

  ~ u 1 u ¼ þ f þ r2 uk1 ak mDt a k m Dt

ð7eÞ

r2 /k ¼ r  u

ð7fÞ

uk ¼ u  r/k

ð7gÞ

(a)

(b)

Fig. 2. One dimensional function u used for the construction of regularized Dirac functions. (a) Function used by Peskin [27], (b) function proposed by Roma et al. [29]. Grid points of the Eulerian grid are located at distances of unity.

3667

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

pk ¼ pk1 þ

/k 1  mr2 /k 2ak Dt 2

ukp ¼ uk1  2ak Dt p

ð7hÞ

Nl X qp qf FðXl ÞDV l þ g mp ðqp  qf Þ l¼1

xkp ¼ xpk1  2ak Dt

! ð7iÞ

Nl   X qp qf Xl  xkp  FðXl ÞDV l Ip ðqp  qf Þ l¼1

ð7jÞ

  þ ak Dt ukp þ upk1 xkp ¼ xk1 p   ukp ¼ uk1 þ ak Dt xkp þ xk1 p p   Ud ðXl Þ ¼ ukp þ xkp  Xl  xkp

ð7kÞ ð7lÞ ð7mÞ

The time advancement is accomplished by an explicit third-order low-storage Runge–Kutta scheme for the convective terms ~ without the momentum source and by a Crank–Nicolson scheme for the viscous terms. First, an intermediate velocity field u ~ is interpolated to the interface C via the regf is computed using fully-explicit viscous terms (7a). In the next step, (7b), u ularized Dirac function (4)–(6) yielding UC. The Lagrangian interface force F then is determined directly at the surface makers by the so-called direct forcing method of Mohd-Yusof [21]. With direct forcing, F is given by the difference of the desired velocity Ud at the interface to the actual velocity UC at the interface, divided by the time step Dt, (7c), with Ud resulting from the equations of motion for the particle according to the previous Runge–Kutta step. Subsequently, F is spread to the grid points, (7d), and the resulting Eulerian force f is introduced in the solution of a Helmholtz Eq. (7e) which is required for the implicit treatment of the viscous terms, yielding the velocity field u⁄. Eq. (7e) is derived by taking the difference between (7a), here with explicit viscous terms, and the same equation but now with semi-implicit viscous terms and IBM force f included. The solution of the pressure Poisson Eq. (7f) and a projection step yields the divergence-free velocity field uk of that Runge–Kutta sub-step (7g). The related pressure field is given by Eq. (7h). The Runge–Kutta coefficients used here are [33]

a1 ¼ 4=15 a2 ¼ 1=15 c1 ¼ 8=15 c2 ¼ 5=12 f1 ¼ 0

a3 ¼ 1=6 c3 ¼ 3=4

ð8Þ

f2 ¼ 17=60 f3 ¼ 5=12:

In contrast to [42], the present implementation makes use of the library hypre [5,6] to solve (7e) and (7f). The second element of the method is constituted by the equations of motion of the particle. Ordinary differential equations are solved for the translation velocity up of the particle (7i) and for its angular velocity xp (7j). Here, mp is the particle mass and Ip the moment of inertia, which for spherical particles is given by Ip ¼ 8pqp R5p =15 with qp and Rp being particle density and radius, respectively. The position of the center of gravity of the particle xp and the rotation angle up are computed by (7k) and (7l), respectively. Finally, the velocity Ud of the surface markers positioned at Xl is found from the superposition of translational and rotational motion (7m). The resulting scheme is versatile and extremely efficient so that simulations with very large numbers of particles are possible as, e.g., performed in [44]. The definition of surface markers and not volume markers all over the particle volume allows for flow also in the interior of the particle since the Navier–Stokes equations are solve throughout on the Eulerian grid as illustrated in Fig. 3b below. This is a well known feature of the method and tolerated as long as the boundary conditions at the particle surface are correctly fulfilled [41,42].

y − periodic 0.8

θ

Ub

y

y / Ly

Outlet

Inlet

0.6 0.4 0.2

x

y − periodic 0

(a)

0.2

0.4

0.6

x / Lx

0.8

1

(b)

Fig. 3. Flow over a two-dimensional stationary cylinder. (a) Sketch of problem and boundary conditions, (b) streamlines for Re = 10 using Grid 1.

3668

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

3. Improved forcing scheme 3.1. Basic method and analysis Numerical experiments in [42] and own simulations showed that the original method of [42] does not exactly impose the desired interface velocity Ud on the fluid. Rather, a velocity error e = jUd  UCj remains, the magnitude of which depends on the configuration. One reason for the velocity error is the explicit forcing procedure (7c) as will be explained in the following. To derive Eq. (7c), the momentum Eq. (1) is considered in dimensionless form and discretized in time with semi-implicit treatment of the viscous term, as used in (7e). The pressure as well as the convective terms are abbreviated here by Ck for simplicity. The provisional velocity at an Eulerian point then is given by

~ ¼ uk þ Dt u

a

k

Re



~ þ uk Þ þ Ck þ f ðiÞ ; r2 ðu

ð9Þ

where is the force applied in the time interval with the upper index referring to the implicit scheme used here. This equation ~ d at the grid point reading can be transformed into an equation for the volume force f(i) required to yield a desired velocity u

f

ðiÞ

¼

 ~ d  uk  a k 2 u ~ þ uk Þ þ Ck : r ðu  Dt Re

ð10Þ

The numerical evaluation of this equation, here referred to as implicit forcing, is demanding as will be detailed in Section 3.2 below. For that reason an explicit forcing scheme is employed in [42]. To derive the explicit scheme, Eq. (1) is now considered with an explicit treatment of the viscous term using the same notation as above. The provisional velocity at the new time level then is given by

~ ¼ uk þ Dt u

  2ak 2 k r u þ Ck þ f ðeÞ ; Re

ð11Þ

~ d to be yielding the volume force f(e) of the explicit scheme needed to impose u

f

ðeÞ

¼

  ~ d  uk 2ak 2 k u r u þ Ck :  Dt Re

ð12Þ

Evaluating (12) at the Lagrangian marker points yields (7c). If f(e) according to (12) is used in (9) instead of (10), a forcing error Df = f(i)  f(e) exists. This error yields an error in the ~ d is not exactly achieved. The difference imposition of the velocity boundary conditions at the interface, so that the velocity u ~ in the resulting equation, the forcing error then Df is given by the difference of Eqs. (10) and (12). Using (11) to replace uk  u is given by



Df ¼ Dt

ak

Re

r2 ðCk þ f ðeÞ Þ þ

2a2k 2

Re



r 2 r 2 uk

 :

ð13Þ

Obviously, Df is a function of the time step Dt and the Reynolds number Re. Hence, even steady state results will be different when the size of the time step is changed, as illustrated below. The issue just described is demonstrated using a two-dimensional test case, the flow around a fixed circular cylinder sketched in Fig. 3. The computational domain X = [0; Lx]  [0; Ly] with Lx = Ly is discretized with Nx  Ny = 128  128 points (Grid 1) or Nx  Ny = 256  256 (Grid 2). On the cylinder surface, NL = 97 marker points are used with Grid 1 and NL = 193 with Grid 2, respectively. The boundary conditions are periodic in y-direction. At the inlet a constant velocity in x-direction with magnitude Ub is prescribed, and a convective boundary condition is used at the outlet. In the simulations L = 1, Ub = 1 so that the results reported can be considered dimensionless. The diameter of the cylinder is Dp = 0.24L. On purpose, the ratio Dp/L was not chosen smaller, i.e. as small as would be used to approximate a cylinder in infinite flow, since this yields a handy test case with low cost perfectly fulfilling the present needs. Now the dependency of the velocity error upon the physical and numerical parameters is investigated. In Fig. 4a, the velocity in streamwise direction along the center line at y = Ly/2 is plotted. Furthermore, the velocity error e along the cylinder surface is shown in Fig. 4b. Various values for the time step are used, here characterized by the CFL number. The Reynolds number is Re = 10. Obviously, the no-slip condition is not exactly imposed. Instead, the error is about 13% at the stagnation point for CFL = 1 and hence is substantial. Fig. 4a shows that this error also affects the wake. The error decreases with the time step, as predicted by Eq. (13). Additional simulations reported in Fig. 5a show the flow with Re = 40 for CFL = 1, 0.6, and 0.1. Comparing both solutions for CFL = 1, the error at the stagnation point is reduced from 11.2 to 8.2% for the larger Reynolds number. This illustrates its decrease with Re. The computation of e for the steady state solution in Figs. 4 and 5 ~ discussed in this involves a large number of timesteps with linear and non-linear operations of which the determination of u section is only one. Hence, while Eqs. (11) and (13) suggest a reduction of e proportional to (Dt)2 and 1/Re, the other components of the algorithm, like projection, the evaluation of the non-linear convective terms, etc., alter the dependency. Important is the demonstration that the error indeed depends on Dt and Re. It is reduced with the time step and with increasing Re (see Fig. 6).

3669

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

1

15 CFL = 1.0 CFL = 0.6 CFL = 0.1

CFL = 1.0 CFL = 0.6 CFL = 0.1

0.6

ε |u| [%]

u / Ub

0.8

0.4

10

5 0.2 0

0

0.2

0.4

0.6

0.8

x / Lx

0

11

0

100

(a)

200

θ [deg]

300

(b)

Fig. 4. Flow over a circular cylinder with Re = 10 using the coarse grid with 1282 points and various CFL numbers. (a) Streamwise velocity component at y = Ly/2 with vertical lines identifying front and back of the cylinder, (b) velocity error at the interface in percentage of Ub.

15

15 CFL = 1.0 CFL = 0.6 CFL = 0.1

10

ε|u| [%]

ε |u| [%]

10

5

5

0

Grid 1 Grid 2

0 0

100

200

θ [deg]

(a)

300

0

100

200

θ [deg]

300

(b)

Fig. 5. Velocity error at the interface for the flow over a circular cylinder with different para-meters. (a) Re = 40 and various time steps. (b) Impact of grid refinement for Re = 10 using Grid 1 with 1282 points and Grid 2 with half the spacing, CFL number 0.6 in both cases.

Fig. 6. Example for the coupling of grid points in the vicinity of the interface with the implicit forcing procedure. Due to the combined interpolation and spreading steps the grid points highlighted by the two circles are now coupled, for example.

Let us now address the dependency of the velocity error at the interface on the spatial discretization. To this end the same configuration is simulated with the grid refined by a factor of two in each spatial direction using Grid 2 instead of Grid 1. The grid and the time step are refined simultaneously, i.e. the condition CFL = 0.6 is maintained. Fig. 5b displays the results of this exercise. The velocity error remains nearly constant for the fine grid, regardless of the smaller time step. Obviously, the error can not be reduced by spatial grid refinement without simultaneous reduction of the time step to maintain the CFL number.

3670

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

Table 1 Performance of the original forcing scheme (nf = 0) and the new forcing scheme (nf > 0), measured by the maximum error of the streamwise velocity on the cylinder surface.

Stationary

Stationary Stationary Accelerated

Grid

Rep

nf

eCFL=0.1 [%]

eCFL=0.6 [%]

eCFL=1.0 [%]

1 1 1 1 1 1 1 2 2 1 1

10 10 10 10 10 10 40 10 10 10 10

0 1 2 3 4 10 0 0 3 0 3

1.270 0.380 0.153 0.068 0.039 0.022 0.856 1.346 0.043 – –

6.981 – – – – – – – – 10.693 0.445

11.219 3.421 1.387 0.608 0.257 0.264 8.238 11.728 0.571 – –

Only if the time step is reduced further yielding a lower CFL number the error is reduced. If the CFL number is increased the error on the fine grid becomes even larger than on the coarse grid. Quantitative data are reported in Table 1. The reason for the observed problem is the inequality between interpolation (7b) and spreading (7d) in the discrete case. To illustrate this issue, the Lagrangian force F is spread to the Eulerian points and interpolated back to the Lagrangian points. Spreading

fðxi;j;k Þ ¼

NL X

FðXl Þdh ðxi;j;k  Xl ÞDV L ;

ð14Þ

l¼1

and interpolation yields the values F⁄ at the Lagrangian points with

F ðXl Þ ¼

Ny X Nz Nx X X i¼1

3

fðxi;j;k Þdh ðxi;j;k  Xl Þh :

ð15Þ

j¼1 k¼1

Inserting (14) into (15) gives



F ðXl Þ ¼

NL X m¼1

FðXm Þ

Ny X Nx X Nz X i¼1

! 3

dh ðxi;j;k  Xm Þdh ðxi;j;k  Xl Þ DV L h :

ð16Þ

j¼1 k¼1

Note that, since the functions dh are regularized, the term in brackets is unequal to dh(Xm  Xl). Consider for illustration a P setting with a single Lagrange point (NL = 1, X1 = X). Then, ijk ðdh ðxi;j;k  XÞÞ2 ¼ 1 would be required for F⁄(X) = F(X), to hold, P but the construction only warrants this relation without the square, i.e. ijk dh ðxi;j;k  XÞ ¼ 1, so that F⁄(X) – F(X).

3.2. Implicit forcing The potentially more accurate implicit forcing (10) might be accomplished together with the solution of the Helmholtz equation for the viscous terms (7e). The still unknown Eulerian force f in (7e) then has to be replaced by the interpolation of the fluid velocity from the Eulerian grid to the particle interface (7b), the direct forcing (7c), and the spreading of the interface force back to the Eulerian grid (7d). The force f in (7e) then can be replaced by

! Np X Ny X NL Nz Nx X X X 1 3 d  fðxi;j;k Þ ¼ u ðxi;j;k Þdh ðxi;j;k  Xl Þh dh ðxi;j;k  Xl ÞDV l : U ðXl Þ  Dt p¼1 l¼1 i¼1 j¼1 k¼1

ð17Þ

where p refers to the index of each of the totally Np particles present in the computational domain with this index being dropped for the quantities in the sum for clarity. Using (7e) in combination with (17) might be a means to achieve higher accuracy for the no-slip condition since the forcing scheme is now consistent with the time-discrete Navier–Stokes equations. Unfortunately, f in (7e) now contains u⁄, so that the operator to be inverted becomes very complex, involving interpolation, forcing and spreading. This rules out a direct inversion as the resulting stencil is extremely large and irregular, changing in each Runge–Kutta step in case the particle moves. An iterative scheme based on the forward evaluation of all these steps, iterating over (17) and (7e), can still be applied, but is not performed here for cost reasons.

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

3671

3.3. New forcing scheme To achieve a numerically efficient and more accurate imposition of the interface boundary conditions an alternative forcing method is proposed here. It is obtained by closely investigating the various stages of the numerical solution procedure (7). In Fig. 7, the error is shown first for the preliminary velocity field without volume force (7a), then after the forcing step and the solution of the Helmholtz Equation (7e), and finally after the solution of the Poisson equation and the projection step (7g). The maximum interface velocity error is about 10% at the end of the time step. Obviously, the interface error is not zero after the solution of the Helmholtz equation. This is due to the forcing error (13) and the inequality between interpolation an spreading described above with the help of (16). A further observation is that the error is not substantially modified due to the projection of the velocity onto a divergence-free field. These results suggest the following road to improvement. Additional forcing loops are performed after the solution of the Helmholtz equation but before the projection step. From a formal point of view, the Helmholtz equation now needs to be solved again in order to account for the interaction of the forcing and the viscous terms. Since the correction is assumed to be small in a Runge–Kutta sub-step this is omitted here. The overall numerical scheme performed in each Runge–Kutta step then is given by 1. 2. 3. 4.

~ without forcing (7a) Computation of the preliminary velocity field u Direct forcing and solution of the Helmholtz Eq. (7e) to obtain u⁄ u(0) = u⁄ Additional forcing loop (the iteration index m is dropped for intermediate quantities and only retained for the unknown u) for m = 1, nf P x PNy PNz ðm1Þ 3 (a) Interpolation UC ðXl Þ ¼ Ni¼1 ðxi;j;k Þdh ðxi;j;k  Xl Þh k¼1 u j¼1 d C (b) Forcing F = (U  U )/Dt PN l (c) Spreading fðxi;j;k Þ ¼ l¼1 FðXl Þdh ðxi;j;k  Xl ÞDV l (d) Velocity correction u(m) = u(m1) + Dt f u ¼ uðnf Þ

5. Solution of the Poisson Eq. (7f), projection step (7g), and pressure correction (7h). In this way the interface velocity is imposed accurately while a divergence-free velocity field is obtained at the end of the time step not requiring additional projection steps which are very costly. This is ensured since the projection step is performed only after the additional forcing loops. The result with the modified scheme is displayed in Fig. 8a using the same configuration as above. Without additional forcing the velocity error at the interface strongly depends on the time step. With the new forcing scheme the solution obtained with CFL = 1 and additional forcing is very close to the one for CFL = 0.1 without additional forcing. The data in Table 1 quantify the improvement achieved. For typical configurations, 2 to 3 additional forcing loops are sufficient to obtain a fairly accurate imposition of the no-slip condition. The computational cost generated by the additional forcing with nf = 3 on an SGI Altix is about 6.7  102 s of the total cost of 7.93  101 s for one time step with Re = 10 on Grid 2, hence 8% in that case. This has to be compared to the cost of the projection step within each time step is 4.31  101 s for this configuration. The modified scheme is now applied to an accelerated cylinder, a case which is highly relevant for the treatment of particle collisions, since during the collision process very large changes in magnitude and direction of the particle velocity have to be coped with. In this simulation the cylinder is first placed in steady flow until convergence, then constantly accelerated from rest to up = Ub/2 over the interval t 2 [0, T] with T = 5.4  103L/Ub. The numerical errors become large with the conventional scheme, while the problem is solved with the new scheme, as shown in Table 1 and Fig. 8b.

40

ε |u| [%]

30

20

10

0

0

100

200

θ [deg]

300

Fig. 7. Velocity error in the flow over a circular cylinder with Re = 10 and CFL = 1.0. : preliminary velocity field without volume force, – – – –: after the forcing step and the solution of the Helmholtz equation, ---: after the solution of the Poisson equation and the projection step.

3672

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

20

20

n f = 0 , up = 0

CFL = 1.0 / nf = 0 CFL = 0.1 / nf = 0

15

10

ε |u| [%]

ε|u| [%]

CFL = 0.1 / nf = 3

nf = 0 , up = -0.5 nf = 3 , up = -0.5

10

5

5

0

n f = 3 , up = 0

15

CFL = 1.0 / nf = 3

0

100

200

θ [deg]

(a)

300

0

0

100

200

θ [deg]

300

(b)

Fig. 8. Performance of the modified forcing scheme with (nf = 3) and without (nf = 0) additional forcing loops. (a) Velocity error at the interface for the flow over stationary cylinder at Re = 10. (b) Velocity error at the interface for the flow over a rapidly accelerated cylinder with Re = 10 at t = 0 and CFL = 0.6 accelerated from up = 0 to up = 0.5 m/s at t = 5.4  103s.

Let us put the present results in relation to the study of Uhlmann [43] where the error in the imposition of the no-slip condition at the particle surface was discussed as well. In that reference, three sources for the velocity error at the interface were named. (1) The force term is determined from a fully explicit discretization but a semi-implicit discretization of the viscous term is used. In that study the error was observed to scale with Dt and 1/Re which confirms the present results in Section 3.1. (2) The spatial interpolation introduces an error of order Dx. This is not confirmed here, since the error is seen to be fairly independent of the resolution of the sphere in terms of Dp/h. (3) An error introduced by the projection step due to the fact that the IBM forces are not solenoidal. This can not be confirmed here. In the present study, the error does practically not increase in the projection step. In [43], volume forcing by positioning Lagrangian marker points over the entire particle volume and not only over its surface was found to reduce the velocity error. This was not considered here as the required number of forcing points then scales as (Dp/h)3 instead of (Dp/h)2. Furthermore, the stability of the method is affected requiring further measures for stabilization [41]. The results presented in Table 1 are obtained without solving the Poisson Eq. (7f) and performing the projection step (7g) for the inner iterations. For the case with Grid 1, CFL = 1.0, and nf = 3 inner iterations, the velocity error is 0.608%, while the divergence of the velocity field in the L2 norm (19) is 1.233  106 and the CPU time is 0.098 s. Inserting a projection step after each iteration on the forces substantially increases the computational cost (a factor of 2 was examined with nf = 3) and is hence not retained in the final method. 3.4. Convergence A refinement study is now performed employing the reference scheme and employing the new scheme to assess its convergence behavior. The setup for the flow around a fixed cylinder described above is used again for this purpose. Spatial and temporal step size are reduced simultaneously. The CFL number is set to 0.75 for all cases. Since an analytical solution is not available, the numerical results are compared to reference data obtained on a very fine grid with Nx  Ny = 512  512 points, so that the following definitions of the L1-error and L2-error will be used:

n o ui;j  uref i;j vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi uPN PN  ref y x u u  u t i¼1 j¼1 i;j i;j ¼ Nx Ny

1 ¼ max

ð18Þ

2

ð19Þ

The order of convergence, s, is estimated by the difference between solutions obtained on systematically refined grids using the expression [7]



log



ð2hÞ =ðhÞ log 2

 ;

ð20Þ

where (h) and (2h) are the errors on a grid with the step size h and 2h, respectively. Table 2 shows that the convergence rate of the original scheme is about 1.5. It is smaller than 2 which is obtained for the basic fluid solver without IBM. This is a consequence of the continuous forcing approach and is a well-known fact [3,18]. A comparison of the results in Tables 2 and 3 reveals that the order of the new scheme is the same.

3673

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

Table 2 Grid refinement study for the flow around a stationary cylinder without additional forcing loops. The errors 1 and 2 are defined by (18) and (19), respectively, and the corresponding order by (20). h

ðhÞ 1

ðhÞ ð2hÞ 1 =1

s1

ðhÞ 2

ðhÞ ð2hÞ 2 =2

s2

1/16 1/32 1/64 1/128 1/256

2.32e01 7.81e02 3.28e02 1.32e02 4.20e03

– 2.97 2.38 2.49 3.14

– 1.57 1.25 1.32 1.65

1.02e01 3.77e02 1.55e02 6.35e03 2.03e03

– 2.70 2.43 2.45 3.13

– 1.43 1.28 1.29 1.64

Table 3 Grid refinement study for the flow around a stationary cylinder with the same parameters as in Table 2 but with nf = 3 additional forcing loops. h

ðhÞ 1

ðhÞ ð2hÞ 1 =1

s1

ðhÞ 2

ðhÞ ð2hÞ 2 =2

s2

1/16 1/32 1/64 1/128 1/256

2.82e01 7.88e02 3.33e02 1.28e02 4.10e03

– 3.58 2.36 2.61 3.12

– 1.84 1.24 1.38 1.64

1.08e01 3.72e02 1.54e02 5.88e03 1.89e03

– 2.90 2.42 2.61 3.10

– 1.54 1.28 1.39 1.63

4. Improved method for a single mobile particle in unbounded fluid While in the previous sections we considered a stationary particle, or one with imposed velocity, we now turn to particles freely moving according to the forces acting upon them. 4.1. Basic method The motion of solid particles (index p) in a flow can be described by the trajectory of the center of mass xp(t) and the orientation, specified by appropriate angles. The equations of motion for an individual particle are then obtained from the conservation of linear momentum

mp

dup ¼ qf dt

I Cp

s  n ds þ V p ðqp  qf Þg

ð21aÞ

and angular momentum

Ip

dxp ¼ qf dt

I

r  ðs  nÞ ds:

ð21bÞ

Cp

Here, Cp is the surface of the particle, g the gravitational acceleration, n the outward-pointing normal vector to the surface Cp, and r = x  xp the position vector of the surface point with respect to the center of mass of the particle. Let us mention that for spherical particles it is not necessary to track the angular position by integrating xp in time. Determining the angular velocity from (21b) is sufficient and is more economic. The Lagrangian markers located at the surface of the spherical particles then do not move with the rotation of the particle, only appropriate velocities are imposed at these marker points. The straightforward evaluation of the surface integrals in (21) requires considerable numerical effort, since the six components of the stress tensor must be interpolated onto the fluid–solid interface and integrated by a quadrature rule. With the present IBM, this approach furthermore is inaccurate as gradients are smeared by the regularized Dirac functions in (7d). Uhlmann [40] therefore proposed another method, which is accurate and numerically more efficient. The momentum balance of the fluid enclosed in an arbitrary volume X with surface C moving with the flow is given by

d dt

Z X

qf udV ¼ qf

Z X

f dV þ qf

I

s  nds;

ð22aÞ

C

where f is the body force acting in the control volume. The corresponding angular momentum balance is

d dt

Z X

qf r  udV ¼

Z X

qf r  fdV þ

I C

qf r  ðs  nÞds:

ð22bÞ

This is applied to a fluid control volume coinciding with the particle volume, i.e. X = Xp with surface C = Cp. The hydrodynamic forces and moments in (21) are now replaced by the corresponding expressions in (22) yielding the modified linear equation of motion

3674

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

dup d ¼ dt dt

mp

Z Xp

qf udV  qf

Z Xp

fdV þ V p ðqp  qf Þg

ð23aÞ

qf r  fdV:

ð23bÞ

and angular equation of motion

Z

dxp d ¼ dt dt

Ip

Xp

qf r  u dV 

Z Xp

Assuming the fluid motion inside the enclosed volume being equal to rigid-body motion yields

d dt

Z Xp

qf u dV ¼ qf V p

dup dt

ð24aÞ

and

d dt

Z Xp

qf r  u dV ¼

qf Ip dxp : qp dt

ð24bÞ

Substituting (24) into (23) gives the following approximate equations of motion for particles [42]

qp qf dup ¼ dt mp ðqp  qf Þ qp qf dxp ¼ dt Ip ðqp  qf Þ

!

Z

f dV þ g

ð25aÞ

Xp

Z

r  f dV:

ð25bÞ

Xp

Since the interpolation kernel in (7b) and (7d) conserves linear and angular momentum [29], the integrals over the Eulerian forces in (25a) and (25b) can be replaced by the corresponding summations over the Lagrangian forces at the marker points,

Z

f dV ¼

Nl X

Xp

FðXl ÞDV l

ð26Þ

Nl   X Xl  xkp  FðXl ÞDV l

ð27Þ

l¼1

and

Z

r  f dV ¼

Xp

l¼1

yielding (7i) and (7j), respectively. The derivation shows that this method can be used for arbitrary particle shapes provided that the correct moment of inertia is used. The rigid body assumption (24) introduces a singularity in the equations of motion for qp/qf = 1. Furthermore, the scheme is shown to be stable only for density ratios qp/qf > 1.2 in [42]. Our tests even yielded a slightly higher limit. Neutrally buoyant particles and particles which are lighter than the surrounding fluid hence can not be computed with the basic method. 4.2. Modified scheme for light particles For cases where the weak, explicit coupling is unstable Uhlmann [41] suggests to use iterations within each Runge–Kutta step to enhance fluid–solid coupling, together with an under-relaxation of all variables. Results were not presented, however. Since the approach involves frequent costly projection steps, a different strategy is proposed here. To enhance the parameter range where the present IBM can be employed the rigid body assumption (24) is replaced by the numerical evaluation of the volume integrals in (23) thus avoiding the problematic factor 1/(qp  qf) in (25). Second-order midpoint quadrature rules are employed reading

Z

u dV ¼

Ny X Nz Nx X X

Xp

Z

i¼1

r  u dV ¼

Xp

ui;j;k

V cell i;j;k ai;j;k

ð28aÞ

j¼1 k¼1 Ny X Nz Nx X X i¼1

ri;j;k  ui;j;k

V cell i;j;k ai;j;k :

ð28bÞ

j¼1 k¼1

Here, a is the volume fraction, defined in each grid cell as the fraction of the particle volume contained in the cell volume, i.e.

ai;j;k ¼

V pi;j;k V cell i;j;k

;

ð29Þ

with V cell i;j;k designating the volume of the Cartesian cell at position i, i, k in the Eulerian grid, with appropriate shifts due to the staggered mesh and V pi;j;k . The intersection of Xp with V cell i;j;k . The rate of change of the linear and the angular momentum of the fluid in the k-th Runge–Kutta sub-step then is given by the finite-difference approximations in time

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

"

d dt

#k

Z

u dV

R

R uk dV  uk1 dV 2ak Dt

¼

Xp

3675

ð30aÞ

and

"

d dt

#k

Z

r  udV

R ¼

Xp

R r  uk dV  r  uk1 dV ; 2ak Dt

ð30bÞ

respectively. Inserting (30) into (23) yields the modified equations of particle motion

0 ! 1 " Z #k Nl X ukp  upk1 2ak q d f @qf gA FðXl ÞDV l þ qf udV þ 1  ¼ dt Xp Dt mp qp l¼1

ð31aÞ

and

xkp  xk1 p Dt

0 " Z #k 1 Nl   X 2ak qf d k @ Xl  xp  FðXl ÞDV l þ r  udV A; ¼ dt Xp Ip l¼1

ð31bÞ

with F according to (7c). 4.3. Computation of the particle volume fraction Since the center coordinates of the particle and the corners of each cubic control volume are given, the particle volume fraction ai,j,k can in principle be determined exactly if the particle surface is given in analytical form. Another approach is a cut-cell approximation, where the particle surface is assumed to be piecewise planar. The intersections of the approximating planes with the edges of the control volume are computed and the type of the cut-cell is identified. Subsequently, the volume fraction is determined. For arbitrarily aligned three-dimensional cells, 64 cut-cell types are found which can be reduced to five generic types [48,24]. One can also employ a finer background mesh for this step, together with a stair-step approximation of first order. In [1], a refinement by a factor of 5 is used so that 53 = 125 values inside each computational cell crossing the interface have to be determined and integrated to compensate the lower order. Since the values of ai,j,k are required in each Runge–Kutta sub-step, and possibly for a very large number of particles if a dense suspension is considered, the two approaches addressed are numerically expensive. A different scheme is hence proposed here which is based on approximating the volume fraction by employing the signed-distance level-set function / of the solid–fluid interface. The sign of / at the corners of a three dimensional control volume with interface intersection is shown in Fig. 9. Inside the solid, / is negative while being positive outside. A simple approximation of the volume fraction then can be obtained by

ai;j;k ¼

P8

 /m Hð/m Þ ; P8 m¼1 j/m j

m¼1

ð32Þ

φ8 > 0 φ >0 6

φ5 > 0

Solid

φ7 > 0 Fluid

r

φ4 < 0 φ1 < 0

φ <0 3 φ2 < 0

Fig. 9. Use of the level-set function / for the determination of the volume fraction.

3676

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

α=1

α = 0 ... 1

α = 0.5

α=0

Fig. 10. Cases where the proposed method to evaluate the volume fraction returns the exact value.

where / is computed at the eight corners of each control volume and where H is the Heaviside function

Hð/Þ ¼

0;

/60

ð33Þ

1; / > 0:

The evaluation of ai,j,k according to (32) is very simple as it only requires sums of / at the corners of the control volume. The determination of / at the corners of the cell depends on the representation of the interface. In the present case with spherical particles it simply amounts to evaluating the Euclidean distance of the grid point from the surface of the particle. This may be generalized to other analytically given shapes. If the interface is not defined analytically, but only by the positions of the Lagrangian marker points, an appropriate scheme to determine / can be derived as well. The above approximation enhances computational speed at the price of reverting from exact integration to numerical integration with the same order of accuracy as the other constituents of the algorithm. Nevertheless, it yields the exact volume fraction in special cases, namely if (a) the interface does not intersect with the cell, i.e. the trivial case, (b) if the interface is plane and perpendicular to one of the coordinate directions of the Eulerian grid, (c) if the interface is plane and contains the midpoint of the cell (see Fig. 10). To assess the numerical error and the convergence behavior of the method, the volume of a sphere and an ellipsoid is computed on a Cartesian grid. The exact solution is given by Vref = 4pabc/3, with a, b and c being the semi-axes of the ellipsoid. The level-set function for spheroidal geometries is

/i;j;k

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðxi;j;k  xp Þ2 ðyi;j;k  yp Þ ðzi;j;k  zp Þ2 ¼ þ þ  1: 2 2 a c2 b

ð34Þ

For validation, the interfaces were placed in the domain X = [0; 2]  [0; 2]  [0; 2]. The diameter of the sphere was set to D = 0.8, while for the ellipsoid a = 0.4, b = 0.6, c = 0.4 was chosen. The center of both particles is xTp ¼ ð1; 1; 1Þ. As shown in Table 4, the numerical error decreases rapidly with increasing resolution of the interface D/h where h is the step size of the Cartesian grid. For D/h = 16 the relative error is about 0.5%. The number of 16 grid points along the interface is a typical value used in production runs and corresponds to approximately 800 Lagrangian marker points on the surface of a sphere. Obviously, second order convergence of the scheme is confirmed. This matches the overall second-order accuracy of the basic fluid solver. Furthermore, the center point of the interfaces was varied in additional computations (not shown here) and second order convergence was confirmed in all cases. If two particles get very close or if the surfaces are in direct contact the scheme proposed in (32) can be used without any modification. Indeed, in that case two interfaces may exist within one cell of the Eulerian grid. But as implied by (31) the computation of the volume integral requires the volume fraction api;j;k of each particle separately and does not require the total particle volume fraction ati;j;k ¼ api;j;k þ aqi;j;k . Hence, the computation of this term is still correct for two interfaces within one cell or for overlapping interfaces. 4.4. Numerical efficiency In this section the numerical efficiency of the modified scheme is investigated for a particle-laden vertical channel flow with Np = 100, 500 and 2000 particles. The computational domain X = [0; 8]  [0; 2]  [0; 1] is discretized with Nx  Ny Table 4 Relative error and convergence of the level-set method for the sample sphere and the ellipsoid. The order of convergence was computed according to Eq. (20). D/h

sphere [%]

8 16 32 64 128 256 512

1.954  100 5.123  101 1.293  101 3.178  102 8.048  103 2.091  103 4.950  104

ssphere

ellipsoid [%]

sellipsoid

1.96 2.00 2.02 1.99 1.97 2.05

2.700  100 6.667  101 1.725  101 3.983  102 1.115  102 2.796  103 8.443  104

2.05 1.97 2.02 1.99 2.02 2.04

3677

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

Table 5 CPU time of fluid solver without the new forcing scheme of Section 3.3 (f) and the time used in particle routines (p). Index 1: scheme of [42], index 2: present scheme (28). Parameters as described in the text. Particles

tf,1 [s]

tp,1 [s]

tf,2 [s]

tp,2 [s]

100 500 2000

4.826 4.830 4.943

0.067 0.130 0.404

4.881 4.894 4.960

0.085 0.194 0.589

 Nz = 512  256  128 points, corresponding to a total number of 16.8 million grid points. The particle surface is discretized with 742 forcing points and the ratio of diameter to the step size of the grid is D/h = 15.4. The simulations were performed on 32 processors of an SGI Altix. The CPU times for the original and the modified scheme are provided in Table 5 confirming the numerical efficiency of the proposed scheme. 4.5. Validation 4.5.1. Sedimenting and buoyant particles in quiescent flow To validate the new scheme the sedimentation of a spherical particle in a fluid at rest is simulated for a configuration used in [42] (Case 2 in that reference), according to experiments of Mordant and Pinton [22]. The experiment was conducted in a large container, while the simulations are conducted with periodic boundary condition with sufficient domain size so as to represent that situation fairly well. Other computations of this case where performed in [32]. The numerical setup was chosen such that the values of the particle Reynolds number

Rep ¼

u p Dp

ð35Þ

m

and the Froude number

up Fr ¼ pffiffiffiffiffiffiffiffiffiffiffi jgjDp

ð36Þ

are identical to those in the experiment. The relevant parameters are assembled in Table 6. The computational domain is X = [0; Lx]  [0; Ly]  [0; Lz] with Lx = 1.28, Ly = 10.24 and Lz = 1.28, covered with a spatially uniform mesh of Nx  Ny  Nz = 128  1048  128 Cartesian grid points. Periodic boundary conditions were applied in all three spatial directions. The fluid–solid interface is discretized with 874 equally-spaced marker points.

Table 6 Characteristics of the simulation of a sedimenting particle. Rep

Fr

(qp/qf)

g

Dp

m

Ly

h

367

1.80

2.56

9.81

0.167

1.04  103

10.24

0.01

0 1

ρp / ρf = 1.0 ρp / ρf = 0.9

-0.5

ρp / ρf = 0.5 ρp / ρf = 0.4

up / uref

up / uref

ρp / ρf = 1.1

-1

-1.5

-2

0.5

ρp / ρf = 0.3

0

0

10

t / tref

(a)

20

30

0

2

4

t / tref

6

(b)

Fig. 11. Settling and buoyant particles in fluid at rest. (a) Comparison between the method in [42] and the present scheme for Rep = 367, : experiment Mordant & Pinton [22], – – – –: method of Uhlmann [42], ---: present scheme, – – – –: present scheme without temporal derivatives (30), (b) Application of the new scheme to particles with various density ratios below qp/qf = 1.4.

3678

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684 Table 7 Setup of the sedimentation experiments given by 10 Cate et al. [35]. Case

Re

mf

qp/qf

1 2

11.6 31.9

1.17e4 6.04e5

1.164 1.166

8

0

Re = 11.6 Re = 31.9

Re = 11.6 Re = 31.9

6

yp / Dp

up [m/s]

-0.05 4

-0.1 2

0

0.5

1

t

1.5

(a)

2

0 0

0.5

1

t

1.5

2

(b)

Fig. 12. Sedimentation of particles with density ratio close to one according the experiment of ten Cate et al. [35] for Re = 11.6 and Re = 31.9. (a) Velocity, (b) trajectories of the particles.

The numerical results of the present method are shown in Fig. 11a in comparison to the experimental and the result pffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffidata ffi of [42]. As in [42], the particle velocity and the time are normalized with uref ¼ Dp jgj and tref ¼ Dp =jgj, respectively. The results show that the computed solution matches the reference data very well. Analogous simulations with Rep = 41 and Rep = 280, (Case 1 and Case 4 in [42], not displayed here) also yielded very good agreement. An additional simulation also reported in Fig. 11a illustrates the effect of neglecting the temporal derivatives (30) in (31). These terms need to be accounted for to correctly represent the unsteady particle motion. As expected, the same particle terminal velocity is reached in all cases, but the transient behavior of the particles is different. The modified scheme was also applied to particles with density ratios qp/qf  1 and to particles which are lighter than the surrounding fluid. The numerical grid and the interface discretization are the same as above. The particles have zero initial velocity and are released in fluid at rest. The results are shown in Fig. 11b, demonstrating that stable time integration is indeed obtained in a wide density range. The case qp/qf = 1 yields a stationary particle, without numerical instability, but will be addressed again with particle motion in the following section. In our experiments the modified scheme is observed to be stable down to density ratios qp/qf  0.3. Excessive grid refinement or time step reduction does not further alter the stability range of the scheme. None of the reported simulations would be possible with the original method but would lead to instability as verified in our own tests. Of particular importance here is the fact that non-buoyant particles qp  qf are now covered as these are relevant for theoretical studies. 4.5.2. Sedimentation of particles with density ratio close to one A detailed validation with such a case is now performed with the experimental data of ten Cate et al. [35], where sedimentation of a single particle with qp/qf  1.16 was investigated. This case was previously simulated by Apte et al. [1]. The computational domain X = [0; 6.67Dp]  [0; 13.33Dp]  [0; 6.67Dp] was discretized with Nx  Ny  Nz = 256  512  256 uniformly distributed grid points. A spherical particle with density qp = 1120 kg/m3 and diameter Dp = 0.015 is released in fluid at rest at a height of yp = 8.5Dp. Two different cases with particle terminal Reynolds numbers of Re = 11.6 and Re = 31.9 are considered here. The simulation parameters are assembled in Table 7. Fig. 12 shows the corresponding results with the present method compared to the experimental data.1 For Re = 31.9 the latter exhibit a shift between the data yp(t) for t < 0.32 and t > 0.38, which possibly is a calibration problem not affecting the velocity data, tough. The agreement between the present results and the experiment is excellent. 4.5.3. Rotation and translation of a sphere in linear shear flow The computation of the rotational and translational motion of a spherical particle with the proposed method is now validated by comparison with the results obtained by the pseudo-spectral method of Bagchi and Balachandar [2]. The first test is performed for the rotational motion of a sphere in linear shear flow given by

uðyÞ ¼ U 0 þ Sy

1

We acknowledge A. ten Cate for proving these data in electronic form.

ð37Þ

3679

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

y

x

u z

Fig. 13. Sketch of the configuration of sphere in linear shear flow.

Table 8 Lift coefficient cL of rotating and non-rotating sphere in linear shear flow with G = 0.2 at three different Reynolds numbers. Rep

cL,nonrotating

20 100 200

cL,rotating

Present

Bagchi [2]

Present

Bagchi [2]

0.012 0.014 0.048

0.011 0.016 0.058

0.053 0.009 0.039

0.047 0.007 0.047

Table 9 Rotation rate of a sphere held by torque-free suspension in linear shear flow with G = 0.2. ReS

xp,3/xf,3 Present

Bagchi [2]

4 20 40

0.694 0.361 0.161

0.695 0.395 0.170

as illustrated in Fig. 13. Here, the position of the particle is fixed, i.e. xp = const. Two situations are considered, the non-rotating case and the case where the particle is free to rotate around the fixed position of its center, experiencing no other torque than by the surrounding fluid. The problem hence is described by two Reynolds numbers, Rep based on the relative fluid velocity at the sphere center Uc, and

ReS ¼

SD2p

m

ð38Þ

where S is the shear rate of the flow. The two Reynolds numbers are related by ReS = G Rep, where

G ¼ SDp =U c

ð39Þ

is the shear parameter. The computational domain is X = [0; 30Dp]  [0; 15Dp]  [0; 15Dp] and is discretized with Nx  Ny  Nz = 512  256  256 equidistant points. The spatial resolution of the sphere is Dp/h  18 with NL = 916 marker points employed over the surface of the sphere. Quantitative results are summarized in Tables 8 and 9. Table 8 reports the lift coefficient

cL ¼ 1

Fy

q U 2 pR2p 2 f c

ð40Þ

for different Reynolds numbers with G = 0.2 in the case of the stationary non-rotating sphere and for the freely rotating sphere. The rotation rate xp,3, scaled with the ambient rotation rate of the fluid xf,3 = S/2, is presented in Table 9. The results of both tables show very good agreement with the reference values and confirm the suitability of the present method for rotating as well as non-rotating particles. In a second test, validation is performed by means of the acceleration rate the particle experiences with respect to translation as well as rotation. This is done for qp/qf = 1.05 and qp/qf = 5. Recall that for the lower value the simulations can not be conducted with the original scheme. For the rotational problem, the sphere is fixed until a fully developed steady flow field is established. Then, the sphere is allowed to rotate freely due to the viscous forces at the surface while its position xp is fixed. Fig. 14a illustrates increase and saturation of the angular velocity in the shear flow for different Reynolds numbers and

3680

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

1 Re = 20 / Re = 20 / Re = 200 / Re = 200 /

0.6

ρ = 1.05 ρ=5 ρ = 1.05 ρ=5

Re = 20 / Re = 20 / Re = 200 / Re = 200 /

ρ = 1.05 ρ=5 ρ = 1.05 ρ=5

0.4

0.6

u / ut

ωp / ω f

0.8

0.4

0.2 0.2 0 0.001

0 0.01

0.1

t*

1

10

0.01

0.1

1

10

100

t*

(a)

(b)

Fig. 14. Transient behavior of the particles released in steady flow, (a) rotation, (b) translation. The reference velocity ut is the terminal linear velocity, xf is the rotation rate of the fluid, and t⁄ is the non-dimensional time with reference time tref = Dp/Uc.

Table 10 Translational relaxation time tr,trans and rotational relaxation time tr,rot defined as the time to reach 63.2% of the asymptotic value.

qp/qf

Rep

1.05 5.00 1.05 5.00

20 20 200 200

t⁄ = tr,rot/tref

t⁄ = tr,trans/tref

Present

Uhlmann [42]

Bagchi [2]

Present

Uhlmann [42]

Bagchi [2]

0.76 1.85 4.70 14.49

– 1.85 – 15.52

0.78 2.10 4.30 13.8

1.81 4.05 5.88 15.4

– 4.03 – 15.45

1.73 4.17 5.58 15.0

density ratios. The relaxation time tr,rot to achieve 63.2% of the asymptotic rotation rate, scaled by tref = Dp/Uc is reported in Table 10. For the translational problem the linear acceleration of the particle in uniform flow is considered, i.e. ReS = 0. Again, the sphere is fixed until steady state is reached. Then, the sphere is released and free translational motion is allowed and the particle is accelerated as shown in Fig. 14b. The relaxation time tr,trans/tref is also reported in Table 10. Both relaxation times are in excellent agreement with the reference data of [42]. The agreement with the experiment is not as good but still acceptable. For the lower Reynolds number and the larger density the deviation is 12% for the rotative acceleration while only 3% for both linear cases. These simulations complete a series of test cases, validating the accurate prediction of the particle terminal velocity as well as the transient behavior characterized by translational and rotational acceleration of the particles. The accurate computation of these phenomena is an essential prerequisite not only for modeling the individual motion of single particles but most importantly for the numerical modeling of particle–particle and particle-wall collisions.

5. Issues related to interface contact and collision modeling 5.1. Problems with the basic scheme In a particle-laden flow, collisions can occur between two particles and between a particle and a wall. In many studies of particle laden flows, a simple repulsive potential is used for collision modeling to prevent particles from overlapping as the one proposed by Glowinski et al. [11]

Fcol n;pq ¼

xp  xq



ðmaxf0; ðfn  SÞgÞ2 ;

ð41Þ

where Fcol n;pq is the repulsive force between particles p and q, fn the normal distance of the two surfaces,  a problem specific model constant and S the range of the repulsive force. The regularized Dirac functions used in Eqs. (7b) and (7d) generate difficulties in contact and collision modeling due to overlapping stencils at low particle distances. Therefore, the value of S characterizing the distance at which the potential sets in is generally chosen equal to twice the step size of the Eulerian grid, so that the surfaces are practically never in direct contact. Eq. (41) constitutes about the simplest possible collision model and bears some arbitrariness in the choice of its parameters. The model was successfully employed for the simulation of dilute suspensions, for example by [44] in the case of a particle-laden flow in a vertical channel with  ¼ 8  104 Dp =ðqf u21 Þ. In other flow configurations, it may be advantageous or even required to use different prefactors [10,25,1].

3681

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

0.02

up [m/s]

0.015

ζn [m]

0.4

CFL = 0.6 CFL = 0.3 CFL = 0.1

0.01

CFL = 0.6 CFL = 0.3 CFL = 0.1

0.2 0

-0.2

0.005

-0.4 0

-0.05

0

t [s]

0.05

0.1

0

-0.05

(a)

t [s]

0.05

0.1

(b)

Fig. 15. Collision of a sedimenting particle with a stationary particle using regularized Dirac functions without modifications. The collision model (41) is used with S ¼ 0. (a) Particle trajectory, (b) particle velocity.

More realistic collision models have to account for tangential forces being transferred during particle contact and lubrication forces generated at instants when surface distances get very small. Modeling these, with part of the phenomena being resolved by the geometric resolution of the particle and other contributions being unresolved, is conceptually sound only when surface distances are allowed to reach values close to or equal to zero. In order to study the behavior of the basic IBM in this situation, the collision of a particle impacting vertically from above onto another particle fixed in space is now simulated with the characteristic length S of the repulsive potential in (41) equal to zero and  = 104. The Stokes number based on the impact velocity is St = 152 and the particle Reynolds number is Rep = 165. The computational domain X = [0; Lx]  [0; Ly]  [0; Lz] with Lx = Ly = Lz = 13.33 Dp was discretized with Nx  Ny  Nz = 256  256  256 points. The spatial resolution of the sphere is Dp/h  20. The surface of the sphere is represented by NL = 1159 marker points. As soon as the surface distance is smaller than 2h the support of the regularized delta functions used for interpolation intersects, as illustrated in Fig. 16. This violates the property of the sum of the shifted regularized delta functions at a given point of the Eulerian grid being equal to one. In mathematical terms the property of the shifted delta functions constituting a partition of unity is lost. As a consequence, an error is introduced whenever the stencils of the particles overlap. This takes place in each time step where the condition is met and hence accumulate over time. When reducing Dt the overall error therefore increases as the number of flawed time steps becomes larger. In Fig. 15 the computed solution for the normal distance and velocity of the upper particle is displayed. Obviously, the rebound trajectory is strongly dependent on the time step chosen. Dramatic is the observation that the error increases when the time step is reduced. In other words, the scheme is inconsistent, which is just intolerable. A further elucidation of the problem by means of the velocity field outside and inside the particle is provided by Fig. 18 discussed in the following section. 5.2. Modified scheme for interface contact If lubrication forces during approach and rebound and tangential collision forces are to be modeled in a conceptually clear way, a remedy for the problem described above has to be provided. Numerous attempts were made during the present study. It was attempted to provide different weights to the stencil such that these add up to values less or equal to 1, even in case of overlapping, for example. Another line of thought was to determine a weight according to the volume fraction of the solid in each cell, etc. None of these strategies yielded satisfactory results, though.

(b) (a)

Fig. 16. Lagrangian markers on the particle surfaces and points in the Eulerian grid where IBM forces are are imposed to couple particles an fluid. (a) Regular calculation of force at Lagrangian surface markers, (b) surface markers surrounding circles show examples for markers removed from the force calculation since they share identical grid points for the regularized delta function.

3682

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

0.02

up [m/s]

0.015

ζn [m]

0.4

CFL = 0.6 CFL = 0.3 CFL = 0.1 CFL = 0.05

0.01

CFL = 0.6 CFL = 0.3 CFL = 0.1 CFL = 0.05

0.2 0

-0.2

0.005

-0.4 0 -0.05

0

0.05

0.1

-0.05

0

0.05

t [s]

t [s]

(a)

(b)

0.1

Fig. 17. Dependency of trajectory on the size of the time step for particle–particle collision. Forcing points with overlapping stencils are removed during interface contact as described in the text. Additional forcing loops are performed to impose the correct velocity at the interface. All other parameters, in particular the collision model, are the same as used in Fig. 15.

Fig. 18. Vertical velocity component v for a vertical particle–particle collision at the end of the phase of direct surface contact. St = 152 and Rep = 165. The definitions for cases a to d are given in Table 11.

The solution finally found is to entirely exclude those surface markers from the forcing scheme whose stencils overlap with the stencil of markers located on the surface of a collision partner. This is illustrated in Fig. 16. With the staggered Eulerian grid employed (Fig. 1a) a Lagrangian point is removed if any of the three grids leads to overlap. In Fig. 16, the two markers in circles would be removed, for example.

3683

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684 Table 11 Parameters for the cases displayed in Fig. 18. Case

Markers removed

CFL

Dt[s]

vmax[m/s]

a b c d

no yes no yes

0.6 0.6 0.1 0.1

1.8e4 1.7e4 1.0e5 1.7e5

1.014 1.012 1.622 0.820

It is important to note that the particle trajectory still can be computed correctly, since the repulsive force during a collision is generated by the collision model which is based on the distance between the surfaces of the particles colliding or the distance of the particle from the wall. The simulations reported in Fig. 17 indeed show that the particle trajectory is independent of the time step. Fig. 18 provides an impression of the secondary flow generated inside the particle with quantitative data reported in Table 11. It is obvious that the magnitude of this flow now is reduced with the time step, while it increases when reducing the time step if the forcing points are not removed. Also note that the discretization is fairly coarse in these simulations. It is to be stressed that the flow inside the particle is not per se a problematic issue, but is accepted with the method if still the boundary conditions on the particle surface are fulfilled as they should. Problematic is the increase with reducing Dt. 6. Conclusions In the present work the immersed boundary method of Uhlmann [42] is enhanced by several modifications. The coupling of solid and fluid phase is improved by an additional forcing scheme. It does not require the solution of any additional elliptic equation so that the computational effort is practically not increased. The additional forcing allows large time steps while the no-slip condition at the interface is still imposed very accurately. The stability range of the method was significantly increased by introducing direct integration of the linear and angular momentum of the fluid inside the particle control volume employing a numerically efficient level set approach. With this modification, no assumption on the motion of the fluid inside the particle as used in [42] is required. That feature is of particular importance for subsequent collision modeling for the following reason. Suppose a model with tangential force is used so that when a particle reaches another surface substantial rotational acceleration can be generated. Angular acceleration impacts on the flow inside the particle which then substantially deviates from rigid body motion since forcing points are only located on the surface of the particle, not inside. As a consequence, the assumption of rigid body motion necessarily has to be replaced by explicit integration. An alternative option would be to distribute Lagrangian points for the interaction between particle and fluid not only along the surface but over the entire volume of the particle. This was considered by Uhlmann [45] but no substantial changes were found in the situations considered, i.e. with repulsive potential and S ¼ 2h. Further information on this option, like accuracy, cost, etc. are currently not available. The present approach is potentially less costly as it requires a substantially smaller number of markers. Furthermore, forcing the entire volume of the particle also can not solve the inconsistency issue for S ¼ 0 in the collision model since the basic problem related to intersecting stencils persists. If interfaces approach closely or are in direct contact, the basic IBM yields an inconsistent numerical scheme since the trajectory of the particles depends dramatically on the time step. We obtained a consistent scheme by excluding all surface marker points from the computation of forces at the particle surface whose stencils overlap with the stencil of a collision partner. The current volume integration employs the analytic definition of the particle geometry. If the geometry is defined locally, e.g. by a triangulated surface, modern algorithms developed for computer graphics can be employed to efficiently determine the level set function, i.e. the distance of an Eulerian grid point from the surface so that the proposed approach can easily be generalized to more complicated shapes and other ways to define the surface. The present paper, concerned with discretization issues for the IBM, lays the ground for advancing the modeling of collisions in the framework of the present type of Euler-Langrange scheme. Such a development, focusing on physical modeling rather than discretization issues has been accomplished and is presented in [14]. Acknowledgments This work was partially funded by DFG trough grant Fr1593/5-1. M. Uhlmann and C. Chan-Braun are acknowledged for stimulating discussions on IBM and particle-laden flows. The computations were performed on the SGI Altix at ZIH Dresden.2 References [1] S.V. Apte, M. Martin, N.A. Patankar, A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows, J. Comput. Phys. 228 (2009) 2712–2738. 2

http://tu-dresden. de/die_tu_dresden/zentrale_einrichtungen/zih/hpc/hochleistungsrechner/altix.

3684

T. Kempe, J. Fröhlich / Journal of Computational Physics 231 (2012) 3663–3684

[2] P. Bagchi, S. Balachandar, Effect of free rotation on the motion of a solid sphere in linear shear flow at moderate Re, Phys. Fluids 14 (2002) 2719–2737. [3] W.P. Breugem, M. Poriquie, B.J. Boersma, Some issues related to the stress IB method for flow around obstacles. in: Proceedings of the Academy colloquium: Immersed Bounday Methods: Current status and future research direction, Amsterdam, The Netherlands, 15–17 June 2009. [4] Y. Cheny, O. Botella, The LS-STAG method: a new immersed boundary/level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties, J. Comput. Phys. 229 (2010) 1043–1076. [5] E. Chow, A.J. Cleary, R.D. Falgout, Design of the hypre preconditioner library, in: SIAM Workshop on Object Oriented Methods for Inter-operable Scientific and Engineering Computing, SIAM, 1998, pp. 106–116. [6] R.D. Falgout, U.M. Yang, hypre: a library of high performance preconditioners, in: International Conference on Computational Science, vol. 3, 2002, pp. 632–641. [7] J.H. Ferzinger, M. Peric, Computational Methods for Fluid Dynamics, Springer-Verlag, 2002. [8] M. Francois, E. Uzgoren, J. Jackson, W. Shyy, Multigrid computations with the immersed boundary technique for multiphase flows, Int. J. Numer. Methods Heat Fluid Flow 14 (2004) 98–115. [9] A. Gilmanov, F. Sotiropoulos, A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies, J. Comput. Phys. 207 (2005) 457–492. [10] R. Glowinski, T.W. Pan, T.I. Hesla, D.D. Joseph, J. Priaux, 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 (2001) 363–426. [11] R. Glowinski, T.-W. Pan, T.I. Hesla, D.D. Joseph, A distributed Lagrange multiplier/fictitious domain method for particulate flows, Int. J. Multiphase Flow 25 (1999) 755–794. [12] F.H. Harlow, J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (1965) 2182– 2189. [13] G. Iaccarino, R. Verzicco, Immersed boundary technique for turbulent flow simulations, Appl. Mech. Rev. 56 (2003) 331–347. [14] T. Kempe, J. Fröhlich, Collision modeling for the interface-resolved simulation of spherical particles in viscous fluids, submitted for publication. [15] D. Kim, H. Choi, Immersed boundary method for flow around an arbitrarily moving body, J. Comput. Phys. 212 (2006) 662–680. [16] D.V. Le, B.C. Khoo, K.M. Lim, An implicit-forcing immersed boundary method for simulating viscous flows in irregular domains, Comput. Methods Appl. Mech. Eng. 197 (25–28) (2008) 2119–2130. [17] P. Leopardi, A partition of the unit sphere into regions of equal area and small diameter, Electron. Trans. Numer. Anal. 25 (2006) 309–327. [18] Z. Li, K. Ito, The immersed interface method: numerical solutions of PDEs involving interfaces and irregular domains, Soc. Indus. Math. (2006). [19] R. Mittal, H. Dong, M. Bozkurttas, F.M. Najjar, A. Vargas, A. von Loebbecke, A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries, J. Comput. Phys. 227 (2008) 4825–4852. [20] R. Mittal, G. Iaccarino, Immersed boundary methods, Ann. Rev. Fluid Mech. 37 (2005) 239–261. [21] J. Mohd-Yusof, Combined immersed boundary/B-Spline method for simulations of flows in complex geometries, Center for Turbulence Research, Annual Research Briefs, NASA Ames/ Stanford University, pp. 317–327, 1997. [22] N. Mordant, J.-F. Pinton, Velocity measurement of a settling sphere, Eur. Phys. J. B 18 (2000) 343–352. [23] H. Nirschl, H.A. Dwyer, V. Denk, Three-dimensional calculations of the simple shear flow around a single particle between two moving walls, J. Fluid Mech. 283 (1995) 273–285. [24] M. Oevermann, C. Scharfenberg, R. Klein, A sharp interface finite volume method for elliptic equations on Cartesian grids, J. Comput. Phys. 228 (2009) 5184–5206. [25] T.-W. Pan, R. Glowinski, Direct simulation of the motion of neutrally buoyant circular cylinders in plane Poiseuille flow, J. Comput. Phys. 181 (2002) 260–279. [26] Ch. S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25 (1977) 220–252. [27] Ch. S. Peskin, The immersed boundary method, Acta Numer. 11 (2003) 479–517. [28] A. Prosperetti, G. Tryggvason, Computational Methods for Multiphase Flow, Cambridge University Press, 2007. [29] A.M. Roma, Ch. S. Peskin, M.J. Berger, An adaptive version of the immersed boundary method, J. Comput. Phys. 153 (1999) 509–534. [30] F. Roman, V. Armenio, J. Fröhlich, A simple wall-layer model for large eddy simulation with immersed boundary method, Phys. Fluids 21 (2009) 101701. [31] F. Roman, E. Napoli, B. Milici, V. Armenio, An improved immersed boundary method for curvilinear grids, Comput. Fluids 38 (2009) 1510. [32] N. Sharma, N.A. Patankar, A fast computation technique for the direct numerical simulation of rigid particulate flows, J. Comput. Phys. 205 (2005) 439– 457. [33] P.R. Spalart, R.D. Moser, M.M. Rogers, Spectral methods for the Navier–Stokes equations with one infinite and two periodic directions, J. Comput. Phys. 96 (1991) 297–324. [34] T. Stoesser, J. Fröhlich, W. Rodi, Turbulent open channel flow over a permeable bed, in: 32th IAHR Congress, Venice, Italy, July 1–6, 2007. [35] A. ten Cate, C.H. Nieuwstad, J.J. Derksen, H.E.A. Van den Akker, Particle imaging velocimetry experiments and lattice-Boltzmann simulations on a single sphere settling under gravity, Phys. Fluids 14 (2002) 4012–4025. [36] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, Y.-J. Jan, A front-tracking method for the computations of multiphase flow, J. Comput. Phys. 169 (2001) 708–759. [37] G. Tryggvason, S. Thomas, J. Lu, B. Aboulhasanzadeh, V. Tsengue, Multiscale computation of multiphase flows, in: 7th International Conference on Multiphase Flows, Tampa, Florida, USA, 2010. [38] H.S. Udaykumar, H.-C. Kan, W. Shyy, R. Tran-Son-Tay, Multiphase dynamics in arbitrary geometries on fixed Cartesian grids, J. Comput. Phys. 137 (1997) 366–405. [39] H.S. Udaykumar, R. Mittal, P. Rampunggoon, A. Khanna, A sharp interface Cartesian grid method for simulating flows with complex moving boundaries, J. Comput. Phys. 174 (2001) 345–380. [40] M. Uhlmann. First experiments with the simulation of particulate flows, Technical Report No. 1020, CIEMAT, Madrid, Spain, 2003, ISSN pp. 1135–9420. [41] M. Uhlmann. New results on the simulation of particulate flows, Technical Report No. 1038, CIEMAT, Madrid, Spain, 2003, ISSN, pp. 1135-9420. [42] M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, J. Comput. Phys. 209 (2005) 448–476. [43] M. Uhlmann, An improved fluid-solid coupling method for DNS of particulate flow on a fixed mesh, in M. Sommerfeld, (Ed.,) Proc. 11th Workshop TwoPhase Flow Predictions, Merseburg, Germany, 2005, Universität Halle. [44] M. Uhlmann, Interface-resolved direct numerical simulation of vertical particulate channel flow in the turbulent regime, Phys. Fluids 20 (2008) 053305. [45] M. Uhlmann, Private communication, 2010. [46] S.O. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, J. Comput. Phys. 100 (1992) 25–37. [47] E. Uzgoren, R. Singh, J. Sim, W. Shyy, Computational modeling for multiphase flows with spacecraft application, Prog. Aerosp. Sci. 43 (2007) 138–192. [48] M. Sint-Annaland van, N.G. Deen, J.A.M. Kuipers, Numerical simulation of gas bubbles behaviour using a three-dimensional volume of fluid method, Chem. Eng. Sci. 60 (2005) 2999–3011. [49] R. Verzicco, J. Mohd-Yusof, P. Orlandi, D. Haworth, LES in complex geometries using boundary body forces, AIAA J. 38 (2000) 427–433.