Численно-аналитические методы решения задач об обтекании препятствий под поверхностью весомой жидкости с образованием солитона Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 5, № 2, 2000

ЧИСЛЕННО-АНАЛИТИЧЕСКИЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ ОБ ОБТЕКАНИИ ПРЕПЯТСТВИЙ ПОД ПОВЕРХНОСТЬЮ ВЕСОМОЙ ЖИДКОСТИ С ОБРАЗОВАНИЕМ

СОЛИТОНА*

В. П. Житников, Н.М. Шерыхалина Уфимский государственный авиационный технический университет, Россия e-mail: pavel@zhitnik.ugatu.ac.ru

The improved solution method for the problem of flow past point obstructions is suggested. This allows to solve more complex problems of flows past obstructions of limiting size. The singularities of the sought function are allowed for by adding terms containing fractional power functions. The results of numerical investigation of the problem of flow past an obstruction shaped as an ellipsoidal semicylinder and located on a channel bottom are shown. The solutions for different ellipsis semiaxes relations including limiting cases such as Stokes wave are obtained.

Одним из наиболее распространенных численно-аналитических методов решения задач со свободными границами, использующих возможности теории функций комплексного переменного, является метод Леви — Чивиты [3]. Усовершенствование этого метода путем учета поведения решения вблизи особых точек было предложено в ряде работ [4-6, 8-10]. В данном сообщении предлагаются приемы, позволяющие сделать следующий шаг в этом направлении, т. е. выделить особенности более высокого порядка. Проведен анализ результатов численного эксперимента, который показал, что учет дополнительных особенностей позволяет существенно увеличить скорость сходимости метода и получить более точные данные при меньших затратах времени.

1. Постановка задачи

В качестве примера рассматривается плоская задача симметричного обтекания донного препятствия ограниченным потоком идеальной несжимаемой весомой жидкости со свободными границами, причем определяются апериодические решения (солитоны). Пусть скорость на бесконечности — V0, асимптотическая толщина струи — h, ускорение силы тяжести g действует вертикально вниз (рис. 1, а).

* Работа выполнена при поддержке федеральной целевой программы “Государственная поддержка интеграции высшего образования и фундаментальной науки на 1997-2000 годы”.

© В. П. Житников, Н. М. Шерыхалина, 2000.

На свободной поверхности ОБ значение модуля вектора скорости жидкости V связано с высотой точки у уравнением Бернулли

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

Обычно такие задачи решаются следующим образом. Область течения физической плоскости £ = х+%у конформно отображается на четверть или половину кольца [6] плоскости изменения параметрического переменного £. На внутреннюю и внешнюю дуги окружностей отображаются границы со сложными условиями. Поскольку областью течения на плоскости комплексного потенциала ю является полоса, то конформное отображение ю (£) определяется через эллиптические функции. Искомая функция £ (С) представляется в виде суммы аналитических функций с известными особенностями и ряда Лорана.

Основная идея модификации метода решения заключается в следующем. Используется решение задачи обтекания диполя. Для этой задачи

где п = N/Щ, N — дипольный момент. Эта же функция применяется для отображения некоторой области плоскости £ на полосу плоскости комплексного потенциала т, которая соответствует области течения в задаче обтекания препятствия. Прообразом данной полосы на плоскости £, определяемым через конформное отображение (1.2), является криволинейное кольцо, в котором внешняя граница — дуга окружности единичного радиуса, а внутренняя — некоторая криволинейная дуга (рис. 1, б). Путем численного исследования было установлено, что эта криволинейная дуга очень мало отличается от дуги окружности. Возникает возможность, не меняя конформного отображения т (£), поставить условие о

а 6 е

(1.1)

(1.2)

/г

£/1Ы±

Б

Ь

А 0

х

-1

0

1 В

Рис. 1. Область течения на различных плоскостях: а — физическая плоскость; б — плоскость изменения параметрического переменного; в — схема расчета обтекания эллиптического цилиндра

с большой горизонтальной полуосью.

форме границы обтекаемого тела на этой криволинейной границе. Это позволяет обойтись минимальными изменениями при переходе от задач обтекания точечных препятствий к препятствиям конечных размеров и построить весьма простые алгоритмы решения задач. Кроме того, поскольку диполем при N < 0 моделируется обтекание тела малого размера, то число добавочных слагаемых решения, необходимых для удовлетворения условию непротекания на границе тела, оказывается небольшим.

