Журнал вычислительной математики и математической физики, 2022, T. 62, № 4, стр. 659-676

Исследование и применение метода декомпозиции области для моделирования тепловыделяющего элемента

М. П. Галанин^1, *, А. С. Родин^1, **

¹ ИПМ им. М.В. Келдыша РАН
125047 Москва, Миусская пл., 4, Россия

^* E-mail: galan@keldysh.ru
^** E-mail: rals@bk.ru

Поступила в редакцию 12.01.2021
После доработки 15.11.2021
Принята к публикации 16.12.2021

EDN: AUBCKK
DOI: 10.31857/S004446692204007X

Аннотация

Предложены вычислительные алгоритмы, основанные на применении двухуровневого аддитивного метода Шварца (вариант метода декомпозиции области с перекрытием) и метода Неймана–Дирихле (вариант метода декомпозиции области без перекрытия) для решения задачи контактного взаимодействия системы тел. Представлены результаты применения алгоритмов для численного моделирования напряженно-деформированного состояния участка тепловыделяющего элемента, включающего в себя от 2 до 100 топливных таблеток. Использовано осесимметричное термоупругое приближение. Выполнено исследование зависимости количества итераций алгоритма, требуемого для достижения заданного уровня точности, от количества тел в системе, шага расчетной сетки и граничных условий, поставленных на контактных поверхностях. Библ. 18. Фиг. 8. Табл. 6.

Ключевые слова: контактная задача, метод декомпозиции области, двухуровневый аддитивный метод Шварца, метод Неймана–Дирихле, тепловыделяющий элемент.

ВВЕДЕНИЕ

При расчете напряженно-деформированного состояния (НДС) реалистичной конструкции часто нужно учитывать контактное взаимодействие ее составных частей. Если отдельных деталей достаточно много, а итоговая конфигурация контактных поверхностей является неизвестной, то стандартные алгоритмы, хорошо работающие для случая контакта двух тел, могут оказаться недостаточно точными и быстрыми. Для эффективного решения общей мультиконтактной задачи могут потребоваться новые итерационные алгоритмы, учитывающие особенности конструкции. В данной работе предложены подобные алгоритмы, основанные на использовании вариантов метода декомпозиции области (МДО) и многосеточного метода.

МДО применяется при решении достаточно широкого класса краевых задач для дифференциальных уравнений в частных производных эллиптического и параболического типа, в том числе для статических задач механики деформируемого твердого тела (см. [1], [2]). Суть метода заключается в сведении численного решения задачи в большой области к решению ряда локальных задач в подобластях меньшего размера в рамках дополнительного (внешнего) итерационного процесса. Для учета контактного взаимодействия тел также используют ряд модернизированных вариантов МДО (см. [3]–[7]). Поскольку МДО является итерационным, то эффективность его применения во многом определяется скоростью сходимости алгоритма. В [1] на примере решения уравнения Пуассона показано, что некоторые варианты МДО обладают привлекательным свойством: количество внешних итераций, необходимое для достижения требуемой точности, практически не зависит от размерности глобальной системы уравнений. Для одноуровневых вариантов МДО (с использованием одной сетки) количество итераций кратно возрастает при увеличении количества подобластей. Для устранения данного недостатка разработаны двухуровневые варианты МДО, в которых кроме основной сетки используется дополнительная грубая сетка (аналог многосеточных методов из [8]–[10]). В [1] показано, что возможно построение двухуровневого метода, для которого количество внешних итераций, необходимое для достижения требуемой точности, практически не зависит от количества вводимых подобластей. Для ряда вариантов МДО решение локальных задач в каждой подобласти может выполняться параллельно. Подобная возможность является критически важной для задач большой размерности.

Разрабатываемые алгоритмы применены для решения задачи контактного взаимодействия тел в тепловыделяющем элементе (“твэл”). В “твэле” обычно находятся сотни топливных таблеток, каждая из которых вступает в контакт с соседними таблетками и с оболочкой, причем конфигурация и положение контактных поверхностей могут значительно меняться в процессе функционирования устройства. В общем случае для моделирования работы “твэла” в режиме нормальной эксплуатации нужно решать трехмерную квазистационарную задачу с учетом ряда нелинейных эффектов: возникновения пластических и радиационных деформаций, деформаций ползучести, растрескивания таблеток и т.д. В данной работе рассмотрена осесимметричная термоупругая статическая постановка задачи.

Целью работы является построение алгоритма на основе двухуровневого аддитивного метода Шварца, для которого количество внешних итераций, требуемое для достижения заданной точности, меняется незначительным образом при измельчении шага сетки и при увеличении количества топливных таблеток в моделируемой области. Алгоритм должен обладать возможностью распараллеливания (для локальных задач) и должен естественным образом обобщаться на трехмерный нестационарный случай с учетом упомянутых нелинейных эффектов.

1. МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ

Рассмотрим в многомерном пространстве конструкцию из $N$ термоупругих тел, занимающую область $G = \bigcup\nolimits_{i = 1}^N {{{G}_{i}}} $ с границей $\partial G$ ($\bar {G} = G \cup \partial G$). Математическая модель статической задачи механики деформируемого твердого тела в рассматриваемой постановке состоит из следующих соотношений (см. [1]):

– уравнение равновесия

(1)

$L{\mathbf{u}} = - \nabla \cdot {\mathbf{\sigma }}\left( {\mathbf{u}} \right) = {\mathbf{f}}({\mathbf{x}}),\quad {\mathbf{x}} \in G,$

– граничные условия (кинематические и силовые)

(2)

${\mathbf{u}}\left( {\mathbf{x}} \right) \cdot {\mathbf{n}} = 0,\quad {\mathbf{x}} \in \partial {{G}_{D}},$

(3)

${\mathbf{\sigma }}\left( {\mathbf{u}} \right) \cdot {\mathbf{n}} = {\mathbf{p}}\left( x \right),\quad x \in \partial {{G}_{N}},$

– соотношение Коши для тензора деформации

(4)

${\mathbf{\varepsilon }}\left( {\mathbf{u}} \right) = \frac{1}{2}\left( {\nabla {\mathbf{u}} + {{{(\nabla {\mathbf{u}})}}^{{\text{т}}}}} \right),$

– закон Гука для тензора напряжений

(5)

${\mathbf{\sigma }}\left( {\mathbf{u}} \right) = {\mathbf{C}}:\left( {{\mathbf{\varepsilon }}\left( {\mathbf{u}} \right) - {{{\mathbf{\varepsilon }}}^{{th}}}\left( {\mathbf{x}} \right)} \right).$

Здесь $u\left( x \right)$ – вектор перемещения точки, определяемой радиус-вектором $x = {{x}_{i}}{{e}_{i}}$, ${\mathbf{f}}\left( x \right) = {{f}_{i}}\left( x \right){{e}_{i}}$ – вектор массовых сил, $\partial {{G}_{D}}$ – участок границы, на котором заданы кинематические условия (условия Дирихле) для нормальной компоненты вектора перемещения, $\partial {{G}_{N}}$ – участок границы, на котором заданы силовые условия (условия Неймана), $p\left( x \right) = p_{i}^{{}}\left( x \right){{e}_{i}}$ – вектор внешней нагрузки, действующей на поверхности $\partial {{G}_{N}}$, ${\mathbf{C}}$ – тензор коэффициентов упругости, ${{{\mathbf{\varepsilon }}}^{{th}}}$ – тензор температурной деформации.

Различные тела контактируют друг с другом, всего в конструкции есть $P$ пар потенциально контактирующих поверхностей $\left( {\Gamma {{{_{{c,1}}^{t}}}_{{}}},\Gamma _{{c,2}}^{t}} \right)$, t = 1, 2, …, P. Введем обозначения

$\Gamma _{{c,1}}^{{}} = \bigcup\limits_{t = 1}^P {\Gamma _{{c,1}}^{t}} ,\quad \Gamma {{_{{c,2}}^{{}}}_{{}}} = \bigcup\limits_{t = 1}^P {\Gamma _{{c,2}}^{t}} ,\quad {{\Gamma }_{c}} = {{\Gamma }_{{c,1}}} \cup {{\Gamma }_{{c,2}}}.$

При использовании на контактных границах условия прилипания дополнительные граничные условия на $\Gamma _{{c,1}}^{t}$ состоят из следующих соотношений (для границы $\Gamma _{{c,2}}^{t}$ можно записать аналогичные соотношения) (см. [3], [4]):

${{\sigma }_{\tau }}({\mathbf{x}}) = {{\sigma }_{\tau }}({\mathbf{\bar {x}}}),$

(6)

$\begin{gathered} {{\sigma }_{n}}({\mathbf{x}}) = {{\sigma }_{n}}({\mathbf{\bar {x}}}) \leqslant 0, \\ {{\operatorname{u} }_{n}}({\mathbf{x}}) + {{\operatorname{u} }_{n}}({\mathbf{\bar {x}}}) \leqslant {{\delta }_{{0n}}}({\mathbf{x}}), \\ \end{gathered} $

${{\sigma }_{\tau }}({\mathbf{x}})\left[ {{{\operatorname{u} }_{\tau }}({\mathbf{x}}) + {{\operatorname{u} }_{\tau }}({\mathbf{\bar {x}}}) - {{\delta }_{{0\tau }}}({\mathbf{x}})} \right] + {{\sigma }_{n}}({\mathbf{x}})\left[ {{{\operatorname{u} }_{n}}({\mathbf{x}}) + {{\operatorname{u} }_{n}}({\mathbf{\bar {x}}}) - {{\delta }_{{0n}}}({\mathbf{x}})} \right] = 0,$

где ${\mathbf{\bar {x}}} \in \Gamma _{{c,2}}^{t}$ – сходственная (ближайшая) точка для ${\mathbf{x}} \in \Gamma _{{c,1}}^{t}$, ${{\delta }_{{0\tau }}}{\mathbf{(x)}} \geqslant 0$, ${{\delta }_{{0n}}}({\mathbf{x}}) \geqslant 0$ – функции, задающие начальный зазор (участки поверхностей в начальный момент могли не соприкасаться друг с другом), ${\mathbf{n}}$ и ${\mathbf{\tau }}$ – вектор внешней нормали и касательный вектор в рассматриваемой точке, ${{\operatorname{u} }_{\tau }} = {\mathbf{u}} \cdot {\mathbf{\tau }}$, ${{\operatorname{u} }_{n}} = {\mathbf{u}} \cdot {\mathbf{n}}$, ${{\sigma }_{\tau }} = \left( {{\mathbf{\sigma }} \cdot {\mathbf{n}}} \right) \cdot {\mathbf{\tau }}$, ${{\sigma }_{n}} = \left( {{\mathbf{\sigma }} \cdot {\mathbf{n}}} \right) \cdot {\mathbf{n}}$.

Последнее равенство в (6) гарантирует, что реализуется один из двух следующих вариантов:

1) рассматриваемая точка ${\mathbf{x}} \in \Gamma _{{c,1}}^{t}$ находится в контакте, и тогда верны соотношения

${{\operatorname{u} }_{\tau }}({\mathbf{x}}) + {{\operatorname{u} }_{\tau }}({\mathbf{\bar {x}}}) = {{\delta }_{{0\tau }}}({\mathbf{x}}),\quad {{\operatorname{u} }_{n}}({\mathbf{x}}) + {{\operatorname{u} }_{n}}({\mathbf{\bar {x}}}) = {{\delta }_{{0n}}}({\mathbf{x}}),\quad {{\sigma }_{n}}({\mathbf{x}}) < 0;$

2) рассматриваемая точка не находится в контакте, и тогда верны соотношения

${{\operatorname{u} }_{n}}({\mathbf{x}}) + {{\operatorname{u} }_{n}}({\mathbf{\bar {x}}}) < {{\delta }_{{0n}}}({\mathbf{x}}),\quad {{\sigma }_{\tau }}({\mathbf{x}}) = 0,\quad {{\sigma }_{n}}({\mathbf{x}}) = 0.$

Если использовать на контактных границах условия скольжения без трения, то дополнительные условия на $\Gamma _{{c,1}}^{t}$ принимают следующий вид (аналогично для границы $\Gamma _{{c,2}}^{t}$) (см. [4]):

${{\sigma }_{\tau }}({\mathbf{x}}) = 0,$

(7)

${{\sigma }_{n}}({\mathbf{x}})\left[ {{{\operatorname{u} }_{n}}({\mathbf{x}}) + {{\operatorname{u} }_{n}}({\mathbf{\bar {x}}}) - {{\delta }_{{0n}}}({\mathbf{x}})} \right] = 0.$

В дальнейшем для обозначения блока граничных условий (6) будем использовать векторное неравенство ${{{\mathbf{g}}}_{1}}\left( {\mathbf{u}} \right) \leqslant 0$, а для обозначения блока граничных условий (7) – ${{{\mathbf{g}}}_{2}}\left( {\mathbf{u}} \right) \leqslant 0$.

2. ВЫЧИСЛИТЕЛЬНАЯ МОДЕЛЬ

Задача (1)–(6) (или (1)–(5), (7)) является нелинейной, поскольку итоговая конфигурация контактирующих поверхностей заранее неизвестна. Для ее решения обычно применяют различные итерационные численные методы, наиболее популярными из которых являются метод множителей Лагранжа, метод штрафных функций или их комбинация (см. [3], [11]).

2.1. Метод множителей Лагранжа

Применение метода множителей Лагранжа позволяет свести исходную задачу к задаче минимизации функционала (см. [3])

(8)

$\Pi = \frac{1}{2}\int\limits_G {{{\sigma }_{{ij}}}{{\varepsilon }_{{ij}}}dG} - \int\limits_G {{{u}_{i}}{{f}_{i}}dG} - \int\limits_{\partial {{G}_{N}}} {{{u}_{i}}{{p}_{i}}d\Gamma } + \int\limits_{{{\Gamma }_{{c,1}}}} {{{\lambda }_{i}}({\mathbf{x}})\left( {{{\operatorname{u} }_{i}}({\mathbf{x}}) + {{\operatorname{u} }_{i}}({\mathbf{\bar {x}}}) - {{\delta }_{{0i}}}({\mathbf{x}})} \right)d\Gamma } ,$

где ${\mathbf{\lambda }} = {{\left( {\begin{array}{*{20}{c}} {{{\lambda }_{\tau }}}&{{{\lambda }_{n}}} \end{array}} \right)}^{{\text{т}}}}$ – множитель Лагранжа.

Для пространственной дискретизации полученных уравнений используем метод конечных элементов (МКЭ) (см. [11]). Построим в расчетной области сетку из четырехугольных конечных элементов, в $i$-м теле сетка состоит из ${{n}_{i}}$ узлов, общее количество узлов $n$.

После применения стандартных процедур МКЭ задачу (8) можно свести к решению следующей системы линейных уравнений:

$\left[ {\begin{array}{*{20}{c}} K&{{{M}^{{(m)}}}} \\ {{{M}^{{(m)T}}}}&0 \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{{U}^{{(m + 1)}}}} \\ {{{\lambda }^{{(m + 1)}}}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} R \\ 0 \end{array}} \right\},$

(9)

$\begin{gathered} \left[ K \right] = {\text{diag}}\left[ {{{K}_{{11}}},\; \ldots ,\;{{K}_{{NN}}}} \right], \\ \left[ {{{K}_{{ii}}}} \right] = \int\limits_{{{G}_{i}}} {{{{\left[ B \right]}}^{{\text{т}}}}} \left[ D \right]\left[ B \right]dG, \\ \end{gathered} $

$\left\{ R \right\} = \int\limits_{\partial G_{N}^{{}}} {{{{\left[ {{{L}_{S}}} \right]}}^{{\text{т}}}}p(x)} dS + \int\limits_{{{G}_{{}}}} {{{{\left[ {{{L}_{G}}} \right]}}^{{\text{т}}}}} f(x)dG + \int\limits_{{{G}_{{}}}} {{{{\left[ B \right]}}^{{\text{т}}}}} \left[ D \right]\left\{ {\varepsilon _{{}}^{{th}}} \right\}dG.$

Здесь $\left[ {{{K}_{{ii}}}} \right]$ – глобальная матрица жесткости размерности $2{{n}_{i}} \times 2{{n}_{i}}$ (для двумерной задачи) для $i$-го тела, $\left[ B \right]$ – матрица связи вектора деформаций с вектором перемещений, $\left[ D \right]$ – матрица коэффициентов упругости, матрица $\left[ M \right]$ размерности $2n \times l$ возникает в результате аппроксимации интеграла по контактной поверхности, $\left\{ U \right\}$ – глобальный вектор узловых перемещений размерности $2n$, $\left\{ \lambda \right\}$ – глобальный вектор множителей Лагранжа размерности $l$, $\left\{ R \right\}$ – глобальный вектор узловых сил размерности $2n$, учитывающий вклад поверхностных сил, массовых сил и температурных деформаций.

Матрица системы (9) имеет нулевой диагональный блок, поэтому для решения системы уравнений нельзя использовать большую часть стандартных итерационных методов. Для решения используют либо прямые методы, либо специальные итерационные методы, предназначенные для решения седловых задач (например, метод Узавы, модифицированный метод симметричной последовательной верхней релаксации, см. [14], [15]).

Матрица является плохообусловленной, поэтому, если в последнем интеграле в (8) учитывать вклад всех поверхностных элементов, лежащих на потенциально контактной поверхности ${{\Gamma }_{{c,1}}}$, то численное решение может не соответствовать граничным условиям (6) или (7). Для получения более точного решения будем проводить итерационный процесс и на (m + 1)-й итерации учитывать в (8) вклад только тех поверхностных элементов, которые лежат на поверхности $\tilde {\Gamma }_{{c,1}}^{{(m)}}$ (объединение участков поверхностей, которые на текущей итерации считаются находящимися в контакте). В конце итерации выполняется вычисление новой контактной конфигурации $\tilde {\Gamma }_{{c,1}}^{{(m + 1)}}$: поверхностные узлы, в которых возникли растягивающие нормальные напряжения, исключаются из $\tilde {\Gamma }_{{c,1}}^{{(m)}}$ (по силовому критерию). Итерационный процесс (в дальнейшем будем называть его внешним процессом) продолжается до тех пор, пока не устанавливается окончательная конфигурация контактной поверхности.

2.2. Метод Неймана–Дирихле

Постановка задачи (1)–(6) во многом аналогична постановке задачи для одного тела, разбитого на подобласти, при использовании варианта МДО без перекрытия. Различия связаны с тем, что при использовании классического МДО условия равенства векторов перемещений и напряжений заданы на всем объединении фиктивных границ между подобластями (которые заранее известны и не меняются в ходе расчета), а в контактной задаче аналогичные равенства будут выполняться только на тех участках поверхности ${{\tilde {\Gamma }}_{c}}$, которые в итоге будут находиться в контакте. Благодаря подобному сходству, для решения контактных задач (с различными условиями на контактирующих поверхностях) применяют модифицированные варианты МДО без перекрытия (см. [1]). Суть этих методов заключается в том, что в рамках одной итерации для каждого тела на контактных поверхностях ставится либо кинематическое условие (условие Дирихле), либо силовое (условие Неймана). Затем полученные на границах значения перемещений или напряжений корректируются с помощью специальных процедур (при необходимости корректируется и конфигурация находящихся в контакте поверхностей). Итерационный процесс продолжается до тех пор, пока оба условия (Дирихле и Неймана) не будут выполнены с достаточной точностью.

В частности, в методе Неймана–Дирихле на одной половине контактных поверхностей $\tilde {\Gamma }_{{c,1}}^{{(m)}}$ всегда ставится условие Неймана, а на противолежащих им поверхностях $\tilde {\Gamma }_{{c,2}}^{{(m)}}$ всегда ставится условие Дирихле (таким образом, чтобы поле перемещений обеспечивало выполнение кинематической части условий (6) или (7)). Для условий скольжения граничные условия на контактных поверхностях на (m + 1)-й итерации принимают вид

(10)

${\mathbf{\sigma }}\left( {{{{\mathbf{u}}}^{{(m + 1)}}}} \right) \cdot {\mathbf{n}} = {\mathbf{p}}_{с}^{{(m)}}\left( x \right),\quad x \in \Gamma _{{c,1}}^{{}},$

