О соответствии требованиям МАГАТЭ кодов RELAP-5 и Athlet
Дата: 16/03/2021
Тема: Атомная наука


Катковский Е.А. Руководитель экспертизы АО «ВО «Безопасность», К.т.н.

Новый документ МАГАТЭ SSG-2 ред.1 [1] «Детерминистический анализ безопасности атомных электростанций» (ДАБ) регулирует детерминистический анализ безопасности для нормальной эксплуатации, ожидаемых при эксплуатации событий, проектных аварий (ПА) и условий расширения проекта, включая тяжелые аварии, определения которых приводятся в документе SSR-2/1 (ред. 1) [2] и в Глоссарии МАГАТЭ по вопросам безопасности, и является важным средством подтверждения адекватности мер по обеспечению безопасности.


В данной работе рассмотрено соответствие требованиям документа [1] кодов RELAP-5 MOD3 (аттестационный паспорт Ростехнадзора №180 от 28.10.04) и Athlet-2.1 (аттестационный паспорт Ростехнадзора № 350 от 17.04.14), используемых в настоящее время в Российской Федерации при проектном обосновании безопасности АЭС.
Следует отметить, что более поздние версии кодов RELAP-5 и Athlet не имеют отличий в методических подходах от вышеуказанных версий и не изменяют выводов данной статьи. Кроме того, всё нижесказанное относится и к версиям кодов RELAP-5 3D и к Athlet CD.

Требования к ДАБ

Из всего объема требований к ДАБ в [1] рассмотрим требования раздела 5. «Использование компьютерных кодов для детерминистического анализа безопасности».

В п. 5.1 раздела 5 [1] указано, что Требование 18 документа МАГАТЭ GSR Часть 4 (ред. 1) [3] предусматривает, что «Любые методы расчета и компьютерные коды, используемые при проведении анализа безопасности, должны проходить верификацию и валидацию».

В п. 5.2 раздела 5 [1] указано, что в отношении выбора компьютерных кодов следует подтвердить, что:

  1. Физические модели, использованные для описания процесса, обоснованы.
  2. Упрощающие допущения, принятые в моделях, обоснованы.
  3. Корреляции, использованные для представления физических процессов, оправданы, а пределы их применимости установлены.
  4. Пределы применения кода установлены. Это важно, когда модель или метод расчета разработаны лишь для моделирования физических процессов в конкретном диапазоне условий, и код не следует применять за пределами этого диапазона.
  5. Численные методы, использованные в коде, точны и надёжны.
  6. Для разработки, написания, испытаний и документирования кода применяется систематический подход.
  7. Соответствие исходного кода его описанию в документации системного кода оценено.

В п.5.3 раздела 5 указано, что оценка точности отдельных компьютерных кодов должна включать ряд шагов:

  1. Определение важных явлений во вспомогательных опытных данных и ожидаемом поведении станции;
  2. Оценка неопределенностей, связанных с численными методами, использованными в коде;
  3. Оценка неопределенностей в основных моделях, использованных в коде;
  4. Выявление чувствительности важных процессов к значениям основных переменных.

В пункте 4.60 документа [3] указано, что верификация компьютерного кода должна включать верификацию модели и системного кода.

В п. 5.15 раздела 5 [1] указано, что в верификацию компьютерного кода следует включить подтверждение того, что код (исходный код и алгоритм) точно отражает математическую модель реальной системы (верификация модели) и соответствует документации кода (верификация системного кода).

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

В п. 5.16 [1] указано, что верификацию компьютерных кодов следует проводить средствами экспертизы, инспекции и аудита. Контрольные листы могут быть представлены для экспертизы и инспекции. Аудиту можно подвергнуть отдельные элементы для обеспечения качества.

В п. 5.17 [1] указано, что верификацию компьютерного кода следует проводить для экспертизы исходного кода на предмет его описания в документации кода. В верификацию следует включить экспертизу концепции проекта, основной логики, блок-схем, алгоритмов и вычислительной среды.

Анализ отступлений от требований [1] к ДАБ.

Сразу следует отметить, что требование 18 документа [3] в отношении верификации кодов не выполнено, и в различных руководствах по этим кодам никаких ссылок на верификацию нет.

Корректность постановки граничных и начальных условий

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

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

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

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

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

Рассмотрим постановку задач по моделированию гидродинамики в известных программах RELAP-5 [4] и Athlet [5] c позиции определений математической физики.

Рассмотрим в [4] раздел 3.1 Field Equations.
В подразделе 3.1.1  Basic Differential Equations приведена основная система дифференциальных уравнений в частных производных 3.1.2 – 3.1.12, которая представляет собой типичный пример уравнений математической физики.

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