Представим границу области на плоскости Z в параметрическом виде Z = r (а) вга. Величина а задается в диапазоне [0, п/2]. Значение r (а) определяется из условия равенства нулю мнимой части w (а)

/ 2r (а) \ n { 1 .Л

arctan ------ sin а —- — r (а) sin а = 0. (1.3)

V1 — r2 (а) ) 4 \r (а) )

Отсюда при а ^ 0 нетрудно определить образ критической точки

С = 5 = \\-----------п/2 . (1.4)

1/2 - п/2+ 2^/1 - п/2

Уравнение (1.3) имеет единственное решение при г Є (0, 1), п < 0 и может быть решено численно методом Ньютона. Краевое условие, определяющее форму границы обтекаемого тела, в этом случае примет вид

/ Iх [г И ,^] ,У [г И ,^]} = 0- (1-5)

Для решения задачи используются методы двух типов — прямого конформного отображения г (£) и с использованием функции Жуковского.

2. Метод прямого конформного отображения

Функция £ (£) является аналитической в полукруге |£ | < 1, > 0 нечетной (в силу сим-

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

£ (С) = к [£о (С) + £1 (С) + £2 (С) + £з (С)] • (2-1)

Функция £0 (£) = — 1п—+7 отображает полуокружность на полосу единичной ширины.

п 1 - С

Функция £1 (£) = £ (С)/к — £0 (С) — £2 (С) — £3 (С) нечетна, имеет чисто действительные значения на действительном диаметре, аналитична внутри полукольца и непрерывна на границе. В связи с этим ее можно аналитически продолжить на все кольцо и искать в виде ряда Лорана с действительными коэффициентами

ГО

£1 (с )= Е с2„,+1С'2т+1- (—•—)

Коэффициенты с2т+1 подбираются так, чтобы выполнялись условия (1.1) и (1.5) (при этом V = —— , х = у = ). При численном решении количество слагаемых суммы

dz

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

Функции г2 (С) и г3 (£) введены для учета особенностей решения г (£) при £ = 1 ив точке излома свободной поверхности и определяются аналогично [8]:

Наличие степенной особенности с показателем 2/3 в (2.4) означает, что свободная граница в точке излома С образует угол, равный 2п/3, что соответствует известной особенности волны Стокса. Эта величина получена в предположении, что предельное значение угла в при подходе к точке С существует и не равно нулю. При нулевом предельном значении в можно ожидать существование других режимов с критической точкой С.

Здесь также учитывалось, что, как и г (£), функции г2 (С) и г3 (£) должны быть нечетными и иметь действительные значения на действительной оси. Тем самым определен вид функции (2.1).

3. Использование функции Жуковского

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

где в — угол наклона вектора скорости к оси х, т = 1п(У/^).

Поскольку на горизонтальном и вертикальном диаметре в = 0, то ш {£) является четной функцией. Представим ш в виде суммы

здесь (С) — функция Жуковского для аналогичной задачи о невесомой жидкости

(2.3)

(2.4)

где а — первый корень трансцендентного уравнения Стокса

п 1 2 = Ё? ’

0 < а < 1.

(2.5)

(3.1)

ш (С) = шо К) + ш1 (С) + ш2 (С) + шз к),

(3.2)

С2 - з2

(3.3)

— ряд Лорана, коэффициенты которого, как и в (2.2), подбираются так, чтобы выполнялись условия (1.1) и (1.5):

(3.4)

(множитель 1 — £2 вводится для того, чтобы и (1) = 0).

Слагаемые и2 (£) и и3 (£), как и £2 (С), Сз (С) в (2.1), введены для учета особенностей решений в точках £ = ±1 и £ = ±* с включением членов более высокого порядка

^2 (С) = ¿В1

1 - с2

+ ¿в

1 - с2

1 - с2

2 \ а+1

(3.5)

^з (С) = 3 1п“+^ + ¿С1

1 + С2 2

2\ в

1

(3.6)

