Mechanical Interaction in Pressurized Pipe Systems: Experiments and Numerical Models

Abstract

The dynamic interaction between the unsteady flow occurrence and the resulting vibration of the pipe are analyzed based on experiments and numerical models. Waterhammer, structural dynamic and fluid–structure interaction (FSI) are the main subjects dealt with in this study. Firstly, a 1D model is developed based on the method of characteristics (MOC) using specific damping coefficients for initial components associated with rheological pipe material behavior, structural and fluid deformation, and type of anchored structural supports. Secondly a 3D coupled complex model based on Computational Fluid Dynamics (CFD), using a Finite Element Method (FEM), is also applied to predict and distinguish the FSI events. Herein, a specific hydrodynamic model of viscosity to replicate the operation of a valve was also developed to minimize the number of mesh elements and the complexity of the system. The importance of integrated analysis of fluid–structure interaction, especially in non-rigidity anchored pipe systems, is equally emphasized. The developed models are validated through experimental tests.

1. Introduction

Pipe systems have frequently changes in flow or pressure, induced by a routine maneuver of hydromechanical equipment, leading to a hydraulic transient event. If the change is gradual, the resulting pressure transient may not be severe; however, if the change is sudden, the transient pressure can cause significant damage to the whole system [1]. The most common causes of rapid flow changes in a pipe system are typically by pump or turbine power failures, pipe breaks or valve maneuvers by opening or closure. However, they can also result from natural causes, equipment malfunction, or even operator error. In addition, resonant vibrations can damage poorly anchored pipes if they are not properly fixed. Pressure transients are common when operating a pipe system, where the challenge lies in keeping the pressure transients and associated phenomena within controllable limits. The damage caused by unexpected violent surge pressures is not always predictable. Often, the consequences, such as leakage occurrence, pipe ruptures and joints displacements, do not become apparent until long after the event. Most of the time, the damage caused by waterhammer generally surpasses the cost of surge control and preventive measures [2].

Therefore, numerical models allow the hydraulic engineer analyzing a potential hydraulic transient event, in order to identify and evaluate alternative solutions for predicting and controlling possible extreme pressures that may occur in the pipe system [3]. In this way, a numerical 1D model based on the method of characteristic (MOC) is used and compared with experimental tests taken in a laboratory network. However, the behavior of various components in these complex systems is typically highly nonlinear, in the geometry, due to large motions and deformations, and in the material rheological behavior laws and pipe support displacements [4]. Studying the behavior of the fluid and the structure, as a whole, enables the response of the structure change by the flow modification and, consequently, on return, the response of the fluid by structure vibration [5]. Thus, in order to simulate fast transient events, appropriate numerical modeling for interactions between different elements is required. In most complex events, simulations based on CFD analyses are used to better understand the dynamic phenomena. When dealing with a joined problem, two domains are created, the fluid-flow and the pipe structure. Hence, specific analyses are developed in order to show possible consequences in terms of displacements when the system is under different hydraulic loads and types of pipe supports. The importance of integrated analysis (fluid–structure interaction) and the role of supports in water infrastructures design is equally emphasized.

In conventional waterhammer 1D analyses, the pipe elasticity is incorporated into the wave speed of the pressure waves [6]. Pipe inertia and axial pipe motion are not taken into account for rigidly anchored pipe systems (e.g., buried pipes). For less restrained systems (e.g., pipe bridges and long pipes), fluid–structure interaction may become very important. In this case, the dynamic behavior of the fluid and the pipe system are treated simultaneously [7].

Moreover, calculations with fluid–structure interaction are suitable in cases with both domains influences, in the development of design criteria and operating rules [8,9,10], being specifically useful in post-accident analyses and interpretations. The results show that the new computational techniques are efficient and can yield accurate evaluation of the FSI in pressurized systems.

2. Experimental Study

Laboratory Facility

The laboratory facility is comprised of an 85 m PVC pipe, a ball valve, and a returned open flow channel. The pipe is braced, through metal clamps, on iron supports approximately 0.5 m in height. The connections are flexible and made through Gibault joints. The maximum discharge is about 25 L/s. From downstream to upstream, 2 pressure transducers (PT2 and PT1) are installed throughout the facility (Figure 1).

Figure 1. Scheme of the laboratory facility.

Figure 1. Scheme of the laboratory facility.

In the configuration for waterhammer tests, the pipe is connected upstream to a constant head of 4 m. Two pressure transducers are installed along the entire pipe system: pressure transducer 1 (PT1) located 2 m from the FCV (Section 1) and pressure transducer 2 (PT2) located 1 m from Section 3, the operating valve. These transducers are connected to an oscilloscope with a bandwidth of 20 MHz for repetitive waveforms, a vertical amplifier frequency response of 5 mV to 5 V with an accuracy of ±3%–5% and a horizontal amplifier frequency response DC (Direct Current) to 2 MHz with an accuracy of ±3%–6%. Accelerometer equipment is used to register the pipe vibrations over time and a manometer.

After creating a steady regime for a flow of 9 L/s, two of the iron supports are removed (in Section 2) and a surge is induced by a quickly closing of the ball valve (Section 3). Repeated tests are performed for the same conditions to verify the consistency of results. Figure 2 shows the data acquisition recorded in the pressure transducers, in the flow meters (located upstream of the PVC pipe), and also the structure vibration in each axis direction (i.e., x, y, z).

Figure 2. Data acquisition in the laboratory facility.

Figure 2. Data acquisition in the laboratory facility.