Однако в разделе 3.1 вообще не встречаются словосочетания «boundary condition» (граничные условия) или «initial condition» (начальные условия) как будто в моделируемых физических объектах этих условий не существует!

Может быть в [5], которое датируется более чем на 10 лет позже, есть корректная постановка задачи.

Рассмотрим в [5] раздел 2 ATHLET Thermo-Fluiddynamics.
В подразделе 2.1 Basic Field Equations приведена основная система дифференциальных уравнений в частных производных 2.1-2.7, которая, как и в случае Relap-5 представляет собой пример уравнений математической физики.

И как в случае RELAP-5 в подразделе 2.1 вообще не встречаются словосочетания «boundary condition» (граничные условия) или «initial condition» (начальные условия)!

Подозрительное сходство в постановке гидродинамической задачи сразу наводит на мысль о подмене общепринятых в математической физике понятий своими, так называемыми «балансовыми соотношениями» (Momentum, Mass and Energy Balances) и далее, в обоих случаях [4 и 5], именно т.н. балансовые соотношения подменяют и граничные и начальные условия.

К чему же приводит игнорирование определений математической физики «граничные условия» в случае рассматриваемых кодов?

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

Авторы этих кодов понимают, что вместо законов сохранения необходимо внедрить в код некие «физические» соображения, иначе решить задачу моделирования не получится.

В подразделе 3.1.7  Volume-Average Velocities [4] читаем почему и как это делается:

«Предыдущая разработка разностных уравнений рассматривала трубу, состоящую из ряда односвязных объемов. В RELAP5/MOD3 каждый объем может иметь ноль, один или несколько переходов, прикрепленных к его входному концу, и ноль, один или несколько переходов, прикрепленных к его выходному концу. Таким образом, члены потока на входном или выходном конце объема состоят из суммирования по всем соединениям, присоединенным к этому концу объема.
Средние по объему скорости необходимы для расчета потока импульса, оценки сил трения стенки, оценки теплопередачи стенки, оценки межфазного теплопередачи и предела временного шага Куранта. В простом проходе с постоянной площадью среднее арифметическое между входом и выходом является удовлетворительным приближением. Однако для ветвящихся объемов с несколькими входами и/или выходами или для объемов с резким изменением площади использование среднего арифметического приводит к нефизическому поведению».

Осознав возникшую проблему с нефизичностью получаемых решений в коде, авторы начали эксперименты с подбором неких форм осреднения параметров в местах сочленения каналов:
переходя от RELAP5/MOD1 к RELAP5/MOD3 методом проб и ошибок разработчики экспериментально подобрали формулу осреднения параметров, которая, как утверждается на стр. 3-52 дала «желаемый результат», причем никакой оценки погрешности, вносимой этим действием в результаты расчетов, не дано.

В коде Athlet для замены реальных граничных условий разработчики пошли другим путем.  В разделе 2.3.1 Mass and Energy Balances [5] приводятся формулы осреднения в месте стыка двух разнородных каналов, отличные от RELAP5/MOD1 или RELAP5/MOD3. Главное отличие, что определение скоростей на границах перехода основано на принципе донорских клеток стр. 2−13 (The determination of the velocities at the junction borders is based on a donor-cell principle). Для модели Т-образной пристыковки разработана модель осреднения, описанная в разделе 5.11 The T-Junction Model. И снова никаких оценок погрешностей от применения формул осреднения вместо законов сохранения, а также влияния этих формул на конечный результат.

Налицо явные отступления от требований п. 5.2. (a, b, d, e) [1].

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

Общеизвестно, что для уравнений с разрывными коэффициентами постановка задачи Коши в общем виде не имеет решения. Никакого анализа корректности постановки задачи Коши в [4,5] нет, что является отступлением от требования п. 5.2 (d, e) [1].

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

  1. От любого численного метода требуется, чтобы он давал приближенное решение задачи с заданной точностью е > 0 за конечное число действий. Кроме того, требуется, чтобы схема была универсальной (пригодной для достаточно, широкого класса задач), однородной, устойчивой и экономичной (точнее, экономичным должен быть вычислительный алгоритм, применяемый для решения разностных уравнений). Все эти требования конкурируют друг с другом. Однородность схемы означает, что разностный оператор определяется одной и той же формулой во всех узлах сетки для любых коэффициентов и правой части уравнения из заданного функционального класса, а также для произвольной сетки.
  2. В качестве одной из априорных характеристик схемы используется погрешность аппроксимации. Максимально возможный порядок аппроксимации на решении оценивается в предположении достаточной гладкости решения. Локальная (в точке) аппроксимация не является, вообще говоря, адекватной характеристикой схемы. Естественной нормой для оценки погрешности аппроксимации является негативная («интегральная») норма. Заметим, что сравнение схем по порядку точности имеет смысл лишь при достаточно малых шагах h и t, при этом более высокая степень шага умножается на максимум производной решения более высокого порядка. Практически же используются весьма крупные сетки, на которых асимптотические оценки могут не работать. Может оказаться, что схема первого порядка точности на реальных сетках точнее схемы второго порядка точности. Таким образом, критерий выбора схем нужно формулировать так: схема должна давать достаточную точность на реальных сотках для рассматриваемого класса задач. С этой точки зрения следует учитывать, помимо математических аргументов (если они имеются), соображения качественного характера.
  3.  Одним из таких качественных требований является консервативность разностной схемы. В чем состоит свойство консервативности схемы? Для определенности будем говорить о задачах механики сплошной среды, описываемых уравнениями в частных производных. Приступая к решению задачи такого типа в некоторой области методом конечных разностей, вводят сетку в области G и заменяют дифференциальные уравнения разностными уравнениями для сеточных функций. В результате получают математическое описание дискретной модели среды. Очевидно, что дискретная модель должна отражать основные черты сплошной среды. Свойства сплошной среды определяются интегральными законами сохранения (количества движения, массы, энергии и т. д.) для любой подобласти G’, содержащейся внутри G. Дифференциальные уравнения суть следствия интегральных законов сохранения. Естественно требовать, чтобы разностные схемы выражали законы сохранения на сетке. Законы сохранения для всей сеточной области должны быть алгебраическим следствием разностных уравнений. Схемы, обладающие этими свойствами, и называются консервативными. Если 15—20 лет назад вопрос о том, должна ли быть схема консервативной («дивергентной»), мог служить предметом дискуссии, то в настоящее время мнение о необходимости требования консервативности является общепринятым.
  4. Наряду с консервативными разностными схемами применяются и консервативные дифференциально-разностные схемы. Для получения консервативных разностных схем можно использовать интегро-интерполяционный метод или метод баланса. Разностный оператор в пространстве сеточных функций должен сохранять основные свойства дифференциального оператора, заданного в пространстве функций непрерывного аргумента. Такими свойствами в линейном случае являются самосопряженность и знакоопределенность оператора. Неконсервативный оператор является и несамосопряженным. Недостатки, которыми обладают неконсервативные схемы, не устранимы. Сгущение сетки, на которое можно было бы рассчитывать для повышения точности, в случае неконсервативной схемы может даже увеличить ошибку схемы.

Никакого анализа таких ошибок в [4,5] не проводится, что является отступлением от требований п. 5.2 (e), 5.3 (b) [1].

В тексте документов [4, 5] нет никакого упоминания об анализе неопределенностей (Uncertainty) и чувствительности (Sensitivity), что также является отступлением от требований п. 5.3 (b, c, d) [1].

Требование п. 5.15 [1] в [4, 5] полностью игнорируется.

Рассмотрим, что предлагается в [4, 5] в качестве постановки начальных условий.
Очевидно, что задание начальных условий при нарушении правил задания граничных условий вызовет значительные проблемы. В такой сложной гидросистеме, как первый и второй контур водоохлаждаемого реактора имеется большое количество мест с т.н. местными сопротивлениями. Это места сужения или расширения потока, места поворота потока, устройства препятствия потоку (например: дистанционирующие решетки в ТВС), трубные пучки и т.п. Часто значения коэффициентов местных сопротивлений известны приближенно и их нельзя получить по формулам или таблицам из справочников. В таких случаях обязательно надо проводить анализ неопределенностей и чувствительности конечных результатов расчета по кодам на соответствие требуемым значениям стационарных условий. Особо надо отметить, что коэффициенты местных сопротивлений не являются постоянными величинами и зависят от многих параметров (геометрия каналов, схема соединения, характеристики потока, параметры среды). Например, коэффициенты местных сопротивлений при сужении или расширении потока (стык двух каналов) применяются при параметрах потока в меньшем сечении в месте соединения [7], а в кодах RELAP и Athlet используются другие параметры, т.к. в местах соединения параметры потоков не рассчитываются. Правомерность такого подхода нигде не обосновывается, что также является отступлением от требований п. 5.2 (a, b, e) [1].

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

А ведь такой важный элемент первого контура, как парогенератор (рис. 2 и 3), имеет протяженные U-образные пучки тонкостенных труб, где поток меняет направление на 180 градусов (Отступление от требований п. 5.2 (a, b, c) [1] ).


Рис. 1. Схема течения жидкости при повороте трубы.