Значение а определяется из решения уравнения (2.5). Коэффициенты В2 и В3 можно найти при подстановке (3.5) в (1.1) и предельном переходе при £ ^ 1:

В2 = 3 В2 ОП. Вз = аВі.

(3.7)

Аналогично, при подстановке (3.6) в (1.1) и предельном переходе при £ ^ * можно получить равенства

Сі

&2т ( —1)

1

т=-го

1п| Пй) + 1п(1 -

П 1

(в + 1) р2 = . (3.8)

Таким образом, найден вид функции (3.2).

а

т

3

4. Численное решение

Численно задача решается методом коллокаций. Для этого применяются алгоритмы, аналогичные описанным в [8]. В бесконечной сумме (2.2) сохраняется конечное число слагаемых, а уравнения (1.1) и (1.5) выполняются в дискретных точках ат = пт/2М, т = 1,..., М, и ат = пт/2М_, т = 0,..., М_, соответственно. Для этого вычисляются значения

Жт = ^ (ат) , Ут = ^ (^т) . (4.1)

Тем самым получается система М + М_ + 1 нелинейных уравнений, которая решается относительно параметров 2т+1, т = — М_,...,М — 2, А, п численно методом Ньютона с регулированием шага. При расчете течений со стоксовской особенностью в число искомых параметров включаются А2, Гг (за счет уменьшения числа искомых коэффициентов 2т+1). При расчетах значения к и -0 принимались равными единице.

Таким же способом задача решается при применении функции (3.2). В этом случае

С

1 С т

Ут = г (С) = е”к> -¿-¿С (4.2)

0

Численное интегрирование производится с помощью квадратурной формулы Симпсона.

Для проверки разработанных методов и программ было рассмотрено обтекание препятствия в виде кругового полуцилиндра.

^ . ~Г~ (ат)

аг

Эта задача решалась многими авторами (см., например, [2, 7]) и служит хорошим тестовым примером. На рис. 2 приведен фрагмент зависимости ус (Гг) в увеличенном масштабе. Тонкими линиями на график нанесены результаты, полученные К. Е. Афанасьевым

[1]. Видно, что отличие результатов весьма мало.

На рис. 3 показаны зависимости числа Гг-1/2 от ус — Я при различных радиусах обтекаемых препятствий. Тонкими линиями на рисунке нанесены результаты приближенного решения, о котором будет сказано ниже.

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

(2.2), (3.4) для Я < 1, как правило, можно было ограничить тремя —шестью при точности порядка 10-6.

Ус с*-

К=

0 0.2 у ' 5.3

У. >

1.2 1.25 1.3 1.35 Гг

Рис. 2. Фрагмент зависимости высоты волны ус от числа Фруда для обтекания кругового

полуцилиндра.

0.1- Т. -1/2 Рг

! >У- 'С

1 ^ - 0.6 -

V- 10

л~ ЯП Рг 1=1 _ - 0.2 -

Г

2 3 ус-я

Рис. 3. Зависимости числа Рг 1/2 от разности ус — Я при обтекании кругового полуцилиндра

с различными радиусами Я.

5. Обтекание препятствия в виде эллипсоидального полуцилиндра

При решении задачи в качестве краевого условия (1.5) на поверхности обтекаемого препятствия используется уравнение эллипса

2 2

X2 у2

а

+ = 1, (5Л)

где а и Ь — полуоси эллипса.

Рассмотрим случай, когда Ь фиксированно, а а изменяется от нуля до бесконечности (рис. 4). Каждая кривая для конкретных значений а и Ь начинается при Рг = ж, что соответствует случаю обтекания невесомой жидкостью, а заканчивается на кривой, отвечающей волне предельной высоты с изломом свободной поверхности на угол 120° (волне Стокса). При а = 0 не требуется никаких изменений алгоритма решения, кроме модификации условия (5.1), которое в этом случае представляется в виде

х (от) =0 , т < М- , у (ом_) = Ь,

где М- — количество членов с отрицательными степенями, Ом_ = п/2.

На рис. 4 кривая, соответствующая а = 0, находится левее и выше всех остальных. При а ^ ж для любой конечной зоны с центром в начале координат граница препятствия спрямляется. Поэтому течение определяется решением задачи о солитоне над ровным дном. Однако при этом весь поток поднят над истинным дном на величину Ь. В связи с этим значения числа Рг и высоты солитона должны быть пересчитаны с учетом изменения масштаба. Для этого наряду с числом Рг введем “местное” значение числа Фруда

Рг! У

где VI и Н\ — скорость течения и глубина потока над препятствием в сечении Б, расположенном достаточно далеко от начала координат, но еще дальше от точки пересечения

Рис. 4. Зависимости числа 1/Рг от высоты волны ус для обтекания эллиптического полуцилиндра с различной величиной горизонтальной полуоси.

эллипса с осью абсцисс (рис. 1, в), например, а1/2). В этом случае зависимость

