Simulation of concrete ageing on dams as illustrated by numerical analysis of Jablanica HPP

The paper presents a model for numerical simulation of interaction between the fluid and the reinforced-concrete structure under seismic load for spatial 3D problems, with an emphasis of concrete dams. The model is based on the so called partition scheme for coupled (multi-field) problems. The model especially focuses on the concrete ageing effects, which has proven to be highly significant in the analysis of older structures. The paper starts with a brief description of the numerical model, which is followed by a more detailed description of the concrete ageing model, and by analysis of Jablanica hydropower plant (HPP) on which the validity of the model and related software was checked.


Introduction
Water engineering structures called dams are built over river valleys or riverbeds to make use of the water mass [1]. By creating an artificial water-storage reservoir or lake, or by temporary retention of water, attempts are made to fulfil the basic purposes of dams: regulation of water flow, water supply, irrigation, generation of electricity, navigation, recreation, etc. Arch gravity dams [1], such as the dam of the Jablanica hydropower plant, are dams with the base ratio of (0.33 -0.65) H, where H stands for dam height. Due to considerable thickness of these dams, a greater part of the load is transferred onto the bottom of the valley, while a smaller part is transferred to abutments [1]. During design and construction of first concrete dams, a relatively small attention was paid to the material itself and to its behaviour over time, although even then the engineers were aware of the downsides of such approach. Papers published by Abrams [2] in 1918 about concrete mix design formulas and water-cement ratios, and advances in 1920s leading to better understanding of the cement hydration process, have resulted in an increased care about the material and its behaviour, namely in response to formation of unwanted cracks that expedite and accelerate dam degradation processes. The design and construction of Hoover dam (then the largest dam in the world) in the 1930s, with enormous quantities of concrete, has led to considerable advances in the concrete mix design, transport, placement, and cooling of placed concrete. The influence of water-storage reservoir must be taken into account in dam analysis. In the beginning of analyse, the influence of storage reservoir was considered as hydrostatic load only. However, dams were then massive earthfill structures, with mass often by far exceeding structural requirements. The period of second industrial revolution, marked by an extensive dam construction spree, was also the period of great scientific discoveries. The first strictly scientific analysis of the dam-storage system behaviour was made by Harold Malcolm Westergaard [3] in1933. In his now already legendary paper, he assumed an infinitely long reservoir and a non-compressible liquid, and proposed an analytic expression for the distribution of hydrodynamic pressures exerted on the rigid dam, with regard to horizontal harmonic oscillations. This expression, with some minor modifications, has been used to this day. The study of fluid-structure interaction (FSI) received a special impetus in the 1960s and 1970s, when many researchers focused their studies on this issue. Over the past three decades, the FSI research has been developing in two directions: the one with a greater emphasis on description of the fluid through development of meshfree methods, where the SPH (Smoothed Particle Hydrodynamics) being the best known, and the other with a greater emphasis on description of the structure. Both of these research alternatives have resulted in a very good description of the problems causing great distortion of grid during simulation, which is very favourable for the fluid flow problems. In addition to the study of fluids, very good results have been achieved in simulation of various problems of plastic flow of materials in continuum mechanics. However, some accompanying problems have also occurred, the greatest of them being the need to use powerful computers and a very long time needed to perform calculations. For tackling real problems, such as those involving analysis of gravity dams with water storage reservoirs, engineers use "older methods" based on Euler's description of the motion of fluids, where fluid is considered based on the volume it occupies, without considering behaviour of individual particles. There are two reasons behind the use of these methods. The first one is that we are generally not interested in fluid behaviour during dam analysis, except for the effects of fluids on the structure, and these can be calculated with sufficient accuracy using much simpler methods. The second reason is that a very dense particle grid would be needed to describe local effects of fluid, such as sloshing, waves, etc., which would greatly increase the cost and time of calculation, while total effects exerted on the dam would be small. This is why models based on Euler's approach are still being intensively developed. A number of papers have been published in this area over the past decade: a numerical model of boundary and finite elements has been developed for the analysis of dam-water storage systems [4], the so called endurance time (ET) method for seismic analysis of concrete gravity dams is being intensively studied [5], seismic behaviour of arch dams is studied with regard to water storage topology [6], etc.

Brief presentation of numerical model
As already emphasized, numerical model for real simulation of the structure in direct contact with the fluid must include simulation of interaction between the fluid and the structure, if their real response is to be obtained. This problem is especially pronounced in dynamic/seismic loading, and can usually be found in literature under the name of "coupled problem" or "multi field" problem. The coupled problem involves two or more fields that either touch or pervade one another, such as the dam with a water storage reservoir. The fields are usually time dependent. The state of one field is continuous bind with the state of another field, and so none of them should be considered separately. The link is established through the equations of state, which describe specific physical phenomena. There are two basic approaches to the resolution (simulation) of coupled problems, i.e. to the solution of problems involving the fluid-structure interaction: -Monolithic approach: requires development of a mathematical model for each physical problem separately and, at that, equations describing the state of pressures in fluid and structure displacements are solved using a unique model/solver; -Partitioned approach: enables modularity of the model and, at that, equations describing the state of pressures in fluid and structure displacements are solved separately, and are linked via interaction limit. GRAĐEVINAR 71 (2019) 9, 749-767 Simulation of concrete ageing on dams as illustrated by numerical analysis of Jablanica HPP The model and calculation software developed for numerical modelling of the fluid-structure interaction problem is briefly described below. It is based on the partitioned approach, with the interaction force transfer in each problem-solving phase. It is favourable for problems in which the fluid has limited displacements, such as those of the dam and water storage behaviour under seismic load. The developed model and software is based on the finite element method (FEM) for space discretisation, and on the finite difference method (FDM) for time discretisation of the system [7][8][9][10][11]. The displacement formulation was used for the structure and soil, while the displacement potential formulation was used for the fluid [9][10][11]. The nonlinear structural analysis model used is presented in [12][13][14], while the eigen problem calculation model is presented in [10].

Numerical models
Fluid is a substance (liquid or gas) that is continuously deformed under the influence of external forces. Fluid can be viscous with friction between particles, or inviscid without friction. Each real fluid has some amount of viscosity, but this viscosity can often be neglected. In this paper, the fluid (water) is considered to be inviscid, and the so called Euler formulation is used [9][10][11]. If the structure and fluid are discretised by the finite element network and if appropriate basic functions are selected for displacement potential in the fluid N Ψ and for structural displacements N u , using the standard Galerkin formulation: (1) where Ψ is the potential of displacement of nodes of the finite element network of the fluid, u are displacements of nodes of the finite element network of the structure, then the differential equation of dynamic equilibrium of the system is obtained: (2) where: M f is the fluid mass matrix; C f is the radiation damping matrix; K f is the fluid stiffness matrix; f f is the vector of external node forces; Q t is the matrix of interaction of the fluid-structure system; ψ is the displacement potential vector; r f is the fluid density; u is the relative vector of displacement of movable limit as related to the base, and d is the base displacement vector. In explicit form, matrices frm equation (1) can be expressed as follows: where: V f is the domain (volume) of the fluid; W r is the surface/ limit of radiation; W f is the surface/boundary of free face of the fluid; W r is the surface/bundary of interaction, i.e. the contact surface between the fluid and the structure. When solving the fluid-structure interaction problem, nonlinearities are usually related to behaviour of the structure, and the fluid is considered to be linear, which implies unlimited negative pressures in the fluid. In reality, air bubbles are formed if the absolute pressure in an area falls below the fluid vapour pressure. This phenomenon is known as cavitation, and it can cause considerable damage to the structure. A simple fluid behaviour model based on the bilinear relationship between the pressure and mass dilatation is presented in Figure 1 [9,10]. Therefore, cavitation occurs if the following expression is satisfied: where: s -he mass dilatation p h -the hydrostatic pressure p a -the atmospheric pressure p v -the evaporation pressure and c is the speed of sound in the fluid.

Numerical model for structure
By spatial discretisation of the structure and by using the FEM, the equation of dynamic equilibrium with unknown node displacements u can be written in the known matrix form [7][8][9][10]: where: M s is the structural mass matrix; C s is the structural damping matrix; K s is the structural stiffness matrix; f s is the vector of external node forces. These matrices and vectors are defined by the expressions: In equations (6) N i are the base (shape) functions, and B is the displacement strain relationship matrix [8,14].
For real structures, the displacement strain relationship matrixis generally nonlinear, i.e.: which represents the so called geometric nonlinearity. In fact, matrix B depends on system displacements, i.e. it is not linear.
The relationship e-u is also known as the geometry model. The stress strain relationship s-e is also generally nonlinear and represents the so called material nonlinearity. The relationship s-e can also be written as follows: where D is the matrix of the stress-strain relationship and, in case of an elastic material, it represents a well known elastic constant matrix [8,14]. The relationship s-e is known as the constitutive law or material model. Material behaviour simulation, especially simulation of elasticbrittle materials such as concrete, in the state of full spatial stress, has not as yet been sufficiently mathematically defined. The situation is further complicated if the model is extended by adding additional effects of material behaviour such as: shrinkage, temperature deformations, creep, ageing, etc.

Numerical 3D model of concrete behaviour
A simple model, based on a small number of easily measurable physical parameters, previously developed at the Faculty of Civil Engineering, Architecture and Geodesy of the University of Split, is presented in this paper [12][13][14]. Poisson ratio, and it can be represented as a multi-surface surface ( Figure 2). The implementation of multi-surface flow surface, for dominant compression, requires separate calculation of stiffness matrix for each sextant but, on the other hand, it enables rapid convergence of the mathematical procedure.
In case of dominant compression, a new flow surface is calculated after liquid limit is exceeded, and the material matrix D is corrected. The procedure is repeated with new material matrix until the convergence is achieved. The traditional Rankine model is used for modelling concrete in tension. In tension, each direction of tension is considered separately. Concrete remains homogeneous and linearly elastic until the tensile strength of f t ' is achieved. The first crack occurs perpendicular to the direction of main stress, and the stiffness reduces, at the moment in which the tensile strength is exceeded. The change in stiffness of uncracked concrete between cracks (Gauss points) is simulated by gradual reduction of the tensile stress component perpendicular to crack surface ( Figure 3).

