Журнал вычислительной математики и математической физики, 2021, T. 61, № 4, стр. 608-624

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

Р. Ф. Марданов^1, *, А. Е. Марданова^1, **

¹ К(П)ФУ
420008 Казань, ул. Кремлевская, 18, Россия

^* E-mail: Renat.Mardanov@kpfu.ru
^** E-mail: AlEMardanova@kpfu.ru

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

DOI: 10.31857/S0044466921040086

Аннотация

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

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

ВВЕДЕНИЕ

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

(0.1)

${{\Delta }^{2}}f = v,$

где $\Delta = \frac{{{{\partial }^{2}}}}{{\partial {{x}^{2}}}} + \frac{{{{\partial }^{2}}}}{{\partial {{y}^{2}}}}$ – оператор Лапласа. Уравнение (0.1) широко применяется в механике сплошных сред: например, в задачах расчета изгиба плоской упругой пластины под действием нагрузки [1] и течения вязкой жидкости при малых числах Рейнольдса в рамках модели Стокса [2]. В приведенных примерах правая часть уравнения (0.1) является известной функцией, и численные методы решения таких задач хорошо развиты [3], [4].

Отдельный класс задач механики сводится к более сложному бигармоническому уравнению, правая часть которого содержит саму искомую функцию и ее производные. Классическим примером является задача о расчете вязкого течения в рамках модели Навье–Стокса с использованием переменных функции тока $\psi $ и завихренности $\omega = - \Delta \psi $. В этом случае функция тока $\psi $ удовлетворяет уравнению [5], [6]

(0.2)

${{\Delta }^{2}}\psi = {\text{Re}}\left( {\frac{{\partial \psi }}{{\partial y}}\frac{{\partial \omega }}{{\partial x}} - \frac{{\partial \psi }}{{\partial x}}\frac{{\partial \omega }}{{\partial y}}} \right),$

где ${\text{Re}}$ – число Рейнольдса. Уравнение (0.2) нелинейное и его решение находится в основном численными итерационными методами [3], [7], [8].

Настоящая работа посвящена решению неоднородного бигармонического уравнения с линейной правой частью, содержащей искомую функцию и ее производные. Предложенный подход основан на методе граничных и областных элементов (МГОЭ) (boundary domain integral method (BDIM)) [5], [9]–[11], позволяющем решить задачу без организации итерационного процесса. Проведена апробация метода для тестовой задачи сравнением с известным аналитическим решением. В качестве примера практического применения выполнено решение задачи о расчете плоского фильтрационного течения в пористой среде с пространственно неоднородным распределением проницаемости. Полученное решение сопоставлено с результатом вычислительного эксперимента по детальному расчету течения вязкой жидкости в межпоровом пространстве модельной пористой среды в рамках модели Стокса. Результаты числовых расчетов для обеих задач подтвердили эффективность и хорошую точность предложенного численного метода.

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

Рассмотрим функцию $f(x,y)$, удовлетворяющую в области $\Omega $, ограниченной контуром $\Gamma $, неоднородному бигармоническому уравнению (0.1). Правая часть $v(x,y)$ в общем случае зависит от искомой функции $f$ и ее производных:

(1.1)

$v = {{v}_{0}} + {{v}_{1}}f + {{v}_{2}}\frac{{\partial f}}{{\partial x}} + {{v}_{3}}\frac{{\partial f}}{{\partial y}} + {{v}_{4}}\Delta f + {{v}_{5}}\frac{{\partial (\Delta f)}}{{\partial x}} + {{v}_{6}}\frac{{\partial (\Delta f)}}{{\partial y}},$

где заданные функции ${{v}_{i}} = {{v}_{i}}(x,y)$, $i = \overline {0,6} $, являются функциями координат $(x,y)$. Граничные условия могут иметь следующий вид:

(1.2)

${{\left. f \right|}_{\Gamma }} = {{h}_{1}}(s),\quad {{\left. {f{\kern 1pt} '{\kern 1pt} } \right|}_{\Gamma }} = {{h}_{2}}(s),\quad {{\left. {\Delta f} \right|}_{\Gamma }} = {{h}_{3}}(s),\quad {{\left. {(\Delta f){\kern 1pt} '{\kern 1pt} } \right|}_{\Gamma }} = {{h}_{4}}(s),$

где ${{h}_{i}}(s)$, $i = \overline {1,4} $ – известные функции, $s$ – дуговая абсцисса контура $\Gamma $, отсчитываемая от некоторой точки так, что область $\Omega $ остается слева, а штрих означает производную по направлению внешней нормали $n$, т.е. $(\,){\kern 1pt} ' = \partial (\,){\text{/}}\partial n$. Граница $\Gamma $ разбита на участки ${{\Gamma }^{{(j)}}}$, $j = \overline {1,N} $, на каждом из которых заданы по два граничных условия вида (1.2), в зависимости от конкретной решаемой задачи.

Требуется определить функцию $f(x,y).$

2. РЕШЕНИЕ

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

(2.1)

$\Delta f = g,\quad \Delta g = v$

для функций $f(x,y)$ и $g(x,y)$. Используя известное интегральное соотношение Рэлея–Грина для бигармонического уравнения [6], запишем эквивалентную пару интегральных уравнений:

(2.2)

