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

Размещено на http://www.Allbest.Ru/

Размещено на http://www.Allbest.Ru/

Размещено на http://www.Allbest.Ru/

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

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

Л.В. Панов, Д.В. Чирков

г. Новосибирск

Введение

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

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

Режим работы гидротурбины по отношению к явлению кавитации определяется кавитационным коэффициентом установки (числом Тома) у. Направим ось z расчетной системы координат вниз (рис. 1). Тогда, согласно стандарту МЭК 60193 [2],

, (1)

где

, (2)

(3)

В (2)-(3) индекс «1» соответствует входному сечению S₁ спиральной камеры, индекс «2» - выходному сечению S₂ отсасывающей трубы (см. рис. 1). Давление p_abs - среднее абсолютное давление в соответствующем сечении, p_V - абсолютное давление насыщенного водяного пара при данной температуре воды. z_r - опорный уровень, за который для радиально-осевых гидротурбин принимается уровень расположения средней линии направляющего аппарата. Величины скорости v в сечениях 1, 2 определяются, согласно стандарту, как средние расходные: , i = 1,2.

Для оценки воздействия кавитации на работу турбомашины проводят серию испытаний модельной турбины на лабораторном стенде, в ходе которых напор H остается постоянным, а коэффициент у варьируется путем создания разрежения в вакуумном баке. В результате испытаний для заданного режима работы строятся зависимости расхода, мощности и КПД от у (рис. 2). При больших и средних значениях у кавитации нет, либо она незначительна, КПД здесь остается практически постоянным. Однако, при уменьшении у ниже некоторого критического (срывного) значения у_s, КПД и расход, пропускаемый турбиной, начинают резко падать. Вид зависимости очень важен при проектировании гидротурбины. Необходимо, чтобы критический коэффициент у_s был с запасом меньше кавитационного коэффициента у_pl, соответствующего условиям работы натурной турбомашины на станции.

Рис. 1. Схема модели гидротурбины



Рис. 2. Пример экспериментальных зависимостей КПД и расхода от у

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

За последние пятнадцать лет в вычислительной гидродинамике предложен ряд моделей для описания двухфазных течений «жидкость-пар» с фазовыми превращениями. Наибольшее распространение получили простые односкоростные модели, в которых уравнения Навье-Стокса для смеси «жидкость-пар» замыкаются уравнением переноса массовой (или объёмной) доли жидкости или пара [5-8]. Эти модели позволяют учесть неравновесность процессов конденсации и парообразования. В [7, 9] подобные модели применены для расчета кавитационных течений в насосах, в [10] - для радиально-осевых турбин различной быстроходности. В [7,9,10] показано достаточно хорошее соответствие расчетов с данными эксперимента.

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

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

1. Основные уравнения

Для моделирования кавитационного течения жидкости в проточном тракте гидротурбины используется квазигомогенная изотермическая модель движения сжимаемой смеси «жидкость-пар», в которой предполагается, что скорости жидкой и паровой фаз совпадают, а распределение объемной доли жидкой фазы_L может быть описано уравнением переноса с источниковыми членами, отвечающими за парообразование и конденсацию:

,(4)

,(5)

(6)

Здесь с - плотность смеси, [кг/м³];

v - скорость, [м/с];

t - время, [с]; ,

p - давление, [Па];

k - кинетическая энергия турбулентных пульсаций, [м²/с²].

При расчёте течения во вращающейся системе координат вектор массовых сил f, помимо силы тяжести, включает центробежную и кориолисову силы:

,

где щ - угловая скорость вращения РК. В (4)-(6) плотность смеси «жидкость-пар» рассчитывается по формуле



где с_L - плотность жидкости,

с_V - плотность пара.

В (6) компоненты тензора вязких напряжений имеют вид:

,

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

,

где м_L, м_V - динамические коэффициенты вязкости жидкости и пара. Для определения турбулентной вязкости смеси система (4)-(6) дополняется стандартной k-е моделью турбулентности.

Динамика объёмной доли жидкой фазы описывается уравнением переноса (6) с источниковыми членами, отвечающими за конденсацию пара (m⁺) и испарение (m-). Вид этих членов для четырёх моделей, рассмотренных в настоящей работе, представлен в таблице 1.

Таблица 1

Источниковые члены для уравнения (6)

Модель	m+	m?	Константы
Модель 1 [5]			Сprod= 80 Cdest = 1
Модель 2 [6]			Сprod= 200 Cdest = 105
Модель 3 [7]			Сprod= 0.137 Cdest= 0.274
Модель 4 [8]			-

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

Для численного интегрирования уравнений (4)-(6) применяется метод, изложенный в предыдущей работе авторов [11] и являющийся распространением метода [4]. В стационарном случае уравнения (4)-(6) переписываются в виде