Figure 3. Concrete behaviour in tension
Even after occurrence of cracks, concrete is considered as a continuum with reduced stiffness. The fixed orthogonal smeared cracks model is used for the occurrence and development GRAĐEVINAR 71 (2019) 9, 749-767 Simulation of concrete ageing on dams as illustrated by numerical analysis of Jablanica HPP of cracks. According to this model, it is assumed that -after occurrence of the first crack -the crack position and direction do not change with the change of load. New cracks may occur only perpendicular to their direction. After occurrence of cracks, concrete becomes anisotropic, and crack directions define main axes of anisotropy. Cracks can be opened and closed during dynamic/variable application or relaxation of load, which can also be monitored by the model. Possible crack occurrence and development patterns are presented in Figure 4.

Numerical interaction model
Numerical model for the fluid-structure interaction problem is based on the so called partitioned solution approach. The term "structure" implies the structure itself and the surrounding soil. This model can be described with two second-order differential equations [9][10][11]. If displacement potential is selected as unknown for the fluid, and if displacements in nodes are selected as unknowns for the structure, then these two equilibrium equations can be written as follows: where (10) of the matrix Q, which is the matrix of interaction between the fluid and the structure ( Figure 5), and is obtained by integration via interaction surface in the direction normal to the surface ( n  ). According to the partitioned approach, each field is solved separately in each increment of the given load and in each iteration, including interaction forces (11) at the fluid and structure contact surface. The corresponding flow chart is given in Figure 6. The standard predictor-corrector procedure is used in this paper for finding solutions. (corrosion) of concrete. As a rule, this phenomenon occurs on the upstream side of the dam, while it is not so pronounced on the downstream side, at least not on the part that is not in contact with water.
Three basic types of concrete degradation in structures can generally be differentiated [15]: sulphate action (sulphate corrosion) which results in chemical and physical degradation of concrete, and is caused by water abounding in sulphates; alkali-aggregate (alkali-silica) reaction which also results in chemical and physical degradation of concrete, and is caused by certain types of aggregate abounding in silicium. In combination with some substances from cement gel, the so called swelling silica gel is formed, which physically destroys the concrete.
degradation by freezing and thawing which is generally physical in nature, but it also alters material properties of concrete.
The fourth degradation mechanism is corrosion of reinforcing (structural) steel. This mechanism primarily affects the structures where reinforcement is essential for mechanical resistance and stability, but the structures themselves are not protected with a sufficient layer of concrete. In any case, when other degradation mechanisms initiate damage to protective layer, then reinforcement corrosion can pose a serious threat to mechanical resistance and stability of the structure.

Change of material properties of concrete and dam safety
One of greatest challenges faced by engineers involves estimation of dam safety and its service life. This is especially difficult for the existing dams already affected by degradation of material, cracking and similar adverse phenomena, which point to the change in their bearing capacity. If the degradation process or the occurrence of cracks have been observed, then the probability of failure gradually increases and reaches the point in which the risk is no longer acceptable. If there is no degradation, the risk generally remains the same, i.e. the design risk is applicable. The risk can even slightly be reduced over time due to positive effects of concrete ageing (increase in compressive strength and elastic modulus). Of course, this statement is valid only for dams that are not exposed to some extraordinary events/effects, such as an earthquake action. A good indication of real condition of dams can be obtained by their proper monitoring, especially with regard to changes in material properties of such structures. In addition, numerical models with optimally created finite element network, with good simulation of material properties, and with accurately set load parameters, may constitute a powerful tool for predicting future behaviour and safety level of dams.

Results obtained by testing changes to concrete properties at US dams
The developed model, incorporated in a computer software, is based on investigations performed on the US dams by the US Department of the Interior, Bureau of Reclamation, Technical Service Center, Denver, Colorado, USA. These results are presented in [2,16]. Although investigations presented in the above literature are more based on finding chemical composition of concrete, and not so much on establishing physicomechanical properties, they still ofer a good indication of change in concrete behaviour over time. These investigations were used in this paper as the basis for predicting time-related changes in compressive strength and elastic modulus.

