A conservative viscous vorticity method for unsteady unidirectional and oscillatory flow past a circular cylinder

Ocean Engineering 191 (2019) 106504

Contents lists available at ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

A conservative viscous vorticity method for unsteady unidirectional and oscillatory flow past a circular cylinder Chunlin Wu 1 , ∗, Spyros A. Kinnas 2 , Zhihao Li 3 , Yunjun Wu 3 Ocean Engineering Group, Department of Civil, Architectural and Environmental Engineering, University of Texas at Austin, 301E E Dean Keeton St C1700, Austin, TX 78712, USA

ARTICLE

INFO

Keywords: Vorticity equation Impulsively started flow Oscillatory flow 2-D unsteady

ABSTRACT A conservative vorticity method for solving the 2-D incompressible unsteady viscous flow in the eulerian frame is proposed. The numerical method implements the vorticity-stream function formulation of the Navier–Stokes (N–S) equations, solved by the cell-centered Finite Volume Method (FVM). A convection limiter, coupling the first order upwind scheme and the Quadratic Upstream Interpolation for Convective Kinematics (QUICK) scheme, is introduced to maintain the accuracy and keep the monotonicity of vorticity interpolations. The vorticity induced stream functions are recovered by a closed-form expression, which also accurately predicts the far-field boundary conditions. The convective flux is evaluated by the difference of stream functions on the neighboring two nodes, which spontaneously eliminates the net mass flux out of a control cell. The conservative FVM and QUICK scheme and the definite evaluations of stream functions and convective flux guarantee the conservation of circulation and mass. The impulsively started flow past a circular cylinder at Reynolds numbers ranging from 100 to 9500 and the inline oscillatory flow around a circular cylinder at Reynolds numbers 100 and 200 are numerically simulated using this method. The current results under the inertial and non-inertial translational frame are well validated by comparisons with other numerical and experimental benchmark tests.

1. Introduction The purpose of this piece of work is to present an eulerian VIScous Vorticity Equation (VISVE4 ) method for 2-D incompressible unsteady viscous flow past a circular cylinder and the numerical results obtained by this method. The vorticity methods evolve from the traditional computational methods for fluid dynamics into a distinguishing branch of approaches. Generally, these methods are classified into two types according to the frame of reference. The first class is known as the mesh-free Lagrangian methods, which are also called the vortex ‘‘blob’’ methods or the discrete vortex methods (Chorin and Bernard, 1973). In these methods, the motions of vortex ‘‘blobs’’ are described by a set of ordinary differential equations via splitting the convection and diffusion (Chorin, 1973) so that the discrete point vortices are individually traced in time. Special treatments, including the vortex merging algorithms (Nielsen and Schwind, 1971; Moore, 1974) or the artificial viscosity corrections (Chorin, 1973; Krasny, 1987), are essential to avoid the singularity of point vortices and to keep the interactions

between neighboring vortices bounded. However, several problems exist within the discrete vortex methods. Firstly, an increasing number of particles is inevitable to capture the small scale structures developed over time, which will impose great challenges on the computational resources. In addition, local refinements are necessary to keep the point vortices densely distributed to maintain the accuracy during longtime simulations (Saffman and Meiron, 1986). The growing spatial resolution will eventually be out of control and limit the affordable simulation time, which triggers the blooms of certain fast algorithms (Lindsay and Krasny, 2001; Ploumhans et al., 2002). Furthermore, although the discrete vortex methods have been extensively used to study the development of vortex structures in the inviscid flow, the incorporation of viscous diffusion is relatively complicated (Leonard, 1980). Last but not least, the distance between two point vortices has to be limited to a certain range to have a smooth velocity field, which could be difficult to accomplish for 3-D extensions of these methods. The other class is named as the mesh based Eulerian methods, in which the governing equations can be numerically solved by methods

∗ Corresponding author. E-mail address: [email protected] (C. Wu). 1 Ph.D. Candidate. 2 Professor. 3 Graduate student. 4 The present method is referred to as VISVE in this study for the sake of compact documentation.

https://doi.org/10.1016/j.oceaneng.2019.106504 Received 9 June 2019; Received in revised form 28 August 2019; Accepted 28 September 2019 Available online 3 October 2019 0029-8018/© 2019 Elsevier Ltd. All rights reserved.

Ocean Engineering 191 (2019) 106504

C. Wu et al.

The impulsively started flow, both unidirectional and oscillatory, past a circular cylinder is a challenging and substantial topic. Especially for numerical simulations, the methods should resolve both the early abrupt initiations and the subsequent separations. The interests in numerical computations were started by Payne (1958), following whom many methods have been proposed to solve this problem. The discrete vortex methods were popular in predicting flow past a cylinder (Chang and Chern, 1991; Smith and Stansby, 1988; Koumoutsakos and Leonard, 1995) for some time, but one major disadvantage of these methods is the ever-increasing number of particles, which limits the affordable simulation time. The N–S methods were also widely used to study these problems (Dennis and Chang, 1970; Nieuwstadt and Keller, 1973; Liu et al., 1998; Wright and Smith, 2001; Calhoun, 2002; Russell and Wang, 2003; Chiu et al., 2010; Pier, 2002; Sipp and Lebedev, 2007; Marquet et al., 2008), whose topics vary among flow regimes, body forces, and wake stability analyses. However, the N–S based schemes usually require very large computational domains to guarantee unperturbed far-field boundary conditions, leading to a large number of cells in the grid. Lately, researchers began to show interests in using the lattice Boltzmann method to model the flow past cylinders (Mei and Shyy, 1998; Guo and Zhao, 2003; Imamura et al., 2005; Dupuis et al., 2008; Cui et al., 2017; Hejranfar and Ezzatneshan, 2014). Most of the work only considers flow past a fixed cylinder at relatively low Reynolds numbers (≤200), yet there are still doubts about their capability of resolving the impulsive start process at higher Reynolds numbers.

such as Finite Difference Method (FDM) or FVM. Whether the velocity or the stream function is used to represent the kinematics of the flow further categorizes these methods into the vorticity–velocity formulation and the vorticity-stream function formulation. The velocity field and the vorticity field are usually associated by a Poisson equation (Wu et al., 1995). Remarkably, Wu and Thompson (1973) proposed an integro-differential scheme that enables calculating the velocity by a direct and explicit integral of the vorticity. Since a known vorticity field uniquely determines the entire flow field (Munson et al., 2014), the simulation of vorticity dynamics and the calculation of induced velocity can be separated into two consecutive steps, that involve only the portion of grids with non-zero vorticity. The biggest advantage of this integro-differential method is that the numerical simulation is confined to only the regions of interests, and meanwhile the number of unknowns remains under control, which further leads to a dramatic drop in the computational time and allows us the luxury of conducting longtime simulations. In the current work, the integro-differential technique is reformulated to fit in with the vorticity-stream function formulation of the N–S equations, whose number of unknown variables reduces to two. The validations of the proposed vorticity method under the inertial frame are placed on predicting the unsteady viscous impulsively started flow past a circular cylinder. The flow past an impulsively started cylinder is a well-established prototype for unsteady flow since it reveals the structural characteristics of time-dependent flow with separations and has triggered more than abundant amounts of academic investigations. This problem imposes great difficulties on not only the experimental studies but also the numerical approaches because of the uniformly sudden beginning of the flow field. The experiments should be properly designed to accommodate for the instantaneous launch of a cylinder in the quiescent flow (Coutanceau and Bouard, 1977a,b). The numerical schemes must resolve the initial vorticity generation to accurately model the sharp start and precisely catch the following development of topological structures in the wake, formed by the interactions of multiple shedding vortices. A considerable number of numerical attempts have been successfully made and yield promising results such as the finite difference scheme of Anderson and Reider (1996), the discrete vortex methods of Koumoutsakos and Leonard (1995) and Subramaniam (1997), the vorticity–velocity formulation of Qian and Vezza (2001), and the lattice-Boltzmann application of Li et al. (2004). The above existing literature provides sufficient diagnostics, including equi-vorticity contours, boundary layer profiles, and body forces, to validate the current method. Additionally, the method is validated by simulating the unsteady viscous oscillatory flow around a circular cylinder in the non-inertial frame of reference. The inline oscillatory flow past a cylinder or the cylinder oscillating in the flow at rest has a variety of engineering applications, especially in the fields of civil and offshore constructions. The interests in this problem featured by complex vortex–body and vortex–vortex interactions arise from predicting the loads on structures subject to periodic waves or currents to optimize the design of structures (Chakrabarti, 1987). Experiments on this subject have been documented by Williamson (1985), Tatsuno and Bearman (1990), and Sarpkaya (1976, 2002), who have conducted comprehensive research about the flow regimes and instantaneous forces. Sustained numerical studies which successfully recovered some of the experimental observations can be found in the publications of Justesen (1991), Badr et al. (1995), Dütsch et al. (1998), Zhou and Graham (2000), Guilmineau and Queutey (2002), Nehari et al. (2004), Elston et al. (2006), and An et al. (2011). The majority of the numerical simulations are based on 2-D models until some recent investigations (Elston et al., 2006; Nehari et al., 2004; An et al., 2011) focused on the 3-D symmetry-breaking and instabilities experimentally visualized by Honji (1981). The results, such as velocity profiles, vorticity contours, and inline forces, from the present method are compared with those of established benchmark tests to assess the accuracy of this method.

In this work, a conservative numerical method that solves for the vorticity in Eulerian coordinates to simulate the incompressible unsteady viscous flow around bluff bodies in 2-D is proposed. This method was originally proposed by Tian and Kinnas (Tian, 2014; Tian and Kinnas, 2015). It is designed to be spatially compact and computationally efficient, and meanwhile capable of modeling the dynamics of vortices under both inertial and non-inertial frame. The integro-differential scheme is improved to be compatible with the vorticity-stream function formulation of the N–S equations, which are numerically solved by the cell-centered FVM. The mass flux of convection is calculated by the direct subtraction of the stream functions on adjacent cell nodes to automatically nullify the net flux of a control cell. Further, the stream function Poisson equation is solved based on the Green’s function, and an analytical expression is derived to compute the vorticity induced stream functions on the cell vertices exactly in 2-D. In addition, the amount of vorticity created is determined by enforcing the non-penetrating and no-slip boundary conditions on the wall. A convection limiter, coupling the first order upwind scheme and the QUICK scheme, to control the interpolations of vorticity is utilized in the current method to maintain the order of accuracy and resolve the sharp gradients of vorticity at the same time. The special treatments of explicitly recovering the stream functions and convective flux and the conservativity of FVM and QUICK schemes guarantee the divergencefree of velocity and preserve the vorticity throughout the numerical calculations. The remainder is organized as follows. The formulations and numerical schemes are explained in Section 2. Section 3 is devoted to showing the numerical results of impulsively started flow past a circular cylinder at Reynolds numbers ranging from 100 to 9500 to validate the proposed method under the inertial frame. Section 4 presents the results of oscillatory flow around a circular cylinder at Reynolds numbers 100 and 200 to validate this method in the non-inertial translating frame of reference. Section 5 gives comprehensive descriptions of the computational efficiency and parallelization. Finally, the achievements and future work are discussed in Section 6. 2

