Журнал вычислительной математики и математической физики. T. 63, Номер 10, 2023

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

Журнал вычислительной математики и математической физики, 2023, T. 63, № 10, стр. 1591-1599

Сеточно-характеристический численный метод на нерегулярной расчетной сетке с расширением шаблона интерполяции

А. В. Васюков 1, *, И. Е. Смирнов 1, **

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

* E-mail: a.vasyukov@phystech.edu
** E-mail: smirnov.ie@phystech.edu

Поступила в редакцию 21.03.2023
После доработки 27.05.2023
Принята к публикации 26.06.2023

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

Аннотация

В работе предложен сеточно-характеристический численный метод для решения многомерного уравнения переноса на неструктурированной расчетной сетке с порядком выше первого без использования вспомогательных точек на ребрах и гранях. Отсутствие вспомогательных точек на ребрах и гранях позволяет упростить топологию расчетной сетки при ее движении, что актуально при решении динамических задач механики деформируемого твердого тела. Для повышения порядка аппроксимации в работе используется аналог расширения сеточного шаблона, реализованный для неструктурированной сетки. В работе приведены результаты тестирования предложенной численной схемы для непрерывно дифференцируемых, непрерывных, разрывных решений. Библ. 14. Фиг. 7. Табл. 5.

Ключевые слова: сеточно-характеристический метод, неструктурированная сетка, уравнение переноса, динамическая задача, механика деформируемого твердого тела.

1. ВВЕДЕНИЕ

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

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

Определяющая система уравнений в частных производных, описывающая динамические процессы в деформируемом твердом теле в трехмерном случае, имеет следующий характерный вид [4]:

(1)
$\frac{{\partial {\mathbf{w}}}}{{\partial t}} + {{A}_{x}} \cdot \frac{{\partial {\mathbf{w}}}}{{\partial x}} + {{A}_{y}} \cdot \frac{{\partial {\mathbf{w}}}}{{\partial y}} + {{A}_{z}} \cdot \frac{{\partial {\mathbf{w}}}}{{\partial z}} = 0;$
здесь ${\mathbf{w}}$ – вектор неизвестных, ${{A}_{x}},\;{{A}_{y}},\;{{A}_{z}}$ – постоянные матрицы. Точный вид ${\mathbf{w}},\;{{A}_{i}}$ зависит от конкретной задачи, в том числе от используемых реологических моделей.

Известно [5], [6], что для многих практически значимых постановок существует аналитическое разложение матриц в следующем виде:

${{A}_{i}} = \Omega _{i}^{{ - 1}} \cdot {{\Lambda }_{i}} \cdot {{\Omega }_{i}}.$

Для решения многомерной системы уравнений можно использовать расщепление по направлениям [7], [8]. Точный вид схемы расщепления, использованной в данной работе, приведен ниже при описании численного метода. Но для любой схемы расщепления для пространственной производной по одному из направлений возникает соотношение:

$\frac{{\partial {\mathbf{w}}}}{{\partial t}} + {{\Omega }^{{ - 1}}} \cdot {{\Lambda }_{x}} \cdot \Omega \cdot \frac{{\partial {\mathbf{w}}}}{{\partial x}} = 0,$
откуда следует

$\Omega \cdot \frac{{\partial {\mathbf{w}}}}{{\partial t}} + {{\Lambda }_{x}} \cdot \Omega \cdot \frac{{\partial {\mathbf{w}}}}{{\partial x}} = 0.$

После введения обозначения ${\mathbf{u}} = \Omega \cdot {\mathbf{w}}$, исходная система распадается на независимые уравнения вида:

$\frac{{\partial {{u}_{i}}}}{{\partial t}} + {{\lambda }_{{{{i}_{x}}}}} \cdot \frac{{\partial {{u}_{i}}}}{{\partial x}} = 0.$

В силу изложенного в рамках данной работы не рассматривается точный вид ${\mathbf{w}},\;{{A}_{i}}$, а численно решается уравнение переноса, которое лежит в основе полной численной схемы для многомерных динамических задач механики деформируемого твердого тела.

При практической реализации численных схем частым их желательным свойством является обеспечение разумно высокого порядка аппроксимации без необходимости использования большого сеточного шаблона. Одним из традиционных подходов к этой задаче является использование компактных продолженных схем [9]. Для сеточно-характеристического метода на структурированных расчетных сетках данный подход также применим [10], [11]. Однако для неструктурированной расчетной сетки такой подход на данный момент не реализован.

