Extension of local domain-free discretization method to simulate 3D flows with complex moving boundaries

Computers & Fluids 64 (2012) 98–107

Contents lists available at SciVerse ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Extension of local domain-free discretization method to simulate 3D flows with complex moving boundaries C.H. Zhou a,⇑, C. Shu b a b

Department of Aerodynamics, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, PR China Department of Mechanical Engineering, National University of Singapore, 119260 Singapore, Singapore

a r t i c l e

i n f o

Article history: Received 31 October 2010 Received in revised form 10 May 2012 Accepted 15 May 2012 Available online 26 May 2012 Keywords: Domain-free discretization Non-boundary-conforming method Moving boundary Carangiform swimming Immersed boundary Unsteady flows

a b s t r a c t This paper is the first endeavor to present the local domain-free discretization (DFD) method for the solution of the three-dimensional Navier–Stokes equations. The computational domain may contain complex moving boundaries. The strategy of DFD is that the discrete form of partial differential equations at an interior point may involve some points outside the solution domain. The functional values at the exterior dependent points are evaluated by the approximate form of solution near the boundary. Compared to the previous work, the tedious task of constructing new interpolation tetrahedrons is eliminated, and this reduces the complexity of DFD implementation. An efficient algorithm for classifying mesh points is also presented. Simulation of flow around a stationary sphere is used to validate the numerical method, and three distinct flow regimes have been obtained with varied Reynolds numbers of up to 300. The ability of the method for flows with complex moving boundary is demonstrated by simulating flows over an undulating fish-like body. The results of force coefficient, structure of wake patterns and propulsive efficiency at critical Strouhal number have been presented. All predictions show a good agreement with the reference data. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction In [1], Shu et al. proposed a domain-free discretization (DFD) method to solve partial differential equations (PDEs) on irregular domains. So far, the DFD method has been successfully applied to simulate two-dimensional flows around moving bodies [2–4] and three-dimensional inviscid compressible flows around a stationary body [5]. This paper is the first endeavor to present the local DFD method for simulation of three-dimensional flows over geometrically complex and moving bodies. In the DFD method, the PDE is discretized at all grid nodes inside the solution domain, but the spatial discretization at an interior node may involve some nodes outside the solution domain. The critical issue for successful implementation of the DFD method is how to evaluate the functional values at the exterior dependent nodes, i.e., to construct some approximate form of the solution in the local region. In the earlier DFD method [1], the approximate form of solution is pursued along the whole mesh line that only involves two boundary points. This way is not applicable for complex domains. To make the method be more general, the local DFD was developed in [2–5]. In the local DFD, the low order schemes are adopted for spatial discretization and also for the approximate form of solution near the wall boundary. The functional values at ⇑ Corresponding author. E-mail address: [email protected] (C.H. Zhou). 0045-7930/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2012.05.012

the exterior dependent points are evaluated by extrapolation along the normal direction to the wall in conjunction with boundary conditions and the simplified momentum equation in the vicinity of the wall boundary. As described above, the DFD method is a kind of non-boundaryconforming methods since a solid boundary can be immersed on a grid. In recent years, these methods are attracting more and more attentions, especially for problems involving complex moving boundaries. Using non-boundary-conforming methods, flow equations can be solved on a fixed grid for moving boundary problems, and some critical issues associated with grid updating at each time step, such as grid quality and interpolation error due to solutiontransferring, can be removed. Non-boundary-conforming methods can be roughly classified into two major categories. The first category consists of methods in which ‘‘continuous forcing’’ is employed and the effect of boundaries is transmitted into the momentum equations as a forcing term. The original immersed boundary method of Peskin [6] and the later versions such as those of Peskin and McQueen [7,8], Goldstein et al. [9] and Saiki and Biringen [10] fall in this category. The advantage of these methods is that they are formulated independent of the spatial discretization. However, the smeared boundary is an inherent drawback of this category of methods. The second category consists of methods that employ ‘‘discrete forcing’’. These forces are placed at points that adjoin the boundary in order to impose the velocity boundary conditions on the solid surface. The

99

C.H. Zhou, C. Shu / Computers & Fluids 64 (2012) 98–107

forcing scheme depends much on the spatial discretization. The obvious advantage of this category of methods is that the immersed boundary can be tracked as a sharp interface. The methods of Mohd-Yusof [11], Fadlun et al. [12], Gilmanov and Sotiropoulos [13], Marella et al.[14], Udaykumar et al. [15] and Mittal et al. [16] fall in the second category. The basic idea of DFD is different from the above methods. It is based on the extension of PDE solution and the boundary conditions are imposed via the approximate form of solution in the local region. In this paper, the local DFD method is extended to solve the three-dimensional, unsteady, incompressible Navier–Stokes equations. The computational domain may contain geometrically complex boundaries moving with prescribed kinematics. Compared to the previous work in [4,5], there is no need to construct new interpolation tetrahedrons (triangles for two-dimensional case) in the procedure of evaluation of flow variables at the exterior dependent points. This improvement is very important for 3D computations since the construction of new interpolation tetrahedrons leads to a large number of possible intersections between the fixed grid and the wall boundary, and various irregular sub-cells are formed. If the boundary is moving, the construction has to be performed at each time step. So, removing the construction of new interpolation tetrahedrons reduces greatly the complexity of DFD implementation and also the computational cost. The present numerical method is validated for stationary immersed boundary problems by simulating the flow around a stationary sphere, and three distinct flow regimes have been obtained with varied Reynolds numbers of up to 300. For moving boundary problems, we test the method by simulating the flow induced by a sphere rotating in a fluid at rest. All predictions show a good agreement with the benchmark numerical results or experimental data available in the literature. Simulation of the fluid-structure interactions that dominate the dynamics of flows, such as flows over swimming and flying animals and flows over micro air vehicles (MAVs), poses a challenge to numerical techniques, and is currently at the forefront of ongoing work in computational fluid dynamics. One can find recent applications of non-boundary-conforming methods in simulation of flows with complex moving boundaries in [13,16–18]. In this paper, we demonstrate the ability of the local DFD method for flows with complex moving boundary by simulating flows over an undulating fish-like body. The results of force coefficient, structure of wake patterns and propulsive efficiency at critical Strouhal number have been presented, which confirm the numerical results of Borazjani and Sotiropoulos [17], elucidate some results of previous experiments and help understand the propulsion of undulatory swimmers. This test simulation underlines the potential of present local DFD method as a powerful simulation tool for bio-fluid problems. For the remainder of this paper, the content is arranged as follows. In Section 2, the governing equations and numerical approximation are described. In Section 3, we discuss the implementation of the local DFD method in three-dimensional computations. In Section 4, an efficient algorithm for classifying the nodes of the Eulerian mesh is presented. Section 5 is devoted to the numerical experiments. Finally, in Section 6, summary and conclusions are given.

