CC BY
47
УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том XIX 19 8 8
№ 4
УДК 533.6.011/533.6.013.42
МНОГОУРОВНЕВЫЙ МЕТОД РАСЧЕТА АЭРОДИНАМИЧЕСКИХ НАГРУЗОК НА ЛЕТАТЕЛЬНОМ АППАРАТЕ
В. Д. Ильичев, В. В. Мозжилкин
Предлагается многоуровневый по степени детализации геометрии летательного аппарата метод расчета аэродинамических нагрузок, который строится на основе панельных методов и метода вложенных сеток Федоренко. Иллюстративные расчеты аэродинамических характеристик компланарных пластин бесконечного размаха, гармонически колеблющихся в потоке несжимаемой жидкости, показывают достаточно высокую точность и эффективность такого подхода.
1. Вычислительные методы и программы расчета аэродинамических нагрузок, ориентированные на параметрические исследования в процессе проектирования летательных аппаратов (ЛА), должны отвечать одновременно требованиям и достоверности, и эффективности, чтобы можно было проводить эти исследования в «темпе» разработки проекта.
Такого рода расчеты проводятся как для ЛА в целом, так и для отдельных его частей (агрегатов), однако во всех случаях размерность задачи, зависящая от числа аэродинамических «элементов» расчетной схемы (панелей, вихрей) ограничена предельной величиной, определяемой возможностями используемых ЭВМ. Совмещение же воедино более подробных схем отдельных агрегатов с менее подробной схемой ЛА часто приводит к сильным искажениям в виде появления колебательности решения: Однако значительного продвижения в эффективности решения задач с большой исходной размерностью можно добиться, применяя так называемые «многосеточные» методы [1], комбинируя их с многоуровневой декомпозицией исходной задачи.
Будем применять классификацию вычислительных методов по двум признакам: уровню агрегирования Л А (Л А в целом, его отдельные агрегаты и даже отдельные элементы агрегатов) и числу вложенных сеток, используемых в расчетной схеме каждого уровня.
1. Одноуровневые односеточные методы, в частности, существующие панельные методы [2, 3]. Они могут служить основой для разработки многоуровневых и (или) многосеточных методов.
2. Одноуровневые многосеточные методы [1]. Здесь используются вложенные друг в друга сетки, причем невязка решения вычисляется на самой крупной сетке. Они могут быть использованы для получения решения одновременно в областях с различными сетками. Процедура вычисления невязки решения на крупной сетке является сглаживающей и осциллирующие искажения решения здесь не возникают. Могут быть использованы сетки различной, в том числе большой размерности при полностью фиксированной внешней геометрии ЛА, вплоть до деталей.
3. Многоуровневые односеточные методы можно рассматривать как обобщенные одноуровневые многосеточные, но с переходом на мелкую сетку отдельно на каждом агрегате. При этом возможен упрощенный подход, удобный в проектировочных расчетах, когда интерференция агрегата и ЛА учитывается на основе упрощенной геометрической схемы ЛА, а детали геометрии учитываются на уровне агрегатов отдельно в каждом агрегате.
4. Многоуровневые многосеточные методы.' В них, кроме расчета ЛА на крупной сетке, осуществляется поагрегатный расчет ЛА на мелкой сетке. В результате получаем решение на мелкой сетке для всего ЛА. Очевидно, что реализация этих методов требует создания единой базы данных, единой математической формулировки и совместимости алгоритмов.
2. Основные положения. Сопоставим теперь возможности перечисленных выше методов расчета нестационарных и стационарных аэродинамических нагрузок с вопросами, которые возникают в связи с решением комплексных задач аэроупругости и анализа состояний ЛА при широкополосном внешнем воздействии на него, от квазистатичес-кого до акустического. Эти задачи, решаемые на основе конечно-эле-ментных упруго-массовых расчетных схем, описываются уравнениями очень большой размерности, так как учет низких частот движения упругого ЛА требует рассмотрения его в целом, а учет высоких частот—■ подробного описания внешней и внутренней пространственной геометрии агрегатов ЛА в расчетной схеме \[4].
Такие комплексные задачи решаются методом математической декомпозиции (ММД), когда система уравнений движения и соответственно расчетная схема ЛА разделяются по уровням агрегирования ЛА и по диапазонам частот внешнего воздействия на частные подзадачи, (имеющие размерности не более предельной, определяемой возможностями ЭВМ) с сохранением связей между ними. Это положение иллюстрируется диаграммой декомпозиции, приведенной в таблице, где в трех левых колонках приведены в однозначное соответствие уровни агрегирования ЛА, уровни характерных частот движения ЛА, его агрегатов, его элементов и уровни диапазонов движения: движение полета ЛА в «большом» движение ЛА в «малом» (упругие деформации), движение в «малоти» агрегатов в системе конструкции ЛА, движение в «малом» элементов в системе конструкций агрегатов.
Если деагрегирование схемы ЛА на агрегаты, агрегатов — на элементы проведено так, что на каждом его уровне — нулевом, первом, втором и третьем — основные собственные частоты движения ЛА и его подконструкции (т. е. агрегатов, элементов) покрывают диапазон частот внешнего воздействия, имеющий номер данного уровня, (см. таблицу), то в диагональных клетках диаграммы оказываются задачи динамики подконструкций, ниже них — задачи, в которых реакцию подконструк-ций на внешнее воздействие можно считать квазистатической, а частные задачи, относящиеся к клеткам выше диагональных, соответст-
Номер диапазона 0 1 1 2 1 3
Уровни агреги- рования Уровни харак- терных частот Диапазон ^ч частот ^ч. воздействия Уровни движения ^ч_^ Управление полетом и статическая аэроупругость Динамическая аэроупругость Акустическая аэроупругость Акустическая прочноость
0 10» Движение полета ЛА в большом Динамика полета ЛА в большом С Р Е з
I 101 Движение ЛА в малом Квазистатистика ЛА в малом Динамика ЛА в малом Ч А С
II 102 Движение 'агрегата ЛА в малом Квазистатика агрегата ЛА Квазистатика агрегата ЛА Динамика агрегата ЛА О Т
III Юз Движение элемента агрегата в малом Квазистатика элемента агрегата Квазистатика элемента агрегата Квазистатика элемента агрегата Динамика элемента агрегата
вующих собственным частотам более высоким, чем основные, как обычно, не рассматриваются («срез» частот в таблице).
Так как предельная размерность задач ограничена возможностями ЭВМ, то для всего ЛА на уровне / диаграммы вместо предельно подробной схемы придется использовать упрощенную, с числом элементов не более предельного, и лишь на втором (или даже на третьем) уровне диаграммы для каждого агрегата (элемента) применять подробную расчетную схему, имеющую сетку с достаточной разрешающей способностью, а иногда и содержащую такие детали геометрии, которых не было в упрощенной схеме уровня I.
Далее, зная координаты и параметры элементов аэродинамической и упруго-массовой схем, можно на основе декомпозиции и перечисленных выше методов аэродинамического расчета определять аэродинамические операторы уравнений движения ЛА и его частей.
В ряде случаев такой расчет на уровне т. е. для ЛА в целом, может быть проведен одноуровневым односеточным методом. На уровень II, т. е. на агрегаты, аэродинамическая нагрузка может быть перенесена либо путем интерполирования распределения давления, либо в форме узловых сосредоточенных сил с помощью интерполяционных соотношений, связывающих перемещения узловых точек элементов схемы уровня I и уровня II из условия равенства виртуальных работ узловых сил уровня / и уровня II (см. [3]). Очевидно, что при этом перенос нагрузки на мелкую сежу уровня II не дает увеличения точности. Кроме того, если геометрия ЛА сложна, то аэродинамическая схема ЛА уровня I, отвечающая требованиям к точности расчета, будет иметь слишком большую размерность. В ряде случаев эту трудность можно преодолевать, применяя на уровне / одноуровневый метод с последующим переносом нагрузки на мелкую сетку уровня II. Алгоритм такого метода более эффективен и может применяться как в конечном, так и в итерационном вариантах. Однако в некоторых случаях погрешности, возникающие при переносе нагрузки с уровня I на уровень II, могут быть неприемлемо велики. Многоуровневый односеточный метод позволяет избавиться от переноса нагрузки, так как аэродинамический расчет проводится как на уровне 1 для ЛА, так и на уровне II — для агрегатов. В последнем случае аэродинамическая интерференция агрегата с ЛА может быть учтена приближенно, на сетке уровня I. Уточнение учета интерференции может быть достигнуто применением многоуровневого многосеточного метода.
Рассмотрим теперь вопрос построения многоуровневых алгоритмов нестационарной аэродинамики.
Алгоритмы методов решения задач нестационарной аэродинамики сводятся к решению матричного уравнения
Л/С=^, (2.1)
где А — матрица коэффициентов, К — искомый вектор, № — известная правая часть.
Одноуровневый многосеточный алгоритм. Пусть на уровнях / (1=1) и II (1 = 2) имеют место системы
А'К^ХГ,, 1=1,2. (2.2)
Решение системы для Кг связано с большой размерностью, превосходящей ресурсы современных ЭВМ. Поэтому применим для ее решения метод Федоренко [1].
Пусть —я-я итерация, а /((«+'>= /С<п> + ^2. Тогда для невязки У2 имеет место уравнение:
А2У2=Ш2-А2К™.
Интерполируя это соотношение на уровне I, имеем:
А, V, ^ 1\А2 У2 = 1\(Ч?% - А2 К™).
Здесь и далее через ^ обозначена матрица интерполяции с матрицы уровня у на уровень а. Отсюда можно получить итерационные формулы для Кг.
АЧИ-1) = К(2п) + Л2 Л Г172 (^2 - ^2 Цп)) (2-3)
или
К^^КЫ + ЦА^^г -1\А.гКМ). (2.4)
Критические по ресурсам памяти умножения А2К2 и интерполяции сводятся к перемножению клеточных матриц и легко реализуемы на ЭВМ. Обратная матрица вычисляется только на уровне I. Это существенно экономит и ресурсы памяти и затраты машинного времени.
Пусть имеется N агрегатов, на каждый агрегат приходится t узлов сетки уровня II. На уровне I тоже имеется ^ узлов. Тогда, если положить, что при интерполяции с сетки II на уровень I требуется т1 операций, а с уровня I на уровень II т2 операций, то на одну итерацию по формуле (2.3) требуется
(2^— 1) + N— 1] + т2) 4- £(2£ — 1) (2.5)
операций, а по формуле (2.4) число операций равно
52 = Л[М(2/ — 1) + ./V — 1] 4- М + Ы(тх + т2)-{- Ц21 — 1) + *. (2.6)
Очевидно, 52<54.
Многоуровневые многосеточные алгоритмы. Без нарушения общности можно считать, что рассчитываемому на уровне II агрегату а соответствует верхняя диагональная клетка матрицы А. Тогда для агрегата получаем соотношения для невязки
А'ъ, 1Л - + 4; V, .'= К -А*.кп - А* К»- |
Здесь а дополнение агрегата до всего ЛА, а мелкая сетка панелей используется только на агрегате уровня II. Поэтому матрица Ааа играет роль интерполятора. Расчет я+1-го приближения осуществляется на агрегате интерполяцией
+ (2-8)
Отметим, что алгоритм (2.3) требует поагрегатного хранения клеточных матриц с числом элементов ЫН2. Если пересчет с крупной сетки на мелкую производить поагрегатно по формулам (2.7), (2.8), то тре-
буется Л/72(3—2/Л/а) элементов. Поэтому алгоритм (2.7) при N<2 требует ресурсов памяти, одинаковых с (2.3). Если Ы>-3, то алгоритм (2.7) более экономичен по ресурсам памяти. Алгоритм (2.7) можно дополнительно упростить, если положить
ЦгАа-аК^А'а-КП, АіаКу~АіаІ\К&.
Тогда формулы (2.7) перепишем в виде
Ка К + ЛК. У1а = Г„ (Г. - А1КУ) - Аа-а №>, Аі V' + АVI = ЧГъ - А і Я А™ — А!_ К™ ■
аа и аа а % аа £ а аа а
(2.9)
Под хранение матриц здесь требуется Л/д^2 элементов. Поэтому данный алгоритм является самым экономичным по используемой памяти ЭВМ.
Алгоритм (2.7) требует на одну итерацию по всем агрегатам ЛА число операций
— 1) 4- + тх + тг — ^8 — ' (2-Ю)
Можно проверить, что 53<52 при 4. При Л/<3 алгоритм (2.5) более быстродействующий.
Отметим, что алгоритмы (2.7) и (2.9) учитывают влияние агрегата уровня II. на его дополнение уровня I.
Алгоритмы (2.3) — (2.9) можно использовать как для одноуровневого расчета на вложенных сетках, так и для многоуровневого. Заметим, что интерполяция с уровня II на уровень I упрощается, если использовать на уровне II вложенные сетки панелей, кратные 3 по отношению к уровню I. В этом случае она сводится к вычислению соответствующих значений в точках коллокации уровня /.
Предлагаемые формулировки многоуровневых многосеточных методов позволяют включать в рассмотрение на уровне агрегата новые геометрические детали и получать многоуровневые односеточные методы. Это осуществляется заменой в формулах (2.7), (2.9) матрицы А'аа на матрицу 1^аа, вычисленную панельным методом с учетом этих деталей.
Аналогично можно по мере уточнения компоновки ЛА вводить новые геометрические детали. Пусть на некотором этапе разработки агрегат не рассматривался. В это время формируется матрица влияния дополнения агрегата По мере конструирования агрегата
формируются матрицы влияния А^а; Аай; Ааа. На первой итерации К}%> можно вычислить по формулам:
АааКЫ + А~ааК^=Ша, после чего можно применить формулы (2.7) и (2.9).
3. Рассмотрим задачу обтекания двух тандемно расположенных профилей, совершающих гармонические колебания в потоке несжимаемой жидкости. Эту задачу можно разбить на две: обтекание двух симметричных профилей под нулевым углом атаки и обтекание гармонически колеблющихся пластин. Первая задача решается в явном виде методом источников. Поагрегатный или многоуровневый подход приводит здесь к интегрированию на более мелкой сетке панелей и поэтому тривиален.
Задача расчета обтекания гармонически колеблющихся пластик сводится к интегральному уравнению
Здесь предполагается, что пластины расположены в областях хє[0; а]; ^є[й; с]; до— скос потока, со — число Струхаля, К(х')—интенсивность распределения диполей на пластинах.
Применяя к (3.1) процедуру панельного метода, получим:
Функции £(«) и {(б) вычисляются по формулам из работы [5].
Рассмотрим многосеточный алгоритм расчета обтекания одной пластинки л:е[0, 1].
(х' — х)2
е-іш(а-4г')
сіх’ +
О а
(3.1)
— 2тКиІ(х і у—і, х, у, х) + К{а)е
X
X [•/{а — х, ш) — 1{Ь — х, <в)] + ^2 ^ К2]І (*2у-і , *2 у. •*) +
/=і
+ К (Ь) е-«“(с-х) J(c — X, ш),
где
К?=Іі2і\ 7 = 1, 1(хи хъ х)= {(х1 — х)(х2 — х)}-1.
Если '•<(), то
/(-с, <в) = е~іап |і- —Іт gі (а>т) — г/(ап)]}.
Если х < 0, то
У(т, <ю) = аш + |~ — іш \g(^o\ т І) -І- і/(ш І Т І )]\ еы^
оо
оо
На рис. 1 приведены данные распределения коэффициента давления Аср при со = 0,1 и колебаниях вида у = хеш. Можно отметить хорошее совпадение результатов многосеточного метода ЛА11 = 12, точного решения при со = 0 (штриховая линия) и одноуровневого расчета при Ы—18 (сплошная линия).
Рассмотрим обтекание двух компланарных пластин
На рис. 2 приведены данные распределения Аср при со = 0,1 и колебаниях вида у = хеш. Результаты расчетов на мелкой сетке при ^ = 6; Ы2=9 (сплошная линия) и многоуровневого метода при —9 (светлые кружки) хорошо согласуются.
Точность многосеточного и многоуровневого алгоритма существенным образом зависит от способа межуровневой интерполяции. В расчетах использовались интерполяционные формулы, учитывающие особенности на кромках профиля:
а) в окрестности передней кромки
у = 0; * £ [0, 1]; *£[1,5; 2].
где | — расстояние до передней кромки;
б) в окрестности задней кромки
К% = Кт (1 + ш 0,5 Л2) + К" (0,5 Л2)3/2,
з
К” = \К}т - К Т (1 + І* 1,5 Аа)] (1,5 Н2) \ Кт = ют а - р/Сед-Г, а = 5,19615{5; р = 0,238314/(4,19615 + ^7.3,29424), а= 1,2.
24 г
20 І
* ё —р-х-ох-9-]Х-С<->Р-^-ох-хох^рх -ргх&і
о
0,2 0,4 0,6 0,8 х 1
3—«Ученые записки» № 4
33
Рис. 2
Здесь Й1 и й2 — хорды панелей уровня / и II на агрегате.
Во внутренних точках используется параболическая интерполяция.
Возможности многоуровневого расчета можно проиллюстрировать на примере обтекания пластин хе[0, 1] и хе[1,5; 1,7] при колебаниях вида у = еш, © = 0,1. На рис. 3 сплошной линией показаны данные многоуровневого многосеточного расчета при Л71 = Л^ = 6. Кружки соответствуют распределению Аср на первой пластине без учета влияния второй, при этом кружки соответствуют одноуровневому многосеточному расчету на первой итерации, точки — данные 6-й поагрегатной итерации. Совпадение результатов свидетельствует о сходимости данного метода.
4. Оценки скорости сходимости метода. Можно показать, что невязка одноуровневого многосеточного метода определяется формулой
У^) = [Е-1\(А,)-Ч\А2\ ум.
Тогда скорость сходимости итераций соответствует убыванию членов геометрической прогрессии со знаменателем
9=||^_/?(Л1)->/>Л2[|.
Аналогично можно показать, что скорость сходимости многоуровневого многосеточного алгоритма определяется величиной
0 Е Е,Е1 ’ где
Е,- Е -11{А'аа)-Ч\А1а>
При проведении расчетов, описанных в п. 3, относительная погрешность К, равная 10-6, достигалась за 10—12 итераций. Аналогичные результаты были получены при пробных расчетах обтекания гармонически колеблющихся крыльев и тандемно расположенных несущих поверхностей.
Предлагаемый подход позволяет в рамках любого численного метода решения интегральных уравнений линейной аэродинамики ЛА создавать многоуровневые одно- и многосеточные алгоритмы. Его рабо-тоспособость проиллюстрирована на примерах.
Практическая реализация многоуровневого метода аэродинамического расчета требует построения базы данных. Она включает в себя поагрегатное представление геометрии ЛА и результаты расчетов агрегатов на различных сетках и геометрических моделях. Система многоуровневых методов, состоящая из одноуровневого односеточного метода для 'инженерных расчетов агрегатов и многоуровневого многосеточного метода для уточнения нагрузок на агрегате и поверочных расчетов всего ЛА функционирует над базой данных. Пользователю должна быть представлена возможность как для итерационного, так и для по-
шагового использования алгоритмов. Эта же база данных может быть основой для разработки и применения многоуровневых методов аэродинамического расчета JIA.
ЛИТЕРАТУРА
1. Федоренко Р. П. О скорости сходимости одного итерационного процесса. — Ж. вычисл. мат. и мат. физ., 1964, т. 4, № 3.
2. 3 а х а р о в А. Г. Численный метод расчета аэродинамических характеристик плоских и неплоских крыльев в сверхзвуковом потоке. — Ученые записки ЦАГИ, 1973, т. 4, № 3.
3. Гостев П. М., Сухарев В. М., Мозжилкин В. В., Шевы-р е в С. П., С а п у н к о в Я- Г. Метод и алгоритм расчета нестационарных аэродинамических характеристик ЛА на дозвуковых скоростях полета. —
Труды ЦАГИ, 1986, вып. 2305.
4. Ильичев В. Д. Матричные методы синтеза динамических и упругих характеристик линейных неконсервативных конструкций.—-Ученые записки ЦАГИ, 1975, т. 6, № 2.
5. D a t R., М а 1 f о i s J.-P. Sur le calcul du noyau de l’equation de la surface portante en e’coulement subsonique instationaire. — La Rech. Acrosp., 1970, N 5.
Рукопись поступила 30/111 1987 г..
