STRUCTURAL ASSESSMENT OF GRAVITY DAMS BY MEANS OF A SEMI-DISCRETE APPROACH

Abstract

A new numerical model for the structural assessment of gravity dams by means of a semi-discrete approach is proposed. Gravity dams are massive structures, which their stability depends on the gravity loads applied into the structure. Mainly, its structural assessment is performed by means of a gravity approach. However, this approach is too conservative and mostly does not reflect the real structural behaviour of the dam. In this context, there is the need of models that are simplified enough to allow a simple and fast parametric analyses. The proposed model idealizes the dam as a set of rigid elements, where the damage and the deformation are concentrated in the contact sides between adjacent elements. Thus, the elements are rigid, but the material is considered as deformable. As the proposed model is semi-discrete, it can detect separation or sliding between elements. However, initial contacts do not change during the analyisis and a relative continuity among elements exists, in order to simplify the computational effort. The effective performance of the proposed model is demonstrated by numerical validation and by comparisons with some numerical models presented in the literature.

Keywords: Gravity dams, Semi-discrete model, non-linear behaviour, rigid body spring model, damage

1. Introduction

Gravity dams are massive structures, which their stability depends on the gravity loads applied into the structure. Mainly, the structural assessment of gravity dams is performed by means of a gravity approach, where the resultant of all forces acting on the dam must lie in the third middle of the base. This analysis must calculate the [1]: a) Position of the resultant force, where the resultant force must lie in the middle third of the base; b) Inclination of the resultant force, in order to evaluate the shear forces and the possible sliding of the dam; and c) Compressive stresses, in order to avoid the crushing of the material. Figure 1 depicts the forces acting on the dam; where P_H is the hydraulic pressure (v=vertical, h=horizontal and e=downstream), P_p is the own weight, P_s is the uplift pressure, P_A is the earth pressure, S is the seismic loads (V=vertical, H=horizontal) and P_O is the hydrodynamic pressure.

Error creating thumbnail: File missing

Figure 1 Forces acting on the dam

However, this approach is too conservative and mostly does not reflect the real structural behaviour of the dam. Although there are several authors whom have proposed different analytical and numerical models [2-7], the study of gravity dams by using simple models is still an open problem.

In this context, there is the need of models that are simplified enough to allow a simple and fast parametric analyses, but they should also take into account the peculiar behavior of the material. Thus, in this paper a new numerical model for the structural assessment of gravity dams by means of a semi–discrete approach is presented.

The proposed model follows the philosophy of the Rigid Body Spring Models (RBSM) [8 - 10]. This philosophy is based on:

• There is a limited number of damage mechanisms. Thus, zones, types and directions of the possible damage can be known a priori.

• Damage is not concentrated in a single point of the structure.

• There are zones with similar behaviour of a rigid body, between them there are damaged zones.

In this study, only the in–plane deformations are considered. The dam is idealized as an assemblage of rigid elements. Damage and deformation are concentrated in the contact sides between adjacent elements. These elements are quadrilateral and have the kinematics of rigid bodies with two linear displacements and one rotation (Fig. 2). Three devices (springs) connect the common side between two rigid elements or the restrained sides. These connections are two axial devices, separated by a distance 2b that take into account a flexural moment, and one shear device at the middle of the side.

Error creating thumbnail: File missing

Figure 2 Forces and displacements

The material is considered deformable but this deformation is concentrated in the connecting devices, while the element is not deformable. Each connecting device is independent of the behavior of the others connecting devices and depends only on the Lagrangian displacements. In other words, the connecting device represents the elastic and post–elastic mechanical characteristics of the material and, at the same time, represents the capacity of the model to take into account the separation or the sliding between elements.

