At < — (191)
This is the so-called CFL condition for stability of the explicit difference scheme (not sufficient for all explicit difference schemes) and as Ä is prescribed as external parameter, Ax is defined by number of computational volumes and pipe length (desired accuracy), it follows that stability restriction suggests the maximal possible time step At. LeVeque [78] showed that the Courant-Friedrichs-Levy condition is a necessary and sufficient condition for the stability of the characteristic upwind difference scheme.
First order upwind scheme for systems with mixed sign characteristics. The characteristic partial differential Eq. (185) refers to the physical model (set of partial differential equations), where the sign of the characteristic velocities varies with position and time. It is necessary to perform the appropriate spatial differencing according to the sign of the characteristic speed Ä of the partial differential equation in order to obtain a useful one-sided scheme. An improvement of the upwind scheme with upwinding principle for systems of equations with characteristics of a mixed sign [49] enables appropriate splitting between the characteristics propagating to the left and to the right (superscripts - and +, respectively):
^^"- (^^^-^^- (A- 1/ 2ßrflS (192)
where A*=diag(A1+,...,ÄM+), l\.~=diagyA1-,...,AM J. Subscripts j, j+1 and j-1 denote the grid points of the
spatial discretisation defined in the middle of the each computational volume. Subscripts j+1/2 and j-1/2 denote the values in the midpoint of two computational volumes, Ax denotes the length of one computational volume, superscripts n and n+1 denote the time levels and At denotes the time step interval between time levels n and n+1. The appropriate splitting between positive and negative waves is given through the application of the correction factors fp:
Ap=\Ap\-fp and ^p=UJ-fp- (193)
- 53 -
where index p is running over M eigenvalues of the system. The correction factors fp read:
f+=max
A
and fp- =min
P %\ P ...... IA
v J
0,A
(194)
The matrix A+ assembles all characteristics with the positive sign that travel to the right (positive characteristics) and the matrix A assembles all characteristics with negative sign that travel to the left (negative characteristics). The full matrix of characteristics A is correlated to partial matrices A+ and A at position j-1/2 and time n through the relationship:
A^1/2= (A+ )n_1/2 + (A- );n_1/2 (195)
Second order Lax-Wendroff scheme. The upwind scheme is first-order accurate and introduces numerical diffusion, yielding poor accuracy and smeared results. The method can be improved by approximating derivatives with 2nd order differences [77,78]. The basic form of these correction terms is motivated by the standard second order accurate Lax-Wendroff method. The Lax-Wendroff method for a system of partial differential equations is based upon the Taylor series expansion [78, p.100]. The first three terms (only) on the right-hand side of the Taylor series expansion where the spatial derivatives are replaced by central finite difference approximations gives the Lax-Wendroff finite difference scheme:
ly"+WMA+ );_1/2KV^)^ + (^/ (196)
However, the second (and higher) order methods are also not acceptable for shock wave simulations because they yield oscillations near discontinuities.
5.3. High resolution finite difference schemes
The term “high resolution” applies to methods that are at least second order accurate on smooth solutions and yet give well resolved, non-oscillatory discontinuities. The idea behind is to use a high order method, but to modify the method in the neighborhood of discontinuities to the monotone first order method that behaves well near discontinuities. For precise theory and derivation of the high resolution method see LeVeque [77,78], Hirsch [57] or Toro [128]. Here, a brief overview is given.
Total variation diminishing schemes. Hirsch [57] showed that the total variation of the nonlinear conservation law in Eq. (186), that is defined as:
?d#w
TV = —5 cfx (197)
A1
is not increasing with time. The total variation of a discrete solution at time level n is defined as:
N
TV (L") = X ^l+1 - 41 (1 98)
j=1
Total variation diminishing schemes are numerical schemes for which the total variation of the numerical solutions is not increasing (LeVeque [77,78]). Finite difference scheme to be TVD:
šr1=Sj-cjMtj-Šj-J+D1/2{fa-$) (199)
- 54 -
is total variation diminishing if the following condition holds for each j:
Cj_1/2 > 0
and
D;+1/2 ^ 0
and
C/-1/2 +^y+1/2 - 1
(200)
The first order upwind scheme can be written in the form of Eq. (199):
«+1
,gn Ul + ^ A^ gn en \ , Ul ~ ^ Af ^ gn n
W W W-1J+ W+1 W
Ax
Ax
(201)
It is obvious that the upwind scheme is total variation diminishing. Analogously, the Lax-Wendroff scheme Eq. (196) written in the form of Eq. (199) yields:
C =^
^Äx + 1)^Af
Ax
en _ en <5j bj-
X&
Ax
1)a
Af
Ax
en en
W+1 W
(202) < 0 .
The Lax-Wendroff scheme is not total variation diminishing because the CFL condition gives D;+1/2
Hirsch [57] showed that:
• The linear total variation diminishing scheme can be only first order accurate.
• The total variation diminishing scheme for which the new value in pointy at time level n+1 is evaluated from three discrete points defined at time level n, can be only first order accurate.
Hirsch therefore showed that second order total variation diminishing schemes are only schemes where more than three points defined at time level n are nonlinearly applied. The possible solution is a combination of the first order upwind scheme and the Lax-Wendroff scheme:
n+1
(203)
where parameters are flux or slope limiters. The difference scheme in Eq. (203) is total variation
diminishing only if the slope limiters are appropriately defined to conform to LeVeque’s conditions. In addition, only an appropriate definition of the slope limiters gives a high resolution scheme with minimized numerical dissipation and minimized oscillations.
Slope limiters. The slope limiters define the share of the first and second order difference schemes. The slope limiters depend on smoothness and gradient of the characteristic variable E, at considered
point (volume). The slope limiter is close to 1 if the solution at time level n in the vicinity of the considered point is smooth and a larger part of the second order Lax-Wendroff difference scheme is applied. In vicinity of the discontinuous solutions, the slope limiter is close to 0 and a larger part of the first order upwind difference scheme is applied.
LeVeque [77,78], Hirsch [57] and Toro [128] gave several functions for the evaluation of the slope limiters and some of the most frequently used slope limiters are defined as:
0 = max(0,min(1,6>))
<» = (M+0)/(M+1)
0 = max(0,min((1 + 0)/2,2,20)) 0 = max(0,min(26>,1),min(6>,2))
Parameter 0 measures the smoothness of the characteristic variable E, near the considered point. The steepest solutions are obtained with the superbee limiter, while the most 'smeared' solutions but still
• Minmod:
• Van Leer:
MC: Superbee:
- 55 -
second-order accurate, are obtained with the minmod limiter. Solutions obtained with the van Leer and MC limiter are between (for example see Figures 27 and 28).
The slope limiter fixed to a constant value transforms the difference scheme in Eq. (203) into one of the following basic linear difference schemes:
• first order upwind scheme: v d At: The source terms are not stiff, the CFL time step At is applied, and the applied difference scheme is explicit given with Eq. (207).
• At > AtTiM > 0.1 At: The source terms are stiff (moderate), the time step AtTm is applied, and the applied difference scheme is explicit given with Eq. (207).
- 59 -
• AtTiM < 0.1 At: The source terms are stiff (significantly), the time step is smaller than 0.1 zlfand is defined manually according to previous numerical trials. An additional implicit iterations are applied.
Implicit iterations. For piping systems with very small cross-section (small moment of inertia lt), or for long piping systems (large Ax) the oscillation time step becomes much smaller than the time step defined by the CFL condition. Implicit iterations are introduced to avoid unreasonable small time steps and therefore long simulation times. For only one iteration, the minimal time step given with AtTm can be more than ten times larger. The best-on-test and therefore applied numerical procedure for implicit iterations due to stiff source terms is implicit iterative method with upwind averaged source terms. The explicit difference scheme given in Eq. (207) is then rewritten in the form:
Y? = Yj -C+ j_1/2 (fj -rU)j^-C]+1/2 (Vj+1 -fj)^-D j-1/2 R j-1/2M-D j+1/2 R j+1/2At = 0 (223)
The first iteration starts with: y/] = Y^ . The predicted variables y/J are accepted and ^"+1 = y/J if:
Yj-Yj
Yj
< L (224)
Otherwise, the solutions are corrected/re-predicted with another iteration, using the relationship y/] = iff**. Matrices C and D are updated for new variables xjf- after each iteration within the time step.
The acceptance criterion i for variables of the Timoshenko beam at the end of each time step checks if the critical variable changes for given value within one time step (currently i = 2). If so, then the time step is repeated. Maximal number of iterations is limited, because it turned out that if the number of iterations exceeds 10, the solution usually diverges. The solution becomes convergent for smaller computational time steps, and for sufficiently small time steps the implicit iterations would not be necessary anymore.
5.5. Numerical errors and difficulties
Numerical dissipation. Tiselj and Petelin [123, 124] stated that the weak side of the first order numerical methods is the numerical dissipation, which tends to smear discontinuities on coarse grids. The numerical dissipation (also numerical diffusion, dispersion) is a consequence of the discretisation of the pure advection equation (which, by definition, is free of dissipation) with spatial and temporal difference schemes that are first order accurate. Numerical dissipation is similar to the physical diffusion (viscosity friction, heat conduction, etc.). It tends to smooth sharp gradients and it stabilizes oscillations of the numerical method. The first error term (also local truncation error) is defined as the difference between the approximate and the true solutions. The leading error term is not included in the upwind difference scheme, therefore, the explicit upwind difference scheme for positive characteristic with the leading error term can be written as:
An+1=An-A^ Ln-Ln -A^(Ax-AAf)^ (225)
j j"f(^-^n)-A2 (AX-AAf)0"
The error term is essentially similar to the dissipation term (second order derivative) with numerical dissipation coefficient e:
A
L = — (Ax - AAf) (226)
- 60 -
A positive dissipation coefficient yields a stable difference scheme. A negative dissipation coefficient increases oscillations and sharpen gradients, thus small disturbances increase. Diffusivity becomes negative if the characteristic Ä is negative (not possible - upwind) or if the time step is too big, hence, the time step must be defined according to the CFL condition in Eq. (191). The dissipation coefficient in Eq. (226) is made up of two terms (Ax and At); the first with stabilizing and the second with destabilizing effect. Both of them summon and minimal numerical dissipation are obtained if the CFL time step At is as close as possible to the ratio Ax /A. For very small time steps the stabilizing term prevails and larger numerical dissipation is introduced into the solution.
Stiff source terms of the Timoshenko beam equations. Simulations show that an increasing number of implicit iterations within one time step introduces additional numerical error, which behaves like numerical dissipation. The number of implicit iterations can be controlled by application of smaller time steps. However, there is an optimal time step for which the numerical dissipation due to small time step in explicit scheme and due to the number of implicit iterations is minimal.
Operator splitting. The operator splitting is a source of a specific non-accuracy that behaves like numerical dissipation, and was discussed by Burman and Sainsaulieu [18], Bereux [11] and Tiselj and Horvat [125]. The problem arises for small relaxation times, where the relaxation time is a time period in which the relaxation quantity approaches to its equilibrium value. With special transformation of the equations and with appropriate complex numerical schemes described by Burman and Sainsaulieu [18], and Bereux [11], the numerical dissipation of the operator splitting method can be avoided. However, these complex procedures are not applied in this thesis.
Numerical evaluation of the eigensystem. A very important step in the characteristic upwind numerical procedure is the diagonalization i.e., the evaluation of the eigenvalues and eigenvectors. The eigensystem is evaluated at each time step for each computational volume. Some simple physical models are analytically solvable, but most of the physical systems discussed in this dissertation are solved numerically. The numerical evaluation of the eigensystem is introduced through the application of the subroutines from EISPACK numerical library. The simulations showed that the results obtained with the same physical model, once solved analytically and then numerically, are identical, which in turn validates the numerical approach to evaluate the eigensystem.
The numerical evaluation of the eigensystem has only one weak side: the eigenvalues and eigenvectors are sorted at output, but for the characteristic upwind numerical method, each characteristics needs to be at the same position in the diagonalized Jacobian matrix C. Eigensystem re-sorting is thus conducted always. Re-sorting is simple when eigenvalues are almost constant with time and position and when the eigenvalues are roughly predictable in advance. However, in two-phase flow, the eigensystem significantly changes during the simulation and the re-sorting becomes very pretentious.
Other sources of numerical errors and difficulties. Various numerical tests performed showed that the influence of other types of numerical errors and difficulties emerging in the code during the simulation is negligibly small:
• The basic variables in fluid equations are written in a non-conservative form. Several tests showed that non-conservation is negligible for simulations of fast transients with short duration.
• The eigensystem is partially evaluated numerically - only slight differences to analytical values are found.
• Artificially introduced relative velocities between phasic velocities and pipe axial velocity in two-phase system enable appropriate sorting of eigenvalues and cause that the system is initially not at rest. The error is small.
• Application of various techniques of extrapolation of the boundary condition influence the solution and stability of the simulation.
- 61 -
5.6. Outlook of the code
The physical model with numerical discretisation was compiled into the computer code using FORTRAN programming language. The FORTRAN files contain over 9000 lines of code. The scheme below describes the basic steps of the algorithm:
Calculate initial values (read input datafile, set initial values (see Section 4.3), set basic parameters)
Time step loop
Calculate boundary conditions (set values in blind volumes, see Section 5.7) Calculate averaged basic variables in accordance with Eq. (209) Calculate characteristic variables and slope limiters (see Section 5.3)
Calculate basic matrices and vectors (A, B, C, S , R, etc.)
Calculate eigensystem (A , L) numerically and/or analytically
Calculate matrices C+, C, D+ and D
Evaluate variables at new time level t with Eq. (223)
Check for stiff source terms of the Timoshenko beam equation (accept solution or iterate Basic loop)
Check for stiff source terms due to the relaxation (two-phase flow only, relax if necessary)
Conclude time step (accept new variables, define new additional variables, outputs, etc.)
End of time step loop
5.7. Boundary conditions
Each simulated piping system or section is divided into N computational volumes. There are two blind volumes at the beginning and at the end of the piping system or section (Figure 5) used for prescription of the boundary conditions. The most important is the first blind volume, while the second blind volume is used only for evaluation of the gradients for calculation of the slope limiters for the second order correction. All boundary conditions fall into one of two fundamental groups of boundary conditions:
• Fixed value; the value at the boundary is fixed at a certain value due to boundary condition (constant tank pressure, closed end) or is prescribed through the junction coupling mechanisms.
• Extrapolation; the value at the boundary is extrapolated from the last few computational volumes near the boundary. Figure 6 shows principles of various techniques.
Transients are generated by time-varying boundary conditions [118]. Table 7 shows boundary conditions for some frequently used types of supports for up to eight equation models.
BLIND VOLUMES
BLIND VOLUMES
PIPE
I i-2 | i-1 l_J_
i i+1 i+2 i+3 ... N-3 N-2 N-1 N
N+1N+2l
._L J
COMPUTATIONAL VOLUMES
Fig. 5: Each pipe is divided into N computational volumes with two blind volumes at both sides of the pipe for definition of the boundary conditions.
jf exponential
linear If straight
N+2
Fig. 6: Various types of extrapolation of the basic variables into blind volumes.
- 62 -
Valve. Boundary conditions in real experimental setups are not abrupt and sharp; the valve closure for example is not instantaneous. Each valve has its own closing function and time. In characteristic upwind numerical method where volume averaged values are used, the flow change function during the valve closure can be described by varying flow, fluid velocity or cross-section. The fluid velocity variation is chosen in the present dissertation and it can change according to the relationship:
u+ v
x0
)T(t)
(227)
where ux is the axial velocity of the pipe, ux0 is the initial velocity of the pipe, vf is the fluid velocity, vf0 is the initial fluid velocity, ± sign depends on the pipe orientation and flow direction, and r(t) is a given empirical closing function of time. Each type of valve has its own closing function. Tijsseling [118] presented the following equation, based on experimental measurements, for non-instantaneous ball valve closing in his experiment:
"T (t)
t
3.53
0.394 1
1.70
for 0\ / y /\\ /a u
\/ A / \\ // \
\l I \ I \\ //' \ / \ /
\KJ V. .v . V x
Centric mass Eccentric mass
W fK\
0.11------ 0.05 ¦ | 01 -0.05 ¦ lil ll r /I fl ft \ i II
T i L
-0.15 ¦ V ni*
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
0.05 0.1 0.15 0.2
Fig. 15: Geometry and lateral Fig. 16: Geometry and lateral
oscillation of the midpoint for the oscillation of the midpoint for the
endpoint oscillation for arbitrarily pipe with centric force and pipe with centric and eccentric
Fig. 17: Geometry and right
equivalent mass.
mass.
supported and loaded pipe.
The additional mass m is included into the simulation as external load through differential term with mT in Eq. (73). If the additional mass m (usually m > mT ) is concentrated in a single volume then the abrupt jump in the mass represents a singularity in the eigensystem that locally increases the numerical dissipation. Table 10 shows that a simulation gives more accurate results for the cases with denser grid and the cases where the mass is distributed over more computational volumes. The analysis was performed with centric mass oscillating pipe problem depicted in Fig. 15.
I
F
^ ^
F
y1
c2
F
y2
K
T
-40 -60
0.05 -0.05
-0.1
-0.15 -0.2
-0.25 -0.3
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Length [m]
0.15 0.1 0.05
-0.05 -0.1 -0.15
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Length [m]
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Length [m] Length [m]
10
0
-10
-20
-30
-40
-50
Fig. 18: Bending, shear force Fig. 19: Bending, shear force Fig. 20: Bending, shear force
and momentum in cantilever pipe and momentum in pipe due to and momentum in cantilever pipe
due to static force at endpoint. static forces and gravity. due to static force at endpoint.
s
L
L
L
0.25
Time [s]
Time [s]
Time [s]
F
q
y
c
c
L
L
L
-1
-2
-3
-4
-5
100
50
50
-50
-100
-50
-100
-150
-200
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Length [m]
- 69 -
Table 10: Accuracy of the simulation depends on the grid and number of computational volumes over which the mass is distributed.
Grid density [comp. volumes] Mass distribution [comp. volumes] Frequency [Hz]
10 1 31.5
30 3 26.3
90 1 36.9
90 3 26.8
90 9 25.1
270 27 24.8
Analytical solution: 24.5
Additional examples are presented in Figures 18, 19, and 20 (diagrams are grouped in columns). The pipes are loaded with forces at various positions (c = c2 = 1.6 m, and c1 = 0.6 m). The external lateral force for the cases of Figs. 18 and 20 is equal to Fy = 100 N and the beam in Fig. 19 is loaded with two eccentric lateral forces, Fy1 = 20 N and Fy2 = 50 N, and distributed mass q = 46.51 N/m. The natural frequency of oscillation of the (unloaded) cantilever pipe in Fig. 18 is f = 16.0 (15.6) Hz, and the natural frequency of the midpoint of the (unloaded) pipe in Fig. 19 oscillates with frequency f = 43.5 (44.1) Hz. Each movement is damped until the equilibrium is reached. Figures 18, 19, and 20 show equilibrium bending, shear force and bending momentum and all simulated variables (blue dots) perfectly match with analytical solution depicted with red continuous line. One should note, that the achievement of the equilibrium is not a consequence of the physical damping, but due to the numerical damping.
Lateral oscillations of a fluid-filled pipe. The planar eight-equation quasi-two-phase flow FSI physical model described in Section 4.2.3 is applied to verify the influence of various fluids on lateral dynamics of the pipe. The geometry and geometric properties are the same as for the structure in Figure 15, the piping system is closed at both ends, pressurized, and filled with four different fluids (properties in Table 8). The system is initially at rest, relative changes in pressure and displacements of the midpoint are traced.
The case where the pipe is not loaded with an external load is considered first. Various fluids have various densities and therefore the total mass of the piping system per unit of length is different for different fluids. At time zero, the gravity starts to act on the piping system and the structure starts to oscillate. Figure 21 shows that for lateral dynamics, the heaviest fluid yields to the largest amplitude and the smallest frequency of oscillation. The structural oscillations through FSI coupling induce oscillations of the pressure in the midpoint however, the pressure changes due to oscillations are actually negligible. It is possible to see that the frequency of the pressure waves is few times larger than the frequency of the lateral movement, which in turn prevents the appearance of significant FSI effects in the lateral direction.
Figure 22 shows a case similar to the previous one with one exception. There is no gravity field, the transient is initiated by the action of the lateral force Fy = 100 N in the midpoint. Then the amplitude is equal for all fluids, because different density of different fluids (mass) affect only the frequency of the oscillations. The heaviest fluid R12 has the largest period of oscillation. The pressure history shows that the frequency of pressure pulsations is much shorter than the frequency of structural dynamics, the heaviest fluid experiences the highest pressure peaks, and the pressure pulsations are actually negligible.
The study of the lateral dynamics in fluid-filled piping system with FSI yields to the following conclusions:
• FSI is very weak in lateral dynamics of structures.
• The most important FSI parameter for oscillation in the lateral direction is the density of the fluid (mass).
• The total pipe volume change is very small for lateral dynamics, therefore, the pressure changes in the fluid are negligible.
• With consideration of the gravity: the denser fluid decreases the natural frequency of the oscillation and increases the amplitude.
- 70 -
Without consideration of the gravity: the denser fluid decreases the natural frequency of the oscillation while the amplitude is not affected.
0 10 20 30 40 50 60
Time [ms]
70
80
90
100
0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6
0
10
20
30
70
80
90
100
40 50 60
Time [ms] Fig. 21: Pressure and oscillation history in the midpoint of the piping system filled with various fluids. The transient is initiated by introduction of the gravity forces.
10
-5 -10
0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6
0 10 20 30 40 50 60
Time [ms]
70
80
90
100
0
10
20
30
70
80
90
100
40 50 60
Time [ms] Fig. 22: Pressure and oscillation history in the midpoint of the piping system filled with various fluids. The structure is loaded with lateral force (no gravity considered).
0
5
0
- 71 -
6.2. Delft hydraulic benchmark problems
All developed FSI computer codes and the underlying mathematical models and numerical methods, have to be extensively tested. For this purpose Tijsseling [113] adopted several benchmark problems, known as Delft hydraulics benchmark (DHB) problems A to F. These benchmark problems are based on an inventory of field cases and were defined in consultation with prof. Wiggert. The DHB problems are:
• DHB problem A: 20 m straight pipe with free massless valve.
• DHB problem B: 20 m straight pipe with fixed valve.
• DHB problem C: 330 m straight pipe with free massless valve.
• DHB problem D: 330 m L-shaped pipeline, free bend 20 m away from fixed valve.
• DHB problem E: 330 m L-shaped pipeline, fixed bend 20 m away from fixed valve.
• DHB problem F: 330 m system with four bends.
6.2.1. Delft hydraulic benchmark problems A and B
TANK
PIPE
VALVE
Initial flow direction
P1
i=1 i+1 i+2 i+3
P2 -p
N-1 NI yc I
777777777?
Table 11: Piping system and water properties for DHB problems A and B.
Fig. 23: Numerical model of the DHB problem A (valve is free) and B (valve is fixed).
Piping system Water
L = 20 m Pi = 1000 kg/m3
R = 0.3985 m K = 2.1 GPa
e = 0.008 m p = ptank = 0 bar
E = 210 GPa v = 1 m/s
pt = 7900 kg/m3 T = 20 °C
v= 0.3
The DHB problems A and B refer to a simple straight tank-pipe-valve system depicted in Fig. 23. The pipe with properties described in Table 11 is filled with water at room temperature. The transient is initiated in the fluid by the instantaneous closure of the valve, which rapidly stops the steady-state water flow. The valve is massless and free for case A and fixed for case B. The results were obtained with application of the linear four-equation physical model for simulations of axial quasi-two-phase flow FSI coupling and are compared to the results obtained with the validated computer code of Tijsseling [113] s1ax01. The s1ax01 is based upon a linear physical model with constant characteristics for axial FSI coupling and is solved with the method of characteristics (for details see Lavooij and Tijsseling [75], Wylie and Streeter [150], or Tijsseling [117, 113]). Recall that the method of characteristics (MOC) actually provides the analytical solutions of FSI problems based on several assumptions that have to be introduced into the physical model: the transient is single-phase, the pipe is uniform, straight, with constant material properties, the state properties of the fluid are constant, all parameters that introduce damping are eliminated, and the source terms are set to zero. These assumptions mean the linearization of the leading partial differential equations. In addition, Wylie and Streeter [150] stated that the numerical error of MOC is partly eliminated by initial tuning of (constant) densities to appropriately adjust wave speeds to conform into equidistantly spaced grid points in distance-time plane (Lavooij and Tijsseling [75, 117]).
To be comparable to the MOC results, all, although non-realistic assumptions were implemented into our physical model and consequently, all advantages of our approach were lost. However, the results enable a direct comparison between the method of characteristics and the characteristic upwind numerical method and thus the validation of the proposed numerical scheme. Figures 24 and 25 show the fluid pressure, the pipe axial force, and the pipe axial velocity histories obtained for the DHB problems A and B at position P2 that is located as close as possible to the valve. The simulation with both approaches yields practically the same results however, detail shows that the results obtained with the method of characteristics are slightly sharper. The differences are attributed to:
• Numerical dissipation. The numerical dissipation that comes with the characteristic upwind method smoothens the discontinuities in the results. The numerical dissipation cannot be
Water
- 72 -
eliminated, but it can be minimized (the applied number of computational volumes N = 400, the optimal time step is reduced by factor 0.2, and the Superbee limiter is applied). The amount of numerical dissipation contributes the largest part into differences between methods. • Position of the measurement. The MOC gives values that are valid for a position exactly at the end of the pipe (point), while the characteristic upwind method, with computational volumes, gives an averaged value over the whole computational volume (volume average). The practical representation of the slightly shifted and averaged measuring position again reduces the sharpness of discontinuities, which is an effect similar to the numerical dissipation.
One can conclude that Figures 24 and 25 validate the characteristic upwind numerical method since the simulation perfectly matches the results obtained with the verified MOC method.
2
1.5 1 0.5
0
-0.5
-1
-1.5
1000
0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
Time [s]
500
-500
0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18
Time [s]
0.2
0.6 0.4
0.2
0
-0.2
-0.4
0
0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
Time [s] Fig. 24: DHB problem A: time history of basic variables at position P2.
Figure 26 shows the pressure history in P2 for the first 4 seconds of the transient without consideration of the damping. The results for the DHB problem B exhibit the secondary frequency of pressure oscillations, which appears due to the Poisson coupling mechanism of the FSI. It is also evident that the maximal pressure peak near the valve is not achieved immediately after the valve closure, but it appears as transient evolves. The maximal pressure and therefore the maximal load of the fluid on the structure are not obtainable without appropriate FSI simulation. It is necessary to stress, that
0
0
- 73 -
simulations were performed without consideration of any of the damping mechanisms. Damping in real experiments is always significant and affects the pressure history.
2 -1.5 -
1 0.5 0 ¦ -0.5 -
-1 "
-1.5 -
-2 -
2
1.5 1
0.5
0
-0.5
-1
-1.5
0
600
400
200
0
-200
0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18
Time [s]
0.2
-400
0
0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18
Time [s] Fig. 25: DHB problem B – time history of basic variables at position P2.
0.2
2 - ll I
lil llllllll A 1.5 - , 1 I I 0.5 I 1 III llllllll
llllllll 'llllllll 0
¦ -0.5 llllllll lllllll III
:ll| 111 v I |||l|l||||l|l||||l|l||||l|l||||l|l -1 ' -1.5 - ' 1 III iflllllll lllllll
-2 - 1 III
1.5 2 2.5
Time [s]
1.5 2 2.5
Time [s]
Fig. 26: DHB problems A and B – Pressure history in P2 for the first 4 seconds of the transient.
Minimization of the numerical dissipation. The numerical dissipation of the applied numerical method is relatively small and manageable i.e. it is possible to control the amount of the numerical dissipation by various parameters of the simulation. One of the possible mechanisms are slope limiters given in Section 5.3. Figures 27 and 28 show that each limiter has a different influence on the sharpness of the result. The sharpest results are obtained with the Superbee limiter, the most diffusive (smearest) solutions, but still 2nd order accurate, are obtained with the Minmod limiter, while the solutions with Van Leer and MC limiters lie between the sharpest and the smearest results. All results obtained with application of limiters are second order accurate on smooth solutions (see first order solution Figs. 27 and 28). The best results were obtained with the Superbee limiter but this limiter, especially in some more complex two-phase simulations, sometimes introduces nonstability into the simulation. The use of the Superbee limiter is recommended in single-phase flow but should be used with caution in two-phase flow.
0
0.5
3
3.5
4
0
0.5
3
3.5
4
- 74 -
2
1.5 1
0.5
0
-0.5
-1
-1.5
0
0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
Time [s] Fig. 27: DHB problem A: pressure time history in P2 with study slope limiters.
2
1.5 1
0.5
0
-0.5
-1
-1.5
0
0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
Time [s] Fig. 28: DHB problem B: pressure time history in P2 with study of slope limiters.
2
1.5 1
0.5 0
-0.5
-1
-1.5
0
0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
Time [s] Fig. 29: DHB problem A: pressure time history in P2 with study of grid refinement.
__ 1.5 ¦
1 f 0.5 0
-0.5
-1
-1.5
0
0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
Time [s] Fig. 30: DHB problem B: pressure time history in P2 with study of grid refinement.
- 75 -
9
The second mechanism that affects numerical dissipation is the nodalization of the computational domain. The importance of this mechanism is shown by a very trivial grid refinement study in Figs. 29 and 30. The sharpest results are obtained with N = 400 computational volumes however, N = 50 computational volumes are sufficient for fast and accurate simulation. The case with 50 computational volumes took one minute of calculation time on Pentium IV 3.0 GHz processor.
FSI coupling mechanisms. There are two major coupling mechanisms for the description of the FSI: the distributed Poisson coupling and the local junction coupling. Figures 31 and 32 show the importance of the coupling mechanisms, because both the Poisson and the junction coupling mechanisms must be accounted for accurate results. The valve in DHB problem B in Fig. 32 is fixed, which prevents junction coupling, thus the case “No Poisson coupling” is equal to the case “No FSI” (valid only for the straight pipe). The DHB problem A pressure history is equal to the DHB problem B pressure history if FSI is not considered. The simulation without the FSI is denominated as standard water hammer simulation, and rather simple box-like pressure histories near the valve are observed due to rapid valve closure (Figures 31, 32, and C-2) in non-cavitating flow. The standard water hammer pressure history at the valve and the first few moments in pressure history at the valve due to the junction coupling mechanism can be reconstructed with Joukowsky’s equation, which defines the pressure rise/drop due to the disturbance in the flow velocity (see also Appendix C):
Ap = PfCfAvf
where
Av
f,before
f ,after
(232)
where cf is the speed of sound, ?f is the density of the fluid, and ?vf is the flow velocity change with vf,before as the fluid velocity before and vf,after as the fluid velocity after the disturbance of the flow. Figures 31 and 32 show that for standard water hammer simulations without consideration of the FSI, the fluid velocity after the valve closure is zero (vf,after = 0) and a typical box-like pressure history is obtained. The continuous (black) line on Figure 31 shows that with consideration of the FSI, just after the valve closure, the fluid pressure rise is smaller than for the no FSI case, while few milliseconds later, the continuous line increases above the standard solution (see the same effects in Figs. 13 and 36). The discrepancies from the standard solution are explained using Joukowsky’s theory: at the beginning, the flow is not entirely stopped although the valve is closed, because the pipe starts to extend. Therefore, vf,after = ux, where ux is the axial velocity of the pipe (valve). Consequently, the pressure near the valve remains lower than that of the fixed valve. The induced axial stress wave travels along the pipe and reflects at the tank as discussed in Section 6.1. When the reflected stress wave reaches the valve again, the valve starts to move in the opposite direction with velocity vf,after = -ux, which corresponds to the case of reversed flow. The fluid is additionally compressed, and the pressure near the valve increases above the no FSI case. The initially simple FSI phenomenon becomes much more complicated as time evolves.
2
1.5
1
0.5
0
-0.5
-1
-1.5
0
0.2
0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18
Time [s] Fig. 31: DHB problem A: pressure time history with study of the FSI coupling mechanisms.
Figure 32 shows DHB case B with a fixed valve where only the Poisson coupling mechanism appears. The pressure gradually increases and later on decreases. This is the effect of the Poisson coupling mechanism. After the valve closure, the pressure rise through the Poisson coupling swells the cross-section of the pipe and generates a new stress wave. The stress wave that travels along the pipe
- 76 -
approximately three times faster than water hammer wave, swells the pipe cross-section and therefore reduces the pressure. This pressure reduction is known as precursor wave. The pressure reduction changes sign at the tank and becomes a pressure increase, and this can be observed on the pressure history near the valve. This effect becomes more intense as time evolves. It is possible to conclude that
simulations of an FSI are not possible without appropriate consideration of the Poisson and junction coupling mechanisms.
1.5 1
0.5
0
-0.5
-1
-1.5
0
0.2
0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18
Time [s] Fig. 32: DHB problem B: pressure time history with study of the FSI coupling mechanisms.
Precursor and successor wave. At the moment of the valve closure in the straight axial system, an axial stress wave and a pressure wave are generated that travel along the pipe at the corresponding characteristic velocity. The axial stress waves travel approximately 3 times faster than pressure waves, and due to the axial extension of the pipe, throughout Poisson effect, cause contraction of the pipe cross-section. Consequently, the pressure is raised. This pressure rise travels along the pipe with the same speed as the axial stress wave in the pipe, so it is faster than other pressure waves in the fluid. This pressure disturbance wave is known as the precursor wave, and can be observed as a small pressure disturbance that travels in front of the pressure wave in the fluid during the transient. It is most evident at the beginning of the FSI transient (Figure 33), while later on, the pressure longitudinal profile becomes rather complicated and the precursor waves are not traceable anymore. The precursor waves are very evident in DHB case B, where only Poisson coupling mechanism is considered. The less famous counterpart of the precursor wave can be observed in the axial force longitudinal cross-section. At the position of the pressure wave, a small disturbance in the axial force is evident. This disturbance appears because the pressure wave swells the pipe and the pipe starts to contract throughout the Poisson coupling mechanism. In analogy to the precursor wave, this wave can be denominated as a successor wave. The successor wave travels in the pipe wall with the characteristic velocity of the pressure waves in the fluid.
Profile at t =0.003 s
Precursor wave
Pressure wave
2
1
Stress wave
/
Successor wave
10 Length [m]
Fig. 33: DHB problem A: Pressure and axial force longitudinal profile at time t = 3 ms with indication of the precursor and successor wave.
0.4
0.2
0
0
2
4
6
8
10
12
14
16
18
20
400
0
0
2
4
6
8
12
14
16
18
20
- 77 -
Common characteristics of the precursor and successor waves are, that they are both observed in one medium, while they are a consequence of the transient in the other medium. The propagation velocity of the transient wave in the first medium is equal to the characteristic velocity of the transient in the second medium. The precursor and successor waves can be described only with appropriate consideration of the Poisson coupling mechanism.
6.2.2. Delft hydraulic benchmark problems C, D and E
The Delft hydraulic benchmark (DHB) problems C, D and E are an extension of the DHB problems A and B. The piping system is longer, and it contains an additional 90° elbow in cases D and E. The material properties are collected in Table 12, the geometrical model for DHB problem C is sketched in Fig. 34 and the model for DHB problems D and E is sketched in Fig. 35. The DHB problems D and E are similar except the elbow that is located 20 meter away from the valve is not supported in case D, while it is fixed in case E. The valve is massless and free for case C and fixed for cases D and E. The piping system is initially filled with water at room temperature. The transient is initiated at time zero when the valve is instantaneously closed and the steady-state water flow is stopped. The flow is single-phase. The fluid thermodynamic state properties are assumed constant. The simulations were performed with the eight-equation linear quasi-two-phase flow model for FSI in arbitrarily shaped planar piping systems.
TANK
PIPE
VALVE
Initial flow direction
P1
i=1 i+1 i+2 i+3
Water
P2 -r
N-1 N I jfč I
Table 12: Piping system and water properties for DHB problems C, D and E.
777777777/
Fig. 34: Numerical model of the DHB problem C.
Piping system Water
L = 330 m Pi = 880 kg/m3
R = 0.1032 m K = 1.55 GPa
e = 0.00635 m p = ptank = 0 bar
E = 210 GPa v = 4 m/s
Pt = 7900 kg/m3 T = 20 °C
v= 0.3
TANK
PIPE
Initial flow direction
P1
i=1 i+1 i+2 i+3
90° ELBOW VALVE P3
Fig. 35: Numerical model of the DHB problem D (free elbow 20 m away from the fixed valve) and E (fixed elbow 20 m away from the fixed valve).
The results for the DHB problem C obtained with the simulations are compared to the results obtained with the verified computer code s1ax01 of Tijsseling [113, 117]. The physical model and assumptions applied in our code have been defined in a way to enable direct comparison of numerical methods (MOC vs. Characteristic upwind) without additional influences. The difference compared to the analysis of Section 6.2.1 is that eight-equation model is applied. Figure 36 shows some most important basic variables and comparisons to the results of the MOC. Numerical dissipation and non-equal measuring position contribute to some small discrepancies. However, the agreement between simulation with our code and verified MOC is perfect, which again, validates the application of the high resolution characteristic upwind numerical method.
Water
- 78 -
2 1.5
1
0.5
0
-0.5
-1
Fig.
1.5 Time [s]
1.5 2
Time [s]
2.5 3
100
0
-100
0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2
Time [s] Time [s]
36: DHB problem C: History of basic variables in position P2 – comparison between
characteristic upwind numerical method.
2.5 3
MOC and
Junction coupling at valve. Figure 37 shows a comparison for three variations of the DHB problem C precisely, the no FSI case is compared to the case with fixed valve (Poisson coupling mechanism only) and the case with free valve (Poisson and junction coupling mechanisms). Figure 37 shows that significantly higher pressures in the considered transient were recorded for the case with free valve, which means that the junction coupling mechanism at the valve plays a very important role. This confirms the generally accepted rule, which states that FSI is more intense in weakly or wrongly supported piping systems. It is possible to conclude that fast operating valves should be adequately supported in order to minimize FSI and transfer of energy from the pipe into the fluid.
10
0
0.5
1
2.5
3
3.5
4
1.5 2
Time [s] Fig. 37: Comparison of pressure histories in straight pipe (DHB C) near the valve for cases without FSI, with FSI with fixed valve and with FSI with free valve.
Junction coupling at an elbow. The simulation results of the DHB problems D (elbow is free) and E (elbow is fixed) are compared with the results of the simulation of the DHB problem C (straight pipe). Figure 38 shows the pressure histories close to the valve. The valve is fixed and closed instantaneously in all considered cases. It is possible to see that the highest pressure peaks (more than 50% higher)
8
0
0.5
2
2.5
3
0
0.5
1.5
300
200
-200
-1.5
0
5
0
- 79 -
were obtained for the case where the elbow was free. This again confirms the general rule, which states that the piping system has to be sufficiently supported to prevent development of the FSI coupling. Every elbow of the piping system has to be adequately supported. Comparison of the case with straight pipe and pipe with fixed elbow gives another interesting conclusion: the pressure history for case with fixed elbow (DHB problem E) is essentially the same as the pressure history for the case with straight pipe (DHB problem C, fixed valve). We can conclude that during the FSI pipe transient the piping system with fixed elbows behaves similarly as the straight piping system.
10 ------Elbow is free Elbow is fixed...........Straight pipe
5 ~ *—\ j-i-Sä. JX-rJÄ \-Aex**.
i 1 v ^Hä i \ ¦ i i \j i \ 7^
0 -5 V i 1 i j ;V i v /J \
Fig
0 0.5 1 1.5 2 2.5 3 3.5 4
Time [s] 38: Comparison of pressure histories in position as close as possible to the valve for cases with free and fixed elbow (DHB D & E) and for case without elbow (DHB C).
Elbows are one of the most common parts of the piping systems, and simulations (and experiences) show that they should be adequately supported. Of course, this is not a trivial task. In the piping systems of a nuclear power plant for example, the thermal extensions of the primary coolant loop piping system lead to up to 20 centimeter extensions/contractions. Fixed supports would impose gigantic loads. The general solution was found in snubbers (mechanical shock arrestors). The snubber allows thermal displacements with velocities lower than threshold velocity (few cm/s), during the water hammer transient occurrences (or plant transients, earthquakes, etc.), velocities of displacements increase over the threshold and the snubber blocks any movement. For example there are approximately 2300 snubbers mounted to support various piping systems in nuclear power plant Krško.
6.3. Valve closure experiment in single elbow pipe
Simpson [103] in 1986 conducted a series of 9 experiments during his Ph.D. research where he experimentally and numerically studied cavitation in single-elbow piping system. The valve closure experiment in the single elbow pipe, also known as Simpson's pipe column separation water hammer experiment [103], became one of the fundamental benchmarks for two-phase flow computer codes because of the simple geometry, initial conditions and clear water hammer initiating mechanism. The valve closure experiment in the single elbow pipe was applied in this dissertation for verification of the quasi-two-phase flow model and to study the influence of the FSI on two-phase flow transients. The simulated results were obtained with the eight-equation quasi-two-phase flow physical model for FSI in planar arbitrarily shaped piping systems, which was solved with the characteristic upwind numerical method. The quasi-two-phase flow model was applied for simulation of the column separation. The results of the simulations were compared to the experimental measurements.
Simpson varied the initial velocity of the steady state fluid to affect the appearance of the two-phase flow and the severity of the water hammer and cavitation. The other initial parameters were held as constant as possible during the execution of the experiments. The matrix with initial conditions for various experiments is presented in Table 13. The transient in initially steady state tank-pipe-valve system (Fig. 40) is induced by instantaneous ball valve closure. The ball valve is closed by hand to avoid the amount of pipe movement and vibration resulting from the fast closure of the valve. The hand operation of the valve leads to longer closing times and smooth and non-standard flow diminishing functions of the ball valve. For the simulation of all experiments we assumed that the fluid temperature is Tf = 297 K, the
- 80 -
valve closing function is defined by Eq. (228), and the applied uniform valve closing time is tc = 10 ms. The piping system in Fig. 40 is 36 meters long with one 90° elbow at 12.5 meter, and the section from 12.5 to 36 meters is slightly inclined (? = 2.44°). The experimental apparatus is fixed with brackets to the wall at approximately 2.5 meters interval and at additional three points: at the tank, at the elbow and at the valve. Simpson [103] claims that two brackets rigidly support the elbow and any movement of the elbow during the transient is expected to be negligible. Because of the fixed supports that are located at the crucial parts of the piping system, it is expected that the junction coupling FSI mechanism will be negligible.
Fig. 39: Plane (left) and side (right) schematic view (not to scale) of experimental apparatus of the valve
closure experiment in the single elbow pipe [103].
VALVE P3
Fig. 40: Numerical model for experimental apparatus of the valve closure experiment in the single elbow pipe.
Table 13: Matrix of initial conditions for experimental runs.
Label Fluid velocity [m/s] Tank and pipe pressure [bar] Valve closing time [ms] Fluid temperature [K] Wall friction coeff. Cavitation
Exp 1 0.239 3.369 29 297.1 0.0325 No
Exp 2 0.332 3.281 35 297.6 0.0315 Yes
Exp 3 0.401 3.281 17 296.5 0.0290 Yes
Exp 4 0.466 3.259 34 296.5 0.0280 Yes
Exp 5 0.507 3.244 25 296.5 0.0280 Yes
Exp 6 0.596 3.265 29 297.1 0.0270 Yes
Exp 7 0.696 3.233 25 297.1 0.0260 Yes
Exp 8 0.938 3.196 35 297.1 0.0240 Yes
Exp 9 1.125 3.118 43 297.1 0.0230 Yes
Table 14: Material properties of experimental apparatus (Tf = 297 °K).
Material Elasticity modulus [Pa] Poisson ratio Density [kg/m3] Inner pipe radius [m] Pipe wall thickness [m] Effective speed of sound in fluid [m/s]
Copper 1.19E11 0.34 8920 9.525E-3 1.588E-3 1355
Alloy 0.75E11 0.25 1280
- 81
18
The material properties of the piping system are described in Table 14. Simpson [103, p.97] stated that the applied pipe was made of copper and consideration of the copper material properties and pipe geometry properties results in the following effective speed of sound in fluid: cf = 1355 m/s. Simpson performed a spectral analysis of all experimental data and as a result he obtained an effective speed of sound equal to cf,real = 1280 m/s, which was actually applied in his numerical studies. Simpson didn’t give any explanation for the discrepancy between the real and theoretical values for effective speed of sound. Kovač [71] showed that the Young modulus and other material properties depend on orientation of the crystals and therefore properties of a monocrystal material significantly changes compared to a policrystal material. It is possible that orientation of the copper crystals was affected during the manufacturing of the pipe. Another possible reason can be, if instead of the pure copper, a copper alloy was used as pipe material. Simulations showed that the properties for the undefined alloy give very good agreement with the experimental results. The alloy elasticity modulus was calculated backwards from the Korteweg’s equation (D-9), and the Poisson ratio was estimated based on comparisons between the experiment and simulation. A material with such properties was not found in the literature; we will denominate it in this dissertation as Alloy.
The fluid pressure was measured along the pipe at three positions that are described in Table 15. In addition, the pressure in the tank was measured (manually) with a Mercury manometer. Other variables were not measured.
Table 15: Position of the measuring equipment and elbow.
Label Position [m] Volume No. Variable Equipment
Tank 0.00 0 Pressure Mercury manometer
P1 9.00 25 Pressure Pressure transducer
Elbow 12.50 35
P2 27.00 75 Pressure Pressure transducer
P3 36.00 100 Pressure Pressure transducer
0.6 0.5
0.4 0.3
0
0.7 Experiment Alloy Alloy, no FSI Copper
0.6 /A If | If * A A &
0.5 i/ j J i / 1 l \ '" *»/ 1 1 I \
0.4 0.3 a \\ i 1 i i / A Vi / il a 1 I \ s\ f' i A \ i J \ \ *-•* /1 \ V ^
0.2 \. Vi \
0.1 V/I W \; 1
0.05
0.1
0.15 0.2
Time [s]
0.25
0.3
0.35
0.4
0.6 0.5
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
Time [s] Time [s]
Fig. 41: Pressure history at positions P1, P2 and P3 for single-phase Exp 1. Comparison between experiment and simulation with consideration of the Copper or Alloy material properties.
0
0.1
0.1
0.4
- 82 -
0
0
1.2 1 0.8 0.6 0.4
0.2
\
Experiment Alloy Alloy, no FSI
Copper
; *
0.05
0.1
0.15
0.2 Time [s]
0.25
0.3
0.35
0.4
0
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
Time [s] Time [s]
Fig. 42: Pressure history at positions P1, P2 and P3 for moderate cavitation Exp 3. Comparison between experiment and simulation with consideration of the Copper or Alloy material properties.
Experiment------Alloy..........Alloy, no FSI
Copper
0.05 0.1 0.15 0.2 0.25
Time [s]
0.3
0.35
0.4
0.15 0.2 0.25
Time [s]
2 1.8 1.6 1.4
1.2
1
0.8 0.6 0.4 0.2 0
0.2 Time [s]
Fig. 43: Pressure history at positions P1, P2 and P3 for severe cavitation Exp 9. Comparison between experiment and simulation with consideration of the Copper or Alloy material properties.
0
0.6
0.6
0.4
0.4
0.4
2
1
0
0
1.2
0.6
0
0.05
0.1
0.3
0.35
0.4
0.05
0.1
0.15
0.25 0.3
0.35
0.4
- 83 -
Figures 41, 42, and 43 show a pressure history comparison between the simulation and the measurement at various locations along the pipe, for three representative experiments (Single phase experiment Exp 1, experiment Exp 3 with moderate cavitation and experiment Exp 9 with severe cavitation). The continuous (red) line was obtained with the assumption of the alloy material properties (Table 14), the dotted (blue) line stands for default copper material properties, and the dash-dotted (green) line stands for Alloy but without consideration of the FSI coupling mechanisms. The simulation with imaginary Alloy material properties perfectly match with the experiment in all cases. It is also possible to see, that the results for Alloy with and without consideration of the FSI are almost identical. This points to the fact, that sufficiently supported piping system, like Simpson pipe is, prevents junction coupling FSI during the column separation water hammer transient. The pipe has a small cross-section and a relatively small radius over thickness ratio (R/e = 5.9) and consequently the Poisson coupling is present but it is so weak, that it can be easily neglected. The simulations show that the influence of the FSI coupling is negligible for the Simpson pipe experiment. These findings are further supported by Figures from 44 to 47.
0.6
0.5 0.4
0.3
0.2
0.1
0
0
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
Time [s]
0.7 0.6
0.5 0.4
0.3
0.2
0.1
0
0
0.05 0.1 0.15 0.2
Time [s]
0.35 0.4
Fig. 44: Pressure history at positions P1 and P3 for Exp 2.
0
0.2
0
0.4
0.2
0 0.15 0.2 0.25 0.3 0.35 0.4 0 0.05 0.1 0.15 0.2
Time [s] Time [s]
Fig. 45: Pressure history at positions P1 and P3 for Exp 4.
0.35 0.4
0.15 0.2 0.25
Time [s]
Experiment Alloy
0.15 0.2
Time [s]
0.35 0.4
Fig. 46: Pressure history at positions P2 and P3 for Exp 6.
0.25
0.3
0
0
0.05
0.1
0.25
0.3
0.8
0.8
0.6
0.6
0
0
0.05
0.1
0.3
0.35
0.4
0
0.05
0.1
0.25
0.3
- 84 -
0.2
1.8 1.6 1.4 1.2 1
0.6 0.4 0.2
0.15 0.2
Time [s]
0.15 0.2
Time [s]
Fig. 47: Pressure history at positions P2 and P3 for Exp 8.
Quasi-two-phase flow model. Gale and Tiselj [43] described their simulations of the Simpson’s pipe experiment with the six-equation two-fluid WAHA code (Tiselj et al. [126]) and showed that one does not need a complete two-fluid model for successful simulation of the column separation water hammer. Especially at low temperatures, flashing and condensation of the steam are not governed by the heat and mass transfer between the phases, but by dynamics of the liquid column. Thus, modeling of the transient in water at room temperature usually does not require energy equations (Bergant et al. [12]) and simplified cavitation models can be applied.
0.3
0.25
0.2
0.15
0.1
Time [s]
Pipe length [m]
Time [s]
Pipe length [m]
Fig. 48: Vapor volume fraction in time-space plane for the case with fixed valve (Exp 3 and 9).
_. r" " " ¦-,
x 10 - r- X X
4^ , - ~~ k - k " I v x v
¦
3- ,,- k " " k ' k - " V- " " k - t -v v-
2" ,- k " - k " , I- - " L- - " L- - " V- " " "j 1 ^ t ^ v v
0.4
•«^fc
0.35 \>
0.3
0.25
0.2
0.15
Time [s]
|8j§g^^r 25
15
Pipe length [m]
j^
0.03.
0.02,
0.01 ,
0.
0.3
0.1
Time [s]
^
Pipe length [m]
Fig. 49: Vapor volume fraction in time-space plane for the case with free valve – with FSI (Exp 3 and 9).
0
0
0
0.05
0.1
0.25
0.3
0.35
0.4
0
0.05
0.1
0.25
0.3
0.35
0.4
0.04
0.03
0.02
0.01
0.3
35
35
30
25
20
15
10
0 0
0 0
0.04
0.4
0.:
35
35
0.2
30
30
25
20
0.1
15
10
0.05
- 85 -
Figures from 41 to 47 show and validate at the same time, that a relatively simple quasi-two-phase flow model describes the column separation type of the two-phase flow in water at room temperature with great accuracy. The duration of the cavitating flow in pressure history is exactly simulated, while it cannot be anticipated, that simulation will catch all small disturbances in a measured pressure history. Figure 48 shows the distribution of the vapor volume fraction in the time-space plane for two cases: moderate cavitation (Exp 3) and severe cavitation (Exp 9). In both cases, the initial cavitating area extends over larger sections of the pipe and later on appears also in the middle of the pipe. Figure 50 shows the vapor volume fraction (VVF) history in the closest volume to the valve. The growth of the vapor bubble is made up of two phases: the first phase where the vapor bubble grows and the second phase, shorter, where the vapor bubble collapses. The simulation shows that the length of the cavitation phase, vapor generation gradient, and amount of the generated vapor linearly depend on the initial fluid velocity. This is actually a property of the chosen quasi-two-phase flow model. The variable initial fluid velocity in lesser extent influences the length of the cavitating area and to a larger extent influences the amount of the generated vapor. The total amount of the vapor is exactly predicted with our code, while it is difficult to discuss or validate the distribution of the vapor without comparison to the experiment.
0.025 0.02
0.015
0.01
0.005
0
0
0.015
0.01
0.005
0
0
2 Exp 3 Exp 7 Exp 9
1.8
1.6 -/ "M /';
/\ ¦
1.4 J
1.2 I—v^ /1
1 0.8 ! "~~1 J ! ,> /i rJ-
0.6 / / v/ i t
0.4 ,
0.2
.1 , \ l f\
0.15 0.2 0.25
Time [ms]
0.15 0.2 0.25
Time [ms]
Fig. 50: VVF and pressure history at position P3 for the case with fixed valve.
2 1.8 1.6 1.4 1.2
1 0.8 0.6 0.4 0.2
0.15 0.2 0.25
Time [ms]
0.15 0.2 0.25
Time [ms]
Fig. 51: VVF and pressure history at position P3 for the case with free valve – with FSI.
FSI coupling mechanisms and cavitation. The Simpson pipe experimental apparatus is relatively strongly supported and thus enables only development of the Poisson coupling FSI mechanism. To include the junction coupling mechanism into the transient, the simulations were repeated with a free valve and released supports in the axial direction. The results with the free valve and the junction coupling are presented in Figures 49 and 51 and some interesting differences in histories of basic variables can be found. Figure 51 shows that junction coupling affects the vapor generation rate. The vapor is generated faster or slower depending on the relative velocity between the pipe (valve) and the fluid. Figures 48 and 49 show the vapor volume fraction in the time-space plane. The vapor is distributed along larger sections of the pipe for the cases with the free valve than for the cases with a fixed valve, especially at the beginning of the transient pipe flow. The reason for the appearance of the additional distributed cavitation is FSI coupling, namely, significant axial stresses are generated at the valve through junction coupling, and as stress waves travel along the pipe, through Poisson coupling
0.035
0.03
0.05
0.1
0.3
0.35
0.4
0
0.05
0.1
0.3
0.35
0.4
0.035
0.03
0.025
0.02
0
0.05
0.1
0.3
0.35
0.4
0
0.05
0.1
0.3
0.35
0.4
- 86 -
generates additional pressure waves in the fluid denoted precursor waves. The precursor waves are added to the basic pressure waves and in case that two ‘negative’ waves are added, the pressure may occasionally drop below saturation and the cavitation starts. Tijsseling [113] denoted this type of cavitation as the Poisson coupling induced cavitation. As distributed cavitation may appear anywhere along the pipe according to the pressure profile, the Poisson coupling induced cavitation appears only as a consequence of the (negative) precursor waves. As time evolves, cavitation generates more and more pressure waves that travel along the piping system and add to the complexity of the pressure profile. In such cases, the precursor waves become more and more hidden and less important. Therefore, the Poisson coupling induced cavitation can be detected probably only at the beginning of the transient.
6.4. Rod impact experiments in hanging piping systems
Vardy and Fan [132] studied Fluid-Structure Interaction without cavitation in liquid-filled piping systems. An extensive series of high-quality measurements was carried out on four different piping systems in a test rig built in the Hydraulics Laboratory at the University of Dundee (UK). The interaction between the fluid and the structure in these piping systems is very strong. The piping systems are suspended by long, thin, vertical steel wires from the ceiling, and can move freely in a nearly horizontal plane. All piping systems are filled with pressurized tap water. The impact rod is used to simultaneously initiate stress waves in the pipe wall and pressure waves in the water by the axial impact of the rod on the left free end of the piping system. The experiments on the wired piping systems were carried out with intention to reduce all influences of the complex initial and boundary conditions, which usually affect the transient. The hanging experimental apparatus actually isolates the effects of the fluid-structure interaction in the case of axial wave propagation. Compared to the standard tank-pipe-valve system, the following unwanted influences were eliminated:
• Influence of the initial steady-state distribution of the fluid properties (pressure, temperature, velocity, density, etc.).
• Influence of the closing valve (closing time, closing function, induced vibrations).
• Influence of the boundary conditions and supports of the piping system.
Therefore, the only measured input in the simulation is the impact rod velocity that is theoretically equal to (with very small deviation in practice):
V°d = 29hr0d (233)
where hrod is the elevation of the rod in raised position.
In his dissertation, Tijsseling [113] studied the cavitation and developed a computer code based on the method of characteristics for simulations of cavitation during FSI occurrence. At the test rig of Vardy and Fan [132] in the Hydraulics Laboratory at the University of Dundee, Tijsseling performed an additional set of experiments with reduced initial pressure in the piping system and thus included cavitation. Tijsseling [39, 113] stated that, at the time of his investigation, the FSI experiments with cavitation were the only well-documented experiments in which FSI and cavitation occur simultaneously and are both significant. Tijsseling provided several experimental data measured during a two-phase rod impact experiment and published them at his web site dedicated to the FSI phenomenon (www.win.tue.nl/fsi/). These experiments quickly became fundamental benchmarks for the FSI two-phase flow codes. The single-phase rod impact experiment in a straight pipe performed by Vardy and Fan and the two-phase flow rod impact experiment in a single elbow pipe performed by Tijsseling are simulated and discussed in the following subsections of this Section.
6.4.1. Rod impact experiment in a straight pipe
The rod impact experiment in the straight piping system consists of a stainless steel pipe closed at both ends with an end plug and an end cap (mass m1 and m2). The piping system presented in Figure 52 is filled with pressurized tap water and suspended by two long (about 3.3 m) thin, vertical steel wires from the ceiling. The impact rod is used to simultaneously initiate a transient by the axial impact of the rod on
- 87 -
the left free end of the piping system. The impact rod and the piping system are separated when the contact force becomes tensile, which happens 2 milliseconds after the strike. The material properties and the initial conditions of the considered piping system, water and impact rod are collected in Table 16. Friction and gravity effects are unimportant and were not considered. The measured initial flexure of the pipe caused by its weight is less than 0.2 mm relative to the suspension points, which are situated about 0.95 m from the pipe ends. The experiment was simulated with the four-equation axial quasi-two-phase flow FSI coupling model.
D2
-ED = pressure transducer ™ = strain ganges
Impact
remote end
P1P2 P3
=*i
P4
P5 P6
i=1 i+ 1 i+2 i+3
Fig. 52: Sketch of experimental apparatus of the rod impact experiment from Tijsseling [113] (above), and nodalization of the piping system for simulation (below).
Vardy and Fan [132] actually performed several single-phase experiments with variable initial pressure and velocity of the impact rod. The pipe in the considered experiment is pressurized with an initial pressure p = 20 bar, which is sufficiently high to prevent cavitation during the transient. Wiggert and Tijsseling [144] stressed that the underlying theory is almost linear for single-phase transients, thus, pressure, velocities and strains are linearly proportional to the impact velocity of the rod. The comparison of the measurement and simulation for one particular impact velocity are considered to be sufficient due to linearity.
Table 16: Material and state properties of the piping system, water and impact rod.
Piping system Water Impact rod
L = 4.502 m Pi = 999 kg/m3 Lrod = 5.006 m
R = 0.02601 m K = 2.14 GPa Rrod = 0.02537 m
e = 0.003945 m p = 20 bar Erod = 200 GPa
E = 168 GPa v = 0 m/s Prod = 7848 kg/m3
Pt = 7985 kg/m3 T = 20 °C v0,rod = 0.739 m/s
v= 0.29 Yrod = 80109.7 kg/s
m-12- 1.2866 / 0.2925 kg
The piping system was extensively instrumented [39, 132, 113]. Table 17 shows a summary of the instrumentation used to obtain data used for comparisons with results of simulation. The axial force was not directly measured but it is related to the averaged axial strain (measured at 4 positions around the circumference) and pressure throughout relationship:
Nx = | EStex + vSt — p
R e
where ?x = (?x,1 + ?x,2 + ?x,3 + ?x,4) / 4
(234)
The experimental data obtained during the rod impact experiments are very valuable, because of the simple and clear geometry and transient initiating mechanism. In comparison to the Simpson pipe experiment, the rod impact experimental data exhibit a higher level of interaction between the fluid and the structure. The numerical reproduction of the results without consideration of the Poisson and junction coupling mechanisms is not possible. These experiments therefore serve as an excellent
Water
N-2N-1 N
- 88 -
benchmark for any numerical simulation that considers the FSI phenomenon. Figure 53 shows a comparison between the experimental measurement and simulation with our code for all available measured variables (pressure in P1, P4, P6, axial force in P5 and axial pipe velocity in P2). The agreement is excellent.
Table 17: Position of the measuring equipment from the impact (left) end.
2.5 <¦ 2 -
1.5 -1 -
0.5 -0 -
2.5 2
1.5 1
0.5 0
Label Position [m] Variable Equipment
P1 0.0195 Pressure Piezoelectric pressure transducer
P2 0.0465 Axial velocity Laser-Doppler vibrometer
P3 0.5740 Axial strain Strain gauges (4 records)
P4 2.2510 Pressure Piezoelectric pressure transducer
P5 3.9440 Axial strain Strain gauges (4 records)
P6 4.5020 Pressure Piezoelectric pressure transducer
Experiment Simulation
1 nA
1 "A-o^^J , r4
i fj m -
Hi i ,j
2.5
2 1.5
1
0.5
0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20
Time [ms] Time [ms]
10
0
-15
2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20
Time [ms] Time [ms]
1.5
0
Time [ms]
Fig. 53: History of measured basic variables (pressure in P1, P4, and P6, axial force in P5, and axial pipe velocity in P2). Comparison between the experiment and simulation.
The good agreement between the numerical simulations and the measurements allows us to display and explain the physical phenomena in the relatively simple experiment with the help of a comprehensive computer output. Figure 54 very illustratively shows the basic variables in the timespace plane. From the aspect of the integrity of the structure, the most important variables are pressure
5
-5
-10
0
0
2
4
6
8
10
12
14
16
18
20
- 89 -
and axial force. Figure 54 shows that the pressure is strongly correlated with the axial pipe velocity. The direction of movement of the pipe actually determine the pressure profile. If the pipe moves to the right, the pressure at the right end is low and the pressure at the left end is high. Therefore, the high pressure always appears in a pair with low pressure (like the children’s swing). The pipe is not loaded with the pressure in the middle of the pipe (axis of the swing). Timing and position of the maximal pressure is not known in advance, all what is certain is, that the maximal pressure will appear in the vicinity of the left or right end. The largest compressive axial force appears just after the strike of the impact rod. When the stress wave travels along the pipe for the first time and reflects from the remote end, it shares energy with the fluid through the junction coupling and the axial force becomes much more moderate. In addition, the energy from the impact partly transforms into kinematic energy (movement) of the piping system. The largest compressive axial force appears at the beginning of the transient along the whole pipe, while the largest tensile axial force appears as transient evolves at arbitrary time and position. Similar analyses can be performed for other parameters of the simulation.
Fig. 54: Basic variables in time – space plane.
Generation of new waves in the piping system. The absolute value of the pressure or axial force at any time and any position is a superposition of all stress and pressure waves that travel along the pipe at current instance. The superposition is initially simple, made up of only few waves (example: after the rod impact, the precursor wave is added to the pressure wave), but it becomes quite complicated as time evolves. Each pressure (stress) wave through the junction coupling at geometric singularities generates new stress (pressure) wave, therefore, the number of stress and pressure waves in the piping systems is rapidly increasing with time. More precisely, the pressure and the stress wave are generated simultaneously at the impact end (left), therefore the faster stress wave is the first that reaches the remote end (right) of the pipe (Figure 54). The stress wave pushes the remote end away from the liquid so that the pressure drops and a new negative pressure wave is generated. The new pressure wave travels backwards and when it meets the forthcoming pressure wave, they are added together (the
- 90 -
precursor wave is also added). Afterwards, the first pressure wave reaches the remote end and reflects from the end cap. The end cap is pushed further away, the pressure is raised and a new tensile stress wave is generated. Only the first reflection of the pressure and the stress wave at remote end was described above but in general, any reflection of any wave through the junction coupling generates new waves.
Linear and nonlinear physical models. The general form of a one dimensional nonlinear system of conservation laws is given as:
M+*(**,') = 0 (235)
dt ds
Linearization of the nonlinear system with application of the Jacobian matrix C yields equation in the following vectorial form:
M + C(^x,0^ = 0 (236)
dt v ds
The system of equations written in form of Eq. (236) is still nonlinear however, LeVeque [77] called it quasilinear because it resembles a linear system. The balance equations for the fluid developed in Chapter 2 are nonlinear. Moreover, every set of partial differential equations considered in this dissertation can be written in vectorial form by equation (236) with corresponding Jacobian matrix. The character (dependency) of the Jacobian matrix is very important for the selection of the numerical method. The following systems are possible:
• Linear system with constant coefficients:
• Linear system with variable coefficients:
• Nonlinear system:
The linear model with constant coefficients is the only acceptable form that can be solved with the method of characteristics, because these systems have constant eigenvalues with time and position, and such systems are mostly used by the FSI community. The model is also known as a standard physical model, and the solution with the method of characteristics is known as a standard FSI procedure. The novelty of our approach introduced in this dissertation is application of the characteristic upwind numerical scheme, which provides numerical solutions of nonlinear systems and linear systems with variable coefficients in space and time. The ability of the characteristic upwind numerical method to solve nonlinear systems is an advantage of great practical importance because it enables application of advanced physical models. Figure 55 shows eigenvalues history in P6 for case with linear constant coefficient system (labeled Linear-const.), which is compared to the system with real water properties (labeled Linear-var.), and to the system with included convective terms (labeled Nonlinear). Figure 56 left shows eigenvalue profiles at time t = 10 ms, where constant coefficient system is compared to the case where masses of the end cap and end plug are taken into account. Figure 56 right shows a similar comparison, where real water properties are used. Figures 55 and 56 show, that the characteristic velocities (eigenvalues) that appear in simulations of experimental (real) cases are not constant. They change with time and position however, Figure 57 shows that for the single-phase rod impact experiment in the straight piping system, consideration of advanced physical models yields relatively small differences to the results obtained with standard linear constant coefficient system. The solution with nonlinear matrix C (labeled Nonlinear) is actually identical to the linear constant matrix C (labeled Linear-const.) due to the acoustic approximation (vf « cf). The advanced physical models become much more important for accurate simulations of the two-phase flow FSI transients, for simulations of FSI transients in piping systems of more complicated geometry or for simulations of transients where the fluid velocity is not small (convection).
dyr dt ds = 0
dp + dt C(x,0 d\j/ ds~
dp + dt C(y/,x ' dS
- 91 -
Improvements of the physical model that affect character of the Jacobian matrix and that were applied in dissertation consist in taking into account the convective terms, exact fluid thermodynamic properties (density, compressibility, speed of sound, etc.), various geometric changes like cross-section changes, variable pipe wall thickness, elbows, flexibility due ovalization, masses and loads, elastic supports, etc. Several examples of use of these improvements with discussion on advantages are given in the following sections.
1363-
---- inear-co nst. ......... Linear-v ar. — Nonlinear fr
I ^1-^-.-
1 1
1
0 2 4 6
-------Linear-const.......... ¦' j Linear-var. - ------Nonlinear ¦
n ! !«(
l!
8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20
Time [ms] Time [ms]
Fig. 55: History of eigenvalues in P6 obtained with various systems of equations.
Profile at t = 10 ms
Profile at t = 10 ms
4.5 4
3.5 3
2.5 2
15 1
0.5 0
Linear-const. Linear-var.
i:
Linear-const. Linear-var.
1.5 2 2.5 3
Length [m]
1.5 2 2.5 3
Linear-const. Linear-var.
2 2.5
Length [m]
Fig. 56: Eigenvalues profiles obtained with various systems of equation.
0
0 2 4 6
8 10 12 14 16 18 20
Time [ms]
0 2 4 6
8 10 12 14 16 18 20
Time [ms]
Fig. 57: Pressure history in P6 and axial force history in P5 obtained with various systems of equations.
Mass of the end cap and end plug. An external load that is mounted on the pipe (end cap, end plug, and other objects with mass) affects inertia of the piping system, and therefore also the eigensystem itself. The eigenvalues of such system are not constant in space (Figure 56 left) but they are still constant in time. Figure 56 shows that the eigenvalues for the rod impact experiment are reduced at the
1365
4618.1
1361
4617.9
1359
1362
1366
1364
1361.5
1362
1361
1360
1360.5
1358
0
0.5
1.5
2.5
3
3.5
4
4.5
0.5
3.5
4
4.5
5000
4618.2
4500
4618.1
4000
4618
3500
4617.9
3000
4617.8
0
0.5
3.5
4
4.5
0.5
1.5
3.5
4
4.5
10
5
-5
-10
- 92 -
beginning and at the end of the experimental apparatus due to the end plug/cap mass, which is numerically distributed over the first/last few volumes. The effect of the mass of the end plug can be physically compared to the mass hanging on the spring. A heavier mass gives a longer period of natural oscillation and smaller natural frequency and the characteristic velocity is therefore reduced to account for this effects. Consideration of the end mass noticeable improves accuracy of results, therefore all results in Figure 53 were obtained with consideration of the additional mass of the piping system (differences are indicated in Figure 57).
Friction and thick-walled model. Figure 58 shows the wall friction forces, obtained with a constant Darcy-Weisbach factor f = 0.01. The influence of the wall friction force is negligible for the considered experiment (small relative velocity between fluid and pipe wall) although it turned out that the friction forces are important for accurate simulation of the Simpson’s pipe experiment (see Table 13). The thick-walled model described by Eq. (159) affects the characteristic velocities of the pressure/stress waves. The characteristic velocity of the pressure wave in the fluid is reduced to cf = 1350.6 (1361.7) m/s and the characteristic velocity of the stress wave in the pipe is reduced to ct,x = 4615.2 (4617.9) m/s. Values in brackets were obtained with the default (thin-walled) model. Figure 59 shows that the reduction of the characteristic velocities of the pressure wave in thick-walled model affects the course of the transient.
2.5
2.5
2 1.5
1 0.5
0 2 4 6 8 10 12 14 16 18
Time [ms]
Fig. 58: Pressure history in P6 – wall friction forces are negligible.
3
2 4 6 8 10 12 14 16 18
Time [ms]
Fig. 59: Pressure history in P6 – influence of the thick-walled model.
Processor time consumption. The number of waves in the piping system is rapidly increasing with time and this may become quite processor demanding for certain numerical methods like method of characteristics where each wave (characteristics) needs to be detected and traced with time. The advantage of the characteristic upwind numerical method is that the processor time consumption during the simulation is constant regardless the number of pressure, stress and other waves in the system. Actually, only the relaxation in the two-phase flow (see WAHA manual [126]) or eventually the implicit iterations for the solution of the stiff source terms can slow down the simulation.
3.5
2.5 2
8 10 12 14 16 18 20
Time [ms]
Fig. 60: Pressure history in P6 – grid refinement study.
20
20
4
1.5
0.5
0
0
2
4
6
- 93 -
Figure 54 shows the results obtained with the numerical model that contains 45 computational volumes. The simulation itself on a Pentium IV 3.0 GHz processor lasts 50 seconds (extensive output to files and screen). The grid refinement study in Figure 60 shows that results are almost independent of the grid. The rough grid with N = 45 volumes is sufficiently dense and accurate since the results obtained on a denser grid with N = 90 and N = 450 computational volumes are not essentially better. The total computational time spent was 30 seconds (reduced output to file) on a P IV 3.0 GHz computer for the rough case and 2 minutes and 50 minutes for the cases with 90 and 450 volumes, respectively.
6.4.2. Rod impact experiment in a single elbow pipe
The rod impact experiment in a single elbow piping system is the second piping system of the experimental work of Vardy and Fan [132] who studied Fluid-Structure Interaction in wire hanging liquid-filled piping systems. Tijsseling [39, 113] improved the set of initially single-phase experiments with two-phase flow experiments. A detailed description of the experiments is given by Tijsseling, Vardy and Fan [114]. Figure 61 shows the geometry of the single elbow piping system that is hanging on three steel wires. The system consists of a stainless steel pipe closed at both ends with an end plug and an end cap (mass m1 and m2) and is filled with pressurized tap water. All geometry and state properties for the piping system, water and impact rod are collected in Table 18. The appearance of the two-phase flow depends on the initial pressure in the system, while the other parameters are held constant as much as possible. The transient starts when the impact rod axially strikes the left end of the piping system.
1210
125
^C^ I
Fig. 61: Sketch of the experimental apparatus from Vardy and Fan [114] (above), and numerical model of the rod impact experiment in single-elbow piping system (below).
Table 18: Material and state properties of the piping system, water and impact rod.
Piping system Water Impact rod
L = 4.51 + 1.34 m v = 0 m/s Lrod = 5.006 m
R = 0.02601 m K = 2.14 GPa Rrod = 0.02537 m
e = 0.003945 m pif = 2.0 MPa Erod = 200 GPa
E = 168 GPa p2F = 0.30, 0.67, 0.87, 1.08, 1.24 MPa Yrod = 80109.7 kg/s
?t = 7985 kg/m3 T = 20 °C v0,rod = 0.809 m/s
v= 0.29 ?f = 999 kg/m3 ?rod = 7848 kg/m3
m12 = 1.312 / 0.3258 kg
Tijsseling [113] and Tijsseling, Vardy and Fan [114] described the extensive instrumentation that was mounted on the piping system to collect data from the experiment. Table 19 shows only a summary of the instrumentation applied in the present dissertation. The axial force is related to the averaged axial strain and pressure, and the bending momentum is related to the top and bottom axial stresses:
- 94 -
Nx = j EStex + vSt
R
and
M
Elt
x,1
+ L,
x,3
R +e
(237)
Table 19: Applied measuring equipment and position measured from the impact end.
0
-1 -1.5
1 0.8 0.6
0
0.25 0.2
0.15 0.1
0
0.05
0
-0.05
-0.1 -0.15
-0.2 -0.25
Label Position [m] Variable Equipment
P1 0.0195 Pressure Piezoelectric pressure transducer
P2 0.0465 Axial velocity Laser-Doppler vibrometer
P3 0.5740 Axial strain Strain gauges (4 records)
P4 4.6400 Pressure Piezoelectric pressure transducer
P5 5.2500 Axial strain Strain gauges (4 records)
P6 5.8500 Pressure Piezoelectric pressure transducer
0.5
0
-0.5
-1
-1.5
-
2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20
Time [ms] Time [ms]
0.5 0.4 0.3 0.2 0.1 0 -0.1
-0.2
-0.3 -0.4
0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20
Time [ms] Time [ms]
Fig. 62: Histories of measured basic variables (pressure in P1 and P6, axial force in P5, and axial pipe velocity in P2 and momentum in P3 and P5). Comparison between the experiment and simulation.
Single-phase transient. A single phase transient with an initial pressure p = 20 bar was considered first. The end masses and the flexibility factor of the elbow improve the accuracy of the simulation and were taken into account. Figure 62 shows a comparison between the experiment and the simulation for all available measured variables. The general agreement is very good for all variables; small
2
10
1.2
5
0
-5
0.4
0.2
-15
- 95 -
discrepancies are attributed to the complex 3D phenomena at the elbow which are linearized and solved with the 1D physical model. Figure 63 gives an additional illustrating overview of timing and position of extreme values of basic variables in the time-space plane.
Fig. 63: Basic variables in time-space plane.
Flexibility factor – ovalization of the elbow. In 1911, Von Karman [133] explained the complex behavior of an elbow under in-plane bending. He identified the ovalization of the circular cross-section of the pipe under the bending. Von Karman showed that using simplifying assumptions, the elbow is much more flexible than an equivalent straight pipe and that more complex stress distribution is induced. He introduced the concept of a flexibility factor and stress intensification factor that compares the flexibility under bending and the maximum stress to those of an equivalent straight pipe. Dodge and Moore [33] reported that ovalization gives the elbows an elastic flexibility that is 5 to 20 times more than the flexibility of the same size straight pipe. The flexibility of the elbow is accompanied by stresses and strains that are typically 3 to 12 times those of a straight pipe under the same load. The Von Karman’s concept essentially remains unchanged in today’s state of the art piping design and analysis software as specified by various design codes and standards [84, 9]. The effect of ovalization is taken into account as a reduction of the moment of inertia of the pipe's cross-section for flexibility factor k. In accordance to the ASME code (ASME, B&PVC, Class 1 components, NB-3686) the flexibility factor yields:
- 96 -
1.65R2
eRp
p 1+ 6
7/3
1/3
(238)
where Rp stands for the radius of curvature of the elbow, R stands for the internal radius of the pipe, and e stands for the thickness of the pipe. In addition, the following conditions must hold: Rp/R < 1.7, centerline length Rp ? > 2 R, and there are no flanges or stiffeners on the elbow. The presented equation for the flexibility factor without part in square brackets is the most commonly used short form; the part in square brackets is a correction for the strengthening effect of the internal pressure.
Figure 64 (upper left) shows the flexibility factor along the piping system for the rod impact experiment. The flexibility factor is different than at the elbow (or any other curved sections of the piping system). The introduction of the flexibility factor influences the eigensystem of the Jacobian matrix. Although the eigenvalues are still constant with time and position, the flexibility factor influences the corresponding eigenvectors L of the Timoshenko beam equations (148):
1
1
0.5
0
0.25 0.2
0.15 0.1
0.05
0
0
-
-0.05 -0.1
-0.15 -0.2
-0.25
LTIM =
11
23
Length [m]
1 0.8 0.6
0.5 0.4 0.3 0.2 0.1 0 -0.1
-0.2
-0.3 -0.4
(239)
0 2 4 6 8 10 12 14 16 18 20
Time [ms]
0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20
Time [ms] Time [ms]
Fig. 64: Flexibility factor profile and histories of some basic variables: comparison with and without
consideration of the flexibility factor.
Figure 64 shows that the history of various variables fits better to the measurement if the flexibility of the elbow is taken into account. The reduction of the flexibility at the elbow decelerates traveling waves in comparison to the simulation without flexibility factor. The ability to take into account different flexibility factors in the same computational section is one of the advantages of the Godunov numerical method. DeYong [29] and Kannapan [69] also utilized flexibility factors to account effects of the elbow ovalization during their simulations of the FSI in piping systems.
k
1.2
3
2.5
2
0.4
4
5
- 97 -
Reinforcement of the elbow. The ASME Boiling and Pressure Vessel code (NB-3641) recommends a minimum wall thickness prior to the bending. For sharp elbows (Rp < 3R) the code suggests an increase of the thickness for 25% compared to the straight section. Introduction of this recommendation into the computational model causes the system to become linear with variable coefficients in space. Figure 65 shows the pipe wall profile with a thicker pipe wall at the elbow and the corresponding eigenvalues profile. Figure 66 shows that the different thickness affects the history of all variables and that the influence on the pressure history is smaller than the influence on the momentum history. The elbow of the experimental apparatus in the rod impact experiment was not reinforced, thus the simulations obtained with the reinforced pipe at the elbow are less accurate compared to the measurement.
2
2.5 2
1.5 1
0.5
0
-0.5
-1
-1.5
4610 1820
Length [m] Length [m]
Fig. 65: Pipe wall profile with thicker wall at elbow and eigenvalues profile.
8 10 12
Time [ms]
14 16 18 20
0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5
0
8 10 12 14 16 18 20
Time [ms]
Fig. 66: Pressure and momentum histories for simulation with normal and reinforced pipe wall.
Water properties. The water properties are evaluated from the water tables in special subroutine, with temperature and pressure as inputs to the subroutine, and density and specific internal energy as outputs of the subroutine. The water properties are evaluated for each computational volume at each time step. Figure 67 shows that the temperature of the fluid affects the pressure and momentum histories because the temperature through the water properties affects the characteristic velocity (speed of sound) of the pressure waves. The influence of the temperature on the basic structural variables (velocity of the pipe, internal forces, momentum, etc.) is indirect, and therefore less evident. Table 20 shows that the speed of sound in a six kelvins colder fluid is 20 m/s lower compared to the speed of sound at room temperature. Small changes in density and temperature have considerable influence on the speed of sound. The speed of sound is fundamental parameter in FSI simulations, therefore it is recommended to perform simulations with exact water properties, and when unavoidable, to be cautious and exact when defining constant water properties in order to estimate the speed of sound as close to reality as possible. Exact temperature should be measured and reported as an important initial condition.
The temperature of the fluid in the (isothermal) eight equation planar quasi-two-phase flow FSI physical model is constant during the simulation (no energy balance equation included), and the experiments presented so far show, that this assumption is accurate for all transients at room temperature. The
1380
1370
1360 4620
4615
15
800
1780
10
5
4587
0
4586
0
2
4
6
2
4
6
- 98 -
temperature changes due to thermal relaxation in such transients are actually small and duration of the transients is so short that exchanges of heat through the pipe walls are negligible.
Table 20: State properties of the fluid (water).
Parameter Constant values T = 293 K T = 287 K
Density [kg/m3] 999.00 999.12 1000.16
Bulk modulus [GPa] 2.1400E9 2.1887E9 2.1401E9
Speed of sound [m/s] 1463.6 1480.1 1462.8
2.5 2
15 1
0.5 0
0.5 0.4 0.3 0.2 0.1 0 -0.1
-0.2
-0.3 -0.4
8 10 12 14 16 18 20
Time [ms]
0 2 4 6
8 10 12 14 16 18 20
Time [ms]
Fig. 67: Pressure and momentum time histories for fluid with variable initial temperature.
Tijsseling et al. [115] reported that they conducted the rod impact experiments at room temperature (it is assumed that room temperature is 7 = 293 K). Tijsseling [113] performed also numerical simulation of the experiment with constant water properties (MOC method), and got a very good agreement between simulation and experiment. The same experiment was simulated with our model, once with Tijsseling’s constant fluid properties and then with accurate water properties. The simulations yield the following conclusions:
• The agreement between the simulation, the experiment and the results of Tijsseling is very good (Fig. 62).
• The results of simulation obtained with accurate and constant water properties are almost identical (Fig. 57, linear-cons. vs. linear-var.). The improvement of real water properties over constant water properties represents a minor improvement for simulations of single-phase transients in cold water. However, it becomes very important for two-phase flow transients or transients in warm and hot water.
• The best agreement with the experiment was obtained for the fluid temperature T = 287 K. This temperature was then applied for all simulations of considered rod impact experiments. The applied temperature is six kelvins lower than the presupposed ‘room’ temperature. However, the decision is justified by the accuracy of the results and by the fact that exact temperature of the water in the rod impact experimental apparatus is not available.
Von Mises stress. The Von Mises stress or simply the Mises stress is a scalar function of the deviatoric components of the stress tensor that gives an appreciation of the overall magnitude of the shear components of the tensor. This allows the onset and amount of plastic deformation under triaxial loading to be predicted from the results of a simple uniaxial tensile test. It is most applicable to ductile materials. In three-dimension, the Mises stress can be expressed as:
CT = '
o-H-tr, + ö", -o-, + \o,-G<
2
(240)
where g, 02, and cr3 are the principal stresses. In one-dimension, this reduces to the uniaxial stress. In terms of a local coordinate system, the Von Mises stress can be expressed as:
2
4 6
- 99 -
°v = J2^°XX ~ °yy ^2 + ^yy ~ °zz ^2 + ^zz " °xx ^2 + 6 (°^ + ^2 + ^) (241)
Von Mises yield criterion for the onset of yield in ductile materials was first formulated by Maxwell [41] in 1865 but is generally attributed to Von Mises in 1913. Von Mises yield criterion can be interpreted physically in terms of the maximum distortion strain energy, which states that yielding in three-dimension occurs when the distortion strain energy reaches that required for yielding in uniaxial loading. Mathematically, this is expressed as: av < ay. In the two-dimensional stress space (shell, pipe wall, a3 = 0) shown in Fig. 68, the yield criterion represents the interior of an ellipse. Stress states a and 02 not touching the boundary of the ellipse produce only elastic deformation. Figure 68 shows also Tresca's maximum shear stress criterion (dashed line labeled Maximum shear). This criterion is more conservative than Von Mises's criterion since it lies inside the von Mises ellipse.
Fig. 68: Von Mises yield criterion for two Fig. 69: Von Mises stress scalar function in time-space dimensional stress space. plane for rod impact experiment in straight piping system.
Fig. 70: Von Mises stress in time-space plane in the upper part of the pipe and the cross-section envelope
for rod impact experiment in single elbow piping system.
Figure 70 shows the appearance of the maximal Mises stress in the upper part of the pipe and the envelope of the Mises stress for the whole cross-section of the piping system. It is evident that the critical section of the pipe with the maximal load is the section in the vicinity of the elbow at the beginning of the transient. The maximal stresses are less than '«I|J ji J
i-n J f 1 Y.
2 4 6 8 10 12 14 16 18 20
Time [ms]
0 0
0 2 4 6 8 10 12 14 16
Time [ms]
2 4 6
0.5 0.4 0.3 0.2 0.1 0 -0.1
-0.2
-0.3 -0.4 -0.5
8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16
Time [ms] Time [ms]
Fig. 72: Histories of basic variables (p2F = 0.87 MPa): measurement vs. simulation.
3.5
3.5
3
3
0.5
0.5
0
18 20
0.5
0
0
18 20
3
2.5
2
0.5
0.5
0
0
18 20
1.2
0.4
0.2
0
0
18 20
- 101 -
1 0.8 0.6
-
-15 ftv
Experiment Simulation
K i*.
'mwm^h
2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16
Time [ms] Time [ms]
Fig. 73: Histories of basic variables (p2F = 1.08 MPa): measurement vs. simulation.
-
0.2
0.15
0.1
0.05
0
-0.05
-0.1
-0.15
2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20
Time [ms] Time [ms]
Fig. 74: Histories of basic variables (p2F = 0.69 MPa): measurement vs. simulation.
1.2 1 0.8 0.6 0.4 0.2 0
1 0.8 0.6
-
2 4 6 8 10 12 14 16 18 20
Time [ms]
0.5
0.4
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
-0.5 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20
Time [ms] Time [ms]
Fig. 75: Histories of basic variables (p2F = 0.30 MPa): measurement vs. simulation.
The simulation was performed with the eight-equation nonlinear model for FSI in a planar arbitrarily shaped piping system with consideration of real water properties and quasi-two-phase flow model. The available experimental data have been applied for comparison with simulation in Figures 71 to 74. The
10
1.2
5
0
-5
0.4
0.2
0
0
18 20
10
5
0
-5
-15
0
1.2
0.4
0.2
10
5
0
-5
-15
0
- 102 -
overall agreement found between the experimental and numerical results is good in view of the complexity of the phenomena. The magnitude and timing of the extreme pressures, the axial velocities, the extreme axial forces and the extreme bending moments are predicted accurately. Discrepancies between numerical and experimental results are attributed to (i) the use of a simple quasi-two-phase flow model of two-phase flow, which is not able to describe cross-section specific localized cavitation at elbow, (ii) the experimental uncertainty (Tijsseling et al. [113] showed that reproducibility of the experiments is generally good, although discrepancies rise for flows with stronger cavitation), (iii) other standard errors, discussed in previous sections, which come with applied numerical method and physical model.
Although an advanced physical models and an advanced numerical method are applied, the overall agreement between the numerical and experimental results is only slightly better compared to the results obtained by Tijsseling et al. [113]. This points to the fact that Tijsseling et al. in their simplified approach, included in physical model all fundamental terms that are essential for the description of the FSI phenomenon coupled with cavitating flow. However, the results presented in this thesis show that all improvements of our physical model and numerical method, give more accurate results and enable significantly more advanced analyses and simulations. While the approach based on the method of characteristics, already reached its limit, our approach enables the inclusion of numerous additional improvements.
Generation and distribution of vapor. Figure 76 shows the comparison of the vapor volume fraction with (left) and without (right) consideration of the FSI for the experiment with moderate cavitation and analogously Figure 77 for the case with severe cavitation. In the cases with FSI, the amount of generated vapor is slightly larger however, the overall generation and distribution with and without FSI are almost the same.
, ~ - t- ^ >..
1.
v. 1 -•
JX
Time [ms] Pipe length [m] Time [ms] Pipe length [m]
Fig. 76: Vapor volume fraction in time-space plane for moderate (p2F = 1.08 MPa) cavitation, with (left) and without (right) consideration of the FSI effects.
5Y
20 \,
x 10 __ -~ "^ .-1
~" r
1 r " k - ~~ "
3-v v r r
0v - " " " jiit*
- K S,
Time [ms]
Pipe length [m]
Time [ms]
Pipe length [m]
Fig. 77: Vapor volume fraction in time-space plane for severe (p2F = 0.3 MPa) cavitation, with (left) and without (right) consideration of the FSI effects.
x 10
x 10
0, 20
20
15
15
10
5
x 10
20
15
15
10
10
5
- 103 -
Von Mises stress in two-phase flow. Comparison of Figures 70 and 78 show that the severity of the cavitation does not influence the overall stress situation in the piping system. The maximal stresses appear at the same position and the temporal distribution is approximately the same. The moderate cavitation case exhibits only one anticipated effect (higher stress) that comes due to different initial pressure. Because all parameters of the simulation, except the initial pressure, are held constant, the maximal stresses are shifted for contribution of the initial pressure. Consequently, the maximal absolute stresses are encountered in single-phase flow, because of the maximal contribution of the initial pressure. It is possible to conclude that the appearance of the two-phase flow, regardless to the intensity of the cavitation, does not reduce the relative maximal stresses in the piping system. The statement that two-phase flow isolates and prevents FSI is wrong. FSI is strong even in compressible two-phase flow.
Fig. 78: Von Mises stress envelope in time-space plane for moderate (left, p2F = 1.08 MPa) and severe (right, p2F = 0.3 MPa) cavitation.
- 104 -
7. Conclusions
Fluid-Structure Interaction (FSI) during a transient in a fluid is an important phenomenon that affects dynamics and integrity of many piping systems. The FSI coupling is strong in soft piping systems conveying less compressible single or two-phase fluids, and is weak in stiff and fully supported piping systems. The integrity of soft piping systems is jeopardized during the FSI occurrence, while the stiff piping systems are costly and are jeopardized by thermal loads. Exchange of energy between the fluid and the structure during the FSI transient occurrence is thus a very important issue and should be predicted (and controlled) as precisely as possible. There are no simple criteria for the estimation of the FSI coupling importance for a general piping system; therefore, several physical and computational approaches have been developed in the past. The approaches that are based on modeling of the FSI coupling phenomena with a single physical model for the thermo-fluid and structural dynamics, and that are solved with an appropriate numerical method, turned out to be the best way for design of new piping systems, for modifications and improvements of the existing piping systems and for analyses of the past accidents. New and more accurate concepts for derivation of advanced physical models, coupling of physical sub-models, and new more accurate numerical algorithm are founded in this dissertation.
The one dimensional general balance equation of thermo-fluid dynamics is derived in arbitrarily shaped, moving and deformable Lagrangian coordinates. The general balance equation is then applied for the derivation of the mass, momentum and energy balance equations of the two-phase flow. The mass, momentum and energy balance equations are used for the construction of the six-equation two-fluid model of the two-phase flow intended for simulations of the two-phase flow transients with FSI. Analogously, the single-phase flow model and the simplified two-phase flow models namely, the homogeneous equilibrium model and the quasi-two-phase flow model, are derived. The Lagrangian derivation yields several new terms in the balance equations compared to the equivalent equations derived in the standard Eulerian coordinate system or compared to the generally accepted fluid dynamics models applied for simulations of the FSI. These terms are terms for local junction coupling FSI mechanism at curvatures and distributed Poisson coupling FSI mechanism, terms with axial force of the pipe and axial and lateral velocity of the pipe, terms accounting for deformations and curvature of the pipe. Accurate thermodynamic states of the fluid are evaluated from the realistic water properties. It is shown that assuming stiff, straight and fixed pipe, the balance equations naturally degenerate into the standard equations derived in the Eulerian frame.
The one dimensional physical models for axial, rotational, lateral and torsional structural dynamics are derived for arbitrarily shaped piping systems described with a general arc length parameter. Several new terms appear compared to the standard physical models of the beam theory, namely, terms accounting for the junction coupling FSI mechanism at curvatures, and distributed Poisson and friction coupling FSI mechanisms. Considering the assumption of a straight piping system, the models for arbitrarily shaped piping systems degenerate into the standard models for straight piping systems. The equations for rotational and lateral movement of the straight piping systems correspond to the equations of the Timoshenko beam theory. The physical models are valid for phenomena in piping systems where longitudinal and lateral frequencies of vibration prevail over circumferential vibration of the pipe wall (low-frequency or long wavelength assumption). The high-frequency phenomena are negligible for the fluid transient and also for the FSI and are thus neglected.
The physical models for thermo-fluid dynamics are grouped with the physical models for structural dynamics into various models of various complexities. They range from the two-equation single-phase flow model to the eighteen-equation model for two-phase flow FSI simulations in arbitrarily shaped and deformable piping systems lying in a 3D space. The fundamental physical model of this dissertation is the eight-equation quasi-two-phase flow FSI model for arbitrarily shaped piping systems in a 2D plane. The equations that make up these physical models are first order nonlinear partial differential equations. The majority of the physical models are hyperbolic by default i.e. the Jacobian matrix of the physical model is diagonalizable, while in order to assure hyperbolicity of all physical models that imply the six-equation two-fluid model of the two-phase flow, a differential virtual mass term is introduced. Analytical solutions for diagonalization exist only for the simple physical models; numerical diagonalization is performed for all other physical models. The eigensystem at output of the numerical diagonalization is
- 105 -
sorted from the largest eigenvalue towards the smallest, therefore additional sorting of the eigensystem is performed at the end of each diagonalization. The sorting is trivial for physical models containing single-phase flow, quasi-two-phase flow, or homogeneous equilibrium two-phase flow, while it becomes quite complicated for the physical models containing the six-equation two-fluid model of a two-phase flow.
The high resolution characteristic upwind finite difference numerical scheme is applied. The term high resolution applies to methods that are at least second order accurate on smooth solutions and yet give well resolved, non-oscillatory discontinuities by modifying the method in the neighborhood of a discontinuity to a monotone first order method. The numerical scheme for the convective terms is explicit: the results exhibit second order accuracy. The stiff source terms of the Timoshenko beam equations are solved by implicit iterations and the stiff relaxation source terms are integrated with the two-step operator splitting approach. The most important property of the numerical method is that it operates with piecewise constant characteristics of the nonlinear physical models. This enables the introduction of numerous improvements like consideration of realistic water properties, consideration of curvatures of the piping system, consideration of geometry changes (cross-section area, pipe wall thickness, material properties, etc.), consideration of the distributed cavitation, consideration of the external loads and masses, friction terms based on Reynolds number, reinforcement of the pipe wall, etc. All these phenomena cannot be addressed with methods dealing with linear equations like the method of characteristics. The weakest point of the characteristic upwind numerical method is the numerical dissipation, which tends to smoothen the solutions at discontinuities. However, the numerical dissipation can be minimized by an appropriate definition of the slope limiter and by increasing the number of computational volumes per unit of length. The numerical method is proclaimed as effective since the average processor time spent for simulation is few minutes, and the most accurate solutions presented in this dissertation took less than one hour on personal computer Pentium IV with 3.0 GHz processor.
The accuracy of the characteristic upwind numerical method is verified by comparison with the results obtained by the method of characteristics. The physical model used was linearized and simplified to match the same assumptions that are used in standard physical models (linear constant coefficient system, no source terms). The agreement is almost perfect, small differences between the characteristic upwind numerical method and the method of characteristics refer to numerical dissipation.
The analytical verification of the results validated the code for the evaluation of the natural frequency of oscillation for arbitrarily supported and loaded piping systems. The code exhibits high accuracy for the non-standard evaluation (simulation with increased damping) of equilibrium values – deflection, internal forces, bending momentum, and other variables, which are used in static calculations in civil engineering. The analysis of the axial and lateral oscillations points out that intensity of the FSI in the axial direction strongly depends on compressibility of the fluid while the intensity of the FSI in the lateral direction is negligible. The most important parameter for lateral oscillations is the mass of the fluid.
Precursor and successor waves are simulated, their origin, appearance and significance for the FSI transient pipe flow is explained. While precursor waves are rather well known and frequently mentioned in the literature, their counterpart, successor waves, are usually neglected; their existence was not traced in the literature by the author yet. Precursor waves under certain circumstances cause the formation of the Poisson coupling induced cavitation. The phenomenon is explained and discussed.
The quasi two-phase flow model is verified by means of experiments for the case with moderate and severe cavitation in a very stiff piping system. The quasi two-phase flow is utilized for the simulation of the rod impact experiment, where FSI coupling is very strong and isolated from the surroundings. Agreement between the experiment and the simulation is remarkable.
The stiffness of a piping system at elbows is reduced due to ovalization of the elbow. The standard model, based on a flexibility factor, is integrated into the computer code. The accuracy of the simulation is increased for the cases with unsupported and flexible elbow.
The 3D stress tensor is represented with scalar values based on the approach established by Von Mises. The Von Mises stress plotted in the time-space plane is one of the most important engineering results of the FSI coupling simulations. It gives timing and position of the maximal stresses in the pipe. The simulations of the pipe with an elbow show that the position and the height of the maximal stresses
- 106 -
are unpredictable without appropriate FSI simulation, and that the duration of the maximal stresses is short – pulsations.
The good agreement between the numerical simulations and the measurements allows explaining the physical phenomena in the relatively simple experimental setups with the help of extensive computer output. Many lessons can be learned from it. It clarifies the phenomenon, when significant cavitation and pipe motion due to water hammer occur simultaneously. Simulations of the FSI coupling transient pipe flow are not possible without consideration of the local junction coupling mechanism and distributed Poisson coupling mechanism. In view of possible transient pipe flows, the junction coupling can be prevented by fixing elbows and/or valves, and the Poisson coupling can be prevented with thick walled pipes. However, none of the solutions to prevent FSI is reasonable. FSI must be allowed to occur, because exchange of energy between the fluid and the structure is in most of the cases beneficial for the integrity of the piping system. Exchange of energy during the FSI must be controlled, and this is not possible without appropriate computer simulations.
The developed advanced physical models and the applied characteristic upwind numerical method have been compiled into the computer code, which is the most prominent practical outcome of the present dissertation. The code is assigned as verified since it was successfully applied for simulation of numerous phenomena related to the FSI and cavitation, and can be now used for simulations of practical cases. There are two possibilities for further research and development. The first is the application of the six-equation two-fluid model of the two-phase flow into the current eight-equation model. The second possibility is a further extension of the equations for 3D space so that it would be possible to simulate arbitrary structures in 3D. It is necessary to stress, that the abovementioned improvements mainly concern computer programming skills, while the firm theoretical background has been established within the scope of the current work. On the other hand, the lessons learned and the advanced physical models derived can be implemented also into the WAHA code, which will improve the WAHA for simulations of the FSI.
- 107 -
References
[I] Abaqus Documentation, Ver. 6.4-1.
[2] Abbot, M.B., Basco, D.R., 1989, Computational Fluid Dynamics: An introduction for Engineers, John Wiley and Sons, Inc., Harlow, England
[3] Allievi, L., 1903, Teoria generale del moto perturbato dell'acqua ani tubi in pressione, Ann. Soc. Ing. Arch. Ithaliana, (French translation by Alleivi, 1904, Revue de mechanique).
[4] Allievi, L., 1913, Teoria del copo d'ariete, Atti Collegio Ing. Arch, (English translation by Halmos, E.E., 1929, The theory of water hammer, Trans. ASME).
[5] Anderson, A., 2000, Celebrations and challenges - water hammer at the start of the 20th and 21st centuries, Proc. of the 8th International Conference on Pressure Surges, BHR Group, 317-322.
[6] Ansys Documentation
[7] Antes, H., 2003, Fundamental solution and integral equations for Timoshenko beams, Computers and Structures, 81: 383-396.
[8] Arscott, F. M., 1955, The oscillation of a heavy spring, The Mathematical Gazette, 39 328: 126-131.
[9] ASME Boiler & Pressure Vessel Code, 1998, Section III, Division 1, Class 1 components.
[10] Barth, T., Ohlberger, M., 2004, Finite Volume Methods: foundations and analysis, Encyclopedia of Computational Mechanics, John Wiley & Sons, Ltd.
[II] Bereux, F., 1996, Zero-relaxation limit versus operator splitting for two-phase fluid flow
computations, Computer Methods in Applied Mechanics and Engineering, 133: 93-124.
[12] Bergant, A., Simpson, A. R., Tijsseling, A. S., 2006, Water hammer with column separation: A historical review, Journal of Fluids and Structures, 22 2: 135-171.
[13] Bestion, D., 1990, The Physical closure laws in the CATHARE code, Nuclear Engineering and Design, 124, pp. 229-245.
[14] Bettinali, F., Molinaro, P., Ciccotelli, M., Micelotta, A., 1991, Transient analysis in piping networks including fluid-structure interaction and cavitation effects, Trans. of SMiRT11, Paper K35/5: 565-570.
[15] Boure, A., 1975, On a Unified Presentation of the Non-Equilibrium Two-Phase Flow Models, Proc. of ASME Symp., New York.
[16] Brennen, C. E., 1995, Cavitation and Bubble Dynamics, Oxford University Press.
[17] Burell, M. J., Enix, D., Lerchel, E., Miro, J., Techendorff, V., Wolferth, K., The Thermal-Hydraulic Code ATHLET for Analysis of PWR and BWR Systems, Fourth Int. Meeting on Nuclear Reactor Therma-Hydraulics, NURETH-4, Proceedings, Vol 2: 1234-1239.
[18] Burman, E., Sainsaulieu, L., 1995, Numerical analysis of two operator splitting methods for and hyperbolic system of conservation laws with stiff relaxation terms, Computer Methods in Applied Mechanics and Engineering, 128: 291-314.
[19] Burmann, W., 1975, “Water Hammer in Coaxial Pipe Systems”, ASCE Journal of the Hydraulics Division, 101, pp. 699-715.
[20] Carlson, K.E., Riemke, R.A., Rouhani, S.Z., Shumway, R.W., Weaver, W.L., 1990, RELAP5/MOD3 Code Manual,Vol. 1-7, NUREG/CR-5535, EG&G Idaho, Idaho Falls.
[21] Carpenter, R. C., 1893, Experiments on Waterhammer, Trans. of ASME, 15.
[22] Casadei, F., Halleux, J. P., Sala, A., Chille, F., 2001, “Transient Fluid-Structure Interaction Algorithms for Large Industrial Applications“, Computer Methods in Applied Mechanics and Engineering, 190, pp. 3081-3110.
[23] Chalabi, A., Qiu, Y., 2000, “Relaxation Schemes for Hyperbolic Conservation Laws With Stiff Source Terms: Application to Reacting Euler Equations”, Journal of Scientific Computating, 15(4), pp. 395-416.
[24] Chalabi, A., 1999, “Convergence of Relaxation Schemes for Hyperbolic Conservation Laws with Stiff Source Terms”, Mathematics of Computation, 68 (227), pp. 955-970.
[25] Clough, R. W., Penzien, J., 2003, Dynamics of Structures, Third edition, Computers & Structures, Inc.
- 108 -
[26] Coutris, N., 1993, Balance equations for fluid lines, sheets, filaments and membranes, International Journal of Multiphase Flow, 19 4: 611-637.
[27] Cowper, G. R., 1966, The shear coefficient in Timoshenko’s beam theory, Journal of Applied Mechanics, 33: 335-340. Timoshenko beam
[28] Davis, J. L., 1988, Wave Propagation in Solids and Fluids, Springer-Verlag, New York Inc.
[29] De Jong, C. A. F., 1994, “Analysis of Pulsations and Vibrations in Fluid-Filled Pipe Systems”, PhD. Thesis, Eindhoven University of Technology, Department of Mechanical Engineering, Eindhoven, The Netherlands.
[30] Debenedetti, P. G., 1996, Metastable liquids, Concepts and Principles, Princeton University Press.
[31] Delhaye, J. M., Achard, J. L., 1976, On the Averaging Operators Introduced in Two-Phase Flow Modeling, Proc. CSNI Specialist Mtg. on Transient Two-Phase Flow, Toronto.
[32] Delhaye, J. M., Giot, M., Riethmuller, M. L., 1981, Thermo hydraulics of two-phase flow systems for industrial design and nuclear engineering, Hemisphere and McGraw-Hill.
[33] Dodge, W.G., Moore, S.E., 1972, “Stress Indices and Flexibility Factors for Moment Loadings on Elbows and Curved Pipes,” Welding Research Council Bulletin, 179.
[34] Downar-Zapolski, P., Bilicki, Z., Bolle, L., Franco, J., 1996, The non-equilibrium relaxation model for one-dimensional flashing liquid flow, International Journal of Multiphase Flow, 22 3: 473-483.
[35] Drew, D. A., Cheng, L., Lahey, R. T., 1979, The analysis of virtual mass effects in two-phase flow, International Journal of Multiphase Flow, 5: 232-242.
[36] Drew, D., Lahey, R. T., 1979, Application of General Constitutive Principles to the Derivation of the Multidimensional Two-Phase Flow Equations, International Journal of Multiphase Flow, 5 (4): 243-264.
[37] Elansary, A. S., Contractor, D. N., 1994, “Valve Closure: Method for Controlling Transients”, Journal of Pressure Vessel Technology, Transactions of the ASME 116 (4), pp. 437-442.
[38] Erath, W., Nowotny, B., Maetz, J., 1999, “Modeling the FSI Produced by a Water Hammer During Shutdown of High-Pressure Pumps”, Nuclear Engineering and Design, 193, pp. 283-296.
[39 Fan, D., Tijsseling, A. S., 1992 “Fluid-Structure Interaction with Cavitation in Transient Pipe Flows”, Journal of Fluids Engineering, Transactions of the ASME, 114, pp. 286-274.
[40] FLUSTRIN, 1990, Manual, ver. 1, Delft Hydraulics.
[41] Ford, H., Alexander, J. M., 1963, Advanced Mechanics of Materials, Longman, London.
[42] Frizell, J. P., 1898, Pressures Resulting From Changes of Velocity of Water in Pipes, Trans. Am. Soc. Civ. Eng., 39: 1-18.
[43] Gale, J., Tiselj, I., 2004, WAHA (WAter HAmmer) Computer Code, Proc. of the 9th International Conference on Pressure Surges.
[44] Gale, J., Tiselj, I., 2005, “Applicability of the Godunov's Method for Fundamental Four-Equation FSI Model”, Proc. of Int. Conf.: Nuclear Energy for New Europe 2005.
[45] Gale, J., Tiselj, I., 2008, Godunov's method for simulations of fluid-structure interaction in piping systems, ASME Journal of Pressure Vessel Technology (final manuscript approved).
[46] Gallouet, T., Massela, J. M., 1996, “A Rough Godunov Scheme”, C.R.A.S. Paris, 323, pp. 77-84.
[47] Ghidaoui, M. S., Zhao, M., McInnis, D. A., Axworthy, D. H., 2005, A Review of Water Hammer Theory and Practice, Applied Mechanics Reviews, 58: 49-76.
[48] Giot, M., Prasser, H. M., Castrillo, F., Dudlik, A., Ezsol, G., Jeschke, J., Lemonnier, H., Rubbers, A., Tiselj, I., Van Hoove, W., Potapov, S., 2003, “Two-Phase Flow Water Hammer Transients and Induced Loads on Materials and Structures of Nuclear Power Plants (WAHALoads)”, FISA 2003, EU research in reactor safety, pp. 523-528.
[49] Godunov, S. K., 1959, Mat. Sb., 47: p. 271.
[50] Goncalves, P. J. P., Brennan, M. J., Elliot, S. J., 2007, Numerical evaluation of high-order modes of vibration in uniform Euler-Bernoulli beams, Journal of Sound and Vibration, 301: 1035-1039.
- 109 -
9931
[51] Harris, C. M., Crede, C. E., 1976, Shock and vibration handbook, McGraw-Hill book company, 2nd edition.
[52] Harvey, A. H., Peskin, A. P., Klein, S. A., 2004, “NIST Standard Reference Database 10,” NIST/ASME Steam properties, Ver. 2.21, US Department of Commerce.
[53] Heinsbroek, A. G. T. J., Kruisbrink, A. C. H., 1993, “Fluid-Structure Interaction in Non-Rigid Pipe Systems - Large Scale Validation Experiments”, Transactions of SMiRT12, pp. 211-216.
[54] Heinsbroek, A. G. T. J., Tijsseling, A. S., 1993, Fluid-Structure Interaction in non-rigid pipeline systems - a numerical investigation II, Proc. of SMiRT 12, J08/2.
[55] Heinsbroek, A. G. T. J., 1997, “Fluid-Structure Interactions in Non-Rigid Pipe Systems”, Nuclear Engineering and Design, 172, pp. 123-135.
[56] Hetsroni, G, 1982, Handbook of multiphase systems, McGraw Hill.
[57] Hirsch, C., 1988, Numerical Computation of Internal and External Flows, Vol 1-2, John Wiley & Sons.
[58] Holmboe, E. L., Rouleau, W. T., 1967, The effect of viscous shear on transients in liquid lines, ASME Journal of Basic Engineering, 89: 174-180.
[59] Hu, C. K., Phillips, J. W., 1981,”Pulse Propagation in Fluid-Filled Elastic Curved Tubes”, ASME Journal of Pressure Vessel Technology, 103, pp. 43-49.
[60] Hutchinson, J. R., 2001, Shear coefficients for Timoshenko beam theory, ASME Journal of Applied Mechanics, 68: 87-92.
[61] Ishii, M., 1975, Thermodynamics of Two-Phase Flow, Eyrolles, Paris.
[62] Ishii, M., Hibiki, T., 2006, Termo-Fluid Dynamics of Two-Phase Flow, Springer Science+Business Media, Inc.
[63] Jaeger, C., 1933, Theorie Generale du Coup de Belier, Dunod, Paris.
[64] Jin, S., Levermore, D., 1996, “Numerical Schemes for Hyperbolic Conservation Laws with Stiff Relaxation Terms”, Journal of computational physics, 126 (1), pp. 449-467.
[65] Jones, S. E., Wood, D. J., 1972, The effect of axial boundary motion on pressure surge generation, ASME Journal of Basic Engineering, 94: 441-446.
[66] Joukowsky, N. E., 1900, On the hydraulic hammer in water supply pipes, Memories of the Imperial Academy Society of St. Petersburg , Series 8, 9(5), (in German) English translation, partly by o Simin: Proc. Am. Water Works Assoc. 24, pp.341-424: 1-71.
[67] Joung, I., Shin, Y. S., 1985, “A New Model on Transient Wave Propagation in Fluid-Filled Tubes”, ASME Proc. of the Pressure Vessels and Piping Conference, 98(7), pp. 275-281.
[68] Kalkwijk, J. P. T., Kranenburg, C., 1971, Cavitation in horizontal pipelines due to water hammer, ASCE Journal of the Hydraulics Division, 97 No. HY10: 1585-1605.
[69] Kannappan, H., 1986, Introduction to pipe stress analysis, John Willey & Sons.
[70] Kladnik, R., 1991, Visokošolska Fizika (in slovene), 1. del, mehanski in toplotni pojavi, DZS, Ljubljana.
[71] Kovač, M., 2004, Influence of microstructure on development of large deformations in reactor pressure vessel steel, Ph.D. thesis, University of Ljubljana, Faculty of mathematics and physics.
[72] Kranenburg, C., 1972, The effect of free gas on cavitation in piplines induced by water hammer, Proc. of the 1st Int. Conf. on Pressure Surges, BHRA, Canterbury, UK, Paper C4: 41-52.
[73] Kries, J. von, 1883, On the relations between pressure and velocity, which exist in the wave-like motion in elastic tubes, Festschrift der 56. Versammlung deutscher Naturforscher und Arzte (in German).
[74] Kucienska, B., 2004, Friction Relaxation Model for Fast Transient Flows, PhD dissertation.
[75] Lavooij, C. S. W., Tijsseling, A. S., 1991, “Fluid-Structure Interaction in Liquid-Filled Piping Systems“, Journal of Fluids and Structures, 5, pp. 573-595.
[76] Le Veque, R. J., Yee. H.C., 1990, “A Study of Numerical Methods for Hyperbolic Conservation Laws with Stiff Source Terms”, Journal of computational physics, 86(1), pp. 187-210.
- 110 -
[77] Le Veque, R. J., 1992, Numerical Methods for Conservation Laws, Second edition, Birkhauser verlag, Basel.
[78] Le Veque, R. J., 2002, Finite volume method for hyperbolic problems, Cambridge University Press, Cambridge.
[79] Lee, U., Kim, J., 1999, “Dynamics of Branched Pipeline Systems Conveying Internal Unsteady Flow”, ASME Journal of Vibration Acoustics, 121, pp. 114-121.
[80] Lemmon, E. W., McLinden, M. O., Friend, D. G., 2003, Thermophysical Properties of Fluid Systems, NIST Chemistry WebBook, NIST Standard Reference Database , 69.
[81] Lemonnier, H., 2002, An attempt to apply the homogeneous relaxation model to the WAHALoads benchmark tests with interaction with the mechanical structure, CEA-T3.3-D61-200302.
[82] Lemonnier, H., 2002, Two-fluid 1D-averaged model equations for a pipe undergoing arbitrary motions, CEA-T3.3-D61-301002.
[83] Leslie, D. J., Vardy, A. E., 2001, Practical guidelines for fluid-structure interaction in pipelines: a review, Proc. of the 10th international meeting of the work group on the behavior of hydraulic machinery under steady oscillatory conditions. See also web site of PHMSA office for pipeline safety: http://ops.dot.gov/stats/LQ_SUM.HTM
[84] Lubis, A., Boyle, J. T., 2004, “The Pressure Reduction Effect in Smooth Piping Elbows – Revisited”, International Journal of Pressure Vessels and Piping, 81, pp. 119-125.
[85] Macchelli, A., Melchiorri, C., 2005, Modeling and control of the Timoshenko beam. The distributed port Hamiltonian approach, SIAM Journal on Control and Optimization, 43 2: 743-767.
[86] Menabrea, L. F., 1885, Note sur les effects de choc de l'eau dans les conduites, C.R. Hebd. Seances Acad. Sci, 47 Jul-Dec: 221-224.
[87] Mendez-Sanches, R. A., Morales, A., Flores, J., 2005, Experimental check on the accuracy of Timoshenko's beam theory, Journal of Sound and Vibration, 279: 508-512.
[88] Mills, A. F., 1999, Heat Transfer, second edition, Prentice-Hall, Inc.
[89] Moody, F. J., 1990, Introduction to Unsteady Thermo fluid Mechanics, John Wiley & Sons, Inc.
[90] MpCCI, Technical reference, Ver. 3.0.
[91] Muller, W. Ch., 1987, Uncoupled and coupled analysis of a large HDR pipe, Transactions of SMiRT9, Vol. F: 31-36.
[92] Naguleswaran, S., 2002, Vibration of an Euler-Bernoulli beam on elastic end supports and with up to three step changes in cross-section, Int. Journal of Mechanical Sciences, 44: 2541-2555.
[93] Naguleswaran, S., 2004, Transverse vibration of an uniform Euler-Bernoulli beam under linearly varying axial force, Journal of Sound and Vibration, 275: 47-57.
[94] Obradović, P., 1990, “Fluid-Structure Interactions: an Accident which has Demonstrated the Necessity for FSI Analysis”, Transactions of the 15th IAHR Symposium on Hydraulic Machinery and Cavitation, Paper J2.
[95] Paidoussis, M.P., 1998, Fluid-Structure Interactions - Slender structures and axial flow, Volume 1, Academic Press, San Diego.
[96] Paidoussis, M.P., 2004, Fluid-Structure Interactions - Slender structures and axial flow, Volume 2, Elsevier Academic Press, London.
[97] Papalexandris, M. V., Leonard, A., Dimotakis, P. E., 1997, “Unsplit Schemes for Hyperbolic Conservation Laws with Source Terms in One Space Dimension”, Journal of computational physics, 134(1), pp. 31-61.
[98] Parmakian, J., 1955, Water-Hammer Analysis, Prentice-Hall Englewood Cliffs, N.J. (Dover Reprint, 1963).
[99] Regetz, J. D., 1960, An experimental determination of the dynamic response of a long hydraulic line, National Aeronautics and Space Administration, Technical note D-576.
[100] Rich, G., 1951, Hydraulic transients, 1st edition, McGraw-Hill, New York (Dover reprint).
[101] Schwarz, W., 1978, “Waterhammer Calculations Taking into Account the Radial and Longitudinal Displacements of the Pipe Wall”, Ph.D thesis, Institut fur Wasserbau, Universitat Stuttgart (in German).
- 111 -
[102] Sha, W. T., Soo, S. L., 1979, On the effect of p?? term in multiphase mechanics, International Journal of Multiphase Flow, 5 2: 153-158.
[103] Simpson, A. R., 1986, Large water hammer pressures due to column separation in sloping pipes, Ph.D thesis, The University of Michigan, Department of Civil Engineering.
[104] Skalak, R., 1956, “An Extension of the Theory of Waterhammer”, Trans. of the ASME, 78, pp. 105-116.
[105] Stadtke, H., Bestion, D., 2002, State-of-the-Art and Limitations of Present Nuclear Thermal Hydraulic Codes and Models, Eurofastnet, Astar, Deliverable D1.
[106] Streeter, V. L., Lai, C., 1963, Waterhammer Analysis Including Fluid Friction, TransTrans. Am. Soc. Civ. Eng., 128: 1491-1524.
[107] Streeter, V. L., 1983, Transient cavitating pipe flow, ASCE Journal of Hydraulic Engineering, 109, No. HY11: 1408-1423.
[108] Svingen, B., 1996, Fluid Structure Interaction in Piping Systems, PhD thesis, The Norwegian University of Science and Technology, Faculty of Mechanical Engineering.
[109] Tang, T., Teng, Z.-H., Wang, J., 2001, “Convergence Analysis of Relaxation Schemes for Conservation Laws with Stiff Source Terms”, Methods and Applications of Analysis, 8(4), pp. 667-680.
[110] Taylor, S. T., Yau, S. C. B., 2003, Boundary control of a rotating Timoshenko beam, ANZIAM -Australian & New Zealand Industrial and Applied Mathematics Journal, 44 E: 143-184.
[111] Tentarelli, S. C., 1990, Propagation of noise and vibration in complex hydraulic tubing systems, Ph.D thesis, Department of Mechanical Engineering, Lehigh University.
[112] Tijsseling, A. S., Fan, D., 1991, The response of liquid-filled pipes to vapor cavity collapse, Transactions of SMiRT11, 183-188.
[113] Tijsseling, A. S., 1993, Fluid-structure interaction in case of waterhammer with cavitation, Ph.D thesis, Delft University of Technology, Faculty of Civil Engineering.
[114] Tijsseling, A. S., Vardy, A. E., Fan, D., 1996, “Fluid-Structure Interaction and Cavitation in Single-Elbow Pipe System”, Journal of Fluids and Structures, 10, pp. 395-420.
[115] Tijsseling, A. S., 1996, “Fluid-Structure Interaction in Liquid-Filled Pipe Systems: A Review“, Journal of Fluids and Structures, 10, pp. 109-146.
[116] Tijsseling, A. S., Lavooij, C. S. W., 1997, “Water hammer with fluid-structure interaction”, Applied Scientific Research, 47, pp. 273-285.
[117] Tijsseling, A. S., 2003, “Exact Solution of Linear Hyperbolic Four-Equation System in Axial Liquid-Pipe Vibration”, Journal of Fluids and Structures, 18, pp. 179-196.
[118] Tijsseling, A. S., 2007, Water hammer with fluid-structure interaction in thick-walled pipes, Computers and Structures, 85: 844-851.
[119] Timoshenko, S. P., Young, D. Y., 1961, Vibration problems in engineering, 3rd edition, D. van Nostrand, New York, 329-331.
[120] Timoshenko, S. P., 1921, On the correction for shear of the differential equation for trasverse vibrations of prismatic bars, Philosophical Magazine, 41: 744-746.
[121] Tiselj, I., Cizelj, L., 1993, Stresses in the steam generator U tubes during large loss of coolant accident, IJS Report, IJS-DP-6766, Confidential.
[122] Tiselj, I., 1997, Sheme drugega reda natančnosti za dvofluidne modele dvofaznega toka (in Slovene), Second order accurate schemes for two-fluid models of two-phase flow, Ph.D disertation, Faculty of mathematics and physics, University of Ljubljana.
[123] Tiselj, I., Petelin, S., 1997, “Modeling of Two-Phase Flow with Second-Order Accurate Scheme”, Journal of computational physics, 136, pp. 503-521.
[124] Tiselj, I., Petelin, S., 1998, First and Second-order Accurate Schemes for Two-Fluid Models, ASME Journal of Fluids Engineering, 120: 363-368.
[125] Tiselj, I., Horvat, A., 2002, “Accuracy of the Operator Splitting Technique for Two-Phase Flow with Stiff Source Terms”, Proceedings of ASME FEDSM02, 2002 ASME FED Summer Meeting.
- 112 -
[126] Tiselj, I., Černe, G., Horvat, A., Gale, J., Parzer, I., Giot, M., Seynhaeve, J. M., Kucienska, B., Lemonnier, H., 2004, “WAHA3 Code Manual“, WAHALoads project deliverable D10, Ljubljana, Download: http://www.rcp.ijs.si/jgale/WAHA.pdf.
[127] Tiselj, I., Gale, J., 2008 (accepted), Numerical integration of unsteady friction models in pipe flow simulations, Journal of Hydraulic Research.
[128] Toro, E. F., 1999, Rieman Solvers and Numerical Methods for Fluid Dynamics, A practical introduction, 2nd edition, Springer.
[129] Trac-M/FORTRAN90 Theory Manual, 2001, NUREG/CR-6724.
[130] Trapp, J. A., Ransom, V. H., 1982, A Choked-Flow Calculation Criterion for Nonhomogeneous, Nonequilibrium, Two-Phase Flow, International Journal of Multiphase Flow, 8 6: 669-681.
[131] Valentin, R. A., Phillips, J. W., Walker, J. S., 1979, “Reflection and Transmission of Fluid Transients at an Elbow”, Transactions of SMiRT5, Paper B 2-6.
[132] Vardy, A. E., Fan, D., 1986, Water hammer in a closed tube, Proc. of the 5th International Conference on Pressure Surges BHRA, 123-137.
[133] Von Karman, 1911, “Ueber die Formanderung dunnwandiger Rohre”, Zeitschrift des Vereines deutscher Ingenieure, 55(2), pp. 1889-1895.
[134] Walker, J. S., Phillips, J. W., 1977, “Pulse Propagation in Fluid-Filled Tubes”, Journal of Applied Mechanics, 44, pp. 31-35.
[135] Warsi, Z. U. A., 1998, Fluid Dynamics, Theoretical and Computational Approaches, 2nd edition, CRC Press.
[136] Weber, W., 1866, Theorie der durch Wasser oder andere incompressible Flussigkeiten in elastischen Rohren fortgepflanzten Wellen (Theory of waves propagation in water or in other incompressible liquids contained in elastic tubes), Berichte uber die Verhandlungen der Koniglichen Sachsischen Gesellschaft der Wissenschaften zu Leipzig, Germany, Mathematical-Physical Section, 18 (in German): 353-357.
[137] Weber, E.-H., 1850, Ueber die Anwendung die Wellenlehre vom Kreislaufe des Blutes und insbesondere auf die Pulslehre (On the application of wave theory to the circulation of blood and particular on the pulse), Berichte uber die Verhandlungen der Koniglichen Sachsischen Gesellschaft der Wissenschaften zu Leipzig, Germany, Mathematical-Physical Section, 2 (in German): 164-204.
[138] Westinghouse Technology Advanced Manual R-504P, Rev 0296
[139] Weston, E. B., 1885, Description of Some Experiments Made on the Providence, RI Water Works to Ascertain the Force of Water Ram in pipes, Trans. Am. Soc. Civ. Eng., 14: 238.
[140] Wiedermann, A. H., 1982, An elastic-plastic pipe response model for small thickness to diameter ratio pipes in water hammer analysis, ASME-PVP, Fluid Transients and Fluid-Structure Interaction, 64: 116-126.
[141] Wiggert, D. C., Otwell, R. S., Hatfield, F. J., 1985, “The Effect of Elbow Restraint on Pressure Transients”, Journal of Fluids Engineering, Transactions of the ASME, 107, pp. 402-406.
[142] Wiggert, D. C., 1986, “Coupled Transient Flow and Structural Motion in Liquid-Filled Piping Systems: A Survey”, Proc. of the ASME Pressure Vessels in Piping Conference, Paper 86-PVP-4.
[143] Wiggert, D. C., Hatfield, F. J., Stuckenbruck, S., 1987, “Analysis of Liquid and Structural Transients in Piping by the Method of Characteristics”, ASME Journal of Fluids Engineering, 109, pp. 161-165.
[144] Wiggert, D. C., Tijsseling, A. S., 2001, “Fluid transients and fluid-structure interaction in flexible liquid-filled piping”, ASME Applied Mechanical Review, 54(5), pp. 455-481.
[145] Wilkinson, D. H., 1978, “Acoustic and Mechanical Vibrations in Liquid-Filled Pipework Systems”, Proc. of the BNES International Conference on Vibration in Nuclear Plant, pp. 863-878.
[146] Wood, F. M., 1937, The Application of the Heaviside’s Operational Calculus to the Solution of Problems in Waterhammer, Trans. ASME, 59: 707-713.
[147] Wood, D. J., 1968, A study of the response of coupled liquid flow-structural systems subjected to periodic disturbances, ASME Journal of Basic Engineering, 90: 532-540.
- 113 -
[148] Wood, D. J., 1969, Influence of line motion on waterhammer pressures, ASCE Journal of the Hydraulics Division, 95: 941-959.
[149] Wood, D. J., Chao, S. P., 1971, Effect of pipeline junctions on waterhammer surges, ASCE Transportation Engineering Journal, 97: 441-456.
[150] Wylie, E. B., Streeter, V. L., 1978, Fluid Transients, McGraw - Hill International Book Company, New York.
[151] Wylie, E. B., Streeter, V. L., 1978, Column separation in horizontal pipelines, Proc. of the Joint Symp. on Design and Operation of Fluid Machinery, IAHR, ASME,ASCE, Colorado State university, Fort Collins, USA, 1: 3-13.
[152] Young, T., 1802, On the velocity of sound, Journal of the Royal Institution of Great Britain, 1: 214-216.
[153] Young, T., 1808, Hydraulic investigations, subservient to an intended Croonian lecture on the motion of the blood, Philosophical Transactions of the Royal Society (London), 98: 164-186.
[154] Youngdahl, C. K., Kot, C. A., Valentin, R. A., 1980, Pressure transient analysis in piping systems including the effects of plastic deformation and cavitation, ASME Journal of Pressure Vessel Technology, 102: 49-55.
[155] Zhang, L., Tijsseling, A. S., Vardy, A. E., 1999, “FSI Analysis of Liquid-Filled Pipes”, Journal of Sound and Vibration, 224(1), pp. 69-99.
[156] Zhang, X., Huang, S., Wang, Y., 1994, “The FEM of Fluid Structure Interaction in Piping Pressure Transients”, Proc. of the First International Conference on Flow Interaction, pp. 532-535.
- 114 -
Appendix A. Gauss theorem and Leibniz rule
A natural set of parameters to describe a pipeline or the fluid domain inside that pipeline consists in a line such as the neutral fiber of the pipe and a cross-section sliding on this line (see Figure A-1). The sliding section on the line generates a domain that Coutris [26] denominates a fluid filament. The fluid filament is thus defined by the inner space of a moving and deforming pipe of circular cross-section. The derivations of the Gauss theorem and the Leibniz rule for a fluid filament of arbitrary cross-section are based on extensive theoretical discussions by Lemonnier [82], Delhaye [32] and Coutris [26] and exceed the scope and purpose of the present work. The interested reader can find the derivation in the above references. The purpose of this section is to define the Leibniz rule and the Gauss theorem for a fluid filament of arbitrary cross-section.
Y
Z /
X
Fig. A-1: A fluid filament generated by the motion of the section S limited by circle 12, the center of which is P and lies on the curve C.
The Gauss theorem for a filament of arbitrary cross-section. Let us described the divergence theorem of vector calculus more commonly known in the older literature as the Gauss theorem. Let V be the region in space with boundary S. Then the volume integral of the divergence V- B of an arbitrary vector or tensor B over V and the surface integral of B over the boundary S of V are related by:
_[VBd\/ = j"BdS
VS
(A-1)
Coutris [26, Eq. C5] mentioned in a footnote that the limiting form of the Gauss theorem holds for arbitrary cross-sections. Lemonnier [82] proved this statement. Following Coutris, the starting point is the Gauss theorem applied to a finite and arbitrary portion of a filament (see Figure A-1) limited by two end sections S? and S2, and the lateral surface of the filament X The Gauss theorem applied to B on V is given by:
_[ V • BdV = JB • nEdS + _[ B • n1dS + _[ B • n2dS
(A-2)
where n is the unit vector normal to \/ pointing outwards the volume and where the subscripts 1, 2, and X refers respectively to surfaces S1, S2, and X. After derivation given by Lemonnier [82], the limiting form of the Gauss theorem, within the Fresnet frame, yields:
_[ VBAdS
S(P,t)
dS,
J BtdS+ j ^^Ädü.
0 S(p,t) n(p,t) nz nn
(A-3)
where Ä is the geometrical factor, and nn is the normal unit vector to the section perimeter and pointing outwards S. This equation is very similar to that of Delhaye [32] related to fixed and straight pipes. The main difference arises from the geometrical factor Ä = 1 - y/R accounting for the lesser weight of points of the cross-section lying closer to the center of the curvature of the neutral fiber of the pipe.
- 115 -
The Leibniz integral rule for a filament of arbitrary cross-section. The Leibniz integral rule gives a formula for the differentiation of a definite integral whose limits are functions of the differential variable:
^^f(x,t)dx=\m^dx + f{b(t),t)^-f{a(t),t)^ (A-4)
dtia(t) Ja(») dt [ dt [ dt
It is sometimes known as the integration under the integral sign. According to the derivation of Lemonnier [82], the Leibniz integral rule for arbitrarily shaped fluid filaments yields (see also [26, Eq. C7]):
f —AdS = — f fAdS +
S(P,t) m m S(P,t)
U C.t j" ac/s + (U Ct)— j fAdS—— j f(U S-t)dS- j f^^Adl
"s0 s(p,o "s0 s(p,f) "s0 s(p,f) n(p,t) ni'nQ
(A-5)
where t is the vector tangent of the Fresnet frame, U2 is the displacement velocity of the lateral surface
(pipe wall), U S is the displacement velocity of the section S and U C is the velocity vector of a point P
attached to the line C. The derivation of the Leibniz rule can be found in Lemonnier [82]. When compared with the analogous equation of Delhaye [32] for fixed and straight pipes, Equation (A-5) shows many new terms. New are the second and the third term on the left hand side and they are related respectively to the stretching of C and the motion of the center of the cross-section S on C. The new term on the right hand side is the second term, which is relative to the motion of the cross-section S. It is naturally encountered when the balance equations are written on a section translating with respect to the frame of reference.
Identity for the pressure. An identity for the pressure is derived from the Gauss theorem. It is a classical transformation in the two-fluid model (Hetsroni, [56]) to collect terms involving the pressure on the interface and pipe boundary and to substitute them with the average of the pressure gradient on the section. This identity is still valid in the case of a moving and deformable pipe. Its derivation is a consequence of the limiting form of the Gauss theorem, that is formally valid also for tensors. Lemonnier [82] derived and proved the following relationship:
škVPk = ^Lškp t + j AnM+ J -ßnM (A-6)
"s0 A n(p,t) nL ' nkn nk(p,t) nL ' nkn
- 116 -
Appendix B. Derivation rules for vectors in Fresnet frame
This appendix gives an overview of the time and space derivatives of vectors in the Fresnet frame (t,n,b) basis. The results are expressed as functions of the derivatives of the vector components. Let
B be an arbitrary vector in the (t,n,b) basis given by:
B = ßft + ßnn + ßöb (B-1)
Time derivatives. The time derivative of B is given by:
B = ^t + ^n + ^b * n 9b
dt dt dt dt ' dt " dt b dt
moreover, since (t,n,b) is a rigid frame:
t = Qxt, and n = Qxn, and b = Qxb (B-3)
dt dt dt
where Q is the angular velocity of rotation of the (t,n,b) frame. As a result:
B = ^t+^n + ^b + RQxt + ßnQxn + ßhQxb (B-4)
dt dt dt dt ' n b
therefore:
9B 38, t dBn - dBb - -. 5
— = —Lt+—n + —b + QxB (B-5)
dt dt dt dt
The components of the time derivative are then given by:
t.3B=^L + t.(QXB) (B-6)
dt dt v
n.^B=^L + n.(öxB) (B-7)
dt dt v
b,dB = dBL b ,ö B
dt dt v / (B8)
Space derivatives. The balance equations encompass derivatives of vectors with respect to the arc length. By using the components of B it yields:
|B = ^(ßft + ßnn + ßöb) (B-9)
Expanding the terms and using the derivative rules for the vectors of the Fresnet frame:
eftn cfn t b cfbn
— = —, — =-------------, and — = — (B-10)
ds Rp ds Rp Tp ds Tp
- 117 -
the space derivative of B yields:
dB rt*°L + B(n
ds
ds
KRPJ
ds
c
v RP TPJ
ds b
-n
JPJ
(B-11)
Collecting the terms together yields:
dB ds
?Bt Bn
-?s R
+ n
PJ
ds Rp Tp
+ b
ds TpJ
(B-12)
Finally the projections of these derivatives are:
t.dB = dBL_^_ ds ds R
(B-13)
9s 9s Rp 7p
(B-14)
3s 3s T
(B-15)
- 118 -
Appendix C. Thermo-fluid dynamics equations in Eulerian coordinates
This appendix gives a short insight into the standard 1D thermo-fluid dynamics models that are based on the equations derived in Eulerian coordinates. The Eulerian coordinates are fixed in space. Typical single and two-phase flow water hammer models in Eulerian coordinates together with the results of a simple water-hammer experiment are given.
C.1. Single-phase flow equations
The research of single-phase water hammer has a long tradition. Blaise Pascal in the 1600s contributed to some of the initial theory to this field, Bernoulli (1738) and Euler established the general equations of hydrodynamics (see also Menabrea [86], Young [152,153], brothers Weber [137,136], Weston [139], Carpenter [21], and Frizell [42]). Joukowsky [66] in 1898 developed the fundamental equation of the single-phase water hammer that relates the pressure changes Ap to the velocity changes Avk in the fluid:
Ap = pkckAvk
where
Av
f,before
f ,after
(C-1)
where pk is the fluid density and ck is the speed of sound in fluid k. The equation was developed from the momentum jump condition (conservation law) under the special case where the flow velocity is negligible compared to the speed of sound. Equation (C-1) is commonly known as the Joukowsky equation, sometimes Joukowsky-Frizell or the Allievi equation. Although the water hammer equations were fully established by the 1960s (Allievi [3,4], Jaeger [63], Wood [146], Rich [100], Parmakian [98], Streeter and Lai [106], Wylie and Streeter [150], etc.), these equations have since been analyzed, discussed, redefined and elucidated in numerous classical texts [47]. The full single-phase water hammer equations with convective terms, gravity and friction terms are defined as:
(C-2)
1 dp pfcf dt 1 pfcf dx dvf dx
dvf dx dx = p
pfgsin(a)-pffw
wf
vr vr
4R
(C-3)
These equations constitute the fundamental equations for 1D water hammer problems and contain all the physics necessary to model wave propagation in complex single-phase fluid-filled piping systems. The standard single-phase water hammer theory correctly predicts extreme pressures and wave periods, but it usually fails in accurately calculating damping and dispersion of wave fronts [150]. The reason is a number of effects that are not taken into account in the standard theory (dissolved and free air, solidified sediment deposit at the pipe walls, unsteady friction and unsteady minor losses, non-elastic behavior of the pipe wall material, etc.). Another omitted effect is FSI in all possible ways.
Tank - Pipe - Valve
V
/-
v = 0.4 m/s, p = 3.42bar, pf = 997.6 kg/m , K = 2.18E9 Pa St = 2.84E-4 m2, E = 0 and 120 GPa, L = 36 m.
Fig. C-1: Geometry of the simple Tank-Pipe-Valve system. The valve is rapidly closed at time zero.
__ 1
0.8 0.6 0.4 0.2
- Initial pre ssure
0.15 Time [s]
Fig. C-2: Pressure history near the valve.
L
-0.2
0
0.05
0.1
0.2
0.25
0.3
- 119 -
Figure C-2 shows a typical single-phase flow pressure history near the valve in a Tank-Pipe-Valve system due to instantaneous valve closure. Equations (C-2) and (C-3) were applied without source terms, thus FSI, damping, body forces, friction, and two-phase flow effects were excluded. The pressure rise above and drop below the initial pressure equals the Joukowsky equation prediction (C-1). The duration of one pressure peak is exactly t1 = 2 L / cf.
C.2. Two-phase flow equations
Physical models and basics. Ishii and Hibiki [62] stressed that two-phase flow thermo-fluid dynamics is an order of a one magnitude more complicated subject than that of the single-phase flow due to the existence of a moving and deformable interface and its interactions with the phases. Significant efforts have been made in recent years to develop accurate general two-phase formulations, mechanistic models for interfacial transfer and interfacial structures, and computational methods to solve these predictive models. Standard two-phase flow models are classified according to the number of partial differential equations that constitute the model. Table C-1 shows some of the most important models.
Table C-1: Typical two-phase flow models.
Number of eqs. Description Vector of basic variables Notes
2 Single-phase water hammer equations y/={p,vf}
3 Homogeneous Equilibrium Model (HEM). Thermal and mechanical equilibrium. W={Pm,PmVm,Pmem} Theoretically important, interesting for FSI.
Inhomogeneous model without heat transfer yr={a,Vg,Vf)
4 Drift-flux model - one phase in saturation (usually vapor). W={Pm,Pg,PmVm,Pmem}
5 Thermal non-equilibrium, mechanical equilibrium. V={Pg,Pf,PmVm,Pfef,Pgeg} Not used in practice
One phase in saturation, the other in non-equilibrium, mechanical non-equilibrium W={Pg,Pf,PgVg,PfVf,Pmem}
6 Thermal and mechanical non-equilibrium. Both pressures equal. Y={Pf,Pg, PfVf ,PgVg, Pfef , Pgeg] Main model in nuclear thermal-hydraulics computer codes.
7 Thermal and mechanical non-equilibrium. Non-equal pressures X = p2 Y={Pf,Pg, PfVf ,PgVg, Pfef ,Pgeg, X] Two-pressure model
Transport equation for interfacial area concentration X = agf
Transport equation for concentration of non-condensable gas X = Og
8 Multi-field models Annular flow: the same phase is modeled with a separate conservation equation for the liquid film at the wall and a separate equation for the droplets in the core.
Multi-group models Bubbly flow: bubble size spectra divided into various classes. Each class of bubbles treated with a separate balance equation.
The most frequently used six-equation two-fluid model is formulated by considering each phase separately. The model is described in every classical work concerning multiphase flows including Ishii and Hibiki [62], Moody [89], Toro [128], Warsi [135], and also researchers like Lemonnier [82] or Tiselj et al. [126] who also integrated the model into the WAHA (WAterHAmmer) computer code. The six-equation two-fluid model of the two-phase flow is derived from Navier-Stokes's system and is expressed in terms of conservation equations governing the balance of mass, momentum, and energy for each phase. The two-phase flow field equations describing the conservation principles require additional constitutive relations or balance equations (also closure relations) for mass, momentum and energy transfers from the kth-phase to the interface. Closure relations encompass the turbulence effects for momentum and energy lost with averaging as well as for interfacial exchanges for mass, momentum and energy transfer. The interfacial transfer rates can be considered as the product of the interfacial flux and the available interfacial area. In the two-phase flow analysis, the void fraction and distribution of the
- 120 -
interface area represent the two fundamental geometrical parameters, which are closely related to the corresponding flow regime. The 1D application of the two-phase models is used for simulations of transients in piping systems whenever pipe’s radius over pipe’s length ratio is small. The closure relations replace various information that are lost with the averaging made to get the 1D model. The computer power still limits the utilization of the full 3D models for simulations of processes. Ishii and Hibiki [62] classify two-phase flows according to the structure of the interface into the following major groups called flow regimes or patterns:
• Separated flows: Film flow, annular flow, jet flow.
• Mixed or transitional flows; Cap, slug or churn turbulent flow, Bubbly annular flow, Droplet annular flow, Bubbly droplet annular flow.
• Dispersed flows: Bubbly flow, Droplet flow, Particulate flow.
Flow regime maps are coded in computer programs with various levels of complexity: from a very detailed flow regime map applied in the RELAP5 code [20] to less detailed applied in the CATHARE code [13]. The standard method using the flow regime transition criteria and the flow regime-dependent closure relations is limited by the fact that closure relations are mainly based on experimental measurements and are thus valid only for certain (steady state) conditions and certain flow regimes. The uncertainty of such correlations and conditions for the transition between flow regimes is even higher during the fast transients, which are considered in this study.
Typical six-equation two-phase flow model. A typical 1D six-equation two-phase flow model in the field of nuclear pipeline thermo-fluid dynamics is given in this section. The described two-fluid model of the two-phase flow belongs to the WAHA [126] code and it similar to the models of RELAP5 [20], Cathare [13], Trac and Trace [129] computer codes. The basic equations are the mass, momentum and energy balances for the liquid and the vapor, with terms for pipe elasticity and without terms for wall-to-fluid heat transfer. The continuity, momentum and internal energy balance equations used in the WAHA code are (equations are derived in Eulerian coordinates):
5 c d c
—akSPk +—akSPkVk
-/'ST
(C-4)
dt
®kSPk Vk +
d „ 2dp -akSPk vk+akS-
iS-CVM + Spj
dak dx
iSF + SF
k,gx
SFk,t
iSTgVi
(C-5)
dt
akSpk ei
tot,k
dx
ak SPk e tot,kVk
d d
+ p—akS +—akSpvk
SQik - /Sr (/,; + v2k /2) + SvkFkgx (C-6)
Where etotk = ek + v2k/2 is the total specific phasic energy, hk = ek + p/pk is the specific phasic
enthalpy, ak is the phasic volume fraction, rg is the vapor generation rate,; is the phase factor, Qik is the phasic volumetric heat flux, F, is the interface drag force, p, is the interfacial pressure term, Fkt is the phasic wall friction, Fkgx is the phasic gravity force in axial direction, and CVM is the virtual mass term. The p, and CVM terms were added to ensure hyperbolicity of the system of equations, which is required by the numerical method. The equations can be applied for the fluid and vapor phase using the definitions in Table C-2. For more details on this model see the WAHA code manual [126].
Table C-2: Application of the WAHA code equations for fluid and vapor phase or mixture.
Property index k vapor volume fraction ak phase factor i
Fluid f 1 - a 1
Vapor g a -1
Figure C-3 shows the geometry and initial conditions for a water hammer transient that is initiated in steady state flow by rapid valve closure. Figure C-4 shows the single and two-phase pressure histories near the valve that are compared to the experiment. Table C-1 contains a brief description of various two-phase flow models. Two of them, the six-equation model of WAHA code and the three-equation
- 121 -
HEM model, are compared in Figure C-5. Figure C-6 shows the importance of an accurate specification of the flow regime and the corresponding flow regime dependent closure terms in two-phase flow codes.
Tank - Pipe - Valve
^>
\
/-
v = 0.4 m/s, p = 3.42bar, ?f = 997.6 kg/m3, K = 2.18E9 Pa St = 2.84E-4 m2, E = 120 GPa, L = 36 m.
Fig. C-3: Geometry of the Simpson pipe [103] experiment; the valve is rapidly closed.
1.2 1
0.8
0.6 0.4
0.2 - 0 --0.2
Experiment
{f\ — Single-phase ------Two-phase
\ |
^/y ,' ü ! l/r
f ^ x 1
Time [ms]
Fig. C-4: Pressure history near the valve: two-phase code accurately describes transient.
1 0.8 0.6 0.4 - 0.2 "
0 -
fvl ------Experiment ..........Six-equation ------HEM
¦> > J L , /! 1] ''*
0.15 Time [ms]
Fig. C-5: Pressure history near the valve:
comparison between six-equation WAHA code and
HEM model.
1 - 0.8p" 0.6 - 0.4 -
0.2- 0 -
------Experiment
'X'l ...........Dispersed flow corr.
------Stratified flow corr.
"\ 1 1
' ^IJ \ ^
1 / "'"mK /\
0.15 Time [ms]
Fig. C-6: Pressure history near the valve: flow regime dependent closure terms are important.
The transition between the single and two-phase flow in the six-equation two-phase flow model starts when the pressure of the single-phase liquid drops below the saturation pressure or when the vapor volume fraction exceeds the numerical uncertainty criterion. The relaxation source terms i.e. heat, mass and energy transfers are calculated only in two-phase flows. The transition back to single-phase flow starts when the vapor volume fraction becomes zero.
The phenomenon of single to two-phase flow transition represents one of the crucial differences between the 'standard' numerical simulations of single-phase flows (CFD) and the numerical simulations of the multi-phase flows (CMFD - Computational Multi-Fluid Dynamics). The transition from single to two-phase flow emerges inside the particular computational volume as a consequence of the thermo-dynamical instability in the single-phase flow or as a consequence of the convection of two-phase flow from a contiguous volume. The Navier-Stokes equations are sufficient to describe single-phase flow with a three-equation model and, when equipped with appropriate closure relations at the interface, also two-phase flow with (usually) six-equation model. But the problem that arises in transition modeling is that the Navier-Stokes equations are not sufficient to describe the transition from single to two-phase flow, where for example bubbles are being generated in the single-phase liquid due to the cavitation. The applied model describes single-phase flow with six equations where the second, non-existing phase is virtually present with a residual volume fraction (?rez = 10-10). As eigenvalues and eigenvectors are defined also in the volumes with single-phase, this approach automatically solves the problem with convection of the two-phase into the volume previously filled with the single-phase.
C.3. Closure relationships for six-equation two-fluid model
The mass, momentum and internal energy balance equations contain several undefined source terms which are given in this section. The modeling of the interface heat, mass and momentum exchange, wall friction forces, body forces, etc. in various physical models relies on the empirically derived, (usually)
L
0
0.05
0.1
0.2
0.25
0.3
0
0.05
0.1
0.2
0.25
0.3
0
0.05
0.1
0.2
0.25
0.3
- 122 -
nondifferential correlations that are usually flow-regime dependent. Our experiences with transient flow modeling show, that a very simple flow regime map with only one flow regime that corresponds to the dispersed bubbly flow is sufficient for accurate simulations of the FSI transient flow. The closure relations were directly extracted from the WAHA code [126] and are briefly described below. More details are given in the WAHA code manual.
Water properties. Whenever realistic water properties are evaluated, an additional equation of state for each phase is needed. The equation of state for phase k is defined by Eq. (111) and Eq. (D-1). The derivatives on the right hand side of the Eq. (111) are determined by the water property subroutines that are based on the ASME steam tables. Water properties are pre-tabulated and saved in ASCII file for 400 pressures: -95 bar < p < 1000 bar and 500 temperatures: 273 K < T < 1638 K. The equation of state relationships are obtained with a three-point interpolation.
Phase-to-interface mass flux. The mass transfer (vapor generation rate rg) is calculated as:
O + O
r
/7.-/7,
(C-7)
where hk* are specific phasic enthalpies (hk = ek + p /pk) and Qik are fluid-to-interface heat fluxes. The specific enthalpies are defined as:
hf=hf and "g g,sat if rg>0
"f = "f,sat and hl=hg if rg<0
(C-8)
Phase-to-interface heat flux. The fluid-to-interface volumetric heat fluxes are calculated as:
Qjk = Hjk(Ts -Tk) k = f, g
(C-9)
where Hik are the heat transfer coefficients. In WAHA, for the dispersed flow the vapor generation rate is calculated with the Homogeneous Relaxation Model (HRM) proposed by Downar-Zapolski et al. [34], and modified by Lemonnier [81]:
-Prr.
x-x
sat
e
(C-10)
where X is the vapor quality, Xsaf is the saturation quality and 0 is the relaxation time correlation:
'108
X
— and Xsat
Pm
"g,saf ~"f,sal
and 0 = max
-4 -0.257
6.51-10"V,
3.84-10-7«-0540
Ap
sat
-2.24
if p 106Pa
( C-11)
Ap
-1.76
Pc-Ps
if p 106Pa
where pc = 221.2 bar is the critical pressure, Ap = max(1000 Pa, \ps - p\), ps is the saturation pressure and ccw = max(106,«).
The vapor heat transfer coefficient Hig is assumed to be extremely large to bring the vapor in equilibrium at a given pressure (similar in the RELAP5 code for dispersed flows):
- 123 -
H ig = min
106
max(«,10 )
(1 + 77- (100 + 25 77)), max
1, r„W-/V
Ts-Tg
max(«,10 9) where r/ = |max(-2,7s -Tg)\. The fluid volumetric heat transfer coefficient is then deduced as:
(C-12)
Hjf
Hig{Ts-Tg)-rg(h;-h*
(Ts-Tf)
(C-13)
In most of the other flow regimes Hif and Hig are calculated first, and /"follows from the heat fluxes Qif and Qig.
Body forces. The body forces per unit volume due to gravity in a pipe with inclination angle ^and considering void fraction, direction and density of the material are defined as:
fluid k:
pipe:
k,gx
Ft,gx
akpkgsin(y) and
ptgsin{y)
and
k,gy
Ft,gy
akpkgcos(r) P,gcos(r)
(C-14) (C-15)
Phase-to-interface momentum flux. The interface friction force is defined as:
Fi =C,\vr\vr+CVM
v - v
(C-16)
In WAHA, for the dispersed flow regime, the interfacial friction coefficient Ci in the momentum equations is calculated from correlations, which are valid for two-component and/or two-phase flow (similar to the RELAP5 model):
1
Ci = 8pfCDag
where CD is the drag coefficient of the slug and agf is the interfacial area concentration: CD = 24(1 + 0.1 Re075)/Re and a
3.6abub/d0
(C-17)
(C-18)
where Re is the Reynolds number, ObUb is the modified vapor volume fraction, and d0 is the average bubble diameter, that are defined as:
(We ¦ a) (1 - or) . _5 (We a) Re =---------j== and abub = max(min(«,0.5),10 ) and d0=2
MfMa Pfvfg
(C-19)
The term ßf is the liquid viscosity, (l/l/e0.6
The applied virtual mass term does not ensure unconditional hyperbolicity. Tiselj [122] reported that complex eigenvalues may appear for very large relative velocities in comparison to the mixture speed of sound (|vr| > 0.3 cm). However these are extremely rare occasions, not relevant in realistic two-phase flows.
The wall friction. The pressure losses due to the wall friction within a given length of a pipe L are defined by the Darcy-Weisbach equation (Wylie and Streeter [150]):
^P = PA
vr vr
4R
(C-23)
where vr = vk - ux. The wall friction force per unit volume is defined as a product of shear stress and contact surface area: Fkt = rktSx/SL and the force balance on a differential pipe section in steady state
flow is defined as: ApS = rktSx . Then the friction force per unit volume is:
F ,— = Pk'k,t
k,t
SL k k,t 4R
Much research is currently focused on transient friction (Tiselj and Gale [127], Kucienska [74]).
(C-24)
- 125 -
Appendix D. Equation of state
The main properties of the fluid are the thermodynamic functions of state and are defined by the equation of state. Consideration of the accurate thermodynamic state properties of the fluid is essential for accurate simulations of single-phase and two-phase flow transients, therefore, the Navier-Stokes system of equations for each phase is supplemented with an additional closure relationship -corresponding equation of state. Fundamentals are briefly described in this appendix.
D.1. Compressibility and speed of sound
According to the set of applied basic variables in this thesis the density of the fluid k is written as a function of pressure and internal energy pk = pk (p, e^. The differential equation of state is:
fyk
fa) dp + [^\ dek
dp
dek
(D-1)
where ek is the phasic specific internal energy, p is the pressure, pk is the density of the fluid.
The single-phase liquid transient pipe flow is under certain circumstances almost isothermal (AT~0) where also Aek ~ 0. Such circumstances are encountered during transient pipe flow with ‘room’ temperature water. This is generally the case in this work. Then the equation of state simplifies into:
*.|Ll*
(D-2)
The compressibility ß is a measure of the relative volume change of a fluid k as a response to a pressure change:
Pt or i
^UpJrors A I 3P Jr or s
(D-3)
where V is the volume, index T stands for the isothermal compressibility (const. temperature T), and S stands for the adiabatic compressibility (const. entropy S). Figure D-1 shows, that differences between isothermal and adiabatic compressibility are small, but not negligible.
Pressure [MPa]
Temperature [K]
Pressure [MPa]
0 700
Temperature [K]
Fig. D-1: Isothermal (left) and adiabatic (right) bulk modulus for water - NIST water properties [80].
The inverse of the compressibility is the bulk modulus K. The bulk modulus of a fluid essentially measures the resistance of the fluid to uniform compression. It is a thermodynamic quantity therefore the bulk modulus varies with variable thermodynamic state (Fig. D-1). The bulk modulus is defined as the pressure increase needed to affect a given relative decrease in volume.
1.5
1.5
0.5
20
20
- 126 -
Kt or s =„------= Pf\ ^n I
A or s KdPkJt ori
dp)
(D-4)
The adiabatic bulk modulus and the density of a material determine the speed of sound c0 of a fluid k:
0,k
Pk
dp) dpkJs
dp
de
Pk (dpk)
(D-5)
Figure D-2 gives the speed of sound in water in a pressure-temperature diagram. The sharp discontinuity appears near the transition from water into steam phase. The speed of sound in single-phase water or steam flow is well defined, while the speed of sound in two-phase flow is a complex flow function. The speed of sound in two-phase water hammer flows can be obtained by substituting the effective bulk modulus of elasticity Ke and the effective density pe in place of K and pk in the equation for the speed of sound (D-5) in single-phase flow. The effective quantities Ke and pe are obtained by the weighted average of the bulk modulus and density of each component, where the partial volumes are the weights. Tiselj [122] derived the following relationship for the speed of sound in two-phase flow, which is valid for the WAHA code described in Appendix C with consideration of the virtual mass term:
PfPg
2F
apt , 1-a)Pg
2 2
CvmPm+aPf+1-a)Pg Cvmp2m + PfPg
(D-6)
where cg and cf are the (effective) speeds of sound in gas and liquid, respectively.
r—
1450
1400 1350 /^ .--; S^^ pL
1300 1250 '/y ""
1200 1150 if/ lil/ Cop per Rod i mpact St eel
1100 1050 if 1 HI Alum nium T = 290 K T = 300 K T = 350 K
ll , 1 , T = 400 K
Temperature [K] Pressure [MPa]
Fig. D-2: The speed of sound water according to the NIST water properties [80].
100 150 200
Elasticity of the pipe [GPa]
Fig. D-3: The effective speed of sound in water conducted by elastic pipe.
Tijsseling [115] in the standard FSI approach applied the following assumption of compressible material with constant density and bulk modulus. The properties of the fluid are evaluated at the beginning of the transient from the initial state. The temporal derivative of density is replaced by the temporal derivative of the pressure given by relationship (D-2):
fyk
(dpk dp
dp
where
fyk dp
Pk
IS
const.
(D-7)
The theory of water hammer based on the assumption of constant fluid properties is sometimes called the elastic water hammer theory according to the analogy with Hooke’s linear elasticity theory for solids.
20
0
50
250
300
700
- 127 -
D.2. Effective speed of sound
The standard water hammer and the FSI theory consider the correction of the speed of sound in a fluid as a consequence of the pipe elasticity. The effective speed of sound in an elastic conduit is given by [66,150]:
1 = 1 + correction = f^l +^^ (D-8)
cf c 2, (dp)s S dp
where variable S stands for the cross-section area and index S stands for the adiabatic process. The speed of sound in a compressible fluid within a rigid pipe is obtained by setting dS/dp = 0, and the speed of sound in an incompressible fluid within an extremely flexible pipe is obtained by setting {dpf/dp)s=0.
Korteweg [150] in 1878 introduced the fluid properties through the state equation (D-4) and used the elastic theory of continuum mechanics to evaluate the correction term. Korteweg developed the relationship for the effective speed of sound and improved it with terms for the Poisson coupling (axial stresses and inertia are considered) and a term to account for the propagation of stress waves along the pipe:
1 (dpf) 2R pt 2R pt 1 1 2R
^ = \iL\ + — y/pf= — + — w>f=J-L where — = — + — w (D-9)
c 2 f [dpjs EeYyf K EeYyf K' K' K EeY
where iff is the multiplication factor that depends on geometry and support conditions. Wylie and Streeter [150] find the following options for the multiplication factor (see Appendix E for details):
• y/=1 - v/2 for a pipe anchored at its upstream end only.
• y/=1 - v2 for a pipe anchored throughout from axial movement.
• y/=1 for a pipe anchored with expansion joints throughout.
Korteweg indicated that his theory is valid for long wavelengths with respect to the pipe diameter. Figure D-3 shows the effective speed of sound as a function of elasticity modulus E of the piping system for various initial temperatures of the fluid according to Korteweg’s equation.
- 128 -
Appendix E. Hooke’s law of elasticity
Hooke’s law for uniaxial strain. Hooke's law of elasticity is an approximation that states that the amount by which a structure is deformed (the strain) is linearly related to the force causing the deformation (the stress). The materials for which Hooke's law is a useful approximation are known as linear-elastic or "Hookean" materials. The most commonly encountered form of the Hooke law is the spring equation, which relates the force exerted by a spring F to the distance w it is stretched by a spring constant k:
F = kw (E-1)
It can be rewritten in a form valid for Hookean materials under certain loading conditions:
a = Ee (E-2)
where E is the modulus of elasticity, zz
0 0 1 + v Gyz axy.
(E-7)
where ^ are the strains and a,, and ^ are the corresponding stresses depicted in Fig. E-1. The stiffness matrix is equal to the inverse of the compliance matrix, and is given by:
yy
yz
xy
(1 + v)(1-2v)
1-v V V 0 0
V 1-v V 0 0
V V 1-v 0 0
0 0 0 1-2v 0
0 0 0 0 1-2v
0 0 0 0 0
0 1 L XX
0 L yy
0 Lzz
0 0 1-2v yz L ZX xy
(E-8)
<7 =
®xx ^ xy *-*xz *-*yx yy yz °zx <7zy °zz
Fig. E-1: Infinitesimal cube with stress vectors and the corresponding stress tensor.
Application of Hooke’s law for piping systems. The wall of the piping system can be treated as a 2D shell structure ( o 2v
dS = 2SdLw =S—1-v2 )dp-S^dNx (E-18)
- 131 -
Appendix F. Timoshenko’s beam theory
Lateral and rotational motion of piping sections can be described with a Timoshenko beam equation (Taylor and Yau [110], Mendez, Sanches, Morales and Flores [87]), which can be decomposed into a system of four linear first order partial differential equations. Timoshenko's beam theory (TBT) constitutes an improvement over the Euler-Bernoulli theory in that it incorporates shear deformation and rotation inertia effects. These improvements of the Euler-Bernoulli theory make the system of four partial differential equations hyperbolic, which is a fundamental condition for the application of the upwind characteristic numerical method. The importance of the shear deformation and rotatory inertia in the description of the response of beams is well documented (Timoshenko and Young [119]) and first improved theory was given by Timoshenko [120] already in 1921. Antes [7] demonstrated that the Euler-Bernoulli theory is a special case of TBT, and that the difference depends on the geometry and is in general very small for static problems. The four Timoshenko beam equations are (see also Section 3.3):
¦
^-^{ir-^w (F-5)
d\-Eltd\-KGSt{d^ dt 2 ' ds 2 { ds
ltPt2t-Elt2-KGSt rf _^ |=0 (F-6)
After elimination of the rotational displacement the fundamental Timoshenko beam equation gives
Elt d4Wy + ltPt d4Wy mT ds 4 /cGSt dt 4 fy{S,t) If, d2 (ltPt + IC ' (mT kGSu fy{s,t) El d4Wy d2Wy ds2dt2 dt2
mT icGStmT dt2 icGS tmT ds2
(F-7)
The left hand part of the equation consists of four terms having the units of force per mass or acceleration. There are terms involving bending momentum, rotational motion, shear force, and
- 132 -
translational motion on the left side of the equation and terms with axially distributed forces on the right side of the equation, respectively. When the shear and rotational terms in (F-7) are small and disregarded, and when there are no external load, the equation becomes the equation of motion of the Euler-Bernoulli beam (Nanguleswaran [92 and 93], Goncalves, Brennan and Elliot [50]):
El. d w d w = ptSt ds 4 d2
Timoshenko’s shear coefficient. The Timoshenko shear coefficient ris defined as the average shear stress on beam cross-section over the shear stress at the neutral axis. The value depends on the shape of the cross-section, and, as pointed out by Cowper [27], on the material’s Poisson ratio v and on the considered frequency range (only for dynamic problems). Hence, different approximations exist that typically give values between 0.5 and 0.6 for a hollow circular cross-section. Hutchinson [60] defined the precise expression of the Timoshenko shear coefficient valid for pipes:
6R2+b2 2 1 + v2
K =------------------------------------------------------------------t------------------------- (F-9)
7R4 +34R2b2 +7b4 +v 12R4 +48R2b2 +12b4 + v2 4R4 +16R2b2 +4b4
where b = R + e is outer radius of the pipe. For a thick walled pipe Cowper [27] proposed:
6 ( 1+v) (+/772 2 1
*- =-------------------1--------'------------- with m = —— (F-10)
( 7 + 6v) (1 + m2) + ( 20 + 12v) m2 1 +
R
For a thin walled pipe m tends to 1 and a simpler and most frequently used definition for the shear coefficient independent of geometry is given by [27]:
2(1 + v)
r =------- (F-11)
4 + 3v
The definitions of the Timoshenko shear coefficient are valid for long wavelengths and low frequencies, since they are based on a quasi-static shear-stress distribution. For short wavelengths, the Timoshenko shear coefficient k becomes frequency dependent [113].
- 133 -
Appendix G. Sklopitev linijske konstrukcije in dvofaznega toka tekočine med hitrimi prehodnimi pojavi - povzetek v slovenskem jeziku
G.1. Uvod
ITK - Interakcija med tekočino in konstrukcijo (ang.: FSI – Fluid-Structure Interaction) je splošen pojem za opis izmenjave (kinetične) energije med gibajočo tekočino in deformabilno konstrukcijo. Količina izmenjane energije je v veliki meri odvisna od gibkosti in deformabilnosti konstrukcije in upora, ki ga le ta nudi toku tekočine (nanj vplivajo geometrijske in materialne lastnosti konstrukcije, podpore, zunanje obremenitve, itd.). Prav tako je količina izmenjane energije odvisna od tekočine, t.j. od gradientov in amplitude tlačnih valov, hitrostnega polja in od lastnosti stanja tekočine, predvsem stisljivosti. ITK se pojavlja v konstrukcijah, ki so potopljene v tekočino, v konstrukcijah v katerih se pretaka tekočina in v kombinaciji obeh vrst konstrukcij. Področje, ki ga zajema ITK je zelo obširno, od aeronavtike, gradbeništva, (procesnega) strojništva, kemične in naftne industrije, in mnogih drugih področij, do končno glasbil in človeškega telesa.
V disertaciji smo se osredotočili na majhen del področja ITK: na dolge (1D) tanke votle konstrukcije v katerih se pretaka tekočina, t.j. cevovode. Cevovodi in cevni sistemi omogočajo transport velikemu spektru tekočin, od vode, kemikalij, nafte, plina, ipd., hkrati pa v nekaterih primerih zagotavljajo varnostno funkcijo (na primer v jedrski elektrarni – hlajenje sredice reaktorja). Tlačni sunki in mehanske vibracije cevovoda med hitrim prehodnim pojavom lahko zmanjšajo varnost in kvaliteto funkcije, ki jo opravlja cevovod. Odpoved cevovoda ima lahko katastrofalne posledice, lahko povzroči poškodbe in smrti ljudi, ekonomsko škodo, okoljske posledice, odpoved ali vibracije cevovoda in hrup. Tijsseling [113] pravi, da ITK ni splošno priznan pojav, čeprav znatno prispeva k določenim odpovedim cevovodov, kot sta na primer utrujanje ali korozija cevovoda.
Raziskave ITK v cevovodih so relativno mlada veja znanosti, še posebno dvosmerna sklopitev konstrukcije in tekočine in raziskave z dvofazno tekočino, čeprav prve hidravlične raziskave vodnega udara segajo daleč nazaj v 18 stoletje (najbrž pa še mnogo dlje, saj so markantne vodne sisteme gradili že Egipčani in Rimljani). Svingen [108] je klasificiral delo, ki se izvaja na področju ITK in delo, ki ga je še potrebno opraviti, na naslednje 4 točke:
• Eksperimentalne preiskave, ki omogočajo širjenje baze podatkov in tako omogočajo izdelavo vedno boljših računalniških programov in povečujejo poznavanje in znanje o pojavu ITK.
• Izpeljava fizikalnih modelov in uporaba naprednih numeričnih shem prevedenih v učinkovite in zanesljive računalniške programe za raziskovalno in praktično inženirsko delo.
• Priprava smernic temelječih na eksperimentih, terenskih izkušnjah in računalniških programih, ki omogočajo prepoznavanje pomembnosti ITK in izdelavo točnejših računskih analiz.
• Določitev projektantskih kriterijev za cevovode, s katerimi se je v večini primerov mogoče na preprost način izogniti nepotrebnim posledicam ITK med hitrimi prehodnimi pojavi.
Raziskovalno delo predstavljeno v pričujoči disertaciji se neposredno nanaša na drugo točko, izpeljava naprednih fizikalnih modelov in vpeljava zmogljivejše numerične sheme. Prva točka (eksperimentalne raziskave) je uporabljena za validacijo računalniškega programa, medtem ko tretja in četrta točka logično sledita iz analize rezultatov.
V nadaljevanju povzetka v slovenskem jeziku, so na kratko navedeni osnovni koraki in najpomembnejši rezultati in dosežki, ki so sicer podrobneje opisani v angleškem delu besedila.
Namen in cilji disertacije. Disertacija je orientirana k modeliranju interakcije med tekočino in konstrukcijo med eno ali dvofaznim hitrim prehodnim pojavom v cevovodu. Obsežne raziskave in publikacije prof. Tijsselinga [112 - 118] predstavljajo tako imenovani stat-of-the-art v svetu modeliranja ITK. Skoraj izključno se uporabljajo linearni hiperbolični sistemi parcialnih diferencialnih enačb prvega
- 134 -
reda, ki v enem modelu opisujejo tako tekočino kot konstrukcijo. Tijsseling in ostala skupnost raziskovalcev uporablja izključno metodo karakteristik (MK). Z MK se lahko rešuje le linearne modele s konstantnimi karakterističnimi hitrostmi, kar že v osnovi izključuje uporabo nelinearnih fizikalnih modelov oziroma izboljšav obstoječih fizikalnih modelov z nelinearnimi členi. Identifikacija omejitev metode karakteristik sugerira osnovni namen pričujoče disertacije, ki je, vpeljati in preveriti numerično metodo, ki omogoča reševanje nelinearnih sistemov diferencialnih enačb. Na podlagi bogatih izkušenj, pridobljenih v času razvoja računalniškega programa WAHA [126] (namenjen je simulacijam termo-hidrodinamike v tekočini med hitrim prehodnim pojavom v popolnoma togem cevovodu brez ITK), in dejstva da se take numerične metode v zadnjem desetletju pogosto uspešno uporabljajo za reševanje inženirskih problemov [44,45,125,127], smo uporabili karakteristično privetrno shemo visoke resolucije drugega reda natančnosti. Drugi pomembni cilji disertacije so (i) izpeljati nelinearne ravnovesne enačbe za tekočino v Lagrangejevem koordinatnem sistemu, (ii) izpeljati enačbe za opis dinamike poljubno ukrivljenega cevovoda, (iii) izpeljati ustrezne fizikalne modele za simulacije ITK, in (iv) izdelati učinkovit in zanesljiv računalniški program za simulacije ITK.
Originalni prispevek. Pričujoča disertacija vsebuje več pomembnih originalnih prispevkov k področju modeliranja interakcije med tekočino in konstrukcijo v cevnih sistemih napolnjenih z eno ali dvofazno tekočino med hitrim prehodnim pojavom. Rezultati raziskav so objavljeni v reviji Journal of Pressure Vessel Technology (Gale in Tiselj [45]). Naše raziskave so privedle do naslednji originalnih prispevkov:
1. Izpeljani so ustrezni nelinearni več-enačbni fizikalni modeli za točnejše simulacije pojava ITK med dvofaznim hitrim prehodnim pojavom v poljubno ukrivljenih in deformabilnih cevovodih.
o Splošna ravnovesna enačba tekočine je izpeljana v premičnem, deformabilnem in poljubno ukrivljenem Lagrangejevem koordinatnem sistemu in uporabljena za izpeljavo ravnovesne enačbe za maso, gibalno količino in energijo. Iz teh enačb so sestavljeni tudi osnovni modeli (termo) dinamike tekočine: šest-enačbni model dvofaznega toka, tri-enačbni model homogenega ravnovesnega toka (HEM model) in dvo-enačbni model kvazi-dvofaznega toka.
o Izpeljan je fizikalni model za opis osne, prečne, rotacijske in torzijske dinamike poljubno ukrivljenega cevovoda napolnjenega s tekočino.
o Izpeljani so različno kompleksni fizikalni modeli za opis ITK, sestavljeni iz fizikalnih modelov za tekočino in konstrukcijo.
2. Uporabljena in preverjena je karakteristična privetrna numerična shema visoke resolucije, ki temelji na metodi končnih razlik, ter izpeljana numerična rešitev za stabilno integracijo togih členov v enačbah Timoshenkovega nosilca.
3. Izpeljan in preverjen je računalniški program za simulacije hitrih prehodnih pojavov v dvofazni tekočini v cevovodu ob upoštevanju prispevka ITK. Različne izboljšave in modeli so integrirani v program in uspešno uporabljeni v simulacijah (točne lastnosti vode, faktor upogljivosti kolena, debelostenski model, zunanje sile, obremenitve, elastične podpore, napetosti po von Misesu, model tanka, ventila, ipd.).
G.2. Ravnovesne enačbe v Lagrangejevem koordinatnem sistemu
Standardne enačbe za opis termodinamike tekočin so izpeljane v Eulerjevem koordinatnem sistemu (dodatek Appendix C). Eulerjeve koordinate so fiksne v prostoru, medtem ko so Lagrangejeve koordinate fiksirane na določen odsek premikajočega se in deformabilnega cevovoda. Lemonnier [82] pravi, da je izpeljava ravnovesnih enačb iz Eulerjevega k.s. v Lagrangejev k.s. obvezna za točen opis termodinamike tekočine v premikajočih se in deformabilnih cevovodih, še posebej, če so premiki cevovoda znatni. Najpomembnejši korak pri transformaciji enačb v Lagrangejev k.s. je povprečenje ravnovesnih enačb preko prečnega preseka cevovoda, ter nato izpeljava in ustrezna uporaba Gausovega (divergenčnega) teorema in Leibnizovega integralskega pravila. Vpeljan je nov splošni parameter s s katerim je enoznačno opisana poljubna krivulja v prostoru. Krivulja predstavlja nevtralno os cevovoda oz. delčka tekočine, ki je predstavljen na sliki G-1.
- 135 -
Y ^
Z
X
Slika G-1: Delec tekočine, ki ga omejujeta dva prečna prereza S? in S2 ter ploskev Z, ki nastane, če
ploskev S z obodom i2 in centralno točko P, drsi po krivulji C. Gibanje oboda i2 generira ploskev Z, ki
predstavlja steno cevovoda. Krivuljo C opišemo s splošnim parametrom krožnega loka s.
Splošna ohranitvena enačba. Vsaka posamezna faza oziroma enofazni delec tekočine omejen z jasnimi mejami oziroma stenami cevovoda, se matematično opiše s sistemom Navier-Stokesovih enačb, ki vsebuje kontinuitetno, gibalno in energijsko enačbo. Navier-Stokesove enačbe so podrobneje obravnavane v mnogih knjigah (Moody [89], Ishii in Hibiki [62], Davis [28], Warsi [135], ipd.). Uporabili smo ravnovesne enačbe v najsplošnejši obliki (Ishii in Hibiki [62]):
dt
PkYk + V • {pkvkVk ) + V-Jk-Pk0k = 0
(G-1)
kjer je pk gostota tekočine k (k = f za kapljevine in k = g za pare), vk je vektor hitrosti tekočine, ostale
spremenljivke pa so definirane v tabeli G-1. S pomočjo tabele G-1 lahko izpeljemo ustrezno kontinuitetno, gibalno in energijsko enačbo iz splošne ravnovesne enačbe.
Tabela G-1: Spremenljivke v splošni ohranitveni enačbi.
Ravnovesje Yk Jk A
Mase 1 0 0
Gibalne količine v* pI-V, g
Totalne energije etot,k=ek + 0.5v2k qk+(pI-Vk)-vk Fg-vk
Splošna 1D ohranitvena enačba v Lagrangejevem koordinatne sistemu. Po izpeljavi, podrobneje opisani v angleškem delu, dobimo naslednjo splošno ravnovesno enačbo v Lagrangejevih koordinatah:
— Sakpky/k + ux— Sakpky/k +^i—Sakpkvky/k+^i—SakJkt-^ Sakpky/k + pi imkWk + J ¦ nk) + pk (J ¦ nk) - S^kPA = 0
(G-2)
Splošna ravnovesna enačba v Lagrangejevih enačbah se lahko uporabi za poljubno tekočino k, pravila za aplikacijo na posamezno fazo pa so definirana v tabeli G-1.
Tabela G-2: Uporaba splošne ravnovesne enačbe v Lagrangejevih koordinatah za poljubno fazo k.
Faza Indeks k Volumski delež pare oik Fazni faktor i
Kapljevina f 1 - a 1
Para g a -1
Kontinuitetna, gibalna in energijska ohranitvena enačba. Če v splošno ohranitveno enačbo vstavimo relacije y/k =1; Jk = 0; 0). Splošna kontinuitetna enačba, gibalna enačba in energijska enačba se poenostavijo:
dt S(XkPk + ds
S, in mT—---------y =-----+— + mTg cos (v)—— (G-10)
^z dt ds Rp Tp V) As
. dt-t (G-14)
Ob predpostavki, da je cevovod raven, odpadejo vsi členi v izvirih, ki imajo v delitelju radij ukrivljenosti cevi Rp ali radij torzijske ukrivljenosti cevi Tp, saj se v takih primerih oba parametra približujeta neskončnosti. Enačbe (G-9) do (G-14) se poenostavijo in postanejo identične standardnim enačbam, ki se uporabljajo v simulacijah ITK (Wiggert, Hatfield in Stuckenbruck [143]). Primerjava pokaže, da so členi z radijem ukrivljenosti in torzijskim radijem ukrivljenosti novi, predstavljajo pa mehanizem vozliščne sklopitve. Simulacije ITK v poljubno ukrivljenih cevovodih niso mogoče brez teh členov. To pomanjkljivost standardnih enačb rešujejo tako, da cevovod razrežejo na odsekoma ravne odseke, kolena in vozliščno sklopitev pa opišejo z dodatnimi enačbami (Tijsseling [115]). Enačbe dinamike cevovoda so uporabne tudi za opis nihanja cevovoda (Lavooij in Tijsseling [75]), izračun lastne frekvence in statičnih sil in premikov cevovoda (z in brez polnila - vode).
G.4. Fizikalni modeli za opis ITK
Za opis ITK lahko uporabimo različne modele za različno kompleksne pojave. Osnovni in tukaj največkrat uporabljeni model za simulacije eksperimentov je osem-enačbni kvazi-dvofazni model ITK za
- 138 -
ravninski cevovod, ki je v nadaljevanju podrobneje predstavljen. Najbolj splošni model sestavlja 18 enačb z naslednjimi osnovnimi spremenljivkami (Tiselj [122] ali Hirsch [57]):
• 6 spremenljivk za tekočino: hitrost vsake faze vg in vf, tlak p, volumski delež pare «in specifična notranja energija vsake faze eg and ef.
• 12 spremenljivk za konstrukcijo: osna hitrost cevi ux in osna sila Nx, prečni hitrosti cevi uy in uz in prečni sili Qy in Qz, rotacijski hitrosti j-1/2
Shema v enačbi (G-28) je 1. reda natančnosti po kraju in času in je stabilna, če je izpolnjen CFL pogoj:
At<
rmax
(G-31)
Karakteristična privetrna shema visoke resolucije. Osnovna ideja privetrnih shem visoke resolucije je, da so rezultati na gladkih rešitvah najmanj drugega reda natančnosti, medtem ko je v območju nezveznih (strmih) rešitev, uporabljena metoda prvega reda z večjo numerično difuzijo (LeVeque [77,78], Hirsch [57], Toro [128]). Po izpeljavi, ki presega okvir tega povzetka, dobimo karakteristično privetrno shemo visoke resolucije, ki je v osnovi enaka enačbi (G-28). Spremenijo se le korekcijski faktorji definirani v enačbi (G-30):
fp+=max
0,A
PlAx
in
fp-=min
k
PlAx
(G-32)
Prvi člen v korekcijskih faktorjih je že znan smerni faktor, drugi člen pa je popravek drugega reda natančnosti določen z omejitvenim faktorjem /2
yy;
+ o-,
yy
+ U>V
+6(4+4+4
(243)
Če je napetost po von Misesu manjša od meje elastičnosti (?v < ?y), potem lahko predpostavimo elastične deformacije. Na primer, v 2D konstrukcijah, kamor prištevamo tudi cevovode, mora biti napetostno stanje določeno s ?1 in ?2 znotraj elipse (slika G-8 levo). Slika kaže tudi konzervativnejši
- 146 -
Tresca kriterij (črtkana črta). Slika G-8 desno kaže položaj in časovni potek napetosti po Misesu (ovojnica največjih napetosti). Kritični del konstrukcije pri obravnavanem prehodnem pojavu je koleno na začetku, ko skozi koleno potuje prvi napetostni val. Največje napetosti so ?v,max = 50 MPa, medtem ko je tipična meja elastičnosti v jeklu ?y = 250 MPa. Trajanje največjih napetosti je kratko – pulziranje.
Dvofazni eksperimenti. Če začetni tlak v cevovodu znižamo, med hitrim prehodnim pojavom pride do kavitacije. Kavitacija je lahko zmerna ali močna, odvisno od začetnega tlaka v sistemu (slike od G-9 do G-11). Glede na uporabljen preprost kvazi-dvofazen model dvofaznega toka in kompleksnost pojava lahko zaključimo, da je ujemanje med meritvijo in eksperimentom zelo dobro. Magnituda in časovni potek največjih napetosti, tlakov in kavitacije sta opisana zelo dobro. Opaziti je možno, da ujemanje nekoliko upada z naraščanjem intenzivnosti kavitacije. Razloge za neujemanje lahko pripišemo (i) uporabi poenostavljenega model dvofaznega toka, (ii) eksperimentalni nezanesljivosti (Tijsseling [113] je pokazal da je ponovljivost relativno dobra, pada pa s povečevanjem intenzivnosti kavitacije), (iii) drugim napakam, kot je na primer numerična disipacija, nepopolni fizikalni modeli, ipd.
3.5 3
2.5 2
1.5 1
0.5
0.5
0.4
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
-0.5 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16
Cas [ms] Cas [ms]
Slika G-9: Časovni potek nekaterih (izbor) osnovnih spremenljivk (p2F = 1.24 MPa).
1 0.8
10
5
0 -5
-15 t|f
Eksperiment Simulacija
2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16
Cas [ms] Cas [ms]
Slika G-10: Časovni potek nekaterih (izbor) osnovnih spremenljivk (p2F = 1.08 MPa).
0 -5
0.5
0.4
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
-0.5 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16
Cas [ms] Cas [ms]
Slika G-11: Časovni potek nekaterih (izbor) osnovnih spremenljivk (p2F = 0.30 MPa).
0
0
18 20
1.2
0.6
0.4
-10
0.2
0
0
18 20
10
5
-10
-15
0
18 20
- 147 -
Volumski delež pare. Slika G-12 kaže volumski delež pare za primer močne kavitacije z in brez upoštevanja ITK mehanizmov. Vidimo, da ITK sicer vpliva na porazdelitev pare vendar minimalno.
10 - - ' r L - "~ " r r u r
10 ^3 ¦ ..—-'¦¦
0 0
Cas [ms]
Slika G-12: Volumski delež pare za primer močne kavitacije (p2F = 03 MPa). Levo z upoštevanjem učinkov
ITK in desno brez upoštevanja.
Von Mises napetosti v dvofaznem toku. Slika G-13 kaže primerjavo napetosti za primer z zmerno in močno kavitacijo, rezultate pa primerjamo tudi s sliko G-8 desno kjer ni kavitacije. Najpomembnejša ugotovitev je, da kavitacija ne vpliva bistveno na napetosti oziroma na sam potek ITK in hitrega prehodnega pojava. Opaziti je le, da so največje napetosti dosežene v enofaznem primeru, najnižje maksimalne napetosti pa v dvofaznem primeru z močno kavitacijo. Razlog je v tem, da je v enofaznem primeru začetni tlak v cevi največji, kar seveda povzroči največje napetosti v stenah cevovoda. Sama kavitacija pa na ITK ne vpliva bistveno.
Slika G-13: Ovojnica napetosti po von Misesu za primer zmerne (levo, p2F = 1.08 MPa) in močne (desno, p2F = 0.3 MPa) kavitacije.
G.7. Zaključek
Obravnavali smo interakcijo med tekočino in konstrukcijo (ITK) med hitrimi prehodnimi pojavi v tekočini. Ugotovili smo, da je ITK intenzivna v mehkih in slabo podprtih cevovodih v katerih se pretaka skoraj nestisljiva tekočina (eno ali dvofazna). Uporaba togega in preveč podprtega cevovoda je dokaj preprosta rešitev, toda v takih cevovodih, še posebno v jedrskih elektrarnah, bi nastopile velike notranje sile zaradi termičnih obremenitev. Izmenjava kinetične energije med tekočino in cevovodom je v osnovi ugoden pojav, ki pa ga je potrebno ustrezno kontrolirati. Kontrola ITK med hitrimi prehodnimi pojavi ni možna brez ustreznih orodij za numerične simulacije. Za opis ITK v cevnih sistemih so najprimernejši enodimenzionalni sistemi parcialnih diferencialnih enačb, ki vsebujejo tako enačbe za tekočino kot enačbe za konstrukcijo, in ki so nato numerično rešeni z enotno numerično proceduro. V tej disertaciji so na novo izpeljane diferencialne enačbe za tekočino in konstrukcijo, uporabljena je nova numerična metoda, vse skupaj pa smo prevedli v računalniški program.
- 148 -
r~" "
^ - - ' ^--1
" r
r '
_ - - " ~~ j
i- -
i- - ~~ ~~
_.
r "
x 10
0 20
fe^
Dolžina cevi [m]
Cas [ms]
Dolžina cevi [m]
15
10
3
3
Izpeljali smo enodimenzionalno splošno ohranitveno enačbo termo-hidrodinamike v poljubno ukrivljenem, premikajočem in deformabilnen cevovodu (Lagrangejev koordinatni sistem). Iz splošne ohranitvene enačbe smo izpeljali kontinuitetno, gibalno in energijsko ravnovesno enačbo. V enačbah se pojavijo novi členi in nove spremenljivke, s katerimi je upoštevana dinamika in ukrivljenost cevovoda. Iz splošnih ohranitvenih enačb je izpeljanih več sistemov enačb, ki so primerni za opis hitrih prehodnih pojavov v različnih tekočinah (izotermični eno in kvazi-dvofazni model, šest-enačbni dvofazni model, tri-enačbni dvofazni model toplotnega in mehanskega ravnovesja). Izpeljali smo tudi 1D fizikalne modele osne, rotacijske, prečne in torzijske dinamike poljubno ukrivljenega cevovoda. Pojavi se več novih členov, vsi pa predstavljajo mehanizem vozliščne sklopitve v kolenih. Različni modeli termo-hidrodinamike in različni fizikalni modeli dinamike cevovoda so nato združeni v več sistemov enačb za opis pojava ITK v cevovodih različne kompleksnosti. Vsi sistemi so sestavljeni iz parcialnih diferencialnih enačb prvega reda, so hiperbolični, lastne vrednosti je možno določiti analitično ali pa numerično. Osnovni fizikalni model ITK v disertaciji je osem-enačbni model sklopitve izotermičnega kvazi-dvofaznega toka in ravninskega cevovoda.
Enačbe smo reševali s karakteristično privetrno shemo končnih razlik visoke resolucije in drugega reda natančnosti. Numerična shema je eksplicitna, občasno pa uporabimo tudi implicitne iteracije, ker so členi izvirov togi in bi bil potreben časovni korak integracije nesprejemljivo majhen. Najpomembnejša prednost uporabljene numerične metode je možnost reševanja nelinearnih sistemov oziroma sistemov enačb s spremenljivimi parametri. Tako je možno uporabiti enačbe stanja tekočine, ki ustrezajo dejanskemu termodinamičnemu stanju tekočine, upoštevanje geometrijskih in materialnih sprememb v cevovodu (prečni prerez, debelina cevi, ukrivljenost cevi, drugačna togost cevi ipd.), možnost upoštevanja dvofaznega toka na poljubnih delih cevovoda, upoštevanje zunanjih obremenitev, togih ali elastičnih podpor, ipd. Numerična metoda vnaša numerično disipacijo, ki pa jo je možno relativno enostavno nadzirati in zmanjševati. Točnost numerične metode smo preverili z uveljavljeno metodo karakteristik, ujemanje rezultatov je zelo dobro, malenkostna odstopanja v nezveznostih lahko pripišemo numerični disipaciji in drugačnemu konceptu računanja.
Rezultate smo preverjali tudi z analitičnimi enačbami, pokazali smo, da je z našim pristopom možno določiti lastno frekvenco nihanja poljubno ukrivljenega cevovoda, prav tako je z veliko natančnostjo možno določiti klasične spremenljivke statike konstrukcij: notranje sile, premike, momente, ipd. Ugotovili smo, da na osno sklopitev cevovoda in tekočine v veliki meri vpliva stisljivost tekočine, medtem ko na prečno sklopitev vpliva predvsem masa tekočine. Pokazali smo obstoj valov imenovanih znanilci in nasledniki, prav tako smo pokazali, da valovi znanilci v določenih primerih lahko preko Poissonovega mehanizma sklopitve povzročijo izredno kavitacijo v sistemu. Izpeljali in preverili smo poenostavljen model dvofaznega toka. S tako imenovanim kvazi-dvofaznim modelom smo uspešno in natančno simulirali več dvofaznih eksperimentov. Vgradili in preverili smo model zmanjšanja togosti kolena zaradi pojava ovalizacije. Uporabili smo tudi debelostenski fizikalni model in pokazali, da upoštevanje debeline cevi za simulacijo praktičnih primerov, ni bistveno. Iz inženirskega stališča je najpomembnejši rezultat napetost izračunana po metodi von Misesa. Največje obremenitve cevovoda so kratkotrajne, pojavljajo pa se na različnih in večinoma v naprej nepredvidljivih delih cevovoda. Mesto in čas največje obremenitve cevovoda je možno napovedati samo z ustrezno računalniško simulacijo.
Izpeljani napredni fizikalni modeli in uporabljena karakteristična privetrna numerična shema so bili prevedeni v računalniški program, ki smo ga obširno in uspešno testirali z obstoječimi eksperimentalnimi in numeričnimi podatki. Vidimo dve možnosti za nadaljevanje raziskav in razvoj programa. Prva je uporaba šest-enačbnega dvotekočinskega modela dvofaznega toka tekočine, druga pa uporaba vseh enačb dinamike konstrukcije, s katerimi bo mogoče opisati dinamiko prostorskega cevovoda. Poudariti je potrebno, da se omenjeni izboljšavi v veliki meri nanašata na programersko delo in iskanje ustreznih računalniških rešitev, medtem ko so teoretične osnove že podane v disertaciji. Obstaja tudi možnost, da se s fizikalnimi modeli in izkušnjami pridobljenimi med izvajanjem raziskav, nadgradi obstoječi termo-hidrodinamični program WAHA, tako, da ga bo mogoče uporabiti tudi za simulacije prehodnih pojavov z upoštevanje interakcije med tekočino in konstrukcijo.
- 149 -
IZJAVA O AVTORSTVU
Spodaj podpisani Janez Gale izjavljam, da predložena disertacija predstavlja rezultate lastnega znastveno raziskovalnega dela.
Janez Gale