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



УДК 532

С. Ю. Константинов, Д. В. Целищев

исследование и совершенствование численных МОДЕЛЕЙ КАВИТАЦИОННОГО МАССОпЕРЕНОСА

Рассмотрены основные аспекты численного моделирования процесса кавитации. Предложена новая математическая модель кавита-ционного массопереноса, учитывающая вязкость жидкости по уравнению Релея-Плессета, а также результаты ее верификации. Верификация проводилась на основании экспериментальных исследований процесса кавитации в струйном элементе типа «сопло-сопло». Кавитация; двухфазная среда; модели кавитации; модели кавитационного массопереноса

ВВЕДЕНИЕ

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

Моделирование кавитации является необходимым и наиболее сложным этапом при проектировании ряда гидравлических машин и устройств, таких как лопастные насосы, струйные элементы, гребные винты, форсунки, эжекторы и т.д. Для моделирования кавитации используют эмпирические, аналитические и численные модели. В последние 10 - 12 лет с развитием информационных технологий наибольшую популярность приобретает численное моделирование кавитации. Однако современные численные модели кавитации неадекватно визуализируют зону кавитации, имеют время расчета одной точки от 10 до 30 ч для достижения точности более 95%, требуют использования физического масштаба по времени для давлений свыше 20 атм. Все это приводит к необходимости исследования и совершенствования современных численных моделей кавитации.

1. СОСТОЯНИЕ ВОПРОСА

При моделировании кавитации выделяют два

Контактная информация: 8 (347) 273-09-44 Исследование выполнено при поддержке Министерства образования и науки Российской Федерации, соглашение 14.B37.21.0337 «Научное обоснование, создание и исследование энергосберегающих вихревых технологий фазоразделения, подогрева и редуцирования при транспортировке попутного и природного газа».

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

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

Динамика роста (замыкания) пузырька считается по упрощенному уравнению Релея:

ЖЯ = 12 рн - р

где Я - радиус кавитационного пузырька; рн - давление внутри пузырька (в модели -давление насыщенных паров); р - давление в жидкости (абсолютное давление решателя); р - плотность жидкости; Ар - действующий на пузырек перепад давлений. Подобное упрощение справедливо только для перепадов давления до 2 МПа [1]. На перепадах давления свыше 2 МПа модель перезаполняет расчетную ячейку паром и процесс расчета становится неустойчивым. Для устранения неустойчивости применяют расчет кавитации с малым числом Куранта в решателе: 25-50, однако данный способ приводит к значительному увеличению времени расчета (в 2-3 раза) и помогает лишь в 30-40% случаев. Аналогичным способом является задание физического масштаба по времени при расчете кавитации в ANSYS CFX [2]. В связи с этим разработчики пакетов вычислительной гидродинамики не гарантируют сходимость расчета для конкретных условий и геометрии [1].

2 , (1) 3 р,

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

4

а = П •— 7lRn

3 n '

(2)

где п - число пузырьков в объеме ячейки; Я0 - начальный радиус кавитационного зародыша, а - объемная доля пара в ячейке. Поскольку количество пузырьков в конкретный момент времени в объеме при моделировании неизвестно, то задачу решают путем введения в расчет некоторого начального радиуса пузырька Я0, что изменяет вид статистического компонента модели

П =

3а

4 Rn3 p

(3)

Singhal [1] предложил использовать статистический компонент в следующем виде:

F max(l, 4к )(а) , ^

г Irvap

а

где F - коэффициент парообразования (конденсации), о - коэффициент поверхностного натяжения, k - коэффициент релаксации по давлению, рг - плотность жидкости, Pvp - плотность пара. Из-за наличия коэффициента k модель обладала плохой сходимостью, что привело к серии работ по модифицированию статистического компонента.

Zwart, Gerber, Belamri [1] предложили следующий вид статистического компонента [1, 2]:

F

3<* nuc (1 -О. vap )Р

vap / г vap

vap

(5)

Rn

где ашс = 5^10-4 - коэффициент связи объемной доли с массовой [2].

Schnerr [1] модифицировал статистический компонент с учетом плотностей жидкости и пара, что позволило считать кавитацию в жидкостях сравнимых по плотности с своими парами:

^^а (1 -а )— Р Rn

(6)

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

2. МОДИФИЦИРОВАННАЯ МОДЕЛЬ КАВИТАЦИОННОГО МАССОПЕРЕНОСА

2.1. Методика модифицирования динамического компонента численной модели кавитационного массопереноса

Методика модифицирования динамического компонента численной модели кавитационного массопереноса состоит из четырех этапов:

1) задание исходных данных; 2) расчет безразмерной скорости роста в зависимости от числа Рейнольдса; 3) аппроксимация безразмерной скорости роста тремя аппроксимирующими функциями; 4) синтез нового динамического компонента численной модели кавитационно-го массопереноса. Методика выполнена в программном коде пакета Maple.

