# Comparative Study on Interface Elements, Thin-Layer Elements, and Contact Analysis Methods in the Analysis of High Concrete-Faced Rockfill Dams

Authors: Xiao-xiang Qian; Hui-na Yuan; Quan-ming Li; Bing-yin Zhang
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

This paper presents a study on the numerical performance of three contact simulation methods, namely, the interface element, thin-layer element, and contact analysis methods, through the analysis of the contact behavior between the concrete face slab and the dam body of a high concrete-faced rockfill dam named Tianshengqiao-I in China. To investigate the accuracy and limitations of each method, the simulation results are compared in terms of the dam deformation, contact stress along the interface, stresses in the concrete face slab, and separation of the concrete face slab from the cushion layer. In particular, the predicted dam deformation and slab separation are compared with the in-situ observation data to classify these methods according to their agreement with the in-situ observations. It is revealed that the interface element and thin-layer element methods have their limitations in predicting contact stress, slab separation, and stresses in the concrete face slab if a large slip occurs. The contact analysis method seems to be the best choice whether the separation is finite or not.


---

## Body

## 1. Introduction

The cracking of the concrete slab is the most important factor affecting the safety of concrete-faced rockfill dams (CFRDs). Accurate computation of stress and deformation in the concrete slab are key issues for slab cracking assessment. Numerical methods can be used to predict the deformation and stress distributions in the concrete face slab, where the behavior of the interface between the concrete face slab and the cushion layer plays a significant role. Because the interface can be treated in different ways, the prediction of displacement and stress distribution around the interface may be different. This study focuses on the comparison of different interface analysis methods through the analysis of stress and displacement distributions near the interface in the Tianshengqiao-I CFRD project. The accuracy and limitations of each method are discussed.Much attention has been paid to numerical treatment of the interfaces in geotechnical problems such as buried structures, jointed rocks, and rockfill dams [1–5]. Interface behavior often involves large relative movement or even debonding [6]. Over the past three decades, three numerical methods have been proposed for simulating the displacement jump along the interface: the interface element, thin-layer element, and contact analysis methods. The interface element method originated from the Goodman joint element approach [2–6]. The basic idea was to introduce a constitutive model for an interface of zero thickness [6]. This constitutive model may be elastic, rigid-plastic, or elastic-plastic [2, 6, 7]. As an alternative, a thin-layer element method [8] was proposed. The thin-layer element method regards joints or interfaces as conventional continuums described by solid elements. However, the material modulus for this thin layer is much lower than that for the intact solid [8–11]. This thin-layer element method has been successfully applied to jointed rock masses [10], buried pipes [8], and the interaction of foundation and soil masses [9, 11]. Either the interface element or thin-layer element is limited to small deformation. Different from the previous two numerical methods, the contact analysis method was proposed to simulate the contact behaviors between the concrete face slab and the cushion layer in the Tianshengqiao-I concrete-faced rockfill dam [12]. In this contact analysis method, the concrete face slab and dam body were regarded as two independent deformable bodies, and the contact interface was treated using contact mechanics [13]. This method allows large relative displacements between the concrete face slab and cushion layer. The physical and mechanical properties of the interface can also be nonlinear or elastic-plastic. In the contact analysis method, the detection of the contact is the key issue. Zhang et al. [12] proposed a local contact detection method at the element level, where the search is localized between two elements and thus needs less time. However, the accuracy of this contact detection method is not acceptable when the mapping function for element geometry is not identical to that for displacement interpolation and when the deformation is large. In this paper, a global contact search method is proposed based on a radial point interpolation method [14, 15]. The accuracy of this global search method is controllable.In this study, the numerical performance of three numerical simulation methods, namely, the interface element, thin-layer element, and contact analysis methods is compared through stress-deformation analysis of a high concrete-faced rockfill dam. In Section2, the fundamentals of the three methods are briefly reviewed. A global search method for contact detection is proposed based on the radial point interpolation method. In Section 3, the constitutive models for the rockfill dam body and the concrete face slab are presented. The Duncan EB model [16] is employed to describe the nonlinearity of rockfill materials, and a linear elastic model is used to describe the mechanical properties of the concrete face slab. In Section 4, the FEM models and material parameters are introduced. Section 5 compares the performance of the three numerical methods using the Tianshengqiao-I CFRD project in China as an example. The separation between the concrete face slab and the cushion layer, stresses in the concrete face slab, contact stress along the interface, displacements along the interface, and deformation of the dam body are compared using the in-situ observations available. Finally, conclusions are drawn in Section 6.


## 2. Fundamentals of Numerical Methods for the Interfaces

### 2.1. The Contact Problem

With reference to Figure1, we consider the contact of two deformable bodies, where the problem domain Ω is divided into two subdomains Ω1 (bounded by Γ1) and Ω2 (bounded by Γ2). The bodies are fixed at Γu=Γ1u∪Γ2u and subjected to boundary traction t at Γt=Γ1t∪Γ2t. Γ1c and Γ2c are the potential contacting boundaries of Ω1 and Ω2, respectively, while Γc denotes the exact contact part on Γ1c and Γ2c.Figure 1
Contact of two deformable bodies.


### 2.2. Interface Element Method

For the interface element method (Figure2), the interface conditions are described by
(1)σ·n|Γ1c∩Γc=σ·n|Γ2c∩Γc[δu]≥0,
where σ is the stress tensor, n is the outward normal, and [δu] denotes the increment of a displacement jump [2, 6]. Such a problem has the following weak form:
(2){∫Ω1{δε}T{σ}dΩ-∫Ω1{δu}T{b}dΩ-∫Γ1t{δu}T{t}dΓ}+{∫Ω2{δε}T{σ}dΩ-∫Ω2{δu}T{b}dΩ-∫Γ2t{δu}T{t}dΓ}+∫Γc{σ}[δu]dΓ=0,
where ε is the strain tensor, u is the displacement, b is the body force, and t is the boundary traction. This weak form is composed of three terms:
(3)π1+π2+πinterface=0,
where π1 denotes the terms in the first bracket to express the potential in Ω1, π2 denotes the terms in the second bracket to express the potential in Ω2, and πinterface denotes the last term to express the potential along the interface Γc. On discretizing the interface term πinterface, the element stiffness is obtained as
(4)Kine=∫SeTTNuT[D-]epNuTdS,
where T and Nu are the transformation matrix and shape function of the interface element Se. The material matrix [D-]ep is defined using the following constitutive law of an interface [6]:
(5){ΔσnΔτ}=[D-]ep{[Δun][Δus]},[D-]ep=[knknsksnks],
where σn, τ are the normal and shear stresses, un, us are the normal and shear displacements, kn, ks are the normal and shear stiffness, and ksn, kns are the coupling stiffnesses between normal and shear deformations.Figure 2
Interface element method.Goodman et al. [2] did not consider the coupling effect between normal and shear deformations. They took the material matrix as
(6)[D-]ep=[kn00ks]
and the shear stiffness ks as
(7)ks=k1γw(σnPa)n1(1-Rf1τσntgϕ)2,
where k1 and n1 are two parameters, σn is the normal stress on the interface, τ is the shear stress along the interface, Pa is the atmospheric pressure, γw is the unit weight of water, Rf1 is the failure ratio, and ϕ is the angle of internal friction. k1, n1, Rf1, and ϕ are the four parameters to be determined from direct shear tests. The normal stiffness kn is usually given a large number when the interface element is in compression and a small number when in tension.


### 2.3. Thin-Layer Element Method

In this method, an interface is treated as a thin-layer solid element (Figure3). This thin layer is given a relatively low modulus and can experience large deformation [8–11]. The problem shown in Figure 1 with a thin layer has the following weak form:
(8)π1+π2+πthin=0,
where the term on the thin layer, πthin, is given by
(9)πthin=∫VL{δε}T{σ}dV,
with VL denoting the domain of interface Γc. If VL has a finite thickness of d, the element stiffness of thin layer element Ve is
(10)Kthe=∫VeBTDBdV≅d≪Sed∫SeBTDBdS,
where B is the strain matrix, D is the material matrix, and Se is the element length. Previous studies revealed that the accuracy of element stiffness is sensitive to the aspect ratio d/Se. When the aspect ratio varies in the range of 0.01–0.1, slippage is modeled quite accurately [8–11].Figure 3
Thin-layer element method.


### 2.4. Contact Analysis Method

#### 2.4.1. Contact of Two Deformable Bodies

As shown in Figure4, the potential contact boundaries are Γ1c in Γ1 and Γ2c in Γ2, while the exact contact boundary is denoted as interface Γc, which is usually unknown beforehand. The weak form of each deformable body is expressed individually as follows.Figure 4
Contact analysis method.For deformable bodyΩ1(11){∫Ω1{δε}T{σ}dΩ-∫Ω1{δu}T{b}dΩ-∫Γ1t{δu}T{t}dΓ}-∫Γc{δu}T{P}dΓ=0.For deformable bodyΩ2(12){∫Ω2{δε}T{σ}dΩ-∫Ω2{δu}T{b}dΩ-∫Γ2t{δu}T{t}dΓ}-∫Γc{δu}T{P}dΓ=0,
where {P} is the interaction force. Upon discretizing the weak forms in (11) and (12), the following discrete system equation is obtained for each deformable body:
(13)K11u1+K12u12+L1P=f1forΩ1,(14)K22u2+K21u21+L2P=f2forΩ2,
where u1 and u2 are the displacement increments in Ω1 and Ω2 whose boundaries exclude the exact contact boundary, u12 is the displacement increment along Γ1c, u21 is the displacement increment along Γ2c, and P is the interaction force along the contact interface Γc. It can be proved that P 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 P≡0. The displacement increments u1 and u12 are solved using (13), while u2 and u21 are determined by (14).When the two bodies are in contact, one deformable body imposes constraints on the other. At this time,u12 and u21 are no longer independent, and P is introduced as an unknown. The contact boundary should satisfy the kinematic and dynamic constraints. As shown in Figure 4, if point A on Γ1c coincides with point B on Γ2c, the kinematic constraint is expressed as [12]
(15)(u12A-u21B)η≤TOL,
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)Pt≤-μ·Pn·ηt,
where Pt is the tangential friction traction force, Pn is the normal traction force, μ is the friction coefficient, and ηt is the tangential vector in the direction of relative velocity. Therefore, the unknowns u1, u12, u2, u21, P, and Γc can be completely solved from (13)–(16).


#### 2.4.2. Strategy of Searching Contact Points

The contact interfaceΓc is the key unknown in the contact problem. Zhang et al. [12] used a typical node-edge contact mode to implement contact detection at the element level. The disadvantage of this node-edge contact mode is that the accuracy is low. This study uses curve fitting; that is, the point interpolation method [14, 15], to detect the exact contact interface Γc. The numerical procedure is as follows.Step 1.
Assume potential contact interfacesΓ1c on Ω1 and Γ2c on Ω2.Step 2.
Locate the nodal points on the interfacesΓ1c and Γ2c. There are M nodes on Γ1c, denoted by x11,x12,…,x1i,…,x1M and N nodes on Γ2c, denoted by x21,x22,…,x2i,…,x2N.Step 3.
Interpolate these nodes to form the boundary lines using the radial point interpolation method [14, 15].
One has(17)x=∑i=1MN1ix1i,x=∑j=1NN2jx2j,
where the shape functions N1i, N2j are determined using point interpolation methods [14, 15].Step 4.
Establish the distance functionδ along either boundary line. The point is not in contact when δ>TOL. Otherwise, the point is in contact with the other boundary. Identify the exact contact points through (17). Iterate the same procedure to find out the entire contact boundary.Step 5.
Iterate FEM computation to satisfy the equilibrium of two deformable bodies and the contact boundary conditions.Step 6.
Update nodal coordinates on the contact boundary. Carry out the next step computation, and return to Step3 for the same search procedure for the contact points.


## 2.1. The Contact Problem

With reference to Figure1, we consider the contact of two deformable bodies, where the problem domain Ω is divided into two subdomains Ω1 (bounded by Γ1) and Ω2 (bounded by Γ2). The bodies are fixed at Γu=Γ1u∪Γ2u and subjected to boundary traction t at Γt=Γ1t∪Γ2t. Γ1c and Γ2c are the potential contacting boundaries of Ω1 and Ω2, respectively, while Γc denotes the exact contact part on Γ1c and Γ2c.Figure 1
Contact of two deformable bodies.


## 2.2. Interface Element Method

