A Model Study of a Homogeneous Light-Water Thorium Reactor

Abstract

This work presents a computational study of a ${}^{232}\mathrm{Th}$-based homogeneous light-water reactor. Thorium reactors have been proposed as an alternative to the uranium fuel cycle since they exploit both the availability of thorium and its ability to afford fissile uranium isotopes by a sequence of neutron captures. Besides ${}^{233}\mathrm{U}$, as a result of the neutron captures, a significant amount of ${}^{234}\mathrm{U}$ (36.3%) and 6.46% of ${}^{235}\mathrm{U}$ are formed in the reactor under study. More importantly, the proposed simulation points out the possibility of a continuous withdrawal of the uranium isotopes without compromising the criticality and the power output of the reactor. This withdrawal affords the fissile material for the startup of reactors other than the first one, which requires a one-time only limited amount of fissile material. The significant molar fraction of the ${}^{234}\mathrm{U}$ (0.17) in the extracted fuel does not pose a limitation on weapon proliferation, as a consequence of its high fission cross section for high-energy neutrons.

1. Introduction

Thorium reactors have been studied as an alternative to uranium-based reactors mainly because thorium is more widely available than uranium [1,2]. In a thorium reactor, the fissile material ${}^{233}\mathrm{U}$ is obtained by the non-fissile ${}^{232}\mathrm{Th}$ by the sequence of neutron (n) capture and ${\beta }^{-}$ decays

$${}^{232}\mathrm{Th}\stackrel{n}{\to }{}^{233}\mathrm{Th}\stackrel{{\beta }^{-},21.83\min }{\to }{}^{233}\mathrm{Pa}\stackrel{{\beta }^{-},26.975\mathrm{d}}{\to }{}^{233}\mathrm{U}$$

(1)

Besides the availability of fuel, other advantages have been suggested, such as the reduced production of heavy actinides, in particular ${}^{239}\mathrm{Pu}$ that carries the possibility of nuclear weapon proliferation. The fuel cycle of a thorium based reactor hinges on the fission of ${}^{233}\mathrm{U}$ and does not contain any ${}^{238}\mathrm{U}$, the precursor of ${}^{239}\mathrm{Pu}$. Furthermore, the occurrence of ${}^{232}\mathrm{U}$ by the reaction $\left(n,2n\right)$ on ${}^{233}\mathrm{U}$ is thought to impede the usage of ${}^{233}\mathrm{U}$ for nuclear weapons [3] because of the relatively short half-life time ($t_{1/2}$) of ${}^{232}\mathrm{U}$ (68.9 yr) that entails the presence of its hard-gamma emitting daughter ${}^{208}\mathrm{Tl}$ ($t_{1/2}=3.053\min$). It is not the purpose of this work to discuss the extensive literature on this topic; interested readers may refer to the excellent review by Ault [4]. Most existing suggestions for the configuration of thorium reactors use molten salts as medium and graphite as a moderator [5,6]. Furthermore, fast reactors simulations have been suggested, either based on molten salts [7,8], lead cooled [9], or water. Among the latter we have a boiling water reactor (BWR) study [10], and pressurized water reactor (PWR) studies based on mixed uranium–plutonium [11] or uranium–thorium [12] fuels. A reduction in power density leading to a reduced temperature of the fuel elements was predicted for a mixed oxide ${\mathrm{UO}}_2$-${\mathrm{ThO}}_2$ configuration compared to pure ${\mathrm{UO}}_2$ [13]. Furthermore, a reduction of burnable neutron poisons was pointed out for Small Modular Reactors operating on ${\mathrm{UO}}_2$-${\mathrm{ThO}}_2$ [14]. However, to the author’s knowledge, no study has been published concerning a light-water PWR based solely on thorium (except for the startup phase, requiring a fissile component). The fissile material is being introduced into the reactor under study only on the first two days, i.e., the startup phase. After that, the reactor produces enough fissile material to support the production of energy. More importantly, the proposed simulation points out the possibility of a continuous withdrawal of the uranium isotopes without compromising neither the criticality nor the power output of the reactor. Thus, the proposed reactor could be exploited as a uranium source for the existing uranium-based reactors in the case the uranium supply would become problematic.

After the few thorium based reactors that have been operating in the past were discontinued, a few programs exist, like the 300 MW Advanced Heavy Water Reactor planned in India, making use of an external feed of plutonium. The THOREX (Thorium–Uranium${}^{233}$ Extraction) reprocessing cycle has been discussed in [15].

This work develops a model study for a homogeneous thorium reactor in light water with an external feed that needs to be operating only for the first two days after startup. The requirement for the startup of a thorium reactor implies the use of fissile material [16], since ${}^{232}\mathrm{Th}$ is not itself fissile. To this purpose, both ${}^{235}\mathrm{U}$ and ${}^{239}\mathrm{Pu}$ have been used. Instead, this work explores the possibility of using the ${}^{233}\mathrm{U}$ bred from a thorium reactor operating with continuous reprocessing that would also avoid the necessity of ${}^{233}\mathrm{Pa}$ extraction [17].

Since the homogeneous reactor would breed a mixture of ${}^{233}\mathrm{U}$, ${}^{234}\mathrm{U}$, and ${}^{235}\mathrm{U}$, it would be possible in principle to extract the bred uranium isotopes to start up other reactors. Although the abundance and role of the bred uranium isotopes is prominent in thorium reactors, it has non been the subject of many investigations so far. This possibility is being investigated in this paper along with the feasibility of maintaining criticality in the presence of a significant sink for all uranium isotopes. The consequence of the composition of the withdrawn uranium on weapon proliferation is briefly examined. The collected fuel could thus be used to start a new thorium reactor without resorting to uranium enriched in the isotope ${}^{235}\mathrm{U}$. In particular, we also want to assess whether the amount of non-fissile ${}^{234}\mathrm{U}$ in the bred fuel would impair criticality during the startup phase. The reactor would operate at a nearly stationary state with a steady input of ${}^{232}\mathrm{Th}$.

In the sequence of neutron captures, beta decays, and fissions taking place in the environment of a reactor core, there are many nuclides representing branching points, because alternative processes to neutron capture, such as the $\left(n,2n\right)$ reaction, are taking place. In this study, we only take into consideration the main sequence from ${}^{232}\mathrm{Th}$ to ${}^{235}\mathrm{U}$. We include the production of ${}^{235}\mathrm{U}$ because we deem it important to quantify its formation in view of the eventuality of nuclear proliferation applications.

2. Definitions and Methodology

2.1. Modeling Nuclear Processes in the Reactor

In order to determine the feasibility of a continuous withdrawal of uranium isotopes from a thorium-based reactor, it is necessary to have a detailed time evolution of the densities of all nuclides. The model has to include both sources and sinks of the uranium isotopes, along with the necessary control procedures for the neutron flux, aimed to keep the reactor critical and affording the desired power output. After establishing the form of all the rate coefficients, the equations for the time evolution can be formulated in terms of all the relevant densities. In particular, the neutron density (or equivalently the flux) needs to be evaluated by detailing all its source, sink, and transport terms. The neutron diffusion term $D\Delta n$ (D is the diffusion coefficient and $\Delta$ the Laplacian operator), describing the changes in neutron density as a function of time and space, can be simplified by approximating the neutron depletion in the reactor volume as homogeneous, and due to the neutron flux across the surface. In this way, the time derivative of the neutron density, determined by the neutron flux at the reactor boundary, can in turn be evaluated by solving the spatial eigenvalue equation for neutron density along with the Neumann boundary condition. This procedure allows the problem to be formulated as a system of ordinary differential equations.

Considering a system of nuclei and neutrons (with energy $\epsilon$, velocity $v_n$, in the center of mass reference frame, and volume V) with reduced mass m and giving a reaction with energy-dependent cross section $\sigma$, the corresponding rate coefficient $\eta$ for the process involving thermal neutrons would have the form $\eta =〈\sigma v_n〉$, that is

$$\begin{array}{ccc}\hfill \eta & =& {\left(\frac{\mathrm{\Lambda }}{h}\right)}^3\frac{1}{V}{\int }_0^{\infty }\sigma v_n\rho \left(\epsilon \right)e^{-\epsilon /kT}\mathrm{d}\epsilon \hfill \\ & =& {\left(\frac{8kT}{\pi m}\right)}^{1/2}{\int }_0^{\infty }\sigma \widehat{\epsilon }{e}^{-\widehat{\epsilon }}\mathrm{d}\widehat{\epsilon },\hfill \end{array}$$

(2)

where the second equality is obtained by defining $\widehat{\epsilon }=\epsilon /kT$ (k being the Boltzmann constant and T the temperature). The quantities $\rho \left(\epsilon \right)$ and $\mathrm{\Lambda }$ represent the density of translational states and the DeBroglie wavelength, respectively,

$$\rho \left(\epsilon \right)=2\pi {\left(\frac{2m}{h^2}\right)}^{3/2}V{\epsilon }^{1/2},\mathrm{\Lambda }=\frac{h}{{\left(2\pi mkT\right)}^{1/2}},$$

(3)

(h is the Planck constant). We take into consideration the following reactions and decays,

$$\begin{array}{c}\hfill {}^{232}\mathrm{Th}+n\stackrel{{\eta }_{c02}}{\to }{}^{233}\mathrm{Th}+\gamma ,\end{array}$$

(4)

$$\begin{array}{c}\hfill {}^{233}\mathrm{Th}\stackrel{{\lambda }_{03}}{\to }{}^{233}\mathrm{Pa}+e^{-}+\widehat{{\nu }_e},\end{array}$$

(5)

$$\begin{array}{c}\hfill {}^{233}\mathrm{Pa}\stackrel{{\lambda }_{13}}{\to }{}^{233}\mathrm{U}+e^{-}+\widehat{{\nu }_e},\end{array}$$

(6)

$$\begin{array}{c}\hfill {}^{233}\mathrm{Pa}+n\stackrel{{\eta }_{c13}}{\to }{}^{234}\mathrm{Pa}+\gamma ,\end{array}$$

(7)

$$\begin{array}{c}\hfill {}^{234}\mathrm{Pa}\stackrel{{\lambda }_{14}}{\to }{}^{234}\mathrm{U}+e^{-}+\widehat{{\nu }_e},\end{array}$$

(8)

$$\begin{array}{c}\hfill {}^{233}\mathrm{U}+n\stackrel{{\eta }_{f23}}{\to }\mathrm{FF}+{\chi }_{1}n,\end{array}$$

(9)

$$\begin{array}{c}\hfill {}^{233}\mathrm{U}+n\stackrel{{\eta }_{c23}}{\to }{}^{234}\mathrm{U}+\gamma ,\end{array}$$

(10)

$$\begin{array}{c}\hfill {}^{234}\mathrm{U}+n\stackrel{{\eta }_{f24}}{\to }\mathrm{FF}+{\chi }_{1}n,\end{array}$$

(11)

$$\begin{array}{c}\hfill {}^{234}\mathrm{U}+n\stackrel{{\eta }_{c24}}{\to }{}^{235}\mathrm{U}+\gamma ,\end{array}$$

(12)

$$\begin{array}{c}\hfill {}^{235}\mathrm{U}+n\stackrel{{\eta }_{f25}}{\to }\mathrm{FF}+{\chi }_{1}n,\end{array}$$

(13)

$$\begin{array}{c}\hfill {}^{113}\mathrm{Cd}+n\stackrel{{\eta }_x}{\to }{}^{114}\mathrm{Cd}+\gamma ,\end{array}$$

(14)

$$\begin{array}{c}\hfill {}^1\mathrm{H}+n\stackrel{{\eta }_{cH}}{\to }{}^2\mathrm{H}+\gamma ,\end{array}$$

(15)

