JakobCHR.com

Quick Navigation:
 

Personal:

 Go to Home
 MS Research
 PhD Research
 Curriculum Vitae

General:

 Linux

Soon to come:

 Matlab
 On-line Stores
 Cycling
 Medicine & Health
 LaTeX
 OOP & C++
 Sony PCM-R500 DAT

Discretization of the multi-group diffusion equations

The multi-group diffusion equations are three dimensional in nature but because of constraints on time regarding this project I was forced to limit the work to the 1-D case. In a large thermal reactor we are able to neglect the flux variation in the radial direction and treat the flux as if it were only varying in the axial z-direction. This is not a bad assumption in the BWR case since the radial peaking factor is fairly low due to greater steam voids in the center fuel assemblies. With this simplification we can rewrite the multi-group diffusion equations (1.10) as

$$\begin{eqnarray*} \lefteqn{-\frac{d}{dz} [ D^g(z) \frac{d\phi^g(z)}{dz} ] + \S... ...f^{g'}(z) \phi^{g'}(z) \right), \quad\quad\quad g=1,2,\ldots,G \end{eqnarray*}$$

$\parbox{1.50cm}{\begin{eqnarray} \end{eqnarray}}$

In order to solve the eigenvalue problem given by the multi-group diffusion equations (1.16) on a computer it is necessary to discretize the governing equations. The first step is to superimpose a grid on the z-interval of interest. We make use of a nonuniform grid, ie with the steplength h_k depending on the position.

We now discretize (1.16) by an integration method, ie we integrate the equation over a control volume as indicated by Figure 1.2. In the discretized case we are seeking the solution only at the discrete points z_k, as indicated in the figure. We assume that we have employed some homogenization procedure on the regions (k-1) and k, ie the intervals [z_k-1;z_k] and [z_k;z_k+1]. The fundamental assumption is now that the physical properties entering in the diffusion equation are constant within each region.

In mathematical terms we assume the following

$$D^g(z) \approx \overline{D}{}_k^g, \quad\quad\quad\quad z \in [z_k;z_{k+1}]$$ (1.12)

$$\Sigma_a^g(z) \approx \overline{\Sigma}{}_{a,k}^g, \quad\quad\quad\quad z \in [z_k;z_{k+1}]$$ (1.13)

$$\Sigma_s^{g \rightarrow g'}(z) \approx \overline{\Sigma}{}_... ...{g \rightarrow g'}, \quad\quad\quad\quad z \in [z_k;z_{k+1}]$$ (1.14)

$$\Sigma_f^g(z) \approx \overline{\Sigma}{}_{f,k}^g, \quad\quad\quad\quad z \in [z_k;z_{k+1}]$$ (1.15)

where the $\;\overline{\rule{0mm}{0.7em}\;\;\;}{}\;$ indicate that we are dealing with a mean value, ie a mean taken over the region^1.13 (section 13.2.3 p. illustrates how the mean value cross sections are evaluated in practice).

We split this control volume integration up in several parts corresponding to the different terms in (1.16).

The first term becomes

$$\int\limits_{V_{CV}} \nabla \cdot [ D^g(z) \nabla\phi^g ] d... ...\limits_{V_{CV_{k-1}}} \nabla \cdot [ D^g(z) \nabla\phi^g ] dV$$ (1.16)

where V_CVk and V_CVk-1 are parts of the control volume which lie within regions k and k-1 respectively.

The volume integrals may be converted to surface integral by means of Gauss' theorem which states the identity

$$\int\limits_V \nabla \cdot \hspace{0.2ex}\underline{F}{}\hs... ...0.15ex} \cdot \hspace{0.2ex}\underline{n}{}\hspace{0.15ex} dA$$ (1.17)

where $\hspace{0.2ex}\underline{F}{}\hspace{0.15ex}$ is a vector function and V is a volume enclosed by a surface S with an outward pointing unit normal vector $\hspace{0.2ex}\underline{n}{}\hspace{0.15ex}$.

Applying (1.22) to the two terms on the RHS of (1.21) reveals that

$$\int\limits_{V_{CV}} \nabla \cdot [ D^g(z) \nabla\phi^g ] d... ...^g \cdot (-\hspace{0.2ex}\underline{k}{}\hspace{0.15ex}) dxdy$$ (1.18)

