A robust implicit compact scheme for two-dimensional unsteady flows with a biharmonic stream function formulation

Computers & Fluids 84 (2013) 141–163

Contents lists available at SciVerse ScienceDirect

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

A robust implicit compact scheme for two-dimensional unsteady flows with a biharmonic stream function formulation Shuvam Sen a,⇑, Jiten C. Kalita b, Murli M. Gupta c a

Department of Mathematical Sciences, Tezpur University, Tezpur 784 028, India Department of Mathematics, Indian Institute of Technology Guwahati, Guwahati 781 039, India c Department of Mathematics, The George Washington University, 2115 G Street, NW, Washington, DC 20052, USA b

a r t i c l e

i n f o

Article history: Received 9 March 2012 Received in revised form 25 December 2012 Accepted 22 May 2013 Available online 4 June 2013 Keywords: Navier–Stokes equations Transient NACA aerofoil Biharmonic

a b s t r a c t In this paper, we develop a compact, implicit second order temporally and spatially accurate finite difference scheme for unsteady Navier–Stokes (N–S) equations for incompressible viscous flows. The scheme is specifically designed for flows in fluid-embedded body interaction as well as curved regions. As the scheme utilizes the 4th order pure stream function form of the N–S equations, the entire flow field can be described in terms of only one equation and thus avoids the difficulties associated with pressure field and nonphysical vorticity boundary conditions. We carry out a spectral analysis as well as von Neumann stability analysis of this scheme and also provide an algorithm for the flow computation. We use the scheme to simulate the time development of 2D viscous flows of varied physical complexities in different geometrical settings: Taylor-Green decaying vortices, constricted channel, flow past rotating and in-line oscillating circular cylinders, and flow past an elliptic cylinder and symmetric aerofoils with various angles of attack. The results obtained are compared with experimental and numerical results available in the literature. Excellent match is obtained in all cases, establishing the efficiency and the accuracy of the proposed scheme. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction The Navier Stokes (N–S) equations governing the flow of fluids have applications in a wide range of physics and engineering problems. Several important numerical techniques have been developed for solving N–S equations where simulations have been carried out for flows as diverse as the flow inside a complex biological region to the interaction between fluid and different types of embedded bodies. Many problems of interest involve dynamics of coupled fluidembedded body interaction. The most common occurrences are associated with fluid flowing over a bluff body or with the movement of a natural or artificial body. Bluff bodies exist in different shapes and sizes and include the flows past an airplane, automobile, submarine and wind blowing past a bridge, a high-rise building, etc. At low Reynolds number the flow over a bluff body is highly viscous; as the Reynolds number increases beyond a critical value, vortex shedding occurs, resulting in a significant pressure drop on the rear surface of the body. Traditionally, the primitive variable and streamfunction– vorticity formulations of the N–S equations have been the most

popular approaches for computing viscous incompressible fluid flows. A review of fundamental formulations of incompressible N–S equations can be found in [1]. For flows in two dimensions, streamfunction–vorticity formulation is more popular for its computational economy; it requires handling of only two unknowns as opposed to three in the case of primitive variables. Furthermore, it ensures exact satisfaction of mass conservation equation. However, in the case of flows in three dimension, the primitive variable formulation is the preferred one. Though three dimensional counterpart of the streamfunction–vorticity formulation also exists in the literature, one has to deal with six unknowns in six equations here. Another approach is to use the fourth order streamfunction vector equation, where each component results in a biharmonic equation and overall, there are only three variables to deal with:

! ! ! ! ! ! @ 1 ðr2 W Þ þ ½ðr  W Þ  rr2 W  ðr2 W  rÞðr  W Þ ¼ r4 W : @t Re ð1Þ !

Here W is known as the three dimensional streamfunction vector. ~ and the velocity ~ v are respectively expressed as: The vorticity x

! ! x ¼ r2 W and ~ v ¼ r  W;

! ⇑ Corresponding author. E-mail addresses: [email protected] (S. Sen), [email protected] (J.C. Kalita), [email protected] (M.M. Gupta). 0045-7930/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2013.05.016

the details of which may be found in [2]. This approach has been used quite successfully in solving the three dimensional cubic

142

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

Fig. 1. The compact computational stencil for the scheme; while all the nodes are used for w, only the nodes denoted by ‘‘’’ are used for wn and wg.

@ 1 ðr2 wÞ þ ½ðr? wÞ  rr2 w ¼ r4 w: @t Re

100 Exact 3D Exact 4D COM 3D COM 4D Central 3D Central 4D

80 60

Λ

40 20 0 −20 −40

0

0.5

1

1.5

2

2.5

3

kh Fig. 2. Comparison of nondimensional characteristics of compact approximations for third and fourth derivatives.

cavity problem [2]. In two dimensions, this system of three biharmonic equations reduces to one biharmonic evolution equation having the scaler streamfunction as the only dependent variable:

ð2Þ

The vorticity is x = r2w and the velocity field is (u, v) = r\w = (wy,  wx). The main advantages of the biharmonic pure stream function formulation are: (i) avoids difficulties associated with primitive variables mainly pressure, (ii) avoids use of artificial vorticity boundary conditions, (iii) Iterations involve only a single variable. The biharmonic pure stream function formulation for the time dependent N–S system in rectangular planar domains has been in existence [3] for almost two decades. Although compact schemes date back to the seventies [4–8] but their application to the biharmonic form of N–S equations is fairly recent [9–13] and started with the pioneering work of Kupferman [9] in 2001. A close look at the literature shows that computations of flows using Eq. (2) has so far been restricted to only rectangular domains. The objective of the present study is to widen the scope of this formulation to non rectangular physical domains and hence test the efficiency of the pure streamfunction formulation. Note that, many challenging and interesting physical problems involve viscous, incompressible fluid flow in geometrically complex regions. Curvilinear coordinate systems generated to maintain coordinate lines coincident with the boundaries provide a key to the development of finite difference solutions of PDEs on regions with arbitrarily shaped boundaries.

Table 1 Problem 1: L1-error, L2-error, L1-error and Rate of convergence for w (dt  h). Mesh

99

Rate

17  17

Rate

33  33

Rate

65  65

t = 0.8 L1 L2 L1

1.142  103 1.297  103 2.020  103

2.92 2.71 2.12

1.507  104 1.980  104 4.662  104

2.93 2.83 2.63

1.979  105 2.778  105 7.552  105

2.97 2.93 2.83

2.534  106 3.656  106 1.060  105

t = 2.0 L1 L2 L1

7.044  104 7.782  104 1.250  103

3.22 3.03 2.51

7.580  105 9.547  105 2.193  104

2.99 2.88 2.64

9.446  106 1.297  105 3.508  105

2.99 2.94 2.84

1.186  106 1.695  106 4.908  106

t = 4.8 L1 L2 L1

1.192  104 1.322  104 1.962  104

3.19 3.01 2.38

1.310  105 1.643  105 3.765  105

3.05 2.92 2.68

1.586  106 2.169  106 5.863  106

3.01 2.95 2.84

1.973  107 2.811  107 8.141  107

143

S. Sen et al. / Computers & Fluids 84 (2013) 141–163 Table 2 Problem 1: L1-error, L2-error, L1-error and Rate of convergence for u and

v components of velocity (dt  h).

Mesh

99

Rate

17  17

Rate

33  33

u t = 2.0 L1 L2 L1

1.6876  103 1.9358  103 3.8506  103

2.10 2.04 2.11

3.9440  104 4.7151  104 8.9154  104

2.06 2.04 1.99

9.4574  105 1.1450  104 2.2439  104

t = 4.8 L1 L2 L1

3.9171  105 5.0489  105 9.9920  105

2.64 2.64 2.26

6.2765  106 8.0912  106 2.0829  105

2.09 2.04 2.02

1.4764  106 1.9635  106 5.1220  106

v t = 2.0 L1 L2 L1

1.6877  103 1.9358  103 3.8500  103

2.10 2.04 2.11

3.9437  104 4.7146  104 8.9156  104

2.06 2.04 1.99

9.4574  105 1.1450  104 2.2439  104

t = 4.8 L1 L2 L1

3.9174  105 5.0489  105 9.9963  105

2.64 2.64 2.26

6.2765  106 8.0912  106 2.0829  105

2.09 2.04 2.02

1.4764  106 1.9635  106 5.1221  106

3

Table 4 Problem 4: Quantitative comparison of the normalized drag coefficient.

2

y

1 Exact 9 9 17 17 33 33 65 65

0

Cases

Experimental CD/CD0, [23]

Numerical CD/CD0, [24]

Numerical CD/CD0, Present

A = 2, ff/f0 = 0.75 A = 2, ff/f0 = 1.5 A = 2, ff/f0 = 3.0

– 1.155 0.869

– 1.245 0.855

2.338 1.077 0.846

-1

ð5Þ

v

ð6Þ

-2 -3 -0.1

-0.05

0

ψ

0.05

0.1

1 ðw xn  wn xg Þ; J g 1 ¼ ðwn yg þ wg yn Þ; J

u¼

where

We assume that the physical (x, y) plane can be transformed into a computational (n, g) plane by a conformal transformation of the form:

z ¼ zðhÞ;

  1 1 2C ; bðn; g; wn ; wg Þ ¼   wg ; ReJ J Re     1 2D 1 E  þ wn ; dðn; g;wn ;wg Þ ¼ þ Cwg  Dwn ; cðn; g;wn ;wg Þ ¼ J Re J Re Jg J Jn J gg C ¼ ; D ¼ ; E ¼ 2C 2 þ 2D2   nn ; J J J J

aðn; gÞ ¼

Fig. 3. Problem 1: grid convergence along the line x = 0.

ð3Þ

z = x + iy and h = n + ig, and Eq. (2) in the computational plane takes the form:

@ Dw ¼ aðn; gÞD2 w þ bðn; g; wn ; wg ÞDwn þ cðn; g; wn ; wg ÞDwg @t þ dðn; g; wn ; wg ÞDw; with the velocities u and

v given by

ð4Þ

J being the Jacobian of the transformation (3) and D  r2. Details of derivation of Eq. (4) are available in the Appendix. Eq. (4) is the fourth order semi-linear stream function form of the N–S equation in the transformed plane. This PDE can be used to simulate fluid flow in any region which is conformally mappable to a rectangular domain. This formulation can also be used in cases where the domain can only be conformally mapped numerically; there exists a large class of domains for which one can numerically

Table 3 Error analysis of the flow through constricted channel.

s = 0.6 Re = 100

s = 0.8 Re = 10

Grid size

Error in L1

[471  21] [706  31] [941  41]

kW2  W1k = 8.388  105

[471  21] [706  31] [941  41]

kW2  W1k = 6.175  106

Rate (a)

Error in L2

2.131 kW3  W2k = 2.780  105

Error in L1

2.178

2.252 kW3  W2k = 5.123  104

kW2  W1k = 3.020  105

kW2  W1k = 3.099  104 2.097

kW3  W2k = 1.020  105

Rate (a)

kW2  W1k = 1.603  103

kW3  W2k = 6.762  105 2.086

kW3  W2k = 2.095  106

Rate (a)

kW2  W1k = 2.060  104

2.112 kW3  W2k = 1.042  104

144

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

(a)

(b) Fig. 4. Problem 2: (a) schematic diagram with s = 0.8, (b) close up view of a typical 941  41 grid.

1 0.5 0 -0.5 -1 -5

0

5

10

15

20

10

15

20

(a) 1 0.5 0 -0.5 -1 -5

0

5

(b) 1

generate a conformable mapping [14]. Eq. (4) contains w as the only dependent variable from which other flow variables can be post processed. The major advantage here is that this formulation bypasses difficulties that often come up while dealing with vorticity in the streamfunction–vorticity formulation or pressure and mass conservation equation in the primitive variable formulation. In this paper, we propose a new compact finite difference scheme for Eq. (4) which is second order accurate both temporally and spatially. Here we not only document the versatility of the fourth order stream function formulation, but also the efficiency of the newly proposed scheme in simulating the dynamics of flow inside curved regions as well as fluid-embedded body interaction. We use the scheme to accurately capture the details of flows consisting of large band of wave-numbers while preserving its physical dispersiveness. We use the new finite difference scheme to simulate time development of 2D viscous flows in different geometrical settings and test it on the following fluid flow problems of varied geometrical settings: (i) Taylor-Green decaying vortices, (ii) constricted channel, (iii) flow past a rotating cylinder, (iv) inline oscillating cylinder in a fluid at rest, (v) flow past an elliptic cylinder with angle of attack, (vi) flow past symmetric aerofoils with angle of attack. The first problem of Taylor-Green decaying vortices has a known closed form solution, which allows us to compare our numerical solutions with the exact ones and establish the rate of convergence of our method. The second problem of flow through constricted channel is an internal flow problem, where we time march to steady state and it serves as a test case for convergence analysis of fluid flow problems without known analytical solutions. Problems 3–6 represent flow over bluff bodies of different curvatures: the flow past circular cylinders having constant curvature, elliptical cylinders having varied curvature and aerofoils having infinite curvature at the trailing edge. For the circular cylinder, two different physical scenarios are considered. We consider the cases of rotating cylinder with constant and oscillatory rotations, and the case of inline oscillatory cylinder. All these rotatory and oscillatory cases symbolize moving boundaries and in each case, we obtain highly accurate solutions. The robustness of our scheme is exemplified by the simulation of flow past elliptical cylinders and symmetric aerofoils with different angles of attack. This is more clearly evident for the aerofoil where we have precisely simulated both uniform and accelerated flows for a wide range of Reynolds numbers and angles of attack. In all cases, our results are validated through comparisons with existing qualitative and quantitative results, both experimental and numerical, and excellent agreement is obtained in all cases. The outline of the paper is as follows. Section 2 deals with discretization procedure and issues related to it, Section 3 with the solution of the algebraic system of equations and the related algorithm, Section 4 with the numerical results and finally in Section 5, we offer a few conclusions.

0.5

2. Discretization procedures

0 -0.5 -1 -5

0

5

10

15

20

(c) Fig. 5. Steady state stream lines for (a) s = 0.6, Re = 100; (b) s = 0.6, Re = 1000 and (c) s = 0.8, Re = 1000.

Amongst the several approaches of discretizing and numerically solving the N–S equations, we are interested in finite difference approximation because of its relative ease in implementation. Our aim now is to discretize Eq. (4) compactly using only the eight neighboring grid points surrounding each node. The main advantages of such an approach will be twofold: (i) the coefficient matrix resulting from the system of linear equations after discretization will have a smaller bandwidth and (ii) no modification will be

145

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

1

1

0.5

0.5

0

0

-0.5

-0.5

-1

-1 -5

0

5

(a)

10

15

20

1

1

0.5

0.5

0

0

-0.5

-0.5

-1

-1 -5

0

5

(c)

10

15

20

1

1

0.5

0.5

0

0

-0.5

-0.5

-1

-5

0

5

-5

0

5

-5

0

5

10

15

20

10

15

20

10

15

20

(b)

(d)

-1 -5

0

5

(e)

10

15

20

(f)

Fig. 6. Time evolution of stream line contours for flow through constricted channel Re = 500, s = 1.0 at (a) t = 1.0, (b) t = 2.0, (c) t = 4.0, (d) t = 6.0, (e) t = 8.0, (f) t = 15.0.

needed at grid points near the boundaries, thus facilitating easier and at times exact implementation and satisfaction of the boundary conditions. The approach here involves discretizing the biharmonic Eq. (4) using not just the grid values of the unknown solution w but also the values of the gradients wn and wg at some grid points inside the computational stencil as shown in Fig. 1. This approach has its roots in the works of Gupta and Manohar [15] and was formulated by Stephenson [16]; later on, significant contributions in this direction were made by Altas et al. [17]. We use the following discretizations for space derivatives D2w and Dw appearing in (4): 2

D2h wi;j ¼ d4n wi;j þ 2d2n d2g wi;j þ d4g wi;j þ Oðh Þ; 2

Dh wi;j ¼ d2n wi;j þ d2g wi;j þ Oðh Þ:

ð7Þ ð8Þ

Here d4n denotes the Stephenson finite difference operator [16] defined as

d4n wi;j 

12  2

h

 dn wni;j  d2n wi;j ;

where d2n and dn are the usual second order central difference operators for second and first order derivatives respectively. Note that the presence of wn and wg in (7) may give the impression of complicating the associated linear system, instead, in the present case, it is an added advantage as these gradients, which are explicitly present in Eq. (4), being already available at all grid points, need not be approximated from the computed values of the solution w. Compatible Padé approximations for wn and wg may be chosen as:

2

!

wni;j ¼

dn wi;j 

h 2 d w 6 n ni;j

wgi;j ¼

dg wi;j 

h 2 d w 6 g ni;j

2

4

ð9Þ

4

ð10Þ

þ Oðh Þ; ! þ Oðh Þ:

Having obtained a second order approximations for space derivatives of Eq. (4), our next step is to discretize time derivative and obtain a stable numerical scheme. We introduce weighted time ðnÞ ðnþ1Þ average parameter k such that t k ¼ ð1  kÞt k þ ktk for 0 6 k 6 1, and approximate Eq. (4) as:

h ðnþ1Þ ðnÞ ðnÞ ðnÞ Dh wi;j ¼Dh wi;j þ dtð1  kÞ ai;j D2h wi;j þ bi;j Dh wni;j þ ci;j Dh wðnÞ gi;j ðnÞ

þ di;j Dh wi;j

i

h ðnþ1Þ ðnþ1Þ þ dtk ai;j D2h wi;j þ bi;j Dh wni;j þ ci;j Dh wðnþ1Þ gi;j

ðnþ1Þ

þ di;j Dh wi;j

i

:

ð11Þ

Varying k provides a class of integrators, for example forward Euler for k = 0, backward Euler for k = 1 and Crank–Nicolson for k ¼ 12. Using the last option, we obtain an O(h2, dt2) accurate scheme for (4). Note that whenever the need for computing the velocities arises, Eqs. (5) and (6) in conjunction with Eqs. (9) and (10) can be used. The finite difference approximation (11) developed here has the following two major benefits as compared to the compact or any other method based on vorticity–streamfunction formulation:

146

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

6

Drag/Lift

4

(b) t=3.5

(a) t=3.5

Drag, α =1 Lift, α =1 Drag, α =3 Lift, α =3

2

0

(d) t=15.0

(c) t=8.0

-2 0

5

10

