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

Д. В. Бережной, М. К. Сагдатуллин, А. А. Саченков

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

Ключевые слова: ортотропные и многослойные оболочки, переходные конечные элементы.

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

Keywords: orthotropic and multilayer shells, transient finite elements.

The technique realizing calculation combined shell of a construction, consisting of various elements of multilayered shells and three-dimensional structure, on stability is constructed. On the basis of a three-dimensional finite element of the continuum the family of the new finite elements is constructed, allowing joining elements of various structures in uniform element model. Calculation of the stabilizer of the plane on durability and stability for various cases of load is lead, comparison of results of calculation with the data received in package ANSYS is lead.

Введение

При построении физических моделей

многослойных ортотропных оболочек, вводятся

всевозможные подходы, основанные как на различных гипотезах для каждого слоя оболочки [1-3], так и на единых гипотезах для всех слоев тонкостенной конструкции [2,4-9]. В первом случае порядок разрешающей системы зависит от количества слоев. Во втором случае порядок системы не зависит от числа слоев, что открывает, в частности возможности, для эффективного применения МКЭ в расчетах

ортотропных многослойных оболочек. Данный ортотропный многослойный конечный элемент (КЭ) получен на основе трехмерного изотропного восьмиузлового КЭ оболочки [10], состоящего в дискретизации трехмерных уравнений теории упругости в криволинейной системе координат и использовании некоторых оболочечных гипотез. Применение подобных подходов к расчету оболочек средней толщины посвящены работы [11-18], введение ортотропности и многослойности в расчетную схему оболочек описано в публикациях [17-21].

Для расчета комбинированных конструкций и конструкций существенно переменной толщины методом конечных элементов целесообразно использовать как трехмерные конечные элементы (КЭ) [22], так и оболочечные. Для их стыковки в единую расчетную модель используются различные способы. В [23] описывается создание и исследование переходных трехмерных КЭ, используемых при решении термоупругих задач, в [24] используются переходные конечные элементы, связывающие в единую расчетную модель квадратичные трехмерные элементы сплошной среды и элементы оболочек типа Тимошенко, в качестве узловых степеней свободы которых выбираются три проекции вектора перемещений и два угла поворота нормали, в [25] для расчета напряженно-деформированного состояния пространственных конструкций применяется комбинированная модель, где кинематические условия упругого сопряжения с оболочечными элементами

учитываются при помощи метода штрафа. В [26] необходимые условия сопряжения по границе

между трехмерными и оболочечными КЭ реализуются путем введения в исходный функционал задачи множителей Лагранжа,

параметры которых исключаются из числа

варьируемых величин на уровне сборки конструкции. Обзор некоторых алгоритмов создания переходных КЭ приведен в [27].

В настоящей работе вводятся

аппроксимации радиус-вектора, ковариантных и контравариантных базисных векторов, метрических тензоров, перемещений, тензоров деформаций и напряжений. Технология использования метода двойной аппроксимации по точкам суперсходимости и методика «понижения порядка аппроксимаций» деформаций поперечного сдвига в трехмерной постановке подробно описаны в [10]. Гипотеза малости напряжений обжатия описана в приближенном виде, «упрощенный закон Гука», связывающий напряжение обжатия с деформацией обжатия, не используется. Определяется матрица упругих констант для ортотропного материала. Упругие характеристики и угол намотки могут быть в общем случае различны для каждого слоя конечного элемента. Описаны соотношения сборки матрицы жесткости для многослойного

ортотропного КЭ.

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

«упрощенный закон Гука», сокращенное

интегрирование и т.д.) на квадратурные точки, которые затем интерполируются в узлы в месте стыковки элемента с оболочечным конечным элементом. Узлы, в которых происходит стыковка переходного конечного элемента с массивными трехмерными элементами, освобождаются от наложения оболочечных гипотез. Таким образом,

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

Критерием оценки напряжений в данной зоне является исследование значений напряжений при разбиении данной зоны различными типами конечных элементов (трехмерными, оболочечными, переходными).

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

1. Кинематические соотношения

