# A Comparative Study on Evaluation Methods of Fluid Forces on Cartesian Grids

Authors: Taku Nonomura; Junya Onishi
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

We investigate the accuracy and the computational efficiency of the numerical schemes for evaluating fluid forces in Cartesian grid systems. A comparison is made between two different types of schemes, namely, polygon-based methods and mesh-based methods, which differ in the discretization of the surface of the object. The present assessment is intended to investigate the effects of the Reynolds number, the object motion, and the complexity of the object surface. The results show that the mesh-based methods work as well as the polygon-based methods, even if the object surface is discretized in a staircase manner. In addition, the results also show that the accuracy of the mesh-based methods is strongly dependent on the evaluation of shear stresses, and thus they must be evaluated by using a reliable method, such as the ghost-cell or ghost-fluid method.


---

## Body

## 1. Introduction

In recent years, Cartesian grids have been widely used in computational fluid dynamics applications. This is not only because of the simplicity of data structure and algorithms, but also because of the fast, automatic, and robust grid generation for complex geometries. The latter reason is particularly important for industrial applications, since grid generation is one of the most time-consuming tasks. Cartesian grids have potential to overcome this difficulty.However, flow simulation using Cartesian grids still has some issues to be addressed before being put into practice. A key issue is the development of accurate methods for evaluating physical quantities and fluxes at the surface of solid objects. This issue comes from the fact that Cartesan grids do not necessarily align with the surface of the objects. As a result, there arises some arbitrariness in the location and the orientation of the object surfaces, and thus the evaluation of physical quantities and fluxes at the object surface is highly dependent on the numerical schemes used. To improve the numerical evaluation described above, many efforts have been already devoted since the pioneering work by Peskin [1]. An extensive review of the methods developed in the literatures is given by Mittal and Iaccarino [2].Despite these efforts, however, there has been no established method that can be applied to any problems. This is partly due to the lack of the systematic investigation on the accuracy of the numerical schemes used in these methods. In this study, we assess and compare the accuracy and the computational efficiency of the numerical schemes for evaluating physical quantities and fluxes in Cartesian grid systems. In particular, we focus on the evaluation of the local fluid stresses at the surface of an object and on the evaluation of the total fluid forces acting on the object body, as the integration of the local stresses.The fluid force acting on a solid body is evaluated by integrating the fluid stresses over the body surface:(1)Fj=∫ΓSσijnidS,where σij is the stress tensor and nj is the unit outward normal of the surface element. The stress tensor σij of Newtonian fluids is given by(2)σij=-pδij+μ∂uj∂xi+∂ui∂xj,where p is the pressure, ui is the velocity, and μ is the viscosity.In numerical computations, the surface integral defined by (1) is evaluated in two steps:(i)
discretization of the surface,(ii)
estimation of the physical quantities at the surface.The numerical procedures to conduct these steps depend on the grid system used in the simulation. In body-conforming curvilinear grid systems, shown in Figure 1(a), both steps can be implemented on the specified grids in a straightforward manner. However, in nonconforming grid systems, shown in Figure 1(b), some special treatment is required since the grids do not exactly match the body surface. To the best of our knowledge, there are two different ways which have been adopted in nonconforming grid systems, namely, the control volume method and the direct calculation.Figure 1
Examples of body-conforming and nonconforming grids.
(a)
Body-conforming grid
(b)
Nonconforming gridThe control volume method introduces a virtual volume which encompasses the body and estimates the fluid forces from the net momentum flux through the surface of this virtual volume. This method was firstly applied to steady forces by Lai and Peskin [3], and later, it was extended to unsteady forces by Balaras [4] and to moving body problems by Shen et al. [5]. The control volume method has an advantage that it is unnecessary to directly handle the complex geometries expressed by the immersed boundaries. On the other hand, it has a severe disadvantage that local stresses cannot be computed; that is, only the total force acting on the body can be evaluated. Moreover, the contributions from pressure and viscous forces cannot be evaluated separately.In the direct calculation, the body surface is usually discretized by a set of polygons, and the local stresses on these surface elements are estimated. Then, the local stresses are integrated over the body surfaces to obtain the total force, as in body-conforming grid systems. However, unlike the cases of body-conforming grid systems, the estimation of the flow variables on the surface elements is not straightforward and requires elaborate interpolation and/or extrapolation methods in order to accurately take into account the location (and orientation if necessary) of the surface. Moreover, the computational cost of the stress integration may become prohibitive when the body shape is complex, since it depends on the number of the polygons used to express the body surface. Such a dependency is undesirable, in particular, in engineering applications. For example, the number of the polygons reaches up to the order of107 for an aerodynamic simulation using a detailed car model with both exterior and interior parts [6].In this study, we discuss numerical schemes that are suitable for evaluating fluid forces in nonconforming grid systems. In particular, we propose a simple and computationally efficient method for calculating the surface integral in (1). The proposed method is a kind of direct calculation described above. Therefore, it has inherited an advantage over the control volume method in that the contributions from pressure and viscous forces can be evaluated locally and separately. Moreover, in the present method, the fluid stresses are integrated over the mesh faces (staircase surfaces) in the grid system used, not over the discretized surface elements of the original body. Therefore, the undesirable dependency on quality or numbers of surface polygons can be removed. Hereafter, the present method is referred to as a “mesh-based scheme (MS),” whereas the method based on the integration over surface polygons is called a “polygon-based scheme (PS).” We present a comparative study of the accuracy of these methods.It is worth noting here that, in both mesh-based and polygon-based schemes, the physical quantities at the discretized surfaces must be estimated by interpolation or extrapolation, and the accuracy of this estimation would influence the final results of the fluid forces. For the PS, we have employed the ghost-cell approach described by Mittal et al. [7]. For the MSs, we have tested three approaches that are often used in the studies of immersed boundary methods: a staircase approach, a ghost-cell approach, and a ghost-fluid approach, the details of which are described in Section 2.The rest of the paper is organized as follows. Section2 describes the details of the computational procedures of a polygon-based scheme and three mesh-based schemes. The computational grid, the discretization of the body surface, and the extrapolation method used in these schemes are fully discussed. In Section 3, both polygon-based and mesh-based schemes are applied to evaluate the drag force acting on a sphere in Stokes flow. The obtained results are compared with the analytical solution. The grid convergence test is also performed to investigate the numerical accuracy of the schemes. Stokes flow is used here because analytical solutions are available not only for the fluid forces but also for the flow fields, that is, pressure and velocity fields. By using the analytical solutions of the flow fields, the error associated with the force evaluation schemes can be assessed separately from the error in the computation of the flow fields themselves. In contrast, in higher Reynolds number flows, the flow fields must be numerically computed. Therefore, these errors are hardly evaluated separately. In Section 4, the extension to the moving body problem and the computational efficiency of the force evaluation schemes are discussed. In addition, to show the applicability of the present mesh-based scheme in higher Reynolds number flows, the results for the case of the Reynolds number of 100 are presented in Section 5. In Section 6, the discussion on the applicability to the complex body problem is conducted. Finally, the last section summarizes this study.


## 2. Numerical Methods

Herein, Section2.1 describes the computational grid and defines the cell types. The polygon-based force evaluation scheme, which seems to be commonly used in the immersed boundary method, and the mesh-based force evaluation scheme proposed in this paper are described in Sections 2.2 and 2.3, respectively.


### 2.1. Computational Grids and Definition of Cell Type