Development of compressive strength and elastic modulus of concrete over time
The compressive strength and elastic modulus of the concrete not degraded by corrosion, were determined experimentally in the study [3] for a number of dams ( Figure 7). The procedure was conducted by in situ sampling, and by subsequent laboratory testing of these samples. Blue markers and blue line represent the trend line of compressive strength results, while rose markers and rose lines represent the trend line of elastic modulus results.
The results obtained during this testing campaign agree quite well for all dams, except for the anomaly marked with a red circle, Figure 7. This anomaly was established at the East Park Dam that was built in 1910. The concrete originating from that period had a much greater water/cement ratio and, as expected, its strength is much lower compared to more recent dams. Simulation of concrete ageing on dams as illustrated by numerical analysis of Jablanica HPP 1938, also suffered considerable degradation of concrete due to corrosion [19]. Perhaps the most extreme example is the American Falls Dam, built in 1927, which was degraded to such an extent that it had to be replaced in 1977.
The aggregate from the Colorado River, used in Hoover Dam construction, is essentially non-reactive, while the aggregate from the Bill Williams River, used in Parker Dam construction, is highly reactive. Comparison of test results clearly reveals that the concrete used in Parker Dam does not show properties similar to concrete used in Hoover Dam, although the same concrete mix and the same cement class IV was used for both dams. Concrete used at Parker Dam achieved only about 60% of the compressive strength and elastic modulus compared to concrete used at Hoover Dam, and the former concrete also exhibits a much smaller increase in strength over time. Development of compressive strength for concrete without corrosion and concrete affected by corrosion is shown in Figure  9. It is interesting to note high compressive strength results that were obtained at Parker Dam 42 years (15330 days) after its construction, which can be attributed to samples taken from bottom parts of the dam (data in rose circle).
Four different types of cement were randomly used during construction of Parker Dam and, at that, one type of cement exhibited low alkalinity. High compressive strength results (ranging from 47.0 to 57.0 MPa) were probably obtained from samples taken from the zone in which this particular type of cement was placed. These results actually reveal the strength that the dam would have had, had this particular type of cement been used for the entire body of the dam. Some other samples exhibit much smaller strength values, which is probably due to the use of other types of cement. The relationship between compressive strength and elastic modulus of samples affected and unaffected by concrete corrosion is shown in Figure 10 for all tested samples. Although the coefficient of correlation is relatively weak, the trend curve clearly points to the difference between these two types of concrete. Individual comparisons (trend curves) of compressive strengths and elastic moduli are much more accurate when specific dams are considered (dams built with the same type of aggregate, dams with one or several types of cement, etc.). Trend curves can be a good indicator about whether a particular dam is affected by corrosion, how deep is the corrosion, and what are the consequences for dam safety.   Figure 11 for concrete without corrosion and concrete with corrosion.

Proposal of curves showing trend of change of material properties of concrete over time
As can be seen in relevant studies [15][16][17][18][19] concrete properties change over time: compressive strength, tensile strength, and elastic modulus. These parameters obviously depend on the composition of concrete (concrete mix) and on chemical composition of water acting on the dam. The above facts must be taken into account in the analysis of dams, older dams in particular.
Based on the presented investigations, a trend of change of these properties over time is proposed.

Trend of change of mechanical properties of concrete over time
Many authors have been studying the trend of change of mechanical properties of concrete over time. Their papers mostly focus on shrinkage and creep, as can be seen for instance in [20,21].
In this paper, the deformation of concrete due to ageing is considered indirectly, i.e. the initial elastic modulus and concrete strength increase over time. In fact, current mechanical properties of concrete and the corresponding material "changed over time" (either strengthened or weakened) are taken into account when choosing the s−e m relationship over a given period of time. The incorporation of concrete ageing for the uniaxial elastoplastic relationship s−e m is schematically presented in Figure 12, where t is the current time, t 0 is the initial time, f cm is the initial design compressive strength of concrete (usually at 28 days), f ctm is the initial design tensile strength of concrete, E cm is the unique initial elastic modulus of concrete in compression and tension, e c,0 is the initial design crushing deformation of concrete in compression, and e t,0 is the initial design crushing deformation of concrete in tension. The increase in strength as a result of concrete ageing is considered in many relevant regulations and, in that respect, [23] the following expression is proposed by CEB: In the above expression and in Figure 12 Simulation of concrete ageing on dams as illustrated by numerical analysis of Jablanica HPP It can be seen that the modified curve (2) follows almost ideally the experiment. An increase in time exponent from -0.5 to -0.2 enables faster increase of curve over time, and the empirical coefficient 1.28 enables better adjustment to real results. proposed in [11] Expression (2) is proposed for elastic modulus of concrete in time t, according to EC-2 [24]: uz (14) where E cm,t denotes the mean value of the so called secant modulus of elasticity (GPa), while f ck,t is the characteristic compressive strength of a cylinder (MPa) in time t. Experimental results from Figure 8, describing change of elastic modulus of concrete over time were also read, converted into MPa, and plotted as shown in Figure 15 (blue curve). If presented results of compressive strength over time are used in the calculation of elastic modulus via equation (14), it can be seen that elastic modulus is also not well approximated by the proposed expression (green curve). That is why the curve of change of elastic modulus over time was also modified; it is proposed in the following form: , All these curves are presented in Figure 14. It can be seen that the proposed curve (4) is too low in the first days, but it increases over time and gets very close to the experiment. An increase in time exponent from -0.5 to -0.3 ) enables faster increase of curve over time, and the empirical coefficient 1.14 enables better adjustment to real results. The coefficient "s" value was adopted as 0.25 in the calculation of compressive strength and elastic modulus, as proposed in literature [23].