$\begin{array}{*{20}{c}} {\chi (x,y)f(x,y) = \int\limits_\Gamma \,[f(s)G_{1}^{'}(x,y,s) - f{\kern 1pt} '(s){{G}_{1}}(x,y,s) + } \\ { + \;g(s)G_{2}^{'}(x,y,s) - g{\kern 1pt} '(s){{G}_{2}}(x,y,s)]ds + \int\limits_\Omega \,v(\xi ,\eta ){{G}_{2}}(x,y,\xi ,\eta )d\xi d\eta ,} \end{array}$

(2.3)

$\chi (x,y)g(x,y) = \int\limits_\Gamma \,[g(s)G_{1}^{'}(x,y,s) - g{\kern 1pt} '(s){{G}_{1}}(x,y,s)]ds + \int\limits_\Omega \,{v}(\xi ,\eta ){{G}_{1}}(x,y,\xi ,\eta )d\xi d\eta ,$

где $\chi (x,y) = 2\pi $ для внутренних точек $(x,y) \in \Omega $, $\chi (x,y) = \beta $ для граничных точек $(x,y) \in \Gamma $ ($\beta $ – внутренний к области $\Omega $ угол в точке на границе $\Gamma $). Функции Грина для бигармонического уравнения имеют вид

(2.4)

$\begin{gathered} {{G}_{1}} = ln\rho ,\quad {{G}_{2}} = \frac{{{{\rho }^{2}}}}{4}(ln\rho - 1), \\ \rho (x,y,s) = \sqrt {{{{({{x}_{1}}(s) - x)}}^{2}} + {{{({{y}_{1}}(s) - y)}}^{2}}} ,\quad \rho (x,y,\xi ,\eta ) = \sqrt {{{{(\xi - x)}}^{2}} + {{{(\eta - y)}}^{2}}} , \\ \end{gathered} $

где $({{x}_{1}},{{y}_{1}})$ – координаты точки интегрирования на границе $\Gamma $ с дуговой абсциссой $s$, а $(\xi ,\eta )$ – координаты точки интегрирования в области $\Omega $.

Функция $v(x,y)$ (1.1), которую с учетом (2.1) перепишем в виде

(2.5)

$v = {{v}_{0}} + {{v}_{1}}f + {{v}_{2}}\frac{{\partial f}}{{\partial x}} + {{v}_{3}}\frac{{\partial f}}{{\partial y}} + {{v}_{4}}g + {{v}_{5}}\frac{{\partial g}}{{\partial x}} + {{v}_{6}}\frac{{\partial g}}{{\partial y}},$

содержит производные неизвестных функций $f$ и $g$. Для получения замкнутой системы интегральных уравнений продифференцируем (2.2) и (2.3) по переменным $x$ и $y$:

(2.6)

$\begin{gathered} \chi (x,y)\frac{{\partial f(x,y)}}{{\partial x}}\int\limits_\Gamma \,\left[ {f(s)\frac{{\partial G_{1}^{'}(x,y,s)}}{{\partial x}} - f{\kern 1pt} '(s)\frac{{\partial {{G}_{1}}(x,y,s)}}{{\partial x}}} \right. + \\ \left. {\mathop + \limits_{_{{_{{_{{_{{_{{_{{}}}}}}}}}}}}} g(s)\frac{{\partial {{G}_{2}}(x,y,s)}}{{\partial x}} - g{\kern 1pt} '(s)\frac{{\partial {{G}_{2}}(x,y,s)}}{{\partial x}}} \right]ds + \int\limits_\Omega \,v(\xi ,\eta )\frac{{\partial {{G}_{2}}(x,y,\xi ,\eta )}}{{\partial x}}d\xi d\eta , \\ \end{gathered} $

(2.7)

$\begin{gathered} \chi (x,y)\frac{{\partial f(x,y)}}{{\partial y}}\int\limits_\Gamma \,\left[ {f(s)\frac{{\partial G_{1}^{'}(x,y,s)}}{{\partial y}} - f{\kern 1pt} '(s)\frac{{\partial {{G}_{1}}(x,y,s)}}{{\partial y}}} \right. + \\ \left. {\mathop + \limits_{_{{_{{_{{_{{_{{_{{_{{}}}}}}}}}}}}}}} g(s)\frac{{\partial G_{2}^{'}(x,y,s)}}{{\partial y}} - g{\kern 1pt} '(s)\frac{{\partial {{G}_{2}}(x,y,s)}}{{\partial y}}} \right]ds + \int\limits_\Omega \,v(\xi ,\eta )\frac{{\partial {{G}_{2}}(x,y,\xi ,\eta )}}{{\partial y}}d\xi d\eta , \\ \end{gathered} $

(2.8)

$\chi (x,y)\frac{{\partial g(x,y)}}{{\partial x}} = \int\limits_\Gamma \,\left[ {g(s)\frac{{\partial G_{1}^{'}(x,y,s)}}{{\partial x}} - g{\kern 1pt} '(s)\frac{{\partial {{G}_{1}}(x,y,s)}}{{\partial x}}} \right]ds + \int\limits_\Omega \,v(\xi ,\eta )\frac{{\partial {{G}_{1}}(x,y,\xi ,\eta )}}{{\partial x}}d\xi d\eta ,$

(2.9)

$\chi (x,y)\frac{{\partial g(x,y)}}{{\partial y}} = \int\limits_\Gamma \,\left[ {g(s)\frac{{\partial G_{1}^{'}(x,y,s)}}{{\partial y}} - g{\kern 1pt} '(s)\frac{{\partial {{G}_{1}}(x,y,s)}}{{\partial y}}} \right]ds + \int\limits_\Omega \,{v}(\xi ,\eta )\frac{{\partial {{G}_{1}}(x,y,\xi ,\eta )}}{{\partial y}}d\xi d\eta .$

Согласно методу граничных элементов, представим границу $\Gamma = \bigcup\nolimits_{j = 1}^n \,{{\Gamma }_{j}}$ в виде набора прямолинейных отрезков ${{\Gamma }_{j}}$, а область $\Omega = \bigcup\nolimits_{k = 1}^m \,{{\Omega }_{k}}$ в виде набора площадных элементов ${{\Omega }_{k}}$ (трех- или четырехугольных). Функции $f(s)$, $f{\kern 1pt} '(s)$, $g(s)$, $g{\kern 1pt} '(s)$ аппроксимируем кусочно-постоянными функциями со значениями ${{f}_{j}}$, $f_{j}^{'}$, ${{g}_{j}}$, $g_{j}^{'}$ на отрезках ${{\Gamma }_{j}}$, а функции ${{v}_{i}}(\xi ,\eta )$, $i = \overline {0,6} $, $f(\xi ,\eta )$, $\partial f{\text{/}}\partial x(\xi ,\eta )$, $\partial f{\text{/}}\partial y(\xi ,\eta )$, $g(\xi ,\eta )$, $\partial g{\text{/}}\partial x(\xi ,\eta )$, $\partial g{\text{/}}\partial y(\xi ,\eta )$ – кусочно-постоянными функциями со значениями ${{v}_{{ik}}}$, $\mathop {\tilde {f}}\nolimits_k $, ${{f}_{{xk}}}$, ${{f}_{{yk}}}$, $\mathop {\tilde {g}}\nolimits_k $, ${{g}_{{xk}}}$, ${{g}_{{yk}}}$ в площадных элементах ${{\Omega }_{k}}$. Тогда уравнения (2.2), (2.3), (2.6)–(2.9) с учетом (2.5) примут вид

(2.10)

$\chi (x,y)f(x,y) = \sum\limits_{j = 1}^n \,\left\{ {{{f}_{j}}\int\limits_{{{\Gamma }_{j}}} \,G_{1}^{'}ds - f_{j}^{'}\int\limits_{{{\Gamma }_{j}}} \,{{G}_{1}}ds + \left. {{{g}_{j}}\int\limits_{{{\Gamma }_{j}}} \,G_{2}^{'}ds - g_{j}^{'}\int\limits_{{{\Gamma }_{j}}} \,{{G}_{2}}ds} \right\} + \sum\limits_{k = 1}^m \,{{{\tilde {v}}}_{k}}\int\limits_{{{\Omega }_{k}}} \,{{G}_{2}}d\xi d\eta ,} \right.$

(2.11)

$\begin{array}{*{20}{c}} {\chi (x,y)g(x,y) = \sum\limits_{j = 1}^n \,\left\{ {{{g}_{j}}\int\limits_{{{\Gamma }_{j}}} \,G_{1}^{'}ds - g_{j}^{'}\int\limits_{{{\Gamma }_{j}}} \,{{G}_{1}}ds} \right\} + \sum\limits_{k = 1}^m \,{{{\tilde {v}}}_{k}}\int\limits_{{{\Omega }_{k}}} \,{{G}_{1}}d\xi d\eta ,} \end{array}$

(2.12)

$\chi (x,y)\frac{{\partial f(x,y)}}{{\partial x}} = \sum\limits_{j = 1}^n \,\left\{ {{{f}_{j}}\int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial G_{1}^{'}}}{{\partial x}}ds - f_{j}^{'}\int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial {{G}_{1}}}}{{\partial x}}ds + \left. {{{g}_{j}}\int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial G_{2}^{'}}}{{\partial x}}ds - g_{j}^{'}\int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial {{G}_{2}}}}{{\partial x}}ds} \right\} + \sum\limits_{k = 1}^m \,{{{\tilde {v}}}_{k}}\int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{2}}}}{{\partial x}}d\xi d\eta ,} \right.$

(2.13)

$\chi (x,y)\frac{{\partial f(x,y)}}{{\partial y}} = \sum\limits_{j = 1}^n \,\left\{ {{{f}_{j}}\int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial G_{1}^{'}}}{{\partial y}}ds - f_{j}^{'}\int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial {{G}_{1}}}}{{\partial y}}ds + } \right.\left. {{{g}_{j}}\int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial G_{2}^{'}}}{{\partial y}}ds - g_{j}^{'}\int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial {{G}_{2}}}}{{\partial y}}ds} \right\} + \sum\limits_{k = 1}^m \,{{\tilde {v}}_{k}}\int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{2}}}}{{\partial y}}d\xi d\eta ,$

(2.14)

$\begin{array}{*{20}{c}} {\chi (x,y)\frac{{\partial g(x,y)}}{{\partial x}} = \sum\limits_{j = 1}^n \,\left\{ {{{g}_{j}}\int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial G_{1}^{'}}}{{\partial x}}ds - g_{j}^{'}\int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial {{G}_{1}}}}{{\partial x}}ds} \right\} + \sum\limits_{k = 1}^m \,{{{\tilde {v}}}_{k}}\int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{1}}}}{{\partial x}}d\xi d\eta ,} \end{array}$

(2.15)

$\begin{array}{*{20}{c}} {\chi (x,y)\frac{{\partial g(x,y)}}{{\partial y}} = \sum\limits_{j = 1}^n \,\left\{ {{{g}_{j}}\int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial G_{1}^{'}}}{{\partial y}}ds - g_{j}^{'}\int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial {{G}_{1}}}}{{\partial y}}ds + } \right\} + \sum\limits_{k = 1}^m \,{{{\tilde {v}}}_{k}}\int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{1}}}}{{\partial y}}d\xi d\eta ,} \end{array}$

где

${{\tilde {v}}_{k}} = {{v}_{{0k}}} + {{v}_{{1k}}}{{\tilde {f}}_{k}} + {{v}_{{2k}}}{{f}_{{xk}}} + {{v}_{{3k}}}{{f}_{{yk}}} + {{v}_{{4k}}}{{\tilde {g}}_{k}} + {{v}_{{5k}}}{{g}_{{xk}}} + {{v}_{{6k}}}{{g}_{{yk}}}.$

Рассмотрев выражения (2.10)–(2.11) в центрах линейных элементов ${{\Gamma }_{i}}$ с координатами $({{x}_{{ci}}},{{y}_{{ci}}})$, запишем

(2.16)