Ocean Engineering 191 (2019) 106504

C. Wu et al.

2. Formulations and numerical schemes 2.1. Governing equations Assuming the non-inertial frame only moves with a time-dependent velocity 𝑼 (𝑡), a translational acceleration term appears in the velocity– pressure formulation of Navier–Stokes equations ∇𝑝 𝜕𝒖 d𝑼 + (𝒖⋅𝛁)𝒖 + =− + 𝜈∇2 𝒖, 𝜕𝑡 d𝑡 𝜌

(1)

𝛁⋅𝒖 = 0,

(2)

Fig. 1. Flowchart of the general solving procedures.

where, 𝒖 = (𝑢, 𝑣) is the velocity under the non-inertial frame, 𝜈 is the kinematic viscosity of the fluids, 𝑝 is the pressure, 𝜌 is the density of the fluids, and 𝑡 is the time. The vorticity in the non-inertial frame of reference is the same as that under the inertial frame since the background flow 𝑼 is homogeneous in space 𝜔=

𝜕𝑣 𝜕𝑢 − . 𝜕𝑥 𝜕𝑦

(3)

It can be shown that by absorbing the transnational acceleration term into the pressure term as a potential in Eq. (1), the curl operation on this equation will result in a vorticity transport equation that takes the same form as the one in the inertial frame of reference 𝜕𝜔 + (𝒖⋅𝛁)𝜔 = 𝜈∇2 𝜔, (4) 𝜕𝑡 with the initial and boundary conditions 𝒖𝑡=0 𝒖 𝒖

= = =

𝑼, 𝑼 𝑏 , on body surface, 𝑼 , at infinity,

Fig. 2. The FVM control cell.

where, 𝛥𝑙𝑓 indicates the edge length, 𝑚𝑓 = ∫𝜕𝐴 (𝒖 ⋅ 𝒏)𝑑𝑙 is the mass flux 𝑓 across a cell edge, 𝐴 is the area of the control cell, 𝜕𝐴 is the boundary enclosing the cell, and 𝒏 is the unit normal vector of the boundary pointing outwards. Under the FVM discretization, boundary conditions are needed for Eq. (7): (a) the far-field boundary condition is 𝜔 = 0; = 0; (c) the wall boundary (b) the out-flow boundary condition is 𝜕𝜔 𝜕𝑛 condition is 𝜕𝜔 = 0. The homogeneous Neumann boundary condition 𝜕𝑛 is recommended at all the boundaries except for the far field, which is exhaustively reasoned by Cottet et al. (2000). The mass flux 𝑚 can be calculated by knowing the stream functions 𝜓 at both ends of the cell edge as given by Eq. (8). Notice that the direction of the edge is determined by the right-hand rule. The most substantial benefit of using Eq. (8) is that the net mass flux out of a control cell is strictly 0 since the evaluation is exact.

(5)

where, 𝑼 𝑏 is the velocity of the wall surface. Since a stationary 2-D cylinder is studied in this work, 𝑼 𝑏 = 𝟎. Therefore, the frame motion only enters into the solution of the problem through the initial and boundary conditions. In 2-D, the velocity 𝒖 can be represented by the gradient of stream function 𝜓 as 𝒖=(

𝜕𝜓 𝜕𝜓 ,− ). 𝜕𝑦 𝜕𝑥

(6)

The solving of the above system can be divided into a kinetic procedure and a kinematic procedure (Wu and Thompson, 1973; White, 1976). The kinetic process convects and diffuses the vorticity, produced from the wall surfaces, in the flow domain by the vorticity equation (4). The kinematic process, described by Eqs. (2) and (3), calculates the velocity (or stream function) corresponding to a given vorticity distribution. Also, one needs to account for appropriate boundary conditions for both steps. Actually, a closed-form integral expression that recovers the velocity field from a known vorticity field can be found (Thompson, 1971). The velocity field is sequentially used to determine the vorticity boundary conditions on the wall through a vorticity creation scheme. In the present method, the stream function is adopted instead of the velocity. The flowchart of general solving procedures is shown in Fig. 1.

𝑚 = 𝜓end − 𝜓start .

In the current study, the Alternating Directional Implicit (ADI) scheme proposed by Peaceman and Rachford (1955) is used to march the vorticity from the current time step 𝑡 to the next time step 𝑡 + 𝛥𝑡, where 𝛥𝑡 is the time step size. Eq. (7) is then split into two sweeps and expressed as [ ( )]−1 [ ( )] 𝐶1 + 𝐷1 1 + 𝛥𝑡 𝐶2 + 𝐷2 𝜔𝑡 , 𝜔𝑡+𝛥𝑡∕2 = 1 − 𝛥𝑡 2 2 (9) [ ( )]−1 [ ( )] 𝜔𝑡+𝛥𝑡 = 1 − 𝛥𝑡 𝐶2 + 𝐷2 1 + 𝛥𝑡 𝐶1 + 𝐷1 𝜔𝑡+𝛥𝑡∕2 , 2 2 where, 𝐶1,2 and 𝐷1,2 are the convective and diffusive coefficients in two directions. The interpolations of convective terms are evaluated by a convection limiter which will be described thoroughly in Section 2.3. The diffusive terms are calculated by a linear interpolation scheme since the grids used in this application are orthogonal.

2.2. Vorticity transport equation While solving the vorticity equation (4), the velocity 𝒖 is treated explicitly since the flow field is uniquely determined for a given vorticity distribution according to the Helmholtz velocity decomposition. Eq. (4) is solved by using the cell-centered FVM. Considering the configuration of a control cell surrounded by 4 edges from the structural mesh as shown in Fig. 2, and applying the Gauss theorem to the conservative form of Eq. (4) over the control cell, we have the discretized vorticity equation as 𝐴

4 4 ∑ 𝜕𝜔 𝜕𝜔 ∑ + 𝑚𝑓 𝜔𝑓 = 𝜈 𝛥𝑙 , 𝜕𝑡 𝑓 =1 𝜕𝑛𝑓 𝑓 𝑓 =1

(8)

2.3. Convection limiter The convection limiter is a hybrid algorithm combined the first order upwind scheme with the QUICK scheme. The QUICK scheme was originally introduced by Leonard (1979), and the formulation for nonuniform grids was derived by Ferziger and Peric (2012). Even though the QUICK scheme is third order in space and extremely conservative,

(7) 3

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Eq. (13) is solved by using Green’s function in 2-D (Thompson, 1971; Wu and Thompson, 1973) 𝜓(𝒙𝑐 ) = −

it encounters considerable difficulties in interpolating values involving sharp gradient changes. Especially in case of vorticity simulations, the vortical flow regions are enclosed by distinct borders, and the QUICK scheme would give nonphysical results near the borders where the impulsive gradient variations exist. Hence, a convection limiter is implemented in the present method. The first order upwind scheme is employed to keep the results monotone for large gradients. The QUICK scheme is used to maintain the order of accuracy for small and mild gradients. The QUICK scheme, corrected by the convection limiter 𝐹 (𝜃) (Woodfield et al., 2004), in terms of variable 𝜉 reads [ 𝐿(𝒙𝑒 , 𝒙𝑈 )𝐿(𝒙𝑒 , 𝒙𝑈 𝑈 ) 𝜉𝑒 = 𝜉𝑈 + 𝐹 (𝜃) (𝜉 − 𝜉𝑈 ) 𝐿(𝒙𝐷 , 𝒙𝑈 )𝐿(𝒙𝐷 , 𝒙𝑈 𝑈 ) 𝐷 ] 𝐿(𝒙𝑒 , 𝒙𝑈 )𝐿(𝒙𝑒 , 𝒙𝐷 ) + (𝜉𝑈 − 𝜉𝑈 𝑈 ) , (10) 𝐿(𝒙𝑈 , 𝒙𝑈 𝑈 )𝐿(𝒙𝐷 , 𝒙𝑈 𝑈 )

The convection and diffusion of vorticity are described by Eq. (4), and meanwhile the vorticity generated from the wall surface should also be included in the flow. The vorticity creation can be modeled as a surface vorticity distribution, or a vortex sheet, on the wall (Chorin, 1978). The strength of the vortex sheet is evaluated by enforcing the non-penetrating and no-slip boundary conditions of the velocity field determined by a known vorticity field (Koumoutsakos and Leonard, 1995). Assuming at time step 𝑡, the correct velocity field 𝒖𝑡 advances the vorticity from the current time step 𝜔𝑡 to the next time step 𝜔̂ 𝑡+𝛥𝑡 . However, 𝜔̂ 𝑡+𝛥𝑡 will induce a velocity field that has normal components 𝑡+𝛥𝑡 𝑡+𝛥𝑡 = 𝜕 𝜓̂𝜕𝑛 on the wall 𝑢̂ 𝑡+𝛥𝑡 = − 𝜕 𝜓̂𝜕𝑠 and tangential components 𝑢̂ 𝑡+𝛥𝑡 𝑠 𝑛 boundaries, where 𝑛 and 𝑠 are defined in Fig. 4. Note that the velocity 𝑢̂ 𝑡+𝛥𝑡 includes the contribution of both the vorticity and the background flow. In order to assure the non-penetrating boundary condition, a Laplace’s equation for the potential 𝜙 is solved (Cottet et al., 2000)

(11)

The convection limiter 𝐹 (𝜃) is ⎧0, ⎪𝜃∕𝛿, ⎪ 𝐹 (𝜃) = ⎨1, ⎪1 − 𝜃∕𝛿, ⎪ ⎩0,

