What is simple immersed boundary method?

denotes the entire fluid domain . Discretization of these equations can be done by assuming an Eulerian grid on the fluid and a separate Lagrangian grid on the fiber. Approximations of the Delta distribution by smoother functions will allow us to interpolate between the two grids.

What is the immersed boundary method for computational fluid dynamics?

The immersed boundary method is a numerical method in computational fluid dynamics where the flow boundary is immersed in the grid that does not conform with the boundary . In the immersed boundary method, special treatment has to be taken at the boundary to incorporate the boundary conditions.

What are the advantages of immersed boundary method?

The Immersed Boundary Method (IBM) has an advantage in simulating fluid–structure interaction, owning to its simplicity, intuitiveness, and ease of handling complex object boundaries . The interpolation function plays a vital role in IBM and it is usually computationally intensive.

What is the immersed boundary force?

The immersed boundary (IB) force density is directly computed by the difference between the unforced fluid velocity and the particles' boundary velocity on the Lagrangian points .

What is a simple boundary condition?

Boundary conditions are constraints necessary for the solution of a boundary value problem . A boundary value problem is a differential equation (or system of differential equations) to be solved in a domain on whose boundary a set of conditions is known.

What is the immersed body method?

Fully resolved simulation techniques for immersed bodies in fluids are those in which fluid flow around each individual body is fully resolved without the use of approximate governing equations or models such as drag laws . The immersed objects can be rigid, elastic, or self-propelling bodies.

What is immersed boundary method software?

Immersed Boundary Method Adaptive Mesh Refinement Software Infrastructure . IBAMR is a distributed-memory parallel implementation of the immersed boundary (IB) method with support for Cartesian grid adaptive mesh refinement (AMR). Support for distributed-memory parallelism is via MPI, the Message Passing Interface.

What is immersed boundary method in foam theory implementation and use?

In this sense, the Immersed Boundary Method is a useful and efficient approach for cheap mesh generation in Computational Fluid Dynamics : it leverages on non- conformal boundary surfaces through a Cartesian mesh, mimicking the presence of the body through the imposition of boundary forces, by including additional 11 ...

(PDF) A simple and efficient direct forcing immersed boundary framework for fluid–structure interactions

Journal of Computational Physics 231 (2012) 5029–5061

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics

journal homepage: www.elsevier .com/locate / jcp

A simple and efficient direct forcing immersed boundary frameworkfor fluid–structure interactions

Jianming Yang ⇑, Frederick SternIIHR – Hydroscience and Engineering, University of Iowa, Iowa City, IA 52242, USA

a r t i c l e i n f o

Article history:Received 29 August 2011Received in revised form 3 April 2012Accepted 4 April 2012Available online 26 April 2012

Keywords:Fluid–structure interactionImmersed boundary methodDirect forcingStrong couplingPredictor–corrector algorithmField extensionMoving boundaryVortex-induced vibrationGallopingFlutteringTumblingVortex shedding

0021-9991/$ - see front matter � 2012 Elsevier Inchttp://dx.doi.org/10.1016/j.jcp.2012.04.012

⇑ Corresponding author. Tel.: +1 319 335 5749; faE-mail address: [emailprotected] (J. Ya

a b s t r a c t

A direct forcing immersed boundary framework is presented for the simple and efficientsimulation of strongly coupled fluid–structure interactions. The immersed boundarymethod developed by Yang and Balaras [J. Yang, E. Balaras, An embedded-boundary formu-lation for large-eddy simulation of turbulent flows interacting with moving boundaries, J.Comput. Phys. 215 (1) (2006) 12–40] is greatly simplified by eliminating several compli-cated geometric procedures without sacrificing the overall accuracy. The fluid–structurecoupling scheme of Yang et al. [J. Yang, S. Preidikman, E. Balaras, A strongly-coupled,embedded-boundary method for fluid–structure interactions of elastically mounted rigidbodies, J. Fluids Struct. 24 (2008) 167–182] is also significantly expedited by moving thefluid solver out of the predictor–corrector iterative loop without altering the strong cou-pling property. Central to these improvements are the reformulation of the field extensionstrategy and the evaluation of fluid force and moment exerted on the immersed bodies, bytaking advantage of the direct forcing idea in a fractional-step method. Several cases withprescribed motions are examined first to validate the simplified field extension approach.Then, a variety of strongly coupled fluid–structure interaction problems, including vortex-induced vibrations of a circular cylinder, transverse and rotational galloping of rectangularbodies, and fluttering and tumbling of rectangular plates, are computed. The excellentagreement between the present results and the reference data from experiments and othersimulations demonstrates the accuracy, simplicity, and efficiency of the new method andits applicability in a wide range of complicated fluid–structure interaction problems.

1. Introduction

The immersed boundary method was originated by Peskin in the last seventies [20,21] for the computational studies offluid–structure interaction (FSI) problems in a human heart. This method and its subsequent developments/extensions havesubstantially expanded the applicability of the traditional finite difference methods relying upon simple Cartesian grids. InPeskin’s method [20,21], the effect of flexible structures, such as heart valves and muscular wall, on the blood flow was rep-resented by an external, singular force field in the Navier–Stokes equations to be solved in a regular domain including boththe flow field and the structures; and these structures were modeled by sets of spring-linked elements using a Lagrangianrepresentation. Discrete delta functions were used to connect the fluid flow and the immersed boundaries, i.e., to interpolatethe fluid velocity from the Eulerian grid points to the Lagrangian elements and to spread out the boundary forces from thelatter to the former. For flexible structures, the boundary forces were readily available using a generalized Hooke’s law fromthe relative positions of the boundary elements depending on the boundary material properties. However, for a solid bound-

x: +1 319 335 5238.ng).

http://dx.doi.org/10.1016/j.jcp.2012.04.012

mailto:[emailprotected]

http://dx.doi.org/10.1016/j.jcp.2012.04.012

http://www.sciencedirect.com/science/journal/00219991

http://www.elsevier.com/locate/jcp

5030 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

ary, there is no relative displacement between neighboring boundary elements and the above method could not be applied ina straightforward manner. Goldstein et al. [11] developed a feedback forcing strategy in a pseudo-spectral method usingideas from control systems theory, in which the local forcing function would adapt itself through a feedback controller tosatisfy the desired velocity boundary conditions. Saiki and Biringen [23] implemented this feedback forcing scheme in afourth-order finite difference method to study low Reynolds number flows past stationary and moving circular cylinders.They used a discrete hat function (area-weighted averaging) to transfer velocity/force information between the backgroundgrid and the immersed boundary. A major disadvantage of the feedback forcing method is the overly restrictive time stepconstraint resulting from the feedback mechanism, which contains two large-magnitude tuneable constants that makethe equations very stiff. This limitation was removed by Mohd-Yusof [18], who derived a direct forcing formulation byimposing the velocity boundary condition at the immersed boundary ‘‘exactly’’ in the discrete-time equations withoutany feedback adjustment. The resulting forcing term is no longer defined in the continuous space and associated with theboundary force density function; hence there is no need to transfer force information from boundary elements to Euleriangrid points and the immersed boundary discretization is no longer required. This approach was implemented in a pseudo-spectral method in [18] and then applied in a finite difference method by Fadlun et al. [8]. With Mohd-Yusof’s derivation, it isessentially a local solution reconstruction procedure to satisfy the desired boundary conditions at the immersed boundaryand an explicit forcing term was not required in the momentum equations at all, as detailed in [8]. On the other hand, thedirect forcing idea was also combined into some approaches [26,37,27,36] immediately related to Peskin’s method, to whicha discrete delta function is central for the information transfer between the Eulerian grid points and the Lagrangian elements.In these approaches, the velocities at surrounding grid points are first interpolated using the discrete delta function to animmersed boundary node; then, instead of using a constitutive law or feedback adjustment, the local forcing is determinedby directly requiring the velocity boundary condition at the boundary node to be satisfied after the forcing is applied; finallythe local forcing at each boundary node is distributed using the discrete delta function to surrounding grid points. In thepresent study, we shall focus on a direct forcing approach making no use of the discrete delta functions. Nonetheless, theextension of our FSI algorithm to those with discrete delta functions is straightforward and requires almost no modificationsto either the original immersed boundary solver or our algorithm, as will be clear later.

The direct forcing immersed boundary approaches have been well received among the developers and practitioners ofnon-boundary conforming methods in the computational fluid dynamics (CFD) field, mainly because of the simplicity ofthe concept and the ease of the implementation of the formulations. Initially, the developments were focused on stationaryimmersed boundaries [8,15,25,2]; and very few were applied to problems with moving boundaries, due to the fact that theimplications of boundary movement on a fixed grid in a time-splitting scheme, such as the fractional-step method, were notsystematically addressed. (Note that this problem is not prominent in approaches using discrete delta functions as the forc-ing is spread out over a few grid points across the immersed boundary.) Then Yang and Balaras [31] demonstrated that, non-physical historical information may enter the flow field in a time step, when some grid points with reconstructed solutions atthe previous time step become normal fluid points, if no treatment is applied to recover the correct historical information atthese points. They proposed a field extension strategy that, at the end of each time step, the flow field is extended into thegrid points with non-physical values near the immersed boundary through extrapolations. A variety of examples rangingfrom laminar to turbulent flows were given to show the accuracy and effectiveness of this approach. They further extendedit in [32] and then [33] to FSI problems with multiple rigid bodies using a strong coupling scheme in which the structuraldynamics was solved via Hamming’s fourth-order predictor–corrector algorithm. In their strong coupling scheme, the fluidand the structure are treated as elements of a single dynamic system, and both sets of governing equations are integratedsimultaneously and interactively in the time domain. It is a very efficient iterative scheme as the number of iterations doesnot change much with increased number of DOFs (degrees of freedom) of the structural part and the convergence of thewhole coupled system usually is reached within ten iterations. In addition, it is not limited to FSI problems with solid bodies,for example, in [28], it was used to study the aerodynamic performance of a flexible hovering wing.

In general, the approaches discussed above for FSI problems (prescribed motion in [31] and predicted motion in [33],respectively) are quite straightforward and efficient. One issue with the field extension strategy is that the definition of pres-sure points requiring the extension operation is not as simple as that for velocity points. Instead of a simple geometric rela-tionship, for instance, the closest grid points to the immersed boundary along grid lines (for velocity components, closestgrid points in the fluid phase and solid phase are defined as forcing points and pseudo forcing points, respectively, in[31]), the status of all surrounding velocity points (four points in two-dimensional (2D) and six points in three-dimensional(3D) cases, respectively) has to be used to find a pressure point that needs the extension operation. And such a point may bein either the fluid or the solid phase, depending on the configuration of the immersed boundary and the grid layout. This isnot particularly convenient in terms of algorithm design and implementation. Furthermore, for grid points in the solidrequiring field extension operation, ambiguities may sometimes exist near sharp corners or under-resolved regions with re-gard to the normal directions to the immersed boundary. To remove these ambiguities, significant amount of coding work isneeded and the clarity and efficiency of the algorithm may be affected.

On the other hand, in [31,33], the fluid force on an immersed boundary was evaluated through a surface force integrationprocedure. Basically, the immersed boundary is first discretized into elements of size similar to the local grid spacing; thenthe pressure value and velocity derivatives at the surface are obtained through extrapolation and one-sided differencing,respectively; with the stress tensor and the geometric information (normal and area) available for each boundary element,the surface force distribution and total force and moment can be evaluated directly. The overall procedure is very generalized

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5031

and not limited to a specific variant of the immersed boundary methods. However, the results from this approach depend onthe resolution of the surface discretization and the locations of the probe points used to obtain flow field information. Also, itis a post-processing step that can only be performed after the whole flow field has been solved. Actually, Lai and Peskin [16]gave several different approaches for evaluating fluid force on an immersed body, which are applicable to different types ofimmersed boundary methods. The most straightforward approach is to integrate the momentum forcing function over thewhole domain, which represents the total effect of immersed boundary on the fluid. Since the forcing function is availablebefore the final solution of the momentum equations, with this pointwise integration approach, in contrast to the surfaceintegration approach, the fluid force can be evaluated far before the whole field is solved. Unfortunately, this approach can-not be directly applied to methods given in [31,33] as the local solution reconstruction and field extension procedures wereapplied to the predicted and final solution fields, respectively. An explicit forcing function was not used therein and indeed itwas cumbersome to define because of the latter procedure.

Another concern with the strong coupling FSI algorithm in [33] is the computational cost. Weak coupling schemes, usu-ally with a time lag of the fluid force/moment on the solid, are very cost-effective for many applications, but they may be-come unstable in numerous cases with very light structures of similar or even lower density compared to that of thesurrounding fluid. For example, in [33], a weak coupling scheme was also implemented and compared with the strong cou-pling scheme using an elastically mounted circular cylinder with two DOFs in a free-stream. The results showed bothschemes gave similar results for the case with a larger mass ratio (solid to liquid). However, for the case with a mass ratiocloser to unity, the weak coupling scheme diverged very soon due to an exponential instability; whereas the strong couplingscheme was still stable, even for cases with lower mass ratios. This desired property of numerical stability comes at the priceof a much higher computational cost due to the interactive manner of solving the coupled systems. In general, the mostexpensive part for simulating FSI problems with an incompressible flow is the Poisson equation solver. In some cases, fastPoisson solvers utilizing FFT (fast Fourier transform) and cyclic reduction can be very useful to reduce the computationaltime. But for complicated cases, such as two-phase flows and adaptive mesh refinement, fast Poisson solvers are not avail-able or straightforwardly applicable and the computational cost from using a strong coupling scheme can be much higherthan that from using a weak coupling one. Interestingly, Kim and Choi [14] developed a strong coupling scheme with theirdirect forcing immersed boundary method without the coupling iterations by using a non-inertial reference frame fixed tothe body. This scheme can be very efficient for some applications, but it is limited to a single solid body in an unconfineddomain and the use of a non-inertial reference frame makes the boundary conditions, especially, inflow and outflow, difficultto handle.

In the present study, we shall present a new formulation that can address the three issues discussed above in a simple andelegant manner. First, we shall reformulate the pressure field extension: instead of the pressure field itself, the field exten-sion operation will be on the three (two for 2D cases) pressure gradient components, as the non-physical historical informa-tion for pressure points enters the system through the pressure gradient term in the momentum equations. On a staggeredgrid, the pressure gradient components are collocated with the velocity components. Therefore, it is possible to utilize theknown stencils for the local reconstruction of velocity components at forcing points, without going through the cumbersomeprocess of identification and extrapolation for the pressure extension points. Then, we shall describe a new strategy consis-tent with our field extension operation for obtaining the momentum forcing term explicitly. Therefore, we can use thestraightforward point integration of the momentum forcing term to evaluate fluid force and moment exerted on the im-mersed bodies instead of the surface integration approach in [31,33]. Finally, by taking advantage of the direct forcing ideain the context of a fractional-step method, the fluid solver can be moved out of the FSI coupling iterative loop without alter-ing the strong coupling property of the overall FSI algorithm. This represents an acceleration of several times or even oneorder of magnitude to the whole computation using the original method in [33].

2. Mathematical formulation

The computational domain for a FSI problem consists of the fluid and structural domains. In this work, a primitive var-iable method is used for the fluid domain, i.e., the velocity u and pressure p are the unknown variables of the fluid flow field.For the structural part, the displacement is chosen as the unknown field variable and, for rigid body motion considered in thepresent study, it can be represented by the linear displacement and orientation of a rigid body. It is obvious that the struc-tural displacement defines the interface of the fluid and structural domain. In an immersed boundary approach, usually aCartesian grid is used to cover the whole domain, which means the structural boundary or the fluid–structure interface,C, has to be given (and tracked for moving boundary problems) on the fixed Cartesian grid as well. On the other hand, sincethe focus of this work is on the incompressible flow interacting with moving bodies, the introduction of body-fixed non-iner-tial reference frames facilitates the solution of rigid body dynamics naturally. Therefore, the reference frames used in thepresent study are discussed first.

2.1. Reference frames

As shown in Fig. 1, two types of reference frames are usually used in FSI problems: an inertial reference frame fixed to ormoving at a constant velocity with respect to the earth, and one or multiple non-inertial reference frames with each attached

Fig. 1. Relationship between the inertial and non-inertial reference frames.

5032 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

to one moving body under consideration. All fluid variables are conveniently defined in the inertial reference frame, whereasfor a moving body, the motion can be described by the linear and angular displacements with respect to the inertial referenceframe, and the linear and angular velocity with respect to the its own body-fixed non-inertial reference frame.

For example, in the inertial reference frame, the translational displacement of body a is given by

da ¼ ðxOa ; yOa; zOa Þ

T; ð1Þ

where the superscript T denotes the transpose of a vector or a matrix. It is the position vector from the origin of the inertialreference frame to the origin of the non-inertial reference frame attached to body a. The latter is chosen to coincide with thecenter of mass of body a. And

Ha ¼ ð/a; ha;waÞT ð2Þ

gives the three Euler angles for defining the orientation of body a. For simplicity, subscript referring to a body will be omittedhereafter. On the other hand, the linear and angular velocities of the body are defined as

v ¼ ðvx;vy;vzÞT; ð3ÞX ¼ ðXx;Xy;XzÞT; ð4Þ

with respect to the non-inertial reference frame, respectively.The Euler angle H can be used to define a rotation matrix R, which can be written as:

R ¼cos w cos h � sin w cos /þ cos w sin h sin / sin w sin /þ cos w sin h cos /

sin w cos h cos w cos /þ sin w sin h sin / � cos w sin /þ sin w sin h cos /

� sin h cos h sin / cos h cos /

0B@

1CA; ð5Þ

according to the xyz-convention in the rotation order [10].Then, a position vector in the inertial reference frame, x, can be transformed into a position vector in the non-inertial ref-

erence frame, xr, through a linear transformation given by the linear translation d and the rotation matrix R:

xr ¼ RTðx� dÞ: ð6Þ

The translational velocity of the body with respect to the inertial reference frame, _d, can be related to the linear velocity inthe non-inertial reference frame using the following equation:

_d ¼ Rv: ð7Þ

The angular velocity X in the non-inertial reference frame can be related to the time rate of change of the Euler angles asfollows

_H ¼ QX; ð8Þ

where Q is the transformation matrix given by

Q ¼1 sin / sin h= cos h cos / sin h= cos h

0 cos / � sin /

0 sin /= cos h cos /= cos h

0B@

1CA: ð9Þ

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5033

It is evident Eq. (9) fails for h = ±90�. This usually is not an issue for many systems as they are not designed to have suchoperational conditions. Nonetheless, two different systems of Euler angles or quaternions can be adopted to avoid thissingularity.

2.2. Incompressible fluid flow

In the inertial reference frame defined above, the Navier–Stokes equations governing the incompressible viscous flowscan be written as:

@u@tþr � uu ¼ �rpþ mr2uþ f; ð10Þ

r � u ¼ 0; ð11Þ

where t is the time, f represents the momentum forcing term due to the immersed body, and m is the kinematic viscosity ofthe fluid. The above equations can be non-dimensionalized using a reference length L and a reference velocity U. The non-dimensionalization results in one single dimensionless parameter, the Reynolds number, defined as Re = UL/m.

2.3. Rigid body dynamics

The dynamics of a rigid body in a fluid–structure coupled system is described by the Newton–Euler equations in the body-fixed non-inertial reference frame:

msdvdtþX� v

� �¼ FNIRF þ ENIRF; ð12Þ

IsdX

dtþX� IsX ¼ TNIRF þ SNIRF; ð13Þ

where ms is the mass of the solid body, Is is the moment of inertia of the body, FNIRF and TNIRF are the fluid force and momentin the non-inertial reference frame, respectively. The force and moment FNIRF and TNIRF can be transformed to those in theinertial reference frame as

FIRF ¼ RFNIRF; ð14ÞTIRF ¼ RTNIRF; ð15Þ

respectively. Here subscripts NIRF and IRF denote quantities in the non-inertial and inertial reference frames, respectively.Likewise, all other sources of force and moment applied on the body in the non-inertial reference frame are represented byENIRF and SNIRF, respectively. The external force ENIRF can be obtained from that in the inertial reference frame using the fol-lowing transformation:

ENIRF ¼ RTEIRF; ð16Þ

noting RT ¼ R�1 as R is an orthogonal matrix. And the external force in inertial reference frame is given by

EIRF ¼ �C _d�Kdþ ðms �mfÞg; ð17Þ

with C the damping matrix, K the stiffness matrix, g the gravitational acceleration, and mf = qfVs the mass of the fluid of den-sity qf displaced by a solid body of volume Vs. Note that the external moment on the body, SNIRF, can be defined in a similarmanner except that the anchor points of the springs may be fixed in the inertial reference frame or following the rigid bodymotion, and the angular velocity X can be related to the time rate of change of the Euler angles _H directly only in one DOFrotational cases.

2.4. Fluid–structure coupling conditions

Both the kinematic and dynamic boundary conditions have to be satisfied at the fluid–structure interface. In the presentwork, non-slip velocity boundary condition is used, and the kinematic condition can be expressed as

uC ¼ Rðv þX� xr;CÞ; ð18Þ

for any point at the interface C with xr,C the position vector of such a point in the non-inertial reference frame attached tothe body. Also, a rigid body motion is imposed on the fluid enclosed by the fluid–structure interface. Hence Eq. (18) is to besatisfied everywhere occupied by the structure.

The dynamic boundary condition is given as

rf ;C � n ¼ rs;C � n; ð19Þ

where r is the stress tensor and n is the interface normal. This condition is automatically satisfied in the present direct forc-ing immersed boundary method as the Navier–Stokes equations are discretized in the whole domain including the solidphase.

5034 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

3. Computational method

The details of many aspects of the computational method in this work have been presented in several previous studies[31,33,34]. Here only a summary of each module is given and the differences from the previous work, especially, the newfluid–structure coupling algorithm, are focused on.

3.1. Basic Navier–Stokes solver

A four-step fractional-step method [6] is employed for velocity–pressure coupling, in which a pressure Poisson equationis solved to enforce the continuity equation. A semi-implicit scheme is used for time advancement with the second-orderAdams–Bashforth scheme for the convective terms and the Crank–Nicolson scheme for the viscous terms. The algorithmcan be written as follows:

u� un�1

Dt¼ �1

2ð3Cn�1 � Cn�2Þ þ 1

2bD þ Dn�1� �

�rpn�1 þ fn; ð20Þ

u� ¼ uþ Dtrpn�1; ð21Þ

r2pn ¼ 1Dtr � u�; ð22Þ

un ¼ u� � Dtrpn; ð23Þ

where superscript n denotes time step, u and u⁄ are the intermediate velocity vectors, C and D are spatial operators contain-ing the convective and viscous terms,r � uu and 1

Rer2u, respectively. The evaluation of the momentum forcing term f will be

discussed next. Note that an explicit time advancement scheme can be applied in a straightforward manner as in [31].The spatial derivatives are discretized using a second-order central difference scheme on a staggered Cartesian grid. The

approximate factorization method [3] is used for the discrete momentum equations and a parallel tridiagonal solver [17] isthen applied to invert the linear systems. The Poisson equation is solved using either a direct solver [24] or a semi-coarseningmultigrid solver [9].

3.2. Momentum forcing

To solve a FSI problem on a Cartesian grid, a geometric description of the immersed body has to be given first for definingthe body position on the grid. In this study, we use a water-tight triangulation of the body surface and track its position in aLagrangian manner. At an instance of time, given an arbitrary position vector on the body surface, i.e., the fluid–structureinterface C, in the body-fixed non-inertial reference frame, the corresponding position vector in the inertial reference framecan be defined by

xC ¼ Rxr;C þ d; ð24Þ

providing the linear translation d and the orientation defined by the rotation matrix R for this body are known.The definition of a moving immersed body with complex geometry on the Cartesian grid is a major task. In this work, a

grid point inside the solid body is identified as a solid point; if a grid point is in the fluid phase and one or more grid linesegments connecting its immediate neighboring grid points are intersected by the fluid–structure interface, then this gridpoint is defined as an interface point; all the rest grid points are defined as fluid points.

In our direct forcing immersed boundary method, a forcing term is imposed on a grid point in the solid point and interfacepoint categories to represent the effect of an immersed rigid body on the fluid flow. Therefore, the grids from both categoriesare collectively defined as forcing points for convenience. Depending on the specific approach used, the evaluation of the forc-ing term can be quite different. Here, a simple and straightforward approach is adopted by directly reconstructing the solu-tion near the immersed boundary such that the boundary conditions on the immersed boundary are satisfied. In the currentsemi-implicit time advancement scheme, following the approach in [15], a preliminary explicit step is performed first:

~u� un�1

Dt¼ � 3

2Cn�1 � 1

2Cn�2

� �þ Dn�1 �rpn�1; ð25Þ

with the viscous terms evaluated using an explicit Euler scheme. The momentum forcing term is constructed on top of theauxiliary velocity field ~u. Essentially, to satisfy the velocity distribution inside and near the immersed body, a correction hasto be applied to ~u,

unf ¼ ~uþ Dt~fn; ð26Þ

such that the effect of the immersed body is included in unf . Therefore, substituting ~u in Eq. (26) with the definition from Eq.

(25), we can obtain the momentum forcing term based on ~u:

~fn ¼ unf � ~uDt

¼ unf � un�1

Dtþ 3

2Cn�1 � 1

2Cn�2

� �� Dn�1 þrpn�1: ð27Þ

Fig. 2. Local reconstruction at interface points (M) using fluid points (h) and body forcing at body points (�) for immersed boundary treatment.

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5035