15

Time

(a) (a) 1.5 (b) t=3.0

y

(a) t=1.5

1 0.5 0

0

Numerical Present Experimental Badr et al. Numerical Chew et al.

1

2

3

4

x

(b) Fig. 7. (a) Evolution of drag and lift coefficients for Re = 1000 and a = 1,3, (b) Comparison between the paths of the first vortex obtained numerically by our computation with experimental [21] and numerical [22] works.

1. Apart from the initial condition, stream function and velocity boundary conditions are sufficient to carry out numerical computation. The formulation entirely avoids the use of artificial vorticity boundary conditions. 2. A single discretized biharmonic equation involving only streamfunction is used here as compared to coupled system of equations involving both streamfunction and vorticity.

2.1. Fourier analysis of finite difference approximations To examine in detail the characteristics of the compact scheme developed above, a spectral analysis of the finite difference approximations of the higher order derivatives viz. third and fourth order derivatives have been carried out. Note that truncation error does not represent all the characteristics of a finite difference approximation and modified wave number analysis provides a tool to compare the resolution of finite difference operator with its analytical counterpart as well as other difference approximations. Here we carry out spectral analysis of the finite difference approximations used for the third and fourth order derivatives found in Eq. (4).  pffiffiffiffiffiffiffi Applying the trial function eIkn I ¼ 1 where k is the wave number, we see that the exact characteristic of the third order derivative wnnn and fourth order derivative wnnnn as

(d) t=5.0

(c) t=5.0

(b) Fig. 8. Computed and experimental [21] stream lines at different times for Re = 1000 (a) a = 1 (b) a = 3.

K3D Exact ¼ Ik

3

4

and K4D Exact ¼ k ;

respectively. These derivatives have been compactly approximated by using the relations d4n wi;j and d2n wni;j respectively. Replacing these difference operators with the corresponding modified wave numbers, the characteristic functions obtained are given by

 6 sinðkhÞðcosðkhÞ  1Þ and 3 2 þ cosðkhÞ h " # 2 12 2  2 cosðkhÞ  sin ðkhÞ ; ¼ 4 2 þ cosðkhÞ h

K3D COM ¼ K4D COM

I



respectively. In order to compare we write below the characteristics of regular wide molecule second order central difference approximations of the third and fourth order derivatives

147

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

6

where a > 0, b, c, d are considered to be independent of w and its derivatives. Note that a(n,g) > 0 in Eq. (4). The scheme (11) with k ¼ 12 when applied to Eq. (12) gives

Ratio=0.75 Ratio=1.5 Ratio=3.0

ðnþ1Þ

ð1  jÞDh wi;j

4

2

ðnþ1Þ

 ah D2h wi;j 2

ðnÞ

ðnþ1Þ

 bDh wni;j

ðnÞ

 cDh wðnþ1Þ gi;j

ðnÞ

Drag/Lift

¼ ð1 þ jÞDh wi;j þ ah D2h wi;j þ bDh wni;j þ cDh wðnÞ gi;j ;

ð13Þ

2 adt bdt cdt with a ¼ 2h 2 ; b ¼ 2 ; c ¼ 2 ;

j ¼ ddt . 2

Theorem 1 (Stability). The finite difference scheme (13) is unconditionally stable, in the von Neumann sense, if d 6 0 and for d > 0 the scheme is stable under sufficient condition dt < 2d.

0

-2 ðnÞ

0

10

20

30

40

50

ðnÞ

Proof. Let wi;j ¼ b eIhn i eIhg j where b(n) is the amplitude at time pffiffiffiffiffiffiffi level n; I ¼ 1 and hn = 2ph/K1, hg = 2ph/K2 are the phase angles

60

Time

with wavelengths K1 and K2 respectively; then from relations (9) and (10) we get

(a)

ðnÞ

wni;j ¼ I

3 sin hn ðnÞ w ; hð2 þ cos hn Þ i;j

ð14Þ

3 sin hg ðnÞ w : hð2 þ cos hg Þ i;j

ð15Þ

and

wðnÞ gi;j ¼ I

(ii)

(i)

If G is the amplification factor associated with the finite difference Eq. (13) then

G¼

where

(iii)

(b)

A ¼ ð2  cos hn  cos hg Þ;

Fig. 9. (a) Evolution of drag () and lift (  ) coefficients for Re = 150, A = 2 with different ff/f0 ratio, (b) computed results (i) are compared with numerical results of Sengupta and Bhumkar [24] (ii) and experimental results [23] (iii) for flow past a circular cylinder executing rotary oscillation for Re = 150, A = 2, ff/f0 = 1.5.

K3D CD ¼ I K4D CD ¼

2 h

3

B ¼ 14  8ðcos hn þ cos hg Þ þ 2 cos hn cos hg ! 2 2 sin hg sin hn 9 þ ; 2 þ cos hn 2 þ cos hg C¼

sinðkhÞðcosðkhÞ  1Þ and

2 h

ð1 þ jÞA  2aB þ IC ; ð1  jÞA þ 2aB  IC

½cosð2khÞ  4 cosðkhÞ þ 3: 4

The non-dimensional characteristics of the derivatives as a function of kh have been depicted in Fig. 2. The figure clearly depicts that in comparison to second order wide molecule central difference scheme our compact formulation possesses better spectral resolution characteristic. Comparison of the resolution of fourth order derivative clearly outlines the fact that the present compact scheme has much smaller dissipative error. Similarly comparison of the characteristics of the third order derivative shows that the scheme has less dispersive error. The above discussion suggests that the compact algorithm developed here possesses better resolution properties and may be best suited for high Reynolds number flow computation. 2.2. Stability analysis We perform a von Neumann stability analysis to the linearized version of Eq. (4)

@ Dw ¼ aD2 w þ bDwn þ cDwg þ dDw; @t

ð12Þ

  sin hg 3 sin hn bþ c A: h 2 þ cos hn 2 þ cos hg

For j 6 0 ) d 6 0,jGj 6 1 and the scheme is stable. For j > 0 ) d > 0 so the fourth term on right hand side of (12) is a growth term and the stability condition requires jGj < 1 + Kdt for some K independent of dt. This is obviously true for j < 1 ) dt < 2d. This completes the proof. h Note that the von Neumann stability analysis presented here is a representative one. Because when applied, it holds for uniform meshes, periodic boundary conditions for linearized systems. In practice, for nonlinear differential equations, this analysis is not valid as many of the assumptions are violated; we still use the analysis on localized linear versions of differential equations in order to obtain some guidance towards the appropriate values of dt. As such the stability limits are much lower. Also bigger time step and grid spacing affect accuracy of the solutions obtained. For example, for the problem of flow past stationary cylinder considered in section Numerical Results, at lower Re, a grid of size 61  101 with dt = 0.01 is sufficient to obtain a stable numerical solution while for Re = 5000 this choice of grid suffers from stability consideration and we are required to chose 181  301 grid with dt = 0.001. Similarly for the problem of flow past elliptic cylinder it will be seen that for the flow with Re = 200 at angle of attack 45°, dt = 0.001 yields very good simulation, but for Re = 1000, when the angle of attack is increased to 60°, dt was reduced to 0.0001. Likewise, for

148

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

4

4

2

2

0

0

-2

-2

-4

-4 0

5

10

15

0

5

10

15

0

5

10

15

0

5

10

15

(a) 4

4

2

2

0

0

-2

-2

-4

-4 0

5

10

15

(b) 4

4

2

2

0

0

-2

-2

-4

-4 0

5

10

15

(c) Fig. 10. Stream lines (left) and vorticity contours (right) for a cylinder performing oscillatory rotation (the dotted lines representing negative values of the contours) at Re = 150, A = 2 and frequency ratio ff/f0= (a) 0.75, (b) 1.5, (c) 3.0.

the flow past NACA 0012 at Re = 3000 and angle of attack 45°, a grid of size 201  121 and dt = 0.005 is effective, but while simulating accelerated flow past NACA 0015 for Re = 5200 at angle of attack 90° we have used dt = 0.0001.

By setting k = 0 in Eq. (11) we get a first order time accurate formula which has the matrix representation

  ðnÞ A4 Wðnþ1Þ ¼ F 4 WðnÞ ; Wn ; WðnÞ g :

ð19Þ

3. The solution of algebraic system of equations

Here A4 = A1 and we have the advantage that W(n+1) can be estimated directly. Thus we have the following algorithm

Let us now discuss the solution of algebraic systems associated with the finite difference approximation (11). As we are interested in second order accuracy in time we chose k ¼ 12. The resulting system of equation in matrix form can be written as:

1. 2. 3. 4.

  ðnÞ ðnþ1Þ : A1 Wðnþ1Þ ¼ F 1 WðnÞ ; Wn ; WðnÞ ; Wðnþ1Þ g ; Wn g

ð16Þ

ðnþ1Þ

Predict Wold using (19). ðnþ1Þ Predict Wnold ; Wðnþ1Þ gold using (17) and (18) respectively. Correct to Wðnþ1Þ new using (16). ðnþ1Þ Correct to Wnnew ; Wðnþ1Þ using (17) and (18) respectively. g new ðnþ1Þ ðnþ1Þ 5. If Wnew  Wold <  then Wðnþ1Þ ¼ Wðnþ1Þ new . ðnþ1Þ