Рис. 2 Схема расположения трубного пучка горизонтального парогенератора ВВЭР



Рис. 3 Схема расположения трубного пучка вертикального парогенератора PWR.

Многочисленные примеры валидации кодов RELAP-5 и Athlet полностью игнорируют силу инерции потока, движущегося по криволинейной траектории и силу инерции, связанную с Кориолисовым ускорением, да и сами экспериментальные установки, предназначенные для имитации первого и второго контура АЭС с водой под давлением (ВВЭР или PWR), не содержали U-образные пучки труб в моделях парогенераторов как в России (для горизонтальных ПГ), так и за рубежом (для вертикальных ПГ).

В нестационарных двухфазных потоках с течением в местах поворота потока вышеуказанные эффекты еще более усиливаются. Соответственно, карты режимов течения двухфазных потоков, которые встроены в коды RELAP-5 и Athlet, неприменимы! А это не только U-образные трубы парогенераторов, а любые места в первом контуре, где есть поворот потока – входы/выходы в корпус реактора, повороты в нижней и верхней камерах смешения, входы/выходы в коллекторах парогенераторов и т.п.

В обширном руководстве РБ-040-09 [8] охвачены практически все случаи теплообмена в ЯЭУ, но полностью отсутствует описание течения и теплообмена в U-образных трубах парогенераторов! Отчего такой пробел, и как практически решается задача расчета теплообмена и течения? Для однофазного потока задача расчета теплообмена и течения решаются методами 3D CFD, для двух- и более -фазных потоков вопрос требует особого исследования. Не взирая на это, кстати необоснованно, некоторые расчетчики, обоснующие безопасность ЯЭУ, пытаются создавать псевдо-трехмерные расчетные схемы конструкциями из одномерных каналов.

В кодах [3? 4] применяется т.н. «метод установления», когда запускается решение соответствующей задачи без возмущений. При этом «замораживаются» некоторые параметры, которые считаются заранее заданными. Например: расходы в петлях контура, напор насоса, тепловая мощность. Считается (причем бездоказательно!) что все параметры в системе приходят к правильному начальному состоянию. Но ведь даже простые изменения местных сопротивлений изменяют пространственное распределение потоков! В активной зоне (в одномерном приближении) распределение параметров потоков в каждой ТВС однозначно не будет соответствовать «желаемому» начальному распределению. Анализ чувствительности распределения потоков в ТВС от задания коэффициентов местных сопротивлений внутри реактора в указанных кодах не проводится – это отступление от требования п. 5.3 (d) [1].

В п. 6.7 [1] указано, что «необходимо определить параметры, к которым результаты анализа имеют наибольшую чувствительность. Анализ чувствительности должен проводиться с систематическим изменением ключевых исходных переменных для определения их влияния на результаты.

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

Важность анализа чувствительности к местным сопротивлениям, находящимся на тракте течения теплоносителя внутри корпуса реактора ВВЭР-1000 показана на диаграмме Рис. 4. Методика и результаты нейросетевого анализа распределения теплогидравлических параметров по каждой ТВС в зависимости от 50-ти местных сопротивлений и от координат конкретного сечения ТВС приведены в [9], там же приведена номенклатура этих сопротивлений и диапазон неопределенности их значений. Сверху вниз, последовательно по перечню местных сопротивлений из [9], показана относительная чувствительность результата распределения параметров течения в активной зоне. Три последних бара – это три координаты в каждой ТВС. Анализ показывает, что некоторые местные сопротивления оказывают значительное влияние на параметры потока в ТВС.

Интересно отметить, что именно анализ повышенной погрешности прогноза параметров течения в начальных и конечных координатах [9] показал, что нейронная сеть «чувствует» ошибочные граничные условия и сигнализирует об этом повышенной погрешностью. Этот факт подтверждает наличие проблем с законами сохранения при соединении или разветвлении однородных каналов в коде Athlet, по которому производилась подготовка обучающего множества нейросети
.
Однако, никакого анализа чувствительности, ни методики проведения такого анализа в документации [4, 5] на коды нет (отступление от требований п. 6.7) [1].


Рис. 4. Диаграмма результатов нейросетевого анализа чувствительности теплогидравлических параметров к местным сопротивлениям внутри корпуса реактора.

При анализе отступления от требований МАГАТЭ к ДАБ необходимо прояснить некоторые особенности численных методов, используемых в [3, 4] для решения систем уравнений гидродинамики.

Полный текст статьи читать здесь.







Это статья PRoAtom
http://www.proatom.ru

URL этой статьи:
http://www.proatom.ru/modules.php?name=News&file=article&sid=9584