При выполнении расчетов для иженерных объектов сложной формы область интегрирования описывается, как правило, именно с использованием неструктурированной расчетной сетки, в двухмерном случае – из треугольников, в трехмерном – из тетраэдров [12]. Таким образом, требуется численно решать уравнения переноса на данной сетке с высоким порядком аппроксимации.

Очевидным образом на треугольниках и тетраэдрах можно обеспечить интерполяцию с первым порядком для значений в произвольной точке по значениям в вершинах. Традиционным способом повышения порядка является введение дополнительных точек на ребрах и гранях элементов сетки [13], [14]. Этот подход хорошо показал себя в решении многих задач. Однако он также содержит и некоторые недостатки. Так, при расчетах задач прочности возникает необходимость описывать значительные деформации объекта, что требует перемещения узлов расчетной сетки. В этом случае использование вспомогательных расчетных узлов, которые жестко зафиксированы на ребрах и гранях, приводит к значительным сложностям. Либо расчет перемещений данных узлов де-факто не выполняется, а их движение интерполируется по узлам в вершинах. Либо при движении узлов на ребрах и гранях существенно нарушается исходная топология элементов расчетной сетки, что резко усложняет построение численной схемы для нее.

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

2. ПОСТАНОВКА ЗАДАЧИ И ЧИСЛЕННЫЙ МЕТОД

Рассматривается численное решение двумерного уравнения переноса на двумерных нерегулярных расчетных сетках. Решаемое уравнение для функции $u(x,y,t)$ в квадрате $[ - 1,1] \times [ - 1,1]$ с периодическими граничными условиями имеет следующий вид:

$\begin{gathered} {{u}_{t}} + {{\lambda }_{x}}{{u}_{x}} + {{\lambda }_{y}}{{u}_{y}} = 0, \\ {{\left. {u(x,y,t)} \right|}_{{t = 0}}} = F(x,y), \\ \end{gathered} $
(2)
$u(a \cdot {{T}_{x}} + x,b \cdot {{T}_{y}} + y,t) = u(x,y,t),$
$\begin{gathered} {{T}_{x}} = {{T}_{y}} = 2;\quad a,b \in Z, \\ x \in [ - 1,1],\quad y \in [ - 1,1]. \\ \end{gathered} $
Здесь $F(x,y)$ – функция начальных условий; ${{T}_{x}},\;{{T}_{y}}$ – период функции $u(x,y,t)$; ${{\lambda }_{x}},\;{{\lambda }_{y}}$ – скорости по соответствующим направлениям.

Для дискретизации уравнения по времени используется равномерная сетка по времени с шагом $\tau $, такая что $T = N \cdot \tau $. Здесь число $N$ – количество слоев по времени. На каждом слое по времени используется неравномерная сетка из треугольников, построенная при помощи алгоритма Делоне, пример расчетной сетки показан на фиг. 1.

Фиг. 1.

Пример расчетной сетки.

При численном решении задача расщепляется по пространственным переменным на два независимых уравнения, решаемых последовательно на каждой временной итерации:

$\begin{gathered} \frac{{\partial u}}{{\partial t}} + {{\lambda }_{x}} \cdot \frac{{\partial u}}{{\partial x}} = 0, \\ \frac{{\partial u}}{{\partial t}} + {{\lambda }_{y}} \cdot \frac{{\partial u}}{{\partial y}} = 0. \\ \end{gathered} $

К значениям на текущем временном слое применяется оператор, соответствующий пространственной производной вдоль оси $OX$, затем к результату применяется оператор, соответствующий пространственной производной вдоль оси $OY$. Результат применения второго оператора является значением на новом временном слое. Данная схема имеет второй порядок аппроксимации.

Для решения каждого из получившихся одномерных уравнений переноса используется сеточно-характеристический метод. Рассмотрим метод на примере уравнения для $x$ координаты относительно $u(x,y,t)$:

$\frac{{\partial u}}{{\partial t}} + {{\lambda }_{x}} \cdot \frac{{\partial u}}{{\partial x}} = 0.$

В рамках данного метода дифференциальное уравнение первого порядка в частных производных сводится к обыкновенному дифференциальному уравнению вдоль характеристики:

$\frac{{dx}}{{dt}} = {{\lambda }_{x}}.$

Тогда верно следующее:

$\frac{{\partial u}}{{\partial t}} + {{\lambda }_{x}} \cdot \frac{{\partial u}}{{\partial x}} = 0 \to \frac{{\partial u}}{{\partial t}} + \frac{{dx}}{{dt}} \cdot \frac{{\partial u}}{{\partial x}} = 0 \to \frac{{dw}}{{dt}} = 0,$
т.е. вдоль характеристики $dx{\text{/}}dt = {{\lambda }_{x}}$ решение не зависит от времени:

${{\left. {u(x,y,t)} \right|}_{{dx/dt = {{\lambda }_{x}}}}} = u(x,y).$

Таким образом, для нахождения значения функции в момент времени ${{t}^{{n + 1}}}$ из заданной точки опускается характеристика (прямая, задаваемая уравнением $x = {{\lambda }_{x}} \cdot t$) на предыдущий слой по времени ${{t}^{n}}$. В точке пересчения $({{x}_{0}},\;{{y}_{0}},\;{{t}^{n}})$ этой прямой с плоскостью $t = {{t}^{n}} = \tau \cdot n = {\text{const}}$ аппроксимируется значение функции, с использованием известных значений в точках на данном слое по времени. Далее это значение переносится в точку $({{x}_{0}} + {{\lambda }_{x}} \cdot \tau ,\;{{y}_{0}},\;{{t}^{{n + 1}}})$. Данный подход проиллюстрирован на фиг. 2.

Фиг. 2.

Перенос значения вдоль характеристики.

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

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

Для получения схемы второго порядка необходимо определить коэффициенты интерполяционного полинома

$P(x,y) = \alpha \cdot {{x}^{2}} + \beta \cdot {{y}^{2}} + \gamma \cdot xy + \delta \cdot x + \epsilon \cdot y + \zeta .$

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

$\begin{gathered} P({{x}_{1}},{{y}_{1}}) = {{u}_{1}},\quad P({{x}_{2}},{{y}_{2}}) = {{u}_{2}},\quad P({{x}_{3}},{{y}_{3}}) = {{u}_{3}}, \\ P({{x}_{4}},{{y}_{4}}) = {{u}_{4}},\quad P({{x}_{5}},{{y}_{5}}) = {{u}_{5}},\quad P({{x}_{6}},{{y}_{6}}) = {{u}_{6}}. \\ \end{gathered} $
Здесь ${{x}_{i}},{{y}_{i}}$ – координаты выбранных ближайших точек, ${{u}_{i}}$ – значения функции в них. Из решения данной системы находятся необходимые коэффициенты полинома $P$, после чего определяется значение функции в требуемой точке.

Для рассматриваемой модельной задачи (2) существует аналитическое решение:

${{u}_{{{\text{analitic}}}}}(x,y,t) = F(x - {{\lambda }_{x}} \cdot t,y - {{\lambda }_{y}} \cdot t).$

Численное решение связано с точным аналитическим решением следующим образом:

${{u}_{{{\text{numeric}}}}} = {{u}_{{{\text{analytic}}}}} + R(h).$
Здесь $R(h)$ – невязка, функция от мелкости пространственной сетки $h$.

Предполагается, что

$R(h) = {{h}^{p}} + o({{h}^{p}}).$
Таким образом, порядок аппроксимации $p$ может быть определен следующим образом:

$r = \left| {{{u}_{{{\text{numeric}}}}} - {{u}_{{{\text{analytic}}}}}} \right| = {{h}^{p}} + o({{h}^{p}}) \to \ln r{\text{/}}\ln h \simeq p.$

С учетом того, что сетка неструктурированная, понятие ее мелкости может быть трактовано разным образом. В данной работе в качестве $h$ используются:

• обратный характерный масштаб сетки ${{h}_{{{\text{scale}}}}} = 1{\text{/scale}}$, где scale – наибольший линейный размер ячейки;

• обратный корень количества точек сетки ${{h}_{{{\text{dots}}}}} = 1{\text{/}}\sqrt {{\text{dots}}\;{\text{number}}} $.

Для определения $r$ предполагается, что значения в точках на каждом слое по времени можно занумеровать, после чего используются следующие нормы:

• ${{e}_{1}} = {{e}_{\infty }} = \max ({\kern 1pt} \left| {{{u}_{{{\text{numeric}}}}}[i] - {{u}_{{{\text{analytic}}}}}[i]} \right|{\kern 1pt} )$,