(11)

$\operatorname{u} _{n}^{{(m + 1)}}({\mathbf{x}}) + \operatorname{u} _{n}^{{(m + 1)}}({\mathbf{\bar {x}}}) - {{\delta }_{{0n}}}({\mathbf{x}}) \leqslant 0,\quad x \in \Gamma _{{c,2}}^{{}},$

где ${\mathbf{p}}_{с}^{{(m)}}\left( x \right)$ – вектор распределенных контактных сил (он равен нулю, если $x \notin \tilde {\Gamma }_{{c,1}}^{{(m)}}$), который корректируется таким образом, чтобы в конце итерационного процесса выполнялись силовые условия в (6).

Рассмотрим пространство $H_{0}^{1}\left( {G,\partial {{G}_{D}},\Gamma _{{c,2}}^{{}}} \right)$ функций ${\mathbf{u}}$, для которых выполнены кинематические условия (2) и (11). Введем билинейную форму для функций ${\mathbf{u}},{\mathbf{v}} \in H_{0}^{1}(G,\partial {{G}_{D}},\Gamma _{{c,2}}^{{}})$:

$a\left( {{\mathbf{u}},{\mathbf{v}}} \right) = \int\limits_G {{\mathbf{C}}\,} {\text{:}}\,{\mathbf{\varepsilon }}\left( {\mathbf{u}} \right)\,{\text{:}}\,{\mathbf{\varepsilon }}\left( {\mathbf{v}} \right)dG,$

и линейные функционалы

$g({\mathbf{v}}) = \int\limits_G {{\mathbf{f}} \cdot {\mathbf{v}}} dG + \int\limits_{\partial {{G}_{N}}} {{\mathbf{p}} \cdot {\mathbf{v}}} dS + \int\limits_G {{\mathbf{C}}\,} {\text{:}}\,{{{\mathbf{\varepsilon }}}^{{th}}}\left( {\mathbf{x}} \right)\,{\text{:}}\,{\mathbf{\varepsilon }}\left( {\mathbf{v}} \right)dG,$

$g_{c}^{{(m)}}({\mathbf{v}}) = \int\limits_{\Gamma _{{c,1}}^{{}}} {{\mathbf{p}}_{c}^{{(m)}} \cdot {\mathbf{v}}} dS.$

Тогда слабая постановка задачи (1)–(5), (10), (11) на (m + 1)-й итерации принимает следующий вид (см. [1], [16]):

(12)

$a\left( {{{{\mathbf{u}}}^{{(m + 1)}}},{\mathbf{v}}} \right) = g({\mathbf{v}}) + g_{c}^{{(m)}}({\mathbf{v}})\quad \forall {\mathbf{v}} \in H_{0}^{1}(G,\partial {{G}_{D}},\Gamma _{{c,2}}^{{}}).$

После использования МКЭ для пространственной дискретизации (12) сводится к системе линейных уравнений

(13)

$\left[ {K_{{}}^{{(m)}}} \right]\left\{ {U_{{}}^{{m + 1}}} \right\} = \left\{ R \right\} + \left\{ {R_{c}^{{(m)}}} \right\},\quad m = 0,\;1,\; \ldots \;.$

Опишем структуру матрицы $\left[ {{{K}^{{(m)}}}} \right]$ для условия скольжения без трения. Множество узлов сетки $J$, которые относятся к области $G$, можно представить в виде объединения подмножеств $J = {{J}^{1}} \cup {{J}^{2}} \cup {{J}^{3}} \cup {{J}^{4}}$, где ${{J}^{1}}$ – внутренние узлы и узлы, относящиеся к $\partial {{G}_{N}}$, ${{J}^{2}}$ – узлы, относящиеся к $\tilde {\Gamma }_{{c,1}}^{{(m)}}$, ${{J}^{3}}$ – узлы, относящиеся к $\partial {{G}_{D}}$, ${{J}^{4}}$ – узлы, относящиеся к $\tilde {\Gamma }_{{c,2}}^{{(m)}}$.

Для узлов $j \in {{J}^{1}} \cup {{J}^{2}}$ строки матрицы $\left[ {{{K}^{{(m)}}}} \right]$ соответствуют строкам стандартной матрицы жесткости. Компоненты вектора узловых сил $\left\{ R \right\}$ формируются из вкладов поверхностных сил, массовых сил и температурной деформации. Компоненты вектора узловых контактных сил $\left\{ {R_{c}^{{(m)}}} \right\}$ отличны от нуля только для узлов $j \in {{J}^{2}}$. Для узлов $j \in {{J}^{3}}$ компонента вектора $\left\{ R \right\}$ равна нулю, а строка матрицы $\left[ {{{K}^{{(m)}}}} \right]$, соответствующая бóльшей (по модулю) компоненте вектора нормали ${{{\mathbf{n}}}_{j}} = {{\left( {{{n}_{r}},{{n}_{z}}} \right)}^{{\text{т}}}}$, имеет вид $\left[ {0,\; \ldots ,\;0,{{k}_{{2j - 1}}},{{k}_{{2j}}},0,\; \ldots ,\;0} \right]$, ${{k}_{{2j - 1}}} = {{n}_{r}}$, ${{k}_{{2j}}} = {{n}_{z}}$. Для узлов $j \in {{J}^{4}}$ компонента вектора $\left\{ R \right\}$ равна $\delta _{{0n}}^{j}$, а строка матрицы $\left[ {{{K}^{{(m)}}}} \right]$, соответствующая бóльшей компоненте вектора нормали, имеет вид (узел $j$ находится напротив поверхностного элемента, в который входят узлы с номерами $p$ и $q$) $\left[ {0,\; \ldots ,\;0,{{k}_{{2j - 1}}},{{k}_{{2j}}},0,\; \ldots ,\;0,{{k}_{{2p - 1}}},{{k}_{{2p}}},0\; \ldots ,\;0,{{k}_{{2q - 1}}},{{k}_{{2q}}},0\; \ldots ,\;0} \right]$, ${{k}_{{2p - 1}}} = - {{n}_{r}}{{L}_{{S,p}}}({{{\mathbf{\bar {x}}}}_{j}})$, ${{k}_{{2p}}} = - {{n}_{z}}{{L}_{{S,p}}}({{{\mathbf{\bar {x}}}}_{j}})$, ${{k}_{{2q - 1}}} = - {{n}_{r}}{{L}_{{S,q}}}({{{\mathbf{\bar {x}}}}_{j}})$, ${{k}_{{2q}}} = - {{n}_{z}}{{L}_{{S,q}}}({{{\mathbf{\bar {x}}}}_{j}})$, где ${{L}_{{S,p}}}$ – базисная функция поверхностного элемента для узла $p$.

В конце итерации выполняется вычисление новой контактной конфигурации $\tilde {\Gamma }_{c}^{{(m + 1)}}$: поверхностные узлы, в которых возникли растягивающие нормальные контактные силы, исключаются из $\tilde {\Gamma }_{c}^{{(m)}}$ (по силовому критерию); узлы, которые в результате сдвига перестали находиться напротив противолежащей поверхности, выходят из контакта по кинематическому критерию; а новые узлы, для которых возник захлест относительно противолежащей поверхности, добавляются в список контактных узлов.

Затем вычисляются новые значения векторов распределенных контактных сил $\left\{ {p_{{c,t}}^{{(m + 1)}}} \right\}$ для каждой поверхности $\Gamma _{{c,1}}^{t}$ (для узлов, которые не находятся в контакте, соответствующие компоненты вектора равны нулю):

(14)

$\left\{ {p_{{c,t}}^{{(m + 1)}}} \right\} = \left( {1 - {{\theta }_{t}}} \right)\left\{ {p_{{c,t}}^{{(m)}}} \right\} - {{\theta }_{t}}\left\{ {\bar {p}_{{c,t}}^{{(m + 1)}}} \right\},$

где $\left\{ {\bar {p}_{{c,t}}^{{(m + 1)}}} \right\}$ – вектор значений распределенных контактных сил в сходственных точках, лежащих на поверхности $\Gamma _{{c,2}}^{t}$, ${{\theta }_{t}}$ – итерационный параметр.

После этого определяется новое приближение для вектора узловых контактных сил:

(15)

$\left\{ {R_{c}^{{(m + 1)}}} \right\} = \int\limits_{\tilde {\Gamma }_{{c,1}}^{{(m + 1)}}} {{{{\left[ {{{L}_{S}}} \right]}}^{{\text{т}}}}p_{c}^{{(m + 1)}}(x)} dS.$

Применяется следующая аппроксимации функции распределенных контактных сил:

$p_{c}^{{(m + 1)}}(x) = \sum\limits_j {{{\chi }_{j}}} {{\left\{ {p_{c}^{{(m + 1)}}} \right\}}_{j}},$

где в качестве базисных функций ${{\chi }_{j}}$ применяются кусочно-постоянные в ячейке Дирихле для узла $j$ функции.

В расчетах с несогласованными сетками более точные результаты получаются, если формулу (11) применять не для ${{\left\{ {p_{{c,t}}^{{(m + 1)}}} \right\}}_{j}}$, а для значения контактных сил ${{\left\{ {F_{{c,t}}^{{(m + 1)}}} \right\}}_{j}} = {{S}_{{D,j}}}{{\left\{ {p_{{c,t}}^{{(m + 1)}}} \right\}}_{j}}$ (${{S}_{{D,j}}}$ – площадь ячейки Дирихле для узла $j$):

(16)

$\left\{ {F_{{c,t}}^{{(m + 1)}}} \right\} = \left( {1 - {{\theta }_{t}}} \right)\left\{ {F_{{c,t}}^{{(m)}}} \right\} - {{\theta }_{t}}\left\{ {\bar {F}_{{c,t}}^{{(m + 1)}}} \right\},$

где ${{\left\{ {\bar {F}_{{c,t}}^{{(m + 1)}}} \right\}}_{j}}$ – значение распределенных контактных сил, проинтегрированных по участку поверхности $\Gamma _{{c,2}}^{t}$, расположенному напротив ячейки Дирихле для узла $j$.

