PWR (quasi) non-linear dynamics solver

Previously i discussed here how to solve linearized lumped equation for PWR core dynamics. However linearized system have some limitation, so here we are, in this work i mantain the non linearity of power equation but reactivity feedback remain linear. For this reason i used the term quasi non-linear (almost non-linear)

Index:

• DAEs system of equations
• Initial conditions
• Matlab solver

DAEs system of equations

Equation are the same as before, 1 Group of precursors and incompressible thermalhydraulics:

Power and Precursors:

\frac{dP}{dt} = \frac{\rho - \beta}{\Lambda}P +\lambda C
\frac{dC}{dt} = \frac{\beta}{\Lambda}P - \lambda C

Fuel and coolant energy balance:

c_fM_f\frac{dT_f}{dt} = P - K(T_f - T_c)
c_cM_c\frac{dT_c}{dt} = K(T_f - T_c) -\Gamma c_c(T_{out} - T_{in})

Reactivity:

Coolant equation can be modified substituting outlet temperature:

c_cM_c\frac{dT_c}{dt} = K(T_f - T_c) - 2\Gamma c_c(T_{c} - T_{in})

Due to the reactivity this is a DAE system (differential-algebraic system of equations)

Initial conditions

To solve this system with matlab it’s important to define initial condition for our variables and also for the derivatives.

For power we exploit the continuity:

To obtain the precursos concentation we exploit precursors equation at steady state:

Initial fuel temperature and mass flow rate are user defined.

Defining initial inlet coolant temperature is possible to obtain the average coolant T:

From energy balance: P(t=0) = \Gamma(t=0^-) C_c (T_{out}(t=0^-) - T{in}(t=0^-))

Then initial coolant temperature is: T_c(t=0^+) = \frac{T_{out}(t=0^-) + T_{in}(t=0^+)}{2}

Finally reactivity initial condition is:

Matlab solver i implemented use ode15i which require both initial condition for state variables and for their derivative, so:

\frac{dP}{dt}(t=0^+) =\frac{\rho (t=0^+)}{\Lambda} P(t=0)

\frac{dC}{dt}(t=0^+) = 0

\frac{dT_f}{dt}(t=0^+) = 0

\frac{d\rho}{dt}(t=0^+) = \alpha_c \frac{dT_c}{dt}(t=0^+) + \alpha_h \frac{dh}{dt}(t=0^+)

Matlab Solver

Definition of symbolic variables and equations:

syms P(t) C(t) T_f(t) T_c(t) rho(t) h(t) Tin(t) G(t)


%% set of equations
eq1 = diff(P) == (rho - b)/L*P + l*C;
eq2 = diff(C) == b/L*P - l*C;
eq3 = M_f*C_f*diff(T_f) == P - k*(T_f - T_c);
eq4 = M_c*C_c*diff(T_c) == k*(T_f - T_c) -2*G*C_c*(T_c - Tin);
eq5 = rho == alpha_f*(T_f - T_f0) + alpha_c*(T_c - T_c0) + alpha_h*(h-h_ref);

%INPUT EQUATIONS _ MODIFY ALSO INITIAL CONDITIONS!
eq6 = Tin == Tin_0 + sin(2*pi*t/T_FIN*30);
initial_Tin = Tin_0;
initial_der_Tin = 2*pi/T_FIN*50;
eq7 = G == G0 ;
initial_G = G0;
initial_der_G = 0;
eq8 = h == h_ref;
initial_h = h_ref;
initial_der_h =0;

As you see in last lines it’s possible to define arbitrary function for inlet temperature, mass flow rate and rod insertion (h) however initial conditions need to be coherent!

Build variable vector and store the dimension:

%% variable definitions
vars = [P(t);C(t);T_f(t);T_c(t);rho(t);Tin(t);G(t);h(t)]; %column vect var
origVars = length(vars);

Reduce order (may be not neccesary)

%% reduce order
%in case order is hugher than 1 must be reduces
[eqns,vars] = reduceDifferentialOrder(eqns,vars);
if isLowIndexDAE(eqns,vars) == 0 %then index must be reduced
    [DAEs,DAEvars] = reduceDAEIndex(eqns,vars);
else
    DAEs =eqns;
    DAEvars = vars;
end

Conversion from symbolic to numerical 

%% conversion to matlab function handle

f = daeFunction(DAEs,DAEvars); %function handle
F = @(t,Y,YP) f(t,Y,YP);
%YP is first derivative of variables Y

Impose initial condition and run decic (which corrects possibles errors)

%% impose initial conditions

%conditions for true variables

y0_guess_true = [P0;C0;T_f0;((Tin_0 + DT_io) + initial_Tin)/2;alpha_h*(initial_h-h_ref);initial_Tin;initial_G;initial_h];

%define all guess initial conditions (also for non physical variables)
y0_guess = zeros(size(DAEvars));
y0_guess(1:length(y0_guess_true)) = y0_guess_true;

%define guess initial condition for derivative
yp0_guess = zeros(size(DAEvars));
yp0_guess(1) = (y0_guess(5) - b)/L*P0;
yp0_guess(4) = (k*(y0_guess(3) - y0_guess(4)) -2*y0_guess(7)*C_c*(y0_guess(4) - y0_guess(6)))/M_c/C_c;
%yp0_guess(4) = 2*G0*(initial_Tin - Tin_0)/M_c;
yp0_guess(5) = alpha_c*yp0_guess(4) + alpha_h*initial_der_h;
yp0_guess(6) = initial_der_Tin;
ypo_guess(7) = initial_der_G;
yp0_guess(8) = initial_der_h;

opt = odeset('RelTol', 10.0^(-7),'AbsTol',10.0^(-7));
%find consistent initial value
[y0,yp0] = decic(F,0,y0_guess,[],yp0_guess,[],opt);
% in some cases if the guessed initial condition is wrong then decic give
% error

Almost done, let’s solve the system with ode15i:

%% solve system

%solve system with ode15i(system, time span, initial variable, initial derivative, options)
[tSol,ySol] = ode15i(F,[T_0 T_FIN],y0,yp0,opt);

And now let’s plot the results:

for k = 1:origVars-3
  
  subplot(origVars-3,1,k)
  plot(tSol,ySol(:,k),'LineWidth',2)
  title(title_vector(k));
  xlabel('time [s]')
  grid on
  
end

figure

for k = 6:origVars
  
  subplot(origVars-5,1,k-5)
  plot(tSol,ySol(:,k),'LineWidth',2)
  title(title_vector(k));
  xlabel('time [s]')
  grid on
    
end

With the selected input (oscillatory inlet temperature) the solution is the one in the picture below:

Note: Precursors unit of measure should be Watt and not Ws as in figures

Another example is the response to a step rod extraction of 10 cm:

Finally it is possible to study the superposition of different effects, for example a linear increase of the flow rate plus an oscillatory behavior of the inlet temperature: