The TIP for multigroup neutron diffusion equation in 3D hexagonal geometry

In 3D Cartesian coordinates (Fig. 1), the multigroup neutron diffusion equation within node k can be written as $\begin{array}{ll}-D_g^k\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+\frac{{\partial }^2}{\partial z^2}\right){\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\hfill & +{\Sigma }_{r\mathrm{,}g}^k{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\hfill \\ \hfill & =Q_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\hfill \end{array}$ (1) where: $Q_g^k\left(x\mathrm{,}y\mathrm{,}z\right)=\left[{\displaystyle \sum _{g^{\prime }\ne g}}{\Sigma }_{g^{\prime }\to g}^k+\frac{{\chi }_g}{k_{\text{eff}}}{\displaystyle \sum _{g^{\prime }=1}^G}\nu {\Sigma }_{f\mathrm{,}{g}^{\prime }}^k\right]{\Phi }_{g^{\prime }}^k\left(x\mathrm{,}y\mathrm{,}z\right)$ (2) Integrating Eq. (1) over the volume of node k leads to the following nodal balance equation: $\begin{array}{ll}\frac{1}{3a}{\displaystyle \sum _{x=x\mathrm{,}u\mathrm{,}\nu }}\left[J_{gx}^k\left(a\right)-J_{gx}^k\left(-a\right)\right]\hfill & +\frac{1}{2a_z^k}\left[J_{gz}^k\left(a_z^k\right)-J_{gz}^k\left(-a_z^k\right)\right]\hfill \\ \hfill & +{\Sigma }_{r\mathrm{,}g}^k{\overline{\Phi }}_g^k={\overline{Q}}_g^k\hfill \end{array}$ (3) where: ${\overline{\Phi }}_g^k=\frac{1}{V_k}{\displaystyle {\int }_{-a_z^k}^{a_z^k}}\text{d}z{\displaystyle {\int }_{-a}^a}\text{d}x{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\text{d}y$ (4) $\begin{array}{c}{\overline{Q}}_g^k=\frac{1}{V_k}{\displaystyle {\int }_{-a_z^k}^{a_z^k}}\text{d}z{\displaystyle {\int }_{-a}^a}\text{d}x{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}{Q}_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\text{d}y\\ V_k={\displaystyle {\int }_{-a_z^k}^{a_z^k}}\text{d}z{\displaystyle {\int }_{-a}^a}\text{d}x{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}\text{d}y=4\sqrt{3}{a}^2{a}_z^k\end{array}$ (5) $y_s\left(x\right)=\frac{1}{\sqrt{3}}(2a-|x|)$ (6) where Vk is the volume of node k; ${\overline{\Phi }}_g^k$ and ${\overline{Q}}_g^k$ are the neutron flux and neutron source averaged over Vk, respectively; a is the lattice half pitch; and $a_z^k$ is the haft height of node k.

Fig. 1Full-screen View

Hexagonal coordinate system

In Eq. (3), $\begin{array}{l}{J}_{gx}^k\left(\pm a\right)=\frac{1}{2a_z^k}{\displaystyle {\int }_{-a_z^k}^{a_z^k}}\text{d}z\\ \cdot {\left[\frac{1}{2y_s\left(x\right)}{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}\left(-D_g^k\frac{\partial }{\partial x}{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\right)\text{d}y\right]}_{x=\pm a}\end{array}$ (7) and $J_{gz}^k\left(\pm a_z^k\right)=\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}\text{d}x{\left[{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}\left(-D_g^k\frac{\partial }{\partial z}{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\right)\text{d}y\right]}_{z=\pm a_z^k}$ (8) are the net neutron currents averaged at both surfaces in the $x(x=x\mathrm{,}u\mathrm{,}\nu )$ and z-directions, respectively.

The transverse integration procedure for Eq. (1) leads to the following equation: $\begin{array}{ll}-D_g^k\frac{{\text{d}}^2}{\text{d}{x}^2}{\Phi }_{gx}^k\left(x\right)+{\Sigma }_{r\mathrm{,}g}^k{\Phi }_{gx}^k\left(x\right)\hfill & =Q_{gx}^k\left(x\right)-L_{gx}^k\left(x\right)\hfill \\ \hfill & \mathrm{,}x(x=x\mathrm{,}u\mathrm{,}\nu )\hfill \end{array}$ (9) where ${\Phi }_{gx}^k\left(x\right)={\displaystyle {\int }_{-a_z^k}^{a_z^k}}\text{d}z{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\text{d}y$ (10) $Q_{gx}^k\left(x\right)={\displaystyle {\int }_{-a_z^k}^{a_z^k}}\text{d}z{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}{Q}_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\text{d}y$ (11) and $J_{gx}^k\left(x\right)={\displaystyle {\int }_{-a_z^k}^{a_z^k}}\text{d}z{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}\left[-D_g^k\frac{\partial }{\partial x}{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\right]\text{d}y$ (12) are the transverse integrated neutron flux, transverse integrated neutron source, and transverse integrated neutron current, respectively. $L_{gx}^k$ is the transverse leakage given by $\begin{array}{c}{L}_{gx}^k\left(x\right)=\frac{2}{\sqrt{3}}\left[J_{gn}^k\left(x\mathrm{,}{y}_s\left(x\right)\right)-J_{gn}^k\left(x\mathrm{,}-y_s\left(x\right)\right)\right]\\ +J_{gz}^k\left(x\mathrm{,}{a}_z^k\right)-J_{gz}^k\left(x\mathrm{,}-a_z^k\right)\\ -D_g^k\frac{sgn\left(x\right)}{\sqrt{3}}\frac{d}{dx}\left[{\Phi }_g^k\left(x\mathrm{,}{y}_s\left(x\right)\right)+{\Phi }_g^k\left(x\mathrm{,}-y_s\left(x\right)\right)\right]\\ -D_g^k\frac{2\delta (x)}{\sqrt{3}}\left[{\Phi }_g^k\left(x\mathrm{,}{y}_s\left(x\right)\right)+{\Phi }_g^k\left(x\mathrm{,}-y_s\left(x\right)\right)\right]\end{array}$ (13) Here, the following relationship between the transverse integrated flux and transverse integrated current is used: $\begin{array}{ll}{J}_{gx}^k\left(x\right)=\hfill & -D_g^k\frac{d}{dx}{\Phi }_{gx}^k\left(x\right)\hfill \\ \hfill & +D_g^k{y^{\prime }}_s\left(x\right)\left[{\Phi }_g^k\left(x\mathrm{,}{y}_s\left(x\right)\right)+{\Phi }_g^k\left(x\mathrm{,}-y_s\left(x\right)\right)\right]\hfill \end{array}$ (14) In Eq. (13) ${\Phi }_g^k\left(x\mathrm{,}\pm y_s\left(x\right)\right)={\displaystyle {\int }_{-a_z^k}^{a_z^k}}{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right){|}_{y=\pm y_s\left(x\right)}\text{d}z$ (15) $J_{gz}^k\left(x\mathrm{,}\pm a_z^k\right)={\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}{\left[-D_g^k\frac{\partial }{\partial z}{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\right]}_{z=\pm a_z^k}\text{d}y$ (16) $J_{gn}^k\left(x\mathrm{,}\pm y_s\left(x\right)\right)={\displaystyle {\int }_{-a_z^k}^{a_z^k}}{\left[-D_g^k\frac{\partial }{\partial n}{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\right]}_{y=\pm y_s\left(x\right)}\text{d}z$ (17) and ${y^{\prime }}_s\left(x\right)=-\frac{1}{\sqrt{3}}sgn(x)\mathrm{,}{{y^{\prime }}^{\prime }}_s(x)=-\frac{2}{\sqrt{3}}\delta (x)$.

Similarly to Eq. (9), the transverse integrated equation along the z-direction is given by $-D_g^k\frac{{\text{d}}^2}{\text{d}{z}^2}{\Phi }_{gz}^k\left(z\right)+{\Sigma }_{r\mathrm{,}g}^k{\Phi }_{gz}^k\left(z\right)=Q_{gz}^k\left(z\right)-L_{gz}^k\left(z\right)$ (18) where, ${\Phi }_{gz}^k\left(z\right)={\displaystyle {\int }_{-a}^a}\text{d}x{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\text{d}y$ (19) $Q_{gz}^k\left(z\right)={\displaystyle {\int }_{-a}^a}\text{d}x{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}{Q}_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\text{d}y$ (20) $L_{gz}^k\left(z\right)={\displaystyle \sum _{x=x\mathrm{,}u\mathrm{,}\nu }}\left(J_{gx}^k\left(a\mathrm{,}z\right)-J_{gx}^k\left(a\mathrm{,}z\right)\right)$ (21) $J_{gx}^k\left(\pm a\mathrm{,}z\right)={\left[{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}\left(-D_g^k\frac{\partial }{\partial x}{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\right)\text{d}y\right]}_{x=\pm a}$ (22) We define the transverse averaged neutron flux, transverse averaged neutron current, and transverse averaged neutron source as follows: $\begin{array}{l}{\overline{\Phi }}_{gx}^k\left(x\right)=\frac{1}{4y_s\left(x\right)a_z^k}{\displaystyle {\int }_{-a_z^k}^{a_z^k}}\text{d}z{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\text{d}y=\frac{1}{4y_s\left(x\right)a_z^k}{\Phi }_{gx}^k\left(x\right)\hfill \end{array}$ (23) $\begin{array}{l}\begin{array}{l}{\overline{J}}_{gx}^k\left(x\right)=\frac{1}{4y_s\left(x\right)a_z^k}{\displaystyle {\int }_{-a_z^k}^{a_z^k}}\text{d}z{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}\left[-D_g^k\frac{\partial }{\partial x}{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\right]\\ \text{d}y=\frac{1}{4y_s\left(x\right)a_z^k}{J}_{gx}^k\left(x\right)\end{array}\hfill \end{array}$ (24) $\begin{array}{l}{\overline{Q}}_{gx}^k\left(x\right)=\frac{1}{4y_s\left(x\right)a_z^k}{\displaystyle {\int }_{-a_z^k}^{a_z^k}}\text{d}z{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}{Q}_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\text{d}y=\frac{1}{4y_s\left(x\right)a_z^k}{Q}_{gx}^k\left(x\right)\hfill \end{array}$ (25) $\begin{array}{l}{\overline{\Phi }}_{gz}^k\left(z\right)=\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}\text{d}x{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\text{d}y=\frac{1}{2\sqrt{3}{a}^2}{\Phi }_{gz}^k\left(z\right)\hfill \end{array}$ (26) $\begin{array}{l}{\overline{Q}}_{gz}^k\left(z\right)=\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}\text{d}x{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}{Q}_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\text{d}y=\frac{1}{2\sqrt{3}{a}^2}{Q}_{gz}^k\left(z\right)\hfill \end{array}$ (27) Substituting Eqs.(23)-(27) into Eqs.(9) and (18), the following equations for the transverse averaged neutron fluxes are obtained $\begin{array}{ll}-D_g^k\frac{{\text{d}}^2}{\text{d}{x}^2}{\overline{\Phi }}_{gx}^k\left(x\right)+{\Sigma }_{r\mathrm{,}g}^k{\overline{\Phi }}_{gx}^k\left(x\right)\hfill & ={\overline{Q}}_{gx}^k\left(x\right)-{\overline{L}}_{gx}^k\left(x\right)\hfill \\ \hfill & \mathrm{,}x(x=x\mathrm{,}u\mathrm{,}\nu )\hfill \end{array}$ (28) $-D_g^k\frac{{\text{d}}^2}{\text{d}{z}^2}{\overline{\Phi }}_{gz}^k\left(z\right)+{\Sigma }_{r\mathrm{,}g}^k{\overline{\Phi }}_{gz}^k\left(z\right)={\overline{Q}}_{gz}^k\left(z\right)-{\overline{L}}_{gz}^k\left(z\right)$ (29) where $\begin{array}{c}{\overline{L}}_{gx}^k\left(x\right)=\frac{1}{\sqrt{3}{y}_s\left(x\right)}\left[{\overline{J}}_{gn}^k\left(x\mathrm{,}{y}_s\left(x\right)\right)-{\overline{J}}_{gn}^k\left(x\mathrm{,}-y_s\left(x\right)\right)\right]\\ +\frac{1}{2a_z^k}\left[{\overline{J}}_{gz}^k\left(x\mathrm{,}{a}_z^k\right)-{\overline{J}}_{gz}^k\left(x\mathrm{,}-a_z^k\right)\right]\\ -D_g^k\frac{sgn\left(x\right)}{2\sqrt{3}{y}_s\left(x\right)}\frac{d}{dx}\left[{\overline{\Phi }}_g^k\left(x\mathrm{,}{y}_s\left(x\right)\right)+{\overline{\Phi }}_g^k\left(x\mathrm{,}-y_s\left(x\right)\right)-4{\overline{\Phi }}_{gx}^k\left(x\right)\right]\\ -D_g^k\frac{sgn\left(x\right)}{2\sqrt{3}{y}_s\left(x\right)}\left[{\overline{\Phi }}_g^k\left(x\mathrm{,}{y}_s\left(x\right)\right)+{\overline{\Phi }}_g^k\left(x\mathrm{,}-y_s\left(x\right)\right)-2{\overline{\Phi }}_{gx}^k\left(x\right)\right]\end{array}$ (30) ${\overline{L}}_{gz}^k\left(z\right)=\frac{1}{3a}{\displaystyle \sum _{x=x\mathrm{,}u\mathrm{,}\nu }}\left[{\overline{J}}_{gx}^k\left(a\mathrm{,}z\right)-{\overline{J}}_{gx}^k\left(-a\mathrm{,}z\right)\right]$ (31) are the transverse averaged leakages along the $x(x=x\mathrm{,}u\mathrm{,}\nu )$ and z-directions, respectively, and $\begin{array}{l}{\overline{\Phi }}_g^k\left(x\mathrm{,}\pm y_s\left(x\right)\right)=\frac{1}{2a_z}{\Phi }_g^k\left(x\mathrm{,}\pm y_s\left(x\right)\right)\mathrm{,}\\ {\overline{J}}_{gn}^k\left(x\mathrm{,}\pm y_s\left(x\right)\right)=\frac{1}{2a_z}{J}_{gn}^k\left(x\mathrm{,}\pm y_s\left(x\right)\right)\mathrm{,}\\ {\overline{J}}_{gz}^k\left(x\mathrm{,}\pm a_z^k\right)=\frac{1}{2y_s\left(x\right)}{J}_{gz}^k\left(x\mathrm{,}\pm a_z^k\right)\mathrm{,}\\ {\overline{J}}_{gx}^k\left(\pm a\mathrm{,}z\right)=\left[\frac{1}{2y_s\left(x\right)}{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}}\left(-D_g^k\frac{\partial }{\partial x}{\Phi }_g^k\left(x\mathrm{,}y\mathrm{,}z\right)\right)\text{d}y\right]\end{array}$ (32)

Solution to the transverse integrated equation using Green’s function under the Neumann boundary condition

Equation (28) can be written as follows (the index g and k and the - symbol are omitted for convenience): $-D\frac{{\text{d}}^2}{\text{d}{x}^2}\Phi (x)+{\Sigma }_r\Phi (x)=Q(x)-L(x)\mathrm{,}x=x\mathrm{,}u\mathrm{,}\nu$ (33) The Green’s function is the solution of the following equation, which corresponds to Eq. (33) $-D\frac{{\text{d}}^2}{\text{d}{x}^2}G\left(x\mathrm{,}{x}_0\right)+{\Sigma }_{r}G\left(x\mathrm{,}{x}_0\right)=\delta \left(x-x_0\right)$ (34) Using Neumann boundary conditions, ${\left[\frac{\text{d}G\left(x\mathrm{,}{x}_0\right)}{\text{d}x}\right]}_{x_0=\pm a}$ and Eqs.(14), (23), and (24), the following equation for the transverse averaged flux is obtained $\begin{array}{c}\Phi (x)={\displaystyle {\int }_{-a}^a}\left(Q\left(x_0\right)-L\left(x_0\right)\right)G\left(x\mathrm{,}{x}_0\right)\text{d}{x}_0-J_x\left(a\right)G(x\mathrm{,}a)\\ +J_x\left(-a\right)G(x\mathrm{,}-a)-\frac{D}{2a}\left[{\Phi }_1+{\Phi }_2-2{\Phi }_x\left(a\right)\right]G(x\mathrm{,}a)\\ +\frac{D}{2a}\left[{\Phi }_4+{\Phi }_5-2{\Phi }_x\left(-a\right)\right]G(x\mathrm{,}-a)\end{array}$ (35) where $\begin{array}{l}{\Phi }_1=\Phi \left(a\mathrm{,}{y}_s\left(a\right)\right)\mathrm{,}{\Phi }_2=\Phi \left(a\mathrm{,}{y}_s\left(a\right)\right)\mathrm{,}\text{ and}\\ {\Phi }_3=\Phi \left(\mathrm{0,}-y_s\left(0\right)\right)\mathrm{,}{\Phi }_4=\Phi \left(-a\mathrm{,}-y_s\left(a\right)\right)\mathrm{,}\\ {\Phi }_5=\Phi \left(-a\mathrm{,}{y}_s\left(a\right)\right)\mathrm{,}{\Phi }_6=\Phi \left(\mathrm{0,}{y}_s\left(0\right)\right)\end{array}$ are the corner point flux values of the hexagonal node.

We denote the transverse leakage by the sum of the individual terms as follows: $L(x)=L_{J\mathrm{,}y}\left(x\right)+L_{J\mathrm{,}z}\left(x\right)+L_{\text{sgn}}\left(x\right)+L_{\delta }\left(x\right)$ (36) where, $\begin{array}{ll}\hfill & L_{J\mathrm{,}y}\left(x\right)=\frac{1}{\sqrt{3}{y}_s\left(x\right)}\left[J_n\left(x\mathrm{,}{y}_s\left(x\right)\right)-J_n\left(x\mathrm{,}-y_s\left(x\right)\right)\right]\mathrm{,}\hfill \\ \hfill & L_{J\mathrm{,}z}\left(x\right)=\frac{1}{2a_z}\left[J_z\left(x\mathrm{,}{a}_z\right)-J_z\left(x\mathrm{,}-a_z\right)\right]\mathrm{,}\hfill \\ \hfill & L_{\text{sgn}}\left(x\right)=-D\frac{sgn(x)}{2\sqrt{3}{y}_s\left(x\right)}\frac{d}{dx}\left[\Phi \left(x\mathrm{,}{y}_s\left(x\right)\right)+\Phi \left(x\mathrm{,}-y_s\left(x\right)\right)-4\Phi (x)\right]\mathrm{,}\hfill \\ \hfill & L_{\delta \left(x\right)}=-D\frac{\delta x}{\sqrt{3}{y}_s\left(x\right)}\left[\Phi \left(x\mathrm{,}{y}_s\left(x\right)\right)+\Phi \left(x\mathrm{,}-y_s\left(x\right)\right)-2\Phi (x)\right]\mathrm{,}\hfill \end{array}$ (37) where $L_{\text{sgn}}\left(x\right)$ and $L_{\delta \left(x\right)}$ are the singular (discontinuous) terms that appear in the hexagonal geometry.

Using the following linear approximations for $\Phi (x)\mathrm{,}\Phi \left(x\mathrm{,}{y}_s\left(x\right)\right)\mathrm{,}\Phi \left(x\mathrm{,}-y_s\left(x\right)\right)$ $\Phi (x)=\frac{\Phi \left(x\mathrm{,}{y}_s\left(x\right)\right)+\Phi \left(x\mathrm{,}-y_s\left(x\right)\right)}{2}$ (38) $\Phi \left(x\mathrm{,}{y}_s\left(x\right)\right)=\{\begin{array}{l}\frac{{\Phi }_1-{\Phi }_6}{a}x+{\Phi }_{6}x>0\hfill \\ \frac{{\Phi }_5-{\Phi }_6}{a}x+{\Phi }_{6}x<0\hfill \end{array}$ (39) $\Phi \left(x\mathrm{,}-y_s\left(x\right)\right)=\{\begin{array}{l}\frac{{\Phi }_2-{\Phi }_3}{a}x+{\Phi }_3\mathrm{,}x>0\hfill \\ \frac{{\Phi }_3-{\Phi }_4}{a}x+{\Phi }_3\mathrm{,}x<0\hfill \end{array}$ (40) in Eq. (35), the integral for the discontinuous terms leads to the following equation $\begin{array}{c}\Phi (x)={\displaystyle {\int }_{-a}^a}\left(Q\left(x_0\right)-L_J\left(x_0\right)\right)G\left(x\mathrm{,}{x}_0\right)\text{d}{x}_0-J_x\left(a\right)G(x\mathrm{,}a)+J_x\left(-a\right)G(x\mathrm{,}-a)\\ -\frac{D}{2a}\left[{\Phi }_1+{\Phi }_2-2{\Phi }_x\left(a\right)\right]G(x\mathrm{,}a)+\frac{D}{2a}\left[{\Phi }_4+{\Phi }_5-2{\Phi }_x\left(-a\right)\right]G(x\mathrm{,}-a)\\ +\frac{D}{2a}\left[{\Phi }_3+{\Phi }_6-2{\Phi }_x\left(0\right)\right]G(x\mathrm{,0})+\frac{D}{2a}\left({\Phi }_3+{\Phi }_6-{\Phi }_1-{\Phi }_2\right){\displaystyle {\int }_0^a}\frac{G\left(x\mathrm{,}{x}_0\right)}{\left(2a-x_0\right)}\text{d}{x}_0\\ +\frac{D}{2a}\left({\Phi }_3+{\Phi }_6-{\Phi }_4-{\Phi }_5\right){\displaystyle {\int }_{-a}^0}\frac{G\left(x\mathrm{,}{x}_0\right)}{\left(2a+x_0\right)}\text{d}{x}_0\end{array}$ (41) where $L_J\left(x\right)=L_{J\mathrm{,}y}\left(x\right)+L_{J\mathrm{,}z}\left(x\right)$ denotes the practical leakage term in the y and z-directions.

Substituting Eq. (38), Eq. (39) and (40) again into Eq. (41), the fourth, fifth, and sixth terms in right side of Eq. (41) become zero, and thus, the following equation for the transverse averaged flux is obtained: $\begin{array}{c}\Phi (x)={\displaystyle {\int }_{-a}^a}\left(Q\left(x_0\right)-L_J\left(x_0\right)\right)G\left(x\mathrm{,}{x}_0\right)\text{d}{x}_0-J_x\left(a\right)G(x\mathrm{,}a)+J_x\left(-a\right)G(x\mathrm{,}-a)\\ +\frac{D}{2a}\left({\Phi }_3+{\Phi }_6-{\Phi }_1-{\Phi }_2\right){\displaystyle {\int }_0^a}\frac{G\left(x\mathrm{,}{x}_0\right)}{\left(2a-x_0\right)}\text{d}{x}_0+\frac{D}{2a}\left({\Phi }_3+{\Phi }_6-{\Phi }_4-{\Phi }_5\right){\displaystyle {\int }_{-a}^0}\frac{G\left(x\mathrm{,}{x}_0\right)}{\left(2a+x_0\right)}\text{d}{x}_0\end{array}$ (42)

Constitution of the closed system of equations

The transverse-averaged flux, source, and leakage terms are expanded using orthogonal polynomials: $\begin{array}{l}\Phi (x)={\displaystyle \sum _{n=0}^{n=2}}{\Phi }_n{\omega }_n\left(x\right)\\ Q(x)={\displaystyle \sum _{n=0}^{n=2}}{Q}_n{\omega }_n\left(x\right)\\ L_J\left(x\right)={\displaystyle \sum _{n=0}^{n=2}}{L}_n{\omega }_n\left(x\right)\end{array}$ (43) where ${\omega }_0\left(x\right)=\mathrm{1,}{\omega }_1\left(x\right)=3\sqrt{\frac{2}{5}}\frac{x}{a}\mathrm{,}{\omega }_2\left(x\right)=\sqrt{\frac{20}{127}}\left(\frac{9x^2}{a^2}-\frac{5}{2}\right)$ (44) are orthogonal [18]: $\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}2y_s\left(x\right){\omega }_m\left(x\right){\omega }_n\left(x\right)\text{d}x={\delta }_{mn}$ (45) where ${\displaystyle {\int }_{-a}^{a}dx};{\displaystyle {\int }_{-y_s\left(x\right)}^{y_s\left(x\right)}dy}={\displaystyle {\int }_{-a}^{a}2y_s}\left(x\right)dx=2\sqrt{3}{a}^2$ denotes the area of a hexagon with lattice pitch $2a$.

From Eq. (45), the expansion coefficients in Eq. (43) are given by $\begin{array}{l}{\Phi }_n=\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}2y_s\left(x\right)\Phi (x){\omega }_n\left(x\right)\text{d}x\\ Q_n=\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}2y_s\left(x\right)Q(x){\omega }_n\left(x\right)\text{d}x\end{array}$ (46) $L_n=\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}2y_s\left(x\right)L(x){\omega }_n\left(x\right)\text{d}x$ (47) where Qn is related to ${\Phi }_n$ via Eq. (2).

The zeroth order moment of the transverse averaged flux equals the node-averaged flux, that is, $\begin{array}{ll}{\Phi }_0\hfill & =\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}2y_s\left(x\right)\Phi (x){\omega }_0\left(x\right)\text{d}x\hfill \\ \hfill & =\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}2y_s\left(x\right)\Phi (x)\text{d}x=\overline{\Phi }\hfill \end{array}$ (48) By contrast, the zeroth order moments of the leakage terms in the y and z coordinate directions are given by $\begin{array}{ll}{L}_{J\mathrm{,}{y}_0}\hfill & =\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}2y_s\left(x\right)L_{J\mathrm{,}y}\left(x\right)\text{d}x\hfill \\ \hfill & =\frac{1}{3a}{\displaystyle \sum _{s=u\mathrm{,}\nu }}\left[J_s\left(a\right)-J_s\left(-a\right)\right]={\overline{L}}_{J\mathrm{,}y}\hfill \end{array}$ (49) $\begin{array}{ll}{L}_{J\mathrm{,}{z}_0}\hfill & =\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}2y_s\left(x\right)L_{J\mathrm{,}z}\left(x\right)\text{d}x\hfill \\ \hfill & =\frac{1}{3a}{\displaystyle \sum _{s=u\mathrm{,}\nu }}\left[J_z\left(a_z\right)-J_z\left(-a_z\right)\right]={\overline{L}}_{J\mathrm{,}z}\hfill \end{array}$ (50) The first- and second-order moments can be obtained using an extrapolation procedure for two adjacent nodes, as in a square node.

The transverse-averaged flux at the boundary $x=\pm a$ is as follows: $\begin{array}{ll}\Phi (\pm a)\hfill & ={\displaystyle {\int }_{-a}^a}(Q\left(x_0\right)-L_J\left(x_0\right)G\left(x_0\mathrm{,}\pm a\right)\text{d}{x}_0\hfill \\ \hfill & -J_x\left(a\right)G(a\mathrm{,}\pm a)+J_x\left(-a\right)G(-a\mathrm{,}\pm a)\hfill \\ \hfill & +\frac{D}{2a}\left({\Phi }_3+{\Phi }_6-{\Phi }_1-{\Phi }_2\right){\displaystyle {\int }_0^a}\frac{G\left(x_0\mathrm{,}\pm a\right)}{2a-x_0}\text{d}{x}_0\hfill \\ \hfill & +\frac{D}{2a}\left({\Phi }_3+{\Phi }_6-{\Phi }_4-{\Phi }_5\right){\displaystyle {\int }_{-a}^0}\frac{G\left(x_0\mathrm{,}\pm a\right)}{2a+x_0}\text{d}{x}_0\hfill \end{array}$ (51) Inserting Eq. (43) for $Q\left(x_0\right)$ and $L_J\left(x_0\right)$ into Eq. (51) leads to the following equation $\begin{array}{l}\Phi (a)=\left[G_I+\right]-J_x\left(a\right)G(a\mathrm{,}a)+J_x\left(-a\right)G(-a\mathrm{,}a)+\left[G_S+\right]\hfill \end{array}$ (52) $\begin{array}{ll}\Phi (-a)\hfill & =[G_I-]-J_x\left(a\right)G(a\mathrm{,}-a)+J_x\left(-a\right)G(-a\mathrm{,}-a)\hfill \\ \hfill & +[G_S-]\hfill \end{array}$ (53) where, $=[G_I\pm ]={\displaystyle \sum _{n=0}^2}\left(Q_n-L_n\right){\left[G_I^{\pm }\right]}_n$ (54) $\begin{array}{l}[G_S\pm ]=\frac{D}{2a}\left({\Phi }_3+{\Phi }_6-{\Phi }_1-{\Phi }_2\right){\displaystyle {\int }_0^a}\frac{G(x\mathrm{,}\pm a)}{(2a-x)}\text{d}x+\frac{D}{2a}\left({\Phi }_3+{\Phi }_6-{\Phi }_4-{\Phi }_5\right){\displaystyle {\int }_{-a}^0}\frac{G(x\mathrm{,}\pm a)}{(2a+x)}\text{d}x\\ =\frac{1}{2a\kappa \sinh 2\kappa a}\left[\left({\Phi }_3+{\Phi }_6-{\Phi }_1-{\Phi }_2\right){\displaystyle {\int }_0^a}\frac{\cosh \kappa (a\pm x)}{(2a-x)}\text{d}x+\frac{D}{2a}\left({\Phi }_3+{\Phi }_6-{\Phi }_4-{\Phi }_5\right){\displaystyle {\int }_{-a}^0}\frac{\cosh \kappa (a\pm x)}{2a+x}\text{d}x\right]\end{array}$ (55) $\begin{array}{l}[G_I^{\pm }]n={\displaystyle {\int }_{-a}^a}{\omega }_n\left(x\right)G(x\mathrm{,}\pm a)\text{d}x\\ =\frac{1}{\kappa D\sinh 2\kappa a}{\displaystyle {\int }_{-a}^a}{\omega }_n\left(x\right)\cosh \kappa (a\pm x)\text{d}x\mathrm{,}\\ n=\mathrm{0,1,2}\end{array}$ (56) where $\kappa =\sqrt{\frac{{\Sigma }_r}{D}}$.

Using the continuity conditions of the transverse-averaged flux and the neutron current at the interfaces between adjacent nodes, $f_{+}^k{\Phi }^k\left(a\right)=f_{-}^{k+1}{\Phi }^{k+1}\left(-a\right)$ (57) $J^k\left(-a\right)=J^{k-1}\left(a\right)\mathrm{,}{J}^{k+1}\left(-a\right)=J^k\left(a\right)\mathrm{,}$ (58) the net neutron current coupling equation can be rewritten as $\begin{array}{l}[G_I^{\pm }]n={\displaystyle {\int }_{-a}^a}{\omega }_n\left(x\right)G(x\mathrm{,}\pm a)\text{d}x\\ =\frac{1}{\kappa D\sinh 2\kappa a}{\displaystyle {\int }_{-a}^a}{\omega }_n\left(x\right)\cosh \kappa (a\pm x)\text{d}x\mathrm{,}\\ n=\mathrm{0,1,2}\end{array}$ (59) where $f_{+}^k$ and $f_{-}^{k+1}$ are the discontinuity factors on the right boundary of node k and left boundary of node $k+1$, respectively. Moreover, $R^k=G^k\left(a\mathrm{,}a\right)\mathrm{,}{T}^k=G^k\left(a\mathrm{,}-a\right)$ are called the reflective and transmission factors, respectively.

We substitute Eq. (43) into Eq. (42) and multiply both sides of Eq. (42) by ωm(x). Using $\frac{2y_s\left(x\right)}{2\sqrt{3}{a}^2}$ as the weighting function and integrating over the interval $[-a\mathrm{,}a]$, the following equation for the expansion coefficients of the transverse-averaged flux is obtained: $\begin{array}{ll}\Phi =\hfill & [G](Q-L)+[G+]J(a)+[G-]J(-a)\hfill \\ \hfill & +\frac{D}{2a}\left({\Phi }_3+{\Phi }_6-{\Phi }_1-{\Phi }_2\right)[G_{IS}+]\hfill \\ \hfill & +\frac{D}{2a}\left({\Phi }_3+{\Phi }_6-{\Phi }_4-{\Phi }_5\right)[G_{IS}-]\hfill \end{array}$ (60) where $[G]$ denotes a 3 × 3 matrix with the following elements: $\begin{array}{ll}{[G]}_{mn}\hfill & =\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}2y_s\left(x\right){\omega }_m\left(x\right)\text{d}x{\displaystyle {\int }_{-a}^a}{\omega }_n\left(x_0\right)G\left(x\mathrm{,}{x}_0\right)\text{d}{x}_0\hfill \\ \hfill & =\frac{1}{3a^2}{\displaystyle {\int }_{-a}^a}(2a-|x|){\omega }_m\left(x\right)\text{d}x[{\displaystyle {\int }_{-a}^x}{\omega }_n\left(x_0\right)\frac{\cosh \kappa (a-x)\cosh \kappa \left(a+x_0\right)}{\kappa D\sinh 2\kappa a}\text{d}{x}_0\hfill \\ \hfill & +{\displaystyle {\int }_x^a}{\omega }_n\left(x_0\right)\frac{\cosh \kappa (a+x)\cosh \kappa \left(a-x_0\right)}{\kappa D\sinh 2\kappa a}\text{d}{x}_0]\mathrm{,}m\mathrm{,}n=\mathrm{0,1,2}\hfill \end{array}$ (61) The elements of the column vectors $[G\pm ]$ and $[G_{IS}\pm ]$ are as follows: $\begin{array}{ll}{[G\pm ]}_n\hfill & =\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}2y_s\left(x\right){\omega }_n\left(x\right)G(x\mathrm{,}\pm a)dx\hfill \\ \hfill & =\frac{1}{3a^2\kappa D\sinh 2\kappa a}{\displaystyle {\int }_{-a}^a}(2a-|x|){\omega }_n\left(x\right)\cosh \kappa (a\pm x)\text{d}x\mathrm{,}n=\mathrm{0,1,2}\hfill \end{array}$ (62) $\begin{array}{ll}{[G_{IS}+]}_n\hfill & =\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}2y_s\left(x\right){\omega }_n\left(x\right){\displaystyle {\int }_0^a}\frac{G(x\mathrm{,}{x}_0)}{2a-x_0}\text{d}{x}_0\hfill \\ \hfill & =\frac{1}{3a^2\kappa D\sinh 2\kappa a}{\displaystyle {\int }_0^a}\frac{1}{2a-x_0}d{x}_0[\cosh \kappa \left(a-x_0\right){\displaystyle {\int }_{-a}^0}(2a+x){\omega }_n\left(x\right)\cosh \kappa (a+x)dx\hfill \\ \hfill & +\cosh \kappa \left(a-x_0\right){\displaystyle {\int }_0^{x_0}}(2a-x){\omega }_n\left(x\right)\cosh \kappa (a+x)\text{d}x\hfill \\ \hfill & +\cosh \kappa \left(a+x_0\right){\displaystyle {\int }_{x_0}^a}(2a-x){\omega }_n\left(x\right)\cosh \kappa (a-x)\text{d}x]\mathrm{,}n=\mathrm{0,1,2}\hfill \end{array}$ (63) $\begin{array}{ll}{\left[G_{IS}-\right]}_n\hfill & =\frac{1}{2\sqrt{3}{a}^2}{\displaystyle {\int }_{-a}^a}2y_s\left(x\right){\omega }_n\left(x\right){\displaystyle {\int }_{-a}^0}\frac{G\left(x\mathrm{,}{x}_0\right)}{2a+x_0}\text{d}{x}_0\hfill \\ \hfill & =\frac{1}{3a^2\kappa D\sinh 2\kappa a}{\displaystyle {\int }_{-a}^0}\frac{1}{2a+x_0}d{x}_0[\cosh \kappa \left(a-x_0\right){\displaystyle {\int }_{-a}^{x_0}}(2a+x){\omega }_n\left(x\right)\cosh \kappa (a+x)\text{d}x\hfill \\ \hfill & +\cosh \kappa \left(a+x_0\right){\displaystyle {\int }_{x_0}^0}(2a+x){\omega }_n\left(x\right)\cosh \kappa (a-x)\text{d}x\hfill \\ \hfill & +\cosh \kappa \left(a+x_0\right){\displaystyle {\int }_0^a}(2a-x){\omega }_n\left(x\right)\cosh \kappa (a-x)\text{d}x]\mathrm{,}n=\mathrm{0,1,2}\hfill \end{array}$ (64) By contrast, the corner point flux values of a hexagonal node are included in Eqs. (55) and (60), which are also unknown quantities. For the interior nodes away from the reactor core boundaries, we approximate the corner point flux in terms of the node-averaged fluxes, surface-averaged fluxes, and diffusion coefficients in three adjacent nodes, as follows [21]: $\begin{array}{ll}{\Phi }_0=\hfill & [2\left(D_1+D_2\right){\Phi }_{1^{\prime }}+2\left(D_2+D_3\right){\Phi }_{2^{\prime }}+2\left(D_3+D_1\right){\Phi }_{3^{\prime }}\hfill \\ \hfill & -{\displaystyle \sum _{i=1}^3}{D}_i{\Phi }_i]/3\left(D_1+D_2+D_3\right)\hfill \end{array}$ (65) where ${\Phi }_i$ and ${\Phi }_{i^{\prime }}\left(i=\mathrm{1,2,3}\right)$ are the node- and surface-averaged fluxes, respectively.