𝜃 ≤ 0; 0 < 𝜃 ≤ 𝛿; 𝛿 < 𝜃 < 1 − 𝛿; 1 − 𝛿 ≤ 𝜃 < 1; ≥ 1;

(12)

∇2 𝜙 = 0,

(15)

with the kinematic boundary condition on the wall surface

where, 𝛿 represents the length of the hybrid region as shown in Fig. 3. 𝛿 should be wisely chosen to strike a good balance between less numerical dissipation and lower risks of failures. In the present study, 𝛿 is set to be 0.2, which is also recommended by Woodfield et al. (2004), from a series of numerical tests. In addition, the QUICK scheme will conflict with the ADI solver. Therefore, the high order terms in Eq. (10) are put into the RHS of the algebraic system, which is subsequently solved by using the deferred-correction technique (Richardson and Gaunt, 1927).

𝜕 𝜓̂ 𝑡+𝛥𝑡 𝜕𝜙 = −𝑢̂ 𝑡+𝛥𝑡 . (16) = 𝑛 𝜕𝑛 𝜕𝑠 In the present method, Eq. (15) is transformed into an integral form by using Green’s second identity in 2-D (Kellogg, 1953) and solved based on the Boundary Element Method (BEM) (Banerjee and Butterfield, 1981) to obtain the potential 𝜙 on the wall. The gradient of 𝜙 together with the tangential velocity 𝑢̂ 𝑡+𝛥𝑡 on the wall is converted 𝑠 into the surface vorticity distribution 𝛾 𝜕𝜙 𝜕 𝜓̂ 𝑡+𝛥𝑡 − . (17) 𝜕𝑠 𝜕𝑛 The surface vorticity 𝛾 shares the same dimension with the velocity (m/s). Actually, 𝛾 is defined as the vortex sheet strength 𝜔ℎ as the height of the sheet ℎ → 0. Notably, the potential 𝜙 only contributes to the surface vorticity 𝛾 on the wall and does not affect the other parts of the domain. The created surface vorticity 𝛾 is successively released into the flow field, which acts as implementing the no-slip boundary condition. Every point on the wall can release the created vorticity 𝜔𝛾 from the 𝛾=−

2.4. Stream function calculation The velocity field and the vorticity field are related through the stream function 𝜓. For a given vorticity field, the stream functions at the nodes of the control cells should be retrieved in order to march equation (7) to the next time step. Considering the relation between 𝒖 and 𝜓 in Eq. (6), Eq. (3) becomes ∇2 𝜓 = −𝜔.

(14)

2.5. Vorticity creation on the wall

where, 𝑒 , 𝑈 𝑈 , 𝑈 , and 𝐷 represent the current edge, the far upwind cell, the near upwind cell, and the downwind cell, respectively. The function 𝐿(𝒙1 , 𝒙2 ) calculates the distance between two locations 𝒙1 and 𝒙2 . 𝜃 is a non-dimensionalized variable that evaluates the gradient changes 𝜉𝑈 − min(𝜉𝑈 𝑈 , 𝜉𝐷 ) . |𝜉𝑈 𝑈 − 𝜉𝐷 |

𝜔(𝒙) ln |𝒙𝑐 − 𝒙| d𝜐 + 𝜓𝑖𝑛 , 2𝜋

here, 𝒙𝑐 is the location where 𝜓 is evaluated, and the dummy variable 𝒙 runs through the integral domain 𝛶 , covering all the vortical flow regions. 𝜓𝑖𝑛 = 𝑦𝑈1 − 𝑥𝑈2 is the stream function of the background flow 𝑼 = (𝑈1 , 𝑈2 ). Supposing that the vorticity is constant within every control cell, Eq. (14) is split into the numerical summation of the influence of all cells where the vorticity exists. An analytical solution of Eq. (14) for a control cell with uniform vorticity distribution is derived in Appendix. Applying Eq. (14) throughout the domain is relatively expensive. The computational complexity is O(𝑁 2 ) for 𝑁 cells in the vortical regions. However, the computation effort is nevertheless under control in 2-D applications, especially with proper parallelization as explained in Section 5. To achieve better performance on a single-core processor, several fast algorithms can be utilized such as the Fast Multipole Method (FMM) of O(𝑁 log(𝑁)) by Lindsay and Krasny (2001) and the FMM of O(𝑁) by Greengard and Rokhlin (1987). Alternatively, certain numerical methods, including the FVM, the FDM, or the Finite Element Method (FEM), for elliptic partial differential equations can be adopted to solve the Poisson equation (13), but one has to treat the non-orthogonality of the grids in this case.

Fig. 3. The convection limiter 𝐹 (𝜃) plotted against the dimensionless parameter 𝜃 (Woodfield et al., 2004).

𝜃=

∫𝛶

(13) 4

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Fig. 4. Coordinate system of the flow past a circular cylinder.

body surface. Koumoutsakos and Leonard (1995) established a relation 𝜕𝜔 between the surface vorticity 𝛾 and the normal vorticity flux 𝜕𝑛𝛾 across the wall 𝜕𝜔𝛾 𝛾 𝜈 =− . (18) 𝜕𝑛 𝛥𝑡 The surface vorticity 𝛾 enters the flow domain by a diffusion equation (19) with the Neumann type boundary condition (18) on the wall, the homogeneous far-field Dirichlet boundary condition 𝜔𝑡+𝛥𝑡 = 0, and 𝛾 the homogeneous initial condition 𝜔𝑡𝛾 = 0. 𝜕𝜔𝛾 𝜕𝑡

Fig. 5. Computational domain and mesh of the flow past a cylinder.

= 𝜈∇2 𝜔𝛾 .

(19)

Eq. (19) is numerically solved by using the FVM in the exact same way of solving equation (4). Eventually, the vorticity at time step 𝑡 + 𝛥𝑡 reads 𝜔𝑡+𝛥𝑡 = 𝜔̂ 𝑡+𝛥𝑡 + 𝜔𝑡+𝛥𝑡 . 𝛾

3. Numerical results of impulsively started flow past a circular cylinder

(20)

The proposed method is applied to the impulsively started flow, initialized as the inviscid flow with a slip boundary condition on the wall, past a circular cylinder to validate its accuracy and robustness in the inertial frame of reference. This type of problem is unanimously believed to be the standard benchmark test for unsteady separated flow and has been widely investigated. In this study, the cylinder with a diameter of 𝐷 locates at the origin (𝑥∕𝐷, 𝑦∕𝐷) = (0, 0) of the global coordinate system, and the uniform inflow comes in the 𝑥-direction. A circular computational domain with a radius of 14.5𝐷 is used for the simulations as shown in Fig. 5. The grid grows uniformly in the circumferential direction 𝒔 and exponentially in the normal direction 𝒏. Although such a large computational domain is overkill for the vorticity methods as indicated by Anderson and Reider (1996), remarkably, the calculations are only conducted in a small fraction of the grids, where the vorticity is non-zero. The Reynolds number 𝑅𝑒 is defined as

2.6. Surface pressure and sheer stress calculation In the viscous flow, the non-penetrating and no-slip boundary conditions require zero velocity on the fixed body surface. So, the unsteady term and the convective term in the momentum equation disappear on the body surface. Considering a local coordinate system on the surface of body (𝒔, 𝒏) as shown in Fig. 4, in which 𝒔 represents the unit vector tangential to the wall, and 𝒏 is the unit vector normal to the wall, the momentum equation is further simplified to 𝜕𝑝 𝜕𝜔 = −𝜇 , (21) 𝜕𝑠 𝜕𝑛 where, 𝑝 represents the pressure, and 𝜇 is the dynamic viscosity of the fluids. is already known from the vorticity creation process, In Eq. (21), 𝜕𝜔 𝜕𝑛 and 𝜕𝑝 describes how the pressure field changes on the body surface. 𝜕𝑠 Therefore, once the pressure at an arbitrary point is given, the pressure along the body surface can be evaluated according to Eq. (21). In the current work, the pressure at the upstream stagnation point is chosen to be the reference and calculated using Bernoulli’s Equation, which is 𝜌𝑈 2 . 2 The shear stress 𝜏 on the body surface is determined by the changing rate of tangential velocity 𝑢𝑠 in the normal direction. Recalling the vorticity defined by Eq. (3) in 2-D, since the normal velocity 𝑢𝑛 is strictly zero on the body surface, it can be shown that the surface shear stress is directly related to the vorticity on the wall

2𝑈 𝑅 , (24) 𝜈 where, 𝑅 = 𝐷∕2 is the radius of the cylinder, 𝑈 is the inflow velocity, and 𝜈 is the kinematic viscosity of the fluids. The Strouhal number 𝑆𝑡 is defined as 2𝑓 𝑅 2𝑅 𝑆𝑡 = i = , (25) 𝑈 𝑈 𝑇i 𝑅𝑒 =

where, 𝑓i and 𝑇i represent the frequency and period of the vortex shedding, respectively. The non-dimensional time 𝑇 is defined as

(22)

𝜏 = −𝜇𝜔.

𝑇 =

The total external force exerting on the body 𝑭 = (𝐹𝐷 , 𝐹𝐿 ) is the summation of pressure force 𝑭 𝑝 and friction force 𝑭 𝜏 , which are calculated by integrating the pressure and the shear stress along the body surface 𝑆𝑏 , respectively 𝑭 = 𝑭𝑝 + 𝑭𝜏 = −

∫𝑆𝑏

𝑝𝒏d𝑙 +

∫𝑆𝑏

𝜏𝒔d𝑙,

𝑡𝑈 . 𝑅

(26)

In this section, the flow development at early flow times during which the wake exhibits rather complex shedding structures for 𝑅𝑒 = 550, 1000, 3000, and 9500 along with the late stage vortex shedding phase for 𝑅𝑒 = 100 and 200 are presented. The descriptions of wake structures in the very early phase are provided for only half of the wake because the flow is symmetric with respect to the flow axis, which is also the 𝑥-axis defined in Fig. 4, before the periodic vortex shedding begins for these Reynolds numbers (Bouard and Coutanceau, 1980).

(23)

