A SIMPLE MOMENTUM FORCING METHOD ON UNSTRUCTURED
CARTESIAN MESH FOR IMMERSED BODIES IN INCOMPRESSIBLE
FLOWS
Dartzi Pan and Tzung-Tza Shen Department of Aeronautics and Astronautics
National Cheng Kung University Tainan, Taiwan, ROC e-mail: dpan@mail.ncku.edu.tw
Key words: Immersed Boundary Method, Direct Momentum Forcing, Unstructured Cartesian
Mesh
Abstract. The incompressible Navier-Stokes equations are solved by an implicit pressure
1. INTRODUCTION
Recently numerical methods for solving incompressible flows on fixed Cartesian grids are gaining popularity for their relative ease in treating complex immersed bodies [1-12]. The differences among these methods lie on the different ways the boundary conditions for the immersed bodied are enforced. The immersed boundary method [1] and the virtual boundary method [2] simulate the no-slip boundary condition on the body surfaces by adding appropriate momentum forcing terms to the appropriate cells. The force added to the flow field is spread over several grid cells using a delta function-like distribution normal to the fluid-solid interface. In the direct-forcing method [4, 5] the no-slip condition at the fluid-solid interface is used to interpolate the velocity distribution for cells around the solid body without actually solving the fluid equations. The fictitious domain method [6] is a finite element method that enforces the fluid-body interface boundary condition in a weak form using Lagrange multipliers. The Lagrange multiplier can be viewed as a momentum-forcing term similar to the term added in the immersed boundary method. In Cartesian cut-cell method [7], the flow variables for the cut cells are solved based on the actual shape of the merged cut cells. In Xiao [8], the solid body is viewed as being made of a material distinct from the surrounding fluid. The motion of the body is identified and tracked by a color function. In Kajishima and Takiguchi [9], the volume fraction of solid is used to couple the fluid velocity and the solid velocity. This coupling is written as a forcing term in the momentum equation, which is implemented by a separate step in a fractional step method.
Very recently, Pan [10,11] viewed the domain inside the solid body as being occupied by the same fluid as outside with a prescribed divergence-free velocity field. In this view a fluid-body interface is similar to a fluid-fluid interface commonly encountered in the Volume of Fluid (VOF) method for the two-fluid flow problems. Thus a Volume of Body (VOB) function analogous to the VOF function can be used to identify and track the presence of the immersed body. For the grid cells containing the fluid-body interface, a mixture of the “two fluids” based on the fluid volume is assumed. This volume averaging of the velocity automatically enforces the no-slip boundary condition inside the interface cell, but it also smears the fluid-body interface to the size of one cell width. In this work, we tried to improve the representation of the fluid-body interface by using a scaled momentum forcing term. Similar to the VOB method, a direct momentum forcing is added to cells inside the body to enforce the no-slip boundary condition. This momentum forcing is further assumed to have an effective range of one cell size normal to the body surface, acting on cells external to the immersed body. For those cells within the active range, the amount of added momentum forcing is linearly scaled by its normal distance from the cell center to the body surface. Test results of this simple scaled momentum forcing (SMF) method are reported for various steady and unsteady flows over a circular cylinder. The flow over an impulsively started circular cylinder and the flow over an ellipse in Figure-8 motion are also computed to verify the capability of the current method to treating moving bodies.
The incompressible Navier-Stoked equations are 0 1 ) ( 0 2 +∇ = ∇ − ⋅ ∇ + ∂ ∂ = ⋅ ∇ P v Re v v t v v v v v v v v v v (1)
where vv and P are Cartesian velocity and pressure, Re is Reynolds number. The backward
time differencing scheme is used to advance the momentum equations with the pressure fixed at time level n : 0 * 1 3 2 * 1 + +∇ = Δ + − − n n n P R t v c v c v c v v v v v (2)
where Δt is the time increment andRv* represents the convection and the viscous fluxes. The constants are c1=1.5, c2=2 and c =0.5 for the second-order accurate backward differencing 3
scheme, and c1=1, c2=1 and c =0 for the first order Euler implicit scheme. The intermediate 3
velocity vv generally does not satisfy the divergence-free condition. It is corrected by the ∗
following correction step:
φ ∇ Δ − = ∗ + v v v v t vn 1 (3)
where φ =Pn+1−Pn is the pressure correction. By requiring vvn+1 be divergence-free, we obtain the Poisson equation for the pressure correction:
t v Δ ⋅ ∇ = ∇2 v* v φ (4)
Equations (2), (3) and (4) constitute the implicit fractional step pressure correction method used in this work.
3. FINITE VOLUME DISCRETIZATION
Equation (2) is the second-order backward differencing scheme in time. The technique of sub-iteration is employed on the right hand side, and an implicit multi-grid method is developed to accelerate the convergence of the inversion process. Also, a similar implicit multi-grid method is developed to solve Eq. (4), or the pressure Poisson equation.
4. SCALED MOMENTUM FORCING
The position and the velocity vvB of the immersed bodies are assumed known from some
appropriate governing equations. The body velocity vvB is assumed to be divergence-free, such that the pressure field inside the body surface obeys the same governing equation as the pressure field outside. Under this assumption, Eq. (4) can be applied to the entire computational domain including the body interior.
The current immersed boundary method simulates the existence of solid body by adding a direct momentum forcing term to the momentum equation:
0 * 1 3 2 * 1 + +∇ + = Δ + − − IBM n n n f P R t v c v c v c v v v v v v (5) The modeling of fvIBM constitutes the major difference among various immersed boundary
methods. In this work, it is assumed that fvIBM depends on the normal distance r of a cell n
center to the body surface, where r > 0 indicates the body exterior, n r < 0 indicates the body n
interior and r = 0 indicates the body surface. n
For cell centers inside the body surface with r < 0, the working of n fvIBM should enforce the no-slip boundary condition, resulting the following equation for each body cell
0 1 1 * 1 = Δ − + t v c v c v vBn (6) where vvBn+1 is the known body velocity at the new time level n+1. This suggests that the direct
forcing added to body cells is
) ( * 1 3 2 1 1 n n n n B IBM R P t v c v c v c f + +∇ Δ + − − = v + v v − v v v (7) This fvIBM is nothing but the discretized Navier-Stokes equations with a negative sign and a
modified acceleration term.
) ( * 1 3 2 1 1 , n p p n p n p n Bp SMF p IBM R P t v c v c v c f + +∇ Δ + − − = − + v v v v v v φ (8) where the subscript “p” indicates the projection point of the cell center under consideration on
the body surface and φSMF is defined as
⎩ ⎨ ⎧ > ≤ ≤ − = 1 / , 0 1 / 0 , / . 1 η η η φ n n n SMF r r r (9)
where η is a preset effective length of fvIBM,p and is chosen to be the minimum cell size lmin of the mesh. Not that the condition of φSMF =1 when rn =0 will enforce the no-slip
boundary condition on the body surface, and the condition of φSMF =0 when rn ≥η indicates
that no forcing is added when the cell is one cell size away from the body surface. Generally, the projection point does not coincide with the cell center and the terms inside the brackets of Eq. (8) need to be estimated using nearby cell center values. In this work, we choose a bilinear interpolation based on the cell center values surrounding the projection point. By the above definition, φSMF =0.5 is a closed contour external to the body surface at a normal distance of 0.5lmin. This is in the middle of the transition zone between the solid body and
flow outside.
Note that here fIBM,p v
is similar to the Lagrangian force developed in the Physical Virtual Model (PVM) of Silva et al. [12]. In PVM method, the Lagrangian force acting on the fluid-structure interface is estimated by some special form of interpolation and then distributed to nearby cells through a discretized Delta function.
6. FORCE CALCULATION
To compute the aerodynamic forces acting on the immersed body, the surface integral of pressure and viscous stress is needed:
) 1 ( v S Re S P f surface Body v v r v v Δ ⋅ ∇ + Δ − =
∑
(9)where the summation is performed on all surface elements of the body surface. The pressure at each surface element center is obtained by a bilinear interpolation based on the vertex values of the cell containing the surface element center. The velocity gradient for a surface element, however, is not estimated at the element center in this work, but at a point having a normal distance 0.5lmin away from the surface element. In other words, the velocity gradient
7. FLOWS OVER A CIRCULAR CYLINDER
The steady and unsteady flows over a circular cylinder of unit diameter at Re=20, 40, 100 and 200 are computed on an unstructured Cartesian grid. The grid shown in Fig. 1(a) is refined around the cylinder to have about 384 cells containing the cylinder surface. The cylinder volume computed by the rn =0 contour is about 0.1% less than the true value. The outer boundaries are 20 diameters away from the cylinder. The uniform flow condition is set to the inflow and the two side boundaries. The downstream boundary follows the upwind differenced equation of (∂vv/∂t)+Un(∂vv/∂x)=0, where U is the normal outflow velocity n
at the boundary.
For the steady case of Re=40, the Euler implicit method is used with the maximum CFL number around 20. The infinity norm of the steady state residual for φSMF =0 cells, or the terms in Eq. (2) without involving Δt, dropped 4 orders of magnitude in 500 steps. The computed streamlines are plotted in Fig. 1(b). For convenience, the contour of rn =0 is displayed as the cylinder wall. The flow is steady with a separation bubble behind the cylinder. The streamlines around the cylinder are smooth, indicating the effects of the interface cells in smoothening the zigzag representation of the cylinder surface. Figure 2(a) shows the computed pressure contours. Note that the pressure inside the cylinder adjusts itself automatically to the pressure field outside. The pressure contours intersect the cylinder wall in a nearly orthogonal manner. The pressure coefficient at the intersection of the cylinder surface and the grid lines is interpolated using the surrounding center values and plotted in Fig. 2(b). The data from Fornberg [13] are also included for comparison. The two results generally agree with each other very well. The current computation predicted a slightly lower surface pressure distribution on the leeward surface. Note that the computed surface pressure has only very minor wiggles around the shoulder region. It would show zigzag oscillation if the fourth order dissipation term was turned off.
For the unsteady case of Re=200, the second order backward difference scheme is used for time integration. The time increment is chosen such that the expected vortex shedding cycle takes about 50 time steps to complete. The instantaneous streamlines at certain instant in the periodic vortex shedding process are plotted in Fig. 3(a). The vorticity contours are shown in Fig. 3(b). The unsteady vortex shedding behind the cylinder is clearly seen.
Table I lists the computed lift coefficient (Cl), drag coefficient (Cd), the wake length (Lw) normalized by the diameter (d), the Strauhal number (St) of unsteady vortex shedding and some published results by others. The comparisons between the present work and the work of others are generally satisfactory.
8. IMPULSIVELY-STARTED CYLINDER
moved to the left at a constant speed. The Reynolds number based on the cylinder diameter D and cylinder speed Ucyl is 550. The mesh is clustered along the path of the cylinder movement such that there are always about 312 interface cells around the cylinder. During the computation of 600 time steps, the cylinder was moved a distance of 3 diameters from its initial position. Figure 4(a) shows the time dependent development of velocity along the symmetry axis in the wake. In this figure the velocity u was measured relative to the moving cylinder and normalized by the cylinder speed Ucyl. The distance x/D along the symmetry
axis was referenced relative to the cylinder center. The symbols in Fig. 4(a) are measured manually from the experimental points in Bouard and Coutanceau [15], while the lines are from our computation. Different symbols correspond to different time ts =tUcyl /D with an increment of 0.5, starting from ts =0.5. The agreement between computation and experiment is generally satisfactory. Figure 4(b) shows the computed instantaneous streamlines at
0 . 3 =
s
t . A pair of isolated secondary vortex between the shoulder and the rear stagnation can be clearly seen. This is in excellent agreement with the experimental observation of Bouard and Coutanceau [15]. This example basically has validated the capability of the present method to treat moving bodies in the flow field.
9. FLOW OVER A FLAPPING ELLIPSE
In Fig. 5(a), an ellipse is in a forced “figure-8” motion. The down-stroke phase is plotted in black while the up-stroke phase is in red. As described in Wang [16], a “figure-8” motion is used to simulate the flapping motion of an insect wing. The movement of the center of the ellipse is governed by A(t)=0.5A0[cos(2πt/T)], where T is the flapping period andA is the 0
amplitude. The stoke plane along which the center moves has an inclination angle β with respect to the horizontal axis. The angle of attack is governed by 0.25π[1−sin(2πt/T+ϕ)], where ϕ is a phase lag. In our choice, the reference speed is 400cm/s, the reference length is 4cm, the major axis of the ellipse is 1cm, and the kinematic viscosity is 2cm2 / . The s
flapping motion has β =60o, A0=2.5cm, T=0.025sec and ϕ =0. The Reynolds number
amounts to Re=800. The time increment is set to 1/400 T. The outer boundary of the grid is 40cm in length. The grid is refined in the region where the figure-8 motion occurs. The computed volume of the ellipse is about 0.25% lower than the analytical value.
periods. Comparing with the Fig. 3 of Wang [16], although many differences are observed, but the magnitude of the aerodynamic forces and the general trend of variation are comparable to the reference.
9. CONCLUSIONS
A simple immersed boundary method on unstructured Cartesian meshes is developed and tested for incompressible flows with immersed bodies. The solid body is modeled by a direct forcing term added to the momentum equations. For cell centers inside the solid body, this forcing term enforces the no-slip boundary condition directly. For cells external to the body within an effective range, this forcing term is linearly scaled by the normal distance from the cell center to the body surface. The same pressure Poisson equation applies to the entire computational domain including the body interior. An implicit pressure correction method and multi-grid methods are used to integrate the governing equations. Steady and unsteady flows over a stationary circular cylinder with Re=20, 40, 100 and 200 are computed to validate the present method. The comparisons of the computed results with published data are generally satisfactory. The flow over an impulsively started cylinder and the flow over an ellipse in Figure-8 motion are computed to verify the capability to treat moving bodies in the flow field.
REFERENCES
[1] Peskin C. S., The Immersed Boundary Method, Acta Numerica, 2002; 1-39.
[2] Goldstein D., Handler R. and Sirovich L., Modeling a No-Slip flow Boundary with an External Force Field, J. Comput. Phys., 1993; 105, 354-366.
[3] Saiki E. M. and Biringen S., Numerical Simulation of a Cylinder in Uniform Flow: Application of a Virtual Boundary Method, J. Comput. Phys, 1996; 123, 450-465.
[4] Fadlun E.A., Verzicco R., Orlandi P. and Mohd-Yusof J., Combined Immersed-Boundary Finite-Difference Methods for Three-Dimensional Complex Flow Simulations, J. Comput. Phys, 2000; 161, 35-60.
[5] Kim J., Kim D. and Choi H., An Immersed-Boundary Finite-Volume Method for Simulations of Flow in Complex Geometries, J. Comput. Phys., 2001; 171, 132-150. [6] Baaijens Frank P.T., A Fictitious Domain/Mortar Element Method for Fluid-Structure
Interaction, Int. J. Numer. Meth. Fluids, 2001; 35, 743-761.
[7] Ye T., Mittal R., Udaykumar H.S. and Shyy W., An Accurate Cartesian Grid Method for Viscous Incompressible Flows with Complex Immersed Boundaries, J. Comput.
Phys., 1999, 156, 209-240.
[8] Xiao F., A Computatioinal Model for Suspended Large Rigid Bodies in 3D Unsteady Viscous Flows, J. Compu. Phys., 1999, 155, 348-379.
[9] Kajishima T. and Takiguchi S., Interaction Between Particle Clusters and Particle-Induced Turbulence, Int. J. Heat and Fluid Flow, 2002, 23, 693-646.
[11] D. Pan, “An Immersed Boundary Method on Unstructured Cartesian Meshes for Incompressible Flows with Heat Transfer”, Numerical Heat Transfer, Part B,
Fundamentals, accepted and in press, 2005.
[12] A.L.F. Lima E Silva, A. Silverira-Neto, and J.J.R. Damasceno, “Numerical Simulation of Two-dimensional Flows over a Circular Cylinder Using the Immersed Boundary Method,” Journal of Computational Physics, Vol. 89, pp. 351-370, 2003.
[13] Fornberg B.,A Numerical Study of Steady viscous Flow Past a Circular Cylinder, J.
fluid Mech., 1980, 98, 819-855.
[14] Kiris C. and Kwak D., Numerical Solution of Incompressible Navier-Stokes Equations Using a Fractional-Step Approach, Comp. & Flu., 2001, 30, 829-851.
[15] Bouard, R. and Coutanceau, M. The Early Stage of Development of the Wake Behind An Impulsively Started Cylinder for 40<Re<104, J. Fluid Mech., 1980, 101, pp. 583-607.
[16] Z. Jane Wang, “Two Dimensional Mechanism of Hovering,” Phys, Rev. Lett, Vol. 85, pp. 2216-2219, 2000. Methods Re Cd Lw/d Cl St Current method 20 2.03 0.90 Fornberg (1980) 20 2.00 0.91 Ye et al. (1999) 20 2.03 0.92 Current method 40 1.51 2.25 Fornberg (1980) 40 1.50 2.24 Ye et al.(1999) 40 1.52 2.27 Current method 100 1.32 ±0.32 0.164 Kim, Kim and Choi (2001) 100 1.33 ±0.32 0.165 Current method 200 1.30±0.04 ±0.61 0.196 Kiris and Kwak (2001) 200 1.27±0.04 ±0.67 0.184
-1 0 1 x -1 0 1 y
Grid For A Stationary Cylinder
x y -1 0 1 2 3 -2 -1 0 1 2 Streamlines, Re=40 Lw/D = 2.245 (a) (b)
Fig. 1 Flows over a stationary circular cylinder, Re=40, (a) grid and (b) streamlines.
x y -1 0 1 2 3 -2 -1 0 1 2 Cp, Re=40 Theta Cp 0 50 100 150 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Fornberg[1980] computed Re=40, Surface Cp (a) (b)
x y -1 0 1 2 3 -2 -1 0 1 2 Streamlines, Re=200 x y 0 5 10 -5 0 5 Re = 200, Stationary Cylinder (a) (b)
Fig. 3 Flow over a stationary cylinder, Re=200, (a) instantaneous streamlines(b) vorticity contours. X/D U /U cyl 0.5 1 1.5 2 2.5 -1.5 -1 -0.5 0 0.5 1
Suddenly Started Cylinder, Re=550
x y 0.3 0.4 0.5 0.6 0.7 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 Instantaneous Streamlines Re = 550, Ts = 3
Fig. 4 Flow over an impulsively started cylinder,
(a) Velocity u along symmetry axis in the wake of an impulsively started cylinder, 550
/ν =
D
Ucyl , Symbols: experimental data of Bouard and Coutanceau (1980) at different time ts =tUcyl /D with an increment of 0.5, starting from ts =0.5 (solid circle), Lines: computed,
(a) (b)
(c) (d)
(e) (f)
