Теплофизика высоких температур, 2022, T. 60, № 6, стр. 906-915

ВВЕДЕНИЕ

В процессе развития тяжелой аварии с расплавлением активной зоны водо-водяного реактора (ВВЭР) расплав активной зоны, содержащий тепловыделяющие продукты деления ядерного топлива, приходит в контакт с корпусом реактора и при отсутствии внешнего охлаждения проплавляет корпус, выходя в бетонную шахту реактора. Удержание расплава внутри корпуса при аварии реактора мощностью 1 ГВт (эл.) и более до сравнительно недавнего времени полагалось нереалистичным, и, хотя концепция удержания уже принята в нескольких новых зарубежных проектах типа AP1000, ее возможность для ныне действующих реакторов (далее речь пойдет о реакторах ВВЭР-1000, не снабженных устройством локализации расплава) пока не обоснована полностью. Проведенные в разное время экспериментальные и численные исследования, аналитические и методами вычислительной гидродинамики [1–6], показывают, что тепловые нагрузки на корпус реактора могут находиться в пределах, допускаемых возможностями внешнего охлаждения и требованиями механической прочности. Но для существующих блоков с ВВЭР-1000 реализация концепции внутрикорпусного удержания расплава требует ряда исследований и конструктивных доработок. В связи с этим актуален вопрос о долговременных последствиях выхода расплава в бетонную шахту и эволюции состояния всей АЭС при тяжелой аварии, что определяет меры по ее локализации.

Проблема взаимодействия расплав–бетон (ВРБ) изучается с 1970-х годов [7]. В мире разработан ряд расчетных средств моделирования ВРБ, но применяемые подходы и методы расчета пока дают ответ на поставленный вопрос со слишком большой погрешностью. Это показал сравнительный анализ, проведенный группой WGAMA [8] на примере двенадцати расчетных кодов, предлагаемых для задач ВРБ. К тому же, вследствие существенных отличий строения помещений и расположения реактора ВВЭР-1000 от зарубежных аналогов, постановка задачи тут имеет ряд особенностей. В частности, расчетное средство должно учитывать возможность распространения расплава за пределы бетонной шахты в соседние помещения, вследствие вероятного проплавления стальной двери, закрывающей вход в шахту [9 ] . Для расчета же долговременной эволюции требования к точности численных моделей должны быть существенно повышены. Не детализируя неточности существующих моделей ВРБ, отметим три основные: недостаточно полные и корректные модели термохимии и теплофизики разложения бетона и его взаимодействия с расплавом (как правило, это – очередность химических реакций по назначенным приоритетам и не вполне корректный учет состава бетона, исходящий из остаточных компонентов после прокаливания); весьма грубые модели продвижения границы эрозии бетона, основанные на геометрических построениях и одномерных аналитических аппроксимациях поля температуры в приграничном слое бетона [8]; неучтенные возможности растекания расплава по большой площади.

В настоящей работе рассмотрены уточняющие подходы к решению первых двух проблем: описаны термодинамическая модель химических реакций компонентов бетона с расплавом, учитывающая также продукты деления ядерного топлива, модель разложения бетона и основанные на этом универсальные процедуры расчета теплоемкости и плотности бетона, не зависящие от его типа. Совокупное применение таких процедур для сквозного расчета ВРБ (включая описанное в [8]) авторам неизвестно. Описана одномерная версия модели продвижения границы эрозии бетона, применимая для случаев плоского бассейна расплава с малым аспектным отношением.

О ЗАДАЧЕ МОДЕЛИРОВАНИЯ ЭРОЗИИ БЕТОНА В КОНТАКТЕ С РАСПЛАВОМ АКТИВНОЙ ЗОНЫ РЕАКТОРА

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