На первом этапе выполняется задание исходных данных, которыми являются:

1) интервал перепадов давлений, для которых выполняется коррекция численной модели кавитационного массопереноса. Интервал задается нижним и верхним значением;

2) шаг по перепаду давлений;

3) физические свойства жидкости, для которой корректируется динамический компонент численной модели кавитационного массопере-носа: вязкость ц и плотность р;

4) безразмерное время окончания роста кави-тационного пузырька от 3 до 10 по рекомендациям Ю. Л. Левковского [3].

Второй этап методики позволяет рассчитать зависимость безразмерной скорости роста кави-тационного пузырька от числа Рейнольдса. Этап состоит из трех частей:

1) решение безразмерного уравнения Релея-Плессета для роста пузырька с учетом вязкости жидкости:

b

d b 3

- + —

2

dx'

db dx

+

1 db Re • b dx

= 1

(7)

где Ь - безразмерный радиус кавитационного пузырька; Яе - число Рейнольдса для кавитационного пузырька; т - безразмерное время роста кавитационного пузырька. Число Рейнольдса, безразмерный радиус и безразмерное время роста пузырька выражаются следующим образом [3]:

2

Re =

Яол/ЛрЧ

i=R

4ц

R

R0'

(8)

т =

Ro\

Ap Р '

где Я0 = 10-6 м - начальный радиус пузырька, выбирается по рекомендациям из [1-3]; Я - текущий радиус пузырька;

2) во второй части результаты решения уравнения (1) комплектуются в единую переменную -двумерный массив;

3) в третьей части результаты расчета выводятся на экран в виде графика.

На этапе аппроксимации зависимости безразмерной скорости роста аппроксимируется тремя функциями. Функция аппроксимации №1:

£ = /(Яе) = tghK • Яе^ ), (9)

где w1 и w2 - коэффициенты аппроксимации. Функция аппроксимации №2:

£ = I(Яе) = ^ • Яе^ (10)

где w1 - коэффициент аппроксимации. Функция аппроксимации №3:

£=f (Re) =J2 (.gh(Wl

Rew2) + w3 ■ Rew4) ,(11)

где w1, w2, w3, w4 - коэффициент аппроксимации.

Аппроксимация выполняется средствами программного пакета Maple.

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

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

(12)

йЯ = йЪ Ар

Ж йх у р

Результатом работы методики является зависимость скорости роста кавитационного пу-

зырька от вязкости и плотности жидкости, которая рассчитывается от перепада давлений Ар и начального радиуса пузырька Я0.

2.2. Расчет нового динамического компонента численной модели кавитационного массопереноса

Синтез нового динамического компонента для численной модели кавитационного массопе-реноса выполняется в интервале давлений от 100 до 106 Па с шагом по давлению 1000 Па, для воды с плотностью 1000 кг/м3 и вязкостью 1004 сСт. Безразмерное время роста пузырька принимается равным х = 5 по рекомендациям [3]. Новый динамический компонент будет использоваться в численной модели кавитационного массопере-носа в пакете ANSYS CFX.

В результате выполнения расчета аппроксимирующие функция №1 принимает вид:

(13)

(14)

(15)

£ = I(Яе) ^3 «1,221 • Яе0 353 ) функция аппроксимации №2:

ЖЪ = I (Яе) = ^3 «1,428 • Яе)

функция аппроксимации №3:

ЖЪ = I (Яе) = ]3 (^(1,734 • Яе 0<457 ) -

- 0,078 • Яе - 0,61 0 ) .

Соответственно погрешность аппроксимации для функции №1 - 6,36%, функции № 2 -42,08%, функции №3 - 0,08%. Наиболее точная аппроксимация у функции №3, однако при внесении в код пакета моделирования ANSYS CFX данная функция является слишком сложной, поэтому выбирается функция №1, погрешность аппроксимации у которой составляет 6,36%.

Динамический компонент для модели кави-тационного массопереноса с учетом выбранной аппроксимирующей функции №1 принимает вид:

di = 3 ^221

(

RoVAp 4m

р

) ■

AP . (16) р

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

t

С учетом (16) и статистического компонента модели Zwart-Gerbert-Belamri (5) новая модель кавитационного массопереноса принимает вид:

ггв

О. Л ' 9И II

^ 3 ^ nuc (1 ^ vap ) Р vap

Р < РН, me = Fvap - X

R

0

(17)

x tgh(1.22 ((Рн - Р)P )035)

4m

2 Рн- p

3 p i

Р > Рн , mc = Fc

3 ^ va^P v

cond

R

x tgh(1.22 ((Р" " Р)P )035)

4m v

2 Р - Рн

3 Pi

