1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 62, № 7, 2022

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

Нелинейный конечно-объемный метод для задачи переноса–сжатия границы раздела сред на неструктурированных адаптивных сетках

Ю. В. Василевский 1, 2, *, К. М. Терехов 1, 3, **

1 ИВМ РАН
119333 Москва, ул. Губкина, 8, Россия

2 Сеченовский университет
119991 Москва, ул. Трубецкая, 8, стр. 2, Россия

3 МФТИ
141701 Долгопрудный, пер. Институтский, 9, Россия

* E-mail: yuri.vassilevski@gmail.com
** E-mail: terekhov@inm.ras.ru

Поступила в редакцию 08.01.2022
После доработки 08.01.2022
Принята к публикации 11.02.2022

Полный текст (PDF)

Аннотация

Данная статья посвящена нелинейному методу конечных объемов для решения задачи отслеживания границы раздела сред на неструктурированных адаптивных сетках. Мы рассматриваем подход, основанный на объеме жидкости. Положение фронта описывается долей жидкости в каждой расчетной ячейке. Распространение границы раздела сред включает одновременное решение задачи переноса фракции и задачи сжатия границы раздела. Задача сжатия решается для восстановления резкости границы, которая теряется вследствие численной диффузии. Для дискретизации задачи используется нелинейный монотонный метод конечных объемов. Мы применяем метод к неструктурированным сеткам с адаптивным локальным уточнением. Библ. 62. Фиг. 16.

Ключевые слова: неявное слежение за фронтом, объем жидкости, сжатие границы, нелинейный метод конечных объемов, монотонный метод.

1. ВВЕДЕНИЕ

Методы отслеживания границы раздела сред востребованы во многих физических приложениях: в задачах течения жидкости со свободной поверхностью [1]–[3], течения многофазных жидкостей [4], течения, вызванного средней кривизной [5]–[8], отвердения сплавов [9], улучшения изображений [10] и других. Многочисленные существующие коды применяют метод объема жидкости для неявного отслеживания границы раздела сред [11]–[15]. В литературе представлены два основных подхода к уменьшению размытия границы: геометрическая реконструкция [16]–[18] и методы сжатия границы [19]–[23]. Последний подход требует решения обратного параболического уравнения в частных производных. Типичной проблемой метода сжатия является искажение границы, которое становится более выраженным при больших шагах по времени [18]. В ряде работ предлагается адаптация параметров сжатия границы для решения этой проблемы [24]–[26]. Популярным альтернативным подходом является метод функции уровня [27]–[30]. Метод функции уровня требует решения задачи реинициализации и не сохраняет массу жидкости, ряд работ направлен на решение этих вопросов [31]–[34]. Мы рассматриваем подход, основанный на объеме жидкости, в сочетании с методом сжатия границы.

В литературе предложен ряд схем высокого разрешения для задачи переноса объема жидкости [35]–[38], в частности, для задачи переноса–сжатия [39], [23], [40]–[42]. В настоящей работе мы используем нелинейный метод конечных объемов, разработанный ранее для задачи анизотропной диффузии [43]–[47], для решения задачи переноса–сжатия. В сочетании с неявной схемой Эйлера первого порядка для производной по времени метод сохраняет физичными границы рассчитанной доли жидкости и допускает шаги по времени с числом Киранта, превышающим единицу. Метод применим к неструктурированным сеткам общего вида, т.е. к согласованным сеткам с произвольными многогранными ячейками.

Статья организована следующим образом. В разд. 2 мы описываем постановку задачи и нелинейный метод конечных объемов. В разд. 3 мы представляем нашу параллельную технологию для адаптивного сгущения неструктурированных сеток общего вида. В разд. 4 мы рассматриваем численные эксперименты, которые демонстрируют возможности метода и основанной на нем технологии.

2. НЕЛИНЕЙНАЯ КОНЕЧНО-ОБЪЕМНАЯ ДИСКРЕТИЗАЦИЯ

Пусть скалярная функция $\psi $ представляет в каждой точке расчетной области $\Omega $ долю объема, занимаемую жидкостью, а ${\mathbf{u}}$ – заданное векторное поле скорости. Перенос неизвестной функции $\psi $ описывается уравнением

(2.1)
$\begin{gathered} \frac{{\partial \psi }}{{\partial t}} + \operatorname{div} \left( {\psi {\mathbf{u}}} \right) = 0\quad {\text{в}}\quad \Omega , \\ {{{\mathbf{n}}}^{{\text{T}}}}\nabla \psi = 0\quad {\text{на}}\quad \partial \Omega , \\ {{\left. \psi \right|}_{{t = 0}}} = {{\psi }^{0}}\quad {\text{в}}\quad \Omega , \\ \end{gathered} $
где ${\mathbf{n}}$ – внешняя нормаль к $\partial \Omega $. На границе мы рассматриваем однородное граничное условие Неймана.

Численное решение уравнения переноса приводит к потере четкости границы раздела сред из-за численной диффузии. Для того чтобы восстановить границу, уравнение переноса (2.1) дополняется членом сжатия границы раздела сред [26], [19]:

(2.2)
$\frac{{\partial \psi }}{{\partial t}} + \operatorname{div} \left( {\psi {\mathbf{u}} + \alpha \psi \left( {1 - \psi } \right)\frac{{\nabla \psi }}{{\left\| {\nabla \psi } \right\|}}} \right) = 0\quad {\text{в}}\quad \Omega ,$
при тех же начальных и неймановских граничных условиях. Параметр $\alpha $ обычно определяется с помощью $\left\| {\mathbf{u}} \right\|$, он контролирует степень сжатия границы раздела сред. В данной работе мы настраиваем $\alpha $ адаптивно, следуя [26]. В знаменателе мы используем регуляризованную норму Фробениуса.

Замечание 2.1. Задача (2.2) является нелинейной задачей переноса–диффузии: $\frac{{\partial \psi }}{{\partial t}} + \operatorname{div} \left( {\psi {\mathbf{u}} - D\nabla \psi } \right) = 0$ с отрицательным коэффициентом диффузии $D = - \alpha \psi (1 - \psi ){{\left\| {\nabla \psi } \right\|}^{{ - 1}}},$ в результате чего (2.2) соответствует обратной (во времени) параболической задаче [48].

2.1. Задача переноса

Пусть $\Omega $ покрыта согласованной многогранной сеткой, которая может адаптивно перестраиваться на разных временных шагах. Для ячейки $V$ обозначим ее центр через ${{{\mathbf{x}}}_{V}}$ и припишем ему степень свободы.

Сначала мы рассмотрим отдельно задачу переноса (2.1). Применяя теорему Остроградского–Гаусса к интегралу по ячейке $V$, получаем

(2.3)
$\int\limits_V {\left( {\frac{{\partial \psi }}{{\partial t}} + \operatorname{div} \left( {\psi {\mathbf{u}}} \right)} \right)dV} = \int\limits_V {\frac{{\partial \psi }}{{\partial t}}dV} + \oint\limits_{\partial V} {\psi d{{{\mathbf{S}}}^{{\text{T}}}}{\mathbf{u}}} = 0,$
которая аппроксимируется со вторым порядком точности следующим образом:
(2.4)
${{\left. {\left| V \right|\frac{{\partial \psi }}{{\partial t}}} \right|}_{{{{{\mathbf{x}}}_{V}}}}} + \sum\limits_{f \in \partial V} |{\kern 1pt} f{\kern 1pt} |{\kern 1pt} {{\beta }_{f}}{{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}} = 0,$
где $\left| V \right|$ – объем ячейки $V$, ограниченной множеством граней $\partial V$, $|{\kern 1pt} f{\kern 1pt} |$ и ${{{\mathbf{x}}}_{f}}$ – площадь и барицентр грани $f$, а ${{\beta }_{f}}$ – проекция скорости на нормаль грани:
(2.5)
${{\beta }_{f}} = {{\left. {{\mathbf{n}}_{f}^{{\text{T}}}{\mathbf{u}}} \right|}_{{{{{\mathbf{x}}}_{f}}}}} \approx \frac{1}{{|{\kern 1pt} f{\kern 1pt} |}}\int\limits_f {d{{{\mathbf{S}}}^{{\text{T}}}}{\mathbf{u}}} .$
Знак ${{\beta }_{f}}$ зависит от ориентации нормали к грани ${{{\mathbf{n}}}_{f}}$.

Чтобы дискретизировать (2.4), мы должны аппроксимировать производную по времени. Для простоты изложения используем неявную схему первого порядка (неявную схему Эйлера):

(2.6)
${{\left. {\frac{{\partial \psi }}{{\partial t}}} \right|}_{{{{{\mathbf{x}}}_{V}}}}} \approx {{\left. {\frac{{{{\psi }^{{n + 1}}} - {{\psi }^{n}}}}{{\Delta t}}} \right|}_{{{{{\mathbf{x}}}_{V}}}}}.$
На практике мы используем схему Кранка–Николсон.

Вектор невязки ${{R}_{a}}({{\psi }^{{n + 1}}})$ для неявной схемы Эйлера (2.4) в ячейке $V$ принимает вид:

(2.7)
${{\left. {{{R}_{a}}({{\psi }^{{n + 1}}})} \right|}_{V}} = {{\left. {|{\kern 1pt} V{\kern 1pt} |{\kern 1pt} \left( {{{\psi }^{{n + 1}}} - {{\psi }^{n}}} \right)} \right|}_{{{{{\mathbf{x}}}_{V}}}}} + \Delta t\sum\limits_{f \in \partial V} |{\kern 1pt} f{\kern 1pt} |{\kern 1pt} \beta _{f}^{{n + 1}}{\kern 1pt} {{\left. {{{\psi }^{{n + 1}}}} \right|}_{{{{{\mathbf{x}}}_{f}}}}},$
что требует приближения потока $q = {{\left. {\beta _{f}^{{n + 1}}{{\psi }^{{n + 1}}}} \right|}_{{{{{\mathbf{x}}}_{f}}}}}$ на каждой грани $f$.

Вектору невязки ${{R}_{a}}({{\psi }^{{n + 1}}})$ соответствует якобиан

${{J}_{a}}({{\psi }^{{n + 1}}}) = \partial {{R}_{a}}({{\psi }^{{n + 1}}}){\text{/}}\partial {{({{\psi }^{{n + 1}}})}^{{\text{T}}}}.$
При дальнейшем описании свойств якобиана мы опускаем индекс временного уровня и подразумеваем $\psi \equiv {{\psi }^{{n + 1}}}$ и ${{\beta }_{f}} \equiv \beta _{f}^{{n + 1}}$.

В данной работе мы предполагаем, что на каждом временном слое на гранях сетки задано дискретное бездивергентное поле скорости ${\mathbf{u}}$

(2.8)
$0 = \sum\limits_{f \in \partial V} \,{{\left. {{\mathbf{n}}_{f}^{{\text{T}}}{\mathbf{u}}} \right|}_{{{{{\mathbf{x}}}_{f}}}}}\left| f \right| \approx \sum\limits_{f \in \partial V} \,\int\limits_f {d{{{\mathbf{S}}}^{{\text{T}}}}{\mathbf{u}}} = \oint\limits_{\partial V} {d{{{\mathbf{S}}}^{{\text{T}}}}{\mathbf{u}}} {\kern 1pt} \; = \int\limits_V {\operatorname{div} ({\mathbf{u}})dV} .$
Поскольку ${{\beta }_{f}} = {{\left. {{\mathbf{n}}_{f}^{{\text{T}}}{\mathbf{u}}} \right|}_{{{{{\mathbf{x}}}_{f}}}}}$ предполагается заданным, ключевым вопросом в приближении потока $q$ является вычисление ${{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}}$.

Рассмотрим внутреннюю грань $f = {{V}_{1}} \cap {{V}_{2}}$ с нормалью ${{{\mathbf{n}}}_{f}}$, направленной наружу ячейки ${{V}_{1}}$ и внутрь ячейки ${{V}_{2}}$. Пусть ${{\beta }_{f}} > 0$, тогда метод против потока аппроксимирует ${{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}} = {{\psi }_{1}}$ с первым порядком точности (см. фиг. 1).

Фиг. 1.

В методе против потока ${{\psi }_{1}}$ используется для аппроксимации $\psi $ в барицентре грани.

При такой аппроксимации ${{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}}$ в (2.7) строка якобиана, соответствующая переносу в ячейке ${{V}_{1}}$, окруженной соседями ${{V}_{k}}$, собирается следующим образом:

(2.9)
${{\left. {{{J}_{a}}} \right|}_{{{{V}_{1}}}}} \leftarrow \Delta t\left( {\sum\limits_{f \in \partial {{V}_{1}}} \frac{{\left| {{{\beta }_{f}}} \right| + {{\beta }_{f}}}}{2}} \right)\partial {{\psi }_{1}} - \Delta t\sum\limits_{f \in \partial {{V}_{1}} \cap {{V}_{k}}} \frac{{\left| {{{\beta }_{f}}} \right| - {{\beta }_{f}}}}{2}\partial {{\psi }_{k}},$
что приводит к якобиану, который является M-матрицей [49]. Действительно, диагональные элементы якобиана положительны, а внедиагональные элементы отрицательны, и в силу условия (2.8) сумма элементов строки равна нулю.

Приближение второго порядка для $q$ достигается антидиффузионной поправкой

(2.10)
${{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}} \approx {{\psi }_{1}} + {{\left( {{{{\mathbf{x}}}_{f}} - {{{\mathbf{x}}}_{1}}} \right)}^{{\text{T}}}}\nabla \psi .$

Сперва рассмотрим антидиффузионную поправку на одномерной сетке с одинаковыми расстояниями $h{\text{/}}2$ между центрами коллокации и гранями. Используя обозначения из фиг. 2, введем два различных приближения для поправки $\vartheta = {{\left( {{{{\mathbf{x}}}_{f}} - {{{\mathbf{x}}}_{1}}} \right)}^{{\text{T}}}}\nabla \psi $ в точке ${{{\mathbf{x}}}_{f}}$:

(2.11)
${{\vartheta }_{1}} = \frac{{{{\psi }_{1}} - {{\psi }_{0}}}}{2},\quad {{\vartheta }_{2}} = \frac{{{{\psi }_{3}} - {{\psi }_{2}}}}{2},$
где ${{\vartheta }_{1}}$ и ${{\vartheta }_{2}}$ – приближения в ячейках ${{V}_{1}}$ и ${{V}_{2}}$ соответственно. С этими поправками мы получаем два приближения второго порядка ${{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}}$:
(2.12)
${{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}} \approx {{\psi }_{1}} + {{\vartheta }_{1}} = \frac{{3{{\psi }_{1}} - {{\psi }_{0}}}}{2},\quad {{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}} \approx {{\psi }_{1}} + {{\vartheta }_{2}} = {{\psi }_{1}} + \frac{{{{\psi }_{3}} - {{\psi }_{2}}}}{2},$
со следующими вариациями:

(2.13)
$\partial {{\psi }_{1}} + \partial {{\vartheta }_{1}} = \frac{3}{2}\partial {{\psi }_{1}} - \frac{1}{2}\partial {{\psi }_{0}},\quad \partial {{\psi }_{1}} + \partial {{\vartheta }_{2}} = \partial {{\psi }_{1}} + \frac{1}{2}\partial {{\psi }_{3}} - \frac{1}{2}\partial {{\psi }_{2}}.$
Фиг. 2.

Обозначения для одномерной антидиффузионной поправки.

Поправки ${{\vartheta }_{1}}$ и ${{\vartheta }_{2}}$ не меняют значения суммы в строке у матрицы Якоби. Приближения второго порядка ${{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}} \approx {{\psi }_{1}} + {{\vartheta }_{1}}$ в ячейке ${{V}_{1}}$ и ${{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}} \approx {{\psi }_{1}} + {{\vartheta }_{2}}$ в ячейке ${{V}_{2}}$ вносят вклад в матрицу Якоби следующим образом:

(2.14)
$\begin{gathered} {{\left. {{{J}_{a}}} \right|}_{{{{V}_{1}}}}} \leftarrow \Delta t\left| f \right|\left( {\frac{{3{{\beta }_{f}}}}{2}\partial {{\psi }_{1}} - \frac{{{{\beta }_{f}}}}{2}\partial {{\psi }_{0}}} \right), \\ {{\left. {{{J}_{a}}} \right|}_{{{{V}_{2}}}}} \leftarrow \Delta t\left| f \right|\left( {\frac{{{{\beta }_{f}}}}{2}\partial {{\psi }_{2}} - {{\beta }_{f}}\partial {{\psi }_{1}} - \frac{{{{\beta }_{f}}}}{2}\partial {{\psi }_{3}}} \right), \\ \end{gathered} $
что сохраняет свойство М-матрицы. Однако использование разных приближений потока $q$ с разных сторон грани $f$ нарушает закон сохранения. Чтобы сделать поток единственным, мы рассматриваем нелинейную выпуклую комбинацию двух приближений с неотрицательными коэффициентами ${{\mu }_{1}} + {{\mu }_{2}} = 1$:

(2.15)
${{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}} \approx {{\psi }_{1}} + {{\mu }_{1}}{{\vartheta }_{1}} + {{\mu }_{2}}{{\vartheta }_{2}}.$

Воспользуемся весами, известными как ограниченные средние Ван Лира [50]:

(2.16)
${{\mu }_{1}} = \frac{{|{\kern 1pt} {{\vartheta }_{2}}{\kern 1pt} |}}{{|{\kern 1pt} {{\vartheta }_{1}}{\kern 1pt} | + |{\kern 1pt} {{\vartheta }_{2}}{\kern 1pt} |}},\quad {{\mu }_{2}} = \frac{{|{\kern 1pt} {{\vartheta }_{1}}{\kern 1pt} |}}{{|{\kern 1pt} {{\vartheta }_{1}}{\kern 1pt} | + |{\kern 1pt} {{\vartheta }_{2}}{\kern 1pt} |}},$
получим ${{\mu }_{1}}{{\vartheta }_{1}} + {{\mu }_{2}}{{\vartheta }_{2}} = 0$ для ${{\vartheta }_{1}}{{\vartheta }_{2}} < 0$ и ${{\mu }_{1}}{{\vartheta }_{1}} + {{\mu }_{2}}{{\vartheta }_{2}} = 2{{\mu }_{1}}{{\vartheta }_{1}} = 2{{\mu }_{2}}{{\vartheta }_{2}}$ для ${{\vartheta }_{1}}{{\vartheta }_{2}} > 0$. Чтобы избежать деления на ноль, оператор взятия модуля в (2.16) регуляризируем следующим образом ${{\left| x \right|}_{\varepsilon }} = \sqrt {{{x}^{2}} + {{\varepsilon }^{2}}} $, где $\varepsilon $ – малая константа.

Пренебрегая вариациями ${{\mu }_{1}}$ и ${{\mu }_{2}}$ и рассматривая в случае ${{\vartheta }_{1}}{{\vartheta }_{2}} > 0$ приближение ${{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}} \approx {{\psi }_{1}} + 2{{\mu }_{1}}{{\vartheta }_{1}}$ для ${{V}_{1}}$ и ${{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}} \approx {{\psi }_{1}} + 2{{\mu }_{2}}{{\vartheta }_{2}}$ для ${{V}_{2}}$, мы получаем следующие вклады приближений (2.15) в якобиан:

(2.17)
$\begin{gathered} {{\left. {{{J}_{a}}} \right|}_{{{{V}_{1}}}}} \leftarrow \Delta t\left| f \right|\left( {\left( {{{\beta }_{f}} + {{\beta }_{f}}{{\mu }_{1}}} \right)\partial {{\psi }_{1}} - {{\beta }_{f}}{{\mu }_{1}}\partial {{\psi }_{0}}} \right), \\ {{\left. {{{J}_{a}}} \right|}_{{{{V}_{2}}}}} \leftarrow \Delta t\left| f \right|\left( {{{\beta }_{f}}{{\mu }_{2}}\partial {{\psi }_{2}} - {{\beta }_{f}}\partial {{\psi }_{1}} - {{\beta }_{f}}{{\mu }_{2}}\partial {{\psi }_{3}}} \right). \\ \end{gathered} $
Данный прием известен как метод Пикара [51], который сохраняет свойство М-матрицы у якобиана на каждой нелинейной итерации и свойство консервативности при сходимости нелинейных итераций. На практике мы используем аппроксимацию (2.15) в рамках метода Ньютона.

Антидиффузионная поправка легко обобщается на неструктурированные многогранные сетки в $d$ измерениях. Опять же, мы предполагаем, что ${{\beta }_{f}} > 0$ и, таким образом, ${{V}_{1}}$ является ячейкой против потока. Найдем две поправки ${{\vartheta }_{1}}$ и ${{\vartheta }_{2}}$, используя комбинации соседних точек коллокации, см. фиг. 3 для двумерного случая. Определим антидиффузионную поправку ${{\vartheta }_{1}} = {{\left( {{{{\mathbf{x}}}_{f}} - {{{\mathbf{x}}}_{1}}} \right)}^{{\text{T}}}}\nabla \psi $ для грани $f$ в ячейке ${{V}_{1}}$ следующим образом.

Фиг. 3.

Обозначения для двумерной антидиффузионной коррекции на многоугольной сетке.

Пусть множество ${{{\mathbf{\sigma }}}_{1}}$ состоит из ячеек ${{V}_{k}}$ с точкой коллокации ${{{\mathbf{x}}}_{k}}$, имеющих хотя бы один общий узел с ячейкой ${{V}_{1}}$, и примыкающих к ${{V}_{1}}$ граничных граней с нормалью ${{{\mathbf{n}}}_{{{{f}_{k}}}}}$. Для каждого элемента ${{e}_{k}}$ из ${{{\mathbf{\sigma }}}_{1}}$ мы определим вектор ${{{\mathbf{v}}}_{k}}$ и значение ${{\psi }_{k}}$:

$\left( {{{{\mathbf{v}}}_{k}},{{\psi }_{k}}} \right) = \left\{ \begin{gathered} \left( {{{{\mathbf{x}}}_{k}} - {{{\mathbf{x}}}_{1}},{{\psi }_{k}}} \right),\quad {\text{если}}\quad {{e}_{k}} - {\text{ячейка,}} \hfill \\ \left( {{{{\mathbf{n}}}_{{{{f}_{k}}}}},{{\psi }_{1}}} \right),\quad {\text{если}}\quad {{e}_{k}} - {\text{грань}}{\text{.}} \hfill \\ \end{gathered} \right.$
Отметим, что скалярное произведение вектора ${{{\mathbf{v}}}_{k}}$ с градиентом $\nabla \psi $ представляет собой конечную разность в направлении ${{{\mathbf{v}}}_{k}}$, ${\mathbf{v}}_{k}^{{\text{T}}}\nabla \psi = {{\psi }_{k}} - {{\psi }_{1}}$, а в случае, когда ${{e}_{k}}$ – граничная грань, ${\mathbf{v}}_{k}^{{\text{T}}}\nabla \psi = {{\psi }_{1}} - {{\psi }_{1}} = 0$ определяет однородное граничное условие Неймана. В множестве ${{{\mathbf{\sigma }}}_{1}}$ выберем $d$ элементов ${{e}_{{{{i}_{k}}}}}$, $k = 1, \ldots ,d$, обеспечивающих ${{{\mathbf{x}}}_{1}} - {{{\mathbf{x}}}_{f}} = \sum\nolimits_{k = 1}^d \,{{\gamma }_{{{{i}_{k}}}}}{{{\mathbf{v}}}_{{{{i}_{k}}}}}$ с неотрицательными весами ${{\gamma }_{{{{i}_{k}}}}} \geqslant 0$, которые минимизируют выражение
(2.18)
$\sum\limits_{k = 1}^d \,{{\gamma }_{{{{i}_{k}}}}}\left| {\frac{{{{{\mathbf{v}}}_{{{{i}_{k}}}}}}}{{\left| {{{{\mathbf{v}}}_{{{{i}_{k}}}}}} \right|}} - \frac{{{{{\mathbf{x}}}_{1}} - {{{\mathbf{x}}}_{f}}}}{{\left| {{{{\mathbf{x}}}_{1}} - {{{\mathbf{x}}}_{f}}} \right|}}} \right|,$
тогда антидиффузионная поправка ${{\vartheta }_{1}}$ представляется в виде
(2.19)
${{\vartheta }_{1}} = \left( {\sum\limits_{k = 1}^d \,{{\gamma }_{{{{i}_{k}}}}}} \right){{\psi }_{1}} - \sum\limits_{k = 1}^d \,{{\gamma }_{{{{i}_{k}}}}}{{\psi }_{{{{i}_{k}}}}},$
где для граничной грани имеем ${{\psi }_{{{{i}_{k}}}}} = {{\psi }_{1}}$. Вклад поправки ${{\vartheta }_{1}}$ сохраняет строчные суммы в якобиане и его свойство М-матрицы:
(2.20)
$\partial {{\vartheta }_{1}} = \left( {\sum\limits_{k = 1}^d \,{{\gamma }_{{{{i}_{k}}}}}} \right)\partial {{\psi }_{1}} - \sum\limits_{k = 1}^d \,{{\gamma }_{{{{i}_{k}}}}}\partial {{\psi }_{{{{i}_{k}}}}}.$
На каждой граничной грани $f$, прилегающей к ячейке ${{V}_{1}}$, единственное приближение с поправкой ${{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}} \approx {{\psi }_{1}} + {{\vartheta }_{1}}$ также сохраняет свойства якобиана как М-матрицы.

Отметим, что изменение направления аппроксимированного вектора ${{{\mathbf{x}}}_{1}} - {{{\mathbf{x}}}_{f}}$ на ${{{\mathbf{x}}}_{f}} - {{{\mathbf{x}}}_{1}}$ требует изменения знаков коэффициентов в (2.19), (2.20).

Аналогично, используя множество ${{{\mathbf{\sigma }}}_{2}}$ для ячейки ${{V}_{2}}$, выберем в множестве ${{{\mathbf{\sigma }}}_{2}}$ $d$ элементов ${{e}_{{{{i}_{k}}}}}$, $k = 1, \ldots ,d$, обеспечивающих ${{{\mathbf{x}}}_{f}} - {{{\mathbf{x}}}_{1}} = \sum\nolimits_{k = 1}^d \,{{\gamma }_{{{{i}_{k}}}}}{{{\mathbf{v}}}_{{{{i}_{k}}}}}$ с неотрицательными весами ${{\gamma }_{{{{i}_{k}}}}} \geqslant 0$, которые минимизируют выражение

$\sum\limits_{k = 1}^d \,{{\gamma }_{{{{i}_{k}}}}}\left| {\frac{{{{{\mathbf{v}}}_{{{{i}_{k}}}}}}}{{\left| {{{{\mathbf{v}}}_{{{{i}_{k}}}}}} \right|}} - \frac{{{{{\mathbf{x}}}_{f}} - {{{\mathbf{x}}}_{1}}}}{{\left| {{{{\mathbf{x}}}_{f}} - {{{\mathbf{x}}}_{1}}} \right|}}} \right|,$
тогда антидиффузионная поправка в ячейке ${{V}_{2}}$ имеет вид

${{\vartheta }_{2}}\;{\kern 1pt} = \sum\limits_{k = 1}^d \,{{\gamma }_{{{{i}_{k}}}}}{{\psi }_{{{{i}_{k}}}}} - \left( {\sum\limits_{k = 1}^d \,{{\gamma }_{{{{i}_{k}}}}}} \right){{\psi }_{2}}.$

Значение ${{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}}$ аппроксимируется нелинейной комбинацией (2.15) с коэффициентами (2.16).

2.2. Сжатие границы раздела сред

Рассмотрим интеграл по ячейке для уравнения сжатия границы раздела сред:

(2.21)
$\int\limits_V {\operatorname{div} \left( {\alpha \psi (1 - \psi )\frac{{\nabla \psi }}{{\left\| {\nabla \psi } \right\|}}} \right)dV} = \oint\limits_{\partial V} {\alpha \psi (1 - \psi )d{{{\mathbf{S}}}^{{\text{T}}}}\frac{{\nabla \psi }}{{\left\| {\nabla \psi } \right\|}}} .$
Невязка неявной схемы Эйлера для конечного объема $V$ (2.2) имеет вид
(2.22)
${{\left. {R(\psi )} \right|}_{V}} = {{\left. {{{R}_{a}}(\psi )} \right|}_{V}} + \Delta t\sum\limits_{f \in \partial V} \left| f \right|{{\left. {\alpha \psi (1 - \psi ){{{\left\| {\nabla \psi } \right\|}}^{{ - 1}}}{\mathbf{n}}_{f}^{{\text{T}}}\nabla \psi } \right|}_{{{{{\mathbf{x}}}_{f}}}}},$
где ${{\left. {{{R}_{a}}(\psi )} \right|}_{V}}$ определяется в (2.7).

В задаче (2.22) необходимо дискретизировать поток сжатия $c$:

(2.23)
$c = {{\left. {\alpha \psi (1 - \psi ){{{\left\| {\nabla \psi } \right\|}}^{{ - 1}}}{\mathbf{n}}_{f}^{{\text{T}}}\nabla \psi } \right|}_{{{{{\mathbf{x}}}_{f}}}}} = {{\left. {{{\varkappa }_{f}}\psi (1 - \psi )} \right|}_{{{{{\mathbf{x}}}_{f}}}}},$
где ${{\varkappa }_{f}} = {{\left. {\alpha {{{\left\| {\nabla \psi } \right\|}}^{{ - 1}}}{\mathbf{n}}_{f}^{{\text{T}}}\nabla \psi } \right|}_{{{{{\mathbf{x}}}_{f}}}}}$ – нормальная составляющая скорости сжатия.

Учитывая нормальную составляющую ${{\beta }_{f}}$ скорости ${\mathbf{u}}$ на грани $f$, общей для ячеек ${{V}_{1}}$ и ${{V}_{2}}$, определим аппроксимацию параметра $\alpha $ в ячейке ${{V}_{j}}$, $j = 1,2$:

(2.24)
${{\alpha }_{{f,j}}} = \left| {{{\beta }_{f}}} \right|\frac{{\left| {{{{({{{\mathbf{x}}}_{2}} - {{{\mathbf{x}}}_{1}})}}^{{\text{T}}}}\nabla {{\psi }_{j}}} \right| + \Delta t\left| {{\mathbf{u}}_{{f,j}}^{{\text{T}}}\nabla {{\psi }_{j}}} \right|}}{{\left\| {{{{\mathbf{x}}}_{2}} - {{{\mathbf{x}}}_{1}}} \right\|\left\| {\nabla {{\psi }_{j}}} \right\| + \Delta t\left\| {{{{\mathbf{u}}}_{{f,j}}}} \right\|\left\| {\nabla {{\psi }_{j}}} \right\|}}.$
Данный параметр [26] позволяет избежать складок на поверхности границы раздела сред [18] и таким образом улучшает общепринятый выбор $\alpha = \left| {{{\beta }_{f}}} \right|$ [24], [19].

В формуле (2.24) используются вектора ${{{\mathbf{u}}}_{{f,j}}}$ и градиенты $\nabla {{\psi }_{j}}$ в ячейке ${{V}_{j}}$, $j = 1,2$. Вектор скорости ${{{\mathbf{u}}}_{{f,j}}}$ в ячейке ${{V}_{j}}$ восстанавливается из нормальных скоростей граней [52]:

(2.25)
${{{\mathbf{u}}}_{{f,j}}} = {{\beta }_{f}}{{{\mathbf{n}}}_{f}} + \left( {\mathbb{I} - {{{\mathbf{n}}}_{f}}{\mathbf{n}}_{f}^{{\text{T}}}} \right)\frac{1}{{|{\kern 1pt} {{V}_{j}}{\kern 1pt} |}}\sum\limits_{{{f}_{k}} \in \partial {{V}_{j}}} \left| {{{f}_{k}}} \right|{{\beta }_{{{{f}_{k}}}}}\left( {{{{\mathbf{x}}}_{{{{f}_{k}}}}} - {{{\mathbf{x}}}_{j}}} \right).$

Градиент $\nabla {{\psi }_{i}}$ в любой ячейке ${{V}_{i}}$ восстанавливается методом наименьших квадратов: множество ${{{\mathbf{\sigma }}}_{i}}$ дает уравнения ${\mathbf{v}}_{k}^{{\text{T}}}\nabla {{\psi }_{i}} = {{\psi }_{k}} - {{\psi }_{i}}$, которые собираются в систему ${{A}_{i}}\nabla {{\psi }_{i}} = {{{\mathbf{r}}}_{i}}$. Система решается методом наименьших квадратов $\nabla {{\psi }_{i}} = (A_{i}^{{\text{T}}}{{A}_{i}}{{)}^{{ - 1}}}A_{i}^{{\text{T}}}{{{\mathbf{r}}}_{i}}$. В дальнейшем мы пренебрегаем вариацией градиента в якобиане, $\partial \left\| {\nabla {{\psi }_{i}}} \right\| = 0$.

Для восстановления нормальной компоненты скорости сжатия ${{\varkappa }_{{f,1}}}$ на грани $f$ ячейки ${{V}_{1}}$ используем множество ${{{\mathbf{\sigma }}}_{1}}$ и определим подмножество ${{{\mathbf{\tilde {\sigma }}}}_{1}}$ допустимых шаблонов, обеспечивающих $\sum\nolimits_{k = 1}^d \,{{\gamma }_{{{{i}_{k}}}}}{{{\mathbf{v}}}_{{{{i}_{k}}}}} = - {{{\mathbf{n}}}_{f}}$ при ${{\gamma }_{{{{i}_{k}}}}} \geqslant 0$, и вычислим их средневзвешенное значение

(2.26)
${{\varkappa }_{{f,1}}} = {{\alpha }_{{f,1}}}{{\left\| {\nabla {{\psi }_{1}}} \right\|}^{{ - 1}}}\sum\limits_{{{e}_{{{{i}_{{_{k}}}}}}} \in {{{{\mathbf{\tilde {\sigma }}}}}_{1}}} \sum\limits_{k = 1}^d \,{{\gamma }_{{{{i}_{k}}}}}\left( {{{\psi }_{1}} - {{\psi }_{{{{i}_{k}}}}}} \right){\text{/}}\sum\limits_{{{e}_{{{{i}_{{_{k}}}}}}} \in {{{{\mathbf{\tilde {\sigma }}}}}_{1}}} {{\left( {\sum\limits_{k = 1}^d \,{{\gamma }_{{{{i}_{k}}}}}} \right)}^{{ - 2}}}.$
Такое усреднение приводит к сглаживанию аппроксимации градиента [19]. Вычисление ${{\varkappa }_{{f,2}}}$ на грани $f$ ячейки ${{V}_{2}}$ производится аналогично. Поскольку ${{\psi }_{j}} \in [0,1]$, выражение ${{\psi }_{j}}(1 - {{\psi }_{j}}) \in [0,\;1{\text{/}}4]$ и, следовательно, вклады ${{\psi }_{j}}(1 - {{\psi }_{j}})\partial {{\varkappa }_{{f,j}}}$ из ${{\varkappa }_{{f,j}}}$ в якобиан дают нулевую строчную сумму и сохраняют свойство М-матрицы, $j = 1,2$.

На граничной грани $f$, смежной с ячейкой ${{V}_{j}}$, скорость сжатия ${{\varkappa }_{{f,j}}}$ равна нулю из-за граничного условия Неймана.

Остается построить дискретизацию ${{\left. {\psi (1 - \psi )} \right|}_{{{{{\mathbf{x}}}_{f}}}}}$ в сжимающем потоке ${{c}_{{f,j}}} = {{\left. {{{\varkappa }_{{f,j}}}\psi (1 - \psi )} \right|}_{{{{{\mathbf{x}}}_{f}}}}}$ на грани $f = {{V}_{1}} \cap {{V}_{2}}$, для каждой ячейки ${{V}_{j}}$, $j = 1,2$. Как только это будет сделано, мы используем их линейную комбинацию для определения единого сжимающего потока ${{c}_{f}}$ на грани $f$:

(2.27)
${{c}_{f}} = {{\mu }_{1}}{{c}_{{f,1}}} + {{\mu }_{2}}{{c}_{{f,2}}},$
где коэффициенты выбираются аналогично (2.16):

(2.28)
${{\mu }_{1}} = \frac{{|{\kern 1pt} {{c}_{{f,2}}}{\kern 1pt} |}}{{|{\kern 1pt} {{c}_{{f,1}}}{\kern 1pt} | + |{\kern 1pt} {{c}_{{f,2}}}{\kern 1pt} |}},\quad {{\mu }_{2}} = \frac{{|{\kern 1pt} {{c}_{{f,1}}}{\kern 1pt} |}}{{|{\kern 1pt} {{c}_{{f,1}}}{\kern 1pt} | + |{\kern 1pt} {{c}_{{f,2}}}{\kern 1pt} |}}.$

Аппроксимация ${{c}_{{f,j}}}$ основана на стратегии аппроксимации против потока, применяемой к нелинейной немонотонной функции $\psi (1 - \psi )$:

(2.29)
${{c}_{{f,j}}} = \left( {{{\nu }_{{1,j}}}{{c}_{{1,j}}} + {{\nu }_{{2,j}}}{{c}_{{2,j}}} + \varepsilon {{c}_{{1/2,j}}}} \right){\text{/}}\left( {{{\nu }_{{1,j}}} + {{\nu }_{{2,j}}} + \varepsilon } \right),$
где $\varepsilon > 0$ – малая константа, а коэффициенты ${{\nu }_{{i,j}}} \geqslant 0$ выбираются по формулам
(2.30)
$\begin{array}{*{20}{l}} {{{\nu }_{{1,j}}}}&{ = \max \left( {{\kern 1pt} \operatorname{sgn} ({{\varkappa }_{{f,j}}})\left( {1 - 2\left( {{{\psi }_{1}} + {{{\bar {\vartheta }}}_{{1,j}}}} \right)} \right),0} \right),} \\ {{{\nu }_{{2,j}}}}&{ = - \min \left( {{\kern 1pt} sgn({{\varkappa }_{{f,j}}})\left( {1 - 2\left( {{{\psi }_{2}} + {{{\bar {\vartheta }}}_{{2,j}}}} \right)} \right),0} \right),} \end{array}$
антидиффузионные поправки в ячейке ${{V}_{j}}$ равны ${{\bar {\vartheta }}_{{i,j}}} = {{\left( {{{{\mathbf{x}}}_{f}} - {{{\mathbf{x}}}_{i}}} \right)}^{{\text{T}}}}\nabla {{\psi }_{j}}$, $i,j = 1,2$.

Стратегия аппроксимации против потока для ненулевой ${{\varkappa }_{{f,j}}}$ приводит к трем возможным случаям: дискретизация с обеих сторон $f$ допустима (${{\nu }_{{1,j}}}$ и ${{\nu }_{{2,j}}}$ оба положительны), дискретизация с одной из сторон допустима (${{\nu }_{{1,j}}}$ или ${{\nu }_{{2,j}}}$ положительно), и ни одна из дискретизаций не допустима (${{\nu }_{{1,j}}}$ и ${{\nu }_{{2,j}}}$ оба равны 0). Последний случай проиллюстрирован на фиг. 4. Максимум функции (красный круг) находится между двумя приближениями ${{\left. \psi \right|}_{{{{{\mathbf{x}}}_{f}}}}}$ (желтые круги). В этом случае допустимая дискретизация получается при рассмотрении аппроксимации в максимуме ${{\left. {\psi (1 - \psi )} \right|}_{{{{{\mathbf{x}}}_{f}}}}} \approx 1{\text{/}}4$, что приводит к

${{c}_{{1/2,j}}} = \frac{{{{\varkappa }_{{f,j}}}}}{4},\quad \partial {{c}_{{1/2,j}}} = \frac{{\partial {{\varkappa }_{{f,j}}}}}{4}.$
В результате мы добавляем $\varepsilon $-взвешенный член в (2.29).

Фиг. 4.

Два случая, когда оба односторонних приближения недопустимы.

Используя разложение в ряд Тейлора второго порядка в ячейке ${{V}_{j}}$ и $\partial {\kern 1pt} [{{\psi }_{i}}(1 - {{\psi }_{i}})] = [1 - 2{{\psi }_{i}}]\partial {{\psi }_{i}}$, мы получаем

(2.31)
${{\left. {\psi (1 - \psi )} \right|}_{{{{{\mathbf{x}}}_{f}}}}} \approx {{\psi }_{i}}(1 - {{\psi }_{i}}) + (1 - 2{{\psi }_{i}})({{{\mathbf{x}}}_{f}} - {{{\mathbf{x}}}_{i}}{{)}^{{\text{T}}}}\nabla {{\psi }_{j}} = {{\psi }_{i}}(1 - {{\psi }_{i}}) + (1 - 2{{\psi }_{i}}){{\bar {\vartheta }}_{{i,j}}}.$
Поскольку $0 \leqslant \psi (1 - \psi ) \leqslant \frac{1}{4}$, антидиффузионная поправка должна быть ограничена множителем ${{\eta }_{{i,j}}} \in [0,1]$, чтобы удовлетворять

(2.32)
$0 \leqslant {{\psi }_{i}}(1 - {{\psi }_{i}}) + (1 - 2{{\psi }_{i}}){{\eta }_{{i,j}}}{{\bar {\vartheta }}_{{i,j}}} \leqslant \frac{1}{4}.$

Из условия (2.32) получим

(2.33)
${{\eta }_{{i,j}}} = \frac{1}{{(1 - 2{{\psi }_{i}}){{{\bar {\vartheta }}}_{{i,j}}}}}\left\{ \begin{gathered} \frac{1}{4} - {{\psi }_{i}}(1 - {{\psi }_{i}}),\quad (1 - 2{{\psi }_{i}}){{{\bar {\vartheta }}}_{{i,j}}} > 0, \hfill \\ - {{\psi }_{i}}(1 - {{\psi }_{i}}),\quad (1 - 2{{\psi }_{i}}){{{\bar {\vartheta }}}_{{i,j}}} < 0. \hfill \\ \end{gathered} \right.$

Используя (2.23), получаем односторонние приближения ${{c}_{{i,j}}}$ сжимающего потока на грани $f$ ячейки ${{V}_{j}}$

(2.34)
${{c}_{{i,j}}} = {{\varkappa }_{{f,j}}}{{\psi }_{i}}(1 - {{\psi }_{i}}) + {{\varkappa }_{{f,j}}}(1 - 2{{\psi }_{i}}){{\eta }_{{i,j}}}{{\bar {\vartheta }}_{{i,j}}}.$

Вариация данного приближения

(2.35)
$\partial {{c}_{{i,j}}} = \left( {{{\psi }_{i}}(1 - {{\psi }_{i}}) + (1 - 2{{\psi }_{i}}){{\eta }_{{i,j}}}{{{\bar {\vartheta }}}_{{i,j}}}} \right)\partial {{\varkappa }_{{f,j}}} + {{\varkappa }_{{f,j}}}\left( {1 - 2\left( {{{\psi }_{i}} + {{\eta }_{{i,j}}}\bar {\vartheta }_{i}^{j}} \right)} \right)\partial {{\psi }_{i}} + {{\varkappa }_{{f,j}}}(1 - 2{{\psi }_{i}}){{\eta }_{{i,j}}}\partial {{\bar {\vartheta }}_{{i,j}}}$
накладывает ограничения на выбор аппроксимации ${{c}_{{i,j}}}$: последняя допустима, если коэффициент при $\partial {{\psi }_{i}}$ в (2.35) имеет знак ${{( - 1)}^{{i + 1}}}$, $i = 1,2$. Данное обстоятельство приводит к правилам дискретизации против потока (2.29)–(2.30). Поправки ${{\bar {\vartheta }}_{{i,j}}}$ восстанавливаются в ячейке ${{V}_{j}}$ по (2.18)–(2.19) с учетом знака ${{\varkappa }_{{f,j}}}(1 - 2{{\psi }_{i}})$ при $\partial {{\bar {\vartheta }}_{{i,j}}}$.

По построению каждый поток ${{c}_{{f,j}}}$ сохраняет знаки в матрице Якоби, когда вклад сжимающего потока ассемблируется в уравнение соответствующей ячейки ${{V}_{j}}$. Однако свойство нулевой суммы по строке матрицы нарушается, поскольку скорость сжатия не является бездивергентной. В результате метод сжатия границы раздела сред может породить новый экстремум решения в промежутке $[0,1]$ даже для неявной схемы Эйлера.

2.3. Итерационное решение нелинейной алгебраической системы

Для решения нелинейной системы применим метод Ньютона. Пусть ${{\xi }^{l}}$ – приближение к ${{\psi }^{{n + 1}}}$ на $l$-й итерации, в качестве начального приближения возьмем решение с предыдущего шага по времени ${{\xi }^{0}} = {{\psi }^{n}}$. На итерации метода Ньютона найдем решение линейной системы

(2.36)
$J({{\xi }^{l}})\Delta \xi = - R({{\xi }^{l}}),$
где $R({{\xi }^{l}})$ – невязка из (2.22), а $J({{\xi }^{l}}) = \partial R({{\xi }^{l}}){\text{/}}\partial {{({{\xi }^{l}})}^{{\text{T}}}}$ – соответствующий якобиан. Следующее приближение получается по формуле
(2.37)
${{\xi }^{{l + 1}}} = {{\xi }^{l}} + \omega {\kern 1pt} \Delta \xi .$
Параметр $\omega $ ограничивает максимальное изменение ${{\xi }^{l}}$ в (2.37) на любой ячейке ${{V}_{k}}$: $\omega \leqslant 0.3/{\kern 1pt} |{\kern 1pt} \Delta {{\xi }_{k}}{\kern 1pt} |$ [53]. Итерации могут стагнировать, поэтому мы применяем следующую эвристику, уменьшая $\omega $ с ростом числа нелинейных итераций $l$:
(2.38)
$\omega = \min \left( {\omega ,\frac{{1 + \exp \left( { - \frac{5}{2}} \right)}}{{1 + \exp \left( {\frac{l}{6} - \frac{5}{2}} \right)}}} \right).$
Параметр $\omega $ является единственным ограничением при вычислении следующего итерационного приближения в методе Ньютона.

Итерации продолжаются до тех пор, пока не выполнится условие $\left\| {R({{\xi }^{l}})} \right\| \leqslant \min \left( {{{\tau }_{{{\text{abs}}}}},{{\tau }_{{{\text{rel}}}}}\left\| {R({{\psi }^{n}})} \right\|} \right)$, где ${{\tau }_{{{\text{abs}}}}}{{ = 10}^{{ - 9}}}$, ${{\tau }_{{{\text{rel}}}}}{{ = 10}^{{ - 2}}}$. Линейная система (2.36) собирается и решается с использованием функционала платформы INMOST [54]–[56].

3. АДАПТАЦИЯ СЕТКИ

Для динамической адаптации неструктурированных сеток на многопроцессорных системах мы используем инструментарий платформы INMOST (www.inmost.org). Подробное описание алгоритмов доступно по ссылкам [54], [57], [58].

Ячейка многогранной сетки локально измельчается путем введения нового узла в центре ячейки и висячих узлов в центрах граней и ребер, и дальнейшего разбиения ячейки, граней и ребер, как показано на фиг. 5. Сгущение является постепенным: две соседние ячейки могут отличаться не более чем на один уровень измельчения.

Фиг. 5.

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

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

Фиг. 6.

Организация множеств элементов в древовидную структуру (а). Исходная неразделенная сетка из трех ячеек (б). После разбиения ячейки ${{C}_{3}}$ к родительскому (корневому) набору присоединяется новый набор (в). Это множество запоминает все мелкие ячейки, образующие исходную ячейку ${{C}_{3}}$.

При проходе по сетке происходит измельчение или огрубление не более чем на один уровень. В параллельной среде, после каждого прохода, автоматически восстанавливается согласованность сетки между процессорами, а также недостающие или избыточные элементы перекрывающихся слоев [58], [54]. В данной работе мы используем пакет Parmetis для разбиения сетки [59], [60].

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

На подготовительном этапе мы вычисляем градиенты для ${{\psi }^{{n + 1}}}$ и ${{\psi }^{n}}$ с помощью метода наименьших квадратов. На этапе измельчения интерполяция со старой ячейки ${{V}_{0}}$ на новые ячейки ${{V}_{i}}$ вычисляется следующим образом:

(3.1)
${{\psi }_{i}} = {{\psi }_{0}} + \eta {{({{{\mathbf{x}}}_{i}} - {{{\mathbf{x}}}_{0}})}^{{\text{T}}}}\nabla {{\psi }_{0}},\quad \nabla {{\psi }_{i}} = \nabla {{\psi }_{0}},$
где $\eta $ выбирается так, чтобы удовлетворять условию для новых ячеек ${{V}_{i}}$
(3.2)
$\mathop {\min }\limits_{{{V}_{j}} \in {{\mathcal{V}}_{n}}({{V}_{0}})} ({{\psi }_{j}}) \leqslant {{\psi }_{0}} + \eta {{({{{\mathbf{x}}}_{i}} - {{{\mathbf{x}}}_{0}})}^{{\text{T}}}}\nabla {{\psi }_{0}} \leqslant \mathop {\max }\limits_{{{V}_{j}} \in {{\mathcal{V}}_{n}}({{V}_{0}})} ({{\psi }_{j}}),$
здесь $\psi $ – ${{\psi }^{{n + 1}}}$ или ${{\psi }^{n}}$, а ${{\mathcal{V}}_{n}}({{V}_{0}})$ – множество ячеек, имеющих хотя бы один общий узел с ячейкой ${{V}_{0}}$. Такая интерполяция консервативна при условии $\sum\nolimits_i \,|{\kern 1pt} {{V}_{i}}{\kern 1pt} |{\kern 1pt} {{{\mathbf{x}}}_{i}} = \,|{\kern 1pt} {{V}_{0}}{\kern 1pt} |{\kern 1pt} {{{\mathbf{x}}}_{0}}$ и монотонна в силу (3.2). При огрублении ${{V}_{i}}$ в ячейку ${{V}_{0}}$ мы используем простое усреднение:
(3.3)
${{\psi }_{0}} = {{\left| {{{V}_{0}}} \right|}^{{ - 1}}}\sum\limits_i \left| {{{V}_{i}}} \right|{{\psi }_{i}},\quad \nabla {{\psi }_{0}} = {{\left| {{{V}_{0}}} \right|}^{{ - 1}}}\sum\limits_i \left| {{{V}_{i}}} \right|\nabla {{\psi }_{i}},$
которое является монотонным и консервативным.

Критерий измельчения основан на оценке изменения решения на гранях сетки. Пусть $f$ – грань, разделяющая ячейки ${{V}_{1}}$ и ${{V}_{2}}$. Обе ячейки будут помечены для измельчения, если максимальный уровень сгущения не будет превышен, и если

(3.4)
$\frac{{|{\kern 1pt} \psi _{2}^{{n + 1}} - \psi _{1}^{{n + 1}}{\kern 1pt} |}}{{|{\kern 1pt} {{{\mathbf{x}}}_{2}} - {{{\mathbf{x}}}_{1}}{\kern 1pt} |}} > 5.$

Ячейка огрубляется, если ни одна из ее граней не указывает на необходимость измельчения, и если минимальный уровень огрубления не достигнут.

Приближенное решение нестационарной задачи (2.2) методом конечных объемов с адаптацией сетки происходит следующим образом.

Шаг 1. Решим уравнение переноса (2.1) для ${{\tilde {\psi }}^{{n + 1}}}$, используя дискретизацию первого порядка против потока и неявную схему Эйлера за одну линейную итерацию, чтобы предсказать положение границы раздела на следующем временном шаге.

Шаг 2. Локально измельчим сетку, используя изменение ${{\tilde {\psi }}^{{n + 1}}}$ на гранях сетки, сбалансируем сетку между процессорами и интерполируем ${{\tilde {\psi }}^{{n + 1}}}$ и ${{\psi }^{n}}$ на измельченную сетку.

Шаг 3. Выполним несколько итераций Ньютона для задачи переноса (2.1), используя нелинейную дискретизацию второго порядка по пространству в сочетании со схемой Кранка–Николсон, чтобы получить лучшее начальное предположение ${{\hat {\psi }}^{{n + 1}}}$ для следующего шага.

Шаг 4. Решим уравнение переноса–сжатия (2.2) с помощью нелинейной дискретизации второго порядка по пространству в сочетании со схемой Кранка–Николсон для получения ${{\psi }^{{n + 1}}}$, при этом ${{\hat {\psi }}^{{n + 1}}}$ используется как начальное приближение для ньютоновских итераций.

Шаг 5. Выполним огрубление сетки, используя изменение ${{\psi }^{{n + 1}}}$ на гранях сетки, сбалансируем сетку между процессорами и интерполируем ${{\psi }^{{n + 1}}}$ и ${{\psi }^{n}}$ на огрубленную сетку.

Функции ${{\tilde {\psi }}^{{n + 1}}}$, ${{\hat {\psi }}^{{n + 1}}}$ и ${{\psi }^{{n + 1}}}$ обеспечивают приближение положения границы раздела сред по изоповерхности значения $\psi = 0.5$ с разной точностью, наилучшее приближение обеспечивает ${{\psi }^{{n + 1}}}$, наихудшее – ${{\tilde {\psi }}^{{n + 1}}}$.

4. ЧИСЛЕННЫЕ ПРИМЕРЫ

4.1. Тест Залесака

Рассмотрим тест Залесака со сферой [61], вращающейся в поле скорости ${\mathbf{u}} = [u,{v},w{{]}^{{\text{T}}}}$

(4.1)
$u = \pi \left( {\frac{1}{2} - y} \right),\quad {v} = \pi \left( {x - \frac{1}{2}} \right),\quad w = 0,$
заданном в единичном кубе $\Omega {{ = [0,1]}^{3}}$.

Вращающийся объект определяется сферой радиусом $r = 0.15$ с центром в точке ${\mathbf{x}}{{ = [0.5,0.75,0.5]}^{{\text{T}}}}$, с выемкой, образованной областью $[0.45,0.55] \otimes [0.6,0.725] \otimes [0,1]$. Проекция скорости на грань ячейки ${{\beta }_{f}}$ и фракции жидкости в ячейках $\psi $ вычисляются путем деления многоугольников и многогранников на треугольники и тетраэдры, соответственно, и интегрированием средних значений с $7$-м порядком точности. Фракция жидкости $\psi $ равна 1 внутри объекта и 0 вне объекта. Задача интегрируется по времени $t \in [0,2]$, за которое объект выполняет один оборот.

Начальная сетка получена из кубической сетки $16 \times 16 \times 16$ с двумя уровнями измельчения к границе раздела сред, см. фиг. 7. Полученная адаптивная сетка с двумя уровнями уточнения является согласованной, если ее кубические ячейки рассматривать как многогранные. Объект отображается с помощью программы визуализации Paraview [62] через изоповерхность $\psi = 0.5$.

Фиг. 7.

Изоповерхность $\psi = 0.5$ и срединный ($z = 0.5$) срез сетки при $t = 0$, раскрашенный согласно палитре $\psi \in [0,1]$. Аналогичная палитра используется на следующих рисунках.

На фиг. 8 мы демонстрируем влияние использования дискретизации второго порядка для решения задачи переноса (2.1) с шагом по времени $\Delta t = 0.01$ и числом Куранта $K \approx 1$. Точность дискретизации как по времени, так и по пространству очень важна. Однако схема является достаточно диссипативной и сильное размазывание границы раздела приводит к огрублению сетки, что еще больше снижает точность восстановленной границы раздела сред.

Фиг. 8.

Изоповерхность $\psi = 0.5$ и срединный ($z = 0.5$) срез сетки при $t = \{ 0.5,\;1,\;1.5,\;2\} $, окрашенный согласно палитре $\psi $. Решение (2.1) противопотоковым методом первого порядка с неявной схемой Эйлера (UPW1-BE), нелинейным методом конечных объемов второго порядка с неявной схемой Эйлера (UPW2-BE), нелинейный метод конечных объемов второго порядка со схемой Кранка–Николсон (UPW2-CN). Шаг по времени $\Delta t = 0.01$ и $K \approx 1$.

На фиг. 9 мы приводим сравнение нелинейного метода конечных объемов для задачи переноса–сжатия (2.2) с дискретизациями по времени неявной схемой Эйлера и схемой Кранка–Николсон с шагом по времени $\Delta t = 0.01$, $K \approx 1$. Метод сжатия значительно улучшает разрешение границы раздела сред. Схема Кранка–Николсон разрешает форму объекта лучше, чем неявная схема Эйлера, сохраняя при этом монотонное решение.

Фиг. 9.

Изоповерхность $\psi = 0.5$ и срединный ($z = 0.5$) срез сетки при $t = \{ 0.5,\;1,\;1.5,\;2\} $, окрашенный согласно палитре $\psi $. Решение (2.2) методом нелинейных конечных объемов второго порядка с неявной схемой Эйлера (CUPW2-BE), нелинейным методом конечных объемов второго порядка со схемой Кранка–Николсон (CUPW2-CN). Шаг по времени $\Delta t = 0.01$ и $K \approx 1$.

На фиг. 10 показана применимость нелинейного метода конечных объемов для решения (2.2) при более высоких числах Куранта $K \approx 2$ ($\Delta t = 0.02$) и $K \approx 4$ ($\Delta t = 0.04$). Для неявной схемы Эйлера при $K \approx 4$ сжатие границы раздела сред недостаточно сильное, чтобы сохранить резкую границу, но $\psi \in [0,1]$ вплоть до машинной точности. Схема Кранка–Николсон не является монотонной при $K > 1$. Поэтому далее мы используем схему Кранка–Николсон с $K \approx 1$, как обеспечивающую наилучшее разрешение границы раздела сред и все еще монотонное решение.

Фиг. 10.

Изоповерхность $\psi = 0.5$ и срединный ($z = 0.5$) срез сетки при $t = \{ 0.5,\;1,\;1.5,\;2\} $, окрашенный согласно палитре $\psi $. Решение (2.2) нелинейным методом конечных объемов второго порядка с неявной схемой Эйлера (BE) и схемой Кранка–Николсон (CN) с шагом по времени $\Delta t = 0.02$, $K \approx 2$ (вверху) и $\Delta t = 0.04$, $K \approx 4$ (внизу). Темно-зеленый и светло-зеленый цвета соответствуют $\psi < 0$ и $\psi > 1$ соответственно.

Влияние измельчения кубической сетки до уровней $L = 3$ и $L = 4$ показано на фиг. 11. Несмотря на то что при более мелкой сетке граница раздела сред разрешается лучше, искажения границы вследствие сжатия более заметны при таком разрешении. На фиг. 12 мы сравниваем эффект от адаптивного выбора $\alpha $ (2.24), предложенного в [26], и общепринятого выбора $\alpha = \left| {{{\beta }_{f}}} \right|$. Артефакты на поверхности раздела менее выражены благодаря адаптивному выбору (2.24). Этот эффект не так заметен на более грубых сетках ($L = 2$, $L = 3$). Отметим, что мы используем гораздо больший шаг по времени, чем обычные схемы.

Фиг. 11.

Изоповерхность $\psi = 0.5$ и срединный ($z = 0.5$) срез сетки при $t = \{ 0.5,1,1.5,2\} $, окрашенный согласно палитре $\psi $. Решение (2.2) нелинейным методом конечных объемов второго порядка со схемой Кранка–Николсон на кубических сетках (C) с $L = 3$ уровнями уточнения с шагом по времени $\Delta t = 0.005$ (вверху) и $L = 4$ уровнями уточнения с $\Delta t = 0.0025$ (внизу). В обоих случаях $K \approx 1$.

Фиг. 12.

Изоповерхность $\psi = 0.5$ при $t = \{ 0.5,\;1,\;1.5,\;2\} $, посчитанная на кубической сетке с $L = 4$ уровнями измельчения. Адаптивный выбор $\alpha $ по формуле (2.24) (а) и общепринятый выбор $\alpha = \left| {{{\beta }_{f}}} \right|$ (б).

Представленный метод реализован в рамках платформы INMOST и, таким образом, применим для решения на параллельных компьютерах. Распределение кубической сетки с $L = 4$ уровнями измельчения между $16$ процессорами показано на фиг. 13.

Фиг. 13.

Срединный ($z = 0.5$) срез кубической сетки с $L = 4$ уровнями измельчения при $t = \{ 0.5,\;1,\;1.5,\;2\} $, окрашенный индексом процессора, $0,\;1,\; \ldots ,\;15$.

Динамика изменения количества ячеек и количества необходимых итераций Ньютона во время моделирования на кубической сетке с $L = 3$ уровнями измельчения показана на фиг. 14.

Фиг. 14.

Изменение числа ячеек (а) и числа итераций Ньютона (б) для решения (2.2) нелинейным методом конечных объемов и схемой Кранка–Николсон на кубической сетке с $L = 3$ уровнями измельчения.

Наконец, мы демонстрируем способность метода работать с многогранными сетками общего вида. Мы рассматриваем шестигранные и треугольные призматические сетки с $L = 2$ и $L = 3$ уровнями измельчения, см. фиг. 15. Сетки содержат $16$ призматических слоев в $z$-направлении, а разрешение грубой сетки в $x$- и $y$-направлениях примерно соответствует разрешению $16 \times 16 \times 16$ кубической сетки.

Фиг. 15.

Изоповерхность $\psi = 0.5$ и срединный ($z = 0.5$) срез сетки при $t = \{ 0.5,\;1,\;1.5,\;2\} $, окрашенный согласно палитре $\psi $. Решение (2.2) нелинейным методом конечных объемов второго порядка со схемой Кранка–Николсон на гексагональной призматической сетке (H) и треугольной призматической сетке (P) с $L = 2$ уровнями измельчения ($\Delta t = 0.01$) и $L = 3$ уровнями измельчения ($\Delta t = 0.005$). В обоих случаях $K \approx 1$.

4.2. Тест Энрайта

Далее мы рассмотрим тест Энрайта [61]. Поле скоростей ${\mathbf{u}} = [u,{v},w{{]}^{{\text{T}}}}$ задано в единичном кубе $\Omega {{ = [0,1]}^{3}}$:

(4.2)
$\begin{gathered} u = 2\sin {{(\pi x)}^{2}}\sin (2\pi y)\sin (2\pi z)\cos (\pi t{\text{/}}3), \\ {v} = - \sin (2\pi x)\sin {{(\pi y)}^{2}}\sin (2\pi z)\cos (\pi t{\text{/}}3), \\ w = - \sin (2\pi x)\sin (2\pi y)\sin {{(\pi z)}^{2}}\cos (\pi t{\text{/}}3). \\ \end{gathered} $
Задача решается нелинейным методом конечных объемов со схемой Кранка–Николсон на временном интервале $t \in [0,1.5]$. В начальный момент времени $t = 0$ переносимый объект определяется сферой с радиусом $r = 0.15$ с центром в ${\mathbf{x}}{{ = [0.35,0.35,0.35]}^{{\text{T}}}}$. Решения на кубических сетках с уровнями измельчения $L = 3$ и $L = 4$ продемонстрированы на фиг. 16. Поверхность объекта остается гладкой, хотя на более мелкой сетке разрешение объекта намного лучше. Разрешение сетки 16-128 (с $L = 3$ уровнями измельчения) является недостаточным для поддержания связанной изоповерхности $\psi = 0.5$.

Фиг. 16.

Тест Энрайта: изоповерхность $\psi = 0.5$ при $T = \{ 0.3,\;0.75,\;1.2,\;1.5\} $. Кубическая сетка с уровнями измельчения $L = 3$ (б) и $L = 4$ (а).

5. ВЫВОДЫ

В работе предложен нелинейный метод конечных объемов для задачи переноса–сжатия границы раздела сред в рамках метода объема жидкости. Предложенная схема сохраняет решение в заданных границах $[0,1]$, применима к неструктурированным адаптивным сеткам общего вида и допускает большие временные шаги. Одним из возможных применений схемы является полностью неявное моделирование течений со свободной поверхностью.

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

Кирилл Терехов благодарит Брэдли Мэллисона и Хамди Челепи за обсуждение применения нелинейного метода конечных объемов к задаче переноса. Яшару Мехмани выражается благодарность за обсуждение проблемы сжатия межфазной границы.

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

  1. Ketabdari M.J. Free surface flow simulation using VOF method // Numerical Simulation: From Brain Imaging to Turbulent Flows. (Ed. Lpez-Ruiz). BoDBooks on Demand. 2016. V. 365.

  2. Nikitin K.D., Olshanskii M.A., Terekhov K.M. et al. An adaptive numerical method for free surface flows passing rigidly mounted obstacles // Comput. Fluids. 2017. V. 148. P. 56–68.

  3. Vassilevski Y.V., Nikitin K., Olshanskii M., Terekhov K. CFD technology for 3D simulation of large-scale hydrodynamic events and disasters // Russ. J. Numer. Anal. Math. Model. 2012. V. 27. № 4. P. 399–412.

  4. Kumar B., Crane M., Delauré Y. On the volume of fluid method for multiphase fluid flow simulation // Int. J. Model. Simul. Sci. Comput. 2013. V. 4. № 02. P. 1350002.

  5. Ruuth S.J., Wetton B.T. A simple scheme for volume-preserving motion by mean curvature // J. Sci. Comput. 2003. V. 19. № 1. P. 373–384.

  6. Qi Y., Lu J., Scardovelli R. et al. Computing curvature for volume of fluid methods using machine learning // J. Comput. Phys. 2019. V. 377. P. 155–161.

  7. Nikitin K.D., Olshanskii M.A., Terekhov K.M., Vassilevski Y.V. A splitting method for numerical simulation of free surface flows of incompressible fluids with surface tension // Comput. Methods Appl. Math. 2015. V. 15. № 1. P. 59–77.

  8. Nikitin K.D., Terekhov K.M., Vassilevski Y.V. Two methods of surface tension treatment in free surface flow simulations // Appl. Math. Lett. 2018. V. 86. P. 236–242.

  9. McFadden S., Browne D. A front-tracking model to predict solidification macrostructures and columnar to equiaxed transitions in alloy castings // Appl. Math. Model. 2009. V. 33. № 3. P. 1397–1416.

  10. Malladi R., Sethian J.A. Level set methods for curvature flow, image enchancement, and shape recovery in medical images // Math. Vis. Springer. 1997. P. 329–345.

  11. Popinet S. Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries // J. Comput. Phys. 2003. V. 190. № 2. P. 572–600.

  12. Lalanne C., Magdelaine Q., Lequien F., Fullana J.-M. Numerical model using a volume-of-fluid method for the study of evaporating sessile droplets in both unpinned and pinned modes // Eur. J. Mech. B Fluids. 2021.

  13. Kunkelmann C., Stephan P. CFD simulation of boiling flows using the volume-of-fluid method within OpenFOAM // Numer. Heat Transf. A. 2009. V. 56. № 8. P. 631–646.

  14. Gamet L., Scala M., Roenby J. et al. Validation of volume-of-fluid OpenFOAM® isoadvector solvers using single bubble benchmarks // Comput. Fluids. 2020. V. 213. P. 104722.

  15. Albadawi A., Donoghue D., Robinson A. et al. On the analysis of bubble growth and detachment at low capillary and bond numbers using volume of fluid and level set methods // Chem. Eng. Sci. 2013. V. 90. P. 77–91.

  16. Puckett E.G., Almgren A.S., Bell J.B. et al. A high-order projection method for tracking fluid interfaces in variable density incompressible flows // J. Comput. Phys. 1997. V. 130. № 2. P. 269–282.

  17. Sussman M., Puckett E.G. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows // J. Comput. Phys. 2000. V. 162. № 2. P. 301–337.

  18. Cifani P., Michalek W., Priems G. et al. A comparison between the surface compression method and an interface reconstruction method for the VOF approach // Comput. Fluids. 2016. V. 136. P. 421–435.

  19. Rusche H. Computational fluid dynamics of dispersed two-phase flows at high phase fractions: Ph.D. thesis / Imperial College London (University of London). 2003.

  20. Okagaki Y., Yonomoto T., Ishigaki M., Hirose Y. Numerical study on an interface compression method for the volume of fluid approach // Fluids. 2021. V. 6. № 2. P. 80.

  21. Aboukhedr M., Georgoulas A., Marengo M. et al. Simulation of micro-flow dynamics at low capillary numbers using adaptive interface compression // Comput. Fluids. 2018. V. 165. P. 13–32.

  22. Patel J.K., Natarajan G. A generic framework for design of interface capturing schemes for multi-fluid flows // Comput. Fluids. 2015. V. 106. P. 108–118.

  23. Arote A., Bade M., Banerjee J. An improved compressive volume of fluid scheme for capturing sharp interfaces using hybridization // Numer. Heat Transf. B: Fundam. 2020. V. 79. № 1. P. 29–53.

  24. Mehmani Y. Wrinkle-free interface compression for two-fluid flows // arXiv preprint arXiv:1811.09744. 2018.

  25. Piro D.J., Maki K. An adaptive interface compression method for water entry and exit: Tech. Rep. 2013-350: University of Michigan. Department of Naval Architecture and Marine Engineering, 2013.

  26. Lee H., Rhee S.H. A dynamic interface compression method for VOF simulations of high-speed planing watercraft // J. Mech. Sci. Technol. 2015. V. 29. № 5. P. 1849–1857.

  27. Sethian J.A., Smereka P. Level set methods for fluid interfaces // Annu. Rev. Fluid Mech. 2003. V. 35. № 1. P. 341–372.

  28. Adalsteinsson D., Sethian J.A. The fast construction of extension velocities in level set methods // J. Comput. Phys. 1999. V. 148. № 1. P. 2–22.

  29. Terekhov K.M., Nikitin K.D., Olshanskii M.A., Vassilevski Y.V. A semi-Lagrangian method on dynamically adapted octree meshes // Russ. J. Numer. Anal. Math. Model. 2015. V. 30. № 6. P. 363–380.

  30. Nikitin K., Olshanskii M., Terekhov K., Vassilevski Yu.V. Preserving distance property of level set function and simulation of free surface flows on adaptive grids // Numerical Geometry, Grid Generation and Scientific Computing (NUMGRID-2010). 2010. P. 25–32.

  31. Ausas R.F., Dari E.A., Buscaglia G.C. A geometric mass-preserving redistancing scheme for the level set function // Int. J. Numer. Methods Fluids. 2011. V. 65. № 8. P. 989–1010.

  32. Ge Z., Loiseau J.-C., Tammisola O., Brandt L. An eficient masspreserving interface-correction level set/ghost fluid method for droplet suspensions under depletion forces // J. Comput. Phys. 2018. V. 353. P. 435–459.

  33. Ningegowda B.M., Ge Z., Lupo G. et al. A mass-preserving interfacecorrection level set/ghost fluid method for modeling of three-dimensional boiling flows // Int. J. Heat Mass Transf. 2020. V. 162. P. 120382.

  34. Guermond J.-L., de Luna M.Q., Thompson T. An conservative antidi fusion technique for the level set method // J. Comput. Appl. Math. 2017. V. 321. P. 448–468.

  35. Leonard B., Mokhtari S. Beyond first-order upwinding: The ultra-sharp alternative for non-oscillatory steady-state simulation of convection // Int. J. Numer. Methods Eng. 1990. V. 30. № 4. P. 729–766.

  36. Silva L., Fontes C., Lage P. Front tracking in recirculating flows: a comparison between the TVD and RCM methods in solving the VOF equation // Braz. J. Chem. Eng. 2005. V. 22. № 1. P. 105–116.

  37. Lu C.-N., Wang R.-Y., Sun J.-S. WENO finite volume method for tracking moving interfaces on unstructured triangle meshes // J. Hohai Univ. 2009. V. 37. № 1. P. 105–109.

  38. Pirozzoli S., Di Giorgio S., Iafrati A. On algebraic TVD-VOF methods for tracking material interfaces // Comput. Fluids. 2019. V. 189. P. 73–81.

  39. Darwish M., Moukalled F. Convective schemes for capturing interfaces of free-surface flows on unstructured grids // Numer. Heat Transf. B: Fundam. 2006. V. 49. № 1. P. 19–42.

  40. Ubbink O., Issa R. A method for capturing sharp fluid interfaces on arbitrary meshes // J. Comput. Phys. 1999. V. 153. № 1. P. 26–50.

  41. Tsui Y.-Y., Lin S.-W., Cheng T.-T., Wu T.-C. Flux-blending schemes for interface capture in two-fluid flows // Int. J. Heat Mass Transf. 2009. V. 52. № 23–24. P. 5547–5556.

  42. Zhang D., Jiang C., Liang D. et al. A refined volume-of-fluid algorithm for capturing sharp fluid interfaces on arbitrary meshes // J. Comput. Phys. 2014. V. 274. P. 709–736.

  43. Bertolazzi E., Manzini G. A second-order maximum principle preserving finite volume method for steady convection-diffusion problems // SIAM J. Numer. Anal. 2005. V. 43. № 5. P. 2172–2199.

  44. Droniou J., Potier C.L. Construction and convergence study of schemes preserving the elliptic local maximum principle // SIAM J. Numer. Anal. 2011. V. 49. № 2. P. 459–490.

  45. Lipnikov K., Svyatskiy D., Vassilevski Y. Minimal stencil finite volume scheme with the discrete maximum principle // Russ. J. Numer. Anal. Math. Model. 2012. V. 27. № 4. P. 369–386.

  46. Chernyshenko A., Vassilevski Y. A finite volume scheme with the discrete maximum principle for diffusion equations on polyhedral meshes // Finite Volumes for Complex Applications VII-Methods and Theoretical Aspects. Springer. 2014. P. 197–205.

  47. Terekhov K.M., Mallison B.T., Tchelepi H.A. Cell-centered nonlinear finite-volume methods for the heterogeneous anisotropic diffusion problem // J. Comput. Phys. 2017. V. 330. P. 245–267.

  48. Lee J., Sheen D. A parallel method for backward parabolic problems based on the Laplace transformation // SIAM J. Numer. Anal. 2006. V. 44. № 4. P. 1466–1486.

  49. Varga R.S. Matrix Iterative Analysis. Prentice-Hall Series in Automatic Computation. Englewood Cliffs: Prentice-Hall. 1962.

  50. Jameson A. Analysis and design of numerical schemes for gas dynamics, 1: artificial diffusion, upwind biasing, limiters and their effect on accuracy and multigrid convergence // Int. J. Comut. Fluid Dyn. 1995. V. 4. № 3–4. P. 171–218.

  51. Lipnikov K., Svyatskiy D., Vassilevski Y.V. Anderson acceleration for nonlinear finite volume scheme for advection-diffusion problems // SIAM J. Sci. Comput. 2013. V. 35. № 2.

  52. Perot B. Conservation properties of unstructured staggered mesh schemes // J. Comput. Phys. 2000. V. 159. № 1. P. 58–89.

  53. Younis R., Tchelepi H.A., Aziz K. Adaptively localized continuation-Newton method–nonlinear solvers that converge all the time // Soc. Pet. Eng. J. 2010. V. 15. № 02. P. 526–544.

  54. Vassilevski Y., Terekhov K., Nikitin K., Kapyrin I. Parallel Finite Volume Computation on General Meshes. Springer Nature. 2020.

  55. Terekhov K. Parallel multilevel linear solver within INMOST platform // Russian Supercomputing Days / Springer. 2020. P. 297–309.

  56. Terekhov K. Greedy dissection method for shared parallelism in incomplete factorization within INMOST platform // Russian Supercomputing Days / Springer. 2021. P. 87–101.

  57. Terekhov K., Vassilevski Y. Mesh modification and adaptation within INMOST programming platform // Numerical Geometry, Grid Generation and Scientific Computing. Springer. 2019. P. 243–255.

  58. Terekhov K. Parallel dynamic mesh adaptation within INMOST platform // Russian Supercomputing Days / Springer. 2019. P. 313–326.

  59. Karypis G., Schloegel K., Kumar V. Parmetis parallel graph partitioning and sparse matrix ordering library: Tech. Rep. 97-060: University of Minnesota. Department of Computer Science and Engineering. 1997.

  60. Karypis G., Kumar V. MeTis: Unstructured Graph Partitioning and Sparse Matrix Ordering System, Version 4.0. http://www.cs.umn.edu/ metis. 2009.

  61. Enright D., Fedkiw R., Ferziger J., Mitchell I. A hybrid particle level set method for improved interface capturing // J. Comput. Phys. 2002. V. 183. № 1. P. 83–116.

  62. Ahrens J., Geveci B., Law C. Paraview: An end-user tool for large data visualization // The Visualization Handbook. 2005. V. 717. № 8.

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

Инструменты

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

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

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