I am trying to solve the neutron diffusion equation to model neutron flux distribution in a one-dimensional two-group setting using the following iterative scheme. The governing equations of the system are:



where $D$=Diffusion Co-efficient, $\phi$=neutron flux, $\mathcal{E}_{a}$=absorption cross-section, $\mathcal{E}_{s}$=scattering cross-section, $\nu$=neutrons per fission, $\mathcal{E}_{f}$=fission cross-section and $k$=multiplication factor. 1 and 2 refers to fast and thermal neutron groups.

I divided the geometry into 420 mesh elements and discretized the equations using Forward Difference Method (FDM). Putting the whole system into matrix form yields:



Where $\Bigl[C\Bigr]$=Co-efficient Matrix and $h$=mesh element length.

Now, for the $n$th iteration, I am supposed to use $\phi_{1}^{n-1}$, $\phi_{2}^{n-1}$ and $k^{n-1}$ to solve the two above equations and calculate $\phi_{1}^{n}$, $\phi_{2}^{n}$ and $k^{n}$ using the iterative scheme:



$$k^n=k^{n-1}\frac{\int dr(\nu_1\mathcal{E}_{f_1}\phi_1^n+\nu_2\mathcal{E}_{f_2}\phi_2^n)}{\int dr(\nu_1\mathcal{E}_{f_1}\phi_1^{n-1}+\nu_2\mathcal{E}_{f_2}\phi_2^{n-1})}$$

I am confused about how to proceed with the last equation. Do I take the differential length $dr$ as mesh element length $h$ and turn the integration into a summation over the geometry? Or do I actually try some kind of integration? I cannot imagine how an integration can be introduced here given that $\phi_1$ and $\phi_2$ are not known as functions but as matrices.

Please help! And thanks for reading through such a long question!


1 Answer 1


You are correct that you have to turn this into a sum over the geometry. You solved for the flux values at the gridpoints and you have the step size between the gridpoints. From this, your choice of integration scheme is limited to those that use uniformly spaced nodes. You could just do a direct sum and multiply by $h$, but in this case you can actually do a bit better by using the trapezoid rule.

  • $\begingroup$ Thanks, I now have a concrete picture in my mind. I'll try to use the trapezoid rule in my code since it looks like this will produce a better approximation. $\endgroup$ Commented Jun 11, 2019 at 4:36
  • $\begingroup$ I would not recommend blindly using a trapezoid rule, you should use an integral that is consistent with your basis functions. If you are assuming that your flux is linear in each node, then use this information to determine the integral. $\endgroup$ Commented Apr 2, 2021 at 18:51

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .