Investigation of Joints Behavior in Seismic Response of Arch Dams

There has been extensive research to examine the effects of contraction joints opening in the nonlinear seismic response of arch dams. However, there has been less attention on the effects of perimetral joint on the response. In this study, a special finite element program is developed, that is described and verified initially. Based on this program, the nonlinear dynamic behavior of a typical thin arch dam is studied. The perimetral joint, as well as, contraction joints are included in the nonlinear cases analyzed. It is shown that the distributions of the maximum tensile stresses are very sensitive to the properties used for the joints. Moreover, the modified joint model proposed for the modeling of perimetral joint is found to be less sensitive and very effective in seismic analysis of concrete arch dams.


Introduction
During the last two decades, many researchers have studied the effects of contraction joints opening in the seismic response of concrete arch dams by different techniques.However, there has been less attention on the effects of perimetral joint on the response.In some of the previous studies such as the work of Fenves et al. [1], they have utilized isoparametric interface elements with a relatively simple separation model relying on the assumption of linear shear stresses.In other studies, e.g., the work of Hohberg [2], more advanced models are developed based on plasticity theories, which are less experimented in the actual applications on earthquake response of concrete arch dams.
In this study, a special finite element program is developed [3] to investigate joints behavior in response of arch dams due to earthquakes.Two nonlinear options are available in the program.A separation alternative, which simply models joint opening and closure while shear stresses in the joint are assumed to behave linearly, and the second option which is separation along with a shear stresses release mechanism.Both models implemented are relatively simple and efficient with no convergence problems.
In the upcoming sections of this paper, the basic concepts and the nonlinear models employed in the program are explained briefly at first.Then, a simple example is considered which helps to illustrate the differences in behavior between these two models, as well as, verifying the developed program.Later on, the nonlinear dynamic behavior of a typical thin arch dam is studied by the application of the models discussed.The perimetral joint, as well as, contraction joints are included in the nonlinear cases analyzed, and parametric study on the properties of both types of joints are performed.e e J J S S E E International

Basic Concepts 2.1 Nonlinear Dynamic Analysis
Consider the equation of dynamic equilibrium at instant n+1: M and C are total mass and damping matrices of the system.Meanwhile, F n+1 and R n+1 in this equation are the vector of nodal forces equivalent to internal stresses and the nodal external forces respectively.
It could be easily shown that by applying Newmark's method, equation (1) can be transformed to the following equation: Based on this equation, the (i+1) th correction for displacements vector is obtained through an iterative process which is explained elsewhere in details [4].K t in this equation is the tangential stiffness matrix of the system.The total tangent stiffness matrix of the system is obtained by assemblage of different elemental stiffness matrices which can be defined for a given discrete crack model.Assuming that non-linearities are limited to contraction and perimetral joints, it is clear that solid elements can be taken as usual linear isoparametric elements whose stiffnesses are calculated once at the beginning, assembled and stored for recursive application throughout the analysis.While interface elements stiffnesses need to be updated at each iteration.

Isoparametric Interface Element
The interface element is utilized to model contraction and perimetral joints in this study.It consists of two spatial 8 noded isoparametric layers A, and B placed originally on top of each other with potential for being separated partially or completely.The tangent stiffness of these elements is obtained by the following integration [4]: In which T, B and G are matrices of rotation, strain interpolation and transformation of displacements to strains (relative displacements of the two surfaces).Matrix D l provides the relation between local stresses and relative displacements in two arbitrary tangential directions and the normal direction to the joint surface.It is a 3 × 3 diagonal matrix with 3 1 j ; k D j l jj − = = k j being the stiffness coefficient in the j th direction.In the finite element program, MAP-76 [3], there are three options available for the above constitutive relations that are explained briefly below:

Linear Model (NCRIT=0)
In this case, the stiffness coefficients k j are set equal to relatively large numbers which are kept constant during the analysis.Therefore, the stiffness matrix would remain constant such as any other linear element.This implies that even though the joints are modeled, but the results in this case should theoretically be similar to a model which joints are not considered.Of course, this would happen if the penalty coefficients (k j ) are chosen correctly.
In general, this option can be used to investigate whether the stiffness coefficients are selected appropriately or not.Meanwhile, this option can be activated in case the natural frequencies of the system are to be calculated.