For a grid of size m  n, the matrix A1 has the dimension mn. Due to the compact nature of our scheme A1 is a banded matrix with nine ðnÞ ðnþ1Þ non zero diagonals. Also Wðnþ1Þ ; WðnÞ ; Wn ; WðnÞ ; Wðnþ1Þ are g ; Wn g (n) all mn component vectors. At any time step, once W has been ðnÞ approximated Wn ; WðnÞ g can be obtained by solving tridiagonal systems ðnÞ

A2 Wn ¼ F 2 ðWðnÞ Þ;

ð17Þ

ðnÞ A3 WðnÞ g ¼ F 3 ðW Þ;

ð18Þ

respectively. Eqs. (17) and (18) are the corresponding matrix forms of the relations (9) and (10). Thus the main objective now is to solve Eq. (16), thereby evaluating unknown vector W(n+1). But a difficulty arises due the presence of (n + 1)th time level gradients of W on the right hand side of Eq. (16) as those quantities will be available only after solving for streamfunction at the (n + 1)th time level. To overcome this difficulty we adopt a predictor–corrector approach.

6. Wold

¼ Wðnþ1Þ new , goto step 3.

Direct solution of any of the above linear system is impractical because of huge size of the coefficient matrix and enormous storage requirements even for moderate values of step length h. On the other hand condition number of the coefficient matrix increases rapidly with reduced step length h and one must be very cautious while attempting to solve such linear systems using iterative solvers. As the coefficient matrix A1 is not diagonally dominant, conventional solvers such as Gauss–Seidel also cannot be used. Therefore all the computations were performed using the biconjugate gradient stabilized (BiCGStab) method without preconditioning, where, thanks to the compact grid, it is easy to implement matrix vector multiplication A1W without the need of storing all the entries of the matrix A1. The convergence criterion for BiCGStab iteration based on norm of residual was set at 108 and the stopping criterion for the

149

1

1

0.5

0.5

y

y

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

0

0

-0.5

-0.5

-1

-1 -1

0

1

-1

(a)

1

1

0.5

0.5

0

1

0

-0.5

-0.5

-1

-1 -1

0

1

u

-1

(b)

1

1

0.5

0.5

0

0

1

v

y

y

0

v

y

y

u

0

-0.5

-0.5

-1

-1 -1

0

1

u

-1

(c)

0

1

v

Fig. 11. Comparison of the velocity components with experimental results [27] at four cross-sections with constant x = 0.6D(, j), 0D(, N), 0.6D(, ), 1.2D(. . ., ). Phase angles (a) 90°, (b) 120°, (c) 240°.

corrector was set at 1012. All our computations were carried out on a Pentium Dual-Core processor based PC with 4 GB RAM using double precision floating point arithmetic. 4. Numerical results 4.1. Taylor-Green decaying vortices: test for accuracy At first we consider the Taylor-Green problem of flow decayed by viscosity. This problem is often used as a model problem for analyzing the spatial accuracy of schemes [11,13,18]. The analytical solution of this problem expressed in terms of stream function is: 2t

w ¼  cosðxÞ cosðyÞeRe : We solve this problem inside the domain D ¼ fðx; yÞ : xþ pffiffiffiffiffiffiffi iy ¼ ð1 þ iÞðn þ igÞ;  p=2 6 n 6 p=2; p=2 6 g 6 p=2; i ¼ 1g,

with Re = 1. The initial and boundary conditions have been derived from the analytical solution. We carry out calculation with four different grids having sizes 9  9, 17  17, 33  33 and 65  65 using dt = 0.4, 0.2, 0.1 and 0.05 respectively. The choice of dt is motivated by our desire to keep dt  h so as to check the second order accuracy of the scheme. Results obtained for t = 0.8, t = 2.0 and t = 4.8 have been tabulated in Tables 1 and 2. In Fig. 3 we present a comparison between our numerical solution and the exact solution along the x = 0 line on four grids of increasing sizes. As can be seen from the figure, grid independence is obtained on a grid of size as coarse as 17  17. The L1, L2 and L1 errors in stream function field have been presented in Table 1. The table shows that the present scheme achieves a better rate than the expected second-order spatial accuracy which may be attributed to the nature of the test problem being solved [11,13]. Since in practical simulation, the order of accuracy for velocity is more important than that for stream function we provide the order of accuracy for u and v components of velocity in Table 2. Here again the scheme shows a second order

150

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

2

2

1

1

5 Re=15, Drag Re=15, Lift Re=30, Drag Re=30, Lift

0

0

-1

-1

-2

-2

-1

0

1

2

-2

-2

-1

0

1

Drag/Lift

4

2

3

2

(a) 2

2

1

1

1

0 0

0

-1

-1

0

2

4

6

8

10

Time

(a) 1.75 -2

-2

-1

0

1

2

-2

-2

-1

0

1

2

(b) 2

1

1

0

0

-1

-1

1.5

Drag/Lift

2

Drag Lift

1.25

1

0.75 -2

-2

-1

0

1

2

-2

-2

-1

0

1

2

(c) 2

2

1

1

0

0

-1

-1

-2

-2

-1

0

1

2

-2

0.5

0

10

20

30

40

50

60

Time

(b) Fig. 13. Evolution of drag and lift for (a) Re = 15 and Re = 30 at k = 0.1, h = 45°; (b) Re = 163 at k = 0.2, h = 45°.

Table 5 Comparison of steady state drag and lift coefficients for elliptic cylinder with h = 45° and k = 0.1.

-2

-1

0

1

2

(d) Fig. 12. Stream lines (left) and vorticity (right) contours for in-line oscillating cylinder (the dotted lines representing negative values of the contours) at Re = 100, KC = 5 for four different phase angles (a) 6°, (b) 102°, (c) 198°, (d) 270°.

accuracy. The problem being symmetric, the errors in u and v show identical behavior. 4.2. The problem of constricted channel We consider this problem because apart from being a representative of flow problems in confined domains this is an idealization of flows in geometries that contains a re-entrant corner [19,20]. At such a corner the flow becomes singular and in order to obtain

Re

Ref. [31]

Ref. [35]

Present

15

CD CL

2.135 1.310

1.865 1.050

2.484 1.316

30

CD CL

1.430 0.935

1.402 0.930

1.829 1.164

accurate results near such a corner, it is necessary to have a highly refined non-uniform grid there. In this problem the flow occurs in a non-uniform channel containing a step-down constriction. Another motive here is to establish rate of convergence for problems having no known analytical solution. The shape of the channel boundary can vary from one having a smooth constriction to one with a very sharp corner. We transform the physical channel into a rectangular computational plane by using the algebraic conformal mapping:

151

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

0.6 0.3 0 -0.3 -0.6

-1

0

1

2

3

(a) 0.6 0.3 0 -0.3 -0.6

-1

0

1

2

3

(b) Fig. 14. Comparison of steady state stream lines from the present computation (left) with those obtained by Dennis and Young [35] (right) (a) for Re = 15, k = 0.1, h = 45°, (b) for Re = 30, k = 0.1, h = 45°.

Table 6 Time evolution of drag and lift coefficients for elliptic cylinder for Re = 200 and h = 45°. (The values in parenthesis are from [31].) k = 0.2 t CD CL k = 0.1 t CD CL

0.84 1.552 (1.565) 1.768 (1.905)

1.585 1.605 (1.735) 1.752 (2.030)

2.395 1.471 (1.508) 1.273 (1.440)

3.205 1.271 (1.350) 0.631 (0.575)

4.42 1.298 (1.340) 0.617 (0.585)

8.9 1.162 (1.230) 0.854 (0.810)

10.15 1.272 (1.505) 1.327 (1.545)

11.2 1.303 (1.565) 0.945 (1.280)

11.6 1.331 (1.550) 0.867 (1.060)

12.05 1.343 (1.485) 0.829 (0.910)

As the positioning of the inlet and outlet boundaries has a direct bearing on the computed flow [19,20], we fix the distances of the inlet and outlet at x  7 and x  20 respectively. A close up view of the grid near the constriction for a typical grid of size 941  41 is shown in Fig. 4b. The ratio between the length and the breadth of the channel needs to be kept sufficiently large so as to allow the flow to be fully developed at the outlet and as such there are many more grid points along ordinate axis than abscissa. The two main controlling parameters for this flow are s and Re, and several combinations of them have been tested in our computation. For all the computations we use dt = 0.001. In Table 3, we present the converged steady state error results for the cases s = 0.6, Re = 100 and s = 0.8, Re = 10 showing the spatial rate of convergence a. This table clearly reveals a convergence rate of order two for the method. Note that in the absence of an analytical solution, the rate of convergence a can be estimated using the relation: a

x ¼ An þ

B ½n sinhð2nÞ  g sinð2gÞ; H

B y ¼ Ag þ ½g sinhð2nÞ þ n sinð2gÞ; H where H ¼ coshð2nÞ þ cosð2gÞ; A ¼ bþa , and B ¼ ba . The constants a 2s 2s and b are the radii of the channel far upstream and far downstream respectively. In our computation, we take a = 1.0 and b = 0.5. In addition, 0 < s 6 1 controls the sharpness of the constriction, with a higher value of s indicating a sharper corner in the domain of the flow. A schematic diagram of the problem is presented in Fig. 4(a). The Reynolds number for this problem is defined as Re ¼ aum0 , where uo is the average velocity at the entrance and m is the kinematic viscosity. A Poiseuille parabolic velocity profile is introduced at the inlet and outlet. These boundary conditions are designed to satisfy conservation of mass in the computational domain. At the walls, a no slip condition is assumed. The boundary conditions can be summarized as:

@w w ¼ 1 and ¼ 0 on the wall g ¼ s; @g @w w ¼ 1 and ¼ 0 on the wall g ¼ s; @g   2 g g @w and 3 2 ¼ 0 as n ! 1: w¼ @n 2s s

a

kW2  W1 k : h1  h2 ¼ kW3  W2 k ha2  ha3 where h1, h2, h3 are mesh sizes in three different uniform grids used and W1, W2, W3 respectively are computed solutions on these grids. We show the streamline contours for Re = 100, s = 0.6 and for Re = 1000, s = 0.6 and 0.8 in Fig. 5. This figure (as well as Fig. 6) clearly indicates that vortices near the constriction are formed only for high values of s . For the case s = 1.0, Re = 500 we present the evolution of streamlines in Fig. 6 where one can see the development of the main vortices at the top and bottom near constrictions. To begin with, there are two small vortices which grow in size as time elapses. As the main vortex acquires strength, smaller vortices are shed. These vortices propagates towards the outlet and diffuse completely before the steady state is reached. 4.3. Flow past a rotating cylinder In this section, we investigate the flow past a circular cylinder which starts translating and rotating impulsively from rest in a viscous fluid. Two different types of oscillations are considered. In the first case the cylinder is considered to be rotating with constant angular velocity while in the second case we consider rotatory oscillation of the cylinder. Both these cases have importance of their own. The first one exhibits the suppression of the separated vortex where studies have shown that for high enough velocity

152

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

2

2

1

1

0

0

-1

-1

-2

0

2

4

6

8

10

-2

(a)

2

2

1

1

0

0

-1

-1

-2

0

2

4

6

8

10

-2

(b)

2

2

1

1

0

0

-1

-1

-2

0

2

4

6

8

10

-2

(c)

2

2

1

1

0

0

-1

-1

-2

0

2

4

6

8

10

-2

(d)

2

2

1

1

0

0

-1

-1

-2

0

2

4

6

8

10

-2

(e)

0

2

4

6

8

10

0

2

4

6

8

10

0

2

4

6

8

10

0

2

4

6

8

10

0

2

4

6

8

10

Fig. 15. Stream lines (left) and vorticity contours (right) depicting the flow past elliptic cylinder with Re = 163, k = 0.2 and h = 45° at (a) t = T, (b) t ¼ T þ T40 , (c) t ¼ T þ T20 , (d) t ¼ T þ 3T40 and (e) t = T + T0.

ratios, steady flows with no vortex shedding is possible even for high Re values [21,22]. The second case has the richness of the vortex structure associated with such flows. Here the rotary oscillation of the cylinder causes separation of the flow from the cylinder surface, interrupting the formation of regular von Kármán vortices in the wake, leading to the creation of small-scale vortices at the cylinder surface [23,24]. It is a true test of efficiency for any numerical scheme to be able to capture these small physical vortices and to record how strongly the vorticity distribution is affected in the near wake by the forcing frequency. We assume the cylinder to be of unit radius placed in an infinite domain. A uniform grid spacing is employed along the cross radial direction and nonuniform grid spacing in the radial direction with clustering around the surface of the cylinder using the transformation

x ¼ eðpnÞ cosðpgÞ;

y ¼ eðpnÞ sinðpgÞ:

The flow variables are non-dimensionalized by using x = x⁄/a, y = y⁄/a, u = u⁄/U1, v = v⁄/U1, t = t⁄U1/D and the Reynolds number Re is defined as Re = DU1/m where D = 2a is the diameter of the cylinder; here the asterisk sign ⁄ denotes the dimensional quantities. The boundary conditions for stream function and its first order derivatives are obtained as: 1. On the surface of the cylinder we impose rotational velocity. 2. In the far upstream u = U1, v = 0 ) wn = pe(pn) sin(pg), wg = pe(pn) cos(pg) and w = e(pn) sin(pg) which corresponds to the potential flow. 3. In the far downstream we use the Neumann condition that @V r t ¼ 0 ¼ @V . @n @n

153

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

0.6

0.6

(b) t=0.5

(a) t=0.25 0.4

0.4

0.2

0.2

0

0

-0.2

-0.2

-0.4

-0.4

-0.6 -0.4

-0.2

0

0.2

0.4

0.6

0.8

-0.6 -0.4

-0.2

0

0.2

0.4

0.6

0.4

0.2

0.2

0

0

-0.2

-0.2

-0.4

-0.4

-0.2

0

0.2

0.4

0.6

0.8

-0.6 -0.4

(i)

-0.2

0

0.2

0.4

0.8

(b) t=22.0

(a) t=20.0 0.4

0.4

0.2

0.2

0

0

-0.2

-0.2

-0.4

-0.4

-0.2

0

0.2

0.4

0.6

0.8

-0.6 -0.4

-0.2

0

0.2

0.4

0.6

0.8

0.6

0.6

(d) t=26.0

(c) t=24.0 0.4

0.4

0.2

0.2

0

0

-0.2

-0.2

-0.4

-0.4

-0.6 -0.4

0.6

0.6

0.6

-0.6 -0.4

0.8

(d) t=1.5

(c) t=1.0 0.4

-0.6 -0.4

0.6

0.6

-0.2

0

0.2

0.4

0.6

0.8

-0.6 -0.4

(ii)

-0.2

0

0.2

0.4

0.6

0.8

Fig. 16. (i) Early evolution of wake and (ii) wake structure during shedding behind elliptic cylinder for Re = 1000 k = 0.1, h = 60°.

4. In the far downstream we use convective boundary condition for w given by @w þ U 1 @w ¼ 0 appropriately translated onto the @t @x n-g plane.

For initial conditions we use w = 0 everywhere except at the outer boundary. For the scheme developed here in the context of Neumann boundary conditions, one sided second order

154

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

0.3 Actual Numerical

0.2

y

0.1 0 -0.1 -0.2 -0.3

0

0.2

0.4

0.6

0.8

1

x

(a) 100 NACA0012 used Joukowski NACA0012 NACA0015 used Joukowski NACA0015

% of error

80

60

40

20

0

0

0.2

0.4

0.6

0.8

1

Aerofoil axis

(b) 0.5

-0.99

0

-0.995

-0.5

different flow fields for Re = 1000 with a = 1 and 3; as pointed out in [21], the two flow fields have distinctly different characters. The evolution of drag and lift coefficients for both the cases are shown in Fig. 7a. From the figure it is clear that for a = 1.0, both drag and lift assume a periodic nature; but for a = 3.0, steady value is reached. To provide a qualitative comparison between our computed and the experimental results of Badr et al. [21] and numerical results of Chew et al. [22], we plot the path of the first vortex centre for Re = 1000, a = 1.0 in Fig. 7b. An excellent comparison can be seen. Note that our computation offer better result vis-a-vis the numerical work of Chew et al. [22] for a large duration of motion. The streamline profiles at different instants of times along with some of the experimental visualizations of [21] for both the cases are provided in Fig. 8. Once again, one finds an excellent match between the computed and experimental streamlines. 4.3.2. Case 2: oscillatory rotation The angular velocity x ¼ h_ of the cylinder about its axis in this case is characterized by two non-dimensional parameters (i) forcing amplitude (A) and (ii) frequency ratio (ff/f0) [23] and is given by h(t) = h0 cos(2pfft). Here ff is called the forcing frequency, f0 is the shedding frequency of the stationary cylinder at same Re and að2pf h Þ A ¼ U1f 0 . Note that some authors prefer to use Strouhal ratio (Stf/St0) as the control parameter [26] instead of the forcing ratio (ff/f0). For oscillatory rotation we consider three different cases of frequency ratio ff/f0 = 0.75, 1.5, 3.0 while keeping a fixed value of Re = 150 and A = 2. The time variation of drag and lift coefficients for all the three cases leading to the periodic state is presented in Fig. 9a. Further, in Fig. 9b, a qualitative comparison of computed vorticity contours for the case Re = 150, A = 2 and ff/f0 = 1.5 with numerical results [24] and experimental results [23] clearly establishes the accuracy of the numerical results obtained by us. The computed streamlines and the post processed vorticity contours at time t = 40.0 for all the three cases are depicted in Fig. 10. The vorticity contours obtained in our simulation matches well with the detailed analysis carried out by Protas and Wesfreid [26]. Next we carry out a quantitative comparison for our cases in Table 4 with the experimental works of Thiria et al. [23] and numerical works of Sengupta and Bhumkar [24]. We have normalized the mean drag coefficient CD of a particular rotary oscillation case by the mean drag coefficient C D0 of a stationary cylinder at the same Reynolds number. Again a strong agreement of the data can be seen here.

-1

4.4. In-line oscillating cylinder in a fluid at rest

-1 -1.005 -0.005

-1.5 -0.5

0

0.5

0

1

0.005

1.5

(c) Fig. 17. (a) Comparison with actual NACA 0012 aerofoil, (b) comparison of error percentage, (c) typical 301  121 grids around NACA 0015 having angle of attack = 90°.

