CC BY
77
УДК 532.529.519
Е. А. Данилкин, А. В. Старченко
Томский государственный университет пр. Ленина, 36, Томск, 634050, Россия
E-mail: [email protected], [email protected]
ПАРАЛЛЕЛЬНАЯ РЕАЛИЗАЦИЯ ЧИСЛЕННОГО МЕТОДА РЕШЕНИЯ СИСТЕМЫ УРАВНЕНИЙ НАВЬЕ-СТОКСА ПРИ МОДЕЛИРОВАНИИ КРУПНЫХ ВИХРЕЙ ТУРБУЛЕНТНЫХ ТЕЧЕНИЙ *
В работе представлена параллельная реализация алгоритма численного решения системы уравнений Навье-Стокса при моделировании турбулентности методом крупных вихрей. Для подсеточного моделирования применялась модель Смагоринского, которая в сочетании со схемой Ван Лира позволила получить хорошее согласование с экспериментальными данными. Исследуются различные способы геометрической декомпозиции при численном решении уравнения переноса, проведен теоретический анализ их эффективности и разработаны рекомендации по их использованию. В работе большое внимание уделяется повышению производительности написанных программ, при этом использовались такие приемы повышения эффективности, как порядок вложенности циклов, учет эффектов кэш-памяти.
Ключевые слова: декомпозиция расчетной области, параллельные вычисления, моделирование турбулентности, вихреразрешающее моделирование.
Введение
Будем предполагать, что турбулентное движение несжимаемой жидкости описывается уравнениями Навье-Стокса
du, dut 1 dp 2 /in
47 + Uj ^ = — + v VU, 1 = ^,3, (1)
dt dXj p dxi
dr=°' (2)
dxi
где ut = ut (X, t) - компоненты мгновенного поля скорости; v - коэффициент кинематической вязкости; p - плотность; p - мгновенное значение давления; V2 = d2/ dXj dxj - оператор Лапласа. По повторяющимся индексам проводится суммирование.
В настоящее время существует три основных и часто используемых подхода моделирования турбулентных течений. Первый подход - это прямое численное моделирование (Direct Numerical Simulation, DNS). Под прямым численным моделированием подразумевается численное решение системы дифференциальных уравнений (1)-(2), аппроксимированных ко-нечными-разностями высокого порядка. DNS обладает тем достоинством, что не требует каких-либо турбулентных замыканий, а результаты, полученные при моделировании, используются подобно экспериментальным данным. Прямое численное моделирование пригодно для исследования внутренних механизмов турбулентности и может выступать как инструмент для детального изучения мелкомасштабной турбулентности [Курбацкий, 2000]. Недостатками прямого численного моделирования являются сложность выбора начальных условий и высокая стоимость проведения расчетов на ЭВМ, т. е. необходимо использовать суперкомпьютеры с массовым параллелизмом, время счета на которых исчисляется десятками часов, а расчеты ограничены течениями при относительно низких числах Рейнольдса.
Один из наиболее дорогих с точки зрения вычислительных затрат расчетов турбулентного потока в канале потребовал приблизительно 120 суток непрерывного счета на 2048 процессорах массивно-параллельного компьютера при числе Рейнольдса равного 2003. За это время было промоделировано 0,8 секунды движения воздуха в канале [Hoyas, Jimenes, 2005].
* Работа выполнена при финансовой поддержке РФФИ (проект № 07-05-01126) и программы «СКИФ_ГРИД».
ISSN 1818-7900. Вестник НГУ. Серия: Информационные технологии. 2009. Том 7, выпуск 2 © Е. А. Данилкин, А. В. Старченко, 2009
Такая высокая затратность прямого численного моделирования является результатом природы турбулентности, поскольку для того, чтобы точно предсказывать турбулентную структуру течения, необходимо моделировать весь спектр масштабов турбулентности Колмогорова [Курбацкий, 2000]. Поэтому маловероятно, что прямое численное моделирование в самом ближайшем будущем будет регулярно применяться в инженерных расчетах.
Второй подход - моделирование крупных вихрей (Large Eddy Simulation, LES) - использует более грубую сетку по сравнению с предыдущей методикой. Он основан на концепции «фильтрования» турбулентности, т. е. осреднении уравнений Навье-Стокса по пространству с некоторым весом - «фильтром» [Там же]. Так как посредством грубой сетки нельзя моделировать наименьшие масштабы турбулентности, то используются модели подсеточных масштабов для учета мелкомасштабной диссипации энергии турбулентности и обратного рассеяния энергии потока от малых к большим масштабам [Турбулентные..., 1982]. Вместе с этим крупномасштабные турбулентные движения рассчитываются явно и не требуют моделирования. Этот метод более точен в сравнении с подходом, основанным на использовании осредненных по времени уравнений Навье-Стокса, в котором моделируется одновременно весь спектр масштабов турбулентности. Вихреразрешающие модели впервые были разработаны для предсказания погодных явлений [Smagorinsky, 1963], однако вскоре повсеместно стали использоваться для моделирования турбулентных течений.
Относительно высокая расчетная стоимость все еще ограничивает применение вихрераз-решающих моделей для предсказания разномасштабных процессов в окружающей среде. Поэтому в рамках этого подхода существует острая необходимость создания высокопроизводительных решателей с привлечением суперкомпьютеров. Неоспоримым достоинством LES-моделей является способность предсказывать нестационарную турбулентную структуру течения и различные турбулентные корреляции, поэтому иногда для получения такой информации данный подход применяют при моделировании сложного турбулентного движения воздушных масс в уличном каньоне [Walton et al., 2002].
Третий подход базируется на осредненных по времени уравнениях Навье-Стокса (Reynolds-Averaged Navier-Stokes, RANS). Здесь турбулентное течение рассматривается как сумма двух компонент: осредненной части и пульсационной [Курбацкий, 2000]. Среднее течение рассчитывается с использованием осредненных уравнений Навье-Стокса, которые включают дополнительные нелинейные члены, содержащие напряжения Рейнольдса u'u'j (u' - это
пульсация скорости, а черта сверху обозначает осреднение Рейнольдса). Для неизвестных напряжений Рейнольдса также может быть выведено транспортное уравнение, однако оно
будет включать дополнительные неизвестные третьего порядка u'U' u'k и т. д.. Эта задача определения напряжений известна в моделировании турбулентности как «проблема замыкания». Существуют различные модели турбулентности, которые аппроксимируют uUuU , начиная от простых алгебраических соотношений и заканчивая теми, которые включают дополнительные транспортные уравнения для каждой из шести независимых компонент тензора напряжений Рейнольдса. Данный подход моделирования турбулентности требует значительно меньше вычислительных ресурсов компьютера по сравнению с двумя предыдущими.
Хотя в изложенном выше материале все подходы моделирования турбулентности были для простоты разбиты на три категории, на самом деле существует множество гибридных методов, в которых, для примера, у твердых стенок используются нестационарные осреднен-ные по времени уравнения Навье-Стокса, а в основном потоке расчеты осуществляются по вихреразрешающей модели. Обсуждение этих гибридных методов (Detached Eddy Simulation (DES), Very Large Eddy Simulation (VLES) и Unsteady RANS (URANS)) приводится в работе [Spalart, 2000].
Постановка задачи
В данной работе моделирование турбулентного течения осуществляется с использованием метода крупных вихрей. В турбулентном потоке образуются вихри различных масштабов, взаимодействующие и обменивающиеся энергией друг с другом, что обусловлено нелинейностью уравнений Навье-Стокса. Крупные вихри являются энергосодержащими, а мелкие
вихри ответственны за процесс вязкой диссипации. При численном решении этих уравнений необходимо описывать весь спектр турбулентности, при этом разрешающая способность модели должна позволять описывать ответственные за диссипацию вихри.
Основная идея метода LES заключается в формальном математическом разделении крупных и мелких структур. Крупные вихри разрешаются явно путем численного моделирования, а мелкомасштабная турбулентность параметризуется, то есть определяется характеристиками крупномасштабных движений. Возможности к этому открывает теория универсального равновесия Колмогорова [Лойцянский, 2003], в соответствии с которой если нет возможности для разрешения всех масштабов турбулентного движения при численном моделировании, то следует моделировать мелкомасштабные вихри как приближенно изотропные структуры по сравнению с анизотропными крупномасштабными вихрями .
Чаще всего для разделении крупных и мелких структур используется операция фильтрации [Турбулентные..., 1982; Wilcox, 1993], согласно которой переменные ячеечного масштаба определяются уравнением
f (x, t) = J G(x, x') f (x', t)dx', J G(x)dx = 1.
2ч
(3)
Здесь G - функция фильтра с характерным масштабом длины А. Отфильтровывая уравнения Навье-Стокса (1), (2), получаем
dU, 5й;й: дт— 1 dp d2и, —- +-1 + —— =---i- +-idt dx— dx— p dx, dx—dx—
^ = 0, (4)
dxi
т p = u,uj - . (5)
Система уравнений (3)-(5) остается незамкнутой, так как (5) содержит исходные нефильтрованные компоненты скорости. Тензор (5) характеризует влияние мелкомасштабных вихрей на эволюцию крупномасштабных вихрей, и его необходимо моделировать как и в подходе Буссинеска через установление его связи со средними скоростями u. Последнее и составляет суть методики подсеточного моделирования.
Выбор подсеточной модели. Существует большое количество подходов моделирования подсеточных масштабов, но наиболее широкое распространение получили модели, основанные на использовании турбулентной вязкости (eddy viscosity models, EVM) для не разрешенных явно турбулентных движений. В таких моделях тензор напряжений вычисляется по формуле
т а = -2 KTSif, S а =
у 1 у7 J 2
1 f a— tt
dui + диу dx. dxf
J 1
Здесь - тензор скорости деформации, построенный по фильтрованному полю скорости
и, Кт = Кт (и, х, t) > 0 - коэффициент турбулентной вязкости, зависящий от решения. Основным свойством ЕУМ является их способность параметризовать прямой энергетический каскад от крубных вихрей к мелкомасштабным, а недостатком является низкая корреляция рассчитанных компонент тензора напряжений с наблюдаемыми значениями этих компонент [БаМта й а1., 1980].
Выбор зависимости Кт = Кт (и, х, t) чрезвычайно разнообразен. Смагоринский [Smagorin-sky, 1963] предложил влияние мелкомасштабных вихрей на эволюцию крупномасштабных вихрей аппроксимировать выражением
т Г =-2С2 ^^^ ,
где А = И - шаг сетки модели (если используется неравномерная сетка, то А = (ИхИуИ2 )1/3), |б| = ^2БуБу - норма тензора скорости деформации, С3- постоянная Смагоринского.
R
R
Классическая модель Смагоринского обладает рядом достоинств, оправдывающих ее применение в LES:
- в работах Лилли была найдена теоретическая оценка параметра Смагоринского Cs ~ 0,17, связанная с постоянной Колмогорова;
- в работе [Bardina et al., 1980] на основании априорного анализа DNS было показано, что модель Смагоринского хорошо аппроксимирует сток энергии в области подсеточных масштабов.
К недостаткам классической модели Смагоринского с постоянным коэффициентом Cs следует отнести:
- завышенную диссипацию в областях с ламинарным потоком и в областях, в которых турбулентность анизотропна;
- низкую корреляцию между наблюдаемым и параметризуемым тензором напряжений, присущую всем моделям, основанным на введении эффективной турбулентной вязкости, вследствие несовпадения ориентации тензора турбулентных напряжений и тензора деформации.
Граничные условия. Поскольку пространственное разрешение не позволяет описывать процессы в вязком приповерхностном слое, то суммарное воздействие мелкомасштабных неоднородностей стенки должно быть учтено при помощи той или иной пристеночной модели.
Большинство пристеночных моделей основано на введении буферного слоя некоторой толщины, в котором применяется RANS модель. Известен ряд работ, в которых удалось достичь правильного поведения средних величин в LES-моделях, ограничившись простейшей пристеночной моделью, не выходящей за пределы первого расчетного слоя.
Чтобы смоделировать граничные условия в пристеночной расчетной ячейке, которая располагается внутри логарифмической области пограничного слоя, используется следующий закон стенки:
U = iln EUn
uT к v
где к - это постоянная Кармана, E = 9,793, uT - динамическая скорость, п - расстояние от стенки.
Аппроксимация и метод решения. Аппроксимация дифференциальной задачи осуществляется на основе метода конечного объема. Основная идея этого метода заключается в разбиении расчетной области на непересекающиеся, граничащие друг с другом конечные объемы так, чтобы один узел расчетной сетки содержался только в своем конечном объеме. Значения компонент скорости определяются на гранях конечных объемов, а скалярные характеристики - в центре. Разбив таким образом расчетную область, интегрируем каждое уравнение по каждому конечному объему. При вычислении интегралов применяется кусочно-полиномиальная интерполяция для зависимых величин. Аппроксимация конвективных членов уравнений выполняется с использованием противопотоковой схемы MLU Ван Лира [Van Leer, 1974]. Для решения уравнений движения применяется явная схема, поэтому результатом приближенного интегрирования уравнения движения по одному конечному объему является готовая для вычислений формула
ap°Ф j = aPi,j,кФП j,к + ae, ,кФ^,,к + an, ¿Ф¡ 1,к + atu к,м +
(6)
+ aw,,j,кфП_1,j,к + as, j кфПj_1,к + ab,,j,кфП,,к-1 + Ъ, j,к.
Здесь Ф j - обобщенная скалярная величина, которая может представлять любую компоненту скорости. Полученная разностная схема имеет первый порядок аппроксимации по времени и второй по пространству, является условно устойчивой и монотонной.
При расчете течений в областях сложной геометрии использовался метод фиктивных областей, суть которого заключается в том, что значения векторных и скалярных величин в области преграды равны нулю и на границах фиктивных конечных объемов отсутствуют диффузионные потоки.
Согласование полей скорости и давления в гидродинамической части модели выполняется с помощью подхода коррекции поля давления. Для решения уравнения Пуассона для поправки давления применяется метод неполной факторизации Булеева [Данилкин, Старченко, 2006].
Параллельная реализация
Сравнение различных способов распараллеливания рассмотренного явного алгоритма решения нестационарных уравнений Навье-Стокса с целью отыскания оптимального с точки зрения минимизации затрат на пересылку проводилось на примере одного адвективно-диффузионного уравнения.
Математическая постановка. Требуется решить трехмерное уравнение переноса в единичном кубе (х,у,z) е (0,1)х (0,1)х (0,1)
дФ д(иФ) д(уф) д(^Ф) д
д дх ду дz дх
Y— д гдф д гдф"
dx + dy _ dy _ + dz _ dz _
+S,t>0
дФ « дФ = 0; дФ
= 0; —
dx x=0 дУ y=0 "дГ z=0
дФ =0; дФ = 0- дФ
dx x=1 ' дУ = 0; y=i "дТ z=1
dx dy dz
с граничными условиями второго рода и заданными начальными условиями
= 0 Ф1 t=0 = 0 = 0.
Функции u, v, w, Г > 0 и S - определены в рассматриваемой области. Эта задача решается численно представленным выше методом.
Методика распараллеливания. В качестве основного подхода распараллеливания выбрана геометрическая декомпозиция сеточной области. В рассматриваемом случае возможны три различных способа разделения значений сеточной функции по вычислительным узлам - одномерная или ленточная схема, двухмерное или блочное разбиение или трехмерное разбиение узлов вычислительной сетки.
Каждому процессорному элементу вместе с выделенной сеточной подобластью распределяются все значения сеточной функции Ф" j k, принадлежащие этой подобласти [Данилкин, 2008].
После этапа декомпозиции, когда производится разделение данных на блоки для построения параллельного алгоритма, переходим к этапу установления связей между блоками, вычисления в которых будут выполняться параллельно, - планированию коммуникаций. Из-за используемого шаблона явной разностной схемы для вычисления очередного приближения в приграничных узлах каждой подобласти необходимо знать значения сеточной функции с соседнего граничащего процессорного элемента. Для этого на каждом вычислительном узле создаются фиктивные ячейки для хранения данных с соседнего вычислительного узла и организуются пересылки этих граничных значений, необходимых для обеспечения однородности вычислений по явным формулам (6) [Там же] (рис. 1). Пересылка данных осуществляется с использованием процедур библиотеки MPI [Беликов и др., 2008].
Обратимся к предварительному теоретическому анализу эффективности различных способов декомпозиции расчетной области для рассматриваемого случая. Будем оценивать время работы параллельной программы как время работы последовательной программы Tcalc, деленное на количество используемых процессоров, плюс время пересылок: Tp = Tcalc / p + Tcom. Время пересылок для различных способов декомпозиции можно приблизительно выразить через количество пересылаемых данных: TlD = t d • 2N2 х 2 , T2D = t d • 2N2 x 4 / p112,
com send ' com send ± ?
где N3 - размерность конечно-разностной задачи, p tsend - время пересылки одного числа.
T3D = t
end 2N2 x 6/ p2'J, (7)
количество вычислительных узлов,
Ш декомпозиция
2Т> декомпозиция
Фштпинын
Организация обменов на примере Ш
декомпозиции
2 пэ
3Б декомпозиция
Рис. 1. Различные способы декомпозиции. Схема организации обменов для 1Б декомпозиции.
Из формулы (7) видно, что для различных способов декомпозиции затраты на пересылку данных можно представить в виде Тсот = tsend ■ 2Ы2 х к(р), где коэффициент пропорциональности к(р) зависит от способа декомпозиции и числа используемых процессорных элементов.
В табл. 1 представлены числовые значения к(р) . Видно, что при р > 5 более эффективными будут 2Б и 3Б декомпозиции, а при р > 11 в случае 3Б декомпозиции потребуется пересылать между процессорными элементами меньшее количество сеточных значений функции Ф,"++1к, и можно ожидать, что в этом случае затрачиваемое на пересылку время будет минимально.
Таблица 1
Зависимость коэффициента пропорциональности к (р) от числа процессорных элементов и способа декомпозиции
Число процессов 3 4 5 6 10 11 12 16 60 120 250
Ш декомпозиция 2 2 2 2 2 2 2 2 2 2 2
2Б декомпозиция 2,31 2,00 1,79 1,63 1,26 1,20 1,15 1,00 0,51 0,36 0,25
3Б декомпозиция 2,88 2,38 2,05 1,82 1,29 1,21 1,14 0,94 0,39 0,24 0,15
Анализ результатов. Расчеты производились на кластерной системе ТГУ СКИФ СуЬейа на сетках 240 х 240 х 240 и 360 х 360 х 360 при использовании до 144 процессоров. Результаты вычислительного эксперимента показали наличие хорошего ускорения при решении задач данного класса. Основное внимание уделялось сравнению времени пересылок и времени вычислений при различных способах декомпозиции.
Сетка 240x240x240 Сетка 360x360x360
0 30 60 90 120 150 0 30 60 90 120 150
Число процессоров Число процессоров
—•—1Р—В-2Р-*-3Р —Ю —■- 2Э - -3Р
Рис. 2. Ускорение для различных способов декомпозиции расчетной области
На первом этапе использовалась одна общая программа, размеры массивов от запуска к запуску не менялись, на каждом процессорном элементе нумерация элементов массивов начиналась с нуля. Несмотря на то, что в соответствии с теоретическим анализом 3Б декомпозиция является оптимальным вариантом распараллеливания, вычислительные эксперименты показали, что лучших результатов можно достичь, используя 2Б декомпозицию при числе используемых процессов от 25 до 144 (рис. 2).
Рис. 3. Время вычислений без учета затрат на пересылку данных для различных способов декомпозиции
Исходя из предварительного теоретического анализа на графиках должны присутствовать следующие закономерности. Время вычислений без учета межпроцессорных коммуникационных затрат при разных способах декомпозиции должно примерно совпадать для одинакового числа процессоров и уменьшаться как 1сс1к / р. В реальности же расчетные данные (рис. 3) указывают на то, что использование 2Б декомпозиции на различных сетках дает минимальные затраты на проведение вычислений и расчетные графики зависимости времени вычислений от числа используемых процессоров размещаются существенно выше, чем ТсаЬ / р .
Для объяснения полученных результатов необходимо обратить внимание на допущения, которые были сделаны при предварительном теоретическом анализе поставленной задачи.
Во-первых, предполагалось, что независимо от способа распределения данных на одном процессорном элементе выполняется одинаковый объем вычислительной работы, который должен приводить к идентичным временным затратам. Во-вторых, принималось, что время, затраченное на межпроцессорную пересылку любой последовательности одинакового количества данных, не зависит от способа их выборки из оперативной памяти.
Для того чтобы разобраться, что же происходит в реальности, был проведен следующий набор тестовых расчетов. Для оценки состоятельности первого допущения был рассмотрен случай, когда программа запускалась в однопроцессорном варианте, и при этом имитировались различные способы геометрической декомпозиции данных при одинаковом объеме вычислений, выполняемом каждым процессором. Было реализовано две версии программы: в первой использовались статические массивы размером N, во второй - «динамические», когда размеры массивов специально подбирались в соответствии со способом декомпозиции для того, чтобы минимизировать временные затраты, связанные с выборкой данных из оперативной памяти.
Таблица 2
Влияние способа декомпозиции сеточной области и размерности используемых массивов на время вычислений на одном вычислительном узле
Разме рность задачи 3603
Статические массивы 360 х 360 х 360 «Динамические» массивы с учетом декомпозиции
Р 1D 2D 3D Р 1D 2D 3D
8 9,3 8,1 10,4 8 8,8 7,5 7,3
24 3,1 2,6 3,7 24 3,0 2,5 2,5
64 1,2 1,0 1,6 64 1,2 0,9 1,0
В табл. 2 представлено время в секундах, которое требуется процессору для проведения вычислений в своей сеточной подобласти, при моделировании различных способов геометрической декомпозиции. От расчета к расчету объем вычислительной работы сохраняется, меняется лишь характер расположения обрабатываемых данных в оперативной памяти.
Проанализировав представленные результаты, можно сделать вывод, что при написании параллельной программы целесообразно использовать «динамические» массивы с подстраиваемыми под выделенное число процессоров размерами. В большей степени это существенно для 3D декомпозиции. Преимущества использования «динамических» массивов, на наш взгляд, определяется меньшими временными затратами на выборку обрабатываемых данных из оперативной памяти и их передачу через кэш-память.
Для оценки реализуемости второго допущения ( Tom = tsend • 2N2 х k(p) ) был рассмотрен тест, в котором измерялось время MPI-пересылки различных сечений трехмерного массива. Заметим, что решение этой задачи (пересылки, скажем, i-j, i-k, j-k - сечений массива, элементы которого имеют индексы i, j, k) в стандарте MPI может осуществляться различными способами в зависимости от того, какие сечения массивов рассматриваются.
Вычислительный эксперимент был организован следующим образом. В параллельной программе многократно выполнялась пересылка сечения массива одного базового типа между двумя процессорными элементами и фиксировалась продолжительность этого процесса. Использовались функции MPI_SEND и MPI_RECV и процедуры создания производных типов MPI_TYPE_VECTOR и MPI_TYPE_INDEXED. В тестовых расчетах рассматривались массивы с количеством элементов 303, 603, 1203. Полученные результаты показали, что и при организации межпроцессорной передачи данных также решительное преимущество имеет использование «динамических» массивов, поскольку в этом случае, допускается применение производного типа MPI_TYPE_VECTOR, который для пересылки сечения j-k в несколько раз по времени предпочтительней, чем более универсальный тип MPI_TYPE_INDEXED.
Таблица 3
Время на пересылку различных сечений трехмерного массива (с)
Количество элементов в массиве / сечение i-j i-k j-k
303 0,00189 0,00192 0,00265
603 0,00402 0,00556 0,00884
1203 0,01601 0,02151 0,03563
В табл. 3 представлено время, затрачиваемое на 100-кратную пересылку данных из разных сечений массива, осуществляемую с помощью производного типа МР1_ТУРБ_УБСТОЯ. Анализируя данные, можно отметить, что время пересылки сечений массивов в МР1-БОЯТКЛК-программе зависит от их типа. Поэтому наименьшие затраты на передачу данных дает Ш декомпозиция, наибольшие - 3Б декомпозиция, причем это отличие становится более существенным при увеличении размера сечения от 30 х 30 до 120 х 120. Последнее связанно с увеличением времени на выборку передаваемых данных из оперативной памяти. Таким образом, в формулах (7) необходимо вводить дополнительные коэффициенты, которые бы учитывали тип пересылаемых сечений массивов сеточных функций.
Сетка 360x360x360
0 30 60 90 120 150
Число процессоров —*— 1D —■— 2D - rtr - 3D
Рис. 4. График ускорения при использовании «динамических» массивов
На рис. 4 представлен график ускорения параллельной программы, использующей «динамические» массивы, для различных способов декомпозиции. Здесь так же, как и на рис. 2, наилучшие показатели демонстрирует 2D декомпозиция, причем имеет место прирост ускорения. Незначительно ей уступает 3D декомпозиция, проблемы которой в маштабируемости параллельной программы связаны, главным образом, с большими, чем в 1D и 2D случаях, затратами на выборку, обработку и пересылку данных из оперативной памяти.
Таким образом, для явных разностных методов решения уравнения переноса возможно применение одномерной, двумерной и трехмерной декомпозиции, однако результаты тестирования программ показали, что 3D декомпозиция не дает выигрыша во времени по сравнению с 2D декомпозицией, по крайней мере, для количества процессоров, не превышающего 250, при этом 3D декомпозиция отличается более трудоемкой программной реализацией, а использование 2D декомпозиции является достаточным для масштабирования данной задачи при рассматриваемом числе вычислительных узлов.
Результаты тестирования численного метода и LES-модели
Рассмотрим турбулентный поток вокруг цилиндра квадратного сечения, который расположен в канале (рис. 5). Выбор области течения и значения числа Рейнольдса Re = UnD/v соответствуют эксперименту [Lyn et al., 1995]. Хотя представленный случай имеет довольно простую геометрию, за цилиндром образуется неустойчивый след. В расчетах использовалась следующая сетка: Nx1 х Nx2 х Nx3 = 132 х 122 х 22.
Рис. 5. Расчетная область
Движение жидкости поддерживалось за счет задания постоянной скорости на входной границе, направленной вдоль оси канала. В качестве начальных данных использовались постоянные значения скорости йх = Um = 0,55 м / с, U2 = U3 = 0. Расчеты (и осреднение величин) проводились до восьми периодов отрыва вихря после времени инициализации, равного 50D / Um .
Расчеты показали хороший уровень согласования с экспериментальными данными [Ibid.]. На рис. 6 приведены результаты сравнения применения различных подсеточных моделей для описания мелкомасштабной турбулентности, взятые из [Nakayama, Vengadesan, 2002] (в этой работе для аппроксимации конвективных слагаемых использовалась центрально-разностная схема). На рис. 7 представлены расчеты, полученные в данной работе, где для моделирования подсеточных масштабов использовалась модель Смагоринского.
Рис. 6. Сравнение различных подсеточных моделей
Рис. 7. Результат, полученный в данной работе
Сравнение результатов (рис. 6 и 7) показывает, что схема Ван Лира вносит свой вклад в подсеточное моделирование и позволяет достичь более точного согласования результатов численного эксперимента с измерениями [Lyn et al., 1995] для средних характеристик потока.
Для оценки пульсационных составляющих использовалась следующая формула:
uf =<ut >2 - < ui > , которая выводится из утверждения < ui >= < ui > + u'.
0,3
0,4
0,2
О Expti-Lyn st äL 1995 Smagorrsky .----NorrtotW -----Dynamic model — One ei|uatiim modsl
/ _n№nDD„ \ ^^^m l -зЛ i-'iüin^-- Ш/f ' Ш/ -f'""* ic Д ff"JTji □ ¡j u Ii В п..........!..,,
-10 12 3 4 5 6 x,fD
Рис. 8. Турбулентные напряжения для пульсаций скорости и1, справа - результат, полученный в данной работе
oje
0,6
0.4
1111 II 111 1 1 | 1 Г Г I ; м 1 1 ; м м | м м 0 Lyn е'. Gl. (*995)
л -Smajorinsky
t ----Nö model
i -----Dynamic
I f -Л — One equaüon
1 А : f 1 ч ч ^ \а \ л» -
: im \ V,
! p 4—
. V D v
' J П
■V' /D . '-r fo
1 JJ : W / /
/
/
✓
в.
И 0 1 г S 4 6 $
Рис. 9. Турбулентные напряжения для пульсаций скорости u2 , справа - результат, полученный в данной работе
Из рис. 8 и 9 видно, что получено достаточно хорошее соответствие с данными эксперимента [Lyn et al., 1995] для относительных пульсаций компонент скорости u1,u2. Отсутствует характерное для расчетов [Nakayama, Vengadesan, 2002] по модели Смагоринского серьезное завышение изменчивости u12 / U2n и недооценка изменчивости u^ / U2n вдоль оси Ох1.
Рис. 10. Ускорение параллельной программы для моделирования крупных вихрей турбулентных течений
На рис. 10 представлено ускорение параллельной программы, полученное при численном решении системы уравнений Навье-Стокса (3)-(4) с использованием подсеточной модели Смагоринского на сетке 132 х 122 х 22 и двумерной декомпозиции области.
Заключение
Результаты численных экспериментов показали, что построенная математическая модель турбулентности способна воспроизводить характерные особенности турбулентного потока. Использование модели Смагоринского в сочетание со схемой Ван Лира позволило получить хорошее согласование средних величин потока с данными эксперимента для исследуемой области.
Применение в расчетах 2D декомпозиции дает 76 % эффективности при использовании 25 вычислительных узлов. При дальнейшем увеличении количества вычислительных узлов до 100 при выбранных размерах сетки, получено характерное для задач данного класса значение эффективности около 50 %.
Список литературы
Беликов Д. А., Говязов И. В., Данилкин Е. А., Лаева В. И., Проханов С. А., Старченко А. В. Высокопроизводительные вычисления на кластерах. Томск: ТГУ, 2008. 198 с.
Данилкин Е. А., Старченко А. В. К выбору способа декомпозиции при численном решении систем связанных дифференциальных уравнений на многопроцессорной технике с распределенной памятью // Третья Сибирская школа-семинар по параллельным вычислениям. Томск: ТГУ, 2006. С. 95-101.
Данилкин Е. А. К вопросу об эффективности 3D декомпозиции при численном решении уравнения переноса с использованием МВС с распределенной памятью // Вестник ТГУ. Томск: ООО «Издательство научно-технической литературы», 2008. № 2. С. 39-46.
Курбацкий А. Ф. Лекции по турбулентности: В 2 ч. Новосибирск, 2000.
Лойцянский Л. Г. Механика жидкости и газа: Учеб. для вузов. 7-е изд., испр. М.: Дрофа. 2003. 840 с.
Турбулентные сдвиговые течения: Пер. с англ.; ред. А. С. Гиневский. М.: Машиностроение, 1982. Т. 1. 432 с.
Bardina J., Ferziger J. H., Reynolds W. C. Improved Subgrid Scale Models for Large-Eddy Simulation // Am. Inst. Astronaut., 1980. P. 80-135.
Hoyas S., Jimenez J. Scaling of the Velocity fluctuations in turbulent Channels up to Re = 2003 // Annual Research Briefs, Center for Turbulence Research, NASA Ames/Stanford Univ. 2005. P.351-356.
Lyn D., Einav S., Rodi W. et al. A Laser-Doppler Velocimetry Study of Ensemble Averaged Characteristics of the Turbulent Near Wake of a Square Cylinder // J. Fluid Mech. 1995. Vol. 304. P. 285-319.
Nakayama A., Vengadesan S. N. On the Influence of Numerical Schemes and Subgrid-Stress Models on Large Eddy Simulation of Turbulent Flow Past a Square Cylinder // Int. J. Numer. Methods Fluids. 2002. Vol. 38. P. 227-253.
Smagorinsky J. General Circulation Experiments with the Primitive Equations. I: The Basic Experiment // Monthly Weather Review. 1963. Vol. 91. No. 3. P. 99-165.
Spalart P. R. Strategies for Turbulence Modelling and Simulations // Int. J. of Heat and Fluid Flow. 2000. Vol. 21. P. 252-263.
Van Leer B. Towards the Ultimate Conservative Difference Scheme. II. Monotonicity and Conservation Combined in a Second Order Scheme // J. of Computational Physics. 1974. Vol. 14. P. 361-370.
Walton A., Cheng A. Y. S., Yeung W. C. Large-Eddy Simulation of Pollution Dispersion in an Urban Street Canyon. Part 1: Comparison with field data // Atmospheric Environment. 2002. Vol. 36. P. 3601-3613.
Wilcox D. C. Turbulence Modeling for CFD: DCW Industries Inc, 1993. P. 460.
Материал поступил в редколлегию 15.02.2009
E. A. Danilkin, A. V. Starchenko
PARALLEL NUMERICAL METHOD OF SOLUTION OF NAVIER-STOKES EQUATIONS FOR LARGE EDDY
SIMULATION OF TURBULENT FLOWS
Parallel implementation of algorithm of numerical solution of Navier-Stokes equations for large eddy simulation (LES) of turbulence is presented in this research. The Smagorinsky model is applied for sub-grid simulation of turbulence. This model with numerical Van Leer's scheme for advection terms provides good agreement with experimental data. Various ways of geometrical decomposition for parallel numerical solution of transport equations are investigated. A theoretical analysis of parallel algorithms effectiveness was performed and recommendations for their application were developed. Special techniques like the order of nested loops and the effect of cash memory are applied to increase performance of the developed parallel programs.
Keywords: domain decomposition, parallel computations, turbulence modelling, LES approach.