where, 𝐹𝐿 is the lift, and 𝐹𝐷 is the drag. They are the components of 𝑭 in the global Cartesian coordinate system. The lift coefficient 𝐶𝐿 and drag coefficient 𝐶𝐷 are non-dimensionalized by 𝜌𝑈 2 𝑅. 5

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Fig. 6. (a) grid dependence study and (b) time step size dependence study on the surface pressure coefficient distributions, 𝐶𝑝 , of an impulsively started cylinder plotted against the angle/𝜋 (rad) for 𝑅𝑒 = 3000 at flow time 𝑇 = 2.0. The present results are compared with those computed by Anderson and Reider (1996).

Fig. 7. Vorticity contours (—— represent contours with positive values, and – – – – indicate contours with negative values) computed by the present method (a) and streamlines comparison between the results of current method (b) and the experimental results of Bouard and Coutanceau (1980) (c) at 𝑅𝑒 = 550 and 𝑇 = 5.

6

Ocean Engineering 191 (2019) 106504

C. Wu et al.

3.1. Grid and time step size dependence study A grid dependence study is conducted for the impulsively started cylinder at Reynolds number 𝑅𝑒 = 3000 with various grid resolutions (circumferential direction × normal direction). Three levels of grid resolution are considered, including the coarse grid (100 × 80), the medium grid (200 × 135), and the fine grid (300 × 190). The surface pressure coefficient distributions at 𝑇 = 2.0, given in Fig. 6(a), verify that the solutions become independent on grids as the resolution increases. Also, the results are found to be in line with those calculated by Anderson and Reider (1996). Additionally, the influence of time step size on the results is investigated. The simulations are performed in the medium grid with a time step size 𝛥𝑇 = 1.0 × 10−2 , 2.0 × 10−3 , and 4.0×10−4 , respectively. It is obvious that the surface pressure coefficient distributions are not dependent on the time step size if 𝛥𝑇 is larger than 2.0 × 10−3 as shown in Fig. 6(b). During the simulations, the absolute total circulation is limited within 1.0 × 10−9 , and the absolute total net mass flux is controlled under 1.0 × 10−14 . 3.2. Impulsively started process for 𝑅𝑒 = 550, 1000, 3000, and 9500 Fig. 8. Radial velocity distributions along the rear centerline of the cylinder: (——), current; (◦), Subramaniam (1997); ( ), Qian and Vezza (2001); (▵), Dupuis et al. (2008) for 𝑅𝑒 = 550 at flow times: 𝑇 = 1.0, 𝑇 = 2.0, 𝑇 = 3.0, 𝑇 = 4.0, 𝑇 = 5.0.

In this part, the initial stages of the impulsively started flow past a cylinder are studied. The results from the present method are compared with other numerical results (Koumoutsakos and Leonard, 1995; Anderson and Reider, 1996; Subramaniam, 1997; Qian and Vezza, 2001) and experimental data (Bouard and Coutanceau, 1980). Koumoutsakos and Leonard (1995) and Subramaniam (1997) studied the early-stage growth of wake behind the cylinder using the mesh-free discrete vortex methods. Similar studies were carried out by Anderson and Reider (1996), who solved the N–S equations using the finite difference schemes of both second order and fourth order. Qian and Vezza (2001) simulated not only the start but also the longtime evolution of the flow past a cylinder based on the vorticity–velocity formulation of the N–S equations. Bouard and Coutanceau (1980) meticulously designed the experiments to achieve the abrupt initiation of the impulsively started cylinder and investigated the wake and separations behind the cylinder by visualizations at Reynolds numbers ranging from 40 to 104 . 3.2.1. 𝑅𝑒 = 550 Firstly, the impulsively started flow at 𝑅𝑒 = 550 is simulated using the proposed method. The grid with a number of elements (200 × 135) and the first layer grid size in the normal direction ℎ0 ∕𝐷 = 0.001 is utilized for the calculation. The time step size to march the flow in time is 𝛥𝑇 = 0.002. Fig. 7 compares the instantaneous streamlines produced by the present method with the visualized flow patterns obtained from the experiments of Bouard and Coutanceau (1980) at 𝑇 = 5.0. After the flow starts from rest, vorticity is created from the body surface and shed into the domain, and a recirculating region is formed by multiple vortical structures in the wake. Particularly in this case, a secondary vortex, bounded by a primary vortex, is observed in Fig. 7. The radial velocity distributions along the rear centerline of the cylinder are compared with those computed by Subramaniam (1997), Qian and Vezza (2001), and Dupuis et al. (2008) at flow times: 𝑇 = 1.0, 2.0, 3.0, 4.0, and 5.0, which are illustrated graphically in Fig. 8. The velocity is negative within the recirculating zone due to the reversed flow and reaches a minimum value at a specific location. The intercepts where the velocity profiles cross the 𝑥-axis (𝑢∕𝑈 = 0) other than the origin also reveal the evolution in time of the wake length. The wake length as well as the minimum velocity and its location can all be used to characterize the wake configuration, and the comparisons show good agreement. In addition, the development in time of drag coefficients calculated by the present method closely matches the predictions of both Koumoutsakos and Leonard (1995) and Qian and Vezza (2001) as given in Fig. 9.

Fig. 9. Comparisons of the drag coefficients evolution in time at 𝑅𝑒 = 550: (——), current; (◦), Koumoutsakos and Leonard (1995); ( ), Qian and Vezza (2001).

3.2.2. 𝑅𝑒 = 1000 Also, the flow at 𝑅𝑒 = 1000 is studied. In this case, the grid with a resolution (200 × 135) and ℎ0 ∕𝐷 = 0.001 is generated for the simulation. The time step size is 𝛥𝑇 = 0.002. The chronological development of the vorticity in the wake is given in Fig. 10 as equal value contours, demonstrating that the present method well captures the phenomena 𝛼 (Bouard and Coutanceau, 1980). It is seen, in particular, that a primary vortex is generated behind the cylinder in half of the recirculating region immediately after the flow starts and grows quickly with time. Later, a secondary vortex begins to develop at 𝑇 = 2.0 due to the separation of the primary vortex from a thin layer into a strong vortex, that rotates in the opposite direction of and is enclosed by the primary vortex. The vorticity contours from the current method are in good agreement with those provided by Qian and Vezza (2001), Anderson and Reider (1996), and Koumoutsakos and Leonard (1995). 7

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Fig. 10. Vorticity contours (—— represent contours with positive values, and – – – – indicate contours with negative values) of an impulsively started cylinder for 𝑅𝑒 = 1000 at flow times: (a), 𝑇 = 1.0; (b), 𝑇 = 2.0; (c), 𝑇 = 3.0; (d), 𝑇 = 4.0; (e), 𝑇 = 5.0; (f), 𝑇 = 6.0.

coefficients are almost the same as the results offered by Koumoutsakos and Leonard (1995) and Qian and Vezza (2001) as shown in Fig. 12. Actually, the variation in drag history is a result of the interactions between the primary and secondary vortices. In the beginning, the drag continuously drops due to the growth of primary vortex. After the primary vortex starts to separate from the cylinder, the formation and development of the secondary vortex raise the drag to a maximum value. Then, further separations of the primary vortex and the decay of the secondary vortex gradually reduce the drag.

The evolution of vortices behind the cylinder is directly manifested in the variation of surface vorticity. Fig. 11 shows the surface vorticity distributions compared with those from Qian and Vezza (2001) at flow times: 𝑇 = 1.0, 2.0, 3.0, 4.0, 5.0, and 6.0, which also show good agreement. The changes in the strength of the primary and secondary vortices in time are clearly observed in Fig. 11. The readers may notice from the curves that a tertiary vortex, hardly identified from the vorticity contours, emerges at later times. This vortex is more distinctly observable for higher Reynolds numbers. Additionally, the present drag 8

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Fig. 11. Comparisons of the surface vorticity distributions, 𝜔, of an impulsively started cylinder plotted against the angle/𝜋 (rad) with the results of Qian and Vezza (2001) for 𝑅𝑒 = 1000 at early flow times.

The streamlines at 𝑇 = 5 are compared with the experimental observations by Bouard and Coutanceau (1980) in Fig. 14, which also offers essential insights into the interactions between the primary vortex, the secondary vortex, and the tertiary vortex. In fact, the vortices in the wake significantly affect the surface vorticity. In Fig. 15, the surface vorticity distributions are presented along with the results from Qian and Vezza (2001) at flow times: 𝑇 = 1.0, 2.0, 3.0, 4.0, 5.0, and 6.0. Strong spatial fluctuations in the surface vorticity, which also exposes the changes of the vortex strength, are observed. The secondary vortex appears at about 𝑇 = 2.0, and its strength keeps increasing. As the secondary vortex grows and attempts to penetrate the primary vortex, the strength of the primary one decreases due to the diffusion and weakened feeding. Later, at 𝑇 = 3.0 another vortex is generated that fights against the secondary vortex from cutting off the feeding of the primary vortex. The interactions between the secondary and tertiary vortices increase the strength of both of them for a short period of time, and the primary vortex gains strength again due to the reinforced feeding. Eventually, the strength of all three vortices decays slowly through diffusion as the wake tends to be steady before the vortex shedding is initiated. Fig. 16 shows the comparisons of radial velocity distributions along the rear flow axis behind the cylinder with those computed by Subramaniam (1997) and Qian and Vezza (2001) at flow times: 𝑇 = 1.0, 2.0, 3.0, 4.0, and 5.0. In general, the results are almost the same.

Fig. 12. Comparisons of the drag coefficients evolution in time at 𝑅𝑒 = 1000: (——), current; (◦), Koumoutsakos and Leonard (1995); ( ), Qian and Vezza (2001).

3.2.3. 𝑅𝑒 = 3000 Further, the Reynolds number is increased to be 𝑅𝑒 = 3000. In this test, the grid with (200 × 135) cells and ℎ0 ∕𝐷 = 0.001 is used together with a time step size 𝛥𝑇 = 0.002 to advance in time. The wake evolution with time is described in equal value contours of vorticity in Fig. 13. At a higher Reynolds number, the primary vortex generates a larger adverse pressure gradient, which induces a stronger secondary vortex that subsequently creates a distinguished tertiary vortex. The tertiary vortex stops the secondary vortex from going through the primary one. These three vortical structures appear in a temporal and logical order. The primary vortex is firstly formed and dominates the flow, and the secondary vortex is induced as the primary vortex begins to isolate from the cylinder. Similarly, the development and separation of the secondary vortex create the tertiary vortex. This remarkable phenomenon has also been detected by Koumoutsakos and Leonard (1995), Anderson and Reider (1996), and Qian and Vezza (2001), which can be reliable sources for the examinations of the current results.

