CC BY
48
Бережной Д.В.1, Габибова А.К.2 ©
'К.ф.-м.н., доцент, кафедра теоретической механики; 2аспирант, кафедра теоретической механики, Казанский (Приволжский) федеральный университет
КОНЕЧНЫЙ ЭЛЕМЕНТ ДЛЯ РАСЧЕТА ОДНОРОДНЫХ, ОРТОТРОПНЫХ И МНОГОСЛОЙНЫХ ПЛАСТИН И ОБОЛОЧЕК СРЕДНЕЙ ТОЛЩИНЫ
В данной работе на основе трехмерных уравнений теории упругости реализован один из алгоритмов построения многослойного ортотропного конечного элемента пластин и оболочек средней толщины. Использована гипотеза малости напряжений обжатия, реализован метод двойной аппроксимации, решен ряд тестовых задач.
Ключевые слова: многослойные оболочки, двойная аппроксимация Keywords: multilayered shells, double approximation
При построении физических моделей многослойных ортотропных оболочек вводятся всевозможные подходы, основанные как на различных гипотезах для каждого слоя оболочки [1, 72], так и на единых гипотезах для всех слоев тонкостенной конструкции [2, 214-219]. В первом случае порядок разрешающей системы зависит от количества слоев. Во втором случае порядок системы не зависит от числа слоев, что открывает, в частности возможности, для эффективного применения метода конечных элементов. При расчете пластин и оболочек средней толщины широкое распространение получили специальные элементы пластин и оболочек, которые имеют в качестве узловых степеней свободы перемещения узлов, расположенных на лицевых поверхностях [3, 9; 4, 98].
Деформацию среды описывают компоненты линейных деформаций e xx ,e уу,e zz и деформации сдвига fху,fyz,fzx , которые выражаются через перемещения известными соотношениями Коши. Напряженное состояние представляет тензор напряжений Коши в виде компонент нормальных (s хх,s уу,s 22) и касательных (tху,tyz,t zx) напряжений.
Деформации и напряжения связаны между собой линейной зависимостью, называемой обобщенным законом Гука, которая для анизотропного материала имеет весьма сложную структуру и представляется в виде компонент матрицы упругих констант размером 6‘‘ 6. Для ортотропных материалов вводят систему ортогональных декартовых координат a , b , f , определяющую так называемые оси ортотропии. Относительно этих осей определяются деформации eaa , e bb , £ ГГ , f ab , f bf , f fa и напряжения s aa , S bb , S ff , 1 ab , T bf , % fa , которые связаны между собой обобщенным законом Гука в виде
Аннотация
Если ввести в рассмотрение направляющие косинусы а х, а у, а г, Ь х, Ь у, Ь г, 1 х, 1 у , 1 г между осями а , Ь ,1 и х,у,7, то справедливы зависимости
© Бережной Д.В., Габибова А.К., 2013 г.
2 2 2 «аа = е хха х +е ууа у + « 22а 2 + 1 хуа ха у + 1 у2А Уа 2 + 1 «а А х ,
еЦ = « ххЬ 2Х + е ууЬ 2 +« аР 2 + 1 хуЬ хЬ у + 1 уЛ уЬ 2 + 1 3 Л х,
2 2 2 £„ = е хх1 х +е уу1 у +е 21 2 + 1 ху1 х1 у + 1 уг11 2 + 1 2,1 21 х ,
1а? = 2 хха хЬ х + 2 ууа уЬ у + 2 2а 2Ь 2 + 1 у (а хЬ у +а уЬ х ) +
+ 1 У2 (а уЬ 2 +а 2Ь у )+1 2х (а 2Ь х +а х^ 2 ) ,
Ьг = 2£ ххЬ х1 х + 2 ууЬ у1 у + 2 22Ь 21 2 +1 у (Ь 11 у + Ь у1 х ) +
+ 1 У2 (Ь у1 2 + Ь 21 у )+ 1 2, (Ь 1 2 + Ь 21 х ) ,
1Ьа = 2 х1 а х + 2 уу1 уа у + 2 221 2а 2 +1 у (1 а у + 1 уа х ) +
+ 1 У2 (1 Уа 2 + 1 2а у )+ 1 2, (1 а 2 +1 2а х ) ,
которые можно переписать в матричном виде
{«а} = [ Т ]{«} ,
(1)
Если обозначить приведенные векторы деформаций в осях x, y, 2 и а , Ь ,1
{«} Т = {е хх,е уу,е 22 ,1 ху ,1 У2 ,1 2х} , {«а} Т = {«аа ,е ЬЬ ,£ 11 ,1 а? ,1 Ьа ,11а]
И приведенные векторы напряжений в тех же осях
{^} = {^ хх ,^ уу ,^ 22,7 ху ^ у2^ 2х} , {^а} = {^аа ,^ЬЬ , ^ 11 ^а$ ^ Ьг ^ Га} ,
то соответствующие соотношения упругости можно записать в виде
{*} = №} , К} = [А]{£а} .
В этом случае матрица упругих постоянных [ я] для заданной матрицы [ Д ] определяется следующим образом
[ я] = [ Т ]Т [ Д][ Т ].
Структурно предлагаемый конечный элемент представляет собой искривленный параллелепипед, состоящий из набора N слоев по толщине, каждый из которых является ортотропным материалом с осями ортотропии а, Ь ,1 . При этом предполагается, что плоскость а , Ь параллельна плоскости локальных координат X Л вдоль любой прямой по толщине, т.е. ось ортотропии 1 параллельна локальной координате { (рисунок 1).
Рис.1. - Многослойный конечный элемент
Для задания ориентации осей а, Ь,1 предусмотрены две возможности. Если конечный элемент плоский и его плоскость совпадает с плоскостью х,у, то определяется угол ф , как угол между осью х и направлением а . Если элемент искривленный, то вводится в рассмотрение система ортогональных ортов Д, Д, р3 таким образом, что Д направлен по касательной к координатной линии X , а р3 - перпендикулярен поверхности X Я , т.е. направлен вдоль координатной линии { . Эту систему удобно определять с помощью следующих соотношений
Д = д г/д{/1д г/Ц I Д = (д г/дХ) ( (д г/дц)/1(д г/дХ) < (д г/дц) I Д = Д ' Д.
Ориентация осей а , $ , лежащих в плоскости р1 - р2, задается в виде угла ф между р1 и осью а . Для построения матрицы жесткости необходимо определить матрицу физических соотношений [ & ] в законе Гука для напряжений и деформаций относительно
глобальных осей х,у,2 через матрицу упругих констант [ &а], определяющую закон Гука в осях ортотропии. Для этого следует ввести в рассмотрение матрицу «поворотов деформаций» на угол ф в плоскости а , $ [ Та \, для которой справедливо выражение
[ т,} =
Если конечный элемент плоский и плоскость а , Ь совпадает с плоскостью х,у
С о 2 sin2 § 0 sinf СОБ§ 0 0 щ
К sin2 § cos2 § 0 - sin§ cosf 0 0 Ї
к 0 0 1 0 0 0 ъ
к С о Х/1 - & ъ
к2Бт§ cosf 2sin § cosf 0 0 0 ъ
к0 0 0 0 cosf - sinf ъ
к ъ
л 0 0 0 0 sinf cosf ы
ж
Щ
.4
з § = а , х ч, то и ш
[ &]=[ Та]Т[ В'][ Та] .
Если конечный элемент искривленный, то необходимо ввести в рассмотрение ‘матрицу поворотов” деформаций из системы х,у,2 с ортами /, j, k в систему с ортами
¥, Щ = к
Эта матрица будет иметь вид
й Рі2 к РІ РІу рі РіхРіу Рі уРіг РігРіх Щ
Ріу рі Рі хРі у СЧ Р 2 Р ъ Р 2 гР2 х ъ
н 2 2 ккк рЗу рі Р ЗхР 3 у Рз уР зг РзгРзх ъ
2 Рі уР 2 у 2 Рі.Рі г РіхРі у + Рі хРі у Рі уРі г + Рі у Ріг РігР2х + РігРіх ъ
X СО 2 2 к 2 Р 2 уР 3 у 2 Р 2 гРз г Рі хРз у + Рз хРі у Рі уР з г + Рз уРі г Рі гРзх + РзгРіх ъ
«Р X СО 2 кк 2 Р3 уРі у 2 Рз гРіг РзхРі у + РіхР 3 у Рз у Ріг + Рі уРзг РзгРіх + РігРзх Ы
Заметим, что матрица ЛТР Щ структурно совпадает с матрицей определяющих пересчет деформаций (1). В результате матрица упругих характеристик [ D] будет вычисляться в виде
Т Т
[ Щ = ЛТр Щ [ Т.\ [ Д\[ Т\ ЙТр Щ.
Величины толщин каждого слоя hs удобно задавать в виде их относительных значений
h N N
А 4 ^, h = е ^, е А , = 1.
h
$= і
«= і
2
Тогда координаты межслойных поверхностей внутри элемента будут определяться в
виде
s -1
z, = - 1, z 2 =- 1 + 2А ,, z з =- 1 + 2А 1 + 2А 2,..., z, = - 1 + 20 А ,Z „, = + 1,
n = 1
а координата середины £ -го слоя определяется как ( £ £ £ +1 ) /2.
В глобальной декартовой системе координат х,у,2 предлагаемый КЭ представляет собой искривленный параллелепипед, у которого верхняя и нижняя поверхности существенно искривлены, а четыре боковые грани являются линейчатыми поверхностями (рисунок 2). Для удобства построения подобного КЭ вводится тройная нумерация узлов в соответствии с рисунком 1.
(3,3,2)
(3,3,1)
(1,3,1)
Рисунок 2 - Тройная нумерация узлов КЭ.
Радиус-вектор аппроксимируется в виде биквадратичного полинома по координатам X , h и линейного — в поперечном направлении, т.е.
N 3 2
r(Xh,Z)= е е е r(x(Ои,(h)Lk(О,
s = 1 m, n = 1 к = 1
где Hn (X ), Hm (h ) и Lk (z ) задаются соотношениями
2
h, (X)= 2 X(X-1), H 2 (X) = l -X2, h з (X)--1 X(x +1), і (Х^-Й2-^, і (0 = /^
2 2 z S + 1 z S z S + 1 z s
Вектор перемещений определяется в виде
N 3 2
V(Xя•{) = е е е '7(Х(~'л'"),("))н.(X)н,(Я)4({■).
£ = 1 т, п = 1 к = 1
что обеспечивает кусочно-линейное изменение перемещений по толщине оболочки. Вводятся ковариантные компоненты деформаций:
3 г . ,3 V 3 г . ,3 V 3 г , ,3 V
Зr.,ЗV Зr .,Зv
Зr .,Зv Зr ..Зі
З r З v З r З
e = _ u_ + _ u_ e = _ u_ + _ u_ с = _ u_ + _ u_
Xh ЗХ Зh Зh ЗХ ’ hz Зh зz її Зh ’ xz ЗХ зz зz зх '
(2)
Связь между приведенными векторами деформаций
{ 0 і1 = { 0 xyx У,с уУУ У»0 z yz y»g xjy y5 g yyz y»g z yxy} , { 0 a } = {c a a , 0 pp ,0 gg ,g ap , g pa , g ga\
записывается аналогично (1).
При вычислении деформаций (2) через перемещения используется метод двойной аппроксимации. В соответствии с этой методикой для каждого * -го слоя деформация е %% вычисляется в точках % = ± і/л/з; Ц==- і, 0, і;
І = І * *+ і, и по этим значениям строится
квадратичная аппроксимация вдоль координаты Ц и линейные вдоль % и І . Деформация вычисляется в точках % = - і,0,і; ц = ± іД/з; І = І *,І * + і, и для нее определяются линейные полиномы по Ц и І , и квадратичные — по % . Деформация обжатия вычисляется в точках % = - і, 0,і; І = І *, что дает квадратичные аппроксимации по % и Ц ,
и постоянное значение по толщине. Деформация 7 вычисляется в точках % = ± іД/з;
І = І * *+ і и определяется в виде линейных аппроксимаций по всем координатам.
Деформация поперечного сдвига 71( вычисляется в точках % = ± іД/з; Ц = - і, 0, і; І = І *, по которым строится квадратичная аппроксимация по Ц , линейная — по % и постоянная — по І . Сдвиговая деформация 7ч( вычисляется при % = - і, 0, і; ц = ± іД/з; І = І * и это дает квадратичный полином по % , линейный по Ц и постоянный вдоль І .
Результатом таких построений будут соотношения
N з 2
= е е е н,(;)н,(Ч)^(і )г(%<ч,ч(-),{<*>)чр(%т,ч")),
*= і і,,,т= і к = і
N з 2
е'„ = е е е н,(і)н„(Ц )ік(І)Г(% 1т\пи,Л,п)ч?(%"),ч ,пл <к)),
*= і і,,,т= і к = і
N з 2
е'і; = е е е Н,(%)н,(Ц)Щ;*)Щ(;)?(гм<")ч?(%«ц1»),
*= і і,, = 11 ,к= і
N з 2
7'1,-- е е е (Р (%) р,Лч) + р, (5) Рт (Ц)) і (І) Г (% ,т).ч(п),(т) ч? (;(і “І№)),
* = і і,,,т = і к = і
ґ Р(%(п>,ц (1>,Іа>) Ч?(% (і),ц(1>,І(к>),
где нт, ( % ) определяется соотношениями
н,!І(Х)= і (і-;'/з) нт (% ) н у (%1,1) + і (і *%ЩН1(1;<1) н у({‘!)),
а Рп, (% ) . — соотношениями
р„(і)== і (і -^Тз) н ((<”) н і (і +^тз) нп(С1))н ї (% <»>).
Потенциальная энергия деформаций вычисляется путем численного интегрирования по всем трем координатам % ,Ц , І в местных координатах ху, у), і), т.е.
Для верификации оболочечного конечного элемента рассмотрим ряд тестовых задач. Используем модели изотропного и ортотропного тела. Для изотропного тела примем модуль Юнга ЕВП200 , коэффициент Пуассона т = 0,3. Для ортотропного тела:
Еі = ЕЩ10 , ЕГ=[і00 , ті2 = т 2з = т зі = 0,0з; Сі2 = Є2з = GJІHdl8,6 .
Приведем сравнение с аналитическим решением и с результатами, полученными в ППП ANSYS. Рассмотрим квадратную пластинку, длина стороны которой ом 0,4 , толщина
Им 2 Чі 0-з . Из условий симметрии рассматривается четверть пластинки. Представим следующие способы нагрузки:
і) Шарнирно опертая пластинка, нагруженная в центре сосредоточенной силой PH 500 . Определяется прогиб № (см) в точке приложения силы Р . Результаты расчета
представлены в таблице і. Сравнение проводится с аналитическим решением [5, 255] и с численным решением, полученным в [з, і5].
Таблица 1
^^теор (см) Сетка 2 ґ 2 з ґ з 4 ґ 4 5 ґ 5 і0ґ і0
Изотропный материал
0.00бзз Рго^гат 0.00645 0.006з9 0.006з6 0.006з4 0.006зз
[з, і5] 0.00бз 0.006з 0.006з 0.006з 0.006зз
ANSYS 0.006586 0.00646і 0.0064і2 0.006з88 0.006з5
Ортотропный материал
0.0і45 Program 0.0і52 0.0і56 0.0і57 0.0і58 0.0і586
[з, і5] 0.0і6з 0.0і67 0.0і69 0.0і7 0.0і7і
ANSYS 0.0і7і88 0.0і6984 0.0і6997 0.0і7027 0.0і7і06
2) Шарнирно опертая квадратная пластинка, нагруженная распределенной нагрузкой н
q = і,25Чі0 —- Определяется прогиб № (см) в центре пластины. Результаты расчета
м
представлены в таблице 2. Сравнение проводится с аналитическим решением [5, 257], с численным решением, полученным в [з, і6], а также с результатами, полученными в ППП ANSYS.
Таблица 2
^^теор (см) Сетка 2 ґ 2 з ґ з 4 ґ 4 5 ґ 5 і0ґ і0
Изотропный материал
0.0087 Program 0.0085 0.0086 0.0087 0.0088 0.0088
[з, і6] 0.0087 0.0088 0.0088 0.0088 0.0089
ANSYS 0.008з78 0.008646 0.008744 0.00879 0.00885і
Ортотропный материал
0.02і4 Program 0.0209 0.02і6 0.02і7 0.02і8 0.02і9
[з, і6] 0.02і6 0.02і8 0.02і9 0.02і9 0.022
ANSYS 0.02і007 0.02і526 0.02і729 0.02і822 0.02і95і
Для определения границ применимости описанной выше модели относительно геометрических размеров и различий жесткостных характеристик слоев, с учетом ортотропной структуры материала проводилось сравнение конечно-элементного решения задачи изгиба трехслойной квадратной шарнирно опертой пластины под синусоидальной нагрузкой со следующими параметрами:
Еі = ЕМ41Ш2тъ , Щ<аШ0ъ G , і2 = ізМП45 Чі0з ,
G23 - 1.38Ч103I I a, m12 = m 23 - m 31 - 0.01, q1 - q0 Sin—Sinpy. a b
Рассматривался вариант укладки слоев трехслойной пластины с углами намотки
у - 0° ,90° ,0° и толщинами И/ 4, И/ 2, И/ 4. Верификация проводилась с другими
приближенными решениями авторов и с использованием ППП ANSYS. В таблице 3 представлены безразмерные величины
w - wQh jq0 И a4, s * - s И2/q0 И2, t * - t h 2!q0 h2,
q ' л 4g12 +(Ei + (i + 2m 12) e2 )/( i -m12m 21)щ/12,
где a - линейный размер пластины, И - толщина пластины.
___________________________________________________________________________Таблица 3
a / И Тип решения * w (a/ 2, a/ 2, И/ 2) * s x (a/ 2, a/ 2, И/ 2) * t xy (0,0, И/ 2)
100 Program 1,0046 0,502 0,0209
[3, 17] 1,06 0,482 0,0196
[6, 100] 1,008 0,539 0,0214
[7, 77] 1,003 0,566 0,0223
Исходя из полученных результатов, следует отметить, что разработанная численная методика исследования напряженно-деформированного ортотропных и многослойных пластин и оболочек дает результаты, хорошо согласующиеся с теоретическими значениями и результатами других авторов. Следовательно, на ее основе можно рассчитывать подобные конструкции и получать достоверные результаты.
Работа выполнена при частичной финансовой поддержке РФФИ в рамках научных проектов № 12-01-00955, № 12-01-97026, № 12-01-31212, № 13-97057, № 13-01-97058.
Литература
1. Соловьев С.С. Конечно-элементная модель многослойной оболочки с анизотропными слоями переменной толщины. // Известия Вузов. Авиационная техника. 1989, №4, с. 71-75.
2. Рикардс Р.Б. Метод конечных элементов в теории оболочек и пластин. - Рига: Знание, 1988, 284с.
3. Бережной Д.В., Сагдатуллин М.К., Голованов А.И. Многослойный ортотропный конечный элемент оболочек средней толщины. // Вестник Саратовского государственного технического университета. - 2011. -№3 (57). Выпуск 1. - С. 9-19.
4. Бурман Я.З., Соловьев С.С. Расчет упругопластического деформирования оболочек на основе теории течения и МКЭ. Исследования по теории пластин и оболочек. - Казань: Изд-во КГУ, 1990, №22, с. 98-107.
5. Шлычков С.В. Оценка точности расчетной модели напряженно- деформированного состояния тонкостенных элементов музыкальных струнных инструментов. // Исследовано в России. - 2000. - С. 245-262.
6. Пагано Н., Хэтфилд С. Упругое поведение многослойного двунаправленного композиционного материала // Ракетная техника и космонавтика. - 1972. - Т. 10. - № 7. - С. 98 - 101.
7. Panda S., Natarajan R. Finite element analysis of laminated composite plates // International Journal for Numerical Methods in Engineering. - 1979. - V. 14. - № 1. - P. 69 - 79.
