АЛГОРИТМЫ ЧИСЛЕННОГО РЕШЕНИЯ УРАВНЕНИЙ НАВЬЕ-СТОКСА ПРИ НАЛИЧИИ КАВИТАЦИИ Текст научной статьи по специальности «Физика»

Электронный журнал «Труды МАИ». Выпуск № 44

www.mai.ru/science/trudy/

УДК 532.528

Алгоритмы численного решения уравнений Навье-Стокса при наличии кавитации 1

Н.Л. Маркина

Аннотация

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

Ключевые слова:

кавитация; двухфазное течение; PISO; уравнения Навье-Стокса; численное моделирование.

Математическая модель течения жидкости при наличии кавитации

С проблемой возникновения и развития кавитации сталкиваются при рассмотрении широкого круга вопросов, связанных с течением жидкости: от проектирования насосов, турбин, клапанов, жиклеров, до исследования тока крови в сосудах. Явление кавитации в гидравлических системах может сопровождаться рядом неблагоприятных эффектов такими как шум, эрозия, вибрация, увеличение потерь энергии, уменьшение КПД и др.

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

1 Работа выполнена при поддержке РФФИ (проект № 08-08-00203-а)

алгоритмов. Кроме того, существенным фактором становится выбор метода связывания давления и плотности парожидкостной среды.

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

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

д(рти) д , 2 дыч д , ды. дР

УИт ' +—(ри2-т—)+—{р„ыу- =—) — (1)

дt дх дх ду ду дх

д(р„У) д , ду. д , 2 ду. дР УИт ' + — (ртыу-и—) + — (рУ- =—)----(2)

•л, -л т ' -л ' -л т ' -л ' -л ^ '

дt дх дх ду ду ду

дР т , д(РтЫ) , д (РтУ)

дt дх ду

= 0 (3)

а+д(аы1+д(а1 = т+-т^ (4)

дt дх ду

Рт = а Р1 + (1 -а )Ру (5)

где t - время, ы и у- компоненты вектора скорости, Р - давление, рт- эффективная

плотность, а - объемная доля жидкой фазы, т - скорость испарения, т+ - скорость конденсации пара, р1 -плотность жидкости, ру - плотность пара, и - динамическая вязкость.

В работах [4-6] источниковые члены т+ и т являются функциями давления и содержат эмпирически определяемые коэффициенты С^ и Сргоа. Также известны подходы к

вычислению коэффициентов С^ иСргоа через нахождение скоростей интерфейсов между

паром и жидкостью [8].

Для связи уравнения неразрывности и уравнения момента использован модифицированный алгоритм PISO [9], в котором дискретизация уравнений производится по методу контрольного объема на разнесенной сетке.

Вывод дискретных аналогов

Основная идея метода PISO [2] заключается в том, что для расчета давления используются два разностных уравнения для поправки поля давления, полученных из дискретных аналогов уравнений моментов и неразрывности. Такой подход связан с тем, что скорректированные первой поправкой скорости, могут не удовлетворять уравнению неразрывности, поэтому вводится второй корректор, который позволяет вычислить скорости и давления, удовлетворяющие линеаризованным уравнениям количества движения и неразрывности. Ниже приведен алгоритм получения уравнений для поправок давлений.

Используя процедуру, подробно описанную в [10], дискретный аналог уравнений (1)-(3) для нахождения скоростей и* и v* в контрольном объеме с гранями Дх и Dy можно записать в виде (предполагается, что известны поля плотности р°т и давления P0; узловые

. 1 . . . 1 точки поля и имеют индексы i + —, j, узловые точки поля v имеют индексы i, j + —, узловые

точки поля P имеют индексы i, j , где i, j е Z ):

* ДхДу п

и 1 (-рт 1 + a 3 + a 1 + a 1 + a 1 ) =

i+2j Dt iV V i-2j i+2j+1 i+2j-1

= ДДДУ Pm" 1 u" 1 +a 3 u 3 +a 1 U\ +a 1 u 1 +a 1 ,u 1 - (PP+1j - ФДУ , (6.1)

Dt i+ 2j i+2j i+ 2j i+2j i-2j i-2j i+2j+1 i+2J+1 i+ 2j-1 i+2J-1

* ДхДу n

V 1(^T" Pm 1 +a , 1 +a , 1 +a 3 +a 1) =

j+2 Д j+2 1+1 j+2 1-1 j+2 i + 2 j-2

ДхДУ n n * * * * /- -r^n t^O \ *

- Pm .. 1V. 1 +a 11V11 +a 11V11 +a 3V* 3 +a 1V* 1 - (P,+1 - P )Дх , (6.2)

» i m.. i .. i .1.1.1.1 .1.1.1.1 .. з .. j ,.i..i

Дt j+ 2 'J+ 2 i+ 2j+ 2 i+ 2j + 2 i 2j+ 2 ' 2J + 2 j+ 2 j+ 2 j 2 j 2

" ij+1 ij

j +— j +— j — j — 2 2 2 2

где коэффициенты a вычисляются по формулам:

ax++3. = Pi+1j I) + max(-^+1j,0), a1 . = D.A(| Pe. |)+max(^j.,0),

i 2j i 2j

ax 1 , = ^ 1 1 A(|Pe 1 1!) + max(-^1 1,0), a 1 , = D , ,A(\P„ , , |)+max(F 1 1,0);

i+—I+1 i+—I +— i+—н— i+—I +— i+—i-1 i+—i — i+—i — i+—i—

2 2 2 2 2 2 2 2 2 2 2 2 2 2

.. Ду .. Ду f 1 1ДХ f 1 1ДХ

D = f+1jДУ D = fjДУ D = i+1 j+2 D = i+2j - 2

i+1 j х х , ij х х , D 1 1 = , D 1 1 = ,

Xi+3 -x.+1 х.+ 1 -х-1 -2j+2 yj+1 -yj i+2j-2 у. -yj-1

2 2 2 2

FMJ = P°mi+lJul+lJAy, Fu = p°miluIlAy, F. i , i =P™+. i , iv.. i , iAx,

miJ ij

i+- J ±- i+- J ±- i+- J ±-

2 2 2 2 2 2

диффузная проводимость и интенсивность конвекции на гранях контрольного объема;

F

i+ij

F,

F , i

i+- j ±

_ P =2JL P = 2

ei+i] D ' eiJ D ' ei+iJ ±A D

DUL 2j 2 LJ. i Л

- числа Пекле; при решении в качестве

yi+ij

i+- J ±-2 2

аппроксимации конвективных членов использовалась схема со степенным законом распределения А(Ре) = тах(0,1 -0.11 Ре |5).

Для конечного объема с центром в точке р можно записать уравнения (6.1) и (6.2) в операторном виде:

uP = H(u) , - D , (VP0), +Afy (PmU)

где u.

u,л

v

V P 0

P ' Af At v"m"/p'

- вектор скорости, H(u), =

H(u), V H(v) P 0

(6.3)

- конечно-разностный оператор,

полученный при дискретизации конвективных и диффузионных потоков,

'D(u), 0 ^

D

V

0 D(v)

матрица коэффициентов, где D(u)P = Ay и D(v)P = AX,

p

AU

AV

(V P),

/ (VP), ^

(VP) P

Vv J, 0

- градиент давления в конечном объеме.

ф * *

Т.к. скорости ы и V не удовлетворяют уравнению неразрывности, то для нового

* ** **

поля давления Р рассчитываются скорректированные скорости ы и V :

** ,АхАу п .

ы 1 (-рт 1 + а 3 + а 1 + а 1 + а 1 ) =

'V А 1+2] 1+2] 1-2] 1+2]+1 1+2]-1

AxAy

Pmn i u" i. + а 3 и з + a i и r + а r u*^ + a r -(P+i. - P*)Ay, (7.i)

At ' "'i+ 2J i+2J i+ 2J i+2J i-2J i-2J i+2J+i i+2J+i i+J i+2J-i

„AxAy

V i(^T" P m i + a , i + a , i + a 3 + a i) =

J+2 At J+2 i+i J + 2 i-^J+2 4+2 J-2

AxAy

Pmn. iv i + a i.iv\ .i + a i.iv\ .i + a 3v* 3 + a iv* i - (P*+i-P*)Ax, (7.2)

» ' Ш .. i . . i .1.1.1.

At J+ 2 J+ 2 i+ 2J + 2 i+ 2J + 2 i 2J+ 2 ' 2J + 2 J+ 2 J+ 2 J 2 J 2

Вычитая из уравнений (7) уравнения (6) получим выражения для коррекции

скоростей:

и 1 = и 1 + и , = и

Ду

г+-/ г+-/ 1+2 Г 2

V 1 = V 1 + V 1 = V 1

К . 1 . , (р+1 у ру).

—] и— ] А 2 2 ^..1

Дх

1+-] 2