The ball valve at Section 3 takes 0.13 s to close completely. After 0.27 s, the pressure variation recorded in PT2 is observed in PT1 (time the wave takes to cross the PVC pipe length until it reaches PT1). Simultaneously, the accelerometer records the pipe vibrations during the transient event.

The pressure wave attains the upstream boundary and reflects to the source perturbation (valve maneuver) in Section 3, T_E = 0.4 s. A sudden change in pressure occurs in PT2, leading to an overpressure after 0.93 s (Figure 2).

Figure 3 shows the first two main overpressures in PT2 and the acceleration induced in the structure. As soon as the ball valve closes (Section 3), i.e., after 0.13 s, the structure starts to move in three directions due to the internal pressure variation created by the transient event. During T_E, the pipe reaches the maximum displacement in the x-axis of 3.14 m/s² (i.e., in longitudinal direction of the pipe) and 2.2 m/s² in the y-axis. In the z-axis (vertical), the structure vibrates but practically maintains in plan XY, due to the roller supports that allow movements in this plan. Apart from the upstream (i.e., Section 1) and downstream (i.e., ball valve in Section 3) sections of the system, that are fixed to concrete supports, all the remaining supports (i.e., roller) allow movements only in the pipe direction (i.e., x or y, depending on the pipe orientation). Thus, the restrictions imposed by the roller supports change according to the direction of the pipes. After T_E, the structure moves significantly in y-axis, reaching its maximum displacement (2.72 m/s²) at t = 0.67 s. Compared to the pressure variation, the accelerometer is affected immediately by the wave pressure coming from the ball valve towards the downstream end of the PVC pipe. At t = T_E, a symmetry can be noticed, for about 0.1 s, in x-axis, coinciding with the instant the pressure wave reaches the downstream end and travels downstream towards the ball valve (elastic reflection time). During the relief wave, the pipe represents a harmonic behavior in y-axis until the end of time.

Figure 3. Pressure variation (m) and pipe vibration (in acceleration units).

Figure 3. Pressure variation (m) and pipe vibration (in acceleration units).

According to Figure 3 and the known frequency of the pressure wave, it is possible to estimate the pipe vibration in terms of displacements. Figure 4 shows the total displacement calculated and the respective maximum displacement achieved in the laboratory facility.

Figure 4. Pipe displacements over time.

Figure 4. Pipe displacements over time.

3. Methods

3.1. Modified Method of Characteristic Model with Non-Elastic Effects

Mathematical models are extensively used as part of routine design projects, as well as for management and control behavior of hydraulic systems, in particular in the prediction of dynamic verification under normal, incidental and emergency operating conditions. In the solution of 1D hydraulic transients, MOC has been used because of its easy programming and efficient results for long pipelines. Considering the partial differential equations, based on continuity and momentum equations:

$$\begin{array}{c}{c}^2\frac{\partial Q}{\partial x}+gS\frac{\partial H}{\partial t}=0\\ \frac{\partial Q}{\partial t}+gS\frac{\partial H}{\partial x}+RQ\left|Q\right|=0\end{array}$$

(1)

where $R=f/2DS$ is the friction term, $f$ the Darcy–Weisbach friction factor, $D$ the inner pipe diameter, $S$ is the cross-sectional area of pipe, $Q$ is the volumetric flow, $H$ is the piezometric head, $c$ is the wave speed, $g$ is the gravitational acceleration, $x$ is the space coordinate, and t is the time.

A linear combination of Equations (1), can be used,

$$\left(\frac{\partial Q}{\partial t}+\lambda c^2\frac{\partial Q}{\partial x}\right)+\lambda gS\left(\frac{\partial H}{\partial t}+\frac{1}{\lambda }\frac{\partial H}{\partial x}\right)+RQ\left|Q\right|=0$$

(2)

If $H=H\left(x,t\right)$ and $Q=Q\left(x,t\right)$, then the total derivatives are expressed by

$$\frac{dQ}{dt}=\frac{\partial Q}{\partial t}+\frac{\partial Q}{\partial x}\frac{dx}{dt}$$

(3)

and

$$\frac{dH}{dt}=\frac{\partial H}{\partial t}+\frac{\partial H}{\partial x}\frac{dx}{dt}$$

(4)

Defining the unknown multiplier $\text{λ}$ as

$$\frac{1}{\text{λ}}=\frac{dx}{dt}=\lambda c^2$$

(5)

with

$$\text{λ}=\pm \frac{1}{c}$$

(6)

and using Equations (3) and (4), Equation (2) can be written as

$$\frac{dQ}{dt}+\frac{gS}{c}\frac{dH}{dt}+RQ\left|Q\right|=0$$

(7)

if

$$\frac{dx}{dt}=c$$

(8)

and

$$\frac{dQ}{dt}-\frac{gS}{c}\frac{dH}{dt}+RQ\left|Q\right|=0$$

(9)

if

$$\frac{dx}{dt}=-c$$

(10)

Equations (7) and (9) are valid when Equations (8) and (10) are also valid [11]. Imposing Equations (8) and (10), the one independent variable, $x$, is eliminated, converting the partial differential equation (Equation (1)) into ordinary differential equations in the independent variable $t$. In the x–t plane, Equations (8) and (10) represent two straight lines with slopes $\pm 1/c$, called characteristic lines. As shown in Figure 5, Equations (7) and (9) represent the C⁺ and C⁻ compatibility relations, respectively.

Figure 5. Characteristic lines in x–t plane.

Figure 5. Characteristic lines in x–t plane.

Multiplying the left-hand side of Equation (7) by $dt$ and integrating generates Equation (11):