where w is the vector of flow variables, fi and si are the convective and viscous flux vectors respectively

2

p

3

2

6u 7 6 17 w ¼ 6 7; 4 u2 5

ui

3

6 u u þ pd 7 1i 7 6 1 i fi ¼ 6 7; 4 u2 ui þ pd2i 5

u3

2 si ¼

0

3

7 1 6 6 si1 7 6 7 4 Re si2 5

ð2Þ

si3

u3 ui þ pd3i

m

and I ¼ diagð 0 1 1 1 Þ is the modified identity matrix annihilating the temporal derivative of pressure from the continuity equation. In (2), ui denotes the Cartesian velocity components, p the pressure, sij the viscous stress tensor, and Re the Reynolds number. We suppose that X  R3 is a connected open set containing a solid body, and denote the boundary of X by C. With h a space discretization step, a tetrahedrization Th of X is introduced. The wall boundary is superimposed upon the computational mesh. The essence of DFD is that the discrete form of governing equations at an interior point can involve some points outside the solution domain. It still needs a numerical approach to do discretization and transfer the differential equation into a discrete form. The spatial discretization employed in this work is the direct three-dimensional extension of the two-dimensional Galerkin finite element approach used in [4], which is second-order accurate. Here, we describe the extension in brief. Let F denote the convective flux tensor, and S the viscous flux tensor. The Navier–Stokes equations can be rewritten as

Im 

@w þrF¼rS @t

ð3Þ

Applying the Galerkin finite element approximation to (3) and employing the concept of a lumped mass matrix, one obtains

Im  Xi

n n X @wi X FA þ FB þ FC 4 e S  DSABC ¼  DSABC  3 @t 3 e¼1 e¼1

In the above equation, the summation is over all the tetrahedrons that contain the vertex i, and Xi represents the volume sum of the tetrahedrons. These tetrahedrons define the influence domain of i. As illustrated in Fig. 1, DSABC represents the directed (normal) triangle area of the face of each tetrahedron e on the outer boundary of the influence domain. FA, FB and FC are the convective fluxes at the three vertices of this triangle, and Se is the constant viscous flux over e. We can see that there is no imposition of wall boundary conditions in the spatial discretization. The wall boundary conditions will be implemented via evaluating the flow variables at exterior dependent points as discussed in the next section. The semi-discrete Eq. (4) are integrated in time via a third-order accurate dual-time-stepping scheme combined with artificial compressibility approach [4,19]. As we know, the standard Galerkin finite element approximation is equivalent to the central difference. In order to prevent odd-even decoupling, a third order artificial dissipation term is

C B

2. Governing equations and numerical approximation We will consider a homogeneous, laminar flow without body force effect. The governing equations are the time-dependent incompressible Navier–Stokes equations in the non-dimensional form

Im 

@w @f i @si þ ¼ ; @t @xi @xi

I SABC

A i ¼ 1; 2; 3

ð4Þ

ð1Þ Fig. 1. Directed area of influence-domain-boundary face of a tetrahedron.

100

C.H. Zhou, C. Shu / Computers & Fluids 64 (2012) 98–107

added to the spatial residual [19]. Following the previous work for two-dimensional flows [4], the artificial dissipation term can be constructed as follows

Di ðwÞ ¼ Pr1 

n X a þ ak 2 ðr wi  r2 wk Þ e i 2 k¼1

h K maxðjV 1 j; jV 2 j; jV 3 jÞ

Fluid

ð5Þ Wall

where Pr is the preconditioning matrix in artificial compressibility approach, e is an empirical coefficient which is taken as 1/128, a is a factor proportional to an isotropic value of the maximum eigenvalue at each node, and n represents the number of edges meeting at the node i. In (5), the undivided Laplacian of w at the node i, can be approximated as the summation of the differences of w along all edges meeting at i. On the outer boundary C, approximate non-reflecting far-field boundary conditions are constructed based on the linearized characteristics approach to improve accuracy and the rate of convergence [19]. To address the issue of ‘‘freshly cleared nodes’’ due to the body motion [4], the physical time step is limited as follows

Dt 6

1

ð6Þ

