CC BY
33
УЧЁНЫЕ ЗАПИСКИ Ц А Г И Том XIX 1 9 88
№ 2
УДК 533.6.011.35 : 629.7.025.73
ПРИМЕНЕНИЕ МНОГОСЕТОЧНОГО МЕТОДА ФЕДОРЕНКО ДЛЯ РАСЧЕТА ТРАНСЗВУКОВОГО ОБТЕКАНИЯ ПРОФИЛЯ
А. П. Аралов, А. А. Шагаев
Излагается алгоритм многосеточного метода Федоренко для численного решения с применением проекционно-сеточной схемы Бубнова—Галер-кина задачи потенциального трансзвукового обтекания профиля. Приводятся результаты расчетов, которые показывают увеличение скорости сходимости многосеточного алгоритма при определенном способе построения проекционно-сеточной схемы на вспомогательных сетках.
Рассматривается задача невязкого безвихревого обтекания профиля трансзвуковым потоком газа. Искусственная вязкость вводйтся посредством модификации плотности р, определяемой обычным изоэнтропическим выражением, по предложенной в работе [1] формуле
= Р — m Да*
Здесь р„ — производная плотности вдоль линии тока, As — величина порядка шага сетки, |i — функция включения. Функция включения задается выражением
[х= [max (0,1 - М^/М2)]1/2,
где М — местное число Маха, а Мс=0,96-^0,99.
Потенциал скорости Ф определяется в результате решения задачи обтекания профиля для уравнения неразрывности с модифицированной плотностью
div (р grad Ф) = 0. (1)
Задача обтекания профиля формулируется обычным образом [2]. На поверхности профиля задается условие непротекания, а на внешней границе области течения — потенциал невозмущенного потока. Для аппроксимации задачи обтекания применяется изложенная в работе [2] процедура проекционно-сеточного метода Бубнова — Галеркина. Основное отличие от работы [2] заключается в применении глобальной криволинейной системы координат (I1, I2) для всей области течения вместо локальных изопарамет-рических преобразований в криволинейные системы координат для каждого конечного элемента (ячейки сетки). Приближенное решение ищется на сетке Qл в расчетной области £2(|\ I2) в виде линейной комбинации билинейных базисных функций
M'i./G1,!2)
Ф« 2Ф* ,ЧГ* ,(6*,6*), k.l
где k, I — индексы узлов сетки, а Ф^; — значения сеточной функции потенциала Ф ‘ в узлах. Используемая сетка и преобразование координат относятся к типу О и описаны в работе [3].
Проекционно-сеточная схема в операторном виде имеет следующий вид
1Н ф*=/* з (2)
где Vх — сеточный оператор схемы, а /Л — сеточная функция, определяемая краевыми условиями на внешней границе области течения. Во внутренних узлах схема (2) имеет вид
1 1
Е X ф*+г.1+л<1, ®*+,.1+,] = о- (3)
Г = —1 5=— 1
Выражение в квадратных скобках означает скалярное произведение базисных функций
а *=1 (3=1 ^ ^
где ^ — метрический тензор криволинейной системы координат (£>, 5а), а §■ = = ёе15«р. Для вычисления скалярных произведений (4) функции а“Р ='р’уг£г£а^ представляются в виде линейных комбинаций билинейных базисных функций
а£А&’^’ (5)
к.1
где — значения функций ав узлах сетки £2л. Следует отметить, что фактически описанная схема отличается от схемы работы [2] лишь разной аппроксимацией функций аа$ и разными квадратурными формулами вычисления интегралов (4).
Уравнение (2) решается при помощи многосеточного метода Федоренко [4]. Возможность применения этого метода к расчету трансзвуковых течений была показана
Брандтом в работах, ссылки на которые можно найти в его обзоре [5]. В работе
Джеймсона [6] метод Федоренко был впервые использован для решения схемы
конечного объема, аппроксимирующей задачу трансзвукового обтекания профиля для полного уравнения потенциала, а в работах [7, 8] он использовался для решения дискретного аналога этой задачи, построенного методом Бубнова — Галеркина. В этих работах не исследовалась зависимость скорости сходимости многосеточного алгоритма от способа аппроксимации задачи на вспомогательных сетках. Для ускорения сходимости многосеточного алгоритма в случае линейных эллиптических задач Хакбушем в работе [9] были рекомендованы некоторые правила построения операторов перехода с сетки на сетку и оператора, аппроксимирующего задачу на вспомогательных сетках. В настоящей статье рассматривается применимость рекомендаций Хакбу-ша для ускорения сходимости многосеточного алгоритма в случае задачи трансзвукового обтекания профиля.
Итерационный процесс решения задачи строится в соответствии с методом искусственной сжимаемости [1]. По значениям потенциала с предыдущего итерационного слоя п вычисляется модифицированная плотность р и определяется оператор Ьк сеточной схемы (2). Схема (2) с определенным таким образом оператором Ьн решается на текущем итерационном слое п+1 многосеточным методом Федоренко [4]. Для этого кроме основной сетки Й/1, состоящей из 128X32 ячеек, используются вспомогательные укрупненные сетки (\'=2, 4, 8). Сетки формируются исключением каждого
второго узла на каждой линии более мелких сеток £2^. В качестве сглаживающего
У
алгоритма применяется алгоритм метода верхней релаксации по точкам.
Расчет одной сглаживающей итерации на сетке £1ь дает приближенное решение Ф'1 системы уравнений (2). Затем находятся сеточные функции поправок ЗФ'* на вспомогательных сетках в результате решения систем уравнений
гл6фЧ/, = _^й (6)
при помощи сглаживающего алгоритма с начальным приближением 8Ф^ = 0. В уравнениях (6) — операторы проекционно-сеточной схемы Бубнова—Га-
леркина на сетках 2чЛ, а <ГЙ— сеточные функции дефекта, которые определяются следующим образом
Й2й = /2* (£* ф" _/*); уЛ V/*
= /** (I г'бФ Г+ 2 ), V > 2,
~2
где — оператор переноса невязок с сетки 2„й на сетку 2„й. Если (т, п) и Т . . Т
(А, I) — индексы общего узла сеток 2чА и 2їЛ і х0 действие оператора /** на се-
Т Т
точную функцию О2 описывается выражением
чк
"Л 1 і
__ //'’Л г» 2 \ _
"т.п ІЛл" ІМ 1 + |г( + 1,| + |г|.|5|.
которое следует из соотношения
УЇІ
ТГЛ & 6*) = у у ЧГ*+г.Н<^»р) (7)
т,і гА г+Т7| + ,*1 + |г1-151 ’
связывающего билинейные базисные функции на сетках 2ЧЙ и 2чЛ. Потенциал
Т
на итерационном слое п + 1 рассчитывается по формуле
Ф*+1 = Фй (5Ф2* + /“ (ЗФ4Й + /84£ ЗФ8Й)),
vЛ
где /Д — оператор билинейной интерполяции для переноса поправок с сетки
на 2УЙ. На сетке 2аЛ вычисляются две, а на сетках 24А и 2^ — десять сгла-~2
живающих итераций.
В качестве оператора А'1* были опробованы операторы 1,\н и Оба этих оператора определяются выражениями, которые получаются из выражений (3) и (4) заменой к на м/г, а их различие обусловлено разным способом вычисления скалярных произведений функций (51, £г). При вычислении скалярных произведений для определения оператора /,]* на сетках используется аппроксимация, аналогичная аппроксимации (5). При определении оператора скалярные произведения выРажаются через скалярные произведения ба-
зисных функций на более мелкой сетке 2„л при помощи представления (7). В
Т
результате для оператора выполняется рекомендованное в работе [9] соот-
ношение
V/* vЛ
тчк _ рЬ. т Т ,~2
2 ^ 2 ЧН ’
2
мй
в котором матрицы операторов и /Д удовлетворяют так называемому усло-
"2
вию Хакбуша
Скорость сходимости многосеточного алгоритма характеризуется зависимостью* величины Я/Яо, равной отношению максимальных невязок системы уравнений (2) на текущем и нулевом итерационном слоях, от объема вычислений N, выраженного в так называемых единицах вычислительной работы. За единицу вычислительной работы принимается объем вычислений, требующийся для расчета одной сглаживающей итерации на сетке Пн. Таким образом, п единиц работы многосеточного алгоритма по затратам машинного времени эквивалентны расчету п итерационных слоев односеточного алгоритма.
На рис. 1 цифрами 1, 2, 3 и 4 «указаны зависимости от объема вычис-
лений N соответственно для односеточного, двухсеточнош, трехсеточного и четырехсеточного алгоритмов с оператором в случае расчета симметричного обтекания профиля NACA 0012 при Мао = 0,8. Как видно, применение многосеточного алгоритма позволяет существенно сократить затраты машинного времени на расчет трансзвукового обтекания профиля. На рис. 2 штриховой и сплошной линиями изображены зависимости \g\RIRo\ от объема вычислений четырехсеточного алгоритма соответственно с операторами и . Применение Оператора увеличивает скорость,
сходимости многосеточного алгоритма примерно на 10%. Аналогичное увеличение скорости сходимости имело место и при расчете других трансзвуковых режимов обтекания профиля ИАСА 0012. Эти данные показывают, что при аппроксимации задачи на вспомогательных сетках в соответствии с рекомендациями Хакбуша [9] скорость сходимости многосеточного алгоритма увеличивается и в случае расчета трансзвуковых течений.
Рис. 3
На рис. 3 приведено сравнение результатов "расчета симметричного обтекания профиля NACA0012 при Моо = 0,8 с данными из работ [2, 3]. Эпюра давления 1 рассчитана при помощи настоящего метода на сетке 128x32, эпюра 2— при помощи консервативной схемы конечного объема [3] на сетке 256X64, а эпюра 3— при помощи проекционно-сеточной схемы Бубнова—Галеркина [2] на сетке 70X15. Приведенное сравнение показывает, что точность настоящего метода лежит в обычном диапазоне (см., например [10]), характерном для численных методов расчета потенциальных трансзвуковых течений.
ЛИТЕРАТУРА
1. Hafez М. М., South I., Murman Е. М. Artificial compressibility methods for numerical solution of transonic full potential equation. —
AIAA J„ 1979, vol. 17, N 8.
2. Habashi W. G., Hafez М. M. Finite element solutions of transonic flow problems. — AIAA J., 1982, vol. 20, N 10.
3. Jameson A. and others. Accelerated finite-volume calculation of transonic potential flows. — In. „Numerical Methods For the computation of inviscid transonic flows with shock waves." Eds. A. Rizzi/H. Viviand,
1981.
4. Федоренко P. П. Релаксационный метод решения разностных эллиптических уравнений. — Ж. вычисл. матем. и матем. физ., 1961, т. 1,
№ 5. ,
5. В г a n d t A. Guide to multigrid development. — Lect. Notes in Matem.— 1981, vol. 960.
6. Jameson A. Acceleration of transonic potential flow calculations on arbitrary meshes by multiple grid method — AIAA Paper 79—1458, 1979,
7. Deconinck H., Hirsch Ch. A multigrid method for the transonic full potential equation discretized with finite elements on an arbitrary body fitted mesh. — J. of Comput. Phys., vol. 48, 1982.
8. В red if M. A fast finite element method for transonic potential flow calculations. — AIAA Paper 83-0507, 1983.
9. Hackbusch W. Introduction to multi-grid Methods for the numerical solution of boundary value probrems. — In „Computational Methods Turbulent, Transonic and Viscous Flows". — Ed. J. A. Essers, 1983.
10. Rizzi A., Viviamd H. Collective comparison of the solutions to the workshop problems. — In ’’Numerical Methods for the Computation of Inviscid Transonic Flows with Shock Waves". Eds. A Rizzi/H. Viviand,
1981.
Рукопись поступила 20/V 1986 г.
