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

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

Том 15, № 5, 2010

Решение стационарной задачи о течении вязкой жидкости в канале, вызванном заданным перепадом давлений, при наличии внутренних источников

H.A. Гейдаров, Ю.Н. Захаров Кемеровский государственный университет, Россия e-mail: naznaz@mail. ru, zyn@kemsu. ru

Ю.И. Шокин

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: shokin@ict.nsc.ru

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

Ключевые слова: внутренние источники, многосвязная область, перепад давления, вязкая жидкость, уравнения Навье—Стокса.

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

Вторая постановка заключается в задании на участках втекания—вытекания не скоростей, а давления (т. е. движение жидкости в области протекания в такой постановке осуществляется за счет разности давлений). Следовательно, в этом случае задачу следует решать не в формулировке "вихрь—функция тока", а только в формулировке "скорость—давление". В работе [2] показано, что для существования и единственности решения нестационарной задачи, кроме задания давления на участках втекания-вытекания, необходимо задать лишь одну из компонент скорости так, чтобы вектор скорости на границе втекания—вытекания был ей перпендикулярен. Однако отсутствие значения второй компоненты вектора скорости не позволяет без дополнительных условий на нее построить процесс численного решения таких задач. Обычно численно решаются задачи протекания в каналах с параллельными осям координат прямыми границами и прямыми углами (см. [3, 4]). В этом случае на входе и выходе задаются

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

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

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

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

и выходе при наличии прямоугольного источника/стока:

V Ди

д(и2) д(иь) др

(1)

vДv

д(иь) д(ь2) др

(2)

дх ду ду

(3)