where Vi(i = 1, 2, 3) are the Cartesian components of the body-velocity vector, K is the order of the approximation of time derivatives, and h is the minimum mesh interval in the vicinity of the solid wall. 3. Implementation of local DFD in three-dimensional computations DFD is based on the extension of PDE solution to exterior of the solution domain and the discrete form at an interior point may involve some points outside the solution domain. Therefore, all the mesh points are classified into three categories: interior points, exterior dependent points and exterior independent points that are blanked out of the computation. The key process of the local DFD method is how to evaluate the functional values at the exterior dependent points, and in this process, the wall boundary conditions will be imposed. The implementation of the local DFD, described as following, is the extension and improvement of the two-dimensional version presented in [4]. Compared to the previous work, the tedious task of constructing new interpolation tetrahedrons (triangles in two-dimensional case) is eliminated, and thus the complexity of implementation is reduced. According to the spatial discretization and the artificial dissipation presented in the last section for three-dimensional case, the calculation of the vector of flow variables at any point inside the solution domain depends on the vector and its undivided Laplacian at each surrounding point connected to this interior point by an edge of tetrahedron. In the practical computations, the coefficient e in Eq. (5) is set to be zero for all the edges intersecting with the solid wall. This treatment means that only evaluation of the vector of flow variables at an exterior dependent point is required in calculations, and there is no need to compute its undivided Laplacian. All our numerical experiments show that this cancelation of the action of artificial dissipation term for the wall-intersected edges does not affect the stability of numerical computation. For an edge intersecting with the wall boundary, the interior end point is a computed node at which all flow variables are obtained by solving the governing equations, and the exterior end point is one of the dependent points of this computed node. All this kind of interior computed nodes near the wall boundary has at least one dependent point outside the solution domain, i.e. the exterior dependent point. An example is shown in Fig. 2. In the local DFD, the values of flow variables at an exterior dependent point are extrapolated from the flow field or deter-

3

2

Solid 4

1: Interior computed node 2, 3, 4: Exterior dependent points

Fig. 2. Interior computed node, exterior dependent points.

mined by the local simplified flow equations, along the normal line passing through the exterior dependent point. Therefore, some points on the normal to the boundary and inside the solution domain should be constructed for the extrapolation and the solution of the local flow equations. These points may not be the mesh points, so we call them fictitious points. For a given exterior dependent point, its mirror image to the wall is chosen initially as the candidate of fictitious point. Denote the tetrahedron in which the mirror image locates by sM. If the four vertices of sM are all inside the solution domain as illustrated in Fig. 3, the mirror image of an exterior dependent point is finally chosen to be its fictitious point. In this case, by the linear interpolation over sM, the values of flow variables at the fictitious point can be obtained. For example, the pressure at the fictitious point is computed by

pf ¼

4 X /k ðxf ; yf ; zf Þpk

ð7Þ

k¼1

where the summation over k refers to the four vertices of tetrahedron sM, and /k(xf, yf, zf) is the linear interpolation function for each vertex at the fictitious point. However, when some of the four vertices of sM are outside the solution domain as shown in Fig. 4, the values of flow variables at any exterior vertex (for example, point 4 in the figure) of sM are not known currently. In this case, the mirror image will not be chosen to be the fictitious point of the given exterior dependent point because flow variables at the mirror image cannot be calculated by the interpolation over sM. As illustrated in Fig. 4, we define the wall-closest intersection point between the normal line and the interior faces to be the fictitious point. Here, an interior face means the tetrahedron face the three vortices of which are all inside the solution domain. The tetrahedron face (a triangle) in which the fictitious point locates is denoted by sF, as shown in Fig. 4. By linear interpolation over sF, the flow variables at the fictitious point can 1

3 2 Fictitious point of D (Mirror image)

Fluid 4

Solid

τM : 1234

D (Exterior dependent point) Fig. 3. Definition of the fictitious point of an exterior dependent point D (Case 1).

C.H. Zhou, C. Shu / Computers & Fluids 64 (2012) 98–107

Du 1 2  n ¼ rp  n þ r un Dt Re

Fictitious point of D

1

3

2

Mirror image of D Mi Fluid

τ F : 123 τ M : 1234 4 Solid

Normal line (dash dot) D (Exterior dependent point)

Fig. 4. Definition of the fictitious point of an exterior dependent point D (Case 2).

be evaluated. For example, the pressure at the fictitious point is computed by

pf ¼

3 X

uk ðxf ; yf ; zf Þpk

ð8Þ

k¼1

where the summation over k refers to the three vertices of triangle sF, and uk(xf, yf, zf) is the linear interpolation function at the fictitious point. The above prescription for fictitious points will avoid the overshoot of an extrapolated velocity. The values of flow variables at the fictitious points, at each time step, will be used in conjunction with the wall boundary conditions and the simplified momentum equation in the direction normal to the wall boundary to update continually the values of flow variables at the exterior dependent points. Now, we give formulations for the evaluation of flow variables at an exterior dependent point with reference to Fig. 5. For viscous flows, the velocity at an exterior dependent point D is determined from the non-slip boundary condition, that is, at the boundary intersection point W, the velocity of fluid is equal to the velocity of body movement

uw ¼ V w

ð9Þ

where uw and Vw are the velocity vectors at W, respectively of fluid flow and body motion. Then, the velocity of flow at the exterior dependent point D can be computed using linear extrapolation between the fictitious point F and the boundary intersection point W

ud ¼

jFDjVw  jWDjuf jFWj

ð10Þ

where ud and uf represent the velocity vectors at D and F respectively, and the absolute value means the distance between two points. To determine the pressure at the exterior dependent point D, one can consider the projection of the momentum equation of a fluid particle in the direction normal to the wall boundary

101

ð11Þ