После вычисления $\left\{ {F_{c}^{{(m + 1)}}} \right\}$ определяется вектор $\left\{ {p_{c}^{{(m + 1)}}} \right\}$ и по формуле (15) находится вектор $\left\{ {R_{c}^{{(m + 1)}}} \right\}$.

Если достигнутая точность не удовлетворяет заданному критерию, то происходит переход к следующей итерации.

2.3. Аддитивный метод Шварца для решения систем линейных уравнений

Для эффективного решения системы линейных уравнений (13) можно использовать аддитивный метод Шварца (вариант МДО с перекрытием подобластей).

Представим область $G$ в виде объединения $G = \bigcup\nolimits_{k = 1}^M {{{\Omega }_{k}}} $ конечного числа своих подобластей ${{\Omega }_{1}}$, …, ${{\Omega }_{M}}$ с границами $\partial {{\Omega }_{1}}$, …, $\partial {{\Omega }_{M}}$ (${{\bar {\Omega }}_{k}} = {{\Omega }_{k}} \cup \partial {{\Omega }_{k}}$). В области ${{\Omega }_{k}}$ находится ${{m}_{k}}$ узлов сетки.

Для решения (13) применим предобусловленный итерационный метод Ричардсона (его итерации будем называть внутренними) (см. [1])

(17)

${{\left\{ {U_{{}}^{{m + 1}}} \right\}}^{{(l + 1)}}} = {{\left\{ {U_{{}}^{{m + 1}}} \right\}}^{{(l)}}} + \alpha {{\left[ {{{B}_{{as}}}} \right]}^{{ - 1}}}{{\left\{ r \right\}}^{{(l)}}},$

где ${{\left\{ r \right\}}^{{(l)}}} = \left\{ R \right\} + \left\{ {R_{c}^{{(m)}}} \right\} - \left[ {K_{{}}^{{(m)}}} \right]{{\left\{ {U_{{}}^{{m + 1}}} \right\}}^{{(l)}}}$ – вектор невязки решения на внутренней итерации с номером $l$, $\alpha $ – некоторый итерационный параметр, существенно влияющий на сходимость данного итерационного процесса, ${{\left[ {{{B}_{{as}}}} \right]}^{{ - 1}}}$ – предобусловливающая матрица для аддитивного метода Шварца. Она задается равенством

(18)

${{\left[ {{{B}_{{as}}}} \right]}^{{ - 1}}} = \sum\limits_{k = 1}^M {\left[ {{{T}_{k}}} \right]} {{\left[ {{{B}_{k}}} \right]}^{{ - 1}}}{{\left[ {{{T}_{k}}} \right]}^{{\text{т}}}},$

где $\left[ {{{T}_{k}}} \right]$ – матрица размерности $2n \times 2{{m}_{k}}$, у которой $j$-й столбец совпадает с $j$-м столбцом единичной матрицы размерности $2n \times 2n$, $\left[ {{{B}_{k}}} \right]$ – матрица размерности $2{{m}_{k}} \times 2{{m}_{k}}$, которая по сути является блоком матрицы $\left[ {K_{{}}^{{(m)}}} \right]$, который, в свою очередь, соответствует узлам сетки, принадлежащим подобласти ${{\Omega }_{k}}$ (вектор $\left\{ {{{r}_{k}}} \right\} = {{\left[ {{{T}_{k}}} \right]}^{{\text{т}}}}{{\left\{ r \right\}}^{{(l)}}}$ состоит из компонент вектора ${{\left\{ r \right\}}^{{(l)}}}$, соответствующих данным узлам сетки).

Как показано в [1], при увеличении количества вводимых подобластей количество внешних итераций в аддитивном методе Шварца, необходимое для достижения заданного уровня точности, может кратно увеличиваться. Подобная зависимость существенно ухудшает эффективность метода. Для преодоления данного недостатка можно использовать аддитивный метод Шварца с грубосеточной коррекцией, предполагающий, что наряду с базовой расчетной сеткой в области задана грубая сетка из ${{N}_{c}}$ узлов, размеры ячеек которой сопоставимы с размерами введенных подобластей. Узлы грубой сетки могут не совпадать с узлами основной сетки, границы расчетной грубой области могут отличаться от границ основной расчетной области. Двухуровневый метод является некоторым аналогом многосеточного метода (V-цикла) (см. [8]–[10]).

В этом случае предобусловливатель задается формулой (см. [1], [17])

(19)

${{\left[ {{{B}_{{as2}}}} \right]}^{{ - 1}}} = \left[ {{{T}_{0}}} \right]{{\left[ {{{B}_{0}}} \right]}^{{ - 1}}}{{\left[ {{{T}_{0}}} \right]}^{{\text{T}}}} + {{\left[ {{{B}_{{as}}}} \right]}^{{ - 1}}},$

где $\left[ {{{T}_{0}}} \right]$ – матрица размерности $2n \times 2{{N}_{c}}$, у которой столбцы, соответствующие узлу грубой сетки $j$, заполнены значениями базисной функции ${{L}_{{сj}}}$ в точках расположения узлов основной сетки, $\left[ {{{B}_{0}}} \right] = {{\left[ {{{T}_{0}}} \right]}^{{\text{т}}}}\left[ {{{K}^{{(m)}}}} \right]\left[ {{{T}_{0}}} \right]$ – матрица размерности $2{{N}_{c}} \times 2{{N}_{c}}$.

С учетом (18) и (19) формулу (17) можно переписать следующим образом:

(20)

${{\left\{ {U_{{}}^{{m + 1}}} \right\}}^{{(l + 1)}}} = {{\left\{ {U_{{}}^{{m + 1}}} \right\}}^{{(l)}}} + \alpha \sum\limits_{k = 0}^M {\left[ {{{T}_{k}}} \right]} {{\left\{ {{{\omega }_{k}}} \right\}}^{{(l)}}},$

где для нахождения векторов ${{\left\{ {{{\omega }_{k}}} \right\}}^{{(l)}}}$ нужно решить ($M + 1$) систему линейных уравнений

(21)

$\left[ {{{B}_{k}}} \right]{{\left\{ {{{\omega }_{k}}} \right\}}^{{(l)}}} = \left\{ {{{r}_{k}}} \right\},\quad k = 0,\; \ldots ,\;M.$

Локальные задачи (21) можно решать независимо друг от друга.

Кроме предобусловленного метода Ричардсона для решения (10) можно использовать предобусловленные итерационные метода крыловского подпространства. Матрица системы является несимметричной, поэтому применим предобусловленный метод сопряженных невязок. Данный метод задается соотношениями (см. [17], [18])

${{\left\{ r \right\}}^{{(0)}}} = \left\{ R \right\} + \left\{ {R_{c}^{{(m)}}} \right\} - \left[ {K_{{}}^{{(m)}}} \right]\left\{ {U_{{}}^{m}} \right\},\quad {{\left\{ p \right\}}^{{(0)}}} = {{\left\{ r \right\}}^{{(0)}}},$

${{\nu }_{l}} = \frac{{{{{\left\{ r \right\}}}^{{(l)T}}}\left[ {{{K}^{{(m)}}}} \right]{{{\left\{ p \right\}}}^{{(l)}}}}}{{{{{\left( {\left[ {{{K}^{{(m)}}}} \right]{{{\left\{ p \right\}}}^{{(l)}}}} \right)}}^{{\text{т}}}}\left[ {{{K}^{{(m)}}}} \right]{{{\left\{ p \right\}}}^{{(l)}}}}},\quad l = 0,\;1,\; \ldots ,$

(22)

$\begin{gathered} {{\left\{ {U_{{}}^{{m + 1}}} \right\}}^{{(l + 1)}}} = {{\left\{ {U_{{}}^{{m + 1}}} \right\}}^{{(l)}}} + {{\nu }_{l}}{{\left\{ p \right\}}^{{(l)}}}, \\ {{\left\{ r \right\}}^{{(l + 1)}}} = {{\left\{ r \right\}}^{{(l)}}} - {{\nu }_{l}}\left[ {{{K}^{{(m)}}}} \right]{{\left\{ p \right\}}^{{(l)}}}, \\ \end{gathered} $

${{\mu }_{{l,s}}} = - \frac{{{{{\left( {\left[ {{{K}^{{(m)}}}} \right]{{{\left\{ p \right\}}}^{{(s)}}}} \right)}}^{{\text{т}}}}\left[ {{{K}^{{(m)}}}} \right]{{{\left[ {{{B}_{{as2}}}} \right]}}^{{ - 1}}}{{{\left\{ r \right\}}}^{{(l + 1)}}}}}{{{{{\left( {\left[ {{{K}^{{(m)}}}} \right]{{{\left\{ p \right\}}}^{{(s)}}}} \right)}}^{{\text{т}}}}\left[ {{{K}^{{(m)}}}} \right]{{{\left\{ p \right\}}}^{{(s)}}}}},\quad s = 0,\; \ldots ,\;l,$

${{\left\{ p \right\}}^{{(l + 1)}}} = {{\left[ {{{B}_{{as2}}}} \right]}^{{ - 1}}}{{\left\{ r \right\}}^{{(l + 1)}}} + \sum\limits_{s = 0}^l {{{\mu }_{{l,s}}}{{{\left\{ p \right\}}}^{{(s)}}}.} $

В формулах (22) вычисление вектора ${{\left[ {{{B}_{{as2}}}} \right]}^{{ - 1}}}{{\left\{ r \right\}}^{{(l + 1)}}}$ также сводится к решению локальных задач (аналогичных (21)) в подобластях и на грубой сетке.

2.4. Двухуровневый аддитивный метод Шварца

Кроме использования аддитивного метода Шварца для эффективного решения систем линейных уравнений можно рассмотреть его применение в качестве итерационного процесса в соответствующих функциональных пространствах (см. [2], [16]).

Введем обозначения $\partial {{\Omega }_{{N,k}}} = \partial {{G}_{N}} \cap \partial {{\Omega }_{k}}$, $\partial {{\Omega }_{{D,k}}} = \partial {{G}_{D}} \cap \partial {{\Omega }_{k}}$, $\Gamma _{{k,j}}^{{}} = \Gamma _{{c,j}}^{{}} \cap {{\Omega }_{k}}$, $\Gamma _{c}^{k} = \Gamma _{{k,1}}^{{}} \cup \Gamma _{{k,2}}^{{}}$.

