Point reactor kinetics—direct and inverse

  • Difficulty: 042/100
  • Author: jeremy theler
  • Keywords: CONST, vecsum, INCLUDE, :=,

These examples deal with the point reactor kinetics equations.

direct.was

INCLUDE parameters.was    # kinetic parameters

# our phase space is flux, precursors and reactivity
PHASE_SPACE phi c rho

end_time = 60

# prescribed reactivity
FUNCTION react(t) INTERPOLATION akima DATA {
0       0
0.5     0
1       0
2       1e-5
3       3e-5
4       7e-5
5       9e-5
6       10e-5
7       10e-5
8       10e-5
10      10e-5
20      10e-5
30      3e-5
40      0
50      0
60      0
}

# initial conditions for the DAE system
rho_0 = 0
phi_0 = 1
c_0(i) = phi_0 * beta(i)/(Lambda*lambda(i))

# DAE system (reactor point kinetics)
rho .= react(t)
phi_dot .= (rho-Beta)/Lambda * phi + vecdot(lambda, c)
c_dot(i) .= beta(i)/Lambda * phi - lambda(i)*c(i)

# write the result to the standard output
PRINT t phi 1e5*rho
$ wasora direct.was > direct.dat
$ qdp direct.dat
$ 

inverse-integral.was

DEFAULT_ARGUMENT_VALUE 1 direct.dat
INCLUDE parameters.was    # kinetic parameters
FUNCTION flux(t) FILE_PATH $1

VAR t'   # dummy integration variable

# define a flux function that allos for negative times
phi(t) := if(t<flux_a, flux(flux_a), flux(t))

# compute the reactivity with the integral formula
rho(t) := { Lambda * derivative(log(phi(t')),t',t) +
  Beta * ( 1 - 1/phi(t) *
   integral(phi(t-t') * sum((lambda(i)*beta(i)/Beta)*exp(-lambda(i)*t'), i, 1, nprec), t', 0, 1e4) ) }

PRINT_FUNCTION rho MIN flux_a MAX flux_b STEP (flux_b-flux_a)/1000
$ 

inverse-dae.was

DEFAULT_ARGUMENT_VALUE 1 direct.dat
INCLUDE parameters.was    # kinetic parameters
FUNCTION flux(t) FILE_PATH $1

# our phase space is flux, precursors and reactivity
PHASE_SPACE phi c rho

end_time = 60

# initial conditions for the DAE system
rho_0 = 0
phi_0 = 1
c_0(i) = phi_0 * beta(i)/(Lambda*lambda(i))

# DAE system (reactor point kinetics with fixed flux)
phi .= flux(t)
phi_dot .= (rho-Beta)/Lambda * phi + vecdot(lambda, c)
c_dot(i) .= beta(i)/Lambda * phi - lambda(i)*c(i)

# write the result to the standard output
PRINT t rho
$ wasora inverse-dae.was > inverse-dae.dat
$ 

pukin.was

INCLUDE parameters.was

end_time = 60
dt = 0.05

# prescribed reactivity
FUNCTION react(t) INTERPOLATION akima DATA {
0       0
0.5     0
1       0
2       1e-5
3       3e-5
4       7e-5
5       9e-5
6       10e-5
7       10e-5
8       10e-5
10      10e-5
20      10e-5
30      3e-5
40      0
50      0
60      0
}


VECTOR cbez SIZE nprec

# condiciones iniciales
pkbez_0 = 1
cbez_0(i) = pkbez_0

cbez(i) = lag(pkbez, 1/lambda(i))
pkbez = lag(vecdot(beta, cbez)/(Beta - react(t)), Lambda / (Beta - react(t)))

PRINT t pkbez 1e5*react(t)
$