INTRODUCCIÓN

En las centrales nucleares se lleva a cabo la producción de energía, la cual se genera a partir de la fisión nuclear. Esta es una reacción muy importante en la física, se origina cuando un neutrón térmico impacta el núcleo de un átomo pesado, siendo el Uranio-235 el material más utilizado, formando un núcleo compuesto en estado excitado Uranio-236, con una energía de activación mayor que la barrera de fisión; del impacto resulta una competencia entre la fuerza de coulomb y la fuerza de superficie descrito en el modelo de la gota liquida con correcciones considerando el modelo de capas, con esta teoría es posible explicar la división del núcleo en dos o más núcleos con mayor probabilidad que la división sea asimétrica, liberando cero o más neutrones, produciendo una reacción en cadena debido a que los neutrones generados pueden impactar a otros núcleos en los elementos de combustible. Este proceso libera una gran cantidad de energía calorífica, mediante transformaciones de energía es distribuida como energía eléctrica (^{Stacey, 2018}). Para controlar con seguridad el reactor nuclear es necesario conocer la reactividad, esto se consigue mediante métodos numéricos, los cuales usados en diferentes áreas como en el comportamiento térmico en régimen transitorio de flujo de gas en redes de ductos, en presencia de varios tipos de componentes como válvulas, compresores, intercambiadores de calor, y bifurcaciones (^{Ortega et al., 2009}), en la propagación de una fisura (^{Palma et al., 2010}), en el proceso de mezclado metano-oxígeno para un sistema confinado con arreglo de chorros 4-Lug-Bolt bajo la influencia de la variación de presión (^{De la Cruz et al., 2015}), para la reactividad nuclear puede ser obtenida discretizando la ecuación de la cinética puntual (^{Shimazu et al., 1987}; ^{Kitano et al., 2000}; ^{Tamura, 2003}), otros métodos discretizan la dependencia integral de la ecuación inversa de la cinética puntual (^{Hoogenboom y Van Der Sluijs,1988}; ^{Binney y Bakir, 1989}; ^{Ansari, 1991}).

Otras investigaciones se han desarrollado con el fin de discretizar la integral reduciendo el costo computacional y aumentando la precisión en el cálculo de la reactividad (^{Malmir y Vosoughi, 2013}) y, un método usando la generalización de los números de Bernoulli (^{Suescún-Díaz et al., 2020}). En este trabajo se discretiza el histórico de la densidad de neutrones empleando la fórmula de Euler-Maclaurin y se hace una aproximación de las derivadas de la respuesta del sistema a un impulso unitario (^{Haykin, 2014}) y una aproximación de los infinitos números de Bernoulli (^{Suescún et al., 2013}). Se simulan las fluctuaciones en la densidad de neutrones suponiendo que sus medidas contienen ruido con una distribución Gaussiana alrededor de un valor medio de la densidad de neutrones y se emplea el filtro paso-bajo de primer orden de atraso (^{Shimazu et al., 1987}) para analizar su comportamiento frente a la reducción de fluctuaciones.

DESARROLLO DE ECUACIONES

Las ecuaciones de la cinética puntual describen la evolución temporal de la distribución de neutrones y la concentración de precursores de neutrones atrasados. Este conjunto de siete ecuaciones fuertemente acopladas modela la dinámica de un reactor nuclear (^{Duderstadt y Hamilton, 1976}):

Las Ecuaciones (1) y (2) están sujetas a las siguientes condiciones iniciales:

siendo P(t) la densidad de neutrones, C _i(t) la concentración del i-ésimo grupo de neutrones atrasados, ρ(t) la reactividad, Λ el tiempo de generación de neutrones instantáneos, β_i la i-ésima fracción de neutrones atrasados, β la fracción total efectiva de neutrones atrasados y λ_i la constante de decaimiento del i-ésimo grupo de neutrones atrasados.

Es posible calcular una expresión para la concentración de precursores al resolver la Ecuación (2) haciendo uso de las condiciones iniciales representadas en las Ecuaciones (3) y (4). Al despejar la reactividad de la Ecuación (1) y remplazar el término que se obtiene para la concentración de precursores, se tiene la ecuación inversa de la cinética puntual:

Reactividad con la aproximación del primer y segundo número de Bernoulli

Para discretizar la dependencia integral de la Ecuación (5) se hace uso de la fórmula de Euler-Maclaurin (^{Kwok., 2010}):

donde los B_p son los números de Bernoulli, F(t) es una función continua con infinitas derivadas continuas. La fórmula de Euler-Maclaurin permite aproximar las integrales de una función continua como una aproximación de la función en valores discretos. Para simplificar el integrando en la Ecuación (6), se hace la siguiente sustitución:

siendo h_i la respuesta del sistema a una función de impulso unitario (^{Haykin, 2014}). Ahora, haciendo el producto entre la respuesta al impulso dado por la Ecuación (7) y la densidad de P(t) se obtiene el integrando de la Ecuación (6):

La versión discreta de la Ecuación (8) será:

Si se trunca la serie hasta p=1 en la Ecuación (6) y al derivar una vez la Ecuación (9) es posible mostrar que la reactividad con el primer número de Bernoulli (Suescún et al., 2020) es:

donde T es el tamaño de paso de tiempo; al truncar la serie en p=2 en la Ecuación (6) y al derivar dos veces la Ecuación (9) es posible mostrar que la reactividad con el segundo número de Bernoulli es la mostrada en la Ecn. (11), siendo ρ_FIR[n] los cuatro primeros términos de la Ecuación (10).

Generalización de los números de Bernoulli

Es posible generalizar el cálculo de la reactividad obteniendo más precisión considerando más números de Bernoulli (p tiende a infinito), sin embargo sería una expresión muy extensa, la idea es simplificar esa expansión encontrando alguna serie que converja, con esto los infinitos números de Bernoulli pueden dar una mayor precisión en el cálculo de la reactividad, para encontrar esa expresión, se hace uso del teorema del binomio, obteniendo derivadas de cualquier orden para la respuesta del sistema a una función de impulso unitario h_i. La forma general de las Ecuaciones (10-11) puede ser escrita,

La Ecuación (12) contiene los términos pares e impares de las derivadas para la respuesta del sistema al impulso unitario h_i(t) y para P(t). De la Ecuación (7) es posible obtener la n-ésima derivada de la respuesta del sistema a un impulso unitario con respecto a t ’ de la siguiente forma:

Remplazando la versión discreta de la Ecuación (13) en la Ecuación (12) se tiene:

Las infinitas derivadas de la densidad de población de neutrones que aparecen en la Ecuación (14), se pueden reducir a solo dos derivadas, haciendo la siguiente suposición (Suescún et al.,2020):

donde a está dado por:

Al remplazar las Ecuaciones (15-17) en la Ecuación (14), después de realizar varias operaciones algebraicas y haciendo Z^2p _i,1=[T(a+λ_i)]^2p y Z^2p _i,2=[T(a-λ_i)]^2p, se puede llegar a una expresión que representa una serie de Laurent de la forma ∑^∞ _{p=1 (-1)} ^p-1(B_pZ^2p)/(2p)!=(z/e^z-1)-1+z/2, es posible obtener la Expresión (18).

La Ecuación (18) representa la reactividad con la generalización de los números de Bernoulli, aparentemente es una expresión complicada, pero depende solo hasta la segunda derivada de la densidad de neutrones. Es una fórmula exacta para las funciones que cumplan el criterio dado por las Ecuaciones (15-17). En la siguiente subsección se propone someter el método presentado en la Ecuación (18) a fluctuaciones las cuales se reducen mediante la aplicación del filtro paso-bajo.

Filtro paso-bajo

Para simular las fluctuaciones, se supone un ruido con una distribución Gaussiana alrededor de un valor medio de la densidad de neutrones, representado de la siguiente forma (^{Kitano et al., 2000}):

donde N es el número de muestras y P_j es un vector de N muestras de la densidad de neutrones para ser filtradas. En este trabajo, para reducir las fluctuaciones en el cálculo de la reactividad, consideramos el filtro paso-bajo de primer orden de atraso (^{Shimazu et al., 1987}), el cual está dado por la siguiente ecuación:

donde τ (tau) es la constante de filtrado, P_i-1 y P_i son las densidades de neutrones filtradas, T es el tamaño de paso como está definido anteriormente.

