Application of the generalized coordinates method to the study of three-dimensional dynamics of nuclear reactors
Abstract Generated abstract
The paper applies the method of generalized coordinates to the three-dimensional space-time dynamics of a cylindrical nuclear reactor in a one-group diffusion approximation. Neutron flux deviations are expanded in eigenfunctions of the spatial operator, while fuel and coolant temperature deviations are represented by Fourier series, reducing the coupled partial differential equations with control-rod perturbations to a system of ordinary differential equations for modal amplitudes. The resulting model was implemented on an EMU-10 analog computer using a truncated set of twelve neutron-flux modes and four temperature modes, which the authors report gave the required accuracy. The approach is presented as a basis for constructing a global reactor model, with possible refinements including xenon and iodine poisoning, multigroup neutron treatment, and more complete delayed-neutron kinetics.
Full Text
PROCEEDINGS
OF THE TOMSK ORDER OF THE RED BANNER OF LABOR POLYTECHNIC
INSTITUTE named after S. M. KIROV
Vol. 168
1969
APPLICATION OF THE METHOD OF GENERALIZED COORDINATES TO THE STUDY OF THREE-DIMENSIONAL DYNAMICS OF NUCLEAR REACTORS
V. A. TARTAKOVSKII, I. G. VINTIZENKO
(Presented by the scientific seminar of the computing laboratory of TPI)
This paper considers the application of expansion in a series in eigenfunctions for solving, on an analog computer, partial differential equations describing the space-time dynamics of a nuclear reactor.
The essence of the method considered is that the sought function is represented by an infinite series in the eigenfunctions of the spatial operator of the original problem. The time dependence is contained in the coefficients of the series—the modes.
1. Let us consider a nuclear reactor that is a cylinder with extrapolated height \(H\) and radius \(R\). This cylinder contains fuel, moderator, and control rods. Cooling, sufficient to compensate for radial and azimuthal nonuniformities of the temperature field, is carried out by gas circulating from bottom to top through channels parallel to the vertical axis (the \(z\)-axis) of the cylinder. The influence of the moderator on thermal processes in the system is assumed to be negligible. We neglect heat diffusion in the fuel and coolant.
The steady state of such a reactor in the one-group diffusion approximation is described by the equations
\[ M^2 \nabla^2 \Phi_c(r,\theta,z) + B_c^2 \Phi_c(r,\theta,z) = 0; \tag{1} \]
\[ \mu U \overline{\Phi}_c(z) - \gamma_{\mathrm{ig}} F_{\mathrm{i}}\,[T_{\mathrm{i}}(z)-T_{\mathrm{g}}(z)] = 0; \tag{2} \]
\[ C_{\mathrm{g}} M_{\mathrm{g}} v_{\mathrm{g}} \frac{dT_{\mathrm{g}}(z)}{dz} - \gamma_{\mathrm{ig}} F_{\mathrm{i}}\,[T_{\mathrm{i}}(z)-T_{\mathrm{g}}(z)] = 0, \tag{3} \]
where
\[ \overline{\Phi}_c(z) = \frac{1}{\pi R^2} \int_{0}^{2\pi}\int_{0}^{R} \Phi_c(r,\theta,z)\, r\,dr\,d\theta, \]
and the boundary conditions for \(\Phi_c\) have the form
\[ \Phi_c(R,\theta,z)=\Phi_c(r,\theta,H)=\Phi_c(r,\theta,0)=0. \]
It is known that equation (1) is an eigenvalue equation. In general form it is written as
\[ \nabla^2 f_{lmn} + L_{lmn}^{2} f_{lmn} = 0. \]
For the cylinder under consideration and with zero boundary conditions, representing \(f_{lmn}\) in the form \(f_{lmn}=Z(z)R(r)S(\theta)\), we find solutions that are finite and continuous in \(U\):
\[ Z(z)=\sin \lambda_n z;\qquad \lambda_n=\frac{\pi(n+1)}{H};\qquad n=0,1,2,3\ldots \]
\[ S_l(\theta)=\sin m\theta;\qquad S_n(\theta)=\cos m\theta;\qquad m=0,1,2\ldots \]
\[ R(r)=J_m\left(\frac{\gamma_{lm}}{R}r\right);\qquad J_m(\gamma_{lm})\equiv 0; \]
\[ \|f_{lmn}\|^2=\frac{\pi HR^2}{2\varepsilon_{lm}^2};\qquad m\ne 0 \]
\[ \|f_{lon}\|^2=\frac{\pi HR^2}{\varepsilon_{lo}^2};\qquad L_{lmn}^2=\frac{\gamma_{lm}^2}{R^2}+\lambda_n^2. \]
If
\(B_c^2=\mathrm{const}\), then, since \(\Phi_c(r,\theta,z)>0\),
\[ \Phi_c=P\sin\frac{\pi}{H}z\cdot J_0\left(\frac{\gamma_{00}}{R}r\right); \tag{4} \]
\[ B_c^2=M^2\left(\frac{\gamma_{00}^2}{R^2}+\frac{\pi^2}{H^2}\right). \]
Here \(P\) is a constant characterizing the power level. Next we find that
\[ T_r=-\frac{\kappa\mu UHP}{\pi C_r M_r v_r}\cos\frac{\pi}{H}z; \tag{5} \]
\[ T_u=\kappa\mu UP\left[ \frac{\sin \frac{\pi}{H}z}{\gamma_{ur}F_u} -\frac{H\cos\frac{\pi}{H}z}{C_r M_r v_r\pi} \right]. \tag{6} \]
The reactor is controlled by five rods located, respectively, at the points
\[ (r',0,z);\quad \left(r',\frac{\pi}{2},z\right);\quad (r',\pi,z);\quad \left(r',\frac{3\pi}{2},z\right);\quad (0,\theta,z). \]
In the stationary state they are at half the height of the reactor. The thickness of the rods is incomparably smaller than the cross section of the reactor; therefore we write the equation for a rod in the following form:
\[ \rho_i=\delta[(r,r_i)(\theta,\theta_i)]\times\left(z-\frac{H}{2}\right). \tag{7} \]
- Let us write the equations for a nonstationary reactor (for deviations from the stationary level).
\[ M^2\nabla^2\varphi(r,\theta,z,t)+B_c^2\varphi(r,\theta,z,t) +\left[\alpha_u\tau_u(z,t)+\alpha_r\tau_r(z,t)+\sum_{i=1}^{5}a_i\rho_i\right]\times \]
\[ \times \Phi_c(r,z)=l^*\frac{\partial\varphi}{\partial t}; \tag{8} \]
\[ C_uM_u\frac{d\tau_u(z,t)}{dt} =\mu U\bar{\varphi}(z,t)-\gamma_{ur}F_u[\tau_u(z,t)-\tau_r(z,t)]; \tag{9} \]
\[ C_{\mathrm{r}}M_{\mathrm{r}}\frac{\partial \tau_{\mathrm{r}}(z,t)}{\partial t} = \gamma_{\mathrm{wr}}F_{\mathrm{u}}\,[\tau_{\mathrm{u}}(z,t)-\tau_{\mathrm{r}}(z,t)] - C_{\mathrm{r}}M_{\mathrm{r}}v_{\mathrm{r}}\frac{\partial \tau_{\mathrm{r}}}{\partial z}(z,t). \tag{10} \]
We represent \(\varphi(r,\theta,z,t)\) by a series in eigenfunctions
\[ \varphi = P \sum_{l,m,n=0}^{\infty} J_m\!\left(\frac{\gamma_{lm}}{R}r\right) \sin \lambda_n z \,[\varphi_{lmn}(t)\cos m\theta+\psi_{lmn}(t)\sin m\theta]. \]
Substitute this expression into equation (8) and take (4), (7) into account. Then multiply what is obtained by some eigenfunction and integrate over the region \(U\); since \(f_{lmn}\) are orthogonal, we obtain:
\[ \dot{\varphi}_{lmn} = \omega_{lmn}^{2}\varphi_{lmn} + \|f_{lmn}\|^{-2} \left\{ \int_U [\alpha_{\mathrm{u}}\tau_{\mathrm{u}}+\alpha_{\mathrm{r}}\tau_{\mathrm{r}}]\, J_0\!\left(\frac{\gamma_{00}}{R}r\right) J_m\!\left(\frac{\gamma_{lm}}{R}\right) \right. \]
\[ \left. {}\times \sin\lambda_0 z\cdot \sin\lambda_n z\cos m\theta\,dv + \sum_{i=1}^{5} a_i r' J_0\!\left(\frac{\gamma_{00}}{R}r'\right) J_m\!\left(\frac{\gamma_{lm}}{R}r'\right) \cos m\theta_i\,h_n^i[z(t)] \right\}. \tag{11} \]
\[ \dot{\psi}_{lmn} = \omega_{lmn}^{2}\psi_{lmn} + \|f_{lmn}\|^{-2} \sum_{i=1}^{5} a_i r' J_0\!\left(\frac{\gamma_{00}}{R}r'\right) J_m\!\left(\frac{\gamma_{lm}}{R}r'\right) \]
\[ {}\times \sin m\theta_i\,h_n^i[z(t)], \tag{12} \]
where
\[ \omega_{lmn}^{2}=M^2\bigl(L_{000}^{2}-L_{lmn}^{2}\bigr), \]
and
\[ h_0^i = \frac{H}{4\pi} \left[ \pi+\sin\frac{2\pi}{H}z(t)-\frac{2\pi}{H}z(t) \right], \]
\[ h_1^i = \frac{H}{6\pi} \left[ 3\sin\frac{\pi}{H}z(t)-\sin\frac{3\pi}{H}z(t)-4 \right]. \]
- From (11), (12) it is evident that the coupling of the neutron flux with the temperature is effected through the integral
\[ I_n = \int_U J_0\!\left(\frac{\gamma_{00}}{R}r\right) J_m\!\left(\frac{\gamma_{lm}}{R}r\right) \sin\lambda_n z\sin\lambda_0 z\cos m\theta \times [\alpha_{\mathrm{u}}\tau_{\mathrm{u}}+\alpha_{\mathrm{r}}\tau_{\mathrm{r}}]\,dv. \]
This expression is nonzero only for \(l=0\) and \(m=0\). Taking this into account, we obtain
\[ I_n = \frac{2\pi R^2}{\varepsilon_{00}^{2}} \int_0^H \sin\lambda_0 z\cdot \sin\lambda_n z [\alpha_{\mathrm{u}}\tau_{\mathrm{u}}+\alpha_{\mathrm{r}}\tau_{\mathrm{r}}]\,dv. \]
Now our task is to eliminate the variable \(z\) from this expression. There are several ways to do this. We shall consider representing \(\tau_{\mathrm{u}}\) and \(\tau_{\mathrm{r}}\) on the interval \(0-H\) by a cosine Fourier series; we shall have
\[ I_n = \frac{R^2H^2}{2\varepsilon_{00}^{2}} \sum_{\kappa=1}^{2} \alpha_{\kappa} \sum_{\nu=0}^{\infty} \tau_{\kappa}^{\nu} \left[ \frac{\sin(n-\nu)\pi}{n-\nu} + \frac{\sin(n+\nu)\pi}{n+\nu} - \frac{\sin(n-\nu-2)\pi}{n-\nu+2} \right]. \]
Substitute the expressions
\[ \tau_{\kappa}(z,t) = \sum_{\nu=0}^{\infty} \tau_{\kappa}^{\nu}(t)\cos\frac{\pi\nu}{H}z; \]
\[ \bar{\varphi}(z,t)=\frac{1}{\pi R^2}\sum_{l,n=0}^{\infty}\varphi_{10n}(t)\sin\lambda_n z \int_0^{2\pi}\int_0^R rJ_0\left(\frac{i_{10}}{R}r\right)\,dr\,d\theta = \]
\[ =\frac{2}{R^2}\sum_{l,n=0}^{\infty} b_l\varphi_{10n}(t)\sin\lambda_n z \]
In equations (9), multiply them by
\[ \left\|\cos\frac{\pi\nu}{H}z\right\|^{-2}\cos\frac{\pi\nu}{H}z \]
and integrate from 0 to \(H\); we obtain:
for \(\nu=0\)
\[ \dot{\tau}_{\mathrm{n}}^0= \frac{\mu U H}{\pi^2 R^2 C_{\mathrm{n}}M_{\mathrm{n}}} \sum_{l,n=0}^{\infty} b_l\varphi_{10n} \frac{G(n)}{n+1} -\frac{\gamma_{\mathrm{ng}}F_{\mathrm{n}}}{C_{\mathrm{n}}M_{\mathrm{n}}} (\tau_{\mathrm{n}}^0-\tau_{\mathrm{g}}^0); \]
\[ \dot{\tau}_{\mathrm{g}}^0= \frac{v_{\mathrm{g}}}{\pi}\sum_{j=0}^{\infty}\tau_{\mathrm{g}}^jY(j) +\frac{\gamma_{\mathrm{ng}}F_{\mathrm{n}}}{C_{\mathrm{g}}M_{\mathrm{g}}} (\tau_{\mathrm{n}}^0-\tau_{\mathrm{g}}^0); \]
for \(\nu\ne0\)
\[ \dot{\tau}_{\mathrm{n}}^\nu= \frac{2\mu U H}{\pi^2 R^2 C_{\mathrm{n}}M_{\mathrm{n}}} \sum_{l,n=0}^{\infty} b_l\varphi_{10n} \left[ \frac{G(n+\nu)}{n+\nu+1} + \frac{G|n-\nu|}{n-\nu+1} \right] -\frac{\gamma_{\mathrm{ng}}F_{\mathrm{n}}}{C_{\mathrm{n}}M_{\mathrm{n}}} (\tau_{\mathrm{n}}^\nu-\tau_{\mathrm{g}}^\nu); \]
\[ \dot{\tau}_{\mathrm{g}}^\nu= \frac{2v_{\mathrm{g}}}{\pi} \sum_{\substack{j=0\\j\ne\nu}}^{\infty} j\tau_{\mathrm{g}}^j \left[ \frac{Y(\nu+j)}{\nu+j} + \frac{Y(j-\nu)}{j-\nu} \right] + \frac{\gamma_{\mathrm{ng}}F_{\mathrm{n}}}{C_{\mathrm{g}}M_{\mathrm{g}}} (\tau_{\mathrm{n}}^\nu-\tau_{\mathrm{g}}^\nu), \]
where
\[ G(n)= \begin{cases} 1, & \text{for even } n,\\ 0, & \text{for odd } n, \end{cases} \]
\[ Y(n)=1-G(n). \]
The equations written above were solved on the “EMU-10” computer. Naturally, only the first terms of the infinite series were taken into account. Twelve “modes” were taken for the neutron flux and four for the temperature. This approximation provided the required accuracy.
It should be noted that the analog scheme obtained was simple; therefore, further refinement of the model is possible. In this respect a nuclear reactor offers considerable possibilities. One way of refinement is to take into account the effect of poisoning by xenon and iodine. The difficulties that arise in this connection are associated with the calculation of rather complicated definite integrals. Next, the distribution of neutrons over energies should be taken into account (at least in a two-group approximation), and the kinetics with delayed neutrons should be more fully allowed for. However, it must be remembered that the discussion concerns small oscillations of the investigated parameters about a stationary level; moreover, these refinements are feasible only if a sufficient amount of precision equipment is available.
Taking into account everything stated above, the method considered makes it possible to construct a global model of a nuclear reactor.
REFERENCES
-
L. V. Kantorovich, V. I. Krylov. Approximate Methods of Higher Analysis. GITTL, Moscow, 1949.
-
D. Jackson. Fourier Series and Orthogonal Polynomials. GIIL, Moscow, 1948.
-
Gian Paolo Caliguri, Herve d’Hoor. Nuovo metodo sperimentale analogico per la simulazione tridimensionale dei reattori nucleari. Rivista di ingegneria. Vol. 12, No. 5, 493–507, 1962.