Численное моделирование турбулентного перемешивания с использованием высокопроизводительных систем

Помимо несомненного интереса, явление турбулентного перемешивания, при развитии неустойчивостей, имеет большое значение для ряда практических задач. Создание новых источников энергии, основанных на процессе управляемого термоядерного синтеза, в частности, лазерного термоядерного синтеза (JITC), проектирование сложных технических конструкций, в которых происходит нестационарный нелинейный перенос вещества среды, требует подробного исследования законов переходных и турбулентных процессов.

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

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

К сожалению, развитие программного обеспечения для суперкомпьютеров отстает от общего развития технологий разработки программного обеспечения для персональных компьютеров[94,97]. Это связано с разницей в алгоритмах решения одних и тех же задач для однопроцессорной машины и многопроцессорного вычислительного комплекса. Поэтому тот огромный потенциал программного обеспечения, созданный для решения задач гидродинамики на персональном компьютере полностью должен быть пересмотрен применительно к многопроцессорным комплексам.

Данная диссертационная работа посвящена разработке и созданию программного комплекса для моделирования гидродинамических неустойчивостей и исследования основных характеристик турбулентного перемешивания на МВС с распределенным видом памяти.

В работе рассмотрено 3 вида неустойчивостей. Неустойчивость Рихтмайера — Мешкова [5,6] развивается на контактной границе 2х сред, когда ударная волна проходит из менее плотной среды в более плотную и наоборот. В [7] на основании численных расчетов было показано, что такого типа неустойчивость развивается и в том случае, когда ударная волна проходит через контактную границу двух равноплотных газов, но с различными показателями адиабаты. Второй тип неустойчивостей — неустойчивость КельвинаГельмгольца, возникающая на границе раздела двух сред разной плотности при их скольжении относительно друг друга. И неустойчивость Рэлея-Тейлора [8−10], которая развивается на контактной границе в том случае, когда более плотная среда помещена в поле тяготения над менее плотной средой, либо когда менее плотным слоем пытаются ускорить более плотный слой.

Работы по исследованию неустойчивости Рихтмайера-Мешкова были продиктованы в 60-е годы прикладными вопросами физики высокой плотности энергии. В 80-е годы начались экспериментальные и численные исследования этой задачи [39,40]. Первые расчеты были сделаны по лагранжевым методикам, и это позволило аккуратно отследить динамику развития контактных границ на линейной и нелинейной стадиях развития, но из-за сильного искажения формы лагранжевых ячеек не удавалось промоделировать переходную стадию развития процесса, наблюдаемую в экспериментах [38]. Лишь только качественный скачок в развитии разностных методов конце 90-х годов позволил в дальнейшем проведение численных экспериментов поздней нелинейной стадии развития неустойчивости вплоть до начала процессов турбулентного перемешивания [36,41,43,44].

Как известно, решение различных задач математической физики, которые описываются гиперболическими системами уравнений, могут быть гладкими в одних областях и разрывными в других [56]. Разрывные решения могут при этом возникать из гладких начальных данных. Такие свойства решений налагают на алгоритмы численного решения гиперболических систем уравнений достаточно противоречивые требования. С одной стороны, численный метод должен уметь сохранять свойство монотонности в тех областях, где искомые решения имеют большие перепады значений. С другой стороны, тот же метод должен обладать высоким порядком точности в тех областях, где решение является гладким. В 1959, в работе [53] Годунов показал, что линейных разностных схем, удовлетворяющих этим двум требованиям одновременно быть не может. В данной работе была сформулирована монотонная схема Годунова, имеющая первый порядок и на каждом шагу которой в соседних ячейках решалась задача Римана о распаде разрыва. В дальнейшем это породило целый класс схем, называемыми схемами типа Годунова. К таким схемам относятся также схемы типа TVD. Подробный обзор работ, посвященных численным методам задач гидродинамики, проведен в [47,48].

Создание TVD схем связано с уточнением понятия монотонности схемы для систем линейных дифференциальных уравнений гиперболического типа. Общие идеи и принципы построения сводятся к следующим: невозрастание полной вариации разностного решения (Total Variation Diminishing) схемы.