Зона прогрева бетона расплавом, определяемая теплопроводностью, узкая на начальной стадии ВРБ, когда мощность велика. Но прогрев может охватить всю оставшуюся толщину перекрытия ~1 м на поздней стадии процесса – через несколько дней или недель. Так как прочность бетона снижается практически до нуля при прогреве до 400°С [8], температура бетона должна определяться с достаточной точностью. Для целей ВРБ численные модели типа конечно-элементных требуют оптимизации, поскольку задача нестационарна, и необходимо разрешение ~1–4 мм, дающее слишком большое количество ячеек, даже в 2D-приближении. Загрубление разрешения до 1–10 см и более дает неконтролируемые искажения. В расчетных кодах, используемых в настоящее время (CORCON [10], MEDICIS [11], CorQuench [12], WECHSL [7]) зона прогрева обычно не моделируется явно, и модель продвижения 2D-границы эрозии основана на балансных соотношениях типа условия Стефана [13]. Новое положение граничной линии при неравномерном распределении граничного потока тепла находится геометрическими построениями, корректность которых оценивается на простых тестах.

Последовательный расчет продвижения границы плавления должен исходить из решения двумерной задачи теплопроводности. 3D-постановки могут быть сведены к 2D. В экономичной сеточной модели эрозии используется то, что зона прогрева сначала довольно узкая. Проблема сеточного разрешения устраняется путем ограничения расчетной сетки на приграничную зону прогрева небольшой ширины и регулярной перестройкой сетки, аналогично процедурам, давно применяемым в вычислительной гидродинамике [14]. Если разрешение сетки неизменно, то в 2D-геометрии количество узлов сетки сначала небольшое и увеличивается пропорционально размеру каверны. Это позволяет ускорить расчет глубокого проплавления при ВРБ на порядок по сравнению с фиксированной сеткой требуемого размера. Процедура разработана в 1D- и 2D-исполнениях. В данной работе использована 1D-модель, достаточно реалистичная с точки зрения большинства, хотя и не всех, вероятных конфигураций бассейна расплава в шахте ВВЭР-1000.

МОДЕЛИРОВАНИЕ ХИМИЧЕСКИХ РЕАКЦИЙ КОМПОНЕНТОВ РАСПЛАВА И БЕТОНА

Моделируемый состав расплава. Расплав на днище корпуса реактора представляет собой многокомпонентный раствор, в который входят компоненты топлива, управляющих элементов активной зоны реактора, стали конструкций, радиоактивные продукты деления (ПД). В процессе ВРБ добавляются компоненты бетона. Элементный состав системы включает U и Zr, входящие в ядерное топливо (UO₂) и оболочки, Fe, Ni, Cr – основные компоненты нержавеющей стали, Al, Ca, Na, K, Si, Mg – основные элементы, продукты декомпозиции бетона, и продукты деления, представленные в модели 18 классами элементов [15]: Cs, I, Mo, Ru, Ba, Sr и др. Сюда входят также вода Н₂О и оксиды CO и CO₂, продуцирующие водород, кислород и углерод.

В модели рассматривается в общем случае трехфазная система, состоящая из двух конденсированных (оксидная и металлическая) и газовой фаз. Металлическая фаза расплава включает компоненты стали и оболочек: Fe, Ni, Cr, Zr, металлический уран, который может появиться в результате восстановления из UO₂, а также металлы Al, Ca, Na, K и Si, которые являются продуктами восстановления из соответствующих оксидных компонентов бетона. Это возможно, когда металлическая фаза содержит значительное количество Zr, напрямую взаимодействующего с бетоном посредством реакций типа Zr + SiO₂ = ZrO₂ + Si. Кроме того, к металлической фазе следует отнести металлические ПД – Mo, Ru, Pd, Rh, Te.

В состав оксидной фазы, помимо топлива UO₂ и окисленной части циркония ZrO₂, входят оксиды компонентов стали (FeO, Fe₂O₃, Cr₂O₃, NiO) и продуктов деления, а также компоненты декомпозиции бетона Al₂O₃, CaO, Na₂O, K₂O, SiO₂, MgO, MnO. Кроме того, в случаях, когда в процессе окисления циркония кислородный потенциал оксидной фазы невысок, в оксидной фазе могут присутствовать в небольших количествах металлы Fe, Cr, Ni и другие.