En la siguiente sección se presentan los resultados obtenidos al someter el método de la generalización de los números de Bernoulli, dado en la Ecuación (18), el ruido con distribución Gaussiana alrededor de un valor medio dado por la Ecuación (19), y el filtro paso-bajo el cual está dado por la Ecuación (20). El porcentaje de reducción se calcula mediante la Ecuación (21).

RESULTADOS

En esta sección se presentan algunos de los resultados más sobresalientes obtenidos para diferentes experimentos numéricos, usando el método propuesto dado por la ecuación (20). Los experimentos numéricos se realizan para la densidad de neutrones de la forma P(t)=exp(ωt) que se cumple en la aproximación dada por las ecuaciones (15-17) con diferentes valores de ω(s^-1), y para la forma P(t)=a+bt³ que no cumple el criterio de las ecuaciones (15-17) donde se fija el valor de a y se varía el valor de b, estas simulaciones numéricas se realizaron para diferentes intervalos en el tiempo de simulación denotado por el tiempo final t_f. La curva de reactividad se comparó con el método de referencia, el cual se obtiene al resolver analíticamente la ecuación (5); el método de referencia permite analizar la precisión del método propuesto ya que es posible obtener una solución exacta de la ecuación de la reactividad presentada en la ecuación (5).

Las constantes físicas empleadas en este trabajo obedecen a la interacción de los neutrones con el material combustible ²³⁵U; estas constantes son la constante de decaimiento λ_i = {0.0127, 0.0317, 0.115, 0.311, 1.4, 3.87} s^-1, la fracción de neutrones atrasados β_i = {0.000266, 0.001491, 0.001316, 0.002849, 0.000896, 0.000182} y el tiempo de generación de neutrones instantáneos Λ=2x10⁵ s. Para simular el ruido, se supone que varía con una distribución Gaussiana alrededor de un valor medio de la densidad de la población de neutrones, para esta suposición se utiliza una semilla generadora de números aleatorios de la forma 2³¹-1 y una desviación estándar que se varía entre σ=0.001 y σ=0.01. Para reducir las fluctuaciones, se emplea el filtro paso-bajo de primer orden de atraso dado por la ecuación (20) con una constante de filtrado τ=0.1 s y τ=2.5 s con un tamaño de paso fijo de T=0.1 s para el cálculo de la reactividad en todos los casos.

Los diferentes resultados numéricos obtenidos con el método propuesto usando la generalización de los números de Bernoulli, denotado como GB, son comparados con el método de Euler Maclaurin con la aproximación del primer número de Bernoulli denotado como EM₁ y con la aproximación del segundo número de Bernoulli denotado EM₂. En las Tablas 1 - 4 se muestran los valores de las máximas diferencias en la reactividad en pcm (partes por cien mil) y del error medio absoluto obtenidos con el método propuesto dado en la ecuación (18) con un ruido de la forma dada en la ecuación (19) y filtrando de la forma dada en la ecuación (20) conocido en este trabajo como filtro paso-bajo de primer orden de atraso (FPB) con una desviación estándar σ=0.001. En la Tabla 1 se observa que es posible una reducción en la fluctuación de la reactividad calculada usando la ecuación (21), considerando la máxima diferencia para un valor muy bajo de ω=0.00243 (s^-1) alcanzando un valor de hasta 7.81%. En la Tabla 2 muestra que se pueden alcanzar reducciones de 81.13%, 35.29% y 34.69% para valores de ω=0.00243 (s^-1), ω=0.006881 (s^-1), y ω=0.01046 (s^-1), respectivamente, con los tres métodos EM₁+FPB, EM₂+FPB y GB+FPB cuando se considera el error medio absoluto. Así mismo, en las Tablas 3-4 se observan reducciones, en el caso cuando se considera la máxima diferencia, se obtienen reducciones de 77.30%, 75.24%, 76.59% y 46.72% y de 75.23%, 71.67%, 63.63% y 37.97% cuando se considera el error medio absoluto, en ambos casos se observa las reducciones para los valores ω≤0.02817(s^-1).

