INTRODUCTION
Arch dams are huge complicated facilities built mainly for water and electric
power supply usually in dry and semidry climate countries such as Iran. Karun1
dam is an aged high elevated concrete arch dam located in north western of Khouzestan
province of Iran country near the MasjedSoleiman city which its building was
completed in 1962. Under the effect of operation and climate conditions as well
as dam body internal chemical reactions such as cement hydration and alkalisilica
reaction, Karun1 aged high elevated concrete arch dam experiences some of changes
in its structural characteristics like concrete modulus of elasticity of dam
body. In recent years, safety control of this dam is highlighted due to the
fact that the electrical energy generation of the dam has been increased thorough
the development of the second phase of its power plant. In order to prevention
of initiation and growth of thermal cracks in dam body, it has been normal procedure
to provide vertical contraction joints in arch dams at approximately 15 to 20
m spacing. In this way, continuous dam body is separated to several vertical
blocks mainly named monoliths in scientific literatures. These contraction joints
were injected with cement grout after the concrete has cooled to mean temperature.
It is worth to mention that adjacent blocks has no any relative shear motions
and even under strong ground shakings some little dislocations were reported
and many studies conducted till now were focused on the investigation about
the effects of contraction joints on earthquake response of the arch dams (Fenves
et al., 1992; Zhang et al., 2000; Xinjia
et al., 2002; Arabshahi and Lotfi, 2009;
Jing and Chen, 2009; Sheng et
al., 2009; Li et al., 2008; Zhao
et al., 2007; Du and Tu, 2007). However,
an arch dam in compare to gravity dam has a wider surface and less thickness,
so, temperature differences between upstream and downstream faces of these dams
can lead to creation of larger stress and strains in relative to their corresponding
values in gravity dams. In a few works has been performed to date on thermal
analysis of arch dams, one dimensional heat transfer across the vertical section
of dam has been studied in which the influence of vertical contraction joints
were not included (Sheibani and Ghaemian, 2006). They
incorporated the thermal stress field in Karaj arch dam in Iran. They contributed
in their study the effect of solar radiations as well as other resources of
heat generation in dam like air and reservoir temperature changes. They concluded
that two dimensional thermal analysis of an arch dam cannot yield accurate results
and 3D numerical simulation is needed. Ardito et al.
(2008) investigated the effect of seasonal thermal loading as well as hydrostatic
pressure on an arch dam in Italy. They used truncated Fourier series for heat
transfer calculation across dam thickness. They computed damage diagnosis thorough
minimization of a batch discrepancy function between measured and computed displacements
as inverse analysis in a linear thermoelasticity context. Leger
and Leclerc (2007) presented frequency domain solution algorithms of one
dimensional transient heat transfer equation to accommodate thermal analysis
of Schlegeis arch dam. They developed two algorithms to solve two kinds of heat
transfer problems: direct problems in which temperature variations are specified
on both upstream and downstream faces of dam and inverse problems where temperatures
are measured with thermometers located inside instrumented dam sections. Leger
and Seydou (2009) presented the hybrid dam displacement model to simulate
the seasonal thermal displacements of gravity concrete dams. They used their
model for predict the behavior of dam under extreme thermal effects which cannot
be happened in real conditions.
According to what mentioned above, in this research a try has been done to study the effect of vertical contraction joints on the Karun1 dam safety certainly under effect of thermal loads by means of a relatively 3D exact simulation of the geometric, material behavior and boundary conditions of the dam using ABAQUS finite elements package. ABAQUS is commercial powerful software for finite elements analysis developed by HKS Inc of Rhode Island, USA and marketed under the SIMULIA brand of Dassault Systemes DS. The ABAQUS product suite consists of three core products: ABAQUS/Standard, ABAQUS/Explicit and ABAQUS/CAE. ABAQUS/Standard is a generalpurpose solver using a traditional implicit integration scheme to solve finite element analyses. ABAQUS/Explicit uses an explicit integration scheme to solve highly nonlinear transient dynamic and quasistatic analyses. ABAQUS/CAE provides an integrated modeling (preprocessing) and visualization (post processing) environment for the analysis products. The ABAQUS is used in the automotive, aerospace and industrial product industries. The product is popular with academic and research institutions due to the wide material modeling capability and the program's ability to be customized. Principal stress tensors, displacement vectors and strain energy were selected as stability indexes safety control and examined. It will be shown that the contraction joints have significant effects on the thermostatic behavior of Karun1 dam.
MATERIALS AND METHODS
As it has been emphasized, in this study we are going to investigate the effect
of vertical contraction joints embedded in a typical arch dam on its thermomechanical
behavior. To accommodate this issue, thorough May to November 2009 the contraction
joints of Karun1 dam were modeled using smallsliding contact feature of wellknown
scientific package ABAQUS. In smallsliding contact theory it is assumed that
two contact surfaces has only relatively small sliding relative to each other
but arbitrary rotation of the bodies is permitted.