$${{\displaystyle \int }}_L^{P}dQ+\frac{gS}{c}{{\displaystyle \int }}_L^{P}dH+R{{\displaystyle \int }}_L^{P}Q\left|Q\right|dt=0$$

(11)

where the use of subscripts L and P indicate the locations in the x–t plane. Since Equation (11) is only valid along the characteristic line LP, the first two integral terms are evaluated easily. However, for the third term, representing the friction losses, a second order approximation of the integral was used, since it yields satisfactory and stable results for engineering applications [12,13,14,15,16,17]:

$$R{{\displaystyle \int }}_L^{P}Q\left|Q\right|dt=R\Delta t\left|Q_L\right|Q_P$$

(12)

Hence Equation (11) becomes

$$Q_P-Q_L+\left(H_P-H_L\right)\frac{gS}{c}+R\Delta t\left|Q_L\right|Q_P=0$$

(13)

and for Equation (9)

$$Q_P-Q_R-\left(H_P-H_R\right)\frac{gS}{c}+R\Delta t\left|Q_R\right|Q_P=0$$

(14)

Thus, the compatibility equations for the positive $C^{+}$ and negative $C^{-}$ characteristic lines can be written as follows [12,18]:

$$\begin{array}{c}{C}^{+}: H_P-H_L+\left(Q_P-Q_L\right)\frac{c}{gS}+R\left|Q_L\right|Q_P=0\\ C^{-}: H_P-H_R-\left(Q_P-Q_R\right)\frac{c}{gS}-R\left|Q_R\right|Q_P=0\\ R=\frac{R^{\prime }}{gS}\Delta x, \Delta x=c\Delta t, R^{\prime }=\frac{f}{2DS}\end{array}$$

(15)

where $H_i$ is the piezometric head at node i, $Q_i$ the volumetric flow at node i, $R$ is the friction term using a second order approximation, $\Delta t$ the time step and $\Delta x$ the space increment.

The characteristics of Equation (15) are solved simultaneously with the conditions imposed by the boundaries. Therefore, an explicit finite-difference method is used, i.e., the partial derivatives are replaced by finite-difference approximations such that the unknown conditions at a point at the end of a time step are expressed in terms of the known conditions at the beginning of the time step [11]. To solve the ordinary differential equations (ODEs), a proper time-integration algorithm is needed [14]. Therefore, combining the MOC model with the conditionally stable Euler forward integration method, results in a system of algebraic equations solved explicitly (i.e., one at a time) [11,16,18]. This simple Euler fixed time step integration algorithm has first-order accuracy [15,19,20]. According to [13], the numerical scheme is stable if $\Delta x\ge c\Delta t$ is verified. This condition is referred to as the Courant–Friedrich–Lewy (CFL) stability condition [13]. The Courant number ($C_r$) is defined as the ratio of the real wave speed, c, and the numerical wave speed ($\Delta x/\Delta t$). Thus, for stability reasons, the computational time step $\Delta t$ and the reach length $\Delta x$ must be selected such that $C_r\le 1$.

According to the generic formulation applicable to any system characteristic, the use of dimensionless parameters of relative head, $h$, relative head losses, $\Delta h_0$, and relative discharge, $q$, is also considered [12,17]:

$$h=\frac{H}{\left(\frac{c Q_0}{gS}\right)}=\frac{H}{\Delta H_J}; \Delta h_0=\frac{\Delta H_0}{\Delta H_J}; q=\frac{Q}{Q_0}$$

(16)

$\Delta H_J$ is the Joukowsky overpressure.

In an elastic behavior, the head variation is given by

$$h=\frac{h_0}{1+h_{0}K\Delta h_0\left(\tau -{\tau }_0\right)}$$

(17)

where $h_0$ is the dimensionless head at time ${\tau }_0$, being $\tau =t/T_E=t/\left(2L/c\right)$.

The modification of the head loss coefficient is obtained by a multiple constant factor. For the description of the non-elastic behavior due to fluid, pipe material (plastic pipe) and supports, two additional parameters (KH and KQ) were included in the MOC equations [17]. Parameter KH stands for the head decay induced by a discharge variation due to the non-elastic behavior of the fluid, the pipe viscoelasticity and resistance of supports. KQ is a raised coefficient in the discharge value that represents the non-elastic response in the recuperation phase of a deformation caused by a head variation, according to [12,20]. Through Equations (16) and (17), the decay coefficients were estimated and included in the MOC as dissipative parameters:

$$\Delta H=KH\frac{c}{gS}\Delta Q-J$$

(18)

$$\Delta Q=KQ\frac{\Delta H-J}{c/\left(gS\right)}$$

(19)

where $J$ is the head loss term, and $\Delta H$ and $\Delta Q$ are the head and discharge variation, respectively.

In characteristic lines, one-sided difference is taken to extract information from the interior of the domain. As for an incoming characteristic line, the characteristic equation is replaced by the boundary condition.

The numerical model was calibrated based on experimental data and benchmark analysis [12,20]. The physical characteristics adapted in the experimental simulations are the ones registered in the pressure transducers. The wave speed is obtained from the physical properties of the pipe material. Table 1 shows characteristic parameters from experimental test and coefficients adjusted for the tested scenario.

Table 1. Physical characteristics of the laboratory facility and model coefficients.

Description	Parameter/Coefficient	Value
Upstream head (m)	H	14.4
Discharge (m3/s)	Q	0.009
Wave speed (m/s)	$c$	350
Time closure of the ball valve (s)	tf	0.20
Total simulation time (s)	tt	3
Head coefficient induced by a discharge variation by non-elastic fluid and pipe deformation (--)	KH	0.32
Discharge coefficient induced by a head variation, due to a non-elastic response in the recuperation phase of the deformation (--)	KQ	2.9

Table 1. Physical characteristics of the laboratory facility and model coefficients.

Description	Parameter/Coefficient	Value
Upstream head (m)	H	14.4
Discharge (m3/s)	Q	0.009
Wave speed (m/s)	$c$	350
Time closure of the ball valve (s)	tf	0.20
Total simulation time (s)	tt	3
Head coefficient induced by a discharge variation by non-elastic fluid and pipe deformation (--)	KH	0.32
Discharge coefficient induced by a head variation, due to a non-elastic response in the recuperation phase of the deformation (--)	KQ	2.9

The celerity is calculated through Equation (20):

$$c=\sqrt{\frac{K/\rho }{1+\left(K/E\right)\left(D/e\right)\psi }}$$

(20)

where K is the volume compressibility modulus; E is the elasticity modulus of the pipe material; D is the internal diameter; e is the thickness; ρ is the fluid density and ψ is the non-dimensional parameter that depends upon the elastic properties of the pipe, calculated for different anchored conditions (Table 2).

Table 2. Non-dimensional parameter associated to the support conditions.

Support Conditions	Non-Dimensional Parameter ($\psi$)
Pipe with frequent expansion joints	$1$
Pipe against longitudinal movement throughout its length	$1-{\text{μ}}^2$
Pipe against longitudinal movement at the upper end	$1-\frac{\text{μ}}{2}$

Table 2. Non-dimensional parameter associated to the support conditions.

Support Conditions	Non-Dimensional Parameter ($\psi$)
Pipe with frequent expansion joints	$1$
Pipe against longitudinal movement throughout its length	$1-{\text{μ}}^2$
Pipe against longitudinal movement at the upper end	$1-\frac{\text{μ}}{2}$

The values obtained for the wave speed for each anchoring condition are presented in Table 3.

Table 3. Pipe characteristics and wave speed values for each anchoring condition.

Material	PVC
Di (m)	0.107
e (m)	0.0035
E (GPa)	2.98
K (N/m2)	2.19 × 109
$\text{ρ}$ (kg/m3)	1000
$\text{υ}$ (--)	0.46
$\text{ψ}$	1	0.79	0.77
c (m/s)	300	350	350
wave speed adopted (m/s)	350

Table 3. Pipe characteristics and wave speed values for each anchoring condition.

Material	PVC
Di (m)	0.107
e (m)	0.0035
E (GPa)	2.98
K (N/m2)	2.19 × 109
$\text{ρ}$ (kg/m3)	1000
$\text{υ}$ (--)	0.46
$\text{ψ}$	1	0.79	0.77
c (m/s)	300	350	350
wave speed adopted (m/s)	350

Hence, the results obtained through 1D model simulations and experimental tests are presented in Figure 6. In the transducer, located near the ball valve (PT2), the maximum overpressure value reached in the second main peak is also achieved by the computational model. However, as soon as the ball valve closes, major discrepancies can be noticed. These differences are due to the pipe displacements and deformations. For the first transducer (PT1), a similar behavior with the MOC model is observed. Herein, it is also verified that the maximum overpressure reached in the computational model is below the experimental one (not common) due to the existence of strong interaction between the fluid and the structure, not adequately taken into account by the MOC.

To achieve the solution of the momentum and continuity equations in the MOC model, the following conditions were presumed, i.e., the fluid is homogeneous and incompressible, the elasticity of the pipeline and fluid follows a linear pattern and the flow is uni-dimensional (1D). Herein, the characteristic and compatibility equations are written for the whole pipe domain, where the entire pipeline is divided into a predefined number of pipe increments (∆x = 5.25 m).

Figure 6. Comparison between method of characteristic (MOC) model results and experimental data.

Figure 6. Comparison between method of characteristic (MOC) model results and experimental data.

3.2. FSI Modeling and Solution

3.2.1. Basic Concepts

In cases of significant interaction between fluid and structure, suitable numerical models are necessary to solve the mass and momentum differential equations [21]. Among different methods, the finite element method (FEM), used in the CFD models, stands out because of its versatility and generality. Its first characteristic is that the continuum field is divided into cells, called elements, which build a mesh. These elements (in 3D) have tetrahedral shape. In FEM, when a structure is submitted to external forces, the objective is to determine in any point of the structure: (i) the displacement μ (vector with 3 components); (ii) the state of deformation ε (tensor with 6 independent components); and (iii) the state of tension σ (tensor with 6 independent components). In each point, it is possible to stabilize 15 equations: displacements and deformations (6 compatibility equations); deformations and tensions (6 constitutive equations); and tensions and mass forces (3 equilibrium equations).

$$\text{Compatibility equations}:\text{ ε}=L u, L-diferential operator$$

(21)

$$\text{Constitutive equations}:\text{ σ}=D \text{ε}, D-elasticity matrix$$

(22)

$$\text{Equilibrium equations}:L^T\text{σ}+f=0, f-mass forces$$

(23)

To summarize, the problem of calculating a structure can be reduced to determining the three components of the displacements at each point [21,22]. Hence, it is necessary to use the Navier equations concerning, in each point, three displacement components, u, with three components of the mass forces, f. The mass forces can be gravitational forces, inertial forces related to seismic accelerations forces and/or pressure waves and pipe deformations (Figure 7).

Figure 7. Navier formulation (adapted from [21,22]).

Figure 7. Navier formulation (adapted from [21,22]).

In practice, the calculation of a structure summed up in solving the differential equations of Navier (24), taking into account the geometry of the structure and the boundary conditions (i.e., the forces applied and its restrictions of the support):

$$L^T\left(DLu\right)+f=0, {\forall }_{\left(x1, x2, x3\right) \in structure}$$

(24)

3.2.2. Geometry and Mesh Adaptation

All calculations have been performed on a PC (Intel 5, CPU 3.90 GHz, RAM 8 GB) with 4 cores and threads running in parallel. Herein, the coupled procedure is taken to solve the flow and the structural equations at the interface, which means the flow changes induce structural deformation and vice versa. For FSI, the fluid is described by an arbitrary Lagrangian–Eulerian (ALE) formulation fully coupled to general solids and structures [23]. The structures can undergo highly nonlinear response due to large deformations and nonlinear material behavior [24]. The FSI problem is solved using a two-way coupling method. Since the structure density is higher than the fluid one and thus the added-mass effect is not significant [25,26,27,28,29], the coupling scheme is stable. In this way, an explicit method can be employed [30]. After discretizing the governing fluid and structural fields, the equilibrium and compatibility conditions must be satisfied at the interface [31,32,33]. The fluid solution is computed using the Navier–Stokes solver on the overlapping domain. The boundary velocity is determined by projecting the fluid velocity to the neighboring overlapping mesh of the structural domain. The structural equations are then solved to yield detailed stress and strain distribution on the solid, which is used to find the resultant force that generates the motion of the structure [24]. The implicit backward differentiation formula (BDF) with adaptive time stepping, which is an extension of the backward Euler method for variable order [33], is used.

In the CFD model, the 3D geometry is based on the Guanajuato lab facility, created in AutoCAD. Inside the model, two types of domains were considered: Fluid (water volume) and Structure (pipe volume, including the wall thickness). In the fluid flow, a turbulent k-ɛ model was selected. For the solid domain, the PVC pipe, and for the supports, some restrictions were imposed, according to their characteristics. Figure 8 represents the 3D model, with the two domains incorporated that are composed of 130,681 and 254,370 tetrahedral elements, respectively, with 188 vertex elements, i.e., 564 degrees of freedom (DOF) and 172,526 triangular elements. The total number of DOF is 1,673,631.

Figure 8. Fluid and structure integrated mesh.

Figure 8. Fluid and structure integrated mesh.

According to the type of supports, apart from the ball valve that is fixed to a concrete support, all the roller supports (Figure 9) allow movements only in the pipe direction, with the exception of two (i.e., first roller support indicated in the legend) that allow movements in every direction, besides z-axis.

Figure 9. Type of supports.

Figure 9. Type of supports.

In the FSI solutions, the ALE formulation is used, to allow reasonable mesh deformations in each solution step [23]. Once the mesh has been significantly distorted, a mesh adaptation is required. The mesh adaptation prevents distorted elements by distributing and sizing them. Appropriate element sizes are firstly estimated and then the decision is made whether to adapt the mesh, according to the quality of mesh achieved. Since the mesh is to be adapted, a new mesh is created and the previous solution is mapped onto that mesh. This mapping also includes the mapping of the boundary conditions and material properties. However, when the fluid element presents a great amount of minimum element sizes, as in boundary layers, and the structural deformations are rather large in an incremental step measured on these fluid element sizes, in such cases, overlapping elements may occur and the fluid element mesh can become distorted. For such problems, the solution is started with a rather coarse finite element mesh, allowing large structural deformations to take place. Then, the fineness of the fluid element mesh is increased. Consequently, since the structure will only adjust its deformations with small changes in the geometry, small fluid element sizes can now be accommodated in the mesh adaptations to obtain adequate accuracy. Finally, the mesh, in its original and final configurations, is represented in Figure 10. After the solution convergence, the maximum displacement achieved in two main directions is 0.087 m.

Figure 10. Mesh in original and final configurations for t = 0.21 s.

Figure 10. Mesh in original and final configurations for t = 0.21 s.

3.2.3. Boundary Conditions

In a 3D model, when simulating a valve maneuver, different methods to describe the opening and closing of the device can be used. It is difficult to model this type of device due to the number of mesh elements and the moving boundary layers, i.e., the part of the geometry that acts as shutter. An interesting alternative is to use a viscous specific element, i.e., a material property that can vary over time. In this case, a very high viscosity (~infinite viscosity) is specified in order to control the flow. To simulate the movement of this specific fluid region with high viscosity, the model includes a logical expression in the domain settings.

Figure 11 shows the model geometry and the operation of the fluid valve. The flow passes through the valve section, where the viscosity wall varies between z₀ and z₁ (see detail in Figure 11).

Figure 11. Geometry and operation of the high viscosity fluid valve.

Figure 11. Geometry and operation of the high viscosity fluid valve.

In the inlet condition, a constant velocity is considered (0.99 m/s) in the moving boundary, according to a function of time and for the outlet condition an open boundary with atmospheric pressure is adopted. The movement of the viscous element is described by an analytic formulation, which returns the value of one in the area corresponding to the boundary (face) valve, and zero elsewhere. Hence, the viscous element is expressed by:

$$={\eta }_0+range·{\eta }_{\infty }$$

(25)

where ${\eta }_0$ is the fluid viscosity, ${\eta }_{\infty }$ is a very high viscosity value (~infinite) and the $range$ is described by:

$$range=range\_x·range\_y·range\_z$$

(26)

where

$$\begin{array}{c}range\_x=(x>x_0)(x<x_1)\\ range\_y=(y>y_0)(y<y_1)\\ range\_z=1-\left(z>z_0\right)+(z>z_1)\end{array}$$

(27)

$x_0, x_1, y_0, y_1, z_0$ are fixed in time and describe the size of the valve boundary and the amplitude with which the boundary moves. $z_1$ is a function that describes the valve boundary amplitude over time, according to:

$$z_1=z_{max}\times int\left(t\right)$$

(28)

where, $z_{max}$ is a constant and corresponds to the valve boundary amplitude, given by the distance between $z_0$ and $z_1$ and $int\left(t\right)$ is a time function that depends on the type of shutter used (i.e., ball valve), given by Figure 12.

Thus, $z_1$ depends on the type of valve given by the respective characteristic curve. Since the shutter is a ball valve, the characteristic curve is represented in Figure 13.

Figure 12. Time function according to the shutter characteristics.

Figure 12. Time function according to the shutter characteristics.

Figure 13. Characteristic curve of the ball valve.

Figure 13. Characteristic curve of the ball valve.

4. Results

4.1. Velocity Profiles

In order to understand the unsteady velocity distribution during the water-hammer wave propagation, three vertical cut lines at three locations along the pipe were chosen to illustrate the numerical results (Figure 14).

Figure 15 shows the velocity distribution according to the cross sections identified in Figure 14, from 0.16 s to 0.5 s, covering the time that the ball valve closes completely. As the valves closes, a reverse in the velocity profile is verified ($0.21 s \le t\le 0.41 s$). After the complete closure of the ball valve, the flow travels to upstream ($t>0.50 s$) (Figure 15).

Figure 14. Identification of the three locations for analysis.

Figure 14. Identification of the three locations for analysis.

Figure 15. Velocity profiles in three pipe locations: upstream, middle, and valve locations.

Figure 15. Velocity profiles in three pipe locations: upstream, middle, and valve locations.

Regarding the velocity profile obtained in the valve section (Figure 16), as the high pressure wave is established at the closing valve (t = 0.21 s), the pressure wave progressively brings the entire velocity profile to be zero. As soon as the valve starts to close, the deceleration of the velocity can be noticed (for 20% opening). However, the process of deceleration is not uniform, showing a near zero mean flow for t = 0.41 s. Figure 16 shows that it is possible to obtain high rates of energy dissipation and shear stress with the inversion of the velocity profile.

According to the RANS, the effect of the boundary layer growth before the closing valve was captured in the CFD model. Through Figure 17, it is possible to identify the region of the flow in which the effect of the viscous shearing forces caused by the fluid viscosity is felt (i.e., the boundary layer). This hypothetical boundary surface divides the flow into two regions: the boundary layer region, in which the viscous effects and the velocity changes are significant, and the irrotational (core) flow region, in which the frictional effects are still negligible and the velocity remains almost constant in the middle of the pipe cross section.

Figure 16. Velocity profile obtained at the valve section (location 1).

Figure 16. Velocity profile obtained at the valve section (location 1).

Figure 17. Development of the velocity boundary layer in the three pipe sections.

Figure 17. Development of the velocity boundary layer in the three pipe sections.

According to Figure 18, before the closure of the valve, the turbulent regime of the flow, due to the disorder velocity vectors, is observed. As soon as the valve starts to close (i.e., t = 0.21 s), the flow inverts its direction, creating vortex near the valve wall. During the closing time, it is verified that the flow deceleration is more significant near the pipe wall by the flow inversion occurrence. After the flow inversion is visible, the vortexes originate a velocity field with a practically null spatial mean velocity. In less than T_E (elastic reflection time), the inside flow is almost stopped and the unsteady friction effect can be neglected. This can also be noticed in Figure 19, where the connection between the closure of the valve and the velocity streamlines are shown. This figure illustrates the velocity vectors in the pipe system from the ball valve to upstream (Section 2 of Figure 14).

Figure 18. Velocity vectors, in the valve section, during the closing time of the shutter.

Figure 18. Velocity vectors, in the valve section, during the closing time of the shutter.

Figure 19. Velocity contours over time.

Figure 19. Velocity contours over time.

4.2. Wave Propagation and Pipe Displacements

The transient pressure wave propagates along the pipe system (Figure 20) at the same time as the structure moves over time (Figure 21). This can be noticed when comparing the pressure variation in the valve section with the internal pressure created by the FSI during the waterhammer event, especially during the first period.

Figure 20. Wave propagation pressure contours on the fluid–structure interface.

Figure 20. Wave propagation pressure contours on the fluid–structure interface.

Figure 21. Fluid velocity (streamlines) and pipe displacement (contours) fields during the first wave period: (a) t = 0.50 s; (b) t = 0.55 s; (c) t = 0.60 s; and (d) t = 0.96 s.

Figure 21. Fluid velocity (streamlines) and pipe displacement (contours) fields during the first wave period: (a) t = 0.50 s; (b) t = 0.55 s; (c) t = 0.60 s; and (d) t = 0.96 s.

Figure 21 represents the velocity streamlines for four instants, i.e., immediately after the valve closes completely. Apart from the streamlines, the structure’s displacement from its initial position, represented by the movable configuration, is verified. The total displacement achieved in the CFD model represents the transference of forces and momentums between the fluid and the pipe wall during the occurrence of the unsteady condition created by the fast closure of the ball valve. This total displacement can be decomposed in deformations in x-axis and y-axis (Figure 22).