approximations have been used [25]. Identical approach have been followed for other problems discussed in this manuscript as well. 4.3.1. Case 1: constant rotation As pointed out in [21] the flow field in the first case depends on three main parameters viz. Reynolds number, velocity ratio and time. The velocity ratio (a) is defined as a ¼ Uax , where x is the 1 angular velocity of the cylinder about its axis. We consider two

The flow due to bodies oscillating in a stationary or unsteady fluid and stationary bodies in an oscillating unsteady flow are considered to be fluid-structure interaction problem and is of immense importance in many fields of science and engineering. Some important experimental and numerical works on such problems can be found in [27–29]. Here we consider the case of an inline oscillating cylinder in a fluid at rest to validate the scheme’s ability in handling moving boundaries. In the present study the translational motion of the cylinder is prescribed by a simple harmonic oscillation

xðtÞ ¼ A cosð2pftÞ;

yðtÞ ¼ constant;

where A is the amplitude of the cylinder motion and f is the characteristic frequency of the oscillation. In this problem distances and velocities have been non-dimensionalized with respect to the diameter D and Umax respectively, where Umax is the maximum velocity of the cylinder motion. The two controlling parameters for this flow D are the Reynolds number Re ¼ Umax and the Keulegan–Carpenter m max number KC ¼ UfD ; m being the kinematic viscosity.

155

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

(a)

(b)

(c) Fig. 18. Comparison of computed stream lines (left) with experimental visualizations [38] (right) for NACA 0012 at Re = 1000, h = 34° (a) t = 1.6, (b) t = 4.0, (c) t = 5.6.

We simulate the flow in a non-inertial frame attached to the cylinder without changing the position of the artificial far-stream boundaries [28,30]. The relevant boundary conditions in non-inertial frame are as given below:

theirs. Note that corresponding to the parameter combination Re = 100, KC = 5 the flow is characterized by stable, symmetric and periodic vortex shedding and has been captured quite effectively in Fig. 12.

1. At the free-stream the time dependent velocity and stream function conditions are (u(t), v(t)) = (Umax sin(2pft), 0) and w = u(t)y respectively. 2. On the surface of the cylinder (u(t), v(t)) = (0, 0) and also w = 0.

4.5. Flow past an elliptic cylinder with angle of attack

The initial value of w has been set as zero. In the present work we have considered Re = 100 and KC = 65. Fig. 11 shows the computed velocity profiles (using lines) in the oscillation direction and transverse direction at four different x locations viz x = 0.6D, 0D, 0.6D and 1.2D each computed at three different phase angles 2pft = 90°, 120° and 240°. This is compared with the experimental results of Dütsch et al. [27] and a good match can be seen. Note that there is a phase difference between our phase angle and that of [27] because of the difference in the expression of x given above. In Fig. 12 we present computed stream function and vorticity contours for four different phase angles 6°, 102°, 198° and 270° for Re = 100, KC = 5. The parameters chosen corresponds to the work of Dütsch et al. [27] and our results match very well with

In this section, we further test our numerical scheme by conducting a series of simulations for the unsteady viscous flow over an inclined elliptic cylinder placed in a uniform stream. The purpose of this test is to investigate the performance of the scheme around bluff bodies with very high curvatures. The last century has seen some outstanding works [31–34] on this problem which continues to excite researchers in this century as well [18,35]. Following [35], we consider the elliptic transformation

x¼

cosh n cosðg þ hÞ ; 2 cosh n0

y¼

sinh n sinðg þ hÞ ; 2 cosh n0

where h, called the angle of attack, is the angle made by the uniform free stream with the positive direction of the major axis and n0 = tanh1(k), k being the ratio of minor and major axes of the cylinder, called the aspect ratio. Note that n = n0 defines the surface of the cylinder with chord length c = 1. The Reynolds number for the flow is defined as Re ¼ U1m c, where U1 is the free stream velocity and m is

156

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

1

1

0.5

0.5

0

0

-0.5

-0.5

-1

0

0.5

1

1.5

2

-1

0

0.5

1

(a) 1

0.5

0.5

0

0

-0.5

-0.5

0

0.5

1

1.5

2

-1

0

0.5

1

(c) 1

0.5

0.5

0

0

-0.5

-0.5

0

0.5

1

1.5

2

1.5

2

(d)

1

-1

2

(b)

1

-1

1.5

1.5

2

(e)

-1

0

0.5

1

(f)

Fig. 19. Evolution of post processed vorticity contours for NACA 0012 at Re = 1000, h = 34° (a) t = 1.6, (b) t = 2.8, (c) t = 3.2, (d) t = 4.0, (e) t = 4.8, (f) t = 5.6.

the kinematic viscosity. Apart from Re the other parameters controlling the flow are k and h. For this problem we consider free stream to be at a distance R1  20c and U1 = 1. The boundary and initial conditions are derived as in Section 4.3 for the flow past a rotating cylinder. The expressions used for drag coefficient (CD) and lift coefficient (CL) are as given in [35]. We start our computation with Re = 15 and 30 at h = 45°, k = 0.1. For such parameters the flow attains steady state [31,35]. We plot the evolution of drag and lift in the Fig. 13a as we time march towards the steady state; we compare CD and CL values obtained with the present scheme with those obtained by Lugt and Haussling [31] (by converting the values using the present definitions), and

Dennis and Young [35] in Table 5. The streamlines at steady state corresponding to both the above flow configurations (Fig. 14) match extremely well with the published results. In Table 6 we compare our computed drag and lift for Re = 200, h = 45° and k = 0.2 and k = 0.1 with those found in [31] for different times at its early stage of evolution and again we find a close match. We next consider the flow past elliptic cylinder with Re = 163, h = 45° and k = 0.2. Here flow field ultimately leads to periodic vortex shedding. This case has been considered by Peng et al. [18]. The time evolution of drag and lift coefficients towards an asymptotic periodic state can be seen in Fig. 13b. We present the temporal evolution of streamlines and post processed vorticity contours over

157

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

0.5

0

-0.5

-1

(a)

0

0.5

1

1.5

2

0

0.5

1

1.5

2

0

0.5

1

1.5

2

0

0.5

1

1.5

2

0.5

0

-0.5

-1

(b)

0.5

0

-0.5

-1

(c)

0.5

0

-0.5

-1

(d)

Fig. 20. Computed stream lines (left) and post processed vorticity contours (right) for NACA 0012 at Re = 3000, h = 45° (a) t = 1.5, (b) t = 2.5, (c) t = 3.5, (d) t = 4.5.

one complete vortex shedding cycle of duration T in Fig. 15. The evolution of von Kármán vortex street is clearly seen in these figures. Note that two eddies are shed just behind the cylinder within each period. Finally we consider the flow field with controlling parameters Re = 1000, h = 60°, k = 0.1. In Fig. 16i we present the early evolu-

tion of wake while the wake structure behind elliptic cylinder during shedding is shown in Fig. 16ii. In Fig. 16i one can see the formation of a main eddy along with a pair of secondary vortices mimicking the a-phenomenon together with the first vortex shedding from the trailing edge which was reported in [33] as well. This a-phenomenon is also seen to appear alternately

158

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

Fig. 21. Comparison of computed vorticity contours (left) with the visualization pictures (right) from Freymuth [36] for Re = 5200, h = 90° at time t = 0.5, 0.9, 1.9 and 2.2.

behind the tips of the elliptical cylinder for fully developed flow (Fig. 16ii). 4.6. Flow past symmetric aerofoils with angle of attack The last example is the determination of unsteady flow past symmetric aerofoils NACA 0012 and NACA 0015. The

literature contains some great classical [36–38] and modern [39,40] experimental results as well as some supporting computational works on incompressible flow past these aerofoils. The first challenge for us in this example is to model the NACA aerofoils and generate the surrounding grid using a conformal transformation. Here we use the conformal transformation

159

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

0

0

-0.5

-0.5

-1

-1

0

0.5

1

1.5

0

0.5

(a)

0

-0.5

-0.5

-1

-1

0.5

1

0

1.5

0.5

0

0

-0.5

-0.5

-1

-1

0.5

1

1.5

1

1.5

(d)

(c)

0

1.5

(b)

0

0

1

1

1.5

(e)

0

0.5

(f)

Fig. 22. Evolution of stream lines for NACA 0015 at Re = 5200, h = 90° (a) t = 0.1, (b) t = 0.5, (c) t = 0.9, (d) t = 1.2, (e) t = 1.9, (f) t = 2.2.

a a2 z ¼ k þ c1 þ c2 2 ; k k

k ¼ aeix ;

1 2

x ¼ n þ ig; a ¼  c1 :

For NACA 0012 we take c1 = 0.2238 and c2 = 0.0165; whereas for NACA 0015, c1 = 0.2165 and c2 = 0.0180. A comparison of the aerofoil NACA 0012 thus generated with the actual ones is presented in Fig. 17a. In Fig. 17b we depict the numerical difference between the actual NACA profiles and the model NACA profiles used in our computation. In the same figure we have also shown the error percentage of profiles generated by using Joukowski transformation. Note that Joukowski profiles have been and still continue to be useful models for the analysis of aerofoils [14]. It is worth mentioning that only at the edge, the percentage error is on a higher side although the absolute error is not very significant. This is primarily

