Разработка кода КАНАЛ

Начиная примерно с 2000 года в США (NRC) с целью экономии сил и средств, вместо разработки двух основных расчетных теплофизических кодов, которыми являются RELAP5 [6] и TRAC [2,7] NRC перешла на развитие лучшего из кодов, которым был признан TRAC. В результате почти десятилетних работ значительно улучшилась достоверность расчетов, поскольку база теплофизических экспериментальных работ, на которые опирается код, выросла с 34 тестов (1997 г. версия TRAC-P [2]) до 500 (2001 г. версия TRAC-М [7]). В настоящее время число успешно описываемых интегральных и локальных экспериментов перевалило за 1000 (2007 г. вышла версия TRACE 5.0, тестируется версия 6.0, хотя детальное описание кода пока недоступно). Новая версия кода получила название TRACE и предназначена для выполнения расчетов по техническому обоснованию безопасности не только водо-водяных реакторов типов PWR, кипящих реакторов BWR, тяжеловодных реакторов типа CANDY, но и реакторов на быстрых нейтронах с натриевым и свинцово-висмутовым теплоносителями.

Код TRACE в процессе тестирования и модернизации силами 2 национальных лабораторий (Айдахо и Лос-Аламос) и 2 университетов получил: переработанную систему замыкающих соотношений, расширил карту режимов кипения, получил замыкающие в области теплообмена к жидким металлам, научился считать задачи водородной безопасности, прошел тестирование на задачах повторного залива, улучшился расчет процессов конденсации и кипения, модифицированы модели межфазного трения, изменены некоторые из моделей, например модель тройника, появились новые модели (струйный насос) и т.д. Кодом стало быстрее пользоваться за счет пре- и постпроцессора SNAP [9], способного переводить исходные файлы других кодов (например, Relap5) в исходный файл кода TRACE. Современные коды вычислительной гидродинамики (CFD) способны рассчитывать сопротивление каждого из компонент из групп последовательных сопротивлений при их взаимном влиянии друг на друга. Отдельные элементы контура с неизвестными коэффициентами сопротивления или двумя близко расположенными и влияющими друг на друга сопротивлениями, для которых нет данных в справочниках, теперь могут рассчитываться с использованием дополнительного CFD кода, с последующим использованием результатов вычислительного эксперимента для задания сопротивлений при моделировании контурной динамики. Это позволяет в перспективе значительно повысить точность расчета сопротивлений первых контуров. Аналогичная работа (проект моделирования ядерных реакторов Nuresim) выполнялась во Франции с кодом CATHARE [8], а CFD код для совместной работы с ним развивался под именем Neptune [10].

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

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

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

Прототипом кода КАНАЛ выбран код TRAC и его модификации [2,7]. В отличие от прототипа, в коде КАНАЛ в замыкающих соотношениях исправлены некоторые формулы, на основании полученных в расчетах расхождений с экспериментальными результатами. Эти исправления проводились по мере обнаружения проблем при решении практических задач. В основном использовались формулы отечественной теплофизической школы взамен устаревших зарубежных корреляций [5]. Используется адаптированная к отечественным энергетическим аппаратам система замыкающих соотношений [2] по трению и теплоотдаче, в частности, для коэффициентов теплоотдачи в режиме пузырькового кипения была применена формула Лабунцова Д.А. Вместо формулы Мартинелли, полученной на базе экспериментов на смеси воздуха и воды, применены отечественные нормативные опытные данные по трению, полученные на базе экспериментов на двухфазных пароводяных потоках [11]. Основной картой режимов кипения является карта, используемая в коде TRAC. Впервые применена универсальная карта режимов кипения для любого наклона труб (вызывается по ключу как вариант).

Код КАНАЛ разрабатывается на базе системы уравнений сохранения кода TRAC [2] и конечно - разностной численной схемы SETS [3] (двухшаговый метод повышенной устойчивости). В ходе разработки кода эта численная схема была переработана для ускорения кода.