The grid system used in this work is illustrated in Figure2. In this figure, the entire domain is divided into equally sized cubic cells. These cells are classified into two groups according to their locations:(i)
fluid cells (cells whose center is in the fluid region),(ii)
solid cells (cells whose center is in the solid region).Within these categories, special cells close to the interface are defined as follows: (i)
boundary cells (fluid cells adjacent to at least one solid cell).Figure 2
Classification of the computational cells.The flow properties used to evaluate (1), namely, the pressure and the velocity, are assumed to be available only at the center of the fluid cells, including the boundary cells. When the values at the solid cells are needed, they must be interpolated or extrapolated by some methods considering a body shape. In the following, we call such solid cells ghost cells.Note that, although we use collocated layout of variables throughout this study, the extension to staggered layout of variables is trivial. The accuracy of mesh-face integration, which is key idea of this paper, does not seem to be affected by the layout of variables. Note also that the above cell classification is not trivial, especially when the solid region has a complicated shape. This issue, however, is beyond the scope of the present study and is not further discussed here.


### 2.2. Polygon-Based Scheme

The polygon-based scheme (PS) appears to be most commonly implemented in the immersed boundary method [8]. This scheme divides three-dimensional body surfaces into sets of polygons or two-dimensional surfaces into series of lines. Hereafter, to avoid confusion, these surface elements will be referred to as “facets.” A schematic of this approach is shown in Figure 3(a).Figure 3
Polygon and mesh-based methods.
(a)
Polygon-based method
(b)
Mesh-based methodIn PS, the surface integral defined by (1) is conveniently evaluated by summing over the facets:(3)Fj=∑facetσijniΔSfacet,where ΔS is the surface area of the facet, and the normal stresses are summed over all facets. The remaining task is to estimate the bracketed quantities at the center of each facet, which requires the geometric properties (the surface normal and the area of the facet) and the physical properties (the pressure and the velocity). The former quantities can be computed from the coordinates of the polygon vertexes. In this paper, the latter are estimated by an approach using the normal, which is described in the following.As shown in Figure3(a), a normal probe is extended from the center of the facet (indicated by I), and two markers (indicated by I1 and I2) are placed at distances d1 and d2, respectively, along the normal probe. These markers are used as reference points for interpolating the value at the center of the facet. The pressure field assumes a linear distribution along the probe:(4)pd=C1d+C0.The coefficients C1 and C0 are determined such that(5)pd1=p1,pd2=p2,where p1 and p2 are the values at the reference points I1 and I2, respectively. From this linear extrapolation, the value at the point I is obtained as(6)pI=p0=p1-d1d2-d1p2-p1.The velocity components are computed from a quadratic function. For example, the x-component of the velocity is fitted to(7)ud=C2d2+C1d+C0,where the coefficients are determined by the values at both reference points and the center of the facet:(8)ud1=u1,ud2=u2,u0=uB,where u1 and u2 are the velocities at the reference points I1 and I2, respectively, and uB is the velocity at the center of the facet. Note that uB is the velocity of the object because no-slip and no-penetration conditions are imposed. Differentiating this obtained function yields the components of the velocity gradient tensor at the point I(9)∂u∂nI=u1d22-u2d12d1d2d2-d1,∂u∂sI=0,∂u∂tI=0,where the coordinates n, s, and t are defined along the normal and the tangential directions at the facet, respectively. The other components are computed similarly.In the above extrapolation, the values at the reference points should be estimated in advance. In this paper, the reference values are computed by trilinear interpolation from the quantities in the surrounding fluid cells. If the reference points are not carefully chosen, the interpolation involves a quantity in a solid cell that must also be extrapolated. To avoid such recursive interpolations and extrapolations, the distances to the reference points are simply selected to be long enough. In this paper,d1 and d2 are set to 1.75Δx and 3.5Δx, respectively. Note that the value 1.75Δx comes from the possible largest distance between neighboring cell centers, namely, 3Δx for a three-dimensional case, and the value 3.5Δx is twice as long as it. This method is called a simplified ghost-cell method in this paper.


### 2.3. Mesh-Based Scheme

The PS described above is simple in algorithm, but it seems to be prohibitive for practical applications where complex geometries must be handled. This is because the PS requires operations ofONS, where NS is the number of the polygons used to represent the geometries. For example, NS is of the order of 107 for a detailed car model with both exterior and interior parts [6]. This problem becomes significant, in particular, on distributed memory platform, since the operations are intensive only in the small area of the computational domain and thus the work load is not balanced well.As an alternative to the PS, we propose mesh-based schemes (MSs), in which the body is simplified by mesh faces between boundary and solid cells, and the surface integral defined by (1) is discretized over the mesh faces. A schematic illustration of this approach is shown in Figure 3(b).In MSs, the surface integral defined by (1) is decomposed into its directional components, each approximated from the mesh face in its corresponding direction:(10)Fx=∫σxxnxdS+∫σyxnydS+∫σzxnzdS≈∑face,xσxxnxΔSface,x+∑face,yσyxnyΔSface,y+∑face,zσzxnzΔSface,z.Here, since the directional components can be independently summed, the surface normal and the area are easily evaluated from the geometrical characteristics of the mesh face. For example, nx=±1 and ΔS=ΔyΔz are held for the mesh faces in the x-direction. The surface normal and the areas of mesh faces in other orientations can be similarly evaluated. Note also that the summations involve only the faces of the boundary cells, which simplifies the computational implementation.The stress at the center point of the mesh face (Ix and Iy in Figure 3(b)) can be computed in several ways. In the following, three methods, which seem to be used widely in the studies of immersed boundary methods, are discussed.


#### 2.3.1. MS1: A Staircase Approach

In MS1, the pressure is simply extrapolated using the cell-centered value:(11)pI,x=pb.c.,pI,y=pb.c.,where pI,x and pI,y are, respectively, the pressures at Ix and Iy and pb.c. is the pressure at the related boundary cell. To approximate the shear stress, the velocities at Ix and Iy are assumed to correspond to that of the object uI=0:(12)∂u∂xI,x=ui,j-uIΔx/2,∂u∂yI,x=ui,j+1-ui,j-12Δx,where the computational stencil used for these derivatives is depicted in Figure 4.Figure 4
Stencils for stress evaluation and extrapolation to the solid cell.
(a)
Stencil for computing the derivatives
(b)
Interpolation to solid cellsAlthough the above staircase approach may be used less often nowadays, it is introduced in order to show its inaccuracy and to clearly demonstrate the important aspects of the other schemes introduced later.


#### 2.3.2. MS2: A Ghost-Cell Approach

In MS2, the assumed pressure atIx is the average of the pressures at the centers of the adjacent boundary cell and the solid cell.(13)pI,x=pi,j+pi-1,j2.The shear stress is approximated by a second-order central difference scheme. The x- and y-derivatives of the x component of the velocity at the point Ix are estimated as follows:(14)∂u∂xI,x=ui,j-ui-1,jΔx,∂u∂yI,x=12ui,j+1-ui,j-12Δx+ui-1,j+1-ui-1,j-12Δx.The computational stencil used to compute these derivatives is also depicted in Figure 4(a).Note that the computation of bothx- and y-derivatives requires quantities inside the solid region, in other words, in the ghost cells. For MS2, these quantities are estimated by the simplified version (avoiding recursive interpolation and extrapolation) of the ghost-cell method [7]. Concretely, the quantities in solid cells are computed by the fitting functions (see (6) and (7)) constructed in the normal probe approach illustrated in Figure 4(b):(15)pS=p-d0,uS=u-d0,where d0 is the distance between the solid cell and the closest point on the body surface.


#### 2.3.3. MS3: A Ghost-Fluid Approach

