Modification of the LS-STAG immersed boundary method for simulating turbulent flows Текст научной статьи по специальности «Медицинские технологии»

DOI: 10.18698/1812-3368-2017-5-19-34

MODIFICATION OF THE LS-STAG IMMERSED BOUNDARY METHOD FOR SIMULATING TURBULENT FLOWS

I.K. Marchevsky V.V. Puzikova

[email protected] [email protected]

Bauman Moscow State Technical University, Moscow, Russian Federation

Abstract

We constructed the LS-STAG discretisation for 2D Reynolds-averaged Navier — Stokes equations, filtered Navier — Stokes equations (as used for large eddy simulation and detached eddy simulation) and equations employed in the Smagorinsky, Spalart — Allmaras, k — e, k — ra and k — ra Menter's Shear Stress Transport turbulence models. We added a fourth grid to the LS-STAG mesh consisting of three staggered grids. We computed the following parameters at the centres of the additional mesh cells: turbulent shear stress and, depending on the turbulence model used, turbulence kinetic energy, turbulent viscosity, and turbulent kinetic energy dissipation rate. We verified the developed numerical method by solving the problem of flow around a circular airfoil when the flow has a high Reynolds number (102 ...107). The obtained results are in good agreement with published experimental data and numerical results of other researchers. Our modification of the LS-STAG immersed boundary method made it possible to model the so-called "drag crisis" phenomenon for a circular airfoil when Re = 105 ...106

Keywords

Immersed boundary method, LS-STAG method, turbulence models, Reynolds-averaged Navier — Stokes equations, large eddy simulation, detached eddy simulation, airfoil

Received 23.01.2017 © BMSTU, 2017

The research is supported by Russian Ministry of Education and Science (proj. 9.2422.2017/PP), Russian Federation President grant for young Russian PhD scientists (proj. MK-7431.2016.8), Russian Foundation for Basic Research (proj. 17-08-01468a)

Introduction. In number of engineering applications, for example in flow simulation around wind turbine rotors, heat exchanger pipes, overhead and underwater cables and pipes, building structures, marine infrastructure elements, etc., it is necessary to solve coupled hydroelastic problems. Such problems are enough difficult for the numerical solution and require high-precision numerical methods usage. There is a special class of numerical methods — the immersed boundary methods — in which the mesh is not connected to the body boundary and is not modified during the entire computation, despite the immersed body movement [1]. These methods involve the rectangular meshes usage. Cells of irregular shape, called the ''cut-cells'', are formed at the intersection of a rectangular mesh with the immersed boundary. One of the most effective methods in this class is the LS-STAG method [2]. This method has not been implemented in both commercial and free software packages.

To obtain accurate quantitative results when simulating unsteady flows characterized by high-speed airfoils movement, and hence by high values of local Reynolds number, strong mesh refinement is need. It leads to a sharp increase in computational cost of direct numerical simulation. The traditional way here is turbulence simulation by using some well-known approaches and turbulence models. However, the corresponding numerical schemes haven't been developed for LS-STAG approach.

In the present research the LS-STAG discretizations for two-dimensional RANS, LES and DES equations and transport equations from Smagorinsky, Spalart — Allmaras, k -s, k -ra and k -ra SST turbulence models are constructed.

Governing equations. The problem is considered for 2D unsteady case when the flow around an airfoil assumed to be viscous and incompressible within the framework of RANS, LES and DES approaches. In contrast to direct numerical simulation (DNS) based on solution of Navier — Stokes equations and resolution of all turbulent movement scales, turbulence models usage involves simulation of turbulence scales contribution to the averaged motion (in case of RANS approach) or simulation of scales that do not exceed the filter width A (in case of LES approach). In case of RANS approach one speaks of the Reynolds stress simulation and in case of LES approach one speaks of the subgrid stress simulation.

The Reynolds-averaged Navier — Stokes equations are being solved in RANS

approach, and the filtered Navier — Stokes equations are being solved in LES approach

instead of the Navier — Stokes equations. DES approach usage means that RANS

equations are being solved in one part of the computational domain, and LES equations

are solved in the other part. It is possible to write down the unified governing equations

in dimensionless variables for all approaches, because the form of LES equations is

similar to the form of RANS equations. So the incompressible flow is described by the

following RANS/LES/DES equations:

Sv . t

V-v = 0, — + (v • V)v = -Vp + vAv + W. (1)

St

Here v = v(x, y, t) = u • ex + v • ey is the dimensionless Reynolds averaged of filtrated velocity, p = p(x, y, t) is the dimensionless Reynolds averaged of filtrated pressure, v = 1/Re is the dimensionless viscosity, xt is the Reynolds or subgrid stresses tensor. The relationship between xt and flow Reynolds averaged or filtrated variables is given by the turbulence model.

Flow around a fixed airfoil in computational domain Q with boundary r = = r1 ur2 ur3 uT4 is considered (Fig. 1). p2p4 In all our simulations, the upstream and

outflow boundaries are set at the distances 8D and 15D, respectively, from the airfoil center, and the blockage ratio is equal to Fig. 1. Computational domain 1/12. The previous studies have shown that

such computational domain is sufficiently wide to obtain results that are don't depend on the domain size.

Fig. 1. Computational domain

We denote the airfoil boundary as K. Then the boundary conditions are the following:

I _ öv

v |Г1 иГ2 иГ3- vœ > ——

ön

-0, v к - vtb.

T4

It is possible to distinguish the linear turbulence scale lturb = kurb (r) for all turbulence models. In the framework of RANS approach this scale lturb is equal to scale Irans = Irans (r), which is determined by the turbulence model (Table 1).

Table 1

Turbulence scale Irans for some turbulence models [3]

Turbulence model Irans Comments

Spalart — Allmaras [4] dw dw is the distance between the field point and the nearest wall

к-e [5] к 3/2e-1 s is the dissipation rate of the turbulent kinetic energy k

к-ю [6], k-ю SST [7] k1/2(ß* ю)-1 ra is the specific dissipation rate of the k, P* = 0.09

In case of LES approach, the scale lturb is equal to subgrid scale:

Iles = Clesa. (2)

Here A = A(r) is the characteristic filter size at the point of computational domain with the radius vector r, and CLES is the empirical constant, which choice depends on the turbulence model and numerical method used to solve the problem (1) in the whole. Within the DES approach the linear turbulence scale lturb is equal to hybrid linear scale

Ides = minims, Cdes A}. (3)

Here CDES is the empirical constant similar to CLES, and the maximum of the mesh steps at the point of computational domain with the radius vector r is used as the characteristic filter size A = A(r). Thus, DES operates as RANS in the domain where the mesh is too coarse and not suitable for resolving turbulent structures, i. e. at CDES A > Irans , and DES operates as subgrid model for LES in the domain where the grid is sufficiently fine [3]. It should be noted that Smagorinsky model is used only within LES approach.

In this paper Eddy Viscosity turbulence models (EVM) are considered. In EVM models the eddy viscosity vt (and the turbulent kinetic energy k in case of two-equation models) is simulated and Reynolds or subgrid stresses are evaluated using the Boussinesq eddy viscosity assumption [3]:

t ^ t du 2, t ^ t dv 2.

Txx = 2vt—+-k; tL = 2vt—+-k; (4)

T Xy = vt

du + dv 8y dx

(5)

Here ■txx and xtyy are normal Reynolds or subgrid stresses and ■txy is shear Reynolds or subgrid stress. In cases of algebraic turbulence models or models with one differential equation the turbulent kinetic energy is assumed to be zero and only eddy viscosity value is computed. For example, in the Smagorinsky model [8] the eddy viscosity is defined by the following formula:

vf =(Cs A)2,

(6)

Here CS is the empirical constant (the Smagorinsky constant). Choice of the CS value depends on the numerical method used to solve the problem, because at LES approach the accuracy of large-scale vortex structures resolution depends not only on the mesh, but also on numerical method properties, in particular, numerical dissipation. If the numerical dissipation is large, it is necessary to choose smaller values of CS, and if numerical dissipation is small, the CS value should be chosen larger.

For EVM with differential equations the governing equations, initial and boundary conditions are given by the turbulence model. In the most general way, they may be written as the following:

Sy

~dt'

rs

y(r, 0) = yo(r), y |r, ur2ur,= V», V k = Vib,

5n r4

Here Prod is the production term which describes the generation of Reynolds or subgrid stresses; Dis is the destruction term; Add is the additional term; y and rv are given by the particular turbulence model (Table 2).

Table 2

Itemization of symbols in (7) and rules for vt computation for some turbulence

models [3]

- (v • V)y = Prod- Dis + V • [(v + )Vy] + Add;

(7)

= 0.

Term Spalart — Allmaras k -e k-ю k-ю SST

V v k e k ю k ю

Add 0 0 0 0 0 0 (1 - F1)DkM

Prod P v Px 1.44ePx k Px 5юРx 9k Px (0.44 + 0.11I\)PX v

Dis D v k3/2 hurb 1.92e2 k k3/2 hurb 3ю2 40 k3/2 / hurb 0.0828 - 0.0078F1

Г„ v vf vf /1.3 vf/2 vf/2 (1 - 0.5!!)у( (0.856 - 0.356Fi)vt

vf v/v1 0.09k2/ e k / ю 0.31k / G1

The following designations are introduced in Table 2: V is the Spalart — Allmaras (S-A) working variable [4]; s is the dissipation rate of the turbulent kinetic energy k; ra is the specific dissipation rate of k;

F - tanh

min

( Г 4k 500v] 3.424k

max<

0.09radw dW® I CDka>dw

4Л

CDk® -max(Dkffl ,10-20); -

1.712Vk -Vra

ю

Pv - 0.1355[1 - ft2]SV; ft2 -1.2e_0.5^2; x -v; S -

V

öu öv öy öx

+f

v2

0.1681/ ,2

turb

fv2-1-"—"; fv1 - ^

1 + X

v1

x3 + 357.911

Dv - (3.2391fw -0.8061ft2 )l —I ; fw -

lturb

; pT - t du t öv t (du öv |_

■; P - Txx~ + xyy + Xxy I "T + T- I;

öx öy ^ öy öx )

s 1/6

65

g6 + 64

; g - r + 0.3(r6 - r);

r - min

Л

-ütt—,10

0.1681Slturb )

(

; G1 - max

0.31ю,

öu öv öy öx

F - tanh

max

l4k 500v

У

0.09ю dw dW ю

Modification of the LS-STAG immersed boundary method. The Cartesian mesh with cells Q,,- = (x,_1, Xi )x (y-_1, y-) is introduced in the rectangular computational domain Q. It is denoted that r,,- is the face of Qi,- cell and xc,- =(xc, ycj) is the center of this cell, which is called "base mesh". Pressure is computed in the center of Q,,-. Unknown components u,,- and v,,- of velocity vector v are computed in the middle of fluid parts of the cell faces. These points are the centers of control volumes Q"- = (xC, x-+1)x (y-_1, y-) (x-mesh) and Qv,- =

= (x,_1, xi) x (yc-, yc-+1) (y-mesh) with faces ru- and rv- and squares Mx and My-, respectively (Fig. 2).

The level-set function 9 = 9(r) = 9(x, y) [9] is used for immersed boundary Fb description [2]:

9(r) <0, r e Qf = Q \ {Qib u rib};

9(r) = 0, r e rib;

9(r)>0, r e Qib.

Vj+1

Vi

г

yj-1

Xi.

Щ-U

hj

Pij

du I öül

Txx\i,j' xyy\i,j

Vi,j-\

Axt

t

-1 '.7

Sul Sol

i

Г".

xi + ]

Fig. 2. Staggered arrangement of the variables on the LS-STAG mesh

The boundary Fb is represented by a line segment on the cut-cell Qi,j (Fig. 3). Location of this segment endpoints is defined by a linear interpolation of the variable ,j = 9(x;, yj). The cell-face fraction ratios 9"j and 9J,j are introduced [2].

yj

yi-\j

Vj- 1

ф«-и>0

Ф,;у>0

P ib

1 Uj

1

%J

ui-\,j

9i-lj-l<0

PiJ

%j-1,

xi- 1

<P;j-1 < 0

^ 0/j-1 ^

6

\ QljAyj

Fig. 3. Example of the cut-cell on the LS-STAG mesh

They take values in interval [0,1] and represent the fluid parts of the east and north faces of r, j, respectively. The cell-face fraction ratios are defined by a one-dimensional linear interpolations of function ^(xi, y) in interval [yj_1, yj ] and 9(x, yj) in interval [xi _1, xi ]:

min(^i, j _i, ^i, j) _ min(^i _i, j, ^i, j)

0u i =-

'' min^, j-1, ф;, j) - max(9i, j-1, ф;, j)

; = —

'' mrn^ -1, j, ф;, j) - max^ -1, j, ф;, j)

In 2D case, the cut-cells can be classified into trapezoidal, triangular and pentagonal cells. Examples of each type cut-cells are shown on Fig. 4.

To preserve the five-point structure of the stencil of the MAC method, which can be considered as some kind of ''predecessor'' of the LS-STAG method, we need to make distinction between the discretization of the normal and shear stresses [2] (Fig. 4). It is conveniently to sample the eddy viscosity vt, the turbulent kinetic

du I 8vI

