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

Журнал вычислительной математики и математической физики, 2019, T. 59, № 12, стр. 2007-2023

Зачем нужны сетки Вороного–Делоне? Основные свойства метода конечных объемов с использованием ячеек Вороного

К. Гертнер 1, *, Л. Каменски 1

1 m4sim GmbH, Seydelstr. 31
10117 Berlin, Germany

* E-mail: info@m4sim.de

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

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

Аннотация

В отличие от схем, которые локально нарушают существенные свойства устойчивости параболических и эллиптических задач, методы конечного объема в мультиматериальных областях с использованием ячеек Вороного и сеток Делоне, согласованных с границей, в которых центры приграничных шаров Делоне принадлежат расчетной области, обеспечивают хорошее приближение геометрии задачи и способны при этом сохранять существенные качественные свойства решения при любом разрешении в пространстве и времени, а также при изменениях временных масштабов на несколько порядков величин. К сожалению, в последнее время эти методы используются все реже и реже, поскольку отсутствуют автоматические генераторы сеток с нужными свойствами. Данная работа дает краткое описание ключевых свойств методов конечного объема с использованием сеток Вороного–Делоне с примерами из практики применения и мотивацию для более широкого использования таких методов и для разработки алгоритмов построения сеток Делоне–Вороного. Библ. 37. Фиг. 9.

Ключевые слова: метод конечных объемов, сетки Вороного–Делоне, сетки Делоне, согласованные с границей.

1. ВВЕДЕНИЕ

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

В дискретном случае это анализ менее очевиден. Зачастую параметры дискретизации во времени и пространстве, $h$ и $\tau $, которые на самом деле являются пределами, рассматриваются как конечные и универсально достижимые. Однако это не всегда бывает в ситуациях с ограниченной гладкостью, например, в состояниях удаленных от термодинамического равновесия или при более высокой плотности свободной энергии, и особенно если требуется приближение негладких границ. Решения могут иметь плато и внутренние или пограничные слои, соединяющие плато и граничные значения. Плотность энергии часто очень велика по сравнению с привычной средой обитания человека. Известным примером является перенос заряда в полупроводниках, где большая разница в силе является результатом того, что разность электростатического потенциала всего лишь в 1 В соответствует разности высот в 4000 км в гравитационном поле Земли (механическое отношение эквивалента тепла на 300); соответственно 100 В, что не так уж много для таких задач, соответствует разнице высот от Земли до места где-то за Луной. Многие задачи из области электрохимии и химических реакций и из других областей науки приводят к такому же типу проблем и требуют сохранения качественных свойств в их дискретной форме.

В таких случаях классическая ошибка дискретизации имеет лишь ограниченную значимость с точки зрения гарантии устойчивости дискретного решения. Часто даже физически не имеет значения, достигается ли стационарное предельное состояние проблемы при наименьшем масштабе времени в 1 с за ${{10}^{{19}}}$ или ${{10}^{{20}}}$ с, поскольку время реакции и продолжительность жизни человека составляют соответственно 1 с и 100 лет ($100\;{\text{лет}} \approx 100 \times 365 \times 24 \times 3600\;{\text{с}} \approx 3e + 9\;{\text{с}}$). Напротив, большой интерес представляют стационарные состояния и качественные характеристики, например, обеспечивает ли модель А лучшее решение в гигагерцовом диапазоне, чем модель B. Поэтому существенным является сохранение структурных свойств задач в пространстве и времени для конечных и не слишком малых параметров дискретизации.

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

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

Идеи устойчивой дискретизации для упрощенных задач восходят к 1950-м годам [3], [19], [29] и в обобщенном виде остаются актуальными и по сей день. Диаграмма Вороного [36] и соответствующая ей триангуляция Делоне [6] использовались уже в [27] (и, например, в классической книге Варги [35]). Двумерные триангуляции Делоне являются оптимальными для различных критериев, например, для максимизации наименьшего угла триангуляции [25], для минимизации квадрата ${{L}^{2}}$-нормы градиента кусочно-линейной интерполяции, и для других критериев (см., например, [31, Sect. 4] и список литературы к статье). Трудно проследить историю развития сеток Делоне в деталях, но основные моменты были связаны с ${{C}^{1}}$ поверхностной интерполяцией и с такими задачами, как перенос нейтронов и перенос зарядов в полупроводниках. Допустимые сетки [7] подводят итог для минимально необходимых свойств сеток, позволяющих доказать сходимость аппроксимаций законов сохранения.

По-видимому, местом рождения термина сетки Делоне, согласованные с границами – boundary conforming Delaunay meshes, была Вена (см., например, [8]). Такие сетки требуются в случае межматериальных границ (интерфейсов) со скачками свойств материала (например, коэффициентов диффузии), нелинейных граничных условий III рода, или в некоторых других особых случаях, чтобы гарантировать, что границы и внутренние интерфейсы области задачи являются частью сетки Вороного–Делоне, в то время как сетки Делоне с ограничениями – constrained Delaunay meshes [5], ставшие популярными в последнее время [30], [32], для этой цели не подходят. Локальные ортогональные координаты были бы лучшим выбором. Несмотря на то что они часто не существуют глобально, т.е. во всей области, локально они могут стать ключом к сочетанию устойчивости, хороших аппроксимирующих свойств и невысокой вычислительной сложности дискретной задачи.

Ниже перечислены некоторые желаемые свойства дискретной задачи по сравнению с дифференциальной постановкой:

• те же качественные свойства устойчивости, что и у точной задачи;

• никаких искусственных предположений о малости или гладкости каких-то величин;

• использование слабой формы дискретной задачи, а также хорошо известных тестовых функций и алгоритмов;

• следовательно, слабый дискретный энергетический функционал;

• слабая сходимость решений, когда правая часть задается как распределение;

• расположение внутренних слоев как конечных разрывов в решении (случай неразрешенных слоев).

Примерами, где список частично выполнен, являются: перенос электронов и дырок в полупроводниках и разделение фаз (не Cahn-Hilliard) [10]–[12], [14], [15], [17].

2. СЕТКИ ДЕЛОНЕ, ЯЧЕЙКИ ВОРОНОГО И СЕТКИ ДЕЛОНЕ, СОГЛАСОВАННЫЕ С ГРАНИЦЕЙ

Рассмотрим дискретизацию (сетку) ${{T}_{h}}(\Omega )$ области $\Omega \subset {{\mathbb{R}}^{d}}$ в предположении, что $\Omega $ – это объединение конечных $d$-мерных подобластей, $\Omega = \bigcup\nolimits_i^{} {{{\Omega }_{i}}} $. Каждая подобласть содержит только один материал и состоит из $d$-мерных симплексов. ${\mathbf{E}}_{j}^{d}$, ${{\Omega }_{i}} = \bigcup\nolimits_j^{} {{\mathbf{E}}_{j}^{d}} $. Вершины симплексов сетки обозначим через ${{{\mathbf{v}}}_{k}}$. Следовательно, все интерфейсы между материалами и границы области составлены из симплексов меньшей размерности ${\mathbf{E}}_{k}^{m}$, $1 \leqslant m < d$. Предполагается, что описание геометрии области содержит все критические линии и точки.

2.1. Триангуляция Делоне

Дискретизация ${{T}_{h}}(\Omega )$ заданной $d$-мерной области $\Omega $ симплексами называется сеткой Делоне или триангуляцией Делоне, если она удовлетворяет свойству пустого шара, то есть если никакая вершина сетки не находится внутри описанной сферы любого симплекса, которая является границей шара, называемого шаром Делоне:

Для заданной сетки, симплекс сетки называется симплексом Делоне, если его шар Делоне пуст (фиг. 1). Ребро сетки называется ребром Делоне, если существует пустая сфера, содержащая это ребро, но не имеющая других вершин сетки внутри себя. Следовательно, симплекс является Делоне тогда и только тогда, когда все его ребра являются ребрами Делоне. Ребра Делоне не обязательно являются самыми короткими из возможных ребер (см. фиг. 1, справа), но двумерная триангуляция Делоне максимизирует минимальный угол триангуляции [25], а также квадрат ${{L}^{2}}$-нормы градиента кусочно-линейной интерполяции [28].

Фиг. 1.

$\Delta {{{\mathbf{v}}}_{1}}{{{\mathbf{v}}}_{2}}{{{\mathbf{v}}}_{3}}$удовлетворяет условию Делоне, потому что его окружность пуста (слева), в то время как $\Delta {{{\mathbf{v}}}_{1}}{{{\mathbf{v}}}_{2}}{{{\mathbf{v}}}_{4}}$ не удовлетворяет условию Делоне, потому что его описанная окружность содержит в себе ${{{\mathbf{v}}}_{3}}$ (в центре). Ребро Делоне не обязательно является самым коротким из возможных: ${{{\mathbf{v}}}_{2}}{{{\mathbf{v}}}_{4}}$ удовлетворяет условию Делоне, но не является самым коротким из возможных (справа).

Триангуляция Делоне однозначна, если вершины сетки находятся в общем положении, т.е. никакие $d + 2$ вершины не лежат на общей пустой $d$-сфере.

2.2. Объемы, поверхности и ячейки Вороного

Объем Вороного ${{V}_{i}}$ вершины ${{{\mathbf{v}}}_{i}}$ – это часть ${{R}^{d}}$, состоящая из всех точек, для которых вершина ${{{\mathbf{v}}}_{i}}$ ближайшая из всех вершин сетки:

${{V}_{i}}: = \left\{ {{\mathbf{x}} \in {{\mathbb{R}}^{d}}:\left\| {{\mathbf{x}} - {{{\mathbf{v}}}_{i}}} \right\| < \left\| {{\mathbf{x}} - {{{\mathbf{v}}}_{j}}} \right\|\forall {{{\mathbf{v}}}_{j}} \in {{\mathcal{T}}_{h}}(\Omega )} \right\}.$
Граница $\partial {{V}_{i}}$ объема Вороного называется соответствующей поверхностью Вороного (фиг. 2, слева).

Фиг. 2.

Объем ${{V}_{4}}$, поверхность $\partial {{V}_{4}}$ (слева), и элемент объема (ячейка) Вороного ${{V}_{{4j}}} = {{V}_{4}} \cap {\mathbf{E}}_{j}^{2}$ (справа).

Разбиение Вороного для заданного набора вершин является двойственной структурой по отношению к триангуляции Делоне. Каждая грань Вороного соответствует ребру триангуляции Делоне.

Пересечение объема Вороного ${{V}_{i}}$ с симплексом ${\mathbf{E}}_{j}^{d}$, ${{V}_{{ij}}} = {{V}_{i}} \cap {\mathbf{E}}_{j}^{d}$, является ячейкой или элементом объема Вороного вершины ${{{\mathbf{v}}}_{i}}$ по отношению к ${\mathbf{E}}_{j}^{d}$ (фиг. 2, справа). Будучи пересечением полупространств, ячейки и грани Вороного всегда выпуклые, и их меры являются непрерывными функциями координат вершин.

2.3. Триангуляции (сетки) Делоне, согласованные с границей

Введение граничных или межфазных условий в сетку может нарушить свойство Делоне. Например, рассмотрим однородные граничные условия Неймана для границы области (внутренней или внешней) фиг. 2 вдоль ребер ${{{\mathbf{v}}}_{1}}{{{\mathbf{v}}}_{3}}$ и ${{{\mathbf{v}}}_{3}}{{{\mathbf{v}}}_{2}}$ (фиг. 3, слева). В схемах, в которых величины хранятся в ячейках, распространенной техникой для граничных условий Неймана является введение дополнительной виртуальной вершины ${\mathbf{v}}_{4}^{*}$ как зеркального отображения ${{{\mathbf{v}}}_{4}}$ через ${{{\mathbf{v}}}_{1}}{{{\mathbf{v}}}_{3}}$. Если виртуальная вершина ${\mathbf{v}}_{4}^{*}$ находится внутри описанного круга для треугольника $\Delta {{{\mathbf{v}}}_{1}}{{{\mathbf{v}}}_{3}}{{{\mathbf{v}}}_{4}}$ (что произойдет тогда и только тогда, если ${{{\mathbf{v}}}_{4}}$ находится внутри описанного круга для ${{{\mathbf{v}}}_{1}}{{{\mathbf{v}}}_{3}}$), то треугольник $\Delta {{{\mathbf{v}}}_{1}}{{{\mathbf{v}}}_{4}}{{{\mathbf{v}}}_{3}}$ нарушит свойство Делоне. Эта проблема присутствует и в более высоких измерениях и мотивирует следующее определение.

Фиг. 3.

Если угол $\angle {{{\mathbf{v}}}_{1}}{{{\mathbf{v}}}_{4}}{{{\mathbf{v}}}_{3}}$ тупой, то виртуальная вершина ${\mathbf{v}}_{4}^{*}$ для граничного условия Неймана нарушит условие Делоне (слева). Вставка новой вершины как ортогональной проекции ${{{\mathbf{v}}}_{4}}$ на ${{{\mathbf{v}}}_{1}}{{{\mathbf{v}}}_{3}}$ восстанавливает граничное соответствие сетки и не нарушает свойства Делоне ребер ${{{\mathbf{v}}}_{1}}{{{\mathbf{v}}}_{4}}$ и ${{{\mathbf{v}}}_{4}}{{{\mathbf{v}}}_{3}}$ (в центре); жирная пунктирная линия соответствует граничному условию, имеющему отношение к ${{{\mathbf{\tilde {v}}}}_{4}}$.

Сетка Делоне называется сеткой Делоне, согласованной с границей, если ни одна вершина сетки не находится внутри самой малой (диаметральной, экваториальной) сферы какого угодно маломерного симплекса границы или внутреннего интерфейса:

Свойство пустого диаметрального шара иногда называют свойством Габриэля [9]. Иными словами, сетка Делоне, согласованная с границей, представляет собой сетку, где каждый симплекс удовлетворяет критерию Делоне, а каждый граничный симплекс – критерию Габриэля.

В ситуации фиг. 3 (слева) согласование с границей можно достичь, разделив $\Delta {{{\mathbf{v}}}_{1}}{{{\mathbf{v}}}_{3}}{{{\mathbf{v}}}_{4}}$ новой вершиной на граничном ребре ${{{\mathbf{v}}}_{1}}{{{\mathbf{v}}}_{3}}$. Пересечение диаметральных дисков ${{v}_{1}}{{v}_{4}}$ и ${{v}_{4}}{{v}_{3}}$ с ребром границы ${{v}_{1}}{{v}_{3}}$ соответствует проекции ${{v}_{4}}$ на ${{v}_{1}}{{v}_{3}}$ (фиг. 3, в центре) и гарантирует восстановление согласования сетки с границей (фиг. 3, справа) так, что при этом все остальные симплексы сетки сохранят свойство Делоне.