Trend of change of mechanical properties of degraded corrosion-affected concrete over time
Compressive strength results for concrete not affected by corrosion and concrete affected by corrosion are shown in Figure  9. Experimental results from that figure were read, converted into MPa and plotted as shown in Figure 15. Experimental results for nondegraded concrete and corrosion-degraded concrete are compared in this figure.
The curve of degraded concrete has the form similar to that of nondegraded concrete. This curve has the following form where coefficient "s" amounts to 0.55.
It can be seen that this curve, although initially "slightly displaced", later on follows the experiment quite closely. As in curve presented in expression (2), an increase in time exponent from -0.5 to -0.17 enables faster increase of curve over time, and the empirical coefficient 1.14 enables better adjustment to real results. Experimentally determined change in tensile strength of concrete over time is shown in Figure 11. As in the case of compressive strength, experimental results from that figure were read, converted into MPa and plotted as shown in Figure  16, where experimental compressive strength results for undegraded concrete and corrosion-degraded concrete are compared. Results obtained by splitting test (ST) were used as reference tensile strength. Expression (17) is proposed for the mean tensile strength of concrete, according to [24]: Straight lines were chosen for interpolation curves as the trend of change in tensile strength is almost linear. The following expression is proposed for interpolation of nondegraded concrete: while the following expression is proposed for interpolation of degraded concrete (19) where t is the time in days and a is the coefficient of concrete degradation, which was adopted as 0.6 in the above expression, as approximately proposed in [17].

Figure 16. Comparison of experimental tensile strength results for nondegraded concrete and corrosion-degraded concrete [17] and curves proposed in [11]
The expressions deduced have been incorporated in the developed program, and are used as basis for determining compressive strength, tensile strength and elastic modulus, for a predefined number of days.

Example: Jablanica Dam
The dam of the Jablanica Hydropower Plant was analysed to illustrate and check the developed model. This arch gravity dam was built in 1954 and has been used ever since for seasonal regulation of the Neretva River. The dam is composed of fourteen blocks interconnected with wide contact areas (measuring approximately 1 m), which were realised for dissipation of concrete temperature. The data about the dam and water storage reservoir of the Jablanica Hydropower Plant ( Figure 17) were obtained from [25].
The construction height of the dam is 85,0 m, the dam length at crest level is 210,0 m, and the dam volume is 130.000,0 m 3 . The installed flow is 180,0 m 3 /s, and the maximum flow of all evacuation spillways is 2915,0 m 3 /s. The installed capacity is 150,00 MW, and the mean annual production is 792,0 GWh.

Description of soil properties and adopted parameters
The length of Jablanica water storage reservoir is 27 km in straight line. The greatest depth is about seventy meters near the dam of the Jablanica Hydropower Plant. The water basin takes up an area of approximately 14,3 km 2 , and the volume of the water storage reservoir is approximately 290.000.000 m 3 . The lake is longitudinally shaped, with the dominantly east-west orientation. It has numerous coves, namely in the central and bottom parts of the lake. An average width of the water storage reservoir varies between 250 and 400 meters [1][2][3].
The Jablanica Lake basin is situated at 270 m above sea level. The biggest watercourse that empties into the water storage reservoir is the Neretva River, with the mean annual flow rate of 53,0 m 3 /s. As to landscape and scenery around the lake, the area abounds in deciduous forests.

Seismic features
According to former classifications, the Jablanica Lake basin belongs to zone 6°MCS, and the epicentral area is situated in Konjic with earthquakes of the order of 4°MCS. There are two epicentral zones to the south-east along the Neretva River valley, with earthquakes of the order of 6°MCS to 7°MCS [1][2][3].
Local seismic activity involves all earthquakes in the vicinity of the Jablanica HPP (43.692° N17.732° E) within the perimeter of 75 km. As many as 5,472 earthquakes were registered in this area in the period from 1386 to 2011. Considering the significance of seismic activity in this zone, a special attention was paid to the earthquakes within the perimeter of 25 km from the Jablanica HPP dam. 430 earthquakes were registered in this perimeter in the period from 1386 to 2011. GRAĐEVINAR 71 (2019) 9, 749-767 Simulation of concrete ageing on dams as illustrated by numerical analysis of Jablanica HPP Maximum horizontal accelerations and earthquake intensities at the dam site are presented in Table 1.

Geological properties of soil
A respectable array of geological investigations and soil testing campaigns was conducted in the surrounding area and in close proximity of the dam. General geological composition of soil in dam area is dominated by a gabbro massif, which is intersected with various dike equivalents, while joints and small cavities are intersected with various mineral parageneses [29]. Table 2 shows soil and rock massif parameters that were used in the simulation of the fluid-structure interaction for the Jablanica dam [4]. Rock massif parameters were determined by individual zones via statistic processing of data obtained following geophysical survey as conducted by IGH Zagreb.

