Incompressible flow computations over moving boundary using a novel upwind method

Computers & Fluids 46 (2011) 348–352

Contents lists available at 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

Incompressible flow computations over moving boundary using a novel upwind method J.C. Mandal ⇑, C.R. Sonawane, A.S. Iyer, S.J. GosaviInamdar Department of Aerospace Engineering, Indian Institute of Technology Bombay, Mumbai 400076, India

a r t i c l e

i n f o

Article history: Available online 22 August 2010 Keywords: Incompressible flow HLLC–AC Arbitrary Lagrangian Eulerian Dual time stepping Radial basis function Moving indentation

a b s t r a c t In this paper a novel method for simulating unsteady incompressible viscous flow over a moving boundary is described. The numerical model is based on a 2D Navier–Stokes incompressible flow in artificial compressibility formulation with Arbitrary Lagrangian Eulerian approach for moving grid and dual time stepping approach for time accurate discretization. A higher order unstructured finite volume scheme, based on a Harten Lax and van Leer with Contact (HLLC) type Riemann solver for convective fluxes, developed for steady incompressible flow in artificial compressibility formulation by Mandal and Iyer (AIAA paper 2009-3541), is extended to solve unsteady flows over moving boundary. Viscous fluxes are discretized in a central differencing manner based on Coirier’s diamond path. An algorithm based on interpolation with radial basis functions is used for grid movements. The present numerical scheme is validated for an unsteady channel flow with a moving indentation. The present numerical results are found to agree well with experimental results reported in literature. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Incompressible flows over moving boundaries are encountered in many practical situations. The main feature of these flows is their unsteadiness with respect to shape of the boundaries and corresponding flow patterns. One such example is the flow over a moving indentation inside a channel [1], which represents flow features of oscillating stenosis of a blood vessel. The flow inside the channel with moving boundary results in transient and complex flow phenomena mainly due to the interaction between the moving boundary and the flowing fluid. In this study, the flow patterns are computed by solving unsteady incompressible Navier– Stokes (NS) equations in artificial compressibility (AC) formulation [2] with Arbitrary Lagrangian Eulerian (ALE) [3] and dual time stepping (DTS) [4] approach using a novel upwind method. One of the most common techniques to solve incompressible Navier–Stokes equations is the AC formulation originally introduced by Chorin [2] for steady state solutions. In recent times, this formulation has drawn considerable attention due to its many desirable properties, such as the possibility of solving the system in a closely coupled manner, superior convergence [5], automatic satisfaction of divergence free condition to the level of residual error at steady state and its ability to adapt advanced numerical techniques developed for compressible flow.

⇑ Corresponding author. Tel.: +91 22 25767129; fax: +91 22 25722602. E-mail address: [email protected] (J.C. Mandal). 0045-7930/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2010.08.011

Mostly two approaches, namely the Jameson method and the Roe method are used for finite volume discretization of the convective fluxes in the AC formulation. The Jameson finite volume method uses central differencing for the convective fluxes, where dissipative terms are needed to stabilize the discretized system. In the Roe’s method, the convective part is linearized and a Riemann problem is solved at a cell interface to evaluate the numerical fluxes. There is yet another method to solve steady incompressible flows in AC formulation, recently introduced by Mandal and Iyer [6], where the non-linear convective part is discretized following an approach similar to the Harten Lax and van Leer with Contact (HLLC) [7] method for compressible flow equations. In the above HLLC type upwind method for AC formulation, henceforth referred to as HLLC–AC method, a three wave structure is considered for solving a Riemann problem arising in the finite volume discretization. In this paper, an extension of the HLLC–AC method is considered in order to incorporate dual time stepping along with the Arbitrary Lagrangian Eulerian (ALE) formulation to compute unsteady flows over moving boundary involving dynamic meshes. The Thin Plate Spline (TPS) based radial basis function (RBF) interpolation method [8] is used for the mesh movements, which is found to maintain better grid quality as compared to linear spring analogy method. The paper is organized in the following manner. After a brief introduction to governing equations in integral form in Section 2, the finite volume discretization is explained in Section 3. Evaluation of ALE flux and application of Geometric Conservation Law (GCL) are also explained in Section 3. Brief description about