where n is the outward normal unit vector. With the present solution-extension procedure, the viscous term in (11) turns out to be zero as all velocity components are linear in the vicinity of the wall. Then, Eq. (11) reduces to be the projection of the inviscid momentum equation in the direction normal to the wall boundary. According to the no-slip boundary condition, the material derivative of velocity at the left hand side of (11) is equal to the acceleration of the material point on the solid boundary which can be calculated directly from the known body movement. Finally, (11) is transformed to be

  @p dV n ¼ @n dt c

ð12Þ

where V is the known velocity vector of body motion and c represents the solid boundary. After approximating the derivative in the left hand side of (12) by the difference of pressure between the fictitious point and the exterior dependent point, we obtain

pd ¼ pf þ jFDj

  dV n dt w

ð13Þ

where pd and pf represent the pressure at the exterior dependent point D and the fictitious point F respectively. In the previous work for two-dimensional problems [4], the fictitious point corresponding to a given exterior dependent point is always defined to be its mirror image to the wall boundary. If the fictitious point locates in a triangle that intersects with the wall, a new interpolation triangle that does not belong to the triangulation of X needs to be constructed locally to evaluate the flow variables at the fictitious point. For moving-boundary problems, this construction must be performed at each time step. In the three-dimensional situation, the local construction of new interpolation tetrahedrons becomes more complicated and computationally costly. In the present work, we do not extend this treatment to the three-dimensional case. As we have done, if the tetrahedron that contains the mirror image intersects with the body surface (some of the four vertices of the tetrahedron are outside the solution domain), another point is chosen to be the corresponding fictitious point. Compared to the mirror image, this point is farther away from the body surface and the velocity-extrapolation stencil becomes a little larger, but the tedious task of constructing new interpolation tetrahedrons is avoided. Numerical experiments show that the effect of the slight enlargement of extrapolation stencil on the accuracy of computed results can be ignored. There exist some special exterior dependent points at which the flow variables may be multi-valued. These points are associated with the thin part of a body, the width of which is smaller than two or one grid interval. The discussion and treatment of these multi-valued points can be found in [3,5]. 4. Algorithm for mesh point classification

WALL

D

W FLOW FIELD

F

n Fig. 5. Local system for evaluation of flow variables at an exterior dependent point.

As done in any sharp-interface method, the first step in the implementation of the DFD method is to classify the nodes of Eulerian mesh into the three categories as mentioned before. For moving boundary problems involving complex three-dimensional bodies, this classification is critical for the overall efficiency of computation as it must be performed at each time step. Here, an efficient classification algorithm is presented. We employ the ray-intersection method [20] to determine whether a mesh point is inside or outside the solution domain. For simplicity, we illustrate this approach for the two-dimensional case. The extension to three-dimensional case is straightforward.

102

C.H. Zhou, C. Shu / Computers & Fluids 64 (2012) 98–107

require that any grid clustering near the body must be maintained to the far-field, and the direct application of Cartesian meshes will lead a very large node number in the calculation of flows around a body. In order to address this problem, we generate a coarse Cartesian mesh at the beginning, and then an initial uniform tetrahedral mesh is constructed by dividing hexahedral cells of this coarse Cartesian mesh. The final tetrahedral mesh used in calculations is obtained by refining the initial tetrahedral mesh gradually from the body surface to the far-field. We flag the tetrahedrons for refinement in terms of their distances to the body surface, and therefore the refined mesh is geometrically adaptive. The details of the local refinement technique can be found in [5].

Fig. 6. Illustration of ray-intersection method for two-dimensional case.

For two-dimensional problems, the boundary of an object is described by a polygon. At first, we define a smallest rectangle that contains all vertices defining the body surface. Any mesh point outside the rectangle is inside the solution domain and consequently it is an interior point. We only need to classify the points inside the rectangle. Then, starting from a given point in the bounding rectangle, cast a random half-infinite ray and count the number of intersections between the half ray and polygon edges. As shown in Fig. 6, if the number is odd then the point is located inside the polygon (outside the solution domain), otherwise it is located outside the polygon (inside the solution domain). After having identified whether a mesh point is an interior or exterior point, it is easy to find the exterior dependent points from all exterior points. If one of the two end points of a tetrahedron edge is interior and another is exterior, this edge must intersect with the wall boundary and the exterior end point must be an exterior dependent point. Till now, by applying the ray-intersection method in the bounding rectangle, we have classified all mesh points into three categories: interior points, exterior dependent points, and exterior independent points. The above classification procedure is performed only at the first time step. From the second time step onwards, the application of the ray-intersection method is not necessary for all points in the bonding rectangle. As presented in Section 2, to address the issue of ‘‘freshly cleared nodes’’, the physical time step is limited by (6). The formula means that the solid object transverses at most one entire mesh interval within one time step for K P 1. As a result, only the status of the end points of an edge that intersects with the immersed boundary at the previous time step may change at the current time step. Therefore, from the second time step onwards, the application of the ray-intersection method is limited to the end points of this kind of edges. 5. Numerical experiments

(a) 1

(b)

Y

