Nowadays, dam safety is becoming more relevant in our society due to the importance of its functions (power generation, water supply, flood control) and the severity of the consequences in case of a serious breakdown. The "safety of hydraulic infrastructures" is one of the priorities according to the RD Spanish State Program which is oriented to social challenges.
Thanks to recent improvements in the modelling of dams as well as the evolution of computational resources, it is possible to perform more detailed analysis. However, most of commercial software is devoted to more common problems (e.g. buildings, bridges) without considering specific aspects of dams: concrete ageing process, joints during the construction process, contact with the ground or uplift pressure.
The development of a specific application for dam engineering, both for construction phase and operating period, is the main objective of this project. The thesis presents the computational tool, which is based on finite element formulations and solves thermomechanical problem using weak coupling, designed for analysing the operating period.
The proposed software is validated with a real case: La Baells dam. The computational results have shown excellent agreement with the obtained data during monitoring process.
This master thesis would not have been the same without the help of many people, thus I would like to dedicate them few words.
Firstly, I would like to thank Prof. Eugenio Oñate for giving me the opportunity of becoming a member of CIMNE. I would also like to thank my supervisors M.Sc. Fernando Salazar and Dr. Antonia Larese for their guidance and support during the writing of this thesis. I would like to make these thanks extensive to all the members of Room C6 and C4, particularly to Ignasi, for their wide advices and support. I would also like to thank to Dr. Josep Maria Carbonell and Dr. Jose Manuel González for his patience in my beginning with KRATOS and all provided geometries of dams. I can not forget to mention my friends and colleagues Tomás, Ruben and Vicente, with whom I spent so many hours.
I would like to also give special emphasis to my family, especially my parents and brother, because despite the distance, have always given me their unconditional support.
Finally, I would like to thank to all people that shared their time in one way or other with me; lifelong friends; Ana, Fernando, Carlos, Leo, Zito, classmates; Samu, Marc, Eric and Marta and flatmates; Jorge and Martin, whom have helped and supported me not only during this work, but also during the whole master course.
This work has been funded by the Ministry of Economy and Competitiveness under Retos Programme through the grant, RTC201537945  ACOMBO.
Water is a basic resource for all humans. Nowadays, the earth is not still affected by a lack of water but, according with some researches, this situation will change in a short period of time [1].
Only 1% of the water in the Earth is freshwater (mainly in rivers, groundwater reserves, lakes and reservoirs, Figure 1 ^{1}). However, only the reservoirs can be handled by humans. For this reason, dams play an important role in our society.
Figure 1: View of Kariba Dam, the world's biggest reservoir. 
Dams are the infrastructures that allow accumulate water with the aim of supply the population, generate energy and control the floods before it arrives to the sea through rivers.
The main objective of a dam is to resist the solicitations that the dammed water generates. This water not only induces a significant mechanical solicitations but also thermal loads due to the gradient of temperatures. Some researches have proved that the influence of the temperature in the global structural response can not be neglected [2].
(^{1}) http://www.worldbank.org/en/region/afr/brief/thekaribadamrehabilitationprojectfactsheet
Dam safety is becoming more relevant due to the importance of its functions (power generation, water supply, flood control) and the severity of the consequences in case of a serious breakdown.
In Spain, the majority of the dams were built during last 50 years [3]. For this reason, companies and institutions that are in charge of these structures are worried about the structural response of these critical infrastructures, which is one of the priorities inside the RD State Program (Programa Estatal de I+D+I Orientada a los Retos de la Sociedad).
These worries are based on past experiences. The Vajont disaster (Italy)(Figure 2 ^{1}) is one of the most important and well known. In 1963 an important landslide produced a tsunami which caused 2,000 casualties, however the body of the dam resisted, and today, continues standing [4].
Figure 2: View of Vajont Dam one day after the landslide. 
The rainfall generated by Typhoon Nina in Zhumadian (China, 1975) caused the worst disaster originated by dam break in the history. The number of casualties reached the amount of 171,000, and 11 millions of people had to leave their homes. The main factor of this disaster was the extreme conditions generated by the collision of Typhoon Nina and a cold front. The project flood was calculated using a returnperiod of 1,000 years but that day the flood was the equivalent to a returnperiod of 2,000 years [5].
During the last 10 years, at least 15 disasters happened around the world, the most important were: Shakidor Dam (Pakistan) in 2005 with 70 casualties, Situ Gintung Dam (Indonesia) in 2009 with 98 casualties, KyzylAgash Dam (Kazakhstan) in 2010 with 43 casualties or Köprü Dam (Turkey) in 2012 reaching 10 casualties.
In this context, the aim of this work is to develop a tool that simulates the real behaviour of dams and predicts its response.
(^{1}) https://en.wikipedia.org/wiki/Vajont_Dam
The main objective of this thesis is to develop a computational tool for structural verification of dams, including the possibility of thermomechanical analysis.
Since the software is specific for dams, the possibility to assign some adhoc boundary conditions is implemented, such as the Bofang temperature conditions [6] and the possibility to vary the water level through an input table data.
Another important objective of this thesis is to check whether accurate results can be obtained with tetrahedral elements. Dam structural calculations are typically performed with hexahedral meshes, which are difficult to generate for complex geometries. The use of tetrahedral meshes avoids these problems.
The present thesis is divided in six parts. The first one is the Introduction where the motivation and the objectives of the thesis were presented.
Chapter 2 is dedicated to literature review. In this chapter, a general guideline for the use of numerical models in dam engineering [7] are presented, focusing in recommendations for new dams and dealing with particular issues on existing dams.
In chapter 3 the proposed methodology for solving the problem is presented. This chapter deals with the mechanical, thermal and thermomechanical formulation. The derivation of the weak formulation is stated as well as the coupling between thermal and mechanical problem.
Chapter 4 is devoted to introduce the KRATOS environment and the Dam Applica tion. KRATOS is a framework for building multidisciplinary finite element programs. It provides several tools for easy implementation of finite element applications and a common platform with effortless interaction between them. On the other hand, Dam Application is a new specific application for structural verification of dams. The main implementations, how it works internally and its dependencies of other applications are presented.
Chapter 5 is dedicated to numerical results. La Baells arch Dam was selected as a case study: a numerical model was built and computed via the new Dam Application, and the results were compared to; a) the actual dam behaviour, as measured by the monitoring devices, and b) numerical results obtained with the software COMET.
In Chapter 6 some conclusions about the work performed as well as future work are discussed.
The use of numerical tools such as finite element, boundary element and finite difference methods, has become standard practice in dam engineering. There is extensive literature about this subject and a specific document published by the International Commission on Large Dams (ICOLD) in his bulletin 155 [7].
The object of the computational activity is often missed during the process. Due to this fact, it is important to emphasise the different purposes of a numerical simulation of damfoundationreservoir:
In the numerical analysis of dams, there are two different scenarios depending on the stage: design of a new structure, or analysis of an existing one. In case of new dams, general recommendations for the process of numerical analysis in this field are presented, in this stage the two first purposes are dealt: prediction of the structural stability and predesign. In the case of existing dams, special attention to the information provided by the monitoring system and the dealing with the last two commented purposes must be taken. The main distinction between modelling new and existing dam lies on the fact that the numerical models can not be calibrated for a new dam since no information on the actual dam response is available.
The basic goal of the numerical analysis is to provide adequate answers to a set of relevant questions related to the performance of a dam. Past incidents and failures have confirmed the necessity of this type of studies in the design phase considering all possible scenarios. Several aspects must be taken into account when considering the application of numerical models for design and analysis dams.
The size and refinement of the finite element grid
It is important to remember that a higher finite element density, and correspondingly finite elements of smaller size, is necessary where higher gradients are expected. These high gradients can appear in areas where there are geometrical singularities (e.g. sharp corners), as well as in relation to the loading condition.
The domain dimension around the dam is another important issue regarding the numerical model. This question takes special relevance when dynamic analysis is performed but can not be ignored for static or quasistatic analysis.
The SaintVenant principle states that a disturbance to a given model has an influence on the stress level in a surrounding domain which has equal to the dimension of the disturbance itself. For static loads and linear analysis, this principle can be applied. Considering the dam as the disturbance, the extension of the surrounding foundation should be at least equal to a characteristic dimension of the structure.
Generally, 2dimensional models are adopted for embankment dams or straight plant masonry/concrete dams; however, the application of this simplification to other typologies as: curved plant masonry/concrete dams or arch dams can lead to problems. In these typologies of dams the use of 3dimensional models are needed.
The inclusion of discontinuities
The consideration of discontinuities (e.g. construction joints, interfaces, cracks) demand a nonlinear finite element modelling. A preliminary linear elastic analysis in order to check the stress field is always a good idea. With this information it is possible to activate the kinematics (opening or sliding) of the discontinuities. In case of concrete arch dams can not avoid the inclusion of discontinuities.
The definition of the initial stressstrain state
This is one of the most complex problems to be solved for existing structures, not only because of the lack of information on construction but also due to the lack of knowledge on the thermoviscomechanical parameters which govern the phenomena.
In traditional numerical modelling, the dead weight of the structure is simulated as an instantly applied load acting on the complete structure to obtain its final conditions at the boundaries. In practice, the dead weight and the stiffness is applied progressively according with the dam construction.
With regards to concrete dams, the consideration of dead weight during the dam construction may be important to evaluate the stressstrain state caused by the generation and dispersion of hydration heat in concrete block masses which are not totally free to contract during the cooling process. In arch dams, the simulation of contraction joints injection could be more significant.
Treatment of consolidation of foundations
The evolution of the effects of this type of process (jetgrouting, cement injection) are quite difficult due to the physicalmechanical characteristics of the rock foundations. Field tests may be useful to evaluate the efficiency of the consolidation works. This information also can be used to prepare the numerical model.
In the design phase numerical analysis plays an important role. The design affects the final cost of the project but also the maintenance and the risk are factors that must be taken into account.
Mainly in the design of dams two tools are used: past dam engineering experience (behaviour of existing dams) and the analytical and numerical models. When these conditions are met computational models provide a safe baseline for design of new dams, based on past experience and knowledge with the support of associated design criteria.
In the design of a new dam, most of the information is generally conventional, there is no much data that can be used as projectspecific reference. Due to this lack of accurate information, this situation leads to an iterative design process with successive definitions of the project characteristics (dimensions) until a satisfactory degree of optimization is reached.
It is recommendable to perform an analysis in the linear elastic range, even for extreme loading conditions. The goal of this analysis is to identify possible damaged zones.
After carrying out linear analysis, it is recommendable to perform a nonlinear study under extreme loading scenarios (severe seismic ground motions) to obtain better approximations of dam behaviour. Nonlinear models are not easy to interpret. The main difficult comes from the stress transmission when the strength of the material is reached. This procedure can be corrupted for several reasons:
Nonlinear finite element models are a useful tool to capture nonelastic mechanical properties.
The appropriate framework for the numerical analysis process consist of:
Same procedure is suitable for both the design of a new structure and the analysis of an existing one.
Numerical models can produce a relevant picture of the state of stresses and the displacement field of the structure. Numerical models are essential to comprehend the behaviour of the dam. These models are widely used during the optimization process.
Figure 3 shows the different steps in the numerical model.
Figure 3: The four steps of a numerical model. 
In many cases this process is iterative. The evaluation of the first case often leads to modify the model, the imposed boundary conditions... when calibrations are possible. In this process the experience of the professionals plays an important role, they must be aware about the assumptions and simplifications used. This is particularly important in nonlinear models since the uniqueness of the solution is not guaranteed.
Preliminary designs are based mainly in past experience of the designer with respect to similar dams built. These designs usually have similar conditions: geological, topographical, hydrological and material parameters.
In case of having similar parameters and conditions a simple model can be suitable. On the other hand, when the characteristics are unprecedented (quite often in arch dams) a too simple approach may be insufficient, and should be considered as only initial step. When more information is available a more realistic model is performed.
The numerical analysis of new projects is typically associated with the following challenges:
In the case of extreme scenarios such as Maximum Credible Earthquake or Probable Maximum Flood, more sophisticated models may be necessary to explain this complex behaviour. Clear examples of this necessity can be transient underground seepage, cooling of new concrete, contact between two blocks or two types of materials and dam behaviour during earthquakes.
The aim of nonlinear finite elements is the representation of nonlinear stress behaviour of soils, concrete, rock foundation, cracks and joints. It is not necessary a global nonlinear analysis, since location of nonlinearities are known (e.g. damfoundation contact interface). Another important consideration in the case of global nonlinear analysis is that the results are nonunique due to they are dependent on the loadhistory, the model, and the numerical solution method.
The definition of acceptable safety criteria always requires a correct selection of design parameters since the solution is sensitive to modelling assumptions. Most of safety guidelines are not addressed to establish acceptable margins through the interpretation of finite element results. These factors of safety are sensitive to the calculation techniques.
The stresses are dependent on model discretization and mesh size. A clear example of that can be a heel of a gravity dam subjected to linear elastic analysis. The whole area related to the heel is mesh dependent. There are different techniques to avoid these problems, for example the use of adaptive techniques with the aim of refining the mesh in the affected area. Another possible solution could be the introduction of nonlinearities in this area.
Following, the particular issues on existing dams modelling are dealt.
Many of the considerations commented in previous section can be applied to this section. Nonetheless, this section is focused on the relation between the behaviour of the structure and the need of periodic assessment with understanding of the nature, source and relevance of the uncertainties involved.
One of the basic aspects related to existing dams is the surveillance. There are three important phases in the surveillance of dam life: during construction, during first impounding of the reservoir and during the operating period. Monitoring activity is an important ingredient to validate design assumptions, and also can be used to calibrate mathematical models or to detect abnormal deviations for longterm behaviours [8].
Construction
All data provided by monitoring devices are very useful for checking the critical design parameters (e.g. the foundation deformability or hydration heat) and to make corrections during construction phase.
First impounding of the reservoir
In this phase some monitored parameters concerning the damfoundationreservoir response to environmental loads with predicted values from design are compared. In case the differences are significant, it is a must the investigation of the origin of these differences and to take corrective measures if needed. This phase is also important in the mathematical model calibration.
Operation
The measurements from the monitoring systems and visual inspections carried out periodically are processed and compared to predicted values. In case of longterm operation, it is also important to check the influence of ageing process in the structure. In the case of detecting an important ageing process, the model must be calibrated. Another important mission of this monitoring process is the detection of behaviour changes and the study of these changes.
The obtained data are crucial in the interpretation of dam behaviour. However, the measurements can be under suspicion due to uncertainties linked to the installation, calibrating process or its conservation. In fact, wrong measurements can lead to incorrect decisions, due to this, some dam engineers can be sceptic about their importance [9].
The combination of statistical models and registered data can lead to better understanding of the dam behaviour.
Dams are designed for long time and their require an important investment not only during the construction period but also during the operating period. The surveillance of dams is crucial to avoid future disasters. A failure of a large dam can generate unacceptable casualties and important damages.
Due to above commented facts, there is a necessity to make forecasting models to predict the dam behaviour. The data processing of these predictions with the real data give the possibility of a quick response in case of issues.
Generally, forecasting models are of the following form [10]:

(1) 
where
The independent variables are:
There are four different ways for computing the forecasting indicators:
These methods are based on mathematical models. The quality of the forecast depends mainly on the capability of the mathematical model to describe the physical reality, the quality of the numerical solution, the knowledge of the parameters describing the materials , the simplifications introduced in the model and the knowledge of the independent variables at the time of observation. They do not need previous measurements, and usually is used for a first reservoir impounding or for the first years of operating.
A set of influence functions corresponding to expected variation of each indicator is chosen. Each component of each function is multiplied by a parameter or unknown coefficient, and these coefficients are defined with previous measurements. It is necessary the use of an a priori analytical form of the indicators (based on experience). The calibration process is performed through the minimization of a given norm. This technique can be applied to different phenomena.
These methods combine the above introduced models. By one hand, the analytical form of the indicators is obtained through the deterministic approach, but each single term of the function is multiplied by an unknown coefficient. The calibration is carried out as in statistical approach.
This methodology [11] is relatively new and is inspired by neural networks. The neural networks are composed of many simple intercommunicating units, or neurons, working in parallel to solve a given problem. It expands as the network learns to operate, based on a set of training data provided by users. Nowadays, several types of neural techniques exist, the most important one in the back analysis is the multilayer feedforward and backpropagation algorithm. These methods have lot of potential.
The surveillance of a dam is mainly performed by visual inspections carried out by qualified personal and by processing data obtained from monitoring devices.
Most of dangerous events as sinkholes, settlements, cracks, concentrated seepages, generally can not be detected by monitoring devices and must be detected through visual inspections. When an abnormal behaviour is detected, its evolution can be surveyed and analysed on the basis of the data collected by a monitoring system suitable installed. The parameters provided by the monitoring devices can be grouped in two categories: environmental actions and the physicalmechanical quantities describing the response of the damfoundationreservoir.
Inside this group it is possible to find the following quantities: reservoir water level, air temperature, water temperatures at different depths in the reservoir, solar radiation, seismic actions, rainfall...
In this group it is necessary to distinguish between different typology of dams:
It is important to remind that the monitoring devices tends to modify and disrupt their local value, especially in quantities like strains and stresses.
Thanks to the important growth of electromechanical devices and advances in topic of the transmission data, this information each time is more relevant and easy to access.
The basic idea of any structural identification method is to make a numerical model behaviour as close as possible to the observed one.
The structural identification methods can be based on static and dynamic monitoring data.
These methods are basically those already shown in previous section 2.2.2: the calibration step usually performed to set up forecasting models may be seen as identification process which allows evaluating the best estimate of the physical/mechanical parameters of the damfoundation system.
Dynamic methods have taken more importance in last decades. These methods are based on dynamic excitation. The aim of the dynamic structural identification consist in the evaluation of the dynamic features (e.g. natural frequencies, modal shapes) according to excitation (input) and response (output) of the system. This process is performed through dynamic campaigns. The structural response can be measured in terms of acceleration, velocity or displacement vs. time. Depending on the available devices, the measurement is carried out in one term or other. According to the measured parameter a different range of frequencies must be used (e.g. velocity measurements requires low frequencies).
The excitation of the structure can be obtained in a different ways depending on the available resources, the importance of the structure...
The selection of a proper method is very important so that no damage are induced on the structure during the test. To this aim is highly recommended a previous evaluation to predict the amplitudes on the model.
The selection of an adequate mathematical algorithm for identification is very important, and the following aspects must be taken into account:
Frequency and/or modal shapes may be used directly for matching and error minimization, as already was mentioned before for time domain identification.
The extension of a dynamic analysis to a nonlinear field must be done carefully and considering the following differences:
Summing up, structural identification not only provides a powerful diagnostic tool but also a refined system to support a deep understanding of the dam behaviour.
In this chapter the used methodology is presented. Due to the nature problem the structure of the chapter is divided in three sections: mechanical, thermal and thermomechanical problem.
The first section deals to the mechanical problem focusing in Linear Elasticity Theory, a branch of General Elasticity Theory. In this section the governing equations of the mechanical problem as well as the weak formulation and the numerical integration techniques are presented. In second section the thermal problem is stated. Same structure than in previous section is followed. Finally, the thermomechanical problem with a weak coupling is introduced.
Many of the concepts discussed below are dealt extensively in the book of Prof. Eugenio Oñate, Structural Analysis with the Finite Element Method. Linear Statics [12] [13],as well as in the book of Prof. Xavier Oliver and Prof. Carlos Agelet de Saracibar, Mecánica de medios continuos para ingenieros [14].
Usually, a mechanical problem can be defined in the following way; given an initial geometry and a loading field, analyse the response according to a stiffness and strength. For solving the problem is necessary to state the basic equations.
Depending on the constitutive equation the solid behaviour is different. The following models describe how a solid responds to an applied force:
This study is focused in Linear elasticity since the approximation that it provides is good enough for the purpose. Most of engineering problems can be solved using this type of hypothesis.
Linear elasticity is the mathematical study of how solid objects deform and become internally stressed due to prescribed loading conditions. Linear elasticity is a simplification of the more general nonlinear theory of elasticity and is a branch of continuum mechanics.
The fundamental "linearizing" assumptions of linear elasticity are[14]:
The solution for elastic problems can be stated through two different approaches depending on the unknown variable;
In this thesis, a displacement approach is used, but in the cases where the contour stresses are known, the use of BeltramiMichell approach is reasonable. Nonetheless, there are mixed formulations that allow to solve the problem without the necessity of knowing all stresses [15].
In direct tensor form, the governing equations for a linear elastic problem are [14]:

(2) 

(3) 

(4) 
The boundary conditions and the initial conditions are

(5) 
Figure 4: Scheme of generic problem in linear elasticity. 
Regarding the analysis type, there are three important groups:

(6) 
Depending on the analysis type, the system to be solved is different. The most common analysis types are: quasistatic and dynamic.
Using a Finite Element Method and a quasistatic analysis, the system reads as

(7) 
where K is the stiffness matrix, u is the vector of displacements and f is the vector forces.
In case of dynamic analysis, the system becomes more complex and reads like

(8) 
where M is the mass matrix, C is the damping matrix and, and are acceleration and velocity vector, respectively.
In this thesis, due to the considered loads are not influenced by the time, the analysis type is quasistatic.
The Finite Element Method (FEM) is a numerical technique for finding approximate solutions to boundary value problems for partial differential equations. It uses subdivisions of a whole problem domain, and variational methods to solve the problem by minimizing an associated error function.
The main concept is connect many simple element equations over many small domains to approximate a more complex equation over a larger domain. Although continuum structures are always threedimensional, if the proper simplification hypothesis are fitted, one can accurately describe the behaviour of a structure by means of uni or bidimensional mathematical models, like in earth dams, tunnels, etc.
Next, the generic stages for the analysis of a structure through the Finite Element Method are presented. Assuming a quasistatic analysis, the equilibrium equation in strong form can be stated as

(9) 
In solid mechanics, the six stress components can be computed through the six components of strain () which are computed from the displacements

(10) 
To obtain the weak form an arbitrary weight function vector 11 is introduced

(11) 
and the equation is integrated over whole domain

(12) 
after applying the divergence theorem is obtained

(13) 
introducing the concept of virtual strains that can be stated as the derivation of virtual displacements

(14) 
which is the threedimensional equivalent virtual work statement.
The virtual displacements and the virtual strains are interpolated in terms of the virtual displacement values in the standard form

(15) 
Substituting and simplifying the virtual displacements

(16) 
where B and N are the strain matrix and the shape function matrix, respectively. Using the stresses expression

(17) 
finally the following equation is obtained

that can also be expressed as

(18) 
where

(19) 
is the elastic stiffness matrix of the element, and

(20) 
the nodal vector force. The nodal vector is composed by the body forces and tractions forces, but may also be other contributions. These concepts and the way to compute the stiffness and the vector force are detailed in Appendix A.
To finish, the global system is constructed

(21) 
and solved for the unknown variable, the nodal displacement vector. Once it is obtained, it is quite straightforward to calculate the strains and stresses at each element as well as the reactions at the nodes with prescribed displacements.
Integration schemes are methods to integrate differential equations. There are two groups of integration schemes; explicit and implicit.
The explicit schemes use the previous step solution () to compute the current step solution (), these types of schemes are conditionally stable and require that the time step size be less or equal than a critical time step. The critical time step depends on element size and maximum wave speed for element material. In some cases, this value has to be too small and requires high number of steps to reach the reality, driving to high computational effort.
The implicit schemes use previous step solution () and current step solution () to compute the current solution. Generally, an implicit solution is unconditionally stable [16], and the used time step can be one or two orders of magnitude larger than in an explicit schemes. However, the accuracy of the implicit schemes deteriorates as the time step increases relative to the period of response of the system.
In this work two well known numerical integration schemes are presented: Newmark method and Bossak method, which is an extension of Newmark method.
Newmark Method This method is widely used in numerical evaluation of the dynamic response of structures and solids [17]. Originally, this method is from 1959 but some improvements have been developed. In matrix form the linear dynamic equilibrium equations can be written as

(22) 
Using truncated Taylor's series as two additional equations expressed in the following form

(23) 
and assuming linear acceleration within the time step

(24) 
it is just necessary to substitute Eq.24 into Eq.23 to produce Newmark’s equations in standard form

(25) 
It is known that values of =0.5 and =0.25 make an unconditionally stable method. The use of values of higher than 0.5 introduces errors associated to numerical damping.
Bossak Method
The Bossak method is the extension of the Newmark method. The acceleration is taken prior to . The method can successfully replace the Newmark method in all cases, and can be stated as

(26) 
The use of =0 leads to Newmark method. The Bossak method is characterized by a good damping in high frequencies range.
The used strategy is typically from nonlinear problems, but also can be applied to linear ones. One of the main features of linear elasticity is the uniqueness of solution, however in nonlinear problems there is a lack of uniqueness of the solution. In fact, in many cases the correct solution is path dependent, so it depends on the path followed.
In matrix form a nonlinear problem (given a prescribed boundary conditions and having applied the particular time discretization) can be stated as

(27) 
where u represents the vector of unknowns, K the 'stiffness matrix', and f the force vector, both K and f depend on the solution u. Because of u is unknown, this system is solved by some iterative algorithm, for instance the NewtonRaphson method.
Figure 5: Example of an incrementaliterative method. 
Figure 5 shows an incremental iterative method. In incremental iterative methods there is a part of prediction and a part of iteration. The iterations continue until reach some prescribed tolerance. There are different types of convergence criteria depending the case.
Before each solution, the NewtonRaphson method evaluates the outofbalance load vector, which is the difference between the restoring forces (the loads corresponding to the element stresses) and the applied loads. Then, a linear solution is performed using the outofbalance loads, and it is checked for a prescribed convergence criterion. If the convergence criterion is not satisfied, the outofbalance load vector is reevaluated, the stiffness matrix is updated, and a new solution is obtained. This iterative procedure continues until the problem converges.
The Picard's method is also a famous iterative method for solving nonlinear problems, it is computationally cheap but the convergence is very slowly, on the other hand, NewtonRaphson requires higher computational cost but the rate of convergence is higher.
The thermal problem is governed by classical heat equation. In dam engineering, the thermal problem takes on special importance due to the temperature variation in different seasons. The temperature field inside of the dam can be predicted solving the diffusion problem. Once the temperature field inside of the dam is obtained, its gradient can be computed.
The energy equation for a solid in a vector form can be stated as [18]

(28) 
where is the density, is the specific heat, is the thermal conductivity and is the rate of heat generation.
Neglecting the heat generation term (since the system does not generate its own heat per unit volume), and arranging terms, the equation reads as

(29) 
the above equation is known as heat equation. The thermal diffusivity () is obtained dividing the thermal conductivity between the density and the specific heat. Finally, the equation reads as

(30) 
Once the governing equation is obtained, the boundary conditions and the initial conditions are introduced. Boundary conditions are given by

(31) 
where is an applied flux, is the convective flux and is the radiative flux, which are dependent on the convective heat transfer coefficient and the effective radiation heat transfer coefficient, respectively

(32) 
The initial condition on the temperature is given by

(33) 
The governing equation of the thermal problem in a strong form has been presented, now, the way to obtain the weak formulation of this equation is introduced. Due to it is a linear problem, the material coefficients does not depend on temperature. These assumptions can be taken into account since the variation of the temperature is supposed not too much high. In other type of problems in which the variation of the temperature is high, the problem must be considered temperature dependent.
Using Eq.28, neglecting the heat generation term, multiplying by a test function , and integrating over the element, the following expression is obtained

(34) 
after, the divergence theorem is applied

(35) 
Applying Fourier's law, the rate of flow of heat energy per unit area through a surface is proportional to the negative temperature gradient across the surface , the following expression is obtained

(36) 
Now, assuming that the time dependence can be separated from the spatial variation, the finite element approximation can be stated as

(37) 
where is the vector of shape functions and is a vector of nodal temperatures. Finally, making the substitution of according with weakform Galerkin, the following expression is obtained

(38) 
in matrix form reads as

(39) 
where C is the specific heat matrix, K is the conductivity matrix, Q is the vector of heat fluxes, is the vector of nodal temperatures and is the time rate of the nodal temperatures.
The ordinary differential Eq.39 should be integrated with respect to time to obtain a transient response. In general, it is not possible to integrate these equations analytically, so they are approximated in time to obtain a set of algebraic equations in terms of the nodal temperature.
The most popular used method for diffusion problems is called family and it consists of the following approximation

(40) 

(41) 
and rearranging terms the equation can be expressed as

(42) 
The choice of produces the forward euler scheme, an explicit method which is conditionally stable

(43) 
and may be rearranged to clearly show its explicit nature

(44) 
Eq.44 requires that matrix C be easily invertible for an effective implementation. The reduction of the matrix to a diagonal or "lumped" form can be accomplished by various procedures, for instance rowsum [19,29]. ForwardEuler method is conditionally stable and the stability condition is governed by the thermal diffusivity of the material, and the finite element space [20]. Euler method is firstorder accuracy in time. The main attraction of an explicit scheme is the fact that a matrix solution is not required and computer storage requirements are minimal.
When , a family of implicit methods are produced. Depending on the used value there are different schemes:
Inside the code is an input parameter which leads to one scheme or other. Computations have been performed using , backward Euler scheme, and reads like

(45) 
Next, the solution strategies for this type of problems are introduced.
The discretization in time and space has already presented, now, some strategies to solve the problem are presented. The linear problems do not require iterative solution. These problems just require a matrix solution.
Nonetheless, the nonlinear problem using a backward Euler is presented, and the system reads like

(46) 
In case of dealing with nonlinear equations (Eq.46), the use of predictorcorrector, extrapolation methods, and quasilinearization can often reduce the computational effort at each step without reducing the accuracy of the method.
A quasilinearization is good choice for mild nonlinearities and modest time steps. The system reads like

(47) 
which is now a singlestep (a non iterative method). Another choice is the use of an extrapolation procedure. Usually, it gives better solution than quasilinearization and allows a larger time step. The extrapolation procedure takes the following form

(48) 
where

(49) 
In case of strong nonlinearities, these techniques do not offer accurate results. In these cases, the coupling of a forward Euler predictor with a corrector is a viable option [21][22].
Two different problems need to be solved: a mechanical problem and a thermal problem. There are different approaches to do so, but the most common are probably two: the fully coupling and the oneway coupling.
The fully coupling is usually reserved to solve problems in which both the mechanical and thermal components can affect each other in a similar way. In order to deal with this kind of problems one needs to solve the two components at the same time, taking into account the interaction between them. Among others, this approach is used in the solution of metal forming deposition and additive manufacturing. In this sort of situations, the material properties change dramatically in time and so it is crucial that both components are fully coupled so as to properly capture the real behaviour of the medium. In the case of additive manufacturing, for instance, the higher temperatures reached during the process of binder generate a totally different material behaviour as compared to that encountered in the cooling process.
On the other hand, the oneway coupling technique is devoted for those problems in which the interaction between both parts exists, but it is not symmetric. In essence, one of the components has a strong influence on the other, but the contrary does not occur. This is the case of dams, where the mechanical problem is really affected by the temperature changes, but the thermal problem is almost insensible to the deformation of the geometry.
Next, the main features of the latter coupling are introduced.
First, the diffusion problem is solved, once it is solved, instead to consider an elastic problem, a thermoelastic problem taking into account the influence of the temperatures through the constitutive law is considered.
The main difference between linear thermoelasticity and linear elasticity is that the process is not considered isothermal any more. The temperature value changes along the time , it means

(50) 
Nonetheless, the hypothesis of adiabatic process is still considered.
The thermal contribution are added to the stresses as

(51) 
where is the constitutive matrix, is the thermal properties matrix and is the reference temperature field. Usually, the difference between temperatures is written as .
Eq.51 can be inverted to obtain the relation of the strains and it reads as

(52) 
The definition of thermal expansion coefficient is

(53) 
In case of isotropic material, the thermal expansion coefficient becomes a scalar.
Taking the assumption of working with an isotropic material (case study), the relations between thermal and nonthermal for stresses and strains are presented

(54) 

(55) 
KRATOS is designed as an OpenSource framework for the implementation of numerical methods for the solution of engineering problems. It is written in C++ and is designed to allow collaborative development by large teams of researchers focusing on modularity as well as on performance [23] [24].
The KRATOS features a "core" and "applications" approach where "standard tools" (such as linear algebra or search structures) come as a part of the core. They are available as building blocks in the development of "applications" which focus on the solution of the problems of interest. Its ultimate goal is to simplify the development of new numerical methods.
Figure 6: KRATOS' Logo. 
A new application called Dam Application () has been created due to an specific necessity. This Application is focused in dam engineering. Its main features as; how the application works or what are the dependencies are described below.
has dependencies of Solid Mechanics Application, ConvectionDiffusion Application and in less measure of Poromechanics Application. Figure 7 shows the basic flowchart of the application.
Convectiondiffusion Application, as its own name shows is an application for solving convectiondiffusion problems, it is in charge of Dr. Pablo Becker. Solid Mechanics Application is devoted to solve all type of problems related solid mechanics, it is in charge of Dr. Josep Maria Carbonell. Poromechanics Application solves porepressure problems and the responsible is Ing. Ignasi de Pouplana. I would like to thank all of them for allowing me to use their code.
Figure 7: Dam Application's dependencies. 
Apart from computational features, a basic but powerful interface where boundary conditions, materials and elements are taken into account has been created. This interface works in GiD, which is a universal, adaptive and userfriendly pre and post processor for numerical simulations in science and engineering [25].
Figure 8 shows the flowchart for solving a generic problem.
Figure 8: Flowchart for solving a generic problem. 
borns with the aim of solving a real necessity in world of dam engineering. It provides userfriendly interface and specific conditions for dam problems. Internally, it solves a thermomechanical problem using a weak coupling.
For solving the thermal problem, the Convectiondiffusion Application is used. Since the convective part is not considered in the problem, all its contributions are set to zero, leading to a diffusive problem. After solving the thermal problem for a given time, the mechanical problem taking into account the influence of thermal contribution through the constitutive law must be solved.
Following lines are devoted to present . The structure is the following; first, the interface is introduced, then the specific boundary conditions of dams are commented, and finally the main files are presented.
The interface is a powerful tool to define all required information for solving the problem. The development of a correct interface leads to a better understanding and use of the software. A simple but powerful interface has been developed.
The interface is based on the classical problemtype of GiD (Figure 9)
Figure 9: GiD interface 
Figure 10 shows the divisions of the main menu. Concretely, the menu is divided in six points:
Figure 10: Dam Application menu interface. 
The displacements and temperatures field can be assigned in this menu. The appearance of menu is shown in Figure 11.
Figure 11: Interface menu for applying displacements and temperatures. 
This section is devoted to the application of forces (Figure 12). Several types of forces have been implemented. The user must be aware that there are some limitations, for instance in a 3D view it is not possible to apply a line normal load since the normal direction can not be computed. Below, the different possibilities of loads that can be applied are introduced.
Figure 12: Interface menu for applying load conditions. 
In this section mechanical and thermal conditions related to the materials can be applied. These conditions are needed due to the Convectiondiffusion Application works with nodal variables. Gravitational or acceleration effects can also be assigned in this menu.
This section is dedicated to assign the discrete elements. For instance, in a 2D case it is possible to use linear triangles and quadrilaterals, and also high order elements. In 3D analysis, the possibilities are linear tetrahedrons and hexahedrons as well as high order elements.
It is worth to remark that the thermomechanical problem only can be solved by triangles or tetrahedrons since the thermal problem is only implemented for this type of elements. Facing a mechanical problem the user can use any type of element.
Material assignation is performed in this section. Figure 13 shows the different possibilities for this assignation. Some default materials are ready to use, although if is convenient, the user has the possibility to create new materials. Through the definition of the constitutive law the type of problem (mechanical, thermal or thermomechanical problem) is specified.
Figure 13: Interface menu for applying materials. 
To finish this section, the different tabs inside of the Problem Parameters menu are presented.
In this menu (Figure 14), the domain size of the problem, the analysis type (quasistatic or dynamic), the time scale and the temporal values are defined.
Figure 14: Interface menu: General Data. 
This menu (Figure 15) is the responsible to keep the conditions constant (Unmodified) along time or update them at some point (Table Interpolation).
Figure 15: Interface menu: Imposed Conditions. 
There are some special conditions as water load and Bofang temperature that are governed by other table Table Evolution Data, this table is explained in the following section.
For computing more than one case, it is possible to use this table (Figure 16). The needed values are the following ones: the month in which the computations are carried out (Bofang formulation), the evolution of the water level along the time, the temperature in the not wet wall upstream and the reference temperature.
Figure 16: Interface menu: Table Evolution Data. 
In this menu (Figure 17), some parameters about the mechanical solver can be set: number of maximum iterations until the convergence is reached, the tolerance of the DOF's and for the residual, the linear solver (direct or iterative) and its type, and finally the number of threads.
Figure 17: Interface menu: Mechanical Solver. 
This menu (Figure 18) is devoted for selecting the output parameters. According to interest outputs one or another are selected.
Figure 18: Interface menu: Postprocess. 
One of the most remarkable features of the is the possibility to apply the following conditions:
These conditions are wide used in dam engineering. In this subsection, how these conditions have been implemented and its relevance are commented.
In the upstream wet wall, it is possible to approximate its temperature to the temperature inside of the reservoir.
Up to some years ago, in many research articles, Stucky formulation [26] was used to provide a temperature field, but this temperature field is not too much realistic. Other approaches combine computational fluid dynamics and heat transfer models [27], however, these approaches are complex and computationally expensive.
Analytical models are other alternative to high computation approaches. These approximations are straightforward and its formulation is quite simple. Inside this analytical models, the Bofang formulation has taken important relevance [2]. The Bofang formulation provides a temperature field inside of the reservoir according to an atmospherics input parameters. It is worth to remark that depth of the dam must be higher than 30m, otherwise, the results can be inappropriate.
The proposal of Bofang formulation Bofang is depending on time and depth

(56) 
where
Once the general form of Bofang formulation Eq.56 has been presented, the contribution of each part is introduced.
The yearly mean water temperature () varies with the and may be expressed by

(57) 
where
The way to obtain these parameters is explained in following lines.
In the case of the surface temperature (), it has to be computed in the following way. In regions where the surface of the reservoir does not freeze in winter, the yearly mean water temperature at the surface is computed by the yearly mean air temperature plus an increment due to solar radiation. In cold regions, the yearly mean air temperature must be modified. The modification consist on set a 0C the temperatures lower than 0C. According with observed data .

(58) 
In the case of the bottom temperature (), these temperatures are computed as the mean of the mean value of the air temperature in the three months in the winter.

(59) 
The amplitude of annual variation of water temperature () is the highest at the surface of reservoir and decreases with the depth of water. It may be expressed as

(60) 
where is the amplitude of variation at the surface of reservoir. Normally, this amplitude can be computed as the mean of maximum () and minimum () monthly mean air temperature. Usually, in the Northern Hemisphere, these values belongs to the mean air temperature of July and January, respectively.

(61) 
In cold regions is computed as

(62) 
where is the increment of amplitude due to solar radiation.
The phase difference of annual variation of water temperature depends on the depth of water. In case of working in months, it is suggested to compute as

(63) 
To understand in a better way how this formulation works, an application example is shown. Figure 19 shows the temperature variation according the period of the year and the depth. The input parameters are described below

Figure 19: Application of Bofang Formulation. 
The temperature at the surface and just below it suffers high variations depending on the period of the year while in deeper levels remains much more constants without barely variations.
The application of this formulation in a real case is presented hereafter. The case study is La Baells dam and the input parameters are

Figure 20 shows the temperature field in two different seasons when the reservoir is full (). In January (Figure 20a) the variation barely reaches 3C while in July (Figure 20b) it is possible to observe that the variation between the surface and the bottom is amount 12C.
(a) January  (b) July 
Figure 20: Application of Bofang Formulation at La Baells dam. 
This condition has been implemented in two different program languages: and . The first language is used to write the input file (), the second one for updating conditions at each time step. There are two different ways to update the conditions: interpolate evolution (given two times, the evolution is computed according an interpolation law) or exact evolution (the condition only changes when the current time arrives to a determined time). The implementation of the developed code in C++ is stated in Appendix C.
A fluid in a resting state generates forces on the walls and bottom parts of the container and also on the surfaces of any object submerged. These generated forces are perpendicular to container's walls or any submerged object's surface regardless of the face's direction.
This pressure depends on the density of fluid and the reference height of the measurement. The absolute pressure can be computed as

(64) 
where
Many times it works with relative pressure, neglecting the atmospheric contribution. The application of this condition is performed by the application of nodal forces according to the water level and the coordinates of the nodes.
Below two application examples are shown. The scene of application is La Baells dam. Figure 21a shows the hydrostatic pressure (Pa) with 40m of water level, and outside of water the applied force is equal to 0. In the case of Figure 21b, the water level is equal to 93m, it means that the reservoir is full.
(a) Water level = 40 m  (b) Water level = 93 m 
Figure 21: Application of Hydrostatic Pressure at La Baells dam. 
The same procedure for coding this boundary condition than in previous case has been followed. The developed code in can be found in Appendix C. As in previous case, there are two ways of updating the conditions; exact or interpolated. It is worth to remark that this file also updates the conditions for uplift pressure.
Uplift pressure is a phenomenon very relevant in gravity dams, in arch dams this phenomenon is less important. In gravity dams, the uplift pressure has special importance due to the length of the base is relevant.
Uplift pressure can be defined as the force that is generated due to the seepage in saturated soils at the foundations. These forces act in the base of the structures. Uplift pressure can play an important role in the stability of the structure because it represents a destabilized force.
There are several approaches to compute the uplift pressure. Usually, to know this force, using the flow it is possible to draw the pressure's diagrams, or it is also possible to compute the force using numerical techniques.
The following approach of the uplift pressure has been implemented. In the upstream part the same force than in the lower part of the reservoir is assigned, and this distribution is decreasing until it arrives to 0 at the end of the length of the base. Figure 22 shows this approximation.
Figure 22: Scheme of uplift pressure approximation. 
The formulation to compute the law is

