Содержание к диссертации
Введение
1. Современные подходы к численному моделированию трансзвуковых течений газа 12
1.1. Исходные положения 12
1.1.1. Уравнения одномерной газовой динамики 12
1.1.2. Общий обзор известных численных схем 13
1.2. Противопоточные схемы 15
1.2.1. Методы расщепления векторов потоков 15
1.2.2. Схема Годунова 19
1.2.3. Схема Ошера 20
1.2.4. Схема Роу 22
1.3. Схемы с искусственной диссипацией. CUSP схема 25
1.3.1. Концепция искусственной диссипации 25
1.3.2. Базовые положения CUSP схемы 26
1.3.3. CUSP схема, обеспечивающая постоянство полной энтальпии в стационарном потоке 32
1.4. Повышение порядка точности схем 33
1.4.1. MUSCL подход и метод экстраполяции потоков 33
1.4.2. Монотонные, TVD и LED схемы 36
1.4.3. JST схема 38
1.4.4. Введение ограничителей. SLIP схема 41
1.4.5. Мягкий ограничитель 43
1.4.6. Обобщение на систему уравнений Эйлера 45
1.5. Методы регуляризации уравнений газовой динамики 46
1.5.1. Предварительные замечания 46
1.5.2. Регуляризация несжимаемых уравнений Эйлера 48
1.5.3. Регуляризация сжимаемых уравнений Эйлера 50
1.5.4. Определение параметра сжимаемости 52
2. Математическая модель и численный метод 55
2.1. Математическая модель 56
2.1.1. Определяющие уравнения 56
2.1.2. Введение обобщенных координат 58
2.1.3. Геометрические характеристики и индексация конечных объемов 60
2.2. Моделирование турбулентности 62
2.2.1. Концепция турбулентной вязкости 62
2.2.2. Модель Спаларта-Аллмараса 63
2.2.3. Высокорейнольдсовая k-є модель 64
2.3. Дискретизация стационарного оператора 66
2.3.1. Дискретизация конвективных потоков 66
2.3.2. Дискретизация вязких и диффузионных членов 68
2.4. Неявная схема и дискретизация стабилизирующего оператора 69
2.4.1. Регуляризация уравнений по методу масштабирования сжимаемости 69
2.4.2. Построение неявной схемы 71
2.4.3. Метод приближенной факторизации 72
2.4.4. Дискретизация конвективной части стабилизирующего оператора по схеме, вводящей скалярную диссипацию 73
2.4.5. Дискретизация конвективной части стабилизирующего оператора по схеме расщепления матрицы коэффициентов 74
2.4.6. Дискретизация диффузионной составляющей стабилизирующего оператора 75
2.4.7. Решение результирующей системы 75
2.4.8. Ускорение сходимости за счёт введения локального шага установления повремени 76
2.5. Использование блочно-структурированных сеток 77
3. Методические и тестовые расчеты 80
3.1. Квазиодномерное невязкое течение в сопле Лаваля 80
3.2. Невязкое двумерное обтекание крылового профиля NACA-0012 81
3.3. Ламинарный сверхзвуковой пограничный слой на продольно обтекаемой пластине 82
3.4. Турбулентный пограничный слой на продольно обтекаемой пластине 85
3.5. Двумерное турбулентное обтекание решетки турбинных лопаток 88
4. Пространственные турбулентные течения в диффузорах и межлопаточных каналах турбомашин 119
4.1. Влияние геометрии поворотного участка и входных условий на течение газа в осерадиальных диффузорах 119
4.2. Течение в осерадиалыюм диффузоре при наличии локального тангенциального вдува 122
4.2.1. Предварительные замечания 122
4.2.2. Описание экспериментов ЦКТИ 124
4.2.3. Численное моделирование 126
4.3. Расчет трехмерного течения и потерь давления в выхлопном патрубке мощной паровой турбины 128
4.3.1. Предварительные замечания 128
4.3.2. Описание экспериментов ЛМЗ 129
4.3.3. Трехмерное численное моделирование 131
4.4. Численное моделирование трехмерного течения в трансзвуковой турбинной решетке 133
4.4.1. Предварительные замечания 133
4.4.2. Описание эксперимента 134
4.4.3. Трехмерное численное моделирование 135
Заключение 164
Литература 167
Приложения 177
- CUSP схема, обеспечивающая постоянство полной энтальпии в стационарном потоке
- Методы регуляризации уравнений газовой динамики
- Геометрические характеристики и индексация конечных объемов
- Дискретизация конвективной части стабилизирующего оператора по схеме расщепления матрицы коэффициентов
Введение к работе
Исследование и инженерный анализ сложных течений вязкого газа методами численного моделирования приобрели в последнее время широкое распространение. Это связано как с разработкой эффективных численных методов, так и с развитием высокопроизводительных ЭВМ.
Одним из важных аспектов численного моделирования является сравнение результатов вычислений с данными эксперимента. Если экспериментальные данные точны, то расхождение с результатами численного моделирования свидетельствует о наличии эффектов, неучтенных или неправильно воспроизведенных в используемой математической модели или же о недостатках/погрешностях численного метода, требующих устранения. Совершенствование численных методов и исследование диапазона их применимости - актуальная задача, успешное решение которой может происходить только во взаимодействии с экспериментом. В ряде случаев численное моделирование способно заменить эксперимент, однако в большинстве случаев они дополняют друг друга. Для исследователей, производящих измерения на опытных установках, численное моделирование позволяет лучше понять изучаемое явление, а для численного моделирования эксперимент помогает определить диапазон применимости модели и находить пути ее улучшения.
Применительно к проточным частям турбомашин проведение измерений связано с определенными трудностями, в частности, с высокими требованиями к разрешающей способности аппаратуры или большими размерами исследуемых моделей, что обуславливает высокую стоимость проведения экспериментов. В связи с этим численное моделирование имеет особо большое значение для практического применения, поскольку позволяет снизить объем дорогостоящих экспериментов, а также предсказать в процессе проектирования новых машин характеристики течения для широкого диапазона режимных параметров.
В проточных частях турбомашин весьма часто реализуются трансзвуковые течения. Особую важность при численном моделировании таких течений приобретает способ аппроксимации конвективных потоков в уравнениях Навье-Стокса. Одно из главных требований, предъявляемых к методам аппроксимации конвективных слагаемых, заключается в обеспечении возможности достаточно точно разрешать ударные волны на небольшом числе внутренних точек. Другим, не менее важным требованием, является отсутствие осцилляции поля течения в окрестности ударных волн.
Указанными свойствами обладают получившие большую популярность за последние десятилетия методы, основанные на использовании решения задачи Римана о распаде разрыва. Среди них особенное распространение получил метод Роу [Roe 1981, 1985, 1986], основанный на решении линеаризованной задачи Римана. Однако при проведении сложных трехмерных расчетов матричные операции, необходимые для вычисления конвективных потоков в схеме Роу становятся весьма затратными. В связи с этим в последние годы появились более экономичные схемы, использующие идеи, заложенные в схеме Роу, и не уступающие ей по качеству получаемых решений. Среди них - так называемая CUSP (Convective Upwind Split Pressure) схема, которую, по существу, можно трактовать, как схему с матричной искусственной диссипацией [Jameson 1995а; Swanson et al. 1997]. В то же время при практической реализации CUSP схемы матричных операций при вычислении конвективного потока удается избежать.
Для получения стационарных решений уравнений Навье-Стокса часто применяется метод установления по псевдовремени в дельта-форме. Известно, что эффективность этого метода снижается в случае существенно дозвуковых течений. Это связано с «жесткостью», которую приобретают уравнения динамики сжимаемого газа при существенно дозвуковых скоростях потока вследствие сильного различия характерных времен процессов конвективного переноса и распространения акустических возмущений. Для построения метода, одинаково хорошо работающего при любых числах Маха, необходимо применять специальные процедуры регуляризации уравнений. Одним из таких методов является метод масштабирования сжимаемости [Стрелец, Шур 1988].
Исходя из изложенных соображений, определены следующие цели и задачи работы:
разработка и реализация трехмерной версии CUSP схемы в сочетании с методом масштабирования сжимаемости для расчетов трансзвуковых течений на основе базового программного комплекса SINF, развиваемого сотрудниками кафедры гидроаэродинамики СПбГПУ;
проведение ряда методических расчетов с целью определения возможностей реализованного алгоритма по разрешению ударных волн, а также ламинарных и турбулентных пограничных слоев;
апробация разработанной расчетной схемы на примере двумерного трансзвукового турбулентного обтекания решетки турбинных лопаток;
проведение параметрических расчетов турбулентного дозвукового течения в осерадиальных диффузорах с целью выяснения влияния геометрии поворотного участка и входных условий на течение газа и коэффициент полных потерь;
проведение численного моделирования течения в осерадиальном диффузоре при наличии локального тангенциального вдува;
расчет трехмерного трансзвукового турбулентного течения газа и потерь давления в проточной части выхлопного патрубка мощной паровой турбины;
численное моделирование трехмерного течения в трансзвуковой турбинной решетке.
Первая глава диссертации посвящена обзору методов численного моделирования трансзвуковых течений. Проводится классификация численных методов для моделирования такого рода течений, приводится детальное описание некоторых расчетных схем для решения одномерных уравнений Эйлера. Рассматривается вопрос повышения порядка аппроксимации и связанная с ним проблема возникновения осцилляции в численном решении в окрестности ударных волн. Также приводится классификация методов регуляризации уравнений с целью ускорения сходимости к стационарному решению. В рамках введенных классификаций излагаются одномерная версия CUSP схемы и метод масштабирования сжимаемости, использующиеся при расчетах в настоящей работе.
Во второй главе дается изложение математической модели и численного метода, на основе которого были проведены расчеты. Проводится обобщение CUSP схемы на случай трехмерных уравнений Навье-Стокса применительно к использованному для расчетов программному "конечно-объемному" комплексу SINF, второго порядка точности. Дается описание неявной итерационной схемы и дискретизация стабилизирующего оператора, рассматриваются некоторые способы ускорения сходимости к сошедшемуся решению.
В третье главе приводятся результаты методических и тестовых расчетов. С целью определения возможностей CUSP схемы по разрешению ударных волн рассматриваются две задачи о течении невязкого газа. Первая из них формулируется для нерасчетного квазиодномерного режима течения в сопле Лаваля. Расчетные распределения газодинамических параметров сравниваются с точным (аналитическим) решением. Вторая задача, рассматриваемая в невязкой постановке, -трансзвуковое обтекание крылового профиля NACA-0012. Результаты расчетов по CUSP схеме сопоставляются с высокоточными результатами численного моделирования, доступными в литературе. Далее в третьей главе, с целью определения точности CUSP схемы при разрешении пограничных слоев, рассматриваются задачи о течении в пограничном слое на продольно обтекаемой теплоизолированной пластине. Расчетные значения коэффициентов пристенного трения и восстановления полной температуры, полученные при моделировании ламинарного сверхзвукового пограничного слоя, сравниваются с результатами, следующими из теории пограничного слоя. Моделирование турбулентного
пограничного слоя проводится с использованием модели турбулентности Спалларта-Аллмараса. Расчетные значения коэффициентов трения и восстановления, а также формпараметра пограничного слоя сравниваются с экспериментальными данными. Наконец, представляются результаты численного моделирования двумерного течения в решетке турбинных направляющих лопаток. Расчеты проводятся в широком диапазоне выходного числа Маха как для ламинарного, так и для турбулентного режимов течения. Расчетные значения статического давления по обводу лопатки, коэффициент потерь энергии и угол выхода потока сравниваются с опытными данными.
В четвертой главе излагаются постановка задачи и результаты расчетов пространственных турбулентных течений в диффузорах и межлопаточных каналах турбомашин. Анализируется влияние геометрии поворотного участка и входных условий, а также локального тангенциального вдува на течение газа в осерадиальных диффузорах. Расчетные значения полных потерь и потерь давления сравниваются с экспериментами, проведенными в ЦКТИ и "Skoda Energo". Представляются данные расчетов трехмерного течения и потерь давления в выхлопном патрубке мощной паровой турбины. Расчетные значения общих потерь сравниваются с экспериментальными значениями, полученными на испытательном стенде ЛМЗ. На основе детального анализа сложного трансзвукового течения в выхлопном патрубке предлагается оптимизировать форму и положение дефлектора с целью уменьшения интенсивности ударных волн и снижения потерь давления. Описываются расчетные данные, полученные для трехмерного трансзвукового турбулентного течения в плоской решетке турбинных профилей. Результаты выполненных расчетов сравниваются с экспериментальными данными и с результатами аналогичных расчетов, проведенных с использованием коммерческого пакета Fluent. Проводится анализ вторичных течений, выясняется влияние вторичных вихревых структур на потери полного давления.
В заключении представлены основные результаты и выводы, полученные в работе.
1. СОВРЕМЕННЫЕ ПОДХОДЫ К ЧИСЛЕННОМУ
МОДЕЛИРОВАНИЮ ТРАНСЗВУКОВЫХ ТЕЧЕНИЙ ГАЗА
1.1. ИСХОДНЫЕ ПОЛОЖЕНИЯ
1.1.1. Уравнения одномерной газовой динамики
В природе и технике подавляющее число течений происходит при больших числах Рейнольдса, вдали от стенок в них преобладают конвективные эффекты, а потому по своим свойствам эти течения близки невязким течениям, описываемым системой уравнений Эйлера. При численном интегрировании уравнений Навье-Стокса особую важность имеет дискретизация конвективных производных, особенно при рассмотрении сверхзвуковых течений и течений с ударными волнами. Дискретизация диффузионной части уравнений не представляет особой сложности, по этой причине в рамках настоящего обзора мы рассмотрим только вопросы численного интегрирования уравнений Эйлера. Обобщение рассматриваемых методов на уравнения Навье-Стокса производится непосредственно и не представляет сложности.
В рамках настоящего обзора будем рассматривать только одномерную систему уравнений Эйлера, которую запишем в консервативной форме
^+Ад*)=о, (l.i.i.i)
где w - вектор консервативных переменных, f{w) - вектор потоков:
w = (p ри рЕ)т, (1.1.1.2)
7 = \ри ри2 + р риН) , (1.1.1.3)
E = cvT + \u2,H = cpT + \u2. (1.1.1.4)
Здесь и - скорость, р - плотность, р- давление, Е - полная энергия, Н - полная
энтальпия, Т - температура, теплоемкости ср, cv и их отношение cp/cv = y
принимаются постоянными. Система (1.1.1.1) замыкается уравнением состояния, которое здесь и далее представляет собой уравнение состояния совершенного газа
Р- = {с -cv)r. (1.1.1.5)
13 1.1.2. Общий обзор известных численных схем
Методы решения системы уравнений Эйлера можно разделить на две группы [Hirsch 1990]. Первая объединяет противопоточные схемы, в которых для аппроксимации пространственных производных используются односторонние разности, расписываемые с учетом направления распространения малых возмущений, тем самым дискретизация производится с учетом физических свойств системы. К этой группе относятся, например, схемы Стегера-Уорминга [Steger, Warming 1981] и Ван-Лира [Van Leer 1982], основанные на расщеплении векторов потоков. Особо важным классом в группе противопоточных схем являются схемы типа Годунова, базирующиеся на решении задачи Римана о распаде разрыва. Учет физических свойств течения в схемах типа Годунова носит более глубокий характер, поскольку дискретизация производится с привлечением частных аналитических решений уравнений Эйлера. Однако решение задачи Римана сопряжено с большими вычислительными издержками, поскольку на каждом шаге по времени требует решения нелинейных алгебраических уравнений для каждой грани расчетной сетки. В связи с этим получили широкое распространение методы приближенного решения задачи Римана, такие как метод Ошера [Osher 1981] и метод Роу [Roe 1981].
Вторая группа - это центрально-разностные методы, в которые для повышения устойчивости тем или иным способом добавляется численная диссипация. Природа численной диссипации может быть различна. Так в схемах, в которых временная и пространственная дискретизации проводятся одновременно, численная диссипация возникает естественный образом за счет дискретизации производной по времени. Примерами таких схем являются схема Лакса-Фридрихса [Lax 1954] и схемы типа Лакса-Вендроффа [Lax, Wendroff 1960], причем в последних величина численной диссипации оказывается недостаточной, что приводит к возникновению осцилляции в областях больших градиентов. Для их подавления в схему вводится специально сконструированный дополнительный диссипативный поток - искусственная диссипация, действие которой аналогично действию физической вязкости. Одним из недостатков схем типа Лакса-Вендроффа является зависимость решения от выбранного шага интегрирования по времени, что является следствием одновременной дискретизации временных и пространственных производных. Отмеченного недостатка лишены схемы, в которых дискретизации пространственных и временных производных производятся независимо. К этому типу относятся, например, схема Бима-Уорминга [Головачев 1997; Флетчер 1991] и JST схема [Jameson et. al 1981]. Здесь, однако, естественный механизм возникновения численной вязкости отсутствует, что приводит к необходимости конструирования искусственных диссипативных слагаемых.
Заметим, что противопоточные методы, например, такие как методы Роу или Ошера, также обладают численной диссипацией, и именно "противопоточность" является естественным механизмом ее возникновения. Схема Роу видится особенно привлекательной, так как она допускает существование идеальных дискретных скачков уплотнений с единственной внутренней точкой. Если рассматривать схему Роу с точки зрения схем, основанных на концепции численной диссипации, то коэффициент численной диффузии для схемы Роу оказывается матричной величиной. Последнее обстоятельство особенно сказывается при проведении трехмерных расчетов, когда применение схемы Роу становится чрезвычайно трудоемким с вычислительной точки зрения.
В последние годы в литературе [Jameson 1995а] появились более экономичные схемы, основанные на схеме Роу и сохраняющие свойство хорошего разрешения ударных волн. Это, в частности, так называемая, CUSP (Convective Upwind Split Pressure) схема и схема со скалярной искусственной диссипацией. CUSP схему, как и схему Роу, можно трактовать, как схему с матричной численной диссипацией, однако при ее практической реализации удается избежать матричных операций при вычислении численной диссипации. В то же время CUSP схема допускает существование ударных волн с единственной внутренней точкой. Схема со скалярной диссипацией является еще более экономичной, чем CUSP схема, однако она не допускает существование идеальных ударных структур.
Другим, не менее важным, чем вопрос о выборе схемы дискретизации, является вопрос повышения порядка точности схем. Известно, что схемы порядка точности выше первого порождают нефизические осцилляции в окрестности ударных волн. Для их устранения на этапе дискретизации приходится вводить дополнительную нелинейность, за счет использования, так называемых, ограничителей, которые локально понижают порядок точности схемы до первого. Недостатком ограничителей является то обстоятельство, что они действуют не только в окрестности скачка, но и на гладких экстремумах, уменьшая локальную точность схемы до первого порядка. Мягкие ограничители позволяют повысить порядок точности в областях гладких экстремумов, путем введения порога в ограничитель, однако в этом случае появляется дополнительный параметр, подбор которого индивидуален для конкретной задачи. Высокого порядка точности иа гладких экстремумах удается достичь в получивших большую популярность в последние годы ENO и WENO схемах [Harten et al. 1987; Shu 1997], однако их использование также связано со значительными вычислительными издержками.
Решение уравнений Эйлера в задачах с малыми числами Маха осложняется значительным возрастанием времени необходимого для достижения установившегося решения, что связано с жесткостью, которую приобретают уравнения при малых
скоростях потока. Для построения метода, одинаково хорошо работающего для любых чисел Маха, необходимо применять специальные процедуры регуляризации уравнений. Суть методов регуляризации состоит в построении системы дифференциальных уравнений, которая, с одной стороны, имеет то же стационарное решение, что и система уравнений Эйлера, а с другой - позволяет устранить жесткость исходной нестационарной дифференциальной задачи при малых значениях чисел Маха.
1.2. ПРОТИВОПОТОЧНЫЕ СХЕМЫ
1.2.1. Методы расщепления векторов потоков
При построении схем, основанных на расщеплении потоков, учитываются направления распространения малых возмущений, переносимых различными типами волн, т.е. характеристические свойства исходных уравнений, благодаря чему удается достичь высокого качества получаемых решений.
Рассмотрим в начале, так называемую, схему расщепления матрицы коэффициентов, наиболее полно сформулированную в работе [Chakravarthy et al. 1980]. Эта схема базируется на неконсервативной форме записи уравнений Эйлера
dw ,dw . — + А— = 0, dt дх
где А - это матрица Якоби, которая имеет вид
(1.2.1.1)
0 (7-І)
(у-ъу
А =
(1.2.1.2)
(Г-З)и yuE уЕ Ъ{у - \)u"
іу-іУ-
\
Р Р 2
Если система (1.2.1.1) стоит только из одного скалярного уравнения
(1.2.1.3)
dw dw Л — + а— = 0, dt dx
то для его дискретизации в точке с индексом j можно применить схему с аппроксимацией производных противопоточными разностями. В том случае, когда значение а имеет произвольный знак, явная противопоточная схема первого порядка точности запишется в виде
wf-wn
(1.2.1.4)
^~ +—(a+A-'-a-A+-Yw"=0,
Ахк h j >
где через
A~Wj = Wj - Wj_x, A+Wj = wJ+l -Wj (1.2.1.5)
обозначены разности первого порядка "назад" и "вперед" соответственно, и
**=(ау±|ву|)/2. (1.2.1.6)
Заметим, что при любых знаках а справедливо равенство
a*+a]=aj. (1.2.1.7)
Обобщение схемы (1.2.1.4), (1.2.1.5) на случай системы нескольких уравнений дается в виде
wTl-w1
1 -w- 1 / \
^ J- +—[A+A--+A-A+-)w"=0, (1.2.1.8)
At Axv J J
где введены обозначения
A±=RA±R~\ А± =diag(A±k). (1.2.1.9)
4=(Л±|Л|)/2, (1.2.1.10)
и Як, к = 1,2,3 - собственные числа матрицы A; R - матрица, составленная из
правых собственных векторов матрицы А как из столбцов.
Важно отметить то обстоятельство, что составляющие расщепленной матрицы А = А++А~ имеют неотрицательные и неположительные собственные числа (Я+ и Л~) соответственно, обеспечивая тем самым учет характеристических свойств системы при аппроксимации пространственных производных. Последнее обеспечивает хорошее качество решений при расчете течений без газодинамических разрывов. Однако схема становится неприменимой для расчета течений с ударными волнами, поскольку базируется на неконсервативной форме записи исходной системы уравнений и не обеспечивает выполнения законов сохранения на разрывах, что даже на сколь угодно мелких сетках приводит к неверным значениям амплитуды изменения параметров потока на скачках. Отмеченный недостаток схемы расщепления матрицы коэффициентов устраняется в близких к ней по идеям схемам расщепления векторов потоков.
Рассмотрим сперва первый вариант этого метода, предложенный Стегером и Уормингом [Steger, Warming 1981]. Этот метод существенно базируется на свойствах собственно уравнений Эйлера, в отличие от метода расщепления матрицы коэффициентов, сформулированного для произвольной гиперболической системы уравнений.
Будем исходить из уравнений Эйлера записанных в консервативной форме (1.1.1.1) и рассматривать вектор потоков как функцию консервативных переменных. Принципиально важно, что вектор газодинамических потоков является однородной
функцией первой степени от компонент вектора консервативных переменных W, то есть что
f = Aw. (1.2.1.11)
Справедливость этого равенства можно проверить прямым перемножением матрицы А на вектор w.
Соотношение (1.2.1.11) позволяет представить вектор потоков / в виде суммы
составляющих /+ и /~, отвечающих положительным и отрицательным собственным числам матрицы А
f = Aw = (A++A~)w = f++f-,me (1.2.1.12)
f±=A±w = RA±R~lw. (1.2.1.13)
Заметим, что Якобианы расщепленных потоков d^/dw, вообще говоря, отличны от
расщепленных Якобианов А±; отличны между собой и их собственные числа. Аналогично схеме (1.2.1.8), можно построить противопоточную схему
л У+т-А-/++Л+Г)>0. (1.2.1.14)
At Ах
Необходимые для построения векторов /* собственные числа матрицы А
имеют вид
Л1 = Л = и, A2,i ~^ =и±с, (1.2.1.15)
где с - Jcp(y-l)T - скорость звука.
Важно отметить тот факт, что в сверхзвуковом течении при и > 0 все собственные числа матрицы А положительны, Ґк = Хк, ?Гк = 0 и, следовательно,
/+=/, /~=0, т.е. в схеме (1.2.1.14) присутствуют только разности назад. Аналогично для сверхзвукового потока, движущегося в противоположном направлении, схема (1.2.1.14) будет содержать только разности вперед, то есть в обоих случаях возмущения по разностной сетке передаются только вниз по потоку. Если течение дозвуковое |w| < с, то схема (1.2.1.14) содержит как разности вперед, так
и разности назад.
Существенный недостаток рассмотренной схемы заключается в недифференцируемости расщепленных потоков в точках, где собственные числа
меняют знак, т. е. в точках и = 0 и и = ±с. Это связано с тем, что потоки j содержат
в качестве множителей величины Л*, определяемые по формулам (1.2.1.10), которые,
в свою очередь являются недифференцируемыми функциями Лк.
Недифференцируемость расщепленных потоков /* в окрестностях точек торможения
и звуковых точек нередко приводит к невозможности получить сошедшееся решение из-за возникновения осцилляции в окрестностях этих точек.
Простой способ устранения осцилляции был предложен в работе [Buning, Steger 1982], путем изменения способа вычисления величин А±, используемых при построении расщепленных потоков /*
3і -Лшсі —
А±^ІА2+є2
(1.2.1.16)
где є - некоторая малая величина. Неравенства А+той > О, Amod < 0, выполнение
которых необходимо для построения схемы расщепления векторов потоков, справедливы по-прежнему и для такого способа расщепления. В то же время
использование (1.2.1.16) приводит к тому, что оба потока /*, не обращаются в ноль ни при каких значениях Ак, в том числе и при \и\ > с, что равносильно появлению
механизма распространения возмущений вверх по потоку в сверхзвуковых областях течения. Кроме того, выбор таких значений параметра є, которые обеспечивали бы одновременно сходимость к стационарному решению и отсутствие пространственных осцилляции решения в областях больших градиентов, это нетривиальная задача, особенно при расчете течений с различными числами Маха.
Указанного недостатка лишен метод расщепления потоков Ван Лира [Van Leer 1982], в котором введено отличное от (1.2.1.10) и (1.2.1.16) расщепление собственных
значений. В этом методе при конструировании потоков /+ и /~, кроме выполнения соотношения (1.2.1.12), требуется выполнение дополнительных условий. В частности требуется, чтобы потоки /* и Якобианы df± /dw были непрерывными функциями числа Маха и являлись полиномами, по возможности, наиболее низкой степени. Кроме того, собственные значения df+/dw должны быть неотрицательны, а
собственные значения df~/dw - неположительны, при этом в дозвуковой области
течения одно из собственных значений Якобианов должно обращаться в нуль. Расщепленные потоки схемы Ван Лира имеют следующий вид
f^±Pf(M±\y
УК 2 ,
(1.2.1.17)
f-x
1±—-м 2 )
а соответствующие собственные числа определяются формулами
y + \
lf^(M + l):
/ + 1
3-M +
(M-l)2
ґ2=^(м+іу
4-3^+^-1(1+^
19 (1.2.1.18)
(1.2.1.19) (1.2.1.20)
Л~(М)=-Л;(-М), / = 1,3. (1.2.1.21)
Очевидно, что расщепленные собственные значения Ван Лира отличаются от расщепленных собственных значений Стегера-Уорминга. Равенство нулю одного из
собственных значений Якобиана df*/дМ при \м\ < 1 позволяет разрешать ударные
волны не более чем на двух ячейках расчетной сетки [Van Leer 1982].
1.2.2. Схема Годунова
Основная идея схемы Годунова [Годунов 1959] состоит в использовании точных решений уравнений газовой динамики при аппроксимации пространственных производных. Для гиперболической системы уравнений Эйлера такие решения распадаются на совокупности независимых и сравнительно просто рассчитываемых деталей - "распадов разрывов".
При реализации численного метода с первым порядком точности по пространству, в каждый момент времени значение газодинамических функций в пределах расчетной ячейки считается постоянным и равным значению в ее центре. Иными словами, во всей расчетной области искомые функции представляют собой кусочно-постоянные распределения, терпящие разрывы на гранях расчетной сетки (которые отмечаются полуцелым индексом). Таким образом, на каждой грани расчетной сетки реализуются условия задачи о распаде произвольного газодинамического разрыва, имеющей аналитическое решение [Годунов, Забродин 1976; Hirsch 1990]. Распад разрыва сопровождается возникновением ударных волн, центрированных волн разрежения и контактных разрывов. При этом, в зависимости от начальных условий, могут возникнуть различные конфигурации течения, параметры которого вычисляются по разным формулам, имеющим обобщенный вид
w = iH&w,,wy+1), (1.2.2.1)
*7+1/2
где автомодельная переменная =
х-х,
t-t'
Из этих формул видно, что значения
консервативных переменных в точке xJ+U2 будут сохраняться равными некоторым постоянным, которые мы обозначим Wjl'U2, до тех пор, пока в эту точку не придут
волны, образовавшиеся при распаде разрывов на соседних гранях. После этого рассмотренную схему течения нужно пересматривать заново. В случае неравномерной сетки моменты перестройки будут различны для разных узлов. В этих условиях процесс распада разрыва рассматривается в течение не слишком большого промежутка времени т
г < minДху /тах(я), (1.2.2.2)
такого, чтобы значения w\+\n, на всех гранях считать неменяющимися.
Распределение газодинамических величин по координате х на момент времени tn+l, построенное с помощью алгоритма расчета распада разрывов на всех гранях xj+i/2 оказывается весьма сложным. В оригинальном методе Годунова оно заменяется кусочно-постоянным распределением, которое получается, если внутри каждой ячейки плотности массы р, импульса ри и энергии рЕ полагаются равными их
і У+1/2 **J ,-,,2
средним значениям по ячейке. Тогда после осреднения величин wjl{/2 в пределах ячейки
/Я(л)(х,Ги+фс. (1.2.2.3)
и интегрирования системы (1.1.1.1) по методу конечного объема, для нахождения этого нового кусочно-постоянного распределения получается соотношение
{wf -*;)дгу + М"[/;+ш -ТНіг\= 0, (1.2.2.4)
где / - численный поток, для схемы Годунова имеющий вид
/#=/И)- о-2-2-5)
Как следствие процедуры осреднения (1.2.2.3), некоторые детали точного решения задачи Римана оказываются потерянными. Так как точное решение задачи Римана требует разрешения нелинейного алгебраического уравнения, что является довольно трудоемкой процедурой, рядом авторов были предложены приближенные методы ее решения. Наиболее интересные аппроксимации решений задачи Римана были разработаны Ошером [Osher 1981; Osher and Solomon 1982] и Poy [Roe 1981].
1.2.3. Схема Ошера
Схема Ошера решения системы уравнений Эйлера основывается на противопоточной схеме, разработанной Энгкуистом и Ошером [Engquist, Osher 1980, 1981] для закона сохранения скалярной величины v
^-Д-0. (1.2.3.1)
dt дх
Численный поток схемы Энгкуиста-Ошера имеет вид
где a(v) = dfjdv - характеристическая скорость распространения возмущений по газу.
Главное достоинство схемы (1.2.3.2) по сравнению с простейшей противопоточной схемой Мурмана-Коле [Murman, Cole 1971]
//+1/2 = ±(/}+i +fj)-i"gn{j+uiyfj+i ~fj) (1-2.3.3)
- простейшей схемой интегрирования скалярного уравнения переноса, заключается в том, что схема (1.2.3.2) удовлетворяет энтропийному условию [Lax 1973], т.е. не допускает возникновения в решении ударных волн разрежения, тогда как напротив, схема Мурмана-Коле может приводить к их возникновению [Hirsch 1990].
Схему Энгкуиста-Ошера (1.2.3.2) можно рассматривать как схему Годунова с приближенным решением задачи Римана. В работе [Van Leer 1984] показано, что единственное различие при решении задачи
Wi+2/З
Л=и-с
Щ+1
Рисунок 1.1. Интегрирование в фазовом пространстве по схеме Ошера для одномерной Римана в двух схемах заключается в том, что в системы уравнений Эйлера схеме Энгкуиста-Ошера скачок уплотнения
заменяется центрированной волной сжатия, что является результатом интегрирования в соотношении (1.2.3.2). Как следствие, схема не допускает возникновения ударных волн разрежения, в то же время разрешая ударные волны не более чем на двух ячейках.
Обобщение схемы Энгкуиста-Ошера на случай решения системы уравнений Эйлера проведено в работах [Osher 1981; Osher, Solomon 1982] и основывается на представлении (1.2.3.2), в котором одномерное интегрирование заменяется интегрированием по фазовому пространству консервативных переменных w, а функция |a(v)j - на абсолютное значение Якобиана \a(wJ . Таким образом, численный поток схемы Ошера имеет вид
/j+i/2 =i(fj +/y+i)"±f 14*)И- (1-2.3.4)
Путь интегрирования в фазовом пространстве от Wj до ivJ+l раскладывается на
совокупности характеристик отвечающих различным собственным значениям системы (1.1.1.1). Очевидно, что подобное разложение не является единственным, поскольку из одной точки фазового пространства может быть выпущено несколько характеристик, отвечающих различным собственным числам. Интегрирование,
изначально предложенное Ошером, производится вдоль характеристик в последовательности, отвечающей убыванию собственных чисел. Это иллюстрируется на рисунке 1.1, где интегрирование сперва производится по характеристике отвечающей собственному числу и + с, затем и и, наконец, и-с. Интегрирование можно проводить и по пути, соответствующему последовательности возрастающих собственных чисел [Hemker, Spekreijse 1986].
Промежуточные точки iVy+i/3 и wJ+2/3 однозначно определяются из условия
постоянства инвариантов Римана вдоль характеристик. Это еще раз подтверждает высказанный выше тезис о том, что в схеме Ошера на этапе решения задачи Римана скачок уплотнения заменяется центрированной волной сжатия.
Примеры практического применения схемы Ошера можно найти в [Chakravarthy, Osher 1983а, 1983b; Osher, Chakravarthy 1983; Hemker, Spekreijse 1986].
1.2.4. Схема Poy
С учетом соотношения (1.2.1.11) противопоточная схема (1.2.1.8) для уравнений Эйлера может быть переписана в виде
-»И+1 -»Л ~7-* "г*
%_ZJH+fj«'2-fj-U2s:0t (12А1)
At Ах
/;./2 =±(+/у+,)+іИУ*у -±HU,*y+I о-2-4-2)
и аналогично для f*_V2 Здесь модуль Якобиана А есть ни что иное, как разность
Якобианов А*
\А\ = А+-А~, (1.2.4.3)
определяемых соотношениями (1.2.1.9), (1.2.1.10). Приближенная для общего случая непостоянства матрицы А замена \Ал и \AJ+l\ на Uy+1/2 дает численный поток вида
fj.m = iifj + /;+і)-іК+і/2|(^+і -*Д (1-2.4.4)
однако при этом теряется свойство консервативности.
Схема Poy [Roe 1981, 1985, 1986] предлагает новый способ расчета матрицы \А\,
обеспечивающий выполнение консервативности. В то же время ее можно трактовать как схему с приближенным решением задачи Римана.
Заметим, что для линейной системы уравнений, когда матрица А всюду является постоянной численный поток (1.2.4.4) будет консервативным, поскольку в этом случае справедливо равенство
/,ч.-Л=Фу+1-*,1 (1-2.4.5)
Основная идея схемы Роу для решения нелинеаризованных уравнений Эйлера
заключается в нахождении матрицы А такой, чтобы она удовлетворяла свойству (1.2.4.5) во всех точках расчетной области, приводя при этом, как и в случае линеаризованных уравнений, к консервативной схеме. Таким образом, требования к
искомой матрице А следующие
1. Для любой пары Wj, wJ+l должно выполняться следующее равенство
7у+1 - fj = 3(wy,wJ+lXwJ+i -Wj) (1.2.4.6)
2. При Wj = w-+l = w матрица А должна совпадать с Якобианом А, т.е.
A(w,w) = A(w) = ^- (1.2.4.7)
3. Матрица А должна иметь только действительные собственные числа, которым
отвечают линейно независимые собственные векторы.
Заметим, что собственные значения матрицы А удовлетворяют соотношению
/y+i-/y=3fc+i-*y)» О-2-4-»)
которое совпадает с соотношением Ренкина-Гюгонио на скачке между состояниями
газа Wj и wJ+l движущемся со скоростью Лг Отсюда прослеживается связь схемы
Роу со схемой Годунова в том смысле, что при решении задачи Римана схема Роу предполагает наличие на каждой грани расчетной ячейки только разрыва состояния газа, тогда как в схеме Годунова задача Римана решается полностью. Здесь уместно вспомнить схему Ошера, в которой скачок уплотнения исключался из этапа решения задачи Римана, отставляя место только гладким ее решениям. Как следствие того, что схема Роу сводит задачу Римана только к разрешению ударных волн, возможно возникновение ударных волн разрежения. Это происходит в тех волнах разрежения, которые содержат звуковые точки.
Что касается структуры матрицы А, то оказывается, что она совпадает с
Л/*
Якобианом вектора потоков А = — , определяемым соотношением (1.2.1.2), если
рассматривать последний как функцию переменных р, и и Н, при условии, что эти переменные усреднены при помощи формул
-I ~ _4p^UJ+\+jPjUJ rr _ jP^Hj+\+jPjHj „ „ л m
Л+1/2 -^PjPj+l ' Uj+\I2 - /= /= ' -"y+1/2 - 7= 7= j(i.^.4.yj
называемых осреднением Роу. Таким образом,
*(*j>*J+ih 4*»ш)> (1.2.4.10)
где w - вектор консервативных переменных, рассматриваемый как функция переменных (1.2.4.9), называемых часто переменными Роу,
Из соотношения (1.2.4.10) следует, что все собственные числа матрицы А -действительные, т.к. они совпадают по форме с собственными числами матрицы А с той лишь разницей, что являются функциями переменных Роу
Л1=и, Л2 = Л+ =и +с , Л3 = Л~ =и —с, (1.2.4.11)
с=^(у-\\Й-\и2) (1.2.4.12)
осредненная скорость звука. Собственным числам (1.2.4.11) отвечают линейно независимые собственные векторы
1 и —
1 =
г2=4г-(і її + с H + ucJ,r3=-^-(\ и-с Н -ucj. (1.2.4.13)
Таким образом, матрица А (1.2.4.10) удовлетворяет всем трем сформулированным ранее требованиям, т.е. обеспечивает получение искомого приближенного решения задачи Римана. В работе [Roe, Pike 1984] показано, что
полученная матрица А является единственной, удовлетворяющей соотношению (1.2.4.6), которое означает, что в каждой ячейке расчетной области исходная консервативная система (1.1.1.1) заменяется на локально линеаризованный ее аналог
(1.2.4.14)
dw ~dw . — +А— = 0, dt дх
что обеспечивает выполнение консервативных свойств. Численный поток схемы Роу аналогичен (1.2.4.4) с той лишь разницей, что матрица А заменяется матрицей А
&1/2 =i(/y + fj+x)'i\^J+mpj+x -). (1.2.4.15)
Для случая нелинейного скалярного уравнения переноса схема Роу совпадает со схемой Мурмана-Коле (1.2.3,3), которая, как следует из предыдущего параграфа, допускает возникновение ударных волн разрежения.
Для предотвращения возникновения ударных волн разрежения в работе [Harten, Hyman 1983] предложено на этапе приближенного решения задачи Римана вводить волну разрежения, когда фиксируется, что эта волна имеет звуковой переход. Это достигается посредством модификации абсолютных значений собственных чисел следующим образом
если
> є
л=\
є + -
2\
если
(1.2.4.16)
25 где параметр є определяется как
є = max(o,Iy+1/2 -ЯуДу+1 -XJ+m). (1.2.4.17)
13. СХЕМЫ С ИСКУССТВЕННОЙ ДИССИПАЦИЕЙ. CUSP СХЕМА
1,3.1. Концепция искусственной диссипации
Концепция схем с искусственной диссипацией появилась как средство подавления нефизических осцилляции в центрально-разностных схемах повышенного порядка аппроксимации. В качестве примера центрально-разностной схемы приведем схему Лакса-Вендроффа [Lax, Wendroff 1960], численный поток которой имеет вид
/Al/2 =\(fj+l+l)-\^J+»l{fM-fj\ (1.3.1.1)
Эта схема имеет второй порядок точности по пространству, вследствие чего в окрестности ударных волн генерируются высокочастотные осцилляции, а второго слагаемого в правой части (1.3.1.1), представляющего собой численную диссипацию, оказывается недостаточно для их подавления. Для устранения нефизических осцилляции Фон Нейман и Рихтмайер [Von Neumann, Richtmyer 1950] ввели понятие искусственной диссипации, которая должна моделировать эффект физической вязкости в окрестности газодинамических разрывов и быть незначительной, т.е. порядка точности схемы или выше, в областях гладкого изменения газодинамических функций. (Здесь и далее под численной диссипацией понимается дополнительная добавка к центрально-разностной дискретизации потоков безотносительно к способу его появления в схеме. Термин искусственная диссипация употребляется для выделения того факта, что численная диссипация появляется умышленным образом, в отличие о случаев ее естественного появления в результате того или иного способа дискретизации). Наличие достаточной численной диссипации также необходимо для выполнения энтропийного условия [Lax 1973], которое гарантирует отсутствие в решении ударных волн разрежения. Рихтмайером и Фон Нейманом был предложен диссипативный поток, который может быть приложен для модификации любой центрально-разностной схемы следующим образом
HAD)* _
/Ж =/,41,2 -*Д*Р
1 Щи^-uj). (1.3.1.2)
\UJ
Лаке и Вендрофф [Lax, Wendroff 1960] предложили обобщенную численную диссипацию в форме
f$Z=fj\m -Ф,.*,Л*/+і ~*Д (1-3.1.3)
где D - любая положительная функция разности wJ+] - Wj, пропорциональная, по
крайней мере, первой степени wJ+l - Wj.
Всевозможные формы диссипативных потоков были предложены различными авторами за прошедшие десятилетия. Так в работе [MacCormak, Baldwin 1975] был использован диссипативный поток третьего порядка точности
(1.3.1.4)
/(^*=/*J^M±
Р В работе [Steger 1978] была предложена численная диссипация четвертого порядка
/<^'=/Ч^| + с)0. (1.3.1.5)
Основным недостатком схем (1.3.1.2), (1.3.1.4), (1.3.1.5) является наличие коэффициента є, оптимальный выбор которого может зависеть от конкретной задачи, а также весьма грубое разрешение ударных волн. Отмеченные схемы редко используются в последнее время, поскольку противопоточные методы обладают существенными преимуществами.
Тем не менее, концепция искусственной диссипации получила удачное развитие в последнее десятилетие, благодаря работам Джемесона. Рассматриваемая ниже CUSP схема сочетает в себе точность, присущую противопоточным схемам, и экономичность схем с численной диссипацией.
1.3.2. Базовые положения CUSP схемы
Отправной точкой в конструировании диссипативных схем Джемесона является тот факт, что противопоточные схемы можно трактовать как схемы с численной диссипацией. Действительно, рассмотренная в предыдущем параграфе схема Роу может трактоваться как схема с численной диссипацией [Jameson 1995с], т.е. схема в которой численный поток представляется в виде суммы центрально-разностного и диссипативного потоков
Л+1/2=|(й+//+іМу+1/2> 0-3.2.1)
где, согласно (1.2.4.15)
dJ+u2 = ipy+]/2|(wy+1 -Wj). (1.3.2.2)
Достоинство схемы Роу состоит в том, что для предотвращения возникновения осцилляции она вводит минимальную численную диссипацию и, что с использованием линеаризации Роу допускается существование дискретной ударной волны с единственной внутренней точкой. Однако схема требует матричных вычислений, что оказывается трудоемким при проведении трехмерных расчетов. Для
27 уменьшения вычислительных затрат в формуле (1.3.2.2) можно заменить AJ+y2 па скалярный коэффициент aJ+y2
dj+уг = aj+y2&wJ+y2, (1.3.2.3)
где aJ+y2 определяется как минимальное по модулю собственное число матрицы
Aj+y2. Схема со скалярным диссипативным потоком разрешает ударные волны
примерно с тремя внутренними точками. Она широко используется при расчетах трансзвуковых течений, так как является экономичной и надежной [Jameson 1995а].
Промежуточный по трудоемкости класс схем может быть сформулирован путем определения диссипативного потока как комбинации разностей консервативных переменных и векторов потоков
dM2 = 2(^4+1/2(^4. -wy)+-/?y+1/2(/;.+1 -/у). (1.3.2.4)
Вектор потоков / далее может быть разложен на конвективную составляющую и составляющую, отвечающую за распространение волн давления
f = uw+fp, (1.3.2.5)
Л =(0 Р ир)т. (1.3.2.6)
Тогда
dJ+u2 = 2^)/+1/2(^41 -^/)+ jfij+uJ^faj+i -uj)+\fpj+l-fPj)l (1.3.2.7)
(«ЮуЧІ/2 = («*4+1/2 + Л+1/2Й/+1/2> (1.3.2.8)
где її и w - средние арифметические
u=i(uJ+l +uj), $=i(wJ+l +Wj). (1.3.2.9)
Схема, использующая диссипацию в форме (1.3.2.4), впервые была предложена Джемесоном [Jameson 1993] и называется CUSP (Convective Upwind Split Pressure) схемой.
Идея расщепления потока на конвективную составляющую и составляющую с давлением восходит к, так называемой, AUSM (Advection Upstream Splitting Method) схеме [Liou, Steffan 1993]. Смысл подобного расщепления состоит в том, что в дозвуковом потоке распространение волн давления осуществляется вверх по потоку, тогда как возмущения переносимые за счет конвекции распространяются вниз по потоку. В AUSM схеме подобный механизм распространения возмущений реализуется за счет аппроксимации соответствующими разностями (вверх или вниз по потоку) различных составляющих потока, полученных в результате расщепления.
Для дальнейшего изложения диссипативный поток удобно представить в обобщенном виде
.(*у+1-*Д (1.3.2.10)
dj+yi - 2Dj+\/2\
в случае CUSP
V+l/2
Для схемы Роу матрица DJ+y2 есть модуль матрицы Роу схемы матрица DJ+y2 имеет вид
DJ+u2 =Pj,mAj+m+{ac)^mI, (1.3.2.11)
что следует из соотношения (1.3.2.4) после подстановки в него линеаризации Роу (1.2.4.6).
Забегая вперед, отметим, что согласно [Swanson et al. 1997] схема с численной диссипацией в форме (1.3.2.10) оказывается TVD схемой при условии справедливости неравенства
/1/2 *
^/+1/2
(1.3.2.12)
Отсюда сразу следует, что схема Роу является TVD. Для CUSP схемы неравенство (1.3.2.12) не выполняется, а потому она не является TVD схемой.
j+1 j+2
Рисунок 1.2. Схема скачка с единственной внутренней точкой
Коэффициенты CUSP схемы определяются из условий существования ударной волны с единственной внутренней точкой (рис. 1.2). Согласно [Jameson 1995b] это возможно, в случае если численная диссипация удовлетворяет следующим двум условиям 1. В сверхзвуковых областях течения численная диссипация такова, что производимый ею вклад в численный поток аналогичен его противопоточному вычислению; 2. На выходе из "области скачка" численная диссипация удовлетворяет обобщенной задаче на собственные значения
{aar+Dar\wr-wa)=0. (1.3.2.13)
Здесь AAR - это матрица Роу, a DAR - матрица, определяющая численную диссипацию на грани AR.
Схема Роу удовлетворяет обоим условиям, тогда как схема со скалярной диссипацией не удовлетворяет первому условию.
Для коэффициентов CUSP схемы из первого условия следует
а*=0, р = ъ%ъ{м), М>\. (1.3.2.14)
Подставляя в соотношение (1.3.2.13) матрицу D, определяемую выражением (1.3.2.11), имеем
*AR
а с т
(wR-wA)=0,
29 (1.3.2.15)
величина
Отсюда следует, что вектор wR - wA является собственным вектором матрицы AAR, а ас
\ + /3
ее собственным значением, т.е.
Л =
а с
(1.3.2.16)
где Л - одно из собственных чисел определяемых по (1.2.4.11). Определяя из последнего соотношения величину а*с и подставляя его в (1.3.2.11) получим
а)
следующее выражение для диссипативной
матрицы D
/+1/2=-^7+1/2^' + ^+1/2^+1/2-^+1/2^-(1-3.2.17)
Для получения положительной диссипации в
случае и > 0 собственное число Л должно быть отрицательным, т.е.
а* с =-(1 + /3)1-. (1.3.2.18)
б)
Аналогично для случая и < О
а'
а*с={\ + р)Л\ (1.3.2.19)
Таким образом, соотношение (1.3.2.16) позволяет построить семейство однопараметрических схем. Для коэффициента /? Джемесоы [Jameson 1995b] предложил следующие значения
в)
+ max
и -Л
-max О,
и -Л+ sgn(M)
>1
Рисунок 1.3. Коэффициенты CUSP схемы
Ограничения /3 > 0 при и > 0 и /? < 0 при и < 0 имеет следствием центрально-разностную аппроксимацию потока / при числах Маха меньших 1/2.
Скачок разрешается с единственной внутренней точкой, если положить
(1.3.2.21)
Таким образом, коэффициенты (1.3.2.20) и (1.3.2.21) обеспечивает полную противопоточность при сверхзвуковом течении (/? = sgn(M), а = 0). Выбор а*с = \и\
при (5 = 0 приводит к непрерывному распределению коэффициента а с в дозвуковой области и низкому уровню диссипации в окрестности точки торможения, что является привлекательным свойством при расчете вязких течений. Коэффициенты а, а и J3 показаны на рисунке 1.3.
Заметим, что при числах Маха меньших 1/2, когда /3 = 0 соотношения (1.3.2.18), (1.3.2.19) не выполняются. Следовательно, разрешение скачка с единственной внутренней точкой и отсутствие осцилляции гарантируется только для скачков, число Маха за которыми не меньше 1/2. Коэффициенты CUSP схемы для расчета сильных скачков, число Маха за которыми не меньше 1/3 могут быть переопределены следующим образом [Swanson et al. 1997]
+ max
її + ^Я" и -Л
J3 = {-max 0, ^
ЇЇ-Л+
sgn(M)
>1
(1.3.2.22)
а с =<
2|н|, Р = 0
-(1 + /?)1~, /?>0, 0<М<1
(i-j3)X+, /?<о, -1<м<о'
(1.3.2.23)
0,
>1
м м
При выборе коэффициентов CUSP схемы возможны некоторые упрощения. Например, в работах [Jameson 1995b; Tatsumi et al. 1995] использовались следующие значения для коэффициента а
>:
М.
а = <
(
є + -
(1.3.2.24)
тогда как fi вычисляется по тем же соотношениям (1.3.2.20). Выбор коэффициента а (1.3.2.24) не удовлетворяет условиям (1.3.2.18), (1.3.2.19), поэтому скачок с единственной внутренней точкой существовать не может.
Идя дальше по пути упрощения расчетных формул, осреднение Роу можно заменить обычным осреднением, однако в этом случае возможно появление осцилляции в окрестности разрывов [Swanson et al. 1997].
Собственные числа матрицы D (1.3.2.11) //,с, //2с, //Зс , где
Mi =
М2 =
М + ll, М^\
м, м<\
a-j3, ^<М<\
>1,
М-1,
(1.3.2.25)
показаны на рисунке 1.4. В области М <-j все собственные числа равны между
собой и CUSP схема становится аналогичной схеме со скалярной диссипацией. При
этом
диссипативныи
коэффициент
пропорционален
, что делает CUSP схему
\
/
/
/
привлекательной для проведения расчетов в существенно дозвуковых областях с применением методов регуляризации исходных уравнений (см.
CUSP схема, обеспечивающая постоянство полной энтальпии в стационарном потоке
В стационарном невязком потоке энтальпия торможения Н должна быть постоянна. Из этого положения следует, в частности, тот факт, что стационарные уравнения сохранения массы и энергии идентичны, если постоянный множитель Н вынесен из-под знаков производных в уравнении энергии. Дискретные и полудискретные схемы не всегда удовлетворяют этому свойству. В случае численного потока в форме (1.3.2.1) допускается решение с постоянной энтальпией Я, если численная диссипация для уравнения энергии приводится к аналогичной численной диссипации для уравнения неразрывности, в котором р заменено рН. Стандартные схемы Роу и CUSP не обладают этим свойством, и в дискретном решении полная энтальпия не постоянна. Источник ошибки в энтальпии торможения находится в расхождении между конвективным членом в векторе потоков который содержит рН и вектором состояния, который содержит рЕ. Это расхождение может быть устранено введением модифицированного вектора состояния Принимая последнее, введем линеаризацию Роу где матрица Ah, отличная от стандартной матрицы Роу, записывается в виде [Jameson 1995b] Заметим, что Л и Л имеют те же знаки, что и+с и и-с , и меняют знак при прохождении звуковых точек и = ±с . Диссипативный поток CUSP схемы теперь выражается следующим образом Численные схемы, рассмотренные в предыдущих параграфах, аппроксимировали исходную систему уравнений с первым порядком точности по пространству.
На практике этого оказывается недостаточно, вследствие чего необходимо переходить к схемам второго, а иногда и более высоких порядков. Однако непосредственное обобщение изложенных методов на случай повышенного порядка аппроксимации приводит к возникновению нефизических осцилляции в окрестности ударных волн, которые не удается устранить даже путем измельчения расчетной сетки. Подобные осцилляции могут отсутствовать на слабых скачках уплотнения [Hirsch 1990], что является следствием нелинейности уравнений Эйлера, однако, в общем случае это не так, поскольку можно показать [Engquist, Osher 1981], что линейные схемы второго порядка всегда генерируют осцилляции. Рассмотрим в начале способы повышения порядка аппроксимации, после чего обратимся к методам борьбы с осцилляциями. Вспомним, что метод Годунова, а также методы, основанные на приближенном решении задачи Римапа, оперировали с осредненными по ячейке параметрами потока, причем последние полагались постоянными в пределах ячейки. Такое кусочно-постоянное распределение газодинамических переменных приводит к первому порядку аппроксимации. Соответственно, если перейти к кусочно-линейному распределению, то результирующая схема будет иметь второй порядок точности. Заметим, что это никак не отразится на методе решения задачи Римана, поскольку эти два этапа схемы являются независимыми друг от друга. В задаче о распаде разрыва будут фигурировать значения переменных, полученные в результате экстраполяции на поверхность разрыва осредненных значений в соседних ячейках. Такой подход, предложенный в работе [Van Leer 1979] и заключающийся в повышении порядка аппроксимации за счет экстраполяции газодинамических переменных получил название MUSCL (Monotone Upstream-centered Schemes for Conservation Laws). Метод конечного объема, на основе которого проводится дискретизация исходных уравнений в данной работе, оперирует с осредненными в пределах расчетной ячейки величинами. Поэтому линейное или квадратичное (в схемах третьего порядка аппроксимации) распределения газодинамических переменных в пределах ячейки также зависят от осредненных параметров, а не от значений в центре расчетной ячейки [Hirsch 1990] - осредненное по ячейке значение функции v(x). Формула (1.4.1.1) получается из разложения в ряд Тейлора функции v(x), в котором принято во внимание следующее соотношение, выражающее связь между средним значением функции v(x) и ее значением в центре ячейки
Последнее соотношение, в свою очередь, получается в результате интегрирования разложения функции v(x) в пределах ячейки. При к = )/ъ соотношение (1.4.1.1) представляет собой разложение с точностью до величин третьего порядка, и эта параболическая интерполяция приводит к третьему порядку точности аппроксимации по пространству. Другие значения к соответствуют различным видам остаточных членов при линейной интерполяции. Выражая первую и вторую производные посредством центрально-разностных формул найдем значения v(x) слева и справа от грани у+ 1/2, которые позднее будут фигурировать в задаче Римана как значения по разные стороны разрыва Параметр введен в интерполяционные формулы (1.4.1.6) для возможности переключения на первый порядок точности в случаях, о которых будет сказано в дальнейшем. Возвращаясь к уравнениям Эйлера, после нахождения новых значений консервативных переменных Wj+y2, wf+y2 по формулам (1.4.1.6), численный поток любой из рассмотренных ранее схем можно записать с повышенным порядком точности следующим образом
Методы регуляризации уравнений газовой динамики
Решение уравнений Эйлера в задачах с малыми числами Маха осложняется значительным возрастанием времени, необходимого для достижения установившегося решения. При использовании метода установления допустимый локальный шаг интегрирования по времени обратно пропорционален максимальному собственному числу системы. В слабо сжимаемых течениях максимальное собственное число приблизительно равно скорости звука, тогда как минимальное -скорости течения. По этой причине за один временной шаг возмущения, переносимые с конвективной скоростью, не могут распространиться далеко, и достижение сходимости требует большого числа итераций. Различные приемы для ускорения сходимости, такие как многосеточные методы или неявные схемы, не решают этой проблемы.
Техника, известная как регуляризация исходных уравнений, помогает разрешить проблему сходимости для слабо сжимаемых течений. В последние годы появилось несколько методов регуляризации [Стрелец, Шур 1988; Choi, Merkle 1993; Turkel 1993; Turkel et al. 1996; Turkel et al. 1997], нацеленных на моделирование слабо сжимаемых течений численными методами, разработанными для существенно сжимаемых потоков. Востребованность методов регуляризации обусловлена двумя важными моментами. Во-первых, существуют течения, в которых одновременно присутствуют области как существенно сжимаемого, так и слабо сжимаемого течения газа. Во-вторых, предпочтительно использовать существующий и уже отлаженный код по решению сжимаемых задач для более широкого круга течений, чем специально разрабатывать и поддерживать отдельный код для решения несжимаемых задач.
За последнее время сложилось два основных способа реализации методов регуляризации. Первый рассматривает регуляризацию только как средство ускорения сходимости, что достигается за счет умножения производных по времени на специальную матрицу. Введение такой операции приводит к понижению максимального собственного числа системы и увеличению допустимого шага интегрирования по времени. При этом никак не затрагивается пространственная дискретизация уравнений. Сошедшееся решение регуляризованной таким образом системы полностью идентично решению исходной системы. Второй способ подразумевает изменение дискретного аналога пространственного оператора исходной системы уравнений, а не только производных по времени. При таком подходе регуляризация вводится на этапе дискретизации, при вычислении численной диссипации, а сошедшееся решение уже не будет идентично решению исходной системы. Обоснование этого подхода можно найти в работе [Turkel et al. 1994], где показано, что стандартные численные схемы для решения сжимаемых уравнений в пределе при стремлении числа Маха к нулю не приводят к схемам решения несжимаемых уравнений, тогда как регуляризованные вторым способом схемы оказываются аналогичными схемам для несжимаемых уравнений.
Будем решать стационарную систему уравнений Эйлера (1.1.1.1) методом установления. Поскольку нас интересует только стационарное решение мы вправе модифицировать производные по времени по своему усмотрению с целью достижения наилучшей скорости сходимости. В работе [Turkel 1985] показано, что наилучшая скорость сходимости достигается, когда число обусловленности системы оказывается порядка единицы, т.е. все собственные числа равны между собой. Метод регуляризация заключается в замене исходной системы (1.1.1.1) регуляризованной системой вида с целью уменьшения числа обусловленности.
В дальнейшем при изучении свойств матрицы Р 1 будем заниматься системой, записанной в неконсервативной форме где q - вектор неконсервативных переменных, A = —L-J- _ матрица Якоби, что, тем не менее, не окажет влияния на способ пространственной дискретизации, которая по-прежнему будет проводиться с использованием консервативной формы уравнений. Единственное ограничение, накладываемое на матрицу регуляризации Р 1 - это существование обратной матрицы. В работе [Turkel 1985] показано, что для того, чтобы регуляризованная задача была корректно поставленной (в частности, в смысле граничных условий), желательно, чтобы система была симметричной гиперболической системой (или приводилась к таковой). Для симметричной гиперболической системы в случае положительной определенности матрицы регуляризации (положительно определенная матрица - матрица, все собственные числа которой больше нуля) гарантируется, что регуляризованная система (1.5.1.1) имеет то же число граничных условий, что и исходная система (1.1.1.1) и, что процесс регуляризации не приводит к обращению времени в отдельных уравнениях системы. Если матрица Р х не является положительно определенной, то изменяется физическая постановка задачи, так, например, матрица Р х =-1 приводит к обращению времени и, следовательно, к изменению граничных условий задачи. Для систем более общего вида необходимо проводить детальный анализ для доказательства идентичности решений двух систем [Turkel 1998].
Геометрические характеристики и индексация конечных объемов
В настоящей работе для численного решения системы уравнений Навье-Стокса (2.1.2.12) используется метод контрольных (конечных) объемов [Hirsch 1990; Флетчер 1991]. Расчетные схемы, построенные на основе этого метода, обладают свойством консервативности, обеспечение которого крайне желательно, а для течений с ударными волнами просто необходимо. При расчете трехмерных течений на структурированных сетках в качестве конечных объемов используются восьмивершинники. Предположим, что в физической области построена структурированная сетка и сеточные линии являются координатными линиями обобщенной системы координат (,, 2, 3) Пронумеруем узлы сетки дробными индексами / +1/2, j +1/2, k + \/2, а центры восьмивершинников, построенных таким образом, чтобы узлы расчетной сетки являлись их вершинами, целыми индексами /, j, к. Заметим, что может реализоваться случай, когда точки, являющиеся вершинами граней восьмивершинника не лежат в одной плоскости (рис. 2.2, точки А, В, С, D). Для реализации метода конечных объемов необходимо знать некоторые геометрические характеристики рассматриваемых восьмивершинников, а именно: его объем и вектора площадей боковых граней. Рассмотрим боковую грань EFGH (рис. 2.3). Определим вектор площади грани SEFGH = S xl2 как вектор равный по модулю ее площади и ортогональный к изоповерхности 2 = const. Вектор площади грани SEFGH складывается из векторов площадей двух треугольников EFG и GHA. Их можно выразить в виде соответствующих векторных произведений В соотношении (2.1.3.1) г - это радиус-вектор соответствующего узла сетки. Аналогичные выражения можно записать для векторов S 2 , S k+ 2 . Объем восьмивершинника найдем как сумму объемов трех пирамид DABFE, DEFGH, DBCGF.
Имеем Заметим, что если обобщенная система координат «привязана» к расчетной сетке так, что в вычислительной области контрольные объемы имеют единичный размер по всем направлениям (А = Д2 = А3 = 1), то якобиан обратного преобразования у-1 вычисленный в центре контрольного объема равен объему ячейки V, а введенные в предыдущем разделе метрические коэффициенты %та имеют наглядную геометрическую интерпретацию: вектор с декартовыми компонентами \%т і К т 2 JX т,з J Л ) рассчитанный для центра грани, есть вектор площади грани AV"1 , m = i,j,k. В этом случае величина G m в соотношении (2.1.2.14) есть ни что иное как компонента вектора массового расхода через грань . В рамках рассматриваемых моделей турбулентности предполагается, что для расчета компонент тензора рейнольдсовых напряжений может быть использована гипотеза Буссинеска [Лойцянский 1987] с изотропным представлением кинематического коэффициента турбулентной вязкости Турбулентная вязкость моделируется по одной из рассматриваемых ниже моделей турбулентности, выражение для ее определения зависит от конкретных приближенного описания процесса перехода от ламинарного к турбулентному режиму течения в заданном месте пристенного слоя, когда внешний поток не турбулентный. Даже при малой степени турбулентности внешнего потока турбулентная вязкость в пристенном слое генерируется автоматически, и в этом случае функции /п и/,2 можно положить равными нулю. Модель Спаларта-Аллмараса сформулирована в низкорейнольдсовой постановке, что означает, что она описывает поведение эффективной вязкости во всей области течения, включая пристенную область и вязкий подслой. На стенке v полагается равной нулю.
Дискретизация конвективной части стабилизирующего оператора по схеме расщепления матрицы коэффициентов
Для ускорения сходимости к установившемуся состоянию в каждой ячейке сетки можно использовать свой шаг по времени установления (псевдовремени), зависящий от локальной скорости распространения возмущений и от размеров ячейки. В каждой точке промежуточное решение будет приближаться к установившемуся решению с максимально возможной для данного числа Куранта скоростью. При этом даже для нерегуляризованных уравнений промежуточное решение теряет смысл как решение физической нестационарной задачи. Для простой противопоточной схемы первого порядка, аппроксимирующей одномерное уравнение Эйлера, анализ устойчивости по методу Неймана приводит к следующему ограничению на шаг по псевдовремени [Флетчер 1991] Максимальное собственное значение матрицы Іг1 получается из соотношений (2.4.1.9) и имеет вид 1 Ограничение на шаг типа (2.4.8.1) приводит к следующей формуле для вычисления локального шага по псевдовремени где схемный параметр CFL - это число Куранта-Фридрихса-Леви (сокращенно -число Куранта), задаваемое пользователем программы на основе опыта аналогичных расчетов. На трехмерный случай формулу (2.4.8.3) можно обобщить разными способами. Наиболее распространенный, и используемый в настоящей работе, основан на непосредственном использовании результатов одномерного анализа. Рассмотрим направление ,. Оценки скорости и шага в окрестности ячейки в направлении , следующие Аналогичные выражения получаются для локальных шагов Ат2, Аг3 при рассмотрении направлений 2, 3 Далее выбирается минимальный шаг Применение соотношения (2.4.8.6), полученного для уравнений Эйлера, при решении уравнений Навье-Стокса может приводить к завышению допустимого временного шага, особенно в тех областях, где вязкие эффекты преобладают над конвективными. Для учета вязких эффектов в локальный шаг по псевдовремени вводится соответствующая модификация [Cabuk и др. 1992] где допустимый временной шаг, обусловленный вязкими эффектами, определяется соотношением Здесь Cvis - задаваемый параметр - аналог числа CFL для учета вклада вязких эффектов.
При вычислении локального шага по псевдовремени на этапе интегрирования уравнений турбулентности используется только одномерный подход, где в качестве собственного числа Ятах берется скорость в соответствующем направлении. Также проводится дополнительная модификация Атт с целью улучшения диагонального преобладания за счет учета источниковых членов уравнений турбулентности [Smirnov 1993]. При моделировании реальных течений в проточных частях турбомашин приходится иметь дело с конфигурациями сложной геометрии, будь то выходной патрубок турбины или лопатка компрессора. Построение расчетной сетки занимает существенную часть времени, отведенного на решение задачи, и требует особой тщательности, поскольку во многом определяет качество полученных результатов. Часто построение единой структурированной расчетной сетки значительно осложнено, а иногда и просто невозможно. Для преодоления указанных трудностей в настоящей работе применяются блочно-структурированные сетки, использование которых значительно облегчает и сокращает процесс построения расчетной сетки. Кроме того, этот метод позволяет экономить машинную память, и эффективен в случае организации параллельных вычислений. Многоблочный подход начал разрабатываться в начале 80-х годов [Lee, Rubbert 1980; Rai 1986] и в настоящее время является весьма популярным.
Основная идея метода состоит в разделении расчетной области на подобласти (блоки), дискретизация в каждом из которых проводится отдельно, причем в общем случае блоки могут частично перекрываться. В случае неперекрывающихся сеток реализация алгоритма существенно зависит от того, продолжаются ли сеточные линии при переходе из блока в блок. В любом случае, в результате дискретизации получается набор структурированных сеток, число которых равно количеству блоков. Вычисления в ходе одной итерации на каждой из сеток ведутся независимо, после каждой итерации производится обмен информацией на межблочных границах. В рамках настоящей работы применялись два из реализованных в ПК SINF многоблочных алгоритмов. Оба алгоритма подразумевают, что стыкуемые блоки не перекрываются. В первом сетки строятся так, что сеточные линии продолжаются при переходе из одного блока в другой (рис. 2.4), во втором допускается обрыв сеточных линий на интерфейсе стыковки, но подразумевается, что к нескольким ячейкам одного блока пристыковывается только одна ячейка другого (рис. 2.5). Индексные координатные линии в стыкуемых блоках могут быть ориентированы произвольным образом. С целью обмена информацией между блоками для каждой стыковки создается отдельный виртуальный блок. Применение виртуального блока позволяет произвести стыковку блоков с примыкающими к границе скошенными ячейками (то есть в том случае, когда сеточные линии на границе двух блоков терпят излом) без нарушения свойства консервативности. Виртуальный блок имеет свою внутреннюю индексную систему координат, его размеры определяются размерами сегментов связи, задаваемых индексными координатами двух противолежащих вершин интерфейса для каждого из двух стыкуемых блоков. Реализованный алгоритм многоблочных вычислений подразумевает, что виртуальный блок состоит из двух слоев приграничных ячеек каждого из стыкуемых блоков (рис. 2.6). В общем случае число требуемых слоев, взятых из обоих блоков, может быть как меньшим, так и большим и зависит от характера вычислений, производимых внутри виртуального блока. Блоку соотносится информация о значениях переносимых величин в используемых ячейках, коэффициентах переноса и вся информация о геометрии ячеек. Таким образом, виртуальный блок представляет собой вполне независимый объект, в котором вычисляются потоки через площадки, составляющие сегмент связи двух реальных блоков, и значения самой переносимой величины на границе. Эти вычисления проводятся по тому же алгоритму, что и внутри стыкуемых блоков. Полученные результаты рассылаются по соответствующим заграничным ячейкам исходных ("реальных") блоков и используются при вычислении невязок на следующей итерации.