$$\begin{array}{c}\hfill \mathrm{FF}\stackrel{{\lambda }_{dn}}{\to }\mathrm{F}+{\chi }_{2}n,\end{array}$$

(16)

plus the elastic scattering of thermal neutrons by ${}^{232}\mathrm{Th}$ and ${}^1\mathrm{H}$. All rate coefficients for ${\beta }^{-}$ decay are indicated with the letter $\lambda$ with either an ordering subscript or the subscript “$dn$” for the emission of delayed neutrons. The rate coefficients for neutron-nuclei interactions are indicated by $\eta$, with the subscript f for “fission”, c for “capture” by heavy nuclei, and x for the neutron capture by either the fission products or the control neutron absorber. We use the reference density of metallic thorium ${\rho }_0=11.7\mathrm{g}{\mathrm{cm}}^{-3}\equiv 3.04\times {10}^{28}{\mathrm{m}}^{-3}$ to define the reduced densities as $a_{ij}={\rho }_{ij}/{\rho }_0$, ${\rho }_{ij}$ being the number density. We follow the convention used in the Manhattan project, indicating the reduced density with the first lowercase letter of the element symbol, its subscript being the last digit of the atomic number (i) and the last digit of the mass number (j) [18]. For example, the reduced density of ${}_{90}^{232}\mathrm{Th}$, ${}_{91}^{233}\mathrm{Pa}$, and ${}_{92}^{235}\mathrm{U}$ are noted as $t_{02}$, $p_{13}$, and $u_{25}$, respectively. The same indexes label the various sources and sinks s. We use for the reduced density of neutrons the same symbol used for neutrons, n, and the fission fragments in Equations (4)–(16) are indicated as FF. Since the rate coefficients in processes (4)–(16) can be calculated from known data, we may write a set of differential equations for the time evolution of the species, once the initial conditions for the number density for all species are specified. In the non-dimensional time unit ${\lambda }_{13}t=2.9741\times {10}^{-7}{\mathrm{s}}^{-1}t$, the evolution of the homogeneous system for reactions (4)–(16) is given by the equations

$$\begin{array}{ccc}\hfill {\dot{t}}_{02}& =& -{\omega }_{c02}{t}_{02}n+s_{02},\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill {\dot{t}}_{03}& =& {\omega }_{c02}{t}_{02}n-{\omega }_{03}{t}_{03},\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill {\dot{p}}_{13}& =& {\omega }_{03}{t}_{03}-p_{13}\left(1+{\omega }_{c13}n\right),\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill {\dot{p}}_{14}& =& {\omega }_{c13}{p}_{13}n-{\omega }_{14}{p}_{14},\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill {\dot{u}}_{23}& =& p_{13}-\left({\omega }_{c23}+{\omega }_{f23}\right)u_{23}n+s_{23},\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill {\dot{u}}_{24}& =& {\omega }_{14}{p}_{14}+{\omega }_{c23}{u}_{23}n-\left({\omega }_{c24}+{\omega }_{f24}\right)u_{24}n+s_{24},\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill {\dot{u}}_{25}& =& \left({\omega }_{c24}{u}_{24}-{\omega }_{f25}{u}_{25}\right)n+s_{25},\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill \dot{n}& =& \kappa n+s_n+{\chi }_2{\omega }_{dn}f+D\Delta n,\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill \dot{x}& =& -2{\omega }_{x}xn+0.074{\mathsf{\Omega }}_{f}n+s_x,\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill \dot{f}& =& {\mathsf{\Omega }}_{f}n-{\omega }_{dn}f,\hfill \end{array}$$

(26)

with the dot indicating the derivative with respect to ${\lambda }_{13}t$, and

$$\kappa =-{\mathsf{\Omega }}_c+{\chi }_1{\mathsf{\Omega }}_f,$$

(27)

$${\mathsf{\Omega }}_c={\omega }_{c02}{t}_{02}+{\omega }_{c13}{p}_{13}+{\omega }_{c23}{u}_{23}+{\omega }_{c24}{u}_{24}+{\omega }_{x}x+{\omega }_{cH}h,$$

(28)

$${\mathsf{\Omega }}_f={\omega }_{f23}{u}_{23}+{\omega }_{f24}{u}_{24}+{\omega }_{f25}{u}_{25},$$

(29)

$s_{ij}=sourc{e}_{ij}/{\lambda }_{13}{\rho }_0$, and the parameters given in

Table 1. Definitions and values of the non-dimensional parameters. The subscripts refer to the processes: c to neutron capture; x to neutron capture by either the fission products or the control neutron absorber; f to fission; $dn$ to delayed neutrons; $el$ to elastic scattering of neutrons.

Parameter	Definition	Value
${\omega }_{03}$	${\lambda }_{03}/{\lambda }_{13}$	$1.779\times {10}^3$
${\omega }_{14}$	${\lambda }_{14}/{\lambda }_{13}$	$96.627$
${\omega }_{c02}$	${\eta }_{c02}{\rho }_0/{\lambda }_{13}$	$1.624\times {10}^{11}$
${\omega }_{c13}$	${\eta }_{c13}{\rho }_0/{\lambda }_{13}$	$8.888\times {10}^{11}$
${\omega }_{c23}$	${\eta }_{c23}{\rho }_0/{\lambda }_{13}$	$1.113\times {10}^{12}$
${\omega }_{c24}$	${\eta }_{c24}{\rho }_0/{\lambda }_{13}$	$2.207\times {10}^{12}$
${\omega }_{cH}$	${\eta }_{cH}{\rho }_0/{\lambda }_{13}$	$7.483\times {10}^9$
${\omega }_x$	${\eta }_x{\rho }_0/{\lambda }_{13}$	$8.927\times {10}^{14}$
${\omega }_{f23}$	${\eta }_{f23}{\rho }_0/{\lambda }_{13}$	$1.187\times {10}^{13}$
${\omega }_{f24}$	${\eta }_{f24}{\rho }_0/{\lambda }_{13}$	$1.467\times {10}^9$
${\omega }_{f25}$	${\eta }_{f25}{\rho }_0/{\lambda }_{13}$	$1.240\times {10}^{13}$
$D{\xi }^2$	Equation (48)	$2.06\times {10}^8$
${\omega }_{dn}$	${\lambda }_{dn}/{\lambda }_{13}$	$1.012\times {10}^6$
${\omega }_{elH}$	${\eta }_{elH}{\rho }_0/{\lambda }_{13}$	$8.393\times {10}^{11}$
${\omega }_{el02}$	${\eta }_{el02}{\rho }_0/{\lambda }_{13}$	$5.431\times {10}^{10}$

