CC BY
49
Вычислительные технологии
Том 13, № 4, 2008
Расчет оболочки вращения при произвольном нагружении с
использованием МКЭ на основе функционала Рейснера
Н.А. Гуреева, Ю.В. Клочков, А. П. Николаев Волгоградская государственная сельскохозяйственная академия, Россия e-mail: natalya-gureeva@yandex. ru
An eight nodal volumetric finite element with nodal displacement and stress as unknown parameters is developed for calculation of an arbitrary loaded shell of rotation. Approximations for stress and displacement were constructed using three linear relations. The modified stiffness matrix is obtained with the help of a variational principle based on the Reysner functional. We present an example, which shows that the proposed finite element method has a good potential for calculations of plates and shells of any thickness.
Введение
Для расчета произвольно нагруженной оболочки вращения с использованием криволинейной системы координат разработан объемный шестигранный конечный элемент в смешанной формулировке, узловыми неизвестными которого являются компоненты вектора перемещений и тензора напряжений. Каждая компонента внутренней точки конечного элемента аппроксимировалась через узловые значения трилинейными соотношениями. При получении матрицы деформирования конечного элемента использовался вариационный принцип на основе функционала Рейснера.
1. Геометрия оболочки вращения.
Радиус-вектор произвольной точки срединной поверхности (рис. 1) определяется выражением
где і, ,і, к — орты декартовой системы координат; х — осевая координата; г — радиус вращений точки срединной поверхности; в — угловая координата.
Базисные векторы в точке срединной поверхности, касательные к ней, определяются дифференцированием (1.1):
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.
R = xi + r sin ej + r cos ek,
(1.1)
a1 = R,s = x,s i + r,s sin ej + r,s cos ek; a2 = R,6> = r cos ej — r sin ek.
(1.2)
Рис. 1
(1.5)
Нормаль к срединной поверхности определяется векторным произведением
a3 = -a1 х a2 = — r,s i + x,s sin 9j + x,s cos 9k. (1.3)
2
Соотношения (1.2) и (1.3) можно представить в матричном виде
{a} = [M]{i}, (1.4)
3x1 3x33x1
где {a}T = {a1, a2, a3} — вектор-строка базисных векторов; {i}T = {i, j, k} — вектор-строка единичных векторов.
Производные базисных векторов произвольной точки срединной поверхности через орты декартовой системы координат имеют вид
ai,s = x,ss i + r,ss sin 9j + r,ss cos 9k; a2,s = r,s cos 9j — r,s sin 9k; a3,s = —r,ss i + x,ss sin 9j + x,ss cos 9k; a1,e = r,s cos 9j — r,s sin 9k; a2,0 = —r sin 9j — r cos 9k; a3,6> = x,s cos 9j — x,s sin 9k, которые можно представить в матричном виде
{a,s} = [N1]{i}; {a,e} = [N2] {i}, (1.6)
где {a s}T = {a1s, a2 s, a3 s} — вектор-строка производных базисных векторов по длине
дуги; {a,#}T = {a1;0, a2;^, a3;^} — вектор-строка производных базисных векторов по уг-
ловой координате.
Введем матричное соотношение, обратное (1.4),
{i} = [M ]-1 {a}. (1.7)
Выражения производных базисных векторов через векторы этого же базиса запишутся в виде
{a,s } = [N1] [M]-1{a} = [m]{a};
{a,0 } = [N2][m]-1{a} = [n]{a}.
(1.8)
2. Перемещения и деформации
Вектор перемещения произвольной точки оболочки выражается через компоненты локального базиса в виде
V = V :аі + ^2а2 + ^3а3, (2.1)
где V1, V2, V3 — проекции вектора перемещения на векторы локального базиса. Производные вектора перемещения по координатам Б, в, і имеют вид
V* = /і1аі + /іа2 + /іаз;
V,в = /2і аі + /2 а2 + /23а3; (2.2)
V,* = Vі а і + ^2а2 + V3 а3,
где /і, /2І (і = 1, 2, 3) — коэффициенты, содержащие компоненты вектора перемещения и их производные.
Радиус-вектор точки, отстоящей на расстояние і от срединной поверхности (рис. 1), определяется выражением
И* = И, + іа, (2.3)
а радиус-вектор этой точки в деформированном состоянии записывается суммой векторов
И** = И* + V. (2.4)
Векторы, касательные к произвольной поверхности в исходном и деформированном состояниях, определяются дифференцированием (2.3) и (2.4):
& = = аі + іа3,і; ( )
g* = и,** = ^ + V,і. ( . )
Компоненты тензора деформаций, определяемые разностью компонент метрических
1 (о*
2
тензоров Є. = 2 (д*. — д.) [1], запишутся соотношениями
Є 11 = g і ■ V,*; Є22 = g2 ■ V,в ; Єзз = gз ■ V,*;
є 12 = 2^ 1 ■ V,в +g2 ■ V,* ) = Є2 1;
1 (2.6) Є13 = 2^і ■ V,* +gз ■ V,*) = Є3 1;
Є23 = 2^2 ■ V,* +gз ■ V,* ) = Є32.
С использованием (1.8) векторы gi (2.5) можно представить выражениями
g 1 = ^11а1 + ^1 2а2 + ^13а3; g2 = ^2 1 а1 + ^22а2 + ^23а3І
gз = а3.
(2.7)
Учитывая (2.7), выразим (2.6) через компоненты вектора перемещения
Є11 = ^ А11 + ^2 А12 + ^3 А13 + V1А14 + ^А15 + 3А16;
Є22 = V,1 А21 + V,2 А22 + V,3 А23 + V1А24 + V2 А25 + V3А26;
Є33 = V,3 Азі;
Є12 = Є21 = ^0 В11 + ^0 В12 + V,3 В13 + ^ В14 + V,2 В15 + V,3 В16 +
+ V 1ві7 + V 2ві8 + ^Ві9;
Є13 = Є31 = 2(^* В21 + В22 + V,3 В23 + V1В24 + V1В25 + ^В26 + ^В27 + ^3 );
Є23 = Є32 = 2(^* В31 + В32 + В33 + v 1В34 + v2^35 + v3^36 + v^ ),
(2.8)
где коэффициенты Ац...Вэб — функции координаты £ компонент метрического тензора и элементов матриц [т], [п]. Например, Аи = (1 + £т31), А12 = г2£т32, А13 = £т33,
А14 = ти(1 + £т31) + £г2т12т32 + ^13т33 и т. д.
Или в матричном виде
1x12
в вышеуказанной точке; [Д] — матрица дифференциальных операций.
6x12
3. Соотношения между деформациями и напряжениями
Закон Гука для изотропной среды в главных осях тензора деформаций и тензора напряжений выражается соотношениями [1]
{Є} = [Д] {V},
(2.9)
6x1 6x1212x1
т
где (е) = |е11е22е332е122е132е23} — вектор-строка деформаций точки, отстоящей на
1x6
расстояние £ от срединной поверхности; (V}т = (V V2v3} — вектор-строка перемещений
(3.1)
Рис. 2
Представим выражение (3.1) в матричном виде
И = ИМ,
6x1 6x66x1
(3.2)
Г тТ г 11 22 33 12 13 23Т и и
где {а) = {ааааа а23) — вектор-строка напряжении в произвольном точке
1x6
оболочки; [С] — матрица упругости материала.
6x6
Функционал РеИснера Пд для оболочки вращения имеет вид [2]
П
{а)Т[Д]{и) - 1 {£}Т{а)
— I {а*)Т ({и) — {и*)) ^з,
(3.3)
где V — объем оболочки; (д*}, (^*) {и*} — заданные перемещения; зи, лами и перемещениями.
заданные поверхностные и граничные силы; части поверхности оболочки с заданными си-
4. Матрица деформирования конечного элемента
Конечный элемент выбран в виде объемного шестигранника с узлами і, ^', к, /, т, п,р, к (рис. 4, а). Для выполнения численного интегрирования глобальные координаты з, 0, і аппроксимировались через локальные координаты куба £, п, С (рис. 4, б) трилинейными соотношениями
А = {/(£,П,С )){Ау ), (4.1)
1x8 8x1
где {Ау}Т = {АгА^АкАгАтАгаАрА^) — вектор-строка узловых значений величины А. Под
1x8
символом А понимаются координаты з,0,і.
Дифференцированием (4.1) определяются производные глобальных координат в локальной системе $,£ £,п $,£ 0,5 0,п 0,с і,5 і,5 і,п і,С и локальных координат в локальной системе £,8 £,# £,* п,5 п,0 п,* С,« С,0 с,*.
Принимая в качестве узловых неизвестных перемещения и напряжения, представим векторы узловых неизвестных в виде
К }т =<! Мт Ют Ют
1x24
yJ I у) I у_ 1x8 1x8 1x8
(V ^17V1кV“V1тV ^ ^1Л ... v3Vpv3h}; т
— < < гг , ( .. , ( а , ( а , ( а , ( а
уууууу
(4.2)
т 11 т 22 т 33 т 12 т 13 т 23 т
(ау} = \(ау } (ау } (ау } (ау } (ау } (ау } / -
1x48
г 11г 11.7 11к 111 11т 11га 11р 11Ь. 23га 23р 23^1
= (а а а а а а а а ••• а а а },
где использованы матричные обозначения
(vy}т = (V ^17 V1к V11V ^ ^ ^1^};
(vy2}т = (V ^ 27 V 2к V21 v2mv2rav 2^};
(vy3}т = (V -V 37 V 3к V31 V 3^ 3> 3^};
33 (ау }
11 т (ау }
22 т (ау }
331 т у
12 т (ау }
13 т (ау }
23 т (ау }
Г_11^ 11.7' 11к 11/ 11т^ 11га ^ 11р^ 11М
(а а а а а а а а }
Г _22^227'22к _221 22т _22га 22п 22^1
(а а а а а а а а }
г 33г 337 33& 331 33т 33га 33р 33Ь,1
(а а а а а а а а }
г 12г 127 12к 121 12т 12га 12р 12^1
(а а а а а а а а }
г 13г 137 13к 131 13т 13га 13р 13^1
(а а а а а а а а }
Г _23г _237'23^ _231 23т _23га 23п 23^1
(а а а а а а а а }.
Используя для аппроксимации перемещений и напряжений внутренней точки конечного элемента трилинейные соотношения (4.1), сформируем матричные соотношения
где
[С]
6x48
(а} = [С] (ау};(и} = [А] (иу}
6x1 6x4848x1 3x1 3x2424x1
(е} = [Д](и} = [ Д ] [А] (иу} = [В] (иу},
6x3 3x1 6x3 3x24 24x1 6x24 24x1
-(£}т
(4.3)
(£}
т
(£}
т
(£}
т
(£}
т
(£}
т
(1 х 8) (1 х 8) (1 х 8) (1 х 8) (1 х 8) (1 х 8)
тт ..................
т
[А] = ... (£}
3x24 ......................... (£}т
(1 х 8) (1 х 8) (1 х 8)
Здесь (£}т = (/(£,п,С)}
1x8 1x8
Т
Для формирования матрицы [В] используем производные компонент вектора перемещений
V,! = / П,а +(/,п }ТП,а +(/,С }ТС,Л ](*£} (2 =1, 2, 3),
где под символом Л понимаются координаты 5, 0, £.
С учетом матричных соотношений (4.3) функционал (3) запишем в виде
Пя = (ау}т / [С]т [В] ^[иу] - 2(ау}т[С] [С] ^(а}-
1x48 и 48x6 6x24 48x1 2 1x48 6x66x48 48x1
V
— (иу}Т / [А] т(9*}^-
1x24 и 24x3 3x1
(4.4)
В (4.4) заданные перемещения приняты равными нулю.
Дифференцируя равенство (4.4) по узловым неизвестным (ау}т и (иу}т, получим систему уравнений
дп
дт-Ят = [Н] (ау} + [3] (иу} = 0; д (ау }
у
дПя
48x48 48x1 48x24 24x1
д(и }т = [3] Т(ау} — (/} = 0
д (иу} 24x48 48x1 24x1
(4.5)
где
[3]
24x48
[С] т [В] ^; [Н] = [С] т [С] [С] ^;
IV 48x6 6x24 (/}
48x48
[А] т д*
IV 48x6 6x6 6x48
24x1 J sa 24x3 3x1
Запишем систему (4.5) в традиционной для МКЭ форме
[к] [Яу] = [^] ,
72x7272x1 72x1
где
[к]
72x72
(Яу }
1x72
т
— [Н] [3]!
48x48 48x24
[3] т [О]
.48x24 24x24.
(ау }Т(иу}Т
1, 1x48 1x24
матрица деформирования конечного элемента;
(^}т = I (0}т (/}т — вектор узловых усилий элемента.
1x72 [1x48 1x24
При формировании матрицы деформирования всей конструкции используется традиционная процедура МКЭ [3].
В качестве примера рассчитана находящаяся под действием равномерно распределенного давления интенсивности д оболочка вращения, меридиан которой описывается параболической зависимостью (рис. 4):
г = ^ х2 + А.
I2
тт
г 4'-
Рис. 3
Приняты следующие исходные данные: А = 0.25 м; Д = 0.02 м; I = 0.5 м; Е =
2 ■ 105МПа; V = 0.3.
даН
Расчет выполнен в двух вариантах. В первом варианте при д = 10---------^ толщина
см2
оболочки принималась равной к = 0.004 м (оболочка считалась тонкой). Численные значения меридиональных а33 и кольцевых иве напряжений в зависимости от количества
( даН \
элементов дискретизации пэ приведены в табл. I в ------— .
V см2 /
Анализ табличных результатов показывает хорошую сходимость вычислительного процесса. Численные значения меридиональных напряжений а33 в точках 3 и 4, близкие к нулю, соответствуют физическому смыслу задачи. Средние значения напряжений в точках I и 2, полученные из условия равновесия оболочки, равны
^ = д(АЗ - д2> =308.07 даНН,
2(А + к/2)к см2
что примерно на I % отличается от табличного значения.
даН
Во втором варианте расчета при д = 100----— толщина оболочки принималась равной
см2
к = 0.07 м (оболочка не являлась тонкой).
Численные значения меридиональных напряжений в зависимости от числа элементов дискретизации приведены в табл. 2.
Т а б л и ц а 1
Пэ = = 56 Пэ = 100 Пэ = 212
№ точки (Ш (Ш &вв
1 305.58 485.78 305.44 485.44 306.11 485.46
2 311.16 480.39 308.22 479.22 307.69 478.88
3 4.6 84.11 0.92 81.4 0.27 81.02
4 1.12 54.69 1.42 53.42 1.65 53.49
Т аблица2
Пэ = 56 Пэ = = 100 Пэ = 212
№ точки & вв &вв & вв &вв & вв &вв
1 138.44 315.72 136.11 315.13 136.14 315.21
2 159.53 248.06 159.08 248.07 158.06 247.95
3 32.44 276.84 2.16 277.76 3.13 278.27
4 5.77 -12.4 0.75 -14.07 0.37 -14.65
Анализ табличных результатов свидетельствует о хорошей сходимости вычислительного процесса.
В обоих вариантах число элементов по толщине оболочки принималось равным четырем. В данной задаче и при двух элементах по толщине оболочки получились практически одинаковые результаты.
Исходя из изложенного можно сделать вывод, что разработанный элемент вполне пригоден для расчета пластин и оболочек произвольной толщины.
Список литературы
[1] Седов Л.И. Механика сплошной среды. М.: Наука, 1976. Т. 1. 536 с.
[2] Галагер Р. Метод конечных элементов. Основы: пер. с англ. М.: Мир, 1984. 428 с.
[3] Постнов В.А., Хархурим И.Я. Метод конечных элементов в расчетах судовых конструкций. Л.: Судостроение, 1974. 344 с.
Поступила в редакцию 14 мая 2007 г., в переработанном виде — 3 марта 2008 г.