Description of properties and structural parameters
According to [29] the mean value of the static modulus of elasticity for concrete amounts to: E cm = 31000 MPa. In addition, according to [29], it is recommended to use in dynamic analysis the dynamic modulus of concrete that can be obtained from the following empirical equation: (20) The rectified elastic moduli were calculated according to (4), and the following values were obtained for the 60-year old concrete: -static modulus of elasticity of E cm = 42 400,0 MPa, -dynamic modulus of elasticity of E cm,D = 58 000,0 MPa.
Both results corresponds quite well with measured results [29]. According to [29], the mean compressive strength of concrete that is considered in calculations amounts to f cm = 30,0 MPa. The rectified compressive strength of concrete was calculated according to (2), and the compressive strength of f cm = 41 MPa was obtained for the 60-year old concrete. According to [29], the mean tensile strength of concrete incorporated in the dam amounts to: f cm = 3,0 MPa. The rectified tensile strength of concrete, calculated according to expression (7), was used in this paper. The tensile strength amounting to f ctm = 4,3 MPa was obtained for the 60-year old concrete.  The Poisson ratio of concrete varies between 0.15 and 0.22. That is why, according to [29], its mean value obtained by laboratory testing is used for static calculations: n c =0, 19.
In dynamic analyses according to [29] the value of this coefficient must be increased by 30%, i.e. to the value of: n c =0,25. Table 3 shows concrete parameters that were used in the simulation of the fluid-structure interaction for the Jablanica dam. Table 4 shows water parameters that were used in the simulation of the fluid-structure interaction for the Jablanica dam.

Model
A finite element network model composed of elements of the dam itself, elements of the surrounding soil, and elements of the water (upstream and downstream storage), was prepared for the study of dam behaviour in the scope of static and dynamic analyses. The finite element network of the overall system is shown in Figure 18, and the finite element network for the dam body is shown is Figure 19. The dam, fluid, and soil were discretised with 240, 2300, and 6920 global finite elements, respectively, which were further divided, if needed, into smaller elements during the analysis. It should be emphasized that temperature effects, which can greatly change stress values and behaviour of the structure, were not considered in the analyses.

Static analysis
In static analysis, the dam is loaded by its own weight and by hydrostatic pressure. The influence of temperature on the dam was not considered in this paper. In the numerical model, the dam water level was set to the spillway level. The results from this analysis were compared with the dam monitoring results, i.e. with the results from the "Technical monitoring report, prepared in 2014, for the dam, structures, and soil in the dam and water storage areas" [25].
In the scope of this monitoring, dam displacement is measured using the inverse plumb lines method. Dam displacements are monitored, using this method, in three measurement sections: IV, VII and X, at three positions for each section. Points along the dam body in which automatic coordimeters were placed are presented in Figure 20. Dam displacements for own-weight load and hydrostatic load, increased by 2000 times, are shown in Figure 21. Dam geometry is shown in light blue colour, and displacement by hydrostatic load at maximum water storage level is presented in dark blue colour. Dam displacements in the central segment (Segment VII) are presented in Figure 22. Dam displacements obtained by numerical model (11.49 mm in the top part of the central segment) correspond very well with displacement results presented in the "Technical monitoring report prepared in 2014 for the dam, structures, and soil in the dam and water storage areas". 11.28 mm [25]. Numerical analysis conducted with initial values of elastic modulus, compressive strength, and tensile strength, gives somewhat greater displacements (15.28 mm at the top part of the central segment), which is not in accordance with measured values. This simple example clearly shows that the influence of ageing, i.e., concrete hardening, is a significant factor in the analysis of older structures.

Analysis of free oscillations of the dam
Free oscillation modes of the dam, i.e. the first five of its own oscillation modes, are presented in Figure 23. In Table 5.the frequencies of the first five forms of dam oscillation are given obtained in this work with those presented in [29]. Table 5 presents comparison of dam frequencies obtained in this paper with those presented in [29]. Dam is shown in blue, and displacement, i.e. modal form, is shown in red. Detailed comparison of results from this paper and the results from [29] (analysed using the DIANA software) shows that the frequencies are somewhat different, which is quite expected considering the size and complexity of the model, and the size of the surrounding terrain and water storage reservoir covered by the model. It is interesting to note that own oscillation modes 1 and 2 are quite similar in both papers, with good correspondence of frequencies. At the same time, the frequency and mode of oscillation of the fourth mode from this paper, 7.895 Hz, is close to the form and frequency of the third mode from [29], which is 8.188 Hz. It can generally be concluded that the results correspond very well with the results given in [29]. The differences can be attributed to the density of the network and to the way in which own values were calculated.

Note on structural damping
The Rayleigh viscous damping is regularly used in the method of direct integration of equations of motion, and the latter is also used in this study. In this viscous damping, the damping matrix C is expressed as a linear combination of matrices M and K.
The damping matrix C is introduced in the context of the use of linear elastic analyses and so, when applying nonlinear dynamic analyses with material models that include internal dissipation of energy, the damping matrix C should cover only that part of the dissipation of energy that is not comprised in material models. An open question still remains: to what extent should the Rayleigh viscous damping be included when nonlinear models of material and geometry are introduced. To take into account this fact, the viscous damping of 0.5% is used in the paper. This damping percent simulates "other nonlinearities" because the nonlinear behaviour of material is included through the previously presented model. Otherwise, as can be seen from the results, greater nonlinearities have not been registered in the dam even in cases of strong earthquake activity.

