CC BY
108
2016, Т. 158, кн. 4 С. 530-543
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА. СЕРИЯ ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
ISSN 1815-6088 (Print) ISSN 2500-2198 (Online)
УДК 539.3
МЕТОД МНОГОСЕТОЧНЫХ КОНЕЧНЫХ ЭЛЕМЕНТОВ В РАСЧЕТАХ ТРЕХМЕРНЫХ ОДНОРОДНЫХ И КОМПОЗИТНЫХ ТЕЛ
А.Д. Матвеев
Институт вычислительного моделирования СО РАН, г. Красноярск, 660036, Россия
Аннотация
В работе предложен метод многосеточных конечных элементов для расчета упругих трехмерных однородных и композитных тел при статическом нагружении. Предлагаемый метод построен на основе алгоритмов метода конечных элементов с применением однородных и композитных трехмерных многосеточных конечных элементов (МнКЭ). Рассмотрены процедуры построения МнКЭ. имеющего форму прямоугольного параллелепипеда и сложную форму. Достоинства МнКЭ состоят в том, что они учитывают по правилам микроподхода неоднородную и микронеоднородную структуры тел, описывают трехмерное напряженно-деформированное состояние (без упрощающих гипотез) в однородных и композитных телах, порождают дискретные модели малой размерности и позволяют получать численные решения с малой погрешностью.
Ключевые слова: композиты, упругость, трехмерные тела, многосеточные конечные элементы, микроподход, малая погрешность
Введение
При анализе напряженно-деформированного состояния тел с неоднородной структурой широко используют микро- и макроподходы [1]. Реализация макроподхода для композитных тел регулярной структуры сводится к сложной проблеме нахождения эффективных модулей упругости. Однако для трехмерных композитов нерегулярной структуры определение эффективных модулей упругости связано с большими трудностями, особенно для композитных тел с малыми коэффициентами наполнения. Расчет композитных тел по методу конечных элементов (МКЭ) [2-4] в трехмерной постановке задачи теории упругости [5] с учетом их сложной неоднородной и микронеоднородной структуры по правилам микроподхода, когда каждая компонента неоднородной структуры тела представляется однородными трехмерными конечными элементами (КЭ), требует построения базовых разбиений высокой размерности. Размерности таких базовых разбиений могут достигать несколько миллиардов (порядка 109 + 1010 узловых неизвестных МКЭ), ширина ленты системы уравнений (СУ) МКЭ - более миллиона. Применение в этом случае программ расчета ЛМ8У8, МЛВТИЛМ [6] и др. затруднительно. Кроме того, решение, полученное для СУ МКЭ высокого порядка, содержит вычислительную погрешность, определить точное значение которой достаточно сложно. Построение аналитических решений для трехмерных тел, имеющих сложную неоднородную и микронеоднородную структуру, связано с большими трудностями. При расчете ряда конструкций и деталей (авиационной и ракетной техники) возникает необходимость строить приближенные решения с малой погрешностью. Это вызвано тем,
что для коэффициентов запаса таких конструкций заданы определенные ограничения. Расчет по МКЭ однородных тел сложной формы, имеющих сложное закрепление и нагружение, с применением уравнений трехмерной теории упругости (без упрощающих гипотез) сводится к построению дискретных моделей высокого порядка, что порождает большие трудности при численной реализации МКЭ.
В настоящей работе предложен метод многосеточных конечных элементов для анализа напряженного состояния трехмерных однородных и композитных тел с различными коэффициентами наполнения при статическом нагружении. Предлагаемый метод реализуется на основе алгоритмов МКЭ с применением однородных и композитных многосеточных конечных элементов (МнКЭ). Отличие МнКЭ от существующих КЭ состоит в следующем. При построении т-сеточного КЭ используются т вложенных сеток. Мелкая сетка порождает мелкое (базовое) разбиение, которое учитывает неоднородную структуру и форму МнКЭ, остальные т — 1 крупная сетка применяются для понижения размерности МнКЭ, причем с увеличением т размерность - общее число узловых неизвестных - МнКЭ уменьшается.
Особенность и достоинство МнКЭ состоят в том, что при построении трехмерных композитных МнКЭ (без увеличения их размерности) можно использовать сколь угодно мелкие базовые разбиения композитных тел, состоящих из односе-точных КЭ 1-го порядка, то есть по сути используется микроподход в конечно-элементной форме [1]. Такие базовые разбиения позволяют учитывать в МнКЭ, следовательно, в базовых дискретных моделях композитных тел, неоднородную и микронеоднородную структуру, сложную форму, сложный характер нагружения и закрепления, и описывать сколь угодно точно напряженно-деформированное состояние уравнениями трехмерной теории упругости без введения дополнительных упрощающих гипотез.
Краткая суть МнКЭ [7, 8] заключается в следующем. На базовом разбиении (на мелкой сетке) т-сеточного КЭ, т > 2, определяем его полную потенциальную энергию как функцию Г многих переменных, которыми являются узловые перемещения мелкой сетки. На остальных т — 1 крупной сетке (вложенных в мелкую сетку) строим согласно алгоритмам МКЭ функции перемещений, которые используем для понижения размерности функции Г, что позволяет получать МнКЭ малой размерности. При построении МнКЭ используются функции перемещений в виде степенных и лагранжевых полиномов различных порядков и уравнения трехмерной задачи теории упругости [5], записанные в локальных декартовых системах координат. Полиномы Лагранжа эффективно используются при построении МнКЭ пластинчатого и балочного типов. Приводятся процедуры построения МнКЭ формы прямоугольного параллелепипеда и сложной формы, пластинчатого и балочного типов. В настоящей работе предлагаются смешанные многосеточные дискретные модели, которые применяются для расчета трехмерных однородных и композитных тел сложной формы.
Перечислим достоинства предлагаемых МнКЭ:
• учитывают неоднородную, микронеоднородную структуру и сложную форму тел;
• описывают трехмерное напряженное состояние в однородных и композитных телах;
• образуют многосеточные дискретные модели трехмерных тел, размерность которых в 104 ^ 105 раз меньше размерностей базовых моделей;
• порождают численные решения с высокой скоростью сходимости к точным решениям, что позволяет строить решения с малой погрешностью.
Расчеты показывают, что реализация МКЭ для многосеточных дискретных моделей требует в 105 + 106 раз меньше объема памяти ЭВМ, чем для базовых моделей, при малых временных затратах. Приведен пример расчета однородной
Рис. 1. Сетки 2сКЭ Уа (Ур)
и многослойной пластин с применением трехсеточных КЭ пластинчатого типа и базовых дискретных моделей, которые имеют 2.2 млрд узловых неизвестных МКЭ.
Отметим, что результаты настоящей работы анонсированы в [9]. В отличие от [9] в настоящей статье более подробно описаны процедуры построения МнКЭ, их свойства и достоинства. Рассматривается также вопрос о сходимости приближенных решений, построенных с применением МнКЭ, к точным решениям (см. замечание 7), изложена процедура построения МнКЭ высокого порядка формы треугольной прямой призмы, представлены расчеты многослойной пластины и однородной пластины. Проведен сравнительный анализ полученных результатов.
1. Первая процедура построения композитных двухсеточных КЭ
Основные положения данной процедуры покажем (не теряя общности суждений) на примере построения композитного двухсеточного КЭ (2сКЭ) Уа 2-го порядка в форме куба размерами 6Н х 6Н х 6Н, рис. 1, Охуг - локальная декартова система координат. Считаем, что между компонентами неоднородной структуры 2сКЭ связи идеальны. Функции перемещений, напряжений и деформаций компонентов (которые являются однородными изотропными телами) удовлетворяют закону Гука и соотношениям Коши, отвечающим трехмерной задачи теории упругости [5], то есть во всей области 2сКЭ реализуется трехмерное напряженное состояние без упрощающих гипотез.
Пусть 2сКЭ моделирует параллепипед, армированный волокнами, параллельными оси Ох, сечения волокон на рисунке закрашены. Для построения 2сКЭ используем две вложенные сетки: мелкую На и крупную На. Равномерная сетка На с шагом Н порождена базовым разбиением Ка 2сКЭ, которое состоит из однородных односеточных КЭ (1сКЭ) Уе 1-го порядка в форме куба со стороной Н [2-4]. Базовое разбиение учитывает согласно правилам микроподхода [1] неоднородную структуру 2сКЭ, его форму, нагружение и закрепление, е = 1,... ,М, М - общее число КЭ Уе, рис. 1, М = 216. На сетке На определяем равномерную крупную сетку На (размерности п х П2 х пз) с шагом 3Н, узлы сетки На на рис. 1 отмечены точками, П1 = П2 = пз = 3. Полную потенциальную энергию иа разбиения Ка 2сКЭ Уа представим в виде
М < 1
а
е=1
и а = £ 2 qe \Ke\q е - Ре , (1)
где [Ке] - матрица жесткости, Ре - вектор узловых сил, q е - вектор узловых перемещений 1сКЭ Уе, верхний индекс Т означает транспонирование.
С помощью полиномов Лагранжа [3] (2-го порядка) на сетке На определяем функции перемещений иа, уа, юа 2сКЭ, которые запишем в форме
п п п
Па = ^ Мр ир, Уа = ^ Мр Ур, Юа = ^ Мр Юр, (2)
р=1 в=1 в=1
где ир, Ур, юр, Мр - перемещения и функция формы в-го узла сетки На, п -общее число узлов сетки На, п = П1П2П3 (для 2сКЭ Уа 2-го порядка п = 27, рис. 1).
Пусть qа = {и1,..., ип, ..., уп, Ю1,..., юп}т - вектор узловых перемещений крупной сетки На, то есть вектор узловых перемещений 2сКЭ. Используя (2), вектор q е выражаем через вектор q а, в результате получим равенство
q е ={Аае^ а, (3)
где {Аа] - прямоугольная матрица, е = 1,..., М.
Подставляя (3) в (1), из условия диа/дqа = 0 получаем для 2сКЭ Уа соотношение {Ка^ а = Еа, отвечающее равновесному состоянию, где
м м
[Ка] ^{Аа]т [Ке][Аа], ¥а = ][>а]Т Ре, (4)
е=1 е=1
[Ка], Еа - матрица жесткости и вектор узловых сил 2сКЭ Уа соответственно.
Порядок 2сКЭ Уа при па = П1 = П2 = пз равен па — 1. Для рис. 1 па = 3.
Замечание 1. Решение, построенное для крупной сетки На 2сКЭ Уа, с помощью формулы (3) проецируем на мелкую сетку На базового разбиения 2сКЭ, что дает возможность вычислять напряжения в 1сКЭ Уе базового разбиения 2сКЭ, то есть определять напряжения в любом компоненте неоднородной структуры тела.
Замечание 2. При построении аппроксимирующих функций перемещений для 2сКЭ на крупной сетке На можно использовать степенные полиномы 1-го, 2-го и 3-го порядка [2-4].
2. Вторая процедура построения композитных двухсеточных КЭ
Вторую процедуру кратко рассмотрим (не теряя общности суждений) на примере построения 2сКЭ Ур в форме куба размерами 6Н х 6к х 6к, который имеет неоднородную структуру и расположен в локальной декартовой системе координат Охуг, как и 2сКЭ Уа (рис. 1). В процедуре используем мелкую сетку На и базовое разбиение Яа 2сКЭ Уа. На базовом разбиении Яа 2сКЭ Ур с помощью метода конденсации [3] строим суперэлемент О. Полную потенциальную энергию ир суперэлемента О запишем в форме
ир = 2 qT [Кд М д — qтg Рд, (5)
где [Кд ], Рд, q д - матрица жесткости, векторы узловых сил и перемещений суперэлемента соответственно.
На границе суперэлемента О определяем крупную сетку Нр , вложенную в мелкую На. Пусть q р - вектор узловых перемещений сетки Нр, то есть узловых перемещений 2сКЭ Ур . Используя функции перемещений, построенные на крупной сетке Нр, между векторами q д, q р установим связь вида
q д = {Ap]q
(6)
Рис. 2. Трехсеточный КЭ Уь
где [Ар] - прямоугольная матрица.
Подставляя (6) в (5), из условия дUp/дqр = 0 получаем для 2сКЭ Ур следующее соотношение [Кр^р = Рр, которое отвечает равновесному состоянию, где
[Кр] = [Ар]т [Кд ][Ар ],
А Р
[АР ]Т Рд,
здесь [Кр] - матрица жесткости, Рр - вектор узловых сил 2сКЭ УР.
Замечание 3. Как показывают расчеты, 2сКЭ Ур (построенные по 2-й процедуре) позволяют получить более точные приближения, чем 2сКЭ Уа (построенные по 1-й процедуре). Однако 2-я процедура включает построение суперэлементов О, связанное с обращением матриц высокого порядка. Это увеличивает временные затраты на построение 2сКЭ Ур.
Замечание 4. В силу (3) размерность вектора qa (размерность 2сКЭ Уа) не зависит от числа М, то есть от размерности разбиения Ка. Поэтому для учета в 2сКЭ Уа сложной неоднородной и микронеоднородной структуры можно использовать сколь угодно мелкие базовые разбиения Ка, состоящие из 1сКЭ Уе (см. п. 1). В этом случае в 1сКЭ Уе сколь угодно точно описывается трехмерное напряженное состояние (без упрощающих гипотез).
При резком увеличении размерностей базовых разбиений 2сКЭ Уа , то есть числа М (размерности суперэлемента О, см. п. 2), резко увеличиваются временные затраты на построение 2сКЭ Уа (Ур). В этом случае целесообразно применять трехсеточные КЭ, построение которых требует меньше временных затрат и которые порождают дискретные модели тел меньшей размерности, чем 2сКЭ.
3. Процедура построения композитных трехсеточных КЭ
Основные положения процедуры рассмотрим на примере построения трехсе-точного КЭ (3сКЭ) Уь 2-го порядка в форме куба размерами 12Н х 12Н х 12Н, который расположен в локальной декартовой системе координат Охуг, рис. 2. Область 3сКЭ состоит из 2сКЭ УП, п = 1,... , N - общее число 2сКЭ (размерами 6Н х 6Н х 6Н), для рис. 2 имеем N = 8. Крупные сетки 2сКЭ образуют мелкую сетку Нь 3сКЭ. На мелкой сетке Нь определяем крупную сетку Нь 3сКЭ размерности То1 х Ш2 х тз . Узлы сетки Нь на рис. 2 отмечены точками (27 узлов), т1 = т-2 = тз = 3.
Функции перемещений пь, Уь , юь 3сКЭ, построенные на сетке Нь с помощью полиномов Лагранжа, запишем в форме
т т т
пь = £ МвУь = £ Nчв, юь = £ Жеч^, (7)
в=1 в=1 в=1
где Чр, Цр, ч^, М в - перемещения и функция формы в -го узла сетки Нь, в = = 1,...,т, т - общее число узлов сетки Нь, т = Ш1Ш2 тз, для рис. 2 имеем т = 27.
При построении 3сКЭ Уь используем три характерные сетки: две сетки 2сКЭ У а и крупную сетку Нь. Полную потенциальную энергию Шь 3сКЭ представим выражением
* /1 \
Шь = ]Т (2(чП)т[ККп - (qan)Tрп] , (8)
где [Ка] - матрица жесткости и q<n, РП - векторы узловых перемещений и сил 2сКЭ Уа.
Пусть q ь - вектор узловых перемещений сетки Нь. Используя (7), выражаем узловые перемещения вектора qа через узловые перемещения вектора q ь крупной сетки Нь :
qan = АП] qь, (9)
где [Аьп] - прямоугольная матрица, п = 1,..., М.
Подставляя (9) в функционал (8) и минимизируя его по перемещениям q ь, для 3сКЭ Уь получаем соотношение [Кь^ь = Рь, которое отвечает равновесному состоянию, где [Кь] - матрица жесткости, р - вектор узловых сил 3сКЭ Уь, определяемые по формулам
N N
[Кь] = ]Г[АП]Т [Кап][АП], = ]Т[АП]Т РП. (10)
=1 =1
Порядок 3сКЭ Уь при пь = т1 = т-2 = тз равен пь — 1. На рис. 2 пь = 3.
В силу (9) размерность вектора q ь (то есть размерность 3сКЭ Уь) не зависит от общего числа М 2сКЭ У а , значит, разбиение 3сКЭ Уь на 2сКЭ У а может быть сколь угодно мелким (то есть размеры 2сКЭ У а могут быть сколь угодно малыми). В этом случае в базовых разбиениях 2сКЭ У а, то есть в 1сКЭ Уе (см. п. 1, 2) сколь угодно точно описывается трехмерное напряженное состояние (без упрощающих гипотез). Определение напряжений в случае 3сКЭ Уь сводится к следующему. Пусть найден вектор qь. По формуле (9) находим векторы qа узловых перемещений 2сКЭ У а и по алгоритмам п. 1, определяем напряжения в случае 1сКЭ Уе базового разбиения 2сКЭ У^, п = 1,..., М.
Замечание 5. С помощью 3сКЭ по процедуре, аналогичной п. 3, проектируются четырехсеточные КЭ и т-сеточные КЭ [10, 11], т > 4. Отметим, что т -сеточные КЭ порождают дискретные модели меньшей размерности, чем (т — 1)-сеточные КЭ. Предлагаемый подход проектирования т-сеточного КЭ можно рассматривать как процедуру последовательных вложений к-сеточных КЭ в (к + 1)-сеточный КЭ, где к = 1,. .. ,т — 1.
Замечание 6. При проектировании 3сКЭ пластинчатого и балочного типов мы используем полиномы Лагранжа, построенные на крупных сетках 3сКЭ. Полиномы Лагранжа для 3сКЭ пластинчатого (балочного) типа имеют по осям Ох,
Оу, лежащих на нижней поверхности пластины (по оси балки Оу), более высокий порядок аппроксимаций, чем по оси Ох; т-1, т-2 > тз (по осям Ох, Ох; т-2 > т-1, тз, см. п. 3).
Двумерные композитные 2сК и 3сКЭ, имеющие форму прямоугольника, проектируются по процедурам, которые аналогичны процедурам п. 1-3. При построении композитных 2сКЭ и 3сКЭ в форме треугольника используются полиномы 1-го, 2-го и 3-го порядка [2-4, 7].
Однородные МнКЭ проектируются по процедурам п. 1-3 в предположении, что все компоненты МнКЭ есть однородные изотропные тела с одинаковыми модулями упругости.
Замечание 7. При использовании мелких разбиений композитные тела представляются однородными МнКЭ. Расчеты показывают, что МнКЭ п-го порядка (п = 1, 2, 3) в форме прямоугольного параллелепипеда (спроектированные с применением мелких разбиений) для однородных тел порождают решения, совпадающие с решениями, построенными с применением однородных односеточных КЭ п -го порядка в форме прямоугольного параллелепипеда (таких же размеров, как МнКЭ), узловые сетки которых совпадают с крупными сетками МнКЭ. Как известно [2-4], однородные односеточные КЭ п-го порядка в форме прямоугольного параллелепипеда порождают решения, которые сходятся (в пределе) к точным. Следовательно, решения, построенные для тел с неоднородной структурой с применением однородных и композитных МнКЭ 1-го, 2-го и 3-го порядка в форме прямоугольного параллелепипеда, при Н^ ^ 0, где Н^ - максимальный характерный размер МнКЭ, сходятся к точным решениям.
4. Смешанные многосеточные дискретные модели композитных тел
Краткая суть совместного применения одно- и многосеточного моделирования процесса деформирования трехмерных упругих тел сложной формы с неоднородной структурой состоит в следующем. Подобласть тела, которая включает границу сложной формы или сложного закрепления тела, представляем (мелким) базовым разбиением, которое состоит из однородных 1сКЭ 1-го порядка и учитывает неоднородную структуру и сложную форму тела. Остальную часть тела представляем крупным разбиением, состоящим из композитных МнКЭ в форме прямоугольного параллелепипеда. Мелкое и крупное разбиения соединяем с помощью связующих МнКЭ. В результате получаем смешанную многосеточную дискретную модель тела, размерность которой меньше размерности базовой модели. Подробно смешанные многосеточные дискретные модели трехмерных композитных тел сложной формы рассмотрены в работе [12]. Смешанные многосеточные дискретные модели эффективно применяются при расчете трехмерных однородных тел сложной формы, имеющих сложное закрепление.
5. Проектирование МнКЭ сложной формы
5.1. Рассмотрим процедуру, которая эффективна при использовании мелких разбиений тел на МнКЭ, на примере двумерного 2сКЭ АБС Б 1-го порядка, рис. 3. Криволинейная граница АБ 2сКЭ совпадает с частью границы тела. Узлы А, Б, С, Б крупной сетки 2сКЭ отмечены точками. При построении функции перемещений на крупной сетке 2сКЭ используем полином Р[(х,у) 1-го порядка вида
Л(х, у) = а1 + а2х + азу + а^ху.
то
■ 1 с
О х
Рис. 4. 2сКЭ 2-го порядка рядка
Сложную форму 2сКЭ учитываем с помощью базового разбиения. На рис. 3 показана сетка базового разбиения 2сКЭ, учитывающая его сложную форму. При построении функций перемещений на крупных сетках МнКЭ используем степенные полиномы [2-4]. При этом часть узлов крупных сеток лежит на криволинейных границах МнКЭ. На рис. 4 показан двумерный 2сКЭ 2-го порядка сложной формы, узлы крупной сетки отмечены точками (8 узлов), 5 узлов лежат на криволинейной границе ABC 2сКЭ, которая есть часть границы тела. При построении функций перемещений на крупной сетке 2сКЭ (рис. 4) используем полином P2(x,y) 2-го порядка
Pi(x,y) = ai + a2x + азу + a4xy + a5x2 + aey2 + a7x2y + a8 xy2.
На рис. 5 криволинейная граница ABCD трехмерного 2сКЭ 1-го порядка (крупная сетка имеет 8 узлов, которые отмечены точками) совпадает с границей тела. Функции перемещений на крупной сетке 2сКЭ (рис. 5) строим с помощью полинома Pi(x,y,z) вида
Pi (x, y, z) = ai + a2x + a¡y + a4z + a^xy + a^xz + a^zy + a%xyz.
Отметим, что 2сКЭ (МнКЭ) эффективно используются при расчете однородных трехмерных тел. С помощью базового разбиения 2сКЭ (МнКЭ) учитывается сложная форма, сложный характер закрепления и нагружения трехмерных однородных тел.
Кроме того, 2сКЭ (МнКЭ) порождают дискретные модели однородных тел, размерности которых меньше размерностей базовых моделей.
5.2. Рассмотрим процедуру построения МнКЭ в форме прямой призмы, основания которых имеют сложную форму [13]. В процедуре используются аппроксимирующие функции перемещений МнКЭ (построенные на крупных сетках МнКЭ), которые представлены в виде произведения степенных и лагранжевых полиномов. Суть процедуры рассмотрим на примере 2сКЭ в форме прямоугольного параллелепипеда с отверстием сложной формы, рис. 6. Сечение отверстия заштриховано, узлы крупной сетки отмечены точками (48 узлов), которая по оси Oy является равномерной.
Базовое разбиение 2сКЭ учитывает его неоднородную структуру и форму. В плоскости Oxz (локальной системы координат Oxyz) на узлах крупной сетки 2сКЭ (12 узлов) определяем с помощью МКЭ аппроксимирующие функции перемещений, используя полином Ps(x, z) 3-го порядка вида:
P^(x, z) = ai + a2x + a3z + a4xz + a5x2+
9 9 90 о о Я / \
+a6z + a^x z + a^xz + agx z + aioxz + aiix + ai2z . (11)
Рис. 5. порядка
Рис. 7. 2сКЭ в форме треугольной прямой призмы
Базисную функцию (х,у,г) для в-го узла крупной сетки 2сКЭ, представляем в виде произведения степенных и лагранжевых полиномов, то есть
(х,у,г) = Мр (х,г)Ьр (у) (12)
где Мр (х,г) и Ьр(у) - базисная функция степенного полинома (11) и полином Лагранжа соответственно, отвечающие в-му узлу крупной сетки, в = 1, . . . , N, N - общее число узлов крупной сетки 2сКЭ. На рис. 6 N = 48, полином Ьр(у) имеет 3-й порядок.
Описанный подход упрощает процедуру построения МнКЭ, имеющего форму треугольной прямой призмы (рис. 7), с применением аппроксимаций высокого порядка [13].
На рис. 7 узлы крупной сетки 2сКЭ отмечены точками (24 узла), которая по оси Оу является равномерной. При построении функций перемещений в плоскости Охх на узлах крупной сетки 2сКЭ, 6 узлов на рис. 7, используем полином Р^(х, г) 2-го порядка (записанный в локальной декартовой системе координат Охуг) вида
Р^(х, г) = а1 + а2х + азг + а4хг + а^х2 + а^г2. (13)
Для рис. 7 полином Лагранжа Ьр(у) имеет 3-й порядок, для (12) имеем в = = 1,..., 24. При построении функций перемещений в треугольном сечении призмы на узлах крупной сетки 2сКЭ можно использовать также степенные полиномы 1-го и 3-го порядка [2-4].
6. Результаты численных экспериментов
6.1. Расчет многослойной пластины. Рассмотрим в декартовой системе координат Охуг модельную задачу изгиба (прямоугольной в плане) шестислойной пластины У\ размерами 216Нх х 432Ну х 6НХ (рис. 8), где Нх = Ну = 3.5; Нг = 6; Н = 6НХ - толщина пластины. При у = 0 выполнены граничные условия: и = V = = и> = 0. Слои являются изотропными однородными телами и имеют толщину
Рис. 8. Размеры пластины Vi (V0)
Нс = Н/6 = Нг = 6. Модули продольной упругости (модули Юнга) слоев, начиная с нижнего, соответственно равны 10, 25, 40, 55, 70, 85. Коэффициент Пуассона
равен 0.3. При у > 324Ну, г = Н пластина нагружена вертикальной равномерной поверхностной нагрузкой Рг = 0.0003.
В расчетах используем базовые дискретные модели Е1,..., Е^ пластины V. Базовая модель ЕЩ состоит из одинаковых 1сКЭ Vекп 1-го порядка размерами НЩ х НП х НП, где е = 1,..., Мп, Мп - общее число КЭ V^n базовой модели ЕП, п = 1,..., 6,
hxn = hx/(2n - 1), hl = hy/(2n - 1), hzn = hz/(2n - 1).
(14)
Сетка базовой модели Rl имеет размерность ml х ml х ml, где
= 216(2n-1)+1,
= 432(2n-1) + 1,
= 6(2n-1) + 1, n = 1,..., 6. (15)
С помощью базовых моделей Rl проектируем многосеточные дискретные модели Rl пластины Vi ( n = 1,..., 6 ). Модели Rl состоят из пластинчатых композитных 3сКЭ Vl в форме прямоугольного параллелепипеда размерами 108hl х х 108hl х 6hl, где i - порядковый номер 3сКЭ (см. п. 3). Крупная сетка 3сКЭ Vl имеет по осям Ox, Oy, Oz соответственно шаги 36hl, 36hl и 6hl, то есть полиномы Лагранжа имеют 3-й порядок по осям Ox, Oy и 1-й - по оси Oz. Трех-сеточный КЭ Vl состоит из 2сКЭ W" 3-го порядка в форме куба размерами 6hl х 6hl х 6hl (см. п. 3), j - порядковый номер 2сКЭ, j = 1,..., 324. Крупная сетка 2сКЭ Wl имеет 32 узла, по осям Ox, Oy, Oz - соответственно шаги 2hxl, 2hl, 2hl. Двухсеточный КЭ строим по процедуре п. 2. На крупной сетке 2сКЭ W" аппроксимирующие функции перемещений определяем с помощью степенных полиномов 3-го порядка [4], записанных в локальной декартовой системе координат.
Результаты расчетов приведены в табл. 1, где wl, al - максимальные прогиб и эквивалентное напряжение дискретной модели Rl пластины V1, n = 1,..., 6. Параметры ôl(%), Al(%) находим по формулам
Sl(%) = 100% • al - a1n_1\/al, Al(%) = 100% • \wl - wl_i\/
n = 2,..., 6.
Напряжения аЩ вычисляем по 4-й теории прочности.
Анализ результатов расчетов показывает (см. табл. 1), что характер изменения параметров дЩ(%), ДЩ(%) имеет высокую скорость равномерной сходимости максимальных напряжений аЩ и перемещений иЩ к точным решениям, п = 1,..., 6. Так как величины Д^ = 0.00038, д\ = 0.015 малы и параметры дЩ(%), ДЩ(%) имеют высокую скорость изменения, то перемещение и^ = 767.838 и напряжение ад = 1.166 определены с малой погрешностью. Тогда значения и\, а1 можно считать искомым решением. Напряжения аЩ возникают в окрестности крепления пластины в верхнем слое.
1
2
3
m
m
m
l
l
l
l
Табл. 1
Максимальные прогибы и эквивалентные напряжения моделей Кп ше-стислойной пластины У1
Rn Ri R2 Rs R4 R5 Re
wi 685.854 754.407 764.893 767.851 768.131 767.838
дП(%) - 9.087 1.371 0.385 0.036 0.038
0.686 0.948 1.071 1.125 1.153 1.166
si (%) - 27.637 11.485 4.800 2.428 1.115
Табл. 2
Максимальные прогибы и эквивалентные напряжения моделей Qn однородной пластины Уо
Q n Qi Q2 Qs Q4 Q5 Qe
wi 702.659 775.147 783.436 785.781 786.068 785.905
Ai(%) - 9.352 1.058 0.298 0.037 0.021
„0 1.305 1.750 1.921 1.989 2.017 2.029
- 25.429 8.902 3.419 1.388 0.591
Базовая дискретная модель R0 пластины имеет 2270396480 узловых неизвестных ( ~ 2.2 • 109 неизвестных МКЭ), ширина ленты системы уравнений (СУ) МКЭ равна 955967. Размерность многосеточной модели Д равна 318384, ширина ленты СУ МКЭ - 14699. Реализация МКЭ для многосеточной модели R6 требует в 463770 раз меньше объема памяти ЭВМ, чем для базовой модели R6 пластины V.
6.2. Расчет однородной пластины. Проведен расчет однородной изотропной пластины V0 (рис. 8), которая имеет размеры и закрепление, как и для пластины V. Модуль Юнга пластины Vo равен 100, коэффициент Пуассона -0.3. При z = h, y > 324hy пластина Vo нагружена равномерной поверхностной нагрузкой Pz = 0.0009. В расчетах используем базовые Qn и многосеточные Qn дискретные модели однородной пластины Vo, n = 1,..., 6, правила измельчения базовых разбиений (14), (15) и однородные 3сКЭ типа Vn размерами 108hn х 108hn х 6hzn, которые состоят из однородных 2сКЭ 3-го порядка в форме куба размерами 6hn х 6hn х 6hn (см. п. 6.1). Результаты расчетов пластины Vo даны в табл. 2, где w'n, о~П - максимальные прогиб и эквивалентное напряжение дискретной модели Qn, n = 1,..., 6. Параметры ôn (%), ДП(%) определяем по формулам
ôn(%) = 100% >n - a0_i\/aOn, An(%) = 100% • W - w°n-i\K, n = 2,..., 6.
Из анализа результатов расчетов следует (см. табл. 2), что характер изменения параметров ôn(%), Д^%) демонстрирует высокую скорость равномерной сходимости максимальных напряжений о~П и перемещений w'n к точным решениям. Так как величины Д° = 0.00021, ô0 = 0.00591 малы и параметры ôn(%), ДП(%) имеют высокую скорость изменения, то максимальные перемещение w0 = 785.905 и эквивалентное напряжение а0 = 2.029 можно считать точным решением. Следует отметить, что поскольку ôn > ôn, n = 2,..., 6, то скорость сходимости напряжений an к точному решению, отвечающих однородной пластине Vo, выше скорости
сходимости напряжений an1 , построенных для композитной пластины V1 (при прочих равных условиях).
Заключение
Предложен метод многосеточных конечных элементов для численного анализа напряженно-деформированного состояния трехмерных однородных и композитных тел с различными коэффициентами наполнения при статическом нагружении. Предлагаемый метод основан на методе конечных элементов с применением однородных и композитных многосеточных конечных элементов. Описаны процедуры построения двухсеточных и трехсеточных конечных элементов. Приведенный пример расчета однородной и композитной пластин с помощью трехсеточных конечных элементов показал высокую эффективность их применения.
Литература
1. Фудзии Т., Дзако М. Механика разрушений композиционных материалов. - М.: Мир, 19S2. - 232 c.
2. Зенкевич О. Метод конечных элементов в технике. - М.: Мир, 1975. - 541 с.
3. Норри Д., де Фриз Ж. Введение в метод конечных элементов. - М.: Мир, 19S1. -304 с.
4. Секулович М. Метод конечных элементов. - М.: Стройиздат, 1993. - 442 с.
5. Самуль В.И. Основы теории упругости и пластичности. - М.: Высш. шк., 19S2. -264 с.
6. Голованов А.И., Тюленева О.Н., Шигабутдинов А.Ф. Метод конечных элементов в статике и динамике тонкостенных конструкций. - М.: ФИЗМАТЛИТ, 2006. - 391 с.
7. Матвеев А.Д. Некоторые подходы проектирования упругих многосеточных конечных элементов. - Красноярск, 2000. - 30 с. - Деп. в ВИНИТИ 24.11.2000, № 2990-В00.
S. Матвеев А.Д. Многосеточное моделирование композитов нерегулярной структуры с малым коэффициентом наполнения ll Прикл. механика и техн. физика. - 2004. -№ 3. - С. 161-171.
9. Matveev A.D. Multigrid finite element method in stress analysis of three-dimensional elastic bodies of heterogeneous structure Ц IGF Conf. Ser.: Mater. Sci. Eng. - 2016. - V. 15S, No 1 . - Art. 012067, F. 1-9.
10. Матвеев А.Д. Построение сложных многосеточных элементов с микронеоднородной структурой ll Тез. докл. XXIII Всерос. конф. «Численные методы решения задач теории упругости и пластичности». - Новосибирск: Параллель, 2013. - С. 142-144.
11. Матвеев А.Д. Построение сложных многосеточных конечных элементов с неоднородной и микронеоднородной структурой ll Изв. Алт. гос. ун-та. Сер. Матем. и механика. - 2014. - Т. S1, № 1-1. - С. S0-S3.
12. Матвеев А.Д. Смешанные дискретные модели в анализе упругих трехмерных неоднородных тел сложной формы Ц Вестн. Перм. нац. исслед. политехн. ун-та. - 2013. -№ 1. - С. 1S2-195.
13. Матвеев А.Д. Расчет трехмерных композитных балок сложной формы с применением двухсеточных конечных элементов Ц Вестн. КрасГАУ. - 2015. - № S. - С. 92-9S.
Поступила в редакцию 1S.10.16
Матвеев Александр Данилович, кандидат физико-математических наук, старший научный сотрудник
Институт вычислительного моделирования СО РАН
Академгородок, д. 50, стр. 44, г. Красноярск, 660036, Россия E-mail: mtv@icm.krasn.ru
ISSN 1815-6088 (Print) ISSN 2500-2198 (Online) UCHENYE ZAPISKI KAZANSKOGO UNIVERSITETA. SERIYA FIZIKO-MATEMATICHESKIE NAUKI (Proceedings of Kazan University. Physics and Mathematics Series)
2016, vol. 158, no. 4, pp. 530-543
Multigrid Finite Element Method in Calculation of 3D Homogeneous and Composite Solids
A.D. Matveev
Institute of Computational Modeling, Siberian Branch, Russian Academy of Sciences, Krasnoyarsk, 660036 Russia E-mail: mtv@icm.krasn.ru
Received October 18, 2016 Abstract
In the present paper, a method of multigrid finite elements to calculate elastic three-dimensional homogeneous and composite solids under static loading has been suggested. The method has been developed based on the finite element method algorithms using homogeneous and composite three-dimensional multigrid finite elements (MFE). The procedures for construction of MFE of both rectangular parallelepiped and complex shapes have been shown. The advantages of MFE are that they take into account, following the rules of the microapproach, heterogeneous and microhomogeneous structures of the bodies, describe the three-dimensional stress-strain state (without any simplifying hypotheses) in homogeneous and composite solids, as well as generate small dimensional discrete models and numerical solutions with a high accuracy.
Keywords: composites, elasticity, three-dimensional solids, multigrid finite elements, mic-roapproach, high accuracy
Figure Captions
Fig. 1. 2gFE grids Va (Vp).
Fig. 2. Three-grid FE Vb.
Fig. 3. Two-dimensional 2gFE of the first order.
Fig. 4. 2gFE of the second order.
Fig. 5. Three-dimensional 2gFE of the first order.
Fig. 6. 2gFE with an outlet.
Fig. 7. 2gFE of triangular right prism shape.
Fig. 8. Vi (Vo) plate sizes.
References
1. Fudzii T. and Dzako M. Fracture Mechanics of Composite Materials. Moscow, Mir, 1982. 232 p. (In Russian)
2. Zenkevich O. Finite Element Method in Technology. Moscow, Mir, 1975. 541 p. (In Russian)
3. Norrie D., de Vries J. An Introduction into Finite Element Method. Moscow, Mir, 1981. 304 p. (In Russian)
4. Sekulovich M. Finite Element Method. Moscow, Stroyizdat, 1993. 664 p. (In Russian)
5. Samul' V.I. Fundamentals of the Theory of Elasticity and Plasticity. Moscow, Vyssh. Shk., 1982. 264 p. (In Russian)
6. Golovanov A.P., Tyuleneva O.N., Shigabutdinov A.F. The Finite Element Method in Statics and Dynamics of Thin-Walled Structures. Moscow, FIZMATLIT, 2006. 391 p. (In Russian)
7. Matveev A.D. Some Approaches for Constructing Elastic Multigrid Finite Elements. Krasnoyarsk, 2000. 30 p. Dep. VINITI, Nov. 24, 2000, no. 2990-V00. 2000. (In Russian)
8. Matveev A.D. Multigrid models of composite materials of irregular structure with a small filling ratio. J. Appl. Mech. Tech. Phys., 2004, vol. 45, no. 3, pp. 440-448. doi: 10.1023/B:JAMT.0000025027.41912.a7.
9. Matveev A.D. Multigrid finite element method in stress analysis of three-dimensional elastic bodies of heterogeneous structure. IOP Conf. Ser.: Mater. Sci. Eng., 2016, vol. 158, no. 1, art. 012067, pp. 1-9.
10. Matveev A.D. Constructing composite multigrid elements with microheterogeneous structures. Tez. dokl. XXIII Vseros. konf. "Chislennye metody resheniya zadach teorii uprugosti i plastich-nosti" [Proc. XXIII All-Russian Conf. Numerical Methods for Solving the Problems of Elasticity and Plasticity Theory]. Novosibirsk, Parallel, 2013, pp. 142-144. (In Russian)
11. Matveev A.D. Constructing the composite multigrid elements of heterogeneous and microhetero-geneous structures. Izv. Altai. Gos. Univ., Ser. Mat. Mekh., 2014, vol. 81, no. 1-1. pp. 80-83. (In Russian)
12. Matveev A.D. Combined discrete models in the three-dimensional elastic inhomogeneous bodies analysis of complex shape. Vestn. Permsk. Nats. Issled. Politekh. Univ. Mekh., 2013, no. 1, pp. 182-195. (In Russian)
13. Matveev A.D. The calculation of the three-dimensional irregular-shaped composite beams using the double-grid finite elements. Vestn. Krasnoyarsk. Gos. Agrar. Univ., 2015, no. 8, pp. 92-98. (In Russian)
/ Для цитирования: Матвеев А.Д. Метод многосеточных конечных элементов I в расчетах трехмерных однородных и композитных тел // Учен. зап. Казан. ун-та. \ Сер. Физ.-матем. науки. - 2016. - Т. 158, кн. 4. - С. 530-543.
/ For citation: Matveev A.D. Multigrid finite element method in calculation of 3D ( homogeneous and composite solids. Uchenye Zapiski Kazanskogo Universiteta. Seriya \ Fiziko-Matematicheskie Nauki, 2016, vol. 158, no. 4, pp. 530-543. (In Russian)