This study validates the current DFD method for three-dimensional flows around stationary boundaries by simulating flows past a sphere immersed in an unbounded fluid. There are distinct regimes in the flow that can be identified by the Reynolds number (based on the sphere diameter and the free stream velocity) [21,22]: steady axisymmetric flow (Re < 200), steady asymmetric flow (210 < Re < 270) and unsteady flow (Re > 280). Four different Reynolds numbers (Re = 50, 150, 250, 300) are considered in this study. The computational domain is 15D  10D  10D, where D is the diameter of sphere. The used mesh contains 518,319 nodes and the mesh spacing in the vicinity of sphere is about 0.015D. The surface of sphere is discretized with a triangular mesh consisting of 16,200 elements. The physical time step is taken to be Dt = 102. A uniform flow is employed as an initial condition for all test simulations. Fig. 7a and b illustrates the streamlines through the symmetry plane (z = 0) for Reynolds numbers of 50 and 150. As expected, the flow has a steady, axisymmetric wake in these two cases. Fig. 8a and b shows streamlines constructed from in-plane velocity vectors in the (x–z) and (x–y) planes for Re = 250. It is clear from Fig. 8a that the flow field is symmetric about the (x–y) plane. Since Fig. 8b corresponds to the (x–y) symmetry plane, the streamlines in this figure are the real ones. It is apparent from Fig. 8b that the toroidal vortex has tilted and the top and bottom of the vortex ring are different from each other. Thus the flow is asymmetric. Fig. 9a– h illustrates the streamlines constructed from in-plane velocity vectors in the (x–y) and (x–z) planes for Re = 300, with time interval of one quarter period. It can be seen clearly from Fig. 9a–d that the flow is unsteady at this Reynolds number. Fig. 9e–h shows that the (x–y) plane remains a symmetry plane of flow throughout the period. The movement of the reversed flow region is evident here. From Fig. 9a–d, it can also be seen that the vortices are always shedding from the upper side of the back of sphere, and this phenomenon corresponds to the non-zero time-averaged lift. Table 1 shows our results for the drag and lift coefficients and Strouhal number (related to vortex shedding frequency) versus referenced numerical solutions and experimental data. Note that CD

Y

In the present work, tetrahedral meshes are obtained by dividing the hexahedral cells of Cartesian meshes. Cartesian meshes

5.1. Flow over a stationary sphere

1

0

-1 -1

0

1

X

2

0

-1 -1

0

1

X

Fig. 7. Axisymmetric streamlines past the sphere: (a) Re = 50; (b) Re = 150.

2

103

(a)1

(b)1

Z

Y

C.H. Zhou, C. Shu / Computers & Fluids 64 (2012) 98–107

0

-1 -1

0

1

0

-1 -1

2

0

1

X

2

X

Fig. 8. Streamlines of projected velocity vectors at Re = 250: (a) the (x–z) plane; (b) the (x–y) plane.

(e) 1

(a) 1

Z

Y

0.5 0

0

-0.5 -1 -1

0

1

-1 -1

2

0

X

1

2

X

(b) 1

(f) 1 Z

Y

0.5 0

0

-0.5 -1 -1

0

1

-1 -1

2

0

X

1

2

X

(c) 1

(g) 1 Z

Y

0.5 0

0

-0.5 -1 -1

0

1

-1 -1

2

0

X

2

(h) 1 Z

Y

(d) 1 0

-1 -1

1

X

0

1

2

X

0

-1 -1

0

1

2

X

Fig. 9. Streamlines of projected velocity vectors at Re = 300 for every quarter period: (a)–(d) the (x–y) plane; (e)–(h) the (x–z) plane.

and CL for Re = 300 are the time-averaged values. As shown in the table, the present method can accurately predict steady and unsteady characteristics of laminar flows around a sphere up to Re = 300.

5.2. Flow over a rotating sphere In this validation case, we simulate the flow induced by a sphere of radius R rotating at constant angular velocity x about

104

C.H. Zhou, C. Shu / Computers & Fluids 64 (2012) 98–107

Table 1 Drag and lift coefficients and Strouhal number for flows over a sphere. References

Clift et al. [22] Johnson and Patel [21] Marella et al. [14] Choi et al. [18] Kim et al. [23] Present work

Re = 50

Re = 150

Re = 250

CD

CD

CD

CL

CD

CL

St

1.57 1.57 1.56

0.89 0.90 0.85

0.70

0.061

0.069

0.70 0.70 0.69

0.052 0.059 0.058

0.656 0.621 0.658 0.657 0.650

0.137 0.133 0.134 0.134 0.131

1.56

0.86

Z

a diameter directed along the z-axis in an unbounded incompressible fluid at rest. The Reynolds number is defined as Re = xR2/m (m is the kinetic viscosity). For the Reynolds number in the range of Re = 1  100, benchmark solutions have been reported by Dennis et al. [24], who solved the steady axisymmetric Navier–Stokes equations in spherical, polar coordinates using a vorticity-streamfunction formulation. Even though the flow exhibits steady and axisymmetric solutions within this range of Reynolds number, we solve the full three-dimensional and unsteady flow problem with the sphere starting to rotate impulsively from rest relative to the stationary and non-boundary-conforming grid. The computational domain is a 10R  10R  10R cube. The mesh contains 369,274 nodes and the spacing in the vicinity of sphere is about 0.015D. The surface of sphere is discretized with a triangular mesh consisting of 16,200 elements. The physical time step is Dt = 102T where T is the rotation period. Calculations have been carried out for three Reynolds numbers, Re = 20, 50, and 100 on the same mesh. As discussed in [24], the flows from the northern and southern hemispheres collide at the equator forming an equatorial jet that is directed radially outward and the streamlines have to close to form two toroidal vortex rings located symmetrically with respect to the equator. The computed steady-state flow pattern at the y-plane of symmetry for Re = 100 is shown in Fig. 10. It is very similar to the result reported by Dennis et al. [24]. To compare the present results quantitatively with

2

Re = 300

0.068 0.067 0.066

