# How to code this iterative scheme for solving 2-group diffusion equation?

I am trying to write a python code to solve the neutron diffusion equation to model neutron flux distribution in a one-dimensional two-group setting. 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.

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[\nu_1\mathcal{E}_{f_1}\Bigr]\Bigl[\phi_1\Bigr]^{n-1}+\frac{1}{k}\Bigl[\nu_2\mathcal{E}_{f_2}\Bigr]\Bigl[\phi_2\Bigr]^{n-1}+\Bigl[\mathcal{E}_{s21}\Bigr]\Bigl[\phi_2\Bigr]^{n-1}$$

$$\Bigl[C_2\Bigr]\Bigl[\phi_2\Bigr]^n=\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 code 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.

• It’s posted on the Mathematics Stack just in case they answer first ... see math.stackexchange.com/q/3258194/532586 Jun 11 '19 at 4:11
• Yes, I posted there seeing as the problem I want help with is related more to mathematics than engineering. But I wanted to post here too because, in the end, it is a part of nuclear engineering first and foremost. So there might be people here who have previous experience of solving this particular problem. Jun 11 '19 at 4:15
• Update, an accepted answer over on Mathematics. Jun 11 '19 at 4:35

You haven't given enough information to answer this exactly, but I am assuming from your equations that you are solving mesh-centered equations and the assumption is that the flux is linear in each cell.

In that case, the average value of the flux in each cell is taken as the midpoint, and your integration just becomes a sum $$k^n=k^{n-1}\frac{\sum h(\nu_1\mathcal{E}_{f_1}\phi_1^n+\nu_2\mathcal{E}_{f_2}\phi_2^n)}{\sum h(\nu_1\mathcal{E}_{f_1}\phi_1^{n-1}+\nu_2\mathcal{E}_{f_2}\phi_2^{n-1})}$$