В пределах каждого элемента вводится локальная система координат #1 , #2 #3, переводящая

искривленный параллелепипед в единичный куб. Определим вектор-функцию произвольной

материальной точки тонкостенной конструкции.

Г = Х (#1,#2,#3 ),

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

X'- (,#2,#3 ) = Х х‘,Ы, (,£,£ ),

/ =1

где Ы, (, #2, #3) - функции формы КЭ оболочки. Определим ковариантные базисные векторы

дх' '

Г* = дтте' = Г е.

д#

Компоненты соответствующего базиса запишем в виде

г' = у х'

* ^ ' д#*,

а ковариантный базис вычислим по формулам

* д#* *'

Г * =-2- = Г* •' е,..

дх' '

Метрические тензоры в различных базисах,

характеризующие геометрию можно представить в

виде следующих соотношений

о. = ГГ. = V ГтГт,

&У ' . ^ ’

т

о"- = Г 'г " = V г',тг.’т.

т

Введем аппроксимации вектора и и градиента вектора и. перемещений и = ит (,£,£ ),

^ т дNs

и' = ет V и. ~д#Г.

£ =1 д^

Компоненты тензора деформаций определим в следующем виде

е.= 2 (г + ги ) = Е V <Е,

где

-г ,т

"V ,

1.

дЫ дЫ дЫ дЫ

+ ■

д? д?

Технология понижения порядка аппроксимации для деформаций поперечного сдвига єі3, где і = 1,2 подробно описана в [10].

Обобщенный закон Гука запишем в следующем виде

с.. = О.. є . (1)

V і.тп тп V '

Введем гипотезу малости напряжений обжатия

Сзз = 0. (2)

Используя соотношение (1), запишем выражение (2) в следующем виде

(3)

(4)

(5)

О33 є = 0.

33тп тп

Представим (3) в виде выражения /0 р + /0 р + /0 р + /0 р + /0 р +

^3311°11 ~-^3312°12 ~-^3313°13 ~-^3321°21 ~-^3322°22 ~

+Н3323р23 + Н3331р31 + Н 3332р32 + Н3333р33 = 0

Выразим из (4) деформацию обжатия

р33 = _(Н3311р11 + Н3312р12 + Н3313р13 + Н3321р21 +

+Н3322р22 + Н3323р23 ++ Н3311р31 + Н3332р32 ) 1 Н3333

Подставляя (5) в (1) получим

= В'"11р11 + В"12р12 + Н"13р13 + Н"21р21 + Н"22р22

+Н"23 р23 + Н"31р31 + Н"32р32 ~ Н"33 3311 р11 +

+Н3312р12 + Н3313р13 + Н3321р21 + Н3322р22 + Н3323р23 + +Н3311р31 + Н3332р32 ) 1 Н3333

Запишем (6) в более простом виде

(6)

с.. = О. є =

V і/тп тп

Оутп

О/33 О33тп ^

О3333

где

О/11 = О/11 -

О..33 О

О/13 = О/13 —

у33^3311

О

і.33 3313

О/12 = 0і/12 —

о..33 О.

3333

о./33 О,

0і/22 0/22

03333 0/3303322

О 21 = О/ 21 —

1.33 3312

О

1.33 3321

3333

О./33 О,

О/31 = О/31 -■

О3333 О/33 О3331

О V 23 = О 23 -

О3333 О..33 О.

О3

О/32 = Оі/32 -

і.33 3323

О3333 Оі/32 О3332

О3

О3333 О3333.

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

2. Ортотропность и многослойность

Среди множества способов расчета слоистых композитных оболочек, особое место занимают численные методы, в частности МКЭ. Главное достоинство подобных алгоритмов состоит в возможности расчета тонкостенных конструкций сложной геометрии с переменными механическими характеристиками под действиями различных нагрузок.

Для анизотропного материала в общем случае соотношения теории упругости имеют

довольно сложную структуру и записываются в виде компонент матрицы упругих констант [ О]

размерностью 6 х 6.

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

относительно которых определяются деформации и напряжения. Обобщенный закон Гука для

ортотропного тела имеет следующий вид

є11 = (с11 — ^12с22 — А13С33 ), У23 = 2є23 = а ,

Е1 а23

1 / \ С12

Є22 = Т^(,С22 — А23С33 — М21С11 ),У12 = 2Є12 = ~р, ,

Е.

^ (33 А32С11 А32С22 ), у31 2є31

а

а,

^3 '-'31

Здесь упругие постоянные выражены, через применяемые в технике показатели: модули упругости Е и сдвига О, а также коэффициента поперечной деформации ц. Таким образом, для установления связей между тензорами напряжений и деформаций необходимо 12 упругих постоянных, из которых 9 постоянно независимы. При этом коэффициенты Пуассона и модули упругости связаны соотношениями

М12 = М21 М23 = М32 М31 = М13

Е1 Ег’ Е2 Е3 ’ Е3 Е1 '

Тогда компоненты матрицы упругости для ортотропного материала запишутся в следующем виде

О1111 =

1

ае2

1

V Е3

Ал

Е2 У

1

ае3

_1_

Е

О1122 О2211

О2233 О3322

ае3

1

О1133 О3311

ае1

1

ае2

^31^23

е2 е1 У

И12И31 , М23

Е3 Е2 )

