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

Вычислительные технологии

Том 9, № 5, 2004

МЕТОД ВИХРЕЙ-В-ЯЧЕЙКАХ ДЛЯ ПЛОСКИХ СЖИМАЕМЫХ ТЕЧЕНИЙ*

Ю. Н. Григорьев Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: [email protected]

New Eulerian-Lagrangian vortex-in-cell algorithms are suggested for 2D weakly compressible medium in the framework of Oberbeck — Boussinesq approximation. General case of compressible ideal gas and flows with relaxation are considered as well.

Введение

Впервые чисто лагранжев метод дискретных (точечных) вихрей был реализован в работе [1] для моделирования неустойчивости вихревой пелены в идеальной несжимаемой жидкости. Смешанный эйлерово-лагранжев алгоритм "вихри-в-ячейках" для плоских течений идеальной несжимаемой жидкости был предложен в начале 70-х годов в [2, 3]. При этом авторы опирались на аналогию динамики вихревых прямолинейных нитей, перпендикулярных к плоскости течения и двумерных движений точечных электрических зарядов, что позволило использовать многие элементы метода "частицы-в-ячейках", успешно развивавшегося в вычислительной физике плазмы. Немалую роль здесь сыграло и то обстоятельство, что системы таких вихрей широко изучались в теоретической гидродинамике. Однако, несмотря на ряд показательных результатов расчетов тонкой структуры вихревых течений [2] и создание известного коммерческого пакета VORTEX [4], методы "вихри-в-ячейках" для задач гидродинамики в последующие годы не получили столь широкого распространения, как, например, методы "частицы-в-ячейках" в математическом моделировании задач физики плазмы. Более того, на волне стремительного прогресса компьютерной техники в вычислительной гидродинамике появились новые чисто лагранжевы алгоритмы, такие как SPH (Smoothed Particle Hydrodynamics) [5, 6], которые стимулировали совершенствование известных лагранжевых вихревых методов. Так, в работе [7] метод дискретных вихрей, восходящий к [8], был обобщен на невязкие плоские течения в слабо сжимаемой жидкости. В последнее время в работах [9, 10] (см. также библиографию в [9]) предложен и апробирован лагранжев алгоритм вихревых частиц для общего случая плоских сжимаемых течений с учетом вязкой диссипации и теплопроводности.

Можно предположить, что застой в развитии метода "вихри-в-ячейках" связан с тем, что долгое время все методы этой группы развивались исключительно в прикладных рабо-

* Работа выполнена при поддержке Российского фонда фундаментальных исследований (грант № 0301-00160), RFOSS (грант № 2314.2003.1), а также интеграционного гранта СО РАН (№ 2).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.

тах, зачастую производя впечатление искусственного know how. В то же время математическое обоснование лагранжевых алгоритмов дискретных вихрей [7, 11] сыграло важную роль в их практических приложениях. Поэтому, хотя трактовку алгоритмов "частицы-в-ячейках" с точки зрения вычислительных методов можно считать завершенной [12], сложившееся у многих вычислителей впечатление все еще сохраняется.

Вместе с тем, как известно [11], лагранжевы алгоритмы вихревых частиц далеко не экономичны. Даже в случае несжимаемой идеальной жидкости соответствующий им объем вычислительной работы оценивается как O(N2) и даже как O(N3) (N — число лагранжевых частиц). Для сжимаемых течений, несмотря на использование в [9] некоторых современных вычислительных приемов из [5, 6], снижающих трудоемкость, эта оценка сохраняется. В то же время метод "вихри-в-ячейках" [3, 4, 11] относится к категории экономичных и имеет оценку объема вычислений порядка O(N)(O(NlnN)).

В некоторых работах по лагранжевым алгоритмам утверждается [7], что смешанные методы "частицы-в-ячейках" имеют более низкий уровень точности. Однако такое утверждение не соответствует практике применения этих методов, например, в вычислительной физике плазмы, где они получили наибольшее развитие. Несложные оценки [12] показывают, что достижение в них второго порядка точности не требует от вычислителя чрезвычайных усилий.

В отличие от идеальной несжимаемой жидкости, где завихренность в плоских течениях только переносится [13], что упрощает построение алгоритмов вихревых частиц, в сжимаемых течениях имеются внутренние механизмы генерации завихренности. В частности, это конвективные движения неоднородных жидкостей в силовых полях, бароклинный механизм возникновения вихрей за фронтом искривленных ударных волн и в процессе развития неустойчивостей контактных границ. Поэтому обобщение известных лагранжевых и смешанных методов вихревых частиц [3, 11] для моделирования сжимаемых потоков является далеко не очевидным. Тем более, что в процессе развития лагранжевы методы частиц все более отличаются от аналогичных эйлерово-лагранжевых. Это связано с тем, что в лагранжевых алгоритмах многие свойства моделируемой среды — вязкую диссипацию, сжимаемость или другую реологию — необходимо отразить непосредственно в "частицах", что приводит к значительному усложнению.

В последние годы появились важные с точки зрения приложений задачи расчета завихренных течений с объемными деформациями, в которых с целью исключить влияние схемной вязкости использовались лагранжевы вихревые методы [9, 10]. Очевидно, что здесь можно было бы успешно применить соответствующий вариант метода "вихри-в-ячейках", альтернативного как конечно-разностным схемам, так и чисто лагранжевым методам частиц.

Нам известна единственная работа [14], где упоминается об использовании метода "вихри-в-ячейках" в сочетании с известными программными комплексами, базирующимися на конечно-разностных схемах, в расчетах неустойчивости Рихтмайера — Мешкова на границе раздела двух идеальных газов. Однако в статье полностью отсутствуют какие-либо детали реализации метода, а ограничения, введенные в исходную постановку задачи, фактически исключают влияние сжимаемости, а вихревая компонента поля течения предполагается малой. Но для математического моделирования наибольший интерес представляет применение метода как раз за рамками сделанных авторами предположений. В связи с этим последовательная разработка алгоритма "вихри-в-ячейках" для сжимаемых течений представляет интерес как с точки зрения вычислительной математики, так и с точки зрения возможных приложений.

1. Метод "вихри-в-ячейках" для слабо сжимаемых течений

Модель слабо сжимаемых течений используется в тех случаях, когда малые изменения плотности в потоке определяются вариацией температуры среды или/и концентрации примеси, а не изменением давления. Такие течения возникают, например, в результате свободной конвекции в поле силы тяжести. В частности, это могут быть конвективные движения в морской воде с неоднородным распределением температуры и солености. Аналогичные потоки, характеризуемые малыми числами Маха, при практически постоянном давлении возникают в неоднородно нагретой атмосфере Земли.

В общем случае подобные течения описываются системой Обербека — Буссинеска [15]. Она состоит из уравнения импульсов

Ро(ddu + uVu) = -yP + - ат(T - To) + ас(C - C°)]g + VS (1)

и двух конвективно-диффузионных уравнений для температуры T и концентрации пассивной примеси C

дТ

— + uVT = кт AT, dC

— + uVC = Dc AC. (2)

Система замыкается уравнением состояния в форме

p = p°[1 - ат(T - T°) + ас(C - C°)]. (3)

Здесь p — плотность жидкости (газа); p° = const — плотность в однородном стационарном состоянии; u = (u,v) — вектор скорости; p — давление, S — тензор вязких напряжений, V, A — операторы градиента и Лапласа на плоскости r = (x, y).

Для дальнейшего целесообразно перейти к упрощенной системе, тем не менее сохраняющей связь с исходными уравнениями Обербека — Буссинеска. Так как перенос пассивного скаляра, которым являются температура и концентрация примеси, легко включается в ла-гранжев этап методов "частицы-в-ячейках", а соответствующие диссипативные процессы рассчитываются на основе конечно-разностных схем на эйлеровой сетке [12], мы не будем рассматривать уравнения (2). Из тех же соображений опустим в уравнении импульсов (1) слагаемое VS, описывающее вязкую диссипацию. Так как динамика температуры и концентрации не рассматриваются, положим также

