Прямое численное моделирование пограничного слоя на плоской пластине с учетом сжимаемости в двумерном случае

При расчете течений воздуха с относительно большими дозвуковыми скоростями возникает необходимость учета сжимаемости. Рассмотрим задачу о пограничном слое на плоской пластине в сжимаемо-вязкой постановке при помощи уравнений Навье-Стокса. Соответствующие уравнения данной задачи приводятся, например, в [1, 2]. Используется численная схема расщепления уравнений Навье-Стокса по физическим процессам. Невязкая сжимаемая подзадача решается с помощью TVD схемы [3, 4]. Вязкая подзадача решается методом, описанным в [5]. Правильный выбор шага интегрирования по времени является важным для этих подзадач. Для невязкой сжимаемой подзадачи шаг по времени выбирался согласно критерию Куранта-Фридрихса-Леви. Для вязкой сжимаемой подзадачи шаг по времени выбирался, как в работе [5]. Так как шаги для данных подзадач не совпадают между собой, использовалась схема интегрирования с разными шагами по времени. В процессе моделирования были получены профили продольной скорости для двух сечений перпендикулярных плоскости пластины для разных чисел Рейнольдса. Результаты сравнивались с профилем Блазиуса для ламинарного случая.

МАТЕМАТИЧЕСКАЯ ФОРМУЛИРОВКА МЕТОДА

Уравнения Навье-Стокса для сжимаемого газа в двумерном случае имеют вид [1]

du    dp  _ d (  du)  d ( (dud p— —---+ 2 — I ц— I +--ц--1I — dt    dx dx V  dx)  dy ( V dy

• 2    d ( (d u d v (

ц|            -I

• 3    d x V V d x d y J_v

dv _— d p d dT” ~d y d x

d u d v d y d x

d ( d v

+ 2 ц — d y V d y

^^^^^^^e

• 2    d ( (du dvH ц1.(1)

• 3    d y V _Vd x d y ))

Для замыкания задачи необходимо добавить уравнения неразрывности dp + d(pu) + d(pv) — 0

dt    dx

и энергии dE +5(Eu) +d(Ev) — div(pv) — 0, (3) dt    dx     dy

где

d z             x d z              x dlv(Pv) — — (PnU + Pnv) + — ^U + p..V) .(4) dx             dy

Здесь

_ d u 2 (d u d v P 11 ^—— P ⁺ 2ц-—-ц — -^т d x 3 V d x d y ^

,

d u 2

d x

du

P12 — P21 — ц , ' , v dy

, d v 2 (d u dv^л

P22 ——P + 2ц-—тц ,   „ dy 3 V dx dy ^ .

Также необходимо уравнение состояния

E — —^7 + ¹ p ( u 2 + v² ) , или p C v T — —^7.(6) к — 1 2 ^{x 7}               к — 1

Коэффициент динамической вязкости для воздуха для (1) и (5) определяется следующим выражением [2]

ц = 1,45 х 10

. -

32 6 T

T + 110

.

c =

Система уравнений (1) решается путем расщепления по физическим процессам. Данная идея не нова и взята из [6], но в данной работе сделан переход от модельной задачи к уравнениям Навье-Стокса.

На этапе конвекции решается задача в невязкой постановке:

5p + 5 ( p u) + 5 ( p v) = 0

51    5 x

5(pu) 5(pu2) 5(puv)

——-+——- + —- = -

Соответствующие правые собственные торы находятся следующим образом [4]

R x (w) =

век-

5 x

5 y

5 x

5 ( p v) 5 ( p vu) 5 ( p v²)     5 p

——- + —-- + ——- = —-

d t

5 x

5 y

,

5 y

u - c

_V H - uc ,

R X (w) =

V ^v J

, R X (w) =

, R 4 (w) =

u

i⁽u² + v^{2 )}

u + c

,

,

5 E + d ((E + p)u) + d ((E + p)v) = 0

_V H + uc ,

I 1

5 x

5 y

с помощью метода Total Variable Diminishing (TVD) [3, 4].

Для удобства изложения схемы метода TVD запишем уравнения Эйлера движения сжимаемого вязкого газа в потоковой форме [1]

d w 5 f_x(w) 5 f (w)

— + ——— + —y--= 0

R y (w) =

u

- c

R y (w) =

u

V н - vc j

R³ y (w) =

,

i⁽u^{2 +} v^{2 )}

,

S t      d x

Здесь

d y

.

V ^u J

w =

p

mx my

V E J

f x^(w) =

mx m2x/ P + P mxmy / P

,

^fy^(w) =

V(E + p)mx / pj , f     my.   1

mxmy / P my/p +p v (E + p)my/ pv

где m_x = p u - импульс по оси x , m_y = p v -импульс по оси y .

Система уравнений (8) замыкается с помощью уравнения состояния для идеального газа, формула (6), где показатель адиабаты для воздуха к = 1.4 .

Собственные значения матриц якобиана [4]

Sf             5f (w)

A x (w) = ^d- x (w ² A y (w) = 4     (11)

dw          dw определяются, как

X _xl = ^{u - c} , X _x 2 = ^u , X x3 = ^u , X x4 = u ⁺ c ,

X _y1 = v - c , X y2 = v , X y3 = v , ^X y4 = ^{v + c} ,(12) где c – скорость звука:

, R⁴ y (w) =

u

_V H + vc ,

, (15)

где H – полная энтальпия:

H = (E + p)/ p = c² /( к- 1) + |(u² + v²) .

Левые собственные векторы-строки Lⁱ , умноженные на вектор столбец конечных разностей переменных (10), можем записать в следующем виде

a i+2,J a2 , .

i+2,j a3  , .

i+2,j a4  , .

i + |.j

= ^Aw i + 2,j , = ^L A w i + 1 J , = ^{L A} w i+1 J , 4

= ^L A w i + lj ,

где левые собственные векторы-строки Lⁱ определяются [4], как

т I / A I 1//U K-1/ 2    2«   1//1 к-1 A

^L x (^w) -I 1> (^{_ +} 7i“(^{u + v} )), ^- )2 (^_+—— ^u), ^-

V c 2c                 c c

к- 1 к- 1 1

2c² , 2c² J

,

T2 / л L к-1, 2   2. к-1 к-1

L x (w) = I 1 -—— (u + v ),^-u,^-v, V 2c            c c

L³ x (w) = ( - v, 0,1,0 ) ,

L>) = 1 Ж"u +K- 1(u² + v²)),^¹--^u), V c 2c                c c

L 1y (w) = I >2 (- + —

V c 2c

к- 1 --Г u, 2c²

^^^^^B

^^^^^B

к- 1 1

c² J

,(17)

к- 1 к- 1 c ^v, c

,

1 к- 1 . к- 1 1

-+—^x) , I c c² 2c 2 J

,

Л к- 1         к- 1 к- 1

L y (w) =1 1 ^- -TT⁽U + v I 2 u 2 v

V 2c            c c

L³ _y (w) = ( - u, 1,0,0 )

к- 1 )

c² J ^,

^L 4x ⁽w) =f Ж^- v + 7-7 ⁽u^{2 +} v²)), V c 2c²

^к- 1     1 ^к- 1

—Г ^u,1 2 ⁽---- ^v)

2c2c c2

к- 1

, 2c2

TVD схема согласно [3, 4] описывается следующими формулами:

wif(t + дО = Wii(t)_~ f*+! -f* I L (19) id V A c/       id        2h Xi +2,j     Xi-1,j fi+ij = 2 ( f(wi,j) + f(wi+1,j) - di+ij) d1+,= —tPtyR^, i+ 2,j        4- "^ i+ 2,j x i + 2,j

A t c k = 1

где

для трех-точечной схемы первого порядка Рое [7]. Данная схема использовалась вблизи границ течения. В случае пяти-точечной схемы второго порядка [3, 4]:

Ptx j = Q^. + Y k 1 X^. - (g kj + g k + 1,j ) . (23) 2 ,                              2 ,                    2 ,                 2 ,                   ,                     ,

В формулах (22) и (23)

v‘ ^A ^ x _t (w„j ,     ⁽²4)

2 A X и ограничителем значения функции в (23):

gkj = sk+1 j max[0, min ( ak+1 j ,ak1 ^^1 j)] / 8 , (25) >J 1+ 2>j                                    i+ 2>j            2>j 1+ 2,j где

S ^k₊, = sign( a^k₊ 1 ) ,            ⁽²⁶⁾

⁺ 2 , ^{j                        +} 2 , ^j

Y k । • _' 11+1-j [(gk+ 1,j -gkj)/ak+ 1,j, ak+ 1,j ^ 0 222 0, ak, _ 0 , (27) I                      , i+ 2,j Q(z) _ z2 + i, (28) wjt + 2a tc) = L*xL*yL*yL>i,j(t), (29)

L* : w*j(t + _At_c) = w, j(t)-—(f*. , . - f* ,

X i, j A c             i, j                               1                ¹

^,j 2hV i ⁺ 2 ^{,j xi} 2 ^,j/ , (30)

t

^Ly ^: W.^t + A t c⁾ = ^wi,j^{(t +} A t c^{) -} — (^fyi + 2 ,j ^{- f}yi - j . (31)

Граничные условия (ГУ) непротекания потока на поверхности пластины для невязкой подзадачи для слоя фиктивных ячеек j = 0, находящихся с «обратной стороны» пластины, сводятся к vi,0 =-vi,1,                 (32)

где ячейка с координатами (i, 1) является прилегающей к поверхности пластины.

Шаг по времени для процесса движения невязкого течения удовлетворяет неравенству atc S 0.48h, c

для критерия Куранта-Фридриха-Леви [3], равного 0,48, который был найден при решении данным методом задачи Рое [7].

Уравнения, описывающие вязкую сжимае- мую подзадачу, имеют вид:

d ( p u) _ 4 дц d u

S t    3 d x d x

d x

. d ²u дцГd u d v)

+ 2ц—₇ + — — + — + d x² d y (d y d x J

+ ц

d u S v |

—⁺ „ „ I ( d y d x d y J

2 Г 5ц d v    Г 5 ² u d ² v )

+ Ц +

3 (d x S y    (d x ² d y d x J_v

d ( p v) _ Зц г S u + d v ' d y     d x ( d y d x ,

+ ц

Г d ²u    d ²v )

----+ —7 + (34)

(S y S x d x² J

4 дц d v _ d ² v + —-— + 2ц—

3 dy dy

2 дц d u

3 ( d y d x

+ ц

d u d v | +2ТГ I

(^dxd y ^d y JJ

dE , d(цu) du , d2u 2 d(цu) (du dv) 2 (d2u d2v = 2^-2+ 2цu2-2l   +|цul + dt      dx  dx      dx2  3  dx  (dx  dyJ  3   (dx2

d(цv) Г du  dv) Г d2u   d2v)  d(цu) Г dv  du) Г d2u

+11| + ц У l1| +11I + Ц U l +I + d x (d y   d x J (d x d y   d x² J   d y (d x   d y J (d y²

, d(цv) dv  , d2v  2 d(цv) Г du

+2      + 2цv— + dy dy dy   3 dy (dx

2 Г d ²u d ²v) - цvl----+

3 (d x d y d y² J

Для конечно-разностной аппроксимации системы уравнений (34, 35) использовались сле- дующие схемы:

d ( p u) _e p j uj-pj^

d            Atd du _ uw2J_Z_ur1J_ dx2

d2u _ uL1,j - 2ut,j + ut+1,j n 2 ^

dx

,

• ^d 2 u ^ u i + 1,j + 1 ^- u i + 1,j - 1 ^- u i - 1,j + 1 ⁺ u i - 1,j - 1

• dxdy              4 a x a y

где A td - шаг по времени для процесса диффузии (вязкой подзадачи), а A x и A y - по пространству. Индекс i отвечает за изменение переменной вдоль оси x, j – вдоль оси y. Остальные производные аппроксимируются аналогичным образом.

Граничные условия прилипания потока на поверхности пластины для слоя фиктивных ячеек j _ 0, находящихся с «обратной стороны» пластины, сводятся к ui,0 _-ui,1  vi,0 _ 0,

где ячейка с координатами (i,1) является прилегающей к поверхности пластины.

В работах [8, 9] было показано, что при выборе шага по времени в виде д td = kh2/ v               (38)

ошибка численного решения уравнения диффузии с помощью схемы «донор-акцептор» и схемы «вперед по времени, центральная по пространству» (ВВЦП) зависит только от константы k и от количества сделанных шагов по времени. Здесь h – шаг ячейки расчетной сетки по пространству h = _д x = _д y . В работе [9] путем серии численных экспериментов было найдено значение к = 0.22 , при котором ошибка решения минимальна. В работе [5] было показано, что данное значение коэффициента k справедливо и для нашей численной схемы.

ТЕСТИРОВАНИЕ ЧИСЛЕННОЙ СХЕМЫ

Рассматривается задача о продольном обтекании плоской пластины и формирования на ней пограничного слоя. Данная задача имеет точное в рамках теории пограничного слоя решение, полученное Блазиусом для ламинарного случая, которое хорошо подтверждается в эксперименте [9]. Известно, что с помощью введения безразмерной нормальной координаты

П = y/' (39) v x и величины продольной компоненты скорости u u = — (40)

профили компоненты скорости u^* в разных сечениях пластины x = const будут совпадать, что является удобным для сравнения получаемых данных.

Сравнение профилей скорости с решением Блазиуса производилось для двух сечений на расстояниях xп » 0.5b и xп » 0.75b от переднего края пластины, где b – длина пластины.

Расчеты проводились в диапазоне чисел Рейнольдса от 10⁴ до 10⁷ . Результаты численного моделирования приведены на рисунках 1-8. Здесь a_x и a_y – размеры расчетной области вдоль координат x и y соответственно, n_x и n_y – число ячеек сетки вдоль координат x и y. Видно, что необходимый шаг расчетной сетки пропорционален квадратному корню из числа Рейнольдса. Это условие вытекает из обратной пропорциональности толщины ламинарного пограничного слоя корню квадратному из числа Рейнольдса.

В результате численного моделирования задачи было выявлено, что в диапазоне чисел Рей- нольдса от 104 до 107 не происходит перехода ламинарного пограничного слоя в турбулентный. При этом никаких ограничений на численную схему не накладывалось. Этот вопрос требует дальнейших исследований. Также оказалось сюрпризом, что TVD схема первого порядка Рое [7] показывает лучшие результаты для ламинарного пограничного слоя по сравнению со схемой второго порядка [4].