Pi,j[ tokj&kj

"i-lj

f

Km

lb. . v? . Tt I . .

hP hP ХУ PJ

8u\ dv\

Sy\ij,Sx\ij

x' I- ■ x' I

du I 8v

,8x\ijMi,r

x' I- • хг I- ■ ■ViJ- 1

k- ■ v'- -x' I- •

hP hP Xy\VJ

8u\ 8vI tykj&hj

"hi

k- ■ V? хг I •

hP hP xy I ''

5и| dv\ ду\и,8х\и

Fig. 4. Location of the variables discretization points in case of base types cells of the LS-STAG

mesh:

a — Cartesian Fluid Cell; b — North Trapezoidal Cell; c — Northwest Pentagonal Cell; d — Northwest

Triangle Cell

energy k, the dissipation rate of the kinetic energy, s and the specific dissipation rate of kinetic energy a at the same points as the shear stresses. Thus, in case of the LS-STAG method usage for RANS-based models the fourth mesh (xy-mesh) with cells =(xC ,xC+1)x (yc, yc+1) is needed. The faces of these cells are r*y (Fig. 2) and their

areas are MXy. If i and j take values from ranges 1,..., N and 1,..., M, respectively,

the base mesh contains E = NM cells, x-mesh contains Ex =(N -1)M cells, y-mesh contains Ey = N (M -1) cells and xy-mesh contains Exy = (N - 1)(M -1) cells. It is possible to assign a weight ai,j to each cell Qi,j of the base mesh:

0, if Qi,j is the solid cell; ai,j = <! 1/ 3, if Qi,j is the triangular cell; 1/4, otherwise.

Then Mxy can be expressed through the area of base mesh cells:

Mxy = a*, j_1Vi> j_i +ai _i, jVi _i, j +a,i, ¡Vi, j +a,-, j+V, j+i.

Here V, j is the area of the cell Q*, j.

Since vt and shear Reynolds or subgrid stresses (5) are sampled at the same points, it follows that

IXy \i,j Vi,j

du dv Л

v°y i, ! __ dx *',j J

(8)

whereas averaged values of turbulent viscosity vtj and the turbulent kinetic energy kj,j should be used for the computation of the normal Reynolds or subgrid stresses

(4):

it \.. = 2V ■ ^

Lxx \t,j , -

dx

., J

2- t . _ t dv

+3 kt,J ; 1 yy \t-,J = 2v,,J —

., j

+1 k. + 3kt ' ;

/itj = at,J (V,J + v.,j+vt_1,J +vt_1,j_j); kt,J = at,j (kt,j + k.,j-i + kt-i,j + k-i,j-i). (9)

Formulae for normal stresses —

Ox

ov

and —

computation are the following:

', j

du Ox

4 jUi, J-ej_i, u-1, J + (0U-1, J _eu J X,

j

Vi,J / Ayj

dv dy

ev, v J -ev, J _ivt, J _i+(0v, J _i-ev, J )vtbj

j

Vi,j / Ах,

formulae for shear stresses —

dy

cell:

', j

. dv and —

dx

computation depend on the type of Ц,

xy

', J

du dy

j

ut, j+i ut, j

(0u ,+iAyj+i +0u j Ay, )/2 du

dy du

dy

, ut,j+i - u(xt, y'j^j+i) .. " 0uj+iay,+i/2 ' _ u(Xi, yb ) - u,,j

', j

eu, j Ay, /2

if eu j * o, eu j+i * 0;

if euj = o, eu,,+i *0;

if eu j * o, eu j+i = 0;

ÔV

dx

>> j

Vi+l,j - Vi,J

(0V+1> j Axi+i +9V, j AXi )/2

">,w 1,. . -MÎv'b

if 0V,j * 0, 0V+1,j * 0;

ÔV

dx

ÔV

dx

i, j

i, j

-Vi+0j - V(xi+l,J, ^ ), if 0Vj =0, 00 * 0;

0V+i, j Axi+i/2 "j i+1, j

V(xfbj, yj ) - Vi, j 0V,jAxt /2 ,

if 0V,j * 0, 0;+,, j = 0.

i+i, j

Components of hydrodynamic force acting on the immersed boundary can be computed as the following:

f

F =

1 xa

X .,(0"-i,j--0Uj)Ayj

cut-cells q'

Pi, j -V

du dx

\

- v Quad;

i, j

i, j

i ,j

du

— I

dy

F =

1 ya ~

cut-cells Qib.

i, J

, ôv

- v Quadibj ' Ô 1 dx

ex nj + (0v,j-i -0 V,j )Axi

(

pi,j -v

d dy

i ,j y

Here Fxa is the drag force, Fya is the lift force, Axt = Xi - Xi_i, Ayj = yj - yj_i, Quadi^j is the quadrature of the shear stresses. Quadibj has to be adapted to each type

of cut-cells. This quadrature is based on the location of point where the shear stresses are sampled in Fig. 4 and the trapezoidal rule.

It is conveniently to sample the linear turbulence scale lturb and the characteristic filter size A for LES and DES at the same points as the vt and k. We recall that the maximum mesh step at the given point of the computational domain is used as a filter size A for DES approach. Since we deal with xy-mesh, the characteristic filter size is defined as a following:

Ai, j = A- =max{Ayx-i, J , Ay^, AyX*, J , Aj Ax%, Ax^},

where

Ayxy =2 (out j- Ayj +eUj+iAyj+i); AxXy = 2 (er j AXi + 0^ Ax,+i). Within LES approach the following filter can also be used on the LS-STAG mesh:

Ai, j = AvoJ =Jm~x;

According to the concept of the LS-STAG method, equations (1), (7) should be written in integral form for cell of base mesh, cell of x-mesh, cell of y-mesh and cell of xy-mesh, respectively:

J v - ndS = 0;

(i0)

li, j

— J udV + J (v • n)udS+ J pex • ndS - J vVu • ndS -

dt nU rU yU rU

"г, j 4 j 4 j 4 j

- J xXxex • ndS - J xXyey • ndS = 0;

yU yU

i, j i, j

d J vdV + J (v • n)vdS+ J pey • ndS - J vVv • ndS -

dt nV yV yV rV

ni, j 4 j 4 j 4 j

(10)

d_ dt

- J xtyyey • ndS - J xtxyex • ndS = 0;

yV yV

i,j i,j

J ydV + J (v • n)ydS = J (v + yv)(Vf n)dS + J (Prod-Dis+ Add)dV. (11)

„xy yxy yxy nxy

, j 4 j j ni, j

In case of fixed immersed boundaries by analogy with the LS-STAG discretization of Navier — Stokes equations [2] the general form of the LS-STAG discretization for (10), (11) can be written as the following:

DxUx + DyUy + U>b =0;

— (MxUx) + CxUx + Gx(P -Txx)- DXTxy -vKxUx + S'b,c-vSXb,v =0; (12) di

— (MyUy) + CyUy + Gy (P - Tyy) - DyTxy - vKyUy + Syb,c - vSyb,v = 0; dt

— (Mxy y) + Cxy [Un ]Vprn + S*f,n [ y», Vb ] - Kxy [G^ ] ]T xvn - , t _

dt (13)

- Siy [Gy, y», y'b ] - MxyPDAn = 0.

Here P is the discrete pressure, Ux and Uy are the discrete components of the velocity vector, Txy is the discrete shear Reynolds or subgrid stresses, Txx and Tyy are

the discrete normal Reynolds or subgrid stresses, y is the discrete y, fv is the

—ib . ... discrete ry; U denotes the mass flux arising in case of vib ^0; Sx,c, Sxb,v, Sy,c,

Syb,v, Slxyc, Sxy are source terms; Kx, Ky and Kxy represent the discretization of the

diffusive terms; Dx, Dy, Dx, Dy are the divergence discrete analogues on the

corresponding meshes; Cx, Cy and Cxy represent the discretization of the convective

terms; Gx = -Dj and Gy = -Dj are the gradient discrete analogues; PDA is the

discrete analog of (Prod - Dis + Add); Gy = Gxyy + Sib,g[y», yib ] is the discrete

analogue of dy/ dx and dy / dy which are computed in the middle of Q,j fluid

faces.

The time integration of the differential algebraic system (12) is performed with a semi-implicit projection method based on the Adams — Bashforth/second-order

backward differentiation formula (AB/BDF 2) scheme. Predictor step leads to discrete analogues of the Helmholtz equation for velocities prediction Ux, Uy at the time tn+1 = (n + 1)At:

U _ ATjn + Tjn-1

Mx —x-x-— + 2CU + 2SX'cn _ CrUr1 _ SXb'cn-1 _

2 At

(14)

-DT (Pn _ Tn) _ Dl'Ty _vKxUx _vS'xb,v = 0;

3Uy _ 4U"n+m _ b , , b ,

My —y—^—y—+2Cnun + 2Sf'cn _ q_ ^1 _ St-1 _

_DT (Pn _Tyy)_ DTy'nTXny _vKyUy _vSyb,v = 0.

Here At is the constant time discretization step. Corrector step is the following:

3 n-n+1 _TT 3 U"+1 _U

-M -x--_DT(Pn+1 _Pn) = 0; -M -—_D (Pn+1 _Pn) = 0;

2 x At x 2 y At y

, , —;b,n+1

DUn+1 + D U + U =0.

xx y y

It leads to the following discrete analogue of Poisson equation for pressure function O = 2At (Pn+1 _ Pn)/ 3:

AO = DxT7X +DyU y +U'M+1, (15)

A = _Dx (Mx )_1(Dx )T _Dy (My )_1(Dy )T. Then flow variables at the time point tn+1 are computed by the following formulae:

3O

Un+1 = Ux + M-_ 1DT O; Ul+1 = Uy + M_ D O; Pn+1 = — + Pn. (16)

2 At

After this, new values of Reynolds or subgrid stresses Txx+1, Tyy+1, Txy+1 are computed according to (8), (9) by solving the discrete analogues of the equations from the used turbulence model:

Mxy ynA_ yn + Cxy[Un ]yn + S%cn , y'b ] _ Kxy[Gy]f nv _ (17)

'v [Gn 'xy I ^y >

-SXy[Gy, , угЬ ]- MxyPDAn =0.

Numerical experiments. The flow past circular airfoil was simulated using the developed modification of the LS-STAG method. The time averaged drag coefficient CD and the Strouhal number St were computed. The coefficient CD was obtained by

FXa (t )

averaging over a large period of time the unsteady load CD (t) =

pVo2/2'

Flow simulation at low Reynolds numbers (Re<103). The LS-STAG method allows to receive the results are in good agreement with the experimental and computational data, even on very coarse meshes (Table 3). No turbulence models have been used.

Table 3

Comparison of CD and St at Re = 100 and Re = 200 with known experimental and numerical results from the literature

Source Re = 100 Re = 200

CD St CD St

Zdravkovich [10] (experiment) 1.21...1.41 0.16.0.17 - -

LS-STAG (present study, 120 x 148) 1.31 0.17 1.33 0.20

LS-STAG (present study, 240 x 204) 1.32 0.17 1.33 0.20

LS-STAG (present study, 480 x 408) 1.32 0.17 1.33 0.20

Cheny [2] (LS-STAG, 550 x 350) 1.32 0.17 1.33 0.20

Henderson [11] 1.35 0.16 1.34 0.20

He [12] 1.35 0.17 1.36 0.20

Flow simulation at medium Reynolds numbers (Re = 103...104). The flow was simulated at the Reynolds numbers Re = 1000 (on non-uniform meshes 120 x148 with At = 5-10"2 and 240 x 296 with At = 10"3) and Re = 3900 (on non-uniform meshes 120x148 with At = 10"3 and 240x296 with At = 5-10"4). These values of the Re were chosen because the experimental data [10, 11] and results of other researchers [14-17] are known for them. Computational results are shown in Table 4 and in Fig. 5. These results are in good agreement with experimental data for simulation on coarse meshes by using the proposed modification of the LS-STAG method. But since the considered models works well only at high Reynolds numbers, it is hardly possible to improve the numerical results, for example, by mesh refinement. In our opinion the wall functions usage can be an efficient solution to this problem.

Table 4

Comparison of CD and St values with known experimental and numerical results

Turbulence model Number of cells Re = 1000 Re = 3900

CD St Cd St

Experiment [10] - 0.98 0.21 0.93 0.22

Experiment [13] - 1.12 - 1.01 -

Smag., LES, Cs = 0.1 [14] 1 103 520 - - 1.08 -

Smag., LES, Cs = 0.1 (spectral) [15] 30 720 - - 1.01 0.22

Smag., LES, CS = 0.1 (finite volume) [15] 855 040 - - 1.07 0.24

k -e, RANS [16] 46 304 1.00 0.15 1.00 0.15

Real k - e, RANS [16] 46 304 - 0.17 - 0.20

k-ra SST, RANS [16] 46 304 - 0.23 - 0.25

k - e, RANS [17], ANSYS 388 550 1.12 - 0.74 -

k -ra SST, RANS [17], ANSYS 388 550 0.99 - 0.62 -

End of Table 4

Turbulence model Number of cells Re = 1000 Re = 3900

cD St cD St

Smag., LES, Cs = 0.1 [17], ANSYS 388 550 1.15 0.21 1.07 -

Smag., LES, A = Amax, CS = 0.2, present study 71 040 1.35 0.24 1.11 0.26

Smag., LES, A = Amax, CS = 0.5, present study 71 040 1.37 0.25 1.10 0.25

S-A, RANS, present study 71 040 1.37 0.25 1.13 0.25

S-A, DES, CS = 0.7, present study 71 040 1.37 0.25 1.11 0.25

k -s, RANS, present study 71 040 1.36 0.25 1.23 0.28

k -s, LES, A = Amax, CS =0.9, present study 71 040 1.37 0.25 1.11 0.25

k -ra, RANS, present study 71 040 1.32 0.24 1.18 0.24

k-ra, DES, CS =1.0, present study 71 040 1.32 0.25 1.00 0,25

k -ra SST, RANS, present study 71 040 1.34 0.25 1.14 0.25

a b

Fig. 5. Computed unsteady load CD (t) and CL (t) (RANS, k -ra model, mesh 240 x 296):

a — Re = 1000 ; b — Re = 3900

Fig. 6. Comparison of the drag coefficient computed values with experimental [18] and computational data on meshes M1-M5 [19] and [20]:

- — Wieselsberger [18]; ------------ — Henderson, 2D calculations [20]; O — mesh M1; +--mesh

M2; □ — mesh M3; x — mesh M4 with Smagorinsky model; A — mesh M4; * — mesh M5; • — LS-STAG (k - ra, RANS, 240x296); O — LS-STAG (480x592)

Flow simulation at high Reynolds numbers (Re = 105 ...107).The flow was simulated at the Reynolds numbers Re = 105 ...107 (on non-uniform meshes 240 x 296 with At = 5-10_4 and 480 x 592 with At = 10_4). Results obtained on mesh 480 x 592 are very close to experimental data [18], see Fig. 6. At Re = 2 -10s, the boundary layer on the cylinder surface undergoes a transition from laminar to turbulent [19]. This transition leads to a delay of the separation of flow from the cylinder surface causing a substantial reduction in the drag force. This is often referred to as ''drag crisis''. This phenomenon was simulated by using modified LS-STAG immersed boundary method (Fig. 6).

Conclusions. The key points of the LS-STAG method [2, 21] extension for RANS/LES/DES turbulence models were described. For the shear Reynolds stresses and for the eddy viscosity an additional mesh (xy-mesh) is introduced. The general approach to the construction of the LS-STAG discretization for differential equations of the EVM models on the additional xy-mesh shown. The Smagorinsky, Spalart — Allmaras, k _ s, k _ ra and k _ ra SST turbulence models are considered. To validate modified LS-STAG immersed boundary method the flow past a circular airfoil at Re = 102 ...107 was simulated. Computational results are in good agreement with established results from the literature. Also, the so-called ''drag crisis'' phenomenon of circular cylinder at Re = 105... 106 was simulated.

REFERENCES

[1] Mittal R., Iaccarino G. Immersed boundary methods. Ann. Rev. Fluid Mech., 2005, vol. 37, pp. 239-261. DOI: 10.1146/annurev.fluid.37.061903.175743

Available at: http://www.annualreviews.org/doi/abs/10.1146/annurev.fluid.37.061903.175743

[2] Cheny Y., Botella O. The LS-STAG method: A new immersed boundary/level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties. J. Comput. Phys., 2010, vol. 229, iss. 4, pp. 1043-1076.

DOI: 10.1016/j.jcp.2009.10.007

[3] Spalart P.R. Strategies for turbulence modelling and simulations. Int. J. Heat and Fluid Flow. 2000, vol. 21, iss. 3, pp. 252-263. DOI: 10.1016/S0142-727X(00)00007-2

[4] Spalart P.R., Allmaras S.R. A one-equation turbulence model for aerodynamic flows. Recherche Aerospatiale, 1994, no. 1, pp. 5-21.

[5] Jones W.P., Launder B.E. The prediction of laminarization with a two-equation model of turbulence. Int. J. Heat Mass Transfer, 1972, vol. 15, iss. 2, pp. 301-314.

DOI: 10.1016/0017-9310(72)90076-2

[6] Wilcox D.C. Reassessment of the scale-determining equation for advanced turbulence models. AIAA Journal, 1988, vol. 26, no. 11, pp. 1299-1310. DOI: 10.2514/3.10041

Available at: https://arc.aiaa.org/doi/abs/10.2514/3.10041?journalCode=aiaaj

[7] Menter F.R. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal, 1994, vol. 32, no. 8, pp. 1598-1605. DOI: 10.2514/3.12149

Available at: https://arc.aiaa.org/doi/abs/10.2514/3.12149?journalCode=aiaaj

[8] Smagorinsky J. General circulation experiments with the primitive equations. I. The basic experiment. Monthly Weather Review, 1963, vol. 91, no. 3, pp. 99-164.

DOI: 10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2 Available at: http://journals.ametsoc.org/doi/abs/10.1175/ 1520-0493%281963%29091%3C0099%3AGCEWTP%3E2.3.C0%3B2

[9] Osher S., Fedkiw R.P. Level set methods and dynamic implicit surfaces. New-York, Springer, 2003. 273 p.

[10] Zdravkovich M.M. Flow around circular cylinders. Oxford University Press, 1997. 694 p.

[11] Henderson R.D. Nonlinear dynamics and pattern formation in turbulent wake transition. J. Fluid Mech., 1997, vol. 352, pp. 65-112. DOI: 10.1017/S0022112097007465

[12] He J.W., Glovinski R., Metcalfe R., Nordlander A., Triaux J.P. Active control and drag optimization for flow past a circular cylinder. Part I: Oscillatory cylinder rotation. J. Comput. Phys., 2000, vol. 163, iss. 1, pp. 87-17. DOI: 10.1006/jcph.2000.6556

[13] Zahm A.F. Flow and drag formulae for simple quadrics. NACA Technical Report 253. Washington, US Government Printing Office, 1927. 29 p.

[14] Breuer M. Large eddy simulation of the subcritical flow past a circular cylinder: Numerical and modelling aspects. Int. J. Numer. Meth. Fluids, 1998, vol. 28, iss. 9, pp. 1281-1302.

DOI: 10.1002/(SICI)1097-0363(19981215)28:9<1281::AID-FLD759>3.0.CO;2-#

[15] Blackburn H.M., Schmidt S. Large eddy simulation of flow past a circular cylinder. Proc. 14th Australasian Fluid Mechanics Conf., Adelaide University, 2001, pp. 689-692.

[16] Rahman M.M., Karim M.M., Alim M.A. Numerical investigation of unsteady flow past a circular cylinder using 2-D finite volume method. J. Naval Arch. and Marine Eng., 2007, vol. 4, no. 1, pp. 27-42. DOI: 10.3329/jname.v4i1.914

[17] Patel Y. Numerical investigation of flow past a circular cylinder and in a staggered tube bundle using various turbulence models. Master's thesis. Lappeenranta University of Technology, 2010. 87 p.

[18] Wieselsberger von C. Neuere feststellungen uber die gesetze des flussigkeits und luftwiderstandes. Phys. Zeit, 1921, vol. 22, pp. 321-328.

[19] Singh S.P., Mittal S. Flow past a cylinder: Shear layer instability and drag crisis. Int. J. Num. Meth. In Fluids, 2005, vol. 47, iss. 1, pp. 75-98. DOI: 10.1002/fld.807

[20] Henderson R.D. Details of the drag curve near the onset of vortex shedding. Physics of Fluids, 1995, no. 7, pp. 2102-2104. DOI: 10.1063/1.868459

Available at: http://aip.scitation.org/doi/10.1063/L868459

[21] Marchevsky I.K., Puzikova V.V. Flow-around simulation of circular cylinder performing rotary oscillations by LS-STAG method. Vestn. Mosk. Gos. Tekh. Univ. im. N.E. Baumana, Estestv. Nauki [Herald of the Bauman Moscow State Tech. Univ., Nat. Sci.], 2014, no. 3, pp. 93-107 (in Russ.).

Marchevsky I.K. — Cand. Sc. (Phys.-Math.), Assoc. Professor of Applied Mathematics Department, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, str. 1, Moscow, 105005 Russian Federation).

Puzikova V.V. — Cand. Sc. (Phys.-Math.), Assist. Professor of Applied Mathematics Department, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, str. 1, Moscow, 105005 Russian Federation).

Please cite this article in English as:

Marchevsky I.K., Puzikova V.V. Modification of the LS-STAG Immersed Boundary Method for Simulating Turbulent Flows. Vestn. Mosk. Gos. Tekh. Univ. im. N.E. Baumana, Estestv. Nauki [Herald of the Bauman Moscow State Tech. Univ., Nat. Sci.], 2017, no. 5, pp. 19-34. DOI: 10.18698/1812-3368-2017-5-19-34

И.В. Блудова, Э.Н. Белянова

Начала топологии в примерах и задачах

киздмшьство

В Издательстве МГТУ им. Н.Э. Баумана вышло в свет учебное пособие авторов И.В. Блудовой, Э.Н. Беляновой

«Начала топологии в примерах и задачах»

Рассмотрены различные классические примеры топологических и метрических пространств и непрерывных отображений, сформулированы все необходимые топологические определения и утверждения. Читателям предложено самостоятельно доказать некоторые свойства указанных выше топологических и метрических пространств, а в случае недостаточной успешности попыток получить эти доказательства—узнать подробные решения предложенных задач.

По вопросам приобретения обращайтесь:

105005, Москва, 2-я Бауманская ул., д. 5, стр. 1

+7 (499) 263-60-45

[email protected]

www.baumanpress.ru