close

Вход

Забыли?

вход по аккаунту

?

Явная конечноразностная схема для решения задачи обтекания методом установления.

код для вставкиСкачать
УЧЕНЫЕ
т о ом
удк
ЗАПИСКИ
ЦАГИ
.мб
1972
l//
533.011.05
ЯВНАЯ КОНЕЧНОРАЗНОСТНАЯ СХЕМА ДЛЯ РЕШЕНИЯ ЗАДАЧИ
ОБТЕКАНИЯ МЕТОДОМ У~ТАНОВЛЕНИЯ
А. Н. Мuнайлос
Предложена двухшаговая схема второго порядка точности по
пространственным переменным. Схема отличается маJlЫМ количеством
арифметических операций при расчете параметров в одном узле
сетки. Она использована ДJlЯ систематических расчетов двумерных
течений.
1. Явные схемы блаrодаря своей простоте широко используются для реше­
ния газодинамических задач. С увеJlичением числа независимых переменных
возрастает объем вычислительной работы, приходящейся на расчет одного УЗJlа
сетки, а также число самих узлов. Это приводит К значительным затратам ма­
шинного времени, поэтому существенной становится задача создания эффектив­
ной конечноразностной схемы. Такая схема должна быть простой, удовлетворяя
при этом требованиям аппроксимации и устойчивости.
'
Одной из простейших является известная схема П. д. Лакса. Ее обобщения,
ослабляющие влияние
рены в работах
в
[1]
и
искусственной вязкости с помощью параметра а, рассмот­
[2J.
Обозначим
такую схему
L [JJ.
Недостатки ее состоят
следующем.
При большом коэффициенте вязкости (а близко нулю, схема близка схеме
П. д. Лакса) решение быстро устанавливается, но содержит большую погреш­
ность, обусловленную наличием вязких членов; при малом коэффициенте вязко­
сти (а близко единице) решение становится неустоЙчивым. При этом рассматри­
вается линейная зависимость между шагами сетки 't и h по временной и прос­
транственным переменным, поскольку при зависимости 't '"-
h 2,
дающей устойчи­
вый счет при а = J, шаг по времени слишком мал и схема экономически невы'­
годна. Счет в работе [1] с шагом 't ' " h при а = 1 был устойчив на сетке с малым
числом ячеек (6Х2) при стабилизирующем влиянии ударной волны. Другая рас­
смотренная в работе [1] схема L - W принадлежит к классу, исследованному
П. д. Лаксом и, Б. Вендровом. Она применялась в работе [3] и дает более точ­
ные результаты, чем схема L. Колебания фронта ударной волны в процессе
установления при применении этой схемы имеют небольшую
довательно, меньше возмущается поле в расчетной области
амплитуду, сле­
и быстрее идет
процесс установления. Однако при расчете параметров в каждом узле сетки по
этой схеме требуется значительно большее число операций, чем по схеме L,
так как вторые производные по времени определяются путем дифференцирова­
ния системы уравнений по всем независимым переменным. Громоздкость выра­
жений возрастает при расчете пространственных течений.
Настоящая статья посвящена выбору такой разностной схемы, которая,
сохранив преимущества схемы
была бы менее громоздкой, а число ариф­
метических операций на каждом шаге в ней было бы существенно меньшим.
L - W,
2. Предлагаемая схема N представляет собой упрощенный Rариант схемы
из работы [4]. Схема в работе [4] - обозначим ее R - представляет собой моди­
фикацию схемы Р. д. Ри хтмаЙера. Это двухшаговая схема типа prcdlctor-corrector
второго порядка точности.
119
Схема N отличается ОТ' схемы R тем, что на втором шаге (corrector) про­
странственная производн'ая с промежуtочного слоя берется только по одному
н:аправлению,
а
по другим
используются
производн ые
со
старого слоя.
Опишем схему N для случая двух пространственных
Систему газодинамических уравнений представим в виде
и;
s, n.
переменных
р(и, и~, и~);
=
(1)
вектор искомых функций, F - оператор правых частей системы, Р,. - разно­
стный оператор, аппроксимирующий оператор f< с порядком аппроксимации 2.
Как и в схеме R, расчет проводится в два этапа.
U-
Сначала рассчитываются значения искомых функций на вспомогательном
Сllое, совпадающем со слоем t
'С. В отличие от схеме R эти значения опреде­
ляются только в двух точках, лежащих на одной координатной линии с цент­
+
ральной
точкой
(i, Л. Таким образом, схема анизотропна, поэтому она теряет
общность в том смысле, что требует для своего применения выбора определен­
ного направления, т. е. некоторых знаний о решении задачи. В задаче обтека­
ния тупого тела точки на промежуточном слое берутся на координатной линии
s = const, т. е. на некотором луче, направленном от тела к ударной волне.
Именно в этом направлении, как показывает опыт, происходят основные коле­
бания решения в процессе установления. Правильность такого выбора под­
тверждается практическими расчетами. Значения во вспомогательных точках
и, j
1/2) и (i, j - 1/2) определяются по схеме L при а;
О:
+
=
И i ,1+ U i,1-1
Иl , 1-1/2 =
2
(2)
иl,1-1/2 = И l , 1-1/2 + ~Р,. (И, И;, И~)i.
Индекс
t
в выражениях опущен.
Формулы для точки
единицы к индексу
t
Затем
определяются
промежуточном
t
+ 'С.
и,
j
+ 1/2)
значения
получаются
Р,.
в
точке
для определения значений на новом слое
иц = т.
; + 'с [~P,. (И.
из
формул
(i, J)
t + 't
U
прибавлением
.старом· слое
и;, И:)l, 1 + (1 - ~)P,. (и. и~. И~)i. J] ;
_ И l • Нl/2 + Ui, j-l/2
[, j -
на
(2)
2
Ui • 1+1/2 - U
(И--)'
1,1 ш=
h
j•
n
t
и на
используется формула
здесь
120
j-l/2 .
.
1-1/2
(3)
в схеме L нужно Обращатьси в блок pac'leTa оператора F h три раза
Тс учетом значений rJi , J-l/2' сохранившихся от расчета предыдущей точки
(/, j - 1)]. В варианте с ~ = О обращаться к блоку расчета оператора Fh в одной
точке
нужно
два
раза.
Ниже приведена сравнительная
таблица
объема вычислительной
работы
в одном узле сетки для различных схем (за единицу принято количество опера­
ций при однократном обращении в блок расчета оператора
Fh):
Таблица
L
I L-W I
R
N
6-8
4-5
2-З
12-14
6-7
2-З
две пространственные производиые
Три
пространственные про извод-
ные
3. Анализ упрощенной системы показывает, что при 't - - h схема N неус­
,0Йчива. Неустоjiчивость вызвана членами, определиющими производные по
направлению s на слое t. для стабилизации решения используем процесс сгла­
живании, эквивалентный введению искусственной вязкости. Отметим, что lIa
шаге predictor метод
П. д. Лакса.
мы
уже
содержит
члены с искусственной
вязкостью в схеме
Рассмотрим влияние различных законов сглаживания на устойчивость схе­
L для простого моделыIгоo уравнении:
ut+au.l:=O (а>О).
Схему
~
представим в виде ит
_
= u -
а"
(4)
2h (ит +! -
и т -l) .
Результаты аllализа спектральным методом даны в табл.
2:
Таблица
а
Сг лажен.ное
1.
2
-:Г
4
5
II
1
I
Условия
устойчивости
h
2"" (ит+l + Um _ 1)
"<а-
"3 (и т + 1 + llm + Um_ 1)
'C<2V 2
h
+ 4и т +llт _ 1 )
,,< -VS
За
h
Ю(U т + 1 +8u m +u m _ 1)
Ц;;:5ti h
О
"3
значение
1
1
6
(ит + 1
1
ит
2
За
З
Неустойчива
с ослаблением сглаживания уменьшаются предельные значения 'С.
Пример во второй' строке таб./l. 2 соответствует сглаживанию по методу
наименьших квадратов при задании линейного сглаживающего ПОЛИlIома, пост­
роенного на трех узлах. Нетрудно убедиться, что все 9ТИ ЛИllейные виды сгла­
живания сохраняют значения, лежащне на прямой линии, и переводят кривую
второго порядка в кривую второго порядка со сдвигом, пропорциональным h'l.
121
Действительно, формула
_
1
U=
1
+ " т + Um_ 1) = " т + 3
(и т + 1
3
переводит линию тh саму в себя, а кривую
Наличие
к кривой
тому
такого сдвига
второго
было
при
[(ит+l - " т )
сглаживании
h2
h2 + 3
в т2
m 2 h2
+ (и т _ 1 -
и т )]
.
формы ударной волны (близкой
порядка) приводит к большому искажению результатов, поэ­
использовано
сглаживание, основанное на
интерполяции
с помощы()
квадратичного полинома. В этом случае метод наименьших квадратов при четы­
рех используемых узлах дает для уравнения (5) формулу сглаживания:
<
Схема устойчива при 't
0,288 h/a. Сглаживание такого типа и было ис­
пользовано для повышения устойчивости схемы N. В случае трех прострапст­
венных переменных сглаживание проводится отдельно в каждой меридиональ­
ной полу плоскости по формуле
_
3
20 {х [и,. Н2 -
"/, j = Uj,j -
+ (1- х) [и/+2, j
где
Ui,j_l -
" 1_ 1, г-
-
3 (ui, j+l - ui , j)]
3 (и1+1,
j -
+
uu)I},
0';;;:1..<1.
СглаЖИВJlНие применяется не на каждом шаге 't, а через некоторое число
шагов
1. При этом условие УСТОЙ'lивости становится более жестким [для
k>
схемы
k
=
(5)
при
Значения
10 -+- 20, Х
4.
Схема
k = 2,
't
0,25
< -аh].
и Х для отдельных вариантов выбираются
k
= 0,5.
N
была использована
при
создании
эмпирически. Обычн(}
совместно с А. П. Базжиным
и С. В. Пироговой программы расчета обтекания затупленного тела сверхзвуко­
вым потоком [5]. Параметры на границах рассчитывались методом характерис­
тик, близким к схеме д. Моретти
тывались по схеме N.
[3].
Внутренние узлы области решения рассчи­
ПО этой программе проведена большая серия 'расчетов двумерных течений
в совершенном газе (с постоянным отношением удельных теплоемкостей) и в рав­
новеснодиссоциирующем
вращения
кругового
и
воздухе. Рассчитывалось
эллиптических
образующей в виде степенного
обтекание сфер, эллипсоидов
цилиндров,
тел
вращения
с
уравнением
одночлена.
Опыт эксплуатации программы [5]-[7] подтвердил ее эффективность при
удовлетворительной точности (ошибки около 2-4%). Ниже представлены неко­
торые
результаты.
На
фиг.
ношением
1
осей
показана
а/Ь,
скорость на поверхности
эллипсоидов
вращения с от­
изменяющимся от нуля (плоский торец) до двух. Скорость
отнесена к предельной скорости и построена в зависимости от длины дуги
эллипсоида, в качестве единицы линейного размера взят радиус миделя эллип­
соида (полуось Ь). Кривая в случае а
О получена методом интегральных соот­
ношений.
На фиг. 2 сопоставлены картины ударных волн, линий тока и линий посто­
=
янных значений числа М при обтекании эллиптических цилиндров с отношением
осей эллипса, равным единице (круговой цилиндр) и полутора. Результаты соот­
ветствуют значениям % = 1,4 и М ro
1,5. Расчет методом установления при ма-
=
лых
сверхзвуковых
(ослабевают
результатов
ших
числах М СО затруднен из-за плохого установления решения
возмущения в поле течения). Кроме этого, известно, что точность
расчета плоских течений ниже, чем осесимметричных, из-за боль­
расстояний от поверхности
тела до ударной
волны,
поэтому фиг.
2
иллю­
стрирует результаты " трудных· расчетов. Отклонение значений интеграла Бер­
нулли от точного значения не превышает для этих результатов 1,5%. Однако
для малых чисел М СО этот интеграл обычно вычисляется с высокой точностью,
и
точность
его
вычисления
не
является
достаточным
условием
точности
резуль­
татов. Ошибки в расходе не превосходят 3%, а скорость движения ударной
волны в конце счета, т. е. степень "установления· варианта, не больше 0,02 Vmax'
122
IZ=(JJJ
11
(.1'- ~7.f - ~/- ~1.1' ~
-А
11
/
~ .......
I Р j ~~ ~
1/ / ~ 1="""
I ~
)1fМ
~~
..--j ~
4'=11
/. 'l /,
,/,
'/ / /;
v. '/ /
/,
V/ L
/ /
L/
/// / / /
.yj '/ / /'
IM 'l' /.: , /
~~ ~ ....... .....
~~
q
-
~
U-rr=
:1 =
~
~
r-
1
-
~
f
t---j(-
i--
x"'I.II-,М_"'i
:,/
~ ~,
I
Фиг.
Именно
критерий
при
малых
-
.r
I
счета
-
точности по скорости
сверхзвуковых
временной переменной в
4-5
1
движения волны требует
скоростях,
раз. В работе
увеличивая
[8],
длительного
количество
шагов
по
результаты которой изображены
на фиг. 2 штриховой линией, !нот критерий, по-видимому, не рассматривался,
течение при М СО = 1,5 ие установилось, полученные результаты содержат значи­
тельные ошибки. Для
сопоставления
достаточно рассмотреть УГОЛ'" между линией тока и звуковой линией в
точке на ударной волне (см. фиг.
)( = ~ ""; м~ ~.f
Полученное
аналитически
2).
точное
значение", ::::: 9Б • результат настоя­
щей работы", = 89°, результат ра­
боты [8] - примерно БОа. Для по­
вышения
О
точности
значения
"',
по­
лученного в настоящей работе, не­
обходим более длительный счет и,
возможно, более мелкое разбиение
расчетной области в направлении,
нормальном
к
телу.
r
p,S
Фиг.
2
Фиг,
z
3
123
В. И. Благосклоновым
проведены
расчеты обтекания
разующей в виде степенного одночлена
тел
r=zm (1/2<.m<. 1).
вращения
На фиг.
поле течения для варианта со значениями т = 0,65; 'Х, = 1,4; М""
= 8.
с об­
показано
3
Характерна
~ольшая по сравнению с более тупыми телами длина линий М = сопs!.
На
J(
-= 1.1f "
М"",,= Q
фиг.
изображена
4
форма
ударной
волны, изомахи и изохоры для тела вращения,
образующая которого задана формулой
-
r=
+
1+
.
Здесь
tg а: V
-=-
( ;)n
1+
~/
1
()n
2
Z
2
-
а1
+
Zl
(+)n
(а2 Z2
+ bz Z + С2) •
1
R2
2
С2= 1 - Т Ь2'
Расчеты течения у тел такой формы про­
вел А. П. Косых. В 9ТИХ расчетах исполь­
зовано
неравномерное разбиение расчетной
области в физической плоскости в направле­
нии вдоль образующей тела (вдоль оси х).
Варианту, представленному
на
фиг. 4,
соответствуют
следующие
ров: для формы тела
=
Фиг.
n = 100; zl 0,5045;
-'Х, = 1,4; М"" = 6.
4
- R1
для
значения
парамет­
= 0,6; R2 = 0,2; а: = 600;
набегающего
потока
ЛИТЕРАТУРА
1.
Косых
А.
П.,
Минайлос А.
Н. О явных схемах метода
установления в задаче сверхзвукового обтекания затупленного тела.
Журн. вычисл. матем. и матем. физ., т. 10, J!ё 2, 1970.
2. V а п L е е r В. stаыliаtiопп of difference schemes for the Equations of lnviscid compressibIe f(ow Ьу arllficial diflusion. J. Compu! Phys.,
У. 3, No 4, 1969.
3. М о r е t t У О., А Ь Ь е t М. А timedependent computatlonal method for bIunt body flows. AIAA J., У. 4. No 12, 1966.
4. В u r s t е i n S. Z. Finile-dlfference caJculations for hydrodinamic
flows сопtаtпlпg dlscontlnuities. J. Comput., Phys., У. 2, No 2, 1967.
5. Б а з ж и н А. П., Б л а г о с к л о н о в В. Н., М и н а й л о с А. Н.,
Пир о г о в а С. В. Обтекание сферы сверхзвуковым потоком совер­
шенного газа .• Ученые записки ЦАГИ', т. 11, N2 3, 1971.
6.
К о с ы х
поверхности
А. П., М и н а й л о с
сверхзвуковым
воздуха .• Ученые записки ЦАГИ', т.
7.
М и н а й л о с
зависимости
А.
потоком
ll,
А. Н. Параметры
осесимметричного
Н. Обтекание сферической
равновеснодиссоциирующего
М
5, 1971.
подобия
сверхзвукового
и
корреляционные
течения
у 9ЛЛИПСО­
ида. Известия АН СССР, МЖГ,
1973.
8. М о r е t I О. Three-dimensional Invlscid flow about supersonic
blunt cones а! angle о! atlack. Part 2, SC - CR - 68 - 3728, 1968.
Рукопись поступила
5// 1972
г.
Документ
Категория
Без категории
Просмотров
9
Размер файла
272 Кб
Теги
решение, методов, обтекании, конечноразностная, явная, задачи, схема, установление
1/--страниц
Пожаловаться на содержимое документа