3. ЭКСПЕРИМЕНТАЛЬНЫЕ ИССЛЕДОВАНИЯ СТРУЙНОГО

ЭЛЕМЕНТА ТИПА «СОПЛО-СОПЛО»

Для верификации численной модели кавитационного массопереноса проводились экспериментальные исследования струйного элемента типа «сопло-сопло» на стенде «Диагностика и идентификация гидросистем» УНИЦ «Гидропневмоавтоматика» и УГАТУ.

Характеристики кавитатора (рис. 1) снимались с давлениями в струйной камере 5, 10, 15 атм. При этом давление на выходе поддерживалось постоянным (давление в баке), а давление на входе переменным в диапазоне от 10 до 100 атм.

Рис. 1. Струйный элемент типа «сопло-сопло» (кавитатор): С = 1,6 мм; С = 1,25; Н = 0,93

После обработки результатов эксперимента были получены три расходно-перепадные характеристики струйного элемента при постоянном давление на выходе и переменном на входе для различных давлений в струйной камере, которые представлены на рис. 2.

/ f

/ -i Hht

-HI«* Lira

129

dp, ич

Рис. 2. Расходно-перепадная характеристика струйного элемента типа «сопло-сопло» при постоянном давлении на выходе и переменном давлении в струйной камере 5, 10, 15 атм

Для различных давлений в струйной камере расходно-перепадные характеристики практически совпадают (расхождение 1-3%) при одинаковом перепаде давлений. Это можно объяснить наличием гидравлического разрыва потока в струйной камере в сочетание с эффектом эжек-ции. Расход на выходе из струйного элемента представляет собой сумму расхода на входе и расхода, поступающего в струйную камеру. Поскольку статическое давление в струе зависит от перепада давлений между входом и выходом элемента, то расход эжекции зависит только от давления в струе и не зависит от давления в струйной камере. Данный феномен объясняется постоянством гидравлического сечения между струей из струйной трубки и приемным соплом. Данное сечение запирает расход эжекции и тем самым обеспечивает постоянный расход на сливе струйного элемента.

4. численное моделирование

КАВИТАЦИИ В СТРУЙНОМ ЭЛЕМЕНТЕ ТИПА «СОПЛО-СОПЛО»

Численное моделирование проводилось в пакете прикладных программ ANSYS CFX для экспериментальной характеристики при давлении в струйной камере 10 атм. Давление на выходе при этом - 1,9 атм. Расход на входе менялся от 7 до 12 л/мин с интервалом 0,5 л/мин. Температура рабочей жидкости - 32 °С. Рабочая жидкость: масло Shell T46 с плотностью 872 кг/м3 и вязкостью 80 сСт. Расчетная сетка имеет число узлов 985 013. Значение невязок (погрешности) при котором расчет считался сошедшимся 10-4.

X

Моделирование выполнялось тремя сериями по 11 расчетов: 1) без учета кавитации; 2) с учетом кавитации моделью Zwart-Gerber-Belamri (ZGB):

3« пис (1 )Р шр

Р < Рн > me = Fva

Р > р , mc = Fc0

2 Рд - Р

3 Р, '

3ДР2 Р - Рш

(18)

4 3 Р;

3) новой моделью кавитационного массопере-носа (17). В качестве модели кавитации используется гомогенная модель кавитации «Эйлера» с моделью турбулентности к-е.

Расходно-перепадные характеристики, полученные в результате моделирования, представлены на рис. 3.

—*—Б« учета шняднн Модна гОВ

Рис. 3. Результаты численного моделирования в сопоставлении с экспериментом

5. ВЕРИФИКАЦИЯ ЧИСЛЕННЫХ МОДЕЛЕЙ КАВИТАЦИИ НА ОСНОВАНИИ НАТУРНОГО ЭКСПЕРИМЕНТА

Для верификации новой численной модели на основании результатов натурного эксперимента были выбраны три критерия: 1) погрешность моделирования; 2) время моделирования; 3) качественная оценка моделирования.

5.1. Оценка погрешности моделирования

Для расчета относительной погрешности по расходу используется зависимость:

6 - 6,

G =

100% ,

(19)

где о - относительная погрешность в %, 6 -расход при моделировании; 6э - расход при эксперименте. На рис. 4 представлена зависимость относительной погрешности от перепада давления на струйном элементе при моделировании для каждой серии расчетов.

Л/», атм

Рис. 4. Изменение погрешности в зависимости от перепада давлений

Средняя относительная погрешность считалась среднеарифметическим от погрешностей для выбранной серии. Максимальная средняя погрешность фиксируется у численного моделирования без учета кавитации: 22,56%; наименьшая погрешность у модели Zwart-Gerber-Belamri -2,11%. Погрешность новой численной модели кавитационного массопереноса составляет 3,99%, что вполне приемлемо для численного моделирования. Однако по сравнению с моделью Zwart-Gerber-Belamri новая модель учитывает влияние сил вязкости, которые оказывают значительное влияние на динамику пузырьков, особенно при нестационарном течении.

