ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2010, том 50, № 12, с. 2275-2284
УДК 519.634
К столетию со дня рождения академика А.А. Дородницына
МОДЕЛИРОВАНИЕ ТЕЧЕНИЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ И СЛАБОСЖИМАЕМОГО ГАЗА НА МНОГОЯДЕРНЫХ ГИБРИДНЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ^
© 2010 г. А. А. Давыдов*, Б. Н. Четверушкин* **, Е. В. Шильников*
(* 125047Москва, Миусская пл., 4, ИПМим. Келдыша, РАН;
** 141700Долгопрудный, М.о., Институтский пер., 9, МФТИ) e-mail: alexander.a.davydov@gmail.com; chetver@imamod.ru; shiva@imamod.ru Поступила в редакцию 12.04.2010 г.
Рассматривается логически простой алгоритм на основе явных схем для моделирования течений несжимаемой жидкости и слабо сжимаемого газа. В качестве исходной математической модели используется гиперболический вариант квазигазодинамической системы. В качестве инструмента расчета используется оригинальный вычислительный кластер на основе графических плат NVIDIA. Библ. 16. Фиг. 4. Табл. 2.
Ключевые слова: кинетические схемы, многоядерные вычислительные системы, графические платы, устойчивость явных схем.
ВВЕДЕНИЕ
В настоящее время исследователям стали доступны вычислительные системы с производительностью более 100 TFlops. Это открывает новые большие возможности в математическом моделировании, позволяя моделировать сложные процессы как можно ближе к их реальному протеканию. В первую очередь, это относится к междисциплинарным индустриальным приложениям, которые требуют высокой точности решения в связи с оптимизацией технологических процессов. Имеются, однако, существенные проблемы, связанные с реальным использованием возможностей таких систем. Эти трудности связаны как со спецификой вычислительных алгоритмов для высокопроизводительных вычислительных систем, так и с проблемами вспомогательного программного инструментария (см. [1]). На решение этих проблем направлены усилия специалистов ведущих индустриальных стран.
В настоящее время усиливается тенденция к использованию вычислительных систем, базирующихся на существенно многоядерных процессорах нетрадиционной архитектуры. Это открывает возможности для заметного снижения стоимости системы на единицу производительности и, что не менее важно, многократно снижает потребление электроэнергии, необходимой для обеспечения работы таких высокопроизводительных компьютеров. В свою очередь, это упрощает решение проблемы охлаждения системы. Однако трудности эффективного использования таких систем значительно возрастают даже по сравнению с непростыми проблемами использования уже существующих высокопроизводительных компьютеров.
Многоядерные многопроцессорные вычислительные системы требуют разработки специального математического обеспечения, учитывающего гибридную структуру памяти — распределенную между процессорами и общую для ядер одного процессора. Не менее важным для таких систем являются особо жесткие требования к вычислительным алгоритмам, которые должны быть как можно проще и прозрачнее с логической точки зрения и вместе с тем достаточно эффективны.
В этом смысле очень привлекательными оказываются явные схемы, которые легко адаптируются к различной архитектуре ЭВМ и позволяют эффективно использовать вычислительные си-
Работа выполнена при финансовой поддержке РФФИ (коды проектов 09-01-00600, 08-01-00644).
2275
стемы, содержащие 103—104 ядер. Однако явные схемы накладывают жесткое ограничение на шаг по времени, связанное с требованием устойчивости, особенно при решении уравнений и систем параболического типа, в связи с чем перспективным направлением является разработка явных схем с как можно более мягким условием устойчивости. Таким свойством обладают кинетические схемы (КС) для решения задач газовой динамики, о которых пойдет речь в данной работе. Кроме того, в работе будет описан оригинальный вычислительный кластер, который позволяет расширить возможности графических плат для решения задач, описываемых уравнениями математической физики.
Одному из авторов (Б.Н. Четверушкину) особенно хотелось обратить внимание на тот интерес, который в конце 80-х годов академик А.А. Дородницын проявил к новому направлению в вычислительной гидродинамике — кинетическим схемам. Его поддержка явилась важным стимулом для дальнейшего развития этого направления.
1. КВАЗИГАЗОДИНАМИЧЕСКАЯ СИСТЕМА И КИНЕТИЧЕСКИЕ СХЕМЫ
Кинетические, или больцмановские, схемы (см. [2]—[5]) в настоящее время широко применяются в вычислительной газовой динамике и являются эффективным подходом к моделированию задач механики сплошной среды на многопроцессорных вычислительных системах. Эти схемы отличаются от других алгоритмов, в первую очередь, тем, что их вывод основывается на некоторой дискретной модели одночастичной функции распределения. Например (см. [2]), предположим, что в некоторый момент времени I = Р функция распределения является локально максвел-ловской/0(Р, х, %) и слабо меняется на расстояниях порядка длины свободного пробега молекул I. Здесь х = {х1, х2, х3} — пространственные координаты, % = £,2, £,3} — компоненты скорости молекулы газа. В течение интервала времени между столкновениями молекул т = Р +1 — Р происходит их бесстолкновительный разлет, который может быть описан соотношением
/(г1+\ X, %) = / (гX - %т, %). (1.1)
В момент времени I = Р +1 происходит мгновенная максвеллизация функции распределения. Принимая во внимание малость изменения^ на расстояниях порядка I, разложим правую часть (1.1) в ряд Тейлора. В результате получим балансное соотношение с точностью до членов второго порядка 0(т2) малости по т:
/ + 1 - /- + = А^к! /. (1.2)
-л * V*/
т дх, о%1 2 дхк
Интегрируя равенство (1.2) по скоростям молекул с сумматорными инвариантами ф(%) = {1, %, %2/2}, используя основное предположение процедуры Чепмэна—Энскога (%)/й% = (%)/0й% и
заменяя разделенную разность в левой части (1.2) на производную по времени, получаем квазигазодинамическую (КГД) систему уравнений. Для трехмерного течения идеального политропно-го газа эта система, записанная в стандартных обозначениях в прямоугольных декартовых координатах, в консервативной форме имеет следующий вид:
д + = о, (1.3)
дг д х1
д , д . , д и т-г /1 „ч
и + + Тр = ^-П, (1.4)
дг дх, дх, дх ,
д . (Е + р\ д „ _ д
д г дх- V р ) дх - - д х ,
Е + + ^ = -^Пи, а5)
Здесь и далее по повторяющимся индексам подразумевается суммирование, Е = р(и2/2 + б) — полная энергия. Поток массы задается как
= Р(и1- Щ), Щ = т(^-Риуи;- + ^-р) . (1.6)
рчдх- д х1 ;
Тензор вязких напряжений П, тепловой поток q и соотношения, замыкающие систему, определяются следующими соотношениями:
тт ттш , ( д ,1 д V с- ( д , д ^ ¡л
Щ. = Щ + хри^ик — и + р —р) + Т.икдХР + УРдХик) , (1/7)
NS , ( д д 2 о д
Щ = +. - з 8.дХкик), (1-8)
NS ( д , д А NS д г Пч
4 = 4 - три;иД —6 + р— -), ъ = -к— Т, (1.9)
д х. р) дх;-
р = рб(у - 1), Т = Р. (1.10)
рЯ
Динамическая вязкость ц, коэффициент теплопроводности к и релаксационный параметр т, имеющий размерность времени, имеют следующий вид:
Ц = Ц»Iт) , к = ^р , т = "Й", (1.11)
V Т») (у - 1) Рг рБс
где у — показатель адиабаты, Рг и 8е — числа Прандтля и Шмидта соответственно. Различные способы вывода КГД-системы уравнений, а также большое количество задач, решенных с применением этой системы, можно найти в [2], [6], [7].
КГД-система (1.3)—(1.5) получена из балансного уравнения (1.2) в предположении малости времени т по сравнению с характерным газодинамическим временем ?газ. Учитывая следующий член разложения функции распределения в ряд Тейлора, приходим к такому варианту системы:
тд!р + дР + А]т1 = 0, (1.12)
2-^2 д/ д х д / ил/^
2
тд ри, д , д . , д дп /л
+ -р и, + —]ти + —р = — П., (1.13)
2 д/ д / дх. дх, дх.
тд2Е д ^ , д . (Е + р^ , д д ,-г 1/|Ч
- т + н;Е + —) + тъ = . (1.14)
2 д/2 д/ дх. ^ р ) дх. дх.
Эта система имеет гиперболический тип (см. [8]). Важным свойством КГД-системы уравнений (1.12)—(1.14) является ее близость к классическим уравнениям Навье—Стокса (НСУ). Связь между этими двумя системами может быть представлена в виде
КГД = НСУ + 0(Кп2), (1.15)
где Кп — число Кнудсена. Следует отметить, что сами уравнения Навье—Стокса получены из уравнения Больцмана с точностью до членов второго порядка малости по числу Кнудсена. Дополнительные диссипативные члены в КГД-системе, содержащие множитель т, обращаются в нуль в тех областях течения, где решение удовлетворяет стационарным уравнениям Эйлера. Заметим также, что наличие вторых производных по времени в уравнениях (1.12)—(1.14) окажется важным фактором при формировании вычислительного алгоритма, описанного в разд. 3.
Конечно-разностная аппроксимация систем (1.3)—(1.5) и (1.12)—(1.14) строится на основе метода контрольных объемов. Все газодинамические параметры относятся к центрам ячеек разностной сетки, а в качестве контрольного объема берется сама ячейка. Интегрируя уравнения (1.3)—(1.5) по объему ячейки, мы получаем законы сохранения в интегральной форме. При этом изменение газодинамических параметров в ячейке определяется суммой потоков консервативных переменных (плотности р, импульса рu и полной энергии Е) через все ее грани. Для аппроксимации пространственных производных, входящих в выражения для потоков, используются центральные разности, а значения газодинамических переменных в центрах граней вычисляются с помощью линейной интерполяции. Для обеспечения устойчивости счета системы (1.3)—(1.5) по
явной схеме к релаксационному параметру т в (1.11) добавляется слагаемое, пропорциональное шагу пространственной сетки:
т = -Н- + а -, (1.16)
pSc c
где c — локальная скорость звука, h = Jax^ + Ax2 + Ax3, , Ax1, Ax2, Ax3 — шаги прямоугольной пространственной сетки, а — числовой параметр порядка единицы, подбираемый экспериментально.
Можно показать (см. [2]), что если истинные вязкие члены малы по сравнению с конвективными, то дополнительные члены дискретного происхождения в (1.3)—(1.5) играют роль искусственной вязкости, используемой при решении уравнений Эйлера для идеального газа. Если же истинные вязкие члены одного порядка с конвективными, то члены сеточного происхождения малы по сравнению и с теми, и с др
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.