Явная конечноразностная схема для решения задачи обтекания методом установления Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ Ц АГ И То м III 197 2

№ 6

УДК 533.011.05

ЯВНАЯ КОНЕЧНОРАЗНОСТНАЯ СХЕМА ДЛЯ РЕШЕНИЯ ЗАДАЧИ ОБТЕКАНИЯ МЕТОДОМ УСТАНОВЛЕНИЯ

А. Н. Минайлос

Предложена двухшаговая схема второго порядка точности по пространственным переменным. Схема отличается малым количеством арифметических операций при расчете параметров в одном узле сетки. Она использована для систематических расчетов двумерных течений.

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

Одной из простейших является известная схема П. Д. Лакса. Ее обобщения, ослабляющие влияние искусственной вязкости с помощью параметра а, рассмотрены в работах [1] и [2]. Обозначим такую схему Ь [1]. Недостатки ее состоят в следующем.

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

числом ячеек (6X2) при стабилизирующем влиянии ударной волны. Другая рассмотренная в работе [1] схема Ь—принадлежит к классу, исследованному П. Д. Лаксом и Б. Вендровом. Она применялась в работе [3] и дает более точные результаты, чем схема Ь. Колебания фронта ударной волны в процессе установления при применении этой схемы имеют небольшую амплитуду, следовательно, меньше возмущается поле в расчетной области и быстрее идет процесс установления. Однако при расчете параметров в каждом узле сетки по этой схеме требуется значительно большее число операций, чем по схеме Ь, так как вторые производные по времени определяются путем дифференцирования системы уравнений по всем независимым переменным. Громоздкость выражений возрастает при расчете пространственных течений.

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

2. Предлагаемая схема N представляет собой упрощенный вариант схемы из работы [4]. Схема в работе [4] — обозначим ее Я — представляет собой модификацию схемы Р. Д. Рихтмайера. Это двухшаговая схема типа рге<11с(ог-соггес-1от второго порядка точности.

Схема N отличаетсй or схемы R тем, что на втором шаге (corrector) пространственная производная с промежуточного слоя берется только по одному направлению, а по другим используются производные со старого слоя.

Опишем схему N для случая двух пространственных переменных s, п. Систему газодинамических уравнений представим в виде

U't = F(U, U's, U'ny,

0)

U — вектор искомых функций, F — оператор правых частей системы, F/,-разностный оператор, аппроксимирующий оператор t с порядком аппроксимации 2. Как и в схеме R, расчет проводится в два этапа.

Сначала рассчитываются значения искомых функций на вспомогательном сдое, совпадающем со слоем В отличие от схеме R эти значения опреде-

ляются только в двух точках, лежащих на одной координатной линии с центральной точкой (г, j). Такйм образом, схема анизотропна, поэтому она теряет общность в том смысле, что требует для своего применения выбора определенного направления, т. е. некоторых знаний о решении задачи. В задаче обтекания тупого тела точки на промежуточном слое берутся на координатной линии s = const, т. е. на некотором луче, направленном от тела к ударной волне. Именно в этом направлении, как показывает опыт, происходят основные колебания решения в процессе установления. Правильность такого выбора подтверждается практическими расчетами. Значения во вспомогательных точках (г, у -J- 1/2) и {I, j — 1/2) определяются по схеме L при а — 0:

U

I. У—1/2'

_ Ui.j + UUj-1 .

/—1/2'і

Ч' _ W + l..J-U l-l.fl + ^Ч-W-l ~ Ul~UJ~l) = -------------------------------------------

4 Л,

Wl,J-\ldn ~~

U\j-1/2

