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

УДК 519.6:532.5

Н.Б. Еремеева, С.С. Куликов

Метод численного решения задачи протекания для системы уравнений неоднородной жидкости

Предлагается метод численного расчета задачи протекания вязкой неоднородной жидкости. Полученная методами конечных элементов и баланса разностная схема, аппроксимирующая систему уравнений Навъе-Стокса в приближении Вуссинеска, решается итерационным методом. I [редставлены различные варианты результатов расчета, демонстрирующие появления зон расслоения течения жидкости при малых числах Фру-да.

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

В двумерном случае, который рассматривается в настоящей работе, при численном исследовании используются уравнения На-вье-Стокса, записанные как относительно физических неременных скорость — давление [6],так и в переменных функция тока — вихрь [7]. Общим недостатком последнего метода является использование в том или ином виде граничного условия для вихря на твердой поверхности, которое отсутствует в физической постановке задачи. Наличие дополнительного итерационного процесса, связанного с указанным граничным условием для вихря, лимитирует скорость сходимости численных алгоритмов. Очевидная ограниченность методов решения в переменных функция тока — вихрь в случае исследования пространственных потоков вязкой жидкости объясняет возросший в последнее время интерес к численному решению уравнений Навье-Стокса, записанных в естествен-

ных переменных скорость — давление:

дI

= -V Р + ^ А и +

+ (и V )и 1

Re

div и = 0;

Гт2

£Р;

(И

Ф дг

+ (ииV)Р = 0.

Система уравнений (1) записана в безразмерных неременных. Здесь компоненты скорости и = (и, V) отнесены к характерной скорости и0 , координаты х, у — к характерной длине Ь, время - к Ций,

2

давление отнесено к р0и0 , где р 0 — характерная плотность; вектор £ = (0,-1);

АР = Ртах - Ртт, ГГ = и„ Д/Ар /р„ - ЧИСЛО Фруда; Re = и0Ь/V — число Рсйнольдса; V — кинематическая вязкость; р — плотность, отнесенная к р; £ - ускорение свободного падения. В случае медленных течений

(и << 1)

,

оолыноп протяженности можно опустить нелинейные инерционные члены [4]. В области

□ = {( х, у):0 < х < 1, ф (х) < у < / (х)ф (х) < / (х)}

(рис.1)

Рис. 1. Область течения

для стационарной системы (1 ) ставится следующая краевая задача:

1

и = 0,и = к (у) > 0,

I Г0 I Г! 1

и г = к2(у) > 0,

I 1 2

4 = 0,^ =Р1(у)(пРи к1 (у)>0), (2) здесь Г0 = ЛБС1)ОЕЕ; Г1 = ЛЕ — участок

втекания; Г- = СО

;

Г = Г„ II Ц II Г2.

Для построения дискретного аналога задачи (1)^(2) в области О, используется сетка

_|(Х, у = ¡ку = ]к ,к = УM,

^ = 14'^ •■-7Ц I ще Мщ

к = (Дх) -<Цх))/ N,0 < г < М,0 < ] < ы\ N - целые положительные числа. Посредством О!1 обозначим множество четырехугольников О. у ,г = 0,М - 1, у = 0, N - 1 с

вершинами (х,, уР ) , (х,+1, уг-+1,}) ,

(у1+1,^ (xг+l, у а Г,.

г+1 'У г,у+1 >

граница

^у. Каждый четырехугольник разбивается на два треугольника посредством соединения верпнти (х 1, у у )и (Х^^ уг+1]+1). Пусть Н

— пространство непрерывных функций, линейных на каждом треугольнике, а Q — множество функций, постоянных па каждом

четырехугольнике ^у. В качестве базисных

Н

«пирамиды» Кураттта., носитель которых изображен тта. рисунке 1 в виде шестиугольника. Множество этих функций

Му(х*, ут) = 8л8т Л к = а М, у, т = 0, N; и являются линейными на каждом треугольнике разбиения области пк . Каждая

функция Щк е Н может быть представлена как

М N

Щк (х, у),

г=0 ]=0

(3)

где Щ значение Жк (х, у) в узле (х ,, у у ),

Значения Ру, г =1,М, ]=0^, (р0у =р 1(у0у ),у = 0, N) будем искать из условия минимума функционала, к построению которого сейчас приступаем. Считая, что выполнено условие и = 0, имеем

у (Р ) = | (иР х + VР у = = | ирёу — урйх »

- ^(ЦД. +и+1,;А+1,; )(у+1,; -у1} ) +

(иг+1,]+Р+1 ]+1 + 4+1,у Рг+1,у )кг+1 +

+(иуР+1 +4+13+Р+ХмЩ+1 у+1]+1) ] +4,у+Рг, + -- (V иР+и + )к+ О^у* + У+щР+ь+1)к1 =