Fig. 1: 
Slave nodes interacting with a two dimensional surface 
In this program two kind of surface is defined: slave surface and master surface
through them we can define the contact physical phenomenon. A kinematic constraint
that the slave surface node can not penetrate the master surface is enforced.
Four internal contact elements designed in the program to handle all kinematic
constraint conditions. The three dimensional contact between a slave node and
a deformable master surface is selected for the proposed study. For better understanding
of relations which will be come into next discussion, the formulation of this
kind of contact is presented as below (ABAQUS, Inc. and Dassault Systemes).
In this case, a two dimensional slave node interacting with a first order master
surface. This formulation can be generalized to second order as well as three
dimensional situations.
As shown in Fig. 1, the nodes 102, 103 and 104 are located on slave surface and nodes 1, 2 and 3 are positioned on master surface. In the first step, unit normal vectors are computed for all the nodes on the master surface. For instance, the unit normal N_{2} is computed by averaging the unit normal vectors of segments 12, 23. The user can also define unit normal vector for each node on master surface. Software assigns additional unit normal vectors for each segment a distance εL from each end of the segment. N_{ε} is the unit normal vector in the middle of segment 12. The unit normals computed are then used to define a smooth varying normal vector, N (X), at any point, X, on the master surface. An ‘anchor’ point on the master surface X_{0} is computed for each slave node so that the vector formed by the slave node and X_{0} coincides with the normal vector N (X_{0}). Assume that the anchor point of slave node 103 is X_{0} on the segment 12. Then we define:
where, X_{1} and X_{2} are the coordinates of nodes 1 and 2, respectively and u_{0} is calculated so that X_{0}X_{103} coincides with N(X_{0}). Furthermore, the contact plane tangent direction, v_{0}, at X_{0} is selected so that it is perpendicular to N(X_{0}); i.e.,
where, T is a (constant) rotation matrix.
The smallsliding contact constraint is achieved by requiring that slave node 103 interact with the tangent plane whose current anchor point coordinates are, at any time, given by:
where, N_{1} (u_{0}) = 1u_{0} and N_{2} (u_{0}) and whose current tangent direction is given by:
where, N_{1}^{u} (u_{0}) = 1 and N_{2}^{u}
(u_{0}). Since, the above expressions for the point x_{0} and
the vector v resulted from barycentric (affine) combination of the points x_{1}
and x_{0} that is (Christopher, 2007):
The contact plane will be mapped properly under affine transformations such
as translation, scaling (stretching) and rotation. At each slave node that can
come into contact with master surface we construct a measure of overclosure
h (penetration of the node into the master surface) and measures of relative
slip s_{i}. These kinematic measures are then used, together with appropriate
Lagrange multiplier techniques, to introduce interaction theories as will be
described later. In three dimensions, the overclosure h along the unit contact
normal n between a slave point x_{N+1} and a master plane p (ζ_{1},
ζ_{2}), where, ζ_{i} parametrize the plane, is determined
by finding the vector (px_{N+1}) from the slave node to the plane that
is perpendicular to the tangent vector v_{1} and v_{2} at p.
Mathematically, we express the required condition as:
When:
If at a given slave node h<0, there is no contact between the surfaces at
that node and no further surface interaction calculations are needed. If h≥0,
the surfaces are in contact. The contact constraint h = 0 is enforced by introducing
a Lagrange multiplier, ,
whose value provides the contact pressure at the point. To enforce the contact
constraint, we need the first variation δh and for the Newton iterations,
we need the second variation dδh. Likewise, if frictional forces are to
be transmitted across the contacting surfaces, the first variation of relative
slip δs_{i} and the second variations dδs_{i} are
needed in the formulation. The three dimensional smallsliding deformable contact
formulations presented here (Chen and Hisada, 2007;
Puso and Laursen, 2002). A point on the contact plane
associated with a slave node x_{N+1} is represented by the vector:
where, the plane’s anchor point x_{0} and its two tangent vectors v_{1} and v_{2} are functions of the current master node coordinates x_{0}, …, x_{N}. Linearization of Eq. 6 yields:
where, δx_{0} = f (δx_{1},..., δx_{N}) and δv_{i} = g_{i} (δx_{0},..., δx_{N}).
Taking the dot production of Eq. 9 with n results in the following expression for δh:
Likewise, taking the dot product of Eq. 9 with t_{1} = v_{1}/v_{1} and setting h = 0 results in the following expression for the variation of the first slip component:
Similarly, taking the dot product of Eq. 9 with t_{2} = v_{2}/v_{2} and setting h = 0 results in the following expression for the variation of the second slip component:
After formulation of smallsliding contact theory, the next step is the surface
interaction description. For modeling the normal and inplane behavior of monoliths
interactions in Karun1 arch dam body, in this research, the Hard contact and
Coulomb contact theory were implemented as constitutive models respectively
(Franco and RoyerCarfagni, 2005). In Hard contact hypothesis,
we have:
where, p is the contact pressure between two surfaces and h is the overclosure at a given point. The contact constraint is enforced with a Lagrange multiplier representing the contact pressure in a mixed formulation. The virtual work contribution is:
And the linearized form of the contribution is:
And in Coulomb constitutive models below relations are used:
where, τ_{eq} is the equivalent shear stress at a point on the contact surface and τ_{1} and τ_{2} are the shear stresses in two orthogonal directions in contact surface at the same point. If the τ_{eq} is less than a critical value τ_{crit} which can be defined based on any constitutive rational computations, no relative motion occurs between two contact surfaces. In this study, the τ_{crit} is defined as:
In which, μ is the friction coefficient that is defined simply as the:
where, Φ is the internal friction angle between two adjacent monoliths
(vertical contraction joints) surfaces. In this research, the μ is set
to be equal 10 to guarantee that no any relative shear movements are possible
at contraction joints. To complete the behavior definition of interaction between
adjacent blocks in our arch dam, we must describe the heat flow across an interface
via conduction. It is assumed that heat transfer can occur only in the normal
direction. Heat conduction across the vertical contraction joints is supposed
as below in the present study (Sheibani and Ghaemian, 2006):
where, q is the heat flux per unit area crossing the joint from point A on one surface to point B on the other, θ_{A} and θ_{B} are the temperatures of the points on the surfaces and k is the gap conductance. The derivatives of q are:
Where:
The contribution to the variational statement of thermal equilibrium is:
where, A is the area. The contribution to the Jacobian matrix for the Newton solution is:
Where:
In this study, we consider the tied thermal contact. This means that the temperature at point A is constrained to have the same temperature as point B. the Lagrange multiplier method is used to impose the constraint by augmenting the thermal equilibrium statement as follows:
where, λ is the Lagrange multiplier. The contribution to Jacobian matrix for Newton solution is:

Fig. 2: 
Numerical solution consequences 
At this stage, it is supposed that the theoretical background necessary for
doing our proposed problem is approximately prepared. Flowchart of numerical
solution based on the mentioned theory is outlined as Fig. 2.
Figure 2 shows that solving this problem encounter nonlinear
coupled displacementthermal boundary value analysis. Solution involves an incremental
procedure and convergence in each increment obtained only after that equivalent
internal nodal force vector is balanced with equivalent external nodal force
vector within desirable tolerance.
RESULTS AND DISCUSSION
Four model of Karun1 arch dam was simulated in order to reach better clarify of the effect of vertical contraction joints in thermostatic dam behavior. Model A considers dam body as a continuous volume without seeing the vertical adjacent block interactions in two loading mode: with (A1) and without (A2) thermal forces. Model B involves the same two loading condition (B1 and B2) with this difference that vertical contraction joints between monoliths have been simulated in the analysis.
In Fig. 3, the 3D finite elements model of dam and its abutments
is shown. The types of the elements used in this mesh were selected as C3D20RT
for dam body and for abutments as C3D10MT.

Fig. 3: 
Finite elements model of Karun1 concrete arch dam 

Fig. 4: 
Vertical contraction joints in Karun1 dam 
Total number of nodes of the mesh was 110253 and corresponding value for the
elements was 47102. Figure 4 indicates the positions of the
contraction joints in Karun1 dam. Upstreamdownstream displacements along the
highest vertical section of Karun1 dam which is named usually as central cantilever
in arch dam engineering literature have been shown in
Fig. 5. The height level of bottom and top of this cantilever are 342 and 542 m from sea level, respectively. Figure 6 shows the trend of movements along central cantilever in gravity direction. Four model results were compared in this picture. Changes in displacements at crest level of dam in river direction have been shown by Fig. 7. In Fig. 7, the horizontal axis indicated the distance of nodes from central cantilever. Variations of arch stresses along crest of the dam can be checked over in Fig. 8. The same variation for cantilever stresses can be found in Fig. 9.
Regard to curves presented in Fig. 5, it can be deduced that
the movements of dam body along the central cantilever in upstreamdownstream
direction have been increased when contraction joints are considered in modeling
in compare to situation which effect of these joints has been neglected in simulation
process. The extent of this increase is about approximately 100% for two mode
of loading (thermostatic and only static loads). Another important result which
can be concluded from this figure is that under the effect of thermal loads
the body of the dam is moving back toward the upstream direction. This effect
as it is obvious from the picture is significant. When thermal forces are included
in analysis, the maximum displacement decreases approximately 50%.

