Loss in reactivity due to Xe and Sm

During the burnup of the fuel we have production of neutron poisons, the most important are Xe and Sm. Xe-135 (Xenon) is produced both from fission and from the decay of I-135 (Iodine), Xe-135 decay \beta^- with an half life of 9 hours. Sm-149 (Samarium) is produced from the decay of Pm-149 (Promethium) and it’s stable.

Knowing the production and decay rates we can build the following lumped system:

\dot{I} = -\lambda_II + \sigma_f\Phi Y_I U^5
\dot{Xe} = \lambda_II - (\sigma_{Xe}\Phi + \lambda_{Xe})Xe + \sigma_f\Phi Y_Xe U^5
\dot{U^5} = - \sigma_a^5 \Phi U^5
\dot{Pm} = -\lambda_{Pm} Pm + \sigma_f\Phi Y_{Pm} U^5
\dot{Sm} = \lambda_{Pm}Pm - \sigma_{Sm}\Phi Sm

Where I, Xe, Pm, Sm, U are the concentration of the relative elements. Y is the Yield factor and λ the decay coefficient and finally σ is the microscopic cross-section (f = fission, a = absorption).

A way to solve this set of equations is discretize the time and perform a cycle, where the equations become:

I(t+dt) = dt[-\lambda_II(t) + \sigma_f\Phi(t) Y_I U^5(t)] + I(t)
Xe(t+dt) = dt[\lambda_II(t) - (\sigma_{Xe}\Phi(t) + \lambda_{Xe})Xe(t) + \sigma_f\Phi(t) Y_Xe U^5(t)] + Xe(t)
U^5(t+dt) = dt[-\sigma_a^5 \Phi(t) U^5(t)] +U^5(t)
Pm(t+dt) = dt[-\lambda_{Pm} Pm(t) + \sigma_f\Phi(t) Y_{Pm} U^5(t)] + Pm(t)
Sm(t+dt) = dt[\lambda_{Pm}Pm(t) - \sigma_{Sm}\Phi(t) Sm(t)] + Sm(t)

Finally we have that reactivity is defined as: |\rho_{POISON}| = f {\Sigma_{POISON} \over \Sigma_f^{FRESH FUEL}}

Show the matlab code