- Email: [email protected]

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 ﬂows 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 ﬂow Particles Immersed boundary method Euler–Lagrange coupling Collision modeling

a b s t r a c t An efﬁcient approach for the simulation of ﬁnite-size particles with interface resolution was presented by Uhlmann [M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate ﬂows, 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 ﬂuid density larger than some critical value above 1, hence excluding, for example, non-buoyant particles. This can be remedied by an efﬁcient integration step for the artiﬁcial ﬂow ﬁeld 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 ﬂows in which solid particles are transported. These range from mechanical and chemical process engineering to energy conversion, biological ﬂows, environmental ﬂows, etc. The accurate and efﬁcient numerical simulation of this type of ﬂow hence is of substantial importance for fundamental research as well as for understanding and optimizing complex ﬂow problems. The three-dimensional, fully coupled and interface-resolving simulation of ﬂows with a huge number of particles of ﬁnite size to date is an area of active research and development. An efﬁcient 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 efﬁcient Cartesian grid for the discretization of the ﬂuid phase and to represent the immersed ﬂuid–solid interface which does not coincide with any grid line by modiﬁcations 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 difﬁcult and tedious for complex conﬁgurations. Moreover, certain geometries such as closely packed spheres, for example, just can not be handled with body-ﬁtted grids, except if the discretization is fully unstructured, while they pose no problem for an IBM [34]. For immobile ﬂuid–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 ﬂow 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 ﬂow 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 ﬁxed 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 ﬂow [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 ﬂuid onto the solid may be generated if the scheme is not sufﬁciently 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 ﬂows, 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 ﬂuctuations 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 ﬂuid, to employ local body-ﬁtted 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 ﬂuid forces onto moving particles in a relatively easy way, since regularization of the forces coupling ﬂuid 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 ﬂow ﬁeld as a constant-density ﬁeld, 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 ﬂow. 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 ﬁrst 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 modiﬁcations 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 ﬂuid 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 efﬁcient 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 ﬂuid 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 ﬂuids 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 ﬂuid 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 ﬂuid, qf the density of the ﬂuid, and t the time. The spatial discretization of (1)–(3) is performed by a second-order ﬁnite-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 ﬂuid solver alone has been throughly validated and convergence of second order in space and time has been checked. Coupling of the ﬂuid 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 ﬂuid–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 deﬁned by the three-point version of Roma et al. [29]

8 qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ > 1 > 5 3krk 3ð1 krkÞ2 þ 1 ; 0:5 6 krk 6 1:5 > >6 < qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 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 efﬁciency 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 ﬁeld 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 ﬁeld 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 ﬁeld uk of that Runge–Kutta sub-step (7g). The related pressure ﬁeld is given by Eq. (7h). The Runge–Kutta coefﬁcients 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 efﬁcient so that simulations with very large numbers of particles are possible as, e.g., performed in [44]. The deﬁnition of surface markers and not volume markers all over the particle volume allows for ﬂow 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 fulﬁlled [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 ﬂuid. Rather, a velocity error e = jUd UCj remains, the magnitude of which depends on the conﬁguration. 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 ﬂow around a ﬁxed 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 inﬁnite ﬂow, since this yields a handy test case with low cost perfectly fulﬁlling 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 ﬂow 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 ﬂow over a circular cylinder with different para-meters. (a) Re = 40 and various time steps. (b) Impact of grid reﬁnement 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 conﬁguration is simulated with the grid reﬁned by a factor of two in each spatial direction using Grid 2 instead of Grid 1. The grid and the time step are reﬁned 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 ﬁne grid, regardless of the smaller time step. Obviously, the error can not be reduced by spatial grid reﬁnement 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 ﬁne 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 ﬂuid 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 efﬁcient 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 ﬁrst for the preliminary velocity ﬁeld without volume force (7a), then after the forcing step and the solution of the Helmholtz Equation (7e), and ﬁnally 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 modiﬁed due to the projection of the velocity onto a divergence-free ﬁeld. 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 ﬁeld 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 ﬁeld 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 modiﬁed scheme is displayed in Fig. 8a using the same conﬁguration 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 conﬁgurations, 2 to 3 additional forcing loops are sufﬁcient 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 conﬁguration. The modiﬁed 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 ﬁrst placed in steady ﬂow 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 ﬂow over a circular cylinder with Re = 10 and CFL = 1.0. : preliminary velocity ﬁeld 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 modiﬁed forcing scheme with (nf = 3) and without (nf = 0) additional forcing loops. (a) Velocity error at the interface for the ﬂow over stationary cylinder at Re = 10. (b) Velocity error at the interface for the ﬂow 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 conﬁrms the present results in Section 3.1. (2) The spatial interpolation introduces an error of order Dx. This is not conﬁrmed 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 conﬁrmed 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 ﬁeld 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 ﬁnal method. 3.4. Convergence A reﬁnement study is now performed employing the reference scheme and employing the new scheme to assess its convergence behavior. The setup for the ﬂow around a ﬁxed 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 ﬁne grid with Nx Ny = 512 512 points, so that the following deﬁnitions of the L1-error and L2-error will be used:

n o ui;j uref i;j vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2ﬃ 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 reﬁned grids using the expression [7]

s¼

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 ﬂuid 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 reﬁnement study for the ﬂow around a stationary cylinder without additional forcing loops. The errors 1 and 2 are deﬁned 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 reﬁnement study for the ﬂow 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 ﬂuid 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 ﬂow can be described by the trajectory of the center of mass xp(t) and the orientation, speciﬁed 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 sufﬁcient 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 ﬂuid–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 efﬁcient. The momentum balance of the ﬂuid enclosed in an arbitrary volume X with surface C moving with the ﬂow 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 ﬂuid 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 modiﬁed 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 ﬂuid 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 ﬂuid hence can not be computed with the basic method. 4.2. Modiﬁed 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 ﬂuid–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, deﬁned 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 ﬂuid in the k-th Runge–Kutta sub-step then is given by the ﬁnite-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 modiﬁed 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 identiﬁed. Subsequently, the volume fraction is determined. For arbitrarily aligned three-dimensional cells, 64 cut-cell types are found which can be reduced to ﬁve generic types [48,24]. One can also employ a ﬁner background mesh for this step, together with a stair-step approximation of ﬁrst order. In [1], a reﬁnement 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–ﬂuid 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 deﬁned 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

sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 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 conﬁrmed. This matches the overall second-order accuracy of the basic ﬂuid solver. Furthermore, the center point of the interfaces was varied in additional computations (not shown here) and second order convergence was conﬁrmed 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 modiﬁcation. 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 efﬁciency In this section the numerical efﬁciency of the modiﬁed scheme is investigated for a particle-laden vertical channel ﬂow 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 ﬂuid 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 modiﬁed scheme are provided in Table 5 conﬁrming the numerical efﬁciency of the proposed scheme. 4.5. Validation 4.5.1. Sedimenting and buoyant particles in quiescent ﬂow To validate the new scheme the sedimentation of a spherical particle in a ﬂuid at rest is simulated for a conﬁguration 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 sufﬁcient 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 ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 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 ﬂuid–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 ﬂuid 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 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃdata ﬃ 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 modiﬁed scheme was also applied to particles with density ratios qp/qf 1 and to particles which are lighter than the surrounding ﬂuid. The numerical grid and the interface discretization are the same as above. The particles have zero initial velocity and are released in ﬂuid 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 modiﬁed scheme is observed to be stable down to density ratios qp/qf 0.3. Excessive grid reﬁnement 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 veriﬁed 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 ﬂuid 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 ﬂow 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 ﬁrst test is performed for the rotational motion of a sphere in linear shear ﬂow 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 conﬁguration of sphere in linear shear ﬂow.

Table 8 Lift coefﬁcient cL of rotating and non-rotating sphere in linear shear ﬂow 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 ﬂow 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 ﬁxed, i.e. xp = const. Two situations are considered, the non-rotating case and the case where the particle is free to rotate around the ﬁxed position of its center, experiencing no other torque than by the surrounding ﬂuid. The problem hence is described by two Reynolds numbers, Rep based on the relative ﬂuid velocity at the sphere center Uc, and

ReS ¼

SD2p

m

ð38Þ

where S is the shear rate of the ﬂow. 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 coefﬁcient

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 ﬂuid xf,3 = S/2, is presented in Table 9. The results of both tables show very good agreement with the reference values and conﬁrm 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 ﬁxed until a fully developed steady ﬂow ﬁeld is established. Then, the sphere is allowed to rotate freely due to the viscous forces at the surface while its position xp is ﬁxed. Fig. 14a illustrates increase and saturation of the angular velocity in the shear ﬂow 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 ﬂow, (a) rotation, (b) translation. The reference velocity ut is the terminal linear velocity, xf is the rotation rate of the ﬂuid, 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 deﬁned 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 ﬂow is considered, i.e. ReS = 0. Again, the sphere is ﬁxed 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 ﬂow, collisions can occur between two particles and between a particle and a wall. In many studies of particle laden ﬂows, 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 speciﬁc model constant and S the range of the repulsive force. The regularized Dirac functions used in Eqs. (7b) and (7d) generate difﬁculties 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 ﬂow in a vertical channel with ¼ 8 104 Dp =ðqf u21 Þ. In other ﬂow conﬁgurations, 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 modiﬁcations. 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 ﬁxed 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 ﬂawed 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 ﬁeld outside and inside the particle is provided by Fig. 18 discussed in the following section. 5.2. Modiﬁed 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 ﬂuid. (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 deﬁnitions for cases a to d are given in Table 11.

The solution ﬁnally 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 ﬂow generated inside the particle with quantitative data reported in Table 11. It is obvious that the magnitude of this ﬂow 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 ﬂow 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 fulﬁlled 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 modiﬁcations. The coupling of solid and ﬂuid 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 signiﬁcantly increased by introducing direct integration of the linear and angular momentum of the ﬂuid inside the particle control volume employing a numerically efﬁcient level set approach. With this modiﬁcation, no assumption on the motion of the ﬂuid 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 ﬂow 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 ﬂuid 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 deﬁnition of the particle geometry. If the geometry is deﬁned locally, e.g. by a triangulated surface, modern algorithms developed for computer graphics can be employed to efﬁciently 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 deﬁne 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 ﬂows. 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-ﬂow interactions in complex ﬂows, 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 ﬂow 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 ﬂow 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 ﬂows 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 Scientiﬁc 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 ﬂows, Int. J. Numer. Methods Heat Fluid Flow 14 (2004) 98–115. [9] A. Gilmanov, F. Sotiropoulos, A hybrid Cartesian/immersed boundary method for simulating ﬂows 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 ﬁctitious domain approach to the direct numerical simulation of incompressible viscous ﬂow past moving rigid bodies: application to particulate ﬂow, J. Comput. Phys. 169 (2001) 363–426. [11] R. Glowinski, T.-W. Pan, T.I. Hesla, D.D. Joseph, A distributed Lagrange multiplier/ﬁctitious domain method for particulate ﬂows, Int. J. Multiphase Flow 25 (1999) 755–794. [12] F.H. Harlow, J.E. Welch, Numerical calculation of time-dependent viscous incompressible ﬂow of ﬂuid with free surface, Phys. Fluids 8 (1965) 2182– 2189. [13] G. Iaccarino, R. Verzicco, Immersed boundary technique for turbulent ﬂow 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 ﬂuids, submitted for publication. [15] D. Kim, H. Choi, Immersed boundary method for ﬂow 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 ﬂows 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 ﬂows 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 ﬂows 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 ﬂow 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 ﬁnite 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 ﬂow, J. Comput. Phys. 181 (2002) 260–279. [26] Ch. S. Peskin, Numerical analysis of blood ﬂow 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 ﬂows, 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 inﬁnite and two periodic directions, J. Comput. Phys. 96 (1991) 297–324. [34] T. Stoesser, J. Fröhlich, W. Rodi, Turbulent open channel ﬂow 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 ﬂow, J. Comput. Phys. 169 (2001) 708–759. [37] G. Tryggvason, S. Thomas, J. Lu, B. Aboulhasanzadeh, V. Tsengue, Multiscale computation of multiphase ﬂows, 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 ﬁxed 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 ﬂows with complex moving boundaries, J. Comput. Phys. 174 (2001) 345–380. [40] M. Uhlmann. First experiments with the simulation of particulate ﬂows, Technical Report No. 1020, CIEMAT, Madrid, Spain, 2003, ISSN pp. 1135–9420. [41] M. Uhlmann. New results on the simulation of particulate ﬂows, 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 ﬂows, J. Comput. Phys. 209 (2005) 448–476. [43] M. Uhlmann, An improved ﬂuid-solid coupling method for DNS of particulate ﬂow on a ﬁxed 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 ﬂow 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-ﬂuid ﬂows, J. Comput. Phys. 100 (1992) 25–37. [47] E. Uzgoren, R. Singh, J. Sim, W. Shyy, Computational modeling for multiphase ﬂows 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 ﬂuid 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.