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

1.1. Гидродинамическая постановка задачи

Рассматривается задача об исследовании движения ограниченной жидкости, когда инерционными силами можно пренебречь (приближение Стокса [11], [32]). Пусть жидкость заключена между двумя цилиндрическими телами с сечением произвольной формы. Тогда задача сводится к краевой задаче для бигармонического уравнения в плоской двусвязной области $D$, граница которой $\partial D$ состоит из внутреннего контура $\partial {{D}_{0}}$ и внешнего $\partial {{D}_{1}}$ (см. фиг. 1).

Точки границы контура движутся со скоростью ${\mathbf{V}}({{V}_{x}},{{V}_{y}})$, компоненты которой могут быть выбраны произвольно. Если контуры движутся как твердые тела, то компоненты скорости граничных точек складываются из компонент поступательной $({{V}_{x}}^{k},{{V}_{y}}^{k})$ и вращательной угловой скорости ${{\omega }_{k}}$:

где $\rho _{x}^{k}$, $\rho _{y}^{k}$ – координаты радиус-вектора точки границы с началом в точке, через которую проходит ось вращения.

Уравнения движения несжимаемой жидкости в приближении Стокса [6], [20] имеют вид:

где $p$ – гидродинамическое давление, $\mu $ – коэффициент вязкости среды.

Для плоской задачи Лагранж ввел функцию тока $\Psi (x,y)$ (см. [24]) и выразил через нее компоненты скорости жидкости, пользуясь условием несжимаемости:

Тогда уравнения движения жидкости (уравнения Стокса) примут вид:

Если теперь продифференцировать первое уравнение (3) по $y$, а второе – по $x$ и записать разность полученных выражений, то видно, что функция тока удовлетворяет бигармоническому уравнению ${{\Delta }^{2}}\Psi = 0$. Это наблюдение использовалось при построении точных решений самим Стоксом [11], Н.Е. Жуковским и С.А. Чаплыгиным [3]–[5] и в дальнейшем устойчиво вошло в практику и широко применяется в классической литературе по гидродинамике вязкой жидкости [32]–[34].

Граничные условия для определения функции тока получаются из условия прилипания, в соответствии с которым на контурах $\partial {{D}_{0}}$ и $\partial {{D}_{1}}$ скорость жидкости совпадает со скоростью точки контуров:

1.2. Граничные условия для функции тока

Общая краевая задача имеет вид

Для компонент нормали имеем соотношения ${{n}_{x}} = \frac{{\partial y}}{{\partial s}},\quad {{n}_{y}} = - \frac{{\partial x}}{{\partial s}}$, где $x(s)$, $y(s)$ – параметрические уравнения границ $\partial {{D}_{k}}$, $k = 0,1$, $s$ – дуговая координата и меняется от 0 до ${{l}_{0}}$ на внутреннем контуре и от ${{l}_{0}}$ до ${{l}_{0}} + {{l}_{1}}$ на внешнем, ${{l}_{0}},{{l}_{1}}$ – длины контуров; обход контура выбирается так, чтобы область $D$ оставалась слева (фиг. 1). С помощью этих соотношений проекции скорости на нормаль ${{V}_{n}}$ и касательную ${{V}_{s}}$ можно выразить через декартовы компоненты скорости

От рассматриваемой гидродинамической задачи можно перейти к задаче об определении функции тока $\Psi $ по известным граничным значениям самой функции и ее нормальной производной. Действительно, функция ${{V}_{n}}(s)$ задана на границе. Поэтому первое граничное условие (5) можно проинтегрировать и тогда краевая задача (5) примет вид

В таком виде обычно формулируется общая краевая задача для бигармонического уравнения и для нее доказана теорема единственности решения (см. [28, c. 399–400]). Заметим, что функция ${{g}^{{(k)}}}(s)$ определена с точностью до произвольной постоянной. Поэтому при заданных функциях ${{g}^{{(k)}}}(s)$ и $g_{n}^{k}(s)$ решение задачи (6) единственно, а при заданных ${{V}_{n}}$ и ${{V}_{s}}$ функция $\Psi $из решения задачи (5) определяется с точностью до произвольной постоянной. Аналогично и при заданных ${{V}_{x}}$ и ${{V}_{y}}$ функция $\Psi $ из решения задачи (4) также определяется с точностью до произвольной постоянной.

