Exponentially convergent numerical scheme for the stream function of potential flow over axisymmetric bodies

Abstract

A boundary element scheme for the problem of potential flow over an axisymmetric toroidal body is considered. An integral equation for the velocity distribution on the body is derived. It is shown that the numerical scheme for solving the considered equation converges exponentially.

Full Text

Численная схема с экспоненциальной сходимостью для функции тока потенциального обтекания тел с осевой симметрией [1]

ВВЕДЕНИЕ

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

В работах [1]–[5] были разработаны численные схемы метода граничных элементов для решения задач потенциального течения идеальной несжимаемой жидкости. Потенциал поля скорости удовлетворяет уравнению Лапласа, и для него выводится линейное интегральное уравнение на граничной поверхности. Оно связывает между собой значение функции и ее нормальную производную на этой поверхности. Для плоской и осесимметричной задач граничную поверхность определяет одномерный контур. Поэтому интегральное уравнение является одномерным, и с помощью удачно подобранных квадратурных формул интегральное уравнение аппроксимируется линейной системой уравнений. Для задачи Дирихле по заданному потенциалу из системы уравнений находится нормальная производная. Для задачи Неймана по заданной нормальной производной потенциала находится потенциал. Потенциал в каждой точке области, ограниченной контуром, линейно выражается через распределение потенциала и его нормальной производной на контуре. В этом и состоит суть метода граничных элементов.

Существует также много более современных работ [6]$–$[10], в которых предлагаются различные численные схемы также для потенциала.

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

Интегральные операторы интегрального уравнения на граничном контуре действуют на периодические функции. Периодом является длина контура. Если контур аналитичен, то коэффициенты n-й гармоники ряда Фурье периодических функций убывают по экспоненте $\exp \mathrm{<mo\; stretchy="false">(</mo>}-cn\mathrm{<mo\; stretchy="false">)</mo>}$. С помощью этого наблюдения в [12], [13] разработаны аппроксимации интегральных уравнений линейной системой, погрешность которых убывает с ростом числа элементов сетки по экспоненте.

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

1. ФУНКЦИЯ ГРИНА

Компоненты скорости $v_r\mathrm{,}$ $v_z$ в цилиндрической системе координат $z\mathrm{,}r\mathrm{,}\phi$ выражаются через функцию тока $\Psi$ следующим образом:

$\begin{array}{c}{v}_r\mathrm{<mo>=</mo>}\frac{1}{r}\frac{\partial \Psi }{\partial z}\mathrm{,}{v}_z\mathrm{<mo>=</mo>}-\frac{1}{r}\frac{\partial \Psi }{\partial r}.\end{array}$

Угловая компонента вектора вихря $\text{r}ot\nu$ имеет вид

$\begin{array}{c}\text{r}o{t}_{\phi }\nu \mathrm{<mo>=</mo>}\frac{\partial v_r}{\partial z}-\frac{\partial v_z}{\partial r}\mathrm{<mo>=</mo>}\frac{\partial }{\partial z}\left(\frac{1}{r}\frac{\partial \Psi }{\partial z}\right)+\frac{\partial }{\partial r}\left(\frac{1}{r}\frac{\partial \Psi }{\partial r}\right)\mathrm{<mo>=</mo>}{D}^2\Psi \mathrm{.}\end{array}$

Рассмотрим потенциальное осесимметричное течение вне осесимметричного тора, $S$ $—$ область, которая получается в результате сечения тора меридианальной плоскостью. Из приведенных равенств получаем уравнение потенциального (безвихревого) течения

$D^2\Psi \mathrm{<mo>=</mo>0.}$

Функцией Грина $G$ осесимметричного потенциального течения будет функция тока вихревого осесимметричного кольца $C$ радиуса $a$. Согласно формуле Био–Савара поле скорости вихревого кольца имеет вид (фиг. 1):

$v\mathrm{<mo>=</mo>}-{\displaystyle \oint }\frac{d{s}^{\prime }\times \mathrm{<mo\; stretchy="false">(</mo>}{r}^{\prime }-r\mathrm{<mo\; stretchy="false">)</mo>}}{\mathrm{<mo>|</mo>}{r}^{\prime }-r{\mathrm{<mo>|</mo>}}^3}\mathrm{.}$

Фиг. 1.Вихревое кольцо

В декартовой системе координат $x\mathrm{,}$ $y\mathrm{,}$ $z$ векторы под интегралом имеют компоненты

$r^{\prime }\mathrm{<mo\; stretchy="false">(</mo>}a\cos \varphi \mathrm{,}a\sin \varphi \mathrm{,0<mo\; stretchy="false">)</mo>,}r\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,0,}z\mathrm{<mo\; stretchy="false">)</mo>,}d{s}^{\prime }\mathrm{<mo\; stretchy="false">(</mo>}-a\sin \varphi \mathrm{,}a\cos \varphi \mathrm{,0<mo\; stretchy="false">)</mo>}d\varphi \mathrm{.}$

Отсюда находим радиальную скорость $v_r$ вихревого кольца

$v_r\mathrm{<mo>=</mo>}{\displaystyle {\int }_0^{2\pi }}\frac{az\cos \varphi d\varphi }{{\mathrm{<mo\; stretchy="false">(</mo>}{a}^2-2ar\cos \varphi +r^2+z^2\mathrm{<mo\; stretchy="false">)</mo>}}^{\mathrm{3<mo>/</mo>2}}}\mathrm{.}$