М12М23 .А '1

Е Е3 У

О2323 а23,

' О1212 а12 ,

О3131

(7)

О3333

ае1

V Е2

Ал Е

где

а =

еЕе (1 — 2^12^23^31 — Е- А321 — Е А122 — Е3 А223). (8)

^3 ^1 2

Наиболее часто в практике расчетов многослойных конструкций встречается преобразование

коэффициентов матрицы упругости ортотропного тела при повороте системы координат вокруг своей оси Ое3

(которая совпадает с нормалью к плоскости слоя) на угол р. Матрица коэффициентов упругости [Б] в осях

слоя Ое1 е'2е3 имеет вид (7), (8). Преобразование коэффициентов упругости осуществляется следующим образом

[ Б**] = [Д]Т [ Б*][Д], (9)

где для случая поворота вокруг оси Ое3 коэффициенты матрицы угла поворота [@] имеют следующий вид:

[Д] =

(10)

с2 s 0 5с 0 0

52 с2 0 - 5с 0 0

0 0 1 0 0 0

-2.se 25с 0 с2 - 52 0

0 0 0 0

0 0 0 0

где с = Сор, 5 = Бтр.

Подставив (10) в (9), получим следующие соотношения для ненулевых компонентов симметричной матрицы [Б]:

Б1112 =\_ Б1111с +(Б1122 + 2Б1212 )(^ — с )— Б22225 ]Х Х5с, Б1111 = Б1111с + Б22225 + 2 (Б1122 + 2Б1212 ) 5 С , Б2212 = |^Б11115 +(Б1122 + 2Б1212 ) (с — 5 )— )2222с ]Х Х5с, Б2222 = Бцц5 + Б^222с + 2^Бц22 + 2Б1212) 5 с , Б1122 = Б1122 (с + 5 ) + (Б1111 + Б2222 — 4Б1212 ) 5 с ,

Б1212 =( Б1111 — 2Б1122 + Б2222 ) 5 с + Б1212 (с — 5 ) , Б1133 = Б1133 с + Б2233 5 , Б3312 = (Б1133 — Б2233 ) 5с,

Б2233 = Б11335 + Б2233с , Б2323 = Б2323с + Б31315 ,