In MS3, the ghost-fluid method described by Gibou et al. [9] is applied to estimate the quantities in the ghost cells.For example, in order to estimate a quantity at the solid celli-1,j in Figure 4(a), a quadratic function along the x-axis is introduced and is fitted to the two points at the fluid cells i,j and i+1,j and the interface point:(16)fx=C2x2+C1x+C0,where f denotes the pressure or a component of the flow velocity, x is the coordinate of the x-axis, and Cii=0,1,2 is the fitting coefficients to be determined. Note that the interface point is defined as the intersection between the solid surface and the line segment connecting i,j and i-1,j. Here, the origin is set to i,j, and the distance from the origin to the interface point is assumed to be θΔx, without loss of generality. For the x-component of the flow velocity, three coefficients in (16) are uniquely determined by solving the following system of equations:(17)u0=ui,j,uΔx=ui+1,j,u-θΔx=uB,where uB is the prescribed boundary value for u. As a result, the target value at the solid cell is estimated as follows:(18)ui-1,j=u-Δx.The values for the other components can be calculated in a similar manner. The pressure can be also estimated in a similar manner, but it is necessary to replace one of the conditions by the Neumann boundary condition:(19)p0=pi,j,pΔx=pi+1,j,∂p∂x-θΔx=∂p∂xB,where ∂p/∂xB is the x-component of the pressure gradient at the boundary, which is assume to be 0 in this study.


## 2.1. Computational Grids and Definition of Cell Type

The grid system used in this work is illustrated in Figure2. In this figure, the entire domain is divided into equally sized cubic cells. These cells are classified into two groups according to their locations:(i)
fluid cells (cells whose center is in the fluid region),(ii)
solid cells (cells whose center is in the solid region).Within these categories, special cells close to the interface are defined as follows: (i)
boundary cells (fluid cells adjacent to at least one solid cell).Figure 2
Classification of the computational cells.The flow properties used to evaluate (1), namely, the pressure and the velocity, are assumed to be available only at the center of the fluid cells, including the boundary cells. When the values at the solid cells are needed, they must be interpolated or extrapolated by some methods considering a body shape. In the following, we call such solid cells ghost cells.Note that, although we use collocated layout of variables throughout this study, the extension to staggered layout of variables is trivial. The accuracy of mesh-face integration, which is key idea of this paper, does not seem to be affected by the layout of variables. Note also that the above cell classification is not trivial, especially when the solid region has a complicated shape. This issue, however, is beyond the scope of the present study and is not further discussed here.


## 2.2. Polygon-Based Scheme

The polygon-based scheme (PS) appears to be most commonly implemented in the immersed boundary method [8]. This scheme divides three-dimensional body surfaces into sets of polygons or two-dimensional surfaces into series of lines. Hereafter, to avoid confusion, these surface elements will be referred to as “facets.” A schematic of this approach is shown in Figure 3(a).Figure 3
Polygon and mesh-based methods.
(a)
Polygon-based method
(b)
Mesh-based methodIn PS, the surface integral defined by (1) is conveniently evaluated by summing over the facets:(3)Fj=∑facetσijniΔSfacet,where ΔS is the surface area of the facet, and the normal stresses are summed over all facets. The remaining task is to estimate the bracketed quantities at the center of each facet, which requires the geometric properties (the surface normal and the area of the facet) and the physical properties (the pressure and the velocity). The former quantities can be computed from the coordinates of the polygon vertexes. In this paper, the latter are estimated by an approach using the normal, which is described in the following.As shown in Figure3(a), a normal probe is extended from the center of the facet (indicated by I), and two markers (indicated by I1 and I2) are placed at distances d1 and d2, respectively, along the normal probe. These markers are used as reference points for interpolating the value at the center of the facet. The pressure field assumes a linear distribution along the probe:(4)pd=C1d+C0.The coefficients C1 and C0 are determined such that(5)pd1=p1,pd2=p2,where p1 and p2 are the values at the reference points I1 and I2, respectively. From this linear extrapolation, the value at the point I is obtained as(6)pI=p0=p1-d1d2-d1p2-p1.The velocity components are computed from a quadratic function. For example, the x-component of the velocity is fitted to(7)ud=C2d2+C1d+C0,where the coefficients are determined by the values at both reference points and the center of the facet:(8)ud1=u1,ud2=u2,u0=uB,where u1 and u2 are the velocities at the reference points I1 and I2, respectively, and uB is the velocity at the center of the facet. Note that uB is the velocity of the object because no-slip and no-penetration conditions are imposed. Differentiating this obtained function yields the components of the velocity gradient tensor at the point I(9)∂u∂nI=u1d22-u2d12d1d2d2-d1,∂u∂sI=0,∂u∂tI=0,where the coordinates n, s, and t are defined along the normal and the tangential directions at the facet, respectively. The other components are computed similarly.In the above extrapolation, the values at the reference points should be estimated in advance. In this paper, the reference values are computed by trilinear interpolation from the quantities in the surrounding fluid cells. If the reference points are not carefully chosen, the interpolation involves a quantity in a solid cell that must also be extrapolated. To avoid such recursive interpolations and extrapolations, the distances to the reference points are simply selected to be long enough. In this paper,d1 and d2 are set to 1.75Δx and 3.5Δx, respectively. Note that the value 1.75Δx comes from the possible largest distance between neighboring cell centers, namely, 3Δx for a three-dimensional case, and the value 3.5Δx is twice as long as it. This method is called a simplified ghost-cell method in this paper.


## 2.3. Mesh-Based Scheme

The PS described above is simple in algorithm, but it seems to be prohibitive for practical applications where complex geometries must be handled. This is because the PS requires operations ofONS, where NS is the number of the polygons used to represent the geometries. For example, NS is of the order of 107 for a detailed car model with both exterior and interior parts [6]. This problem becomes significant, in particular, on distributed memory platform, since the operations are intensive only in the small area of the computational domain and thus the work load is not balanced well.As an alternative to the PS, we propose mesh-based schemes (MSs), in which the body is simplified by mesh faces between boundary and solid cells, and the surface integral defined by (1) is discretized over the mesh faces. A schematic illustration of this approach is shown in Figure 3(b).In MSs, the surface integral defined by (1) is decomposed into its directional components, each approximated from the mesh face in its corresponding direction:(10)Fx=∫σxxnxdS+∫σyxnydS+∫σzxnzdS≈∑face,xσxxnxΔSface,x+∑face,yσyxnyΔSface,y+∑face,zσzxnzΔSface,z.Here, since the directional components can be independently summed, the surface normal and the area are easily evaluated from the geometrical characteristics of the mesh face. For example, nx=±1 and ΔS=ΔyΔz are held for the mesh faces in the x-direction. The surface normal and the areas of mesh faces in other orientations can be similarly evaluated. Note also that the summations involve only the faces of the boundary cells, which simplifies the computational implementation.The stress at the center point of the mesh face (Ix and Iy in Figure 3(b)) can be computed in several ways. In the following, three methods, which seem to be used widely in the studies of immersed boundary methods, are discussed.


### 2.3.1. MS1: A Staircase Approach

In MS1, the pressure is simply extrapolated using the cell-centered value:(11)pI,x=pb.c.,pI,y=pb.c.,where pI,x and pI,y are, respectively, the pressures at Ix and Iy and pb.c. is the pressure at the related boundary cell. To approximate the shear stress, the velocities at Ix and Iy are assumed to correspond to that of the object uI=0:(12)∂u∂xI,x=ui,j-uIΔx/2,∂u∂yI,x=ui,j+1-ui,j-12Δx,where the computational stencil used for these derivatives is depicted in Figure 4.Figure 4
Stencils for stress evaluation and extrapolation to the solid cell.
(a)
Stencil for computing the derivatives
(b)
Interpolation to solid cellsAlthough the above staircase approach may be used less often nowadays, it is introduced in order to show its inaccuracy and to clearly demonstrate the important aspects of the other schemes introduced later.


### 2.3.2. MS2: A Ghost-Cell Approach

