Model for integral codes for calculating containment thermohydraulics in severe accidents at nuclear power plants with a water coolant

Texto integral

Введение

Для моделирования тяжелых аварий (ТА) на АЭС с РУ ВВЭР в России разработан код СОКРАТ в версиях В1/В2 и В3 [1], широко используемый в задачах анализа безопасности АЭС. При авариях с потерей теплоносителя динамика протекания внекорпусной, а при больших течах и внутрикорпусной фазы ТА, существенно зависят от параметров в ЗО.

Исторически теплогидравлика ЗО моделировалась с помощью отдельных контейнментных кодов АНГАР и КУПОЛ. На начальном этапе развития кода СОКРАТ давление и состав газовой среды в ЗО, как граничные условия для кода СОКРАТ, рассчитывались итерационно путем последовательного расчета источника массы и энергии в ЗО по коду СОКРАТ и запуска автономного расчета по контейнментному коду с целью получения граничных условий приемлемого уровня точности. Позднее контейнментые коды стали связываться с кодом СОКРАТ через универсальный интерфейс, что позволило выполнять самосогласованный расчет. Однако такой подход все же обладал рядом недостатков как технического (например, использование различных методик расчета свойств парогазовой среды, согласование шагов интегрирования), так и административного характера (например, направленность валидации контейнментных кодов на проектные режимы, рассогласованные сроки действия аттестационных паспортов и фактическая зависимость использования кода СОКРАТ от решений организаций-разработчиков контейнментных кодов по повторной аттестации), что указывало на необходимость наличия в коде СОКРАТ собственной модели для расчета параметров в ЗО. В окончательной формулировке задача была расширена на возможность моделировать медленное горение водорода и перенос радиоактивных веществ в ЗО для выполнения расчетов в обоснование безопасности АЭС, включая детерминистические расчеты в поддержку ВАБ-2. Последнее накладывает требования высокого быстродействия и возможностей описывать весь спектр систем безопасности в ЗО, используемых на действующих и проектируемых АЭС для управления и смягчения последствий ТА.

Ниже приведено описание модели, реализованной в виде модуля CONT_TH для кода СОКРАТ. Текущее модульное наполнение кода СОКРАТ позволяет описывать спектр тяжелых аварий при автоматическом учете взаимосвязи процессов в РУ и ЗО.

Общая характеристика модели

Модель предназначена для расчета давления, температуры и состава парогазовой среды, температур воды и поверхности стен в каждом из моделируемых помещений ЗО, массы воды в приямках, скорости парогазовой среды в гидравлических связях между помещениями и в местах утечек из помещений ЗО в окружающую среду в каждый момент времени.

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

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

•  истечение двухфазного многокомпонентного теплоносителя из контуров РУ в помещения и обратное течение газовой фазы;
•  перенос парогазовой среды между помещениями;
•  конденсацию водяного пара на стенках и в объеме;
•  теплообмен между стенами и парогазовой средой;
•  теплопроводность стен, отвод тепла в окружающую среду через стены;
•  работу спринклерной системы (тепло и массообмен среды с каплями);
•  работу рекомбинаторов водорода;
•  поступление конденсата в приямки и его испарение, а также перетечки между приямками.

В основе моделей лежат следующие допущения:

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

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

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

Применение модуля CONT_TH по свойствам жидкой и газовой фазы ограничивается следующим диапазоном термодинамических параметров:

•  давление от 0.01 до 27.5 МПа;
•  массовая концентрация неконденсируемого газа от 0 до 1;
•  температура газовой фазы от 270 до 3300 К;
•  температура жидкой фазы от 273.16 до температуры насыщения +10 К (но не выше 635 К).

Применение модуля для задач поверхностного кипения воды и конденсации пара ограничивается по давлению диапазоном от 0.01 до 22 МПа.

Конечно-разностная линеаризованная система алгебраических уравнений решается методом Гаусса.

Численное интегрирование уравнения теплопроводности для стен проводится по неявной конечно-разностной схеме с использованием метода прогонки. Уравнения тепло- и массообмена спринклерных капель, записанные относительно координаты вдоль линии движения, линеаризуются и решаются в явном виде с приращением скорости, температуры или массы не более чем 0.5% за внутренний шаг интегрирования.

В составе ПрЭВМ СОКРАТ-В1/В2 модуль CONT_TH аттестован для расчета давления и мольных долей неконденсируемых газов в помещениях герметичных ограждений АЭС с РУ ВВЭР [2].

Методика построения расчетных схем

В модуле CONT_TH имеются следующие модели элементов и систем:

•  элементы гидравлической сети:
  
  •  помещение (ROOM) – гидравлический элемент с сосредоточенными параметрами:
  •  источник массы и энергии (INJECTION) – моделирует объемный источник/сток массы и соответствующей ей энергии в помещениях;
  •  тепловыделение (HEAT_SOURCE) – моделирует объемный источник тепла в помещениях;
  •  спринклер (SPRAY) – моделирует работу спринклерной системы;
  •  рекомбинатор водорода (RECOMB) – моделирует работу ПАР;
  •  связь (CONNECTION) – соединение двух помещений по газовой фазе;
  •  вентилятор (FAN) – обеспечивает заданный расход среды в связи;
  •  клапан (VALVE) – используется в связи для изменения проходного сечения;
  •  дренаж (DRAIN) – моделирует сток воды за счет действия гравитационных сил;
•  условия на границах системы (CTH_BOUNDARY):
  
  •  вывод/ввод теплоносителя с различным расходом;
  •  заданное давление, температура и состав газовой фазы;
•  тепловые элементы (WALL):
  
  •  плоская стенка с теплообменом с теплоносителем или с заданными параметрами на двух границах;
  •  система контроля и управления:
  •  датчик;
  •  триггер;
  •  регулятор;
  •  контрольный блок.

Элемент “помещение” (ROOM) используется для описания поведения теплоносителя в одном или группе смежных помещений защитной оболочки. На рис. 1 приведен пример нодализационной схемы, в которой помещение А разделено по высоте на два расчетных объема А1 и А2 для учета возможной стратификации газовой фазы. Помещение характеризуется объемом, относительной высотой центра масс (в этой точке определяется давление), характерной длиной, гидравлическим диаметром и площадью проходного сечения. В каждом помещении пользователем задается давление и состав газовой фазы, ее объемная доля, а также температуры газовой и жидкой фаз. В дальнейшем они рассчитываются на каждом шаге и могут быть использованы в системе контроля и управления. Характерным размером помещения в общем случае следует выбирать высоту, а проходным сечением – отношение объема к гидравлическому диаметру. В помещении может задаваться один или несколько источников INJECTION двухфазного теплоносителя (из контуров РУ или других систем). Для источника пара и воды используется энтальпия, а для неконденсируемого газа – температура. Для передачи импульса в газовую фазу используется проходное сечение источника (с учетом возможного обжатия струи). При помощи источника возможно независимо отбирать воду, а также каждый из компонентов газовой фазы, при этом удельная энтальпия отбираемого компонента будет соответствовать текущему значению в данном помещении. Для изменения энергосодержания каждой из фаз можно использовать объемное тепловыделение (HEAT_SOURCE).

Рис. 1. Пример нодализационной схемы для модуля CONT_TH.

Спринклерная система может работать в нескольких помещениях (на рис. 1 это А1 и А2), она характеризуется:

•  количеством форсунок;
•  высотой установки форсунки;
•  углом распыла;
•  эквивалентным диаметром сопла (или диаметром сопла и коэффициентом истечения);
•  расходом;
•  эквивалентным диаметром капель (диаметр Саутера или D₃₂);
•  для каждого помещения высотой пролета капли и долей расхода, попадающей на стены.

