CC BY
152
Вычислительные технологии
Том 7, № 6, 2002
КИНЕТИЧЕСКИЕ СХЕМЫ И ВЫСОКОПРОИЗВОДИТЕЛЬНЫЕ МНОГОПРОЦЕССОРНЫЕ ВЫЧИСЛЕНИЯ В ГАЗОВОЙ ДИНАМИКЕ
Б. Н. Четверушкин Институт математического моделирования РАН
Москва, Россия e-mail: buzz@imamod.ru
Some problems connected to the use of multiprocessor computer systems for numerical simulation of gas dynamic problems are discussed on the basis of original approach, i.e. kinetically coherent finite difference schemes. The problems of algorithms adaptation to the multiprocessor systems with distributed memory architecture are considered. The difficulties connecting to the developing of modern parallel software, are discussed. The results of numerical simulation of some actual gas dynamic problems using high performance computer systems are presented.
Введение
В настоящей работе обсуждаются две на первый взгляд слабо связанные научные проблемы. С одной стороны, это развитие нетрадиционного подхода в вычислительной газовой динамике, основанного на использовании тесной связи между кинетическим и макроскопическим описаниями сплошной среды — кинетически согласованных разностных схем (КСРС) [1]1. Такой подход в известном смысле приводит к переформулировке базовых уравнений механики сплошной среды. Другая проблема, с которой пришлось столкнуться в последние годы, вызвана трудностями при использовании многопроцессорных суперЭВМ. Эта проблема в первую очередь связана со сложностью адаптации алгоритмов решения задач, описываемых уравнениями математической физики, к архитектуре многопроцессорных вычислительных систем с распределенной памятью. Такие системы наиболее распространены [2] в настоящее время. Имеются также серьезные трудности при создании математического обеспечения для таких многопроцессорных систем.
Многие из проблем, связанных с использованием многопроцессорных вычислительных систем для решения задач газовой динамики, могут быть решены с помощью применения кинетически согласованных схем, опирающихся на представления статистической механики и механики сплошных сред. Демонстрация возможностей такого подхода для высокопроизводительных вычислений представляет одну из целей данной работы. Более того, эффективные многопроцессорные алгоритмы, построенные на основе КСРС, позволяющие
© Б.Н. Четверушкин, 2002.
ХВ зарубежной литературе такие схемы часто называют больцмановскими.
проводить массовые расчеты с использованием систем, которые имеют быстродействие 102 Гигафлоп (1011 операций с плавающей запятой в секунду) и выше, дают возможность для развития, казалось бы, далеких непосредственно от проблем механики жидкости и газа параллельного математического обеспечения (программ многопроцессорной визуализации, динамической перезагрузки процессоров при расчете и т.д.).
1. Кинетически согласованные разностные схемы
1. Как известно, уравнения газовой динамики — Эйлера или Навье — Стокса — могут быть получены с помощью процедуры Чепмена — Энского из более сложного кинетического уравнения Больцмана [3, 4] для одночастичной функции распределения f:
f + Zif = // (f(x,C',t)(x, CV) - f(x,С,t)f (x, с!,*)) |gIdChdbde.
Здесь (i — скорость молекулы в направлении оси Oxi; Z', — скорости, которые принимают молекулы со скоростями £ и после столкновения; |g| = — (^1|; b — прицельное расстояние, на котором происходит столкновение молекул; е — угол в плоскости перпендикулярной к вектору g, 0 < е < 2п.
Зная функцию распределения f, можно вычислить все остальные газодинамические параметры [1, 3, 4]:
Р = J fd(, pu = j f Сd(, Рг] = j fCiCjd(,
3 г с2 с C2
-kT = j f-2d(, q = J fC^d(, C = С — u,
где p — плотность; u — вектор скорости; T — температура; Pj — тензор напряжений; q — вектор теплового потока.
Однако обычно связь между кинетическим и макроскопическим описаниями не используется при построении вычислительных алгоритмов газовой динамики. Алгоритмы для решения систем уравнений Эйлера и Навье — Стокса конструируются без учета их органической связи с уравнением Больцмана.
В кинетических схемах в качестве основы принимается дискретная модель одночас-тичной функции распределения, осредняя которую по скоростям молекул, можно сразу получить разностные уравнения для газодинамических параметров. Здесь уже на алгоритмическом уровне происходит согласование кинетического и газодинамического описаний сплошной среды.
Продемонстрируем процедуру получения КСРС на примере одномерного по пространству случая. Простейшая разностная дискретизация уравнения Больцмана, схема с направленными разностями, примет вид
fi — fim I fi+1m — fi-1m 1 Ax fi+1m — 2fim + fi-1m / f / (л i\
-л--+ Z-:-=--:—~-, —oo < Z < oo. (1.1)
At s 2Ax 2 Ax2 ' " s " v ;
При этом сделано естественное для газодинамического описания предположение, что на j-м слое по времени функция распределения постоянна на отрезке x £ [xi,xi+1] и принимает локально-максвелловский вид2
fi = i = JPL3/2 e^. (1.2)
(2nRTi)3/2
2Заметим, что на плодотворность использования кусочно-постоянного представления функции распределения обращалось внимание еще в монографии [5].
Обратим внимание на отсутствие интеграла столкновений в правой части (1.1), который для максвелловской функции распределения обращается в нуль.
Последовательно умножая (1.1) на сумматорные инварианты 1, £, С2/2 и интегрируя по скоростям молекул, получим разностную схему для газодинамических параметров [6]:
puerf (s) + 1 exp (—s2)
(pu)j+1 — (pu)j ( 2 ) . Ax 2 ) .. . pu exp (—s2) 4 '% 4 '% + (pu2 + p) x =- (pu2 + p) erf (s) + ' K J
At v ; 2
Ej+1 — E%% . ^ m, . Ax
"^"Ai + (u (E + P)) x=-
u (E + p) erf (в) + —П (e + 2) exp (—s2)
pu2 1
ние; E — полная энергия; e — внутренняя энергия; E =--+ e; в = , ; s = ^u; erf (s)
в^ 1 2 ' ^ v " 7 • (L3) Здесь использовались дополнительные обозначения: p = pRT — газокинетическое давле
pu2 1
--+ e; в = . ;
2 ^ \plRT
x
2 С
— функция ошибок erf (x) = exp (—z ) dz; ax — центральная разность для аппрок-
у/п J 0
симации первой производной, axx — разностная аппроксимация второй производной.
Обратим внимание на члены в правых частях системы (1.3), которые являются членами искусственной диссипации. В отличие от классической искусственной вязкости фон Неймана, применяемой в вычислительной газовой динамике, вид этих членов согласован с диссипативными членами в разностном аналоге кинетического уравнения (1.1). Это послужило поводом дать схемам название кинетически согласованных.
2. Дадим краткий обзор работ по кинетическим схемам. В [7] на основе дискретной модели для функции распределения построена разностная схема для задач газовой динамики в лагранжевой системе координат3. В [8] строилась схема, использующая принципы, достаточно близкие к тем, что применялись в работе [6]. Однако диссипативные члены в этом алгоритме в явном виде выражены не были. В работах [9, 10] газодинамические параметры непосредственно определялись интегрированием конечно-разностного решения уравнения Больцмана.
Получаемые на основе упомянутых работ данные, с одной стороны, показывали их применимость для нахождения газодинамических параметров, что, учитывая связь уравнения Больцмана и уравнений газовой динамики, неудивительно. С другой стороны, качество расчетов с учетом вычислительных затрат уступало существующим на начало 80-х годов традиционным методам решения уравнений Эйлера. В то же время полученный в [6] (1983 г.) алгоритм сразу стал конкурировать с традиционными подходами и немедленно вызвал цепь дальнейших публикаций [11-18]. Следует отметить, что совпадающие с КСРС [6] схемы были получены в [19-20] (1986 г.) и также послужили основой исследований течений невязкого газа, подробно описанных в обзоре [21]. Появившиеся в конце 80-х годов Lattice Boltzmann Schemes (LBS)4 в настоящее время широко используются для расчета целого ряда газодинамических течений [22]. В этих схемах вместо уравнений газовой
3Эта работа выполнена в качестве дипломной под руководством автора настоящей статьи. К глубокому сожалению, вскоре после ее завершения В. Поткин трагически погиб.
4Решетчатые больцмановские схемы.
динамики решается кинетическое уравнение для ограниченного набора лучей, направления которых жестко согласованы с пространственной сеткой. По своей сути LBS весьма близки к КСРС.
Начиная с 1985 г. интересы автора статьи были сосредоточены на использовании КСРС для моделирования течений вязкого газа, включая моделирование неустановившихся течений, турбулентных потоков, течений химически реагирующих газов, решения задач аэроакустики и аэроупругости. С 1988 г. большинство расчетов проводились с использованием многопроцессорных вычислительных систем. По указанному направлению имеются более ста публикаций, включая монографию [1].
Естественно, что моделирование течений вязкого газа с использованием КСРС невозможно без понимания связи уравнений, базирующихся на их основе, с традиционными физическими моделями, применяемыми для описания таких течений. Частично эти вопросы будут рассмотрены ниже. Среди работ, имеющих определенную связь с обсуждаемым подходом к моделированию течений вязкого газа, укажем работы [23-26].
В рамках данной статьи автор не ставит своей задачей обсуждение достоинств и недостатков указанных работ. Отметим лишь, что эти работы отличает желание по-новому взглянуть на уравнения, описывающие сплошную среду. В этом контексте особо хочется отметить работы [25, 26], в которых так же, как и КСРС, вычислительный аспект и формулировка физической модели органически связаны.
3. Сформулируем в самом общем виде способ получения КСРС через баланс функции распределения. Предположим, что в момент времени t = tj в некоторой ячейке I (рис. 1), имеющей произвольную форму, величина f может быть представлена в виде максвеллов-ской функции (1.2) f/m, постоянной внутри ячейки. В течение времени At = tj+1 — tj имеет место бесстолкновительный разлет молекул. При t = tj+1 происходит мгновенная максвел-лизация молекул, причем выполняется основное положение процедуры Чемпена — Энского
[3, 4]:
У fj+V (Z) dZ = J fm+1 V (Z) d(, (1.4)
где v (Z) — сумматорные инварианты; fj+1 — значение функции распределения перед максвеллизацией. В итоге балансное соотношение для функции распределения примет вид
f J+1 _ fj i k г fj I fj i k r
i=1 c i=1 с
Рис. 1.
Здесь mes (VI) — объем (площадь для двумерной задачи) центральной ячейки; иц — нормаль к K-й грани (ребру) ячейки.
Чтобы получить КСРС на произвольной сетке, достаточно проинтегрировать (1.4) с сумматорными инвариантами ^ (Z) = 1, 2/2.
Замечание 1. В случае одномерной задачи, интегрируя выражение (1.5), придем к системе разностных уравнений (1.3).
Замечание 2. Полученные таким образом КСРС соответствуют первому члену разложения Чепмена — Энского разностного аналога уравнения Больцмана [1, 27].
Замечание 3. Рассмотренная процедура построения КСРС позволяет легко строить схемы до определения газодинамических параметров на любой последовательности сеток, включая приобретающие большое значение при моделировании на высокопроизводительных системах, неструктурированных сетках.
Балансное уравнение (1.3) можно упростить, введя кинетическое время т, за которое частицы в среднем пересекают расчетную ячейку. Для членов, описывающих аппроксимацию по пространству, используя технику, развитую в [28], можно получить дифференциальное приближение. При этом выведем дифференциально-разностное уравнение, описывающее изменение функции распределения:
Р+1 _ fj Ji Ji
fim
At
+ div (C, /
j \ -m,
2 ltl/
(1.6)
где £/ = С, (/).
Выражение (1.6) является основой для получения квазигазодинамической системы уравнений (КГУ) [16], используемой в большинстве расчетов. Чтобы получить КГУ в пространственно-двумерном случае, проинтегрируем (1.6) по скоростям молекул с весами
р (С) = 1, с,<2/2:
Pj+1 - Pj , д M
+ дХ +
дтд , 2 Л дтд , 2 Л дтд.
1 дХ 2 дх^ + р + дУ 2 ду^ + р + дХду(puv),
( pu)j+1 — ( pu)j д д
+ ^(pu2 + р) + ду (puv) =
At
дтд , о ч дтд , 2
дХ 2 дх^ +3pU + ду 2 ду^
(pv)j+1 — (pv)j д
+ pu) + дхду (pu2v + pv)
д
At
+ дх (puv) + ^ (pv2+p
дУ
д т д ( 2 дх 2 дх ^
дтд д тд
+ H + ду 2 ду (pv° +3pv) + дХду (pv2u + ,
E^1 - Ej д д
___ + _u (E + p) + „ (E + p) =
u
дтд дх 2 дх
д трдр
5
E +2 Р
12 + ^ pv
+
5
E +2 Р
+ ^ pu
Y
Y — 1
+ 1 д т д /р
дх 2 рдх Pr дх 2р-х \p
дтд
ду 2 ду
дтрдр 1 дт д/р ду 2 p-у Pr ду 2Рду \p
+
+
дд +—т—uv (E + р) дх ду
(1.7)
(1.8)
2
v
(7 — показатель адиабаты; Рг — число Прандтля). Значение т выбиралось в виде т = (2^)/р, что соответствует времени между столкновениями молекул, полученному в рамках элементарной кинетической теории [3, 4]: т = 1/е.
Единственным физическим допущением, используемым при выводе КГУ, было предположение о слабом изменении максвелловской функции на расстоянии длины свободного пробега и мгновенной максвеллизации при £ = £?+1.
4. Что же дает возможность использовать КГУ, внешне отличающейся диссипатив-ными членами от уравнений Навье — Стокса, для моделирования течений вязкого газа? Оказывается, диссипативные члены, стоящие в правых частях КГУ, можно разделить на две составляющие — члены естественной вязкости и естественной теплопроводности и остальные дополнительные члены. Как показывает анализ, оставшиеся члены всегда малы (имеют порядок малости О(Кп), безразмерное число Кнудсена определяется как Кп = 1/Ь, где Ь — характерный размер задачи) по сравнению с конвективными членами. Естественно, они имеют тот же порядок малости, что и члены естественной диссипации
5
там, где последние сопоставимы с конвективными, например в пограничном слое5.
Тем самым решение, полученное с помощью КГУ, не отличается от решения, полученного с помощью уравнений Навье — Стокса, где последние справедливы. Это подтверждается многочисленными сопоставлениями данных расчетов. Возникает главный вопрос: в чем же тогда смысл использования КСРС и КГУ для моделирования течений вязкого газа?
Обратим внимание на одно противоречие, возникающее при выводе уравнений газовой динамики. Как известно, традиционный вывод этих уравнений заключается в рассмотрении баланса массы, импульса и энергии для конечного объема У0. Далее, заменяя поверхностные интегралы на объемные и полагая объем У произвольным, придходим к известным уравнениям газовой динамики.
Зададимся вопросом, не существует ли каких-либо ограничений на объем У? С точки зрения кинетики газодинамическое приближение основано на близости функции распределения / к ее равновесному максвелловскому значению /т. В этом случае функция распределения / должна слабо меняться на расстоянии длины свободного пробега /. Поэтому нет смысла рассматривать объемы, линейный размер которых по порядку величины меньше /.
Во вновь полученных уравнениях, включая уравнение неразрывности, должна учитываться конечность объема У0. С другой стороны, эти уравнения не должны давать сколько-нибудь заметных отличий по сравнению с результатами, полученными на основе использования уравнений Навье — Стокса, которые прошли полуторавековую проверку индустриальной практикой.
Кинетически согласованные разностные схемы и КГУ удовлетворяют двум этим условиям. Во-первых, в силу существования максвелловского распределения в ячейке ее размеры не должны превышать длину свободного пробега (для КГУ т должно быть порядка времени между столкновениями молекул). Во-вторых, практически нет различий между решениями, полученными на их основе, и решениями, полученными с помощью уравнений Навье — Стокса.
5Последнее замечание представляется крайне важным, ведь отличие решения уравнения Эйлера от уравнений Навье — Стокса объясняется не наличием членов порядка О(Кп). В задачах газовой динамики число Кнудсена мало (Кп < 0.001). Это отличие определяется тем, что формально малые члены становятся сопоставимыми с конвективными (пограничный слой). Дополнительные диссипативные члены КГУ эффекта пограничного слоя не вызывают.
Полезность применения КСРС и КГУ заключается в следующем. Эти системы корректны по своему физическому построению, так как существует минимальный пространственный (временной) размер, на котором могут происходить сильные изменения параметров. Напомним, что для уравнений Навье — Стокса существуют хорошо известные трудности обоснования их корректности [29, 30]. Кроме того, иная математическая форма записи по сравнению с уравнениями Навье — Стокса дает дополнительную свободу в построении вычислительных алгоритмов.
По всей видимости, наличие внутренней корректности и обусловило успех применения кинетических схем для моделирования наиболее сложных течений вязкого газа: неустановившихся течений, задач аэроупругости и аэроакустики [31-38]. Важно отметить, что с вычислительной точки зрения дополнительные диссипативные члены можно трактовать как эффективные регуляризаторы, позволяющие развиваться естественным неустойчивос-тям и сглаживающие неустойчивости чисто счетного характера6.
Особенно важным для параллельных вычислений является то, что при моделировании многих типов вязких течений эти регуляризаторы позволяют вести счет по явным схемам с условием устойчивости типа Куранта:
т < К. (1.9)
Этот результат представляется особенно ценным, так как обычно явные схемы для данных задач требуют гораздо более жесткого ограничения
т < К2. (1.10)
5. В основе построения таких регуляризаторов лежат скорее физические, чем математические представления. Это дает возможность построения их аналогов не только при газодинамическом описании среды.
Рассмотрим, например, разностную схему (1.1) для уравнения Больцмана. Ее дифференциальным приближением в общем случае при отказе от обязательного наличия мак-велловского распределения будет уравнение
д/+0/ = // (// - /1®|+ах/т, * = 1,2' з. (1.11)
В уравнении (1.11) в отличие от уравнения Больцмана появились дополнительные диффузионные члены, пропорциональные — шагу пространственной сетки.
С увеличением степени детальности описания (уменьшения Ах) вклад этих дополнительных членов будет уменьшаться. Могут ли на сверхподробных сетках величины Ах и соответственно дополнительных членов быть сколь угодно малыми?
По построению роль этих членов заключается в сглаживании решения на ячейке Х £ [ж;,ж;+1]. Функция распределения — вероятностная величина, и для ее определения необходимо, по крайней мере, несколько десятков молекул. Эти молекулы не локализованы в одной точке пространства, а занимают объем, линейный размер которого равен /1. Для большинства газов, включая воздух нормальной плотности, имеет место соотношение /0 ^ 11 ^ /, где /0 — размер молекулы.
6Так как главной функцией регуляризаторов является сглаживание решения на расстоянии 1, при условии выполнения этой функции имеется определенный произвол в их выборе.
Таким образом, для описания функции распределения можно использовать уравнение
[39]
Ж + СЖ41 - /л) + £'1|с| Ж• (1Л2)
Так как /1 существенно меньше длины свободного пробега /, основные следствия из уравнения Больцмана, включая величины транспортных коэффициентов, останутся практически неизменными. Дополнительные диффузионные члены в (1.12) выступают, как и в КГУ, в качестве регуляризаторов, дающих большую свободу в применении численных методов и снимающих вопросы корректности для уравнения Больцмана.
Аналогичные регуляризаторы можно построить и для задач фильтрации [40], которые обычно формулируются следующим образом:
др
— + ^у (ри) = 0, и = -к grad (р) , р = ро + (р - ро) вж, (1.13)
где коэффициент к зависит как от свойств жидкости, так и от конкретной породы; коэффициент вж, как правило, велик для жидкостей, поскольку небольшое изменение их плотности приводит к сильному изменению давления.
Наблюдаемое движение фильтрующейся жидкости есть результат осреднения на масштабе порядка сотен зерен породы. Это подтверждается тем, что экспериментальный закон Дарси (1857 г.) является следствием осредненного на этом масштабе газодинамического уравнения переноса импульса [41]. Однако в уравнении неразрывности (1.13) факт осреднения непосредственно не отражен.
По аналогии с КГУ изменим уравнение (1.13), записав его в виде
I + -<-) = ¿Ш (114)
где I' — размер, сравнимый с суммарным диаметром сотен зерен породы; с — скорость звука в фильтрующейся жидкости.
Как показывает опыт расчетов [40], дополнительный диссипативный член в (1.14) не меняет результатов расчетов по сравнению с исходными уравнениями (1.13), однако его наличие позволяет эффективно сглаживать счетные осцилляции.
Система КГУ (1.7), (1.8) применялась для моделирования трансзвуковых и умеренно сверхзвуковых течений. В последнее время, используя кинетическое уравнение Энского [42], удалось построить аналог КГУ, применяемый для моделирования течений газа с небольшими скоростями, а также течения вязкой жидкости [43]. Этот алгоритм успешно поддается параллелизации [44].
До последнего времени дискуссионным оставался вопрос о применимости КГУ для моделирования течений с большими гиперзвуковыми скоростями. Связано это с сильной нелинейной зависимостью от скорости дополнительных диссипативных членов, которые при больших скоростях становятся велики и уже не выполняют свою роль регуляризато-ров. В рамках исходных КСРС эту трудность обходили эффективным увеличением длины свободного пробега в направлении большой макроскопической скорости [45]:
'эф = V + («2/с2)-
В работе [46] удалось обобщить этот подход и в рамках КГУ.
Итак, следует отметить, что построение регуляризаторов для уравнений механики сплошной среды, которые опираются на использование физически обоснованных мини-
мальных масштабов, оказывается чрезвычайно полезным при проведении высокопроизводительных параллельных вычислений.
2. Адаптация КСРС к архитектуре многопроцессорных систем
1. Рассмотрим общие проблемы, возникающие при решении задач, описываемых уравнениями математической физики, на параллельных вычислительных системах. Одним из важнейших условий эффективности такого подхода является равномерная загрузка процессоров. В противном случае наиболее загруженный процессор, медленнее всего выполняющий свою работу, будет определять скорость расчета в целом.
Вторым требованием для систем, построенных на базе архитектуры распределенной памяти, является сведение к минимуму числа обменов между процессорами. Оно связано с относительно низкой пропускной способностью каналов связи между процессорами. Если объем информации велик, то процессор может не успеть вовремя получить необходимые для продолжения расчета данные от других процессоров. Это требование наиболее трудновыполнимо. Именно оно, как правило, вызывает основные проблемы при адаптации алгоритмов к параллельным вычислительным системам.
Третье требование — экономичность по ресурсам — обусловлено значительным объемом информации, перерабатываемым многопроцессорной системой, который с трудом поддается дальнейшей обработке на обычном персональном компьютере.
Четвертым важным требованием является логическая простота алгоритма. Важность этого требования определяется быстрым развитием и сменой поколений многопроцессорной вычислительной техники, сопровождаемыми неизбежным в той или иной степени изменением математического обеспечения, которое в свою очередь заставляет постоянно переделывать действующие программы.
Обратим внимание еще на одну проблему, характерную для газодинамических расчетов на высокопроизводительных и не обязательно многопроцессорных системах. Очень часто приходится сталкиваться с такой ситуацией, когда результаты расчетов на грубых сетках более или менее соответствуют данным натурного эксперимента. Однако на подробных сетках, использование которых допускают вычислительные системы большой производительности, вместо ожидаемого улучшения согласия с экспериментальными данными наблюдается обратный эффект.
2. Рассмотрим теперь, как алгоритмы, построенные на базе КСРС и КГУ, адаптируются к архитектуре многопроцессорных вычислительных систем с распределенной памятью [1].
Одними из наиболее сложных для моделирования являются течения, характеризующиеся взаимодействием вязкой и невязкой частей потока. Как правило, расчет в каждой из частей ведется по своим собственным алгоритмам. Однако с помощью КГУ расчет в обеих частях можно вести по единой схеме. Для этого диссипативные члены, за исключением членов, описывающих естественную сдвиговую вязкость и теплопроводность (1.7), (1.8), масштабируются пропорционально шагу пространственной сетки. При этом в невязкой части потока дополнительные члены играют роль искусственной вязкости фон Неймана, а в вязкой части эти члены оказываются малыми по сравнению с членами естественной диссипации и играют роль регуляризаторов'.
7Данная схема получила название КСРС с коррекцией и начиная с 1985 г. использовалась для расчета различных задач динамики вязкого газа [1, 47, 48]
Однородные схемы оказываются наиболее удобными для равномерной загрузки процессоров. Достаточно разбить всю расчетную область на число подобластей, равных числу процессоров, и в каждой подобласти взять одинаковое число счетных узлов. Ввиду однородности схемы, когда расчет в каждом узле ведется по единому алгоритму, автоматически обеспечивается равномерная загрузка процессоров.
Схема КСРС с коррекцией обладает логической простотой. Начиная с 1988 г. эта схема применялась для проведения расчетов на восьми различных многопроцессорных системах. При этом соответствующие программы легко адаптировались и использовались на новых типах вычислительных систем.
Наиболее простыми с точки зрения минимизации обменов между процессорами являются явные схемы. В них значения величины Uj+1 на новом слое t = tj+1 по времени целиком определяются через ее значения на предыдущем слое t = tj с помощью простой вычислительной процедуры:
Nk
Uj+1 = £ агк U_1- (2.1)
k=0
В неявных схемах, когда функция Uj+1 определяется с помощью значений на старом и новом слоях по времени, для решения соответствующих систем в общем случае приходится применять достаточно сложные прямые и итерационные методы [49].
Для явных схем происходят обмены только значениями функций, лежащих в приграничных для каждой подобласти узлах, что делает их идеальными для адаптации к многопроцессорным системам с распределенной памятью.
Однако серьезным препятствием является, как правило, жесткое ограничение сверху на допустимый шаг по времени, вытекающее из требований устойчивости расчета. Например, для задач динамики вязкого газа обычным условием является условие типа (1.10). Этот факт делает практически бессмысленным расчет на высокопроизводительных многопроцессорных системах, резко увеличивая расчетное время при использовании большого числа пространственных точек.
С другой стороны для нахождения решения альтернативных неявных схем
Nk
U'+1 = £ агкU-1 + £ bipU+ (2.2)
к=0
требуются достаточно сложные алгоритмы уже при использовании однопроцессорных систем. В многопроцессорном варианте эти трудности значительно возрастают. В самом деле, для итерационного алгоритма решения системы (2.2) требуются не только быстрота сходимости, но и минимизация обменов между процессорами. Особые сложности возникают при использовании нерегулярных сеток, которые широко применяются для моделирования технологических процессов. В этом случае матрица, описывающая изменение функции Uj+1, меняет свою ленточную структуру от точки к точке, что создает дополнительные сложности при построении итерационных алгоритмов.
Для КСРС с коррекцией условием устойчивости является гораздо более мягкое, чем (1.10), соотношение (1.9). Здесь линейное уменьшение шага по времени At вместе с h представляется естественным, так как обычно большая детальность описания по пространству должна сопровождаться более подробным описанием по времени. Простота адаптации КСРС к архитектуре многопроцессорных систем дала возможность провести первые двумерные по пространству расчеты неустановившихся течений вязкого газа на решетке транспьютеров T-414 еще в 1988 г.
Как уже отмечалось в п. 1, решение, полученное с помощью КСРС и КГУ, является сглаженным на расстоянии порядка длины свободного пробега /. По всей видимости, эта корректная с физической точки зрения постановка и является причиной того, что при использовании подробных сеток происходит лишь ожидаемое уточнение картины течения. Примеры сопоставления решений, полученных на различных сетках и подтверждающие это утверждение, приведены ниже.
Здесь не обсуждается вопрос о несоответствии большего объема информации, обрабатываемого многопроцессорной системой, и возможностей персональных компьютеров. Решение этого вопроса не связано непосредственно с КСРС. Однако проведение с их помощью большого числа разнообразных многопроцессорных расчетов поставило исследователей лицом к лицу с рядом проблем параллельного математического обеспечения и заставило находить пути их решения.
3. Моделирование течений вязкого газа
Рассмотрим несколько примеров, иллюстрирующих возможности КСРС для моделирования на высокопроизводительных параллельных системах задач, представляющих особый интерес в современной гидро- и газовой динамике. Рассмотрение начнем с неустановившихся течений вязкого газа, характеризующихся тем, что даже при постоянных внешних условиях наблюдаемая картина не выходит на стационарный режим, реализуя незатухающие осцилляции газодинамических параметров. Помимо чисто практической стороны, связанной с получением конкретных амплитуд и частот осцилляций давления, моделирование этих течений имеет общетеоретический смысл, связанный с реализацией при ряде стационарных значений параметров колебательной ветви решения.
1. Моделирование неустановившихся течений вязкого сжимаемого газа было начато с помощью КСРС в 1986 г. еще на однопроцессорных ЭВМ. Рассматривались двумерные задачи, описывающие колебательные режимы около цилиндрических тел с выступающими иглами [31]. Появление многопроцессорных вычислительных систем дало возможность моделирования трехмерных задач с числом пространственных узлов 106 ... 108.
Ниже будут рассмотрены результаты моделирования неустановившихся течений вблизи прямоугольной трехмерной каверны [33, 50, 51]. Данная задача при всей внешней простоте является хорошим объектом для моделирования возникновения колебательных режимов в различных ситуациях, например, в отсутствие элемента теплозащитного покрытия космического корабля многоразового пользования, при раскрытии различных емкостей во время полета и т. д. Для расчета использовались как вычислительные мощности ИММ РАН на базе 44-процессорной системы Рагэу1ее СС (более 10 Гигафлоп), так и мощности межведомственного суперкомпьютерного центра на базе отечественной системы МВС-1000М с общей производительностью более 0.5 Терафлоп.
Характерная схема обтекаемой каверны изображена на рис. 2. Набегающий (слева) поток составляет угол ^ с линией ВС, который менялся в диапазоне 0 < ^ < 6°. Профиль потока соответствовал ламинарному пограничному слою
и 3 /г\ 1 /г
и" = 2 Ы - 2 и
т. е. менялся от нуля на поверхности до значения невозмущенного потока на границе пограничного слоя толщиной
Течение представляется в виде периодически возникающих и исчезающих крупных вихрей как в плоскости, параллельной передней грани каверны, так и в плоскости, к ней перпендикулярной. Также наблюдаются жгуты вдоль bf — кромок каверны. Заметим, что нестационарные режимы возникают лишь при определенных соотношениях геометрических размеров каверны (l, h, w) и скоростях набегающего потока. Опишем некоторые результаты расчетов. На рис. 3 изображены мгновенные распределения давления для ^ = 0. При числе Маха (отношение скорости невозмущенного потока к скорости звука) M^ = 1.35 существуют незатухающие колебания (рис. 3, а). С увеличением числа Маха амплитуда незатухающих колебаний уменьшается. При M^ > 3 колебания затухают (рис. 3, б). Эта картина хорошо соответствует данным натурных экспериментов [52, 53].
Для ^ = 0 течение в центральной плоскости каверны близко к двумерному. Для двумерного случая был предпринят вычислительный эксперимент, направленный на исследование численных решений, полученных при сгущении пространственных сеток. Интерес к исследованиям подобного рода определяется тем, что при сгущении сеток очень часто данные расчетов перестают соответствовать экспериментальным данным [54]8.
На рис. 4 изображен фрагмент поля скоростей вблизи отрывной зоны, полученный на подробной сетке из 700 х 800 расчетных точек. Видно образование мелкой вихревой структуры, которая впоследствии создает вихревой подслой отрывного слоя. Важно отметить, что крупная вихревая структура сохраняет те же параметры, что и при моделировании на грубой сетке, в том числе амплитуды и частоты давлений. Единственное, что происходит со спектром, так это его обогащение частотами колебаний меньшей интенсивности, генерируемых мелкомасштабной вихревой структурой. Та же тенденция повторяемости результатов при смещении сеток имела место в недавних трехмерных расчетах на системе МВС-1000М.
На рис. 5 приведен спектр колебаний давления Spb в точке f при различных углах
Видно, что при небольших углах (^ < 4°) течение сохраняет те же черты, что и при ^ = 0, это еще раз косвенно подтверждает корректность использованного алгоритма. При дальнейшем увеличении угла вид течения радикально меняется.
В связи с многочисленными расчетами неустановившихся режимов встал вопрос о визуализации данных многопроцессорных вычислений. При этом возникает целый ряд проблем, таких как выделение небольших фрагментов поля течения, где находятся интере-
8Одной из возможных причин подобной ситуации является следующая. На грубых сетках существует заметная схемная диффузия, эффективно сглаживающая решение, которая значительно уменьшается при сгущении сеток, что приводит к развитию счетных неустойчивостей, в свою очередь причиной которых является исходная некорректность алгоритма.
Рис. 3.
сующие особенности9, визуализация частиц газа, переходящих из подобласти одного процессора в подобласть другого.
Однако наиболее сложная проблема заключается в следующем. Огромный объем информации, обработанный многопроцессорной системой, хранится на сервере, из которого с помощью обычного персонального компьютера (РС) выбираются данные для построения конкретного изображения. Несоответствие между возможностями РС, низкой пропуск-
9Выделение мелких фрагментов в общем поле возможно при использовании индикаторов, указывающих на появление тех или иных особенностей. Например, функция rot U будет индикатором для выделения вихревой структуры. Приведенный пример типичен для современного состояния математического обеспечения многопроцессорных систем, которое может успешно развиваться лишь в тесном контакте с моделированием конкретных задач.
Рис. 4.
Рис. 5.
Рис. 6.
ной способностью обычных каналов связи и количеством информации, из которой надо выбирать исходные данные, приводит к крайне медленному построению изображения.
Выход лежит в применении технологии "клиент — сервер". В этом случае выборку информации из хранилища осуществляет в параллельном режиме высокопроизводительная система, а отобранный небольшой объем информации, необходимый для построения изображения, передается по относительно медленным каналам связи на РС. С помощью программ, построенных на этом принципе, была проведена визуализация большинства задач моделирования на основе КСРС [55].
Замечание 4. Описанная технология визуализации может быть использована для работы с любыми базами данных сверхбольшого объема.
Вернемся вновь к моделированию задач динамики вязкого газа на многопроцессорных системах. На рис. 6 изображен результат расчета слоя смешения, включая его турбу-лизацию. В принципе, с учетом внутренней корректности КСРС возможен прямой подход к моделированию турбулентности, когда все неустойчивости получаются в процессе расчета.
Однако существующие в настоящее время мощности многопроцессорных систем позволяют моделировать лишь неустойчивости крупных масштабов или в лучшем случае — (см. рис. 4) средних масштабов. Неустойчивости более мелких масштабов, характерные
для развитой турбулентности, остаются на существующих сетках неучтенными.
Для моделирования такого рода течений использовался подход, недавно предложенный А. Бегу1еих (ШША, Франция) и И. В. Абалакиным (ИММ РАН). Этот подход является развитием идеи, высказанной в [56]. Суть метода заключается в следующем. На подробной сетке, допустимой с точки зрения производительности системы, проводится расчет поля течения. При этом напрямую определяются крупномасштабные неустойчивости. Оставшийся вклад турбулентности более мелких масштабов определяется с помощью (к — е)-модели турбулентности [57]10.
Интересным представляется использование КСРС для моделирования явлений аэроакустики, весьма важных и в то же время трудных задач современной вычислительной газовой динамики.
Перенос звуковых возмущений моделируется с помощью линеаризованных КГУ (1.7), (1.8). Линеаризованные уравнения, как и исходная КГУ, позволяют вести расчеты по явной схеме с условием устойчивости (1.9), что дает возможность легко адаптировать алгоритм к архитектуре многопроцессорных систем [1, 36, 38]. Тем самым проблема моделирования звуковых возмущений в рамках КСРС решается достаточно просто. На рис. 7, а, б приведены спектры акустического давления БР1, фиксируемого в точке С обтекающего каверну потока (рис. 7, в), при задании случайной по времени частоты искусственного источника. Расчеты проводились на различных сетках. Данные, представленные на рис. 7, б, получены на в два раза более подробной сетке, чем на рис. 7, а. Интересно обратить внимание на совпадение частот и амплитуд основных колебаний.
2. Приведенные выше расчеты, представляющие достаточно большой физический и газодинамический интерес, опирались на простые геометрические конструкции. Однако в реальных приложениях чаще используют конфигурации сложной формы. Около таких тел газодинамические поля хорошо аппроксимируются с помощью неструктурированных сеток, треугольных для двумерного случая и тетраэдных для трехмерного. Аппроксимация с помощью неструктурированных сеток, как правило, уступает по точности варианту с использованием сеток регулярных. Однако высокопроизводительные системы дают возможность применения большого числа расчетных точек, что компенсирует недостаток точности неструктурированных сеток.
На рис. 8 изображен фрагмент двукрылового профиля. Формирование треугольной сетки проводилось на основе метода продвинутого фронта [58]. Вблизи поверхности крыла и щели сетка заметно сгущается.
Как уже отмечалось, на основе выражения для баланса функции распределения (1.4) легко построить разностную аппроксимацию для газодинамических полей на любой сетке. Построенные таким образом явные схемы обладают достаточно мягким условием устойчивости, что делает их удобными для использования на многопроцессорных системах с распределенной памятью [1, 59].
Однако при этом возникает новая проблема, ранее не встречавшаяся при многопроцессорном моделировании на регулярных сетках. Нерегулярные сетки в отличие от регулярных имеют только одноиндексную нумерацию. Поэтому тривиальное разбиение с равным числом узлов в каждой подобласти, когда, например, узлы с нумерацией от 1 до N обрабатывает первый процессор, узлы от N + 1 до 2N — второй процессор и т.д., приводит к образованию подобластей с большой протяженностью внутренних границ. Причем мо-
10Чем меньше масштаб турбулентности, тем более обоснованно использование тех или иных моделей. Связано это с тем, что для мелкомасштабной турбулентности хорошо выполняются приближения однородности и общие положения классической теории.
в
Рис. 7.
гут образовываться несвязные подобласти. Это приводит к увеличению обменов между процессорами и уменьшению эффективности параллелизации.
Рациональному (с минимальным числом связей) разбиению на подобласти при использовании неструктурированных сеток уделяется серьезное внимание. Примененный автором с коллегами подход [55] опирался на метод спектрального разбиения, в основе которого лежит использование собственных векторов матрицы Лапласа для заданного графа [60]. На рис. 9 приведен пример такого разбиения для течения вблизи двукрылового профиля при использовании разного числа процессоров (а — 8 процессоров, б — 32 процессора).
3. Одной из важных задач, требующих для своего решения систем сверхвысокой производительности, является моделирование течений химически реагирующего многокомпонентного газа. Подобной постановкой, в частности, описывается горение природного газа метана СН4. Важной проблемой при этом является уменьшение выброса токсических компонентов.
Главным отличием (1.7), (1.8) от исходной системы КГУ в этом случае является наличие уравнений неразрывности для каждого из химических компонентов:
^ +1 ) + ду <* v ) =
0 Т 0 0 Т 0 0 0
= ах2 ах(р1,2 + р> + дУ2 ду^ + + ахта/*,/у +- • * =1 ■ (31)
Рис. 8.
Здесь г — номер компонента; ( — скорости образования и потерь ¿-го вещества в процессе химических реакций, которые в свою очередь нелинейным образом зависят от температуры Т и концентраций р^
Другим отличием от исходной КГУ является наличие в уравнении энергии (1.8) членов,
Рис. 9.
а
описывающих нагрев или охлаждение за счет химических реакций.
Замечание. По своей структуре система КГУ хорошо приспособлена для моделирования течений химически реагирующих газов, так как в правой части (3.1) уже присутствуют члены, описывающие самодиффузию различных компонентов.
Система уравнений, описывающих течение химически реагирующего вязкого газа, решается методом суммарной аппроксимации (расщеплением по физическим процессам). На первом этапе решается система уравнений газовой динамики без членов Ьг — ргQг, описывающих химические реакции, на втором — система
^ = Ь — р^г. (3.2)
В качестве примера рассмотрим задачу о горении в воздухе метанового факела, вырывающегося из плоской щели. Учитывалось взаимодействие в процессе горения следующих веществ: СН4, СН3, СН30, СН20, СНО, СО, С02, Н2, 02, Н, ОН, Н02, Н20, N2, N N0, ^0, участвующих в 65 прямых и обратных реакциях [1, 61, 62]. На рис. 10 приведены концентрации метана СН4, угарного газа С0 и одного из оксидов азота — N0. Видно (рис. 10, а),
а б в
Рис. 10.
как струя метана, имеющая начальную скорость 1.2 • 102 м/с, тормозится в воздухе. Вследствие подтока воздуха из приземных слоев происходят уменьшение концентрации метана и образование вихревых структур. Из рис. 10, б, в видно, что образование токсических компонентов происходит в узкой зоне на периферии факела. В этой зоне повышенной температуры происходит подавляющее большинство других химических реакций.
Моделирование процесса горения осуществлялось с помощью высокопроизводительных многопроцессорных систем Рагяу1ее СС и МВС-1000М. Обратим внимание на одно важное для многопроцессорной реализации обстоятельство.
Система уравнений, описывающих горение, требует для своего решения гораздо большего компьютерного времени по сравнению с системой, описывающей однокомпонентную газовую динамику, хотя бы из-за большого количества уравнений неразрывности (3.1). Однако основной объем вычислений (для данной задачи более 90 %) приходится на второй этап, описываемый внешне простой системой (3.2). Связано это с тем, что ряд из 65 учитываемых в модели химических реакций являются быстротекущими с характерными
временами на много порядков меньше газодинамического шага используемого на первом этапе. Зависимость этих реакций от температуры имеет пороговый характер: если температура меньше критического значения Т0, то реакции можно практически исключить из рассмотрения.
На рис. 11 приведено распределение вычислительной нагрузки между процессорами при условии равенства пространственных узлов в подобласти каждого из них. Видно, что на процессор 3 приходится большая часть нагрузки. Именно он будет реально снижать эффективность работы многопроцессорной системы в целом.
В этом случае необходимо применять специально разработанные подпрограммы динамической перезагрузки процессоров, суть которых заключается в следующем [2, 55].
На некотором этапе расчета проводятся замеры объема вычислений, совершаемых процессорами, и определяются наиболее загруженные из них. Затем часть расчетных точек из них (процессор 4 на рис. 12) передаются соседним (3 и 5). Из процессоров 3 и 5 с тем, чтобы они в свою очередь не снижали эффективность вычислений, расчетные точки, но уже в меньшем количестве передаются более удаленным процессорам. Данная процедура выполняется до тех пор, пока загрузка процессоров почти не сравняется. Повторение цикла перезагрузки происходит через определенное время, достаточное для заметного изменения картины течения. Следует обратить внимание на то, что алгоритмы динамической перезагрузки требуют для своего непосредственного применения достаточно подробной информации о сути протекающего физического процесса.
4. Рассмотренные в данном параграфе примеры опирались на явную схему КСРС с коррекцией, применимую для моделирования течений с большими скоростями. В рамках дальнейших исследований КСРС удалось построить вычислительный алгоритм, по своей логической структуре совпадающий с явной схемой, удобный для использования на многопроцессорных системах [43]. Этот алгоритм показал свою эффективность при моделировании течений газа с М > 0.1.
Рис. 11.
При моделировании течений по схеме [43] с меньшими числами Маха использовался прием, основанный на практической несжимаемости газов при М < 0.2, для которых характерным параметром является число Рейнольдса И,е = иЬ/д/р. Например, для расчета
Рис. 12.
конкретного течения с М = 0.005 и физическим числом И,е* используется параллельная версия алгоритма с М = 0.1 (газ фактически, как и при М = 0.005, несжимаем) и И,е = И,е*.
На рис. 13 приведены результаты расчета температурного поля в задаче моделирования синтеза газа с помощью платинового катализатора. Использование такого катализатора позволяет заметно уменьшить до 1200 ... 1300° температуры, при которых происходят химические реакции. Важно, что при этих температурах практически не образуются токсичные оксиды азота. Обратим внимание на несимметрию поля температур, возникающую из-за влияния силы тяжести на течение с небольшой скоростью.
Рис. 13.
Заключение
В данной работе автор обратил внимание лишь на небольшую часть возможных применений многопроцессорных вычислительных систем. В последние годы в ИММ РАН на этих системах решались и другие задачи, связанные с моделированием добычи нефти, переноса загрязнений подземными водами, процессов, происходящих при лазерном термоядерном синтезе, моделированием полупроводниковых приборов. Характерным для успеха многопроцессорного моделирования этих процессов явилось использование нетривиальных подходов, опирающихся на информацию практически из всех областей знаний — от построения математических моделей, допускающих параллельную реализацию, до специальных вопросов системного программирования, тесно связанных с решением конкретных задач.
Реальные возможности такого моделирования сегодня определяются не только наличием высокопроизводительной техники, а в первую очередь оригинальными решениями, полученными на базе фундаментальных исследований. Именно фундаментальный базис
высокопроизводительных вычислений ставит на повестку дня не менее актуальную проблему подготовки в ведущих вузах России специалистов подобного профиля.
Автор искренне благодарен своим коллегам С. В. Полякову, М. В. Якобовскому, Т. Г. Елизаровой, И. В. Абалакину, А. В. Александрову, Л. В. Дородницыну, А. В. Жоховой, Т. К. Козубской, М. А. Корнилиной, Т. А. Кудряшовой, Г. Н. Никишиной, М. А. Трапезниковой, Н. Ю. Романюхе, Н. Г. Чурбановой, Е. В. Шильникову, М. А. Шумкову, Д. А. Мачину, без активного творческого участия которых было бы невозможно получение результатов, лежащих в основе этой работы.
Список литературы
[1] ЧЕТВЕРУШКИН Б. Н. Кинетически согласованные схемы в газовой динамике. М.: МГУ, 1999. 232 с.
[2] ЧЕТВЕРУШКИН Б. Н. Проблемы эффективного использования многопроцессорных вычислительных систем // Информационные технологии и вычисл. системы. 2000. №2. С. 22-34.
[3] КОГАН М. Н. Динамика разреженного газа. М.: Наука, 1967. 440 с.
[4] ЧЕПМЕН С., КАУЛИНГ Т. Математическая теория неоднородных газов. М.: Иностр. лит-ра, 1960. 550 с.
[5] ВЛАСОВ А. А. Статистические функции распределения. М.: Наука, 1966. 440 с.
[6] ВолчинСКАЯ М. И., ПАВЛОВ А. Н., ЧЕТВЕРУШКИН Б. Н. Об одной схеме интегрирования уравнений газовой динамики. М., 1983 (Препр. / АН СССР. ИПМ им. М. В. Келдыша; №113).
[7] ПОТКИН В. А. Кинетический анализ разностных схем для газовой динамики // Журн. вычисл. математики и мат. физики. 1975. Т. 15, №6. С. 1492-1498.
[8] PULLIN D.I. Direct simulation methods for compressible in viscous ideal gas //J. Comput. Phys. Vol. 34, No. 2. P. 231-244.
[9] REITZ R. D. One-dimensional compressible gas dynamic calculations using the Boltzmann equation // Ibid. 1981. Vol. 42, No. 1. P. 108-132.
[10] Аристов В. В., Черемисин Ф. Г. Решение уравнений Эйлера и Навье —Стокса на основе операторного расщепления кинетического уравнени // Докл. АН СССР. 1983. Т. 272, №3. С. 555-559.
[11] Елизарова Т. Г., ЧЕТВЕРУШКИН Б.Н. Использование кинетических моделей для расчета газодинамических течений. М., 1984 (Препр. / АН СССР. ИПМ им. М. В. Келдыша; №165).
[12] Елизарова Т. Г., ЧЕТВЕРУШКИН Б.Н. Об одном вычислительном алгоритме для расчета газодинамических течений // Докл. АН СССР. 1984. Т. 279, №1. С. 80-83.
[13] Елизарова Т. Г., Павлов А. Н., ЧЕТВЕРУШКИН Б. Н. Применение кинетического алгоритма для расчета газодинамических течений // Дифференц. уравнения. 1985. Т. 21, №7. С. 1179-1185.
[14] Граур И. А., Елизарова Т. Г., Четверушкин Б. Н. Использование кинетических алгоритмов для расчета газодинамических задач, моделирующих вязкие течения. М., 1985 (Препр. / АН СССР. ИПМ им. М. В. Келдыша; №165).
[15] Елизарова Т. Г., Четверушкин Б. Н. Кинетические алгоритмы для расчета газодинамических течений // Журн. вычисл. математики и мат. физики. 1985. Т. 25, №10. С. 1526-1533.
[16] Елизарова Т. Г., Четверушкин Б. Н. Использование кинетических моделей для расчета газодинамических течений. Математическое моделирование. Процессы в нелинейных средах. М.: Наука, 1986. С. 261-278.
[17] Четверушкин Б.Н., Чурбанова Н. Г. Консервативные кинетические схемы для решения задач газовой динамики. М., 1986 (Препр. / АН СССР. ИПМ им. М.В. Келдыша; №78).
[18] Граур И. А., Елизарова Т. Г., Четверушкин Б. Н. Моделирование сложных газодинамических течений на основе кинетических алгоритмов // Дифференц. уравнения. 1986. Т. 22, №7. С. 1173-1180.
[19] DESHPANDE S. M. On the maxwellian distribution symmetric form and entropy conservation for Euler equations // NASA Technical Paper. 1986. No. 2583.
[20] DESHPANDE S. M. Kinetic theory based new upwind methods for in viscous compressible flow // AIAA Paper. 1986. No. 86-0275.
[21] DESHPANDE S. M. Kinetic flux splitting schemes. Comput-Dynamic Review, Wiley, Chichester, 1995. P. 161-181.
[22] SAURO S. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond. Oxford Sci. Publ. Clarendon Press, 2001. 288 p.
[23] Климонтович Ю. Л. Статистическая теория открытых систем. М.: Янус, 1995.
[24] Алексеев Б. В. Обобщенная больцмановская физическая кинетика // Теплофизика высоких температур. 1997. Т. 35, №1. С. 129-146.
[25] ONATE E. Derivation of stabilized equations for numerical solution of advective-diffusive transport and fliud flow problems // Comput. Math. Appl. Mech. Eng. 1998. Vol. 151. P. 233-265.
[26] Onate E., Manzan M. Stabilization Techniques for Finite Element Analysis of Convection-diffusion Problems // Publication CIMNE, N183, Barcelona, Feb. 2000.
[27] Лукшин А. В., Четверушкин Б. Н. К теории кинетически согласованных разностных схем // Мат. моделирование. 1995. Т. 7, №11. С. 109-125.
[28] ВЛАДИМИРОВ В. С. Математические задачи односкоростной теории переноса частиц // Тр. Мат. ин-та АН СССР. 1961. Т. 61. 158 с.
[29] ЛАДЫЖЕНСКАЯ О. А. Математические вопросы динамики вязкой несжимаемой жидкости. М.: Наука, 1970. 288 с.
[30] ТЕМАМ R. Navier — Stokes Equations. Theory and Numerical Analysis. Amsterdam, North-Holland, 1974.
[31] Елизарова Т. Г., Павлов А. Н., Четверушкин Б. Н. Использование квазигазодинамической системы уравнений для расчета обтекания тела с иглой // Докл. АН СССР. 1987. Т. 292, №2. С. 327-331.
[32] Антонов А. Н., Елизарова Т. Г., Четверушкин Б.Н., Шеретов Ю.В. Численное моделирование пульсационных режимов при сверхзвуковом обтекании полого цилиндра // Журн. вычисл. математики и мат. физики. 1990. Т. 30, №4. С. 548-556.
[33] Антонов М.А., Граур И. А., Косарев Л. В., Четверушкин Б. Н. Численное моделирование пульсаций давления в трехмерных выемках // Мат. моделирование. 1996. Т. 8, №5. С. 76-90.
[34] Антонов А. Н., Антонов М.А., Граур И. А. и др. Численное влияние формы выемки на характеристики пульсационного течения газа // Инженерно-физ. журн. 1998. Т. 71, №3. С. 447-483.
[35] Антонов А. Н., Антонов М. А., Граур И. А. и др. Взаимодействие сверхзвукового потока с упругой поверхностью обтекаемого тела // Мат. моделирование. 1998. Т. 10, №2. С. 15-24.
[36] ANTONOV A.N., KOZUBSKAYA Т.К. Acoustic noise simulation for supersonic viscous compressible gas flows // Comput. Fluid' Dynamics'98, CEDS. K.D. Papailioni. Elsevier, Amsterdam, 1998. P. 145-200.
[37] Дородницын Л. В., Четверушкин Б. Н. Кинетически согласованные схемы для моделирования течений вязкого газа // Журн. вычисл. математики и мат. физики. 2000. Т. 40, №12. С. 1875-1889.
[38] Alexandrov A. B., KOZUBSKAYA Т.К., ChetveruSHKIN B.N. Numerical investigation of viscous compressible gas flows by means of flow fields exposure to acoustic radiation // Parallel Computation Fluids Dynamics. Trends and Application / Ed. C. B. Jenssen. Elseveir, 2001. P. 267-273.
[39] Абалакин И. В., Четверушкин Б. Н. Кинетически согласованные разностные схемы как модель для описания течений умеренно разреженных газов // Мат. моделирование. 1993. Т. 5, №5.
[40] Белоцерковская М. С., Трапезникова М. А., Четверушкин Б. Н. Аналог кинетически согласованных схем для моделирования задач фильтрации // Мат. моделирование. 2002. Т. 14, №10.
[41] Бахвалов Н. С., Панасенко Г. П. Осреднение процессов в периодических средах. М.: Наука, 1984. 352 с.
[42] Гиршфельдер Д. Ж., Кертис Ч., Берд Р. Молекулярная теория газов и жидкостей. М.: Иностр. лит-ра, 1961. 930 с.
[43] Дородницын Л. В., Четверушкин Б.Н., Чурбанова Н. Г. Кинетически согласованные разностные схемы и квазигазодинамическая система для моделирования течений плотных газов и жидкостей // Мат. моделирование. 2001. Т. 13, №4. С. 45-57.
[44] CHETVERUSKIN B. N., CHURBANOVA N. G. Simulation of low Mach number flows using the equations dynamic system // XV Intern. Conf. on Chem. Reactors (Chem Reactor-15), Helsinki, Finland, 2001. Abstracts. P. 319-321.
[45] Абалакин И. В., Четверушкин Б. Н. О расширении возможности газодинамического описания с помощью кинетически согласованных разностных схем // Мат. моделирование. 1994. №7. С. 3-14.
[46] Дородницын А. В., Четверушкин Б. Н. Об одном варианте квазигазодинамической системы уравнений // Там же. 2002. Т. 14 (в печати).
[47] Граур И. А., Дородницын А. В., Елизарова Т. Г., Четверушкин Б. Н. Кинетически согласованные схемы газовой динамики с неполной коррекцией. М., 1987 (Препр. / АН СССР. ИПМ им. М. В. Келдыша; №5).
[48] Елизарова Т. Г., Четверушкин Б. Н. Кинетически согласованные разностные схемы для моделирования течений вязкого газа // Журн. вычисл. математики и мат. физики. 1988. Т. 28, №11. С. 695-710.
[49] Самарский А. А., Николаев е. С. Решение систем сеточных уравнений. М.: Наука, 1976.
[50] Шильников Е. В., Шумков М.А. Моделирование трехмерных нестационарных течений газа на МВС с распределенной памятью // Мат. моделирование. 2001. Т. 13, №4. С. 35-46.
[51] CHETVERUSKIN B.N., Shilnikov E.V., SHOUMKOV M. A. Using massively parallel computer systems for numerical simulation of 3D viscous gas flows. Parallel Computional Fluids Dynamics. Trends an Application / Ed. C. B. Jenssen et al. 2001. Elsevier. P. 509-515.
[52] Антонов А. Н., Вишняков А. Н., Шагаев С. Г. Экспериментальное исследование пульсаций давления в выемке, обтекаемой дозвуковым или сверхзвуковым потоком газа // ПМТФ. 1981. №2. С. 89-98.
[53] Антонов А. Н., Купцов В. М., Комаров В. В. Пульсации давления при струйных и отрывных течениях. М: Машиностроение, 1990. 271 с.
[54] Arnal M., LAUREL O., Liek Z. Peric M in AGAMM Workshop // Notes on Numer. Fluid Mech. Vol. 36 / Eds M. Devill, T-Hle, Yu. Morchois (Braunschwieg).
[55] ЯкобовскиЙ М. В. Распределенные системы и сети. М.: МГТУ "Станкин", 2000. 118 с.
[56] Белоцерковский О. М. Прямое численное моделирование свободно развитой турбулентности // Журн. вычисл. математики и мат. физики. 1985. Т. 25. С. 1857-1882.
[57] Ландау Л. Д., Лифшиц е. М. Гидродинамика. Т. VI: Теоретическая физика. М.: Наука, 1988.
[58] Тхир А. В. Метод продвинутого фронта для построения двумерных неструктурированных сеток // Численные методы и приложения / Под ред. Ю. А. Кузнецова. М.: ИВМ РАН, 1995. С. 151-160. Деп. в ВИНИТИ №394-395.
[59] Абалакин И. В., Жохова А. В., Четверушкин Б. Н. Кинетически согласованные схемы повышенного порядка точности // Мат. моделирование. 2001. Т. 13, №5. С. 53-61.
[60] Fiedler M. A property of eigenvectors of nonnegative symmetric matrices and its application to graph theory // Chechoslovak. Math. J. 1975. Vol. 25. P. 619-633.
[61] CHETVERUSHKIN B. N. Kinetic shemes and their application for simulation of viscous gas dynamic problems // Proc. Europ. Conf. of Comput. Methods in Applied Sci. and Eng. (ECCOMAS) CD-ROM.
[62] Дородницын А. В., Корнилина М.А., Четверушкин Б.Н., Якобов-ский М. В. Моделирование газовых течений при наличии химически активных компонентов // Журн. физ. химии. 1997. Т. 71, №12. С. 2275-2281.
Поступила в редакцию 14 января 2002 г.