And as discussed in Section 2.4, for a rigid body motion, the velocity at an arbitrary point on the body surface or inside thesolid is given by

uns ¼ Rðvn þXn � xrÞ: ð28Þ

For a solid point, it is evident that unf ¼ un

s . For an interface point, in general it will not locate exactly on the body surface anda local reconstruction is necessary to obtain un

f . We use a linear interpolation stencil, which include one point on the im-mersed boundary ðun

s Þ and three (two for 2D cases) fluid points near the interface point, as shown in Fig. 2, to reconstructthe expected un

f using the auxiliary velocity field ~u. For details, the reader is referred to [31,33]. Note that the body positiond, orientation H, and velocity un

f , are unknown in a FSI problem and have to be obtained by solving the dynamic equationsfor rigid body motion.

3.3. Field extension reformulation

In [31], a field extension strategy was proposed to recover the physical values of pressure/velocity derivatives at the fluidpoints which were interface points at the previous time step. This field extension can be implemented by extrapolating thephysical information from the flow field into the regions containing non-physical information before the immersed bound-ary is relocated at each time step. For velocity components, the non-physical region can be easily defined as the grid pointscovered by the immersed body (for under-resolved grids or geometries with lower dimensions than the grids, the definitionwill involve those neighboring points of an interface point that are in the fluid too). However, the extension for pressure fieldis performed on points that are selected based on their association with non-physical velocity points due to the staggering ofthe grid. The reader is referred to [31] for more details. The overall field extension strategy is quite successful in eliminatingthe spurious forces otherwise generated near the immersed boundary, as shown in the results given in [31] for prescribedmotions and in [33] for predicted motions. Note that FSI applications are very sensitive to spurious force oscillation, espe-cially for cases with light structures of density close to or less than the fluid density. The robustness of the strong coupling FSIalgorithm in [33] also shows the effectiveness of the field extension strategy.

On the other hand, the straightforward implementation of field extension by modifying the solution field itself makes thedefinition of the momentum forcing term complicated. In [31,33], explicit third-order Runge–Kutta and second-orderAdams–Bashforth schemes were used for time advancement, respectively. Therefore, the direct forcing immersed boundarytreatment was fulfilled by directly modifying the predicted velocity values at forcing points and a distinct forcing term wasnot readily available for each point. In general, an explicitly defined forcing term is not required for performing the immersedboundary treatment. Actually, in [8], a semi-implicit time advancement scheme was used and the coefficients of the tridi-agonal systems resulting from the approximate factorization of the momentum equation were directly modified to respectthe boundary conditions at the immersed interfaces. In any case, the force/moment can be evaluated by a convenient surfaceintegration utilizing the interpolation stencils established for the local reconstruction at the interface points as in [31,33].However, this type of force/moment calculation can be only performed when the complete solution of a time step is avail-able, i.e., the velocity prediction, the pressure Poisson equation, and the following velocity correction steps have to be fin-

5036 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

ished. By contrast, in most direct forcing immersed boundary methods, the momentum forcing term is available before solv-ing the Poisson equation. As given in Eq. (36), the integration of the forcing function at all forcing points contains the force tothe fluid flow from the immersed boundary, since the forcing term represents the effect of the immersed body on the fluidflow [16]. Here, we shall also explore this fact to substantially accelerate the strong coupling FSI scheme proposed in [33]. Tothis end, an explicitly defined forcing term has to be given.

Interestingly, our approach to solve the above issues is quite simple and effective. In [31], the boundary condition forpressure,

@p@n

� ��C

¼ � DuDt� n

� ��C

; ð29Þ

was used to extend physical pressure information into the regions with non-physical pressure values. Notice that the pur-pose of this pressure extension is to fix the non-physical pressure gradient components at those velocity forcing points, wecan actually fix components of the pressure gradient at the velocity interface points directly. At the immersed boundary, theboundary condition for components of the pressure gradient can be written as

@p@xi

� ��C

¼ � Dui

� ��C

; ð30Þ

where xi (i = 1, 2 for 2D or i = 1, 2, 3 for 3D problems) is the ith Cartesian coordinate and ui is the corresponding Cartesianvelocity component. This equation gives the Dirichlet boundary condition for the ith component of the pressure gradient,which is also collocated with the ith velocity component. And the acceleration vector for a point at the fluid–structure inter-face is

DuC

Dt¼ d

dtðRv þ RX� xr;CÞ: ð31Þ

Therefore, the reconstruction stencil for the velocity component can be directly used to interpolate the corresponding com-ponent of the pressure gradient at the same location. Compared to the extrapolation of pressure itself in [31], this step rep-resents a great simplification of the whole procedure for the immersed boundary treatment. The pressure points requiringfield extension operation are no longer necessary to be defined and all associated data structures and computational over-head are eliminated. The only extra expense is some temporary memory space to store the pressure gradient, which has to beevaluated in the original approach anyway.

On the other hand, the pressure gradient extension is only meaningful for those fluid points at the current time step thatwere interface points at the previous time step; also, it relies on velocity reconstruction stencils for those interface points.Therefore, in order to reuse the data structures for these stencils from the previous time step, in the actual implementation,the extension was applied to all velocity interface points defined in the previous time step before we start solving the fluidflow for the current time step, as it involves repositioning the immersed body and redefining all the forcing points. Appar-ently, for the former interface points that remain their status or change to solid points at the current time step, an unclear adhoc impact will be brought into the system by the pressure gradient extension on these points. More specifically, step two(Eq. (21) in the four-step fractional-step method) removes the pressure gradient from the predicted velocity field; hence theoperation to the pressure gradient will affect the velocity field in an unpredictable manner. In fact, similar to imposing rigidbody motion on all solid points, we specify

@p@xi

� �f¼ 0; ð32Þ

for all current forcing points. It should be noted that it is applied to both the solid and interface points, even the latter are inthe fluid phase. Ideally, if the true velocity field is used together with the velocity boundary conditions at the immersedboundary to reconstruct the velocity solution at the interface points, the velocity values at these points (and the solid points,on which the solid body velocity is imposed) should not be changed in the final velocity corrector step. Therefore, all theforcing points would be functioned as velocity boundary conditions for the immersed boundary, and the treatment aboveis a practical implementation for applying the hom*ogeneous Neumann condition for pressure at these points. Of course,within the present direct forcing framework, a predicted velocity field is used in the local velocity solution reconstruction;also the hom*ogeneous Neumann condition for the pressure Poisson equation is not explicitly implemented. Nonetheless, as apart of the field extension, the step imposing Eq. (32) is essential for a much improved accuracy and the avoidance of any adhoc treatments.

Central to our reformulated field extension for the pressure gradient is the adoption of the four-step fractional-step meth-od. In the pressure increment fractional-step method used in the original approach [31,33], both the pressure and its gradi-ent, as well as the pressure increment and its gradient, appear in the whole solution process. The operation on the pressuregradient (i.e., Eq. (30)) alone will result in decoupling of the pressure field and its gradient in the system. On the other hand,it should be noted that, unlike the present method that requires the step discussed above for imposing Eq. (32), the pressurefield extension in the previous approach [31] would not bring any ad hoc impact into the solution. Therein, the modifiedpressure field is only used for calculating the pressure gradient that enters the system through the predicted velocity field.

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5037

For the former interface points that remain their status or change into solid points, their velocity values are locally recon-structed (interface points) or directly imposed (solid points) anyway, and the values of the pressure gradient at these pointsbecome irrelevant to this process. Then the gradient of the pressure increment (i.e., solution of the Poisson equation) is usedto correct the predicted velocity field, again without involving pressure itself.

In the previous field extension formulation [31], the velocity field is extended into the solid in order to obtain physicalvelocity derivatives at the interface points for the use at the next time step. The grid points requiring extension operation,as defined in [31], are relatively easier to identify than those for pressure field extension. But similar to the pressure points,extra steps are necessary for the fulfillment of field extension at grid points near sharp corners or under-resolved thin struc-tures. It is instinctive to ask if the velocity field extension could be simplified too. Actually, in our numerical experiments, wefound that this procedure could be removed without much negative effect on the results when the grid resolution is reason-able, although a slight improvement was observed when it was used on a very coarse grid. This is quite different from thecase of pressure field extension, which is essential to moving boundary problems in our direct forcing approach. There aretwo major reasons for the differences. First, for an incompressible viscous flow with immersed solid bodies, the velocity fieldis continuous across the fluid–solid interface when a no-slip boundary condition at the interface and a rigid body motioninside the solid are imposed; but the pressure field has a jump condition at the interface [30]. Second, in our direct forcingapproach, the momentum forcing function is not regularized to have a smooth transition near the interface, thus the pres-sure field retains its rapid change across the interface; on the other hand, we evaluate the momentum forcing function uti-lizing the predicted velocity field and do not impose the pressure jump conditions in the Poisson equation, hence as a resultof the current fractional-step type of velocity–pressure coupling, the sharp variation of the velocity field across the interfaceis mitigated to some extent. Nevertheless, we shall demonstrate in the results section that, although the whole field exten-sion process is greatly simplified, the reformulation of the pressure field extension and the elimination of the velocity fieldextension discussed here do not alter the overall accuracy or the numerical stability property of the original approach. Actu-ally, larger time steps can be used with the present method as a semi-implicit time-advancement scheme is adopted.

3.4. Force and moment evaluation

For rigid body dynamics considered here, essentially, the integration of the forcing term f in the momentum equation foran immersed body,

Rqf f dV , represents the total effect of this body on the fluid flow, including the force, i.e., mf

€d, required toimpose a rigid body motion of the portion of fluid enclosed by the immersed boundary. For a forcing point, either an interfacepoint or a solid point, the forcing term is given by Eq. (27). As discussed in the previous part, Eq. (32) is applied to all forcingpoints. Therefore, the actual forcing term is

fn ¼ unf � un�1

Dtþ 3

2Cn�1 � 1

2Cn�2

� �� Dn�1; ð33Þ

or,

fn ¼ unf � ~uDt

�rpn�1; ð34Þ

with the auxiliary velocity field ~u still evaluated using Eq. (25) and the pressure field extension fulfilled by interpolating thepressure gradient for all interface points at the previous time step using Eq. (30). In a strong coupling FSI scheme, the posi-tions of the immersed bodies have to be adjusted multiple times in one time step, hence the forcing points for a body have tobe re-defined each time when the position of a body is changed. It is particularly convenient to adopt Eq. (34) for evaluatingthe momentum forcing term, as the pressure field extension procedure is significantly simplified by applying Eq. (32) onlyonce before Eq. (21).

In addition, a comparison of Eqs. (20) and (25) shows that the pressure gradient term offsets for all current forcing points.Also, the desired velocity un

f , either reconstructed using the predicted velocity field at the surrounding fluid points for theinterface points, or directly given by the predicted rigid body motion for the solid points, is not affected by the pressure gra-dient term at these points. Therefore, the value of the pressure gradient term at a forcing point is not relevant to the actualforcing term added to the momentum equation, which means it is totally acceptable to apply the first step of the pressurefield extension given in Eq. (30) to all previous interface points, regardless of their status changes later. Together with thesecond step in Eq. (32), it is apparent that the operation on the pressure gradient does not affect the forcing term at all,as the changes added to the system through the pressure gradient in step one (i.e., Eq. (20)) will be cancelled out in steptwo (i.e., Eq. (21)) entirely!

Therefore, the force exerted on the body by the fluid can be given as

FIRF ¼ �Z

qf f dV �mf€d

� �; ð35Þ

in the inertial reference frame. For solving the Newton–Euler equations, we need the fluid force and moment on the body inthe non-inertial reference frame. They can be written as

5038 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

FNIRF ¼ �RTZ

qf f dV �mf€d

� �¼ � RT

Zqf f dV �mf

dvdtþX� v

� �� ; ð36Þ

TNIRF ¼ � RTZðx� dÞ � qf f dV � If

dtþX� IfX

� �� ; ð37Þ

respectively. If is the moment of inertia of the portion of fluid enclosed in the volume Vs defined by the immersed boundary.

3.5. Predictor–corrector scheme for rigid body dynamics equations

As in [33], Hamming’s 4th-order predictor–corrector scheme is used to solve the governing equations for rigid bodydynamics. First, the equations for rigid body dynamics are rewritten as follows

_vn ¼ m�1s Fn

NIRF þ EnNIRF �msX

n � vn

; ð38Þ_Xn ¼ I�1

s TnNIRF þ Sn

NIRF �Xn � IsXn : ð39Þ

We also rewrite Eqs. (7) and (8) as

_dn ¼ Rnvn; ð40Þ_Hn ¼ Q

nXn; ð41Þ

respectively. Then, we can solve a system of first-order ordinary differential equations (ODE) with a numerical integrationscheme. In general, one equation from this system can be given as

dqdt¼ _q; ð42Þ

where q represents the displacement or velocity in a certain direction. With Hamming’s scheme [13], we have three steps:

(i) predictor step: the solution is explicitly predicted and then modified using the error estimated from the previous timestep

