ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2008, том 48, < 2, с. 309-320
УДК 519.634
ВЫЧИСЛИТЕЛЬНЫЙ И ПРОГРАММНЫЙ ИНСТРУМЕНТАРИЙ ДЛЯ МОДЕЛИРОВАНИЯ ТРЕХМЕРНЫХ ТЕЧЕНИЙ ВЯЗКОГО ГАЗА НА МНОГОПРОЦЕССОРНЫХ СИСТЕМАХ1)
© 2008 г. Б. Н. Четверушкин, Е. В. Шильников
(125047 Москва, Миусская пл., 4А, ИММ РАН) e-mail: shiva@imamod.ru, chetver@imamod.ru Поступила в редакцию 20.07.2007 г.
Предлагается параллельный программный комплекс, предназначенный для численного моделирования трехмерных задач вязкого газа. В качестве основы вычислительного алгоритма используются кинетически согласованные разностные схемы на локально сгущающихся сетках. Программный комплекс апробирован на ряде модельных задач о сверх- и дозвуковом обтекании различных тел. Отмечена возможность прямого моделирования турбулентных течений. Исследована эффективность параллельной реализации в зависимости от размера задачи и от числа используемых процессов. Библ. 16. Фиг. 7. Табл. 1.
Ключевые слова: кинетические схемы, газовая динамика, локально сгущающиеся сетки, параллельные вычисления, программный комплекс.
ВВЕДЕНИЕ
Существенно нестационарные и турбулентные режимы течения вязкого сжимаемого газа привлекают все более пристальное внимание исследователей, обусловленное, в частности, важностью нестационарных процессов, возникающих при решении широкого спектра индустриальных задач. В случае симметричной геометрии при не слишком больших числах Рейнольдса вполне успешным оказывается численное моделирование таких задач на основе двумерных нестационарных уравнений Навье-Стокса. При больших числах Рейнольдса, которые, как правило, и реализуются на практике, на квазипериодическое двумерное нестационарное течение накладываются существенно трехмерные стохастические турбулентные флуктуации. Таким образом, численное моделирование должно проводиться в трехмерной постановке даже для простой симметричной геометрии течения. Моделирование детальной структуры трехмерных течений вязкого газа является сложным и дорогостоящим и возможно только с применением современных многопроцессорных вычислительных комплексов, бурный рост производительности которых отмечается в последние годы. Производительность самой мощной в настоящее время вычислительной системы составляет 360 TFlops. В свою очередь, производительность систем, находящихся вблизи нижней границы списка 500 самых высокопроизводительных систем в мире (ТОР-500) превосходит 2 TFlops. Имеются, однако, существенные проблемы, связанные с реальным использованием возможностей таких систем. Эти трудности связаны как со спецификой вычислительных алгоритмов для высокопроизводительных вычислительных систем, так и с проблемами вспомогательного программного инструментария (см. [1]). На решение этих проблем направлены усилия специалистов ведущих индустриальных стран.
В данной работе рассматривается один из возможных подходов к решению пространственно трехмерных задач динамики вязкого газа. В качестве основы вычислительного алгоритма используются кинетические схемы (см. [1]), которые вместе с родственными им Lattice Boltzmann-схемами (LBS) (см. [2], [3]) становятся одним из эффективных механизмов моделирования задач механики сплошной среды на многопроцессорных вычислительных системах. Для пространственной аппроксимации используются структурированные локально сгущающиеся сетки. При этом в значительной мере использовались наработки, полученные при решении двумерных и трехмерных задач (см. [4]-[6]). Следует отметить, что структурированные сетки вполне конкурентоспособны при моделировании ряда задач с неструктурированными сетками, которые наряду с преимуществами при описании сложных конфигураций требуют решения специальных проблем, связанных с их построением и аппроксимацией на них соответствующих уравнений (см. [7]-[8]).
1)Работа выполнена при финансовой поддержке РФФИ (коды проектов 06-01-00187, 05-07-90230, 05-01-00510).
1. ОБЩИЕ СВЕДЕНИЯ О ПРОГРАММНОМ КОМПЛЕКСЕ
Разработка предлагаемого комплекса велась на алгоритмическом языке Fortran-77. Параллельные модули ППП реализованы на основе функций из MPI-библиотек (MPICH 1.2), что обеспечивает высокую переносимость комплекса. Распараллеливание исходной задачи производится на основе принципа геометрического параллелизма, суть которого для задач математической физики заключается в следующем. Область, в которой ищется численное решение задачи, разбивается на подобласти с примерно одинаковым числом счетных сеточных узлов в каждой, а общее количество подобластей равно числу процессоров. Если используемые схемы однородны, то этот подход автоматически обеспечивает равномерную загрузку процессоров.
Прозрачность логики программы весьма желательна из-за намного большей трудности обнаружения и анализа логических ошибок в распределенных программах по сравнению с последовательными. Кроме того, логическая простота программы значительно облегчает ее перенос с одной многопроцессорной ЭВМ на другую. Этот фактор наиболее заметен в случае, если установки имеют не полностью идентичную архитектуру.
Главный принцип, принятый в качестве основополагающего, базировался на идее, что все сложные, но не обязательно параллельные обработки следует вынести в предваряющие (подготовительные) последовательные программы. В связи с этим автоматически определилась структура комплекса. Процесс обработки различных данных, требующихся для моделирования изучаемой задачи на многопроцессорной ЭВМ, состоит из четырех относительно самостоятельных по своим функциям стадий:
1) описание геометрии исследуемого объекта, построение разностной сетки, задание граничных условий;
2) разбиение сеточной области на подобласти, каждая из которых будет рассчитываться на отдельном процессоре, окончательная подготовка к расчету;
3) собственно расчет на многопроцессорной вычислительной установке;
4) визуализация полученных результатов.
При этом с необходимостью параллельной является только третья стадия. Две первые вполне могут быть выполнены на однопроцессорной ЭВМ. Что касается четвертой стадии, то здесь ситуация менее однозначна. Построению визуализационной картинки на экране персонального компьютера предшествует выбор соответствующих данных из хранилища информации, полученной в результате высокопроизводительных вычислений. Этот выбор может осуществляться в многопроцессорном режиме, основываясь на принципе геометрического параллелизма. При этом каждый процессор обрабатывает ту часть данных, которая была им получена в процессе расчетов. Такой подход к визуализации данных высокопроизводительных вычислений на основе технологии "клиент - сервер" в настоящее время достаточно широко распространен (см. [9]-[10]). Однако если объемы данных не слишком велики, а скорость построения картинки не является критичной для нас, то и с этой задачей может справиться однопроцессорное приложение.
Несмотря на относительную самостоятельность указанных выше задач, они все-таки взаимосвязаны, причем определяющей является третья, а точнее - тип выбранного численного метода. Поэтому в первую очередь остановимся на этом этапе.
2. ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ УРАВНЕНИЙ ГАЗОВОЙ ДИНАМИКИ
В качестве основы для численного моделирования использовалась квазигазодинамическая система уравнений (КГУ), которая может быть интерпретирована как дифференциальное приближение кинетически согласованных разностных схем (КСРС, см. [1]). При этом основой для получения системы КГУ является балансное уравнение для одночастичной функции распределения f(t, x, X):
f1 + 1 - f M + = д £,£kAtдf m (2 1)
A t ^ д xi д xi 2 д xk' (
Здесь x = {x1, x2, x3} - пространственные координаты, X = - скорость молекулы газа,
t - время, At - шаг по времени, j - номер шага по времени, fM - локально-максвелловская функция
распределения
-.2
Р( x, t) f fo-M^O]'
[2пRT(x, t )] 3/2 pI 2 RT(x, t)
/м = - - - - t - , 3/2 exp[ -1 „ZZ':' I. (2.2)
В свою очередь в (2.2) использовались стандартные обозначения: р - плотность, T - температура, u = {m1, u2, u3} - скорость газовой среды.
Система КГУ получается из уравнения (2.1) путем интегрирования последнего по скоростям молекул с компонентами вектора столкновений j(X) = {1, X, X2/2} и с использованием основного предположения процедуры Чепмена-Энскога
J j (X)/dX = J j (X) /м<%. (2.3)
КСРС по своему способу получения близки к Lattice Boltzmann-схемам, которые в качестве своей основы опираются на дискретные модели для функции /. Наиболее популярной для получения LBS является модель Бхатнагара-Гросса-Крука
= v(/м - /), (2.4)
где v - частота столкновений молекул. В свою очередь, балансное уравнение (2.1) может быть получено непосредственно из (2.4) (см. [11]), что еще более подчеркивает тесную связь между КСРС, КГУ и LBS. Связь уравнений (2.1) и (2.4) с физической точки зрения представляется естественной, так как оба уравнения получены в предположении, что после столкновений функция распределения мгновенно становится локально максвелловской. Наличие КГУ системы дает значительные преимущества КСРС по сравнению с LBS. Дифференциальная форма КГУ предоставляет возможность конструирования более широкого спектра алгоритмов.
При At = т, где т - время между столкновениями молекул, получаемая из (2.1) система КГУ связана с уравнениями Навье-Стокса (N-S) следующим соотношением (см. [1]):
КГУ = N-S + 0(т2). (2.5)
Следует отметить, что сами уравнения Навье-Стокса получены из уравнения Больцмана с точностью до членов второго порядка малости. Полученная таким образом система КГУ непосредственно использовалась для моделирования вязких течений. Однако для моделирования сверх- и трансзвуковых течений хорошо зарекомендовали себя кинетические схемы с коррекцией, являющиеся упрощено аппроксимацией исходной системы КГУ. Эти схемы явились основой вычислительного алгоритма, используемого в описываемом программном комплексе. На прямоугольной сетке для случая 3Б-течения идеального политропного газа кинетическая схема с коррекцией выглядит следующим образом:
^Ат + (рu)j ~ 'т-(р"2
+ (рUj)°, = (Tj(рUj + p); )j (2.6)
P Ui P Ui + (pUiUj + p Sj) j = (V Ui)-( 1 + + (Tj (р u i u2
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.