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$ Jun 11 '19 at 4:15
  • $\begingroup$ Update, an accepted answer over on Mathematics. $\endgroup$
    – Solar Mike
    Jun 11 '19 at 4:35
1
$\begingroup$

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})}$$

$\endgroup$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.