. The variable h indicates the reduced density of ${}^1\mathrm{H}$ and the variable f indicates the reduced number density of the fission fragments, while the same letter as a subscript indicates fission.

The parameters ${\chi }_1=1.27$ and ${\chi }_2=0.00624$ represent the net number of the prompt and delayed neutrons emitted after fission, respectively [19]. In Equation (25) we summarize the effect of the neutron absorber ${}^{113}\mathrm{Cd}$ (obeying the equation $\dot{c}=-{\omega }_{x}cn+s_c$, the source term $s_c$ being subject to a feedback control to keep the power output at the desired level) with the neutron poison fission products y, exhibiting a high cross section for neutron capture (as ${}^{135}\mathrm{Xe}$ and ${}^{149}\mathrm{Sm}$, having total a fission yield of 0.074 and obeying the equation $\dot{y}=-{\omega }_{x}yn+0.074{\mathsf{\Omega }}_{f}n$). Representing all neutron absorbers by the single variable x and the absorption cross section ${\sigma }_x=3.12\times {10}^4\mathrm{b}$, and summing these two equations we obtain Equation (25).

In the conventional time and density units, the diffusion coefficient $D^{\prime }\left({\mathrm{length}}^2{\mathrm{time}}^{-1}\right)$ for neutrons may be expressed through the diffusion length $L_d$ as

$$D^{\prime }=\frac{1}{3}v{L}_d,$$

(30)

with

$${\rho }_0{L}_d={\left({\Sigma }_c+{\Sigma }_f+{\Sigma }_{el}\right)}^{-1},$$

(31)

$${\Sigma }_c={\sigma }_{c02}{t}_{02}+{\sigma }_{c13}{p}_{13}+{\sigma }_{c23}{u}_{23}+{\sigma }_{c24}{u}_{24}+{\sigma }_{cH}h+{\sigma }_{x}x,$$

(32)

$${\Sigma }_f={\sigma }_{f23}{u}_{23}+{\sigma }_{f24}{u}_{24}+{\sigma }_{f25}{u}_{25},{\Sigma }_{el}={\sigma }_{elH}h+{\sigma }_{el02}{t}_{02}.$$

(33)

The subscript $el$ refers to the elastic scattering of neutrons. Defining $L=v/{\lambda }_{13}$, and expressing $L_d$ in terms of the reduced parameters we have

$$L_d=\frac{L}{{\mathsf{\Omega }}_c+{\mathsf{\Omega }}_f+{\mathsf{\Omega }}_{el}},\mathrm{with}{\mathsf{\Omega }}_{el}={\omega }_{elH}h+{\omega }_{el02}{t}_{02}.$$

(34)

The reduced diffusion coefficient $D\left({\mathrm{length}}^2\right)$ thus is

$$\begin{array}{c}\hfill D=\frac{D^{\prime }}{{\lambda }_{13}}=\frac{1}{3}L{L}_d.\end{array}$$

(35)

Neglecting the contribution to the neutron flux of both the source $s_n$ and the delayed neutrons in Equation (24), we have:

$$\dot{n}=\kappa n+D\Delta n.$$

(36)

Assuming constant coefficients, a solution of the form $n={\psi }_t{\psi }_r$ [20], and separating the variables, we obtain the time equation

$${\dot{\psi }}_t=\left(\kappa -D{\xi }^2\right){\psi }_t,$$

(37)

with solution

$${\psi }_t=e^{\left(\kappa -D{\xi }^2\right)t},$$

(38)

the spatial eigenvalue equation

$$\Delta {\psi }_r=-{\xi }^2{\psi }_r,$$

(39)

with the Neumann boundary condition [21]

$$\frac{\mathrm{d}{\psi }_r}{\mathrm{d}r}=-\frac{3}{2L_t}{\psi }_r,$$

(40)

and solution

$${\psi }_r=C\frac{sin\xi r}{r}.$$

(41)

The quantity $L_t$ in Equation (40) is the transport length given by

$${\rho }_0{L}_t={\left({\Sigma }_c+{\Sigma }_f\right)}^{-1}.$$

(42)

The constant C was chosen to match the initial conditions for a homogeneous neutron number density $n_0$ in a spherical reactor of radius R

$$n_0\frac{4}{3}\pi R^3=C{\int }_0^R\frac{sin\xi r}{r}4\pi r^2\mathrm{d}r,$$

(43)

that gives

$$C=\frac{1}{3}\frac{{\xi }^2{R}^3}{sin\xi R-\xi Rcos\xi R}{n}_0.$$

(44)

The complete solution is

$$n=\frac{n_0}{3}\frac{{\xi }^2{R}^3}{sin\xi R-\xi Rcos\xi R}\frac{sin\xi r}{r}{e}^{\left(\kappa -D{\xi }^2\right)t},$$

(45)

and the Neumann boundary condition (40) becomes

$$tan\xi R=\frac{\xi R}{1-\frac{3}{2}\frac{R}{L_t}}.$$

(46)

If we approximate the time derivative of the neutron density in the whole volume in terms of the instantaneous flux at the boundary ${\phi }_R$, as $4{\pi }_R{R}^2{\phi }_R/V=3{\phi }_R/R$, and express ${\phi }_R$ according to Equation (45) as

