Аппроксимация уравнений квазигазодинамики на треугольных сетках Текст научной статьи по специальности «Математика»

УДК 517.958:533.7

АППРОКСИМАЦИЯ УРАВНЕНИЙ КВАЗИГАЗОДИНАМИКИ НА ТРЕУГОЛЬНЫХ СЕТКАХ

Т.Г. Елизарова*), В.В. Серегин

(.кафедра математики) E-mail: seregin@afrodita.phys.msu.ru

Рассматривается метод решения квазигазодинамических уравнений, описывающих течение вязкого газа в двумерном случае, в областях сложной формы с использованием неструктурированных пространственных сеток. Для построения разностной схемы используется интегро-ин-терполяционный метод (метод конечного объема). Аппроксимация системы квазигазодинамических уравнений производится явной по времени разностной схемой на треугольной сетке, удовлетворяющей принципу триангуляции Делоне.

Система квазигазодинамических (КГД) уравнений в безразмерных переменных имеет вид [1, 2]

др

dt дри

дЕ

at '

СНул = 0, (1)

+ ® и) + \7р = <Иу П, (2)

+ = (з)

Здесь р — плотность среды; и — газодинамическая скорость; р — давление; Е — полная энергия; j — вектор плотности потока массы; П — тензор вязких напряжений КГД-еиетемы; q — вектор теплового потока КГД-еиетемы.

В отличие от уравнений Навье-Стокса в КГД системе присутствуют дополнительные дивергентные слагаемые с параметром, имеющим размерность времени. Дополнительная диссипация, присутствующая в КГД-уравнениях, обеспечивает эффективность численных алгоритмов, построенных на их основе.

Перепишем систему уравнений (1)-(3) в векторном виде

Ш

— сИу = 0, (4)

где векторы ииШ введены следующим образом:

U =

w, -j;

t р\

PUx

PUy

\EJ

,W =