349

J.C. Mandal et al. / Computers & Fluids 46 (2011) 348–352

inviscid flux evaluation using HLLC–AC [6] method, and viscous flux computation are given in Section 4. Grid movement strategy by RBF based interpolation is presented briefly in Section 5. Description of the moving indentation channel problem [1] and its solution with the present method are provided in Section 6 in order to demonstrate the potential of the present approach.

k

n i

2. Governing equations

L R

l

j

The integral form of 2D unsteady incompressible Navier–Stokes equations in AC formulation, where DTS approach is considered for time accurate solution and ALE formulation for mesh movements, can be written as [14]

Z Z

@W dX þIM X @s

Z Z

@W dX þ HM X @t

I

q

½ðEc Ev Þlx þðGc Gv Þly dA ¼ 0

A

ð1Þ T

where W = (p, u, v) is the field vector; E = (U, uU + p, vU)T, Gc = (V, uV, vV + p)T are the convective fluxes; Ev = (0, rxx, rxy)T, Gv = (0, ryx, ryy)T are the viscous fluxes; lx, ly are direction cosines for the face area A, which is enclosing the cell volume X. Here, the convective velocity components in the referential frame are U = u  xt and V = v  yt with xt and yt as the velocities of the grid in x and y directions respectively. The variables t and s are real time and pseudo-time respectively. A steady state solution in pseudotime corresponds to an instantaneous unsteady solution in real time. The expressions for the stress components used here are

   1 @u 2 @u @ v 2  þ rxx ¼ Re @x 3 @x @y    1 @ v 2 @u @ v 2 þ ryy ¼  Re @y 3 @x @y   1 @u @ v rxy ¼ ryx ¼ þ Re @y @x 0 2 1 0 0 0 b B B C M I ¼ @0 1 0A H ¼ @ 0 0 0 1 0 0

0

0

1

C 1 0A

In the finite volume method, a computational domain is generally divided into structured/unstructured grid consisting of quadrilateral/triangular cells in 2D. Approximating the @/@t operator by a second order implicit backward difference formula and splitting the convective fluxes (Ec and Gc) into stationary reference fluxes denoted by subscript s and ALE fluxes denoted by subscript a, the finite volume discretization of the Eq. (1) can be written for cell i (see Fig. 1) as

X

" # nþ1 Wnþ1  4Xni Wni þ Xn1 Wn1 @Wi M 3Xi i i i þI @s 2 Dt

þ HM

M X 

ðEcs  Ev Þlx þ ðGcs  Gv Þly

nþ1 m

DAm

M X  c nþ1 Ea lx þ Gca ly m DAm ¼ 0

X

dxdy þ

I A

½xt dy  yt dx ¼ 0

which for the cell i becomes M X 3Xnþ1  4Xni þ Xn1 i i ½xt Dy  yt Dxnþ1 ¼0  m 2 Dt m¼1

ð4Þ

Thus, the moving mesh effect can be taken care by subtracting the volumetric increment DXm along the face m from continuity equation and subtracting the velocity component multiplied by the volumetric increment from the respective momentum equations. It is clear from Eq. (2) that the numerical discretization of the stationary reference convective fluxes can be treated in a conventional manner as in steady problem. At the mth face, the normal components of the stationary convective flux F ij ¼ ðEcs lx þ Gcs ly Þ at cell interface between two cells i and j is evaluated based on a new upwind method called HLLC–AC [6], which is described next. The ordinary differential equation in pseudo-time (3) is integrated using a five stage Runge–Kutta scheme [10]. Note that the required real time accurate solution at time level (n + 1) satisfies R ðWnþ1 Þ ¼ 0 and this is found i by marching Eq. (3) to a steady state in pseudo-time. 4. Computation of convective and viscous fluxes