p°[1 - ат(T - T°) + ас(C - C°)]g = pF,

где F = (fi, f2) — некоторая заданная массовая сила, которую, не сужая рамок исходной системы, можно положить потенциальной (консервативной), так что

F = Vф, V х F = rotF = 0.

Вместо уравнения состояния (3), связывающего плотность с вариациями температуры и концентрации, для плотности p введем уравнение неразрывности в обычной форме [13] с условием несжимаемости Vu = 0.

В результате получаем систему уравнений в виде

др + иУр = 0; (4)

ди _ Ур р . .

— + иУи =--+ Е. (5)

дЬ ро ро

С точностью до постоянного множителя р0 в последнем слагаемом правой части (5) полученная система описывает развитие конвективной неустойчивости в идеальной несжимаемой неоднородной жидкости из первоначального состояния покоя [7, 15]. Эта система близка также к полной системе уравнений гидродинамики идеальной неоднородной жидкости [16].

Проведенные преобразования позволят при построении алгоритма "вихри-в-ячейках" не отвлекаться на второстепенные детали, одновременно очерчивая достаточно широкий круг задач для его приложения.

Для дальнейшего необходимо из уравнения (5) получить уравнение динамики завихренности. Так как рассматриваются плоские течения, то завихренность имеет одну компоненту

дь ди дх ду

перпендикулярную к плоскости течения.

С учетом несжимаемости жидкости и консервативности силы Е из уравнения (5) имеем

и = Ух и = го1и = — — — = ьх — иу, (6)

+ иУи = Рх/2 — ру ¡1- (7)

При этом поле скорости выражается через функцию тока ф, которая связана с завихренностью уравнением Пуассона

Дф = —и, и = фх, ь = —фу. (8)

Для системы (4), (7), (8) на временном шаге т от момента времени Ь введем следующую схему расщепления

+ У(р1и 1) = 0, + У^и 1) = 0; (9)

Р1(Ь) = р(Ь), и (Ь) = и(Ь), й 1(Ь) = и1(Ь); (10)

др2 П ди2 ~ г ~ г пг\

-дц~ = 0, "дТ = й1х/2— й1у¡1; (11)

Р2(Ь) = й1^ + т), и2(Ь) = й1(Ь + т); (12)

Дф2 = — £й2, ^ |ж = (иБ)|дп, (13)

дп

ий2 = ф2х, й = —ф2у -

Здесь п, б — единичные векторы внешней нормали и касательной к границе области течения.

Система (9), (10) описывает лагранжев этап переноса [12], моделируемый на лагранже-вой сетке вихревых частиц. Система (11)—(13) соответствует эйлерову этапу, на котором она интегрируется на эйлеровой сетке с помощью какой-либо разностной схемы.

Введенное разбиение исходной системы можно рассматривать как расщепление по физическим процессам. Система (9), (10) отражает перенос вещества и завихренности полем скорости, индуцированным распределением завихренности в поле течения. Такое поле по аналогии с внутренними электромагнитными полями в плазме [12] можно назвать самосогласованным. Соответственно, система (11), (12) описывает генерацию завихренности в неоднородной жидкости в поле массовой силы Е. Завихренность определяет самосогласованное поле скорости через уравнения (13).

Опишем кратко реализацию лагранжева этапа, используя стандартную схему метода "частицы-в-ячейках" [12]. Система (9) в векторной форме записывается как

^ + У(ия) = 0. (14)

Решение представляется в виде следующей интерполяционной формулы:

N

я(г) = £ О,Л(г, г,- (*)) (15)

3 = 1

на лагранжевой сетке частиц, которые в данном случае мы будем называть вихрями.

Здесь О, = (т,, Г,) — вектор признаков, переносимых ]-м вихрем, компоненты которого т, и Г, есть масса и циркуляция данной частицы соответственно. Функция Л (г, г, (¿)) — ядро модельной частицы, обладающая следующими универсальными свойствами