Separation With Linear Shear Stresses Behavior Model (NCRIT=1)
This is the behavior considered in most of the previous studies in this subject, such as the work of Fenves et al. [1].In this option, the stiffness coefficients k j are set to relatively large numbers.However, in case the normal stress ( σ l 3 ) at a Gauss point is greater than a specified limit σ * (tensile strength of the joint), perfect tensile softening is assumed, and the element opens up.This means k 3 =0 and the tensile strength at that point is reduced to zero.Of course, it regains its stiffness if that Gauss point tries to act in compression at a later stage.During this whole process, stiffness coefficients corresponding to tangential directions are kept constant.This can be presumed as a separation along with linear shear stresses behavior law, which is a realistic model for contraction joints when shear keys are anticipated.Of course, assuming that shear keys are not damaged, and the amounts of separation are not larger than their heights.

Separation With Nonlinear Shear Stresses Behavior Model (NCRIT=2)
In this case, the stiffness coefficients k j are set initially equal to relatively large numbers as usual.However, in case the normal stress ( σ l 3 ) at a Gauss point is greater than a specified limit σ * (tensile strength of the joint), perfect tensile and shear softening are assumed, and the element opens up.This mean k j =0; j=1-3, and the tensile and shear strengths at that point are reduced to zero.Of course, it regains its stiffnesses in all three directions if the normal stress at the Gauss point tries to move into the compression range at a later stage.Therefore, the separation is occurring along with shear stresses release mechanism.
It should be mentioned that, this concept can be perceived as utilizing Mohr-Coulomb law for shear strength of the joint.Assuming that cohesion factor and friction angle are taken as a large number and 90 degrees respectively, and as the tensile damage occurs, the cohesion factor is reduced to zero.This is a more realistic model for contraction joints without shear keys, or in case of perimetral joints.Of course, it is noted that, this can also be seen as an approximate slippage model.That is, it allows unrestrained relative tangential displacement after a tensile damage.While, if the normal stress at the joint goes into the compression range, the shear strength would be infinite, having in mind that the friction angle is assumed 90 degrees.Meanwhile, it also regains its shear stiffnesses similar to normal direction.

Verification Example
As mentioned in the previous section, two nonlinear models are available for the interface elements.A separation alternative (NCRIT=1) which simply models joint opening and closure while shear stresses are behaving linearly, and the second option which is separation with nonlinear shear stresses behavior model (NCRIT=2).For an initial verification of the developed program and to illustrate differences in behavior between these two nonlinear models, the following example has been defined.Consider a 16 noded isoparametric interface element with dimensions of 2× 2 meters (Fig. 1a) acting based on the two above mentioned nonlinear material behavior options.The basic material data for this element include the normal stiffness coefficient k n =1500 GPa/m, tangential stiffness coefficients k s,t =k n /2, and the initial tensile strength of the joint which is selected as The nodes of the lower plane (1-8) are kept fixed, and specified displacements are enforced for the nodes of the upper plane (9-16) in horizontal and vertical directions (x, z) while they are fixed in the y-direction.The applied control displacements path are illustrated in Fig. 1b for e e J J S S E E International node 13, and the other nodes of the upper plane are in similar conditions.It is noted that vertical displacements are enforced initially (Path O-A), then horizontal displacements are applied with vertical displacements being kept constant (Path A-B).The controlled displacements are continued in a decreasing manner till it reaches point C, which has a zero vertical displacement or the two planes are in contact again.In the next step, displacements are varied with the same rate down to point D. The history of normal and shear stresses are plotted in Fig. 2 for the two above mentioned nonlinear material behavior options (NCRIT=1 and 2).
In the first case, separation with linear shear stresses behavior model (NCRIT=1), it is observed that normal stress increases up to 2 MPa, and suddenly tensile crack occurs and the normal stress drops to zero afterward.This stress remains zero during the path A-B-C as the joint is open.In path C-D, as the joint closes and the two surfaces compress against each other, the normal stress moves into the compression range till it reaches a maximum compressive value of -6 MPa.As it is observed, the shear stress variation is linear in this case during the whole path O-A-B-C-D.
In the second case which is separation with nonlinear shear stresses behavior model (NCRIT=2), it is observed that history of normal stress is similar to the previous case.However, shear stress behavior is nonlinear in this case.As in path O-A, the shear stress is zero, due to displacements being merely imposed in vertical direction.Meanwhile, the shear strength is set to zero due to tensile crack.Therefore, during path A-B-C, shear stress behaves similar to normal stress and remains zero, as the joint is open in this range.In path C-D, the shear stress increases linearly to -2 MPa, since the joint is in a closed state.
The normal and shear stresses calculated in this example could be easily checked by hand calculations, knowing the joint stiffness coefficients and specified displacements.This assures the fundamental verification of the developed program and the material models implemented for interface elements in this study.

Application on Concrete Arch Dams
In this section, the nonlinear behavior of Shahid Rajaee arch dam is being studied by application of the models discussed above.The dam is 130 m high, with the crest length of 420 m and it is located in north of Iran in the seismically active foothills of Alborz Mountains, near the city of Sari.

Finite Element Idealization
Selected Models Five cases are considered.A linear model, case J1 and four discrete crack cases which are referred to with their symbolic names J6, J7, J8 and J9.The finite element mesh for discrete crack cases (Fig. 3) include 76 isoparametric 20 noded solid elements, and 84 isoparametric 16 noded interface elements.This model consists of a total of 873 nodes and 2316 degrees of freedom.The mesh for the linear case (J1) is similar to the mesh for discrete crack cases excluding interface elements.
Although the geometry of the dam is not exactly symmetric, but it was decided to idealize it as a symmetric model based on average parameters of the two sides.It is clear in this case, discretization might be performed on one half of the dam and conditions of symmetry or anti-symmetry would be applied in the mid-plane for different components of earthquake excitations.However, this was disregarded and the whole dam is discretized in this study.Therefore, symmetric results are expected throughout the analyses knowing that cross canyon excitation is also neglected, and this fact is used as an additional check for the accuracy of obtained results.Dam-reservoir interaction is considered based on modified Westergaard's method [5].
Meanwhile, the foundation is taken as rigid.Both of these idealizations were decided to reduce computational efforts.Of course, it is clear that rigid foundation assumption would create large tensile stresses near the boundaries in the linear case.However, this was seen to be a more challenging test for the nonlinear models in the attempt to reduce these high stresses.
In the analyses carried out, the Rayleigh damping matrix is applied and the corresponding coefficients are determined such that equivalent damping for frequencies close to the first and sixth modes of vibration would be 10% of the critical damping.

Basic Parameters
The concrete is assumed to have the following basic properties: Elastic modulus E c = 30.0Gpa Poisson's ratio ν c = 0.18 Unit weight γ c = 24.0kN/m 3   Interface elements utilized for discrete crack models are applied with different stiffness coefficients and nonlinear behavior laws (NCRIT) mentioned in Table 1 below.The water is taken as incompressible, invisid fluid, with weight density of 10.0 kN/m 3 and the water level to be at elevation 485.0 masl (h=122.0m).

Loading
It should be mentioned that static loads (weight, hydrostatic pressures) are visualized as being incrementally increasing in time until they reach their full magnitude.Therefore, the same time step of 0.01 second, which is chosen in dynamic analysis, is also considered as time increment of static loads application.It is noted that time for static analysis is just a convenient tool for applying the load incrementally, but it is obvious that inertia and damping effects are disregarded in the process.In this respect, the dead load is applied in ten increments and hydrostatic pressures thereafter in another ten increments at negative range of time.At time zero, the actual nonlinear dynamic analysis begins with the static displacements and stresses being applied as initial conditions.
The dynamic excitations include two components of Friuli-Tolmezzo earthquake records (cross-canyon excitation neglected to have a symmetric loading) normalized based on the frequency content for MDE condition with a peak ground acceleration of 0.56g in the horizontal direction.The time duration considered is 12 seconds in each case.