$\begin{gathered} \sum\limits_{j = 1}^n \,\{ A_{{ij}}^{{(1)}}{{f}_{j}}\, + \,B_{{ij}}^{{(1)}}f_{j}^{'}\, + \,C_{{ij}}^{{(1)}}{{g}_{j}}\, + \,D_{{ij}}^{{(1)}}g_{j}^{'}\} \, + \,\sum\limits_{k = 1}^m \,\{ E_{{ik1}}^{{(1)}}{{{\tilde {f}}}_{k}}\, + \,E_{{ik2}}^{{(1)}}{{f}_{{xk}}}\, + \,E_{{ik3}}^{{(1)}}{{f}_{{yk}}}\, + \,E_{{ik4}}^{{(1)}}{{{\tilde {g}}}_{k}}\, + \,E_{{ik5}}^{{(1)}}{{g}_{{xk}}}\, + \,E_{{ik6}}^{{(1)}}{{g}_{{yk}}}\} \, = \,{{a}_{i}}, \\ \sum\limits_{j = 1}^n \,\{ C_{{ij}}^{{(2)}}{{g}_{j}} + D_{{ij}}^{{(2)}}g_{j}^{'}\} + \sum\limits_{k = 1}^m \,\{ E_{{ik1}}^{{(2)}}{{{\tilde {f}}}_{k}} + E_{{ik2}}^{{(2)}}{{f}_{{xk}}} + E_{{ik3}}^{{(2)}}{{f}_{{yk}}} + E_{{ik4}}^{{(2)}}{{{\tilde {g}}}_{k}} + E_{{ik5}}^{{(2)}}{{g}_{{xk}}} + E_{{ik6}}^{{(2)}}{{g}_{{yk}}}\} = {{b}_{i}}, \\ \end{gathered} $

где $i = \overline {1,n} $. Рассмотрев выражения (2.10)–(2.15) в центрах площадных элементов ${{\Omega }_{l}}$ с координатами $({{\tilde {x}}_{{cl}}},{{\tilde {y}}_{{cl}}})$, запишем

$\begin{gathered} \sum\limits_{j = 1}^n \,\{ A_{{lj}}^{{(3)}}{{f}_{j}} + B_{{lj}}^{{(3)}}f_{j}^{'} + C_{{lj}}^{{(3)}}{{g}_{j}} + D_{{lj}}^{{(3)}}g_{j}^{'}\} + \\ \, + \sum\limits_{k = 1}^m \,\{ E_{{lk1}}^{{(3)}}{{{\tilde {f}}}_{k}} + E_{{lk2}}^{{(3)}}{{f}_{{xk}}} + E_{{lk3}}^{{(3)}}{{f}_{{yk}}} + E_{{lk4}}^{{(3)}}{{{\tilde {g}}}_{k}} + E_{{lk5}}^{{(3)}}{{g}_{{xk}}} + E_{{lk6}}^{{(3)}}{{g}_{{yk}}}\} = {{c}_{l}}, \\ \sum\limits_{j = 1}^n \,\{ C_{{lj}}^{{(4)}}{{g}_{j}} + D_{{lj}}^{{(4)}}g_{j}^{'}\} + \sum\limits_{k = 1}^m \,\{ E_{{lk1}}^{{(4)}}{{{\tilde {f}}}_{k}} + E_{{lk2}}^{{(4)}}{{f}_{{xk}}} + E_{{lk3}}^{{(4)}}{{f}_{{yk}}} + E_{{lk4}}^{{(4)}}{{{\tilde {g}}}_{k}} + E_{{lk5}}^{{(4)}}{{g}_{{xk}}} + E_{{lk6}}^{{(4)}}{{g}_{{yk}}}\} = {{d}_{l}}, \\ \end{gathered} $

(2.17)

$\begin{gathered} \sum\limits_{j = 1}^n \,\{ A_{{lj}}^{{(5)}}{{f}_{j}} + B_{{lj}}^{{(5)}}f_{j}^{'} + C_{{lj}}^{{(5)}}{{g}_{j}} + D_{{lj}}^{{(5)}}g_{j}^{'}\} + \\ \, + \sum\limits_{k = 1}^m \,\{ E_{{lk1}}^{{(5)}}{{{\tilde {f}}}_{k}} + E_{{lk2}}^{{(5)}}{{f}_{{xk}}} + E_{{lk3}}^{{(5)}}{{f}_{{yk}}} + E_{{lk4}}^{{(5)}}{{{\tilde {g}}}_{k}} + E_{{lk5}}^{{(5)}}{{g}_{{xk}}} + E_{{lk6}}^{{(5)}}{{g}_{{yk}}}\} = {{e}_{l}}, \\ \sum\limits_{j = 1}^n \,\{ A_{{lj}}^{{(6)}}{{f}_{j}} + B_{{lj}}^{{(6)}}f_{j}^{'} + C_{{lj}}^{{(6)}}{{g}_{j}} + D_{{lj}}^{{(6)}}g_{j}^{'}\} + \\ \end{gathered} $

$\begin{gathered} \, + \sum\limits_{k = 1}^m \,\{ E_{{lk1}}^{{(6)}}{{{\tilde {f}}}_{k}} + E_{{lk2}}^{{(6)}}{{f}_{{xk}}} + E_{{lk3}}^{{(6)}}{{f}_{{yk}}} + E_{{lk4}}^{{(6)}}{{{\tilde {g}}}_{k}} + E_{{lk5}}^{{(6)}}{{g}_{{xk}}} + E_{{lk6}}^{{(6)}}{{g}_{{yk}}}\} = {{h}_{l}}, \\ \sum\limits_{j = 1}^n \,\{ C_{{lj}}^{{(7)}}{{g}_{j}} + D_{{lj}}^{{(7)}}g_{j}^{'}\} + \sum\limits_{k = 1}^m \,\{ E_{{lk1}}^{{(7)}}{{{\tilde {f}}}_{k}} + E_{{lk2}}^{{(7)}}{{f}_{{xk}}} + E_{{lk3}}^{{(7)}}{{f}_{{yk}}} + E_{{lk4}}^{{(7)}}{{{\tilde {g}}}_{k}} + E_{{lk5}}^{{(7)}}{{g}_{{xk}}} + E_{{lk6}}^{{(7)}}{{g}_{{yk}}}\} = {{q}_{l}}, \\ \sum\limits_{j = 1}^n \,\{ C_{{lj}}^{{(8)}}{{g}_{j}} + D_{{lj}}^{{(8)}}g_{j}^{'}\} + \sum\limits_{k = 1}^m \,\{ E_{{lk1}}^{{(8)}}{{{\tilde {f}}}_{k}} + E_{{lk2}}^{{(8)}}{{f}_{{xk}}} + E_{{lk3}}^{{(8)}}{{f}_{{yk}}} + E_{{lk4}}^{{(8)}}{{{\tilde {g}}}_{k}} + E_{{lk5}}^{{(8)}}{{g}_{{xk}}} + E_{{lk6}}^{{(8)}}{{g}_{{yk}}}\} = {{t}_{l}}, \\ \end{gathered} $

где $l = \overline {1,m} $. Таким образом, соотношения (2.16), (2.17) представляют собой систему $2n + 6m$ линейных алгебраических уравнений (СЛАУ) относительно $4n + 6m$ неизвестных: ${{f}_{j}}$, $f_{j}^{'}$, ${{g}_{j}}$, $g_{j}^{'}$, $j = \overline {1,n} $; ${{\tilde {f}}_{k}}$, ${{f}_{{xk}}}$, ${{f}_{{yk}}}$, ${{\tilde {g}}_{k}}$, ${{g}_{{xk}}}$, ${{g}_{{yk}}}$, $k = \overline {1,m} $. Коэффициенты СЛАУ определены следующими выражениями:

$A_{{ij}}^{{(1)}} = \int\limits_{{{\Gamma }_{j}}} \,G_{1}^{'}({{x}_{{ci}}},{{y}_{{ci}}},s)ds - {{\beta }_{i}}{{\delta }_{{ij}}},\quad B_{{ij}}^{{(1)}} = - \int\limits_{{{\Gamma }_{j}}} \,{{G}_{1}}({{x}_{{ci}}},{{y}_{{ci}}},s)ds,$

$C_{{ij}}^{{(1)}} = \int\limits_{{{\Gamma }_{j}}} \,G_{2}^{'}({{x}_{{ci}}},{{y}_{{ci}}},s)ds,\quad D_{{ij}}^{{(1)}} = - \int\limits_{{{\Gamma }_{j}}} \,{{G}_{2}}({{x}_{{ci}}},{{y}_{{ci}}},s)ds,$

$E_{{ikp}}^{{(1)}} = {{{v}}_{{pk}}}\int\limits_{{{\Omega }_{k}}} \,{{G}_{2}}({{x}_{{ci}}},{{y}_{{ci}}},\xi ,\eta )d\xi d\eta ,\quad {{a}_{i}} = - \sum\limits_{k = 1}^m \,{{{v}}_{{0k}}}\int\limits_{{{\Omega }_{k}}} \,{{G}_{2}}({{x}_{{ci}}},{{y}_{{ci}}},\xi ,\eta )d\xi d\eta ,$

$C_{{ij}}^{{(2)}} = A_{{ij}}^{{(1)}},\quad D_{{ij}}^{{(2)}} = B_{{ij}}^{{(1)}},$

$E_{{ikp}}^{{(2)}} = {{{v}}_{{pk}}}\int\limits_{{{\Omega }_{k}}} \,{{G}_{1}}({{x}_{{ci}}},{{y}_{{ci}}},\xi ,\eta )d\xi d\eta ,\quad {{b}_{i}} = - \sum\limits_{k = 1}^m \,{{{v}}_{{0k}}}\int\limits_{{{\Omega }_{k}}} \,{{G}_{1}}({{x}_{{ci}}},{{y}_{{ci}}},\xi ,\eta )d\xi d\eta ,$

$A_{{lj}}^{{(3)}} = \int\limits_{{{\Gamma }_{j}}} \,G_{1}^{'}({{\tilde {x}}_{{cl}}},{{\tilde {y}}_{{cl}}},s)ds,\quad B_{{lj}}^{{(3)}} = - \int\limits_{{{\Gamma }_{j}}} \,{{G}_{1}}({{\tilde {x}}_{{cl}}},{{\tilde {y}}_{{cl}}},s)ds,$

$C_{{lj}}^{{(3)}} = \int\limits_{{{\Gamma }_{j}}} \,G_{2}^{'}({{\tilde {x}}_{{cl}}},{{\tilde {y}}_{{cl}}},s)ds,\quad D_{{lj}}^{{(3)}} = - \int\limits_{{{\Gamma }_{j}}} \,{{G}_{2}}({{\tilde {x}}_{{cl}}},{{\tilde {y}}_{{cl}}},s)ds,$

$E_{{lkp}}^{{(3)}} = {{{v}}_{{pk}}}\int\limits_{{{\Omega }_{k}}} \,{{G}_{2}}({{\tilde {x}}_{{cl}}},{{\tilde {y}}_{{cl}}},\xi ,\eta )d\xi d\eta - 2\pi {{\delta }_{{lk}}}{{\delta }_{{1p}}},\quad {{c}_{l}} = - \sum\limits_{k = 1}^m \,{{{v}}_{{0k}}}\int\limits_{{{\Omega }_{k}}} \,{{G}_{2}}({{\tilde {x}}_{{cl}}},{{\tilde {y}}_{{cl}}},\xi ,\eta )d\xi d\eta ,$

$C_{{lj}}^{{(4)}} = A_{{lj}}^{{(3)}},\quad D_{{lj}}^{{(4)}} = B_{{lj}}^{{(3)}},$

$E_{{lkp}}^{{(4)}} = {{{v}}_{{pk}}}\int\limits_{{{\Omega }_{k}}} \,{{G}_{1}}({{\tilde {x}}_{{cl}}},{{\tilde {y}}_{{cl}}},\xi ,\eta )d\xi d\eta - 2\pi {{\delta }_{{lk}}}{{\delta }_{{4p}}},\quad {{d}_{l}} = - \sum\limits_{k = 1}^m \,{{{v}}_{{0k}}}\int\limits_{{{\Omega }_{k}}} \,{{G}_{1}}({{\tilde {x}}_{{cl}}},{{\tilde {y}}_{{cl}}},\xi ,\eta )d\xi d\eta ,$

$A_{{lj}}^{{(5)}} = \int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial G_{1}^{'}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},s)}}{{\partial x}}ds,\quad B_{{lj}}^{{(5)}} = - \int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial {{G}_{1}}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},s)}}{{\partial x}}ds,$

$C_{{lj}}^{{(5)}} = \int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial G_{2}^{'}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},s)}}{{\partial x}}ds,\quad D_{{lj}}^{{(5)}} = - \int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial {{G}_{2}}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},s)}}{{\partial x}}ds,$

$E_{{lkp}}^{{(5)}} = {{{v}}_{{pk}}}\int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{2}}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},\xi ,\eta )}}{{\partial x}}d\xi d\eta - 2\pi {{\delta }_{{lk}}}{{\delta }_{{2p}}},\quad {{e}_{l}} = - \sum\limits_{k = 1}^m \,{{{v}}_{{0k}}}\int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{2}}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},\xi ,\eta )}}{{\partial x}}d\xi d\eta ,$

$A_{{lj}}^{{(6)}} = \int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial G_{1}^{'}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},s)}}{{\partial y}}ds,\quad B_{{lj}}^{{(6)}} = - \int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial {{G}_{1}}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},s)}}{{\partial y}}ds,$

$C_{{lj}}^{{(6)}} = \int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial G_{2}^{'}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},s)}}{{\partial y}}ds,\quad D_{{lj}}^{{(6)}} = - \int\limits_{{{\Gamma }_{j}}} \,\frac{{\partial {{G}_{2}}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},s)}}{{\partial y}}ds,$

$E_{{lkp}}^{{(6)}} = {{{v}}_{{pk}}}\int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{2}}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},\xi ,\eta )}}{{\partial y}}d\xi d\eta - 2\pi {{\delta }_{{lk}}}{{\delta }_{{3p}}},\quad {{h}_{l}} = - \sum\limits_{k = 1}^m \,{{{v}}_{{0k}}}\int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{2}}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},\xi ,\eta )}}{{\partial y}}d\xi d\eta ,$