58,59,62−64] и неувеличение числа локальных экстремумов при ограниченном росте полной вариации (Essentional Non-Oscillatory) [77]. TVD свойство является более слабым требованием, чем монотонность, однако, за счет ослабления требований к схеме, удается построить схемы более высокого порядка, не осциллирующих на разрывах решения. Развитие TVD схем связано с именем Хартена. В работах [60,61,67,69] Хартен ввел класс TVD схем и изучил их свойства, и кроме того, предложил схему 2-го порядка. Для повышения порядка схемы и удовлетворению TVD условию развита специальная техника кусочно-линейной реконструкции сеточных функций. Наклоны кусочно-линейных сеточных функций внутри ячеек для выполнения TVD условия ограничиваются специальными ограничителями (limiters). Они действуют на конечные разности сеточных функций. В работах [73−75] дан анализ современному состоянию дел в области построения ограничителей. Во многих TVD схемах применяется так же и метод энтропийной коррекции, представленный в работе Хартена [68]. В последнее время практически все схемы переменного порядка точности основаны на кусочно-полиномиальной реконструкции дискретных сеточных функций, удовлетворяющих TVD-свойству или свойствам его обобщающим. Это такие схемы как UNO[66,72], TVB [65], ENO [66−70], и др.

Многие схемы не являются в общем случае TVD схемами, так же как и TVD схемы не являются монотонными, а имеют более слабые свойства.

Создание многочисленных разновидностей методик привело к значительному повышению качества получаемых численных решений, и даже вошло в обращение понятие схемы высокого разрешения (high resolution schemes) [71,73−76].

Другим перспективным направлением решения задач гидродинамики является развитие экономичных численных методов. При решении задач математической физики часто приходится иметь дело с изолированными областями высоких градиентов, в которых сложно аппроксимировать решение и которые, однако, могут существенным образом влиять на весь дальнейший процесс. Примерами тому служат задачи динамики жидкости и газа со сложными конфигурациями ударных волн и контактных разрывовзадачи лазерного термоядерного синтеза [52], в которых характерны значительные (на порядки) изменения масштабов длины за время сжатия мишени и т. д. В настоящее время существует достаточно много вычислительных методик для решения задач указанного класса. Среди них можно выделить такие направления как многосеточные итерационные методы, методы подвижных сеток и адаптивные методы. Перспективным и интенсивно развивающимся за последние десять лет является математическое моделирование на основе иерархических локально — измельчающихся адаптивных сеток (adaptive mesh refinement — AMR). Преимущество данного метода перед другими методами заключается в возможности получить решение требуемой точности за минимальные вычислительные затраты. Кроме того, иерархическая структура построения и хранения расчетных сеток предоставляет возможность эффективного использования высокопроизводительных машин с распределенной памятью. Однако, динамическая природа адаптивности при моделировании течений, зависящих от времени, делает достаточно сложным и трудоемким приложение методологии данного типа для современных параллельных суперкомпьютеров. Подробное описание метода AMR было предложено в работе М. Berger и J. 01iger[138]. Развитие метода и приложение его для решения систем газовой динамики описано в работе М. Berger и P. Colella [139] и обобщено на трехмерный случай Bell и др. [140−145]. -В последнее время появилось огромное количество работ, использующих схемы сквозного счета на движущихся адаптивных сетках [83−88]. Использование схем сквозного счета на движущихся адаптивных сетках позволяет придавать численным методам сквозного счета элементы не только высокого, но и практически точного разрешения особенностей течения среды.

Огромное количество работ посвящено экспериментальным и численным исследованиям турбулентного перемешивания при развитии гравитационной неустойчивости. Развитие малых начальных возмущений на границе раздела двух жидкостей приводит к образованию зоны турбулентного перемешивания. Со временем процесс развития такой зоны должен выйти на автомодельный режим, когда характерные размеры структур течения будут пропорциональны ширине самой зоны, а она будет нарастать по закону L (t)=aAgt2, A=(p2-pi)/(Рг+Рь) — число Атвуда р2, рьплотности тяжелой и легкой жидкостей, g — ускорение, t — время [11] Относительно недавно были сформулированы условия реализации автомодельного режима турбулентного перемешивания (Анучина и др. 1978; Youngs 1984 [12,13]) — (I) полное «забывание» начальных возмущений, (II) отсутствие в задаче макромасштаба размерности длины (то есть отсутствие влияния границ области и т. д.), и (III) пренебрежимая малость вклада физических механизмов, ответственных за диссипацию энергии. Удовлетворение всех трех этих требований оказалось чрезвычайно трудным, и экспериментально и в вычислительном отношении.

