ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2007, том 47, < 8, с. 1387-1401
УДК 519.6:531.33
СОСТАВНЫЕ КОМПАКТНЫЕ СХЕМЫ ВЫСОКОГО ПОРЯДКА ДЛЯ МОДЕЛИРОВАНИЯ ТЕЧЕНИЙ ВЯЗКОГО ГАЗА1}
© 2007 г. А. Д. Савельев
(119991 Москва, ул. Вавилова, 40, ВЦ РАН) e-mail: savel@ccas.ru
Поступила в редакцию 15.11.2004 г.
Переработанный вариант 02.10.2006 г.
Рассмотрены дифференциальные схемы до 7-го порядка включительно для численного описания законов сохранения на криволинейных сетках. Они представляют собой комбинацию симметричных компактных разностей и ориентированных в соответствии с направлением характеристик операторов диффузного типа. Рассмотрены спектральные свойства схем, проведены расчеты ряда модельных задач. Разработанный метод численного интегрирования уравнений Навье-Стокса с применением двухпараметрической модели турбулентности применяется для моделирования двумерных течений вязкого газа. Библ. 21. Фиг. 11.
Ключевые слова: компактные аппроксимации, схемы высокого порядка, криволинейные координаты, турбулентные течения вязкого газа.
ВВЕДЕНИЕ
Постоянно повышающиеся требования к уровню научно-технических разработок делают актуальной задачу совершенствования используемых численных методов. В настоящее время они применяются при исследовании проблем турбулентности, аэродинамики, теплообмена, акустики и т.д. Улучшить точность проводимых расчетов призваны методы высокого порядка аппроксимации. Поскольку повышение порядка неизбежно ведет к укрупнению сеточного шаблона, одним из перспективных направлений разработки численных алгоритмов является использование компактных аппроксимаций. Для задач аэрогазодинамики ориентированные против потока компактные разности третьего порядка успешно применяются с середины 70-х годов (см. [1]).
В настоящее время существуют два подхода к применению компактных разностей высокого порядка на неразнесенных расчетных сетках. В первом случае это симметричные разности высокого порядка, дополненные различными фильтрами для подавления схемных осцилляций (см. [2]-[4]). Несмотря на кажущуюся простоту, реализация данного подхода при решении задач с достаточно сильными градиентами параметров потока может представлять серьезную проблему. Выбор механизма демпфирования паразитных осцилляций в ряде случаев составляет предмет отдельного исследования.
Другой подход состоит в использовании несимметричных компактных аппроксимаций, ориентированных против потока (см. [1], [5], [6]). Разработанный изначально для схем третьего порядка, он применяется для получения схем произвольного порядка на основе мультиоператорного принципа. В схемах этого класса содержится естественный диссипативный механизм подавления схемных осцилляций. К сожалению, использование подобных схем в областях с криволинейными границами вызывает вопросы, связанные, например, с точностью сохранения невозмущенного потока. Альтернативный вариант применения несимметричных компактных разностей предложен в [7] как консервативный метод конечных объемов высокого порядка аппроксимации.
Большого различия между двумя упомянутыми подходами нет. В [8] убедительно показано, что известные разности, ориентированные против потока, легко могут быть получены из центральных путем добавления операторов, учитывающих направление характеристик. В [9] такие схемы на основе компактных аппроксимаций высокого порядка были построены. При этом учитывались следующие аспекты моделирования течений вязкого сжимаемого газа компактными аппроксимациями.
1) Работа выполнена при финансовой поддержке РФФИ (код проекта 05-01-00584).
Преодоление сеточных негладкостей. Узлы разностных сеток, адаптированных к геометрии обтекаемого тела, как правило, распределяются в физической плоскости неравномерно. Зачастую нельзя обойтись без их сильного сгущения в направлении твердой поверхности, кривизны и скоса координатных линий. Влияние данных факторов на разностное решение сильно зависит от способа вычисления метрических коэффициентов. Другой аспект проведения расчетов на криволинейных сетках состоит в проблеме сохранения невозмущенного потока. Появление даже очень слабых возмущений параметров течения в сверхзвуковом набегающем потоке следует признать очень неприятным явлением. Как показано в [4], для преодоления упомянутых трудностей необходимо при определении метрических коэффициентов расчетной области использовать те же дифференциальные операторы, что и для аппроксимации уравнений законов сохранения. В трехмерном случае этого недостаточно, однако и здесь существуют подходы, позволяющие устранить нежелательные моменты, вызванные криволинейными координатами.
Эффективность диссипативного механизма. Симметричные разностные операторы не содержат никакого механизма подавления высокочастотных гармоник, не имеющих физического смысла. Введение в уравнения членов типа искусственной вязкости не всегда удобно, поскольку не позволяет проводить сквозной счет разрывных решений. С другой стороны, использование ограничителей потока вида [10] может приводить к искажению результатов, например, к отсутствию отрыва пограничного слоя там, где эксперимент фиксирует его наличие. В этом плане схемы вида [5], [6], [11], обладающие естественными диссипативными свойствами за счет ориентации разностей против потока, имеют заметные преимущества. Применяющийся в [7] механизм монотонизации решения также представляет собой ориентацию конечных разностей с учетом направления характеристик. Можно сделать вывод о преимуществах схем, ориентированных против потока, хотя способы реализации данного подхода могут отличаться.
Разностное представление диффузных членов. С повышением порядка аппроксимации схем, описывающих конвективные составляющие уравнений, встает вопрос о более точном представлении вязких членов. Привлекательный на первый взгляд алгоритм вычисления вязких членов с использованием симметричных компактных разностей успешно применялся в [4]. Однако опыт [9] показал, что это оправдывает себя лишь при сглаживании полей физических переменных по методу искусственной вязкости. В данном случае был использован традиционный подход, когда необходимые производные параметров потока определяются в полуцелых узлах сетки. Применение интерполяционных формул высокого порядка дало хорошие результаты.
Ниже рассмотрены разностные схемы, являющиеся модификацией [9], и представлены результаты расчетов ряда течений вязкого газа.
1. ПОСТАНОВКА ЗАДАЧИ
Рассматривается задача обтекания двумерного (плоского) тела турбулентным потоком сжимаемого газа с постоянным отношением удельных теплоемкостей. Решение определяется путем численного интегрирования осредненных по Рейнольдсу нестационарных двумерных уравнений Навье-Стокса. Задание турбулентных составляющих осуществляется в терминах добавочной турбулентной вязкости. Обезразмеренные традиционным способом по параметрам набегающего потока на бесконечности и характерному линейному размеру Н, они имеют следующий вид в декартовой системе координат (х, у):
ди(ф) + ЭЕ( ф) + ЭС( ф) = Ке_1 УЭ[Е 1 (ф, фх) + Г2( ф, фу)] + Э[ С 1 (ф, фх ) + С2( ф, фу)] Л ^ ^ д( Эх Эу V Эх Эу /
р ри р V р
и = р и , Е = р и2 + р , С = р иv 2 р V + р , ф = и
Р V ри V V
е (е + р) и (е + р) V
Г =
0
(4/3 )цих
^ X
(4/3 иих + Vvvx + цРГ1 кх
Г =
0
-(2/3 Vy
Vиу
(2/3 uvy + ^ иу
С, =
0 V V х -(2/3 их Vuvx - (2/3их
С9 =
0
V их
(4/3 ^ Vу V иих + (4/3 )Vvvy + V Рг-1к
Р = р[(У -1 )У-1 к + (2/3) д2 ], е = р[у-1к + 0.5( и2 + V2) + д2 ],
V = V; + Vt, VPг-1 = V; Рг-1
V, Рг-1,
-1
Re = р^и^Я^.
Здесь , - время, р - плотность, и и V - компоненты вектора скорости в направлениях х и у, р -давление, к - энтальпия, е - полная энергия, д = к1/2 (к - кинетическая энергия турбулентности), у - отношение удельных теплоемкостей, V; и V, - коэффициенты молекулярной и вихревой вяз-костей, Рг; = 0.72 и Pгr = 0.9 - ламинарное и турбулентное числа Прандтля, Re - число Рейнольдса. Система дополняется зависимостью молекулярной вязкости от энтальпии в виде формулы Сатерленда (см. [12]).
Для расчета турбулентных характеристик используется дифференциальная двухпараметри-ческая модель (см. [13]). Ее уравнения имеют вид
Эи(ф, ф,(ф, ф,)^ЭСЛФ, ф,) = Ке-1 (д¥п(ф,х) , дС,2(фу)
д ,
дх
ду
дх
ду
и,
и, =
рд , г, = рид , с, = р д , ф, = д
р uv р VV V
(1.2)
=
+ ^Гц1) Цх
(V; + V, Рг;1 )v х
С,9 =
+ V, РгЦ) Цу (V; + V, Рг^ К
И
Я
Я ^
где V = (дЬ 1)1/2, Ь - масштаб турбулентности. Вектор источниковых членов И имеет компоненты
Я = рд(V5 2 В
2 I 2
2 V V
Я = Я V = Т
а( Сё-2в)-pv2
Здесь
5 = 2и\ + (иУ + V2) + 2v2y--D2, В = их + Vy, а = а/ + «2, = 1 - ехр(^Дерд^ 1),
п - расстояние по нормали от поверхности, а С = 0.09, а1 = 0.555, а2 = 0.065, в = 0.92, ^ = 0.0075, Ргд = 1 и Prv = 1.3 - константы. Коэффициент турбулентной вязкости определяется по формуле
V, = ReCv/vP(VJ .
Граничные условия задаются следующим образом. На входной границе при сверхзвуковой скорости набегающего потока фиксируются параметры течения на бесконечности. При дозвуковой скорости течения задаются полное давление, температура торможения и угол вектора скорости, а четвертый параметр определяется на основе характеристических соотношений для одномерного нестационарного течения невязкого газа. На входной границе также фиксируются значения турбулентных параметров дм и V«,. На выходных границах ставятся условия гладкого
вытекания; если же течение дозвуковое, то задается статическое давление. На твердых поверхностях плотность р„ определяется из уравнения неразрывности, ставятся условия прилипания им> = ^^ = 0 и значение энтальпии (температуры) поверхности Для турбулентных параметров используются условия qw = 0 и дv/дn = 0. Также возможно находить vw путем решения на границе второго уравнения системы (1.2).
Для расчета обтекания тел с криволинейной поверхностью на основе исходных уравнений в строгой консервативной форме осуществляется переход к обобщенной системе координат (см. [14]). Область течения в физической плоскости (х, у) проецируется на единичный квадрат расчетной плоскости п) с помощью отображения общего вида \ = ^(х,
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.