Подробности о сетках Делоне, согласованных с границей, можно найти в работе [33] и в ссылках в ней.

3. ОСНОВНЫЕ СВОЙСТВА МЕТОДА КОНЕЧНЫХ ОБЪЕМОВ С ИСПОЛЬЗОВАНИЕМ ЯЧЕЕК ВОРОНОГО (VFVM)

Если применить теорему Гаусса к ячейке Вороного или ее пересечению с $\partial \Omega $, то каждая грань двойственна ребру и, следовательно, ее можно представить как набор симплексов. Более удобно ассемблировать уравнения конечных объемов с помощью симплексов двойственной сетки Делоне, в частности, гораздо удобнее определять свойства материалов на симплексах, чем на объемах Вороного.

3.1. Уравнения конечных объемов с использованием ячеек Вороного

Матрица смежности ${{\tilde {G}}_{i}}$ для $i$-мерного симплекса, отображающая вершины на ребра может быть определена рекурсивно через

${{\tilde {G}}_{{i + 1}}} = \left( {\begin{array}{*{20}{c}} {{{{\tilde {G}}}_{i}}}&{{{{\mathbf{0}}}_{{\tfrac{{i(i - 1)}}{2}}}}} \\ { - {{\mathbb{I}}_{{i + 1}}}}&{{{{\mathbf{1}}}_{{i + 1}}}} \end{array}} \right),\quad i = 1,2, \ldots ,$
где ${{\mathbb{I}}_{j}}$ обозначает единичную матрицу $j \times j$, ${{{\mathbf{0}}}_{j}} = {{(\overbrace {0, \ldots ,0}^j)}^{{\text{T}}}}$, ${{{\mathbf{1}}}_{j}} = {{(\overbrace {1, \ldots ,1}^j)}^{{\text{T}}}}$, а ${{\tilde {G}}_{1}} = ( - 1,1)$ определяет разность на интервале,
${{\tilde {G}}_{1}}\left( {\begin{array}{*{20}{c}} {{{{\mathbf{x}}}_{1}}} \\ {{{{\mathbf{x}}}_{2}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} { - 1}&1 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{{{\mathbf{v}}}_{1}}} \\ {{{{\mathbf{v}}}_{2}}} \end{array}} \right) = {{{\mathbf{v}}}_{2}} - {{{\mathbf{v}}}_{1}}.$
Получаем

