Abstract
We consider a dam-water system modeled as a fluid-structure interaction, specifically, a coupled hyperbolic second-order problem, formulated in terms of the displacement of the structure and the fluid pressure. Firstly, we investigate the well posedness of the corresponding variational formulation using Galerkin approximations, energy estimates, and mollification. Then, we apply the finite element method along with the state-space representation of the discrete problem in order to perform a 3D numerical simulation of Cabril arch dam (Zêzere river, Portugal). The numerical model is validated by comparison with available experimental data from a monitoring vibration system installed in Cabril dam.
1. Introduction
The aim of this paper is the analysis and numerical simulation of the dynamics of a dam-water system, in particular, its response to earthquake actions.
Since the core engineering problem lies on the safety and performance evaluation of the structure alone, many related problems have been mostly addressed from a solid mechanics perspective. A simple approach to incorporate the water effects in the dynamic behavior of the dam is to consider an extended Lagrangian model consisting merely of displacements variables in both the reservoir and the dam as in [1, 2] and references therein. However, to take into account the sensitivity of the dam body and foundations to water pressure, the techniques from solid mechanics alone may not be the most adequate. Fluid mechanics plays an important role as well, and appropriate equations describing the water’s behavior and its interaction with the structure should to be taken into account for a more realistic modeling of the system. For this purpose, a pressure-elastic displacement formulation is presented in [3], which is obtained by simplifying the equations of motion of the fluid, so that the velocity variable is eliminated. In [4–6], the authors perform a dynamic analysis of the two-dimensional problem. An advantage of the formulation based on the fluid pressure, a scalar unknown, is the reduction of the computational effort in the numerical simulations. Furthermore, the use of the pressure variable is advantageous due to the possibility of comparison of values from the numerical simulations of the mathematical model with real measurements of the water pressure alongside the dam-water wall. It is common that, in order to check if there exist cracks forming in the concrete and for general safety control [7–10], such structures are kept under constant monitoring, especially in seismic hazardous regions.
To the best of our knowledge, the model we use in this paper for the simulation of dam-water systems, and used in similar engineering studies by other authors, e.g., [11–13], has not been analysed from the mathematical point of view. Therefore, our first aim is the well posedness of the fluid-structure interaction problem. The transmission condition on the fluid-structure interface causes new difficulties because the unbalanced order of the time derivatives of the pressure and the structure’s displacement do not allow a direct application of available mathematical results for linear hyperbolic systems. As a consequence, the regularity we will obtain for the two unknowns is not the same. If we had considered a fluid model formulated in terms of a velocity potential, then the problem would be a hyperbolic coupling similar to the problem treated in [14–17]. In these references, the fluid domain was unbounded and, in order to find the fluid velocity and the displacement field in an elastic body as a result of an incident acoustic wave, the authors used the Laplace transform with respect to the time variable, which is suitable for the application of discretization schemes based on the boundary element method and convolution quadrature method (see also [18]).
In a second stage of our study, the discretization in space by the finite element method of the fluid-structure interaction model is followed by a state space representation and a diagonalization process that lead to a system of first-order linear differential equations. This approach is preferred over numerical time stepping schemes with less computational cost, as suggested in [3], because a modal analysis allows for a direct comparison with natural frequencies experimental data. In practice, such data may correlate eventual deterioration processes with changes in modal parameters over time [8]. Moreover, in the context of seismic responses, the modal analysis plays an important role, as it is crucial that a structure’s natural frequency does not match a frequency of an eventual earthquake; otherwise it may continue to resonate and experience severe structural damage.
The plan of the paper is the following. The mathematical model is presented in Section 2. In Section 3, we derive the variational formulation of the continuous problem in appropriate function spaces and investigate the a priori mathematical properties of solutions; for this, a pressure potential is introduced so that the classical techniques based on Galerkin approximations for second-order hyperbolic equations can be applied to the modified coupled system. The main results on the well posedness of the formulation in pressure-elastic displacement are stated in Theorem 3. We proceed, in Section 4, with the numerical solution of the three-dimensional problem: the discretization in space by the finite element method leads to a state-space representation of the problem (50), (47), and (48), followed by a suitable integration over time (58)-(59). The use of the modal superposition technique allows us to study, in Section 5, the main natural frequencies of the system and corresponding modal configurations. The comparison of these results with available experimental data from a monitoring vibration system of Cabril arch dam (Zêzere river, Portugal) will validate the numerical model. Finally, the seismic response of Cabril dam will be presented by means of displacement fields.
2. Analytical Modeling of the System
The region occupied by the structure (dam body and foundation) is denoted by , while the fluid domain (reservoir) is denoted by . Under the assumptions (i) small displacements and deformations in the structure, (ii) fluid displacement remains small; while interaction is substantial, the domains can be considered fixed, constant in time, and the Eulerian description can be adopted for both the fluid and the structure motions. For the boundaries of and , we have the following: is the interface between and ; represents the ground or rock mass under the fluid; is the air-fluid interface; is an artificial ’wall’, a delimitation of the fluid domain’s extension; represents the rock mass under the structure; and finally, is constituted by the structure walls with contact with the air. The time interval domain during which the dynamic fluid-structure interaction occurs is , .
The formulation of the fluid-structure interaction problem in elastic displacement-pressure variables is the following (see [3], pgs.634-637, and [11]): given the seismic activity and initial conditions , , , and , find and such thatThe gravitational acceleration and the acceleration of the ground due to earthquake vibrations are included in the Navier equation, where is the structure density and is a positive constant representing a damping effect. The Cauchy stress tensor is given by , where and the parameters and are the Lamé constants. They relate to other elastic moduli constants, namely the Young’s modulus and the Poisson ratio , by: , The structure domain boundary is and is the outward unit normal vector at . The viscous effects in the fluid have been neglected and it was assumed that the fluid density varies very little, so that it may be considered constant, . The constant is the speed of sound in the fluid and is the acceleration due to gravity. The fluid domain boundary is . In the air-water interface , a wavelike motion is observed, and is permeable in the sense that the solution of (1b) in should be composed only of outgoing waves on , as no input from this boundary portion exists [3].
3. Variational Formulation. Existence and Uniqueness of Weak Solutions
3.1. Notation and Auxiliary Results
Let be a bounded domain and . We will use the classical notations and results for Lebesgue spaces , , , Sobolev spaces , , and their dual spaces . By and we denote the norms of and , respectively. We will use the notation for vector valued functions, and when and are tensor-valued functions, we will write , where represents contraction of tensors. Analogous meaning holds for .
Suppose that has nonzero -dimensional measure and let Recall that the Korn inequality is valid in : . Therefore is a Hilbert space with respect to the inner product and the associated norm is equivalent to the norm induced by .
If is a Banach space based on the domain , we will write for its dual space and for the evaluation of at . The Bochner spaces associated with will be denoted by , , will denote the space of weakly continuous functions with values in and will be a space of distributions.
We recall the method of mollifiers, which provides approximation by smooth functions. Here we will consider and mollification in the variable. Let the function satisfy the following properties: The mollifier , , is defined by (). For a Banach space , a smoothing operator can be defined by convolution with in the following way: given and , set for all , and let be defined by . Then Since convolution commutes with differentiation, we have and if , where and is a Banach space. Another property of mollification is + .
We will also use the following classical embedding result (see [19], pg. 392, Lemma 11.9).
Lemma 1. Let be Hilbert spaces such that . Then
3.2. Variational Formulation
In order to derive the variational formulation of Problems (1a)-(1j), we consider the function spaces where is equipped with inner product and is equipped with inner product .
Suppose and are smooth solutions of Problem (1a)-(1j). Multiplying both sides of (1b) by and integrating by parts over yields where by (1d), (1e), (1f), and (1c)Therefore Analogously, dot multiplying both sides of (1a) by , integrating by parts over , and using yields Imposing the boundary conditions (1g) and (1h), we obtain
Hence the following variational problem is associated with (1a)-(1j): find and such thatand , , , and , in some sense, to be specified in terms of continuity properties of the solution. In (11a)-(11b) we have used the following notation (see [3]): In (12), and are the so-called mass operators, and are drag operators, and are stiffness operators, and and are interaction operators for the structure and the fluid, respectively. Note that , but the interaction operators appear in (11a)-(11b) in an unbalanced way in terms of the order of the time derivatives, which makes it difficult to obtain good a priori estimates to use in the mathematical analysis of the model.
3.3. Some Considerations on the Regularity of the Pressure: Formal Energy Equation and a Priori Estimates
Let us assume that a solution of (1a)-(1j) exists and is sufficiently regular. If we try to obtain a (formal) energy estimate for such a solution by replacing in (11a) by and by in (11b), followed by time integration in , we end up with the identity
In order to “skew-symmetrize” the bilinear form associated with the equations on , so that the unbalanced terms are not present in the energy equation, it is convenient to introduce a pressure potential for (with respect to the variable): , .
Note that taking into account the relation , in , between the fluid velocity and pressure in this model (see [3], pg. 634), the variable can be seen as a velocity potential for the eliminated unknown in the fluid equations. This is consistent with the assumption of irrotational flow. Now, replacing in (11a) yields while the same replacement in (11b) produces a higher order ordinary differential equation Then where is a constant (distribution) and, more precisely, . If , , and are regular functions with , then the chosen value for is just the value of at . Hence, we consider the new equation for the unknown which is the result of formally integrating (11b) in . The same procedure that leads to (13) allows us to derive a formal energy equation for and (18), in turn, yields the following basic a priori estimates for :where the constants , , depend on the data. Relation (18) suggests that the pair satisfies , and , and . The above considerations also indicate that the time regularity we can expect of the pressure is less than that of the structure displacement, more precisely, with . Moreover, if and satisfies and in , it follows that and . Note that (21) is obtained from (17) taking test functions .
Remark 2. The variational formulations (14) and (17) are actually a formulation in displacement-velocity potential. From and the relation , in , between the fluid velocity and the pressure in this model (see [3], pg. 634), we get .
3.4. Definition of Weak Solution and Well Posedness of the Mathematical Model
Based on the considerations made in the previous section, we shall say that a weak solution to Problems (1a)-(1j) is a pair satisfying (11a)-(11b) andThe main result concerning the well posedness of the mathematical models (1a)-(1j) is as follows.
Theorem 3. Let be locally Lipschitz bounded domains with the common boundary , as described in Section 2. Assume that , , , and . Then there exists a unique solution to Problems (11a)-(11b).
Theorem 3 will be proved along several steps in the next subsections. The a priori estimates (19a)-(19c) for will be used to construct a weak solution for the problem, starting from a generalization of the classical Galerkin method for second-order hyperbolic problems (see [19, 20]). The method used for studying the transmission problem in [16] might also be adapted to the equations satisfied by .
3.5. Galerkin Approximations
Let and be basis of smooth functions for and , respectively. Let and . We seek approximating solutions in the form and with the coefficients and () obtained fromfor , and the initial conditions where and are projection operators onto and , respectively. This is a system of second-order ODEs, which has a unique solution. Indeed, (23a)-(23b) can be written in the form
where we have used the notation and The matrices and are symmetric, positive definite and since the matrix is nonsingular.
Letso that and . Then the approximating functions satisfyNow we multiply (29a) by and (29b) by and sum from to to conclude that and satisfy relations similar to (14) and (17) in the spaces and and, consequently, the a priori estimates (19a)-(19c).
Lemma 4. For the Galerkin approximations it holds that
The uniform bounds collected in Lemma 4 allow us to obtain a pair as a weak limit of the Galerkin approximations such that together with the identitiesfor all and and all null in . We now take and in (32a)-(32b) and obtain From these two relations, it follows and . By Lemma 1 and the well-known density result , we obtain the following.
Lemma 5. and .
3.6. Strong Continuity, Energy Equality and Uniqueness
Let have the weak continuity properties established in the previous section and consider the energy function defined by We shall prove that and use this result to prove the strong continuity of . To accomplish this, we will resort to the following intermediate result.
Lemma 6. If is the weak solution obtained in the previous section then and is an absolutely continuous function.
Proof. Once (35) is proved, since we immediately conclude that is an absolutely continuous function.
To prove (35), we use the mollification method, as described in Section 3.1. Consider , some fixed with , , and let . Inserting the special test functions and , with , in (32a)-(32b) and a simple calculation yieldsNow, since and , we can take and in (37a)-(37b), thus obtaining We let and use the convergence properties (3)-(4) to obtain (35) in . Since are arbitrary, we conclude that (35) holds a.e. in .
The result of Lemma 5 can be improved to
Lemma 7. The pair constructed in the previous subsection satisfies and the energy equality (18).
Proof. Using the weak continuity properties of and the continuity of from Lemma 6, we find that
Uniqueness for the constructed is consequence of the linearity of the problem (14) and (17) and the validity of the energy equality (18).
Concerning uniqueness for the original , it suffices to show that the only weak solution (in the sense of the definition that supports Theorem 3) with null data is . Again, the idea is to introduce a primitive for the original pressure . Since and another primitive of will differ from by a constant, it follows that .
4. Discretization of the Equations
4.1. Semidiscretization in Space
Once a mesh and the corresponding basis functions, and , have been generated, and are approximated by , , and Problems (11a)-(11b) are approximated by (see (26)) Matrices , , and are dimensional, , , and are dimensional, is dimensional, and and are, respectively, - and -dimensional vectors. The matrices and are symmetric, positive definite and the matrix is nonsingular, as already seen in the previous section. It will be convenient to introduce additional notationso that (41) can be written as
4.2. State-Space Representation of the Problem
We proceed with the so-called state-space representation of Problem (41) with subsequent diagonalization, so we end up with a system of independent first-order linear differential equations for the discretized problem. As already mentioned in the Introduction, this is preferred over a direct numerical time stepping scheme applied to (44) because it allows for a comparison with the measured natural frequencies for different values of water level.
Starting from (44) and following [21, 22], a new variable is introduced so that (44) is equivalent to which in turn can be written as The matrix is nonsingular because . Defining and putting , , , we can write (41) as From now on, we consider the system of ODEs with the obvious identifications as concerns . In order to solve (50), the solution of the homogeneous equation can be sought in the form , yielding the generalized eigenvalue problem Let . If all the eigenvalues are distinct, the matrix pencil is diagonalizable; that is, there exist matrices and such that The matrices and contain the left and right eigenvectors of , that is and In practice, since the response of the mechanical system is typically dominated by a relatively small number of the lowest modes, it may not be necessary to compute the complete eigensystem of . Other approaches [12, 13], based on pseudo-symmetric techniques, can be used for computing the mode shapes and natural frequencies of this fluid-structure model.
Using the transformation for modal coordinates and multiplying both sides of (50) by , yields the modal system with and solution of the linear system . From (51) we get or so that the calculation of only requires direct matrix multiplications. Now we have a system of independent first-order linear ordinary differential equations, which is easy to solve. Each equation of the system is of the form where is the -th component of the vector , which in practice is only known in a set of discrete points.
4.3. Time Discretization of the Equations
For the discretization in time of the diagonalized system (52), we take the solution of each equation (53) in its integral form where are the state eigenvalues and , . Each function is known in discrete instants of time with , , and therefore will approximate by a constant in each interval : From the relation and the approximation (55), we get In vector form, the algorithm for computing reads with . Once is known, will provide the solution for the original variable.
5. Numerical Simulations
The objective of this section is to illustrate how the displacement-pressure model can be used for studying the dynamic behavior of a dam. For this purpose, several numerical simulations of Cabril dam (shown in Figure 1(a)) were performed, considering different water levels. This is a common procedure [23–25] as the natural frequencies of the system depend on the reservoir water level. Cabril dam, built in 1954, is the highest Portuguese arch dam. It presents some cracking near the crest which was considered in the computational geometry of the model (these are horizontal cracks, developed since the first filling of the reservoir due to a particular design shape of the dam crest). The numerical results will be compared with available experimental data from a permanent monitoring vibration system installed by LNEC in 2008. More precisely, the computed main natural frequencies (module of the state eigenvalues introduced in Section 4) and the corresponding modal configurations (obtained from the components of the state eigenvectors ) will be compared with experimental data identified from acceleration records. A fast Fourier transform (FFT) technique was used in [26] to obtain the natural frequencies from these acceleration records.
(a) Located at Zêzere river, 132 m high
(b) Finite-element meshes
5.1. Material Parameters and Finite Element Modeling of Cabril Arch Dam
A 3D finite element model of the dam-water system was developed, using meshes composed of serendipity isoparametric finite elements of cubic type, with 20 nodes in each element (see [27]) in the structure and fluid domains.
Concerning the concrete and foundations parameters, Young’s modulus is GPa, the Poisson’s ratio is , and the concrete specific mass is tn/m. Although the damping effect on the structure is well defined through the drag operator , the so-called Rayleigh’s hypothesis was used to model material damping. This assumes as a linear combination of the mass and stiffness matrices, , and we took, as in previous studies [26, 28], and .
The finite element algorithm described in Section 4, based on the geometry shown in Figure 5(b), was implemented in a Matlab program which allows to compute the modal characteristics of the whole system (static analysis) and the corresponding displacements and stress histories under prescribed forces (dynamic analysis), in particular, under seismic loads.
5.2. Comparison between Numerical and Experimental Results of Cabril Dam. Natural Frequencies and Modal Configurations
We recall that the approximation process to simulate the behavior over time of the structure and the fluid ended up in the form the original variable being recovered from a linear combination of the right eigenvectors of the state-space pencil To each vibration mode, i.e., each eigenvector , , there corresponds an eigenmode of the state-space system which is characterized by the way of responding to external forces, i.e., by vibrating with a certain frequency , called a natural frequency of the system, and a certain relative damping coefficient .
Mode shapes and natural frequencies for full reservoir are presented in Figure 2. In the simulations, 202 serendipity finite elements with 20 nodes each were used at the dam body and foundations and 424 elements of the same type were defined at the fluid domain. Here only the displacement of the structure is shown for the first 8 modes. Concerning the graphical representation of the computed mode shapes, it should be noted that each mode corresponds to an oscillatory motion at the frequency which is represented by waves at the nodal points (three waves at each point) whose amplitude and phase are calculated from the complex components of the eigenvector .
In Figure 3 we present a comparison between numerical and experimental values of the main natural frequencies of the dam-reservoir system for different values of the water level. The experimental values of the natural frequencies were obtained during the year 2014 by the analysis of several acceleration records from a permanent vibration monitoring system installed in Cabril dam [26]. The monitoring system includes 16 acceleration sensors (uniaxial sensors for measuring radial accelerations) located in the dam body, near the top: 9 at the crest gallery and 7 in a gallery under the cracked zone. Acceleration files were recorded every hour, during several months, in order to experimentally characterize the influence of the water level on the main natural frequencies of the system. The experimental natural frequencies were identified as abscises of spectral peaks obtained using a FFT technique to analyze the time series of measured accelerations.
The results presented in Figure 3 show a good agreement between identified and computed natural frequencies for different values of water level, meaning that the implemented numerical model is reliable and well calibrated. Note that, for each mode, the natural frequency increases as the water level decreases, as usually calculated with classical approaches based on the simplified model of associated water masses proposed by Westergaard [23] for the simulation of hydrodynamic pressure on the upstream dam face.
5.3. Simulation of Seismic Behavior of Cabril Dam Using the Calibrated Model
We used the accelerogram of the 921 Earthquake, also known as Chichi Earthquake, as input data for ground motions . This accelerogram presents a peak ground acceleration of 0.1g (=1 m/), as shown in Figure 4. Such intensity corresponds to a strong perceived shaking when the earthquake occurs [29] and is in agreement with the characteristics of the seismic zone where the dam is located [30].
(a) Radial displacement history at the top (crest central point)
(b) Displacement field at the instant when the maximum displacement towards upstream occurs (approx. 15 mm)
The corresponding time response was computed considering the diagonalization technique for the time integration described in Section 4, which lead to algorithms (60)-(61). A reduced number of eigenvectors, more precisely 160 vibration modes, were fixed after several numerical tests have been made with a larger number of modes and response values with the same precision were obtained. Figure 5 shows the radial displacement response over time at the top of the central cantilever, with a maximum displacement at the top of the central section of about 15 mm.
The maximum stress values occur at the crest and are of about 1.6 MPa (tension and compression). In conclusion, Cabril dam withstands against the applied earthquake.
6. Conclusions
We studied a fluid-structure interaction model in pressure-displacement formulation which was derived from the Euler and Navier equations, as explained in [3].
The difficulty in the theoretical study of the variational formulation of the coupled equations was caused by the unbalanced order of the time derivatives of the pressure and the structure’s displacement in the interface conditions. We introduced a pressure potential so that the classical techniques based on Galerkin approximations for second-order hyperbolic equations could be applied to investigate the well posedness of the fluid-structure interaction problem.
Concerning the discretization of the problem, the application of the Finite Element Method (FEM) to the space variables was followed by a discretization in time associated with a state space approach appropriate for the dynamic analysis of arch dams. This numerical formulation was tested through a comparative study of the dynamic behavior of Cabril dam, namely, the simulation of dam’s natural frequencies and vibration mode shapes, for different reservoir water levels. The good agreement observed between modal identification outputs (from real measurements) and numerical results shows the reliability of the state space formulation implemented in the developed 3-D FEM program.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
Ana L. Silvestre would like to acknowledge the partial financial support received from CEMAT through Project FCT-UID/Multi/04621/2013.