р \лБ = Р1(у^ р\сю = P2(y),

и \бс = 0, илд = 0, \(и, п)\ = |и\ на границе АВ, \(и, п)\ = \и\ па границе ОБ,

(4)

(5)

(6) (7)

и |Г1 = и1 (у)•

и|г2 =

и2(х), и |гз = из (у),

и|Г4 =

и4(х),

'ийв

= 0,

(8) (9)

вед(А> В' С' Б')

где и = (и, V) — вектор скорости, и = и(х,у), V = V(х,у); р = р(х,у) — давление; V > 0 — коэффициент кинематической вязкости; & — область решения (рис. 1); Г = Г1 и Г2 и Г3 и Г4 — граница источ ника Примем, что об ласти & и&8 имеют форму прямоугольников, стенки которых расположены параллельно осям координат, как показано на рис. 1 (в ином случае задача усложнится лишь технически, но не изменится содержательно).

Равенства (6), (7) означают, что вектор скорости и перпендикулярен границам АВ и ОБ, а условие (9) означает, что суммарный поток скорости через границу внутреннего источника равен нулю. Подобное условие используется в работе [7] как достаточное условие существования глобального решения стационарной задачи о течении вязкой жидкости в канале с заданным вектором скорости на границе.

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

Далее будем считать, что задача (1)-(9) имеет хотя бы одно решение.

2. Для численного решения полученной дифференциальной задачи (1)-(9) построим в ограниченной области & прямоугольную и в общем случае неравномерную сетку согласованную с границей Г = АВ и ВО и ОБ и АБ и Г8.

На узлах (х^, у?) Е & введем сеточные функции и?, vij1 р? и аппроксимируем систему дифференциальных уравнений (1) разностной схемой

V Ан и?

VАнVij -

22

кхл+1

и

з+1 — из -^¿з -1 — р^13

+ К

к

У?+1

+ к

УЗ

vi+1j щ+13 — vi-1j ui-1j vij+1 — vij -1

к

Р?+1

xi+1 + кхЛ

Рз-1

^+1 + кя

к

У3+1

к

"/ • 1.У — "/ 1.У '"/.У • 1 — '"/.У 1 _ ^

^+1 + кя

к

У3+1

+ к

+ к ■ + кУ3

, (ъ,3) Е 31

У3+1

+ к

УЗ

(10)

УЗ

У

В

А

а г2

^ о, г3

г4

С

В

О

Рис. 1. Область течения. Канал АБС Б с внутренним источником О.

где </1 — множество индексов (г, 3), отвечающих внутренним точкам. Разностный оператор Д^ аппроксимирует оператор Лапласа:

«¿+1.7 «у Му «г—у «¿7+1 «у «у — 1

^ _-1—Гл--1--1—и,-'

2 2

Вместе с аппроксимацией краевых условий (4)-(9) получим разностную задачу относительно вектора и = («у, ^ ), (г, 3) £ /ь Введем также множество индексов /2, которые соответствуют тем граничным узлам, где нет краевых условий на компоненты «у вектора решения.

Число уравнений полученной системы меньше числа неизвестных в силу того, что, во-первых, отсутствуют условия на горизонтальную компоненту вектора скорости на границах протекания (на вертикальную компоненту скорости поставлено условие (6), (7), означающее V = 0 на границах АВ и СБ), во-вторых, нет условий па давление на части границ области.

Для предположения существования и единственности решения полученной разностной задачи необходимо каким-либо образом на входе и выходе определить значение горизонтальной составляющей вектора скорости, а также величину давления на твердых стенках, но так, чтобы эти условия были следствиями условий (1)-(9). В противном случае задача может не иметь решения.

Для вычисления «у на множестве узлов /2 рассмотрим следующее. Так как на АВ и СБ эти переменные не заданы, то их можно определить из уравнения (1) путем аппроксимации его внутрь области решения, заменяя производные односторонними разностями первого или второго порядка аппроксимации. На твердых стенках ВС и АБ, а также на границе источника Г краевые условия на давление отсутствуют. Считая, что уравнение неразрывности (3) выполнено и на данных границах области, также аппроксимируем его внутрь области решения, В итоге аппроксимации замыкают разностную схему (10), и в этом случае число уравнений системы будет равно числу неизвестных.

Таким образом, полученную разностную задачу можно записать в виде системы билинейных уравнений размерности т

А(и)М = /, (11)

где = А1(м)г; + А2г;. Здесь А1(м) — матрица, полученная в результате аппрокси-

мации нелинейной части исходной дифференциальной задачи, А1 — числовая матрица,

/

лучаемая из краевых условий. Элементы матрицы А1 линейно зависят от компонент вектора М = (м1, ..., мт)т.

Поскольку указанная задача (11) включает в себя аппроксимацию уравнения неразрывности (3), то при сходимости какого-либо применяемого для ее решения итерационного метода получаемый вектор решения будет удовлетворять разностной схеме, аппроксимирующей уравнение сИум = 0,

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

«п+1/2 = - Гга+1[А(и>га - /], (12)

мп+1 = мп+1/2 + ага+1гга, п = 0,1, 2,..., (13)

где гга = (г1,..., гт)т — некоторый ненулевой вектор; и0 — произвольное начальное приближение из области определения оператора А; гп+1 _ либо числовой итерационный параметр, либо квадратная диагональная матрица; ага+1 — диагональная матрица итерационных параметров, В работе [8] показано: если элементы матрицы ага+1 выбирать из условия минимума последовательности невязок, то норма невязки ||гга|| = ||А(ига)ига — / || монотонно убывает и при некоторых ограничениях на оператор А(и) стремится к нулю. Это обстоятельство позволяет находить хотя бы одно из решений разностной задачи (подробнее о методе и способах выбора ага+1 см, в [8]),

Рассмотрим первый шаг схемы (12), (13) в случае, если тга+1 является диагональной матрицей. Для этого выпишем систему (11) в виде

ап(м)м1 + Й12(м)м2 + ... + а1т(«)«т = /1 «21(и)и1 + а22(и)«2 + ... + а2т(и)«т = /2

ат1(и)«1 + ат2(и)«2 + ... +

Далее будем полагать, что агг(и) = 0 для всех и, Пусть элементы матрицы этой системы зависят от и линейно, т. е, представляют собой линейные комбинации компонент вектора и

%(и) = X]Х' + Й°>,

причем

А1(и) = (Йу (и)), (и) = «у,

А 2 = ), А(и) = А1(и) + А2,

где икр — одна го компонент вектора и = (и1,... , ит)т, 1 < к?, < ш, 1 < г,^, р < т, , «0,- _ некоторые числа.

Запишем матрицу А(и) системы (11) как сумму двух матриц А(и) = 5(и) + Т(и),

где

а и 1

5(и) = (и)) , «г, (и) = Е , > ^

I

Т(и) = (и)) , ¿г,"(и) = Е , < ^

г

Йг, (и) = (и) + ¿г,(и).

Выпишем разложение матриц 5(и) и Т(и) на нижнетреугольную, диагональную и верхнетреугольную части соответственно:

5 (и) = $ь(и) + 5Л(и) + (и),

Т (и) = Ть(и) + ТЛ(и) + Ти (и).

Введем следующие матрицы:

£(и) = 5^(и) + Ть(и), Л(и) = 5л (и) + Тл(и), и (и) = (и) + Ти (и).

Получим

А(и) = Ь(и) + Л(и) + и (и) = = Бь (и) + Ть(й) + Бл(и) + ТА(и) + Би (и) + Ти (и).

Тогда первый шаг схемы (12), (13), где тп+1 = сИа§(тп+1,..., тт+1), можно записать в виде

[БА(ип) + ТА(ип+1/2)] ип+1/2 = [БА(ип) + Тл(ип+1/2)] ип+к/(2т)-

-т к+1п+(к+1)/(2ш) - - т П+1/2

1 п+1 ' к+1 ек+1 ... I п+1'т ет,