Ее можно выразить через функцию тока $G$:

$\begin{array}{c}\begin{array}{c}{v}_r\mathrm{<mo>=</mo>}\frac{1}{r}\frac{\partial G}{\partial z}\mathrm{,}\\ G\left(r\mathrm{,}z\mathrm{,}a\right)\mathrm{<mo>=</mo>}-{\displaystyle {\int }_0^{2\pi }}\frac{ar\cos \varphi d\varphi }{{\left(a^2-2ar\cos \varphi +r^2+z^2\right)}^{\mathrm{1<mo>/</mo>2}}}\mathrm{.}\end{array}\end{array}$ (1)

Перейдем теперь к выводу функции Грина. Пусть $M\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,}z\mathrm{<mo\; stretchy="false">)</mo>}$ и $M\mathrm{<mo\; stretchy="false">(</mo>}{r}^{\prime }\mathrm{,}{z}^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ две произвольные точки пространства. Функция Грина $G\mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{,}{M}^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}$ осесимметричного течения строится следующим образом. Проводим контур через точку $M\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,}z\mathrm{<mo\; stretchy="false">)</mo>}$ и вычисляем значение функции тока вихревого кольца в точке $M^{\prime }\mathrm{<mo\; stretchy="false">(</mo>}{r}^{\prime }\mathrm{,}{z}^{\prime }\mathrm{<mo\; stretchy="false">)</mo>.}$ Найденное значение будет равно значению функции Грина в точках $M$ и $M^{\prime }\mathrm{.}$ По построению значение $G\mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{,}{M}^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}$ будет равно функции тока (1), в которой первый аргумент равен $r^{\prime }\mathrm{,}$ второй $–$ ( $z^{\prime }-z$ ) и третий $a\mathrm{<mo>=</mo>}r$. Таким образом, функция Грина такова:

$G\left(M^{\prime }\mathrm{,}M\right)\mathrm{<mo>=</mo>}-{\displaystyle {\int }_0^{2\pi }}\frac{r^{\prime }r\cos \varphi d\varphi }{{\left[{\left(r^{\prime }\right)}^2-2r^{\prime }r\cos \varphi +r^2+{\left(z-z^{\prime }\right)}^2\right]}^{\mathrm{1<mo>/</mo>2}}}\mathrm{.}$

Функция Грина симметрична относительно перестановки точек $M$ и $M^{\prime }\mathrm{.}$ Поэтому для определения $G\mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{,}{M}^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}$ контур можно также провести через точку $M^{\prime }$ и вычислить функцию тока в точке $M\mathrm{.}$ Полученный интеграл $G\mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{,}{M}^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}$ можно выразить через эллиптические интегралы (см. [14, т. 1, с. 202])

$G\left(M\mathrm{,}{M}^{\prime }\right)\mathrm{<mo>=</mo>}-2\left(r_2+r_1\right)\left[K\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}-E\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}\right]\mathrm{,}$ (2)

где приняты следующие обозначения:

$\begin{array}{c}{r}_1\mathrm{<mo>=</mo>}\sqrt{{\left(z^{\prime }-z\right)}^2+{\left(r^{\prime }-r\right)}^2}\mathrm{,}{r}_2\mathrm{<mo>=</mo>}\sqrt{{\left(z^{\prime }-z\right)}^2+{\left(r^{\prime }+r\right)}^2}\mathrm{,}k\mathrm{<mo>=</mo>}\frac{r_2-r_1}{r_2+r_1}\mathrm{,}\\ K\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle {\int }_0^{\pi \mathrm{<mo>/</mo>2}}}{\left(1-k^2{\sin }^2{x}^2\right)}^{-\mathrm{1<mo>/</mo>2}}dx\mathrm{,}E\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle {\int }_0^{\pi \mathrm{<mo>/</mo>2}}}{\left(1-k^2{\sin }^2{x}^2\right)}^{\mathrm{1<mo>/</mo>2}}dx\mathrm{.}\end{array}$

Для функции $K\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}-E\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}$ полезна аппроксимация

$\begin{array}{c}K\left(k\right)-E\left(k\right)\mathrm{<mo>=</mo>}P\left(m_1\right)+\ln \left(\mathrm{1<mo>/</mo>}{m}_1\right)Q\left(m_1\right)\mathrm{,}\\ P\left(m_1\right)\approx 0.3862943611-0.3468568385m_1-\\ -0.0305489156m_1^2-0.0088886006m_1^3\mathrm{,}\\ Q\left(m_1\right)\approx 0.5-0.1249477647m_1-\\ -0.02138998364m_1^2-0.00288878875m_1^3\mathrm{,}\\ m_1\mathrm{<mo>=</mo>1}-k^2\mathrm{<mo>=</mo>}\frac{4r_1{r}_2}{{\left(r_1+r_2\right)}^2}\mathrm{.}\end{array}$

Максимальное отклонение аппроксимации от точного значения в области определения $\mathrm{0<mo><</mo>}k\mathrm{<mo><</mo>1}$ не превышает величину $1.33\times {10}^{-7}$.

2. ИНТЕГРАЛЬНОЕ УРАВНЕНИЕ ДЛЯ ФУНКЦИИ ТОКА