those in [24,13], we compute the radial velocity profiles along the equator and the angle hs between the z-axis and the radius passes through the center of the toroidal vortex ring (see Fig. 10 for definition). The comparison for hs is shown in Table 2, and the comparison for the velocity profiles is shown in Fig. 11. The three numerical solutions are in very good agreement. 5.3. Simulation of fish-like swimming In this subsection, we seek to demonstrate the applicability of the local DFD method to simulate flows with complex, flexible immersed boundaries and underline its potential as a powerful simulation tool for bio-fluid problems. We simulate the flow induced by a fish-like, flexible object of length L moving in a straight line with constant axial velocity U and undulating in the lateral direction with a characteristic frequency f. The geometry of the three-dimensional body is provided by Sotiropoulos and Borazjani [25], which is based on measurements obtained from a real mackerel fish. Since the emphasis in numerical simulation is on the body/caudal fin mode of aquatic swimming, only the caudal fin is retained in the model fish. The body surface is discretized with triangular elements. The threedimensional view of the model fish is shown in Fig. 12. In the following descriptions, lengths are scaled with L, velocities are scaled with U, and time is scaled with f1. The kinematics for the swimmer is taken from Borazjani and Sotiropoulos [17]. The swimming motion is prescribed by specifying the lateral displacement of the fish backbone as a function of time, which is given in terms of a backward-traveling wave of varying amplitude as follows

hðx1 ; tÞ ¼ aðx1 Þ sinðkx1  xtÞ 1

-3

-2

-1

0

θS 0

ð14Þ

3

1

2

X 3

2.5

Re=100

2

u Re 1/ 2

-1

-2

1.5

Re=50

1 Fig. 10. Steady-state streamlines at the y = 0 plane for Re = 100. Re=20

0.5 Table 2 Comparison of the angle hs.

0

References

Re = 20

Re = 50

Re = 100

Dennis et al. [24] Gilmanov and Sotiropoulos [13] Present work

62.6 61.6 62.0

69.4 68.2 68.8

73.8 72.1 73.0

0

2

4

6

8

10

12

(x/R − 1) Re 1/ 2 Fig. 11. Radial velocity profiles along the equator. Solid lines: present computations; dashed lines: numerical results of Gilmanov and Sotiropoulos [13]; filled squares: benchmark data of Dennis et al. [24].

105

C.H. Zhou, C. Shu / Computers & Fluids 64 (2012) 98–107

aðx1 Þ ¼ a0 þ a1 x1 þ a2 x21

Z

X

where the coefficients a0 = 0.02, a1 = 0.08, and a2 = 0.16. With this amplitude envelope, the maximum displacement at the tail will be amax = 0.1, which gives hmax = 0.1. The computational domain is a 8L  3L  2L cuboid, discretized with a grid containing 772,172 nodes and 4,550,644 tetrahedral cells. The mesh size near the body surface is about 0.008L. The tail beat period T is divided into 500 time steps, i.e. the physical time step Dt = T/500. The component of the instantaneous hydrodynamic force along the swimming direction can be obtained by integrating the pressure and the viscous forces acting on the body as follows

Y

Fig. 12. Three-dimensional view of mackerel body (the geometry is supplied by Sotiropoulos and Borazjani [25]).

FðtÞ ¼ 

g ¼ TU=ðTU þ P L Þ

St=0.3

2

-1

St=0.5 Borazjani et al. Present

1.5 1 0.5 0 -0.5 -1

-1.5 -2 11

-1.5 11.5

12

12.5

-2 11

13

11.5

t/T

3

St=0.6

Borazjani et al. Present

2.5

1.5 1 0.5 0 -0.5 -1

13

Borazjani et al. Present

St=0.7

2 1.5 1 0.5 0 -0.5

-1.5 -2 11

12.5

3.5

Force coefficient

Force coefficient

2

12

t/T

3 2.5

ð18Þ

where T is the time-averaged thrust force, and PL is the time-averaged power loss due to lateral undulations. In our applications, Eq. (18) can only be used to compute the propulsive efficiency at the critical Strouhal number St for which the time-averaged force acting on the fish body is zero and constant-speed inline swimming is possible. Following the work of Borazjani and Sotiropoulos [17], the instantaneous net force F(t) is decomposed into the thrust contribution T(t) and the drag contribution D(t) as follows 3

-0.5

ð17Þ

It is useful to define a Froude propulsive efficiency based on the thrust force for constant-speed inline swimming as follows (see [17,26])

2.5

0

ð16Þ

C F ¼ FðtÞ=qU 2 L2

1.5

0.5

ðpn1 þ s1j nj ÞdA

where nj is the jth component of the unit normal vector on dA. The force coefficient is defined as

2

Borazjani et al. Present

Z A

Force coefficient

Force coefficient

In the above equation, x1 is the axial (flow) direction measured from the tip of the head, a(x1) is the wave amplitude that is assumed to vary nonlinearly along the fish body, k is the wave number of the undulations that corresponds to a wave length k, and x is the angular frequency. The following four similarity parameters must be specified in the fish-like swimming simulations: (1) the Reynolds number, which is based on the fish length L, the swimming speed U, and the kinetic viscosity m: Re = LU/m; (2) the Strouhal number, which is based on the maximum lateral excursion of the tail 2hmax and the tail beat frequency f: St = 2fhmax/U; (3) the non-dimensional wave length k; (4) the non-dimensional amplitude a(x1). In all simulations, the shape parameters, i.e. the wave length k and the parameters in the amplitude envelope a(x1), are kept the same. The non-dimensional wavelength k is taken to be 0.95. The wave amplitude is a quadratic function of x1

1

ð15Þ

-1 11.5

12

t/T

12.5

13

-1.5 11

11.5

12

t/T

Fig. 13. Time history of force coefficient for different Strouhal numbers at Re = 4000.

12.5

13

106

C.H. Zhou, C. Shu / Computers & Fluids 64 (2012) 98–107