qnp ¼ qn�4 þ 4

3Dtð2 _qn�1 � _qn�2 þ 2 _qn�3Þ; ð43Þ

qk¼0c ¼ qn

p �112121ðqn�1

p � qn�1c Þ; ð44Þ

where subscripts p and c represent the predicted and corrected values, respectively, superscript k is the iteration index forthe corrector step, and the last term on the right hand side of Eq. (44) is the truncation error estimated from the previoustime step.

(ii) corrector step: the solution is corrected iteratively

for k ¼ 1; . . . ;m; . . .

qmc ¼

18ð9qn�1 � qn�3Þ þ 3

8Dtð2 _qm�1

c þ 2 _qn�1 � _qn�2Þ

if jqmc � qm�1

c j < �; qnc ¼ qm

c ; exit: ð45Þ

where the tolerance � = 10�8 is used in the present study. And this convergence criterion has to be met for all equations inthe system before the iteration can be terminated.

(iii) finalizer step: the solution is updated using the estimated truncation error

qn ¼ qnc þ

9121ðqn

p � qnc Þ: ð46Þ

This scheme is very robust and efficient as shown in [33].

3.6. Fluid–structure coupling procedure

The fluid–structure coupling procedure in one time step can be given as follows:

(i) Perform the reformulated field extension for pressure by reconstructing components of the pressure gradient at theforcing points from the previous time step using Eq. (30).

(ii) Perform the preliminary velocity predictor step using Eq. (25).(iii) Perform the actual velocity predictor step without the forcing term.(iv) This step has three branches:

� If it is the first time reach this step, perform the structural predictor step.

Fig. 3.dotted

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5039

� If it is not the first time and the ODEs are not converged, perform the structural corrector step and checkconvergence.

� If the ODEs are converged or the maximum number of iterations is reached, perform the structural finalizer step.

(v) With the structural displacement and velocity from the above step, redefine the structural position and the fluid/solid/

interface status of all grid points, and set up the local reconstruction and field extension stencils.(vi) Define and integrate the momentum forcing term to obtain the fluid force and moment. Invert the governing equa-

tions for rigid body motion to acquire the acceleration of the structures.(vii) If this step follows the structural finalizer in step (iv), then go to the next step; otherwise go back to step (iv).

(viii) Solve the tridiagonal systems to obtain the predicted velocity field with the forcing term added (the effect of immersedboundaries is included here).

(ix) Solve the Poisson equation to obtain the pressure field for the current time step.(x) Perform the velocity corrector to obtain the velocity field for the current time step.

In the above procedure, the fluid solver part is essentially excluded from the iterative loop, which starts at Step (iv) and exitat Step (vii). It only includes the solution of the ODE system, the setup of immersed boundaries, and the definition of themomentum forcing term. Compared to the original approach in [33], this represents a substantial computational accelerationof several times or one order of magnitude, depending on the number of iterations it may take. This is because the fluid sol-ver with the expensive Poisson equation inversion is only performed once at each time step.

Furthermore, the approach proposed here can be applied to all other immersed boundary and immersed interface meth-ods using a momentum forcing function, either continuous or discrete. In addition, multiple structures or flexible structurescan be treated in a straightforward manner.

4. Results

4.1. Prescribed motions

4.1.1. In-line oscillation of a circular cylinderIn this case, a circular cylinder is given a prescribed translational motion in a fluid at rest. The equation of motion is given

by x(t) = �Asin(2pft), where A is the amplitude of the oscillation. Two key parameters are the Reynolds number, Re = UmaxD/m, and the Keulegan–Carpenter number, KC = Umax/fD, where Umax is the maximum velocity of the cylinder, D is the cylinderdiameter, and f is the characteristic frequency of the oscillation. The parameters from the experiments and numerical sim-ulations in [7] are used here, i.e., Re = 100 and KC = 5. The size of computational domain is 40D � 24D in the x and y direc-tions, respectively, with the cylinder located at the center initially. Three different grids of 240 � 120, 480 � 240, and960 � 480 are used for this case and around the cylinder the grids are approximately uniform in both directions with spacingof 0.04D, 0.02D, and 0.01D, respectively. A constant time step of 0.005D/Umax is used on all three grids. Fig. 3 shows the timehistory of the in-line force, Fx(t), acting on the cylinder. Results from all three grids are in very good agreement with the ref-erence data from [7]. This validates the accuracy of our simplified approach for applying the field extension.

Fig. 4 shows the pressure and vorticity contours at four different phase-angles from the fine grid simulation. The resultsare in very good agreement with the corresponding results reported in [7,31].

4.1.2. Transversely oscillating cylinder in a free-streamTo further demonstrate that the present simplified field extension strategy does not alter the overall accuracy of the pre-

vious approach [31], the case of a transversely oscillating circular cylinder in a free-stream, which was considered in [31], is

/(ρD

U 2 )

Time history of the in-line force acting on a circular cylinder oscillating in a fluid at rest at Re = 100 and KC = 5. Symbols: experimental data [7];line: grid 240 � 120; dashed line: grid 480 � 240; solid line: grid 960 � 480.

Fig. 4. Pressure (left) and vorticity (right) contours at four different phase angles for an in-line oscillating circular cylinder in a fluid at rest at Re = 100 andKC = 5: (a) / = 0�; (b) / = 96�; (c) / = 192�; and (d) / = 288�.

5040 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

also examined here. The Reynolds number based on the free-stream velocity U and the cylinder diameter D is Re = UD/m. Thetransverse oscillation of the cylinder is given by y(t) = Acos(2pfet), with y(t) the transverse position of the cylinder center, Athe amplitude of the oscillation, and fe the excitation frequency. The parameters are chosen to match those used in [12,31],i.e., Re = 185, A = 0.2D, and 0.8 6 fe/f0 6 1.2, with f0 the natural vortex shedding frequency for a stationary circular cylinder.

The computational domain is [�10D,40D] � [�15D,15D] and the grid is 400 � 300 points in the streamwise and trans-verse directions, respectively. The grid is the same as the medium grid used in [31] and the local grid spacing around thecylinder is 0.01D. At the inlet Dirichlet velocity boundary condition is applied. A convective boundary condition is used atthe outlet. Slip-wall condition is used for far-field boundaries. Four cases, i.e., fe/f0 = 0.8, 1.0, 1.1, and 1.2, are considered.For each case, a constant time step between 0.004 � 0.005D/U (the specific value depends on the excitation frequency fe)is used.

Fig. 5 shows the instantaneous vorticity fields when the cylinder is at the positive peak of the oscillation. Compared tocases fe/f0 = 0.8 and 1.0, a topological change of the vorticity distribution in the near wake for case fe/f0 = 1.1 is evident.The overall patterns for these cases are very similar to those reported in [31].

Fig. 6 shows the drag and lift coefficients as functions of time for all four cases. The drag and lift coefficients are defined asCD ¼ fx

12 qf DU2� �.

and CL ¼ fy12 qf DU2� �.

with fx and fy the instantaneous drag and lift forces on the cylinder, respectively. In

Fig. 5. Snapshots of the vorticity field for the transversely oscillating circular cylinder in a free-stream at Re = 185: (a) fe/f0 = 0.8; (b) fe/f0 = 1.0; (c) fe/f0 = 1.1;and (d) fe/f0 = 1.2.

Fig. 6. Drag and lift coefficients (upper and lower curves in each sub-figure, respectively) as functions of time for the transversely oscillating circularcylinder in a free-stream at Re = 185: (a) fe/f0 = 0.8; (b) fe/f0 = 1.0; (c) fe/f0 = 1.1; and (d) fe/f0 = 1.2. Solid line: previous immersed boundary method [31];dashed line: present immersed boundary method.

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5041

order to have a direct comparison, the present results are integrated from the local surface force distribution, same as theapproach used to obtain the reference results on the same grid in [31]. The curves of time evolution from the present methodand the method in [31] are almost on top of each other, except for some deviation due to different transient phases for casefe/f0 = 1.2. The present mean drag coefficients are 3% higher than their corresponding values in [31], which might be an effectof the different treatments of the pressure field extension.

Fig. 7 shows the distribution of the pressure and skin friction coefficients on the cylinder surface at the instant of the po-sitive peak of the cylinder oscillation. Both the results from a body-fitted method in [12] and the previous immersed bound-ary method in [31] are given for reference. The results from [31] were obtained on the same 400 � 320 grid used here. The

-32

-3

Fig. 7. Distribution of the pressure and skin friction coefficients on the cylinder surface at the instant of positive oscillation peak: (a) fe/f0 = 0.8; (b) fe/f0 = 1.0; (c) fe/f0 = 1.1; and (d) fe/f0 = 1.2. Symbols: body-fitted method in a non-inertial reference frame [12]; solid line: previous immersed boundarymethod on the 400 � 320 grid in [31]; dashed line: present immersed boundary method on the 400 � 320 grid.

5042 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

pressure coefficient is not very sensitive to grid resolution and the curves from the Cartesian grid methods match those froma boundary-conforming method very well. The skin friction coefficient, proportional to the surface vorticity distribution inthe present 2D case, is much sensitive to grid resolution; and the local grid spacing of 0.01D in the Cartesian grid is not fineenough to capture the peaks of its variation on the cylinder surface. In general, the present results can hardly be distin-guished from those in [31] in terms of accuracy when compared with the results in [12] using a body-fitted method. Thisshows the present simplified method does not change the overall accuracy of the method in [31].

4.1.3. A hovering wingThe two cases discussed above only contain translational motions of a circular cylinder in one direction. The fluid–struc-

ture interaction becomes more complicated when a case combines rotation and translational motions in more than onedirection. Wang [29] studied the 2D mechanism for insect hovering using a single wing with a prescribed motion. The wingis a 2D ellipse with chord length C and the aspect ratio between the major and minor axes is K. The hovering flight has aflapping period T and a flapping distance A0 for the wing center translation. The stroke plane is inclined at an angle b.The prescribed motion can be defined by giving the center coordinates (x,y) and the angle of attack for the wing as follows

xðtÞ ¼ A0

2cos

2ptT

� �þ 1

� �cos b; ð47Þ

yðtÞ ¼ A0

2cos

2ptT

� �þ 1

� �sin b; ð48Þ

aðtÞ ¼ a0 1� sin2pt

Tþ /

� �� ; ð49Þ

where / is the phase difference. The maximum speed along the flapping path, U = pA0/T, is used as the velocity scale. And theReynolds number is defined as Re = pA0C/Tm. The parameters used here are chosen to be the same as those used in [29,30]:K = 4, A0 = 2.5C, a0 = p/4, and Re = 157.

A square computational domain of [�24C,24C] � [�24C,24C] is used with the wing center initially located at (0,0). A non-uniform grid of 576 � 576 covers the whole domain and the grid spacing around the wing is 0.03C. This is very close to thegrid resolution used by Xu and Wang [30] in their simulation with an immersed interface method, in which a uniform grid of512 � 512 was used for a much smaller domain of [�10C,6C] � [�10C,6C]. A constant time step, Dt = 0.008C/U is used here.

Fig. 8 shows four snapshots of the vorticity field around the hovering wing during one flapping period. The vortex dipoleformed and then shed during each period is the key to the lift generation in the 2D hovering flight. These snapshots are verysimilar to those given in [29,30]. In the latter, there was a clockwise flow at the lower right part of the plots due to the small

Fig. 8. Snapshots of the vorticity field around a hovering wing of elliptical shape at Re = 157: (a) t/T = 0.25, during the downstroke, a pair of counter-rotatingvortices is generated; (b) t/T = 0.44, the pair of vortices sheds as the wing rotates; (c) t/T = 0.74, during the upstroke, the pair of vortices becomes a vortexdipole; and (d) t/T = 0.99, the vortex dipole moves downward.

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5043

size of the computational domain. The vortex dipoles shed from the wing were driven toward the left side by this flow. Thisphenomenon does not happen in the present simulation as a much larger domain is adopted to minimize the blockage effectfrom the domain boundary.

Time series of drag and lift coefficients are given in Fig. 9. Computed results from a conformal mapping finite differencemethod in a wing-fixed non-inertial reference frame [29] and a Cartesian grid immersed interface method [30] are also givenfor comparison. The results from the present direct forcing immersed boundary method on a grid with similar resolution tothat in [30] match the data from the immersed interface method very well. In addition, as mentioned above, the present sim-ulation is performed with a larger domain, thus the resulting drag and lift coefficients have better periodicity than those in[30].

4.2. Vortex-induced vibration of a circular cylinder

The FSI problem of an elastically mounted circular cylinder in a free-stream is considered to show the accuracy of our newformulation. The cylinder can be modeled as a mass-spring-damper system, and the governing equations for the cylindermotion are as follows:

€X þ 2f2pU�

� �_X þ 2p

U�

� �2

X ¼ 2pm�

CD; ð50Þ

€Y þ 2f2pU�

