0
$\begingroup$

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.

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

$\endgroup$
3
  • $\begingroup$ It’s posted on the Mathematics Stack just in case they answer first ... see math.stackexchange.com/q/3258194/532586 $\endgroup$ – Solar Mike Jun 11 '19 at 4:11
  • $\begingroup$ 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. $\endgroup$ – Amit Hasan Arpon Jun 11 '19 at 4:15
  • $\begingroup$ Update, an accepted answer over on Mathematics. $\endgroup$ – Solar Mike Jun 11 '19 at 4:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.