Выберем некоторую функцию ${{{\mathbf{u}}}^{{(0)}}}$, удовлетворяющую граничному условию (2). Построим последовательность функций $\left\{ {{{{\mathbf{u}}}^{{(m)}}}} \right\}$. На (m + 1)-й итерации сначала нужно решить локальные контактные задачи (с использованием метода Неймана–Дирихле) в подобластях ${{\Omega }_{k}}$ ($k = 1,\; \ldots ,\;M$):

$L{\mathbf{u}}_{k}^{{(m + 1)}} = {\mathbf{f}}({\mathbf{x}}),\quad {\mathbf{x}} \in {{\Omega }_{k}}{{\backslash }}\Gamma _{c}^{k},$

${\mathbf{\sigma }} \cdot {\mathbf{n}} = {\mathbf{p}}\left( {\mathbf{x}} \right),\quad {\mathbf{x}} \in \partial {{\Omega }_{{N,k}}},$

(23)

${\mathbf{\sigma }}\left( {{\mathbf{u}}_{k}^{{(m + 1)}}} \right) \cdot {\mathbf{n}} = {\mathbf{p}}_{с}^{{(m + 1)}}\left( x \right),\quad {\mathbf{x}} \in \Gamma _{{k,1}}^{{}},$

$\operatorname{u} _{{k,n}}^{{(m + 1)}}({\mathbf{x}}) + \operatorname{u} _{{k,n}}^{{(m + 1)}}({\mathbf{\bar {x}}}) - {{\delta }_{{0n}}}({\mathbf{x}}) \leqslant 0,\quad x \in \Gamma _{{k,2}}^{{}},$

${\mathbf{u}}_{k}^{{(m + 1)}}({\mathbf{x}}) = {{{\mathbf{u}}}^{{(m)}}}({\mathbf{x}}),\quad {\mathbf{x}} \in \bar {G}{{\backslash }}({{\Omega }_{k}} \cup \partial {{\Omega }_{{N,k}}}).$

Выделим в пространстве $H_{0}^{1}(G,\partial {{G}_{D}},\Gamma _{{c,2}}^{{}})$ замкнутые подпространства

${{H}_{k}} = \left\{ {{\mathbf{v}} \in H_{0}^{1}(G,\partial {{G}_{D}},\Gamma _{{c,2}}^{{}}):{\mathbf{v}}({\mathbf{x}}) = {\mathbf{0}},\;{\mathbf{x}} \in \bar {G}{{\backslash }}({{\Omega }_{k}} \cup \partial {{\Omega }_{{N,k}}})} \right\}$ ${{H}_{k}} = \left\{ {{\mathbf{v}} \in H_{0}^{1}(G,\partial {{G}_{D}},\Gamma _{{c,2}}^{{}}):{\mathbf{v}}({\mathbf{x}}) = {\mathbf{0}},\;{\mathbf{x}} \in \bar {G}{{\backslash }}({{\Omega }_{k}} \cup \partial {{\Omega }_{{N,k}}})} \right\}$

(см. [2]). Функция ${\mathbf{w}}_{k}^{{(m + 1)}} = {\mathbf{u}}_{k}^{{(m + 1)}} - {{{\mathbf{u}}}^{{(m)}}}$ принадлежит подпространству ${{H}_{k}}$. Тогда слабая постановка для локальной задачи (23) принимает вид

(24)

$a\left( {{\mathbf{w}}_{k}^{{(m + 1)}},{{{\mathbf{v}}}_{k}}} \right) = g({{{\mathbf{v}}}_{k}}) + g_{c}^{{(m + 1)}}({{{\mathbf{v}}}_{k}}) - a\left( {{\mathbf{u}}_{{}}^{{(m)}},{{{\mathbf{v}}}_{k}}} \right)\quad \forall {{{\mathbf{v}}}_{k}} \in {{H}_{k}}.$

В конце каждой итерации итоговое решение вычисляется по формуле

${{{\mathbf{u}}}^{{(m + 1)}}} = {{{\mathbf{u}}}^{{(m)}}} + \alpha \sum\limits_{k = 1}^M {{\mathbf{w}}_{k}^{{(m + 1)}}.} $

В уравнении (24) в правой части стоит слагаемое с интегралом от распределенных контактных сил ${\mathbf{p}}_{с}^{{(m + 1)}}$, которые не известны. Поэтому для решения нужно проводить дополнительный внутренний итерационный процесс. После использования МКЭ для пространственной дискретизации (24) сводится к следующей линейной системе уравнений на (l + 1)-й внутренней итерации:

(25)

$\left[ {K_{k}^{{(m)}}} \right]{{\left\{ {{{w}_{k}}} \right\}}^{{(l + 1)}}} = \left\{ {{{r}_{k}}} \right\} + {{\left\{ {\Delta {{R}_{{k,c}}}} \right\}}^{{(l)}}},\quad k = 1,\; \ldots ,\;M.$

Здесь $\left[ {K_{k}^{{(m)}}} \right]$ – матрица размерности $2{{m}_{k}} \times 2{{m}_{k}}$, которая является блоком матрицы $\left[ {K_{{}}^{{(m)}}} \right]$, $\left\{ {r_{k}^{{(m)}}} \right\} = \left\{ {{{R}_{k}}} \right\} + \left\{ {R_{{k,c}}^{{(m)}}} \right\} - \left[ {K_{k}^{{(m)}}} \right]\left\{ {U_{k}^{{(m)}}} \right\}$ – вектор невязки на предыдущей внешней итерации, ${{\left\{ {\Delta {{R}_{{k,c}}}} \right\}}^{{(l)}}} = {{\left\{ {R_{{k,c}}^{{(m + 1)}}} \right\}}^{{(l)}}} - \left\{ {R_{{k,c}}^{{(m)}}} \right\}$, где вектор узловых контактных сил ${{\left\{ {R_{{k,c}}^{{(m + 1)}}} \right\}}^{{(l)}}}$определяется по формулам, аналогичным (15), (16).

Так же, как делалось выше, можно вместе с базовой расчетной сеткой использовать дополнительную грубую сетку из ${{N}_{c}}$ узлов. В полученном двухуровневом аддитивном методе Шварца глобальный вектор перемещений вычисляется по формуле

$\left\{ {U_{{}}^{{(m + 1)}}} \right\} = \left\{ {U_{{}}^{{(m)}}} \right\} + \alpha \sum\limits_{k = 0}^M {\left\{ {w_{k}^{{(m + 1)}}} \right\}} .$

Вспомогательная задача на грубой сетке имеет вид

$\left[ {{{K}_{0}}} \right]\left\{ {\Delta U_{0}^{{(m + 1)}}} \right\} = \left\{ {R_{0}^{{(m)}}} \right\}.$

Здесь $\left[ {{{K}_{0}}} \right]$ – матрица размерности $2{{N}_{c}} \times 2{{N}_{c}}$, соответствующая матрице жесткости МКЭ на грубой сетке для эффективных материалов, $\left\{ {\Delta U_{0}^{{(m + 1)}}} \right\}$ – вектор размерности $2{{N}_{c}}$, $\left\{ {w_{0}^{{(m + 1)}}} \right\} = \left[ {{{T}_{0}}} \right]\left\{ {\Delta U_{0}^{{(m + 1)}}} \right\}$, компоненты вектора $\left\{ {R_{0}^{{(m)}}} \right\}$ совпадают с компонентами вектора ${{\left[ {{{T}_{0}}} \right]}^{{\text{т}}}}\left\{ {{{r}^{{(m)}}}} \right\}$ за исключением узлов, в которых заданы кинематические условия.

3. РЕЗУЛЬТАТЫ РАСЧЕТОВ С ПРИМЕНЕНИЕМ МЕТОДА ДЕКОМПОЗИЦИИ ОБЛАСТИ

В описанных ниже расчетах итерационный процесс продолжается до тех пор, пока не выполняется соотношение

${{E}_{r}} = {{{{{\left\| {\left\{ {r_{{}}^{{(m + 1)}}} \right\}} \right\|}}_{2}}} \mathord{\left/ {\vphantom {{{{{\left\| {\left\{ {r_{{}}^{{(m + 1)}}} \right\}} \right\|}}_{2}}} {{{{\left\| {\left\{ {r_{{}}^{{(0)}}} \right\}} \right\|}}_{2}}}}} \right. \kern-0em} {{{{\left\| {\left\{ {r_{{}}^{{(0)}}} \right\}} \right\|}}_{2}}}} \leqslant {{\varepsilon }_{0}}.$

Для расчетов, если не указано иное, использовались конечные элементы первого порядка на четырехугольной сетке (см. [11]), выбирались следующие значения входных параметров: $\left\{ {{{U}^{0}}} \right\} = 0$, ${{\varepsilon }_{0}} = {{10}^{{ - 6}}}$, $\alpha $ = 0.5.

Рассмотренные алгоритмы протестированы на решении осесимметричной задачи о контактном взаимодействии $N$ однотипных участков упругой трубы ($N$ = 2, 4, 8), поставленных в направлении оси Oz друг на друга в виде столба, под действием постоянных давлений: ${{p}_{1}}$, приложенного к верхнему торцу верхнего тела, и ${{p}_{2}}$, приложенного к внешним боковым поверхностям тел.

В такой постановке контактная задача эквивалентна задаче для единого тела, для которого известны аналитические решения для компонент тензора напряжений (см. [13]). Конфигурация контактных поверхностей не изменяется, задача фактически является линейной, поэтому “mortar-метод” (вариант метода множителей Лагранжа) дает корректное решение за одну итерацию. Метод Неймана–Дирихле, если в качестве начального приближения для контактного давления на всех нижних торцах выбрать величину ${{p}_{1}}$, тоже сходится за одну итерацию; если выбрать нулевое приближение, то при увеличении количества тел в $M$ раз, количество требуемых итераций будет увеличиваться приблизительно в ${{M}^{{1/2}}}$ раз.

