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

СИСТЕМЫ И

ПРОЦЕССЫ

УПРАВЛЕНИЯ

УДК517.9 : 532.5

ПРИМЕНЕНИЕ ИТЕРАЦИОННЫХ МЕТОДОВ К РЕШЕНИЮ ВНЕШНИХ ЗАДАЧ ГИДРОДИНАМИКИ

КОЛОСОВА С.В., ЛАМТЮГОВА С.Н.,

СИДОРОВ М.В.____________________________

Рассматривается применение методов R-функций, последовательных приближений и метода Г алеркина-Петрова к расчету осесимметричных стационарных течений вязкой несжимаемой жидкости (обтекание конечных тел вращения).

Введение

Актуальность исследования. Многие явления, наблюдаемые в атмосфере и океане, а также проблемы гидроаэродинамики, теплоэнергетики, химической кинетики, биомедицины можно изучать в рамках модели несжимаемой вязкой жидкости. Задачи, представляющие практический интерес, как правило, описываются уравнениями Навье-Стокса [1 - 3], существенной особенностью которых является нелинейность, а также наличие малого параметра при старшей производной (величина, обратная числу Рейнольдса). Кроме того, задачи для уравнений Навье-Стокса часто приходится решать в областях сложной геометрии, а область может быть и бесконечной (задачи обтекания тел, течения жидкости в трубах и пр.). В большинстве случаев при численном решении задач обтекания условия на бесконечности сносятся на некоторый контур, расположенный достаточно далеко от обтекаемого тела, что приводит к чрезмерным затратам ресурсов ЭВМ.

Существует множество подходов к расчету вязких течений. В основном эти подходы используют метод конечных разностей и метод конечных элементов [4 -12]. Эти методы просты в реализации, но не обладают необходимым свойством универсальности - при переходе к новой области (особенно неклассической геометрии) необходимо генерировать новую сетку, а часто и заменять сложные участки границы простыми, составленными, например, из отрезков прямых. Точно учесть геометрию области, а также краевые условия (в том числе и условие на бесконечности) можно, воспользовавшись конструктивным аппаратом теории R-функций акад. НАН Украины В.Л. Рва-чева [13]. Метод R-функций в задачах гидродинамики использовался в работах [14 - 18], но рассматривались задачи расчета течений идеальной жидкости [14] или же вязкой в ограниченных областях [15 - 17], или при наличии винтовой симметрии [18].

РИ, 2012, № 3

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

Цели и задачи исследования. Целью данной работы является создание современного и эффективного метода математического моделирования стационарных задач обтекания тел вращения. Здесь не обсуждается степень строгости, условия применимости использованных уравнений движения жидкости, они рассматриваются как математические модели, подлежащие численной алгоритмизации. Для достижения поставленной цели необходимо решить следующие задачи: на основании методов теории R-функций построить полную структуру решения краевой задачи для функции тока; заменить исходную нелинейную задачу последовательностью линейных краевых задач; для решения линейных задач на каждом шаге итерационного процесса разработать численный алгоритм на основании метода Г алеркина-Петрова.

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

Рассмотрим стационарное обтекание тела вращения потоком вязкой несжимаемой жидкости. Будем считать, что в пространстве введена декартова система координат (x, y,z), а обтекаемое тело образовано вращением вокруг оси Oz фигуры Q, лежащей в плоскости Oxz (фигура Q односвязная, конечная и симметричная относительно оси Oz ). Кроме того, предположим, что поток жидкости равномерный, его скорость равна UM и он сонаправлен с осью Oz . Такие течения удобно рассматривать в сферической

системе координат (r, 9, ф):

x = rsin 9 cos ф, y = rsin9sinф , z = rcos9,

0 < r <+да , 0 <9<л, 0 <ф< 2п .

Стационарные уравнения Навье-Стокса в сферической системе координат имеют вид [3]:

dvr v9 5vr ^ 5vr v9 + Уф 1 dp

dr r 59 rsin9 5ф r p 5r

+v

Avr -

2vr

2

r2 59

2v9 ctg9 2

r2 r2 sin 9 Зф ^

2

5v9+ v9 5у9 + Уф 5уе + уг vQ v9ctg9

r 5r r 59 rsin9 5ф r r

115p ( 2 5vr

=-----— + v| Av9 +—-- r

p r 59

r2 59

13

Vq 2cos9

r2 sin2 9 r2 sin2 9 дф ) ’

(1)

4.=о, ду

Sn

= 0.

5Q

(4)

^ф V9 vф vф vr

vr—— + ——— + — -----+—----+

dr r 59 rsin9 5ф r

v9 vф ctgQ

1 1 dp

p rsin9 5ф

(

+ v

v,

ф

ЛVФ 2 . 2,

v r sin I

+ 2 •

5vr + 2cos9 5vq

r2 sin

r“ sin 9 дф

d-(f2Vrsin9)+)Г sin9)+ ( ^

2 9 дф

= о

5r' '59 5ф

Здесь vr, vq , Vф - радиальная, угловая и осевая компоненты скорости жидкости соответственно; p -давление; v - кинематический коэффициент вязкости; p - плотность жидкости,

2

А = 4Щ r2 А| +

r2 5r ^ 5r

1

д

sin 9— | +

r2sin 9 д9 V дві -2

1

д

r2 sin2 9 дф2

Совокупность уравнений (1) представляет собой систему нелинейных уравнений в частных производных относительно неизвестных функций vr , Vq , Vф , p . Ввиду сложности системы (1) ее прямой численный анализ затруднен.

В сделанных выше предположениях относительно обтекаемого тела осевая компонента скорости Vф равна нулю, vr, vq и p являются функциями только от r и 9 , а уравнение неразрывности (четвертое уравнение в (1)) интегрируется введением функции тока у по формулам

Vr =-

1 ду

' sin 9 д9

v9 =-

1 ду rsin9 дr

(2)

Исключив из оставшихся двух первых уравнений перекрестным дифференцированием давление, для функции тока у = у(r, 9) получим нелинейное уравнение четвертого порядка [19]:

^2 1 Г ду дБу дудБу|| +

V у = r2 sin 9І5в^Т-^Г ~S^J +

r2 sin 9

„ _ ду 2 ду'і

2ctg9^T - 7 791 бу вне Q, (3)

r

где n - внешняя к д. нормаль.

Условие на бесконечности имеет вид

у ~ — ида r2 sin2 9 при r ^да (5)

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

Соответствующая (3) - (5) линейная задача (приближение Стокса) была рассмотрена в [20].

Итак, для расчета течения около рассматриваемого тела вращения нужно решить краевую задачу (3) - (5).

2. Построение структуры решения

Для решения задачи (3) - (5) воспользуемся методом R-функций: с помощью конструктивных средств теории R-функций построим структуру решения краевой задачи (3) - (5), т.е. пучок функций, точно удовлетворяющий краевым условиям на д. и условию при r ^да.

Пусть вне Q известна достаточно гладкая функция ю (r, 9), обладающая следующими свойствами:

1) ю (r, 9) > 0 вне Q ;

2) ю (r, 9) = 0 на д.;

.ч дю (r,9) „

3) —^-^ = -1 на д.,

дп

где п - вектор внешней нормали к д..

Для областей произвольной формы, ограниченных кусочно-гладким контуром, такая функция ю может быть построена в виде единого аналитического выражения благодаря использованию R-функций [13].

Введем в рассмотрение достаточно гладкую функцию У = fM (x) [21], удовлетворяющую следующим требованиям:

а) fM(0) = 0;

б) fM (0) = 1;

в) fM (x) > 0 Vx > 0;

где Бу

д2у sin 9 д ( 1 ду Л д9 r2 дві, sin 9 д9

Б2у = Б (Бу).

Уравнение (3) следует дополнить условиями на д. и на бесконечности (при r ).

Если граница обтекаемого тела неподвижна и непроницаема, то из условий прилипания следуют такие краевые условия:

г) fM (x) = 1 Vx > M (M = const > 0).

Условиям а) - г) удовлетворяет, например, функция

fM(x) =

1 - exp Mx , 0 < x < M; x - M

1, x > M.

Кроме того, очевидно, что такая fM (x) є Сда [0, + да).

14

РИ, 2012, № 3

Обозначим

raM(r, 0) = fM[ra(r, 0)] . (6)

Легко проверить, что функция юм(г, 0) удовлетворяет условиям 1) - 3). Кроме того,

Юм(г, 0) = 1, если ю (r, 0) > M.

Заметим, что это условие означает, что если функция ю (r, 0) монотонно возрастает при удалении от dQ, то функция юм(г, 0) вида (6) отлична от единицы лишь в некоторой кольцеобразной области

vE2u = Б(уо + u) -vE2yо вне Q , (9)

I n du uL„ = 0 , —