• ${{e}_{2}} = \sum\nolimits_{i = 0}^N \left| {{{u}_{{{\text{numeric}}}}}[i] - {{u}_{{{\text{analytic}}}}}[i]{\kern 1pt} } \right|{\text{/}}N$,

• ${{e}_{3}} = \sqrt {\sum\nolimits_{i = 0}^N \,{{{({{u}_{{{\text{numeric}}}}}[i] - {{u}_{{{\text{analytic}}}}}[i])}}^{2}}} {\text{/}}N$.

Для определения фактического порядка сходимости решения строилось аналитическое решение в узлах сетки, затем вычислялись численные значения в этих же узлах. После этого рассчитывался вектор ошибки (невязки). По результатам расчета норм векторов ошибок на нескольких сетках с известными величинами шага по пространству $h$ с помощью метода наименьших квадратов строится прямая в координатах $(\ln r,\ln h)$, по ее наклону определяется порядок аппроксимации $p$.

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

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

Во всех экспериментах расчетная область представляла собой квадрат $[ - 1,1] \times [ - 1,1]$ с периодическими граничными условями. Количество cлоев по времени $N = 51$, безразмерное время расчета $T = 1$, постоянный шаг по времени $\tau = T{\text{/}}(N - 1) = 0.02$. Скорости распространения возмущений по осям $OX$ и $OY$ соответственно равны: ${{\lambda }_{x}} = - 2,$ ${{\lambda }_{y}} = 5$.

3.1. Гладкое решение

Начальное условие было задано в виде:

(3)
$F(x,y) = {{\cos }^{4}}(x \cdot \pi {\text{/}}2) \cdot {{\cos }^{4}}(y \cdot \pi {\text{/}}2).$

На фиг. 3 приведены результаты: слева показан общий вид решения, справа – график определения фактического порядка сходимости. Количественные значения для порядка аппроксимации представлены в табл. 1. Приведены результаты всех вычисляемых норм. В первом столбце в качестве шага по пространству используется ${{h}_{{{\text{dots}}}}}$, во втором ${{h}_{{{\text{scale}}}}}$.

Фиг. 3.

Результаты для начальных условий (3).

Таблица 1.  

Порядок аппроксимации для начальных условий (3)

e ${{h}_{{{\text{dots}}}}}$ ${{h}_{{{\text{scale}}}}}$
${{e}_{1}}$ 1.914 1.864
${{e}_{2}}$ 2.687 2.610
${{e}_{3}}$ 3.581 3.479

3.2. Быстро спадающая экспонента

Начальное условие было задано в виде:

(4)

На фиг. 4 приведены результаты: слева показан общий вид решения, справа – график определения фактического порядка сходимости. Количественные значения для порядка аппроксимации представлены в табл. 2.

Фиг. 4.

Результаты для начальных условий (4).

Таблица 2.  

Порядок аппроксимации для начальных условий (4)

e ${{h}_{{{\text{dots}}}}}$ ${{h}_{{{\text{scale}}}}}$
${{e}_{1}}$ 1.242 1.203
${{e}_{2}}$ 1.899 1.842
${{e}_{3}}$ 2.615 2.537

3.3. Конус

Начальное условие было задано в виде:

(5)

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

Фиг. 5.

Результаты для начальных условий (5).

Таблица 3.  

Порядок аппроксимации для начальных условий (5)

e ${{h}_{{{\text{dots}}}}}$ ${{h}_{{{\text{scale}}}}}$
${{e}_{1}}$ 0.888 0.861
${{e}_{2}}$ 1.793 1.739
${{e}_{3}}$ 2.482 2.408

3.4. Корень

Начальное условие было задано в виде:

(6)

На фиг. 6 приведены результаты: слева показан общий вид решения, справа – график определения фактического порядка сходимости. Количественные значения для порядка аппроксимации представлены в табл. 4.

Фиг. 6.

Результаты для начальных условий (6).

Таблица 4.  

Порядок аппроксимации для начальных условий (6)

e ${{h}_{{{\text{dots}}}}}$ ${{h}_{{{\text{scale}}}}}$
${{e}_{1}}$ 0.500 0.487
${{e}_{2}}$ 1.439 1.397
${{e}_{3}}$ 1.993 1.935

3.5. Ступенька

Начальное условие было задано в виде:

(7)
$F(x,y) = \left( \begin{gathered} 1,\quad \max ({\kern 1pt} {\text{|}}x{\text{|}},\;{\text{|}}{\kern 1pt} y{\kern 1pt} {\text{|}}) \leqslant 0.5, \hfill \\ 0,\quad \max ({\kern 1pt} {\text{|}}x{\kern 1pt} {\text{|}},\;{\text{|}}y{\kern 1pt} {\kern 1pt} {\text{|}}) > 0.5. \hfill \\ \end{gathered} \right.$

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

Фиг. 7.

Результаты для начальных условий (7).

Таблица 5.  

Порядок аппроксимации для начальных условий (7)

e ${{h}_{{{\text{dots}}}}}$ ${{h}_{{{\text{scale}}}}}$
${{e}_{1}}$ –0.003 –0.002
${{e}_{2}}$ 0.956 0.929
${{e}_{3}}$ 1.477 1.435

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

В настоящей статье предложена версия сеточно-характеристического численного метода для уравнения ${{\partial }_{t}}u + {{\lambda }_{x}}{{\partial }_{x}}u + {{\lambda }_{y}}{{\partial }_{y}}u = 0$ на нерегулярной расчетной сетке. Новизна подхода заключается в том, что для повышения порядка аппроксимации используется аналог расширения сеточного шаблона, реализованный в данном случае для неструктурированной сетки. Предложенная схема обеспечивает порядок аппроксимации выше первого без использования вспомогательных точек на ребрах и гранях элементов сетки.

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

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

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

  1. Беклемышева К.А., Васюков А.В., Голубев В.И., Петров И.Б. Численное моделирование воздействия сейсмической активности на подводный композитный трубопровод // Матем. моделирование. 2019. Т. 31. № 1. С. 103–113.

  2. Беклемышева К.А., Петров И.Б. Моделирование разрушения гибридных композитов под действием низкоскоростного удара // Матем. моделирование. 2018. Т. 30. № 11. С. 27–43.

  3. Беклемышева К.А., Васюков А.В., Петров И.Б. Численное моделирование динамических процессов в биомеханике сеточно-характеристическим методом // Ж. вычисл. матем. и матем. физ. 2015. Т. 55. № 8. С. 1380–1390.

  4. Магомедов К.М., Холодов А.С. Сеточно-характеристические численные методы: учебное пособие для бакалавриата и магистратуры. М.: Издательство Юрайт, 2019. 313 с.

  5. Челноков Ф.Б. Явное представление сеточно-характеристических схем для уравнений упругости в двумерном и трехмерном пространствах // Матем. моделирование. 2006. Т. 18. № 6. С. 96–108.

  6. Челноков Ф.Б. Численное моделирование деформационных процессов в средах со сложной структурой. Дис. … канд. физ.-матем. наук. М.: МФТИ, 2005.

  7. Федоренко Р.П. Введение в вычислительную физику. М.: Издательство МФТИ, 1994. 528 с.

  8. Петров И.Б., Холодов А.С. Численное исследование некоторых динамических задач механики деформируемого твердого тела сеточно-характеристическим методом // Ж. вычисл. матем. и матем. физ. 1984. Т. 24. № 5. С. 722–739.

  9. Рогов Б. В., Михайловская М. Н. Монотонные бикомпактные схемы для линейного уравнения переноса // Матем. моделирование. 2011. Т. 23. № 6. С. 98–110.

  10. Голубев В.И., Петров И.Б., Хохлов Н.И. Компактные сеточно-характеристические схемы повышенного порядка точности для трехмерного линейного уравнения переноса // Матем. моделирование. 2016. Т. 28. № 2. С. 123–132.

  11. Khokhlov N.I., Petrov I.B. On one class of high-order compact grid-characteristic schemes for linear advection // Russian Journal of Numerical Analysis and Mathematical Modelling. 2016. T. 31. № 6. C. 355–368.

  12. Васюков А.В., Петров И.Б. Использование сеточно-характеристического метода на неструктурированных сетках из тетраэдров с большими топологическими неоднородностями // Ж. вычисл. матем. и матем. физ. 2018. Т. 58. № 8. С. 62–72.

  13. Агапов П.И., Челноков Ф.Б. Сравнительный анализ разностных схем для численного решения двумерных задач механики деформируемого твердого тела // Моделирование и обработка информации. М.: МФТИ. 2003. С. 19–27.

  14. Петров И.Б., Фаворская А.В. Библиотека по интерполяции высоких порядков на неструктурированных треугольных и тетраэдральных сетках // Информационные технологии. 2011. № 9. С. 30–32.

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

Инструменты

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

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

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