В настоящее время[14] нет полного согласия относительно выполнения условия (I): можно ли когда-либо полностью «забывать» о начальных условиях (Cook & Dimotakis 2001; Clark 2003) [15,16]. Тем не менее, многие исследователи сходятся во мнении, что начальные условия чувствуются до тех пор, пока ширина зоны перемешивания не становится больше по сравнению с наибольшей длиной волны, присутствующей в начальных возмущениях (Glimm. 2001; Dimonte, 2004) [17,18]. Следовательно, для выполнения этого условия необходимо задавать мелкомасштабные возмущения. В численном моделировании, минимальная длина волны возмущения ограничена шагом по сетке, и разрешающей способностью вычислительной разностной схемой. Относительно условия (II) известно, что процесс 'слипания пузырьков' прекращается, когда доминирующая длина волны, или 'пузырьковый диаметр', вырастает до ширины области. Кроме того, как только количество мод в области уменьшается, статистические измерения ширины зоны перемешивания нельзя считать истинными. Следовательно, в экспериментах желательно делать размер области как можно больше. Относительно (III): можно пренебречь явлениями диффузионных процессов, если происходят расчеты неустойчивости в случае несмешивающихся жидкостей. В работе [14] показано, что, хотя процессы перемешивания смешивающихся жидкостей происходят в масштабах, намного меньших чем длины волн, присутствующие в возмущениях, «то есть пузырьковый диаметр» намного больше, чем Колмогоровский масштаб, тем не менее диффузионные процессы оказывают влияние на развитие крупномасштабных структур. Наконец, для смешивающихся жидкостей диффузионными эффектами нельзя пренебрегать ни для каких значений числа Рейнольдса. A Youngs (1991) [19] отметил, что в экспериментах, использующих несмешивающиеся жидкости (Read 1984) [20] и проводимых для больших чисел Рейнольдса, возможно, полученная оценка значения, а является завышенной применительно для смешивающихся жидкостей.

Многочисленные экспериментальные и численные исследования турбулентного перемешивания посвящены проверке квадратичного закона роста зоны перемешивания и количественной оценке величины а. Определение величины, а имеет важное значение для полуэмпирических моделей развития турбулентности. В различных натурных экспериментах и численных расчетах это значение меняется в довольно широком диапазоне [17,20,22−25].

Эксперименты обладают значительным преимуществом перед численным моделированием в достижение большого числа Рейнольдса. Cook, и Zhou (2002)[21] выполнили прямое численное моделирование РТН используя 535 миллионов ячеек, однако, их числа Рейнольдса (5500) было недостаточно, чтобы выявить значительный инерционный интервал в энергетическом спектре. Однако, намного проще варьировать начальными возмущениями в численной модели, нежели порой даже определить в эксперименте. Эксперименты, включающие ускорение корпуса (КучеренкоRead 1984; Dimonte и Schneider 2000) [26,20,27], обычно не характеризовали начальные возмущения, а эксперименты, включающие выдергивание пластины на границе раздела разноплотных газов (ЗайцевDuff, Harlow и Hirt 1962;

Dalziel, Linden и Youngs 1999) [28- 30], как известно, вводят нежелательные крупномасштабные нарушения течения. Так же и другие условия натурных экспериментов далеки от идеальных условий * автомодельного режима. В них, в частности, присутствует параметр размерности длины, определяемый конечными размерами ускоряемой т капсулы и ограничено время наблюдения развития неустойчивости.

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

