W-модификация метода Годунова и ее приложения в моделировании газодинамических течений с ударными волнами

СОДЕРЖАНИЕ
ВВЕДЕНИЕ ................................................................. 4
1. ТЕОРЕТИЧЕСКИЕ ОСТОВЫ И СОДЕРЖАНИЕ Ш-МОДИФИКАЦИИ МЕТОДА ГОДУНОВА ДЛЯ ГИПЕРБОЛИЧЕСКИХ СИСТЕМ УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ ......................................... 10
1.1. Общие понятия и определения. Задача Римаиа. Формулировка схемы Годунова. Подход к повышению точности схемы Годунова .................. 10
1.2. Линейные и нелинейные модификации схемы Годунова второго порядка для модельного уравнения переноса. Условие монотонности.
Второе дифференциальное приближение. Примеры функций-лимитеров ., 15
1.3. Нелинейная модификация схемы Годунова для двумерного уравнения переноса. Ш-шаблон вычисления поправок. Теорема о монотонности
в двумерном случае. Примеры функций осреднения гш^а, 6) ........ 23
1.4. Обобщение \У-модификации для системы уравнений в одномерном
и двумержхм случаях .............................................. 27
2. И'-МОДИФИКАЦИЯ МЕТОДА ГОДУНОВА НА ПОДВИЖНЫХ КРИВОЛИНЕЙНЫХ СЕТКАХ ПРИМЕНИТЕЛЬНО К ДВУМЕРНЫМ НЕСТАЦИОНАРНЫМ ТЕЧЕНИЯМ ИДЕАЛЬНОГО ГАЗА...............................29
2.1. Дифференциальные и интегральные формы уравнений движения.
Расчетные формулы для подвижной четырехугольной ячейки ........... 29
2.2. Процедура ускорения численного решения задачи Р и мал а
о распаде произвольного разрыва....................................33
2.3. Вычисление ^’-поправок к методу Годунова для уравнений Эйлера на подвижных криволинейных сетках. Тестовые расчеты и оценка порядка точности \¥-метода ........................................... 35
2.4. Перемещение границы подвижной сетки, связанной с ударной волной --44
3. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ АВТОМОДЕЛЬНЫХ ТЕЧЕНИЙ
С УДАРНЫМИ ВОЛНАМИ ................................................... 47
3.1. Истечение сверхзвуковой струи из расширяющейся щели. Влияние скорости расширения щели и начального перепада давлений на характер течения. Периодические пульсации течения при медленном расширении щели .................................................... 49
3.2. Распад комбинированного разрыва около прямого угла. Режим
течения с периодическими пульсациями ........................... 79
3.3. Маховское отражение сильной ударной волны от клина. Различные случаи отр?-жений. Пульсации спутной струйки в двойном маховском отражении ............................................................ 85
3.4. Маховское отражение слабой ударной волны от клина в условиях
парадокса Неймана, Методика численного моделирования и результаты
вычислений. Четырехволновая схема течения. Сопоставление
численных, теоретических и экспериментальных результатов.........110
4. ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ СТАЦИОНАРНОГО ОТРАЖЕНИЯ КОСОЙ УДАРНОЙ ВОЛНЫ ОТ ПЛОСКОЙ ПОВЕРХНОСТИ ........................ 1*27
4.1. Некоторые особенности численного моделирования ................ 130
4.2. Влияние подветренного угла порождающего клина и внешнего противодавления ................................................... 1-34
4.3. Численное исследование неоднозначности стационарного отражения
косой ударной волны. Эффекты гистерезиса ................... 152
4.4. Эффект гистерезиса при вариации внешнего противодавления ...... 162
5. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ НЕСТАЦИОНАРНЫХ ТЕЧЕНИЙ ЗАПЫЛЕННОГО ГАЗА .................................................. 171
5.1. Уравнения движения запыленного газа. Трехэтапная разностная схема численного решения ................................................. 173
5.2. Численное моделирование взаимодействия плоской ударной волны
с пылевым облаком .............................................. 178
5.3. Течение запыленного газа в сверхзвуковых соплах и струях. Постановка задачи. Методика построения сетки а реализация граничных условий. Результаты численного моделирования ................................ 185
ЗАКЛЮЧЕНИЕ..........................................................;.... 199
СПИСОК ЛИТЕРАТУРЫ...................................................... 202
3
ВВЕДЕНИЕ
Исследование различных течений жидкости и газа с помощью моделирования на ЭВМ стало в настоящее время признанным и быстроразвивающимся направлением пауки, получившим название вычислительная аэрогидродинамика. Эта дисциплина состоит из теории механики жидкости и газа, вычислительных методов, математического и программного обеспечения ЭВМ.
Численное моделирование аэрогидродинамических явлений относится к разряду наиболее ресурсоёмких компьютерных задач. Стремительный рост быстродействия и памяти ЭВМ не может в полной мере удовлетворить потребностям практики и научных исследований, опережающий рост запросов которых во многом стимулируется именно прогрессом вычислительной техники. По этой причине совершенствование имеющихся и создание новых экономичных и быстродействующих численных методов неизменно является актуальной и первостепенной задачей вычислительной аэрогидродинамики.
Для серийных практических расчетов наиболее важными требованиями являются относительная простота и надежность численного метода. Для газодинамических течений с ударными волнами первым методом, удовлетворяющим указанным критериям, по праву считается метод, основанный на "расиадной” схеме С.К. Годунова [34]. Гибкость метода, сочетающая устойчивый сквозной счет с возможностью использования подвижных сеток и простых процедур выделения разрывов, определила его широкое распространение не только в России. Значительную роль в популяризации метода сыграла вышедшая в 1976г. первая в своем роде коллективная монография [36].
Последующее развитие явной схемы Годунова связано с появлением монотонной модификации Колгана [43], имеющей второй порядок по пространственной переменной, идеи которой далее были развиты и обобщены в целом ряде отечественных [44, 47, 26, 27, 64, 65, 57, 33, 4] и зарубежных работ [162, 163, 105, 92, 157, 174, 106, 175, 93].
Монотонные схемы для одномерной нестационарной газодинамики изучены достаточно полно. К 90-ым годам работы по монотонным методам повышенной точности переместились в актуальную плоскость создания и исследования схем для нестационарных многомерных течений однофазных, двухфазных и вязких сред. С повышением сложности математических моделей неизбежно возрастает и сложность разностных схем. В связи с этим повышается роль тестирования методов и их апробации на примере классических двумерных нестационарных задач, в которых еше имеется целый спектр собственных нерешенных проблем (например, задача о маховском отражении ударных волн).
Целями настоящей работы являются:
— создание и тестирование эффективных полностью двумерных монотонных разностных схем повышенной точности для расчетов нестационарных плоских к осесимметричных течений однофазных и двухфазных сжимаемых сред
— исследование пульсирующей неустойчивости автомодельных решений в некоторых течениях струйного типа с ударными волнами
— комплексное численное исследование нестационарного маховского отражения сильных ударных волн от плоского клина
— численное моделирование с высокой точностью маховского отражения слабой ударной волны от плоского клина и создание обоснованных теоретических схем течения около тройной точки в условиях парадокса Неймана
— численное моделирование стационарного отражения косой ударной волны от плоскости и исследование гистерезисных эффектов
— численное моделирование и исследование течений запыленного газа в соплах и перерасширенных струях, включая нестационарную фазу формирования течения
Разработка новых численных методик осуществляется на базе известного метода С.К. Годунова. При исследовании свойств разностных схем используется метод дифференциального приближения [69]. Применяется теоретическое обоснование монотонности разностной схемы для модельного двумерного уравнения переноса. Численные расчеты проводятся с помощью созданной W-модификапии метода Годунова на криволинейных подвижных сетках с применением как известной [48], так и модифицированной процедуры выделения ударных волн.
Структура и объем работы. Диссертация состоит из введения, пяти глав, заключения и списка литературы; содержит 213 стр., включая 86 стр. с рисунками и 12 стр. списка литературы. В работе 106 рисунков и 178 библиографических ссылок.
— В главе 1 введеды-основные понятия для гиперболических систем уравнений в частных производных. Сформулирована задача Римака в виде частного случая задачи Коши с кусочно постоянными начальными условиями и выписаны ее решения в частных случаях. Сформулирована разностная схема С.К. Годунова в двумерном случае, как явная двухслойная схема, использующая для определения числовых потоков решение задачи Римана о распаде разрыва на границах ячеек. Далее излагается подход к повышению порядка точности схемы Годунова путем коррекции числовых потоков с помощью введения аддитивных поправок к исходным аргументам задачи Римана. Выписаны выражения для поправок и обоснован второй порядок аппроксимации на гладких решениях. На примере модельного уравнения переноса рассматриваются различные способы вычисления поправок в виде нелинейных средних конечных разностей на трехточечном шаблоне. Вводится функция выбора среднего mid(a,6), линейным разновидностям которой соответствуют различные классические немонотонные схемы второго порядка. Выписано достаточно общее условие на функцию mid(a,6) и компактно доказана теорема, о том, что это условие обеспечивает монотонность. Для функции осреднения вводится порождающая функция одного аргумента (лимитер v>(0 )> для которой выведены условия монотонности, аналогичные условиям Roe и Sweby. В общих предположениях гладкости ср{£) выведено вто]х>е дифференциальное. приближение для схемы с поправками. Рассмотрен вопрос, о влиянии лимитера на стационарные решения второго дифференциального приближения и показано преимущество лимитеров, обладающих свойством ¥>"(1) > 0. Предлагаются новые разновидности гладких адаптивных лимитеров для схем второго и третьего порядков.
На примере двумерного уравнения переноса конструируется полностью двумерная монотонная схема второго порядка с вычислением поправок на переменном \¥-шаблоне, ориентация которого зависит от знаков коэффициентов уравнения (направление скорости переноса). Выписано достаточное условие на функцию осреднения т»с1(а,Ь) к доказана, теорема о монотонности в двумерном случае. Приводится обобщение \У-модификации на случай гиперболических систем уравнений с компактной вектор но-матричной формой записи алгоритма вычисления поправок.
Глава 2 посвящена подробному описанию \У-модификации метода Годунова для нестационарных уравнений газовой динамики и ее тестированию. На основе интегральных законов сохранения выводятся формулы для пересчета параметров в отдельной подвижной расчетной ячейке. Выписаны беэытерационные формулы для осесимметричных течений в случае неявной аппроксимации Источниковых членов. Приведены формулы для вычисления поправок, обеспечивающих повышение порядка точности для общего случая криволинейных подвижных сеток. Алгоритм поправок увеличивает затраты компьютерного времени на практических задачах в среднем на 45% яо сравнению с исходным методом Годунова. Предложена, ускоряющая процедура приближенного решения задачи о распаде произвольного разрыва с контролем точности, которая сокращает общие затраты времени СР11 на решение задачи Римапа в среднем в 2.5-3 раза, при этом суммарные затраты времени всей программы сокращаются на 25-20%, 'го есть компенсируется около половины затрат, идущих на повышение порядка. Представлены результаты расчетов двух тестовых задач, которые демонстрируют эффективность метода на примере сложного разрывного течения и второй порядок точности »га примере осесиметричного непрерывного течения. В последнем разделе главы описаны две процедуры выделения ударных волн, используемые в последующих расчетах в случае подвижных сеток, одна из которых основана на принципе Гюйгенса (методика Крайко, Макарова и Тилляевой), а другая является ее улучшенной модификацией.
В главе 3 технология подвижных сеток с выделением наиболее важных разрывов используется для численного решения некоторых двумерных нестационарных задач, допускающих автомодельное решение. В первом разделе исследуется задача об истечении сверхзвуковой струи газа из одного полупространства в другое через щель, размер которой увеличивается с постоянной скоростью. Здесь при некоторых исходных данных обнаружена неустойчивость автомодельного течения. Вместо автомодельного формируется пульсирующее течение с постоянным периодом в логарифмической шкале времени. Пульсации аналогичной природы обнаружены также в течении, возникающем при распаде комбинированного разрыва около кромки прямого угла. Результаты численного решения этой задачи представлены во втором разделе главы. Третий и четвертый раздел посвящены явлению нестационарного маховского отражения ударной волны от клина. Исследование проводятся с целью ответить на ряд вопросов, обозначенных в проблемной статье Веп-Иог к Такауата (1992). Приведены результаты численного решения более 50-ти вариантов в диапазоне сильных ударных
волн. Для четырех вариантов проведено сравнение результатов численного моделирования с экспериментальными интерферограммами. Выявлены некотрые закономерности для интенсивности стебля Маха и угла траектории тройной точки. Обнаружена неустойчивость тангенциального разрыва и струйки в двойном м&ховском отражении и установлена определяющая роль струйки в раскачке течения. С помощью процедуры принудительной стабилизации течения выяснен характер взаимодействия второго стебля Маха с тангенциальным разрывом в двойном маховском отражении. Обнаружен переходный режим между двойным и тройным маховским отражением. В последнем разделе изучается отражение слабых ударных волн в диапазоне парадокса Неймана. С помощью уникальной численной методики удалось выяснить природу парадокса и состроить четырехволяовую схему течения в тройной точке с центрированной волной разрежения за отраженным фронтом, которая органически дополняет трехударную схему Неймана и снимает теоретический аспект парадокса. Показано, что в диапазоне слабых взаимодействий отраженный фронт имеет логарифмическую особенность. Приводится сравнение численных, теоретических и экспериментальных данных.
Глава 4 посвящена численному моделированию стационарного отражения косой ударной волны от плоской поверхности. Исследовано влияние подветренного угла порождающего клина и внешнего противодавления за клином на характер отражения падающего скачка и формирующееся за ним течение. Исследован и подтвержден эффект гистерезиса по углу падающей волны. Показано, что выделение фронта падающей волны создает побочный численный эффект более раннего перехода регулярного отражения в маховское и сужает область гистерезиса. С помощью методики выделения линии срыва рассчитаны зависимости для минимальной высоты канала, при котором возможно формирование стационарного течения. Исследовано влияние противодавления на переходные процессы и обнаружен эффект вынужденного гистерезиса в переходах между регулярным отражением и инверсным маховским отражением. Найдены две разновидности переходных процессов ЛИ.—»1пМИ: с промежуточной стадией в виде переходного регулярного отражения (ТЛИ) и без промежуточной стадии. Определены границы между этими разновидностями. Установлено, что переходы в обоих направлениях сопровождаются скачкообразной перестройкой структуры течения.
Глава 5 посвящена развитию методов расчета и численному моделированию двухфазных течений в соплах и струях. В первом разделе в рамках единого подхода повышения точности с использованием \\г-шаблон а построена трехэтапная разностная схема второго порядка для двумерных нестационарных течений в модели запыленного газа. Во втором разделе обсуждаются возможности метода на примере численного решения задачи о взаимодействии ударной волны со сферическим облаком пыли. Проводится анализ эволюции облака. Третий раздел посвящен численному моделированию формирования двухфазного течения в сопле и в перерасширенной струе за срезом сопла. Обсуждается предыстория вопроса. Изложены детали численной методики, касающиеся построения сетки и реализации граничных условий. Изучено влияние размера частиц, начальной загрузки и коэффициента сопротивления на характеристики
стационарного течения в сопле и на конфигурацию скачков в перерасширенной струе. Обсуждается влияние частиц на нестационарную фазу формирования течения в струе и в дозвуковой секции сопла.
Научная новизна результатов диссертации состоит в следующем:
— предложена новая полностью двумерная \\>-модификация схемы Годунова второго порядка точпости по пространству и времени, на базе которой построены аффективные методы рачета нестационарных однофазных и двухфазных течений сжимаемого газа
— предложены: новые виды гладких функций-лимитеров, удовлетворяющие условиям монотонности и значительно повышающие степень локализации слабых ударных волн и контактных разрывов; процедура быстрого решения задачи Рим ада с контролем точности; модифицированная процедура выделения ударной волны, учитывающая скорость распространения возмущений за фронтом волны
— исследована новая задача об истечении газа из расширяющейся щели, моделирующая начальную стадию раздвижения диафрагмы в ударной трубе, обнаружено, что для режимов медленного раздвижения имеет место неустойчивость автомодельных решений* вместо которых наблюдаются пульсации течения с периодом равным времени 4-х — 5-ти кратного увеличения размера щели
— обнаружены аналогичные пульсации в течении, возникающем при распаде комбинированного разрыва с сектором плотного газа около кромки прямого угла
— в нестационарном маховском отражении установлены некоторые закономерности для зависимостей интенсивности стебля Маха и угла траектории тройной точки от исходных данных, обнаружена неустойчивость тангенциального разрыва и спутнон струйки в двойном маховском отражении
— для двойного маховского отражения впервые обнаружен каскадный тип взаимсь действия второго стебля Маха с тангенциальным разрывом, который характеризуется чередующейся последовательностью волн сжатия и разрежения с убывающей интенсивностью
— в области двузначности решения обнаружены новая разновидность двойного маховского отражения и переходный вид между двойным и тройным маховским отражением
— обнаружено присутствие центрированной сверхзвуковой волны разрежения за отражепным фронтом для слабого маховского отражения, построена новая четы-рехволнозая теория течения в окрестности тройной точки, которая дополняет трехударную теорию и устраняет теоретический аспект парадокса Неймана
— выведено уравнение для формы отраженного фронта вблизи тройной точки и показано, что отраженный фронт имеет логарифмическую особенность
— классифицированы три вида отражений слабых ударных волн с особенностью: слабое неймановское, сильное неймановское и слабое махов ское, причем показано, что для двух последних характерный размер особенности более чем на порядок меньше пределов разрешимости эксперимен тальных измерений в настоящее время
— установлено, что предельные конфигурации ударных волн и параметры течения непосредственно в тройной точке для всех видов отражений описываются трехударной и четырех волновой теориями
— для стационарного маховского отражения в канале подтвержден эффект гистерезиса по углу падающей волны, обнаружен и исследован эффект гистерезиса в переходах между регулярным и инверсным маховским отражением при вариации противодавления, найдены две разновидности переходных процессов и определены границы между ними
— изучено влияние размера частиц и начальной загрузки на формирование течения запыленного газа в сопле и в перер&сширенной струе за срезом сопла, установлено, что влияние частиц 4 малого” и "большого" размеров приводят к смешению диска Маха в струе в противоположных направлениях, обнаружены пульсации в струе из-за неустойчивого положения диска Маха, показано, что присутствие пыли стабилизирует течение в струе и уменьшает время формирования течения внутри сопла.
Научная и практическая значимость работы состоит в создании эффективных численных методов для нестационарных однофазных и двухфазных течений с ударными волнами. Изложенные подходи распространимы на более сложные модели и на трехмерные течения.
Некоторые особенности рассмотренных в диссертации газодинамических течений имеют место в важных практических приложениях. Так, характер и механизм обнаруженных пульсаций автомодельных течений имеет общую природу с пульсациями, возникающими при сверхзвуковом обтекании тел со струей, направленной навстречу потоку. Исследованные гистерезисные явления по противодавлению, приводящие к резкой перестройке структуры течения, могут возникать в перерасширенных струях за соплами с коническим насадком. Стабилизирующее влияние второй фазы в перерас-ширенной струе весьма существенно для коротких сверхзвуковых сопел.
Таким образом, результаты диссертации могут быть использованы для оптимизации работы сопловых устройств реактивных двигателей на различных эксплуатационных режимах.
Теоретические результаты диссертации входят в спецкурс "Численные методы для гиперболических систем”, читаемый автором на математическом факультете Волгоградского университета.
Глада 1
ТЕОРЕТИЧЕСКИЕ ОСНОВЫ И СОДЕРЖАНИЕ \У-МОДИФИКАЦИИ # МЕТОДА ГОДУНОВА ДЛЯ ГИПЕРБОЛИЧЕСКИХ СИСТЕМ
УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ
1.1* Общие понятия и определения. Формулировка схемы Годунова для линейных систем и нелинейного уравнения. Подход к повышению точности схемы Годунова
Широкий круг задач в механике сплошных сред и, в частности, в газовой динамике описывается системами уравнений гиперболического типа. Рассмотрим систему квазилинейных дифференциальных уравнений относительно т-мерной вектор-функции и(*,Т,у) = (щ,и2у...,Х1т)Т
дп .дп , пди , .
dt* дх ду~ ' ( }
где вещественные (го х тп) матрицы Л, В и правая часть h в общем случае зависят от iy X, у и и. Если h(u) линейна по и, и матрицы Ау В не зависят от и, то имеем линейную систему с переменными коэффициентами, а в случае А, В — const — линейную систему с постоянными коэффициентами.
Систему (1.1) с А} В = const, будем называть гиперболической (35, 56), если собственные вектора матрицы А (а также и В) образуют базис в Rm. Или другими словами существуют невырожденные вещественные матрицы Ra и Rb такие, что
A = RaARa\ B = RbMR£\ (1.2)
где Л = (Аь...,Ат) и М = (ді,...,дт) - диагональные матрицы собственных значений матриц А и В соответственно.
Из (1.2) видно, что столбцы матриц RA и Rb являются собственными векторами матриц А и Д. В отличие от других данное определение не требует, чтобы собственные значения были различными.
Для линейных систем с переменными коэффициентами свойство гиперболичности зависит от точки (і,Хуу)у а для квазилинейных систем также и от решения u(t, ж,у).
В ряде случаев система (1.1) может быть записана в дивергентной форме
Sr+T+T-1- <‘J>
Такие системы называются консервативными. В зарубежной литературе принят термин conservation laws (законы сохранения), так как именно в подобном виде в газовой динамике записываются законы сохранения массы, импульса и энергии среды. Вектор-функции f(u) и g(u) в (1.3) называются функциями потоков, а правая часть h(ii) — источшїковьім членом.
Интегральным аналогом системы (1.3) является уравнение
j>(\inK + fnx ~ yhrfv, (1.4)
av v
10
где V - произвольный гомеоморфаый сфере объем в трехмерном пространстве (*, я, у), (п^п^Пу) - компоненты внешней нормали к поверхности объема дУ. Негладкие или разрывные решения интегрального уравнения (1.4) называются слабыми решениями системы (1.3).
Основной идеей схемы С.К. Годунова [34, 36] является непосредственное использование в разностной схеме точного решения задами Римана о распаде одномерного разрыва. Задача Римана - это задача Коши для (1.3) с кусочно-постоянными начальными условиями в виде скачка
где /(ас,у) = 0 - уравнение положения прямолинейного фронта скачка в момент t = 0. Очевидно, что в такой постановке под решением подразумеваются слабые решения (1.3), т.е. удовлетворяющие уравнению в интегральном смысле (1.4), При отсутствии Источниковых членов задачу (1.3) с начальными условиями (1.5) можно свести к одномерной задаче, выбрав пространственную координату перпендикулярно фронту разрыва.
Запишем задачу Римана для одномерного случая
ди д((и) (и1, х<х*
В скалярном случае ( т = 1, /(и) = /(«) ), а также в случае линейной системы с постоянными коэффициентами ( /(и) = Ли ) решение задачи (1.6) можно довольно просто представить в виде одного или нескольких эволюционных разрывов, между которыми решение постоянно. При этом решение в точке разрыва х* при < > 0 выглядит следующим образом
если т = 1, то
,,// гЧ - + цД л.*т( Я“*)-/(ц*)\
«(<>*)- 2 + 31811 ^ иг,_ик ^ 2 .
если /(и) = Ли, то =----------+ з1§п(Л)-----------,
здесь функция одп от матрицы Л : эдп(Л) = в!до(Л)Лд\ где sigи(Л) - диаго-нальная матрица знаков собственных значений матрицы Л. Отметим, что в нелинейном случае решение (1.6) может быть не единственным и для устойчивости решения в виде одиаочного эволюционного разрыва иЦ^х) = п(х — ИЬ) необходимо выполнение энтропийного условия [144, 132].
В более общем случае (1.6) по аналогии с предыдущими можно выписать линейное приближение для и(/, х*) в виде
~ и*1 — ил -*'•
«(!,*•) = —2— + »*«"(/-)--------2---» где /« = ^|и=1(„Ч.и »)•
Схема С.К. Годунова [36] принадлежит к классу явных двухслойных разностных схем. Пусть сетка выбрана равномерной по координатам (х, у), причем, направление роста индексов ячеек (у,^) совпадает с направлением осей координат. Обозначим через
11
й,ч вектор параметров в ячейке і) в момент Ь 4- Аі. В этом случае схема Годунова первого порядка для (1.3) запишется в виде
Щі-Uji *(4,44) - g(ut+p - S(u.-^)
At Ax Ay
- h(u^) = 0. (1.8)
Величины Ци^), £»(и«±£) называются числовыми потоками. Аргументы и и^т. определяются из решения задачи Римана на границах ячейки. В линейном приближении (1.7) они имеют вид:
Основная идея, большинства двухслойных схем повышенной точности заключается в консервативной поправке числовых потоков для какой-либо известной схемы первого порядка (в том числе и для схемы Годунова). Нелинейный характер поправок, а также переменный шаблон их вычисления позволяет обойти известную теорему Годунова [34] и совместить второй порядок точности схемы с монотонностью.
Условно можно выделить три основных различных подхода., применяемых при построении схем высокого разрешения:
1) использование поправок для числовых потоков, определяемых как разности потоков для монотонной схемы первого порядка и схемы более высокого порядка, с последующим их ограничением из соображений сохранения монотонности [84, 162, 177].
2) использование поправок к исходным данным для задачи Римана, определяемых путем более точного представления распределения параметров внутри ячеек на каждом временном шаге (43, 163, 47], а также использование обобщенной задачи Римана с линейным профилем исходных данных [75].
3) использование аддитивных поправок к потоковым функциям [105] или к аргументам потоковых функций [93] в решаемых уравнениях с последующим применением схемы первого порядка.
Подробный обзор работ о монотонных схемах в хронологической последовательности до 1990 года, описание, и классификация наиболее известных разностных схем приведены в [129, 130].
Неоспоримые преимущества нелинейных монотонных схем убедительно продемонстрированы для одномерных законов сохранения. Наиболее распространенный подход в повышении точности методов для многомерных нестационарных задач основан па сочетании схем расщепления и одномерных схем повышенной точности. Методов, не использующих расщепления (unsplit schemes) или так называемых полностью многомерных (full multidimensional) методов повышенной точности мало и их возможности еще не раскрыты.
В данной работе выбрана линия разработки полностью многомерной монотонной модификации метода Годунова повышенной точности. Идея модификации состоит в комбинации подходов [163] и [105] и в некоторых чертах схожа с [93].
Схема Годунова первого порядка (1.8-1.9), в.которой используется второй порядок аппроксимации источниковых членов, будучи примененной к модифицированному
(1.9а)
(1.96)
уравнению
у + а.Пи + а) + |_8(и + 0) = ь> (М0)
даст приближенное решение исходного уравнения (1,3) со вторым порядком точности по пространству и времени, если поправки
Ах . Дг Л Дд у V Д<
а = —+ уи(, $ = — 8!ёо(8и)и, + —и,. (1.11)
Покажем это подробнее. Разностный оператор схемы Годунова для уравнения (110) со вторым порядком аппроксимации источниковых членов запишется в виде:
, г(имНК-|) , 8(Ц.+р-8(^-р 1і(ц,.)+Ь(а,ч)
У }~ Д* Ах Ау 2 К *
где величины и вычисляются по формулам (1-9), но с соответствующими поправками:
и,±, = + ± 8І£п(г,)(»+а>* (1ЛЗв)
„а1, 1 ±8Іеп(^)(^к^^, (1.Ш)
і А
Пусть и(£,ж,у) и 1(и) - гладкие вектор-функции. Тогда
=г“(“-ти>-*■)(“*•* ■ и>ч) + °(дг3)' <114)
в свою очередь
2и* + иі+1>і + - иі+1,і
и?+^ + и;_^ =---------2 + 2-
+ ^ + аі^£+_°±±і + 8;8п(£и)д^Д1±М =
I* &
= [2 и - ет§п(£,) Дхих]д- + [2а — зійп{Ги) ДхосхЬ% + 0(Дх2).
Последнее, с учетом выражений для поправок а и (3 (1.11), представляется в виде
и-ц + и-_^ = [2 и + АЬи^г + 0(Дх2) = + 0( Дх2, Д*3). (1.15)
Аналогично для разностей:
«», - а,., - ■»; “'-»,<№|Цац-
.~ , 3|6а(Г^2оГл1> ~ ,Д:;" Т =
О ^
д 2 д 2
= [Дхиг - + [Ах<х.х - зівд(См) “ + 0(Ах%
Последнее также с учетом (1.11) представляется в виде
и,+ £ - = Дх[иг + + <?(Дя3) = Дх[и*]я[+^ + 0(Ат3, ДхД*2). (1.16)
13
Объединяя (1.14), (1.15) н (1.16) получим, что
і{иМ)-ЦпН) = (г(ц)і]..|(^ + 0(Дх2) д<3) (1Л7)
Аналогично для потоков по у:
— к(иі— і) г / Ч 1 I гл/ л 2 л 2\
—у — = + 0(Ау\ А*2). (1.18)
Если к этому добавить очевидные разложения:
= + ---у—- = |ь(и)Ь.|,^ + 0(ЛП
то суммируя их с (1.17-1.18) получим, что разностный оператор (1.12) на гладкой функции и (і, х, у) представим в виде:
1(и) = [и, + ^и)* + 8(и)„ - Ь(и))я|1+^, + 0(Лх2, Ау2, Д<3),
т.е. обращается в нуль на решении дифференциального уравнения (1,3) с точностью до величин второго порядка малости. Следовательно, схема (1.12-1.13) имеет второй порядок аппроксимации.
Заметим, что для обоснованности вывода необязательно точное выполнение соотношений (1.11), достаточно, чтобы эти условия выполнялись со вторым порядком по Д* и Ау для ас, /9 и для их производных. Последнее посуществу предполагает определенную гладкость поправок и может быть заменено более сильным требованием выполнимости (1.11) с точностью до (Дт3,Ду3).
Для практического вычисления поправок удобно заменить и« пространственными производными из уравнения (1.3), тогда (1.11) преобразуется к виду:
а = ^зщп(Г„}и„ - ^(и)* - ^8(и)„ + у^Ь(и),
(1.19)
. Ау , Ді Д* Д<
0 = —- уг(и)х - уЄ(и)ї + у-Ь(и).
Создать процедуру приближенного вычисления выражений в правых частях (1.19) с помощью конечных разностей возможно различными способами. Дело осложняется тем, что обычные способы, использующие линейные конечные разности приводят в итоге к немонотонным схемам второго порядка. Далее в этой главе рассматривается нелинейный способ аппроксимация величин сх и (3 на переменном двумерном шаблоне (\\'-шаблои), ориентация которого зависит от направления характеристических скоростей (знаков собственных чисел матриц ^ и и выводятся достаточные условия, которые в случае скалярного линейного уравнения обеспечивает сохранение монотонности решения при том же числе Куранта, что а метод Годунова.
14
1.2. Линейные и нелинейные модификации схемы Годунова второго порядка для модельного уравнения переноса. Условие монотонности. Второе дифференциальное приближение. Примеры функций-лимитеров.
В данном разделе рассмотрим вопросы повышения точности схемы Годунова на примере простейшего одномерного линейного уравнения переноса с достоянным коэффициентом
Эи Эи
+ а — = 05 а- = const. (1.20)
Решения задачи Копти этого уравнения хорошо известны.
Пусть сетка выбрана равномерной, причем направление роста номеров ячеек (j) совпадает с направлением оси координат х . Обозначим через Uj и Uj значения и(а?) в ячейке j в моменты t и < + At соответственно, а через j— и j-f номера ячеек, находящихся соответственно против потока и по потоку по отношению к ячейке j :
j- = j - sign(a), = j -f sign(a).
В этих обозначениях разностная схема Годунова для (1.20) запишется в виде
S, = u, - v(iij - Uj.), где v = |а|^. (1.21)
Известно, что схема (1.21) устойчива и монотонна при выполнении следующих условий
«/>0, 0, . (1.22)
первое из которых обеспечено автоматически в силу использования разностей против потока, которые получаются как результат применения решения задачи Римана для (1.20).
Модифицированную схему в соответствии с подходом, изложенным в предыдущем разделе запишем как
Щ = uj - v{(u + а)} -(и 4- а),-), (1.23)
где поправка а$ согласно формулам (1.19) должна аппроксимировать в j-ой ячейке выражение
^-=-^siga(a)Axg. (1.24)
Будем рассматривать возможность аппроксимации (1/24) на трехточечном шаблоне по значениям щ-и щ, При различных способах линейной аппроксимации будем получать разные известные классические схемы. Так, при использовании разностей вперед по потоку получим схему Лакса-Вендроффа [131] и Мак-Кормака [135] (для линейного уравпения переноса (1.20) они эквивалентны), при использовании разностей против потока получим схему Бима-Уорминга [169, 170], при использовании центральных разностей получим схему Фромма [99].
Для каждой ячейки j введем обозначения для разностей против потока pj
Рз = Uj- и,- (1.25)
15
Из очевидного тождества (уЧ)— = 0’")+ = ^ следует, что («,+ — и}) = р,+ , т.е. разность вперед по потоку в ячейке ) равна разности против потока в ячейке ^Ч , которая находится впереди по потоку по отношению к ячейке } . С учетом этого замечания линейные модификации схемы (1/23), упомянутые в предыдущем абзаце, можно записать следующим образом
1
сх} = -(1 — и) Р}+ схема Лакса-Всуцюффа,
а> = -(1 — у) Р] схема Вшла-Уор.минга.
£
1 2?* 4г л •
от, = -(1 — и) 3 ^ схема Фромма,
а,- = ^(1 - и) ((1 — 0)ру + &р3+) обобщенная схема Фромма 0 < в < 1.
Второе дифференциальное приближение в П-форме [69, 1] для обобщенной схемы Фромма имеет вид
д 2
щ + аих = - *№ ~1'~ зб)~
|а|Ат
з
(1 - о)(2 - Zv + vlArAvO — Щите* + 0(Дя4), (1.26)
о
откуда видно, что при 9 = (2 —1/)/3 дисперсионный член с и„х исчезнет и обобщенная схема Фромма будет иметь третий порядок аппроксимации [99]. Однако все упомяпутые схемы не обладают монотонностью.
Для выяснения свойств монотонности рассмотрим более общий вариант схемы
<х, = i(l - v)mid(pj,pi+), (1.27)
где mid(a, b) - некоторая функция осреднения а и b. Тогда имеет место следующая Теорема 1. Если функция mid(a, Ь) такова, что для любых а, //, с
|mid(a, 6) - mid(c,а) - —Г* ^ Г- \а\у (1.28)
г> = rh' Г2=Ь (L29)
то при условии (1.22) разностная схема (1.23, 1.25, 1/27) устойчива и монотонна.
Доказательство. Неравенство (1.28), имеющее смысл лишь при условии (1/22), можно переписать в виде равенства с некоторым неопределенным коэффициентом 7
mid(a, 6) - mid(c, a) - Га a = 7— — a, —1 ^ 7 ^ 1.
и ь
Тогда для разности поправок из (1.27) с учетом того, что p(j+)- = Pj получим
otj-aj- = l(l-j/)(mid(p,-,p,+)-mid(pj_,p3))= l(l-i/)(''2 ~- + T"-g Г‘)
Последнее соотношение позволяет переписать схему (1.23) в виде, аналогичном схеме Годуиова (1.21):
щ = щ - v(uj - «j_), где v = i/[l + ~(1 - v){r2 - n + 7(^*2 + П))].
16
После элементарных преобразовании с учетом (1.29), а также с учетом диапазона изменения у получим следующие оценки ДЛЯ V
1Ч" у
V = —J1 0 0^1.
Следовательно схема (1.23, 1.25, 1.27) устойчива и монотонна. Теорема доказана. ♦
В соответствии со смыслом функции осреднения потребуем от m>d(a, 6), чтобы она была однородной по своим аргументам и тривиальной в случае их равенства:
mid(fca, kb) = к mid(a, 6), mid(a, a) = a.
Однородность позволяет перейти к порождающей функции одного переменного
mid(o,6) — amid(l,&/a) = a<p(b/a).
Тогда из условия (1.28) следует, что порождающая фзгнкция для любых £ и 77 должна удовлетворять неравенству
“ri ^ р(() “ ^¥>0?) < гз (1-30)
Лалее ограничимся рассмотрением таких <р(£), которые неотрицательны при ( ^ 0 и обращаются в нуль при £ < 0. Последнее условие означает, что функция осреднения обращается в нуль, если аргументы разного знака. Сказанное позволяет из (1.30) выписать следующее условие на функцию <р(£)
Г 0 <¥>({) ^minfa&fa), £>0 . .
Wtf) = o, *<о. (L31)
Впервые подобную функцию практически одновременно ввели в рассмотрение Sweby Р. [157] (с ограничением при гг = г2 = 2) и Roe Р. [149]. Область монотонности (т е. условие (1.31)) и некоторые примеры для изображены на рис. 1.2, откуда, в частности, видно, что линейные функции, кроме <р(£) г 0, условию монотонности не удовлетворяют. Ограничения на функцию (р(() накладывают соответствующие ограничения и на функцию осреднения mid(e, 6). В литературе для у>(£) часто используется термин функция-лимитер.
Далее рассмотрим вопрос о влиянии поведения функции <р(£) на точность схемы и, в частности, на ее дифференциальдое приближение. Так как поправка о (1.27) должна аппроксимировать выражение (1.24), то необходимо, чтобы mid(a,a) = а . Это условие уже было упомянуто. Из него следует у?(1) = 1. Предположим, что функция <?(£) гладкая, и введем обозначения
О' = v'«)|e=l, (С = V>"(0j£=1. = v>'"«)|(=1- (1.32)
Тогда разложение но степеням Ах в точке х, в предположении, что а > 0, приводит к следующим соотношениям
mid(pj,pi+) = Ат(их + ($'- \)Ах ихх + | Ах2 иххх + \ 0”Ах2 ulJux
+(n*'~ £)A*3tw + (^'+ \nb**»lJ< + 0(Ax%
17
«“!<%-, Р>) = Ах(и, + ($'- |)Джи<са. + (| -0')Ахгиххх+ зв"Д**о|./««
+(&*- § )л*3 «„** + (! о" + \(Г)^<К
« -0"Д:С3
«х*«ссл* /«* + 0(Дг4)),
откуда разность
т^(Р;>Р;+) - хгаа(р,--,р>) = Дх2^** + (0' - 1 )Ахи1ХТ + ( £ - £0')Да?а«мад;
- *0"Дх2 и3х/и2 + 0"Дх2 «„и,,*/«. + 0(Дх5)).
Объединим два последних слагаемых в предыдущем выражении и выпишем разложение для разности поправок
<*; “ <*$- =
^-у^Дх2^ + (Р - 1)Ахиххх + (у - ~)Дх2«**** + ^Лх2 (^)#+0(Да?3)^.
Видим, что но сравнению с линейным случаем (0" — 0"' = 0) имеет место еще один член такого же порядка малости как и член с производной , который в свою очередь
является наименьшим учитываемым членом во втором дифференциальном приближении (1.26). Следовательно этот добавочный член войдет просто дополнительным слагаемым в дифференциальное приближение (1.26). Для обобщения на случай а < О нужно только заменить Дх на — Дх, В окончательном виде второе дифференциальное приближение для нелинейного случая с условиями (1.32) на (р(() будет иметь вид (в том числе и для а < 0):
щ + аи9 = ^гЦ'1 - *')(2 -V- 30')«*** ~ (1 - |/)0"(—)
О 4 ' их 1 *
- ^—(1 - *)((1 - И(2 -у)- 2«'(1 - 2*/)) «„„+ 0(Л*4). (1.33)
Видим, что для любых конечных значений 9' к в” модифицированная схема имеет второй порядок аппроксимации. Третий порядок также как для схемы Фромма реализуется при (У = (2 — г/)/3. Четвертый порядок в общем случае недостижим, однако существует целый класс ограниченных решений исходного уравнения, (1.20) на кото-+ рых достигается четвертый порядок аппроксимации.
Дифференциальное приближение для схемы третьего порядка (О' = (2 — 1/)/3) имеет вид
Ъ + а»ж = -Щ£(1-у)9"(Щ -Щ£з(1-у)П1-в’)иг*** + 0(Ах4) (1.34)
4 4 их ' * о
Отбрасывая члены 0(Дх4) и переходя в систему координат, движущуюся со скоростью а, получим уравнение
Ъ / ,,2 \ о0"
и, = -||.|Д*^(1 - УЩ1 - в1) («2* + и,*.) . где 8=ЗЙ,(1_Й,). (1.35)
Функции и(х), обнуляющие правую часть уравнения (1.35), не будут изменять свою форму, т.е. будут язляться стационарными решениями. Эти функции удовлетворяют
18
уравнению 4-го порядка и, следовательно, в общем случае зависят от четырех кон-стает. С другой стороны, прямой подстановкой легко проверяется, что если а(аг) стационарное решение то с-\ и(с2Х + сз) + с4 тоже является стационарным решением, т.е. решение инвариантно относительно растяжений и переносов вдоль осей координат. Таким образом, задача нахождения общего стационарного решения сводится к нахождению частных решений.
В случаях 5 = —1 и 8 = 0 стационарные решения легко находятся в конечном виде. В первом случае это сумма линейной и тригонометрической функции, а во втором — кубические полиномы. Из них можно скомбинировать функцию в виде уединенной волны с гладкой переходной областью. Так, для волны с амплитудой 2 и углом паклопа в пептре переходной области 45°, получим:
На таких функциях схема 3-го порядка с соответствующим значением в" будет иметь четвертый порядок точности (кроме точек сопряжения). Стационарные решения можно рассматривать как предел эволюции в процессе численного решения близких к ним функций, например, волны с линейным начальным профилем.
Для других значений з поиск предельных профилей сводится к специальным функциям. Для некоторых значений з предельные профили изображены на рис. 1.1. Все профили имеют непрерывную первую производную. Вторя производная непрерывна при $ < 0 и разрывна в точках сопряжения при $ > 0, причем, при з > О она обращается в бесконечность. Из графиков видно, что размазывание начального линейного профиля в угловых точках значительно меньше, если з > 0, то есть, если вторая производная 9" функции-лимитера в точке £ = 1 положительна. Этот факт важен для практических расчетов, так как позволяет уменьшить эффект размазывания на границах волн разрежения и обеспечивает отсутствие роста ширины зоны размазывания контактных разрывов. В то же время при слишком больших $ возможен и противоположный эффект генерирования разрывов первой производной решения в гладкой переходной области, например, в окрестности экстремумов.
Тестированию различных монотонизированяых методов повышенною порядка посвящено очень большое количество работ. Однако многие по существу различные методы применительно к линейному скалярному уравнению выглядят одинаково, и их сравнение в этом случае превращается в сравнение свойств применяемых лимитеров [128, 1.40, 153, 178). Большинство описанных в литературе лимитеров используют трехточечный шаблон и представимы в виде (1.27). Более сложный вариант предложен Harten & Osher в так называемой UNO-схеме [106], лимитер в которой использует пятиточечный шаблон.
Наиболее распространенные лимитеры, основанные на трехточечном шаблоне, приведены на рис. 1.2. Самые простые - это кусочно-линейные лимитеры, они легко
-2 < * < 2, х ^ 2,
2,
при s = 0, у(#) = 1 + х(4х2127 — 1),
О,
я ^ -1.5, —1.5 < х < 1.5, х ^ 1.5,