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-1APPENDIX 3.6BPIPE WHIP ANALYSIS(ETR 1002, Appendix C)
WSES-FSAR-UNIT-3 C-2APPENDIX C of ETR 1002PIPE WHIP ANALYSIS1.METHODS OF ANALYSIS1.1Dynamic AnalysisThe 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 oftime 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 dependentthrust 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.2Equivalent Static AnalysisBecause of the relatively large number of breaks which must beanalyzed, 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 initialrestraint system design may be based on an equivalent static analysis.Following such static analysis, the designed restraining system will beverified 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-3C-3Revision 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)
F es = K 1 K 2 PA(1.1)(DRN 00-1032)
Where: A is the break area
K 2 is a thrust coefficient which has the value of 1.26 for steam andsaturated water, and 2.0 for subcooled water, and K 1 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 K 1 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 K
- 1. Our experience to date, although limited, points a value of 3.0 as sufficient for Ebasco's typical piping restraintconfigurations wing gaps below 4 inches.It should be stressed, however, that the choice of K 1 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 K 1 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 operatingpressure in the line.
For feedwater lines, however, the subcooleddecompression 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 K 1 2.OP o A at time zero, decreasing linearly to K 11.26P sat 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 P sat the saturation pressure corresponding to the maximum operating
temperature in the line.(DRN 00-1032)
WSES-FSAR-UNIT-3 C-42.DESCRIPTION OF PLAST2.1Analytic DescriptionPLAST is a FORTRAN program, based on the "matrix-displacement method"1/, for three dimensional space frames in which sections of pipecomprise 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 system's 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 ondamping. Upon defining a time dependent forcing function, an elasto-plastic constitutive law and constraints on the system's boundaries, a set of ordinary differential equations in the time domain, are delineated in a form amenable to solution. The details of the aboveoutline will now be presented.2.1.1The Element Stiffness MatrixA typical straight pipe element, Figure 2.1.1 is shown oriented to aset 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 positivedisplacements, 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:
{}xxyzxyzxyzxyz = [() I () ]
T J{}))]f = [(f I (f TxxJffMMMffMMMyzxyzyzxyzwhere lower case subscripts refer to local element axes along which vectorcomponents 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) isthen found by substituting the first two of equations (2.1.1.2) into (2.1.1.1) and premultiplying both sides by [T]
T. Thensince[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.2The Element Mass MatrixThe element mass matrix referred to local axes is a 12 X 12 diagonalmatrix 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.3Damping MatrixPLAST utilizes a damping matrix that is a linear combination of thestiffness 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.4Elasto Plastic Constitutive LawsConstitutive laws governing the elasto-plastic behavior of threedimensional 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 andyielding 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 directioncosines with all principal stress axes. 5,6/ The surface is defined mathematically as;(DRN 00-1032)
F('ij)-e 2/3 = 0(2.1.4.1)
The yield condition then states;
If F('ij)-e 2/3>0and F('ij) > 0 then yielding occurs If F('ij)-e 2/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
= (fi P)(2.1.4.4) fi P (i = x, y, z) = generalized moments and forces.
This is the equation of a unit sphere in a generalized force space ofgeneralized forces. The increment in generalized plastic displacements
of a beam section undergoing yielding is a vector normal to and
pointing outward from the yield su rface 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 hardening5/ 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 plastichinges 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 onedimensional 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 2p x p e(2.1.5.2)(DRN 00-1032) as the effective stress, and the increment in effective plas tic 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 inte rval t n , 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 MethodThe equations of motion for the system of lumped masses tied to eachother 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
displacementsderives from the observation that the force displacement relationship;[K] {X} = {F}
applies only to elastic media. By artificially subdividing,displacementsinto "elastic" and "plastic" component s, a plastic"correction" force arises; ie,[K] {X - X p} = {F}
and(DRN 00-1032)[K] {X} = {F} + [K] {X p} + {F} {F p}(2.1.6.2)
Hence equation (2.1.6.1) may be rewritten as;
[M] {X} + [C] {X} + [K] {X}
e = {F} + {F p}(2.1.6.3)where {F} , {F p} are the externally applied force vector and theplastic 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 Newmark7/
Appendix C. An internal procedure of checking convergence and stabilityis 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.(DRN 00-1032)
Boundary ConditionsIn 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 anoverspecified 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 moredisplacement 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 t
- n. Then if 0.~_~nor(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 n f
&~where , N designates linear and nonlinear functions respectively of their arguments; and the dot over the " n~ " argument designates "the increment of n ".(DRN 00-1032)
WSES-FSAR-UNIT-3C-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 displacemen t 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/(DRN 00-1032) 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 significanteffect 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 ofmagnitude 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 theserotations must be of the order of 2(.001)1/2 radians or
+/-3.62 degrees, before second order terms begin to become significant.(DRN 00-1032)
The pipe system is monitored by PLAST to determine if rigid bodyrotations 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 wouldoccur with respect to monitonically increasing deformations.
WSES-FSAR-UNIT-3C-11 Revision 11-A (02/02)2.2Discussion and Summary of Program Features and CharacteristicsPLAST solves initial value problems in the stress analysis of piping systems thatare 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 pipingsystem 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 andis 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 no nincreasing 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 exceededv/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 reviewingthe 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-3C-12 Revision 11-A (02/02)element. The updated maximum parameters mentioned in the preceding paragraphare 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.3Typical Pipe Whip Dynamic Analysis 2.3.1 Steam Line with Yielding RestraintsThe 43 in. steam line shown in Figure (2.3-1) was analyzed for a guillotine break atnode 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 therestraint 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
pastblowdown initiation. The restraint is seen to permanently deform to about 8.5in. 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 occursin 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 #GapElastic Spring "K"(DRN 00-1032)26."8.536 X 106 #/in.46.""
56.""
85.5""
95.5""
114.""
WSES-FSAR-UNIT-3C-13 Revision 11-A (02/02)3.1PLAST 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:2222 x y x y t)+ m )= p (x,t)
(EI 22 (3.1.1)and its solution if found to be (see Reference 9/ for example) yxt(,)=(x) q(t) nn=1 n(3.1.1.2)where q n , n are the mode displacements and mode shapes respectively and are shown below as; q n n=2 L (+1-(-1)(t-n) Mo nsin 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)(DRN 00-1032)
WSES-FSAR-UNIT-3C-14 Revision 11-A (02/02)
The solution to this problem is displayed in Figures (3.1.2 A&B) as displacementand 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.2E lastic 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 ofapplication whereas Reference10/
does not; hence increasingdeviations 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 to10, 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 10lumped masses. However, in the case the masses were arbitrarily
multiplied by (25) 2 which results in a minimum system period 25 timesgreater 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-3C-15 Revision 11-A (02/02)(DRN 00-1032)
Comparisons are shown in Figures (3.1.5A&B) for the same problem runthrough ANSYS 11/ using zero strain hardening and FAP 12
/ using 10percent 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.(DRN 00-1032)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. Thedeviation 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 addi tional 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 animposed velocity Vo at time t o = 0 uniformly distributed along its length and the second is a half sine curve distribution of an imposed velocity U* at time t o = 0. In both cases, the imposed velocity drops to zero for all times t > t o 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:
U i 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 U i and U*i areinitially 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-16V* = 4 V oNon-dimensional plots of displacement and velocity responses of the beam at itscenter 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 M o = fully plastic moment in beamm = beam mass/unit lengthE,I = Young's 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 ascalculated 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 tocontinue to oscillate past this point. The authors' contention that solutiondifferences decrease with time or are non-increasing is clearly demonstrated.Figure (3.1.6B) compares central deflection of the above beam versus time forreference 13 / and PLAST, subject to the sinusoidal initial velocity distribution andelastic-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 RELAP4.1GeneralRELAP-3 14 / is a FORTRAN program which describes the behavior of hydraulicsystems under transient conditions. The program itself, which has been welldocumented elsewhere 14,15,16 /, was developed specifically to describe thebehavior 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 controlvolumes connected by flow paths. Check valves, valves, pumps can be inserted at any flow path.
WSES-FSAR-UNIT-3C-17 Revision 11-A (02/02)
To simulate a break, a leak can be opened either instantaneously, or as a functionof 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 equationa - 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)(DRN 00-1032)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-3C-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-=vi V i V d dt d (4.5)But i i V dv$ is just the total mass in volume i = M i and!!"$dS V i S 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 andpotential 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 i V)q+v h (-=dV)q+v (V-=dV d t u)d(i i S V(4.8)i ij j ij Q+h W=dt du i (4.9)The momentum equation restricted to one dimensional flow and integratedover 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)(DRN 00-1032)
WSES-FSAR-UNIT-3C-19 Revision 11-A (02/02)
For the friction term, the pressure losses are assumed to be proportional to v 2,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 A w/ 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 dW i i l , 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 A j j k 2 c g 144 k'=(4.15)P i and P i+1 are the inlet and outlet pressure of a control volume in psi, and Z i+1 and Z i , are the elevation of two contiguous controlvolumes 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 iandi + i, is added to the momentum equation(DRN 00-1032)
WSES-FSAR-UNIT-3C-20 Revision 11-A (02/02)(DRN 00-1032) j pj 1+i i j 1 144 1+P+P-P=144 1j 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)(DRN 00-1032)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.2Choice 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)(DRN 00-1032)
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
volumet W j M j(4.19)
WSES-FSAR-UNIT-3C-21 Revision 11-A (02/02) andt h j U (W j 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 2x (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 2P 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)(DRN 00-1032)
WSES-FSAR-UNIT-3C-22 Revision 11-A (02/02)
Thus(DRN 00-1032) 0=x v+x v+t P 1 2C (4.28)0=x 1+x v v+t v(4.29)Ift P<<x v (4.30)andt v<<x v v (4.31)then 0=x v+t P 1 2C (4.32)(DRN 00-1032)v t=1P x=0 (4.33)Thetwo equations above can be reduced to the familiar wave equation2 22 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 obtainedfrom 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)(DRN 00-1032)
WSES-FSAR-UNIT-3C-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)(DRN 00-1032)
Where S being the total surface of the fluid element, can be divided into the wetted surface, S w and the surfaces across which flow occurs, S 1 and S 2.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 F ext 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 thewalls 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!!!!!$$$"$$(DRN 00-1032)
(4.38)The last term in equation (4.38) is small in blowdown studies and can thus beneglected. 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 S 1$ 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-3C-24 Revision 11-A (02/02)appear at the exit plane. If there is a discontinuity at the exit plane, as is the caseforchocked 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 (-d t d-=i i S V ds V d F(4.39)4 4 43 4 4 42 1 chocked is flow if only Pnds s Pnds s 2 1$$(DRN 00-1032)
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 V i 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 T m=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 gt+t cl (4.42)(DRN 01-1032)
If one assumes that a break develops instantaneo usly, 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.(DRN 01-1032)
WSES-FSAR-UNIT-3C-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)(DRN 00-1032)Herein P i stands for the pressure inside the control volume (fluid pressure) and P i+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-P a)A, does then indeed equate the acceleration component T a.Unbalanced external forces at the break area occur at times other than zero, whenthe 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 P e 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 Ag W t jc jt=W+ 144(Pe- Pa) A+(W g j j jjt+t i c 2)l (4.45)The RELAP program 14/ has been modified to solve for the three components andthe total thrust force. Three additional minor edit variables have been added to theprintout 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 CalculationsThe RELAP program chooses as flow through a junction, the smaller of the inertial flow or choked flow.
WSES-FSAR-UNIT-3 C-26The limiting mass flow is determined by use of Moody's two phase choked flowmodel 22 /. This model also allows computation of the throat pressure (taken asthe exit plane pressure) as a function of stagnation pressure and enthalpy.Simon 23 / has discussed the many models which may be used to predictblowdown flow rates of two phase fluids, such as would result following an initialsubcooled blowdown. His comparisons among limited test data show that for sufficiently large diameter pipes, and for breaks located a sufficient distance from the reservoir, Moody's 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 initiallyfilled 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 RELAPwith 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.5Typical ResultsThere are unfortunately no large scale blowdown thrust values reported in theliterature. 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 tothe thrust measured by Hansen 26 / for the apparatus shown in Figure 4-6. Theagreement is rather good.Figure 4-7 and 4-8 also compare calculated thrusts with Hansen's measurement.In these cases an orifice was present in the discharge pipe of the apparatusreducing 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-27Thrust forces at the broken pipe break plane have been calculated by themodified 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 tothe number of volumes and junctions chosen for the calculation. The long term blowdown thrust (quasi-steady) is essentially insensitive to the number of nodeschosen. The initial thrust, however, is quite sensitive, and several runs may beindicated 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 thepredicted initial thrust to the choice of the j for the break control volume ischosen to the nearest elbow or change in area.The RELAP predicted quasi-steady thrust compare very well with the theoreticalvalues of Moody 25/ with the proper friction factor.
WSES-FSAR-UNIT-3C-28 Revision 11-A (02/02)
NOMENCLATURE FOR PLASTGeneral 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 f x ,f y ,f z = forces along x,y,z and X,Y,Z axes respectively.
F X ,F Y ,F ZMMMxyz ,,= moments about s,y,z and X,Y,Z axes respectively.
M X ,M Y ,M Zx ,y ,z = Displacements along x,y,z and X,Y,Z, axesX ,Y ,Zx ,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
I xx , I yy , I zz = Area moments of inertia of pipe about respective local axes.
g = acceleration of gravity
A = cross-sectional area of pipe steelmIIIxxyyzz,,, = translational and rotational inertia terms along local principal pipe axes.
WSES-FSAR-UNIT-3C-29 Revision 11-A (02/02)IIIIIIxyxzyzyxzxzy
,,,,,, = 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 J 2=Second invariant of stress deviatorij=i,j component of strain tensorE=Young's modulus G=Shear modulus
&=Poisson's ratios=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.(DRN 00-1032)
=yield surface function in force space.MMMFxyzx ,,, = Internal moments and axial force on pipe element node referred to local axes.MMMFxyzx,,, = Internal fully plastic moments and axial force for pioe element referred to local axes.
m i (i= x,y,z),f x = generalized moments and axial force referred to local element axes.=fluid density v!=local fluid velocity vector WSES-FSAR-UNIT-3 C-30Nomenclature for RELAPe=u + v 2 2 + = total specific energy q=heat flux p=pressure=viscous stress tensoru=thermodynamic internal energy=potential energy functiong=gravitational accelerationV=volume of control volume S=surface of control volume M i=total mass in a control volume W ij=mass flow rate to control volume i through flow pathh=specific enthalpyU=thermodynamic energy h ij=enthalpy constant of fluid moving through flow path Q i=rate of energy transferred through all the surfaces of control volume i f=Fanning friction factor for single pipe A w=wetted wall areaA=cross-sectional
g c=gravitational conversion constant n=normal unit vector=body force/unit mass=length of control volume WSES-FSAR-UNIT-3 C-31REFERENCES(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. andDesign, 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)User's Manual-PLAST -2267 - Ebasco Publication.
(9)"Vibration Theory and Applications," W.T. Thompson, pp. 294-299 Prentice HallInc, 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 andEnvironmental 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 aLOCA," - 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-32REFERENCES (Cont'd)(17)Transport Phenomena, R.B. Bird, W.W. Stewart, E.N. Lightfoor; John Wiley andSons, 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, IntersciencePublishers 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, USAECDiv. 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 Methodand Experimental Confirmation," G.H. Hanson; TID-4500, USAEC, June 1970.