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

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

Том 16, № 4, 2011

Численные алгоритмы моделирования кавитационных течений вязкой жидкости*

Л .В. Панов, Д.В. Чирков, С.Г. Чёрный Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: [email protected]

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

Ключевые слова: численное моделирование, кавитационные течения жидкости, баротропная модель.

Введение

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

Первый подход будем называть методом выделения границы жидкость —пар. Метод состоит в нахождении границы, разделяющей области жидкости и пара и в решении уравнений движения несжимаемой жидкости в области, занятой жидкостью [1, 2]. Метод выделения границы жидкость —пар неплохо подходит для описания пленочной кавитации вблизи тела, однако имеет ряд недостатков, среди которых — необходимость перестроения сетки на каждом шаге, сложность распространения на задачи моделирования трехмерных течений [2] и течений с отрывом паровых каверн.

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

* Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00475).

объемной (или массовой) доли одной из фаз, В баротропных моделях [3-5] система Навье — Стокса сжимаемой смеси замыкается с помощью уравнения состояния

Р = pip). (!)

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

В последнее десятилетие широко применяются модели с уравнением переноса фазы (УПФ), В них для замыкания системы Навье —Стокса используется дополнительное уравнение переноса объемной (или массовой) доли жидкой (или паровой) фазы, а плотность смеси восстанавливается по формуле

Р = aLpL + (1 - aL)pv, (2)

где aL — объемная доля жидкой фазы, pL — плотность жидкости, pv — плотность пара. Уравнение переноса имеет вид

= »

где источники m+,m- описывают соответственно процесс конденсации пара и испарения жидкости. Таким образом, модели с УПФ учитывают неравновесность процессов испарения и конденсации и влияние инерционных сил на динамику паровой каверны, В литературе предложено несколько подобных моделей [6-9], различающихся лишь видом источпиковых слагаемых m+ и m-, В [10] приводится подробный сравнительный анализ моделей с УПФ и делается вывод, что все модели дают качественно схожие результаты, Недостатки моделей с УПФ — наличие эмпирических констант в выражениях для m+ и m- и существенное возрастание источникового члена при pv/pL << 1, что является причиной неустойчивости расчетов,

В большинстве практических задач область кавитации мала по сравнению со всей рассматриваемой областью течения жидкости. Поэтому для решения уравнений кавитационных течений адаптируют численные методы, хорошо зарекомендовавшие себя при расчете течений несжимаемой жидкости, такие как SIMPLE [6, 8, 9], либо основанные на методе искусственной сжимаемости [7],

В настоящей работе на основе численного алгоритма из [11], базирующегося на неявном методе искусственной сжимаемости, предложены численные алгоритмы расчета пространственных кавитационных течений для баротропной модели и для модели с УПФ, За счет совместного решения уравнений движения и УПФ, специальной явно-неявной аппроксимации псточнпкового члена, а также введения релаксации удалось

pv

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

1. Основные уравнения 1.1. Уравнения движения

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

^+сИу(ру) = 0, (4)

др\

+ сИу(ру 0 V) + \7р = сИу(т) + pf. (5)

2

Здесь р — плотность смеси, кг/м3; V — скорость, м/с; Ь — время, с; р = р + -рк, р —

3

давление, Па; к — кинетическая энергия турбулентных пульсаций, м2/с2; ¥ — вектор

т^ — тензор вязких напряжений Тц = {р + Рт)

дпг дпЛ 2 дпк

~ -Ог-

дху ' дхг) 3 %3 дхк

Рт Р

вязкости смеси, вычисляемый по формуле

р = аьрь + {1 - аь)ру,

в которой рь,ру — динамические коэффициенты вязкости жидкости и насыщенного пара соответственно. При Т =17 °С для воды р^ = 10-3 Па • с; для насыщенного пара ру = 10-5 П а •с.

Коэффициент турбулентной вязкости вычисляется по формуле рт = С^рк2/е, где параметры к и е рассчитываются по стандартной к — е-модели турбулентности с логарифмическими пристеночными функциями,

1.2. Модель кавитации

Для замыкания системы уравнений (4) - (5) необходимо задать уравнение для опреде-

р

1.2.1. Варотропная модель

В баротропных моделях предполагается, что плотность смеси может быть найдена по

р{ р)

использовать функцию, имеющую ¿"-образный график с максимальным наклоном в

р = р у

др _ 1

где Ст;п скорость звука смеси, В [5] рекомендуется брать Стп = 1 ^ 2 м/с, В настоящей работе предлагается следующий баротропный закон:

[ ру' ф< —п/2, р — Ру

Р(Р) = { (рь-ру)/28тф+(рь + ру)/2, -тг/2 < ф < тг/2, ф= у РУ (6)

[ рь, ф>п/2. СтЛрЬ — ру)/2

Конкретный вид этого закона определяется параметром Стп (рис, 1),

Рис. 1. Зависимость баротропного закона (6) от C,

1.2.2. Модели с УПФ

Более полными являются модели, основанные па уравнении переноса массовой или объемной доли одной из фаз, В литературе встречается около 10 различных моделей, записанных дня массовой или объемной доли жидкости или пара 16—101, Для всех этих моделей плотность смеси р восстанавливается по одной из формул:

р = рь аь + ру ау, аь + ау = 1 (7)

или

- = - + -, Д + /у = 1, (8)

р рь ру

где аь — объемная доля жидкости; ау — объемная доля пара; ¡ь — массовая доля жидкости; ¡у — массовая доля пара. Объемные и массовые доли, как легко видеть, связаны следующим образом:

, Р , Р

®ь = }ь —, «у =/у—• рь ру

Эти модели путем несложных преобразований могут быть записаны в универсальной форме (3) для аь и различаться только видом источниковых членов т+ и т-, отвечающих за конденсацию и испарение, В настоящей работе используются четыре наиболее распространенные модели, источниковые члены т+ и т- которых представ-

Дх

~йоо

ная длина, и,, — характерная скорость на большом удалении от тела). Заметим, что во всех моделях выполняются неравенства

т+ = т+(р, аь) > 0, т- = т-(р, аь) < 0.

лены в таблице (здесь too = r*dp — характерный временной масштаб, Lxap — характер-

Эмпирические коистанты С7рГОс1 и имеют смысл соответственно скорости продуцирования жидкости (конденсации) и скорости разрушения жидкости (испарения).

Источниковые члены для моделей с УПФ

Модель т+ т Константы

Модель 1 [6] Сргоатах[р - ру, 0](1 - аь)рь С^тт[0,р -ру]р\аь ^ргО (1 = 80 Cdest = 1

¿со(р^/2) ру^Рьи^/2)

Модель 2 [7] СрГОа(1 - аь)а2ьРУ ¿со рь САещаьтт[0,р - ру]ру ipLU^0/2)toopL Срго а = 100 Cdest = Ю0

Модель 3 [8] 1 ^ тт \2р-ру] 2 (1 -аь)ру (-/ргск1'^оо Рь о 3 рь Р 1 ^ ГТ \2р-рууаьрь ^ев^оо РЬРУ „ 3 Рь Р Срго а =0.137 =0.274

Модель 4 [9] рьтах[р — ру, 0](1 — аь) р^тт[0,р — ру]аь _

tooi.PL - Ру№,п - Уу,п)2 pvtooipL ~ РУ№,п ~ Уу,п)2

В модели 4 [9] VI — интерфейеа, Уу,п — скорость паровой фазы по нор-

мали к интерфейсу, = и • п, п = Для нестационарной кавитации скорость

интерфейса вычисляется в процессе счета, для стационарной постановки в работе [9]

Уу,п ~ —Уь,п

предлагается формула V/ п =--, где Уьп — скорость жидкой фазы в направ-

Ру

лении нормали к интерфейсу. Для ее вычисления предлагается формула Уь,п = / ■ УУ,п где / = 0.9, Заметим, что модель 1 [6] переходит в модель 4 [9], если произвести замену

^ргоё 1 ^^^ 1

->■ 7-ттт7-т^^т, —77^"7 ->■

рьи1 /2 (рь - ру)(Ут,п - Уу,п)2 рьи1 /2 (рь - ру)(Ут,п - Уу,п)2

2. Численный метод

В этом разделе вначале рассмотрим численный метод решения уравнений движения (4)-(5) и модели кавитации с УПФ (3), а затем численный метод для баротропной модели (4)-(5) и (1),

2.1. Численный метод для моделей с УПФ 2.1.1. Метод искусственной сжимаемости

Используя уравнение переноса (3) и соотношение (2), уравнение неразрывности (4) преобразуем к виду

сИу(У) = (т+ + т~) ( —---— ]. (9)

\Рь РУ )

Решение уравнений (9), (5) и (3) основано на методе искусственной сжимаемости, согласно которому в каждое из этих уравнений добавляется производная неизвестной функции (р, или а^) по псевдовремени т и полученная система решается методом установления по т:

+ = + (10)

рдт дху \рь ру

дрпг дрпг дрпгщ

дт