(7)

В системе (7) в соответствии с подходом искусственной сжимаемости введено псевдовремя, по которому осуществляется установление, - параметр искусственной сжимаемости. В соответствии с [12], в последнее уравнение системы внесен дополнительный член , делающий уравнение переноса _L совместным с уравнением неразрывности. - вектора невязких и вязких потоков в направлении x_i. Уравнения системы (7) решаются совместно, поскольку предварительное тестирование показало, что совместное решение существенно ускоряет сходимость. Система (7) аппроксимируется неявным методом конечных объёмов с противопотоковой MUSCL-схемой для невязких членов. Для нелинейных источниковых членов используется специальная явно-неявная аппроксимация [11]. Полученная в результате дискретизации система линейных алгебраических уравнений на приращение неизвестных Q при переходе со слоя s на слой s+1 по псевдовремени решается приближенно, путем замены матрицы А на произведение матриц A₁A₂, каждая из которых обращается бегущим счетом [11].

В отличие от [11], в настоящей работе введено ограничение снизу на объёмную долю жидкости в ячейке . Как показали численные эксперименты, такое ограничение не привело к изменению расчётных значений интегральных величин момента на рабочем колесе и расхода. Кроме того, при расчёте невязких потоков через грани ячейки по схеме MUSCL

внесён специальный ограничитель

как в работе [12], за счёт которого порядок аппроксимации снижается до первого в областях больших градиентов б_L. Перечисленные модификации позволили существенно улучшить сходимость метода и решить проблему сходимости при малых плотностях пара, которая имела место в [11].

Шаг по псевдовремени вычислялся локально в каждой ячейке по числу Куранта CFL=1. Такое число Куранта обеспечивает скорость сходимости кавитационных расчётов, равную скорости сходимости расчета по модели несжимаемой жидкости, а значит, получение решения за приемлемое время.

Расчет течения в гидротурбине осуществлялся в области, включающей направляющий аппарат, рабочее колесо и отсасывающую трубу. При этом использовалась циклическая постановка, т.е. течение моделировалось только в одном межлопаточном канале НА и одном межлопастном канале РК, в предположении, что течение в остальных каналах циклически повторяется. При передаче данных между областями НА и РК, РК и ОТ осуществлялось их усреднение в окружном направлении [4].

пар кавитация турбулентность энергетический гидротурбина

3. Задание граничных условий

Введем обозначение для полной энергии потока в выходном сечении ОТ:

, где (8)

Тогда из (1)-(3), с учетом того, что, получим:

(9)

Таким образом, при заданном значении у и известном напоре H остается фиксированной интегральная величина полной энергии потока E_2,IEC в выходном сечении ОТ. При этом расход жидкости и давление в точке 2 по-отдельности априори не известны. Таким образом, при численном решении уравнений (4)-(6) в выходном сечении ОТ держалась полная энергия E_2,IEC, заданная по формуле (9). Давление p_abs вычислялось как среднее по сечению.

На входе в направляющий аппарат (НА) держался угол входа потока и величина полной энергии, равная

, (10)

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

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

Разработанный метод применен для расчета кавитационных течений в проточном тракте моделей (D₁=0.46 м) двух гидротурбин с параметрами быстроходности , равными 313 и 240. Параметры турбин и выбранных режимов указаны в таблице 2. Во всех расчётах учитывалась сила тяжести g =9.81 м/с².

Таблица 2

Параметры моделей гидротурбин

	Режим	Открытие НА a0, мм	Частота вращения РК n, об/мин	Напор, H, м
Турбина 1, ns = 313	Режим оптим. КПД (1)	28	797.6	20.5
	Режим макс. мощности (2)	36	978.9	20.5
Турбина 2, ns = 240	Режим макс. мощности	34	793	24.5

Полученное в расчёте значение расхода Q_c корректировалось по формуле . КПД вычислялся через момент M на валу РК и расход

, (10)

где з_vol - объёмный КПД,

з_mech - механический КПД, которые оценивались по инженерным формулам.

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

Методические расчеты турбины 1

На примере режима максимальной мощности для турбины 1 исследованы сеточная сходимость и влияние плотности пара с_V на кавитационные характеристики.

Базовая сетка состояла из 62-25-25 ячеек НА, 90-25-25 ячеек РК и 57930 ячеек ОТ. В подробной сетке число узлов в РК и ОТ увеличено в 1.5 раза по каждому пространственному направлению. На рис. 3 показано влияние сетки на кавитационные кривые , и при использовании модели кавитации 1 [5]. Видно, что на базовой сетке для бескавитационного режима (= 0.38) КПД на 1% ниже, чем в эксперименте. На подробной сетке абсолютные значения всех интегральных параметров ближе к эксперименту. На подробной сетке при = 0.25 стационарное решение получить не удалось, диапазон колебаний величин M, Q и казан на рис. 3 символами «».