Газовая фаза включает неконденсируемые H₂O, H₂, O₂, CO₂, входящие в состав атмосферы над расплавом (азот не учитывается) или газовых пузырей, возникающих при взаимодействии расплава с бетоном, и около ста газообразных компонентов ПД, расплава и бетона.

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

• кинетике химических реакций ${{\tau }^{{{\text{kin}}}}}$,

• выравниванию состава за счет массопереноса (конвекции, диффузии) ${{\tau }^{{{\text{mass}}}}}$,

• радиоактивным превращениям ${{\tau }^{{{\text{rad}}}}}$.

Газовыделение при взаимодействии расплава с бетоном эффективно выравнивает состав, и характерное время массопереноса ${{\tau }^{{{\text{mass}}}}}$, оцененное как стократное время всплытия пузыря, не превышает ~100 с. При температурах расплава T ~ ~ 1700−2200 К скорости химических реакций достаточно велики, чтобы выполнялось условие ${{\tau }^{{{\text{kin}}}}} \ll $ ${{\tau }^{{{\text{mass}}}}}$, и на интервале времени ~${{\tau }^{{{\text{mass}}}}}$ система находилась в состоянии локального термохимического равновесия. Радиоактивный распад ПД, являющихся компонентами расплава, приводит, кроме тепловыделения, к изменению элементного и, следовательно, молекулярного состава системы. Как показано в [15], времена жизни большинства нуклидов, ответственных за тепловыделение, отвечают неравенству ${{\tau }^{{{\text{mass}}}}} \ll {{\tau }^{{{\text{rad}}}}}$, и на временах ~${{\tau }^{{{\text{mass}}}}}$ в системе устанавливается термодинамически равновесное состояние, которое на масштабе ${{\tau }^{{{\text{rad}}}}}$ меняется квазистатически, по мере переноса компонентов расплава и изменений в изотопном и элементном составах, связанных с радиоактивным распадом.

В описании равновесия многокомпонентной, многофазной системы далее используются следующие обозначения: ${{n}_{e}}$ – полное число элементов (одноатомные компоненты) в системе; ${{n}_{p}}$ – число фаз в системе; ${{n}_{s}}$ – полное число компонентов в системе; верхний индекс $p$ задает номер фазы; нижние индексы нумеруют компоненты; $N_{i}^{{(p)}}$ – число молей компонента $i$ в фазе p; ${{N}^{{(p)}}} = \sum\nolimits_{i \in {{I}_{p}}} {N_{i}^{{(p)}}} $ – суммарное количество молей в фазе $p;$ $x_{i}^{{(p)}} = {{N_{i}^{{(p)}}} \mathord{\left/ {\vphantom {{N_{i}^{{(p)}}} {{{N}^{{(p)}}}}}} \right. \kern-0em} {{{N}^{{(p)}}}}}$ – мольная доля $i$-го компонента в $p$-й фазе; ${{b}_{{ji}}}$ – стехиометрические коэффициенты, задающие число атомов i в молекуле j.

Для расчета параметров равновесного состояния {N^(eq)} системы находится абсолютный минимум свободной энергии Гиббса G_tot, взятой как функция количества молей компонентов N_i, при заданных температуре T, давлении P_tot и элементном составе системы, который связан с компонентным составом соотношениями

Очевидное дополнительное условие состоит в положительности

Минимизация G_tot сводится к решению системы нелинейных алгебраических уравнений

где A_i – активность i-го вещества. Для компонентов конденсированной фазы активность пропорциональна мольной доле вещества где ${{\gamma }_{i}}$ – коэффициент активности. В случае газовой фазы активность компонента A_i совпадает с его парциальным давлением, выраженным в атмосферах: где p₀ = 1.01325 × 10⁵ Па, P_tot – полное давление в газовой фазе. В (3) использовано уравнение состояния идеального газа. Коэффициенты активности в уравнении (2) в общем случае являются функциями состава расплава и температуры. В работе используется упрощенная модель – коэффициенты ${{\gamma }_{i}}$ = 1 для всех компонентов, за исключением основных ПД: Sr, Ba, Ce, La, Mo, Ru, которые описываются моделью неидеального раствора [15].