3.2.4. 𝑅𝑒 = 9500 The impulsively started flow at 𝑅𝑒 = 9500 is critical for the numerical methods because of the complex wake structures and separation details. For this simulation, the grid with a number of (300 × 200) cells and ℎ0 ∕𝐷 = 0.0005 is used in conjunction with a time step size 𝛥𝑇 = 0.002. Fig. 17 displays how the wake advances in time in terms of equi-vorticity contours and streamlines and also gives the experimental results from Bouard and Coutanceau (1980) as validations of the present work. The phenomenon 𝛽 is evidently seen during the development of flow (Bouard and Coutanceau, 1980). In the beginning, a thin layer of vortex is formed and soon grows to be the primary vortex. After a small lapse of time, this strong primary vortex splits the wake into 2 regions: a pair of secondary counter-rotating vortices appears near the front separation point, and the tertiary vortex separates from the primary one in the vicinity of the rear reattachment point. The contours from a high-resolution simulation are accessible in the publication of Koumoutsakos and Leonard (1995). 9

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Fig. 13. Vorticity contours (—— represent contours with positive values, and – – – – indicate contours with negative values) of an impulsively started cylinder for 𝑅𝑒 = 3000 at flow times: (a), 𝑇 = 1.0; (b), 𝑇 = 2.0; (c), 𝑇 = 3.0; (d), 𝑇 = 4.0; (e), 𝑇 = 5.0; (f), 𝑇 = 6.0.

The comparisons with lower Reynolds numbers confirm that the decrease in viscosity results in the growth of small scale structures in the wake as the Reynolds number increases. The same trend can also be observed from the rising oscillations of surface vorticity. The surface vorticity variations with time as computed by the current method, Koumoutsakos and Leonard (1995), and Qian and Vezza (2001) are shown in Fig. 18 at flow times: 𝑇 = 0.4, 1.0, 1.2, 1.8, and 2.0. The decrease in viscosity at higher Reynolds numbers leads to weaker diffusion. As a result, the wake tends to be unstable and vulnerable

to perturbations. Compared with the previous cases, the strength of the inner secondary vortex is much larger and increases quickly as the vortex grows under thin diffusion. This strong secondary vortex penetrates and splits the primary vortex rapidly since the viscosity is no longer effective enough to hold the vortices together. The primary vortex then starts to absorb the tertiary vortex and detach from the cylinder. Additionally, the radial velocity distributions along the rear centerline of the cylinder are compared with the results obtained by Subramaniam (1997) and Qian and Vezza (2001) at flow times: 𝑇 = 10

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Fig. 14. Streamlines comparison between (a) VISVE and (b) Bouard and Coutanceau (1980) at 𝑅𝑒 = 3000 and 𝑇 = 5.

Fig. 15. Comparisons of the surface vorticity distributions, 𝜔, of an impulsively started cylinder plotted against the angle/𝜋 (rad) with the results of Qian and Vezza (2001) for 𝑅𝑒 = 3000 at early flow times.

0.5, 1.0, 1.5, 2.0, and 2.5. The corresponding graphic representations are given in Fig. 19.

scheme. Therefore, the vortex shedding is triggered manually by a perturbation at a relatively early time. The impulsively started flow past a cylinder at 𝑅𝑒 = 100 and 𝑅𝑒 = 200 is simulated using the present method, respectively. The evolution of lift coefficient 𝐶𝐿 and drag coefficient 𝐶𝐷 eventually develops into steady sinusoidal fluctuations due to the periodic vortex shedding during the longtime simulation as plotted in Fig. 20. The lift coefficient, drag coefficient, and Strouhal number, predicted by the current method during the established shedding phase, are compared with the results of other latest studies as shown in Table 1. The present results agree well with them quantitatively. The amplitude of lift & drag coefficient and the Strouhal number at 𝑅𝑒 = 200 are larger than those at 𝑅𝑒 = 100 during the vortex shedding phase, which is caused by the drop of viscosity. Fig. 21 displays a snapshot of the alternating vortices

3.3. Vortex shedding in longtime simulation for 𝑅𝑒 = 100 and 200 In this part, the longtime evolution of impulsively started flow past a circular cylinder at 𝑅𝑒 = 100 and 200 is presented to assess the present method in capturing the periodic vortex shedding. The mesh with (200 × 210) elements and the first layer grid size in the normal direction ℎ0 ∕𝐷 = 0.001 is adopted with a time step size 𝛥𝑇 = 0.002. The radius of the computational domain is increased to be 17.0𝐷. The vortex shedding phase will not be automatically activated by the numerical errors even when 𝑇 approaches 200 because of the highly conservative 11

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Table 1 Comparison of lift coefficient 𝐶𝐿 , drag coefficient 𝐶𝐷 , and Strouhal number 𝑆𝑡 at 𝑅𝑒 = 100 and 𝑅𝑒 = 200 during the shedding phase. Reynolds number

Author(s)

𝐶𝐷

𝐶𝐿

𝑆𝑡

100

Liu et al. (1998) Calhoun (2002) Russell and Wang (2003) Ding et al. (2004) Choi et al. (2007) Chiu et al. (2010) Hejranfar and Ezzatneshan (2014) Cui et al. (2017) Present

1.35 ± 0.012 1.330 ± 0.014 1.38 ± 0.007 1.325 ± 0.008 1.34 ± 0.011 1.35 ± 0.012 1.3360 ± 0.0114 1.358 ± 0.013 1.314 ± 0.009

±0.339 ±0.298 ±0.300 ±0.28 ±0.315 ±0.303 ±0.3354 ±0.342 ±0.308

0.165 0.175 0.169 0.164 0.164 0.167 0.1626 0.167 0.163

200

Liu et al. (1998) Calhoun (2002) Wright and Smith (2001) Russell and Wang (2003) Ding et al. (2004) Choi et al. (2007) Chiu et al. (2010) Hejranfar and Ezzatneshan (2014) Cui et al. (2017) Present

1.31 ± 0.049 1.172 ± 0.058 1.33 ± 0.04 1.29 ± 0.022 1.327 ± 0.045 1.36 ± 0.048 1.37 ± 0.051 1.329 ± 0.0439 1.340 ± 0.049 1.321 ± 0.045

±0.69 ±0.668 ±0.68 ±0.50 ±0.60 ±0.64 ±0.71 ±0.6814 ±0.655 ±0.669

0.192 0.202 0.196 0.195 0.196 0.191 0.198 0.1941 0.196 0.195

where, 𝑈o is the maximum velocity of the oscillatory inflow, and 𝑓o is the frequency of the velocity oscillation. The oscillatory flow is characterized by the Keulegan–Carpenter number 𝐾𝐶 and the Reynolds number 𝑅𝑒. The Reynolds number of the oscillatory flow is defined as 𝑈o 𝐷 , (28) 𝜈 where, 𝜈 is the kinematic viscosity of the fluids. The Keulegan– Carpenter number is defined as

𝑅𝑒 =

𝑈o 𝑇 o , (29) 𝐷 where, 𝑇o = 1∕𝑓o is the period of the oscillatory flow. Further, the ratio of 𝑅𝑒 to 𝐾𝐶 gives the frequency parameter or Stokes number 𝛽, which is expressed as 𝐾𝐶 =

𝛽=

(30)

2-D oscillatory flow past a circular cylinder only exists at low 𝐾𝐶 and 𝑅𝑒 numbers (Tatsuno and Bearman, 1990). Within the scope of current work, the oscillatory flow at 𝑅𝑒 = 100 & 𝐾𝐶 = 5 is firstly studied in this section to validate the proposed method in the non-inertial translating frame of reference. The present results are checked against other numerical schemes and experimental visualizations. Dütsch et al. (1998) numerically reproduced three regimes found by Tatsuno and Bearman (1990) in 2-D, and in particular provided experimental measurements as validations. Guilmineau and Queutey (2002) carried out similar studies on a 2-D circular cylinder subject to inline oscillations with extensions to cross-flow oscillations. An et al. (2011) conducted direct numerical simulations to investigate the Honji vortices (Honji, 1981) in oscillatory flow past a 3-D circular cylinder, and discovered that the instabilities are directly related to 𝛽. Upon validations, the results of oscillatory flow at 𝑅𝑒 = 200 & 𝐾𝐶 = 1, 2, and 4 are further presented to investigate the influence of 𝐾𝐶 number on the inline forces exerting on the cylinder. The 𝐾𝐶 and 𝑅𝑒 numbers used here are small such that the flow is symmetric and subcritical, and laminar boundary layers can be established (Sarpkaya, 1986).

Fig. 16. Radial velocity distributions along the rear centerline of the cylinder: (——), current; (◦), Subramaniam (1997); ( ), Qian and Vezza (2001) for 𝑅𝑒 = 3000 at flow times: 𝑇 = 1.0, 𝑇 = 2.0, 𝑇 = 3.0, 𝑇 = 4.0, 𝑇 = 5.0.

shed behind the cylinder in terms of vorticity contours. A much shorter distance between two neighboring vortices is observed at 𝑅𝑒 = 200 than that at 𝑅𝑒 = 100. Also, the vortices are dissipated much faster by the larger viscosity at 𝑅𝑒 = 100. 4. Numerical results of oscillatory flow past a circular cylinder The inline oscillatory flow past a stationary circular cylinder is simulated to show the capability and accuracy of the proposed method under the non-inertial frame in translation. This problem is of diverse interests to both hydrodynamic and structural academia (Blevins, 1977). In addition, the accurate predictions of time-dependent loads on structures matter enormously in the civil and offshore industry (Naudascher and Rockwell, 2012). Similar to the unidirectional flow, the cylinder with a diameter of 𝐷 is fixed at the origin (𝑥∕𝐷, 𝑦∕𝐷) = (0, 0) of the global coordinate system. The circular computational domain has a radius of 10.0𝐷, and the flow oscillates along the 𝑥-direction with a sinusoidal velocity 𝑈 = 𝑈o cos(2𝜋𝑓o 𝑡),

𝑅𝑒 𝐷2 = . 𝐾𝐶 𝜈𝑇o