Многие вопросы течения турбулентного перемешивания при развитии неустойчивости РэлеяТейлора являются все еще открытыми, хотя значительный скачек в исследованиях в последнее время связан с возможностью моделирования данного явления на высокопроизводительной технике. Подробный обзор ранних работ, посвященных численным и экспериментальным исследованиям, проведен Sharp (1984) [31]. В работе Linden, Redondo и Youngs (1994) [32], сочетающей и экспериментальные и трехмерные численные модели развития неустойчивости, использующие растворимые жидкости, было проведено всестороннее сравнение и получено качественное подобие, но столь же хорошее количественное согласие получить не удалось. Позже, в работе Dalziel, Linden и Youngs (1999) [30] также были проведены численные и натурные эксперименты турбулентного перемешивания между двумя растворимыми жидкостями, порожденного неустойчивостью Релея-Тейлора. В этой работе обсуждаются вопросы улучшения диагностики эксперимента и обеспечения более высокого уровня взаимодействия между экспериментами и численным моделированием. Показано, что на рост зоны перемешивания оказывает влияние в большей степени пространственная структура начальных возмущений, чем их амплитуда. В данной работе обсуждаются различные постановки экспериментов.

В работе Cook, Cabot & Miller (2004) [14] разрабона крупно-вихревая методология моделирования неустойчивости Рэлея-Тейлора, выполнены расчеты РТН на сетке с числом ячеек 1152. Результаты подтверждают наличие автомодельной стадии турбулентности. Скорость J роста зоны перемешивания подчиняется известному закону вне зависимости от начальных условий. Значения коэффициента, а коррелируют с 1 результатами Glirnm и др. (2001) [17]. Эти значения ниже описанных в литературе, но соответствуют общей убывающей тенденции при увеличивающейся разрешающей способности сетки в численном моделировании (Dimonte и др. 2004) [18]. В работе предложено объяснение почему эксперименты с несмешивающимися жидкостями (Read 1984; Dimonte и Schneider 2000) [20,27], и численные модели (Glimm и др. 2001)[17] показывают более быстрые темпы роста зоны, чем случаи, учитывающие диффузию (Dimonte и др. 2004) [18].

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

Несомненно, предпосылки появления многопроцессорной вычислительной техники появились намного раньше её создания. Ярким примером параллельной обработки данных ещё в докомпьютерную эру, является расчет эволюции взрывной волны, произведенный Александром Андреевичем Самарским в начале 50-х годов. Несколько десятков человек производили расчеты на арифмометрах и передавали результаты друг другу на словах. Таким образом была реализована основная идея, положенная в основу параллельных вычислительных систем. Не случайно, так же, история создания параллельной вычислительной техники в работе [4] начинается с момента создания первых ЭВМ (конец 40 — начало 50-х годов XX столетия).

Так с появлением первых компьютеров основным вопросом стало повышение их производительности, которое достигалось за счет совмещения во времени различных этапов выполнения соседних команд. j Следующими этапами увеличения производительности стали реализация принципов конвеерности, векторности и независимости функциональных устройств. В дальнейшем, увеличение производительности было связано в первую очередь с принципом параллелизма, но это был внутренний параллелизм, предполагающий, что процессор содержит несколько функциональных устройств, работающих независимо друг от друга. Следующим шагом в увеличении производительности вычислительной машины является создание многопроцессорной конфигурации. Архитектура параллельных компьютеров развивается невероятно быстрыми темпами в самых различных направлениях [90]. Но, тем не менее, основополагающими являются 2 класса. Первый — машины с общей памятью (SMP), когда в системе присутствуют несколько равноправных процессоров, которые имеют одинаковый доступ к единой памяти. И компьютеры с распределенной памятью (МРР). Где по сути каждый узел является полноценным компьютером со своим процессором, памятью, системой ввода/вывода, операционной системой. Кроме того, существует огромное количество способов организации параллельных систем. В настоящее время ещё не существует полной классификации параллельных машин, отражающей все основные характеристики систем. Хотя попыток создания классификации было множество: в работе[4] рассмотрены наиболее известные из них. Наиболее ранней и общепринятой в настоящее время является классификация М. Флина. В её основу положено понятие потока, (т. е. последовательность команд и данных, обрабатываемых процессором). Выделяется четыре класса архйтектур.

SISD — (одиночный поток команд — одиночный поток данных) SIMD — (одиночный поток команд — множественный поток данных) MISD — (множественный поток команд — одиночный поток данных) MIMD — (множественный поток команд — множественный поток данных).

