- Email: [email protected]
Inverse and direct transfer functions for the fatigue follow-up of piping systems submitted to stratification
ELSEVIER
Nuclear Engineering and Design 153 (1995) 183 195
Nuclear Engineering and Design
Inverse and direct transfer functions for the fatigue follow-up of piping systems submitted to stratification M. Guyette, M. De Smet Tractebel Energy Engineering, Avenue Ariane 7, 1200 Brussels, Belgium
Abstract In this paper we outline a methodology to assess the fatigue induced in piping systems submitted to thermal stratification. More specifically, the transformation from the measured outer wall temperature time histories to stress time histories in any point of the line is treated. By means of inverse transfer functions, the fluid temperature distribution is calculated from the outside wall temperatures measured in a limited number of temperature sections. Using direct transfer functions, the local stresses due to stratification may be determined as well as the pipe free curvatures and the pipe free axial strains. Using a finite beam element model of the line, the global response of the line (in terms of displacements or stresses) due to the applied curvatures, axial strains, end point displacements, internal pressure and possible contacts with the pipe environment may be determined. The method is illustrated for the surge lines of the Doel 2 and Doel 4 nuclear power plants. An excellent correlation is found between measured and calculated displacements. Typical stress time histories are shown for a plant cool down.
1. Introduction Thermal stratification has been observed in several pressurized water reactor systems for a couple of years. Thermal stratification typically occurs in the surge line and the main feed water lines but other lines such as the pressurizer spray line may also exhibit the phenomenon. The large number of thermal cycles of great amplitude produced by stratification has raised some concern about the damage induced by fatigue in these lines. The piping systems affected by stratification have been equipped with thermocouples mounted on the outer pipe wall and with displacement
sensors. Temperature and displacement time histories have been recorded during several heat-ups and cool-downs. In view of a fatigue analysis one has to pass from the recorded outer wall temperature time histories to the stress time histories at the most highly loaded points in the piping. The procedure which has been used is outlined in the following sections.
2. Main problems to be solved The physical parameters which influence the fatigue of a highly loaded point are usually very
0029-5493/95/$9.50 © 1995 Elsevier Science S.A. All rights reserved SSDI 0029-5493(94) 00826-K
184
M. Guyette, M. De Smet / Nuch, ar Engineering and Design 153 (1995) 183 195
limited in number. In most cases, the knowledge of the local pressure and of a few values of fluid temperature (usually less than three) is sufficient to allow the evaluation of the local stresses by use of direct transfer functions. The problem of investigating the thermal stratification in a pipe is notably more complex since, as a result of the non-linearity of the temperature distribution over a cross-section, even in an unrestrained pipe, stratification induces local stresses and a pipe curvature. The stresses at a given location in a restrained pipe submitted to stratified flow are therefore a function not only of the local temperature distribution but also of the curvatures induced by the stratification at all other locations. Usually the stratification pattern is variable along the length of the line. The local temperature distribution as well as the distribution of the curvature along the line (or the linear thermal gradient, which is proportional to the curvature) must therefore be gained from the monitoring sensors. To reach a sufficient accuracy, a minimum of five measurement points per section are required and a minimum of four or five measurement sections along the length of the pipe has to be contemplated. Ideally, the measurement should give the fluid temperature which is difficult, the measurements of the pipe outer wall temperature being easier by far to perform. This leads either to the assumption that the external wall temperatures are a fair representation of the water temperatures (which appear to be inaccurate and non-conservative particularly for transients) or to the development of an algorithm allowing the calculation of the fluid temperature from the external wall values. The development of this algorithm requires the knowledge of the film coefficient between the fluid and the wall, the value of which is usually not known with good accuracy. Large inaccuracies on the fluid temperatures could therefore be generated by such a process. In the sequel, one is, however, not directly interested by the fluid temperature but rather by the temperature distribution in the metal. Provided that use is made, in the direct transfer functions (allowing calculation of the temperatures, stresses and strains in the
metal from the fluid temperatures), of the same film coefficient values, the inaccuracies mentioned above are fully compensated. The fluid temperature must be regarded as an intermediate result without direct interest in the calculation chain. This point will be commented on in more detail in Section 4. From a physical point of view, this deduction may appear as relatively evident since the whole process consists of starting from temperatures measured in some particular points in the metal to deduce other temperatures in the same piece of metal. Once the local curvatures are known, they must be introduced in a global model of the line which will be used to determine the forces and moments at all locations of interest. The problem is complicated by the fact that the gaps at the pipe whip restraints are not always sufficient to accommodate the large displacements caused by the stratification. Contact forces with these restraints take place at some locations and completely disturb the distribution of moments. The main problems to be solved in view of a fatigue analysis based on measured outer wall temperature time histories are therefore --- the development of a practical way of coping with the complex temperature and stress patterns both in a cross-section and axially, - - the calculation of the fluid temperature distribution from the measured outside wall temperature, - - the evaluation of the global behaviour of the line (influenced by the local conditions in each cross-section), and the presence of gap contact situations at some pipe whip restraints. The way by which these problems have been solved is described in detail in the next sections.
3. C o m p l e x i t y o f the temperature and stress fields
The temperature field in the pipe presents important variations in azimuthal and radial directions while the variations are much smoother in the axial direction. In addition, the temperature field is strongly variable in time. Based on these assumptions and neglecting the small variations in
M. Guyette, M. De Smet /Nuclear Engineering and Design 153 (1995) 183-195
axial direction, the temperature field will satisfy the differential equation a2T 1 a T 1 a2T ~ar + -r ~ r + r 2 ~02
1 aT K at
(1)
where T is the metal temperature (see also nomenclature in Appendix A) and K is the metal thermal diffusivity. The following boundary conditions will be applied:
185
eigenmodes and eigenfrequencies in structural dynamics will be used for each harmonic. The time history response of the temperature at a given time t for the point of coordinates r and 0 in the pipe will be
Tk(r't)=[M°k(r)--~M'k(r)lTfk(t),=l 4-
Y, K4;,,M~,(r) /=l
2(ar~ --
k~rJ
= htT(r,, O, t) - Tf(O, t)]
(2)
....
1
x
-_o r=r e
which assume a heat exchange, with a film coefficient h, between the fluid and the metal at the inner surface and a perfect insulation at the outer surface. The initial temperature is assumed to be uniformly equal to 0. This last assumption could be considered as an oversimplification of the problem. Taking, however, into account that the transients that will be followed are extremely long compared with the period needed for the initial conditions to die out, this assumption has a negligible effect on the computed results. In each cross-section, the temperature presents a mirror symmetry with respect to the vertical plane. Based on these observations, it has been decided to express the temperature field in the pipe as a Fourier series of the azimuthal angle around the pipe axis. The temperature will thus be expressed as
T(r, O, t) = ~ T,(r, t) cos kO
(4)
k=0
In order to keep the equations linear, the thermal as well as the mechanical properties of the pipe material such as the thermal conductivity, the thermal capacity, the Young's modulus and the thermal expansion coefficient will be considered to be independent of the temperature. The thermal conduction equations in transient regime will be solved for each harmonic of the Fourier series. A technique of eigenvalues and eigenfunctions (see Boley and Weiner (1985) or Carslaw and Jaeger (1959)) similar to that of the
e x p [ - K~b~k(t - v)lTfk(r) dr
(5)
where Mik(r) is the eigenform of index l for the harmonic k, ~bCkthe eigenvalue of index l for the harmonic k, ~ the thermal diffusivity of the pipe material, Tfk(t) the component of harmonic k of the fluid temperature and Mok(r) the steady state temperature distribution of harmonic k for a fluid temperature equal to unity in the same harmonic. The eigenforms M~k(r) are normalized in such a way that
mok(r) = ~ mkl(r)
(6)
/=1
The eigenforms are combinations of first- and second-kind Bessel functions of the same order as the harmonic to which they belong. The derivation of the stresses caused by the temperature field is performed by bringing the value of this temperature field into the thermoelastic equations. Use is again made of the mirror symmetry with respect to the vertical plane. In a first step, a strict plane strain hypothesis is used for the determination of the stress field. In a second step, this plane strain hypothesis is transformed into a generalized plane strain hypothesis, under the assumption that the net axial force and the bending moment perpendicular to the vertical symmetry plane are both equal to zero. It is assumed that the thermoelastic equations are not explicitly functions of the time, which corresponds physically to neglect of the inertia forces caused by the thermal expansion. Because of the slowness of the temperature changes on the one hand and the smallness of the thermal dis-
M. Guyette, M. De Smet / Nuclear Engineering and Design 153 (1995) 183-195
186
placements on the other hand, these inertia forces are indeed negligible. The solution technique is based on the theory of the thermoelastic potential (see Timoshenko and Goodier (1970)) and on a Fourier development of all the functions (displacements, stresses, strains) intervening in the equations. Beyond the generalized plane strain conditions mentioned above, the boundary conditions imposed on the stress solution are simply that the radial normal stress and the shear stress a,.o vanish on both the inner and outer surfaces of the pipe. For each Mik(r ) eigenform (l being the mode subscript and k the harmonic subscript), the thermoelastic theory will give stress functions Stk(r) for each stress component which, after combination, gives the stress response of harmonic k:
4. Calculation of the azimuthal fluid temperature distribution from the outer wall distribution: inverse transfer functions
The equations above may only be used if the fluid temperature is known. In fact as already explained, the outside wall temperatures are, in most cases, the only ones available. A formulation has therefore to be developed to obtain the fluid temperature from the outside wall temperatures. In such a process, the transfer functions used may be called "inverse" since they allow us to obtain the cause of the phenomenon from the knowledge of its effects in opposition to normal or "direct" transfer functions which permit us to calculate the effects of a phenomenon from its causes. Each of Eqs. (5) relative to a particular harmonic may be written for the external radius r e and put in the form
Tk(re, t) = Ck(re)Tfk(t ) +
Gk(r e, t -- r)Tfk(r) dr
+ ~ ~'O~kS,k(r)
(8)
I-I
x
;o
exp[ - ~CC~k(t -- r)]Tfk(r) dr
(7)
where S~k(r) is the modal stress caused by the application of temperature distribution M~k(r) and Soa.(r) is the steady state stress caused by the temperature field Mok(r ). The last equation is valid for any one of the stress components. A similar formulation may also be used for the displacements, strains etc. The stresses are first calculated under the assumption of plane strain conditions. Under these conditions the harmonic 0 usually causes an axial force. In a second step, to obtain the pipe free expansion, an axial force of opposite sign is superimposed. The same reasoning may be used concerning the bending moment caused by the harmonic I. The superposition of a moment of opposite sign causes a curvature of the pipe called hereafter the free curvature. The above formulation allows, for any section where the fluid temperature is known, the local stresses, the free axial expansion and the free curvature to be calculated.
with Ck(re) = Mok(rc)- ~ Mlk(re)
(9)
I-I
Gk(re, t) = ~ ~Cb21Mlk(re) exp(-~c~b2/t)
(10)
l-I
r e being a constant, Ck a constant and Gk a function of t only. Eq. (8) is the basis to obtain the inverse transfer functions. In the case where Tk(re, t) is known and Tfk(t) unknown, it may be considered as an integral equation in Tfk, the solution of which may be written as
Tfk(t ) = Qk(re)Tk(re, t) +
f+
Uk(re, ti~)Tk(re, r) dr
(11)
The mathematical developments leading to the values of the constant Qk and to the expression of Uk are beyond the scope of this presentation. They are described in detail in Guyette (1994). The Uk function is the inverse transfer function. As this function decreases sharply as t tends to
M. Guyette, M. De Srnet / Nuclear Engineering and Design 153 (1995) 183-195
+ oo or - oo, the final expression of Tfk(t) may be put in the form Tfk(t ) ~ Ok(re)Tk(re, t) -Jr
f
t+A
Uk(Fe, t -- r)Tk(re, z) dr
(12)
,,)t-A-
where A - and A ÷ are chosen in such a way as to give a good engineering accuracy to the expression of Tfk. A function Uk and a constant Qk are obtained for each harmonic in turn. The inversion process is therefore executed in the following steps: (1) perform a Fourier analysis of the outer wall temperature to obtain the Fourier coefficients of the external wall temperature; (2) knowing the time history of each of these outer wall Fourier coefficients, obtain, with the help of Eq. (12), the response in time of the Fourier coefficient for the fluid temperature; (3) to obtain the fluid temperature at a given azimuthal position, assess the value of the Fourier series at this point. The inversion process is therefore fully two dimensional. The formulation which is followed warrants the uniqueness of the solution. Contrary to what is usually given in the literature on the inverse problem, the solution presented is based on the use of future times: the integration domain in Eq. (12) extends up to t + A +. This has the advantage that the solution is very stable and accurate (which is not the reputation of usual inverse problems) but it is necessary to retain measurements until the time t + A + before the temperatures and stresses at time t can be evaluated. A + being usually of the order of a few hundred seconds, this may not be considered as a true drawback for the type of problem which is solved here. From the theoretical point of view, the main approximation is the limitation of the number of eigenfunctions taken into account in the inversion process of each harmonic. It is, however, shown, in Guyette (1994), that, from a practical point of view, the accuracy on the results of the inversion process, i.e. the bulk fluid temperatures, is mainly limited by the fluctuations of the outer
187
wall temperatures caused by the measurement devices: beyond systematic errors, the measurement chains usually produce erratic fluctuations which are of the order of a few tenths of a kelvin (for temperatures between approximately 150 and 350 °C). These erratic fluctuations induce a numeric noise which is of the order of 0.5 K on the results of the inversion. An important point to be commented on is the influence of the film heat transfer coefficient on the results of the direct and inverse transfer functions. This film coefficient is taken into account in both types of functions and any uncertainty on its values may strongly affect the results. However, when use is concurrently made of inverse and direct transfer functions and provided that the same value of the film coefficient is used for the two types of functions, a compensation of the error occurs. For instance, if a very large film coefficient is used, the inverse transfer functions will predict the inside wall temperature in place of the bulk fluid temperature. However, if the same value of the film coeffÉcient is used in the direct transfer functions, the calculated variables, which may be temperatures, displacements, stresses or strains, will assume that the starting temperature is the inside wall temperature and the calculated values will be correct. As a matter of fact, the process described above, combining inverse and direct transfer functions, may be seen as a way of evaluation of the heat flux exchanged with the fluid at the inner wall, from which the temperature at any point in the metal may be calculated. This process does not attempt to make predictions of what occurs in the fluid itself. This conclusion is of course limited to the simulation of the phenomena which have a measurable impact on the outer wall temperatures. Rapid fluctuations of the fluid temperatures which are undetectable at the outer wall surface, such as what is known in the literature as the striping phenomenon, are not accounted for in the abovedescribed process. Their contribution to the usage factor of the equipment has been shown to be rather limited. The expression (12) is particularly useful when the outer wall temperatures are sampled at regular time intervals. If one makes the hypothesis either
188
M. Guyette, M. De Srnet / N u c l e a r Engineering and Design 153 (1995) 183-.-195
of a linear variation in the outer wall temperature (Co continuity) or a more elaborate C, continuity (using third-degree polynomials), the expression (12) may be transformed into T,-,(t) = ~, w,T~(t + i~)
(13)
i
where ~ is the sampling time interval, the index i being such as to cover the time period between t - A - and t + A +. The w~ are weights depending on the value of the constant Q~ and of the function U~. As the latter are unique for a given geometry and material, the w~may be calculated once for all, which allows a fast calculation of the fluid temperatures from the outer wall temperatures. The following example illustrates the good accuracy obtained in the use of the inverse transfer functions. Figs. 1 and 2 show the measured outside wall temperatures recorded during a typical stratification transient in a cross-section of the surge line of the Doel 2 power plant. The measured values indicated by diamond symbols on the figures correspond to five thermocouples located 45 ° apart on the outside wall of the piping. A fitting of the
220 200 180 ~,~:
"---222. '2~
0 16o = 140 ~ 120 100 80 60
i
0
i
i
150
300
i
I
i
i
450
600
750
900
Time (s) Fig. 1. Accuracy check of the inverse transfer functions: O, measured outside temperatures; the fine curves are the calculated inside temperatures (obtained by the use of inverse transfer functions) while the thick curves are the calculated outside wall temperatures (derived with direct transfer functions).
220
, •,,, .v,~ ~,~,"'~ <'" <"""" 2
200 180 L
~ / / ~ ' " ~"
~" o 160 140
K 120 ~,,
•
lOO 80 60 900
i
i
i
1050
1200
1350
i
i
1500
i
i
i
1650
1800
Time (s) Fig. 2. Accuracy check of the inverse transfer functions (continued). See comments from Fig. 1.
outside wall temperatures by a Fourier series is first performed for all measured times. With the help of inverse transfer functions, applied on each term of the Fourier series, the development of the inner wall temperature is obtained. Assessing the series for the five azimuthal positions of the thermocouples, one obtains the inside wall temperatures just in front of the thermocouples. These are shown as fine line curves on Figs. 1 and 2. Using now direct transfer functions, the outside wall temperatures may be recalculated. The values corresponding to the five thermocouple positions are plotted as thick line curves on Figs. 1 and 2. Any error or inaccuracy in the inner wall temperatures would bring a discrepancy between the starting values of the outer wall temperatures and the final calculated values. The figures show a very good agreement with a discrepancy well below 1 K (thick lines marked by diamonds) except at the onset of the sharp transient where differences are about 2 K. The figures also clearly show that the inner wall temperature curves anticipate with a sharper slope the trends seen for the outside wall temperature, which conforms to the physics of the phenomenon. Examination of Figs. 1 and 2 shows that the inside wall temperature may be as much as 40 K
M. Guyette, M. De Smet / Nuclear Engineering and Design 153 (1995) 183-195
higher than the outer wall at some moments of the transient. This illustrates the errors which could have been made if the outer wall temperature had been used in place of the inside wall temperature.
189
I Pressu- I I Section 1 Section 2
Section3
Section 4
Section5
I influence influence I function function
influence function
influence function
influence function
.'--'--<.T_-:::-:----: 0
0
0
:-:..:
=t'3 0
~_1
5. Global behaviour of the line
Fig. 3 shows the layout of the surge line of the Doel 2 power plant with the location of the provisional measurement devices (thermocouples and displacement sensors) used to monitor some heat-ups and cool-downs. This provisional measuring installation is typical (number and position of the captors) of what may be used for permanent monitoring. The use of the inverse transfer functions combined with the stress and strain analysis described here above allows us to calculate the local stresses, the axial free expansion and free curvature of the pipe in each measurement section. To assess the global behaviour of the line, a model of the latter is built; this model may be made with the help of either a specific piping
T1.2 TI.1
T2.2~.
TI.S \
Section3 ~ s~nsor3
I
T4.2[ ~
4
Fig. 3. Layout of the Doel 2 surge line showing the locations of the temperature and displacement sensors.
Fig. 4. Schematic representation of the influence functions of the various sections.
calculation program provided that it is able to assess the effects of transverse thermal gradients (causing a curvature of the pipe elements) or a general-purpose finite element program. In the second case, care must be taken that the elements representing the pipe elbows and other specific piping fittings reflect the true flexibility. It is assumed that the axial free expansion and the free curvature vary linearly between measurement sections. Influence functions are defined for each section as shown in Fig. 4. The function relative to a particular section varies from 1 in the considered section to 0 in the adjacent sections. A series of responses of interest are defined for which the time evolution during the transients is wanted. These responses may be displacements, stresses, strains or forces at particular locations in the piping system. These responses are assumed to be stored according to a well-defined order in a vector v. For each measurement section, load cases corresponding to unit curvature and to unit average temperature, distributed in the line according to the above-mentioned influence functions, are calculated. The values of the responses of interest are extracted from the results of these load cases and stored column by column in the matrix V~ for what concerns the unit curvatures and in V¢ for the mean unit temperatures. Load cases corresponding to unit end displacement in each of the coordinate directions are also calculated and the responses of interest are stored columnwise in the V s matrix. As, from the previous calculations, the mean temperatures and the curvatures are known for each time step in the transient, the values of the
M. Guyette, M. De Smet / Nuclear Engineering and Design 153 (1995) 183 195
190
responses of interest may simply be calculated for each step by
Meas. - dp. 2
(14)
v = Vp[~ + V i . T + V a 6 + Vpp
where [i and I" are vectors containing respectively the curvature and the mean temperature in the measurement sections, 6 is the vector containing the end displacements for the considered time, Vp is the vector of the responses of interest for a unit internal pressure and p is the internal pressure. The accuracy of this global model may easily be checked by performing a comparison between the measured values of the displacements and the calculated response of the model. As shown in Fig. 3, in addition to the five measurement sections equipped with thermocouples, the vertical displacement of the line is measured at three locations. These displacements are compared with the calculated displacements in Figs. 5 and 6 for a typical transient the recorded temperatures of which are given in Fig. 7 for what concerns Section 1 and in Fig. 8 for Section 5 (the temperature variations in Sections 2 - 4 are similar to those of Section 1 and are not shown). Examination of Figs. 5 and 6 shows a very good agreement be-
Calc. - - - - - Meas. - . . . . . . . Calc. dp. 2 dp. 3 dp, 3
7O
60
50 4o
.~ 30
ilo 2°
-10 0
'
,
,
,
3600
7200
10800
14400
Calc. dp. 1
,
Time (s) Fig. 6. C o m p a r i s o n of measured and calculated vertical displacements at sensors 2 and 3 (Doel 2).
-
-
T1.1
- -
-
T1,2
.........
Tt.3
.
.
.
.
T1.4
-
-
T1.5
220
200
/
,s0 Meas. dp. 1
,
18000 21600
I j,/
i140
70
12o
6O
(~ 100
f/;:? t,/' i~ii :
i,,
J.: i ! f
,
t
,;
ii ,'
Z" 50 60
~ 40
,
64800
68400
,
,
,'.-r.-,'
,
,
7 2 0 0 0 7 5 6 0 0 79200 82800
,
86400
r~
Time (s) "~ 30 Fig. 7. Evolution in time of the measured outside temperatures in Section 1.
} 2o b2
10 J
L
3600
7200
i
10800
i
14400
i
18000
i
21600
Time (s) Fig. 5. C o m p a r i s o n of measured and calculated vertical displacements at sensor 1 (Doel 2).
tween the calculated and measured values with an excellent synchronization in time. This validates the model of the line and thus the forces and moments which are obtained by the combination of the unit load cases. The stresses at the most highly loaded point (which for Doel 2 is located in the elbow under
M. Guyette, M. De Smet / Nuclear Engineering and Design 153 (1995) 183-195 -
-
T5.1
. . . . .
T5.2
.........
T5.3
. . . . . . .
T5.4
.......
T5.5
191
160
220
12°
200 80 ~180 ~ 160 140
0
~ 120
,~ , "
< -40
o~
~ 100 80 6O 64800
.,~"
-80 :-,.:/s
-120
. . . . .
3600
68400 72000 75600 79200 82800 86400
i
. . . . .
t
25200
the pressurizer) may then be calculated by summing the local stresses caused by the stratification with the stresses resulting from the forces and moments in the line. In the calculation of the stresses, use is made of stress intensification factors either directly extracted from the code or defined with the help of a detailed finite element model of the pipe particular discontinuity. In addition, the stresses are multiplied by an uncertainty factor of 1.2 to cover conservatively all the uncertainties and approximations introduced in the calculations. Fig. 9 shows the time evolution of the axial stress at this highest load point for a cool-down (the plot covers a little more than 24 h) while Fig. 10 is relative to the hoop stress at the same point.
6. Contacts with the line environment For each potential point of contact, one defines a local axis going through the axis of the pipe and normal to the contact surface. The origin of this axis is on the position of the pipe axis before any loading is exerted on it. Denoting by d~ the displacement along the local axis of the contact point i, three possible situations may occur:
i
,
,
~
i
68400
i
i
i
i
90000
Time (s)
Time (s) Fig. 8. Evolution in time of the measured outside temperatures in Section 5 (limited to the first five thermocouples).
. . . . .
46800
Fig. 9. Axial stress distribution in the elbow below the pressurizer during a cool-down. 140 10o
"~ ~
20 -20 -60 -100 -140
. . . . .
3600
'
'
25200
'
~
'
'
L
~
J
~
46800
~
,
I
L
68400
~
,
,
L
t
~
~
90000
Time (s) Fig. 10. Hoop stress distribution in the elbow below the pressurizer during a cool-down.
(1) gF < d i < g~-, where gF and gT are the limits between which the pipe is free and the interference force Fi therefore equals 0; (2) di > g i +, for which the interference force is given by F i =Kg(g~" - d i )
(15)
where Kg is the stiffness of the contact surface;
M. Guyette, M. De Smet /Nuclear Engineering and Design 153 (1995) 183 195
192
(3) d~ < g [ , for which the interference force is given by F, = g,.(g i - d , )
(16)
Similarly to what has been done for the responses of interest, the displacements in the direction of the local axis may be extracted from the above-mentioned unit load cases and the free displacements (displacements in absence of any contact) at the potential contact points may be calculated for each time in the transient, in a manner similar to what has been done in Eq. (14). These free displacements may therefore be put in the form flrree = D~[i + Di-~" + D~6
(17)
with a similar definition for Dp etc. as for Vp etc. Each potential contact point is assumed to be in one of the three states previously defined. The forces at the contact points and the displacement in the direction of the local axis will be calculated on the basis of this assumption. From the values of the forces and displacements it will be possible to verify whether the assumed state is correct. If not, another assumption has to be made on the state of the contact points and the process must be iterated. One will first describe how the forces and displacements are obtained. Eqs. (15) and (16) for all the potential contacts which are assumed to be in states (2) or (3) are grouped in F* = K * ( g * - d*)
(18)
The asterisk is used to indicate that the various vectors and matrices are not the full quantities (for all potential contacts) but are limited to the contacts in state (2) or (3). The total displacements in the contact local coordinates may be written as d
= dfree -f-
DvF
(19)
where the DFO is the displacement in local coordinates at contact point i for a unit force at contact j. These values may be obtained from the global model of the line by calculating a series of load cases corresponding to a unit force applied at each potential contact in the direction of the local axis.
Bringing Eq. (19) into Eq. (18), we obtain alter transformations F*(K* ' + D * )
l(g, - d f, ~ )
(20)
The displacements may be calculated by Eq. (19). The assumptions made on the potential contact states are correct if: • for a state (1) contact point g70 In case one or several of the states are not correct, a new trial has to be done. Thanks to the small number of potential contact points (usually less than five) a very crude way of iterating (assuming successively all possible situations for each potential contact point) may be performed, the converged situation of a calculation time being assumed as first trial for the next time. Finally the responses of interest are calculated by an equation similar to Eq. (14), taking the interference forces into account: v = Vpli + V+T + V~ii + VvF + Vpp
(21)
Fig. 11 shows the layout of the Doel 4 surge line. Several contacts with pipe whip restraints have been observed on this line. The most important is the restraint WP6 (see Fig. 11) which presents a very stiff contact with a downwards clearance of 31 mm. Fig. 12 shows the comparison of the measured and calculated vertical displacements at sensors 1 and 3 during a typical transient, while Fig. 13 gives the same comparison for the displacements at sensors 2 and 4. The comparison of the calculated and measured displacements is good. The differences are probably due to other contact points which were neglected in the calculation. Fig. 14 shows the evolution of the contact force during the transient. The maximum force of more than 50 kN will obviously cause significant moments in the line. A calculation of the line in absence of the restraint shows that the deflection
M. Guyette, M. De Smet / Nuclear Engineering and Design 153 (1995) 183-195
-
Calc. - . . . . dp. 4
-
Calc. dp. 2
193
Meas. - - - - - Meas. dp. 4 dp. 2
40 30 T1.1
."0"-
T1.4~-~"
~e">
T2.3 N~ ~ " / ~ , ~
L -'v
~ 2o 10
Displ.
o .~ -10 b2 -20
-30 -40 0
i
i
1200
2400
i
i
i
i
3600
4800
6000
7200
Time (s) Fig. 13. C o m p a r i s o n of measured and calculated displacements at sensors 2 and 4 (Doel 4). Fig. 1 I. Layout of the Doel 4 surge line showing the locations of the temperature and displacement sensors.
60000
50000 -
Calc. - . . . . dp. 1
-
Calc. dp. 3
Meas. - - - - - Meas. dp. 1 dp. 3
40000 z
10 ~
30000
5 U 20000
~ o k.,-~----------~-~,
~,~,e "-- - . ~ _ 10000
i
~ -10
1200
i
i
2400
i
3600
4800
i
I
6000
7200
Time (s)
~ -15
Fig. 14. Force at the gap contact point WP6 (Doel 4). -20
~
0
i
i
1200
2100
i
i
i
i
i
3600
4800
6000
7200
Time (s) Fig. 12. C o m p a r i s o n of measured and calculated displacements at sensors I and 3 (Doel 4).
at the location of the restraint would be 2.5 times as large as what it is with the restraint present.
7. Fatigue analysis Once the local stresses and the force and moment distributions are calculated, the peak stresses at the points o f interest may easily be evaluated leading to time history curves o f the various stress components. The fatigue analysis based on the
194
M. Guyette, M. De Smet / Nuclear Engineering and Design 153 (1995) 183 195
rainflow method is then straightforward. The following steps are performed in turn: • identify first the extremes of the stress diagram; • identify the smallest cycles first, count them and eliminate them from the time history diagram and proceed further to finally take into account the largest cycles.
Electrabel kindly allowed the publication of measurements made in their Doel 2 and Doel 4 power plants.
A p p e n d i x A: N o m e n c l a t u r e
c~ d 8. C o n c l u s i o n s d free
• Provisional measurement devices were installed on the surge line of the Doel 2 and the Doel 4 power plants to monitor the stratification in the line. From measurements obtained during the heat-up and cool-down periods, it has been shown that a detailed fatigue monitoring of the actual transients occurring in the line is feasible. The main problems which were solved are as follows: --the complexity of the temperature and stress fields in space as well as in time; --the determination of the temperature and stress fields from the knowledge of the outside wall temperatures; - - the calculation of the global line behaviour and its addition to the local effects; - - the effect of intermittent contacts with pipe whip restraints. • The implementation of a permanent monitoring system may be contemplated as an alternative for the justification of the additional fatigue induced by the stratification on the surge lines. • The usage factors obtained by this method are notably smaller than those calculated by classical design analysis. The period of observation, limited to a few heat-ups and cool-downs, is, however, too short to draw a final conclusion on this last point.
d~ D
F
F, g gi g~ Gk h k K
K, l m
Mok(r)
Mlk (r)
P Ok Acknowledgements
r
Sik(r) The authors want to thank Dr. C. Leblois and Dr. L. Geraets, both from Tractebel, for the fruitful discussions they had with them on the subject.
Sok(r) T(r, O, t)
function defined by Eq. (9) vector containing the displacements of the line at the contact points vector containing the free displacements of the line at the potential contact points displacement of contact i in vector d matrix containing displacements (in local axes) at the contact points, due to unit loads at all potential contacts vector containing the forces acting on the pipe at contact point i force acting on the pipe at contact point i vector containing clearance values at the contact points clearance between pipe and restraint i on the negative side of local axis clearance between pipe and restraint i on the positive side of local axis function defined by Eq. (10) heat transfer coefficient current harmonic index stiffness matrix at the potential contact points stiffness of the contact point i current eigenform index maximum number of eigenforms taken into account steady state temperature distribution in the harmonic k for a unit fluid temperature in same harmonic eigenform of index ! for the harmonic k internal pressure constant defined by Eq. (11) current radial coordinate stress caused by the application of the temperature distribution Mik(r) steady state stress caused by the temperature field Mok(r) metal temperature
M. Guyette, M. De Smet / Nuclear Engineering and Design 153 (1995) 183-195 t
Tk(r, t)
Tfk(t)
Uk(r, t) v
V
wi
current time vector containing the mean temperature in each measurement section coefficient of the kth harmonic in the Fourier development of the metal temperature component of harmonic k of the fluid temperature inverse transfer function of harmonic k vector containing response of interest matrix containing response in the points of interest due to unit loads at the measurement sections weights for the calculation of the fluid temperature
Greek letters vector containing the free curvature at each measurement section
0 K 2 ak(r, t)
q51k
195
vector containing the end displacements in each coordinate direction current azimuthal coordinate metal thermal diffusivity metal thermal conductivity sampling time interval coefficient of harmonic k of the stress distribution (represent any stress component) eigenvalue of index l for harmonic k
References B.A. Boley and J.H. Weiner, Theory of Thermal Stresses, Krieger, Malabar, FL, 1985. H.S. Carslaw and J.C. Jaeger, Conduction of Heat in Solids, Clarendon, Oxford, 1959. M. Guyette, Prediction of fluid temperature from measurements of outside wall temperatures in pipes, ASME J. Pressure Vessel Technol., 116 (1994) 179 187. S.P. Timoshenko and J.N. Goodier, Theory of Elasticity, McGraw-Hill, Kogakusha, Tokyo, 1970.
Nuclear Engineering and Design 153 (1995) 183 195
Nuclear Engineering and Design
Inverse and direct transfer functions for the fatigue follow-up of piping systems submitted to stratification M. Guyette, M. De Smet Tractebel Energy Engineering, Avenue Ariane 7, 1200 Brussels, Belgium
Abstract In this paper we outline a methodology to assess the fatigue induced in piping systems submitted to thermal stratification. More specifically, the transformation from the measured outer wall temperature time histories to stress time histories in any point of the line is treated. By means of inverse transfer functions, the fluid temperature distribution is calculated from the outside wall temperatures measured in a limited number of temperature sections. Using direct transfer functions, the local stresses due to stratification may be determined as well as the pipe free curvatures and the pipe free axial strains. Using a finite beam element model of the line, the global response of the line (in terms of displacements or stresses) due to the applied curvatures, axial strains, end point displacements, internal pressure and possible contacts with the pipe environment may be determined. The method is illustrated for the surge lines of the Doel 2 and Doel 4 nuclear power plants. An excellent correlation is found between measured and calculated displacements. Typical stress time histories are shown for a plant cool down.
1. Introduction Thermal stratification has been observed in several pressurized water reactor systems for a couple of years. Thermal stratification typically occurs in the surge line and the main feed water lines but other lines such as the pressurizer spray line may also exhibit the phenomenon. The large number of thermal cycles of great amplitude produced by stratification has raised some concern about the damage induced by fatigue in these lines. The piping systems affected by stratification have been equipped with thermocouples mounted on the outer pipe wall and with displacement
sensors. Temperature and displacement time histories have been recorded during several heat-ups and cool-downs. In view of a fatigue analysis one has to pass from the recorded outer wall temperature time histories to the stress time histories at the most highly loaded points in the piping. The procedure which has been used is outlined in the following sections.
2. Main problems to be solved The physical parameters which influence the fatigue of a highly loaded point are usually very
0029-5493/95/$9.50 © 1995 Elsevier Science S.A. All rights reserved SSDI 0029-5493(94) 00826-K
184
M. Guyette, M. De Smet / Nuch, ar Engineering and Design 153 (1995) 183 195
limited in number. In most cases, the knowledge of the local pressure and of a few values of fluid temperature (usually less than three) is sufficient to allow the evaluation of the local stresses by use of direct transfer functions. The problem of investigating the thermal stratification in a pipe is notably more complex since, as a result of the non-linearity of the temperature distribution over a cross-section, even in an unrestrained pipe, stratification induces local stresses and a pipe curvature. The stresses at a given location in a restrained pipe submitted to stratified flow are therefore a function not only of the local temperature distribution but also of the curvatures induced by the stratification at all other locations. Usually the stratification pattern is variable along the length of the line. The local temperature distribution as well as the distribution of the curvature along the line (or the linear thermal gradient, which is proportional to the curvature) must therefore be gained from the monitoring sensors. To reach a sufficient accuracy, a minimum of five measurement points per section are required and a minimum of four or five measurement sections along the length of the pipe has to be contemplated. Ideally, the measurement should give the fluid temperature which is difficult, the measurements of the pipe outer wall temperature being easier by far to perform. This leads either to the assumption that the external wall temperatures are a fair representation of the water temperatures (which appear to be inaccurate and non-conservative particularly for transients) or to the development of an algorithm allowing the calculation of the fluid temperature from the external wall values. The development of this algorithm requires the knowledge of the film coefficient between the fluid and the wall, the value of which is usually not known with good accuracy. Large inaccuracies on the fluid temperatures could therefore be generated by such a process. In the sequel, one is, however, not directly interested by the fluid temperature but rather by the temperature distribution in the metal. Provided that use is made, in the direct transfer functions (allowing calculation of the temperatures, stresses and strains in the
metal from the fluid temperatures), of the same film coefficient values, the inaccuracies mentioned above are fully compensated. The fluid temperature must be regarded as an intermediate result without direct interest in the calculation chain. This point will be commented on in more detail in Section 4. From a physical point of view, this deduction may appear as relatively evident since the whole process consists of starting from temperatures measured in some particular points in the metal to deduce other temperatures in the same piece of metal. Once the local curvatures are known, they must be introduced in a global model of the line which will be used to determine the forces and moments at all locations of interest. The problem is complicated by the fact that the gaps at the pipe whip restraints are not always sufficient to accommodate the large displacements caused by the stratification. Contact forces with these restraints take place at some locations and completely disturb the distribution of moments. The main problems to be solved in view of a fatigue analysis based on measured outer wall temperature time histories are therefore --- the development of a practical way of coping with the complex temperature and stress patterns both in a cross-section and axially, - - the calculation of the fluid temperature distribution from the measured outside wall temperature, - - the evaluation of the global behaviour of the line (influenced by the local conditions in each cross-section), and the presence of gap contact situations at some pipe whip restraints. The way by which these problems have been solved is described in detail in the next sections.
3. C o m p l e x i t y o f the temperature and stress fields
The temperature field in the pipe presents important variations in azimuthal and radial directions while the variations are much smoother in the axial direction. In addition, the temperature field is strongly variable in time. Based on these assumptions and neglecting the small variations in
M. Guyette, M. De Smet /Nuclear Engineering and Design 153 (1995) 183-195
axial direction, the temperature field will satisfy the differential equation a2T 1 a T 1 a2T ~ar + -r ~ r + r 2 ~02
1 aT K at
(1)
where T is the metal temperature (see also nomenclature in Appendix A) and K is the metal thermal diffusivity. The following boundary conditions will be applied:
185
eigenmodes and eigenfrequencies in structural dynamics will be used for each harmonic. The time history response of the temperature at a given time t for the point of coordinates r and 0 in the pipe will be
Tk(r't)=[M°k(r)--~M'k(r)lTfk(t),=l 4-
Y, K4;,,M~,(r) /=l
2(ar~ --
k~rJ
= htT(r,, O, t) - Tf(O, t)]
(2)
....
1
x
-_o r=r e
which assume a heat exchange, with a film coefficient h, between the fluid and the metal at the inner surface and a perfect insulation at the outer surface. The initial temperature is assumed to be uniformly equal to 0. This last assumption could be considered as an oversimplification of the problem. Taking, however, into account that the transients that will be followed are extremely long compared with the period needed for the initial conditions to die out, this assumption has a negligible effect on the computed results. In each cross-section, the temperature presents a mirror symmetry with respect to the vertical plane. Based on these observations, it has been decided to express the temperature field in the pipe as a Fourier series of the azimuthal angle around the pipe axis. The temperature will thus be expressed as
T(r, O, t) = ~ T,(r, t) cos kO
(4)
k=0
In order to keep the equations linear, the thermal as well as the mechanical properties of the pipe material such as the thermal conductivity, the thermal capacity, the Young's modulus and the thermal expansion coefficient will be considered to be independent of the temperature. The thermal conduction equations in transient regime will be solved for each harmonic of the Fourier series. A technique of eigenvalues and eigenfunctions (see Boley and Weiner (1985) or Carslaw and Jaeger (1959)) similar to that of the
e x p [ - K~b~k(t - v)lTfk(r) dr
(5)
where Mik(r) is the eigenform of index l for the harmonic k, ~bCkthe eigenvalue of index l for the harmonic k, ~ the thermal diffusivity of the pipe material, Tfk(t) the component of harmonic k of the fluid temperature and Mok(r) the steady state temperature distribution of harmonic k for a fluid temperature equal to unity in the same harmonic. The eigenforms M~k(r) are normalized in such a way that
mok(r) = ~ mkl(r)
(6)
/=1
The eigenforms are combinations of first- and second-kind Bessel functions of the same order as the harmonic to which they belong. The derivation of the stresses caused by the temperature field is performed by bringing the value of this temperature field into the thermoelastic equations. Use is again made of the mirror symmetry with respect to the vertical plane. In a first step, a strict plane strain hypothesis is used for the determination of the stress field. In a second step, this plane strain hypothesis is transformed into a generalized plane strain hypothesis, under the assumption that the net axial force and the bending moment perpendicular to the vertical symmetry plane are both equal to zero. It is assumed that the thermoelastic equations are not explicitly functions of the time, which corresponds physically to neglect of the inertia forces caused by the thermal expansion. Because of the slowness of the temperature changes on the one hand and the smallness of the thermal dis-
M. Guyette, M. De Smet / Nuclear Engineering and Design 153 (1995) 183-195
186
placements on the other hand, these inertia forces are indeed negligible. The solution technique is based on the theory of the thermoelastic potential (see Timoshenko and Goodier (1970)) and on a Fourier development of all the functions (displacements, stresses, strains) intervening in the equations. Beyond the generalized plane strain conditions mentioned above, the boundary conditions imposed on the stress solution are simply that the radial normal stress and the shear stress a,.o vanish on both the inner and outer surfaces of the pipe. For each Mik(r ) eigenform (l being the mode subscript and k the harmonic subscript), the thermoelastic theory will give stress functions Stk(r) for each stress component which, after combination, gives the stress response of harmonic k:
4. Calculation of the azimuthal fluid temperature distribution from the outer wall distribution: inverse transfer functions
The equations above may only be used if the fluid temperature is known. In fact as already explained, the outside wall temperatures are, in most cases, the only ones available. A formulation has therefore to be developed to obtain the fluid temperature from the outside wall temperatures. In such a process, the transfer functions used may be called "inverse" since they allow us to obtain the cause of the phenomenon from the knowledge of its effects in opposition to normal or "direct" transfer functions which permit us to calculate the effects of a phenomenon from its causes. Each of Eqs. (5) relative to a particular harmonic may be written for the external radius r e and put in the form
Tk(re, t) = Ck(re)Tfk(t ) +
Gk(r e, t -- r)Tfk(r) dr
+ ~ ~'O~kS,k(r)
(8)
I-I
x
;o
exp[ - ~CC~k(t -- r)]Tfk(r) dr
(7)
where S~k(r) is the modal stress caused by the application of temperature distribution M~k(r) and Soa.(r) is the steady state stress caused by the temperature field Mok(r ). The last equation is valid for any one of the stress components. A similar formulation may also be used for the displacements, strains etc. The stresses are first calculated under the assumption of plane strain conditions. Under these conditions the harmonic 0 usually causes an axial force. In a second step, to obtain the pipe free expansion, an axial force of opposite sign is superimposed. The same reasoning may be used concerning the bending moment caused by the harmonic I. The superposition of a moment of opposite sign causes a curvature of the pipe called hereafter the free curvature. The above formulation allows, for any section where the fluid temperature is known, the local stresses, the free axial expansion and the free curvature to be calculated.
with Ck(re) = Mok(rc)- ~ Mlk(re)
(9)
I-I
Gk(re, t) = ~ ~Cb21Mlk(re) exp(-~c~b2/t)
(10)
l-I
r e being a constant, Ck a constant and Gk a function of t only. Eq. (8) is the basis to obtain the inverse transfer functions. In the case where Tk(re, t) is known and Tfk(t) unknown, it may be considered as an integral equation in Tfk, the solution of which may be written as
Tfk(t ) = Qk(re)Tk(re, t) +
f+
Uk(re, ti~)Tk(re, r) dr
(11)
The mathematical developments leading to the values of the constant Qk and to the expression of Uk are beyond the scope of this presentation. They are described in detail in Guyette (1994). The Uk function is the inverse transfer function. As this function decreases sharply as t tends to
M. Guyette, M. De Srnet / Nuclear Engineering and Design 153 (1995) 183-195
+ oo or - oo, the final expression of Tfk(t) may be put in the form Tfk(t ) ~ Ok(re)Tk(re, t) -Jr
f
t+A
Uk(Fe, t -- r)Tk(re, z) dr
(12)
,,)t-A-
where A - and A ÷ are chosen in such a way as to give a good engineering accuracy to the expression of Tfk. A function Uk and a constant Qk are obtained for each harmonic in turn. The inversion process is therefore executed in the following steps: (1) perform a Fourier analysis of the outer wall temperature to obtain the Fourier coefficients of the external wall temperature; (2) knowing the time history of each of these outer wall Fourier coefficients, obtain, with the help of Eq. (12), the response in time of the Fourier coefficient for the fluid temperature; (3) to obtain the fluid temperature at a given azimuthal position, assess the value of the Fourier series at this point. The inversion process is therefore fully two dimensional. The formulation which is followed warrants the uniqueness of the solution. Contrary to what is usually given in the literature on the inverse problem, the solution presented is based on the use of future times: the integration domain in Eq. (12) extends up to t + A +. This has the advantage that the solution is very stable and accurate (which is not the reputation of usual inverse problems) but it is necessary to retain measurements until the time t + A + before the temperatures and stresses at time t can be evaluated. A + being usually of the order of a few hundred seconds, this may not be considered as a true drawback for the type of problem which is solved here. From the theoretical point of view, the main approximation is the limitation of the number of eigenfunctions taken into account in the inversion process of each harmonic. It is, however, shown, in Guyette (1994), that, from a practical point of view, the accuracy on the results of the inversion process, i.e. the bulk fluid temperatures, is mainly limited by the fluctuations of the outer
187
wall temperatures caused by the measurement devices: beyond systematic errors, the measurement chains usually produce erratic fluctuations which are of the order of a few tenths of a kelvin (for temperatures between approximately 150 and 350 °C). These erratic fluctuations induce a numeric noise which is of the order of 0.5 K on the results of the inversion. An important point to be commented on is the influence of the film heat transfer coefficient on the results of the direct and inverse transfer functions. This film coefficient is taken into account in both types of functions and any uncertainty on its values may strongly affect the results. However, when use is concurrently made of inverse and direct transfer functions and provided that the same value of the film coefficient is used for the two types of functions, a compensation of the error occurs. For instance, if a very large film coefficient is used, the inverse transfer functions will predict the inside wall temperature in place of the bulk fluid temperature. However, if the same value of the film coeffÉcient is used in the direct transfer functions, the calculated variables, which may be temperatures, displacements, stresses or strains, will assume that the starting temperature is the inside wall temperature and the calculated values will be correct. As a matter of fact, the process described above, combining inverse and direct transfer functions, may be seen as a way of evaluation of the heat flux exchanged with the fluid at the inner wall, from which the temperature at any point in the metal may be calculated. This process does not attempt to make predictions of what occurs in the fluid itself. This conclusion is of course limited to the simulation of the phenomena which have a measurable impact on the outer wall temperatures. Rapid fluctuations of the fluid temperatures which are undetectable at the outer wall surface, such as what is known in the literature as the striping phenomenon, are not accounted for in the abovedescribed process. Their contribution to the usage factor of the equipment has been shown to be rather limited. The expression (12) is particularly useful when the outer wall temperatures are sampled at regular time intervals. If one makes the hypothesis either
188
M. Guyette, M. De Srnet / N u c l e a r Engineering and Design 153 (1995) 183-.-195
of a linear variation in the outer wall temperature (Co continuity) or a more elaborate C, continuity (using third-degree polynomials), the expression (12) may be transformed into T,-,(t) = ~, w,T~(t + i~)
(13)
i
where ~ is the sampling time interval, the index i being such as to cover the time period between t - A - and t + A +. The w~ are weights depending on the value of the constant Q~ and of the function U~. As the latter are unique for a given geometry and material, the w~may be calculated once for all, which allows a fast calculation of the fluid temperatures from the outer wall temperatures. The following example illustrates the good accuracy obtained in the use of the inverse transfer functions. Figs. 1 and 2 show the measured outside wall temperatures recorded during a typical stratification transient in a cross-section of the surge line of the Doel 2 power plant. The measured values indicated by diamond symbols on the figures correspond to five thermocouples located 45 ° apart on the outside wall of the piping. A fitting of the
220 200 180 ~,~:
"---222. '2~
0 16o = 140 ~ 120 100 80 60
i
0
i
i
150
300
i
I
i
i
450
600
750
900
Time (s) Fig. 1. Accuracy check of the inverse transfer functions: O, measured outside temperatures; the fine curves are the calculated inside temperatures (obtained by the use of inverse transfer functions) while the thick curves are the calculated outside wall temperatures (derived with direct transfer functions).
220
, •,,, .v,~ ~,~,"'~ <'" <"""" 2
200 180 L
~ / / ~ ' " ~"
~" o 160 140
K 120 ~,,
•
lOO 80 60 900
i
i
i
1050
1200
1350
i
i
1500
i
i
i
1650
1800
Time (s) Fig. 2. Accuracy check of the inverse transfer functions (continued). See comments from Fig. 1.
outside wall temperatures by a Fourier series is first performed for all measured times. With the help of inverse transfer functions, applied on each term of the Fourier series, the development of the inner wall temperature is obtained. Assessing the series for the five azimuthal positions of the thermocouples, one obtains the inside wall temperatures just in front of the thermocouples. These are shown as fine line curves on Figs. 1 and 2. Using now direct transfer functions, the outside wall temperatures may be recalculated. The values corresponding to the five thermocouple positions are plotted as thick line curves on Figs. 1 and 2. Any error or inaccuracy in the inner wall temperatures would bring a discrepancy between the starting values of the outer wall temperatures and the final calculated values. The figures show a very good agreement with a discrepancy well below 1 K (thick lines marked by diamonds) except at the onset of the sharp transient where differences are about 2 K. The figures also clearly show that the inner wall temperature curves anticipate with a sharper slope the trends seen for the outside wall temperature, which conforms to the physics of the phenomenon. Examination of Figs. 1 and 2 shows that the inside wall temperature may be as much as 40 K
M. Guyette, M. De Smet / Nuclear Engineering and Design 153 (1995) 183-195
higher than the outer wall at some moments of the transient. This illustrates the errors which could have been made if the outer wall temperature had been used in place of the inside wall temperature.
189
I Pressu- I I Section 1 Section 2
Section3
Section 4
Section5
I influence influence I function function
influence function
influence function
influence function
.'--'--<.T_-:::-:----: 0
0
0
:-:..:
=t'3 0
~_1
5. Global behaviour of the line
Fig. 3 shows the layout of the surge line of the Doel 2 power plant with the location of the provisional measurement devices (thermocouples and displacement sensors) used to monitor some heat-ups and cool-downs. This provisional measuring installation is typical (number and position of the captors) of what may be used for permanent monitoring. The use of the inverse transfer functions combined with the stress and strain analysis described here above allows us to calculate the local stresses, the axial free expansion and free curvature of the pipe in each measurement section. To assess the global behaviour of the line, a model of the latter is built; this model may be made with the help of either a specific piping
T1.2 TI.1
T2.2~.
TI.S \
Section3 ~ s~nsor3
I
T4.2[ ~
4
Fig. 3. Layout of the Doel 2 surge line showing the locations of the temperature and displacement sensors.
Fig. 4. Schematic representation of the influence functions of the various sections.
calculation program provided that it is able to assess the effects of transverse thermal gradients (causing a curvature of the pipe elements) or a general-purpose finite element program. In the second case, care must be taken that the elements representing the pipe elbows and other specific piping fittings reflect the true flexibility. It is assumed that the axial free expansion and the free curvature vary linearly between measurement sections. Influence functions are defined for each section as shown in Fig. 4. The function relative to a particular section varies from 1 in the considered section to 0 in the adjacent sections. A series of responses of interest are defined for which the time evolution during the transients is wanted. These responses may be displacements, stresses, strains or forces at particular locations in the piping system. These responses are assumed to be stored according to a well-defined order in a vector v. For each measurement section, load cases corresponding to unit curvature and to unit average temperature, distributed in the line according to the above-mentioned influence functions, are calculated. The values of the responses of interest are extracted from the results of these load cases and stored column by column in the matrix V~ for what concerns the unit curvatures and in V¢ for the mean unit temperatures. Load cases corresponding to unit end displacement in each of the coordinate directions are also calculated and the responses of interest are stored columnwise in the V s matrix. As, from the previous calculations, the mean temperatures and the curvatures are known for each time step in the transient, the values of the
M. Guyette, M. De Smet / Nuclear Engineering and Design 153 (1995) 183 195
190
responses of interest may simply be calculated for each step by
Meas. - dp. 2
(14)
v = Vp[~ + V i . T + V a 6 + Vpp
where [i and I" are vectors containing respectively the curvature and the mean temperature in the measurement sections, 6 is the vector containing the end displacements for the considered time, Vp is the vector of the responses of interest for a unit internal pressure and p is the internal pressure. The accuracy of this global model may easily be checked by performing a comparison between the measured values of the displacements and the calculated response of the model. As shown in Fig. 3, in addition to the five measurement sections equipped with thermocouples, the vertical displacement of the line is measured at three locations. These displacements are compared with the calculated displacements in Figs. 5 and 6 for a typical transient the recorded temperatures of which are given in Fig. 7 for what concerns Section 1 and in Fig. 8 for Section 5 (the temperature variations in Sections 2 - 4 are similar to those of Section 1 and are not shown). Examination of Figs. 5 and 6 shows a very good agreement be-
Calc. - - - - - Meas. - . . . . . . . Calc. dp. 2 dp. 3 dp, 3
7O
60
50 4o
.~ 30
ilo 2°
-10 0
'
,
,
,
3600
7200
10800
14400
Calc. dp. 1
,
Time (s) Fig. 6. C o m p a r i s o n of measured and calculated vertical displacements at sensors 2 and 3 (Doel 2).
-
-
T1.1
- -
-
T1,2
.........
Tt.3
.
.
.
.
T1.4
-
-
T1.5
220
200
/
,s0 Meas. dp. 1
,
18000 21600
I j,/
i140
70
12o
6O
(~ 100
f/;:? t,/' i~ii :
i,,
J.: i ! f
,
t
,;
ii ,'
Z" 50 60
~ 40
,
64800
68400
,
,
,'.-r.-,'
,
,
7 2 0 0 0 7 5 6 0 0 79200 82800
,
86400
r~
Time (s) "~ 30 Fig. 7. Evolution in time of the measured outside temperatures in Section 1.
} 2o b2
10 J
L
3600
7200
i
10800
i
14400
i
18000
i
21600
Time (s) Fig. 5. C o m p a r i s o n of measured and calculated vertical displacements at sensor 1 (Doel 2).
tween the calculated and measured values with an excellent synchronization in time. This validates the model of the line and thus the forces and moments which are obtained by the combination of the unit load cases. The stresses at the most highly loaded point (which for Doel 2 is located in the elbow under
M. Guyette, M. De Smet / Nuclear Engineering and Design 153 (1995) 183-195 -
-
T5.1
. . . . .
T5.2
.........
T5.3
. . . . . . .
T5.4
.......
T5.5
191
160
220
12°
200 80 ~180 ~ 160 140
0
~ 120
,~ , "
< -40
o~
~ 100 80 6O 64800
.,~"
-80 :-,.:/s
-120
. . . . .
3600
68400 72000 75600 79200 82800 86400
i
. . . . .
t
25200
the pressurizer) may then be calculated by summing the local stresses caused by the stratification with the stresses resulting from the forces and moments in the line. In the calculation of the stresses, use is made of stress intensification factors either directly extracted from the code or defined with the help of a detailed finite element model of the pipe particular discontinuity. In addition, the stresses are multiplied by an uncertainty factor of 1.2 to cover conservatively all the uncertainties and approximations introduced in the calculations. Fig. 9 shows the time evolution of the axial stress at this highest load point for a cool-down (the plot covers a little more than 24 h) while Fig. 10 is relative to the hoop stress at the same point.
6. Contacts with the line environment For each potential point of contact, one defines a local axis going through the axis of the pipe and normal to the contact surface. The origin of this axis is on the position of the pipe axis before any loading is exerted on it. Denoting by d~ the displacement along the local axis of the contact point i, three possible situations may occur:
i
,
,
~
i
68400
i
i
i
i
90000
Time (s)
Time (s) Fig. 8. Evolution in time of the measured outside temperatures in Section 5 (limited to the first five thermocouples).
. . . . .
46800
Fig. 9. Axial stress distribution in the elbow below the pressurizer during a cool-down. 140 10o
"~ ~
20 -20 -60 -100 -140
. . . . .
3600
'
'
25200
'
~
'
'
L
~
J
~
46800
~
,
I
L
68400
~
,
,
L
t
~
~
90000
Time (s) Fig. 10. Hoop stress distribution in the elbow below the pressurizer during a cool-down.
(1) gF < d i < g~-, where gF and gT are the limits between which the pipe is free and the interference force Fi therefore equals 0; (2) di > g i +, for which the interference force is given by F i =Kg(g~" - d i )
(15)
where Kg is the stiffness of the contact surface;
M. Guyette, M. De Smet /Nuclear Engineering and Design 153 (1995) 183 195
192
(3) d~ < g [ , for which the interference force is given by F, = g,.(g i - d , )
(16)
Similarly to what has been done for the responses of interest, the displacements in the direction of the local axis may be extracted from the above-mentioned unit load cases and the free displacements (displacements in absence of any contact) at the potential contact points may be calculated for each time in the transient, in a manner similar to what has been done in Eq. (14). These free displacements may therefore be put in the form flrree = D~[i + Di-~" + D~6
(17)
with a similar definition for Dp etc. as for Vp etc. Each potential contact point is assumed to be in one of the three states previously defined. The forces at the contact points and the displacement in the direction of the local axis will be calculated on the basis of this assumption. From the values of the forces and displacements it will be possible to verify whether the assumed state is correct. If not, another assumption has to be made on the state of the contact points and the process must be iterated. One will first describe how the forces and displacements are obtained. Eqs. (15) and (16) for all the potential contacts which are assumed to be in states (2) or (3) are grouped in F* = K * ( g * - d*)
(18)
The asterisk is used to indicate that the various vectors and matrices are not the full quantities (for all potential contacts) but are limited to the contacts in state (2) or (3). The total displacements in the contact local coordinates may be written as d
= dfree -f-
DvF
(19)
where the DFO is the displacement in local coordinates at contact point i for a unit force at contact j. These values may be obtained from the global model of the line by calculating a series of load cases corresponding to a unit force applied at each potential contact in the direction of the local axis.
Bringing Eq. (19) into Eq. (18), we obtain alter transformations F*(K* ' + D * )
l(g, - d f, ~ )
(20)
The displacements may be calculated by Eq. (19). The assumptions made on the potential contact states are correct if: • for a state (1) contact point g7
(21)
Fig. 11 shows the layout of the Doel 4 surge line. Several contacts with pipe whip restraints have been observed on this line. The most important is the restraint WP6 (see Fig. 11) which presents a very stiff contact with a downwards clearance of 31 mm. Fig. 12 shows the comparison of the measured and calculated vertical displacements at sensors 1 and 3 during a typical transient, while Fig. 13 gives the same comparison for the displacements at sensors 2 and 4. The comparison of the calculated and measured displacements is good. The differences are probably due to other contact points which were neglected in the calculation. Fig. 14 shows the evolution of the contact force during the transient. The maximum force of more than 50 kN will obviously cause significant moments in the line. A calculation of the line in absence of the restraint shows that the deflection
M. Guyette, M. De Smet / Nuclear Engineering and Design 153 (1995) 183-195
-
Calc. - . . . . dp. 4
-
Calc. dp. 2
193
Meas. - - - - - Meas. dp. 4 dp. 2
40 30 T1.1
."0"-
T1.4~-~"
~e">
T2.3 N~ ~ " / ~ , ~
L -'v
~ 2o 10
Displ.
o .~ -10 b2 -20
-30 -40 0
i
i
1200
2400
i
i
i
i
3600
4800
6000
7200
Time (s) Fig. 13. C o m p a r i s o n of measured and calculated displacements at sensors 2 and 4 (Doel 4). Fig. 1 I. Layout of the Doel 4 surge line showing the locations of the temperature and displacement sensors.
60000
50000 -
Calc. - . . . . dp. 1
-
Calc. dp. 3
Meas. - - - - - Meas. dp. 1 dp. 3
40000 z
10 ~
30000
5 U 20000
~ o k.,-~----------~-~,
~,~,e "-- - . ~ _ 10000
i
~ -10
1200
i
i
2400
i
3600
4800
i
I
6000
7200
Time (s)
~ -15
Fig. 14. Force at the gap contact point WP6 (Doel 4). -20
~
0
i
i
1200
2100
i
i
i
i
i
3600
4800
6000
7200
Time (s) Fig. 12. C o m p a r i s o n of measured and calculated displacements at sensors I and 3 (Doel 4).
at the location of the restraint would be 2.5 times as large as what it is with the restraint present.
7. Fatigue analysis Once the local stresses and the force and moment distributions are calculated, the peak stresses at the points o f interest may easily be evaluated leading to time history curves o f the various stress components. The fatigue analysis based on the
194
M. Guyette, M. De Smet / Nuclear Engineering and Design 153 (1995) 183 195
rainflow method is then straightforward. The following steps are performed in turn: • identify first the extremes of the stress diagram; • identify the smallest cycles first, count them and eliminate them from the time history diagram and proceed further to finally take into account the largest cycles.
Electrabel kindly allowed the publication of measurements made in their Doel 2 and Doel 4 power plants.
A p p e n d i x A: N o m e n c l a t u r e
c~ d 8. C o n c l u s i o n s d free
• Provisional measurement devices were installed on the surge line of the Doel 2 and the Doel 4 power plants to monitor the stratification in the line. From measurements obtained during the heat-up and cool-down periods, it has been shown that a detailed fatigue monitoring of the actual transients occurring in the line is feasible. The main problems which were solved are as follows: --the complexity of the temperature and stress fields in space as well as in time; --the determination of the temperature and stress fields from the knowledge of the outside wall temperatures; - - the calculation of the global line behaviour and its addition to the local effects; - - the effect of intermittent contacts with pipe whip restraints. • The implementation of a permanent monitoring system may be contemplated as an alternative for the justification of the additional fatigue induced by the stratification on the surge lines. • The usage factors obtained by this method are notably smaller than those calculated by classical design analysis. The period of observation, limited to a few heat-ups and cool-downs, is, however, too short to draw a final conclusion on this last point.
d~ D
F
F, g gi g~ Gk h k K
K, l m
Mok(r)
Mlk (r)
P Ok Acknowledgements
r
Sik(r) The authors want to thank Dr. C. Leblois and Dr. L. Geraets, both from Tractebel, for the fruitful discussions they had with them on the subject.
Sok(r) T(r, O, t)
function defined by Eq. (9) vector containing the displacements of the line at the contact points vector containing the free displacements of the line at the potential contact points displacement of contact i in vector d matrix containing displacements (in local axes) at the contact points, due to unit loads at all potential contacts vector containing the forces acting on the pipe at contact point i force acting on the pipe at contact point i vector containing clearance values at the contact points clearance between pipe and restraint i on the negative side of local axis clearance between pipe and restraint i on the positive side of local axis function defined by Eq. (10) heat transfer coefficient current harmonic index stiffness matrix at the potential contact points stiffness of the contact point i current eigenform index maximum number of eigenforms taken into account steady state temperature distribution in the harmonic k for a unit fluid temperature in same harmonic eigenform of index ! for the harmonic k internal pressure constant defined by Eq. (11) current radial coordinate stress caused by the application of the temperature distribution Mik(r) steady state stress caused by the temperature field Mok(r) metal temperature
M. Guyette, M. De Smet / Nuclear Engineering and Design 153 (1995) 183-195 t
Tk(r, t)
Tfk(t)
Uk(r, t) v
V
wi
current time vector containing the mean temperature in each measurement section coefficient of the kth harmonic in the Fourier development of the metal temperature component of harmonic k of the fluid temperature inverse transfer function of harmonic k vector containing response of interest matrix containing response in the points of interest due to unit loads at the measurement sections weights for the calculation of the fluid temperature
Greek letters vector containing the free curvature at each measurement section
0 K 2 ak(r, t)
q51k
195
vector containing the end displacements in each coordinate direction current azimuthal coordinate metal thermal diffusivity metal thermal conductivity sampling time interval coefficient of harmonic k of the stress distribution (represent any stress component) eigenvalue of index l for harmonic k
References B.A. Boley and J.H. Weiner, Theory of Thermal Stresses, Krieger, Malabar, FL, 1985. H.S. Carslaw and J.C. Jaeger, Conduction of Heat in Solids, Clarendon, Oxford, 1959. M. Guyette, Prediction of fluid temperature from measurements of outside wall temperatures in pipes, ASME J. Pressure Vessel Technol., 116 (1994) 179 187. S.P. Timoshenko and J.N. Goodier, Theory of Elasticity, McGraw-Hill, Kogakusha, Tokyo, 1970.