4.1. 𝑅𝑒 = 100 And 𝐾𝐶 = 5 In this case, the grid with a resolution (300 × 150) and ℎ0 ∕𝐷 = 0.001 is built for the simulation. The time step size is 𝛥𝑇 = 0.002. According to the classifications of Tatsuno and Bearman (1990), a combination of 𝐾𝐶 = 5 and 𝑅𝑒 = 100 falls into the category of regime A, of which the flow is 2-D and symmetric with respect to the 𝑥-axis, and separations occur periodically. Fig. 22 displays the development of vorticity in half a cycle when the cylinder begins to move in the positive

(27) 12

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Fig. 17. Vorticity contours (left column), streamlines (middle column), and experimental results from Bouard and Coutanceau (1980) (right column) of an impulsively started cylinder for 𝑅𝑒 = 9500 at early flow times. 1 from thin layers into a pair of dominant causes the growth of vortices ⃝ 2 which gradually vortices, accompanied by the separation of vortices ⃝, 3 due to diffusion. Eventually, the feeding fade together with vortices ⃝

𝑥-direction from its maximum negative location. At first, two thin layers 1 partially bounded by a pair of vortices ⃝ 2 and a pair of vorticity ⃝ 3 are observed. The motion of the cylinder of fully detached vortices ⃝ 13

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Fig. 18. Comparisons of the surface vorticity distributions, 𝜔, of an impulsively started cylinder plotted against the angle/𝜋 (rad) with the results of Qian and Vezza (2001) and Koumoutsakos and Leonard (1995) for 𝑅𝑒 = 9500 at early flow times.

present method, which provides a different perspective to view the same problem under the non-inertial frame. It is noteworthy to point out that similar deviations in velocity profiles also exist between the experimental data of Dütsch et al. (1998) and the numerical results of the mentioned three studies (Dütsch et al., 1998; Guilmineau and Queutey, 2002; An et al., 2011). Finally, the variations of inline force coefficients, non-dimensionalized by 𝜌𝑈o2 𝑅, in time are reported in Fig. 25. 4.2. 𝑅𝑒 = 200 And 𝐾𝐶 = 1, 2, 4 In this part, the grid and time step size are the same as the previous case. The Reynolds number is raised to be 200. The flow at a 𝐾𝐶 number of 1, 2, and 4 is simulated using the current method, respectively. The vorticity contours are shown in Fig. 26 at a phase angle of 180◦ for 𝐾𝐶 = 1, 2, and 4. It is seen that the vorticity is attached to the cylinder at lower 𝐾𝐶 numbers, but separations will occur for the flow at larger 𝐾𝐶 numbers. The developments of inline force coefficients, non-dimensionalized by 𝜌𝑈o2 𝑅, in time are displayed in Fig. 27. Similar to the case at 𝑅𝑒 = 100 and 𝐾𝐶 = 5, the pressure part contributes much more to the total force than the friction part does for all 𝐾𝐶 numbers presented in this work. The authors also notice a drop in the total force, mostly resulting from the reduction of pressure force, when the 𝐾𝐶 number increases. The reason for this decline is that a higher 𝐾𝐶 number means a larger period of the oscillation and a lower acceleration, which generates a milder adverse pressure gradient on the down-stream side of the cylinder. A closer look at the friction force is given in Fig. 28, demonstrating the decrease of friction force as the 𝐾𝐶 number rises, which shares the same trend of how the pressure force changes with the 𝐾𝐶 number. According to Eq. (22), the friction force is proportional to the intensity of vorticity, also the velocity gradient, on the wall. The vorticity has enough time to convect away and detach from the cylinder for higher 𝐾𝐶 numbers while in contrast the vorticity is largely confined to the regions near the cylinder for lower 𝐾𝐶 numbers as shown in Fig. 26. Therefore, the flow at smaller 𝐾𝐶 numbers is characterized with larger vorticity magnitude in the vicinity of the cylinder, resulting in bigger friction force.

Fig. 19. Radial velocity distributions along the rear centerline of the cylinder: (——), current; (◦), Subramaniam (1997); ( ), Qian and Vezza (2001) for 𝑅𝑒 = 9500 at flow times: 𝑇 = 0.5, 𝑇 = 1.0, 𝑇 = 1.5, 𝑇 = 2.0, 𝑇 = 2.5.

1 are cut off by the formation of two new thin layers links of vortices ⃝ 4 on the wall surface. of vorticity ⃝ Fig. 23 shows the comparisons of velocity profiles at different locations and phase angles (2𝜋𝑓o 𝑡) between the current solutions and the experimental results of Dütsch et al. (1998). Notably, the measurements by Dütsch et al. (1998) were made with respect to the fluids. On the contrary, the current numerical simulations adopted a cylinder fixed global coordinate system. Therefore, coordinates transformations were conducted to produce reasonable comparisons. As is seen in Fig. 23, the results share the same trend in general but still exhibit some visual discrepancies. To further demonstrate the robustness of the proposed method, the velocity profiles are compared with the numerical results predicted by Dütsch et al. (1998), Guilmineau and Queutey (2002) as well as An et al. (2011) in Fig. 24, which shows good agreement. In the simulations of Dütsch et al. (1998) and Guilmineau and Queutey (2002), the cylinder is forced to oscillate in the quiescent flow. In contrast, the flow field alternates around the fixed cylinder in the

5. Computational efficiency and parallelization To show the efficiency of the proposed vorticity method, an analysis of the computational complexity is conducted in the simulation of an 14

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Fig. 20. Longtime evolution of drag (——) and lift (– – – –) coefficients of the impulsively started cylinder at: (a) 𝑅𝑒 = 100 and (b) 𝑅𝑒 = 200.

Fig. 21. Vorticity contours (—— represent contours with positive values, and – – – – indicate contours with negative values) of a longtime simulation at: (a) 𝑅𝑒 = 100 and (b) 𝑅𝑒 = 200.

to quantify as indicated by Kevlahan and Vasilyev (2005). Additionally, the size of the computational domain and the number of cells in the grid are compared between the current method and the N–S methods as shown in Table 2 at a Reynolds number around 200. As is seen, the present method has a much smaller computational domain and a less number of cells. The authors also remark that the actual number of occupied cells involved in the computations is only a small fraction of the total number. It, therefore, is obvious that the present method is more efficient compared with the other methods in terms of computational complexity. The coarse grid with much fewer cells produces comparably fine results of the other methods. In order to further improve the computational efficiency, the present numerical scheme is paralleled. The most time-consuming part, which turns out to be the stream function calculations and takes 99% of the total time cost, is located by profiling. Since the data of spatially

impulsively started cylinder at 𝑅𝑒 = 3000 and 𝑅𝑒 = 200. The cells actively involved in the calculations are defined as the occupied cells, whose vorticity magnitude is larger than a predefined value 10−6 in this study. Firstly, at 𝑅𝑒 = 3000, the number of occupied cells used in the present method is compared, in particular, with those of Koumoutsakos and Leonard (1995) and Kevlahan and Vasilyev (2005), who put forward a wavelet method to calculate the flow around a cylinder. At flow time 𝑇 = 6.0, to resolve the fine wake comprised of several eddies, 3.0 × 105 particles are employed by Koumoutsakos and Leonard (1995) in the discrete vortex method, as opposed to only 1.3 × 104 , which is about 23 times smaller, occupied cells utilized in the current method. In the work of Kevlahan and Vasilyev (2005), 4.3 × 104 , which is 3.3 times larger than the number of occupied cells in the present method, wavelets are necessary to capture the wake structures. However, the computational cost of one wavelet element over a vorticity cell is hard 15

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Fig. 22. Vorticity contours (—— represent contours with positive values, and – – – – indicate contours with negative values) of alternating flows around a cylinder for 𝐾𝐶 = 5 and 𝑅𝑒 = 100 (𝛽 = 20) at phase angles: (a) 90◦ ; (b) 126◦ ; (c) 162◦ ; (d) 198◦ ; (e) 234◦ ; (f) 270◦ .

Table 2 Comparison of computational domain size and number of cells in grid. Author(s)

Method

Domain size

Number of mesh cells

Calhoun (2002) Russell and Wang (2003) Linnick and Fasel (2005) Taira and Colonius (2007) Chiu et al. (2010) Present

N–S N–S N–S N–S N–S VISVE

Rectangular 32𝐷 × 16𝐷 Rectangular 32𝐷 × 16𝐷 Rectangular 46𝐷 × 16𝐷 Rectangular 30𝐷 × 30𝐷 Rectangular 40𝐷 × 20𝐷 Circular (Radius) 17𝐷

204,800 204,800 180,121 106,250 245,000 42,000

connected cells is not necessarily stored in the memory continuously, the vectorization, which is the vector calculation of a processor, does not function effectively in case of random memory access. However, the evaluations of Eq. (14) are perfect for parallel computing, in which independent calculations are performed by multiple computer processors simultaneously. Because the computations of inductions using Eq. (14) are not dependent on each other, the work can be distributed to available processors and conducted at the same time. The speed-up of parallelization in terms of strong scaling (Amdahl, 1967) and the number of processors used exhibit an almost ideally linear relationship

in Fig. 29, in which the speed-up is defined as 𝑡1 ∕𝑡𝑛 where 𝑡𝑛 is the time cost using 𝑛 processors, and 𝑡1 is the time cost using 1 processor. 6. Conclusions and future work A conservative vorticity method that solves the vorticity-stream function formulation of Navier–Stokes equations in the Eulerian coordinates is proposed to simulate the incompressible unsteady viscous flow past bluff bodies in 2-D. The numerical scheme is based on the modified integro-differential technique, which permits the simulation 16

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Fig. 23. Comparisons of the velocity profiles (left column 𝑢; right column 𝑣) at various 𝑥∕𝐷 locations for 𝐾𝐶 = 5 and 𝑅𝑒 = 100 (𝛽 = 20) between present results and experimental measurements by Dütsch et al. (1998) at three phase angles: (a) 180◦ ; (b) 210◦ ; (c) 330◦ .

17

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Fig. 24. Velocity components profiles (a) 𝑢 and (b) 𝑣 at phase angle 180◦ and location 𝑥∕𝐷 = −0.6 for 𝐾𝐶 = 5 and 𝑅𝑒 = 100 (𝛽 = 20). The present results are compared with the numerical results computed by Dütsch et al. (1998), Guilmineau and Queutey (2002) and An et al. (2011).

