I am currently trying to understand the TUBRNP model which is used to calculate the isotope compositions and the radial power profile evolution with the burnup. I am talking about the basic model (the one from 1994 with only 6 isotopes taken into account) in the case of a LWR with UO2.

In order to understand it, I try to create a simple algorithm doing the same thing. i will be the burnup step and n the inner iteration index for convergence.

1) From an initial composition vector N_0 (determined with the enrichment), calculate the parameter A.

2) Determine the matrix M so that I can calculate the next composition N(n)_i+1 with N(n)_i+1*M=b with b depending on the old concentration N_old (see below the code for a given position):

Code:

```
function [M,b]=create_matrices(delta_BU,A,sigma_abs,sigma_capt,radial_pos,shape_function,N_average,N_old)
f_r=shape_function(radial_pos);
M=zeros(6);
%U235
M(1,1)=1+sigma_abs(1)*A*delta_BU;
%U238
M(2,2)=1;
%Pu239
M(3,3)=1+sigma_abs(3)*A*delta_BU;
%Pu240
M(4,3)=-sigma_capt(3)*A*delta_BU;
M(4,4)=1+sigma_abs(4)*A*delta_BU;
%Pu241
M(5,4)=-sigma_capt(4)*A*delta_BU;
M(5,5)=1+sigma_abs(5)*A*delta_BU;
%Pu242
M(6,5)=-sigma_capt(5)*A*delta_BU;
M(6,6)=1+sigma_abs(6)*A*delta_BU;
b=N_old;
b(2)=b(2)-N_average(2)*sigma_abs(2)*f_r*A*delta_BU;
b(3)=b(3)+N_average(2)*sigma_abs(2)*f_r*A*delta_BU;
```

4) Calculate the power profile from the average power that I want and the flux.

5) Determine the burnup increment from the power, the initial uranium mass and the simulated time increment

6) Compare the compositions N(n)_i+1 and N(n-1)_i+1. If the difference is small enough, go to next burnup step, otherwise take N(n-1)_i+1 to calculate a new A and a new matrix M, ...

First, does this reproduce properly TUBRNP? I took the cross sections and the p_i parameters for the shape function from the FRAPCON user manual.

Then I have several issues with the shape function f(r). As there is a volumetric normalization step, the outer radius value matters. Thus, the unit of this parameter influences the results and unfortunately in the publication, no unit is given. If I try to reproduce the Plutonium radial profiles given in this paper with cm as units, the shape is not peak enough, but if I try with mm, it is not working either.

Finally, I face some problems calculating the conversion factor alpha. From what I understand, it is a factor created to convert time to burnup units. As I am working with MWd/kg_HM as burnup unit and cm as spatial unit, we have : [BU]=[alpha]*(fission/cm^3)/(kg/cm^3) =MWd/kg_HM

So for me : [alpha]=MWd/fission

alpha=200MeV/reac*1.602e-19/(24*3600)=3.7088e-22

If my calculation is correct, then why in the Lassman's paper the value of alpha is 3.35e-16 in the case of MWd/t_HM ? This 1e6 factor does not make sense to me.

I would really appreciate your help on this ! And if you have questions or if I have not been clear enough on a point, please ask me.