Derivation and validation of a novel implicit second-order accurate immersed boundary method

https://doi.org/10.1016/j.jcp.2008.03.031Get rights and content

Abstract

A novel implicit second-order accurate immersed boundary method (IBM) for simulating the flow around arbitrary stationary bodies is developed, implemented and validated in this paper.

The IBM is used to efficiently take into account the existence of bodies within the fluid domain. The flow domain consist of simple Cartesian cells whereas the body can be arbitrary. At the triangulated interface of the body and the fluid, the immersed boundary, the coefficients obtained from discretizing the Navier–Stokes equations are closed with a second-order accurate interpolation arising from the immersed boundary condition employed at the interface. Two different conditions are developed in this paper and it is shown that for the mirroring method the resulting coefficients lead to a well-posed and diagonally dominant system which can be efficiently solved with a preconditioned Krylov sub-space solver. The immersed boundary condition generates a fictitious reversed velocity field inside the immersed boundary, which is excluded from the continuity equation to account for the presence of the IB in the pressure correction equation, resulting in no mass flux over the IB. The force acting on the object from the fluid is determined by integrating the pressure and the viscous forces over the object.

The method is validated by simulating the flow around a sphere for a range of Re numbers. It is shown that the drag is very well in agreement with experimental data. Accuracy and convergence studies are employed, proving the second-order accuracy of the method and showing the superiority in convergence rate compared to other IBM. Finally the drag force of a cluster of non-spherical particles is employed to show the generality and potential of the method.

Introduction

Although dispersed multiphase flows are common in both nature and industry, their fundamental behavior is still poorly understood. In recent years, multiphase computational fluid dynamics (CFD) has become an important tool to model and gain more insight into such flows. However, to successfully predict and simulate multiphase flows at a large scale, the fundamental behavior at small scale needs to be understood. The majority of multiphase flows are fluid–solid flows, in which a dispersed solid phase is present in a continuous fluid phase.

Almost all multiphase CFD today is directly aimed at resolving large scales. One of the key issues in successfully describing fluid–solid flows is the physics of the interaction of the two phases. This interaction is a complex process and, to successfully capture it, all the important length and time scales must be resolved; hence true direct numerical simulations (DNS) are required. In true DNS, the full Navier–Stokes equations are solved, taking into account the particle surface, thereby resolving the boundary layer and wake of each particle. In this article, we present a second-order accurate and relatively inexpensive true DNS method to accurately resolve the flow around individual particles.

One of the first methods developed to accurately simulate fluid–solid flows is the arbitrary Lagrangian–Eulerian (ALE) method, proposed by Hu [11], [12]. In this method, an unstructured grid is created around the particles and adapted as the particles move. This requires remeshing each time step. Although the method shows good results, the remeshing is computationally very expensive and makes it difficult to provide accurate results. To decrease the computational cost involved by remeshing and to potentially increase the accuracy, three methods without the necessity of remeshing have been employed in the literature: Cartesian, Lagrange Multiplier and Immersed Boundary methods.

In the Cartesian method, or the cell-splitting method [26], [33], [13], [17], the domain mainly consists of Cartesian grid cells. Where a Cartesian cell is cut by the surface segment of a particle, the fluid cells are split and locally unstructured cells are formed. The flow is resolved for the new grid with irregular boundary conditions. The drawback is that the irregular boundary conditions are hard to generalize and small cells can lead to problems with the conservation and stability in the solver.

The Lagrange multiplier method [7], [8], [6] is a finite-element method where the coupling between the boundary of the particle and the fluid is introduced by a Lagrange multiplier in the integral formulation of the Navier–Stokes equations. The weak coupling is implicitly implemented and, by minimizing the combined integral formulation, the solution of the flow field is found. The method is second-order accurate in space. The drawback with the method is that it is relatively expensive and only preserves the mass and momentum globally, not locally.

The immersed boundary (IB) method resolves these difficulties by using a Cartesian grid in the whole domain. Where the particle surface segments (or immersed boundaries) cross the Cartesian grid cells, the flow variables are directly modified, without the necessity of rearranging the grid. This modification of the flow variables, representing the coupling between the two phases, can be implemented in several different ways. One common way is to exert a Lagrangian force at the immersed boundary onto the fluid phase. The Navier–Stokes equations are then solved on the Cartesian grid, including a number of point forces, representing the influence of the immersed boundaries on the fluid flow.

There exist in general two methods for calculating the forces to represent the immersed boundaries. In the first approach, the forces are calculated from the momentum equations and distributed on the Cartesian mesh by employing a distribution function, first proposed by Peskin [29]; this is called the distributive IB method. The second approach is a non-distributive method, where the Eulerian force is calculated from the momentum equations directly or by employing a boundary condition exactly at the immersed boundary.

The distributive IB method [29], [18] employs a function, based on the discrete Dirac delta function, to distribute the Lagrangian forces on the Cartesian mesh. As a result, the Lagrangian forces are smeared out over several grid cells, leading to a smeared representation of the immersed boundary and to only first-order accuracy. In this method, the Lagrangian force is determined explicitly from the generalized Hooke’s law.

Goldstein et al. [9] also employs the distributive IB method, but calculates the Lagrangian force by constraining the fluid to the no-slip boundary condition on the immersed boundary. The force proposed by Goldstein et al. [9] contains two unphysical parameters which are dependent on the time scales of the flow. These constants tend to become relatively large, which gives stiff equations to solve. Due to this the method only works well for small time steps. Silva and co-workers [21], [22], [27] improved the calculation of the Lagrangian forces by making the forces originate from a physical point of view, generally increasing the robustness. The method is explicit and first-order accurate.