$C_{{lj}}^{{(7)}} = A_{{lj}}^{{(5)}},\quad D_{{lj}}^{{(7)}} = B_{{lj}}^{{(5)}},$

$E_{{lkp}}^{{(7)}} = {{{v}}_{{pk}}}\int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{1}}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},\xi ,\eta )}}{{\partial x}}d\xi d\eta - 2\pi {{\delta }_{{lk}}}{{\delta }_{{5p}}},\quad {{q}_{l}} = - \sum\limits_{k = 1}^m \,{{{v}}_{{0k}}}\int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{1}}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},\xi ,\eta )}}{{\partial x}}d\xi d\eta ,$

$C_{{lj}}^{{(8)}} = A_{{lj}}^{{(6)}},\quad D_{{lj}}^{{(8)}} = B_{{lj}}^{{(6)}},$

$E_{{lkp}}^{{(8)}} = {{{v}}_{{pk}}}\int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{1}}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},\xi ,\eta )}}{{\partial y}}d\xi d\eta - 2\pi {{\delta }_{{lk}}}{{\delta }_{{6p}}},\quad {{t}_{l}} = - \sum\limits_{k = 1}^m \,{{{v}}_{{0k}}}\int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{1}}({{{\tilde {x}}}_{{cl}}},{{{\tilde {y}}}_{{cl}}},\xi ,\eta )}}{{\partial y}}d\xi d\eta ,$

где ${{\delta }_{{ij}}}$, ${{\delta }_{{lk}}}$, ${{\delta }_{{rp}}}$, $r = \overline {1,6} $ – символы Кронекера, а ${{v}_{{pk}}} = {{v}_{p}}({{\tilde {x}}_{{ck}}},{{\tilde {y}}_{{ck}}})$, $p = \overline {0,6} $. Для замыкания системы уравнений необходимо добавить к ней $2n$ соотношений из двух граничных условий вида (1.2), записанных в точках $({{x}_{{ci}}},{{y}_{{ci}}})$.

3. ВЫЧИСЛЕНИЕ ИНТЕГРАЛОВ

При численной реализации предложенного метода решения основную сложность и расчетное время как при составлении матрицы СЛАУ (2.16), (2.17), так и при нахождении искомых функций в расчетной точке с координатами $(x,y)$ по формулам (2.10)–(2.15) занимает вычисление интегралов от функций ${{G}_{1}}$, $G_{1}^{'}$, ${{G}_{2}}$, $G_{2}^{'}$, $\partial {{G}_{1}}{\text{/}}\partial x$, $\partial {{G}_{1}}{\text{/}}\partial y$, $\partial G_{1}^{'}{\text{/}}\partial x$, $\partial G_{1}^{'}{\text{/}}\partial y$, $\partial {{G}_{2}}{\text{/}}\partial x$, $\partial {{G}_{2}}{\text{/}}\partial y$, $\partial G_{2}^{'}{\text{/}}\partial x$, $\partial G_{2}^{'}{\text{/}}\partial y$ по линейным элементам ${{\Gamma }_{j}}$ и областным элементам ${{\Omega }_{k}}$. Аналитические формулы вычисления интегралов по линейным элементам получены в работе [12]. Для областных интегралов введем следующие обозначения:

$\begin{gathered} {{I}_{1}} = \int\limits_{{{\Omega }_{k}}} \,{{G}_{1}}d\xi d\eta ,\quad {{I}_{2}} = \int\limits_{{{\Omega }_{k}}} \,{{G}_{2}}d\xi d\eta ,\quad {{I}_{3}} = \int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{1}}}}{{\partial x}}d\xi d\eta , \\ {{I}_{4}} = \int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{1}}}}{{\partial y}}d\xi d\eta ,\quad {{I}_{5}} = \int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{2}}}}{{\partial x}}d\xi d\eta ,\quad {{I}_{6}} = \int\limits_{{{\Omega }_{k}}} \,\frac{{\partial {{G}_{2}}}}{{\partial y}}d\xi d\eta . \\ \end{gathered} $

Аналитическая формула для вычисления ${{I}_{1}}$ в случае, когда областным элементом является треугольник, получена в работе [10]. На основе подходов, использованных в этих работах, выведем аналитические формулы для вычисления интегралов ${{I}_{2}}$–${{I}_{6}}$ также для случая треугольного областного элемента.

Воспользуемся обозначениями, введенными в работе [10]. Обозначим комплексную координату расчетной точки $(x,y)$ через $z = x + iy$, комплексные координаты вершин треугольника ${{\Omega }_{k}}$ через $({{z}_{1}},{{z}_{2}},{{z}_{3}})$, комплексную координату текущей точки интегрирования $(\xi ,\eta )$ в треугольнике – $z{\kern 1pt} * = \xi + i\eta $. Порядок индексации вершин треугольника должен быть таков, чтобы при обходе вершин с возрастанием индекса внутренность треугольника оставалась слева. Для сокращения записи введем новую комплексную координату $\zeta = z{\kern 1pt} * - \;z$, т.е. перейдем в новую комплексную плоскость, в которой начало координат $\zeta = 0$ совпадает с расчетной точкой $z$. Тогда новые координаты вершин треугольника запишутся в виде (фиг. 1a)

${{\zeta }_{1}} = {{z}_{1}} - z,\quad {{\zeta }_{2}} = {{z}_{2}} - z,\quad {{\zeta }_{3}} = {{z}_{3}} - z.$

Так же обозначим

(3.1)

${{\zeta }_{{21}}} = {{\zeta }_{2}} - {{\zeta }_{1}} = {{l}_{1}}{{e}^{{i{{\gamma }_{1}}}}},\quad {{\zeta }_{{32}}} = {{\zeta }_{3}} - {{\zeta }_{2}} = {{l}_{2}}{{e}^{{i{{\gamma }_{2}}}}},\quad {{\zeta }_{{13}}} = {{\zeta }_{1}} - {{\zeta }_{3}} = {{l}_{3}}{{e}^{{i{{\gamma }_{3}}}}},$

где ${{l}_{1}}$, ${{l}_{2}}$, ${{l}_{3}}$, ${{\gamma }_{1}}$, ${{\gamma }_{2}}$, ${{\gamma }_{3}}$ – длины сторон треугольника и углы, которые они образуют с горизонталью соответственно. Во введенных обозначениях выражение для вычисления интеграла ${{I}_{1}}$ в случае, когда $z$ располагается снаружи треугольника, следующее [10]:

(3.2)

${{I}_{1}} = - \frac{1}{4}{\text{Re}}\left\{ {\left| {sin({{\gamma }_{2}} - {{\gamma }_{1}})} \right|{{e}^{{ - i({{\gamma }_{2}} + {{\gamma }_{1}})}}}\left( {\frac{{{{\zeta }_{{21}}}}}{{{{\zeta }_{{13}}}}}\mathop {[{{\zeta }^{2}}(2ln\zeta - 3)]}\nolimits_{{{\zeta }_{1}}}^{{{\zeta }_{3}}} + \mathop {[{{\zeta }^{2}}(2ln\zeta - 3)]}\nolimits_{{{\zeta }_{1}}}^{{{\zeta }_{2}}} } \right)} \right\},$

а когда $z$ совпадает с одним из углов, например, с ${{z}_{1}}$ имеем

(3.3)

${{I}_{1}} = \frac{1}{4}{\text{Re}}\left\{ {\left| {sin({{\gamma }_{2}} - {{\gamma }_{1}})} \right|{{e}^{{ - i({{\gamma }_{2}} + {{\gamma }_{1}})}}}{{\zeta }_{2}}\mathop {[\zeta (2ln\zeta - 3)]}\nolimits_{{{\zeta }_{2}}}^{{{\zeta }_{3}}} } \right\}.$

Здесь и ниже использовано обозначение, введенное в [10], [12]

$[f(\zeta )]_{{{{t}_{1}}}}^{{{{t}_{2}}}} = f({{t}_{2}}) - f({{t}_{1}}).$

Фиг. 1.

Интегрирование по треугольнику.

Вычислим интеграл ${{I}_{2}}$. Рассмотрим сначала случай, когда расчетная точка $z$ находится снаружи треугольника. Для выполнения интегрирования введем локальные оси координат: ось $\sigma $, направленную вдоль стороны ${{\zeta }_{1}}{{\zeta }_{2}}$ треугольника, и ось $\tau $, направленную параллельно стороне ${{\zeta }_{2}}{{\zeta }_{3}}$ треугольника, но проходящую через текущую точку интегрирования $\zeta $. Элементарной площадью интегрирования является параллелограмм со сторонами $d\sigma $ и $d\tau $ и углом $({{\gamma }_{2}} - {{\gamma }_{1}})$ между ними (фиг. 1б), тогда площадь этого элемента равна

(3.4)

$d\Omega = \left| {sin({{\gamma }_{2}} - {{\gamma }_{1}})} \right|d\tau d\sigma .$

Из (2.4) имеем

${{G}_{2}} = \frac{1}{4}\operatorname{Re} \{ \zeta \bar {\zeta }(\ln \zeta - 1)\} .$

Здесь и далее черта сверху означает комплексное сопряжение. Вывод аналитического выражения для вычисления интеграла ${{I}_{2}}$ проведем, следуя методике вычисления интеграла ${{I}_{1}}$ в [10]. Подставив последнее выражение в ${{I}_{2}}$ с учетом (3.4), запишем

(3.5)

${{I}_{2}} = \frac{1}{4}\left| {sin({{\gamma }_{2}} - {{\gamma }_{1}})} \right|\operatorname{Re} J,\quad J = \int\limits_0^{{{l}_{1}}} \int\limits_0^l \,\zeta \bar {\zeta }(\ln \zeta - 1)d\tau d\sigma ,$

где $l$ – длина отрезка интегрирования вдоль оси $\tau $. Так как при интегрировании вдоль оси $\tau $ дифференциал $d\zeta = {{e}^{{i{{\gamma }_{2}}}}}d\tau $, то перейдем во внутреннем интеграле к интегрированию по комплексной переменной $\zeta $:

$J = \int\limits_0^{{{l}_{1}}} \,{{e}^{{ - i{{\gamma }_{2}}}}}{{J}_{1}}d\sigma ,\quad {{J}_{1}} = \int\limits_{{{t}_{1}}}^{{{t}_{2}}} \,\zeta \bar {\zeta }(ln\zeta - 1)d\zeta .$

Подынтегральная функция в ${{J}_{1}}$ не является аналитической, поэтому, воспользовавшись способом вычисления подобного интеграла по частям, приведенным в [12], найдем

${{J}_{1}} = \frac{1}{4}\mathop {\left[ {{{\zeta }^{2}}\bar {\zeta }(2ln\zeta - 3) - \frac{{{{\zeta }^{3}}{{e}^{{ - 2i{{\gamma }_{2}}}}}}}{9}(6ln\zeta - 11)} \right]}\nolimits_{{{t}_{1}}}^{{{t}_{2}}} .$

Для комплексных координат ${{t}_{1}}$ и ${{t}_{2}}$ (начальной и конечной точек интегрирования на оси $\tau $) имеет место соотношение

${{t}_{1}} = {{\zeta }_{1}} + \sigma {{e}^{{i{{\gamma }_{1}}}}},\quad {{t}_{2}} = {{t}_{1}} + \frac{{{{\zeta }_{{32}}}}}{{{{\zeta }_{{21}}}}}({{t}_{1}} - {{\zeta }_{1}}),$

откуда

$d{{t}_{1}} = {{e}^{{i{{\gamma }_{1}}}}}d\sigma ,\quad d{{t}_{2}} = - \frac{{{{\zeta }_{{13}}}}}{{{{\zeta }_{{21}}}}}{{e}^{{i{{\gamma }_{1}}}}}d\sigma .$

При вычислении интеграла по переменной $\sigma $ перейдем к интегрированию по комплексным переменным ${{t}_{1}}$ и ${{t}_{2}}$. С учетом последних соотношений подставим найденное значение интеграла ${{J}_{1}}$ в $J$:

(3.6)

$J = - \frac{{{{e}^{{ - i({{\gamma }_{2}} + {{\gamma }_{1}})}}}}}{4}\left[ {\frac{{{{\zeta }_{{21}}}}}{{{{\zeta }_{{13}}}}}\left( {{{J}_{2}} - \frac{{{{e}^{{ - 2i{{\gamma }_{2}}}}}}}{9}{{J}_{3}}} \right) + {{J}_{4}} - \frac{{{{e}^{{ - 2i{{\gamma }_{2}}}}}}}{9}{{J}_{5}}} \right],$

где

${{J}_{2}} = \int\limits_{{{\zeta }_{1}}}^{{{\zeta }_{3}}} \,t_{2}^{2}{{\bar {t}}_{2}}(2ln{{t}_{2}} - 3)d{{t}_{2}},\quad {{J}_{3}} = \int\limits_{{{\zeta }_{1}}}^{{{\zeta }_{3}}} \,t_{2}^{3}(6ln{{t}_{2}} - 11)d{{t}_{2}},$

${{J}_{4}} = \int\limits_{{{\zeta }_{1}}}^{{{\zeta }_{2}}} \,t_{1}^{2}{{\bar {t}}_{1}}(2ln{{t}_{1}} - 3)d{{t}_{1}},\quad {{J}_{5}} = \int\limits_{{{\zeta }_{1}}}^{{{\zeta }_{2}}} \,t_{1}^{3}(6ln{{t}_{1}} - 11)d{{t}_{1}}.$

Вычислив эти интегралы

${{J}_{2}} = \frac{1}{9}\mathop {\left[ {{{\zeta }^{3}}\bar {\zeta }(6ln\zeta - 11) - \frac{{{{e}^{{ - 2i{{\gamma }_{3}}}}}}}{8}{{\zeta }^{4}}(12ln\zeta - 25)} \right]}\nolimits_{{{\zeta }_{1}}}^{{{\zeta }_{3}}} ,$

${{J}_{4}} = \frac{1}{9}\mathop {\left[ {{{\zeta }^{3}}\bar {\zeta }(6ln\zeta - 11) - \frac{{{{e}^{{ - 2i{{\gamma }_{1}}}}}}}{8}{{\zeta }^{4}}(12ln\zeta - 25)} \right]}\nolimits_{{{\zeta }_{1}}}^{{{\zeta }_{2}}} ,$

${{J}_{3}} = \frac{1}{8}[{{\zeta }^{4}}(12ln\zeta - 25)]_{{{{\zeta }_{1}}}}^{{{{\zeta }_{3}}}},\quad {{J}_{5}} = \frac{1}{8}[{{\zeta }^{4}}(12ln\zeta - 25)]_{{{{\zeta }_{1}}}}^{{{{\zeta }_{2}}}}$

и подставив их в (3.6), а затем полученное выражение в (3.5), в итоге найдем

(3.7)

$\begin{gathered} {{I}_{2}} = - \frac{1}{{144}}\left| {sin({{\gamma }_{2}} - {{\gamma }_{1}})} \right|\operatorname{Re} \left\{ {{{e}^{{ - i({{\gamma }_{2}} + {{\gamma }_{1}})}}}\left( {\frac{{{{\zeta }_{{21}}}}}{{{{\zeta }_{{13}}}}}{{J}_{6}} + {{J}_{7}}} \right)} \right\}, \\ {{J}_{6}} = [{{\zeta }^{3}}\bar {\zeta }(6ln\zeta - 11) - F({{\zeta }_{3}}){{\zeta }^{4}}(12ln\zeta - 25)]_{{{{\zeta }_{1}}}}^{{{{\zeta }_{3}}}}, \\ {{J}_{7}} = [{{\zeta }^{3}}\bar {\zeta }(6ln\zeta - 11) - F({{\zeta }_{2}}){{\zeta }^{4}}(12ln\zeta - 25)]_{{{{\zeta }_{1}}}}^{{{{\zeta }_{2}}}}, \\ \end{gathered} $

где

(3.8)

$F({{\zeta }_{1}}) = \frac{{{{e}^{{ - 2i{{\gamma }_{1}}}}} + {{e}^{{ - 2i{{\gamma }_{3}}}}}}}{8},\quad F({{\zeta }_{2}}) = \frac{{{{e}^{{ - 2i{{\gamma }_{2}}}}} + {{e}^{{ - 2i{{\gamma }_{1}}}}}}}{8},\quad F({{\zeta }_{3}}) = \frac{{{{e}^{{ - 2i{{\gamma }_{3}}}}} + {{e}^{{ - 2i{{\gamma }_{2}}}}}}}{8}.$

Если текущая точка $z$ совпадает с одним из углов треугольника, то формула вычисления интеграла упрощается. Для случая, когда $z$ совпадает с вершиной ${{z}_{1}}$, с учетом (3.1) имеем

(3.9)

${{\zeta }_{1}} = 0,\quad {{\zeta }_{{21}}} = {{\zeta }_{2}},\quad {{\zeta }_{{13}}} = - {{\zeta }_{3}}.$

Подставив эти значения в (3.7) и упростив выражение, с учетом (3.8) получим

(3.10)

${{I}_{2}} = \frac{1}{{144}}\left| {sin({{\gamma }_{2}} - {{\gamma }_{1}})} \right|\operatorname{Re} \left\{ {{{e}^{{ - i({{\gamma }_{2}} + {{\gamma }_{1}})}}}{{\zeta }_{2}}[{{\zeta }^{2}}\bar {\zeta }(6ln\zeta - 11) - {{\zeta }^{3}}(12ln\zeta - 25)F(\zeta )]_{{{{\zeta }_{2}}}}^{{{{\zeta }_{3}}}}} \right\}.$

Для применения этой формулы в случае совпадения $z$ с ${{z}_{2}}$ или с ${{z}_{3}}$ достаточно циклически переиндексировать вершины треугольника. Если текущая точка $z$ располагается на одной из сторон треугольника или внутри него, то этот треугольник необходимо разбить на два или три треугольника (см. фиг. 1в) и вычислить интеграл как сумму двух или трех интегралов с учетом формулы (3.10).

Для вычисления интегралов ${{I}_{3}}$–${{I}_{6}}$ продифференцируем формулы (3.2), (3.3), (3.7), (3.10) по $x$ и $y$ с учетом соотношений

$\frac{{d\zeta }}{{dx}} = - 1,\quad \frac{{d\overline \zeta }}{{dx}} = - 1,\quad \frac{{d\zeta }}{{dy}} = - i,\quad \frac{{d\overline \zeta }}{{dy}} = i.$

В итоге для случая, когда $z$ расположена снаружи треугольника, получим

${{I}_{3}} = \operatorname{Re} {{J}_{8}},\quad {{I}_{4}} = - \operatorname{Im} {{J}_{8}},$

${{I}_{5}} = \frac{1}{{144}}\left| {sin({{\gamma }_{2}} - {{\gamma }_{1}})} \right|{\text{Re}}\left\{ {{{e}^{{ - i({{\gamma }_{2}} + {{\gamma }_{1}})}}}\left( {\frac{{{{\zeta }_{{21}}}}}{{{{\zeta }_{{13}}}}}((1 - 8F({{\zeta }_{3}})){{J}_{9}} + {{J}_{{10}}}) + \left. {(1 - 8F({{\zeta }_{2}})){{J}_{{11}}} + {{J}_{{12}}}} \right)} \right)} \right\},$

${{I}_{6}} = - \frac{1}{{144}}\left| {sin({{\gamma }_{2}} - {{\gamma }_{1}})} \right|{\text{Im}}\left\{ {{{e}^{{ - i({{\gamma }_{2}} + {{\gamma }_{1}})}}}\left( {\frac{{{{\zeta }_{{21}}}}}{{{{\zeta }_{{13}}}}}((1 + 8F({{\zeta }_{3}})){{J}_{9}} - {{J}_{{10}}}) + (1 + 8F({{\zeta }_{2}})){{J}_{{11}}} - {{J}_{{12}}}} \right)} \right\},$

где

${{J}_{8}} = \left| {sin({{\gamma }_{2}} - {{\gamma }_{1}})} \right|{{e}^{{ - i({{\gamma }_{2}} + {{\gamma }_{1}})}}}\left( {\frac{{{{\zeta }_{{21}}}}}{{{{\zeta }_{{13}}}}}\mathop {\left[ {\zeta (ln\zeta - 1)} \right]}\nolimits_{{{\zeta }_{1}}}^{{{\zeta }_{3}}} + \mathop {\left[ {\zeta (ln\zeta - 1)} \right]}\nolimits_{{{\zeta }_{1}}}^{{{\zeta }_{2}}} } \right),$

${{J}_{9}} = [{{\zeta }^{3}}(6ln\zeta - 11)]_{{{{\zeta }_{1}}}}^{{{{\zeta }_{3}}}},\quad {{J}_{{10}}} = 9[{{\zeta }^{2}}\bar {\zeta }(2ln\zeta - 3)]_{{{{\zeta }_{1}}}}^{{{{\zeta }_{3}}}},$

${{J}_{{11}}} = [{{\zeta }^{3}}(6ln\zeta - 11)]_{{{{\zeta }_{1}}}}^{{{{\zeta }_{2}}}},\quad {{J}_{{12}}} = 9[{{\zeta }^{2}}\bar {\zeta }(2ln\zeta - 3)]_{{{{\zeta }_{1}}}}^{{{{\zeta }_{2}}}}.$

Для случая, когда $z$ совпадает с ${{z}_{1}}$, с учетом (3.9) имеем

${{I}_{3}} = - \operatorname{Re} {{J}_{{13}}},\quad {{I}_{4}} = \operatorname{Im} {{J}_{{13}}},$

${{I}_{5}} = - \frac{1}{{144}}\left| {sin({{\gamma }_{2}} - {{\gamma }_{1}})} \right|{\text{Re}}\left\{ {{{e}^{{ - i({{\gamma }_{2}} + {{\gamma }_{1}})}}}{{\zeta }_{2}}({{J}_{{14}}} + {{J}_{{15}}} - {{J}_{{16}}})} \right\},$

${{I}_{6}} = \frac{1}{{144}}\left| {sin({{\gamma }_{2}} - {{\gamma }_{1}})} \right|{\text{Im}}\left\{ {{{e}^{{ - i({{\gamma }_{2}} + {{\gamma }_{1}})}}}{{\zeta }_{2}}( - {{J}_{{14}}} + {{J}_{{15}}} - {{J}_{{16}}})} \right\},$

где

${{J}_{{13}}} = \left| {sin({{\gamma }_{2}} - {{\gamma }_{1}})} \right|{{e}^{{ - i({{\gamma }_{2}} + {{\gamma }_{1}})}}}{{\zeta }_{2}}[ln\zeta ]_{{{{\zeta }_{2}}}}^{{{{\zeta }_{3}}}},$

${{J}_{{14}}} = [{{\zeta }^{2}}(6ln\zeta - 11)]_{{{{\zeta }_{2}}}}^{{{{\zeta }_{3}}}},\quad {{J}_{{15}}} = 9[\zeta \bar {\zeta }(2ln\zeta - 3)]_{{{{\zeta }_{2}}}}^{{{{\zeta }_{3}}}},$

${{J}_{{16}}} = 8[{{\zeta }^{2}}(6ln\zeta - 11)F(\zeta )]_{{{{\zeta }_{2}}}}^{{{{\zeta }_{3}}}}.$

Для вычисления интегралов ${{I}_{1}}{\kern 1pt} - {\kern 1pt} {{I}_{6}}$ по четырехугольному элементу достаточно разбить его на два треугольника по любой из диагоналей.

4. ТЕСТОВЫЙ РАСЧЕТ

В тестовом расчете в качестве области $\Omega = \{ 0 \leqslant x \leqslant 1;0 \leqslant y \leqslant 1\} $ возьмем квадрат с единичной стороной, левый нижний угол которого совпадает с началом координат. Каждую из сторон квадрата разобьем на равное количество ${{n}_{1}}$ линейных элементов одинаковой длины, тогда общее число линейных элементов будет равно $n = 4{{n}_{1}}$. Расчетную область покроем равномерной квадратной сеткой размерности $({{n}_{1}} \times {{n}_{1}})$, состоящей из одинаковых четырехугольных элементов. Таким образом, общее количество площадных элементов равно $m = n_{1}^{2}$. Для оценки точности введем величины абсолютной и относительной погрешностей:

(4.1)

$\varepsilon _{a}^{f} = \mathop {max}\limits_k |f({{x}_{k}},{{y}_{k}}) - {{f}_{a}}({{x}_{k}},{{y}_{k}})|,$

(4.2)

$\varepsilon _{e}^{f} = \frac{1}{K}\sum\limits_k \,\left| {\frac{{f({{x}_{k}},{{y}_{k}}) - {{f}_{a}}({{x}_{k}},{{y}_{k}})}}{{{{f}_{{max}}}}}} \right|,$

где ${{f}_{a}}(x,y)$ – аналитическое решение, а ${{f}_{{max}}} = \mathop {max}\limits_{(x,y) \in \Omega } {{f}_{a}}(x,y)$. В качестве набора контрольных точек $({{x}_{k}},{{y}_{k}})$, $k = \overline {1,K} $, выберем точки, лежащие в углах областных элементов и образующие равномерную квадратную сетку размерности $(({{n}_{1}} + 1) \times ({{n}_{1}} + 1))$, покрывающую область $\Omega $. Аналогично введем абсолютную $\varepsilon _{a}^{g}$ и относительную $\varepsilon _{e}^{g}$ погрешности для функции $g(x,y)$.

В тестовой серии расчетов проведено решение неоднородного бигармонического уравнения (0.1) с правой частью (1.1), в которой функции ${{v}_{i}} = {{v}_{i}}(x,y)$, $i = \overline {0,6} $, приняты равными

$\begin{gathered} {{v}_{0}} = {{x}^{2}}{{y}^{2}}cosx,\quad {{v}_{1}} = x,\quad {{v}_{2}} = sinx,\quad {{v}_{3}} = y, \\ {{v}_{4}} = cosx,\quad {{v}_{5}} = {{x}^{2}}y,\quad {{v}_{6}} = xy \\ \end{gathered} $

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

${{f}_{a}}(x,y) = - {{g}_{a}}(x,y) = ysinx.$

В качестве граничных зададим следующие условия:

${{\left. f \right|}_{\Gamma }} = {{f}_{a}}({{x}_{1}}(s),{{y}_{1}}(s)),\quad {{\left. {\Delta f} \right|}_{\Gamma }} = {{g}_{a}}({{x}_{1}}(s),{{y}_{1}}(s)).$

В ходе тестовых расчетов определялись значения погрешностей $\varepsilon _{a}^{f}$, $\varepsilon _{a}^{g}$, $\varepsilon _{e}^{f}$, $\varepsilon _{e}^{g}$ в зависимости от параметра ${{n}_{1}}$, определяющего детальность разбиения границы $\Gamma $ и области $\Omega $ на линейные и областные элементы. На фиг. 2 представлены результаты вычислений для значений ${{n}_{1}} = 10,\;20,\;30,\;40,\;50$. Как видно из фигуры, абсолютная и относительная погрешности для функций $f(x,y)$ и $g(x,y)$ не велики и уменьшаются с ростом числа разбиений, что свидетельствует о хорошей точности и сеточной сходимости разработанного метода.

Фиг. 2.

Зависимости погрешностей в тестовом расчете от ${{n}_{1}}$.

5. РАСЧЕТ ФИЛЬТРАЦИОННОГО ТЕЧЕНИЯ В РАМКАХ МОДЕЛИ БРИНКМАНА В НЕОДНОРОДНОЙ ПОРИСТОЙ СРЕДЕ

Продемонстрируем возможности предложенного метода на примере двумерной задачи о фильтрации жидкости в пористой области с неоднородным распределением проницаемости в рамках модели Бринкмана [13]

(5.1)

$ - \nabla p - \frac{\mu }{k}{\mathbf{u}} + {{\mu }_{b}}\Delta {\mathbf{u}} = 0,$

(5.2)

$\operatorname{div} {\mathbf{u}} = 0,$

где ${\mathbf{u}} = ({{u}_{x}},{{u}_{y}})$ – вектор скорости фильтрации жидкости, $p$ – давление, $\mu $ – вязкость жидкости, ${{\mu }_{b}} = {{\mu }_{b}}(x,y)$ – эффективная вязкость фильтрационного потока, $k = k(x,y)$ – проницаемость пористой среды. В качестве расчетной выберем прямоугольную область $\Omega = \{ 0 \leqslant x \leqslant L;0 \leqslant y \leqslant H\} $ (фиг. 3a). Будем считать, что жидкость течет в области $\Omega $ слева направо, поступая через отрезок $AD$ и выходя через отрезок $BC$, а горизонтальные участки границы $AB$ и $CD$ являются линиями тока. Расход потока задан и равен $Q$.

Фиг. 3.

Расчетные области $\Omega $ (a) и ${{\Omega }_{{v}}}$ (б).

На основании уравнения неразрывности (5.2) введем функцию тока $\psi $ и завихренности $\omega $ равенствами

(5.3)

${{u}_{x}} = \frac{{\partial \psi }}{{\partial y}},\quad {{u}_{y}} = - \frac{{\partial \psi }}{{\partial x}},\quad \omega = \frac{{\partial {{u}_{y}}}}{{\partial x}} - \frac{{\partial {{u}_{x}}}}{{\partial y}} = - \Delta \psi .$

Выбрав в качестве характерных размерных величин высоту $H$ области и среднюю скорость потока $U = Q{\text{/}}H$, перейдем к безразмерным переменным

(5.4)

$\begin{gathered} \bar {x} = \frac{x}{H},\quad \bar {y} = \frac{y}{H},\quad \bar {u} = \frac{u}{U},\quad \bar {\psi } = \frac{\psi }{{UH}},\quad \bar {\omega } = \frac{H}{U}\omega , \\ \bar {p} = \frac{H}{{U\mu }}p,\quad \mathop {\bar {\mu }}\nolimits_b = \frac{{{{\mu }_{b}}}}{\mu },\quad \bar {k} = \frac{k}{{{{H}^{2}}}}. \\ \end{gathered} $

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

Перепишем закон Бринкмана (5.1) с учетом (5.3) в виде

$\frac{{\partial p}}{{\partial x}} = - \frac{1}{k}\frac{{\partial \psi }}{{\partial y}} - {{\mu }_{b}}\frac{{\partial \omega }}{{\partial y}},\quad \frac{{\partial p}}{{\partial y}} = \frac{1}{k}\frac{{\partial \psi }}{{\partial x}} + {{\mu }_{b}}\frac{{\partial \omega }}{{\partial x}}.$

Исключив из этих уравнений давление $p$ и выразив член со старшей производной, получим

(5.5)