Для двухуровневого аддитивного метода Шварца расчетная область, соответствующая продольному сечению рассматриваемой конструкции, разбивалась на $M = N + 1$ подобластей. Проведенные расчеты показали, что количество итераций, требуемое для достижения заданной точности, незначительно отличалось для различного количества тел и различных значений шага сетки (рассмотрены три последовательно сгущающиеся сетки).

Рассмотрим задачу, позволяющую моделировать ряд термомеханических эффектов, происходящих в “твэле” (см. [7], [15]). В конструкцию “твэла” входят цилиндрическая оболочка и столб из топливных таблеток (их количество достигает нескольких сотен). Оболочка имеет следующие характерные размеры: длина – 3–4 м, диаметр – 1 см, толщина – 1 мм. Продольное сечение топливной таблетки, которая является цилиндром с центральным отверстием и фасками на обоих торцах, приведено на фиг. 1а. Расчетная область для участка “твэла” с двумя таблетками и оболочкой показана на фиг. 1б. В ходе работы «твэла» происходит нагрев таблеток до температуры около 1600 K. Контактное взаимодействие таблеток друг с другом и оболочкой оказывает существенное влияние на НДС “твэла”.

Фиг. 1.

(а) – Продольное сечение топливной таблетки с фасками и центральным отверстием; (б) – cхема расчетной области с разбиением на подобласти (для случая N = 2, N = 3).

Выполним численное моделирование выхода участка «твэла» на номинальную мощность в осесимметричном термоупругом приближении. Свойства материала таблетки соответствуют диоксиду урана, свойства материала оболочки соответствуют сплаву циркония, теплофизические характеристики нелинейно зависят от температуры.

Расчетная область включала в себя $N$ таблеток и участок оболочки. Сначала во всей области решалась динамическая температурная задача с заданной постоянной (по пространству) плотностью тепловыделения в таблетках, температура всей оболочки считалась равной 623 K, на боковых поверхностях таблеток и внутренней поверхности оболочки ставилось граничное условие третьего рода.

После выхода температуры на установившийся режим вычисленное температурное поле применялось для решения механической задачи со следующими граничными условиями: нижние торцы первой таблетки и оболочки закреплены в направлении оси Oz, на верхнем торце верхней таблетки задано давление${{p}_{1}}$ = 50 МПа, на внешней поверхности оболочки задано давление ${{p}_{2}}$ = = 10 МПа. На поверхности контактирующих тел поставлено условие скольжения без трения.

Отметим следующие характерные особенности рассмотренной мультиконтактной задачи. Начальный зазор для всех контактных пар равен нулю. Для каждой контактной пары таблетка/таблетка в результате нагрева происходит выход из контакта значительного количества узлов сетки, расположенных ближе к фаске. Из-за большой протяженности конструкции вдоль оси Oz топливные таблетки, особенно в верхней части столба, смещались на значительное расстояние относительно своего начального положения (и относительно оболочки), поэтому поверхностные сетки в контактных парах таблетка/оболочка всегда не согласованы. Например, для 100 таблеток итоговое смещение сопоставимо с высотой одной таблетки.

В табл. 1 приведены количество итераций и время расчетов с использованием “mortar-метода”, в которых количество таблеток менялось от 2 до 50. Для интерполяции перемещений и множителей Лагранжа использовались КЭ второго порядка, в одной таблетке сетка состояла из 10 ячеек в радиальном и 10 ячеек в осевом направлениях (${{h}_{{1r}}} = 0.308$, ${{h}_{{1z}}} = 1.0$). Для решения системы линейных уравнений (9) применена процедура SparseLU из библиотеки Eigen, предназначенная для работы с разреженными матрицами. Итерационный процесс продолжался до тех пор, пока не устанавливалась окончательная конфигурация контактных поверхностей. Поскольку на каждой итерации матрица менялась, то на каждой итерации LU-разложение делалось заново.

Таблица 1.

Количество итераций и время расчета для “mortar-метода”

$N$	2	5	10	25	50
Количество итераций	12	37	59	140	236
Время, с	2	20	181	3274	36 756

Из табл. 1 видно, что при увеличении количества таблеток количество итераций возрастает по зависимости, близкой к линейной ($O(N)$), а время расчета увеличивается по зависимости, близкой к кубической ($O({{N}^{3}})$). Поэтому данный алгоритм не является эффективным при моделировании большого количества таблеток.

В [15] для решения уравнения (9) применен модифицированный метод симметричной последовательной верхней релаксации (MSSOR) с предобусловливателем, который позволил существенно сократить время счета. Выполнено сравнение результатов расчета участка “твэла” из 10 таблеток, полученных с помощью “mortar-метода” и двухуровневого аддитивного метода Шварца.

В методе Неймана–Дирихле граничные условия на контактных поверхностях поставлены следующим образом: в контактных парах таблетка/таблетка на торце нижней таблетки задано силовое условие (10), а на торце верхней таблетки – кинематическое (11). В контактных парах таблетка/оболочка на внешней боковой поверхности таблетки задано силовое условие, а на внутренней боковой поверхности оболочки – кинематическое. В качестве начального приближения на первой итерации выбрано нулевое контактное давление. Для всех контактных пар итерационный параметр ${{\theta }_{j}}$ = 0.5.

Для решения на каждой итерации системы линейных уравнений (13) могут быть использованы разные методы. Сначала рассмотрим последовательный алгоритм, учитывающий, что при таком выборе граничных условий система уравнений по сути принимает блочный вид: в первую очередь решается блок системы, относящийся к первой таблетке (он не зависит от переменных, соответствующих другим телам). Затем решается блок системы, относящийся ко второй таблетке. Строки матрицы, соответствующие узлам нижнего торца, в которых заданы кинематические условия, переписываются следующим образом: все величины, относящиеся к первой таблетке (они уже известны), переносятся в правую часть локальной системы уравнений. Аналогичным образом происходит последовательное вычисление всех блоков системы, соответствующих таблеткам. В конце выполняется решение блока системы уравнений для оболочки.

В использованном варианте метода Неймана–Дирихле после установления итоговой конфигурации контактных поверхностей для узлов, относящихся к поверхностям $\Gamma _{{c,2}}^{t}$, всегда выполнена кинематическая часть граничных условий (6) или (7) (последнее неравенство и последнее равенство). Для корректности решения всей контактной задачи нужно выполнение силовых условий (первые два равенства в (6) или (7)) для узлов, относящихся к поверхностям $\Gamma _{{c,1}}^{t}$. Для оценки текущей погрешности для поверхности $\Gamma _{{c,1}}^{t}$ будем использовать следующую величину (для условия скольжения):

$E_{{F,t}}^{{(m + 1)}} = \sqrt {\frac{{\sum\limits_{j = 1}^{{{l}_{t}}} {{{{\left[ {\left( {{{{\left\{ {F_{{c,t}}^{{(m + 1)}}} \right\}}}_{{2j - 1}}} + {{{\left\{ {\bar {F}_{{c,t}}^{{(m + 1)}}} \right\}}}_{{2j - 1}}}} \right){{n}_{{r,j}}} + \left( {{{{\left\{ {F_{{c,t}}^{{(m + 1)}}} \right\}}}_{{2j}}} + {{{\left\{ {\bar {F}_{{c,t}}^{{(m + 1)}}} \right\}}}_{{2j}}}} \right){{n}_{{z,j}}}} \right]}}^{2}}} }}{{\sum\limits_{j = 1}^{{{l}_{t}}} {{{{\left[ {{{{\left\{ {F_{{c,t}}^{{(m + 1)}}} \right\}}}_{{2j - 1}}}{{n}_{{r,j}}} + {{{\left\{ {F_{{c,t}}^{{(m + 1)}}} \right\}}}_{{2j}}}{{n}_{{z,j}}}} \right]}}^{2}}} }}} ,$

где $\left\{ {F_{{c,t}}^{{(m + 1)}}} \right\}$ – вектор контактных сил в узлах $\Gamma _{{c,1}}^{t}$ (распределенные силы, умноженные на площади ячеек Дирихле) и $\left\{ {\bar {F}_{{c,t}}^{{(m + 1)}}} \right\}$ – вектор контактных сил, вычисленных на участках поверхности $\Gamma _{{c,2}}^{t}$, расположенных напротив соответствующих ячеек Дирихле. В дальнейшем будем приводить значение величины $E_{F}^{c} = \mathop {\max }\limits_t E_{{F,t}}^{{(m + 1)}}$ – максимум среди относительных погрешностей вектора контактных сил по всем парам таблетка/оболочка.

В табл. 2 для метода Неймана–Дирихле приведены время расчета ${\text{iter}}\_{\text{out}}$, количество итераций, в течение которого происходит выход из контакта узлов в парах таблетка/таблетка (по силовому критерию), а также количество итераций, необходимое для достижения различных значений погрешности невязки ${{E}_{r}}$ и погрешности $E_{F}^{c}$ (приведено в скобках, если заданное значение не достигнуто за проведенное количество итераций, то стоит символ “–” ). Расчеты проводились для следующих значений шагов сетки: в таблетке h₁ = 0.25, в оболочке h₂ = 0.165. В процессе расчета в контактных парах таблетка/таблетка из контакта выходило 5/6 от общего количества элементов на поверхности (10 элементов из 12). Количество таблеток варьировалось от 2 до 100.

Таблица 2.

Количество итераций для метода Неймана–Дирихле

$N$	Время, с	${\text{iter\_out}}$	${\text{iter}}$ ${{E}_{r}}$ = $E_{F}^{c}$ = 10^–3	${\text{iter}}$ ${{E}_{r}}$ = $E_{F}^{c}$ = 10^–4	${\text{iter}}$ ${{E}_{r}}$ = $E_{F}^{c}$ = 10^–5	${\text{iter}}$ ${{E}_{r}}$ = $E_{F}^{c}$ = 10^–6
2	8	15	20 (29)	25 (35)	31 (40)	36 (–)
4	18	21	28 (39)	34 (44)	39 (50)	45 (–)
10	55	40	46 (58)	53 (63)	57 (69)	63 (–)
25	235	80	87 (95)	92 (101)	96 (106)	101 (–)
50	799	138	142 (154)	147 (159)	154 (165)	161 (–)
100	3372	248	261 (272)	265 (277)	270 (283)	278 (–)

