go to main safpwr site

SAFPWR program validation.

The validation process involves the following tasks:

a) Checking that the coding correctly reflects the mathematical model.
b) Checking its model on reference problems whose exact mathematical solution is available, or against well-defined benchmarks whose reference solution has been obtained by "reference" programs.
c) As explained in the introduction, SAFPWR vocation is not to carry best estimate evaluations, but to insure that the bundle of model assumptions and correlations incorporated or introduced in the program leads to a pessimistic evaluation of safety margins.
This can only be achieved by performing sensitivity calculations allowing to ascertain that all model "parameters" are on the safe side.
Checking closure of field balance equations

As the SAFPWR model has no rescourse to any linear approximation of the original model, all the field variables balances should be closed exactly for each component of the system.
Checking global closure on a collection of components is a also a powerful debugging method because it allows also checking the connection flows between the components.

A special coding has been developed for testing global balances on the "0" sub-system including the primary system without the pressu "oxp" sub-subsystem (the Outlet + eXpansion line + Pressu vessel) which is tested separately.
Checking calculation is triggered by inserting the keyword primary_bal in the ini, step and do or redo sequences.[?] for an example.
The results are edited on output file
ss_m2_0 (1)) represents the eos mass of O calculated from global balance.
wyl are the flows injected into l; weo + wss is the net flow to oxp.
Likewise (2) is the global enthalpy balance;
q2c is the core contribution, ql is the net heat rate collected by l from pipe wall, steam-generator and pump.
(3) is the boron balance.
If the coding is correct, the values should be identical (except for rounding errors) to the summation s_m2_0, of 0 subsystem components values (4).
The same comparison is made for the system oxp.

Note that the mass accumulation in the spray line is neglected: wss represents the spray flow into the pressu as well as the spray flow taken at outlet of loop 1.
The oxp balances have been calculated for the drift pressu model only. v2xj is the eos volume of xj element and v2p the eos volume of the pressu region p.

Heat balance checking of fuel rods and pressu wall and sg are inherent to each individual models.

Test vs Nordheim asymptotic solution

Input data: [05_cntpt.dat] (download)
Point kinetics is simulated by setting
k9c=i9c=1 and zero leakage (clogbc=cloghc=0).
Core_config: 1 configuration, 1 density state.
Nordheim formula requests the delayed precursor data under Lstgg and
Λ= svog/snog= 24e-6 (typical for PWR core) under Lstsxvg.
sdog is not relevant (no diffusion in point kinetics).
sknog is arbitrary: it allows calculating f2c from the flux.
Uniform core reactivity is conveniently inserted here as negative values in the Lstcfb/rau(sec) interpolator [Reactivity;Arret;Urgence].
A 600e-5 reactivity is inserted as a ramp terminated in .1 s.
Steeper ramps are never encountered in practical problems.
oqc= 1 allows by-passing heat conduction in fuel rods (q2c=f2c).
Prompt-jump is active whenever ypjc < epsf3c.
In the present test no flux iterations are required (it9fc=1) because there is no feedback, and
dsec is imposed (it9dsc=0).
The first part of the application uses small dsec to cover the fast ramp reactivity release.
After .2 s, dsec is relaxed to pursue the case till asymptotic condition.

Chart1:2 gives the Nordheim plots omega(rho) for positive and negative reactivity.

A first test application simulates a positive reactivity injection of 600e-5 in .1 s.
The corresponding Nodheim omega value is 1.17.
Chart 3 compares results for 5 values of parameters dsec and epsf3c.

Case1: dsec=.05 and epsf3c=.1: prompt-jump approximation is enabled (when ypjc > epsf3c) for step sec=.1 only (graph 4).
Chart 3 shows that the reference omfc of 1.17 is correctly predicted after about 4 sec.
Case 4: (dsec= .05 / epsf3c=0): If prompt jump is deactivated all the time (epsf3c=0), it is necessary (chart 5) to reduce dsec down to .01 to get the same asymptotic omfc and a non-physical undershoot is observed around .1 s.
Chart 6 shows a larger scale picture.
Case 2 (chart 3) using larger dsec (.5) from .2 s and epsf3=.1 leads to a lower asymptotic omfc.

We can therefore observe that activating PJA allows getting best results with relatively large dsec.
This conclusion is valid thanks to the very short neutron life time and the results should not be sensitive to that values as long as it remains in the same order.
This is of course an instructive observation insofar as an accurate estimation of lambda is difficult.

Large Negative reactivity insertion test
05a_cntpt.dat (download) presents a point kinetic case with a step reactivity decrease of -0.16 in .1 s.
Here point kinetics is simulated by setting i9c=k9c=1 and by a linear drop of configuration ig=2 of reactivity weight -.16 in .1 s.
For such a large negative values omfc (chart 2) reaches a constant value of .0128 which is actually the decay constant amjg(1) of the longest life time decay precursor.
chart 7 shows that with dsec= .1 up to sec=.1 and dsec=10 later, omfc actually converges towards .0128.
Immediately after a step drop the decay precursors contribution to flux keep their initial value, and as the PJ approximation is valid all the time for negative reactivity insertion, it is easily shown that the initial value of the ratio f2c0/f2c should be close to (be-rc)/be =(760-600)/760 =22.05.
As, in this application, the drop lasts .1 s, we must extend the fitting of 1/f2c after the drop (.1 s) and over a time interval (.1,.2) where dsec has been kept at the low value .01 s.
chart 8 confirms that this fitting intersects well the y axix at 22.05.
This observation suggests a convenient method for measuring the drop reactivity.

Weak Negative step reactivity drop

The application is now repeated for a drop of 101 e-5 in .1 s.
dsec is kept at .01 s up to .2 s, and relaxed to 10 s thereafter.
Nordheim omfc(-101e-5)= -0.0075 (chart 2).
Here, (chart 9) omfc reaches only .0071 after 1000 s.
With dsec diminished to 1 s the asymptotic omfc reaches .0074, just below the reference. This is because the shorter half-life decay products still contribute for this minor reactivity insertion.
chart 10 shows that the linear fitting method on post-drop (.01, .02) interval leads to a value 1.13017 fairly close to the theoretical one (760+101)/760= 1.1328.

Step insertion of 600e-5 by space kinetics


(cos_600_rau.dat) (download). The same core now is represented in spatial kinetics with 14 equal meshes.
Boundary conditions are zero current (clogbc=0) at bottom and zero flux (cloghc=1e24) at top of the core, so that the upper half only is actually simulated.
For uniform critical (λ=1) core the planar reactivity r0ci, calculated by ini, core should be identical to (dΔφ)/(fφ) in order to compensate neutron axial leakage.
For an uniform core the theoretical flux distribution is
cos(πz/azc) where z is taken from mid plane and azc is the total core height. Thus,
rci= -sd(π/azc)2/sn=
1.4/.03(π/365.001)2= 345.7e-5.
A cos f2ci distribution is imposed at input.
Uniform reactivity insertion is simulated by -rau(sec) like for case above.
Chart 11 confirms that the cos flux shape 14/f2c/fci remains identical to the fundamental shape in the course of the transient, and that the calculated mesh reactivity rci calculated at ini phase is very close to the theoretical value of 345.7e-5.
The later deviations are probably due to the limited precision in entering the f2ci.
As anticipated, f2c and omfc evolution are identical as for the point kinetics representation.

.1 sec step insertion 600e-5 by group 1 withdrawal

Here (cos_600_zgri.dat) (download)the step insertion of 600e-5 is simulated by true linear withdrawal, in .1 s, of a 600e-5 reactivity configuration 1 from the critical configuration 2.
After a small flux profile deformation during the group 1 displacement (chart 11), the profiles recovers rapidly is initial fundamental shapes because the core becomes again uniform.
Likewise, omfc and f2c deviate from the uniform case predictions solely in the course of group 1 withdrawal, and omfc recovers (chart 12) its Norheim value rapidly later as the higher flux harmonics vanish.
However, as the result of the initial perturbation lasting about .2 s, f2c remains permanently about 2% higher later on.

Conclusion: the tests show that SAFPR reproduces fairly well the theoretical predictions.
Best results are always obtained by enabling prompt-jump approximation, which allow fairly good agreements, even with relatively large dsec of .05 .