Теория интегральных уравнений для плоских задач уравнения Лапласа переносится на осесимметричный случай следующим образом. Оператору Лапласа ставится в соответствие оператор $D^2\text{}\mathrm{,}$ а гармонической функции $\Phi$ — функция тока $\Psi \mathrm{.}$ Кроме того, в интегралах нужно учесть весовой множитель $\mathrm{1<mo>/</mo>}{r}^{\prime }\mathrm{,}$ где $r^{\prime }$ $—$ координата точки интегрирования $M^{\prime }\mathrm{.}$

Для функции $\Psi \mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{<mo\; stretchy="false">)</mo>,}$ удовлетворяющей в области $M\in S$, ограниченной замкнутым гладким контуром $\partial S$, введем следующий интеграл:

$\begin{array}{c}I\left(M\right)\mathrm{<mo>=</mo>}\tilde{A}V\left(M\right)+\tilde{B}\frac{\Psi }{r}\left(M\right)\mathrm{,}V\left(M\right)\mathrm{<mo>=</mo>}\frac{1}{r}\frac{\partial \Psi }{\partial n}\left(M\right)\mathrm{,}\end{array}$ (3)

где $\tilde{A}$ и $\tilde{B}$ $—$ два интегральных оператора:

$\begin{array}{c}\tilde{A}V\left(M\right)\mathrm{<mo>=</mo>}-{\displaystyle \underset{0}{\overset{l}{\int }}}G\left(M\mathrm{,}{M}^{\prime }\right)V\left(M^{\prime }\right)d{s}^{\prime }\mathrm{,}\\ \tilde{B}\Phi \left(M\right)\mathrm{<mo>=</mo>}{\displaystyle \underset{0}{\overset{l}{\int }}}\frac{\partial G\left(M\mathrm{,}{M}^{\prime }\right)}{\partial n^{\prime }}\frac{\Psi }{r}\left(M^{\prime }\right)d{s}^{\prime }\mathrm{,}\end{array}$ (4)

Интегралы берутся по контуру $\partial S$. Контур $\partial S$ разбивает плоскость на две области: внутреннюю $S\mathrm{,}$ ограниченную контуром $\partial S$ снаружи, и внешнюю $\overline{S}\mathrm{,}$ ограниченную контуром изнутри, $M$ $–$ произвольная точка на плоскости, $M^{\prime }\in \partial S$ $–$ точка интегрирования, $d{s}^{\prime }$ $–$ элемент дуги в точке $M^{\prime }$, $l$ $–$ длина контура $\partial S$, $n$ $–$ вектор нормали к контуру $\partial S\mathrm{,}$ внешней по отношению к $S$, $\frac{\partial }{\partial n}$ и $\frac{\partial }{\partial n^{\prime }}$ $—$ производные по направлению внешней нормали в точках $M$ и $M^{\prime }$ соответственно, $G\mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{,}{M}^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ функция Грина, $r\mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{,}{M}^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ расстояние между точками $M$ и $M^{\prime }$. (Обозначения показаны на фиг. 2).

Выражения $\tilde{A}V\mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{<mo\; stretchy="false">)</mo>}$ и $\tilde{B}\Phi \mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{<mo\; stretchy="false">)</mo>}$ $—$ называются потенциалами простого и двойного слоев соответственно. Известно [15, с. 241$–$242], что потенциал простого слоя непрерывен на границе $\partial S\mathrm{,}$ а потенциал двойного слоя имеет скачок. Значение потенциала двойного слоя в граничной точке $M_0\in \partial S$ равно среднеарифметическому двух предельных значений потенциала при $M\to M_0\mathrm{,}$ точек $M$, находящихся внутри области $S$ ( $M\in S$ ) и снаружи – ( $M\in \overline{S}$ ).

Фиг. 2.Обозначения к основному тождеству

Пользуясь свойствами потенциала двойного слоя, можно вывести следующие тождества для интеграла (3):

$I\left(M\right)\mathrm{<mo>=</mo>}\left\{\begin{array}{ll}4\pi \Psi \left(M\right)\mathrm{,}\hfill & M\in S\mathrm{,}\hfill \\ 2\pi \Psi \left(M\right)\mathrm{,}\hfill & M\in \partial S\mathrm{,}\hfill \\ \mathrm{0,}\hfill & M\in S.\hfill \end{array}\right.$ (5)

Отсюда для значений функции тока на границе получим

$\begin{array}{c}2\pi \Psi \left(M\right)\mathrm{<mo>=</mo>}\tilde{A}V\left(M\right)+\tilde{B}\frac{\Psi }{r}\left(M\right)\mathrm{.}\end{array}$ (6)

Оно представляет линейную связь между значениями функции $\Psi$ и ее нормальной производной в граничных точках $M\in \partial S\mathrm{.}$ Так же как и в плоской задаче, тождество (7) можно рассматривать как интегральное уравнение относительно $\Psi \mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{<mo\; stretchy="false">)</mo>}$ при заданной нормальной производной или относительно $\partial \Psi \mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{<mo\; stretchy="false">)</mo><mo>/</mo>}\partial n$ при заданном значении $\Psi \mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{<mo\; stretchy="false">)</mo>}$ на границе. После того как $\partial \Psi \mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{<mo\; stretchy="false">)</mo><mo>/</mo>}\partial n$ и $\Psi \mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{<mo\; stretchy="false">)</mo>}$ будут найдены, можно вычислить функцию $\Psi \mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{<mo\; stretchy="false">)</mo>}$ во внутренних точках $M\in S$:

 $\begin{array}{c}4\pi \Psi \left(M\right)\mathrm{<mo>=</mo>}{\displaystyle \underset{\partial S}{\oint }}\frac{1}{r^{\prime }}\left[-G\left(M\mathrm{,}{M}^{\prime }\right)\frac{\partial \Psi }{\partial n}\left(M^{\prime }\right)+\Psi \left(M^{\prime }\right)\frac{\partial G\left(M\mathrm{,}{M}^{\prime }\right)}{\partial n}\right]d{s}^{\prime }\mathrm{.}\end{array}$ (7)

