Расчет устойчивости изгибно-крутильных колебаний крыла в дозвуковом потоке методом конечных элементов Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ И А Г И

Том Ш 1972 № 1

УДК 629.735 33.015.4:533.6.013.42

РАСЧЕТ УСТОЙЧИВОСТИ ИЗГИ БИО-КРУТИЛ ЬН ЫХ КОЛЕБАНИИ КРЫЛА В ДОЗВУКОВОМ ПОТОКЕ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ

В. П. Кандидов, Л. П. Ким

Рассматривается применение метода конечных элементов для расчета устойчивости изгибно-крутильных колебаний крыла самолета в дозвуковом потоке. Излагается способ, позволяющий повысить точность в описании упругих и инерционных свойств отдельного неоднородного элемента без увеличения числа степеней свободы. Рассмотрен пример.

В последнее время интенсивное развитие получил метод конечных элементов — эффективное средство расчета статики и динамики непрерывных упругих систем [I]. В литературе уделяется большое внимание вопросу сходимости [2) и оценке погрешности метода [3], (4]. Для повышения точности расчета поперечных колебаний стержней в работах |5] и [6] предлагаются неоднородные конечные элементы. В работе [7] для этих целей применяется искусственный прием, практическое использование которого, по-видимому, довольно ограничено.

В настоящей статье излагается способ, позволяющий повышать точность в описании упругих и инерционных свойств отдельного неоднородного элемента. Он используется при построении конечных элементов для расчета совместных изгибно-крутильных колебаний консольного стержня, а также для получения его матриц аэродинамической жесткости и демпфирования при исследовании колебаний в дозвуковом потоке.

!. Основные положения метода конечных элементов. В методе конечных элементов исследуемая система разбивается на части, каждая из которых заменяется некоторой моделью, называемой конечным элементом (КЭ). Функция деформации КЭ и {г, /) определяется в пространстве базисных функций:

и (г, I) =<)/ (г)•?(/),

где г—пространственная переменная, область изменения которой /?; у(/)—зек-тор координат функции деформации в пространстве базисных функций;

—►

ф (г) — вектор, компоненты которого являются базисными функциями.

Введя матрицу преобразования координат а, удобно перейти к вектору

<1 (/), координаты которого имеют простой физический смысл линейных и угловых

смещений узловых точек конструкции. Тогда функция деформации в новой

системе координат примет вид

и (г, 0 = <|| (г) а 1} (/). (1.1)

III

Существенно, что КЭ с функцией деформации (1.1), в отличие от соответствующего участка непрерывной системы, имеет конечное число координат, определяющих его положение и деформацию. Благодаря этому уравнения движения элемента можно представить в матричной форме, введя матрицы инерции и жесткости. Распределенные силы реакции, действующие на контуре КЭ, сводятся к эквивалентной системе сосредоточенных сил взаимодействия между элементами. Такая процедура может быть осуществлена либо простым усреднением реакции на отрезках контура [8]. либо с помощью вариационных методов [ 11, [9]. Последний способ, являющийся наиболее общим, используется в настоящей работе.

Используя принцип виртуальных перемещений, получим уравнение движения н-го КЭ в виде

Мп ч„ (0 ~ Кп д„ (/) - ?„ (/) - 0„ (/) = О, где М„ — матрица инерции,

М„ = ( |1 (г) а|/ (г) у (г) а' Лг\

Я

Кп — матрица жесткости.

Кп = [ е„(/-)<*г;

/?

Рп(1) — система сосредоточенных сил, эквивалентная распределенной нагрузке плотности /(г, О.

МО- \пг, 1)аЬ(г)<1г, (1.5)

I

0л(0— вектор обобщенных сил реакции, соответствующий вектору обобщенных координат <?(/);

¡л (г) — погонная инерционная плотность; е„(г)—матрица плотности потенциальной энергии.

Для конкретного КЭ с произвольными параметрами матрицы К„. Мп и вектор Р„(0 вычисляются по (1.3) —(1.5) методом численного интегрирования. Однако, аппроксимируя параметры КЭ в области /? какими-либо несложными функциями, можно получить приближенные аналитические выражения для элементов матриц Кп, Мп и вектора р„ (<) КЭ в некотором пространстве <|/ (г).

В этом случае вычисление К,„ Мп и /г„(/) значительно упрощается.

Обычно принято |1] мри определении матриц К„ и Мп использовать одно

и то же пространство ф(г). Это является неоправданным ограничением. Действительно, в каждое из выражений (1.3) и (1.4) может быть подставлена своя,

независимая система базисных функций. Обозначим через (г) пространство

функций для определения Кп и через Ь'* (г) — для М„. В этом случае функция

деформации и* (г, () и перемещения и1 (г, <) не совпадают. Однако вектор <?„(<) должен быть одним и тем же. Используя разные системы базисных функций

(г) и (/')* можно независимо варьировать точность описания упругих и инерционных свойств отдельного КЭ с помощью матриц Кп и Мп. Выбор пространства базисных функций для каждой из матриц определяется, с одной стороны, стремлением достичь более высокой точности, с другой, — уменьшить число степеней свободы КЭ для сокращения объема вычислений.

Для соседних элементов в общих узловых точках имеют место условия

сопряжения: геометрические, связывающие координаты вектора q„ (/), и динамические. накладываемые на вектор обобщенных сил РпМ). В узлах, лежащих на границе системы, накладываются условия, соответствующие конкретному виду закрепления ее краев.

Совокупность уравнений движения (1.2), записанных для всех КЭ, и условий сопряжения во всех узлах образует систему уравнений, которая описывает поведение исходной системы, представленной набором из КЭ. В общем случае она состоит из алгебраических и обыкновенных дифференциальных уравнений второго порядка. При этом число дифференциальных уравнений равно общему

(1.2)

(1.3)

(1.4)

числу степеней свободы набора из КЭ. Точность описания исходной системы возрастает как при увеличении числа КЭ в наборе |4), так и при улучшении аппроксимации свойств отдельного участка системы соответствующим КЭ [5].

2. Конечный элемент для расчета изгибно-крутильных колебаний. Задача изгибно-крутильных колебаний стержня рассматривается при обычных предположениях, когда можно считать, что упругая связь между деформациями изгиба и кручения отсутствует и совместно с колебаниями кручения возникают поперечные колебания только в плоскости м>0; (фиг. I). Смешение в этой плоскости Судет

и («, У, О = «(«, О -I- У9 (5. /). (2.1)

Положение отдельного КЭ определим шестью координатами, в качестве которых удобно взять следующие величины:

?’(0 =

<?о(0. МО. «МО. т (0. 14™

<7;

с-.о

(О

. 1~\ (0І. (2.2) дг 1с=/ |

Здесь индексом .0“ обозначены значения, взятые на левом конце элемента при 5 = 0, индексом / —на правом при г; — 1.

йГ У Середина нрыла \№> <№>■*> А,*

Щ

0 Ып / 7*’”"’"

* 1 / 1

Ось местности

Фиг.

Вектор обобщенных сил реакции КЭ, соответствующий введенному вектору (2.2), будет иметь вид

ма т м, (/)

МО. Г, и). Р0(Ц. Р,(0.

(0 =

I

I

(2.3)

где Т{1)—закручивающий момент относительно О;, Р{1) — перерезывающая сила, М (/) — изгибающий момент относительно ОУ.

Рассмотрим деформацию КЭ в пространстве четырех базисных функций:

— / £ с* и ч

*'(«)-(і. -Ь 1Г- 1г)- (2-4)

Пусть ?(;, /) линейно зависит от ;, тогда матрицы преобразования координат будут:

(2.5)

Записав выражение для матрицы плотности потенциальной энергии изгиб-но-крутильиой деформации, найдем по (1.4) матрицу Кп• В частности, если жесткость на изгиб является линейной функцией /:./(;) —/3 ^/.| + 7_, _2_ ^ , то

' 0 0 1 0 0 0 1 0 0 0 0 0 '

0 0 0 0 1 0 -1 1 0 0 0 0

. а9 =

0 0 -3 3 -2 -1 0 0 0 0 0 0

. 0 0 2 -2 1 1 ! - . 0 0 и 0 0 0.

Кп =

Я

Симметрично

0 0

0 0

Ч- 2 їй 6 Ч 4 *

— 2 7.2 —6 Х| — 4 7.

"Ь У-2 2 7., 7-2

1 -г 3 7.

(2.6)

О/(:)</: — средняя жесткость на кручение.

Равенство нулю части элементов матрицы является следствием сделанных предположений об отсутствии упругой свнзи между изгибом и кручением.

Для получения матрицы Мп в (1.3) подставим (2.4) и матрицу преобразования, которая имеет вид [см. (2.1)]

а а“ +_у а'р.

При плотности, постоянной вдоль длины элемента ц (5, у) = [д(у), получим 1/3/ 1/6/ 7/20 5 3/20 ^ 1/20 ї —1/30 і

1/3 і 3/20« 7/20 і

13/35 сі 9/70 сі

13/35 сі

Ж„=/

Симметрично

1/30 ї 11/210 гі 13/420 сі 1/105 г/

—1/20 і —13/420 сі

— 11/210 сі

— 1/140 (і 1/105 сі

(2.7)

где с — ц (у)уг с1у— погонный момент инерции;

О

ь

5 = ^ (у) уЛу — погонный статический момент;

о

ь

^ = ] Р (У) ¿У — погонная плотность.

о

КЭ с матрицей М„ (2.7) имеет шесть колебательных степеней свободы. Можно сократить их до четырех, если для построения матрицы инерции взять следующие базисные функции

(2.8)

Полагая в этом случае линейное изменение параметров /(5) = /, +/2£//, 5 (?) = .9, -)-*2 ;//, с1 (5) = </, + с1.2 £//, найдем

/,4-1/4 /, 1/2/,+ 1/4/3 5,-1-1/4

/| --г 3/4 /о 1/2 5,4 1 /4 5о

г/, -)- I /4 с1г

М„ =

Симметрично

1/2»,+ 1/4«, 0 0

+ 3/4 0 0

1/2 ¿, + 1/4 0 0

а, + з/4 оо

о о

о

(2.9)

Четыре степени свободы имеет также КЭ, у которого инерционные параметры являются сосредоточенными. В этом случае

Мп — і

0

п

-’о

0

"о

Симметрично

0 0 0 "

0 0

0 0 0

сі/ 0 0

0 0

0

(2.10)

Модель из КЭ с матрицей жесткости (2.6) при *, = 0 и матрицей, инерции (2.10) является аналогом модели Хольцера — Майклстеда (10|. Эту модель можно наглядно представить в виде системы точечных масс, соединенных невесомыми однородными стержнями.

При сокращении числа степеней свободы КЭ безусловно уменьшается точность в аппроксимации динамических свойств. Однако, увеличивая число разбиений, можно достичь первоначальной или более высокой точности при одинаковом числе степеней свободы всего набора элементов в целом. Используя полученные матрицы жесткости и инерции для расчета конкретных систем, можно оценить точность той или иной аппроксимации.

В случае стационарного приближения теории несущей полосы [12] плотность нормальной аэродинамической нагрузки и аэродинамического момента определяется по известным формулам

/«-(?, 0--рЮ(5)йв>|-»«. 0+И<р(5. О 6(€)[1/2 —я(5)]у(с. /)} . (2.11)

Р (=, t)=-t>Vb*(Z) + b(l)

+

-ю(<. 0+^(5. /) +

1

2 0 (;)

а («)] 4(5) {

Н*. + 0.

(2.12)

где р— плотность воздуха, с® (;)—производная по о от коэффициента подъемной силы.

Полная система сил, эквивалентная всей аэродинамической нагрузке, действующей на элемент, будет иметь вид

Fn (') = - VD„ q„ (0 - V: В„ q„ (i).

(2.13)

где й„ — матрица аэродинамической жесткости, D,, — матрица аэродинамического демпфирования.

Если для получения аналитических выражений элементов матриц Вп и Dn взять вектор обобщенных координат (2 2) и систему базисных функций (2.8), то при линейном изменении параметров и с£(;) = const получим

0000 0000

Я„=2 cf, pt

-I

5

Ь1

20

Ь0 1 bt Ьй bi

4 "Г 12 12 12

Ьо , Ь[ bo , b[

~i2 ' 12 12 ■ 4

^ (<, + *,) + bj L U bf _La. Л a. + t'iA (a,+a„)-f -La. 30 °г 20 "зо 20

— з2 ^.ai + b0bl(, +x )+^_i 20 2 30 1 ’ -26“ 5

0 0

0 0

0000

0000

(2.14)

Л1 2 т "0. *2 * 2

Матрица й,, имеет похожую структуру. Полагая и Ь/ = Ь0, легко

перейти к однородному элементу, параметры которого постоянны.

3. Некоторые результаты исследований. В качестве примера рассматривались изгибно-крутильные колебания прямого крыла большого удлинения, параметры которого представлены на фиг. 2. Сначала были рассмотрены собственные колебания крыла в пустоте. Чтобы оценить сходимость метода при расчете кои-

кретных систем, крыло разбивалось на 6, 8, 12, 16 и 18 элементов. Решение матричных уравнений бсуществлялось на ЭЦВМ М-20 с использованием стандартных программ |11) для вычисления собственных значений матриц.

Результаты расчета представлены в таблице. Видно, что с увеличением числа степеней свободы имеет место сходимость, причем для однородного КЭ снизу, а для линейного сверху. Как было указано ранее, целесообразно применение матриц инерции с уменьшенным числом степеней свободы, поскольку, увеличивая N. можно получить такую же или даже большую точность (см. таблицу). Сравнение линейной и однородной аппроксимации параметров—в пользу неоднородных элементов, так как при том же числе степеней свободы они дают большую точность. Это позволяет с помощью неоднородных КЭ уменьшить без потери точности число степеней свободы при расчете и тем самым уменьшить время счета на ЭЦВМ.

Элемент Матрица инерции Число элементов Число степеней свободы Ш| ' 0>і “5 <■>0

(2-7) 6 18 31,00 119,9 156.9 310.5 352.7 604,6

(2-9) 8 16 30,88 120,7 156,3 300,3 354.6 547,0

Однородный (2-9) 12 24 31,11 121,3 157.2 305,1 354,9 553,5

(2-9) 16 32 31,45 122.4 157,3 310,7 355,0 558.8

(2-9) 18 36 31.52 122.6 157,4 312,3 356.2 559,5

(2-9) 6 12 32.34 130,2 159.8 335.5 370,5 592,2

(2-9) 8 16 31,94 124,0 158.4 323,6 362,1 517,9

Линейный (2-9) 12 24 31,70 123,2 157,8 318,2 358,5 558,1

(2-9) 16 32 31,59 123.0 157,8 317.9 357,9 563,9

(2-9) 18 36 31,59 122,9 157.7 316,7 357.7 560,4

При исследовании устойчивости изгибно-крутильиых колебаний крыла было использовано разбиение на 8 и 12 элементов. В первом случае модель имела 16 степеней свободы, во втором —24. Изменения параметров на элементе аппроксимировались линейными функциями.

На фиг. 3 представлены траектории корней для первых трех тонов при разбиении крыла на 8 и 12 элементов. Скорость потока V изменялась через 50 м/сек. Видно, что траектории корней на первом тоне для N =■ 8 и 12 практи-

— #=г — п I кручение со М [Тёк] ио-

20 300 501 Г

УСТ

Мл И

300_ . —

ш - 10 Ж изгиб

т 10 1-

7-

ш юо

I изгиб -

м

Л т -»І- 500

300 «00 500 V[м/сек]

Фиг. 3

Фиг. 4

чески совпадают. При скорости V, равной примерно 320 м/сек, происходит вырождение первого тона, и при И=460 м/сек наступает статическая неустойчивость типа дивергенции. Следовательно, в рассматриваемой области дозвуковых скоростей крыло остается устойчивым Некоторое различие в поведении траекторий корней на втором и третьем тонах при N = 8 и 12 можно объяснить тем, что модель из двенадцати КЭ является более точной.

Представляет интерес сравнить с приведенными данными результаты расчета по методу Бубнова — Галеркииа (13). При этом в качестве системы координатных функций используются собственные формы крыла » пустоте, полученные методом конечных элементов. Как и следовало ожидать, результат существенно зависел от того, какие формы взяты в качестве координатных функций. На фиг. 4 приведена зависимость частоты и> и коэффициента затухания 6 от скорости потока V. Если координатными функциями являются первая и третья формы или первая, вторая и третья формы, найденные при N — 12. то «(К) и о(У) совпадают с полученными по методу конечных элементов при том же N (на фиг. 4 кривые, отмеченные точками). Если в качестве координатных функций используется первая форма или первая и вторая формы, то результат значительно отличается от предыдущего (на фиг. 4 кривые, отмеченные крестиками). Такая же картина наблюдается, если за координатные функции взяты собственные формы крыла в пустоте при А/=8.

Из этого исследования можно сделать вывод, что статическая неустойчивость на первом топе возникает в результате взаимодействия первого изгибного и первого крутильного тонов колебаний. При расчете по методу Бубнова — Га-леркина в систему координатных функций должны быть включены формы колебаний этих тонов.

При расчете собственных колебаний одномерных систем методом конечных элементов с помощью неоднородных КЭ можно без потери точности уменьшить число степеней свободы системы и тем самым сократить объем вычислительной работы.

Метод конечных элементов может с успехом использоваться для исследования неконсервативных задач аэроупругости. Описанная методика позволяет легко получить матрицы аэродинамической жесткости и аэродинамического демпфирования консольного крыла по заданным распределениям погонной плотности нормальной нагрузки и аэродинамического момента.

Метод Бубнова—Галеркииа. где в качестве координатных функций взяты формы собственных колебаний в пустоте, полученные методом конечных элементов, дает такие же результаты, что и метод конечных элементов, существенно уменьшая время счета на ЭЦВМ.

ЛИТЕРАТУРА

I. Zienkiewicz О. С. The finite element method in structural and continum mechanics. London, New York, 1967.

2. Eduardo R. de Or antes e Oliveira. Theoretical foundations of the finite element method. Ini. J of Solids and Structures, v. 4,

No 10. 1968.

3. L e с k I • F. A., L 1 n d b erg Q. M. The effect of lumped parameters on beam frequencies. The Aeron. Quart., 14, No 3, 1963.

4. Кандидов В. П., К и м Л. П. Дискретная модель для исследования колебаний балок. Вестник МГУ, сер. .Физика, астрономия*,

Л* 5, 1968

5. Krishna Murty А. V., Prabhakaran К R- Vibration of tapered cantilever beams and shafts. The Aeron. Quart., 20, part 2, 1969.

6. Kan war K. Kapur. Vibration of a Timoshenko beam using finite—element approach. J. of the Acousl. Soc. America, v. 40, No 5, 1966.

7. Cohen E., Me Callion H. Improved deieimination functions

for the finite elements analysis of beam systems. Int. J. of Numerical

methods for engrs., v. 1, No 2, 1969.

8. Кандидов В. П. Приближенный расчет неоднородных пластин путем разбиения на элементы. Вестник МГУ, сер. .Математика. механика“, № 4, Г964.

9. Me lose h R. J. A stiffness matrix for the analysis of thin plates

in bending. JASS, v 28. No I, 1961.

JO. Mlnhinick i. T. The theoretical determination of normal modes and frequencies of vibration. ARC, RaeM, No 3039. 1957.

11. Буньков В. Г. Программа расчета собственных колебаний неконсервативных систем. Труды ЦАГИ, вып. 977. 1965.

12. Бисплннгхофф Р. Л., Эшли X., Халфмэн Р. Л. Аэроупругость. М., Изд иностр. лит., 1958.

13. Стрелков С. П. Введение в теорию колебаний М.. .Наука",

1969.

Рукопись поступила 15/ VI 1971 г.