A Fluid-Structure Interaction Model for Dam-Water Systems : Analytical Study and Application to Seismic Behavior

We consider a dam-water systemmodeled 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 usingGalerkin approximations, energy estimates, andmollification.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.


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][5][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][8][9][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][12][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 2 Advances in Mathematical Physics 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][15][16][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.

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 Ω 푓 , Ω 푠 ⊂ R 3 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: Γ 1 is the interface between Ω 푠 and Ω 푓 ; Γ 2 represents the ground or rock mass under the fluid; Γ 3 is the air-fluid interface; Γ 4 is an artificial 'wall' , a delimitation of the fluid domain's extension; Γ 5 represents the rock mass under the structure; and finally, Γ 6 is constituted by the structure walls with contact with the air.The time interval domain during which the dynamic fluid-structure interaction occurs is [0, ],  > 0.
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  →  0 ,  → V 0 ,  0 , and  0 , find  →  =  →  (, ) and  = (, ) such that The gravitational acceleration  →  and the acceleration of the Advances in Mathematical Physics 3 The structure domain boundary is 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  fl |  →  | is the acceleration due to gravity.The fluid domain boundary is In the air-water interface Γ 3 , a wavelike motion is observed, and Γ 4 is permeable in the sense that the solution of (1b) in Ω 푓 should be composed only of outgoing waves on Γ 4 , as no input from this boundary portion exists [3].
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  푟 (, ; ), 1 ≤  ≤ ∞,  푤 ([, ]; ) will denote the space of weakly continuous functions with values in  and D 耠 (0, ; ) will be a space of distributions.
Another property of mollification is supp  훿 ⊂ supp J 훿 + supp .

Advances in Mathematical Physics
Therefore Analogously, dot multiplying both sides of (1a) by  →  ∈  푠 , integrating by parts over Ω 푠 , and using (  →  ) : Imposing the boundary conditions (1g) and (1h), we obtain Hence the following variational problem is associated with =  0 , and  푡 (0) =  0 , 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), M and S are the so-called mass operators, C 푠 and C 푓 are drag operators, K and H are stiffness operators, and Q 푠 and Q 푓 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.in (11b), followed by time integration in [0, ], we end up with the identity

Some Considerations on the Regularity of the
In order to "skew-symmetrize" the bilinear form associated with the equations on Γ 1 ×]0, ], so that the unbalanced terms in the energy equation, it is convenient to introduce a pressure potential for  (with respect to the  variable):  ∈ D 耠 (0, ;  푓 ),  =  푡 .
Note that taking into account the relation 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  → V 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, which is the result of formally integrating (11b) in [0, ].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  푖 ,  = 1, 2, 3, depend on the data.Relation (18) suggests that the pair and  푡  →  ∈ ([0, ];  푠 ).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 (0) =  0 .Moreover, if  =  푡  and satisfies Advances in Mathematical Physics and ∇ 푥 (0) = 0 in Ω 푓 , it follows that  푡  ∈ ([0, ];  −1 (Ω 푓 )) and (0) =  0 in  −1 (Ω 푓 ).Note that ( 21) is obtained from (17) taking test functions  →  ∈  1 0 (Ω 푓 ).
The main result concerning the well posedness of the mathematical models (1a)-(1j) is as follows.
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 E(  →  훿 ,  훿 ) 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 (  →  , ) = (  → 0 , 0).Again, the idea is to introduce a primitive  ∈ D 耠 (0, ;  푓 ) for the original pressure  =  푡 .Since  = 0 and another primitive of  will differ from  by a constant, it follows that  = 0.

Discretization of the Equations
Matrices M, C 푠 , and K are (3 푠 ) × (3 푠 ) dimensional, S, C 푓 , and dimensional, and f and s are, respectively,  푠 -and  푓dimensional vectors.The matrices M and S are symmetric, positive definite and the matrix is nonsingular, as already seen in the previous section.It will be convenient to introduce additional notation so 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 z = ẏ is introduced so that (44) is equivalent to which in turn can be written as The matrix and putting k fl u , q fl ṗ , x fl [y z] ⊤ = [u p k q] ⊤ , we can write (41) as From now on, we consider the system of ODEs B ẋ () = Ax () + F () ,  ∈ ]0, ] with the obvious identifications as concerns x 0 .In order to solve (50), the solution of the homogeneous equation B ẋ () = Ax(),  ∈]0, ] can be sought in the form x() = b 휆푡 , yielding the generalized eigenvalue problem Ab − Bb = 0. Let  fl dim(A) = dim(B) = 2(3 푠 +  푓 ).If all the eigenvalues are distinct, the matrix pencil (A, B) is diagonalizable; that is, there exist  ×  matrices W and V such that The matrices W and V contain the left and right eigenvectors of (A, B), that is  푖 = w * 푖 Ak 푖 and  푖 = w * 푖 Bk 푖 .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 (A, B).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 x = Vz and multiplying both sides of (50) by W * , yields the modal system with C = W * F and z 0 solution of the linear system Vz 0 = x 0 .From (51) we get z 0 = Ω −1 퐴 W * Ax 0 or z 0 = Ω −1 퐵 W * Bx 0 so that the calculation of z 0 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 C, which in practice is only known in a set of discrete points.

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  푖 fl  푖 / 푖 are the state eigenvalues and c푖 =  푖 / 푖 ,  = 1 : .Each function  푖 is known in discrete instants of time 0 <  1 < ⋅ ⋅ ⋅ <  푛−1 <  푛 =  with Δ fl  푘+1 −  푘 ,  = 1 :  − 1, and therefore will approximate  푖 by a constant in each interval [ 푘 ,  푘+1 ]: From the relation and the approximation (55), we get In vector form, the algorithm for computing z reads with Λ = diag( 1 , . . .,  푚 ).Once z is known, will provide the solution for the original variable.

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][24][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 u 푖 of the state eigenvectors k 푖 ) 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.

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  = 32.5GPa, the Poisson's ratio is ] = 0.2, and the concrete specific mass is 2.45tn/m 3 .Although the damping effect on the structure is well defined through the drag operator C 푠 , the so-called Rayleigh's hypothesis was used to model material damping.This assumes C 푠 as a linear combination of the mass and stiffness matrices, C 푠 = M + K, and we took, as in previous studies [26,28],  = 0.5 and  = 0.0025.
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.

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 x being recovered from a linear combination of the right eigenvectors of the state-space pencil (A, B) To each vibration mode, i.e., each eigenvector k 푖 ,  ∈ {1, . . ., }, 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  푖 = −Re( 푖 )/ 푖 .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 u 푖 .
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.

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/ s 2 ), 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].
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.

Conclusions
We studied a fluid-structure interaction model in pressuredisplacement 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 secondorder 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.
we obtain the following.