: Ui, у—1/2 + tFh (U > Vs. U'Xi-

r—1/2 •

(2)

Индекс < в выражениях опущен.

Формулы для точки (/, /-(-1/2) получаются из формул (2) прибавлением единицы к индексу /.

Затем определяются значения Р/, в точке («, ]) на .старом" слое £ и на промежуточном < + т.

Для определения значений на новом слое < + г используется формула

Ulj= Ui, і + т [pf* (U, U's, U'n)i, j + (1 - ?)Fa(U, U's. Un)i, j) ;

(3)

здесь

Tl - j +1'2 u‘. /-1'2 •

ui, і--------------о---------- ’

at \'—n, V — ^i+iJ ^i-Uj ■

Wi, A - (Ui. л------------2X--------- ’

Ui,j+H2 — Ui, y—1/2 .

(Ui, j)w =

Ui. /4-і — Uі, /-1 2 hh

В схеме /. нужно обращаться в блок расчета оператора Р^ три раза {с учетом значений £/г< у_1/2, сохранившихся от расчета предыдущей точки

У— !)]• В варианте с р = 0 обращаться к блоку расчета оператора Рн в одной точке нужно два раза.

Ниже приведена сравнительная таблица объема вычислительной работы л одном узле сетки для различных схем (за единицу принято количество операций при однократном обращении в блок расчета оператора Р^):

Таблица 1

L L-W R N

Две пространственные производные 1 6-8 4-5 2-3

Три пространственные производные 1 12-14 6-7 2—3

3. Анализ упрощенной системы показывает, что при х — Л схема N неустойчива. Неустойчивость вызвана членами, определяющими производные по направлению 5 на слое t. Для стабилизации решения используем процесс сглаживания, эквивалентный введению искусственной вязкости. Отметим, что на шаге predictor метод уже содержит члены с искусственной вязкостью в схеме П. Д. Лакса.

Рассмотрим влияние различных законов сглаживания на устойчивость схемы L для простого модельного уравнения:

ut + aux — 0 (я>0). (4)

— az

Схему представим в виде ихт = и - (um+i — ит-\) ■

Результаты анализа спектральным методом даны в табл. 2:

Таблица 2

а Сглаженное значение И Условия устойчивости

1 А

0 T(“m+l + «™-l)

1 . 3 1 3 (ит+1 + ит + ит-\) '<*V л О CL

2 3 1 -6"(e«+l + 4“* + “m-l) 3 а

4 5 + 8ия* + мт-і) 3 Т<

1 ит Неустойчива

С ослаблением сглаживания уменьшаются предельные значения т.

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

и — - д- (**„,_!_! + Кщ + Чт-\) = ат + "з- [(“(Я+1 ит) + (мт—1 — ит)]

Л2

переводит линию /пЛ саму в себя, а кривую /н2Л2 в т2^2-)--^--

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

~ 3

и = ит~- ~20 1нт+2 — “т-1 ~ 3 (“т+1 “ ат) 1 •

Схема устойчива при т <; 0,288 Л/я. Сглаживание такого типа и было использовано для повышения устойчивости схемы N. В случае трех пространственных переменных сглаживание проводится отдельно в каждой меридиональной полуплоскости по формуле

~ 3

а1, ] = “/,/ - 20 К /+2 “ «и-1 - 3 (“*. у+1 - “л у)1 +

+ (1 - X) К+2> / - “г-1. ; ~ 3 (“;+1, у - «/, ;)1}.

где 0<х<1*

Сглаживание применяется не на каждом шаге т, а через некоторое число шагов &>1. При этом условие устойчивости становится более жестким [для

0,25

схемы (5) при «=2, т<—а~~

Значения к и х для отдельных вариантов выбираются эмпирически. Обычно к = 10 -ь 20, х = 0,5.

4. Схема N была использована при создании совместно с А. П. Базжиным и С. В. Пироговой программы расчета обтекания затупленного тела сверхзвуковым потоком [5]. Параметры на границах рассчитывались методом характеристик, близким к схеме Д. Моретти [3]. Внутренние узлы области решения рассчитывались по схеме N.

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

Опыт эксплуатации программы [5]—[7] подтвердил ее эффективность при удовлетворительной точности (ошибки около 2—4%). Ниже представлены некоторые результаты.

На фиг. 1 показана скорость на поверхности эллипсоидов вращения с отношением осей я/£>, изменяющимся от нуля (плоский торец) до двух. Скорость отнесена к предельной скорости и построена в зависимости от длины дуги эллипсоида, в качестве единицы линейного размера взят радиус миделя эллипсоида (полуось Ь). Кривая в случае а = 0 получена методом интегральных соотношений.

На фиг. 2 сопоставлены картины ударных волн, линий тока и линий постоянных значений числа М при обтекании эллиптических цилиндров с отношением осей эллипса, равным единице (круговой цилиндр) и полутора. Результаты соответствуют значениям г. = 1,4 и Мпо=1,5. Расчет методом установления при малых сверхзвуковых числах затруднен из-за плохого установления решения (ослабевают возмущения в поле течения). Кроме этого, известно, что точность результатов расчета плоских течений ниже, чем осесимметричных, из-за больших расстояний от поверхности тела до ударной волны, поэтому фиг. 2 иллюстрирует результаты „трудных" расчетов. Отклонение значений интеграла Бернулли от точного значения не превышает для этих результатов 1,5%. Однако для малых чисел Мот этот интеграл обычно вычисляется с высокой точностью, и точность его вычисления не является достаточным условием точности результатов. Ошибки в расходе не превосходят 3%, а скорость движения ударной волны в конце счета, т. е. степень „установления" варианта, не больше 0,02 Ктах-

и

Фиг. 1

Именно критерий точности по скорости движения волны требует длительного счета при малых сверхзвуковых скоростях, увеличивая количество шагов по временной переменной в 4—5 раз. В работе [8], результаты которой изображены на фиг. 2 штриховой линией, этот критерий, по-видимому, не рассматривался, течение при Мсо= 1,5 не установилось, полученные результаты содержат значительные ошибки. Для сопоставления достаточно рассмотреть угол со между линией тока и звуковой линией в точке на ударной волне (см. фиг. 2). Полученное аналитически точное значение ш я: 96°. результат настоящей работы <о = 89°, результат работы [8] — примерно 60°. Для повышения точности значения <о, полученного в настоящей работе, необходим более длительный счет и, возможно, более мелкое разбиение расчетной области в направлении, нормальном к телу.

В. И. Благосклоновым проведены расчеты обтекания тел вращения с образующей в виде степенного одночлена r = zm (1/2<.т<; 1). На фиг. 3 показано поле течения для варианта со значениями т — 0,65; *=1,4; Мот = 8. Характерна большая по сравнению с более тупыми телами длина линий М = const.

На фиг. 4 изображена форма ударной волны, изомахи и изохоры для тела вращения, образующая которого задана формулой

н =£*; Мт=6

Г =

1 +

Ш'

tg а V Z2 —а\ +

+

Здесь

Ф*

(я2 ^ bnZ + с2).

а\

Ri

tg2 а ’

h — -Fr +

а3 — — гI tg а

2 Я,

Фиг. 4

с - 1 - — 62 с2 — I — 2 2 •

Расчеты течения у тел такой формы провел А. П. Косых. В этих расчетах исполь-

зовано неравномерное разбиение расчетной области в физической плоскости в направлении вдоль образующей тела (вдоль оси х).

Варианту, представленному на фиг. 4, соответствуют следующие значения параметров: для формы тела — = 0,6; = 0,2; а = 60°;

п = 100; 21 = 0,5045; для набегающего потока

— % = 1,4; М = 6.

ЛИТЕРАТУРА

1. Косых А. П., Минайлос А. Н. О явных схемах метода установления в задаче сверхзвукового обтекания затупленного тела. Журн. вычисл. матем. и матем. физ., т. 10, № 2, 1970.

2. Van Leer В. Stabilization of difference schemes for the Equations of Inviscid compressible flow by artificial difluslon. J. Comput Phys., v. 3, No 4, 1969.

3. M о r e 11 у G., A b b e t M. A timedependent computational method for blunt body flows. AIAA J., v. 4, No 12, 1966.

4. Bur stein S. Z. Finite-dlfference calculations for hydrodinamic flows containing discontinuities. J. Comput., Phys., v. 2, No 2, 1967.

5. Базжин А. П., Благосклонов В. H„ Минайлос А. Н., Пирогова С. В. Обтекание сферы сверхзвуковым потоком совершенного газа. .Ученые записки ЦАГИ", т. И, № 3, 1971.

6. Косых А. П., М и н а й л о с А Н. Обтекание сферической поверхности сверхзвуковым потоком равновеснодиссоциирующего воздуха. .Ученые записки ЦАГИ", т. II, № 5, 1971.

7. Минайлос А. Н. Параметры подобия и корреляционные зависимости осесимметричного сверхзвукового течения у эллипсоида. Известия АН СССР, МЖГ, 1973.

8. More И G. Three-dimensional Inviscid flow about supersonic blunt cones at angle of attack. Part 2, SC — CR — 68 — 3728, 1968.

Рукопись поступила 5/1 1972 г.