m¼1

 HM

Z Z

½Eca lx þ Gca ly m DAm ¼ ðDX; uDX; v DXÞTm

0 1

3. Finite volume discretization

nþ1 i

@ @t

The second term in Eq. (4) represents the rate of face swept area, where Dx and Dy are the projections of a face on x and y co-ordinates respectively. In Eqs. (2) and (4) the summations are over all the cell faces corresponding to the index m = 1, 2, . . ., M that constitute the boundary of the control volume Xi. In Eq. (2), Wi is the finite volume average value of the field vector at cell center; Ecs ¼ ðu; u2 þ p; uv ÞT and Gcs ¼ ðv ; uv ; v 2 þ pÞT are the stationary reference convective fluxes; Eca ¼ ðxt ; uxt ; v xt ÞT ; Gca ¼ ðyt ; uyt ; v yt ÞT are the ALE convective fluxes. It can be shown, the effect of ALE flux is the volumetric increment along the face m and is given as

The matrices IM, HM in (1) are given as M

Fig. 1. Sketch of a typical finite volume grid.

c

ð2Þ

4.1. Stationary reference convective flux computation

m¼1

or

Xnþ1 i

@Wi þ R ðWnþ1 Þ¼0 i @s

ð3Þ

where Xnþ1 is obtained by solving the discretized form of Geometric i Conservation Law (GCL) [9]

The normal components of the convective flux vector F ij ¼ ðEcs lx þ Gcs ly Þ at the interface between two cells i and j (as seen in Fig. 1) are evaluated based on a new upwind method for pseudocompressibility formulation, HLLC–AC [6]. The interface stationary reference convective flux is evaluated based on a Riemann solver using three wave structure (Fig. 2) as

350

J.C. Mandal et al. / Computers & Fluids 46 (2011) 348–352

t

* S

SL v*n L

v*n R

4.1.1. Interface values of primitive variables (WL and WR) The left and right state of primitive variables are evaluated at the interface between the cells i and j from cell average values using a piecewise linear reconstruction of Barth and Jespersen [11] as

S R

WL ¼ Wi þ wi ðrWi  ~ r L Þ;

p* u*n p W = un L vn L

WR ¼ Wj þ wj ðrWj  ~ rR Þ

where rWi is the gradient of W at the cell center i. The function w is Barth limiter [11] and ~ r is the displacement vector from cell center i to the interface. The gradients required in the above equation is evaluated using a least squares technique proposed by Barth and Jespersen [11].

p W = un R v n R n

4.2. Viscous flux computation Fig. 2. Three wave system for the present formulation.

8 FL > > > < F ¼ F þ S ðW  W Þ L L L L L Fi;j ¼   > F ¼ F þ S ðW  W R R RÞ > R R > : FR

0 6 SL SL 6 0 6 S S 6 0 6 SR

ð5Þ

0 P SR

where SL, SR and S* are the estimated wave speeds. The vectors, FL and FR, are convective fluxes on the left and right of the interface between cells i and j. The primitive variable vectors, WL and WR, are evaluated at left and right states using a piecewise linear reconstruction presented by Barth and Jespersen [11]. Realizing that the flux Jacobian for the above system has eigenvalues un, un  a, and un + a, the fastest left and right running waves and the middle wave are estimated as

SL ¼ min½ðun  aÞL ; ðun  aÞR  SR ¼ max½ðun þ aÞL ; ðun þ aÞR  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S ¼ un with a ¼ u2n þ b2 Note that the subscripts L and R denote left and right states at the interface respectively. The unknown vectors defining the star quantities

WL ¼ ½p un

v n TL

and

WR ¼ ½p un

v n TR