to be limited to only the regions with non-zero vorticity. The computation is then split into two consecutive steps, that are the solving of vorticity transport equation using the cell-centered FVM and the recovery of stream functions using an analytical solution. The exact evaluations of mass flux from the explicitly calculated stream functions further secure the continuity of velocity field. To ensure the accuracy and monotonicity of vorticity interpolations, a convection limiter coupling the first order upwind scheme and the QUICK scheme is implemented. In addition, a vorticity creation scheme, effectively resolving the impulsive initiation and boundary layer distribution, is employed to determine the amount of vorticity released from the wall by imposing the boundary conditions on the body surface. The conservative FVM and QUICK scheme and the exact evaluations of stream functions and convective flux together construct a numerical scheme that ensures the conservation of circulation and mass. The proposed vorticity method has been validated by simulating the impulsively started flow past a circular cylinder at Reynolds numbers ranging from 100 to 9500 in the inertial frame of reference and the oscillatory flow around a circular cylinder at 𝑅𝑒 = 100 & 𝐾𝐶 = 5 as well as 𝑅𝑒 = 200 & 𝐾𝐶 = 1, 2, 4 in the non-inertial translational frame of reference. The present results coincide well with those from other established numerical and experimental benchmark tests, demonstrating that the current method produces high-quality results for both earlystage evolution and longtime vortex shedding phase of the flow past a circular cylinder even in comparatively coarse grids. The current method can simulate 2-D laminar flow past solid bodies in both inertial and non-inertial translational frame of reference in an efficient way. Even though only the results of the flow past a circular cylinder are shown in this work, the method can be applied to lifting or bluff bodies of arbitrary shapes. The method can also handle turbulent flow past a 2-D hydrofoil (Yao and Kinnas, 2019) with appropriate modifications or cavitating flow (Xing et al., 2018). Also, this method is being extended to 3-D, and a preliminary study about the flow past a 3-D hydrofoil has been shown by Wu and Kinnas (2019). Due to the advantage of promising efficiency, this method has the potential for simulating the flow around 3-D long cylinders, such as the risers connecting a subsea oil well to an offshore platform.

Fig. 25. Chronological evolution of inline force coefficients for 𝐾𝐶 = 5 and 𝑅𝑒 = 100 (𝛽 = 20): (——), total force; (– – – –), pressure force; (– . – . – .), friction force.

on Cavitation Performance of High Speed Propulsors’’. The authors are grateful for the inspiring ideas provided by Dr. Ye Tian to complete this method and finish this work. The authors also acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. Appendix. Vorticity induced stream function The induction from a control cell 𝐴 with constant vorticity distribution 𝜔 to a free space point with a relative location vector 𝒓 = (𝑥, 𝑦) reads 𝜔 ln 𝑟d𝑠, 2𝜋 ∫𝐴 √ where, 𝑟 = 𝑥2 + 𝑦2 . Considering the following relationship 𝜓 =−

ln 𝑟 = Acknowledgments

𝛁⋅(𝑥 ln 𝑟2 − 1, 𝑦 ln 𝑟2 − 1) . 4

(A.1)

(A.2)

Eq. (A.1) is rewritten as a line integral form based on the Gauss theorem

Support for this research was provided partly by the U.S. Office of Naval Research (Grant Nos. N00014-14-1-0303 and N00014-18-12276; Dr. Ki-Han Kim) and by Phases VII and VIII of the ‘‘Consortium

𝜓 =− 18

𝜔 (𝑥 ln 𝑟2 − 1)d𝑦 − (𝑦 ln 𝑟2 − 1)d𝑥. 8𝜋 ∫𝜕𝐴

(A.3)

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Fig. 26. Vorticity contours (—— represent contours with positive values, and – – – – indicate contours with negative values) of alternating flows around a cylinder at 𝑅𝑒 = 200 with: (a) 𝐾𝐶 = 1; (b) 𝐾𝐶 = 2; (c) 𝐾𝐶 = 4 for phase angle 180◦ .

Fig. 28. Inline pressure force coefficients at 𝑅𝑒 = 200 with: (——), 𝐾𝐶 = 1; (– – – –), 𝑡𝑈 𝐾𝐶 = 2; (– . – . – .), 𝐾𝐶 = 4. The non-dimensionalized time 𝑇 = 𝑅𝑜 .

Fig. 29. The speed-up in terms of strong scaling for 𝑅𝑒 = 3000 in the medium grid.

Now, for an edge with a starting point 𝒓0 = (𝑥0 , 𝑦0 ) and an end point 𝒓1 = (𝑥1 , 𝑦1 ) in Fig. A.30, we have

Fig. 27. Chronological evolution of inline force coefficients at 𝑅𝑒 = 200 with (a) 𝐾𝐶 = 1, (b) 𝐾𝐶 = 2, and (c) 𝐾𝐶 = 4: (——), total force; (– – – –), pressure force; (– . 𝑡𝑈 – . – .), friction force. The non-dimensionalized time 𝑇 = 𝑅𝑜 .

𝑥 𝑦 19

= =

𝜅1 𝜒 + 𝑥0 , 𝜅2 𝜒 + 𝑦0 ,

(A.4)

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Choi, J.-I., Oberoi, R.C., Edwards, J.R., Rosati, J.A., 2007. An immersed boundary method for complex incompressible flows. J. Comput. Phys. 224 (2), 757–784. Chorin, A.J., 1973. Numerical study of slightly viscous flow. J. Fluid Mech. 57 (4), 785–796. Chorin, A.J., 1978. Vortex sheet approximation of boundary layers. J. Comput. Phys. 27 (3), 428–442. Chorin, A.J., Bernard, P.S., 1973. Discretization of a vortex sheet, with an example of roll-up. J. Comput. Phys. 13 (3), 423–429. Cottet, G.-H., Koumoutsakos, P.D., Petros, D., Koumoutsakos, P.D., et al., 2000. Vortex Methods: Theory and Practice. Cambridge university press. Coutanceau, M., Bouard, R., 1977a. Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform translation. Part 1. Steady flow. J. Fluid Mech. 79 (2), 231–256. Coutanceau, M., Bouard, R., 1977b. Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform translation. Part 2. Unsteady flow. J. Fluid Mech. 79 (2), 257–272. http://dx.doi.org/10.1017/ S0022112077000147. Cui, X., Yao, X., Wang, Z., Liu, M., 2017. A hybrid wavelet-based adaptive immersed boundary finite-difference lattice Boltzmann method for two-dimensional fluid–structure interaction. J. Comput. Phys. 333, 24–48. Dennis, S., Chang, G.-Z., 1970. Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100. J. Fluid Mech. 42 (3), 471–489. Ding, H., Shu, C., Yeo, K., Xu, D., 2004. Simulation of incompressible viscous flows past a circular cylinder by hybrid FD scheme and meshless least square-based finite difference method. Comput. Methods Appl. Mech. Engrg. 193 (9–11), 727–744. Dupuis, A., Chatelain, P., Koumoutsakos, P., 2008. An immersed boundary–latticeBoltzmann method for the simulation of the flow past an impulsively started cylinder. J. Comput. Phys. 227 (9), 4486–4498. Dütsch, H., Durst, F., Becker, S., Lienhart, H., 1998. Low-Reynolds-number flow around an oscillating circular cylinder at low Keulegan–Carpenter numbers. J. Fluid Mech. 360, 249–271. Elston, J.R., Blackburn, H.M., Sheridan, J., 2006. The primary and secondary instabilities of flow generated by an oscillating circular cylinder. J. Fluid Mech. 550, 359–389. Ferziger, J.H., Peric, M., 2012. Computational Methods for Fluid Dynamics. Springer Science & Business Media. Greengard, L., Rokhlin, V., 1987. A fast algorithm for particle simulations. J. Comput. Phys. 73 (2), 325–348. Guilmineau, E., Queutey, P., 2002. A numerical simulation of vortex shedding from an oscillating circular cylinder. J. Fluids Struct. 16 (6), 773–794. Guo, Z., Zhao, T.-S., 2003. Explicit finite-difference lattice Boltzmann method for curvilinear coordinates. Phys. Rev. E 67 (6), 066709. Hejranfar, K., Ezzatneshan, E., 2014. Implementation of a high-order compact finitedifference lattice Boltzmann method in generalized curvilinear coordinates. J. Comput. Phys. 267, 28–49. Honji, H., 1981. Streaked flow around an oscillating circular cylinder. J. Fluid Mech. 107, 509–520. Imamura, T., Suzuki, K., Nakamura, T., Yoshida, M., 2005. Acceleration of steady-state lattice Boltzmann simulations on non-uniform mesh using local time step method. J. Comput. Phys. 202 (2), 645–663. Justesen, P., 1991. A numerical study of oscillating flow around a circular cylinder. J. Fluid Mech. 222, 157–196. Kellogg, O.D., 1953. Foundations of Potential Theory, Vol. 31. Courier Corporation. Kevlahan, N.K., Vasilyev, O.V., 2005. An adaptive wavelet collocation method for fluid-structure interaction at high Reynolds numbers. SIAM J. Sci. Comput. 26 (6), 1894–1915. Koumoutsakos, P., Leonard, A., 1995. High-resolution simulations of the flow around an impulsively started cylinder using vortex methods. J. Fluid Mech. 296, 1–38. Krasny, R., 1987. Computation of vortex sheet roll-up in the Trefftz plane. J. Fluid Mech. 184, 123–155. Leonard, B.P., 1979. A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Comput. Methods Appl. Mech. Engrg. 19 (1), 59–98. Leonard, A., 1980. Vortex methods for flow simulation. J. Comput. Phys. 37 (3), 289–335. Li, Y., Shock, R., Zhang, R., Chen, H., 2004. Numerical study of flow past an impulsively started cylinder by the lattice-Boltzmann method. J. Fluid Mech. 519, 273–300. Lindsay, K., Krasny, R., 2001. A particle method and adaptive treecode for vortex sheet motion in three-dimensional flow. J. Comput. Phys. 172 (2), 879–907. Linnick, M.N., Fasel, H.F., 2005. A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains. J. Comput. Phys. 204 (1), 157–192. Liu, C., Zheng, X., Sung, C., 1998. Preconditioned multigrid methods for unsteady incompressible flows. J. Comput. Phys. 139 (1), 35–57. Marquet, O., Sipp, D., Jacquin, L., 2008. Sensitivity analysis and passive control of cylinder flow. J. Fluid Mech. 615, 221–252. Mei, R., Shyy, W., 1998. On the finite difference-based lattice Boltzmann method in curvilinear coordinates. J. Comput. Phys. 143 (2), 426–448. Moore, D., 1974. A numerical study of the roll-up of a finite vortex sheet. J. Fluid Mech. 63 (2), 225–235.