Как видно из табл. 2, значение ${\text{iter\_out}}$ значительно возрастает при увеличении количества таблеток. При этом после того, как в контактных парах таблетка/таблетка устанавливается окончательная конфигурация, последующее количество итераций, необходимое для достижения фиксированных значений ${{E}_{r}}$ и $E_{F}^{c}$, для разных расчетов меняется незначительным образом. Для большого количества таблеток время расчета увеличивается по зависимости, близкой к квадратичной ($O({{N}^{2}})$).

Поскольку для расчетов с большим количеством таблеток происходит значительное смещение узлов сетки на поверхности оболочки относительно поверхностных элементов таблеток, а в местах контакта углов фасок с оболочкой возникают концентраторы напряжений, то для стабилизации сходимости итерационного процесса в определенный момент расчета выполнялась фиксация положения сходственных точек для узлов, находящихся в контакте, и соответствующих им ячеек Дирихле.

Для большей наглядности на фиг. 2 приведены графики зависимости ${\text{lg}}({{E}_{r}})$ (фиг. 2а) и ${\text{lg}}(E_{F}^{c})$ (фиг. 2б) от номера итерации ${\text{iter}}$ для расчетов с $N$ = 10, 50, 100.

Фиг. 2.

Графики зависимости от номера итерации ${\text{lg}}({{E}_{r}})$ (а) и ${\text{lg}}(E_{F}^{c})$ (б).

Как уже упоминалось, описанный выше алгоритм решения системы уравнений (13) является строго последовательным. Проведены дополнительные расчеты, в которых для решения (13) применялись два варианта алгоритма с предобусловливателем $\left[ {{{B}_{{as2}}}} \right]$, основанном на использовании аддитивного метода Шварца с грубосеточной коррекцией: метод Ричардсона (20) для $\alpha = 0.4$ и метод сопряженных невязок (22).

В аддитивном методе Шварца расчетная область разбита на подобласти без перекрытия ${{\bar {\Omega }}_{1}}$, …, ${{\bar {\Omega }}_{M}}$, $M = N + 1$. На фиг. 1б показана расчетная область с разбиением на подобласти для случая $N = 2$, $M = 3$. Подобласть ${{\bar {\Omega }}_{1}}$ является объединением нижней половины нижней таблетки и соответствующего ей участка оболочки, ${{\bar {\Omega }}_{M}}$ является объединением верхней половины верхней таблетки и соответствующего ей участка оболочки, для других случаев ${{\bar {\Omega }}_{k}}$, $k = 2,\; \ldots ,\;M - 1$, является объединением верхней половины $k$-й таблетки, нижней половины (k + 1)-й таблетки и соответствующего им участка оболочки. Подобласть ${{\Omega }_{k}}$ получается объединением подобласти ${{\bar {\Omega }}_{k}}$ и дополнительных участков соседних подобластей ${{\bar {\Omega }}_{{k - 1}}}$ и ${{\bar {\Omega }}_{{k + 1}}}$ с размерами, задаваемыми выбранным коэффициентом относительного перекрытия (обычно 0.2). Грубая сетка строилась для двух тел: одно тело соответствовало объединению всех таблеток без учета их фасок, второе тело совпадало с оболочкой. В радиальном направлении помещалась одна ячейка, в осевом направлении размер ячеек брался равным половине высоты таблетки.

На фиг. 3 показаны графики изменения десятичного логарифма величины $E_{r}^{{{\text{in}}}} = {{{{{\left\| {{{{\left\{ r \right\}}}^{{({\text{iter}})}}}} \right\|}}_{2}}} \mathord{\left/ {\vphantom {{{{{\left\| {{{{\left\{ r \right\}}}^{{({\text{iter}})}}}} \right\|}}_{2}}} {{{{\left\| {{{{\left\{ r \right\}}}^{{(0)}}}} \right\|}}_{2}}}}} \right. \kern-0em} {{{{\left\| {{{{\left\{ r \right\}}}^{{(0)}}}} \right\|}}_{2}}}}$ в зависимости от номера внутренней итерации ${\text{iter}}$ для первой внешней итерации ($m = 0$) в расчетах с количеством таблеток 2, 10, 25.

Фиг. 3.

Графики зависимости ${\text{lg}}(E_{r}^{{{\text{in}}}})$ от номера итерации: (а) – метод сопряженных невязок, (б) – метод Ричардсона.

Как видно из графиков, благодаря грубосеточной коррекции, количество внутренних итераций, необходимых для решения системы уравнений, полученной для различного количества таблеток, отличается незначительно. При этом локальные задачи для каждой подобласти (и для грубой сетки) можно решать независимо друг от друга. Поэтому данный подход может быть более эффективным при использовании параллельных вычислений.

Рассмотрим применение двухуровневого аддитивного метода Шварца. В этом случае в каждой из введенных подобластей нужно решать нелинейную задачу. Для подобластей ${{\Omega }_{1}}$, ${{\Omega }_{M}}$ решаются задачи контакта двух термоупругих тел (контактная пара таблетка/оболочка), для остальных подобластей – задачи контакта трех термоупругих тел (одна контактная пара таблетка/таблетка и две контактные пары таблетка/оболочка). Все локальные задачи решаются независимо друг от друга. Для всех контактных пар итерационный параметр ${{\theta }_{j}}$ = 0.5. При решении локальных задач число внутренних итераций всегда равнялось двум.

Корректировка конфигурации контактных поверхностей проводилась только в конце внешних итераций, на внутренних итерациях конфигурация считалась неизменной. Также в конце внешней итерации учитывалось изменение участков оболочки, которые входят в подобласти ${{\Omega }_{k}}$ (они всегда должны быть расположены напротив соответствующих участков топливных таблеток).

На грубой сетке для каждого из двух тел отдельно решалась линейная задача упругости, значения модуля Юнга и коэффициента Пуассона для эффективного материала определялись путем усреднения значений соответствующих характеристик в элементах базовой сетки, которые попадали в рассматриваемую ячейку грубой сетки. Нижние торцы тел считались закрепленными по оси. На внешней поверхности первого тела и на внутренней поверхности второго тела ставилось условие Дирихле для радиального перемещения. Значения приращений в узлах грубой сетки получались путем усреднения приращений радиальных перемещений с предыдущей итерации в узлах основной сетки, расположенных на контактной границе таблетка/оболочка и относящихся к соответствующим ячейкам грубой сетки.

В табл. 3 приведены время расчетов и количество внешних итераций (аналогично табл. 2) для различного количества таблеток в расчетной области. Из табл. 3 видно, что при увеличении количества таблеток от 2 до 10 наблюдается увеличение количества итераций, но в дальнейшем процесс стабилизируется: увеличение количества таблеток вплоть до 100 не приводит к существенному изменению количества итераций. Время расчета увеличивается по зависимости, близкой к линейной ($O(N)$). Поэтому, хотя для моделирования нескольких таблеток данный алгоритм требует существенно бóльших затрат, чем «mortar-метод» и метод Неймана–Дирихле, но для расчета большого количества таблеток он может оказаться существенно более эффективным.

Таблица 3.

Количество итераций для двухуровневого метода Шварца

$N$	Время, с	${\text{iter\_out}}$	${\text{iter}}$ ${{E}_{r}}$ = $E_{F}^{c}$ = 10^–3	${\text{iter}}$ ${{E}_{r}}$ = $E_{F}^{c}$ = 10^–4	${\text{iter}}$ ${{E}_{r}}$ = $E_{F}^{c}$ = 10^–5	${\text{iter}}$ ${{E}_{r}}$ = $E_{F}^{c}$ = 10^–6
2	65	22	45 (74)	70 (99)	95 (124)	120 (–)
4	141	26	58 (89)	86 (119)	117 (150)	148 (–)
10	351	28	68 (92)	93 (125)	127 (158)	159 (–)
25	910	28	71 (93)	93 (126)	126 (159)	159 (–)
50	1865	28	78 (91)	91 (124)	123 (158)	155 (–)
100	4127	28	77 (93)	87 (126)	119 (160)	152 (–)

Для большей наглядности на фиг. 4 приведены графики зависимости ${\text{lg}}({{E}_{r}})$ (фиг. 2а) и ${\text{lg}}(E_{F}^{c})$ (фиг. 2б) от номера итерации ${\text{iter}}$ для расчетов с $N$ = 2, 10, 25.

Фиг. 4.

Графики зависимости от номера итерации ${\text{lg}}({{E}_{r}})$ (а) и ${\text{lg}}(E_{F}^{c})$ (б).

На фиг. 5 приведены графики зависимости ${\text{lg}}({{E}_{r}})$ от номера итерации ${\text{iter}}$ в расчетах тремя методами: “mortar-методом” (calc 1), методом Неймана–Дирихле (calc 2) и двухуровневым аддитивным методом Шварца (calc 3). На фиг. 5а представлены графики для 10 таблеток, а фиг. 5б – для 25 таблеток. Для первого метода на каждой итерации решается система уравнений (9). Несмотря на то что для решения используется прямой метод, в силу плохой обусловленности матрицы системы без процедуры корректировки конфигурации контактных поверхностей полученное численное решение имеет весьма низкую точность. Для второго метода на каждой итерации нужно решать систему уравнений (13), это можно сделать, в том числе, с помощью предобусловливателя, использующего аддитивный метод Шварца. Для третьего метода на каждой итерации выполняются по два решения локальной задачи (25) в каждой подобласти и решение задачи на грубой сеткe; все эти задачи могут быть выполнены независимо друг от друга. В расчетах для 10 и 25 таблеток для первых двух методов количество итераций меньше, чем для третьего. Но для большего количества таблеток для первых двух методов требуемое количество итераций продолжит значительно увеличиваться, а для третьего метода останется практически неизменным.

Фиг. 5.

Графики зависимости ${\text{lg}}({{E}_{r}})$ от номера итерации: (а) – $N$ = 10, (б) –$N$ = 25.

В табл. 4 приведено количество внешних итераций двухуровневого аддитивного метода Шварца для случая, когда в расчетной области находятся 10 таблеток и используются сетки с разными шагами.

Таблица 4.