Рг1 (^) <5/2>

определяется из решения задачи о солитоне над ровной поверхностью [8]. Необходимо связать значения к1 и У1 с к и У0. Для этого, аналогично [7], запишем уравнение сохранения

расхода жидкости и уравнение Бернулли для двух сечений (Б и бесконечно удаленного):

Последнее уравнение, имеющее единственное решение при ц > 1, может быть разрешено относительно ц методом Ньютона. Общий порядок решения при этом следующий. Задается число Рг1. Множество значений Рг1 ограничено с одной стороны предельным решением типа солитона Стокса (Рг1 = 1.290890455 [8]), с другой — солитоном, выродившимся в

Далее находим Ь/к = уц + в и Рг из (5.6). Кривая, соответствующая этому решению, обозначена на рис. 4 1.

Однако равномерный поток (в этом случае у = (ус — Ь)/к1 = 1) может существовать не только при Рг1 = 1, но и при любых значениях Рг1. Это приводит к возможности продолжения предельного решения, которое получается путем изменения Рг1 от единицы до бесконечности, решении уравнения (5.7) относительно ц и определении Ь/к = ц + в и Рг из (5.6). Данная ветвь решения обозначена на рис. 4 2.

Как видно из рис. 4, решения задачи с а = 2.5 и а = 5 проходят очень близко от кривых 1 и 2. Расчеты показывают, что кривые с а > 10 практически совпадают с кривыми 1 и 2. Поэтому их можно считать предельным решением задачи об обтекании эллипса. Тем самым оправдывается употребляемая в [1, 2] терминология “решения, ответвляющиеся от равномерного потока или от солитона”.

При решении задачи с 0 < а < 10 количество коэффициентов с отрицательными степенями в сумме (2.2), как правило, можно было ограничить числом 12-18 при точности порядка 10-5.

Подобные приближенные решения имеет смысл применить к другим предельным случаям, когда какой-нибудь размер обтекаемого препятствия уходит в бесконечность. Рассмотрим, например, случай а = Ь ^ ж, т. е. обтекание кругового цилиндра большого радиуса. Как показал численный эксперимент, прямое применение формул (5.3)-(5.7) к такой задаче дает довольно грубые результаты. Значительного улучшения удается достичь, если

(5.3)

Отсюда нетрудно получить уравнение

(5.4)

(5.5)

VVі/й

С учетом этого равенство (5.5) приводится к виду

(5.6)

(5.7)

равномерный поток (Рг1 = 1). Затем находим ц из (5.7) и у = (ус — Ь)/к1, обращая (5.2).

учесть, что при движении вдоль нормали к направляющей кругового цилиндра скорость

в сечении струи не постоянна, а меняется в соответствии с условием сохранения нулевой

циркуляции:

= VI Я + к1 , 0 < х < к1,

х 1 Я + х ~

где VX — скорость жидкости на расстоянии х от поверхности кругового цилиндра, V, к1 —

скорость на свободной поверхности и толщина струи вблизи верхней точки цилиндра. В

этом случае уравнение сохранения расхода примет вид

Щ = V! (Я + к1)1пЛ+ ^ , (5.8)

R

а местное значение числа Фруда

Т7_ _ V _ ГГ „ _ h1 .. _ R ^ ГЛ

Fri !—— , . , n ^ . (5-9)