In MS2, the assumed pressure atIx is the average of the pressures at the centers of the adjacent boundary cell and the solid cell.(13)pI,x=pi,j+pi-1,j2.The shear stress is approximated by a second-order central difference scheme. The x- and y-derivatives of the x component of the velocity at the point Ix are estimated as follows:(14)∂u∂xI,x=ui,j-ui-1,jΔx,∂u∂yI,x=12ui,j+1-ui,j-12Δx+ui-1,j+1-ui-1,j-12Δx.The computational stencil used to compute these derivatives is also depicted in Figure 4(a).Note that the computation of bothx- and y-derivatives requires quantities inside the solid region, in other words, in the ghost cells. For MS2, these quantities are estimated by the simplified version (avoiding recursive interpolation and extrapolation) of the ghost-cell method [7]. Concretely, the quantities in solid cells are computed by the fitting functions (see (6) and (7)) constructed in the normal probe approach illustrated in Figure 4(b):(15)pS=p-d0,uS=u-d0,where d0 is the distance between the solid cell and the closest point on the body surface.


### 2.3.3. MS3: A Ghost-Fluid Approach

In MS3, the ghost-fluid method described by Gibou et al. [9] is applied to estimate the quantities in the ghost cells.For example, in order to estimate a quantity at the solid celli-1,j in Figure 4(a), a quadratic function along the x-axis is introduced and is fitted to the two points at the fluid cells i,j and i+1,j and the interface point:(16)fx=C2x2+C1x+C0,where f denotes the pressure or a component of the flow velocity, x is the coordinate of the x-axis, and Cii=0,1,2 is the fitting coefficients to be determined. Note that the interface point is defined as the intersection between the solid surface and the line segment connecting i,j and i-1,j. Here, the origin is set to i,j, and the distance from the origin to the interface point is assumed to be θΔx, without loss of generality. For the x-component of the flow velocity, three coefficients in (16) are uniquely determined by solving the following system of equations:(17)u0=ui,j,uΔx=ui+1,j,u-θΔx=uB,where uB is the prescribed boundary value for u. As a result, the target value at the solid cell is estimated as follows:(18)ui-1,j=u-Δx.The values for the other components can be calculated in a similar manner. The pressure can be also estimated in a similar manner, but it is necessary to replace one of the conditions by the Neumann boundary condition:(19)p0=pi,j,pΔx=pi+1,j,∂p∂x-θΔx=∂p∂xB,where ∂p/∂xB is the x-component of the pressure gradient at the boundary, which is assume to be 0 in this study.


## 2.3.1. MS1: A Staircase Approach

In MS1, the pressure is simply extrapolated using the cell-centered value:(11)pI,x=pb.c.,pI,y=pb.c.,where pI,x and pI,y are, respectively, the pressures at Ix and Iy and pb.c. is the pressure at the related boundary cell. To approximate the shear stress, the velocities at Ix and Iy are assumed to correspond to that of the object uI=0:(12)∂u∂xI,x=ui,j-uIΔx/2,∂u∂yI,x=ui,j+1-ui,j-12Δx,where the computational stencil used for these derivatives is depicted in Figure 4.Figure 4
Stencils for stress evaluation and extrapolation to the solid cell.
(a)
Stencil for computing the derivatives
(b)
Interpolation to solid cellsAlthough the above staircase approach may be used less often nowadays, it is introduced in order to show its inaccuracy and to clearly demonstrate the important aspects of the other schemes introduced later.


## 2.3.2. MS2: A Ghost-Cell Approach

In MS2, the assumed pressure atIx is the average of the pressures at the centers of the adjacent boundary cell and the solid cell.(13)pI,x=pi,j+pi-1,j2.The shear stress is approximated by a second-order central difference scheme. The x- and y-derivatives of the x component of the velocity at the point Ix are estimated as follows:(14)∂u∂xI,x=ui,j-ui-1,jΔx,∂u∂yI,x=12ui,j+1-ui,j-12Δx+ui-1,j+1-ui-1,j-12Δx.The computational stencil used to compute these derivatives is also depicted in Figure 4(a).Note that the computation of bothx- and y-derivatives requires quantities inside the solid region, in other words, in the ghost cells. For MS2, these quantities are estimated by the simplified version (avoiding recursive interpolation and extrapolation) of the ghost-cell method [7]. Concretely, the quantities in solid cells are computed by the fitting functions (see (6) and (7)) constructed in the normal probe approach illustrated in Figure 4(b):(15)pS=p-d0,uS=u-d0,where d0 is the distance between the solid cell and the closest point on the body surface.


## 2.3.3. MS3: A Ghost-Fluid Approach

In MS3, the ghost-fluid method described by Gibou et al. [9] is applied to estimate the quantities in the ghost cells.For example, in order to estimate a quantity at the solid celli-1,j in Figure 4(a), a quadratic function along the x-axis is introduced and is fitted to the two points at the fluid cells i,j and i+1,j and the interface point:(16)fx=C2x2+C1x+C0,where f denotes the pressure or a component of the flow velocity, x is the coordinate of the x-axis, and Cii=0,1,2 is the fitting coefficients to be determined. Note that the interface point is defined as the intersection between the solid surface and the line segment connecting i,j and i-1,j. Here, the origin is set to i,j, and the distance from the origin to the interface point is assumed to be θΔx, without loss of generality. For the x-component of the flow velocity, three coefficients in (16) are uniquely determined by solving the following system of equations:(17)u0=ui,j,uΔx=ui+1,j,u-θΔx=uB,where uB is the prescribed boundary value for u. As a result, the target value at the solid cell is estimated as follows:(18)ui-1,j=u-Δx.The values for the other components can be calculated in a similar manner. The pressure can be also estimated in a similar manner, but it is necessary to replace one of the conditions by the Neumann boundary condition:(19)p0=pi,j,pΔx=pi+1,j,∂p∂x-θΔx=∂p∂xB,where ∂p/∂xB is the x-component of the pressure gradient at the boundary, which is assume to be 0 in this study.


## 3. Force Estimation for Stokes Flows

In this section, the accuracies of PS, MS1, MS2, and MS3 are assessed by comparison with analytical solutions of Stokes flow past a sphere.


### 3.1. Stokes Flow Past a Sphere

Analytical solutions to the pressure and velocity of Stokes flow around a sphere are given by [10]:(20)p=p∞-32μρu∞axr3,u=u∞-14u∞ar3+a2r2-34u∞ax2r31-a2r2,v=-34u∞axyr31-a2r2,w=-34u∞axzr31-a2r2,where a is the radius of the sphere, ρ is the mass density, μ is the viscosity, u∞ is the freestream velocity, p∞ is the reference pressure, x, y, and z are the Cartesian coordinates, and r is the distance from the origin (where the origin is the center of the sphere).The fluid is assumed to flow parallel to thex-axis; therefore, the drag force acting on the sphere is the x-component of the force defined in (1). Substituting the analytical solutions (20) into this definition, the pressure drag and the frictional drag are, respectively, obtained as (21)CDp=Fxp1/2ρu∞2A=2πaμu∞ρu∞2/2πa2=8μρu∞D=8Re,CDf=Fxv1/2ρu∞2A=4πaμu∞ρu∞2/2πa2=16μρu∞D=16Re,where CDp and CDf are the coefficients of the pressure drag and frictional drag, respectively, and Re=ρu∞D/μ is the Reynolds number. D=2a and A=πa2 define the diameter and the projected area of the sphere, respectively.


### 3.2. Computational Settings and Data

The drag force is now computed by the numerical methods described in Section2. To compute the drag forces, a sphere of the diameter D is centralized within a cubic domain of the side length L. The value of L is arbitrary but must be sufficiently large to encompass the sphere. The present study selects L=2D. The computational cells are generated by dividing the cubic domain into 2N equally sized cells in each direction. The fineness of the sphere is altered by setting the number of grid points N occupied by the sphere in each direction to N=4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048. The solid cells constituting the sphere of grid points less than 128 are illustrated in Figure 5. The PS requires an additional polygon mesh. The force computed on each Cartesian mesh properly converges with that on a spherical equally distributed polygon mesh with 3,200×3,200 points. Therefore, this polygon mesh is adopted in the remainder of the study. The analytical solutions (20) are assumed to be available only at the centers of the fluid cells, including the boundary cells.Figure 5
Spheres constructed on meshes with different resolutions. From left to right, a sphere of diameterD is constructed on mesh resolved to Δx=a/4, x=a/8,a/16,a/32,a/64, and a/128.