1]1]+ 2 2

?/+- У+- А У 2 2 Л 1 1]+-2

(Ру +1 - Р])

(8.1)

(8.2)

где

А 1 =

1+—] 2

ДхДу

Рт. 1 .+ а 3. + а 1 .+ а 1., + а 1

' т. 1 . . .з . .1. . 1 . , . 1 . ,

Д/ 1 +— / 1 +— / 1--/ 1+— /+1 1+— /-1

V ^ 2 2 2 2 2 0

. А 1 =

ГАхАу п

——Рт 1 +а , 1 +а 1 + а 3 +а 1

Д/ 1/ +— 1+1 / +— 1-1 /+— 174— 1]--

V ^ 2 ■'2 2 2 2 0

или в операторном виде:

и** = и; -Ър(УР')р . (8.3)

Связывание полей плотности и давления производится через введение аналога скорости звука Ср по формуле:

Рт* Рт0 + Рт] РУ + С оУр/ ■

О

(9)

ту ' ту ' т1] Гц р/ 1]

Уравнение неразрывности для новых полей скорости и плотности принимает вид:

Р^ т1] Р^ т1] А А * ** А * ** А * ** А * ** А /"\

-ДхДу + р 1 и 1 Ду -р 1 и 1 Ду + р 1V 1 Дх -р 1V 1 Дх = О.

* **

* **

Д/

1+— / 1+— /

2 2

1—] 1—] 2 2

1 / +— 1 / +— 2 2

(10)

1 — п— 2 2

Подстановка в уравнение неразрывности (10) выражений (9) и (8) приводит к уравнению для нахождения поправок давления Р :

р\ДхДУс +р0 Ду2

¿Л а^ ^рп^Гт. 1

Ду2

Дх

Дх

Д/

Я. 1 . л

1+-] А 2 Л 1 1 +- ] 2

+ Рт 1 + Рт" -+ Рт" -+

1 ] А , у^ А , у^ А

2' 1 .

1-2 ]

2 1

У+-2

2"-.. 1 1]-2

+Ср х тах(и 1 .0)Ду + С 1 тах(-и 1 .0)Ду + С 1 max(v х.0)Дх +

1'+-у 1'+-у 1—у 1—у у+- у+-

2 2 2 2 2 2

+Ср - тах(^* -. 0)Дх) + Р1+1] (р

Ду2

я - „ СР 1 тах(-и 1 .0)Ду) +

1+2]А 1 Р2]

2

!+- У 2

+Р-1У (-Рт0 1

Ду2

Дх2

1 - 2 У А 1 1-2 ]

Ср 1 тах(и 1 .0)Ду) + Ру+Х(-Рт - ---Ср - тах(-у -.0)Дх) +

1-2 ]

1-2 ]

1] +- А -1 2 Л 1 У+

У+-Г,

У +-г,

Лт2

+Р--Сг 1 тах(у* 1,0)Лх)= -

о

Ртп -Р

ту г ту

"■■ 1 л

1 - 2 Л 1

} " 2

2

2

Л/

ЛхЛу - р 1 и* 1 Лу -

'+-} '+-} 2 2

-Рш" 1 1 ЛУ + Рш° 1V* 1Лт РШ 1V* 1Лт)

1--} '--}

2 2

}+- У+-2 2

'}--} -

2 2

(11.1)

или в операторной форме: ЛтЛу

о

Л/

(СР)р + Л(СрРи)р-Л(РШЩУР)•=)

Ртр р

тр г~тр

Л/

ЛтЛу -Л(рри *) р. (11.2)

При выводе уравнения (11) слагаемыми вида (СрРЩ(УР ) • 8) пренебрегается. Следует отметить, что для дискретизации членов Л(СрРи*) необходимо применение

противопоточной схемы. В работе использовалась схема против потока первого порядка.

При вычитании из уравнения (7) уравнения (6) слагаемые, содержащие оператор Н, сокращаются, что приводит к необходимости ввода второго корректора для получения полей

*** *** Т~»**А *-» *-» / 1 \

и , V и Р . Аналогично уравнению (6) записывается дискретный аналог уравнений (1)-

(3) для полей и*** и V***:

*** , ЛТЛу п ч

и , (-р , + а , + а , + а , + а , ) =

. 1 Л л , ' т. 1 . . 3 . . 1 . . 1 . , .1.,''

1+2} Л/ V V !-2} 1+2 }+1 !+2}-1

ЛтЛу Л/

Ртп! и"! + а з и"з + а х и** х + а ^ и* 1 + а г - (Р+1. - Р**)Лу, (12.1)

1+^1 1+^1 1+^1 1+^1 1-^1 1-^1 1+^.+1 1+—}+1 1+-}-1 1+-}-1

*** .ЛХЛу " ч

1(^Т" Рт 1 +а , 1 +а , 1 +а 3 +а 1):

1+2 Л/ 1+2 ^^ 1 Ч+2 1-2

ЛтЛу

Рт" 1 ^ 1 + а 1 . 1 Л. 1 +а 1. 1 Л. 1 +а 3 ^'з + а 1 1 - (р1+1 - РГ)Лт . (12.2)

» /т., 1..1 .1.1.1.1 .1.1.1.1 ,.з..з ..1..1 ^ У+1 У

Л/ 1+ 2 11+ 2 1+ 21 + 2 1+ 2]+2 1 21 + 2 1 21+2 1 + 2 1+ 2 1 2 1 2

Вычитая из уравнений (12) уравнения (7) получим выражения для коррекции

скоростей:

и

1 = и\ +и и 1 = и*1 ^ ^^- 4)-d 1(р++11- р;1)

_ т ^J__т ^J__т ^J__т /I ^^ ^J__т ^ ^

2 2 2

21 Л 1

}

2

1+- } 2

V 1 = V 1 + V , = V , +

2 2 2

1-' 1 ' X (v*¡ - ^пк ) - d 1(P,'+1 - P,").

2 1+2 Л 1 11+2 1 1 '

(13.1)

(13.2)

Лу Ах

где d 1 = —, d 1 =-

1+-> Л 11+— Л 2 ^ 2 1

В операторном виде выражения (13.1) и (13.2) записываются как и7 = ир + Н(и * - и*) ; - Б ; (УР ) ; .

Коррекция полей плотности и давления производится

Рту Рту + Рт1у р* + СрУр/ .

Р] 1У

(13.3) по формулам: (14)

Р* = Р/+ Р/.

(15)

Подставляя в уравнение неразрывности

Рту Рту_ ДхДу + р**- и*** Ду - р**1 и™- Ду + р*'*1 V™ Дх -р**1 V™ Дх = 0

Д/

1+— * 1+— * 2 2

1—1 1 —] 2 2

1/ +— 1/ +— 2 2

1— 1/ -2 2

(16)

выражения для поправки плотности (14) и скоростей (13). получим уравнения для вычисления поправки давления Р'':

р(ДДуСру +Рт*+- Ду + Рт*-- аДу + Рт*/+-] ДХ + р], ]Лх +

Д/ 1+21 1+21 1 21 1 21 1]+2 1]+2 1 2 1 2

+Ср , тах(и , .0)Ду + Ср , тах(-и , .0)Ду + Ср , max(v ,.0)Дх + Ср , тах(-V ,.0)Дх) +

1'+-у 1'+_ у 1—у 1—у у+— у+- у— у—

2 Г Г Г 2 2 2 2

+Р++1 ] ( рр 1 а - Ду - Ср - тах(-и**- . 0)Ду) + р- (-рр - а - Ду-Ср - тах(и** - .0)Ду) +

1+— / 1+— / 1+— / 1+— / 1— / 1— / г— / 1— /

2 2 2 2 2 2 2 2

+рт'+1 ( рр 1 а - Дх - Ср 1тах(-т**1.0)Дх) + р-^-р - а - Дх - Ср - тах(у* -.0) Дх) =

1*+- 1*+- 1]'+- ^ — У — 11 — —

2 2 2 2 2 2 2 2

рт] рт1]'

ДхДу - (рт 1 и 1 Ду-рт 1 и 1 Ду + рт 1V 1 Дх - рт 1 1ДхУ

.......■ " 2 У-2

Д/ 1+2у 1+2у 1* 1-^] *у+~ у-~ *

1—* 1—* 2 2

У'+- У+-2 2

. Ду * / ** * ч Ду * х-™1 / ** * \

-рт 1 Xа , (и - и ) ^^-рт 1 Xа , (и - и ) +

А , 1+-* 1+-у А , — - 1 ■

, 1+-] •+-]

1 2 2 1+^ У 2

1-21

1--1 У

2 2

Дх * ^ т , ** *. Дх * ^ т , ** *.,

+ ~-рт.. 1 X а 1(V - V) ^-р». - X а - V ))

А 1 У+;

А. -

У-7

(17.1)

у-;;

или в операторной форме:

(СрР ) ; + Д(СрР"и**); -Д(р:Б(УР") •=); Д

рт; р:; ДхДу-Д(р:Ц**)

-Д(р:н(и - и*) • 8); (17.2)

Уравнение переноса (4) дискретизировалось в недивергентной форме с использованием противопоточной схемы первого порядка точности:

ак+1 -ак ак+1 -ак+1 ак+1 -ак+1 ак+1 -ак+1

а а + тах^+1, 0) ^ а + тах(-ик+1, 0)а-+ тах(^+1, 0)1 а +

Л/

хг+1 х

х Х1-1

ак+1 - ак+1 д

+тах(-^к+1,0) _ + ак Ч

у I - у 1-1 Эх

к+1

Эv +—

Эу

к+1

= ) ш+ - т- .

Чтобы предотвратить схемные осцилляции, слагаемое

