CC BY
36
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «ИННОВАЦИОННАЯ НАУКА» №6/2016 ISSN 2410-6070_
УДК 629.7.036.54
К.П. Баслык
к.т.н, доцент кафедры «Космические аппараты и ракеты-носители», Московский государственный технический университет им. Н.Э. Баумана
Д.Х. Зарипов студент группы СМ6-122 Факультет специального машиностроения, Московский государственный технический университет им. Н.Э. Баумана,
г. Москва, Российская Федерация CA. Чурилин инженер отдела СМ1-1 Научно-исследовательский институт специального машиностроения Московский государственный технический университет им. Н.Э. Баумана,
г. Москва, Российская Федерация
ОСОБЕННОСТИ ЧИСЛЕННОЙ РЕАЛИЗАЦИИ АЛГОРИТМА ПОСТРОЕНИЯ КОНЕЧНОГО ЭЛЕМЕНТА КОМПОЗИТНОГО ТЕЛА ВРАЩЕНИЯ ДЛЯ РЕШЕНИЯ НЕЛИНЕЙНОЙ ЗАДАЧИ С ИСПОЛЬЗОВАНИЕМ МОДИФИЦИРОВАННОГО ПОДХОДА ЛАГРАНЖА
Аннотация
Рассмотрен алгоритм построения конечного элемента тела вращения из композиционных материалов для решения нелинейной задачи в осесимметричной постановке. Используется модифицированный подход Лагранжа. Приведены зависимости для определения коэффициентов термоупругости в различных системах координат. Представлено решение тестовых примеров: определение критического внешнего давления сферической оболочки и определение напряженно-деформированного состояния ортотропной цилиндрической оболочки.
Ключевые слова
Метод конечных элементов (МКЭ), конечный элемент (КЭ); тело вращения; композиционный материал (КМ); нелинейная задача, модифицированный подход Лагранжа.
Используемые сокращения:
СК - система координат; ГСК - глобальная система координат; ЛСК - локальная система координат; МЖЭ - матрица жесткости элемента; ВУС - вектор приведенных узловых сил; МЖК - матрица жесткости конструкции; ЦТ - центр тяжести; ТКН - температурный коэффициент напряжения;
КЛТР - коэффициент линейного термического расширения.
Модифицированный подход Лагранжа (Updated Lagrangian) [1,3,8] является одним из методов решения задач теории упругости в геометрически и физически нелинейной постановке. Достоинство подхода состоит в том, что в его формулировке присутствуют достаточно простые деформационные соотношения -модифицированные приращения деформаций Грина, а также используются коэффициенты упругости, соответствующие текущему деформированному состоянию. При получении результата не требуется решать обширную систему нелинейных уравнений.
Отмечается (см. [8]), что данный подход может быть использован, когда деформации не превосходят 0,04. Предельные деформации для композиционных материалов меньше этого значения, поэтому использование модифицированного подхода Лагранжа представляется вполне корректным при нелинейном анализе таких конструкций.
МЖЭ и ВУС четырехузлового изопараметрического КЭ тела вращения с 8-ю степенями свободы вычисляются в ГСК. Для аппроксимации приращения узловых степенй свободы используются нормированные координаты четырехугольника [6]. Данный КЭ может применяться для решения геометрически и физически нелинейных задач тел вращения в осесимметричной постановке.
МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «ИННОВАЦИОННАЯ НАУКА» №6/2016 ISSN 2410-6070 1. Построение МЖЭ и ВУС четырехузлового изопараметрического КЭ тела вращения
Рассмотрим квазидвумерный КЭ тела вращения с четырехугольным сечением (рис. 1).
Рисунок 1 - Четырехугольный КЭ тела вращения: а - в различных СК, соответствующих разным деформированным состояниям; Ь - в нормированной СК
Для описания геометрии, кинематики деформирования и деформационных соотношений введем несколько ортогональных криволинейных, вообще говоря, систем координат, используя обозначения, принятые в [1]:
ОХ1Х2 - ГСК, соответствующая деформированному состоянию П(0); ОХ1Х2 - ГСК, соответствующая деформированному состоянию О^-1; 0717 - ГСК, соответствующая деформированному состоянию О(№1);
О] - ЛСК, связанная с узлами конечного элемента в каждом из деформированных состояний;
(^И-1,1].
Ось («)з - это суть окружная координата, она дополняет указанные системы координат до правой тройки; ось вращения ориентирована горизонтально.
Будем считать, что переход от деформированного состояния к деформированному состоянию осуществляется на шаге нагружения номер «N+1».
1.1. Геометрические соотношения В деформированных состояниях П(0), и О(№1) координаты любой точки изопараметрического четырехугольного КЭ (см. рис. 1) определяются как
4 , л 4
х=х / (ы*1(г); *2=1 / (#,т)*2(0, (хХ,7)- (1)
'=1 ¡=1
Здесь (х(^)), (х(^Х()), (У\У^')), (' = 1,..,4) - координаты узлов КЭ в каждом из перечисленных деформированных состояний;
= (1 -£)(! -]) ; _(1 + £)(1 -]); =(1 + ^)(1 + ]); =(1 -^)(1 + ]) (2)
/ 4 ' 2 4 ' /3 4 ' 4 4
- функции формы, записанные относительно нормированных координат Е,и т] (см. [6]). Выражения (1) для удобства следует переписать в матричном виде. Так для переменной Х1 получаем:
X! = N • ; X = N • ; У1 = N • У1, (1,2), (3)
где N = [/ / /з /4 ] (4)
- матрица функций формы аппроксимации геометрии;
X =[х(1) х<2) х<3) х<4) ]", (х,Х,¥); (ХХ,7); (1,2) (5)
- векторы и соответствующие координаты.
ÖX1
Частные производные вида-, (1,2); вычисляются с использованием соотношений (2) и (3), а
именно:
ÖXi = дк X(1)+к X (2)+к X(з)+к X (4) =ÖN. X c^ c^ 1 c^ 1 c^ 1 c^ 1 c^ 1:
(j
(6)
где ^ = 1 [_ 1 + л 1 1 + ^ -1 — = 1 [-1 + £ -1 1 + £ 1 _£]. (7)
Тогда компоненты матрицы Якоби для деформированных состояний О(0), О1^1 и могут быть
вычислены как
Ш
(0) J =
ÖXj Öx2 ÖN
ÖF "dT
ÖXj Öx2 ÖN
ö j Ö j Ö j
cn
' Х1 — ~ * Хо
ÖN
■* x,
"ÖXj ÖX2'
; (N) J = ÖX Ö j ÖX2 Öj =
ÖN. x^ — - X 2
ÖN,x dN.x
dj 1 dj 2
(N+1)
dN * Y dN * Y
d^ 1 d^ 2
dN * Y dN * Y
д j д j
Теперь могут быть определены частные производные
"ÖY1 дУ2'
5Y, дУ2
Ö j Öj
(8)
dYt дХ,
■, (/,/=1,2). Опуская промежуточные
выкладки, получаем:
"Öy " "ÖN * y 1
dX = ( (N) J у * 1
dYj ön*Y
ÖX2 Ö j
ÖY2 dXj ÖY2 X
:( (N ) J У *
dN * y2
2
dN * Y
д j
(9)
Соотношения (9), очевидно, позволяют вычислить все требуемые частные производные в любой точке КЭ.
1.2. Кинематика деформирования Вектор приращений перемещений точки КЭ с координатами (Х1Х2) на шаге нагружения имеет вид:
ли = [А^ ди2 ]т. (10)
Все используемые в дальнейшем переменные могут зависеть от двух координат: Х1 и Х2, (хХ,У), индексы которых указываться не будут. Индексы (М\и) и (^+1)(») при проведении выкладок также, кроме некоторых отдельных случаев, не указываются. При этом (см. [1]):
X = хг + и; ^ = Х1+Ди1; Х2 = х2 + и2; 72 =Х2+Ди2. (11)
Здесь и = [и и2 ]т, - вектор декартовых перемещений, переводящий КЭ из состояния О(0) в
деформированное состояние О^; ли = [ли1 ди2 ]т - вектор приращений перемещений, переводящий
КЭ из состояния в состояние О(№1).
1.3. Деформационные соотношения, модифицированный вектор приращения деформаций Грина
Для компонент модифицированного вектора приращений деформаций Грина запишем, как указано в [1, 8]:
Ае„ = Деп +Д j =
ÖAU 1 ■ + ■
ÖX 2
(ял тт \2 ^ ял tt \2
dAU
v^1 у
+
cau2 v^Ty
A*e22 = A*£22 + A* j22 =
ÖAU, 1 -2 + —
ÖX2 2
ÖAU
л2 r
KÖX2 у
+
ÖAU,
v dX 2 у
(12) (13)
x
2
2
2Д*е12 = 2Д*£12 + 2Д*^ =
дДЦ дДи,
дХ,
дХ1
■ +
дДЦ дДЦ дДЦ dAU2 дХ дХ2 +
дХ дХ
л * л * л * Д U^O 1
Д^3 = Д^33 +Д^33 + -
л, 2
ДЦ v Х 2 у
(14)
(15)
Тогда можно сформировать векторы:
Д*е=[Д*еп А*в22 Жв12 Д*езз]т; (16)
Д*Е=[Д*£п ДД£22 2Д£ Д£зз]т; (17)
Д*П=[Д>1 Д'т 2Д-ГЦ2 Д>з]т, (18)
где Д*е - модифицированный вектор приращений деформаций Грина; Д*Е - линеаризованная составляющая модифицированного вектора приращений деформаций Грина; Д*л - нелинейная составляющая модифицированного вектора приращений деформаций Грина.
1.4. Соотношения упругости
Связь между модифицированным вектором приращений напряжений Кирхгофа и модифицированным вектором приращений деформаций Грина после линеаризации записывается как (см. [1])
Д*о = С* • Д*Е - р*ДГ, (19)
I # * * * Т
где До = [Д ах 15 Д а22, Д аи, Д а33\ (20)
- модифицированный вектор приращений напряжений Кирхгофа; ДТ - изменение температуры на шаге нагружения.
Для случая изотропного материала матрица коэффициентов упругости и вектор ТКН имеют вид:
C* =
с11 * С12 0 * С13
* С21 * С22 0 * С23
0 0 * С44 0
* С31 * С32 0 * С33
р*=д Д Д; Д;]1.
(21)
как
Здесь коэффициенты термоупругости (индекс «*» в дальнейшем указываться не будет) вычисляются (1 -у)Е уЕ
С11 С22 С33
(1 + v)(1 - 2v)'
(1 + v)(1 - 2v)
E
ß11 ß22 ß33
E-a
ßl2 = 0 ,
(22)
2(1 + у)' "" (1 - 2у)
где Е, у, а- модуль Юнга, коэффициент Пуассона и КЛТР материала.
Также считаем, что матрица С и вектор Р - это коэффициенты термоупругости материала в текущем
деформированном и температурном состоянии
1.6. Аппроксимация перемещений и необходимых производных
Для аппроксимации вектора приращений обобщенных перемещений Ди (см. формулу (11)) введем вектор приращений узловых степеней свободы КЭ:
Дд* = [Дц(1) Ди® ДЦ(2) Ди(2) Ди<3) Ди<3) ДЦ(4) ди<4) ]
Тогда Ди =
ДЦ ДЦ
Ф - Дц
2
/ 0 /2 0 /3 0 /4 0
0 /1 0 /2 0 /30 /4
Дц.
(23)
(24)
При построении КЭ также необходимо располагать значениями первых производных от функций
сц dfi
формы по координатамХ1 иХ2: -¿±-; -,(/'=1,...,4)). Тогда целесообразно сформировать матрицу вида:
сх сх0
2
Г Sfi Sf2 S/з s/4
SX, SX, SX, SX
( dФ ^ S/i S/2 s/з S/4 . (25)
LdXi2 J SX2 SX2 SX2 SX2
fi /2 /з /4
L X 2 X 2 X 2 X 2
Для этого сначала вычислим производные от функций формы по нормированным координатам:
S/ _ S/ SX Sf SX . — 1 ; f = S/ SXi + S/ SX2 , (Z=1,^,4). (26)
s^ sx s^ SX 2 S^ SX, s^ sX2 s^
Теперь аналогично соотношению (12) с использованием формулы (26) получим:
Sf Sf Sf f • (((N) J)- 7 . (27)
SX, SX2 [S^ S^ \\ ' /
Этих соотношений, очевидно, достаточно для вычисления компонент матрицы (28) в любой точке нормированной области
1.7. Вариационное уравнение равновесия изолированного КЭ и аппроксимация его составляющих
Теперь необходимо задать аппроксимации слагаемых из используемой вариационной формулировки (см. [1,3]), где за исходное принимается деформированное состояние Для рассматриваемого здесь КЭ тела вращения она имеет вид:
11 11
к*ъ(тл.((тЛ1К1, г г ел* т
•о я •(N) г
det ((N ) J )d£d^
JjSA * E T C* •А * E^(N) г det ((N) i)d&V + \\öA * n -i-i -i-i
i i
JJSA * E T •ß * •Ar •(N) г det ((N) -i -i
ii ii
J J SAU T Ag^(N) г det((N) J- jj SAU T g(N) г det((N) J+
-i-i
-i-i
i i
+ JjSA * E T о я •(N) г det ((N) J )d^-SAq
-i -i
(N) г • (N+i) At '(i) At (i)
(N) г • (N+i) At
г(2) At (2) (N) г • (N+i) At
г(3) At (3)
(N) г • (N+i) At
г(4) At (4)
- SAq
(N)
г • t
r(i) 1 (i)
(N)
(N)
(N)
(2) (2)
t
t
(3) (3)
t
г
(N)
-SAq
T (N )i
' Ca
2
(N+i)
AP(Sim + S )
(N+i)
AP(Sim + Sin ) (N+i) Ap(Sim + S)
(N)
-SAq
T (N)
Ca
2
(4) (4)
P(Sim + Sin )
P(Sim + Sin )
P(Sim + Sin )
P(Sim + Sin )
=0.
(28)
Здесь 0е - вектор напряжений Коши-Эйлера, накопленный по результатам предыдущих шагов нагружения; Ag=[Agl Ag2]т - приращение вектора объемных сил на шаге нагружения «N+1»; g=[Agl Ag2]т -вектор объемных сил, накопленных к началу этого шага нагружения; (№1)Д1(г-), (/—1,.. .,4) - приращение вектора узловых реакций, приложенных к узлу номер «/» КЭ; 1®, (/=1,..,4) - вектор узловых реакций, накопленных по результатам предыдущих шагов нагружения и приложенных к узлу номер «/» КЭ; Др=[Др1 Др2]т -приращение вектора распределенных сил в ГСК, приложенного к узлу номер «/» КЭ; р©=[р1 Р2]т - вектор распределенной нагрузки, накопленной по результатам предыдущих шагов нагружения и преобразованный к ГСК. Вектор р(/), действует в узле номер «/» КЭ; т и п - локальные номера узлов конечного элемента, идентифицирующие нагруженную сторону; между ними действуют распределенные силы Др, р;
T
Smn - символ Кронекера; Lo и reo - длина и радиус окружности, проходящей через середину нагруженной стороны КЭ, взятые в деформированном состоянии Q(N).
Отметим, что векторы (N+1)At(;> и t(¿) содержат погонные, то есть распределенные по окружностям
величины, и поэтому умножаются на величины соответствующих радиусов вращения (28)).
АЕ = 'Ли = ФАЯ = в АЯ.
(N)
Г^ (см. формулу
Здесь
B = L Ф
а 0 а 0 т
ах, ах2 Г/1 0 /2 0 /3 0 /4 0
0 а ах2 а äX" 1 X2. 0 /1 0 /2 0 /3 0 /4
Г / dX, 0 d/2 dX 0 d/3 dXx 0 d/4 dX, 0
= 0 df dX, 0 d/2 dX, 0 #3 0 d/4 dX2
f df d/2 d/2 d/3 d/3 d/4 d/4
dX 2 dX, dX2 dX, dX 2 dXx dX 2 dXx
0 f1 0 /2 0 /3 0 /4
X2 X 2 X 2 X 2
С учетом формул (29) и (30) получаем:
ЗА * Е т С* -А * Е = З(В -Ац)т • С* (В -Ая) = ЗАя т • В т СВ Ац .
ЗА * Е т -в *АТ = З(В -Ац)т -в *АТ = ЗАц т В т в *АТ.
ЗА* Е т о Е = З(В -Ац)т о Е =ЗАц т В т о Е.
Для составляющих из формулы (33), где присутствует вектор А*^г, имеем:
ЗА* пт оЕ = ЗА* г]п - сгЕ + ЗА* щ2 - аЕ2 + 2ЗА* - аЕ + ЗА * щъ - а3Е3.
£4 * Е
Для слагаемого ЗА Т]п - СТ11 получим:
3A Ли - О = -
ац ац ац2 ац2 ах ах
ах ах1
^ 1 "ац ац2" "ац ац"
J"2 SXj SXj ' ах
ац d/1 0 d/2 0 d/3 0 d/4 0
ах 1 dX, dX, dX, dX,
ац 2 0 d/1 0 d/2 0 d/3 0 /4
а^. dX, dX, dX, dX,
Тогда
3A П11 -ofi T
dX 0
Aq.
dX,
dX 0
dX,
¿X 0
dX,
dXj 0
dX
x
d/1 0 #2 0 #3 0 /4 0
dX, dX dX dX
0 d/1 0 0 #3 0 d/4
dX dX dX dXt
(29)
(30)
(31)
(32)
(33)
(34)
(35)
(36)
Aq О
E —
1
= ^ 3Aq T A T A oE Aq = ^ 3Aq T A OiE Aq;
1
T
T
2
3
4
0
0
0
0
X
2
3
4
МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «ИННОВАЦИОННАЯ НАУКА» №6/2016 ISSN 2410-6070
Аналогично формулы (37) получаем:
SA • a22 —
su su su su2 SX 2 sx 2 + sx 2 sx 2 J
1 SU SU 2 su SU 2
2 SX i , SX 2 SX 2 , SX 2
= i ^SAq T 2
d/i 0 d/2 0 d/3 0 d/4 0
dXn dX2 dX2 dX2
0 d/i 0 d/2 0 d/3 0 d/4
dX2 dX2 dX2 dX2
d/i 0 df2 0 d/3 0 d/4 0
dX2 dX2 dX2 dX2
0 d/i 0 df2 0 d/3 0 d/4
dX2 dX2 dX2 dX2
Aq •a
=^ SAqT A 2T A 2a2E2Aq — isAqT A 22a2E2Aq;
И, наконец, для «сдвиговой» компоненты: Г
.E
'i2
2SA*^i2 а1г —
sau sau sau2 sau2
vSXi SX 2
SXi SX 2 ,
a
-E = 12=
SAU SAU
sx sx
SAU sau2
sx SX2
+
SAU SAU
SX, SX,
SAU SAU
SX SX
тЛ
a
.E = i2
(38)
"SAqT (AiTA2 + A2TAi— ^SAqT(Ai2 + A2i)ai2Aq .
l2 1 "^2 iri 2^4 ^ ^ V^ i 2 Для окружных деформаций и напряжений следует выражение:
i Г au V i
2
(
SA *tf33aE3 —
2
v X 2 J
Г AU ^ T AU
Л
v X 2 J
X
1 л T
2 Aq
0 а
X,
0 а X
А
X,
(39)
f X,
0 а 0
X.
Ä X.
0 f 0 ii-XX
Здесь
Ai i —{4,};
л —
a22 — {5,};
A 12+A2 i— {p,};
в —
i
A33 — {p, );
P —
bkb,, i — 2k - i, j — 2/ - i; 0, i — 2k - i, j — 2/; bkbt, i — 2k, j — 2/; 0, — 2k, j — 2/ - ;
akb + atbk, i — 2k - i, j — 2/ - i;
0, i — 2k - i, j — 2/;
akbi + a/bk, — 2k, j — 2/; 0, i — 2k, j — 2/ - ;
0, i — 2k -1;
P —
|ЧС/, i — 2k, j — 2/; 0, i — 2k, j — 2/ - ;
Aq • a3E3 =iSAqT A3T A3aEAq — iSAqT A3Aq . (40)
akat, i — 2k - i, j — 2/ - i; 0, i — 2k - i, j — 2/; a a , i — 2k, j — 2/; 0, i — 2k, j — 2/ - ;
T
T
x
x
T
T
0
0
x
x
f ; h ; . -h
а ; \ = ^; ^ = ^, (у=1,.,8); (£,/=1,..,4).
к ах к ах к х 2
1.7.1. Аппроксимация слагаемых - элементарной работы объемных и поверхностных сил
Если конечные элементы достаточно малы, то значения функций формы (см. формулу (2)) можно взять в ЦТ КЭ: /=0,25; (/=1,..,4). Тогда
ЗАи т Аg = ЗАц т Ф т Аg = ЗАц т1 ^ А^2 А^ А^2 ^ Аg2 Аg1 Аg2 ]т. (41)
Объемные силы Аg и g сразу задаются в ГСК, поэтому никаких преобразований производить не требуется. При этом для перехода к шагу нагружения номер «N+2» необходимо положить
^= §(ю +д§(^. (42)
Заслуживает рассмотрения задача получения аппроксимации ЗАи Ар и ЗАиТ р (см. уравнение
(33)), где Ар, р - приращение вектора поверхностных распределенных сил на шаге нагружения номер «N+1» вектор поверхностных распределенных сил, накопленных по результатам предыдущих шагов нагружения (рис. 2).
Вектор приращений распределенных сил (вектор Ар*), показанный на рисунке 2, задается в ЛСК, связанной с ориентацией нагруженной стороны в деформированном состоянии Следовательно, для подстановки в вариационное уравнение (1.32) он должен быть преобразован к декартовой ГСК по
формулам: {АЛ =^ - ^ - Ар- ^ (43)
1АР =АЛ - 12 + АР* - 11.
(N) у(п) _(Мт) (Ып) (N)\г(т)
Здесь / = —Х-; / = Х2 ~ х2
11 (Ы)£ ' 12 (Ы)т
тп тп
- направляющие косинусы, вычисляемые по координатам узлов КЭ в деформированном состоянии П^); (N) ^ - длина нагруженной стороны.
Рисунок 2 - Компоненты вектора распределенной нагрузки в ЛСК
1.7.2. Формирование векторов узловых реакций конечного элемента
Как следует из формулы (28), на шаге нагружения «N+1» необходимо располагать значениями вектора узловых реакций, накопленного по результатам предыдущих шагов нагружения: (-[^-Уо)-^) ^-Ур)-^) (Ы)Г(3)Л(3)
№Г(4)-1(4)Г).
Отметим, что его целесообразно вычислять сразу же в ГСК, чтобы избежать лишних преобразований. Рассмотрим, как формируется этот вектор на примере компоненты ¿Х1, совпадающей по направлению с осью ОХ1. Очевидно, он должен быть равен отношению суммы всех осевых сил, приложенных по предыдущим шагам нагружения (7=1Х), к текущей длине окружности,
1 N 1 N
а именно: , = - V 2л-(/ -1) г-(/ )А^ = - У (/-1)г-( 1 )А^ . (44)
2л-( N) г £ 1 (N)г 1=1 1
Тогда получаем окончательное выражение в виде:
( N )
N
r-txl (1-1)r•(1)Atxl
i=1
(45)
В качестве особенности формирования (№1)А1 - приращения вектора узловых реакций на шаге нагружения отметим, что составляющую (№1)А1 также целесообразно задавать в ГСК.
1.8. Общая структура МЖЭ и ВУС конечного элемента
Выполнив преобразования слагаемых вариационного уравнения (28), получим составляющие МЖЭ и ВУС конечного элемента, а именно:
K ^¿A*ETC* • A*E(N); K ^ SA cE
r
с
Ap ^ SA * ET • ß *AT •(N) rc 0SAUT Ag •(N) 0 SAq
( N ),
t T (N ),
(N)
L
' Ca
(N+1) Ap • S + S )■ (N+1) Ap • S +*1„) (N+1) Ap • S +Si„)
(46)
(47)
AP, ^
( N )
SAU T g •(N) Ф SAqT
L„
-SA * E
* _ t
r Ca
c E ^N) rr
p- (S1m +S1„) p- (S1m +S1„) p- (S1m +S1„).
•det ((N) J )•1 + SAqT
(N) (N) (N)
Г • t
'(1) 1 (1)
Г • t
'(2) 1(2)
Г • t
'(3) 1 (3)
(48)
Матрицы и векторы из формул (46) - (48) формируют уравнение равновесия изолированного КЭ в локальной СК при выполнении шага нагружения номер «N+1», а именно:
K • Aq - AP =
(N) r • (N+1)at
'(i) al (i)
(n)r •(N+1)at
'(2) al (2) (n)' •(N+1)at
'(3) at (3)
(49)
где
К = К + К ; АР = АР +АР - МЖЭ и ВУС конечного элемента.
£ } £ }
1.9. Вычисление модифицированного вектора приращений напряжений Кирхгофа, модифицированного вектора напряжений Кирхгофа и преобразование этого вектора к напряжениям
Коши-Эйлера
Перед выполнением шага нагружения номер «N+1», то есть в деформированном состоянии ОЦ мы располагаем значением вектора напряжений Коши-Эйлера (см. формулы (28)) вида:
аЕ =к о* оЕзГ • (50)
После того, как будет вычислено приращение вектора узловых степеней свободы по результатам шага нагружения номер «N+1», следует определить линеаризованную составляющую модифицированного вектора приращения деформаций Грина по формуле (29): АЕ = Ц • Аи = Ц ФАя = В • Ая,
а затем по формуле (19) - приращение модифицированного вектора напряжений Кирхгофа:
А*о = С* • А*Е -А*от •
Для перехода к следующему шагу нагружения (номер «N+2») необходимо преобразовать модифицированный вектор напряжений Кирхгофа стя+А*о к вектору напряжений Коши-Эйлера О+АО. Для этого необходимо выполнить преобразование (см. [1]):
1 дУ д¥,( ч (51)
a +AaE =
dfr,Y2,Y3) x dX,
ö(X1, X 2, X 3 )
2
2
где
д(—1 ЛЛ) _ D
(N+1)
э(х, x 2, x з) d
( N )
= det
д—; д—; д— д— д—; 0
дХ дХ2 дХ дХ, дХ2
д—2 дХ, д—2 дХ2 д—2 дХ = det д—2 дХ, д—2 дХ2 0
д—3 д—з д—3 0 0 д—3
дХ, дХ2 дХ3 дХ3
(52)
Производные
д— ЭХ,
(// из определителя (52) вычисляются по формулам (9).
Частную производную
дХ
будем вычислять как отношение радиусов вращения окружностей,
проходящих через ЦТ КЭ в деформированном состоянии П(№1) и соответственно. Тогда
дУ3 дХ3
Окончательно формулы (56) принимают вид: +А°п =
(N + 1)
r
(53)
(N )„
- D
(N)
D (N + 1) E
§; )+(дХ: )**+a-)+ 2 £ J—:&+*-..)
&2i2 + A^22 =
D
(N)
D
(N+1)
ЧдХ 1 У
'д—2л2
+ A*0"n)+ —^ (^2E2 +A*^22)+ 2
vdX 2 y
д- д-
2 ^ 2 ( E , л*. 22} 1 2 - 12
ЭХ дХ2
^2 + A*^12 )
°"12 +A^12
D
( N )
д-2
D
\ 2
( N +1)
■X
+A*a„)+I^m ^+a*^22)
чдХ 2
E А * 22 +Aa22 ' +
д-j д-2 д—! д-2 U E А.
чдХ, дХ2 дХ2 дХ
+A*^12 )
^ \ Л^-E D
iT33 + =
(N) f
33
D
( N+1)
д— väX3
j ^ +A*^33 ) .
(54)
Таким образом, определяются компоненты вектора напряжений Коши-Эйлера в деформированном состоянии 0(№1) для выполнения следующего шага нагружения.
2. Соотношения упругости анизотропного тела вращения для осесимметричной задачи
Рассмотрим тело вращения (рис. 3), имеющее слоистую структуру, и обладающего свойством осевой симметрии.
Тело вращения (см. рис. 3,а) образовано слоями композиционного материала с модулями упругости и коэффициентами Пуассона: Е, О,/ , И/, (/,/=1,2,3), уложенного с углами армирования ±у, как показано на рисунке 3,Ь.
Введем системы координат, связанные с геометрией тела вращения (см. рис. 3,Ь): СК 0'а1'а2'аэ' ориентирована в соответствии с направлением армирования слоя: 0'а1' - по волокну; 0'аз' - по нормали к поверхности слоя; 0'а2' - дополняет СК до правой тройки. Криволинейная СК 0'а^аз ориентирована в соответствии с геометрией тела вращения: 0'а1 - направлении меридиана; 0'аз - по нормали к поверхности слоя; 0'а2 - дополняет СК до правой тройки.
с
2
X
Рисунок 3 - Тело вращения: а - слоистая структура материала; Ь - геометрия укладки армирующих слоев
Тогда коэффициенты матрицы С'={сг/}, (/г/=1,..,6) и вектора ТКН Р' из соотношений упругости
Ао' = С' • Ас'-р' АT, где Ао' = [АопАо22\ Ао33АоАо31Ао Т;
А8' =[А£п\А£22 '' А£33'. АГ23'. АГзЛ АУп'Г ' вычисляют, как указано в работе [5]:
c11 =
1 - ^32 Е
c22
A A
с ' = G
c55 G13;
11;
E22; c23
c12
^21 + E A '
22 ;
c13
^31 + Ц E
A '
33;
^32 + У3К2 E
A
33
c33
1 E c G .
д E33; c44 _ g23 ;
c '_ G •
c66 g12 ;
c
cü'; A _ (1 - ^23 - 2l2l12 )(1 + ^23 );
(55)
(56)
(57)
(58)
"55 ^13' ^66 ^у ' ^ V 23 1 23/
Для компонент вектора ТКН, если слой представляет собой трансверсально изотропный материал, справедливо (см. [4])
Р11 = 011 «АХ ^ 012 «ТК ^ C13 «ТЯ; Р22 = C12 «АХ ^ 022 «ТЯ ^ C23 «ТЯ;
Р33 = C13 «АХ ^ C22 «ТК ^ C33 «ТЯ; Р23 = Р31 = Р12 = 0 • (59)
Здесь «ах - КЛТР слоя в направлении армирования; «гя - КЛТР слоя в трансверсальном направлении. Преобразование, переводящее тензор напряжений Тст '= {Аог>' } в тензор напряжений Тст = {Ао^},
вычисленный в СК О'ша^аз (см. рис. 4), определяется как [2]: Т = а ' Та 'T, (60)
о s 0"
где матрица преобразования « имеет вид:
а _
- ^c 0 0 0 1
(61)
Здесь с=соэ(у); Б=8т(у); у-угол армирования (см. рис. 4,Ь).
Следствием преобразования (60) являются формулы для вычисления коэффициентов упругости и температурных коэффициентов напряжений. Необходимо только учитывать, что слои - перекрестно армированные под углом ±у. С учетом этого получим:
^ 2 2 ^ 2 2 С11 = С11 0 + 2012 3 0 + 022 3
+ 4обб'30 ; 01 3 = 013 0 + 023 3 ;
c
4 + 2c12's 2c2 + c22'c4 + 4c66's lc2; c
2 2
22
"44
i 2 i 2
c44 c + c55 S ;
c33 c33 ;
c12 (c11 +c22 2(c12 +2c66 ))c S + c12 ;
c
23
? 2 ? 2 ? 2
c13 s + c23 c ; c55 _ c44 s
+ c55' c2 ; c66 _(c11 '+c22 '-2(c12 '+2c66 '))c2s2 + c
66 ;
ß11 _ cß11 }+S^ß22 ; ß22 _ sß11'+c ß22 ; ß33 _ ß33*; ß12 _ß23 _ ß31 _ 0 •
Затем необходимо преобразовать коэффициенты термоупругости от СК О'а^аз, связанной с геометрией поверхности тела вращения, к цилиндрической СК OYZX (рис. 5). Здесь л-ß - угол между соответствующими осями СК.
Рисунок 4 - Взаимное расположение СК, связанной с поверхностью тела вращения, и глобальной
цилиндрической СК
Они могут быть получены аналогично соотношениям (62). Если
äff 11 äff12 äff13 ÄffYY ÄffYZ ÄffYX
T„ = а^12 äff22 äff23 ; rp (YZX) _ а ÄffYZ ÄffZZ ÄffZX
_äff13 äff23 äff11 _ ÄffYX ÄffZX ÄffXX
(63)
- приращения тензоров напряжений в «старой» СК О'а^аз и «новой» СК 0У2Х, то матрица преобразования имеет вид:
а
а
а
а
Y " sin (ß) 0 cos(ß)" Y " s 0 c
a = Z 0 1 0 =Z 0 1 0
X - cos(ß) 0 sin (ß)_ X - c 0 s
, (//=1,..,6) вычисляются как (см. [5])
Тогда коэффициенты матрицы упругости C^YZX ^ = B
Бп = С S + 2(Сз + 2c55 )s2c2 + c33c4;
B13 = (с11 + С33 - 2(с13 + 2c55))sC + С13; B15 = (c33c' - С11^ + (с13 + 2С55)(^ - С'))sC ;
B22 = С22; ; B25 = (c23 - C12)sc ; B33 = cnc4 + C3S4 + 2(c3 + 2c55>V ; (65)
B35 = -SC(c11c — c33s — С13С — 2С55С + (c13 + 2С55)s ); B44 = C44S + С66С ;
: (c44 - c66 )sC ; B55 = (c11 + С33 - 2(c13 + 2С55))s'С' + С55 ; B66 = С44С' + C66S' '
B12 C12S + С23С :
B
а сама матрица коэффициентов упругости и вектор ТКН имеют структуру:
C
(YZX) '
SYY SZZ SXX Yzx Yyx Yyz
" B11 B12 B13 0 B15 0
ffZZ B12 B22 B23 0 B25 0
ffXX B13 B23 B33 0 B35 0
ffZX 0 0 0 B 44 0 B46
B15 B25 B35 0 B55 0
0 0 0 B 46 0 B66
-'(YZX) '
ßYY 'ßns2 +ß33C2"
ßzz ß 22
ßxx ß11c 2 + ß33s 2
ßzx 0
ßYX (ß33 -ßjsc
ßyz 0 _
(64)
(66)
Перейдем к векторам напряжений и деформаций, формируемых согласно формул (18) - (20). В глобальной системе координат 0ХХ2 для матрицы коэффициентов упругости и вектора температурных коэффициентов напряжений будем иметь:
а2 а3
C
(XYZ)
£XX SYY Yxy £ZZ
^xx B33 B13 B35 B23
= ^ YY B13 Bn B15 B12
°XY B35 B\5 B55 B25
°ZZ B23 B12 B25 B22
ß
(XYZ)
ßxx Pyy
ßYX
ßzz
ßnc2 +ß33s
ßnS2 +ß33c
(ß33 -ßn)
sc
ß.
22
(67)
Процедура вычисления угла Р (см. рис. 5) проиллюстрирована на рисунке 5. Очевидно, что этот угол может быть вычислен как
р = (у. - у )Дх - )).
(68)
О
X
-1
Рисунок 5 - К определению Р - угла между осями систем координат
3. Тестирование
На алгоритмическом языке PASCAL была разработана оригинальная программа, реализующая МКЭ. В статье приводятся результаты двух тестовых примеров из числа решенных при отработке программы. 2.1. Определение критического внешнего давления сферической оболочки
Рассматривается задача определения критического внешнего давления сферической оболочки. Расчетная схема задачи представлена на рисунке 6.
Приняты следующие геометрические размеры оболочки и механические характеристики материала: ^=100 мм - радиус серединной поверхности оболочки; й=0,5-мм - толщина; E=20 МПа - модуль упругости материала; v=0,3 - коэффициент Пуассона.
Оболочка разбивалась на 16000 квазидвумерных КЭ с четырехугольным сечением, 8 конечных элементов по толщине оболочки.
Если внешняя нагрузка считается приложенной к серединной поверхности оболочки, то решение может быть получено аналитически [7], а именно:
h2
Рр = 1,21E— = 1,21 • 20 -IQ6
R
25-10-0,12
= 605 Па.
(69)
Рисунок 6 - Сферическая оболочка под внешним давлением: a - расчетная схема задачи; Ь - форма потери устойчивости
2
2
8
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «ИННОВАЦИОННАЯ НАУКА» №6/2016 ISSN 2410-6070_
Алгоритм определения критического давления при решении задачи статики в геометрически нелинейной постановке сводится к следующему: реализуется процесс шагового нагружения оболочки внешним давлением, и отслеживается приращение максимального отклонения деформированного контура (то есть в нашем случае нормального прогиба) на шаге нагружения:
Max(AU2 + AV2 )2 . (70)
Если величина из формулы (70) на каком-либо шаге нагружения становится много большей, чем на предыдущих шагах, это означает, что оболочка потеряла устойчивость.
В первом случае при решении задачи МКЭ давление прикладывали к серединной поверхности оболочки. Критическое давление, определенное численно, составляет 613 Па, что хорошо соотносится с аналитическим решением (см. формулу (69)). Осесимметричная форма потери устойчивости представлена на рисунке 6,b, а зависимость «максимальный прогиб - нагрузка» - на рисунке 7.
Рисунок 7 - Значения максимального прогиба оболочки в зависимости от величины нагрузки, приложенной к серединной поверхности
2.2. Расчет ортотропной цилиндрической оболочки под действием
внутреннего давления
На рисунке 8 показана короткая ортотропная цилиндрическая оболочка, жестко заделанная по торцам и нагруженная внутренним давлением величины _ро=108 Н/м2 (см. [3]).
Рисунок 8 - Ортотропная цилиндрическая оболочка под внутренним давлением
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «ИННОВАЦИОННАЯ НАУКА» №6/2016 ISSN 2410-6070_
Заданы следующие геометрические размеры и механические характеристики материала оболочки: L=0,2 м - длина; R=0,2 м - радиус; h=10-2 м - толщина; £¡=75 ГПа; £2=£3=20 ГПа; Gi2=12,5 ГПа; G23=G3i=6,3 ГПа - модули упругости материала; Vi2= V23=0,25; V3i=0,067 - коэффициенты Пуассона. Оболочка разбивалась на 5760 четырехугольных КЭ: 32 КЭ по толщине оболочки и 180 КЭ по ее длине.
На рисунке 9 представлена зависимость величины максимального прогиба в оболочке от значения приложенного внутреннего давления. Произведено сравнение с аналогичными результатами из работы [3]. Можно отметить их хорошее соответствие: максимум прогиба, рассчитанный по оригинальной программе, составляет 13,78-10-2 м.
Анализ эпюр напряжений, представленных на рисунке 10 a,b, позволяет сделать вывод о непротиворечивом характере их распределения: меридиональные напряжения имеют максимум в центре свободной поверхности оболочки, если не считать участков закрепления, (см. рис. 10,a); для нормальных напряжений (см. рис. 10,b) выполняются силовые граничные условия.
(а) (Ь)
Рисунок 9 - Зависимость величины максимального прогиба в оболочке от значения приложенной нагрузки:
а - приведенная в работе [3] (перемещения - в «см»); Ь - рассчитанная по оригинальной программе (перемещения - в «м»)
(а) 0>)
Рисунок 10 - Напряжения Коши-Эйлера (Па) в оболочке: a - меридиональные; b - нормальные
4. Заключение
В статье приведен алгоритм построения четырехугольного изопараметрического КЭ тела вращения с 8-ю степенями свободы с использованием модифицированного подхода Лагранжа: указаны функции для аппроксимации основных неизвестных, вариационное уравнение равновесия изолированного КЭ, а также аппроксимации входящих в него слагаемых. Выведены формулы преобразования модифицированного вектора напряжений Кирхгофа к вектору напряжений Коши-Эйлера применительно к данной задаче.
Подробно рассмотрен вопрос определения термоупругих характеристик слоев композиционного материала, составляющих тело вращения. При получении этих зависимостей авторы использовали материалы, представленные в работе [5].
Предложенная методика построения конечного элемента тела вращения может быть использована для построения КЭ другого вида: как более простых изопараметрических треугольных, так и элементов более высокого порядка - субпараметрических.
Корректность разработанного алгоритма подтверждается решением нескольких тестовых примеров. Список использованной литературы:
1. Васидзу К. Вариационные методы в теории упругости и пластичности. - М.: Мир, 1987. - 542 с.
2. Гафинов В.И. Основы теории упругости и пластичности. Часть I. - М.: МВТУ, 1977. - 58 с.
3. Голованов А.И., Тюленева О.Н., Шигабутдинов А.Ф. Метод конечных элементов в статике и динамике тонкостенных конструкций. - М.: Физматлит, 2006. - 392 с.
4. Кристенсен Р. Введение в механику композитов. - М.: Мир, 1982. - 334 с.
5. Носатенко П.Я. Исследование геометрически нелинейного напряженно-деформированного состояния анизотропных оболочек вращения методом конечных элементов. Диссертация на соискание ученой степени к.ф-м.н. М.: МАМИ, 1984. - 164 с.
6. Секулович М. Метод конечных элементов. - М.: Стройиздат, 1993. - 664 с.
7. Биргер И.А., Шорр Б.Ф., Иосилевич Г.Б. Расчет на прочность деталей машин. - М.: Машиностроение, 1993. - 640 с.
8. Bathe K.J. Finite Element Procedures in Engineering Analysis. - Prentice? Hall, Inc., 1982. - 735 p.
© Баслык К.П., Зарипов Д.Х., Чурилин С.А., 2016
УДК 62-112.83
Н.С. Васильев
Студент 6 курса, кафедра ПТМ С.Л. Заярный к.т.н., доцент, кафедра ПТМ КФ МГТУ им. Баумана Г. Калуга, Российская Федерация
МОДЕЛЬ ФРИКЦИОННОГО БОЛТОВОГО СОЕДИНЕНИЯ, С УЧЕТОМ ПОДАТЛИВОСТИ СТЫКА
Аннотация
В статье приводится расчетная модель фрикционного болтового соединения, нагруженного сдвигающими силами, учитывающая податливость стыка и материала сопрягаемых пластин.
Ключевые слова
Фрикционное болтовое соединение, контактный слой, контактное взаимодействие, податливость стыка.
Получение расчетных зависимостей и моделей работы фрикционного болтового соединения, в его до сдвиговой и сдвиговой стадиях, представляет интерес для оценки его повреждаемости [1-3, 5].
Болтовое фрикционное соединение (рис. 1, а) имеет двухмерную циклическую структуру (рис. 1, б), болтовым суперэлементом (БСЭ) которой, с характерным размером, определяемым площадь опорного конуса, является отдельное болтовое скрепление (рис. 1, в). С точки зрения обеспечения межэлементных взаимодействий БСЭ является сложной функциональной структурой (рис. 2 а,б), которая может быть аппроксимирована трехэлементной схемой (рис. 2 в) в которой учитываются обобщенные взаимодействия: силовое замыкание элементов распределенной нагрузкой ; ретрансляция силового потока (элементы 1, 2);
