Численная стабилизация с границы решений модельного одномерного реактора типа РБМК Текст научной статьи по специальности «Математика»

УДК 519.6

ЧИСЛЕННАЯ СТАБИЛИЗАЦИЯ С ГРАНИЦЫ РЕШЕНИИ МОДЕЛЬНОГО ОДНОМЕРНОГО РЕАКТОРА ТИПА РБМК

А. А. Корнев

i

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

Ключевые слова: численная стабилизация с границы, нейтронные поля, моделирование реактора типа РБМК.

The problem of construction of first king boundary conditions providing an asymptotic change of the trivial solution of a model one-dimensional RBMK reactor to the required stationary state is numerically studied according to specific features of this model. Results of calculations are presented for different admissible modes. The principal feasibility of efficient stabilization of the dynamics of occurring processes by boundary control of fast and slow neutrons is shown as well as its essential slow-down in the control of only fast neutrons.

Key words: numerical stabilization from the boudary, neutron fields, simulation of RBMK-reactor.

Постановка задачи. Рассмотрим следующую систему уравнений:

приближенно описывающую процесс эволюции нейтронных полей в модельном одномерном ядерном реакторе (см. [1]). Здесь искомыми являются функции Ф1^,:?;), Ф2(£,х) и C1(í,ж), задающие распределение быстрых, медленных и запаздывающих нейтронов в одномерной области со = [жь^г]-Отметим, что для физически корректных параметров модели однозначная разрешимость задачи для произвольных начальных условий из пространства ^(ш) на любом требуемом интервале времени следует из общей теории линейных уравнений параболического типа. В данном случае одномерность задачи, одна группа запаздывающих нейтронов С1 вместо стандартных шести и краевые условия первого рода выбраны для наглядности. При этом все полученные результаты могут быть однозначно адаптированы к реалистичным многомерным многокомпонентным моделям.

Далее мы будем использовать каноническую форму записи исходной системы (1) в терминах вектор-функции 17 = (Ф^Ф 2,С1)Т

(1)

с начально-краевыми условиями

Фй(0,ж) = Ф§(ж), С\0, ж) = СЦх) {<$>l{t,x),<$>2{t,x))\x=xk =¡i[k){t), Hk)(t) = (nlk)(t),txfk)(t)), к = 1,2,

(2)

— U = AVU, U\an =Kt) = 0«(i),/i(2)), U(0,x) = U0(x), полученную в результате переноса коэффициентов Vip в правую часть.

д

1 Корнев Андрей Алексеевич — доктор физ.-мат. наук, проф. каф. вычислительной математики мех.-мат. ф-та МГУ, e-mail: kornevQmech.math.msu.su.

Задача (1), (2) рассматривается в области со = [330,660] при следующих типичных (см. [1,2]) значениях параметров:

Vi = 2 • 107, V2 = 2, 2 • 105, Ai = 0, 012, Dx = 0, 4, D2 = 0,1, /?i = 0, 215 • Ю-3,

An = 0,015, A i2 = 0,149, A21 = 0,014, Ai = 0,009, A2 = 0,15,

а коэффициент A22 выбирается равным A22 = 0,138 25 5 205 475 8 44 из условия вырожденности оператора Ау правой части, что является необходимым для корректности рассматриваемой модели. Действительно, пусть /ü(iy2)(t) = 0. Тогда решение задачи на собственные значения Ау£ = не вызывает трудностей, так как все собственные функции имеют вид

£(n,i) = [k¡n,íy,kfníiy,kfníi)}T sin (^^ззд330^) ,П = 1,2,..., i = 1,2,3.

Следовательно, для определения п-й группы из трех собственных чисел г = 1,2,3, и коэффи-

12 3

циентов к,' \ для соответствующих собственных векторов получаем задачу на собственные значения

(га,г; ' '

с (3 х 3)-матрицей. В табл. 1 приведены приближенные собственные значения ^ для п = 1,..., 6,

12 3

г = 1,2,3. При этом вектор коэффициентов к^ отвечающий нулевому собственному числу Адд), имеет вид

[0, 994 912 667 312 951; 0,100 740 252 696 917; 0, 000 431169 097 015]т.

При увеличении п значения изменяются следующим образом: все собственные числа остаются действительными и отрицательными; первое собственное число п-й группы стремится к значению —Ai = —0, 012, которое является точкой сгущения спектра; два других монотонно убывают к минус бесконечности. При этом система собственных функций образует базис и как следствие решение U(t, х) задачи (1), (2) можно представить в виде