� �_Y þ 2p

U�

� �2

Y ¼ 2pm�

CL; ð51Þ

where the mass ratio is m⁄ = ms/mf, the reduced velocity is U⁄ = U/fND. The natural vibration frequency of the structure andthe damping ratio are given by fN ¼ 1

ffiffiffiffiffiffiffiffiffiffiffik=ms

pand f ¼ c=ð2

ffiffiffiffiffiffiffiffiffikms

pÞ, respectively, where k is the spring constant and c is the

damping coefficient. Also, X = x/D and Y = y/D are the non-dimensionalized horizontal and vertical displacements with xand y the streamwise and transverse displacements of the cylinder center, respectively. Note that the same length and veloc-ity reference scales as in the Navier–Stokes equations (the cylinder diameter, D, and the free-stream velocity, U) were usedhere.

The parameters are selected to match the spectral element simulation in a non-inertial reference frame in [5]. The Rey-nolds number is Re = 200, the reduced velocity is U⁄ = 5.0, the damping ratio is f = 0.01, and the mass ratio is m⁄ = 4/p. In oursimulation, the inflow of a uniform velocity U is located 10D upstream of the cylinder. The outflow boundary is located 30Ddownstream of the cylinder and a convective boundary condition is applied. Both top and bottom boundaries are located 10Daway from the cylinder and a free-stream condition is used. Three different grids are used to cover the whole domain and thenumber of grid points for the three grids are 160 � 120, 320 � 240, and 640 � 480 in the streamwise and transverse direc-tions, respectively. The grids are stretched in both directions. Near the cylinder, the grids are uniform with local spacing ofsize 0.04D � 0.04D, 0.02D � 0.02D, and 0.01D � 0.01D, respectively. A constant time step of 0.004D/U is used for all threegrids.

Fig. 9. Drag and lift coefficients as functions of time for the 2D hovering flight of an elliptical wing at Re = 157. Dashed line: conformal mapping finitedifference method in a non-inertial reference frame [29]; solid line: immersed interface method [30]; symbols: present direct forcing immersed boundarymethod.

Fig. 10. Centerline displacement and displacement–velocity phase plots for a circular cylinder freely vibrating in a free-stream at Re = 200 and m⁄ = 4/p. �:spectral element simulation in a non-inertial reference frame [5]; dotted line: previous immersed boundary simulation on grid 640 � 480 [33]; solid line:present method on grid 160 � 120; dashed line: present method on grid 320 � 240; dot-dashed line: present method on grid 640 � 480.

5044 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

On all three grids, the simulations are carried out first for a stationary cylinder. After the periodic vortex shedding patternis established, the cylinder is released to allow two DOF translational motions in both streamwise and transverse directions.The frequency of the cylinder oscillation is 0.187, same as the frequency of the vortex shedding, which means the presentcase is in the synchronization regime of vortex-induced vibration. Fig. 10 shows the centerline displacement and displace-

Fig. 11. Instantaneous vorticity field for a circular cylinder freely vibrating in a free-stream at Re = 200 and m⁄ = 4/p.

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5045

ment–velocity phase plots. The grid convergence is evident from the phase plots. The maximum amplitude of the cylinderoscillation on the fine grid is 0.603D. The centerline displacement data from [5] using a high order spectral element methodin a non-inertial reference frame and from a previous immersed boundary simulation [33] are also given in the figure. Thecenterline displacement data have been shifted to have the ‘‘figure of eight’’ plots from all cases centered at (0,0). Our finegrid results match the reference data very well and there is a clear trend of grid convergence of the solution, which can beseen from the displacement–velocity phase plots as well. On the fine grid, the center of cylinder oscillation is at x = 0.651D,which is slightly further downstream than the value 0.62D from the spectral element simulation in [5]. This gives a 5% dif-ference of the cylinder displacement in the streamwise direction, which is due to a difference of the same percentage in thecomputed mean drag forces between the two simulations, as evident from Eq. (51) that the spring extension is given by thecylinder displacement and the mean drag force is balance by the mean spring force in the streamwise direction. Fig. 11shows the instantaneous vorticity contours from the free oscillating cylinder, which gives a clear 2S (two single vortices shedduring each oscillation period) vortex pattern that has been observed in many experiments and simulations.

To further demonstrate the robustness of the present strongly coupled FSI algorithm, the vortex-induced vibration casefor a circular cylinder is performed for a much lower mass ratio, m⁄ = 0.5. The computational domain, initial and boundary

Fig. 12. Centerline displacement and displacement–velocity phase plots for a circular cylinder freely vibrating in a free-stream at Re = 200. Solid line:m⁄ = 4/p; dashed line: m⁄ = 0.5.

Fig. 13. Instantaneous vorticity field for a circular cylinder freely vibrating in a free-stream at Re = 200 and m⁄ = 0.5.

5046 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

conditions, and other parameters are the same as used in the m⁄ = 4/p case on the 640 � 480 grid: Re = 200, U⁄ = 5.0, f = 0.01,and Dt = 0.004D/U. As shown in Fig. 12, after the transient the cylinder retains a very stable oscillation and exhibits a ‘‘figure-of-eight’’ trajectory. Compared to the m⁄ = 4/p case, the cylinder oscillates around a center further downstream and at a high-er amplitude (0.668D). This case is also in the synchronization regime of vortex-induced vibration as both the oscillation andvortex shedding frequencies are 0.175. Fig. 13 shows an instantaneous vorticity field. Due to a lower shedding frequency,there are less vortices in the wake compared to the m⁄ = 4/p case.

4.3. Galloping rectangular bodies

Vortex-induced vibration of a structure, as one type of flow-induced vibration, is characterized by the lock-in or synchro-nization phenomenon, in which the structure oscillates at a frequency similar to that of the vortex shedding behind thestructure. Another type of flow-induced vibration happens at frequencies much lower than the wake vortex shedding fre-quencies, usually associated with a non-circular structure in a flow with a large incoming velocity. If the structure is limitedto have one DOF motions, i.e., transverse or rotational oscillations, then this type of vibrations is called galloping. Robertsonet al. [22] studied transversely and rotationally galloping rectangular bodies using a high-order spectral element method in anon-inertial reference frame fixed to the body. In this part, two cases in [22], one transversely galloping body and one rota-tionally galloping body in a free-stream, are selected to validate the present FSI algorithm. For rotational galloping, the casewithout damping is also presented. A case with both a transversely galloping body and a rotationally galloping body is stud-ied to demonstrate the applicability of our algorithm to FSI problems with multiple bodies. For a body with a rectangularcross-section of thickness D and width L, a width-to-thickness ratio can be defined as K = L/D. This ratio is one of the majorparameters for determining a rectangular body is susceptible to galloping motion or not.

4.3.1. Transverse gallopingFor the transverse galloping, the parameters are selected to match the spectral element simulation of a rigid body with a

square cross-section in [22], i.e., K = 1. The governing equation for the transversely galloping motion of a rectangular body issimilar to Eq. (51), but with a different right-hand-side, i.e., CL/(2m⁄), instead of CL=

12 pm�

, as the displaced fluid volume perunit length is qfDL for a body with a rectangular cross-section instead of qfpD2/4 for that with a circular cross-section. TheReynolds number based on the free-stream velocity U and the body thickness D is Re = 250. Parameters in the structural gov-erning equation are: the reduced velocity U�y ¼ 40; the damping ratio fy = 0.0037; and the mass ratio m⁄ = 20. The computa-tional domain is [�10D,30D] � [�11D,11D]. Boundary conditions are the same as given in the previous case. A non-uniformgrid of 480 � 480 is used with a uniform grid spacing of 0.015D around the square body. The constant time step isDt = 0.006D/U.

A simulation with the body fixed (centered at (0,0)) for the same Reynolds number is performed first. Then the FSI sim-ulation is restarted from the solution field with the fixed body by allowing the body to have a one DOF transversely heavingmotion. Fig. 14 shows the transverse response of the body in the cross-flow. As a direct comparison, the data from Robertsonet al. [22] using a body-fitted spectral element method in a non-inertial reference frame are also given. The excellent agree-ment of the results from two totally different approaches is remarkable. In our case, the maximum response of the body is1.146D, and the frequency of the body oscillations is 0.0238, which is close to the natural frequency of the structure (0.952fy,where the natural frequency fy = 0.025). To further verify the accuracy of the present algorithm, two more grids are used forthis case. The fine one has 720 � 720 points with a uniform grid spacing of 0.01D around the body. The coarse grid is ob-tained from the fine grid by taking every other point, thus it has 360 � 360 points with a uniform grid spacing of 0.02Daround the body. And the constant time steps for these two grids are Dt = 0.004D/U and 0.008D/U, respectively. All otherparameters in the simulations are the same as the medium grid (480 � 480). The maximum responses from the coarseand the fine grids are 1.073D and 1.144D, respectively; whereas the oscillation frequency remains the same on all three gridsused. This demonstrates that the accuracy of the present FSI algorithm as the results from different grids with different timesteps are very close to each other (the difference between the maximum responses from the coarse and the fine grids is about6%).

Fig. 14. Transverse response as a function of time for a square body galloping in a cross-flow at Re ¼ 250; U�y ¼ 40; fy ¼ 0:0037, and m⁄ = 20. Solid line:spectral element simulation in a non-inertial reference frame [22]; dashed line: present immersed boundary simulation.

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5047

Fig. 15(a) gives an instantaneous vorticity field for the square body galloping in the cross-flow. Due to the rectangularshape of the body, boundary layer separations at the sharp corners are evident. As a comparison, a snapshot of the vortexshedding from the fixed square body is shown in Fig. 15(b), which has a frequency of 0.144. The spacing and strength ofthe shed vortices from both cases are very similar. For the former case, due to the high amplitude, low frequency galloping

Fig. 15. Comparison of vortex shedding patterns behind a body with a square cross-section in a free-stream at Re = 250: (a) transverse galloping atU�y ¼ 40; fy ¼ 0:0037, and m⁄ = 20; (b) fixed; and (c) vortex-induced vibration at U�y ¼ 7; fy ¼ 0:0037, and m⁄ = 20.

Fig. 16. Transverse response as a function of time for the vortex-induced vibration of a square body in a cross-flow at Re ¼ 250; U�y ¼ 7; fy ¼ 0:0037, andm⁄ = 20.

5048 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

response, the vortex shedding forms an oblique wake structure (deviates from the centerline and toward the upper rear endof the domain at the given instant).

To further illustrate the difference between galloping and vortex-induced vibration, a case with a much lower reducedvelocity, U�y ¼ 7, is performed with all other parameters the same as the galloping case. The vortex shedding, featured bythe sharp corner separations as shown in Fig. 15(c), is similar to the case of fixed body. However, the spacing between vor-tices in the streamwise direction is reduced due to the higher vortex shedding frequency at 0.174, compared to the fixedbody case. The time series of transverse response is shown in Fig. 16. The vortex-induced vibration has a much lower ampli-tude than that of the galloping motion at the much higher reduced velocity. Also, the oscillation frequency is the same as thevortex frequency, which is the typical lock-in phenomenon in vortex-induced vibration.

4.3.2. Rotational gallopingThe governing equation for the one DOF torsional model of a rectangular body can be given as follows

Fig. 17I⁄ = 400

€hþ 2fh2pU�h

� �_hþ 2p

U�h

� �2

h ¼ CT

2I�; ð52Þ

where h is the rotational angle of the body about its center, the moment of inertia ratio is I⁄ = Is/(qs D4) withIs ¼ 1

12 qsDLðD2 þ L2Þ the moment of inertia of the body per unit length, the rotational moment coefficient is

CT ¼ T 12 qf D

2U2� �.

with T the rotational moment on the body per unit length, the reduced velocity is U�h ¼ U=ðfhDÞ. The nat-

ural vibration frequency of the structure and the damping ratio are given by fh ¼ 12p

ffiffiffiffiffiffiffiffiffiffikh=Is

pand fh ¼ ch=ð2

ffiffiffiffiffiffiffiffikhIs

pÞ, with kh the

torsional spring constant and ch the torsional damping coefficient, respectively.The simulation parameters are selected to match the spectral element simulation of a rigid body with a rectangular cross-

section of K = 4 in [22]. The Reynolds number based on the free-stream velocity U and the body thickness D is Re = 250. Otherparameters in the structural governing equation are: the reduced velocity U�h ¼ 40; the damping ratio fh = 0.25; and the massmoment of inertia ratio I⁄ = 400. The computational domain is [�11D,50D] � [�21D,21D]. Boundary conditions are the sameas given in the previous case. A non-uniform grid of 320 � 400 is used with a uniform grid spacing 0.03D around the rect-angular body. The constant time step is Dt = 0.006D/U.

. Rotational galloping response as a function of time for a rectangular body freely rotating in a cross-flow at Re ¼ 250; U�h ¼ 40; fh ¼ 0:25, and. Solid line: spectral element simulation in a non-inertial reference frame [22]; dashed line: present immersed boundary simulation.