Fig. A.30. Evaluation of the induced stream function due to a polygonal cell with a constant vorticity distribution.

where, 𝜒 indicates the small increment along the edge, and 𝜅1 & 𝜅2 are constants defined as 𝜅1

=

√

𝜅2

=

√

𝑥1 −𝑥0

(𝑥1 −𝑥0 )2 +(𝑦1 −𝑦0 )2 𝑦1 −𝑦0 (𝑥1 −𝑥0 )2 +(𝑦1 −𝑦0 )2

, .

(A.5)

Therefore, applying Eq. (A.3) to the edge 𝒓

𝜔 ∫𝒓 1 (𝑥 ln 𝑟2 − 1)d𝑦 − (𝑦 ln 𝑟2 − 1)d𝑥 − 8𝜋 0

=

𝛥𝑙

𝜔 ∫0 [(𝜅2 𝑥0 − 𝜅1 𝑦0 ) ln 𝑟2 − (𝜅2 − 𝜅1 )]d𝜒 − 8𝜋

(A.6) { [ 𝜅2 𝑥0 −𝜅1 𝑦0 = (𝜅2 𝑥0 − 𝜅1 𝑦0 ) (2(𝜅1 𝑦0 − 𝜅2 𝑥0 ) arctan 𝜅 𝑥 +𝜅 𝑦 +𝜒 ) 1 0 2 0 } ]|𝛥𝑙 +(𝜅1 𝑥0 + 𝜅2 𝑦0 + 𝜒) ln 𝑟2 − 2𝜒 | + (𝜅1 − 𝜅2 )𝛥𝑙 , |0 √ 2 2 where, 𝛥𝑙 = (𝑥1 − 𝑥0 ) + (𝑦1 − 𝑦0 ) is the length of the edge. Note that Eq. (A.6) also works for general non-orthogonal grids with arbitrary polygonal cells. The summation of contributions from all the edges of a control cell gives the induced stream function. 𝜔 − 8𝜋

References Amdahl, G.M., 1967. Validity of the single processor approach to achieving large scale computing capabilities. In: Proceedings of the April 18-20, 1967, Spring Joint Computer Conference. ACM, pp. 483–485. An, H., Cheng, L., Zhao, M., 2011. Direct numerical simulation of oscillatory flow around a circular cylinder at low Keulegan–Carpenter number. J. Fluid Mech. 666, 77–103. Anderson, C.R., Reider, M.B., 1996. A high order explicit method for the computation of flow about a circular cylinder. J. Comput. Phys. 125 (1), 207–224. Badr, H., Dennis, S., Kocabiyik, S., Nguyen, P., 1995. Viscous oscillatory flow about a circular cylinder at small to moderate Strouhal number. J. Fluid Mech. 303, 215–232. Banerjee, P.K., Butterfield, R., 1981. Boundary Element Methods in Engineering Science, Vol. 17. McGraw-Hill London. Blevins, R.D., 1977. Flow-induced vibration. Van Nostrand Reinhold Co., New York, p. 377. Bouard, R., Coutanceau, M., 1980. The early stage of development of the wake behind an impulsively started cylinder for 40< Re< 10 4. J. Fluid Mech. 101 (3), 583–607. Calhoun, D., 2002. A Cartesian grid method for solving the two-dimensional streamfunction-vorticity equations in irregular regions. J. Comput. Phys. 176 (2), 231–275. Chakrabarti, S.K., 1987. Hydrodynamics of Offshore Structures. WIT press. Chang, C.-C., Chern, R.-L., 1991. A numerical study of flow around an impulsively started circular cylinder by a deterministic vortex method. J. Fluid Mech. 233, 243–263. Chiu, P.-H., Lin, R.-K., Sheu, T.W., 2010. A differentially interpolated direct forcing immersed boundary method for predicting incompressible Navier–Stokes equations in time-varying complex geometries. J. Comput. Phys. 229 (12), 4476–4500. 20

Ocean Engineering 191 (2019) 106504

C. Wu et al.

Subramaniam, S., A new mesh-free vortex method.. Taira, K., Colonius, T., 2007. The immersed boundary method: a projection approach. J. Comput. Phys. 225 (2), 2118–2137. Tatsuno, M., Bearman, P.W., 1990. A visual study of the flow around an oscillating circular cylinder at low Keulegan–Carpenter numbers and low Stokes numbers. J. Fluid Mech. 211, 157–182. http://dx.doi.org/10.1017/S0022112090001537. Thompson, J.F., 1971. Two Approaches to the Three-Dimensional Jet-in-Cross-Wind Problem: a Vortex Lattice Model and a Numerical Solution of the Navier-Stokes Equations (Ph.D thesis). Georgia Institute of Technology. Tian, Y., 2014. Leading edge vortex modeling and its effect on propulsor performance. (Ph.D thesis). University of Texas at Austin. Tian, Y., Kinnas, S.A., 2015. A viscous vorticity method for propeller tip flows and leading edge vortex. In: Fourth international symposium on marine propulsors (smp’15), Austin, Texas, USA. White, D.J.S., 1976. Numerical boundary conditions for viscous flow problems. AIAA J. 14 (8), 1042–1049. Williamson, C., 1985. Sinusoidal flow relative to circular cylinders. J. Fluid Mech. 155, 141–174. Woodfield, P.L., Suzuki, K., Nakabe, K., 2004. A simple strategy for constructing bounded convection schemes for unstructured grids. Internat. J. Numer. Methods Fluids 46 (10), 1007–1024. Wright, J.A., Smith, R.W., 2001. An edge-based method for the incompressible Navier–Stokes equations on polygonal meshes. J. Comput. Phys. 169 (1), 24–43. Wu, C., Kinnas, S.A., 2019. A distributed vorticity model for flow past a 3-D hydrofoil. In: SNAME 24th Offshore Symposium. The Society of Naval Architects and Marine Engineers. Wu, J.C., Thompson, J.F., 1973. Numerical solutions of time-dependent incompressible Navier-Stokes equations using an integro-differential formulation. Comput. & Fluids 1 (2), 197–215. Wu, X., Wu, J., Wu, J., 1995. Effective vorticity-velocity formulations for three-dimensional incompressible viscous flows. J. Comput. Phys. 122 (1), 68–82. Xing, L., Wu, C., Kinnas, S.A., 2018. Visve, a vorticity based model applied to cavitating flow around a 2-d hydrofoil. In: Proceedings of the 10th International Symposium on Cavitation (CAV2018). ASME Press. Yao, H., Kinnas, S.A., 2019. Coupling viscous vorticity equation (VISVE) method with openFOAM to predict turbulent flow around 2-D hydrofoils and cylinders. In: The 29th International Ocean and Polar Engineering Conference. International Society of Offshore and Polar Engineers. Zhou, C., Graham, J., 2000. A numerical study of cylinders in waves and currents. J. Fluids Struct. 14 (3), 403–428.

Munson, B.R., Okiishi, T.H., Rothmayer, A.P., Huebsch, W.W., 2014. Fundamentals of fluid mechanics. John Wiley & Sons. Naudascher, E., Rockwell, D., 2012. Flow-Induced Vibrations: an Engineering Guide. Courier Corporation. Nehari, D., Armenio, V., Ballio, F., 2004. Three-dimensional analysis of the unidirectional oscillatory flow around a circular cylinder at low Keulegan–Carpenter and 𝛽 numbers. J. Fluid Mech. 520, 157–186. http://dx.doi.org/10.1017/ S002211200400134X. Nielsen, J., Schwind, R.G., 1971. Decay of a vortex pair behind an aircraft. In: Aircraft Wake Turbulence and its Detection. Springer, pp. 413–454. Nieuwstadt, F., Keller, H., 1973. Viscous flow past circular cylinders. Comput. & Fluids 1 (1), 59–71. Payne, R., 1958. Calculation of viscous unsteady flows past a circular cylinder. J. Fluid Mech. 4, 81–86. Peaceman, D.W., Rachford, H.H., 1955. The numerical solution of parabolic and elliptic differential equations. J. Soc. Ind. Appl. Math. 3 (1), 28–41. Pier, B., 2002. On the frequency selection of finite-amplitude vortex shedding in the cylinder wake. J. Fluid Mech. 458, 407–417. Ploumhans, P., Winckelmans, G., Salmon, J.K., Leonard, A., Warren, M., 2002. Vortex methods for direct numerical simulation of three-dimensional bluff body flows: application to the sphere at re= 300, 500, and 1000. J. Comput. Phys. 178 (2), 427–463. Qian, L., Vezza, M., 2001. A vorticity-based method for incompressible unsteady viscous flows. J. Comput. Phys. 172 (2), 515–542. Richardson, L.F., Gaunt, J., 1927. VIII. The deferred approach to the limit. Phil. Trans. R. Soc. A, contain. pap. math. phys. character 226 (636–646), 299–361. Russell, D., Wang, Z.J., 2003. A Cartesian grid method for modeling multiple moving objects in 2D incompressible viscous flow. J. Comput. Phys. 191 (1), 177–205. Saffman, P.G., Meiron, D., 1986. Difficulties with three-dimensional weak solutions for inviscid incompressible flow. Phys. Fluids 29 (8), 2373–2375. Sarpkaya, T., 1976. Vortex shedding and resistance in harmonic flow about smooth and rough circular cylinders at high Reynolds numbers, Tech. Rep., Monterey, California. Naval Postgraduate School. Sarpkaya, T., 1986. Force on a circular cylinder in viscous oscillatory flow at low Keulegan—Carpenter numbers. J. Fluid Mech. 165, 61–71. Sarpkaya, T., 2002. Experiments on the stability of sinusoidal flow over a circular cylinder. J. Fluid Mech. 457, 157–180. Sipp, D., Lebedev, A., 2007. Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows. J. Fluid Mech. 593, 333–358. Smith, P.A.t., Stansby, P.K., 1988. Impulsively started flow around a circular cylinder by the vortex method. J. Fluid Mech. 194, 45–77.

21