РОССИЙСКАЯ АКАДЕМИЯ НАУК

ИНСТИТУТ ПРОБЛЕМ БЕЗОПАСНОГО РАЗВИТИЯ

АТОМНОЙ ЭНЕРГЕТИКИ

(ИБРАЭ РАН)

На правах рукописи

Кондаков Василий Гаврильевич

Обобщение схемы КАБАРЕ на многомерные

уравнения задач газовой динамики

специальность: 05.13.18 – «Математическое моделирование, численные методы и комплексы программ»

Диссертация на соискание ученой степени кандидата физико-математических наук

Научный руководитель:

д.ф.-м.н. С.А. Карабасов Москва - Содержание Введение

Глава 1. Схема КАБАРЕ для уравнений одномерной газовой динамики....... 11  ОПИСАНИЕ СХЕМЫ КАБАРЕ

Сеточные множества и сеточные функции

Аппроксимация системы дивергентных уравнений

Вычисление новых значений потоковых переменных

Звуковые точки

Вариант № 1. Ускорение потока влево

Вариант № 2. Торможение потока вправо

Вариант № 3. Торможение потока влево

Вариант № 4. Ускорение потока вправо

Учет особенностей уравнения состояния

Граничные условия

СВОЙСТВА СХЕМЫ КАБАРЕ

Аппроксимация

Другая форма представления алгоритма

Свойства акустического приближения

Некоторые вычислительные свойства нового алгоритма

ПРИМЕРЫ РАСЧЕТОВ

Задача Сода

Задача Шу-Ошера

Глава 2. Обобщение схемы Кабаре на двумерные ортогональные сетки....... 36  Уравнения Эйлера и инварианты Римана

Схема КАБАРЕ

Звуковые точки

Вариант №1. Ускорение потока вправо

Вариант №2. Ускорение потока влево

Вариант №3. Торможение потока, двигающегося вправо

Вариант №4. Торможение потока, двигающегося влево

Граничные условия

Выбор шага по времени

Постановка тестовых примеров

Задача №1. Изолированный Q-вихрь.

Задача №2. Двойное маховское отражение

Глава 3. Схема КАБАРЕ в трехмерной области на регулярной сетке............. 69  Разностная схема

Вычисление промежуточных значений консервативных величин........... 73  Вычисление потоковых переменных вдоль оси x

Вычисление потоковых переменных вдоль оси y

Вычисление потоковых переменных вдоль оси z

Вычисление новых значений консервативных величин.

Аппроксимация тензора вязкости

Выбор шага по времени

Особенности распараллеливания

Обтекание обратного уступа

Постановка задачи и примеры расчетов в литературе

Постановка задачи

Результаты расчета

Заключение к Главе 3

Глава 4. Схема КАБАРЕ в трехмерной области на произвольной сетке с гексагональными ячейками.

Разностная схема

Идеальный газ

Слабосжимаемый газ

Тензор вязких напряжений

Вычисление тензора вязких напряжений

Вычисление шага по времени

Глава 5. Моделирование высокоскоростной турбулентной струи из одноконтурного сопла, эксперимент JEAN.

Постановка задачи и примеры расчетов в литературе

Результаты расчетов по методу Кабаре

Заключение к Главе 5

Защищаемые положения

Литература

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

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

С одной стороны, метод должен обладать высоким порядком точности в тех областях, где решение является гладким. Даже в этих областях получение решения требуемой точности в прикладных задачах газовой динамики представляет известные трудности. В частности, известными недостатками классических конечно-разностных методов для задач газовой динамики в эйлеровых переменных являются большие диссипативные и дисперсионные ошибки [2-4]. Одним из часто используемых подходов для улучшения диссипативных и дисперсионных свойств классичеcких схем являются подходы, основанные на применении центральныx конечноразностныx схем повышенного порядка аппроксимации. В таких подходах за основу берутся уравнения переноса, записанные в неконсервативной форме, и проводится оптимизация коэффициентов разностного шаблона уравнения линейного переноса с целью минимизации дисперсионной ошибки при заданном порядке аппроксимации [5-7]. Оптимизированные конечноразностные схемы были разработаны как способ перенесения замечательных свойств спектральных методов в пространственно-временную область для снятия проблем, возникающих при постановке граничных условий в областях конечного размера. Оптимизированные конечно-разностные схемы эффективны для задач линейного переноса на слабо неоднородных сетках при наличии медленно меняющегося фонового поля на сетке: уже 4-го порядка схемы демонстрируют рекордную точность при переносе быстро осциллирующих решений по медленно меняющемуся среднему полю. К другому классу методов относятся так называемые методы высoкого разрешения, которые основаны на решении уравнений переноса в консервативной форме [8]. Для улучшения разностных свойств используются формулы высокого порядка аналитического восполнения (реконструкции) функции потоков или переменных. Для обеспечения баланса между дисперсионными и диссипативными ошибками, возникающими в решении при наличии плохо разрешённых градиентов, используются подходы на основе введения дополнительных диффузионных или антидиффузионных членов, улучшаюших разностные свойства изначальной аппроксимации [9].

С другой стороны, метод решения должен уметь сохранять свойство монотонности в областях с большими перепадами значений. Согласно теореме Годунова (1959) эти два требования одновременно не выполняются в рамках линейных разностных схем. Для преодоления этих трудностей могут применяться разностные методы с выделением разрывов. Эти методы основаны на прямом выделении разрывов, возникающих в численном решении. Выделение достигается путем построения дискретной сетки, связанной с разрывами. В частности может быть использован метод характеристик [10]. Для того чтобы вблизи разрывов не возникали нефизические осцилляции, требуется применение либо монотонных схем первого порядка точности, либо искусственной вязкости. B областях гладкости решения требуется использование схем более высокого порядка точности. Для этих целей были применены т.н. гибридные разностные схемы, или, по-другому, разностные схемы переменного порядка точности.

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

Борис и др.[13-15] развили гибридный метод, который позволяет повысить порядок точности численных расчетов путем применения специальной процедуры коррекции потоков – FCT (Flux Corrected Transport).

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

Заметим, что эти условия эквивалентны условию неувеличения полной вариации численного решения, или условию TVD (Total Variation Diminishing). Для того чтобы решение удовлетворяло TVD-условию, развита специальная техника кусочно-линейной реконструкции сеточных функций.

Эта техника впервые была предложена в работах Колгана [16] и, позднее, развита в работах Ван Лира [17]. Наклоны кусочно-линейных распределений сеточных функций внутри дискретных ячеек для выполнения TVD-условия ограничиваются специальными ограничителями для обеспечения свойства минимальной вариации решения [18, 19].

Практически все схемы переменного порядка точности основаны на кусочно-полиномиальной реконструкции дискретных сеточных функций, удовлетворяющей TVD-свойству. Бурное развитие TVD-алгоритмов привело к созданию ряда их многочисленных модификаций, например TVB [20] (Shu, 1987), AUSM [21] (Liou, 1996), WENO[22] (Liu, Osher, Chan, 1994), WAF[23] (Toro, 1989). Создание этих методик привело к значительному повышению качества получаемых численных решений по сравнению с классическими разностными схемами фиксированного порядка точности. В частности, cхемы ТВД второго-третьего порядка хорошо зарекомендовали себя при расчете ударных волн. Однако для задач, связанных с линейным переносом, повышенного (начиная с 5-го) порядка аппроксимации, сочетающиx конечно-элементные методы высoкого порядка с Римановскими солверами (Разрывный Галёркин) [24] и использующих вместо лимитерной фунцкции переменный шаблон (ЕНО/ВЕНО) [25], этот недостаток практически устранен.

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

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

Схема Кабаре была впервые представлена для линейного уравнения переноса в трехслойном виде в работах Aйзерлиса [26] и Самарского и Головизнина [27, 28] в двухслойном виде, с разнесенными переменными по пространству и по времени, в работе Головизнина и др. [29]. Консервативная нелинейная коррекция схемы Кабаре была впервые рассмотрена в работе Головизнина и Карабасова [30]. Более совершенные варианты алгоритма нелинейной коррекции для уравнений в одномерном случае были предложены в работах Головизнина и др. [29] [31]. В неконсервативном виде схема Кабаре была обобщена для двумерных уравнений адвекции в работах Роу [32], Трана и Шорера [33] и Кима[34]. В компактной форме предикторкорректор, соответствующей сеточному шаблону в одну пространственновременную ячейку, схема Кабаре для гиперболических законов сохранения была представлена в работах Головизнина [35] и Карабасова и Головизнина [36]. Для уравнений двумерного линейного переноса консервативное обобщение схемы Кабаре с учётом нелинейной коррекции было предложено Головизниным и др. в работе [37]. Для двумерных уравнений Эйлера схема Кабаре была представлена Карабасовым и Головизниным [38] а также Головизниным и др. в работе [39]. В последней работе, так же, как и в работе Фараносова и др.[40], посвященной oбобщению схемы Кабаре на трехмерные уравнения Навье-Стокса с учетом параллельной реализации алгоритма, содержится cущественный вклад автора настоящей диссертации.

B рамках настоящей диссертации предложено oбобщение схемы Кабаре на уравнения газовой динамики в двумерном и трехмерном случае с улучшенной процедурой реконструкции потоков в случае звуковой точки.

Проведена серия двумерных расчетов изолированных вихрей и их взаимодействия с ударными волнами, на которых подтверждены такие положительные качества алгоритмов Кабаре, как малая диссипативность и отсутствие паразитных осцилляций при наличие больших градиентов решения. На основе схемы Кабаре разработаны параллельные алгоритмы для решения трехмерных уравнений Навье-Стокса. Эти алгоритмы иcпользованы для решения задачи о моделировании обтекания обратного уступа турбулентым потоком и истечения высокоскоростной затопленной турбулентной струи (эксперимент JEAN[41]) в рамках подхода крупных вихрей. Проведена статистическая обработка решения для нескольких сеток, и показано хорошее совпадение расчета с данными эксперимента.

Equation Chapter 1 Section Глава 1. Схема КАБАРЕ для уравнений одномерной газовой динамики.

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

здесь – плотность, u – скорость, p – давление, – удельная внутренняя энергия, E – полная энергия. Уравнение состояния удовлетворяет равенству:

где c – скорость звука.

Свойство гиперболичности этой системы состоит в том, что в областях гладкости она приводится к т.н. характеристической форме:

где Вычислительные алгоритмы, базирующиеся на дивергентной форме записи уравнений (1.1) и сохраняющие свойство дивергентности на разностном уровне называются консервативными [42]. Консервативные разностные схемы, дополненные процедурами монотонизации численного нелинейная коррекция потоков), составляют основу современного парка алгоритмов сквозного счета газодинамических задач [43].

Вычислительные схемы, опирающиеся на представление уравнений газовой динамики в форме (1.3), получили название характеристических [44], или сеточно-характеристических [45]. Характеристические алгоритмы по точности заметно превосходят схемы сквозного счета в областях гладкости решения, однако при возникновении сильных разрывов теряют однородность обстоятельство обуславливает относительную редкость их практического использования, несмотря на неослабевающий интерес к ним со стороны математического сообщества.

В работе [46] предложен новый вычислительный метод, объединяющий достоинства характеристических и консервативных разностных схем. Его отличительной особенностью является наличие двух видов переменных, т.н.

«консервативных» и «потоковых» величин. Специализация переменных позволяет использовать дивергентную форму записи уравнений для вычисления «консервативных» переменных и характеристическую – для определения потоков. Новый метод, названный «баланснохарактеристическим» заметно превосходит по всем параметрам известные TVD-, TVB- и ENO – алгоритмы [43].

