Теплоэнергетика, 2022, № 12, стр. 40-53

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

Рассматривается процесс пленочной конденсации пара на поверхности гладкой горизонтальной трубы внешним диаметром D_out и внутренним диаметром ${{D}_{{in}}}$ (рис. 1). Поток пара, не содержащий примеси, набегает со скоростью U₀ на цилиндр под углом ${{{{\varphi }}}_{0}}$ относительно направления вектора ускорения свободного падения ${\mathbf{g}}.$ На внутренней поверхности трубы задаются необходимые граничные условия – как правило, средняя температура охлаждающей воды ${{T}_{{cw,ave}}}$ и коэффициент теплоотдачи от внутренней стенки к охлаждающей воде ${{{{\alpha }}}_{{cw}}}.$

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Одномерное уравнение сохранения массы для пленки имеет вид [11]

где δ – толщина пленки; t – время; $\xi $ – координата вдоль поверхности трубы, отсчитываемая от верхней образующей в сторону увеличения угла ${{\varphi ;}}$ ${{u}_{f}}$ – средняя скорость движения пленки (f от англ. fluid или film); ρ – плотность конденсата; ${{m}_{f}}$ – плотность потока массы из паровой фазы в жидкую; ${{m}_{{out}}}$ – плотность потока массы отрывающегося от пленки конденсата.

Выражение для связи между объемным расходом конденсата на единицу ширины пленки Г и толщиной пленки можно записать в виде

где ${{{{\tau }}}_{v}}$ – касательное напряжение на межфазной границе (v от англ. vapor); $\mu $ – динамический коэффициент вязкости конденсата; ${{K}_{1}},{{K}_{2}}$ – коэффициенты, учитывающие массообмен [11], в настоящей работе принимается ${{K}_{1}},$ ${{K}_{2}} = 1$ (отличные от единицы значения характерны для жидких металлов [11]); ${{g}_{\xi }}$ – проекция вектора ускорения свободного падения на ось $\xi {\text{.}}$

Плотность потока из паровой фазы в жидкую рассчитывали по соотношению

где ${{\lambda }}$ – коэффициент теплопроводности конденсата; ${{T}_{{sat}}}$ – температура насыщения; ${{T}_{{sur}}}$ – температура внешней поверхности стенки трубы; ${{h}_{{vf}}}$ – теплота фазового перехода (vf от англ. vapor-fluid).

Для определения температуры поверхности стенки трубы решали одномерное (по угловой координате) уравнение для характерной температуры стенки ${{T}_{{w,a}}}$ (w, a от англ. wall average).

где ${{R}_{{in}}} = {{{{D}_{{in}}}} \mathord{\left/ {\vphantom {{{{D}_{{in}}}} 2}} \right. \kern-0em} 2};$ ${{R}_{{out}}} = {{{{D}_{{out}}}} \mathord{\left/ {\vphantom {{{{D}_{{out}}}} 2}} \right. \kern-0em} 2};$ ${{R}_{{med}}} = {{\left( {{{R}_{{in}}} + {{R}_{{out}}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{R}_{{in}}} + {{R}_{{out}}}} \right)} 2}} \right. \kern-0em} 2};$ ${{\lambda }_{w}}$ – коэффициент теплопроводности материала стенки; ${{\delta }_{w}} = \left( {{{R}_{{out}}} - {{R}_{{in}}}} \right);$

– плотность теплового потока на внешней поверхности стенки трубы (от пара к стенке);

– плотность теплового потока на внутренней поверхности стенки; ${{k}_{{out}}}$ – коэффициент теплопередачи от пара к стенке; ${{k}_{{in}}}$ – коэффициент теплопередачи от стенки к воде; ${{{{\alpha }}}_{{cw}}}$ – (cw от англ. cooling water) коэффициент теплоотдачи (определяли по формуле Петухова – Кириллова [12], а температуру воды – по уравнению теплового баланса исходя из заданной температуры воды на входе в охлаждающую трубу).