due to narrow thickness at the edge leading to extremely small numerical values. Almost at every point along the length of both the aerofoils our models show much better approximation vis-avis Joukowski’s transformation. The computational grids for NACA 0015 at angle of attack 90° is shown in Fig. 17c. A zoomed view of the trailing edge has also been shown in the same figure from which the orthogonality of the grid lines in the physical plane is obvious. We take the chord length c = 1 with its tip at the origin. We consider two different combinations of flow control parameters viz. Reynolds number (Re) and angle of attack (h). For NACA 0012 aerofoils, they are (a) Re = 1000, h = 34°; (b) Re = 3000, h = 45°. In case of NACA 0015 aerofoils, we further consider two accelerated flow fields (c) Re = 5200, h = 90° and (d) Re = 8000, h = 30°. These choices were inspired by the existence of numerous

160

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

0.2

0

-0.2

-0.4 0

0.2

0.4

0.6

0.8

(a)

0.2

0

-0.2

-0.4 0

0.2

0.4

0.6

0.8

0

0.2

0.4

0.6

0.8

(b)

0.2

0

-0.2

-0.4

(c) Fig. 23. Comparison of computed nondimensional velocity and vorticity (left) with the experimental visualization [40] for Re = 8000, h = 30°: (a) t = 1.0, (b) t = 2.0, (c) t = 3.0.

numerical results and experimental visualizations in the literature. The Reynolds number is based on the chord length c of the aerofoil and the free-stream speed U1. For cases (a) and (b) we consider the flow to start impulsively with U1 = 1 and R1  12c, while for cases (c) and (d) we consider uniform acceleration where nondimensional velocity of unity is attained, starting with zero velocity, at times 0.0521 and 2.5 respectively. As an initial condition we take w = 0 = wn = wg. We also consider potential boundary condition at upstream and convective boundary condition for w at the downstream of the flow. In Fig. 18 we present comparisons of the streamlines between the computed results and those reported in [38] for Re = 1000, h = 34° showing excellent agreement. The time evolution of the flow structure is correctly reproduced and all the main and secondary vortices detected by experimental visualization are captured

by our numerical simulation. To illustrate the vorticity field, the post processed vorticity contours corresponding to different times are presented in Fig. 19. For Re = 3000, h = 45° the stream lines and the associated vortex structures can be seen in Fig. 20. As has been pointed out in [37], it is clear from the computed flow field that the initial wake follows an unsteady process consisting of the development of a large-scale leading-edge vortex over the upper surface and subsequent detachment of that vortex, followed by the establishment of reverse streamlines on the upper-surface and the subsequent development and shedding of the smaller trailing-edge vortex. The generation of leading-edge vortex with clockwise rotation and the trailing-edge vortex with counterclockwise rotation can also be seen in Fig. 20. Freymuth [36] experimentally studied the effect of uniform acceleration on flow past NACA 0015 aerofoil for different

161

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

0.5

0.5

0

0

-0.5

-0.5

0

0.5

1

0

(a) 0.5

0

0

-0.5

-0.5

0.5

0

1

0.5

0.5

0

0

-0.5

-0.5

0.5

0.5

1

(d)

(c)

0

1

(b)

0.5

0

0.5

1

(e)

0

0.5

1

(f)

Fig. 24. Evolution of stream lines for NACA 0015 at Re = 8000, h = 30° (a) t = 0.5, (b) t = 1.0, (c) t = 2.0, (d) t = 3.0, (e) t = 4.0, (f) t = 4.5.

Reynolds numbers (Re) and angles of attack (h). It is well known that for a fixed Re, as h increases the shedding of vortices from the leading and trailing edges gets more pronounced. Note that in the above mentioned work, the Reynolds number has been defined in terms of the uniform acceleration and the chord length of the aerofoil. As pointed out by Sengupta et al. [39] the Reynolds number based on instantaneous speed, as defined here will vary linearly with time. We follow the convention adopted in [39] and choose the combination Re = 5200 and h = 90° to carry out computations in order to compare our results with the strong vortical structures that were seen in the experiments of Freymuth [36]. In Fig. 21 the vorticity contours are compared with the flow visualization experiment performed in [36]. The match between

the computation and experiment is striking for the vortices that emanate from both leading and trailing edges. Our computation clearly captures the instability of the vortices shed from the top and bottom resulting in a train of smaller vortices. In Fig. 22 we present the evolution of streamlines, where as time progresses, the growth of two dissimilar bubbles, with the larger one on the trailing edge can be observed. Another case is inspired by the works of Soria et al. [40] where the NACA 0015 aerofoil was held at 30° angle of attack. The flow here is accelerated uniformly from rest upto a desired free stream speed and once the free stream speed is attained it has been maintained thereafter. In Fig. 23, we compare our results with the first of the two experimental results reported in

162

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

[40] for Re = 8000. It shows the temporal evolution of velocity and vorticity as the flow accelerates and shortly afterwards. From Fig. 23a it is clear that the leading edge separation starts after t = 1.0 which is very accurately captured by our numerical simulation in Fig. 23b at time t = 2.0 as reported in [40]. Fig. 23b also confirms that the boundary layer remains attached to the downstream side of the aerofoil. As the acceleration phase gets over at time t = 2.5, the velocity stays constant and flow becomes exceedingly complex; this can be seen in Fig. 23c which shows the development of strong shear layer at the leading edge at t = 3.0. Note that our computed figures on the left of Fig. 23 match very well with the experimental visualization of [40] on the right. The evolution of streamlines for the NACA 0015 aerofoil at h = 30° for Re = 8000 can be seen in Fig. 24. Because of the insufficiency of grid resolution some fine structures that were reported in the experiments have not been found in the numerical simulation. Nevertheless our scheme accurately captures the physical phenomena.

5. Conclusion Numerical schemes based on pure streamfunction formulation of the transient Navier–Stokes equations have so far been developed only on cartesian coordinates for rectangular domains. Here we develop an implicit, temporally and spatially second order accurate compact finite difference scheme for such a formulation that allows us to effectively simulate unsteady viscous incompressible flows in irregular geometries that can be conformably mapped onto rectangular domains. This formulation has been found to be very effective in capturing flow physics within curved regions as also the dynamics of a fluid-embedded body interaction. Spectral analysis shows that the scheme produces accurate solutions in terms of phase and amplitude errors. We perform numerical experiments on six different problems of varying physical complexities in order to extensively investigate the efficiency of the present scheme. In the process, we also carry out a stability and convergence analysis; the first problem considered is that of decaying vortices with known closed formed solution, where the rate of convergence of the present method is established. The problem of flow through constricted channel, where we time march to steady state, is an internal flow problem and serves as a suitable test case for convergence analysis of problems having no analytical solution. The next series of numerical test cases represent flow over bluff bodies of different curvatures: the flow past circular cylinders having constant curvature, elliptical cylinders having varied curvature and aerofoils having infinite curvature at the trailing edge. For the circular cylinder, two different physical scenarios are considered. The cases of rotating cylinder with constant and oscillatory rotations are considered along with the case of inline oscillatory cylinder is considered. All these rotatory and oscillatory cases symbolize moving boundaries. The robustness of our scheme is exemplified by the simulation of flow past elliptical cylinders and symmetric aerofoils with different angles of attack; more so for the aerofoils where both uniform and accelerated flows are precisely simulated for a wide range of Reynolds numbers. In all cases, our results are validated by comparison with existing qualitative and quantitative results, both experimental and numerical; besides excellent comparison with the experimental results, our results are closer to the experimental ones than some previously published numerical results. This shows clearly that our scheme is successful in capturing two-dimensional incompressible viscous flows in non-rectangular physical domains.

6. Appendix The N–S equations for unsteady 2D incompressible viscous flows together with the equation of continuity in non dimensional form can be written as:

ux þ v y ¼ 0; @ 1 u þ uux þ v uy ¼ px þ ðuxx þ uyy Þ; @t Re @ 1 v þ uv x þ vv y ¼ py þ ðv xx þ v yy Þ: @t Re

ð20Þ ð21Þ ð22Þ

Here Re ¼ UL m is the Reynolds number based on characteristic length L and characteristic velocity U of the flow. Introducing streamfunction w and vorticity x we have u ¼ @w ; v ¼  @w ; x ¼ @@xv  @u and @y @x @y the above formulation may be rewritten as:

wxx þ wyy ¼ x;

ð23Þ

@ 1 x ¼ ðxxx þ xyy Þ  ðuxx þ v xy Þ: @t Re

ð24Þ

Let the physical (x, y) plane be transformed into a computational (n, g) plane using the mapping:

x ¼ xðn; gÞ;

y ¼ yðn; gÞ:

ð25Þ

Under this transformation Eqs. (23) and (24) in the computational plane becomes:

a1 wnn þ e1 wng þ b1 wgg þ c1 wn þ d1 wg ¼ f1 ;

ð26Þ

@ x ¼ a2 xnn þ e2 xng þ b2 xgg þ c2 xn þ d2 xg ; @t

ð27Þ

where

a1 ¼

c1 ¼

d1 ¼

1 J

2

b1 ¼

1

J2

 x2n þ y2n ;

1h

 

  yg x2g þ y2g xnn  2ðyg yn þ xg xn Þxng þ x2n þ y2n xgg J  

 i þ xg x2g þ y2g ynn  2ðyg yn þ xg xn Þyng þ x2n þ y2n ygg ; 3



  1 h  2 yn xg þ y2g xnn  2ðyg yn þ xg xn Þxng þ x2n þ y2n xgg 3 J  

   leftxn x2g þ y2g ynn  2ðyg yn þ xg xn Þyng þ x2n þ y2n ygg 