=ау Ру +ЬуР+,] +C]Р+lh]+1 +d]Рl]+1, г = 0,М-1у = 0, N-1.

М -1 N-1

Пусть I (Р ) = 2 XX ■

•=0 ]=0

Минимум этого функционала будем искать как решение системы д1 (Р)

РР]

-■0г=1, М,у = 0, N

что приводит к следующей системе алгебраических уравнений:

] +ь_1 у1- _1,] + с -1 ,]Ч 1 -1 +

+]] = 0 (4)

г=1Ми=0, N.

В системе (4) полагаем I у = 0, если выполнено одно из условий:

г<0, г>М-1,у<0,у>^1.

Приближенным решением задачи (1) - (2) назовем совокупность функций

{ик, Ук, Рк, Рк }е Н х Н х Н х

удовлетворяющих системе (4) и соотноше-

ниям

1

Я7е

(УГ, VЩк)+(Рк ,У Щк) +

+((икЧ) ик ,жк)=к ,Щк);

(ё1у ик ,1) ^ =| uкdxdy = 0.

] П]

г = 0, М -1, ] = 0, N -1,

(5)

(6)

где ик = (ик, Vк); Щк =М у (х, у), Мт (х, у)};

г,к=1М-1; ],т=1^-1;

(ик, V1) =| ик (х, y)vк (х, y)dxdy - скалярное

г

произведение в пространстве Ь2(0.к ) — функций, интегрируемых с квадратом. Систему (4)-(6) решаем итерационным методом по следующей схеме (опущен индекс к):

pn

-—(Vun+1, V W) + (Pn+1, VW) + Re

+(ГУ) un WW)=-^(gpn ,W); Fr

'+1| = Pn\ -t Ь' divun+1|

ri ri » I г

(7)

(8)

где

div un

=-1— [ (u + v )dxdy =

mesпД x У

1-1 udy - vdx =

mes D.

1

2mesQ, ^ + u«)("У) +

+ (ui+1, j+1 + ui+1, j )h+1 + +(ui ,j+1 + u+1,j+1)(y! ,j+1 - У+1, j+1) -- (uu + u,j+1)h -

- (V+1, j + V, j )h+(v, j+1 + v+1, j+1)h].

Система (7)-(8) относительно и"+1,Р"+1 решается (Р,х)-методом, подробное изложение

[]

тема (4) имеет положительно определенную

симметричную матрицу, свойства которой

[]

этой системы применялся попеременно-[ ].

приведены изображения линий уровня функции тока х,у) течений в прямоугольной области. Здесь Re=1,

к( у) = 4у(1 - у),0 < у < 1,

И2( у) = 256у(0.25 - у),0 < у < 0.25,

р 1(у) = 1 — у/ (0, у). Значения чисел Фруда указаны под рисунками. Приведенные результаты демонстрируют влияние величины плотностного числа Фруда на характер течения.

Литература

1. Yin C.S. Slraüfied flows. New York, 1980.

2. Васильев О.Ф., Квон В.И., Лыткин Ю.М., Розовский И.Л. Стратифицированные течения // Итоги науки и техники. Гидродинамика. Т. 8. М., 1975.

3. Белолииецкий В. М., Костюк В. Ю., Шокии Ю. М. Математическое моделирование стратифицированной жидкости. Новосибирск, 1991.

4. Белолииецкий В.М., Гекова С. II., Туг овиков В. Б., Шокии Ю.М. Численное моделирование задач гидроледотермики водотоков. Новосибирск, 1994.

5. Антонцев С. II., Кажихов А. В., Монахов В.II. Краевые задачи механики неоднородной жидкости. Новосибирск, 1983.

6. Белоцерковекий О. М. Численное моделирование

в механике сплошных сред. М., 1984.

7. Кускова Т. В. Численное исследование двумерных течений вязкой несжимаемой жидкости // Некоторые применения метода сеток в газовой динамике. Выи. 3. М., 1972.

8. Кобельков Г. М. О численных методах решения уравнений Навье-Стокса в цеременных скорость-давление // Вычислительные процессы и системы. Выи. 8. М., 1991.

9. Кузнков С. С., Семенов С. П. Метод численного расчета задач цротекания стратифицированной жидкости // Вычислительные технологии. Т. 4. ц12. 1995.

К). Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М,, 1978.

п. 0

ал

о ъ

О 2

20 О.ЗО

ОЧО 0.50 О,

Рис. 2. Изображение изолиний функции тока при РгЮ.1

о.8

ал

о .о

Рис. 3, Изображение изолиний функции тока при Рт=0.05

О £>

Рис. 4. Изображение изолиний функции тока при Рг=0.04