CC BY
62
Том XII
УЧЕНЫЕ ЗАПИСКИ Ц А Г И
1981
№ 4
УДК 629.735.33.015.4:533.6.013.422
РАСЧЕТ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ АЭРОДИНАМИЧЕСКИХ СИЛ НА КОЛЕБЛЮЩЕМСЯ КРЫЛЕ В СВЕРХЗВУКОВОМ ПОТОКЕ
: . V ; ’
Деформации тонкого крыла выражаются через смещения ряда узлов, образующих множество трапециевидных и треугольных конечных элементов. Нестационарные аэродинамические нагрузки на упругом колеблющемся крыле определяются как сосредоточенные силы в узлах, связанные с перемещениями в этих узлах через комплексную матрицу аэродинамических коэффициентов влияния. Интеграл, связывающий потенциал скорости со скосом потока на крыле и пелене возмущений, вычисляется с помощью специальной интегрирующей сетки с переменным шагом. Исследуется влияние, на флаттерные характеристики крыла чисел М и БЬ. Метод реализован для крыла со сверхзвуковыми задними кромками.
1. Используется известный интеграл [1], связывающий потенциал скорости <Р со скосом ж.
где 5 — зона влияния, включающая часть крыла ^ и, как правило, часть пелены возмущений 5П (рис. 1). Здесь применены следующие обозначения: х, г—безразмерные координаты крыла (х — ось симметрии или бортовая хорда); ч\ — применяются вместо х, г при интегрировании, в этом случае х, г — координаты точки; х', г'— размерные координаты, связанные с безразмерными через характерный размер Ь—он же габарит крыла в направлении оси х\ ЭИ — число Струхаля, юЬ/У (ш — частота колебаний, V — скорость); р = = 8ЬМ/(М3—1) —число Струхаля с учетом сжимаемости; М — число М, М=1//1/зв (1/3в — скорость звука); х0 = £ — т], г0 = г — -ц — коорди-наты, связанные с зоной интегрирования; /? = V— Р2 2$ — гиперболический радиус (р3 = М2 — 1).
В. Г. Буньков
ср = — — [ \ но ехр (—Шрх0) сое (рЯ)1Я (14 с1г[, (1)
те с
Безразмерный скос на крыле равен:
уу — ~ 1БЬ/,
дх
где/—заданная форма колебаний для деформации у=Ь$(х, г)ехр(гшг!).
Скос вне крыла неизвестен и определяется из интегрального уравнения (1) при условии обращения в нуль потенциала 9 впереди и сбоку крыла. Пелена сзади крыла не участвует в расчете, так как мы ограничиваемся случаем сверхзвуковых задних кромок.
Коэффициент давления с'(с'=ср14) равен:
— + IБЬ ср. дх
(3)
Он связывает разность давлений Др с плотностью воздуха и скоростью следующим образом:
'*с'реш.
<Хж\
/ ‘Г,
г
г (х)
(4)
Шзона
Рис.
В Зоны нрыла.(1) Рис. 2
2. Идея разбивки крыла и пелены на конечные элементы заимствована из работы [2], однако вместо треугольного элемента за основу взят трапециевидный. Мы ограничиваемся лишь уравнением для потенциала скоростей (1) в отличие от работы [2], в которой используется также и уравнение, связывающее возмущение скорости с давлением на крыле. Предлагается следующая схема крыла.
Консоль крыла состоит из I трапециевидных зон Каждая зона может состоять из трех участков: переднего, среднего и заднего. Обязательным является только передний участок, а средний и задний могут отсутствовать. Каждая /г-я зона делится на гк вертикальных полос, а полосы делятся на элементы. Полосы в средних и задних участках делятся просто: tl элементов во всех средних участках и соответственно — в задних. Передние же участки могут состоять из смеси трапециевидных и треугольных элементов, для чего на левом краю к-го переднего участка задается элементов, а на правом 8к. При = получаются только трапециевидные элементы. При сначала определяется
число узлов на промежуточных вертикалях участка (приближенно по линейному закону), а затем все узлы соединяются перемычками, причем перемычки от лишних узлов стягиваются в один передний узел. Соответственно у передней кромки крыла образуются один или несколько треугольных элементов (рис. 2).
Шаг по ширине полос и по расстоянию между узлами на вертикалях может варьироваться по геометрической прогрессии.
Задние участки предусмотрены на случай задания более мелких элементов на элеронах крыла.
Пелена возмущений разбивается на конечные элементы по тому же принципу' но с сохранением лишь передних участков. Кроме того, для одиночного крыла, когда не нужна пелена между крылом и оперением, пелену оказалось удобнее разбивать на полосы, параллельные линиям Маха (см. рис. 2).
Потенциал ср и деформации / на поверхности трапециевидного конечного элемента (рис. 3) предлагается определять по их значениям в углах элемента на основе следующих гипотез.
Рз
Значения потенциала в четырех углах обозначим <р2, <?в, <р4. Потенциал на любой стороне меняется по линейному закону. Также по линейному закону изменяется потенциал на любой хорде Д, параллельной основаниям и Д2. Такому же закону подчиняется и деформация / в элементе.
При такой схематизации давление р (разность давлений) окажется линейно меняющимся на каждой хорде элемента, а стационарная часть (вещественная составляющая) — постоянной. Чтобы привести получившееся давление к сосредоточенным силам Ри Р2, Р3, Р4 в углах элемента, рассуждаем следующим образом. Физический смысл заданного закона деформации соответствует жестким хордам и отсутствию жесткости на кручение элемента. В таком случае все давление в окрестности каждой хорды передается в виде равнодействующих сил на передний и задний концы хорды, а получающиеся при этом погонные нагрузки на передней и задней стороне элемента дают силы в углах. Связь сил в углах с потенциалом в этих углах можно записать в матричном виде:
-22-1 1 -2 2—11 -11—2 2
-2 2
/2 1 2 1 12 12 2 16 3 12 3 6
(Л /
р2 й \
р3 ■ 12
Р4 \
+
і<і БИ
1 1
72
6 3 2 Г
3 6 1 2
2 1 2 1
1 2 1 2
+
(5)
где (I— ширина конечного элемента.
3. Интеграл (1) вычисляется с помощью интегрирующей сетки с заранее затабулированными значениями интеграла в нескольких тысячах узлов (см. рис. 1).
Построим прямоугольную сетку в характеристических координатах: и = л;0 — $г0, V = х0 + р20, для чего зададим положение узлов сетки по формуле:
(6)
где 0-<а<Л —коэффициент, регулирующий сгущение сетки вблизи границы с особенностями (и = 0; ю = 0); Л^= 40-ь60 — число узлов сетки на отрезке Узлов должно быть взято достаточно,
чтобы покрыть треугольник « + ^<2, т. е. Л^12^2Л''.
В интеграле (1) в подынтегральной части выделим гиперболический радиус К = угиъ, после чего интеграл приобретает вид
где функция /(и, V) включает в себя и скос т, и часть ядра, и множители, в частности множитель 1/2(3 при переходе от скй-ц к йи (IV.
Будем полагать, что функция /(и, V) задана в узлах сетки: Лк~7(Щ> г"» к = 0, 1, 2, . . а в каждой прямоугольной клетке
меняется по линейному закону как вдоль оси и, так и вдоль оси V. Тогда интеграл (7), если область интегрирования условно дополнить до прямоугольника с пустыми узлами, берется в явном виде:
Расшифровывая величину /ш получаем требуемую формулу, выражающую потенциал через значение скосов чю1к в узлах г, к сетки:
Важно отметить три особенности практической реализации этой формулы:
1) для симметричного крыла в расчете должна участвовать только половина крыла (например, правая консоль). Весь интеграл представляется в виде суммы двух частей: одна часть — из точки (я, г) с обрезанием зоны интегрирования слева от оси х, другая часть — из симметричной точки (л:, —£), также с обрезанием слева от оси х (рис. 4). Для симметричных деформаций эти части складываются, а для антисимметричных вычитаются.
*=я
У иу
(7)
т п
(8)
1=0 к = 0
где
с0=4 Уиг;з,
2ик+\ У чк+1 + (ик-Зик+1) Уик1
1
“А + 1 - “к
(10)
I, к
где
Щ ехр (-№рх0) соэ (рЯ).
2) Суммирование узлов следует вести вдоль ЛИНИЙ V = const, начиная с -у = 0. При этом, если линия v = const отсекается снизу границей зоны интегрирования при 0 (ось х на рис. 4), то необходимо учитывать дробную долю для нижнего узла, т. е. определять, какая доля отсеклась от первой клетки, точнее, от отрезка линии.
3) Потенциал от крыла и от пелены вычисляется отдельно, независимо от того, где лежит центр (х, г) интегрирования: на крыле или на пелене. При этом у передней и боковой кромки крыла образуются „рваные“ клетки. В таких клетках также учитывается дробная доля пропорционально доле, отсекаемой на очередном отрезке линии ■иссопе! границей зоны интегрирования (кромкой крыла). Причем, если рассекается отрезок между и0 = 0 и и1г т. е. отрезок, содержащий интеграл с особенностью, то учет доли должен быть особый, например, для ./N/ = 40 и а = 0,5 получается поправочный множитель й0 = 1,51/5 для нулевого узла и ^о = 2,21 — — 1,811/5 для 1-го узла (т. е. со стороны крыла и со стороны пелены), где 5 = «/«! — отсекаемая доля.
4. Задача о расчете потенциалов методом конечных элементов сводится к отысканию потенциала <р в п точках, являющихся углами конечных элементов. Обозначим неизвестные потенциалы в этих п точках через вектор Попутно приходится определять неизвестные скосы в т' = т — т0 точках пелены, где т — общее число точек на пелене, а т0 — число точек на передней кромке пелены, где скос нулевой. Скосы в т' точках пелены обозначим вектором УУт', а потенциал ут>.
После интегрирования описанным выше способом получаем систему уравнений:
?п = АУп + В]Ут•; |
<?т, = СУп + В\Ут' = 0, |
гХе А, И — квадратные комплексные матрицы, а В, С — прямоугольные. На крыле вместо скосов оказалось удобнее подставлять непосредственно вектор деформаций Уп в соответствии с формулой (2) и свойствами конечных элементов.
Решением системы являются матрицы Р и
(7* = Л+.£/=■, (12)
связывающие потенциал и скосы УУт' с деформациями крыла Уп:
<?n = G9Yn, Wmi =FYn.
(13)
Ряд строк в матрицах А, В, а следовательно, и в С?9, заведомо
равны нулю. Это строки, соответствующие точкам, лежащим на
передней и боковой кромках крыла. Из матрицы (39 может быть получена так называемая матрица аэродинамических коэффициентов влияния С#, связывающая сосредоченные силы в п точках с деформациями в этих точках:
Я=ОкУп. (14)
Для получения необходимо сначала выразить через срп;
(15)
где матрица перехода составляется на основе формулы (5), применяемой поэлементно, а далее:
= Ов, . (16)
Аналогично получается матрица давлений.
Матрица влияния может быть непосредственно использована в задачах аэроупругости, решаемых методом конечных элементов при условии, что аэродинамические элементы совпадут с элементами упругой схемы. При использовании других методов, например, метода многочленов [3], так называемые матрицы аэродинамической жесткости В и аэродинамического демпфирования И получаются по следующим формулам:
5 = _ТХТСЛ1А; Я = -Т'ЛГ1 <?*,*; ч = 2РУЧ; т'=у§^ (17)
где Оц 1, Ода—вещественная и мнимая части матрицы вц, а ЛГ—матрица перехода от, вектора обобщенных координат к вектору деформаций:
У а = XV. (18)
5. Решение описываемой задачи на ЭВМ БЭСМ-6 было разделено на ряд отдельных блоков (стандартных программ): автоматическая разбивка крыла и пелены на конечные элементы, сортировка узлов, расчет таблицы для интегрирующей сетки, составление основных матриц системы, перемножение и обращение комплексных матриц большого порядка и др. Наибольшее число узлов на крыле равно «=160, на пелене также т.' =160 (с пустыми узлами т = т' + т0 = = 192). Время расчета для п—36, 94, 159 получилось соответственно 1, 5, 23 мин.
Приведем характерные результаты.
Крыло специфической формы (правая консоль показана на рис. 5). Сравниваются точные и расчетные распределения давления в трех сечениях при М = 'У"2. В соответствии с особенностью закладываемых в расчет конечных элементов расчетные эпюры давления оказываются ломаными. Эти ломаные хорошо накладываются на точные кривые. В этом расчете принимали А/—40 (N=60 практически улучшений не дает), а = 0,5, л =121, т' — 66. Время расчета 8 мин.
Прямоугольное крыло с удлинением Х = 3. Исследовалась возможность появления отрицательного демпфирования при вращении вокруг оси, параллельной передней кромке. Не занимая места рисунком, отметим, что отрицательное демпфирование получилось для БЬ = 0; 0,2; 0,4 соответственно при М<1,38, 1,36, 1,28, а при вращении вокруг передней кромки: М<1,3, 1,26, 1,21. При 511>0,45 оно пропадает совсем.
Как показывают расчеты, отрицательное демпфирование имеет место для достаточно малых чисел Струхаля (5Ь<;0,5). Это согласуется с результатами расчетов для колеблющегося профиля по нелинейной теории [4]. Диапазон чисел М, в котором возможно отрицательное демпфирование, сокращается с уменьшением удлинения крыла. Так, для треугольных крыльев со стреловидностью больше 60° или прямоугольных с удлинением, меньшим >. = 1,5, отрицательное демпфирование по расчету пропадает. Однако это справедливо только для недеформируе-мых крыльев, при колебаниях же упругих крыльев картина может измениться, так как отдельные участки крыла при своей деформации могут имитировать крыло большого удлинения, например,элероны при их отклонении.
' Исследовалось влияние чисел М и БИ на флаттерные характеристики. В качестве примера была рассчитана без-рулевая симметричная форма флаттера самолета с крылом малого удлинения.
На рис. 6 показана зависимость критической скорости флаттера от М (скорость выражена в индикаторном виде, т. е. 1/инд= 1/звМ/Д, где Д —получившаяся относительная плотность при
Рис. 6
заданной истинной скорости V=V3BM-, Vwm характеризует скоростной напор). По числу Струхаля было взято два варианта: Sh = 0 (0,001) nSh=3(4TO в среднем соответствует истинному числу Струхаля). Эти два варианта оказались качественно противоположными: при Sh=0 скорость флаттера падает с уменьшением числа М и становится нулевой при М= 1,1 (отрицательное демпфирование), а при Sh=3 она постоянна и даже несколько подрастает. Показан разброс результатов при разном числе точек; так, результат при п — 94, т' = 92 й п =159, т! ==155 отличается на 4% по 1/„нд. ф. Интересно, что для получения Су прямоугольного крыла с точностью 1% достаточно взять я = 36.
ЛИТЕРАТУРА
1. Красильщикова Е. А. Крыло конечного размаха в сжимаемом потоке. М. —Л., ГИТТЛ, 1952.
2. Hashimoto М., Washizy К. Application of the finite element technique combined with the collocation method to low supersonic lifting surface problems. International Journal Numerical Methods in Engineering, vol. 11, N 10, 1977.
3. Буньков В. Г. Учет деформаций сдвига при расчете колебаний крыла малого удлинения методом многочленов. „Ученые записки ЦАГИ“, т. Щ, № 4, 1972.
.4. Кузьмина С. И. Численное исследование колебаний профиля в трансзвуковом потоке газа. Труды ЦАГИ, вып. 1862, 1977.
Рукопись потупила 31\Ш 1980 г.