0.8

simple way to calculate the efficiency. The Froude efficiency based on EBT for steady swimming is given as follows

0.6

Mean force coefficient

0.4

1 2

gEBT ¼ ð1 þ bÞ

0.2 0

where b = U/V is the slip velocity defined as the ratio of the swimming speed U to the wave speed V of undulatory wave along the body. An improved EBT efficiency formula (denoted herein as EBT-2) that takes into account the slope of the fish tail a = (k/ 2p)  [h0 (L)/h(L)] is given by Chen and Blickhan [28], as follows

-0.2 -0.4 -0.6 -0.8 -1

Borazjani et al. Re=300 Borazjani et al. Re=4000 Present Re=300 Present Re=4000

-1.2 -1.4 0

0.2

0.4

0.6

0.8

1

gEBT2

1.2

St Fig. 14. Variation of mean force coefficient with Strouhal numbers.

Z    pn1 dA þ  pn1 dA A A  Z Z   1  þ s1j nj dA þ  s1j nj dA 2 A A

1 2

TðtÞ ¼

DðtÞ ¼

Z

Z    pn1 dA þ  pn1 dA A A Z   Z   1 þ  s1j nj dA þ  s1j nj dA 2 A A

1 2

ð19Þ

Z

ð20Þ

The power loss due to lateral undulations of the fish body is calculated as follows

PL ¼

Z A

ð22Þ

_ ðpn2 þ s2j nj ÞhdA

ð21Þ

where h_ is the time derivative of the lateral displacement, i.e., the velocity of the lateral undulations. Lighthill’s elongated body theory (EBT) [27] is based on the assumption of slender body and two-dimensional inviscid flow. Although EBT overestimates the efficiency [17,28], it provides a

! 1 1 2 b2 ¼ ð1 þ bÞ  a 2 2 1þb

ð23Þ

We will compare the numerical results obtained by the present local DFD method with those of Borazjani and Sotiropoulos obtained by hybrid Cartesian/immersed boundary (HCIB) method [17]. In Fig. 13, the time history of the instantaneous hydrodynamics force coefficient CF for Re = 4000 and various St is shown. Positive and negative values indicate that the net force is of thrust- and drag- type, respectively. The values of force coefficient in Fig. 13 and the subsequently presented figures have been scaled with the force coefficient calculated for the rigid body (St = 0). It can be seen from Fig. 13 that the present computed results are smoother than those of Borazjani and Sotiropoulos. In Fig. 14, we show the variation of the time-averaged force coefficient C F with St for Re = 300 and 4000. From Figs. 13 and 14, it can be seen that the numerical results obtained by DFD and HCIB are in agreement. With the local refinement technique, the mesh used in our work has the same spacing near the body surface as that used by Borazjani and Sotiropoulos. However, because of the limitation of computational resources, our mesh is globally coarser than the mesh of Borazjani and Sotiropoulos (with 5.5  106 nodes). So, in Fig. 14, the disparity in the result for Re = 4000 is greater than that in the result for Re = 300 as the higher Reynolds number flow generally requires a finer mesh. In both results for Re = 300 and Re = 4000 shown in Fig. 14, the disparity is greater for lower Strouhal numbers. In [17], Borazjani and Sotiropoulos took the time step Dt = T/120. It may be too large for smaller St (smaller St corresponds to larger T).

Fig. 15. Instantaneous streamlines and vorticity contours at mid-plane. (A) Re = 4000, St = 0.2; (B) Re = 4000, St = 0.3; (C) Re = 4000, St = 0.7; (D) Re = 300, St = 0.3.

C.H. Zhou, C. Shu / Computers & Fluids 64 (2012) 98–107 Table 3 Froude efficiency at critical Strouhal number. Re

300 4000

107

Acknowledgments

Froude efficiency (%) EBT

EBT-2

CFD [17]

CFD (present)

59.3 67.25

59.5 66.67

18.86 22.95

16.28 26.93

This work was supported by National Science Foundation of China under Grant 11072113. We are grateful to F. Sotiropoulos and I. Borazjani for providing us with the geometry of the fish-like body. References

The simulated velocity and vorticity fields in mid-plane (x3 = 0) are shown in Fig. 15. The penetration of streamlines in the figure is due the fish motion. Different wake patterns depending on Re and St have been revealed. A single row regular Karman vortex street in Fig. 15A (Re = 4000 and St = 0.2) is of drag type. The wake pattern in Fig. 15B (Re = 4000 and St = 0.3) is also of single row type but more disorganized. It is still of drag type. The double row reverse Karman vortex street shown in Fig. 15C (Re = 4000 and St = 0.7) is of thrust type. The wake in Fig. 15D (Re = 300 and St = 0.3) is of drag type with a single row of vortices, and the pattern is more organized than that in Fig. 15B (higher Re case). In Fig. 15, the flow patterns near body are very similar to those obtained in the computation of Borazjani and Sotiropoulos [17]. There exist some slight differences in the far-field where our mesh is coarser and the vortices may be dissipated more. Finally, in Table 3, the propulsive efficiency values (at the critical Strouhal number St) obtained in this work are presented and compared with the CFD results of Borazjani and Sotiropoulos [17] and the results obtained by EBT and EBT-2. Obviously, the agreement between two CFD results is good. 6. Summary and conclusions In this paper, we extend the local DFD method to simulate three-dimensional viscous flows with geometrically complex boundaries moving with prescribed kinematics. The immersed boundary is treated as a sharp interface and the body surface is discretized with an unstructured triangular mesh. The discrete form of governing equations at an interior point may involve some points outside the solution domain. The functional values of flow variables at the exterior dependent points are updated at each time step by the approximate form of solution near the wall boundary, which is constructed through linear extrapolation along the normal direction to the boundary. Compared to the previous work, the tedious task of constructing new interpolation tetrahedrons is eliminated, and thus the complexity of DFD implementation is reduced greatly. An efficient algorithm for classifying mesh points is proposed to improve the overall efficiency of three-dimensional computations involving moving boundaries. The local DFD method for three-dimensional viscous flows was validated by its application to simulate flow past a stationary sphere at varied Reynolds numbers up to 300 as well as flow induced by a sphere that is rotating in an unbounded fluid at rest. For both cases, good agreement with referenced numerical solutions or experimental data was achieved. To demonstrate the applicability of the method to flows with complex, flexible immersed boundaries, we simulated the flow past an undulating fish-like body. The results were compared with those obtained by Borazjani and Sotiropoulos using HCIB method [17]. Both numerical solutions, for time history of force, time-averaged force, structure of wake patterns and propulsive efficiency, are in acceptable agreement even though the mesh used in this work is relatively coarser. The last test underlines the potential of present local DFD method as a powerful simulation tool for bio-fluid problems.