Особую точку интегрального уравнения (7) можно устранить так же, как в плоской задаче

$\begin{array}{c}AV\left(s\right)+B\frac{\Psi }{r}\left(s\right)\mathrm{<mo>=</mo>0,}\\ B\Phi \left(M\right)\mathrm{<mo>=</mo>}{\displaystyle \underset{0}{\overset{l}{\int }}}\frac{\partial G\left(M\mathrm{,}{M}^{\prime }\right)}{\partial n^{\prime }}\left(\frac{\Psi }{r}\left(M^{\prime }\right)-\frac{\Psi }{r}\left(M\right)\right)d{s}^{\prime }\end{array}$. (8)

Здесь оператор $\tilde{B}$ заменился на $B$, а оператор $A$ такой же, как $\tilde{A}$ в формуле (4).

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

3. УРАВНЕНИЯ ДЛЯ ФУНКЦИИ ТОКА ПРИ ОБТЕКАНИИ ТОРА

Пусть $\Psi \mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,}z\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ функция тока потенциального течения идеальной жидкости, обтекающего тор $S$ со скоростью $U$ на бесконечности с циркуляцией $\Gamma \mathrm{.}$ На бесконечности функция тока имеет следующее асимптотическое разложение:

$\Psi \mathrm{<mo>=</mo>}\frac{1}{2}{r}^{2}U+C+\mathcal{O}\left(\mathrm{1<mo>/</mo>}\sqrt{r^2+z^2}\right)$. (9)

Тогда рассмотрим область ${\Omega }_R\mathrm{,}$ ограниченную изнутри контуром $\partial S\mathrm{,}$ а снаружи $–$ границей полукруга $C_R\mathrm{,}$ состоящей из полуокружности и диаметра $\mathrm{<mo\; stretchy="false">(</mo>}-R\mathrm{,}R\mathrm{<mo\; stretchy="false">)</mo>}$ на оси симметрии (фиг. 3). Функция тока в полукруге $C_R$ удовлетворяет уравнению $D^2\Psi \mathrm{<mo>=</mo>0}$ и можно применить тождество (7). Граница области $\partial \Omega \mathrm{<mo>=</mo>}\partial S\cup C_R$ состоит из двух контуров, поэтому в интегральное тождество (7) добавится второй интеграл по границе полукруга $C_R\mathrm{.}$ В интеграле по контуру $\partial S$ следует учесть изменение знака, так как нормаль, внешняя к ${\Omega }_R$, будет внутренней к $S\mathrm{.}$ В результате для точек на контуре $M\in \partial S$ получим уравнение

$-2\pi \Psi \left(M\right)\mathrm{<mo>=</mo>}\tilde{A}V\left(M\right)+\tilde{B}\frac{\Psi }{r}\left(M\right)+I\mathrm{,}$ (10)

$I\mathrm{<mo>=</mo>}{\displaystyle \underset{C_R}{\int }}\frac{1}{r^{\prime }}\left[G\left(M\mathrm{,}{M}^{\prime }\right)\frac{\partial \Psi }{\partial n}\left(M^{\prime }\right)-\Psi \left(M^{\prime }\right)\frac{\partial G\left(M\mathrm{,}{M}^{\prime }\right)}{\partial n}\right]ds\mathrm{.}$

Фиг. 3.Обозначения к внешней задаче

Если же точка лежит вне контура $\partial S\mathrm{,}$ то согласно (5) будем иметь равенство

$-4\pi \Psi \left(M\right)\mathrm{<mo>=</mo>}\tilde{A}V\left(M\right)+\tilde{B}\frac{\Psi }{r}\left(M\right)+I\mathrm{.}$ (11)

Перейдем к вычислению интеграла $I$ по бесконечно удаленному контуру $C_R\mathrm{.}$ Интеграл $I$ состоит из интеграла $I_1$ по полуокружности и интеграла по отрезку прямой $I_2\mathrm{.}$ Интеграл $I_1$ вычисляется с помощью асимптотик функции Грина. Пусть $M^{\prime }\mathrm{<mo\; stretchy="false">(</mo>}R\sin \theta \mathrm{,}R\cos \theta \mathrm{<mo\; stretchy="false">)</mo>}$ $–$ точка полукруга и $M\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,0<mo\; stretchy="false">)</mo>}$ $–$ точка наблюдения. Проведем вихревой контур через точку $M\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,0<mo\; stretchy="false">)</mo>}$. Тогда для функции тока (9) имеем асимптотики

$\begin{array}{c}\Psi \left(M^{\prime }\right)\mathrm{<mo>=</mo>}\frac{1}{2}U{R}^2{{\displaystyle \sin }}^2\theta +C+\mathcal{O}\left(\mathrm{1<mo>/</mo>}R\right)\mathrm{,}\\ \frac{\partial \Psi \left(M^{\prime }\right)}{\partial n}\mathrm{<mo>=</mo>}UR{{\displaystyle \sin }}^2\theta +\mathcal{O}\left(\mathrm{1<mo>/</mo>}{R}^2\right)\end{array}$ (12)