Для рекомбинаторов водорода задается производительность в виде функции необходимых параметров системы контроля и управления или выбирается тип, для которого производительность уже определена.

Помещения по газовой фазе соединяются при помощи связей (CONNECTION), не содержащих массы и характеризующихся (кроме проходного сечения) относительной высотой для расчета гидростатического напора со стороны каждого из помещений: длиной, гидравлическим диаметром и локальным сопротивлением для расчета сил трения. Скорость течения газа в связи ограничена скоростью звука в среде донорного помещения. Все достаточно большие проемы, в которых может существовать противоточное течение газа, рекомендуется задавать двойной связью с сохранением всех характеристик кроме уменьшенного вдвое проходного сечения и разнонаправленности. Разнонаправленность вертикальных связей (между помещениями А1 и А2 на рис. 1) позволяет через поправки к уравнению импульса моделировать обмен массой между помещениями в случае, если внизу находится среда с меньшей плотностью или расчетная среднеобъемная скорость хотя бы в одном из помещений отлична от 0. В случае, когда температура в помещении А1 больше, чем в помещении А2, автоматически рассчитывается тепловой поток между ними (без массопереноса). Массообмен между помещениями А1 и В через двойную связь будет осуществляться только в случае разных плотностей газовой фазы между помещениями за счет гидростатического напора. Отметим, что в случае задания одной связи (при подключении граничных условий) противоток среды смоделировать невозможно из-за модельных допущений об одинаковой скорости для каждого компонента газовой фазы. Относительные высотные отметки двойных горизонтальных связей для прямоугольного проема рекомендуется устанавливать на ¼ и ¾ высоты, а для более общего случая – на высотах, ниже которых расположены ¼ и ¾ проходного сечения.

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

Интенсивность дренирования конденсата (DRAIN) осуществляется с постоянным расходом или может быть задана, например, как зависимость от уровня воды в помещении. Под конденсатом понимается вода, образовавшаяся за счет процессов объемной или поверхностной конденсации, а также поступившая из источника. Перенос спринклерных капель между помещениями рассчитывается в модели SPRAY. Уровень воды в приямке (помещение А1 на рис. 1) рассчитывается пользователем на основании массы воды и поперечного сечения помещения.

Система уравнений для теплоносителя

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

В основе модели лежит подход, схожий с используемым в модели контурной теплогидравлики ПрЭВМ СОКРАТ и многих других расчетных средств, аналогичных по функционалу. Система уравнений (1)÷(6) описывает термодинамические параметры двухфазной среды с примесью неконденсируемых газов и представляет собой систему из 5+N дифференциальных уравнений в частных производных с алгебраическими замыкающими соотношениями. Основными уравнениями модели являются уравнения неразрывности и энергии, для газовой фазы используется уравнение движения. Для численного решения в качестве основных независимых переменных выбраны следующие величины: ${\alpha }_g,P,h_g,h_f,V_g,X_1\mathrm{,..,}{X}_N$.

Уравнение неразрывности газовой фазы:

$\frac{\partial }{\partial t}\left({\alpha }_g{\rho }_g\right)+\frac{1}{Vol}\sum _m{\left(A{\rho }_g{V}_g\right)}_m={\Gamma }_i+{\Gamma }_w+{\Gamma }_{spr}+\sum _{n=1}^N{S}_n+S_v$. (1)

Уравнения неразрывности неконденсируемых газов (n = 1…N):

$\frac{\partial }{\partial t}\left({\alpha }_g{\rho }_g{X}_n\right)+\frac{1}{Vol}\sum _m{\left(A{\rho }_g{X}_n{V}_g\right)}_m=S_n$. (2)

Уравнение неразрывности жидкой фазы:

$\frac{\partial }{\partial t}\left({\alpha }_f{\rho }_f\right)=S_f-{\Gamma }_i-{\Gamma }_w-{\Gamma }_{spr}$. (3)

Уравнение энергии газовой фазы записано относительно энтальпии (с учетом уравнений массы и импульса как предложено, например, в [3]):

$\begin{array}{c}{\alpha }_g{\rho }_g\frac{\partial h_g}{\partial t}+\sum _m{\left(\frac{A{\rho }_g{h}_g{V}_g}{Vol}\right)}_m-\frac{h_g}{Vol}\underset{m}{}{\left(A{\rho }_g{V}_g\right)}_m-{\alpha }_g\frac{\partial P}{\partial t} =\\ =({\Gamma }_i+{\Gamma }_{spr})\left(h_{iv}-h_g\right)+{\Gamma }_w\left(h_{wv}-h_g\right)+\left({\chi }_w+{\chi }_{local}\right)V_g +\\ +Q_{iv}+Q_{gf}+Q_{wg}+Q_g+Q_{rec}+Q_{spr}+\sum _{n=1}^N{S}_n\left(h_n^{\text{'}}-h_g\right)+S_v\left(h_v^{\text{'}}-h_g\right)\end{array}$. (4)