Б2331 = (—(2323 + Б3131 )5с, Б3131 = Б23235 + Б3131с ,

Б3333 = Б3333 .

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

А N N

А* = Т, А = ЁА*, Ёд* =1

А *=1 *=1

где N — количество слоев. Если толщина каждого слоя одинакова, то А * = 1/N. Межслойные координаты по толщине внутри КЭ можно определить следующим образом £ = —1, * = —1 + 2А1, * = —1 + 2А1 + 2А2,

.., N =—1 + 2§ А *,..., £ N+1 = 1.

*=1

При вычислении матрицы жесткости применяется квадратурная формула Гаусса-Лежандра порядка 2 х 2 х N, где N — количество слоев. Таким образом, координаты квадратурных точек в направлении оси N определяются в середине толщины каждого слоя

= 2 (N* +£+1).

Тогда компоненты матрицы жесткости можно вычислить следующим образом

БгГь = 2А*,

где г, с = 1,8, а, Ь,',", т, п, о, р = 1,3.

Компоненты напряжений вычисляются в четырех квадратурных точках, определенных в

середине каждого слоя и далее интерполируются в узлы. Так как при расчете НДС тонкостенных конструкций используется однослойная

аппроксимация, то приемлемая точность достигается уже при количестве слоев N > 3 .

3. Линеаризованная задача устойчивости

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

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

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

Исходными уравнениями для этой задачи являются уравнения виртуальных работ для основного и возможного смежного равновесного состояний Ц] = Л!¥5иЩ + ЦРЯиЖ, , (11)

V V sa

в которых для деформаций используются геометрически нелинейные соотношения. Например, исходное деформированное состояние характеризуется тензором

1

* 2

Л

дг ди + дг дц + ди ди

дх. дх. дх. дх. дх. дх.

V 1 і і1 . і у

а смежное - тензором

Г Зи ди* ^

дх . дх.

1 і У

* ^ ди ди* ^

------1----

.дх . дх .

V . 1У

дг

дх,.

ди ди

--------1--------

дх, дх.

*\

Г ди ди

+ 1--------1-----

V дхі дхі

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