### 3.3. Results

#### 3.3.1. Total Drag Forces

The forces estimated by PS, MS1, MS2, and MS3 at different mesh resolutions are plotted in Figure6. Here, the Reynolds number is Re=0.01 for all the cases. Correspondingly, the pressure and friction drags are, respectively, solved as CDp,exact=800 and CDp,exact=1600, using the analytical solutions. The result shows that the convergence rate for the pressure drag is first order in MS1 and second order in PS and MS2 and that the convergence rate for the friction drag is zeroth order in MS1 and first to second order in PS, MS2, and MS3. These results indicate that the proposed MS2 and MS3 can evaluate the fluid force as accurately as PS, but at much smaller computational cost and using a much simpler computational code. Thus, fluid forces can be simply and efficiently evaluated by MS2 and MS3 in the immersed boundary method. Comparing MS2 and MS3, it interestingly shows that MS2 has fewer errors than MS3. This clarifies that the MS2 based on the ghost-cell method is better for the force estimation than MS1 based on the ghost-fluid method. Apart from the success of MS2 and MS3, MS1 is shown to be unacceptable for practical use since the friction drag is not converged even when a very fine mesh is used. At Δx=0.06, sudden decrease in friction drag of PS is observed and this might be caused by canceling out two different error sources which is not clarified well in this study.Figure 6
Estimated drag in Stokes flow around a sphere as functions of mesh spacing. Analytical pressure and friction drag areCDp,exact=800 and CDp,exact=1600, respectively.
(a)
Pressure drag
(b)
Friction dragIt is worth noting that the PS results are sensitive to the density of the polygon mesh at coarser polygon mesh resolutions, whereas we use a very fine polygon mesh in this study to eliminate the error introduced by coarse polygon mesh. This implies that the polygon density requires special attention in the PS method.


#### 3.3.2. Force Distributions

Then, we visualized the force distributions on the each cell projected onx, y, and z plane. The partial drag coefficients for each projected cell are defined as follows:(22)Cdp,x=∫cellpdydz1/2ρu∞2ΔyΔzCdf,y=∫cellτxy+τyzdxdz1/2ρu∞2ΔxΔz⋮.Here, integration is conducted for the surface in each projected cell. In the case of a sphere, there are two projections: front and rear projections. In this definition, the reference solution can be obtained with very fine discretization. On the other hand, MSs evaluate them in simplified ways as in the inside summation of (10). The distributions of reference solutions and MSs are shown in Figures 7 and 8. Here, the reference solution is obtained with the ten times finer submesh. The results of MS1 have strong oscillation in friction drag because this method does not compute shear stress correctly. On the other hand, MS2 and MS3 can predict them quantitatively, though they have weak oscillation. This behavior is illustrated in Figures 9 and 10, which shows the partial drag on the line indicated in Figures 7 and 8. These figures show the same trend more quantitatively. Figure 8 shows that the MS2 predict partial friction drag more accurately than MS3, as discussed above for total friction drag.Figure 7
Distribution of partial pressure drag for each cell front-projected onx plane. Contour range is set to be -0.1<Cdp,x/CDp,exact<0.8. Analytical pressure drag is CDp,exact=800.
(a)
Reference
(b)
MS1
(c)
MS2
(d)
MS3Figure 8
Distribution of partial friction drag for each cell front-projected ony plane. Contour range is set to be -0.1<Cdf,y/CDf,exact<0.4. Analytical pressure drag is CDf,exact=1600.
(a)
Reference
(b)
MS1
(c)
MS2
(d)
MS3Figure 9
Distribution of partial pressure drag for each cell front-projected onx plane on line A in Figure 7. Analytical pressure drag is CDp,exact=800.Figure 10
Distribution of partial pressure drag for each cell front-projected onx plane on lines A and B in Figure 8. Analytical friction drag is CDf,exact=1600.
(a)
Line A
(b)
Line B


#### 3.3.3. Effects of Estimation Methods for Ghost Cells

In MS, the error in the friction drag appears to arise from two sources:(1) the estimated shear stress and (2) the estimated surface area vector. These error sources are distinguished for clarity.The error in the estimated shear stress at the interface of the staircase body can be completely removed by analytically computing the shear stress in the Stokes flow, rather than approximating it by a finite difference scheme. The estimated drag with analytical shear stress at the boundary is plotted in Figure11. This solution is labeled as “MS with the analytical shear stress” in the figure. Clearly, the analytical shear stress yields much more accurate results than MS1, MS2, and MS3, in which the shear stress is numerically computed. This result illustrates that (1) the error in the estimated shear stress is significant and prevents the numerically computed force from converging, and (2) the error in the surface area vector is probably insignificant.Figure 11
Estimated drag in Stokes flow around a sphere as functions of mesh spacing using a mesh-base scheme with analytical shear stress. Analytical estimates of pressure and friction drag areCDp,exact=800 and CDp,exact=1600, respectively.
(a)
Pressure drag
(b)
Friction dragAlso note that, according to Figure11, the error in shear stress for MS2 depends on the accuracy of the ghost-cell method. Therefore, the order of accuracy in MS2 (shown in Figure 6) seems to be determined by the formal order of accuracy in the ghost-cell method.


#### 3.3.4. Surface Area Vector Estimation

We now discuss the error in the estimated surface area vector. The surface area vectors (x>0, y>0, and z>0 of the staircase body represented by the mesh faces) are summed in Figure 12. The sum converges to the exact integral over the surface area as the mesh density increases. It should be noted that the surfaces of the staircase body never exactly sum to the area of the body surface, even on very fine meshes. The estimated force acting on the body depends on the surface area vector, rather than on the surface area itself. This may explain why the staircase body approximation for integration does not matter.Figure 12
Evaluation of summation of surface area vectors and surface areas in MS.
(a)
x component of surface area vector
(b)
Surface area


#### 3.3.5. Computational Costs

The computational costs for PS and MS2 forN=4 to 128 are discussed in this subsection. Here, costs for MS1 and MS3 are almost the same as MS2. The computation is measured by a core of the computer process unit of Xeon X5450 3.0 GHz. For the PS computation, computational costs strongly depend on the number of polygon meshes. In this study, numbers of polygon mesh for each grid resolutions are determined with increasing it twice until the drag coefficient is firstly converged with 10% error. The results are shown in Table 1. The number of polygon meshes for convergence is larger than we expected, and resulting computational costs of PS are much higher than that of MS2. The computational cost of MS2 is 5–30 times lower than that of PS, and table shows that the MS2 is much more efficient than PS. This is because MS2 only requires the mesh-face loop to compute the body force while PS requires the polygon mesh loop, the number of which is much larger than that of the mesh face.Table 1
Computational costs for body force estimation of Stokes flow around a sphere.
N
PS
MS2
Computational costs [s]
Number of polygons
Computational costs [s]
4
0.0006043
64
×
64
0.0001197
8
0.002353
128
×
128
0.0003260
16
0.04947
512
×
512
0.001319
32
0.1955
1024
×
1024
0.00677
64
0.6897
2048
×
2048
0.05297
128
2.939
4096
×
4096
0.3567


## 3.1. Stokes Flow Past a Sphere