[1] Shu C, Fan LF. A new discretization method and its application to solve incompressible Navier–Stokes equations. Comput Mech 2001;27:292–301. [2] Shu C, Wu YL. Adaptive mesh refinement-enhanced local DFD method and its application to solve Navier–Stokes equations. Int J Numer Meth Fluids 2006;51:897–912. [3] Zhou CH, Shu C, Wu YZ. Extension of domain-free discretization method to simulate compressible flows over fixed and moving bodies. Int J Numer Meth Fluids 2007;53:175–99. [4] Zhou CH, Shu C. A local domain-free discretization method for simulation of incompressible flows over moving bodies. Int J Numer Meth Fluids 2011;66:162–82. [5] Zhou CH, Shu C. A local domain-free discretization method to simulate threedimensional compressible inviscid flows. Int J Numer Meth Fluids 2009;61:970–86. [6] Peskin CS. Flow patterns around heart valves: a numerical method. J Comput Phys 1972;10:220–52. [7] Peskin CY, McQueen DM. A three-dimensional computational method for flood flows in the heart: I. Immersed elastic fibers in a viscous incompressible fluid. J Comput Phys 1989;81:372–405. [8] Peskin CY, McQueen DM. Cardiac fluid dynamics. Crit Rev Biomed Eng 1992;20:451–9. [9] Goldstein D, Handler R, Sirovich L. Modeling a no-slip flow boundary with an external force field. J Comput Phys 1993;105:354–66. [10] Saiki EM, Biringen S. Numerical simulation of a cylinder in uniform flow: application of a virtual boundary method. J Comput Phys 1996;123:450–65. [11] Mohd-Yusof J. Combined immersed boundaries/B-splines methods for simulations of flows in complex geometries. CTR Annual Research Briefs, NASA Ames/Stanford University; 1997. p. 317–27. [12] Fadlun EA, Verzicco R, Orlandi P, Mohd-Yusof J. Combined immersedboundary/finite-difference methods for three-dimensional complex flow simulations. J Comput Phys 2000;161:35–60. [13] Gilmanov A, Sotiropoulos F. A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies. J Comput Phys 2005;207:457–92. [14] Marella S, Krishnan S, Liu H, Udaykumar HS. Sharp interface Cartesian grid method I: an easily implemented technique for 3D moving boundary computations. J Comput Phys 2005;210:1–31. [15] Udaykumar HS, Mittal R, Rampunggoon P, Khanna A. A sharp interface Cartesian grid method for simulating flows with complex moving boundaries. J Comput Phys 2001;174:345–80. [16] Mittal R, Dong H, Bozkurttas M, Najjar FM, Vargas A, von Loebbecke A. A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. J Comput Phys 2008;227:4825–52. [17] Borazjani I, Sotiropoulos F. Numerical investigation of the hydrodynamics of Carangiform swimming in the transitional and inertial flow regimes. J Exp Biol 2008;211:1541–58. [18] Choi Jung-I1, Oberoi Roshan C, Edwards Jack R, Rosati Jacky A. An immersed boundary method for complex incompressible flows. J Comput Phys 2007;224:757–84. [19] Belov A, Martnelli L, Jameson A. Three-dimensional unsteady incompressible flow calculations using multigrid. AIAA-1997-0443; 1997. [20] Huang CW, Shih TY. On the complexity of point-in-polygon algorithms. Comput Geosci 1997;23:109–18. [21] Johnson TA, Patel VC. Flow past a sphere up to a Reynolds number of 300. J Fluid Mech 1999;378:19–70. [22] Clift R, Grace JR, Weber ME. Bubbles, drops and particles. New York: Academic Press; 1978. [23] Kim J, Kim D, Choi H. An immersed boundary finite volume method for simulations of flow in complex geometries. J Comput Phys 2001;171:132–50. [24] Dennis SCR, Singh SN, Ingham DB. The steady flow due to a rotating sphere at low and moderate Reynolds numbers. J Fluid Mech 1980;101:257–79. [25] Sotiropoulos F, Borazjani I. Private, communication; 2009. [26] Tytell ED, Lauder GV. The hydrodynamics of eel swimming, 1. Wake structure. J Exp Biol 2004;207:1825–41. [27] Lighthill MJ. Aquatic animal propulsion of high hydromechanical efficiency. J Fluid Mech 1970;44:265–301. [28] Chen JY, Blickhan R. Note on the calculation of propeller efficiency using elongated body theory. J Exp Biol 1994;192:169–77.