Содержание к диссертации
Введение
1. Математические модели и методика численного моделирования. 11
1.1. Математические модели движения сжимаемого газа. 13
1.2. Повышение порядка точности схемы Годунова. 16
1.3. Методика ускорения решения вспомогательной задачи Римана .
1.4. Сеточные методы погашения нефизических эффектов численного решения за сильными малоподвижными ударными волнами .
1.4.1. Метод искусственного излома сеточных линий. 32
1.4.2. Метод вибрирующей сетки. 41
2. Численное моделирование истечения перерасширенной струи газа из короткого осесимметричного сопла
2.1. Постановка задачи и методика численного решения 51
2.2. Результаты расчетов. 55
2.2.1. Результаты расчетов в модели идеального газа . 55
2.2.2. Результаты расчетов в модели вязкого теплопроводного газа.
3. Численное моделирование распространения ударной волны в канале с препятствием .
3.1. Постановка задачи и методика решения. 95
3.2. Результаты расчетов. 100
3.2.1. Начальный этап формирования течения без учета вязких эффектов .
3.2.2. Результаты расчетов в модели Навье-Стокса с учетом эффектов турбулентности.
3.2.3. Результаты расчетов в модели идеального газа. 115
Заключение 122
Список литературы
- Методика ускорения решения вспомогательной задачи Римана
- Сеточные методы погашения нефизических эффектов численного решения за сильными малоподвижными ударными волнами
- Результаты расчетов в модели идеального газа
- Начальный этап формирования течения без учета вязких эффектов
Введение к работе
Исследование различных газодинамических течений с помощью моделирования на ЭВМ приобрело в последнее время широкое распространение. Это связано с появлением как простых и надежных численных методов, так и доступных производительных ЭВМ. За последние три десятилетия разработаны новые методы численного моделирования и существенно улучшены старые. Однако, несмотря на бурный рост вычислительной мощности ЭВМ и прогресс в развитии методов компьютерного моделирования газодинамических течений, задача совершенствования методов на настоящий момент остается актуальной. Так, для большинства методов второго порядка сужение зон локализации ударных волн в численном решении порождает побочные флуктуации решения за сильными ударными волнами. Этой проблеме за последние десятилетия посвящено большое количество публикаций [36, 39, 45, 53, 56, 57,74,75,81,87,88,92].
Одним из важных применений численного моделирования является проверка полноты понимания физического явления путем сравнения результатов вычислений с результатами эксперимента. Если экспериментальные данные точны, то расхождение с результатами численного моделирования свидетельствует либо о неучтенных в данной модели эффектах, либо о погрешностях в методике численного расчета. Совершенствование математических моделей и исследование диапазона их применимости -актуальная задача и это направление активно может продвигаться лишь во взаимодействии с экспериментом. В ряде случаев численное моделирование способно заменить эксперимент, однако в большинстве случаев они дополняют друг друга. Для экспериментаторов численное моделирование позволяет лучше понять исследуемое явление, а для
численного моделирования эксперимент позволяет выявлять диапазон применимости модели и изыскивать направления улучшения модели.
Актуальность вышеупомянутых задач совершенствования численных методов и математических моделей в последнее время возросла в связи с широким внедрением в практическое использование пакетов для инженерно-прикладного моделирования, таких как ANSYS, LS-DYNA, STAR-CD. Одним из основных требований, предъявляемых этими программными продуктами к численным методам и математическим моделям, является высокая надежность. Для численных методов это означает надежность вычисления, а для математических моделей - надежность диапазона их применения.
Целями данной работы являются:
создание эффективной методики устранения побочных флуктуации численного решения за фронтом сильной малоподвижной относительно расчетной сетки ударной волной в одномерном и двумерном случае;
численное моделирование и исследование с помощью созданной методики пульсирующих течений перерасширенной струи из короткого осесимметричного сопла;
создание на основе сравнения с результатами эксперимента адекватной численной модели и моделирование течения, возникающего при распространении ударной волны в длинном канале с препятствием.
Исследования проводятся с помощью привлечения нескольких моделей для нестационарных течений сжимаемого газа, а именно: модель идеального газа Эйлера, модель Навье-Стокса вязкого теплопроводного газа и модель Навье-Стокса с учетом эффектов турбулентности. В качестве базового численного метода использовалась хорошо зарекомендовавшая себя W-модификация метода Годунова второго порядка точности, разработанная Васильевым Е.И. [6]. При анализе достоверности расчетов
привлекаются результаты экспериментов, проводимых в рамках совместных исследований с группой ученых университета им. Бен-Гуриона (Израиль, Беер-Шева)
Структура и объем работы. Диссертация состоит из введения, трех глав, заключения и списка литературы; содержит 130 стр., включая 34 стр. с рисунками и 7 стр. списка литературы. В работе 93 библиографических ссылки.
Методика ускорения решения вспомогательной задачи Римана
Базовой моделью для численного исследования газодинамических течений с ударными волнами традиционно считается модель идеального совершенного газа с постоянными теплоємкостями (у= const). В ряде случаев, особенно, когда проводится сравнение с экспериментами, эта модель нуждается в обобщении, которое зависит от особенностей исследуемого течения. Так, например, при высоких температурах за ударными волнами вводят зависимость теплоємкостей от температуры.
Наличие в области течения активных локальных зон возвратного циркуляционного течения приводит к необходимости учета эффектов турбулентности в окрестности этих зон. Для течений, рассматриваемых в данной работе, возникает именно такая ситуация. Учет эффектов турбулентности осуществляется добавлением к базовой модели (А )-модели турбулентности. Промежуточной по сложности численной реализации между этими двумя моделям является модель вязкого газа. Она также используется в данной работе, однако, корректное ее применение возможно лишь при достаточно низких числах Рейнольдса.
Численная реализация для всех трех перечисленных моделей строится на базе W-модификации метода Годунова С.К. [12, 13], развитой Васильевым Е.И. [5, 6], и имеющей второй порядок аппроксимации по пространству и времени. Основные моменты W-метода и его обобщения для используемых моделей изложены во втором параграфе данной главы.
Третий и четвертый параграфы посвящены новым усовершенствованиям для методов годуновского типа. В третьем параграфе предложено ускорение процедуры численного решения задачи о распаде произвольного разрыва (задача Римана). Эта процедура является массовой в методах годуновского типа, так как используется при вычислении потоков через границы расчетных ячеек [13]. Предлагаемая модификация алгоритма приближенного решения задачи Римана позволяет значительно сократить затраты машинного времени на вычисление потоков.
Четвертый параграф посвящен способам устранения нефизических численных флуктуации, возникающих в численных решениях за сильными ударными волнами. При расчетах слабо-нестационарных течений с "медленными" ударными волнами схемы сквозного счета генерируют волнообразные осцилляции за фронтами ударных волн. Эти осцилляции, порождаемые уже схемами первого порядка, становятся более выраженными для схем высокого порядка из-за их низкой диссипации. Хотя такие осцилляции не велики и часто могут быть проигнорированы, они создают помехи, которые весьма нежелательны. Дополнительные побочные численные эффекты возникают при решении двумерных по пространству уравнений Эйлера с применением разностей против потока. Например, известный "карбункул-эффект" и "слоистый шум" за сильными стационарными ударными волнами. В данной работе рассмотрен новый подход к погашению негативных флуктуации, главная роль в котором отводится расчетной сетке.
Для моделирования нестационарных течений совершенного газа в данной работе использовались три модели: модель идеального газа, не учитывающая вязкость и теплопроводность (модель Эйлера), модель вязкого газа (модель Навье-Стокса) и модель вязкого газа с дополнительным учетом (&-е)-модели турбулентности (модель Навье - Стокса + (к-є) ).
Модель Эйлера. Для идеального газа используем дифференциальные уравнения относительно полей плотности, скорости и давления газа в виде законов сохранения массы, импульса и энергии: где плотность среды p, давление /?, полная энергия единицы объема е и вектор скорости для двумерных течений V = (и, v). Отношение теплоємкостей газа принимается равным у= 1.4 . Модель Навье-Стокса. В модели вязкого газа в правой части второго и третьего уравнения системы (1.1) возникают члены с молекулярной вязкостью и теплопроводностью:
Сеточные методы погашения нефизических эффектов численного решения за сильными малоподвижными ударными волнами
При численном моделировании нестационарных газодинамических течений характерно присутствие в решении большого количества разры вов. Наиболее простыми в реализации для численного моделирования этих течений являются методы сквозного счета. Эти методы легко обобщаются на многомерные случаи. Разработаны как методы покоординатного рас щеплення, так и методы «истинно» многомерные. Методы сквозного счета несложно адаптировать к криволинейным и подвижным сеткам. Все это способствовало широкому распространению методов такого типа. В методах сквозного счета производные аппроксимируются через разрывы. Разрыв при этом «размазывается» на некотором отрезке, ширина которого уменьшается с увеличением порядка точности разностной схемы.
Однако такой подход приводит к тому, что схемы сквозного счета, в том числе метод Годунова, в некоторых случаях генерируют нефизические флуктуации численного решения. В частности, флуктуации возникают за фронтом сильных ударных волн, стационарных или медленно движущихся относительно расчетной сетки. В качестве метода подавления флуктуации изначально использовалась искусственная вязкость (линейная или квадратичная). Модификациям, исследованию и применению этого метода посвящено большое количество работ, в частности [26, 28, 33-35, 37, 71, 86, 90]. Обобщение метода на многомерные случаи приведено в обзорной работе [91]. Однако использование метода искусственной вязкости может существенно исказить результаты численного решения [33, 65]. Это вызвано сильным воздействием метода на ширину фронтов слабых ударных волн и скачков, и существенно более слабым - на сильные ударные волны, генерирующие флуктуации. Поэтому результаты, полученные с применением этого подхода, требуют дополнительной проверки, несмотря на появление методов исправления вводимых искусственной вязкостью погрешностей, например [40, 44, 70, 73]. Также для устранения флуктуации достаточно широко используются линейные и нелинейные подавляющие фильтры [23, 30, 38, 65] и специальные добавочные члены [15, 36]. Для устранения этих осцилляции возможно использование гибридных схем, снижающие порядок точности схемы до первого при вычислениях на разрывах. Однако удовлетворительного решения этой задачи в настоящее время не существует.
Максимальная амплитуда этих флуктуации согласно работе [39] не превышает 5% от величины разрыва для схем первого порядка. При повышении порядка точности схемы диссипация схемы снижается, что приводит к увеличению амплитуды флуктуации. Методы сквозного счета воспринимают такие «паразитические» разрывы как элементы численного решения, что в некоторых типах задач может привести к существенному искажению результатов вычислений. Помимо метода искусственной вязкости, к настоящему времени разработаны другие подходы. Так, используются различного рода поправки для численных потоков [39, 53, 75, 81, 82], как непосредственно при их вычислении, так и при решении задачи Римана. Согласно результатам исследований в публикациях [39, 77, 78], причиной возникновения флуктуации численного решения за сильной малоподвижной ударной волной является результат решения задачи Римана, приводящий к дисбалансу между кинетической и внутренней энергиями.
С ростом доступных вычислительных мощностей повсеместно стали решаться многомерные газодинамические задачи на криволинейных сетках. В результатах численного решения некоторых задач исследователями был обнаружен еще один присущий методам сквозного счета эффект нефизического порядка - так называемый «карбункул»-эффект [77]. Этот эффект заключается в характерном искажении фронтов ударных волн и вызывается возникновением во фронте ударной волны слабых флуктуации. Эти флуктуации распространяются вдоль фронта волны. При распространении ударной волны флуктуации могут усиливаться, искажая как фронт ударной волны, так и численное решение за ним. Непосредственной причиной возникновения этого эффекта может являться неоднородность сетки, породившая начальные флуктуации во фронте, или ошибки округления при вычислении координат узлов сетки в случае очень длительных расчетов. Для его погашения, как и в случае с флуктуациями за фронтом ударной волны, также используются поправки [53, 74, 75, 77-79].
В данном параграфе предлагается новый подход к уменьшению степени проявления нежелательных эффектов нефизического характера в результатах численного моделирования нестационарных течений. Основная идея изложенных ниже методов заключается в изменении расчетной сетки во время вычисления потоков через границы ячейки.
Результаты расчетов в модели идеального газа
В данной главе рассматривается течение совершенного газа в перерасширенной сверхзвуковой струе за коротким осесимметричным соплом. Особенностью течения внутри осесимметричного сопла с прямолинейными образующими является наличие так называемого висячего скачка. Этот скачок формируется вблизи критического сечения из-за переразгона потока около стенки сопла. Действительно, при аналитическом контуре сопла поток продолжал бы разгоняться. Однако наличие конического участка является своего рода препятствием. Интенсивность разгона потока уменьшается, а при определенных условиях возможно и торможение потока. Течение торможения вызывает возникновение висячей ударной волны. Вблизи стенки скачок является слабым, но по мере приближения к оси симметрии его интенсивность увеличивается, что может приводить к существенной неоднородности течения внутри сверхзвуковой части сопла [31]. В точке взаимодействия (отражения) висячего скачка с осью симметрии имеет место локальное, но достаточно существенное повышение давления. При этом формируется уходящий отраженный скачок, который ослабевает по мере удаления от оси симметрии.
Для коротких осесимметричных сопел, наиболее привлекательных на практике с технологической точки зрения, положение локального всплеска давления на оси симметрии может выходить за срез сопла и попадать в область струи. В сверхзвуковой струе за соплом также присутствует серия ударных волн (боковые скачки, диск Маха и др.), причем в существенно перерасширенной струе диск Маха может располагаться достаточно близко к срезу сопла. Комбинация выше упомянутых факторов может приводить к взаимодействию диска Маха с приосевой частью висячего скачка и порождать новые сложные волновые конфигурации.
Впервые возможность существования таких конфигураций обнаружена с помощью численного моделирования в работах [7, 42]. Там же отмечен факт пульсирующей нестабильности волновой конфигурации, проявлявшийся в том, что течение в численном моделировании не выходило на стационарный режим. Однако детального изучения закономерностей течения в упомянутых работах не проводилось.
В данной главе в рамках трех моделей: идеального, вязкого и вязкого газа с к-є моделью турбулентности проводится численное исследование закономерностей взаимодействия фронта Маха и висячего скачка для сопла. Геометрические характеристики сопла взяты с точностью до подобия из работ [42, 60, 63] и идентичны использованным в работе [54]. Используемое сопло является осесимметричным JPL - соплом, приведенным в работе [72], с измененным углом наклона сужающейся стенки сопла. Параметры в камере высокого давления взяты из работ [7, 42]. Противодавление в пространстве варьируется в диапазоне от 0.8 до 1.5 атм. Во всем этом диапазоне струя за соплом является перерасширенной (в выходном сечении сопла давление потока лежит в пределах от 0.25 атм (на оси) до 0.5 атм на стенке).
Решение осуществляется численно с помощью описанных в главе 1 методик. Для повышения качества расчетов в области пульсирующих ударных волн применялась вибрирующая сетка (см. 1.4.2). По результатам большой серии расчетов подробно изучается эволюция волновых конфигураций при циклическом изменении противодавления в диапазоне 1.5— 0.8— 1.5 атм. Обнаружены близкие к периодическим пульсации волновых конфигураций и эффект гистерезиса, связанный с их перестройкой. Эффект гистерезиса численного решения по величине давления не нов [41]. Исследуется влияние моделей течения среды (идеальная, вязкая, к-є турбулентность) на величину и характер обнаруженных эффектов.
Постановка задачи и методика численного решения. Рассматривается течение в осесимметричном сопле, камера высокого давления которого (КВД) представляет собой полубесконечный цилиндрический канал, а расширяющаяся часть (камера низкого давления - КНД) имеет выход в неограниченное пространство. Перед стартом газ внутри сопла покоится, причем газы в КВД и КНД имеют существенно различное давление и отделены диафрагмой, которая может быть установлена, например, в критическом сечении (рис. 2.1). При разрыве диафрагмы в некоторый момент времени to начинается так называемый процесс запуска сопла: влево вовнутрь КВД распространяется волна разрежения, а вправо по расширяющейся части движется ударная волна, которая затем выходит в пространство. Через какое-то время внутри сопла формируется стационарное течение, а за соплом формируется сверхзвуковая струя.
Для описания течения будем использовать различные модели, дифференциальные уравнения которых подробно описаны в главе 1. При использовании модели Эйлера на стенках сопла выставляются условия непротекания. В модели Навье-Стокса условие непротекания заменялось условием прилипания, то есть нулевой скоростью на стенке, и дополнительно задавалось условие отсутствия теплообмена газа со стенками (дТ/дп = 0 ).
При использовании к-є модели турбулентности начальное значение кинетической энергии турбулентных пульсаций к принималось равным к = 2-Ю"4 Яоо2, где иг,» - скорость звука в пространстве. Начальное значение скорости диссипации кинетической энергии є подбиралось таким образом, чтобы турбулентная вязкость была в два раза меньше молекулярной вязкости. В этой модели использовалась упрощенная постановка: считалось, что турбулентность внутри сопла заморожена. Тем самым выставление граничных условий для к и є на стенке не требовалось.
Начальный этап формирования течения без учета вязких эффектов
Влияние граничных условий на результаты численного решения. Рассмотренные выше результаты численного решения задачи были получены в постановке уравнений Навье-Стокса с учетом эффекта турбулентности (k-s модель турбулентности). Упомянутые выше вычисления проводились с выставлением на границах канала условий прилипания. Однако ограниченность разрешения расчетной сетки может приводить к неправильному описанию пограничного слоя и его отрыва. Чтобы выяснить роль отрыва пограничного слоя на нижней стенке, были проведены расчеты с граничными условиями скольжения, которые исключают возникновение отрыва на нижней стенке. На рис. 3.8 приведены в альбомной проекции результаты такого численного решения на различные моменты времени. Основным отличием от результатов, приведенных на рис. 3.4-3.6 является отсутствие эффекта отрыва струи J от нижней стенки канала. Действительно, на кадре t\ рис. 3.8 отрыва струи не наблюдается, в то время как в случае выставления условий прилипания отрыв струи происходит (см. кадр t(, рис. 3.5). В результате отсутствия отрыва картина течения существенно меняется. В застойной зоне над струей происходит накопление турбулентности. С течением времени эта зона увеличивается и препятствует образованию нестационарных вихрей, наблюдаемых в предыдущем случае. Вместо этого в струе можно наблюдать (см. кадры t\i) семейство ударных волн. В области течения также можно выделить зоны А и В. Однако область А, в отличие от предыдущей постановки задачи, явно разделяется на две подобласти: струю J и застойную область над ней. Так же, как и в предыдущем случае (см. кадр Ц рис. 3.6), отразившись от правого закрытого торца канала, ударная волна W попадает в область неоднородного течения, о чем свидетельствует вытягивание фронта ударной волны (см. кадр tu).
Показания датчиков при измененном граничном условии отличаются от приведенных на рис. 3.7. Диаграммы давления, построенные по результатам численного решения, приведены на рис. 3.9. Из анализа развития
течения по рис. 3.8 следует, что основные отличия, наблюдаемые в головной части струи, не должны проявляться в показаниях датчика A i и Л з, так как они удалены от этого участка струи. Это подтверждается сравнением эпюр давления на рис. 3.7 и 3.9. Для датчика N2 отличия более значительны.
После прихода фронта волны W на датчики N2 и вплоть до возвращения ударной волны R показания датчиков согласуются с экспериментальными данными вполне удовлетворительно, а для верхней стенки даже лучше, чем приведенные на рис. 3.7 результаты вычислений с прежними граничными условиями. Это следует из падения давления в показаниях датчиков Ni после прохождения фронта W, как и в экспериментальных данных, хотя и более резкое, в то время как в результатах вычислений на рис. 3.7 снижения давления практически не происходит. Для показаний нижнего датчика N2 существенным отличием от приведенных на рис. 3.7 графиков являются флуктуации давления с большой амплитудой непосредственно при прохождении фронта волны R. Эти всплески объясняются общей картиной течения в струе. В силу отсутствия отрыва струи и вызванной им турбулентности в струе возникают малоподвижные скачки, разделенные областями с очень сильными градиентами давления. Давление повышается в скачке, находящемся в струе, и затем достаточно резко падает вдоль «бочки» вплоть до следующего скачка, которым замыкается «бочка». Струя в целом нестабильна и упомянутые скачки испытывают небольшие колебательные перемещения вдоль оси х. Если позиция датчика попадает в область смещения скачка, датчик регистрирует пилообразные всплески давления с течением времени. В расчетах позиция датчика N2 не попадает в такую область. Однако при приходе ударной волны R волновая конфигурация в струе деформируется. Датчик регистрирует колеблющееся поведение одного из скачков сжатия внутри струи, что видно на графике рис. 3.9 при t = 30 единиц безразмерного времени.
Сравнение этих расчетов с экспериментальными эпюрами (см. рис. 3.7) для датчика Л на оси симметрии, а также для других сечений внутри струи, которые здесь не приводятся, свидетельствует об отсутствии в реальном течении такой ярко выраженной последовательности скачков подобной той, что наблюдается на рис. 3.8, кадры tyU.
В тоже время, сравнивая показания датчика Л на стенке с экспериментальными данными, видим, что при условиях скольжения согласование с экспериментом лучше. Это, по-видимому, свидетельствует о том, что в реальном течении эффект отрыва на нижней стенке слабее и оттесняемая отрывом струя гораздо меньше воздействует на верхнюю стенку, где расположен датчик Л .
Таким образом, сравнение расчетных эпюр давления с двумя видами граничных условий с экспериментальными данными в трех сечениях канала не позволило однозначно отдать предпочтение тому или иному виду граничных условий. При условии прилипания лучше согласуются с экспериментом датчики на нижней стенке канала, а при условии скольжения - на верхней стенке.
На взгляд автора более правдоподобны условия прилипания, хотя для более определенного ответа нужна экспериментальная фотосъемка картины течения. На рис. 3.10 в альбомной проекции представлены результаты расчетов течения в модели Эйлера в виде изолиний плотности на различные моменты времени. На серии последовательных кадров изображена часть канала от препятствия до торца (собственно правый торец канала на рисунках не показан и находится за правой границей).
На кадре t\ изображена ударная волна W, прошедшая через препятствие. Стрелкой показано направление распространения ударной волны по потоку. Отраженная от препятствия и ушедшая налево против потока ударная волна на рисунке отсутствует.