для интеграла $I_1$ по дуге круга радиуса $R$ в пределе при $R\to \infty$ получим

$\begin{array}{c}{I}_1\mathrm{<mo>=</mo>}{\displaystyle \underset{0}{\overset{\pi }{\int }}}\frac{Rd\theta }{R\sin \theta }\left[\begin{array}{c}-\pi \frac{r^2}{R}{{\displaystyle \sin }}^2\theta \cdot UR{{\displaystyle \sin }}^2\theta -\\ -\frac{1}{2}U{R}^2{{\displaystyle \sin }}^2\theta \cdot \pi \frac{r^2}{R^2}{{\displaystyle \sin }}^2\theta \end{array}\right]\mathrm{<mo>=</mo>}\\ \mathrm{<mo>=</mo>}-\frac{3}{2}\pi U{r}^2{\displaystyle \underset{0}{\overset{\pi }{\int }}}{{\displaystyle \sin }}^3\theta d\theta \mathrm{<mo>=</mo>}-2\pi U{r}^2\mathrm{,}\end{array}$

$r^{\prime }\mathrm{<mo>=</mo>}R\sin \theta \mathrm{,}{z}^{\prime }\mathrm{<mo>=</mo>}R\cos \theta \mathrm{,}ds\mathrm{<mo>=</mo>}Rd\theta \mathrm{.}$

В интеграле $I_2$ при $r^{\prime }\to 0$ получаем

$\begin{array}{c}\begin{array}{c}G\left(M\mathrm{,}{M}^{\prime }\right)\simeq -\frac{\pi r^2{\left(r^{\prime }\right)}^2}{{\left(r^2+{\left(z^{\prime }\right)}^2\right)}^{\mathrm{3<mo>/</mo>2}}}\mathrm{,}\\ \frac{1}{r^{\prime }}\frac{\partial G}{\partial r^{\prime }}\simeq -\frac{2\pi r^2}{{\left(r^2+{\left(z^{\prime }\right)}^2\right)}^{\mathrm{3<mo>/</mo>2}}}\mathrm{,}\end{array}\\ \Psi \left(M^{\prime }\right)\mathrm{<mo>=</mo>}C+\mathcal{O}{\left(r^{\prime }\right)}^2\mathrm{.}\end{array}$

Подставляя эти выражения в интеграл $I_2$ в пределе $R\to \infty$, имеем

$I_2\mathrm{<mo>=</mo>}-{\displaystyle \underset{\infty }{\overset{\infty }{\int }}}\frac{C2\pi r^{2}d{z}^{\prime }}{{\left(r^2+{\left(z^{\prime }\right)}^2\right)}^{\mathrm{3<mo>/</mo>2}}}\mathrm{<mo>=</mo>}-4\pi C\mathrm{.}$

Подставляя в (10) вычисленное значение $I\mathrm{<mo>=</mo>}{I}_1+I_2$ и учитывая, что $\Psi \mathrm{<mo>=</mo>0}$ на границе $\partial S\mathrm{,}$ получим следующее уравнение:

$AV\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>4}\pi \left(\frac{1}{2}U{r}^2+C\right)\text{}\mathrm{.}$ (13)

К нему нужно присоединить условие на циркуляцию

${\displaystyle \underset{0}{\overset{l}{\int }}}V\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}ds\mathrm{<mo>=</mo>}-\Gamma \mathrm{.}$ (14)

Из системы уравнений (13) и (14) находим постоянную $C$ и распределение скорости на профиле $V\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>.}$ Функцию тока $\Psi \mathrm{<mo\; stretchy="false">(</mo>}M\mathrm{<mo\; stretchy="false">)</mo>}$ вне контура $\partial S$ найдем, пользуясь (11):

 $\Psi \left(M\right)\mathrm{<mo>=</mo>}\frac{1}{2}U{r}^2+C-\frac{1}{4\pi }\tilde{A}V\left(M\right)\mathrm{.}$ (15)

4. ЧИСЛЕННАЯ СХЕМА РЕШЕНИЯ ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ

Пусть $\partial S$ граница сечения тора меридиональной плоскостью. Так же как и для плоской задачи, для проведения численных расчетов вводится дискретизация контура $\partial S$ конечным числом точек $M_1\mathrm{,}{M}_2\mathrm{,}\dots \mathrm{,}{M}_N\mathrm{<mo>=</mo>}{M}_0\mathrm{.}$ Координаты точек $M$ задаются параметром $\zeta \in \mathrm{<mo\; stretchy="false">(</mo>0,1<mo\; stretchy="false">)</mo>,}$ так чтобы точке $M_i$ соответствовало значение $\zeta \mathrm{<mo>=</mo>}{\zeta }_i\mathrm{<mo>=</mo>}i\mathrm{<mo>/</mo>}N\mathrm{,}$ $i\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}N\mathrm{.}$ Параметр $\zeta$ и длина дуги на контуре $s$ связаны дифференциальным соотношением $ds\mathrm{<mo>=</mo>}f\mathrm{<mo\; stretchy="false">(</mo>}\zeta \mathrm{<mo\; stretchy="false">)</mo>}d\zeta$.