(65) 
The proposal formulation is related with the water level, it means that in the part that are in contact with the reservoir, the value of the uplift pressure is maximum and it decreases until arrive to 0. It is a must to provide the coordinates in the uplift direction of the upstream to perform the computations.
In arch dams this phenomenon does not take special relevance, and due to the arch geometry, it is a must the modification of the reference coordinate to apply the law properly making harder the application of this type of boundary condition for the user.
Figure 23 shows the application of hydrostatic and uplift pressure in a gravity dam. In the bottom of the upstream there is a matching in values between hydrostatic and uplift pressure.
Figure 23: Hydrostatic and Uplift pressure in a gravity dam. 
As it has been already commented in previous section, the implemented code can be found in Appendix C.
The or also known as is the file that contains all the information related to geometry and boundary conditions of our problem; the mesh (nodes and elements), the applied boundary conditions, the assigned materials, submeshes for updating conditions, and tables. On the other hand the file has information about General Data; domain size, , type of solver...
The file is written by file. This file is the responsible to translate all applied specifications during the preprocess in a script file. A small sample of the file is shown below.
Begin ModelPartData //VARIABLE_NAME value End ModelPartData Begin Table 1 TIME DISPLACEMENT_X 0.00000e+00 0.00000e+00 End Table Begin Properties 1 CONSTITUTIVE_LAW_NAME LinearElastic3D DENSITY 2.50000e+03 YOUNG_MODULUS 3.50000e+10 POISSON_RATIO 2.00000e01 End Properties Begin Nodes 1 0.00000e+00 4.00000e+01 1.00000e+01 2 2.00000e+00 4.00000e+01 1.00000e+01 3 0.00000e+00 3.80000e+01 1.00000e+01 4 0.00000e+00 4.00000e+01 8.00000e+00 End Nodes Begin Elements SmallDisplacementThermoMechanicElement3D4N 1 1 12658 12754 12598 13227 2 1 16536 16594 16471 16462 3 1 15255 15426 15541 15161 4 1 12848 12341 12013 12122 End Elements
In the case of the file, the responsible of the generation of the input file is . A small part of the file is detailed below.
## General Data –––––––––––––––––– domain_size = *GenData(Domain_Size,INT) NumberofThreads = *GenData(Number_of_threads,INT) time_scale = "*GenData(Time_Scale)" evolution_type = "*GenData(Evolution_Type)" delta_time = *GenData(Delta_Time) ending_time = *GenData(Ending_Time)
These files are not included in Appendix C due to its extension. In case of interest, these files can be found in https://kratos.cimne.upc.es/projects/kratos/repository/show/kratos/applications/DamApplicationKRATOS repository.
Once the input data have been generated, these are read by the main script. This script () is composed by different stages for solving the problem. In the following lines the different stages are presented.
The time is initialized and all necessary modules are imported, also the mechanical and thermal solver are defined.
The time parameters, scale of times and output variables are defined.
The file (that previously has been generated) is called, and the unknown variables are added to the solver, according the inputs. The is read and all information about mesh and conditions are extracted. The degrees of freedom are added to the solvers (thermal and mechanical) and the process info is finally set. The previous steps are filled with the boundary conditions, if necessary.
All functions that are involved in the solution of the problem are initialized.
Once everything has been initialized, inside of the temporal loop the resolution of the problem starts. First of all, the time variables and boundary conditions are updated, after the thermal problem is solved, then the mechanical problem with the thermal contribution is solved, and finally the results according the selected parameters in file are written.
Finally, all functions are closed and the elapsed time is printed.
The implemented code can be found in Appendix C.
This chapter is dedicated to numerical results. La Baells arch dam (Figure 24 ^{1}) was selected as a case study: a numerical model was built and computed via the new Dam Application, and the results were compared to: a) the actual dam behaviour, as measured by the monitoring devices, and b) those obtained with the software COMET.
The data used for the study correspond to La Baells dam. It is a double curvature arch dam, with a height of 102 m, which entered into service in 1976. The monitoring system records the main indicators of the dam performance: displacement, temperature, stress, strain and leakage. The data were provided by the Catalan Water Agency (Agencia Catalana de l’Aigua, ACA), the dam owner, for research purposes.
Figure 24: View of La Baells dam. 
(^{1}) http://www.bergactual.com/wpcontent/uploads/2009/04/373363.jpg
Among the available records, the study focuses on 10 variables: 6 correspond to displacements measured by pendulums (radial), and four to temperature.
The location of these pendulums is shown in Figure 25. The exact position of the pendulums is detailed in Table1. The origin of the coordinate system is the sea level and the generatrix of the archs.
Figure 25: Position of pendulums at La Baells dam. 
Coordinates (m)  
X  Y  Z  
P1DR1  31.60  135.10  610.00 
P1DR4  31.80  136.10  590.00 
P2IR1  31.60  135.10  610.00 
P2IR4  31.80  136.10  590.00 
P5DR1  90.00  101.27  628.00 
P5DR3  90.40  101.80  610.00 
Figure 26 shows the position of the thermometers and Table 2 provides the exact position of them.
Figure 26: Position of thermometers at La Baells dam. 
Coordinates (m)  
X  Y  Z  
6196D  105.95  90.03  619.00 
6193D  60.62  120.76  619.00 
5812I  28.78  140.87  581.00 
5611I  45.72  130.57  561.00 
The analysis is performed from 1994 until 2007. The information of the water level and the air mean temperature during these years are shown in Figure 27 and 28, respectively.
Figure 27: Hydrostatic load evolution 19942007. 
Figure 28: Air temperature evolution 19942007. 
The problem is solved in the following way. Firstly, the thermal problem according with prescribed boundary conditions (Bofang and the air temperature) is solved. The used time step for solving the thermal problem is month. Then, the mechanical problem taking into account the contributions of the thermal problem is solved.
Due to the thermal effects play an important role in the problem it is considered appropriate to start the resolution of the thermal problem in 1990 with the aim of obtaining a more realistic temperature field.
In all the simulations two different materials were used, concrete in the body dam, and foundation in the base. The material specifications are detailed in Table 3
Property  Concrete  Foundation 
Density  2,400  3,000 
Specific heat  982  950 
Conductivity  2.43  2.2 
Young Modulus  4.67e10  3.1e10 
Poisson  0.25  0.25 
Thermal expansion  1e05  1e05 
To simplify the input parameters, it was decided to compute the Bofang parameters for the period of study. To do that, a monthly mean during the study time is computed. Once the monthly mean is obtained, the input parameters to apply the formulation are computed. In Table 4, the monthly mean values of air temperature during 19942007 are shown.
Month  Temperature 
January  4.27 
February  5.22 
March  7.22 
April  10.42 
May  13.00 
June  18.27 
July  21.50 
August  21.92 
September  19.13 
October  15.60 
November  10.82 
December  5.38 
The input parameters for this period of study are

Two different Dirichlet boundary conditions need to be applied on the upstream dam face. The first condition that is applied is Bofang at the wet wall, but when the reservoir is not totally full, it is necessary to assign a value at the dry area. In this case, this temperature has been approximated adding 2 to the air temperature.
Downstream is also necessary to assign a temperature value. There is a vast literature about how to compute these values [2]. The used approximation was to add 2 to the air temperature. Another important data that must be introduced is the reference temperature. It is considered appropriate to set this value as the joints injection temperature, .
Figure 29 shows a scheme with all thermal and mechanical conditions.
Figure 29: Thermal and mechanical boundary conditions on the dam body 
(a) P1DR1  (b) P1DR4 
(c) P2IR1  (d) P2IR4 
Figure 30: Convergence analysis at La Baells dam 
A mesh of 86,000 nodes is used (Figure 31), although the mesh has converged before arriving to 86,000 nodes, the interest problem is a thermomechanical one and this problem is more restrictive with the element size in convergence criteria. The mesh is composed by 86,993 nodes and 323,955 tetrahedrons.
Figure 31: Tetrahedral mesh. 
The contact between both volumes (dam and foundation) is carried out just by sharing nodes.
The temperature analysis was performed using the information provided by the four thermometers previously presented. Figure 32 shows the comparison of temperatures during 2007 year between real data and the values obtained by Dam Application.
(a) 6196D  (b) 6193D 
(c) 5812I  (d) 5611I 
Figure 32: Temperature analysis at La Baells dam 
Figures 32a and 32b show the evolution of temperature during 2007 in the higher parts of the dam. The temperature variation in these thermometers between winter and summer is quite high due to its position. As was commented in Chapter 4, upper positions are highly influenced by climatic conditions. The opposite occurs in thermometers 5812I and 5611I shown in Figure 32c and 32d, respectively. In these parts the variation is quite low.
In general terms, the temperature field is well approximated by the used formulation, the greatest differences are located in summer months.
The displacement measurement is provided by the six pendulums previously presented. Dam displacements are measured in local axes: tangent and perpendicular to the dam axis at each location. They are termed "tangential" and "radial" displacements, respectively. In this work the radial displacements were considered, since their magnitude is greater and thus the effect of the measurement error is lower.
Since the numerical model results are provided in global Cartesian axes, they need to be rotated to be compared to the observed values. The angle of this rotation depends on the location of each pendulum. Moreover, the results suggest that some discrepancies might exist between the theoretical and the actual position of some local axes, as shown below.
Finally, a constant value was added to each radial displacement to account the displacement previous to the installation of the measurement device.
Taking all above mentioned into account the transformation of the displacement field is performed according the next formulation:

(66) 
In this work two different ways to compute the angle are presented. The first approach to compute the angle is to get this parameter using geometric relations provided by the model. The second approach consists in the minimization of the error, the idea is to find the angle that offers the best approximation.
Figures 33 and 34 show the and field in December of 2007.
Figure 33: xdisplacement field, 200712 
Figure 34: ydisplacement field, 200712 
Next the obtained results using each approach at each pendulum are presented.
Figure 35: Radial displacement at P1DR1 pendulum. 
Figure 36: Radial displacement at P1DR4 pendulum. 
Figure 37: Radial displacement at P2IR1 pendulum. 
Figure 38: Radial displacement at P2IR4 pendulum. 
Figure 39: Radial displacement at P5DR1 pendulum. 
Figure 40: Radial displacement at P5DR3 pendulum. 
Table 5 provides a summary of the angles at each pendulum according to each approach.
Id  Geometric  Optimized 
P1DR1  13  41 
P1DR4  16  38 
P2IR1  13  39 
P2IR4  16  36 
P5DR1  35  34 
P5DR3  39  41 
The obtained values using a geometric value of provides accurate results, concentrating the error at peaks.
In the second approach, the optimization process, the behaviour of dam is captured more accurately. The obtained values at the central pendulums keep the symmetry between them, and provide good results although the value is higher than geometric values. The obtained values at the exterior pendulums (close to the abutments) have similar value than geometric ones.
In Table 6 the error depending the used approach is shown. Two different types of errors have been performed, the first one, the mean absolute error which is computed as

(67) 
and the mean relative error which is computed as