Константы равновесия ${{K}_{i}}$ в уравнениях (1) зависят только от температуры и выражаются через изменение стандартной энергии Гиббса реакции образования вещества i:

где $\mu _{i}^{0}(T)$ – химический потенциал, или энергия Гиббса одного моля чистого вещества j в стандартном состоянии, $R$ – универсальная газовая постоянная. Данные о химических потенциалах чистых компонентов взяты из термодинамической базы ИВТАНТЕРМО [16, 17] и справочных изданий [18, 19].

РАСЧЕТ ПОЛНОЙ ТЕПЛОЕМКОСТИ БЕТОНА

Модель фактического состава бетона. Бетон – это неоднородный, многокомпонентный и многофазный материал, состав которого существенно изменяется при нагревании. Состав бетона обычно дается в виде набора компонентов – продуктов его теплового разложения, называемых далее остаточными компонентами. Но исходный состав бетона включает соединения, отличные от остаточных компонентов. Основные физико-химические процессы, сопровождающие нагрев бетона и превращения соединений, следующие [8, 20]:

– испарение содержащейся в порах, связанной и адсорбированной воды при температуре около 400 К. Cюда же можно отнести дегидратацию гидроалюминатов, равно как и распад и перекристаллизацию гидросульфоалюминатов кальция (3CaSO₄⋅3СаО⋅Al₂O₃⋅3Н₂О) в диапазоне температур 473–673 К;

– разложение гидроксидов кальция и магния

при температурах около 760 и 520 К соответственно с образованием водяного пара;

– разложение карбонатов кальция и магния по реакции

при температурах ~1100 и ~560 К соответственно с образованием углекислого газа;

– оставшаяся оксидная матрица в контакте с расплавом плавится при температуре в пределах 1500–2400 К [8] для бетонов, используемых при строительстве АЭС. В контакте с расплавом плавление этих оксидов происходит растворением в оксидной фазе расплава.

Теплоемкость бетона связана с взаимодействием и поглощением тепла на трех уровнях структуры: 1) внутри остаточных компонентов; 2) внутри молекул соединений, при разложении которых образовались эти остаточные компоненты; 3) взаимодействие этих молекул, объединенных в агломераты не вполне определенной структуры, внутри которых находятся еще молекулы воды. Расчет теплоемкости на основе данной структуры практически невозможен. Расчет только по остаточному составу заведомо некорректен. Приближенный метод, взятый в настоящей работе, базируется на следующих предположениях:

• бетон полагается состоящим из соединений уровня 2 – например, для Ca это молекулы CaO, Ca(OH)₂, CaCO₃;

• взаимодействие на уровне 3 не учитывается, вода  рассматривается как отдельное соединение, основные энергозатраты связаны с ее испарением;

• теплоемкость соединений уровня 2 включает вклады остаточных компонентов, рассчитываемые по справочным данным в предположении аддитивности вкладов.

Это позволяeт построить модель и ввести ее константы, практически не обращаясь к эмпирическим данным по теплоемкостям бетонов. Теплоты реакций разложения известны. При известном полном количестве остаточной воды неопределенным параметром остается количество химически несвязанной (“свободной”) воды, которое должно быть задано извне.

Пусть начальный состав бетона задан массовыми долями компонентов ${{w}_{j}}$, в частности, заданы доли оксидов кальция и магния, “свободной” и “связанной” (H₂O^(bond)) воды и газообразного диоксида углерода: $w_{{{\text{CaO}}}}^{{\left( 0 \right)}}$, $w_{{{\text{MgO}}}}^{{\left( 0 \right)}}$, $w_{{{{{\text{H}}}_{2}}{\text{O}}}}^{{\left( {{\text{free}}} \right)}}$, $w_{{{{{\text{H}}}_{2}}{\text{O}}}}^{{\left( {{\text{bond}}} \right)}}$, $w_{{{\text{C}}{{{\text{O}}}_{2}}}}^{{\left( 0 \right)}}$. В исходном состоянии связанная вода содержится в бетоне в виде гидроксидов, а СО₂ есть продукт разложения карбонатов кальция и магния. Поэтому полагается, что при температуре ${{T}_{0}}$ = 298.15 К оксиды кальция и магния входят в состав бетона в виде шести соединений: (Ca, Mg)(OH)₂, (Ca, Mg)CO₃ и (Ca, Mg)O. Весовые доли первых двух можно определить по начальным весовым долям связанной воды и диоксида углерода, полагая, что CO₂ и H₂O^(bond) входят как соединения с кальцием и магнием:

где символ Me обозначает Ca или Mg, параметр ${{\alpha }_{{{\text{Ca}}}}} = {{w_{{{\text{CaO}}}}^{{(0)}}} \mathord{\left/ {\vphantom {{w_{{{\text{CaO}}}}^{{(0)}}} {\left( {w_{{{\text{CaO}}}}^{{(0)}} + w_{{{\text{MgO}}}}^{{(0)}}} \right)}}} \right. \kern-0em} {\left( {w_{{{\text{CaO}}}}^{{(0)}} + w_{{{\text{MgO}}}}^{{(0)}}} \right)}} = 1 - {{\alpha }_{{{\text{Mg}}}}}$, ${{M}_{j}}$ – молярные массы компонентов (кг/моль). Весовые доли CaO и MgO в модельном составе бетона ${{w}_{{{\text{CaO}}}}}$ и ${{w}_{{{\text{MgO}}}}}$ определяются из условия сохранения массы Ca и Mg:

Интервалы температур и состав бетона при разложении. При нагреве бетона сначала испаряется свободная вода, потом начинают разлагаться гидраты кальция и магния с выходом воды. Далее разлагаются карбонаты кальция/магния, выходит углекислота. Затем начинается плавление (растворение) остаточных компонентов, в расплав поступают оксиды, включая окислы железа и кремния (SiO₂). В модели разложения бетона рассматриваются пять температурных интервалов, за которые состав претерпевает изменения. Это – интервал температур, в котором испаряется “свободная” вода, и четыре интервала разложения гидроксидов и карбонатов магния и кальция, параметры которых приведены в табл. 1. Теплоты разложения рассчитаны по данным ИВТАНТЕРМО [16, 17].

Полагается, что свободная вода испаряется в некотором интервале температур. Масса свободной воды $m_{{{{{\text{H}}}_{{\text{2}}}}{\text{O}}}}^{{{\text{(free)}}}}$ зависит от температуры по линейному закону

где по умолчанию ${{T}_{0}}$ = 298.15 К, ${{T}_{1}}$ = 400, δT₁ = 10 и введена функция θ(x) – функция Хэвисайда. Здесь и дальше используется обозначение $T_{i}^{ \pm } = {{T}_{i}} \pm \delta {{T}_{i}}$.

В остальных температурных интервалах $T_{j}^{ \pm }$, j = = 2–5, перечисленных в табл. 1 (нумерация соответствует компонентам Mg(OH)₂, MgCO₃, Ca(OH)₂, CaCO₃), происходит разложение гидроксидов и карбонатов по реакциям (4), (5). При этом масса соединения ${{m}_{j}}$ линейно убывает, а масса соответствующего оксида m_k(j) (k(2) = k(3) = MgO, k(4) = k(5) = = CaO) линейно растет с ростом температуры:

Соотношения (6)–(8) позволяют рассчитать температурную зависимость состава бетона, выраженного в форме массовых долей компонентов ${{w}_{j}}\left( T \right)$, и параметр $\mu \left( T \right) = \frac{{{{m}_{{{\text{conc}}}}}\left( T \right)}}{{{{m}_{{{\text{conc}}}}}\left( {{{T}_{0}}} \right)}}$, описывающий изменение массы бетона с температурой.

Плотность и полная теплоемкость бетона. Теплоемкость бетона вычисляется по аддитивной модели (см., например, [10]):

где суммирование выполняется по всем компонентам бетона с учетом описанной выше модификации исходного состава с весовыми долями ${{w}_{j}}\left( T \right)$. Теплоемкости компонентов ${{c}_{j}}\left( T \right)$ рассчитываются по данным, приведенным в ИВТАНТЕРМО [16, 17].