Таким образом, функция Грина от двух точек $M$ и $M^{\prime }$ на контуре $\partial S$ превратится в функцию двух аргументов $G\mathrm{<mo\; stretchy="false">(</mo>}\zeta \mathrm{,}{\zeta }^{\prime }\mathrm{<mo\; stretchy="false">)</mo>.}$ Она имеет логарифмическую особенность при ${\zeta }^{\prime }\to \zeta$ 

$G\left(\zeta \mathrm{,}{\zeta }^{\prime }\right)\mathrm{<mo>=</mo>2}r\ln \left(\frac{lf}{8\pi r}\sin \pi \left|\zeta -{\zeta }^{\prime }\right|\right)+4r+g\left(\zeta \mathrm{,}{\zeta }^{\prime }\right)\mathrm{,}$ (16)

где $g\mathrm{<mo\; stretchy="false">(</mo>}\zeta \mathrm{,}{\zeta }^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ аналитическая функция аргумента ${\zeta }^{\prime }$ периода 1.

Для периодических функций с такой периодической особенностью в работе [11] выведены специальные квадратурные формулы (формулы (17) и (18)). Для них в [12] доказана экспоненциальная сходимость. С их помощью получаем следующую аппроксимацию для матрицы $A$:

$\begin{array}{c}AV\left({\zeta }_i\right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^N}{A}_{ij}{V}_j\mathrm{,}\\ A_{ij}\mathrm{<mo>=</mo>}-\frac{l{f}_j}{N}\left(2r_i\beta \left(\left|i-j\right|\right)+G_{ij}\right)\text{}\mathrm{,}{G}_{ij}\mathrm{<mo>=</mo>}\left\{\begin{array}{c}G\left({\zeta }_i\mathrm{,}{\zeta }_j\right)\mathrm{,}i\ne j\mathrm{,}\\ 2r_i\ln \frac{l{f}_i}{8\pi r_i}+4r_i\mathrm{,}\end{array}\right.\end{array}$ (17)

и систему линейных уравнений, аппроксимирующих интегральные уравнения (13) и (14):

$\begin{array}{c}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^N}{A}_{ij}{V}_j-4\pi C\mathrm{<mo>=</mo>2}\pi r_i^{2}U\mathrm{,}\\ {\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^N}\frac{l{f}_j}{N}{V}_j\mathrm{<mo>=</mo>}-\Gamma \mathrm{.}\end{array}$ (18)

5. РАСЧЕТЫ ОБТЕКАНИЯ ТОРООБРАЗНЫХ ТЕЛ С ОСЕВОЙ СИММЕТРИЕЙ

Общее распределение скорости $V$ на границе осесимметричного тела представляется в виде линейной комбинации $V\mathrm{<mo>=</mo>}U{V}_0+\Gamma V_{\Gamma }$ двух фундаментальных распределений $V_0$ и $V_{\Gamma }$.

Распределение скорости $V_{\Gamma }$ на границе профиля и константа $C_{\Gamma }$ для единичной циркуляции и нулевой скорости потока на бесконечности находятся из системы уравнений, которая получается подстановкой в (18) $\Gamma \mathrm{<mo>=</mo>1,}U\mathrm{<mo>=</mo>0}$ 

${\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^N}{A}_{ij}{V}_j-4\pi C\mathrm{<mo>=</mo>0,}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^N}\frac{l{f}_j}{N}{V}_j\mathrm{<mo>=</mo>}-1.$ (19)

Аналогично находится распределение скорости $V_0$ на границе профиля и константа $C_0$ при $\Gamma \mathrm{<mo>=</mo>0}$ и $U\mathrm{<mo>=</mo>1}$

${\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^N}{A}_{ij}{V}_j-4\pi C\mathrm{<mo>=</mo>2}\pi r_i^2\mathrm{,}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^N}\frac{l{f}_j}{N}{V}_j\mathrm{<mo>=</mo>0.}$ (20)

В табл. 1 демонстрируется, как увеличивается точность расчетов с увеличением числа точек $N$ на контуре тора. Точки границы тора $\partial S$ задаются уравнением $r\mathrm{<mo>=</mo>2}+\sin t\mathrm{,}$ $z\mathrm{<mo>=</mo>0.5}\cos t\mathrm{.}$ В первой колонке приведены значения параметра $t\mathrm{,}$ в следующих трех колонках приведены расчеты значений скоростей $V_{\Gamma }\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ для $N\mathrm{<mo>=</mo>16,}$ 32 и 64 при обтекании контура с циркуляцией $\Gamma \mathrm{<mo>=</mo>1}$ и скоростью потока $U\mathrm{<mo>=</mo>0.}$ В последних трех колонках приведены расчеты значений скоростей $V_0\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ для $N\mathrm{<mo>=</mo>16,}$ 32 и 64 с циркуляцией $\Gamma \mathrm{<mo>=</mo>0}$ и $U\mathrm{<mo>=</mo>1.}$ При $N\mathrm{<mo>=</mo>16,}$ 32 и 64 достигается точность 3, 4 и 6 знаков соответственно.

Taблица 1

$\frac{t}{2\pi }$	Γ =1, U = 0			Γ = 0, U =1
	VΓ (N =16)	VΓ (N =32)	VΓ (N =64)	V0 (N =16)	V0 (N =32)	V0 (N =64)
