ML16256A181
| ML16256A181 | |
| Person / Time | |
|---|---|
| Site: | Waterford  |
| Issue date: | 08/25/2016 |
| From: | Entergy Operations |
| To: | Office of Nuclear Reactor Regulation |
| Shared Package | |
| ML16256A115 | List:
|
| References | |
| W3F1-2016-0053 | |
| Download: ML16256A181 (32) | |
Text
3.6B-1 APPENDIX 3.6B PIPE WHIP ANALYSIS (ETR 1002, Appendix C)
WSES-FSAR-UNIT-3 C-2 APPENDIX C of ETR 1002 PIPE WHIP ANALYSIS 1.
METHODS OF ANALYSIS 1.1 Dynamic Analysis The dynamic method of analysis employed to describe the effects (motion, loads, etc.) of pipe whip relies on the execution of two computer programs: RELAP-3 (or RELAP-4) and PLAST, described hereinafter in more detail.
The RELAP program calculates thrust forces generated as a function of time throughout a complex piping system following a break in such systems. For many cases the thrust force will only be required at the break plane, since for the majority of whipping pipe target geometries, the impacting event (on a structure, restraint, etc.) will have taken place before additional thrust forces, differently directed, caused by pressure unbalances in the piping system, have been generated.
The exit plane time dependent thrust forces, or the time dependent thrust forces at various points in the piping system (i.e. elbows, changes in area, etc.), as the case may be, are then used as input to the PLAST program, which analyzes the response of the piping and restraint system to such loads.
1.2 Equivalent Static Analysis Because of the relatively large number of breaks which must be analyzed, as well as the non-uniqueness of any restraining system, which may be utilized, it would be extremely time consuming to dynamically analyze every break with every possible restraint system configuration and permutation.
In order to aid in the design process, portions of the initial restraint system design may be based on an equivalent static analysis.
Following such static analysis, the designed restraining system will be verified for adequacy by dynamically analyzing its response to postulated major breaks with the methods of section 1.1 of this Appendix. If the system is not adequate, appropriate design changes will be made, and the dynamic analyses repeated until adequacy is demonstrated.
In this fashion an adequate, rather than optimum, restraining system will be designed, but in a much shorter time.
WSES-FSAR-UNIT-3 C-3 Revision 11-A (02/02)
The static loads applied in the initial design of the restraints and restrain bearing structures will be in the shape of a step function of amplitude given by
¨(DRN 00-1032)
Fes = K1 K2 PA (1.1)
Where:
A is the break area K2 is a thrust coefficient which has the value of 1.26 for steam and saturated water, and 2.0 for subcooled water, and K1 is a dynamic load factor accounting for the dynamic response of the target to the suddenly applied load and the potential increase in load caused by the acceleration of the pipe through any gap while the escaping fluid thrust force is still acting.
A proper value of K1 can, of course, be obtained only after exhaustive dynamic analysis, and is bound to vary from configuration to config-uration. After a sufficient number of typical cases have been dynamically analyzed, it will be possible to point out an upper bound to the value of K1. Our experience to date, although limited, points a value of 3.0 as sufficient for Ebasco's typical piping restraint configurations wing gaps below 4 inches.
It should be stressed, however, that the choice of K1 would affect period of design rather than the design itself, in that a dynamic verification is always required. Thus an incorrect (too small) value of K1 would be ascertained during the dynamic analysis, and the design process would have to be repeated.
¨(DRN 00-1032)
The pressure in equation (1.1) is normally the maximum operating pressure in the line. For feedwater lines, however, the subcooled decompression is generally very short (of the order of milliseconds or less) because of the large speed of sound. For such lines therefore use of the maximum operating pressure in equation (1.1) would result in a very conservative estimate of the thrust force. More properly for these lines an equivalent static force will be used with a magnitude equal to K12.OPoA at time zero, decreasing linearly to K11.26 Psat A at a time, t = L/C, where L is the length of the line from the break to the reservoir or pump, C is the speed of sound in the fluid, and Psat the saturation pressure corresponding to the maximum operating temperature in the line.
WSES-FSAR-UNIT-3 C-4 2.
DESCRIPTION OF PLAST 2.1 Analytic Description PLAST is a FORTRAN program, based on the "matrix-displacement method" 1/, for three dimensional space frames in which sections of pipe comprise the structural elements. In this method, the overall piping system is subdivided into substructures, called "elements" and these elements are joined to each other at "nodes". The systems mass, consisting of linear and rotary inertia terms, is lumped and located at these nodes. An "assembled" stiffness matrix is then generated by first referring each of the element matrices to a system of "global axes" and then adding appropriate terms of the matrices. A similar method is also applied to the generation of an "assembled" mass matrix.
Following this, a damping matrix may be generated with upper bounds on damping. Upon defining a time dependent forcing function, an elasto-plastic constitutive law and constraints on the systems boundaries, a set of ordinary differential equations in the time domain, are delineated in a form amenable to solution. The details of the above outline will now be presented.
2.1.1 The Element Stiffness Matrix A typical straight pipe element, Figure 2.1.1 is shown oriented to a set of global axes, (X,Y,Z). Local axes, (x,y,z) are indicated at each end node I and J such that the x axis is positive from J to I.
The positive direction of both sets of axes are used to define positive displacements, forces, accelerations and velocities. The flexibility matrix is commonly used to derive a stiffness matrix. 1.2/ This stiffness matrix is square, symmetric and 12 X 12; its dimensionality corresponding to the 12 degrees of freedom available to each pipe element. If the transposed displacement and force vectors for an element are defined as:
{ }
x x
y z
x y
z x
y z
x y
z
= [(
) I (
) ]
T J
{ }
)
) ]
f = [(f I (f T
x x
J f f M M M f f M M M y
z x
y z
y z
x y
z where lower case subscripts refer to local element axes along which vector components are directed, then the following relationship holds for [k], {x} and {f}.
[k] {x} = {f}
(2.1.1.1)
The orthogonal transformation matrix [T] that rotates the (X,Y,Z) axes into the (x,y,z) axes also rotates displacements, forces, velocities, and accelerations.
WSES-FSAR-UNIT-3 C-5
{x} = [ T ] {X}
{f} = [ T ] {f}
{}
x = [ T ] {X}
{}
x = [ T ] {X}
(2.1.1.2)
The transformed stiffness matrix (referred to global axes) is then found by substituting the first two of equations (2.1.1.2) into (2.1.1.1) and premultiplying both sides by [T]T. Then since
[T]T [T] = [ I ] (the identity matrix)
([ T ]T [ k ] [T]) [F]
(2.1.1.3) and the transformed stiffness matrix is defined as:
[ k ] = [ T ]T [ k ] [ T ]
2.1.2 The Element Mass Matrix The element mass matrix referred to local axes is a 12 X 12 diagonal matrix in which are defined the linear and rotary inertia terms; m = w L / 2g (2.1.2.1) xx I
m A
xx I
zz I
A yy I
m yy I
=
=
3/
2.1.3 Damping Matrix PLAST utilizes a damping matrix that is a linear combination of the stiffness and mass matrixes. In doing so, constants may be determined such that a specified damping factor is bounded within a required frequency range. This under estimate of the system damping is conservative in regard to estimating system responses to pipe rupture blowdown forces.
2.1.4 Elasto Plastic Constitutive Laws Constitutive laws governing the elasto-plastic behavior of three dimensional space frames were developed in reference 4/ for elastic-perfectly plastic materials based on the Von Mises yield criteria. 5/
WSES-FSAR-UNIT-3 C-6 Revision 11-A (02/02)
As in two and three dimensional continua, a "yield surface" (or "interaction surface") in force space is developed to define conditions under which yielding may take place. Because of the "nonuniqueness" of elasto-plastic solutions; ie, because of solution dependency on path, and "incremental" constitutive law, relating plastic displacement increments to the internal force states (luring yielding, is applicable. 5/ Incremental constitutive laws have the further property that the load sequence of yielding, unloading, rigid body rotation and yielding again result in continuous states of plastic deformation. 5/
The yield surface for Von Mises materials is defined by a right circular cylinder in "stress-space," with axes having equal direction cosines with all principal stress axes. 5,6/ The surface is defined mathematically as;
¨(DRN 00-1032)
F('ij)-e2/3 = 0 (2.1.4.1)
The yield condition then states; If F('ij)-e2/3>0and F('ij) > 0 then yielding occurs If F('ij)-e2/30and F('ij) < 0 then action is elastic If yielding does occur then plastic strain increment is given by:
ij ij s
)
(
F
=
p ij (2.1.4.2)
(DRN 00-1032) and this strain vector is oriented outward and normal to the yield surface; ie, s is always positive.
In a similar manner a yield surface in force space is represented as:
- 1 = 0 (2.1.4.3) where
= (fiP)
(2.1.4.4) fiP (i = x, y, z) = generalized moments and forces.
This is the equation of a unit sphere in a generalized force space of generalized forces. The increment in generalized plastic displacements of a beam section undergoing yielding is a vector normal to and pointing outward from the yield surface in generalized force space.
{
}
f X
f P
=
(2.1.4.5)
WSES-FSAR-UNIT-3 C-7 Revision 11-A (02/02) 2.1.5 Strain Hardening Of the available strain hardening theories, that of isotropic hardening 5/ is the simplest mathematically and requires the least data on material properties. In essence, this theory predicts that the yield surface in stress space expands isotropically about its origin as a material yields. In order to develop strain hardening parameters for this approach, the following assumptions are made:
a -
Yielding occurs simultaneously throughout an entire cross section or not at all.
b -
Small deformations occur.
c -
Because beam elements can be loaded only at their nodes plastic hinges occur only at nodes-the exception being when both nodes yield in such fashion as to cause the entire beam to yield.
As in elastic analyses, material properties are obtained from one dimensional tensile/compression tests. Referring to Figures 2.1.2 and 2.1.3, an "effective-stress-strain" curve can be developed by plotting the plastic strain component, (p) against stress,. By defining; 2
e J
3
=
(2.1.5.1)
}
)
(
1/2
+
)
(
{
2/3
=
2 p
x 2
p x
p e
(2.1.5.2)
¨(DRN 00-1032) as the effective stress, and the increment in effective plastic strain respectively, it maybe observed that there is a one to one correspondence between ((e, e) for the effective-stress-strain curve and (, p) for the uniaxial case. An expression for f at some node may now be found for a time interval tn, using the stress state at the beginning of the interval, the displacement incurrent and the effective stress strain curve.
(DRN 00-1032) 2.1.6 Equations of Motion - Solution Method The equations of motion for the system of lumped masses tied to each other by a flexural piping system is written below for elastic displacements;
{F}
=
{X}
[K]
+
{X}
]
C
[
+
{X}
]
M
[
e (2.1.6.1)
WSES-FSAR-UNIT-3 C-8 Revision 11-A (02/02)
The superscript "e" is applied only to displacements to indicate "elastic displacements" are intended. This "elastic" designation for displacements derives from the observation that the force displacement relationship;
[K] {X} = {F}
applies only to elastic media. By artificially subdividing, displacements into "elastic" and "plastic" components, a plastic "correction" force arises; ie,
[K] {X - Xp} = {F}
and
¨(DRN 00-1032)
[K] {X} = {F} + [K] {Xp} + {F} {Fp}
(2.1.6.2)
Hence equation (2.1.6.1) may be rewritten as;
[M] {X} + [C] {X} + [K] {X}e = {F} + {Fp}
(2.1.6.3) where {F}, {Fp} are the externally applied force vector and the plastic correction force vector respectively. The system of second order ordinary differential equations in the time domain is solved by a step-by-step, numerical integration method developed by Newmark 7/
Appendix C. An internal procedure of checking convergence and stability is used in conjunction with a variable integration step. If instability is detected, the program automatically aborts and prints diagnostic messages. A brief algorithmic flow chart, Figure 2.1.5, is included to indicate the sequence of operations towards reaching a numerical solution.
Boundary Conditions In a well posed boundary and initial value problem, displacement and force specifications at boundaries may be imposed subject to the following restrictions; ie, - displacement and force specifications may not be made for the same boundary point, in the same direction simultaneously. Violation of this restriction results in an overspecified problem. The pipe whip problem involves two general types of displacement and one type of force boundary specification.
In the first type of displacement boundary condition, one or more displacement components vanish at a boundary point. This corresponds to an anchor or restraint fixed for all times. The second type of displacement boundary condition is displacement dependent, as indicated below.
WSES-FSAR-UNIT-3 C-9 Revision 11-A (02/02)
Let
~ be a prescribed displacement vector from some point on a boundary in its undisturbed, initial state (Figure 2-14); let
~n be the subsequent displacement of that point at time tn. Then if 0
~
~
n or (2.1.6.3)
~
0
~ n and
~
~
n the displacement of the boundary point is unconstrained.
If
~
> 0
~n
¨(DRN 00-1032) and
~
~
n (2.1.6.4)
(DRN 00-1032) then one of the following constraints may be imposed on the displacement of the boundary point; (a)
~
=
~
~
~n
(b)
~
~
~n
=
( )
n
~
(2.1.6.5)
¨(DRN 00-1032)
(c)
)
~
(
N
=
~
~
p n
e n
n nf
~
where, N designates linear and nonlinear functions respectively of their arguments; and the dot over the "
n
~ " argument designates "the increment of n ".
WSES-FSAR-UNIT-3 C-10 Revision 11-A (02/02)
The first constraint condition (a) represents a rigid restraint, with gap; condition (b) represents a linear spring restraint, with gap and condition (c) represents an elasto-plastic, strain hardening restraint, with gap. All restraints are aligned with the vector ~ which represents the initial gap.
The effect of such displacement dependent constraints is to impose a sudden impulse on a boundary consistent with the constitutive laws of the restraint materials. The resultant change in momentum of a boundary mass accounts for its instantaneous rebound velocity as well as the peaks in the restraint reactions.
Obviously, gap size affects the boundary mass velocity just prior to impact, while energy dissipation due to damping and yielding in both the pipe and restraints affects the momentum exchange, mass rebound velocity, and peak restraint reactions. The force boundary condition consists of the time history of the blowdown force at the pipe rupture. This force is always kept tangent to the pipe axis for guillotine breaks and normal to it for slot breaks, as the pipe deforms and rotates.
2.1.7 Kinematic Considerations The small deformation theory, used in PLAST, is based on the following considerations with regard to the kinematics of flexible piping systems;
¨(DRN 00-1032)
The axial strain in a pipe undergoing finite deformations is given as:
2, 12/
2 z
2 y
2 x
2 2
2 2
z x
x u
+
x u
+
x u
1/2
+
)
u r(
+
x u
=
x u
x y
xx (2.1.7.1)
¨(DRN 00-1032)
This expression differs from that for small strains only in the last, bracketed term. For flexible piping systems, only the last two of the squared terms in brackets need to be considered as having a significant effect on strain. In order to differentiate between "small" and "large" strains it is observed that for large strains the bracketed, second order terms in equation (2.1.7.1) are of the same order of magnitude as the sum of the remaining terms. Now the significant second order terms, above, represent rotations about both bending axes.
Thus, for example, if a pipe yields at (xx = +/-.001) then these rotations must be of the order of 2(.001)1/2 radians or +/-3.62 degrees, before second order terms begin to become significant.
The pipe system is monitored by PLAST to determine if rigid body rotations of the nodes exceed some preset magnitude. This occurrence would imply that a stage has been reached in the dynamic analyses in which the dropping of second order terms from the strain tensor would result in detectable error. At this point a slow divergence in solutions obtained, using small and large deformation assumptions would occur with respect to monitonically increasing deformations.
WSES-FSAR-UNIT-3 C-11 Revision 11-A (02/02) 2.2 Discussion and Summary of Program Features and Characteristics PLAST solves initial value problems in the stress analysis of piping systems that are nonlinear in materials and boundary conditions. This is accomplished by application of a numerical integration method in the time domain to a system discretized in space. Speed and accuracy of solutions depend in large part on the fineness of discretization, intrinsic physical properties of segments of the piping system and on PLAST itself. It is noted in Appendix C that finer discretization and larger maximum stiffness lead to shorter time steps and longer computer execution time. This is characteristic of all numerical integration methods used in the solution of initial value problems. However, this attribute is not significant, except with regard to the time interval over which solutions are of interest.
The time interval of solution interest varies with each problem and is not sharply defined. The upper bound of this interval is determined by a set of criteria denoting some sequence of events that must occur. The program user in part retains the option of judging whether or not these criteria are satisfactorily met.
For example, in executing a dynamic pipe whip analysis with pipe whip restraints, the following criteria are used to denote solution termination;
¨(DRN 00-1032)
(DRN 00-1032) a -
The system response reaches steady state. This is defined by the motion of some characteristic point of the system, such as the point of load application or a restraint point.
When this point oscillates with nonincreasing amplitudes about some median displacement, a steady state condition has occurred.
b -
Some point of maximum strain in the piping system or pipe whip restraint has exceeded v/2.
c -
At a given point, allowable deflection is exceeded and inter-ference with some safety related object occurs.
d -
At the location of some component, essential to safe shutdown such as an isolation valve, rotations, deflections and internal forces have exceeded their allowable values.
The user exercises his option of continuing or terminating the program by reviewing the final printout. Maximum strains, rotations, deflection at applied load, reactions at restraints, deflections at restraints and their times of occurrence are monitored by PLAST and printed out at each print interval as an aid to the user. If the upper bound of the time interval of solution needs to be extended, the user exercises a "restart" option and execution continues to the new specified upper bound in time.
At regular, specified intervals the entire problem solution is printed out. This includes time, displacement, velocity, acceleration, external reactions, internal forces, elastic and plastic strains and stresses corresponding to each degree of freedom for each node and each
WSES-FSAR-UNIT-3 C-12 Revision 11-A (02/02) element. The updated maximum parameters mentioned in the preceding paragraph are also included. Externally applied forces at a node, and external spring forces applied to a node appear in the column titled "REACTION" in the printout. Sundry details regarding I.etc, are described in the user's manual 8 / for PLAST and will not be described here.
2.3 Typical Pipe Whip Dynamic Analysis 2.3.1 Steam Line with Yielding Restraints The 43 in. steam line shown in Figure (2.3-1) was analyzed for a guillotine break at node 16 and the blowdown force of Figure (2.3.4). Effective gaps, restraint orientations, and spring constants are indicated on Figures (2.3.2 and 2.3.3) for the restraint located at node 14. The restraint deflection (Figure 2.3.5) is taken as an indicator of steady state behavior and this is seen to occur at about 0.13 seconds past blowdown initiation. The restraint is seen to permanently deform to about 8.5 in. and to oscillate about this point. No yielding occurs in the pipe, however, the continuance of blowdown beyond 0.3 seconds keeps the pipe in oscillatory conduct with the restraint. (Figure 2.3.6) 2.3.2 Steam Line with Elastic Restraint
¨(DRN 00-1032)
The 34 in. line in Figure (2.3.7) is shown with six restraint locations. A guillotine break occurs at node 12 and the applied blowdown force is given in Figure (2.3.10).
Oscillatory motion of the load point, Figure (2.3.8) indicates onset of steady state conditions and starts at 0.1106 seconds past blowdown initiation. Yielding occurs in the pipe at the point of impact (node 9). Restraint reactions form sharp spikes and unload rapidly. In this case pipe impact occurs in two directions (Figure 2.3.9) at the same restraint. Spring constants and gap data are given below for each restraint in all four spring directions.
Node #
Gap Elastic Spring "K" (DRN 00-1032) 2 6."
8.536 X 106 #/in.
4 6."
5 6."
8 5.5" 9
5.5" 11 4."
WSES-FSAR-UNIT-3 C-13 Revision 11-A (02/02) 3.1 PLAST Verification of PLAST was made by comparing solutions by PLAST against analytic and numerical results from other methods and programs. These include; (a) an elastic simple beam dynamic analysis by modal superposition (b) plane frame, elastic and elastic-plastic numerical results (c) a simple beam elastic-plastic analysis.
3.1.1 Simple Uniformly Loaded Beam The equations of motion for the beam illustrated in Figure (3.1.1) are:
2 2
2 2
x y
x y
t
) + m
) = p (x,t)
(EI 2
2 (3.1.1) and its solution if found to be (see Reference 9/ for example) y x t
(, ) =
(x) q (t) n n=1 n
(3.1.1.2) where qn, n are the mode displacements and mode shapes respectively and are shown below as; qn n
=
2 L (+1 - (-1)
(t -
n )
Mo n
sin t
n
)
2 (3.1.1.3)
¨(DRN 00-1032)
)
L x
n
(
sin 2
=
n (DRN 00-1032) and
n
=
(n
)
EI Mo L 3 2
¨(DRN 00-1032)
(The values of E, I, L, and Mo are given in Figure (3.1.1)
WSES-FSAR-UNIT-3 C-14 Revision 11-A (02/02)
The solution to this problem is displayed in Figures (3.1.2 A&B) as displacement and acceleration time histories respectively. These solutions representing the superposition of ten non-null mode shapes, are compared with the numerical solution obtained by PLAST, plotted on the same axes. It may be noted, that of the two solutions, that of PLAST is the smoother since it contains rotational as well as displacement codes. In most instances, comparisons of numerical and analytic solutions are in excellent agreement in (displacements and in poorer agreement with time derivatives of displacement as these derivatives approach higher orders. The convergence criteria used in PLAST have, however, resulted in excellent agreement of accelerations; ie, second order derivations of displacement with respect to time.
3.1.2 Elastic Plane Frame Analysis Figure (3.1.4A) and (3.1.4B) display time histories of displacements of points on the frame shown in Figure (3.1.3). These displacements were obtained by PLAST using 5 and 10 lumped masses and compared with similar data obtained from the program of Reference 10 / (which used 10 lumped masses). In comparing these results, the following should be noted;
¨(DRN 00-1032) a -
At a time +/-0.1 seconds past initiation of load application, rotation of the frame at the point of application of horizontal load exceeds 2(.001)1/2 radians. In view of the remarks made in Section 2.1.8 of this report, it would appear that omission of second order strain terms for solutions beyond this time leads to solution divergence for both PLAST and Reference 10 /
as deformations increase.
(DRN 00-1032) b -
PLAST permits the applied loads to rotate with the point of application whereas Reference 10 / does not; hence increasing deviations in the results of these two programs appear with increasing magnitudes of rotation.
¨(DRN 00-1032) c -
An increase in the number of lumped masses in PLAST, from 5 to 10, increases the number of frequency modes available to the discrete system. Although the solution accuracy theoretically improves, there appears to be very little difference between solutions up to +/- 0.125 seconds. Subsequent differences are then due to the nonlinear effects of load rotation.
(DRN 00-1032) 3.1.3 Elasto-Plastic Plane Frame Analysis The same frames shown in Figures (3.1.4A) and 3.1.4B) was analyzed elasto-plastically under the same load for 7 and 10 lumped masses. However, in the case the masses were arbitrarily multiplied by (25)2 which results in a minimum system period 25 times greater than that of the case discussed above. The effect is to filter out system responses to high frequencies in the forcing functions.
WSES-FSAR-UNIT-3 C-15 Revision 11-A (02/02)
¨(DRN 00-1032)
Comparisons are shown in Figures (3.1.5A&B) for the same problem run through ANSYS 11 / using zero strain hardening and FAP 12 / using 10 percent hardening. Again, it is noted that the frame rotations at the point of application of the horizontal load exceeds 2(.001)1/2 at a time
+/- 0.42 seconds from initiation of load application and that the omission of second order strain terms from this time on leads to increases in solution divergence. Nevertheless, solution agreement among these three programs remains remarkably close, well beyond this time.
The effect of the strain hardening is seen to be important only in the large deformation range, and has the effect of limiting the maximum frame deflection after the load has reached steady state. The deviation noted between PLAST and ANSYS for zero hardening occurs in the large deformation range where the effect of applied load rotation (in PLAST) is noticeable. The small variation in deflection curves for PLAST resulting from the increase from 7 to 10 nodes is due to the larger number of available mode shapes in the elastic range and to the formation of an additional hinge (in the horizontal member) in the plastic range. Physically, the additional hinge indicates the spread of yielding along the beam from the vertical load outward.
3.1.4 Elastic-Plastic Beam Subject To An Initial Velocity
¨(DRN 00-1032)
The beam in Figure (3.1.6A) is shown under two condition of loading; the first is an imposed velocity Vo at time to = 0 uniformly distributed along its length and the second is a half sine curve distribution of an imposed velocity U* at time to = 0. In both cases, the imposed velocity drops to zero for all times t > to to while the beam velocity assumes some distribution and shape in accordance with the beam's response.
In Reference 13 /, the author have show that if:
Ui and U i are two initial velocity fields and the stresses and displacements are zero initially for both cases and a quantity o remains shall compared to the initial kinetic energy associated with U i then:
a -
the difference in solutions corresponding to Ui and U i are initially small and cannot increase. Further, if plastic deformation takes place then; (DRN 00-1032) b -
the difference in the above solutions will decrease.
The quantity o, above is the difference in the initial kinetic energies associated with each velocity field.
It is shown in Reference 13 / that for o to be a minimum for the two velocity fields described above,
WSES-FSAR-UNIT-3 C-16 V* = 4 Vo
Non-dimensional plots of displacement and velocity responses of the beam at its center line for each load condition are shown in Figures (3.1.6A) and (3.1.6B) corresponding to the dimensionless parameter, where,
= m EI (V*)
M 2
o 2
Mo = fully plastic moment in beam m = beam mass/unit length E,I = Youngs modulus and area moment of inertia of beam Figure (3.1.6A) compares the central velocities of the simple beam versus time for the two load conditions above as given in Reference 13 / for a rigid plastic beam and as calculated by PLAST for an elastic-perfectly plastic beam. Note that for rigid plastic materials, beam motion ceases after particle velocity reaches zero. In elasto-plastic materials some residual elastic energy is available for the beam to continue to oscillate past this point. The authors contention that solution differences decrease with time or are non-increasing is clearly demonstrated.
Figure (3.1.6B) compares central deflection of the above beam versus time for reference 13 / and PLAST, subject to the sinusoidal initial velocity distribution and elastic-perfectly plastic material properties. The two plots virtually coincide over the region indicated and differ very slightly in peak values of U.
4.
DESCRIPTION OF ANALYTIC METHODS OF RELAP 4.1 General RELAP-3 14 / is a FORTRAN program which describes the behavior of hydraulic systems under transient conditions. The program itself, which has been well documented elsewhere 14,15,16 /, was developed specifically to describe the behavior of water-cooled nuclear reactors during postulated accidents, such as loss of coolant, pump failure, or power transients. RELAP is however sufficiently versatile to describe other, simpler hydraulic systems, not necessarily confined to water fluids.
The system being studied is modeled as an assembly of user defined control volumes connected by flow paths. Check valves, valves, pumps can be inserted at any flow path.
WSES-FSAR-UNIT-3 C-17 Revision 11-A (02/02)
To simulate a break, a leak can be opened either instantaneously, or as a function of time, at the appropriate flow path in the system. The transient energy, momentum, and state equations are then solved for all the control volumes and flow paths to determine the fluid properties anywhere in the system including at the break.
4.1.1 Conservation Equations The following is a brief description of the fluid equations used in RELAP-3. A more extensive discussion can be found in several excellent references which are available on thermal hydraulics 17,18,19 /.
In addition to solving the equation of state (steam tables), RELAP makes use of the mass, momentum, and energy conservation equation a -
Mass Conservation
¨(DRN 00-1032)
)
v
(
=
t (4.1) b -
Energy Conservation
)
v
(
)
v (p
q
)
v
(
=
)
(
e t
e (4.2) c -
Momentum Conservation
=
g
+
)
v v
(
t
)
(
v (4.3)
Under certain general assumptions, which obviously delineate the limitations of the RELAP3 program as presently written, the conservation equations can be further simplified.
WSES-FSAR-UNIT-3 C-18 Revision 11-A (02/02)
The assumptions made in RELAP3 are:
a - Stationary control volumes b - One-dimensional flow c - Negligible body forces (except gravity)
The mass equation can be integrated over a fixed control volume i, thus
¨(DRN 00-1032) dV
)
v
(
=
dv V
i V
i t
(4.4)
Gauss Divergence Theorem further simplifies the above equation to
dS v
=
v
i V
i V
d dt d
(4.5)
But i
i V
dv
is just the total mass in volume i = Mi and
dS V
iS is the total mass flow rate in and out of volume i, across its surface, i.e. its flow paths, then
j ij W
=
t
)
(
d M
d i
(4.6)
The energy equation is simplified ignoring frictional, kinetic and potential energy effects. Thus
)
q
+
v h
(
V
=
t
)
(
d u
d (4.7)
Integration over control volume i, results in the following
S d
h iV
)
q
+
v h
(
=
dV
)
q
+
v
(
V
=
dV dt u) d(
i i
S V
(4.8) i ij j
ij Q
+
h W
=
dt dui (4.9)
The momentum equation restricted to one dimensional flow and integrated over the control volume i, yields with the aid of the Divergence Theorem dV g
+
dS dS
)
S d
v
(
v
=
V d
v i
i i
i i
V S
S S
V
F dt d
(4.10)
WSES-FSAR-UNIT-3 C-19 Revision 11-A (02/02)
For the friction term, the pressure losses are assumed to be proportional to v2, ie:
¨(DRN 00-1032)
A
)
(
)
(
k'
=
)
dS v
(
v
+
.dS F
j j
j j
j j
j j
S A
W W
k v
v i
=
(4.11) where j
W
=
j) v
(
is the average flow rate, and k is a friction coefficient which is directly related to the Fanning friction factor for a single pipe:
(DRN 00-1032) k = f Aw / 2A (4.12)
Momentum flux terms are included as part of the friction term k. This simplification is correct only where the flow is unidirectional during the transient 20/. RELAP will be further modified to properly account for these momentum flux terms 21/ (RELAP-4).
¨(DRN 00-1032)
With this simplification the first term in the momentum equation reduces to dt dWi il
, thus dZ gA
+
144g A
)
P (P
W dt dW j
1 j
c j
1
+
i i
j j
j
+
=
i Z
Z j
j j
i A
W k
l (4.13)
The last term account for gravity effects in vertical control volumes.
Hence the momentum equation reduces to Z
144 1
+
)
P (P
dt W
d 144g i
Z 1
1
+
i i
j j
c j
d W
W k
A j
Z j
j j
j
+
=
l (4.14)
Therein Aj jk 2
c g
144 k'
=
(4.15)
Pi and Pi+1 are the inlet and outlet pressure of a control volume in psi, and Zi+1 and Zi, are the elevation of two contiguous control volumes center of gravities.
Since it is possible to have a pump at a flow path, a change in pressure at a flow path j separating two contiguous control volumes i and i + i, is added to the momentum equation (DRN 00-1032)
WSES-FSAR-UNIT-3 C-20 Revision 11-A (02/02)
¨(DRN 00-1032) j pj 1
+
i i
j 1
144 1
+
P
+
P P
=
144 1
j j
j j
Z j
c W
W k
dz dt W
d A
g i
Z i
+
l (4.16)
The junction inertia for a flow path j is given by
)
A 1
+
i
+
i
(
2 1
=
)
A
(
j l
l l
j (4.17)
The friction coefficient can either be computed prior to starting the program, or can be calculated internally by the program if flow exists at the flow paths, from the steady state momentum equation.
Equations (4.6) (4.9) and (4.16) together with the equation of state are solved using a forward finite difference approximation, with each conservation equation being advanced independently for a small time step.
4.1.2 Choice of Time Step Several considerations enter in the choice of suitable time steps chosen for the forward finite difference solution of the conservation equations, namely numerical stability of the solution, potential instabilities associated with paradoxical emptying a volume of all material or overfilling a volume, and ability of approximating the wave equation describing the preparation of acoustic disturbances throughout the system.
Numerical stability of the solution is normally assured by choosing the time step to be a fraction of the natural oscillation period of the control volume having the smallest volume and flow path inertia. The natural period of such volume is given by 2 C
V A
l where C is the speed of sound in the medium. Thus
¨(DRN 00-1032)
A V
C 2
t l
(4.18)
In addition, either under conditions of very high flow and small mass content or energy content within a volume, to avoid complete replacement of fluid within the volume t
Wj M
j
(4.19)
WSES-FSAR-UNIT-3 C-21 Revision 11-A (02/02) and t
hj U
(Wj j
)
(4.20)
For a one dimensional system the continuity equation may be written as
¨(DRN 00-1032) x v
x v
=
t (4.21)
If the gravity and viscous effects are ignored, the momentum equation can be written as x
v x
v v
2
=
t v
+
t v
2
x (4.22)
Assuming that for each time step the pressure is a function of density and enthalpy, P = P (,S), ie: esentropic process involving small pressure changes within the time step, then (DRN 00-1032) x C
=
x
=
x 2
s
(4.23) and t
C
=
t
=
t 2
s
(4.24)
Substitution of the above in the continuity and momentum equations yield:
¨(DRN 00-1032) 0
=
x v
+
x v
+
t 1
2
P C
(4.25) and 0
=
+
x v
+
x v
v 2
+
t v
+
t P
v 2
2 x
C
(4.26)
Multiplying the first equation by v and subtracting the resultant equation from the second yields 0
=
t v
+
x
+
x v
v
(4.27)
WSES-FSAR-UNIT-3 C-22 Revision 11-A (02/02)
Thus
¨(DRN 00-1032) 0
=
x v
+
x v
+
t P
1 2
C (4.28) 0
=
x 1
+
x v
v
+
t v
(4.29)
If
t P
x v
(4.30) and
t v
x v
v (4.31) then 0
=
x v
+
t P
1 2
C (4.32)
v t
= 1 P
x
= 0 (4.33)
The two equations above can be reduced to the familiar wave equation
2 2
2 P
t P
x
= C 2 2
(4.34) describing pressure disturbances in a medium, by differentiating the first equation with respect to the time variable, the second equation with respect to the space variable, and subtracting the two. The two inequality conditions above, provide a guide for the choice of time steps necessary to describe acoustical phenomena in the fluid system. Normally it is necessary to execute a problem to verify that the inequalities are satisfied.
4.1.3 Thrust Force Equation The general analytical expression for the force on a fluid container can be obtained from the general vector form of momentum equation 18 /. For the control volume depicted in Figure 4.1, the momentum theorem states that they hydraulic force is given by
¨(DRN 00-1032)
)
v
(
v
+
v
+
F
+
=
S V
ext S
ds dV t
dV ds P
F V
=
(4.35)
WSES-FSAR-UNIT-3 C-23 Revision 11-A (02/02)
¨(DRN 00-1032)
In the above equation the term ds P
S
represents surface forces, and hence it can be written as ds Pnds Pnds Pnds ds P
S
w w
2 1
S S
S S
=
(4.36)
Where S being the total surface of the fluid element, can be divided into the wetted surface, Sw and the surfaces across which flow occurs, S1 and S2.
The term ds w
S
is due to friction between the wetted surface and the fluid, with
representing the pressure for the wall due to the shear stress in the fluid.
The unit vector n is taken normal to the surface considered.
The term Fext represents any external force acting on the control volume, and the term dv V
=
represents body forces such as due to gravity.
Forces on the wetted surface are the only mechanism for exerting forces on the walls of the container. Thus the force on a fluid container is given by:
ds Pnds F
w w
S S
+
=
(4.37)
Substitution of (4.37) in (4.36) and (4.35) yields
¨(DRN 00-1032) v
+
F
+
)
v
(
v v
t
=
V ext S
S S
2 1
d x
P Pnds Pnds ds dV F
V
(4.38)
The last term in equation (4.38) is small in blowdown studies and can thus be neglected. The remainder of the equation can further be simplified, if the quantity of interest is the fluid thrust at the break plane of a broken piping system, rather than the combined thrust for an overall system.
This is generally the case for pipe whip studies where the quantity of interest is the fluid thrust normal to the broken pipe end.
In this case the choice of a suitable control volume surrounding the break plane results in the two terms Pnds S1
and Pnds S 2
to approximately cancel each other out, since the normal vector is in opposite directions, the cross-sectional areas are the same, and the static pressure nearly the same, as long as discontinuities do not
WSES-FSAR-UNIT-3 C-24 Revision 11-A (02/02) appear at the exit plane. If there is a discontinuity at the exit plane, as is the case for chocked flow, then the pressure terms will indeed contribute to the force, as will be shown later.
The break plane control volume for circumferential and longitudinal breaks are shown in figures 4-2 and 4-3 respectively.
If the momentum theorem 18 / is applied to an approximate control volume surrounding the break, then the thrust at the break plane is evaluated from:
¨(DRN 00-1032)
v
)
v
(
dt d
=
i i
S V
ds V
d F
(4.39) 4 4
4 3
4 4
4 2
1 chocked is flow if only Pnds s
Pnds s
2 1
The minus sings signify that the thrust on the pipe is in a direction opposite to the flow out of the break.
Therein F is the force exerted on the pipe.
RELAP being a dimensional program, the momentum theorem in one dimension reduces to pnds Pnds ds dV g
v dt d
F c
Vi 2
i i
S S
S c
2 g
v
+
=
(4.41)
The second term represents the outflow of momentum out of the control volume, and for a break area A, it is simply given by Tm =
Av g
=
W Ag 2
c 2
c
(4.41)
The first term is an acceleration (or deceleration) term which expresses the rate of increase (or decrease) in mass flow rate within the control volume (unsteady flow).
This term is approximated in the RELAP program by the following T
W t
a t
= W g
t +
t c
l (4.42)
¨(DRN 01-1032)
If one assumes that a break develops instantaneously, then at time t = 0, there is an unbalanced pressure force at the break plane due to the difference between the fluid pressure and atmospheric pressure.
This unbalanced force must be converted to a fluid acceleration overcoming the fluid inertia, so that fluid begins flowing out of the break.
WSES-FSAR-UNIT-3 C-25 Revision 11-A (02/02)
The RELAP momentum equation for the condition of no flow (t=o), and ignoring gravity effects can be written as
¨(DRN 00-1032) 144
)
(P
=
)
W
(
g 1
i t
+
t c
j
+
i t
j P
j t
W A
l (4.43)
Herein Pi stands for the pressure inside the control volume (fluid pressure) and Pi+i the pressure in the next control volume, in this case none, thus atmospheric pressure.
The thrust force at time t=o, when there is no flow, and which is given by 144 (P-Pa)A, does then indeed equate the acceleration component Ta.
Unbalanced external forces at the break area occur at times other than zero, when the flow out of the break becomes choked. With choked flow there exists at the break plane a discontinuity in pressure which leads to a thrust force component equal to A
144
)
P (P
=
Pnds
+
Pnds
=
a e
S S
2 i
T P
(4.44) where Pe is the exit plane pressure in psi.
The total thrust force at a break location (flow path) is given by the sum of the three components T
A g W
t j
c jt
= W
+ 144 (Pe - Pa) A + (W g
j j
j jt +
t i
c 2
) l (4.45)
The RELAP program 14 / has been modified to solve for the three components and the total thrust force. Three additional minor edit variables have been added to the printout capability described in Reference 14 /. These minor edit variables correspond to the three components of the thrust; the leak force variable LF, is the total thrust force.
The program may also be used to predict overall thrust forces in a complex system, by properly adding the thrust component as calculated from pressures in the system. These pressures are determined by the program. The RELAP-3 program is currently being modified to perform the thrust component addition in the three directions. For the time, however, this is done manually.
4.1.4 Choked Flow Calculations The RELAP program chooses as flow through a junction, the smaller of the inertial flow or choked flow.
WSES-FSAR-UNIT-3 C-26 The limiting mass flow is determined by use of Moodys two phase choked flow model 22 /. This model also allows computation of the throat pressure (taken as the exit plane pressure) as a function of stagnation pressure and enthalpy.
Simon 23 / has discussed the many models which may be used to predict blowdown flow rates of two phase fluids, such as would result following an initial subcooled blowdown. His comparisons among limited test data show that for sufficiently large diameter pipes, and for breaks located a sufficient distance from the reservoir, Moodys modeling may be expected to predict conservative flow rates. For breaks at very short distances from the reservoir, and small diameters, the model might predict too low flow rates.
Comparison of RELAP calculations of pressure variations within a vessel initially filled with subcooled liquid water, following blowdown 24 /,
show that a Moody flow multiplier of 0.6 must be used to predict the correct pressure behavior. The particular experiment reported in 24 /
consisted of a small vessel blowing down through a relatively large diameter and relatively long pipe. It was to be expected therefore that critical flow rates predicted by Moody would be excessive.
Figure 4-4 shows the good agreement of vessel pressure calculated by RELAP with a Moody flow multiplier of 0.6 and experimental vessel pressure 24 /.
Since for different reservoir piping configurations, it is not always possible to ascertain a prior whether a 0.6 flow multiplier or the unmodified flow will result in a higher thrust, RELAP is normally run in both cases, and the higher thrust force then chosen.
4.1.5 Typical Results There are unfortunately no large scale blowdown thrust values reported in the literature. Thus results from RELAP have been compared to those published for small scale blowdown experiments and to theoretical values predicted by Moody 25 /.
Figure 4-5 shows a comparison of thrust calculated using results from RELAP to the thrust measured by Hansen 26 / for the apparatus shown in Figure 4-6. The agreement is rather good.
Figure 4-7 and 4-8 also compare calculated thrusts with Hansens measurement.
In these cases an orifice was present in the discharge pipe of the apparatus reducing flow area to 30 percent and 10 percent of full area respectively. Again the agreement is quite good for the 30 percent orifice, but not so good for the 10% orifice after the initial peak.
WSES-FSAR-UNIT-3 C-27 Thrust forces at the broken pipe break plane have been calculated by the modified RELAP code for typical GE, B&W, W, and CE steam and feedwater lines.
Typical results are shown in Figures 4-9 through 4-19.
Figures 4-9 through 4-12 illustrate the sensitivity of the calculated thrust force to the number of volumes and junctions chosen for the calculation. The long term blowdown thrust (quasi-steady) is essentially insensitive to the number of nodes chosen. The initial thrust, however, is quite sensitive, and several runs may be indicated prior to establishing a suitable minimum number of nodes for the particular vessel-piping configuration being studied.
Figures 4-9 and 4-13 among others, show also the extreme sensitivity of the predicted initial thrust to the choice of the j for the break control volume is chosen to the nearest elbow or change in area.
The RELAP predicted quasi-steady thrust compare very well with the theoretical values of Moody 25 / with the proper friction factor.
WSES-FSAR-UNIT-3 C-28 Revision 11-A (02/02)
NOMENCLATURE FOR PLAST General Notes: Lower case letters refer to variable referred to local coordinate axes, except as noted. Upper case letters refer to variables referred to global coordinate aces, except as noted.
{V} or V
~ designates the vector "V"
¨(DRN 00-1032)
{M} or M
~ designates the matrix "M" (DRN 00-1032)
[M]T = Transpose of M
{V}T = Transpose of V fx,fy,fz = forces along x,y,z and X,Y,Z axes respectively.
FX,FY,FZ M
M M
x y
z
= moments about s,y,z and X,Y,Z axes respectively.
MX,MY,MZ x,y,z = Displacements along x,y,z and X,Y,Z, axes X,Y,Z x,y,z = rotations about s,y,z and X,Y,Z, axes.
X,Y,Z
[k], [K] = stiffness matrices
{x}, {X} = displacement vectors
[T] = orthogonal transformation matrix w = weight of pipe/unit length L = length of pipe Ixx, Iyy, Izz = Area moments of inertia of pipe about respective local axes.
g = acceleration of gravity A = cross-sectional area of pipe steel m I I
I xx yy zz
= translational and rotational inertia terms along local principal pipe axes.
WSES-FSAR-UNIT-3 C-29 Revision 11-A (02/02)
I I
I I
I I
xy xz yz yx zx zy
, = product of inertia terms in rotational inertia tensor.
=
critical damping
=
circular frequency
[m],[M]=
mass matricas
[C]
=
damping matrix
'ij
=
i,j component of stress tensor J2
=
Second invariant of stress deviator ij
=
i,j component of strain tensor E
=
Young's modulus G
=
Shear modulus
=
Poisson's ratio s
=
Stress and plastic strain dependent coefficient for yield surface gradient in stress space.
f
=
Force dependent and plastic displacement dependent coefficient for yield surface gradient in force space.
¨(DRN 00-1032) f(ij)
=
yield surface function in stress space.
=
yield surface function in force space.
M M
M F
x y
z x
= Internal moments and axial force on pipe element node referred to local axes.
M M
M F
x y
z x
= Internal fully plastic moments and axial force for pioe element referred to local axes.
mi (i= x,y,z),fx = generalized moments and axial force referred to local element axes.
=
fluid density v
=
local fluid velocity vector
WSES-FSAR-UNIT-3 C-30 Nomenclature for RELAP e
=
u + v2 2
+ = total specific energy q
=
heat flux p
=
pressure
=
viscous stress tensor u
=
thermodynamic internal energy
=
potential energy function g
=
gravitational acceleration V
=
volume of control volume S
=
surface of control volume Mi
=
total mass in a control volume Wij
=
mass flow rate to control volume i through flow path h
=
specific enthalpy U
=
thermodynamic energy hij
=
enthalpy constant of fluid moving through flow path Qi
=
rate of energy transferred through all the surfaces of control volume i f
=
Fanning friction factor for single pipe Aw
=
wetted wall area A
=
cross-sectional gc
=
gravitational conversion constant n
=
normal unit vector
=
body force/unit mass
=
length of control volume
WSES-FSAR-UNIT-3 C-31 REFERENCES (1)
"Theory of Matrix Structural Analysis" - J.S. Prezemieniecki; Chapter 6, McGraw Hill, 1968.
(2)
"Survey of Advanced Structural Design Techniques," Z. Zudan; Nuc. Eng. and Design, 1969; pp. 400-440.
(3)
"Dynamics of Structures," W.C. Hurty, M.F. Rubenstein; Chap. 7, Prentice Hall, 1964 (4)
"A General Procedure for the Analysis of Elastic and Plastic Frameworks,"
G.A. Morris, S.J. Fenves; Office of Naval Research, Department of the Navy Contract None. 1834 (03), Project NR-064-183-(Univ. of Illinois Struct. Res.
Series 325.)
(5)
"The Mathematic Theory of Plasticity," Hill, R. Oxford University Press, 1950.
(6)
"Plasticity; Theory and Application," - Mendelson. A MacMillan, 1968.
(7)
"A Method of Computation for Structural Dynamics," N.M. Newmark, J. Eng. Mech. Proc. ASCE July 1959; EM3 pp. 67-94.
(8)
Users Manual-PLAST -2267 - Ebasco Publication.
(9)
"Vibration Theory and Applications," W.T. Thompson, pp. 294-299 Prentice Hall Inc, 1965.
(10)
WANGS "ELASTIC FRAME ANALYSIS PROGRAM", (George Washington University)
(11)
"ANSYS" by Swanson Analysis Systems, Inc. (SASI).
(12)
"FAP" (Frame Analysis Program); Department of Civil Mechanical and Environmental Engineering, George Washington University.
(13)
"Approximate Solutions for Impulsively Loaded Elastic-Plastic Beams," J.B.
Martin, LSS Lee, J. Appl. Mech. (Paper #68 WA/APM-23).
(14)
"RELAP-3... A Computer Program for Reactor Blowdown Analysis," W.H. Retting, IN-1321, 1970.
(15)
"Nodal Sensitivity in Modeling A Large Pressurized Water Reactor System for a LOCA," - C.E. Slater, R.D. Hentzen; CONF-710302, Vol. I USAEC Division of Technical Information, March 1971, p. 391.
(16)
"Effects of State Properties on High Pressure Compressible Hydrodynamics,"
D.V. Moore, G.E. Gruen; CONF-710302, Vol. 1, USAEC Division of Tech.
Information, March 1971, p. 426.
WSES-FSAR-UNIT-3 C-32 REFERENCES (Contd)
(17)
Transport Phenomena, R.B. Bird, W.W. Stewart, E.N. Lightfoor; John Wiley and Sons, Inc., New York, 1960.
(18)
The Dynamics and Thermodynamics of Compressible Fluid Flow, A.H. Shapiro; Vol.
1, Ronald Press, New York, 1953.
(19)
Supersonic Flow and Shock Waves, R. Courart and K.O. Friedricks, Interscience Publishers Inc., New York, 1967.
(20)
"Momentum Flux Terms in Transient Hydraulic Codes," K.V. Moore, C.E. Slater, L.J. Ybarrondo, G.E. Gruen; CONF-730304, USAEC Div. of Tech. Information, March 1973, p. 522.
(21)
"RELAP-4, A Computer Program for Transient Thermal-Hydraulic Analysis," ANCR-1127, USAEC 1973, D.V. Moore, W.H. Rettig.
(22)
"Maximum Flow Rate of a Single Component, Two Phase Mixture," F.J. Moody, Journal of Heat Transfer, ASME Trans 87. 134-142, Feb. 1965.
(23)
"Blowdown Rates of Initially Saturated Water," U. Simon; CONF-730304, USAEC Div. of Tech. Information, March 1973, P. 172.
(24)
USAEC National Reactor Testing Station Quarterly Technical Reports, IDO-17266 (Feb. 1969).
(25)
"Time Dependent Pipe Forces Caused by Blowdown and Flow Storage," F.J. Moody, Journal of Heat Transfer, 73-FE-23.
(26)
"Subcooled-Blowdown Forces on Reactor-System Components; Calculational Method and Experimental Confirmation," G.H. Hanson; TID-4500, USAEC, June 1970.