- Email: [email protected]
An advanced method for determination of loss of coolant accident in nuclear power plants
Nuclear Engineering and Design 241 (2011) 2013–2019
Contents lists available at ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
An advanced method for determination of loss of coolant accident in nuclear power plants R. Mahmoodi, M. Shahriari ∗ , A. Zolfaghari, A. Minuchehr Department of Engineering, Shahid Beheshti University, GC, Evin, Tehran, IR, Iran
a r t i c l e
i n f o
Article history: Received 24 March 2010 Received in revised form 17 December 2010 Accepted 9 January 2011
a b s t r a c t A major objective in reactor design is to provide the capability to withstand a wide range of postulated events without exceeding specified safety limits. Assessment of the consequence of hypothetical loss of coolant accident (LOCA) in primary circuit is an essential element to address fulfilment of acceptance criteria. In addition, finding the position of rupture, one could manage accident in a right direction. In this work, the transient vibration signal from a pipe rupture is used to determine the position of LOCA. A finite element formulation (Galerkin Method) is implemented to include the effect of fluid-structure interaction (FSI). The coupled equations of fluid motion and pipe displacement are solved. The obtained results are in good agreement with published data. Fast Fourier transform (FFT) provides an alternate way of representing data. Instead of representing vibration signal amplitude as a function of time, the signal is represented by the amount of information, which is contained at different frequencies. The most of frequencies of structure and fluid coupled are presented in the FFT of structural response and through it the dominant frequency of excitation is obtained. Furthermore, the power spectral density (PSD), a measurement of energy at various frequencies is worked out. MATLAB software is used to convert signals from the time to frequency domain and to obtain PSD of signals. © 2011 Published by Elsevier B.V.
1. Introduction The loss of coolant accident (LOCA) is one of the incidents that occurred in nuclear power plants possibly that if uncontrolled can be disastrous. Analysis of such LOCA along other phenomena which are occurred in a nuclear plant through using safety systems are issue of scientific and practical research in the field of nuclear energy. When in a fluid-filled piping system under stationary conditions at some location the fluid velocity changes, pressure surges are induced and propagate at speed of sound through the system. Spectacular pressure surges, referred to as water hammer, which can be caused by rapid closing or opening of valve, stopping or starting of pump, pipe rupture, loud rejection of turbine, seismic excitation and so forth. Tijsseling (1996) carried out a very detailed review of transient phenomena in liquid-filled pipe systems. He dealt with water hammer, cavitations, structural dynamics and fluid-structure inter-
∗ Corresponding author at: Engineering Department of Shahid Beheshti University, GC. Evin, Tehran, IR, Post Box: 1985963113, Iran. Tel.: +98 21 22431596; fax: +98 21 22431780. E-mail address: [email protected] (M. Shahriari). 0029-5493/$ – see front matter © 2011 Published by Elsevier B.V. doi:10.1016/j.nucengdes.2011.01.009
action (FSI). His main focus was on the history of FSI research in time domain. In conventional water hammer analyses, pipe elasticity is incorporated in the propagation speed of the pressure waves. Contrary to conventional water hammer treatment analysis and uncoupled analysis, the effects of FSI are problem dependent. Treatment including fluid-structure interaction may lead to higher or lower extreme pressure and stresses, change in natural frequencies of the system and more damping and dispersion in the pressure and stress history in comparison to analysis, which the coupling is not considered. Wiggert et al. (1987) used the MOC to study transients in pipeline systems. They identified seven wave components, coupled axial compression of liquid and pipe material, along transverse shear and bending of the pipe elements in two principal direction and torsion of the pipe wall. Heinsbroek (1996) reported an application of FSI in the nuclear industry. His analysis was based on a combination of MOC and FEM. He showed the superiority of MOC for axial dynamics. In this paper, finite elements analysis of the coupled system is carried out through using Galerkin approach, which leads to a matrix problem. Using finite element technique for solving structural vibration equations has more flexibility. Information about vibration of structure is contained in the peak amplitude and fre-
2014
R. Mahmoodi et al. / Nuclear Engineering and Design 241 (2011) 2013–2019
The behavior of fluid is governed by extended equations of momentum and mass conservation while, the dynamics of the pipe is modeled by standard beam theory. Water hammer theory has been presented by several authors from the 19th century onwards (e.g. Tijsseling and Lavooij (1990) and among them, Joukowsky (1990) conducted a systematic study of the water distribution system in Moscow and expressed: P = cV
Fig. 1. Hermite interpolation functions.
quency content of the initial wave and frequency content of each signal is considered in dominant frequency of signal instead of using the fluid and structure frequencies. The dominant frequency of combined system, coupled fluid and structure, is obtained by fast Fourier transfer as a method of signal processing. Li et al. (2002) used frequency domain analysis of FSI in liquidfilled pipe systems by transfer matrix method, along taking the Laplace transform. They implemented both harmonic and frequency response analysis. The methods are useful for obtaining the solution at any point of system. They also showed how the frequencies and mode shapes change with the radius and material properties of the piping system. The main goal of the present work is to describe a new method for determination of small or medium LOCA in a nuclear power plant. The method is based on frequency analysis of mechanical vibration signals due to water hammer phenomenon, which is generated during LOCA. Furthermore it is also utilized to find the position of pipe rupture following LOCA accident. 2. Governing equations for fluid and structure The transient behavior of fluid-filled pipeline system is governed by acoustic waves propagating in both fluid and pipe. Propagation of pressure waves in fluid, water hammer causes axial, lateral and torsional stress waves in the pipe wall. These waves may have very steep fronts when are generated by, for example, closure of lightly damped check valves, collapse of liquid column separations, cavitations or other phenomena. Steep-fronted waves are prone to excite the structural system and result pipes vibration. Pipe motion generates water hammer and invokes fluidstructure interaction. The numerical simulation of water-hammer events including fluid-structure interaction is usually based on onedimensional mathematical models.
Fig. 2. Tijsseling and Lavooij (1990) geometry and parameters.
This formula relates pressure fluctuations P to velocity changes V by a constant factor c, where is the mass density of the fluid and c the velocity of sound in the fluid. The expression for c was first derived to predict standing waves in musical instruments and pulsatile flows in blood vessels. For a compressible fluid in an elastic tube, c depends to the bulk elastic modulus of the fluid, the elastic modulus of the pipe, the inner radius of the pipe and to wall thickness. The water hammer equations are another version of the compressible fluid flow equations. The choice of version is problem-dependent. In basic water hammer approach friction and damping mechanisms are neglected and in classic water hammer fluid wall friction is taken into account. Finally extended water hammer allows pipe motion. The set of pipe dynamic equations suggested by Wiggert et al. (1987) can be written as: EAp u − mu¨ + 2AP = 0
(1)
w V˙ + P = 0
(2)
∗2
p˙ + w a V = 0
(3)
where a∗2 =
Kf /w , (1 + Kf D)/E ∗ t
E∗ =
E 1 − 2
where Kf is the fluid bulk modulus, mp , E, , Ap , D, w , t, u, P and V, are the mass per unit length of pipe, Young’s modulus of elasticity, Poisson’s ratio, the cross-sectional area, the inner diameter, density of fluid, the thickness of the pipeline, displacement of the pipe in x-direction, pressure and velocity of flow, respectively. The first term at the left hand side of Eq. (1) represents the axial stiffness of the pipe. The second term accounts for the axial inertia of the pipe. The second and third equations describe the compressible fluid flow equations. Taking derivative of Eq. (2) with respect to the axial direction and Eq. (3) with respect to time, respectively, one could eliminate one of the variables to obtain two wave equations, which are either in terms of pressure or in terms of velocity. In the present study the equations involving with velocity are used. 2.1. Computational method The wave equations obtained are elliptical in nature and suitable for FEM. The finite element method can be regarded as a generalization of the finite difference method. The regular mesh of the latter method is replaced generally in the former method by an arbitrary mesh. Several authors have illustrated this generalization by many examples. Finite difference methods for differential equations are based usually on a Taylor power series expansion for a function. The derivatives of a function at a point are expressed approximately in terms of the values of the function at nodes of the mesh, which are in the neighborhood of the point. The use of a local polynomial approximation is usual for finite elements method, but more general forms for a local approximation are sometimes used. Examples are given for problems in variety of fields. In finite elements consideration has to be given to the size, shape and interface continuity conditions for the elements. The completeness of the local polynomial for the element is important for the selection of an element.
R. Mahmoodi et al. / Nuclear Engineering and Design 241 (2011) 2013–2019
2015
Fig. 3. Pressure history near the valve from present study.
The elements must be capable of representing closely systems of complex shape. For two dimensions, arbitrary triangular and rectangular elements are in common use. Curve boundaries can be approximated closely with triangles. In three dimensions, arbitrary tetrahedral and triangular prisms are used. The former permits a very faithful representation of the three-dimensional system. It is a good idea to utilize finite element method as a solving procedure due to the structural vibration equations with standard beam theory. Since the boundary condition for the valve opening or closing event is in terms of flow velocity, it is preferred to use the wave equation in terms of flow velocity, which is given by: 1 ∂2 V ∂3 u ∂2 V − 2 − 2 2 = 0 2 2 a∗ ∂t ∂x ∂t ∂x
(4)
The relation between pressure and velocity given by Eq. (3) is used to obtain the pressure from velocity. On using the fact that effect of pipe displacement on the fluid is negligible, one can ignore the third term of Eq. (4) to reduce it to a simple form. The finite element form of wave equation, Eq. (4), and pressure equation, Eq. (3), using Galerkin technique is obtained as: [G]{V¨ } + [H]{V } = {0}
(5)
˙ + [B]{V } = {0} [A]{p}
(6)
For modeling the pipe, 3D beam element with six degree-offreedom per node is used. With the use of Eq. (1), one obtains: ¨ + [K]{u} − [S]{P} = 0 [M]{u}
(e)
[G] = [H] =
1 a∗2
l
l
NVT NV d 0
N V NV d T
0 l
NPT NP d
[A] =
0
[B] = w a∗2
[S] = 2
l
l
NPT NV d
0
NST NP d
0
where NS , NP and NV are shape functions of structure of pipe, pressure, and velocity of fluid, respectively. The Runge–Kutta fourth order integration scheme is used to evaluate the transient response for the valve opening or closure event. The fully coupled equations are given by:
⎡ ⎤
⎡
0 u˙ ⎢ u¨ ⎥ ⎢ −[M]−1 [K] ⎢˙ ⎥ ⎢ 0 ⎢V ⎥ = ⎢ ⎣ V¨ ⎦ ⎣ 0 P˙ 0
[I] 0 0 0 0
0 0 0 −[G]−1 [H] −[A]−1 [B]
⎤
⎧u⎫ 0 0 ⎪ ⎪ ⎨ u˙ ⎪ ⎬ 0 [M]−1 [S] ⎥ ⎪ ⎥ [I] 0 ⎥ V ⎪ ⎦⎪ ⎪ 0 0 ⎩ V˙ ⎪ ⎭ P 0 0 (8)
(7)
where [M] and [K] are the mass and stiffness matrix of the pipe and [S] is a vector which represents the interaction of pressure with the structure due to the Poisson coupling. [G], [H], [A], [B] and [S] are obtained through using Galerkin method along implementing Hermity equations as shape functions. The shape functions for elements are defined as: u(e) () = HN0 ()uN + HN1 ()
Using above functions yield the matrix coefficient of Eqs. (5)–(7) as:
∂u ∂
(e)
+ HN2 ()
N
∂2 u ∂ 2
(e)
The huge size of matrix causes an overflow error and even reducing spatial meshes or increasing the time steps cannot improve the required memory in PC computers. To overcome this difficulty, the modal reduction technique is implemented to decrease the size of both structural and fluid matrices (Jayaraj et al., 2003). According to modal reduction mode shapes are used with respect to natural frequencies. 2.2. Natural frequencies and mode shapes
N
which for a node (e refers the nodes) they are: H10 = 1 − 10 3 + 15 4 − 6 5 H20 = 10 3 − 15 4 + 6 5 H11 = − 6 3 + 8 4 − 3 5 H21 = −4 3 + 7 4 − 3 5 H12 = 0.5( 2 − 3 3 + 3 4 − 5 ) H22 = 0.5( 3 − 2 4 + 5 ) where HN0 (), HN1 (), and HN2 () are shape functions and illustrated in Fig. 1.
It is important to use dynamic displacement of structure equation along fluid flow according to FSI to find natural frequencies and mode shapes accurately. Paidoussis and Li (1993), reported very details of dynamics response of pipe vibrations with fluid content. They expressed the equation as: EI
∂2 u ∂2 u ∂2 u ∂4 u − (Ap + Ai ) 2 = 0 + w Ai U 2 2 + 2w Ai U ∂x∂t ∂x4 ∂x ∂t
(9)
where I, , Ai , and U are the cross-sectional moment of inertia, the density of pipe, and the velocity of fluid, respectively. Implementing Galerkin method, one obtains M u¨ + Cc u˙ + Ku = 0
2016
R. Mahmoodi et al. / Nuclear Engineering and Design 241 (2011) 2013–2019
Fig. 4. The velocity of fluid near the valve from present study.
where Cc is the accelerating Crevlious matrix. The equations can be cast:
X˙
=
[0]n×n −M −1 K
[I]n×n −M −1 Cc
X
where,
X
J=
=
u u˙
[0]n×n −M −1 K
[I]n×n −M −1 Cc
Fig. 5. Time domain plots of the displacement of pipe.
(10)
R. Mahmoodi et al. / Nuclear Engineering and Design 241 (2011) 2013–2019
Fig. 6. The maximum displacement of pipe.
2017
Fig. 8. Mode shape of pipe system in domain frequency (109 Hz).
The eigenvalues and eigenvectors of J matrix are the structural natural frequencies and mode shapes of pipe with constant fluid flow, respectively. The first few mode shapes of the structure [ϕs ] as well as the fluid [ϕf ] are used for transforming the variables by substituting: {u} = [ϕS ]{x}
(11)
{V } = [ϕf ]{Vm }
(12)
The mode shapes of structure and fluid are obtained using dynamic vibration of equation with constant flow fluid interaction. The mode shapes corresponding to fluid and structural frequencies with the same range are used. When the frequencies of the fluid are
much higher than those of the fundamental frequency of the structure, one has to include a few mode shapes of the structure having frequencies in the range of fluid frequencies. After substitution and multiplying through by [ϕs ]T and [ϕf ]T respectively, one gets: T
T
T
T
T
[f ] [G][f ]{V¨ m } + [f ] [H][f ]{Vm } = 0
(14)
˙ + [B][f ]{Vm } = 0 [A]{P}
(15)
[S ] [M][S ]{¨x} + [S ] [K][S ]{x} − [S ] [S2 ]{p} = 0
(13)
The values of the structural variables as well as of the fluid flow velocity in the modal coordinates are transformed to the nodal coordinates by multiplying to [ϕS ] and [ϕf ]
Fig. 7. Frequency domain of vibration signals with FFT method.
2018
R. Mahmoodi et al. / Nuclear Engineering and Design 241 (2011) 2013–2019
Fig. 9. The power spectral density (PSD) at different values of l.
respectively. The pressure values can be used directly without transformation.
For investigation of the modal reduction, the fluid frequencies are evaluated from the finite element form of Eq. (5) and Fig. 4 shows the velocity of fluid near to the point.
2.3. Assumptions 3.1. Initial and boundary conditions The formulation is based on using of linear elasticity, no buckling, negligible radial inertia, low Mach number. Furthermore, it is understood that instantaneous pressure will remain above vapor pressure so that no cavitations will occur. 3. Validation study Tijsseling and Lavooij (1990) used the water hammer theory for the fluid coupled with beam theory for pipe to model FSI problems in non-rigid pipeline systems. They used simple model, which represented by four equations. The system analyzed consists of a pipe with length of 20 m. The other parameters of the reservoir, pipe and valve system have been illustrated in Fig. 2. In a steady state situation fluid is flowing from reservoir to toward a valve. The structural boundary conditions for the pipeline system were taken with no displacements at the valve as well as at the upstream reservoir end. At t = 0 ms, the valve is closed instantaneously. The pressure rise propagates with a sound speed toward the reservoir. Fig. 3 shows the pressure history near the value due to valve closure and furthermore illustrates a comparison of the results from the four model formulation of Tijsseling and Lavooij (1990). Table 1 shows the fluid and structural frequencies of this geometry.
For complete mathematical description of the problem, determination of initial and boundary conditions are essential. Denoting
Table 1 Structural and fluid frequencies, in Hz, for Tijsseling and Lavooij (1990), geometry. Serial no.
Structural frequency
Fluid frequency
1 2 3 4 5
43.96 113.76 221.67 365.32 544.8
167.28 334.43 667.8 833.77 999
Table 2 Frequencies in Hz of pipe system and their amplitude at l = 20 cm. Serial no.
Frequency (Hz)
Amplitude (m)
1 2 3 4 5 6
190 112 115 120 248 252
2.33E-5 1.532E-5 1.656E-5 6.941E-6 1.386E-6 6.19E-6
R. Mahmoodi et al. / Nuclear Engineering and Design 241 (2011) 2013–2019
to Fig. 2 and giving attention to the restrain points in the pipe ends, one may specify the structural boundary conditions. In the steady state condition fluid is flowing from reservoir to valve and the pressure and velocity of fluid boundary conditions are known in advance. Furthermore, initial conditions are taken from steady state situation of system too, but it is assumed that valve acting instantaneously and after 0.1 s the velocity and pressure of fluid retain steady state situation linearity.
2019
4.2. Power spectral density It is important to find severity of the pressure surge impacts that clarify the intensity of LOCA. The power spectral density (PSD), a measurement of the energy at various frequencies, is computed. Fig. 9 shows the PSD of the vibration signals at the different position of pipe length. It is obvious that the energy of the vibration signals reduce along increasing length of pipe, which has been shown in Fig. 9.
4. Processing analysis 5. Conclusion 4.1. Frequency domain analysis In this work a software with a graphically interface under the MATLAB environment, was developed to illustrate the signal processing techniques used. In order to validate the used technique, the outcome of our modeling is compared with the results of Tijsseling and Lavooij (1990). To prevent large amount of computation and analyzing data, we extended our test case problem for a pipe with 2 m length and the same parameters illustrated in the geometry. Fig. 5 shows time domain plot of pipe displacement for different values of l (the position of a rupture from a transient vibration signal). Clearly, as l increases (from the valve) the displacement of pipe decreases. Fig. 5 shows this decrease during the time vs. displacement of the pipe. FFTs and power spectrum are useful for measuring the frequency content of vibration transient signals. FFTs produce the average frequency content of signal over the entire time that a signal is acquired. For this reason, to use FFT for the vibration signals analysis, one needs to include the energy at each frequency line for power spectrum computation. In this simulation the sampling frequency and sampling rate are taken 500 kHz and 2.5 kHz, respectively. Furthermore, to verify the major frequency component of excitation, the FFT of the structural response including the effect of vibration of the fluid on the structure is shown in Fig. 6. In this case, some of frequencies are suppressed and some are slightly deviated from their original values. The dominant frequency of excitation of the system is 109 Hz. The frequencies of pipe system are dependent to the length of pipe. It is concluded from Fig. 6 that maximum displacement of the pipe decreases with increasing l along damping oscillation trend. Fig. 7 shows the amplitude of vibration signals in FFT plots is decreased. Fig. 8 also illustrates decreasing damped trend, which represents the mode shape of structure at the dominant frequency 109 Hz. The natural frequency of the structure and fluid are found separately without coupling the structure and fluid for this case as well, and it is shown in Table 1. The first sixth frequencies of pipe system are given in Table 2. The obtained amplitudes are at those frequencies and l = 20 cm.
It was obtained, theoretically, that in a single straight pipe with open–closed and fixed–fixed boundaries for liquid and pipe, respectively, FSI affects the induced system vibration. In order to treat vibration of a pipe, a finite element formulation for the fully coupled pipe dynamic theory was implemented and water hammer phenomenon which occurs due to sudden valve closure or opening modeled by using a new structural response based finite element formulation. The comparison of the outcome of the pressure propagation with a benchmark problem validates the present formulation. On using the above method, the rupture of pipe, which causes water hammer phenomenon and loss of coolant accident in nuclear power plants can be recognized. Signal processing methods, FFTs, are effective procedures to process the vibration signals following a LOCA for finding frequencies and mode shapes in pipeline systems .The dominant frequency of the system, coupled fluid and structure, can be obtained by fast Fourier transfer. The severity of each signal or in other word the energy of each signal can be obtained through the PSD analysis. Although we focused on using the proposed method for nuclear power plants but it can be use for all power plants. References Heinsbroek, A., 1996. Fluid-structure interaction in non-rigid pipeline system. J. Nucl. Eng. Des. 172, 123–135. Jayaraj, K., Ganesan, N., Padmanabham, C., 2003. Model reduction for parametric instability analysis of shells conveying fluid. J. Sound Vib. 262, 633–649. Joukowsky, N., 1990. Uber den hydraulischen stoss in wasserleitungsrohren. Meˇımoires de l’Acadeˇımie Impeˇıriale des Sciences de St. Peˇıtersbourg Series 8 (9), 5. Li, Q.S., Yang, K., Zhang, L., Zhang, N., 2002. Frequency domain analysis of fluid structure interaction in liquid-filled pipe systems by transfer matrix method. Int. J. Mech. Sci. 44, 2067–2087. Paidoussis, M.P., Li, G.X., 1993. Pipe conveying fluid: a model dynamical problem. J. Fluid Struct. 7, 137–204. Tijsseling, A.S., 1996. Fluid-structure interaction in liquid-filled pipe systems: a review. J. Fluids Struct. 10, 109–146. Tijsseling, A.S., Lavooij, C.S.W., 1990. Water hammer with fluid-structure interaction. J. Appl. Sci. Res. 47, 273–285. Wiggert, D.C., Hatfield, F.J., Struckenbruck, S., 1987. Analysis of liquid and structural transients in piping by the method of characteristics. ASME. J. Fluid Eng. 109, 161–165.
Contents lists available at ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
An advanced method for determination of loss of coolant accident in nuclear power plants R. Mahmoodi, M. Shahriari ∗ , A. Zolfaghari, A. Minuchehr Department of Engineering, Shahid Beheshti University, GC, Evin, Tehran, IR, Iran
a r t i c l e
i n f o
Article history: Received 24 March 2010 Received in revised form 17 December 2010 Accepted 9 January 2011
a b s t r a c t A major objective in reactor design is to provide the capability to withstand a wide range of postulated events without exceeding specified safety limits. Assessment of the consequence of hypothetical loss of coolant accident (LOCA) in primary circuit is an essential element to address fulfilment of acceptance criteria. In addition, finding the position of rupture, one could manage accident in a right direction. In this work, the transient vibration signal from a pipe rupture is used to determine the position of LOCA. A finite element formulation (Galerkin Method) is implemented to include the effect of fluid-structure interaction (FSI). The coupled equations of fluid motion and pipe displacement are solved. The obtained results are in good agreement with published data. Fast Fourier transform (FFT) provides an alternate way of representing data. Instead of representing vibration signal amplitude as a function of time, the signal is represented by the amount of information, which is contained at different frequencies. The most of frequencies of structure and fluid coupled are presented in the FFT of structural response and through it the dominant frequency of excitation is obtained. Furthermore, the power spectral density (PSD), a measurement of energy at various frequencies is worked out. MATLAB software is used to convert signals from the time to frequency domain and to obtain PSD of signals. © 2011 Published by Elsevier B.V.
1. Introduction The loss of coolant accident (LOCA) is one of the incidents that occurred in nuclear power plants possibly that if uncontrolled can be disastrous. Analysis of such LOCA along other phenomena which are occurred in a nuclear plant through using safety systems are issue of scientific and practical research in the field of nuclear energy. When in a fluid-filled piping system under stationary conditions at some location the fluid velocity changes, pressure surges are induced and propagate at speed of sound through the system. Spectacular pressure surges, referred to as water hammer, which can be caused by rapid closing or opening of valve, stopping or starting of pump, pipe rupture, loud rejection of turbine, seismic excitation and so forth. Tijsseling (1996) carried out a very detailed review of transient phenomena in liquid-filled pipe systems. He dealt with water hammer, cavitations, structural dynamics and fluid-structure inter-
∗ Corresponding author at: Engineering Department of Shahid Beheshti University, GC. Evin, Tehran, IR, Post Box: 1985963113, Iran. Tel.: +98 21 22431596; fax: +98 21 22431780. E-mail address: [email protected] (M. Shahriari). 0029-5493/$ – see front matter © 2011 Published by Elsevier B.V. doi:10.1016/j.nucengdes.2011.01.009
action (FSI). His main focus was on the history of FSI research in time domain. In conventional water hammer analyses, pipe elasticity is incorporated in the propagation speed of the pressure waves. Contrary to conventional water hammer treatment analysis and uncoupled analysis, the effects of FSI are problem dependent. Treatment including fluid-structure interaction may lead to higher or lower extreme pressure and stresses, change in natural frequencies of the system and more damping and dispersion in the pressure and stress history in comparison to analysis, which the coupling is not considered. Wiggert et al. (1987) used the MOC to study transients in pipeline systems. They identified seven wave components, coupled axial compression of liquid and pipe material, along transverse shear and bending of the pipe elements in two principal direction and torsion of the pipe wall. Heinsbroek (1996) reported an application of FSI in the nuclear industry. His analysis was based on a combination of MOC and FEM. He showed the superiority of MOC for axial dynamics. In this paper, finite elements analysis of the coupled system is carried out through using Galerkin approach, which leads to a matrix problem. Using finite element technique for solving structural vibration equations has more flexibility. Information about vibration of structure is contained in the peak amplitude and fre-
2014
R. Mahmoodi et al. / Nuclear Engineering and Design 241 (2011) 2013–2019
The behavior of fluid is governed by extended equations of momentum and mass conservation while, the dynamics of the pipe is modeled by standard beam theory. Water hammer theory has been presented by several authors from the 19th century onwards (e.g. Tijsseling and Lavooij (1990) and among them, Joukowsky (1990) conducted a systematic study of the water distribution system in Moscow and expressed: P = cV
Fig. 1. Hermite interpolation functions.
quency content of the initial wave and frequency content of each signal is considered in dominant frequency of signal instead of using the fluid and structure frequencies. The dominant frequency of combined system, coupled fluid and structure, is obtained by fast Fourier transfer as a method of signal processing. Li et al. (2002) used frequency domain analysis of FSI in liquidfilled pipe systems by transfer matrix method, along taking the Laplace transform. They implemented both harmonic and frequency response analysis. The methods are useful for obtaining the solution at any point of system. They also showed how the frequencies and mode shapes change with the radius and material properties of the piping system. The main goal of the present work is to describe a new method for determination of small or medium LOCA in a nuclear power plant. The method is based on frequency analysis of mechanical vibration signals due to water hammer phenomenon, which is generated during LOCA. Furthermore it is also utilized to find the position of pipe rupture following LOCA accident. 2. Governing equations for fluid and structure The transient behavior of fluid-filled pipeline system is governed by acoustic waves propagating in both fluid and pipe. Propagation of pressure waves in fluid, water hammer causes axial, lateral and torsional stress waves in the pipe wall. These waves may have very steep fronts when are generated by, for example, closure of lightly damped check valves, collapse of liquid column separations, cavitations or other phenomena. Steep-fronted waves are prone to excite the structural system and result pipes vibration. Pipe motion generates water hammer and invokes fluidstructure interaction. The numerical simulation of water-hammer events including fluid-structure interaction is usually based on onedimensional mathematical models.
Fig. 2. Tijsseling and Lavooij (1990) geometry and parameters.
This formula relates pressure fluctuations P to velocity changes V by a constant factor c, where is the mass density of the fluid and c the velocity of sound in the fluid. The expression for c was first derived to predict standing waves in musical instruments and pulsatile flows in blood vessels. For a compressible fluid in an elastic tube, c depends to the bulk elastic modulus of the fluid, the elastic modulus of the pipe, the inner radius of the pipe and to wall thickness. The water hammer equations are another version of the compressible fluid flow equations. The choice of version is problem-dependent. In basic water hammer approach friction and damping mechanisms are neglected and in classic water hammer fluid wall friction is taken into account. Finally extended water hammer allows pipe motion. The set of pipe dynamic equations suggested by Wiggert et al. (1987) can be written as: EAp u − mu¨ + 2AP = 0
(1)
w V˙ + P = 0
(2)
∗2
p˙ + w a V = 0
(3)
where a∗2 =
Kf /w , (1 + Kf D)/E ∗ t
E∗ =
E 1 − 2
where Kf is the fluid bulk modulus, mp , E, , Ap , D, w , t, u, P and V, are the mass per unit length of pipe, Young’s modulus of elasticity, Poisson’s ratio, the cross-sectional area, the inner diameter, density of fluid, the thickness of the pipeline, displacement of the pipe in x-direction, pressure and velocity of flow, respectively. The first term at the left hand side of Eq. (1) represents the axial stiffness of the pipe. The second term accounts for the axial inertia of the pipe. The second and third equations describe the compressible fluid flow equations. Taking derivative of Eq. (2) with respect to the axial direction and Eq. (3) with respect to time, respectively, one could eliminate one of the variables to obtain two wave equations, which are either in terms of pressure or in terms of velocity. In the present study the equations involving with velocity are used. 2.1. Computational method The wave equations obtained are elliptical in nature and suitable for FEM. The finite element method can be regarded as a generalization of the finite difference method. The regular mesh of the latter method is replaced generally in the former method by an arbitrary mesh. Several authors have illustrated this generalization by many examples. Finite difference methods for differential equations are based usually on a Taylor power series expansion for a function. The derivatives of a function at a point are expressed approximately in terms of the values of the function at nodes of the mesh, which are in the neighborhood of the point. The use of a local polynomial approximation is usual for finite elements method, but more general forms for a local approximation are sometimes used. Examples are given for problems in variety of fields. In finite elements consideration has to be given to the size, shape and interface continuity conditions for the elements. The completeness of the local polynomial for the element is important for the selection of an element.
R. Mahmoodi et al. / Nuclear Engineering and Design 241 (2011) 2013–2019
2015
Fig. 3. Pressure history near the valve from present study.
The elements must be capable of representing closely systems of complex shape. For two dimensions, arbitrary triangular and rectangular elements are in common use. Curve boundaries can be approximated closely with triangles. In three dimensions, arbitrary tetrahedral and triangular prisms are used. The former permits a very faithful representation of the three-dimensional system. It is a good idea to utilize finite element method as a solving procedure due to the structural vibration equations with standard beam theory. Since the boundary condition for the valve opening or closing event is in terms of flow velocity, it is preferred to use the wave equation in terms of flow velocity, which is given by: 1 ∂2 V ∂3 u ∂2 V − 2 − 2 2 = 0 2 2 a∗ ∂t ∂x ∂t ∂x
(4)
The relation between pressure and velocity given by Eq. (3) is used to obtain the pressure from velocity. On using the fact that effect of pipe displacement on the fluid is negligible, one can ignore the third term of Eq. (4) to reduce it to a simple form. The finite element form of wave equation, Eq. (4), and pressure equation, Eq. (3), using Galerkin technique is obtained as: [G]{V¨ } + [H]{V } = {0}
(5)
˙ + [B]{V } = {0} [A]{p}
(6)
For modeling the pipe, 3D beam element with six degree-offreedom per node is used. With the use of Eq. (1), one obtains: ¨ + [K]{u} − [S]{P} = 0 [M]{u}
(e)
[G] = [H] =
1 a∗2
l
l
NVT NV d 0
N V NV d T
0 l
NPT NP d
[A] =
0
[B] = w a∗2
[S] = 2
l
l
NPT NV d
0
NST NP d
0
where NS , NP and NV are shape functions of structure of pipe, pressure, and velocity of fluid, respectively. The Runge–Kutta fourth order integration scheme is used to evaluate the transient response for the valve opening or closure event. The fully coupled equations are given by:
⎡ ⎤
⎡
0 u˙ ⎢ u¨ ⎥ ⎢ −[M]−1 [K] ⎢˙ ⎥ ⎢ 0 ⎢V ⎥ = ⎢ ⎣ V¨ ⎦ ⎣ 0 P˙ 0
[I] 0 0 0 0
0 0 0 −[G]−1 [H] −[A]−1 [B]
⎤
⎧u⎫ 0 0 ⎪ ⎪ ⎨ u˙ ⎪ ⎬ 0 [M]−1 [S] ⎥ ⎪ ⎥ [I] 0 ⎥ V ⎪ ⎦⎪ ⎪ 0 0 ⎩ V˙ ⎪ ⎭ P 0 0 (8)
(7)
where [M] and [K] are the mass and stiffness matrix of the pipe and [S] is a vector which represents the interaction of pressure with the structure due to the Poisson coupling. [G], [H], [A], [B] and [S] are obtained through using Galerkin method along implementing Hermity equations as shape functions. The shape functions for elements are defined as: u(e) () = HN0 ()uN + HN1 ()
Using above functions yield the matrix coefficient of Eqs. (5)–(7) as:
∂u ∂
(e)
+ HN2 ()
N
∂2 u ∂ 2
(e)
The huge size of matrix causes an overflow error and even reducing spatial meshes or increasing the time steps cannot improve the required memory in PC computers. To overcome this difficulty, the modal reduction technique is implemented to decrease the size of both structural and fluid matrices (Jayaraj et al., 2003). According to modal reduction mode shapes are used with respect to natural frequencies. 2.2. Natural frequencies and mode shapes
N
which for a node (e refers the nodes) they are: H10 = 1 − 10 3 + 15 4 − 6 5 H20 = 10 3 − 15 4 + 6 5 H11 = − 6 3 + 8 4 − 3 5 H21 = −4 3 + 7 4 − 3 5 H12 = 0.5( 2 − 3 3 + 3 4 − 5 ) H22 = 0.5( 3 − 2 4 + 5 ) where HN0 (), HN1 (), and HN2 () are shape functions and illustrated in Fig. 1.
It is important to use dynamic displacement of structure equation along fluid flow according to FSI to find natural frequencies and mode shapes accurately. Paidoussis and Li (1993), reported very details of dynamics response of pipe vibrations with fluid content. They expressed the equation as: EI
∂2 u ∂2 u ∂2 u ∂4 u − (Ap + Ai ) 2 = 0 + w Ai U 2 2 + 2w Ai U ∂x∂t ∂x4 ∂x ∂t
(9)
where I, , Ai , and U are the cross-sectional moment of inertia, the density of pipe, and the velocity of fluid, respectively. Implementing Galerkin method, one obtains M u¨ + Cc u˙ + Ku = 0
2016
R. Mahmoodi et al. / Nuclear Engineering and Design 241 (2011) 2013–2019
Fig. 4. The velocity of fluid near the valve from present study.
where Cc is the accelerating Crevlious matrix. The equations can be cast:
X˙
=
[0]n×n −M −1 K
[I]n×n −M −1 Cc
X
where,
X
J=
=
u u˙
[0]n×n −M −1 K
[I]n×n −M −1 Cc
Fig. 5. Time domain plots of the displacement of pipe.
(10)
R. Mahmoodi et al. / Nuclear Engineering and Design 241 (2011) 2013–2019
Fig. 6. The maximum displacement of pipe.
2017
Fig. 8. Mode shape of pipe system in domain frequency (109 Hz).
The eigenvalues and eigenvectors of J matrix are the structural natural frequencies and mode shapes of pipe with constant fluid flow, respectively. The first few mode shapes of the structure [ϕs ] as well as the fluid [ϕf ] are used for transforming the variables by substituting: {u} = [ϕS ]{x}
(11)
{V } = [ϕf ]{Vm }
(12)
The mode shapes of structure and fluid are obtained using dynamic vibration of equation with constant flow fluid interaction. The mode shapes corresponding to fluid and structural frequencies with the same range are used. When the frequencies of the fluid are
much higher than those of the fundamental frequency of the structure, one has to include a few mode shapes of the structure having frequencies in the range of fluid frequencies. After substitution and multiplying through by [ϕs ]T and [ϕf ]T respectively, one gets: T
T
T
T
T
[f ] [G][f ]{V¨ m } + [f ] [H][f ]{Vm } = 0
(14)
˙ + [B][f ]{Vm } = 0 [A]{P}
(15)
[S ] [M][S ]{¨x} + [S ] [K][S ]{x} − [S ] [S2 ]{p} = 0
(13)
The values of the structural variables as well as of the fluid flow velocity in the modal coordinates are transformed to the nodal coordinates by multiplying to [ϕS ] and [ϕf ]
Fig. 7. Frequency domain of vibration signals with FFT method.
2018
R. Mahmoodi et al. / Nuclear Engineering and Design 241 (2011) 2013–2019
Fig. 9. The power spectral density (PSD) at different values of l.
respectively. The pressure values can be used directly without transformation.
For investigation of the modal reduction, the fluid frequencies are evaluated from the finite element form of Eq. (5) and Fig. 4 shows the velocity of fluid near to the point.
2.3. Assumptions 3.1. Initial and boundary conditions The formulation is based on using of linear elasticity, no buckling, negligible radial inertia, low Mach number. Furthermore, it is understood that instantaneous pressure will remain above vapor pressure so that no cavitations will occur. 3. Validation study Tijsseling and Lavooij (1990) used the water hammer theory for the fluid coupled with beam theory for pipe to model FSI problems in non-rigid pipeline systems. They used simple model, which represented by four equations. The system analyzed consists of a pipe with length of 20 m. The other parameters of the reservoir, pipe and valve system have been illustrated in Fig. 2. In a steady state situation fluid is flowing from reservoir to toward a valve. The structural boundary conditions for the pipeline system were taken with no displacements at the valve as well as at the upstream reservoir end. At t = 0 ms, the valve is closed instantaneously. The pressure rise propagates with a sound speed toward the reservoir. Fig. 3 shows the pressure history near the value due to valve closure and furthermore illustrates a comparison of the results from the four model formulation of Tijsseling and Lavooij (1990). Table 1 shows the fluid and structural frequencies of this geometry.
For complete mathematical description of the problem, determination of initial and boundary conditions are essential. Denoting
Table 1 Structural and fluid frequencies, in Hz, for Tijsseling and Lavooij (1990), geometry. Serial no.
Structural frequency
Fluid frequency
1 2 3 4 5
43.96 113.76 221.67 365.32 544.8
167.28 334.43 667.8 833.77 999
Table 2 Frequencies in Hz of pipe system and their amplitude at l = 20 cm. Serial no.
Frequency (Hz)
Amplitude (m)
1 2 3 4 5 6
190 112 115 120 248 252
2.33E-5 1.532E-5 1.656E-5 6.941E-6 1.386E-6 6.19E-6
R. Mahmoodi et al. / Nuclear Engineering and Design 241 (2011) 2013–2019
to Fig. 2 and giving attention to the restrain points in the pipe ends, one may specify the structural boundary conditions. In the steady state condition fluid is flowing from reservoir to valve and the pressure and velocity of fluid boundary conditions are known in advance. Furthermore, initial conditions are taken from steady state situation of system too, but it is assumed that valve acting instantaneously and after 0.1 s the velocity and pressure of fluid retain steady state situation linearity.
2019
4.2. Power spectral density It is important to find severity of the pressure surge impacts that clarify the intensity of LOCA. The power spectral density (PSD), a measurement of the energy at various frequencies, is computed. Fig. 9 shows the PSD of the vibration signals at the different position of pipe length. It is obvious that the energy of the vibration signals reduce along increasing length of pipe, which has been shown in Fig. 9.
4. Processing analysis 5. Conclusion 4.1. Frequency domain analysis In this work a software with a graphically interface under the MATLAB environment, was developed to illustrate the signal processing techniques used. In order to validate the used technique, the outcome of our modeling is compared with the results of Tijsseling and Lavooij (1990). To prevent large amount of computation and analyzing data, we extended our test case problem for a pipe with 2 m length and the same parameters illustrated in the geometry. Fig. 5 shows time domain plot of pipe displacement for different values of l (the position of a rupture from a transient vibration signal). Clearly, as l increases (from the valve) the displacement of pipe decreases. Fig. 5 shows this decrease during the time vs. displacement of the pipe. FFTs and power spectrum are useful for measuring the frequency content of vibration transient signals. FFTs produce the average frequency content of signal over the entire time that a signal is acquired. For this reason, to use FFT for the vibration signals analysis, one needs to include the energy at each frequency line for power spectrum computation. In this simulation the sampling frequency and sampling rate are taken 500 kHz and 2.5 kHz, respectively. Furthermore, to verify the major frequency component of excitation, the FFT of the structural response including the effect of vibration of the fluid on the structure is shown in Fig. 6. In this case, some of frequencies are suppressed and some are slightly deviated from their original values. The dominant frequency of excitation of the system is 109 Hz. The frequencies of pipe system are dependent to the length of pipe. It is concluded from Fig. 6 that maximum displacement of the pipe decreases with increasing l along damping oscillation trend. Fig. 7 shows the amplitude of vibration signals in FFT plots is decreased. Fig. 8 also illustrates decreasing damped trend, which represents the mode shape of structure at the dominant frequency 109 Hz. The natural frequency of the structure and fluid are found separately without coupling the structure and fluid for this case as well, and it is shown in Table 1. The first sixth frequencies of pipe system are given in Table 2. The obtained amplitudes are at those frequencies and l = 20 cm.
It was obtained, theoretically, that in a single straight pipe with open–closed and fixed–fixed boundaries for liquid and pipe, respectively, FSI affects the induced system vibration. In order to treat vibration of a pipe, a finite element formulation for the fully coupled pipe dynamic theory was implemented and water hammer phenomenon which occurs due to sudden valve closure or opening modeled by using a new structural response based finite element formulation. The comparison of the outcome of the pressure propagation with a benchmark problem validates the present formulation. On using the above method, the rupture of pipe, which causes water hammer phenomenon and loss of coolant accident in nuclear power plants can be recognized. Signal processing methods, FFTs, are effective procedures to process the vibration signals following a LOCA for finding frequencies and mode shapes in pipeline systems .The dominant frequency of the system, coupled fluid and structure, can be obtained by fast Fourier transfer. The severity of each signal or in other word the energy of each signal can be obtained through the PSD analysis. Although we focused on using the proposed method for nuclear power plants but it can be use for all power plants. References Heinsbroek, A., 1996. Fluid-structure interaction in non-rigid pipeline system. J. Nucl. Eng. Des. 172, 123–135. Jayaraj, K., Ganesan, N., Padmanabham, C., 2003. Model reduction for parametric instability analysis of shells conveying fluid. J. Sound Vib. 262, 633–649. Joukowsky, N., 1990. Uber den hydraulischen stoss in wasserleitungsrohren. Meˇımoires de l’Acadeˇımie Impeˇıriale des Sciences de St. Peˇıtersbourg Series 8 (9), 5. Li, Q.S., Yang, K., Zhang, L., Zhang, N., 2002. Frequency domain analysis of fluid structure interaction in liquid-filled pipe systems by transfer matrix method. Int. J. Mech. Sci. 44, 2067–2087. Paidoussis, M.P., Li, G.X., 1993. Pipe conveying fluid: a model dynamical problem. J. Fluid Struct. 7, 137–204. Tijsseling, A.S., 1996. Fluid-structure interaction in liquid-filled pipe systems: a review. J. Fluids Struct. 10, 109–146. Tijsseling, A.S., Lavooij, C.S.W., 1990. Water hammer with fluid-structure interaction. J. Appl. Sci. Res. 47, 273–285. Wiggert, D.C., Hatfield, F.J., Struckenbruck, S., 1987. Analysis of liquid and structural transients in piping by the method of characteristics. ASME. J. Fluid Eng. 109, 161–165.