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

Вестник КемГУ

№ 4

2008

Математика

УДК 519.6

РЕШЕНИЕ ЗАДАЧИ О ПРОТЕКАНИИ ВЯЗКОЙ ЖИДКОСТИ ПРИ ЗАДАННОМ ПЕРЕПАДЕ ДАВЛЕНИЯ

Н. А. Гейдаров, Ю. Н. Захаров

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

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

Вторая постановка заключается в задании на участках втекания-вытекания давления, а не скоростей. Т. е. движение жидкости в области протекания осуществляется за счет разности давлений. На участках протекания границы отсутствуют условия для компонент скорости, что затрудняет процесс решения. В работах [2], [3] показано, что для существования и единственности решения нестационарной задачи, кроме задания давления на участках втекания-вытекания границы, необходимо задать и значения одной из компонент скорости. Отсутствие значения второй компоненты вектора скорости не позволяет без дополнительных условий на эту компоненту построить численный процесс решения таких задач. Обычно численно решаются задачи протекания в каналах с параллельными прямыми границами и прямыми углами. В этом случае задаются «естественные» краевые условия на скорости вида ди ду п

— = 0, — = 0, у = 0 и т. п. (см. [2, 3]), которые по-

дх ду

зволяют построить алгоритмически замкнутый численный метод решения. Однако такое «простое» определение скоростей на границе не всегда возможно, т.к. оно может не вытекать из характера течения на границе. Например, если на КЬ и МИ (см. рис. 1) задать вертикальную составляющую скорости равную нулю (в соответствии с [1]), то получающееся решение не является «естественным», т. к. вблизи точек М и N жидкость будет вытекать под прямым углом к границе. Т. о., задание всех компонент вектора скорости на участках протекания позволяет получить приемлемые решения лишь в областях сравнительно простой формы: с параллельными твердыми стенками и перпендикулярными им границам входа и выхода.

Существует несколько подходов к решению стационарной системы уравнений Навье-Стокса. Наиболее широко используется метод стациониро-вания: решение исходной задачи находится как предел решения нестационарной задачи. Второй подход состоит в решении каким-либо итерационным методом непосредственно системы нелинейных алгебраических уравнений, аппроксимирующей исходную стационарную задачу. В первом случае обычно используют неявные схемы решения, которые требуют задания краевых условий для всех граничных компонент вектора неизвестных. Еще одной особенностью такого подхода является требование соле-ноидальности поля скоростей на каждом временном слое, удовлетворение которого также является технически сложной задачей. Во втором случае не возникает трудностей как с выполнением условия со-леноидальности поля скоростей (соответствующее уравнение является частью решаемой системы), так и с устойчивостью применяемого итерационного метода.

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

1. Рассмотрим двумерную задачу о протекании вязкой однородной несжимаемой жидкости в криволинейном канале с заданным перепадом давления на входе и выходе.

8(и2) 8(иу~) 8р

дх р р х

д( иу) д( у2 ) 8р

р р р у

Р Р у /■ N ^

-----1----= 0, (х, у) е О.

р р

vAu — vAу —

р

КЬ'

МИ

P2, Р1 > P2,

и ,,, = 0, и = 0,

\ЬМ ’ 1км ’

(1)

(2.1)

(2.2)

(2.3)

|(«| кь ,^^ = |и\ кь |, |(^|ми,п^ = |и\ми|,

где и = ( и , у) - вектор скорости, р - давление. V > 0 - коэффициент кинематической вязкости.

| Вестник КемГУ № 4 І008 Математика

и = и(x, у) , у = у(x, у), О - конечный криволинейный канал, являющийся областью решения (см. рис. 1).

м

уАи и , и

п^, + п.

К+і + Кі

Пхі+1 + Пхі

_ У+і іиі+іі -У-і Іиі-1 і _ У,і+1 -Уі _ Ріі+1 - Ріі-1 (3)

Пхі+1 + Пхі

Пхі +1 + Пхі Пу]+1 + Пу]

Куі+1 + Куі Куі+1 + Куі

■■ 0, (і, і) є J,

Ы- , ■ - Ы- Ы - - и. , •

і+1 і і] і і-1]

иї+1- иї иї- иї-1

