Fluid–structure interaction (FSI) means that neither a fluid nor a structural system can be solved independently due to the unknown forces in the interface region. The solution for such problems is based on relations of continuum mean mechanics, solved based on finite element methods, being a computational challenge due to its complex geometry, physics of fluids and the interaction between two physic means (i.e., fluid and structure). FSI simulations are mainly divided into simultaneous (direct) and partitioned (iterative) solution procedures. While in the simultaneous procedure a two-way coupled model was used, in the partitioned procedure the one-way coupled model was chosen to work with. This paper analyses an accident that occurred in a main water supply system, in Lisbon, where due to pressure variation, the pipe system suffered successive movements provoking major displacements and leading to rupture of concrete support blocks. According to the diagnosis made in situ after the accident, several analyses were developed addressing different procedures to simulate the accident and to understand the main causes behind this event.
velocity strain tensor
vector of body forces
solid shear modulus
reference configuration of the material deformation
outer normal surface vector of the fluid
outer normal surface vector of the structure
displacement in x axis direction
displacement in y axis direction
fluid stress tensor
Cauchy Stress tensor
velocity of the moving ALE frame
velocity vector of structural displacements
divergence and gradient operator
When dynamic forces cause the displacement of pipes and fittings, there is a significant fluid–structure interaction (FSI), which requires an integrated analysis, considering the liquid and the solid as a whole. Moreover, the interaction between the fluid and structure can be solved using different approaches depending on the characteristics, extension and complexity of the problem. Thus, a finite element method (FEM) of FSI consisting of different types and numbers of elements at the fluid structure (FS) interface can be considered in a partitioned or in a simultaneous solver procedure. A partitioned procedure is more accurate in a large-scale model, whereas a simultaneous approach is required for a smaller scale (Park 1980).
Hence, two coupling methodologies, with different solvers, were used in this study: a partitioned procedure for a one-way coupled and a simultaneous procedure for a two-way coupled. These two procedures were applied in order to assess which best represents the incident that occurred in a real water supply system (WSS). This is a case of interaction that can arise in a typical water pipe system during normal operation. Thus, the phenomenon behind the transference of forces and momentums between the fluid and the pipe wall during the occurrence of flow and structural unsteady conditions is directly related to the safety, reliability and performance of each system type whether WSS, hydraulic circuits of power plants, or industrial pipeline delivery systems (Kratz & Ungar 2003).
This research relates to an accident that occurred in a real WSS, caused by a slow closure manoeuvre of a valve during maintenance operations, which led to a significant movement and the consequent rupture of concrete support blocks. Detailed analyses were also carried out based on the diagnosis made in situ, and simulations performed in the affected area of the system.
The challenge of understanding and foreseeing problems arising from FSI emerges as an interesting topic in engineering and in the science field as a whole. Valve manoeuvres in a pipe system induce pressure variations that propagate and reflect waves in liquid-filled tubes with steep fronts (Koetzier et al. 1986), or the collapse of liquid column separations, leading to dangerous consequences (Bergant et al. 2006). Steep-fronted waves are prone to excite the structural system resulting in pipe movements. On the other hand, displacements generate pressure oscillations inducing FSI (Tijsseling 2007), which are of great relevance when considering hydraulic systems. The timescale of the pipe is lower than the timescale of the fluid and greater than its own excitation, i.e., the propagation velocities and pressure waves are quite different among system components, requiring an integrated analysis that includes the behaviour of different elements of the system (Tijsseling 1996). The interaction manifests itself as vibration behaviour and in disturbances in the velocity and pressure variations of the fluid. The additional pressure loads are transmitted to the mechanical supports and may damage them, inducing dangerous dynamic instabilities (Liu et al. 2012). Moreover, the fluid flow in a pipe network poses a problem of nonlinear nature, due to the fact that the dependence of the flow rate and pressure is nonlinear (Gabrielaitiene et al. 2003).
Hence, the analysis of pressure variations is very important in the design stage of any hydraulic system, as well as in the selection of suitable pipe material and resistances, wall thickness, design of protection devices and operation rules, as well as in the diagnosis and monitoring of anomalies in existing systems (Almeida & Ramos 2010). At an industrial scale, significant levels of noise and vibrations, or transient events, or even phenomena of cavitation can be generated by singularities in pipe systems (e.g., connections, junctions, divergents/convergents, bends and valves) or external loads, inducing instability in the system operation and even damages in the structure components (Ramos 2003; Moussou et al. 2004).
A situation that conveys such events is described by the accident that occurred in the hydroelectric power station of Sayano–Shushenskaya, in Khakassia, Russia in 2009. The report states the accident was primarily caused by turbine vibrations which led to the fatigue and damage of penstock and turbine supports. After the accident, 41 recovered bolts presented fatigue cracks. This violent accident resulted in the flood of the whole turbine engine room and in the collapse of the ceiling of the power station. The entire plant output, totalling 6,400 MW and a significant proportion of the supply to the local electric grid, was lost, leading to widespread power failure in the local area.
This accident was an example of many accidents that can occur in real systems, causing the shutdown of hydraulic systems and compromising support structures. Accidents such as this can be associated with FSIs, which, in recent years, have been the subject of several application domains, since they can arise in many different forms. The complexity of the problem relies on creating and developing new studies of concepts and designs, which require complex geometry, time dependence and multiphysics predictive simulation methods. Most of the problems involving FSI are actively approached in different research fields, such as wind turbines (Bazilevs et al. 2011, 2012), the improvement of parachute systems performance (Kalro & Tezduyar 2000; Takizawa et al. 2011) and even in medicine, namely in the risk assessment of the rupture of abdominal aortic aneurysms (Chandra et al. 2013), as well as in determining the wall tension in cerebral artery aneurysms (Isaksen et al. 2008; Takizawa & Tezduyar 2011) and artery blood flow (Bazilevs et al. 2006; Tezduyar et al. 2010), where identifying and understanding circumstantial events (e.g., strokes and heart attacks) is a challenge.
In the governing equations of the study of FSI, to model the behaviour of a structure, the Lagrangian formulation of motion (i.e., particles are followed in their movement) is required, whereas for the fluid flow the Eulerian formulation is usually applied in order to understand the behaviour of the fluid at a particular spatial position. However, in the presence of a FSI, the fluid domain changes as a function of time, and so an arbitrary Lagrangian Eulerian (ALE) description of motion is needed. This description is a combination between Lagrangian and Eulerian formulations (Donea 1983). The subscripts (s) and (f) are used to represent a quantity that either belongs to the solid or to the fluid, respectively. The fluid is considered incompressible, whereas the solid material is elastic-linear. The domains are described as Ωf and Ωs.
Fluid flow equations
If , the ALE frame (or mesh in a FE discretisation) is not moving, thus the Euler formulation is used. If , the ALE frame is moving with the particles, so the Lagrangian formulation of motion is used.
Coupled fluid and structure equations
In a fluid–structure interface (ΓI), equilibrium and compatibility equations are conditions that must be satisfied. The equilibrium condition, , ensures that the forces are in equilibrium at the interface, where ns and nf are the unit outer normal surface vectors of the fluid and structural domains, respectively. As for the compatibility condition at the interface, us = uf, it ensures the material will not overlap and no gaps will be formed, us and uf being the vectors of displacements of the structure and fluid, respectively.
Furthermore, to discretise a FSI problem, it is convenient to apply different meshes for each system, such as in fluid flows, a finer mesh is needed, when employing the Navier–Stokes or Euler equations for the fluid domain, unlike the structure mesh. To guarantee that equilibrium is satisfied, the forces imposed by the fluid onto the structure need to be calculated and used to construct the fully coupled approach. First the stress tensor of the fluid at the interface is calculated, and from the traction values, structural nodal forces imposed by the fluid onto the structure are then obtained. Therefore, as velocities and displacements at the interface are related by the time integration scheme, all the interface degrees of freedom (DOF) can be expressed as displacements (Bathe 1998). Using the interpolation functions of the structure, the DOF corresponding to the fluid can be expressed as a function of the structural displacements, allowing a compatibility condition between FS.
FEM can be examined as a powerful tool for the approximate solution of mass and momentum differential equations, allowing the description of different physical processes and being extensively applied in engineering practice to obtain the response of flows and structures behaviour. FEM is generally required for a structure solver and also as a fluid solver by the formulation of anALE. The ALE–FEM is widely used in moving boundary, large deformation or in interface problems (Braess & Wriggers 2000).
The ALE method allows arbitrary motion of grid/mesh points to their frame of reference by taking the convection of mesh points into account (Hirt et al. 1974; Hughes et al. 1981). This is the most well-known method used to represent a FSI problem, having the advantage of low computational costs. However, for large movements and rotations of solid or inhomogeneous movements of the grid/mesh points, fluid elements tend to become ill-shaped, which reflects on the accuracy of the solution (Van Loon et al. 2007). Thus, an advantage can be taken, for large-scale systems with the one-way coupled model, by solving the problem sequentially. The ﬂuid solution is computed on the overlapping domain using the Navier–Stokes solver, with a moving grid. The boundary velocity of the overlapping ﬂuid domain is determined by projecting the ﬂuid velocity, obtained from the underlying ﬂuid domain, to the neighbouring 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 ﬁnd the resultant force that generates the motion of the structure (Bathe 1998).
According to Rugonyi & Bathe (2000), partitioned procedures are very attractive since already developed finite element (FE) codes for the fluid flow and structure can be employed. Although powerful, partitioned procedures are not very effective in the solution of FSI problems in which the structure is very flexible, requiring a strong coupling between different means. On the other hand, the simultaneous procedure requires the coupling and solution of all equations. The effectiveness of this procedure can be enhanced by defining the structural DOF prior to the solution. In the simultaneous solution procedure the ﬂuid and solid equations are established and solved together. In a partitioned procedure, each ﬁeld (ﬂuid and solid) is solved separately and solution variables (i.e., forces, or velocities and displacements at the interface) are transferred iteratively from one ﬁeld to the other until convergence is achieved (for each time or load step). In structural and fluid dynamics each equation of motion can be discretised using FE procedures and formulated in terms of the unknown displacement, u. Since the fluid flow equations are expressed using the ALE formulation of motion, the unknown variables include the displacements of the FE mesh nodes (Bathe 1996).
The ALE method is characterised by a part of the calculation domain that moves and deforms according to the fluid motion. This movement can be represented by the moving mesh. To ensure that the flow has a physical consistency, it is necessary that the calculation algorithm meets the geometric law of conservation, which requires that, regardless of the movement of the mesh, it is able to maintain a uniform flow of the solution after the deformation. The ALE method rewrites the equations of fluid motion, in order to understand the differences in geometry and therefore the mesh. This reformulation is performed due to a remapping of the domain, starting from an initial configuration to a configuration that occurred in a given time.
SYSTEM CHARACTERISTICS AND SIMULATION MODELS
Description of the accident
An accident that occurred in a pumping station of a water distribution system in Lisbon (Portugal) served as a case study for this research. The pumping station layout (Figure 1) comprises three sets of parallel pumps which are supplied by a large reservoir of water through two steel pipes: the ramification pipe on the left E, CRE; the ramification pipe on the right D, CPD. These pipes are connected to the horizontal main steel suction pipe placed on concrete supports (labelled CPC). This pipe is connected to the pumps and has two isolating butterfly valves (V58 and V59) and a third valve is placed at the upper pipe (V53) of the pipe branch E.
Owing to routine maintenance and valve replacement, pumps for Sintra were shut down and valve V59 began to be closed by a slow manoeuvre in order to isolate the pipe branch E. During the final instants of the manoeuvre (in V59) the CPC pipe moved, valve V59 moved 0.065 m towards valve V58, joint JP2 opened almost completely (0.07 m) and the pipe branch E moved 0.015 m in the opposite direction.
According to Cesteiro & Ramos (2009, 2010) and Almeida & Ramos (2010), several displacements were measured in the pipe system including pipe CRE and pipe CRD, as well as the pipes connected to the pumps (Figure 1). After this accident, valve V59 was slowly opened and the system remained without any movement. Figure 1 shows the support block rupture identified in situ and the different displacements measured in the pumping station, respectively.
In the one-way coupled method, the problem is partitioned into fluid and structural parts. Data exchanges on the interface (boundary between fluid and solid) are shared between the fluid and the structure solvers. Here, only the fluid pressure acting on the structure is transferred to the structure solver. Then, the received fluid forces are applied as boundary conditions for the structural sub-problem and the structural displacements are obtained. The data exchange on the interface is done once per time-step. The new structural position is considered on the next time-step. Then, the fluid domain is modified to fit its new boundaries and the corresponding flow field is found.
In a two-way coupled algorithm, both parts of the FSI problem are solved simultaneously. For this purpose, one system of equations is created after discretising the governing fluid and structural equations and taking into account the boundary conditions on the interface. Hence, the whole FSI problem is solved at once using a simultaneous scheme.
The geometry system between the V59 and the upper zone (Figure 1) is defined in a 3D software, in order to minimise the number of mesh elements associated with the entire system. Here, two domains are created describing the interior zone (i.e., fluid) and the external volume defined by the pipe thickness (i.e., solid). Both models were simulated in a PC (Intel 5, CPU 3.90 GHz, RAM 8 GB) with four cores and threads running in parallel.
Different geometries for the fluid flow and the structural mechanics are used, where two separate models are set up: one is solved for the turbulent flow within the pipe, using a turbulent flow, k–ɛ interface. The fluid domain is formed by 15,673 tetrahedral elements and 4,734 triangular elements (Figure 2), with 61,221 DOF, corresponding to an average element quality of 75%. The fluid inside the pipe is water, with its own characteristics. As for the structure domain, the total amount of tetrahedral and triangular elements is 57,086 and 37,502 (283,764 DOF), respectively, with an average element quality of 70%. The material applied was steel with a Young modulus of 205 GPa and density of 7,850 kg/m3. In the structural model an identity mapping coupling operator was used to couple the total fluid stress on the pipe surfaces (from the flow simulation). A partial differential equation interface was used to restrict the evaluation of the coupling operator to a single instant. The main advantage lies in plotting the stress components of the fluid and applying them to the structural model without re-evaluating the operator coupling. The boundary conditions applied in both domains are described in Figure 2.
In the fluid flow simulation, the goal was to estimate the initial conditions of the system, before the accident occurred. Thus, a pressure inlet was introduced in the P0 faces as well as a pressure outlet, indicated by P1 in Figure 2. In the pressure inlet, atmospheric pressure was the input (P0). As for the outlet boundary, a function was introduced in order to describe a slow manoeuvre of a valve V59. The condition given by the closing valve is displayed in Figure 2.
For the structural mechanics simulation, it was assumed that all pipe walls were fixed, preventing them from moving in the z-axis. Also, fixed constraints were established in the two faces represented by Face 1 and Face 2 (Figure 2). For Face 3, a free displacement in x-axis was also prescribed, and as for Face 4 a spring-damper was introduced in order to simulate the real support conditions in situ.
After the ﬂuid ﬁeld was solved at a given time instance with an assumed interface location, the resulting ﬂuid pressure and stress were then applied to the structure as external forces.
Since the moving ﬂuid particles interact with the solid boundaries inducing the deformation of the structure, which, in turn, affects the ﬂow motion, the problem is therefore fully coupled (Marti et al. 2006). Usually an important ingredient of the fully coupled analysis is the fact that the fluid domain can be meshed with a much finer discretisation than the structure. Completely different meshes are thus usually employed for the fluid and the structure, and they must adjoin each other in a compatible manner. For the fluid and structure domains, 59,284 and 49,558 tetrahedral elements were needed to perform the mesh (Figure 2), respectively, with 61 vertex elements (183 DOF) and 33,512 triangular elements. The total number of DOF is about 427,062 with an average element quality of 76%.
In this model the major differences, apart from how it is structured, lie in the solver simulation. As previous mentioned, the coupled model was solved simultaneously for a time-dependent study. For each instant, pressure and velocity values were reached and mutually affected by the structure displacement.
One-way coupled model
The pressure variation along the CPC pipe length (x-axis) and along the CRE pipe (y-axis) located 25 m from the main pipe (CPC) is described in Figures 3 and 4. Each line represents the pressure variation for each time instant. In Figure 3, a pressure variation over the time is the input of the upstream condition (Face 3, Figure 2), simulating the effect induced by closing the valve V59. Concerning the Michaud formulation, used for slow manoeuvres, the maximum differential of pressure attained in the simulations is about 70 m, for a time of 2.0 s. Notice that the time operation of the V59 took 10 s until the valve was fully closed. However for the first 2 s, the valve was practically closed but with flow still passing.
As soon as the valve starts to close, the pressure at the downstream end of the pipeline rises changing the pressure condition of the system. At time 2.0 s, a pulse of high pressure in the section immediately downstream of the union of the CPC pipe and CRE pipe is observed. The reservoir is at atmospheric pressure, inducing an unbalanced pressure condition along the CRE pipe. The pressure wave is then reflected in the CRD pipe. A relief wave travels towards the CRD pipe, where the pressure pulse differs from the CRE pipe, due to differences in the pipe section between the two pipes. Here, in CRD pipe, there is exactly the same condition described for the CRE pipe, but with lesser magnitude, until the end of the main pipe (CPC) is reached.
Figure 4 shows the pressure variation along the CRE pipe. Up to 2.0 s the pressure is always rising, as soon as the valve is fully closed. Between 3.5 and 4.0 s the pressure oscillates due to the propagation waves that travel along the pipes. According to the Michaud formulation, the differential pressure obtained here is about 34 m, which is similar to the value obtained through this simulation (Figure 4).
In the case of slow manoeuvres (i.e., with closure time longer than elastic time T > 2L/c) depressure waves near the valve are still being generated when the reflected waves downstream of the reservoir reach the valve, attenuating the first wave. The same behaviour is seen in Figure 4, immediately upstream, in the intersection between CRD and CPC pipes. It can be assumed that the surrounding pressure varies linearly along the entire length of the CRE pipe (Figure 4).
Regarding the solid displacements, after the fluid simulation, the structure movements (Figure 5) and stresses (Figure 6) associated with the fluid friction are obtained for the last instant, through the coupling operator.
Two-way coupled model
The maximum pressure achieved in this model simulation, for t = 4 s, is higher than the one achieved in the one-way model. After t = 3 s the pressure is always higher with a pressure difference of 20 m. Observing the values achieved in the CRE pipe (Figure 8), the pressure, at the end of the closing time of the valve, begins to decline along the waveguide bend until it reaches the intersection with the main pipe (CPC). Comparing it with the fluid simulation, this decline occurs due to the two-coupled simulation, since in each time-step the pressure variation is also influenced by the pipe displacements.
In Figure 9, the displacement is presented for the CRE pipe, relative to the two main axes (x and y). As can be seen from Figure 9(b), the displacements are zero in the fixed supports (Figure 2), and as they move forward on the pipe length, they begin to be felt due to the stresses induced by the fluid in the structure. In order to observe the effect of the spring and damping over the time closure, the simulation was set up to 6 s. The spring-damper is able to represent contact force and elastic deformation during the impact process. In addition, this system describes movements in an opposite direction to the resonance frequency oscillations of the structure, so that when unwanted vibrations occur (like the accident), the causing resonance can be reduced and/or absorbed avoiding the complete destruction of the structure. Thus, due to the characteristics of the system it seems adequate to use a spring-damper at the final end of the CPC pipe (Face 4, Figure 2), allowing movements along the x-axis and decreasing the vibration intensity, while still supporting the pipe.
The maximum displacements achieved in the principal y plan directions, for a time-step of 0.1 s, are almost the same in the x-axis for the two sections described in Figure 9 (ux = −0.045 m). On the other hand, in y-axis, the maximum differs in value and in time. In section 5 the maximum is achieved after 0.1 s (∼uy = 8 × 10−4 m), whereas for section 29, it is reached after 2 s (∼uy = 6 × 10−4 m). Looking at Figure 10(a), it is possible to observe the spring and the damping effect after the closing time of the valve V59. In the main pipe CPC (Figure 10), the maximum displacements registered were about −0.06 m in the x-axis (section 56) and −0.011 m in the y-axis direction.
Comparing simulation results to those obtained in situ, it is possible to conclude that in each model the displacements achieved are similar. Comparing the two models’ results, the pressure difference is about 20 m (i.e., in the two-way coupled model). This difference is associated with the procedure used, since when solving a nonlinear-coupled problem using the one-way model, the fluid solution obtained at the beginning (i.e., the iterative procedure began with the fluid flow equations) may be far from the actual coupled solution, and as a consequence, the iteration procedure may diverge, leading to unreliable results. Also, the displacements estimated in both models differ because in the one-way model the excitation of the structure was calculated only in the final solution time of the fluid component, while in the two-way model the structure and the fluid interact simultaneously in each time-step, inducing oscillations movement. Comparing each model, it is possible to verify that for the one-way model the solution time is faster to obtain, with a quality mesh of 72% for 72,759 finite elements. This is due to the fact that in this model, the simulation time is subdivided into two phases, the fluid and structure, where the first solution (i.e., fluid) is obtained after 14 min while the second (i.e., solid) immediately after 65 s. The latter, despite having a larger number of elements over the fluid domain, it is only simulated for the last instant. Regarding the two-way model, since it is a simultaneous solution, the three displacement components at each point of the structure are calculated. Therefore, this solution takes about 2 h, where in each time-step, the displacements, strains and stresses are obtained.
Even though the displacements achieved in each model are almost equal, the most accurate model to describe the real accident is the two-way coupled, where all the domains interact with each other, as time varies. The overpressure caused in valve V59 triggered successive disorderly movements of the particles of the fluid, causing the misalignment of the conduit due to the mass inertia at valve V59.
Figure 11 shows the comparison of vorticity between one-way and two-way models, respectively. The vorticity becomes stronger along the 20 m from valve section V59, at the maximum velocity. For the same instant, the mean maximum velocity magnitude achieved is different in the one-way and two-way models, respectively. This means that starting from t = 3.5 s, in the two-way model, the structure displacements are felt, providing an increase in the velocity magnitude. For the highest velocity, the vorticity is stronger particularly near the pipe walls of CRE and CRD pipes and downstream of valve V59 (Figure 11). In the two-way coupled model (Figure 11(a)), the distribution of velocity vectors presents the highest intensity and magnitude in CRE and CRD pipes and at the entrance of CPC pipe, compared to the one-way coupled model.
Both models change the dynamic response of the structure and induce perturbation in the flow. However, the two-way coupled model shows, in an integrated way, how the instability of the structure is affected by the state of fluid flow.
The coupling fluid–structure phenomenon is seen in this research as a great approach to describe an accident induced by pressure variation and significant displacements which put the Lisbon WSS and neighbouring municipal water services out of operation. The hydrodynamic pressure tends to be amplified due to the fluid motion created by perturbations, which, in turn, lead to additional structure deformations. Consequently, these deformations will also modify the hydrodynamic pressure. Another important remark is related to the increasing turbulence detected inside the pipe, being an important key factor for the instability of the whole set (e.g., structure behaviour). Even though modelling in 3D formulation takes time and requires complex analysis, it provides more interpretative data allowing a better understanding of diagnosis analyses and problem management.
The two model approaches together with diagnosis were used in the field to explore and detect the real causes related to the origin of the problem that occurred in the Sintra water conveyance pipe system. In the two-way coupled model, significantly more iterations were required at each time to converge to the solution, when compared to the one-way model. However, although a greater number of DOF were required in the two-way model in order to accurately solve for the fluid flow, when compared to the structural response, iterations were only applied for the nonlinearities of the problem, thereby reducing the calculation time.
Since in this research the problem is compliant, or complex, and the number of structural DOF is significant when compared to the fluid DOF, it is more suitable to solve the coupled equations through the two-way model. If the structure is very flexible, a greater dependency between the two domains occurs, and the solution is coupled. If the structure is rigid and barely deforms, the solution adopted can be the one-way coupled model. In addition, the number of elements has implications on the type of model (e.g., geometry size, singularities and structure's configuration), since with greater complexity comes greater simulation time. Therefore, it is necessary to verify how to simulate the geometry and what type of models should be used to represent the reality.
Hence, from the two models presented, the one that best describes the accident is the two-way model, due to the system characteristics, such as the geometry, the relation between the DOF in the fluid–structure system and the dynamic behaviour of the fluid–structure in an integrated way (i.e., the influence of the fluid movements in the structure and consequently the displacements of the structure in the fluid oscillations).
This research was funded by the Portuguese Foundation for Science and Technology (FCT) through the Doctoral Grant (SFRH/BD/68293/2010) to M. S. Also thanks to CEHIDRO, Hydro Systems Research Centre from the Department of Civil Engineering, Architecture and Georesources (DECivil), at Instituto Superior Técnico, Lisbon.