CC BY
8
Illl
УДК 691.33
А.А. Абдюшев - старший научный сотрудник лаборатории нелинейной механики оболочек Института механики и машиностроения Казанского научного центра Российской академии наук (ЛНМО ИММ КНЦ РАН), старший преподаватель кафедры мостов и транспортных тоннелей Казанского государственного архитектурно-строительного университета (КазГАСУ)
МОДЕЛИРОВАНИЕ ПОДКРЕПЛЕННОЙ ОБОЛОЧКИ В АРС ЭРА-ПК2000
АННОТАЦИЯ
Применение метода конечных элементов к расчету подкрепленных оболочек при анализе общей прочности летательных аппаратов предъявляет определенные требования к моделированию поведения отдельных её элементов и корректного учета их взаимодействия. Эта проблема, в определенной степени, была решена в Автоматизированной расчетной системе (АРС) ЭРА-ПК2000. Работа посвящена математическому выводу матриц жесткости конечных элементов, моделирующих изгибаемые подкрепления и обшивку оболочки, реализованных в АРС. Предложен простой алгоритм учета влияния обшивки на напряженное состояние подкрепляющих элементов.
A.A. Abdiushev - senior scientific researcher of the Nonlinear Mechanical Engineers of Shells of Institute Mechanical Engineers and Machine Building Laboratory of Kazan Scientific Centre of Russia Academy of Sciences (LNMS IMM KSC RAS), senior lecturer of the Department of Bridge and Transport Subways of Kazan State University of Architecture and Engineering (KSUAE)
MODELING OF SUPPORTED SHELL IN ESA-PC2000
ABSTRACT
The application of a method of finite elements to the calculation of the supported shells in the analysis of general durability of flying devices determines the certain requirements to modeling behaviour of its separate elements and the correct account of their interaction. This problem, in the certain degree, has been solved in the Automatic Settlement System (ASS) ESA-PC2000. The article is devoted to a mathematical conclusion of matrixes of rigidity of final elements. These matrixes are simulating the bent reinforcements and the shell cover, realized in ASS. The simple algorithm of the account of influence of a cover on the intense condition of supporting elements is described.
Расчет подкрепленных оболочек является неотъемлемой частью анализа общей прочности летательных аппаратов. В середине XX века практически одновременно, за рубежом (Эбнер) и в отечественной практике (Беляев), была предложена расчетная модель работы тонкостенной подкрепленной оболочки типа моноккок, суть которой заключалась в “предписывании” разным элементам оболочки определенного характера работы. А именно: поперечные подкрепляющие элементы (шпангоуты) работают на изгиб в своей плоскости (плоскость перпендикулярная оси оболочки), продольные подкрепляющие элементы (стрингеры или эффективные пояса) воспринимают растяжение-сжатие, а обшивка подкрепленной оболочки воспринимает “чистый” сдвиг. С появлением вычислительной техники и развитием численных методов расчета, усилиями Аргириса подобная модель, при расчете фюзеляжа летательного аппарата (ЛА), была реализована Методом конечных элементов (МКЭ) в форме метода сил, на основе стационарности функционала Кастильяно [1]. Реализованная по МКЭ модель позволяла проводить анализ общей прочности элементов фюзеляжа с учетом их взаимного влияния. При этом поперечное сечение фюзеляжа было представлено в форме полигона. Расчеты Аргириса и его последователей [2] показывали хорошие соответствия с экспериментальными данными. Тем не менее, формулировка МКЭ в форме метода перемещений, основанная на минимизации функционала Лагранжа, оказалась, в силу своей алгоритмичности, более привлекательной. Поскольку модель, предложенная Аргирисом, зарекомендовала себя хорошо, возникла потребность “переноса” всех положительных качеств модели на МКЭ в форме метода перемещений. Эту работу можно условно разбить на две части, определяемые подходом к проблеме: формальный - создание конечноэлементной базы, позволяющей реализовать модель и неформальный -создание методики анализа результатов, позволяющей правильно оценить поведение конструкции. Ниже приводятся материалы, позволившие, в определенной степени, решить эти задачи в рамках Автоматизированной расчетной системы
ТЕОРИЯ ИНЖЕНЕРНЫХ СООРУЖЕНИЙ И СТРОИТЕЛЬНЫЕ КОНСТРУКЦИИ
ЭРА-ПК2000, созданной в отраслевой научно-исследовательской лаборатории САПР фюзеляжа вертолета при кафедре строительной механики КазГАСУ (КИСИ) [3].
Для рассмотрения более широкого класса задач при переходе к МКЭ в форме метода перемещений было принято решение рассматривать шесть естественных степеней свободы конструкции в узле при расчете подкрепленных оболочек. Это позволило, в частности, безболезненно включать в конструкцию поперечные и продольные изгибаемые элементы, ориентированные произвольно относительно глобальной системы координат. Поэтому в качестве базового изгибаемого элемента был принят прямолинейный двухузловой конечный элемент (КЭ) пространственного бруса с шестью степенями свободы в узле [5].
Изгибаемые подкрепляющие элементы
Особенностью рассматриваемых оболочек является наличие эксцентриситета положения нейтральной оси изгибаемых элементов и узлов расчетной схемы (рис. 1). Как правило, узлы расчетной схемы обозначают нейтральную ось тонкой безмоментной обшивки. Математическая процедура приведения элементов обшивки к нейтральной оси подкрепляющих изгибаемых элементов не сложнее, чем приведение изгибающих элементов к обшивке. Однако, второй вариант предпочтительнее, поскольку позволяет строить более точные модели. При формировании глобальной матрицы жесткости (МЖ) проблема несовпадения узлов нейтральной оси и узлов расчетной схемы решается статически эквивалентным “переносом” элементной МЖ из одних узлов в другие.
Рис.1. Схема эксцентриситетов
Для малых перемещений имеем:
Р = С ■ Р
ге е г
иГе = Се ■ иг
(1)
где Рг, и г - векторы узловых сил и перемещений в глобальной системе координат, соответствующие точкам на гтральной оси элемента;
Рге ,иге - векторы узловых сил и перемещений в глобальной системе координат, соответствующие узлам расчетной схемы;
С е - матрица “эксцентриситетов” - трансформирующая матрица.
Се =
1
0 1 0
0 0 1
0 Д2І -Ду 1
-Д2І 0 Дх 0 1
Ду - Дх 0 0 0
1
0 1
0 0
0 Д2]
-Д2] 0
ДУз - Дх
1
-Ду,
Дх,
0
1
0 1 0 0 1
(2)
0
0
0
Где Ах = х1 - хґ; Ау = у1 - уґ; Аг = - іґ.
Поскольку
се • Се = Е (3)
тогда
кге • иге = Рге ^ Кге • (Се • иг) = Се • Рг ^ се • Кге • (Се • иг) = се • с
откуда
Ке = Се • К • се (5)
В ЭРА- ПК2000 вычисление эксцентриситетов производится при помощи сервисной программы, исходными данными для которой являются расстояния от наружной полки сечения элемента до нейтральной оси и угол наклона узлового сечения относительно вертикальной оси.
В конструкции фюзеляжа имеются изгибаемые элементы, к которым обшивка примыкает как по “нижней”, так и “верхней” полкам. В этом случае в ЭРА-ПК2000 предлагается применять составной элемент (макроэлемент) балка-стенка. Этот элемент образуется двумя изгибаемыми поясами и стенкой. В качестве поясов выступают рассмотренные выше балочные элементы с эксцентриситетом. Стенка составляется из трехузловых КЭ оболочки, полученных суперпозицией тонкой изгибаемой треугольной плиты по гипотезе Кирхгофа и мембранного трехузлового элемента Бергана-Фелиппа [3]. Элемент Бергана имеет три степени свободы в узле и был отобран в результате тестирования при моделировании продольного изгиба. Геометрические характеристики компонент составного сечения определяются из соотношений, записанных для заданного и эквивалентного сечений (по рис. 2).
Рис. 2. Тестовое сечение изгибаемого элемента и эквивалентного ему сечения макроэлемента балка-стенка
F = F + - (y -y )b; F = F- - (y - y i )b;
в w max sc' ’ н c s mm / ’
Sв = s + - (ymax - Ус )2 b /2; Sh = S - + (Ус - ym,n )2 b /2; s / s /
Ув = s/f .Ун = H/F ,Iye = IyH = Iy/2,IKpe = IKPh = (IKP -bh/3)/2, (6)
вн
Ize = Izh = (Iz - Fe • Ув2 - Fh • yH - b ■ h 3/12 - bh(yc - h/2)2)/2-
(4)
e
Здесь: Е - общая площадь сечения; I - момент инерции сечения относительно оси 7; I у - момент
нейтральной оси части сечения выше этой оси; р и $ - площади и статические моменты инерции относительно нейтральной оси части сечения ниже этой оси; Рн, Ев - площади поясов; $н, $ - статические
моменты поясов эквивалентного сечения относительно нейтральной оси; I , I , Iгв, I гн, Iкрн, 1крв -
геометрические характеристики поясов.
Для уточнения модели изгибаемых подкреплений в АРС ЭРА-ПК2000 предлагается КЭ пространственного бруса переменного сечения. Вывод жесткости элемента переменного по длине сечения также использует допущение о том, что вертикальная ось сечения является главной центральной. Для начального (с индексом
1) и конечного (с индексом 2) сечений элемента составляются соотношения (6). Исходя из эквивалентного представления сечений, можно положить основные характеристики линейно изменяющимися по длине, например:
длиной I (рис. 3) по формуле Мора строится матрица податливости § размерностью 6 X 6 с ненулевыми элементами:
Для получения полной матрицы жесткости элемента переменного сечения формируем матрицу, уравновешивающую узловые силы первого узла относительно 2-го. Эта матрица образуется из матрицы равновесия относительно 2-го узла изменением знака на противоположный у блока, относящегося к 1-му
инерции сечения относительно оси у; Iкр - момент инерции сечения при кручении; Xу , Xг - коэффициенты работы сечения на сдвиг по оси у и 7; Е + и Б + - площади и статические моменты инерции относительно
Рв П) = Рв1 + Р - Рв1 Л , Рн (П) = Рн1 + (Рн2 - Рн1 )п, Н(п) = Н + (¿2- Юп , Ь(п) = Ь + (Ь2- ьі )п
и так далеедля ^в(П), ^н(ПК Ьв(ПX ^н(П), Ув(П), Ун(П^Ьрв(ПX Ьрн(П).
(7)
0
(8)
После обращения матрицы податливости, получаем жесткость узла 2:
(9)
узлу.
Окончательно получаем:
К = В‘ • К22 • В .
(10)
В =
-1 0 0 0 0 0 1 0 0 0 0 0
0 -1 0 0 0 0 0 1 0 0 0 0
0 0 -1 0 0 0 0 0 1 0 0 0
0 0 0 -1 0 0 0 0 0 1 0 0
0 0 -1 0 -1 0 0 0 0 0 1 0
0 1 0 0 0 -1 0 0 0 0 0 1
Рис. 3. Элемент переменного сечения
Элементы, моделирующие обшивку
В качестве базовых элементов в АРС ЭРА- ПК2000 применяются плоские трех- и четырехузловые мембранные, реализующие плоско-напряженное состояние [5], и сдвиговые, реализующие чистый сдвиг в некоторой прямоугольной системе координат, конечные элементы.
Для вывода матрицы жесткости семейства сдвиговых элементов была принята гибридная формулировка на базе обобщенного функционала Кастильяно [4]. Матрица жесткости для произвольного выпуклого п -узлового многоугольника (рис. 4), реализующего чистый сдвиг в плоской декартовой системе координат, выражается следующим образом:
К = Г-И --Т (11)
здесь И - матрица податливости выпуклого п - угольника, Т - трансформирующая матрица.
Рис. 4. Плоский выпуклый п -узловой многоугольник
И -1 = £
2Бп-(1 + V)-1
(11)
Здесь Е - модуль упругости материала, V - коэффициент Пуассона, I - толщина элемента, 25п - удвоенная площадь выпуклого плоского п -угольника.
п—1
25п = -Хп 'У\ - ХГ Уп + X (-X ■ У +1 + X +1 ■ У )
(12)
2
Х2п ; Уп2 ; Х31; у13 ;••• 'Х(г +1)(г-1); У(г-1)(г'+1) ;•• • Х1(п-1); У(п-1)1.
(13)
(14)
где X г, у, - координаты г -ого узла, Х ц - X - .■; Уу - у, - У}- •
Величина сдвига в элементе определяется по вычисленным перемещениям:
т - Н- Т^и Здесь и - вектор перемещений узлов в плоскости
В реальных задачах узлы четырехугольных элементов обшивки часто не являются компланарными Игнорировать этот факт нельзя, так как после получения матрицы жесткости плоского элемента, образованного проекциями узлов модели на некоторую плоскость и формального добавления этой жесткости к жесткости расчетных узлов, нарушаются соотношения равновесия в глобальной матрице. Для того чтобы не нарушать равновесия, в ЭРА-ПК2000 предлагается два варианта трансформирования матрицы жесткости плоского четырехузлового элемента, оба основанные на принципе наименьших квадратов. Укажем один из них, более простой в реализации, тем более, что оба приводят к одинаковому результату
На рис 5 показаны четырехузловые элементы: а) - плоский четырехузловой элемент в пространстве, б) - его прототип Для обоих элементов рассмотрено некоторое одинаковое единичное состояние и показаны реакции этого состояния. Фактически это некоторый столбец (строка) матрицы жесткости Запишем функционал метода наименьших квадратов в виде:
4 4 4
р -X (р - о2 +Х р - о2+Х (р - р;у (15)
г-1 г-1 г-1
С помощью неопределенных множителей Лагранжа добавим в функционал (15) шесть уравнений равновесия для реакций закрученного элемента и получим:
Рис 5^ Исходный а) и закрученный б) четырехузловые элементы в пространстве
р ** - ¿(р - к?+Х(р - р;,)1+Х(р - р' )2+лх р„+к X р. +я X р+
г-1
г-1
г-1
г-1
г-1
г-1
4 4 \ 4 4 \ 4 4 л
14 X р* ■У, -£ РУ, ■* + Я5 X р. ■г . -X р г гг ■Хг + Яб м У Хг -£ рхг ■Уг
и-1 г-1 0 У1-1 г-1 0 ^■-1 г-1 0
Здесь х, , у, , - координаты узлов закрученного элемент
Вариации (16) по р.,, Ру,, РЕ,, я, приводят к системе линейных уравнений в матричном виде:
(16)
-1
$0 • Я = Я'
' Р : 1
: Я' =
1 _ 0
Я
Р = {Р РРРРРРРР
1 V х1 1 у1 1 г\ 1 х 2 1 у 2 1 г2 1 х3 1 у3 ^
у1
у 2
$ о
Р' у2 р' 1 г 2 Р' х3
Е2 0 0 0
0 Е2 0 0
0 0 Е2 0
0 0 0 Е2
Е Е Е Е
$1 $ 2 $ 3 $ 4
у3
у3
о
о
г3
Е $ Е $ 2 Е $3 Е $
0
0
(17)
15 1б):
Р х4 РУ 4
Рх4 Р4
Р 4): Р^):
(18)
Е - единичная матрица 3X3, Е2 = Е • 2;
Уг 0 - г.
- х 0
0 -У Хг
$
Алгоритм “закрутки” будет выглядеть следующим образом:
1) формируется К' - матрица жесткости плоского элемента в глобальной системе координат,
2) вычисляется $01 и из нее выбирается подматрица $ [12 X12]
$ $я
$1 $11
(19)
3) окончательно вычисляется
К = • К'-$ (20)
В продольном направлении в качестве подкрепления используются те же изгибаемые элементы и двухузловые ферменные КЭ.
Анализ взаимодействия подкрепляющих элементов и обшивки
Для более корректного анализа, при получении усилий в двухузловых подкрепляющих элементах, предлагается определять эти усилия по вычисленным перемещениям, используя совокупную жесткость двухузловых элементов и КЭ обшивки. В этом случае предполагается, что дополнительные узловые усилия от элементов обшивки в узлах подкрепляющих элементов будут уравновешены некоторым сдвиговым потоком от панелей обшивки (рис. 6). Такой подход позволяет максимально приблизиться к модели Аргириса [1]. При использовании вместо сдвиговых элементов КЭ мембраны необходимо пересчитывать геометрические характеристики сечений, а для изгибаемых и положение нейтральной оси сечения.
В АРС ЭРА-ПК2000 эта процедура производится автоматически.
Такой подход позволяет получить гладкую и достоверную картину напряженного состояния в любой цепочке подкрепляющих элементов.
Рис. 6. Представление действия обшивки как локальной нагрузки на двухузловой элемент
Литература
1. Аргирис Дж. Современные достижения в методах расчета конструкций с применением матриц. - М.: Изд-во лит-ры по стр-ву, 1968. - 241 с.
2. Бурман З.И., Лукашенко В.И., Тимофеев М.Т. Расчет тонкостенных подкрепленных оболочек методом конечных элементов с применением ЭЦВМ. - Казань: Изд-во Казанского ун-та, 1973. - 569 с.
3. Лукашенко В.И., Абдюшев А.А., Доронин М.М., Нуриева Д.М., Сладков А.В. Экспертиза, расчет, анализ пространственных конструкций. АРС ЭРА-ПК2000. - Казань: Изд-во КГАСУ, 2006. - 320 с.
4. Пиан Т. Вывод соотношений для матриц жесткости элемента, основанный на выборе закона распределения напряжений. // Ракетная техника и космонавтика, 1964, N° 7. - С. 219-222.
5. Przemieniecki J.S. Theory ofMatrix Structural Analysis. - N.Y.: Mc. Graw-Hill Book Company,1968. - 246 p.
