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

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

То м XV 19 84

№ 5

УДК 629.735.33.018.4 : 620.178.3

УЧЕТ НАЧАЛЬНЫХ УСИЛИЙ В СТАТИЧЕСКИХ И ДИНАМИЧЕСКИХ РАСЧЕТАХ КОНСТРУКЦИЙ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ

В. П. Агапов, В. Д. Ильичев, В. А. Коротков, А. В. Стрелин

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

Явление статической потери устойчивости конструкций (в том числе и авиационных) и их элементов достаточно хорошо изучено (см., например, [1—3]). При использовании метода конечных элементов критический параметр нагрузки находится из известного уравнения (см. И)*:

[/С-ХК,]{£/}=0,

где [/(] — матрица жесткости, X [Ка\ — матрица начальных напряжений**, причем [АТ®] соответствует единичному значению параметра нагрузки X, {6'} — вектор узловых перемещений, определяющий форму потери устойчивости.

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

* Здесь и далее матрицы обозначаются квадратыми скобками, а векторы — фигурными. ,й

** Эту матрицу называют также матрицей геометрической [21 или дифференциальной [4] жесткости.

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

Рассматриваются следущие задачи.

1. Статический расчет конструкций с учетом начальных усилий с использованием уравнения [6]

[К-'кК,]{Ъ\ = {Р), (1)

где {8} —вектор узловых перемещений.

Предполагается, что матрица [ Ка ] не зависит от вектора нагрузки {Р}, стоящего в правой части уравнения (1). Такая ситуация возникает при продольно-поперечном изгибе стрежней, пластин и в других аналогичных случаях.

2. Определение частот и форм собственных колебаний с учетом начальных усилий с помощью уравнения

[/С-Ж.]{84} = 2?[Л1]{8/}, /= 1, 2, . . ., V,

где [М] — инерционная матрица, йг — частота /-го тона собственных колебаний, [6;] — вектор узловых перемещений, характеризующий форму ¿-го тона собственных колебаний.

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

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

[/С-Х^]{Й}-62[^]{5} = {Я|,

где 0 — частота вынуждающих сил, {Р} и {6} — векторы амплитуд узловых сил и перемещений соответственно.

Для решения перечисленных задач составлена программа для ЭВМ. Программа написана на языке ФОРТРАН и названа ПУСК-3 (Прочность, Устойчивость и Колебания). Реализованный в программе метод расчета на устойчивость описан в работе [5].

В библиотеку конечных элементов программы включены: 1) стержень (элемент фермы); 2) балка (элемент рамы); 3) плоский несогласованный элемент комбинированного типа (комбинируются плосконапряженное и изгибное состояния); 4) плоский согласованный элемент комбинированного типа.

Отметим, что в теории метода конечных элементов принято называть согласованными такие плоские конечные элементы, которые при изгибе пластины обеспечивают совместность как прогибов, так и углов поворота для соседних элементов не только в узловых точках, но и вдоль границ. Однако построены и применяются на практике плоские элементы, которые не обеспечивают совместность углов поворота вдоль границ. Такие элементы принято называть несогласованными [2]. В работе '[2] показано, что плоские несогласованные конечные элементы в некоторых случаях обеспечивают более высокую точность расчетов, чем согласованные, поэтому в библиотеке конечных элементов целесообразно иметь элементы обоих типов.

Отметим также, что все плоские элементы с различными порядками аппроксимаций мембранных и изгибных перемещений в расчетах оболочек являются несогласованными ¡[2] как по перемещениям, так и по поворотам вдоль границ. Совместность перемещений и углов поворота в узловых точках обеспечивается, однако, в любом случае.

Для вычисления матриц жесткости и массы элементов использовались известные алгоритмы (см., например, {2, 4]). Матрица начальных напряжений \К1\ для всех элементов получалась с помощью формулы

(2)

где V—объем конечного элемента.

Матрица [б] в формуле (2) связывает вектор поворотов {ш} и вектор узловых перемещений {8е} элемента по формуле

{«•>}-= [О] {*}•

В свою очередь, вектор поворотов состоит из компонентов <ах, о>у и о)г, которые определяются следующим образом [7]:

__ 1 ( дт дк \ _ _ 1 I ди дни \ __ 1 / дг» ди \

<°х 2 \ ду дг I ’ 2 \ дг дх ) ’ (°г 2 \ дх ду }'

Матрица [5] в формуле (2) состоит из компонентов тензора начальных напряжений и имеет следующий вид:

[5] =

/ 0 (<3уу

о \

°ху

/О . 0 \

\?хх -Г °гг)

О

— ауг

/ О (°хх-

0 \

аУУ/

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

Рис. 1

Схема элемента приведена на рис. 1. При получении всех матриц элемента аппроксимирующие функции для мембранных (и, г>) и изгибных (до) перемещений приняты в виде

(3)

И = а 1 + а2 + а3у,

V — а4 а5 Л + а6_у,

® = ?1+?2.* + ?33' + ?4^ + Яъ ху + дву2 + д, хъ + ?8 ху"- + у3. (4)

Компоненты напряжений в начальном положении

о

Зхх '■

-О

М°гг

о _~о К 2

-- Зуу ■ ^

■>ху ■

''ху

К*

6 — «Ученые записки» № 5

8]

о о о п

О 22 <Зуг — <Згх - и,

где <з°х, о°уу, — мембранные напряжения; М% — изгибаю-

щие и крутящие моменты, а /—момент инерции пластинки. Следовательно, матрица [5] в данном случае

[5]

-о Ь М%г

Оуу I ) 1 а*У I

-о М%г\ (5°

аху / 1 ахх I

О

О

О

М°хг

I -уу I ,

Подставив (5) в (2) и выполнив интегрирование по г с учетом того, что повороты со*, щ, <мг, а следовательно, и элементы матрицы [й] от 2 не зависят, получим

К Л

(5)

[*■'] = а|[0]т[5] \G\dF,

(6)

где Л и ственно,

толщина и площадь конечного элемента соответ-

Г Оуу -о аху 0

[5] = -о аху ахх 0 (7)

0 0 “О , “0 ахх ~Г ауу .

Рассмотрим отдельно получение матрицы [Ке„\ для изгибных и мембранных движений.

При изгибных движениях перемещение т определяется формулой (4), а перемещения и и V [8]

Тогда

дно дни

и = г~дх ’ V = ■ — 2~л— ■ ду

1 ( дт дь \ дт

х 5=1 2 \"37 ~ дг 1 ~~~д^у

1 ( ди дш> \ дт

У = 2 \ дг дх / дх

1 1 ди ди ) 1 = 0.

г = 2 \ дх ду )

(8)

Используя (4) и (8), находим

М = 1С„]{?},

где {?} = |_<7,<72 • . . ?9_Г*>

О 1 0 2 х у

[Оа

0 0 1 Ю О О

О

0x2у ООО

З*2

о

о

у2 0 '

Злу 3у* О 0 .

(9)

* Здесь и далее угловыми скобками обозначается вектор-строка, а символ вверху означает операцию транспонирования'.

Вектор обобщенных координат {<7} элемента связан с вектором узловых перемещений {8и} зависимостью

{8и} = [Л„]{<7}, (И)

где

{8и} = [»1 ср,, <р, 1 ю2 ?х2<?у2 <?х з <ру з]т, (12)

а [Ли] есть числовая матрица перехода от обобщенных координат к физическим [4]. В формуле (12) ъ<?х1, <?у1 представляют собой прогиб и углы поворота вокруг осей х и у I-го узла конечного элемента (¿=1, 2, 3). Из (11) следует, что

{9} = ИиГ‘{85}. (13)

Подставляя (13) в (9), получаем

М = [0„]Ии1~1{8и}. (14)

На основании (6) и (14)

\&] = ({АИ]-'У (й | [ОиГ [5] [Ои](15)

Обозначим

1^,9] = л|[ОиГ[5][Ои]^. (16)

V

Подставляя в формулу (16) соотношения (7) и (10) и выполняя интегрирование, получаем симметрическую матрицу начальных напряжений для изгибных перемещений треугольного элемента в обобщенных координатах (выписаны лишь главная диагональ и расположенные выше нее члены), Причем Ох = Охх, Оу = Оуу, X — Оху'-

[Ка, ,

о о

ах 100

0 0

Чо 2°* 1\ 0

о <> о 2т/10

43г/20

0 0 0 0 0 -

ах А>1 + + ^10 2 т/0, Зах1ы 3Л^02 + + 2-с/ц 3-/02

. *А>1 + + °у Ло 2ву/01 Зт/20 ХА)2 + +2зу /и ЗЗу /02

2(°*/ц+ 4~ х4о) 4-с/ц КАо 2(аЛА2 + +2х/21) 6^12

аХ Лй + +2х/ц + + °у^20 2(*/<в + +ау Л1) 3(°л-/21 + +т^зо) А)3 + + 31/12+' +2йу 1п з (т/оз + + °у Л 2)

4зу /02 6^21 2 (т/оз Ч- + 2 оу/12) ^03

1 ■<Ь 9®* Ло 3(зл/22 + +2т/31) °ЛГ Ли + + 4т/13+ + 4ау /22 9х/22 3 (1/04 + + 2®у/13) 9зу /04 _

(17)

В матрице (17)

Переход к физическим координатам (узловым перемещениям) осуществляется на основании (15) и (16) с помощью соотношения

[К”1 = [Л^1]1 <?]

Если движения носят мембранный характер, то на основании (3)

1 I дтю

2 \ ду

1 ( ди

2 \ дг

1 / ди ~2 \1ЬГ

дно

~дх

ди

~(>У

= о,

1

(«5 - Кз)-

(18)

Соотношения (18) в матричной форме записываются в виде

Н = [0«]{«}, (19)

где {а} —вектор обобщенных координат, равный в данном случае

{«} =[<*! а2 а3 а4 а5 а6 ]т,

а

*0 0 0 0 0 0"

[С?м] = 0 0 О ООО. (20)

1.0 0 -0,5 0 0,5 0.

Вектор обобщенных координат {а} элемента связан с вектором узловых перемещений {8^} зависимостью

{&м}=Им]{а}, (21)

где

{8м} = [«1 г>1 м2 Щ «з ‘»зГ. (22)

а Им] — числовая матрица перехода от обобщенных координат к физическим [2]. В формуле (22) щ, V* представляют собой перемещения

г-го узла конечного элемента (¿=1, 2, 3) в направлении осей х и у

соответственно.

Из (21) получаем:

{«} = Им]-Ч8м}.

Подставляя (22) в (19), находим

К} = [0„1ИМГ ФЬ

На основании (6) и (24)

И = (Им]-1)7 (а | [Ом]т [5] [Ом] ^ ] [Д,]'1-