in Eq. (5) are evaluated in the following manner. As shown in Fig. 2, there exists jumps in all primitive variables across the fastest left and right running waves (SL, SR). Based on generalized Riemann Invariant analysis similar to that in Ref. [7], it is found that there is a jump in the face tangential velocity across the intermediate wave with wave speed S*, while the pressure and the normal component of velocity remain invariant across the middle wave. Therefore, the intermediate (star) quantities for pressure and normal velocity can be obtained from the two wave system as

SR pR  SL pL þ b2 ðun ÞL  b2 ðun ÞR SR  SL SR ðun ÞR  SL ðun ÞL þ ðu2n ÞL þ pL  ðu2n ÞR  pR  un ¼ SR  SL p ¼

To evaluate the tangential velocity components in the intermediate (star) region, i.e. v L ; v R , the Rankine Hugnoit conditions are solved across the left and right running waves leading to

v nL ¼

v nL ½SL  ðun ÞL  SL  ðun ÞL

v nR ¼

Discretization of viscous fluxes is done in a central differencing manner, where the calculation of spatial derivatives of velocity components at cell interfaces are performed by applying the Green-Gauss divergence theorem along the Courier diamond path [12].

v nR ½SR  ðun ÞR  SR  ðun ÞR

The interface stationary reference convective flux is thus evaluated for each edge of a cell using Eq. (5) once WL and WR are computed as above.

5. Grid movement using radial basis function interpolation Given displacement of boundary nodes dðxb j Þ at discrete grid points j = 1, . . ., nb, the displacements of the interior nodes are derived using radial basis function interpolation [8]. The interpolation function ‘‘s” describing the displacement in the whole domain, is approximated by a sum of radial basis functions

sðxÞ ¼

nb X

aj /ðkx  xbj kÞ þ pðxÞ

ð6Þ

j¼1

where xbj = [xbj,ybj] are the co-ordinates of boundary nodes with known displacement values. Here, p(x) = c1 + c2x + c3y is polynomial, nb is number of boundary nodes and / is a given Thin Plate Spline (TPS) radial basis function with respect to the Euclidean distance kxk and given as

/ðkx  xbj kÞ ¼ ðkx  xbj kÞ2 lnkðx  xbj Þk The coefficients aj and c1, c2, c3 from the polynomial p(x) are found using the interpolation conditions

sðxbj Þ ¼ dðxbj Þ

ð7Þ

and an additional requirements nb X

aj qðxbj Þ ¼ 0

ð8Þ

j¼1

for all first degree polynomials q. The above Eqs. (7) and (8) lead to system of equations for the coefficient vectors a = [a1,   , anb]T and c = [c1, c2, c3]T of the form

"

M b;b

Pb

P Tb

0

# 

 

a d ¼ c 0

ð9Þ

where, Mb,b is nb  nb matrix containing evaluation of the basis function /bi ;bj ¼ /ðkxbi  xbj kÞ; Pb is nb  3 matrix with row j given by [1, xbj, ybj] and d ¼ ½dðxb1 Þ;    ; dðxbn b ÞT . The values for the displacement in the interior grid points d(xin), can now be derived by evaluating the interpolation function (6) as d(xin) = s(xin). 6. Numerical results The present numerical scheme, which is an extension of HLLC–AC to include DTS and ALE formulation, is applied to solve

351

J.C. Mandal et al. / Computers & Fluids 46 (2011) 348–352

0.3

0.25

0.15

0.1

0.05

0

f ðx; tÞ ¼ gðxÞ  hðtÞ

-10

-5

0

5

10

X coordinate

where h(t) is given as,

1 ½1  cos2pt=T; t=T 2 ½0; 1 28 0 for x1 > > > > 1 >  ½1 þ tanhb ðx  x Þ for x2 > a dyn <2 gðxÞ ¼  for x3 > > > for x5 > 12 ½1 þ tanhbdyn  > > :  for x5

dp calculated dp Theory

0.2

dp