Изложенная модель пленки не описывает формирование и срыв капель конденсата, поэтому для учета данного эффекта использовали дополнительную упрощенную модель. Считали, что при превышении некоторой критической толщины пленки ${{\delta }_{{cr}}}$ начинался унос конденсата с поверхности трубы. Для моделирования этого процесса использовали линеаризованную форму искусственного источника (стока)

где

– некоторая “большая” величина, отличная от нуля при $\delta > {{\delta }_{{cr}}};$ $\Delta t$ – шаг по времени.

Указанный источник вычисляли во всех контрольных объемах, однако активировали только в контрольном объеме пленки, толщина которой превышала ${{\delta }_{{cr}}} = 10{{\delta }_{0}}$ (здесь ${{\delta }_{0}}$ – характерная толщина пленки из решения Нуссельта). Значение константы 10 было выбрано на основании предварительных тестов. Следует отметить, что изменение этой константы в широком диапазоне (10–25) не оказывало существенного влияния на среднюю теплоотдачу и распределение толщины пленки по периметру трубы.

Для определения касательного напряжения на межфазной границе ${{{{\tau }}}_{v}}$ решали уравнение для безразмерной скорости $F$ [13], к которому может быть сведена стационарная система уравнений пограничного слоя для пара, обтекающего проницаемый цилиндр, с произвольным распределением интенсивности отсоса пара на его поверхности с граничными условиями $F(\xi ,0) = 0$ и $F(\xi ,\infty ) = 1$:

где $F(\xi ,\eta ) = {{u(\xi ,\eta )} \mathord{\left/ {\vphantom {{u(\xi ,\eta )} {{{u}_{\infty }}\left( \xi \right)}}} \right. \kern-0em} {{{u}_{\infty }}\left( \xi \right)}};$ u – скорость движения пара; $\eta $ – расстояние от стенки трубы; ${{u}_{\infty }}\left( \xi \right)$ – заданное распределение скорости вдоль внешней границы пограничного слоя.

Преобразованные координаты $\hat {\eta }{\text{,}}$ $\hat {\xi }$ и параметр $\beta $ рассчитывали по формулам

где ${{\nu }_{v}}$ – кинематический коэффициент вязкости пара.

Безразмерную функцию тока $f$ вычисляли по формуле

где

– значение функции тока на межфазной границе; ${{u}_{{{{\eta }} = 0,v}}}$ – нормальная скорость пара на поверхности цилиндра.

Координату $\xi $ в уравнении (3) отсчитывали от угла ${{\varphi }} = {{{{\varphi }}}_{0}}.$ Считали, что относительно указанного начала отсчета на поверхности цилиндра формировались два пограничных слоя – один в направлении отсчета угла по часовой стрелке, второй при отсчете угла в противоположном направлении. Характеристики каждого слоя рассчитывали из решения (3), при этом необходимое в (4) распределение ${{u}_{{{{\eta }} = 0,v}}}$ определяли из текущего распределения интенсивности конденсации на поверхности трубы.

Распределение скорости на внешней границе пограничного слоя ${{u}_{\infty }}$ вычисляли по соотношению для потенциального течения

где $\Phi = {{\varphi }} - {{{{\varphi }}}_{0}}.$

Для проверки влияния профиля ${{u}_{\infty }}$ на результаты расчета использовали также распределение скорости из работы [14]

ЧИСЛЕННАЯ РЕАЛИЗАЦИЯ

Для решения уравнения (1) периметр трубы разбивали на одномерные контрольные объемы (КО). На концах отрезка (0–2π) задавали периодические граничные условия и использовали явную численную схему. На дискретном уровне уравнение (1) записывается в виде

где

верхним индексом 0 обозначены величины с предыдущего шага по времени, индексом Р – параметры в центре КО.

Переписывая (7) относительно толщины жидкой пленки в центре (узле) КО, можно получить выражение

