Теплоэнергетика, 2021, № 8, стр. 17-25

ПОСТАНОВКА ЗАДАЧИ

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

Математическое описание включает в себя систему уравнений, состоящую из уравнений Навье–Стокса для сжимаемого газа, усредненных по Рейнольдсу и Фавру [1], уравнения состояния, а также уравнения турбулентности для замыкания системы. Ударные волны не моделируются, поэтому объемная вязкость в уравнениях движения не учитывается [2].

При получении уравнений для сжимаемого газа плотность, давление и тепловой поток осредняются по Рейнольдсу следующим образом:

где $\bar {\rho },$ $\bar {p}$ и $\overline {{{q}_{i}}} $ – средние значения плотности, давления и удельного теплового потока; $\rho {\kern 1pt} ',$ $p{\kern 1pt} '$ и $q_{i}^{'}$ – их пульсации; индекс i обозначает проекцию вектора на ось координат.

Скорость, энтальпия, внутренняя энергия и температура осредняются по Фавру:

где $\widetilde {{{u}_{i}}},$ $\tilde {h},$ $\tilde {e}$ и $\tilde {T}$ – средние значения скорости потока, энтальпии, внутренней энергии и температуры; и – их пульсации.

В конечном итоге получается следующая система уравнений для сжимаемого газа, включающая в себя уравнения неразрывности, движения и энергии:

Здесь $\overline {{{{{\tau }}}_{{ji}}}} $ – тензор вязких напряжений; $\overline {{{{{\tau }}}_{{t,ji}}}} $ = =  – тензор турбулентных напряжений; – турбулентный тепловой поток; k – кинетическая энергия турбулентности. Слагаемое в уравнении энергии учитывает молекулярный и турбулентный перенос кинетической энергии турбулентности. При записи системы уравнений (1) компоненты тензоров записываются с помощью индексов i и j, при этом используется метод немого суммирования по дважды повторяющемуся индексу. Индекс t указывает на турбулентные составляющие.

Следующим важным шагом является описание взаимосвязи между термодинамическими характеристиками рабочего тела – плотностью, давлением и температурой, для чего используется уравнение состояния реальных газов Редлиха –Квонга [3], которое имеет следующий вид:

где R – универсальная газовая постоянная; V_m – молярный объем газа; a и b – константы, зависящие от конкретного вещества и определяемые по критическим значениям температуры ${{T}_{{cr}}}$ и давления ${{p}_{{cr}}}{\text{:}}$

Если ввести коэффициент сжимаемости $Z(p,T) = {{p{{V}_{m}}} \mathord{\left/ {\vphantom {{p{{V}_{m}}} {\left( {RT} \right)}}} \right. \kern-0em} {\left( {RT} \right)}},$ то уравнение состояния (2) преобразуется в кубическое уравнение относительно Z:

где коэффициенты A и B рассчитываются как

Учитывая, что молярный объем можно представить как ${{V}_{m}} = {M \mathord{\left/ {\vphantom {M {{{\rho }}\left( {p,T} \right)}}} \right. \kern-0em} {{{\rho }}\left( {p,T} \right)}}$ (М – молярная масса), итоговую зависимость плотности от давления и температуры можно записать в виде

где Z(p,T) – единственный в области применимости уравнения (2) положительный действительный корень уравнения (4).

Уравнение Редлиха – Квонга (2) может использоваться в интервале давлений и температур, для которых справедливо неравенство ${p \mathord{\left/ {\vphantom {p {{{p}_{{cr}}}}}} \right. \kern-0em} {{{p}_{{cr}}}}} < 0.5{T \mathord{\left/ {\vphantom {T {{{T}_{{cr}}}}}} \right. \kern-0em} {{{T}_{{cr}}}}}.$ В данной работе расчеты проводятся для давлений и температур, находящихся в области применимости этого уравнения.

В силу того что связь между турбулентными составляющими тензора напряжений ${{\tau }_{t}}$ и вектора плотности теплового потока ${{\vec {q}}_{t}}$ с параметрами осредненного течения неизвестна и должна быть определена, система (1) оказывается незамкнутой. Чтобы ее замкнуть, необходимо использовать одну из моделей турбулентности. В настоящей работе применяется обновленная k–ω-модель турбулентности, так как она показала наилучшую устойчивость при расчетах сжимаемых течений [4]. Уравнения указанной модели турбулентности в общем виде записываются следующим образом:

где ω – скорость диссипации энергии; G_k, G_ω – члены, отвечающие за генерацию k и ω; Γ_k, Γ_ω – коэффициенты диффузии; Y_ω, Y_k – коэффициенты затухания. В настоящей работе из-за громоздкости изложения выкладки для членов G_k, G_ω, Γ_k, Γ_ω и Y_ω, Y_k не приводятся.

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

где h_хол – холодопроизводительность; h_s – располагаемый теплоперепад; ${{h}_{{\text{н}}}}$ и ${{h}_{{\text{к}}}}$ – значения энтальпии в начальной и конечной точках процесса расширения в h, s-диаграмме.

Эффективность ступени турбодетандерного агрегата наглядно можно показать на h, s-диаграмме процесса детандирования, для построения которой, а также для определения изоэнтропийного КПД необходимо предварительно вывести зависимости энтропии и энтальпии от давления и температуры. Поскольку в настоящей работе в качестве рабочего тела рассматривается газообразный гелий (ТУ 0271-135-31323949-2005), все зависимости будут приведены для этого вещества. Предварительные расчеты показали, что при использовании соотношений для идеального газа возникают серьезные отличия расчетных значений от экспериментальных (по температуре на выходе разница составила около 3 К), поэтому необходимо применить зависимости для реального газа. Также следует отметить, что все расчетные величины (КПД, энтальпии) определяются по статическим параметрам потока на выходе из ступени проточной части турбодетандера.

Энтропия в общем виде может быть представлена в виде зависимости

где s₀(T) – энтропия вещества в идеально-газовом состоянии [5]; p – давление, МПа.

Для одноатомного газа теплоемкость с_p0 в идеально-газовом состоянии не зависит от температуры и определяется по формуле

Это существенно упрощает расчет энтропии в идеально-газовом состоянии:

где ${{s}_{{00}}}\left( {{{T}_{{01}}}} \right)$ – значение энтропии при температуре начала отсчета Т₀₁.

В настоящей работе при температуре Т₀₁ = 20 К принимается значение энтропии равное ${{s}_{{00}}}$ = = 17.462 кДж/(моль ⋅ К).

Формулу для вычисления энтальпии для газообразного гелия можно представить в виде [5]

где h₀ – энтальпия в идеально-газовом состоянии; А₃ – эмпирический коэффициент.

Определение энтальпии в идеально-газовом состоянии также упрощается потому, что теплоемкость ${{с}_{{p0}}}$ не зависит от температуры, т.е.

где ${{h}_{{00}}}$ – энтальпия при температуре ${{T}_{{02}}}$ = 100 К, а $h_{0}^{0}$ – теплота сублимации гелия при Т = 0 К.

Значение теплоты сублимации гелия принято на основании данных [5] и равно $h_{0}^{0}$ = = 59.00 кДж/кмоль = 14.7404 кДж/кг. Значение энтальпии при температуре начала отсчета ${{T}_{{02}}}$ = = 100 К равно ${{h}_{{00}}}$ = 2.0785 кДж/кг.

Коэффициент A₃ получен усреднением серий равноточных уравнений состояния для соответствующих интервалов температур [5]:

где b – матрица коэффициентов, полученных методом наименьших квадратов; m и n – размерность матрицы b; $w = {{{{\rho }}(p,T)} \mathord{\left/ {\vphantom {{{{\rho }}(p,T)} {{{{{\rho }}}_{{cr}}}}}} \right. \kern-0em} {{{{{\rho }}}_{{cr}}}}}$ – приведенная плотность; $\tilde {T} = {T \mathord{\left/ {\vphantom {T {{{T}_{{cr}}}}}} \right. \kern-0em} {{{T}_{{cr}}}}}$ – приведенная температура.

При расчетах приняты следующие значения критических параметров: Т_cr = 5.19 К, p_cr = 222 740 Па, ρ_cr = 69.64 кг/м³.

На рис. 1 показана h, s-диаграмма, построенная на основе представленных выше формул, в диапазоне давлений от 0.02 до 1.5 МПа и температуры от 20 до 300 К.

Изотермы на диаграмме представляют собой горизонтальные линии, что соответствует идеальному газу, энтальпия которого зависит только от температуры. Однако в действительности для давлений выше 0.12 МПа наблюдается отклонение изотерм от горизонтальных линий, т.е. появляется зависимость энтальпии от давления, что дополнительно иллюстрируется увеличенным фрагментом на рис. 1.

Трехмерная модель ступени турбодетандерного агрегата приведена на рис. 2. Принцип ее работы довольно хорошо известен: газ через входные отверстия 1 поступает на направляющий аппарат 2, после чего, расширяясь, совершает работу, вращая рабочее колесо с определенной частотой, а затем через диффузор 4 попадает в выходное отверстие 5 и далее в магистраль (не показана).

На основе трехмерной модели строится модель проточной части, которая в дальнейшем разбивается на элементарные объемы, из которых формируется объемная расчетная сетка. В отличие от представленных в литературе расчетных исследований [6], в которых каждая часть турбинной ступени рассчитывается отдельно, в данной работе используется разбиение проточной части на зоны с общими скользящими интерфейсами, что позволяет провести расчеты полной модели. Таким образом, итоговая расчетная модель состоит из четырех зон: направляющий аппарат (НА), рабочее колесо (РК), зазор между диффузором и рабочим колесом и диффузор (Д).

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

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

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

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

Места контакта расчетных зон турбинной ступени НА – РК и РК – Д заменялись скользящими интерфейсами, что позволило провести комплексный расчет агрегата в сборе, а не по частям.

На входных отверстиях 1 (см. рис. 3) задается условие pressure-inlet, т.е. задаются давление ${{p}_{{\text{н}}}}$ и температура на входе ${{T}_{{\text{н}}}}.$ На выходное отверстие 2 накладывается условие pressure-outlet, т.е. задается давление в выходной магистрали ${{p}_{{\text{к}}}},$ при этом температура на выходе из диффузора ${{T}_{{\text{к}}}}$ является расчетным параметром. Для каждой расчетной точки задается частота вращения рабочего колеса Ω в соответствии с протоколом испытаний. Рассматривается однофазное и однокомпонентное рабочее тело без возможных примесей, поэтому фазовые переходы отсутствуют.

ЭКСПЕРИМЕНТАЛЬНЫЕ ДАННЫЕ

Экспериментальные данные предоставлены авторам сторонней организацией. На рис. 4 в упрощенном виде показана принципиальная схема экспериментальной установки. Были проведены статические обкаточные испытания на холостом ходу, т.е. без имитации нагрузки.

Из баллона 1 через клапан 2 рабочее тело (гелий) поступает во входную группу турбодетандерного агрегата (НА), затем, вращая рабочее колесо, расширяется и через диффузор попадает в магистраль через отсечной кран 10. Частота вращения рабочего колеса определяется датчиком 7. Параметры рабочего тела на входе регистрируются датчиками давления и температуры 3, а на выходе – датчиками 8. Регулировка расхода и давления осуществляется клапаном 2.

По значениям давления и температуры определялись термодинамические параметры рабочего процесса, на основании которых по формуле (8) проводился расчет изоэнтропийного КПД. Для удобства восприятия в табл. 1 представлены исходные параметры, которые задавались в эксперименте и затем использовались для проведения CFD-расчетов.

РЕЗУЛЬТАТЫ РАСЧЕТОВ И ИХ ОБСУЖДЕНИЕ

В ходе исследования была произведена серия расчетов, результаты которых представлены в табл. 2 и 3. Начальные условия для проведения расчетов выбирались из условий выполнения эксперимента (см. табл. 1). Как видно из табл. 3, различие расчетных и экспериментальных данных не превышает 3%.

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

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

На рис. 5 представлена упрощенная диаграмма процесса детандирования в точке, соответствующей опыту № 7 (см. табл. 1–3), построенная по результатам расчета варианта с наивысшим из представленных значений изоэнтропийного КПД. Расчетное исследование позволяет построить полную диаграмму процесса расширения газа, так как распределения макропараметров известны в любой точке расчетной сетки на каждом шаге по времени, однако результаты экспериментов не содержат данных для промежуточных точек, ввиду чего сравнение диаграмм не представляется возможным.

Методы конечных объемов, в отличие от аналитических методов решения, дают существенно бóльшую наглядность результатов для инженеров, чья специализация не связана с проектированием турбомашин. В частности, интерес вызывают, например, распределения температуры, скоростей и давлений в сечении для установившегося течения, представленные на рис. 6 и 7.

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

В распределении скорости (см. рис. 6, б) также видна неравномерность: поток газа на выходе из рабочего колеса сильно закручен, что приводит к образованию области с высокой скоростью на стенках диффузора, затем поток постепенно замедляется.

По распределению скорости и давления (см. рис. 7) видно, что гелий поступает в направляющий аппарат с высоким давлением и низкой скоростью, однако на выходе из ТДА ситуация меняется: давление газа уменьшается, а скорость многократно возрастает, что доказывает физичность постановки задачи и адекватность расчетной методики.

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

ВЫВОДЫ

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

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

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

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

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

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