к = 1,... ,т, п = 0, 1, 2,..., [Бл(ип) + Тл(ип+1/2)} ип+к/(2т) = [Бл(ип) + Тл(ип+1/2)} ип-

-т 1 _п+1/(2т) е - 2 гп+2/(2т) - - к _п+к/(2т)р тп+1'1 е1 тп+1' 2 е2 ... тп+1'к ек,

где и0 — заданный произвольный вектор, ек = (0,..., 0, 1, 0,..., 0)т — вектор с едини-

„ 7 п+к/(2т) 1

цеи на месте с номером к гк _ к-я компонента вектора невязки

гп+к/(2т) = [Бь(ип) + Ть(ип+1/2) + Бл(ип)+

+Тл(ип+1/2) + Би(ип) + Ти(ип+1/2)]ип+к/(2т) - ¡.

Обозначим

Ь(ип, ип+1/2) = Бь(ип) + Ть(ип+1/2), Л(ип, ип+1/2) = Бл(ип) + Тл(ип+1/2), и (ип,ип+1/2) = Би (ип) + Ти (ип+1/2), А(ип,и п+1/2) = Ь(ип,ип+1/2) + Л(ип,ип+1/2) + и (ип,ип+1/2).

Зависимость матриц от двух параметров показывает, что некоторые компоненты данных матриц содержат ип, а другие ип+1/^. Как следует из обозначений, Ь(ип,ип) = Ь(ип), Л(ип,ип) = Л(ип), и(ип,ип) = и(ип), А(ип, ип) = А(ип). Тогда

Л(ип,ип+1/2)ип+1/2 = Л(ип, ип+1/2)ип - тп+1(1(ип,ип+1/2)ип+1/2+ + (Л(ип, ип+1/2) + и(ип,ип+1/2))ип - f).

Отсюда получаем следующий вид (12) в случае, когда тп+1 является квадратной диагональной матрицей итерационных параметров:

т"+\ (Л(ип, ип +1/2) + тЬ(ип, ип+1/2)) (ип +1/2 - ип) +

+ (А(ип, ип+1/2)ип - f) = 0. ^

т 1+1

||гп+1/(2т) тп2+1 — го условия минимума нормы ||гп+2/(2т) ||, тк+1 — из условия минимума нормы ||гп+к/(2т)||. Как показано в [8], эти параметры вычисляются как корпи кубического уравнения по явным формулам Кардано,

Таким образом, схема (12), записанная в виде (14), является методом последовательной верхней релаксации решения систем билинейных уравнений с покомпонентной вариационной оптимизацией параметров. Подробнее о методе см, [9].

4. В данном разделе приведены результаты расчетов течения вязкой несжимаемой жидкости, вызванного перепадом давления, в прямоугольном канале при наличии прямоугольного источника/стока как с условием па скорость (7), так и без пего.

Отметим, что кроме схемы (10), построенной дня "совмещенной" сетки (в одном и том же узле сетки рассчитываются все сеточные функции), использовалась и схема па "разнесенной" сетке (см, [10]), когда каждая из переменных и, V, Р рассчитывается на своем наборе узлов.

На рис, 2-5 приводятся линии тока (оттенками серого цвета показано распределение функции давления),

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

1) возможно ли эффективно применить предлагаемый в работе итерационный метод (12) дня решения поставленной задачи;

2) правомерно ли практическое использование краевого условия (6), (7) па входе и выходе капана;

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

Расчеты проводились для различных значений вязкости V, шага сетки Н, перепада давления. Итерационный процесс останавливайся при выполнении условия

11 У

невязка па итерационном шаге п. На-

< 10 4, где г°

чапьное приближение и0 = 0, на левой границе выполнено условие (6) (т. е. V = 0),

Рис. 2. Течение с условием на касательную составляющую вектора скорости V \оо = 0| на правой границе канала; р1 = 0.1 р2 = 0; условия на границах источника: и\гх = -0.05, и\г2 = -0.05, V \гз = 0.1 V\г4 = 0.1 V\гх = V \г2 = и\гз = и \ Г4 =0

Рис. 3. Течение без условия на касательную составляющую вектора скорости на правой границе; р 1 =0.1 Р2 = 0; условия на границах источника: и \ гх = -0.05, и \ г2 = -0.05, V \ Гз = 0.1 V \ г4 = 0.1 V \ п = V \ Г2 = и \ Гз = и \ г 4 =0

Результаты расчетов показали, что итерационный процесс (12), (13) сходится дня всех используемых значений Н, V и для всех решаемых задач. Отметим, что закон сохранения массы в каждом расчете выполняется с точностью порядка 0,01 %. При этом максимум модуля элементов вектора невязки находился в диапазоне 10-4-10-3, а максимум модуля элементов вектора дивергенции не превышал 10-5,