Для решения уравнения сохранения массы требуется аппроксимация пространственной производной от объемного расхода на единицу ширины пленки. Наиболее устойчивой схемой является схема с разностями “против потока”:

где индексами w, e обозначены параметры на левой и правой гранях КО с узловой точкой P; индексами W, E – значения в узловых точках левого и правого соседних КО; max – оператор, возвращающий максимальное из двух значений.

Результаты тестовых расчетов показали, что данная схема приводит к появлению заметных разрывов в пленке в областях ее “захлебывания”, т.е. в окрестности расчетных точек одномерной сетки, где пленка, увлекаемая внешним потоком пара, начинает двигаться против направления силы тяжести. На грубых сетках появление разрыва можно предотвратить применением схем высокого порядка точности (противопоточная схема второго порядка точности, центрально-разностная схема), однако при измельчении сетки проблема возникает снова, т.е. использование схем высокого порядка не обеспечивает сеточной сходимости решения. Анализ показал, что это связано с особенностями профиля скорости в окрестности данных точек – направления средней скорости движения пленки и скорости движения поверхности пленки при определенных параметрах могут иметь разные знаки.

Выражение для профиля скорости можно представить в виде [11]

где $A = \frac{{{{\tau }_{v}}{{K}_{2}}\delta }}{{{{K}_{1}}\mu {{u}_{s}}}};$ $\bar {\eta } = \frac{\eta }{\delta };$ ${{u}_{s}}$ – скорость на межфазной поверхности.

Координата ${{\eta }_{0}},$ соответствующая u = 0 вне стенки, определяется уравнением

При конденсации следует рассматривать область от стенки до внешней поверхности пленки, т.е. $0 < {{\bar {\eta }}_{0}} < 1.$ Только при $A$ > 2 возможно u = 0 (рис. 2).

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

где ${{\Gamma }_{1}}$ – объемный расход в сечении 0 ≤ ${{\bar {\eta }}}$ < ${{{{\bar {\eta }}}}_{0}};$ Γ₂ – объемный расход в сечении ${{\bar {\eta }}_{0}}$ ≤ $\bar {\eta }$ ≤ 1.

При численной реализации для исключения особенности в точке $A = 4$ в (10) и (11) лучше перейти от средней скорости пленки ${{u}_{f}}$ к скорости на границе ${{u}_{s}},$ используя соотношение

Показатель А в соответствующих узловых точках E, P, W одномерной сетки для пленки можно обозначить символами ${{A}_{E}},{{A}_{P}},{{A}_{W}}.$ Схема “против потока” была модифицирована следующим образом. Если каждое из значений ${{A}_{E}},{{A}_{P}},{{A}_{W}}$ не превышало 2.0, использовали стандартную схему “против потока” (8), (9) для аппроксимации потоков на границе. Если значение одного из этих параметров превышало 2.0, схему изменяли. При ${{A}_{E}} \geqslant 2$ по уравнениям (10), (11) вычисляли ${{\Gamma }_{{1,E}}}$ и ${{\Gamma }_{{2,E}}}$ и расчет ${{\left( {\delta {{u}_{f}}} \right)}_{e}}$ проводили по соотношению

На рис. 3 показаны схема “против потока” и модифицированная схема “против потока” при давлении насыщения ${{p}_{{sat}}} = 0.1\,\,{\text{МПа,}}$ ${{U}_{0}} = 1.97\,\,{\text{м/с,}}$ температуре внешней поверхности стенки трубы 64°C, ${{{{\varphi }}}_{0}} = 60^\circ $ и числе узлов сетки равном 720. Профили скорости вблизи точки растекания конденсата, полученные с помощью модифицированной схемы, показаны на рис. 4.