0.125	–0.057078	–0.057069	–0.057068	1.36138	1.36147	1.36148
0.25	–0.065127	–0.065111	–0.065109	2.79678	2.79736	2.79743
0.375	–0.057078	–0.057069	–0.057068	1.36138	1.36147	1.36148
0.5	–0.100577	–0.100602	–0.100606	0.22108	0.22067	0.22062
0.625	–0.32877	–0.328961	–0.328984	–1.28767	1.28858	–1.28869
0.75	–0.859781	–0.859018	–0.858927	–3.91833	–3.91419	–3.91369
0.875	–0.32877	–0.328961	–0.328984	–1.28767	–1.28858	–1.28869
0.0	–0.100577	–0.100602	–0.100606	0.22108	0.22067	0.220625

Зная два фундаментальных распределения скоростей, можно построить общее распределение скорости на границе обтекаемого контура с произвольно заданными значениями $U$ и $\Gamma$: $V\mathrm{<mo>=</mo>}U{V}_0+\Gamma V_{\Gamma }$ и определить константу $C\mathrm{<mo>=</mo>}U{C}_0+\Gamma C_{\Gamma }$.

Функция тока вне контура $\partial S$ находится квадратурой (15). С помощью квадратурной формулы для периодических функций (17) получаем формулу, аналогичную для плоской задачи [11]:

$\Psi \left(r\mathrm{,}z\right)\mathrm{<mo>=</mo>}\frac{1}{2}U{r}^2+C+\frac{l}{4\pi N}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^N}G\left(r\mathrm{,}z\mathrm{,}{r}_j\mathrm{,}{z}_j\right)V_j{f}_j\mathrm{.}$ (21)

Для контуров типа крыла циркуляция находится с помощью постулата Чаплыгина-Жуковского $–$ линия тока сходит с острой кромки. Острую кромку удобно сгладить, так чтобы кривизна на ней была достаточно большой, но конечной. Тогда для выполнения постулата Чаплыгина-Жуковского достаточно потребовать, чтобы скорость в точке максимальной кривизны профиля была нулевой.

На фиг. 4 изображены линии тока обтекания торов, образованных вращениями профилей Жуковского. Профили направлены под различными углами атаки к оси $z$ — штрихпунктирная линия. Ось профиля – прямая, соединяющая вершину и острую кромку, направлена под разными углами к оси вращения (углы атаки). Углы атаки имеют значения $\mathrm{0,}\mathrm{2,}\mathrm{4,}\mathrm{6,}8$ и $10$ градусов. Вершины всех профилей находятся на одинаковом единичном расстоянии от оси вращения (штриховой отрезок).

Фиг. 4.Линии тока при обтекании торов, образованных профилями Жуковского, направленных под различными углами атаки.

Профили Жуковского определяются с помощью конформного отображения [16, с. 190].

$\frac{z-z_0}{z+z_0}\mathrm{<mo>=</mo>}{\left(\frac{Z-Z_0}{Z+Z_0}\right)}^{z_0\mathrm{<mo>/</mo>}{Z}_0}\mathrm{.}$

Оно отображает плоскость комплексного переменного $Z$ на физическую область комплексного переменного $z$. Окружность $Z\mathrm{<mo>=</mo>}{Z}_0-a{e}^{-i\beta }+a^{-i\mathrm{<mo\; stretchy="false">(</mo>}t-\beta \mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}t\in \mathrm{<mo\; stretchy="false">(</mo>0,2}\pi \mathrm{<mo\; stretchy="false">)</mo>}$, на плоскости $Z$ отображается на профиль Жуковского на плоскости $z$. Острая кромка профиля соответствует значениям $t\mathrm{<mo>=</mo>0}$ и $t\mathrm{<mo>=</mo>2}\pi$. Угол при острой комке равен $\tau \mathrm{<mo>=</mo>}\pi \mathrm{<mo\; stretchy="false">(</mo>2}-z_0\mathrm{<mo>/</mo>}{Z}_0\mathrm{<mo\; stretchy="false">)</mo>}$. Параметр $\beta$ является мерой изогнутости профиля. При $\beta \mathrm{<mo>=</mo>0}$ получается симметричный профиль. На фиг. 4 изображены профили Жуковского с параметрами $a\mathrm{<mo>=</mo>1.05,}{Z}_0\mathrm{<mo>=</mo>1.108,}\tau \mathrm{<mo>=</mo>0.1}$ и $\beta \mathrm{<mo>=</mo>}-0.1$. Параметр $z_0$ вычисляется по формуле $z_0\mathrm{<mo>=</mo>2}{Z}_0+\tau \mathrm{<mo>/</mo>}\pi$.

Циркуляцию можно найти с помощью постулата Чаплыгина-Жуковского описанным выше способом. Расход однородного потока через сечение радиуса $R$ равен $Q_0\mathrm{<mo>=</mo>}U\pi R^2$. Внутреннюю границу тора можно рассматривать как сопло, через которое течет жидкость. По мере увеличения угла атаки увеличивается циркуляция. С увеличением угла атаки при постоянной скорости потока $U\mathrm{<mo>=</mo>1}$ циркуляция растет и, как следствие, растет расход жидкости через тор (сопло) $Q$.

Функция тока на границе профиля равна нулю, а на оси вращения значение функции тока равно $-C$. Соответственно, расход жидкости через сопло - тор равен $Q\mathrm{<mo>=</mo>}-2\pi C$ и для относительного расхода получаем $Q\mathrm{<mo>/</mo>}{Q}_0\mathrm{<mo>=</mo>}-2C$.

На фиг. 5 представлены зависимости относительного расхода $Q\mathrm{<mo>/</mo>}{Q}_0$ через сопло от угла атаки, при значениях параметра $\beta \mathrm{<mo>=</mo>}-\mathrm{0.1<mo>;</mo>}-\mathrm{0.15<mo>;</mo>}-\mathrm{0.2<mo>;</mo>}-0.25$. Для каждого значения $\beta$ вычисляются расходы для углов атаки $\theta {\mathrm{<mo>=</mo>0}}^{\circ }\mathrm{<mo>;</mo>}{2}^{\circ }{4}^{\circ }{6}^{\circ }{8}^{\circ }$ и $10$ градусов. Через них проводится кривая на фиг. 5.

Фиг. 5.Зависимости расхода через сопло от угла атаки, при значениях параметра.

Программа реализована в пакете Wolfram Mathematica [17]. Точность 5 знаков достигается при N = 64. Для расчета обтекания всех шести профилей и построения графика зависимости $Q\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo><mo>/</mo>}{Q}_0$ требуется менее секунды расчетного времени.