${{\Delta }^{2}}\psi = - \frac{1}{{{{k}^{2}}{{\mu }_{b}}}}\left( {\frac{{\partial k}}{{\partial x}}\frac{{\partial \psi }}{{\partial x}} + \frac{{\partial k}}{{\partial y}}\frac{{\partial \psi }}{{\partial y}}} \right) - \frac{1}{{k{{\mu }_{b}}}}\omega + \frac{1}{{{{\mu }_{b}}}}\left( {\frac{{\partial {{\mu }_{b}}}}{{\partial x}}\frac{{\partial \omega }}{{\partial x}} + \frac{{\partial {{\mu }_{b}}}}{{\partial y}}\frac{{\partial \omega }}{{\partial y}}} \right).$

Это уравнение является неоднородным бигармоническим уравнением вида (0.1), где с учетом (1.1)

$\begin{gathered} f = \psi ,\quad g = - \omega ,\quad {{v}_{0}} = {{v}_{1}} = 0,\quad {{v}_{2}} = - \frac{1}{{{{k}^{2}}{{\mu }_{b}}}}\frac{{\partial k}}{{\partial x}},\quad {{v}_{3}} = - \frac{1}{{{{k}^{2}}{{\mu }_{b}}}}\frac{{\partial k}}{{\partial y}}, \\ {{v}_{4}} = \frac{1}{{k{{\mu }_{b}}}},\quad {{v}_{5}} = - \frac{1}{{{{\mu }_{b}}}}\frac{{\partial {{\mu }_{b}}}}{{\partial x}},\quad {{v}_{6}} = - \frac{1}{{{{\mu }_{b}}}}\frac{{\partial {{\mu }_{b}}}}{{\partial y}}. \\ \end{gathered} $

Отметим, что если считать пористую среду однородной, то ${{v}_{2}} = {{v}_{3}} = {{v}_{5}} = {{v}_{6}} \equiv 0$ и уравнение (5.5) перейдет в известное уравнение Бринкмана

${{\Delta }^{2}}\psi - {{S}^{2}}\Delta \psi = 0,$

где ${{S}^{2}} = {{(k{{\mu }_{b}})}^{{ - 1}}}$.

Искомая функция $\psi (x,y)$ удовлетворяет в области $\Omega $ уравнению (5.5), а на ее границе следующим граничным условиям. На входном участке $AD$ границы

(5.6)

$\psi = q(y),\quad \frac{{\partial \psi }}{{\partial n}} = 0,$

что соответствует заданному профилю скорости потока ${{u}_{x}}(0,y) = \frac{{dq}}{{dy}}$, ${{u}_{y}}(0,y) \equiv 0$, где $q(y)$ – известная монотонно возрастающая функция, удовлетворяющая условиям $q(0) = 0$, $q(1) = 1$. На выходном участке $BC$ границы зададим мягкие граничные условия

(5.7)

$\frac{{\partial \psi }}{{\partial n}} = 0,\quad \frac{{\partial \omega }}{{\partial n}} = 0.$

На участках $AB$ и $CD$ потребуем выполнения условий непротекания

(5.8)

$\psi = 0,\quad \omega = 0,$

(5.9)

$\psi = 1,\quad \omega = 0$

соответственно. Для решения краевой задачи (5.5)–(5.7) используем метод, описанный в разд. 2.

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

(5.10)

$ - \nabla p + \Delta {\mathbf{v}} = 0,$

(5.11)

$\operatorname{div} {\mathbf{v}} = 0,$

где ${\mathbf{v}} = ({{v}_{x}},{{v}_{y}})$ – истинная скорость вязкой жидкости. На основании равенства (5.11), записав соотношения, аналогичные (5.3), введем функцию тока ${{\psi }_{s}}$ и завихренности ${{\omega }_{s}}$ вязкого течения. Из (5.10) следует, что функция тока удовлетворяет однородному бигармоническому уравнению

(5.12)

${{\Delta }^{2}}{{\psi }_{s}} = 0.$

Рассмотрим модельную пористую среду, представляющую собой периодическую структуру продублированных в вертикальном и горизонтальном направлениях элементарных ячеек (ЭЯ), каждая из которых содержит твердое включение. Для настоящего расчета будем полагать, что ЭЯ имеет форму квадрата со стороной $h$, в центре которого расположено круговое включение радиуса $r$. Такая пористая среда аналогична модельной пористой среде конфигурации S1 в работе [14], однако, радиусы круговых включений в ЭЯ будем полагать различными, что позволит моделировать заданное неоднородное распределение проницаемости $k(x,y)$ пористой среды. Таким образом, расчетная область $\Omega $ полностью покрыта квадратной сеткой ЭЯ размерности $({{n}_{x}} \times {{n}_{y}})$, где ${{n}_{x}} = L{\text{/}}h$, ${{n}_{y}} = H{\text{/}}h$ (фиг. 3б). Для каждой ЭЯ локальное значение объемной концентрации твердых включений определим по формуле

(5.13)

$\phi = \frac{{\pi {{r}^{2}}}}{{{{h}^{2}}}},$

тогда локальное значение пористости $m = 1 - \phi $.

Дадим следующую математическую постановку задачи вычислительного эксперимента. Требуется определить функцию тока ${{\psi }_{s}}(x,y)$, удовлетворяющую уравнению (5.12) в области ${{\Omega }_{v}}$ – межпоровом пространстве пористой среды в области $\Omega $. На участках внешней границы области ${{\Omega }_{v}}$ потребуем выполнения граничных условий (5.6), (5.7). На поверхности каждого твердого включения ${{\ell }_{k}}$, $k = \overline {1,{{n}_{x}}{{n}_{y}}} $, зададим условия прилипания

${{\psi }_{s}} = {{\psi }_{k}},\quad \frac{{\partial {{\psi }_{s}}}}{{\partial n}} = 0,$

где ${{\psi }_{k}}$ – неизвестное заранее значение функции тока на поверхности $k$-го включения. Для решения этой задачи используем метод граничных элементов (МГЭ) по аналогии с тем, как это было сделано в работах [12], [14]. Для твердых включений малого размера, для которых вычисленное по формуле (5.13) значение концентрации $\phi < 0.01$, воспользуемся подходом МГЭ с представлением неизвестных функций на границе включений усеченным рядом Фурье, описанным в [15]. Это позволит снизить количество искомых неизвестных без потери точности расчета.

По рассчитанной в области ${{\Omega }_{{v}}}$ функции ${{\psi }_{s}}(x,y)$ построим осредненную функцию ${{\tilde {\psi }}_{s}}(x,y)$, определенную во всех точках области $\Omega $, способом, описанным в [15]. Значения функции ${{\tilde {\psi }}_{s}}(x,y)$ возьмем равными значениям ${{\psi }_{s}}(x,y)$ в угловых точках и серединах сторон всех ЭЯ, а в их центрах – соответствующим значениям ${{\psi }_{k}}$ на поверхностях твердых включений, содержащихся в каждой ЭЯ. Перечисленные точки покрывают всю область $\Omega $ равномерной квадратной сеткой. В остальных точках области $\Omega $ осредненная функция тока строится билинейной интерполяцией по найденным значениям в узлах этой сетки. Будем называть решение в рамках микроскопического подхода вычислительным экспериментом.

Для сравнения решим также поставленную задачу в рамках модели Дарси. Уравнение для функции тока ${{\psi }_{d}}$ фильтрационного течения по этой модели в случае неоднородного распределения проницаемости приведено в работе [11] и во введенных безразмерных переменных (5.4) имеет вид

(5.14)

$\Delta {{\psi }_{d}} = \frac{1}{k}\left( {\frac{{\partial k}}{{\partial x}}\frac{{\partial {{\psi }_{d}}}}{{\partial x}} + \frac{{\partial k}}{{\partial y}}\frac{{\partial {{\psi }_{d}}}}{{\partial y}}} \right).$

Математическая постановка задачи следующая. Требуется определить функцию тока ${{\psi }_{d}}(x,y)$, удовлетворяющую в области $\Omega $ уравнению (5.14), а на ее границе граничным условиям:

$AB\,:\;\;{{\psi }_{d}} = 0,\quad BC\,:\;\;\frac{{\partial {{\psi }_{d}}}}{{\partial n}} = 0,\quad CD\,:\;\;{{\psi }_{d}} = 1,\quad AD\,:\;\;{{\psi }_{d}} = q(y).$

Уравнение (5.14) представляет собой уравнение Пуассона с неизвестными функциями в правой части. Для его решения также воспользуемся МГОЭ, изложение которого для решения этого уравнения приведено в [11].

Зададим неоднородное распределение проницаемости в виде следующей зависимости:

(5.15)

$k(x,y) = {{k}_{{min}}} + {{k}_{0}}\left( {\frac{x}{L},\frac{y}{H}} \right)({{k}_{{max}}} - {{k}_{{min}}}),$

где ${{k}_{{min}}}$ и ${{k}_{{max}}}$ – заданные минимальное и максимальное значения проницаемости в области $\Omega $, а

${{k}_{0}}(\xi ,\eta ) = 1 - \varphi (\xi )\varphi (\eta ),\quad \varphi (\xi ) = \left\{ \begin{gathered} {{2}^{{p - 1}}}{{\xi }^{p}},\quad 0 \leqslant \xi \leqslant 0.5, \hfill \\ 1 - {{2}^{{p - 1}}}{{(1 - \xi )}^{p}},\quad 0.5 < \xi \leqslant 1. \hfill \\ \end{gathered} \right.$

График функции $\varphi (\xi )$ для $p = 4$ и изолинии функции ${{k}_{0}}(\xi ,\eta )$ приведены на фиг. 4.

Фиг. 4.

График функции $\varphi (\xi )$ и изолинии функции ${{k}_{0}}(\xi ,\eta )$ для $p = 4$.

Для проведения расчетов в рамках микроскопического подхода необходимо задать радиус твердого включения в каждой ЭЯ, исходя из локального значения проницаемости $k({{x}_{i}},{{y}_{j}})$, вычисленного по формуле (5.15) в центре ЭЯ с координатами ${{x}_{i}} = (i - 0.5)h$, $i = \overline {1,{{n}_{x}}} $, ${{y}_{j}} = (j - 0.5)h$, $j = \overline {1,{{n}_{y}}} $. Для высокопористых сред, пористость $m$ которых близка к единице, значение проницаемости с большой степенью точности можно определить по формуле из работы [16], которая во введенных безразмерных переменных с учетом (5.13) имеет вид

(5.16)

$k(\phi ) = \frac{1}{{8\pi {{\kappa }^{2}}}}\left( { - ln\phi - 1.476 + 2\phi - 1.774{{\phi }^{2}} + 4.076{{\phi }^{3}}} \right),$

где $\kappa = H{\text{/}}h$ – безразмерный геометрический параметр, определяющий отношение характерного линейного размера задачи к размеру ЭЯ. Из выражения (5.13) получим формулу для вычисления радиуса твердого включения

$r(k) = h\sqrt {\frac{{\phi (k)}}{\pi }} ,$

где $\phi (k)$ – функция, обратная к (5.16).

Для вычисления эффективной вязкости воспользуемся формулой

${{\mu }_{b}}(k) = \mathop {\left[ {1 + \mathop {\left( {\frac{{k{{\kappa }^{2}}}}{{\lambda _{*}^{2}}}} \right)}\nolimits^{ - B} } \right]}\nolimits^{1/B} ,$

полученной в работе [14], где для рассматриваемой конфигурации модельной пористой среды $B = 2.46$, ${{\lambda }_{*}} = 0.375$. Тогда имеем

$\frac{{\partial {{\mu }_{b}}}}{{\partial x}} = \frac{{d{{\mu }_{b}}}}{{dk}}\frac{{\partial k}}{{\partial x}},\quad \frac{{\partial {{\mu }_{b}}}}{{\partial y}} = \frac{{d{{\mu }_{b}}}}{{dk}}\frac{{\partial k}}{{\partial y}}.$

При проведении расчетов положим $L = 1$, т.е. ${{n}_{x}} = {{n}_{y}}$. Сетка граничных и областных элементов для всех трех моделей построим аналогично тому, как это было сделано в тестовых расчетах в разд. 4, положив ${{n}_{1}} = 20,\;40,\;60$. В вычислительном эксперименте каждую сторону внешней границы разобьем на ${{n}_{1}} = 100$ линейных элементов, а границы крупных твердых включений, расположенных в ЭЯ с локальным значением концентрации $\phi \geqslant 0.01$, аппроксимируем ${{n}_{2}} = 30$ линейными элементами. Значения ${{k}_{{min}}}$ и ${{k}_{{max}}}$ определим по формуле (5.16) для фиксированных предельных значений концентраций ${{\phi }_{{min}}} = {{10}^{{ - 6}}}$, ${{\phi }_{{max}}} = 0.1$ и различных $\kappa = {{n}_{y}}$.

Для оценки точности расчетов найдем погрешности $\varepsilon _{a}^{\psi }$, $\varepsilon _{a}^{{\psi d}}$ $\varepsilon _{e}^{\psi }$ и $\varepsilon _{e}^{{\psi d}}$ по формулам (4.1), (4.2), где в качестве численного решения положим $f = \psi (x,y)$ (модель Бринкмана) и $f = {{\psi }_{d}}(x,y)$ (модель Дарси) соответственно. В качестве тестового решения примем ${{f}_{a}} = \mathop {\tilde {\psi }}\nolimits_s (x,y)$, полученную из вычислительного эксперимента. При этом максимальное значение функции тока для всех моделей в силу постановки задачи ${{\psi }_{{max}}} = 1$.

На фиг. 5 приведены результаты расчетов для ${{n}_{y}} = 5$ (сверху) и ${{n}_{y}} = 10$ (снизу) в случае ${{n}_{1}} = 40$. На левых графиках представлены результаты вычислительных экспериментов. Сплошными линиями показаны линии тока детального течения вязкой жидкости в межпоровом пространстве, а штриховыми линиями – изолинии осредненной функции тока $\mathop {\tilde {\psi }}\nolimits_s (x,y)$. Твердые включения изображены серым цветом, а их положения в случае их малого размера помечены символом “$ \times $”. Из фигур видно, что изолинии осредненной функции тока повторяют общее поведение линий тока детального расчета, а в областях с высокой проницаемостью (т.е. с низкой концентрацией твердых включений) практически совпадают с ними. На правых графиках приведено сравнение изолиний осредненной функции тока $\mathop {\tilde {\psi }}\nolimits_s (x,y)$ (сплошная линия) с линиями тока фильтрационных течений по моделям Бринкмана (штриховые линии) и Дарси (пунктирные линии), построенных как изолинии функций $\psi (x,y)$ и ${{\psi }_{d}}(x,y)$ соответственно. Из фигур видно, что результаты расчета по модели Бринкмана с большей степенью точности совпадают с результатами вычислительного эксперимента, чем результаты по модели Дарси. Отметим, что линии тока, построенные по модели Бринкмана, визуально совпадают с изолиниями функции $\mathop {\tilde {\psi }}\nolimits_s $ практически всюду за исключением небольшой области больших градиентов функции $k(x,y)$ (см. фиг. 4) при ${{n}_{y}} = 5$, а при ${{n}_{y}} = 10$ – во всей расчетной области. Линии тока по модели Дарси не совпадают с изолиниями $\mathop {\tilde {\psi }}\nolimits_s $ во всей области $\Omega $, однако, различие становится меньше c увеличением ${{n}_{y}}$.

Фиг. 5.

Сравнение линий тока течения по модели Дарси и Бринкмана с результатами вычислительного эксперимента для ${{n}_{y}} = 5,\;10$: 1 – вычислительный эксперимент, 2 – изолинии осредненной функции тока в вычислительном эксперименте, 3 – модель Бринкмана, 4 – модель Дарси.

В табл. 1 представлены значения погрешностей $\varepsilon _{a}^{\psi }$, $\varepsilon _{a}^{{\psi d}}$ $\varepsilon _{e}^{\psi }$ и $\varepsilon _{e}^{{\psi d}}$ для значений ${{n}_{y}} = 5,\;10,\;15,\;20$. Видно, что с увеличением ${{n}_{y}}$ (т.е. с увеличением $\kappa $) погрешности по моделям Дарси и Бринкмана уменьшаются. При этом погрешность модели Бринкмана во всех расчетах в несколько раз меньше погрешности модели Дарси. Аналогичные выводы были сделаны и в работе [14] для случая однородной пористой среды. Так же отметим, что при увеличении разбиений ${{n}_{1}}$ погрешность модели Бринкмана уменьшается, что свидетельствует о сеточной сходимости разработанного метода, а погрешность модели Дарси увеличивается, что свидетельствует о том, что эта погрешность является погрешностью модели, а не вычислительной погрешностью.

Таблица 1.

Погрешности $\varepsilon _{a}^{\psi }$, $\varepsilon _{a}^{{\psi d}}$, $\varepsilon _{e}^{\psi }$ и $\varepsilon _{e}^{{\psi d}}$

${{n}_{1}}$	${{n}_{y}}$	$\varepsilon _{a}^{\psi }$	$\varepsilon _{a}^{{\psi d}}$	$\varepsilon _{e}^{\psi }$	$\varepsilon _{e}^{{\psi d}}$
20	5	0.819e–2	0.281e–1	0.257e–2	0.155e–1
20	10	0.535e–2	0.111e–1	0.139e–2	0.466e–2
20	15	0.401e–2	0.801e–2	0.956e–3	0.253e–2
20	20	0.360e–2	0.647e–2	0.818e–3	0.164e–2
40	5	0.731e–2	0.288e–1	0.217e–2	0.160e–1
40	10	0.373e–2	0.112e–1	0.927e–3	0.500e–2
40	15	0.222e–2	0.808e–2	0.541e–3	0.280e–2
40	20	0.191e–2	0.657e–2	0.451e–3	0.186e–2
60	5	0.707e–2	0.290e–1	0.206e–2	0.161e–1
60	10	0.329e–2	0.113e–1	0.803e–3	0.508e–2
60	15	0.170e–2	0.809e–2	0.417e–3	0.287e–2
60	20	0.137e–2	0.659e–2	0.331e–3	0.193e–2

ЗАКЛЮЧЕНИЕ

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

Авторы благодарят Ш.Х. Зарипова и В.Ф. Шарафутдинова за полезные замечания.

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

Тимошенко С.П., Войновский-Кригер С. Пластины и оболочки. М.: Наука, 1966. 636 с.
Хаппель Дж., Бреннер Г. Гидродинамика при малых числах Рейнольдса. М.: Мир, 1976. 630 с.
Флетчер К., Бреннер Г. Вычислительные методы в динамике жидкостей: в 2 т. М.: Мир, 1991. Т. 2. 504 с.
Алгазин С.Д. Численные алгоритмы классической математической физики. М.: Диагол–МИФИ, 2010. 240 с.
Wu J. Problem of General Visous Flow. Developments in BEM. London: Elsevier Applied Science Publication, 1982. V. 2.
Camp C.V., Gipson G.S. Boundary element analysis of nonhomogeneous biharmonic phenomena. Berlin, Heidelberg: Springer, 1992. 268 p.
Роуч П. Вычислительная гидродинамика. М.: Мир, 1980. 618 с.
Андерсон Д., Таннехил Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен: в 2 т. М.: Мир, 1990. 384 с.
Skerget L., Hribersek M., Kuhn G. Computational fluid dynamics by boundary–domain integral method // International Journal for Numerical Methods in Engng. 1999. V. 46. 8. P. 1291–1311.
Mardanov R.F., Zaripov S.K. Solution of Nonhomogeneous Helmholtz Equation with Variable Coefficient Using Boundary Domain Integral Method // Lobachevskii Journal of Math. 2018. V. 39. 6. P. 783–793.
Mardanov R.F., Sharafutdinov V.F., Ibragimov I.Z., Zaripov S.K., Baganina A.E. Solving the fluid flow problems with boundary domain integral method // J. Phys.: Conf. Ser. 2019. V. 1158. 032027.
Mardanov R.F., Dunnett S.J., Zaripov S.K. Modeling of fluid flow in periodic cell with porous cylinder using a boundary element method // Engng Analysis with Boundary Elements. 2016. V. 68. P. 54–62.
Brinkman H.C. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles // Appl. Sci. Res. 1947. V. A1. 1. P. 27–34.
Zaripov S.K., Mardanov R.F., Sharafutdinov V.F. Determination of Brinkman model parameters using Stokes flow model // Transport in Porous Media. 2019. V. 130. 2. P. 529–557.
Mardanov R.F., Zaripov S.K., Maklakov D.V. Two-dimensional Stokes flows in porous medium composed of a large number of circular inclusions // Engng Analysis with Boundary Elements. 2020. V. 113. P. 204–218.
Sangani A.S., Acrivos A. Slow flow through a periodic array of spheres // Internat. Journal of Multiphase Flow. 1982. V. 8. 4. P. 343–360.

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

Инструменты

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

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

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