CC BY
23
УДК 532.526
МЕХАНИКА И МАШИНОСТРОЕНИЕ
СХЕМА РАСЧЕТА СКОРОСТИ ДЛЯ МЕТОДА "ВИХРЬ В ЯЧЕЙКЕ" ПРИМЕНИТЕЛЬНО К МОДЕЛИРОВАНИЮ ДВУМЕРНОГО ЛАМИНАРНОГО ПОГРАНИЧНОГО СЛОЯ
© 2005 В.В. Никонов, В.Г. Шахов Самарский государственный аэрокосмический университет
Предлагается новая схема вычисления скорости в методе "вихрь в ячейке". На каждом шаге по времени процессы конвекции, диффузии в свободном потоке и с поверхности обтекаемого тела рассматриваются отдельно. Основная идея предлагаемой схемы, названной как метод расщепления завихренности (МРЗ), заключается в рассмотрении независимо друг от друга частных производных скорости, входящих в определение завихренности. Кроме того, применяется интегрирование с разными шагами по времени по причине разности скоростей протекания процессов диффузии и конвекции. В качестве тестовой задачи используется прямое численное моделирование ламинарного пограничного слоя на плоской пластине в широком диапазоне чисел Рейнольдса. Профили скорости, полученные с помощью новой схемы, сравниваются с решением Блазиуса и результатами других авторов. Показано, что схема МРЗ позволяет получать результаты с достаточной точностью в широком диапазоне изменения чисел Рейнольдса.
Математическая формулировка метода "вихрь в ячейке"
Метод "вихрь в ячейке" (ВЯ) часто используется для расчета отрывных течений в рамках модели вязкой несжимаемой жидкости. Данный метод относится группе методов дробных шагов для решения уравнения Навье-Стокса в безразмерных переменных завихренность-скорость, которое для двумерного случая имеет вид
— + (и • У)ш = уАш д I
где ш - завихренность
ш = ■
д V д и
д х д у
t - время, и - вектор скорости, V - безразмерная кинематическая вязкость v=1/Re, где Яе - число Рейнольдса, у - оператор Гамильтона, А - оператор Лапласа.
В численной схеме метода ВЯ вместо завихренности используется вихревая интенсивность, по размерности совпадающая с циркуляцией
х1 +h /2 Уj +^2
= | | ш (х,у)^у. (3)
х1 ^/2 у] -Ы2
Поле скорости в вихревых методах представляется в виде суммы постоянной и ш и
соленоидальной и у составляющих
и = и ш + и у, (4)
где и ш - скорость набегающего потока, а и у
определяется как сумма полей скорости, индуцированных потенциальными вихрями (закон Био-Савара)
(1)
(2)
и (х) ^ X (X -У
и У (Х) Л Ъ
2 п j=l 1:
(х-У2
. (5)
Число необходимых операций при применении схемы (5) составляет порядок N 2, где N = пх пу. Для сокращения времени вычисления скорости по (5) применяется кластеризация и метод мультипольного разложения [2].
При решении уравнения (1) методом ВЯ процессы конвекции, диффузии завихренности в свободном потоке и с поверхности тела на каждом шаге по времени рассматриваются отдельно [1, 4].
Схема метода расщепления завихренности
При разработке данной схемы мы исходили из того, что для заданной величины за-
вихренности может существовать бесконеч-
д V д и
ное количество пар значений дХ и ду , при
которых равенство (2) будет выполняться. По этой причине в методе расщепления завихренности (МРЗ) указанные частные производные рассматриваются отдельно и вместо циркуляции (3) в численной схеме используются величины
х; +И/2 Уj +И/2
А иу(Х1,^) = / /
х; -h/2 yj-И/2 х; + h/2 у j+h/2
^(ЗД) = / /
х; - h/2 у j - h/2
ди(х,у)
д у
д v(x,y)
дХ
dxdy
dxdy
. (6)
Вычисление скорости
и моделирование процесса конвекции
при использовании схемы МРЗ
Для определения скорости данный метод использует вместо (5) схему интегрирования вдоль координатных линий, при этом двумерная задача сводится к нескольким одномерным. Если известна величина скорости с одной стороны ячейки, а также величины изменения скорости при проходе через ячейку в данном направлении (6), то скорость с противоположной стороны ячейки определяется как
и(У + 0,5) = и(У - 0,5) + Аиу(У)^х, V(I + 0,5,я) = V(I-0.5,:)е + Аvx(i,j)/hy . (7)
Скорость в центрах ячеек сетки рассчитывается следующим образом:
и(У) = и(У - 0,5) + 0,5А uy(I,j)/hx,
v(I,j) = V(I - 0,5, Я) + 0,5Ата(У)/Иу. (8)
Формулы (7) и (8) используются при интегрировании в положительном направлении осей ОХ и ОУ. Интегрирование в обратном направлении проводится аналогичным образом с учетом смены знака. Число необходимых операций для схемы МРЗ имеет порядок
Из граничных условий (ГУ) задач внешнего обтекания тел следует, что величина скорости на внешней границе (невозмущенном
потоке) равна и ш, тогда скорость внутри вяз-
кого течения удобно находить с помощью (7)-(8) интегрированием от внешней границы течения до другой границы течения или поверхности тела.
В схеме метода МРЗ частицы движутся с потоком и переносят величины А и А .
Г УХ иу
После расчета поля скоростей течения новые координаты частиц получаются аналогично методу ВЯ численным интегрированием системы обыкновенных дифференциальных уравнений методом Эйлера. Новое местоположение частиц не обязательно совпадет с координатами расчетной сетки, поэтому процедура перераспределения используется отдельно для величин А и А :
¥Х иу
Аиу(х^ Уя) = Аиу(хк,у1)Л(х1 -Хк)Л(Уя -у1),
А vx(Хi,Уj) = А ух(хк,у1)Л (х - хк) Л(у - у^ . (9) В качестве функции перераспределения Л в (9) используется формула "облако в ячейке" [8, 9]
Л (7) =
)1 - * I 0,
0 < 7* < 1
* 1
7 > 1
здесь 7* = |*|/Ь.
Шаг по времени для процесса конвекции определялся с помощью неравенства
И
и ^ (10)
А^ < 1.5-
соответствующего критерию Курант-Фрид-рих-Леви [10] с величиной коэффициента 1,5. При этом он не обязательно совпадет с оптимальным шагом для расчета процесса диффузии. По этой причине предлагается применять метод интегрирования с раздельными шагами по времени для процессов диффузии и конвекции.
Диффузия в свободном течении
Диффузия в свободном потоке для схемы МРЗ рассчитывается с использованием метода Д-А аналогично [1, 5], но отдельно для А и А ,
^ иу ¥Х5
А иуО + А ^У^ =А иу^^У^ +
+ X (Диу^у^'^у)-Аиу^Х^^у)) ,
А vx(t + А t,ХI,УJ) = А „Д^Уя) +
+ X (A„t(t,xq;yJ)Gi;(x) - Avx (t,xq,yj)G*(x)) ,
i- nd < q<i + nd
(11)
где коэффициенты G*jq определяются по формуле
Gpq(z) =
erf
Ч + h/2 - Zq^
>/4vAt
-erf
( z - h/2 - z ^
Здесь z - xp или yp - координаты центра ячейки p-го вихря, erf - интеграл вероятности (функция ошибок) [6]. В работе [7] было показано, что для достижения заданной точности шаг по времени для метода Д-А определяется следующим соотношением:
A t = kdh2/ V , (12)
где kd - константа, зависящая только от радиуса nd "диффузионной молекулы", и для nd = 1 удовлетворяет условию 0,2 < kd < 0,21. Здесь принимается kd = 0,21, так как в данном случае [4] ошибки методов Д-А и моделирования процесса конвекции будут иметь разные знаки и компенсировать друг друга.
Диффузия с поверхности обтекаемого тела
В схеме МРЗ с поверхности тела вместо завихренности диффундируют величины Auy и Avx. Для удобства сначала определим Aun и Avs в системе координат, связанной с панелью, а затем пересчитаем их в глобальную систему координат (СК). Вследствие ГУ непротекания и прилипания жидкости получается
A un(Si,nj) = u(Si,+0)h
(
erf
n. + h/2
Л
V4VJ
(
- erf
n. - h/2
V4VAt
A vs(Si,nj) = 0,
(13)
где величина скорости жидкости у поверхности тела определяется следующим образом:
и(^, +0) = limu(si,n) (14)
n^+0
ности тела и самой поверхности тела. Однако скорость последнего равна нулю так как рассматривается обращенное движение.
Пересчет в глобальную СК выполняется по формулам
А иу(Х^У,) = Лип^,п,)Уп,
Ата(х^) = Лип^)Хп, (15)
где п = {хп,уп} - единичный вектор, нормальный к рассматриваемой панели, записанный в глобальной СК, (si, п,) и (xi, у,) - координаты центра рассматриваемой ячейки в локальной и глобальной СК соответственно.
В случае произвольного контура тела координаты (х, у,) в (15) не обязательно совпадут с ячейками глобальной сетки, по этой причине для значений DDuy и DDvx применяется процедура перераспределения (9).
Коррекция поля скорости для удовлетворения уравнению неразрывности
В случае применения МРЗ в схему метода необходимо включить корректировку поля скорости для учета уравнения неразрывности. Заметим, что в обычной схеме метода ВЯ уравнение неразрывности выполняется автоматически вследствие локально вращательного представления движения. Для каждой жидкой частицы вниз по потоку профиль скорости будет менее заполненным, чем для предыдущей частицы, вследствие развития пограничного слоя. Это приводит к тому, что в каждую жидкую ячейку в продольном направлении втекает больше жидкости, чем вытекает из нее. Поэтому предпагается, что избыточное количество жидкости вытекает в поперечном направлении, вызывая тем самым появление вертикальной компоненты скорости течения.
Следуя уравнению неразрывности для двумерного случая, можно записать
здесь И - размер ячейки локальной сетки в направлении координаты б. В общем случае И приблизительно равен шагу сетки И, так как длина панели 1 (на которые разбивается тело) не всегда кратна И. В выражении (13) вместо величины скорости (14) должна стоять разность скоростей жидкости у поверх-
д v д y
д u
д x '
(16)
Тогда, согласно (16) и ГУ непротекания на поверхности пластины (V = 0), вертикальный компонент скорости для ячеек, прилегающих к пластине сверху, находится следую-
щим образом:
j-W fx
(17)
Для остальных ячеек, расположенных над пластиной, вертикальная компонента скорости определяется как
(/ „ \ / „ \ >
. (18)
Vi,j= Vi,j-i - 0,5hy
fu
,f x.
/_fu ^
\f x/
V 1 ' 1,j * ' ч-1 у
Для аппроксимации пространственных производных в правых частях выражений (17)-(18) использовалась центральная конечно-разностная (ЦКР) схема.
Предварительные расчеты обтекания плоской пластины показали, что ЦКР схема позволяет получать хорошие результаты во всем рассматриваемом диапазоне чисел Рей-нольдса (10 < Яе < 10 6), в отличие от схемы вперед по потоку 1-го порядка. Конечно-разностная (КР) схема вперед по потоку 2-го порядка имеет одинаковый с ЦКР порядок погрешности, но является более трудоемкой и менее компактной, так как является пятиточечной схемой.
Формулы для ячеек, расположенных с нижней стороны пластины, получаются зеркальным отражением (17), (18) относительно оси абсцисс ОХ.
Алгоритм метода МРЗ
Краткая форма алгоритма схемы МРЗ имеет следующий вид:
1) Определение величины разности скорости на поверхности тела с помощью (7).
2) Расчет диффузии в свободном потоке (11).
3) Расчет диффузии с поверхности тела
(13), (15).
4) Расчет поля скорости по (8).
5) Коррекция поля скорости согласно
(17), (18).
6) Конвекция жидких частиц и перераспределение их в ячейки сетки (9).
7) Сохранение результатов расчета и переход к следующему расчетному шагу.
Данный алгоритм записан для случая равенства шагов по времени = А^ для расчета процессов диффузии и конвекции. Если из условий (12) и (10) получается, что
^ А^, то применяется схема интегрирования с разными шагами по времени. Шаги 1 - 3 алгоритма выполняются, если текущее расчетное время достигает величины t = ^ + где ^ - предыдущий момент времени расчета диффузии. Аналогично для процесса конвекции шаги 4 - 6 алгоритма выполняются, если текущее расчетное время достигает величины t = ^ + А^, где ^ - предыдущий момент времени расчета конвекции. Сохранение результатов (шаг 7 алгоритма) выполняется через указанное пользователем число расчетных шагов в файлах данных.
Прямое численное моделирование ламинарного пограничного слоя на плоской пластине
Рассматривается задача о продольном обтекании плоской пластины для того, чтобы проверить диффузию с поверхности тела совместно с моделированием процессов диффузии в свободном потоке и конвекции в методах ВЯ и МРЗ. Данная задача имеет приближенное решение, полученное Блазиусом для ламинарного случая, которое хорошо подтверждается в эксперименте [11]. Известно, что с помощью введения безразмерной нормальной координаты
Л = У.
и величины вертикального компонента скорости
v u^
V
(19)
профили компонентов скорости u и vf в разных сечениях пластины x = const будут совпадать, что является удобным для сравнения получаемых данных.
Сравнение профилей скорости с решением Блазиуса производилось для четырех сечений на расстояниях x1 0,25, X 0,5, х 0,75 и x4 0,9 от переднего края пластины.
С помощью предлагаемой здесь схемы МРЗ удалось получить хорошие результаты в диапазоне чисел Рейнольдса от 10 до 10 6. Параметры численного моделирования при-
u „ х
vf =
ведены в табл. 1. Здесь t - конечное безразмерное время счета, Tm - затраты машинного времени. Видно, что необходимый шаг расчетной сетки пропорционален квадратному корню из безразмерной кинематической вязкости. Это условие вытекает из обратной пропорциональности толщины ламинарного пограничного слоя корню квадратному из числа Рейнольдса. Также следует отметить существенный рост потребного машинного времени при увеличении числа Рейнольдса (табл. 1), так для числа Re = 10 6 расчет для схемы МРЗ длился более 5 суток. Для метода ВЯ такой расчет занял бы в log 2(N)/2 раз больше времени, что не представляет практичес-
кого интереса.
Для проверки сходимости метода МРЗ были проведены расчеты на вдвое более мелких сетках.
Для числа Re = 10 3 необходимый шаг сетки равнялся И = 0,01. Остальные параметры моделирования для удобства сравнения результатов были выбраны равными параметрам расчета с помощью вихревых методов.
Результаты прямого нестационарного моделирования пограничного слоя несжимаемой жидкости в сравнении с аналитическим решением Блазиуса для двух моментов времени приводятся на рис. 1-2.
При сравнении полученного решения с
Таблица 1. Параметры численного моделирования пограничного слоя в зависимости от числа Рейнольдса
Re h L x H nx x Пу Atd Atc t [чч:мм:сс]
10 0,1 4x8 40x80 0,021 0,105 16,8 00:00:05
0,05 4x8 80x160 0,00525 0,0525 16,8 00:00:51
0,025 4x8 160x320 0,0013125 0,0328125 16,8 00:12:12
100 0,03125 4x8 128x256 0,0205078125 0,0205078125 16,8 00:01:03
0,015625 4x4 256x256 0,0051269531 0,0205078125 16,8 00:08:28
10 3 0,01 4x2 400x200 0,021 0,021 16,8 00:03:11
0,005 4x2 800x400 0,00525 0,00525 16,8 00:45:18
10 4 0,003125 4x1 1280x320 0,0205078125 0,0041015625 16.8 01:35:15
0,0015625 2x0.5 1280x320 0,005126953125 0,001708984375 16,8 04:49:06
10 5 0,001 2x0.5 2000x500 0,021 0,0015 8,4 02:23:40
5-10-4 2x0.1 4000x200 0,00525 7,5-10-4 8,4 10:35:53
10 6 3,125-10-4 2x0.1 6400x320 0,0205078125 4,660867-10-4 8,4 36:57:44
1,563-10-4 2x0.1 12800x640 0,005126953125 2,330433-10-4 4,2 134:51:09
u
1
0.8 0.6 0.4 0.2
0
v
0 1 2 3 4 5
6
Л
Рис. 1. Профиль продольной скорости ламинарного пограничного слоя в сравнении с результатами других авторов (Re = 10 3, h = 0,01, t = 16,8):
_- профиль Блазиуса [11], - Wu [12]
(x = 0,5)______.- Ota [13] (x = 0,9),
.......- Nakamura [3] (x = 0,5),
настоящая работа: - x = 0,5, x - x = 0,9
f 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0 1 2 3 4 5
6
Л
Рис. 2. Распределениевертикальной компоненты скорости V внутри пограничного слоя в сравнении с аналитическим решением Блазиуса ^е = 10 3, 11 = 0,01, t = 16,8): численное решение: ^ - х = 0,25, + - х = 0,5, □ - х = 0,75, х - х = 0,9;_- Блазиус [11]
данными других авторов: [3, 12-13] (рис. 1), видно, что наилучшие результаты дает предлагаемая в данной работе схема МРЗ. Для профиля вертикальной скорости vf (19) схема МРЗ позволяет также получать хорошие результаты (рис. 2). Заметим, что данные о полученных значениях вертикальной скорости в работах [3, 12-13] не приводятся.
Оценка точности полученных с помощью схемы МРЗ результатов для профилей продольного и вертикального компонентов скорости в сечении, проходящем через центр пластины, приводится в табл. 2. Здесь 5 u и 5v - величины относительной погрешности продольного и вертикального компонентов скорости соответственно. Видно, что при необходимом размере ячейки сетки h = 3,125-10-4 для числа Рейнольдса 10 6 не удалось получить удовлетворительных результатов для профиля вертикальной скорости v. Однако на вдвое более меткой сетке удалось достичь точности, достаточной для инженерных целей. Интервал времени моделирования был сокращен с t = 8,4 до 4,2 из-за большой трудоемкости вычислений. Предыдущим расчетом на более грубой сетке было показано, что к этому моменту времени решение по всем рассматриваемым сечениям х = const принимает установившийся характер. Полученные профили скорости приводятся на рис. 3-4. Заметим, что данные для чисел Рейнольдса Re > 10 3 в работах [3, 12-13] не приводятся.
Следует отметить, что точность определения продольного компонента скорости u
Таблица 2. Ошибка численного решения при моделировании пограничного слоя на плоской пластине методом МРЗ для профилей скорости в центральном сечении x 0,5
U
Re h t 5u 5v
10 0,1 16,8 0,04 0,14
0,05 16,8 0,02 0,06
0,025 16,8 0,01 0,05
100 0,03125 16,8 0,02 0,11
0,015625 16,8 0,01 0,01
10 3 0,01 16,8 0,01 0,10
0,005 16,8 0,007 0,02
10 4 0,003125 16,8 0,01 0,09
0,0015625 16,8 0,006 0,02
10 5 0,001 8,4 0,03 0,13
5-10-4 8,4 0,009 0,03
10 6 3,125-10-4 8,4 0,05 0,62
1,563-10-4 4,2 0,01 0,04
1
0.8 0.6 0.4 0.2 0
0 1 2 3 4 5
6
Л
Рис. 3. Распределение продольной компоненты скорости и внутри пограничного слоя в сравнении с аналитическим решением Блазиуса ^е = 10 6, 11 = 1,5625ЧЧ10-4, t = 4,2): численное решение: ф - х = 0.25, + - х = 0,5, □ - х = 0,75, X - х = 0,9;_- Блазиус [11]
Vf
f 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0 1 2 3 4 5 6
Л
Рис. 4. Распределение вертикальной компоненты скорости V внутри пограничного слоя в сравнении с аналитическим решением Блазиуса ^е = 10 6, 1 = 1,5625 10-4, t = 4,2): численное решение: ^ - х = 0,25, + - х = 0,5, □ - х = 0,75, х - х = 0,9;_- Блазиус [11]
выше, чем V. Это объясняется тем, что для определения компонента V используется конечно-разностная схема, а для компонента и применяется интегральный метод.
Результаты численного моделирования показали, что разработанная численная схема метода с расщеплением завихренности позволяет получать хорошие результаты в диапазоне чисел Рейнольдса 10 < Яе < 10 6 в отличие от вихревых методов. Кроме того, предложенная схема МРЗ является значительно более быстрой (порядок числа операций равен 2^, чем самая быстрая реализация
метода ВЯ (порядок числа вычислений
N log 2(N)).
СПИСОК ЛИТЕРАТУРЫ
1. Taranov A., Kornev N., Leder A. Development of the Computational Vortex Method for Calculation of Two-Dimensional Ship Sections with Flow Separation. Schiflbauforschung, 39, 2000, 2.
2. Kornev N., Leder A., Mazaev K. Comparison of two fast algorithms for the calculation of flow velocities induced by a three-dimensional vortex field. Schiflbauforschung, 40, 2001, 1.
3. Nakamura H., Kamemoto K. Igarashi T. Analysis of unsteady heat transfer in the wake behind a circular cylinder in a uniform flow by a vortex and heat element method// Proceedings of the Second International Conference on Vortex Methods. Turkey, 2001.
4. Никонов В.В., Шахов В.Г. Модификация схемы "донор-акцептор" для расчета диффузии завихренности и ее применение в методе "вихрь в ячейке"// Вестник СГАУ 2003. № 1 (3).
5. Basin M., Kornev N. Beruecksichtigung der Reibung in der Wirbelmethode. ZAMM, 78, 1998, 5.
6. Справочник по прикладной статистике в 2-х т./ Под ред. Э. Ллойда, У. Ледермана. М.: Финансы и статистика, 1989.
7. Nikonov V., Kornev N., Leder A. The Ratio between Spatial and Time Resolutions for the Diffusion Substep in 2D Computational Vortex Methods. Schiflbauforschung, 2002, vol. 41, N 3/4.
8. Cottet G.-H., Koumoutsakos P. Vortex methods: theory and practice. Cambridge University Press, 2000.
9. Григорьев Ю.Н., Вшивков В.А. Численные
методы "частицы-в-ячейках", Новосибирск: Наука, Сибирская издательская фирма РАН, 2000.
10. Ferziger J., Peric M. Computational methods for fluid dynamics // 3 rev. ed., SpringerVerlag. 2002.
11. Шлихтинг Г., Теория пограничного слоя, пер. с нем. Г.А. Вольперта, под. ред. Л.Г. Лойцянского. М.: Наука, 1974.
12. Wu J.C. Numerical boundary conditions for viscous flow problems // AIAA Journal, 14, 1976.
13. Ota S., Kamemoto K. Study on higher resolution of vorticity layer over a solid boundary for vortex methods // Proc. of The Second Intern. Conf on Vortex Methods, Istanbul, Turkey, 2001.
VELOCITY CALCULATION SCHEME IN THE 2D COMPUTATIONAL VORTEX-IN-CELL METHOD APPLIED TO THE SIMULATION OF THE LAMINAR BOUNDARY LAYER
© 2005 VV. Nikonov, V.G. Shakhov
Samara State Aerospace University
The new velocity calculation scheme in the Vortex-in-Cell (VC) method is proposed. Within each time step three different kinds of mechanisms: convection, diffusion and vorticity creation at boundaries are treated separately. The proposed scheme named as a Vorticity Splitting Method (VSM) and its main idea consist in separately consideration of velocities partial derivatives, which are included in the vorticity. Moreover Multi-Time-Step simulation is used in the numerical algorithm because the differences in the rates of convection and diffusion processes. Direct Numerical Simulation of the laminar boundary layer on a flat plate in a wide range of Reynolds numbers is used for a test purpose. Velocity profiles obtained with the ordinary and new schemes are compared to Blasius solution and results of other authors. VSM scheme is allowed results obtaining with sufficient accuracy in the wide range of Reynolds numbers.