Analysis Results
As mentioned in the previous section, five cases are considered, the linear model (LN), and four discrete crack (DC) nonlinear models.In each case, envelopes of maximum tensile and compressive principal stresses are obtained, and the results of maximum stresses occurring on both faces are summarized on Furthermore, the envelopes of maximum tensile stresses are illustrated in Figs. 4 and 5 for cases J1 and J6.For the linear case (J1), it is observed that maximum tensile stresses for the spillway and abutments regions reach to maximum values of 8.9 MPa and 15.72 MPa, respectively (Fig. 4).The maximum tensile principal stresses of these zones occur in the arch direction and perpendicular to the abutments as expected.The maximum compressive stresses of this case are -15.24MPa and -11.91 MPa for the spillway and abutments regions, respectively.
For the case J6, which is a discrete crack model, it is noted that maximum tensile stress for the whole analysis is limited to 5.85 MPa, which is occurring approximately 25 meters below the spillway on the downstream face (Fig. 5).It is clear that high tensile stresses in the spillway and abutments regions observed in the linear case, are released in this case by opening of the joints.
However, as the central cantilevers move toward upstream, a maximum tensile stress of 5.85 MPa develops in the cantilever direction on downstream face.Of course, this limited stress is still higher than the tensile strength of concrete lifts or mass concrete itself.At the base and abutments, the maximum tensile stress observed is a much lower value of 2.61 MPa.
In regard to compressive stresses (Table 2), it is noted that maximum value on the upstream face of the spillway region is decreased, while the value on the abutments at downstream face is increased slightly in comparison to the linear case.
The history of displacements at the right quarter point and the central point at the dam crest are also monitored through time and the maximum values are summarized in Table 3.It is noted that maximum value for the case J6 is for the central point along the stream direction (12.02 cm), and it is increased slightly in comparison to the linear case J1.Furthermore, the value of maximum joint openings at the crest for the middle and 1/4 point's joints can be found in Table 4.The maximum value for the case J6, occur at the downstream side of the middle joint with a magnitude of 2.34 cm.

e e J J S S E E
International Let us now concentrate on the results of case J7, which the parameters and nonlinear behavior model used for the joints are similar to the case J6, except for the shear stiffness coefficients.The ratio of k s,t /E=2 is chosen in this case (Table 1), that is one half of the value used for the ratio of normal stiffness coefficient which is compatible to Hohberg's suggestion [2].
The envelope of maximum tensile stresses for the case J7 are presented on Fig. 6.It is observed that overall distribution of stresses and maximum values in most regions are close to the case J6.However, very high tensile stresses (σ 1 =8.50 MPa) are noticed at the base and abutments intersections on the upstream face.This phenomenon, which is referred to as joint locking effect is due to the fact that shear stiffness coefficients are too high, and the perimetral joint can not easily open on these corners.This becomes more clear in the case of 90 degrees kinks.Because, in that case it is easy to see that opening of the joint on one face corresponds to tangential displacement on the adjacent perpendicular interface element.Therefore, for high shear stiffness coefficient, the joint is behaving similar to being locked.In the spillway and middle portion of the dam, it is noticed that contraction joints are behaving appropriately.The high arch stresses underneath the spillway are released and the maximum tensile stress in this region is 6.63 MPa, which occurs on the downstream face in the cantilever direction.
As for the compressive stresses, displacements and joint openings at the crest (Tables 2-4), the behavior are similar to case J6, and the maximum values are not changed significantly.
The next case is J8, which is similar to J7 except that stiffness coefficients are decreased by a factor of 100 (Table 1).The results are summarized on Tables 2-4, and the envelopes of maximum tensile stresses are presented on Fig. 7.It is noticed that significant differences exist as far as the maximum tensile, compressive stresses, displacements and joint openings in comparison with the case J6.
To see the differences more vividly, the history of displacement at mid-crest along the stream direction and the adjacent joint opening are compared in Fig. 8 for the two cases of J6 and J8.It is obvious that the joint stiffness coefficients are two low in the case of J8, such that there are large negative joint openings, which are not realistic.Meanwhile, the maximum displacement of mid-crest along the stream direction show an increase of about 40% for the case of J8 in comparison with J6, and reaches a high value of 16.83 cm.
The final considered case is J9, which has parameters similar to the case J7.That is both normal and shear stiffness coefficients are high values.However, in the recent case, the joint modeling behavior for the perimetral joint is changed to a case of separation with nonlinear shear stresses behavior option (NCRIT=2).
The results are summarized on Tables 2-4, and the envelopes of maximum tensile stresses are presented on Fig. 9.It is noticed that the distribution of tensile stresses and the maximum values are very similar to the case J6.Although a high value is chosen for the shear stiffness coefficients, but there are no sign of joint locking phenomena.This is mainly due to the improved joint modeling (NCRIT=2) employed in this case.
As for the maximum compressive stresses, displacements and joint openings (Table 2-4), the results are also similar to the case J6.For a better evaluation of the results, it was decided to compare the cases from natural frequency's point of view.Therefore, an Eigen-value problem is solved in each case based on the initial stiffness of the models, or it could perceived that linear option (NCRIT=0) is used temporarily for discrete crack models as far as these computations are concerned.Clearly, these values are valid as long as there are no joint separation in the real cases.Meanwhile, the natural frequencies of nonlinear cases should theoretically be very close to the ones related to the linear case, if the joint stiffness coefficients are appropriate.
The two initial natural frequencies of different cases are presented in Table 5 below, as well as the percentage difference of the first natural frequency of each case in comparison to the linear case.It is observed that in the case of J6, the percentage difference in fundamental natural frequency is only 2.5%.This comparison and the fact that there is not large negative joint opening show that joint stiffness coefficient values are selected in a suitable range in this case.Meanwhile, due to low value of shear stiffness coefficient, there was no joint locking problem.
As for the case J7, it is noticed that the percentage difference in fundamental frequency is even decreased to a low value of 0.6%.However, as it was seen, joint locking phenomena originated in this case due to large value of shear stiffness coefficients, and lack of shear stress release in zones having tensile cracks.This caused very high tensile stresses in the upstream face near boundary kinks, even though far from these corners, the maximum tensile stress values are close to the previous case (J6).
In the case J8, the percentage difference in natural frequency is a large value of 22%.This is consistent with unrealistic results observed in that case, such as the large negative joint openings mentioned.Therefore, the percentage difference in fundamental frequency of the system in comparison to the linear case could be used as a basic criteria to evaluate appropriate selection of joint stiffness coefficients.However, it should be noticed that joint locking phenomena due to its local effects cannot be predicted by these values.Furthermore, it is observed that joint locking effects can be avoided by selecting a relatively low value of shear stiffness coefficient, or a more rigorous solution is to utilize a modified joint model (NCRIT=2) for perimetral joints.In the second alternative, the results are even less sensitive, and high values of shear stiffness coefficients can be applied without creating any joint locking.Meanwhile, it improves the value of percentage difference in fundamental frequency, as it was obtained in the case of J9 (0.6%).

