CC BY
11
УДК 519.63:532.5
ТЕНЗОРНО-ИНДЕКСНОЕ ПРЕДСТАВЛЕНИЕ И РАЗНОСТНАЯ АППРОКСИМАЦИЯ ДЛЯ КВАЗИГАЗОДИНАМИЧЕСКОЙ
СИСТЕМЫ
Д. В. Тихомиров
(.кафедра математики)
Исследована квазигазодинамическая система (КГД-система) на сложных расчетных областях. Построен общий вид КГД-системы в тензорно-индексном представлении в произвольных координатах, а также разностная аппроксимация полученных уравнений в потоковой форме в виде законов сохранения.
Введение
Истоки научного направления исследования задач газодинамики на основе системы квазигазодинамических уравнений (КГД-системы), рассматриваемого в настоящей работе, принадлежат Б.Н. Четве-рушкину и Т. Г. Елизаровой [1]. С момента основания нового подхода был проведен ряд аналитических и численных исследований [2]. Так, с помощью КГД-системы проведены численные расчеты течения перед уступом [3]; течения в окрестности плоской пластины [4]; колебательных течений вблизи цилиндрического тела с выступающей иглой и вблизи полого цилиндра [5] и др.
Настоящая работа продолжает исследования КГД-системы и расширяет класс существующих задач на сложные расчетные области. В данной работе построен общий вид КГД-системы в тензорно-ин-дексном представлении в произвольных координатах, удовлетворяющих требуемым условиям гладкости; построена разностная аппроксимация в потоковой форме в соответствии с записью КГД-уравнений в виде законов сохранения [6, 7].
1. Постановка задачи
Рассмотрим задачу обтекания твердого тела в приближении квазигазодинамической системы. В общем случае задача формулируется в безразмерных переменных в виде системы дифференциальных уравнений в частных производных, начально-граничных условий и условий замыкания:
др
Ы (фи) т
дЕ_
т
Л |
сН
Ш) + Ур = ШуП, ■р)] + <1г^ = Шу(Пи).
(1)
Соответствующие газодинамические величины, входящие в систему, выражаются следующими формулами:
3 = р(
\¥= -(сИри®и) + V?), р
II II,,,+ -ио [р(и\7)и + Ур] + + г/[(иУ)р + 7рШуи], = р [(V ® и) + (V ® и)т — 2/3/сНуи
р 7
т =
Рг 7 — три М
Ые•Эс р
1
VI?
р
7 ~~
!/-(2 —(_/„,
-(и\7) ( -
v р
/и2
1 »-V
7-1 р
Здесь р — плотность среды; и — газодинамическая скорость; Е — полная энергия; р — давление среды; J — вектор плотности потока массы; — КГД-добавка вектора плотности потока массы; П — тензор вязких напряжений КГД-системы; —
тензор вязких напряжений системы Навье-Стокса; Ся — скорость звука в среде; q —вектор теплового потока КГД-системы; т-параметр КГД-системы; -у — показатель адиабаты; р — коэффициент динамической вязкости; Ые — число Рейнольдса; Рг — число Прандтля; Эс — число Шмидта; М — число Маха; I — единичный тензор-инвариант 2-го ранга.
Начальные, а также граничные условия на область входящего потока, твердого тела и свободной границы задачи обтекания для КГД-системы сформулируем следующим образом соответственно:
р к=о= 1,
и |4=0= -и0,
[р |*=о= 1/7,
^ = 0 <Уп
< и = 6>,
др
' др
дп
= 0,
дп
дУх
дп
др
= 0,
= 0,
^ = 0, ^ = 0,
<Эп
<Эп
дп
= 0.
Здесь вектор п — внешняя нормаль к соответствующей границе.
2. Общие положения тензорного анализа
Рассмотрим некоторый базис ej в евклидовом пространстве [8], любой вектор в котором задается координатами (ж1,ж2,ж3). Определим переход к новому базису Гг с помощью следующего преобразования (здесь и далее предполагается суммирование по повторяющимся индексам):
гг = Це0, Щ =
дх3 дх1
:г = хг (
Ж-*- /у» /у»
> Ои > Ои
!)
где Щ — матрица перехода от базиса ej к новому базису Гг, любой вектор в котором задается координатами (ж1,ж2,ж3). При этом метрический тензор, отражающий дифференциальные свойства геометрии области имеет вид
дхк дхк дх1 дхз
(2)
Рассмотрим некоторый вектор и = игег в базисе е^, где и% — соответствующие координаты вектора в этом базисе. При переходе к новому нормированному базису (ёг = у/\дц\) координаты данного вектора преобразуются по следующей формуле:
С/-7 = иЩ \/
431'
иг =
и*
>/¡9.
-Л,
(3)
зз\
где Ъ = Ь-1 — матрица обратного перехода от базиса Гг к базису ej.
Представим основные формулы тензорного анализа, необходимые для построения КГД-еиетемы в тензорно-индексном представлении.
Формула вычисления скалярного произведения векторов:
(АВ )=дыАкВ1. (4)
Формула вычисления ковариантной и контравари-антной производной скалярной функции:
ял
™ т-- =
(5)
Формула связи координат вектора в нормированном базисе ëj (11г) и в базисе г^ (11г):
иг =
хр
\f\9i-
(6)
Формула вычисления ковариантной производной тензора первого и второго ранга (с последующей сверткой):
■ д111 3 дхз
итт
т,р
а дгз
дх%
Гк = 1 и (д9Ц , д9й ддц * 'Г \ дх* дх3 дх1
(7)
(8) (9)
где Г^- — символы Кристоффеля. Формула вычисления дивергенции вектора —
^ = (9 = йегШ\). (10)
3. Приведение КГД-системы
к тензорно-индексному представлению
Последовательно преобразуем все уравнения системы (1) к тензорно-индексному представлению, используя формулы раздела 2. Для начала рассмотрим первое уравнение системы вместе с величинами участвующими в данном уравнении. Преобразование уравнения и соответствующих величин проводится по формулам (5), (6), (8)—(10):
др Ж
д
\Мдх
I )
РГ = -
Т3с
д
дхз II1
и3
= 0,
и1
У= ( ^
уУ
и* 4 ит
ттР
из
а 9Р
9x3
Тензорно-индексное представление второго уравнения системы (1), с помощью формул (4)-(8), вместе с величинами и П принимает вид:
и1
д_
г-7 I Л1 П% у/Ы
9x3 \ у/\д1(
Гг I Р
тз{
ит
а др дхз
м
п" = п£
из
дх з
йк
Г3^
д
гг грт ГП] '
и1
др
ит>
=г:
+тд
,31
ик
Л
\П9г.
др 1
+7Р
тк
Лк
дхк
д
Ък
■¿к
дхк
д
у/ЫдхГ'
ш
ит
\f\9r,
и1
лк
д
дхк \V\m\ и*
ит
>/\9 тт
=г;
тк
г-7
тк
2/3 д
:ЗК
1__д_
у/Ы-
ит
\f\9n
Используя формулы (4)—(10), выразим третье уравнение КГД-еиетемы и q в тензорно-индекеном представлении:
дЕ_
т
_1__д_
Е■
Р
_1__д_
\f\g\dxi
Г (У№)
р
д
Е =
9ы ик
и1
уН/
2 Vы
р ,
т Р.
9> =
Р
7 ^
д
из
Рг 7 — 1 дхк
1 ик д
-тр
\/У
и?-1 у\д,
'кк
дхк
Р
Р! Х/Й
+р
ик д
дхк
Р)\
С учетом формулы (3) начальные условия в тен-зорно-индексном представлении имеют вид
' р |*=о= 1, [р |*=о= 1/7-
Граничные условия на область входящего потока, твердого тела и свободной границы с помощью формул (3)-(6) преобразуются к следующему виду соответственно:
Р= 1, _ <Р = 1/7,
7=0,
др
дх Ц* = 0
др дх
7=0,
7=0,
др
дх' _д_
др
и1
= 0,
. дх
т=0.
Здесь индекс г соответствует индексу координатной оси вдоль нормали к области.
4. Разностная аппроксимация КГД-системы в тензорно-индекеном представлении
Покроем трехмерную расчетную область равномерной сеткой с шагами Ь,\, /13 • Символом ш обозначим множество узлов, являющихся серединами ячеек этой сетки:
ш = {(ж1,ж2,ж3): х] = ¿1 + (г — 0.5)ДХ; х] = 12 + - 0.5)Л2; х3к = Ь3 + (к^ 0.5)Л3;
* = 1,
j = l,...,N2^, к = 1,
Здесь ¿1,1*2, ¿з-координаты начала границы расчетной области.
Все газодинамические величины будем относить к узлам из и. Значения произвольной функции ф из р^и1,112,173,р в средних точках вычислим по формулам:
^п±1=0.5(^„±1 + ^„), (П)
Для остальных ф = ф(р, II1, II2, II3,р):
(12)
Фгзк = Ф(Ргзк, (и^ук, (и2)ук, (и3)ф,Рф). (13)
Приведем результат построения разностной аппроксимации дифференциальных уравнений КГД-системы в тензорно-индекеном представлении, используя формулы предыдущего пункта и (11)—(13).
Уравнение баланса массы:
д*
|9
1 I тта
\9,
т
2>
К
тш __/
<У(%г±|) ~ Р(Ьа± |) ^
177:
(%г±±)
к
гого(1т±|) I
(%г±|)
Т(А
ИТГ
(%г±|)
РЦшЦ)
ттп
х 1 л == х
1»пп(гт±|)(г„ + |)1 7 -Р(«г)(«„-§)
ттгп
1$тт(гт±±)(г„ + §)1
2 '
19пп(гт ± 4 ) (¿п — 4 ) I
5'тт(гт±|)(гп-|)
¿7.
•Г"
ЫатЦ)Р(гтЦ)
(Ьг,М)
ТГТП
Цш+1)
№
р т
кп(гт±\)Р{ь-п±Ь2)
тт(цп±|)1 у 15тт(гт±|) I
'П
(¿т +1)
|5и(гт±А)1 Л/ Ьпп(гт±|)
Уравнение баланса импульса:
ттпт _ттпт
Р\ V
+ т?пи1п + ттжк
^Аят,
тп ттт, тп
Цп + 1) Цп+1) Цп-1)
1 Л™
1
'тт(1п +1) I
'тт(г„-|)1/ "
ги1ит - т^пг-ик
9т,т,
\f\9mm\g
т,к
П"1
■та
_дпт
гк
9кк
ТТП
пп(г„±|) I
ик.
РЦп±\) и™ , г
ит
'тт(г„±|)(гк + |) I 9т,т,{г„±^){гк-^) I ^ {ьМ)
И1,. . х
2 рт
(¿«±1)1
+ 5
• Т/,
2 ' ш
'(¿»±§)5(«„±±) Х
ик.
Ьк
'(¿»±5)
А
0П(г„±±)(гг + ±)
\
0П(г„±±)(гг-±)
п
та __I п1
ттт (¿п±|)Й + |)
,тт(г„±|)(гг + |)1
ттт
'тт(г„ ± |) (г; — |)
ик.
ы
±¿1
рт
ТТП
(¿п±|)Й + |)
'пп(г,г±|)(гг + |)1
кк^п^-^) I
ТТП,
'««(^„¿^(гг-!)^ 1
ик
2 ! рП
-2/3
„та %г±|)
\
\
Уравнение баланса энергии:
_
(г„±|)(^-|) I ^
/<; /<; + л/1
\
рг пк
\
УиЦп-Ь
т пк
(Ьг 2) г) / к.
у/Ш
\9(г„
ЕапЦ) +РапЦ)
4)1
Еа -1
(1-п 2
Л)
■2' \ ТП
Тп I -1.
1
ы
7
%„±1) =--т,- ,9« ^
Ра,г±Ык-Ь)\ 1
т
1
' Т(г»±|)Р(г„±|)
пп(г,г±|)
7-1
кк(г„М) I +
ик
Р(гп±\Ш-Ь)\ 1 , "(¿»±5)
2 2 *
г»±|)1
1
1
^»±1) ^(¿»±1)
1
Для обеспечения устойчивости разностной схемы при больших числах Рейнольдса, введем специальный регуляризатор:
= М 1 2 Ь_ Т Re.Scр 8 С/
Здесь а — параметр регуляризации (числовой коэффициент); Л — размер ячейки при разностной аппроксимации (Н = \/{Н\)2 + (Лг)2 + (^з)2)-
Алгоритм нахождения плотности, компонент скорости и давления на каждом временном слое состоит из двух этапов. Сначала заполняются фиктивные ячейки на основании граничных условий, далее вычисляются значения искомых величин во внутренних точках области. Критерием окончания итерационного процесса вычислений является выполнение следующего условия:
Рг Рг
тах
А t
< е,
где е — параметр точности.
Замечание 1. Метрический тензор и символы Кристоффеля характеризуют геометрию области и не зависят от времени, что позволяет при программной реализации проводить вычисление этих величин только один раз и на всех шагах по времени использовать уже вычисленные значения.
Замечание 2. При программной реализации возможно как аналитическое задание перехода к координатам (с помощью задания явных функций), так и численно, используя таблицу соответствия координат. В случае численного задания значения метрического тензора и символов Кристоффеля вычисляются по разностным аналогам формул (2) и (9).
Замечание 3. Для устранения особенностей в точках на ребрах и углах расчетной области необходимо разбить область на подобласти, окружающие тело, и проводить расчет последовательно в каждой из полученных областей, как предложено в [6].
Замечание 4. При расчетах в координатах, обладающих свойством периодичности (например, угловая координата ф в цилиндрических координатах), необходимо проводить операцию сшивания значений в фиктивных точках на границе смены периода.
Автор выражает глубокую признательность профессору Т. Г. Елизаровой за постоянную поддержку и помощь в работе.
Литература
1. Елизарова Т.Г., Четверушкин Б.Н. // ДАН СССР. 1984. 279, №1. С. 80.
2. Шеретов Ю.В. Математическое моделирование течений жидкости и газа на основе квазигазодинамических и квазигидродинамических уравнений. Тверь, 2000.
3. Граур И.А., Елизарова Т.Г., Четверушкин Б.Н. // Дифф. уравнения. 1986. 22, №7. С. 1173.
4. Elizarova Т.О., Graur I.A., Lengrand J.C., Chpoun А. // AIAA J. 1995. 33, N 12. Р. 2316.
5. Антонов А.Н., Елизарова Т.Г., Четверушкин Б.Н., Шеретов Ю.В. // Журн. вычисл. матем. и матем. физ. 1990. 30, №4. С. 548.
6. Шеретов Ю.В. Применение функционального анализа в теории приближений. Тверь, 2001.
7. Елизарова Т.Г., Соколова М.Е. // Вестн. Моск. ун-та. Физ. Аастрон. 2003. №5. С. 19.
8. Дубровин Б.А., Новиков С.П., Фоменко А.Т. Современная геометрия. Методы и приложения. Т. 1, 2. М., 1998.
Поступила в редакцию 21.01.05