U(t,x) = ^ (¿)£(ra,í)(X>.

га,г

Отсюда в том числе следует, что исходная задача является асимптотически устойчивой и для произвольных начальных данных IJq{x) = U(0,x) решение U(t,x) в отсутствие ошибок округления стремится к стационарному состоянию Z(x) = (0)<^(1д)(ж) при t —> оо. Это соответствует общей физической картине исследуемых процессов. При этом, согласно структуре спектра, характерная асимптотическая скорость выхода решения U на стационарный режим Z для почти всех начальных условий Uo(x) не выше, чем Ce-0'012í. Это означает, что интегральная норма разности U — Z уменьшается приблизительно в три раза за 100 единиц времени.

Рассмотрим следующую задачу: найти такие краевые условия /i(i;2)(t), что решение исходной системы (1), (2) для IJq{x) = 0 асимптотически стремится к стационарному решению Z(x) = £(1,1)(ж)-По сути, нас интересует математический вопрос о возможности управления по краевым условиям решением приведенной системы уравнений для физически разумных параметров в классе положительных функций (так как только такие решения имеют в данном случае физический смысл), т.е. теоретическая возможность перевода модельного реактора из нулевого состояния в заданный стационарный режим посредством изменения значений для (Ф1, Ф2) только в крайних точках х = х\,х2. Отметим, что в работе речь идет об имитации узкого класса процессов, близких к стационарным, происходящих в реакторах типа РБМК. При этом исследуются только математические аспекты указанной задачи.

Метод решения. В соответствии с разработанным в [3-5] подходом искомые граничные функции будем строить следующим образом. Формально продолжим исходную задачу (1), (2) в расширенную область oj так, чтобы стационарное решение Z(x) расширенной задачи совпало на подобласти ш с интересующим нас стационарным решением Z(x) исходной системы. Далее, начальное условие IJq{x) = 0 системы (1), (2) продолжим в область Со до функции IJq{x) так, чтобы решение U(t, х) расширенной задачи сходилось при t —> оо к стационару Z(x). Тогда след функции U (t, х) на границе со дает искомое граничное управление 1,2) С^-)? обеспечивающее сходимость решения U(t,x) исходной задачи (1), (2) к требуемому стационару Z(x) при нулевой начальной функции.

11 ВМУ, математика, механика, №3

Таблица1 Старшие собственные числа матрицы Av

i Номер группы п

1 2 3 4 5 6

1 ю-'' -0,01146 -0,01179 -0,01189 -0,01193 -0,01195

2 -9 -212 -546 -1003 -1573 -2244

3 -332 303 -334 281 -337 582 -348188 -348188 -355 515

В данном случае удобно выбрать со = [0, 990] с нулевыми краевыми условиями

(ФЧ*,ж),Ф2(*,ж))|*=0>990 = 0.

Тогда все собственные функции оператора Ау расширенной задачи будут иметь вид

~ гГ1 у 2 ГЗ 1Т • {7ТТ1Х\

Цп,г) = [к(п,г) 5 (га,г)! к(п,г)\ 8111 ^"999") '

При этом все собственные числа оператора А-у действительные, а система собственных функций полна в пространстве решений. В табл 2. представлены приближенные собственные значения г = 1, 2, 3, п = 1,... , 6,15. Отметим, что собственная функция £(зд)(ж), соответствующая собственному числу А(3д) = 0, имеет вид