Особый интерес представляет частный случай движения твердых поверхностей. Найдем краевые условия (5) на поверхностях $\partial {{D}_{k}},\;k = 0,1\;$, движущихся, как твердые тела. Подставляя в (6) выражения (1), найдем функции ${{g}^{{(k)}}}(s)$ и $g_{n}^{k}(s)$, входящие в граничные условия (5):

1.3. Представление Гурса и определение констант

Справедливо следующее утверждение (теорема Гурса): бигармоническая в плоской области $D$ функция $\Psi $ может быть представлена через две независимые гармонические в этой области функции $\varphi $ и $\psi $:

Бигармоничность такого представления функции $\Psi $ очевидна. Достаточные условия существования функций $\varphi (x,y)$ и $\psi (x,y)$ для односвязных областей сформулированы, например, в [28, с. 400–401]. Ниже будет дан алгоритм проверки существования этих функций для рассматриваемых областей $D$ конкретного вида.

Запишем граничные условия (5), (7) для функции вида (8):

Условия (9) – граничные условия для гармонических функций $\varphi $ и $\psi $.

Из формулы Гаусса–Остроградского (см., например, [36]) следует, что гармонические в плоской области $D$ функции $\varphi $ и $\psi $ удовлетворяют условиям:

где направление обхода контуров выбирается так, чтобы область $D$ оставалась слева.

Уравнения (10) могут быть использованы как условия для вычисления неопределенных постоянных ${{C}_{0}}$, ${{C}_{1}}$, входящих в граничные условия (7).

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

Для некоторой функции $\varphi (s)$, гармонической в односвязной области $D$, ограниченной замкнутым гладким контуром $\partial D$, справедливо тождество:

где $A$ и $B$ – линейные операторы:

Интегральное тождество для гармонической функции в таком виде, насколько нам известно, впервые было использовано в [22]. Его отличительной особенностью является то, что подынтегральные функции оператора $B$ не содержат особенностей.

Тождество (11) распространяется и на двухсвязную область $D$, ограниченную контурами $\partial {{D}_{0}}$ и $\partial {{D}_{1}}$, изображенную на фиг. 1, если $\varphi (s)$ гармоническая и однозначная в области $D$ функция. Для вывода этого факта можно сделать разрез $L$, приведенный на фиг. 1 штриховой линией. Тогда область, ограниченная контуром $\partial {{D}_{0}}$, правым берегом разреза $L$, контуром $\partial {{D}_{1}}$ и противоположным берегом разреза $L$, становится односвязной. Сумма интегралов $G(s,s{\kern 1pt} '){{\varphi }_{n}}(s{\kern 1pt} ')$ и ${{G}_{n}}(s,s{\kern 1pt} ')\left( {\varphi (s{\kern 1pt} ') - \varphi (s)} \right)$ по берегам разреза $L$ равна нулю в силу противоположности нормалей по берегам и однозначности подынтегральных функций.

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

Таким образом, каждая из двух гармонических функций $\varphi $ и $\psi $, входящих в (8), удовлетворяет на гладкой границе $\partial D = \partial {{D}_{0}} \cup \partial {{D}_{1}}$ двухсвязной области $D$ интегральному тождеству

где $s$ меняется от 0 до ${{l}_{k}}$, $A$ и $B$ – линейные операторы: Тогда уравнения (9), (10) и (12) можно рассматривать как систему двух алгебраических и четырех интегральных уравнений для четырех граничных функций $\varphi (s)$, $\psi (s)$, ${{\varphi }_{n}}(s)$, ${{\psi }_{n}}(s)$ и двух констант ${{C}_{1}}$, ${{C}_{2}}$, входящих в (7). Алгоритм вычисления этих функций состоит из следующих процедур.

1. Составляется сетка на граничных контурах, состоящая из $N$ узловых точек.

2. Для каждой из четырех функций $\varphi (s)$, $\psi (s)$, ${{\varphi }_{n}}(s)$, ${{\psi }_{n}}(s)$ в узлах сетки определяются $N$значений $\varphi ({{s}_{i}})$, $\psi ({{s}_{i}})$, ${{\varphi }_{n}}({{s}_{i}})$, ${{\psi }_{n}}({{s}_{i}}),$ $i = 1,\;...,\;N$.

3. Аппроксимируя интегральные уравнения (12), получим две линейные системы 2N уравнений, а аппроксимация интегральных условий (10) даст еще 2 уравнения. Равенство соотношений (9) в узлах сетки даст еще 2N уравнений. Итого мы получим 4N + 2 уравнений относительно 4N неизвестных значений $\varphi ({{s}_{i}})$, $\psi ({{s}_{i}})$, ${{\varphi }_{n}}({{s}_{i}})$, ${{\psi }_{n}}({{s}_{i}}),$$i = 1,\;...,\;N$, и двух констант ${{C}_{1}}$, ${{C}_{2}}$. От точности аппроксимаций зависит точность решения задачи. Ниже будут применяться квадратурные формулы без насыщения, которые обеспечат ненасыщаемость схемы вычисления искомых функций с экспоненциальным уменьшением погрешности.

2.3. Квадратурные формулы и дискретизация системы уравнений

Ранее границы области $\partial {{D}_{k}}$ были заданы параметрически уравнениями $x = {{x}_{k}}(s)$, $y = {{y}_{k}}(s)$ ($k = 0,1$), имеющими естественный период ${{l}_{k}}$, где ${{l}_{k}}$ – длина дуги контура $\partial {{D}_{k}}$. Далее удобно ввести новый параметр $\zeta $ и задавать контуры области параметрическими уравнениями:

где ${{x}_{k}}(\zeta )$, ${{y}_{k}}(\zeta )$ ($k = 0,1$) – функции параметра $\zeta $ с единичным периодом.

Вводится дискретизация контуров конечным числом точек $M_{i}^{{(0)}}$ ($i = \overline {1,{{N}_{0}}} $), $M_{i}^{{(1)}}$ ($i = \overline {1,{{N}_{1}}} $) так, что точке $M_{i}^{{(k)}}$ соответствует значение $\zeta _{i}^{{(k)}} = \frac{i}{{{{N}_{k}}}}$. Параметр $\zeta $ и координата $s$ на каждом контуре связаны соотношением $ds = {{l}_{k}}{{f}_{k}}(\zeta )d\zeta $ ($k = 0,1$), ${{f}_{k}}(\zeta )$ – плотность точек на контуре $\partial {{D}_{k}}$. В частности, если ${{f}_{k}}(\zeta ) \equiv 1$, то точки на соответствующем контуре распределены равномерно.

Значения функций $\varphi (s)$, $\psi (s)$, $x(s)$, ${{\varphi }_{n}}(s)$, ${{\psi }_{n}}(s)$, ${{x}_{n}}(s) = \partial x{\text{/}}\partial n = \partial y{\text{/}}\partial s$ в точках $M_{i}^{{(0)}}$ ($i = \overline {1,{{N}_{0}}} $), $M_{i}^{{(1)}}$ ($i = \overline {1,{{N}_{1}}} $) обозначим через ${{\Phi }_{q}}$, ${{\Psi }_{q}}$, ${{X}_{q}}$, $\Phi _{q}^{n}$, $\Psi _{q}^{n}$, $X_{q}^{n}$, $q = \overline {1,{{N}_{0}} + {{N}_{1}}} $. Первые ${{N}_{0}}$ индексов относятся к точкам $M_{i}^{{(0)}}$ на контуре $\partial {{D}_{0}}$, а вторые ${{N}_{1}}$ индексов – к точкам $M_{i}^{{(1)}}$на контуре $\partial {{D}_{1}}$.

Запишем алгебраическую систему уравнений (9) после дискретизации

Для аппроксимации интегралов используем квадратурные формулы В статьях [24], [27] доказано, что эти формулы для аналитических функций $f(x)$ периода $l$ являются квадратурными формулами без насыщения и имеют порядок погрешности выше любой степени по шагу сетки $h$. Формула (15) применяется для оператора $B$ и интегралов (10), подынтегральные функции которых не содержат особенностей. Формула (16) применяется для оператора $A$, который содержит логарифмическую особенность.

Система интегральных уравнений (10) с помощью квадратурной формулы без насыщения (15) запишется в виде

где

Наконец, запишем систему уравнений, аппроксимирующую уравнения (12):

где элементы блочных матриц ${\mathbf{A}}$ и ${\mathbf{B}}$ выражаются через значения функции Грина и ее производной по квадратурным формулам без насыщения (15) и (16).

Таким образом, имеем $4({{N}_{0}} + {{N}_{1}})$ значений ${{\Phi }_{q}},\;{{\Psi }_{q}},\;\Phi _{q}^{n},\;\Psi _{q}^{n},$ $q = \overline {1,{{N}_{0}} + {{N}_{1}}} $, и две константы ${{C}_{0}}$, ${{C}_{1}}$. Для их определения получена система линейных уравнений из $4({{N}_{0}} + {{N}_{1}}) + 2$ уравнений. В нее входят $2({{N}_{0}} + {{N}_{1}})$ уравнений (14), 2 уравнения (17) и $2({{N}_{0}} + {{N}_{1}})$ уравнений (18).

Окончательно систему линейных алгебраических уравнений (14), (17), (18) можно записать в самом общем матричном виде:

где ${\mathbf{Z}}$ – столбец неизвестных значений функций $\varphi (s),\;\psi (s),\;{{\varphi }_{n}}(s),\;{{\psi }_{n}}(s)$ и констант ${{C}_{0}}$, ${{C}_{1}}$, ${\mathbf{g}}$ – столбец правых частей уравнений, а матрица ${\mathbf{P}}$ определяется только формой области $D$ и не зависит от вида граничных условий задачи.

2.4. Определение функции тока и скорости жидкости внутри области

Решив систему уравнений (19), определим неизвестные дискретные значения функций $\varphi (s),\;\psi (s),\;{{\varphi }_{n}}(s),\;{{\psi }_{n}}(s)$ на границах $\partial {{D}_{0}}$ и $\partial {{D}_{1}}$ области $D$, зная которые можно приближенно определить значение искомой бигармонической функции тока $\Psi $ в любой точке, лежащей внутри области $D$.

Для гармонических функций $\varphi $ и $\psi $ во внутренней точке области $D$ справедливы интегральные выражения:

Интегралы (20) не имеют особенностей, и для функций $\varphi $ и $\psi $ внутри области $D$ с помощью квадратурной формулы без насыщения (15) получаются следующие приближенные выражения: где ${{N}_{{0,1}}} = {{N}_{0}}$ при $q = 1,2,\;...,\;{{N}_{0}}$ и ${{N}_{{0,1}}} = {{N}_{1}}$ при $q > {{N}_{0}}$.

Тогда окончательное выражение для искомой бигармонической функции тока $\Psi $ внутри области $D$ примет вид

Для полного решения необходимо также определить поле скоростей вязкой жидкости внутри области $D$ по формулам (2), т.е. следует найти частные производные найденной функции тока (22). Поскольку эти функции заданы явно как функции переменных $x$ и $y$, то дифференцирование можно провести аналитически. Тогда для компонент скорости произвольной внутренней точки $(x,y)$ двусвязной области $D$ получим

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

3.1. Проверка существования представления Гурса и сходимости расчетной схемы

Пусть бигармоническая в двухсвязной области $D$ функция $\Psi $ имеет представление Гурса $\Psi = \varphi (x,y) + x \cdot \psi (x,y)$ через 2 гармонические функции $\varphi (x,y)$ и $\psi (x,y)$. Тогда, как показано выше, ${\mathbf{Z}}$ – столбец размером ${{N}_{0}} + {{N}_{1}} + 2$, в который входят неизвестные значения функций $\varphi (s),\;\psi (s),\;{{\varphi }_{n}}(s),\;{{\psi }_{n}}(s)$ на контурах $\partial {{D}_{0}}$ и $\partial {{D}_{1}}$ и две константы ${{C}_{0}}$, ${{C}_{1}}$, находятся из системы линейных алгебраических уравнений (19)

где ${\mathbf{g}}$ – заданный столбец правых частей уравнений размером ${{N}_{0}} + {{N}_{1}} + 2$, а матрица ${\mathbf{P}}$, размером $({{N}_{0}} + {{N}_{1}} + 2) \times ({{N}_{0}} + {{N}_{1}} + 2)$, определяется только формой области $D$ и не зависит от вида граничных условий задачи. Столбец Z можно найти, обращая матрицу ${\mathbf{P}}$ в виде ${\mathbf{Z}} = {{{\mathbf{P}}}^{{ - {\mathbf{1}}}}}{\mathbf{g}}$, где обратная матрица ${{{\mathbf{P}}}^{{ - {\mathbf{1}}}}}$ тоже не зависит от вида граничных условий задачи. Если существует обратная матрица ${{{\mathbf{P}}}^{{ - {\mathbf{1}}}}}$, то это означает существование столбца ${\mathbf{Z}}$ и значений гармонических функций $\varphi (s),\;\psi (s),\;{{\varphi }_{n}}(s),\;{{\psi }_{n}}(s)$ на контурах $\partial {{D}_{0}}$ и $\partial {{D}_{1}}$. Следовательно, для существования представления Гурса нужно доказать существование обратной матрицы ${{{\mathbf{P}}}^{{ - {\mathbf{1}}}}}$, что можно выполнить вычислительным путем.

Проверить сходимость схемы без насыщения для любой двухсвязной области $D$ можно с помощью построения краевой задачи (6) с точным аналитическим решением. Для этого можно выбрать бигармоническую функцию $u(x,y)$, например в виде кубического полинома от переменных $x,y$, вычислить ее граничные значения ${{\left. u \right|}_{{\partial {{D}_{k}}}}} = {{g}^{{(k)}}}(s)$ и ее нормальное производное ${{\left. {\frac{{\partial u}}{{\partial n}}} \right|}_{{\partial {{D}_{k}}}}} = g_{n}^{k}(s)$. Функция $\Psi (x,y)$ будет очевидно являться точным решением для этой задачи. После этого вычисляем по предложенной выше численной схеме граничные значения ${{\Phi }_{q}},\;{{\Psi }_{q}},\;\Phi _{q}^{n},\;\Psi _{q}^{n},\;{{X}_{q}},$ $q = \overline {1,{{N}_{0}} + {{N}_{1}}} $ и проверяем сходимость численных значений ${{\Phi }_{q}} + {{X}_{q}}\,{{\Psi }_{q}}$ к их точным граничным значениям $u({{X}_{q}},{{Y}_{q}})$. В любой внутренней точке $(x,y) \in D$ приближенное значение решения вычисляется по формуле (22). Его можно сравнить с точным значением $u(x,y)$ и найти погрешность.

Пусть рассматривается тестовая функция $u\left( {x,y} \right) = {{x}^{3}} + {{y}^{3}}$ – бигармоническая в двусвязной области, ограниченной окружностями радиусов ${{\rho }_{1}} = 1$ и ${{\rho }_{0}} = 0.3$, расстояние между центрами которых $\Delta \,x = 0.35$ (фиг. 2).

Результаты численного решения задачи при ${{N}_{0}} = 32$, ${{N}_{1}} = 64$ представлены на фиг. 3: серыми сплошными линиями изображены графики значений тестовой функции ${{u}_{1}}(\zeta )$, ${{u}_{2}}(\zeta )$, ${{u}_{3}}(\zeta )$, $\zeta \in (0,1)$ на трех окружностях 1, 2 и 3, показанных на фиг. 2; черными штриховыми линиями – соответствующие графики значений приближенной функции ${{\tilde {u}}_{1}}(\zeta )$, ${{\tilde {u}}_{2}}(\zeta )$, ${{\tilde {u}}_{3}}(\zeta )$, найденной численно с помощью формул (21) и (22). В табл. 1 приведены значения относительной погрешности функции $u$ на этих окружностях для различных значений ${{N}_{0}}$$\left( {{{N}_{1}} = 2{{N}_{0}}} \right)$. Указанная погрешность вычисляется по формуле

где ${{u}_{i}}$ – точное значение исследуемой величины (соответствует аналитическому решению), а ${{\tilde {u}}_{i}}$ – приближенное ее значение, полученное численно.

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