e1 ¼ 

a2 ¼

 x2g þ y2g ;

2 J2

ðyg yn þ xg xn Þ;

1  ReJ

2

 x2g þ y2g ;

f 1 ¼ x;

b2 ¼

1

ReJ2

 x2n þ y2n ;

  1 1 h c2 ¼ ðuyg  v xg Þ þ yg x2g þ y2g xnn 3 J ReJ

   2ðyg yn þ xg xn Þxng þ x2n þ y2n xgg  

 i þ xg x2g þ y2g ynn  2ðyg yn þ xg xn Þyng þ x2n þ y2n ygg ;  1 1 h  2 2 xnn  2ðyg yn þ xg xn Þxng d2 ¼ ðuyn þ v xn Þ þ y x þ y n g g J ReJ3   

  þ x2n þ y2n xgg  xn x2g þ y2g ynn  2ðyg yn þ xg xn Þyng

 i þ x2n þ y2n ygg ; e2 ¼ 

2 ReJ 2

ðyg yn þ xg xn Þ;

S. Sen et al. / Computers & Fluids 84 (2013) 141–163

where J = xnyg  ynxg is the jacobian of the transformation with u and v as given in (5) and (6) respectively. If the transformation (25) is a conformal one i.e. of the form:

z ¼ zðhÞ;

ð28Þ

z = x + iy and h = n + ig; then Eqs. (26) and (27) with the help of expressions (5) and (6) simplify to:

wnn þ wgg ¼ J x;

ð29Þ

@ 1 1 x¼ ðxnn þ xgg Þ  ðxn wg  xg wn Þ: @t ReJ J

ð30Þ

Note that in this case the jacobian of the transformation is, J ¼ xn yg  yn xg ¼ x2n þ y2n ¼ x2g þ y2g . Eliminating x from Eqs. (29) and (30), we obtain the following form of the N–S equations:

ReJ

@ ðw þ wgg Þ ¼ðwnnnn þ 2wnngg þ wgggg Þ  ð2C þ Rewg Þ @t nn  ðwnnn þ wngg Þ  ð2D  Rewn Þðwnng þ wggg Þ þ ðE þ CRewg  DRewn Þðwnn þ wgg Þ

ð31Þ

where C = Jn/J, D = Jg/J, and E = 2C2 + 2D2  Jgg/J  Jnn/J. Eq. (31) is of form:

@ 2 r w ¼ aðn; gÞr4 w þ bðn; g; wn ; wg Þr2 wn @t þ cðn; g; wn ; wg Þr2 wg þ dðn; g; wn ; wg Þr2 w:

ð32Þ

References [1] Gresho PM. Incompressible fluid dynamics: some fundamental formulation issues. Annu Rev Fluid Mech 1991;23:413–53. [2] Chung TJ. Computational fluid dynamics. Cambridge: CUP; 2002. [3] Goodrich JW, Sox WY. Time-dependent viscous incompressible Navier–Stokes equations: The finite difference Galerkin formulation and streamfunction algorithms. J Comput Phys 1989;84:207–41. [4] Hirsh RS. Higher order accurate difference solutions of fluid mechanics problems by a compact differencing technique. J Comput Phys 1975;19:90–109. [5] Krause E, Hirschel EH, Kordulla W. Fourth order Mehrstellen integration for 3D turbulent boundary layer. Comput Fluids 1976;4:77–92. [6] Adam Y. Highly accurate compact implicit methods and boundary conditions. J Comput Phys 1977;24:10–22. [7] Lecointe Y, Piquet J. On the use of several compact methods for the study of unsteady incompressible viscous flow round a circular cylinder. Comput Fluids 1984;12:255–80. [8] Schiestel R, Viazzo S. A Hermitian-Fourier numerical method for solving the incompressible Navier–Stokes equations. Comput Fluids 1995;24:739–52. [9] Kupferman R. A central difference scheme for a pure streamfunction formulation of incompressible viscous flow. SIAM J Sci Comput 2001;23:1–18. [10] Fishelov D, Ben-Artzi M, Croisille JP. A compact scheme for the streamfunction formulation of Navier–Stokes equations. Notes Comput Sci 2003(266):809–17. Springer–Verlag. [11] Ben-Artzi M, Croisille JP, Fishelov D, Trachtenberg S. A pure-compact scheme for the streamfunction formulation of Navier–Stokes equations.. J Comput Phys 2005;205:640–64.

163

[12] Gupta MM, Kalita JC. A new paradigm for solving NavierStokes equations: streamfunction–velocity formulation. J Comput Phys 2005;207:52–68. [13] Kalita JC, Gupta MM. A streamfunction–velocity approach for the 2D transient incompressible viscous flows. Int J Numer Method Fluids 2010;62:237–66. [14] Papamichael N, Stylianopoulos N. Numerical conformal mapping domain decomposition and the mapping of quadrilaterals. World Scientific; 2010. [15] Gupta MM, Manohar RP. Direct solution of the biharmonic equation using noncoupled approach. J Comput Phys 1979;33:236–48. [16] Stephenson JW. Single cell discretization of order two and four for biharmonic problems. J Comput Phys 1984;55:65–80. [17] Altas I, Dym J, Gupta MM, Manohar RP. Multigrid solution of automatically generated high-order discretizations for the biharmonic equation. SIAM J Sci Comput 1998;19:1575–85. [18] Peng YF, Mittal R, Sau A, Hwang RR. Nested Cartesian grid method in incompressible viscous fluid flow. J Comput Phys 2010;229:7072–101. [19] Pandit SK, Kalita JC, Dalal DC. A transient higher order compact scheme for incompressible viscous flows on geometries beyond rectangular. J Comput Phys 2007;225:1100–24. [20] Mancera PFA. A study of numerical solution of the steady two dimensional Navier–Stokes equations in a constricted channel problem by a compact fourth order method. Appl Math Comput 2003;146:771–90. [21] Badr HM, Coutanceau M, Dennis SCR, Menard C. Unsteady flow past a rotating circular cylinder at Reynolds numbers 103 and 104. J Fluid Mech 1990;220:459–84. [22] Chew YT, Cheng M, Luo SC. A numerical study of flow past a rotating circular cylinder using a hybrid vortex scheme. J Fluid Mech 1995;299:35–71. [23] Thiria B, Durand SG, Wesfreid JE. The wake of a cylinder performing rotary oscillations. J Fluid Mech 2006;560:123–47. [24] Sengupta TK, Bhumkar YG. New explicit two-dimensional higher order filters. Comput Fluids 2010;39:1848–63. [25] Kalita JC, Sen S. Triggering asymmetry for flow past circular cylinder at low Reynolds numbers. Comput Fluids 2012;59:44–60. [26] Protas B, Wesfreid JE. Drag force in the open-loop control of the cylinder wake in the laminar regime. Phys Fluids 2002;14:810–26. [27] Dütsch H, Durst F, Becker S, Lienhart H. Low-Reynolds-number flow around an oscillating circular cylinder at low Keulegan–Carpenter numbers. J Fluid Mech 1998;360:249–71. [28] Blackburn HM, Henderson RD. A study of two-dimensional flow past an oscillating cylinder. J Fluid Mech 1999;385:255–85. [29] Choi JI, Oberoi RC, Edwards JR, Rosati JA. An immersed boundary method for complex incompressible flows. J Comput Phys 2007;224:757–84. [30] Pope SB. Turbulent flows. Cambridge: CUP; 2000. [31] Lugt HJ, Haussling HJ. Laminar flow past an abruptly accelerated elliptic cylinder at 45° incidence. J Fluid Mech 1974;65:711–34. [32] Feng J, Hu HH, Joseph DD. Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid. Part 1: Sedimentation. J Fluid Mech 1994;261:95–134. [33] Chou MH, Huang W. Numerical study of high-Reynolds-number flow past a bluff object. Int J Numer Method Fluids 1996;23:711–32. [34] Nair MT, Sengupta TK. Unsteady flow past elliptic cylinders. J Fluid Struct 1997;11:555–95. [35] Dennis SCR, Young PJS. Steady flow past an elliptic cylinder inclined to the stream. J Eng Math 2003;47:101–20. [36] Freymuth P. The vortex patterns of dynamic separation: a parametric and comparative study. Prog Aerospace Sci 1985;22:161–208. [37] Ohmi K, Coutanceau M, Daube O, Loc TP. Further experiments on vortex formation around an oscillating and translating airfoil at large incidences. J Fluid Mech 1991;225:607–30. [38] Shen WZ, Loc TP. A coupling finite difference/particle method for the resolution of 2D Navier–Stokes equations in velocity–vorticity form. Aerospace Sci Technol 1997;2:97–109. [39] Sengupta TK, Lim TT, Sajjan SV, Ganesh S, Soria J. Accelerated flow past a symmetric aerofoil: experiments and computations. J Fluid Mech 2007;591:255–88. [40] Soria J, New TH, Lim TT, Parker K. Multigrid CCDPIV measurements of accelerated flow past an airfoil at an angle of attack of 30°. Exp Thermal Fluid Sci 2003;27:667–76.