В расчетах на подробной сетке на кривых и наблюдается провал в диапазоне = 0.25-0.27, который отсутствует в эксперименте. Замечено, что получаемая в расчетах на этих режимах паровая каверна на тыльной стороне лопасти схлопывается вблизи ее выходной кромки. В то же время, в области схлопывания паровой каверны, как и в стационарных расчетах кавитационного обтекания затупленных тел [11,12], происходит значительное искривление линий тока с образованием зоны возвратного течения. Эта особенность потока приводит к отклонению угла схода потока с выходной кромки лопасти, и, как следствие, к изменению величин М и Q. С другой стороны, в эксперименте паровая каверна «дышит», она зарождается и отрывается от поверхности лопасти выше по потоку (см. рис. 9), не возмущая угол схода потока с лопасти. Таким образом, мы считаем, что вероятной причиной провала при = 0.25 - 0.27 в стационарных расчетах является искажение потока вблизи выходной кромки лопасти. Эта гипотеза требует проверки с использованием нестационарной модели течения.

Вне диапазона = 0.25-0.27 качественное поведение , , и расчетные значения _s на базовой и подробной сетках совпадают. Поэтому для дальнейших расчетов использована базовая сетка.

Рис. 3. Влияние сетки на момент, расход и КПД турбины. с_V=1 кг/м³

Для воды при T = 17°C плотность пара с_V = 0.01 кг/м³. Расчеты с таким малым с_V сопряжены с вычислительными трудностями. Однако, как показано в [11,12] в практических расчетах величина с_V может быть увеличена. На рис. 4. показано влияние с_V на интегральные характеристики работы турбины 1. Значения момента, расхода и КПД приведены по отношению к их значениям для бескавитационного режима. Видно, что при с_V 10 кг/м³ результаты практически не меняются. Поэтому дальнейшие кавитационные расчеты проведены с с_V = 1 кг/м³.

Известно, что выбор модели турбулетности сильно влияет на картину кавитационного течения, поэтому нами проведено сравнение трёх известных двухпараметрические моделей - стандартной k-е, RNG k-е и k-е Кима-Чена.

Модели RNG и Кима-Чена существенно хуже сходились - удалось получить результаты только для режима оптимального КПД. На этом режиме (рис. 6) в диапазоне = 0.15-0.20 у моделей RNG и Кима-Чена наблюдается существенное увеличение момента и КПД, которое не наблюдается в эксперименте и которого нет на графике для стандартной k-е модели.

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

Рис. 4. Влияние с_V на момент, расход и КПД турбины



Рис. 5. Сравнение различных моделей турбулентности для АКО постановки. Богучанская. Режим оптимального КПД

Моделирование течения в гидротурбине 1 (n_s = 313)

На рис. 6 изображены зависимости момента, расхода и КПД установки от у для оптимального режима работы турбины 1, полученные по моделям 1-4 (см. табл. 1). Видно, что качественное поведение всех трех величин, а также наклон графиков при малых у близки к полученным в эксперименте. Однако в расчетах кавитационный срыв наступает раньше, при на 0.05 больше, чем в эксперименте.

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

Важно отметить, что все четыре рассмотренные модели парообразования и конденсации дают очень близкие количественные результаты, несмотря на разные предположения, принятые при их выводе. Аналогичный вывод сделан в [11] на задаче обтекания затупленного цилиндра. Этот факт позволяет утверждать, что проблема обнаруженного на рис. 6 и 7 сдвига кавитационных кривых не связана с выбором констант моделей.

Рис. 6. Влияние у на момент, расход и КПД турбины. с_V = 1 кг/м³. Сравнение моделей кавитации 1-4. Режим оптимального КПД

Рис. 7. Влияние у на момент, расход и КПД турбины. с_V = 1 кг/м³. Сравнение моделей кавитации 1 - 4. Режим максимальной мощности

На рис. 8 для режима оптимального КПД представлено распределение пара и модуль вектора скорости на тыльной стороне лопасти для нескольких значений коэффициента у: от у = 0.20, когда кавитации практически нет, до у = 0.09, когда интенсивная кавитация дает в расчете падение КПД более чем на 5%. На рис. 9 приведено сравнение фотоснимка эксперимента при у = 0.175 и результата численного моделирования. В эксперименте кавитация на этом режиме проявляется в виде отдельных больших пустот, периодически формирующихся на тыльной поверхности лопасти и схлопывающихся ниже по потоку. Таким образом, одной из возможных причин завышения значения у_s в расчете является использование стационарной постановки.

у = 0.20 у = 0.17 у = 0.15 у = 0.09

Рис. 8. Распределение объемной доли пара и модуля вектора скорости на тыльной стороне лопасти. Пунктиром показана граница, на которой p = p_V.

Рис. 9. Фотоснимок течения в эксперименте (слева) и изоповерхность у_L=0.3 в расчете (справа). у = 0.175. Пунктиром показана область, где пара 1%

Моделирование течения в гидротурбине 2 (n_s = 240)

На рис. 10 представлено сравнение момента, расхода и КПД с графиками кавитационных испытаний для турбины 2. Здесь также все четыре рассмотренные модели дают близкие, качественно соответствующие эксперименту результаты. В эксперименте наблюдается увеличение КПД на 0.5% перед срывом, в расчётах этого увеличения нет. Для данной турбины отклонение рассчитанного значения у₁ от экспериментального у₁ = 0.073 существенно меньше, около 0.02.

Рис. 10. Влияние у на момент, расход и КПД турбины PO910. с_V = 1 кг/м³. Сравнение моделей кавитации 1-4

Заключение

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

Показано, что для получения интегральных параметров потока, не зависящих от значения плотности пара, достаточно положить с_V = 1 кг/м³. Проведено сравнение четырех моделей парообразования и конденсации. Показано, что все модели, включая модель [8], не имеющую эмпирических констант, дают сопоставимые результаты. Для турбин различной быстроходности проведено сравнение кривых , и с экспериментальными данными. Показано хорошее качественное соответствие. Однако значение критического кавитационного коэффициента в расчетах завышено, особенно для быстроходной турбины. Возможной причиной подобного расхождения является использование стационарной постановки. Кроме этого, в дальнейшем представляет интерес учесть в модели кавитационного течения нерастворимый газ, который в небольшой концентрации присутствует в потоке.

Предложенная методика позволяет за приемлемое время (характерное время счёта - 2 часа на ПК с частотой процессора 3 ГГц) проводить кавитационные расчёты потока и сравнительную оценку кавитационных характеристик при проектировании проточных частей гидротурбин.

Литература

1. Avellan F. Introduction to cavitation in hydraulic machinery // The 6th Intern. Conf. on Hydraulic Machinery and Hydrodynamics, Timisoara, Romania, 21-22 October - 2004 - P. 11-22.

2. Hydraulic turbines, storage pumps and pump-turbines - Model acceptance tests. IEC Standard 60193, International Electrotechnical Commission, 1999.

3. Топаж Г.И. Расчет интегральных гидравлических показателей гидромашин. - Л.: Изд-во Ленинградского университета, 1989. - 208 с.

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

5. Singhal A.K., Vaidya N., Leonard A.D. Multi-dimensional simulation of cavitating flows using a PDF model for phase change // ASME Fluids Engineering Division Summer Meeting, ASME Paper FEDSM97-3272. - 1997.

6. Lindau J.W., Kunz R.F., Venkateswaran S., Boger D.A. Application of preconditioned, multiple-species, Navier-Stokes models to cavitating flows // 4th International Symposium on Cavitation, Pasadena, California, 20-23 June. - 2001 - 14 p.

7. Athavale M.M., Singhal A.K. at el. Application of the Full Cavitation Model to Pumps and Inducers // Intern. J. of Rotating Machinery. - 2002. - Vol. 8. -№ 1. - P. 45-56.

8. Senocak I., Shyy W. Evaluation of cavitation models for Navier-Stokes computations// FEDSM 02. ASME. Fluid Engineering Division Summer Meeting, Montreal, Canada. - 2002.

9. Bouziad Y.A. et al. Experimental and Numerical Cavitation Flow Analysis of an Industrial Inducer // 22nd IAHR Symposium on Hydraulic Machinery and Systems, Stockholm, Sweden, June 29 - July 2. - 2004 - 9 p.

10. Kurosawa S., Lim S. M. and Enomoto Y. Virtual model test for a Francis turbine / 25th IAHR Symposium on Hydraulic Machinery and Systems, 2010 10 p.

11. Панов Л.В., Чирков Д.В., Чёрный С.Г. Численные алгоритмы моделирования кавитационных течений вязкой жидкости // Вычислительные технологии. - 2011. - Т. 16, №4. - С. 96-113.

12. Kunz R.F., Boger D.A., Stinebring D.A. et al. A preconditioned Navier-Stokes method for two-phase flows with application to cavitation prediction // Computers & Fluids - 2000 - №29 - P. 849-875.

Размещено на allbest.ru