Об одном методе численного решения уравнений Навье-Стокса сжимаемого газа Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ Ц А Г И

Т о м III 197 2 Мб

УДК 532.516

ОБ ОДНОМ МЕТОДЕ ЧИСЛЕННОГО РЕШЕНИЯ УРАВНЕНИЙ НАВЬЕ—СТОКСА СЖИМАЕМОГО ГАЗА

А. И. Толстых

Описывается метод численного решения уравнений Навье— Стокса для сжимаемого газа, основанный на растяжении областей с большими градиентами искомых функций. В качестве модельных рассматриваются задачи о течении в ударной волне и об обтекании затупленных тел при достаточно больших числах Рейнольдса. Для последнего случая приводится разностная схема повышенного порядка точности.

1. При решении многих задач аэродинамики часто возникает необходимость использовать уравнения Навье —Стокса. Это происходит либо при малых числах Рейнольдса (например, в случае течений газа малой плотности), либо при больших числах Рейнольдса, когда теория пограничного слоя не дает эффективного описания течения во всей рассматриваемой области. К последнему классу задач следует отнести задачи о взаимодействии скачка уплотнения с пограничным слоем, о следе за телом и т. д.

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

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

Применимость существующих разностных схем для полных уравнений течения вязкого сжимаемого газа [1] —[4] (помимо пере-

численных в обзорах [1] —[2], отметим также работы [3] и [4]), по-видимому, ограничена областью сравнительно небольших значений чисел Re. Ниже описывается метод получения численных решений уравнений Навье —Стокса для различных задач аэродинамики в широком диапазоне чисел Re. Приводятся примеры расчетов, проведенных с целью выяснения эффективности метода в относительно простых случаях.

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

Рассматривая для простоты плоский или осесимметричный случай, выберем ортогональную систему координат s, п таким образом, чтобы линии разрыва, в которые переходят эти зоны при

Re -» оо, пересекали линии s •-= const под ненулевыми углами.

Пусть в общем случае требуется построить решение краевой задачи для уравнений Навье—Стокса в области

0<s<s0, п0 (s)п пг (s) оо, О

где t—время, пй (s)—поверхность тела, a nr(s)—некоторая „внешняя“ граница. Вместо переменных s, п введем новые независимые переменные s, С

s = s\ С = £(и, v, s, п), (1)

где и и v — составляющие скорости, соответствующие координатам sun.

Потребуем, чтобы функция С удовлетворяла следующим условиям:

С |я=я0 — 0; С |я=лг = 1;

У = дудп > 0 («о < п < яг).

Первое условие является условием нормировки, второе означает, что преобразование взаимно однозначно.

Кроме того, выберем функцию С так, чтобы она менялась на свой порядок в областях с большими градиентами, т. е. чтобы она приблизительно повторяла изменение вдоль координаты составляющих скоростей, терпящих разрыв при Re -* оо .

Можно ожидать, что решение исходных дифференциальных уравнений, записанное в координатах (1), будет содержать сравнительно медленно меняющиеся функции переменной С, так что в плоскости s, С исчезнут области с малыми характерными размерами, стремящимися к нулю при возрастании числа Re. В этом случае погрешности, возникающие при замене дифференциальных уравнений разностными, вследствие сравнительно небольших значений высших производных по С от искомых функций должны значительно уменьшиться по сравнению с аналогичными погрешностями при аппроксимации в физической плоскости s, п. Можно предположить также уменьшение немонотонности решений разностных уравнений в случае схем высокого порядка точности.

Отметим, что в отличие от обычного сгущения узлов разностной сетки в областях с большими градиентами функций (напри-

мер, [5]), преобразование (1) не требует априорного знания расположения этих областей; при равномерной сетке вдоль координаты С оно реализует это сгущение „автоматически" в процессе получения решения.

После перехода к координатам (1) в исходных уравнениях появляется новая функция х(я. С) — дС1дп; систему замыкает второе из равенств (1), являющееся определением функции С. Входящая в (1) и в коэффициенты Лямэ для криволинейной системы координат я, п переменная п с достаточной степенью точности может быть определена с помощью квадратуры

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

Примером, когда выбор функции С особенно прост, является задача о структуре прямолинейной ударной волны. Отметим, что эта классическая задача часто является критерием применимости той или иной разностной схемы (например, [6]).