ldQ ’ dn

= 0,

dQ

(10)

lim u = 0. (11)

r^+ro

Заметим, что в силу свойств функции у0 вида (7) в области {ю (r, 0) > M} имеем

E2y0 = 0 ,

{0 <ю(г, 0) <M},

которая содержится во внешности Q и прилегает к dQ.

Известно [19], что уравнение (3) имеет частное решение вида

у = у(г, 0) = (r2 + c2r-1 ^ sin2 0 ,

где c1 , c2 - произвольные постоянные, причем такая у удовлетворяет и уравнению Еу = 0. Тогда функция

У0 =У0(г,0) = юМ(г,0):2uot r2sin2 0 (7)

удовлетворяет краевым условиям (4) и условию (5). Кроме того, на функции у0 вида (7) в области {ю(r, 0) > M} уравнение (3) обращается в тождество.

Из сказанного выше следует

Теорема. При любом выборе достаточно гладких функций Ф1 и Ф2 (Ф1 ^ 0 при r ^ +оо) краевым условиям (4) и условию на бесконечности (5) удовлетворяет пучок функций

у = у0 +юМФ1 +юМ(1 -юм)ф2 , (8)

где функция У0 имеет вид (7).

Таким образом, формула (8) задает структуру решения краевой задачи (3) - (5).

3. Построение итерационного процесса

Обозначим нелинейный оператор в правой части уравнения (3) буквой Б:

Б 1 (ду 9Еу ду Жу!

^ r2 sin 0 V д0 dr dr д0 )

r2 sin(

*<*>£ - 2 *1 Еу.

В задаче (3) - (5) сделаем замену у = у 0 + u, где u -

новая неизвестная функция, а у0 - функция вида (7).

Тогда для функции u получим краевую задачу с однородными краевыми условиями

РИ, 2012, № 3

Б(у0 + u) = Bu + Uro cos 0--Uro

dr

sin 0 dEu r d0 .

Для решения задачи (9) - (11) воспользуемся методом последовательных приближений. Пусть начальное приближение u(0) задано. Например, можно взять

u(0) = 0.

Если k -е приближение u(k) построено, то новое

(k +1) -е приближение u(k+1) находим как решение линейной задачи

vE2u(k+1) = Б(у0 + u(k)) - vE2у0 вне Q , (12)

(k+1)

= 0.

dQ

du(k+1)

dn

dQ

= 0,

(13)

lim u(k+1) = 0 (14)

В соответствии с теоремой структура решения задачи (12) - (14) имеет вид

u(k+1) = ю2 ф(k+1) + 2 (1 -ю )ф(k+1) u =юМ Ф1 +юМ(1 -юМ)ф2 .

Для аппроксимации неопределенных компонент

Ф(k+1), Ф^1 воспользуемся методом Галеркина-

Петрова.

Известно [22, 23], что общее решение уравнения E^ = 0 при отсутствии в физической постановке сингулярностей может быть записано в виде

го /

у (r, 0) = Z (An rn + Бп r1-n + Cn rn+2 +

n=2 '

+Dn r3-n) Jn (cos0), (15)

где An , Бп , Cn , Dn - произвольные постоянные;

Jn(Z) - функции Гегенбауэра первого рода. Представлением (15) воспользуемся для выбора координатных последовательностей.

15

Для аппроксимации неопределенной компоненты Ф(к+1) воспользуемся функциями системы

{r-1J2(cos9), r-1J4(cos9), r-2 J3(cos9), r-2 J5(cos9),

r-n Jn+1 (cos 9), r-n Jn+3 (cos 9),...}, (16)

Поскольку функции Tn - частные решения уравнения E2u = 0 , то E2u(Nk+1) = 0 в области {со (r, 9) > M} . Функция F в этой области имеет вид

F = Bu(k)

+ U0 cos 9

5Eu(k)

dr

- U

sin 9 5Eu(k)

r 59

а для аппроксимации неопределенной компоненты

(k

ф21і+1) воспользуемся функциями системы

{j3(cos9),r J2(cos9),r2 J2(cos9),r4 J2(cos9), r3 J3(cos9), (cos9),..., rn Jn(cos9), rn+2 Jn(cos9),...} . (17)

r? J3(

Итак, функции Ф^1 и ФІУ1 представим в виде

Ф^4 и ф11"+1) = ank+1)Tn

1 1, mi ^ n n

1 n=1

ф

(k+1) и ф(k+1) =

m2

,(k+1).

Ф2,m2 S an+m1 Ln+m1 , n=1

где ть

Tm1+1,

(17).

Тогда

Ш - первые Ш1 функций системы (16), а

Ш1+Ш2

- первые m2 функций системы

N

,(k+1) и u(k+1) = N a(k+1)

N

= Zank+14

n=1

(18)

и допускает оценку

F = O

при r ^ +00 .

Таким образом, в области {со (r, 9) > М} невязку можно сделать сколь угодно малой за счет подходящего выбора M .

Для нахождения коэффициентов a(k+1), ..., а^+^ воспользуемся методом Галеркина-Петрова [24], взяв в качестве проекционной систему функций {fi}, где

fi = юМ Ті . Тогда a(k+1), ..., a^+9 найдем из условия

ортогональности невязки Rn (20) элементам f1, ..., fN проекционной последов ательности:

(Rn, fi) = 0, і = UN, (21)

причем в силу сделанных выше замечаний интегрирование в (21) при вычислении скалярных произведений можно производить только по конечной области

{0 <ю (r, 9) < М} .

где N = Ш1 + ш2 ,

2 2 Ф1 = юм т1, Фш1 = юМ ТШ1,

2

фШ1+1 = юМ (1 — юМ)тШ1+1 , •••, фN = юМ (1 — юМ)ТШ1+ш2 •

Таким образом, построенные функции фп образуют координатную последовательность.

Обозначим

F = B(^0 + u(k)) - vE2y0 . (19)

При подстановке u(Nk+1) из (18) в уравнение (12) получим невязку

N

Rn = vE2uNk+1) - F = v£ ank+^E2Фп - F . (20)

n=1

Подставляя (20) в (21), получаем, что a(k+1), ..., a(k+1)

- это решение системы линейных алгебраических уравнений:

N (k+1) 2

S a(k+1) (vE2фn, fi) = (F, fi), i = 1, 2, ..., N .(22)

П=1

Решив систему (22), получим новое приближение u(k+1) . Итерации следует прекратить, когда

u(k+1) - u(k) І <є, где є > 0 - малое число.

Предварительные результаты работы были доложены на пятнадцатой Всеукраинской (десятой международной) студенческой научной конференции по прикладной математике и информатике СНКПМИ-2012 [25].

Выводы

Изучим поведение невязки Rn (20) и функции F (19) в области {со (r, 9) > М} . Заметим, что в этой области

u(k+1) = Ф^1 и S ank+1)тп ,

n=1

u(k) = Ф«

Ш1

San

(k).

n=1

Впервые разработан численный метод расчета внешних течений вязкой несжимаемой жидкости, основанный на совместном применении методов последовательных приближений, R-функций и Галеркина-Петрова. Разработанный метод отличается от известных универсальностью (алгоритм не изменяется при изменении геометрии области) и тем, что структура решения точно учитывает как краевые условия на границе обтекаемого тела, так и условие на бесконеч-

16

РИ, 2012, № 3

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

Литература: 1. Ламб Г. Гидродинамика. М.: РХД, 2003. Т. 1. 452 с.; Т. 2. 452 с. 2. Ландау Л.Ф., Лифшиц Е.М. Теоретическая физика. В 10 т. Т. VI. Гидродинамика. М.: Физматлит, 2003. 736 с. 3. Лойцянский Л.Г. Механика жидкости и газа. М.: Дрофа, 2003. 840 с. 4. Андерсон Д., Таннехил Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. М.: Мир, 1990. Т. 1. 384 с.; Т. 2. 392 с. 5. Гущин В.А., Матюшин П.В. Математическое моделирование пространственных течений несжимаемой жидкости // Мат. моделирование. 2006. 18, № 5. С. 5-20. 6. Приходько А. А., Редчиц Д. А. Численное моделирование нестационарного течения в следе за цилиндром на основе уравнений Навье-Стокса // Прикладна гідромеханіка. 2005. 7 (79), № 1. С. 56-71. 7. Рябенький В.С., Торгашов В.А. Безытерационный способ решения неявной разностной схемы для уравнений Навье-Стокса в переменных: завихренность и функция тока // Мат. моделирование. 1996. 8, № 10. С. 100-112. 8. РоучП. Вычислительная гидродинамика. М.: Мир, 1980. 616 с. 9. Флетчер К. Вычислительные методы в динамике жидкостей. М.: Мир, 1991. Т.1. 504 с; Т. 2. 552 с. 10. Chung T. J. Computational fluid dynamics. UK: CUP, 2002. 1012 p. 11. Liu J.-G., Wienan E. Simple Finite Element Method in vorticity formulation for incompressible flow // Math. of Computation. 2003. 10, № 2. P. 1130 - 1145. 12. Pozrikidis C. Fluid dynamics: theory, computation, and numerical simulation. USA: Kluwer academic publishers, 2001. 557 p. 13. Рвачев В.Л. Теория R-функций и некоторые ее приложения. К.: Наук. думка, 1982. 552 с. 14. Колосова С.В. Применение проекционных методов и метода R-функций к решению краевых задач в бесконечных областях. Дисс. ... к.ф.-м.н.: 01.01.07 Вычислительная математика. Харьков: ХИРЭ, 1972. 85 с. 15. Колосова С.В., Сидоров М.В. Применение метода R-функций к расчету плоских течений вязкой жидкости // Вісн. ХНУ. Сер. Прикл. матем. мех. 2003. № 602. С. 61 - 67. 16. Суворова И.Г. Компьютерное моделирование осесимметричных течений жидкости в каналах сложной формы // Вестн. НТУ ХПИ. Харьков, 2004. № 31. С. 141 - 148. 17. Тевяшев А.Д., Гибкина Н.В., Сидоров М.В. Об одном подходе к математическому моделированию плоских стационарных течений вязкой несжимаемой жидкости в конечных односвязных областях // Радиоэлектроника и информатика. 2007. № 2. С. 50-57. 18. Максименко-Шейко К.В. Математическое моделирование теплообмена при движении жидкости по

каналам с винтовым типом симметрии методом R-функций // Доп. НАН України. 2005. № 9. С. 41 - 46. 19. Полянин А.Д., Зайцев В. Ф. Справочник по нелинейным уравнениям математической физики: Точные решения. М.: ФИЗМАТЛИТ, 2002. 432 с. 20. Колосова С.В., Ламтюгова С.М., Сидоров М.В. Про один метод розв’язання зовнішніх задач гідродинаміки в’язкої рідини у наближенні Стокса // Тринадцята міжнародна наукова конференція ім. акад. М. Кравчука: Матеріали конференції. Т. 2. К.: НТУУ, 2010. С. 150. 21. Стрельченко А.Й., Колосова С.В., Рвачов В.Л. Про один метод розв’язування крайових задач // Доп. АН УРСР, сер. А. № 9. 1972. С. 837 - 839.

22. Полянин А.Д. Справочник по линейным уравнениям математической физики. М.: ФИЗМАТЛИТ, 2001. 576 с.

23. Хаппель Дж., Бреннер Г. Гидродинамика при малых числах Рейнольдса. М.: Мир, 1976. - 630 с. 24. Красносельский М. А., Вайникко Г. М., Забрейко П. П., Рутицкий Я. Б., Стеценко В. Я. Приближенное решение операторных уравнений. М.: Наука, 1969. 420 с. 25. Ламтюгова С.М. Застосування методів послідовних наближень та R-функцій до розрахунку зовнішніх вісесиметричних в’ яз-ких течій // П’ ятнадцята всеукраїнська (десята міжнародна) студентська наукова конференція з прикладної математики та інформатики СНКПМІ-2012: Тези доповідей. Львів: ЛНУ, 2012. С. 229 - 231.

Поступила в редколлегию 10.09.2012

Рецензент: д-р физ.-мат. наук, проф. Колосов А.И.

Колосова Светлана Васильевна, канд. физ.-мат. наук, проф. кафедры прикладной математики ХНУРЭ. Научные интересы: математическое моделирование, численные методы математической физики. Увлечения и хобби: театр, искусство и литература. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. (057) 7021436.

Ламтюгова Светлана Николаевна, аспирантка кафедры прикладной математики ХНУРЭ, ассистент кафедры высшей математики ХНАГХ. Научные интересы: математическое моделирование, численные методы, математическая физика, теория R-функций и её приложения. Увлечения и хобби: искусство и литература. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. (057) 7021436.

Сидоров Максим Викторович, канд. физ.-мат. наук, доц. кафедры прикладной математики ХНУРЭ. Научные интересы: математическое моделирование, численные методы, математическая физика, теория R-функций и её приложения, стохастический анализ и его приложения. Увлечения и хобби: всемирная история, история искусств. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. (057) 7021436.

РИ, 2012, № 3

17