CC BY
63
МЕХАНИКА
УДК 532.5
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ НЕСТАЦИОНАРНЫХ ТЕПЛОГИДРАВЛИЧЕСКИХ ПРОЦЕССОВ В КОНТУРАХ ЦИРКУЛЯЦИИ ВОДЯНОГО ТЕПЛОНОСИТЕЛЯ ПЕРСПЕКТИВНОЙ АЭС
© 2013 г. В.И. Будников1, О.Г. Савихин2, А.С. Чистов1
'НИИ механики Нижегородского госуниверситета им. Н.И. Лобачевского 2Нижегородский госуниверситет им. Н.И. Лобачевского
savihin@list.ru
Поступила в редакцию 23.05.2012
Рассматриваются разностные аппроксимации, основанные на явной и полунеявной схемах, системы уравнений сохранения, описывающей нестационарные теплогидравлические процессы в контуре циркуляции АЭС с водяным теплоносителем. Приводятся результаты тестовых расчетов применительно к перспективной АЭС с паровым энергетическим реактором.
Ключевые слова: численное моделирование, теплогидравлические процессы, контур циркуляции, водяной теплоноситель, одномерные законы сохранения, разностная схема.
1. Введение
Одной из главных задач, стоящих перед атомной энергетикой России на ближайший период, является обеспечение высокого уровня безопасности АЭС. Для обоснования безопасности, выполняемого средствами математического моделирования, обязательной компонентой является расчет нестационарных процессов, возникающих как в режимах нормальной эксплуатации, так и в различных аварийных ситуациях. С этой целью используются системные расчетные коды улучшенной оценки, которые характеризуются как широтой охвата моделируемого оборудования, так и совершенством и обоснованностью реализованных в них физикоматематических моделей. Наиболее известными зарубежными системными расчетными кодами, предназначенными для анализа переходных и аварийных режимов в АЭС с водяным теплоносителем, являются: TRAC, RELAP5 (США), CATHARE (Франция), ATHLET (Германия) [1]. Среди отечественных программ такого уровня можно отметить РАСНАР, БАГИРА, ТРАП, РАДУГА, КОРСАР [2,3]. Повышение точности при расчетном обосновании безопасности действующих АЭС в процессе их реконструкции и продления ресурса достигается ценой увеличения сложности математических моделей, объема программного кода и, как следствие, времени для получения результата.
В настоящее время в мире не вводятся в эксплуатацию АЭС с новыми типами реакторов. В то же время не прекращается разработка перспективных реакторных установок, которые используют принципиально новые решения, направленные на повышение безопасности и коэффициента полезного действия.
На начальном этапе проектирования, когда конструктивные решения и параметры АЭС постоянно корректируются, нецелесообразно использовать сложный программный код. Для получения консервативной оценки безопасности водо-водяных энергетических реакторов (ВВЭР) достаточно использовать гомогенные модели двухфазного течения теплоносителя. Это снижает зависимость программы от экспериментальных данных, делает ее более гибкой и позволяет быстрее получить ответ на конкретный вопрос при проектировании АЭС.
Целью данной работы является построение математической модели теплогидравлических процессов в контурах циркуляции АЭС с водяным теплоносителем, разработка численных методов ее решения, основанных на явной и полунеявной схемах, и проведение тестовых расчетов. В качестве модельных задач были выбраны следующие. Расчет волны сжатия воды в одиночном канале, а также расчет переходного процесса, вызванного ступенчатым возмущением теплового потока в активной зоне, и аварийного процесса, обусловленного
разрывом трубопровода, в перспективной АЭС с паровым энергетическим реактором [4].
2. Математическая модель теплогидравлических процессов
Разветвленный контур циркуляции представим в виде отдельных каналов, связанных между собой в точках смешения теплоносителя. Для описания теплогидравлических процессов в обогреваемом канале используем систему уравнений, содержащую одномерные законы сохранения энергии, массы, импульса [5]:
и2 и2
ф(8 + —) дG(I + —)
^--------— +----------— = QS,
Далее преобразуем систему уравнений сохранения к форме, удобной для разностной аппроксимации:
Sp— + О — - S(— + и —) = QS + ^ви2,
дt
дх
дt
дх
дt дх
х = о,
дt дх
+ ™ +г,ои + X дР = о,
дt дх дх
О = ирБ,
в= / - Р.
Р
Здесь приняты следующие обозначения: X -площадь сечения канала, р - плотность, е -внутренняя энергия, и - скорость вещества в канале, О - расход, I - энтальпия, Q - тепловой поток, t - время, х - пространственная координата, Р - давление, Е, - коэффициент гидродинамического сопротивления.
В качестве теплоносителя используется вода в различных фазовых состояниях. Это связано с тем, что для оценки безопасности необходимо проводить расчеты широкого круга глубоких аварийных процессов. Такие теплогидравлические процессы могут сопровождаться изменением давления от сверхкритического до атмосферного, а также переходом теплоносителя через все фазовые состояния: вода, пароводяная смесь и пар. Для описания пароводяной смеси применяется гомогенная равновесная модель [6, 7]. В результате дифференциальные уравнения для всех фазовых состояний имеют один и тот же вид, а все различие проявляется в уравнении состояния. Так как температура пароводяной смеси при нагреве не изменяется и равна температуре насыщения, то в качестве независимых термодинамических переменных выберем энтальпию и давление. Уравнение состояния имеет следующий вид:
Р = Р(Л1).
Поскольку в стационарном состоянии расход по контуру постоянен, его тоже удобно взять за независимую переменную.
,д1 , дР. дО
Х (р' ¥+ р Р ¥ > + ¥= °’ дО+дОи+ои+X дР = о,
дt дх дх
и = О.
рХ
Система уравнений дополняется начальными и граничными условиями.
В качестве начальных условий используется решение стационарной задачи:
1 (х,0) = 1 о (х), Р( х,0) = Ро (х),
О( х,о) = Оо( х).
Граничные условия задаются на входе в канал по энтальпии, а также по давлению либо расходу теплоносителя. На выходе канала задается давление либо расход:
1 (о, г) = 1х (О, О(о, Г) = Ох (О,
Р(1, t) = Рых (t).
Здесь I - длина канала, индекс «вх» относится к данным на входе в канал, а индекс «вых» - к данным на выходе из канала.
3. Метод численного решения 3.1. Уравнение состояния
Уравнение состояния строится на основе аппроксимации данных скелетных таблиц для воды и пара [8]. С целью повышения точности аппроксимации уравнение состояния задается в виде зависимости удельного объема V от давления и энтальпии:
V = V (Р, I):
1
где V = —.
Р
В области пароводяной смеси уравнение состояния используется в следующей форме [7]:
V"(Р) - V'(Р)
V = V'(Р) + У (Р) У (Р) (I -1'(Р))
I' (Р) -1' (Р)
при I' < I < I". Здесь переменные с одним штрихом относятся к воде на линии насыщения, а с двумя - к пару на линии насыщения.
Поскольку скелетные таблицы используют прямоугольную сетку по температуре и давлению, потребовалось с помощью линейных методов преобразовать ее в сетку по энтальпии и давлению. Для нахождения удельного объема и температуры по произвольным значениям энтальпии и давления применялась билинейная аппроксимация. Детальное исследование трех-
мерных поверхностей, построенных на основе скелетных таблиц, позволило сделать следующий вывод. Несмотря на то, что поверхности являются гладкими, производные плотности по энтальпии и давлению имеют разрывы. Эти разрывы носят не физический характер. Они являются следствием аппроксимации, использованной при составлении скелетных таблиц, поскольку ее целью являлось получение наиболее точных значений термодинамических функций, а не их производных.
В качестве одного из способов преодоления трудностей при численном моделировании было выбрано использование скелетных таблиц для самих производных термодинамических функций. Такие таблицы были построены на основе исходных данных с помощью вычисления соответствующих разностей в соседних точках [9].
3.2. Явная разностная схема
Для численного решения системы уравнений была выбрана аппроксимация при помощи явной разностной схемы «крест» [5]. В этой схеме используется левая разность по координате на п-м временном слое для аппроксимации производной по расходу и правая разность на (п+1)-м слое для аппроксимации производной по давлению:
РП (1 +
(р I) п
ТИ+1 _ тп к к к Г'П (Т т2\п
рП (рР)
А/
= С *оп (и2) п +
Рп _ Рп
гк гк-1
Ах
гп - гп оп - Оп
- оп * 1к-1 к к-1 ,
Рп+1 _ рп
- (РР )к к к
А/
Ах
оп - ои
Як Ах
+ (Р I)
тп+1 тп Г \п 1к - 1к
А/
Оп+1 Оп ^~пттп Г2.п ТТп
А/
к + опип - о;_1u;_1 +
Ах
Т)п+\ рп + 1
+ Сои + як к+1 ~ к = о.
Ах
Последовательность шагов при вычислении состоит в нахождении значения энтальпии на верхнем временном слое из уравнения энергии и подстановки его в уравнение неразрывности. Аналогично из уравнения неразрывности находится значение давления на верхнем временном слое. Оно используется в уравнении импульса для нахождения расхода.
Условием устойчивости явной разностной схемы является условие Куранта, связанное со скоростью звука в теплоносителе.
3.3. Полунеявная разностная схема
С целью расширения области устойчивости конвективные слагаемые, связанные с перено-
сом вещества, аппроксимируем по явной схеме, а слагаемые, связанные с быстрыми акустическими процессами, - по неявной.
Уравнение тепловой энергии запишем в форме, которая используется в интегральной модели количества движения [7]:
дI дI дР
Яр— + О— = я— + QS.
д/ дх д/
Это позволяет применить для решения разностных уравнений скалярную прогонку вместо векторной. Интегральная модель служит для расчета процессов, характерное время которых определяется временем прохода канала теплоносителем. Основное упрощающее предположение интегральной модели заключается в том, что перепад давления мал по сравнению с величиной самого давления в канале. Поэтому в контурах циркуляции с большой разностью температуры теплоносителя вкладом в тепловую энергию слагаемого, связанного с перепадом давления, можно пренебречь:
РкЯк
тп+1 тп тп тп
к к , г^п к к-1 о глп . гг
------------+ О----------------------= яМк + \
А/
Ах
+1 +1
гк гк-1
А/
(рР )к к
Рп+1 -рп тп+\ тп /"'тп+1 /"'тп + 1
-* /, - Р , 9 \
+ (рI)к
А/
+1 +1 +1
Г \п -^к - 1к , Ок - Ок-1 = о
А/
оп+1 - оп + оип - о;_lU;_l +
А/
Ах
Т)п+\ рп + 1
+ Сои + к+1 ~ к = о.
Ах
Выразим расход на верхнем временном слое из уравнения сохранения количества движения:
оп+1 = Ак (Рп+11 - рп+1) + рщ ,
А/
где Ак = як~,
Ах
А/
Fkn = оп - А/с ои-—(ои - о^и^У
Ах
Подставим полученное выражение, а также производную энтальпии по времени из уравнения энергии в уравнение сохранения массы и сгруппируем слагаемые:
АкРп++ + Щ - 2Ак )Ркп+1 + АР++11 = ЩРп + Fkn_l -
- F; + (рЖ (оп (I; - I;_l) - опяк Ах),
Рк
Ах
где щп = ^~ як (р'р )к.
А/ р
Приведенное выражение является линейным уравнением относительно давления с трехдиагональной матрицей. Для его решения используется метод скалярной прогонки.
Рис. 1. Вид переходного процесса при расчете по явной Рис. 2. Вид переходного процесса при расчете по не-схеме явной схеме
О
пг
щн
Рис. 3. Схема двухконтурной АЭС
4. Тестовые расчеты
4.1. Одиночный канал
Рассматривается труба без обогрева, заполненная водой при закритическом давлении. На входе трубы задано постоянное во времени возмущение по давлению, равное одному проценту от стационарного значения. На выходе трубы задано граничное значение скорости теплоносителя. Нулевое значение скорости соответствует закрытой части трубы.
Расчеты производились с различными шагами по времени, а также с различным числом расчетных точек по координате. Через равные интервалы времени отслеживалось распределение давления по координате. Вид переходного процесса при расчете по явной схеме показан на рисунке 1, а по неявной - на рисунке 2.
Из приведенных рисунков следует, что форму волны две разностные схемы передают примерно одинаково. При использовании явной разностной схемы фронт волны получается более крутой, но за ним следуют связанные с разностной аппроксимацией колебания большей амплитуды. В случае неявной схемы увеличение шага интегрирования по времени за границу области устойчивости явной схемы приводит к сильному сглаживанию фронта волны. Несмотря на то, что численный метод, ос-
нованный на полунеявной схеме, остается устойчивым, расчет такого вида переходных процессов при больших шагах интегрирования по времени теряет смысл.
4.2. Расчет переходного процесса в перспективной АЭС с паровым энергетическим реактором
Описание перспективной АЭС с паровым энергетическим реактором приводится в статье [4]. В данной работе используется ее упрощенная расчетная схема, представленная на рисунке 3. Здесь приняты следующие обозначения: АЗ -активная зона, ПГ - парогенератор, ГЦН -главный циркуляционный насос.
В расчетной схеме двухконтурная АЭС представлена шестью каналами. Второй контур представлен одним каналом парогенератора.
Для данной схемы произведен расчет переходного процесса, вызванного ступенчатым возмущением теплового потока в активной зоне. Возмущения такого рода используются при исследовании нейтронно-физической устойчивости стационарного режима. На рисунках 4, 5 приведены результаты расчетов по явной и по-лунеявной схемам. Синим цветом на графиках отмечены данные, полученные по явной схеме, а сиреневым - по полунеявной схеме.
Рис. 4. Энтальпия на выходе активной зоны Рис. 5. Энтальпия на выходе парогенератора
Рис. 6. Зависимость расхода теплоносителя от време- Рис. 7. Зависимость температуры теплоносителя от
ни в различных точках контура циркуляции времени в различных точках контура циркуляции
Из рисунков следует, что форма кривых для двух вариантов одинакова. Происходит переход на одно и то же новое стационарное значение. Результаты расчета с помощью полунеявной разностной схемы имеют опережающий по времени характер. Это, по-видимому, связано со свойством неявных схем увеличивать скорость передачи возмущения за счет уменьшения крутизны фронта. График энтальпии на выходе парогенератора демонстрирует различие в величине запаздывания, связанного с проходом теплоносителя по контуру циркуляции. Последнее обусловлено принятым упрощающим предположением при записи уравнения энергии, которое используется для полунеявной схемы. В целом, разность значений данных, полученных с помощью двух методов, находится в пределах одного процента.
4.3. Расчет аварии с разрывом трубопровода
Для того чтобы продемонстрировать возможность применения приведенной тепло-
гидравлической модели для расчета быстрых аварийных процессов, была рассмотрена авария с разрывом трубопровода на входе в активную зону (рис. 3). Для оценки допустимого времени задержки срабатывания систем аварийной защиты реактора необходимо знать такие параметры переходного процесса, как расход и скорость увеличения температуры теплоносителя в активной зоне.
В работе рассмотрен случай постепенного разрыва трубопровода, что соответствует появлению мелкой трещины и нарастанию её со временем. В математической модели в качестве граничного условия в точке разрыва задано давление, которое уменьшается со временем по линейному закону до атмосферного. В некоторый момент времени скорость истечения теплоносителя достигает скорости звука. В этом случае осуществляется переход к граничному условию на равенство расхода критическому значению.
Авария с разрывом трубопровода на входе в активную зону сопровождается уменьшением расхода теплоносителя в активной зоне до ну-
левого значения, а затем происходит опрокидывание циркуляции. В связи с этим для обеспечения устойчивости численного метода в разностную схему были внесены изменения. При изменении знака расхода в некоторой точке во всех конвективных слагаемых разностной системы уравнений левая разность по координате менялась на правую. Поскольку целью расчета являлось определение параметров протекания быстрых процессов на начальной стадии аварии, для вычисления использовалась явная разностная схема.
На рисунках 6, 7 представлены результаты расчета аварийного процесса в виде зависимости расхода и температуры теплоносителя от времени в различных точках активной зоны и примыкающего к ней трубопровода. Здесь синим цветом обозначены параметры в середине трубопровода, фиолетовым - в трубопроводе перед точкой разрыва, желтым - в начале активной зоны, голубым - в середине активной зоны, коричневым - в конце активной зоны.
Проведенные расчеты показывают, что численный метод правильно отражает процесс опрокидывания циркуляции и позволяет определить скорость роста температуры теплоносителя в активной зоне. Полученные результаты носят демонстрационный характер. Для точной оценки времени достижения критических значений необходимо дополнить математическую модель системой уравнений нейтронной кинетики и уравнением теплопроводности для тепловыделяющего элемента. Только в этом случае можно будет делать вывод о достаточности времени для срабатывания системы защиты реактора.
5. Заключение
В работе рассмотрены две разностные аппроксимации, основанные на явной и полунеяв-ной схемах, системы уравнений сохранения, описывающей теплогидравлические процессы в контуре циркуляции водяного теплоносителя.
Проведенные тестовые расчеты показали работоспособность численных методов, а также близость получаемых результатов. Из вычислительных экспериментов следует, что две схемы имеют примерно одинаковые затраты машинного времени на один шаг интегрирования по времени. При этом полунеявная схема имеет более широкую область устойчивости. Граница устойчивости определяется условием Куранта, связанным со скоростью теплоносителя, а не со скоростью звука в теплоносителе.
Список литературы
1. Девкин A.C., Мелихов О.И., Москалев A.M. и др. // Зарубежные теплогидравлические коды улучшенной оценки. Опыт разработки, создание и применение. М.: ОЦРК при Минатоме РФ, 2000. 176 с.
2. Мигров Ю.А., Волкова C.H., Юдов Ю.В. и др. КОРСАР - теплогидравлический расчетный код нового поколения для обоснования безопасности АЭС с ВВЭР // Теплоэнергетика. 2001. № 9. С. 36-43.
3. Нигматулин Б.И., Мелихов О.И., Соловьев С.Л. Состояние и развитие отечественных системных теплогидравлических кодов для моделирования аварийных и нестационарных процессов на АЭС с ВВЭР // Теплоэнергетика. 2001. № 3. С. 17-20.
4. Семченков Ю.М., Сидоренко В.А. Перспективы развития АЭС с ВВЭР // Теплоэнергетика. 2011. № 5. С. 1-8.
5. Рождественский Б.Л., Яненко Н.Н. Системы квазилинейных уравнений и их приложения к газовой динамике. 2-е изд. М.: Наука, 1978. 592 с.
6. Кузнецов Ю.Н. Теплообмен в проблеме безопасности ядерных реакторов. М.: Энергоатомиздат, 1989. 296 с.
7. Сабаев Е.Ф. Системы сравнения для нелинейных дифференциальных уравнений и их приложения в динамике реакторов. М.: Атомиздат, 1980. 192 с.
8. Риквин С.Л., Александров А.А. Термодинамические свойства воды и водяного пара: Справочник, 2-е изд. М.: Энергоатомиздат, 1984. 80 с.
9. Савихин О.Г., Линник С.В., Савихин А.О. Аппроксимация производных термодинамических функций для воды, пара и пароводяной смеси // Вестник Нижегородского университета им. Н.И. Лобачевского. 2011. № 4 (3). С. 1090-1091.
NUMERICAL SIMULATION OF NONSTATIONARY THERMAL-HYDRAULIC PROCESSES IN THE WATER COOLANT CIRCULATION LOOPS OF AN ADVANCED NUCLEAR POWER PLANT
V.I. Budnikov, O.G. Savikhin, A S. Chistov
The article considers explicit and semi-implicit finite-difference approximation schemes of one-dimensional conservation laws describing nonstationary thermal-hydraulic processes in the water coolant circulation loops of a nuclear power plant. The results of test calculations for an advanced nuclear power plant with a steam power reactor are presented.
Keywords: numerical simulation, thermal-hydraulic processes, circulation path, water coolant, one-dimensional conservation laws, difference scheme.