Класс MIMD чрезвычайно широк, поэтому многими авторами уточнялось понятие этого класса[85]. Понимание типа архитектуры вычислительной машины необходимо для написания эффективной программы, т.к. для каждого типа параллельных машин характерна своя модель программирования. Так для МРР и компьютеров с кластерной архитектурой стандартной моделью является программирование с использованием сообщений. Сообщение как средство межпроцессорного взаимодействия впервые было введено Хоаром в работе [93]. На основе сообщений было создано много систем параллельного программирования, но самым распространенным на сегодняшний день является MPI [84,95−97]. Интерфейс MPI поддерживает создание параллельных программ в стиле MIMD, что подразумевает объединение процессов с различными исходными текстами. Однако на практике гораздо чаще используется S PMD модель, в рамках которой для всех параллельных процессов используется один и тот же код.

В свою очередь применение в вычислениях многопроцессорных машин требует совершенствования численных алгоритмов и создания параллельных программ для многопроцессорных комплексов. Если для оценки алгоритма для ПК достаточно было трех характеристик — это точность алгоритма, объем требуемой памяти и количество операций, то для создания алгоритма для многопроцессорной системы помимо этих характеристик необходимо учитывать выбор эффективного способа распараллеливания. И, вопрос эффективности базового алгоритма является ключевым в выборе численного метода. Так теоретическим и практическим вопросам программирования применительно к параллельным системам посвящена работа[57], анализ создания эффективных параллельных алгоритмов проведен в работах[4, 87,88]. Наиболее известна работа, посвященная параллельным алгоритмам линейной алгебры[127]. Параллельные алгоритмы решения задач гидродинамики были рассмотрены в работе[99]. И в последнее время это направление активно развивается[102−104].

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

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

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

Нередко выполнение квадратичного закона роста ширины зоны и временное подобие осредненных профилей плотности в области перемешивания на некотором временном интервале считаются достаточными для того, чтобы рассматривать течение, как автомодельное. Для получения более обоснованного вывода об автомодельном характере турбулентного перемешивания на каком-то отрезке времени необходимо анализировать внутренние характеристики потока, базирующиеся на спектральном представлении газодинамических полей. В частности, наряду с установлением пространственно-временного подобия функционалов течения, необходимо проверять наличие инерционного интервала волновых чисел к, на котором реализуются закон Колмогорова-Обухова [33−35] распределения турбулентной кинетической энергии E (k)~k" 5/3. Для анализа результатов расчетов процесса турбулентного перемешивания необходимы эффективные методики исследования основныех характеристик возникающей турбулентности. Ширину зоны перемешивания возможно оценивать различными способами. Например, задав начальное положение границы в дальнейшем рассчитывать её эволюцию. Используют при этом как алгоритмы явного выделения границы [17], так и «лагранжевые» метки -«маркеры"[36]. Известно определение зоны по средним концентрациям [37]. В работе [ 14] вводится понятие зоны «смешанности», используя молярную массу вещества, по ней всесторонне исследуется перемешивающийся слой.

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

Основными целями диссертации являются:

• разработка и создание пакета прикладных программ для решения задач гидродинамики и исследования основных характеристик турбулентности.

• исследование гидродинамических неустойчивостей на многопроцессорных вычислительных системах.

• исследование тонкой структуры турбулентности.

Диссертация состоит из введения и трех глав, заключения и списка литературы.

Основные результаты диссертации:

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

2. Проведено прямое моделирование турбулентного перемешивания на примере двух задач: задачи о развитии неустойчивости РелеяТейлора и задачи о многократном прохождении ударной волны контактной границы элегаз/воздух. В расчетах было дано теоретическое объяснение эффекту уменьшения скорости роста зоны перемешивания при переходе к расчетным сеткам, содержащим большое число узлов. (порядкаЮ). Это объяснение позже получило экспериментальное подтверждение. Получена картина тонкой структуры турбулентности.

3. Для решения задач газодинамических задач с учетом теплопроводности на основе концепции адаптивных локально — измельчающихся сеток (в нашем случае на основе четырехугольных ячеек) были разработаны алгоритмы решения и создан пакет прикладных программ ADAPT 2D.