Fig. 18. Snapshots of the vorticity field around a rectangular body freely rotating in a free-stream at Re ¼ 250; U�h ¼ 40; fh ¼ 0:25, and I⁄ = 400: (a)h = 15.6�; (b) h = 1.7�; (c) h = �14.8�.

Fig. 19. Rotational galloping response as a function of time for a rectangular body freely rotating in a cross-flow at Re ¼ 250; U�h ¼ 40; fh ¼ 0, and I⁄ = 400.

Fig. 20. Snapshots of the vorticity field around a rectangular body freely rotating in a free-stream at Re ¼ 250; U�h ¼ 40; fh ¼ 0, and I⁄ = 400: (a)h = �121.1�; (b) h = 16.1�; (c) h = 122.7�.

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5049

Similar to the case of transverse galloping, the simulation is started with the body fixed. After a stable vortex shedding isestablished the body is allowed to have a one DOF rotational motion. Fig. 17 shows the time series of the rotational gallopingresponse of the body in the free-stream. As a direct comparison, the data from Robertson et al. [22] using a body-fitted spec-tral element method in a non-inertial reference frame are also given. The maximum rotational responses from the presentsimulation and [22] are 15.7�and about 15� (estimated from the solid line in Fig. 17), respectively. The difference is within5%. Also the present simulation gives a higher frequency at 0.0198, comparing to that in [22] (0.0191). The overall agree-ment between these two totally different approaches is very good.

Fig. 18 gives three snapshots of the instantaneous vorticity field around the body during a clockwise rotation from a po-sitive peak to the following negative peak. For a large rotational angle, the boundary layer separates at the leading sharpcorner of a leeward surface; the flow reattaches the surface near the rear end and then separates at the trailing edge again.But for a very small rotational angle, the flow does not reattach and features a large angle separation due to the clockwiserotation of the body.

Robertson et al. [22] also demonstrated a simulation of the rectangular body but with a zero damping ratio (fh = 0). How-ever, due to the restrictions of their algorithm from using a non-inertial reference frame, their simulation was limited to amaximum rotational response much less than 90� and had to be terminated before a stable galloping pattern could be estab-lished. On the other hand, the present algorithm does not have such type of restrictions and the mentioned case can be

5050 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

performed in a straightforward manner, just like other cases. Here the computational setup and all parameters are the sameas listed above except that fh = 0 and Dt = 0.004D/U due to the violent rotational galloping response as shown in Fig. 19. Themaximum galloping amplitude is 123�, about 8 times of that for the fh = 0.25 case. The galloping frequency is 0.0244, whichis also higher than that of the fh = 0.25 case.

Fig. 20 shows three snapshots of the vorticity field around the body during an anti-clockwise rotation from the negativeresponse peak to the positive peak position. Near a maximum response position, the body faces the free-stream with a largeamplitude angle of attack; the flow features rapid vortex shedding from the trailing edge of the windward surface and largerecirculation zone with several major inter-interacting vortices behind the body. Near a zero response position, due to themuch reduced blockage of the body, these vortices are washed down in the wake.

4.3.3. Two galloping rectangular bodiesFor a FSI algorithm relying on a non-inertial reference frame, besides the large angle rotation limitation shown in the zero

damping ratio case in the previous section, another major disadvantage is that it is limited to cases with a single body in anunbounded domain. Here a case with multiple bodies is studied to demonstrate the applicability of the present algorithm.

The transversely galloping body and the rotationally galloping body studied in the previous sections are placed in a free-stream in tandem arrangement. The same Reynolds number is used, Re = UD/m = 250. For the transversely galloping body of asquare cross-section (K = 1), the structural parameters are: U�y ¼ 40; fy ¼ 0:0037, and m⁄ = 20. The rotationally gallopingbody of a rectangular cross-section (K = 4) is placed behind the square body with an initial distance 4.5D between thetwo centroids and has the following structural parameters: U�h ¼ 40; fh ¼ 0:25, and I⁄ = 400. The computational domain is[�11D,52D] � [�21D,21D]. Boundary conditions are the same as given in the previous cases. A non-uniform grid of480 � 320 is used with a uniform grid spacing 0.03D around the two bodies. The constant time step is Dt = 0.006D/U.

This case is a nice illustrative example of body–body interactions through vortex shedding and wake flows. Fig. 21 showsthe comparison of the galloping responses of the body of a square cross-section in the single-body and the two-body con-figurations. Overall the upstream body is barely affected by the downstream body. The galloping amplitude (about 1.27D) isa little higher than that in the single-body case; whereas the response frequency (0.0235) is slightly decreased. In general,the shed vortices from the upstream body are forced to circumvent the downstream body due to the blockage of the latter,

tU/D

y/D

0 200 400 600 800 1000-2

-1

Fig. 21. Responses as functions of time for the transversely galloping body of a square cross-section at Re ¼ 250; U�y ¼ 40; fy ¼ 0:0037, and m⁄ = 20. Dashedline: a single body in a free-stream; solid line: two bodies in tandem arrangement in a free-stream.

tU/D

θ (d

egre

0 200 400 600 800 1000-20

-10

Fig. 22. Responses as functions of time for the rotationally galloping body of a rectangular cross-section at Re ¼ 250; U�h ¼ 40; fh ¼ 0:25, and I⁄ = 400.Dashed line: a single body in a free-stream; solid line: two bodies in tandem arrangement in a free-stream.

Fig. 23. Snapshots of the vorticity field around two rectangular bodies in tandem arrangement freely galloping in a free-stream at Re = 250: (a) y/D = 1.27;h = �6.1�; (b) y/D = 0.75; h = �9.4�; (c) y/D = �0.06; h = �8.2�; (d) y/D = �1.11; h = 0.3�; (e) y/D = �1.30; h = 6.4�; (f) y/D = �0.83; h = 8.7�. For the upstreambody: U�y ¼ 40; fy ¼ 0:0037, and m⁄ = 20. For the downstream body: U�h ¼ 40; fh ¼ 0:25, and I⁄ = 400.

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5051

which provides extra lift force to the upstream body and thus increases its transverse galloping response. And the period fora complete galloping cycle becomes slightly longer too.

Fig. 22 shows the comparison of the galloping responses of the body of K = 4 in the single-body and the two-body con-figurations. It is evident that the galloping response of the downstream body is greatly affected by the upstream body. Itexhibits a much reduced rotational amplitude (about 10�) and the response frequency turns out to be 0.0235, exactly thesame as that of the upstream body and much higher than that in the single body configuration. The phase angle betweenthe responses of these two body is 121�, which apparently is affected by the gap width between the bodies. With the currentnarrow gap (2D), the rotational response of the downstream body is significantly delimitated by the bypassing vorticesformed behind the upstream body and the galloping frequency is modified to follow that of the upstream body.

Fig. 23 shows several snapshots of the vorticity field around the bodies in the course of the upstream body heaving fromthe positive translational peak (a) to the negative peak position (e) and the downstream body pitching from the negativerotational peak (b) to the positive rotational peak position (f). The deflection of the latter to the vortices shed from the formeris evident. These deflected and elongated vortices interact with the vortices shed from the leading and trailing edges of thedownstream body as they bypass it. Also their alternating circumvention above or below the downstream body restricts itsgalloping amplitude, and modulates its response frequency to be the same as that of the upstream body. These vortices thenjoin those shed from the downstream body and are discharged in the combined wake.

Apparently those FSI algorithms relying on a body-fixed non-inertial reference are not applicable to this type of complexbody–body interaction problems. These problems also pose major challenges to algorithms utilizing deforming meshes,especially, when the gaps between the bodies vary in a wide range of length scale. On the other hand, this case demonstratesthe simplicity and efficiency of the present direct forcing immersed boundary approach: compared with the single-bodyproblems in the previous sections, this case only combines the input parameters of the two bodies and then runs in astraightforward manner.

4.4. Fluttering and tumbling rectangular plates

Andersen et al. [1] studied the quasi-2D aerodynamics of freely falling rectangular plates at intermediate Reynolds num-bers in a water tank. In this section, two cases from [1], one fluttering plate at Re = 1147 and one tumbling plate at Re = 837,

5052 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

are selected to validate our FSI algorithm. For a freely falling plate in 2D, three dimensional parameters have to be specified,i.e., the plate thickness, H, the plate width, L, and the plate density, qs. In the experiments, aluminum plates were used andqs = 2700 kg/m3. The width-to-thickness ratio, K = L/H, is one major non-dimensional parameter for the plate. The momentof inertia per unit length of the plate is given by Is = qsHL(H2 + L2)/12. In [1], it was non-dimensionalized using the moment ofinertia per unit length about the axis of a circular cylinder of density qf and diameter L, i.e., pqfL

4/32. The density and thekinematic viscosity of the fluid (water) are qf = 997 kg/m3 and m = 8.8967 � 10�7 m2/s, respectively. The gravitational accel-eration is g = 9.8 m/s2. An estimate of the average descend speed using the buoyancy-corrected gravity can be used as thevelocity scale, together with the plate width L as the length scale, in the definition of the Reynolds number. But the measuredaverage value was used in [1]. Depending on these three non-dimensional parameters and initial conditions, a freely fallingplate may experience fluttering (swinging from side to side), tumbling (rotating and drifting obliquely), or apparently chaoticmotion with mixed fluttering and tumbling.

4.4.1. A fluttering rectangular plateFor the fluttering plate at Re = 1147, the plate thickness is H = 8.1 � 10�4 m, the width-to-thickness ratio K = 14. The

dimensionless moment of inertia is I⁄ = 0.16. The computational domain is [�0.13 m,0.07 m] � [�0.31 m,0.01 m] in the hor-izontal and vertical directions, respectively. The plate is centered at (0,0) and given no initial velocity at time t = 0 s. The ini-tial release angle was not reported in [1]. Here, an angle of 60� with the horizontal direction is chosen in the simulation. Forsimplicity, a uniform grid of 2500 � 4000 is adopted and the size of a grid cell is about 0.1H � 0.1H. A constant time stepDt = 3 � 10�5 s is used, which gives a maximum CFL number around 0.25. It should be noted that the resulting grid resolu-tion is only marginal for resolving the thin plate used in the experiment (about 10 cells across the thickness direction) andfor resolving the boundary layers at such a Reynolds number. Nevertheless, the results are satisfactory and the essential vor-tex wake structures are well captured, as will be shown below.

The plate develops a periodic fluttering motion very soon after it is released with the current initial angle. Fig. 24 showsthe trajectory of the fluttering plate. The measured trajectory from [1] is shown in (a) for comparison. Except the initial tran-sient, the overall agreement is very good. The plate swings from side to side during the free fall under gravity. It glides from aturning point on one side with low angle of attack (the angle between the direction along the plate width and its transla-tional velocity vector at the center of mass). Near the end of its gliding, the center of mass starts to raise until the plate

Fig. 24. Trajectory of a freely falling rectangular plate fluttering in a fluid at Re = 1147: (a) the simulated and measured [1] trajectories, shown as the platecross-sections, are given in green with time intervals of 0.03 s and in black with time intervals of 0.04 s, respectively; (b) the simulated trajectory of thecenter of mass and the overall vortex wake structure at t = 2.4 s. (For interpretation of the references to colour in this figure legend, the reader is referred tothe web version of this article.)

Fig. 25. Simulated (dashed line) and measured (solid line) velocities [1] as functions of time for a freely falling rectangular plate fluttering in a fluid atRe = 1147: (a) the horizontal velocity component ux; (b) the vertical velocity component uy; and (c) the angular velocity x.

Table 1Average translational velocities h uxi and huyi, angular velocity hxi for the fluttering plate withK = 14 at Re = 1147.

huxi(m/s) h uyi(m/s) hxi (rad/s)

Experiment [1] 0.198 ± 0.003 �0.090 ± 0.002 6.8 ± 0.1Simulation 0.195 �0.093 6.9

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5053

reaches the turning point on the other side. Near the turning point, the magnitude of the angular velocity of the plate is largeand the plate shows a rapid rotation until reaching the maximum angle of attack.

As shown in Fig. 25, the linear and angular velocities of the plate from the simulation match the experimental dataremarkably well. It is interesting that in such a complex fluttering motion the vertical velocity can be described very wellas a sinusoidal curve. Also, its oscillatory component has a frequency of twice that of other two velocities due to theside-to-side periodic motion pattern. The average velocities from the experiment and the present simulation are listed inTable 1. The simulation data are averaged from three complete periods after the initial transient (t > 1 s). For the horizontaland angular velocities, the data listed here are the averages of their absolute values, i.e., as defined in [1]. In general, the pres-ent data are within the experimental deviations. From Fig. 25, it is evident that, first, the center of mass raises to the turningpoint, at which the vertical velocity is zero; then the plate starts to descend again, but its center of mass keeps the sidewardmotion before the horizontal velocity reaches zero; after the horizontal velocity reverses its sign, the plate begins its swingmotion toward the other side, whereas the angle of attack is still increasing during the above process until the angular veloc-ity also changes its sign at the position of the maximum angle of attack.

Fig. 24(b) shows the overall vortex wake structure at t = 2.4 s. The vortices formed from the wake instability following theplate is very organized. Away from the plate, especially near the turning points, the vortices interact with each other andproduce a very complex pattern of vortex wake structure. Furthermore, a series of detailed vorticity fields around the flut-tering plate are given in Fig. 26 to illustrate the vortex dynamics of the complex fluttering motion of a rectangular plate. Itshould be noted that in [1] they only considered a plate with an elliptical cross-section at a much lower Reynolds number intheir 2D conformal mapping calculation, and the vortex wake structures from such conditions are not the same as the actual

Fig. 26. Snapshots of the vorticity field for a fluttering plate at Re = 1147: (a) t = 1.500 s; (b) t = 1.530 s; (c) t = 1.560 s; (d) t = 1.590 s; (e) t = 1.605 s; (f)t = 1.620 s; (g) t = 1.635 s; (h) t = 1.650 s; (i) t = 1.665 s; (j) t = 1.680 s; (k) t = 1.695 s; (l) t = 1.710 s; (m) t = 1.740 s; (n) t = 1.770 s; (o) t = 1.800 s. And, (a)downward gliding with wake instability and vortex formation; (b) and (c) the plate pitches up; (d)–(f) fast rotation and shedding of the leading and trailingedge vortices; (g)–(h) the plate is around the turning point; (i)–(k) a vortex dipole is generated and the plate resumes gliding; and (l)–(o) a new period withwake instability and vortex shedding.

5054 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

phenomena in the experiments and what revealed here using the present approach with the same experimental parameters.Here, the series starts from the downward gliding motion in Fig. 26(a). In general, during gliding, the boundary layer devel-oped along the lower surface of the plate is thinner than that on the upper surface due to the downward motion from thegravitational acceleration. The vortices shed from the lower surface is much stronger than those from the upper surface. Theimbalance of the vortex shedding produces a net moment that pitches up the plate during the gliding motion. The alternatingvortex shedding from the two plate surfaces generates an oscillation of fluid moment on the plate, which is clearly illustratedby the oscillations of the angular velocity before its maximum magnitudes are reached as shown in Fig. 25(c). In addition,

Fig. 27. Trajectory of a freely falling rectangular plate tumbling in a fluid at Re = 837: (a) the simulated and measured [1] plate trajectories, shown as theplate cross-sections, are given in green (with time intervals of 0.03 s) and in black, respectively, and (b) the simulated trajectory of the center of mass andthe overall vortex wake structure at t = 1.77 s. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version ofthis article.)

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5055