ак. +1( ^

1 Эх

к+1

Эv +—

Эу

а^(-ик+1 Р

р к+1 ( и эх

к+1

к+1

заменялось эквивалентным

ту

- V,,

к+1

ЭРш

Эу

к+1 _ к+1

Рту -Р

ту

Л/

полученным из уравнения неразрывности.

Алгоритм расчета

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

1. Принимая Р0 = Р", и° = и", V0 = V", гШ = Рш" , а = а" вычислить поля скоростей и* и V* по формуле (6.3) с использованием неявной схемы.

2. Найти величину аналога скорости звука Ср по формуле Ср = С(1 - а* + е), где С и е - константы.

3. Вычислить поле поправок давления Р, решая уравнение (11.2) с использованием неявной схемы.

4. Скорректировать поля скоростей используя (8.3) и найти поле давления Р* = Р0 + Р^ .

5. Решая уравнение

а^л а^ц ** 1 а^ц ** ^^а оС^_11 ** ^^а+1 'ц

+ тах(и„, 0) —--— + тах(-и„, 0) —-- + max(vг; ,0)

Л/

+

хг+1 Х1

Х1 Х1-1

« а1 -а 1 а.. „ Эр +тах(-V* ,0) —-^+ 1-и; Эр

ЭРш

Эу

0 „ 0

Рту_ Р ту

=Л/

у1 у1 -1 Р ту

вычислить значение объемной доли жидкой фазы а**.

6. Найти поле плотности жидкости рт* = а**р1 + (1 -а**)ру.

7. Рассчитать величину аналога скорости звука Ср = С(1 - а** +е) .

у+1- у о ш+- ш,

8. Вычислить поле поправок давления Р . решая уравнение (17.2) с использованием неявной схемы.

9. Скорректировать поля скоростей по формуле (13.3) и найти поле давления

Р** = Р* + Р".

V V V

10. Из уравнения

+тафГ. 0)а+1] -а + тах(-и^. 0) +тах(1. 0) +

Д/ х1+1-х х -х-1 у*+1 -уу

+тах(-]. 0)а* 1 + 1 (-] ^

у* - у*-1 р:* дх

т др:

V ——

1

* _ п

^ ту ' ту \ + -- ^ :+- Р-

Д/

ду

найти значение объемной доли жидкой фазы а***.

11. Вычислить поле плотности жидкости р:** =а***р1 + (1 -а***)ру.

^0 х-|** * *** * *** 0 ** * ***

12. Полагая Р = Р . и = и . V = V . рт = рт . а =а . повторить итерации 2 -11.

Процесс массообмена. происходящий в областях жидкости с давлением близким к

критическому значению давления Рат. моделируется членами т + и т- уравнения (4).

представляющими скорость конденсации и испарения жидкости. В данной работе скорость массообмена вычислялась по формулам:

:+ = С;гоа тах( Р - Pcav .0)(1 -а) m_=-Cdesrmin(р-Pаv,0а (0.5рги¥)/т ' ру (0.5рги¥)/ш '

где - средняя скорость течения в невозмущенном потоке; /¥ = ^^ - характерный

временной масштаб; ЬсЬ -длина тракта; Сае!1(. Сргоа - эмпирически определяемые

коэффициенты. Подробный анализ кавитационных моделей и влияния величины коэффициентов Сае!1( и Сргоа на характеристики течения приведен в [1].

Для связывания полей плотности и давления использовалось соотношение

р: = С(1 - а + е)Р. где С - константа. величина которой не влияет на результат

вычислительного эксперимента. но выбор большого значения С может дестабилизировать расчет; при малом значении С значительно замедляется скорость сходимости итерационного процесса. В данной работе С = 2.

Результаты тестовых расчетов

Для верификации описанного алгоритма рассматривается задача течения жидкости с кинематической вязкостью 10-6 м2/с в плоском тракте длиной 23 мм. Жидкость из входного сечения шириной 32 мм поступает в узкое сечение тракта шириной 4 мм и длиной 16 мм. Так как течение имеет вертикальную ось симметрии, в качестве расчетной области рассматривается половина тракта. На стенках тракта ставятся условия прилипания, на входе в тракт задается профиль скорости и . Истечение жидкости производится в среду с атмосферным давлением, критическое давление Рсск равняется 3 КПа. Плотность пара ру

считается равной 1 кг/м3. В численных экспериментах применяется ламинарная модель течения. На рисунках 1 и 2 приведены расчетные мгновенные поля распределения гидродинамических параметров для режима суперкавитации.

10'

15'

-2

т

У

1000 900 800 700 600 500 400 300 200 100 0

104

15-

-2 -1 0 У

3.0Е+05 2.8Е+05 2.6Е+05 2.4Е+05 2.2Е+05 2.0Е+05 1.8Е+05 1.6Е+05 1.4Е+05 1.2Е+05 1.0Е+05 8.0Е+04 6.0Е+04 4.0Е+04 2.0Е+04 0.0Е+00

5-

10'

15-

-2 -1 0

_У_

а) Фотография течения [10]

б) Поле плотности, кг/м3

в) Поле давления, Па

