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

МОДЕЛЬ, МЕТОД РЕШЕНИЯ И ПРОГРАММНАЯ РЕАЛИЗАЦИЯ В КОДЕ

В основу модели радиолиза воды в расчетном модуле AERMOD положена схема химических реакций в воде и ее растворах с примесями Н₂, Н₂О₂ и O₂, инициируемых воздействием ионизирующего излучения. Модель относится к гомогенному этапу процесса радиолиза, который наступает примерно спустя 0.1 мкс после начала облучения [1, 2]. К этому моменту времени области с первичными ионами и радикалами в воде (“шпоры”) приходят к относительному равновесию по составу ионов и радикалов и равномерно распределяются в объеме облучаемого раствора. Предполагается наличие однофазного замкнутого объема жидкости с одинаковой мощностью поглощенной дозы во всех точках.

Используемый набор реакций между продуктами радиолиза представлен в табл. 1 и соответствует данным из работы [3]. Схема реакций применима для воды с показателем кислотности pH в диапазоне от 1 до 10. Этот диапазон pH может быть расширен до 13 при добавлении реакций из табл. 2. Данные наборы сопровождаются константами, используемыми для вычисления скорости реакции в зависимости от температуры среды.

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

где ${{C}_{i}},$ ${{C}_{p}},$ $\;{{C}_{q}}$ – концентрация компонентов (здесь и далее частиц) сортов i, p, q, моль/дм³; τ – время, с; ${{G}_{i}}$ – выход частиц i-го сорта, моль/Дж; ${{I}_{\gamma }}$ – мощность поглощенной в воде дозы излучения, Дж/(кг ∙ c); ρ – плотность воды, кг/дм³; ${{\upsilon }_{{pq}}},$ ${{\upsilon }_{{pi}}}$ – константы скоростей возникновения и убыли частиц i-го сорта в реакциях между частицами p, q и p, i соответственно, дм³/(моль ∙ c). В зависимости от комбинации индексов i, p, q константы скоростей возникновения ${{\upsilon }_{{pq}}}$ и убыли ${{\upsilon }_{{pi}}}$ частиц приобретают значения равные константам скоростей прямой ${{\upsilon }_{f}}$ либо обратной ${{\upsilon }_{b}}$ реакции (см. табл. 1 и 2).

Данная система дифференциальных уравнений принадлежит классу “жестких” систем, для которых характерна возможность наличия одновременно для разных компонентов существенно различных временных масштабов изменения их концентрации. Это явление обусловлено большим различием в константах скоростей химических реакций и концентрациях компонентов. При использовании явных разностных схем для решения такой системы уравнений может потребоваться неприемлемо малый шаг по времени. Поэтому для решения был выбран неявный метод Розенброка 3-го порядка [4]. Этот метод является сравнительно быстрым, одношаговым, и для него в поставленной задаче существует возможность получить якобиан системы уравнений, что имеет решающее значение для эффективности и надежности данного расчетного метода.

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

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

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

Согласно модели из [5], формулы для потоков массы растворимой примеси через поверхность раздела фаз, кг/с, могут быть записаны следующим образом:

где ${{A}_{{int}}}$ – площадь межфазной поверхности, м²; ${{k}_{g}},{{k}_{l}}$ – коэффициенты массопереноса в газе и жидкости, м/с; ${{C}_{g}},{{C}_{l}}$ – концентрация примеси в газе и жидкости, кг/м³; $C_{g}^{{int}},C_{l}^{{int}}$ – концентрация примеси на межфазной поверхности в газе и жидкости, кг/м³; $E$ – скорость испарения воды, кг/с; $\rho _{g}^{{sat}}({{T}_{l}})$ – плотность насыщенного водяного пара при температуре воды ${{T}_{l}},$ кг/м³; ${{H}_{{cc}}}$ – коэффициент распределения примеси между водой и газом в отсутствие испарения (безразмерная величина, определяемая постоянной Генри). Второе слагаемое в первом уравнении отражает влияние испарения на межфазный массоперенос.

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

Очевидно, что в отсутствие испарения $(E = 0)$ система упрощается и переходит в классический вариант уравнений межфазного массообмена.

Коэффициенты массообмена получены в [5] при условии доминирования естественной конвекции вблизи поверхности воды. В этом случае коэффициент массопереноса для газовой примеси можно оценить из соотношения подобия с переносом молекул воды в паровой фазе:

где ${{D}_{g}}$ ‒ коэффициент диффузии данной примеси в газовой среде, м²/с; g – ускорение свободного падения, м/с²; ${{\rho }_{g}}({{T}_{g}})$ – плотность пара при его текущей температуре, кг/м³; ${{\mu }_{g}}$ ‒ динамический коэффициент вязкости водяного пара, Па ∙ с; $Sc = {{{{\nu }_{l}}} \mathord{\left/ {\vphantom {{{{\nu }_{l}}} {{{D}_{l}}}}} \right. \kern-0em} {{{D}_{l}}}}$ – число Шмидта (${{\nu }_{l}}$ – кинематический коэффициент вязкости жидкости, м²/с; ${{D}_{l}}$ – коэффициент диффузии жидкости, м²/с); ${{U}_{{surf}}}$ – скорость газа над поверхностью воды, м/с.

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

Следует отметить, что наиболее известными в России аналогичными программными средствами для расчета радиолиза воды при воздействии ионизирующего излучения являются код МОРАВА-H2 [6] и йодный модуль в составе кода КУПОЛ-М [7].

Код МОРАВА-Н2 предназначен для описания процесса образования и распределения водорода и других газов, а также технологических добавок теплоносителя в оборудовании I контура реакторных установок с ВВЭР в проектных режимах эксплуатации, включая проектные аварии. Программное средство может применяться также для анализа химического состава теплоносителя при разработке и обосновании норм водно-химического режима I контура и при расследовании инцидентов, связанных с его нарушениями. В расчетах учитываются следующие явления: радиолиз жидкой фазы теплоносителя, радиолиз водяного пара, пароциркониевая реакция и межфазный перенос веществ. Решение системы дифференциальных уравнений реализовано с помощью неявного метода Гира. При решении уравнений массопереноса используется метод Рунге–Кутты 4-го порядка. Валидацию модели проводили как на данных лабораторных экспериментов, так и на результатах измерений на действующих энергоблоках с ВВЭР.

Стоит отметить, что используемый набор реакций и реагирующих компонентов в коде фиксирован, и возможность его расширения или изменения без участия разработчика не предусмотрена. Базовые модули программы реализованы на алгоритмических языках Pascal, Fortran и Assembler. Погрешность расчета концентраций компонентов в элементах оборудования по программе МОРАВА-Н2 в сравнении с экспериментальными данными может достигать 50%.

Цель разработки йодного модуля в составе кода КУПОЛ-М – прогнозирование радиационных последствий, вызванных выбросами радиоактивного йода в помещениях защитной оболочки ВВЭР. Соответственно в коде моделируется большой перечень физико-химических явлений, связанных с поведением йода в водной и паровоздушной атмосфере. В коде реализована модель радиолиза водного теплоносителя, необходимая для правильного описания химических взаимодействий растворенного йода. Также большой блок кода посвящен расчету показателя кислотности pH теплоносителя. В коде реализована модель межфазного массопереноса для водорода, летучих соединений йода и других примесей. Для решения системы дифференциальных уравнений, так же как и в коде МОРАВА-Н2, используется метод Гира. Изначально модуль был создан с применением пакета MATLAB, однако для его интеграции в код КУПОЛ-М он был переписан с использованием языка Fortran 77. При этом за основу метода решения системы дифференциальных уравнений из свободной математической библиотеки SLATEC был взят пакет DASSL. Погрешность при расчете изменения показателя кислотности воды в сравнении с экспериментальными данными составляет 5%, а при расчете концентрации соединений йода варьируется в диапазоне от 10 до 35%.

Таким образом, каждый из указанных кодов имеет определенные преимущества и недостатки. Код МОРАВА-Н2 широко используется на практике и обладает большой валидационной базой в части механизмов радиолиза воды. При этом он является автономным с фиксированной системой химических реакций, что ограничивает область его применения для отличных от рассматриваемых в коде условий. Йодный модуль кода КУПОЛ-М используется в задачах оценки радиационной безопасности в помещениях станции. Его отличительная особенность – детальная модель физико-химического поведения соединений с йодом, не затрагивающая других продуктов деления, поэтому применение данного кода ограничено расчетами для условий в помещениях защитной оболочки водо-водяного реактора.

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

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

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

В соответствии с рекомендациями РБ-166-20 оценка погрешностей результатов расчетов, полученных по программе для ЭВМ, основывается на методике, описанной в стандарте ASME V&V 20 – 2009 [8], а также на методиках, разработанных в ИБРАЭ РАН [9]. При валидации для каждого эксперимента проводилось 100 многовариантных расчетов (данное число удовлетворяет требованиям указанных методик), в которых основные параметры варьировались случайным образом в диапазоне, определяемом с учетом экспериментальных погрешностей. Распределение параметров по интервалу предполагалось гауссовым. Данное число расчетов обеспечивает покрытие двустороннего интервала с вероятностью выше 0.9. В результате проведения многовариантных расчетов были получены коэффициенты чувствительности многомерной линейной регрессии к каждому из варьируемых параметров и затем представлены в виде гистограмм. При последующем статистическом анализе вычислялись минимальные, средние и максимальные значения моделируемых параметров.

