CC BY
42
УДК 532.526
О ВЛИЯНИИ ДИАГОНАЛЬНЫХ КОМПОНЕНТОВ ТЕНЗОРА СКОРОСТЕЙ ДЕФОРМАЦИИ НА ЛАМИНАРНЫЙ ПОГРАНИЧНЫЙ СЛОЙ НА ПЛОСКОЙ ПЛАСТИНЕ
© 2009 В.В. Никонов, В.Г. Шахов
Самарский государственный аэрокосмический университет
Поступила в редакцию 07.05.2008
Рассматривается прямое численное моделирование ламинарного пограничного слоя на плоской пластине в широком диапазоне чисел Рейнольдса. Используется метод конечных разностей. Применяется схема интегрирования с разными шагами по времени из-за разности скоростей протекания процессов диффузии и конвекции. Исследуются два варианта постановки задачи с упрощенным и полным представлением члена, отвечающего за вязкую диффузию. Полученные профили скорости сравниваются с решением Блазиуса и результатами схемы метода расщепления завихренности (МРЗ). Отмечается, что, несмотря на одинаковый порядок точности решения, метод конечных разностей в 20 раз более медленный, чем МРЗ.
Ключевые слова: Пограничный слой, плоская пластина, численный метод, число Рейнольдса, диффузия, конвекция, тензор скоростей деформации, решение Блазиуса
1. МАТЕМАТИЧЕСКАЯ ФОРМУЛИРОВКА МЕТОДА
Математическая формулировка задачи описывается уравнениями Навье-Стокса. Давление в ламинарном пограничном слое на плоской пластине мало меняется, поэтому будем считать его постоянным. Рассматривается два варианта постановки задачи с упрощенным
д и д и д и д 2и
--+ и--+ V— = V—2
д! дх ду ду
д v д v д v д 2v
--+ и--+ v— = v—2
д t дх д y дх
(1)
и полным д и
ди ди
--+ и--+ v— = v
дt дх ду
{ д 2и д 2и ^ +
дх ду
д v дv дv
--+ и--+ v— = v
дt дх ду
{ д 2v дМ дх2 + ду2
(2)
ты, ! = Л/и ш - безразмерное время, где Ь -длина пластины, и ш - скорость набегающего потока. Упрощение системы (2) заключается в отбрасывании диагональных компонентов тензора скоростей деформации.
Для замыкания системы уравнений необходимо добавить уравнение неразрывности, которое для несжимаемой жидкости имеет вид
д и + д v = о дх ду
(3)
представлением члена, отвечающего за вязкую диффузию. Здесь используются безразмерные переменные: и = и*/и , v = v*/uш - компоненты скорости, v = v*/uшb = l/Re - безразмерный коэффициент кинематической вязкости, х = х*/Ь , у = у*/Ь - безразмерные координа-
Никонов Валерий Владимирович, кандидат технических наук инженер НТП "Авиатехнокон". E-mail: [email protected].
Шахов Валентин Гаврилович, кандидат технических наук, профессор, заведующий кафедрой аэрогидродинамики. E-mail: [email protected].
Переопределенность системы (1)-(3) в численном методе преодолевается следующим образом. Сначала решаются уравнения (1) или (2). Затем компонента V скорости корректируется так, чтобы поле скорости удовлетворяло уравнению неразрывности (3).
По аналогии с методом расщепления завихренности [1] изменение скорости в (1) или (2) вследствие конвекции и диффузии рассматривается отдельно, т.е. производится расщепление по физическим процессам.
На этапе конвекции решается следующая конечно-разностная система уравнений:
- v
Д tc
vt+41 - vt.
-Ы = -ut
Д tc
I.J
д v ,д х.
д v
- v-
'Л д у/ц
■ (4)
где д1с - шаг по времени для процесса конвекции, а для аппроксимации частных производных по
пространству применяется конечно-разностная схема "вперед по потоку" 2-го порядка
(дх).. = ( -4и1-и + ^)' иы/у >0и'-1/2]>0 (|)ц = 21(-^ ^0'и-1/2^0 (д^х) . = 21(-и,+Я] +3и,+Ц -3и"и + и'-2>])' и,+"Я] <0'Ц-1/2"] >0,
(дх).. =-2К( -4и'+ц+и'+2^)' и'+1/2^ <0^ -0 (5>
где величины компонент скорости и . и и 2.
определяются как
и
_ 'О "-'О ,, _ 'О
1+1/2^ " 2 ' 2^ " 2 ' (6)
а Ь - размер ячейки однородной расчетной сетки. Для ячеек вблизи плоской пластины и внешних границ расчетной области применялась конечно-разностная схема "вперед по потоку" 1-го порядка
ди
\д Х1 ди
* д х.
и'+у - иу
и ; > 0
и ; < 0'
(7)
и'+ 4 4 - и', ( и
-— = v
4 Ъ
■у+1
V
г+Л1 г Уц 4 - УЦ
д!а
' у!-,- 4+ ^ Л
(9)
где - шаг по времени, необходимый для про-
цесса диффузии. Для системы (1) схема (9) соответствующим образом упрощается.
Граничные условия (ГУ) прилипания пото-
ка на поверхности пластины на этапе диффузии
Граничные условия (ГУ) прилипания потока на поверхности пластины на этапе конвекции для слоя фиктивных ячеек .-1, находящихся с "обратной стороны" пластины, сводятся к
и,н = 0, (8)
где ячейка с координатами (1, . является прилегающей к поверхности пластины.
На этапе диффузии используется схема ВВЦП [2], которая для системы (2) будет иметь следующий вид
и, 2и',,+ иЦ, и, - 2и',,+ и' Л
для слоя фиктивных ячеек .-1, находящихся с "обратной стороны" пластины, сводятся к
и',Н =-ич , (10)
где ячейка с координатами (1, . является прилегающей к поверхности пластины.
В работе [3] было показано, что для достижения заданной точности шаг по времени для схемы ВВЦП определяется следующим соотношением
д! = М2/V, (11)
где ка - константа, от которой зависит ошибка численного решения. Минимальная ошибка получается при ка = 0,22.
Шаг по времени для процесса конвекции удовлетворяет неравенству
4!с - 0.111, (12)
для критерия Куранта-Фридриха-Леви [4], равного 0,1, который получен для схемы (4), (5) серией численных расчетов. При этом он не обязательно совпадет с оптимальным шагом для расчета процесса диффузии. По этой причине предлагается применять метод интегрирования с разными шагами по времени (ИРШ) для процессов диффузии и конвекции.
Корректировка поля скорости для удовлетворения уравнения неразрывности (3) производится на этапе конвекции перед расчетом (4). Для каждой ячейки сетки вниз по потоку профиль скорости будет менее заполненным, чем для предыдущей ячейки, вследствие развития пограничного слоя. Это приводит к тому, что в каждую ячейку сетки в продольном направлении втекает больше жидкости, чем вытекает из нее. Поэтому предполагается, что избыточное количество жидкости вытекает в поперечном направлении, вызывая тем самым появление вертикальной компоненты скорости течения.
Из (3) можно записать
д у = ди
ду д х . (13)
Тогда, согласно (13) и ГУ непротекания на поверхности пластины (у = 0), вертикальная компонента скорости для ячеек, прилегающих к пластине сверху, находится следующим образом:
У1, =-0.51 / ^
4 Л дх
(14)
Для остальных ячеек, расположенных над пластиной, вертикальная компонента скорости определяется как
\ \ Л
У',, = У',, -1 - 05Ьу
. (15)
, -1 У
Для аппроксимации пространственных производных в правых частях выражений (14), (15), как и в [1], использовалась центральная конечно-разностная схема. Формулы для ячеек, расположенных с нижней стороны пластины, получаются зеркальным отражением (14), (15) относительно оси абсцисс ОХ.
2. ПРЯМОЕ ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЛАМИНАРНОГО ПОГРАНИЧНОГО СЛОЯ НА ПЛОСКОЙ ПЛАСТИНЕ
Рассматривается задача о продольном обтекании плоской пластины для того, чтобы проверить диффузию с поверхности тела совместно с моделированием процессов диффузии в свободном потоке и конвекции. Данная задача имеет точное в рамках теории пограничного слоя решение, полученное Блазиусом для ламинарного случая, которое хорошо подтверждается в эксперименте [5]. Известно, что с помощью введения безразмерной нормальной координаты
П = У
u„
* * v x
(16)
и величины вертикальной компоненты скорости
v u ^
u ^ x
v
(17)
профили компонент скорости u и vf в разных сечениях пластины x = const будут совпадать, что является удобным для сравнения получаемых данных.
Сравнение профилей скорости с решением Блазиуса производилось для четырех сечений на расстояниях x « 0.25, x « 0.5, x « 0.75 и x « 0.9
L п 7 п 7 п п
от переднего края пластины.
Расчеты проводились в диапазоне чисел Рейнольдса от 10 до 10 5. Параметры численного мо-
делирования приведены в табл. 1. Здесь te - конечное безразмерное время счета, Tm - затраты машинного времени в часах, минутах и секундах для решения систем уравнений (1) и (2), L и H -размеры расчетной области вдоль координат x и y соответственно, n и n - число ячеек сетки
J x y
вдоль координат x и y.
Видно, что необходимый шаг расчетной сетки пропорционален квадратному корню из безразмерной кинематической вязкости n. Это условие вытекает из обратной пропорциональности толщины ламинарного пограничного слоя корню квадратному из числа Рейнольдса. Также следует отметить существенный рост потребного машинного времени при увеличении числа Рейнольдса (табл. 1), так для числа Re = 10 5 расчеты длились более 2 суток. Заметим, что результаты данной работы получены на компьютере с частотой процессора 3.2 ГГц и 1 Гб оперативной памяти, который является в 2.6 более производительным, чем компьютер, который использовался в работе [1]. Сравнивая время счета для Re = 1000 и h = 0,01, приведенное в табл. 1, со временем счета для МРЗ в работе [1], можно заключить, что рассматриваемый в данной работе метод является в 20 раз более медленным, чем МРЗ. Это происходит потому, что на этапе конвекции в рассматриваемом методе приходится использовать в 15 раз более мелкий шаг (12), чем в схеме МРЗ.
В результате моделирования задачи было выявлено, что в диапазоне чисел Рейнольдса от 100 до 10 5 для систем (1) и (2) получаются одинаковые результаты (таблица 2), за исключением профиля вертикальной компоненты скорости для сечения хп « 0.9 при Re = 100 (см. рис. 1). Заметим, что для вертикальной компоненты скорости удается достичь достаточной для инженерных целей точности лишь на более мелких сетках (табл. 2). При числе Re = 10 для системы (1),
Таблица 1. Параметры численного моделирования пограничного слоя в зависимости от числа Рейнольдса
vf =
Re h LxH nx X ny Atd ¿tC te Tm, [чч:мм:сс]
(1) (2)
10 0.1 4x8 40x80 0.022 7.333333-10"3 16.8 00:00:08 00:00:12
0.05 4x8 80x160 5.5-10"3 2.75-10"3 16.8 00:01:28 00:02:22
0.025 4x8 160x320 1.375-10 3 1.375-10 3 16.8 00:13:57 00:15:32
100 0.03125 4x4 128x128 0.021484 3.069196-10"3 16.8 00:01:36 00:01:38
0.015625 4x4 256x256 5.371094-10"3 1.342773-10"3 16.8 00:32:38 00:37:20
103 0.01 4x2 400x200 0.022 0.001 16.8 00:25:12 00:29:34
0.005 4x2 800x400 5.5-10"3 5.0-10"4 8.4 02:51:08 02:54:33
104 0.003125 2x0.5 640x160 0.021484375 3.069196-10"4 8.4 01:24:47 01:35:28
0.0015625 2x0.5 1280x320 5.371094-10"3 1.534598-10"4 8.4 21:30:10 21:29:51
105 0.001 2x0.1 2000x100 0.022 1.0-10"4 4.2 06:47:54 06:45:09
5-10"4 2x0.1 4000x200 5.5-10"3 5.0-10"5 4.2 61:49:26 64:02:59
Таблица 2. Средняя ошибка численного решения при моделировании пограничного слоя на плоской пластине для профилей скорости в центральном сечении х « 0.5
Яе ь 5и
МРЗ (1) (2) МРЗ (1) (2)
10 0.1 0.04 0.04 0.50 0.14 0.12 0.86
0.05 0.02 0.02 0.38 0.06 0.06 0.74
0.025 0.01 0.01 0.34 0.05 0.05 0.72
100 0.03125 0.02 0.02 0.02 0.11 0.12 0.10
0.015625 0.01 0.006 0.01 0.01 0.03 0.03
103 0.01 0.01 0.01 0.009 0.10 0.10 0.10
0.005 0.007 0.001 0.0008 0.02 0.01 0.01
104 0.003125 0.01 0.007 0.007 0.09 0.08 0.08
0.0015625 0.006 0.004 0.004 0.02 0.01 0.01
105 0.001 0.03 0.006 0.006 0.13 0.20 0.17
5-10"4 0.009 0.005 0.005 0.03 0.01 0.01
Рис. 1. Профиль вертикальной скорости на плоской пластине, Ие = 100, = 5.371094-10-3, 1 = 1.342773-10-3, 1 = 16.8, Ь = 0.015625, ЬхН = 4x4;
' 4 а ' 4 с ' ' ' '
_- профиль Блазиуса [5],
настоящая работа: ^ - хп = 0.25, + - хп = 0.5, □ - хп = 0.75, х - хп = 0.9
Рис. 2. Профиль продольной скорости на плоской пластине, Ие = 10, = 1 = 1.375-10-3, 1 = 16.8, Ь = 0.025, ЬхН = 4х8;
4 а 4 с
_- профиль Блазиуса [5],
настоящая работа: ф - хп = 0.25, + - хп = 0.5, □ - хп = 0.75, х - хп = 0.9
как и для МРЗ в [1] (МРЗ просто другой метод решения системы (1)), удается получить с точностью, достаточной для инженерных целей,
профиль Блазиуса, см. рис. 2.
Погрешность решения, приведенная в табл. 2, определялась следующим образом:
8U = max
U 0<r|<6
U ; - U
Ue
8v = max
0 <n<6
v,j- ve
V
где иех и Уех - поле скорости точного решения Блазиуса.
Для системы (2) при Ие = 10 наблюдается сильное отличие численного решения от решения Блазиуса (рис. 2), что объясняется влиянием диагональных компонентов тензора скоростей деформации. Для окончательного вывода не хватает данных эксперимента о продольном обтекании плоской пластины для данного числа Рейнольдса.
СПИСОК ЛИТЕРАТУРЫ
■ (17) 1. Никонов В.В., Шахов В.Г, Схема расчета скорости для метода "вихрь в ячейке" применительно к моделированию двумерного ламинарного пограничного слоя // Известия СНЦ РАН. 2005. Т. 7. № 2. С. 392-398.
2. Флетчер К. Вычислительные методы в динамике жидкостей / Пер. с англ. М.: Мир. Т. 1, 1991. 504 с.
3. Никонов В.В. О выборе шага по времени в схеме ВВЦП при расчете процесса диффузии // Сборник трудов 13-го Всеросс. семинара по управлению движением и навигации ЛА. Самара: СГАУ, 2007. Ч. 2. С. 55-57.
4. Ferziger J., Peric M. Computational methods for fluid dynamics // 3 rev. ed., Springer-Verlag. 2002. 423 p.
5. Шлихтинг Г. Теория пограничного слоя / Пер. с нем. Г.А. Вольперта, под. ред. Л.Г.Лойцянского. М.: Наука, 1974. 712 с.
ABOUT INFLUENCE OF DIAGONAL COMPONENTS OF A STRAIN RATE TENSOR IN LAMINAR INCOMPRESSIBLE BOUNDARY LAYER ON A FLAT PLATE
© 2009 V.V. Nikonov, V.G. Shakhov
Samara State Aerospace University,
Direct numerical simulation of the laminar boundary layer on the flat plate in the wide range of Reynolds numbers is considered. Finite difference method is used. The Multi-Time-Step Scheme is used in the numerical algorithm because the differences in the rates of convection and diffusion processes. Two variants of the problem with the simplified and full representation of the viscous term are investigated. Obtained profiles are compared with Blasius solution and results of the Vorticity Splitting Method (VSM). It should be noted, that the finite difference method is slow in 20 times than VSM despite of same solution accuracy. Key words: boundary layer, flat plate, numerical method, Reynolds number, diffusion, convection, strain rate tensor, Blasius solution
Valery Nikonov, Candidate of Technics, engineer NTP "Aviatechnokon". E-mail: [email protected]. Valentin Shakhov, Candidate of Technics, Professor, Head of Aerohydrodynamic Department. E-mail: [email protected].