vghi (^ + n)ln(l + h h

Уравнение Бернулли (5.3) принимает вид

Fri ((^ + n)2 ln^ 1 + ^ + 1^ - 2^ П + 1 = 0. (5.10)

Из рис. 3 видно, что тонкие кривые, соответствующие этому приближению (в предположении, что связь относительной высоты струи в точке C y = (yc — R)/h1 с местным значением числа Фруда Fr1 определяется из (5.2) так же, как и для солитона над ровным дном), проходят довольно близко от расчетных кривых. Это особенно проявляется на части кривых, соответствующих равномерному потоку y = 1 (до точечной кривой Fri = 1). На второй части кривых значения числа Fr совпадают достаточно хорошо, но предельная высота солитона у приближенного решения несколько больше расчетной. Это объясняется различием эпюр скорости у солитонов, расположенных над ровным дном и над криволинейной поверхностью.

Рассмотрим другой предельный случай, когда a = 0, а b изменяется от нуля до бесконечности (обтекание пластины). Очевидно, что вблизи верхней точки такого препятствия при b ^ ж образуется некоторое предельное течение. Задача состоит в том, чтобы найти, при каких b/h характеристики течения перестают зависеть от этого отношения. На рис. 5 даны зависимости местного числа 1/Fr1 от относительного расстояния h' между точкой C и верхней точкой препятствия:

Fri = . Vi , h' = (Ус —Ь) Vi, (5.11)

Vd (Ус — b) hV0

где Vi — скорость на свободной поверхности в точке, ордината которой совпадает с максимальной ординатой обтекаемого препятствия b. Связь между этими величинами и общими параметрами потока Fr, n = (yc — b)/h , в = b/h определяется с помощью уравнения Бернулли

V 2 V 2 тт + gh = -2- + gb,

откуда имеем

Vi \2 2

vi =!— FT2(в — 1>- <5-12>

а = о

1/РГ1

5 Ю

0.3

0.9

1.5

2.1

Ь'

Рис. 5. Зависимости местного числа 1/Рг1 от относительной высоты волны Ы для обтекания

Из рис. 5 видно, что все кривые с в > 5 практически совпадают. Это означает, что течение в верхней части практически не зависит от нижней.

Таким образом, с помощью несложного видоизменения методов решения задачи об обтекании диполя удалось решить задачи более сложного вида: обтекания неровности на дне в виде полукруга и полуэллипса. Расчеты показали, что в весьма широком диапазоне изменения длины полуосей эллипса при точности порядка 10-5 количество добавочных слагаемых в решении не превышает 12-18.

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

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

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

[1] АФАНАСЬЕВ К. Е. Приближение нелинейной уединенной волны Тр. VI науч. школы “Гидродинамика больших скоростей”. Чувашский гос. ун-т, Чебоксары, 1996, 3-10.

[2] ГУЗЕВСКИЙ Л. Г. Обтекание препятствий потоком тяжелой жидкости конечной глубины. В “Динамика сплошной среды с границами раздела”. Чувашский гос. ун-т, Чебоксары, 1982, 61-69.

[3] Гуревич М. И. Теория струй идеальной жидкости. Наука, М., 1979.

пластины различной высоты.

Тогда

(5.13)

[4] Житников В. П. Гравитационные волны на ограниченном участке поверхности жидкости. ПМТФ, 37, №2, 1996, 83-89.

[5] Житников В. П., Терентьев А. Г. Струйное обтекание гибкой оболочки потоком идеальной жидкости. Изв. АН СССР. МЖГ, №6, 1982, 43-48.

[6] Киселев О. М., Котляр Л. М. Нелинейные задачи теории струйных течений тяжелой жидкости. Казанский гос. ун-т, 1978.

[7] Маклаков Д. В. Нелинейные задачи гидродинамики потенциальных течений с неизвестными границами. Янус-К, М., 1997.

[8] Шерыхалина Н. М. Разработка численных алгоритмов решения задач гидродинамики с особыми точками на свободной поверхности и экспериментальное исследование скорости их сходимости. Деп. в ВИНИТИ, №2550-В95, М., 1995.

[9] HüNTER J. K., Vanden-Broek J.-M. Accurate computations for steep solitary waves. J. F. M, 136, 1983, 63-71.

[10] WILLIAMS J. M. Limiting gravity waves in water of finite depth. Phil. Trans. Roy. Soc. bond, A302, No. 1466, 1981, 139-188.

Поступила в редакцию 5 декабря 1998 г.