(

дг ди

----------1---------

дх дх,.

д(8и)

дх,.

(

дг ди

-------1------

дх. дх.

V 1 1У

д (8ц)

дх,

йУ =

(12)

= [[[ Р 8 ийУ0 + Л Р8ийБс.

В смежном равновесном состоянии варьируется

у_, *

приращение перемещений и , и уравнение виртуальных работ имеет вид

[

с +с

\1( дг ди ди

Л1 дх.. дх.. дх..

д(8и )

дх.

+!_* +^и +зці і38:) } йу=

V дх. дх. дх! у дхі }

= [[[ Р8ийУ0 + [[ Р8цй»с,

(13)

где с. - приращение напряжений, вызванных

*

перемещением и .

Вычитая из уравнения (13) уравнение (12) с учетом равенства (в силу их произвольности) вариаций 8ц и 8ц* , получим уравнение нейтрального равновесия

(

дг ди Зи

--------1---------1-------

дх. дх. дх.

Л

д (8и*)

3х.

дг ди ди* ^

--------1-------1----

дх,. дх,. дх,.

Vі 1 1 У

д(8и*)

дх..

(14)

+с,,

Гди* 3(8*) + ди* 3(8*)^

дх.. дх,.

дх,. дх..

■йУ = 0.

] 1

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

и условием

• <<-

дг

ди ^ ди ди

сХ( дх( дх( дх(

что означает близость смежного равновесного состояния - так называемая задача устойчивости в малом. Второе предположение огранивает

перемещение и условием ди дг

дх дх( ’

что эквивалентно линейности докритического состояния. С учетом этих неравенств уравнения (12), (14) сводятся к линейным уравнениям (11) и

Г *8 * ди* Зи* V г,

с,,8є,, + с,,----------йУ = 0 .

11 11 11 дх. дх,.

V 11 У

(15)

+

*

<Т

где

Л

дг ди + дг ди

дх. дх. дх. дх.

V ' 1 1 'У

(16)

р=-2 . 2

л

дг ди + дг ди дх. дх. дх. дх.

V ' 1 1 'У

(17)

х дх дх 1 1

Использование конечно-элементной техники для уравнения (11) приводит к известной матричной задаче

[К ]Ы = Н, (18)

где [К] - матрица жесткости, {д} - вектор узловых перемещений, {Р} - вектор внешних сил. Из этого

уравнения определяется первоначальное НДС. Уравнение (15) приводит к однородной системе уравнений

([ К ] + [ Ка ])*} = 0, (19)

где {д*} - вектор приращений перемещений, [Ко] -

матрица геометрической жесткости.

Появление в уравнении (19) той же самой матрицы жесткости [К], что и в (18), обусловлено двумя,

факторами: во-первых, физические соотношения в докритическом состоянии и при выпучивании одинаковы, что можно записать в виде

А1 = |о'118р11ё¥ = |{р}' [Б]{дг}Ж,

V V

А2 = ! = ! {р*}' [Б] {&*} ,

V V

во-вторых, вследствие (16), (17) отличие между деформациями {р} и {р} лишь в физическом смысле

*

векторов и и и и при одинаковой аппроксимации для

них, т.е.

{р = [ир{д}, р} = [ир{д*},

получим

А = {и}' ![ирТ [Б][и^{дд} = {д}' [К]{дд} ,

V

А2 = {д*}' № [Б][ив]^{<5д*} = {д*}" [К]{д*} .

V

Матрица геометрической жесткости [Кст] является

функцией решения задачи (18), т.к. в нее входят

значения {ст.}, характеризующие докритическое

напряженное состояние. В зависимости от величины этих напряжений возможны такие ситуации, когда матрица системы уравнений (19) вырождается и эта система допускает ненулевое решение. Физически это означает, что возможно появление нового равновесного состояния, отличного от найденного из

(18), с полем перемещений и на величину и*.

Математически условие вырождения системы

(19) имеет вид

Л* ([ К ] + [ Кт]) = 0.

Используя линейность всех зависимостей, считаем, что вектор нагрузки пропорционален некоторому параметру Я, т.е.

{РМ{Р,}, (20)

где {Р0} - вектор единичной нагрузки,

характеризующий ее распределение по области;

- из решения задачи (18) определяем т. = ХоТ, что

позволяет записать матрицу геометрической жесткости в виде

[Ко] = 4КО] ; (21)

- подставляя (21) в (19), получаем обобщенную задачу на собственные значения

[ К ]{д* } + Х[ КО ]{д* } = 0; (22)

- решая эту задачу по приведенной далее схеме, находим минимальное собственное значение ХШ1П (или спектр собственных значений) и

соответствующий ему собственный вектор {дШ1п}.

Наименьшая критическая нагрузка, в соответствии с

(20), будет определяться как

{Рсп1 } = Хш1п {р0 } ,

а вектор {дШ1п} характеризует форму потери

устойчивости с точностью до постоянного множителя. Минимальное собственное число и соответствующий ему собственный вектор задачи (22) определят наименьшую критическую нагрузку и соответствующую ей форму потери устойчивости (до постоянного множителя).

4. Метод итерации подпространств

Задача определения всего спектра собственных частот задачи (22) является излишней и вполне достаточно ограничиться определением лишь части спектра, которая соответствует минимальным Х:г" . Для решения этой частичной проблемы собственных значений разработано семейство методов, известных в литературе как методы итераций подпространств. Используется одна из простейших схем такого метода.

1. Имеем некоторую систему линейно независимых векторов

{^0°}, ' = 1, •••,т , (23)

размерность которых совпадает с размерностью задачи (22) N, и число которых значительно меньше этой размерности (т << N).

2. Будем искать решение задачи (22) в виде линейной комбинации (23), т.е.

т

{д} = ЕА{°} .

(24)

3. Составляется отношение Релея, определяющее квадрат частоты колебаний, наиболее близкое по форме к приближению (24)

мд) = —

М' [к]{д} _ Ъ?'^1

(25)

>, 1=1

где

к-К0 }' }, (Ко% =К)}' КН*.}.

4. Параметры разложения (24) Д. находятся из условия минимальности отношения Релея (25), т.е.

'=1

дХ

дД

( ,

из условия минимальности по величине определяемой частоты собственных колебаний

Г( т Л( т \

ък д, ъдк).д

IV1=1 У V1

¥ т Л1

Ъ(к'.), д. Ъдкд Г

V.=1________У V д.=1____У_ = 0

( т Л2 = .

!Д(КЛ Д

V д1=1

Из условия равенства нулю числителя с учетом обозначения (25) получаем алгебраическую задачу

т т

ъ кд,+^Ъ(ко); д=^ (26)

1=1 1=1

которая по своей природе тоже является обобщенной задачей на собственные значения. Из ее решения необходимо найти спектр собственных значений Х(п), п = 1, • , т и соответствующих собственных векторов

{Дп>}'={д(п), Д п}, •, дд п)}.

Число параметров Д существенно меньше неизвестных в системе уравнений (22).

5. Приведенная собственная проблема (26) решается методом вращений Якоби (т.к. матрицы [К] и [Ко]

являются симметричными, а матрица [ К ] к тому же

является положительно определенной.) и находятся

Х(п) и {Д(п)}. Найденные значения Х(п),

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

соответствующие формы потери устойчивости определяются как разложения

т

к }=ъд п) {с}.

(27)

Точность решения (27) целиком и полностью зависит от качества выбора векторов {*(')}, т.е. от

базиса в N -мерном пространстве векторов, в котором ищется решение задачи (22). Как правило, весьма затруднительно сразу определить этот базис (совокупность векторов {*('}}) наилучшим образом.

Поэтому предусматривается процедура изменения этого базиса.

6. Определяется новая совокупность векторов {*}, как решения системы уравнений

[ К ]{) } = —Х(') [ Ко]{д>}.

Далее, повторяются все действия, начиная с пункта 2. Критерием остановки процесса является проверка вычисления величин Х(п) на уточнение их значений на двух соседних итерациях.

5. Расчет стабилизатора и руля высоты легкомоторного самолета на прочность и устойчивость

В качестве примера приведен расчет напряженно-деформированного состояния руля

высоты и стабилизатора легкомоторного самолета.

Стабилизатор и руль высоты представляют собой сложную комбинированную оболочечную конструкцию (рис. 1), усиленную корневыми

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

Рис. 1 - Стабилизатор и руль высоты (вид сверху)

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

0.511

Л

Схема нагружения 1

ПП_ЕГТ~

Схема нагружения 2

Т-------1—------

111

Схема нагружения .1

Т~Г"п—г-—

- Т| и‘1'

Рис. 2 - Схемы нагружения

Схема крепления стабилизатора к фюзеляжу приведена на рисунке 3 (руль высоты крепится к стабилизатору). Исходя из симметрии конструкции, приведена только половина стабилизатора. Крепление к фюзеляжу (для половины стабилизатора) производится в трех точках А , В и С (рис.3). В точке А реализован скользящий шарнир, в точке В - обыкновенный шарнир, в точке С стабилизатор крепится к килю с помощью расчалки.

и2=о г

А -----1

иу—и2 с

их=о

ть=о

и2=о

Рис. 3 - Схема кинематического закрепления стабилизатора

Для верификации программы было проведено сравнение решения задачи с

'=1

результатами, полученными в ППП А№У8 на основе конечного элемента 8ИБЬЬ181. В таблице 1 приведено сравнение реакций Ях, Яу и Я2, возникающих в точках крепления стабилизатора к фюзеляжу в точках А и В .

Таблица 1 - Реакции (Н) в точках крепления А и В , передаваемые на фюзеляж со стабилизатора для трех случаев нагружения

№ Ях

рі^г ansys progr ansys progr ansys

1 А - - -5841 -5062 63,2 -9,2

В -524 -443 2164 1926 26,3 0

2 А - - -9124 -8615 -21 -16

В -349 -337 11016 10497 -14 -3

3 А - - -7811 -7223 -18 -13,8

В -414 -399 5642 5236 -37 -2

Анализ результатов расчета показывает, что разница в результатах расчета силовых факторов при статическом расчете сложной конструкции между полученными по разработанной методике и в ППП А№У8 не превышает 15%.

Для иллюстрации результатов расчета на рисунке 4 приведено распределение прогибов для половины конструкции стабилизатора Ж (в метрах) для схемы нагружения 2.

•.«$5514 -.Й«55 -.#Э#Э*7 -.#171» -.М5.'7*

-.ММ» -.«24117 -.•1155% . 1##Е*#«

Рис. 4 - Общая картина прогибов половины конструкции стабилизатора для схемы нагружения 2

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

Таблица 2 - Отношения первых четырех значений величин критических нагрузок к статическим.

Номер критической нагрузки 1 2 3 4

Схема нагружения 1 1,3931 1,4027 1,6569 1,6576

Схема нагружения 2 1,0251 1,0334 1,0991 1,1079

Схема нагружения 3 1,2361 1,2489 1,3136 1,3263

Можно отметить, что наиболее опасным с точки зрения расчета на устойчивость является вариант нагружения по схеме 2. В этом случае все четыре минимальных критических нагрузки превышают расчетную не более, чем на 10%.

Заключение

Исходя из полученных результатов, следует отметить, что разработанная численная методика исследования напряженно-деформированного

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

Для результатов, полученных в ППП А^У8 использовалась лицензионная версия (l1cense-ansys_9-13-2012-6.57.22), установленная в Казанском (Приволжском) федеральном университете.

Работа выполнена при частичной финансовой поддержке РФФИ в рамках научных проектов № 12-01-00955, № 12-01-97026.

Литература

1. Болотин, В.В. Механика многослойных конструкций /

В.В. Болотин, Ю.Н. Новичков. - М.: Машиностроение, 1980. - 375с.

2. Григолюк, Э.И. Многослойные армированные оболочки. Расчет пневматических шин / Э.И. Григолюк, Г.М. Куликов. - М.: Машиностроение, 1988. - 288с.

3. Соловьев, С. С. Конечноэлементная модель многослойной оболочки с анизотропными слоями переменной толщины. / С.С. Соловьев // Известия Вузов. Авиационная техника. - 1989. - №4. - С. 71-75.

4. Алфутов, Н.А. Расчет многослойных пластин и

оболочек из композиционных материалов / Н.А. Алфутов, П. А. Зиновьев, Б.Г. Попов. - М.:

Машиностроение, 1984. - 264с.

5. Васильев, В.В. Механика конструкций из композиционных материалов / В.В. Васильев. - М.: Машиностроение, 1988. - 270с.

6. Григоренко, Я.М. Статика анизотропных оболочек с конечной сдвиговой жесткостью / Я.М. Григоренко, А.Т. Василенко, Г.П. Голуб. - Киев: Наукова Думка, 1987. -216с.

7. Пискунов, В.Г. Расчет неоднородных оболочек и пластин методом конечных элементов / В.Г. Пискунов. -Киев: Вища Школа, 1987. - 200с.

8. Рикардс, Р.Б. Метод конечных элементов в теории оболочек и пластин / Р.Б. Рикардс. - Рига: Знание, 1988.

- 284с.

9. Рассказов, А.О. Теория и расчет слоистых ортотропных пластин и оболочек / А.О. Рассказов, И.И. Соколовская, Н.А. Шульга. - Киев: Вища школа, 1986. - 191с.

10. Голованов, А.И. Трехмерный конечный элемент для расчета тонкостенных конструкций / А.И. Голованов, М.К. Сагдатуллин // Ученые записки Казанского государственного университета. Серия физ. -мат. наук.

- Казань 2009. - Т. 151, кн. 3. - С. 121-129.

11. Перелыгин, О. А. Исследование прочности цилиндрических оболочек с вмятинами в области радиальных соединений / О.А. Перелыгин, Ш.Ш. Галявиев, Р.Х. Зайнуллин, Д.В. Бережной // Вестник КГТУ, Казань, изд-во КГТУ. - 2001. - №1-2. - С.75-77.

12. Перелыгин, О. А. Исследование прочности цилиндрических оболочек при наличии увода или смещения кромок сварных швов / О.А. Перелыгин, Н.М. Туйкин, Д.В. Бережной, М.Н. Серазутдинов // Вестник КГТУ, Казань, изд-во КГТУ. - 2001. - №1-2. - С.77-79.

13. Бережной, Д.В. Многослойный ортотропный конечный элемент оболочек средней толщины / Д.В. Бережной, М.К. Сагдатуллин, А.И. Голованов // Вестник Саратовского государственного технического ун-та. - 2011. -№3 (57). Выпуск 1. - С. 9-19.

14. Бурман, Я.З. Расчет упругопластического

деформирования оболочек на основе теории течения и МКЭ. Исследования по теории пластин и оболочек / Я.З. Бурман, С.С. Соловьев. - Казань: Изд-во КГУ. - 1990. -№22. - С. 98-107.

15. Parish, H.A. Critical survey of the J-node degenerated shell element with special emphasis on thin shell application and reduced integration / H.A. Parish // Computer Methods in Applied Mechanics and Engineering. - 1979. - vol. 20. - №3. -p. 323-350.

16. Zienkiewicz, O.C. Reduced integration technique in general analysis of plates and shells / O.C. Zienkiewicz, R.L. Taylor, J.M. Too // Int. J. Num. Meth. Eng. - 1971. - v.3. - p. 275-290.

17. Голованов, А.И. Современные конечно-элементные модели и методы исследования тонкостенных конструкций / А.И. Голованов, А.В. Песошин, О.Н. Тюленева - Казань: КГУ, 2005. - 442с.

18. Голованов, А.И. Метод конечных элементов в статике и динамике тонкостенных конструкций / А.И. Голованов, О.Н. Тюленева, А.Ф. Шигабутдинов. - М.: ФИЗМАТЛИТ, 2006. - 392с.

19. Голованов, А.И. Метод конечных элементов в механике деформируемых твердых тел / А.И. Голованов, Д.В. Бережной. - Казань: Издательство «ДАС», 2001. -301с.

20. Попов, Б.Г. Расчет многослойных конструкций вариационно-матричными методами: Учебное пособие. / Б.Г. Попов- М.: Изд-во МГТУ, 1993. - 294с.

21. Хечумов, Р.А. Применение метода конечных элементов

к расчету конструкций: Учебное пособие для

технических вузов / Р.А. Хечумов, Х. Кеплер, В.И. Прокопьев. - М.: Изд-во Ассоциации строительных вузов, 1994. - 353с.

22. Бережной, Д.В. Расчет комбинированных конструкций

методом конечных элементов / Д. В. Бережной, М. К. Сагдатуллин, А.А. Саченков // Научно-технический вестник Поволжья. - Казань: Научно-технический

вестник Поволжья, 2012. - № 4.- С.13-16.

23. Surana, Karan S. Three dimensional solid-shell transition finite elements for heat conduction / Karan S. Surana // Comput. and Struct. - 1987. - V. 26. - № 6. - P. 941-950.

24. Бережной, Д.В. Конечные элементы для расчета

конструкций существенно переменной толщины / Д.В. Бережной, И.Ю. Красновский. // Статика и динамика элементов конструкций сложной формы: межвуз.

сборник. - Наб. Челны, тез. докл. респ. науч.-техн. конф. КамАЗ-КамПИ.- Наб.Челны : КамПИ,1990. - С.30-36.

25. Савула, Я.Г. Применение комбинированной модели

для расчета напряженно-деформированного состояния пространственных конструкций / Я.Г. Савула,

И.И. Дыяк, А.В. Дубовик // Прикладная механика. -1989. - Т. 25. - № 9. - С. 62-67.

26. Адясова, Н.М. Некоторые вопросы расчета нелинейных составных конструкций / Н.М. Адясова, С.А. Капустин, Л.С. Яблонко // Прикладные проблемы прочности и пластичности. - Горький, 1975. - Вып. I -

С. 124-135.

27. Surana, Karan S. Transition finite elements for threedimensional stress analysis / Karan S. Surana // Int. J. Numer. Meth. Eng. - 1980. - V. 15. - № 7. - P. 991-1020.

© Д. В. Бережной - канд. физ.-мат. наук, доц. каф. теоретической механики КФУ, Dmitri.Berezhnoi@mail.ru; М. К. Сагдатуллин - канд. физ.-мат. наук, ст. препод. каф. теоретической механики и сопротивления материалов КНИТУ, ssmarat@mail.ru; А. А. Саченков - канд. физ.-мат. наук, доц. каф. теоретической механики КФУ.