a moving indentation in a channel flow. Numerical results are validated with experimental results of Pedley and Stephanoff [1]. The present results are also compared with computational results of Ralph and Padley [13] and with Zhao and Forhad [14], who use a pseudo-compressibility formulation similar to the present one. The moving indentation is an interesting problem, as the literature results are not consistent and the flow features reported vary across results produced by different researchers. The geometry and dynamics of the moving indentation channel chosen here is identical to the experimental work of Pedley and Stephanoff [1]. The aspect ratio of the channel is 10, and features of the flow are observed to be uniform in the span-wise direction over most of the flow cycle, so that a two-dimensional simulation is appropriate. The schematic of flow problem is shown in Fig. 3. The shape of indentation varies as a function of space and time f(x, t) as,

Fig. 4. Comparison of analytical and computational streamwise pressure variations.

hðtÞ ¼

9 < x 6 x2 > > > > > < x 6 x3 > = < x 6 x4 > > < x 6 x5 > > > > ; < x 6 x6

The parameters, x1 = 13.75, x2 = 11.75, x3 = 9.25, x4 = 1.25, x5 = 1.25, x6 = 13.25,  = 0.38, and bdyn = 4.14, similar to those of Zhao and Forhad [14] have been chosen. No slip boundary condition is applied on the walls, which is specified by imposing (u, v)wall = (0, 0) on static walls, and (u, v)wall = (xt, yt)wall on moving wall. At inlet, a Poiseuille velocity profile is specified as boundary condition. At the outlet, back pressure pb is specified and velocities are extrapolated. The inlet pressure is computed from the interior. The computed streamwise pressure variation in the undisturbed channel is compared with the analytical solution given by expression (10).

dp 2umax l ¼ 2 dx b

ð10Þ

where, umax is the maximum axial velocity along the centerline, l is dynamic coefficient of viscosity, b is half width of channel. The pressure is always normalised by density. Fig. 4 shows comparison of streamwise computed pressure variation (dp = p  pb) with analytical solution for the steady state (corresponding to t = 0), where pb is back pressure, which is set to be 1 atm. For the given shape and amplitude of the wall indentation, the problem is specified by two non-dimensional parameters, Reynolds number Re = bUref/m and Strouhal number St = b/UrefT. Here, b is the undisturbed channel width, Uref is the average velocity of the flow upstream, m is the kinematic viscosity of the fluid and T is the oscillation period. The present calculation are carried out for Re = 507 and St = 0.037. Although quadrilateral grids are used,

Fig. 3. Schematic of the problem setup.

unstructured data structure is used for flow computations. Based on grid independence tests a mesh size with 61,812 nodes is chosen for the flow computation. The movement of grids is achieved using radial basis function based interpolation explained in Section 5. As the indentation starts moving from flat position, a separation vortex ‘A’ generates on its downstream side (see Fig. 5). With the passage of time, a second vortex ‘B’ of opposite sign is found to form on opposite wall, and then a third vortex ‘C’ further downstream of the indentation is generated and convected downstream. As the time increases, the number of vortices and the amplitude of the core waviness increases. Fig. 5 shows the instantaneous streamlines showing the generation of vortices produced by present method and that by other researchers reported in literature [13,14]. The present computations are able to capture the desired flow features very well. A quantitative comparison between the present results and literature results [1,13,14] in terms of time evolution of shed vortices ‘B’, ‘C’, ‘D’ are presented graphically in Fig. 6, which are found to be in good agreement. The discrepancies seen here are within the range of experimental scatter. The discrepancies are also due to expected lower accuracy associated with finding locations of vortices, which are based on derived quantity (vorticity) involving differentiation. 7. Conclusion A dual time stepping method for incompressible flow computation on moving grid using pseudo-compressibility approach with Arbitrary Lagrangian Eulerian formulation is presented. The nonlinear convective fluxes are discretized with a novel upwind method HLLC–AC scheme. Neither external numerical dissipation for stabilization of the numerical scheme, nor linearization of the convective terms are required in the present approach. Since, the eigenvalues of Jacobian matrices in pseudo-compressibility formulation is similar to subsonic like case of compressible flow, no discontinuous solution is expected in the flow near pseudo-steady state. Consequently, as the solution approaches a pseudo-steady state, the wave speed can be estimated very accurately purely based on the eigenvalues. Therefore, the new upwind method is expected to be an accurate modeling for pseudo-compressibility formulation. Here, the viscous fluxes are discretized in a conventional central differencing manner. The grid movements are