The proposed model was developed as a semi–discrete element model (SDEM). Therefore, it can detect separation and sliding of the elements. However, initial contacts do not change during the analysis and a relative continuity among elements exists, in order to simplify the computational effort. Thus, overlapping, separation or sliding between adjacent elements can occur. Numerically, these mean compression, tension or shear in the connecting devices. The semi–discrete model can be though as an analysis technique that combines the advantages of the discrete analysis techniques (e.g. it considers the relative motion among elements) with the computational advantages of the continuous analysis technique (e.g. no new contacts update is necessary).

2. Mathematical Formulation

The dam is considered as two–dimensional plane solid model Ω, partitioned into m quadrilateral elements ωi such that no vertex of one quadrilateral lies on the edge of another quadrilateral. The global Cartesian coordinate frame {O, x, y} is placed in order to have the gravity acceleration g applied in the negative y–axis direction. A local reference frame {o_i, ξ_i, η_i}, whose axes are initially parallel to the global reference frame, is fixed in each element’s barycenter o_i. These elements are rigid, so the displaced configuration of the discrete model is described by the position of these local reference frames, as shown in Figure 2b. Given the local coordinate of a point (ξ, η), the displacements (Δx, Δy) in the x − y plane are evaluated as follows:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://test.scipedia.com:8081/localhost/v1/":): {\textstyle \left\{ \begin{matrix}\Delta x\\\Delta y\end{matrix}\right\} =} Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://test.scipedia.com:8081/localhost/v1/":): \left[ \begin{matrix}1&0&\eta \\0&1&-\xi \end{matrix}\right] \left\{ \begin{matrix}{u}_{e}\\{v}_{e}\\{\psi }_{e}\end{matrix}\right\}

(1)

The translation components u_i, v_i and the rotation angle ψ_i associated with each element, are collected into the vector of Lagrangian coordinates {u}. The loads are condensed into three resultants associated with each element: the forces p_i and q_i applied to the element centroid considering the initial undeformed geometry, and the couple μ_i. They are assembled into the vector of external loads {f_e} which is conjugated in virtual work with {u} as follows:

{u}^T = {u₁, v₁, ψ₁, u₂, v₂, ψ₂, ..., u_m, v_m, ψ_m} (2)

{f_e}^T = {p₁, q₁, μ₁, p₂, q₂, μ₂, ..., p_m, q_m, μ_m} (3)

The elements are interconnected by connecting devices (line springs) placed along each side, in correspondence of three points named P, Q and R, as shown in Figure 3. Three average strain measures are associated with these connecting devices: the axial strains, ε^P and ε^R are associated with the volumes of pertinence V^P and V^R, while the shear strain ε^Q is associated with the volume V^Q = V^P + V^R. Considering a discrete model with r sides which connect all the elements (interfaces), the vector of generalized strain {ε} and the diagonal matrix of volumes of pertinence [V] (Fig. 4) are defined as follows:

{ε}^T = {ε^P₁, ε^Q₁, ε^R₁, ε^P₂, ε^Q₂, ε^R₂, . . . , ε^P_r , ε^Q_r, ε^R_r} (4)

[V] = Diag{V^P₁ , V^Q₁, V^R₁, V^P₂, V^Q₂, V^R₂, ..., V^P_r, V^Q_r, V^R_r} (5)

Error creating thumbnail: File missing

Figure 3 Assembly of rigid elements

a)

Error creating thumbnail: File missing

b)

Error creating thumbnail: File missing

c)

Error creating thumbnail: File missing

d)

Error creating thumbnail: File missing

Figure 4 Volume of pertinence: a) Tensile volume, b) compressive volume, c) shear volume, d) sliding volume

Under small rotation assumption, the strain-displacement relation can be expressed by considering a 3r x 3m matrix [B] as follows:

{ε} = [B]{u} (6)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://test.scipedia.com:8081/localhost/v1/":): {\textstyle \left[ B\right] =\frac{1}{{h}_{i}+{h}_{j}}\left[ \begin{matrix}\mathrm{cos}\,{\alpha }_{i}&\mathrm{sin}\,{\alpha }_{i}&\left[ \mathrm{sin}\,\left( {\alpha }_{i}-\vartheta \right) {d}_{i}+b\right] \quad \\-\mathrm{sin}\,{\alpha }_{i}&\mathrm{cos}\,{\alpha }_{i}&\left[ \mathrm{cos}\,\left( {\alpha }_{i}-\vartheta \right) {d}_{i}\right] \\\mathrm{cos}\,{\alpha }_{i}&\mathrm{sin}\,{\alpha }_{i}&\left[ \mathrm{sin}\,\left( {\alpha }_{i}-\vartheta \right) {d}_{i}-b\right] \end{matrix}\begin{matrix}-\mathrm{cos}\,{\alpha }_{i}&-\mathrm{sin}\,{\alpha }_{i}&-\left[ \mathrm{sin}\,\left( {\alpha }_{i}-\vartheta \right) {d}_{i}+b\right] \\\mathrm{sin}\,{\alpha }_{i}&-\mathrm{cos}\,{\alpha }_{i}&-\left[ \mathrm{cos}\,\left( {\alpha }_{i}-\vartheta \right) {d}_{i}\right] \\-\mathrm{cos}\,{\alpha }_{i}&-\mathrm{sin}\,{\alpha }_{i}&-\left[ \mathrm{sin}\,\left( {\alpha }_{i}-\vartheta \right) {d}_{i}-b\right] \end{matrix}\right]}

(7)

Where α_i is the angle of the connection side of element i referring to ξ–axis and ϑ is called distortion angle. A measure of stress, work–conjugated to the strain, is assigned to each connecting device, and is assembled into the vector {σ} as follows:

{ σ }^T = { σ ^P₁, σ ^Q₁, σ ^R₁, σ ^P₂, σ ^Q₂, σ ^R₂, . . . , σ^P_r , σ^Q_r, σ ^R_r} (8)

Where σ^P and σ^R are the axial stresses in the connection point P and R, and σ^Q is the shear stress in Q. Forces are related by:

{f_e} = [B]^T{σ} (9)

The constitutive law correlates the strains and stresses:

{σ} = [D]{ε} (10)

Where [D] is the tangential stiffness matrix of the connection side:

[D] =Diag [k^P, k^Q, k^R] (11)

Replacing Equation 6 in Equation 10 and this in Equation 9, it obtains:

{f_e} = [B]^T[D][B]{u}=[K]{u} (12)

3. Mechanical Characteristics of the Interfaces

3.1 Elastic properties

The elastic characteristics of the connecting devices are assigned with the criterion of approximating the strain energy of the corresponding volumes of pertinence in the case of simple deformation. For an orthotropic material in plane deformation, the matrix of elasticity is given by:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://test.scipedia.com:8081/localhost/v1/":): {\textstyle \left[ A\right] =\left[ \begin{matrix}{A}_{11}&{A}_{12}&0\\{A}_{21}&{A}_{22}&0\\0&0&{A}_{33}\end{matrix}\right]}

(13)

Where A₁₁=A₂₂=E/(1-v²), A₁₂=A₂₁=vE/(1-v²), A₃₃=2G=E/(1+v²); E is the Young modulus, v is the Poisson’s coefficient and G is the shear modulus.

On the other hand, the stress Σ and the strain Η vectors are:

{Σ}^T = {Σ₁₁, Σ₂₂, Σ₁₂} (14)

{Η}^T = {Η₁₁, Η₂₂, Η₁₂} (15)

The siffness of the elastic devices is obtained by equating the strain energy of the material Π_m and the strain energy of the connections Π_c:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://test.scipedia.com:8081/localhost/v1/":): {\textstyle {\Pi }_{m}=\frac{1}{2}{\left\{ \epsilon \right\} }^{T}\left[ A\right] \left\{ \epsilon \right\} Vol=} Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://test.scipedia.com:8081/localhost/v1/":): \frac{1}{2}{\left\{ u\right\} }^{T}\left[ k\right] \left\{ u\right\} ={\Pi }_{c}

(16)

Thus, the axial and shear stiffness are:

k^P = k^R =A₁₁ (17)

k^Q = A₃₃ (18)

In addition, the two axial devices are separated from the middle point of the side by a length b in order to take into account the bending momento, where b=l/(2√3).

3.2 Strength and Plastic properties

The monotic constitutive laws are assigned to the connecting devices adopting a phenomenological approach. These laws are based on experimental monotonic tests currently available in literature. Different rules are assumed for the axial devices and for the shear device, as sketched in Figure 5. For the axial spring, the skeleton curve under compression is given by:

a)

Error creating thumbnail: File missing

b)

Error creating thumbnail: File missing

Figure 5 Constitutive laws: a) axial, b) shear

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://test.scipedia.com:8081/localhost/v1/":): {\textstyle \sigma =\, {E}_{0}exp\left( \frac{-\epsilon }{{\epsilon }_{c}}\right)}

(19)

Where E₀ is the initial elastic modulus and ε_c is the strain at the peak compression strength σ_c. Along this skeleton curve, the spring stiffness (k^P, k^R) for compression loading is:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://test.scipedia.com:8081/localhost/v1/":): {\textstyle {k}^{P}={k}^{R}={E}_{0}\left( 1-\frac{-\epsilon }{{\epsilon }_{c}}\right) \mathrm{exp}\,\left( \frac{-\epsilon }{{\epsilon }_{c}}\right)}

(20)

The tensile axial response is defined by a tri-linear skeleton curve identified by the couples of points (σ_t, ε_t) and (σ_r, ε_r) which correspond to the peak and residual strengths. The plastic response of each axial connection is independent from the behaviour of any other connection device.

Symmetric stiffness and strength have been attributed to the shear connections. The skeleton curve is tri-linear, defined by four parameters: the initial shear stiffness G, the softening stiffness G_r, the maximum shear strength τ and the residual shear strength τ_r. The shear strength is related to the stresses of the axial connections according with Mohr-Coulomb criterion:

τ = c - σ tan(ϕ) (21)

where c is the cohesion, σ is the axial stress and ϕ is the internal friction angle.

4. Validation

The model was validated by using the discrete element model of the Guil1hofrei dam (Portugal) performed by Bretas et al. [7]. This is a masonry gravity dam, built in 1938, with a maximum height of 39 m and a total length of 190 m. The soil foundation of the dam is a granitic rock mass, of good quality (Fig. 6). Table 1 shows the mechanical properties of the masonry and the soil foundation. For details of the dam characteristics please refer to [7, 11]. Figure 7 shows the SDEM and the DEM meshes for the studied dam.

Error creating thumbnail: File missing

Error creating thumbnail: File missing

Figure 6 Guilhofrei dam [11]

a)

Error creating thumbnail: File missing

b)

Error creating thumbnail: File missing

Figure 7 Numerical models: a) SDEM, b) DEM [11]

Table 1: Mechanical properties of the materials [7]

Property	Dam	Foundation
Volumetric weight [kN/m3]	24	25
Elasticity modulus [GPa]	10	10
Shear modulus [GPa]	4	4
Poisson coefficient	0.2	0.2
Compressive strength [MPa]	10	Elastic
Tensile strength [MPa]	1	Elastic
Cohesion [MPa]	1.58	Elastic
Friction angle [grad]	55	Elastic

Four different load cases were considered:

• Load Case 1: Own weight (PoPo)

• Load Case 2: Hydrostatic pressure (PH)

• Load Case 3: Own weight plus hydrostatic pressure (PoPo + PH)

• Load Case 4: Own weight plus hydrostatic pressure plus uplift pressure (PoPo + PH + Ps)

4.1 Load Case 1: Own weight

For the own weight analysis, only the volumetric weight of the dam was taken into account. Table 2 shows the results obtained by DEM [7] and the proposed model. It can be observed that the results are practically the same for both models. The own weight and the compression stress have error percentages are less than 10%. Although the error for displacements are around 15%, the displacements are in mm. Thus, small variations give great errors.

