where푢 1 and푢 2 are the displacement increments inΩ 1 and
Ω 2 whose boundaries exclude the exact contact boundary,
푢 12 is the displacement increment alongΓ1푐,푢 21 is the
displacement increment alongΓ2푐,and푃is the interaction
forcealongthecontactinterfaceΓ푐.Itcanbeprovedthat푃
is equivalent to the Lagrange multiplier [ 17 , 18 ].
When the two bodies are not in contact, one body
imposes no constraints on the other, and thus ( 13 )and( 14 )
are independent of each other and푃≡0. The displacement
increments푢 1 and푢 12 aresolvedusing( 13 ), while푢 2 and푢 21
are determined by ( 14 ).
When the two bodies are in contact, one deformable
body imposes constraints on the other. At this time,푢 12 and
푢 21 are no longer independent, and푃is introduced as an
unknown. The contact boundary should satisfy the kinematic
and dynamic constraints. As shown in Figure 4 ,ifpointAon
Γ1푐coincides with point B onΓ2푐, the kinematic constraint is
expressed as [ 12 ]
(푢A 12 −푢B 21 )휂≤TOL, (15)
where휂is the directional cosine at the contact point and
TOL is the closure distance or contact tolerance. The dynamic
condition is Coulomb’s friction law in our computation:
푃푡≤−휇⋅푃푛⋅휂푡, (16)where푃푡is the tangential friction traction force,푃푛is the
normal traction force,휇is the friction coefficient, and휂푡
isthetangentialvectorinthedirectionofrelativevelocity.
Therefore, the unknowns푢 1 ,푢 12 ,푢 2 ,푢 21 ,푃,andΓ푐can be
completely solved from ( 13 )–( 16 ).
2.4.2. Strategy of Searching Contact Points.The contact inter-
faceΓ푐is the key unknown in the contact problem. Zhang et
al. [ 12 ]usedatypicalnode-edgecontactmodetoimplement
contact detection at the element level. The disadvantage of
thisnode-edgecontactmodeisthattheaccuracyislow.
This study uses curve fitting; that is, the point interpolation
method [ 14 , 15 ], to detect the exact contact interfaceΓ푐.The
numerical procedure is as follows.
Step 1.Assume potential contact interfacesΓ1푐onΩ 1 andΓ2푐
onΩ 2.
Step 2.Locate the nodal points on the interfaces Γ1푐
andΓ2푐.Thereare푀nodes onΓ1푐, denoted by푥 11 ,푥 12 ,
...,푥1푖,...,푥1푀and푁nodes onΓ2푐, denoted by푥 21 ,푥 22 ,
...,푥2푖,...,푥2푁.
Step 3.Interpolate these nodes to form the boundary lines
using the radial point interpolation method [ 14 , 15 ].
One has
푥=
푀
∑
푖=1푁1푖푥1푖,푥=
푁
∑
푗=1푁2푗푥2푗, (17)
where the shape functions푁1푖,푁2푗are determined using
point interpolation methods [ 14 , 15 ].
Step 4.Establish the distance function훿along either bound-
ary line. The point is not in contact when훿>TOL. Other-
wise, the point is in contact with the other boundary. Identify
the exact contact points through ( 17 ). Iterate the same proce-
duretofindouttheentirecontactboundary.Step 5.Iterate FEM computation to satisfy the equilibrium of
twodeformablebodiesandthecontactboundaryconditions.Step 6.Update nodal coordinates on the contact boundary.
Carry out the next step computation, and return to Step 3 for
the same search procedure for the contact points.3. Constitutive Models for Dam Materials
3.1. EB Model for Rockfill Materials.Rockfill materials and
soil masses behave with strong nonlinearity because of the
high stress levels in dams. This nonlinearity is described by
the following incremental Hooke’s law:{
{
{
푑휎푥
푑휎푦
푑휏푥푦
}
}
}
=[퐷]
{
{
{
푑휀푥
푑휀푦
푑훾푥푦
}
}
}
=
[
[[
[
[
[
[
퐵푡+
4
3
퐺푡 퐵푡−
2
3
퐺푡 0
퐵푡−
2
3
퐺푡 퐵푡+
4
3
퐺푡 0
002퐺푡
]
]]
]
]
]
]
{
{
{
푑휀푥
푑휀푦
푑훾푥푦
}
}
}
,
(18)
where퐵푡is the bulk modulus,퐺푡 =3퐵푡퐸푡/(9퐵푡−퐸푡)is
the shear modulus, and퐸푡is the deformation modulus. The
Duncan EB model [ 16 ] gives the deformation modulus퐸푡as
follows:퐸푡=푘⋅푃푎(
휎 3
푃푎
)
푛
[1−푅푓(1 −sin휙)(휎 1 −휎 3 )
2푐 ⋅cos휙+2휎 3 sin휙]
2
, (19)where(휎 1 −휎 3 )is the deviatoric stress,휎 3 is the confining
pressure,푐is the cohesion intercept,휙is the angle of internal
friction,푅푓is the failure ratio,푃푎is the atmospheric pressure,
and푘and푛are constants.
In the computation, the rockfill material has푐=0and a
variable angle of internal friction휙휙=휙 0 −Δ휙log(휎 3
푃푎
), (20)
where휙 0 andΔ휙are two constants. Another parameter, bulk
modulus퐵푡, is assumed to be퐵푡=푘푏푃푎(
휎 3
푃푎
)
푚
, (21)where푘푏and푚are constants.3.2.LinearElasticModelfortheConcreteFaceSlab.Alinear
elastic model with Young’s modulus퐸and Poisson ratio]is
used to describe the mechanical properties of the concrete
face slab. No failure is allowed.