Построенные зависимости показывают монотонное увеличение расхода при росте угла атаки. Как видно из графиков, при увеличении параметра $\beta$ расход тоже увеличивается от значения, равного 1 при нулевом угле атаки, до значения, большего 4-х при угле атаки ${10}^{\circ }$.

6. ОБТЕКАНИЕ ОСЕСИММЕТРИЧНОГО ТЕЛА

Пусть $\partial S$ сечение односвязной осесимметричной поверхности меридиональной плоскостью. Специфика этой задачи состоит в том, что контур $L\mathrm{,}$ по которому будет вестись интегрирование, не будет замкнутым (фиг. 6). Для этого случая циркуляция $\Gamma$ и постоянная $C$ будут равны нулю. Распределение скорости на контуре $L$ находится из уравнения (13), а функция тока вне контура $–$ из (21). Поскольку контур $L$ незамкнут, то подынтегральная функция непериодическая. Используются по-прежнему линейные уравнения (18). Однако остаточный член квадратурной формулы (17) для матрицы $A$ имеет не экспоненциальную, а степенную оценку $R_N\sim \frac{1}{12N^2}\mathrm{,}$ где $N$ $–$ число точек на контуре.

Фиг. 6.Обозначения

Проверка на известных точных решениях показывает также достаточно быструю сходимость метода и в этом случае. В табл. 2 приводятся сравнения точного решения при отношении осей $\chi \mathrm{<mo>=</mo>1<mo>/</mo>5}$ с расчетами при различном числе точек на границе эллипсоида. Точки эллипсоида с отношением осей $\chi$ определяются из уравнений $r\mathrm{<mo>=</mo>}\chi \sin t\mathrm{,}$ $z\mathrm{<mo>=</mo>}\cos t\mathrm{.}$ Точное решение имеет вид

$\begin{array}{c}V\mathrm{<mo>=</mo>}\frac{2}{2-A\mathrm{<mo\; stretchy="false">(</mo>}\chi \mathrm{<mo\; stretchy="false">)</mo>}}\frac{r}{\sqrt{r^2+{\chi }^4{z}^2}}\mathrm{,}\\ A\mathrm{<mo\; stretchy="false">(</mo>}\chi \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>2}\frac{{\chi }^2}{1-{\chi }^2}\left(\frac{1}{2\sqrt{{\chi }^2-1}}\ln \frac{1+\sqrt{1-{\chi }^2}}{1-\sqrt{1-{\chi }^2}}-1\right)\text{}\mathrm{,}\chi \mathrm{<mo><</mo>1,}\\ A\mathrm{<mo\; stretchy="false">(</mo>}\chi \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>2}\frac{{\chi }^2}{1-{\chi }^2}\left(\frac{\sqrt{{\chi }^2-1}}{\sqrt{{\chi }^2-1}}-1\right)\text{}\mathrm{,}\chi \mathrm{<mo>></mo>1.}\end{array}$

Taблица 2

$t/\left(2\mathrm{\pi }\right)$	Vexact	V (N = 16)	V (N = 32)	V (N = 64)	V (N = 128)
0.125	0.953 762	0.953 939	0.953 616	0.953 733	0.953 753
0.250	1.033 554	1.04132	1.033 337	1.033 595	1.033 559
0.375	1.055 505	1.053 652	1.055 391	1.055 553	1.055 511
0.500	1.059 121	1.062 335	1.059 515	1.05917	1.059 127

При $N\mathrm{<mo>=</mo>16,}$ 32, 64 и 128 достигается точность 3, 4, 5 и 6 знаков соответственно.

ВЫВОДЫ

Построенная численная схема для решения задачи обтекания осесимметричных торов аналогична схеме плоской задачи обтекания профилей [11]. Она сводится к решению системы линейных уравнений для значений скорости в точках сетки на границе обтекаемого тела. Матрица линейной системы находится по квадратурной формуле с погрешностью, которая убывает в зависимости от числа элементов сетки $N$ по экспоненте $e^{-CN}$. Доказательство экспоненциального убывания погрешности квадратурной формулы с периодической особенностью доказывается в [12]. Функция тока в численной схеме метода граничных элементов ранее не использовалась, хотя она обладает многими преимуществами перед численной схемой для потенциала. Во-первых, схема значительно проще существующих, во-вторых, точность ее значительно выше и, в-третьих, схема удобна для построения линий тока.

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

[1] Работа выполнена при финансовой поддержке РНФ, проект $№$ 124012500442-3.