To obtain second-order accuracy, Mohd-Yusof [24], [25] proposed a non-distributive method, the momentum forcing method. In this method, the no-slip boundary condition is enforced directly at the IB. This is done by introducing a smooth velocity gradient over the sharp IB. To determine the smooth gradient, the velocity field is mirrored over the IB. An explicit momentum force is applied to the interior boundary velocities, such that a linear or bilinear interpolation of the velocity exactly fulfills the no-slip boundary condition at the IB. Due to the reversed velocity field, problems with the mass conservation arise in the boundary cell. The method is employed by, among others, Fadlun et al. [3] and Kang et al. [14].

To deal with the mass conservation errors obtained with the Mohd-Yusof non-distributive method, Kim and Choi [16], [15] introduces a mass source/sink in the cells surrounding the IB. This new explicit method shows promising and robust results.

Gilmanov and co-workers [5], [4] have also proposed a second-order hybrid Cartesian and IB method, where the IB is triangulated. The external boundary velocities are set by a Dirichlet boundary condition such that an interpolation of the velocity along the normal of the IB fulfills the no-slip boundary condition. The Neumann pressure boundary condition is set in a similar way. The boundary condition for the pressure is applied explicitly, and the boundary condition for the velocity is discretized with first-order accuracy, directly on the exterior grid points. The method shows promising results for low Reynolds numbers.

Majumdar et al. [23] have extended Mohd-Yusof’s momentum forcing method to a boundary condition at the IB. To constrain the velocity of the fluid at the IB, a boundary condition is employed exactly at the IB. This boundary condition is formulated either explicitly or implicitly, is second-order accurate in space, and employs three different interpolation schemes. The method has been implemented for stationary IBs along with RANS models to describe turbulent flows. The method shows promising results, but has potential problems with the weighting coefficients in the boundary condition which may result in wiggles in the resulting solution. Moreover, the method generates an unphysical mass flux over the IB.

In [32], Tseng and Ferzinger extended Mohd Yusof’s and Fadluns ideas to a ghost-cell immersed boundary method, which extrapolates the fluid flow on a point inside the IB such that the velocity is constrained at the IB, with second-order accuracy. The method shows good results for complex geometries, but an unphysical mass flux over the IB is generated by the method.

The purpose of this paper is to construct, validate and compare two novel stable second-order implicit IB methods. The arbitrary IB is triangulated and in the first method a subset of the triangle vertices are used as control points in which the constraints are applied. A trilinear interpolation is applied to ensure the correct velocity at the IB. Mohd-Yusof’s momentum force [25] is rewritten to a fully implicit immersed boundary condition and the generated internal fictitious velocity field is excluded from the continuity equation to account for the presence of the IB in the pressure correction equation, resulting in no mass flux over the IB.

The second method mirrors the velocity along the normal of the triangulated IB segment such that it becomes exactly defined at the IB. Trilinear interpolation is employed to interpolate the velocity to the mirrored velocity point outside the IB.

The present implementation of the methods is done with a staggered variable configuration using a segregated finite volume solver based upon the pressure correction method SIMPLEC [2]. The methods are validated and compared by determining the drag force over a sphere for several different Reynolds numbers. Additionally, the drag force of a cluster of non-spherical particles is simulated to show the generality and potential of the methods.

In the paper, the novel IBM is proven to be second-order accurate in space and to have a superior convergence rate compared to earlier IBMs. The convergence rate is increased because there exists no mass flux over the IB.

Section snippets

Governing equations

The equations governing the flow around the immersed boundaries are the continuity and momentum equations for incompressible flow, the Navier–Stokes equations,ujxj=0,ρt(ui)+ρujuixj=-pxi+xjμuixj+fi,where ρ is the density of the flow, ui is the velocity field of the flow, p is the pressure of the flow, μ represents the viscosity of the fluid, and fi represents any external source terms. This set of equations is solved together with the immersed boundary condition,ui=uiibon the IB,

The immersed boundary condition

The methods developed and employed in this paper are three-dimensional and second-order accurate immersed boundary methods. The methods constrain the velocity of the fluid to a specified velocity exactly at the surface of any body immersed in the fluid flow. The velocity is constrained by an implicitly formulated immersed boundary condition (IBC). In this work, two different discretization schemes of the IBC have been developed. As a result of the application of the IBC, a fictitious velocity

Calculation of the surface forces

In most IB methods, the fluid velocity and pressure at the IB are controlled by distributive forces. In such cases, the total force on the fluid as a result of the IB is simply the sum of the distributive forces applied to each IB element. Since in the current method the flow properties are controlled in the discretization itself, such a force is not explicitly calculated, but often this force is required. The total force is given by the surface integral of the stress tensor over the IB,Fi=IB(-

Numerical simulation of a sphere immersed in a fluid

To validate the derived IB methods, the flow around a unit sphere at various Reynolds numbers is simulated. For the flow over a sphere at low Reynolds numbers Re1.0, the Stokes drag force is given byfd=fvisc+fpress=4πμRv+2πμRv,where R is the radius of the sphere, μ the viscosity of the fluid and v is the free stream velocity. From this, the Stokes drag coefficient can be derived,Cd=24Re.At higher Reynolds numbers, semi-empirical equations are used to validate the method, Table 1 [30].

The

Conclusions

In this paper, a novel implicit second-order accurate immersed boundary method for simulating the flow around stationary arbitrary bodies is developed, implemented and validated. The velocity of the fluid is constrained by an implicit immersed boundary condition such that it exactly follows the immersed boundary. Two such conditions are developed, validated and compared, the vertex-constraining (VCIB) method which constrains the velocity at vertex control points and the mirroring (MIB) method

Acknowledgment

This work was financed by the Swedish Research Council (VR), Grant No. 621-2004-2008.

References (33)

Cited by (0)

View full text