(68) 
First approach provides an average error of while the second approach just a of error.
Geometric  Optimized  
Abs. Er. (mm)  Rel. Er. ()  Abs. Er. (mm)  Rel. Er. ()  
P1DR1  2.19  11.5  0.90  4.8 
P1DR4  1.21  8.7  0.65  4.6 
P2IR1  2.31  11.4  1.28  6.3 
P2IR4  1.37  9.3  0.96  6.6 
P5DR1  0.77  5.4  0.76  5.4 
P5DR3  0.39  4.8  0.37  4.6 
Average  1.38  8.5  0.82  5.4 
The provided results by as well as the performed transformation to radial displacements approximate accurately the real data provided by the pendulums. The obtained values in the second approach have coherence, keep the symmetry, and the numeric value is totally reasonable taking into account that pendulums can suffer distortions by different factors. The proposed numerical model captures the real behaviour of the dam accurately.
In this section a comparison between the obtained results using the with a tetrahedral mesh and the results obtained using COMET with hexahedrons of 20 nodes (provided by M.Sc. Fernando Salazar) was performed. The aim of this study is to prove that the use of tetrahedral elements can also lead to accurate solutions.
The boundary conditions applied in Software COMET are: a constant temperature value (provided by Bofang formulation at the middle point in the wall) at upstream dam face below water level, the air temperature + 2 at downstream and the reference temperature equal to .
The comparison was performed using the second approach for both cases; COMET and . Figures 41 and 42 provide the comparison through the radial displacement field at each pendulum.
(a) P1DR1  (b) P1DR4 
(c) P2IR1  (d) P2IR4 
Figure 41: Radial displacement comparison at central pendulums 
(a) P5DR1  (b) P5DR3 
Figure 42: Radial displacement comparison at exterior pendulums 
Table 7 shows the obtained results through the mean absolute error for both approaches.
Geometric  Optimized  
Dam App  COMET  Dam App  COMET  
P1DR1  2.19  1.19  0.90  0.96 
P1DR4  1.21  0.64  0.65  0.60 
P2IR1  2.31  1.44  1.28  1.28 
P2IR4  1.37  0.97  0.96  0.95 
P5DR1  0.77  1.56  0.76  0.62 
P5DR3  0.39  0.89  0.37  0.70 
Average  1.37  1.12  0.82  0.85 
Same comparison but this time through the mean relative error is shown in Table 8.
Geometric  Optimized  
Dam App  COMET  Dam App  COMET  
P1DR1  11.5  6.2  4.8  5.1 
P1DR4  8.7  4.6  4.6  4.3 
P2IR1  11.4  7.1  6.3  6.3 
P2IR4  9.3  6.6  6.6  6.5 
P5DR1  5.4  10.9  5.4  4.4 
P5DR3  4.8  10.9  4.6  8.6 
Average  8.5  7.7  5.4  5.9 
According with the obtained results (Table 7 and 8) the use of tetrahedral elements leads to accurate solutions. These results also prove the correct working of the new Application ().
This thesis presents a new application inside the KRATOS environment for elastic structural verification of dams.
The thermomechanical coupling was solved by introducing a thermal component into the mechanical constitutive law (one way coupling), and based on the obtained results, this modified equation is able to approach the physics of the real problem.
The developed graphical interface satisfies the basic user’s requirements. The application has a user friendly design and offers the possibility of selfcustomization through the use of input data tables.
Regarding the new implemented boundary conditions, the Bofang constrain provides an accurate temperature field in the upstream submerged part and its annual variation is in good agreement with the test field measurements. The limited range of temperatures given in the cases analysed, permits simplifying the application of the Bofang law. The uplift condition was implemented without taking into account the influence of drains.
On the other hand, with this thesis it was proved that the use of tetrahedral elements can be a valid option for these type of simulations. In dam engineering, the use of hexahedral elements in numerical analysis is very common due to its good response under bending loads. However, hexahedral elements become too rigid when it is required to deal with nonregular geometries, e.g. when taking into account the influence of a spillway. For a complex geometry, the use of tetrahedral elements presents a good option. Regarding the computational cost, the hexahedrons is expensive compared to the tetrahedrons.
As a final conclusion, the software Dam Application (DamApp) satisfies all the requested objectives, providing accurate results within a userfriendly interface.
As a future works, some improvements can be implemented in order to improve the robustness of this application:
As long term goal, it can be mentioned: the development of a module for simulating the solidification process during the construction phase, the introduction of nonlinear constitutive laws, the coupling between fluid structure and foundation (e.g. in a rock joint), all of them with the aim of analyzing extreme conditions.
This appendix is an overview of Finite Element Method (FEM) applied to the mechanical problem. Mainly, is based in the book of Prof. Eugenio Oñate, Structural Analysis with the Finite Element Method. Linear Statics [12] [13], where all these concepts are dealt extensively and in the book of Prof. Xavier Oliver and Prof. Carlos Agelet de Saracibar, Mecánica de medios continuos para ingenieros [14].
There are a wide number of structures of practical interest which can be analysed following the assumptions of 2D elasticity. These types of structure have a sort of prismatic geometry. Depending on the relative dimensions of the prism and the loading type, two categories can be distinguished;
Figure 43: Plate under plane stress 
Figure 44: Section under plane strain 
Displacements, strains and stresses fields
Plane stress and plane strain assumptions imply that transversal sections to the prismatic deform in the same way, and the displacement along z axis is also negligible. Applying this simplification only a generic 2D transverse section in the plane needs to be considered for the analysis. The displacement vector field of a point is

(69) 
Strains can be derived by the displacement field 69

(70) 
In plane strain case, the longitudinal strain () is assumed to be zero, but not for plane stress problem where is considered zero. The strain vector can be described as

(71) 
From the equations of the strain field 70 is deduced that and are zero. As it was explained, the longitudinal stress () does not contribute to the internal work. The stress vector can be described as

(72) 
According with 3D elasticity theory and applying the stated assumptions, the following matrix relationship can be obtained

(73) 
where C is the elastic material matrix or constitutive matrix

(74) 
Depending on the type of problem, the components of matrix take different expression. An usual way to define these coefficients is as function of Young Modulus and Poisson's ratio. According with MaxwellBetti theorem, is always symmetrical, so .
In the case of isotropic elasticity, the components of 74 are

(75) 
For an orthotropic material with principal orthotropy directions along the 1, 2, 3 axis (where 3 is outofplane direction), the matrix takes the following expression for plane stress and plane strain, respectively.

(76) 

(77) 
where


The symmetry of requires (MaxwellBetti Theorem)

(78) 
In the case that the solid is subjected to initial stress, the constitutive relationship must be modified. The total strain () can be divided as the sum of elastic () and initial () strains. The constitutive equation reads as

(79) 
Thermal strain contributions
One of the possible causes of initial strains 79 is the thermal effect. In the case of isotropic material, the expression of the vector can read as

(80) 
where is the thermal expansion coefficient and T is the temperature variation at each point. As it is shown in the thermal strain vector, the temperature variation does not generate shear strains.
Initial stresses
The solid can also be subjected to initial stress (). This initial stress can come from different sources. The total stress in the new equilibrium configuration is just the addition of initial stress to the Eq.79

(81) 
where can be decomposed as

(82) 
Virtual work expression
The Principle of Virtual Work (PVW) for 2D elasticity problem can be stated as [29]

(83) 
The integral in the l.h.s represents the work performed by the stresses (, , ) over the virtual strains (, , ). The terms in r.h.s represent the virtual work of the body forces (, ); the surface tractions (, ); and the external points loads (, ). is the area and is the boundary of the transverse section. For a plane stress problem is the thickness of the solid, and for plane strain problem is equal to one.
The PVW shows us the continuity requirements. Derivatives of the displacements are needed, so continuity is needed.
The shape function is a fundamental concept to understand the Finite Element Method (FEM). Shape functions () are defined at each node of a finite element, and allow us to obtain the value of a nodal variable at any point of the element through interpolation.
Before introducing the shape functions of the linear triangular elements and the fournoded Lagrangian elements, the concept of natural coordinates is presented.
The natural coordinates (, ) are normalized coordinates (range from 1 to 1). From Figure 45 can be deduced that

(84) 

Figure 45: Natural coordinate system for rectangular element 
where and are the coordinates of the element centroid. The differentials of area in Cartesian and natural systems are related by

(85) 
The integration of a function () over a rectangular element can be expressed in natural system as

(86) 
The shape functions can only reproduce exactly a polynomial solution of order equal or less than a polynomial contained in the shape functions. Pascal's triangle (Figure 46) shows the relation between the polynomial degree and the number of needed terms. It can be deduced that the higher is the order of that complete polynomial, the more accurate is the finite element solution.
Figure 46: Pascal's triangle in two dimensions 
A 2D complete polynomial of degree can be written as

(87) 
where the number of terms is

(88) 
The shape functions of triangles and tetrahedrons are formed by complete polynomials, but quadrilaterals and hexahedrons contain incomplete polynomial terms.
Shape functions must satisfy the same requirements in natural coordinates and in Cartesian coordinates. These conditions are:

(89) 

(90) 
Only two basic types of elements are presented: the triangular element of three nodes [30] and the rectangular Lagrangian element of four nodes.
The fournoded Lagrangian element (Figure 47) is the simplest element of the Lagrangian family and it coincides with the element developed by Argyris and Kelsey [31]. The way to obtain the shape functions is presented hereafter. Considering a generic node , the 1D shape functions in local directions (,) can be expressed as

(91) 
In this case, values of and take a value equal to 1 since is a linear element. To obtain the shape functions for a 4noded Lagrangian element only is necessary to multiply together, obtaining

(92) 
Figure 47: Fournoded Lagrangian element 
The shape functions of the linear triangular element in natural coordinates (Figure 48) are presented below

(93) 
Figure 48: Natural coordinates for a triangular linear element 
For a generic twodimensional element of nodes, the displacement field can be defined as

(94) 
where , are the horizontal and vertical displacements and is the shape function of node . In matrix form can be rewritten as

(95) 
or

(96) 
where u is the displacement vector of a point. The shape functions matrix of the element and the th node are defined as

(97) 
The nodal displacement vector of the element and the th node are defined as

(98) 
In the case of the strain field, it can be written as


(99) 

which in matrix form is written as

(100) 
or

(101) 
where

(102) 
is the element strain matrix and

(103) 
is the strain matrix of node .
To obtain the expression of the stresses field, it is just needed the substitution of Eq.101 into Eq.73, obtaining

(104) 
If initial strains and stresses are considered, from Eq.81 can be deduced that