Рис.1. Распределения гидродинамических параметров при а = 0.65 .

0

0

0

5

24 22 20 18 16 14 12 10 8 6 4 2 0 -2 -4

10-

_ И11111 п -2 -1 0

У

-20

01

5-

10-

15-

±1±

-2 -1 0 У

а) Продольные скорости и. м/с

б) Поперечные скорости V. м/с

Рис.2. Распределения скоростей при а = 0.65 .

5

На рисунке 3 показана зависимость относительной величины кавитационной зоны

Ьст = . где Ьст - абсолютная длина кавитацтонной зоны. - длина узкой части тракта. от

величины перепада давления аР в тракте. Наблюдается хорошее согласование рассчитанной и экспериментальной [11] кривых.

Рис.3. Зависимость величины кавитационной зоны от перепада давления в тракте. 1-

эксперимент; 2-расчет.

Выводы

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

Библиографический список

1. Frikha S., Coutier-Delgosha O., Astolfi J. A. Influence of the Cavitation Model on the Simulation of Cloud Cavitation on 2D Foil Section. //International Journal of Rotating Machinery. - 2008. - pp. 1-12.

2. Issa R.I. Solution of the implicitly discretised fluid flow equations by operator-splitting. //Journal of Computational Physics. - 1985. - vol. 62. - pp. 40-65.

3. Wang G., Ostoja-Starzewski M. Large eddy simulation of a sheet/cloud cavitation on a NASA0015 hydrofoil. //Applied Mathematical Modeling. - 2007. - vol. 31. - pp. 417-447.

4. Merkle C.L., Feng J.Z., Buelow P.E.O. Computational modeling of the dynamics of sheet cavitation. //Third International Symposium on Cavitation. - 1998. - pp. 307-311.

5. Kunz R.F., Boger D.A., Stinebring D.R., Chyczewski T.S., Lindau J.W., Gibeling H.J., Venkateswaran S., Govindan T.R. A preconditioned Navier-Stokes method for two-phase flows with application to cavitation prediction. //Computers and Fluids. - 2000. - vol. 29. - pp. 849-875.

6. Singhal A.K., Vaidya N. and Leonard A.D. Multi-dimensional Simulation of Cavitating Flows Using a PDF Model for Phase Change. ASME FED Meeting. Paper No. FEDSM'97-3272. - 1997.

7. Lee C., Tae-Seong R. Flow instability due to cryogenic cavitation in the downstream of orifice. //Journal of Mechanical Science and Technology. - 2009. - vol. 23. - pp. 623-649.

8. Senocak I., Shyy W. Interfacial Dynamics-Based Modeling of Turbulent Cavitating Flows, Part-1: Model development and steady-state computations. //Int. J. for Num. Methods in Fluids. - 2004. - vol. 44. - pp. 975-995.

9. Senocak I., Shyy W. Interfacial dynamics-based modelling of turbulent cavitating fows, Part-2: Time-dependent computations. //Int. J. for Num. Methods in Fluids. - 2004. - vol. 44. - pp. 9971016.

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

11. Sou A., Tomiyama A., Hosokawa S., Nigorikawa S., Maeda T. Cavitation in a two-dimensional nozzle and liquid j et atomization (LDV Measurement of l iquid velocity i n a nozzle). //JSME International Journal. - 2006. - vol. 49, No. 4. - pp. 1253-1259.

12. Сиов Б.Н. Истечение жидкости через насадки. - М.: Машиностроение, 1968. - 140 с.

Сведения об авторах

Маркина Надежда Леонидовна, аспирант Московского авиационного института (государственного технического университета).

Нижняя Первомайская ул., 50, кв.79, Москва, 105203; тел.: 8-910-472-34-39; е-mail: markina_n@list.ru