$$\begin{array}{c}\hfill {\phi }_R=-D^{\prime }{\left|{\partial }_{r}n\right|}_R=\frac{R}{3}{D}^{\prime }{\xi }^{2}n,\end{array}$$

(47)

where ${\partial }_r\equiv \partial /\partial r$, the the time derivative of the neutron reduced density becomes $D^{\prime }{\xi }^{2}n$, and, in non-dimensional units.

$$D\Delta n\simeq D{\xi }^{2}n.$$

(48)

Solving Equation (46) for $\xi$, we may evaluate the expression (48). Within this approximation, Equation (24) can be written

$$\dot{n}=\left(\kappa -D{\xi }^2\right)n+s_n+{\chi }_2{\omega }_{dn}f,$$

(49)

and the system of Equations (17)–(26) becomes a system of ordinary differential equations.

2.2. Simulation Details

The actual value of $D{\xi }^2$ is obtained iteratively by solving Equation (46) for $\xi$ at the end of a simulation, and use this value back in Equation (49) until self consistency is reached. To validate the approximation given by Equation (48), we report in

Table 2. Quantities involved in the approximation given by Equation (48) during the 5-yr period operating with a uranium sink, with the reactor radius $R=1\mathrm{m}$, power density $50\mathrm{MW}{\mathrm{m}}^{-3}$, $t_{02}=0.05$ (≡$581.12\mathrm{kg}{\mathrm{m}}^{-3}$), $s_{02}=3.07\times {10}^{-4}$, and $s_{23}+s_{24}+s_{25}=-1.53\times {10}^{-4}$.

t	$\mathit{\xi }\left[{\mathbf{m}}^{-1}\right]$	${\mathit{R}}^{-1}sin\mathit{\xi }\mathit{R}\left[{\mathbf{m}}^{-1}\right]$	$\mathit{D}{\mathit{\xi }}^2(\times {10}^8)$
30 d	2.31	0.78	2.05
0.25 yr	2.22	0.80	2.00
0.50 yr	2.23	0.79	2.02
1.00 yr	2.25	0.78	2.04
5.00 yr	2.26	0.77	2.06

the fairly stable values of $\xi$, $sin\xi R/R$, and $D{\xi }^2$ at different times during a simulation.

For $s_{23}=s_{24}=s_{25}=0$, we can look for a stationary solution of the system of the Equations (17)–(26). The stationary reduced densities of the system (17)–(26) (denoted with a bar) are as follows

$$\begin{array}{ccc}\hfill \overline{n}& =& \frac{{\overline{s}}_{02}}{{\omega }_{c02}{\overline{t}}_{02}},\hfill \end{array}$$

(50)

$$\begin{array}{ccc}\hfill {\overline{t}}_{03}& =& \frac{{\overline{s}}_{02}}{{\omega }_{03}},\hfill \end{array}$$

(51)

$$\begin{array}{ccc}\hfill {\overline{p}}_{13}& =& {\overline{s}}_{02}A,\hfill \end{array}$$

(52)

$$\begin{array}{ccc}\hfill {\overline{p}}_{14}& =& \frac{{\overline{s}}_{02}}{{\omega }_{14}}\left(1-A\right),\hfill \end{array}$$

(53)

$$\begin{array}{ccc}\hfill {\overline{u}}_{23}& =& A{B}_{23},\hfill \end{array}$$

(54)

$$\begin{array}{ccc}\hfill {\overline{u}}_{24}& =& B_{24}\left(1-\frac{{\omega }_{f23}}{{\omega }_{c23}+{\omega }_{f23}}A\right),\hfill \end{array}$$

(55)

$$\begin{array}{ccc}\hfill {\overline{u}}_{25}& =& \frac{{\omega }_{c24}}{{\omega }_{f25}}{\overline{u}}_{24},\hfill \end{array}$$

(56)

with

$$A=\frac{{\omega }_{c02}{\overline{t}}_{02}}{{\omega }_{c02}{\overline{t}}_{02}+{\omega }_{c13}{\overline{s}}_{02}},B_{23}=\frac{{\omega }_{c02}{\overline{t}}_{02}}{{\omega }_{c23}+{\omega }_{f23}},B_{24}=\frac{{\omega }_{c02}{\overline{t}}_{02}}{{\omega }_{c24}+{\omega }_{f24}}.$$

(57)

Defining the ratio z between the power density at time t and the required value for the power density ${\rho }_w$,

$$z=\frac{{\lambda }_{13}{\rho }_0{E}_f}{{\rho }_w}{\mathsf{\Omega }}_{f}n,$$

(58)

(the average recoverable energy per fission being $E_f=170\mathrm{MeV}$ [19]), the condition $z=1$ at ${\omega }_{c02}{t}_{02}=const$ gives

$${\overline{s}}_{02}=\frac{{\rho }_w}{{\lambda }_{13}{\rho }_0{E}_f}.$$

(59)

We thus have stationary values for all state variables (except x and $s_x$) in terms of the value of ${\overline{t}}_{02}$ and ${\rho }_w$. To obtain $\overline{x}$, we equate the expression for $\kappa$ in Equations (49) and (27) and solve for x. Since $\kappa$ is also a function of x through the term ${\omega }_l$, that in turns depends on x because of $L_t$, ${\mathsf{\Omega }}_c$, ${\mathsf{\Omega }}_f$, and ${\mathsf{\Omega }}_{el}$, this operation must be done iteratively. Once $\overline{x}$ is known, Equation (25) affords the value of ${\overline{s}}_x$. We actually do not find positive solutions for both $\overline{x}$ and ${\overline{s}}_x$ unless R is very large ($R>15\mathrm{m}$), but nevertheless the above expressions turn out to be useful approximations to a true stationary state.

The values of all state variables are given in

Table 3. Values of the variables for a homogeneous reactor at pseudo-steady state with the reactor radius $R=1\mathrm{m}$, power density $50\mathrm{MW}{\mathrm{m}}^{-3}$, ${\overline{t}}_{02}=0.05$ (≡$581.12\mathrm{kg}{\mathrm{m}}^{-3}$), and ${\overline{s}}_{02}=2.03\times {10}^{-4}$ (≡$22.32\mathrm{kg}{\mathrm{m}}^{-3}{\mathrm{yr}}^{-1}$). The second row is in $\mathrm{kg}{\mathrm{m}}^{-3}$.

${\overline{\mathit{t}}}_{03}$	${\overline{\mathit{p}}}_{13}$	${\overline{\mathit{p}}}_{14}$	${\overline{\mathit{u}}}_{23}$	${\overline{\mathit{u}}}_{24}$	${\overline{\mathit{u}}}_{25}$
$1.14\times {10}^{-7}$	$1.98\times {10}^{-4}$	$4.56\times {10}^{-8}$	$6.12\times {10}^{-4}$	$3.88\times {10}^{-4}$	$6.91\times {10}^{-5}$
$1.34\times {10}^{-3}$	2.34	$5.39\times {10}^{-4}$	7.21	4.59	0.82

for a reactor with a stationary source of ${}^{232}\mathrm{Th}$${\overline{s}}_{02}=2.03\times {10}^{-4}$ ($\equiv 22.13\mathrm{kg}{\mathrm{m}}^{-3}{\mathrm{yr}}^{-1}$); the corresponding burnup rates and fractional power are given in

Table 4. Burnup rates (in $\mathrm{kg}{\mathrm{m}}^{-3}{\mathrm{yr}}^{-1}$) of ${}^{232}\mathrm{Th}$, ${}^{233}\mathrm{U}$, ${}^{234}\mathrm{U}$, and ${}^{235}\mathrm{U}$ and fractional power output of a homogeneous reactor at pseudo-steady state with ${\overline{t}}_{02}=0.05$ (≡$581.12\mathrm{kg}{\mathrm{m}}^{-3}$), the reactor radius $R=1\mathrm{m}$, and power density $50\mathrm{MW}{\mathrm{m}}^{-3}$.

Nuclide	${}^{232}\mathbf{Th}$	${}^{233}\mathbf{U}$	${}^{234}\mathbf{U}$	${}^{235}\mathbf{U}$
Burning rate	22.32	20.67	$1.23\times {10}^{-3}$	1.71
Fractional power	-	0.924	$5.51\times {10}^{-5}$	$7.63\times {10}^{-2}$

. The processed uranium with the composition given in Table 3 corresponds to the molar fractions $0.5724$, $0.3629$, and $0.06463$ for ${}^{233}\mathrm{U}$, ${}^{234}\mathrm{U}$, and ${}^{235}\mathrm{U}$, respectively.

In reality, a proper stationary-state does not exist for realistic values of R. However, a pseudo-steady state can be attained: i.e., a state where all the values of the state variables change on a time scale much longer than the typical reactor operation time. For example, 5 years after startup, we have ${\dot{t}}_{02}/t_{02}=-1.44\times {10}^{-4}{\mathrm{yr}}^{-1}$, ${\dot{u}}_{23}/u_{23}=-3.50\times {10}^{-4}{\mathrm{yr}}^{-1}$, ${\dot{u}}_{24}/u_{24}=2.39\times {10}^{-2}{\mathrm{yr}}^{-1}$, and ${\dot{u}}_{25}/u_{25}=2.96\times {10}^{-2}{\mathrm{yr}}^{-1}$.

We now address the issue of providing the fissile material to start a reactor. We begin by observing that, by including a sink of uranium of the form $s_{ij}=-s_f{\overline{s}}_{02}{\widehat{u}}_{ij}/\widehat{S}$, where $0\le s_f\le 1$ and $\widehat{S}={\widehat{u}}_{23}+{\widehat{u}}_{24}+{\widehat{u}}_{25}$, the stationary point of the system (17)–(26) (denoted by $\widehat{a}$) satisfies the nonlinear system of equations

$$\begin{array}{ccc}\hfill {\widehat{u}}_{23}\left(\widehat{S}-s_f{B}_{23}\right)-A{B}_{23}\widehat{S}& =& 0,\hfill \end{array}$$

(60)

$$\begin{array}{ccc}\hfill {\widehat{u}}_{24}\left(\widehat{S}-s_f{B}_{24}{\widehat{s}}_{02}\right)-B_{24}\widehat{S}\left(1-\frac{{\omega }_{f23}}{{\omega }_{c23}+{\omega }_{f23}}A\right)-B_{24}\frac{{\omega }_{c23}}{{\omega }_{c23}+{\omega }_{f23}}{s}_f{\widehat{u}}_{23}& =& 0,\hfill \end{array}$$

(61)

$$\begin{array}{ccc}\hfill {\widehat{u}}_{25}\left(\widehat{S}-s_f\frac{{\omega }_{c02}}{{\omega }_{f25}}{\widehat{t}}_{02}\right)-\widehat{S}{\widehat{u}}_{24}\frac{{\omega }_{c24}}{{\omega }_{f25}}& =& 0,\hfill \end{array}$$

(62)

with the corresponding expression for the stationary ${}^{232}\mathrm{Th}$ source being

$${\widehat{s}}_{02}={\overline{s}}_{02}\frac{{\omega }_{c02}{\widehat{t}}_{02}}{{\widehat{\mathsf{\Omega }}}_f}.$$

(63)

Since we do not have the option to actually operate the sink at the true stationary point, the actual expression used for the sink is

$$s_{ij}=-\frac{s_f{\overline{s}}_{02}{u}_{ij}}{u_{23}+u_{24}+u_{25}},$$

(64)

and the composition of the extracted uranium fuel is given by integration over time of the sink term

$${\tilde{u}}_{ij}=-{\int }_0^t{s}_{ij}\mathrm{d}\tau .$$

(65)

The solution of the system of Equations (60)–(62) is obtained by iteration until the composition given by Equation (65) converges to a stable value. We thus have the option of running the reactor while withdrawing, for example by continuous processing, the fissile material necessary to start a new reactor, with fuel of composition given by Equation (65).