и значения коэффициентов те же, что и у соответствующей собственной функции £(1,1) (ж) исходной системы. Таким образом, стационарным решением расширенной задачи является произвольная функция вида [То (ж) = й(зд) х £(зд)(ж)) совпадающая при й(зд) = —1 с требуемым стационарным решением Z(x) на области со. Однако, как следует из табл. 2, полученная расширенная задача неустойчива. И если начальная функция [7(0, ж) содержит в разложении первую или вторую собственную функцию с ненулевыми коэффициентами, то норма решения неограниченно растет при Ь —У оо. При этом имеет место следующая

Лемма. Пусть полная система {(т(х) = т = 1,2,...} собственных функций

оператора расширенной задачи упорядочена по убыванию значений собственных чисел Х\ > Аг > Аз = 0 > А4 > ..а начальное условие имеет вид

оо

¿7(0,ж) = -|(3д)(ж) + ^ ат£т(ж), г0 ^ 0, «¿0+1 / 0.

т=г о+1

Тогда решение [7(¿, ж) стремится к функции ¿(х) при £ —> оо, если и только если ¿о ^ 3. При этом верна асимптотическая оценка

\\2-и\\ьф) ^ С • е-л>о+1*.

Таким образом, если начальное условие [То(ж) расширенной задачи имеет указанный в лемме вид, то соответствующее решение асимптотически сходится к требуемому стационару. Поэтому если искомые граничные условия для исходной задачи положить равными = и{Щэш, то соответствующее таким краевым условиям решение [7(£, ж) будет совпадать с [7(£,ж) на со, т.е. удовлетворять условию стабилизации.

Таблица2

Старшие собственные числа матрицы А\-

г Номер группы п

1 2 3 4 5 6 15

1 51,2 28,4 -4 • 10-1и -0,01017 -0,01112 -0,01146 -0,01193

2 -0.014 -0.016 -9 -62 -130 -213 -1573

3 -331718 -331937 -332 303 -332 816 -333 475 -334 281 -348188

Отсюда следует, что для построения управляющих краевых условий fx(t), обеспечивающих стабилизацию решения, достаточно продолжить нулевую начальную функцию из области со в область со = со/со до функции [То(ж) так, чтобы выполнялось указанное в лемме условие. Для этого выберем подпространство допустимых смещений С,0 = span < ¿¿(ж), % = 1,..., ¿о >• где

li(x) = W-1

0, ж € со;

До 1£г(х), X € со, г = 1,..., ¿0-

а вектор-функция грг(х) = А0 1^г(х), фí = [ф!, , по определению является решением следующей задачи:

,/V

-jgj- = X G ш, ФЦэо, = о, j = 1,2,3.

го

И далее будем искать начальную функцию в виде Uo(x) = ^^c^(.t) из условия вложения

г= 1

(Z(x) — Uo(x)) G Н- = span < im(x),m = ¿о + 1, го + 2,... > . Так как собственные векторы £т(х) рассматриваемой задачи неортогональны, то для нахождения

го

коэффициентов {ci,i = 1,.

, го} достаточно решить систему уравнений (Z(x) — ^^ Cili(x),f)j) = О,

г=1

j = 1,...,г'о, где Н^ = span < 'f]i(x),i = 1,...,г'о >-L Н- . Здесь f]i(x) являются старшими собственными функциями оператора А у, сопряженного к оператору А.у исходной задачи. Детальные пояснения к данному алгоритму и его применение, в том числе для нелинейных задач, можно найти в работах |6, 7|.

Приведем результаты расчетов для различных чисел ¿о- Так как для параметров рассматриваемой модели система уравнений (1) является специфически жесткой, то, следуя |2|, для построения численного решения будем использовать полунеявную по всем компонентам вектора U схему с весами при в = 0, 55. В данном случае методы типа расщепления требуют исключительно тонкой настройки соотношения шагов по пространству и времени и поэтому не подходят для наших целей. Расчеты проводились с шагом по пространству hx = 1, что, согласно |2|, на порядок меньше стандартного. Шаг по времени т варьировался в зависимости от поведения решения от Ю-7 на первом этапе стабилизации до Ю-4 на завершающем отрезке по времени. Более того, при расчетах расширенной задачи требовалось каждые 1... 10 шагов заново решать задачу проектирования (Z(x) — U(t,x)) G Н-, т.е. осуществлять стабилизацию с обратной связью. Данная процедура является допустимой и обязательной, так как обеспечивает правильное формирование искомых краевых условий. Отметим, что согласно специфике уравнений граничное управление и получаемое решение для х G со в процессе стабилизации должны оставаться положительными. Для решения промежуточных подзадач в процессе расчетов использовались пакеты ARPACK и SuiteSparse.

На рисунке изображены найденные управляющие условия /.t^(t), т.е. функция Ф1^, .т)|.т=ззо, для г'о = 3,10,

о 12

15. Графики для управляющих условий и /л^ не

приводятся, так как /.t^(i) ~ 0, l/.t^(t), а в силу сим-

1 2

метрии исходной задачи управления /.t^(i) при х = 660 имеют аналогичный вид.

В табл. 3 для типичных величин ¿о указаны приближенное значение для максимума управляющих функций ||/.1^||с и точность стабилизации

20 40 60 80 t

Вид управляющих функций fJ,\(t): 1 - го = 3, 2 - 10, 3 - 15

Т а б л и ц а 3

Точность стабилизации

Егг|Те = || Z-U(Te

\Ьф]

'¿0 Характеристики управления

IKDIIC Егг|те=50 Ursolic Err те=юо 11^ = 100 Не

3 0,14 6,54 0,49 3,34 0,74

10 0,66 6,28 0,51 3,07 0,75

15 90 5,52 0,57 2,62 0,79

с найденными управлениями для моментов времени Те = 50,100. Отметим, что при "выключении" управления с момента Те, т.е. при /л = 0 для £ ^ Те, решение исходной задачи асимптотически стремится к функции 11те(х) = с>'(1д)(Ге)^1д)(х). Соответствующее значение ||?7уе||с = °'(1Д) 1)

также приведено в табл. 3.

12 ВМУ, математика, .механика, .4« 3

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

Помимо представленных было подробно исследовано большое количество других решений, полученных для различных параметров алгоритма. При этом все управления имели следующую структуру (см. рисунок, а, Ь). На первом временном отрезке длины порядка 0(10-1) происходит быстрый рост функции U(t,x) в области ш за счет краевых условий, растущих в пределах от нуля до 0(10° ■■■2) в зависимости от режима. На следующем временном отрезке на порядок большей длины краевые условия убывают до величины порядка 0(10-1). А далее следует заключительный этап на 3-5 порядков большей длины, в течение которого краевое условие изменяется в пределах от 0(10-1) до 0(10"3-"5). Отметим, что на данном временном участке решение U(t,x) внешне меняется незначительно, однако именно на этом этапе коэффициент a(1;1)(i) при собственной функции £(1,1)(ж) в разложении U(t,x) достигает требуемого значения с приемлемой точностью, а коэффициенты при остальных гармониках уменьшаются до несущественных величин. Если процесс стабилизации остановить в некоторый момент t = Те до завершения последнего этапа, формально положив граничные условия равными нулю при t ^ Те, то решение за время 0(10-1) выйдет на стационар,

определяемый текущим значением коэффициента При этом, возможно, норма решения

уменьшится в несколько раз.

Как отмечалось, характерная асимптотическая скорость выхода на стационарный режим для типичных начальных условий определяется величиной порядка С • е—°>012i _ Поэтому в общем случае в рамках данного подхода для приведенных значений параметров теоретически невозможна стабилизация в интегральной норме к требуемому стационару со скоростью выше, чем указанная. Если же в расчетах ¿о ^ Ю, то, согласно табл. 2 и рисунку, Ь, данная скорость достигается практически с 5-процентной точностью.

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

Замечание. Хотя функции Ф2(£, х) и /х^ 2^(t) в процессе расчетов принимают на порядок меньшие значения, чем $1(i, х) и стабилизация посредством функций вида (//^ 2)(i),0) требует существенно больше времени, т.е. вклад /х^ 2)(i) в процесс управления имеет важное значение.

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

Работа выполнена при поддержке гранта РФФИ № 15-01-08023.

СПИСОК ЛИТЕРАТУРЫ

1. Reed Wm.H., Hansen К.F. Alternating direction methods for the reactor kinetics équations // Nuclear Sci. and Eng. 1970. N 41. 431-442.

2. Страховская Л.Г., Федоренко Р.П. Ядерный реактор — жесткая система. III. Квазиасимптотический метод интегрирования и главная спектральная задача (вычислительный эксперимент). Препринт N 44 Ин-та прикладной математики им. М.В. Келдыша РАН. М., 1996.

3. Fursikov А. V. Stabilizability of two-dimensional Navier-Stokes équations with help of boundary feedback control // J. Math. Fluid Mech. 2001. 3. 259-301.

4. Chizhonkov E. V. Numerical aspects of one stabilization method // Russ. J. Numer. Anal. Math. Modelling.

2003. 18, N 5. 363-376.

5. Ведерникова Э.Ю., Корнев A.A. К задаче о нагреве стержня // Вестн. Моск. ун-та. Матем. Механ. 2014. № 6. 10-15.

6. Корнев А. А. Об итерационном методе построения "усов Адамара" // Жури, вычисл. матем и матем. физ.

2004. 44, № 8. 1346-1355.

7. Fursikov A.V., К orne v A.A. Feedback stabilization for Navier-Stokes équations: theory and calculations // Mathematical Aspects of Fluid Mechanics. Lect. Notes Ser. Cambridge: Cambridge University Press, 2012. 130-172.

Поступила в редакцию 27.li.2015