Температурная зависимость плотности бетона ${{\rho }_{{{\text{conc}}}}}\left( T \right)$ описывается уравнением

Данная формула дает изменение начальной плотности ${{\rho }_{{{\text{conc}}}}}\left( {{{T}_{0}}} \right)$, связанное с уходом газов, при сохранении объема, занимаемого бетоном. Энтальпия бетона определяется формулой

где первое слагаемое соответствует аддитивной модели; ${{m}_{j}}\left( T \right)$ – массы компонентов бетона, зависимости которых от температуры даются уравнениями (6)–(9); ${{h}_{j}}\left( T \right)$ – удельные энтальпии компонентов, вычисленные по данным [16, 17] и удовлетворяющие условию ${{c}_{j}}\left( T \right) = {{d{{h}_{j}}\left( T \right)} \mathord{\left/ {\vphantom {{d{{h}_{j}}\left( T \right)} {dT}}} \right. \kern-0em} {dT}}$. Второе слагаемое – это сумма теплот испарения свободной воды и реакций разложения гидроксидов и карбонатов кальция и магния где индексы k = 1, …, 5 соответствуют H₂O^free, Mg(OH)₂, MgCO₃, Ca(OH)₂, CaCO₃; $\Delta {{h}_{k}}$ – удельные теплоты разложения, значения которых представлены в табл. 1. Следует отметить, что, вообще говоря, эти теплоты зависят от температуры и давления в газовой фазе, но температуру можно не учитывать, поскольку в интервалах разложения относительные изменения каждой из $\Delta {{h}_{k}}$ не превышают 1%. В табл. 1 приведены значения $\Delta {{h}_{k}}$ и характерных температур при ${{P}_{{{\text{tot}}}}}$ = 1 атм. С ростом давления температуры разложения немного растут, удельные теплоты слабо убывают (~1% в диапазоне вариации p ~ 1–4 атм, рассматриваемом для ВРБ). Зависимость плотности от температуры необходима для расчета объемной теплоемкости бетона и энтальпии, используемой в численном расчете и для контроля баланса тепла.

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

В качестве проверки даны результаты расчета теплоемкости известкового бетона (“LimeStone/CommonSand”: L/S [10]) с составом, приведенным в табл. 2.

На рис. 2 показано сравнение рассчитанной температурной зависимости теплоемкости L/S-бетона по модели (10) с теплоемкостью, вычисленной по линейной аппроксимации из обзора [20]:

Результаты удовлетворительно согласуются до температуры ~1100 К. Расхождение при росте температуры, видимо, связано с выходом за границы применимости линейной зависимости (13) при T > 1000 К, когда бетон уже практически разложился.

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

Неопределенными параметрами теплофизики разложения остаются теплоты растворения компонентов в расплаве, которые зависят также и от состава расплава и должны определяться по экспериментальным данным.

МОДЕЛИРОВАНИЕ ПРОДВИЖЕНИЯ ГРАНИЦЫ ПЛАВЛЕНИЯ

Учет теплофизики нагрева и разложения бетона. Описанные выше превращения в бетоне протекают с поглощением тепла на некотором интервале температуры. Связь поглощенного тепла с температурой учитывается моделью эффективной теплоемкости, вводимой с помощью энтальпии материала H(T). Как и в (12), H(T) представляется в виде суммы двух слагаемых, одно из которых включает температурную зависимость удельной теплоемкости, вычисляемой по (10), второе – эф-фективную теплоемкость на температурном интервале данного превращения, учитывающую поглощение теплот превращений. В расчете формула (12) используется для единицы объема и учитывает также полную теплоту плавления:

Здесь $\rho (T)$ – текущая плотность бетона; $H_{L}^{{(k)}}$, k = 1, …, 6 – теплоты превращений при соответственно испарении воды, разложении Ca(OH)₂, Mg(OH)₂, CaCO₃, MgCO₃ и при плавлении/растворении всех компонентов в расплаве; T_ref  – температура отсчета энтальпии. Функции $\tilde {c}_{{(k)}}^{{}}$ заданы каждая на своем температурном интервале, вне которого они нулевые, а внутри имеют вид косинусоиды на отрезке [–π/2, π/2] с нормировкой под длину интервала $\Delta T_{L}^{{(k)}}$, соответствующего данной температуре превращения $T_{L}^{{(k)}}$. Интеграл от $\tilde {c}_{{(k)}}^{{}}$ по интервалу нормирован на единицу:

Плотность рассчитывается по (11) с учетом ухода газов, которые служат одним из источников массы в модели термохимии расплава. Процесс плавления/растворения компонентов в расплаве сложен и происходит на широком температурном интервале (до ΔT$_{L}^{{(6)}}$ = 1000 К [8]). Данные по температурам и энтальпиям растворения компонентов бетона в расплаве определенного состава редки, и модель, в которой плавление/растворение трактуется как один процесс со средним интервалом и теплотой превращения, широко применяется. Средние теплоты и интервалы плавления оцениваются по экспериментам.

Одномерная модель перемещаемой границы ВРБ. Расплав в контакте с бетоном прогревает его с границы. При постоянном потоке тепла от расплава распределение температуры в зоне прогрева бетона стабилизируется, и граница взаимодействия продвигается вглубь бетона с некоторой установившейся скоростью u_L. Суммируя затраты тепла на нагрев, разложение и плавление бетона от начальной температуры до полного плавления в некоторую эффективную теплоту плавления $H_{L}^{{{\text{eff}}}}$, вводя среднюю плотность бетона ${{\rho }_{{{\text{sol}}}}}$ и поток тепла f_in от расплава к бетону (Вт/м²), можно выразить поток через среднюю скорость продвижения границы:

При большой величине f_in ширина зоны прогрева расплава пренебрежимо мала относительно толщины бетонной стены или перекрытия, и соотношение (16) применимо. Но по мере снижения f_in, а также в случае большой теплопроводности твердого материала (стальная стенка шахты реактора ВВЭР) зона прогрева может быть широкой и неоднородной вдоль границы. Для учета этого решается полная задача теплопроводности в бетонной или стальной стенке, в которой заданный граничный поток тепла выражается в виде условия второго рода на границе расплав–стенка ∂Ω_m:

Здесь ${{\lambda }_{{{\text{sol}}}}}$– теплопроводность твердой фазы на границе контакта расплава с бетоном, взятая при некоторой эффективной температуре растворения T_liq. Задача решается методом конечных элементов. Затраты тепла на нагрев, химическое разложение, испарение воды и плавление бетона учтены в вычисляемой эффективной теплоемкости, фигурирующей в (14). Формула (16) используется для контроля баланса тепла при вычисленной скорости продвижения границы взаимодействия.

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

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

РАСЧЕТ МОДЕЛЬНОЙ ЗАДАЧИ ВРБ

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

– Одномоментное поступление 50 т железа в смеси с 20 т циркония при температуре Т = 2000 К, что моделирует проплавление корпуса ВВЭР и вытекание металлического слоя расплава (см. примеры расчетов в [6, 15]).

– Мощность остаточного тепловыделения ПД постоянна и равна 15 МВт.

– Состав бетона включает компоненты (в %): Fe₂O₃ (4.31), K₂O (3.74), Na₂O (1.40), CaO (14.36), MgO (3.85), SiO₂ (56.06), Al₂O₃ (6.98), CO₂ (1.37), H₂O (7.93).

– Площадь бассейна расплава соответствует площади пола бетонной шахты – 31.5 м².

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

Результаты приведены на рис. 6–10. Рис. 6 показывает полный источник тепла в расплаве, включающий источник от химических реакций и остаточное тепловыделение. В данном расчете постепенность исчерпания реагентов не учитывалась, и поступающие в расплав порции газа от разложения бетона (пар и углекислота) реагировали практически полностью в соответствии с активностями Zr и Si. Два участка роста и спада мощности скоррелированы с поэтапным исчерпанием Zr и Si, как это видно из графиков рис. 8. К моменту времени t = 12 000 с Zr и Si в расплаве окислены и источник тепла на графике близок к задаваемой мощности остаточного тепловыделения – 482 кВт/м² × 31.5 м² = 15 180 кВт.