Las Tablas 5 - 6 muestran los resultados obtenidos para una densidad de neutrones de la forma P(t)=a+bt³ donde se ha dejado fijo el valor de a en la unidad, se varía el valor de b (^{Suescún et al., 2007}) y el tiempo de cálculo se hace hasta t_f = 3950 s, con una desviación standard σ=0.01 y una constante de filtrado τ = 0.1 s el método de la generalización de los números de Bernoulli y el filtro paso-bajo de primer orden de atraso logra reducir más las fluctuaciones para un valor de b=(0.0127)⁵/9.

La Figura 1 presenta la curva de reactividad para una densidad de neutrones de forma cúbica con a=1 y b=(0.0127)⁵/9, para un tiempo total de simulación de hasta t_f = 10000 s. Es posible apreciar que el método propuesto usando la Ecuación (18) y el filtro paso-bajo dado por la Ecuación (20), reduce las fluctuaciones significativamente acercándose más al método de referencia el cual se obtiene al resolver analíticamente la Ecuación (5), lo que no se logra con los métodos usando el primer número de Bernoulli con el filtro paso-bajo EM₁+FPB y con el segundo número de Bernoulli con el filtro paso-bajo EM₂+FPB que son prácticamente iguales, llegando a tener una máxima diferencia de aproximadamente 25 pcm, se necesitan más números de Bernoulli para aumentar principalmente la precisión y reducir las fluctuaciones como el indicado usando la convergencia de los infinitos números de Bernoulli representado en la Ecuación (18), aunque está ecuación requiere cumplir la condición de convergencia establecida principalmente en la Ecuación (17) presenta mejores resultados en la reducción de las fluctuaciones que los métodos con uno o dos números de Bernoulli.

DISCUSIÓN FINAL

Se presentó un método para el cálculo de la reactividad que reduce las fluctuaciones, empleando los diferentes números de Bernoulli con su generalización y el filtro paso-bajo de primer orden de atraso. El ruido se simuló alrededor de un valor medio de la densidad de neutrones con una distribución Gaussiana y una semilla generadora de números aleatorios de 2³¹-1. Los diferentes experimentos numéricos muestran que, para una selección de la forma de la densidad de población de neutrones exponencial, es posible la reducción de las fluctuaciones representadas por las máximas diferencias y el error medio absoluto, mejores resultados son obtenidos para valores de ω=0.00243 (s^-1) obteniendo reducciones de hasta 81.13%. Cuando se cambia la forma de la densidad de la población de neutrones a una de forma cúbica, es posible la reducción de las fluctuaciones bajo ciertas condiciones reportadas en este trabajo, con una constante de filtrado de τ=0.1 s obteniendo reducciones de hasta el 79.16% con el método propuesto. Cuando el valor de σ aumenta es posible observar que la reducción en las fluctuaciones resulta ser apropiada cuando se usa el método propuesto. En general, las diferencias principalmente presentadas se deben al filtro de repuesta al impulso finito FIR, el cual es el primer término y el dominante en el cálculo de la reactividad.

CONCLUSIONES

De acuerdo con el trabajo presentado y los resultados obtenidos, se pueden plantear las siguientes conclusiones principales: 1) En el cálculo de la reactividad se ve que es posible reducir las fluctuaciones usando el método propuesto con la convergencia de los infinitos números de Bernoulli representada por medio de la ecuación (18) y el filtro paso-bajo de primer orden de atraso representado en la ecuación (20) con diferentes constantes de filtrado, suponiendo que la población de neutrones tiene una fluctuación con una distribución Gaussiana alrededor de su valor medio indicado en la ecuación (19); (2) En la reducción de fluctuaciones se ve que se alcanzan valores hasta un 80% para el caso que cumplan o no el criterio de la reducción de la derivadas pares e impares dada por las ecuaciones (15-17) usando la generalización de los números de Bernoulli; y (3) Para formas exponenciales de la densidad de neutrones es suficiente con usar solamente el primer número de Bernoulli, sin embargo, para la forma cubica que no cumplen las propiedades de la convergencia dada en las ecuaciones (15-17), se ve que el primer y segundo número de Bernoulli presenta dificultades al inicio de la simulación numérica siendo mejor el desempeño cuando se usa la generalización de los números de Bernoulli.