Dynamic analysis of Jablanica HPP dam
In order to determine response to dynamic excitation (earthquake), the dam was analysed using four earthquake records: Displacement spectra and acceleration spectra for all analysed earthquakes, with 0.5 percent of damping, are shown in Figure  24 and 25, respectively. The analysis of these spectra offers a good idea of behaviour of the structure, even before the analysis is conducted, as will be shown later in the paper. Accelerograms of all analysed earthquakes, for earthquake simulation: with the return period of 475 years, scaled to 0.23g; with the return period of 10000 years, scaled to 0.54g.
The dynamic response of the dam was analysed for all four earthquakes for the time step of Dt=0,01 s, with the total of one thousand steps, i.e. the dynamic response of the dam was monitored over a period of ten seconds. Simulation of concrete ageing on dams as illustrated by numerical analysis of Jablanica HPP Figure 26 shows comparative dam crest displacements over time for all analysed earthquakes for a maximum acceleration of 0.23g, when the storage reservoir is full, and Figure 27 shows comparative dam crest displacements over time for all analysed earthquakes for a maximum acceleration of 0.54g, when the storage reservoir is full. Comparative dam crest accelerations over time, also for the full storage reservoir are shown in Figure  28 for all analysed earthquakes for a maximum acceleration of 0.23g, while the same is shown in Figure 29 for all analysed earthquakes for a maximum acceleration of 0.54g.   Figure 32 shows comparative total hydrodynamic forces acting on dam body over time, for all analysed earthquakes, and for a maximum acceleration of 0.23g, when the storage reservoir is full, and Figure 33 shows comparative total hydrodynamic forces acting on dam body over time, for all analysed earthquakes, and for a maximum acceleration of 0.54g, when the storage reservoir is full.

Analysis of results
Accelerograms, spectra and results obtained during analysis clearly show that the earthquakes selected for analysis greatly differ according to their periods, way in which they act, and according to the total effect they have of the dam. Banja Luka is an impact earthquake whose dominant activity lasts for about three seconds, and then the earthquake suddenly calms down for the next two seconds, only to continue with a low-intensity shaking until the end of seismic record. As a result, hydrodynamic pressures for this earthquake are highly variable and at some moments (in some parts of the dam) they attain very high values but, overall, the hydrodynamic force imposed on the dam is not the greatest. The Northridge earthquake, which is more oscillatory compared to the Banja Luka earthquake, is characterised by much smaller displacements and accelerations compared to other earthquakes, but has the greatest hydrodynamic force. It can be seen that the accelerogram of the Northridge earthquake covers very well the free period of the storage reservoir in the time segment between 2.5 and 3.0 seconds, as it then imposes the greatest total hydrodynamic force on the dam. However, as can be seen from the displacement diagram, the dam withstands this force quite well with the maximum relative displacement (measured from hydrostatic state) of approximately 6.5 mm (for the earthquake scaled to 0.54g).
Kobe and Ston earthquakes are somewhere in between these two earthquakes, as can be seen from the response spectrum and the corresponding results. The results of dam behaviour when subjected to the Banja Luka earthquake load are very interesting. Although smaller total hydrodynamic forces imposed on the dam were obtained for this earthquake, for the return period of 475 years (0.23g) the earthquake provoked displacement of approximately 7.0 mm (which is about 60% of the hydrostatic load, with the dam crest acceleration somewhat exceeding 9.81 m/s 2 ). For the return period of 10000 years, the earthquake causes displacement of approximately 12.4 mm (which is even slightly more compared ot the hydrostatic load, with the dam crest acceleration of about 24.6 m/s 2 , which amounts to almost 2.5 g). This earthquake, together with the Kobe earthquake, exerts the greatest influence of the dam. It is interesting that the Northridge earthquake causes twice as low acceleration of dam crest for the same acceleration amplitude. For the return period of 10000 years (0.54g), the Kobe earthquake provides the displacement of approximately 12.0 mm, with a slightly smaller crest acceleration of approximately 20.0 m/s 2 . However, the Kobe earthquake provides the greatest displacements at the dam crest. It is obvious that the dam, due to its considerable thickness, is more sensitive to rapid seismic impacts, while earthquakes with greater periods, despite causing greater hydrodynamic pressures and hydrodynamic forces, actually exert a smaller impact on the dam. Relatively great hydrodynamic pressures occur at the bottom of the dam and so, for the Banja Luka earthquake, the pressure is almost 1000 kPa, with the tension of approximately 600 kPa, for the return period of 10000 years (0.54g). The hydrostatic pressure at the bottom of the dam is approximately 760 kPa, and so the cavitation does not occur. According to numerical analysis results, cavitation occurred at several points at 2/3 of the dam height, but the values were very low and it was observed in a very limited area, which is why it can be attributed to numeric error rather than to an acutal physical occurrence. In addition, it should be noted that these pressures are the so called peak values and that the mean pressure is much lower, i.e. it amounts to approximately ± 200 kPa. As it could be said that the Banja Luka earthquake is, so to speak, a "domestic earthquake", i.e. an earthquake that could most probably be expected in the wider area of the dam, further discussion will actually concentrate on this earthquake, and selected results of dam behaviour will be presented in discrete time segments. Dam displacements caused by the Banja Luka earthquake (0.23g) are presented in Figure 34 for the fully impounded reservoir and for the earthquake duration of 2.25 seconds, when the greatest displacement of 0.68 cm was registered in the negative or downstream direction (positive direction of coordinate system is positioned toward the water storage reservoir). GRAĐEVINAR 71 (2019) 9, 749-767 Simulation of concrete ageing on dams as illustrated by numerical analysis of Jablanica HPP Dam displacements caused by the Banja Luka earthquake (0.54g) are presented in Figure 35 for the fully impounded reservoir and for the earthquake duration of 2.25 seconds, when the greatest displacement of 1.24 cm was registered in the negative or downstream direction, i.e. toward to Neretva River. In this figure, the basic geometry of the dam is marked in gray, while the displaced dam is marked in purple. Displacements are increased by the coefficient of 2000. The formation and propagation of cracks on the dam subjected to the Banja Luka earthquake (0.23g) is shown in Figure 36 for the fully impounded reservoir. The first crack occurred on the dam extrados 0.02 seconds after activation of the earthquake (Figure 36.a). The situation did not change until 2.18 seconds after activation of the earthquake, when two new cracks occurred, one on the extrados and the other on the intrados of the dam (Figure 36.b). Thus, a total of three cracks were registered, two at the extrados and one at the intrados of the dam. The formation and of cracks on the dam subjected to the Banja Luka earthquake (0.54g) is shown in Figure 37 for the fully impounded reservoir. The first crack occurred on the dam extrados 0.02 seconds after activation of the earthquake (Figure 37.a). The situation did not change until 2.16 seconds after activation of the earthquake, when three new cracks occurred at dam crest (Figure 37.b). Four new cracks occurred at 2.17 seconds after activation of the earthquake, two at the extrados slightly below the dam crest and two at the extrados at dam abutment (Figure 37.c). One new crack occurred at 2.23 seconds after activation of the earthquake at the extrados slightly below the dam crest (Figure 37.d). Two new cracks occurred at 2.42 seconds after activation of the earthquake at the extrados at both dam abutments (Figure 37.e). Two new cracks occurred at 2.43 seconds after activation of the earthquake, one at the extrados at dam abutment, and the other at the intrados at dam abutment ( Figure 37.c). Three new cracks occurred at 2.63 seconds after activation of the earthquake, one at the extrados at dam abutment, and the other two at the intrados at both dam abutments ( Figure  37-g).. No new cracks appeared at 2.63 seconds after activation of the earthquake. Thus, a total of sixteen cracks occurred, ten at the extrados, three at the intrados and three at dam crest. It is interesting to note that the final number of cracks (after the earthquake action) was relatively low. This can be explained as follows: First, although according to some classifications this dam can be grouped under arch dams, based on its relatively low height and cross sectional dimensions, it can actually be classified among gravity arch dams. In addition, initial stresses in the dam are relatively small for hydrostatic state so that an earthquake, even the one for the return period of 10000 years, would not generally cause extreme stress levels and cracking..

Some conclusions relating to dam analysis
After analysis of the previously conducted numerical simulations, it can be concluded that Jablanica dam is a globally safe and stable structure, even after sixty years of intensive use. Some cracks that have occurred at inspection galleries and at other places along the dam are more due to local stress concentration around the opening, change in temperature, and shrinkage, than to global instability of the dam. On the other hand, a better model, which is to be developed in the future, and is to include parallelisation of code and analysis of a much greater number of finite elements, could also simulate inspection galleries and shafts, and a greater number of elements along the dam thickness, and is expected to provide an even more accurate answer about the real state of stress in the dam. In any case, the dam simulation across its thickness should be conducted with at least three elements, and the size of a finite element should not exceed 1.0 m.
What should definitely be included in some subsequent phases of analysis of this structure is certainly a better simulation of the surrounding soil, and this primarily from the standpoint of possible occurrence of faults, connection of different lithological environments, etc. Furthermore, the program should certainly contain contact elements that could better simulate the link between the dam and the bedrock. The results of the analysis conducted with concrete parameters that include ageing correspond very well with the dam monitoring results [25], and with the results obtained in [29]. It can be said that the developed model is a "powerful" tool for the analysis of concrete dams, and that it is especially appropriate for older concrete dams.

Conclusion
The problem considered in this study, i.e. simulation of dynamic interaction of the fluid-structure system is highly complex, as already presented in [1], because it encompasses all aspects of fluid mechanics and mechanics of deformable solids. The problem of numerical simulation of the fluid-structure interaction still remains an insufficiently investigated area, although it is currently being intensively developed. But not without a good reason. The growing need for drinking water and water for irrigation, as well as the need for clean water-generated energy, has intensified the efforts to design and build new dams, and to properly maintain the existing ones. An intensive construction of dams for generation of hydroelectric power was initiated in the early 20th century. The dams built in that period are now more than one hundred years old. The issue involving their safety is currently highly topical not only in professional circles, but also in political and safety-concerned circles. The model presented in this paper, which monitors development and pattern of change of mechanical properties of concrete over time (concrete ageing), has proven to be especially significant in the analysis of older structures (dams) The analysed example, involving simulation of behaviour of a real-life dam under seismic conditions, is the proof of the validity and efficiency of the developed model and of the corresponding computation software. Practical application of such models for simulating dynamic fluidstructure interactions has greatly contributed to the study of real behaviour and safety of new and existing dams under dynamic/ seismic conditions.