352

J.C. Mandal et al. / Computers & Fluids 46 (2011) 348–352

Location of shed vortices

Fig. 5. Instantaneous streamlines for moving indentation; (I) present, (II) Ralph and Pedley [13], and (III) Zhao and Forhad [14].

10

D Pedley & stephanhoff Expt Present Computations Ralph & Pedley computation Zhao et al Computation

8

C

6 B

4 2 0.4

0.5

0.6

0.7

0.8

Time (t/T) Fig. 6. Time evolutions of shed vortices ‘B’,‘C’, and ‘D’.

achieved using RBF interpolations, which are found to maintain good quality of the grid. The present method has been applied to solve a 2D channel with a moving indentation for validation. The results obtained from the present method show very good agreement with the experimental and the numerical results reported in literature.

Acknowledgement This work has been partially supported by Aeronautical Research and Development Board, India.

References [1] Pedley TJ, Stephanoff KD. Flow along a channel with a time dependent indentation in one wall: the generation of vorticity waves. J Fluid Mech 1985;160:337–67. [2] Chorin AJ. A numerical method for solving incompressible viscous flow problems. J Comput Phys 1967;2:12–26. [3] Donea J, Huerta A, Ponthot JP, Ferran AR. Arbitrary lagrangian-eulerian methods. John Wiley & Sons Ltd.; 2004. p. 1–25. [chapter 14]. [4] Jameson A. Time dependent Calculations Using Multigrid, with applications to unsteady flows past Airfoils and wings. AIAA paper June 1991; 91:1596. [5] Tamamidis P, Zhang G, Assanis DN. Comparison of pressure based and artificial compressibility methods for solving 3D steady incompressible viscous flows. J Comput Phys 1996;124:1–13. [6] Mandal JC, Iyer AS. An upwind method for incompressible flow computations using pseudo-compressibility approach. AIAA paper no. AIAA-2009-3541, In: 19th computational fluid dynamics conference, Jun 2009, San Antonio, USA; 2009. [7] Toro EF, Spruce M, Speares W. Restoration of the contact surface in the HLL Riemann solver. Shock Waves 1994;4(25):25–34. [8] de Boer A, van der Schoot MS, Bijl H. Mesh deformation based on radial basis function interpolation. Comput Struct 2007;85:784–95. [9] Thomas PD, Lombard CK. Geometric conservation law and its application to flow computations on moving grids. AIAA J 1979;17:1030–7. [10] Jameson A, Schmidta W, Turkel E. Numerical solution of Euler equations by finite volume methods using Runge–Kutta time-stepping scheme AIAA-811259. In: AIAA 14th fluid and plasma dynamics conference, June 23–25, 1981, Palo Alto, California; 1981. [11] Barth TJ, Jespersen DC. The design and application of upwind schemes on unstructured meshes. AIAA paper 89-0366, 1989. [12] Coirier J. An adaptively-refined, Cartesian, cell-based scheme for the Euler and Navier–Stokes equations. Technical report TM-106754, NASA, 1994. [13] Ralph ME, Pedley TJ. Viscous and inviscid flows in a channel with a moving indentations. J Fluid Mech 1989;209:543–66. [14] Zhao Y, Forhad A. A general method for simulation of fluid flows with moving and compliant boundaries on unstructured grids. Comput Methods Appl Mech Eng 2003;192:4439–66.