( Wi\

w2

W3 \w 4,/

= П

' u^pe,

Wi = Пи — q

E + p.

-J-

P

Построим треугольную сетку таким образом, чтобы ее узлы попадали на границу области расчета,

число узлов в области расчета было достаточным для описания течения с большими градиентами, а также сетка удовлетворяла бы принципу триангуляции Делоне [3], в соответствии с которым требуется, чтобы в круг, описанный около любого треугольника, не попадало ни одной узловой точки, отличной от вершин указанного треугольника. Фрагмент такой треугольной сетки приведен на рис. 1.

Рис. 1. Фрагмент треугольной сетки

Сетка представляет собой систему треугольников с вершинами Мо, Мх,..., М{, Мп (рис. 2). Газодинамические величины и,р,р,Е заданы в узлах сетки Л/;. Для аппроксимации исходных уравнений потребуются значения газодинамических величин в центрах треугольников. В качестве центров примем точку пересечения медиан. Эти точки обозначим Ра, Ръ ■ ■ ■! Рк: Рк (см. рис. 2). Газодинамические величины в центре треугольников определяются как среднее арифметическое величин в вершинах треугольников.

о

Институт математического моделирования РАН.

Рис. 2. Контрольный объем

Для построения разностной схемы воспользуемся интегро-интерполяционным методом (методом конечного объема) [4]. Для каждого узла Л/; треугольной сетки строится контур, состоящий из точек Р'к пересечения медиан треугольников, содержащий данный узел. Обозначим число узлов контура Ь через К. Область, ограниченная этим контуром, представляет собой расчетную ячейку — контрольный объем. На рис. 2 приведен контрольный объем с центром в точке Мо и границей Ь. В соответствии с методом конечного объема проинтегрируем систему (4) по контрольному объему. Для интеграла от воспользуемся формулой Оетроградекого-Га-усса, связывающая интеграл по объему с интегралом по поверхности или формулой Грина в двумерном случае:

д_ т

и ¿в — £ (• п) ¿1 = 0.

(5)

Здесь Б — площадь контрольной ячейки, Ь — контур ячейки, п — внешняя нормаль к контуру.

Далее к первому слагаемому в (5) применим формулу среднего значения, получим

(6)

1

где и = — ф иёБ — среднее значение скорости и

в центре расчетной ячейки Б. Далее обозначим и через и.

Аппроксимируем контурный интеграл следующим образом:

Н к (7)

+^у(Рк+1/2)пу(Рк+1/2)) ьк.

Здесь значения потоков Шх, Шу и проекций нормалей пх,Пу вычисляются в точках Рк+1/>2 отрезков Ьк = Ь(Рк, Рк+г), составляющих контур Ь.

Итак, получим разностную формулу для метода конечного объема:

О О ^ У (№х(Рк+1/2)пх(Рк+1/2)+

Б к (8)

+Щ(Рк+1/2)Пу(Рк+1/2))

При расчете по формуле (8) потребуются значения частных производных по ж и у от плотности, скорости, давления в точках Рк+1/>2 отрезков Ьк. Значения газодинамических величин в серединах отрезков определяются как средние между значениями на концах. Рассмотрим способы для выражения производных. Обратимся к рис. 3 и рассмотрим, например, четырехугольник Ма,Рх,М2,Р2. Рассмотрим величину /, определенную в вершинах четырехугольника и найдем ее частные производные Ш- в точке Р-.

ду

3/2 •

У

2(Р2)

1(М2)

0(Р1)

~ 1<Ъ "

3(Мо)

х

Рис. 3. Контур для вычисления частных производных

Частные производные можно искать двумя способами: используя формулы Грина (9)—(10), т. е через интеграл по контуру 0123,

дв

д£ ду

£ 5

/% Л,

(9) (10)

и используя производные по направлениям 1\,12 (11)—(12):

21

дк 91 т2

дf дх дх д1\ д/ дх дх д12

<9/ ду ду дЬ' <9/ ду ду д12'

(Н) (12)

Раскрывая интегралы в формулах (9)—(10) по аналогии с (7) и решая систему (11)-(12) отноеи-можно показать, что эти способы

тельно

й£ дх'

й£

ду

эквивалентны. В результате получим формулу, выражающую производную в Р3/»2 через координаты узлов контура 0123 и значения функции в этих

узлах:

д/ (/2 /о)(уз 2/1) "Ь (/1 /з)(|/2 2/о)

<Эж

2А

= (/2 - /о)(ж! - ж3) + (Д - /з)(ж0 - ж2) ду 2А

(13)

(14)

где А — площадь четырехугольника 0123.

Далее, выразив компоненты вектора через газодинамические величины и частные производные, можем решить систему (8) относительно р,и,р. Соответствующая разностная схема приведена в приложении.

В качестве примера расчета рассмотрим задачу о распаде сильного разрыва в двумерном случае. Слева от разрыва р = 8, р = 480, справа р = 1, р = 1. Значения других величин: и = 0, 7 = 5/3, Рг = Эс = 1, а = 0.5 (см. приложение). Шаг по времени Д£ = 0.002, графики приведены для момента времени £ = 4. Для сравнения с одномерным решением этой задачи будем рассматривать сечение вдоль оси х. На рис. 4, 5 изображены графики плотности и давления.

80 100 120 140 ^160

Рис. 4. Распределение плотности

80 100 120 140 ж 160

Рис. 5. Распределение давления

Проведенные расчеты показывают, что при сгущении сетки решение двумерной задачи стремится к автомодельному решению.

Приложение

Распишем разностную схему (8):

Д*.

Л - Рг ~ X] {зк+1/2пк+1/2 + ^'1+1/21/2)

^ - сГХ! (|Щч-1/2 ^ Зк+1/2ик+1/2 ~Р\ Х

Рг

ХПк+1/2 + [П*;+1/2 ~ И+1/2ик+1/2\ Пк+1/2) Ьк

/2 Зк+1/2ик+1/2\ Пк+1/2~

[П1+1/2 Зк+1/2ик+1/2 р] П1+1/г)

ТТХХ X ^т\ху иу

11к+1/2ик+1/2 ^ 11к+1/2ик+1/2

х Ек+1/2+Рк+1/2 ,х

Чк+1/2 3 к+1/2

Рк+1/2

Щ+1/2ик+1/2 + Щ+1/2ик+1/2 ~ Чк+1/2

Ек+1/2 +Рк+1/2 .у

ьк+1/2"

Рк+1/2

'Зк+1/2

Пк+1/2

Рг = (7 - 1) ¿г - Л

(й?)2 + (йГ)2

Здесь индекс г — номер узла сетки, к — номер узла контура Ь,

т-^К т / ч -1/2

,Ък=0 Ьк,к+1 Рк \

Тк =а-

К

Рк

_ Рк+ Рк+1 _ Рк+1/2 ~ -^-' Рк+1/2 ~ 2

Рк=Рк &СТк, Рк +Рк+1

Рк+1/2 ~

Щ+1/2 ~

Рк + Рк+1

Пк + Щ+1

Е,

к+1/2 ~ Рк+1/2

тк+1/2 ~ (ик+1/2У + (ик+1/2

Тк + Тк+1

2 > 2

2 I !,.У \2

1

7-1

Рк+1/2,

и!

к+1/2

_ тк+1/2 Рк+1/2 д

го

к+1/2

П+1/2

рихиу

к+1/2 к+1/2 \9х/ к+1/2

д_

дх

рихиу

Рк+1/2 / й;+1/2

% ) к+1/2 \9У/ к+1/2.

Зтк+1/2 = Рк+1/2(и- = Рк+1/2 К -№?)

П:

пэ к+1/2 — Рк+1/2

4 _ 2 /диу\

з V дх ) к+1/2 3 V ду ) к+1/2