Обозначим

•*и

[л“а]=,л|[0мг[5][0м]^.

(23)

(24)

(25)

(26)

Подставляя в (26) соотношения (7) и (20) и выполняя интегрирование, находим матрицу начальных напряжений для мембранных перемещений треугольного элемента в обобщенных координатах

0 0 о

о о

0,25®*+;%)

0

0

0

0

о

о

- 0,25(5“*

0,25 (с

°°уу)

-О ч

вуу)

о

о

О

О

О

О

А/\

Переход к физическим координатам осуществляется на основании (24) с помощью матрицы (Им]~Т:

[к“]=ал,]-т”«] и«]

1-і

Матрица начальных напряжений для комбинированного элемента пластинки находится суммированием матриц [А'“] и [АГ“]. Матрицы начальных напряжений для остальных типов конечных элементов были получены аналогичным образом.

Перейдем к рассмотрению результатов расчета подкрепленных тонкостенных конструкций по программе ПУСК-3.

Значительное число задач было решено для проверки точности и оценки сходимости результатов, получаемых с использованием различных конечных элементов. Приведем некоторые ИЗ НИХ.

Для оценки результатов, получаемых с помощью согласованных и несогласованных конечных элементов при расчете оболочечных конструкций, была рассчитана хорошо изученная в теоретическом отношении (см., например, [9]) замкнутая круговая цилиндрическая оболочка, испытывающая действие внешнего давления (рис. 2, а). Исследовалась форма собственных колебаний с шестью волнами по окружности. Это позволило ограничиться в расчетах рассмотрением четвертой части оболочки.

В

Расчетная схема оболочки приведена на рис. 2, б. Вдоль образующей АВ исключены перемещения по оси г и повороты вокруг оси у, а вдоль образующей СО — перемещения по оси х и повороты вокруг оси у. Расчеты проводились при сетке узлов 4X8 и 9X14 при использовании как согласованных, так и несогласованных элементов.

Типичная для данной оболочки зависимость частоты О основного тона собственных колебаний от параметра нагрузки д в диапазоне от нуля до минимального критического значения, полученная с использованием несогласованых треугольных элементов при сетке узлов 9X14, приведена ниже на рис. 2,0. По программе ПУСК-3 найдена частота собственных колебаний ненагруженной оболочки 2? = 1622 с-1. По графику рис. 2, в находим верхнюю критическую нагрузку ¿7В= 1,33 мПа. Теоретические значения Я° и qъ равны* 1633 с-1 и 1,42 мПа соответ ственно. Расхождение составляет 0,6% для частоты собственых колебаний и 7,04%! для критического параметра. При сетке узлов 4X8 несогласованные элементы дали такие результаты: 2? = 1870 с-1; qъ =

= 1,65 мПа.

При использовании согласованных конечных элементов были получены следующие результаты: при сетке узлов 4X8 2? = 1959 с-1, дв=1,95 мПа и при сетке узлов 9X14 2? =1708 с-1, ^в=1,61 мПа. Следует отметить, что несогласованные элементы дают оценку частот собственных колебаний и критического параметра снизу, в то время как согласованные — сверху.

Оболочка, показанная на рис. 2, а, рассчитывалась на вынужденные колебания при действии вибрационной силы приклады-

ваемой в дополнение к статическому внешнему давлению (рис. 3, а). На рис. 3, б показаны зависимости амплитуды перемещения их точки приложения а вибрационной силы в направлении оси х от частоты 0 этой силы при различных значениях внешнего давления

По программе ПУСК-3 исследовалось изменение частот и форм первых пяти тонов собственных колебаний лопасти, изготовленной из стеклопластика, в зависимости от скорости ее вращения. Общий вид лопасти с указанием основных размеров приведен на рис. 4, а. Значения частот собственных колебаний лопасти при различных скоростях вращения приведены в табл. 1. Из табл. 1 видно, что все частоты монотонно возрастают при увеличении скорости вращения лопасти, причем наиболее интенсивно увеличивается частота первого тона собственных колебаний. В табл. 1 приведены также для сравнения частоты собственных колебаний лопасти, найденные экспериментальным путем (указаны в круглых скобках). Наибольшее расхождение (~13%) наблюдается для четвертого тона собственных колебаний.

Анализ форм собственных колебаний лопасти показал, что в исследованном диапазоне скоростей вращения эти формы практически не изменяются. Узловые линии форм собственных колебаний для неподвижной лопасти приведены на рис. 4,6.

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

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

* Теоретические значения определялись по формулам, приведенным в работе [9], при параметрах оболочки, указанных на рис. 2, а.

N4)

