Содержание к диссертации
Введение
1 Современное состояние исследований в области турбулентного перемешивания 19
1.1 Три типа неустойчивостей 19
1.2 Экспериментальная информация 29
1.3 Обзор методик расчёта турбу лентного перемешивания 38
1.4 Проблемы численного моделирования 45
Выводы по главе 1 53
2 Методика расчёта турбулентного перемешивания в одномерных течениях (методика вихрь) 55
2.1 Уравнения модели 55
2.2 Расщепление по физическим процессам 62
2.3 Уравнения газовой динамики 64
2.4 Расчёт термодинамических функций 76
2.5 Расчёт лучистой теплопроводности 80
2.6 Разностные схемы для уравнений ту рбулентного перемешивания 90
2.7 Пакет программ ВИХРЬ
2.7.1 Основы организации 94
2.7.2 Расчёт масс намешанных веществ 95
2.7.3 Определение ширины зоны турбулентного перемешивания 97
2.7.4 Сервисные возможности 98
2.8 Примеры расчётов газодинамических течений 102
2.8.1 Сферическая задача (задача Ноха) 102
2.8.2 Столкновение ударных волн разной интенсивности 105
2.8.3 Расчёт распространения N-волны 109
2.9 Технология проведения расчётов турбулентного
перемешивания в одномерных течениях 112
2.9.1 Моделирование развития неустойчивости на контактной границе 115
2.9.2 Инициализация турбулентного перемешивания 121
2.9.3 О сходимости разностного решения уравнений турбулентного перемешивания 1 2.9.3.1 Безударные течения 124
2.9.3.2 Течения с ударными волнами 129
2.9.4 Определение центра конечно-разностной ударной волны 136
2.9.5 Тактика проведения расчётов 141
2.10 Примеры расчётов турбулентного перемешивания в модельных опытах 143
2.10.1 Опыт Е. Е. Мешкова в цилиндрической геометрии 145
2.10.2 Смыкание зон турбулентного перемешивания (опыт Е. Е. Мешкова в плоской геометрии) 148
2.10.3 Опыты Димонте 150
2.10.4 Опыт Барра 164
2.10.5 Опыт Poggi et al 171
Выводы по главе 2 180
3 Разработка и обоснование методики расчёта турбулентного перемешивания в двумерных течениях 181
3.1 Исходные уравнения 182
3.2 Отработка двумерной методики
3.2.1 Аналитические решения уравнения диффузии 188
3.2.2 Взаимодействие ударной волны с наклонной контактной границей 194
3.2.3 Сравнение с расчётами по другим методикам 211
3.2.4 Плоское двумерное течение 225
3.2.5 Сравнение с экспериментальными данными 245
3.3 Важные направления дальнейших исследований 266
3.3.1 Экспериментальные исследования 266
3.3.2 Численное моделирование 274
Выводы по главе 3 284
Заключение 285
Список сокращений и условных обозначений 288
Список литературы
- Проблемы численного моделирования
- Расчёт термодинамических функций
- О сходимости разностного решения уравнений турбулентного перемешивания 1 2.9.3.1 Безударные течения
- Экспериментальные исследования
Введение к работе
Актуальность работы. В исследованиях по физике высоких
плотностей энергии (например, в задачах по управляемому термоядерному
синтезу и ряде других научно-технических задач) существует проблема
устойчивого сжатия слоистых систем: в процессе сжатия перемешивание
слоёв препятствует достижению максимального энерговыделения
разрабатываемого устройства. В исследованиях подобных задач всё большую роль играют методы математического моделирования. Чрезвычайная сложность математического моделирования движения слоистых систем под действием ударных волн (УВ) обусловлена нелинейностью решаемых систем уравнений и значительной деформацией слоёв в процессе сжатия, что приводит к сильному искажению расчётной сетки и, как следствие, к потере точности вычислений. При моделировании рассматриваемого класса задач имеется ряд фундаментальных проблем, которые на сегодняшний день не имеют решения. Понимание закономерностей развития неустойчивости на контактной границе (КГ) и условий перехода к турбулентному перемешиванию (ТП) веществ является также одной из сложнейших проблем фундаментальной физики. Без решения существующих проблем результаты вычислительного эксперимента не позволяют спрогнозировать работу изучаемых систем.
Степень разработанности темы. Разработке физических и
математических моделей для описания турбулентного перемешивания,
возникающего на контактных границах при сжатии слоистых систем под
действием ударных волн, посвящено большое число работ. Значительный
вклад в развитие этого направления исследований внесли Андронов В.А.,
Анучина Н.Н., Бельков С.А., Демьянов А.Ю., Долголева Г.В., Зайцев С.Г.,
Змитренко Н.В., Иногамов Н.А., Кучеренко Ю.А., Лебо И.Г., Мешков Е.Е.,
Мхитарьян Л.С., Невмержицкий Н.В., Неуважаев В.Е., Никифоров В.В.,
Раевский В.А., Розанов Б.В., Стаценко В.П., Тишкин В.Ф., Янилкин Ю.В. и
многие другие специалисты как у нас в стране, так и за рубежом. Тем не
менее, при численном моделировании подобных задач с использованием
уравнений Эйлера (расчёт средних газодинамических величин) и
полуэмпирических моделей турбулентности (расчёт турбулентных
характеристик) существует ряд нерешённых проблем.
Первая проблема связана с моделированием перехода от ламинарного
течения к турбулентному. Для расчёта скорости роста начальных
возмущений КГ широко используются прямое численное моделирование и
эмпирические модели. Для малоамплитудных одномодовых начальных
возмущений КГ оба подхода с удовлетворительной точностью воспроизводят
имеющуюся экспериментальную информацию по скорости роста
возмущений (роста начальной амплитуды) на линейной и нелинейной стадиях развития возмущений. Стадия перехода от неустойчивости к ТП с использованием прямого численного моделирования не моделируется, а при
моделировании ТП с использованием полуэмпирических моделей
турбулентности возникает проблема определения момента перехода от этапа развития неустойчивости к ТП и задание начальных данных для решения уравнений ТП.
Вторая фундаментальная проблема связана с отсутствием сходимости разностного решения уравнений Эйлера на фронте конечно-разностной ударной волны. В модели Никифорова ширина зоны ТП в численном расчёте уменьшается при измельчении размера счётной ячейки разностной сетки, что является следствием неограниченного роста градиентов газодинамических величин на фронте скачка.
Третья фундаментальная проблема возникает при моделировании взаимодействия ударной волны с турбулентным потоком. Установлено, что при прохождении УВ по турбулентному потоку однородного вещества или смеси разноплотных веществ внутренняя структура скачка меняется. Если пульсации газодинамических величин перед фронтом превышают некоторый порог, фронт скачка деформируется под действием пульсаций и вихрей и становится неустойчивым.
Таким образом, сложность решаемых практических задач и нерешённые математические проблемы, присущие современным численным методикам, выдвигают на первый план задачи повышения точности вычислений и обоснования результатов математического моделирования. Полномасштабную проверку адекватности используемых физических и математических моделей и подходов, используемых при описании развития неустойчивости и ТП на контактных границах многослойных систем, можно осуществить при сопоставлении результатов вычислительного эксперимента с экспериментальными данными.
На сегодняшний день в литературе наиболее широко представлены
результаты экспериментальных исследований эволюции зоны ТП,
возникающей на КГ двухслойных газовых систем при доминирующей роли
неустойчивости Рихтмайера-Мешкова. Для методик, ориентированных на
расчёты развития ТП в многослойных системах, имеющейся
экспериментальной информации недостаточно. Для детального тестирования двумерных численных методик нужны эксперименты с трёхслойными системами. В опытах с трёхслойными системами (по сравнению с двухслойными системами) ТП развивается под действием ударных волн и (или) волн разрежения, циркулирующих в центральном слое, что обычно характерно для слоистых систем, используемых в устройствах по получению высоких концентраций энергии. Следовательно, информация, полученная в опытах с трёхслойными системами, наиболее точно отражает характер течения, формирующегося в многослойных системах под влиянием УВ.
Целью работы является создание новых математических моделей, алгоритмов и комплексов программ для численного моделирования турбулентного перемешивания в устройствах по получению высоких плотностей энергии, которые в отличие от имеющихся аналогов, адекватно
учитывают физику перемешивания и с высокой точностью описывают экспериментальные данные. Для достижения этой цели решались задачи, которые формулируются следующим образом:
-
анализ современного состояния экспериментальных, теоретических и численных исследований в области турбулентного перемешивания;
-
разработка трёхпараметрической разностной схемы для решения уравнений газовой динамики; определение оптимального (с точки зрения точности и монотонности численного решения) набора параметров разностной схемы;
-
разработка алгоритмов построения сходящегося решения разностных уравнений как для безударных течений, так и течений с ударными волнами;
-
повышение точности расчёта генерационных членов уравнений модели Никифорова на фронте ударной волны, движущейся по турбулентному полю; разработка алгоритма для описания взаимодействия ударной волны с турбулентным потоком;
-
повышение точности вычисления генерационных членов уравнений турбулентного перемешивания модели Никифорова за счёт уменьшения интенсивности энтропийных следов и паразитических осцилляций численного решения;
-
построение алгоритма проведения вычислительного эксперимента, включающего моделирование этапа развития неустойчивости (рост амплитуды возмущений на контактной границе), выбор критерия перехода к перемешиванию и построение начальных данных для решения уравнений модели Никифорова;
-
построение аналитических решений ряда модельных уравнений с целью обоснования и тестирования математической модели; математическая формулировка постановок опытов с целью получения новой информации по развитию турбулентного перемешивания на контактных границах трёхслойных газовых систем.
Научная новизна результатов.
-
разработана трёхпараметрическая разностная схема для решения уравнений газовой динамики; на основе анализа диссипативных и дисперсионных погрешностей разностной схемы предложен оптимальный набор параметров;
-
впервые предложены методы построения сходящегося решения разностных уравнений как для безударных течений, так и течений с УВ, движущимися по зоне ТП;
-
впервые разработаны математические модели расчёта взаимодействия УВ с турбулентным потоком разноплотных газов, учитывающие деформацию фронта скачка;
-
сформулирован ряд новых тестовых задач для оценки точности вычислений, в том числе, с использованием аналитических решений модельных уравнений;
5. на основе аналитического решения задачи о взаимодействии УВ с наклонной КГ и численных экспериментов предложены математические постановки новых опытов по изучению развития ТП на контактных границах трёхслойных газовых систем (одна из контактных границ располагалась параллельно фронту УВ, вторая – под различными углами).
Теоретическая значимость работы.
Новые результаты представленных в работе исследований и разработанные на их основе подходы для решения фундаментальных проблем численного моделирования турбулентного перемешивания вносят значительный вклад в повышение точности математического моделирования процессов, происходящих в устройствах по получению высоких плотностей энергии, и способствуют лучшему пониманию закономерностей развития неустойчивостей и перемешивания на КГ под действием ударных волн.
Практическое значение работ.
Разработанные в рамках методики ВИХРЬ [1а, 2а] программы моделирования ТП внедрены в методику СНД-ТУР [3а, 4а] (расчёты ТП в одномерном приближении). Для методики РАМЗЕС (КОРОНА) разработаны программы учёта ТП по модели Никифорова (расчёты в двумерном приближении). Предложенные тестовые задачи, в том числе, на основе найденных аналитических решений ряда модельных уравнений и экспериментальная информация используются для отработки технологии счёта практических задач и повышения точности двумерных численных методик (ЛЭГАК-ВКЛ, РАМЗЕС, МИД (МИМОЗА)); программные комплексы используются для моделирования процессов, происходящих в устройствах по получению высоких плотностей энергии.
Методы исследования.
Для решения поставленных задач используются математическое моделирование, численные методы и комплексы программ. Система исходных дифференциальных уравнений аппроксимируется разностными уравнениями, которые решаются итерационными методами. Для оценки точности вычислений построены аналитические решения ряда модельных уравнений. С целью проверки адекватности математической модели исследуемым объектам привлекается экспериментальная информация. Разработанные математические методы и алгоритмы реализованы в виде комплексов программ.
На защиту выносятся следующие результаты, полученные автором
лично или вклад автора в которые был определяющим:
1. математическая модель и комплекс программ ВИХРЬ, предназначенные
для расчёта одномерных течений в лагранжевых переменных,
учитывающие процессы газодинамики, детонации, лучистой
теплопроводности, турбулентного перемешивания по модели Никифорова;
-
модуль расчёта турбулентного перемешивания, обобщённый на случай неравновесной нестационарной многокомпонентной плазмы в рамках одномерной программы радиационной газовой динамики СНД-ТУР;
-
математическая модель расчёта турбулентного перемешивания, разработанная для двумерной методики РАМЗЕС, основанная на эйлерово-лагражевом подходе для моделирования нестационарных течений в криволинейной системе координат;
-
алгоритм проведения вычислительного эксперимента по исследованию ТП, включающий метод инициализации перемешивания, основанный на расчёте этапа развития неустойчивости (расчёт скорости роста пузырей и струй), предложенного критерия перехода к перемешиванию и разработанного алгоритма построения начальных данных для решения уравнений модели Никифорова;
-
алгоритмы построения сходящегося решения разностных уравнений как для безударных течений, так и течений с ударными волнами;
-
алгоритм для описания взаимодействия ударной волны с турбулентным потоком, учитывающий деформацию фронта скачка;
-
математическая формулировка постановок опытов с целью получения новой информации по развитию турбулентного перемешивания на контактных границах трёхслойных газовых систем с наклонной контактной границей и взаимодействия ударной волны с турбулентным потоком.
Достоверность и апробация результатов. Достоверность полученных результатов подтверждается результатами сравнения с тестовыми задачами, экспериментальной информацией, с данными, полученными по другим методикам. Результаты работ, направленных на создание, обоснование и применение методики, изложены в 113 отчётах ВНИИЭФ. Основные результаты работ по теме диссертации:
опубликованы в журналах, рекомендованных ВАК РФ [1а-17а], две работы [18а, 19а] входят в систему цитирования SCOPUS;
изложены в монографии “Моделирование неустойчивости и турбулентного перемешивания в слоистых системах”, Саров - 2010, 414c.
опубликованы в других периодических журналах ([1в-8в]);
изданы в виде отчёта в USA, LANL, LA - 12896, 1995, 206с. ([25d]);
изложены в 113 научно-исследовательских отчётах ВНИИЭФ;
изложены в трудах международных конференций и семинаров ([1d-28d]):
1) конференция по физике сжимаемого турбулентного перемешивания
(France, Marseille - 1997; USA, Pasadena - 2001; UK, Cambridge - 2004;
France, Paris - 2006; Россия, Москва - 2010);
-
Международные Забабахинские научные чтения (Снежинск - 2001);
-
Международные Харитоновские тематические научные чтения (Саров -2007, Саров - 2009, Саров - 2011).
Соответствие паспорту специальности. Диссертационная работа
соответствует формуле научной специальности 05.13.18 – “Математическое
моделирование, численные методы и комплексы программ” (физико-
математические науки) в пунктах 1, 3, 4, 5, 6, 7.
Проблемы численного моделирования
В теоретических и численных исследованиях развития неустойчивости форма начальной поверхности (начального возмущения) представляется в виде конечного ряда Фурье в котором am, k = 2л/Х, m, фm, X - соответственно амплитуда возмущения, волновое число, номер моды (гармоники), фазовый угол и длина волны. Строго говоря, в практических приложениях форма возмущённой поверхности нам не известна, в некоторых ситуациях можно лишь оценить начальную шероховатость -LQ. С целью получения качественных оценок по развитию одномодовых начальных возмущений в теоретических исследованиях форму начального возмущения достаточно задать в виде синусоиды. Для анализа конкуренции (взаимодействия) различных мод начальное возмущение необходимо представлять в виде многомодового ряда.
На рис. 1.2 приведены две из используемых в теоретических и численных исследованиях формы первоначально возмущённой поверхности. Пунктирная кривая задана уравнением для одномодового возмущения y = a1cos(k1 x), a=0,05, kх=2, где х и у - пространственные координаты (у направлена вертикально вверх, х слева направо), величина ал _ называется амплитудой начального возмущения; кх = 2п/\1 есть волновое число; X і - длина волны возмущения. Сплошная кривая соответствует начальному возмущению КГ, заданному в виде суммы двух мод, j = 21cos( 1x) + 22cos( 2x); ufj=0,05, kx=2; а2=0,02, k2=\0.
При заданных условиях длины волн гармоник равны =71 и Л,2= 0,27Г. Существуют и другие способы задания возмущений на поверхности раздела жидкостей. Например, в численных расчётах в окрестности КГ с помощью датчика случайных чисел можно сформировать пульсации газодинамических величин.
Неустойчивость Рэлея-Тейлора. Рост начальной шероховатости границы (амплитуды возмущённой поверхности), разделяющей слой “тяжёлой” жидкости, поддерживаемый “лёгкой” жидкостью, известен как РТ неустойчивость. Такая неустойчивость появляется под влиянием земной гравитации, если система (тяжёлая жидкость сверху, лёгкая снизу) находится в покое, или при ускорении всей жидкой системы в направлении от лёгкой жидкости к тяжёлой (градиенты давления и плотности направлены в разные стороны). Неустойчивость РТ развивается в окрестности шероховатой (неровной) КГ разноплотных жидкостей. В процессе развития неустойчивости амплитуда шероховатости начинает расти: тяжёлая жидкость проваливается в лёгкую в виде струй, лёгкая жидкость поднимается в тяжёлую в виде пузырей. Пузырями принято называть порции лёгкой жидкости, проникающие в тяжёлую, струями называются порции тяжёлой жидкости, проникающие в лёгкую жидкость. Когда ускорение направлено от тяжёлой жидкости к лёгкой, граница будет вести себя устойчиво (начальная шероховатость КГ затухает со временем). В процессе развития неустойчивости амплитуда начальной шероховатости увеличивается, затем КГ разрушается и между тяжёлой и лёгкой жидкостями образуется слой (зона) перемешанных веществ.
Отметим некоторые факторы, влияющие на развитие неустойчивости КГ: отношение плотностей на КГ, поверхностное натяжение, вязкость и сжимаемость жидкостей, трёхмерные эффекты, вид ускорения (постоянное, возрастающее и т. д.), степень гетерогенности смеси (см. ниже), геометрия явления (плоская, цилиндрическая, сферическая). В практических приложениях существует ещё ряд других факторов, влияющих на развитие неустойчивости, например, вид уравнения состояния (УРС) контактирующих жидкостей, процессы теплопроводности, изменение фазового состояния контактирующих жидкостей и т. д.
В случае многослойных систем устойчивость движения необходимо исследовать на всех КГ системы. При этом разные КГ системы могут двигаться с различными ускорениями, зависящими от времени, а слои с течением времени могут сжиматься или расширяться. Следует также иметь в виду, что отдельные зоны перемешивания, возникшие на неустойчивых КГ, могут со временем смыкаться между собой, в результате чего число перемешанных веществ в зоне может быть больше двух.
Эксперименты по исследованию неустойчивости проводятся, в основном, на КГ жидкость-жидкость. Для предотвращения размывания КГ до начала опыта (вследствие взаимодиффузии контактирующих смешиваемых жидкостей) в некоторых экспериментах используются несмешиваемые жидкости. Развитие малоамплитудных возмущений на КГ принято делить на четыре стадии. Следует отметить, что такое деление является достаточно условным. Стадия 1 - линейная (рис. 1.3). Если амплитуда начального периодического одномодового возмущения значительно меньше длины волны, то амплитуда возмущения растёт экспоненциально по времени. На этой стадии форма начального возмущения остается отклонение линейной теории развития возмущений от экспериментальных данных. симметричной относительно среднего положения (у = 0). Линейная стадия продолжается до тех пор, пока амплитуда связана с длиной волны приближённым неравенством а 0,1 -0,4 X. В дальнейшем наблюдается
Эта стадия роста возмущений может быть теоретически проанализирована с помощью линеаризованных уравнений движения. Для расчёта роста возмущений также широко используется прямое численное моделирование на основе уравнений Эйлера или Навье - Стокса (без использования моделей турбулентности). Расчёты с одномодовым начальным возмущением, по-видимому, следует считать познавательными, поскольку в природе не существует идеально возмущённой поверхности (если специально не создавать такое возмущение в экспериментах). Расчёты, в которых начальное возмущение содержит несколько мод, ближе к реальности и дают больше информации исследователю. В частности, такие расчёты позволяют проследить за процессом взаимодействия (конкуренции) пузырей, имеющих различные радиусы (предполагаем, что пузыри имеют сферическую форму).
Линейная стадия развития возмущений, при t = 0 задано одномодовое периодическое возмущение Стадия 2 – нелинейная (рис. 1.4). На этой стадии начальные амплитуды возмущений растут нелинейно до размеров h1 ,h2 (h1 – глубина проникновения лёгкой жидкости в тяжёлую; h2 – глубина проникновения тяжёлой жидкости в лёгкую). Развиваются гармоники с различной длиной и амплитудой волны. Формируются пузыри и струи различных размеров.
Расчёт термодинамических функций
Модификация модели Никифорова. В последние годы совершенствованием модели Никифорова занимается В. И. Козлов. Наиболее сложные с точки зрения практической реализации изменения модели Никифорова, предложенные В. И. Козловым [74], связаны с решением уравнения для скорости диссипации Q кинетической энергии турбулентности: на фронте УВ в модифицированной модели интегрируется уравнение для скорости диссипации, отличное от уравнения Никифорова. Модифицированное уравнение для скорости диссипации, записанное в лагранжевых массовых переменных, имеет вид dQ д
В случае принадлежности счётной ячейки фронту конечно-разностной УВ функция у/ полагается равной 0 (в противном случае у/=\). Функция у/ определяется соотношением 1, x x О, X X e3/2divU sw где X = W = XSW – задаваемый параметр. SW QCS При вычислении “адиабатического градиента” плотности в модифицированной модели (как и при решении уравнения для скорости диссипации) учитывается принадлежность ячейки разностной схемы фронту УВ: As = y/ дР дг р\ + PCv Л; dp р дг В модифицированной модели изменены формулы расчёта коэффициентов диффузии D и Д: вместо двух коэффициентов модели Никифорова в модифицированной модели рассчитывается один коэффициент - Д, а второй полагается равным D = 2 Dt. Коэффициент турбулентной диффузии в модифицированной модели имеет вид:
Важное преимущество данного способа вычисления f1 перед определением f1 в модели Никифорова состоит в монотонности получаемого профиля f1. Монотонность профиля f1 обеспечивается монотонностью коэффициента турбулентной диффузии Dt (в модели Никифорова в профилях коэффициентов диффузии содержались паразитические осцилляции, так как рассчитывались с использованием градиентов газодинамических величин).
В уравнениях диффузии для Ск (концентраций компонент смеси), Q, е2 (поперечной кинетической энергии турбулентности), Е (удельной внутренней энергии) уточнены константы при Dt.
Обобщено выражение для функции у, характеризующей скорость обмена энергией между продольными и поперечными пульсациями скорости: На фронте УВ используется ограничение на скорость обмена энергией между продольной и поперечной компонентами турбулентных пульсаций скорости где Mt - турбулентное число Маха, Сs - скорость звука. Рекомендованные В. И. Козловым значения констант: 2 = 5; у2 = 0,25; Уз = 0,5. Параметр, характеризующий скорость турбулентной диффузии, как и для К- модели принимается равным а = 0,09.
При численном моделировании различного рода задач широкое распространение получил метод расщепления по физическим процессам. Такой подход к решению систем уравнений, моделирующих несколько физических процессов, имеет ряд достоинств. Во-первых, он позволяет легко варьировать набором физических процессов, которые в силу специфики задачи следует учитывать при решении той или иной задачи и, во-вторых, применять для расчёта конкретного процесса хорошо разработанные разностные методы. Так, например, система уравнений (2.1) используется для расчёта следующих физических процессов: перенос количества движения и тепла, их диффузию и турбулентное перемешивание. Метод расщепления по физическим процессам использовался и при разработке методики ВИХРЬ, предназначенной для расчёта одномерных течений с учётом лучистой теплопроводности и турбулентного перемешивания.
В пакете программ ВИХРЬ система уравнений (2.1) решается методом расщепления по физическим процессам. На первом этапе интегрируются уравнения газовой динамики, записанные в лагранжевых массовых координатах с учётом турбулентного давления Pt : dU д(Р +Pt+g ) дг =-S , — = U, pV=m, dt dm dt дЕ , ч dSU дР , Ч2 dSU ,„„, dt dm dt dm В приведённых уравнениях m - массовая лагранжева переменная (масса контрольного объёма), которая связана с плотностью соотношением m = jp(r,t)Sdr; S - площадь поверхности, ограничивающей контрольный объём V; с - скорость звука; g - математическая вязкость. На втором этапе решаются уравнения для турбулентных величин eh е2, F, Q, R, W, В без учёта диффузионных членов. На этом же этапе пересчитываются скорость и удельная внутренняя энергия за счёт турбулентных слагаемых:
Важное требование, предъявляемое к разностным схемам, аппроксимирующим систему уравнений газовой динамики (2.2), - монотонность профилей газодинамических величин. Это требование не менее важное, чем требование к точности аппроксимации, так как в уравнения для турбулентных величин входят градиенты газодинамических величин, а наличие паразитических осцилляций в решении (например, в окрестности фронта ударной волны) может приводить к физически необоснованному развитию ТП. При разработке пакета программ ВИХРЬ было опробовано несколько явных (и неявных) разностных схем для интегрирования уравнений газовой динамики. В результате в пакете программ ВИХРЬ было реализовано две разностные схемы: явная схема с улучшенными диссипативными свойствами и неявная, которая используются в методике ВИХРЬ для расчёта газодинамических течений с учётом ТП. Рассмотрим построение трёхпараметрической неявной разностной схемы.
Интегрирование системы уравнений газовой динамики выполняется по трёхпараметрической разностной схеме, являющейся модификацией схемы [65]. Ниже используются следующие обозначения: т - шаг интегрирования по времени, п - номер временного шага. Целые нижние индексы у величин указывают на принадлежность величины к узлам разностной сетки (боковым границам ячеек), полуцелые - центрам ячеек. В методике ВИХРЬ пространственные координаты и массовые скорости относятся к узлам, прочие газодинамические величины центрам ячеек. Система обозначений демонстрируется на рис. 2.1.
О сходимости разностного решения уравнений турбулентного перемешивания 1 2.9.3.1 Безударные течения
Наличие в разностной схеме трёх параметров (Зі, (Зг, (Зз позволяет управлять до некоторой степени не только погрешностью аппроксимации, но и диссипативными и дисперсионными свойствами схемы.
Для полной постановки разностной задачи в области искомого решения задаются начальные условия, а на левой и правой границах задачи - граничные условия. В методике ВИХРЬ реализованы следующие граничные условия 1. U = const - задаётся постоянное значение скорости, 2. U = U(t) - величина скорости меняется во времени, 3. Р = const - задаётся постоянное значение давления, 4. Р = P(t) - величина давления меняется во времени, „ сгТ4 5. Р = . Схема с улучшенными диссипативными свойствами. Для расчёта распространения УВ в атмосфере на большие расстояния Ю.А. Бондаренко разработал схему с улучшенными диссипативными свойствами, которая реализована в методике ВИХРЬ. Приведем расчётные формулы схемы.
В работе Шульца [133] для более точного определения схождения ударных волн к центру в сферическом одномерном случае или к оси симметрии в одномерном цилиндрическом случае предложено вязкость в уравнении для скорости учитывать неконсервативным способом. По этой причине закон сохранения внутренней энергии в методе Шульца записывается в иной форме по сравнению с (2.2). Таким образом, система дифференциальных уравнений в методе Шульца отличается от системы (2.2) уравнениями для скорости и внутренней энергии и приобретает следующий вид: du Эр d(S-g) dr — =-S— , — = u, dt 3m 3m dt (2.12) dE d(S-u) n du — =-p -g-S , pv= const. dt 5m 5m
В результате такой записи уравнения энергии вклад вязкости в повышение энтропии около оси или центра симметрии получается меньше, чем при обычной записи уравнения для скорости и уравнения для внутренней энергии, а в плоском случае никаких отличий не имеется. К тому же, на численных примерах установлено, что система (2.12) позволяет более точно рассчитывать значение давления. Учитывая сказанное, для расчёта одномерных течений газа будем использовать запись дифференциальных уравнений в форме (2.12). К этим уравнениям присоединяется уравнение состояния, которое будем считать заданным в виде: р = р (р, Е) .
При построении разностной схемы как и в первой схеме, величины S, г, и будем относить к узлам сетки, а величины р, Р, Е, с, V - к серединам интервалов. Узлы сетки обозначим целыми индексами, центрам ячеек припишем полуцелые индексы.
Для интегрирования системы уравнений (2.12) используется явная разностная схема типа “крест“. Счётные формулы разностной схемы приведём в том порядке, в котором они используются на одном шаге по времени в программе. Вначале приведем разностные уравнения для координат и плотности, затем формулы для счётных вязкостей, за ними следуют разностные уравнения для энергии и для скорости.
Расчёт начинается с вычисления допустимого счётного шага по времени, после чего определяются новые координаты узлов разностной сетки по формуле п+1 п+1/2 П+1/2 п+1 n+l/2 л. n+l/2 n+l ,n = u , где At = t -1 . n+l/2 j At г;+1-г; После вычисления новых координат с их помощью вычисляются из закона сохранения массы новые плотности по формулам р" = ї-
Разностная схема “крест” является бездиссипативной, зато обладает дисперсией. В связи с этим на разрывных решениях схема не монотонна. Для подавления осцилляций в разностные уравнения добавляется искусственная вязкость. В данной работе вязкость вычисляется через градиент скорости и равна сумме квадратичной и линейной вязкостей. При сжатии работает квадратичная вязкость (безразмерный коэффициент вязкости задаётся пользователем в начальных данных). Для монотонизации профилей на ударных волнах и на передних фронтах волн разрежения имеется возможность к квадратичной вязкости добавлять линейную вязкость с разными безразмерными коэффициентами при сжатии ячейки и при её разрежении.
Вычисление значения вязкости начинается с определения градиента скорости по формуле п+1 п+1 п+1 п+1/2 п+1/2 о П+1 о П+1 п+1/ _ j+l j uj+l uj ди yj+l/2 2 mJ+l/2 pdr
В этой формуле S есть площадь поверхности, ограничивающей ячейку сетки в направлении координаты г.
Окончательное значение искусственной вязкости в ячейке сетки вычисляется по формулам, которые записываются в виде Здесь C:+1/2( Ры/г Е" 2 ] изоэнтропическая скорость звука, рассчитываемая по уравнению состояния в данной ячейке сетки, Q и L - безразмерные коэффициенты квадратичной и линейной вязкостей соответственно. Отметим, что величина %, определённая в ячейках сетки и вычисляемая формулой (2.13), используется при вычислении потоков энергии в искусственной теплопроводности. Эта величина имеет смысл размерного коэффициента искусственной теплопроводности, который определён в центрах ячеек сетки.
Безразмерные коэффициенты квадратичной и линейной вязкостей в данной ячейке в данный момент времени определяются формулами При таком выборе коэффициента Q квадратичная вязкость работает только при сжатии, так как задаваемый коэффициент квадратичной вязкости отличен от нуля только в окрестности фронта ударной волны. В то же время линейная вязкость может работать как при сжатии, так и при разрежении, если AV и АР заданы отличными от нуля. Безразмерный коэффициент квадратичной вязкости Q задаётся пользователем в начальных данных. Безразмерные коэффициенты линейной вязкости для сжатия АР и для разрежения AV также определены в начальных данных.
Экспериментальные исследования
В области 1 (0 г г і) находится газ, плотность которого pi = 0,1663. Область 2 (гі г г2) заполнена газом с плотностью р2 = 1,205. На КГ число Атвуда составляет А = 0,757. В расчётах используется уравнение состояния, имитирующее несжимаемую жидкость. Это уравнение имеет следующий вид:
В уравнении состояния через с\ и с2 обозначены массовые концентрации соответствующих газов. Представленные ниже результаты моделирования задачи получены при следующих значениях параметров: г і = 20, г2 = 60, g = 1 - постоянное ускорение.
Расчёты на равномерной сетке. Инициализация ТП выполнялась по модели Никифорова [5] с R0 = 0,5 - начальное значение пульсаций плотности, которое 125 задавалось в одной ячейке слева и в одной ячейке справа от КГ (задание профиля пульсаций плотности в форме прямоугольника).
На рис. 2.36 приведены r диаграммы границ зоны перемешивания, полученные в расчётах на равномерных сетках с различным числом ячеек. Границы зоны ТП определялись по уровню массовых концентраций C = 0,99 и C = 0,01. На рис. 2.36 приняты следующие обозначения: M1 – число ячеек в первой (левой) области, М2 – число ячеек во второй (правой) области. Отметим, что в приведённых расчётах изменяемыми параметрами являются размер счётной ячейки и, как следствие, начальная ширина зоны ТП L0 = h1 + h2 , где h1 – размер счётной ячейки слева от контактной границы, h2 – размер счётной ячейки справа от контактной границы. Для самой грубой сетки (M1 = 100, M2 = 440) линейные размеры ячеек, прилегающих слева и справа к КГ соответственно равны: h1= 0,2, h2= 0,09. Неравенство h1 h2 указывает на то, что при задании ширины начальной зоны учитывалось наличие асимметрии зоны на момент включение уравнений ТП в расчёт: струи тяжёлой жидкости проникают в лёгкую жидкость на большую глубину, чем пузыри лёгкой жидкости в тяжёлую. Из рис. 2.36 следует, что при измельчении сетки ширина зоны ТП заметно уменьшается.
На рис. 2.37 показана зависимость L(S) (ширины зоны ТП) от пройденного пути S = gt2 /2. Данные рис. 2.37 демонстрируют заметную зависимость скорости роста зоны на начальном этапе развития неустойчивости от ширины начальной зоны ТП L0 и уменьшение этой зависимости с течением времени (при выходе решения на автомодельный режим). В выполненных расчётах начальные значения пульсаций плотности задавались в одной ячейке слева и в одной ячейке справа от контактной границы. При таком способе инициализации ТП (метод Никифорова) измельчение сетки приводит к уменьшению начальной ширины зоны ТП L0 = h1 + h2 и увеличению градиента плотности в окрестности КГ. В пределе при L0 = 0, что соответствует заданию нулевых значений пульсаций плотности, уравнения ТП модели Никифорова имеют тривиальное решение.
Чтобы определить какую роль в представленных выше расчётах играет размер счётной ячейки и начальная ширина зоны ТП L0, были выполнены расчёты, в которых менялся алгоритм задания начальных значений пульсаций плотности.
Как и прежде, на самой грубой сетке значения пульсаций плотности задавались только в одной ячейке слева и в одной ячейке справа от контактной границы. При каждом следующем измельчении сетки в 2 раза увеличивали в 2 раза число ячеек слева и справа от контактной границы, в которых задавались значения пульсаций плотности. Тем самым фиксировалась начальная ширина зоны ТП L0, однако в начальной зоне число счётных ячеек не сохранялось (следовательно, градиенты газодинамических величин на КГ зависели от размера ячейки).
Расчёты на неравномерной сетке. Изучим поведение разностного решения уравнений ТП при сохранении как начальной ширины зоны ТП, так и числа счётных ячеек в этой зоне с целью сохранения градиента плотности в окрестности КГ. Для сохранения интегрального значения начальных пульсаций плотности (начальной ширины зоны ТП) будем использовать следующий приём. Размеры счётной ячейки слева и справа от контактной границы сохраняются постоянными во всех расчётах. На других участках областей 1 и 2 расстановка ячеек меняется от расчёта к расчёту. Данный способ задания начальных значений пульсаций плотности отличается от рассмотренного выше приёма тем, что в этом случае одновременно сохраняются и начальная ширина зоны ТП, и число счётных ячеек в начальной зоне (как следствие, сохранение градиентов средних газодинамических величин на КГ).
Первый расчёт проведён на самой подробной равномерной сетке (как и в предыдущем подразделе). Во втором расчёте в областях 1 и 2 сетка сгущена в окрестности контактной границы по закону геометрической прогрессии. В отличие от второго расчёта в третьем расчёте увеличен знаменатель геометрической прогрессии при расстановке сетки (q1 – знаменатель геометрической прогрессии для области 1, q2 – для области 2). При этом во всех расчётах размер ячейки слева и размер ячейки справа от контактной границы сохранены.
Из представленных результатов следует, что при проведении расчётов на различных сетках требуется сохранять ширину начальной зоны, в которой задаются начальные данные для решения уравнений ТП. А чтобы решение задачи не зависело от размера счётной ячейки следует выполнить расчёты на сходимость решения (в том числе и по градиентам средних газодинамических величин) на последовательности измельчающихся сеток.