${{\tilde {G}}_{2}} = \left( {\begin{array}{*{20}{c}} {{{{\tilde {G}}}_{1}}}&{{{{\mathbf{0}}}_{1}}} \\ { - {{\mathbb{I}}_{2}}}&{{{{\mathbf{1}}}_{2}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} { - 1}&1&\vline & 0 \\ \hline { - 1}&0&\vline & 1 \\ 0&{ - 1}&\vline & 1 \end{array}} \right),$
${{\tilde {G}}_{3}} = \left( {\begin{array}{*{20}{c}} {{{{\tilde {G}}}_{2}}}&{{{{\mathbf{0}}}_{3}}} \\ { - {{\mathbb{I}}_{3}}}&{{{{\mathbf{1}}}_{3}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} { - 1}&1&0&\vline & 0 \\ { - 1}&0&1&\vline & 0 \\ 0&{ - 1}&1&\vline & 0 \\ \hline { - 1}&0&0&\vline & 1 \\ 0&{ - 1}&0&\vline & 1 \\ 0&0&{ - 1}&\vline & 1 \end{array}} \right),$
${{\tilde {G}}_{4}} = \left( {\begin{array}{*{20}{c}} {{{{\tilde {G}}}_{3}}}&{{{{\mathbf{0}}}_{6}}} \\ { - {{\mathbb{I}}_{4}}}&{{{{\mathbf{1}}}_{4}}} \end{array}} \right),\;...$

Матрица ${{\tilde {G}}_{{3d}}}$ отображает вершины на ребра. Соответственно $\tilde {G}_{d}^{{\text{T}}}$ отображает ребра на вершины, и $\tilde {G}_{d}^{{\text{T}}}{{\tilde {G}}_{d}}$ является матрицей Кирхгофа (Лапласианом графа) $d$-мерного симплекса.

Обозначим длину ребра ${\mathbf{E}}_{i}^{1}$ как $\left| {{{{\tilde {G}}}_{1}}{\mathbf{E}}_{i}^{1}} \right|$ и диагональную матрицу длин всех ребер симплекса ${\mathbf{E}}_{i}^{d}$ как $\left[ {{\text{|}}{{{\tilde {G}}}_{d}}{\mathbf{E}}_{i}^{d}{\text{|}}} \right]$. Тогда

$\mathop {\left[ {{\text{|}}{{{\tilde {G}}}_{d}}{\mathbf{E}}_{i}^{d}{\text{|}}} \right]}\nolimits^{ - 1} {{\tilde {G}}_{d}}{\mathbf{u}}$
является приближением значения $\nabla u$ в направлении $\partial u{\text{/}}\partial {{n}_{k}}$, перпендикулярном грани ячейки Вороного, относительно симплекса ${\mathbf{E}}_{i}^{d}$, которое преобразовывает заданные линейные функции на каждом ребре в корректную постоянную. Далее,
$\tilde {G}_{d}^{{\text{T}}}[{{\sigma }_{{{\mathbf{E}}_{i}^{d}}}}] = \tilde {G}_{d}^{{\text{T}}}{{\left[ {{\text{|}}{{{\tilde {G}}}_{d}}{\mathbf{E}}_{i}^{d}{\text{|}}} \right]}^{{ - 1}}}\mathop {\left\{ {\left[ {{\text{|}}{{{\tilde {G}}}_{d}}{\mathbf{E}}_{i}^{d}{\text{|}}} \right]\left[ {{{\sigma }_{{{\mathbf{E}}_{i}^{d}}}}} \right]} \right\}}\nolimits^{\tfrac{1}{2}} $
является симметричной формой отображения с ребер на вершины, включающей геометрический вес грани Вороного, пересекаемой ребром, умноженным на длину ребра. Следовательно, метрический фактор $\{ \cdot \} $ имеет смысл удвоенного детерминанта (длина ребра в два раза больше высоты $d$-мерного симплекса, образованного пересечением симплекса и ${{\sigma }_{i}}$), а не объема, пересеченного с соответствующей частью каждого ребра симплекса ${\mathbf{E}}_{i}^{d}$. Это правило надо последовательно соблюдать во всех частях оператора и граничных условий. Следовательно,
$\varepsilon \mathop \smallint \limits_{{\mathbf{E}}_{i}^{d}} \nabla \cdot \nabla u{\text{d}}V = \oint\limits_{\Sigma } {\nabla u{\text{d}}{\mathbf{S}}} \approx \varepsilon \tilde {G}_{d}^{{\text{T}}}[{{\sigma }_{i}}]{{\left[ {\left| {{{{\tilde {G}}}_{d}}{\mathbf{E}}_{i}^{d}} \right|} \right]}^{{ - 1}}}{{\tilde {G}}_{d}}{\mathbf{u}}: = \varepsilon G_{d}^{{\text{T}}}{{G}_{d}}{\mathbf{u}}$
является приближением оператора Лапласа, проинтегрированного по объему на ${\mathbf{E}}_{i}^{d}$, где $\varepsilon = {\text{const}}$ на ${\mathbf{E}}_{i}^{d}$ и ${{\sigma }_{j}} \geqslant 0$ для каждого ребра $j$.

Граничное условие

$\alpha (u)u + \beta (u)\frac{{\partial u}}{{\partial n}} + \gamma (u) = 0$
интегрируется на гранях Вороного, принадлежащих к поверхностным вершинам; грань определяет нормальное направление, которое совпадает с направлением ребра, трансверсального к границе, в случае сетки Делоне, согласованной с границей:
$\varepsilon \frac{{\partial u}}{{\partial n}}\partial {{V}_{{s,i}}} = - \partial {{V}_{{s,i}}}\frac{\varepsilon }{{\beta (u)}}\left( {{{u}_{i}}\alpha ({{u}_{i}}) + \gamma ({{u}_{i}})} \right).$
Нормальная производная вызывает изменение соответствующего диагонального элемента и правой части дискретных уравнений. Условие Дирихле получается для $\beta \approx \varepsilon _{{{\text{machine}}}}^{2}$ и легко гомогенизируется.

Сетка Делоне, согласованная с границей, позволяет корректно разрезать ячейки Вороного вдоль границы, задаваемой симплексами, и ввести то же поверхностное интегрирование на $d - 1$-мерных ячейках Вороного. В случае пограничных слоев вершины на границе могут быть носителями экстремумов, введенных граничным условием. Это будет не так, если для приближения подобласти $\Omega _{i}^{d}$ использовать ячейки Вороного.

Основным условием для метода конечных объемов является справедливость теоремы Гаусса. Наибольшее функциональное пространство, где действует теорема Гаусса – это пространство функций ограниченной вариации. Это очень слабое условие на класс функций.

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

3.2. Механизм компенсации и свойство S-матрицы

Основное преимущество метода конечных объемов с использованием ячеек Вороного (VFVM) состоит в том, что матрица системы является невырожденной M-матрицей. M-матрица (матрица Минковского) – это Z-матрица (матрица с неположительными внедиагональными элементами) с положительной диагональю, сумма всех столбцов которой неотрицательна. Если хотя бы одна столбцовая сумма М-матрицы положительна, то эта матрица обратима и обратная матрица является положительной. Симметричная М-матрица называется S-матрицей (матрица Стильтеса).

Рассмотрим случай ${{\mathbb{R}}^{3}}$ и два допустимых тетраэдра Делоне (фиг. 4). Для удобства изложения скажем, что свободная четвертая вершина, противолежащая каждой треугольной грани, определяет четыре условия Делоне. На фиг. 4 (слева) все центры описанных окружностей (далее $C$ с индексом соответствующего, возможно низкоразмерного симплекса) находятся внутри соответствующего определяющего симплекса (соответственно одно-, двух- или трехмерного).

Фиг. 4.

Слева: все центры окружностей находятся внутри их определяющих $d$-симплексов, $d = 1,2,3$. Справа: центр $C$ окружности тетраэдра и центр окружности треугольника $\Delta 234$ выходят за пределы их определяющих симплексов.

На фиг. 4 (справа) центр ${{C}_{{1234}}}$ окружности тетраэдра и центр ${{C}_{{234}}}$ окружности треугольника $\Delta 234$ находятся вне их определяющих симплексов. Каждая из граней Вороного определяется двумя треугольниками, например: грань, определяемая вершинами 1 и 2, состоит из $\Delta {{C}_{{12}}}{{C}_{{123}}}{{C}_{{1234}}}$ и $\Delta {{C}_{{21}}}{{C}_{{124}}}{{C}_{{1234}}}$, которые перпендикулярны ребру 12 и изменяют знак в соответствии с положением соответствующих центров окружности треугольника и тетраэдра и имеют положительную проекцию на ребра 12 и 21 для внутренних центров окружностей. Соседний тетраэдр, который тоже имеет $\Delta 124$ в качестве грани, имеет те же центры описанных окружностей общих ребер и то же направление нормали к $\Delta 124$. Грань Вороного ограничена полигоном, соединяющим все центры окружностей тетраэдров, разделяющих одно ребро и расположенных в плоскости, перпендикулярной этому ребру. Суммирование всех вкладов тетраэдров, разделяющих общее ребро, дает площадь грани Вороного (со знаком), относящуюся к этому ребру. Условие Делоне гарантирует, что расстояние между центрами окружностей двух соседних тетраэдров является неотрицательным. Таким образом, значение площади грани Вороного неотрицательно.

Этот эффект называется компенсацией: отрицательный вклад компенсируется положительным вкладом соседнего тетраэдра. Он присутствует в любом измерении и гарантирует, что матрица в методе конечного объема является S-матрицей или M-матрицей.

Если на общей $d$-сфере лежат более $m \geqslant d + 2$ вершин, то центры $m - d$ окружностей будут совпадать. Сетка Делоне с ограничениями гарантирует, что все центры окружностей находятся внутри соответствующих определяющих симплексов, части ячеек и поверхностей Делоне являются неотрицательными, и сохраняются ориентированные объемы подобласти ${{\Omega }_{i}}$ и площади поверхностей. Это позволяет легко справляться со всеми граничными условиями III рода.

Для P1 (линейного) метода конечных элементов (МКЭ) в $d$-мерном пространстве матрица Дирихле является проекцией каждой из $d - 1$-мерной грани на все остальные,

${{P}^{{ - {\text{T}}}}}{{P}^{{ - 1}}},\quad P = {{{\mathbf{v}}}_{j}} - {{{\mathbf{v}}}_{0}},\quad j = 1, \ldots ,d.$
В двух измерениях она идентична матрице $G_{d}^{{\text{T}}}{{G}_{d}}$ для VFVM.

В трех- и более измерениях это уже не так (за исключением особых случаев) и компенсация негативных частей работает по-другому, так что граничные сетки Делоне не являются ни достаточными, ни необходимыми для получения S-матрицы для МКЭ (подробнее см. [4], [21], [24], [37]) и ссылки в них).

Так называемое условие остроугольности (все углы между граничными нормалями являются острыми), гарантирующее свойство S-матрицы для матрицы Дирихле для МКЭ является слишком ограничительным и сложным для использования на практике.

3.3. Принцип максимума и диссипативность

Решение VFVM удовлетворяет слабому дискретному принципу максимума. Пусть ${\mathbf{u}}$ удовлетворяет $G_{d}^{d}{{G}_{d}}{\mathbf{u}} = {\mathbf{0}}$ с заданным граничным условием Дирихле ${{{\mathbf{u}}}_{D}}$, $\hat {u} = max{{u}_{D}}$, $\mathop u\limits^ = min{{u}_{D}}$. Предположим, что $maxu = U > \hat {u}$. Проверка уравнения с помощью ${{({\mathbf{u}} - \hat {u})}^{{ + {\text{T}}}}}$ дает ${{({\mathbf{u}} - \hat {u})}^{{ + {\text{T}}}}}G_{d}^{{\text{T}}}{{G}_{d}}{\mathbf{u}} > 0$, потому что ${{G}_{d}}{\mathbf{u}}$ и ${{G}_{d}}\mathop {({\mathbf{u}} - \hat {u})}\nolimits^ + $ имеют либо одинаковый знак (и, как минимум, одно значение отлично от нуля) или ${{G}_{n}}\mathop {({\mathbf{u}} - \hat {u})}\nolimits^ + = {\mathbf{0}}$, что приводит к противоречию.

Диссипативность вытекает из этого же аргумента. В дискретном случае принцип максимума и диссипативность в этом смысле эквивалентны.

Замечание и диссипативности членов нулевого порядка: объемы Вороного образуют диагональную матрицу масс $[V]$, больший носитель невозможно получить без допущения малости вариации решения на носителе.

3.4. Невыполнения условия Делоне для оператора Лапласа

Так как ячейки Вороного имеют неотрицательные площади поверхности, ассемблирование

$A = G_{d}^{{\text{T}}}[\varepsilon ]{{G}_{d}}$
приводит к S-матрице: ${{A}_{{ii}}} > 0$, ${{A}_{{ij}}} \leqslant 0$, $i \ne j$, и ${\mathbf{1}}_{d}^{{\text{T}}}G_{d}^{{\text{T}}} = {{0}_{d}}$ (сумма в столбцах и строках равна нулю для вершин без граничных условий). Очевидно, что матрица ${{A}^{{ - 1}}}$ поэлементно положительна для неприводимой $A$ или связных сеток. Конечное значение $\varepsilon $ и граничные условия добавляют только положительные значения к диагонали по крайней мере для одной вершины. Следовательно, решение ограничено значениями Дирихле и положительно для нулевых значений Дирихле и положительной правой стороны уравнения.

Если заданная сетка не является сеткой Делоне, суммирование ориентированной площади поверхностей для каждого ребра может привести к образованию отрицательных поверхностей вокруг этого ребра. Это может быть произвольное ребро между двумя вершинами, не соответствующими грани Вороного (описанная сфера одного симплекса содержит центр другой). Далее следует ${{A}_{{ij}}} > 0$ и позитивность больше не гарантируется, по крайней мере, может случиться, что появятся локальные экстремумы в вершинах $i$ и $j$ для значительных отклонений от условия Делоне. По построению, ${{G}_{d}}{{{\mathbf{1}}}_{d}} = {{{\mathbf{0}}}_{d}}$. Введение ${{A}_{{ij}}} = 0$ для соответствия с ${{A}_{{ii}}}$ эквивалентно введению локального пограничного условия Неймана и несогласованной дискретизации.

Предположим, что решение $u$ имеет неразрешенный движущийся внутренний слой c параметрами

${{u}_{1}} = 1 - \varepsilon ,\quad {{u}_{j}} = {{\varepsilon }_{j}},\quad {\text{с}}\;0 < {{\varepsilon }_{j}} \leqslant \varepsilon \ll 1,\quad j = i,2, \ldots ,k,$
задающими порядок изменения $u$ за пределами скачка (фиг. 5). Здесь индексы – это просто номера точек сетки, в которых задана скалярная функция $u$. Пусть ${{A}_{{{\text{non}} - {\text{Del}}}}}$ отвечает типу ${{\tilde {G}}^{{\text{T}}}}[W]\tilde {G}$, a сумма берется по всем симплексам, относящимся к вершине $i$ (знак суммирования не указан, порядок сложения по ребрам или симплексам заданного типа может быть выбран соответственно каждой цели, и контекст определяет локальное или глобальное значение).

Фиг. 5.

Скачок решения вдоль ребра ${{{\mathbf{v}}}_{i}}{{{\mathbf{v}}}_{1}}$ от $u({{{\mathbf{v}}}_{i}}) = {{\varepsilon }_{i}}$ до $u({{{\mathbf{v}}}_{1}}) = 1 - \varepsilon $.

Пусть ${{l}_{{ij}}}$ будет длиной ребра для всех соседей ${{{\mathbf{v}}}_{j}}$ вершины ${{v}_{i}}$ и

$\begin{gathered} {{W}_{{ij}}} = \tfrac{{{{\sigma }_{{ij}}}}}{{{{l}_{{ij}}}}} = 1,\quad j = 2, \ldots ,k, \hfill \\ {{W}_{{i1}}} = \tfrac{{{{\sigma }_{{i1}}}}}{{{{l}_{{i1}}}}} = - 1, \hfill \\ \end{gathered} $
т.е. ребро ${{{\mathbf{v}}}_{i}}{{{\mathbf{v}}}_{1}}$ не удовлетворяет условию Делоне. Таким образом, ${{A}_{{{\text{non - del}}}}}{\mathbf{u}}{{{\text{|}}}_{i}}$ при вершине $i$ приводит к тому, что
${{A}_{{{\text{non}} - {\text{Del}}}}}{\mathbf{u}}{{{\text{|}}}_{i}} = \sum\limits_{j = 2}^k {({{u}_{i}} - {{u}_{j}}) - ({{u}_{i}} - {{u}_{1}})} = \sum\limits_{j = 2}^k {({{\varepsilon }_{i}} - {{\varepsilon }_{j}}) - ({{\varepsilon }_{i}} - (1 - \varepsilon ))} \approx 1 > 0,$
в то время как в случае Делоне
${{A}_{{{\text{Del}}}}}{\mathbf{u}}{{{\text{|}}}_{i}} = \sum\limits_{j = 2}^k {({{u}_{i}} - {{u}_{j}}) + ({{u}_{i}} - {{u}_{1}})} = \sum\limits_{j = 2}^k {({{\varepsilon }_{i}} - {{\varepsilon }_{j}}) + ({{\varepsilon }_{i}} - (1 - \varepsilon ))} \approx - 1 < 0.$
Таким образом, правая часть для $i$ велика по сравнению с правой частью для вершин, находящихся на плато уровня $\varepsilon $ или $1 - \varepsilon $. Поскольку $A{\mathbf{u}}$ в методе Ньютона (или другом методе линейной коррекции невязки) является частью функции $f$, и ${{A}_{{{\text{Del}}}}}$ вносит свой вклад в матрицу Якоби ($A$), то
$\delta \tilde {u}{{{\text{|}}}_{i}} = ( - A_{{{\text{Del}}}}^{{ - 1}}f){{{\text{|}}}_{i}} > 0$
и поэтому ${{u}^{ + }} = u + \delta u$ (следующее приближение ${{u}_{i}}$ в методе Ньютона) увеличивается диффузионной частью матрицы Якоби до среднего значения ближайших соседей.

В случае неудовлетворения условия Делоне, коэффициент $ - 1$ и огромная (по абсолютному значению) правая сторона эквивалентна сильному источнику с неправильным знаком. Кроме того, ${{A}_{{{\text{non}} - {\text{Del}}}}}$ может быть неограниченным. Очевидно, постоянные $A$ не могут быть причиной скачка, но другие составляющие могут создать скачки в случае небольшой диффузии. Небольшая ${{A}_{{{\text{non}} - {\text{Del}}}}}$ может быть уже достаточной для того, чтобы нарушить предельные значения $u$ и таким образом привести к серьезным проблемам, которые не могут быть решены простыми ограничителями для $u$, не создавая при этом новых проблем в зависимости от специфики задачи.

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

3.5. Непрерывная зависимость от положения вершин

Рассмотрим одномерную краевую задачу:

$\begin{gathered} u{\kern 1pt} '{\kern 1pt} '(y) + cu(y) = 0,\quad y \in (0,1),\quad c = {\text{const}}, \hfill \\ u(0) = 1, \hfill \\ u(1) = 0. \hfill \\ \end{gathered} $
Эту задачу и ее дискретизацию на заданной одномерной сетке по $y$ можно расширить на двумерную ортогональную сетку по $x$, используя произвольные значения шага сетки и однородные граничные условия Неймана. При этом сохраняется профиль $u$ по $y$ (двумерное решение будет инвариантным в направлении по $x$).

Триангуляция Делоне для получившегося регулярного прямоугольного распределения узлов не однозначна, потому что каждая из четырех соседних вершин одной ячейки в направлениях $x$ и $y$ лежат на одной окружности (фиг. 6, слева: оба направления диагональных ребер удовлетворяют условию Делоне). В данном случае однозначное положение диагональных ребер может получиться только при использовании специального правила, например, “снизу слева – наверх справа”. Использование точной арифметики для теста на выполнение условия Делоне в этом случае не поможет, так как он все равно приведет к произвольному выбору ребра в случае если четыре узла имеют точные координаты и лежат на одной окружности или по построению должны быть на одной окружности, но были сохранены с округлением в представлении с плавающей точкой. Например, поворот сетки в фиг. 6 на $\pi {\text{/}}4$ приведет к тому, что четыре узла, лежавшие вначале на одной окружности, в новой сетке не будут больше лежать на одной окружности, т.к. $sin\pi {\text{/}}4 = cos\pi {\text{/}}4 = \sqrt 2 {\text{/}}2$ не может быть представлен точно числами в арифметике с плавающей запятой без округления.

Фиг. 6.

Выбор диагональных ребер для сеток Делоне для регулярной прямоугольной сетки не однозначен, так как четыре вершины прямоугольника принадлежат к одной окружности. Жесткое правило для прямоугольной сетки, например “снизу слева – наверх справа”, позволяет простроить трансляционно-инвариантную сетку (слева). Некоторые диагональные ребра выбраны по-другому (справа), что приводит к так называемой перекрестной схеме вокруг ${{{\mathbf{v}}}_{1}}$.

Для VFVM выбор диагональных ребер неважен, потому что грань Вороного, соответствующая диагональному ребру, имеет площадь, равную нулю. Поверхности и объемы Вороного масштабируются в соответствии с размером шага по $x$. Таким образом, сохраняется неизменность решения в направлении $x$. Кроме того, небольшие искажения координат узлов приведут к соответственно небольшому изменению площадей граней (объемы ячеек и граней Вороного являются непрерывными функциями координат узлов). То же справедливо и для матриц типа ${{G}^{{\text{T}}}}G$ в любом $d$-мерном пространстве (например, матрицы масс и матрицы Дирихле). Следовательно, очень небольшое изменение положения узлов сетки приведет к очень небольшому изменению дискретного решения.

Для метода конечных элементов (МКЭ) ситуация иная, поскольку здесь топология сетки имеет значение. С одной стороны, для сеток на фиг. 6, матрица жесткости для дискретизации с помощью P1 (линейных) конечных элементов инвариантна по отношению к выбору диагональных ребер, так как в двухмерном случае матрица жесткости МКЭ идентична соответствующей ей матрице метода конечных объемов, которая зависит только от координат узлов. С другой стороны, матрица масс МКЭ зависит от топологии сетки, например,

${{M}_{{ii}}} = \frac{{2\left| {{{\omega }_{i}}} \right|}}{{(d + 1)(d + 2)}},$
где ${{\omega }_{i}}$ – подобласть, относящаяся к узлу ${{{\mathbf{v}}}_{i}}$, а $\left| {{{\omega }_{i}}} \right|$ – ее объем (площадь). Сетка на фиг. 6 (слева) – трансляционно-инвариантна (исключая первый и последний узел), и каждая подобласть ${{\omega }_{i}}$ имеет одинаковый размер. Следовательно, справедливо
${{M}_{{11}}} = {{M}_{{22}}} = \frac{{2 \times \left( {6 \times \tfrac{{{{h}^{2}}}}{2}} \right)}}{{12}} = \frac{1}{2}{{h}^{2}}.$
Сетка на фиг. 6 (справа), хотя и удовлетворяет условию Делоне, но некоторые диагональные ребра были выбраны по другому. Подобласти, принадлежащие ${{{\mathbf{v}}}_{1}}$ и ${{{\mathbf{v}}}_{2}}$, имеют разные размеры, так же, как и соответствующие им элементы в матрице масс:
${{M}_{{11}}} = \frac{{2 \times \left( {8 \times \tfrac{{{{h}^{2}}}}{2}} \right)}}{{12}} = \frac{2}{3}{{h}^{2}},\quad {\text{тогда}}\;{\text{как}}\quad {{M}_{{22}}} = \frac{{2 \times \left( {4 \times \tfrac{{{{h}^{2}}}}{2}} \right)}}{{12}} = \frac{1}{3}{{h}^{2}},$
так что матрица масс не будет трансляционно-инвариантной, в отличие от матрицы жесткости. Как следствие, численное решение потеряет инвариантность в направлении $x$. Более того, несогласованность матрицы масс также может привести и к более низкому, чем ожидаемому, порядку сходимости. Одноточечная квадратура для построения матрицы массы не исправит этот недостаток в общем, потому что распределение массы по узлам зависит от топологии сетки и, таким образом, уже является несогласованной еще до того, как будет проведено вычисление матрицы масс. Подробное обсуждение потери точности линейных МКЭ, а также конечных элементов более высокого порядка в зависимости от топологии сеток содержится в недавних работах Коптевой [22], [23].

4. ПРИМЕРЫ ПРИМЕНЕНИЯ

4.1. Разделение фаз

Рассмотрим диффундирующие черные и белые частицы, сумма которых везде равна 1 (статистика Ферми), короткодействующий потенциал взаимодействия (черное притягивает черное) и диффузию, действующую как расталкивающая сила (такова качественная интерпретация известного соотношения Эйнштейна). Градиенты можно увеличить за счет увеличения потенциала притягивания по отношению к диффузии.

Концентрации ${{u}_{b}}$ черных и ${{u}_{w}}$ белых частиц удовлетворяют соотношению

$1 > 1 - \varepsilon > {{u}_{b}} > \varepsilon > 0,\quad {{u}_{w}} = 1 - {{u}_{b}}.$
Если доминирует диффузионное расталкивание, то решением является идеальное распределение серого. В противном случае, черные и белые частицы будут разделены минимальной поверхностью.

Ниже приводится краткое математическое описание этой нелокальной задачи фазовой сегрегации; подробное обсуждение дано в статье [10] и в литературе к ней. Задача описывается следующей системой уравнений:

(1)
$\begin{gathered} \tfrac{{\partial u}}{{\partial t}} - \nabla \cdot \left( {f\left( {{\text{|}}\nabla v{\text{|}}} \right)\left( {\nabla u + \tfrac{{\nabla w}}{{\phi {\kern 1pt} '{\kern 1pt} '(u)}}} \right)} \right) = 0\quad {\text{в}}\quad \Omega , \hfill \\ u(0, \cdot ) = {{u}_{0}}( \cdot ), \hfill \\ v = \phi {\kern 1pt} '(u) + w, \hfill \\ w(t,{\mathbf{x}}) = \int\limits_\Omega {\mathcal{K}\left( {{\text{|}}{\mathbf{x}} - {\mathbf{y}}{\text{|}}} \right)\left( {1 - 2u(t,{\mathbf{y}})} \right)d{\mathbf{y}}} , \hfill \\ \end{gathered} $
где $\Omega \subset {{\mathbb{R}}^{N}}$, $1 \leqslant N \leqslant 3$ – ограниченная Липшицева область, $\phi $ – выпуклая функция, ядро $\mathcal{K}$ представляет нелокальные силы притяжения, а $w$ и $v$ – потенциал взаимодействия и химический потенциал соответственно. Первоначальное значение ${{u}_{0}}$ удовлетворяет неравенству $0 < {{u}_{0}} < 1$. Система дополняется однородными граничными условиями Неймана.

Эта система была выведена в [16] для случая постоянной функции $f$ и может быть рассмотрена как нелокальный вариант уравнения Кан-Хиллиарда, связанного с локальной свободной энергией Гинзбурга-Ландау

${{F}_{{gl}}}(u) = \int\limits_\Omega {\left( {\phi (u) + \kappa u\left( {1 - u} \right) + \frac{\lambda }{2}{{{\left| {\nabla u} \right|}}^{2}}} \right)} \,d{\mathbf{x}}.$
Функционалом Ляпунова для (1) является
$F(u) = \int\limits_\Omega {\left( {\phi (u) + u\int\limits_\Omega {\mathcal{K}\left( {\left| {x - y} \right|} \right)\left( {1 - u(y)} \right)dy} } \right)\,} d{\mathbf{x}}.$
Рассмотрим специальную модель, в которой функция $f$ постоянна, $\phi (u) = ulnu + (1 - u)ln(1 - u)$, и $\mathcal{K}$ – это функция Грина для эллиптической краевой задачи
$\begin{gathered} - {{\sigma }^{2}}\nabla \cdot \nabla w + w = m\left( {1 - 2u} \right)\quad {\text{в}}\quad \Omega , \hfill \\ \nu \cdot \nabla w = 0\quad {\text{на}}\quad \partial \Omega \hfill \\ \end{gathered} $
с константами $\sigma ,m > 0$. При расчете на малые времена эта модель также используется для реконструкции изображений и подавления шума, и не ограничивается только двумя видами частиц [11], [12], [18].

Подходящим численным методом решения этой проблемы является схема

$\begin{gathered} ({{G}^{T}}{{\sigma }^{2}}G + [V]){{{\mathbf{w}}}^{ + }} = [V]m({\mathbf{1}} - 2{{{\mathbf{u}}}^{ + }}), \\ [V]\frac{{{\mathbf{u}} - {{{\mathbf{u}}}^{ + }}}}{\tau } = {{G}^{T}}f\left( {G{{u}^{ + }} + \left[ {\frac{1}{{\Phi _{a}^{{''}}}}} \right]\frac{1}{2}G({{{\mathbf{w}}}^{ + }} + {\mathbf{w}})} \right), \\ \end{gathered} $
Доказано, что данная схема является диссипативной [10, теорема 2.7]. Более того, если ${{{\mathbf{v}}}^{ + }} \ne c$, $0 < u < 1$, то неявный метод Эйлера сохраняет свойство $0 < {{u}^{ + }} < 1$ [10, утверждение 2.9].

Контроль размера шага по времени в вычислениях основан на разностях типа предиктор-корректор для свободной энергии и дискретной скорости диссипации. Матрицы Якоби выводятся при помощи аналитического дифференцирования на каждом симплексе, а метод Ньютона применяется без демпфирования. В качестве начального решения берется

${{u}_{0}}({\mathbf{x}}) = \varepsilon + \frac{{{{x}_{1}}}}{{128\left( {1 + 2\varepsilon } \right)}},\quad \Omega = (0,128) \times (0,128),\quad m = 8,\quad \varepsilon = {{10}^{{ - 7}}},\quad \sigma = \sqrt f = 2.$

Диссипативная схема, как и ожидалось, выдает приближение минимальной поверхности. После ${{10}^{{20}}}$ шагов по времени исходное распределение серого преображается, проходя через узор из полосок с пузырьками, в черные и белые половинки квадрата, разделенные прямой вертикальной линией посередине квадрата. Фигуры 7 и 8 показывают типичный процесс эволюции: тонкие полосы (фиг. 7, слева) распадаются на пузырьки с диаметром, связанным с шириной полосы (фиг. 7, в центре). Симметрия нарушается из-за ошибок округления, небольших неточностей в исходных данных, а также из-за того, что область $\Omega $ имеет большой размер по сравнению с $\sigma $ и не совместима с произвольной шириной полос. Самая длинная шкала времени определяется временем, необходимым для выравнивания границы между этими половинами до прямой линии.

Фиг. 7.

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

Фиг. 8.

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

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

4.2. DePFET детекторы с активными пикселями

Детекторы DePFET с активными пикселями – это кремниевые сенсоры на базе объединенных p-канальных полевых транзисторов (DePFETs) [20], [26], которые представляют собой полевые транзисторы (FETs), интегрированные в объединенный кремниевый объем. Входящие частицы (фотоны) генерируют электроны в объеме, которые затем собираются на внутреннем затворе под транзисторным каналом. Это увеличивает проводимость транзистора пропорционально количеству накопленного заряда, обеспечивая этим измерение энергии фотонов. Уникальными преимуществами детекторов, основанных на этой технологии, являются усиление и хранение заряда сигнала непосредственно внутри пикселей, неразрушающее считывание, очень низкий уровень шума и высокая энергоэффективность, что делает их идеальными для рентгеновской астрономии и физики частиц [26], [34].

Примером такого детектора является широкоугольная камера (Wide Field Imager, WFI) для усовершенствованного телескопа для высокоэнергетической астрофизики (Advanced Telescope for High-Energy Astrophysics, Athena) [1], [2] Европейского Космического Агентства. WFI – это твердотельный рентгеновский спектрометр для рентгеновских снимков в полосе 0.1–15 кВ с одновременным спектральным и временным разрешением фотонов. Он разработан в совместной работе Лаборатории Полупроводников Общества Макса-Планка (MPG HLL), Эрлангенского Центра Физики Астрочастиц, Института Астрономии и Астрофизики Тюбингена, Дрезденского Технического Университета, PNSensor GmbH, Лейстерского Университета и Научно-исследовательского Института Астрофизики и Планетологии в Тулузе [2]. Численное моделирование проводилось в сотрудничестве MPG HLL с Институтом Вейерштрасса по Прикладному Анализу и Стохастике.

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

$\begin{gathered} - \nabla \cdot \varepsilon \nabla w = C - n - p, \\ \frac{{\partial n}}{{\partial t}} - \nabla \cdot {{n}_{i}}{{\mu }_{n}}{{e}^{w}}\nabla {{e}^{{ - {{\phi }_{n}}}}} = R, \\ \frac{{\partial p}}{{\partial t}} - \nabla \cdot {{n}_{i}}{{\mu }_{p}}{{e}^{{ - w}}}\nabla {{e}^{{{{\phi }_{p}}}}} = R, \\ \end{gathered} $
где $w$ – электростатический потенциал, $n = {{n}_{i}}{{e}^{{w - {{\phi }_{n}}}}}$ и $p = {{n}_{i}}{{e}^{{{{\phi }_{p}} - w}}}$ – плотности электронов и дырок, $C$ – легирование, ${{\phi }_{n}}$ и ${{\phi }_{p}}$ – квазиуровни Ферми, и $R = r({\mathbf{x}},n,p)(n_{i}^{2} - np)$ описывает рекомбинацию/генерацию ($r(x,n,p) > 0$).

Для численного моделирования важно, чтобы дискретная система сохраняла существенные свойства аналитической системы, такие как принцип максимума, положительную концентрацию носителей и сохранение тока. VFVM позволяет перенести эти свойства в дискретную модель. Численные расчеты производились по хорошо известной схеме Шарфеттера-Гуммеля (подробное описание находится, например, в [14]). Фигура 9 показывает моделирование распространения фронта деплеции в пикселе детектора. Ввиду симметричности области, вычисления выполнены на половине пикселя.

Фиг. 9.

Детектор Athena WFI: моделирование распространения фронта деплеции.

5. ЗАКЛЮЧЕНИЕ

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

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

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

  1. Athena – Advanced Telescope for High-Energy Astrophysics. https://sci.esa.int/athena

  2. Athena. The extremes of the universe: from black holes to large-scale structure. Assessment study report ESA/SRE (2011)17, European Space Agency, 2012. http://sci.esa.int/jump.cfm?oidI835

  3. Allen D., Southwell R. Relaxation methods applied to determine the motion, in two dimensions, of a viscous fluid past a fixed cylinder // Quart. J. Mech. and Appl. Math. 1955. V. 8. P. 129–145.

  4. Bank R.E., Rose D.J. Some error estimates for the box method // SIAM J. Numer. Anal. 1987. V. 24 (4). P. 777–787.

  5. Chew P.L. Constrained Delaunay triangulations // Algorithmica. 1989. V. 4 (1). P. 97–108.

  6. Delaunay B. Sur la sphére vide // Izvestia Akademii Nauk SSSR. Otd. Matem. i Estestv. Nauk. 1934. V. 6. P. 793–800.

  7. Eymard R., Gallouet T., Herbin R. Finite volume methods. P. G. Ciarlet and J. L. Lions, Ed., Handbook of Numerical Analysis. Elsevier Science B.V., Amsterdam, 1997.

  8. Fleischmann P. Mesh Generation for Technology CAD in Three Dimensions. Dissertation, Technische Universität Wien, 1999.

  9. Gabriel K.R., Sokal R.R. A new statistical approach to geographic variation analysis // Syst. Biol. 1969. V. 18 (3). P. 259–278.

  10. Gajewski H., Gärtner K. A dissipative discretization scheme for a nonlocal phase segregation model // Z. Angew. Math. Mech. 2005. V. 85. P. 815–822.

  11. Gajewski H., Gärtner K. On a nonlocal model of image segmentation // Z. Angew. Math. Phys. 2005. V. 56. P. 572–591.

  12. Gajewski H., Griepentrog J.A. A descent method for the free energy of multicomponent systems // Discrete Contin. Dyn. Syst. 2006. V. 15 (2). P. 505–528.

  13. Garanzha V., Kudryavtseva L., Tsvetkova V. Structured orthogonal near-boundary Voronoi mesh layers for planar domains. In Numerical Geometry, Grid Generation and Scientific Computing, Lect. Notes Comput. Sci. Eng. 131. Springer International Publishing, 2019. Proceedings of the 9th International Conference, NUMGRID 2018 / Voronoi 150, Celebrating the 150th Anniversary of G.F. Voronoi, Moscow, Russia, December 2018.

  14. Gärtner K. Existence of bounded discrete steady-state solutions of the van Roosbroeck system on boundary conforming Delaunay grids // SIAM J. Sci. Comput. 2009. V. 31 (2). P. 1347–1362.

  15. Gärtner K. Existence of bounded discrete steady state solutions of the van Roosbroeck system with monotone Fermi–Dirac statistic functions // J. Comput. Electron. 2015. V. 14 (3).

  16. Giacomin G., Lebowitz J.L. Phase segregation dynamics in particle systems with long range interactions. I. Macroscopic limits // J. Statist. Phys. 1997. V. 87 (1–2). P. 37–61.

  17. Glitzky A., Gärtner K. Energy estimates for continuous and discretized electro-reaction-diffusion systems // Nonlinear Anal. 2009. V. 70. P. 788–805.

  18. Griepentrog J.A. On the unique solvability of a non-local phase separation problem for multicomponent systems. In Nonlocal elliptic and parabolic problems, Banach Center Publ. Polish Acad. Sci. Inst. Math., Warsaw, 2004. V. 66. P. 153–164.

  19. Il’in A.M. A difference scheme for a differential equation with a small parameter multiplying the second derivative // Mat. Zametki. 1969. V. 6. P. 237–248.

  20. Kemmer J., Lutz G. New detector concepts // Nucl. Instrum. Methods Phys. Res. A. 1987. V. 253 (3). P. 365–377.

  21. Kerkhoven T. Piecewise linear Petrov-Galerkin error estimates for the box method // SIAM J. Numer. Anal. 1996. 33 (5). P. 1864–1884.

  22. Kopteva N. Linear finite elements may be only first-order pointwise accurate on anisotropic triangulations // Math. Comp. 2014. V. 83 (289). P. 2061–2070.

  23. Kopteva N. How accurate are finite elements on anisotropic triangulations in the maximum norm? // J. Comput. Appl. Math. 2020. V. 364. P. 112316.

  24. Kosik R., Fleischmann P., Haindl B., Pietra P., Selberherr S. On the interplay between meshing and discretization in three-dimensional diffusion simulation // IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst. 2000. V. 19 (11). P. 1233–1240.

  25. Lawson C.L. Software for ${{C}^{1}}$ surface interpolation. In Mathematical Software III. 1977. P. 161–194.

  26. Lutz G., Andricek L., Eckardt R., Hälker O., Hermann S., Lechner P., Richter R., Schaller G., Schopper F., Soltau H., Strüder L., Treis J., Wölfl S., Zhang C. DePFET–detectors: New developments // Nucl. Instrum. Methods Phys. Res. A. 2007. V. 572 (1). P. 311–315.

  27. MacNeal R.H. An asymmetrical finite difference network // Quart. Math. Appl. 1953. V. 11. P. 295–310.

  28. Rippa S. Minimal roughness property of the Delaunay triangulation // Comput. Aided Geom. Design. 1990. V. 7 (6). P. 489–497.

  29. Scharfetter D.L., Gummel H.K. Large-signal analysis of a silicon read diode oscillator // IEEE Trans. Electr. Dev. 1969. V. 16. P. 64–77.

  30. Shewchuk J.R. Delaunay refinement algorithms for triangular mesh generation // Comput. Geom. 2002. V. 22 (1–3). P. 21–74.

  31. Shewchuk J.R. General-dimensional constrained Delaunay and constrained regular triangulations. I. Combinatorial properties // Discrete Comput. Geom. 2008. V. 39 (1–3). P. 580–637.

  32. Si H. TetGen, a Delaunay-based quality tetrahedral mesh generator // ACM Trans. Math. Softw. 2015. V. 41 (2). P. 11:1–11:36.

  33. Si H., Gärtner K., Fuhrmann J. Boundary conforming Delaunay mesh generation // Comput. Math. Math. Phys. 2010. V. 50 (1). P. 38–53.

  34. Treberspurg W., Andritschke R., Hauser G., Lechner P., Meidinger N., Ninkovic J., Müller-Seidlitz J., Schopper F. Measurement results of different options for spectroscopic X-ray DePFET sensors // J. Instrum. 2018. V. 13 (9). P. P09014–P09014.

  35. Varga R.S. Matrix Iterative Analysis. Prentice-Hall, 1962. Englewood Cliffs, N. J.

  36. Voronoi G. Nouvelles applications des paramètres continus à la théorie des formes quadratiques // Reine Angew. Math. 1907. V. 133. P. 97–178.

  37. Xu J., Zikatanov L. A monotone finite element scheme for convection-diffusion equations // Math. Comp. 1999. V. 68 (228). P. 1429–1446.

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

Инструменты

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

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

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