Уравнение энергии жидкой фазы записано аналогично: ${\alpha }_f{\rho }_f\frac{\partial h_f}{\partial t}-{\alpha }_f\frac{\partial P}{\partial t}=-({\Gamma }_i+{\Gamma }_{spr})\left(h_{if}-h_f\right)-{\Gamma }_w\left(h_{wf}-h_f\right)+ S_f\left(h_f^{\text{'}}-h_f\right)+$

+ Q_if – Q_gf + Q_wf + Q_f – Q_spr. (5)

Уравнение движения газовой фазы:

${\left({\rho }_g\right)}_m\frac{\partial {\left(V_g\right)}_m}{\partial t}+\frac{\Delta P}{\Delta z}={\chi }_w+{\chi }_{local}-g{\left({\rho }_g\frac{\Delta H}{\Delta z}\right)}_m$. (6)

В общем случае система содержит (5+N)·max(K_j,K_m) уравнений. Здесь K_j – число расчетных помещений включая граничные условия (ГУ), а K_m – число связей между помещениями или помещением и ГУ.

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

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

•  модель теплоотдачи от парогазовой среды к стенкам для расчета Q_wg, Q_wf;
•  модель объемной конденсации пара для расчета Г_i;
•  модель поверхностной конденсация пара в присутствии неконденсируемых газов для расчета Г_w;
•  модель спринклерной системы для расчета Q_spr, Г_spr;
•  модель рекомбинатора водорода для расчета интенсивности рекомбинации и Q_rec (может использоваться для расчета как рекомбинаторов, так и горения водорода в адиабатическом изохорном приближении);
•  модель функционирования вентиляционной системы;
•  модель источников массы и энергии;
•  модель перетечек жидкой фазы.

Линеаризация по времени

Для аппроксимации (линеаризации) уравнений течения теплоносителя (1)÷(6) по времени используется разностная схема, в основе которой положен метод SIMPLE (Semi-Implicit Method for Pressure Linked Equations) [4], когда одновременно решается система для поля давлений и поля скоростей. Эта система изначально предназначена для снижения затрат процессорного времени для CFD кодов. Впоследствии она получила достаточно широкое распространение. Для одномерного случая течения газового теплоносителя этот метод достаточно подробно изложен в работе [5] по описанию полунеявной одномерной модели течения сжимаемого газа. Для описываемой модели, в отличие от работы [5], используется энтальпия вместо температуры (как и в работе [4]), что обусловлено необходимостью учета метастабильных состояний. Кроме этого, добавлены уравнения неразрывности неконденсируемых газов:

${\left({\alpha }_g{\rho }_g\right)}^{*}-\left({\alpha }_g{\rho }_g\right)+\tau \frac{1}{Vol}\sum _m{\left(A{\rho }_g{V}_g^{*}\right)}_m=\tau {\Gamma }_i^{*}+\tau {\Gamma }_w+\tau {\Gamma }_{spr}+\tau \sum _{n=1}^N{S}_n+\tau S_v$. (7)

${\left({\alpha }_g{\rho }_g{X}_n\right)}^{\text{*}}-\left({\alpha }_g{\rho }_g{X}_n\right)+\tau \frac{1}{Vol}\sum _m{\left(A{\rho }_g{V}_g^{*}{X}_n\right)}_m=\tau \sum _{n=1}^N{S}_n$. (8)

${\left({\alpha }_f{\rho }_f\right)}^{*}-\left({\alpha }_f{\rho }_f\right)=\tau {\Gamma }_i^{*}+\tau S_f-\tau {\Gamma }_w-\tau {\Gamma }_{spr}$. (9)

 $\begin{array}{c}\frac{{\alpha }_g{\rho }_g}{\tau }\left(h_g^{*}-h_g\right)+\sum _m{\left(\frac{A{\rho }_g{h}_g^{*}{V}_g}{Vol}\right)}_m-\frac{h_g^{*}}{Vol}\sum _m{\left(A{\rho }_g{V}_g\right)}_m -{\alpha }_g\left(\frac{P^{*}-P}{\tau }+\sum _m{\left(\frac{A{P}^{*}{V}_g}{Vol}\right)}_m\right)=\\ =\left({\Gamma }_i^{*}+{\Gamma }_w+{\Gamma }_{spr}\right)\left(h_{iv}-h_g\right)+\sum _{n=1}^N{S}_n\left(h_n^{\text{'}}-h_g^{*}\right)+S_v\left(h_v^{\text{'}}-h_g^{*}\right)+\frac{P_v}{P}{\eta }_{iv}\left(T_{sv}^{*}-T_g^{*}\right)+\\ +\frac{P-P_v}{P}{\eta }_{gf}\left(T_f^{*}-T_g^{*}\right)+\left(k_w+k_{local}\right)\left|V_g\right|V_g^{*}+Q_{wg}+Q_g+Q_{rec}+Q_{spr}\end{array}$(10)

$\frac{{\alpha }_f{\rho }_f}{\tau }\left(h_f^{*}-h_f\right)-{\alpha }_f\left(\frac{P^{*}-P}{\tau }\right)=-\left({\Gamma }_i^{*}+{\Gamma }_w+{\Gamma }_{spr}\right)\left(h_{if}-h_f\right)+$

$+{\eta }_{if}\left(T_{sv}^{*}-T_f^{*}\right)+S_f\left(h_f^{\text{'}}-h_g^{*}\right)+\frac{P-P_v}{P}{\eta }_{gf}\left(T_g^{*}-T_f^{*}\right)+Q_{wf}+Q_f-Q_{spr}$ (11)

$\frac{{\rho }_g}{\tau }\left(V_g^{\text{*}}-V_g\right)+\frac{\partial }{\partial z}{P}_g^{\text{*}}=\left(k_w+k_{local}\right)V_g^{*}-{\rho }_{g}g\frac{\partial H_m}{\partial z}.$ (12)

Нелинейные члены в уравнениях неразрывности (8) аппроксимируются следующим образом:

$X_n^{*}{\alpha }_g^{*}=X_n{\alpha }_g^{*}+X_n^{*}{\alpha }_g-X_n{\alpha }_g$. (13)

Интенсивность массообмена на верхнем временном слое представляется как:

 ${\Gamma }_i^{*}={\Gamma }_i+\frac{\partial {\Gamma }_i}{\partial F}\left(F^{*}-F\right),$(14)

где $F=\left(P,{\alpha }_g,h_g,V_g,X_1,\dots X_N\right)$  вектор-функция независимых переменных.

Плотность газа и температуры на верхнем временном слое определяются через независимые переменные:

${\rho }_g^{*}={\rho }_g+\frac{\partial {\rho }_g}{\partial P}\left(P^{*}-P\right)+\frac{\partial {\rho }_g}{\partial h_g}\left(h_g^{*}-h_g\right)+\sum _{n=1}^N\frac{\partial {\rho }_g}{\partial X_n}\left(X_n^{*}-X_n\right)$, (15)

$T_g^{*}=T_g+\frac{\partial T_g}{\partial P}\left(P^{*}-P\right)+\frac{\partial T_g}{\partial h_g}\left(h_g^{*}-h_g\right)+\sum _{n=1}^N\frac{\partial T_g}{\partial X_n}\left(X_n^{*}-X_n\right)$, (16)

$T_{vs}^{*}=T_{vs}+\frac{\partial T_{vs}}{\partial P}\left(P^{*}-P\right)$, (17)

$T_f^{*}=T_f+\frac{\partial T_f}{\partial P}\left(P^{*}-P\right)+\frac{\partial T_f}{\partial h_f}\left(h_f^{*}-h_f\right).$ (18)

В выражениях (15)÷(18) температуры и их производные определяются из уравнений состояния. Для воды ${\text{T}}_{\text{f}}={\text{T}}_{\text{f}}\left(\text{P},{\text{h}}_{\text{f}}\right)$, расчет параметров газовой фазы (ρ_g, T_g, ρ_v, P_v) описан ниже.

Пространственная аппроксимация

Пространственная аппроксимация уравнений неразрывности и энергии основана на ”донорном“ методе определения направления потока по состоянию на начало расчетного шага. В результате члены уравнений (7)÷(8) описывающие перетечки среды с учетом (13), преобразуются к виду:

 $\frac{1}{Vol}\sum _m{\left(A{\rho }_g{V}_g^{*}\right)}_m\to \frac{1}{Vol}\sum _m{A}_m{\dot{\rho }}_{gm}\left(V_{gm}^{*}+\left({\dot{\dot{\alpha }}}_{gm}^{*}/{\dot{\alpha }}_{gm}-1\right)V_{gm}\right)$

$\frac{1}{Vol}\sum _m{\left(A{\rho }_g{V}_g^{*}{X}_n\right)}_m\to \frac{1}{Vol}\sum _m{A}_m{\dot{\rho }}_{gm}{\dot{X}}_{nm}\left(V_{gm}^{*}+\left({\dot{\dot{\alpha }}}_{gm}^{*}/{\dot{\alpha }}_{gm}-1\right)V_{gm}\right)$(19)

В левой части уравнении сохранения энергии газовой фазы (10), записанного для помещения j индекс k соседнего помещения определяется из топологии связи m и является зависимой величиной: k=k(m). Суммирование проведено отдельно для втекающих и вытекающих потоков газа:

$\begin{array}{c}\frac{{\alpha }_g{\rho }_g}{\tau }\left(h_g^{*}-h_g\right)+\sum _m{\left(\frac{A{\rho }_g{h}_g^{*}{V}_g}{Vol}\right)}_m-\frac{h_g^{*}}{Vol}\sum _m{\left(A{\rho }_g{V}_g\right)}_m-{\alpha }_g\left(\frac{P^{*}-P}{\tau }+\sum _m{\left(\frac{A{P}^{*}{V}_g}{Vol}\right)}_m\right)\to \\ \to \frac{{\alpha }_{gj}{\rho }_{gj}}{\tau }\left(h_{gj}^{*}-h_{gj}\right)-{\alpha }_{gj}\left(\frac{P_j^{*}-P_j}{\tau }\right)+\sum _m{\left.\frac{A_m{V}_{gm}({\rho }_{gk}{h}_{gk}^{*}-P_k^{*})}{Vo{l}_k}\right|}_{V_{gm\left(k\to j\right)>0}}+\\ +\sum _m{\left.\frac{A_m{V}_{gm}({\rho }_{gj}{h}_{gj}^{*}-P_j^{*})}{Vo{l}_j}\right|}_{V_{gm\left(k\to j\right)<0}}\end{array}$. (20)

Уравнение сохранения импульса (12) в связи m между помещениями j и k принимает вид:

$\begin{array}{c}\frac{{\rho }_j\Delta z_k+{\rho }_k\Delta z_j}{\Delta z_k+\Delta z_j}\left(V_{gm}^{\text{*}}-V_{gm}\right)+\tau \frac{P_k^{*}-P_j^{*}}{\Delta z_k+\Delta z_j}+\tau \frac{\Delta P_{mom}}{\Delta z_k+\Delta z_j}+\tau \frac{\Delta P_{den}}{\Delta z_k+\Delta z_j}=\\ =\tau \left(k_w+k_{local}\right)V_g^{*}-\tau g\frac{\left({\rho }_{gk}\left(H_k-H_m\right)+{\rho }_{gj}\left(H_m-H_j\right)\right)}{\Delta z_k+\Delta z_j}\end{array}$. (21)

Величина $\Delta P_{mom}$ пропорциональна ${\rho }_g{V}_g^2$ и представляет собой движущую силу, обеспечивающую турбулентное перемешивание среды между двумя расчетными помещениями. Для величины $\Delta P_{den}$, представляющей собой движущую силу всплытия легкого газа и пропорциональной $g\left({\rho }_{gj}-{\rho }_{gk}\right)\left(H_j-H_k\right)$, используются только положительные значения. При наличии двух разнонаправленных связей между помещениями (задаются пользователем) $\Delta P_{mom}$ и/или $\Delta P_{den}$ обеспечивают перемешивание газовой фазы между этими помещениями. Выше приведен пример построения расчетной схемы, учитывающей перемешивание среды за счет указанных выше сил даже в случае гидростатического равновесия.

Перенос тепла в конструкционных элементах

Одномерное уравнение теплопроводности (22) аппроксимируется на плоской одномерной сетке:

${с}_w{\rho }_w\frac{\partial T}{\partial t}=\frac{\partial }{\partial x}\left({\lambda }_w\frac{\partial T}{\partial x}\right)+Q\left(x,t\right)$ . (22)

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

${с}_i{\rho }_i\frac{\partial T}{\partial t}=\frac{2}{x_{i+1}-x_{i-1}}\left({\lambda }_{i+\frac{1}{2}}\frac{T_{i+1}-T_i}{x_{i+1}-x_i}-{\lambda }_{i-\frac{1}{2}}\frac{T_i-T_{i-1}}{x_i-x_{i-1}}\right)+Q\left(x_i,t\right)$ . (23)

Для граничных ячеек запись уравнения зависит от используемого условия (задана температура, плотность потока тепла или граничное условие 3-го рода). Плотность является постоянной величиной, соответственно тепловой элемент с изменением температуры не меняет свой размер. Значения плотности, теплоемкости и удельного энерговыделения внутренних точек сетки определяются по величинам в центрах ячеек:

${С}_i=\frac{1}{2}\left({С}_{i-\frac{1}{2}}\left(x_{i-1}+x_i\right)+{С}_{i+\frac{1}{2}}\left(x_{i+1}+x_i\right)\right)$,

${\rho }_i=\frac{1}{2}\left({\rho }_{i-1/2}\left(x_{i-1}+x_i\right)+{\rho }_{i+1/2}\left(x_{i+1}+x_i\right)\right)$,

$Q\left(x_i,t\right)=\frac{1}{2}\left(Q\left(x_{i-1/2},t\right)\left(x_{i-1}+x_i\right)+Q\left(x_{i+1/2},t\right)\left(x_{i+1}+x_i\right)\right).$ (24)

Для аппроксимации по времени используется неявная четырехточечная разностная схема (температуры в правой части уравнения (23) берутся с верхнего временного слоя).

Система из I линейных разностных уравнений теплопроводности представляются в виде линеаризованных уравнений с коэффициентами A, B, C и D:

$B_1{T}_1^{*}+{С}_1{T}_2^{*}=D_1$,

$A_i{T}_{i-1}^{*}+B_i{T}_i^{*}+{С}_i{T}_{i+1}^{*}=D_i$,

$A_I{T}_{I-1}^{*}+B_I{T}_I^{*}=D_I$ (25)

и решается методом прогонки [6].

Замыкающие соотношения и модели оборудования

Теплообмен теплоносителя со стенами

Модель конвективной теплоотдачи от парогазовой среды к стенам предназначена для описания теплообмена между парогазовой средой в помещениях и поверхностью стен помещений и конструкционных элементов и основана на корреляциях из работы [7]. Результатами расчета данной модели являются: объемный источник/сток тепла в атмосферу помещения и плотность теплового потока q на границе теплового элемента. Коэффициенты теплоотдачи вычисляются по нескольким эмпирическим моделям. В качестве гидравлического диаметра D используется характерный размер поверхности или высота элемента. Для расчета режимов теплообмена среды с поверхностями используется пять чисел подобия:

$Pr=\frac{\mu \cdot c}{\lambda }$, $Re=\frac{G\cdot D}{\mu }$, $Gr=\frac{{\rho }^2\cdot g\cdot \beta \cdot \left|T_w-T_g\right|\cdot D^3}{ {\mu }^2}$,

$Ra=Gr\cdot Pr$, $Ri=\frac{Gr}{R{e}^2}$,

где $Pr$ – число Прандтля; $Re$ – число Рейнольдса; $Gr$ – число Грасгофа; $Ra$ – число Релея; $Ri$ – число Ричардсона. Для вычисления плотности массового расхода газа G используется среднеобъемная скорость течения газовой фазы, определяемая как максимум из скорости течения через гидравлические связи данного помещения $\frac{1}{2}\frac{D}{Vol}\underset{m}{}{A}_m\left|V_{gm}\right|$, скорость, определяемую по наличию источника $\frac{D}{Vol}{G}_{inj}{A}_{inj}{k}_{inj}$ и скоростью газа, вызванной передачей части импульса спринклерных капель. Для газа температурный коэффициент объемного расширения β принимается равным 1/T. Множитель $k_{inj}=1+\xi \frac{G_{inj,f}}{G_{inj,g}}$ учитывает передачу части импульса от жидкой к газовой фазе в случае источника двухфазной среды. Множитель $\xi$ ограничивает величину $k_{inj}$ сверху таким образом, чтобы скорость $\frac{D}{Vol}{G}_{inj}{A}_{inj}{k}_{inj}$ была меньше скорости истекающей воды.

Согласно [8] в качестве критерия для определения режима теплообмена (вынужденная или естественная конвекция) используется число Ричардсона, граничные значения которого приведены в табл. 1.

Таблица 1. Карта режимов теплообмена при совместном действии вынужденной и естественной конвекции для газов [8]

В модели для смешанной конвекции при расчете $N{u}_{mixed}$ используется подход, изложенный в [9] в области $Ri\cong 1$:

$N{u}_{mixed}= \phi \cdot N{u}_{forced}+\left(1-\phi \right)\cdot N{u}_{free}, где 0\le \phi =\frac{10-Ri}{10-1}\le 1$. (26)

В этом случае $N{u}_{mixed}$ получается на 2÷5% меньше, чем максимум из $N{u}_{forced} и N{u}_{free}$ и $N{u}_{forced} и N{u}_{free}$. Оценка показывает, что при характерном размере помещения 10 м, скорости газа 1 м/с и указанном выше допущении о линейной зависимости плотности газа от температуры получим, что с учетом $Ri=\frac{g\cdot \left|T_w-T_g\right|D}{ V_g^{2}T}$, при температурном напоре менее 1 К преобладает вынужденная конвекция, а для перепадов температур выше 50 К – естественная. Дополнительная погрешность модели за счет интерполяции в этом случае не превысит одного градуса.

Для помещений ЗО среднеобъемная скорость, используемая для расчета $\Delta P_{mom}$ в уравнении (21), рассчитывается из теплового баланса:

$Q_{g,j}=A_{w,j}{\eta }_{conv,j}\left(T_w-T_g\right)=2{\rho }_g{C}_{p,g}\left(T_w-T_g\right)V_{g,j}\frac{Vol}{D}.$ (27)

Она берется как максимум из индивидуальных скоростей газовой фазы вдоль вертикальных поверхностей теплообмена (стен):

$V_g=\underset{j}{\max }2\frac{{\eta }_{conv,j}}{{\rho }_g{C}_{p,g}}\frac{A_{w,j}D}{Vol}.$ (28)

Множитель 2 получен в предположении, что проходное сечение помещения $\frac{Vol}{D}$ поровну распределено между прямым и возвратным течением в пределах одного объема.

Теплообмен между газом и поверхностью за счет излучения основан на эмпирическом соотношении для поглощения излучения в водяном паре, полученном автором путем аппроксимации данных ВТИ [10, 11]:

${\eta }_{rad,j}=\sigma \cdot \text{min}\left(\frac{24}{T_g}\sqrt{{10}^{-5} P_{v}D }, \frac{1}{2}\right)\left({T_{wj}}^2+{T_g}^2\right)\left(T_{wj}+T_g\right)$. (29)

Суммарный коэффициент теплоотдачи .

Конвективный объемный источник тепла в газовую фазу помещения Q_wg получается суммированием потоков тепла по всем поверхностям j:

$Q_{wg}=\frac{1}{{\alpha }_g\cdot Vol}\underset{j}{}{A}_{w,j}\cdot {\eta }_{wg,j}\left({T_w}_j-T_g\right)$. (30)

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

Число Нуссельта при вынужденном течении вычисляется следующим образом:

$N{u}_{forced}=0.664\cdot R{e}^{0.5}\cdot P{r}^{0.33}$, ламинарное [12], (31)

$N{u}_{forced}=0.036\cdot R{e}^{0.8}\cdot P{r}^{0.33}$, турбулентное, $k_{inj}=1$ [12], (32)

$N{u}_{forced}=0.023\cdot R{e}^{0.8}\cdot P{r}^{0.33}$, турбулентное, $k_{inj}>36/23$ [9]. (33)

Число Нуссельта для вынужденного турбулентного течения интерполируется по $k_{inj}$, описанному выше. Отметим, что в выражениях (32) и (33) число Прандтля возводится в степень 0.33 вместо часто используемой 0.4. Для газа это приводит к отличию числа Нуссельта на 2.5%, но существенно упрощает интерполяцию между режимами при изменении состава газовой фазы.

Для естественной конвекции, учитывая диапазон чисел $Ra$ для газа для помещений ЗО от 10⁶ до 10¹⁵ (при наличии градиента температур между газом и стенкой 0.1 К или больше) [12]:

$N{u}_{free}=0.59\cdot R{a}^{0.25}$, ламинарная, (34)

$N{u}_{free}=0.13\cdot R{a}^{0.33}$, турбулентная. (35)

В [12] предлагается осуществлять переход от корреляции (34), описывающей область ламинарной конвекции и часть переходного режима к (35), описывающей турбулентную область, по граничному числу Ra = 10⁹. Для гладкого перехода граничное значение Ra принято равным 1.63·10⁸. Коэффициент конвективной теплоотдачи выбирается на основании числа Ричардсона согласно (26):

${\eta }_{conv}=\frac{\lambda }{D}\cdot F\left(N{u}_{forced}, N{u}_{free},\text{Ri}\right).$ (36)

Отметим, что при больших числах Ra (как правило преобладающем режиме теплообмена для помещений ЗО) коэффициент теплоотдачи не зависит от характерного размера в силу линейной зависимости числа Нуссельта от диаметра в (35).

Для горизонтальных поверхностей характерный размер D принимается равным отношению площади поверхности к ее периметру [13]. Коэффициент теплоотдачи при вынужденной конвекции рассчитывается аналогично изложенному выше алгоритму вычисления теплоотдачи к вертикальной стене. Вычисление коэффициента теплоотдачи к полу или потолку зависит от направления теплообмена. Для свободной конвекции используются соотношения, приведенные в табл. 2, они справедливы для Ra > 10⁵, при меньших числах Рэлея, как правило, преобладает вынужденная конвекция.

Таблица 2. Корреляции для расчета естественной конвекции для горизонтальных поверхностей

Модель конвективной теплоотдачи от воды в приямке к полу предназначена для описания теплообмена между приямками в помещениях и поверхностью подлежащего пола в предположении неподвижной жидкой фазы. Результатом расчета данной модели является объемный источник/сток тепла Q_wf в воду приямка и плотность теплового потока q на границе теплового элемента. Число Нуссельта для естественной конвекции:

${\text{Ra}}_f=\frac{{{\rho }_f}^{2}g{\beta }_f\left|T_w-T_f\right|{{\delta }_f}^3}{ {{\mu }_f}^2}P{r}_f, N{u}_{free}=0.54{{\text{Ra}}_f}^{1/4}$. (37)

Коэффициент конвективной теплоотдачи:

${\eta }_{wf}=\frac{{\lambda }_f}{{\delta }_f} \text{max}(N{u}_{free}, 0.55).$ (38)

В приведенных выше числах подобия β_f означает коэффициент теплового расширения жидкой фазы, толщина слоя воды δ_f ограничена снизу величиной 10^-3 м, перенос энергии за счет излучения не учитывается.

Объемный источник тепла в воду приямка:

$Q_{wf}=\frac{A_w{\eta }_{wf}}{{\alpha }_{f}Vol}\left(T_w-T_f\right).$ (39)

Интенсивность массообмена

Интенсивность объемной конденсации пара в помещениях ЗО считается пропорциональной его переохлаждению относительно линии насыщения:

${\Gamma }_i=-\frac{k \left(T_{sv}-T_g\right)C_{p,v} {\rho }_v}{h_v-h_{fs}}, k=\frac{\left|V_g\right|}{10}, 1\le k<10.$ (40)

Параметр k с размерностью 1/с зависит от интенсивности течения газа в помещении. Для неявной аппроксимации интенсивности объемной конденсации на новый временной слой используется выражение (14) с функцией $F=\left(P,h_g,X_1,\dots X_N\right).$

Интенсивность конденсации пара на поверхности, при наличии неконденсируемых газов, пропорциональна коэффициенту диффузии водяного пара D_v в атмосфере помещения. Поток массы конденсата в присутствии неконденсируемого газа рассчитывается следующим образом [9]:

${\Gamma }_w=D_v{\eta }_{conv}\frac{{\rho }_{vs}}{\lambda }\text{ ln}\frac{\left(1-P_{vs}/P\right)}{\left(1-P_v/P\right)}$. (41)

Плотность водяного пара на линии насыщения ${\rho }_{vs}$ определяется по парциальному давлению пара в объеме помещения $P_v$, давление насыщенного водяного пара $P_{vs}$ определяется по температуре стенки $T_w$. Коэффициент теплоотдачи для случая чистого пара или незначительного количества неконденсируемого газа ${\eta }_{conv}$ ограничивается сверху величиной ${\eta }_{conv,max}=\frac{{\lambda }_f}{\delta }$, что интерпретируется как диффузионное сопротивление пленки конденсата толщиной δ, которая зависит от ориентации поверхности. Критерии Нуссельта вычисляются в блоке конвективного теплообмена. В качестве гидравлического диаметра используется характерный размер помещения.

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

$D_i=\left(1-c_i\right)\cdot {\left(\underset{j=1}{\overset{i-1}{}}\frac{c_i}{D_{ij}}+\underset{j=i+1}{\overset{N}{}}\frac{c_i}{D_{ij}}\right)}^{-1},$ (42)

где $D_i$ – коэффициент диффузии компонента $i$ в смеси, м²∙с^-1; $D_{ij}$ – коэффициент бинарной диффузии газа $i$ в газе $j$, м²∙с^-1.

Коэффициент бинарной диффузии согласно [16] аппроксимирован выражением:

$D_{ij}=1.858\cdot {10}^{-12}\cdot \frac{{\text{T}}^{1.5}}{P{\Omega }_D{\sigma }_{ij}^2}\sqrt{\frac{1}{M_i}+\frac{1}{M_j}}$. (43)

Характеристическое расстояние и интеграл столкновений для потенциала Леннарда-Джонса вычисляются по формулам из работы [16].

В табл. 3 приведены значения характеристических энергий и расстояний для стандартных газов модуля CONT_TH. Как правило, для этих величин приводятся 3 значащие цифры, но в разных источниках обычно совпадают не более двух. Значения в таблице выбраны таким образом, чтобы с учетом приближений (42) и (43) для смесей H₂-H₂O-N₂-O₂ наилучшим образом соответствовать экспериментальным данным [17, 18, 19].

Таблица 3. Параметры потенциала Леннарда-Джонса

Модель спринклерной системы

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

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

Импульс, предаваемый в газовую фазу, позволяет получить среднеобъемную скорость движения газа в помещении и, по аналогии с выражением (28), использовать ее при расчете перемешивания газовой фазы за счет члена с ∆P_mom в уравнении (21). За базовую переменную взято время движения капли, изменяющееся от нуля до момента соприкосновения с полом последнего помещения по пути движения.

Система уравнений (44) связывает импульс капли, координату, массу и полную энтальпию капли (рис. 2):

Рис. 2. Движение спринклерной капли.

$\frac{d\left(M_d{V}_d\right)}{dt}=F_{grav}-F_{fr}, \frac{dz}{dt}=V_d\cos \frac{\phi }{2}$,

$\frac{d\left(M_d\right)}{dt}=G_{surf}, \frac{d\left(H_d\right)}{dt}=Q_{surf}+Q_{conv}.$ (44)

Капля, в силу малости размера, представляется в виде шара, взаимодействие между каплями не учитывается. Связь радиуса капли с ее объемом Ω, площадью поверхности S, поперечным сечением, массой и полной энтальпией выражается следующим образом:

 ${\Omega }_d=\frac{4}{3}\pi R_d^3,S_d=4\pi R_d^2,A_d=\pi R_d^2,M_d={\rho }_d{\Omega }_d,H_d={\rho }_d {\text{h}}_d.$(45)

Сила тяжести, действующая на каплю (в большинстве случаев угол φ достаточно мал и для упрощения дальнейших уравнений будем полагать, что действие силы тяжести всегда направлено вдоль линии движения капли):

$F_{grav}=\left({\rho }_d-{\rho }_g\right)\cdot g\cdot {\Omega }_d$ (46)

Сила трения о газ для случая витания (движение в турбулентной среде), может быть представлена в виде:

 $F_{fr}=C_d\cdot A_d\cdot \frac{1}{2} {\rho }_g{\left(V_d-V_{g,spr}\right)}^2,$(47)

где $V_{g,spr}$ – скорость газа в спринклерном конусе, м∙с^-1; $C_d$ – коэффициент сопротивления, значение которого зависит от числа Рейнольдса для капли $R{e}_d=2R_d\left|V_d-V_{g,spr}\right|$, и рассчитывается по зависимостям, предложенным в [20]. Число Шмидта для капли определяется через число Нуссельта $\frac{S{h}_d}{N{u}_d}={\left(\frac{{\lambda }_g}{D_v{C}_{p,g}}\right)}^{1/3}$, которое, согласно [21], можно выразить следующим образом:

 $N{u}_d=2.0+0.47R{e}_d^{0.5} .$(48)

Поток массы пара за счет конденсации или испарения вычисляется:

$G_{surf}=-D_v\left({\rho }_{vs,d}-{\rho }_v\right)S{h}_d\frac{2A_d}{R_d} ,$ (49)

где плотность пара на линии насыщения ${\rho }_{vs,d}$ определяется по температуре поверхности капли $T_d$.

Вычисление коэффициента диффузии водяного пара D_v аналогично приведенному выше случаю пристеночной конденсации. Соответствующие потоки тепла к капле за счет потока массы и конвективного теплообмена с газом соответственно определены как:

$Q_{surf}=G_{surf}{h}_v, Q_{conv}={\lambda }_g\left(T_g-T_d\right)N{u}_d\frac{2A_d}{R_d} .$ (50)

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

Скорость изменения радиуса капли во времени за счет массообмена рассчитывается из потока массы пара:

$\frac{d{R}_d}{dt}=-2D_v\cdot S{h}_d\cdot \frac{{\rho }_{vs,d}-{\rho }_v}{{\rho }_d{R}_d} .$ (51)

В [22] предлагается для средней температуры капли брать полусумму между температурой в центре капли и на ее поверхности. Но более точным приближением будет значение температуры в точке $\sqrt[3]{0.5}$ радиуса, отделяющей половину объема капли. В предположении малости внутренних конвективных потоков в капле из-за ее размера будем полагать, что перенос тепла в жидкости молекулярный и по балансу тепловых потоков за счет конвекции и массопереноса для средней температуры капли  получим:

${\lambda }_f\frac{\left(\overline{T_d}-T_d\right)}{1-\sqrt[3]{\mathrm{0,5}}}=-\frac{Q_{surf}+Q_{conv}}{2A_d}$. (52)

Для расчета средней скорости газа $V_{g,spr}$ внутри спринклерного конуса сделаем предположение, что весь импульс капли передается газовой фазе, тогда из закона сохранения импульса следует:

$G_{spr\mathrm{,0}}{V}_{d\mathrm{,0}}-G_{spr}{V}_d=\pi {\left(z\tan \frac{\phi }{2}\right)}^2\cdot {\rho }_g\cdot {V_{g,spr}}^2,$ (53)

где V_d,0 – начальная скорость капли, определяемая по расходу, диаметру сопла и коэффициенту истечения ξ:

 $V_{d\mathrm{,0}}=\frac{G_{spr}}{{\rho }_d}\cdot \frac{4}{\pi {D_{nozzle}}^2}\cdot \frac{1}{\xi } .$(54)

Приток массы газа в конус распыла спринклера согласно (53) может ограничивать поступление пара к поверхности капель и, соответственно, вводить ограничение на скорость его конденсации (49) через эффективное уменьшение коэффициента диффузии пара.

Система уравнений (44) интегрируется по времени, ограничением на шаг интегрирования являются: предельные изменения скорости капли, ее массы или температуры. Значение скорости газа в спринклерном конусе берется в явном виде согласно (53). Интегральные значения потоков массы и энергии в помещениях определяются на основании диапазона высот и количества капель, пролетающих эти помещения в единицу времени.

Если необходимо учитывать взаимодействие спринклера со стенкой, то в каждом помещении выделяется одна тепловая структура, на которую попадает заданная пользователем доля от всех капель θ. При соударении со стенкой в помещении j (предполагается, что это происходит в самой нижней точке помещения) капли принимают температуру стенки и мгновенно перемещаются в приямок. Соответствующий тепловой поток в стенку составит:

$Q_{wj,spr}={\theta }_j{G}_{spr}\left(z\right)\left(h_d-h_{f,wj}\right).$ (55)

Отметим, что в этом случае поток спринклерных капель в (53) зависит не только от изменения их массы, но и от их количества. Среднеобъемная скорость газа за счет передачи импульса от спрея составит:

$V_{g,j}=2 \pi {\left(z_j\tan \frac{\phi }{2}\right)}^2\cdot {\rho }_{gj}\cdot V_{g,spr,j}{\left(\frac{D}{Vol}\right)}_j.$ (56)

Моделирование других систем

Для моделирования взаимодействия водорода с кислородом используются следующие балансные соотношения:

$\frac{d{M}_{H2}}{dt}=-G_{rec} , \frac{d{M}_{O2}}{dt}=-8G_{rec} , \frac{d{M}_v}{dt}=9G_{rec}$. (57)

Скорость рекомбинации G_rec задается пользователем и может зависеть от произвольного числа параметров. Это позволяет использовать модель и для расчета горения водорода в адиабатическом изохорном приближении. Тепловыделение в газовую фазу принято равным 120 МДж на 1 кг водорода.

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

${\dot{\rho }}_m{V}_{gm}^{\text{*}}{A}_m=G_m.$ (58)

Для ограничения расхода парогазовой смеси при достижении потоком скорости звука в уравнение (21) вводится эффективное сопротивление $k_{crit}>k_w+k_{local}$. Массовая скорость $J_{g,crit}$ для смеси газов вычисляется с использованием эквивалентного показателя адиабаты γ с использованием донорных концентраций компонентов газовой фазы:

$J_{g,crit}=\psi {\left(\frac{2}{\gamma +1}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\gamma -1$}\right.}\sqrt{\frac{2\gamma }{\gamma +1}\dot{P}{\dot{\rho }}_g}, \gamma =\frac{\frac{{\dot{c}}_v{\gamma }_v}{{\gamma }_v-1}+{{\displaystyle \sum }}_{n=1}^N\frac{{\dot{c}}_n{\gamma }_n}{{\gamma }_n-1}}{\frac{{\dot{c}}_v}{{\gamma }_v-1}+{{\displaystyle \sum }}_{n=1}^N\frac{{\dot{c}}_n}{{\gamma }_n-1}} .$ (59)

Коэффициент истечения ψ для отверстий принимается равным 0.82 (острая кромка), а для более длинных каналов – учитывает трение о стенку. Критический расход наступает при условии перепада давлений между помещениями более $\frac{{J_{g,crit}}^2}{2{\dot{\rho }}_g}$, уравнение сохранения импульса для связи m записывается следующим образом:

$\tau \frac{{\rho }_j\Delta z_k+{\rho }_k\Delta z_j}{\Delta z_k+\Delta z_j}\left(V_{gm}^{\text{*}}-V_{gm}\right)+\tau \frac{P_k^{*}-P_j^{*}}{\Delta z_k+\Delta z_j}+\tau \frac{\Delta P_{mom}^{}}{\Delta z_k+\Delta z_j}+\tau \frac{\Delta P_{den}^{}}{\Delta z_k+\Delta z_j}=$

$=\tau \left(k_{crit}\right)V_g^{*}-\tau g\frac{\left({\rho }_{gk}\left(H_k-H_m\right)+{\rho }_{gj}\left(H_m-H_j\right)\right)}{\Delta z_k+\Delta z_j}$ ,

$k_{crit}=\frac{\left|P_k-P_j\right|{\dot{\rho }}_g}{J_{g,crit,m} }\frac{2}{ \Delta z_k+\Delta z_j} .$(60)

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

$\frac{d{M}_f}{dt}=G_{f,src} , \frac{d{h}_f}{dt}=h_f+\frac{G_{f,src}{h}_{f,src}}{M_f} .$ (61)

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

Используемые уравнения состояния воды и водяного пара вместе с расчетом производных по определяющим величинам приведены в таблицах 4 и 5.

Таблица 4. Уравнения состояния водяного пара

Таблица 5. Уравнения состояния воды

Для расчета поверхностного натяжения воды, вязкости и теплопроводности используются полиномиальные зависимости, отклонения от данных [24] в стабильной области не превышают 1÷2%.

В модели неконденсируемые газы удовлетворяют уравнению состояния идеального газа:

$\rho \left(P,h\right)=\frac{\gamma }{\gamma -1}\cdot \frac{P}{h} , h\left(T\right)=C_{p}T, P\left(\rho ,T\right)=C_p\frac{\gamma -1}{\gamma }\rho T$. (62)

Параметры уравнений состояния для стандартных газов модуля приведены в табл. 6. Отметим, что коэффициент адиабатического сжатия не является теоретической величиной. Он получен автором из наилучшего соответствия расчетной плотности (62) экспериментальным значениям в области давлений 10⁵÷10⁶ Па и температур 250÷1250 К в предположении постоянной теплоемкости. Теплоемкость соответствует среднему значению в интервале 300÷600 К. Отметим, что из-за неидеальности газов показатель адиабаты может быть как выше, так и ниже теоретического значения.

Таблица 6. Параметры уравнений состояния идеального газа

Для расчета вязкости идеальных газов используются полиномиальные зависимости второго порядка по температуре, для расчета тепловодности – степенная зависимость с показателем от 0.67 до 0.88 в зависимости от газа.

Для определения ρ_g, T_g, ρ_v, P_v и их производных, через P, h_g, X₁...X_N используется система:

${\rho }_g={\rho }_v+{\rho }_g{\sum }_{n=1}^N{X}_n;h_g{\rho }_g={\rho }_v{h}_v+{\rho }_g{\sum }_{n=1}^N{X}_n{h}_n; P=P_v+{\sum }_{n=1}^N{P}_n$

${\rho }_v={\rho }_v\left(P_v,h_v\right); h_v=h_v\left(P_v,T_g\right); P_n=h_n{\rho }_g{X}_n\frac{{\gamma }_n}{{\gamma }_n-1}; h_n=C_{Pn}{T}_g$. (63)

Здесь C_Pn теплоемкость неконденсируемого газа n. Система линеаризуется с использованием производных плотности водяного пара по давлению и энтальпии, энтальпии водяного пара по давлению и температуре и аппроксимации величин на новый временной слой:

${\rho }_g^{*}{h}_n^{*}={\rho }_g{h}_n^{*}+{\rho }_g^{*}{h}_n-{\rho }_g{h}_n ,$  ${\rho }_g^{*}{T}_g^{*}={\rho }_g{T}_g^{*}+{\rho }_g^{*}{T}_g-{\rho }_g{T}_g .$ (64)

Теплоемкость смеси пара с неконденсируемыми газами определяется через массовые доли компонентов и их теплоемкости:

$C_p=C_{p,v}\left(1-{\sum }_{n=1}^N{X}_n\right)+ {\sum }_{n=1}^N{C}_{p,n}{X}_n .$ (65)

Динамическая вязкость смеси вычисляется по формуле, которая является аппроксимацией точного решения на основе строгой кинетической теории Чэпмэна-Энскога для случая низких давлений [16], параметры Ф_nm вычисляются на основе аппроксимации Вильке:

$\mu =\underset{n=0}{\overset{N}{}}\frac{c_n{\mu }_n}{{\displaystyle {\sum }_{m=0}^N{c}_n{\text{Ô}}_{nm}}}, {\text{Ô}}_{nm}=\frac{{\left[1+\sqrt{\raisebox{1ex}{${\mu }_n$}\!\left/ \!\raisebox{-1ex}{${\mu }_m$}\right.}{\left(\raisebox{1ex}{$M_m$}\!\left/ \!\raisebox{-1ex}{$M_n$}\right.\right)}^{0.25}\right]}^2}{\sqrt{8\left(1+\raisebox{1ex}{$M_n$}\!\left/ \!\raisebox{-1ex}{$M_m$}\right.\right)}} .$ (66)

где n=0 и m=0 для упрощения записи соответствует водяному пару, Mn – молярная масса n-ого компонента, кг/моль. Отклонение вязкости смеси, рассчитанной по этой формуле, от измерений составляет меньше 2% даже для смесей, содержащих полярные газы.

Теплопроводность смеси при низких давлениях вычисляется по формуле Васильевой [16], A_ij вычисляется на основе модификации Линдсея и Бромли:

$\lambda =\underset{n=0}{\overset{N}{}}\frac{c_n{\lambda }_n}{{\displaystyle {\sum }_{m=0}^N{c}_n{A}_{nm}}}, {\text{A}}_{nm}=\left(1+\sqrt{\frac{{\mu }_n}{{\mu }_m}{\left(\frac{M_m}{M_n}\right)}^{3/4}\frac{T+S_n}{T+S_m}}\right)\frac{T+S_{nm}}{4\left(T+S_n\right)},$ (67)

где $S_n$ – постоянная Сюзерленда для n-ого компонента, которая может быть оценена по эмпирической формуле $S_n=1.5\cdot T_{b,n}$, где $T_{b,n}$ – температура кипения n-ого компонента при нормальных условиях, К; $S_{nm}=C_S\sqrt{S_n\cdot S_m}$ – постоянная взаимодействия Сюзерленда, $C_S=1$ для случая, когда оба компонента являются неполярными газами и $C_S=0.733$ если один из двух компонентов является полярным газом [25]. Данная формула обеспечивает среднее отклонение от измерений ~3% для неполярных смесей и ~5% для случая полярно-неполярных смесей.

Решение системы уравнений и расчетный шаг

Система уравнений (7)–(12) с учетом пространственной аппроксимации (19)–(21) является простой и легко разрешимой системой разностных уравнений. Но она является неконсервативной с точки зрения сохранения массы в силу того, что с одной стороны плотность газовой фазы определяется по уравнению (15), с другой стороны является функцией базовых переменных – давления, энтальпии и состава газовой фазы:

${\rho }_g^{*}={\rho }_g\left(P_g^{*},h_g^{*},X_1^{*}\dots X_N^{*}\right)$,

${\rho }_g^{*} \ne {\rho }_g+\frac{\partial {\rho }_g}{\partial P}\left(P^{*}-P\right)+\frac{\partial {\rho }_g}{\partial h_g}\left(h_g^{*}-h_g\right)+\underset{n=1}{\overset{N}{}}\frac{\partial {\rho }_g}{\partial X_n}\left(X_n^{*}-X_n\right)$. (68)

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

$\begin{array}{c}{\Delta }_g={\alpha }_g^{*}\left({\rho }_g-{\rho }_g\left(P_g^{*},h_g^{*},X_1^{*}\dots X_N^{*}\right)\right)-\\ -{\alpha }_g\left(\frac{\partial {\rho }_g}{\partial P}\left(P^{*}-P\right)+\frac{\partial {\rho }_g}{\partial h_g}\left(h_g^{*}-h_g\right)+\underset{n=1}{\overset{N}{}}\frac{\partial {\rho }_g}{\partial X_n}\left(X_n^{*}-X_n\right)\right)\end{array}$. (69)

Эта погрешность затем используется а качестве дополнительного источника массы S(∆g) в уравнении (7). Относительная погрешность массы газовой фазы:

${\delta }_g=\underset{k}{\max }\raisebox{1ex}{$2\left|{}_{\Delta gk}\right|$}\!\left/ \!\raisebox{-1ex}{$\left({\rho }_{gk}+{\rho }_{gk}^{*}\right)$}\right.$ (70)

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

Для установления границ возможного шага интегрирования используется относительная погрешность массы газовой фазы δ = max(δ_g, δ_f, δ₁...δ_N), шаг уменьшается если δ > δ_max и может быть увеличен если δ < δ_min.

На рис. 3 показаны основные этапы расчетного шага интегрирования модуля CONT_TH. Расчет датчиков системы контроля и управления обеспечивает не только данные для граничных условий, но и необходимую пользователю расчетную информацию. Взаимодействие с другими модулями кода СОКРАТ осуществляется в явном виде по состоянию на начало шага. Начальное значение шага интегрирования получается из результатов расчета предыдущего шага с учетом задаваемых пользователем границ. Оно может быть скорректировано в случае выхода базовых переменных за рамки допустимых значений или резкого изменения их значений, нефизичных или резко изменившихся значений вторичных параметров (температуры и плотности фаз).

Рис. 3. Расчетный шаг модуля CONT_TH.

Заключение

Модель реализована в составе ПрЭВМ СОКРАТ-В1/В2 в виде модуля CONT_TH и аттестована для расчета противодавления в ЗО при моделировании протекания тяжелой аварии на РУ ВВЭР в связанной постановке, а также для расчета распределения водородосодержащих смесей по помещениям ЗО для оценок предельного давления в условиях их адиабатического сгорания. Совместно с модулем CONT_FP кода СОКРАТ-В3 обеспечивает расчет динамики поведения аэрозолей в ЗО.

К преимуществам модуля CONT_TH по сравнению с другими расчетными кодами с сосредоточенными параметрами можно отнести:

•  единый алгоритм расчета свойств воды, водяного пара и неконденсируемых газов с моделью контурной гидравлики ПрЭВМ СОКРАТ;
•  общий язык описания входных наборов;
•  преимущественное использование автомодельных или слабо зависящих от характерного размера замыкающих соотношений, что упрощает построение расчетных схем;
•  возможность моделировать различное оборудование (в том числе пассивные системы безопасности);
•  низкие затраты процессорного времени.

Количественные возможности модели, в том числе при работе систем безопасности ЗО, используемых для смягчения последствий ТА, будут обсуждены в следующей статье автора.

Модель при незначительной адаптации может быть применена для расчета параметров ЗО РУ с водяным теплоносителем в составе других интегральных кодов.

Список обозначений:

Sobre autores

D. Tomashchik

Autor responsável pela correspondência
Email: tdyu@ibrae.ac.ru
Rússia, Moscow

Bibliografia

Arquivos suplementares

Arquivos suplementares

Ação

1. JATS XML

Baixar

2. Fig. 1. Example of a nodalization scheme for the CONT_TH module.

Baixar (28KB)

Metadados

3. Fig. 2. Movement of a sprinkler drop.

Baixar (6KB)

Metadados

4. Fig. 3. Calculation step of the CONT_TH module.

Baixar (34KB)

Metadados