Conclusions
The nonlinear dynamic behavior of Shahid Rajaee concrete arch dam is subjected to the Friuli-Tolmezzo earthquake records.Five cases are analyzed, a linear model and four discrete crack models based on different joint models and joint stiffness properties.Overall, the main conclusions obtained can be listed as follows: • Although, the criteria suggested in literature [1 and 2] for joint normal stiffness coefficient is suitable, similar values for joint shear stiffness coefficients are not e e J J S S E E International appropriate in cases where perimetral joint is also modeled.In these cases, joint locking phenomena are observed which causes very high tensile stresses at the base and abutment intersections on the upstream face.
• Joint locking distort the results in a local manner, and far from the boundary kinks, the effects are negligible.

•
In the usual joint models (NCRIT=1), if the perimetral joint is considered, the shear stiffness coefficients should be decreased to a value of 0.025 times what has been suggested by Fenves et al. [1], to avoid joint locking phenomena.
• A high value of shear stiffness coefficient can be applied in the modified joint model (NCRIT=2), without any joint locking effects.

•
The percentage difference in fundamental frequency of the system in comparison to the linear case, could be used as a basic criteria to evaluate appropriate selection of joint stiffness coefficients.However, joint locking phenomenon due to its local effects could not be predicted by this parameter.

Fig. 8b
Fig. 8b Comparison of joint openings at upstream face of mid-crest for cases J6 and J8.

Fig. 9a
Fig. 9a Envelope of maximum tensile stresses (MPa) for the case J9 (upstream view)

Table 1
Discrete crack cases analyzed and the properties used

Table 2
Maximum principal stresses (MPa) for different cases

Table 3
Maximum displacements (cm) at crest of the dam

Table 4
Maximum joint openings (cm) at crest of the dam

Table 5
Natural frequencies (Hz) of different cases