CC BY
35
УДК 539.3 АД. Матвеев
МЕТОД МНОГОСЕТОЧНЫХ КОНЕЧНЫХ ЭЛЕМЕНТОВ
A.D. Matveev
MULTIGRID FINITE ELEMENT METHOD
Матвеев А.Д. - канд. физ.-мат. наук, доц., ст. науч. сотр. Института вычислительного моделирования СО РАН, г. Красноярск. E-mail: [email protected]
Для решения ряда важных физических краевых задач (решения уравнений которых эквивалентны нахождению минимума соответствующие функционалов) предлагается метод многосеточных конечных элементов (ММКЭ), который реализуется на основе соотношений и алгоритмов метода конечных элементов (МКЭ) в форме метода Ритца с применением многосеточных конечных элементов (МнКЭ). При построении n-сеточного конечного элемента (КЭ) используем n вложенных сеток. Мелкая сетка порождена базовым разбиением тела, которое учитывает его сложную форму и физические особенности краевой задачи (например, неоднородную структуру упругого тела). Остальные n -1 сетки применяем для понижения размерности МнКЭ (причем с увеличением n размерность МнКЭ уменьшается). Суть МнКЭ заключается в следующем. На базовом разбиении n-сеточного КЭ, n > 2, которое состоит из известных односеточных КЭ, определяем функционал краевой задачи как функцию F многих переменных, которыми являются значения искомой функции в узлах мелкой сетки. На остальных n-1 сетках строим аппроксимирующие функции, которые используем для понижения размерности функции F , что позволяет проектировать МнКЭ малой размерности. Проектирование n-сеточного КЭ проводится по единой матричной процедуре. Основные отличия ММКЭ от МКЭ состоят в следующем. Во-первых, в ММКЭ можно применять сколь угодно мелкие базовые разбиения тел, что позволяет сколь угодно точно учитывать их сложную форму, неоднородную и микронеоднородную структуру упругих тел (без увеличения размерностей многосеточных дискретных моделей). В МКЭ невозможно использовать сколь угодно мелкие разбиения тел, так как ресурсы ЭВМ ограничены, т.е. ММКЭ более эффективный, чем
Matveev A.D. - Cand. Phys. and Math. Sci., Assoc. Prof., Senior Staff Scientist, Institute of Computing Modeling, SB RAS, Krasnoyarsk. E-mail: [email protected]
МКЭ. Во-вторых, реализация ММКЭ на основе базовых моделей тел требует меньше памяти ЭВМ и временных затрат, чем реализация МКЭ для базовых моделей, т.е. ММКЭ более экономичный, чем МКЭ. В-третьих, в ММКЭ применяем упругие однородные и неоднородные МнКЭ, при построении которых используем системы вложенных сеток, что расширяет область применения ММКЭ. В МКЭ применяют однородные односеточные КЭ. Поэтому можно считать, что ММКЭ есть обобщение МКЭ, т.е. МКЭ - частный случай ММКЭ. Изложены процедуры построения МнКЭ различной формы. Предложена верхняя оценка погрешностей приближенных решений.
Ключевые слова: физические краевые задачи, однородные и неоднородные тела, многосеточные конечные элементы, малая погрешность.
To solve a number of important physical boundary value problems (which solutions of the equations being equivalent to finding the minimum of corresponding functional) the multigrid finite element method (MFEM) which is realized on the basis of ratios and algorithms of the method of final elements (MFE) in the form of Ritz method with application of multigrid final elements (MFEM) was proposed. To construct a n -grid finite element (FE), the n -nested grids were used. A finite grid is generated by basic body partition taking into account its irregular shape and physical features of the boundary value problem (e.g. the in-homogeneous structure of elastic body). Other grids n -1 were used to reduce MFE dimension, and with increasing MFE dimension decreases. The essence of MFE is as follows: at a basic partition of n > 2 grid FE, consisting of known single-grid FE, the functional of boundary value problem was determined as a F function of many variables, being the values of the required function at the nodes of a fine grid. On the
other n -1 grids some approximating functions were used for the decrease of the dimension of F function, allowing one to develop small dimensional MFE. The developing n -grid FE is carried out according to a single matrix procedure There are some essential differences between MFEM and FEM. First, in regard to MFEM, some arbitrarily fine base body partitions can be applied, which makes it possible to take into account their irregular shape heterogeneous and microheterogeneous structure (without increasing the dimensions of the multigrid discrete models). As to FEM, it is impossible to use any arbitrarily fine partitions of the bodies, as the computer resources are limited, i.e. MFEM is more efficient than FEM. Secondly, the implementation of MFEM based on the essential models of bodies takes less computer memory and span time than that of FEM for essential models, i.e. MFEM is more time and memory-saving than FEM. Thirdly, in MFEM some elastic homogeneous and inhomogeneous MFE are applied, using the nested grids to construct, significantly expanding the scope of MFEM. Therefore, MFEM can be assumed to be a generalization of FEM, i.e. FEM is a special case of MFEM. The procedures of developing MFE of various shapes were presented. The top assessment of errors of approximate decisions is offered.
Keywords: physical boundary value problems, homogeneous and inhomogeneous bodies, multigrid finite elements, small error.
Введение. Предлагаемый в данной работе метод многосеточных конечных элементов (ММКЭ) реализуется на основе соотношений и алгоритмов стандартного МКЭ [1, 2] (в форме метода Ритца) с применением МнКЭ. Вначале кратко рассмотрим основные соотношения МКЭ на примере важнейших физических краевых задач, таких как, например, распространение тепла в трехмерной области, кручение упругого стержня некругового сечения, задача о безвихревом течении жидкости [1]. Каждая из этих физических задач описывается дифференциальным уравнением в частных производных вида [1]
Dx
К д
Kxx ~дХ ГдУ
K
+-
д_ dz
др ду
K Q = 0
zz дz ) Q
\
+
с граничными условиями на границе S
Р = Рв
(2)
и (или) на границе S
Kxx ^ + Ky, ^ + Kzz ^ +
xx дx x yy ду y дz z +q + k(p-px ) = 0,
(3)
где ^ + £2 = £ , £ - граница области V тела;
р = р (х, у, х) - искомая функция;
Кх = Кхх (х, у, х); Куу = Куу (х, у, х);
Кгг = Кгг (х, у, х) и д = д (х, у, х) не зависят от
функции р; 1х, 1у, /г - направляющие косинусы
вектора нормали к поверхности £2.
Уравнение (1) применимо как к изотропным телам, так и к анизотропным (неоднородным). Координатные оси должны быть параллельны главным осям инерции в анизотропных областях. Кратко сформулируем задачу (1)-(3) в стандартной конечноэлементной формулировке (в форме метода Ритца) [1]. Как известно, с вариационной точки зрения решение уравнения (1) с граничными условиями (2), (3) эквивалентно отысканию минимума функционала [1]
х =
f 1 [ Kx
дx I
^ л2
+ K
уу
др ду
+ K,
Dz
20р] dV \яр + 0,5Нрр-рх . (4)
£
Согласно МКЭ, область тела V представляем КЭ V(е), е = \,...,Е, Е - общее число КЭ, т.е. V = V1 + V(2) +...+V(Е). Для каждого конечного элемента V(e) вместо искомой функции р введем аппроксимирующую функцию р(е) вида
Pe) = [N (е)]{ф(е)},
(5)
где [N(e)] = [Ne,...,Ner ] - матрица функций фор-
(1)
мы КЭ V(е); Ые - базисная функция КЭ V(е) для
узла г, г = 1,...,г; {Ф(е)} = {Фе,...,Фе>г - вектор узловых значений искомой функции р; г - общее число узлов КЭ V(e\ е = 1,...,Е; т- транспонирование.
Минимизацию функционала (4) проводим на множестве узловых значений функции р, кото-
рые определяются на всей области V. Будем минимизировать функционал перед вычислением интегралов. Такой подход позволяет построить характеристики КЭ, которые отвечают каждой конкретной задаче. В связи с этим преобразуем функционал (4). Введем две матрицы [1]:
(2)}т =
»
дф-е) дф(е) дф дх ду дг
(6)
£(е) =
К(е)
хх 0
0 0
0 0
ку 0
К
(7)
Интегралы в (4) разбиваем на интегралы по
отдельным КЭ ^е), е = 1,...,Е. Тогда, учитывая (6), (7), функционал (4) запишем в форме
Е I 1
% = £ | / ^^[°{е)]- | ф(е)£(е)^ + ¡у(е)д(е)ёБ +
е=1 \ъг (е)
+
. 1 Е
| - Н(е)ф(еф(е) - 2ф(е)фэ+фэ2]ёБ = £%(е),
Ме) 2 е=1
Б 2 У
(8)
где %(е) - вклад элемента V(е) в функционал %; + Se2 = Б(е) - область КЭ V(е). Минимизация функционала % сводится к выполнению соотношения
Е Л„Ле)
д%
д% _ _
д{Ф} £ д{Ф(е)}
= 0,
где {Ф} = {{Ф(1)},...,{Ф(Е)}}т, {Ф} - вектор узловых значений у всей области V.
Используя (5), (6), построим соотношения
{£(е)} = [В(е)]{Ф(е)}
где {Ф(е)} - вектор узловых значений функции у, отвечающий КЭ V(е); [В(е)] - матрица дифференцирования вида
(9)
[ В(е)] =
дЩ дЩ
дЩ
дх дх дх
дЩ дЩ дЩ
ду ' ду ду
дЩ дЩ дЩ
дг дг
дг
(11)
(10)
Используя (10), (11), запишем интегралы (8) по элементам V(е) в виде [1]
%(е) = | -{Ф(е)}т [в(е)]т р(е)] [В(е)]{Ф(е)}й^ - |б(е)[N(е)]{Ф(е)}й^ + | д(е)[N(е)] {Ф(е)} ёБ +
Г/И 2
+
| -{ф«}т[Щ(е)]т[Щ(2)]{Ф(2)} ёБ- | Н(е)ф[Щ(е)]{Ф(е)}ёБ+- |. (12)
Вклад отдельного элемента д%/д{Ф(е)} в
Е
общую сумму £ д%(е)/д{Ф(е)} представим в
д%(е) д{Ф(е)}
= [ К (е)]{Ф(е)} + { /(е)}
(e)^
(13)
е=1
следующей компактной форме:
где [К(е)] - матрица элемента Vм;; {/(е)} -вектор нагрузки элемента V(е) [1],
[К(е)] = | [В(е)]т [0(е)] [В(е)]^ + | Н [ Щ(2)]т [ Щ(е)] ёБ;
V ( е )
а(е)
Б
(е)
Б
Б
Б
{/е)} = _ |(е)]тёУ + |q(e)[N(e)]гdS Н(е)р[^е)]^£
гАе) ?(е) s(e)
(15)
Используя матрицы (14) и векторы нагрузок (15) элементов V(е) по алгоритмам МКЭ [1, 2], строим глобальные матрицу и вектор нагрузки для дискретной области V, т.е. определяем систему уравнений МКЭ для области V.
Цель исследований. Анализ метода многосеточных конечных элементов.
Методы исследований. В данной работе для решения ряда широко известных краевых задач предложен ММКЭ, который базируется на алгоритмах МКЭ в форме метода Ритца с применением МнКЭ. При построении МнКЭ используем мелкие (базовые) разбиения, которые учитывают сложную форму МнКЭ и особенности краевой задачи. Например, для неоднородных упругих тел в МнКЭ учитывается (в рамках микроподхода [3]) их структура. При построении МнКЭ используем аппроксимирующие функции в виде степенных и лагранжевых полиномов различных порядков и функционал, отвечающий краевой задаче и представленный в локальной декартовой системе координат. Достоинства МнКЭ состоят в том, что они учитывают сложную форму тел (в упругих телах сколь угодно точно описывают трехмерное напряженное состояние, учитывают в рамках микроподхода неоднородную и микронеоднород-
ную структуру без увеличения размерностей МнКЭ) образуют дискретные модели, размерность которых в 103 ^ 105 раз меньше размерностей базовых, и порождают приближенные решения с малой погрешностью. Изложены процедуры построения МнКЭ различной формы.
Результаты исследований
1. Процедуры построения двухсеточных конечных элементов. Первую процедуру построения двухсеточного КЭ (2сКЭ) V формы
прямоугольного параллелепипеда (треугольника, прямоугольника) изложим на примере краевой задачи (1)-(3). При построении 2сКЭ используем две вложенные сетки: мелкую и крупную. Мелкая сетка Иа порождена базовым разбиением Яа области 2сКЭ V, которое состоит из односеточных
КЭ (1сКЭ) V(е) первого порядка формы куба и учитывает сложную форму 2сКЭ, физические особенности тела краевой задачи. Например, учитывает неоднородную структуру тела в рамках микроподхода [3], е = 1,...,м, м - общее число
КЭ V(е). На базовом разбиении Яа определяем
функционал Ха, который с учетом обозначений
(12)-(15) для КЭ V(е) запишем в виде
м м Л >
х.=Е Х(е) =£ 11{Ф(е)}т [^(е)]{Ф(е)}+{Ф(е)}т {/(е)}
е=1 е=1 V2 У
(16)
На мелкой сетке \ определяем крупную сетку На 2сКЭ. С помощью полиномов Лагранжа [2] (или степенных полиномов [1, 2]) на сетке Н определяем аппроксимирующую функцию Фа, которую запишем в форме
(17)
,=1
где N¡3, ф - базисная функция и значения функции р в , -м узле сетки На; п - общее число узлов крупной сетки.
Пусть {Фа} - вектор узловых значений р крупной сетки На. Используя (17), вектор {Ф(е)} выражаем через вектор {Ф } . В результате получим равенство
{Ф(е)}=[ л: ]{Ф.},
(18)
где [Л. ] - прямоугольная матрица, е = 1,...,м.
Подставляя (18) в (16), из условия д%а /д{Фа} = 0 получаем для 2сКЭ матричное
соотношение [Ка ]{Ф.} = /}, где
м м
[ ка ] = £[ А ]Т [К (е)][ а; ], ш = -£[ а; ]т {/(е)}
(19)
е=1
е=1
где [К], {/} - матрица и вектор нагрузки 2сКЭ
V; -
Вторую процедуру рассмотрим на примере построения 2сКЭ V, который имеет такую же
форму и размеры, как и 2сКЭ V• В процедуре используем мелкую сетку Иа и базовое разбиение Я 2сКЭ V. С помощью метода конденсации [2] на разбиении Яа строим суперэлемент О. Функционал %3 для 2сКЭ V запишем в виде
1
^ = ^ {Ф3 } [К ]{Ф3} + {Ф3 }т /} , (20)
где [К5 ] - матрица, отвечающая суперэлементу О, {/5}, {Ф33}- векторы нагрузки и узловых значений функции р суперэлемента О -
На границе суперэлемента О определяем крупную сетку Н3, вложенную в мелкую сетку ка -
Пусть {Ф} - вектор узловых значений функции р сетки Н, т.е. узловых значений р 2сКЭ V -Используя аппроксимирующую функцию, построенную на сетке Н, между векторами {Ф33}, {Ф5} установим связь вида
Замечание 1. В силу (18) размерность вектора {Фа} (размерность 2сКЭ V) не зависит от числа
М, т.е. от размерности разбиения Яа. Поэтому
для сколь угодно точного учета, например сложной неоднородной и микронеоднородной структуры или формы 2сКЭ V можно использовать сколь угодно мелкие базовые разбиения Я , состоящие из 1сКЭ У(е\ без увелечения размерности 2сКЭ V
2. Процедура построения п-сеточных конечных элементов. Рассмотрим процедуру построения п -сеточного КЭ УlП формы прямоугольного параллелепипеда (треугольника, прямоугольника), п > 2 . Область КЭ у состоит из (п - 1)-сеточных КЭ yn-1 такой же формы, как КЭ V", 1 = 1,...,М ,, М , - общее число КЭ у-1.
и ' ? ? п—1 п—1 I 1
Крупные сетки Ип-1 КЭ Уln 1 (1 = 1,...МИ_1) образуют сетку И, на которой определяем крупную сетку Ни для КЭ у (сетка Ни вложена в сетку И). Для п -сеточного КЭ у функцию Фп, построенную на крупной сетке Н с помощью полиномов Лагранжа (или степенных полиномов [1, 2]), представим
{ф3}=[ а; ]{ф5},
(21)
где [ А; ] - прямоугольная матрица.
Подставляем (21) в (20). Из условия минимизации функционала %3, т.е. из выполнения
д%3 / д{Ф<.} = 0, следует соотношение [ К5] {Ф5 } = Ж5, где
г5 = - [ а; ]т {/5}, [к ] = [ а; ]т [к5 ][ а; ], здесь [К] - матрица, - вектор узловой нагрузки 2сКЭ V.
Фп = 1 к ф; ,
(22)
;=1
где Ыпр - базисная функция ;-го узла сетки Ни; Ф пр - значение функции Ф в ; -м узле сетки Ни; т - общее число узлов сетки Нп.
Функционал хпа л-сеточного КЭ у представляем как сумму функционалов, построенных для КЭ У-1, т. е. имеем
Мп-1( 1
%а = Е( т{Ф(п-1)}т [к/п-1)]{Ф
1=1 V 2
(п-1)
}+{Ф( п-1)}т {/(п-1)}
(23)
где [К(п 1)], /(п 1), Ф(п-1) - матрица, векторы нагруз- ки и узловых значений искомой функции р КЭ V"1 -
94
Обозначим через Ф^ вектор узловых значений функции р крупной сетки Яи КЭ V'П. Используя (22), вектор Ф(пА) выражаем через вектор Ф(«и), т.е.
ф(«-1) = [ап. ] Ф^, (24)
где [К(п)], /(п) - матрица и вектор нагрузки п-сеточного КЭ V.
Замечание 2. Построение матриц и векторов нагрузок п -сеточных КЭ У^ проводится по единой процедуре (по формулам (25)) с применением матриц [К(п-1)\ /(п-1) элемента уп-1
(г = \...,М„-1) и матриц [Л"], [АП]т, построенных с помощью функции (22) п -сеточного КЭ У". Следует отметить, что КЭ У" порождают дискретные модели тел меньшей размерности, чем КЭ уп 1.
Замечание 3. При построении 2сКЭ пластинчатого (стержневого) типа используем полиномы Ла-гранжа, имеющие по осям Ох, Оу (по оси Оу) более высокий порядок аппроксимации, чем по оси Ох (по осям Ох, Ох).
Замечание 4. При использовании мелких разбиений композитные тела представляются однородными МнКЭ. Расчеты показывают, что МнКЭ п -го порядка (п = 1,2,3) формы прямоугольного параллелепипеда (треугольника, прямоугольника), спроектированные на мелких разбиениях, для однородных тел порождают решения, совпадающие с решениями, построенными с применением
Вторую процедуру покажем на примере построения упругого 2сКЭ У формы прямоугольного параллелепипеда с отверстием сложной формы (рис. 2), расположенного в локальной декартовой системе координат Охух . Область 2сКЭ Уъ
представляем базовым разбиением Яа, которое
учитывает сложную форму, неоднородную структуру 2сКЭ и порождает мелкую сетку. На рисунке
где [Лп] - прямоугольная матрица, г = 1,...,М^.
Подставляя (24) в (23), из условия дхпа / сФ[п = 0 получаем для КЭ У^ матричное соотношение вида [К(п)] Ф(п) = /(п), где
(25)
однородных односеточных КЭ п -го порядка формы прямоугольного параллелепипеда (тетраэдра, треугольника, прямоугольника) таких же размеров, как МнКЭ, узловые сетки которых совпадают с крупными сетками МнКЭ. Как известно [1, 2], однородные односеточные КЭ п -го порядка формы прямоугольного параллелепипеда (тетраэдра, треугольника, прямоугольника) порождают решения, которые при измельчении КЭ (т.е. в пределе) сходятся к точным. Следовательно, решения, построенные для упругих неоднородных тел с применением однородных и композитных МнКЭ п -го порядка, формы прямоугольного параллелепипеда (треугольника, прямоугольника) в пределе сходятся к точным.
3. Процедуры построения двухсеточных конечных элементов сложной формы. Первую процедуру рассмотрим на примере упругого двумерного 2сКЭ у сложной формы (рис. 1). Криволинейная граница ЛСВ совпадает с границей тела. На рисунке 1 показано базовое разбиение 2сКЭ у, которое учитывает его форму. Узлы
крупной сетки отмечены точками. При построении функции перемещений на крупной сетке используем полином р (х, х) вида
(26)
2 в плоскости Охх показано базовое разбиение на 1сКЭ. На мелкой сетке определяем крупную сетку, узлы отмечены точками (60 узлов). Основание 2сКЭ У , которое лежит в плоскости Охх и имеет 12 узлов крупной сетки, рассматриваем квадратным КЭ АВСБ 3-го порядка с вырезом (рис. 2).
Мп-1 Мп-1
^ = й А ]Т [ Кп-1)][ лп ] , /« = - £[ Л" ]т /,(п-1),
г=1 г=1
Р (х, х) = а + а2х + а2 + а4хх + а5х2 + ах2 + а7хх2 + а&х2 х, аг = сотХ.
Рис. 2. 2сКЭ К
Для аппроксимации перемещений КЭ abcd используем полином р(x,z) вида
Pb (x, z) = a + a2x + a3z + a4xz +
2 2 2 2 + a^x + a6z + ax z + asxz +
+a)xzi + a0zx + aux3+^z3, a = const. (27)
Базисную функцию щ (x, y, z) для узла i, j,k (рис. 2) крупной сетки Яа 2сКЭ V представляем в следующей форме:
¥jk (x У, z) = Nk(x, z)Lj(У) <
(28)
где (х, г) - базисные функции перемещений КЭ АБСБ, отвечающие полиному вида (27);
L (у) - полином Лагранжа 4-го порядка [2], i, k = 1,...,4, j = 1,...,5.
Рассмотрим 2сКЭ V формы прямой треугольной призмы с отверстием, которое закрашено (рис. 3). Процедура построения 2сКЭ V аналогична 2-й процедуре п. 3. На рисунке 3 показаны мелкая и крупная сетки 2сКЭ, узлы крупной отмечены точками (24 узла). Основание 2сКЭ, которое лежит в плоскости Oxz и имеет 6 узлов крупной сетки, рассматриваем треугольным КЭ 2-го порядка. Для аппроксимации перемещений КЭ V используем полином вида P (x, z) = a + a2x + a3z + a4xz + ax2 + a6z2, a = const. В формуле (28) применяем полином Лагранжа L (у) 3-го порядка, j = 1,...,4.
Рис. 3. 2сКЭ V
4. Построение смешанных многосеточных дискретных моделей для тел сложной
формы. Пусть трехмерное тело V имеет сложную форму. Отметим, что в этом случае примене-
ние МнКЭ простой формы (формы куба, прямоугольного параллелепипеда, квадрата, прямоугольника) затруднительно. Подобласть у с У тела, которая имеет сложную форму, представляем мелким (базовым) разбиением, которое состоит из 1сКЭ первого порядка (формы куба, тетраэдра, треугольной прямой призмы [1, 2]) и учитывает сложную форму области У . Остальную
представляем круп-
часть тела У
У + У2 = У,
ным разбиением, состоящим из МнКЭ формы куба (прямоугольного параллелепипеда, квадрата, прямоугольника). Мелкое и крупное разбиения соединяем с помощью связующих МнКЭ. В результате получаем смешанную многосеточную дискретную модель, размерность которой меньше размерности базовой дискретной модели. Подробно для упругих тел смешанные многосеточные модели рассмотрены в работах [4, 5]. В [4-20] изложены процедуры построения МнКЭ и применения ММКЭ для расчета однородных и композитных тел.
5. Построение верхней оценки погрешности приближенных решений
5.1. Пусть для заданной последовательности разбиений (многосеточных дискретных моделей) (Я }гп=1 тела У построена по ММКЭ последовательность решений (р}п=1 краевой задачи, Р =р(х,у,х), ^ ^ 0 при г ^да, где - характерный максимальный размер МнКЭ, отвечающий дискретной модели Я, г = 1,...,п, п > 3. Введем относительные погрешности 8, е вида
8 =
\ р0 -Р I I I \Р\ I
е =
| \Рг -Рг-х\ |
I |Р,| I :
(29)
где Сг =
| \Ро -Рг-1 | | | \Рг -Рг-1| |
г = 2.....п.
(32)
Введем величину С* по формуле
С*= | \Рп -Рг-2 | |
' | \Рг-1 -Р-2 | |
г = 3,...,п . (33)
Используя в (31) С* вместо С, получаем верхнюю оценку для относительной погрешности 8 решения р, которое отвечает модели Я, т.е.
8<( 1+С*) е, г = 3,...,п.
(34)
где е, С*, согласно (29), (33), вычисляются с помощью решений р_1, р, р_2, р.
Отметим, что верхняя оценка (34) является приближенной (в (31) вместо Сг применяем С*)
и применяется для мелких разбиений тел.
5.2. При расчете упругих конструкций на прочность определяют максимальные эквивалентные напряжения, с помощью которых находят коэффициенты запаса конструкций. Обозначим: ° -максимальное эквивалентное напряжение, отвечающее многосеточной дискретной модели Я тела У; а0 - максимальное эквивалентное напряжение тела, которое отвечает точному решению, отметим, что °0,8 > 0, г = 1,...,п, п > 3. По
аналогии с (29) введем погрешности 8°, вида
где р - решение, отвечающее модели Я; р0 -точное решение, г = 1,...,п.
Согласно неравенству, треугольника имеем
\ \ро -Рг \ \ < \ к-Х -Рг \ \ , \ \ ро - Рг-1 \ \
I I рг I I
Рг
г = 2,...,п.
Рг
Используя (29), (30), получаем
8 < ( 1 + С ) 8г,
(30)
(31)
-К, = -1 , г = 2,...,п , (35)
где максимальные эквивалентные напряжения сгг._!, ° моделей , Я.
Рассуждая аналогично п. 5.1, для относительной погрешности 8° напряжения ° получаем
верхнюю оценку в форме
8° < Р° = ( 1 + С°)
(36)
а
гпо Г? 1ап -?г -2 I - , „
где С1 =-, 1 = з,...,п -
? -1 -2 1
(37)
Обозначим: ж - максимальное перемещение (некоторой точки), возникающее в модели Я тела V , ж > 0 , 1 = 1,...,п . С помощью максимальных перемещений щ._2, , ж, соответственно моделей Я 2, Я_1, Я для относительной погрешности 8™ перемещения ж получаем верхнюю оценку РЩ вида
8ж < Рж = ( 1 + Сж) еж,
(38)
5.3. Практическое применение полученных верхних оценок рассмотрим на примере расчета шестислойной консольной квадратной (в плане) пластины при локальном нагружении с применением трехсеточных КЭ, описанного в работе [18]-В данном случае 3сКЭ оболочечного типа имеет те же размеры, как в [18], но по осям Ох, Оу использовались полиномы 2-го порядка. Результаты расчетов даны в таблице 1, где ж, ? - максимальные прогибы и эквивалентное напряжение многосеточной модели Я пластины, 1 = 1,...,7 (п = 7). Значения е? (%), еж(%) (при 1 > 2) находим по формулам:
где 8Ж =
1 жо - Ж 1
еЖ =
1Ж-1 - Ж 1
ж
ж
(39)
е?(%) = 100% х |?-?-1|/?; еж(%) = 100% х | щ - | /ж -
ГЖ_ | жп - Ж-2| - , и
С =-, 1 = 3,...,п,
| ж-1 - ж-2 |
(40)
п > 3, ж - максимальное перемещение тела, отвечающее точному решению.
Характер изменения величин е?(%), еж(%) (1 = 2,...,7) демонстрирует быструю сходимость ? и ж соответственно к точным решениям ?,
ж0.
Таблица 1
Я Я1 Я2 Яз Я4 Я5 Яб Я]
ж 399,711 479,121 489,484 492,701 494,109 495,042 495,580
еж (%) - 16,57 2,12 0,65 0,28 0,19 0,11
?1 1,346 2,488 2,789 2,928 3,038 3,108 3,158
е?(%) - 45,90 10,79 4,75 3,62 2,25 1,58
Учитывая, что (см. табл. 1)
ж >... > ж+1 > ж >... > ж, ?0 > ... >?+1 >?> ... >?1
решением, т.е. принимаем
ж0 = ж7 ■
?0 =?7
(43)
(41)
согласно (35), (39), перемещения ж0 и напряжения ? отвечают соотношениям
Ж = (1 + 8К, ?0 = (1 + 8?)?, 1 = 2,...,7 - (42)
Так как величины е^ = 0,0011, е? = 0,0158 малы, то перемещение ж7 = 495,580 и напряжение ?= 3,158 (п = 7 ) можно считать точным
Оценки рж, р? находим, используя данные таблицы 1 и формулы (35)-(40). Для нахождения максимальных перемещений жр и эквивалентных напряжений ?р модели Я с помощью верхних оценок рж, р? в формулах (42) вместо точных значений погрешностей 8™, 8? соответственно используем верхние оценки р™, р?, т.е. имеем
у = (1+Ру)у, а* = (1+Р?)?, г = 3,...,7 . (44) це 2, где параметры а(%), <(%) определяются
по формулам:
Результаты расчетов представлены в табли-
< (%) = 100% х | у* - / <; аа(%) = 100% х | а*/а*, г = 3,...,7 . (45)
Таблица 2
Я3 *4 Рб
у* 512,358 501,027 498,185 497,883 496,966
аГ (%) 4,46 1,66 0,82 0,57 0,28
а* 3,568 3,376 3,440 3,324 3,294
аа(%) 21,82 13,28 11,69 6,51 4,12
рГ (%) 4,67 1,69 0,82 0,57 -
57 (%) 1,245 0,584 0,297 0,109 -
р? (%) 27,92 15,31 13,23 6,96 -
5 а (%) 13,230 7,855 3,950 1,608 -
Согласно расчетам, имеем ау а, ^ 0 при г ^да, что в силу (45) означает: у* ^ у, а* при г ^да. Учитывая (41)-(43), точные значения относительных погрешностей 5у, 5°, которые отвечают соответственно перемещению у и напряжению , модели р, г = 1,...,б, находим по формулам:
5 = У0 - У
5а = а0 -а,
У
а
(46)
Анализ расчетов показывает, что у* > у,
а* > а,, Ру > 5у, ра > 5а, где г = 3,...,6, причем, Ру , Р? ^5? , т. е. Ру ^ 0, ра ^ 0 при г ^да.
6. Многосеточные КЭ оболочек вращения и двоякой кривизны. Рассмотрим кратко построение МнКЭ для композитных оболочек вращения и двоякой кривизны на примере упругого криволинейного двухсеточного КЭ (2сКЭ) V (рис. 4). Для 2сКЭ V введем локальные системы координат: декартовую Охух и криволинейную
О^г/^. Плоскость Ог/£ проходит через ось ей оболочки. Ось Оу лежит в плоскости , причем Оу 1О£. Базовое разбиение 2сКЭ V учитывает его неоднородную структуру и порождает мелкую сетку ^ , на которой определяем для
2сКЭ у крупную криволинейную сетку Я^, узлы которой отмечены точками (36 узлов, рис. 4). Для узлов сетки Я введем криволинейную целочисленную систему координат фк (рис. 4), ф, к = 1,...,4 , г = 1,...,3 , для узла р: г = 2 , ф, к = 1. Сечение 2сКЭ V, лежащее в плоскости рассматриваем криволинейным односе-точным КЭ (1сКЭ) у 3-го порядка, который имеет 12 узлов крупной сетки Я^ (рис. 5). Область 2сКЭ V образуется вращением плоского 1сКЭ У/ вокруг оси ей (оси оболочки) на угол а0, где а - угол раствора 2сКЭ у. 1сКЭ у будем называть образующим КЭ. Радиусы р, Р3 (Я2, Я4) описывают нижние (верхние) границы боковых граней 2сКЭ V (рис. 4).
Рис. 4. Оболочечный 2сКЭ V
Рис. 5. Сечение 2сКЭ Vd, 1сКЭ Vda
При малых углах а0 для точки М(х,у, г) помощью которого аппроксимируются перемеще-
можно принять х = £, г = £ . При построении ния °бразующег° К^ т.а 1сКЭ V; (см. рис. 5).
базисных функций перемещений для узлов круп- Полином Рй(у,С) для 1сКЭ VI имеет вид ной сетки Н используем полиномы Лагранжа
Ц(х) и интерполяционный полином Рй(у,С), с
Pd (y,C) = a1 + a2 У + a3C + a4 yC + a5 У2 + a6C + a7 У + a8 yC2 +
+aX3 + аю y3C+^y3 + ai2^3-
(47)
Li(x) = П (x")/(Xi "), i = 1,...,3,
Базисную функцию w (x, y, C) для узла i, j, k крупной сетки Я 2сКЭ Vd (см. рис. 4) пред- p=hрФ{
ставляем в форме \yijk (x, y,C) = Nß (y,C)Lt(x), где x (xi) - координата узла p (узла i ), лежа-где Njk(y,C) - базисные функции 1сКЭ V/, ко- щего на оси В силу малости угла а0 раствора торые отвечают полиному Pd(y,C) вида (47), 2сКЭ Vd имеем j, k = 1,...,4; Ц (x) - полином Лагранжа 2-го порядка, имеющий для точки M вид
x = aRi, xi = о-Д, xp = apR, (48)
где a - угол раствора (точки M) плоскости, которая проходит через ось cd и включает точку M; a (a) - угол раствора узла p (узла i ), лежащего на оси O% (см. рис. 4).
Учитывая представления (48), получаем
Wyk (a, y,C) = Njk (y,C)L (a),
где Li (a) = ^
a-a
, i = 1,...,3
p=1, p^i ai ap
(49)
Для каждого узла 1, у, к сетки На определим целое число ;, введем обозначение Ыр , ; = 1,...,3б и функции перемещений иа, ус{ , ж 2сКЭ V представим в виде
36
= 1N
ß=1
ßUß
36
d = 1N ß=1
ß vß
w
36
d = XN
ß=1
ß wß=
(50)
U
d
где N, пр, \р, Ур - базисная функция и перемещения р -го узла сетки Я^.
Используя (50), по известным алгоритмам [11-13, 20] определяем матрицу жесткости и вектор узловых сил 2сКЭ V. С помощью 2сКЭ строим 3-сеточные КЭ. При расчете оболочек вращения сложной формы, прямых круговых конических и цилиндрических оболочек постоянной толщины целесообразно при построении криволинейных 2сКЭ в качестве образующих КЭ соответственно
использовать 1сКЭ формы (криволинейного) треугольника (рис. 6, а), параллелограмма (рис. 6, б) и прямоугольника (рис. 6, в) 1-го, 2-го и 3-го порядка [1, 2]. На рисунке 6 показаны образующие 1сКЭ 3-го порядка. На рисунке 6, б показан образующий 1сКЭ, применяемый для расчета конических оболочек с углом 2р при вершине. При построении функций перемещений для образующих 1сКЭ, показанных на рисунке 6, б, в, применяем полиномы вида (47).
б
Рис. 6. Образующие 1сКЭ
Предлагаемый подход эффективен при построении криволинейных МнКЭ, которые применяются в расчетах оболочек двоякой кривизны. Для верификации предлагаемых оболочечных МнКЭ можно использовать известный численный метод [20].
Заключение. В данной работе для решения ряда краевых задач предложен метод многосеточных конечных элементов, который базируется на алгоритмах метода конечных элементов в форме метода Ритца с применением многосеточных конечных элементов. Достоинства многосеточных конечных элементов состоят в том, что они учитывают сложную форму тел (в упругих телах сколь угодно точно описывают трехмерное напряженное состояние, учитывают в рамках микроподхода неоднородную и микронеоднородную структуру без увеличения размерностей МнКЭ), образуют многосеточные дискретные модели, размерность которых в 103 ^ 105 раз меньше размерностей базовых моделей, и порождают приближенные решения с малой погрешностью.
Литература
1. Сегерлинд Л. Применение метода конечных элементов. - М.: Мир, 1979. - 392 с.
Норри Д., Фриз Ж. де. Введение в метод конечных элементов. - М.: Мир, 1981. - 304 с. Фудзии Т., Дзако М. Механика разрушений композиционных материалов. - М.: Мир, 1982. - 232 с.
Матвеев А.Д. Смешанные дискретные модели в анализе упругих трехмерных неоднородных тел сложной формы // Вестн. Перм. нац. исслед. политехн. ун-та. - 2013. - № 1. -С. 182-195.
Матвеев А.Д., Гришанов А.Н. Смешанные многосеточные дискретные модели трехмерных цилиндрических композитных панелей и оболочек сложной формы // Сб. ст. XIX зимней школы по механике сплошных сред. -Пермь, 2015. - С. 198-211. Матвеев А.Д. Многосеточное моделирование композитов нерегулярной структуры с малым коэффициентом наполнения // Прикладная механика и техническая физика. -2004. - № 3. - С. 161-171. Матвеев А.Д. Построение сложных многосеточных элементов с микронеоднородной структурой // Численные методы решения задач теории упругости и пластичности: тез. докл. XXIII Всерос. конф. (Барнаул, 2013 г.). -Новосибирск: Параллель, 2013. - С. 142-144. Матвеев А.Д. Построение многосеточных конечных элементов сложной формы с при-
а
в
менением локальных аппроксимаций // Вестник КрасГАУ. - 2013. - № 1. - С. 28-34.
9. Матвеев А.Д. Расчет тонких пластин и оболочек с применением многосеточных конечных элементов со свободными границами // Вестник КрасГАУ. - 2014. - № 3. - С. 44-47.
10. Матвеев А.Д. Построение сложных многосеточных конечных элементов с неоднородной и микронеоднородной структурой // Известия АлтГУ. Серия: Математика и механика -2014. - № 1/1. - С. 80-83.
11. Матвеев А.Д. Применение граничных двух-сеточных элементов в расчетах трехмерных композитных балок // Вестник КрасГАУ. -2014. - № 5. - С. 44-49.
12. Матвеев А.Д., Гришанов А.Н. Одно- и двух-сеточные криволинейные элементы трехмерных цилиндрических панелей и оболочек // Известия АлтГУ. Серия: Математика и механика. - 2014. - № 1/1. - С. 84-89.
13. Матвеев А.Д., Гришанов А.Н. Многосеточные криволинейные элементы в трехмерном анализе цилиндрических композитных панелей с полостями и отверстиями // Ученые записки Казан. ун-та. Серия: Физико-математические науки. - 2014. - Т. 156. -№ 4. - С. 47-59.
14. Матвеев А.Д., Гришанов А.Н. Многосеточные лагранжевые криволинейные элементы в трехмерном анализе композитных цилиндрических панелей и оболочек // Вестник КрасГАУ. - 2015. - № 2. - С. 75-85.
15. Матвеев А.Д. Расчет трехмерных композитных балок сложной формы с применением двухсеточных конечных элементов // Вестник КрасГАУ. - 2015. - № 8. - С. 92-98.
16. Матвеев А.Д. Расчет композитных пластин и балок с учетом их структуры с применением сложных многосеточных конечных элементов // Вестник КрасГАУ. - 2015. - № 9. - С. 100107.
17. Матвеев А.Д. Метод многосеточных конечных элементов в расчетах композитных пластин и балок // Вестник КрасГАУ. - 2016. -№ 12. - С. 93-100.
18. Матвеев А.Д. Метод многосеточных конечных элементов в расчетах трехмерных однородных и композитных тел // Учен. записки Казан. ун-та. Серия: Физ.-мат. науки. - 2016. - Т. 158. - № 4. - С. 530-543.
19. Matveev A.D. Multigrid finite element method in stress of three-dimensionalelastic bodies of het-
erogeneous structure // IOP Conf. Ser.: Mater. Sci. Eng. - 2016. - V. 158, № 1. - Art. 012067.
- P. 1-9.
20. Матвеев А.Д., Гришанов А.Н. Трехмерные композитные многосеточные конечные элементы оболочечного типа // Изв. АлтГУ. Сер. физико-математические науки. - 2017. - № 4.
- C. 120-125.
Literatura
1. Segerlind L. Primenenie metoda konechnyh jelementov. - M.: Mir, 1979. - 392 s.
2. Norri D., Friz Zh. de. Vvedenie v metod konechnyh jelementov. - M.: Mir, 1981. - 304 s.
3. Fudzii T., Dzako M. Mehanika razrushenij kompozicionnyh materialov. - M.: Mir, 1982. -232 s.
4. Matveev A.D. Smeshannye diskretnye modeli v analize uprugih trehmernyh neodnorodnyh tel slozhnoj formy // Vestn. Perm. nac. issled. politehn. un-ta. - 2013. - № 1. - S. 182-195.
5. Matveev A.D., Grishanov A.N. Smeshannye mnogosetochnye diskretnye modeli trehmernyh cilindricheskih kompozitnyh panelej i obolochek slozhnoj formy // Sb. st. XIX zimnej shkoly po mehanike sploshnyh sred. - Perm', 2015. -S. 198-211.
6. Matveev A.D. Mnogosetochnoe modelirovanie kompozitov nereguljarnoj struk-tury s malym kojefficientom napolnenija // Prikladnaja mehanika i tehnicheskaja fizika. - 2004. - № 3.
- S. 161-171.
7. Matveev A.D. Postroenie slozhnyh mnogosetochnyh jelementov s mikroneodno-rodnoj strukturoj // Chislennye metody reshenija zadach teorii uprugosti i plastichnosti: tez. dokl. XXIII Vseros. konf. (Barnaul, 2013 g.). - Novosibirsk: Parallel', 2013. - S. 142-144.
8. Matveev A.D. Postroenie mnogosetochnyh konechnyh jelementov slozhnoj for-my s primeneniem lokal'nyh approksimacij // Vestnik KrasGAU. - 2013. - № 1. - S. 28-34.
9 . Matveev A.D. Raschet tonkih plastin i obolochek s primeneniem mnogosetoch-nyh konechnyh jelementov so svobodnymi granicami // Vestnik KrasGAU. - 2014. - № 3. - S. 44-47.
10. Matveev A.D. Postroenie slozhnyh mnogosetochnyh konechnyh jelementov s neodnorodnoj i mikroneodnorodnoj strukturoj // Izvestija AltGU. Serija: Matematika i mehanika.
- 2014. - № 1/1. - S. 80-83.
11. Matveev A.D. Primenenie granichnyh dvuhsetochnyh jelementov v raschetah trehmernyh kompozitnyh balok // Vestnik KrasGAU. - 2014. - № 5. - S. 44-49.
12. Matveev A.D., Grishanov A.N. Odno- i dvuhsetochnye krivolinejnye jelementy trehmernyh cilindricheskih panelej i obolochek // Izvestija AltGU. Serija: Matematika i mehanika.
- 2014. - № 1/1. - S. 84-89.
13. Matveev A.D., Grishanov A.N. Mnogosetochnye krivolinejnye jelementy v trehmernom analize cilindricheskih kompozitnyh panelej s polostjami i otverstijami // Uchenye zapiski Kazan. un-ta. Serija: Fiziko-matematicheskie nauki. - 2014. -T. 156. - № 4. - S. 47-59.
14. Matveev A.D., Grishanov A.N. Mnogosetochnye lagranzhevye krivolinejnye jelementy v trehmernom analize kompozitnyh cilindricheskih panelej i obolochek // Vestnik KrasGAU. - 2015.
- № 2. - S. 75-85.
15. Matveev A.D. Raschet trehmernyh kompozitnyh balok slozhnoj formy s primeneniem dvuhsetochnyh konechnyh jelementov // Vestnik KrasGAU. - 2015. - № 8. - S. 92-98.
16. Matveev A.D. Raschet kompozitnyh plastin i balok s uchetom ih struktury s primeneniem slozhnyh mnogosetochnyh konechnyh jelementov // Vestnik Kras-GAU. - 2015. - № 9.
- S. 100-107.
17. Matveev A.D. Metod mnogosetochnyh konechnyh jelementov v raschetah kompozitnyh plastin i balok // Vestnik KrasGAU. - 2016. -№ 12. - S. 93-100.
18. Matveev A.D. Metod mnogosetochnyh konechnyh jelementov v raschetah trehmernyh odnorodnyh i kompozitnyh tel // Uchen. zapiski Kazan. un-ta. Serija: Fiz.-mat. nauki. - 2016. -T. 158. - № 4. - S. 530-543.
19. Matveev A.D. Multigrid finite element method in stress of three-dimensionalelastic bodies of heterogeneous structure // IOP Conf. Ser.: Mater. Sci. Eng. - 2016. - V. 158, № 1. - Art. 012067.
- P. 1-9.
20. Matveev A.D., Grishanov A.N. Trehmernye kompozitnye mnogosetochnye konechnye jele-menty obolochechnogo tipa // Izv. AltGU. Ser. fiziko-matematicheskie nauki. - 2017. - № 4. -C. 120-125.