Пусть все функции зависят лишь от одной координаты п, и при п -* оотечение является сверхзвуковым. Отнеся скорость V и плотность р соответственно к скорости Коо и плотности роо невозмущенного потока, давление р и энтальпию к соответственно к величинам р*, 1/^ и I/£,, запишем граничные условия задачи в виде

где Моо — число М невозмущенного потока, f — показатель адиабаты.

Положив С = (у — v0)/(vao — v0), найдем, что функция С удовлетворяет всем сформулированным выше требованиям. Уравнения течения газа в одномерной ударной волне запишутся в виде

где Рг — число Прандтля, <о — постоянная (0,5 < ш < 1); (л — коэффициент вязкости, отнесенный к коэффициенту вязкости р.* при Л=1.

Введем разностную сетку С( = г'Д; (/ = 0, 1,. . . , 1/ДС) и запишем для уравнений (4) заведомо грубую разностную схему 1-го порядка

(2>

о

(3>

при я

Я — ОО ,

(4)

* dn Роо 1/(50 ’ Р

Л F-*

К — 1г. L

_-рЛ; Р==ЛШ,

аппроксимации. Для этого заменим первые производные от всех функций, кроме давления, односторонними разностными отношениями; производные dpjdi будем аппроксимировать с помощью центральных разностей. Аппроксимацию производной [d (pxdfi/d(,)]ldZ запишем с точностью до членов порядка О(ДС), не центрируя коэффициента [V/ относительно узлов, в которых вычисляется производная dhjdС. Система разностных уравнений для системы (4) с граничными условиями (3) и очевидным условием -/(1) = 0 решалась методом последовательных приближений с введением релаксационного параметра; в течение одной итерации каждое уравнение рассматривалось как линейное.

Расчеты проводились при Моо=10, Рг = 0,75, 4—1,4, ш = 0,5. На фиг. 1 приведены зависимости давления, энтальпии и величины % от координаты С при различном числе узлов сетки (N= 5, 10 и 20). Видно, что газодинамические функции всюду имеют сравнительно малую кривизну, так что любая разностная аппроксимация инерционных членов не должна привести к большим погрешностям. Это подтверждается тем, что кривые для функций hup даже в неблагоприятном случае грубой схемы при N=5, 10 и 20 сравнительно мало отличаются друг от друга. Некоторое различие кривых для функции у можно объяснить первым порядком аппроксимации члена с коэффициентом вязкости во втором уравнении системы (4).

Для перехода в физическую плоскость поместим начало координат в точку С = 1/2 и будем вычислять определенные интегралы вида

Ч

2,..., (7V-1)]. (5)

0,5

Ввиду того, что функция х имеет нули при С = 0 и 1, для вычисления (5) строились специальные квадратурные формулы. Результаты интегрирования приведены на фиг. 3, где изображены профили скорости и давления (точки) вместе с точным решением

Фиг. 1

6—Ученые записки ЦАГИ № 6

Фиг. 2

81

(сплошная линия), соответствующим случаю Рг = 0,75, <о=1/2,

Моо =10. ,

4. Рассмотрим теперь осесимметричное течение, в котором при больших числах Ие появляется пограничный слой. Этот случай удобен тем, что полученные результаты можно сравнить с известными данными для невязкого течения и пограничного слоя.

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

к поверхности. Будем искать численное решение полных стационарных уравнений Навье —Стокса в полосе 0 <;га-<гаг («). 0 < в < в0, где гаг (в)—фронт отошедшей ударной волны, которую будем считать поверхностью разрыва. В качестве граничных условий используем условия прилипания и условие для температуры при га = О, условия симметрии При 5 = 0 и соотношения Гюгонио при га = гаг (в). Строго говоря, условия ;/на ударной волне являются недостаточными в случае вязкого газа. Однако при записи разностной схемы будем всегда аппроксимировать уравнения на линии га = гаг с помощью граничных и внутренних узлов области, формально получая замкнутую систему. Основания для этого следующие: во-первых, для наших целей неточности в граничных условиях являются несущественными; во-вторых, при больших числах Ие течение около волны практически является невязким, так что способ аппроксимации членов с коэффициентом вязкости в окрестности линии га = гаг не должен сильно влиять на получаемое решение. На линии 5 = «0 также не будем ставить дополнительных условий; при достаточно большой величине Я,) возмущения, распространяющиеся вверх по потоку от линии «=50> должны быстро затухать.

Введем функцию С следующим образом:

О

0,25 рЛ»,и Фиг. 3

с =

(и/иг) + Сп

1 Спг

(6)

где С — постоянная, которая выбирается достаточно большой, чтобы гарантировать монотонность функции С (га).

: Надлежащий выбор величины С позволяет применять преобра-

зование (5) для широкого класса задач (течение в срывной зоне, в следе, в сопле и т. д.).

В системе Навье—Стокса, записанной в координатах п [7], перейдем к независимым переменным 5, С. Для численного решения полученных уравнений можно использовать известные разностные схемы. Однако расчеты проводились с помощью специальной схемы, имеющей повышенный порядок аппроксимации по С; опишем эту схему в общих чертах.

Уравнения для нормальной составляющей импульса, энергии и неразрывности представим соответственно в виде

Т---1

pv(v-Gu)-j-J----p/z —

б) Ц7(рЯ)+^-[р£(^-Си)-^]=^ 7V,

в) Щ 4- p (® - Gh) = 0.

(7)

Здесь И/' — дифференциальный оператор, содержащий только конвективные члены уравнений и не содержащий производных по С; £ — кривизна контура; Н = 1 + Ля;

Е = h +

и1 + г»2

G

дп (s, С )/ds И

хпп и Яп~ соответственно компонент тензора напряжений и нормальная составляющая теплового потока; Ие — число Рейнольдса, вычисленное по параметрам невозмущенного потока, вязкости при А=1 й радиусу кривизны контура в критической точке /?0: все линейные размеры отнесены к /?0, кривизна — к величине 1//?0; Тх и Т2 — все остальные члены уравнений, содержащие коэффициент вязкости. _

На разностной сетке = г‘Д$, С, = /ДС (г = 0, 1... 50/Д5, / = — 0, 1, ..., 1./ДС) введем операторы А+ и Л_, действующие по формулам

И+/)у+1/2 = -у?г 8/‘7+1 —/‘7+2) ;

М-/)„

