МАТЕМАТИЧЕСКАЯ МОДЕЛЬ РАСПРОСТРАНЕНИЯ ПЛЕНОЧНОГО РЕЖИМА КИПЕНИЯ | Опубликовать статью ВАК, elibrary (НЭБ)

Вдовин С.И.¹, Харитонова И.Ю.²

¹Кандидат технических наук, доцент, Нижегородский государственный технический университет им. Р.Е.Алексеева; ²кандидат технических наук, доцент, Нижегородский государственный технический университет им. Р.Е.Алексеева

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ РАСПРОСТРАНЕНИЯ ПЛЕНОЧНОГО РЕЖИМА КИПЕНИЯ

Аннотация

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

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

Vdovin S.I.¹, Kharitonova I.Yu.²

¹Candidate of Sciences in Technics, assosiate professor, Nizhny Novgorod State Technical University n.a. R.E. Alekseev, ²Candidate of Sciences in Technics, assosiate professor, Nizhny Novgorod State Technical University n.a. R.E. Alekseev

MATHEMATICAL MODEL OF THE FILM BOILING DISTRIBUTION

Abstract

In the article a mathematical model is constructed for calculation of the running heat wave that transforms the bubble boiling into the  film boiling on the surface of a fuel element. Control calculations are carried out that allow to set limits on the minimum step in the spatial coordinate.

Keywords: mathematical modeling, boiling regime, numerical algorithm.

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

В работе [2] при изучении явления перехода пузырькового режима кипения в пленочный (на нагреваемых током нитях и стержнях) было замечено, что смена происходит путем возникновения бегущей тепловой волны, меняющей режим кипения на поверхности тепловыделяющего элемента. Экспериментальный и расчетный анализ изменения теплофизических параметров нити накаливания и устойчивости к локальным возмущениям проводился в [3–6].

В настоящей работе, используя разработанный численный алгоритм, исследуется явление формирования и распространения волны пленочного режима кипения применительно к реальным топливным элементам. Анализируется изменение температурного режима и влияние некоторых конструктивных и режимных параметров на скорость рассматриваемых процессов. Эти вопросы представляют интерес при изучении начальных этапов развития аварий, связанных с перегревом твэлов. Сложность расчета обусловлена сильным разнесением характерных пространственных (температурный фронт ~ 10^–3 м, высота АЗ ~ 1 м) и временных (постоянная прогрева горючего ~ 10 с, элемента оболочки твэла ~ 10^–2 с) масштабов.

Описание математической модели и алгоритма расчета

Уравнения теплового баланса для горючего и оболочки твэла в одноточечном приближении их описания в поперечном сечении z запишем в виде:

           (1)

      (2)

где T, λ, ρ, c, V – соответственно температура, теплопроводность, плотность, теплоемкость и объем; индексы _r и _ст соответствуют горючему и оболочке твэла; ψ(z) – пространственная форма тепловыделения; n – амплитуда первой гармоники потока нейтронов тепловой группы; S – площадь поверхности; K – коэффициент теплопередачи от горючего к оболочке твэла; q – плотность теплового потока.

Уравнение (1) не учитывает теплопередачу вдоль горючего. Это оправдано тем, что топливные таблетки разделены зазорами и длина обогреваемой части в 10² ÷ 10³ раз превышает поперечный размер твэла.

Зависимость плотности теплового потока q от разности температур поверхности нагрева и воды ΔT по опытам С.С.Кутателадзе [7] показана на рис.1. Влияние давления, расхода и паросодержания на значения первой и второй критической плотности теплового потока (q_кр1, q_кр2) рассмотрено в работах [8, 9]. Отметим некоторые характерные черты функциональной обусловленности теплового потока. В области пузырькового режима кипения интенсивность теплообмена сильно зависит от температурного напора (ΔT), т.е. использование той либо иной функции от теплофизических параметров теплоносителя не приводит к существенному изменению температурного режима теплоотдающей поверхности. Другая характерная черта состоит в том, что в широком диапазоне изменения давления и расхода (принудительная циркуляция) отношение q_кр1/q_кр2 имеет примерно то же значение, что и при кипении в большом объеме [9]. Это обусловлено гидродинамической природой кризиса теплообмена первого рода. С учетом отмеченных фактов далее при анализе температурного режима в условиях фиксированных значений давления, расхода и энтальпии воды на входе в АЗ будем использовать зависимость теплового потока, показанную на рис. 1.

Рис.1 – Зависимость плотности теплового потока ( q ) от температурного напора (ΔT = T_ст – T): сплошная линия – зависимость по опытным данным С.С.Кутателадзе [7]; пунктирная линия – измененная в диапазоне пленочного кипения зависимость.

 

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

            (3)

u(0, t) = u₀;                    J(0, t) = J₀;                             z ∈ [0, l];           (4)

где u, J, v – скорость, энтальпия и удельный объем теплоносителя;

v’, v’’ – удельный объем воды и пара на линии насыщения;

J’, J’’ – энтальпия воды и пара на линии насыщения;

u₀, J₀, v’, v’’, J’, J’’  – const;

Q = qS/V – тепловой поток в единицу объема теплоносителя. Зависимость температуры воды от энтальпии будем аппроксимировать кусочно-линейной функцией

где a, b – const, определяемые из работы [10] .

Исследование распространения волны пленочного режима кипения проводилось численным методом на основе математической модели (1) – (4). В области резкого (слабого) изменения температуры теплоотдающей поверхности предусматривалось дробление (объединение) шага сетки по пространственной координате. Величина шага ограничивалась сверху и снизу заданными величинами. Для расчета температуры горючего и теплоносителя использовали равномерные сетки, содержащие соответственно n₁ и n₂ узлов. Шаг по времени выбирался с учетом характерной постоянной времени исследуемого переходного процесса, т.е. имел различные значения для горючего, оболочки твэла и теплоносителя.

Для проверки точности определения температуры стенки были проведены контрольные расчеты распространения температурного фронта при постоянном и переменном шаге по пространственной координате. Вычисления проводились на элементе малой длины (l = 0.1 м), температура горючего и теплоносителя не изменялась. В качестве исходного возмущения принималось ступенчатое увеличение температуры оболочки на части рассматриваемого элемента. Величина возмущения задавалась достаточной для образования бегущей тепловой волны при перепаде температуры на стенке в стационарном режиме ΔT* = 22.8^oС. Скорость распространения температурного фронта определялась по смещению точки с фиксированной температурой, несколько превышающей (≈ 40^oC) температуру поверхности, соответствующую пузырьковому режиму кипения.

Показано, что уменьшение постоянного шага разностной сетки с 10^–3 м до 5·10^–4 м увеличивает расчетную скорость распространения температурного фронта с 5.8·10^–3 м/с до 6·10^–3 м/с, т.е. вычисления с шагом менее 10^–3 м приемлемы для контрольного сравнения. При вычислении с переменным шагом аналогичное снижение ограничения на минимальный шаг разностной сетки уменьшает относительную погрешность с 12% (5.3·10^–3 м/с) до значений менее 1% (6·10^–3 м/с). Количество узлов разностной сетки в этом варианте сокращается примерно в 10 раз. Таким образом требуемая точность расчетов изменения температуры оболочки достигается с помощью надлежащего выбора ограничения на минимальный шаг по пространственной координате.