CC BY
130
________УЧЕНЫЕ ЗАПИСКИ ЦАГИ______________
Том XXX ' 1999 ' №3—4
УДК 629.735.33.015.4:533.6.013.422
РАСЧЕТНЫЕ ИССЛЕДОВАНИЯ ТРАНСЗВУКОВОГО
ФЛАТТЕРА
Ф. 3. Ишмуратов, С. И. Кузьмина, В. А. Мосунов
Разработан один из вариантов математической модели движения упругого самолета с органами управления в трансзвуковом потоке. Математическая модель основана на полиномиальном методе Ритца (программа КС-2).
Проблема исследования нестационарного обтекания упругих деформируемых крыльев в околозвуковом диапазоне чисел Маха решается при помощи конечно-разностного метода Годунова интегрирования уравнений Эйлера (программа TRAN).
В частотной области применяется итерационный метод расчета трансзвукового флаттера. При вычислении обобщенных аэродинамических сил используется принцип гармонической линеаризации.
Получены уравнения движения упругого самолета с органами управления с учетом нестационарных нелинейных аэродинамических сил при трансзвуковом обтекании. Разработаны метод и алгоритм решения уравнений во временнбй области (программа TRTDR).
Предложенные методы иллюстрируются расчетом характеристик флаттера двух крыльев в зависимости от различных параметров. Приведены сравнительные результаты расчетов линейной, линеаризованной и нелинейной систем.
1. Традиционный подход к расчетным исследованиям флаттера основан на линейной аэродинамической теории. Этот подход позволяет изучать подавляющее большинство встречающихся форм и типов флаттера. Большие трудности возникают при расчете флаттера современных самолетов в околозвуковом диапазоне чисел Маха. Эксперименты показывают, что в этом диапазоне происходит существенное уменьшение критического скоростного напора, обусловленное движением скачков уплотнения по поверхности крыла. Кроме того, некоторые типы самолетов сталкиваются в околозвуковом полете с явлением, известным как предельный цикл колебаний. Это тип флаттера, где аэродинамические силы существенно нелинейны, и для исследования которого невозможно использование методов классического линейного анализа аэроупругих колебаний [1].
Решение задачи аэроупругой устойчивости в трансзвуковом потоке непосредственно связано с изучением особенностей обтекания колеблющихся стреловидных крыльев. Для расчета течений с местными сверхзвуковыми зонами широко используются численные методы. Определение нестационарных аэродинамических характеристик составляет самую сложную и трудоемкую часть задачи о флаттере.
Применение нелинейной аэродинамической теории, позволяющей, в отличие от линейной, учитывать толщину профиля, амплитуду колебаний, движение скачков уплотнения, привело к созданию новых методов решения задачи о флаттере.
В данной статье излагаются два подхода к исследованию трансзвукового флаттера: в частотной и временной областях. Аэродинамические нагрузки вычисляются конечно-разностным методом Годунова интегрирования нелинейных нестационарных уравнений Эйлера для идеального газа (программа TRAN).
В частотной области применяется итерационный метод исследования, который основывается на классическом линейном подходе [2], [3] к определению деформаций упругих поверхностей (программа КС-2). На каждой итерации аэродинамические характеристики определяются путем решения нелинейных уравнений Эйлера для исследуемого распределения амплитуды колебаний [4].
Анализ аэроупругих колебаний проводится также во временной области. Применяется метод совместного численного интегрирования уравнений, описывающих перемещения несущей поверхности, и аэродинамических уравнений трансзвукового обтекания (программа TRTDR). Разработанная математическая модель и созданная на ее основе программа расчета предназначены в основном для исследования неоднородных задач динамической аэроупругости, включая задачи аэросервоупругости. В дальнейшем предполагается использовать данный подход для численного моделирования трансзвукового обтекания упругого самолета с учетом функционирования систем управления.
2. Уравнения движения. При выводе уравнений движения упругого самолета будем исходить из методов и алгоритмов, принятых в комплексе программ КС-2 [4]. Используется метод Ритца, когда деформации упругих поверхностей (УП) представляются в виде степенных функций от пространственных координат. При этом вся конструкция моделируется набором тонких, первоначально плоских УП, которые могут быть расположены в пространстве произвольным образом. Для каждой из УП задается распределение масс и жесткостей. Используются различные типы элементов для схематизации реальной конструкции: сосредоточенные массы, балки, работающие на изгиб и кручение, пластины и т. п. Локальные системы координат для каждой упругой поверхности выбираются таким образом, чтобы плоскость xOz совпадала с плоскостью УП. В единую расчетную модель упругие поверхности объединяются с помощью упругих связей, которые позволяют моделировать различные условия закрепления в местах соединения УП между собой.
N
IV(х, 2, /) = ]Г(х, г) и*(О,
*=1
где
/к{х,г) = хщ гПк, т,п = 0,1,—
Для каждой УП может быть выбран свой полином. Коэффициенты щ (/) являются обобщенными координатами метода полиномов.
На УП в потоке воздействуют аэродинамические нагрузки, которые могут быть вычислены по какой-либо аэродинамической теории или взяты из эксперимента.
Для расчета аэродинамических сил на УП вычисляются перемещения и углы наклона в точках аэродинамической расчетной сетки. Иногда бывает необходимо использовать формы деформации упругой поверхности, заданные в виде перемещений в М точках, не совпадающих с точками аэродинамической расчетной сетки, когда перемещения получены, например, расчетом по методу конечного элемента или из эксперимента. В этом случае используется интерполяционная формула [5]:
м
= +Ъ\ +1>1Х + Ы12,
(=1
где
Неизвестные коэффициенты а, и Ъ, определяются из заданных перемещений в Л/точках х„ г, и дополнительных условий:
ММ м
Хл =0; =0-
1=1 1=1 г=1
Таким образом получаются матрицы уравнений движения для объединенного вектора метода полиномов для всей системы упругих поверхностей. Для решения задач динамической аэроупругости (флаттер, динамическая реакция, аэросервоупругость) уравнения приводятся к обобщенным координатам, соответствующим формам собственных колебаний вне потока. В матричной форме они могут быть записаны в виде
Сц + 1)д^ + Сг([ = 0° + Л6 г + р1', у = Н0д + Нхд + Н2<1,
где • ,
0 = со\(дг,де,5) — вектор обобщенных координат, описывающий движение самолета как твердого тела qr, его упругие деформации де, отклонения органов управления 8;
С, 2>о, С — матрицы инерции, конструкционного демпфирования и конструкционной жесткости;
— вектор обобщенных аэродинамических сил;
8, — вектор отклонений штоков рулевых приводов;
/? — матрица эффективности воздействия рулевых приводов;
^— вектор внешних сосредоточенных сил;
у — вектор выходных параметров, включающий смещения, углы поворота, ускорения и угловые скорости в местах установки датчиков, а также силовые факторы в выбранных сечениях конструкции;
Но, Н\,Н2 — матрицы перехода от обобщенных координат, их скоростей и ускорений к физическим координатам.
Аэродинамические силы в линейном случае определяются в частотной области (для гармонического движения). В дозвуковом потоке используется метод дискретных вихрей [2], в сверхзвуковом — панельный метод [3]. Обобщенные силы во временной области могут быть определены при этом с применением гипотезы гармоничности. В этом случае:
0° =-рУОд-рУ2Вд + р¥ И™
где
О, В — матрицы аэродинамического демпфирования и аэродинамической жесткости, вычисленные для заданного числа Струхаля;
р — плотность воздуха;
V — скорость полета;
2)" — вектор эффективности воздействия порыва (в предположении мгновенного охвата конструкции порывом);
уу — интенсивность воздушного порыва.
Решение однородных и неоднородных уравнений с линейной аэродинамикой в частотной и временной областях рассмотрено в работах [2], [3], [6], [7].
Расчет трансзвукового обтекания крыла проводится путем решения нелинейных уравнений, поэтому его результат зависит как от формы деформаций крыла, так и от их амплитуды. Отметим, что в расчетах по линейной теории формы деформаций могут определяться с точностью до произвольного множителя, в то время как при расчетах трансзвукового обтекания они должны иметь определенную амплитуду, от которой зависит распределение разности давлений на поверхности крыла Ар(х, г, ¿).
При решении задачи о флаттере в частотной области используется линеаризованное распределение разности давления. Предусмотрены два подхода к линеаризации, предложенные В. Г-. Буньковым.
В первом подходе, известном как гипотеза одномерной стационарности (ГОС), предполагается, что распределение разности давления вблизи заданного угла атаки можно представить в виде
где
а{х, г, I) ----------------местный динамический угол атаки;
дх V
Ср (х, г) — производная коэффициента давления по углу атаки; определяется численным дифференцированием по результатам двух расчетов статических распределений давления Арх и Ар2 для двух близких значений углов атаки а і и а2:
Тогда на основании ГОС могут быть получены коэффициенты аэродинамических матриц В и Б.
где & — площадь поверхности крыла.
Такой подход позволяет выполнять расчет на флаттер с учетом особенностей трансзвукового распределения давления при заданном угле атаки. Однако при этом не учитывается нестационарность аэродинамических сил. Кроме того, применение метода, основанного на ГОС, для расчета рулевых форм флаттера может привести к большим погрешностям, так как в нем не учитываются особенности обтекания в окрестности элеронов (рулей) при их отклонении.
Поэтому в дальнейшем была реализована идея использования вместо стационарного нестационарное трансзвуковое решение именно для тех форм возмущения потока, которые участвуют в возникновении флаттера. Был разработан новый метод, позволяющий преодолеть указанные трудности и соединить вместе имеющиеся два алгоритма: расчет флаттера по линейной теории и расчет трансзвукового обтекания упругого крыла, колеблющегося с заданной частотой и формой (амплитудой).
Ср(х,г) =
2 Ар1(х,г)-Ар2(х,2)
р Vі щ-а2
1Г{х,г,Ъ = Р(х,г)е1м,
где
Р(х, г) = — комплексная форма колебаний (у = д/-1);
со — заданная круговая частота колебаний.
Для такого гармонического движения вычисляется распределение давления, из которого выделяются коэффициенты основной гармоники разложения Фурье:
Ар(х, г, 0 »[Ад 0,2) + у'Арг (х, £)\е}
'Ш
(2)
Далее вводится гипотеза одномерной динамичности (ГОД), согласно которой гармонически линеаризованное давление определяется комплексным местным углом атаки а = щ+ (при этом давление в общем случае не совпадает по фазе с а):
р V _ _
Др(*,г,/) = — (с“ +ус“ )(а! +уа2).
(3)
Динамический угол атаки при гармонических колебаниях представляется в виде
а(х,г, 0 = -
д1_
дх Ъ
(4)
где
к~а>Ь /V — число Струхаля;
Ь — характерный размер.
Из уравнений (2) — (4) можно найти коэффициенты нестационарного давления:
&Р\
к л
Сп =■
Р\
а,1“»* +А*
дх Ь
л
У.
рУ2
8Р1 к дх Ь
л2
+
дх Ь
с“ =
Р2
. (дрг к -Ар, —+ тЧ | + Дй --ТЪ
дх Ъ
Для коэффициентов аэродинамических матриц В и I) могут быть получены следующие выражения:
Таким образом может быть организован итерационный процесс, в котором расчет аэродинамических сил основан на методе вынужденных колебаний и гармонической линеаризации. Так как для получения Др(х,л,/) необходимо предварительно задать деформации крыла и частоту колебаний при флаттере, то в качестве исходного приближения используются результаты расчета флаттера для случая, когда аэродинамические матрицы вычисляются по какой-либо линейной теории. Затем по нелинейной трансзвуковой теории рассчитывается нестационарное обтекание крыла. Из вычисленного на этом этапе динамического давления выделяются синфазные и квадратурные составляющие основной гармоники, вычисляются коэффициенты аэродинамических матриц, и уравнение колебаний в потоке решается заново. Уравнения в этом случае выглядят так же, как при использовании линейной аэродинамики. При этом получаются новые частота и форма флаттера. Далее процесс повторяется до сходимости (например, по частоте флаттера). В том случае, когда по линейной теории флаттер отсутствует, в качестве начального приближения может быть использована одна из форм собственных колебаний в пустоте.
Применение конечно-разностного метода Годунова интегрирования уравнений Эйлера для расчета нестационарного решения требует дополнительных пояснений. Существует принципиальное отличие между такими понятиями, как стационарное и нестационарное решения, полученные методом установления, хотя в обоих случаях формально производится одна и та же процедура численного интегрирования. Параметры потока, соответствующие каждой ячейке трехмерной расчетной сетки, являются ступенчатыми функциями времени. Шаг по времени определяется из условий устойчивости численной схемы. Найти стационарное решение — это значит вычислить установившееся по времени поле течения вокруг неподвижной несущей поверхности, внезапно помещенной в невозмущенный поток. Требуется достаточно большое число шагов по времени, чтобы параметры потока в каждой ячейке установились с заданной точностью. Нестационарное решение соответствует расчету обтекания колеблющейся с некоторой частотой поверхности. В качестве начальных условий можно взять либо невозмущенный поток, либо стационарное решение. Для определения нестационарного решения применяется подвижная расчетная сетка. Ско-
• дЖ 81Г '
роста границ ячеек Рки их углы наклонов-------и-----определяются коор-
■ дх дг
динатными функциями перемещений/{х, г) и частотой со. Найти нестацио-
нарное решение (в случае гармонических колебаний) — это значит вычислить установившееся периодическое решение, т. е. параметры потока в каждой ячейке должны быть периодическими функциями. Интегрирование при этом проводится в течение нескольких периодов колебаний.
Вопрос о том, что такое нестационарное решение, полученное методом Годунова при произвольных движениях поверхности, достаточно сложен и требует детального рассмотрения в каждом конкретном случае.
В задачах аэроупругости движение деформируемой поверхности, как правило, описывается временной функцией . Как показали расчеты
и проведенные сравнения с известными экспериментальными данными, нестационарные решения колебательного характера, полученные численным методом установления, соответствуют реальным физическим процессам. Термин «установление» при решении таких нестационарных задач указывает на то, что процесс может идти только в направлении увеличения временного параметра где г — число пройденных шагов по
времени.
При исследовании трансзвуковых явлений аэроупругости во вре-меннбй области выполняется совместное численное интегрирование уравнений (1) и уравнений Эйлера трансзвукового обтекания. На каждом шаге интегрирования вычисляется все поле скоростей, давления и плотности методом Годунова [4]. Граничные условия в контрольных точках аэроди-
намическои сетки движущейся упругой поверхности Ш, -—, ---------, вычис-
. дх дг
ляются через векторы обобщенных координат и скоростей:
ттг —,, г' / \ чгХтт 81^
W = XUq -w(x,z,t); —- = Х Иц\ ----------X Щ,
дх дг :
где .
X, Xх, X2 — соответствующие полиномиальные матрицы перехода;
и—модальная матрица;
м>(х, г, 0 — распределение скоростей воздушного порыва.
По полученному распределению разности давления Ар{х,г^) определяются обобщенные аэродинамические силы:
=ХтиТ8Ар,
где 51 — диагональная матрица площадей аэродинамических клеток.
Таким образом, для вычисления правой части уравнения (1) в заданный момент времени используется все пространственное поле параметров течения, зависящее от предыстории процесса, толщины профиля, амплитуды колебаний и движения скачков уплотнения. Интегрирование уравнений выполняется методом Эйлера 1-го порядка. Для выбора оптимального шага интегрирования предварительно исследуется система с линейной аэро-
o,s
M
M
• Эксперимент
-----TRAN, Г ^
-----TRAN, 6’
-----Линейная теория, K.C-2
CA М
В,7 в, 8 *>-» М 11
динамикой. Обычно величина шага ограничивается устойчивостью разностной схемы Годунова, а не упругими колебаниями.
Ниже рассмотрены примеры, иллюстрирующие возможности предложенных методов расчета трансзвукового флаттера в частотной и временной областях.
3. Крыло AGARD 445.6. Данное крыло часто используется для «сертификации» и сравнения расчетных методов. Сравнительные результаты приведены на рис. 1 для варианта «2,5 foot weakened model 3» [8]. В расчете были использованы первые 4 тона колебаний. Их частоты перечислены в таблице.
Значения безразмерного скоростного напора флаттера qj-
для каждого числа Маха вычислены по линейной нестационарной дозвуковой и сверхзвуковой аэродинамической теории (штрих-пунктирная линия на рис. 1).
В применяемой трансзвуковой теории давление на крыле существенно зависит от амплитуды колебаний конструкции. В данном примере амплитуда колебаний характеризуется углом закручивания конца крыла ао.
На рис. 1 представлены результаты расчетов по программе TRAN для двух значений амплитуды колебаний а0=1° (штриховая линия) и а0 = 6° (сплошная линия).
Как видно из рис. 1, изменение амплитуды колебания влечет за собой изменение скорости флаттера, так как при больших дозвуковых скоростях полета аэродинамические характеристики становятся существенно нелинейными из-за учета движения скачков уплотнения и изменения их интенсивности как за счет перемещения, так и за счет скорости перемещения.
В трансзвуковом диапазоне чисел М линейная аэродинамическая теория дает значения qj-, превышающие экспериментальные, в то время
как в чисто дозвуковой и сверхзвуковой областях согласование с экспериментом вполне удовлетворительное. Применение итерационного TRAN-метода, основанного на нелинейной аэродинамической теории, позволяет значительно уточнить результаты. Отличие между экспериментальными данными и результатами расчетов по нелинейной трансзвуковой теории в
Рис. 1. Зависимость относительного скоростного напора флаттера от числа Маха для крыла 445.6 .
№ ¿Гц NASA / Гц ЦАГИ Название тона
1 9,6 9,56 Изгиб 1-го тона
2 38,1 38,09 Кручение 1-го тона
3 50,7 48,15 Изгиб 2-го тона
4 98,5 92,04 Кручение 2-го тона
а) возмущение в вкае силы на конце крыла
б ) возмущение в виде начального условия 4|=(0)=0,1
х TRAN
Рис. 2. Временные зависимости для крыла АвАИО 445.6
основном обусловлено, по-видимому, тем, что в ТЯА1Ч-методе интегрируются уравнения Эйлера, не учитывающие вязкость. В работе [8], из которой взяты результаты эксперимента, не указаны конкретные амплитуды флаттерных колебаний. Результаты расчета при больших амплитудах лучше согласуются с экспериментом, так как начало флаттера в эксперименте обычно определяют по изменению декремента колебаний при заметной величине амплитуды.
Полученные в результате расчетов временные процессы при двух разных типах начального возбуждения для числа М = 0,96, плотности и скорости звука, соответствующего началу флаттера в эксперименте, приведены на рис. 2. В случае а) сила -Г=2,5 кг кратковременно прикладывается на конце крыла, в случае 6) начальное условие представляет собой форму первого собственного тона колебаний.
В обоих случаях временной процесс, вычисленный по линейной теории, является устойчивым (положительное демпфирование), в то время как нелинейная теория дает неустойчивый характер колебаний (антидемпфирование); частота колебаний 14,5 Гц. Результаты расчета флаттера во временной и частотной областях хорошо согласуются между собой.
4. Автоколебания руля направления. На рис. 3 показана форма киля 10%-ной относительной толщины с рулем направления. Жесткость на изгиб и кручение киля были заданы достаточно высокими, чтобы конструкцию можно было рассматривать как систему с одной степенью свободы. Известно, что трансзвуковые автоколебания органов управления (buzz) возникают при отрицательном аэродинамическом демпфировании [9]. Исследовалось влияние ам-
0,3
0,2
0,1
0
- \
\ я L-J
М= 0,95
со = 30 Гц
4" 8" 12“ \lV 20° 24° S
-0,3-
Рис. 3. Зависимость производной шарнирного момента т
от амплитуды колебаний
плитуды колебаний руля 5 на величину производной шарнирного момента
с •
по угловой скорости 5 . Предельный цикл реализуется при амплитуде
8о, когда аэродинамическое демпфирование меняет знак (рис. 3). Особенность исследуемого режима автоколебаний состоит в том, что колебания руля связаны с перемещением скачков уплотнения, причем существует определенное фазовое запаздывание между ними.
На рис. 4 показано сравнение результатов расчетов во временной области по линейной и нелинейной аэродинамическим теориям для двух значений начального отклонения руля. Использовались два подхода для учета трансзвуковых явлений во временных процессах. В одном из них обобщенные аэродинамические силы определялись на каждом временном шаге совместно с уравнениями колебаний, в другом временной процесс получен на основании принципа гармонической линеаризации. Оба подхода дают качественно одинаковые результаты. При начальной амплитуде, меньшей, чем амплитуда предельного цикла, она увеличивается со временем, стремясь достигнуть предельной величины. Когда временной процесс стартует с амплитуды, большей, чем предельная, то он является затухающим при том же числе М = 0,95 и частоте колебаний руля около 30 Гц. Линейная аэродинамическая теория (метод дискретных диполей) дает затухающий по времени процесс при любых начальных условиях. Полученные в расчетах во временной и частотной областях оценки величины амплитуды предельного цикла удовлетворительно согласуются между собой, если интегриро-
8, град
н— нелинейная трансзвуковая аэродинамика; х — гармоническая линеаризация;
• — метод дискретных диполей
вание нестационарных уравнений Эйлера для определения аэродинамических нагрузок в обоих случаях производится на одной и той же расчетной сетке и, кроме того, шаг по времени достаточно мелкий, чтобы обеспечить устойчивость конечно-разностной схемы и корректность расчета временного процесса. Однако параметры аэродинамической сетки и шаг интегрирования могут влиять на величину амплитуды предельного цикла. Вопрос о сходимости по параметрам сетки требует дополнительного исследования на более мощных компьютерах. Проведенные на компьютерах Pentium Pro 200 расчеты типичных переходных процессов занимают около 20 ч при сетке 60 х 15x30; дальнейшее увеличение размерности задачи затруднительно без увеличения мощности компьютеров.
ЛИТЕРАТУРА
1. Charles М., Denegrí. Correlation of classical flutter analyses and nonlinear flutter response characteristics.— International Forum on Aeroelasticity and Structural Dynamics.— 1997, Rome, Italy.
2. Мосунов В. A., Набиуллин Э.Н. Определение аэродинамических сил, действующих в дозвуковом потоке на упругие колеблющиеся поверхности, расположенные в разных плоскостях // Труды ЦАГИ.— 1981.
Вып. 2118.
3. Буньков В.Г., Набиуллин Э.Н. Сравнительный анализ методов расчета аэродинамических сил на колеблющемся крыле в сверхзвуковом потоке // Труды ЦАГИ.— 1985. Вып.2281.
4. Kouzmina S., Mosunov V., Karkle P. Iterative method for transonic flutter calculation.— International Forum on Aeroelasticity and Structural Dynamics.— 1997, Rome, Italy.
5. Harder R.L., Desmarias R.L. Interpolation using surface splines// AIAA Journal.— 1972.
6. Ishmuratov F., Popovsky V., Karkle P. TsAGI experience and investigations in aeroelasticity of flying vehicles.— International Forum on Aeroelasticity and Structural Dynamics.— 1997, Rome, Italy.
7. Ишмуратов Ф. 3., Поповский В. H. Объединенная математическая модель летательного аппарата с системой управления для исследования аэроупругого взаимодействия// ТВФ.— 1997. T. LXXI, №2 (624).
8. Yates С., Jr. AGARD Standard Aeroelastic Configurations for Dynamic Response: I — Wing 445.6 // AGARD Report.— 1988, N 765.
9. Кузьмина С.И. Численное исследование автоколебаний профиля с элероном в трансзвуковом потоке // Труды ЦАГИ.— 1985. Вып. 2281.
Рукопись поступила 3/VI1998 г.