12

(5/17 + 1 + 8fij — /,7-і).

(8)

Разностную схему для уравнения (7, в) запишем в виде

1(Р®)у+1 - (Р®М/ДС + И^р)г/+і/2 = 0;

v = v — Ом,

(9)

аппроксимировав производные по s, входящие в оператор, по трехточечной схеме с использованием значений функций на предыдущих лучах st = const (t >-2). Легко убедиться, что погрешность аппроксимации в точке (i, j + 1/2) имеет вид

-w(#+4+1/!+°<№+4s,>-

так что схема (9) в точке (i, j+ 1/2) аппроксимирует уравнение (8) на его решении с точностью О (ДС3 + Дя2).

ku2

Применяя операторы А+ и Л_ к выражениям Wpv—Pjjy > W (?Е),

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

части (7, а, б) со вторым порядком точности по всем переменным, будет иметь вид шах [О (Д'3 + Ля'-), О (Д£2/Ие -(- Дз2)]. Очевидно, что в случае больших чисел Ие аппроксимирующие свойства рассматриваемой схемы не ухудшаются при переходе к уравнениям (7, а, б) по крайней мере в области небольших градиентов искомых функций.

Заметив, что операторы А+ и Л_ определяют квадратурные формулы вида

: сУ+1

: | /Л = (Лт/)у+1/2 ДС Ь0(ДС4),

V

легко убедиться, что при соответствующей аппроксимации членов, содержащих коэффициент вязкости, схемы для уравнений (7, а, б) являются консервативными и могут быть получены с помощью законов сохранения для ячеек вида 5/_1/2 < 5 <5г+1/2, С/ < С < Су+ь Выбор операторов А+ или Л_ определялся из условий того, чтобы, во-первых, используемые узлы не выходили за границу области

и, во-вторых, чтобы в линейном приближении все собственные числа соответствующих операторов перехода по модулю не превосходили единицы. Выполнение второго условия было связано с величинами и направлением скоростей. _ .

Уравнение импульса в проекции на ось 5 после некоторых преобразований и замены производных конечными разностями записывалось для луча 5 = в следующем виде:

= аи + -£- и ~ 0, 1, ... , 1/ДС), (10)

где величины а.у и Ьг) не содержат значений функции у на 1-м луче, а (ду-х/д'^ц — разностная аппроксимация производной дрх/д С. К уравнению вида (11) можно прийти, используя законы сохранения касательного импульса для элементарной ячейки.

Уравнение (10) для функции у играет принципиальную роль в полученной системе. В частности, можно показать, что в зависимости от знаков коэффициентов а^ и Ьг;- функция у может соответствовать либо безотрывным, либо отрывным течениям.

Общий ход решения системы алгебраических уравнений итерационным методом состоял в следующем. В течение одной итерации происходил последовательный переход от г-го луча к (г-|-1)-у по неявной схеме, причем в случае г >2 использовались найденные значения на двух предыдущих лучах. Значения функций при 5>-5,+2, а также некоторые другие члены уравнений (7, а—в; брались из предыдущей итерации. Специфика решения уравнений вдоль каждого луча состояла в том, что скорость определялась не из уравнения (7, а), а из уравнения неразрывности (7, в); наоборот, плотность ргу- находилась из уравнения для нормального импульса (7, а)

Для определения кц из уравнения (7,6) использовался метод прогонки. Решение vpaвнeния ; 10) после замены производной в левой части двухточечной схемой сводилось к последовательному вычислению значений /. в точках у=1, 2..., 1/ДС. Граничные значения Хго определялись из условия, что скорость ■Ум/дс, найденная

из уравнений (9), равна скорости на скачке vT- Процесс определения величины ул о осуществлялся с помощью метода Ньютона, причем требуемая сходимость достигалась в результате одной-двух итераций.

Найденные значения Хц на луче s = si(i = 0,l, ... , s0/ks) использовались для определения координаты пузлов в физической плоскости, затем из конечных соотношений вида (6) вычислялись значения касательной скорости tfy.

После завершения очередной итерации, т. е. определения параметров на всех лучах $ = 5г(г = 0, 1,..., s0/As), вычислялись углы наклона ударной волны и, следовательно, значения всех функций на волне для следующей итерации. В окончательном виде любая из полученных величин f™ записывалась в виде

7” = Vv + (l-b)/v_1,

где т — номер итерации; X — релаксационный параметр. Даже без оптимального выбора параметра X (во всех расчетах полагалось Х^;0,1) сходимость до трех-четырех знаков, как правило, достигалась после 200—300 итераций при задании нулевого приближения с помощью постоянных и линейных функций С.

5. Остановимся на результатах расчетов, проводившихся при числах Моз = 10, 7 = 1,4 и Рг = 0,72 для случаев охлажденной сферы

[^^=^00=- 1/(т — 1)М|о и со=0,5]. При числе узлов сетки yV=-^<50,

S —

N[ = ~-^8 величина s0<l,2 варьировалась и не влияла заметным

ъЛ о

образом на приводимые данные.

На фиг. 3 представлены профили скорости и и v, энтальпии h, давления р и плотности р на луче s = 0,6 при Re—100, 500 и 5000 (соответствующие числа Рейнольдса, вычисленные по параметрам невозмущенного потока, приблизительно равны соответственно 640, 3200 и 32000). Для случаев Re = 500 и 5000 на фиг. 4 приведены данные, полученные при N—2b (точки) и N — 50 (сплошные линии); эти результаты почти не различаются. По-видимому, достаточная точность расчетов может достигаться с относительно большим шагом ДС разностной сетки. Отметим также, что при увеличении числа Re зависимости газодинамических параметров от координаты С остаются сравнительно консервативными и не имеют резко выраженных градиентов. Основные изменения, происходящие но мере утоныпения пограничного слоя в физической плоскости, сосредоточены в функции x(s, характеризующей завихренность потока; близкая к постоянной во „внешней области", эта функция начинает быстро возрастать, принимая на поверхности тела тем большие значения, чем больше число Re. Однако ширина области, соответствующей пограничному слою, в плоскости s, С остается конечной и не стремится к нулю при Re -* оо .

_ Изменения газодинамических параметров вдоль того же луча s = 0,6 в физической плоскости приведены на фиг. 4. Там же для сравнения изображены результаты расчетов невязкого обтекания [8] (Мсо = 10, 7=1,4), соответствующие тому же сечению s = const.

Как видно на фиг. 4, профили газодинамических функций, полученные с помощью уравнений Навье—Стокса, везде, кроме пристеночной области, близки к профилям течения невязкого газа; основные изменения этих функций происходят в узкой зоне около

0,2* Р&Ь>а

Фиг. 4

1.0

0,5

1-Яе~100

2-Не=500

< /Ум т

— —неузкое N

Ч

< ~ч

N

/ 2.

/

10

0,5 Фиг. 5

5 II?

0,2

0,1

0.05

0025

0,0125

\

<

\ О

\

\

\ с»

,*Л. \

*0 \

\

сл (лограяачмь/й у слой, о)-0,5) \

Ч

ч ч

101 яР Ж? 10*

Фиг. 6

стенки, где, например, плотность изменяется приблизительна в 20 раз.

На фиг. 5 приведены распределения вдоль поверхности тела коэффициента трения Су = (р/Не) (х^и/Л)|:=о/рсо I/». теплового потока ?«г=(р/КеРг)(хдА/дС)|с=о/роо Уж, давления рхг = р{$, 0), а также величины отхода ударной волны г — пт. _

Там же пунктирными линиями указана зависимость и

в = пт(э) для случая идеального обтекания. Меньшее значение расстояния отхода ударной волны для вязкого газа соответствует

отрицательной величине толщины вытеснения пограничного слоя в рассматриваемом случае сильноохлажденной поверхности (hw — = /too = 0,025). Наконец, на фиг. 6 в логарифмическом масштабе приведено изменение производной (dcf/ds)j„о и коэффициента теплопередачи ch = qw(0)j(ho — hw), где h0 — энтальпия торможения, с изменением числа Re. Линейный характер этих функций при Re>103 и угол наклона прямых указывают на обратно пропорциональную зависимость коэффициентов трения и теплопередачи от величины l/Re. Для сравнения на фиг. 6 изображено изменение коэффициента ch с числом Re согласно теории пограничного слоя [9] при ш = 0,5. Таким образом, результаты расчетов, проведенных на основе полных уравнений Навье—Стокса, при больших числах Re хорошо согласуются с данными для соответствующих невязкого течения и пограничного слоя.

В заключение отметим, что преобразование (1) можно построить и для нерассматриваемого здесь случая, когда целесообразно „растягивать" области как пограничного слоя, так и ударной волны.

ЛИТЕРАТУРА

1. Браиловская И. Ю., Кускова Т. В., Чудов Л. А. Разностные методы решения уравнения Навье—Стокса (Обзор). Вычислительные методы и программирование. Изд. МГУ, вып. XI, 1968.

2. П а в л о в Б. М. Методы теоретического исследования сверхзвукового обтекания затупленных тел потоком вязкого газа. Вычислительные методы и программирование. Изд. МГУ, вып. XV, 1970.

3. Молодцов В. К. Численный расчет сверхзвукового обтекания сферы потоком вязкого совершенного газа. Журн. вычислит, матем. и матем. физ., вып. 9, № 5, 1969.

4. А 11 е n J. S., Cheng S. I. Numerical solution of the near wake of supersonic flow with boundary layer over a sharp corner. The Phys. of Fluids, 1970, v. 13, No 1.

5. Crocco L. A suggestion for the numerical solution of the steady Navier—Stokes equations, A1AA Journal, v. 3, No 10, 1965.

6. M о r e 11 i G., Salas M. D. The blunt body problem for a viscous rarefied gas flow, A1AA Paper, No 69—139, 1969.

7. High speed aerodynamics and jet propulsion, Princeton-New Jersey,

1954, v. 3, Fundamentals of gas dynamics, 1958.

8. Белоцерковский О. М. Расчет обтекания осесимметричных тел с отошедшей ударной волной (расчетные формулы и таблицы полей течения). М., Вычислительный центр АН СССР, 1961.

9. F а у J. A., Riddel F. R. Theory of stagnation point heat transfer on dissociation air, JAS, v. 25, No 2, 1958.

Рукопись поступила 23/11 1972 г.