Уравнение теплопроводности (2) дискретизировали с использованием центрально-разностной схемы второго порядка. Полученную систему линейных алгебраических уравнений на каждом шаге по времени решали методом прогонки. Решение уравнения (3) можно свести к последовательному решению системы обыкновенных дифференциальных уравнений с помощью встроенной функции библиотеки SciPy [15]. Для нахождения производных по координате $\xi $ использовали аппроксимацию с разностями “назад” с адаптивным шагом, сгущающимся при приближении к точке отрыва пограничного слоя. Уравнение (3) решали на каждом n-м шаге по времени для пленки. При некотором значении n полученное квазистационарное решение переставало зависеть от него.

ТЕСТИРОВАНИЕ МОДЕЛИ ВНЕШНЕГО ПОТОКА

Для оценки точности определения распределения касательных напряжений, получаемых при решении (3), была выполнена серия тестовых CFD-расчетов обтекания цилиндра с заданной скоростью отсоса пара на его поверхности. Схема расчетной области и фрагмент численной сетки изображены на рис. 5. С использованием СFD-кода ANES [16] решали двумерную нестационарную систему уравнений Навье – Стокса, осредненную по Рейнольдсу. Для расчета турбулентных характеристик потока пара использовали двухпараметрическую SST-модель c универсальными пристеночными функциями Ментера [17]. Касательная компонента скорости на поверхности цилиндра была равна нулю, нормальную компоненту скорости задавали с использованием параметра $j = {{{{u}_{{{{\eta }} = 0,v}}}} \mathord{\left/ {\vphantom {{{{u}_{{{{\eta }} = 0,v}}}} {{{U}_{0}}}}} \right. \kern-0em} {{{U}_{0}}}}.$

На рис. 6 представлены результаты тестовых CFD-расчетов и экспериментальные данные [18] по распределению безразмерного градиента скорости пара на стенке до точки отрыва пограничного слоя, на рис. 7 – по распределению коэффициента давления ${{С}_{p}} = 2{{\left( {p - {{p}_{0}}} \right)} \mathord{\left/ {\vphantom {{\left( {p - {{p}_{0}}} \right)} {\left( {{{\rho }_{v}}U_{o}^{2}} \right)}}} \right. \kern-0em} {\left( {{{\rho }_{v}}U_{o}^{2}} \right)}}$ по поверхности цилиндра [19]. Распределения расчетных величин, представленные на рис. 6, 7, получены путем осреднения по времени данных нестационарного CFD-расчета.

На рис. 8 приведены результаты CFD-моделирования и расчета с использованием уравнения (3) распределений ${{c}_{f}}\operatorname{Re} _{v}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}$ (здесь c_f – коэффициент трения) по поверхности цилиндра при задании скорости на внешней границе пограничного слоя по уравнению (5). Видно, что с ростом параметра  j согласование результатов расчетов в приближении пограничного слоя и CFD-моделирования становится лучше. При низкой интенсивности отсоса пара лучшее согласие достигается при использовании экспериментально определенного распределения скорости (6), что вполне ожидаемо, так как оно было получено при  j = 0.

Следует отметить, что при обдуве цилиндра сверху именно при $\Phi $ < 60° интенсивность конденсации (плотность отсоса пара) максимальна. В этом диапазоне $\Phi $ значения u_∞, рассчитанные по (5) и (6), практически одинаковы и решения уравнения пограничного слоя (3) хорошо воспроизводят результаты CFD-моделирования.

Интерес представляют режимы обтекания цилиндра при неоднородном по периметру отсосе газа, что значительно лучше (особенно в зоне отрыва пограничного слоя) отвечает привлекаемой аналогии между отсосом газа на поверхности цилиндра и массовым потоком пара к поверхности жидкой пленки вследствие конденсации. На рис. 9 представлены результаты решения уравнений пограничного слоя и CFD-моделирования для некоторого “предельного” варианта неоднородного отсоса пара.

В целом можно сделать вывод, что расчет по (3) дает достаточно разумные результаты.

Описанная математическая модель схожа с предложенной в работе [20]. Главное различие моделей заключается в том, что авторы [20] не учитывали влияние массообмена на касательные напряжения и решали задачу в стационарной постановке, что привело к необходимости использовать метод “проб и ошибок” для подбора толщины пленки в произвольной точке цилиндра. Предложенная в настоящей работе модель пленки является нестационарной и для рассмотренных режимов обеспечивает выход на квазистационарный режим конденсации, при котором интеграл по поверхности пленки от правой части уравнения (1) дает значение, близкое к нулю, при этом интегралы от отдельных составляющих правой части практически не изменяют своего значения во времени.

РЕЗУЛЬТАТЫ ВАЛИДАЦИИ МОДЕЛИ

Данные, отобранные для валидации, представлены в табл. 1. Все расчеты были выполнены с использованием (5). В отдельных случаях для дополнительного тестирования причин рассогласования были проведены расчеты с использованием (6) или упрощенной модели (1), интегрированной в CFD-код.

На рис. 10 представлена зависимость экспериментального теплового потока ${{q}_{{w,exp}}}$ от рассчитанного ${{q}_{{w,calc}}}$ (здесь q_w – тепловой поток на стенке трубы) для пара, движущегося сверху вниз. В табл. 2 приведена информация по отклонению ${{q}_{{w,calc}}}$ от данных [21, 22] и среднему отклонению ${{\Delta }_{{mean}}} = \frac{1}{N}\sum\nolimits_N {{{\left| {\left( {{{q}_{{w,exp}}} - {{q}_{{w,calc}}}} \right)} \right|} \mathord{\left/ {\vphantom {{\left| {\left( {{{q}_{{w,exp}}} - {{q}_{{w,calc}}}} \right)} \right|} {{{q}_{{w,exp}}}}}} \right. \kern-0em} {{{q}_{{w,exp}}}}}} $ (здесь N – количество точек).

Была выполнена серия расчетов конденсации пара без учета влияния скорости, т.е. для случая ${{{{\tau }}}_{v}} = 0.$ На рис. 11 представлена зависимость отношения экспериментального коэффициента теплопередачи α_exp к расчетному α_calc от динамического напора ${{{{\rho }_{v}}U_{0}^{2}} \mathord{\left/ {\vphantom {{{{\rho }_{v}}U_{0}^{2}} 2}} \right. \kern-0em} 2}$ (здесь ρ_v – плотность пара) для части данных из работы [22], для которых была приведена информация о температуре стенки. Коэффициент теплоотдачи вычисляли по формуле

где $\left\langle \ldots \right\rangle $ означает температуру поверхности, осредненную по периметру трубы.

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

Авторы [23] предоставили детальную информацию для тринадцати режимов конденсации на внешней поверхности трубы движущегося горизонтально водяного пара. Данные представлены группами режимов с примерно близкими параметрами и одним режимом с низкой скоростью пара (табл. 3). На рис. 12 показана зависимость отношения рассчитанной плотности теплового потока на стенке ${{q}_{{w,calc}}}$ к экспериментально определенному значению ${{q}_{{w,exp}}}$ от динамического напора. Наибольшее рассогласование между экспериментальными и расчетными данными наблюдается для режимов из первой группы, отличающихся наибольшей скоростью набегающего потока.

В [23] представлены данные о распределении температуры стенки по периметру для всех режимов, при этом все распределения для одной группы обозначены одинаковыми символами, поэтому выделить данные о распределении температуры стенки для конкретного режима в каждой группе не представляется возможным. На рис. 13 для каждой группы показано распределение температуры θ, которую вычисляли по формуле