Table 2: Results of the own weight analysis

Result	DEM	SDEM	Error [%]
Own weight [kN]	9,700	9,500	2.0
Maximum compressive stress [MPa]	0.84	0.87	7.4
Maximum horizontal displacement [mm]	2.5	2.1	14.5

a)

Error creating thumbnail: File missing

b)

Error creating thumbnail: File missing

Figure 8 Own weight deformed: a) SDEM, b) DEM [11]

a)

Error creating thumbnail: File missing

b)

Error creating thumbnail: File missing

Figure 9 Vertical axial stress maps due to own weight: a) SDEM, b) DEM [11]

Figure 8 shows the deformation produced by the dam’s own weight, in which it can be observed that it deforms slightly upstream. This coincides with the real phenomenon, since most of the mass is on this side of the dam. So that when the hydrostatic pressure will apply, the forces remain in equilibrium. Figure 9 shows the stress maps for vertical axial stresses. It canbe seen that the maximum compressive stress is located upstream at the foot of the dam, for both models.

4.2 Load Case 2: Hydrostatic Pressure

In this analysis, only the hydrostatic pressure was considered without taking into account the dam’s own weight. The stability of the dam is due to the tensile strength of the material. The resulting hydrostatic pressure was approximately to 5,000 kN/m. Table 3 shows the results obtained, it can be observed that the results are acceptable with a minimum error. Figure 10 shows the deformation produced by hydrostatic pressure. It can be seen that the dam rotates slightly downstream. There is a slightly tensile damage in the base of the dam upstream (Fig. 11), since the own weight of the dam is not considered.

Table 3: Results of the hydrostatic pressure analysis

Result	DEM	SDEM	Error [%]
Maximum tensile stress [MPa]	0.84	0.89	6.0
Maximum compressive stress [MPa]	0.77	0.76	1.0
Maximum horizontal displacement [mm]	6.5	5.7	12.3

a)

Error creating thumbnail: File missing

b)

Error creating thumbnail: File missing

Figure 10 Deformed mesh due to hydrostatic pressure: a) SDEM, b) DEM [11]

a)

Error creating thumbnail: File missing

b)

Error creating thumbnail: File missing

Figure 11 Vertical axial stress maps due to hydrostatic pressure: a) SDEM, b) DEM [11]

4.3 Load Case 3: Own Weight plus Hydrostatic Pressure

In this analysis, the own weight of the dam is first applied and then the hydrostatic pressure load, since the analysis is non–linear. Table 4 shows the obtained results. It can see that no tensile stresses are in the dam due to the own weight load (Fig. 12). In this context, the maximum horizontal displacement is lesser than the load case 2 (when the own weight is not considered, see Fig. 13). This means that the own weight contributes to the stability of the dam.

Table 4: Results of the load combination 3

Result	DEM	SDEM	Error [%]
Maximum tensile stress [MPa]	0.29	0.33	12.1
Maximum compressive stress [MPa]	0.96	1.10	14.5
Maximum horizontal displacement [mm]	4.0	3.4	15.0

a)

Error creating thumbnail: File missing

b)

Error creating thumbnail: File missing

Figure 12 Vertical axial stress maps due to load case 3: a) SDEM, b) DEM [11]

a)

Error creating thumbnail: File missing

b)

Error creating thumbnail: File missing

Figure 13 Deformed mesh due to load case 3: a) SDEM, b) DEM [11]

The error percentages are around 15%. This can be explained in terms of the fundamental assumptions implied. The DEM model used by Bretas et al. [7] is a particular type of DEM with deformable blocks. In this case, the dam block is discretized into a mesh of 4-node elastic finite elements. Only the dam-rock joint is nonlinear. The errors reported are thus expectable, and can be considered acceptable. The issue is that the REM model has computational advantages over this deformable block DEM, since its performance is not as good.