Таким образом, вопрос о существовании представления Гурса (8) может быть решен для каждой конкретной области с помощью составления матрицы ${\mathbf{P}}$. Если матрица ${\mathbf{P}}$ имеет обратную, то для любой правой части ${\mathbf{g}}$ системы (19) в данной области могут быть определены гармонические функции $\varphi $ и $\psi $. Это будет доказывать и существование представления Гурса для данной области и экспоненциальную сходимость указанного метода. В частности, для всех рассмотренных в данной статье областей указанное условие выполняется.

3.2. Движение вязкой жидкости между двумя вращающимися концентрическими круговыми цилиндрами

Рассматривается задача о движении вязкой жидкости между двумя вращающимися концентрично расположенными круговыми цилиндрами. Пусть ${{\rho }_{0}}$ – радиус внутреннего цилиндра (вращается с угловой скоростью ${{\omega }_{0}}$), ${{\rho }_{1}}$ – радиус ${{\omega }_{0}}$ внешнего цилиндра (вращается с угловой скоростью ${{\omega }_{1}}$). Задача рассматривается в полярных координатах ($\rho ,\theta $). Граничные условия:

где $V$ – алгебраическое значение линейной скорости точек жидкости.

Решение такой задачи описано, например, в [22], где показано, что скорость $V$ зависит только от $\rho $ и имеет вид

где следующие коэффициенты:

Для численного решения задачи от граничных условий (24) необходимо перейти к граничным условиям для функции тока. В случае такой постановки, если оси вращения цилиндров проходят через начало координат, то граничные условия (5), (7) для функции тока, с учетом неопределенности констант ${{C}_{k}}$, запишутся в виде

Результаты численного решения задачи для некоторых различных значений входных параметров представлены на фиг. 4: красной сплошной линией изображены графики точного значения модуля скорости в зависимости от $\rho $, а синей штриховой линией – графики соответствующего модуля скорости, найденного численно при ${{N}_{0}} = 200$ и ${{N}_{1}} = 400$. В табл. 2 для каждого случая, представленного на графиках фиг. 2, для различных ${{N}_{0}}$ (${{N}_{1}} = 2{{N}_{0}}$) приведены значения относительной погрешности модуля скорости, которая вычисляется по формуле, аналогичной формуле (23).

Следует также заметить, что на серединной окружности рассматриваемых областей $\left( {\rho = \frac{{{{\rho }_{0}} + {{\rho }_{1}}}}{2}} \right)$ относительная погрешность практически нулевая (в зависимости от параметров области и числа граничных элементов она принимает значения порядка ${{10}^{{ - 12}}}{\kern 1pt} - {\kern 1pt} {{10}^{{ - 16}}}$). Из-за того, что потенциал двойного слоя имеет скачок на границе, то при приближении к ней на расстояние, сравнимое с расстоянием между выбранными точками границы, наблюдается увеличение погрешности. С увеличением количества точек на границе область расхождения уменьшается. Для снижения погрешности в “узких” участках области удобно применять сгущение точек, что достигается за счет функции ${{f}_{k}}\left( \zeta \right)$ ($k = 0,\;1$). Этот прием проиллюстрирован ниже на тестовой задаче.

3.3. Движение вязкой жидкости между двумя вращающимися эксцентрическими круговыми цилиндрами

Рассмотрим теперь течение вязкой жидкости между двумя вращающимися эксцентрично расположенными круговыми цилиндрами. Пусть ${{\rho }_{0}}$ – радиус внутреннего цилиндра (вращается с угловой скоростью ${{\omega }_{0}}$), ${{\rho }_{1}}$ – радиус внешнего цилиндра (вращается с угловой скоростью ${{\omega }_{1}}$), расстояние между центрами цилиндров $\Delta {\kern 1pt} x$. Положительным является направление вращения против часовой стрелки.

По аналогии с результатами, полученными в [1] для случая вращения только внутреннего цилиндра, в [5] получено аналитическое решение рассматриваемой задачи.

Численное решение задачи проведено по описанному в данной статье алгоритму. При этом применено неравномерное распределение точек на контурах, что позволяет улучшить точность без увеличения количества точек на границах. В качестве функций ${{f}_{k}}(\zeta )$, которые характеризуют плотность распределения точек на контуре $\partial {{D}_{k}}$ ($k = 0,1$), используются функции

В частности, при ${{\rho }_{1}} = 1$, ${{\rho }_{0}} = 4$ для различных значений $\Delta x$ получаются следующие картины распределения точек на контурах (${{N}_{0}} = 20$, ${{N}_{1}} = 48$), представленные на фиг. 5.

Если ось вращения внешнего цилиндра проходит через начало координат, граничные условия (5), (7) для функции $\Psi $ примут вид

На фиг. 6, 7 показано сравнение результатов численного и аналитического решения для различных значений параметров ${{\rho }_{0}}$, ${{\rho }_{1}}$, ${{\omega }_{0}}$, ${{\omega }_{1}}$, $\Delta x$ (${{N}_{0}} = 40$, ${{N}_{1}} = 100$).

При наложении фигур 6а и 6б и фигур 7а и 7б линии тока совпадают, что говорит о достаточной точности численного решения.

Для более точной количественной оценки погрешности проведем детальный анализ для еще одного случая: ${{\rho }_{0}} = 0.3$, ${{\rho }_{1}} = 1$, ${{\omega }_{0}} = - 8$, ${{\omega }_{1}} = 4$, $\Delta {\kern 1pt} x = 0.35$ (параметры области такие же, как в п. 3.1). На фиг. 8 представлены результаты численного решения задачи для указанных значений входных параметров при ${{N}_{0}} = 60$, ${{N}_{1}} = 120$: на фиг. 8а приведены линии тока $\Psi = {\text{const}}$ и для каждой линии указано числовое значение функции тока, которому она соответствует; на фиг. 8б представлено векторное поле скоростей точек жидкости (длина вектора соответствует его величине).

Для оценки погрешности полученного решения проведем сначала сравнение функции тока, полученной численно и аналитически. Из точного решения может быть найдена точная функция тока $\Psi (x)$ в зависимости от координаты $x$ при $y = 0$. Соответствующая функция $\tilde {\Psi }(x)$ находится из приближенного решения с помощью (22). В табл. 3 приведены значения относительной погрешности $\delta \left( {\Psi (x)} \right) = \frac{{\max \left| {\Psi (x) - \tilde {\Psi }(x)} \right|}}{{\max \left| {\Psi (x)} \right|}}$ для различных значений ${{N}_{0}}$ (${{N}_{1}} = 2{{N}_{0}}$) на интервале $x \in \left[ { - 0.95, - 0.7} \right] \cup \left[ {0,0.95} \right]$.

Далее выясним вопрос о распределении погрешности модуля скорости внутри области течения. На фиг. 9 приведены графики относительной погрешности модуля скорости на внутренней окружности 2, изображенной на фиг. 2:

где ${{v}_{2}}(\zeta )$ – точное значение модуля скорости в точке $\zeta $, а ${{\tilde {v}}_{2}}(\zeta )$ – приближенное значение, полученное численно. Окружность симметрична относительно оси Ox, поэтому графики приводим для верхней полуокружности (параметр $\zeta \in \left[ {0,0.5} \right]$ по оси абсцисс) для различных значений количества граничных элементов ${{N}_{0}}\;\left( {{{N}_{1}} = 2{{N}_{0}}} \right)$.

Представленные на фиг. 9 графики показывают, как зависит относительная погрешность модуля скорости на дуге одной из окружностей внутри расчетной области от параметра $\zeta $ при различных значениях числа узлов сетки ${{N}_{0}}$. Из графиков видно, что погрешности малы и почти хаотичны в расчетной области. Поэтому более информативной представляется табл. 4, в которой приведены относительные погрешности модуля скорости, вычисленные по формуле (23), вдоль окружностей внутри области, изображенных на фиг. 2, для различных значений ${{N}_{0}}$ (${{N}_{1}} = 2{{N}_{0}}$).