(105) 
From the Principle of Virtual Works expression applied at the equilibrium of an element is possible to obtain the expressions of the stiffness matrix K and the force vector f.
Let us suppose that the uniformly distributed forces per unit area act over the body of the element (mass forces b), and uniformly distributed forces per unit length act over one of its sides (surface forces t). Moreover, supposing that the equilibrium of the element is achieved at the nodes, it is possible to define punctual forces acting at the nodes (nodal forces of equilibrium q) that must balance the forces that appear due to the element deformation and the rest of applied forces.
Taking into account these assumptions and using Eq.83, the Principle of Virtual Works applied to the element is written as

(106) 
In plane stress problems, is the real thickness of the structure, while in plane strain usually takes the value of 1.
From Eq.96 and Eq.101, it is possible to write

(107) 
Substituting these equations in Eq.106, the next expression is obtained

(108) 
Due to the arbitrariness of the virtual displacements , the expression can be written as

(109) 
Using the stress expression of Eq.105, the following expression is obtained

(110) 
after rearranging terms, the equation reads as

(111) 
it can also be expressed as

(112) 
where

(113) 
is the elastic stiffness matrix of the element, and

(114) 
the equivalent nodal force vector for the element, formed by initial strains, initial stresses, body forces and surface tractions. These vectors can be defined as

(115) 

(116) 

(117) 

(118) 
The global equilibrium equation of the mesh is obtained just imposing the sum of nodal forces of equilibrium at each node must be equal to the external forces

(119) 
where the left part of the expression represents the sum of contributions of the vectors of nodal forces of equilibrium of the different elements that are sharing the global node j, and represents the vector of external punctual forces acting on such node. Therefore, the global equilibrium equation of the mesh can be obtained by assembling the contributions of stiffness matrices and the equivalent nodal force vectors of each element. The global matrix equation can be written as

(120) 
where K, a and f are respectively the stiffness matrix, the vector of nodal displacements, and the vector of equivalent nodal forces of the whole mesh.
There are structures that due to geometrical, mechanical or loading aspects are not possible to compute using plane stress, plane strain or axisymmetric assumptions. The only alternative is performing a full three dimensional (3D) analysis using the general 3D elasticity theory.
Despite its apparent complexity, the analysis of a 3D solid with FEM does not introduce major conceptual issues. 3D elasticity theory is a straightforward extension of the 2D case. Some examples of solids with irregular shapes are shown in Figure 49.
] 
Figure 49: Structures which require a 3D analysis: a) Double arch dam including foundation effects. b) Pressure vessel. Imagen taken from [12] 
The structure of this section is the same as previous one.
Displacements, strains and stresses fields
The displacement vector field of a point is defined by three components

(121) 
where , , and are the displacements of a point in the direction of cartesian axes (, , ), respectively.
The strain field is defined by the six standard components of 3D elasticity [32]. The strain vector is written as

(122) 
The normal strain () and the tangential strains (, ) are not equal to 0 anymore, and can be computed as

(123) 
The stress field is defined by the six stress components which are the conjugate of the six nonzero strains Eq.123. The stress vector is

(124) 
According with 3D elasticity theory (Eq.73) and assuming an isotropic material, the constitutive matrix reads as

(125) 
In the case of isotropic materials only are required two parameters: the Young modulus () and the Poisson's ratio (). In the general case, anisotropic elasticity, the constitutive matrix is composed by 21 independent parameters.
Initial strains and stresses
Taking into account the contribution of initial strains and stresses, it is possible to write

(126) 
Thermal strain contribution
The initial strain vector due to thermal effects can be defined as

(127) 
Virtual work principle
The PVW expression for solids in three dimensions is

(128) 
Eq.128 is just an extension of 2D solids Eq.83. It is worth to remark that is written in vectorial form.
Figure 50: Natural coordinates system , , in a tetrahedron 
For a tetrahedron (Figure 50) with right edges a,b,c is defined

(129) 

The differential of volume can be expressed as

(130) 
The integral of a function () over the element is an extension of Eq.86 and can be written as

(131) 
The shape functions for a linear tetrahedron can be expressed in terms of natural coordinates (,,) as

(132) 
The discretization of the displacement field in a 3D case follows the same criterion than in 2D. The only change is the addition of a new component. The displacements can be written in matrix form as

(133) 
In the case of strains, according to Eq.101, it is just necessary to adapt B to new requirements

(134) 
To obtain the expression of the stress field, it is only necessary to substitute Eq.101 into Eq.73, and consider the initials strains and/or stresses if necessary.

(135) 
From the Principle of Virtual Works expression applied at the equilibrium of an element is possible to obtain the expressions of the stiffness matrix K and the force vector f.
Taking into account the same assumptions than in 2D case and using the Eq.128, the Principle of Virtual Works applied of an element can be written as

(136) 
The virtual displacements and the virtual strains are interpolated in terms of the virtual displacement values in the standard form

(137) 
Substituting the Eq.15 into Eq.136 and simplifying the virtual displacements, the following expression is obtained

(138) 
Using the stress expression of Eq.105, the following expression is obtained

(139) 
it can also be expressed as

(140) 
where

(141) 
is the elastic stiffness matrix of the element, and

(142) 
the equivalent nodal force vector for the element, formed by initial strains, initial stresses, body forces and surface tractions. These vectors can be defined as

(143) 

(144) 

(145) 

(146) 
Finally the global matrix equation can be written as

(147) 
where K, a and f are respectively, the stiffness matrix, the vector of nodal displacements, and the vector of equivalent nodal forces of the whole mesh.
At this point, the discretization of equilibrium equations have been already obtained, but a convenient numerical technique for integrating the expressions of the stiffness matrix and the vectors of equivalent nodal forces of the elements has not been presented.
The integration of the expressions of the stiffness matrix and the vectors of equivalent nodal forces is performed by means of GaussLegendre quadratures. This technique allows us to integrate any function over a normalized domain using the tabulated Gauss points (coordinates and weights). However, it is necessary to transform the integrals over the element domain into integrals over the normalized space.
To carry out this transformation, it is necessary to introduce the concept of the isoparametric interpolation. The concept of isoparametric interpolation means that; the displacement shape functions are used to interpolate the element geometry in terms of the nodal coordinates. In this section, the isoparametric formulation for 2D and 3D solids is presented.
In a 2D case, the coordinates of a point within element are expressed in isoparametric form as

(148) 
where are the element shape functions. These equations relate the Cartesian coordinates of a point with the natural coordinates, in this case and . To formulate the change of variable from Cartesian to natural coordinates is necessary to compute the determinant of the Jacobian matrix.
To obtain the Jacobian matrix, it is necessary to derive the shape functions respect to and . Considering and apply the chain rule.

(149) 
or in matrix form

(150) 
where is the Jacobian matrix of the transformation of the derivatives of in the natural global axes. The Eq.150 can also be expressed as

(151) 
where is the determinant of the Jacobian matrix ("the Jacobian"). The determinant relates the differential of area between two coordinates systems

(152) 
Using an isoparametric transformation, the Jacobian can be obtained in terms of

(153) 
and

(154) 
Subtituting Eq.151 into Eq.103 the strain matrix in terms of natural coordinates is obtained

(155) 
where

(156) 
The components of J are detailed in the .
Once the strain matrix in natural coordinates is presented, the element stiffness matrix can also be computed in the normalized natural coordinate space. For a quadrilateral domain (Figure 45) it reads as

(157) 
in the case of triangular domain (Figure 48) it can be stated as

(158) 
Same procedure must be followed for computing the load vector.
Up to here the way to write the stiffness matrix and the force load vector in the natural coordinate space have been presented. To carry out these computations it is necessary to integrate the corresponding functions. Below, the numerical integration by a Gauss quadrature is presented.
The integral () over the normalized isoparametric quadrilateral domain can be evaluated using 2D Gauss quadrature

(159) 
where and are the number of integration points along each natural coordinate (, ), respectively; and are the natural coordinates of the integration point and the corresponding weights are defined by and . The exact evaluation of a fournoded rectangular domain requires a 2x2 quadrature.
A three noded triangular element only requires one integration point [29] [33] and can be written as

(160) 
So to obtain stiffness matrix in a quadrilateral domain, it is just necessary the substitution of Eq.159 into Eq. 157, obtaining

(161) 
Applying the same criterion for triangular domain, Eq.160 into Eq.158, the following expression is obtained

(162) 
To compute the stiffness matrix is necessary to evaluate the Jacobian , the deformation matrix , the constitutive matrix , and the thickness at each integration point.
To compute any of the vectors of equivalent nodal forces same procedure must be followed. For a quadrilateral domain the body force vector is

(163) 
and for a triangular domain

(164) 
To finish this section it is worth to remark that the computation of the vector of surfaces force is different. This difference comes from the integral since it is performed over the contour element. In case of 2D solids, this contour usually represents a straight line. These computations and extensive explanations can be found in [12].
This section is devoted to 3D solids, a generalization of 2D solids. The formulation of a tetrahedral element is presented below. In case of having interest in other type of elements as hexahedrons, the following book of Prof. Eugenio Oñate can be consulted [12].
The coordinates of a point within element are expressed in isoparametric form as

(165) 
The derivatives of the shape functions in Cartesian coordinates are computed using the chain rule
