CC BY
31
УЧЕНЫЕ ЗАПИСКИ ЦАГИ Т о м XVI 19 8 5
№ 4
УДК 533.6.011.35
УСКОРЕНИЕ СХОДИМОСТИ МЕТОДОВ РАСЧЕТА ПЛОСКОГО И ПРОСТРАНСТВЕННОГО ТРАНСЗВУКОВОГО ОБТЕКАНИЯ ТЕЛ В НЕОГРАНИЧЕННОМ ПОТОКЕ
Н. А. Владимирова, В. В. Вышинский, С. В. Ляпунов, Я■ М.- Серебрийский
Показано, что основной причиной медленной сходимости конечноразностных методов, использующих отображение внешности обтекаемого тела на конечную область, является наличие бесконечно удаленной точки в расчетной области. С целью ускорения сходимости применен алгоритм «замораживания» дальнего поля при расчетах на сетках после удвоения. Применение алгоритма «замораживания» при расчете плоских и пространственных трансзвуковых течений идеального газа позволяет сократить время вычисления в 4—6 раз.
Успехи в вычислительной трансзвуковой аэродинамике в последние годы в значительной степени связаны с использованием метода релаксации для расчета потенциальных трансзвуковых течений. Предложенный первоначально для решения трансзвукового уравнения малых возмущений в задаче об обтекании профиля [1] метод нашел широкое применение для решения полного уравнения относительно потенциала скорости в задачах об околозвуковом обтекании профилей [2, 3], тел вращения [4, 5], пространственных тел типа фюзеляжа самолета [6] и изолированных крыльев [7, 8].
Несмотря на свою эффективность, методы релаксации характеризуются довольно низкой скоростью сходимости. Задачи оптимизации формы тел, учета вязкости в приближении пограничного слоя, задачи проектирования, опирающиеся на многочисленные параметрические исследования, связаны с многократным расчетом поля течения. Эти задачи стимулируют дальнейшие работы по сокращению времени вычисления.
В последнее время интенсивно разрабатываются методы и приемы, ускоряющие сходимость релаксационных методов. В работе [9] предложена схема экстраполяции решения совместно с методом релаксации от поверхности тела в поле потока, что позволяет ускорить сходимость в 2—4 раза. В работе [10] предложен смешанный метод, использующий прямые способы решения уравнения Пуассона. Эффективными оказались методы факторизации, приближенной факторизации [11, 12] и ме-
тод многих сеток [13], позволяющие в ряде случаев на порядок сократить время расчета околозвукового обтекания. Все эти методы, за исключением первого, требуют фактически разработки новых программ расчета.
В данной работе показано, что скорость сходимости по мере удаления от тела падает. Исключение области дальнего поля из зоны расчета и постановка на ее границе краевого условия Дирихле для потенциала скорости, полученного из расчета на грубой сетке, позволяют существенно ускорить сходимость итерационного процесса. Предлагается схема проведения расчетов обтекания с «замораживанием» потенциала в дальнем поле на окончательной расчетной сетке после последнего ее удвоения. Предлагаемый метод не исключает возможности сочетания его с методами [9—13] для дальнейшего сокращения времени расчета.
1. В работах [2] и [5] для расчета плоского и осесимметричного обтекания используется конформное отображение на круг единичного радиуса внешности профиля или внешности продольного меридионального сечения тела вращения. При этом бесконечно удаленная точка переходит в центр круга. Уравнение для ф, регулярной при г->0 части полного потенциала, в плоскости круга в полярной системе координат г, 0 для плоского ^ = 0 и осесимметричного у=1 случаев имеет вид:
где а — скорость звука; и и V-—компоненты скорости в плоскости круга; хю — модуль Якобиана отображающей функции; х, у— исходные прямоугольные декартовые координаты, где ось х связана с хордой профиля или осью симметрии тела вращения.
При малых скоростях а2^>и2 + и2 уравнение (1) переходит в уравнение Лапласа:
В эллиптической области течения вторые производные аппроксимируются конечно-разностными выражениями
где индексы означают номер узла разностной сетки по координате 0 и г
соответственно, знак Н---величину потенциала на текущей итерации, в
противном случае ср берется на предыдущей итерации, со — параметр релаксации (равный 1,4 в [2] и 1,6 в [5]). Первые производные аппроксимируются значениями ср на предыдущей итерации центральными разностями.
Наличие в (2) члена г2фгг приводит к вырождению его в бесконечно удаленной точке при г->0.
(а2 — и2) <р99 — 2гиг»<р,в + г2 (а2 — V2) — 2иъ<рд +
4- г2 (а- + и2 — 2г;2) <?г + ~ (и2 + V2) (ита0 + гтюг) -г
0)
(2)
В соответствии с методом Неймана анализа устойчивости разностных схем рассматривается поведение Фурье-гармоники
<fkj = eiAB eiBr, }
1 (4)
и определяется коэффициент нарастанция гармоники X. Подстановка (4) в (3), а затем в (2) и пренебрежение квадратичными по А г членами приводит к выражению:
еы-----X — (2-------4- е~иХ + г2 р2 [2г sin рДгХ' — 4sin2 p/2X] -f-
-f- rp0r sin р (1-~хь) + V—Д0г sin а = О,
где а = ЛД0, р = 5Дг, р = -|£-.
В случае sin р Ф 0 при г 0 для определения X получается дифференциальное уравнение вида
г2)/ + GX = О,
где G = (e~ia — 2/<o)/(2i sin рДгр2). Решение его X = e°lr, Re G =
= — sina/(2sin (ЗДгр2)<;0 для всех гармоник О О О, 0 < и при
г 0 они затухают.
При р = тг в случае г-»- 0 несложно получить
, I - <а/2 (2 — еи)
1 — со/2-^—
причем максимум Ке^=1 достигается при а = 0.
В случае а = 0, Р = я, что соответствует высокочастотным осцилляциям вдоль координаты г, для % получается выражение
1
X =
4г2 р2
1 + 2/со — 1
которое при 0 (0<со<2) стремится к 1 и, следовательно, скорость сходимости рассматриваемых возмущений замедляется.
Наличие окрестности точки г = 0 в расчетной области является, таким образом, причиной медленной сходимости релаксационного процесса. С другой стороны, решение в этой области может быть оценено из известного асимптотического поведения потенциала на бесконечности или путем расчета значений потенциала на грубой разностной сетке. В случае расчета обтекания близких тел дальнее поле может быть взято из расчета другого тела при тех же значениях числа Маха набегающего потока и угла атаки. Последующее «замораживание» величины потенциала в этой области и исключение ее из релаксационной процедуры значительно ускоряют сходимость процесса. Ускорение имеет место также и за счет уменьшения порядка разностных уравнений (сокращения количества узлов разностной сетки и, следовательно, времени одной релаксации).
С целью сохранения непрерывности нормальной производной на границе замораживаемой области через 10 релаксаций с «замораживанием» проводится одна релаксация всего поля. Этот прием, опирающийся на метод Шварца [14], гарантирует сходимость метода, во всяком случае при дозвуковых скоростях набегающего потока, как в плоском, так и в пространственном случаях.
Аналогичный прием ускорения сходимости итерационного процесса был применен к задачам расчета потенциального обтекания изолированных фюзеляжа [6] и крыла [8] самолета. В этих задачах методом релаксации решается полное уравнение для потенциала скорости течения. При этом бесконечная область физического пространства около обтекаемого тела отображается во внутренность конечного параллелепипеда. В расчетных координатах |г]£ уравнение для потенциала возмущения ф(|, г), £) записывается в виде:
+ ВЧ-т + + Я “ °- <5>
Коэффициенты при вторых производных и свободный член уравнения (5) являются функциями младшйх производных <р£, <рг, <рс, коэффициентов матрицы преобразования переменных, координат и угла атаки. Уравнение (5) аппроксимируется в узлах равномерной разностной сетки разностным оператором. Аппроксимация старших производных (используется следящая разностная схема) производится в зависимости от типа уравнения: в сверхзвуковых узлах используются гиперболические разностные операторы первого порядка точности, в дозвуковых — центральные разности:
А<2 Ы=?£ у* - + ‘РГ-и*;
Д-']2 [?,,] = <#_!* - 2/<о?+* - (2 - 2/ш) Т|/А + ?и+1к;
№ [?к] = Т|%_1 - 2Н& - (2 - 2/со) + ?./А+1;
4Д? Дт) [<р£т)] = ®.+1;.+1А *Рг4-1у—1* ~ <Р/_1/4-и + сР/+_1/_ш
4Д| ДС = ?(-+1^ + 1 — 'Рг—1 у*4-1
4Дт! ДС [^с] = ср..+1А+1 - <р+_1й+1 - •
Здесь ср+й и Ф,7А, как и выше, значения <р в узлах сетки на текущей и предыдущей итерациях, Д&, Дтг), ДС — размеры ячеек сетки по соответствующим осям. При малых скоростях (и2 + V2 + ®>2<С1а2, где и, V, IV — компоненты скорости) уравнение (5) имеет эллиптический тип и может быть записано в виде
9% + + Н2 <р,, = 0. (7)
Коэффициент Н2 при производной сра равен якобиану преобразования координат и неограниченно возрастает при удалении от тела. Гармонический анализ уравнения (7) позволяет получить оценку скорости затухания высокочастотных гармоник. Подставляя выражения (6) в (7) и представляя ф в виде
= е‘ (“Е+Рч+тО; ф+к = д (с, г\) е1 (“НР’1+тС)}
пренебрегая членами более высокого порядка, чем А£2 +Ат]2+А£2, можно получить выражение для А, характеризующее скорость убывания амплитуды гармоник по координате т):
1 2/ю — 3
1 + № 0з 2/ю — 1
1 + 1 2/со 1
где 0=(Аг)/А^), откуда для параметра релаксации 1<ш<2 и конечных значений Я следует, что А<1, и рассматриваемые гармоники затухают.
При неограниченном возрастании Н, соответствующем удалению от поверхности тела, Я->1, и сходимость, как и в двумерном случае, замедляется.
2. Аппробация алгоритма «замораживания» в случае плоского обтекания приведена на примерах расчетов обтекания профиля КАСА0012 при угле атаки а=Г и того же профиля при заданном коэффициенте подъемной силы су = 0,5. Расчеты проводились для случаев до-критического (число Маха набегающего потока Моо = 0,4) и закритиче-ского (Моо = 0,8) обтекания на сетке с числом узлов после удвоения 80x16 (80 точек на контуре профиля). Круг «замораживания» содержит до восьми слоев разностной сетки по г, что соответствует в физи-
N АС А ООП <х=1°;М„=0,й
ческой плоскости «замораживанию» поля течения от бесконечности до расстояния от профиля не менее величины его хорды. Расчет заканчивается при выполнении 200 итераций на удвоенной сетке, либо когда значение максимальной поправки к потенциалу на итерации Д<р становится меньше 10~в.
На рисунке вверху представлена сходимость итерационного процесса в виде зависимости логарифма максимальной поправки к потенциалу на итерации Дер от номера итерации Митер. Приведенные данные соответствуют расчету обтекания профиля ЙАСА 0012 при угле атаки а=1° и числе М набегающего потока Моо = 0,8. На верхней поверхности профиля сверхзвуковая зона заканчивается сильным скачком уплотнения. Цифры 0, 4 и 8 на рисунке означают количество слоев разностной сетки с «замороженным» значением потенциала. Наблюдается заметное ускорение сходимости, которое увеличивается с ростом области «замораживания». На рисунке внизу приведено изменение по итерациям коэффициента подъемной силы су. С увеличением количества «замораживаемых» слоев значение Су быстрее сходится к постоянной величине.
Численный эксперимент позволил выявить ошибки, вносимые в результаты при использовании «замораживания». Приведенные ошибки (табл. 1) представляют собой разности соответствующих величин при расчете без «замораживания» при выполнении 200 итераций и с «замораживанием» восьми слоев удвоенной сетки с прекращением расчета
NACA 0012 а = 1° = 0,5
Мае 0,4 0,8 0,4 0,8
Асх в 0,8-10-4 2,0-10-4 1,3-10-* 1,5-10-4
1,2% 0,6%
0,2-Ю-з 5-10-3
0,2И 2,1%
Да — 0,03° 0,04°
0,8% 1,6%
min 0,002 0,002 0,002 0,002
при достижении условия Дф=10 6. Ошибки в определении сопротивления сх в не превышают 2-10'4, что для рассмотренных случаев закри-тического обтекания составляет около 1%. Ошибки в расчете коэффициента подъемной силы су, угла атаки а и минимального значения коэффициента давления на профиле также невелики.
Алгоритм «замораживания» применен к расчету околозвукового обтекания тел вращения (используется метод [5] на сетке с числом узлов после удвоения 41X30) и пространственных удлиненных тел (используется метод [6] на сетке после удвоения 41X30X13). Поверхности тела принадлежит в осесимметричном случае 41 узел, в пространственном 41X13 узлов (используется 13 продольных меридиональных сечений, связанных с наибольшей хордой поверхности). Производится «замораживание» N слоев дальнего поля в осесимметричном случае по полярной переменной г, в пространственном •— по соответствующей цилиндрической переменной. В последнем случае продольная ось цилиндрической системы координат связана с наибольшей хордой поверхности (более подробно см. [6]). Расчет прекращается при выполнении 200 итераций или при уменьшении поправки к потенциалу на итерации ниже 10"5.
Рассчитывалось обтекание следующих тел, взятых из работ [5, 6]: сферы, незамкнутого полубесконечного тела вращения (тело I) и пространственного тела типа фюзеляжа самолета (тело II). Использовались значения параметра N, равные 0 (контрольный расчет), 5, 10, 15 («замораживание» половины пространства) и 20.
В табл. 2 приведена полученная при расчете обтекания сферы зависимость относительной погрешности в определении максимального значения числа Маха на теле Мтах при расчете с «замораживанием» N слоев по сравнению с контрольными результатами, N—0.
В табл. 3 результаты для сферы представлены в случае N=20 при различных числах М,*,. Относительные погрешности в определении волнового сопротивления схв представлены в табл. 4.
Таким образом, использование N<20 в случае не очень сильных
СКаЧКОВ уПЛОТНеНИЯ (Мтах<1,7) Сохраняет ТОЧНОСТЬ Определения Мщах до 0,2%, Схв —до 1 т-2%.
Таблица 2
\ N 5 10 15 20
Дмтах М ’ •lvimax 0,3 0.6 0,08 0,18 0,09 0,18 0,14 0,18 0,12 0,15
Аналогичные результаты дают расчеты обтекания тела II. В табл. 5 для случая N=15 представлены значения погрешностей в определении Мтах и схв. Для сравнения следует указать, что вообще точность инте-
Таблица 3
Моо 0,3 0,6 0,7 0.8 0.85 0,9
ДМгаах мтах’ /0 0.12 0,15 0,10 1,5 0,21 -1,4
Таблица 4
МсХ. 5 10 15 20
0,7 0,26 0,30 0,40 0,63
^ % 0,8 0.76 0.92 1,18 1,86
сх в 0,85 0,96 1,15 1,49 2.34
0,9 -0,76 —0,85 -1,15 -1,80
Таблица 5
0,6 0,85 0.89
АМтах М * ™ 1Птах —0,12 —0,14 -0,38
&сх в 1 О о 7,7-10-4 -0,5-10-4
% в 2,8 4,4 -0,13
Таблица 6
Таблица 7
N 0 5 10 15 20
м* 0,759 0,758 0,757 0,756 0,757
00 о 1! 8 5 М = 0,85 00
1,7 3,5
Су
Д Сг
г , % *•X 2,9 0
грирования схв (интеграл нормальных сил, отнесенный к скоростному напору и площади миделя) порядка 10~3.
Последовательный расчет обтекания при различных значениях числа Моо позволяет определить величину критического числа М*. В табл. 6 для тела I приведены значения М* при различных значениях параметра N.
Число сверхзвуковых точек в случае замкнутого тела существенно не меняется, если сверхзвуковая зона не попадает в область «замораживания», граница которой расположена от тела на расстоянии порядка его длины при N=15.
Результаты численных экспериментов для тел вращения и пространственных удлиненных тел позволяет выбрать в качестве параметра
«замораживания» в реализованном алгоритме Л/=15, что приводит к сокращению времени вычисления в 5—6 раз в осесимметричном и в 3—4 раза в пространственном случаях по сравнению с обычным расчетом без «замораживания» дальнего поля.
Алгоритм «замораживания» применен к расчету околозвукового обтекания изолированных крыльев. Как и в исходной версии [8], используется двукратное удвоение сетки. Окончательная сетка после удвоений имеет 192X24X32 узла. При расчете на ней вырезается внутренний параллелепипед, включающий 150Х 10x32 узла. В нем осуществляется итерационный процесс. В физическом пространстве выделяемый внутренний параллелепипед соответствует области вблизи крыла величиной порядка его наибольшей хорды. В узлах сетки, расположенных вне выделенной области и определяющих поле вдали от обтекаемого крыла, фиксируются значения потенциала, полученные из расчета на предыдущей сетке. В направлении размаха крыла «замораживание» решения не производится. Применение такого подхода приводит к сокращению времени расчета в 3—5 раз.
Анализ, проведенный для крыла с углом стреловидности % = 30° при обтекании его на околокритическом (Моо = 0,8) и закритическом (Моо = = 0,85) режимах при коэффициенте подъемной силы cy~0,5, показывает, что использование алгоритма «замораживания» с данным выбором параметров практически не приводит к потере точности результатов. В табл. 7 представлены относительные погрешности в определении подъемной силы и сопротивления крыльев при расчете с «замораживанием» по сравнению с результатами, рассчитанными по исходной версии программы [8]. Погрешности в определении коэффициентов давления на поверхности крыльев не превышают при этом 1%.
Таким образом, в численных экспериментах выявлена допустимая область «замораживания» при сохранении приемлемой точности расчета. Применение данного алгоритма не требует существенной перестройки программ расчета и позволяет в 4—6 раз сократить время вычисления при сохранении значений интегральных аэродинамических характеристик с точностью до 1—2%.
В заключение авторы выражают благодарность Ю. Б. Лифшицу за полезное обсуждение данной работы.
ЛИТЕРАТУРА
1. Murman Е. М., С о 1 е J. D. Calculation of plane steady transonic flows. AIAA J., 1971, vol. 9, N 1.
2. Bauer F., Garabedian P., Korn D., Jameson A. Supercritical wing sections II. — Lect. Notes in Economics and Math. Syst., 1975.
3. Лифшиц Ю. Б. К теории трансзвуковых течений около профиля.— Ученые записки ЦАГИ, 1973, т. IV, № 5.
4. Jameson A,, S о u t h J. С. Relaxation solulions for inviscid axisvmmetric transonic flow over blunt or pointed bodies. — AIAA Comp.
Fluid Dyn. Conf., 1973.
5. Вышинский В. В. Расчет околозвукового осесимметричного обтекания тел вращения.—Ученые записки ЦАГИ, 1982, т. XIII, № 5.
6. Вышинский В. В. К расчету пространственного околозвукового обтекания удлиненных тел. — Ж. вычисл. матем. и матем. физ., 1983, т. 23, № 1.
7. Jameson A.,. Caughey D. A. Numerical calculation of the transonic flow past a swept wing. NASA—CR—153297, 1977.
8. Владимирова H. А. Исследование обтекания прямых и стреловидных крыльев большого удлинения при околозвуковых скоростях.— Ученые записки ЦАГИ, 1983, т. XIV, № 4.
3—«Ученые записки» № 4
33
9. Y u N. J., R u b b e r t P. E. Acceleration schemes for transonic potential flow calculations. — AIAA Paper N 80—0338, 1980.
10. J a m e s о n A. Transonic potential flow calculation using conservation form. Proc. — AIAA 2nd Comp. Fluid Dyn. Conf., 1975.
M. BallhausW. F., Jameson A., Albert J. Implicit Approximate — factorization schemes for the efficient solution of steady transonic flow problems. — Proc. AIAA 3rd Comp. Fluid. Dyn. Conf., 1977.
12. .Holst T. L. A fast, conservative algorithm for solving the transonic full potential equation.—AIAA Paper N 79—1456, 1979.
13. Jameson A. Acceleration of transonic potential flow calculations on arbitrary meshes by the multiple grid method. — AIAA Paper
N 79—1458, 1979.
14. К а н т о p о в и ч JI. В., К р ы л о в В. И. Приближенные методы
высшего анализа. — М. — Л.: Гос. изд-во физ.-мат. литературы, 1962.
Рукопись поступила 17jVII 1984