А ниї

к

к.

к

V]+1

к

К.

п..

Полученная система да нелинейных алгебраических уравнений является билинейной:

Л(й)й = / , (4)

где Л(и)у = Л1 (и)у + Л2у .

Л1(и ) - матрица, элементы которой зависят от компонент вектора и , Л2 - числовая матрица.

Считаем, что (4) имеет хотя бы одно решение.

3. Для решения системы (4) использовался явный

N

Отметим, что равенство (2.3) означает, что вектор скорости и перпендикулярен границам КЬ и МЫ.

Наряду с задачей (1), (2.1), (2.2), (2.3) рассматривалась задача (1), (2.1), (2.2), т. е. без условия на касательную составляющую вектора скорости на участках протекания.

2. Для решения полученной дифференциальной задачи (1), (2) введем в ограниченной области О прямоугольную и в общем случае неравномерную сетку Оп, согласующуюся с границей Г .

На узлах (хі, у.) є О введем сеточные функции и.., V. , р. и аппроксимируем систему диф-

У У и

ференциальных уравнений (1) разностной схемой: и \ . - и и , . и .V.. , - и.. ,У. , р , - р ,

і +1 ] і-1 ] і] + 1 ії + 1 і]-1 і]-1 -Гі + 1 і ■Гі-1 і

(5)

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

и+1/2 = и — ти+1[Л(ив К — /]

ип+1 = ип+1/2 +ап+1 хп, п = 0,1,2,..., где 2п = (21,...., )Г - некоторый вектор, и° - про-

извольное начальное приближение из области определения оператора Л , тп+1 - итерационный параметр, ап+1 - диагональная матрица итерационных параметров.

Опишем способ нахождения матрицы параметров ап+1.

Пусть ап+1 - квадратная матрица с т ненулевыми элементами а1к , 1, у = 1,2,..., да

произвольные целые числа от 1 до т. (5) может быть записано в виде:

_уи+1/і + /(р-1) +

I

а

(и+1) [

р _1, і,...,т,

(6)

п + 1/І _ п+1/І

где у( р-1) _ и

,,п+1/і „ п+1/і

У(0) _и

Р-1

-Иа1

і _1

, - вектор с одной ненулевой

к -ой компонентой - Щ _ (0,...,0,гк^, 0,..., 0)г . Введем обозначение:

п+1/і

Ч і)

а

(п+1) +

1,к,-

-к

Vя+1/і

Лі-1)

где Ь - множество индексов (1 , у), отвечающих внутренним точкам. Разностный оператор Аь аппроксимирует оператор Лапласа:

+а

(Я + 1) п

у к, ) пк.

) - /, і = 1,і,..., т.

„п+1/І _ п+1

Очевидно, что Г(т) _ г

_ А(ип

') - /

- невязка схемы (3.3). Далее:

„п + 1/і „п + 1/і

г(0 _ г(.--1)

а

(п+1) ТгЧК ) I ґгг(п+1ПІ Т7(.К )

К, + (а,, ) ^іп , (7)

I к 1 п

п+1/ І

Вместе с аппроксимацией краевых условий (2) где мы получим разностную задачу относительно вектора и = и, уу, ру), 1, ] е./ размерности да . Для

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

) _ А(у”+ І, ^ , --+ І

КІ} _ А111, г” ).

Из (7) имеем:

II г(п;1/і IIі = || г"+1/І IIі +

'(і-1)

I («і(к+1})х €] _

_И К

п+1/І IIі -ф (а("+1))

(і -1) II п+1\^ /?

/_ 1,і,

Є1^ ) _ і(гп+1/і,к

()) (;-1) 5 -* 1п

Здесь:

■ ) _ І(г(

вп _И К!,] IIі +і(Г:”+]1)/і, }),

^з(„к) _і(Кіпк),Кі(лк)), ) _||К

(к,) ||І І п II

+

| Вестник КемГУ № 4 2008 Математика

Так какв^' > 0 , то график зависимости

имеет вид, показанный на

и п+1/2 м2 а(п+1)

11 г(і) 11 от акк.

рис. 2.

\У

’.'її 41

J~l.il

/

ИЙ»

£*£■

Лп + 1/2)