Tixy _

linsk+l/2 — lJ'k+1/2

~^Vnsk+1/2 = ßk+1/2

= Pfe+1/2

fdux\ _ (диУ\

К дх j k+1ß \ дх / k+1/2_

YdvT\ _ /диГ\

\dx)k+1/2 \dy)k+1/2

3 V % ) k+1/2 3 \a®/fc+i/2.

n^^fc+l^ = Щ^+1/2 +Tfe+l/2 jwilt-l/2 Pfc+1/2X

X ("^/2 fc+1/2 + «Ï+1/2 (^J fc+i/2) +

^ (DJ-- d),+ l/2 +

/диП \i + 7Pfc+l/2 ^ ^ gx J ^^ + ^ dy J ^J j ,

X

n^fc+1/2 = Щ»к+1/2 + Tk+1/2 \uXk+1/2 pk+1/2 X

\%Л+1/2

Pfe+1/2X

n^fc+i^ = n»ïfc+1/2+rfc+1/2|«»+1/2

\dxj k+1ß

k+1/2 = 4ïk+1/2 + Tk+1/2 |«|+1/2 Pfe+1/2 X ( x (диУ\ У (duV\ \

x Vk+1/2 4+1/2 K^J k+1/2 J

+iPk+i/2 \\ dx J k+1/2 +\ dy J k+1/2/

%+l/2

Pfe+1/2 7 Pr 7 — 1

\dx p)t

+u

Tk+l/2ßk+l/2Uk+1ß у i д p\

L7

fc+1/2 \dypJk+1/2

k+1/2 j_( x (9_P

<¡1+1/2

Pfe+1/2 7 Pr 7 —

^+1/2Р/г+1/2%+1/2

9 p 1 V%P

1 О '(dp

k+1/2

7-1

fc+l/2

dx p ; k+1/2 8 l\

+u|+1/2 llpJfc+i/2J+Pfc+1/4<+1/2 I^pJ

fc+l/2

+lt

fc+1/2

(--) )

Р/ k+1/2/

k+1/2

Здесь К — число узлов в контуре L, Рг — число Прандтля, Sc — число Шмидта, 7 — показатель адиабаты, а — регуля-ризатор [0,1].

Авторы благодарны И. В. Попову за сотрудничество.

Литература

1. Шеретов Ю.В. Математическое моделирование течений жидкости и газа на основе квазигидродинамических и квазигазодинамических уравнений. Тверь, 2000.

2. Елизарова Т.Г., Шеретов Ю.В. // Журн. вычисл. матем. и матем. физ. 2001. 41, №2. С. 239.

3. Попов И.В., Поляков C.B. // Матем. моделирование. 2002. 14, №6. С. 25.

4. Самарский A.A. Теория разностных схем. М., 1989.

5. Кудрявцев Л.Д. Курс математического анализа. Т. 2, М., 1981.

6. Рождественский Б.Л., Яненко H.H. Системы квазилинейных уравнений. М., 1978.

Поступила в редакцию 16.06.04