CC BY
15
-о- Hb Hb
Численные методы расчета конструкций
РАСЧЕТ ПРОИЗВОЛЬНЫХ ОБОЛОЧЕК НА ОСНОВЕ МКЭ В СМЕШАННОЙ ФОРМУЛИРОВКЕ С ИСПОЛЬЗОВАНИЕМ АППРОКСИМАЦИИ ТЕНЗОРНЫХ ПОЛЕЙ
Н.А. ГУРЕЕВА, канд. техн. наук, ст.преподаватель
Ю.В. КЛОЧКОВ, д-р техн. наук, профессор
А.П. НИКОЛАЕВ, д-р техн. наук, профессор
Волгоградская государственная сельскохозяйственная академия
На основе функционала Рейсснера с использованием соотношений теории упругости без упрощающих гипотез разработан алгоритм формирования смешанной матрицы связи между усилиями и перемещениями для шестигранного конечного элемента с узловыми неизвестными в виде перемещений и напряжений. Для выражения перемещений и напряжений через узловые неизвестные использовался способ аппроксимации векторных и тензорных полей. Показана возможность использования разработанного конечного элемента в расчете тонких оболочек.
КЛЮЧЕВЫЕ СЛОВА: МКЭ, шестигранный КЭ, функционал Рейсснера
Обычно при расчете оболочек на основе МКЭ используется теория поверхностей с гипотезами Кирхгофа или Тимошенко о деформировании нормали к срединной поверхности оболочки. В настоящей работе изложен алгоритм расчета произвольной оболочки на основе МКЭ в смешанной формулировке. Разработан объемный конечный элемент в форме шестигранника с узловыми неиз-
вестными в виде компонент вектора перемещении и компонент тензора напряжений. Для выражения компонент вектора перемещения внутренней точки конечного элемента и компонент ее тензора напряжений через узловые неизвестные использовался способ интерполяции векторных и тензорных полей трилинейными соотношениями, что позволило в неявном виде учесть смещение конечного элемента как жесткого целого.
На примере показана эффективность предложенного способа аппроксимации искомых величин как векторных и тензорных полей по сравнению с традиционным способом аппроксимации искомых величин как скалярных полей.
1. Геометрия оболочки. Положение произвольной точки М отсчетной поверхности определяется радиусом - вектором
Я = хк [в1,в2 Ук, (1.1)
где ва - криволинейные координаты; X - декартовы координаты точки; ¡к - орты декартовой системы координат. Здесь и ниже греческие индексы принимают значения 1, 2, а латинские индексы - 1, 2, 3.
Векторы, касательные к отсчетной поверхности, определяются дифференцированием (1.1) по криволинейным координатам ва
аа = Х,аК . (1.2)
Орт нормали к отсчетной поверхности оболочки определится векторным произведением
- а х а
а = -г——-^г. (1.3)
ах а2
Базис векторов {а} = {ааа} в точке М можно представить в матричной форме {а}=[М ]{-}, (1.4)
где {- }Т = {-/2/3} - орты декартовой системы координат; [М] - матрица коэффициентов. Дифференцированием (1.4) определяются производные базисных векторов произвольной точки отсчетной поверхности оболочки
{а,Л=К]{1 (1.5)
Используя (1.4), производные (1.5) можно представить разложением по векторам базиса рассматриваемой точки
{а,а} = [х,а][м]' {а} = [ша ]{а}. (1.6)
2. Перемещения и деформации. Положение точки М, отстоящей на расстоянии Г от отсчетной поверхности, определяется радиусом-вектором
Я' = Я + Га . (2.1)
Базисные векторы точки М определяются дифференцированием (2.1)
ёа = аа + Га,а ; ёз = Я = а . (2.2)
Вектор перемещения произвольной точки М выражается в локальном базисе
V = Vрар + уа = {а}Т {V}, (2.3)
где {у}Т = {у^у2у} - вектор-строка компонент вектора перемещений точки М.
Производные вектора перемещения по координатам а , Г с учетом (1.6) определяются выражениями
V ,а=fарар + аа-Уп=у,р ар + у„а, (2.4)
где /а,/а - функции компонент вектора перемещения и их производных.
Положение точки М ' в деформированном состоянии оболочки определяется радиусом-вектором
Л'* = Я + 'а + V = Л* + У . (2.5)
Базисные векторы М' в деформированном состоянии оболочки определяются дифференцированием (2.5)
£ = Яа = Яа +У,„ = Еа + /Р5р + /„5 ; £ = Д Г = 3 + V,* . (2.6)
Деформации в произвольной точке определяются как разности компонент метрических тензоров деформированного и исходного состояний [1]
*Ц =(?* - §у )/2 ^У7, у + ё У, 1+К у} )/2. (2.7)
Используя (2.2) и (2.4), линейную часть соотношений (2.7) можно представить в матричном виде {е}= [ь]^}, (2.8)
6x1 6x3 3x1
где {е}г = {е11е22е 33 2е122е132е23 }- вектор-строка деформаций в точке М; [Ь] -матрица дифференциальных операций.
3. Закон Гука и функционал Рейсснера. Соотношения между деформациями и напряжениями в произвольной системе координат между компонентами тензоров деформаций и напряжений имеют вид [1]
=-^1 +[ Е+Е ,
И=№},
(3.1)
(3.2)
или в матричной форме
{М 1 11 22 23 12 13 23 ( « «
ст} = р ст ст ст ст ст вектор-строка напряжений произвольной точки оболочки; J1 (ст) = ё тпсттп - первый инвариант тензора напряжений; Е - модуль упругости материала; ц - коэффициент Пуассона.
Функционал Рейсснера для оболочки, загруженной поверхностной нагрузкой, имеет вид [2]
Пя = { {ст}Т [Ь]{М -1{ {ст}Т [в]{а^У -1 Мт {^ ,
2 1
|{у}ОУ - -] {ст} [и\{ст}ау - ] {V} {д}ая, (3.3)
У 2 У я
где У - объем оболочки; - поверхность оболочки с заданными внешними силами; {д}Т = {д'д2 д3}- вектор поверхностных нагрузок.
4. Смешанная матрица деформирования конечного элемента. Объемный конечный элемент выбран в форме шестигранника произвольного слоя оболочки с узлами 1, у, k, I нижней грани и узлами т, п, р, h верхней грани (рис. 1, а). Для выполнения численного интегрирования по объему шестигранника глобальные координаты 6р и ' аппроксимировались трилинейными соотношениями через локальные координаты куба (рис.1, б), изменяющиеся в пределах
-1 < < 1,
А = {/ (£,^)}Т {Лу }, (4.1)
где
{Яу } ={ХА]^2!АтАпАрАк}
- вектор-строка узловых значений величины X. Под символом X понимаются
координаты вр,' . Производные глобальных коор-
У
Рис. 1
динат в локальной системе 0," ,0,^,0,",, 1,ц, и локальные координаты в
глобальной системе определяются диффе-
ренцированием (4.1). Векторные поля перемещений и тензорные поля напряжений внутренней точки конечного элемента аппроксимируются через узловые неизвестные трилинейными соотношениями
V = {/ )}Т V ¡1 С = {/ (е,П,С )}Т Су }, (4.2)
где V } = } - строка векторов перемещений узловых точек;
(—V
1 1 _ ___ __о о _ __1 о _ __1 о _ __ЛЛ _ __I I /* Ч
а = а а1а1 + а а2а2 + а аа + 2а а1а2 + 2а а1а + 2а а2а = <а> {а} -тензор напряжений внутренней точки конечного элемента; а =а а1 а +С а2а2 +а а а + 2а а, а2 + 2а а а + 2а а2а = < а > {а } - тензор напряжений узловой точки конечного элемента (" = 7, к, I, т, п, р, И);
I ~ |Т 7 ~ I ~к ~I ~т ~п ~ р ~И \ «
Су] = аа а а а а аа ^-строка тензоров напряжений в узловых точках.
Базисные векторы узловых точек {а"} = ,а"} с учетом (1.4) выражаются через базисные векторы произвольной внутренней точки {а}т ={а1,а2,а} конечного элемента
{а"} = [М"][М]_1 {а} = [^]{а} , (" = 7,1, к, I, т, п,р, И). (4.3)
3x1 3x3 3x1
С использованием (4.3) диадные произведения базисных векторов узловой точки (а"}можно выразить через диадные произведения базисных векторов произвольной точки конечного элемента (а} в матричной форме
{а-} = [а"]{а} , (м> = 7,1, к, I, т, п, р, И). (4.4)
6x1 6x6 6x1
Столбец узловых векторов перемещений может быть представлен матричным соотношением
{V }=[Т" ]Ь}
8x24 24x1
(4.5)
где элементами матрицы [Т " ] являются базисные векторы узловых точек, а
{"у }Т = "..Vй }.
1x24
Векторы перемещения внутренней точки конечного элемента и его производные, учитывая (4.5) и (4.3), можно представить матричными выражениями
= {/ ^Л {Vy }= {а}
I ^х8 Л 8x1 1x3
V,а = {/
^1x8
где
8x1
Т
/1 [N ]Т /2 [N ]Т ... / [N]
3x3
{"у }; (4.6)
/1 [N] /1 [& ] ... /,1 [^] {"у},
{¥ у }= {а}Т
8x1 1x8
3x3
[Т
8x1
Г
Г
3x3
3x3
с*
Принимая во внимание (2.3) и (2.4), из (4.6) можно получить аппроксимирующие матричные выражения для компонент вектора перемещения и их производных
\Т ( 1 р (г ЛТ ( ^ (-, лТ
{V} = [}, V, р = {'„}Т {vy};V,; = {/„} {vy}; V, „ = {Ъа}Т {vy} . (4.7)
3x1 3x24 24x1
Столбец узловых тензоров напряжений {г ; }= {~ ' ст }ст 1сст ' ст т ст пст р ст к }
можно записать в матричном виде
8x1
{о; }= н; { ст;},
8x48
!стУ
48x1
(4.8)
где {ст;} = |{ст''} {ст {ст•}' {ст'} {стт} {стп} {стр}'{ст'} ,
1x48
а элементами матрицы [Н ] являются диадные произведения базисных векторов узловых точек. Используя (4.8) и (4.4), аппроксимацию тензора напряжений (4.2) можно представить выражением
ст =
{5}Т /1 ]Т/2 ]Т .../ (£,Ч,С)[0"]Т { ст;} . (4.9)
6x6
1x6
Принимая во внимание диадное представление тензора о = {а} { ст}, ком' 6x1 1x6
поненты тензора напряжений внутренней точки элемента на основании (4.9) можно записать произведением матриц
{ст}=[5 ]{ст у }. (4.10)
6x1 6x48 48x1
При использовании аппроксимирующих выражений (4.7), (4.10) функционал (3.3) принимает вид
ПК = {сту}Т |ИТ [5]ауК¡-2{сту}Т {ИТ ШаУ{сту}-{vy}T |Д4 {д}ая, (4.11)
1x48 У48x6 6x24 24x1 2 1x48 У48x6 6x66x48 48x1 1><24 я 24x3 3x1
где [в] = [Ь][А\ .
6x24 6x3 3x24
Минимизируя функционал (4.12) по узловым неизвестным {ст у }, {уу } и приравнивая производные нулю, можно получить систему уравнений П ( ) г„1 ( ) „ дП„
Я - И К }- [Н ]{рг }= 0. [о]Т К}- / } = 0
48x24 24^ 48x48
у ) 48x24 24x1 Т
48x1
у J
24x48 48x1 24x1
(4.12)
где О] = {МТ {в}аУ; [н ] = { |>}Т {ДО} аУ, {/ } = | Г{ А} {д} а5.
48x24 у 48x6 6x24 48x48 У 48x6 6x6 6x48 24x1 ^ 24x3 3к1
Систему уравнений (4.12) можно представить в традиционной для МКЭ форме
^}=/ }, (4.13)
72x72
где
[с ] =
-[н] О]
48x48 48x24
[о]Т [о]
24x48 24x24
72x1
72x1
- матрица деформирования конечного элемента;
1x6
1x6
6x6
6x6
Т
} = {{ау }Т {"у }Т } - вектор узловых неизвестных конечного элемента;
1x72 [ 1x48 1x24 ]
{/у }={ {о}Т {/ }Т } - вектор узловых усилий конечного элемента.
72x1
О
^1x48 1x24
Матрица деформирования всей конструкции формируется с использованием традиционной процедуры МКЭ [3].
5. Оболочка вращения. При расчете оболочек вращения радиуса г = г(х) в
качестве криволинейных координат 61,02 принимаются осевая координата х и окружная координата 0 . Радиус-вектор произвольной точки отсчетной поверхности определяется выражением
Я = хг + rsin6 + гcos0k . (5.1)
Дальнейший расчет оболочки вращения выполняется по алгоритму, изложенному в разд. 1-4.
Пример. Определено напряженное состояние усеченного параболоида вращения, внутренняя образующая которого описывается функцией
2 =
I- а 2 —:—х + а . 12
Были приняты следующие исходные данные (рис. 2): а = 0,2м, d = 0,02м, I = 0,5м, И = 0,02м, Е = 2-105 МПа, л = 0,3, q = 10 МПа.В узловых точках элементов на левом краю оболочки установлены пружинные опоры (рис. 2), при изменении жесткости которых узловые точки на левом краю получают перемещения
в см. Достаточная сходимость вычислительного процесса получена при дискретизации конструкции на 288 элементов с двумя элементами по толщине оболочки. В табл. 1 приведены численные значения меридиональных напряжений <ти и кольцевых напряжений а22 в точках 1-4. Как видно, в точках 3 и 4 меридиональные напряжения равны нулю, а на правом краю среднее меридиональное напряжение равно 465,93 МПа.
Среднее напряжение, подсчитанное из условия равновесия равно
я(а2 -С12)
аф = q
2п(а + И/2)И
= 47,1 МПа.
Рис. 2
Как видно расхождение в результатах не превышает 1,5%.
Таблица 1
№ точки Напряжения МПа
= 0,0000 = 0,0393 = 0,1969
а11 а22 а11 а22 а11 а22
2 3 4 5 6 7
1 45,370 89,756 49,313 93,249 49,446 93,363
2 47,804 82,258 43,872 83,127 43,763 83,156
3 0,037 25,434 0,038 25,434 0,038 25,434
4 0,004 5,420 0,007 5,421 0,007 5,421
= 19,7047 = 197,205 = 1980,13
а11 а22 а11 а11 а22 а11
8 9 10 8 9 10
1 49,481 93,392 49,483 49,481 93,392 49,483
2 43,732 83,164 43,730 43,732 83,164 43,730
3 0,038 25,434 0,038 0,038 25,434 0,038
4 0,007 5,421 0,007 0,007 5,421 0,007
Анализ результатов показывает, что численные значения напряжений при значительных смещениях как жесткого целого V1 остаются практически неизменными. Расчеты при использовании скалярной аппроксимации компонент вектора перемещения и компонент тензора напряжений показали, что численные значения напряжений получились такими же, как в колонках 2,3 табл. 1. только при наличии жестких опор на левом краю оболочки. Если оболочка получала даже весьма малое смещение V1, результаты расчета оказывались просто неприемлемыми. Таким образом, предложенная аппроксимация векторных и тензорных полей искомых величин внутренней точки конечного элемента позволяет учет смещения конструкции как жесткого целого.
Конечный элемент на основе разработанного способа аппроксимации векторов перемещений и тензоров напряжений может найти применение в программных комплексах для расчета прочности оболочек произвольной толщины в трехмерной постановке.
Л и т е р а т у р а
1. СедовЛ.И. Механика сплошной среды, т.1. - М.: Наука, 1976. - 536 с.
2. Галагер Р. Метод конечных элементов. Основы. - М.: Мир, 1984. - 428 с.
3. Постное В.А., Хархурим И.Я. Метод конечных элементов в расчетах судовых конструкций. - Л.: Судостроение, 1974. - 344 с.
' ' нь -о- -о-
ЧИСЛЕННЫЙ АНАЛИЗ УСТОЙЧИВОСТИ РАВНОВЕСИЯ ПРОСТРАНСТВЕННЫХ ФЕРМ В ГЕОМЕТРИЧЕСКИ НЕЛИНЕЙНОЙ ПОСТАНОВКЕ
В.В. ГАЛИТТТНИКОВА, канд. техн. наук, доцент Российский университет дружбы народов, Москва
В работе излагаются алгоритмы вычисления критических состояний пространственной фермы и продолжения решения в сингулярных точках, реализующие методики, разработанные автором. Приводится пример расчета сетчатой оболочки на устойчивость равновесия по разработанной программе.
КЛЮЧЕВЫЕ СЛОВА: устойчивость, пространственная ферма, сетчатая оболочка
Вычисление точек бифуркации гибких, подверженных большим деформациям конструкций, а также продолжений траекторий их нагружения после прохождения этих точек остается одной из сложнейших задач вычислительной механики. Очевидно, что задача устойчивости равновесия таких конструкций может быть достоверно решена лишь при учете предшествующей истории нагру-жения, то есть в геометрически нелинейной постановке. Однако, несмотря на то, что численные методы геометрически нелинейного анализа хорошо развиты, в окрестностях точек бифуркации возникают вычислительные трудности, препятствующие созданию надежно работающих алгоритмов продолжения решения в этих точках. Формулировка и проверка достоверности методов геометрически нелинейного анализа конструкций существенно затрудняется недостатком точных аналитических решений задач для конструкций, подвергающихся значительным деформациям перед потерей устойчивости. Достоверность и точность методик, предлагаемых в данной работе, оценивалась на тестовых примерах, для которых автором получены аналитические решения [1].
Для выполнения геометрически нелинейного расчета пространственных ферм на деформации по методу конечных элементов автором разработана методика, основанная на модификации метода постоянных дуг (Arc Length Method), в которой используется процедура итерационного уточнения секущей матрицы жесткости конструкции [2]. В качестве характеристической кривой решения используется нормализованная траектория нагружения. Алгоритмы, описываемые в настоящей работе, являются частью обобщенного метода геометрически нелинейного расчета пространственных ферм, включающего расчет на деформации и потерю устойчивости.
