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

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

А. В. Сетуха^1, *, А. В. Семенова^1, **

¹ МГУ
119234 Москва, Ленинские горы, Россия

^* E-mail: setuhaav@rambler.ru
^** E-mail: anastasia.semenova@cs.msu.ru

Поступила в редакцию 30.05.2018
После доработки 23.01.2019
Принята к публикации 08.02.2018

DOI: 10.1134/S0044466919060127

Аннотация

В статье рассматривается линейное гиперсингулярное интегральное уравнение на поверхности (замкнутой или разомкнутой с краем), возникающее при решении краевой задачи Неймана для уравнения Лапласа методом граничных интегральных уравнений, если представить решение в виде потенциала двойного слоя. Для такого уравнения строится численная схема, основанная на триангуляции поверхности, аппроксимации решения кусочно-линейной функцией и применении метода коллокации в вершинах треугольников, аппроксимирующих поверхность. В результате возникает система линейных уравнений, коэффициенты которой выражаются через интегралы по ячейками разбиения, содержащие произведения базисных функций на ядро с сильной особенностью. В статье получены аналитические формулы для нахождения этих коэффициентов. Это потребовало вычисления указанных интегралов, выполнив для каждого интеграла обход окрестности особой точки так, чтобы в результате система линейных уравнений аппроксимировала интегралы от неизвестной функции в точках коллокации в смысле конечного значения по Адамару. Проведено тестирование метода на модельных задачах. Библ. 17. Фиг. 9.

Ключевые слова: численные методы, интегральные уравнения, гиперсингулярные интегралы, метод граничных элементов, теория потенциала.

1. ВВЕДЕНИЕ

Численные методы, основанные на решении такого уравнения, возникают, в частности, при решении задач аэродинамики панельными (вихревыми) методами. В так называемом методе “вихревых рамок”, возникшем первоначально из эмпирических соображений, при учете граничного условия по существу решается рассматриваемое интегральное уравнение на границе течения с применением кусочно-постоянной аппроксимации неизвестной функции и принципа коллокации ([1], [2, c. 473–477]). В дальнейшем методы такого типа были распространены на другие классы задач математической физики, в которых возникают аналогичные гиперсингулярные интегральные уравнения на граничной поверхности с сильной особенностью, раскрываемой в смысле конечного значения по Адамару (см. [3]–[5]).

Ряд результатов по сходимости метода кусочно-постоянных аппроксимаций и коллокаций для гиперсингулярных интегральных уравнений рассматриваемого типа был получен в работах [6], [7], однако только для случая, когда уравнение решается на плоском экране и используется равномерное разбиение на прямоугольные ячейки.

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

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

Однако вопросы конструктивной реализации разработанного метода в работе [12] не рассматривались. Этим вопросам посвящена настоящая статья. Ниже получены аналитические формулы для коэффициентов возникающей системы линейных уравнений. Эти коэффициенты выражаются через интегралы по ячейками разбиения, содержащие произведения базисных функций на ядро с сильной особенностью. При этом каждый интеграл вычисляется с обходом окрестности особой точки так, чтобы в результате система линейных уравнений аппроксимировала интегралы от неизвестной функции в точках коллокации в смысле конечного значения по Адамару. Произведено тестирование численного метода на модельных примерах.

2. ВИД РАССМАТРИВАЕМОГО УРАВНЕНИЯ

Пусть $\Sigma $ – простая гладкая поверхность, которая может быть замкнутой или разомкнутой с краем $\partial \Sigma $ (край есть замкнутая кусочно-гладкая кривая). Рассматривается интегральное уравнение относительно функции $g(y)$:

(1)

$(Ig)(x) = \int\limits_\Sigma {g(y)\frac{{{{\partial }^{2}}F(x - y)}}{{\partial {{n}_{x}}\partial {{n}_{y}}}}d{{\sigma }_{y}} = f(x)} ,\quad x \in \Sigma ,\quad F(x - y) = \frac{1}{{4\pi {\text{|}}x - y{\text{|}}}},$

где $d{{\sigma }_{y}}$ – элемент площади, ${\mathbf{n}}(x)$ – единичная нормаль к поверхности $\Sigma $ в точке $x$, а $\partial {\text{/}}\partial {{n}_{x}}$ – производная по направлению вектора ${\mathbf{n}}(x)$, которая вычисляется как частная производная функции, зависящей от переменной $x$, при фиксированном значении переменной $y$. Подынтегральное выражение имеет особенность порядка ${\text{|}}x - y{{{\text{|}}}^{{ - 3}}}$ при стремлении $x$ к $y$. Интеграл в уравнении (1) понимается в смысле конечного значения по Адамару [6, с. 131]:

(2)

$\int\limits_\Sigma {\frac{{{{\partial }^{2}}F(x - y)}}{{\partial {{n}_{x}}\partial {{n}_{y}}}}} g(y)d{{\sigma }_{y}} = \mathop {lim}\limits_{\varepsilon \to 0} \left\{ {\int\limits_{\Sigma \backslash {{U}_{\varepsilon }}(x)} {\frac{{{{\partial }^{2}}F(x - y)}}{{\partial {{n}_{x}}\partial {{n}_{y}}}}} g(y)d{{\sigma }_{y}} - \frac{{g(x)}}{{2\varepsilon }}} \right\},$

где ${{U}_{\varepsilon }}(x)$ – трехмерная окрестность точки $x$ радиуса $\varepsilon $.

В случае разомкнутой поверхности на ее краю ставится условие

(3)

$g(y) = 0,\quad y \in \partial \Sigma .$

Записанное интегральное уравнение возникает при решении методом граничных интегральных уравнений краевой задачи Неймана для уравнения Лапласа в области $\Omega $, границей которой является поверхность $\Sigma $:

(4)

$\Delta u = 0\quad {\text{в }}\quad \Omega ,\quad \partial u{\text{/}}\partial n = f\quad {\text{н а }}\quad \Sigma .$

В случае, если поверхность $\Sigma $ замкнутая, считается, что $\Omega $ есть одна из областей (внешняя ${{\Omega }^{ + }}$ или внутренняя ${{\Omega }^{ - }}$), на которые поверхность $\Sigma $ делит трехмерное пространство. В случае, когда $\Sigma $ есть разомкнутая поверхность, область $\Omega $ есть все пространство вне этой поверхности (здесь мы считаем, что поверхность содержит свой край и поэтому является замкнутым множеством в трехмерном пространстве), граничное условие ставится на обеих сторонах поверхности $\Sigma $ (рассматривается задача на экране) за исключением точек края, функция $f$ одна и та же для обеих сторон поверхности.

В случае внешней задачи для замкнутой поверхности ($\Omega = {{\Omega }^{ + }}$), а также в случае разомкнутой поверхности, ставится условие

$u(x) \to 0\quad {\text{п р и }}\quad \left| x \right| \to \infty .$

В каждой точке поверхности $\Sigma $, за исключением точек края в случае замкнутой поверхности, должны существовать краевые значения у функции $u$ и у ее первых производных (со стороны области $\Omega $ в случае замкнутой поверхности $\Sigma $ и на обеих сторонах в случае разомкнутой поверхности). В случае разомкнутой поверхности ставится требование ограниченности функции $u$ в окрестности края.

Будем искать функцию $u$ в виде потенциала двойного слоя:

(5)

$u(x) = \int\limits_\Sigma {g(y)\frac{{\partial F(x - y)}}{{\partial {{n}_{y}}}}} d{{\sigma }_{y}} \equiv U[\Sigma ,g](x),$

где $g$ – неизвестная плотность.

Пусть ${{H}^{\mu }}(\Sigma ),\;\mu \in (0,1]$ есть пространство функций, непрерывных по Гёльдеру на поверхности $\Sigma $ с показателем $\mu $, ${{H}^{{1,\mu }}}(\Sigma ),\mu \in (0,1]$ – пространство функций, непрерывно диффернцируемых на поверхности $\Sigma $ и у которых поверхностный градиент $Gradg$ непрерывен по Гёльдеру с показателем $\mu \in (0,1]$ на поверхности $\Sigma $ (понятие поверхностного градиента – см. [13, c. 45]).

Если выполнено условие $g \in {{H}^{{1,\mu }}}(\Sigma )$, то в каждой точке $x \in \Sigma $, не лежащей на краю в случае разомкнутой поверхности, существуют краевые значения производной функции $u$, для которых справедливо равенство (см. [6, c. 125–129]):

(6)

$\mathop {\left( {\frac{{\partial u}}{{\partial n}}} \right)}\nolimits^ \pm (x) = \int\limits_\Sigma {g(y)\frac{{{{\partial }^{2}}F(x - y)}}{{\partial {{n}_{x}}\partial {{n}_{y}}}}} d{{\sigma }_{y}},\quad x \in \Sigma ,$

