%% Antes de processar este arquivo LaTeX (LaTeX2e) deve %% verificar que o arquivo TEMA.cls estah no mesmo %% diretorio. O arquivo TEMA.cls pode ser obtido do %% endereco www.sbmac.org.br/tema. \documentclass{TEMA} %\usepackage[brazil]{babel} % para texto em Português \usepackage[english]{babel} % para texto em Inglês \usepackage[latin1]{inputenc} % para acentuação em Português %\input{P-margin.inf} %\usepackage[dvips]{graphics} \usepackage{subfigure} \usepackage{graphicx} \usepackage{epsfig} \newcommand{\B}{{\tt\symbol{92}}} \newcommand{\til}{{\tt\symbol{126}}} \newcommand{\chap}{{\tt\symbol{94}}} \newcommand{\agud}{{\tt\symbol{13}}} \newcommand{\crav}{{\tt\symbol{18}}} \begin{document} %******************************************************** \title {Finite Element Method with Spectral Green's Function in Slab Geometry for Neutron Diffusion in Multiplying Means and One Energy Group\thanks{Authors acknowledge to Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq, BRAZIL) and Fundação de Amparo à Pesquisa do Estado da Bahia (FAPESB, BRAZIL) for the financial support provided.}} \author { R. V. M. ROCHA\thanks{rvmrocha@uesc.br; Corresponding author.}\,, D. S. DOMINGUEZ\thanks{dsdominguez@gmail.com}\,, S. M. IGLESIAS\thanks{smiglesias@uesc.br}\,, Programa de Pós-Graduação em Modelagem Computacional, Departamento de Ciências Exatas e Tecnológicas, UESC - Universidade Estadual de Santa Cruz, 45662-900 Ilhéus, BA, Brasil \\ \\ R. C. BARROS% \thanks{rcbarros@pesquisador.cnpq.br}\,, UERJ - Universidade do Estado do Rio de Janeiro, Instituto Politécnico, 28625-570 Nova Friburgo, RJ, Brasil.} \criartitulo \runningheads {Rocha, Dominguez, Iglesias and Barros}{FEM-SGF for Neutron Diffusion in Multiplying Means} \begin{abstract} {\bf Abstract}. The physical phenomenon of neutrons transport associated with eigenvalue problems appears in the criticality calculations of nuclear reactors and can be treated as a diffusion process. This paper presents a new method to solve eigenvalue problems of neutron diffusion in slab geometry and one energy group. This formulation combines the Finite Element Method, considered an intermediate mesh method, with the Spectral Green's Function Method, which is free of truncation errors, and it is considered a coarse mesh method. The novelty of this formulation is to approach the spatial moments of the neutron flux distribution by the first-order polynomials obtained from the spectral analysis of diffusion equation. The approximations provided by the new formulation allow obtaining accurate results in coarse mesh calculations. To validate the method, we compare the results obtained with the methods described in the literature, specifically the Diamond Difference method. The accuracy and the computational performance of the proposed formulation were characterized by solving benchmarks problems with a high degree of heterogeneity. {\bf Keywords}. Eigenvalue problems, Neutron diffusion equation, Spectral Green's Function. \end{abstract} %******************************************************** \newsec{\label{Introduction} Introduction} Population growth and economic development requires the continuing increase in the capacity of electricity generation. Meet this demand without aggravating environmental problems associated with the greenhouse effect is a major challenge for scientists, economists and politicians today. In this context, nuclear energy appears as one of the generation alternatives that do not emit carbon dioxide. In Brazil, the energy development plan provides, in the next years, an increase of nuclear power generation in the energy matrix \cite{pne2030}. Among the main concerns related to the use of nuclear energy is safety and containment of radioactive materials generated in the core fission reactor which is associated to the monitoring and control of the population of neutrons in the reactor. For this reason, during the design and operation of nuclear plants, we need accurate and efficient numerical methods to assess the neutron flux distribution and the effective multiplication coefficient, respectively linked to eigenfunctions and eigenvalues of the system. This problems are associated with criticality calculations of nuclear reactors and them are known as eigenvalue problems \cite{lamarshBaratta2001}. The physical modeling of the phenomenon of neutrons transport consists in to describe the migration of neutrons in a host medium and the probability of interaction between neutrons and atoms of the medium \cite{bellGlasstone1970}. Many scientists treat the physical phenomenon of neutrons transport as a diffusion process, thus, one of the mathematical models that allow us to describe this physical phenomenon is the neutron diffusion equation (DE). A simple model that provides accurate results for the distribution of neutron flux and the effective multiplication coefficient \cite{duderstadtHamilton1975}. For computational performance reasons it is desirable to obtain accurate results for the distribution of neutron flux and the effective multiplication coefficient in a coarse mesh. Among the methods of coarse mesh, we can highlight the Spectral Green's Function Nodal Methods (SGF) \cite{barrosLarsen1991, dominguez2007}. Which had its origins in the Spectral-Nodal Method with constant approximation for discrete ordinates model in source-fixed problems \cite{barrosLarsen1992}. Since then, several others formulations has been developed, for example the Spectral Nodal Methods to discrete ordinates \cite{snd2003} and diffusion \cite{dominguezHernandezBarros2010} problems in Cartesian geometry for neutron multiplying systems. Due to high algebraic cost and to complexity of the computational algorithm of SGFs methods, the Finite Element Method (FEM) has been used for solve nuclear engineering problems and it has offered good results for fine and intermediate meshes \cite{zienkiewicz1971, stacey2007}. The formulation proposed in this work combines the linear approximation of the FEM with an approach quasi-analytic of the SGF. The combination of these methods, aims to obtain accurate results with high computational performance, preserving the analytical spectral solutions inside the nodes. We verified this new formulation by solving benchmark problems and comparing the obtained results with other methods. In next section, we discussed the mathematical model used for describe the physical phenomenon and the discretization scheme. In the section 3, we present a new hybrid formulation. The numerical results are offered in section 4, where we used a benchmark problem with high degree of heterogeneity to validate the new method. Finally, in section~5, we give the final comments and suggestions for future works. %******************************************************** \newsec{\label{TheMathematicalModel} The Mathematical Model} The deterministic mathematical modeling for the eigenvalue problem of transport of neutrons in this paper is based on DE. We consider a static model in slab geometry, disregarding the energy dependence. This model \cite{lamarshBaratta2001, stacey2007} is a system composed for two equations, the balance equation \begin{equation} \label{Eq:eigenvalue_diffusion} \frac{d J(x)}{d x}+\Sigma_a(x)\phi(x) = \frac{1}{k_{eff}}\nu\Sigma_f(x)\phi(x), \end{equation} and, the Fick's law, \begin{equation} \label{Eq:fick_law} J(x) = -D(x)\frac{d \phi(x)}{d x}, \hspace{10mm} 0 \leq x \leq L, \end{equation} The equation~(\ref{Eq:eigenvalue_diffusion}) represents the mono-energetic diffusion equation, which describes the gain and loss of neutrons, while (\ref{Eq:fick_law}) refers to Fick's law, which describes the neutron's motion, that is, the neutrons migrate from regions of high to low density \cite{duderstadtHamilton1975}. In~(\ref{Eq:eigenvalue_diffusion}) and~(\ref{Eq:fick_law}) $J(x)$ is the neutron current, $\phi(x)$ is the neutron flux, $\Sigma_a(x)$ is the macroscopic absorption cross section, $k_{eff}$ is the effective multiplication coefficient, $\nu$ is the average number of neutrons issued by fission event, $\Sigma_f(x)$ is the macroscopic fission cross section, $L$ is the longitude of the one-dimensional domain. The boundary conditions considered in this work was Albedo boundary conditions \cite{lewisMiller1993}, which are represented by \begin{equation} \label{Eq:boundary_condition} J(0) = -\alpha_L\phi(0), \hspace{10mm} J(L) = \alpha_R\phi(L). \end{equation} the terms $\alpha_{L,R}$ are constants of proportionality of Albedo, that is, in Albedo boundary conditions the incident flux is proportional to emergent flux. %-------------------------------------------------------- \subsection{Spatial Discretization Scheme} \label{TNDE:spatial_discretization} In general, the derivation of the numerical solution for differential equations systems is associated to discretization of the domain. Let us define a spatial mesh with a total length of $L$ cm, with $I$ nodes called $\Omega_i$, $i = 1:I$, each one with $h_i$ of length. The representation of this spatial mesh can be seen in Fig.~\ref{Fig:spatial_mesh}. \begin{figure}[h] \centering % \epsfig{file=spatialMesh.eps,height=10cm,width=12cm} \scalebox{0.25}{\includegraphics{spatialMesh.eps}} \caption{Spatial mesh in the domain of length $L$.} \label{Fig:spatial_mesh} \end{figure} Now, we will write the equations~(\ref{Eq:eigenvalue_diffusion}) and~(\ref{Eq:fick_law}) in discretized form for each node $\Omega_i$ in the domain, \begin{equation} \label{Eq:eigenvalue_diffusion_discretized} \frac{d J(x)}{d x}+\Sigma_{ai}\phi(x) = \frac{1}{k_{eff}}\nu\Sigma_{fi}\phi(x), \hspace{30mm} \end{equation} \begin{equation} \label{Eq:fick_law_discretized} J(x) = -D_i\frac{d \phi(x)}{d x}, \hspace{7mm} i = 1:I, \hspace{7mm} x_{i-1/2} \leq x \leq x_{i+1/2}, \end{equation} the Albedo boundary conditions, equations~(\ref{Eq:boundary_condition}), appear in the following discretized form, \[ \label{Eq:boundary_condition_discretized} J_{1/2} = -\alpha_L\phi_{1/2}, \hspace{10mm} J_{I+1/2} = \alpha_R\phi_{I+1/2}. \] In this paper, the computational modeling requires that we construct balance equations for the moments of zero and first order for each node $\Omega_i$ in the domain. For we obtain this balance equations, we need to apply the operator, \[ \label{Eq:integral_operator_legendre_polynomial} \frac{2l + 1}{h_i} \int_{x_{i-1/2}}^{x_{i+1/2}} \bullet P_l \left( \frac{2(x - x_i)}{h_i} \right) dx, \] in~(\ref{Eq:eigenvalue_diffusion_discretized}) and~(\ref{Eq:fick_law_discretized}). This operator is based in Legendre polynomials \cite{attar2009}, where $P_l$ represents the Legendre polynomial of degree $l$. For the zero order balance equations, $l = 0$, the Legendre polynomial is a constant with value $1$, that is, $P_0(x) = 1$. For first order, $l = 1$, we have a linear polynomial, $P_1(x) = x$. The balance equations for the moment of zero order appear in the following form, \[ \label{Eq:balance_equation_zero_order_1} \frac{1}{h_i} (J_{i+1/2} - J_{i-1/2}) + \Sigma_{ai} \overline{\phi_i} = \frac{1}{k_{eff}}\nu\Sigma_{fi} \overline{\phi_i}, \] \[ \label{Eq:balance_equation_zero_order_2} \overline{J_i} = -\frac{D_i}{h_i} (\phi_{i+1/2} - \phi_{i-1/2}), \] in which, we defined the average flux into node as, \begin{equation} \label{Eq:average_flux} \overline{\phi_i} = \frac{1}{h_i}\int_{x_{i-1/2}}^{x_{i+1/2}} \phi(x) dx, \end{equation} and the average current into node as, \[ \label{Eq:average_current} \overline{J_i} = \frac{1}{h_i}\int_{x_{i-1/2}}^{x_{i+1/2}} J(x) dx. \] For the first order balance equations we have the following representation, \[ \label{Eq:balance_equation_first_order_1} \frac{3}{h_i} (J_{i+1/2} + J_{i-1/2} - 2\overline{J_i}) + \Sigma_{ai} \widehat{\phi_i} = \frac{1}{k_{eff}}\nu\Sigma_{fi} \widehat{\phi_i}, \] \[ \label{Eq:balance_equation_first_order_2} \widehat{J_i} = -\frac{3D_i}{h_i} (\phi_{i+1/2} + \phi_{i-1/2} - 2\overline{\phi_i}), \] where, the first order moment of flux is, \begin{equation} \label{Eq:first_order_flux_moment} \widehat{\phi_i} = \frac{3}{h_i}\int_{x_{i-1/2}}^{x_{i+1/2}} \left( \frac{2(x - x_i)}{h_i} \right) \phi(x) dx, \end{equation} and the first order moment for the current is defined as, \[ \label{Eq:first_order_current_moment} \widehat{J_i} = \frac{3}{h_i}\int_{x_{i-1/2}}^{x_{i+1/2}} \left( \frac{2(x - x_i)}{h_i} \right) J(x) dx. \] %******************************************************** \newsec{\label{TheFEMSGFFormulation} The FEM-SGF Formulation} The new formulation proposed in this paper is a hybrid formulation that combines two methods found in the literature. One of the methods is the FEM, which obtain the solution of the diffusion model for eigenvalue problems by approximation functions. In this work, we choose two linear approximation functions, one for neutron flux and one for neutron current, such functions satisfy the conditions described by integral equations in the domain. The other method that compose the hybrid formulation is the one-dimensional SGF~\cite{barrosLarsen1991}, this method approach the solution of diffusion equation with two spectral para\-meters $\alpha_i$ and $\beta_i$, obtained by spectral analysis of the mathematical model. Thus, the hybrid formulation proposed is based on the combination of linear approximation of the FEM with the approach quasi-analytic of the SGF. The approximation functions of the FEM-SGF appear in the following form, \begin{equation} \label{Eq:hybrid_auxiliar_equation_flux} \phi (x) = \alpha_i \overline{\phi_i} + \frac{2 \beta_i}{h_i} (x - x_i) \widehat{\phi_i}, \end{equation} \begin{equation} \label{Eq:hybrid_auxiliar_equation_current} J (x) = \alpha_i \overline{J_i} + \frac{2 \beta_i}{h_i} (x - x_i) \widehat{J_i}. \end{equation} where, we couple the spectral parameters of the SGF in the linear functions of the FEM. To obtain the discretization of the approximation functions for the FEM-SGF, we need to evaluate this functions in the node edges. For the nodes belonging to first half of domain, we evaluate the functions~(\ref{Eq:hybrid_auxiliar_equation_flux}) and~(\ref{Eq:hybrid_auxiliar_equation_current}) in the right edge of the node, $x_{i+1/2}$, this functions discretized appear as, \begin{equation} \label{Eq:left_domain_hybrid_auxiliar_equation_flux_discretized} \phi_{i+1/2} = \alpha_i \overline{\phi_i} + \beta_i \widehat{\phi_i}, \end{equation} \begin{equation} \label{Eq:left_domain_hybrid_auxiliar_equation_current_discretized} J_{i+1/2} = \alpha_i \overline{J_i} + \beta_i \widehat{J_i}. \end{equation} For the nodes in the second half of domain, the functions~(\ref{Eq:hybrid_auxiliar_equation_flux}) and~(\ref{Eq:hybrid_auxiliar_equation_current}) are evaluated in the left edge of the node and they are represented for, \begin{equation} \label{Eq:right_domain_hybrid_auxiliar_equation_flux_discretized} \phi_{i-1/2} = \alpha_i \overline{\phi_i} - \beta_i \widehat{\phi_i}, \end{equation} \begin{equation} \label{Eq:right_domain_hybrid_auxiliar_equation_current_discretized} J_{i-1/2} = \alpha_i \overline{J_i} - \beta_i \widehat{J_i}. \end{equation} In the case in which the mesh consider an odd number of nodes, we must to equal the average flux, equations~(\ref{Eq:left_domain_hybrid_auxiliar_equation_flux_discretized}) and~(\ref{Eq:right_domain_hybrid_auxiliar_equation_flux_discretized}), and average current, equations~(\ref{Eq:left_domain_hybrid_auxiliar_equation_current_discretized}) and~(\ref{Eq:right_domain_hybrid_auxiliar_equation_current_discretized}), inside central node to ensure the symmetry conditions in the domain. The resultant equations appear as, \[ \label{Eq:central_node_hybrid_auxiliar_equation_flux_discretized} 2 \beta_i \widehat{\phi_i} - \phi_{i+1/2} + \phi_{i-1/2} = 0, \] \[ \label{Eq:central_node_hybrid_auxiliar_equation_current_discretized} 2 \beta_i \widehat{J_i} - J_{i+1/2} + J_{i-1/2} = 0. \] Now, we calculate the spectral parameters, $\alpha_i$ and $\beta_i$, of the approximation functions of the FEM-SGF. For that, we need to obtain general analytic solutions for diffusion equation, which can be obtained by the spectral analysis technical proposed by Case and Zweifel~\cite{caseZweifel1967}. The general analytic solutions of the DE must be preserved inside each spatial node $\Omega_i$, and appears in the following form, \begin{eqnarray} \label{Eq:flux_general_analytical_solution} \phi(x) = k_1 e^{x/\mu_i} + k_2 e^{-x/\mu_i}, \end{eqnarray} \[ \label{Eq:current_general_analytical_solution} J(x) = k_1 \frac{-D_i}{\mu_i} e^{x/\mu_i} + k_2 \frac{D_i}{\mu_i} e^{-x/\mu_i}, \] where $k_1$ and $k_2$ are arbitrary constants, and $\mu_i$ are the eigenvalues of the problem. Next, we choose, arbitrarily, any of the functions of the general analytic solutions, because both functions carry the arbitrary constants $k_1$ and $k_2$, by the way, we choose the function $\phi(x)$, equation~(\ref{Eq:flux_general_analytical_solution}). The next step is to obtain the flux $\phi_{i \pm 1/2}$, by evaluating (\ref{Eq:flux_general_analytical_solution}) in the node edges ($x_{i \pm 1/2}$), which appears as \begin{eqnarray} \label{Eq:flux_general_analytical_solution_discretized} \phi_{i \pm 1/2} = k_1 e^{x_i/\mu_i} e^{\pm h_i/2\mu_i} + k_2 e^{-x_i/\mu_i} e^{\mp h_i/2\mu_i}, \end{eqnarray} next, we calculate the average flux $\overline{\phi_i}$, by substituting~(\ref{Eq:flux_general_analytical_solution}) in the definition~(\ref{Eq:average_flux}), the result appears as \begin{eqnarray} \label{Eq:general_average_flux_discretized} \overline{\phi_i} = k_1 \frac{\mu_i e^{x_i/\mu_i}}{h_i} \left( e^{h_i/2\mu_i} - e^{-h_i/2\mu_i} \right) + k_2 \frac{\mu_i e^{-x_i/\mu_i}}{h_i} \left( e^{h_i/2\mu_i} - e^{-h_i/2\mu_i} \right), \end{eqnarray} and last, calculate the first order moment of neutron flux $\widehat{\phi_i}$, by substituting~(\ref{Eq:flux_general_analytical_solution}) in the definition~(\ref{Eq:first_order_flux_moment}), the result appears in the following form \begin{eqnarray} \label{Eq:general_first_order_flux_moment_discretized} \widehat{\phi_i} = k_1 \frac{3 \mu_i e^{x_i/\mu_i}}{h_i} \left( e^{h_i/2\mu_i} + e^{-h_i/2\mu_i} \right) - k_1 \frac{6 \mu_i^2 e^{x_i/\mu_i}}{h_i^2} \left( e^{h_i/2\mu_i} - e^{-h_i/2\mu_i} \right) \\ \nonumber - k_2 \frac{3 \mu_i e^{-x_i/\mu_i}}{h_i} \left( e^{h_i/2\mu_i} + e^{-h_i/2\mu_i} \right) + k_2 \frac{6 \mu_i^2 e^{-x_i/\mu_i}}{h_i^2} \left( e^{h_i/2\mu_i} - e^{-h_i/2\mu_i} \right). \end{eqnarray} Now, we can substitute the expressions (\ref{Eq:flux_general_analytical_solution_discretized}), (\ref{Eq:general_average_flux_discretized}), (\ref{Eq:general_first_order_flux_moment_discretized}) obtained respectively for $\phi_{i \pm 1/2}$, $\overline{\phi_i}$ and $\widehat{\phi_i}$, in~(\ref{Eq:left_domain_hybrid_auxiliar_equation_flux_discretized}) or~(\ref{Eq:right_domain_hybrid_auxiliar_equation_flux_discretized}) and organize these equations in the following form \begin{eqnarray} \label{Eq:general_hybrid_auxiliar_equation_flux} k_1 (A - C) + k_2 (B - D) = 0, \end{eqnarray} where, \[ \label{Eq:a} A = e^{x_i/\mu_i} e^{\pm h_i/2\mu_i}, \] \[ \label{Eq:b} B = e^{- x_i/\mu_i} e^{\mp h_i/2\mu_i}, \] \[ \label{Eq:c} C = \alpha_i \frac{\mu_i e^{x_i/\mu_i}}{h_i} \left( e^{h_i/2\mu_i} - e^{-h_i/2\mu_i} \right) + \beta_i \frac{3 \mu_i e^{x_i/\mu_i}}{h_i} \left( e^{h_i/2\mu_i} + e^{-h_i/2\mu_i} \right) \] \[ - \beta_i \frac{6 \mu_i^2 e^{x_i/\mu_i}}{h_i^2} \left( e^{h_i/2\mu_i} - e^{-h_i/2\mu_i} \right), \] \[ \label{Eq:d} D = \alpha_i \frac{\mu_i e^{-x_i/\mu_i}}{h_i} \left( e^{h_i/2\mu_i} - e^{-h_i/2\mu_i} \right) - \beta_i \frac{3 \mu_i e^{-x_i/\mu_i}}{h_i} \left( e^{h_i/2\mu_i} + e^{-h_i/2\mu_i} \right) \] \[ + \beta_i \frac{6 \mu_i^2 e^{-x_i/\mu_i}}{h_i^2} \left( e^{h_i/2\mu_i} - e^{-h_i/2\mu_i} \right). \] But, we want the general analytic solutions be preserved for all value of the $k_1$ and $k_2$, thus, (\ref{Eq:general_hybrid_auxiliar_equation_flux}) will be satisfied if $A = C$ and $B = D$. These equalities provide a linear system that lets us calculate the parameters $\alpha_i$ and $\beta_i$. Depending on the real or imaginary feature of the eigenvalue problem $\mu_i$, we can have two alternatives, \begin{flushleft} \begin{description} \item[(\textit{i}).] If $\Sigma_{ai} > \displaystyle\frac{1}{k_{eff}}\nu\Sigma_{fi}$ \\ \begin{eqnarray} \label{Eq:real_parameters_alfa_beta} \alpha_i = \gamma_i \coth{\left( \gamma_i \right)}, \hspace{10mm} \displaystyle \beta_i = \frac{\sinh{\left( \gamma_i \right)}}{\displaystyle\frac{3 \cosh{\left( \gamma_i \right)}}{\gamma_i} - \frac{3 \sinh{\left( \gamma_i \right)}}{\gamma_i^2}}; \end{eqnarray} \item[(\textit{ii}).] If $\Sigma_{ai} < \displaystyle\frac{1}{k_{eff}}\nu\Sigma_{fi}$ \\ \begin{eqnarray} \label{Eq:imaginary_parameters_alfa_beta} \alpha_i = \gamma_i \cot{\left( \gamma_i \right)}, \hspace{10mm} \displaystyle \beta_i = - \frac{\sin{\left( \gamma_i \right)}}{\displaystyle\frac{3 \cos{\left( \gamma_i \right)}}{\gamma_i} - \frac{3 \sin{\left( \gamma_i \right)}}{\gamma_i^2}}; \end{eqnarray} \end{description} \end{flushleft} where $\gamma_i = h_i / 2\mu_i$. The equations~(\ref{Eq:real_parameters_alfa_beta}) and~(\ref{Eq:imaginary_parameters_alfa_beta}) ensure that the general analytic solutions be preserved inside each spatial node $\Omega_i$. The algebraic linear system obtained by FEM-SGF is composed by ($6I + 2$) variables and the same number of discretized equations. This system can be reduced to an equivalent system of order ($2I + 2$), what results in the elimination of the average flux and current ($\overline{\phi_i}$ and $\overline{J_i}$), and of the first order moments of the flux and current ($\widehat{\phi_i}$~and~$\widehat{J_i}$). Then, this reduced algebraic system would be solve using Gaussian Elimination with Scaled Partial Pivoting \cite{burdenFaires2011} to obtain neutron flux distribution, and combined with Power Method \cite{burdenFaires2011} to estimate an effective multiplication coefficient. %******************************************************** \newsec{\label{ResultsDiscussion} Results and Discussion} In order to quantify the accuracy and computational performance of the new method developed we compare the proposed formulation with conventional methods, specifi\-cally the Diamond Difference method (DD) \cite{lewisMiller1993}, which is considered a method of fine mesh and it has a low computational performance. The benchmark chosen for the computational experiments has a domain with $150$~$cm$ of length, divided in $6$~regions, whose materials and geometric parameters are shown in Table~\ref{Tab:materials_parameters}. The boundary conditions have $\alpha_L = 0$ for the left boundary, and $\alpha_R = 0.5$ for the right boundary. The convergence criteria used for stop iterative process was $10^{-4}\%$ in the effective multiplication coefficient, and the nominal power considered in the reactor is $10$~$MW$. \begin{table}[h] \caption{Materials and geometric parameters of benchmark problem.} \label{Tab:materials_parameters} \begin{center} \begin{tabular}{|c|c|c|c|c|} \hline & $D$ ($cm^{-1}$) & $\Sigma_a$ ($cm^{-1}$) & $\nu\Sigma_f$ ($cm^{-1}$) & $length$ $(cm)$ \\ \hline \hline \textbf{Region 1} & $1$ & $0,24$ & $0,2733$ & $25$ \\ \hline \textbf{Region 2} & $1,\overline{3}$ & $0,22$ & $0,2565$ & $25$ \\ \hline \textbf{Region 3} & $1$ & $0,44$ & $0,17$ & $30$ \\ \hline \textbf{Region 4} & $1,\overline{3}$ & $0,1726$ & $0,23$ & $15$ \\ \hline \textbf{Region 5} & $2,\overline{7}$ & $0,046$ & $0,015$ & $30$ \\ \hline \textbf{Region 6} & $1$ & $0,3$ & $0,3526$ & $25$ \\ \hline \end{tabular} \end{center} \end{table} The numerical results for neutron flux distribution can be seen in Figs.~\ref{Fig:flux} and~\ref{Fig:relative_error_flux}. Specifically in Fig.~\ref{Fig:flux}, we plot the behavior of the neutron flux in the domain. The reference solution was generated using a finer mesh composed by $960$~nodes with the FEM method. Note that there are strong gradients of the neutron flux in some regions. This situation is so complicated to be solved by the most of the numerical methods. \begin{figure}[h] \centering % \epsfig{file=flux.eps,height=10cm,width=12cm} \scalebox{0.4}{\includegraphics{flux.eps}} \caption{Neutron flux distribution in the domain for DD, FEM and FEM-SGF methods using 30 nodes in spatial grid.} \label{Fig:flux} \end{figure} Fig.~\ref{Fig:flux} shows us that the FEM-SGF method, with a coarse mesh~(30~nodes), keep the neutron flux close to reference solution, while the FEM, with the same mesh, moves away from the reference solutions in the regions with strong gradients. On the other hand, the DD method generates meaningless results for a coarse mesh calculations. The differences between numerical results obtained by the three methods can be more noticeable in Fig.~\ref{Fig:relative_error_flux}, which shows the relative deviations of the methods with respect to reference solution in all domain. In Fig.~\ref{Fig:relative_error_flux}, we can see that the strongest variation of the neutron flux occurs in the region between $x = 50$~and~$x = 75$, where the numerical methods find the major difficulty in to approximate the neutron flux. The FEM-SGF presents the smallest relative deviations, $10^{-1}\%$ appro\-ximately, while in the FEM the relative errors vary between $1\%$~to~$10^{2}\%$. The relative errors associated to DD are so higher, $10^{2}\%$ approximately. \begin{figure}[h] \centering % \epsfig{file=relative_error_flux.eps,height=10cm,width=12cm} \scalebox{0.4}{\includegraphics{relative_error_flux.eps}} \caption{Relative deviations in neutron flux distribution for DD, FEM and FEM-SGF methods in coarse mesh solution.} \label{Fig:relative_error_flux} \end{figure} For calculate the power density in the reactor, we choose the region between $x = 60$ and $x = 100$ $cm$, which presents a strong variation in the neutron flux and a high density of neutrons. Table~\ref{Tab:power_density} brings the results for the power density, which shows the power in $MW$ for each method and its respective relative deviations. \newcommand{\cc}{\centering} \newcommand{\rr}{\raggedright} \newcommand{\tn}{\tabularnewline} \renewcommand{\arraystretch}{1.5} \begin{table}[h] \caption{Numerical results for power density in the region ($60 \leq x \leq 100$) for DD, FEM and FEM-SGF methods in coarse mesh solution.} %title of the table \label{Tab:power_density} \begin{center} \begin{tabular}{|c|c|c|} \hline \cc \textbf{Method (nodes)} &\cc \textbf{Power (MW)} &\cc \textbf{Rel. deviation (\%)} \tn \hline \hline \rr \textbf{DD (30)} &\cc 6.208733E-02 &\cc 9.769013E+01 \tn \hline \rr \textbf{FEM (30)} &\cc 2.523533E+00 &\cc 6.115466E+00 \tn \hline \rr \textbf{FEM-SGF (30)} &\cc 2.692115E+00 &\cc 1.563635E-01 \tn \hline \rr \textbf{Reference: FEM (960)} &\cc 2.687912E+00 &\cc .. \tn \hline \end{tabular} \end{center} \end{table} The power density calculated with the distribution of neutron flux obtained by FEM-SGF presented a good result when we compare to reference solution, its relative deviation is about~$0.1\%$, while the FEM shown a relative deviation of~$6\%$, approximately. The DD method was more sensitive to coarse mesh and the strong variation of the flux in the region, it presented a relative deviation of~$97.6\%$, approximately. The numerical results for the effective multiplication coefficient ($k_{eff}$) appears in Table~\ref{Tab:effective_multiplication_coefficient}, which shows the value of the $k_{eff}$ for each method and its respective relative deviations. The value of the $k_{eff}$ obtained by FEM-SGF have the smallest relative deviation, $10^{-5}\%$, approximately. The results for the conventional methods, FEM and DD, suffer a bigger influence of the spatial grid, which present relative errors of $10^{-3}\%$~and~$10^{-2}\%$, respectively. \begin{table}[h] \caption{Numerical results for effective multiplication coefficient for DD, FEM and FEM-SGF methods in coarse mesh solution.} %title of the table \label{Tab:effective_multiplication_coefficient} \begin{center} \begin{tabular}{|c|c|c|} \hline \cc \textbf{Method (nodes)} &\cc $k_{eff}$ &\cc \textbf{Relative Error (\%)} \tn \hline \hline \rr \textbf{DD (30)} &\cc 1.140051E+00 &\cc 3.218296E-02 \tn \hline \rr \textbf{FEM (30)} &\cc 1.140309E+00 &\cc 9.573804E-03 \tn \hline \rr \textbf{FEM-SGF (30)} &\cc 1.140419E+00 &\cc 7.634042E-05 \tn \hline \rr \textbf{Reference: FEM (960)} &\cc 1.140418E+00 &\cc .. \tn \hline \end{tabular} \end{center} \end{table} To characterize the computational performance of the FEM-SGF, we consider the dominant operation for the implementations of the numerical formulations analyzed. Dominant operation is associated at the solution of the algebraic linear system generated by the discretization procedure in each method. Due the specific characteristic of the FEM-SGF, where the matrix coefficients depend on problem eigenvalue and must be update in each iteration, the implementation of the hybrid formulation use the direct method of Gaussian Elimination with Scaled Partial Pivoting \cite{burdenFaires2011} for solve the algebraic linear system. The computational cost of this method, in asymptotic notation, is~$O(I^3)$, thus, the computational cost of the FEM-SGF is $O(M) O(I^3)$, which $M$ is iterations number of power method to converge. On the other hand, the computational cost of the conventional methods, FEM and DD, in asymptotic notation, is~$O(I^3)$~and~$O(M)O(I)$, respectively. It is noteworthy that the DD generate an algebraic linear system with~($I + 1$) variables and equations, and the order of algebraic linear system generated by FEM is equivalent to FEM-SGF. Fig.~\ref{Fig:time_precision} shows us the behavior of the computational performance of the formulations studied in this work. In this figure, we measure the processing time according to precision in calculations of the power density into region between $x = 60$~and~$x = 100$~$cm$. We can see that the FEM-SGF obtained a good performance, with a processing time lower than the conventional methods, when the precision was~$0.1\%$. However, when the precision increase, i.e., the number of nodes increase, the performance of the new formulation decreases in relation to conventional methods. This performance loss is due to the computational cost of the Gaussian Elimination method for solve the algebraic linear system when the precision is increased. \begin{figure}[h] \centering % \epsfig{file=time_precision.eps,height=10cm,width=12cm} \scalebox{0.4}{\includegraphics{time_precision.eps}} \caption{Processing time versus on relative deviation in calculating power for DD, FEM and FEM-SGF methods.} \label{Fig:time_precision} \end{figure} %******************************************************** \newsec{\label{ConclusionsFutureWorks} Conclusions and Future Works} We presented in this work a new formulation for solve the one-dimensional neutron diffusion equation in multiplying means and one energy group. The character hybrid of this formulation is in the combination of the linear functions of the FEM with spectral parameters of the SGF methods which let us preserve the general analytic solutions inside node. This formulation let us to obtain numerical results with accuracy in a coarse mesh, what suggest a high computational performance. The FEM-SGF generated a satisfactory numerical results for the neutron flux distribution and the effective multiplication coefficient in a coarse spatial grid. It was superior to conventional methods, FEM and DD, relative to accuracy and computational performance. The computational experiments with coarse meshes showed that in regions with strong variations in the neutron flux, the FEM-SGF keep the numerical results close to reference solution. This is in contrast to the FEM and DD, which generate results so far of reference solution. In relation to computational performance, the hybrid formulation achieved good results for a coarse mesh. We suggest, in future works, to reduce the algebraic linear system for~($I + 1$) order and to explore its numerical resolution with iterative methods. We suggest too, to utilize the hybrid formulation for solve one-dimensional problems for two energy groups and in multi-dimensional geometries. \vspace*{.3cm} \noindent {\bf {\large Acknowledgments}}\ \\ Authors acknowledge the State University of Santa Cruz (UESC, BRAZIL) for the support to this paper and to Núcleo de Biologia Computacional e Gestão de Informações Biotecnológicas (NBCGIB, BRAZIL) for the infrastructure offered for the development of computational experiments. \begin{abstract} {\bf Resumo}. A utilização de centrais nucleares para geração de energia elétrica é uma das principais alternativas para atender a demanda energética nos próximos anos sem agravar os problemas ambientais associados ao aquecimento global. O uso de métodos e técnicas de simulação que estimem a população de nêutrons no núcleo é fundamental para garantir a operação segura e confiável do reator. Neste trabalho apresentamos uma nova formulação híbrida para resolver problemas de autovalor de difusão de nêutrons em domínios unidimensionais e aproximação de uma velocidade. Esta formulação combina a simplicidade do Método de Elementos Finitos com uma aproximação Espectro-Nodal e permite obter resultados precisos em cálculos de malha grossa. A precisão e o desempenho computacional do método são caracterizados mediante a solução de um problema modelo com alto grau de heterogeneidade. \end{abstract} \begin{thebibliography}{8} \bibitem{attar2009} R. E. Attar, ``Legendre Polynomials and Functions'', CreateSpace, 2009. \bibitem{barrosLarsen1991} R. C. Barros, E. W. Larsen, A Numerical Method for Multigroup Slab-Geometry Discrete Ordinates Problems with No Spatial Truncation Error, {\em Transport Theory and Statistical Physics}, {\bf 20}, pp. 441, (1991). \bibitem{barrosLarsen1992} R. C. Barros, E. W. Larsen, A spectral nodal method for one-group X,Y-geometry discrete ordinates problems, {\em Nuclear Science and Engineering}, {\bf 111}, No. 1 (1992), 34--45. \bibitem{snd2003} R. C. Barros, H. A. Filho, E. T. V. Orellana, F. C. Silva, N. Couto, D. S. Dominguez, C. R. G. Hernandez, The Application of Spectral Nodal Methods to Discrete Ordinates and Diffusion Problems in Cartesina Geometry for Neutron Multiplying Systems, {\em Progress in Nuclear Energy}, {\bf 42}, (2003), 385-426. \bibitem{bellGlasstone1970} G. I. Bell, S. Glasstone, ``Nuclear reactor theory'', Van Nostrand Reinhold Co., 1970. \bibitem{burdenFaires2011} R. L. Burden, D. J. Faires, ``Numerical Analysis'', Ninth Edition, BROOKS/COLE - CENGAGE Learning, 2011. \bibitem{caseZweifel1967} K. M. Case, P. F. Zweifel, ``Linear transport theory'', Addison-Wesley Pub. Co., 1967. \bibitem{dominguez2007} D. S. Dominguez, R. C. Barros, The spectral Green's function linear-nodal method for one-speed X,Y-geometry discrete ordinates deep penetration problems, {\em Annals of Nuclear Energy}, {\bf 34}, pp. 958, (2007). \bibitem{dominguezHernandezBarros2010} D. S. Dominguez, C. R. G. Hernandez, R. C. Barros, Spectral nodal method for numerically solving two-energy group X,Y geometry neutron diffusion eigenvalue problems, {\em International Journal of Nuclear Energy, Science and Technology (Print)}, {\bf 5}, (2010), 66. \bibitem{duderstadtHamilton1975} J. J. Duderstadt, L. J. Hamilton, ``Nuclear Reactor Analysis'', John Wiley \& Sons Inc, 1975. \bibitem{pne2030} Empresa de Pesquisa Energética, ``Plano Nacional de Energia 2030'', Ministério de Minas e Energia, Rio de Janeiro, Brasil, 2007. \bibitem{lamarshBaratta2001} J. R. Lamarsh, A. J. Baratta, ``Introduction Nuclear Engineering'', Prentice Hall, 2001. \bibitem{lewisMiller1993} E. E. Lewis, W. F. Miller Jr, ``Computational Methods of Neutron Transport'', American Nuclear Society, Illinois, USA, 1993. \bibitem{stacey2007} W. M. Stacey, ``Nuclear Reactor Physics'', Wiley-VCH, 2007. \bibitem{zienkiewicz1971} O. C. Zienkiewicz, ``The Finite Element Methods in Engineering Science'', McGraw-Hill, 1971. \end{thebibliography} \end{document} \newpage $ \ \ $ \thispagestyle{myheadings} \markboth{ }{ }