Fig. 5: 
Movements along the Karun1 central cantilever in river direction 

Fig. 6: 
Movements along the Karun1 central cantilever in gravity
direction 

Fig. 7: 
Crest displacements in river direction 

Fig. 8: 
Arch stresses variations along the dam crest 

Fig. 9: 
Cantilever stresses alterations along central cantilever 
Furthermore, an interesting phenomenon which can be seen is that the effect
of vertical contraction joints and thermal forces in Karun1 arch dam is about
the same, first in increasing and the second in decreasing the displacements
in river direction (the curves belong to model A1 and B2 are very close to each
other). So, we can say that in the analysis of Karun1 dam if any of two important
factors e.g., contraction joints or thermal forces are not considered in data
processing, the significant error in computing river direction displacements
can be occurred. Another trivial result which can be observed from Fig.
5 is that from the base up to crest of the dam the displacements in river
direction are increasing and maximum movements are happened in dam crest. When
thermal loads are added to other existing static loads, e.g., gravity and hydrostatic
loads, the central cantilever of the dam is moving upwards opposite to gravitational
direction. This concept can be deduced from Fig. 6. In other
words, this means that under the effect of heat transfer from downstream to
upstream face of the dam in critical situation in summer season, the dam is
turning toward the upstream direction about the basis. On the other hand, contraction
joints cause the vertical displacements of dam become greater.
Figure 7 shows that displacements in river direction are
symmetric about central cantilever axis of the dam. This is a logical result
because Karun1 dam has symmetrical plan and the loads applied to it including
thermal forces are also symmetric to the mentioned axis. Vertical contraction
joints have not meaningful effect on the values of arch stresses in Karun1
dam. This can be shown through the examination of Fig. 8.
This is a rational conclusion because the joints are orthogonal to the arch
direction and under effect of thermostatic loads the arches would be compressed.
So, theses joints have not valuable effects on arch stresses. Another interesting
result obtained from that picture could be the arch stress concentrations near
the abutments in the models which consider the heat transfer e.g., model A2
and B2. This can be interpreted in this way: when heat transfer is applied to
the dam, the nonhomogenity which exists between the dam and abutment in thermal
conduction properties cause to the concentration of arch stresses in vicinity
the damabutment line. This stress concentration also occurs for cantilever
stresses due to the same reason (Fig. 9). As it can be seen
from both Fig. 8 and 9, the arch and cantilever
stresses along the crest and central cantilever respectively all have compressive
nature. After surveying the principal stresses inside the dam (Table
1), it was obvious that under the effect of both contraction joints and
thermal forces, the tensile stresses have been increased in dam body whereas
near the abutments the mentioned stresses have been decreased. Not that the
negative values in Table 1 present the compressive stresses
whereas the positives reflect the tensile stresses. So, it is necessary that
for having more reliable assessment of tensile cracking potential in Karun1
dam body, the heat transfer and contraction joints are being simultaneously
considered in the stability analysis conjunction with hydrostatic and gravity
loads.
Table 1: 
Comparison between the maximum and minimum stresses 

With respect to the results obtained in this study, the contribution of thermal
loads in displacements and stresses is as the contribution of hydrostatic plus
gravity loads. This result is in contrast to the conclusion which obtained by
the Sheibani and Ghaemian (2006) for the Karaj arch
dam.
CONCLUSION
The model which has been adopted via using ABAQUS package for Karun1 dam has good potential for considering two important factor in stability evaluation of the mentioned dam; Contraction vertical joints and heat transfer in summer season. Because the resulted displacements from the model are coincidence with the values read from pendulum instruments located in central block of the dam. If one of these two factors is neglected in displacement or stress analysis, the significant errors would result in dam safety qualification.