Из проведенного в этом разделе сравнения численного и аналитического решений видно, что численные результаты в расчетной области мало отличаются уже при небольших значениях ${{N}_{0}}$ и ${{N}_{1}}$. Точность может нарушаться лишь вблизи границ, однако с увеличением значений ${{N}_{0}}$ и ${{N}_{1}}$ область расхождения уменьшается. Кроме того, улучшение точности при одном и том же числе граничных элементов может быть достигнуто за счет неравномерного распределения точек на контурах. При увеличении числа элементов сетки по арифметической прогрессии погрешность убывает по геометрической прогрессии. Таким образом, можно сделать вывод об эффективности предложенного алгоритма для расчета течения вязкой жидкости, заключенной между двумя произвольно движущимися цилиндрами произвольного сечения.

4.1. Расчет модели миксера

Рассмотрим течение вязкой жидкости, заключенной между вращающейся с угловой скоростью $\omega $ прямоугольной лопаткой с линейными размерами ${{a}_{0}}$ и ${{b}_{0}}$ (внутренний контур) и круговым цилиндром радиуса ${{\rho }_{1}}$ (внешний контур), расстояние между центрами круга и прямоугольника равно $\Delta {\kern 1pt} x$ (модель миксера). Пусть выбрана система координат, в которой внутренний контур неподвижен, тогда, очевидно, внешний круговой контур вращается в ней с угловой скоростью $ - \omega $. В такой системе координат течение будет стационарно. Пусть начало координат расположено в центре внешнего кругового контура, стороны прямоугольника параллельны осям координат, центр прямоугольника лежит на оси $Ox$.

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

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

На фиг. 10 представлены линии тока $\Psi = {\text{const}}$, полученные в результате численного расчета для двусвязной области описанного выше вида при $N = 20$ и при количестве граничных элементов ${{N}_{0}} = 80$, ${{N}_{1}} = 200$ для следующих значений входных параметров: ${{\rho }_{1}} = 10$, ${{a}_{0}} = 1$, ${{b}_{0}} = 5$, $\hat {x} = 3$, $\omega = 1$. Формулы (25) соответствуют вертикальному положению лопатки с центром в начале координат, преобразования координат, необходимые для получения контуров, изображенных на фиг. 6, очевидны.

На фиг. 11 изображены линии тока, полученные в результате численного расчета в случае лопатки более сложной формы: параметрические уравнения рассматриваемого контура совпадают с уравнениями (25), если положить в них

Расчеты проведены при ${{N}_{0}} = 120$, ${{N}_{1}} = 180$ (фиг. 10а: $\omega = 4$, фиг. 10б: $\omega = 1$).

4.2. Вращение кругового цилиндра внутри некругового

Приведем результаты расчета течения вязкой жидкости, заключенной между вращающимся с угловой скоростью $\omega $ круговым цилиндром радиуса ${{\rho }_{0}}$ (внутренний контур) и внешним цилиндром эллиптического и сглаженного прямоугольного (описывается формулами (23)) сечения. Расстояние между центрами внешнего и внутреннего контуров равно $\Delta {\kern 1pt} x$. Течение в этом случае также стационарно, поэтому течение жидкости может быть описано с помощью линий тока. Пусть начало координат расположено в центре внешнего контура, центр внутреннего контура лежит на оси $Ox$.

На фиг. 12 изображены линии тока $\Psi = {\text{const}}$, полученные в результате численного расчета при ${{N}_{0}} = 60$, ${{N}_{1}} = 180$ для $\omega = 1$ для внешнего контура в форме: (а) – эллипса с полуосями ${{a}_{1}} = 1$ и ${{b}_{1}} = 0.6$ при ${{\rho }_{0}} = 0.2$, $\Delta {\kern 1pt} x = 0.3$; (б) – скругленного прямоугольника со сторонами ${{a}_{1}} = 1$ и ${{b}_{1}} = 0.6$ при ${{\rho }_{0}} = 0.15$, $\Delta {\kern 1pt} x = 0.3$; (в) – эллипса с полуосями ${{a}_{1}} = 1$ и ${{b}_{1}} = 0.4$ при ${{\rho }_{0}} = 0.1$, $\hat {x} = 0.04$. Значения функции тока на приведенных линиях уровня представлены в табл. 5. Эти значения несут информацию о величине скорости на линиях тока: она, как известно, пропорциональна изменению значений функции тока.

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