Для решения разреженных систем линейных уравнений использовался пакет с открытым исходным кодом superLU [4], вслед за американцами, применившими пакет для нового кода TRACE. Этот пакет позволяет выполнять распределенные вычисления на современных многоядерных компьютерах.

В основе кода КАНАЛ лежат усредненные нестационарные уравнения для средних по объему ячейки величин, в которых для описания процессов отдельно рассматриваем законы сохранения для газопаровой и жидкой среды. Численная схема SETS используется для записи конечно - разностной схемы методом контрольного объема полностью неравновесной модели (6 уравнений), дополнительно включаются уравнения для неконденсируемого идеального газа и растворенной примеси (ортоборной кислоты).

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

Используемая система уравнений имеет вид.

Уравнение сохранения полной энергии смеси:

, (1)

здесь ? объемная доля газовой фазы или истинное объёмное паросодержание;

индекс g означают газовую смесь как фазу в целом;

индекс l означает жидкость как фазу;

? плотность газа;

? плотность жидкости;

e ? удельная внутренняя энергия;

P ? давление;

? вектор скорости жидкости;

? то же для парогазовой смеси;

? мощность, поступившая непосредственно в газовую или жидкую среду (без процесса переноса тепла теплопроводностью);

? удельная, на единицу объема среды скорость теплопередачи от стенки, отдельно к жидкости и газу.

Уравнение энергии газа (газовой фазы):

, (2)

здесь индекс v означает «перегретый» или «сухой насыщенный» пар;

qig ? межфазный теплоперенос за счет теплоотдачи между фазами

qgl ? тепло фазового перехода, за счет скрытой теплоты парообразования;

Г ? скорость межфазного массопереноса;

h'v ? удельная энтальпия пара в объеме (если пар конденсируется) или удельная энтальпия пара на линии насыщения (если жидкость испаряется).

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

, (3)

здесь использовано общепринятое соглашение о знаке скорости межфазного перехода: для уравнения (6) она положительна, а для уравнения (5) отрицательна.

Уравнение закона сохранения массы газовой фазы (пары воды и неконденсируемый газ):

. (4)

Уравнение закона сохранения импульса жидкой фазы:

, (5)

здесь Г- ? минимум между Г и 0 (защита от нефизического решения);

ci ? межфазный коэффициент трения;

cwl ? коэффициент трения о стенку (индекс l означает трение для жидкой фазы);

? вектор силы тяжести.

Уравнение закона сохранения импульса газовой фазы:

, (6)

здесь Г+ ? максимум между Г и 0 (защита от нефизического решения);

cwg ? коэффициент трения о стенку парогазовой фазы (индекс g ).

Система из шести уравнений (1)-(6) имеет шесть независимых переменных: давление (), истинное объемное паросодержание (), температуры газовой и паровой фаз (Tg, Tl), и скорости фаз ( ).

Уравнение сохранения полной массы парогазожидкостной смеси имеет вид:

. (7)

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

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

Закон сохранения массы неконденсируемого газа

. (8)

При описании парогазовой смеси принимается, что полная плотность газовой смеси и энергия являются суммой соответствующих компонент: водяного пара (v) и неконденсируемого газа (a).

, (8)

где ? плотность пара;

? плотность неконденсируемого газа;

. (9)

Предполагается, что справедлив закон Дальтона для смесей газов, поэтому

, (10)

где Pv ? парциальное давление пара, а величина Pa (давление неконденсируемого газа) фактически задает состав парогазовой смеси (объёмными долями).

Наличие растворимой примеси (ортоборная кислота) в коде учитывается введением дополнительного уравнения закона сохранения массы растворимой примеси, движущейся с потоком жидкой фазы. Обычно в первом контуре реакторов примесью является ортоборная кислота (orthoboric acid). Считается, что растворенная примесь никак не влияет на гидродинамику.

Уравнение сохранения растворимой примеси имеет вид:

, (11)

здесь переменная m - это величина концентрации примеси в жидкой фазе (масса примеси, отнесенная к массе воды).

Базовым элементом двухфазного кода КАНАЛ является модель трубы, которую необходимо тестировать.

На Рис. 1 показан тест разгона первоначально неподвижной двухфазной смеси в прямой горизонтальной трубе при заданном перепаде давления на концах трубы. Показано, что модель трубы кода КАНАЛ правильно описывает начальную стадию разгона трубы, когда ускорение двухфазной смеси в трубе соответствует аналитическому решению, следующему из закона Ньютона. Финальная стадия решения также правильно описывается кодом, о чем свиделельствует совпадение с аналитическим решением, полученным из уравнения Дарси-Вейсбаха при заданных коэффициентах сопротивления для каждой из фаз.

Размещено на http://www.allbest.ru/

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

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

Тестирование выполнялось для течения двухфазной смеси. Изотермическое течение было заменено расчетом теплоотдачи в трубе при фиксированных коэффициентах теплоотдачи. Следующим этапом было подключение карты режимов и замыкающих соотношений по теплоотдаче на основе системы замыкающих кода TRAC. Подключено 7-е уравнение, выполнена проверка решения для парогазовоздушной смеси. Подключено 8-е уравнение, выполнена проверка решения для растворимой примеси (борной кислоты).

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

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

Код КАНАЛ написан на языке С++ (код прототип TRAC - на языке Фортран). Код использует конечно - разностную численную схему SETS [3] (двухшаговый метод повышенной устойчивости). В ходе разработки кода было обнаружено, что при отключении второго шага эта численная схема работает в два раза быстрее при примерно тех же результатах расчетов.

Для решения разреженных систем линейных уравнений использовался пакет с открытым исходным кодом superLU [4], выполненная работа подтвердила работоспособность пакета superLU.

На Рис. 2 представлен вид одной (из серии) номограм, входящих в Нормативные методы расчета объектов котлонадзора. Набор этих номограм был отцифрован и переведен в форму, пригодную для интерполяции по трем независимым переменным - давление, степень сухости, удельный массовый расход (P,x,W). Коэффициент сопротивления двухфазного потока в трубе получается в этой методике расчета путем умножения трения по формуле Филоненко для на поправочную функцию, вид которой представлен на номограммах. Этот способ расчета основан на экспериментальных измерениях пароводяных смесей, а не водовоздушных, как имеет место при использовании формулы Мартинелли, широко использующейся в западных кодах.

Размещено на http://www.allbest.ru/

Рис. 2. Номограмма для поправочных коэффициентов для трения при течении двухфазной смеси по отношению к трению однофазной смеси рассчитанной по формуле Филоненко.

В настоящее время продолжается работа над дальнейшим развитием кода КАНАЛ.

В заключение автору хотелось бы поблагодарить людей, так или иначе повлиявших на облик кода КАНАЛ, в частности коллектив фирмы ЭНИМЦ «Моделирующие системы», г. Обнинск, где выполнялась данная работа, его директора к.т.н. Левченко В.А., н.с. Позднякова А.Н. (классы на С++ для термодинамических и транспортных свойств теплоносителей), Леонова И.Н. (подключение пакета superLU), коллег из ГНЦ РФ ФЭИ н.с. Сергеева В.В. (тестирование скорости работы кодов TRAC, Relap5, Корсар) и к.т.н. Зайцева А.А. (тестирование трубы кода КАНАЛ на обобщенных функциях), коллег из каф. ИТФ МЭИ, д.т.н. проф. Ягова В.В. (обсуждение замыкающих соотношений).

Список литературы

1. Сергеев В.В., Евстифеев В.В., Казанцев А.А., Поздняков А.Н., Трунов Н.Б., Фильчуков А.И., Ситник Ю.К. Верификация кода TRAC на гидродинамической модели трубного пучка ПГВ-1500 // Сб. докл. VП Международного семинара по горизонтальным парогенераторам, Подольск, 3 - 5 октября 2006 г., - ФГУП ОКБ "Гидропресс». - на 13 с.