Решение задачи расчета реактивности ядерного реактора с учетом ее некорректности

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

Обработка результатов измерений, основанная на решении уравнений нейтронной кинетики относительно реактивности к (0, связана с численным дифференцированием функции плотности нейтронов и (г), заданной значениями щ. Решение этой задачи оказывается неустойчивым. Попытка

создания компьютерных систем определения реактивности, без учета некорректности задачи, привела к сбоям работы ЭВМ из-за случайных погрешностей в измерениях и, [1]. Для борьбы с помехами в дополнение к аппаратным средствам применялся программный способ защиты, который по сути дела бесконтрольно корректировал исходные данные. Упрощения [2], например, предположение о том, что вклад скорости изменения плотности нейтронов в реактивность незначителен, оказываются грубыми для нестационарных режимов.

Для повышения устойчивости счета необхо

димо использование разностных схем, учитывающих некорректность задачи. Имеющиеся числен

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

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

Уравнения нейтронной кинетики имеют вид:

б

dt dC, dt

L

Начальные условия для Н)

*4 Cv^ = ^JXvLy0,

где C_v(f) - концентрация ядер - предшественников запаздывающих нейтронов v-й группы; Ел, P_v, X_v, L - константы.

Определение реактивности заключается в восстановлении К, в точках /„ 0 <  i < m, образую

щих разбиение t = 0 <  t_x< t₂< ... < t_m с учетом измеренных значений плотности нейтронов п, в тех же точках t;. При этом известно, что

(лЮ-д! <8,, (4) где 5, - уклонение измеренной плотности нейтронов П; от точной величины n(tj).

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

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

Через вектор, составляющими которого являются все измерения плотности нейтронов, поступившие до момента времени /,₊₁, можно записать оценку вектора реактивности Кг ср (и"¹). В момент времени /,₊₁ появляется информация, составляющая вектор измерений «^,и = (и⁽¹⁾, и^^), по которому определяется оценка K_l+1= Ф (и⁽'⁺¹⁾)- Рекуррентный подход заключается в построении алгоритма вида

К₁₊₁ =\у(К,,п,_+х) (5)

Система уравнений (1>—(3) разрешима с помощью одношаговых (Эйлера, Рунге-Кутта) или многошаговых (Адамса, Гира) методов параллельного типа: каждое уравнение системы на каждом шаге решается отдельно. Так, концентрация ядер -предшественников запаздывающих нейтронов, рассчитанная для z+1-го шага, в счет плотности нейтронов на этом шаге не входит и применяется предполагаемое значение данной величины. Параллельные методы дают возможность проводить вычисления на нескольких процессорах, однако они не пригодны для решения обратных задач.

Для этой цели необходимы последовательные методы численного решения систем дифференциальных уравнений. На основе метода работы [4] можно построить новую разностную схему для

системы (1НЗ), предусматривающую использование результатов вычисления по каждому предыдущему уравнению во всех последующих уравнениях и позволяющую на ее основе перейти к решению обратной задачи:

= С,

к - 2Z И' ^+1;

«,+1

т

^-^(Су.+т^) . (6)

21и,₊1 L

"“VJ + l

«,+1         )

■(7)

где т - шаг;

■                       -ру 6■

F₀ =--- + £>_VC_V, ;

LV=1

Fv=^n,-XvCv.. '

Формула (7) содержит операцию численного дифференцирования функции »,(/,), заданной приближенно. Разность и,+) - И; при т —> 0 может стремиться к нулю «несогласованно» с т. Тогда возможно, что отношение (и,+] - и,) / т —> со при т -> 0.

Следовательно, решение, полученное по формулам (6) - (7), может быть неустойчивым. Для построения класса приближенных, но устойчивых решений заменим в формуле (7) численное дифференцирование функции и„ полученной в результате измерений, расчетом производной некоторой непрерывной функции n(t)>, обеспечивающей наилучшее среднеквадратичное приближение к измеренным значениям п,. При этом потребуем выполнения условия (4). Путем дифференцирования многочленов наилучшего среднеквадратичного приближения могут быть получены следующие формулы [5]:

«,У ”     (3«,+, +1 Он, -18и,ч + 6и,_2 - и,„3); (8)

12т

«А 1(^-2 «4+2 ) ~ $(«4-1 ~«А+1ХВ       (9)

12и

” 1^^КП/,+3 ~^6пр^^{2 +18}Vi -¹⁰”д -³Vi) ,(1°) где т ,h, q- интервал дискретности; i, к, р - номер интервала.

В нестационарных режимах, особенно при прохождении экстремумов, результат измерения п_т может значительно уклониться от кривой среднеквадратичного приближения, построенной по предыдущим точкам. Соответственно производная, вычисленная, например, по формуле (8), будет существенно отличаться от точной величины (и,₊₁ - н,) / т. Потребуем, чтобы это отличие не вносило в расчет погрешность, превышающую погрешность измерения плотности нейтронов:

[и,+TH,^Z -и_/+1|<6,₊₁.            (И)

Алгоритм расчета производной с требуемой точностью построим на основе деления шага дис кретности т и последовательном применении формул (8)—(10). Деление шага дискретности позволяет ввести на интервале опроса измерительных устройств [;, г-1] дополнительно несколько фиктивных точек, имитирующих измерения. Введение имитаций измерения дает возможность усилить влияние замера плотности нейтронов в /+1-й точке на расчет производной. При этом используется особенность каждой из формул (8)(10), заключающаяся в различном влиянии на значение производной точки слева и справа от расчетной. Так, если принять т = 2й и в формуле (9) совместить индексы опроса г = к, г +1 = к + 2, а имитацию измерений выполнить в виде

«г-1 + п,

■      ^-2 = «,-1; «4-1 = —-—; «а = «г;

и, + и,.,

«А+1 =---"---: «А+2 = «,+1 ’         О2)

то с применением формулы (9) влияние замера в z+1-й точке на расчет производной усиливается по сравнению с формулой (8) с помощью имитации измерения п^х. Соответственно при т=3д и использовании формулы (10) с имитацией замеров пР-\ = п,ч+2В (п, -   ); пр=п,;

«р^М^'Ч^-^У        (13)

«р+2 = «, + 2 / 3 («,+] - п,); пр+3 = ni+v, дополнительно к одному реальному измерению вводятся две имитации измерения на интервале опроса датчиков. При этом вклад результата измерения пк+\ в формирование производной возрастает еще больше по сравнению с формулами (8) и (9).

В условиях физического эксперимента реактивность имеет конечную производную. Разностный аналог реактивности может быть представлен в виде к1+1 = к,+тк;.           (14)

На основе алгоритма общего вида (5) применение формулы (8) в выражении (14) дает рекуррентную формулу для расчета реактивности:

^/+1=^ +

+ ™(3^,+1 +Ю^, +18^.] +6^_2 -^_з)> (15)

где k_i+x определяется из формул (6), (7) исходя из и,+1 и производной п' по методу алгоритмической регуляризации:

шаг 1 - по замерам плотности нейтронов и формуле (8) рассчитывается п- Если найденная величина удовлетворяет условию (11), то реактивность вычисляется по формулам (6), (7) и (15). Если условие (11) нарушается, то влияние точки на расчет производной усиливается последовательным применением формул (9) и (10);

шаг 2 - путем интерполяции с помощью формулы (12) вводятся две имитации измерений в серединах интервалов [г-1, г] и [г, г+1]. Производная вычисляется по формуле (9). Если полученное

Теплоэнергетика

значение удовлетворяет условию (11), то разрешается переход к счету реактивности по формулам (6), (7) и (15). В противном случае следует шаг 3 - интервалы измерения [z-1, Z] и [/, /+1] делятся на три равные части. По формулам (13) имитируются три замера плотности нейтронов, необходимые для расчета производной по формуле (10). Если значение производной удовлетворяет условию (11), то выполняется переход к формулам (6), (7) и (15). При невыполнении условия (11) следует шаг 4 - производная п- рассчитывается по формуле

и;=(и,+1+5,+1-и,)/т, реализующей наибольшее допустимое сглаживание. Знак перед 5,+] выбирается из условия соответствия наименьшему абсолютному значению производной. Реактивность восстанавливается по формулам (6), (7) и (15)

Для предложенного алгоритма определения реактивности затруднительно провести априорный анализ устойчивости и точности. Исходные данные в виде результатов замера плотности нейтронов задаются приближенно. На каждом шаге измерения вводится лишь одно случайное приближение плотности нейтронов при неопределенности точного значения. Действительная погрешность измерения на данном шаге остается неизвестной. Однако именно эта погрешность наряду с качеством разностной схемы определяет точность и устойчивость счета реактивности. Поэтому предложенный метод алгоритмической регуляризации исследовался на основе модельного эксперимента, позволяющего имитировать различные режимы ядерного реактора, в том числе и аварийные. При этом оказывается возможным и введение погрешности в модель процесса измерения в целях оценки зависимости погрешности определения реактивности от точности измерения плотности нейтронов. На первом этапе модельного эксперимента решается прямая задача для уравнений (1)—(3): по заданной функции к (t) находят п (t). Для решения прямой задачи применяется разностная схема, описанная в работе [4]. На втором этапе в решение прямой задачи - найденную функцию п (t) вносится погрешность заданного уровня

«,(?,) = и(?,)(1+7®)> где у - уровень погрешности; 0 - случайные числа из интеграла [-1; 1], распределенные по нормальному закону На третьем этапе решается обратная задача для уравнений нейтронной кинетики: по заданным значениям и, /, восстанавливается оценка реактивности K,(t).

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

Тестирование, выполненное по описанной схеме, показало, что предложенный метод определения реактивности обладает регуляризирующими свойствами. По каждой паре л,_+], 8,₊₁ находят определенную оценку К_1+ь такую, что К,-, -^ к (й+i) при 8,₊₁ -> 0 В отличие от рекуррентного метода Калмана, основанного на оптимальном оценивании, метод алгоритмической регуляризации позволяет получить оценку, согласованную по точности с погрешностью измерения.

Предложенный метод реализован в устройстве для измерения реактивности ядерного реактива.