дЬ

дхп

др дхг

дт

г]

дхп

+ р!г, г =1, 2, 3,

(11)

/аь\ др даь

даь даьпп 1 . +

1 = — (га+ га ).

дт дЬ дх] рь

(12)

Заметим, что изначально на каждой итерации совместно решались четыре уравнения (10)-(11), а затем решалось уравнение (12) (без слагаемого ( — ) Однако

\рр) дт

такой численный алгоритм имел низкую скорость сходимости. Существенно улучшить сходимость позволило совместное решение (10) - (12), Первое слагаемое в уравнении переноса аь (12) гарантирует совместность (непротиворечивость) уравнений (10) и (12) при расчете предельного случая течения жидкости (аь = 1, т+ = 0 т_ = 0), при котором (10) и (12) совпадают [7].

Систему уравнений (10)-(12) можно переписать в векторной форме

М

дт дЬ дх1

( 1/р

аь \ р/3

н,

( р \

рп\

, Q = рп2

рпз

аь /

(13)

(

впг

рщп + 8ир — ти

рП2Пг + б2гр — т2г рПзПг + 5згр — тзг \ аьпг

н

в

\

,Рь

1

рь

1А / +

— (га + га )

ру

рк

р/2 р/з

{т+ + т )

И* = сИщ{0,1,1,1,1).

/

В (13) можно перейти к простейшим переменным Q = {р,п1,п2,п3,аь)

т

дт дЬ Л дхг

г=1 г

н,

(14)

Р

_1

М

Щ дЦ

/ 1

р

Оь \ р/3

\

р {рь — ру )пг

р {рь — ру )п2 р {рь — ру )пз 1

/

1

1

1

1

2.1.2. Дискретизация неявным методом конечных объемов

Дискретизация (14) осуществляется согласно методике, предложенной в работе [11]. Для построения консервативной разностной схемы система (14) представляется в форме интегральных законов сохранения

V V

К ■ Ж

ШУ,,

(15)

дV

V

где К = (¥1, Р3), дУ — замкнутая поверхность произвольного фиксированного объема V, ¿Б = п ■ ¿Б = (¿Б1, ¿Б2, ¿Б3) — элемент поверхности Б, умноженный на единичную внешнюю нормаль п к ней. Вектор потоков представляется как сумма невязкого и вязкого потоков

К ■ ¿Б = (Кт + Ктз) ■ ¿Б,

Кг

/ ви1 ви2 виз \ / 0 0 0 \

рщ2 + р ри1 и2 ри1и3 т11 т12 т13

ри2и1 ри22 + р ри2и3 , Кт8 = _ т21 т22 т23

ри3и1 ри3и2 риз2 + р т31 т32 т33

V аьи1 аьЩ аьиз ) \ 0 0 0 /

(16)

(17)

Введем следующие обозначения для векторов ((, (( и И, усредненных по ячейке (г,,;, к) с объем ом У^к:

Qгjk

1

V

^к

(¿V, ((ук

1

V

У

((¿У, И

1

]к

ук

ук

V

У

ШУ,

]к

ук

V

]к

и отнесем их к центру ячеики.

Пусть и, в — номера слоя соответственно по времени и псевдовремени. Аппроксимируем интегральное уравнение (15) полностью неявной по псевдовремени разностной схемой

(дп+1).+ 1 _ ((п+1)8

+И

Р-1((Цп+1У)

г3(((п+1)8+1 _ 4((п

Дт

+

п- 1

2Дг

У

ук

(ИИБ^1)

П+1\ 8 + 1

(18)

где Дт — шаг по псевдовремени, ДЬ — шаг по времени. При сходимости итераций по псевдовремени получим полностью неявную схему, аппроксимирующую исходные уравнения (9), (5) и (3) со вторым порядком по времени ¿. Порядок аппроксимации по пространству определяется способом вычисления потоков правой части ИИБ, которая имеет вид

ИИБ = _((КБ) +1/2_(К^Б)-1/2+(К^)>Н1/2_(К^8)^1/2+(К^8)^1/2_(К^)к.1/2)+ИУ, (19)

где (К■ Б) ++1/2, (К■ Б^/2, (К■ Б)к+1/2 представляют собой разностные потоки (невязкие и вязкие) через грани (г + 1/2,;, к), (%;] + 1/2, к), (г,,, к + 1/2) ячейки с номером (г;],к) и объемом У^к. Способ аппроксимации и алгоритм вычисления вязких потоков описан в [11] и заимствован из этой работы. Рассмотрим подробно аппроксимацию невязких потоков.

2.1.3. Аппроксимация невязких потоков

Для вычисления невязких потоков будем использовать схему \П'Ж'!.. В [12] показано, что при решении предобуеловленных уравнений для повышения точности решения следует модифицировать диссипативиый член. Модификация состоит в замене матрицы |А| на матрицу Р-1|РА|

(Кт • 3)т+1/2

К^Яь) + К^Яь)) 8т+1/2 - Рт1+1/2|РА|т+1/2(ЧД

(Яп - Яь) , (20)

РА

/ 0

РА(Я) = Р

<9(Кт(Я) • 8) Щ

врБх:

врБу

— и + и\Бх и\Бу р

Б

— и2Бх II + и2Бч Р

Б,

врБ, ицБх

щБ,

0 0

0

р 0

щБх щБу и + щБг 0

0 0 0 и /

и = V • Я = и1Бх + и2Бу + щБх. При нахождении |РА| = (РА)+ — (РА)- применяется расщепление по знакам собствен-

ных значении

1

(РА) = ИБ Ь, Б^Ша^А* А± А±А±), А* =-(Аг ± |Аг|).

РА

л1,2 = и, лз,4 = и± \/и2 + Р\Б\2, л5 = и.

РА Бх = 0

И

0 0 врС врс 0

—Б, — Бу и(и + с) + вБх и(и - с) + вБх 0

0 Бх ь(и + с)+ вБу ь(и - с) + вБу 0

Бх 0 w(U + с)+ вБ, w(U — с) + вБ, 0

0 0 0 0 1

ь

Б, и — Uw+вБ,

Бу (Uw+вБ,)

(Бхи+Буу)Ц+(Бх2 + Бу2)/3 Х

БхС2

Б, (иь+вБу)

БхС2

рБхС2 с2

Бу и — |Я|2ь иь+вБу (Бхи+Б^)и +(Бх2 + Б, 2)в

рБхс2 с2 Бх с2 Бх с2

и-с Бх Бу Б2

2/3 рс2 2 с2 2 с2 2 с2

и + с Бх Бу &

2/3 рс2 2 с2 2 с2 2с2

0 0 0 0

0

При Бх = 0 необходимо модифицировать первые два столбца матрицы И и первые две строки матрицы Ь, Однако в окончательном алгоритме И и Ь в явном виде не используются, вместо них применяются формулы рационального вычисления потоков, о чем будет сказано ниже,

В формуле (20) (РА)т+1/2 = РА(((ь + Ця)/2), а векторы ((я и ((ь определяются при помощи интерполяции повышенного порядка

(ь = (т + [(1 _ е)Дт-1/2( + (1 + 9)Дт+1/2(] /4,

Ця = (т+1 _ [(1 _ 0)Дт+1/2( + (1 + 9)Дт+3/2(,\ /4,

что при 9 = 1/3 дает третий порядок аппроксимации невязких потоков на гладких решениях.

При вычислении потоков по формуле (20) требуется перемножение матриц И Ю±, Ь, отнимающее значительную часть времени счета, В работе [11] при построении схемы для уравнений Навье — Стокса несжимаемой жидкости получены формулы рационального вычисления потоков через грани ячейки, не требующие дорогостоящих перемножений матриц. Использование этих формул позволило сократить время счета на 30% и упростить программную реализацию. Оказывается, в случае системы с переменной плотностью также удается получить формулы рационального вычисления потоков в (20), Так, в случае и > 0 поток (20) равен

/

(ЮП ■ Б)т+1 /2

виь N

Рьиьиь + РьБх + Ь(ия _ иь) РьУьиь + РьБу + Ь(Уя _ Уь) рьШьиь + РьБг + Ь('Шя _ Ыь)

\

(аь)ь иь

/

/ _вс \

р (и(и _ с) + вБх) р (у (и _ с) + вБу) р (ы(и _ с) + вБг) \ _аьс )

Ь = \(рп ~ РьШп - иь), а = ± (^2{РКр Рь) +(и- с)(ия - и^ .

Значения переменных, у которых индексы ^и Ь не указаны, вычисляются как среднее арифметическое < = (<ь + <я)/2.

2.1.4• Явно-неявная аппроксимация источникового члена

И

И

К1 К1+ К1- 0

рЬ 0 0 р/1

р/2 = 0 + 0 + р/2

р/з 0 0 р/з

V Ка / V к ^ \К ) \ 0 )

ч И+ И- Н

где введены обозначения

К

±

+

-т

рь

К1±

в

,Рь

ру,

±

т

1

Источник Н3+1 можно аппроксимировать по-разному, в общем случае в виде

Н3+1 = И5 + О3(Я3+1 — Я3), (21)

где полностью явная аппроксимация соответствует О = 0, а полностью неявная — выбору

( д{к1 + Ю др

О

^(Н+ + Н~) <9(3

д{К1 + К) \ дат,

\

д(К + К) др

0

д(ь+

ь-)

дат

(22)

/

Известно [11], что для увеличения запаса устойчивости численного алгоритма все слагаемые в источнике, имеющие отрицательные коэффициенты при искомых переменных, нужно аппроксимировать неявно, В случае нелинейного источникового члена необходимо анализировать знак частной производной. Численные эксперименты с моделями 1 [6] и 2 [7] показали, что применительно к специфике нелинейных источниковых членов моделей кавитации справедливы следующие два правила явно-неявной аппроксимации источника, позволяющие значительно повысить устойчивость численного алгоритма (см, ниже рис, 6):

1) отрицательные частные производные источника, находящиеся на главной диагонали матрицы С, необходимо аппроксимировать неявно;

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

Применяя эти правила при аппроксимации источника модели 1 [6] с учетом того, что для нее

дЬ+ дЬ- , „ дЬ+ , „ дЬ-

< 0,

< 0,

< 0,

< 0,

др др дат дат

получим, что для данной модели в (21) в качестве матрицы О5 достаточно взять

ь-) а

О

др

дат

По этим же правилам аппроксимируются источииковые члены для остальных моделей с УПФ.

2.1.5. Реализация

С использованием подходов [11] схема (18) линеаризуется по методу Ньютона, затем производится Ы1-факторизация неявного оператора и схема разрешается за два дробных шага бегущим счетом.

Численный алгоритм модели с УПФ плохо сходится при малых значениях плотности пара ру. Для улучшения свойств сходимости проводится релаксация источникового Н1 ат

(И+)3 = 71 • (Н+)3 + (1 — 71) • (И+)-1, а8ь = 72 • аь3 + (1 — 72) • аЬ-1 с коэффициентами 71 = 0.001 72 = 0.2,

2.2. Численный метод для баротропной модели

При построении численного алгоритма для баротропной модели используются те же методы и подходы, что и для модели с УПФ, Отметим основные особенности, имеющие место в этом случае.

2.2.1. Метод искусственной сжимаемости. Дискретизация

Согласно методу искусственной сжимаемости уравнения движения (4), (5) и баротропной модели (1) преобразуются к виду

др др дври, дт дЬ дх,

0,

(23)

дриг дриг дриги, др дт, дт дЬ дх, дхг дх

+ = «=1,2,3,

р = р(р).

(24)

(25)

Систему уравнений (23) - (25) перепишем в векторной форме и перейдем к простейшим переменным (( = (р,и,у,ы)т

дт дЬ дхг

г=1 г

(26)

1 ( р \ ( вриг \

Р-1 = р р , ( = ри1 ри2 Е- = , Е г рщиг + 5ир _ ти рщиг + 52гр _ т2г ,И

\

р

\ риз /

0 р/1 р/2 V р/3 )

Дискретизуем систему (26) неявным методом конечных объемов. Формулы (15), (18) распространяются па случай баротропной модели с заменой матрицы И* на единичную матрицу:

р-1^ I цёУ + I цёУ + / к • = у ШУ,

V V дУ V

р-1тп+1у)

(цп+1у+1 _ 3(ЦП+1У+1 - Щп + С^"1

ДТ + 2АЬ

Угзк = (ИИБп+1Г+1.

Общий вид правой части (19) не изменится. Невязкий поток в данном случае отличается от невязкого потока (17) только в первом столбце:

К*

ври1 ри12 + р ри2и1 \ ризщ

ври2 ри1и2 ри22 + р ризи2

вриз \ ри1и3 рщиз риз2 + р )

(27)

2.2.2. Аппроксимация невязких потоков

Для баротропной модели невязкие потоки также будем аппроксимировать по схеме МиЯСЬ (20), Заметим, что при вычислении матрицы Якоби РА потока Кт • Я, где Кт имеет вид (27), учтена зависимость р = р(р) баротроиного закона. Предварительные расчеты показали, что это позволило существенно повысить устойчивость и точность

численного решения по сравнению с вариантом, где при дифференцировании Кт • Я р

РА

( вси

РА

Сии + Бх

р

С11уЛ - Бу

р

Сию4 -Б2

р

врБх и + иБх врБу иБу врБ2 \ иБ2

уБх и + уБу уБ2

п)Бх юБу и + юБ2 у

Эта матрица обладает действительными собственными значениями

= и + ир±(1, (1 = у/и2 + Щ + /ЗБ2, ир = /ЗСи/ 2, С

Л

и, Л,

др др

Были получены матрицы правых и левых собственных векторов с учетом зависимости Р = р(р) баротропного закона

И

вр(ир + д)

( 0 0 вр(ир - д)

-Б2 -Бу и(ир + и + д)+ рБх и(ир + и - д) + рБх

0 -Бу у(ир + и + д)+ вБу у(ир + и - д)+ рБу

\-Бх 0 т(ир + и + д)+ вБ? т(ир + и - д) + вБ?

/ Б2и~ Я 2гу иги+рБ2 Бу(ии)+рБ2) (Бхи+Буу)и+(Бх2 + Бу2)р\

рБхс2 с2 Бхс2 Бхс2

Буи-\$\2ь Цу+рБу {Бхи+Бгю)и+{Бх2 + Б2)(3 Б^иу+рБу)

рБхС2 с2 Бх С2 БхС2

и{(1-ир)-с2 Бх(с1 - ир) Бу(<1 - ир) БМ-ир)

2/Зфс2 2 ¿с2 2 ¿с2 2 ¿с2

и(с1 + ир) + с2 Бх(с1 + 11р) Бу{(1+ир) БМ+ир)

V 2/Зфс2 2 ¿с2 2 ¿с2 Ыс2 )

При расчете невязких потоков по схеме (20) необходимо вычислять комплексы Р-1КЮ-Ь, Р-1КЮ+Ь путем перемножения матриц. Однако явное перемножение матриц удается исключить, рассмотрев четыре следующих случая:

и > 0, вБ2 > 2иир; и > 0, вБ2 < 2иир; и < 0, вБ2 > 2иир; и < 0, вБ2 < 2иир.

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

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

Изложенные алгоритмы построены для численного моделирования трехмерных кавита-циоппых течений в турбомашипах, по в данной работе приводится их тестирование па модельной задаче стационарного обтекания затупленного сферой цилиндра турбулентным потоком вязкой жидкости с Яе = 1.36 • 105, Для этой задачи известны экспериментальные данные по распределению коэффициента Ср вдоль поверхности обтекаемого цилиндра |13|

Р - Рс

Сп

Рь иЛ2'

(28)

На рис, 2 представлены расчетная область, построенная сетка размером 120 х 60 и граничные условия па твердой стенке, па большом удалении от цилиндра и па выходе. Давление на выходной границе рс определялось через кавитационное число

а

Р с

ру

Рь иС/2

(29)

которое в эксперименте варьировалось от 0.5 (слабая кавитация) до 0.3 (интенсивная кавитация),

Если дополнительно не указано, то расчет проводился на сетке 120 х 60 с параметрами а = 0.4, ру = 10 △т = 0.001 и в качестве модели транспортного уравнения использовалась модель 1 [6]; в баротропной модели принималось Ст;п = 0.2, Расчеты проводились в два этапа: первый состоя,:: из 1000 итераций расчета течения несжимаемой жидкости и необходим дня получения начального приближения ноля гидродинамических величин, второй — расчет течения сжимаемой смеси — проводился до сходимости итераций,

и= 1

V =

X

0.5

-0.5

- 0.2

- 0.3

- 0.4

- 0.5

- 1.0 Эксперимент

\ \ 1

\

V \

\\\ \

%ч /

■

1 I

0.5 1 1.5

2.5

Рис. 2. Расчетная область, сетка и граничные условия

Рис. 3. Влияние параметра Стп в баротропной модели на распределение Ср вдоль поверхности цилиндра

3.1. Выбор параметра в баротропной модели

Представляет интерес выяснить, как влияет параметр Ст;п баротропного закона (6) на точность решения, В работе [5] предлагается задавать Ст;п = 1 — 2, На рис, 3 приведены распределения Ср вдоль поверхности цилиндра в зависимости от Ст;п баротропного закона (6), Видно, что большие значения Ст;п (в частноети Ст;п = 1) дают неправильное распределение Ср; при Ст;п ~ 0.2 — 0.3 результаты расчетов уже достаточно близки к эксперименту, поэтому следует принять Ст;п ~ 0.2 — 0.3,

3.2. Влияние способа явно-неявной аппроксимации источника на сходимость

Дня подтверждения правил явно-неявной аппроксимации нелинейных иеточпиковых членов (см, раздел 2,1,4) проведено сравнение различных аппроксимаций источника дня модели 1 |6|, Данные рис, 4 подверждают первое правило о том, что отрицательные частные производные от источника на главной диагонали матрицы G следует аппроксимировать неявно. Рисунок 5 подтверждает второе правило — производные па побочной диагонали можно аппроксимировать произвольно, сходимость при этом не меняется,

3.3. Сравнение сходимости численных алгоритмов

В алгоритмах дня моделей с УПФ используется релаксация ноля источника и ноля аь (см, раздел 2,1,5), Влияние релаксации на скорость сходимости показано на рис, 6, Можно заключить, что использование сильной нижней релаксации с коэффициентами 71 = 0.001 72 = 0.2 замедляет сходимость. Однако без релаксации не удается проводить расчеты с малыми плотностями пара ру < 50, а ее применение позволяет снизить этот параметр до ру = 3,

На рис, 7 представлено сравнение сходимости численных алгоритмов дня баротропной модели и моделей с УПФ, Видно, что скорость сходимости алгоритма дня баротроп-

Рис. 4. Учет компонент матрицы О: 1 — со- Рис. 5. Учет компонент матрицы О, содержа-держит диагональные элементы, 2 нулевая щей диаххжальные члены: члены на побочной

диах'онали матрицы присутствуют (1), нулевые (2)

Рис. 6. Сравнение историй еходимоетей чис-ленших) а.;п'оритма для модели 1 с УПФ с релаксацией источника и поля аь (1) и без нее (2) (ру = 50)

Рис. 7. Сравнение историй еходимоетей численных а;п'оритмов для баротропной модели и моделей с УПФ: 1 модель 1 [6], 2 модель 2 [7], 3 модель 3 [8], 4 модель 4 [9], 5 баротропная модель

ной модели выше таковой дня моделей 1, 3, 4 |6, 8, 9|, так как в последних используется релаксация. Модель 2 |7| сходится и без релаксации, поэтому ее скорость сходимости близка к таковой дня баротропной модели. Заметим также, что модель 4 |9| в строгом смысле не сходится, однако невязка устанавливается на уровне 10-3, а этой точности вполне достаточно дня получения решения.

3.4. Сравнение с экспериментом

На рис. 8 представлены результаты сравнения Ср, полученной экспериментально и с применением численных алгоритмов па основе баротропной модели и модели 1 с УПФ. Видно, что численный алгоритм дня модели с УПФ имеет очень близкий к эксперименту профиль Ср, численные результаты на основе баротропной модели немного хуже. Из рисунка также следует, что с увеличением в расчетной области объема пара (т. е. с уменьшением а) точность алгоритмов снижается.

Необходимо отметить, что построенный алгоритм неустойчив при ру < 1 и тем самым не позволяет провести расчет с реальной плотностью пара ру ~ 0.01 кг/м3. Поэтому представляет интерес сравнить, как влияет выбор этой константы па решение. Из рис. 9, где изображены поля аь (а) и давления р (б), видно, что форма каверны и размеры вихря существенно зависят от ру. Заметим, что эти параметры при ру =

100 качественно похожи на полученные в [9], а при малой плотности пара (ру = 10)

Ср

ру

Ср ру = 10 ру = 100

уменьшение ру в модели с транспортным уравнением (до ру = 3) не дает заметного улучшения. Таким образом, выбор ру < 10 в приложении к данной задаче оправдан.

Рис. 8. Сравнение величин Ср, полученных экспериментально и с применением численных алгоритмов, для баротропной модели (а) и для модели 1 [6] (б) при разных а

Рис. 9. Влияние константы ру в модели 1 с УПФ [6] на распределение аь (а) и давления р {б) в расчетной области (сверху вниз — ру = 100, 50, 10)

3.5. Сравнение различных моделей с уравнением переноса фазы

На рис. 11 приведено распределение величины Ср, полученное по каждой из моделей таблицы (см. с. 100). Видно, что результаты, полученные но моделям 1 и 4, очень близки к эксперименту, но модели 3 — несколько хуже, а модель 2 дает исиравилыюс распре-

С

и Саев^ что профиль Ср будет очень близок к экспериментальным точкам. Например,

С = 5 С = 200

Модель 1 [6] Модель 2 [7] Модель 3 [8] Модель 4 [9] Эксперимент

Рис. 10. Влияние константы ру в модели 1 с Рис. 11. Влияние моделей с УПФ на распре-

УПФ [61 на Cp

деление Cp по поверхности цилиндра

Рис. 12. Влияние моделей транспортного уравнения на распределение переменной аь в расчетной области

представлено распределение аь в расчетной области, видно, что в целом формы и размеры каверны и вихря дня моделей 1, 3 и 4 похожи, а в соответствии с расчетами по модели 2 каверна вообще пе образовалась, что объясняется неверным в данном случае СС

с подобранными константами) дают достаточно хорошее совпадение с экспериментом.

Заключение

В работе предложены экономичные численные методы расчета трехмерных кавитаци-оппых течений вязкой жидкости па основе баротроппой модели и моделей с уравнением переноса объемной доли жидкой фазы. Построенные алгоритмы являются развитием метода искусственной сжимаемости |11|, хорошо зарекомендовавшего себя при моделировании трехмерных стационарных течений несжимаемой жидкости. Предложен способ рационального вычисления певязких потоков по MUSCL-схеме, позволяющий па 25% сократить время расчета. В алгоритме дня моделей с УПФ за счет применения

релаксации и явно-неявной аппроксимации источпикового члена частично решена проблема сходимости при малых плотностях пара pv, Построенный алгоритм сходится при Pv/pl > 0.003, На модельной задаче показано, что близкие к эксперименту результаты получаются уже при pv = 0.01pL, дальнейшее уменьшение плотности пара на решение не влияет. Исследование влияния на решение вида зависимости p = p(p) в баротроп-ной модели позволило определить значение параметра Cmin = 0.2, На примере задачи обтекания затупленного цилиндра проведено сравнение баротропной модели и четырех моделей с УПФ, Выявлено, что все модели с УПФ дают примерно одинаковые результаты и в целом несколько точнее баротропной модели, В настоящее время модели с УПФ применяются авторами для изучения кавитационных течений воды в рабочем колесе и в отсасывающей трубе гидравлической турбины.

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

[1] Рождественский В.В. Кавитация. Л.: Судостроение, 1977. 247 с.

[2] HlRSCHl R. Centrifugal pump performance drop due to leading edge cavitation: Numerical predictions compared with model tests //J. Fluids Eng. 1998. Vol. 120. P. 705-711.

[31 Hosangadi А., Ани j a V., Arunajatesan S. A generalized compressible cavitation model // CAV 2001: Fourth Intern. Svmp. on Cavitation. California Institute of Technology. Pasadena, CA USA, 2001.

[4] Delannoy Y., Kueny J.L. Two phase flow approach in unsteady cavitation modeling // ASME. Cavitation and Multi-phase Flow Forum. 1990. Vol. 109. P. 153-159.

[5] Coutier-Delgosha O., Morel P., Fortes-Patella R., Reboud JL. Numerical simulation of turbopump inducer cavitating behavior // Intern. J. of Rotating Machinery. 2005. Vol. 2005, iss. 2. P. 135-142.

[61 Singhal A.K., Vaidya N., Leonard A.D. Multi-dimensional simulation of cavitating flows using a pdf model for phase change // Proc. of ASME Fluids Eng. Division Summer Meeting. 1997. P. 1-8.

[71 Kunz R.F., boger D.A., Stinebring D.A. et al A preconditioned Xavicr Stokes method for two-phase flows with application to cavitation prediction // Comput. k, Fluids. 2000. Vol. 29. P. 849-875.

[81 Athavale M.M., Li H.Y., Jiang Yu, Singhal A.K. Application of the full cavitation model to pumps and inducers // Intern. J. of Rotating Machinery. 2002. Vol. 8, No. 1. P. 45-56.

[91 Senocak I., Shyy W. Evaluation of cavitation models for Xavicr Stokes computations // Proc. of ASME Fluids Eng. Division Summer Meeting. 2002.

[101 Frikha S., Coutier-Delgosha O., Astolfi J.A. Influence of the cavitation model on the simulation of cloud cavitation on 2D foil section // Intern. J. of Rotating Machinery. 2008. Vol. 2008. Article ID 146234. 12 p. doi:10.1155/2008/146234.

[Ill Черный С.Г., Чирков Д.В., Лапин В.Н. Численное моделирование течений в турбо-машинах. Новосибирск: Наука, 2006. 202 с.

[121 VAN LEER- В., Lee W.T., Roe P.L. Characteristic Time-Stepping or Local Preconditioning of the Euler Equations. AIAA Pap. 91-1552, 1991.

[131 R-ouse H., McNown J.S. Cavitation and pressure distribution, head forms at zero angle of yaw // Studies in Engineering. Bulletin 32. State Univ. of Iowa, 1948.

Поступила в редакцию 10 ноября 2010 г., с доработки — 4 февраля 2011 г.