При решении (1)-(9) (т. е, вместе с условием (7) о равенстве пуню касательной составляющей вектора скорости) жидкость вытекает из капана иод прямым углом, без учета влияния источника па течение вблизи границы, даже если он находится достаточно близко от нее, В этом случае линии тока могут резко менять направление (рис, 2), Если же не задавать условие (7), то для замыкания разностной задачи (1)-(6), (8), (9) аппроксимируем исходную систему па правой границе внутрь области решения так же, как это было сделано для определения компонент и р^- вектора решения на соответствующих границах (см, разде,:: 2), При этом получается течение (рис, 3), в котором линии тока подходят к правой границе капана иод разными углами.

Таким образом, если не задавать на правой границе условие V = 0, то жидкость вытекает из канала, но нашему мнению, более естественно.

Поскольку жидкость может вытекать из источника с произвольной заданной скоростью независимо от заданного па границах области решения давления, то течение на левой границе может быть направлено как внутрь, так и вне капана таким образом, что жидкость будет вытекать через границу ЛБ, Пример такого расчета показан на рис, 4 (жидкость вытекает через часть левой границы). При этом заранее предположить, куда будет направлен вектор скорости на границе, довольно трудно. Чтобы в этом убедиться, достаточно немного увеличить длину канала — в данном случае жидкость начнет втекать в область решения через левую границу (рис, 5, а).

Если дня канала, показанного на рис, 5, а, перепад давления уменьшить до значения 0,04, то жидкость снова будет вытекать через левую границу канала (рис, 5,6"), Таким образом, способ замыкания системы уравнений с помощью аппроксимации уравнений

Рис. 4. Жидкость вытекает через часть левой границы ЛБСБ; р\ = 0.1 р2 = 0; условия на границах источника: и |Г1 = -0.12, и |г2 = -0.12, V |Гз = 0.01 V |г4 = 0.01 V |г1 = V |г2 = и |гз = и |г4 =0

Рис. 5. Течение в удлиненном канале; p2 = 0; условия на границах источника: u |гх = -0.12, u |г2 = -0.12, v |г3 = 0.01 v |г4 = 0.01 v |Г1 = v |г2 = u |Гз = u |г4 = 0; а — p1 = 0.1, б — pi = 0.04

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

Из приведенных расчетов следуют выводы:

1) итерационным методом (12), (13) можно достаточно успешно решать задачу о стационарном течении вязкой несжимаемой жидкости в каналах, когда источником движения является перепад давлений на концах отрезка, и при наличии источников/стоков внутри капана;

2) использование аппроксимации самой системы уравнений на входе и выходе канала позволяет находить решения задачи (1)-(9), (11)—(13), не задавая значений одной из компонент вектора скорости.

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

[1| Роуч. П. Вычислительная гидродинамика. М.: Мир, 1980. 612 с.

[2| Рагулип В.В. К задаче о протекании вязкой жидкости сквозь ограниченную область при заданном перепаде давления или напора /7 Динамика сплошной среды: Сб. науч. тр. Ин-та гидродинамики СО АН СССР. Новосибирск, 1976. Выи. 27. С. 78 92.

[3] Кузнецов Б.Г., Мошкип Н.П., Смагулов Ш. Численное исследование течения вязкой несжимаемой жидкости в каналах при заданных перепадах давлений /7 Численные методы динамики вязкой жидкости: Сб. науч. тр. НТПМ СО АН СССР. Новосибирск, 1983. С. 203 207.

[4] Moshkix N., Yambaxgwi D. Steady viscous incompressible flow driven by a pressure difference in a planar T-junction channel // Intern. J. Comput. Fluid Dvn. 2009. Vol. 23, No. 3. P 259 270.

[51 Метод минимальных невязок решения одного класса нелинейных уравнений / Ю.Н. Захаров, Е.Ф. Егорова, М.А. Толстых, Ю.Н. Шокин. Красноярск, 1991 (Пренр. ВЦ СО АН СССР; № 9-91).

[6] Geidarov N.A., Zakharov Y.N. Stability of solution of stationary viscous incompressible fluid flow produced by a given pressure drop problem /7 Proc. of Intern. Conf. "Mathematical and Informational Technologies MIT-20097 Kopaonik, Serbia, 2009. P. 140 143.

[7] Leray J. Etude de diverses équations intégrais non lineaires et de quelques problèmes que pose l'hydrodynamique // J. Math. Pures Appl. 1933. Ser. 9, 12. P. 1-82.

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

[9] Гейдаров Н.А., Захаров Ю.Н. О треугольных методах решения систем линейных и нелинейных уравнений с вариационной оптимизацией параметров // Вестник КемГУ. 2009. № 2(38)'. С. 34-39.

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

Поступила в редакцию 16 июля 2010 г.