Л(г, р) = Л(р, г), ^ = - ^ (16)

с условием нормировки в области течения П

J Л(г, р)вг = 1.

п

Пусть <^(г) — произвольная гладкая финитная функция в П. Подстановка (15) в (14) и интегрирование с весом <^(г) дают

/ <?(г)( + У(ияП ¿г =

,5 О* | !гг - «и^иаг

Е Ц ^"Г J Л(г, г,)У^(г)вг - I Л(г, г,)иУ^(г)вг| = 0. (17)

7 п п

Здесь при преобразованиях были использованы свойства симметрии (16) ядра Л(г, г,) и очевидные соотношения

5Л(г, г,) 5Л(г, г,) ^г, 5Л(г, г,)

дг, ^ 5г ^

Сделав в последних интегралах замену переменных г ^ г,, получаем

N " '"¿г,-

¿0,/ Л(г, г, )У^(г,)

,

,=1 п

ж - и(г' >

¿г, = 0.

В силу произвольности функций ^(г) и Д(г, Гj) интегралы в последнем равенстве должны обращаться в нуль тождественно. Для этого необходимо выполнение условий

¿г ■

- и (г,) = 0, з = !,..., N. (18)

Эта динамическая система описывает движение модельных вихрей. Интегрируя ее на шаге т, мы в силу (15), (17) приближенно воспроизводим пространственно-временную эволюцию решения системы (9), (10).

Очевидно, что в процессе перемещения частиц переносимые ими признаки {Oj}^=1 сохраняются. После окончания лагранжева этапа необходимо проинтерполировать пере-

г р+1т N -

несенные признаки из новых местоположений частиц {г^ }-=1 в узлы эйлеровой сетки. Пусть {5а}М=1 — регулярное декартово разбиение области течения П на ячейки, центры которых образуют эйлерову сетку узлов. При интерполировании необходимо выполнить интегральные законы сохранения массы и завихренности, которые записываются в виде

N „ N М

/ п(г,ГjМг= £д, = £q«lU (19)

j=1 П j=1 а=1

Здесь — площадь декартовой ячейки номера а. В переходе к последнему равенству подразумевается, что для вычисления интеграла по П используется квадратурная формула "средних" прямоугольников для двумерного случая. При этом величины qa = {ра ,Соа} имеют смысл кусочно-постоянных плотностей массы и завихренности в расчете на единицу площади. Очевидно, что законы сохранения будут выполнены, если значения сеточных плотностей вычислять по формулам

1 г , 1 N - N

ISalJ ^ \Sa\ U

X .— 1

aa

х J — 1 s 3 —1

оа оа

qdr = —J2 Qj R(r, ri )dr = Q. R(r> r. ), a = 1,...,M. (20)

_1 J A_1

Здесь

R(r'r.) = \5~\ ' R(r' r.)dr

a

Sa

— сеточное ядро [12]. Суммирование в (20) формально распространяется на все частицы, но используемые в практических расчетах ядра обычно имеют конечный носитель, так что каждая частица дает вклад лишь в несколько ближайших ячеек, с которыми ее носитель имеет непустое пересечение. Например, часто применяемое в расчетах динамики плазмы ядро PIC интерполирует в четыре ближайших узла.

Детали процесса интерполяции "частицы-сетка" для конкретных сеточных ядер можно найти в [12]. В результате мы получаем набор сеточных функций на (p + 1)-м временном слое qa+1 = Ж1,^:1}.

Как следует из системы (11), (12), плотность на эйлеровом этапе остается постоянной. Для эйлерова этапа распространенным является использование явных конечно-разностных схем. По найденной таким образом сеточной функции завихренности {ш2а} на основе некоторой разностной аппроксимации краевой задачи (13) для уравнения Пуассона

= Wf^K = (21)

a

вычисляется сеточная функция "2+1. Для нахождения компонент скоростей {и2+1,'/,р+1} обычно используются центральные разности. В результате получаем набор сеточных функций

= К+1,<\<+1,<+1}, (22)

характеризующих состояние жидкости (газа) на (р + 1)-м временном слое.

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

м

Qj = {ш,, Г,-} = я«|№(г«, г,-), з = (23)

а=1

где интерполирующая функция выражается через сеточное ядро

5(г„, г,) = ^ Я(г, г,)аг. (24)

Непосредственно проверяется, что при этом интегральные масса и завихренность также сохраняются.

Так как сеточная функция скорости не является плотностной характеристикой, ее пересчет выполняется по формулам

м

и, = и (га, г,), з = 1,...,Ж. (25)

а=1

Следует отметить, что в прямой интерполяции по формулам (20) и в обратной интерполяции (23), (25) можно использовать разные сеточные ядра.

2. Метод "вихри-в-ячейках" в общем случае плоских сжимаемых течений

При построении алгоритмов "частицы-в-ячейках" важно выбрать необходимую форму исходной системы уравнений, которая определяет схему расщепления исходной задачи на эйлерову и лагранжеву подсистемы и набор признаков, переносимых модельными частицами.

Будем исходить из системы уравнений газовой динамики в виде [13]

др + У(ри) = 0; (26)

ди и2 1

^ + У(у) - и и = - - Ур; (27)

дг

— + У(еи) = -(7 - 1)еУи. (28)

Система (26)-(28) замыкается уравнением состояния идеального газа

p = (Y - 1)е. (29)

Кроме переменных, описанных в разделе 1, здесь е = pcvT — внутренняя энергия единицы объема газа, 7 — показатель адиабаты. Так как мы рассматриваем плоские течения, вектор завихренности ш = V х u = rotu по-прежнему имеет единственную ненулевую компоненту, перпендикулярную к плоскости течения (см. уравнение (6)). Однако здесь в отличие от случая слабосжимаемых течений Vu = divu = 0. Поэтому аналогично подходу, предложенному для чисто лагранжева метода вихревых частиц [9], от уравнения импульсов (27), записанного в форме Лэмба, удобно перейти к двум уравнениям для скалярных величин ш и в = Vu = divu. Функцию в будем называть дилатацией или расхождением поля скорости.

Применяя к уравнению (27) операцию rot с использованием известных формул векторного анализа [13], приходим к уравнению

дш 1

^ + (uVV = -шв + — Vp х Vp. (30)

It p2

Соответственно, применение к (27) операции div дает уравнение

Iе + (uV)e + (Vu) : (V)T = -^Ap + 1 VpVp, (31)

It p p2

где (Vu) : (V)T = uX + + 2uyvx — тензорная свертка. Так как вектор завихренности нормален к плоскости течения, в уравнении (30) отсутствует характерное для трехмерного случая слагаемое (wV)u, которое описывает эффект растяжения вихревых нитей.

Слагаемое — Vp х Vp является источником генерации завихренности по бароклинному p2

механизму и может быть существенным в ударных волнах с искривленным фронтом [13] или при развитии неустойчивости типа Рихтмайера — Мешкова. В баротропных течениях с уравнением состояния p = p(p) градиент давления

^ dp_, Vp = -f- Vp

dp

необходимо коллинеарен градиенту плотности и данное слагаемое зануляется. При этом уравнение (30) превращается в уравнение Гельмгольца [13] в плоском случае, где генерация завихренности происходит только в результате сжимаемости.

В сжимаемых течениях поле скорости разлагается на соленоидальную и потенциальную компоненты [13]. В нашем случае

u(x,y) = us(x,y) + ug (x,y). (32)

Соленоидальная составляющая us вычисляется через вектор-потенциал A = {0,0, единственная ненулевая компонента которого есть функция тока ^(x,y). Потенциальная составляющая ug выражается через скалярный потенциал <^(x,y). Имеем соответственно

u = V х A = {^y, -^x}, ug = V<£.

С учетом этих выражений, применив к (32) операции div и rot, для функций ip(x,y) и ф(х,у) получаем следующие краевые задачи:

д<р дф Д^ = 0, — = un|dn, Дф = —и, — = us|dQ. (33)

дп дп

Преобразуем левые части в уравнениях (30), (31) к дивергентной форме, а в слагаемых справа выразим давление p из уравнения состояния (29). В результате приходим к следующей системе уравнений:

ж + V<pu> = °

ди _, ч y — 1 / ч

— + V(uu) = (px£y — Py£x),

(34)

д0 Y — 1 Y — 1 — + V(0u) = —2(uy Vx — UxVy) + (Px£x + Py £y)--— Д£,

^ + V(£u) = —(7 — 1)£0.

К полученной системе общая схема построения метода "частицы-в-ячейках" [12] применяется очевидным образом. Вводится расщепление на временном шаге т. Система уравнений лагранжева этапа имеет вид (сравни (9),(10))

др

~8t

дй\

~дГ

двг dt

д1 + V(p iü1) = ü, р 1 (t) = p(t), + V(^iüi) = ü, £i(t) = u(t),

+ V(0iüi ) = ü, h (t) = 9(t), (35)

^ + V(£iüi) = 0, Zi(t) = e(t), üi(i) = üi(t).

Соответственно, для эйлерова этапа (сравни (11), (12))

^ = 0, P2(t)= Р l(t + Т),

д£2 Y — 1

-df = - Pix е iy — р iy £lx, U2(t) = £i(t + Т),

дв2 Y — 1

— = — 2(Uiy V ix — UixV iy) +---z^- (p ixSix + (36)

Y — 1

+P iy z iy) — -Цг— Ае i, 9 2 (t) = 9 i (t + Т); P i

—2 = —(y — 1)Zi9 i, -2(t) = Zi(t + т ). t

Уравнения эйлерова этапа дополняются системой уравнений для вихревой и потенциальной компонент поля скорости

д ~ Z дФ2 ,

А< 2 = 92, -щ = ün|öQ, Ug2 = <2xVg2 = <2y, дф

Аф 2 = —£>2, = üs|sn, Us2 = ф 2y ,Vs2 = — Ф2x■ (37)

Как и выше, введенное расщепление выделяет две группы физических процессов. Система (35) описывает перенос признаков, характеризующих локальное состояние сжимаемого газа: плотности, завихренности, дилатации и внутренней энергии суммарным полем скорости (32). Система (36) отражает процессы локальной генерации переносимых признаков. В частности, завихренность порождается бароклиническим слагаемым, которое аналогично правой части уравнения (11) для завихренности в неоднородной жидкости. Локальное изменение дилатации в и внутренней энергии е определяется сжимаемостью газа. Кроме того, можно отметить, что в уравнении для в слагаемое с производными скоростей по структуре аналогично члену (шУ)и в уравнении для вихря в пространственном случае, где он соответствует деформации вихревых трубок полем скорости.

Видно, что система (35) уравнений лагранжева этапа в векторной форме имеет стандартный вид (14). Модельные частицы-вихри вводятся интерполяционной формулой (15), в которой векторы д, переносимых частицами признаков определяются как

д, = {ш,,г,в,-}, з = 1,...,ж (38)

Здесь в число переносимых признаков, кроме массы ш, и циркуляции Г,, включены ди-латация частицы в, и внутренняя энергия Е,.

Реализация лагранжева этапа осуществляется по описанной выше схеме. Модельные частицы-вихри перемещаются в соответствии с уравнениями (18). После окончания лаг-ранжева этапа сеточные плотности переносимых признаков вычисляются по формулам (20). В результате получаем набор значений сеточных функций на (р + 1)-м временном слое

Я« = {Р+1 в?+1е1+1}, а = 1,...,М, (39)

которые являются начальными данными для эйлерова этапа, реализуемого на основе некоторой конечно-разностной аппроксимации системы (36).

Как и выше, из (36) следует, что сеточная плотность газа {р2а} на эйлеровом этапе остается постоянной. В отличие от системы (11), (12) здесь при вычислении дилатации необходимо знать сеточные функции скорости. В частности, можно использовать значения скоростей, рассчитанных в заключение эйлерова этапа на предыдущем временном шаге. Другой подход состоит в интерполяции на эйлерову сетку скоростей модельных частиц после окончания лагранжева этапа на данном временном слое. При этом, однако, процедура интерполяции должна быть согласована с интерполированной сеточной плотностью {р1+1}<м=1 так, чтобы удовлетворялась некоторая конечно-разностная аппроксимация уравнения неразрывности

Р+1 _ р

^-^ + ^({р 1«и1«) = 0,

т

где — разностный оператор дивергенции. Вариант такой консервативной процедуры изложен в [17].

Обычно на эйлеровом этапе используют явные по времени схемы. В данном случае при решении уравнений для дилатации в2 и внутренней энергии е2 можно дополнительно применить зейделевскую стратегию, подставляя в одно из уравнений значение другой величины, вычисленной на данном шаге.

По найденным сеточным функциям завихренности и дилатации на основе разностных аппроксимаций краевых задач (38)

А^^1 = в2+1, = ,

= -шР+\ ^Й1 =

находятся сеточные функции потенциала и тока "О^1}, по которым вычисляются

потенциальная и вихревая компоненты сеточной функции скорости и0,+1. При этом целесообразно использовать центральные разности

~Р+1 ~Р+1 7р+1 7р+1

"Р+1 = "Р+1 = "Р+1 + "Р+1 = - ^г-Ц + "2,^+1 - "2^-1

"а — "г.? + "ад + 2Л, '

~Р+1 ~Р+1 7р+1 7р+1

?,Р+1 — "Р+1 = „Р+1 + ?,Р+1 = ^2,г,^+1 - ^2,г,^-1 + "2,г+1,^ - "2,г-1„7 ^а "г^ + Чед 2^, + 2^ '

Переход на следующий временной слой завершается обратной интерполяцией рассчитанных значений сеточных функций (39) на частицы по формулам (23) и вычислением скоростей частиц по соотношениям (25).

3. Метод "вихри-в-ячейках" для сжимаемых течений с релаксацией

В заключение сделаем замечание об обобщении алгоритма, изложенного в предыдущем разделе, на случай сжимаемых течений, сопровождаемых некоторым релаксационным процессом. Это могут быть химические реакции, диссоциация-рекомбинация, термали-зация внутренних степеней свободы в молекулярных газах. Аналогичными уравнениями описываются течения жидкостей с пузырьками в квазигомогенном приближении и некоторых сред со сложной реологией. Во всех этих случаях релаксация приводит к относительно слабому диссипативному эффекту. Однако его точная количественная оценка может представлять значительный практический интерес. В этой связи использование конечно-разностных методов приводит к нежелательному затушевыванию эффекта из-за соизмеримого диссипативного влияния схемной вязкости. Естественной альтернативой, оптимальной с вычислительной точки зрения, здесь могли бы служить методы "частицы-в-ячейках".

Рассмотрим систему уравнений двухтемпературной газовой динамики [18], описывающей течения молекулярного газа с возбужденными колебательными модами, релаксиру-ющими к термическому равновесию. В ней уравнение неразрывности для плотности р, а также уравнения для завихренности ш и дилатации в имеют тот же вид, что и в системе (34). В уравнение энергии входит дополнительное слагаемое характерного вида

| + У(еи) = ^ - (7 - 1)е*

Собственно релаксационный процесс описывается уравнением Ландау — Теллера [18]

дЦ + ^ и) = ^ • №

Здесь тг — характерное время релаксационного процесса; е = с^ рТ — удельная энергия квазиравновесных поступательных и вращательных степеней свободы газовых молекул, характеризуемых статической температурой потока Т. Соответственно, е^ = с^ьрТ — энергия колебательных мод с текущей колебательной температурой Т; е^ — энергия колебательных мод в термическом равновесии.

Таким образом, отличие от системы (34) заключается в уравнениях для энергий е и е^. В этом случае система уравнений лагранжева этапа дополняется уравнением переноса для колебательной энергии

^е^ + У(е^1й!) = 0, ё„1(г) = е^ (г).

В систему уравнений эйлерова этапа (36) включаются уравнения энергии в следующем виде:

^ = - (7 - 1)еД, е2(г)= е~1(г + т),

дг тг

де ^2 е,и1 е VI ~ /,\ ~ /, . \

=-, е„2(г) = е„цг + т).

дг тг

В вектор переносимых модельными вихрями признаков (38) добавляется колебательная энергия Е^,

д, = {т,, Г,, в,-, Е, ,Е„,}, 3 = 1,...,Ж

Расчетная схема совпадает с описанной в предыдущем разделе. На эйлеровом этапе, как и для системы (36), можно использовать зейделевскую стратегию.

Список литературы

[1] RosENHEAD L. The formation of vortices from a surface of discontinuity // Proc. Roy. Soc. London A. 1931. Vol. 134. P. 170-192.

[2] Christiansen J.P. Numerical simulation of Hydrodynamics by the Method of Point Vortices // J. Comp. Phys. 1973. Vol. 13, N 3. P. 363-379.

[3] Baker G.R. The "Cloud in Cell" Technique Applied to the Roll Up of Vortex Sheets //J. Comp. Phys. 1979. Vol. 31, N 1. P. 76-95.

[4] Christiansen J.P. VORTEX. Two-dimensional Hydrodynamics Simulation Code. L.: Culham. Lab. Rep., CLM-106. HMSO, 1970.

[5] Gingold R.A., Monoghan J.J. Smoothed particle hydrodynamics: theory and application to nonspherical stars // Mon. Not. R. Astron. Soc. 1970. Vol. 181. P. 375-387.

[6] Monoghan J.J. Particle methods for hydrodynamics // Comp. Phys. Rep. 1985. N 3. P. 71-85.

[7] Anderson C.R. A vortex method for flows with slight density variations //J. Comp. Phys. 1985. Vol. 61. P. 417-444.

[8] Chorin A.J. Numerical study of slightly viscous flow //J. Fl. Mech. 1973. Vol. 57. P. 785-796.

[9] Eldredge J.D., Colonius T., Leonard A. A vortex particle method for two-dimensional compressible flow // J. Comp. Phys. 2002. Vol. 179. P. 371-399.

[10] Eldredge J.D., Colonius T., Leonard A. A dilating vortex particle method for compressible flow //J. Turb. 2002. N 3. P. 1-10.

[11] Leonard A. Vortex methods for flow simulation //J. Comp. Phys. 1980. Vol. 37. P. 289-335.

[12] Григорьев Ю.Н., Вшивков В.А. Численные методы "частицы-в-ячейках". Новосибирск: Наука, 2000.

[13] Кочин Н.Е., Кивель И.А., Розе Н.В. Теоретическая гидромеханика. Ч. II. М.: ФМ, 1963.

[14] Бахрах С.М., Волков С.Г., Куратов С.Е. и др. Использование алгоритма "вихри-в-ячейках" при численном моделировании гидродинамической неустойчивости // Тр. Забаба-хинских науч. чтений-2003. (Электронная публикация www.vniitf.ru/rig/konfer/)

[15] Джозеф Д. Устойчивость движений жидкости. М.: Мир, 1981.

[16] Ландау Л.Д., Лившиц Е.М. Теоретическая физика. Т. VI. Гидродинамика. М.: Наука, 1986.

[17] Grigoryey Yu.N., Vshiykoy V.A., Fedoruk M.P. Numerical "Partical-in-Cell" Methods. Theory and Applications. Utreht-Boston: VSP, 2002.

[18] Лосев С.А. Газодинамические лазеры. М.: Наука, 1977.

Поступила в редакцию 19 июля 2004 г.