Figure 22. Deformed pipe at unsteady states: (a) in x-axis and (b) in y-axis.

Figure 22. Deformed pipe at unsteady states: (a) in x-axis and (b) in y-axis.

4.3. Deformation Gradient

The deformation gradient (F) describes the deformation of a line element in the reference configuration into the current configuration. It is the fundamental measure of deformation in continuum mechanics. It is the second order tensor that maps line elements in the reference configuration into the current configuration. Consider a line element $dX$ emanating from position $X$ in the reference configuration that becomes dx in the current configuration (Figure 23).

Figure 23. The deformation gradient acting on a line element.

Figure 23. The deformation gradient acting on a line element.

Thus, the deformation gradient can be expressed by

$$F=\frac{\partial x}{\partial X}=grad x$$

(29)

The stretch ratio, λ, is defined as the ratio of the length of a deformed line element to the length of the corresponding undeformed element:

$$\text{λ}=\frac{\left|dx\right|}{\left|dX\right|}$$

(30)

If the line element is extended, upstretched or compressed, then $\text{λ}>1,\text{ λ}=1 \text{or} \text{λ}<1$, respectively.

When the material deforms in a pure stretch, it changes length with no change in the angles between them. If the stretches in these directions are described as ${\text{λ}}_1$, ${\text{λ}}_2$ and ${\text{λ}}_3$, then

$$x_1={\text{λ}}_1{X}_1,x_2={\text{λ}}_2{X}_2,x_3={\text{λ}}_3{X}_3,$$

(31)

and the deformation gradient has only diagonal elements in its matrix:

$$F=\left[\begin{array}{ccc}{\text{λ}}_1& 0& 0\\ 0& {\text{λ}}_2& 0\\ 0& 0& {\text{λ}}_3\end{array}\right],F_{ij}={\text{λ}}_i{\delta }_{ij}$$

(32)

For any deformation, there are always three mutually orthogonal directions along which material undergoes a pure stretch. These directions are called the principal axes of the material and the associated stretches are called the principal stretches. Figure 24 represents the deformation gradient and the related pure stretch of a cross section of the PVC pipe in the two principal directions, respectively.

Figure 24. Maximum principal stretches (undeformed, on the left) and strains (deformed, on the right): (a) vertical direction and (b) horizontal direction.

Figure 24. Maximum principal stretches (undeformed, on the left) and strains (deformed, on the right): (a) vertical direction and (b) horizontal direction.

4.4. Stress/Strain Response

The potential range of the stress/strain response of a material is bounded by two extremes. At one extreme, the response can be perfectly elastic, in conformity to Hook’s law, whereby the magnitude of strain is always proportional to the magnitude of the applied stress. Elastic deformation is instantaneous, which means that total deformation (or strain) occurs at the instant the stress is applied. Upon the release of the external stress, the deformation is instantaneously and totally recovered (Figure 25a). The stress/strain relationship is independent of duration of load application. At the other extreme, deformation caused by the application of a stress is neither instantaneous nor proportional to the stress (i.e., viscous behavior). Deformation is delayed and the rate and the final extent of deformation are dependent on the magnitude and the duration of the applied stress. In addition, the deformation that occurs is not reversible after the stress is released (Figure 25b).

Figure 25. Different strain responses to a constant load: (a) Hook’s law (instantaneous elastic deformation); (b) viscous behavior (deformation is delayed); and (c) real elastic recovery (time-dependent strain).

Figure 25. Different strain responses to a constant load: (a) Hook’s law (instantaneous elastic deformation); (b) viscous behavior (deformation is delayed); and (c) real elastic recovery (time-dependent strain).

The imposition of a stress results in a small instantaneous elastic strain that is then followed by a time-dependent strain. Upon removal of the stress, there is a small elastic recovery of strain that is then followed by a time-dependent recovery. This time dependent recovery occurs more quickly for lower values of initial strain and more slowly for an initially larger strain. While the strain recovery may be total, there is almost always some remaining permanent deformation, which, again, is larger for an initially larger deformation (Figure 25c). Thus, due to its viscoelastic nature, PVC does not exhibit a true elastic behavior (Figure 26).

In Figure 26, the pressure fluctuations and the expansion of the pipe take a nearly synchronous course. During the steady state, for a constant pressure of 4 m, the pipe is no longer in absolute rest, presenting a strain of $5\times {10}^{-5}$. Under the transient regime, pressure oscillations occur within the pipe, causing the extension in the pipe axis and consequently the displacement of the system.

Figure 26. Pressure and first principal strain variation over time.

Figure 26. Pressure and first principal strain variation over time.

5. Discussion

Comparing the two numerical models (MOC and CFD), the pressure variation obtained shows a good agreement with the experimental tests (Figure 27). In both models, the pressure fluctuations are rapidly attenuated due to the mechanical proprieties of the pipe system.

Figure 27. Pressure variation between experimental tests and simulations: (a) for PT2 and (b) for PT1.

Figure 27. Pressure variation between experimental tests and simulations: (a) for PT2 and (b) for PT1.

Considering the Joukowsky overpressure, the results shown in Figure 27 are lower, due to the energy that is transferred from the fluid pulsation to the pipe vibration [23]. Besides the agreement (i.e., low relative error) between the two models and the experimental tests, the modified MOC with damping coefficients prediction assumes the whole system is fixed, otherwise the 3D model with FSI prediction (Table 4—CFD model), where the structure behavior is evaluated in terms of displacements, stress, strain, expansion and compression and internal loads. The size of the pipe system, and the devices installed along the facility (i.e., fittings, bends, elbows, valves) represent a significant component of the head losses. There is a higher energy dissipation due to friction and pipe rheological effects along the experimental pipe system. This energy loss is represented through the decay parameters used in the 1D model, creating an overpressure below the Joukowsky overpressure.

Another difference between the two models (e.g., between 0.5 s and 1.0 s in Figure 27a) relies on the solution of the numerical problem and the spatial discretization. The discretization in the 1D model is based on Courant–Friedrichs–Lewy (CFL) condition of the finite difference method of the MOC [34]. For the FEM, approximation functions are used, based on the division of the function domain into elements such as triangles and tetrahedral elements and on finding nodes within them where the numerical solution is determined. These nodal points are then multiplied by basic interpolation functions.

Table 4. Comparison between experimental data and computational models.

Parameter	Experimental Data	MOC Model	CFD Model
$P_{int}^{max}\left(PT1\right)$ (m)	18.2	16.0	17.9
$P_{int}^{max}\left(PT2\right)$ (m)	17.5	17.0	19.0
${\delta }_t^{max}$ (m)	0.063	--	0.087
${\delta }_{ay}^{max}$ (m)	0.037	--	0.051
${\delta }_{ax}^{max}$ (m)	0.053	--	0.071
Relative error (of pressure)		12%	8%

Table 4. Comparison between experimental data and computational models.

Parameter	Experimental Data	MOC Model	CFD Model
$P_{int}^{max}\left(PT1\right)$ (m)	18.2	16.0	17.9
$P_{int}^{max}\left(PT2\right)$ (m)	17.5	17.0	19.0
${\delta }_t^{max}$ (m)	0.063	--	0.087
${\delta }_{ay}^{max}$ (m)	0.037	--	0.051
${\delta }_{ax}^{max}$ (m)	0.053	--	0.071
Relative error (of pressure)		12%	8%

The maximum amplitude presents a relative error between 8% (CFD model) and 12% (MOC model). In Table 4, the effect of pipe motion, by axial stress waves, is visible in the recorded pressures. The pressure variations can be generated by the lateral vibration of the straight liquid-filled pipe. These can be explained by the axial-lateral coupling mechanism that is evident in the initially curved pipes, in combination with Poisson and Bourdon coupling [35,36].

Both models describe a good agreement with the pressure variation registered in the transducers. However, only in the CFD model is it possible to observe not only the fluid flow simulation over time, but also the displacements achieved by the structure under unsteady conditions. With the use of the finite element method to solve the three-dimensional equations for both fluids and solids it is possible to predict in detail the propagating pressure waves and pipe-response effects [37,38].

In the developed analyses, it was found that the non-conventional dynamic effects present a major role, in particular, in the design of pipe systems under transient pressures. Pipe systems, especially installed above ground, during the occurrence of transients (water-hammer) are under relevant dynamic forces (Figure 28). When these forces are associated to the system movement, a higher FSI is generated, which means that the liquid and the pipe must be analyzed together. The transit pressures and the dynamic forces generated by the water-hammer brings to the system new loads such as internal and external pressures, the latter being ones created from the soil and/or through the pipe supports.

Figure 28. Relevance of the system type and dynamic forces.

Figure 28. Relevance of the system type and dynamic forces.

However, in the case of buried pipes (Figure 28), analyses including water-hammer can be performed using conventional methods, such as the MOC model with special damping coefficients, to account the non-elastic behavior of the fluid, the pipe viscoelasticity and resistance of supports and in the recuperation phase of a deformation. The buried pipelines are generally not very sensitive to the phenomenon of FSI, which when occurs a pump shuts down or by a closing of a valve (internal cause). However, in some circumstances, the supports may be able to generate internal waves to the structure (external cause). In such cases, the design standard procedures of a hydraulic circuit under pressure should include an assessment of the increase loads in the pipes due to pressure waves induced by the FSI effect.

6. Conclusions

Hydraulic transient phenomena associated with a sudden valve closure were investigated by simulation and results analysis using different numerical methods. The current research presents experimental and numerical analyses of fluid structure interaction during a water-hammer event. Transient pressure data collected in a laboratory facility were used to calibrate and validate the transient models (i.e., MOC and CFD).

Concerning the effect of FSI, it is important to highlight the type of supports in the induced movements in pipes during the occurrence of transient phenomena. This analysis is only achieved by applying models that consider various support systems (flexible or rigid), since it has repercussions on transient pressures. Thus, the importance of an integrated analysis using CFD models that can predict and simulate the FSI, accounting for the supports, the geometry and the rigidity of pipes, which influence the resulted unbalanced pressures and dynamic forces in the whole system, is clear. Beyond this, the attenuation or damping effect of the elastic wave appears as a major phenomenon. This damping depends on the type of material and the support properties.

To summarize, when dynamic forces cause the displacement of pipes and fittings, there is an important FSI effect, which implies the need to consider the liquid and solid phases as a whole. Thus, the FSI phenomenon assumes a great importance when the time scale of the structural behavior is less than the fluid behavior and less than the excitation, i.e., the propagation of the pressure waves are quite different from the system components, requiring integrated analyses that include the different behaviors of the different constituent elements of the system.

This work highlights the importance of integrated analysis that depends on the type of system. The specific hydrodynamic 3D model of high viscosity to replicate the operation of a ball valve resulted in an interesting approximation to the experimental tests. Both 3D CFD simulation, using integrated mesh and coupling operator, and 1D model, with specific damping coefficients, can be seen as possible models for solving practical engineering applications, depending on the type of loads and pipe constraints.