where $\hspace{0.2ex}\underline{k}{}\hspace{0.15ex}$ is a unit vector pointing in the z direction and S_CVk and S_CVk-1 are parts of the control surface lying in the kth and $(\mbox{k}-1)$th region respectively. We have not written out terms belonging to the sides pointing in a direction parallel to $\hspace{0.2ex}\underline{i}{}\hspace{0.15ex}$ (x-axis) and $\hspace{0.2ex}\underline{j}{}\hspace{0.15ex}$ (y-axis) since we are assuming that $\phi^g$ only varies in the z direction. Consequently the gradient vector has no component parallel to $\hspace{0.2ex}\underline{i}{}\hspace{0.15ex}$ and $\hspace{0.2ex}\underline{j}{}\hspace{0.15ex}$ directions and the results of the scalar products analogous to the ones for the $\pm\mbox{z}$ direction are all zero. It is important to note that because of the material interface lying at z=z_k care is necessary in the evaluation of the two integrals on the right hand side of (1.23). Evaluating the right hand side integrals we can write

$$\begin{eqnarray*} \int\limits_{S_{CV}} D^g(z) \nabla\phi^g \cdot \hspace{0.2ex}... ...eft.[D^g(z)\nabla\phi^g]\right\vert _{z=z_k^-} \right) \right] \end{eqnarray*}$$

$\parbox{1.50cm}{\begin{eqnarray} \end{eqnarray}}$

where the superscripts + and - indicate right-hand limits and left-hand limits respectively.

Now, by using the continuity condition (1.11) and the fact that the group flux only varies with z, we are able to write the final expression for the volume integral:

$$\int\limits_{V_{CV}} \nabla \cdot [ D^g(z) \nabla\phi^g ] d... ...t.\frac{d\phi^g}{dz}\right\vert _{z=z_{k-\frac{1}{2}}} \right]$$ (1.19)

Returning to (1.16), integration of the second term leads to

$$\begin{eqnarray*} \int\limits_{V_{CV}} \Sigma_a^g(z)\phi^g(z) dV &=& \Delta x \... ...^g \int\limits_{z_k}^{z_{k+\frac{1}{2}}} \phi^g(z) dz \right] \end{eqnarray*}$$

$\parbox{1.50cm}{\begin{eqnarray} \end{eqnarray}}$

where we have used the homogenized (mean value) cross sections which per definition are assumed to be constant inside every region (see (1.18)).

Analogously for the next three terms we may write the resulting expressions as

$$\int\limits_{V_{CV}} \sum\limits_{g'=1}^{G} \Sigma_s^{g \ri... ...\} \int\limits_{z_k}^{z_{k+\frac{1}{2}}} \phi^g(z) dz \right]$$ (1.20)

$$\int\limits_{V_{CV}} \sum\limits_{g'=1}^{G} \Sigma_s^{g' \r... ...its_{z_k}^{z_{k+\frac{1}{2}}} \phi^{g'}(z) dz \right\} \right]$$ (1.21)

$$\begin{eqnarray*} \int\limits_{V_{CV}} \frac{\chi^g}{\lambda_0} \sum\limits_{g... ...its_{z_k}^{z_{k+\frac{1}{2}}} \phi^{g'}(z) dz \right\} \right] \end{eqnarray*}$$

$\parbox{1.50cm}{\begin{eqnarray} \end{eqnarray}}$

It may be useful to get an overview over the form of the total equation. Assuming unit lengths in x and y directions ( $\Delta x = \Delta y = 1$) and adding the four above mentioned expressions we get

$$\begin{eqnarray*} \lefteqn{ D^g(z_{k+\frac{1}{2}}) \left.\frac{d\phi^g}{dz}\ri... ...z_k}^{z_{k+\frac{1}{2}}} \phi^{g'}(z) dz \right\} \right] = 0 \end{eqnarray*}$$

$\parbox{1.50cm}{\begin{eqnarray} \end{eqnarray}}$

It is now time to introduce suitable approximations in regard firstly to $d\phi^g/dz$ and secondly to the integrals involving $\phi^g(z)$. In this situation a suitable approximation is one that only uses values of $\phi^g$ at grid points since we are only seeking the solution at the discrete locations given by $z_k , k=\{1,2,3,\ldots,K\}$. It is symptomatic for the vast majority of texts treating reactor physics, the author has knowledge of, that they just state the approximations without commenting the origin or more importantly the error of the approximations. It is, however, the author's belief that the best way of introducing approximations is to make a derivation of those and in the way of deriving them making it absolutely clear the error accompanying them. In this text we only consider finite difference approximations (in regard to $d\phi^g/dz$) because of their simplicity.

Let us first derive an approximation for the derivative of $\phi^g$ evaluated at point $z_{k-\frac{1}{2}}$--this may be accomplished by considering the Taylor series of $\phi^g$ about point $z_{k-\frac{1}{2}}$ with steps of h_k-1/2 and -h_k-1/2. These Taylor series may be written (assuming $\phi^g \in C^4 \mbox{ on } ]z_{k-1};z_k[$)

$$\begin{eqnarray*} \phi^g(z_k) \equiv \phi^g(z_{k-\frac{1}{2}}+ \frac{h_{k-1}}{2... ...g}{dz^4}\right\vert _{z_{k-\frac{1}{2}}} + {\cal O}(h_{k-1}^5) \end{eqnarray*}$$

$\parbox{1.50cm}{\begin{eqnarray} \end{eqnarray}}$

and

$$\begin{eqnarray*} \phi^g(z_{k-1}) \equiv \phi^g(z_{k-\frac{1}{2}}- \frac{h_{k-1... ...g}{dz^4}\right\vert _{z_{k-\frac{1}{2}}} + {\cal O}(h_{k-1}^5) \end{eqnarray*}$$

$\parbox{1.50cm}{\begin{eqnarray} \end{eqnarray}}$

It is important to note that since we have assumed material interfaces to lie only on grid points the above mentioned Taylor series are only valid inside a region^1.14.

Subtracting (1.32) from (1.31) and dividing by h_k-1 yields

$$\begin{eqnarray*} \frac{\phi^g(z_k) - \phi^g(z_{k-1})}{h_{k-1}} &=& \left.\fra... ...}{dz^3}\right\vert _{z_{k-\frac{1}{2}}} + {\cal O}(h_{k-1}^4) \end{eqnarray*}$$

$\parbox{1.50cm}{\begin{eqnarray} \end{eqnarray}}$

Isolating the wanted derivative reveals

$$\left.\frac{d\phi^g}{dz}\right\vert _{z_{k-\frac{1}{2}}} = ... ...}{dz^3}\right\vert _{z_{k-\frac{1}{2}}} + {\cal O}(h_{k-1}^4)$$ (1.22)

Similarly,

$$\left.\frac{d\phi^g}{dz}\right\vert _{z_{k+\frac{1}{2}}} = ... ...hi^g}{dz^3}\right\vert _{z_{k+\frac{1}{2}}} + {\cal O}(h_k^4)$$ (1.23)

Having treated a way of approximating the derivatives we now turn to ways of approximating the integrals in (1.30). Using a Taylor series of $\phi^g$ about point z=z_k together with Taylor's theorem we are able to write $\phi^g(z)$ as (assuming $\phi^g \in C^3 \mbox{ on } [z_k-\frac{h_{k-1}}{2};z_k^-]$)

$$\begin{eqnarray*} \phi^g(z) &=& \phi^g(z_k) + (z-z_k)\left.\frac{d\phi^g}{dz}\... ...\quad\quad\quad\quad\quad z \in [z_k-\frac{h_{k-1}}{2};z_k^-] \end{eqnarray*}$$

$\parbox{1.50cm}{\begin{eqnarray} \end{eqnarray}}$

for some $\zeta \in [z_k-\frac{h_{k-1}}{2};z_k^-]$^1.15.

Using (1.36) we are now able to write the first integral appearing in (1.30) as follows:

$$\begin{eqnarray*} \int\limits_{z_{k-\frac{1}{2}}}^{z_k} \phi^g(z) dz &\equiv& ... ...ac{d^2\phi^g}{dz^2}\right\vert _{z_k^-} + {\cal O}(h_{k-1}^4) \end{eqnarray*}$$

$\parbox{1.50cm}{\begin{eqnarray} \end{eqnarray}}$

where $\xi=(z-z_k)$ and we have assumed $\phi^g \in C^3 \mbox{ on } [z_k - \frac{h_{k-1}}{2};z_k^-]$. Actually we have again made use of our fundamental assumption that material interfaces are aligned with grid points because it is this assumption that ensures $\phi^g$ to be a smooth function in the interval. Accordingly we are able to write the integral as

$$\int\limits_{z_{k-\frac{1}{2}}}^{z_k} \phi^g(z) dz \equiv ... ....\frac{d\phi^g}{dz}\right\vert _{z_k^-} + {\cal O}(h_{k-1}^3)$$ (1.24)

Analogously we are able to write the second integral appearing in (1.30) as

$$\int\limits_{z_k}^{z_{k+\frac{1}{2}}} \phi^g(z) dz \equiv ... ...left.\frac{d\phi^g}{dz}\right\vert _{z_k^+} + {\cal O}(h_k^3)$$ (1.25)

where we assume $\phi^g \in C^3 \mbox{ on } [z_k^+;z_k + \frac{h_k}{2}]$. Before making the final choice of approximation we write the final exact expression for the multi-group diffusion equation by using the proposed expressions for the derivatives and integrals in (1.30)

$$\begin{eqnarray*} \lefteqn{ \overline{D}{}_k^g \left\{ \frac{\phi^g(z_{k+1}) -... ...\vert _{z_k^+} + {\cal O}(h_k^3) \right) \right\} \right] = 0 \end{eqnarray*}$$

$\parbox{1.50cm}{\begin{eqnarray} \end{eqnarray}}$

We now approximate (1.40) by truncating the Taylor series after the 1st term, and using the grid depicted in Figure 1.3 we obtain

$$\begin{eqnarray*} \lefteqn{\overline{D}{}_k^g \left( \frac{\phi_{k+1}^g - \phi_... ...= 0, \quad\quad g=\{1,2,\ldots,G\} , k = \{2,3,\ldots,(K-1)\} \end{eqnarray*}$$

$\parbox{1.50cm}{\begin{eqnarray} \end{eqnarray}}$

where we have used the definition

$$\phi_k \;\hbox{$=$\kern-0.68em\raise1.1ex \hbox{$\scriptscriptstyle\triangle$}}\;\phi(z_k)$$ (1.26)

Collecting terms and moving the term involving $\lambda_0$ to the right hand side yields

$$\begin{eqnarray*} \lefteqn{ \frac{1}{h_{k-1}} \overline{D}{}_{k-1}^g \phi_{k-1}... ... & \hspace{15em} k=2,3,\ldots,(K-1) \quad,\quad g=1,2,\ldots,G \end{eqnarray*}$$

$\parbox{1.50cm}{\begin{eqnarray} \end{eqnarray}}$

In regard to boundary conditions we only consider the zero flux condition in this text. Accordingly we may state

$$\phi_1^g = 0 \qquad \wedge \qquad \phi_K^g = 0, \qquad g=1,2,\ldots,G$$ (1.27)

where we assume that points 1 and K lie on the outer surface of the reflectors.

At this point it may be useful one more to state the nature of the system of (G x (K-2)) equations. The task is to find the $\phi$-values and the eigenvalue $\lambda_0$, ie we are dealing with an eigenvalue problem^1.16.

The principal truncation error associated with the approximation (1.43) of the integrated multi-group diffusion equation (1.30), $\widetilde{\mbox{PTE}}_k$^1.17, at point z=z_k associated with the above mentioned approximation (1.43) is

$$\begin{eqnarray*} \widetilde{\mbox{PTE}}_k &\;\hbox{$=$\kern-0.68em\raise1.1ex ... ...z_k^+} \right\} \\ &=& {\cal O}(h_{k-1}^2) + {\cal O}(h_k^2) \end{eqnarray*}$$

$\parbox{1.50cm}{\begin{eqnarray} \end{eqnarray}}$

Note that our difference scheme is efficient in the sense that it ensures a matching accuracy of individual terms in the principal part of the truncation error.

It should, however, be noted that the principal truncation error defined with respect to the original multi-group diffusion equation (1.16), ${\mbox{PTE}}_k$, can be stated as

$${\mbox{PTE}}_k = {\cal O}(h_{k-1} + h_k)$$ (1.28)

since the integration of the multi-group diffusion equation introduces a factor $\frac{1}{2} (h_{k-1} + h_k )$ in the truncation error defined by (1.45). It is ${\mbox{PTE}}_k$ which is important in practice since it in most cases of practical interest is such that the discretization error^1.18 (or global truncation error) associated with the discrete scheme is of the same order as the local truncation error when we define the local truncation error with respect to the original equation.

If we consider the case of an uniform mesh ( $h_k = h, \quad \forall k$) and non-varying physical properties we obtain a local truncation error with respect to the original multi-group diffusion equation (1.16) of size

$${\mbox{PTE}}_k = {\cal O}(h^2)$$ (1.29)

For this case it would be appropriate to consider a better approximation of $\left.\frac{d\phi}{dz}\right\vert _{z_k \pm \frac{1}{2}}$. However, since we in general will have varying physical properties and step lengths we will not go into that discussion.

  Contents   Index  
 
 
 

Top Quality Developed with Danish Brain Power