due to the sharp plate corner, the lower surface boundary layer starts to separate very soon, as evident in Fig. 26(b), and thena shear layer instability develops into vortex shedding (Fig. 26(c)). As the angle of attack increases, the upper surface bound-ary layer starts to separate too (Fig. 26(d)). At this time, the plate has only very small translational velocities but a very large

5056 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

angular velocity, hence the vortices generated during the upward motion start to detach from the lower surface (Fig. 26(e)).These vortices become totally detached very soon and the upper surface boundary layer becomes fully separated due to therotational motion (Fig. 26(f)). Then the plate starts to descend and reaches the maximum angle of attack; in the mean time, avortex dipole originated from the lower surface sheds at the upper end of the plate (Fig. 26(g)–(k)). Finally, the plate resumesits gliding motion and repeats the above procedure, though with reversed directions of the horizontal and angular velocities(Fig. 26(l)–(o)).

This simulation was run on an IBM cluster equipped with 4.7 GHz Power6 processors using one 32-processor node. Withthe simultaneous multi-threading (SMT) feature enabled 64 processes were actually used. In average each time step needed1.157 s of CPU time, in which the four major steps of the basic fluid flow solver, i.e., the right-hand-side computation for themomentum equations, the parallel tridiagonal solver, the pressure Poisson solver, and the velocity corrector, used 0.021,0.080, 0.721, and 0.040 s, respectively. The latter three steps were included within the predictor–corrector loop in the pre-vious approach [33], but are moved out of the loop in the present approach and only executed once just like the first step. Thetotal wall time for this particular simulation (up to t = 3 s) was 33 h. With the present FSI algorithm around 5 corrector stepswere used in each time step. If the previous strong coupling algorithm is used with the same computational hardware, thetotal wall time will be more than 170 h (more than 1 week), since those three steps have to be carried out 6 more timesincluding the predictor and finalizer steps in the predictor–corrector loop. Therefore, the savings of computational cost byusing the present FSI algorithm is remarkable and substantial. On the other hand, it should be noted that in the present pre-dictor–corrector algorithm several steps of the immersed boundary treatment are carried out only once too. These steps,

Table 2Average linear velocities huxi and huyi, angular velocity h xi for the tumbling plate with K = 8 at Re = 837.

huxi(m/s) h uyi(m/s) hxi (rad/s) Descent angle (degree)

Experiment [1] 0.159 ± 0.003 �0.115 ± 0.005 14.5 ± 0.3 35.8 ± 1.3Simulation 0.149 �0.104 15.4 34.9

Fig. 28. Simulated velocities as functions of time for a freely falling rectangular plate tumbling in a fluid at Re = 837: (a) the horizontal velocity componentux; (b) the vertical velocity component uy; and (c) the angular velocity x.

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5057

including the preliminary velocity predictor, the fulfillment of the forcing term in the momentum equations, and the fieldextension operation, took about 0.013 s. However, the positions of the immersed bodies still have to adjusted multiple timesin one time step until the fluid–structure coupling process is converged. The re-positioning of the immersed bodies and thefollowing immersed boundary setup procedure (grid point classification, closest surface point computation, interpolationstencil construction, etc.) could generate major computational cost. Even for the simple 2D geometry in this case, 0.192 sof CPU time was spent on this part since the relocation/setup procedure for the immersed bodies was repeated 7 times.Of course, the momentum forcing construction step was also repeated 7 times and took 0.088 s in total, whereas the ODE

Fig. 29. Snapshots of the vorticity field for a tumbling plate at Re = 837: (a) t = 1.170 s; (b) t = 1.194 s; (c) t = 1.206 s; (d) t = 1.218 s; (e) t = 1.230 s; (f)t = 1.242 s; (g) t = 1.254 s; (h) t = 1.266 s; (i) t = 1.278 s; (j) t = 1.290 s; (k) t = 1.302 s; (l) t = 1.314 s; (m) t = 1.326 s; (n) t = 1.338 s; (o) t = 1.362 s. And, (a) and(b) downward gliding with wake instability and vortex formation; (c) and (d) the plate pitches up; (e) and (g) fast rotation and shedding of the leading andtrailing edge vortices; (h)–(j) the plate is completing its 180� turn; (k)–(m) a vortex dipole is generated and the plate resumes gliding; and (n) and (o) a newperiod with wake instability and vortex shedding.

5058 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

solver only took 0.001 s that is negligible compared with other parts. It should be noted that the simplified field extensionalso contributes to the computational cost reduction too. For instance, the pressure and velocity points requiring field exten-sion operation in the previous approach are not needed any more here; thus the CPU time to be spent on the point identi-fication and the extrapolation stencil construction can be saved. In this 2D case, this portion of CPU time was about 0.002 s;but it could take a larger portion of the total computational cost for 3D complex geometries and high-resolution grids. More-over, the surface element integration approach would be more expensive than the present point integration approach forforce and moment evaluation. Nevertheless, a robust and efficient setup procedure for complex immersed geometries is cru-cial for the overall efficiency of FSI simulations using the immersed boundary methods [35].

4.4.2. A tumbling rectangular plateThe tumbling plate at Re = 837 in [1] has the same plate thickness H = 8.1 � 10�4 m as the fluttering plate, but a much

smaller width and the width-to-thickness ratio is K = 8. The dimensionless moment of inertia is I⁄ = 0.29. The computationaldomain is [�0.03 m,0.24136 m] � [�0.23 m,0.01 m] in the horizontal and vertical directions, respectively. A uniform grid of3392 � 3000 is adopted and the size of a grid cell is about 0.1H � 0.1H, same as the fluttering plate case. The plate is centeredat (0,0) and no initial velocity is given at time t = 0 s. The initial release angle is the same as reported in [1]: 45� with thehorizontal direction. A constant time step Dt = 3 � 10�5 s is used, which gives a maximum CFL number about 0.27. Althoughthe Reynolds number is a little lower than that of the fluttering plate, the magnitude of the instantaneous velocity vector ishigher in this case as evident by the maximum CFL number. Therefore, the grid resolution for this case is also marginal as thefluttering plate case.

Fig. 27 shows the trajectory of the tumbling plate. After released from its initial position, the plate starts to descend androtate; soon it pitches up just like the fluttering plate. Due to the large angular velocity, the plate rotates more than 135�clockwise from its initial release angle. Because of a large restoring moment, the plate does not glide toward the left sidethough. Instead, it slightly rotates anti-clockwise and dives almost vertically. Along the dive, the plate drifts toward the rightside and finally pitches up as in the fluttering case. A hook-like trajectory is formed in the above procedure. At the end of the‘‘hook’’, the center of mass elevates and reaches the turning point on the right side, similar to the fluttering plate. However,near a turning point, the tumbling plate has such a larger angular velocity that it makes a complete 180� rotation and thelower surface during the gliding becomes upward facing. The plate then resumes its sideward gliding and starts a new periodof the tumbling motion. The measured trajectory from [1] is shown in Fig. 27(a) for comparison. In the simulation, the trans-lational velocities are under-predicted and the plate has a smaller angle of descent and shorter periods for the first few cy-cles. After that, the plate gains velocities very close to these from the experiment and the displaced trajectory from thesimulation matches the experimental trajectory very well.

The average velocities from the experiment and the present simulation are listed in Table 2. The simulation data are aver-aged from three complete periods after the initial transient (t > 1 s). The average translational velocities from the simulationare slightly lower than those from the experiment, whereas the average angular velocity is slightly higher than the exper-imental value. However, the angle of descent is very close to the experimental one and within the experimental deviation.The three velocities from the simulation are shown in Fig. 28. Unlike the fluttering plate, all three velocities have the samefrequency during the tumbling motion. Also, the tumbling plate rotates much faster near the turning points with angularvelocity values around 30 rad/s.

Fig. 27(b) also shows the overall vortex wake structure at t = 1.77 s. Due to the oblique motion, the vortex–vortex inter-action is not as strong as that in the fluttering plate case. Similarly, a series of detailed vorticity fields around the tumblingplate are given in Fig. 29. In [1] they only considered a plate with an elliptical cross-section at a different Reynolds number. Ingeneral, the vortex shedding during gliding is very similar to that of the fluttering plate. However, the vortices generatedfrom the wake instability in the tumbling plate case are weaker in strength and less in number than those from the flutteringplate due to a shorter gliding path and a smaller gliding speed at each gliding cycle.

5. Conclusions and future work

A direct forcing immersed boundary framework has been presented for simulating strongly coupled FSI problems. Thisframework is based on the immersed boundary method developed by Yang and Balaras for moving boundary problems in[31] and the strong coupling scheme for FSI problems by Yang et al. in [33]. However, the field extension strategy in [31]has been greatly simplified in this study. The complicated geometric procedures for identifying pressure and velocity pointsthat require field extension operation are totally eliminated. Based on a four-step fractional-step method, the present fieldextension is centered on tackling the pressure gradient components that are collocated with the staggered velocity compo-nents. The reformulated field extension approach for moving boundary problems is very simple, robust, and efficient, but theoverall accuracy of the original method is not altered. Furthermore, the strong coupling scheme in [33] has been significantlyexpedited by several times or one order of magnitude. With the reformulated field extension procedure, the fluid force andmoment exerted on the immersed bodies can be evaluated by a simple pointwise integration of the momentum forcing term.And this step usually is applicable before solving the expensive pressure Poisson equation in a direct forcing framework.Therefore, in this study the fluid solver has been moved out of the predictor–corrector iterative loop, which now contains

Fig. 30. Simulated (dashed line) and measured (solid line) velocities in the gravitational direction for a settling sphere released with zero initial velocity in afluid at rest. The computational domain is [�10D,54D] � [�16D,16D] � [�16D,16D] (in the gravitational direction, and two horizontal directions,respectively) and the grid is 480 � 192 � 192 with uniform grid spacing of 0.03D around the sphere. The experiment was performed in a water tank with asteel sphere of 1 mm diameter (qs = 7850 kg/m3) by Mordant and Pinton [19]. The average terminal speed from the experiment is 0.383 m/s and thecorresponding Reynolds number is about 430. The present simulation gives an average terminal speed of 0.386 m/s, which is 1% higher than theexperimental result. The faster acceleration before reaching the terminal speed was also observed in the immersed boundary simulation of Kim and Choi[14] on a Cartesian grid with a resolution around the sphere similar to the present case.