4.4 Load Case 4: Own Weight plus Hydrostatic Pressure plus Uplift Pressure

This load case considers additionally to the own weight and the hydrostatic pressure, the uplift pressure and a flood of 5 m over the crown of the dam (failure load). The resultant of the uplift pressure load was equal to 1,015 kN. This type of combination loads are similar to the failure loads of the dam. As the previous case, the own weight of the dam is first applied and then the hydrostatic pressure and uplift pressure loads are applied, since the analysis is non–linear.

The maximum compressive stress at the base of the dam downstream was equal to 1.81 MPa for the DEM [7] and 2.1 for the SDEM (16% of error). Figure 14 shows the deformed mesh and the failure mechanism of the dam. The dam overturns downstream, since the tensile stresses at the base of the dam are overpassed.

a)

Error creating thumbnail: File missing

b)

Error creating thumbnail: File missing

c)

Error creating thumbnail: File missing

d)

Error creating thumbnail: File missing

Figure 14 Deformed mesh a) SDEM, b) DEM [11]; and failure mechanism c) SDEM, d) DEM due to load case 4.

5. Final Remarks

In this paper, a new model for the structural assessment of gravity dams by means of a semi–discrete approach is presented. This model can detect sliding, separation, overturning, crushing, tensile and shear damage. Thus, the proposed model can detect the different collapse mechanism of the dams, mainly: overturning and sliding.

It was validated by comparing with a discrete element model of a dam. The validation of the model was taking into account different load cases. One advantage of the proposed model is that it is no necessary the mechanical properties of the interfaces required in a discrete element model. The mechanical properties of the material are concentrated in the connecting devices between adjacent elements.

The tensile strength of the material is an important parameter in the structural assessment of gravity dams. This parameter is usually neglected, as it cannot be relied upon, given the uncertainties about the contact between the two materials. Therefore, it is required that the dam be stable even without tensile strength along the base joint.

References

[1] W. Creager, Engineering for masonry dams, John Wiley & Sons, Inc., New York, 1917.

[2] S. Bennati, M. Lucchesi, The minimal section of a triangular masonry dam, Meccanica 23 (1988) 221–225.

[3] A. Carpinteri, S. Valente, G. Ferrara, L. Imperato, Fracture mechanics of concrete structures, Taylor & Francis, London, 1992, Ch. Experimental and numerical fracture modeling of a gravity dam, pp. 351–360.

[4] S. Bhattacharjee, P. Léger, Application of NLFM models to predict cracking in concrete gravity dams, Journal of Structural Engineering 120 (4) (1994) 1255–1271.

[5] M. Leclerc, P. Léger, R. Tinawi, Computer aided stability analysis of gravity dams – CADAM, Advances in Engineering Software 34 (2004) 403–420.

[6] M. Ftima, P. Léger, Seismic stability of cracked concrete dams using rigid block models, Computers & Structures 84 (2006) 1802–1814.

[7] E. Bretas, J. Lemos, P. Lourenço, A DEM based tool for the safety analysis of masonry gravity dams, Engineering Structures 58 (2014) 248–260. doi:10.1016/j.engstruct.2013.10.044.

[8] S. Casolo, Modelling in-plane micro-structure of masonry walls by rigid elements, International Journal of Solids and Structures 13 (41) (2004) 3625–3641.

[9] S. Casolo, F. Peña, Modelo de elementos rígidos para el análisis de estructuras de mampostería, Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería 21 (2) (2005) 193–211.

[10] S. Casolo, F. Peña, Rigid element model for in-plane dynamics of masonry walls considering hysteretic behaviour and damage, Earthquake Engineering and Structural Dynamics 36 (8) (2007) 1029–1048. doi:10.1002/eqe.670.

[11] E. Bretas, Desenvolvimento de um modelo de elementos discretos para o estudo de barragens gravidade em alvanaria, Ph.D. thesis, Universidade do Minho, Portugal, 2012, in Portuguese.