Количество итераций аддитивного метода Шварца для разных шагов сеток

Шаги сетки	${\text{iter\_out}}$	${{E}_{r}}$ = $E_{F}^{c}$ = 10^–3	${{E}_{r}}$ = $E_{F}^{c}$ = 10^–4	${{E}_{r}}$ = $E_{F}^{c}$ = 10^–5	${{E}_{r}}$ = $E_{F}^{c}$ = 10^–6
h₁ = 0.25, h₂ = 0.165	28	68 (92)	93 (125)	127 (158)	159 (–)
h₁ =0.125, h₂ = 0.0825	30	75 (96)	96 (131)	131 (166)	166 (–)
h₁ =0.06125, h₂ = 0.04125	32	79 (99)	100 (136)	137 (173)	173 (–)

Из табл. 4 видно, что уменьшение шагов сетки не приводит к существенному изменению количества внешних итераций.

В табл. 5 приведено количество внешних итераций для 10 таблеток в расчетах, для которых в контактных парах таблетка/таблетка запрещен выход из контакта вне определенной зоны даже при возникновении растягивающих напряжений (по силовому критерию), но конфигурация контактных пар таблетка/оболочка по-прежнему менялась. Таблица еще раз наглядно демонстрирует, что чем больше меняется конфигурация контактных поверхностей, тем больше требуется итераций для достижения требуемой точности.

Таблица 5.

Количество итераций аддитивного метода Шварца для различных конфигураций контакта таблетка/таблетка (10 таблеток)

Доля поверхностных элементов, вышедших из контакта	${\text{iter\_out}}$	${{E}_{r}}$ = $E_{F}^{c}$ = 10^–3	${{E}_{r}}$ = $E_{F}^{c}$ = 10^–4	${{E}_{r}}$ = $E_{F}^{c}$ = 10^–5	${{E}_{r}}$ = $E_{F}^{c}$ = 10^–6
0	0	16 (20)	21 (25)	25 (29)	30 (36)
1/2	10	36 (44)	46 (59)	59 (70)	71 (–)
5/6	28	68 (92)	93 (125)	127 (158)	159 (–)

Проведена также серия расчетов, в которых на всех контактных поверхностях использованы условия прилипания. В этом случае грубая сетка строилась в прямоугольной области, включающей в себя и таблетки, и оболочку. В радиальном направлении помещалась одна ячейка, в осевом направлении размер ячеек брался равным половине высоты таблетки. Разделение расчетной области на подобласти проводилось идентично предыдущим случаям, количество внутренних итераций в локальных задачах не превышало пяти. Значения итерационных параметров в методе Неймана–Дирихле ${{\theta }_{j}}$ = 0.25. В контактных парах таблетка/таблетка происходил выход из контакта примерно половины элементов поверхностной сетки, но расхождение поверхностей было существенно меньше, чем в расчетах с условием скольжения. Конфигурация контактных пар таблетка/оболочка в процессе расчета не меняется. В табл. 6 приведено количество внешних итераций для различного количества таблеток в расчетной области. Расчеты проводились для фиксированных шагов сетки (в таблетке – 0.25, в оболочке – 0.125) и фиксированного коэффициента относительного перекрытия (0.1).

Таблица 6.

Количество итераций аддитивного метода Шварца для различного количества таблеток, условия прилипания

$N$	${\text{iter\_out}}$	${{E}_{r}}$ = $E_{F}^{c}$ = 10^–3	${{E}_{r}}$ = $E_{F}^{c}$ = 10^–4	${{E}_{r}}$ = $E_{F}^{c}$ = 10^–5	${{E}_{r}}$ = $E_{F}^{c}$ = 10^–6
2	8	15 (21)	21 (26)	26 (32)	32 (38)
10	9	15 (21)	21 (27)	27 (33)	33 (39)
25	9	17 (24)	23 (30)	29 (36)	35 (42)
100	9	20 (30)	29 (34)	34 (37)	40 (–)

В качестве иллюстрации полученных результатов приведем некоторые графики для расчетов с 10 таблетками (h₁ = 0.125, h₂ = 0.08) для условия скольжения (графики без маркера) и условия прилипания (графики с маркерами). Координаты и перемещения на графиках нормированы на 1 мм, а напряжения – на 1 МПа. На фиг. 6 приведены графики распределения радиального перемещения, на фиг. 7 – осевого перемещения, а на фиг. 8 – контактного давления в узлах, расположенных на внешней поверхности таблеток. Аналогичные графики для узлов, расположенных на внутренней поверхности оболочки, визуально неотличимы от показанных (за исключением участков между фасками, где контакта не происходит).

Фиг. 6.

Радиальные перемещения в узлах на внешней поверхности таблеток.

Фиг. 7.

Осевые перемещения в узлах на внешней поверхности таблеток.

Фиг. 8.

Контактные давления в узлах на внешней поверхности таблеток.

Из графиков следует, что в расчете с условием прилипания в среднем радиальные перемещения на 10% больше, максимальное осевое перемещение на 30% меньше, а контактное давление на 40% больше, чем соответствующие величины в расчете с условием скольжения без трения. Вблизи углов фасок находятся концентраторы напряжений. Значения давлений в этих зонах могут на порядок превосходить значения давлений в центральной части контактной поверхности, они увеличиваются при измельчении сетки, но интегральная контактная сила, действующая на поверхность, остается практически неизменной.

ЗАКЛЮЧЕНИЕ

Представлены результаты применения различных алгоритмов для численного моделирования термомеханического состояния участка “твэла” с учетом контактного взаимодействия топливных таблеток друг с другом и с оболочкой. Задача решена в термоупругом осесимметричном приближении, количество таблеток $N$ в расчетной области изменялось от 2 до 100, на контактных поверхностях использованы либо условия прилипания, либо условия скольжения без трения. В результате нагрева топливного столба происходит существенное изменение конфигурации контактных поверхностей. Анализ полученных результатов позволил сделать вывод, что для большого количества таблеток при использовании “mortar-метода” в сочетании с прямым методом решения системы линейных уравнений время расчета увеличивается по зависимости, близкой к $O({{N}^{3}})$, в методе Неймана–Дирихле – по зависимости, близкой к $O({{N}^{2}})$, а для двухуровневого аддитивного метода Шварца – по зависимости, близкой к $O(N)$. В то же время для расчетов с небольшим количеством таблеток первые два алгоритма являются гораздо более быстрыми, чем третий. Для рассмотренной задачи при определенном выборе входных параметров (шаг грубой сетки, величина захлеста, значение итерационных параметров) для двухуровневого аддитивного метода Шварца количество внешних итераций, необходимое для достижения заданной точности, незначительно меняется при использовании различных шагов расчетной сетки и различного количества топливных таблеток. Решение локальных задач в каждой подобласти может быть выполнено независимым (в том числе и параллельным) способом. Указанные свойства делают данный алгоритм эффективным инструментом для решения аналогичных мультиконтактных задач с подробными расчетными сетками.

Авторы выражают глубокую благодарность П.С. Аронову за проведенные расчеты контактных задач с использованием mortar-метода.

Список литературы

Toselli A., Widlund O. Domain decomposition methods – algorithms and theory. Berlin-Heidelberg: Springer-Verlag, 2005. 450 p.
Марчук Г.И. Методы вычислительной математики. М.: Наука, 1989. 608 с.
Wriggers P. Computational contact mechanics. Berlin-Heidelberg: Springer-Verlag, 2006. 520 p.
Eck C., Wohlmuth B. Convergence of a contact – Neumann iteration for the solution of two-body contact problems // Math. Model. and Meth. Appl. Sci. 2003. V. 13. № 8. P. 1103–1118.
Bayada G., Sabil J., Sassi T. Neumann-Neumann domain decomposition algorithm for the Signorini problem // Appl. Math. Lett. 2004. V. 17. P. 1153–1159.
Цвик Л.Б. Принцип поочередности в задачах о сопряжении и контакте твердых деформируемых тел // Прикл. механ. 1980. Т. 16. № 1. С. 13–18.
Галанин М.П., Крупкин А.В., Кузнецов В.И., Лукин В.В., Новиков В.В., Родин А.С., Станкевич И.В. Математическое моделирование контактного взаимодействия системы термоупругих тел методом Шварца для многомерного случая // Изв. вузов. Машиностр. 2016. № 12 (681). С. 9–20.
Федоренко Р.П. Введение в вычислительную физику. Долгопрудный: Изд-кий Дом “Интеллект”, 2008. 504 с.
Мартыненко С.И. Универсальная многосеточная технология. М.: ИПМ им. М.В. Келдыша, 2013. 243 с.
Галанин М.П., Cавенков Е.Б. Методы численного анализа математических моделей. М.: Изд-во МГТУ им. Н.Э. Баумана, 2010. 591 с.
Бате К.-Ю. Методы конечных элементов. М.: Физматлит, 2010. 1024 с.
Bayada G., Sabil J., Sassi T. Convergence of a Neumann-Dirichlet algorithm for two-body contact problems with non local Coulomb’s friction law // Math. Model. and Numer. Analys. 2008. V. 42. P. 243–262.
Тимошенко С. П., Гудьер Дж. Теория упругости. М.: Наука, 1975. 576 с.
Быченков Ю.В., Чижонков Е.В. Итерационные методы решения седловых задач. М.: БИНОМ, 2010. 349 с.
Аронов П.С., Галанин М.П., Родин А.С. Численное решение задачи контактного взаимодействия элементов твэла с помощью mortar-метода и метода декомпозиции области // Вестн. МГТУ им. Н.Э. Баумана. Сер. Естествен. науки. 2021. № 3. С. 4–22.
Korneev V.G., Langer U. Dirichlet–Dirichlet domain decomposition methods for elliptic problems. Word Sci. Publ., 2015. 463 p.
Ильин В.П., Гурьева Я.Л. О некоторых параллельных методах и технологии декомпозиции областей // Зап. научн. семинаров ПОМИ. 2014. Т. 428. С. 89–106.
Ильин В.П., Ицкович Е.А. О методах полусопряженных направлений с динамическим предобуславливателем // Сиб. журнал индустр. матем. 2007. Т. 10. № 4. С. 41–54.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал