Simulator Development
Simulator Development
Simulator Development
12
Workshop material
© IAEA, 2001
Printed by the IAEA in Austria
April 2001
FOREWORD
The use of particular designations of countries or territories does not imply any
judgement by the publisher, the IAEA, as to the legal status of such countries or
territories, of their authorities and institutions or of the delimitation of their
boundaries.
The mention of names of specific companies or products (whether or not
indicated as registered) does not imply any intention to infringe proprietary
rights, nor should it be construed as an endorsement or recommendation on the
part of the IAEA.
CONTENTS
2.1 Introduction................................................................................................................... 23
2.2 Basic Equations ............................................................................................................ 24
2.3 Mass and Component Balance Equations..................................................................... 24
2.4 Rate of Change of Pressure in a Control Volume......................................................... 26
2.5 Momentum Equation for Links Joining Control Volumes ........................................... 27
2.6 Linearized Momentum Equation .................................................................................. 27
2.7 Pressure Drop Calculation ............................................................................................ 28
2.8 Pump Flow Equation .................................................................................................... 28
2.9 Flow Equation in Link with Pump and Valve .............................................................. 30
2.10 Pump Cavitation Factor .............................................................................................. 30
2.11 Flow Equation to Match Specific Pump Curve and Incompressible Flow................. 30
2.12 Energy Balance Equations .......................................................................................... 30
2.13 Rate of Change of Fluid Temperature in Control Volume ......................................... 31
2.14 External Heat Transfer Rate in Control Volume ........................................................ 32
2.15 Equations of State ....................................................................................................... 32
2.16 Two Phase Systems (Simple Equilibrium Model)...................................................... 32
3.1 Introduction................................................................................................................... 35
3.2 Explicit Euler Integration ............................................................................................. 35
3.3 Implicit Euler Integration.............................................................................................. 38
3.4 System of Equations ..................................................................................................... 40
3.5 Other Important Model Considerations ........................................................................ 42
3.5.1 Range of Validity ............................................................................................... 42
3.5.2 Reliability........................................................................................................... 43
3.5.3 Real Time Simulation ........................................................................................ 44
4.0 CassimÔ Based Simulation Methodology and Simple Modeling Exercises .................... 46
In 1968, General Electric equipped its training center with the first full-
scope training simulator. That was followed by all the major manufacturers
of light water reactors in the U.S., and subsequently, by several utilities
and training institutes, both in the U.S. and overseas. Nowadays, power
plant simulators have become necessities among nuclear utilities
worldwide and have become much more than a training tool. In addition to
imparting to the operator an extensive “hands-on” experience in startups,
power maneuvering and shutdown operations, they give them the
opportunity to face more emergency and abnormal situations. They also
constitute a valuable tool during licensing examination and in the
development of operating procedures. As well, they have been used in
plant design validation as engineering simulators.
Over the span of the last 30 years, the simulator technologies have
changed significantly due to the rapid advances in computer hardware and
software. 15 years ago, computers that were used to simulate a complete
nuclear power plant down to the last electrical fuse occupied a full room
and cost millions of dollars. Nowadays, the same full scope simulation can
be accommodated by a high-end multi-Pentium chips PC sitting on your
desktop, and it only costs few thousands of dollars.
1
cover various operating conditions. It is important that the model will
be generic and flexible enough to facilitate easy customization of the
component for specific application. For example, to change a valve
model from a first order linear valve to a quick-opening non-linear
valve only requires changing certain parameters in the model. Some
vendors called these basic units “generic algorithms”; some call them
“standard model handlers”; others call them “model objects” (a more
modern terminology). In essence, they all amount to the same thing —
that is, they are all compiled into object files and stored in a library that
is readily called upon by the main simulation program.
5. The fifth step is testing to ensure that the subsystem model satisfies the
modeling requirements specification. Such requirements will include:
range of operating conditions; model fidelity ( % deviation from
operating data); model performance on selected malfunctions etc.
2
6. As the several interrelated subsystem models are completed, they will
be integrated to form a large size system model, followed by testing
and model tuning etc. leading to the final plant model.
The process for building generic algorithm library and the large scale
model is illustrated in Figure 1.
Figure 1.
3
1.2 Integral Components of a Modular Modeling System
1.3 Simulation Development System for Control Room Replica Full Scope
Simulators
(a) ROSE (Real time Object oriented Simulation Environment ) from
Canadian Aviation Electronics (CAE), Montreal, Canada.
4
schematics and generate simulation model codes in C and FORTRAN.
Specialized Code Generators include: (a) Hydraulic and Electrical Network
Solvers (2) Relay Logic Model Generator (3) Sequential Logic Generator
ABB’s CETRAN has been used for real time full scope simulators
application for nuclear and fossil power plants. CETRAN software is block
structured. In block structured approach, the model development engineer
creates blocks from an algorithm library, enters physical and
thermodynamic coefficients to make the model match the physical plant,
and then links the blocks together interactively and on-line. A CETRAN
Model is created by interconnecting pre-programmed blocks (e.g. heat
exchanger, time delay etc.). The CETRAN Executive controls all of the
simulation software modules and performs the following simulation
controls : (a) creates the necessary initialization before simulation runs (b)
controls the timing based on the real time clock (c) Passes necessary
information to other modules (d) activates/deactivates the other modules.
The CETRAN Graphics Package generates a variety of displays for the
instructor and operator stations. This software allows the user to build
tabular, graphic, or trend displays, and to make the displays interact with
any simulated variables of the model. CETRAN Network Solvers are
provided for Thermal hydraulic (THLF) Network Models and Electrical
(ELNET) Network Models.
5
(b) The APROS — APMS (Advanced PROcess Simulator — Advanced
Paper and Pulp Mill Simulator) from Technical Research Institute of
Finland
6
inter-task communication within Microsoft Windows operating system,
CASSENG supports Dynamic Data Exchange (DDE), so that CASSBASE
can view block's details during debugging, and LabView can represent the
plant's simulated data via its graphical user interface. (3) LABVIEW:
LABVIEW for Windows is the user interface development environment
from National Instruments of Texas, U.S.A., with many value-added
features provided by C.T.I. for the purpose of developing user-friendly
graphical interfaces for simulator applications. Graphical screens with
buttons, icon symbols, pop-up dialogs, etc. are developed using LabView
for Windows. User interface symbols, and special-purpose modules are
provided in the form of LabView libraries. CTI has the custom LabView
libraries which provide power plant simulation related symbols, modules
for pop-up's, for performing Microsoft Windows Dynamic Data Exchange
(DDE) and TCP/IP communication with CASSENG, and more. Through
either DDE or TCP/IP, user interface screens request data from the
simulation running in CASSENG for display and monitoring. Likewise,
user input interactions such as specific control actions and setpoint changes
can be transmitted via button status changes sent to the running simulation
in CASSENG.
7
1.6 Simulation Fundamentals
1.6.1 Dynamic Simulation versus Steady–State Simulation
Process simulation is the use of mathematical models to develop a
realistic representation of the real process behavior, as measured by the
dynamic responses of the monitored process variables such as levels,
temperature, pressures, and materials compositions etc. Process simulation
can be divided into two categories: Steady–State and Dynamic.
8
example, the technique which applies to a single phase heat exchanger
differs considerably from those which apply to a two phase heat exchanger.
Hence, correct simulation problem definition is the first and the most
critical step in simulation. This is commonly known as "Model
Specification".
The simulation analyst must identify all of the "situations" which are to
be "simulated" by his/her model over the entire operating range, including
"situations" mitigated by component malfunctions (valves, pumps etc.). As
well, it is important to identify simulated conditions (single phase; two
phases) in an equipment such as heat exchanger so that the model will give
correct dynamic responses. However, as a caution, one must be able to
balance the needs for accurate models versus expensive model
development cost, and NOT to "over-specify" the requirements for the
model, resulting in paying for an expensive model feature for insignificant
training (or design) value.
9
· Depth of the details — specify whether every process, every equipment
(pump, valve, etc.), every instrument (pressure transmitter, temperature
transmitter etc.) shall be simulated. If simulation is required, should it
be a detailed simulation (down to the single electrical fuse); or
functional simulation is adequate?
10
Figure 2 — Typical Block Diagram of a control system associated with a
plant process
11
transmit, convert, and dissipate energy. Typical control volumes are
valves, pipes, generators, motors, vessels, and heat exchangers.
12
3
13
The following Table 1 summarizes the Quantity Variables, Flow Variables,
and Potential Variables used for pneumatic, hydraulic, thermal,
mechanical, and electrical systems. The units for these variables are listed
below:
1. The sum of all flow variables for any Junction is equal zero.
2. The sum of all potential variables for any closed loop in the system is
equal to zero.
For Junction n,
.....................(1.1)
where Fni is the Flow Variable from ith Control Volume going into (or
away) from Junction n.
.....................(1.2)
14
5
F1 - F2 - F3 = 0
F2 + F3 - F4 = 0
F4 - F1 = 0
Applying the Second Law for Inner and Outer Loops, respectively,
P1 + P2 + P4 = 0
P1 + P3 + P4 = 0
Note that the direction of flow should be taken into account in the flow
equations. In the potential equations, the potential variable is assumed to
be measured in the direction of the flow variable. Thus P1 is the potential
variable measured across junctions C to A.
15
Type of System First Law: Flow Second Law: Potential
Variable connected to Variable around a
a Junction loop
Pneumatic
(flow rates ) = 0 (pressures) = 0
Hydraulic
(flow rates) = 0 (heads) =0
Thermal
(enthalpy flows)= (temperatures)=0
0
Mechanical
(forces) = 0 (velocities)=0
Electrical
(currents) =0 (voltages) =0
· the relationship between the Flow Variable and the Potential Variable
Q = Quantity Variable
F = Flow Variable
P = Potential Variable
R = Resistance
C = Capacitance
16
The Resistance of a Process Control Volume is defined as the rate of
change of the Potential Variable with respect to the Flow Variable, and is
given by:
........................(1.3)
.........................(1.4)
..........(1.5)
The Resistance and Capacitance for the various physical systems are
introduced below as illustration of the concepts, without the detailed
derivation of the appropriate equations.
...............(1.6)
17
D = Inside Diameter of pipe (ft)
...........(1.7)
and .............(1.8)
Where V is the head; f is the friction factor; D is the inside diameter of the
orifice, and L is the equivalent length.
......................(1.9)
T = Temperature in (Degrees C)
Hydraulic systems also have two types of flow: turbulent flow where the
Reynolds number ( a function of geometry, flow rate and viscosity) is
greater than 4000; and laminar flow if Reynolds number is less than 2000.
The Quantity Variable of Flow Rate; Potential Variable is pressure.
18
..............(1.10)
...................(1.11)
and ......................(1.12)
P = Head, in (ft)
Q = Volume (ft**3)
....................(1.13)
19
Variable is Enthalpy, or commonly known as heat; the Flow Variable is
Enthalpy-Flow, or Heat Flow; the Potential Variable is Temperature.
...............(1.14)
............(1.15)
deg)
...................(1.16)
20
Variable is Force; Potential Variable is Velocity. According to Hooke's
Law, the spring constant is the inverse of the derivative of displacement
with respect to force.
........(1.17)
F = Force (lb)
........(1.18)
I = Current ( Amp)
L = dV/dI.......................(1.21)
21
(7) Summary
Turbulent flow
Turbulent flow
Forced Convection
22
2.0 Overview of Dynamic Simulation for Power Plant Process
2.1 Introduction
Formulating real time simulations of physical processes is, in some
respects, more difficult than formulating other quantitative models. This is
so because as of to date, there is no cogently expressed situation-
independent theory of real time simulation. Over the years, various
modeling techniques have evolved; the more successful survived and the
less successful have not. What is presented here then is a compendium of
techniques which are of necessity situation-oriented.
23
· Indications associated with water chemistry, which are not influenced
by operator actions may be made parameters adjustable by instructor
inputs.
· Momentum Balance
Where:
24
Normally, the level of liquid in the tank is the desired variable to be
calculated:
Where
L = level
= liquid density
where
Expanding the left hand side of equation 2.4, and substituting equation
2.1. After re-arranging, we have:
25
2.4 Rate of Change of Pressure in a Control Volume
Equation (2.1) can be re-written as:
or
26
Note that is a measure of fluid compressibility. For instance, for
Where:
If the pipe has a "check valve" ( i.e. only allows flow in one direction),
the flow rate will become zero whenever downstream pressure is greater
than upstream pressure. On the other hand, if flow reversal is possible for a
particular situation, equation (2.13) becomes:
If the fluid is gas, or vapor, the flow in the link connecting the two
Control Volumes also depends on density. Hence the following equation
may be used:
27
Where A is defined as Admittance and is expressed as:
Note that as the pressure difference goes to zero in equation (2.17), the
Admittance, A, will go to infinity. This can be prevented by limiting the
value of A to some maximum value when the pressure difference goes
below a predetermined valve.
Where:
= (P2 - P1)
Note that the pump's actual static head is related to the Design Static head
as follows:
28
Actual Static Head = Design Static Head * (Pump Speed)2 *
.............. (2.20)
Thus at 100 % rated pump speed, if there is no fluid density change, and
if the pump runs at 100 % efficiency (i.e. no degradation), then the Actual
Static Head is equal to Pump Design Static Head.
Figure 6
With this model, it can be seen that as the pump is started, the dynamic
head = (discharge pressure — suction pressure) is at its maximum. As
pump flow is increased, the dynamic head decreases accordingly.
29
2.9 Flow Equation in Link with Pump and Valve
Consider a "link" connecting two Control Volumes 1 & 2 has a pump
and a valve in series, the flow equation for this link can be expressed as:
................(2.23)
Where:
Where N.P.S.H = Net Positive Suction Head required for the pump to
function.
30
Rate of energy inflow from the incoming streams — Rate of energy outflow
from the outgoing streams + (or - ) Rate of external energy gain (or loss)
Mathematically,
Where:
Expanding the left side and substituting for dM/dt from equation (2.1),
after re-arranging, equation (2.25) will become:
31
2.14 External Heat Transfer Rate in Control Volume
The external heat transfer, Q, is normally calculated from:
for example:
32
same time, the condensate is drained out from the vessel at a certain rate. If
the drainage rate is less than the steam condensation rate, condensate will
accumulate inside the vessel to form level L.
Symbols Definition:
= Density of condensate
L = Level of Condensate
33
Where C is the Capacitance of the Vessel Steam Space. Wc is the steam
condensation rate.
It can also be shown that by using Ideal Gas Law, assuming constant
volume and constant temperature (recall Chapter #1 - on pneumatic
system)
34
3.0 Numerical Methods & Other Important Model Considerations
3.1 Introduction
The real time solution of a set of ordinary differential equations (ODE's)
in a digital computer based simulator requires the application of
appropriate numerical integration techniques to provide numerically stable
solutions, given a specific integration time step.
· Range of Validity
· Realism
· Reliability
35
A numerical solution of equation (3.1) can be obtained by
approximating the time derivative of X by a "backward difference" over
the time interval :
where:
X= X(t)
36
Again substituting (3.3) into (3.7), after re-arranging, we have obtained
a relationship between the error term at the current time and that at the past
time step:
(a)
(b)
· Valve actuators
· Thermocouples
37
Where X0 = initial value of X
= time constant
and
or,
38
Therefore, the Implicit Euler integration technique is unconditionally
stable, as long as the partial derivative of f is negative, which is always
met, as long as the equation is formulated properly.
For linear equations, the implicit integration can be applied easily, for
example, consider the following differential equation:
re-arranging:
Let
39
or,
It can also be shown that the linearization of f(X) does not affect the
stability conditions.
Where:
40
The same techniques for Explicit Euler and Implicit Euler can be
applied to the solution of the system of ODE's.
In the linear case, or in the non linear case after linearization, a typical
equation can be written as:
or,
Where:
41
The matrix solution to equation (3.37) becomes:
It should noted that if the size of matrix A is large, and if it contains very
large and very small values, (such as the case for solving mass and
momentum conservation equations (pressure and flow have small time
constants) and thermodynamics (temperature has large time constant)), the
matrix A is called "stiff" matrix.
Where:
M = mass of coolant
42
UA = overall heat transfer coefficient, multiplied by the average heat
transfer area
Where
This shows that for low flow, m could be small and less than 1/2. In
such case, the right hand side of equation (3.39) is negative, which means
the ratio of temperature differences is negative. This implies either T0 is
less than Ts, or Ti is less than Ts. In either case, this is physically
impossible.
(B) Realism
3.5.2 Reliability
The effectiveness of a simulator can be greatly reduced, if it often
interrupted by a failure of the simulator, either in the hardware, or in the
software model. An example of model failure is the case of division - a
number divided by another number, which suddenly becomes zero, causing
computer to "crash". To obtain a highly reliable software model,
unconditional stability of the numerical techniques may have to used to
solve ODE's, and the implicit algebraic equations. As well, model may
43
have to undergo "stress testing" to ensure proper behavior outside expected
range covered by the model. For example, in addition to testing a boiler
model under all operating conditions, it may have to be tested under
extreme conditions such as all water inventory evaporated, and see whether
the model can survive and still gives reasonable results.
· When all ODE's have been integrated once with a fixed time step ;
as well all other equations — algebraic and logic, are solved once, this
process is then called one ITERATION of the model equations.
Digital computer via a special "simulation executive" will execute the
model equations iteration by iteration continuously in time, unless
stopped by human intervention (in this case, we call it simulation
"freeze").
· The time interval between the start of two consecutive model iterations
is call "Iteration Interval ".
Should the current iteration take longer than sec. to complete, i.e. >
then the simulation executive will wait until the current iteration is
finished and will immediately execute the next iteration.
Now we can discuss the meaning of real time and other terms:
· REAL TIME:
44
· FASTER THAN REAL TIME:
AND
OR
Despite that fact that we set the Iteration Interval to be equal to the
integration time step (i.e. = ), but > .
45
Ô based Simulation Methodology & Simple Modeling Exercises
4.0 CASSIMÔ
· One final note, the completed solution to the simulation model you are
about to build is included to help you as a reference guide. It is
installed in the directory called 2TANKS.SOL. Please feel free to
reference this model solution in the CASSIM™ CD.
46
In block oriented modeling, a simulation model is described as a
sequence of interconnecting blocks, where each block represents a
component in the simulation. The output(s) of one block feeds into one or
more subsequent block(s) downstream. When the simulation runs, each
block is executed in sequence, producing output which is passed onto the
next connected block(s). This next block in turn operates based on its
inputs (and coefficients), and produce outputs for subsequent blocks
connected to its outputs. When all the blocks in the model has executed
once in sequence, the simulation has completed one iteration.
CASSIMÔ models are real time dynamic models which are based on
systematic application of first principles of mathematical modeling of plant
process systems.
47
LOGIC MODEL
Given:
gate 1
gate 4
gate 3
gate 2
MODEL
(COLLECTION OF BLOCKS CONNECTED IN A MEANINGFUL WAY)
ALGORITHM LIBRARY
TRANSMITTER CONTROLLER
PROCESS MODEL
Given:
VALVE 1
INFLOW1
VALVE 2
OUTFLOW1
TANK1
MODEL
(COLLECTION OF BLOCKS CONNECTED IN A MEANINGFUL WAY)
BLOCK BLOCK 3
1 Valve 1 (Inflow1) inflow 1 flowrate
port-area
(Valve 1)
BLOCK 5
(Tank1)
BLOCK 2 BLOCK 4 outflow 1 flowrate
(Valve 2) Valve 2 (Outflow1)
port-area
ALGORITHM LIBRARY
48
INPUTS, OUTPUTS & COEFFICIENTS
Block name: CTRL1_SIG Block name: VALVE1 Algorithm name: 108) CV_SIM Block name: INFLOW1
Cofs:
Stroking time 10
Plug type(1=linear;2=equal %) 1
Block name: CTRL2_SIG Block name: VALVE2 Algorithm name: 108) CV_SIM Block name: INFLOW2
Cofs:
Stroking time 5
Plug type(1=linear;2=equal %) 2
TAGNAME TRANSFORM
BLOCK NAME
Quick Quiz
3) Define iteration.
49
The seamless integration of these three tools provides a powerful
Microsoft Windows-based dynamic real time model development
environment, and a user interface development environment for producing
high fidelity, user-friendly simulation applications. The CASSIM
development system is also capable of producing stand-alone simulation
run time applications.
50
In summary:
51
· LabVIEW is used for user interface development. The user interface
screens communicate with the simulation model running in CASSENG
to allow user interactions with the simulation. LabVIEW Application
Builder is used to produce a stand-alone executable application
consisting of the user interface screens.
In the same plant, there may be many other valves similar to it. From a
simulation development standpoint, it would be cumbersome to model
each and every individual valve separately in the plant when they all have
similar characteristic. It would be more efficient to model a generic valve
only once, with some configurable parameters (or attributes as described
earlier), and store it as a template in a library. Subsequently, every valve
created will be based on the generic template stored in the library, with the
configurable parameters filled in to describe the particular instance of the
valve being modeled.
52
A CASSBASE simulation model consists of blocks which are
customized versions of generic algorithms from the algorithm library. We
are going to illustrate CASSIM™ based modeling using a simple two tank
process. Before we start, we should first describe the process.
· Each menu item will open up a drop-down menu. The first three drop-
down menu: Model, Algorithm, Block are used for managing models,
algorithms, and blocks respectively. The rest of the drop-down menus
are productivity tools to assist you in simulation modeling. The drop-
down menu items will be presented in more detail later. You will also
be using quite a few of these items to build the two tanks system.
· Let's open up a simulation model and begin building the two tanks
system. A model has already been started for you. From the top menu
53
bar, select the menu item Model. The drop-down menu list looks like
this:
· Select Open (hold down the left mouse button, drag the mouse pointer
onto Open, and let go of the left mouse button). This will open up the
file dialog for accessing a model.
· Using the directory list on the right hand side of the file dialog,
position yourself into the C:\CASSDEMO\2TANKS directory. You
should see a list of files ending with .dbf in the file list on the left hand
side of the file dialog.
· Click on any one of these files, for instance, blk.dbf, and press the OK
button in the top right hand corner of the file dialog.
· The file dialog will disappear. A few windows will flash by and you
should see the message:
54
· Select Open Algorithm (hold down the left mouse button, drag the
mouse pointer onto Open Algorithm and let go of the left mouse
button). This will open up a list of algorithms. The dialog looks like the
following:
55
4.6.1 The NON-MOVEABLE BLOCKS Algorithm
This algorithm is a CTI Generic Algorithm that comes with the full
CASSIM Simulation Development System. It is a very simple algorithm.
Its purpose is to bring together outputs (signals) from specific blocks in the
model to this one block for displaying the simulation values. You will see
later on when debugging the simulation that this block is very useful for
monitoring the running simulation values.
Click on the Edit Inps button to open the table for editing the input
descriptions. In the demo version of CASSBASE, you can not change the
contents of the table. In the first column (labeled PINNUM) are the input
numbers. In the second column (labeled DESC) are the descriptions for
each input. (When we developed this algorithm, there really wasn't much
to say about each input. Note that we simply use the input number as the
description). The last column is the data type of each input. The letter A
means analog, which is the same for all the inputs of this algorithm. An
56
analog value is a real numeric value. The data type can be digital, in which
case, the letter D is used. A digital value can have two possible values:
TRUE or FALSE, represented by 1.0, and 0.0 respectively.
Next press Edit Outs to open the table for browsing the output
descriptions. The columns are similar except there is one more column
called TRANSFORM. This column is where you would give short names
to the outputs. You will see later that these short names are used when
connecting blocks. Press Close Outs to close the table.
The Alg File Name field specifies where the FORTRAN source code for
this algorithm is stored. Note: If the Algorithm is CTI's generic algorithm,
source code file is not included. For this reason, if you click Edit Alg
button, you will not be able to see the source file.
Given an opening demand, the valve port area will open quickly and will
slow down when the opening approaches the demand. This behavior is
different than that simulated by algorithm #856 DEMO LINEAR VALVE.
The linear valve algorithm simulates a valve whose port area opens
linearly based on the demand signal.
57
The DEMO 1ST ORDER VALVE algorithm has one input, three outputs,
and one coefficient.
· Press Edit Inps to open the input description table. The one input is
Port Area Opening Demand. It is an analog value specified as a
percentage (decimal value between 0 and 1) opening.
· Press Edit Outs to close the input description table and open the output
description table. There are three outputs: Effective Port Area, Valve
Fully Opened, Valve Fully Closed.
1. Effective Port Area is the actual port area opening of the valve.
Although the demanded port area opening may be a certain value, it
takes some time before the valve drives to the desired opening.
This output simulates over time the action of driving the valve.
· Press Edit Cofs to close the output description and open the coefficient
description table. There is only one coefficient: Stroking Time. This is
use to set the stroking time characteristic of this particular valve. The
stroking time is defined as the amount of time it takes to drive the
valve from fully closed to 90% opened.
· Press Cancel to close the algorithm template.
The CTI Generic algorithm #108 Simple Control Valve incorporates
both of these valve port area opening types in addition to an Equal
Percentage opening type. The selection of which valve type to simulate is
done by specifying one of the coefficients in that algorithm. Feel free to
examine this and other algorithms in the algorithm template library.
Let's have look at how the two algorithm templates: 1) 452 NON-
MOVEABLE BLOCKS ALGORITHM, and 2) DEMO 1ST ORDER
VALVE are used as blocks in a model.
58
4.7 Opening a Block in the Model
One block of the two tanks system has been created for you.
· To open it and see its content, use the Block menu. From the top menu
bar, select the menu item Block. The drop-down menu list appears as
shown below:
· Select Open Block. This will open up a list of blocks defined in the
model in a block selection window.
· The block selection window appears as shown below:
· This window is used for browsing the list of blocks in the model.
Scroll up and down the list by pressing the PgUp and PgDn buttons.
Use the Top, Middle, Bottom buttons to jump to the respective section
of the model. The Search and Repeat buttons are used to locate to a
block beginning with a specific string. We will not be needing the
Search and Repeat functions.
· The blocks in the list are listed in execution sequence order (i.e.,
ordered by the numbers under the SEQ column). To change the display
order, press the Reorder button to display the blocks in order by the
59
numbers under the BLK column. Usually, the sequence numbers and
the block numbers are the same. After moving and deleting blocks in
the model, the numbers in the two columns may no longer match. The
simulation always execute in the order by sequence numbers.
· Since there is only one block in the model, the cursor is on this block.
Press Select to open this block in a block edit window for browsing.
· The block number and the sequence number are both 1. The name of
the block is DISPLAY. There is currently no block description text for
this block. The block is active. Inactive blocks will not be executed by
the simulation engine.
· The next four items concern the algorithm template which this block is
based on. The box following the title Alg Name is used for selecting
which algorithm template to base this block on.
· This block is currently based on algorithm #452 NON-MOVEABLE
BLOCKS ALGORITHM. The block inherits the characteristics of
algorithm 452 and thus has: 99 inputs, 99 outputs, and no coefficients.
· Press Edit Inps to see the input connection table. Many of the entries
have been filled in. Use the scroll bar to the right of the table to scroll
up and down to see the various connections. There are connections
which connect to other blocks' outputs that currently do not exist in the
model. You will be building them shortly.
· You can also move the highlighted cell within the table. Note that as
the highlighted cell moves from one row to another, the gray area
below the input table shows the algorithm description of the input that
the highlighted cell is on. Recall that the algorithm input descriptions
for the algorithm that this block is based on is just the input number.
60
· Press Edit Outs to open up the output table for browsing. Scroll up and
down the table to see other outputs.
· Press the Cancel button to close this block. When prompted whether to
save the changes or not, press No.
· You are now back at the block selection window. Press Close to close
this window.
61
· A line is provided for pumping liquid from the larger tank to the
smaller tank. Likewise, a similar line is available for transferring liquid
from the smaller tank to the larger tank. This is to facilitate the storing
of as much reserve liquid as possible utilizing the larger tank to meet
the demands of when the external sources are not available.
· The delivery of the liquid out of the smaller tank is controlled by a
level controller monitoring the tank level. The level controller adjusts
the valve on the outlet pipe to let liquid out so as to maintain tank level
at or below the controller set point.
· The following are some technical specifications of the system
described above. These specifications will be referenced later as you
build the simulation model of the system.
Tanks Tank 1 Tank2
Tank height 10 m 20 m
62
Physical Properties
Flow Rates Flow into Valve 1 Flow into Valve Flow into Valve
2 3
Control Specifications
· Valve 1's port area restricts the flow through the cross-section area in
the pipe and can be adjusted by controlling the stem position of the
valve to the desired position of opening. This control action is achieved
by manual control.
· A simple PID loop is used to control the outflow from tank 1 through
valve 3 by opening or closing valve 3 to maintain a tank level of
maximum 8 m. The PID loop can be put into manual mode to allow for
manual adjustment of the valve 3 port area opening.
63
4.8.3 Initial Conditions
· Tank 1 initial level = 5 m.
· Tank 2 initially empty.
· Valve 1, 2, 3 closed.
· Pump 1 & 2 stopped.
4.8.4 Flow Calculation Method
In this lesson, a flow is represented as a block using algorithm #860. In
fact, it should be computed as proportional to the square root of the
pressure drop between two points — part of the Hydraulic Network
Solution schema provided by CASSIM™. Since it will involve details
beyond the scope of this lesson, we are going to take a simpler approach
here. However, future Lesson will involve detailed Hydraulic Network
Solutions.
Therefore, in this lesson here, each flow block represents a flow path
which can carry a specified maximum flow rate. The same flow block also
allows the inclusion of one valve and one pump on the flow path to affect
the effective flow rate through the path. As depicted above, there are five
flow paths in the two tanks system.
Algorithm #860 calculates the effective flow in a path using the following
equation:
where
· Maximum Flow is in kg/s
64
· Valve Opening is % port area opening of the valve (in the range 0.0 to
1.0)
· Pump Power is % of full power (in the range 0.0 to 1.0)
You can follow this block diagram design to build a CASSBASE model
representing the two tanks system.
65
Press OK to create the new block.
· You are now looking at a new block in a window similar to the one you
used when browsing block 1. The difference is, this one is empty. You
will fill in the block entries, select an algorithm to base this block on,
and fill in the inputs, outputs, and coefficients.
· In the Block Name entry box, type in VALVE1 and press Enter.
· The next item is the Active item. This specifies whether or not this
block is executed when the simulation runs. By default, it is marked
with an "X" to indicate active. You can toggle the active/inactive
selection by clicking in the box.
· The next item is the algorithm selection for this block. By default, the
first algorithm in the list is selected. For this block, you want the block
to be based on algorithm #857 DEMO 1ST ORDER VALVE. To select
algorithm #857, click on the downward pointing arrow to drop-down a
list of algorithms. Scroll down and select the entry 857)
DEMO_VALVE.
· Press Edit Inps to open the input table for editing. There is only one
input for this block, the port area opening demand.
· For this simulation, the input comes from the user interface, where an
operator would control the valve by typing in a port area opening
demand. This demand value will be passed to the running simulation.
· The block where the input is passed to is the DISPLAY block. If you
had noticed before, input #3 of the DISPLAY block was reserved for
taking the port area opening demand for valve 1. Knowing that the
DISPLAY block transfers the value coming into input #3 to output #3,
VALVE1 block's port area opening demand input should be connected
to output #3 of the DISPLAY block.
66
· In general, an input can be specified in two different ways:
· Instead of typing in the output tag name to connect to, use the Connect
button to the right of the input table. For the VALVE1 block's port area
opening demand (input 1), move the highlighted cell inside the input
table to the cell to the TAGNAME column corresponding to the row
with pinnum = 1.
· Note that as the highlighted cell moves within the table, the gray area
below the input edit table show the algorithm description of the input
that the highlighted cell is on.
· Click the Connect button to the right of the input table. A block
selection window opens to allow for selection of the block to connect
to. Select block #1, the DISPLAY block and click the Select button.
This opens up a Select Transform dialog listing the output transforms
for the DISPLAY block. Select output #3 (XX03) and click the Select
button.
67
· Because this input gets its value from another block's output, there is
no need to enter any value into the DATAVALUE column.
· Press Edit Outs to open the output table. Values can be entered into the
DATAVALUE column to initialize the outputs produced by this block
before the simulation begins to run. In this case, the valve is initially
closed, therefore you have to configure the VALVE FULLY CLOSED
output to be TRUE (1.0).
· Press Edit Cofs to open the coefficient table. Position the highlighted
cell to enter the DATAVALUE for VALVE1's Stroking Time.
According to the specification given for VALVE1, the stroking time is
5 seconds.
· You are now presented with the Create New Block dialog. Press OK to
create block #3.
68
· Press Edit Inps to open the input table. There are two inputs.
· For input #1, the port area opening of the valve on the flow path, the
input comes from output #1 of the VALVE1 block. Position the
highlighted cell on the TAGNAME column for output #1. Press the
Connect button to the right of the input table. Select the VALVE1
block from the block selection window. Select the Effective Port Area
output (A01) transform.
· For input #2, the pump power of the pump on the flow path, there is no
pump on this flow path. To make the flow calculation work, enter a
constant value of 1.0 for this input. To do this, enter a constant tag
name for input #2 in the TAGNAME column. Type in $NO_PUMP.
Also enter 1.0 into the DATAVALUE column for input #2.
· Press Edit Outs to open the output table. Initially, this calculation
should show no flow. Therefore, leave the ACTUAL FLOW output's
DATAVALUE as 0.0.
· You are now presented with the Create New Block dialog. Press OK to
create block #4.
69
· If you refer to the block diagram earlier, you will notice that before we
can create a block for tank 1, there are other flow paths that are input to
tank 1.
· Make block #4 represent pump 1. In the empty block edit window, type
in PUMP1 as the block name. Enter a block description. The block
should be active. The algorithm that this block is based on is algorithm
#855 DEMO PUMP.
· There is only one coefficient for PUMP1: the pump run up time. Once
the pump is turned on, it takes some time for it to run up to full power.
The pump run up time is reflected in the pump power output.
According to the specifications given, PUMP1 runs up in 5 seconds.
Enter 5.0 in the DATAVALUE column.
· You are now presented with the Create New Block dialog. Press OK to
create block #5.
70
4.9.4 Creating Block #5— FLOW4
· The pump power output of the PUMP1 block is used in the calculation
of the flow path called FLOW4. Make block #5 this flow calculation.
Enter FLOW4 as the block name. Enter a block description. The block
should be active. The algorithm this block is based on is algorithm
#860 DEMO FLOW CALCULATION.
· Press Edit Inps to open the input table. There are two inputs.
· For input #1, the port area opening of the valve on the flow path, there
is no valve on this flow path. To make the flow calculation work, enter
a constant value of 1.0 for this input. To do this, enter a constant tag
name for input #2 in the TAGNAME column. Type in $NO_VALVE.
Also enter 1.0 into the DATAVALUE column for input #1.
· For input #2, the pump power of the pump on the flow path, the input
comes from output #1 of the PUMP1 block. Position the highlighted
cell on the TAGNAME column for output #1. Press the Connect button
to the right of the input table. Select the PUMP1 block from the block
selection window. Select the Pump Power output (A01) transform.
· Press Edit Outs to open the output table. Initially, this calculation
should show no flow. Therefore, leave the ACTUAL FLOW output's
DATAVALUE as 0.0.
· You are now presented with the Create New Block dialog. Press OK to
create block #6.
4.9.5 Creating Block #6 — PID Controller
· Let's choose another flow path that leads into tank 1 to model.
· Make block #6 the PID level controller. A PID controller has four
inputs:
71
1. Controller setpoint
2. Present value of the process variable being controlled
3. AUTO/MANUAL mode request
4. Manual control signal
· In the two tanks system, the PID controller is used as a tank level
controller to maintain the level of the liquid in tank 1 such that it does
not exceed a specific height. The specified height is used as the
controller setpoint. The PID level controller's output control signal
drives valve #3 which will open to relief some liquid inside tank 1 if
the level rises above the height setpoint.
· In the empty block edit window, enter PID as the block name, Enter a
block description. The block should be active. The algorithm this block
is based on is algorithm #859 DEMO PID CONTROLLER.
· Input #2, the present value of the process variable being controlled, is
the tank level output of tank 1. The Connect button can not be used for
specifying this connection because the block representing tank 1 has
not been defined yet. However, you can still specify the connection if
you know what the tag name of tank 1's level output is. Since you did
the analysis and design, a likely name for the tank 1 block is TANK1.
And knowing that the transform for the tank level output is A01, you
can enter into the input table's TAGNAME column TANK1/A01 for
input #2.
72
· Press Edit Outs to open the output table. The default values of 0.0 for
all of these outputs are just fine. Therefore, leave them as they are.
· Press Edit Cofs to open the coefficient table. There are four
coefficients. These are values for tuning the controller so that it makes
its automatic adjustments smoothly. Enter the values as shown below:
· You are now presented with the Create New Block dialog. Press OK to
create block #7.
· VALVE 3's port area opening demand signal comes from the output
control signal of the PID block. Press Connect and make the
appropriate selections to form the tag name PID/A01 for input 1.
· Press Edit Outs to open the output table. VALVE3 is initially closed,
therefore you have to configure the VALVE FULLY CLOSED output
to be TRUE (1.0). Leave the other output DATAVALUES.
73
· Press Edit Cofs to open the coefficient table. According to the
specification given for VALVE3, the stroking time is 5 seconds. Enter
this into the DATAVALUE column.
· You are now presented with the Create New Block dialog. Press OK to
create block #8.
· Enter FLOW3 as the block name. Enter a block description. The block
should be active. The algorithm this block is based on is algorithm
#860 DEMO FLOW CALCULATION.
· The inputs are: VALVE3/A01 for input 1, the port area opening of the
valve on the flow path, and $NO_PUMP with a DATAVALUE of 1.0
for input 2, the power of the pump on the flow path.
· The output is just fine. Press Edit Outs to check the table and press
Close Outs to close it.
· You are now presented with the Create New Block dialog. Press OK to
create block #9.
74
· Note that VALVE2 is a LINEAR valve. Choose algorithm #856
DEMO LINEAR VALVE. The inputs, outputs, and coefficients are the
same as a 1ST ORDER valve. The difference is the behavior of the
effective port area opening.
· For each new block, make sure you press the buttons: Edit Inps, Edit
Outs, and Edit Cofs at least once. CASSBASE checks this to try and
ensure you accept the contents of the tables before saving.
· The following tables list the values you should have in the blocks, use
these tables as reference in the unlikely event that you get stuck.
75
· When you are done saving block 12, you are presented with the Create
New Block dialog. Press OK to create block #13.
· The TANK1 block has four inputs. The tank in this demo has two
inlets and two outlets. The algorithm takes into account the flow rates
entering and leaving the inlets and outlets and adjusts the tank level
accordingly.
· Based on the P&ID diagram with the flow rates given, the flow rates
which are input to TANK1 are FLOW5/A01 and FLOW1/A01.
Likewise, the flow rates which deplete from TANK1 are FLOW4/A01
and FLOW3/A01. Make these connections. Connect FLOW5/A01 to
input #1, FLOW1/A01 to input #2, FLOW3/A01 to input #3, and
FLOW4/A01 to input #4.
· TANK1 has two outputs: tank level and tank overflow indicator. At the
start of the simulation, TANK1 should be set to a level of 5 meters. At
a level of 5 meters, there is no overflow. Therefore, the tank overflow
indicator having a 0.0 (FALSE) value is fine.
· Press Edit Cofs to open the coefficient table. The four coefficients for
TANK1 are: tank height (10 m), tank radius (0.25 m), liquid density
(1000 kg/m3), and height of outlet from the bottom of the tank (0.5 m).
Enter these values into the DATAVALUE column.
· You are now presented with the Create New Block dialog. Press OK to
create block #14.
· Set the initial TANK2 level to 0.0. Also note that there is only one
outlet in TANK2.
76
· The following tables list the values you should have in the TANK2
block, use this table as reference.
4.10 Hands-on Testing and Debugging the Simple Two Tanks Process
If you are unfortunate enough to run into some errors, the listing of
errors are captured in a text file. The Verify process will display a message
notifying you that there are errors and will ask you if you want to see the
error listing text file. Press OK to open this text file using MS Write. The
first MS Write dialog will prompt you whether you want to convert the text
file to Write format. Press No Conversion.
The most common error are hanging tag errors. A Hanging Tag is an
erroneous connection (hanging connection) specified under the
77
TAGNAME column in the input table. That is, the TAGNAME entry for
the input is suppose to be the output tag name of a block output to connect
to, but instead is a tag name which does not identify any output tag name in
the model. This usually occurs as a result of typing errors. If you used the
Connect button to specify connections, then you would not have hanging
tag problems.
Once you have corrected the errors (if any) ... YOU ARE DONE! All
the detailed work is done. Now comes the fun part of running the
simulation!
4.10.2 Exporting Model Data File for CASSIM™ Simulation Engine —
CASSENG
The model you have finished building is now ready for testing. To test
the simulation, you need the Simulation Engine — CASSENG and the
Debugger.
· To produce this file, select the menu item Export from the Model drop-
down menu. You are immediately presented with a file dialog to enter
the name of the compact data file to export to. Using the directory
selection list on the right hand side of the dialog, position yourself to
the \CASSDEMO\2TANKS directory. Then, in the entry field under
the File name:, type in the filename TANKS.DAT. Press OK.
· Let's start up CASSENG and open the compact data file. First
minimize CASSBASE. In your Program Manager, double-click on the
icon that looks like a train engine.
78
· These controls are not very exciting since you cannot see the model
values change as the simulation runs. There is a debug tool in
CASSIM™ that allows you to watch the inputs, outputs, and
coefficients of any block while the simulation is running.
The above only demonstrates that CASSENG and the model data file
can run by itself, but you do not see much going on. In order to visualize
the simulation data, you need to invoke the CASSIM Debugger:
· Click on the Model folder (yellow color). A "Model Open" Dialog box
shows up.
· Press "Model" button. A dialog box will appears, requesting your input
for the directory where the Model resides. Select
CASSBASE\CASSDEMO\2TANKS directory, and click on any .dbf
file in the directory. The model will be loaded.
79
· Press "Data" Button. A dialog box will appear, requesting your input
for model data file. Click "Tanks.dat" to select the model data file.
· Press "Done" Button. You will see the Debugger is loading up the
Model.
· Use the mouse to click on "Engine Control" and "Block List" Menu
items. (Note if you do not know where they are, move the mouse
cursor to any of the menu items, you will see a "tail-tips" behind it).
You will see the following dialog boxes.
80
· Double-click on the first block on the list "DISPLAY". The DISPLAY
block will be presented.
81
· The debug block window communicates with CASSENG via
Microsoft Window's Dynamic Data Exchange (DDE) facility. The
debug block window requests data from CASSENG and displays the
returned data. Likewise, input from the user is sent by the debug block
window to CASSENG and these inputs are integrated into the running
simulation by CASSENG.
Exercises:
If the simulation is frozen (i.e., counter not counting), press Run now.
The Iteration count will count up. Press Freeze to stop the simulation.
Press Step Iteration once and notice the iteration count go up by one.
82
Exercise:
TANK1/A01 .........
· This selects the input and it will be highlighted in "blue". This input is
connected to the block TANK1's A01 output. Try to display the debug
window for block TANK1. Loading Initial Conditions
· Note that the 2 tanks process that you are simulating can be at different
process conditions — that is, some pumps may be turned on and some
valves may be opened. So during debugging, you can store and reload a
particular process condition of interest. These initial conditions are
stored as IC files with *.ic file extension.
83
· Select 2 tanks_ic1.ic. This is a pre-stored IC for a normal condition for
the process.
· You will see the IC file is loaded via the "progress indicator".
Exercise:
· Let's tune some of the inputs of TANK1. First bring up the debug
window for TANK1.
· First run the simulation by pressing Run on "Engine Control". Note
that there is an inflow of 15 kg/s (input 1 connected to FLOW5/A01)
and an outflow of 5 kg/s (input 3 connected to FLOW4/A01). The tank
level output of the tank (TANK1/A01) is increasing because there is
more flow coming in than going out.
· Let's turn off the inflow at input 1. Double-click on the input 1 line.
This will bring up a Block Inputs Editor dialog as follows:
84
· Press Output Tag Name, you will see the Input Tagname FLOW5/A01.
· Put a $ sign in front of FLOW5/A01. In other words, you have
disconnected the output connection from Block FLow5 output A01,
and in its place, you have put a constant tag name called $FLOW5/A01
as input.
· Press "Constant Value" in the Block Input Editor dialog box. This will
switch the Data Value area to allow entering a constant value. Enter 0.0
in the entry box, and press Tune. The input is now tuned to have a
constant inflow of 0 kg/s.
85
· You will see "red" arrow pointing input $FLOW5/A01 in TANK1
Block. As well, you will see "red arrow" pointing the block TANK1,
indicating tuning changes are made to this block.
· Now press Run and notice that the tank level output is decreasing
because there is no inflow and still an outflow of 5 kg/s.
· Repeat the steps similar to above to tune the outflow input currently
connected to FLOW4/A01 to have a constant outflow of 0 kg/s. Run
the simulation and notice that the tank level output is constant because
there is no inflow and no outflow.
86
· All tuned values are used in the running simulation only and are not
stored in the model database if not saved. To save the tuned values into
the model database, single click anywhere in a block debug window
(but not on any input, output or coefficient line); then click the "right"
mouse button. A drop-down menu will appear, then select "save
block".
1. Before you save the block changes into the model, make sure the
changes are final and not temporary ones, because after you have
"saved" the block changes, you have to re-export the model data file
TANKS.DAT.
2. If the change is temporary, it is recommended to save the block tuning
changes in a new IC file. You can do that by clicking on the Debugger
Menu "Model" and select "Save IC" in the drop-down menu. You will
be asked to enter a new IC file name. As debugging progresses, you
can continue to make incremental "save" of the tuning changes into this
IC file, or if you wish, you can save into another new file, as the
situation requires. If there are numerous IC files, it is recommended
that you should label the IC files succinctly so that you may easily
recall what tuning changes those files have had, by looking at their
names.
3. After numerous testing, and if you are satisfied that the tuning changes
are final, you can select to update the model database with your IC file.
The way to do it is to first load the IC file as you have done previously
in "Loading Initial Conditions". After the simulation engine is loaded
with IC data, you select Debugger menu item " Save Model with
engine data", the tuning changes will be saved into the model database.
You have to re-export the model data file TANKS.DAT.
Exercise:
· Open TANK1 block and tune the inlet flow to be 10 Kg/s. As well,
tune the outlet flow to be 5 Kg/s.
· Save an IC file for this condition, and call this file tanks_ic2.ic
· Exit the Debugger, and make sure to exit the CASSENG as well.
· Now re-invoke the Debugger, load the original Model Data file
TANKS.DAT.
· Load the IC file tanks_ic.ic
· Restore all the tuning changes that you have made previously.
87
4.11 Running the Simulation Model with LabVIEW® Screen
C:\CASSBASE\CASSDEMO\2TANKS\TANKS.DAT.
· Open Program Manager and double-click on the icon that has two
tanks:
· This program sets the CASSENG simulation into run mode. The screen
appears as shown on the next page.
88
· At this point, you will notice there are THREE tasks running in
Windows: (1) CASSENG (2) Debugger (3) LabVIEW graphics
executable.
· The equipment that you modeled in CASSBASE are presented
graphically in the program window. There are two tanks, with pumps
on the pipes connecting them. Valves 1 and 2 control the inlets to tanks
1 and 2 respectively. There is a level controller on tank 1 controlling
the outflow through valve 3.
· Valves are green in color when they are fully closed, and red when they
are fully opened.
· Turn on valve 2 now by clicking the On/Off switch besides the valve.
Notice the pipe flowing into tank 2 fill up and the level of tank 2 rises.
The level of tank 2 is shown along the scale to the left of the tank and
as a numeric display to the left of the tank.
· Recall that tank 2 is the reservoir tank. Let's try to fill this tank's level
up to 10 m. To speed things along, it may help to use tank 1 to fill up
tank 2. First open the flow from tank 1 to tank 2 by turning on pump 1.
Turn on pump 1 now by clicking the On/Off switch besides the pump.
Notice the pipe connecting pump 1 to tank 2 fill up. This flow into tank
2 will speed up the fill of tank 2.
· At the moment there is only outflow from tank 1 and no inflow.
Eventually, this will deplete the liquid from tank 1. To add some
inflow, turn on valve 1 controlling the flow from an external source.
Turn on valve 1 now to fully opened by typing in 100.0 in the numeric
control besides valve 1 and press Enter. Recall that valve 1 can be
controlled to be partially opened.
89
· A couple of things to take note of here
1. Notice that although you typed in 100.0% opening, the valve still takes
time to stroke. In fact, the valve was modeled to take 5 seconds to
stroke for 0 to 90% opened. The 100.0% entered is the demand
opening only. There is no numeric display shown of the actual valve
opening. The color of the valve helps in this case. When the valve was
fully closed, the color was all green. When the valve is stroking, the
color is half red and half green. When the valve reaches 100.0%
opened, the valve color changes to all red.
2. Also notice that once the valve has completely opened to 100.0%, the
total inflow rate (with pump 2 off) is 5 kg/s as modeled. The total
outflow rate (with valve 3 off) is also 5 kg/s. This effectively makes the
level in tank 1 remain constant.
3. Is tank 2's level filled up to 10 m yet? Once the level gets close to 10
m, turn off first pump 1, then valve 2. This shuts off all flow into tank
2.
4. Time to squeeze some juice out of valve 3! That is, let's use the level
controller to operate valve 3 to let out the liquid in tank 1.
5. All along, valve 3 is being controlled by the level controller which is in
automatic mode. The setpoint is at 8 m. That means as long as the level
in tank 1 is below 8 m, then the controller will keep valve 3 close.
Once the level goes above 8 m, the controller will open up valve 3
accordingly to let out some flow from tank 1.
6. The following diagrams show how to use the level controller in
automatic mode and in manual mode.
· If the level in tank 1 has not reached 8 m yet, wait for it. If you are
impatient, you can lower the setpoint to a lower level. Try typing 7.0
into the setpoint (SP) numeric control. The controller will adjust its
control on valve 3 accordingly.
· Once the tank level has reached the setpoint, valve 3 opens. Turn off
valve 1 now to shut all flow into tank 1. Watch the controller slowly
close valve 3 as the level in tank 1 comes below the setpoint.
· Another method of letting some liquid out of tank 1 is to manually
operate valve 3. To do this, switch the controller to manual mode.
Switch to manual mode now by pressing on the MAN button. The
pointer will switch to point to MAN. You can now control valve 3 by
entering the desired opening in the numeric control under the label V3
Out (%). Open valve 3 completely now to 100.0% by typing in the
numeric control. Notice the flow out of valve 3 increase up to the
modeled maximum flow of 20 kg/s.
· At this point, you have interacted with most of the equipment modeled
using CASSBASE. Feel free to continue to run the simulation, trying
out different operating scenarios.
90
· NOW maximize Debugger, practice some tuning changes and observe
that the changes can be seen in the screen — e.g. change the tank
diameter.
91
momentum conservation equations (flow and pressure) are generally very
large while the eigenvalues caused by thermodynamics (temperature) are
small. Therefore a system of coupled differential equations describing
mass, momentum and energy dynamics of a large process could contain
wide spread of eigenvalues, and thus is usually referred to as “stiff” set of
equations. Stiff equations are usually solved by implicit methods due to
their inherent stability. However, standard implicit methods usually require
finding the Jacobian for the equation set during each iteration (Newton’s
method), as well as solving the resulting set of linear algebraic equations.
These solution techniques can provide solutions with high degree of
accuracy and stability; however, they cannot be solved “economically” in
real time.
Figure 8 illustrates the mass and momentum equations set for a typical
flow path containing equipments such as pumps and valves, and
demonstrates how these equations can be linearized to be solved by matrix
methods. By specifying nodes and links, the flow network problem can be
defined using modular block structured modeling, that is, each node block
and link block can be created by using generic node and link algorithm
respectively.
Hydraulic properties of the nodes and links blocks (e.g. node size, valve
size, pump dynamic head etc.) are either treated as coefficients or inputs to
these blocks. CASSBASE, the database engine of CASSIMÔ, has a
computerized network generation program which will automatically create
matrix elements for the network defined by the connected nodes and links
blocks.
92
HYDRAULIC NETWORK FLOW DIAGRAM
NHZ_X01 NHZ_N01 NHZ_N02 NHZ_X02
Tank 2
(X2)
NHZN01_N0
NHZX01_N01 NHZN02_X02
2
(N2) Valve 1
Tank 1 (N1)
EXTERNAL NODE
(X1) Pump
INTERNAL NODE
XXX_YYY NETWORK BLOCK NAME
Assuming incompressible flow
W X1N1 = K X1N1 (PX1 - PN1)1/2
W N1N2 = K N1N2 (PN1 + PDYH - PN2)1/2 Legend
W N2X2 = K N2X2VA (PN2 - PELV - PX2)1/2
W ij = flow from node i to j
Pi = pressure in node i
Linearizing equation by defining admittance VA = valve port area
AX1N1 = W X1N1 / (PX1 - PN1) PDYH = pump dynamic head
AN1N2 = W N1N2 / (PN1 + P DYH - PN2) PELV = elevation head
AN2X2 = W N2X2 / (PN2 - PELV - PX2) Aij = admittance between node i and j
Kij = conductance between node i and j
Ci = capacitance of node i
Performing mass balance on N1 and N2; assuming constant density ,t = time increment
CN1(dP N1/dt) = W X1N1 - W N1N2 P' i = Pi at previous iteration
CN2(dP N2/dt) = W N1N2 - W N2X2 Mij = matrix element ij
Si = source element i
Using backward Euler; and substitute admittance
CN1((PN1 - P'N1)/, t) = AX1N1(PX1 - PN1) - AN1N2(PN1 + P DYH - PN2)
CN2((PN2 - P'N2)/, t) = AN1N2(PN1 + PDYH - PN2) - AN2X2(PN2 - PELV - PX2)
where
( M11
M21
M12
M22) • ( ) =( )PN1
PN2
S1
S2
93
CASSIM™ Network naming convention, we give this network a 3 letter
pre-fix NHZ. Hence, all nodes and branches block name will have this pre-
fix identifier. For instance, N1 will be labeled as NHZ_N1; the branch
connecting Tank #1 and Node N1, is labeled as NHZX1_N1.
(a) The flow in each link is computed using momentum balance for each
link.
(b) The pressure in each node is computed using mass balance equation for
each node.
(1) Users create node blocks and link blocks. The node block output
contains node pressure values; whereas the link block output contains
link flow values.
(2) After users have constructed all the node and link blocks with the
equipment (pumps and valves), and with proper initializations, users
"generate" the network using CASSIM.
(3) After "network generation", the CASSIM utility has already built the
matrix equation similar to that shown in the above diagram.
(4) In the first iteration, for each link, given the initial flow, initial pressure
in upstream node i and downstream node j respectively, admittance is
calculated in each link block (see equation above); and hence all the
M's and S's in the matrix elements can be calculated. When these
calculations are all done, the Network Solver inside CASSIM solves
the matrix equation in real time to give the pressure values for all the
nodes. Once the pressures of the nodes are solved, they update all the
node blocks' pressure output.
(5) After the pressure update in the node blocks is completed, flow can
then be calculated from executing the link flow block algorithm (see
equation above).
(6) Now it is ready for the next iteration with exactly the same iteration
scheme as described in (4) above.
Note:
(1) The equations formulated below for the hydraulic network pressure and
flow calculations are only good for incompressible flow. For
compressible flow (steam, air systems etc.), users have to consider the
use of "variable" conductance in the link, in order to use the same set
of incompressible flow equations. See explanation below.
94
heat exchanger) is done outside the hydraulic network blocks
boundary. It can be achieved by using "mixing node" algorithm, in
which the output temperature or enthalpy of the mixing node is
calculated from the inlet streams' temperatures or enthalpies, assuming
perfect mixing. The effect of temperature change can lead to a change
in the flow in the link due to the change in the fluid density (input # 7
of the link algorithm).
0.5
Flow = K c × V1 × Pup 2 - Pdown 2
Flow = K c × V1 × Pup
where K c = link conductance; V1 = Valve port area in the link between the
upstream and downstream nodes; Pup , Pdown = upstream and
downstream node pressure respectively.
95
is then connected to the link block as input. In other words, this
technique “allows” the link momentum equation to be an
“incompressible” flow type; however in this instance, the link
conductance should not be fixed, but rather a “variable” one and is
calculated in each iteration from the compressible momentum
equations using back-substitution. The application of this technique
allows the same numerical solution methodology formulated for
“incompressible” systems (Fig. 8) to be applicable for “compressible”
flow systems, with the additional computation of “variable link
conductance”.
· Internal Nodes = Control volumes whose pressures will vary with time.
(A) Pressure drop can be calculated using Darcy’s Law of the exact
piping configuration is known.
· Sharing the overall pressure drop between the branches and the
valves.
96
4) Create Network Control Block with Blockname NH?_CTRL and move
it ahead of the network marker block NH?_END
Why?
Hint
Why?
If a variable has not been assigned a value, the memory space occupied
by it contains garbage data. When you reference this undefined
variable, the garbage data will be interpreted as legitimate data and will
cause unexpected results.
97
TRUE but does nothing if the condition is FALSE. A reference to that
variable just after the IF-THEN-ELSE structure will return the assigned
value if the condition was TRUE and will return garbage data if the
condition was FALSE.
Examples:
123.764 Incorrect
.4352344 Incorrect
1423.34E12 Incorrect
+345.E-4 Incorrect
-.4565788E3 Incorrect
2E6 Incorrect
1234. Incorrect
123.764D0 Correct
1423.34D12 Correct
+345.D-4 Correct
-.4565788D3 Correct
2D6 Correct
1234.D0 Correct
Examples:
TLOAD = 0 Incorrect
TLOAD = 0.0 Incorrect
TLOAD = 0.0E0 Incorrect
98
TLOAD = 0.0D0 Correct
5. When writing a mathematical equation, Use the same data type in all
parts of the calculation to ensure that precision is not lost.
dY/dt = AY + C
99
5.0 Reactor Simulation Model Development
· Macroscopic Cross-section:
The macroscopic cross-section 5 is also proportional to nuclei density
of the medium; hence
5 = N • I ............................................(5.2)
where N = number of nuclei per cm
I = microscopic cross-section for the medium in cm2.
· Fission Cross-Section:
Since different reactions occur with different probabilities, they will be
associated with different cross-sections:
(a) fission cross-section å f ;
(b) absorption cross-section å a ;
(c) elastic scattering cross-section å s ;
(d) inelastic scattering cross-section å i .
100
Since we are primarily dealing with nuclear fission which causes
nuclear chain reaction; hence the number of nuclear fissions per unit
volume per sec is B • å f .
· Reactor Power:
Given a reactor with fissile material with volume V cm3, fission cross
section å f (fission reactions/(neutron.cm), neutron flux B
(neutrons/(cm2.sec)), the number of fissions occurring per sec in such
reactor is V • å f • B .
Note: the average energy released per fission is 200 MeV. Note 1 eV =
1.6 ´ 10 -19 watt sec (Joule), so 200 MeV = 3.2 ´ 10 -17 watt sec.
Of the (typically) 2.5 neutrons produced at fission, only one goes through
this cycle to maintain chain reaction. The other 0.5 are lost by capture or
leakage. We start with N thermal neutrons in the reactor, and go round the
cycle as follows:
2. Not all the neutrons absorbed by the fuel will cause fission, because
some will be lost by radiative capture. Let a fraction “a” cause fission
in U-235, the remaining fraction (1 - a) being lost by radiative capture.
Therefore, we get afN fissions.
3. Each fission will produce n fast neutrons (note for U-235, n = 2.43).
Therefore, an fN fast neutrons are produced by the original N thermal
101
neutrons. The product an is usually replaced by one symbol h so that
h fN fast neutrons are produced by U-235 fissions. The factor h is
called the thermal fission factor. It represents the number of fast
neutrons produced from each thermal neutron captured in the fuel;
whereas n represents the number of neutrons produced from each
thermal neuron actually causing fission.
7. Some thermal neutrons will be leaked out from the reactor. Let Pth be
the thermal neutrons non-leakage probability, then the number of
thermal neutrons which do not leak from the reactor is
A h f Pf p Pth N
· For very large reactors, leakage is very small, hence Pf ~1 and Pth ~ 1,
we get
K¥ = A h f p ....................(5.6)
102
K ¥ is called the infinite multiplication factor, meaning that it is the factor
that would apply if the core in question were made so large that the
neutron leakage became negligibly small. Hence, this called the “four-
factor formula”
· Loss rate due to leakage: Fick’s law states that assuming that scattering
is isotropic; neutron flux is a slowly varying function of position; no
source in the medium, then neutrons flow from regions of higher
densities to lower densities, i.e. the neutron current (number of
neutrons crossing unit area in unit time in the normal direction is
related to the neutron flux -
L
J = - D × gradB = - D × Ñ 2B ............(5.7)
dn
= 0 = K eff × B × S f + S - ( - D × Ñ 2B + B × S a ) , or
dt
103
· The diffusion equation is the basic equation which will be used by
reactor engineers in designing the reactor. However, since diffusion
theory makes the basic assumption that all neutrons have the same
thermal energy, the multiplication factor or critical size calculations are
therefore not very accurate. Application of the slowing down theory
improves the result somewhat since the neutrons are classified into
epithermal (slowing down region) and thermal (diffusion region). Even
this two energy group representation might not be adequate for actual
design of power plant. In order to circumvent this problem, multigroup
diffusion theory equations have been developed. In practice the total
neutron energy spread from 10 MeV to thermal energy (0.025 eV) may
be divided into 100 groups. To account for leakage and other spatial
effects, one has to divide the reactor into some spatial mesh points at
which the diffusion equations are solved by applying the boundary
conditions. For example, if there are 100 spatial points in each
direction, then there are 106 points in the reactor at which one has to
solve the diffusion equations. Since we are primarily interested in the
dynamic behavior of the reactor, we conclude the steady state design
topic discussion here.
la 1
l0 = = ................(5.9)
v vS a
· For a finite reactor, the mean neutron lifetime l = l0 × Pth , where Pth is
the thermal neutrons non-leakage probability.
l × N = L × K eff × N
l
or L= ..................(5.10)
K eff
104
· The reactivity, denoted by ,K or H , is defined in some literature by
the equation
K e - 1 K ex
H = ,K = = ...............(5.12)
Ke Ke
dN ( K e - 1) × N N × DK
= = .............. (5.13)
dt l L
Since neutron flux and reactor power are both proportional to neutron
density, and they will have the similar exponential behavior as the neutron
population growth equation above.; that is
æ DK ö
ç ×t ÷
è L ø
· Neutron density : nt = n0 × e ..........(5.15)
æ DK ö
ç ×t ÷
· Neutron flux: Ft = F0 × eè L ø
........(5.16)
æ DK ö
ç ×t ÷
è L ø
· Neutron power : Pt = P0 × e ........(5.17)
5.1.7 The Reactor Period for Reactor with Prompt Neutrons Only
If we define the reactor period T as the length of time required for the
power (or flux or neutron density) to change by a factor e. That is:
L
T= ...................(5.18)
DK
105
The previous equations now become:
ætö
ç ÷
nt = n 0 × e è T ø
ætö
ç ÷
èTø
Ft = F0 × e ...............(5.19)
ætö
ç ÷
èTø
Pt = P0 × e
You can see that the shorter the reactor period, the faster the power
changes will be. Conversely, should the reactor period be infinite, the
reactor would be critical.
L 0.001
T= = = 2 sec
DK 0.0005
0.001
T= = 0.2 sec.
0.005
· These two examples show how rapid the power increases would be -
even for small reactivity changes of a order of a mk - if all the
neutrons were prompt neutrons.
106
this weighted average lifetime plus the slowing-down time and the
diffusion time in the moderator. The latter is typically 0.001 sec., and so
we arrive at J = 0.085 + 0.001 = 0.086. Hence, when compared with
lifetime of prompt neutrons of 0.001 sec., delayed neutrons increase the
average lifetime of all neutrons by 86 times.
J 0.086
T= = = 17.2
DK 0.005
1
Pt
The relative power increase per second is = e 17.2 = 105986
. ; that is 5.9
P0
% per sec. corresponding to +5 mk increase in reactivity. Therefore, the
power increase is much less rapid than with prompt neutrons alone.
Total = 0.65 %
107
5.2 Point Kinetic Model — Delayed Neutrons
Given total m groups of fission precursors, the decay of the ith group
precursor is described by the following differential equation:
dCi n
= bi × K e × - li × Ci i = 1..... m .......(5.20)
dt l
n
>i K e = Neutrons born to the i-th group precursor
l
dN
= Rate of birth of Prompt Neutrons - Rate of loss of Neutrons + Rate
dt
of birth of Delayed Neutrons ........(5.21)
n
Rate of birth of Prompt Neutrons = K e (1 - b ) .......(5.22)
l
N
Rate of loss of Neutrons = ..........(5.23)
l
m
Rate of birth of Delayed Neutrons = S li × Ci ........(5.24)
i =1
dn n m
= ( K e (1 - b ) - 1) × + S li × Ci ............(5.25)
dt l i =1
dn n m
or = ( K ex - K e b ) × + S li × Ci ...........(5.26)
dt l i=1
dn DK - b m
or = × n + S li × Ci .......................(5.27)
dt L i =1
108
dn DK - b m
= × n + S li × Ci ................(5.28)
dt L i =1
dCi n
= bi × - li × Ci i = 1...... m .........(5.29)
dt L
dn DK - b
= n + l × C ...............(5.30)
dt L
dC b
= n - l ×C ................(5.31)
dt L
DK - b
(s - ) × N ( s ) = l × C( s ) + N 0 .........(5.32)
L
b
( s + l ) × C( s ) = × N ( s ) + C0 ...........(5.33)
L
Rearranging the second equation and substituting expression for C(s) into
the first equation, we have:
DK - b l ×b l × C0 b/L
(s - - ) × N ( s) = N 0 + = N 0 (1 + ) ....(5.34)
L L( s + l ) s+l s+l
dC b
because = 0 for t < 0, hence C0 = N .......(5.35)
dt l ×L 0
b
s+l +
Hence N ( s ) = N 0 L ........(5.36)
DK - b l × DK
s + (l -
2
)s -
L L
l × L - DK + b l × DK
The characteristic equation is: s2 + ( )s - =0
L L
.....(5.37)
109
If l × L << > , the characteristic equation for single precursor model
reduces to:
b - DK l × DK
s2 + ( )s - = 0 ................(5.38)
L L
1 ( b - DK ) é 4LlDK ù
s1, 2 = ê- 1 ± 1 + ú .......(5.39)
2 L ë ( b - D K )2 û
1 1 2 1× 3 3
(1 + x ) = 1 + x- x + x +......
2 2×4 2×4×6
We can keep the first two terms as approximation if x <<1. For small ,K
and the fact that l × L << > , this approximation is valid. Hence, the roots
equation can be approximated as:
1 ( b - DK ) é 2LlDK ù
s1, 2 = ê - 1 ± (1 +
( b - DK )2 úû
) ........(5.40)
2 L ë
lDK
that is s1 = ; ..........(5.41)
( b - DK )
1 é ( b - DK )2 + LlDK ù
s2 = - ê ú ........(5.42)
Lë ( b - DK ) û
( > - DK )
»- since l × L << > ........(5.43)
L
110
L
TP = ...............(5.45)
DK
Case (1) if ,K << > , the reactor period is dominated by first exponential
term :
( b - DK ) 1 b
T= » ...........(5.46)
lDK DK l
T = 28 sec.
L 1 b
<< .............(5.47)
DK DK l
1 m
b
ål
i
T= ............(5.48)
DK i =1 i
L L
T= » ..........(5.49)
( > - DK ) DK
This is the same as the prompt neutron model. This implies that if
,K >> > , the beneficial effect of the delayed neutrons has been lost, and
prompt neutrons dominate the response.
111
2.2
,k=0.003
1.8
N/No
1.6
,k=0.002
1.4
,k=0.001
1.2
1
0 0.2 0.4 0.6 0.8 1 1.2 1.4
Time
112
Dollars: one dollar is defined as the magnitude of ,K which is equal to
> . Thus for Uranium 235 fueled reactor, one dollar is equal to 0.0065,or
+6.5 mK. When the reactivity is one dollar, the reactor is critical on
prompt neutrons alone (since dn/dt = 0; if ,K = > ). Normally reactor
operates near delayed critical with ,K near zero.
1b 0.0065
DK » = = 0.0000234 = 0.023mK
T l 3600 × 0.077
Note: for inherent safety of the reactor, if the period is smaller than 10
seconds, the reactor control is designed to automatically shutdown the
reactor.
dCi bi * N FLUX
= - li * Ci
dt N LIFE
Data:
b1 = 0.000329
b2 = 0.001003
b3 = 0.000856
b4 = 0.002319
b5 = 0.000587
113
b6 = 0.000206
l1 = 0.01261 sec-1
l2 = 0.03051 sec-1
l3 = 0.12 sec-1
l4 = 0.3209 sec-1
l5 = 1.262 sec-1
l6 = 3.239 sec-1
114
7. Copy the new CASSIM.DLL to the CASSENG directory.
10. Create a block for calculating the reactor neutron flux using the custom
algorithm ALG802. Remember to initialize the neutron concentrations
by depositing their values in the output section.
12. Try to maneuver the reactor power to 80% by depositing the liquid
zone reactivity change to no more that -0.1mK.
13. Shutdown the reactor with large negative reactivity, and observe the
reactor power decay.
dCi bi * N FLUX
= - li * Ci
dt N LIFE
115
Use Backward Euler's to approximate the above equations:
dCi bi * N FLUX
= - li * Ci
dt N LIFE
Ci - Ci bi * N FLUX
= - li * C i
Dt N LIFE
re - arranging C i
bi * N FLUX * Dt
Ci * (1 + li * Dt ) = + Ci
N LIFE
bi * N FLUX * Dt
+ Ci
Ci = N LIFE
(1 + li * Dt )
6
( Dk - å bi )* N FLUX 6
dN FLUX
= i =1
+ å li * Ci
dt N LIFE i=1
6
( DK - Sbi ) × N FLUX
N FLUX - N FLUX i =1
6
= + Sl i × Ci
Dt N LIFE i =1
6
( DK - S bi ) × N FLUX é bi × N FLUX × Dt ù
i =1 ê + Ci ú
N FLUX - N FLUX 6 N LIFE
= + Sl i × ê ú
Dt N LIFE i =1 ê
(1 + l i × Dt ) ú
êë úû
116
Multiply both sides by ,t ,
6
( DK - S bi ) × N FLUX × Dt é bi × N FLUX × Dt ù
i =1 ê + Ci ú
6 N LIFE
N FLUX - N FLUX = + Sl i × Dt × ê ú
N LIFE i =1 ê (1 + l i × Dt ) ú
êë úû
defining Ai and Bi
li * Dt * Ci
Ai =
(1+ li * Dt )
li * bi * Dt
Bi =
(1+ li * Dt )* N LIFE
6
( Dk - å bi ) * N FLUX * Dt 6
N FLUX - N FLUX =
i =1
+ å ( Bi * N FLUX * Dt ) + Ai
N LIFE i =1
6
( Dk - å bi )* Dt 6 6
N FLUX * 1 - i =1
- å ( Bi * Dt ) = N FLUX + å Ai
N LIFE i =1 i =1
6
N FLUX + å Ai
N FLUX = 6
i=1
( Dk - å bi ) 6
1 - Dt *
i =1
+ å ( Bi )
N LIFE i=1
117
C
C ALGORITHM NAME: SIMPLE POINT REACTOR KINETICS &
C REACTIVITY
C
C AUTHOR DATE
C
C-----------------------------------------------------------
C
C THIS ALGORITHM CALCULATES CANDU REACTOR KINETICS
&
C REACTIVITY
C
C-----------------------------------------------------------
C
C INPUTS: 2
C
C INP(1) = DKXE, reactivity change due to xenon(K)
C INP(2) = DKLZ, reactivity change due to liquid zone(K)
C------------------------------------------------------------
C
C OUTPUTS: 9
C OUT(1) = NFLUX, neutron flux (normalized)
C OUT(2) = DK, overall reactivity change (K)
C OUT(3) = BETA1, Total delayed neutron fraction
C OUT(4) = C(1), neutron concentration, group 1(normalized)
C OUT(5) = C(2), neutron concentration, group 2(normalized)
C OUT(6) = C(3), neutron concentration, group 3(normalized)
C OUT(7) = C(4), neutron concentration, group 4(normalized)
C OUT(8) = C(5), neutron concentration, group 5(normalized)
C OUT(9) = C(6), neutron concentration, group 6(normalized)
C------------------------------------------------------------
C
C COEFFICIENTS: 13
C COF(1) = TNEUTRON, mean neutron lifetime(sec)
C COF(2) = BETA(1), delayed neutron fraction, group 1(normalized)
C COF(3) = BETA(2), delayed neutron fraction, group 2(normalized)
C COF(4) = BETA(3), delayed neutron fraction, group 3(normalized)
C COF(5) = BETA(4), delayed neutron fraction, group 4(normalized)
C COF(6) = BETA(5), delayed neutron fraction, group 5(normalized)
C COF(7) = BETA(6), delayed neutron fraction, group 6(normalized)
C COF(8) = LAMDA(1), delayed neutron decay constant, group 1 (1/sec)
C COF(9) = LAMDA(2), delayed neutron decay constant, group 2 (1/sec)
C COF(10) = LAMDA(3), delayed neutron decay constant, group 3 (1/sec)
C COF(11) = LAMDA(4), delayed neutron decay constant, group 4 (1/sec)
C COF(12) = LAMDA(5), delayed neutron decay constant, group 5 (1/sec)
C COF(13) = LAMDA(6), delayed neutron decay constant, group 6 (1/sec)
C------------------------------------------------------------
C
C PROGRAM DECLARATIONS
118
C
IMPLICIT NONE
C
INTEGER I
DOUBLE PRECISION INP(*), OUT(*), COF(*), DT
C
DOUBLE PRECISION DKXE,DKLZ
DOUBLE PRECISION NFLUX,DK,C(6), BETA1
DOUBLE PRECISION BETA(6),LAMDA(6), TNEUTRON
DOUBLE PRECISION A(6), B(6)
C
C INPUTS, OUTPUTS, AND COFS SETUP
C
DKXE = INP(1)
DKLZ = INP(2)
C
NFLUX = OUT(1)
DK = OUT(2)
BETA1 = OUT(3)
C(1) = OUT(4)
C(2) = OUT(5)
C(3) = OUT(6)
C(4) = OUT(7)
C(5) = OUT(8)
C(6) = OUT(9)
C
TNEUTRON = COF(1)
BETA(1) = COF(2)
BETA(2) = COF(3)
BETA(3) = COF(4)
BETA(4) = COF(5)
BETA(5) = COF(6)
BETA(6) = COF(7)
LAMDA(1)= COF(8)
LAMDA(2)= COF(9)
LAMDA(3)= COF(10)
LAMDA(4)= COF(11)
LAMDA(5)= COF(12)
LAMDA(6)= COF(13)
C
C-------- TOTAL REACTIVITY CHANGE (K)
C
DK = DKLZ + DKXE
C
C
C-------- TOTAL DELAYED NEUTRON FRACTION
C
119
BETA1 = BETA(1) + BETA(2) + BETA(3) +
& BETA(4) + BETA(5) + BETA(6)
C
C--------DEFINE Ai,
C
DO 1 I = 1, 6
A(I) = LAMDA(I) * C(I) * DT / ( 1D0 + DT * LAMDA(I))
1 CONTINUE
C
C--------DEFINE Bi
C
DO 2 I = 1, 6
B(I) = LAMDA(I) * DT * BETA(I) /
& (TNEUTRON * ( 1D0 + DT * LAMDA(I)))
2 CONTINUE
C
C-------- NEUTRON FLUX CALCULATION (NORMALIZED)
C
NFLUX =DMAX1(((NFLUX + A(1) + A(2) + A(3) + A(4) +
& A(5) + A(6)) /
& (1D0 - DT * ((DK -BETA1)/TNEUTRON +
& B(1) + B(2) + B(3) + B(4) + B(5) + B(6)))),0.0D0)
C
C--------DELAYED GROUPS CONCENTRATION
C
DO 3 I = 1, 6
C(I) = (C(I) + (DT * BETA(I) * NFLUX/TNEUTRON)) /
& ( 1D0 + DT * LAMDA(I))
3 CONTINUE
C
OUT(1) = NFLUX
OUT(2) = DK
OUT(3) = BETA1
OUT(4) = C(1)
OUT(5) = C(2)
OUT(6) = C(3)
OUT(7) = C(4)
OUT(8) = C(5)
OUT(9) = C(6)
C
RETURN
END
120
5.4 The Effects of Neutron Sources on Reactor Kinetics
The point kinetic reactor model equations in the previous chapter are all
based on the assumption that there are no sources of neutrons in the reactor
other than the (induced) fission. This assumption is not entirely true. In this
chapter, we shall consider the effects that neutrons sources other than
fission have on power changes.
1 H 2 + C = 1 H 1 + 0 n1 ................. (5.50)
This photoneutrons source will persist after the reactor has been
shutdown and the prompt and delayed neutrons have died down.
121
S 3 = ( S 0 LK e + S 0 L ) K e + S 0 L = S 0 L(1 + K e + K e2 ) ........... (5.51)
Hence after mth generation, the total neutron population due to the source
will be:
S 0 L(1 - K em )
or , Sm = ........................(5.53)
(1 - K e )
t
since m = , therefore the number of source neutrons at time t is given
L
by:
t
S0 L(1 - K ) L
St = e
..................(5.54)
(1 - K e )
S0 L
St = ..............(5.55)
(1 - K e )
Case (2) K e =1: In this case, S 0 L neutrons will be added to the neutron
population in every generation, so that after a time t equivalent
t
to m = generations, the number of source neutrons will be:
L
t
St = S0 L = S0t ...........(5.56)
L
Case (3) K e > 1: For a supercritical reactor, the neutron population will
increase very rapidly. The fission neutrons already there will
increase at rate governed by the reactor period and the source
neutrons will multiply as given by equation (5.54)
122
total flux will be constant, so will the neutron density. Suppose there
are n neutrons in any generation, they will related by :
nK e + 10-6 n = n ....................(5.55)
S0 L
St = .........(5.55)
(1 - K e )
We can therefore say that the final power level, P¥ , will be given by:
Ps
P¥ = ......................(5.56)
1 - Ke
They will decay with a half life of 18 days and will therefore set limit
on the length of time during which the reactor may still be restarted
without special instrumentation.
123
Table 4 - Photoneutron Groups for a typical HWR.
Total=28 x10-3
%
dCi n
= b i × - l i × Ci for i = 1, ... 15
dt L
Data:
124
b9 = 0.000013 “
b10 = 0.0000111 “
b11 = 0.0000061 “
b12 = 0.0000016 “
b13 = 0.0000008 “
b14 = 0.0000002 “
b15 = 0.0000002 “
l1 = 0.01261 sec-1 Delayed Neutron Group Decay Constant
l2 = 0.03051 sec-1 “
l3 = 0.12 sec-1 “
l4 = 0.3209 sec-1 “
l5 = 1.262 sec-1 “
l6 = 3.239 sec-1 “
l7 = 0.2778 sec-1 Photoneutron Group Decay Constant
l8 = 0.01754 sec-1 “
l9 = 0.00556 sec-1 “
125
1. Modify the Algorithm # 802 to incorporate these changes. Be sure to
include additional inputs, outputs and coefficients
1. Long term effects which consist of fuel burn up, conversion of fertile
material into fuel material (build up), and the accumulation of fission
products. Time scale is in months.
3. Short term effects which arise due to temperature changes in fuel, heat
transport fluid (void) and moderators, The time scale for these short
term effects are of the order of minutes.
126
· The rate at which this occurs depends on the neutron flux, because the
rate dC/dt at which U-235 is destroyed by neutron capture is given by :
dC
= - CI aB .....................(5.57)
dt
B = neutron flux
U 238 + n ® U 239 + C
¯
93 Np 239
¯
94 Pu239
· The conversion ratio is the number of Pu-239 nuclei formed for each
U-235 nuclei consumed. In typical CANDU Heavy Water reactor, it is
about 0.8.
· The rate of change in the concentration Ci of the ith isotope in the fuel
in a flux of B is in general given by:
dCi
= Ci -1s i -1a i-1f + C j l j - Cis if - Ci l i .................(5.58)
dt
The first term Ci-1I i-1= i-1B is the production of the ith isotope by
neutron capture in the (i-1)th isotope.
127
The second term C j l j is the production of ith isotope from radioactive
decay of a parent jth nuclide.
The last two terms are the losses due to neutron absorption and
radioactive decay of the ith isotope respectively.
· The composition of the reactor fuel after a given time t can be obtained
by solving the set of simultaneous equations (5.58) and inserting as a
boundary condition the initial fuel composition at t = 0. This change in
fuel composition will affect both thermal fission factor D , and the
thermal utilization factor, f , and hence D f . This will be given by:
C5 × s × f 5 × n5 + C9 × s × f 9 × n9 + C1 × s × f1 × n1
h f = ........(5.59)
f m
C5s 5 + C8s 8 + C9s 9 + C0s 0 + C1s1 + ( Cms m )
f f
Here the subscript 5 and 8 refer to U-235 and U-238, and subscript 9, 0
and 1 refer to Pu-239, Pu-240 and Pu-241. The absorption in materials
with concentration Cm in the fuel other than fuel isotopes is included in
the last term of the denominator.
· In general, it has been found that the changes in the composition of the
fuel will not change the value of the fast fission factor A , or the
resonance capture probability p, since both of these are dependent on
the U-238 concentration in the fuel, and this will not change
significantly with time.
128
lifetimes are increased by a factor of 85. However, when fuel are burnt,
Pu-239 are formed. Plutonium-239 has lower delayed neutron yields,
and the effect is to reduce the mean lifetime and so make reactor
behavior more dynamic. For example, for equilibrium fuel with
roughly equal U-235 and U-239 content, the neutron generation
lifetime will be around 0.06 seconds. A + 1mk insertion will therefore
cause a more rapid power increase for equilibrium fuel than for fresh
fuel. Therefore, in the design of the reactor control system, an
allowance must be made for the decrease in reactor period as the
plutonium concentration increases.
· The burnup of the fuel will produce change in the neutron flux
distribution across the reactor. This occurs because the burnup at any
point in the core is proportional to the flux at that point. Since the flux
is higher in the central region of the reactor, the burnup proceeds at a
greater rate there. The gradual reduction of reactivity in these regions
will cause the flux to be depressed there relative to what it was before.
This changing flux shape across the core affects refuelling scheme to
be used if maximum burnup is to be obtained from the fuel.
Because Tellurium decays very fast into idoine I-135, we may assume
that iodine is produced directly from fission. Hence the Iodine and xenon
concentrations are given by:
dX
= g X S f f + l I I - l X X - s X f X ...................(5.60)
dt
dI
= g I S f f - l I I ...............................................(5.61)
dt
129
B = neutron flux in neutrons/cm2.sec
The terms on the right hand side of the ODEs can be described as follows:
g IS ff 0
I0 = ........................(5.62)
l I
(l I I0 + g X S f f 0 )
X0 = ...............(5.63)
(l X + s Xf 0 )
(g I + g X ) S f f 0
X0 = ................(5.64)
lX + s X f 0
(g I + g X ) S f f o
X0 » therefore proportional to flux ..........(5.65)
lX
(g I + g X ) S f
X0 » , therefore constant...........(5.66)
sX
130
So for high power reactor, the steady state xenon concentration is
independent of power.
Sa Sa
Let f = without xenon, and f ¢ = with xenon
Sa + Sm Sa + Sm + S X
Note that S m is the absorption cross section of all materials in the reactor
that is not fuel or xenon.
Now
K e¢ - K e f¢- f f S + Sm + S X SX
= = 1- =1- a =-
K e¢ f¢ f¢ Sa + Sm Sa + Sm
.....(5.67)
SX Sa
=- × = -R × f
Sa Sa + Sm
SX X ×I
At equilibrium, = 0 X ..................................(5.68)
Sa Sa
S X s X × (g I + g X ) S f f 0
R= = ...............(5.69)
Sa ( lX + s X f 0 ) × S a
When the reactor is critical with its xenon load (i.e. K e¢ = 1), and - Rf
represents the reactivity loss due to the xenon. In other words, the xenon
load is the product of R and f.
5f
; I X = 3.5 x 10-18 cm2 ; C X =0.003; C I =0.056; l X =2.1x10-5; =0.54
5a
(for natural uranium), substituting these numbers into (5.69), we have R =
0.028. Given f = 0.9, the xenon load R will be 0.0252. This means that this
reactor would suffer a reactivity loss of 25.2 mk caused by equilibrium
xenon poisoning.
131
The presence of an equilibrium xenon load means that enough excess
reactivity must be built into the reactor so that the regulating system can
compensate for the loss in reactivity as the xenon builds up.
Let at t=0 the reactor be shutdown suddenly from a steady state flux B
with corresponding xenon and iodine concentrations X 0, I 0 . Putting B = 0
in equation (5.60), (5.61), we have:
dX
= - lX X + l I I
dt
.......................(5.70)
dI
= -l I I
dt
I ( t ) = I 0 e - lI t
lI ...........(5.71)
X (t ) = X 0 e -l X t + I 0 ( e -lI t - e -l X t )
l X -l I
It can be seen that for initial steady flux of B = 1014 n.cm-2.sec-1, the
xenon reactivity load rises from an equilibrium value of - 28 mk at t = 0 to
a maximum value of -132 mk somewhere between 7 and 12 hours after
shutdown. Note that the peak xenon load level depends proportionally on
the reactor flux level before shutdown. Physically, after shutdown, there
are no neutrons to convert the xenon, so it disappears only by slow natural
decay. But xenon is being produced at a faster rate by the decay of iodine,
so its level increases. Due to the decay of iodine, a maximum is reached,
and the xenon decays thereafter.
To find the peak xenon load, we differentiate (5.71) with time and put
dX/dt = 0, we have
dX l I I0
= ( l X e -l X t - l I e -lI t ) - l X X 0 e -l X t = 0 .......(5.72)
dt l X -l I
1 l I l X ( I0 + X 0 ) - l X
2
X0
t= ln( ) ...........(5.73)
l X -l I l I 2 I0
X0 C
X max = ( I 0 + AX 0 )( B + AB ) .......................(5.74)
I0
132
where
lX lX lX
A = 1- ;B = ;C = ............(5.75)
lI lI l I -l X
The consequences for peaking xenon load are far reaching - it is not
economical to provide the extra fuel required to override such large
negative reactivities. The net result is that after a shutdown, the reactor
can only be restarted within about 40 minutes, before the xenon has grown
too much. If this is not possible, it is necessary to wait about 40 hours for
the xenon to decay naturally to a sufficiently low level.
Xenon effects also accompany both power reductions and increases. For
a power reduction, the effect is essentially as above, with the maximum
xenon reactivity load depending on the magnitude and the rate of the
power reduction. With a power increase, the effect is reversed, and a strong
dip occurs in the xenon concentration. Its depth again depends on the size
and the rate of the increase. To prevent a sharp increase of reactor power
due to this dip, it is necessary to introduce neutron absorbing reactivity
control devices to compensate for it.
With temperature changes, the nuclear cross sections will change due to
changes in mean energies, and density variations will affect the neutron
mean free path and the leakage probabilities; as well, there is the
broadening of the resonance cross section due to an increase in the
temperature in the fuel, known as the “Doppler effect”.
133
where = T is the temperature coefficient ( -15.06 ´ 10-6 /deg. C for fresh
fuel)
For example, in order to bring the reactor power from cold 25 deg. C to hot
condition at 271 deg, C, the reactivity needed to compensate for reactivity
loss due to temperature effects is
It has been shown that reactor core voiding (due to boiling or in extreme
case, loss of coolant accident) will lead to positive reactivity excursion. Is
generally expressed as:
where ,KV max =the maximum reactivity from coolant voiding ( typically 2
mk)
? = Coolant quality
134
dCI
= 9. 445 E - 3 * N - 2. 8717 E - 5 * CI
dt
dCXE
= 9.167 E - 4 * N + 2. 8717 E - 5 * CI - ( 2.1E - 5 + 3. 5 E - 4 * N ) * CXE
dt
CI = concentration of iodine(mK)
135
9) Create a block for calculating the xenon/iodine concentration transient
using the custom algorithm ALG801. Remember to initialize the xenon
and iodine concentration by depositing their values in the output
section.
1) TIMEKEEPER/IT1
14) Use the above model and create transients for the cases of a step
change of:
15) Use a speedup factor of 600. Why is it better to use a smaller speedup
factor, i.e. a smaller integration time-step?
16) Using the model and with the help of the LABVIEW executable,
CURVES.EXE, can you find the following:
136
a) If a reactor core is tripped from full power, what will be the poison-
out period? When must the reactor be re-started before a 'poison-
out' occurs? And if a 'poison-out' did occur, how long before the
reactor can be re-started?
137
6.0 Reactor Control Simulation Model Development
(Courtesy of ABB)
BWR is the abbreviation for the Boiling Water Reactor. These reactors
were originally designed by Allis-Chambers and General Electric (GE).
The General Electric design has survived, whereas all Allis-Chambers
units are now shutdown. The first GE US commercial plant was at
Humboldt Bay (near Eureka) in California. Other suppliers of the BWR
design world-wide include - ASEA-Atom, Kraftwerk Union, Hitachi.
Commercial BWR reactors may be found in Finland, Germany, India,
Japan, Mexico, Netherlands, Spain, Sweden, Switzerland, and Taiwan.
138
(Courtesy of ABB)
In the figure above, water is circulated through the Reactor Core picking
up heat as the water moves past the fuel assemblies. The water eventually
is heated enough to convert to steam. Steam separators in the upper part of
the reactor remove water from the steam.
The steam then passes through the Main Steam Lines to the Turbine-
Generators. The steam typically goes first to a smaller High Pressure (HP)
Turbine, then passes to Moisture Separators (not shown), then to the 2 or 3
larger Low Pressure (LP) Turbines. In the sketch above there are 3 low
pressure turbines, as is common for 1000 MW(e) plant. The turbines are
connected to each other and to the Generator by a long shaft (not one
piece).
The steam, after passing through the turbines, then condenses in the
Condenser, which is at a vacuum and is cooled by ocean, sea, lake, or river
water. The condensed steam then is pumped to Low Pressure Feedwater
Heaters (shown but not identified). The water then passes to the Feedwater
Pumps which in turn, pump the water to the reactor and start the cycle all
over again.
The BWR is unique in that the Control Rods, used to shutdown the
reactor and maintain an uniform power distribution across the reactor,
139
are inserted from the bottom by a high pressure hydraulically operated
system. The BWR also has a Torus or a Suppression Pool. The torus or
suppression pool is used to remove heat released if an event occurs in
which large quantities of steam are released from the reactor or the Reactor
Recirculation System, used to circulate water through the reactor.
The BWR reactor typically allows bulk boiling of the water in the
reactor. The operating temperature of the reactor is approximately 570F
producing steam at a pressure of about 1000 pounds per square inch.
A BWR plant has numerous control rods installed inside the core. Each
control rod is provided with a drive system and hydraulic control system.
Rod position is adjusted by drawing out or insertion by means of the
hydraulic pressure that is remotely controlled in the control room. Only
one rod is operative at a time.
When the power output is changed by using a control rod, the changing
rate is approximately 2 % power per minute. In the event of emergency
shutdown of the reactor, all control rods are promptly inserted at once
(scram) by the reactor trip system.
140
6.2 Brief Review of PWR Reactor Power Control
141
Feedwater System and passes on the outside of those steam generator
tubes, is heated and converted to steam. The steam then passes through the
a Main Steam Line to the Turbine, which is connected to and turns the
Generator.
Light water PWR plants have four main control systems: control rod
control (reactor power control), steam generator level control, pressurizer
pressure control and pressurizer level control.
The control rod control is used for reactor power control (Figure 10).
The control rods are moved up and down when a deviation between
“primary power”, PPRI, and the “reference power”, PREF, formed by the
turbine first-stage pressure (secondary power) exceeds the predetermined
setpoint.
Figure 10 - Control Block Diagram for Reactor Control for a typical PWR.
(Courtesy of IAEA)
142
The boron concentration control system (which requires manual action)
is used for the relatively long-term and slow core reactivity control,
whereas the reactivity control (including reactor trip) responsive to load
changes during normal operation is performed by controlling the rod
positions
143
graphite as its moderator. It is very different from most other power reactor
designs as it was intended and used for both plutonium and power
production. The combination of graphite moderator and water coolant is
found in no other power reactors. The design characteristics of the reactor
were shown, in the Chernobyl accident, to cause instability when low
power. This was due primarily to control rod design and a positive void
coefficient. A number of significant design changes have now been made to
address these problems.
144
2. Pressure tubes -within the reactor each fuel assembly is positioned in
its own pressure tube or channel. Each channel is individually cooled
by pressurized water.
3. Graphite moderator - a series of graphite blocks surround, and hence
separate, the pressure tubes. They act as a moderator to slow down the
neutrons released during fission. This is necessary for continuous
fission to be maintained. Conductance of heat between the blocks is
enhanced by a mixture of helium and nitrogen gas.
4. Control rods - boron carbide control rods absorb neutrons to control
the rate of fission. A few short rods, inserted upwards from the bottom
of the core, even the distribution of power across the reactor. The main
control rods are inserted from the top down and provide automatic,
manual or emergency control. The automatic rods are regulated by
feedback from in-core detectors. If there is a deviation from normal
operating parameters (e.g. increased reactor power level), the rods
can be dropped into the core to reduce or stop reactor activity. A
number of rods normally remain in the core during operation.
5. Coolant - two separate water coolant systems each with four pumps
circulate water through the pressure tubes. Ninety-five percent of the
heat from fission is transferred to the coolant. There is also an
emergency core cooling system which will come into operation if
either coolant circuit is interrupted.
6. Steam separator - in the RBMK design, boiling occurs. The steam
produced passes to the Steam Separator which separates water from the
steam. The steam then passes to the Turbine as in the Boiling Water
Reactor design. Similar to the BWR case, the steam is radioactive,
however, the steam separator introduces a delay time so radiation
levels near the turbine may not be as high as in the BWR case.
7. Containment - the reactor core is located in a concrete lined cavity that
acts as a radiation shield. The upper shield or pile cap above the core,
is made of steel and supports the fuel assemblies. The steam separators
of the coolant systems are housed in their own concrete shields.
In the RBMK design, water coolant slows the reaction down, meaning
that the water coolant absorbs the neutrons, a necessary part of the
reaction. In other words, the water coolant additionally serves to control
the reaction. The lack of water increases the reaction in the RBMK design,
due to positive void reactivity feedback. This is one chief difference
between the power plants in the United States, known as the LWR (Light
Water Reactors), and the RBMK type. In the LWR the lack of water kills
the reaction, opposite of a RBMK.
145
6.4 Pressurized Heavy Water Reactor
(Courtesy of AECL)
146
Each pressure tube is isolated and insulated from the heavy water
moderator by a concentric calandria tube and a gas annulus. Consequently
the moderator system is operated at low temperature and pressure. The
reactivity control and shutdown mechanism reside in the low-pressure
moderator, thus simplifying their design, construction and maintenance
and eliminating the possibility of their ejection in an accident situation. As
well, this cool moderator can act as a heat sink under certain accident
conditions.
The use of natural uranium fuel in an optimized lattice, and heavy water
as moderator and coolant, combined with the capability to refuel the
reactor while at full power, gives the CANDU reactor its good economy
and low excess reactivity. This results in a power with very low fuel costs.
The only reactivity feedback effects that are of any consequence are the
coolant density (or void), fuel temperature and xenon effects. Since the
moderator temperature is controlled independently of the primary coolant,
it does not contribute any significant reactivity feedback. The fuel
temperature reactivity coefficient is negative and therefore stabilizing
whereas the coolant void reactivity coefficient is positive and therefore
destabilizing. The sum of these, and lesser reactivity feedback effects, is a
power reactivity coefficient that is near zero under normal operating
conditions.
(B) Surge Tank Level - it is controlled by the feed and bleed valves shown
on the circuit.
In the boiler, the energy is transferred to ordinary light water steam, and
the steam drives a conventional turbine-generator set. The feedwater
valves, pumps and preheaters complete the chain from the condenser back
to the boiler. Steam can also be discharged directly to the condenser via the
discharge valves with a capacity of 68 % of full power steam flow, or to
the atmosphere via atmospheric steam discharge valves with a capacity of
147
10 %. The need for these valves is related to xenon poison effects
explained earlier.
The boiler level control system (not indicated on the diagram) control
the light water level in the boiler drum via feedwater valves.
ALL SYSTEMS
COMMON
MODERATOR FEEDWATER
SERVICES
CONTAINMENT
A simplified diagram of a generic CANDU nuclear power plant and its digital
Figure 11-
computer control system (Courtesy of AECL)
148
The three main control programs are:
149
3. 2 banks of mechanical Absorber rods (total # of rods = 4) - these
Absorber rods decrease the reactivity inside the reactor by inserting
themselves into the core. Being neutron absorbers, insertion into
the core means more neutron absorption, thereby reducing the
reactivity and the neutron power of the reactor. The total reactivity
worth of both fully inserted absorber rods is approximately - 9mk.
Each absorber bank is assumed to have the same reactivity worth.
150
Figure 13 - Reactivity Mechanism Desk - showing the arrangement of
Adjuster rods, Control Absorbers, Shutdown rods, and Zone Controllers.
(Courtesy of AECL)
151
6.4.2 Moderator Liquid Poison Addition and Removal System
Poison is removed by the ion exchange purification system. But the time
constant of the removal is of the order of hours. Operation of poison
addition and removal is controlled manually.
152
6.5 PHW Reactor Regulating System Brief Overview
The overview of the PHW Reactor Regulating System is shown in the
following block diagram (figure 14):
(Courtesy of AECL)
153
· The Power Measurement & Calibration Program determines average
reactor power and its distribution through the reactor.
· The Demand Power Routine sets the desired reactor power which
depends on the mode of reactor operation. In the “Normal” mode of
operation, reactor power is set to follow turbine generator power such
that boiler pressure is kept constant under all conditions. In this mode,
the setpoint demand signal comes from the Boiler Pressure Controller
(BPC). In the upset conditions, the mode of plant operation is switched
to “Alternate” Mode, in which case the reactor setpoint is manually
adjusted or adjusted via a Reactor Power Setback Routine.
154
4. The adjusters and mechanical control absorbers are driven at a speed
proportional to power error to minimize, at low power errors, the shim
reactivity rate which must be canceled by zone controllers.
(Courtesy of AECL)
155
6.6 Generic CANDU Simulator - RRS Familiarization
Summary of Simulator Features.
156
6.6.1 Familiarization with CANDU Generic Simulator
(1) First load the CANDU generic Simulator up as per instructions given
in the class.
· the window displaying ‘CASSIM’ will be green and the counter under
it will not be incrementing when the simulator is frozen (i.e. the model
programs are not executing), and will turn red and the counter will
increment when the simulator is running;
· to stop (freeze) LabVIEW click once on the ‘STOP’ sign at the top left
hand corner; to restart ‘LabVIEW’ click on the Þ symbol at the top
left hand corner;
· to start the simulation click on ‘Run’ at the bottom right hand corner;
to ‘Stop’ the simulation click on ‘Freeze’ at the bottom right hand
corner;
· the bottom of the screen shows the values of the following major plant
parameters:
· the bottom left hand corner allows the initiation of two major plant
events:
1. ‘Reactor Trip’
157
2. ‘Turbine Trip’
· the box above the Trip buttons shows the display currently selected
(i.e. ‘Plant Overview’); by clicking and holding on the arrow in this
box the titles of the other displays will be shown, and a new one can be
selected by highlighting it;
· the remaining buttons in the bottom right hand corner allow control of
the simulation one iteration at a time (‘Iterate’); the selection of
initialization points (‘IC’); insertion of malfunctions (‘Malf’); and
calling up the ‘Help’ screen.
· The screen shows the status of Shutdown System (SDS) #1, as well as
the reactivity contributions of each device and physical phenomenon
that is relevant to reactor operations.
· The positions of each of the two SDS1 SHUTDOWN ROD banks are
shown relative to their normal (fully withdrawn) position.
· Note also that when the reactor is critical the Total Reactivity must be
zero.
· This screen shows the Limit Control Diagram, and the status of the
three reactivity control devices that are under the control of the Reactor
Regulating System (RRS).
158
· If power error is negative, more (positive) reactivity is needed, hence
liquid zone level will decrease and if this is insufficient, absorber rods
and adjuster rods will be withdrawn from the reactor.
· The ABSORBERS are moved in two banks, and are normally outside
the core. They can be moved by RRS if AUTO is selected, or manually
if they are placed in MANUAL mode.
· The ADJUSTERS are moved in eight banks, and are normally fully
inserted into the core. They can be moved by RRS if AUTO is selected,
or manually if they are placed in MANUAL mode.
· This screen permits control of reactor power setpoint and its rate of
change while under Reactor Regulating System (RRS) control, i.e. in
‘alternate’ mode. Several of the parameters key to RRS operation are
displayed on this page.
159
· TRIP status is indicated by YES or NO; trip is initiated by the
Shutdown Systems, if the condition clears, it can be reset from here.
Note however, that the tripped SDS#1 must also be reset before RRS
will pull out the shutdown rods, this must be done on the Shutdown
Rods Page
· Key components of RRS and DPR control algorithm are also shown on
this screen.
· HOLD POWER 'On' will select ‘alternate mode’ and stops any
requested changes in DEMANDED POWER SETPOINT.
· It is important to note that Reactor Power Setpoint Target and Rate are
specified on the Simulator in terms of %FP and %FP/sec, i.e. as linear
measurements, instead of the logarithmic values used in practice. The
requested rate of change should be no greater than 7% of actual power
per second in order to avoid a reactor LOG RATE trip. This is readily
achieved in the 'at-power' range (above 15%FP), but only very small
rates should be used at low reactor power levels (below 1%FP), such as
encountered after a reactor trip.
· This screen permits control of station load setpoint and its rate of
change while under Unit Power Regulator (UPR) control, i.e. ‘normal’
mode. Control of the Main Steam Header Pressure is also through this
screen, but this is not usually changed under normal operating
conditions.
160
· MAIN STEAM HEADER PRESSURE SETPOINT (MPa) - alters the
setpoint of the Steam Generator Pressure Controller, which is rarely
done during power operation. Caution must be exercised when using
this feature on the Simulator, since the requested change takes place in
a step fashion as soon as the change is made; changes should be made
in increments of 0.1 MPa.
161
· return reactor power to 100%FP at 0.5%/sec, freeze as soon as
Reactor Neutron Power reaches 100% and record parameter
values
· unfreeze and let parameters stabilize, record parameter values
· go to the Reactivity Control page, unfreeze, and note the path of the
operating point on the attached diagram, until power error has
stabilized at or near zero (about 3 - 4 minutes)
· confirm the value of average zone level on the Plant Overview page
100
90
80
70
AVE
60
ZONE
LEVEL 50
(%)
40
30
20
10
0
-7.0 -6.0 -5.0 -4.0 -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0
POWER ERROR (%FP) = ACTUAL - DEMANDED
(error amplitude and rate)
Question - Why is the final zone level higher than the original zone
level?
162
(3) REACTOR AND RRS RESPONSE TO POWER MANEUVER
· go to the Reactivity Control page, unfreeze, and note the path of the
operating point on the attached diagram, until at least one Adjuster
Rod bank is out of the reactor (about 20 minutes) - once the first
Adjuster Bank is more than 50% withdrawn, place Absorbers on
Manual and drive them fully OUT.
· note the time taken from the start of lowering reactor power until
steady operation within the specified error limits is achieved
163
· on the Reactivity Control page place the Absorbers and Adjusters
onto Manual
· manually drive both banks of Absorbers into the core, and by the
simultaneous withdraw of one or more the Adjuster banks keep
reactor neutron power within 2% of 100%FP
· once all the Absorbers are fully in the core, attempt to drive all
remaining Adjusters simultaneously out of the core
100
90
80
70
AVE
60
ZONE
LEVEL 50
(%)
40
30
20
10
0
-7.0 -6.0 -5.0 -4.0 -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0
POWER ERROR (%FP) = ACTUAL - DEMANDED
(error amplitude and rate)
164
· note OUC mode and reactor power
· once the Absorbers have fully withdrawn from the reactor, raise
reactor power to a level dependent on the number of Absorber
banks out of the reactor: for each bank partially or fully out, reactor
power is limited by 5% (i.e. one bank - 95%FP, two banks -
90%FP, etc.)
· wait until Generator power is zero and reactor neuron power less
than 0.1%
· record the time (using the display under the chart recorders) needed
to withdraw all shutdown rods
· wait one hour before resetting Reactor Trip and SDS#1 (use
initialization point given by instructor).
· after the shutdown rods have been withdrawn observe the status of
the reactivity control devices
· note the reactivity changes that have taken place, in particular note
the magnitude and estimate the rate of change of Xenon reactivity
buildup
165
(10) RESPONSE TO SDS#1 REACTOR TRIP AND POISON OVERRIDE
· using the data from the previous two exercises, estimate the time
available to the operator from the initiation of the reactor trip until
the trip must be reset to avoid a poison outage: the desired end state
of this exercise is reactor power at 60%FP and less than one bank
of adjuster rods left in the core (i.e. the last bank is partially
withdrawn)
· note the final state of the adjuster rods and Average Liquid Zone
level
166
6.7 Modeling Exercise - Reactivity Control Mechanism: Liquid Zone
Problem Description:
Create a process model as presented in the figure below.
Level controller
167
1. Identify the process equipment which will be required to be modeled,
and find the corresponding generic algorithms to match the equipment.
2. Identify which method of flow calculation will be used for each flow
paths in the system.
EXTERNAL NODE INTERNAL NODE INTERNAL NODE INTERNAL NODE EXTERNAL NODE
Equipment locations:
NHAX01_N01: ELEV HEAD OF 5 metres
NHAN01_N02: PUMP P1
NHAN02_N03: VALVE V1
6. The flow between TANK1 and TANK2 will be calculated using the
network method, and the return flow from TANK2 to TANK1 will be done
outside of the network.
168
Equipment/function Block name Generic Connection
algorithm used comments
Pump 1 Breaker AAA_P1_S ALG133, The inputs to
BREAKER control the breaker
will be constant
tags.
Pump 1 Motor AAA_P1_MTR ALG134, MOTOR The inputs to
control the motor
will be from the
breaker status from
the breaker block.
Pump 1 Hydraulics AAA_P1 ALG102, The main inputs
PUMP_SIM will be the motor
speed from the
motor block, and
the flow from the
network block
Valve 1 AAA_V1 ALG108, CV_SIM The inputs to the
valve position will
be from the
controller block.
Valve 2 AAA_V2 ALG106, The inputs to
VLV_ONOFF_SI control the valve
M position will be a
constant tag.
Conductance of AAA_V2_COND ALG301, Valve 2
TANK2 return flow 2_MULTI normalized port-
area is multiplied
together will
TANK2 level
factor, such that if
TANK2 runs dry,
it will stop the
return flow.
TANK 2 return flow AAA_RETURN ALG140, The input will the
calculation CONDUCT_FI TANK2 return
flow conductance
TANK1 total inlet AAA_RET_TOT ALG_311, 2_SUM TANK1 total inlet
flow flow is sum of the
return flow from
TANK2 and the
overflow from
TANK2
TANK1 AAA_TANK1 ALG128, TANK The inlet flow is
calculated from
TANK1 total inlet
flow block and the
outflow from the
network
TANK2 AAA_TANK2 ALG128, TANK The inlet flow is
calculated from the
network, and the
outlet flow from the
TANK2 return flow
calculation
169
8. Create a NMB block(ALG452) for communication with LABVIEW.
9. Use CASSBASE and create the blocks for modeling the process
equipment, and non-network type flow calculations.
· Input connections can be done after all the blocks are created.
10. Use CASSBASE and create the blocks related to the hydraulic
network. Remember to create marker blocks for enclosing the network
related blocks.
13. Use CASSBASE to EDIT the blocks and make the connections to the
block inputs.
14. Make sure the connections between the blocks are made correctly, use
the VERIFY command to locate any hanging tags, and clear these
hanging tags by reconnecting the correct output tags to the inputs in
question.
15. Connect the variables you would like to trend to block NMB (up to 8 is
available).
· Observe the transient when you try to set the level setpoint of
TANK2 to 1.5m.
170
16. Operate set the level setpoint back to 2m.
liquid REACTOR
Delay Tank Liquid zone
Liquid zone
zone inlet reactor flux
pump control
valve
PID
demanded setpoint
reactor power
controller
Rate limit
to 0.8%
actual setpoint
Create a model to control the reactor power using the liquid zone.
171
· the liquid zone hydraulic model, IAEA3.
5. Connect the reactivity change due to liquid zone level to the reactor
block.
9. Select the parameters you wish to trend (e.g. reactor flux, TANK2
level.) and connect it to the NMB block. Use the LABVIEW
executable TRENDS.EXE to trend the variables.
· If you are successful, can you raise the reactor power back to
100%? If not why?
11. Use CASSBASE to extract the whole model out, and call it IAEA5.
172
7.0 Development of Adjuster, Absorber, Shutdown Rods Model
where
systems (K)
where
173
· The reactivity change worth of the adjuster banks and absorbers are
assumed to be:
DKCA = åi =ABSR
2
iWORTH * ( MAX ( -1, MIN ( 0 ,
1
Where
174
ABSRiWORTH = total reactivity change worth of absorber bank
i(K)(i=1,2)
· The total reactivity worth for the 2 absorber banks in core is -9 mk; the
total reactivity worth for all 8 banks of adjuster out-of-core is + 17 mk.
· At full speed, an adjuster rod takes 60 seconds for full travel; at full
speed, an absorber rod takes 150 seconds for full travel.
inserted(K)
· The total reactivity worth for the 2 shutdown rods is - 68 mk, when
fully inserted in the core. The maximum time for full rod drop should
be less than 2 sec.
· Design a reactor tripping logic. Make this trip is a latch in once set.
175
Can you model this using generic algorithm?
· Now test the model by manually tripping the reactor at full power.
Monitor xenon-iodine poison buildup. Try to speed up the xenon
poison buildup, if you can.
· Use the liquid zone to raise reactor power. If liquid zone cannot
recover full power, use adjuster rod to help. Do it slowly, and make
IC points on the way.
176
8.0 Reactor Startup and Approach to Criticality
1. The available reactivity is at its maximum, since the fuel has suffered
no burnup and fission product poisons do not exist yet.
4. The critical value of the variable that is used for the startup is
unknown, For instance, if the approach is used by gradual removal of
poison, then how much of poison removal is enough for criticality is
unknown.
177
S 0 L(1 - K em )
St = ...........................(8.1)
(1 - K e )
178
For a subcritical reactor, as Ke is increased, the power levels off more
slowly. As Ke is increased further step by step, a point will finally be
reached when the power will increase steadily without leveling off at all,
and this is the point where criticality (sustained nuclear fission chain
reaction) has been reached i.e. Ke =1. This is the basic method we use for
approaching criticality, by gradually increasing Ke. First let us review the
reactor startup instrumentation which makes this startup possible.
· The ion chamber system extends from 10-7 to 1.5 of full power.
Two sets of neutron detectors (BF3 counters) are used covering the range
from 10-14 of full power up to 10-6 of full power. One set, located in-core,
covers the range from spontaneous source level (10-14 full power) to 10-9
full power. The in-core detectors are installed via the vertical viewing port
penetration from the reactivity deck. This instrumentation is removed after
the low power commissioning period.
179
Figure 16 - Range of sensitivity for different reactor instrumentation.
180
The other set of neutron detectors is located out of core in the three spare
cavities of the shutdown system number 2 ion chamber housings.
181
Figure 17 - Arrangement of Startup Instrumentation.
182
8.2.2 Ion-Chamber Instrumentation
The ion chamber system is a part of the on-line power measurement
equipment of the reactor control system. This instrumentation is on-scale
during normal operation, and remains on-scale during a normal shutdown.
If the shutdown lasts more than 40 to 70 days (depending on reactor
operating history and ion chamber neutron and gamma sensitivity), the flux
decays below 10-7 of the full power value and startup neutron detectors are
needed.
Ion chambers measure only the average flux, as seen at the side of the
core, but are sensitive to very low flux levels, for monitoring sub-critical as
well as full power behavior. Because they are outside the calandria, the
neutron flux they detect has been attenuated by the moderator and any
poison dissolved in the moderator.
The control system flux detectors are of two types. One type has an
Inconel emitter and is used for the zone control system. The other type has
a vanadium emitter, and is used for the flux mapping system.
Both types of control system detectors, along with the Inconel detectors
used for the shutdown system number 1 high neutron power trip, are
contained in vertical flux detector assemblies located throughout the core.
These assemblies extend from the reactivity mechanisms deck to the
bottom of the reactor.
183
· Zone Control Flux Detectors (Inconel)
184
Reactors which do not have a natural photoneutron sources like CANDUs
are usually equipped with an artificial source, in the sense that some
Beryllium is included in the core which undergoes (g, n) reactions in much
the same way as deuterium.
1
= constant x ( 1 - Ke) ..........................(8.2)
count rate
As the moderator level is raised in small steps, the count rates will
increase, and the reciprocal count rates are plotted against moderator level.
A plot of inverse count rate versus moderator level will look like this:
Reciprocal
Count Rate
(c.p.s) Critical Moderator Level
1
= constant x ( 1 - Ke)
count rate
1
= A × DK
count rate Ke
= A × B × ( C - Ccritical ) .............(8.3)
185
C - Boron concentration at time t.
The procedures used are as follows: the count rate is measured for the
initial boron concentration, and the moderator purification circuit is then
valves in. Boron is removed for an interval of time ~ 20 minutes, and the
count rate were recorded again. This procedure is repeated a number of
times and the inverse count rate is plotted against total boron removal time.
Reciprocal
Count Rate
Critical Boron Concentration
(c.p.s)
at this moment
The theory is simple: if all the adjuster rods are in “Manual”, and the
liquid zone control is in “Auto”. As Ke gradually approaches 1 by slowly
removing Boron, liquid zone controllers always control water level at
setpoint, thus making zero reactivity contribution from liquid zone. As Ke
is very close to 1, and if Boron removal is stopped at that time, and one
makes a small step change in reactor power increase by RRS, the average
liquid zone level will drops (for power increase), if in fact Ke =1, after the
desired power level is reached, the liquid zone level will eventually return
to its original level setpoint after a few minutes. However, if in fact Ke is
less than 1, as we explained earlier, the neutron power after initial increase
will level off quickly. Because reactor power cannot reach setpoint, liquid
zone level cannot recover quickly to setpoint, and in fact it will continue to
drop in order to contribute reactivity for power increase. This is the
indication of a subcritical reactor, so Boron removal is resumed for some
time, and then the small step change in reactor power is repeated and the
liquid zone level response will be assessed again to determine criticality. In
the simulated reactor startup that we are going to practice, we will use the
liquid zone response technique to determine reactor criticality.
186
8.4 Practice Reactor Startup and Approach to Criticality Using
Simulator
1. First load the simulator up and initialize the simulator with the IC =
GSS_SD.ic (that is, Reactor in Guaranteed Shutdown State with
moderator overpoisoned).
2. In the Plant Overview screen, you will see the overall plant condition
from the annunciators: Turbine trip; Primary Heat Transport System is
depressurized at low pressure < 300 KPa, at temperature 54 deg. C;
PHT pumps shutdown.
(a) Special Liquid Zone Controls under the box titled “LIQUID ZONE
CONTROL MODE”. A control push-button is provided for
switching the Liquid Zone Control between “Standard” and
“Special Shutdown” Mode. Some explanations are necessary for
these two modes - under normal operating conditions, that is when
reactor power is > 0.1 %, the liquid zone controller is in “standard”
mode, that is, the zone level controller signal is proportional to the
combination of reactor power errors and flux tilt errors
components. However, when reactor power is less than 0.1 % and
at such low power, flux tilt control is not necessary and operator
will put the zone controller in “ Special Shutdown” Mode, meaning
that the zone level controller now responds to local level setpoint as
inputted by the operator. The zone level controller is now
proportional to the level errors only. As you will see later, as we
approach criticality, the liquid zone controller will be “kicked” out
of “Special Shutdown Mode” into “Standard” mode, when power
level is > 0.1 % (-3 decades) and there is a positive power error
present.
(b) The reactor startup instrumentations under the box titled “POWER
LEVEL READINGS”. The instrumentation simulated here are :
· as well, the Ion Chambers and Flux Detectors are simulated, but
they are all off-scale, meaning the reactor power is below -6
decades.
(c) The reactor power maneuvering rates - PWR LOG rate (DEC/sec);
and power error (DEC) are also indicated here for monitoring
187
purposes. The Xenon Load (mk) normally would not be displayed
in real situation, but is displayed here for our study and analysis.
(d) SDS1 Reactor trip Setpoint can be set here under titled box “SDS1
Trip Setpoint”. For safety concern, you always set trip setpoint just
slightly above the current the reactor power level, to prevent the
unlikely situation of reactor power running away. You can set the
setpoint in % full power or in decades. Currently the SDS1 setpoint
is set at -2 decades. Remember to raise the SDS1 trip setpoint as
the reactor power level reaches -2 decades. In real situation, there
will be a warning message informing a “pending” reactor NOP
(neutron overpower) trip, when the trip setpoint is approaching.
(g) The graphical trends are provided here for count rate (counts per
sec); power log rate (dec/s); ion-chamber detectors (dec); boron
poison load (mk); power error (dec); average zone level (%).
Note:
(1) For our purpose, we can ignore the first step regarding operations
in the PHT regarding the warming up of the pressurizer liquid, and
go to Step #2 right away.
188
(2) The currently simulated startup counters will saturate at 6000 cps at
the “Long” position. One should take the approach of plotting the
inverse counts rate versus total boron removal time in order to
arrive at the extrapolated “critical” boron concentration. But to save
time for our exercise, we are using the method of liquid zone
controller response as an approximate indication of “criticality”.
Operating Phases Unit Load Reactor RRS PHT SDC SGPC Check/Monitor
1. Initial State: De- Þ Liquid zones Þ HTS in a cold Þ Both side of Þ In ‘Hold’
pressurized control in de-pressurized the SDC is in mode with
cooldown state: Special S/D state: service under a SG
HTS temperature < mode with a ‘SDC pump’ pressure
54°C and pressure level setpoint Þ HTS mode. setpoint of
< 300 KPa. Reactor of 40%. temperature < 461 KPa
has been down for 4 54°C Þ For SDC #2,
weeks. Þ Adjusters the flow path
control in Þ HTS pressure is MV23,
manual < 300 KPa MV45, SDC
P2, MV31
Þ Absorber Þ HTS inventory (bypass MV18
control in control is in is closed),
Þ Reactor auto. ‘SOLID’ mode Bleed cooler,
power as both the MV20, MV28,
level Þ Ion chamber pressurizer and MV29 to RIH
depends will be off- the bleed #1 & #3.
on how scale if the condenser are
long the reactor power isolated Þ Some bleed to
reactor < -6.25 MV10 (10%
has been Decades Þ Bleed flow is open) &
‘down’ controlled by CV15(Bleed
for. The CV15, the condenser
last group bleed LCV) to
of photo- condenser Purification.
neutrons level control The tempering
modeled valves and the valve MV19
has a half pressure should be
life of the letdown valve closed.
order of MV10 which
days. is opened at
10%.
1. Continued: Initial Þ For SDC#1,
State: De- the flow path
pressurized is MV12,
cooldown state: MV46, SDC
HTS temperature < P1, MV2, SDC
54°C and pressure HX, MV8,
< 300 KPa. Reactor MV9 to RIH
has been down for 4 #2 & #4. The
weeks. bypass line
valve MV4 is
closed and the
Tempering line
MV5 is closed.
189
2. Reactor 2. Set RRS 1. Turn all Þ Circuit #1 & 2 Þ In ‘Hold’ Þ SDS1 poised
Approaching Startup flag pressurizer both in service. mode with with NOP set at
Criticality - Phase1 to ‘ON’ . heaters to a SG LOW POWER
(from screen ‘auto’ and set pressure TRIP= 1 %
Reactor the solid setpoint of
SU/SD mode 461 KPa Þ Moderator
Controls) pressurizer Poison addition
pressure available.
Þ Startup flag setpoint to 7.3
will ensure MPa Þ Note Power
that the RRS (from screen Level and Zone
use a PHT Pressure Level
minimum Control)
reactor power
value of -5.9 Þ The
decades in its pressurizer
calc. may take a
long time to
3. Remove warmup
Poison by IX depending on
columns its initial
until Ion condition.
Chambers (warmup time
are on scale can take up to
and flux 9 hours)
detector’s
response to Þ Speedup
poison control is
removal is available to
significant artificially
(power speedup
changes warmup
within 1-2 process
minutes).
(from screen
Reactor
SU/SD
Controls)
5. Place all
adjuster
rods on
Manual and
fully insert
them into
the core
(from screen
Reactivity
Control)
Þ LZCS
remains in
Special SD
with a level
setpoint of
40%
190
3. Reactor 6. Stop poison Þ Monitor time to
Approaching removal double the
Criticality- Phase 2 when Ion count, which
Chamber gives some
indicates indication of
reactor proximity to
power > -6 Criticality.
decades. Adjust range if
(from screen count saturates
Reactor
SU/SD Þ LZC will be
Controls) "kickout" of
Special S/D
7. Raise power Mode on
by inputting positive power
a reactor error.
power
setpoint
which is 0.3
decade
above the
present
power. (from
screen
Reactor
SU/SD
Controls)
8. Monitor the
average
zone level. If
the average
level
3. Continued: Reactor should drop
Approaching below 25 %,
Criticality- Phase 2 press
the ‘HOLD
POWER’
button.
(from screen
Reactor
SU/SD
Controls)
Þ Note: the
average zone
level is
available as
one of the 6
trends on the
screen
‘Reactor
SU/SD
Controls’
9. Resume
poison
removal so
that the
zones must
rise to hold
power.
Monitor the
zone level
rises and the
reactor
power. As
191
3. Continued: Reactor the zones
Approaching rises to 70
Criticality- Phase 2 %, stop
poison
removal.
Note reactor
power and
zone level.
(from screen
Reactor
SU/SD
Controls)
10. Repeat Steps
7, 8, 9 until
after a
reactor
power SP
increase, the
average
zone level
drops, and
returns to
within 5%
of the
original
level after a
few minutes.
(from screen
Reactor
SU/SD
Controls)
12. Resume
Poison
removal
until the
average
zone level is
around 60
%, and
power
holding.
(from screen
Reactor
SU/SD
Controls)
192
9.0 From Simple Point Kinetic Reactor Model to High Fidelity Model
For Zone i,
dN 14 14 6
li i = (1 )å Kij N j Ni å K ij å l mCmj ....... ..(9.1)
- b - +
dt j=1 j =1 m=1
i = 1, 2, 3......14
dCmj
= b m N j - l mCmj ..................(9.2)
dt
j = 1, 2,......14
m = 1, 2,..... 6
where:
193
bm = Delayed neutron fraction of the mth group
6
={(1 - b) Kii - 1} i + ii å l mCmi +
dNi N K
dt li li m =1
(1 - b) 14 1 14 6
å ij j l j =å
K
li j =1, j ¹i
N + K
1, j ¹i
ij å l m Cmj
m =1
i
i = 1, 2, 3......14 ...............................................(9.3)
(a)
{(1 - >) K - 1}
Ni
ii
li
is the rate of neutronic flux changes in zone i due to the its zone reactivity.
(b)
6
Kii
li
ål m Cmi
m =1
(c)
(1 - >) 14
å Kij N j
li j =1, j ¹i
is the rate of neutronic flux changes in zone i due to the coupling effects of
the neutronic fluxes in the other 13 zones.
194
(d)
1 14 6
å ij å l m Cmj
K
li j =1, j ¹i m=1
is the rate of neutronic flux changes in zone i due to coupling effects from
the concentration of delayed neutron groups in the other 13 zones.
dN i ( DKi - b)
6 14
= N i + å l mCmi + a i å Kij N j +
*
dt Li m =1 j =1, j ¹i
14 6
1
å Kij å l m Cmj
li j =1, j ¹i m=1 .....(9.4)
Where
li
Li =
Kii
(1 - b)
ai =
li
l
l*m = m
li
Equation (9.4) is almost identical to point kinetic model for reactor zone i,
with the exception of an extra zone coupling source terms:
1 14 6
å ij å l m Cmj
14
=i åK
j =1, j ¹i
ij Nj K
li j =1, j ¹i m=1
and
If we absorb the "zone coupling effects" into "the zone i reactivity term
Drii" and rewrite the point kinetic equations for zone i:
14
dNi
( Drii + å Dr ij - b) 6
Ni + å l m Cmi
j =1, j ¹i
= *
dt Li m=1 ............(9.5)
Where
Drij is the reactivity change for zone i due to coupling effects in zone j
195
Equation (9.4) and (9.5) will be identical if
Dr × N1 1 6
( K12 å l mCm )
12
For zone 1 & 2, = a 1 K12 N 2 +
L1 l1 ZONE 2
m=1
Dr × N1 1 6
= a 1 K13 N 3 + ( K13 å l m Cm )
13
For zone 1 & 3,
L1 l1 ZONE 3
m=1
Dr × N1 1 6
= a 1 K14 N 4 + ( K14 å l m Cm )
14
For zone 1 & 4,
L1 l1 ZONE 4
m=1
Dr × N1 1 6
( K1,14 å l m Cm )
1,14
For zone 1 & 14, = a 1 K1,14 N14 +
L1 l1 ZONE14
m=1
Therefore the equation for Drij, the reactivity change for zone i due to
coupling effects in zone j is:
æ ö6
ç ål
Cm ÷ m
ç N j 1 ZONEj
m=1
÷
Dr ij = L i K ij ç a i + ÷ ...........(9.6)
ç N i li Ni ÷
ç ÷
è ø
òi å
n f
ò f (r ) n å f jth( r ) dr
*
th (r ) dr f
i
=
f f
k ij *
òj nå f òi f (r ) nå f th(r ) dr
*
th (r ) dr f
f f
....................(9.7)
where:
å f
= fission neutron production cross section
196
Bth(r) = real thermal flux at position r
Kij involves the computation of the distribution of the real thermal flux
and adjoint fast fluxes which usually involves a lot of CPU time. It would
be impractical to compute these fluxes in real-time environment. To deal
with this problem, an approximate method (3) is implemented to compute
the coupling coefficients in real time using a perturbation form (relative to
the equilibrium values):
K ij = K ij0 + g ij ..................(9.8)
in node j
14
gij = åP K
k =1
k ijk ................(9.9)
· K ij0 can be computed off-line for a given (nominal) core condition and
configuration by using equation (9.7), for all i, j.
197
calculate the new coupling coefficients Kij off-line, for all i, j. With the
use of equation (9.10), K ij1 can be obtained for all i, j. Repeat the same
for other zone to obtain K ijk , for k= 2, ...14.
in turn inputs into the reactivity change “input” of the affected reactor
zone.
· Each point kinetic model will calculate the neutron power based on 6
different neutron delay groups and 9 photoneutron groups, and the
overall change in reactivity.
· The reactor fuel elements are divided into 4 lumped channels. The total
heat generated from the fuel as calculated using the average reactor
flux is assumed to be divided equally among the 4 lumped channels.
· Each channel is assumed to have its own coolant flow, and its own
lumped fuel element. The temperature of the lumped fuel element in
each channel is calculated, and the lumped fuel sheath temperature will
be used in the coolant heat transfer calculation.
198
2. 2 independent lumped absorber banks, there is an equal reactivity
split between 14 zones.
3. liquid zones; the individual zone level will affect the respective
reactivity change in its designated zone reactor.
RRS
Reactivity Change due to
adjusters, absorbers,
Power Error shutdown rods and fuel
burnup
Average Zone Level
Average Reactor Flux
Zone 1 CV Position
Zone Controller 1 Reactor Zone 1I
Zone 2 CV Position
Zone Controller 2 Reactor Zone 2I
Average Reactor
Flux Calculation
To Decay Heat,
Zone 3 to 13 Zone 3 to 13 RRS
Average Zone
Level Calculation
To RRS
Flux Mapping
To Display
I reactivity changes due to temperature change, xenon poisoning and voiding are
within each reactor zone
z coupling is modelled between each neighbouring zones according to prescribed formula
199
Fig.19 Reactor Model Zones Configuration
10
13 8
11
14 9
12
3
1 6
4
2 7
5
(A) For each Zone Reactor, the point kinetic equations and the governing
equations for various reactivity feedbacks are:
dCi bi * NFLUX
= - liCi
dt TNEUTRON i = 1,....15 ..............(9.3-2)
200
3. The rate of change of neutron flux in a point kinetic model can be
expressed as:
NFLUX' + S i =1 Ai
15
NFLUX =
D K - b1
1 - Dt * ( + S i15=1 Bi )
TNEUTRON
where:
li * Ci * Dt
Ai =
1 + li * D t
li * bi * Dt
Bi =
(1 + li * Dt ) * TNEUTRON .................(9.3-4)
.........................................(9.3-5)
201
DKV = overall neutron reactivity change due to channel voiding (K)
systems (K)
zones (K)
202
QR = heat transfer to coolant(norm)
9. Assuming there are I coolant channels in the reactor zone, the reactivity
change due to channel voiding is assumed to be:
( S i =1>Li )
I
DKV = DKVO * N *
I
where
>Li = XEi - XEFP ...............(9.3-9)
dCXE
= 9.167 E - 4 * NFLUX + 2. 8717 E - 5 * CI -
dt
( 2.1E - 5 + 3. 5 E - 4 * NFLUX ) * CXE .......(9.3-11)
CI = iodine concentration
203
12. The reactivity change due to Zone Coupling effects:
æ 6
ö
ç ål Cm ÷
m
ç N j 1 ZONEj
m=1
÷
Dr ij = L i K ij ç a i + ÷ ............(9.3-12)
ç N i li Ni ÷
ç ÷
è ø
Where
14
K ij = K + å Pk × K ijk ............................(9.3-13)
0
ij
k =1
li
Li =
Kii
(1 - b)
ai =
li
l
l*m = m
li
· The reactivity change due to fuel burnup are calculated for one whole
reactor, and split up among the 14 zones according to certain criteria -
e.g. zones near the center of the reactor would have more contribution
than the outer zones.
204
· The total neutron power from the 14 zones reactor (each normalized)
will be summed up and then divided by 14 to get the normalized
overall reactor power. As well, each zonal power will provide into to
Flux Mapping routine for display.
· All the zone levels will be summed up and divided by 14 to get the
average zone level to RRS.
However, under reactor trip, stepback, or if the shutdown rods are fully
inserted, the power error is always set to 0.05.
1. If the absorber banks control is set to auto mode, then the absorber
banks will move according to the power error versus zone level
plot below.
205
ABSORBER BANKS AUTO MODE CONTROL
100
LIQUID ZONE LEVEL(%)
DRIVE 1
BANK IN
80
75
DRIVE 2
BANKS IN
DRIVE 1
DRIVE 2 BANK OUT
BANKS OUT
50
-6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6
POWER ERROR(%)
1. If the adjuster banks control is set to auto mode, then the adjuster
banks will move according to the power error versus zone level
plot below.
2. In the 'DRIVE 1 BANK IN' region, the adjuster bank will be driven
in an 'ascending bank number' manner, i.e. only one bank will be
traveling at a time, and the next bank will start to drive in when the
last bank is completely driven in, and the rest of the banks will be
driven in this sequential order.
206
ADJUSTER BANKS AUTO MODE CONTROL
100
LIQUID ZONE LEVEL(%)
DRIVE 1
BANK IN
75
DRIVE 2
BANKS OUT
50
15
DRIVE 1
BANK OUT
-6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6
POWER ERROR(%)
100
ABSORBER/ADJUSTER
TRAVEL SPEED(%)
50
40
-6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6
POWER ERROR(%)
207
· The reactivity change worth of the adjuster and the absorber banks are
assumed to be:
208
9.4 Coupled Reactor Kinetics Reference
1. Luxat, J.C., “ The Potential of a Generalized Modal Analysis of
Postulated Loss-of-Coolant Accidents in CANDU Reactors”, Nuclear
Technology, Vol. 46, Mid-Dec 1979.
209
Typical Workshop Agenda
Session 1 -
Ô Simulation Fundamentals
3. CASSIMÔ
Session 2 -
210
2.2 Photoneutron Sources in Heavy Water Reactors
2.4 Discussion
3.4 Discussion
Session 3 -
Session 4 -
211
3. Model Integration for Simple Reactor
Session 5 -
212
APPENDIX
Selected CASSIMÔ Algorithms Descriptions
Algorithm # 102: PUMP_SIM, Simple Centrifugal Pump
The pump model is a simple centrifugal pump dynamic model which utilizes
the pump's manufacturer performance data for the following computation:
1. The pump static head is calculated as proportional to the product of the
design static head and the square of normalized pump speed. Other
factors such as pump degradation and fluid density changes relative to the
design density are taken into account in the calculation. That means as the
pump is starting up, the pump static head will be increasing and will be on
a quadratic rise to the design static head as the speed increases.
2. The pump power is computed as quadratic function of pump mass flow,
using curve-fit technique on the pump performance curve to obtain
respective constants (KP1, KP2, KP3) for the quadratic function.
3. Heat dissipated by the pump frictional loss is calculated as a fraction of
pump power
4. Cavitation factor is calculated as (suction pressure - saturation pressure) /
N.P.S.H.
Scope of Application:
Physical effects simulated:
· Static head calculation as function of speed, fluid density changes and
pump efficiency.
· Pump power as a function of mass flow.
213
· Heat dissipation as pump frictional loss.
· Cavitation factor is calculated.
Validity check:
· Cavitation factor should be checked if it is near or equal to 1.
Range of Operation:
· Full range of operation with the exception of pump cavitation.
214
Coefficient Specifications (ALG102):
215
Algorithm # 106: VLV_ONOFF_SM, Simple ON-OFF Valve
Scope of Application:
Physical effects simulated:
· Motive force (e.g. instrument air signal) for the valve.
· 3 types of valve can be specified: linear, equal %, quick opening
· Valve can fail open, closed, or failed to desired position.
Limitations or effects not simulated:
· Valve cavitation is not calculated in the model
· If the valve characteristics is known, and it does not fit the linear, equal
% or the quick opening, an extra block(ALG152-CURVEFIT) can be
created to calculate the effective port-area, using the stem position
output(OUT(2)) from this block as an input to the CURVEFIT block.
Validity check:
· Valve will not operate if motive force is not provided.
Range of Operation:
· Full range of valve operation, provided motive force is kept constant and
is greater than the motive force required to fully open and close the
valve.
216
Input Specifications(ALG106):
Output Specifications(ALG106):
217
Coefficient Specifications(ALG106):
218
Algorithm # 108: CV_SIM, Simple Control Valve
Scope of Application:
Physical effects simulated:
· Motive force (e.g. instrument air signal) for the valve.
· 3 types of valve can be specified: linear, equal %, quick opening
· Valve can fail open, closed, or failed to desired position.
Validity check:
· Valve will not operate if motive force is not provided.
Range of Operation:
· Full range of valve operation, provided motive force is kept constant and
is greater than the motive force required to fully open and close the
valve.
219
Input Specifications(ALG108):
Output Specifications(ALG108):
220
Coefficient Specifications(ALG108):
221
Algorithm #128: TANK, Simple Storage Tank with Heating
3. If the level in the tank is greater than the overflow height, overflow of the
liquid will result. A simple first order lag is used to model the overflow.
4. When the level in the tank falls below the elevation at which the outlet
pipe is connected to the tank(KEMPTY), the outlet flow should decrease.
The outlet flow is assumed to be decreased to zero, when the level falls to
20% of KEMPTY. The tank level factor is used as a measure to reduce the
outlet flow when the level in the tank falls below KEMPTY.
Scope of Application:
Physical effects simulated:
· Overall mass balance to calculate the liquid level inside the tank
· Overall energy balance to calculate the liquid temperature inside the tank
· Liquid overflow will occur, when the liquid level inside the tank exceeds
the overflow limit
· Tank level factor reflecting the required reduction in outflow when the
tank liquid level falls below the elevation at which the outlet pipe is
connected to the tank.
Validity check:
· If reverse flows are simulated, additional blocks must be used to calculate
overall input and output flowrates with respect to the tank.
· Check for two phase fluid.
Range of Operation:
222
Special Requirements & Assumptions:
· All fluids are in single phase, no boiling or condensation takes place in
the tank.
· The properties Cp of both fluids are assumed constant within the range of
their temperature range.
· Only 1 overall inlet flow and 1 overall outlet flow is modelled. If more
than one inlet flow or outlet flow exists, additional blocks will be
required to calculate the overall inlet and outlet flows.
· The energy balance for the tank assumed perfect mixing for the tank.
Input Specifications(ALG128):
Output Specifications(ALG128):
223
Coefficient Specifications(ALG128):
224
Algorithm # 133: BREAKER, Breaker
2. The inputs for the breaker algorithm includes close and open command,
trip command, trip relays, and malfunctions.
Scope of Application:
Physical effects simulated:
· Breaker status
· Breaker alarm indications
225
Range of Operation:
Special Requirements & Assumptions:
226
Algorithm # 134: MOTOR, Simple motor model
3. The total current can be calculated based on the motor start-up current and
the motor running current.
5. The motor start-up current is calculated based on the speed of the motor,
and the initial start-up current when the speed is zero.
Scope of Application:
Physical effects simulated:
· Motor speed
· Current used by the motor
Validity check:
Range of Operation:
227
Input Specifications(ALG134):
Output Specifications(ALG134):
Coefficient Specifications(ALG134):
228
Algorithm # 140: CONDUCT_FI, incompressible flow conductance
W = K*VAL*(Pup-Pdown)0.5
Scope of Application:
Physical effects simulated:
· The flow through the link as predicted by incompressible flow equation
Validity check:
Range of Operation:
229
Input Specifications(ALG140):
Output Specifications(ALG140):
Coefficient Specifications(ALG140):
230
Algorithm #223: RSFF, Reset — Set (RS) Flip Flop
The truth table below represents a typical sequence of events for the output
logic states evolution corresponding to Set Signal & Reset Signal with
Coefficient = 0
Scope of Application:
Physical effects simulated: Not applicable.
231
Range of Operation:
· Input values are either 0 (logical false) or 1 (logical true). No
intermediate values between 0 and 1 are allowed.
Special Requirements:
· None.
232
Algorithm #232: Comparator, Single Analog Comparator with Deadband
Scope of Application:
Physical effects simulated:
· The effects of a comparator switch (e.g. level switch etc.)
Range of Operation:
· Input values can assume any values from -99999999.99999999 to
+99999999.99999999
Special Requirements:
· None
233
Input Specifications (ALG # 232):
234
Algorithm #236: DELTA_INP, , Input = (Current Input — Previous Input)
Scope of Application:
Physical effects simulated: Not applicable
Range of Operation:
· Input values can assume any values from -99999999.99999999 to
+99999999.99999999
Special Requirements:
· None
235
Input Specifications (ALG #236):
236
Algorithm #241: SPEEDER_GEAR, Steam Turbine speed/load changer gear
Scope of Application:
Physical effects simulated:
· Simple turbine speed/load changer gear.
Validity check:
· None
Range of Operation:
· All digital input signals assume value either "0" or "1".
· All analog input signals assume normalized value.
Special Requirements:
· None
237
Input Specifications (ALG # 241):
238
Algorithm #312: 5_SUM, Sum of 5 Inputs
Scope of Application:
Physical effects simulated: Not applicable
239
Input Specifications (ALG # 312):
240
Algorithm #332: PID_MALF, Proportional + Integral + Derivative Controller
with Controller Malfunction
241
When RSP mode requested is accepted, Output #8 will indicate that the
controller is in RSP Mode.
In RSP Mode, the Setpoint signal comes from the Model via Input #17.
(d) Tripped-To-Manual Mode - When the controller is in any of the
above control modes, if Input #5 is "true" (= 1), the controller will be
"tripped to Manual". Output #7 will indicate that the controller is in
Manual; whereas Output #27 will indicate that the controller has been
"tripped to manual". Logic Model may send a Manual Control signal
at Input #21, in order to set the controller output to a desired value,
should the controller be "tripped to Manual".
The " Tripped to Manual" mode will be reset when Input #27 is "true".
When that happens, Output #27 will be reset to "0", and the controller
will remain in "Manual" mode ((b) above).
5. There is a provision to "track" a controller when Input #14 (Track
Command) is "true" ( = 1). When that happens, regardless of what mode
the controller is in, the controller output #1 will follow the signal provided
at Input #15 (Track Value). Output #9 will indicate that the controller is in
"TRACK".
6. The controller can be locked in a particular mode:
Input #11 = Lock Auto
Input #12 = Lock Manual
Input #13 = Lock RSP
It is important for the user to make sure that the"lock" signals for the 3
modes (Input #11, #12, #13) cannot be applied simultaneously.
When the mode "lock" signal is accepted, Output #10 will indicate that the
controller is in "lock".
7. The controller can be inhibited from "Auto", "Manual" or "RSP" mode via
respective Input #8, Input #9 and Input #10. The "Inhibit" function will
prevent the inhibited mode from being set. When a mode inhibit function
occurs, Output #12, #13, #14 will display the mode inhibited function
respectively.
8. The controller accepts "feedforward" bias via Input #18 for MMI, or via
Input #19 for Model input.
9. Usually the controller tuning constants: Proportional Gain, Reset Time, Set
Time are set by the respective coefficient #1 , #2, #3. However, this
algorithm provides a tuning provision with which the user can tune the
controller "on-line" by sending a "controller tune" request at Input #22.
When that happens, Proportional Gain can be provided at Input #23; Reset
Time at Input #24; Set Time at Input #25.
10. The controller output can have a low and high limits specified by the user
respectively at the Coefficient #5 and #6.
11. Note that for the purpose of interfacing to Labview GUI, a special
algorithm #457 has been designed to interface with this algorithm. See
details in ALG # 457.
12. This algorithm accepts controller malfunction request via Input #28 (failed
controller flag); and via Input #29 (failed controller output position) If the
controller has failed, then the controller output will be set to equal Input
#29.
242
Scope of Application:
Physical effects simulated: Not applicable
243
Input Specifications (ALG # 332):
244
27. MMI to Reset D Nil See 4(d) above
"Tripped to
Manual"
28. Controller Failed D Nil See 12 Above
Flag
29. Controller Output A Normalized See 12 Above
Failed Position
245
21.Controller Output A Normalized For algorithm internal use
Low Limit (A9)
22.Controller Output A Normalized For algorithm internal use
High Limit (A10)
23.Controller Setpoint A to follow SP For Labview pop-up interface, see connection details
Registered (A11) in ALG # 457
24.Controller A Normalized For Labview pop-up interface, see connection details
Feedforward Bias in ALG # 457
Registered (A12)
25. MMI Manual D Nil For Labview pop-up interface, see connection details
Control Signal in ALG # 457
Update Request (D10)
26. MMI Setpoint D Nil For Labview pop-up interface, see connection details
Update Request (D11) in ALG # 457
27. Controller Tripped D Nil For Labview pop-up interface, see connection details
to in ALG # 457
Manual (D12)
246
Algorithm #398: MARKER, Dummy Marker for Model Identification
Scope of Application:
Physical effects simulated: Not applicable
247
Algorithm #399: DUMMY_DISPLAY, Dummy Display of 20 Output Tags
Scope of Application:
Physical effects simulated: Not applicable
No Outputs.
No Coefficients.
248
Algorithm #401: H_NET_INT, Internal Node
5. Any flows to and from the node that are not related to the network link can
be inputted to INP(1) from the node.
Scope of Application:
Physical effects simulated:
· Pressure in the node will respond to a mismatch between the inflow and
outflow of the node.
Range of Operation:
· network flows remain in a single phase, for 2 phase systems, a 2 phase
friction factor has to be calculated outside the network blocks.
249
Input Specifications(ALG # 401):
250
Algorithm #402: H_NET_EXT, External Node
4. This block can be seen as the interface block to represent the existence of
an external node.
Scope of Application:
Physical effects simulated: None
Range of Operation:
· See CASSIM USER MANUAL section 2.3.4 and CASSIM TUTORIAL
section 2.3 on hydraulic flow networks for further explanations.
251
Input Specifications(ALG # 402):
252
Algorithm #403: H_NET_LINK, Hydraulic Network Link
3. With the pressure difference between the two connected nodes, the link
flow can be calculated using a linearized form of the momentum equation.
W = KCOND*TANK_LVL_FAC*VAL1*VAL2*(Pup+Pelev +Pstatic-
Pdown)0.5
where:
KCOND, the link conductance is calculated based on the link mass flow
and the link pressure drop at design conditions. For cases, where the
momentum equation for the link does not follow the above form, e.g.
compressible flow, or pump head/flow curve, an external conductance
has to be calculated external to this algorithm, and connected to INP(9)
of this algorithm. COF(8) has to be set to one, if external conductance is
used.
VAL1, VAL2, is the normalized valve port-area for the valves in the link
253
Pelev is the elevation head of the piping arrangement in the link
connecting the two nodes.
Pstatic is the static head for the pump located in the link
6. If there are valves located on this link, the normalized valve port-areas
must be connected.(INP(3), INP(4)). These valve port-areas will reduce the
conductance of the link thereby reducing the link flow.
Scope of Application:
Physical effects simulated:
· Flow calculation based on incompressible flow
· Reverse flows can be simulated.
Range of Operation:
· network flows remain in a single phase, for 2 phase systems, a 2 phase
friction factor has to be calculated outside the network blocks.
254
Input Specifications(ALG # 403):
255
Coefficient Specifications(ALG # 403):
256
Algorithm #451: TIMEKEEPER
· Simulation Time
The simulation time is represented as an eight-digit double precision number.
Starting from the left, digits 1 and 2 represent the hour, digits 3 and 4 represent the
minute, digits 5 and 6 represent the second, and digits 7 and 8 represent half
second.
There are two ways to use the simulation time with TIMEKEEPER:
1) Use the PC's current clock time as the initial simulation start time
Code the data value at output 1 as zero.
2) Specify a particular initial simulation start time
Code a non-zero time as the data value at output 1. Use the eight-digit format
as described above.
Each iteration increments the simulation time by dt (the iteration time step
specified in the model).
Scope of Application:
Physical effects simulated: Not Applicable.
257
Input Specifications (ALG # 451):
258
Algorithm #452: NON-MOVEABLE, Non-moveable Block NMB
Model Model
Scope of Application:
Physical effects simulated: Not Applicable.
259
Input Specifications (ALG # 452):
260
01-00344