For the interface element method (Figure2), the interface conditions are described by
(1)σ·n|Γ1c∩Γc=σ·n|Γ2c∩Γc[δu]≥0,
where σ is the stress tensor, n is the outward normal, and [δu] denotes the increment of a displacement jump [2, 6]. Such a problem has the following weak form:
(2){∫Ω1{δε}T{σ}dΩ-∫Ω1{δu}T{b}dΩ-∫Γ1t{δu}T{t}dΓ}+{∫Ω2{δε}T{σ}dΩ-∫Ω2{δu}T{b}dΩ-∫Γ2t{δu}T{t}dΓ}+∫Γc{σ}[δu]dΓ=0,
where ε is the strain tensor, u is the displacement, b is the body force, and t is the boundary traction. This weak form is composed of three terms:
(3)π1+π2+πinterface=0,
where π1 denotes the terms in the first bracket to express the potential in Ω1, π2 denotes the terms in the second bracket to express the potential in Ω2, and πinterface denotes the last term to express the potential along the interface Γc. On discretizing the interface term πinterface, the element stiffness is obtained as
(4)Kine=∫SeTTNuT[D-]epNuTdS,
where T and Nu are the transformation matrix and shape function of the interface element Se. The material matrix [D-]ep is defined using the following constitutive law of an interface [6]:
(5){ΔσnΔτ}=[D-]ep{[Δun][Δus]},[D-]ep=[knknsksnks],
where σn, τ are the normal and shear stresses, un, us are the normal and shear displacements, kn, ks are the normal and shear stiffness, and ksn, kns are the coupling stiffnesses between normal and shear deformations.Figure 2
Interface element method.Goodman et al. [2] did not consider the coupling effect between normal and shear deformations. They took the material matrix as
(6)[D-]ep=[kn00ks]
and the shear stiffness ks as
(7)ks=k1γw(σnPa)n1(1-Rf1τσntgϕ)2,
where k1 and n1 are two parameters, σn is the normal stress on the interface, τ is the shear stress along the interface, Pa is the atmospheric pressure, γw is the unit weight of water, Rf1 is the failure ratio, and ϕ is the angle of internal friction. k1, n1, Rf1, and ϕ are the four parameters to be determined from direct shear tests. The normal stiffness kn is usually given a large number when the interface element is in compression and a small number when in tension.


## 2.3. Thin-Layer Element Method

In this method, an interface is treated as a thin-layer solid element (Figure3). This thin layer is given a relatively low modulus and can experience large deformation [8–11]. The problem shown in Figure 1 with a thin layer has the following weak form:
(8)π1+π2+πthin=0,
where the term on the thin layer, πthin, is given by
(9)πthin=∫VL{δε}T{σ}dV,
with VL denoting the domain of interface Γc. If VL has a finite thickness of d, the element stiffness of thin layer element Ve is
(10)Kthe=∫VeBTDBdV≅d≪Sed∫SeBTDBdS,
where B is the strain matrix, D is the material matrix, and Se is the element length. Previous studies revealed that the accuracy of element stiffness is sensitive to the aspect ratio d/Se. When the aspect ratio varies in the range of 0.01–0.1, slippage is modeled quite accurately [8–11].Figure 3
Thin-layer element method.


## 2.4. Contact Analysis Method

### 2.4.1. Contact of Two Deformable Bodies

As shown in Figure4, the potential contact boundaries are Γ1c in Γ1 and Γ2c in Γ2, while the exact contact boundary is denoted as interface Γc, which is usually unknown beforehand. The weak form of each deformable body is expressed individually as follows.Figure 4
Contact analysis method.For deformable bodyΩ1(11){∫Ω1{δε}T{σ}dΩ-∫Ω1{δu}T{b}dΩ-∫Γ1t{δu}T{t}dΓ}-∫Γc{δu}T{P}dΓ=0.For deformable bodyΩ2(12){∫Ω2{δε}T{σ}dΩ-∫Ω2{δu}T{b}dΩ-∫Γ2t{δu}T{t}dΓ}-∫Γc{δu}T{P}dΓ=0,
where {P} is the interaction force. Upon discretizing the weak forms in (11) and (12), the following discrete system equation is obtained for each deformable body:
(13)K11u1+K12u12+L1P=f1forΩ1,(14)K22u2+K21u21+L2P=f2forΩ2,
where u1 and u2 are the displacement increments in Ω1 and Ω2 whose boundaries exclude the exact contact boundary, u12 is the displacement increment along Γ1c, u21 is the displacement increment along Γ2c, and P is the interaction force along the contact interface Γc. It can be proved that P 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 P≡0. The displacement increments u1 and u12 are solved using (13), while u2 and u21 are determined by (14).When the two bodies are in contact, one deformable body imposes constraints on the other. At this time,u12 and u21 are no longer independent, and P is introduced as an unknown. The contact boundary should satisfy the kinematic and dynamic constraints. As shown in Figure 4, if point A on Γ1c coincides with point B on Γ2c, the kinematic constraint is expressed as [12]
(15)(u12A-u21B)η≤TOL,
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)Pt≤-μ·Pn·ηt,
where Pt is the tangential friction traction force, Pn is the normal traction force, μ is the friction coefficient, and ηt is the tangential vector in the direction of relative velocity. Therefore, the unknowns u1, u12, u2, u21, P, and Γc can be completely solved from (13)–(16).


### 2.4.2. Strategy of Searching Contact Points

The contact interfaceΓc is the key unknown in the contact problem. Zhang et al. [12] used a typical node-edge contact mode to implement contact detection at the element level. The disadvantage of this node-edge contact mode is that the accuracy is low. This study uses curve fitting; that is, the point interpolation method [14, 15], to detect the exact contact interface Γc. The numerical procedure is as follows.Step 1.
Assume potential contact interfacesΓ1c on Ω1 and Γ2c on Ω2.Step 2.
Locate the nodal points on the interfacesΓ1c and Γ2c. There are M nodes on Γ1c, denoted by x11,x12,…,x1i,…,x1M and N nodes on Γ2c, denoted by x21,x22,…,x2i,…,x2N.Step 3.
Interpolate these nodes to form the boundary lines using the radial point interpolation method [14, 15].
One has(17)x=∑i=1MN1ix1i,x=∑j=1NN2jx2j,
where the shape functions N1i, N2j are determined using point interpolation methods [14, 15].Step 4.
Establish the distance functionδ along either boundary line. The point is not in contact when δ>TOL. Otherwise, the point is in contact with the other boundary. Identify the exact contact points through (17). Iterate the same procedure to find out the entire contact boundary.Step 5.
Iterate FEM computation to satisfy the equilibrium of two deformable bodies and the contact boundary conditions.Step 6.
Update nodal coordinates on the contact boundary. Carry out the next step computation, and return to Step3 for the same search procedure for the contact points.


## 2.4.1. Contact of Two Deformable Bodies

As shown in Figure4, the potential contact boundaries are Γ1c in Γ1 and Γ2c in Γ2, while the exact contact boundary is denoted as interface Γc, which is usually unknown beforehand. The weak form of each deformable body is expressed individually as follows.Figure 4
Contact analysis method.For deformable bodyΩ1(11){∫Ω1{δε}T{σ}dΩ-∫Ω1{δu}T{b}dΩ-∫Γ1t{δu}T{t}dΓ}-∫Γc{δu}T{P}dΓ=0.For deformable bodyΩ2(12){∫Ω2{δε}T{σ}dΩ-∫Ω2{δu}T{b}dΩ-∫Γ2t{δu}T{t}dΓ}-∫Γc{δu}T{P}dΓ=0,
where {P} is the interaction force. Upon discretizing the weak forms in (11) and (12), the following discrete system equation is obtained for each deformable body:
(13)K11u1+K12u12+L1P=f1forΩ1,(14)K22u2+K21u21+L2P=f2forΩ2,
where u1 and u2 are the displacement increments in Ω1 and Ω2 whose boundaries exclude the exact contact boundary, u12 is the displacement increment along Γ1c, u21 is the displacement increment along Γ2c, and P is the interaction force along the contact interface Γc. It can be proved that P 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 P≡0. The displacement increments u1 and u12 are solved using (13), while u2 and u21 are determined by (14).When the two bodies are in contact, one deformable body imposes constraints on the other. At this time,u12 and u21 are no longer independent, and P is introduced as an unknown. The contact boundary should satisfy the kinematic and dynamic constraints. As shown in Figure 4, if point A on Γ1c coincides with point B on Γ2c, the kinematic constraint is expressed as [12]
(15)(u12A-u21B)η≤TOL,
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)Pt≤-μ·Pn·ηt,
where Pt is the tangential friction traction force, Pn is the normal traction force, μ is the friction coefficient, and ηt is the tangential vector in the direction of relative velocity. Therefore, the unknowns u1, u12, u2, u21, P, and Γc can be completely solved from (13)–(16).


## 2.4.2. Strategy of Searching Contact Points

The contact interfaceΓc is the key unknown in the contact problem. Zhang et al. [12] used a typical node-edge contact mode to implement contact detection at the element level. The disadvantage of this node-edge contact mode is that the accuracy is low. This study uses curve fitting; that is, the point interpolation method [14, 15], to detect the exact contact interface Γc. The numerical procedure is as follows.Step 1.
Assume potential contact interfacesΓ1c on Ω1 and Γ2c on Ω2.Step 2.
Locate the nodal points on the interfacesΓ1c and Γ2c. There are M nodes on Γ1c, denoted by x11,x12,…,x1i,…,x1M and N nodes on Γ2c, denoted by x21,x22,…,x2i,…,x2N.Step 3.
Interpolate these nodes to form the boundary lines using the radial point interpolation method [14, 15].
One has(17)x=∑i=1MN1ix1i,x=∑j=1NN2jx2j,
where the shape functions N1i, N2j are determined using point interpolation methods [14, 15].Step 4.
Establish the distance functionδ along either boundary line. The point is not in contact when δ>TOL. Otherwise, the point is in contact with the other boundary. Identify the exact contact points through (17). Iterate the same procedure to find out the entire contact boundary.Step 5.
Iterate FEM computation to satisfy the equilibrium of two deformable bodies and the contact boundary conditions.Step 6.
Update nodal coordinates on the contact boundary. Carry out the next step computation, and return to Step3 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:(18){dσxdσydτxy}=[D]{dεxdεydγxy}=[Bt+43GtBt-23Gt0Bt-23GtBt+43Gt0002Gt]{dεxdεydγxy},
where Bt is the bulk modulus, Gt=3BtEt/(9Bt-Et) is the shear modulus, and Et is the deformation modulus. The Duncan EB model [16] gives the deformation modulus Et as follows:
(19)Et=k·Pa(σ3Pa)n[1-Rf(1-sinϕ)(σ1-σ3)2c·cosϕ+2σ3sinϕ]2,
where (σ1-σ3) is the deviatoric stress, σ3 is the confining pressure, c is the cohesion intercept, ϕ is the angle of internal friction, Rf is the failure ratio, Pa is the atmospheric pressure, and k and n are constants.In the computation, the rockfill material hasc=0 and a variable angle of internal friction ϕ(20)ϕ=ϕ0-Δϕlog(σ3Pa),
where ϕ0 and Δϕ are two constants. Another parameter, bulk modulus Bt, is assumed to be
(21)Bt=kbPa(σ3Pa)m,
where kb and m are constants.


### 3.2. Linear Elastic Model for the Concrete Face Slab

A linear elastic model with Young’s modulusE and Poisson ratio ν is used to describe the mechanical properties of the concrete face slab. No failure is allowed.


## 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:(18){dσxdσydτxy}=[D]{dεxdεydγxy}=[Bt+43GtBt-23Gt0Bt-23GtBt+43Gt0002Gt]{dεxdεydγxy},
where Bt is the bulk modulus, Gt=3BtEt/(9Bt-Et) is the shear modulus, and Et is the deformation modulus. The Duncan EB model [16] gives the deformation modulus Et as follows:
(19)Et=k·Pa(σ3Pa)n[1-Rf(1-sinϕ)(σ1-σ3)2c·cosϕ+2σ3sinϕ]2,
where (σ1-σ3) is the deviatoric stress, σ3 is the confining pressure, c is the cohesion intercept, ϕ is the angle of internal friction, Rf is the failure ratio, Pa is the atmospheric pressure, and k and n are constants.In the computation, the rockfill material hasc=0 and a variable angle of internal friction ϕ(20)ϕ=ϕ0-Δϕlog(σ3Pa),
where ϕ0 and Δϕ are two constants. Another parameter, bulk modulus Bt, is assumed to be
(21)Bt=kbPa(σ3Pa)m,
where kb and m are constants.


