CC BY
8
УДК 532.5.032
ПРЯМОЕ ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПОГРАНИЧНОГО СЛОЯ НА ПЛОСКОЙ ПЛАСТИНЕ С УЧЕТОМ СЖИМАЕМОСТИ В ДВУМЕРНОМ СЛУЧАЕ
© 2021 В.В. Никонов
Самарский национальный исследовательский университет имени академика С.П. Королёва,
Самара, Россия
Статья поступила в редакцию 22.09.2021
В статье моделируется пограничный слой на плоской пластине с учетом сжимаемости воздуха. Учет сжимаемости производился с помощью ТУБ схемы. Был получен профиль продольной скорости для двух сечений перпендикулярных плоскости пластины для разных чисел Рейнольдса. Результаты сравнивались с профилем Блазиуса для ламинарного случая.
Ключевые слова: воздух, сжимаемость, уравнения Навье-Стокса, пограничный слой, прямое моделирование, численное моделирование, двумерное течение. Б01: 10.37313/1990-5378-2021-23-6-94-102
При расчете течений воздуха с относительно большими дозвуковыми скоростями возникает необходимость учета сжимаемости. Рассмотрим задачу о пограничном слое на плоской пластине в сжимаемо-вязкой постановке при помощи уравнений Навье-Стокса. Соответствующие уравнения данной задачи приводятся, например, в [1, 2]. Используется численная схема расщепления уравнений Навье-Стокса по физическим процессам. Невязкая сжимаемая подзадача решается с помощью ТУБ схемы [3, 4]. Вязкая подзадача решается методом, описанным в [5]. Правильный выбор шага интегрирования по времени является важным для этих подзадач. Для невязкой сжимаемой подзадачи шаг по времени выбирался согласно критерию Куранта-Фридрихса-Леви. Для вязкой сжимаемой подзадачи шаг по времени выбирался, как в работе [5]. Так как шаги для данных подзадач не совпадают между собой, использовалась схема интегрирования с разными шагами по времени. В процессе моделирования были получены профили продольной скорости для двух сечений перпендикулярных плоскости пластины для разных чисел Рейнольдса. Результаты сравнивались с профилем Блазиуса для ламинарного случая.
МАТЕМАТИЧЕСКАЯ ФОРМУЛИРОВКА МЕТОДА
Уравнения Навье-Стокса для сжимаемого
газа в двумерном случае имеют вид [1]
( (
Д
du Эр „ 8 ( 8u ^ 8 р— = —^ + 2 — 1 д— | + — dt 8x 8x V 8x J 8y
V V
8u 8v — + —
8y 8x
JJ
2 5
3 dx
Д
dv = 8p +8 dt 8y 8x
v
f r
д
V V
3u dv --1--
cx dy
8 u + 8v 8y 8x
8
+ 2— 8y
dv
2 д( ( dv ^ ---ц--1--
3 dy I I 3x dy
(1)
V ^/у
Для замыкания задачи необходимо добавить уравнения неразрывности ф 5(ри) д^)
- + -8t 8x
8y
= 0
(2)
и энергии
8E 8 (Eu) 8(Ev)
8t 8x
8y
- div(Pv) = 0, (3)
где
8 8
div(Pv) = — (p„u + Pi2v) + — (P2iu + P22v) -(4)
8x 8y
Здесь
P11 =-p + 2д —--д 8x 3
8u 2 (8u 8v ^
Pl2 = P21 = Д
8x 8y J
'8u + 8vч
V8y 8x J' 8v 2 (8u 8v ^
(5)
P22 =-P + 2Д^Т-ТД д. д.
8y 3 V8x 8y J
Также необходимо уравнение состояния
E=¿1+2 p(u2 + v2).
или pCvT =
Никонов Валерий Владимирович, кандидат технических наук, доцент кафедры конструкции и проектирования летательных аппаратов. E-mail: v_nikonov@mail.ru
".(6)
К-1 2 4 К-1
Коэффициент динамической вязкости для воздуха для (1) и (5) определяется следующим выражением [2]
^ = 1,45 х10-
T
3/2
(7)
T +110
Система уравнений (1) решается путем расщепления по физическим процессам. Данная идея не нова и взята из [6], но в данной работе сделан переход от модельной задачи к уравнениям Навье-Стокса.
На этапе конвекции решается задача в невязкой постановке:
Ф | д(pu) + d(pv) = 0 д t дх
д y
д(pu) | д(pu2) + д(puv)
д t
д x
д y
д p
д x
д (pv) + 5(pvu) + д (pv2)
д p
д y
дt дх ду дЕ д((Е + p)u) д ((E + p)v)
= , (8)
д t
д х
д у
= 0
5t - -i-- --1— дx
Здесь
Г p'
m
w= x fx(w) =
m y
V E y
ду
= 0
(9)
fy(w) =
mx
mX/p+p
mxmy/ p
(E + p)mx/ p
л
mxmy/ p m2/p+p
(E + p)my/ p
импульс по оси у.
Система уравнений (8) замыкается с помощью уравнения состояния для идеального газа, формула (6), где показатель адиабаты для воздуха к = 1.4.
Собственные значения матриц якобиана [4]
= —^ Ау^) = (11)
определяются, как
Хх1 = и - с, Хх2 = и, ХХ3 = и , Хх4 = и + с,
Xу1 = V - с, Xу2 = V, Xуз = V, ху4 = v + с ,(12)
где с - скорость звука:
c=fp •
(13)
Соответствующие правые собственные векторы находятся следующим образом [4]
( 1 ^ г 1 ^
RX(w) =
RX(w) =
u - c
v
H - uc
Г 0 ^ 0
, RX(w) =
y
1
V v y
, R4x(w) =
с помощью метода Total Variable Diminishing (TVD) [3, 4].
Для удобства изложения схемы метода TVD запишем уравнения Эйлера движения сжимаемого вязкого газа в потоковой форме [1]
dw dfx(w) dfy(w)
Ry(w) =
R3(w) =
Г 1 ^
u
v - c H - vc
Г о ^ i
u
v
vi(u2+v2)
Г 1 ^
u + c
v
V H + uc y Г 1
y
(14)
Ry (w) =
y
0
v u y
. Ry(w) =
u
v
V!2(u2 + v2)j
Г 1 ^
u
v + c H + vc
, (15)
y
где Н - полная энтальпия: Н = (Е + р)/р = с2/( к-1) + ^(и2 + V2).
Левые собственные векторы-строки П, умноженные на вектор столбец конечных разностей переменных (10), можем записать в следующем виде
a1. 1 . = 1 = L . w. 1 .
i+ij л i+i2,j
(10) a2. 1 . = £ 2
i+ij A i+t,j
a3 1 = £ 3
A i+t,j
pv - a4. 1 . = £ 4
i+ij A i+2,J
(16)
где левые собственные векторы-строки П определяются [4], как
5
Т2, ч (л К-К 2 . К-1 К-1 К-1 ^
(ш) = 11 - ~2с2~(и + ^'"С"»'"С"С2")
L3x(w) = (- v, 0,1,0),
(17)
, u к 1 . л л.. í.1 к 1 , к 1 к 1 i
L4x(w) = i 1(— (u2 + v2)), 1---2u), - —v —
tW -> I i/Л К - 1 2 2 К-1 1 К-1 К-1 ^
Ly(w) +—(u + v W,-^,-^ + —v)T
\ I 1 л К 1 К 1 К 1
L,(w) = I1"^-T(u2 + ---
2c - - -
L3y(w) = (-u, 1,0,0)
2h 3 k
dl+2, = "Г :1,jRXl+|,j '
Д1 c k=1
где
ß\ = pl+ij
V
1+2,j
ak j .
1+i,
В формулах (22) и (23)
V k = ,„
1+2'j 2 x x 1+
^xk(w,+1j) ,
g. =Si+1jmax[0,m1n
где
alk+1 . 1+1 j
s. j = s1gn(a. j ),
1+-2,j ö v 1+ijy
^1+2.j =
(gk+1 ■- gl,)/ak+, ■ , ak+, 0
vo1+1,j с»1,,/ 1+1, j ' 1+1, j
0, ak+, . = 0
' 1+i,j
= + i
Q(z) = z2 + ^ ,
:w1,(t +Дtc) = w1,(t)-A(( -f* , .)
Д '"0
tc < 0.48-,
(33)
(18)
Lt(w) = (X(-I + |-I(u2 + v2))-- ^T1£
TVD схема согласно [3, 4] описывается следующими формулами:
W1,j(t +д tc) = W1,j(t)(-ГJ(19)
fj 2(f(wiJKf(w1+1,j)-d1+1 j) , (20)
(21)
для критерия Куранта-Фридриха-Леви [3], равного 0,48, который был найден при решении данным методом задачи Рое [7].
Уравнения, описывающие вязкую сжимаемую подзадачу, имеют вид:
ö(pu) 4 öu öu „ ö2u öu/öu öv i =--tH— + 2u—- + — + — +
dt 3 öx öx öx öy ^öy öxy
f a2
ö2u ö2v
2
öy öxöy
2föuöv (ö2u ö2v I
+ -
ö(pv) öuYöu öv i ( ö2u ö2v ^ v - -+— + u-+
(22)
для трех-точечной схемы первого порядка Рое [7]. Данная схема использовалась вблизи границ течения. В случае пяти-точечной схемы второго порядка [3, 4]:
- =Q(vk+ к+ 1,-). (23)
öy öx 4 öu öv
öx öy ^öx öyöx j j + (34)
öy öx j ^öyöx öx j
„ ö2v 2
+ —:--+ 2u—T -
3 öy öy öy 3
öu öu
öy öx
f л2 ~a W
ö u ö vi
öxöy öy2
jj
(24)
öE „ö(uu) öu „ ö2u 2 ö(uu) (öu öv i 2 (ö2u ö2v | — = 2^^-—+ 2цц—2---—+— I—ul —г+-1+
ö öx öx öx 3 öx ^öx öyj 3 ^öx öxöyj
+ ö(|iv)(öu + öv |+ ( ö2u + ö2v |+ö(uu)(öv + öu |+ | ö2u + ö2v |+
öx ^öy öx j löxöy öx2) öy ^öx öyj l öy2 öxöy)
(35)
и ограничителем значения функции в (23):
-k =k , ak1|sk+1|)]/s, (25)
2 2'J /
. ö(uv) öv . ö2v 2 ö(uv)) +2 — + 2uv—г —
öy öy öy 3 öy löx öy
öu + öv
\ 1 (л- л- \
2 ö u ö v — uv -+ —2
3 löxöy öy
(26)
,(27)
(28)
Для конечно-разностной аппроксимации системы уравнений (34, 35) использовались следующие схемы:
ч t +. tt +. t t t
Р1,4и1,4
öt
Д td
öu Ä u1+1,j — u1-1,j
öx
2 Д x
(36)
wM(t + 2д tc) = L*xL*y L*yL*x w1, j (t), (29)
ö2u ut-!,j - 2utJ+ ut+1,j
öx2
x
,(30)
L*í: w;:;(t + 4У = w^(t + 41с) - А-. (31)
Граничные условия (ГУ) непротекания потока на поверхности пластины для невязкой подзадачи для слоя фиктивных ячеек j = 0, находящихся с «обратной стороны» пластины, сводятся к
=-^1, (32)
где ячейка с координатами (1, 1) является прилегающей к поверхности пластины.
Шаг по времени для процесса движения невязкого течения удовлетворяет неравенству
ö 2u
ut+1,j+1 - u1+1,j-1 - u1-1,,+1 + u!
1,j-1
5х5у 4 4 X 4 у ,
где 4 ^ - шаг по времени для процесса диффузии (вязкой подзадачи), а4 х и4 у - по пространству. Индекс 1 отвечает за изменение переменной вдоль оси х, )' - вдоль оси у. Остальные производные аппроксимируются аналогичным образом.
Граничные условия прилипания потока на поверхности пластины для слоя фиктивных ячеек j = 0, находящихся с «обратной стороны» пластины, сводятся к
и;0 =-и;1 v;0 = 0, (37)
c
где ячейка с координатами (1,1) является прилегающей к поверхности пластины.
В работах [8, 9] было показано, что при выборе шага по времени в виде
Д ^ = №2/ V (38)
ошибка численного решения уравнения диффузии с помощью схемы «донор-акцептор» и схемы «вперед по времени, центральная по пространству» (ВВЦП) зависит только от константы к и от количества сделанных шагов по времени. Здесь Ь - шаг ячейки расчетной сетки по пространству h = Д х = Д у . В работе [9] путем серии численных экспериментов было найдено значение к = 0.22, при котором ошибка решения минимальна. В работе [5] было показано, что данное значение коэффициента к справедливо и для нашей численной схемы.
ТЕСТИРОВАНИЕ ЧИСЛЕННОЙ СХЕМЫ
Рассматривается задача о продольном обтекании плоской пластины и формирования на ней пограничного слоя. Данная задача имеет точное в рамках теории пограничного слоя решение, полученное Блазиусом для ламинарного случая, которое хорошо подтверждается в эксперименте [9]. Известно, что с помощью введения безразмерной нормальной координаты
Л = У
(39)
и величины продольной компоненты скорости
u
*
u =
(40)
u
профили компоненты скорости u в разных сечениях пластины x = const будут совпадать, что является удобным для сравнения получаемых данных.
Сравнение профилей скорости с решением Блазиуса производилось для двух сечений на расстояниях хп » 0.5b и хп » 0.75b от переднего края пластины, где b - длина пластины.
Расчеты проводились в диапазоне чисел
Рейнольдса от 104 до 107. Результаты численного моделирования приведены на рисунках 1-8. Здесь ax и ay - размеры расчетной области вдоль координат х и y соответственно, nx и ny - число ячеек сетки вдоль координат х и у. Видно, что необходимый шаг расчетной сетки пропорционален квадратному корню из числа Рейнольдса. Это условие вытекает из обратной пропорциональности толщины ламинарного пограничного слоя корню квадратному из числа Рейнольдса.
В результате численного моделирования задачи было выявлено, что в диапазоне чисел Рей-
нольдса от 104 до 107 не происходит перехода ламинарного пограничного слоя в турбулентный. При этом никаких ограничений на численную схему не накладывалось. Этот вопрос требует дальнейших исследований. Также оказалось сюрпризом, что TVD схема первого порядка Рое [7] показывает лучшие результаты для ламинарного пограничного слоя по сравнению со схемой второго порядка [4].
СПИСОК ЛИТЕРАТУРЫ
1. Лойцянский, Л.Г. Механика жидкости и газа / Л.Г. Лойцянский. - М.- Л.: Гос. изд. технико-теоретической литературы. 1950. - 676 с.
2. Себиси, Т. Конвективный теплообмен / Т. Себиси, П. Брэдшоу. - Пер. с англ. С.С. Ченцова и В.А. Хохрякова. - Под. ред. Пирумова У.Г. - М.: Мир, 1987. - 593 с.
3. Harten, A. High resolution schemes for hyperbolic conservation laws / A. Harten // J. of Comp. Phys. - V. 49. - 1983. - P. 357-393.
4. Carofano, G.C. Blast computation using Harten's total variation diminishing scheme / Garry C. Carofano. Technical report ARLCB-TR-84029. US Army Armament Research and Development Center. 1984.
5. Никонов, В.В. О тестировании конечно-разностной схемы для моделирования процесса вязкой диффузии с учетом сжимаемости газа в двумерном случае / В.В. Никонов // Известия Самарского научного центра РАН. 2020. - Т. 22. - № 5. - С. 128-131.
6. Самарский, А. А. Аддитивные схемы для задач математической физики / А.А. Самарский, П.Н. Ва-бищевич. - М.: Наука, 2001. - 319 с.
7. Roe, P.L. Approximate Riemann solvers, parameter vectors, and difference schemes / P.L. Roe // J. of Comp. Phys. - V. 135. - 1997. - P. 250-258.
8. Nikonov, V. The Ratio between Spatial and Time Resolutions for the Diffusion Substep in 2D Computational Vortex Methods / V. Nikonov, N. Kornev, A. Leder // Schiffbauforschung. - 2002. - Vol. 41. - N 3/4. - Pp. 5-12.
9. Никонов, В.В. О выборе шага по времени в схеме ВВЦП при расчете процесса диффузии / В.В. Никонов // Сборник трудов 13-го Всероссийского семинара по управлению движением и навигации летательных аппаратов. - Самара: СГАУ, 2007. - Ч. 2. С. 55-57.
10. Шлихтинг, Г. Теория пограничного слоя / Г. Шлих-тинг // Пер. с нем. Г.А. Вольперта. Под. ред. Лой-цянского Л.Г. - М.: Наука. - 1974. - 712 с.
Известия Самарского научного центра Российской академии наук, т. 23, № 6, 2021 и*
1.2 1.0 0.8 0.6 04 0.2 0.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0
Рис. 1. Распределение продольной безразмерной компоненты скорости и* в пограничном слое над поверхностью пластины в сравнении с решением Блазиуса, схема [7]: Ь = 0.1 м, Яе = 104, ах хау = 0.2х0.14 м, пх хпу = 400х280; и„ = 14.63394141 м/с,
= 10132.5 па, = 0.1225 кг/м3, <И = 7.038175965-Ю-7 с, <ИЛ = 3.758386-Ю-4 с, численное решение: О - сечение X = 0.5Ь, I = 40; А - сечение X = 0.75Ь, I = 40; X - сечение X = 0.75Ь, I = ; — - решение Блазиуса
и*
1.2 1.0 0.8 0.6 0.4 0.2 0.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0
Рис. 2. Распределение продольной безразмерной компоненты скорости и* в пограничном слое над поверхностью пластины в сравнении с решением Блазиуса, схема [4]: Ь = 0.1 м, Яе = 104, ах хау = 0.2 х0.14 м, пх хпу = 400х280; и„ = 14.63394141 м/с,
= 10132.5 па, = 0.1225 кг/м3, <Ц = 7.03817596510-7 с, ^ = 3.758386Ю-4 с, численное решение: О - сечение X = 0.5Ь, I = 40 ; А - сечение X = 0.75Ь, I = 40dtd ; Х- сечение X = 0.75Ь, I = ;--решение Блазиуса
—.-«г „ж 1 К Ж> Жо Ж
/ / к / ч у
ь/
£ / У
/
4 ж о *
о о л ж О ^ ж 0&
X о ж 1
и*
■яг ж А ж
-- •
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0
Рис. 3. Распределение продольной безразмерной компоненты скорости и* в пограничном слое над поверхностью пластины в сравнении с решением Блазиуса, схема [7] :
b = 0.1 M, Re = 105,а Xа = 0.2X0.04 м, n Xп = 800Xl60;u = 14.63394141 м/с,
' ' X у ' X у 'го '
= 101325.0 Па, = 1.225 кг/м3, dt = 3.52436793510-7 с, dtd = 9.395965Ю-4 с,
численное решение:0 - сечение X = 0.5Ь, t = 10 dt, ;Л - сечение X =
X - сечение X = - решение Блазиуса
о & £ Й
О/
к/
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0
Рис. 4. Распределение продольной безразмерной компоненты скорости и* в пограничном слое над поверхностью пластины в сравнении с решением Блазиуса, схема [4], Ь = 0.1 м, Яе = 105, ах хау = 0.2 х0.04 м, пх хпу = 800х160; = 14.63394141 м/с,
= 101325.0 па, = 1.225 кг/м3, dt = 3.524367935-Ю-7 с, = 9.395965Ю-4 с, численное решение: О- сечение X = 0.5Ь, I = 10; Д- сечение X = 0.75Ь, I = 10; Х- сечение X = 0.75Ь, I = 20dtd ; — - решение Блазиуса
1.2 -
1.0 -
0.8 -
0.6 -
0.4 -
0.2 -
о.о -,
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0
Рис. 5. Распределение продольной безразмерной компоненты скорости и* в пограничном
слое над поверхностью пластины в сравнении с решением Блазиуса, схема [7]: Ь = 0.1 м, Яе = 106, ах хау = 0.2х0.02 м, пх хпу = 1600х160; = 146.3394141 м/с,
= 101325.0 Па, = 1.225 кг/м5, Л = 1.760862990-10-7 с, й^ = 2.348991-10-4 с, численное решение: О - сечение х = 0.5Ь, I = 5сЙл ; Л - сечение х = 0.75Ь, I = ; X - сечение х = 0.75Ь, I = 1; — - решение Блазиуса
""<>1 * о & &
/
/
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0
Рис. 6. Распределение продольной безразмерной компоненты скорости и* в пограничном слое над поверхностью пластины в сравнении с решением Блазиуса, схема [4]: Ь = 0.1 м, Яе = 106, ах хау = 0.2х0.02 м, пх хпу = 1600х160; и„ = 146.3394141 м/с,
= 101325.0 па, = 1.225 кг/м3, Л = 1.760862990-10-7 с, й^ = 2.348991-10-4 с, численное решение: О - сечение X = 0.5Ь, I = 5; Л - сечение X = 0.75Ь, I = 5; X - сечение X = 0.75Ь, I = 1; — - решение Блазиуса
V
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0
Рис. 7. Распределение продольной безразмерной компоненты скорости и* в пограничном слое над поверхностью пластины в сравнении с решением Блазиуса, схема [7]:
Ь = 1.0 М, Яе = 107, ахXау = 2.0X0.05 м, пхXпу = 4800Х120; = 146.3394141 м/с, = 101325.0 Па, = 1.225 кг/м3, dt = 5.875709711 -Ю-7 с, = 2.6099903-Ю-3 с, численное решение:0 - сечение X = 0.5Ь, I = 5;А - сечение X = 0.75Ь , I = ; Х- сечение X = 0.75Ь, I = 10dtd ; - решение Блазиуса
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0
Рис. 8. Распределение продольной безразмерной компоненты скорости и* в пограничном слое над поверхностью пластины в сравнении с решением Блазиуса, схема [4]: Ь = 1.0 м, Яе = 107, ах X ау = 2.0 X 0.05 м, пх X пу = 4800 Х120; и „ = 146.3394141 м/с,
= 101325.0 Па, = 1.225 кг/м3, dt = 5.875709711 -10-7 с, = 2.6099903-Ю-3 с, численное решение: О - сечение х = 0.5Ь, I = 5dtd ; А - сечение х = 0.75Ь, I = ; Х- сечение х = 0.75Ь, I = ; __ - решение Блазиуса
2D DIRECT NUMERICAL SIMULATION OF A COMPRESSIBLE BOUNDARY LAYER ON A FLAT PLATE
© 2021 V.V. Nikonov
Samara National Research University named after Academician S.P. Korolyov
In a paper a compressible air boundary layer on a flat plate is modeled. Total Variation Deminition (TVD) scheme is used for compressibility calculation. The longitudinal velocity profile was obtained for two sections perpendicular to the plate plane for different Reynolds numbers. The numerical results are compared with the Blasius' solution for the laminar case.
Keywords: air, compressibility, Navier-Stokes equations, boundary layer, direct simulation, numerical
simulation, two-dimensional flow.
DOI: 10.37313/1990-5378-2021-23-6-94-102
Valeriy Nikonov, PhD, Associate Professor of Aircraft Construction and Design Department. E-mail: v_nikonov@mail.ru