ОПИСАНИЕ СХЕМЫ КАБАРЕ

Сеточные множества и сеточные функции узлов этой сетки обозначим как множество центров ячеек – как. Центры ячеек будем отмечать полуцелыми индексами k+1/2. Множества сеточных функций, определенных на и будем обозначать соответственно как и Для объединения балансного и характеристического подходов будем использовать как множество узлов, так и множество ячеек расчетной сетки.

Введем два типа дискретных величин: т.н. «консервативные» величины давления, k 1/2 - удельные внутренние энергии, Ek 1/2 - полные энергии ячеек, uk - потоковые скорости, pk - потоковые давления, k - потоковые плотности.

Таким образом, количество неизвестных сеточных функций в новом подходе в случае одного пространственного измерения в два раза превышает их число в традиционных подходах [43] и [2].

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

где чертой сверху обозначается осреднение по двум временным слоям:

Вычисление новых значений потоковых переменных Для того, чтобы разрешить систему (1.5) относительно новых значений консервативным переменных, определим неизвестные значения потоковых величин, используя характеристическую форму записи исходных уравнений (1.3).

Вначале вычислим значения консервативных величин на момент времени tn1/2 tn 0,5 n1/2 по разностной схеме аппроксимирующей дивергентную форму дифференциальных уравнений (1.1) со вторым порядком точности по пространству и с первым порядком по консервативных переменных вычисляются при этом со вторым порядком систему (1.3) системой линейных уравнений:

где in1/ Используя обозначения (1.7) можно записать в виде:

Величины r, q, s будем называть «римановыми квазиинвариантами», или просто «квазиинвариантами».

По известным значениям потоковых переменных в момент времени tn и консервативных переменных в момент tn1/2, вычислим соответствующие им значения Римановых квазиинвариантов:

Если в ячейке Din1/2 квазиинварианты являются гладкими функциями, то с точностью O n1/2, hi 1/2 должны выполняться соотношения:

связывающие неизвестные значения квазиинвариантов на новом временном слое с их известными значениями посредством линейной экстраполяции.

Примем эти соотношения в качестве определяющих.

Из уравнений (1.7) следует, что величины квазиинвариантов на новом временном слое должны удовлетворять неравенствам:

где Значения квазиинвариантов на новом временном слое, вычисленные по формулам линейной экстраполяции (1.12) не всегда будут удовлетворять этим неравенствам, поэтому их следует откорректировать в соответствии с (1.13).

В случае дискретных сеточных функций заменим ограничения (1.14) их приближенными сеточными аналогами:

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

Таким образом, для каждого узла с индексами i, n 1 мы получили использовать для нахождения новых значений потоковых переменных, зависит от того, с какой стороны соответствующая характеристика приходит в узел i, n 1. Вычислим «прикидочные» значения uin1 и новой скорости звука cin1 со вторым порядком точности по формулам:

и определим соответствующие характеристические скорости как:

Следует отметить, что на этом этапе нам потребуется только информация о направлении прихода характеристик в данный узел, поэтому вопрос о точности вычисления характеристических скоростей не является критическим.

Возможны всего четыре комбинации прихода характеристик в узел (i,n+1) Новые значения потоковых переменных в первом случае находятся из системы уравнений:

откуда в первом случае получаем:

для сверхзвукового течения влево:

На этом этап вычисления новых значений потоковых переменных во внутренних узлах расчетной сетки завершается.

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

Звуковые точки Рассмотрим все возможные случаи перехода от дозвуковых течений к сверхзвуковым. Для начала определим критерий наличия звуковой точки в узле сетки. Пусть нам известны значения чисел Маха в смежных ячейках на момент времени tn+0.5dt, обозначим их как ML, MR, соответственно в левой и в правой ячейках от узла. Критерием наличия звуковой точки может служить