## 3.2. Linear Elastic Model for the Concrete Face Slab

A linear elastic model with Young’s modulusE and Poisson ratio ν is used to describe the mechanical properties of the concrete face slab. No failure is allowed.


## 4. Computation Models and Parameters

### 4.1. Tianshengqiao-I Concrete-Faced Rockfill Dam Project

The Tianshengqiao-I hydropower project is on the Nanpan River in southwestern China [12]. Its water retaining structure is a concrete-faced rockfill dam, 178 m high and 1104 m long. The rockfill volume of the dam body is about 18 million m3, and the area of the concrete face is 173,000 m2. A surface chute spillway on the right bank allows a maximum discharge of 19,450 m3/s. The tunnel in the right abutment is used for emptying the reservoir during operation. The left abutment has four power tunnels and a surface powerhouse with a total capacity of 1,200 MW. Material zoning and construction stages are shown in Figure 5. The design parameters of the dam materials are listed in Table 1, and the details of each construction stage are given in Table 2.Table 1
Design parameters of dam materials.
Mat. number
Mat. description
Max. particle size (cm)
Dry unit weight (KN/m3)
Void ratio (%)
IIA
Processed limestone
8
22.0
19
IIIA
Limestone
30
21.5
21
IIIB
Limestone
80
21.2
22
IIIC
Mudstone and sandstone
80
21.5
22
IIID
Limestone
160
20.5
24Table 2
Construction stages and time.
Filling step
Time
Remark
①
1996.01–1996.06
Fill dam body
②
1996.07–1997.02
Fill dam body
③
1997.03–1997.05
Cast Phase 1 concrete slab
④ and ⑥
1997.02–1997.10
Fill dam body
⑤
1997.05-1997.05
Water level rises
⑦
1997.06–1997.10
Water level fluctuation
⑧ and ⑨
1997.11–1998.01
Fill dam body
⑩
1997.12–1998.05
Cast Phase 2 concrete slab
⑪
1997.11-1997.12
Water level rises
⑫
1998.02–1998.08
Fill dam body
⑬
1998.06-1998.07
Water level rises
⑭
1998.08–1999.01
Fill dam body
⑮
1999.01–1999.05
Cast Phase 3 concrete slab
⑯
1999.06–1999.09
Store waterMaterial zones and construction stages of Tianshengqiao-I CFRD.
(a)
Material zoning
(b)
Construction stages


### 4.2. Computation Section, Procedure, and Material Parameters

A two-dimensional finite element analysis was performed [19]. The maximum cross-section (section 0+630 m), which is in the middle of the riverbed, was taken for computation. Figure 6(a) shows the finite element mesh for the contact analysis method. It has a total of 402 four-node elements in the dam body and 46 four-node elements in the concrete face slab (the concrete face slab is divided into two layers of elements). The mesh for the interface element method is shown in Figure 6(b), where a row of interface elements is placed along the interface between the concrete face slab and the cushion layer. This mesh model has 23 additional interface elements compared to the mesh for the contact analysis model. If the interface elements in Figure 6(b) are assigned a thickness of 0.3 m, the finite element mesh for the thin-layer element method is obtained. Because the length of each element is 12 m, the thin-layer elements have an aspect ratio of 0.025, in the range of 0.01–0.1 [8–11]. The previous mesh models show that the dam body and concrete face slab can be meshed independently for the contact analysis method. This may produce nonmatching nodes on both sides of the interface [15]. However, the thin-layer element and interface element methods usually require matching nodes on both sides of the interface. This model sets zero displacements along the rock base [12].Two-dimensional finite element mesh.
(a)
For the contact analysis method
(b)
For the interface element methodThe computational procedure follows exactly the construction stages shown in Figure5(b). First, blocks ① and ② of the dam body were built up to El.682 m. In each block, layer-by-layer elements were activated to simulate the construction process, and the midpoint stiffness [20] was used for the nonlinear constitutive model. Before placement of the freshly cast stage-I slab, the calculated displacements of the dam body were set to zero, and the calculated stresses were retained. The elements of stage-I slab ③ were then activated, and dam construction continued. The impounding process was simulated by increasing the water level by 10 m in each increment. The same procedure was repeated until completion of the whole dam body.The concrete face slab had an elastic modulus of3×104 MPa and a Poisson’s ratio of 0.2. Table 3 gives the computational parameters of the rockfill materials for the EB model. An elastic modulus of 6 MPa and Poisson’s ratio of 0.2 were used for the materials in the thin-layer elements. The computational parameters for the Goodman interface model are listed in Table 4.Table 3
Computational parameters for the rockfill materials.
Mat. number
Density (kg/m3)
ϕ
′  [°]
Δ
ϕ  [°]
k
n
R
f
k
b
m
IIA
2200
50.6
7.0
1000
0.35
0.71
450
0.24
IIIA
2100
52.5
8.0
900
0.36
0.76
400
0.19
IIIB
2100
51.0
13.0
564
0.35
0.85
204
0.18
IIID
2050
51.0
13.5
432
0.30
0.80
300
−0.18
IIIC
2150
45.0
10.0
250
0.25
0.73
125
0.00Table 4
Parameters of the Goodman interface model.
ϕ  [°]
k
1
n
1
R
f
1
k
n (MPa)
Compression
Tension
30
1000
0.3
1
10000
1


## 4.1. Tianshengqiao-I Concrete-Faced Rockfill Dam Project

The Tianshengqiao-I hydropower project is on the Nanpan River in southwestern China [12]. Its water retaining structure is a concrete-faced rockfill dam, 178 m high and 1104 m long. The rockfill volume of the dam body is about 18 million m3, and the area of the concrete face is 173,000 m2. A surface chute spillway on the right bank allows a maximum discharge of 19,450 m3/s. The tunnel in the right abutment is used for emptying the reservoir during operation. The left abutment has four power tunnels and a surface powerhouse with a total capacity of 1,200 MW. Material zoning and construction stages are shown in Figure 5. The design parameters of the dam materials are listed in Table 1, and the details of each construction stage are given in Table 2.Table 1
Design parameters of dam materials.
Mat. number
Mat. description
Max. particle size (cm)
Dry unit weight (KN/m3)
Void ratio (%)
IIA
Processed limestone
8
22.0
19
IIIA
Limestone
30
21.5
21
IIIB
Limestone
80
21.2
22
IIIC
Mudstone and sandstone
80
21.5
22
IIID
Limestone
160
20.5
24Table 2
Construction stages and time.
Filling step
Time
Remark
①
1996.01–1996.06
Fill dam body
②
1996.07–1997.02
Fill dam body
③
1997.03–1997.05
Cast Phase 1 concrete slab
④ and ⑥
1997.02–1997.10
Fill dam body
⑤
1997.05-1997.05
Water level rises
⑦
1997.06–1997.10
Water level fluctuation
⑧ and ⑨
1997.11–1998.01
Fill dam body
⑩
1997.12–1998.05
Cast Phase 2 concrete slab
⑪
1997.11-1997.12
Water level rises
⑫
1998.02–1998.08
Fill dam body
⑬
1998.06-1998.07
Water level rises
⑭
1998.08–1999.01
Fill dam body
⑮
1999.01–1999.05
Cast Phase 3 concrete slab
⑯
1999.06–1999.09
Store waterMaterial zones and construction stages of Tianshengqiao-I CFRD.
(a)
Material zoning
(b)
Construction stages


## 4.2. Computation Section, Procedure, and Material Parameters

A two-dimensional finite element analysis was performed [19]. The maximum cross-section (section 0+630 m), which is in the middle of the riverbed, was taken for computation. Figure 6(a) shows the finite element mesh for the contact analysis method. It has a total of 402 four-node elements in the dam body and 46 four-node elements in the concrete face slab (the concrete face slab is divided into two layers of elements). The mesh for the interface element method is shown in Figure 6(b), where a row of interface elements is placed along the interface between the concrete face slab and the cushion layer. This mesh model has 23 additional interface elements compared to the mesh for the contact analysis model. If the interface elements in Figure 6(b) are assigned a thickness of 0.3 m, the finite element mesh for the thin-layer element method is obtained. Because the length of each element is 12 m, the thin-layer elements have an aspect ratio of 0.025, in the range of 0.01–0.1 [8–11]. The previous mesh models show that the dam body and concrete face slab can be meshed independently for the contact analysis method. This may produce nonmatching nodes on both sides of the interface [15]. However, the thin-layer element and interface element methods usually require matching nodes on both sides of the interface. This model sets zero displacements along the rock base [12].Two-dimensional finite element mesh.
(a)
For the contact analysis method
(b)
For the interface element methodThe computational procedure follows exactly the construction stages shown in Figure5(b). First, blocks ① and ② of the dam body were built up to El.682 m. In each block, layer-by-layer elements were activated to simulate the construction process, and the midpoint stiffness [20] was used for the nonlinear constitutive model. Before placement of the freshly cast stage-I slab, the calculated displacements of the dam body were set to zero, and the calculated stresses were retained. The elements of stage-I slab ③ were then activated, and dam construction continued. The impounding process was simulated by increasing the water level by 10 m in each increment. The same procedure was repeated until completion of the whole dam body.The concrete face slab had an elastic modulus of3×104 MPa and a Poisson’s ratio of 0.2. Table 3 gives the computational parameters of the rockfill materials for the EB model. An elastic modulus of 6 MPa and Poisson’s ratio of 0.2 were used for the materials in the thin-layer elements. The computational parameters for the Goodman interface model are listed in Table 4.Table 3
Computational parameters for the rockfill materials.
Mat. number
Density (kg/m3)
ϕ
′  [°]
Δ
ϕ  [°]
k
n
R
f
k
b
m
IIA
2200
50.6
7.0
1000
0.35
0.71
450
0.24
IIIA
2100
52.5
8.0
900
0.36
0.76
400
0.19
IIIB
2100
51.0
13.0
564
0.35
0.85
204
0.18
IIID
2050
51.0
13.5
432
0.30
0.80
300
−0.18
IIIC
2150
45.0
10.0
250
0.25
0.73
125
0.00Table 4
Parameters of the Goodman interface model.
ϕ  [°]
k
1
n
1
R
f
1
k
n (MPa)
Compression
Tension
30
1000
0.3
1
10000
1


## 5. Comparison of the Three Methods

### 5.1. Deformation of the Dam Body in August 1999

The deformation of the dam body in August 1999 (water level: 768 m) was predicted by the previous three numerical methods. Figure7 compares the contours of the predicted settlement using these numerical methods with in-situ observations. The in-situ observation data used in this study were provided by the HydroChina Kunming Engineering Corporation [21]. Horizontal displacements were measured using indium steel wire alignment horizontal displacement meters, and settlements were measured using water level settlement gauges. As shown in Figure 7, the three numerical methods provided almost identical results and agreed reasonably with the in-situ observation data.Contours of settlement in the dam body in August 1999 (unit: cm).
(a)
Comparison of the interface element, contact analysis, and thin-layer element methods
(b)
Comparison of the contact analysis method and in-situ observationFigure8 shows the locations of the observation points along the interface, where C1-H5, C2-H6, C3-H4, and C4-H2 are the horizontal displacement measurement points and C3-V7 and C4-V4 are the settlement measurement points. The settlement-time curves and horizontal displacement-time curves at typical observation points are displayed in Figures 9 and 10, respectively. The in-situ observations are also plotted for comparison. The three numerical methods predicted almost the same settlements and were in reasonable agreement with the in-situ observations. The horizontal displacements predicted by the three numerical methods were also similar and agreed reasonably with the in-situ observations.Figure 8
In-situ observation points along the interface at0+630 section.Settlement of the dam body along the interface.
(a)
C3-V7 point
(b)
C4-V4 pointHorizontal displacement of the dam body along the interface.
(a)
C1-H5 point
(b)
C2-H6 point
(c)
C3-H4 point
(d)
C4-H2 point


### 5.2. Separation of the Concrete Face Slab from the Cushion Layer

Figures11 and 12 show the separation of the concrete face slab from the cushion layer at different construction stages, predicted by the interface element and contact analysis methods, respectively. Table 5 compares the maximum opening width and depth predicted by the three numerical methods with the in-situ observations. The opening width was measured using a TSJ displacement meter, and the depth was measured manually using a ruler.Table 5
Comparison of the maximum openings.
Stage-I slab
Stage-II slab
Width (m)
Depth (m)
Width (m)
Depth (m)
In-situ observation
0.15
7.2
0.10
5.0
Contact analysis method
0.13
8.0
0.40
14.0
Thin-layer element method*
0
0
0.05
13.2
Interface element method*
0
0
0.08
26.4
*The depth of the tensile stress zone in the interface/thin-layer element is taken as the opening depth, and the relative displacement is taken as the opening width.Separation of the slab from the cushion layer at different stages (interface element method).
(a)
Separation of the stage-I slab
(b)
Separation of the stage-II slabSeparation of the slab from the cushion layer at different stages (contact analysis method).
(a)
Separation of the stage-I slab
(b)
Separation of the stage-II slabThe contact analysis method predicted a maximum opening width of 0.13 m and a depth of 8.0 m for the stage-I slab, which were in good agreement with the in-situ observations. The thin-layer element and interface element methods predicted no opening for the stage-I slab. At the completion of dam body construction, the contact analysis method predicted a maximum opening width of 0.40 m and a depth of 14.0 m for the stage-II slab, while the in-situ observations were much smaller with an opening width of 0.1 m and an opening depth of 5.0 m. The opening widths predicted using the thin-layer element and interface element methods were closer to the in-situ observations. However, the interface element method predicted a much larger opening depth.As shown in Figure11, the opening width and depth were mesh-size dependent for both the thin-layer element and interface element methods because they used element information to determine the separation. The opening depth was the depth of the tensile stress zone, and the opening width was the relative displacement. Therefore, the opening width and depth obtained were used only for reference. Conversely, the contact analysis method regarded the concrete face slab and dam body as independent deformable bodies, and thus the separation could be directly calculated and was independent of mesh size as shown in Figure 12. Therefore, it was concluded that the contact analysis method was reliable and accurate in the prediction of the opening width and depth. In summary, the contact analysis method was a better choice for simulating the separation (opening width and depth) of the concrete face slab from the cushion layer.


### 5.3. Normal Contact Stress along the Interface

The normal contact stress along the interface is compared in Figure13 for the three numerical methods. Figure 13(a) shows the contact stress immediately before casting the stage-II slab and Figure 13(b) at the completion of dam body construction. As shown in Figure 13(a), the maximum normal stress predicted by the thin-layer element method occurs at the middle of the interface between the stage-I slab and the cushion layer, which is not reasonable because the self-weight of the stage-I slab and water pressure should produce a larger normal stress at the bottom as predicted by the contact analysis method. At this stage, the thin-layer element method failed to predict any separation. Furthermore, thin-layer element method predicted a tensile stress zone at the top of the stage-II slab after completion of the dam body construction (Figure 13(b)). Physically, no tensile stress should exist if separation of the two materials occurs. Because the thin-layer element was basically a solid element, it was unsuitable for separation simulation [9]. The interface element method predicted oscillatory normal contact stress at both stages, and the elimination of such oscillation was difficult [3, 22]. In addition, the interface element method could not predict the separation before casting the stage-II slab, and the opening depth was mesh-size dependent. Therefore, both the thin-layer element and interface element methods could not correctly compute the contact stress or the separation.Comparison of normal contact stress along the interface.
(a)
Before casting the stage-II slab
(b)
Completion of the dam body construction


### 5.4. Stresses in the Concrete Face Slab

The stress distribution in the concrete face slab, which was complex because of the deflection of the concrete face slab, was important to the development of cracks. The shear and normal stresses in the concrete face slab at the completion of dam body construction predicted by the three numerical methods, are compared in Figure14. Both normal and shear stresses predicted by the interface element method were oscillatory and nonzero at the top of the slab. The thin-layer element method predicted less oscillatory stresses; however, its normal and shear stresses were also nonzero at the top of the concrete face slab. The magnitude of the stresses predicted by the contact analysis method was much lower than the other two methods, and the normal and shear stresses were zero at the top of the slab. Moreover, the stress distributions for the concrete face slab looked reasonable.Comparison of stresses in the concrete face slab at completion of the dam body construction.
(a)
Shear stress
(b)
Normal stress


## 5.1. Deformation of the Dam Body in August 1999

The deformation of the dam body in August 1999 (water level: 768 m) was predicted by the previous three numerical methods. Figure7 compares the contours of the predicted settlement using these numerical methods with in-situ observations. The in-situ observation data used in this study were provided by the HydroChina Kunming Engineering Corporation [21]. Horizontal displacements were measured using indium steel wire alignment horizontal displacement meters, and settlements were measured using water level settlement gauges. As shown in Figure 7, the three numerical methods provided almost identical results and agreed reasonably with the in-situ observation data.Contours of settlement in the dam body in August 1999 (unit: cm).
(a)
Comparison of the interface element, contact analysis, and thin-layer element methods
(b)
Comparison of the contact analysis method and in-situ observationFigure8 shows the locations of the observation points along the interface, where C1-H5, C2-H6, C3-H4, and C4-H2 are the horizontal displacement measurement points and C3-V7 and C4-V4 are the settlement measurement points. The settlement-time curves and horizontal displacement-time curves at typical observation points are displayed in Figures 9 and 10, respectively. The in-situ observations are also plotted for comparison. The three numerical methods predicted almost the same settlements and were in reasonable agreement with the in-situ observations. The horizontal displacements predicted by the three numerical methods were also similar and agreed reasonably with the in-situ observations.Figure 8
In-situ observation points along the interface at0+630 section.Settlement of the dam body along the interface.
(a)
C3-V7 point
(b)
C4-V4 pointHorizontal displacement of the dam body along the interface.
(a)
C1-H5 point
(b)
C2-H6 point
(c)
C3-H4 point
(d)
C4-H2 point


## 5.2. Separation of the Concrete Face Slab from the Cushion Layer

