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

Информационные технологии Вестник Нижегородского университета им. Н.И. Лобачевского, 2013, № 1 (3), с. 301-305

УДК 541.186

РАЗРАБОТКА ДВУМЕРНОГО CFD-КОДА ДЛЯ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ГАЗОДИНАМИЧЕСКИХ ТЕЧЕНИЙ

© 2013 г. Р.В. Жалнин, О.Е. Каледин, Е.Е. Пескова

Мордовский госуниверситет им. Н.П. Огарева, Саранск

kaledinoe@mrsu.ru

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

Рассматривается решение двумерной системы уравнений Навье-Стокса методом конечных объемов. Приведен алгоритм для численного решения системы уравнений Навье-Стокса, алгоритм декомпозиции расчетной области для параллельной реализации. Работа алгоритма демонстрируется на примере обтекания профиля КЛСЛ0012.

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

Введение

Развитие средств вычислительной газовой динамики и компьютерной техники делает возможным разработку и реализацию методов расчета нестационарных течений жидкости и газа в пространственных областях сложной конфигурации. При использовании неструктурированных сеток отпадает необходимость применения многоблочных сеток в многосвязных областях, а также появляется возможность добавления новых узлов без перестройки всей сетки, что требуется для точного расчета пограничных слоев, скачков уплотнения и вихревых структур [1, 2].

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

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

ния профиля КЛСЛ0012 потоком вязкого газа [3].

1. Дискретизация уравнения Навье-Стокса

Рассмотрим двумерную систему уравнений Навье-Стокса, записанную в консервативной форме

д,и+ У-Ре (и)-Ру (и,^и)=0,

&и д¥е (и) -V- (и, Уи)

— +-------—-----------—------ +

& дх

дОе (и) -дОу (и, Уи) п +------—---------—------- = 0,

ду

(1)

и дополненную подходящими начальнокраевыми условиями. Консервативные переменные и и декартовы компоненты (и), (и)

невязкой (эйлеровой) потоковой функции ¥е (и) заданы в виде

' р ры и = - !>,

ру Е

(2)

где р — плотность жидкости; и, V — компоненты скорости, р — давление и

Е = е + [ы2 + V2) /2) — полная энергия.

ры ру

=- ры2 + р , =- рыу

рыу ру1 + р

(Е + р)ы (Е + р)у

Работа выполнена при поддержке Федеральной целевой программы «Научные и научно-педагогические кадры инновационной России» на 2009-2013 гг. (контракт № 14.B37.21.0176).

Полагая, что жидкость подчиняется произвольным уравнениям состояния, р можно записать как р = р(р, е) . В настоящей работе используется уравнение состояния идеального газа

р = е (у-1).

Декартовы компоненты функции (и, Уи),

(и, Уи) вязкого потока (и, Уи) заданы в

виде

Fv (u, Vu) =

Gv (u, Vu) =

О

т

Tx

uT + vt О

Tyx Tyy ut + VT

yx

xy

(3)

Тхх =(Л~ 2И/3)(Ых + ^ ) + 2^Ых ,

ТУУ =(Л~ 2^/3)(Ых + ^ ) + 2ИЫу , (4)

т = т = и(ы + V ),

ху ух ' \ у х/

где л — коэффициент динамической вязкости, X — объёмная вязкость, ц- (/, - = х, у) — компоненты тензора вязких напряжений. В настоящей работе не рассматривается влияние теплопроводности, к = 0. Проинтегрировав по области О, мы получим окончательную формулировку:

J д,udQ +j^V • F(u,Vu)dQ = :Х|7пд( ud Q + ^V^ F(u, Vu)d Q

= О,

(5)

Vv,

где Р(и, Уи) = ^ (и) - ^ (и, Уи) и интегралы по области О были преобразованы в сумму интегралов над совокупностью непересекающихся элементов {Т}, которые будем считать треугольниками. Интеграл по границе контрольного объема вычисляется как сумма произведений значений вектора потока в центрах граней контрольного объема на площади его граней.

Проинтегрировав по частям каждую элементарную {Т} уравнения (5), которое содержит дивергенцию функции потока Навье-Стокса, мы получим ослабленную формулировку

2 [| г(и,Уи)-ш/о--

т

-{ГУ-Г(11,У11>Ю] = 0, (6)

Уу,

где 0Т обозначает границу элемента Т.

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

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

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

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

C помощью программы Triangle сгенерирована треугольная сетка [4]. Имеется 3 файла: *.node - содержит узлы сетки, *.ele - содержит номера треугольников и координаты их вершин; *.neigh - содержит номера треугольников и их соседей. Для параллельной реализации алгоритма необходимо разбить построенную сетку на N частей (N - количество процессов). Этот вопрос можно свести к задаче разбиения графа. Целью разбиения является получение независимых подмножеств вершин исходного графа с минимальным количеством связей между этими подмножествами.

Построим матрицу смежности для исходной сетки. Единицы используем, если ячейки являются соседними; если ячейки не связаны - то ставим ноль. Получаем сильно разреженную матрицу (табл. 1). Шобходимо так перенумеровать её узлы, чтобы как можно ближе к главной диагонали концентрировались ненулевые элементы (табл. 2). Таким образом, при разбиении

сетки мы уменьшим количество связей между полученными подмножествами. Данные преобразования проводились с помощью алгоритма Катхилла-Макки [5].

Таблица 1

0 І 2 З 4 5 б n-І n

0 І 0 0 І 0 0 0 І І І

І 0 І 0 І 0 0 І І 0 І

2 І 0 І 0 І 0 0 І І 0

З І 0 І І 0 І 0 І 0 0

4 0 І 0 І І 0 0 І 0 І

5 І 0 І 0 І І 0 І 0 0

б 0 І 0 І 0 0 І І 0 І

0 0 І 0 І 0 І І 0 І

n-І І 0 0 І 0 І 0 І І 0

n 0 І 0 0 І 0 І І 0 І

0 І 2 З 4 5 б Таблица 2 n-І n

0 І І І І І 0 0 0 0 0

І І І І І І 0 0 0 0 0

2 0 І І І І І 0 0 0 0

З 0 0 І І І І І 0 0 0

4 0 0 0 І І І І І 0 0

5 0 0 0 І І І І І 0 0

б 0 0 0 І І І І І 0 0

0 0 0 0 І І І І І 0

n-І 0 0 0 0 0 І І І І І

n 0 0 0 0 0 І І І І І

После перенумерации узлов было проведено разбиение полученной сетки на N частей (табл. 3). Для каждой части созданы файлы с расширениями*.еіе, *.пеі§Ь, нумерация ячеек

вкоторых начинается с нуля. Также для каждого подмножества были созданы файлы: *^е^ -содержит номера ячеек, значения которых необходимо отправить, номер процесса-получателя, глобальный номер ячейки; *.гееу -содержит номера ячеек, значения в которых получает текущий процесс, номер процесса-отправителя, глобальный номер ячейки.

Таблица З

0 І 2 З 4 5 б n-І n

0 І І І І І 0 0 0 0 0

І І І І І І 0 0 0 0 0

2 0 І І І І І 0 0 0 0

З 0 0 І І І І І 0 0 0

4 0 0 0 І І І І 1 І 0 0

5 0 0 0 І І І І І 0 0

б 0 0 0 І І І І І 0 0

0 0 0 0 І І І І І 0

n-І 0 0 0 0 0 І І І І І

n 0 0 0 0 0 І І І І І

На рис. 1 изображен результат работы модуля декомпозиции для N=8.

Рис. І. Декомпозиция расчетной области

З.Обтекание профиля NACA0012

Для проверки полученного разбиения построен параллельный алгоритм решения уравнений газовой динамики на треугольной сетке. Работа программного комплекса тестировалась на кластере со следующими характеристиками: 4 вычислительных узла, 16 Gb Ram, 2*Intel Xeon E5410 общей производительностью 0,3 терафлопс.

На рис. 2-5 приведены результаты работы программы для задачи обтекания профиля NACA0012 с начальными данными [3, 6]:

1. M = 1 - число Маха,

2. T = 273 К - температура,

3. P = 101325 Па - давление,

4. Alpha = 0 и 5° - угол атаки,

5. f = 0,0000185 Па/c - коэффициент динамической вязкости.

Рис. 2. Поле давления, Alpha = 0, t = 0,0З c

Рис. 4. Поле давления, Alpha = 5°, t = 0,03 c

Рис. 5. Магнитуда скорости, Alpha = 5°, t = 0,03 c

Зависимость ускорения от количества процессоров приведена в таблице 4.

Таблица 1

Результаты вычислительных экспериментов

Максимальное ускорение более чем в 20 раз было получено для параллельного алгоритма при размере сетки 171406 ячеек.

Разработки по проекту соответствуют выполняемым в последнее время научнообразовательным программам по развитию су-перкомпьютерного образования в стране [7-10].

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

1. Волков К.Н., Емельянов В.Н. Течения и теплообмен в каналах и вращающихся плоскостях. М.:Физматлит. 2010. 486 с.

2. Куликовский А.Г., Погорелов Н.В., Семенов

A.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит. 2001. 608 с.

3. AGARD Fluid Dynamics Panel. Test cases for inviscid flow field methods // AGARD Advisory Reports. 1986. № AR-211.

4. A two-dimensional quality mesh generator and delaunay triangulator. [Электронный ресурс] режим доступа: http://www. cs. cmu. edu/~quake/triangle.html

5. Дом проектов SkyFox. Блог Дмитрия Тимофеева [Электронный ресурс] Режим доступа: http://digital-space.ru/fizmat/algoritm-katxilla-makki/

6. Holst T.L. Computational Fluid Dynamics. Drag Prediction-Results fromthe Viscous Transonic Airfoil Workshop. Ames Research Center, Moffett Field, California, 1988.

7. Гергель В.П., Стронгин Р.Г. Опыт Нижегородского университета по подготовке специалистов в области суперкомпьютерных технологий // Вестник Нижегородского университета им. Н.И. Лобачевского. 2010. № 3-1. С. 191-199.

8. Воеводин В.В., Гергель В.П. Суперкомпью-терное образование: третья составляющая суперком-пьютерных технологий // Вычислительные методы и программирование: новые вычислительные технологии. 2010. Т. 11, № 2. С. 117-122.

9. Об опыте проведения программ повышения квалификации профессорско-преподавательского состава по направлению «Высокопроизводительные вычисления» / В.П. Гергель, А.В. Линев, И.Б. Мее-ров, А.В. Сысоев // Открытое и дистанционное образование. 2010. № 3. С. 15-20.

10. Развитие системы суперкомпьютерного образования в России: текущие результаты и перспективы / В.В. Воеводин, В.П. Гергель, Л.Б. Соколинский,

B.П. Демкин, Н.Н. Попова, А.В. Бухановский // Вестник Нижегородского университета им. Н.И. Лобачевского. 2012. № 4. С. 268-274.

Размер сетки hQ S та р во ог дл сПо ны Параллельный алгоритм

8 процессоров 32 процессора

Время Ускорение Время Ускорение

98389 ячеек 0,589 0,076 7,75 0,031 19

171406 ячеек 1,092 0,144 7,58 0,049 22,28

DEVELOPMENT OF TWO-DIMENSIONAL CFD-CODE FOR NUMERICAL MODELING

OF GAS-DYNAMIC

R. V. Zhalnin, O.E. Kaledin, E.E. Peskova

In this paper we consider the solution of two-dimensional Navier-Stokes equations by finite volume. Describes an algorithm for the numerical solution of the Navier-Stokes equations, shows a diagram of the computational domain decomposition for parallel implementation. The algorithm is demonstrated on the example of the flow airfoil NACA0012.

Keywords: Navier-Stokes equations, the computational domain, decomposition, finite volume method.