CC BY
59
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
Том 151, кн. 3
Физико-математические пауки
2009
УДК 539.3
ТРЕХМЕРНЫЙ КОНЕЧНЫЙ ЭЛЕМЕНТ ДЛЯ РАСЧЕТА ТОНКОСТЕННЫХ КОНСТРУКЦИЙ
А.И. Голованов, М.К. Сагдатуллип
Аннотация
Работа посвящена построению нового конечного элемента для расчета средней толщины оболочки па основе модификации трехмерного изопараметрического восьмиузлового элемента путем введения гипотезы малости напряжений обжатия и использования техники понижения порядка аппроксимации, при этом использован метод двойной аппроксимации по точкам суперсходимости. На численных примерах показана эффективность данного подхода.
Ключевые слова: оболочечпый конечный элемент, упругие деформации, метрический тензор, метод двойной аппроксимации, гипотеза малости напряжений обжатия.
Введение
Цель настоящей работы разработка такой методики модификации трехмерного изопараметрического восьмиузлового конечного элемента (КЭ) сплошной среды, чтобы было возможно моделирование оболочек средней толщины при однослойной аппроксимации по толщине.
Обычно подобные элементы строятся с использованием степеней свободы, определенных на срединной поверхности и включающих углы поворота нормального волокна. Обзор таких КЭ представлен в работах [1 3]. Как правило. КЭ с угловыми степенями свободы демонстрируют хорошую точность и весьма эффективны в расчетах оболочек малой и средней толщин. Однако их использование весьма затруднительно при моделировании сопряжений оболочек с массивными телами, так как необходимо выражать узловые перемещения трехмерных элементов через углы поворота оболочечных КЭ. Поэтому получили распространение специальные элементы оболочек, которые имеют в качестве узловых степеней свободы перемещения узлов, расположенных на лицевых поверхностях. Фактически речь идет о модели оболочки средней толщины с учетом обжатия. Примеры подобных КЭ и их использования описаны в работах [3 12]. Отметим, что идейно близкие построения предложены в работах [4. 5. 8]. Описанная методика частично использовалась в работах [6. 7].
Основные этапы построения предлагаемого КЭ состоят в следующем. Фрагмент оболочки представляется как трехмерное тело. Вводятся изопараметрические полилинейные аппроксимации геометрии и неизвестных компонент перемещений. Узловые точки располагаются в вершинах выделенного параллелепипеда. Использование аппроксимации перемещений по толщине эквивалентно введению кинематической гипотезы об их линейном распределении по толщине с учетом обжатия оболочки.
Вариация потенциальной энергии деформации (работа внутренних напряжений на возможных деформациях) вычисляется численно по квадратурной формуле Гаусса Лежандра (два узла по каждому координатному направлению всего
Рис. 1
2 х 2 х 2 = 8 квадратурных узлов на КЭ). В каждом квадратурном узле вводится статическая гипотеза, обычно используемая для оболочек средней толщины, об отсутствии влияния поперечного напряжения на мембранные, изгибные и деформации поперечного сдвига. Фактически задача поперечного обжатия отделяется от классической задачи теории оболочек с учетом поперечного сдвига.
Главный источник ошибок, возникающих при использовании подобных КЭ. которые делают их практически непригодными. появление «ложных» деформаций поперечного сдвига. В настоящей работе предполагается использование метода двойной аппроксимации по точкам суперсходимости [1. 2].
1. Построение трехмерного КЭ оболочки средней толщины
Рассмотрим восьмиузловой трехмерный КЭ. изображенный на рис. 1. Криволинейная Лагранжева система координат £1, £2, £з становится прямолинейной, но в общем случае неортогональной.
Введем в рассмотрение радиус-вектор произвольной материальной точки в виде
К = X®(£1, £2, £3) е
Для проекций X1 соответствующие аппроксимации запишем в виде
X i(£1,£2,£3)^ X Ní(£1j£2,£3),
¿1 ¿2 ¿3\ =
4=1
где полилинейные функции формы определяются по формулам
ад, £2, £з) = 8(! + £1£1) (! + £?£2) (1 + £3£3),
£(, £2, £3 _ локальные координаты узлов, такие что ££ = ±1 в зависимости от номера узла (см. рис. 1). Соотношение
дХ4
определяет ковариантный базис
= = Rke¿
R= ЕXdN = 8Еxi£k (! + £t£') (1 + £Г£т):
t t
где номера к, 1, т удовлетворяют следующим условиям:
к, 1, т =1, 2, 3; к = 1 = т = к.
к
к д£к 1 к к • дХ4 2л/С
где € ктп — символы Леви - Чевита. Значение л/О = И1 • [112 х Из] задает величину элементарного объема ¿V. Компоненты метрического тензора в различных базисах вычисляются по формулам
_^ дХт дХт _л
= = Е ~о£Т = Е йтдт,
тт
С3 = 114113 = V ^^ = V Д4'тД3'т, ^ дХт дХт ^
тт
/о _ \ л г-чтптуг г,
= О ПтПп-
Вектор перемещений
т
т,п
3т
и = ит (£1,£2,£3)
аппроксимируется следующим образом:
и4 (£1,£2,£3) = Е и; N (£1,£2,£3).
Г = 1
Градиенты компонент перемещений представлены в виде
и- = е V ^^ иг =&т иг д£4 .
; = 1
Ковариантиые компоненты тензора бесконечно малых деформаций определяются следующим образом:
1 „ „ тт ч 1 (дит дХт дХт дит
е, = - (И; • К,- + К. • и,- ) = - > -:--- +--:---
43 2К 3 3 2 т=1 V д£4 д£з д£4 д£з
= 1 ЕЕ(итХт + иГХЛЦ. а)
т=1 г,4=1 £ £
Кратко опишем процесс построения матрицы жесткости. Примем в качестве базового вариационное уравнение
J а3 ¿V = ^ Г • Ш ¿V + У <;п • Ш
где V - физический объем, ограниченный поверхностью 5 = 5й и 5" П = 0,
п
части поверхности на которой определены силовые граничные условия. Кинематические условия на части Би поверхности 5 выполняются за счет специальным образом определенных аппроксимаций.
Введем систему квадратурных узлов, которые в локальной системе координат обозначим через , Щ, ^ и соответствующий весовой множитель . Если предполагается использование упомянутой выше восьми точечной формулы Гаусса
Лежандра, то / =1, 2,..., 8; ^ = ±^3' = 1
Далее интеграл в левой части уравнения виртуальных работ заменяется конечной суммой, то есть
J стг]5ег] вУ = JJJ стг]¿е^-л/С^1^2^3 = ^ атп5етп^/С,
V /
где С = 0.
Если ввести обозначение
Ег,т = 1Л хт (дN дN + дN \ т
г] 2 = 1 V д¥ дС + д¥ ) ' 1 ;
то компоненты деформации можно записать в виде
3
=
ЕЕитЕГ. (з)
е.
г= 1 т= 1
Используя закон Гука и учитывая, что в криволинейном базисе роль единичного тензора играет метрический тензор, имеем:
стг] = 2Мег] + АСг] /1,, (4)
где первый инвариант тензора деформации имеет вид:
/ _ стп е
11е С етп •
Отсюда получаем классическое уравнение
стг] = (2^СгтС]п + АСг]Стп) еШп.
Рассмотрим подынтегральное выражение левой части уравнения виртуальных
()
стг]5ег] = {2^0гт0]п + АСг]Стп) етп6ег], (5)
Используя приведенное выше представление (5). получим
3 8
к ХТТ1 т?г,т 173,1
Ш = стг] ¿-г] = £ ]Г икк ] Щ
к,= 1 г,з=1
3 8
х (2^агта]п + АСг]атп^у/а = Е ик °гы дц1з,
к,1= 1 г ,з=1
где = (2^0гт0]п + АСг]Стп) Е^Е^и^О.
Суммируя значения величин , вычисленных в системе квадратурных точек £4, получаем блоки соответствующей матрицы жесткости трехмерного восьмиуз-лового КЭ линейной задачи теории упругости.
Рассмотрим способ учета малости напряжений обжатия. Запишем эту гипотезу в виде равенства
<33 = 0. (6)
Подставим в уравнение (6) выражения компонент тензора напряжений а43 = = С4к акт в виде (4) и получим
<33 = 2^ £33 + АС33/1е =
= 2М£33 + АС33 (С43£4, + 2С43£43 + С33£33) = 0, (7) где г, j = 1, 2. Выразим го уравнения (7) деформацию обжатия £33:
£33 = - 2^А§33О3 (О43£43 +2О43£43) . (8)
Теперь преобразуем соотношения упругости для ковариантных компонент напряжений отдельно для напряжений а^, г, j = 1, 2, характеризующих мембранное и изгибное напряженное состояния, и а43, описывающих напряжения поперечного сдвига.
Из (4) с учетом (8) получаем:
73 = 2м£®3 + АС^ (Стп£тп + 2Ст3£т3 + С33£3^ =
= + АС^ (^*СтП£тп + 2^*Ст3£т^ ,
<43 = 2^£43 + АС43 (^*СтП£тп + 2^*Ст3£т^ , гдеМ = 1 - т,п = 1,2
Так как в рассматриваемом КЭ фигурируют узловые степени свободы на лицевых поверхностях, то фактически в кинематике учитываются возможные обжатия оболочки (изменение ее толщины). Следовательно, в расчетную схему необходимо ввести упрощенный закон упругости, связывающий напряжение обжатия 733 с
£33
<33 = ^3£33,
где Е3 - модуль жесткости на обжатие. В частности, можно принять, что
Е М (3А + 2М)
Е3 = —г-:-.
А + м
Запишем в общем виде закон Гука:
3
^ ^ £mn, j 1, 3
т,п= 1
где компоненты тензора упругости получены на основе приведенных выше рассуждений:
= 2^а7 ^ + АСав
^аз7й = Дза7й = А Саз^*С7<5, ^а373 = Бза13 = Оазз7 = Озаз7 = 2^5а1 + АСаз^*С7?, ^зз7Й = Д7ЙЗЗ = Б1ззз = Бзз1з = 0, Озззз = Е?
где а, в, 7, ^ = 1, 2.
Для контравариантных компонент напряжений получим:
стк1 _ скгс^З _ скгс^З П е
С С г] С С Ог]тп е тп •
Матрица жесткости при этом собирается следующим образом:
^ = С°г СР] А]тп Щтп Е^р^С,
где г, с =1, 8; а, 6, г, т, п, о, р = 1, 3.
Теоретические исследования и опыт применения подобных КЭ свидетельствуют о хорошем влиянии на точность результатов понижения аппроксимаций деформа-
ег3
аппроксимаций [1. 2].
е13
13 128 1
1 еГе? (1+е2е2) (1 + е?е?) х
х (1 + е!еО (1 + е2е2) + е?е1 (1 + е!еЧ (1 + е2е2) (1+е2е2) (1+е?е?)
Первая редукция аппроксимации состоит в исключении переменности этой деформации по толщине:
туг,т _ 1 ^гт
Е13 = 128 Х
е1е? (1+е2е2) (1 + йе1) (1 + е2е2) +
+е?е1 (1 + е^1) (1+е2е2) (1+е2е2)
Вторая редукция направлена на исключение переменности этой деформации е1
1
трт,т _ - у-
128 -е 1 + (е2+е2) е2+е2е2 (е2)21 (е1е? + е?е1).
Третья модификация состоит в нахождении этой деформации в виде линейной аппроксимации по е2 по двум точкам, соответствующим е2 = ±1- После несложных преобразований получим:
Щ1зт = ^ хт [(1 + е2е2) + (е2 + е2) е2] (е1е?+е?е1).
По аналогии можем записать
Е2зт = ^ хт [(1 + + (е1 + е1) е1] (е2е?+е?е2).
Рис. 2
Табл. 1
Wt Сетка 2 х 2 Сетка 3 х 3 Сетка 4 х 4 Сетка 5 х 5
W^ WMOA W,» WMOA WKJ, WMOA W™ WMOA
q 0.68796 0.0046 0.661834 0.00972 0.677329 0.01669 0.683437 0.025391 0.68631
F 0.61152 0.0042 0.5297 0.009 0.5743 0.0154 0.5909 0.0233 0.5987
2. Тестовые задачи
В настоящем разделе рассматриваются примеры численной реализации решения линейных задач. Рассмотрены две тестовые задачи для апробации методики, представленной в предыдущем разделе, приводится сравнение с известными результатами [1, 2. 13].
Задача 1. Рассматривается процесс деформирования пластины (см. рис. 2) под действием равномерного давления q = 1 кГ/см2 или сосредоточенной силы F = 2000 кГ, приложенной в центре пластины. Пусть пластина квадратная со стороной a = 100 см и толщи ной h = 1 см, со следующими механическими свойствами: E = 2 • 106 кГ/см2, коэффициент Пуассона ^ = 0.3. Пластина имеет жесткое защемление по всем боковым граням. В силу симметрии была рассмотрена четверть пластины, использовались различные конечно-элементные сетки. Результаты решения данной задачи приведены в табл. 1, где Wt — максимальный прогиб пластины (см), полученный из приближенного решения [13].
Задача 2. Деформирование под действием собственного веса цилиндрической панели (см. рис. 3), шарнирно опертой по криволинейным границам и свободными прямолинейными гранями со следующими параметрами: L = 1524 см ; R = 762 см; h = 7.62 см (толщина панели); E = 2.1 • 105 кГ/см2; ^ = 0; 7h = 0.044 кГ/см2 (удельный вес).
В силу наличия двух плоскостей симметрии была рассмотрена лишь четверть цилиндрической панели, использовались различные конечно-элементные сетки. Скорость сходимости решения этой задачи, полученная путем использования построенного настоящего конечного элемента, в сравнениях с теоретическим решением н решением обычным классическим трехмерным нзопараметрнческнм вось-миузловым элементом представлена в табл. 2, где Wt — максимальный прогиб цилиндрической панели (см), полученный из приближенного решения [1, 2].
Рис. 3
Табл. 2
Wt Сетка 4 х 4 Сетка 5 х 5 Сетка 10 х 10 Сетка 13 х 13 Сетка 20 х 20
W^ WMOA WKJ, WMOA WKJ, WMOA WKJ, WMOA WKJ, WMOA
9.24 0.5758 8.5715 0.7215 8.6834 1.4159 8.9584 1.8847 9.0108 3.0776 9.0595
Заключение
Предложенная в настоящей работе методика модификации трехмерного восьми-узлового изоиараметрического КЭ теории упругости позволяет получить специальный КЭ, при помощи которого возможно рассчитывать оболочки средней толщины с использованием однослойной аппроксимации по толщине. Приведенные числовые примеры демонстрируют существенное улучшение точности по сравнению с традиционным КЭ трехмерной теории упругости.
Summary
A.I. Golovanov, М.К. Sagdatullin. Three-dimensional Finite Element for Analysis of Tliin-sliell Constructions.
The article regards the construction of a new finite element for calculating the middle thickness of shells on the basis of a modification of three-dimensional isoparametric 8-node element through introduction of a hypothesis of infinit.esimalit.y of compression strains and usage of the technique of approximation order reduction. The method of double approximations on superconvergence points has been applied. Effectiveness of the given approach is shown on numerical examples.
Key words: shell finite element., elastic strains, metric tensor, method of double approximation, hypothesis of infinitesimalit.y of compression strains.
Литература
1. Голованов А.И., Песогиин А.В., Тюлеиева О.Н. Современные конечно-элементные модели и методы исследования тонкостенных конструкций. Казань: Казан, гос. ун-т. 2005. 442 с.
2. Голованов А.И., Тюлеиева О.Н., Шигабутдинов А.Ф. Метод конечных элементов в статике и динамике тонкостенных конструкций. М.: ФИЗМАТЛИТ. 2006. 392 с.
3. Yang H.T.Y., Saigal S., Masud A., Kapania R.K. A survey of recent, shell finite elements // Int. J. for numerical methods in engineering. 2000. V. 47. P. 101 127.
4. Сахаров А.С., Кислоокий B.H., Киричевский В.В., Альте.нбах И., Габберт У., Даи-ке.рт Ю., Кепплер X., Кочы/к 3. Метод конечных элементов в механике твердых тел. Киев: Вища шк., 1982. 480 с.
5. Бережной Д.В. Искривленный конечный элемент пластин и оболочек средней толщины с учетом обжатия // Труды XVII междупар. копф. по теории оболочек и пластин. Казань: Казан, гос. уп-т, 1996. Т. 2. С. 94 99.
6. Гуриелидзе М.Г., Голованов А.И. Расчет толстостенных оболочек с учетом больших деформаций // Труды XVII междупар. копф. по теории оболочек и пластин. Казань: Казан, гос. уп-т, 1996. Т. 2. С. 118 123.
7. Голованов А.И., Гуриелидзе М.Г. Пошаговая постановка решения геометрически нелинейной задачи МКЭ // Динамические и технологические проблемы механики конструкций и сплошных сред. М., 1998. С. 82 87.
8. Баженов В.А., Сахаров А.С., Цыхаповекий В.К. Момептпая схема метода конечных элементов в задачах нелинейной механики сплошной среды // Прикл. механика. 2002. Т. 38, № 6. С. 24 63.
9. Конюхов А.В., Коноплев Ю.Г. Модель термогиперупругостии ее применение к исследованию потери устойчивости раздуваемых пластин // Изв. вузов. Авиац. техника. 2006. № 4. С. 7 13.
10. Kara N.. Kumbasar N. Three dimensional finite element for thick shells of general shape // Int. J. for Physical and Engineering Sciences. 2001. V. 52. P. 1 7.
11. de Sousa R.J. A., Cardoso R.P.R., Valente R.A.F., Yoon J.-W., Gracio J. J., Jorge R.M.N. A new one-point, quadrature enhanced assumed strain (EAS) solid-shell element, with multiple integration points along thickness: Part. I geometrically linear applications // Int.. J. for Numerical Methods in Engineering. 2005. V. 62, No 7. . P. 952 977.
12. Sze K. Y. Three-dimensional continuum finite element, models for plate/shell analysis // Prog. Struct.. Engng Mater. 2002. V. 4. P. 400 407.
13. Тимошенко С.П., Войновекмй-Кргш'.р С. Пластипки и оболочки М.: Наука, 1966. 636 с.
Поступила в редакцию 24.02.09
Голованов Александр Иванович доктор физико-математических паук, профессор кафедры теоретической механики Казанского государственного университета.
E-mail: Alexandr.GolovanovQksu.ru
Сагдатуллин Марат Камилевич аспирант кафедры теоретической механики Казанского государственного университета.