у

А т 1073 1200

Сечение по А-к

а) „

3-и. тон 2-и тон

5-й тон 4-й тон

б)

Рис. 4

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

та б л и ц а 1

Число оббротов лопасти в минуту Частоты собственных колебаний, Гц

Л /, Л Л

0 55,8 125,8 165,0 195,4 304,2

— (54,0) — (169,0) (225,0) (320,0)

700 57,3 >■ 126,5 166,5 196,1 306,2

1300 60.6 128,2 170,5 197,7 310,8

1900 70.5 132,2 180,3 203,1 323,7

А

а)

— расчет ,

— эксперимент

1-и то/г

2-й то/г

3-й тон

5-й тон

ростей вращения от нуля до 5000 оборотов в минуту. Вид лопасти в плане показан на рис. 5,а, а формы собственных колебаний приведены на рис. 5, б (четвертый тон соответствует колебаниям в плоскости наибольшей жесткости и на рис. 5,6 не -показан). В табл. 2 приведены частоты собственных колебаний лопасти в исследуемом диапа-

Таблица 2

Источник данных Число оборотов лопасти в минуту Частоты собственных колебаний, Гц

А А /з Л Л

Эксперимент 0 45,0 200.0 400,0 500,0 800,0

Расчет с помощью согласованных треугольных элементов 0 1000 3000 5000 44.8 64.8 91,6 111,8 207,8 216,5 231,1 243,3 470,1 490,0 527.5 562.5 617,9 619,6 623.0 626.0 799.3 812.3 837,5 859,7

Расчет с помощью 0 41,4 197,1 437,7 617,6 745,7

несогласованных треугольных элементов 1000 61.1 204,1 456,7 619,5 759,2

3000 87,2 215,1 492,8 622,5 783,6

5000 106,2 224,3 526,5 625,6 804,9

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

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

ЛИТЕРАТУРА

1. Вольмир А. С. Устойчивость деформируемых систем. — М.: Наука, 1967.

2. Зенкевич О. Метод конечных элементов в технике.—М.: Мир,

1976.

3. Белоус А. А. Устойчивость овальных и рамных шпангоутов. —

Труды ЦАГИ, 1937, вып. 334.

4. The NASTRAN theoretical manual.—NASA, Washington, 1972.

5. А г а п о в В. П. Итерационный метод расчета упругих систем на устойчивость. — Ученые записки ЦАГИ, 1983, № 6, т. XIV.

6. Агапов В. П. Учет геометрической нелинейности в статических и динамических расчетах конструкций методом конечных элементов. Печатается в настоящем выпуске Ученых записок ЦАГИ.

7. Н овожилов В. В. Основы нелинейной теории упругости. —

М.: ГИТТЛ., 1948.

8. Б е з у х о в Н. И. Основы теории упругости, пластичности и ползучести.—М.: Высшая школа, 1961.

9. Г р и г о л ю к Э. И., Кабанов В. В. Устойчивость оболочек. —

М.: Наука, 1978.

Рукопись поступила ЗОЦХ 1982 г.