где интеграл в правой части определен в смысле конечного значения по Адамару (2).

Таким образом, краевая задача Неймана для уравнения Лапласа может быть сведена к гиперсингулярному интегральному уравнению (1).

Необходимость рассмотрения уравнения (1) с сильной особенностью возникла при решении краевой задачи Неймана (4) в случае, когда поверхность $\Sigma $ является тонким экраном. Такая задача не может быть сведена к интегральному уравнению Фредгольма II рода с интегрируемым ядром, как это делается в случае классической задачи Неймана в области, ограниченной замкнутой поверхностью. Первоначально основным инструментом теоретического анализа уравнений вида (1) была теория псевдо-дифференциальных операторов, действующих в пространствах Соболева. Так, в работе [14] доказана разрешимость уравнения (1) в некотором классе обобщенных функций и предложена схема типа Галеркина для его приближенного решения.

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

В случае замкнутой поверхности ищется решение уравнения (1) $g \in {{H}^{{1,\mu }}}(\Sigma )$. В статье [8] показано, что при $f \in {{H}^{{1,\mu }}}(\Sigma )$, $\mu \in (0,1)$, следующее условие:

(7)

$\int\limits_\Sigma {f(y)d{{\sigma }_{y}}} = 0,$

является необходимым и достаточным условием разрешимости уравнения (1). При этом общее решение имеет вид $g(x) = {{g}_{0}}(x) + A,$ где $A$ – произвольная константа, ${{g}_{0}} \in {{H}^{{1,\mu }}}$, причем функцию ${{g}_{0}}$ можно выбрать удовлетворяющей условию (7).

В случае разомкнутой поверхности ищется решение в классе функций $g \in C(\Sigma )$, удовлетворяющих условию (3) и таких, что на любой поверхности $\Sigma {\kern 1pt} ' \subset \Sigma $ – части поверхности $\Sigma $, отделенной от ее края положительным расстоянием, выполнено условие $g \in {{H}^{{1,\mu }}}(\Sigma {\kern 1pt} ')$ с некоторым $\mu \in (0,1]$ (легко показать, что для функций $g$ такого вида формула (6) по-прежнему верна во всех точках $x$).

В случае, когда поверхность $\Sigma $ есть плоский экран, являющийся выпуклым множеством на плоскости с кусочно-гладкой границей $\partial \Sigma $, существование и единственность такого решения для случая $f \in {{H}^{\mu }}(\Sigma )$ следуют из результатов, полученных в статье [15].

3. КВАДРАТУРНАЯ ФОРМУЛА И СИСТЕМА ЛИНЕЙНЫХ УРАВНЕНИЙ

Опишем предлагаемый численный метод решения уравнения (1) и приведем основные результаты, полученные в статье [12].

Осуществим триангуляцию исходной поверхности $\Sigma $ таким образом, что полученная приближенная простая поверхность $\tilde {\Sigma }$ является объединением треугольников ${{\sigma }_{i}}$, $i = 1,\; \ldots ,\;M$, вершины которых лежат на поверхности $\Sigma $, и треугольники образуют конформную сетку на поверхности $\tilde {\Sigma }$. Последнее требование означает, что множеством пересечения любых двух различных треугольников являются либо их общая целая сторона, либо одна общая вершина, либо пустое множество. В случае замкнутой исходной поверхности $\Sigma $ приближенная поверхность $\tilde {\Sigma }$ также должна быть замкнута. В случае разомкнутой поверхности $\Sigma $ приближенная поверхность $\tilde {\Sigma }$ также есть разомкнутая поверхность с краем $\partial{ \tilde {\Sigma }}$, который является замкнутой ломаной линией.

Обозначим $X = \{ {{x}_{j}} \in \Sigma ,j = 1,\; \ldots ,\;N\} $ множество всех различных вершин треугольников ${{\sigma }_{i}}$ в случае замкнутой поверхности и множество различных вершин треугольников ${{\sigma }_{i}}$, не лежащих на краю $\partial{ \tilde {\Sigma }}$ в случае разомкнутой поверхности. Будем называть такие вершины внутренними.

Сформируем для каждой внутренней вершины ${{x}_{j}} \in X$ множество треугольных ячеек $\sigma _{k}^{j}$ таких, что ${{x}_{j}} \in \sigma _{k}^{j}$, $k = 1,\; \ldots ,\;{{K}_{j}}$, и пусть ${{\tilde {\Sigma }}_{j}} = \cup _{{k = 1}}^{{{{K}_{j}}}}\sigma _{k}^{j}$ (${{\tilde {\Sigma }}_{j}}$ есть поверхность, объединяющая все ячейки, для которых точка ${{x}_{j}}$ является одной из вершин).

Построим на приближенной поверхности $\tilde {\Sigma }$ систему базисных функций ${{\varphi }_{j}}$, $j = 1,\; \ldots ,\;N$, пирамидального вида, каждая из которых является линейной функций на каждом из треугольников ${{\sigma }_{m}}$, $m = 1,\; \ldots ,\;M$, и удовлетворяет условиям:

(8)

${{\varphi }_{j}}({{x}_{i}}) = \left\{ \begin{gathered} 1,\quad i = j, \hfill \\ 0,\quad i \ne j,\quad i = 1,\;...,\;N. \hfill \\ \end{gathered} \right.$

На каждом треугольнике $\sigma _{k}^{j} \in {{\tilde {\Sigma }}_{j}}$, $k = 1,\; \ldots ,\;{{K}_{j}}$, функция ${{\varphi }_{j}}(x)$ принимает значение 1 в точке ${{x}_{j}}$ и значение $0$ в двух других вершинах. На ячейках, не лежащих на поверхности ${{\tilde {\Sigma }}_{j}}$, выполнено равенство ${{\varphi }_{j}}(x) \equiv 0$.

Пусть $g$ – искомое решение уравнения (1). Рассмотрим на поверхности $\tilde {\Sigma }$ линейную на каждой из ячеек аппроксимации ${{\sigma }_{i}}$, $i = 1,\;...,\;M$, функцию

(9)

$\tilde {g}(x) = \sum\limits_{j = 1}^N \,g({{x}_{j}}){{\varphi }_{j}}(x).$

Функция $\tilde {g}(x)$ является непрерывной на поверхности $\tilde {\Sigma }$ и удовлетворяет условию $\tilde {g}({{x}_{i}}) = g({{x}_{i}})$, $i = 1,\;...,\;N$. В статье [12] предложена квадратурная формула для значений интеграла $(Ig)$ в точках ${{x}_{i}} \in X$, основанная на замене поверхности интегрирования $\Sigma $ на поверхность $\tilde {\Sigma }$ и функции $g(x)$ на функцию $\tilde {g}(x)$ в подынтегральном выражении.

Поскольку функция $\tilde {g}(x)$ не дифференцируема в точках ${{x}_{i}}$, простой такой замены недостаточно. Для каждой точки ${{x}_{i}} \in X$ при построении квадратурной формулы для значения интеграла в этой точке строится ее трехмерная окрестность ${{U}_{\varepsilon }}({{x}_{i}})$ радиуса $\varepsilon > 0$, на поверхности $\tilde {\Sigma }$ выделяется участок ${{\tilde {\Sigma }}_{\varepsilon }}({{x}_{i}}) = \tilde {\Sigma } \cap {{U}_{\varepsilon }}({{x}_{i}})$, на котором функция $\tilde {g}$ всюду полагается равной $g({{x}_{i}})$. Введем обозначение ${{\tilde {\Sigma }}_{j}}({{x}_{i}},\varepsilon ) = {{\tilde {\Sigma }}_{j}}{\backslash }{{U}_{\varepsilon }}({{x}_{i}})$ для множества ячеек, содержащих вершину ${{x}_{j}}$ и получающихся из множества ячеек ${{\tilde {\Sigma }}_{j}}$ в результате вырезания окрестности ${{U}_{\varepsilon }}({{x}_{i}})$.

Предположение 1. Параметр $\varepsilon > 0$ выбирается так, что ${{U}_{\varepsilon }}({{x}_{i}}) \cap {{\sigma }_{j}} = {\not {0}}$ при ${{x}_{i}} \notin {{\sigma }_{j}}$.

Обозначим:

(10)

${{J}_{\varepsilon }}({{x}_{i}}) = \mathop {\left. {\left( {\frac{\partial }{{\partial {{n}_{x}}}}\int\limits_{{{{\tilde {\Sigma }}}_{\varepsilon }}({{x}_{i}})} {\frac{{\partial F(x - y)}}{{\partial {{n}_{y}}}}d{{\sigma }_{y}}} } \right)} \right|}\nolimits_{x = {{x}_{i}}} .$

Отметим, что градиент интеграла, стоящего в круглых скобках, имеет на поверхности $\tilde {\Sigma }$ в окрестности рассматриваемой точки ${{x}_{i}}$ устранимую особенность [12], и в данной формуле понимается значение производной от записанного интеграла при $x \to {{x}_{i}}$.

В результате получим следующую квадратурную формулу:

(11)

$\tilde {I}({{x}_{i}}) = \sum\limits_{j = 1}^N \,g({{x}_{j}})\int\limits_{{{{\tilde {\Sigma }}}_{j}}({{x}_{i}},\varepsilon )} {{{\varphi }_{j}}(y)} \mathop {\left. {\frac{{{{\partial }^{2}}F(x - y)}}{{\partial {{n}_{x}}\partial {{n}_{y}}}}d{{\sigma }_{y}}{\kern 1pt} } \right|}\nolimits_{x = {{x}_{i}}} + g({{x}_{i}}){{J}_{\varepsilon }}({{x}_{i}}).$

Записывая интегральное уравнение (1) в узлах коллокации ${{x}_{i}} \in X$ и заменяя интегральный оператор его приближением по формуле (11), получим систему линейных уравнений относительно неизвестных ${{g}_{i}}$, аппроксимирующих значения функции $g$ в точках ${{x}_{i}}$, $i = 1,\;...,\;N$:

(12)

$\begin{gathered} \sum\limits_{j = 1}^N \,{{a}_{{ij}}}g({{x}_{j}}) = f({{x}_{i}}),\quad i = 1, \ldots ,N, \\ {{a}_{{ij}}} = a_{{ij}}^{0} + \delta _{i}^{j}{{J}_{\varepsilon }}({{x}_{i}}),a_{{ij}}^{0} = \int\limits_{\mathop {\tilde {\Sigma }}\nolimits_j ({{x}_{i}},\varepsilon )} {{{\varphi }_{j}}(y)} \mathop {\left. {\frac{{{{\partial }^{2}}F(x - y)}}{{\partial {{n}_{x}}\partial {{n}_{y}}}}d{{\sigma }_{y}}{\kern 1pt} } \right|}\nolimits_{x = {{x}_{i}}} ,\quad \delta _{i}^{i} = 1,\quad \delta _{i}^{j} = 0,\quad i \ne j. \\ \end{gathered} $

В случае разомкнутой поверхности будем решать систему (12).

В случае замкнутой поверхности данная система является вырожденной (сумма элементов каждой строки матрицы ее системы равна $0$ [12]). Это отражает и свойства исходного уравнения (1): для него есть условие разрешимости на правую часть (7) и решение не единственно. Поэтому решение, лежащее в классе функций, удовлетворяющих условию (7), ищем из системы линейных уравнений (см. [12]):

(13)

$\sum\limits_{j = 1}^N \,{{a}_{{ij}}}{{g}_{j}} + \lambda = {{f}_{i}},\quad \sum\limits_{j = 1}^N \,{{g}_{j}}{{s}_{j}} = 0,$

относительно неизвестных ${{g}_{1}},\; \ldots ,\;{{g}_{N}},\lambda $, где $\lambda $ есть регуляризирующая переменная,

${{s}_{j}} = \int\limits_{\tilde {\Sigma }} {{{\varphi }_{j}}(y)d{{\sigma }_{y}}} .$

В статье [12] доказано, что если решение уравнения (1) удовлетворяет условию $g \in {{H}^{{1,\mu }}}(\Sigma )$, $\mu \in (0,1)$, то при определенных предположениях о свойствах замкнутой поверхности $\Sigma $ и при условии, что все углы треугольников разбиения не могут быть меньше некотрого ${{\alpha }_{0}} > 0$, найдется константа $C = C(g)$, не зависящая от разбиения такая, что:

${\text{|}}g({{x}_{i}}) - {{g}_{i}}{\text{|}} \leqslant C\Phi (h),\quad i = 1,\;...,\;N,\quad {\text{|}}\lambda {\text{|}} \leqslant C\Phi (h),\quad \Phi (h) = {{h}^{\mu }}({\text{|}}ln\varepsilon {\text{|}} + 1).$

4. ВЫЧИСЛЕНИЕ КОЭФФИЦИЕНТОВ

Перейдем к рассмотрению вопроса о вычислении коэффициентов ${{a}_{{ij}}}$ системы (12), для которых получим аналитические выражения. Выражение для коэффициентов $a_{{ij}}^{0}$ перепишем в виде:

(14)

$a_{{ij}}^{0} = \sum\limits_{k = 1}^{{{K}_{j}}} \,a_{{ij}}^{{0k}},\quad a_{{ij}}^{{0k}} = \int\limits_{\sigma _{{k,\varepsilon }}^{j}({{x}_{i}})} {{{\varphi }_{j}}(y)} \mathop {\left. {\frac{{{{\partial }^{2}}F(x - y)}}{{\partial {{n}_{x}}\partial {{n}_{y}}}}d{{\sigma }_{y}}{\kern 1pt} } \right|}\nolimits_{x = {{x}_{i}}} ,\quad \sigma _{{k,\varepsilon }}^{j}({{x}_{i}}) = \sigma _{k}^{j}{\backslash }{{U}_{\varepsilon }}({{x}_{i}}).$

Обозначим также, $S_{{k,\varepsilon }}^{j}({{x}_{i}}) = \sigma _{k}^{j} \cap {{U}_{\varepsilon }}({{x}_{i}})$ – сектор с центром в точке ${{x}_{i}}$, вырезаемый на ячейке $\sigma _{k}^{j}$ $\varepsilon $-окрестностью этой точки, $l_{{k,\varepsilon }}^{j}({{x}_{i}}) = \sigma _{k}^{j} \cap \partial {{U}_{o}}n({{x}_{i}})$ – дугу сектора $S_{{k,\varepsilon }}^{j}({{x}_{i}})$, $k = 1,\;...,\;{{K}_{j}}$ – см. фиг. 1. При этом в случае, когда точка ${{x}_{i}}$ не является вершиной ячейки $\sigma _{k}^{j}$, множества точек $S_{{k,\varepsilon }}^{j}({{x}_{i}})$ и $l_{{k,\varepsilon }}^{j}({{x}_{i}})$ являются пустыми. На фиг. 1 серым цветом выделено множество ${{\tilde {\Sigma }}_{j}}({{x}_{i}},\varepsilon )$ – носитель функции ${{\varphi }_{j}}$ за вычетом точек $\varepsilon $-окрестности точки ${{x}_{i}}$.

Фиг. 1.

Структура области интегрирования при вычислении коэффициентов ${{a}_{{ij}}}$: (a) точки ${{x}_{i}}$ и ${{x}_{j}}$ не являются вершинами одного треугольника, (б) существует треугольник с вершинами в точках ${{x}_{i}}$ и ${{x}_{j}}$, $i \ne j$, (в) $i = j.$

Таким образом, необходимо вычислить интегралы, входящие в формулу (14), каждый из которых берется по поверхности $\hat {\sigma } = \sigma _{{k,\varepsilon }}^{j}({{x}_{i}})$, которая есть плоская поверхность, не содержащая точку ${{x}_{i}}$, одного из следующих видов (см. фиг. 1):

– $\hat {\sigma }$ есть треугольник $\sigma _{k}^{j}$, (в случае, когда точка ${{x}_{i}}$ не является вершиной треугольника $\sigma _{k}^{j}$) – ячейки на фиг. 1a, а также некоторые из ячеек на фиг. 1б;

– $\hat {\sigma }$ есть поверхность $\sigma _{{k,\varepsilon }}^{j}({{x}_{i}})$ – треугольник $\sigma _{k}^{j}$ с вершиной в точке ${{x}_{i}}$, из которого вырезан сектор $S_{{k,\varepsilon }}^{j}({{x}_{i}})$ радиуса $\varepsilon $ с центром в точке ${{x}_{i}}$ – случай (б), если при этом $i \ne j$, или случай (в), если при этом $i = j$.

Кроме того, при вычислении коэффициента ${{a}_{{ij}}}$ при $i = j$ (случай (в)) вычисляется интеграл (10).

Для нахождения интегралов по поверхности $\hat {\sigma }$ будем использовать формулу

$\int\limits_{\hat {\sigma }} {{{\varphi }_{j}}(y)\frac{{{{\partial }^{2}}F({{x}_{i}} - y)}}{{\partial {{n}_{{{{x}_{i}}}}}\partial {{n}_{y}}}}} d{{\sigma }_{y}} = {\mathbf{W}}[\hat {\sigma },{{\varphi }_{j}}](x){\mathbf{n}}({{x}_{i}}),$

где мы ввели обозначение

${\mathbf{W}}[S,f](x) = \operatorname{grad} U[S,f](x) \equiv \operatorname{grad} \int\limits_S {f(y)\frac{{\partial F(x - y)}}{{\partial {{n}_{y}}}}} d{{\sigma }_{y}},$

$U[S,f]$ – потенциал двойного слоя с плостностью $f$ вида (5), размещенный на поверхности $S$.

Сначала проведем некоторые вспомогательные рассуждения.

Пусть $S$ – простая гладкая поверхность класса ${{C}^{3}}$ с краем $\partial S$, $f \in {{H}^{{1,\mu }}}(S)$. Тогда для точек $x$, не лежащих на поверхности $S$, справедлива формула (см. [13, с. 68–69]):

(15)

$\begin{gathered} {\mathbf{W}}[S,f](x) = \int\limits_S {\left[ {[Gradf(y) \times {{{\mathbf{n}}}_{y}}] \times gra{{d}_{x}}F(x - y)} \right]} {\kern 1pt} d{{\sigma }_{y}} + \\ + \;\int\limits_{\partial S} {f(y)({\mathbf{\tau }}(y)} \times gra{{d}_{x}}F(x - y))d{{s}_{y}},\quad x \in {{R}^{3}}{\backslash }\Sigma , \\ \end{gathered} $

где $d{{s}_{y}}$ – дифференциал дуги, ${\mathbf{\tau }}(y)$ – орт вектора касательной к контуру $\partial S$, выбираемый так, чтобы при движении вдоль контура в направлении вектора ${\mathbf{\tau }}(y)$ при нормали к поверхности, направленной вверх, поверхность оставалась слева. При этом градиент функции $u$ имеет в каждой точке $x \in S{\backslash }\partial S$ краевые значения, для которых справедлива формула

${{(gradu)}^{ \pm }} = gradu \pm Gradf{\text{/}}2,$

где через $gradu$ обозначено прямое значение, получаемое при непосредственной подстановке точки $x$ в правую часть формулы (15), первый интеграл понимается в смысле главного значения. Доказательство формулы (15) для случая замкнутой поверхности приведено в [13], и в этом случае в ней отсутствует последний интеграл по контуру $\partial S$. Для разомкнутой поверхности с краем формула доказывается аналогичным образом, используя небольшое уточнение в самом конце рассуждения из [13, с. 69].

Пусть теперь $\hat {\sigma }$ – ограниченная плоская поверхность с краем, $\varphi $ – линейная на поверхности $\hat {\sigma }$ функция, точка $x$ лежит вне поверхности $\hat {\sigma }$. Используя формулу (15), представим векторное поле ${\mathbf{W}}[\hat {\sigma },\varphi ](x)$ в виде суммы интегралов:

${\mathbf{W}}[\hat {\sigma },\varphi ](x) = \int\limits_{\hat {\sigma }} {\left[ {[Grad\varphi \times {\mathbf{n}}] \times gra{{d}_{x}}F(x - y)} \right]} {\kern 1pt} d{{\sigma }_{y}} + \int\limits_{\partial{ \hat {\sigma }}} {\varphi (y){\mathbf{\tau }}(y)} \times gra{{d}_{x}}F(x - y)d{{s}_{y}}.$

Последнюю формулу запишем в виде

(16)

${\mathbf{W}}[\hat {\sigma },\varphi ](x) = \gamma \times \mathcal{I}[\hat {\sigma }](x) + \mathcal{J}[\varphi ,\partial{ \hat {\sigma }}](x),$

(17)

$\begin{gathered} \mathcal{I}[\hat {\sigma }](x) = \int\limits_{\hat {\sigma }} {{\mathbf{V}}(x - y)} d{{\sigma }_{y}},\quad \mathcal{J}[\varphi ,L](x) = \int\limits_L {\varphi (y){\mathbf{\tau }}(y)} \times {\mathbf{V}}(x - y)d{{s}_{y}}, \\ {\mathbf{\gamma }} = Grad\varphi \times {\mathbf{n}},\quad {\mathbf{V}}(x - y) = gra{{d}_{x}}F(x - y) = - \frac{{x - y}}{{4\pi {\text{|}}x - y{{{\text{|}}}^{3}}}}, \\ \end{gathered} $

где ${\mathbf{n}}$ – постоянный вектор нормали к плоской ячейки $\hat {\sigma }$, $Grad\varphi $ – поверхностный градиент функции $\varphi $, причем вектор $Grad\varphi $ является постоянным на ячейке $\hat {\sigma }$ в силу линейности функции $\varphi $. В определении интеграла $\mathcal{J}[\varphi ,L]$ кривая $L$ есть кусочно гладкая кривая, замкнутая или разомкнутая.

Рассмотрим интеграл $\mathcal{I}[\hat {\sigma }]$.

Если функция $f$ определена и непрерывно дифференцируема в трехмерной окрестности поверхности (включая саму поверхность), то на поверхности определен поверхностный градиент, который связан с обычным (пространственным) градиентом формулой:

(18)

$Gradf = gradf - {\mathbf{n}}\frac{{\partial f}}{{\partial n}}.$

Пусть $S$ – простая гладкая поверхность с краем $\partial S$ класса ${{C}^{2}}$, функция $\varphi $ является непрерывно дифференцируемой на поверхности $S$ (включая край). Тогда справедлива формула (см. [13, с. 46]):

(19)

$\int\limits_S {Grad\varphi d\sigma } = \int\limits_{\partial S} {\varphi {\mathbf{\nu }}ds} - 2\int\limits_S {\varphi H{\mathbf{n}}d\sigma } ,$

здесь ${\mathbf{\nu }}$ – единичный вектор нормали к граничной кривой $\partial S$, который перпендикулярен нормали ${\mathbf{n}}$ к поверхности и направлен наружу от поверхности $S$, $H$ – средняя кривизна поверхности.

Так как рассматриваемая поверхность $\hat {\sigma }$ плоская, на ней функция $H$ – средняя кривизна, равна нулю, а ${\mathbf{n}}(y) = {\mathbf{n}}$, где через ${\mathbf{n}}$ мы обозначили нормаль к плоскости, в которой лежит ячейка. Используя равенство $gra{{d}_{x}}F(x - y) = - {{\operatorname{grad} }_{y}}F(x - y)$ и формулы (18), (19), получаем

$\begin{gathered} \mathcal{I}[\hat {\sigma }](x) = - \frac{1}{{4\pi }}\int\limits_{\hat {\sigma }} {gra{{d}_{y}}} \frac{1}{{{\text{|}}x - y{\text{|}}}}d{{\sigma }_{y}} = \\ = \; - {\kern 1pt} \frac{1}{{4\pi }}\int\limits_{\partial \mathop {\hat {\sigma }}\nolimits_y } {\frac{{{\mathbf{\nu }}(y)}}{{{\text{|}}x - y{\text{|}}}}} ds - \frac{1}{{4\pi }}\int\limits_{\hat {\sigma }} {{\mathbf{n}}\frac{\partial }{{\partial {{n}_{y}}}}\frac{1}{{{\text{|}}x - y{\text{|}}}}} d{{\sigma }_{y}} = - {{\mathcal{I}}_{1}}[\hat {\sigma }](x) - {{\mathcal{I}}_{2}}[\hat {\sigma }](x). \\ \end{gathered} $

Для интеграла ${{\mathcal{I}}_{2}}$, с учетом того, что ячейка $\hat {\sigma }$ плоская, сразу можем записать:

(20)

${{\mathcal{I}}_{2}}[\hat {\sigma }](x) = {\mathbf{n}}\int\limits_{\hat {\sigma }} {\frac{{\partial F(x - y)}}{{\partial {{n}_{y}}}}} d{{\sigma }_{y}} = \frac{1}{{4\pi }}{\mathbf{n}}\Theta (\hat {\sigma },x),$

где $\Theta (\hat {\sigma },x)$ – ориентированный телесный угол, под которым видна ячейка $\hat {\sigma }$ из точки $x$ (этот угол считается положительным, если из точки $x$ видна положительная сторона поверхности, и отрицательным в противоположном случае). Заметим, что когда поверхность представляет $\hat {\sigma }$ собой треугольник с вырезанным сектором вокруг точки $x$, являющейся вершиной треугольника, телесный угол равен $0$. Алгоритм расчета телесного угла, под которым виден из заданной точки треугольник описан, например, в [16].

Рассмотрим интеграл ${{\mathcal{I}}_{1}}$. Этот интеграл нам необходимо вычислить для двух видов поверхности $\hat {\sigma }$: треугольник и плоская поверхность, возникающая при вырезании из треугольника сектора вокруг одной из вершин, причем радиус этого сектора меньше, чем расстояние от его центра до других вершин. Этот интеграл разбивается на интегралы двух видов.

Пусть $AB$ – дуга окружности с центром ${{x}_{i}}$ и радиусом $\varepsilon $, ${{{\mathbf{\nu }}}_{{AB}}}(y)$ – направленный к точке ${{x}_{i}}$ единичный вектор нормали к дуге $AB$ в лежащей на ней точке $y$. Рассмотрим интеграл

${{{\mathbf{I}}}_{\varepsilon }}[{{x}_{i}},AB] \equiv \int\limits_{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{A} B} {\frac{{{{{\mathbf{\nu }}}_{{AB}}}}}{{{\text{|}}{{x}_{i}} - y{\text{|}}}}} ds = \frac{1}{{{{\varepsilon }^{2}}}}\int\limits_{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{A} B} {({{x}_{i}} - y)ds} .$

Введем систему координат с началом координат в точке ${{x}_{i}}$ и базисными векторами ${{{\mathbf{e}}}_{1}} = \tfrac{{{{x}_{i}} - A}}{{{\text{|}}{{x}_{i}} - A{\text{|}}}}$, ${{{\mathbf{e}}}_{2}} = {\mathbf{n}} \times {{{\mathbf{e}}}_{1}}$. Здесь ${\mathbf{n}}$ – это нормаль к сектору $A{{x}_{i}}B$. Тогда в этой системе координат получаем

(21)

${{{\mathbf{I}}}_{\varepsilon }}[{{x}_{i}},AB] = \frac{1}{{{{\varepsilon }^{2}}}}\int\limits_A^B {( - {{\xi }_{1}}, - {{\xi }_{2}})ds = \frac{1}{{{{\varepsilon }^{2}}}}} \int\limits_0^\psi {( - \varepsilon cos\alpha , - \varepsilon sin\alpha )\varepsilon d\alpha } = - sin\psi {{{\mathbf{e}}}_{1}} + (cos\psi - 1){{{\mathbf{e}}}_{2}},$

где $\psi $ – угол раствора сектора. Заметим, что данный интеграл не зависит от радиуса окружности.

Также, для вычисления интеграла ${{\mathcal{I}}_{1}}$, рассмотрим следующий интеграл по отрезку $AB$:

(22)

$\begin{gathered} {{k}_{{AB}}}(x) \equiv \int\limits_A^B {\frac{{d{{s}_{y}}}}{{{\text{|}}x - y{\text{|}}}}} = \{ y = A + t(B - A)\} = {\text{|}}B - A{\kern 1pt} {\text{|}}\,\int\limits_0^1 {\frac{{dt}}{{{\text{|}}x - A - t(B - A){\text{|}}}}} = \\ = \;\left\{ {\alpha = {\text{|}}B - A{\text{|}}{\kern 1pt} ,\;\beta = \frac{{(B - A,x - A)}}{\alpha },\;\gamma = - {{\beta }^{2}}\; + \;{\text{|}}x - A{\kern 1pt} {{{\text{|}}}^{2}},\;z = \alpha t - \beta } \right\} = \int\limits_{ - \beta }^{\alpha - \beta } {\frac{{dz}}{{\sqrt {{{z}^{2}} + {{\gamma }^{2}}} }}} , \\ \end{gathered} $

(23)

${{k}_{{AB}}}(x) = ln\left| {\frac{{(B - A,B - x)\; + \;{\text{|}}x - B{\text{||}}B - A{\text{|}}}}{{(B - A,A - x)\; + \;{\text{|}}x - A{\text{||}}B - A{\text{|}}}}} \right|.$

Заметим, что полученное выражение, является перестановочным по отношению к точкам $A$ и $B$. Перестановочность важна, так как по построению триангуляции точка $x$ обязательно лежит вне отрезка $AB$, однако может попасть на его продолжение, что может привести к вырождению выражения (23). Если точка $x$ принадлежит прямой, проходящей через точки $A$ и $B$ и, кроме того, $(B - A,A - x) < 0$, то выражение (23) вырождается, однако, в силу перестановочности можно использовать выражение

(24)

${{k}_{{AB}}}(x) = ln\left| {\frac{{(A - B,A - x)\; + \;{\text{|}}x - A{\text{||}}A - B{\text{|}}}}{{(A - B,B - x)\; + \;{\text{|}}x - B{\text{||}}A - B{\text{|}}}}} \right|.$

Заметим, что случай одновременного вырождения обоих выражений в рассматриваемой численной схеме не возникает, так как в вычисляемых в ней интегралах точка $x = {{x}_{i}}$ не может находиться на отрезке $AB$ – стороне ячейки $\hat {\sigma }$.

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

(25)

${\text{Е с л и }}\quad \left( {B - A,\frac{{B + A}}{2} - x} \right) \geqslant 0,\quad {\text{т о }}\;{\text{и с п о л ь з у е т с я }}\;{\text{ф о р м у л а }}\;\quad\left( {23} \right),\;{\text{и н а ч е }}\;--\;\left( {24} \right).$

Теперь мы можем записать следующие выражения для интеграла ${{\mathcal{I}}_{1}}$:

(i) если поверхность $\hat {\sigma }$ есть треугольник $ABC$ ($\hat {\sigma } = \sigma _{k}^{j}$ в случае ${{x}_{i}} \notin \sigma _{k}^{j}$, см. случаи (а) и (б) на фиг. 1), то

(26)

${{\mathcal{I}}_{1}}[\hat {\sigma }]({{x}_{i}}) = \frac{1}{{4\pi }}\left\{ {{{{\mathbf{\nu }}}_{{AB}}}{{k}_{{AB}}}({{x}_{i}}) + {{{\mathbf{\nu }}}_{{BC}}}{{k}_{{BC}}}({{x}_{i}}) + {{{\mathbf{\nu }}}_{{CA}}}{{k}_{{CA}}}({{x}_{i}})} \right\};$

(ii) если поверхность $\hat {\sigma }$ есть треугольник $ABC$ с вырезанным сектором ($\hat {\sigma } = \sigma _{{k,\varepsilon }}^{j}({{x}_{i}})$ в случаях (б) или (в), см. фиг. 1), причем вершина ${{x}_{i}}$, вокруг которой вырезается сектор, соответствует вершине $A$, а дуга $l_{{k,\varepsilon }}^{j}({{x}_{i}})$ имеет концы $D \in [AC]$ и $E \in [AB]$, то

(27)

${{\mathcal{I}}_{1}}[\hat {\sigma }]({{x}_{i}}) = \frac{1}{{4\pi }}\left\{ {{{{\mathbf{\nu }}}_{{EB}}}{{k}_{{EB}}}({{x}_{i}}) + {{{\mathbf{\nu }}}_{{BC}}}{{k}_{{BC}}}({{x}_{i}}) + {{{\mathbf{\nu }}}_{{CD}}}{{k}_{{CD}}}({{x}_{i}}) + {{{\mathbf{I}}}_{\varepsilon }}[{{x}_{i}},DE]} \right\},$

где для вычисления ${{k}_{{AB}}}({{x}_{i}})$ используются формулы (23)–(25), интеграл ${{{\mathbf{I}}}_{\varepsilon }}$ вычисляется в соответствии с формулой (21), ${{{\mathbf{\nu }}}_{{AB}}}$ – вектор нормали к стороне $AB$, лежащий в плоскости рассматриваемой ячейки и направленный наружу от ячейки.

Теперь рассмотрим вектор ${\mathbf{\gamma }}$ из формулы (16). Для нахождения этого вектора необходимо вычислить вектор $Grad\varphi $ – поверхностный градиент линейной функции $Grad\varphi (y)$, который является постоянным вектором на плоской ячейке.

Рассмотрим базисную функцию $\varphi (y) = {{\varphi }_{j}}(y)$ на треугольнике $ABC = \sigma _{k}^{j}$, который является одним из треугольников разбиения, на котором эта базисная функия не равна нулю тождественно. Считаем, что вершины треугольника обозначены так, что $\varphi (A) = 1,\;\varphi (B) = 0,\;\varphi (C) = 0$. Пусть ${\mathbf{c}} = {\mathbf{AB}}$, ${\mathbf{b}} = {\mathbf{AC}}$. В силу линейности функции $\varphi (x)$, справедливо равенство: $\varphi (x) = \varphi (A) + (x - A) \cdot Grad\varphi $, $x \in ABC$. Тогда вектор $Grad\varphi $ ищем из системы уравненний:

$\begin{gathered} {\mathbf{c}} \cdot Grad\varphi = - 1, \\ {\mathbf{b}} \cdot Grad\varphi = - 1, \\ {\mathbf{n}} \cdot Grad\varphi = 0. \\ \end{gathered} $

Решением этой системы является вектор

(28)

$Grad\varphi = \alpha {\mathbf{c}} + \beta {\mathbf{b}},\quad \alpha = \tfrac{{({\mathbf{c}},{\mathbf{b}}) - ({\mathbf{b}},{\mathbf{b}})}}{\Delta },\quad \beta = \tfrac{{({\mathbf{c}},{\mathbf{b}}) - ({\mathbf{c}},{\mathbf{c}})}}{\Delta },\quad \Delta = ({\mathbf{c}},{\mathbf{c}})({\mathbf{b}},{\mathbf{b}}) - {{({\mathbf{c}},{\mathbf{b}})}^{2}}.$

Перейдем к рассмотрению контурного интеграла $\mathcal{J}$ в формуле (16).

Напомним, что формула (16) применяется к каждой ячейке разбиения. При вычислении коэффициента ${{a}_{{ij}}}$ интеграл $\mathcal{J}$ вычисляется для функции ${{\varphi }_{j}}$ по краю $\partial \sigma _{{k,\varepsilon }}^{j}({{x}_{i}})$ каждой ячейки $\sigma _{{k,\varepsilon }}^{j}({{x}_{i}})$, $k = 1,\;...,\;{{K}_{j}}$, на которой функция ${{\varphi }_{j}}$ не равна нулю тождественно, причем по сторонам ячеек, выходящим из точки ${{x}_{j}}$, мы посчитаем этот интеграл дважды для двух соседних ячеек в разных направлениях. На сторонах, не содержащих точку ${{x}_{j}}$, функция ${{\varphi }_{j}}$ обращается в нуль. Поэтому справедлива формула:

$\sum\limits_{k = 1}^{{{K}_{j}}} \,\mathcal{J}[{{\varphi }_{j}},\partial \sigma _{{k,\varepsilon }}^{j}({{x}_{i}})]({{x}_{i}}) = \sum\limits_{k = 1}^{{{K}_{j}}} \,\mathcal{J}[{{\varphi }_{j}},l_{{k,\varepsilon }}^{j}({{x}_{i}})]({{x}_{i}}),$

означающая, что вклад в коэффициент ${{a}_{{ij}}}$ вносят только интегралы по дугам $l_{{k,\varepsilon }}^{j}({{x}_{i}})$. Рассмотрим вычисление таких интегралов.

Пусть треугольник $ABC$ есть ячейка триангуляции, содержащая одновременно узлы ${{x}_{i}}$ и ${{x}_{j}}$ (эти узлы могут как совпадать, так и быть соседними вершинами), причем $A = {{x}_{i}}$ – вершина, вокруг которой вырезается сектор $AB{\kern 1pt} {\text{'}}C{\kern 1pt} {\text{'}}$, где точки $B{\kern 1pt} {\text{'}} \in [AB]$ и $C{\kern 1pt} {\text{'}} \in [AC]$ есть концы дуги окружности, ограничивающей этот сектор. Обозначим $\alpha ,\beta ,\gamma $ углы у вершин $A,B$ и $C$ соответственно.

Рассмотрим интеграл $\mathcal{J}[\varphi ,L]$ по дуге $L = C{\kern 1pt} {\text{'}}B{\kern 1pt} {\text{'}}$, где $\varphi $ – линейная функция, принимающая значение 1 в одной из вершин треугольника и значение 0 в двух других.

В случае, когда $\varphi (A) = 1$ (случай $i = j$, фиг. 2a), из теоремы синусов для треугольника $ABM$ получаем

$\varphi (y) = 1 - \frac{\varepsilon }{{{\text{|}}AM{\text{|}}}} = 1 - sin(\psi + \beta )\frac{\varepsilon }{{{\text{|}}AB{\text{|}}sin\beta }},$

где $y \in L$, $\psi $ – угол между лучами $[AB)$ и $[Ay)$, $M = [Ay) \cap [B,C]$. Тогда получаем формулу

(29)

$\mathcal{J}[\varphi ,L]({{x}_{i}}) = - \frac{1}{{4\pi {{\varepsilon }^{2}}}}{\mathbf{n}}\int\limits_{C'B'} {\varphi (y)ds} = \frac{{\mathbf{n}}}{{4\pi \varepsilon }}\left( {\alpha - \frac{\varepsilon }{{ABsin\beta }}(cos\gamma + cos\beta )} \right),$

где ${\mathbf{n}}$ – нормаль к плоскости треугольника.

Фиг. 2.

(a) $i = j$, $\varphi (A) = 1$, (б) $i \ne j$, $\varphi (A) = 0.$

Eсли $\varphi (A) = 0$ (случай $i \ne j$, фиг. 2б), $\varphi (C) = 1$, аналогичным образом имеем

(30)

$\mathcal{J}[\varphi ,L]({{x}_{i}}) = - {\mathbf{n}}\left( {\frac{1}{{4\pi ACsin\alpha }}(1 - cos\alpha )} \right).$

Перейдем к рассмотрению интеграла ${{J}_{\varepsilon }}({{x}_{i}})$, определяемого формулой (10). Применяя формулу (15) c $f = 1$ к каждому из элементов $S_{{i,\varepsilon }}^{k}({{x}_{i}})$, $k = 1,\;...,\;{{K}_{i}}$ и учитывая, что интегралы по сторонам ячеек $S_{{k,\varepsilon }}^{i}({{x}_{i}})$, которые лежат на лучах, выходящих из вершины ${{x}_{i}}$, взаимно сокращаются для соседних ячеек, получаем (см. фиг. 1)

${{J}_{\varepsilon }}({{x}_{i}}) = - \sum\limits_{k = 1}^{{{K}_{i}}} \int\limits_{l_{{k,\varepsilon }}^{i}({{x}_{i}})} {\frac{{{\mathbf{n}}({{x}_{i}}) \cdot ({\mathbf{\tau }}(y) \times ({{x}_{i}} - y))}}{{4\pi |{{x}_{i}} - y{{|}^{3}}}}} ds = - \frac{1}{{4\pi {{\varepsilon }^{2}}}}\sum\limits_{k = 1}^{{{K}_{i}}} \int\limits_{l_{{k,\varepsilon }}^{i}({{x}_{i}})} {ds} .$

Из последнего выражения имеем

(31)

${{J}_{\varepsilon }}({{x}_{i}}) = - \frac{1}{{4\pi \varepsilon }}\sum\limits_{k = 1}^{{{K}_{i}}} \,{{\psi }^{i}},$

где ${{\psi }^{i}}$ – угол раствора сектора $S_{{k,\varepsilon }}^{i}$.

Таким образом, мы получили следующие формулы для вычисления рассматриваемых коэффициентов матрицы:

(32)

$\begin{gathered} {{a}_{{ij}}} = \sum\limits_{k = 1}^{{{K}_{j}}} \,\tilde {a}_{{ij}}^{{0k}} + \delta _{i}^{j}{{J}_{\varepsilon }}({{x}_{i}}),\quad \delta _{i}^{i} = 1,\quad \delta _{i}^{j} = 0,\quad i \ne j, \\ \tilde {a}_{{ij}}^{{0k}} = {\mathbf{n}}({{x}_{i}}) \cdot \left( {Grad{{\varphi }_{j}} \times {\mathbf{n}}(\sigma _{k}^{j}) \times [{{\mathcal{I}}_{1}}[\hat {\sigma }]({{x}_{i}}) + {{\mathcal{I}}_{2}}[\hat {\sigma }]({{x}_{i}})] + \mathcal{J}[{{\varphi }_{j}},l_{{k,\varepsilon }}^{j}({{x}_{i}})]} \right),\quad t\sigma = \sigma _{{k,\varepsilon }}^{j}({{x}_{i}}), \\ \end{gathered} $

где ${\mathbf{n}}({{x}_{i}})$ – вектор нормали к поверхности $\Sigma $ в точке ${{x}_{i}}$, вектор $Grad{{\varphi }_{j}}$ – поверхностный градиент функции $\varphi $ на ячейке $\sigma _{k}^{j}$, являющийся постоянным вектором, рассчитывается по формуле (28), ${\mathbf{n}}(\sigma _{k}^{j})$ – нормаль к ячейке $\sigma _{k}^{j}$. Интеграл ${{\mathcal{I}}_{1}}$ вычисляется по формулам (26), (27) в зависимости от вида ячейки $\hat {\sigma }$, интеграл ${{\mathcal{I}}_{2}}$ вычисляется по формуле (20). Слагаемое $\mathcal{J}$ рассчитывается по формуле (29), если $i = j$, и по формуле (30) если $i \ne j$, причем этот интеграл отсутствует (полагается равным 0), если $\sigma _{{k,\varepsilon }}^{j}({{x}_{i}}) = \sigma _{k}^{j}$. Для вычисления интеграла ${{J}_{\varepsilon }}({{x}_{i}})$ используется формула (31).

Замечание. При практических расчетах, приведенных далее, вектор нормали ${\mathbf{n}}({{x}_{i}})$ в вершине ${{x}_{i}}$ рассчитывался как сумма нормалей всех треугольников, которым принадлежит точка ${{x}_{i}}$, осредненная по площадям этих треугольников.

5. ЧИСЛЕННЫЙ ЭКСПЕРИМЕНТ: ПРИМЕРЫ РАСЧЕТОВ И ОБСУЖДЕНИЕ

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

5.1. Уравнение на сфере

Рассмотрим задачу о потенциальном обтекании сферы радиуса $R = 1$ потоком идеальной несжимаемой жидкости, имеющим на бесконечности постоянную скорость ${{{\mathbf{w}}}_{\infty }}$ такую, что ${\text{|}}{{{\mathbf{w}}}_{\infty }}{\text{|}} = 1$. Пусть $\Sigma $ – поверхность сферы. Задача состоит в отыскании поля скоростей жидкости ${\mathbf{w}} = {{{\mathbf{w}}}_{\infty }} + gradu$, в котором потенциал $u$ является решением краевой задачи (4) в области ${{\Omega }^{ + }}$ вне сферы с $f = - ({\mathbf{n}},{{{\mathbf{w}}}_{\infty }})$ (такое условие равносильно требованию ${\mathbf{wn}} = 0$ на поверхности сферы). Как упоминалось выше, данная краевая задача сводится к интегральному уравнению (1), если искать потенциал $u$ в виде потенциала двойного слоя (5).

Для сформулированной задачи известно аналитическое выражение для поля скоростей ${\mathbf{w}}$, которое позволяет выписать аналитическое решение возникающего на поверхности сферы $\Sigma $ уравнения (1). В частности (см. [17, c. 42]), на поверхности сферы выражение для вектора скорости ${\mathbf{w}}$ имеет вид

(33)

${\mathbf{w}}(x) = w{\mathbf{\tau }}(x),\quad w \equiv w(\theta ) = \frac{3}{2}sin\theta ,\quad {\mathbf{\tau }}(x) = {{{\mathbf{r}}}_{x}} \times {\mathbf{N}},$

где $x$ – точка на поверхности сферы, ${{{\mathbf{r}}}_{x}}$ – радиус-вектор точки $x$ относительно центра сферы, $\theta = \theta (x)$ – угол между векторами ${{{\mathbf{r}}}_{x}}$ и вектором $ - {{{\mathbf{w}}}_{\infty }}$, ${\mathbf{N}}$ – единичный вектор нормали к плоскости, образованной векторами ${{{\mathbf{r}}}_{x}}$ и ${{{\mathbf{w}}}_{\infty }}$ (при этом ${\mathbf{\tau }}$ есть единичнй вектор касательной к дуге большой окружности, проведенной через точку $x$ параллельно вектору ${{{\mathbf{w}}}_{\infty }}$, если угол $\theta $ равен нулю, вектор ${\mathbf{N}}$ не определен, в этом случае ${\mathbf{w}}(x) = 0$ по непрерывности).

Чтобы отсюда найти функцию $g$ – решение уравнения (1), заметим, что если функция $g$ есть решение этого уравнения, то потенциал двойного слоя (5) является решением одновременно и внутренней, и внешней задач (4) в силу формулы (6). Тогда функция $U(x) = {{{\mathbf{w}}}_{\infty }}{{{\mathbf{r}}}_{x}} + u(x)$, которая в области вне сферы есть потенциал полного течения (${\mathbf{w}} = \operatorname{grad} U$), определена данным выражением и внутри сферы. При этом в области ${{\Omega }^{ - }}$ внутри сферы потенциал $U(x)$ удовлетворяет уравнению Лапласа, а на ее границе $\Sigma $ условию ${{(\partial U(x){\text{/}}\partial n)}^{ - }} = 0$. Значит, $U(x) = {\text{const}}$ в области ${{\Omega }^{ - }}$. Поскольку краевые значения потенциала двойного слоя $u$ на поверхности $\Sigma $ удовлетворяют условию ${{u}^{ + }} - {{u}^{ - }} = g$, заключаем, что на поверхности сферы выполнено условие $g - {{U}^{ + }} = {\text{const}}$.

Теперь возьмем точку $P$ – полюс сферы, обращенный к набегающему потоку, которой соответствует полярный угол $\theta (P) = 0$. Пусть $x \in \Sigma $ – произвольная точка на поверхности сферы, $L(x)$ – дуга большой окружности с началом в точке $P$ и концом в точке $x$. Тогда

$g(x) - g(P) = \int\limits_{L(x)} {{\mathbf{w\tau }}ds} = \int\limits_0^\theta {w(\theta {\kern 1pt} ')d\theta {\kern 1pt} '} = \frac{3}{2}(1 - cos\theta ),$

где $\theta = \theta (x)$, ${\mathbf{\tau }} = {\mathbf{\tau }}(y)$, $y \in L(x)$ и $w = w(\theta )$ – те же, что и в формуле (33). При этом решение, удовлетворяющее условию (7), имеет вид

$g(x) = - \frac{3}{2}cos\theta .$

Была проведена серия расчетов, в которых решалось уравнение (1) с применением триангуляции поверхности на различное количество ячеек (см. фиг. 3). Использовались разбиения, где в качестве начальной поверхности брался правильный октаэдр с вершинами на координатных осях декартовой системы координат $O{{x}_{1}}{{x}_{2}}{{x}_{3}}$ с началом в центре сферы. Далее последовательно строились новые разбиения, получаемые добавлением вершин на лучах, выходящих из центра сферы в центр каждого ребра предыдущего разбиения. На фиг. 3 для некоторых из разбиений, используемых в расчетах, изображено сечение плоскостью ${{x}_{2}} = {{x}_{3}}$, на котором осуществлялось сравнение полученных численно значений функции $g$ и ее поверхностного градиента $Gradg$ с точными значениями. При этом вектор скорости набегающего потока ${{{\mathbf{w}}}_{\infty }}$ направлен вдоль оси $O{{x}_{1}}$.

Фиг. 3.

Разбиения поверхности сферы.

На графиках на фиг. 4 приведены зависимости значения искомой функции $g$ от сферической координаты $\theta $ вдоль линии, лежащей в плоскости ${{x}_{2}} = {{x}_{3}}$ (cм. фиг. 3 ), для различных разбиений (на рисунке указано количество уравнений в системе (13)). На фиг. 5 приведены аналогичные зависимости для величины ${\text{|}}Gradg{\text{|}}$, а также показано поле скоростей жидкости, полученное численно в одном из расчетов. С физической точки зрения ${\text{|}}Gradg{\text{|}} = w$, где $w$ – значение модуля касательной скорости на поверхности сферы (см. формулу (33)). Значения функции $g$ вычислялись в центрах ячеек, попавших на сечение, по формуле (9), в которой за значения $g({{x}_{i}})$, $i = 1,\;...,\;N$, принимались значения ${{g}_{i}}$, полученные при решении системы (13). Поверхностный градиент функции $g$ вычислялся по формуле:

$Grad\tilde {g}(x) = \sum\limits_{j = 1}^N \,{{g}_{j}}Grad{{\varphi }_{j}}(x),$

следующей из формулы (9), градиент базисной функции $Grad{{\varphi }_{j}}(x)$, являющийся постоянным вектором на соответствующей ячейке, вычислялся по формуле (28).

Фиг. 4.

Метод кусочно-линейных аппроксимаций: значение функции $g$ на сечении.

Фиг. 5.

Метод кусочно-линейных аппроксимаций: значение $\left| {Gradg} \right|$ на сечении.

Наблюдается практическая сходимость метода, как в отношении вычисления функции $g$, так и в отношении ее поверхностного градиента.

Следует отметить, что важное значение имеет выбор параметра $\varepsilon $ – радиуса регуляризирующей окрестности. Тестирование было проведено для различного выбора радиуса окрестности области регуляризации, в том числе “динамического”, т.е. подбираемого для каждой вершины пропорционально минимальной высоте, выходящей из этой вершины на стороны содержащих ее треугольников. Однако наилучшая сходимость наблюдалась при фиксированном значении параметра $\varepsilon $, общем для всех точек и близком к максимально возможному в соответствии с предположением 1: в приведенных расчетах мы полагали $\varepsilon = 0.99{{h}_{{{\text{min}}}}}$, где ${{h}_{{{\text{min}}}}}$ – минимальная из всех высот всех треугольников разбиения.

Для сравнения аналогичные расчеты были проведены с применением метода кусочно-постоянных аппроксимаций и коллокаций (метод вихревых рамок), описанного в [2], с. 473–475. При сравнении с предлагаемым в настоящей статье методом вычислялись зачения функции $g$ и ее поверхностного градиента $Gradg$ в центрах ячеек, лежащих на сечении (фиг. 3), при использовании тех же разбиений поверхности сферы (за значение в центре ячейки принималось значение кусочно-постоянной функции $g$ на этой ячейке).

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

$Gradg = w(\theta ){\mathbf{\tau }},\quad w(\theta ) \approx \frac{1}{2}\left( {\frac{{g(\theta ) - g(\theta {\kern 1pt} ')}}{{\theta - \theta {\kern 1pt} '}} + \frac{{g(\theta {\kern 1pt} '') - g(\theta )}}{{\theta {\kern 1pt} {\text{''}} - \theta }}} \right),$

где $\theta $ – значение полярного угла в рассматриваемой точке, попавшей на сечение, $\theta {\kern 1pt} '$ и $\theta {\kern 1pt} ''$ – значения полярного угла в соседних точках на этом сечении (такой способ естественен для точек имеенно выбранных сечений), ${\mathbf{\tau }}$ – тот же, что и в формуле (33), $g(\theta )$ – значение функции $g$, соответствующее точке $x$ с угловым положением $\theta $ на данном сечении.

На фиг. 6 приведены зависимости максимальной абсолютной погрешности для значений функции $g$ и вектора $Gradg$ в центрах ячеек рассматриваемых сечений, в зависимости от числа уравнений при использовании кусочно-линейного и кусочно-постоянного методов аппроксимации. Здесь следует заметить, что при одном и том же разбиении число возникающих дискретных уравнений для сравниваемых методов различно: для замкнутого тела оно на 1 больше числа вершин в случае кусочно-линейного метода и числа ячеек в случае кусочно-постоянного. По оси абсцисс отложено число уравнений, поскольку именно оно, а не число ячеек, определяет вычислительные затраты. Видно, что точность метода кусочно-линейных аппроксимаций выше, особенно при вычислении поверхностного градиента.

Фиг. 6.

Максимальная абсолютная ошибка на сечении: (a) функции $g$, (б) вектора $Gradg.$

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

Фиг. 7.

Средняя ошибка (a) функции $g$, (б) вектора $Gradg.$

5.2. Уравнение на пластине

В качестве следующего теста была рассмотрена задача о потенциальном обтекании прямоугольной пластины потоком идеальной несжимаемой жидкости, имеющим на бесконечности постоянную скорость ${{{\mathbf{w}}}_{\infty }}$, направленную перпендикулярно поверхности пластины. Была взята пластина в форме прямоугольника в плоскости $O{{x}_{2}}{{x}_{3}}$, определяемого соотношениями ${{x}_{2}} \in [ - 0.5,0.5]$, ${{x}_{3}} \in [ - 1,1]$, вектор скорости набегающего потока брался в виде ${{{\mathbf{w}}}_{\infty }} = (1,0,0)$. Так же, как и задача об обтекании сферы, эта задача сводится к уравнению (1) на поверхности пластины $\Sigma $ с правой частью $f = - {\mathbf{(n}},{{{\mathbf{w}}}_{\infty }})$. При этом поле скоростей опять имеет вид ${{w}_{\infty }} + gradu$, где $u$ есть потенциал двойного слоя (5). Отметим, что здесь решается интгеральное уравнение (1) на разомкнутой поверхности.

На фиг. 8 приведены зависимости неизвестной функции $g$ от координаты ${{x}_{3}}$ вдоль линии на пластине, заданной уравнением ${{x}_{2}} = 0$, полученные для различных разбиений пластины – пластина сначала разбивалась на квадратные ячейки, затем каждая квадратная ячейка разбивалась диагоналями на 4 треугольника. На рисунке в формате $n \times m$ указано число квадратов вдоль сторон пластины. Здесь также видна тенденция к сеточной сходимости, причем наибольшие различия между кривыми наблюдаются в точке максимума. Также показано поле скоростей жидкости, полученное численно в одном из расчетов, где интегральное уравнение решалось методом кусочно-линейных аппроксимаций.

Фиг. 8.

Значение функции $g$ на отрезке – сечении ${{x}_{2}} = 0$, метод кусочно-линейных аппроксимаций.

На графике фиг. 9 приведена зависимость максимального значения функции $g$ от числа уравнений. Для сравнения, приведена аналогичная зависимость, полученная методом кусочно-постоянных аппроксимаций на тех же разбиениях поверхности пластины (при одном и том же разбиении число уравнений для методов кусочно-линейных и кусочно-постоянных аппроксимаций различно).

Фиг. 9.

Максимальное значение функции g.

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

В статье разработана численная схема решения гиперсингулярного интегрального уравнения (1) на основе метода кусочно-линейных аппроксимаций и коллокаций. Уравнение сведено к системе линейных алгебраических уравнений, для коэффициентов которой получены аналитические выражения. Произведено тестирование разработанного метода на примере решения интегральных уравнений, возникающих при решении задач о потенциальном обтекании идеальной несжимаемой жидкостью сферы и прямоугольной пластины. Сравнение с известным методом кусочно-постоянных аппроксимаций (методом вихревых рамок), широко используемым в вычислительной аэродинамике, показало более высокую скорость сходимости разработанного метода при вычислении неизвестной функции и, особенно, при вычислении ее поверхностного градиента. Последнее важно с точки приложений. В частности, применяемые в аэродинамике панельные методы, основаны на решении уравнений такого типа. При этом поверхностный градиент плотности потенциала двойного слоя, относительно которой решается уравнение, необходим для нахождения касательных скоростей на поверхности тела [16].

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

Апаринов В.А., Белоцерковский С.М., Лифанов И.К., Михайлов А.А. Расчет нестационарных аэродинамических характеристик тел при отрывном обтекании // Ж. вычисл. матем. и матем. физ. 1988. Т. 2. № 16. С. 1558–1566.
Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент. М.: ТОО Янус, 1995. 520 с.
Гутников В.А., Кирякин В.Ю., Лифанов И.К., Сетуха А.В., Ставцев С.Л. О численном решении двумерного гиперсингулярного интегрального уравнения и о распространении звука в городской застройке // Ж. вычисл. матем. и матем. физ. 2007. Т. 47. № 12. С. 2088–2100.
Daeva S.G., Setukha A.V. Numerical simulation of scattering of acoustic waves by inelastic bodies using hypersingular boundary integral equation // AIP Conference Proceedings. V. 1648. AIP Publishing, 2015. P. 390004-1–390004-4.
Lifanov I.K., Stavtsev S.L., Piven V.F. Mathematical Modelling of the Three-Dimensional Boundary Value Problem of the Discharge of the Well System in a Homogeneous Layer // Russ. J. Numer. Anal. Math. Modelling. 2002. V. 17. № 1. P. 99–111.
Вайникко Г.М., Лифанов И.К., Полтавский Л.Н. Численные методы в гиперсингулярных интегральных уравнениях и их приложения. М.: Янус-К, 2001. 508 с.
Лебедева С.Г., Сетуха А.В. О численном решении полного двумерного гиперсингулярного уравнения методом дискретных особенностей // Дифференц. ур-ния. 2013. Т. 49. № 2. С. 223–233.
Рыжаков С.Г., Сетуха А.В. О сходимости численного метода решения некоторого гиперсингулярного интегрального уравнения на замкнутой поверхности // Дифференц. ур-ния. 2010. Т. 46. № 9. С. 1343–1353.
Рыжаков С.Г., Сетуха А.В. О сходимости метода вихревых рамок с регуляризацией для краевой задачи Неймана на плоском экране // Дифференц. ур-ния. 2011. Т. 47. № 9. С. 1352–1358.
Рыжаков С.Г., Сетуха А.В. Ссатуф И. О численном решении краевой задачи Неймана на экране методом вихревых рамок с использованием неравномерной сетки // Ученые записки Орловского гос. ун-та. 2010. № 4. С. 12–19.
Сетуха А.В. Сходимость численного метода решения гиперсингулярного интегрального уравнения на отрезке с применением кусочно-линейных аппроксимаций на неравномерной сетке // Дифференц. ур-ния. 2017. Т. 53. № 2. С. 237–249.
Сетуха А.В., Семенова А.В. Сходимость метода кусочно-линейных аппроксимаций и коллокаций для некоторого гиперсингулярного интегрального уравнения на замкнутой поверхности // Дифференц. ур-ния. 2017. Т. 53. № 9. С. 1265–1280.
Колтон Д., Кресс Р. Методы интегральных уравнений в теории рассеяния: Пер. с англ. М.: Мир, 1987. 311 с.
Stephan E.P. Boundary Integral Equations for Screen Problem in ${{R}^{3}}$ // Integral Equations and Operator Theory. 1987. V. 10. № 2. P. 236–257.
Сетуха А.В. Трехмерная краевая задача Неймана с обобщенными граничными условиями и уравнение Прандтля // Дифференц. ур-ния. 2003. Т. 39. № 9. С. 1188–1200.
Гутников В.А., Лифанов И.К., Сетуха А.В. О моделировании зданий и сооружений методом дискретных вихревых рамок // Изв. РАН. Механ. жидкости и газа. 2006. № 4. С. 78–92.
Ландау Л.Д., Лифшиц Е.М. Теоретическая физика: Учебное пособие. Т. 6. Гидродинамика. 3-е изд. М.: Наука, 1986. 736 с.

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

Инструменты

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

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

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