При этом считали, что рассчитанная по уравнению (2) характерная температура стенки ${{T}_{{w,a}}}$ и температура в местах закладки термопар совпадают вследствие малого термического сопротивления стенки. Угол ${{\gamma }}$ определяли как ${{\gamma }} = - \left( {{{\varphi }} - {{{{\varphi }}}_{0}}} \right)$ в диапазоне от 0 до 360°. Расчетные распределения для всех режимов одной группы практически не различимы между собой, и поэтому на рис. 13 показана одна представительная линия. Наблюдается заметное рассогласование экспериментальных и расчетных данных при ${{\gamma }}$ = 120–240°. К сожалению, сравнение размерных локальных температур не представляется возможным из-за того, что в работе [23] отсутствует информация о соответствии между конкретными распределениями безразмерных температур и номерами режимов внутри каждой группы экспериментальных данных. Следует отметить, что рассогласование безразмерной температуры при некотором угле ${{\gamma }}$ не означает несовпадения размерных температур в этой же точке, так как распределения безразмерной температуры (13) зависят от значения знаменателя в (13).

Для дополнительной валидации были взяты экспериментальные данные (${{{{\varphi }}}_{0}} = 90^\circ $) из работ [20, 24]. К сожалению, в указанных работах не приводится подробная информация об условиях охлаждения трубы, а представленные данные содержат ряд противоречий. В частности, значение ${{q}_{{w,exp}}}$ из баланса тепла охлаждающей воды, используемое при определении экспериментального значения коэффициента теплоотдачи (12), существенно отличается от значения ${{q}_{{w,exp}}},$ полученного осреднением локальной плотности теплового потока в [24]. Определенное авторами [20, 24] значение коэффициента теплоотдачи (12) при значительно неоднородной температуре стенки существенно отличается от значений, рассчитанных по толщине пленки $\alpha {\kern 1pt} {\text{'}} = \left\langle {{\lambda \mathord{\left/ {\vphantom {\lambda \delta }} \right. \kern-0em} \delta }} \right\rangle ,$ при этом с ростом скорости пара данное отличие усиливается. Поэтому сравнение выполняли только с первичными экспериментальными данными о распределении температуры стенки. При этом по аналогии с расчетами авторов [20] температуру охлаждающей воды принимали равной 15°С, а коэффициент теплоотдачи к охлаждающей воде подбирали таким образом, чтобы обеспечить экспериментально полученный перепад между температурой насыщения и температурой стенки в передней критической точке. Считали также, что рассчитанная характерная температура стенки ${{T}_{{w,a}}}$ и температура в местах закладки термопар совпадают.

Анализировали результаты трех вариантов использования упрощенной модели, отличающиеся расчетом касательных напряжений на поверхности жидкой пленки: из решения уравнения (3) с граничными условиями (5) (вариант а) и (6) (вариант б) и из решения двумерных уравнений гидродинамики с использованием CFD-кода ANES (вариант в).

На рис. 14 представлено распределение безразмерной температуры стенки $\vartheta $ = (T_w,a – T_c,w)/(T_sat – – T_cw) по угловой координате. Можно видеть хорошее согласие рассчитанных и экспериментальных профилей [20, 24], при этом способ расчета касательных напряжений не приводит к существенным различиям результатов. Из-за того что параметр  j для рассматриваемых режимов мал, вариант (б) дает несколько лучшее соответствие экспериментальным данным по сравнению с вариантом (а).

Для некоторых режимов, полученных с использованием CFD-моделирования, можно отметить появление немонотонностей в распределении ϑ в окрестности ${{\gamma }}$ = 180°. Данная особенность связана с тем, что вихри в ближнем следе за цилиндром (рис. 15) начинают активно взаимодействовать с пленкой конденсата. Для указанных режимов часто используемое в различных моделях [20, 24–26] предположение о том, что после отрыва пограничного слоя пленка движется под действием только силы тяжести, оказывается несправедливым.

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

ВЫВОДЫ

1. Результаты валидации позволяют рекомендовать разработанную авторами математическую модель для применения в CFD-кодах для расчета процессов конденсации пара на поверхности горизонтальных труб при произвольном направлении движения внешнего потока пара и умеренных динамических напорах.

2. Предложенные аппроксимации объемного потока на гранях одномерных контрольных объемов с учетом особенностей профиля скорости в пленке позволяют избежать разрывов в толщине пленки в окрестности точки “захлебывания”.

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

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

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