A uranium sink with value $1.53\times {10}^{-4}$ (≡$17.03\mathrm{kg}{\mathrm{m}}^{-3}{\mathrm{yr}}^{-1}$), corresponding to $1/2$ of the source of ${}^{232}\mathrm{Th}$, gives the uranium state variables in

Table 5. Values of the reduced densities for the isotopes of uranium in a homogeneous reactor at pseudo-steady state running with a uranium sink, the reactor radius $R=1\mathrm{m}$, power density $50\mathrm{MW}{\mathrm{m}}^{-3}$, ${\widehat{t}}_{02}=0.05$ (≡$581.12\mathrm{kg}{\mathrm{m}}^{-3}$), ${\widehat{s}}_{02}=3.07\times {10}^{-4}$ (≡$33.48\mathrm{kg}{\mathrm{m}}^{-3}{\mathrm{yr}}^{-1}$), and $s_{23}+s_{24}+s_{25}=-1.53\times {10}^{-4}$ (≡$-17.03\mathrm{kg}{\mathrm{m}}^{-3}{\mathrm{yr}}^{-1}$). The second row of data is in $\mathrm{kg}{\mathrm{m}}^{-3}$.

${\widehat{\mathit{u}}}_{23}$	${\widehat{\mathit{u}}}_{24}$	${\widehat{\mathit{u}}}_{25}$
$4.26\times {10}^{-4}$	$2.15\times {10}^{-4}$	$2.56\times {10}^{-5}$
5.01	2.54	0.30

; the relevant burnup rates (mass per unit volume and unit time) and fractional power are given in

Table 6. Burnup rates (in $\mathrm{kg}{\mathrm{m}}^{-3}{\mathrm{yr}}^{-1}$) of ${}^{232}\mathrm{Th}$, ${}^{233}\mathrm{U}$, ${}^{234}\mathrm{U}$, and ${}^{235}\mathrm{U}$ and fractional power output of a homogeneous reactor at pseudo-steady state with uranium sink, the reactor radius $R=1\mathrm{m}$, and power density $50\mathrm{MW}{\mathrm{m}}^{-3}$.

Nuclide	${}^{232}\mathbf{Th}$	${}^{233}\mathbf{U}$	${}^{234}\mathbf{U}$	${}^{235}\mathbf{U}$
Burning rate	$33.48$	$21.88$	$5.45\times {10}^{-4}$	$0.50$
Fractional power	-	$0.977$	$2.44\times {10}^{-5}$	$2.25\times {10}^{-2}$

.

3. Results and Discussion

We begin detailing the results of this computational study [22] with a simulation of the reactor with a startup phase of $2.0\mathrm{d}$, where each individual source of uranium is set to the value $s_{ij}=38.92{\overline{u}}_{ij}{\tilde{u}}_{ij}/\tilde{S}$, matching the composition given by Equation (65) and totaling $s_{23}+s_{24}+s_{25}=2.02\times {10}^{-2}\equiv 6.13\mathrm{kg}{\mathrm{m}}^{-3}{\mathrm{d}}^{-1}$. At the end of the startup phase, the amount of uranium put into the reactor is $12.26\mathrm{kg}{\mathrm{m}}^{-3}$, and the uranium source is switched off. The ratio of the reduced density of the uranium species to the corresponding pseudo-stationary density is plotted in Figure 1 during the startup phase.

The corresponding reduced density of ${}^{232}\mathrm{Th}$ equals its initial value ${\overline{t}}_{02}=0.05$, equivalent to $581.12\mathrm{kg}{\mathrm{m}}^{-3}$. This value corresponds to the concentration of $2.53\mathrm{mol}/\mathrm{L}$, within the limits of thorium nitrate solubility in water given in [23] as $S_{Th}=2.615+0.010\left(T-25\right)\mathrm{mol}/\mathrm{L}$ for $5\le T\le 60$ °C. The reduced density of hydrogen, assumed to be constant, is equal to 1.07. All values of the state variables are given as ratios to the pseudo-stationary state densities with $s_{ij}=0$ at the power density ${\rho }_w=50\mathrm{MW}{\mathrm{m}}^{-3}$. In the two-day startup run with ${\rho }_w=50\mathrm{MW}{\mathrm{m}}^{-3}$, the maximum value of $s_x$ ($2.72\times {10}^{-4}$) corresponds to $3.99\times {10}^{-2}\mathrm{kg}{\mathrm{m}}^{-3}{\mathrm{d}}^{-1}$. The growth of the reduced density of all uranium isotopes is linear, an indication that the source term is always dominant with respect to the fission rates.

After the two-day startup the uranium source was switched off and the value of the source of ${}^{113}\mathrm{Cd}$ was dynamically assigned based on the ratio z of the current power density to the set power density, i.e., $s_x=0.01\left(z-1\right)$. The ratios of the reduced densities of the uranium species to the corresponding pseudo-steady-state density are reported in Figure 2 for a period of 5 years, and the value of the average neutron flux in this time period is $2.42\times {10}^{18}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}$.

The reduced density of ${}^{233}\mathrm{U}$ at the time of switch-off of the uranium source is ≈$20\%$ above its stationary value, while the other isotopes, ${}^{234}\mathrm{U}$ and ${}^{235}\mathrm{U}$, reach the equilibrium value just at the end of the 5-yr time period.

Operating the uranium sink for a period of 5 years the system stays critical by increasing the neutron flux to $4.18\times {10}^{18}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}$ and lowering the uranium state variables to the values given in

Table 7. Reduced number density of the uranium isotopes for a homogeneous reactor at the end a 5-yr period operating with a uranium sink, the reactor radius $R=1\mathrm{m}$, power density $50\mathrm{MW}{\mathrm{m}}^{-3}$, $t_{02}=0.05$ (≡$581.12\mathrm{kg}{\mathrm{m}}^{-3}$), $s_{02}=3.07\times {10}^{-4}$ (≡$33.48\mathrm{kg}{\mathrm{m}}^{-3}{\mathrm{yr}}^{-1}$), and $s_{23}+s_{24}+s_{25}=-1.53\times {10}^{-4}$ (≡$-17.03\mathrm{kg}{\mathrm{m}}^{-3}{\mathrm{yr}}^{-1}$). The second row of data is in $\mathrm{kg}{\mathrm{m}}^{-3}$.

${\mathit{u}}_{23}$	${\mathit{u}}_{24}$	${\mathit{u}}_{25}$
$3.66\times {10}^{-4}$	$7.31\times {10}^{-5}$	$8.11\times {10}^{-6}$
4.31	0.87	$9.63\times {10}^{-2}$

. Plots of the ratios of the reduced densities of the various species to the corresponding pseudo-steady-state density are reported in Figure 3. The withdrawn uranium has the composition $u_{23}=5.89\times {10}^{-3},u_{24}=1.19\times {10}^{-3},\mathrm{and}{u}_{25}=1.29\times {10}^{-4}$, corresponding to the molar fractions $\chi =0.82,0.17$, and $0.018$ for ${}^{233}\mathrm{U}$, ${}^{234}\mathrm{U}$, and ${}^{235}\mathrm{U}$, respectively. At these molar ratios and the number density of metallic uranium ${\rho }_U=4.83\times {10}^{28}{\mathrm{m}}^{-3}$, the fission rate per fast neutron would be $\eta {\rho }_U\chi =1.32\times {10}^8,1.92\times {10}^7$, and $1.01\times {10}^6{\mathrm{s}}^{-1}$ for ${}^{233}\mathrm{U}$, ${}^{234}\mathrm{U}$, and ${}^{235}\mathrm{U}$, respectively. All these values are able to sustain a fission chain reaction in a metallic uranium sphere with the above isotopic composition. The corresponding ratios between the reduced densities at the end of the five years and the pseudo-steady state reduced densities of uranium are $u_{23}/{\widehat{u}}_{23}=0.86$, $u_{24}/{\widehat{u}}_{24}=0.34$, and $u_{25}/{\widehat{u}}_{25}=0.32$. The withdrawal of the uranium isotopes reduces the reduced density of ${}^{232}\mathrm{Th}$ to $95.9\%$ of its initial value.

Based on a cross section of $1.80\mathrm{mb}$ for the $(n,2n)$ reaction on ${}^{233}\mathrm{U}$ by $6\mathrm{MeV}$ neutrons [3], the formation of ${}^{232}\mathrm{U}$ in a period of five years through the $\left(n,2n\right)$ channel gives $u_{22}/\left(u_{23}+u_{24}+u_{25}\right)=2.70\times {10}^{-4}$, corresponding to a specific activity of $2.18\times {10}^{11}\mathrm{Bq}{\mathrm{kg}}^{-1}$ in the processed uranium with the composition of the withdrawn fuel. Kang and von Hippel estimated the exposure at a distance of one meter from a 5 kg uranium sphere with 0.4 % of ${}^{232}\mathrm{U}$ in secular equilibrium with ${}^{208}\mathrm{Tl}$ to be $0.076\mathrm{Sv}{\mathrm{h}}^{-1}{\mathrm{kg}}^{-1}$. The other uranium isotopes would contribute a specific activity of $2.85\times {10}^{11}$, $3.74\times {10}^{10}$, and $1.41\times {10}^6\mathrm{Bq}{\mathrm{kg}}^{-1}$ for ${}^{233}\mathrm{U}$, ${}^{234}\mathrm{U}$, and ${}^{235}\mathrm{U}$ respectively. For comparison, the specific activity of ${}^{238}\mathrm{U}$ is $1.24\times {10}^7\mathrm{Bq}{\mathrm{kg}}^{-1}$.

The ratio ${\left(s_{23}+s_{24}+s_{25}\right)}_{\mathrm{startup}}t/{\left(s_{23}+s_{24}+s_{25}\right)}_{\mathrm{sink}}=2.02\times {10}^{-2}t/1.53\times {10}^{-4}$ of the total uranium input at the end of the two-day startup period to the uranium sink rate of a reactor with $s_f=0.5$ at the end of the five-year period following startup is $0.72\mathrm{yr}$, that is, a reactor running with a uranium sink for $0.72\mathrm{yr}$ would breed enough uranium to start a new reactor of the same size in 2 days.

We conclude by comparing the yield of this simulated reactor with the more common reactors based on the uranium fuel cycle. The ratio of the fissioned mass to the energy produced for the homogeneous reactor described above as $1.97\times {10}^{-13}\mathrm{kg}{\mathrm{J}}^{-1}$, while the corresponding value for the fission of ${}^{235}\mathrm{U}$ enriched at $3\%$ is $4.20\times {10}^{-13}\mathrm{kg}{\mathrm{J}}^{-1}$.

4. Conclusions

This work explores the possibility of energy production by a fission reactor in the absence of a continuous supply of fissile isotopes, since the dynamics of a homogeneous light-water thorium reactor can be exploited for breeding the uranium isotopes necessary for fission. In fact, after a short startup period of two days requiring fissile uranium isotopes (for a total inventory of $12.26\mathrm{kg}{\mathrm{m}}^{-3}$), the modeled reactor can operate at a nearly steady state with a constant source of ${}^{232}\mathrm{Th}$, both with and without a uranium sink, i.e., a steady withdrawal of all uranium isotopes. After the fissile ${}^{233}\mathrm{U}$ or ${}^{235}\mathrm{U}$ has been used for the startup of the first reactor, its operation with a continuous uranium sink for 0.72 yr can afford the necessary fissile uranium isotopes to start a new reactor, allowing the continuous production of energy from only ${}^{232}\mathrm{Th}$. However, the limitation to nuclear weapon proliferation does not appear to be ensured by the molar fraction of ${}^{234}\mathrm{U}$ in the extracted uranium fuel.