Analytical solutions to the pressure and velocity of Stokes flow around a sphere are given by [10]:(20)p=p∞-32μρu∞axr3,u=u∞-14u∞ar3+a2r2-34u∞ax2r31-a2r2,v=-34u∞axyr31-a2r2,w=-34u∞axzr31-a2r2,where a is the radius of the sphere, ρ is the mass density, μ is the viscosity, u∞ is the freestream velocity, p∞ is the reference pressure, x, y, and z are the Cartesian coordinates, and r is the distance from the origin (where the origin is the center of the sphere).The fluid is assumed to flow parallel to thex-axis; therefore, the drag force acting on the sphere is the x-component of the force defined in (1). Substituting the analytical solutions (20) into this definition, the pressure drag and the frictional drag are, respectively, obtained as (21)CDp=Fxp1/2ρu∞2A=2πaμu∞ρu∞2/2πa2=8μρu∞D=8Re,CDf=Fxv1/2ρu∞2A=4πaμu∞ρu∞2/2πa2=16μρu∞D=16Re,where CDp and CDf are the coefficients of the pressure drag and frictional drag, respectively, and Re=ρu∞D/μ is the Reynolds number. D=2a and A=πa2 define the diameter and the projected area of the sphere, respectively.


## 3.2. Computational Settings and Data

The drag force is now computed by the numerical methods described in Section2. To compute the drag forces, a sphere of the diameter D is centralized within a cubic domain of the side length L. The value of L is arbitrary but must be sufficiently large to encompass the sphere. The present study selects L=2D. The computational cells are generated by dividing the cubic domain into 2N equally sized cells in each direction. The fineness of the sphere is altered by setting the number of grid points N occupied by the sphere in each direction to N=4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048. The solid cells constituting the sphere of grid points less than 128 are illustrated in Figure 5. The PS requires an additional polygon mesh. The force computed on each Cartesian mesh properly converges with that on a spherical equally distributed polygon mesh with 3,200×3,200 points. Therefore, this polygon mesh is adopted in the remainder of the study. The analytical solutions (20) are assumed to be available only at the centers of the fluid cells, including the boundary cells.Figure 5
Spheres constructed on meshes with different resolutions. From left to right, a sphere of diameterD is constructed on mesh resolved to Δx=a/4, x=a/8,a/16,a/32,a/64, and a/128.


## 3.3. Results

### 3.3.1. Total Drag Forces

The forces estimated by PS, MS1, MS2, and MS3 at different mesh resolutions are plotted in Figure6. Here, the Reynolds number is Re=0.01 for all the cases. Correspondingly, the pressure and friction drags are, respectively, solved as CDp,exact=800 and CDp,exact=1600, using the analytical solutions. The result shows that the convergence rate for the pressure drag is first order in MS1 and second order in PS and MS2 and that the convergence rate for the friction drag is zeroth order in MS1 and first to second order in PS, MS2, and MS3. These results indicate that the proposed MS2 and MS3 can evaluate the fluid force as accurately as PS, but at much smaller computational cost and using a much simpler computational code. Thus, fluid forces can be simply and efficiently evaluated by MS2 and MS3 in the immersed boundary method. Comparing MS2 and MS3, it interestingly shows that MS2 has fewer errors than MS3. This clarifies that the MS2 based on the ghost-cell method is better for the force estimation than MS1 based on the ghost-fluid method. Apart from the success of MS2 and MS3, MS1 is shown to be unacceptable for practical use since the friction drag is not converged even when a very fine mesh is used. At Δx=0.06, sudden decrease in friction drag of PS is observed and this might be caused by canceling out two different error sources which is not clarified well in this study.Figure 6
Estimated drag in Stokes flow around a sphere as functions of mesh spacing. Analytical pressure and friction drag areCDp,exact=800 and CDp,exact=1600, respectively.
(a)
Pressure drag
(b)
Friction dragIt is worth noting that the PS results are sensitive to the density of the polygon mesh at coarser polygon mesh resolutions, whereas we use a very fine polygon mesh in this study to eliminate the error introduced by coarse polygon mesh. This implies that the polygon density requires special attention in the PS method.


### 3.3.2. Force Distributions

Then, we visualized the force distributions on the each cell projected onx, y, and z plane. The partial drag coefficients for each projected cell are defined as follows:(22)Cdp,x=∫cellpdydz1/2ρu∞2ΔyΔzCdf,y=∫cellτxy+τyzdxdz1/2ρu∞2ΔxΔz⋮.Here, integration is conducted for the surface in each projected cell. In the case of a sphere, there are two projections: front and rear projections. In this definition, the reference solution can be obtained with very fine discretization. On the other hand, MSs evaluate them in simplified ways as in the inside summation of (10). The distributions of reference solutions and MSs are shown in Figures 7 and 8. Here, the reference solution is obtained with the ten times finer submesh. The results of MS1 have strong oscillation in friction drag because this method does not compute shear stress correctly. On the other hand, MS2 and MS3 can predict them quantitatively, though they have weak oscillation. This behavior is illustrated in Figures 9 and 10, which shows the partial drag on the line indicated in Figures 7 and 8. These figures show the same trend more quantitatively. Figure 8 shows that the MS2 predict partial friction drag more accurately than MS3, as discussed above for total friction drag.Figure 7
Distribution of partial pressure drag for each cell front-projected onx plane. Contour range is set to be -0.1<Cdp,x/CDp,exact<0.8. Analytical pressure drag is CDp,exact=800.
(a)
Reference
(b)
MS1
(c)
MS2
(d)
MS3Figure 8
Distribution of partial friction drag for each cell front-projected ony plane. Contour range is set to be -0.1<Cdf,y/CDf,exact<0.4. Analytical pressure drag is CDf,exact=1600.
(a)
Reference
(b)
MS1
(c)
MS2
(d)
MS3Figure 9
Distribution of partial pressure drag for each cell front-projected onx plane on line A in Figure 7. Analytical pressure drag is CDp,exact=800.Figure 10
Distribution of partial pressure drag for each cell front-projected onx plane on lines A and B in Figure 8. Analytical friction drag is CDf,exact=1600.
(a)
Line A
(b)
Line B


### 3.3.3. Effects of Estimation Methods for Ghost Cells

In MS, the error in the friction drag appears to arise from two sources:(1) the estimated shear stress and (2) the estimated surface area vector. These error sources are distinguished for clarity.The error in the estimated shear stress at the interface of the staircase body can be completely removed by analytically computing the shear stress in the Stokes flow, rather than approximating it by a finite difference scheme. The estimated drag with analytical shear stress at the boundary is plotted in Figure11. This solution is labeled as “MS with the analytical shear stress” in the figure. Clearly, the analytical shear stress yields much more accurate results than MS1, MS2, and MS3, in which the shear stress is numerically computed. This result illustrates that (1) the error in the estimated shear stress is significant and prevents the numerically computed force from converging, and (2) the error in the surface area vector is probably insignificant.Figure 11
Estimated drag in Stokes flow around a sphere as functions of mesh spacing using a mesh-base scheme with analytical shear stress. Analytical estimates of pressure and friction drag areCDp,exact=800 and CDp,exact=1600, respectively.
(a)
Pressure drag
(b)
Friction dragAlso note that, according to Figure11, the error in shear stress for MS2 depends on the accuracy of the ghost-cell method. Therefore, the order of accuracy in MS2 (shown in Figure 6) seems to be determined by the formal order of accuracy in the ghost-cell method.


### 3.3.4. Surface Area Vector Estimation

We now discuss the error in the estimated surface area vector. The surface area vectors (x>0, y>0, and z>0 of the staircase body represented by the mesh faces) are summed in Figure 12. The sum converges to the exact integral over the surface area as the mesh density increases. It should be noted that the surfaces of the staircase body never exactly sum to the area of the body surface, even on very fine meshes. The estimated force acting on the body depends on the surface area vector, rather than on the surface area itself. This may explain why the staircase body approximation for integration does not matter.Figure 12
Evaluation of summation of surface area vectors and surface areas in MS.
(a)
x component of surface area vector
(b)
Surface area


### 3.3.5. Computational Costs

The computational costs for PS and MS2 forN=4 to 128 are discussed in this subsection. Here, costs for MS1 and MS3 are almost the same as MS2. The computation is measured by a core of the computer process unit of Xeon X5450 3.0 GHz. For the PS computation, computational costs strongly depend on the number of polygon meshes. In this study, numbers of polygon mesh for each grid resolutions are determined with increasing it twice until the drag coefficient is firstly converged with 10% error. The results are shown in Table 1. The number of polygon meshes for convergence is larger than we expected, and resulting computational costs of PS are much higher than that of MS2. The computational cost of MS2 is 5–30 times lower than that of PS, and table shows that the MS2 is much more efficient than PS. This is because MS2 only requires the mesh-face loop to compute the body force while PS requires the polygon mesh loop, the number of which is much larger than that of the mesh face.Table 1
Computational costs for body force estimation of Stokes flow around a sphere.
N
PS
MS2
Computational costs [s]
Number of polygons
Computational costs [s]
4
0.0006043
64
×
64
0.0001197
8
0.002353
128
×
128
0.0003260
16
0.04947
512
×
512
0.001319
32
0.1955
1024
×
1024
0.00677
64
0.6897
2048
×
2048
0.05297
128
2.939
4096
×
4096
0.3567


## 3.3.1. Total Drag Forces

The forces estimated by PS, MS1, MS2, and MS3 at different mesh resolutions are plotted in Figure6. Here, the Reynolds number is Re=0.01 for all the cases. Correspondingly, the pressure and friction drags are, respectively, solved as CDp,exact=800 and CDp,exact=1600, using the analytical solutions. The result shows that the convergence rate for the pressure drag is first order in MS1 and second order in PS and MS2 and that the convergence rate for the friction drag is zeroth order in MS1 and first to second order in PS, MS2, and MS3. These results indicate that the proposed MS2 and MS3 can evaluate the fluid force as accurately as PS, but at much smaller computational cost and using a much simpler computational code. Thus, fluid forces can be simply and efficiently evaluated by MS2 and MS3 in the immersed boundary method. Comparing MS2 and MS3, it interestingly shows that MS2 has fewer errors than MS3. This clarifies that the MS2 based on the ghost-cell method is better for the force estimation than MS1 based on the ghost-fluid method. Apart from the success of MS2 and MS3, MS1 is shown to be unacceptable for practical use since the friction drag is not converged even when a very fine mesh is used. At Δx=0.06, sudden decrease in friction drag of PS is observed and this might be caused by canceling out two different error sources which is not clarified well in this study.Figure 6
Estimated drag in Stokes flow around a sphere as functions of mesh spacing. Analytical pressure and friction drag areCDp,exact=800 and CDp,exact=1600, respectively.
(a)
Pressure drag
(b)
Friction dragIt is worth noting that the PS results are sensitive to the density of the polygon mesh at coarser polygon mesh resolutions, whereas we use a very fine polygon mesh in this study to eliminate the error introduced by coarse polygon mesh. This implies that the polygon density requires special attention in the PS method.


## 3.3.2. Force Distributions

Then, we visualized the force distributions on the each cell projected onx, y, and z plane. The partial drag coefficients for each projected cell are defined as follows:(22)Cdp,x=∫cellpdydz1/2ρu∞2ΔyΔzCdf,y=∫cellτxy+τyzdxdz1/2ρu∞2ΔxΔz⋮.Here, integration is conducted for the surface in each projected cell. In the case of a sphere, there are two projections: front and rear projections. In this definition, the reference solution can be obtained with very fine discretization. On the other hand, MSs evaluate them in simplified ways as in the inside summation of (10). The distributions of reference solutions and MSs are shown in Figures 7 and 8. Here, the reference solution is obtained with the ten times finer submesh. The results of MS1 have strong oscillation in friction drag because this method does not compute shear stress correctly. On the other hand, MS2 and MS3 can predict them quantitatively, though they have weak oscillation. This behavior is illustrated in Figures 9 and 10, which shows the partial drag on the line indicated in Figures 7 and 8. These figures show the same trend more quantitatively. Figure 8 shows that the MS2 predict partial friction drag more accurately than MS3, as discussed above for total friction drag.Figure 7
Distribution of partial pressure drag for each cell front-projected onx plane. Contour range is set to be -0.1<Cdp,x/CDp,exact<0.8. Analytical pressure drag is CDp,exact=800.
(a)
Reference
(b)
MS1
(c)
MS2
(d)
MS3Figure 8
Distribution of partial friction drag for each cell front-projected ony plane. Contour range is set to be -0.1<Cdf,y/CDf,exact<0.4. Analytical pressure drag is CDf,exact=1600.
(a)
Reference
(b)
MS1
(c)
MS2
(d)
MS3Figure 9
Distribution of partial pressure drag for each cell front-projected onx plane on line A in Figure 7. Analytical pressure drag is CDp,exact=800.Figure 10
Distribution of partial pressure drag for each cell front-projected onx plane on lines A and B in Figure 8. Analytical friction drag is CDf,exact=1600.
(a)
Line A
(b)
Line B


## 3.3.3. Effects of Estimation Methods for Ghost Cells

In MS, the error in the friction drag appears to arise from two sources:(1) the estimated shear stress and (2) the estimated surface area vector. These error sources are distinguished for clarity.The error in the estimated shear stress at the interface of the staircase body can be completely removed by analytically computing the shear stress in the Stokes flow, rather than approximating it by a finite difference scheme. The estimated drag with analytical shear stress at the boundary is plotted in Figure11. This solution is labeled as “MS with the analytical shear stress” in the figure. Clearly, the analytical shear stress yields much more accurate results than MS1, MS2, and MS3, in which the shear stress is numerically computed. This result illustrates that (1) the error in the estimated shear stress is significant and prevents the numerically computed force from converging, and (2) the error in the surface area vector is probably insignificant.Figure 11
Estimated drag in Stokes flow around a sphere as functions of mesh spacing using a mesh-base scheme with analytical shear stress. Analytical estimates of pressure and friction drag areCDp,exact=800 and CDp,exact=1600, respectively.
(a)
Pressure drag
(b)
Friction dragAlso note that, according to Figure11, the error in shear stress for MS2 depends on the accuracy of the ghost-cell method. Therefore, the order of accuracy in MS2 (shown in Figure 6) seems to be determined by the formal order of accuracy in the ghost-cell method.


## 3.3.4. Surface Area Vector Estimation

We now discuss the error in the estimated surface area vector. The surface area vectors (x>0, y>0, and z>0 of the staircase body represented by the mesh faces) are summed in Figure 12. The sum converges to the exact integral over the surface area as the mesh density increases. It should be noted that the surfaces of the staircase body never exactly sum to the area of the body surface, even on very fine meshes. The estimated force acting on the body depends on the surface area vector, rather than on the surface area itself. This may explain why the staircase body approximation for integration does not matter.Figure 12
Evaluation of summation of surface area vectors and surface areas in MS.
(a)
x component of surface area vector
(b)
Surface area


## 3.3.5. Computational Costs

The computational costs for PS and MS2 forN=4 to 128 are discussed in this subsection. Here, costs for MS1 and MS3 are almost the same as MS2. The computation is measured by a core of the computer process unit of Xeon X5450 3.0 GHz. For the PS computation, computational costs strongly depend on the number of polygon meshes. In this study, numbers of polygon mesh for each grid resolutions are determined with increasing it twice until the drag coefficient is firstly converged with 10% error. The results are shown in Table 1. The number of polygon meshes for convergence is larger than we expected, and resulting computational costs of PS are much higher than that of MS2. The computational cost of MS2 is 5–30 times lower than that of PS, and table shows that the MS2 is much more efficient than PS. This is because MS2 only requires the mesh-face loop to compute the body force while PS requires the polygon mesh loop, the number of which is much larger than that of the mesh face.Table 1
Computational costs for body force estimation of Stokes flow around a sphere.
N
PS
MS2
Computational costs [s]
Number of polygons
Computational costs [s]
4
0.0006043
64
×
64
0.0001197
8
0.002353
128
×
128
0.0003260
16
0.04947
512
×
512
0.001319
32
0.1955
1024
×
1024
0.00677
64
0.6897
2048
×
2048
0.05297
128
2.939
4096
×
4096
0.3567


## 4. Evaluation of Force on a Moving Sphere: Various Alignment of Body and Grid

Then, the Stokes flow around a moving sphere is considered because the immersed boundary method is often used for the moving body problems. In the moving body problems, we should consider following additional two points:(i)
change in the evaluation of the shear stress,(ii)
the various locations of the body compared with the grid.If we consider the exact solution of the Stokes flow around a moving sphere, it is exactly corresponding to the superimposition of the exact solution of the Stokes flow and the uniform velocity of moving body. If we use the MS1, MS2, MS3, and PS for the evaluation of the shear stress, it exactly corresponds to that in the static case because the uniform velocity added in the flow fields is perfectly canceled out by the boundary velocity imposed on the moving body. Therefore, the first point, that is, change in evaluation of the shear stress, does not appear for MS1, MS2, MS3, and PS, in the case we considered. On the other hand, the second point seems to be more important for moving body problems. The smaller force variation due to the alignment of the body and the grid is preferred. In this section, force variation due to the alignment of the body and the grid for each method is investigated.In this test case, the center of the sphere is shifted bycxΔx,cyΔy,czΔz from the original location used in the previous section, where -0.5<cx≤0.5, -0.5<cy≤0.5, and -0.5<cz≤0.5. We set cy=cz=0 for simplicity and force is evaluated with changing cx, assuming flow around the moving sphere through the quiescent air.The force is evaluated with the grid resolutions ofN=4, 8, 16, 32, and 64 in this section. The results of N=16 are shown in Figure 13. The results show that the friction drag of MS1 has large oscillations depending on the alignment of the body and the grid, while the other results have much smaller oscillations. Interestingly, the oscillation is significantly suppressed by PS, though the accuracy of its averaged value is worse than MS2 in this case.Figure 13
Effects of alignment of body and grid on estimated drag forN=16. Analytical estimates of pressure and friction drag are CDp,exact=800 and CDp,exact=1600, respectively.
(a)
Pressure drag
(b)
Friction dragIn order to quantitatively evaluate the force oscillation by the alignment of the body and the grid, the root mean square value is introduced in Figure14, where the bracket represents the averaging operator for -0.5<cx≤0.5. The root mean square of the force oscillation normalized by the analytical force for each case is shown by the thin lines in Figure 14, together with the error of averaged force compared with the analytical force by the thick lines. Here, MS1 has very large oscillations due to the alignment as well as very large errors in averaged value as discussed in the previous section. MS2 and MS3 have similar or smaller oscillations than the errors in averaged value and their oscillations are considered to be acceptable. Besides, PS has very small oscillation and shows the best characteristics for the grid alignment.Figure 14
Estimated averaged drag and oscillation due to alignment of body and grid in Stokes flow around a moving sphere as functions of mesh spacing. Analytical estimates of pressure and friction drag areCDp,exact=800 and CDp,exact=1600, respectively.
(a)
Pressure drag
(b)
Friction dragThe test case shown in this section illustrates that the proposed MS2 and MS3 methods have similar or smaller oscillation than error and can be used for the moving body problems, though PS shows the best characteristics for the grid alignment of the body and the grid.


## 5. Estimating the Force in the Flow around a Sphere at the Reynolds Number 100

To evaluate MS at the higher Reynolds number, the force was estimated in fluid flowing around a sphere at the Reynolds number 100. Since no analytical solution was available for such higher Reynolds number flows, numerical solution was obtained by using an incompressible flow solver developed by Onishi et al. [11]. This solver adopts a simple marker and cell algorithm with collocated variable layout. The grid system is based on the building cube method proposed by Nakahashi [12]. Space is discretized by a second-order finite difference scheme [13] and the simple immersed boundary method based on the ghost-fluid based method proposed by Sato et al. [14]. This test case includes the errors in the computed flow fields, and therefore the result does not directly indicate the performance of the force estimation method.The forces evaluated by PS, MS1, MS2, and MS3 are shown in Figure15. The result shows that the forces estimated by PS, MS2, and MS3 seem to converge to 1.09, which is consistent with the existing data obtained by the empirical relation derived from a number of experiments [15] and the numerical computations [7, 16].Figure 15
Estimated drags in flow around a sphere at the Reynolds number 100 as functions of mesh spacing.Interestingly, MS3 has less error in a coarse grid. This might be because fluid dynamic computation and force evaluation methods are consistent. On the other hand, the magnitudes of the errors in PS and MS2 are almost identical on the coarser and finer meshes. This result illustrates that MS2 and MS3 are adequate methods for evaluating forces acting on flows around objects in the immersed boundary method, although only the information of mesh and flow variables is used in the flow and boundary cell computations.


## 6. Applicability to More Complex Problems: Body with Sharp Edges, Multibody Treatment, and Compressibility

In this paper, we introduced the mesh-based force evaluation methods in the immersed boundary methods. When it is applied to the complex body problems, there are the sharp edge problems or other multiple body problems. However, this method can be applied to such complex body problems, as far as the base immersed boundary method can treat such bodies. This is because the proposed method just uses the cell-face value which is given by the base immersed boundary method, regardless of whichever method, such as a ghost-cell method or a ghost-fluid method, is used as the base immersed boundary method. Therefore, the applicability of this method to the complex body problem is highly depending on the base immersed boundary methods and the further discussion is left for the other papers which discuss applicability of the base immersed boundary method to the complex body problems.In addition, it is worth emphasizing that the present method can be applied to compressible flows as well as incompressible flows. This is because the present method does not rely on the details of fluid flows (i.e., the divergence of velocity field). Although all the numerical results shown in this paper have been obtained by incompressible flows, several studies can be found in the literature on this subject [17, 18].


## 7. Conclusions

We have proposed a simple method for evaluating the forces acting on flows around bodies in the immersed boundary scenario. This method has been developed by employing a novel mesh-face integration method and an extrapolation method for evaluating pressure and shear stresses at the mesh faces, such as the first-order, ghost-cell, or ghost-fluid methods. The present method is, in principle, advantageous over the conventional methods based on control volumes in that pressure and shear stress can be evaluated separately.Moreover, we have applied the present method to the computation of the drag force acting on a sphere in Stokes flow and have investigated the effects of grid spacing and extrapolation methods on the errors originating from the present force estimation method by using the existing analytical solutions. In addition, we have addressed the computational costs. As a result, the accuracy of the proposed mesh-based scheme has been proven to be comparable to that of the polygon-based scheme, which is commonly adopted in straightforward implementation. This indicates that the proposed scheme works better than the polygon-based one when complex geometries are involved, since its implementation is simple and its computational cost is low.The error sources in the proposed implementation are sourced from(1) the surface area vector of the staircase body shape and (2) the approximated shear stress. Of these, error in the evaluated shear stress dominates and is significant. If the shear stress is appropriately evaluated, the fluid force can be accurately obtained by summing over the mesh faces, because the surface area vector components converge with increasing grid density while the surface area does not. The shear stress is adequately evaluated by the second-order finite differencing scheme with the ghost-cell or ghost-fluid method. Sometimes, it is difficult to estimate the shear stress accurately with this method by its complex shape. It should be noted that this difficulty is caused by the immersed boundary methods themselves and the present idea using the staircase integration does not have difficulty.



---

*Source: 8314615-2017-07-19.xml*