Основным результатом проведенного соответственно методике [8, 9] анализа погрешности расчетных результатов является определение интервала полученных погрешностей [ε ‒ u_tot, ε + u_tot], где ε – отклонение результата расчета по программе для ЭВМ от экспериментальных данных, использованных для ее валидации; u_tot – неопределенность валидации, зависящая от неопределенности исходных данных, неопределенности измерений и неопределенности численных методов, реализованных в программе для ЭВМ.

Главной целью проведения эксперимента [10] являлось исследование влияния растворенных в воде водорода и кислорода на процессы образования и разложения пероксида водорода при гамма-радиолизе раствора воды. Для этого в эксперименте с помощью источника гамма-квантов ⁶⁰Co в специальной камере облучали герметичные стеклянные ампулы, содержавшие различные растворы вышеуказанных веществ. Проводился промежуточный анализ состава растворов на заданных промежутках времени. На рис. 2 показано изменение концентрации пероксида водорода ${{C}_{{{{H}_{2}}{{O}_{2}}}}}$ в растворе в зависимости от времени облучения $\tau $ при начальных концентрациях кислорода и водорода в растворе – 25.8 и 700 мкмоль/дм³ соответственно.

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

Изначально находящийся в растворе водород реагирует с образующимися из-за радиолиза OH-радикалами по реакции (1) с образованием атомов H. Затем, из-за наличия в растворе кислорода вследствие реакции (2) и начальной его добавки, атомарный водород преимущественно реагирует по (3) с образованием пероксильного радикала, поскольку скорость реакции (3) на два порядка больше, чем реакции (4). Как только весь кислород прореагировал, водород начинает интенсивно разрушать накопившийся пероксид, запуская цепные реакции (4) и (1).

Максимальное отклонение результатов расчетов от экспериментальных данных для максимальной концентрации пероксида водорода не превысило 5%, а для момента достижения максимальной концентрации отклонение составило 25%. Значения неопределенности валидации для указанных величин не превышают 19% (табл. 3).

Цель следующего эксперимента [11] состояла в исследовании радиолиза водных растворов, содержащих примеси пероксида водорода, водорода и кислорода при воздействии гамма-излучения. Для этого с помощью ⁶⁰Co облучали герметичные стеклянные ампулы с различными растворами указанных веществ. В процессе облучения проводился промежуточный анализ состава растворов на заданных промежутках времени. На рис. 3 для сравнения показаны расчетные и экспериментальные данные изменения концентрации веществ в системе в зависимости от времени облучения $\tau $ при следующих начальных значениях концентраций: ${{C}_{{{{{\left( {{{H}_{2}}{{O}_{2}}} \right)}}_{0}}}}}$ = 2.87 × × 10^‒3 моль/дм³, ${{C}_{{{{{\left( {{{H}_{2}}} \right)}}_{0}}}}}$ = 0.78 × 10^‒3 моль/дм³.

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

Максимальные отклонения результатов расчетов от экспериментальных данных для концентраций H₂O₂, H₂ и O₂ составили 2, 1 и 15% соответственно. Неопределенность валидации находилась в пределах 30%. Конкретные значения по каждому из проведенных экспериментов указаны в табл. 4.

Для проверки модели дегазации совместно с радиолизом использовались результаты исследования [12]. Мощность поглощенной дозы гамма-излучения в эксперименте составляла 2.5 Дж/(кг · с) во всей облучаемой камере. Облучение раствора проводилось в стеклянных ампулах объемом 20 мл. Площадь поверхности межфазного переноса газ – жидкость равнялась площади поперечного сечения ампулы и составляла 3.14 см². Ампулы закрывались алюминиевыми обжимными колпачками с силиконовыми мембранами для обеспечения герметичного уплотнения и предотвращения потерь газообразных веществ.

В представленном эксперименте атмосфера ампулы заполнялась аргоном, а вода деаэрировалась с заданным уровнем кислотности pH 10.6. В процессе измерений фиксировались концентрация водорода в газовом объеме над раствором и концентрация пероксида в жидкости. Для моделирования использовался полный набор реакций из табл. 1 и 2. На рис. 4 представлены результаты сравнения измеренной и рассчитанной концентраций водорода в газовом объеме в зависимости от времени облучения $\tau $ при соотношении объемов газовой и жидкой фаз 1 : 1.

Отклонение расчетных значений концентрации газообразного водорода от измеренных при данных условиях составило – 10%, а неопределенность валидации равнялась 31%.

Следует отметить, что после дополнительных исследований сходимости результатов расчетов при изменении временной и пространственной сеток было установлено, что неопределенность, связанная с использованием различных пространственных разбиений моделируемой области (одна, две, четыре ячейки), и неопределенность, обусловленная выбором максимального шага по времени (0.010, 0.005, 0.001 с), являются незначительными.

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

ВЫВОДЫ

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

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