Figures11 and 12 show the separation of the concrete face slab from the cushion layer at different construction stages, predicted by the interface element and contact analysis methods, respectively. Table 5 compares the maximum opening width and depth predicted by the three numerical methods with the in-situ observations. The opening width was measured using a TSJ displacement meter, and the depth was measured manually using a ruler.Table 5
Comparison of the maximum openings.
Stage-I slab
Stage-II slab
Width (m)
Depth (m)
Width (m)
Depth (m)
In-situ observation
0.15
7.2
0.10
5.0
Contact analysis method
0.13
8.0
0.40
14.0
Thin-layer element method*
0
0
0.05
13.2
Interface element method*
0
0
0.08
26.4
*The depth of the tensile stress zone in the interface/thin-layer element is taken as the opening depth, and the relative displacement is taken as the opening width.Separation of the slab from the cushion layer at different stages (interface element method).
(a)
Separation of the stage-I slab
(b)
Separation of the stage-II slabSeparation of the slab from the cushion layer at different stages (contact analysis method).
(a)
Separation of the stage-I slab
(b)
Separation of the stage-II slabThe contact analysis method predicted a maximum opening width of 0.13 m and a depth of 8.0 m for the stage-I slab, which were in good agreement with the in-situ observations. The thin-layer element and interface element methods predicted no opening for the stage-I slab. At the completion of dam body construction, the contact analysis method predicted a maximum opening width of 0.40 m and a depth of 14.0 m for the stage-II slab, while the in-situ observations were much smaller with an opening width of 0.1 m and an opening depth of 5.0 m. The opening widths predicted using the thin-layer element and interface element methods were closer to the in-situ observations. However, the interface element method predicted a much larger opening depth.As shown in Figure11, the opening width and depth were mesh-size dependent for both the thin-layer element and interface element methods because they used element information to determine the separation. The opening depth was the depth of the tensile stress zone, and the opening width was the relative displacement. Therefore, the opening width and depth obtained were used only for reference. Conversely, the contact analysis method regarded the concrete face slab and dam body as independent deformable bodies, and thus the separation could be directly calculated and was independent of mesh size as shown in Figure 12. Therefore, it was concluded that the contact analysis method was reliable and accurate in the prediction of the opening width and depth. In summary, the contact analysis method was a better choice for simulating the separation (opening width and depth) of the concrete face slab from the cushion layer.


## 5.3. Normal Contact Stress along the Interface

The normal contact stress along the interface is compared in Figure13 for the three numerical methods. Figure 13(a) shows the contact stress immediately before casting the stage-II slab and Figure 13(b) at the completion of dam body construction. As shown in Figure 13(a), the maximum normal stress predicted by the thin-layer element method occurs at the middle of the interface between the stage-I slab and the cushion layer, which is not reasonable because the self-weight of the stage-I slab and water pressure should produce a larger normal stress at the bottom as predicted by the contact analysis method. At this stage, the thin-layer element method failed to predict any separation. Furthermore, thin-layer element method predicted a tensile stress zone at the top of the stage-II slab after completion of the dam body construction (Figure 13(b)). Physically, no tensile stress should exist if separation of the two materials occurs. Because the thin-layer element was basically a solid element, it was unsuitable for separation simulation [9]. The interface element method predicted oscillatory normal contact stress at both stages, and the elimination of such oscillation was difficult [3, 22]. In addition, the interface element method could not predict the separation before casting the stage-II slab, and the opening depth was mesh-size dependent. Therefore, both the thin-layer element and interface element methods could not correctly compute the contact stress or the separation.Comparison of normal contact stress along the interface.
(a)
Before casting the stage-II slab
(b)
Completion of the dam body construction


## 5.4. Stresses in the Concrete Face Slab

The stress distribution in the concrete face slab, which was complex because of the deflection of the concrete face slab, was important to the development of cracks. The shear and normal stresses in the concrete face slab at the completion of dam body construction predicted by the three numerical methods, are compared in Figure14. Both normal and shear stresses predicted by the interface element method were oscillatory and nonzero at the top of the slab. The thin-layer element method predicted less oscillatory stresses; however, its normal and shear stresses were also nonzero at the top of the concrete face slab. The magnitude of the stresses predicted by the contact analysis method was much lower than the other two methods, and the normal and shear stresses were zero at the top of the slab. Moreover, the stress distributions for the concrete face slab looked reasonable.Comparison of stresses in the concrete face slab at completion of the dam body construction.
(a)
Shear stress
(b)
Normal stress


## 6. Conclusions

This study compared the interface treatments in the interface element, thin-layer element, and contact analysis methods, and their numerical performance in predicting deformation, slab separation, contact stress along the interface, and stresses in the concrete face slab in the Tianshengqiao-I concrete-faced rockfill dam, through two-dimensional finite element analysis. Numerical results were also compared with the in-situ observations available. Based on these comparisons, the following conclusions and understanding can be drawn.First, the three numerical methods predicted almost the same settlement and similar horizontal displacement, and the predicted deformation was in good agreement with the in-situ observation data. This indicated that the Duncan EB model used can correctly describe the nonlinearity of this high concrete-faced rockfill dam.Second, interface element method cannot correctly simulate the slab separation. The predicted normal stress along the interface, and stresses in the concrete face slab were oscillatory and not accurate enough for cracking assessment. The thin-layer element method could reasonably predict the normal stress along the interface in some circumstances. However, because solid elements were used, there were intrinsic difficulties in simulating slab separation, and this often led to inaccurate stress distribution in the concrete slab.Third, the contact analysis method could physically and quantitatively simulate the slab separation at different construction stages of the Tianshenqiao-I high CFRD dam. The predicted opening width and depth were in reasonable agreement with the in-situ observations. The normal contact stress along the interface and the stresses in the concrete face slab were reasonable. Furthermore, because no elements were used along the interface, the contact analysis method allowed nonmatching nodes on both sides of the interface and could incorporate complex physical and geometrical properties. The stress distributions obtained could be used for the evaluation of potential cracking risk in CFRDs.The previous discussion indicates that, for contact problems involving large separation or slipping, the contact analysis method (as the most physically realistic approach) is the best numerical method, while the interface element and thin-layer element methods (as simplified contact treatments) are not applicable. Although the performance of these two methods can be largely improved through using more sophisticated constitutive models, applying a tension cut-off criterion, or allowing node-to-node contact, their intrinsic limitations (e.g., contact description based on fixed node pairs) make it difficult for them to obtain satisfactory results for complex contact problems. However, the contact analysis method is a relatively new approach for engineering applications and further studies should be conducted to improve its computational efficiency and stability.



---

*Source: 320890-2013-09-26.xml*