Fig. 31. Instantaneous vortical structures for a settling sphere released with zero initial velocity in a fluid at rest: (a) t = 0.15 s; (b) t = 0.30 s; (c) t = 0.45 s;(d) t = 0.60 s; (e) t = 0.75 s; and (f) t = 0.90 s.

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5059

only a few modules related to the immersed boundary setup procedure. This step greatly simplifies the code structure andalleviates the memory requirement, but the strong coupling property of the FSI algorithm is not weakened at all.

To verify the reformulated field extension procedure, two cases used in [31], the in-line oscillation of a circular cylinder ina fluid at rest and a transversely oscillating circular cylinder in a free-stream, have been performed. Then the simulation of ahovering elliptical wing with a complicated motion has been carried out too. The results of these cases match the referencedata very well, which shows the accuracy and effectiveness of the present simplified field extension strategy.

Furthermore, a variety of strongly coupled FSI problems have been simulated with the present direct forcing immersedboundary framework. The vortex-induced vibrations of a circular cylinder have been studied and, especially, one very lowmass ratio case has been demonstrated to show the desirable property of numerical stability of the present FSI algorithm.In addition, as a different type of flow-induced vibrations, galloping of rectangular bodies has been studied too. These casesinclude a transversely galloping body of a square cross-section in a free-stream, a rotationally galloping body of a rectangularcross-section in a free-stream, and a multi-body case with both bodies in tandem arrangement. The simulation data from

5060 J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061

these fluid-induced vibration cases have been compared with the results from high-order spectral element simulations uti-lizing a body-fitted non-inertial reference frame. The excellent agreement between the present results on non-boundary-conforming Cartesian grids and those on body-fitted grids is truly remarkable.

On the other hand, the performance of the present algorithm in some real world challenges is even more noteworthy. Inthis study, fluttering and tumbling rectangular plates have been simulated and compared with the experimental data fromwater tank tests. The fluttering and tumbling motions of thin plates (less than 1 mm thick) at moderate Reynolds numbers of103 are very sensitive to initial conditions and small geometric features. And surprisingly, the trajectories (including theplate orientation) of the plate on Cartesian grids of marginal resolution match those from the experiments very well. None-theless, the present simulations have also provided detailed vortex wake structures, which can be used to further study theFSI mechanisms in these complicated fluttering or tumbling motions.

In summary, these FSI cases of both prescribed and predicted motions demonstrate the accuracy, simplicity, and effi-ciency of the new method and its applicability in a wide range of complicated FSI problems. However, there are still severalissues to be addressed in the future. First, the present immersed boundary method employs a linear local reconstruction thatlimits its applications in high Reynolds number flows due to the prohibitive computational cost of resolving a high Reynoldsnumber turbulent boundary layer with a Cartesian grid. Therefore, it is imperative to develop immersed boundary based wallmodels for the local reconstruction in high Reynolds number applications. Some preliminary work toward this direction hasbeen reported in [4] and further work will be conducted to develop a robust and accurate wall model for immersed boundarymethods. Interestingly, the present FSI algorithm can be combined with a wall model based on the direct forcing approach ina straightforward manner. This promises the non-compromised efficiency of the present algorithm in high Reynolds numberFSI problems. Second, the single block grid used for the flow solver restricts the grid resolution available near the immersedboundaries. Actually, this is a limitation for the Navier–Stokes solver, instead of one for the present FSI algorithm. Local gridrefinement techniques can be used to reduce the overall computational cost and increase the resolution near the immersedboundaries. For instance, in the simulations of fluttering and tumbling plates, a grid can be locally refined along the path ofthe plate and even be coarsened at a region away from the plate and its wake without decreasing the overall accuracy.

In addition, for applications with a single body in an (approximately) unbounded domain, the adoption of a non-inertialreference frame that partially or completely (one up to six DOFs) follows the body motion can greatly reduce the size of thecomputational domain and the total cost of the computation. For example, a settling sphere travels a very long distance be-fore it reaches the terminal settling speed. Apparently the computation using the present method will be very expensive.Actually, a simply implemented option that allows the grid to follow the translational motion of the sphere has been addedin the code. This modification is similar to the immersed boundary method used in [14] and some moving grid methods byadding a grid velocity in the convection terms. The speed evolution of the sphere case are shown in Fig. 30 and the instan-taneous vortical structures at several instances are shown in Fig. 31. It turns out that the sphere travels a distance of morethan 40D before it is about to settle at its terminal speed. And in order to obtain a stable estimate the simulation has to beperformed up to t = 0.60 s with a total distance nearly 200D. The complete run will be prohibitively expensive for a simula-tion using a stationary grid even with a local grid refinement technique. Therefore, this type of moving grid formulation canbe very useful in more complicated cases such as ship hydrodynamics, where a ship sustains motions of multiple DOFs undera variety of environmental conditions due to the currents, waves, and winds around the ship. Obviously a stationary gridcannot satisfy the requirement of this type of dynamic simulations. However, if the grid follows all the six DOFs, it will alsobe difficult to guarantee good grid resolution, for instance, for the air–water interface when there is a large angle pitchingmotion. Thus the grid can follow the surge (forward/backward translation), sway (sideward translation), and yaw (rotationabout the vertical axis) of the ship, and maintain a good resolution for the ship boundary layers, air–water interface, and theship wake at the same time.

Acknowledgments

This work was sponsored by the Office of Naval Research (ONR) under Grant N000141-01-00-1-7, with Dr. Patrick Purtellas the program manager. The simulations presented in this paper were performed at the Department of Defense (DoD)Supercomputing Resource Centers (DSRCs) through the High Performance Computing Modernization Program (HPCMP).The authors are grateful to Dr. K.W. Cheng for proofreading the initial manuscript.

References

[1] A. Andersen, U. Pesavento, Z.J. Wang, Unsteady aerodynamics of fluttering and tumbling plates, Journal of Fluid Mechanics 541 (2005) 65–90.[2] E. Balaras, Modeling complex boundaries using an external force field on fixed cartesian grids in large-eddy simulations, Computers & Fluids 33 (3)

(2004) 375–404.[3] R.M. Beam, R.F. Warming, An implicit finite-difference algorithm for hyperbolic systems in conservation-law form, Journal of Computational Physics 22

(1) (1976) 87–110.[4] S. Bhushan, P. Carrica, J. Yang, F. Stern, Scalability studies and large grid computations for surface combatant using CFDShip-Iowa, International Journal

of High Performance Computing Applications 25 (4) (2011) 466–487.[5] H.M. Blackburn, G.E. Karniadakis, Two and three-dimensional simulations of vortex-induced vibration of a circular cylinder, in: Third International

Offshore and Polar Engineering Conference, 1993, pp. 715–720.[6] H. Choi, P. Moin, Effects of the computational time step on numerical solutions of turbulent flow, Journal of Computational Physics 113 (1) (1994) 1–4.

J. Yang, F. Stern / Journal of Computational Physics 231 (2012) 5029–5061 5061

[7] H. Dütsch, F. Durst, S. Becker, H. Lienhart, Low-Reynolds-number flow around an oscillating circular cylinder at low Keulegan–Carpenter numbers,Journal of Fluid Mechanics 360 (1998) 249–271.

[8] E.A. Fadlun, R. Verzicco, P. Orlandi, J. Mohd-Yusof, Combined immersed-boundary finite-difference methods for three-dimensional complex flowsimulations, Journal of Computational Physics 161 (1) (2000) 35–60.

[9] R. Falgout, J. Jones, U. Yang, The design and implementation of hypre, a library of parallel high performance preconditioners, in: A. Bruaset, A. Tveito(Eds.), Numerical Solution of Partial Differential Equations on Parallel Computers, vol. 51, Springer-Verlag, 2006, pp. 267–294 (UCRL-JRNL-205459).

[10] T.I. Fossen, Guidance and Control of Ocean Vehicles, John Wiley & Sons, New York, NY, 1994.[11] D. Goldstein, R. Handler, L. Sirovich, Modeling a no-slip flow boundary with an external force field, Journal of Computational Physics 105 (2) (1993)

354–366.[12] E. Guilmineau, P. Queutey, A numerical simulation of vortex shedding from an oscillating circular cylinder, Journal of Fluids and Structures 16 (6)

(2002) 773–794.[13] R.W. Hamming, Stable predictor–corrector methods for ordinary differential equations, Journal of the ACM 6 (1959) 37–47.[14] D. Kim, H. Choi, Immersed boundary method for flow around an arbitrarily moving body, Journal of Computational Physics 212 (2) (2006) 662–680.[15] J. Kim, D. Kim, H. Choi, An immersed-boundary finite-volume method for simulations of flow in complex geometries, Journal of Computational Physics

171 (1) (2001) 132–150.[16] M.-C. Lai, C.S. Peskin, An immersed boundary method with formal second-order accuracy and reduced numerical viscosity, Journal of Computational

Physics 160 (2) (2000) 705–719.[17] N. Mattor, T.J. Williams, D.W. Hewett, Algorithm for solving tridiagonal matrix problems in parallel, Parallel Computing 21 (11) (1995) 1769–1782.[18] J. Mohd-Yusof, Combined immersed-boundary/B-spline methods for simulations of flow in complex geometries, in: Annual Research Briefs, Center for

Turbulence Research, Stanford University, Stanford, CA, 1997, pp. 317–327.[19] N. Mordant, J.-F. Pinton, Velocity measurement of a settling sphere, The European Physical Journal B: Condensed Matter and Complex Systems 18

(2000) 343–352.[20] C.S. Peskin, Flow patterns around heart valves: a numerical method, Journal of Computational Physics 10 (2) (1972) 252–271.[21] C.S. Peskin, Numerical analysis of blood flow in the heart, Journal of Computational Physics 25 (3) (1977) 220–252.[22] I. Robertson, L. Li, S.J. Sherwin, P.W. Bearman, A numerical study of rotational and transverse galloping rectangular bodies, Journal of Fluids and

Structures 17 (5) (2003) 681–699.[23] E.M. Saiki, S. Biringen, Numerical simulation of a cylinder in uniform flow: application of a virtual boundary method, Journal of Computational Physics

123 (2) (1996) 450–465.[24] P.N. Swarztrauber, A direct method for the discrete solution of separable elliptic equations, SIAM Journal on Numerical Analysis 11 (6) (1974) 1136–

1150.[25] Y.-H. Tseng, J.H. Ferziger, A ghost-cell immersed boundary method for flow in complex geometry, Journal of Computational Physics 192 (2) (2003)

593–623.[26] M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, Journal of Computational Physics 209 (2)

(2005) 448–476.[27] M. Vanella, E. Balaras, A moving-least-squares reconstruction for embedded-boundary formulations, Journal of Computational Physics 228 (18) (2009)

6617–6628.[28] M. Vanella, T. Fitzgerald, S. Preidikman, E. Balaras, B. Balachandran, Influence of flexibility on the aerodynamic performance of a hovering wing, Journal

of Experimental Biology 212 (1) (2009) 95–105.[29] Z.J. Wang, Two dimensional mechanism for insect hovering, Physical Review Letters 85 (10) (2000) 2216–2219.[30] S. Xu, Z.J. Wang, An immersed interface method for simulating the interaction of a fluid with moving boundaries, Journal of Computational Physics 216

(2) (2006) 454–493.[31] J. Yang, E. Balaras, An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries, Journal of

Computational Physics 215 (1) (2006) 12–40.[32] J. Yang, S. Preidikman, E. Balaras, A strong coupling scheme for fluid–structure interaction problems in viscous incompressible flows, in: Proceedings of

the International Conference on Computational Methods for Coupled Problems in Science and Engineering, Santorini Island, Greece, May 2005.[33] J. Yang, S. Preidikman, E. Balaras, A strong-coupled, embedded-boundary method for fluid-structure interactions of elastically mounted rigid bodies,

Journal of Fluids and Structures 24 (2) (2008) 167–182.[34] J. Yang, F. Stern, Sharp interface immersed-boundary/level-set method for wave-body interactions, Journal of Computational Physics 228 (17) (2009)

6590–6616.[35] J. Yang, F. Stern, Robust and efficient setup procedure for complex triangulations in immersed boundary simulations, in: Proceedings of the ASME 2012

Fluids Engineering Summer Meeting, Rio Grande, Puerto Rico, July 2012.[36] X. Yang, X. Zhang, Z. Li, G.-W. He, A smoothing technique for discrete delta functions with application to immersed boundary method in moving

boundary simulations, Journal of Computational Physics 228 (20) (2009) 7821–7836.[37] N. Zhang, Z. Zheng, An improved direct-forcing immersed-boundary method for finite difference applications, Journal of Computational Physics 221 (1)

(2007) 250–268.

(PDF) A simple and efficient direct forcing immersed boundary framework for fluid–structure interactions - DOKUMEN.TIPS (2024)

FAQs

What is simple immersed boundary method? ›