SAFPWR: pwr safety and transients computer program general description.
Aims: analyse core safety and protection of Westinghouse-type nuclear Pressurized Water Reactors
SAFPWR program modeling assumptions
Global pressure approximation: amounts to evaluating local water massic volume from its local enthalpy, but from a single uniform global pressure. This approximation has the effect of "erasing" pressure waves, of relieving numerical stiffness and allowing larger time-steps.
In the fluid nodal energy balance equations, by dropping the kinetic energy and gravity terms.
In the course of the accidents and transients, the axial shape of the core power distribution may deform rapidly, which has a major impact on the core safety margins.
This allows conceiving the spatial kinetics representation of the core as resulting from a formal axial collapsing of the general 3D-kinetics model.
As the radial distribution is not explicitly calculated, a "hot channel" must be attached to the core in order to evaluate conservative values of peak power density, fuel temperature, departure of nucleate boiling.
At the occasion of 3D-kinetic benchmarks it will indeed be shown that our axial collapsing model reproduces surprisingly well the results obtained by reference 3D_ kinetics codes.
Let us again stress that the goal of the program is not to reach "best estimate", but rather "conservative" values for the safety-related variables, in conformity with the "licensing" approach.
Except for the above simplifications, the SAFPWR mathematical model is basically the same as in the conventional system programs: conservation equations of mass, energy and kinetic momentum on homogenous nodes
For the sake of maximum portability, SAFPWR was totally programmed in modern FORTRAN 95 language.
The field equations integration method (for water mass, enthalpy, boron, kinetic momentum, for internal energy of metallic parts and for neutron density, delayed neutron precursor density, residual heat fission products density...) aims at solving exactly the original set of those highly non linear constitutive equations.
Numerical stability is achieved by resorting, whenever feasible, to totally implicit (hence stable) numerical schemes, where feedback effects are applied at end-of-step.
The solution scheme also guarantees strict conservation of all extensive (mass, energy, volume, neutron density,...) original field values, even for large time-steps.
This is of course a necessary, but not sufficient condition: it must in addition be demonstrated that the time integration of the field equation is correct.
This can be verified by refining the time stepping and by comparing to theoretical and benchmark solutions.
In addition the strongly implicit scheme alleviates the "Courant stability limitation" which, in some system codes, necessitates reducing node volumes when reducing time-step duration.
New thermal-hydraulic developments.
A non-equilibrium, Lagrangian, pressurizer model has been developed for making it possible simulation of pressurizer and expansion line voiding.
A rigorous "Lagrangian" solution scheme for water transport and heat exchange in the primary loops and down-comer was developed with the purpose of quantifying the effect of "forward numerical diffusion " encountered in the classical eulerian scheme, on the response of over-power and over-temperature protections systems installed on the primary branches of steam-generator.
Although the primary goal the program of educational nature, it should also be able, if fed by proper system and neutronic data derived from static design, to be used as well for practical safety evaluations.
Prerequisites for using SAFPWR computer program
A basic knowledge of PWR system and core kinetics is required for a knowledgeable use of the program. .
First SAFPWR accident application exemple
The first accident application exemple analyses an elementary water coolant transport problem over a linear string of nodes.
The totally implicit eulerian scheme solution depicts the typical unrealistic forward diffusion effect characteristic of the eulerian integration scheme.
Is numerical diffusion effect observed for Euler solution a critical safety issue?
The consequences of numerical diffusion deformation on the response of the temperature sensors in the core "over-temperature" and "over-power" protection systems will be analysed in depth in the chapter on protections.
Solving balance equations
As explained earlier, thanks to the implicit solution scheme, the severe restrictions imposed by the Courant numerical stability criterion are avoided, which makes it feasible to accept large time-steps and volumes without jeopardizing numerical stability and also allows evaluating solution accuracy by simply decreasing node volumes.
The steam-generator deserves a reasonably realistic modeling in view of the important potential reactivity release, caused by a cooling accident, at end of cycle, when the temperature coefficient is large.
It features modelling water natural recirculation, level regulation and steam-water separation in the dome. This should provide a fair prediction of mass and enthalpy inventories and heat-exchange through the steam-generator tubes.
Safety analysis on simple PWR complete system
SAFPWR computer program description
Properties of water and steam.
Properties of sub-cooled water
The coefficients are obtained by least-square fitting on VDI tables over a range of volumic mass, enthalpy and pressure which should bracket the range of values encountered during the analyzed transients.
The function temperature versus enthalpy and pressure is represented by the Kaisermann correlation , which is quite accurate.
Properties of over-heated steam
Nodal mass, enthalpy and boron balances
For Eulerian balances, the control volume is the node geometrical volume .
For the purpose of setting up the balance relations, the node can be viewed as a perfect mixing cell whose content remains homogeneous.
Heat conduction model
Eulerian biphasic drift model
This model accounts for the relative Celerity between the vapor bubbles and the surrounding water
The installed drift model is intended to answer the questions:
has drift a sensible effect on the recirculation flow, on the mass and energy inventory in the steam-generator core, on the secondary heat transfer, on the time-history of an accident such as the Main Steam Line Break?
In other words, is it a safety issue?
Lagrangian homogeneous model
As explained in the introduction, a special lagrangian balances model had to been developed in view of overcoming the "forward diffusion" effect associated with the Euler balances scheme.
Biphasic drift lagrangian balances
In pressurizer system, for non equilibrium condition, the eulerian drift model developed for the steam-generator core does not easily apply because the flow distribution may rapidly change in value and even direction in the course of a time-step.
Furthermore, the vapor region must be explicitly represented, as well as the bubble rise from liquid region.
Therefore, a lagrangian drift model is more appropriate.
Kinetic momentum balances
It is a single loop case with thermal power generated by a heat-generator
The initial flow is imposed as volumetric flow.
Pump manometric head is cut to zero from 20 s to simulate flow slowing down to natural circulation.
Steam generator is simulated through an uniform secondary temperature field
After flow slowdown, natural circulation will be activated by generating a power of 2*10e6 from 400 s.
Fuel rod thermal model
Thermal conduction model in pipes
The model is basically similar as that of the fuel rod, without heat source in conductor
Reactor kinetics model of SAFPWR computer program
Multi-group 3-D kinetics
We start with the multi-group kinetic equations in diffusion approximation, which are written as (1) in order to single-out the neutron axial leakage.
models: neutron density, scalar flux, fission of fissile isotope, delayed neutron emitters, prompt neutrons life-time
Prompt neutrons are emitted into a "prompt" velocity distribution spectrum
The delayed neutrons are emitted into another spectrum softer than prompt neutron emission spectrum.
In summary, neutron density balance:
neutron density time derivaive = prompt source - removal - axial leakage + source by delayed emitters decay;
The first term: production rate of emitters
Second term: removal rate by radio-active decay
Setting to 0 the time derivatives leads to the classical form of balance, where the core eigen-value lambda needs to be inserted to insure core criticality and defines, for fissile isotope i, the stationary, "composite", average emission spectrum for prompt and delayed neutrons. This spectrum is part of the fuel assembly neutronic library.
2-group axial reduction
The neutron importance of family j is defined .
In analogy with the definition of static core reactivity, it is possible to define a planar dynamic reactivity
where it is reminded that the removal term includes the radial leakage but excludes axial leakage.
The obvious interest of adjoint flux weighing is to select, as reduction invariant, the most significant parameter of core physical state, namely the static reactivity in each of its planes.
In terms of static core reactivity we come back to the same Nordheim relation as in point kinetics, because static flux verifies eigen value equation.
Numerical Integration of reactor kinetic equations
The delayed emitters balance equation can be integrated analytically.
The "Prompt-jump approximation" can be accepted whenever the contribution of the term can be neglected.
Solving the flux equation
The method for obtaining a feed back converged solution for the flux equation involves 4 successive calculation stages which must be repeated until convergence.
A control of acceptability of the prompt-jump approximation.
This approximation is generally always valid, except, in fast reactivity transients, during the very short periods when the core approaches or crosses prompt criticality.
Whenever acceptable the Prompt Jump approximation is very effective in so far as it allows, without impairing the solution quality, to tolerate time steps dsec very large compared to prompt neutron life-time. The described algorithm accommodates for marked deformation of the axial profile in the course of the time-step.
It provides an accurate, converged solution at eos, even if the water reactivity temperature coefficient is positive.
Tabulation of neutronic properties
The neutronic properties (reactivity, macroscopic cross-sections, delayed neutron beta's) must be assigned to each element of the partition
It aims at generating a complete initial steady-state condition consistent with a minimal set of design or operational data, namely:
By back-solving the conservation equations, the program then computes successively:
For the flows in loops and in reactor vessel, a multiplying factor correcting the friction pressurizer loss coefficients in such a way that, with the input manometric head, we reproduce the input volumetric flow.
The initial mass, enthalpy and boron of each node are in turn calculated.
For the pressurizer, the power of electric heaters necessary to close the pressurizer enthalpy balance, assuming that core outlet receives saturated water from pressurizer expansion line.
In the core, the neutronic tables correlate, for each core configuration, the neutronic properties with the physical effects actually modelled in SAFPWR (coolant density, boron, effective fuel temperature ).
However, the slow "historical" effects, such as fuel depletion, Xenon and Samarium build-up, which are supposed invariant in the course of the transient, do not need to be explicitly represented:
Assume that, at initial steady state, is available a reference solution accounting for all the effects not explicitly modelled by the program.
This reference solution must satisfy the stationary form of neutronic balance, where lambda is an eigen-value accounting for the fact that the core model selected for generating the neutronic tables is not necessarily critical.
safpwr has now the capability to solve the initial condition eigen value problem
Conceiving the power decay model as a strict analogy to the delayed neutron model allows taking profit of the existing coding. However, as power production is no longer synchronous whith neutron flux, modeling the interaction between neutronic and power effects (for example, the power feedback to reactivity) is more complicated
As explained in the introduction, a hot channel model has been attached to the core average channel, in order to assess commonly used safety quantities such as maximum fuel temperature, fuel heat content, maximum clad temperature, dnbr (departure from nucleate boiling ratios,...), used in safety evaluations.
Departure from Nucleate Boiling Ratio
It is the heat flux at onset of boiling crisis from which clad temperature suddenly increases by degradation of film heat-transfer.
This value, calculated by empirical correlations, depends on the local thermal properties and their upstream distribution.
The pressurizer system includes the core outlet, the expansion line, the pressurizer vessel itself and the Spray system.
General description of multi-regions pressurizer model
For severe transient like the Main Steam Line Break (MSLB), the strong water out surge from the pressurizer system caused by the rapid cooling of primary water may result in complete collapsing of the liquid region in the pressurizer and expansion line.
The water in dome may even enter boiling because its temperature remains warmer and the dome, hence, behaves then as a second pressurizer!
For non-isolated MSLB the core bottom may also experiment moderate bulk boiling.
In order to tackle those extreme situations, a multi-nodal, multi-elements, non-equilibrium, Lagrangian model had to be developed
In the course of SLB transients, the possible return to power causes a rapid water in surge into the emptied pressurizer.
The pressurizer behaves then more as a steam gas pressurizer and this will result in a rapid pressurizerre surge, which is beneficial for opposing to boiling crisis, but at the same time, reduces the boron flow injected by Safety Injection water.
As spray is not active (pressure still too low) at time of water surge, the main mechanism for checking the pressure surge is vapor condensation to wall.
But this can only take place in so far as the wall edge temperature remains below the steam temperature and the condensation latent heat can be absorbed by the wall.
Safety Injection System
Reactor protection by rod dropping
Reactor protection trip (dropping all the rods) is triggered whenever some core protection variable values, exceed their trip threshold level
An example of the trip logics is illustrated on fig for a rod withdrawal accident.
Actually f2c is the nuclear power as monitored by neutron flux detectors. This signal gives only a relative estimation of nuclear power and needs therefore to be rescaled to a "reference" thermal power calculated by core heat balance.
"dphdt" trip protects the core in case of an abnormal flux time derivative caused, for example by an accidental, undetected, drop of one or several rods. This trip protects the core against the effects of an undetected rod drop.
Core thermal protection
The high flux trip is fast, but not accurate, because the measured flux provides only a relative estimation of power, which must be recalibrated on the thermal power taken from the system heat balance.
It protects the fuel against excessive temperature in case of fast reactivity transients such as rod ejection, but does not allow assessing the coolant heat removal capability.
Designing core thermal limits starts (core thermal design) with representing , of core inlet temperature versus core thermal power, the limit lines corresponding to each phenomena which could endanger safe core operation.
A first limit would cover core outlet boiling, a second, overheating of fuel at hot spot due to overpower, a third, the loss of coolant heat removal capability by boiling crisis (loss of Departure from Nucleate Boiling), ....
These core protection "physical limits" depend also on additional parameters such as pressurizer, core flow, core enthalpy form factors, core axial power profile, .....
Consequently, the only procedure for assessing the protections is by submitting the core to set of transients and accidetns sensitive to the dynamic effects and evaluating the protection margins directly in terms of the real safety limits such as max fuel temperature at hot spot, DNBR margins,....
Note that SAFPWR allows nevertheless editing the margins against physical limits which would be calculated on a thermal core model (multi channel,....) more realistic than in our model.
However core thermal design is usually done in steady state condition, with enthalpy distributions which may markedly differ from the transient ones.
In order to investigate this effect, the program may also calculates, on demand, the stationary condition of the core, for some "state points" (thermal power, inlet flow,...) of transient. The trick is to insert pseudo cases at chosen steps.
The drift process facilitates bubbles extraction but the resulting decrease of vapor fraction reduces natural recirculation, so that it is not clear to ascertain whether the correct representation of the drift process is critical.
The steam-generator dome region is modelled as an "ideal" representation including an ideal, zero-volume, separator returning saturated water to the annular region releasing saturated vapor to the dome.
processing starts with calculating the total flow through the pressure relieve valves, by means of the flow Moody critical discharge correlation, using the discharge coefficients.
wzpdn applies for the pressure relief devices, and wypdn simulates the discharge though the break in case of MLSB accident.
This representation avoids modeling the vapor collector and turbine feed valve systems. It amounts to neglecting the effect of vapor collector on pressure variation at turbine entry.
The feed water flow is imposed as input if the steam-generator is free (unregulated) or adjusted by the steam-generator water level regulator
instead of modeling a real regulator, we have rescourse to an ideal regulator based on the same Newton method as used for heat-exchange stabilization.
In spite of its relative simplicity, our "ideal" steam-generator model is believed to be adequate for assessing the effects of steam generator secondary side behavior on core transients.
It is however valid as long as the steam-generator operation remains "normal" ie when the core cn contains saturated vapor and natural recirculation is active.
In extreme cases (like MSLB) some steam-generator may become thermally isolated from the primary, and recirculation will cease because secondary water becomes warmer than that of the primary side.
In such a case, we have rescourse to an elementary "hot water bottle" representation of the whole steam-generator as a single homogeneous volume ".
In addition, the friction coefficients input data along the recirculation loop must be adjusted to reproduce the target recirculation data with the water level in an at its target value.
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.
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.
Test versus Nordheim asymptotic solution
Conclusions: the tests show that SAFPWR 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 time-step of .05 .
ANL one-dimensional kinetik benchmarks
In conclusion the agreement with the anl1 benchmark is quite good, which brings, particularly, also an validation of the data reduction procedure to one-group approximation.
This test illustrates the good stability and accuracy of the program, even for large steps.
The ability of correctly falling back on the good trend, even after missing the top peaks, is doubtless thanks to the second degree scheme used for integrating the neutron kinetics and to the analytical scheme for integrating the delayed precursors equations.
NEA 3-D Kinetics Group Withdrawal benchmark
The reactivity ramp before power peak is about 17e-5/s. For such a slow rate, a first "Doppler" feedback-controlled power peak is observed at about 2 s after core reactivity has reached (at 70 s) its peak, just below beta= 760e-5, but at a power level still under the indicated power trip level of .35 nominal power.
Thus, the power keeps increasing, but at a slower rate due to the Doppler feedback, and reaches eventually the trip level of .35 at 81.3 sec.
The rods start dropping dsecau=.6 s later, which causes immediate reduction of all safety related items.
Occurrence of this singular double peak pattern was not anticipated, but revealed itself as a severe test for the benchmarking exercise.
Direct comparison with the published benchmark results indicates that SAFPWR results are well within the dispersion range of the best benchmark contributions, except however for a tendency to overestimation, mainly for hot channel heating items, which was expected.
SAFPWR application exemples
Cause of accident: enforced core coolant flow reduction resulting from an accidental drop of electric supply grid frequency.
The interest of calculating dnbr on stationary core is to make it possible to apply licensed thermal design methodology, SAFPWR being used to evaluate the min dnbr time and check that the initial distribution of heat flux and vapor quality still apply.
The "dphdt" protection system aims at detecting the accidental fall of one or several rods.
"dphdt" protection system acts as a backup for the other drop detection systems by end-of-stroke detection and by monitoring the radial unbalance of the response 4 fission chambers, appearing in cases of centrally unsymmetrical position of the falling rods.
Testing ODT protection on a power ramp
The "Group Withdrawal accident at Power" is the typical accident protected by this system.
Control Rods Group Withdrawal Accident at Power
Static verification of dnbr at mindnbr state point
It may be advisable to confirm SAFPWR transient dnbr predictions by means of a licensed thermal hydraulic design program like COBRA, equipped with the fuel-dependent dnb correlations and allowing for a better inter-channel mixing grids representation.
Main Steam Line Break accident
The main steam line break (MSLB) accident is initiated by a full size break of one of the main steam line at outlet of one steam generator.
The accident is supposed to occur at end of cycle with the core at hot zero power condition.
All the control rods are then fully inserted, except for the most reactive one, which is supposed to remain stuck at its position at nominal core operating condition.
Core nuclear design must insure that, at accident time, the core be subcritical by at least 2%
At end of cycle, the boron concentration is close to zero, so that the moderator density reactivity coefficient is strongly positive, so that the rapid fall of core moderator temperature caused by the liquid discharge in the faulted steam-generator risks to use up the subcriticalty margin and bring the core again at power.
The core safety margins must be verified in that adverse situation
The accident is particularly of concern for several reason:
The main process for checking the accident is the counter-reactivity developed by the intrinsic core reactivity feedback's as core returns to power.
Boron injected by the safety injection system triggered at the accident onset appears to comes too late for checking the reactivity release.
Steam generators abnormal working conditions: in the faulted steam-generator, modeling the heat exchange, as the tubes start to be uncovered by saturated water, is difficult because the weak heat capacity of dry steam causes rapid steam temperature increase and the uncoupled heat exchange model fails.
Modeling steam discharge through the break is not obvious.
Normally, the unaffected steam-generator will be isolated at the accident beginning by automatic closure of the fast isolating valves, and these steam-generator ceases to behave as a natural recirculation steam-generator because the secondary water becomes warmer than the primary water.
Modeling the large temperature and flow unbalance effects through the vessel downcomer, bottom, core entrance and exit is difficult.
In the core itself, the neutronic-hydraulic interactions will cause rapid spatial power radial and axial redistributions which challenge the core calculation model.
In case primary water shrinks too rapidly before being reheated by the possible core return to power, the pressurizer becomes voided and boiling can occur in the primary system.
Boiling will generally start in the vessel dome, which tends to retain its initial temperature.
Boiling in the primary system mitigates the pressure decrease but challenges the calculation if core is still at power.
When pressurizer refills after complete voiding of the faulted steam-generator, its behavior departs totally from normal operation.
In order to tackle the problem, we will, as usual, try uncoupling the effects in order to identify the critical ones and adopt a bounding analysis approach when modeling difficulties are encountered.
We will start analyzing the accident for the highly improbable case that the isolating valves of the intact steam-generator's fail to operate, so that all steam-generator loose their water through the break.
In that case we will, however, not retain the, equally highly improbable, stuck rod assumption.
Because of the paramount importance of spatial redistribution effects in the core, it is not relevant, to start analyzing the transient with a point kinetics representation of the core.
Furthermore, in the situations where boiling appears in the core, the program execution may crash because of the strong pressure-reactivity coupling oscillations. Therefore, all the applications need to be carried out with a fine mesh representation of the core.
This will allow gaining a first feeling of the accident trend and on the importance of parameters such as core initial power, effective fraction of delayed neutrons, mox core effects, neutronic reactivity feedback's, steam discharge rate, steam-generator feed water flow, steam-generator and pressurizer initial water inventory,...
Furthermore this will provide envelope case for the normally protected MSLB.
Next, the classical unsymmetrical, single faulted steam-generator, accident event will be investigated with a 2 loop representation: loop1 feeding the faulted steam-generator and loop2 duplicated, unsymmetrical core, with, effects of isolating valves closure time, discharge rate, flow mixing parameters, ... can thus be examined.
Because of the simplification of the various input data, it is obvious that the findings gained from our first preliminary investigations must be confirmed whith more realistic and more complete data sets.
For some cases at least, it will be observed that, at most dangerous conditions, the core and system are quasi stationary, which makes it feasible checking the state point by means of core and steam-generator stationary calculations with normal 3-d core and steam-generator design tools.