Проверки модели продвижения границы проводились по выполнению соотношения (15), в котором полная теплота превращения (14) вычисляется в начале расчета, а средняя скорость продвижения границы определяется по положению крайнего узла сетки на границе взаимодействия. Уравнения (14)–(16) связывают работу трех вычислительных процедур:

– расчет теплоемкости и интегральной теплоты плавления бетона;

– расчет температуры в бетоне с использованием полученной теплоемкости;

– расчет продвижения границы при нагреве до критериальной температуры, равной 1400 К.

Степень согласованности работы этих процедур демонстрирует рис. 7, на котором изображена зависимость от времени скорости границы бетона, вычисленная двумя способами: по уравнению (15) через поток тепла в бетон и интегральную теплоту (14) и по движению границы в конечно-элементном расчете. Рост скорости связан с химическим тепловыделением в расплаве: тепловыделение в реакциях сначала циркония, затем кремния несколько ускоряет прогрев, что отчасти убыстряет плавление. В результате скорость продвижения границы медленно растет. Источник, связанный с реакцией с цирконием, более мощный, и, когда Zr исчерпан, поток тепла резко спадает, что замедляет и плавление. Тут вступает менее мощный источник, связанный с реакцией с кремнием, и картина повторяется, но с меньшим ростом скорости, до исчерпания кремния при t = 12 000 c. Мелкие колебания связаны с дискретностью аппроксимации положения границы. После начального установления профиля температуры (около 300 с) вычисленные скорости совпадают с достаточной точностью.

Последовательность реакций иллюстрируют рис. 8–10, на которых даны изменения во времени полных масс компонентов расплава (кмоли) во всей шахте. Для газообразных компонентов (водород, углекислота) указана их наработка. На рис. 8 представлены зависимости от времени содержаний компонентов системы, связанных с основными реакциями окисления и восстановления в расплаве, в которых участвуют вода, цирконий, кремний. До окисления циркония весь кислород воды уходит на реакцию с ним – высвобождается водород. Цирконий забирает кислород также у оксида кремния, и в расплаве появляется чистый кремний. После окисления циркония выходящий пар окисляет кремний, после окисления кремния – железо.

Вследствие существенно меньшей активности железа, чем циркония и кремния, его окисление в расплаве до исчерпания Zr и Si нулевое, а поступающий с бетоном окисел железа Fe₂O₃ полностью восстанавливается цирконием и кремнием, и до исчерпания кремния (t = 12 000 c) масса железа даже растет, а масса окисла FeO почти нулевая (рис. 9).

Что касается генерации водорода, то до полного окисления циркония и кремния весь пар отдает свой кислород, но, когда из металлов в расплаве остается только железо (менее активное), на выходе появляется вода, а темп генерации водорода снижается (рис. 8). Окисление железа паром – это единственная реакция, продуцирующая водород в заметном количестве на больших временах при ВРБ. Как видно из рис. 10, углерод, выходящий в виде двуокиси, в расплаве в основном переходит в монооксид СО. Водород и монооксид углерода – это все горючие газы, продуцируемые при ВРБ.

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

ЗАКЛЮЧЕНИЕ

Разработаны три модели для корректного расчета взаимодействия расплав–бетон:

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

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

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

Модель термохимии и теплофизики бетона исходит из естественных предположений о составе бетона и термодинамических данных по простым веществам, вследствие чего является универсальной и в рамках сделанных предположений не зависит от состава бетона. Температура и энтальпия растворения остаточных компонентов бетона в расплаве требуют дополнительных данных. Модель продвижения границы бетона вкупе с процедурой перестройки сетки позволяет проводить расчеты длительного взаимодействия расплав–бетон с малыми затратами процессорного времени (в приведенном случае – около 1 мин). Эти модели пригодны в качестве ядра уточненной интегральной методики расчета процессов при ВРБ.