5.2. Оценка времени моделирования

Для оценки продолжительности моделирования оценивается число итерации, необходимое для полной сходимости решения в зависимости от перепада давления с учетом того, что все расчёты производятся на одной расчетной сетке. Зависимость числа итераций для новой численной модели кавитации и численной модели Zwart-Gerber-Belamri от перепада давлений представлена на рис. 5.

.■V

-*-Можль£«1 —»-Нммрмнгпь ' ЭТМ

Рис. 5. Зависимость числа итераций, необходимых

для получения решения, для новой численной модели кавитации и модели Zwart-Gerber-Belamri

Полное число итераций для получения серии расчетов с моделью Zwart-Gerber-Belamri равно 3617, новой численной модели кавитации - 1562, что дает экономию во времени 56,8%.

Данный результат делает новую численную модель кавитации более конкурентоспособной. Следует отметить, что с ростом перепада давлений растет и время расчета, что отмечается и другими исследователями [4, 5], однако на больших перепадах давления на сходимость существенно влияет качество сетки [5], с улучшением сетки, как правило, уменьшается число итераций, однако общее время расчёта может расти [6].

5.3. Качественная оценка результатов моделирования

Качественная оценка численного моделирования проводится с целью определения адекватности полученной кавитационной картины. Для качественной оценки моделирования сравнивается визуализация объемной доли пара в расчетной области для обеих численных моделей кавита-ционного массопереноса. Визуализация объемной доли пара для перепада давлений 39,5 атм. представлена на рис. 6.

Рис. 6. Визуализация объемной доли пара для перепада давления 39,5 атм: а - модель Zwart-Gerber-Ве1атп; б - новая численная модель кавитации (17)

Сравнение результатов изображенных на рис. 6, а - б показывает, что зона кавитации отображается моделью Zwart-Gerber-Be1amri

неверно: в осесимметричном канале не может возникать несимметричный отрыв потока показанный для наглядности стрелкой. Также следует считать ошибочным искажение зоны кавитации после диффузора в проточной части кавитатора, что отмечено стрелкой. В то же время у новой модели кавитационного массоперно-са отсутствует подобное искажение. Причина искажений кроется в максимальной объемной доле пара. Так, в случае использования модели Zwart-Gerber-Be1amri она составляет 0,75, а в новой модели - 0,67. Данный факт показывает, что новая модель уменьшает объемную долю пара, тем самым увеличивая расход жидкости через проточную часть и выправляя поле давлений и скоростей, что ведет к получению более реалистичной картины кавитационной зоны.

ВЫВОДЫ

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

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

3. Верификация результатов численного моделирования кавитации и эксперимента показала, что погрешность при расчете без кавитации составляет 22%, с учетом кавитации посредством модели «Zwart-Gerber-Be1amri» не более 2,5%, а с помощью новой модели не более 4%, что находится в пределах допускаемой погрешности. При этом новая модель экономит до 57% машинного времени по сравнению с моделью «Zwart-Gerber-Be1amri», точнее визуализирует течение и учитывает дополнительный параметр, такой как вязкость жидкости.

СПИСОК ЛИТЕРАТУРЫ

1. ANSYS FLUENT 12.0 Theory Guide. April 2009. ANSYS Inc.

2. ANSYS CFX 12.0 Theory Guide. April 2009. ANSYS Inc.

3. Левковский Ю. Л. Структура кавитационных течений. Л.; Судостроение, 1978. 224 с.

4. Vortmann C. Untersuchen zur Thermodynamik des Phasenubergangs bei der numerischen Berechnung kavitierender Dusenstromungen. Karlsruhe, 2001. 132 s.

5. Пермяков Г. С., Целищев Д. В. Исследование эффекта стабилизации расхода в струйных элементах // Вестник УГАТУ, 2010, Т. 14, № 2 (37). С. 21-29.

6. Ахметов Ю. М., Калимуллин Р. Р. и др. Ис-

следование гидродинамических и термодинамических процессов высоконапорного многофазного вихревого течения жидкости // Вестник УГАТУ, 2012. Т. 16, № 2 (47). С. 163-168.

ОБ АВТОРАХ

Константинов Сергей Юрьевич, аспирант каф. прикл. гидро-мех. Дипл. магистра (УГАТУ 2012) Иссл. в обл. математического моделирования кавитационных течений.

Целищев Дмитрий Владимирович, доцент той же каф. Канд. техн. наук по гидравл. машинам (УГАТУ, 2010). Иссл. в обл. электрогидравл. рулевых приводов для систем упр-я летательн. аппаратами.