# How to make use of this iterative scheme for solving the neutron diffusion equation?

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:

$$-D_1\nabla^2\phi_1+(\mathcal{E}_{a_1}+\mathcal{E}_{s12})\phi_1=\frac{1}{k}[\nu_1\mathcal{E}_{f_1}\phi_1+\nu_2\mathcal{E}_{f_2}\phi_2]+\mathcal{E}_{s21}\phi_2$$

$$-D_2\nabla^2\phi_2+(\mathcal{E}_{a_2}+\mathcal{E}_{s21})\phi_2=\mathcal{E}_{s12}\phi_1$$

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:

$$\Bigl[C_1\Bigr]\Bigl[\phi_1\Bigr]=\frac{1}{k}\Bigl[h\Bigr]\Bigl[\nu_1\mathcal{E}_{f_1}\Bigr]\Bigl[\phi_1\Bigr]+\frac{1}{k}\Bigl[h\Bigr]\Bigl[\nu_2\mathcal{E}_{f_2}\Bigr]\Bigl[\phi_2\Bigr]+\Bigl[h\Bigr]\Bigl[\mathcal{E}_{s21}\Bigr]\Bigl[\phi_2\Bigr]$$

$$\Bigl[C_2\Bigr]\Bigl[\phi_2\Bigr]=\Bigl[h\Bigr]\Bigl[\mathcal{E}_{s12}\Bigr]\Bigl[\phi_1\Bigr]$$

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:

$$\Bigl[C_1\Bigr]\Bigl[\phi_1\Bigr]^n=\frac{1}{k}\Bigl[h\Bigr]\Bigl[\nu_1\mathcal{E}_{f_1}\Bigr]\Bigl[\phi_1\Bigr]^{n-1}+\frac{1}{k}\Bigl[h\Bigr]\Bigl[\nu_2\mathcal{E}_{f_2}\Bigr]\Bigl[\phi_2\Bigr]^{n-1}+\Bigl[h\Bigr]\Bigl[\mathcal{E}_{s21}\Bigr]\Bigl[\phi_2\Bigr]^{n-1}$$

$$\Bigl[C_2\Bigr]\Bigl[\phi_2\Bigr]^n=\Bigl[h\Bigr]\Bigl[\mathcal{E}_{s12}\Bigr]\Bigl[\phi_1\Bigr]^n$$

$$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.

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.