CC BY
26
УДК 539.3
М. К. Сагдатуллин
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ КОНЕЧНОГО ЭЛЕМЕНТА
ДЛЯ РАСЧЕТА ТОНКОСТЕННЫХ КОНСТРУКЦИЙ
Ключевые слова: метод конечных элементов, упругие деформации, метод двойной аппроксимации, метрический тензор,
гипотеза малости напряжений обжатия.
Работа посвящена построению нового конечного элемента для расчета средней толщины оболочки на основе модификации трехмерного изопараметрического восьмиузлового элемента путем введения гипотезы малости напряжений обжатия и использования техники понижения порядка аппроксимации, при этом использован метод двойной аппроксимации по точкам суперсходимости. На численных примерах показана эффективность данного подхода.
Keywords: a method offinite elements, elastic strains, method of double approximation, metric tensor, hypothesis of infinitesimality of
compression strains.
The article regards the construction of a new finite element for calculating the middle thickness of shells on the basic of a modification of three-dimensional isoparametric 8-node element through introduction of a hypothesis of infinitesimality of compression strains and usage of the technique of approximation order reduction. The method of double approximations on super convergence points has been applied. Effectiveness of the given approach is shown on numerical examples.
Введение
Данная работа посвящена построению трехмерного конечного элемента оболочки, способного к определению НДС тонкостенных конструкций. Основные этапы построения предлагаемого КЭ состоят в следующем. Фрагмент оболочки представляется как трехмерное тело. Вводятся изопараметрические полилинейные аппроксимации геометрии и неизвестных компонент перемещений. Узловые точки располагаются в вершинах выделенного параллелепипеда.
Использование аппроксимаций перемещений по толщине эквивалентно введению кинематической гипотезы об их линейном распределении по толщине с учетом обжатия оболочки. Вариация потенциальной энергии деформации (работа внутренних сил напряжений на возможных деформациях) вычисляется численно по квадратурной формуле Гаусса-Лежандра второго порядка. В каждом квадратурном узле вводится статическая гипотеза, обычно используемая для оболочек средней толщины, об отсутствии влияния поперечного напряжения на мембранные, изгибные и деформации поперечного сдвига. Фактически задача поперечного обжатия отделяется от классической задачи теории оболочек с учетом поперечного сдвига. Главный источник ошибок, возникающих при использовании подобных КЭ, которые делают их практически непригодными - появление «ложных» деформаций поперечного сдвига. Предлагается использование метода двойной аппроксимации по точкам суперсходимости в трехмерной постановке.
1. Геометрические параметры
Рассмотрим восьмиузловой трехмерный конечный элемент в виде параллелепипеда.
Введем в рассмотрение радиус-вектор, определяющий геометрию поверхности в виде
Я) = X' )ё,.
Для проекций X' соответствующие аппроксимации запишем в виде
X' ))Х'М, ),
г=1
где введены полилинейные функции формы
Ж, )=1 (1)(1)(1+г),
Соотношение, определяющее ковариантные базисные вектора определим в виде
Я = X = п'ё
Наряду с базисными векторами Як вводится сопряженный (контравариантный) базис
* _ _ 1
Rk e =.
=kmnRm x Rn = Rk%
дХ' ' 2^0
Скалярные произведения базисных векторов определяют ковариантные и контравариантные компоненты метрических тензоров, которые характеризуют так называемую внутреннюю геометрию поверхности (площади, длины дуг, углы между касательными).
8 дХт дХт 8
О = щ = у дХ_ = у ятят,
т=1 т=1
О'1 = Я'Я1 = у = у Я'-тЯ1т,
т=1 дХт дХт £1 '
=У ОтпЯЯ .
т,п
2. Метод двойной аппроксимации
Для борьбы с явлением заклинивания получили весьма широкое распространение различные модификации так называемого метода двойной аппроксимации. Суть его в следующем: деформации внутри конечного элемента
аппроксимируются самостоятельно, отдельно от перемещений. При этом если правильно выбрать степени полиномов перемещений и деформаций, можно весьма существенно понизить величину погрешности.
Впервые подобная технология была предложена А.С. Сахаровым в его известной моментной схеме конечных элементов (МСКЭ), как способ более правильного представления смещений, как твердого целого.
Связь между вектором деформаций е и вектором перемещений и в линейной теории упругости задается линейным кинематическим уравнением
е = Еи, (1)
где Е - матрица дифференциального оператора.
Если ввести обозначение
Е™ = 1XХГ & N + ^ ), (2)
] 2 £ ' {д? д? д? д? )
то компоненты тензора деформаций можно представить в следующем виде
е = Х X и7Е™. (3)
г=1 т=1
Теоретические исследования и опыт построения подобных конечных элементов свидетельствуют о «благотворном» влиянии на точность результатов понижения аппроксимаций деформаций поперечного сдвига е3. Опишем технику введения «понижения порядка аппроксимаций» для трехмерного конечного элемента.
Рассмотрим деформации поперечного сдвига е13, е23. Из соотношений выше имеем
е™ = 128 хт [? ? (1+?? )(+?3? )• •(1+??)(+??2)+?? (1+?? )•
•(1+??2)+??2 )(+??)].
Е2™ = 118X™ [? ? (1 + ?? ) (1 + ?? ) ^ •(1+?? )+??)+?? (1+?? )•
•(1+??)+?? )+?? )].
Первая «редукция аппроксимации» состоит в исключении переменности этих деформаций по толщине. Получим
Е™ = 128 хт [? ? (1+??2 )(1+?? )•
•(1+??2)+?? (1+?? )(1+?? )(1+??)].
Е2т = 128 хт [? ? (1+?1?) (+??) •
•(1+??2)+?? (1+??1 )(1+?? )(1+??)].
Вторая «редукция» направлена на исключение переменности вдоль координат ? , ? деформаций е13, е23 соответственно.
1 + (? +?)) + ?? (? )2 "•
•?1? + ??).
1 + (? + ?)? + ?? (?)2"• •(?? + ???).
Третья модификация предполагает определить деформации е13, е23 в виде линейной
аппроксимации по ? и ? , по двум точкам,
? = +1 и ? = +1 соответственно. В итоге получим
е™ = 128 хт [(1+?? )+(? +?)) ]• •(?? +???).
Е2т=118 хт [(1+??)+(? +?)? ]• •(?? +?3?).
Следует заметить, что брать меньшие степени в выражениях для деформаций, чем те которые представлены выше, не рекомендуется, так как это приведет к пониженному рангу матрицы жесткости, или, иначе говоря, к дополнительным ложным жестким смещениям.
Некоторым недостатком описанной выше схемы, является нарушение вариационного принципа, в результате чего для усеченных деформаций не выполняются условия совместности деформаций, что приводит к немонотонной сходимости. Тем не менее, выбор таких аппроксимаций является целесообразным и оправданным с позиций МСКЭ [3].
3. Гипотеза малости напряжений обжатия
Используем представление закона Гука для изотропного материала, учитывая, что в криволинейном базисе роль единичного тензора играет метрический тензор.
] 2^]+^?%, (4)
где первый инвариант тензора деформации имеет вид
1и=(в )--(е) = От"епт.
Отсюда получаем «классическое уравнение» теории упругости
а'] =( 2/иО'тО'" + АО']От" )етп.
Егт =
13
128
-хт
77 гт _ А
Е23 =-Л,
128 '
Рассмотрим процедуру учета малости напряжений обжатия. Полученные, таким образом конечно-элементные модели утрачивают характер объемного поведения КЭ и соответствуют специальной теории оболочек. Формально запишем эту гипотезу в виде равенства
= Яз -И-Яз = 0.
Подставим в это уравнение выражения тензора напряжений
И1
^ = ИОО
получим
Из = 2ие33 +Л°зз Л* =
= 2/ие33 +Л033 (0% + 20'3е33 + О33е33) = 0,
где ', 1 = 1,2. Выразим из этого уравнения деформацию обжатия е33.
Л033
-(О% + 20%).
33 2и+Ю3303
Теперь преобразуем соотношения упругости для ковариантных компонент напряжений отдельно для напряжений ац (', 1 = 1,2), характеризующих мембранное и изгибное напряженное состояние и а3, описывающих напряжения поперечного сдвига.
Итак, получаем
а.
1 = 2ив„ +Л0'/(0тП£тп + 20т3^т3 + 033^33 ) = = 2иеу (и'О^тп + 2и*0а£'3),
а, = 2ива +Я0'3 ("Хп + 2и'0т3ет3),
где
и = 1 -
АО33О
2и + Ю3303
т, п = 1,2.
Так как в качестве степеней свободы в рассматриваемом КЭ фигурируют узловые степени свободы на лицевых поверхностях, то фактически в кинематике учитываются возможные обжатия оболочки (изменение ее толщины). Следовательно, в расчетную схему необходимо ввести «упрощенный» закон упругости, связывающий напряжение обжатия а33 с деформацией обжатия е33. Для этого введем следующее соотношение
И33 = Е'^
где Е' - модуль жесткости на обжатие (в общем случае, он может быть определен из экспериментальных данных). В частности, можно принять
Е* = и((-+2и) А. + и
Запишем в общем виде соотношения упругости
а = Бе.
В дискретном виде соотношение выглядит следующим образом
3
а.. = У Б., е , ', 1 = 1,3
'1 '1тп тп
где
Баруз = 2и8ау8рз +10а,и'075,
Ба.ру3 = Ба.р3у = Л^0арU''0УЪ, Ба3у8 = Б3аЗ = а3ии ,
Ба3у3 = Б3ау3 = Ба33у = Б3а3у = 2и8ау +^0а3и0\
Б33 8 = Б 833 = Б 333 = Б33 3 = 0,
33у8 у533 у 333 33у 3 >
Б3333 = Е
где а, р, у,8 = 1,2 .
Для контравариантных компонент напряжений получим
И = О^ОИ: = 0к'0,1Б: е .
'1 '1тп тп
Анализ результатов решения обширного класса задач теории тонких оболочек и оболочек средней толщины [1, 2, 4-8] показывает, что применение линейной аппроксимации сдвига по толщине стенки, с учетом гипотезы малости напряжений обжатия, приводит к достаточно точным результатам.
4. Составление матрицы жесткости
Примем в качестве базового соотношения вариационное уравнение принципа виртуальных работ для однородного и изотропного материала
| а'18е1 =| / - 8ШУ + | !п8ШБ, (5)
V V Ба
где V - физический объем, ограниченный
поверхностью £ = Б" и , Б" п Ба = 0, У - вектор
заданных внешних объемных сил, - вектор
заданных напряжений на части поверхности Ба, на которой определены силовые граничные условия. Кинематические условия на части поверхности Б" выполняются за счет специальным образом определенных аппроксимаций.
Потенциальная энергия деформации (работа внутренних напряжений на возможных деформациях) вычисляется численно по квадратурной формуле Гаусса - Лежандра. Введем систему квадратурных точек, которые в локальной системе координат обозначим ^, ^ и
соответствующий весовой множитель ау. Если
предполагается использование упомянутой выше восьми точечной формулы Гаусса - Лежандра (второго порядка по каждой из координат), то
"73
У = 1,8; = ±—= ; а у = 1. Далее интеграл в левой
части уравнения виртуальных работ заменяется конечной суммой, то есть
¡а^де^йУ = Щаде1]40й?с1?й? = ^атде„я№, г
V V /
де ? = ?.
Рассмотрим один из вариантов свертки в подынтегральном выражении левой части уравнения виртуальных работ (5)
дЖ = аде у = {2/в'тОу" + Ш'уОт" )етпдеу (6) Используя представление (6), получим
дЖ = адеи = £ £ икгди[Е™Е1 •
к ,1 =1 г, í=1
•(2/0'т0у" +ЯО'уОт")) = ]Г икБ^дП!,
к ,1=1 г, í=1
где
Б" = {2/в'тОу" +АОиОт" )ЕгтпЕтт^О.
Суммируя значения величин Б« , вычисленных в
системе квадратурных точек ?, получаем блоки
соответствующей модифицированной матрицы жесткости трехмерного восьми узлового КЭ линейной задачи теории упругости. Зная компоненты матрицы жесткости всех конечных элементов тонкостенной конструкции и решив систему линейных алгебраических уравнений, можно определить компоненты перемещений, а, следовательно, и деформированное состояние всей конструкции.
5. Числовые примеры
Настоящий раздел посвящен примерам численной реализации решения линейных задач. Рассмотрено несколько тестовых примеров для апробации методики, представленной в предыдущих разделах, приводятся сравнения с аналитическими решениями задач, численными решениями ППП АШУЯ.
Задача 1. Деформирование пластины (см. рис. 1) под действием равномерного давления д = 1 кГ / см2 и сосредоточенной силы Е = 2000 кГ . Пластина -квадратная со стороной а = 100 см и толщиной И = 1 см , со следующими механическими свойствами: модуль Юнга Е = 2 • 106 кГ / см2, коэффициент Пуассона / = 0.3. Пластина имеет жесткое
защемление по всем боковым граням. В силу симметрии была рассмотрена четверть пластины, использовались различные конечно - элементные сетки. Результаты решения данной задачи приведены в таблице 1, где Жт - максимальный прогиб пластины, полученный из приближенного решения [8].
Рис. 1
Таблица 1
Задача 2. Одной из самых распространенных задач является задача об изгибе замкнутой круговой цилиндрической оболочки со свободными торцами под действием самоуравновешенной системы двух сосредоточенных сил (см. рис. 2). Особенностью этой задачи является близость напряженного состояния к чистому изгибу, что представляет определенные трудности при использовании МКЭ. Наиболее обширные численные данные накоплены для оболочки со следующими параметрами:
Ь = 26.2 см Я = 12.5 см И = 0.238 см
Е = 7.4-105 кГ-
см
/ = 0.3125 Р = 45.3 кГ
Аналитическое решение, полученное в предположении нерастяжимости оболочки [8], дает значение прогиба под силой Жт = -0.275 см. Конечно - элементное решение дает несколько большее значение [1,2]. В силу наличия трех плоскостей симметрии была рассмотрена восьмая часть цилиндрической оболочки. Результаты решения данной задачи на различных конечно -элементных сетках рассмотрены в таблице 2.
Таблица 2
Жт 10х10 Ж класс 10х10 Ж " мод 15х15 Ж класс 15х15 Ж мод 20х20 Ж класс 20х20 Ж мод
0.275 0.0132 0.2773 0.0267 0.2826 0.0427 0.2842
3х3 3х3 4х4 4х4 5х5 5х5
Жт Ж класс Ж мод Ж класс Ж мод Ж класс Ж мод
д 0,6879 0,0079 0,6751 0,01669 0,6834 0,0254 0,68631
Е 0,6115 0,009 0,5743 0,0154 0,5909 0,0233 0,5987
Заключение
Предложенная в настоящей главе, методика модификации трехмерного восьмиузлового изопараметрического конечного элемента теории упругости позволяет получить специальный КЭ, при помощи которого вполне реально рассчитывать НДС оболочек средней толщины сложной геометрии с использованием однослойной аппроксимации по толщине. Хочется отметить, что «классический» трехмерный элемент не дает удовлетворительного решения при однослойной аппроксимации, даже с использованием густой сетки. Приведенные числовые примеры демонстрируют существенное улучшение точности по сравнению с традиционным КЭ трехмерной теории упругости.
Литература
1. Голованов А.И., Песошин А.В., Тюленева О.Н. Современные конечно-элементные модели и методы исследования тонкостенных конструкций. Казань: Казан. гос. ун-т, 2005. - 442 с.
2. Голованов А.И., Тюленева О.Н., Шигабутдинов А.Ф. Метод конечных элементов в статике и динамике тонкостенных конструкций. М.:ФИЗМАТЛИТ, 2006. -392 с.
3. Сахаров А.С., Кислоокий В.Н., Киричевский В.В., Альтенбах И., Габберт У., Данкерт Ю., Кепплер Х., Кочык З. Метод конечных элементов в механике твердых тел. Киев: Вища школа, 1982. - 480 с.
4. Голованов А.И., Бережной Д.В. Метод конечных элементов в механике деформируемых твердых тел. Казань, 2001, 300 с.
5. Бережной Д.В., Сагдатуллин М.К., Саченков А.А. Универсальный конечный элемент для расчета комбинированных конструкций // Вестник Казанского государственного технологического университета. -2012. - №17. - С.150-157.
6. Бережной Д.В., Сагдатуллин М.К. Трехмерный конечный элемент для расчета оболочек средней толщины // Вестник Казанского государственного технологического университета. - 2013. - №9. - С.256-261.
7. Баженов В.А., Сахаров А.С., Цыхановский В.К. Моментная схема метода конечных элементов в задачах нелинейной механики сплошной среды // Прикл. механика. - 2002. - Т. 38, № 6. - С. 24-63.
8. Тимошенко С.П., Войновский - Кригер С. Пластинки и оболочки. М.: Наука. - 1966. - 636 с.
© М. К. Сагдатуллин, канд. физ.-мат. наук, доцент каф. теоретической механики и сопротивления материалов, Казанский национальный исследовательский технологический университет, [email protected].
© M. K. Sagdatullin, candidate of physico-mathematical sciences, associate professor, Kazan National Research Technological University, department of theoretical mechanics and strength of materials, [email protected].