не превосходит г-'і^2'1 . Выберем а

„(п+1' 1,к,

г(п+1/2) . Тогда аУпк+1)

уравнения:

I -1 вук' _ 0

ч _і

(9)

„*(п+1)

по явным формулам Кардано. Пусть а^ - реше-

ние (7), отвечающее глобальному минимуму (8). Подставив этот корень в (8), получим:

II гп+1/2 ||2 =|| гп+1/2 ||2 —ф(0 = 2

II 'а) N II'(,■ -1) II ^п+1 на

„п+1/2 м2 г(.--і) 11 ,

о2 _1 -Ф(і) /II гп+1/І ||2 Н(і) 1 ^п+і 1 II '(і-1) II ’

где ф « _ф п+а+і)) > 0. Таким образом,

II її 1/2||<|К-1)/2||< ... <||гс

(10)

п+1/І п-1/І

_ г

а

і ? і :

;( п+1)

Фп^, і _ 1, і,. ., т была наибольшей. При этом, в

)(і' і _

■п + 1 5 ‘

случае, например

к1 _ к2.

11 _ 1, і,..., т, мы будем изменять оба раза одну и ту

и *#п+1

же лі - компоненту вектора и .

Рассуждая аналогично, найдем параметр Тп+1 первого шага схемы, при выводе вместо набора век-

,,п+1/І п+1/І ~п

торов Уі используем и , а вместо е -

п+1 / 2 т-*

г . В итоге ^п+1 находится как решение уравнения IчТп+іГЧп _0, где в1 п _ 2(гn+1/І,^ '

ч_1

віп _||К1п||2 +І(Гk;1/2, К2п ),

_ 2(К,, Кп),

в4п _І|Кіп IIі, а

К1п _ А(11/2, ги+1/2) + А (гп 1 /2, ип+1/2),

К _ А, (г

2 п IV

п+1/2 и+1/2

).

Таким образом, у полинома (8) или два, или один минимум. При этом глобальный минимум величины

из условия минимума нормы вектора невязки

находится как решение

^|| Г(,-1) И-... -Ц-0 ..............

1 = 1,2,..., да, п = 0,1,2,..

и, в силу ограниченности || Г0 ||,

Пт || гп ||=|| гда ||< да.

п^да

Следует отметить, что = 0 только в случае а это возможно, если

С) = 2(г:п—1)/2, ^ = 0. (11)

Следовательно, для того чтобы последовательность {|| гп ||} была убывающей, номера

I, к, 1 = 1,2,..., да, в котором стоит элемент

Подробнее о методе см. [4].

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

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

Наряду с постановкой граничных условий, показанной выше, рассматривалась задача (1), (2), но без условия (2.3). Для замыкания разностной задачи (3),

(2.1), (2.2) мы аппроксимируем исходную систему на границах внутрь области решения.

На первом этапе были изучены течения жидкости в каналах наиболее простой формы: с параллельными или перпендикулярными стенками. Далее были проведены численные эксперименты в каналах «сложной» формы (в качестве такого канала выбран канал-диффузор) и влияние наличия или отсутствия условия (2.3) на получаемое решение.

Во всех расчетах давление р1 = 0.1, р2 = 0 . Параметр вязкости V = 0.01.

Кроме схемы (3), построенной для «совмещенной» сетки (в одном и том же узле сетки рассчитываются все сеточные функции), использовалась и схема, на «разнесенной» сетке (см. [5]), когда каждая из переменных и , у и р рассчитывается на своем наборе узлов.

необходимо выбирать так, чтобы величина

независимо от

Рис. 3. Течение в Т-образном канале. Задача (3), (2.1), (2.2). Без (2.3):

| Вестник КемГУ № 4 2008 Математика

р _ 0.1, р _ 0, р _ 0, у_ 0.01

г верх ' г лев * г прав '

2 -

1 — Ц^ЕЕ

7/. і . . . г 'Л . і . . . Г Л .

X

Рис. 4. Течение в Т-образном канале.

Задача (3), (2.1), (2.2). Без (2.3):

р _ 1, р _ 0, р _ 0, у_ 0.01

г верх ' г лев ' г прав '

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

(2.1), (2.2) (без (2.3)) дает картину течения, совпадающую с течениями, приведенными в [2, 6]. При увеличении перепада давления, вызывающего движение в Т-образном канале, появляются вихри на верхних стенках (см. рис. 3, 4), что также соответствует результатам указанных работ. Отличие от решения соответствующей задачи (2.1), (2.2), (2.3) менее 1 %.

0.8

>“

0 6 0.4 0.2

0 0.5 1 1.5

X

Рис. 5. Течение в канале с траншеей.

Задача (3), (2.1), (2.2). Без (2.3):

рвход _ 0.1, реыгход _ Л У_ 001

Полученное решение краевой задачи (2.1), (2.2), без использования (2.3), в канале-траншее (рис. 5) согласуется с решениями аналогичной задачи в работе [3], где в качестве граничных условий задан вектор скорости на границе.

При малом угле канала-диффузора течения для разных постановок краевых условий качественно совпадают, однако значения вертикальной компоненты вектора скорости для задачи (2.1), (2.2) отличны от 0 на границе вытекания, что можно считать «естественным» для данной задачи.

Рис. 6. Течение в канале-диффузоре с продолжением.

Линии тока. Задача (2.1), (2.2), без (2.3):

рход = О-1, рвыход = ^ ^ °.01

Рис. 7. Течение в канале-диффузоре. Линии тока. Задача (2.1), (2.2), без (2.3):

рвход = О.и рвыход = ^ ^ °.01

В случае течения в расширяющемся канале со стенками, параллельными вблизи выхода, получаем течение, не зависящее от наличия или отсутствия условия (2.3) (рис. 8). В [7] показано течение в канале-диффузоре при заданном поле скоростей на входе в канал. Показанный в этой работе отрыв от верхней стенки обнаруживается и при решении задачи (1), (2.1), (2.2) (см. рис. 7). Однако задание на правой границе нормальной составляющей скорости приводит к течению, отличному от приведенного в

Рис. 8. Течение в Т-образном канале с расширением. Линии тока. Задача (2.1), (2.2), без (2.3):

рверх _ 1, р _ 0.5, р _ 0, у_ 0.01

± верх ’ Г лев ’ г прав >

|| Вестник КемГУ

Дополнив канал с параллельными и перпендикулярными стенками каналом с расширением (см. рис. 9) обнаружим отрывное течение аналогичное показанному на рис. 7.

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

Литература

1. Рагулин, В. В. О гидродинамическом моделировании мембранного фильтрования вязкой жидкости /В. В. Рагулин // Вестник КемГУ. - Кемерово, 2005. - Вып. 4(24). - С. 133 - 137.

2. Кузнецов, Б. Г. Численное исследование течения вязкой несжимаемой жидкости в каналах при заданных перепадах давлений / Б. Г. Кузнецов, Н. П. Мошкин, Ш. Смагулов // Численные методы динамики вязкой жидкости. - Новосибирск, 1983. -С. 203 - 207.

Математика

3. Захаров, Ю. Н. Градиентные итерационные методы решения задач гидродинамики / Ю. Н. Захаров. - Новосибирск: Наука, 2004. - 239 с.

4. Захаров, Ю. Н. Итерационная схема минимальных невязок решения стационарной системы уравнений Навье-Стокса / Ю. Н. Захаров, О. Н. На-горнова // Проблемы динамики вязкой жидкости. -Новосибирск, 1985. - С. 156 - 159.

5. Патанкар, С. Численные методы решения задач тепломеханики и динамики жидкости / С. Патанкар. - М.: Энергоатомиздат, 1984. - 149 с.

6. Дорфман, А. Л. Численное решение задачи о течении вязкой несжимаемой жидкости в разветвляющемся канале при задании перепадов давления между ответвлениями / А. Л. Дорфман // Численные методы механики сплошной среды. - Новосибирск, 1979. - Т. 10. - № 7. - С. 47 - 53.

7. Карякин, В. Е. Расчет течений вязкой жидкости в плоских каналах произвольной формы / В. Е. Карякин, Ю. Е. Карякин // Численные методы механики сплошной среды. - Новосибирск, 1986. -Т. 17. - № 5. - С. 91 - 100.

№ 4

2008