An interfacial layer element for finite element analysis of arch dams

Engineering Structures 128 (2016) 400–414

Contents lists available at ScienceDirect

Engineering Structures journal homepage: www.elsevier.com/locate/engstruct

An interfacial layer element for finite element analysis of arch dams Danni Luo a,b, Yu Hu b, Qingbin Li b,⇑ a b

College of Civil Engineering and Architecture, Guangxi University, Nanning, Guangxi 530004, China State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 15 February 2016 Revised 21 September 2016 Accepted 22 September 2016

An interfacial layer element is developed for finite element analysis of engineering structures with interfaces/joints. In order to simulate complex material property variations and high deformation gradients across thin interface/joint regions using coarse meshes, new shape functions and special formulation of stiffness matrix are developed for interfacial layer elements. Two numerical cases are analyzed using the interfacial layer element method in comparison to analytical solutions and the standard element method using fine meshes. The results show that the use of interfacial layer elements substantially reduces the complexity of mesh refinement around interface/joint regions without compromising on the result accuracy. Furthermore, the structural stability of a real arch dam is analyzed by comparing the numerical results obtained by using interfacial layer elements with field monitoring results, and it is found that the interfacial layer element method has large potential in practical engineering applications. Ó 2016 Elsevier Ltd. All rights reserved.

Keywords: Interface/joint Dam Finite element Shape function

1. Introduction Normally, a dam is constructed in its entirety by combining concrete blocks using different types of joints such as transverse joints, longitudinal joints, and horizontal construction joints. The joints are often the most vulnerable area of the entire structure since they contain many micro-defects in the thin regions, and these dam joints have been an important factor influencing the dam safety [1–3]. Therefore, the research on its local deformation and stress is very essential, and it is important to build a reasonable model for the joints that represents their geometric features, material properties and mechanical behavior. The abovementioned problem is also faced by other similar composite materials or structures [4–6], e.g., masonry structures, repaired concrete structures, rock mass structures or welded components. From this point of view, the research on the modeling for structural joint or bimaterial interface has great scientific and application value in dam engineering or other civil engineering. In the last few decades, some experimental analyses of composite materials have been carried out to explore the mechanical properties of bonded interfaces/joints. Some major characteristics of interfaces/joints are as follows: (1) The interfaces/joints are thin regions with very small thickness. This has been verified by the ⇑ Corresponding author. E-mail addresses: [email protected] (Y. Hu), [email protected] (Q. Li).

(D. Luo),

http://dx.doi.org/10.1016/j.engstruct.2016.09.048 0141-0296/Ó 2016 Elsevier Ltd. All rights reserved.

[email protected]

experimental observations of an alumina/niobium system and an epoxy bonded bimaterial using a high-frequency grating and by high-resolution electron microscopy [7–9]. (2) An interfacial effect region (IER) does exist, and the material properties vary significantly and continuously in this region. Based on the measurement results of the deformation in the local region of bimaterial interfaces, Delale and Erdogan [10] and Kang et al. [9] indicated that complex and continuous material property variations exist in IERs and four different types of distributions are proposed for the bonded bimaterials. (3) High deformation gradients and stress concentrations occur around the narrow interfacial regions. By using some microscopy techniques, the whole-field displacement of bimaterial plates was obtained, then strain and stress were calculated [9,11]. The results show that severe deformation gradients and stress concentrations occur around the interfacial regions under complicated loads. Thus, interfacial damage or debonding is a common failure mode for composite materials or structures. Theoretical analyses have also been carried out to investigate interface/joint problems, and some notable results have been obtained [12,13]. The researches have mainly used two kinds of mechanical models for interfaces/joints: the zero-thickness mechanical model and thin-layer mechanical model [10]. The former assumes that material properties have a sudden jump in value at the interface/joint plane, whereas the latter considers interfaces/ joints as thin layers whose material properties vary significantly and continuously. Although the zero-thickness mechanical model can be used to easily obtain analytical solutions, the effect of the

D. Luo et al. / Engineering Structures 128 (2016) 400–414

microstructure and material properties of interfacial regions on interfacial bond properties cannot be studied by this model, and the theoretical analysis of interfacial cracks carried out using this model shows some unreasonable phenomena (e.g., oscillatory stress singularities and interpenetration of crack faces). In contrast, the thin-layer mechanical model considers a greater number of characteristics of interfaces/joints (e.g., the thickness of interfacial regions and the high material property gradient in IERs), and the aforementioned unreasonable phenomena do not appear when using the thin-layer mechanical model. It is indicated that the geometric features and material properties of interfaces/joints should be represented in modeling interfaces/joints. However, only simple boundary value problems involving interfaces/joints can be solved by theoretical analysis, and theoretical analysis can hardly handle complicated engineering problems. For engineering structures with interfaces/joints, the finite element method (FEM) can be used as an effective means of numerical simulation to perform structural analysis because of its applicability and reliability. In FEM analysis, bimaterial interfaces or structural joints are modeled and simulated by using some special elements. At present, according to Karabatakis and Hatzigogos [14], the simulation of interfaces/joints can be mainly performed by three types of elements: (1) linkage elements [15,16]; (2) interface elements, e.g., zero-thickness interface elements [17,18] and thin-layer interface elements [19,20]; and (3) thin continuum finite elements with either minor or no modification [21,22]. Linkage elements and zero-thickness interface elements can predict the relative movements of interfaces/joints and have the advantage of simple implementation. However, they have some serious deficiencies such as spurious stress oscillations due to inappropriate quadrature schemes, significant errors in the normal stress due to large stiffness in the normal direction, and the arbitrary determination of stiffness parameters without considering the geometric features and material properties of interfaces/joints [23]. Some researchers [19,24] have advocated that thin-layer interface elements are computationally more reliable than zerothickness interface elements because test results indicate that interfaces/joints have some thickness, and thin-layer interface elements can give the distribution of stress and strain and present their concentrations within a finite thin region. In addition, an unreasonable steep stress gradient and inaccurate normal stress

401

do not appear when using the thin-layer elements. However, thin-layer interface elements also present a few challenges, e.g., the numerical ill-conditioning due to low thickness/length ratios and the lack of shear and normal coupling behavior. Thin continuum finite elements can give the complete three-dimensional (3D) state of deformation and stress and can couple shear and normal behavior. However, they considerably raise the complexity of mesh refinement and computational cost [22,25]. Overall, in the FEM analysis of structures with interfaces/joints, the aforementioned elements have exhibited numerical effectiveness and conceptual rationality. However, some major characteristics of bonded interfaces/joints, e.g., thin thickness, complex material property variations of IERs, and severe deformation gradients across interfacial regions, are not fully considered. Moreover, the delicate deformation field and complete stress state cannot be obtained through simple mesh refinement around interfaces/ joints, which affects the evaluation of interface/joint failure or raise the complexity of mesh refinement. To meet the needs of engineering structures, it is important to develop a new special element that is capable of finely simulating complex material property variations and high deformation gradients across thin interfaces/joints regions in structural analysis without mesh refinement around interface/joint regions. Thus the local deformation and stress of interfacial layers and its influence on the behavior of overall structures can be well analyzed. Accordingly, the objective of this study is to develop an interfacial layer element (ILE) for the FEM analysis of the structures with interfaces/joints with low element cost and high calculation accuracy. Furthermore, the new ILE method is used for the performance analysis of real super-high arch dams during their impounding period. 2. Interfacial layer elements 2.1. Modeling method For structures with bonded interfaces/joints, the interfacial region is a very thin layer compared to structure size, and the material properties vary significantly and continuously in the IERs. Thus, high deformation gradients and stress concentrations occur in these regions under complicated loads. The thin-layer mechanical model for bimaterial interfaces/joints is depicted in Fig. 1 (left).

Fig. 1. Thin-layer mechanical model and special numerical model for interfaces/joints.

402

D. Luo et al. / Engineering Structures 128 (2016) 400–414

The accurate simulation of complex material property variations and high deformation gradients across the thin interface/joint regions in FEM analysis is a great challenge. If standard elements are used, a fine mesh resolution is required around interfaces/ joints. In this cases, special enriched elements are an attractive alternative [26,27]. The interface/joint regions and the wide adjacent regions on their two sides are simulated by large-sized twolayer special composite elements, avoiding specially generating fine meshes around interfaces/joints. The special composite elements are placed along the interfaces/joints, i.e., the first layer of elements closest to the interfaces/joints, as illustrated in Fig. 1 (right). Such special composite elements possess continuous material property variations in individual elements and high deformation gradients in the sub-domain of individual elements. Numerous studies have been conducted on composite elements that incorporate nonhomogeneous materials [28–30], e.g., graded finite elements, micro-elements, composite elements, and virtual laminated elements. However, these elements cannot fully meet the specific design requirements for the interfaces/joints described above. For example, the complex piecewise distribution of material parameters cannot be analyzed using graded finite elements, the continuity of material property variations cannot be sufficiently satisfied using micro-elements, high deformation gradients across interfaces/joints cannot be simulated using virtual laminated elements, and the degrees of freedom increase excessively because of complex material property variations when using composite elements. Therefore, in this study, new special composite elements called ILEs are developed for the modeling of interfaces/joints using coarse meshes, which well considers the physico-mechanical characteristics of interfacial layers, especially the local complex inhomogeneity of material property across thin interface/joint regions. The ILEs possess two main characteristics: (1) New shape functions: The shape functions in ILEs account for the rapid variation in the displacement field inside IERs and the smooth variation outside IERs by using a piecewise interpolating field that incorporates special logarithmic and standard polynomial behavior, respectively. (2) Special formulation of the stiffness matrix and material property matrix: To account for complex material property variations at the element level, dense sub-elements are set up in an element according to the micro-element method.

However, to show the continuity of actual material property variations, the determination of the material properties is modified to direct sample at the Gauss points of sub-elements. The element stiffness matrix is obtained by adding the transformed stiffness matrices of all sub-elements. 2.2. Shape functions of interfacial layer elements At the element level, the accuracy of numerical approximation is dependent on the compatibility of the assumed shape functions with the exact displacement field. In this section, to obtain accurate numerical results for the structures with nonhomogeneous interfaces/joints using large-sized two-layer ILEs, the shape functions in the ILEs are derived from the exact displacement solution of plates with spatially varying material property fields. 2.2.1. Some analytical solutions for nonhomogeneous elasticity Boundary value problems involving nonhomogeneous materials are difficult to solve, so only the analytical solutions for plates with infinite length and finite width subjected to certain loads are available in the literature [28]. The loading conditions include remote tension, bending, and a fixed grip in the direction perpendicular to the material property gradient, as shown in Fig. 2a–c. Both exponential and linear material variations are considered. Remarkably, a new analytical solution is derived for plates with linear material variations subjected to tension in the direction parallel to the material property gradient (Fig. 2d) in this study. In this study, the material property variations across interfaces/ joints are assumed to be piecewise linear (Fig. 1). Therefore, the displacement field for plates with linear material variations is mainly analyzed here. The analytical solutions for the four analysis cases (Fig. 2a–d) with linear material variations are expressed as follows: Case an and case b

ux ðx; yÞ ¼ m


 A 2 A x þ Bx  y2 2 2

ð1Þ

uy ðx; yÞ ¼ ðAx þ BÞy

Fig. 2. Analysis cases of plates with material property variations: (a) tension loading in y-axis; (b) bending loading in y-axis; (c) fixed grip in y-axis; and (d) tension loading in x-axis.

403

D. Luo et al. / Engineering Structures 128 (2016) 400–414

Case c

ux ðx; yÞ ¼ me0 x

ð2Þ

uy ðx; yÞ ¼ e0 y

2.2.2. Special shape functions of interfacial layer elements The shape functions in ILEs are constructed based on the analytical solutions. A 3-D ILE has twelve nodes and three degrees of freedom per node, as shown in Fig. 3a. The nodal displacement vector is given by

fdg ¼ ½ u1

Case d ux ðx;yÞ ¼ 8 rt ðx  x1 Þ > E1 > > h  i > >  E0 E1 ðx3 x2 Þ ðE1 x3 E0 x2 Þ rt > > < E1 ðx2  x1 Þ þ ðE0 E1 Þ rt ln  E1 ðx3 x2 Þ x þ E1 ðx3 x2 Þ      E0  ðx3 x4 Þ  ðE0 E2 Þ ðx3 x4 Þ E2 x3 E0 x4  rt > þ ðx  x Þ þ r ln r ln x þ     > 2 1 t t E E Þ E E Þ ðx x Þ E ðx x Þ ðE ðE E > 0 1 1 0 2 0 3 4 0 3 4 > > 1     > > : rt ðx2  x1 Þ þ ðx3 x2 Þ rt ln E0  þ ðx3 x4 Þ rt ln E2  þ rt ðx  x4 Þ E1 E1 E0 E2 ðE0 E1 Þ ðE0 E2 Þ

x1 6 x 6 x2 x2 6 x 6 x3 x3 6 x 6 x4 x4 6 x 6 x5

uy ðx;yÞ ¼ 0 ð3Þ

where ux and uy are the displacements in the x- and y-directions, respectively; E and m are the elastic modulus and Poisson’s ratio, 1

respectively; A is a constant with the unit ½length ; and B is a dimensionless constant. The parameters A and B are determined by the equations expressed in terms of the material properties, plate dimensions, and loads [28]. Based on the above analytical solutions, two characteristics of the displacement field can be described as follows: (1) Based on Eq. (3) for case d, different behaviors of the displacement field are represented inside ½x2 ; x4  (i.e., IER) and outside the regions. Apparently, the variation in the displacement field is very rapid in IERs because of a high material property gradient, whereas the variation in the displacement field is smooth outside the regions. Thus, a piecewise interpolating field can be adopted in an individual element to numerically approximate the distribution of the displacement field. (2) For the nonhomogeneous regions with linear material variations (i.e., IER), a polynomial behavior of the displacement field is represented under the loading in the direction perpendicular to the material property gradient (Eqs. (1) and (2) for cases a–c), whereas a logarithmic behavior of the displacement field is represented under the loading in the direction parallel to the material property gradient (Eq. (3) for case d). It is suggested that the special log-linear mixed behavior should be incorporated into shape functions to modify the numerical approximation. In the same way, the characteristics of the displacement field for exponential material variations in IERs can also be obtained, and an exponential-linear mixed behavior of the displacement field is represented. However, this study focuses on the interfaces/joints with linear material variations in IERs.

v1

w1

   u12

v 12

w12 T

ð4Þ

Four nodes are added to the four hexahedral edges parallel to the material property gradient to approximate the complex displacement field in this direction. The position of these added nodes is determined by the boundary of the IERs, indicating that these nodes are adaptive to different mesh sizes and thicknesses of IERs. To clarify this, a two-dimensional (2-D) ILE for the left zones of interfaces/ joints is illustrated in Fig. 3b. The 2-D ILE is divided into parts inside and outside the IERs by the line from node 5 to node 6. According to the analysis of the displacement field for nonhomogeneous materials in Section 2.2.1, the shape functions in ILEs can be derived from the interpolating field as follows: (1) A piecewise interpolating field is adopted in the direction parallel to the material property gradient (n-axis). The linear representation and log-linear mixed representation of the interpolating field are respectively adopted inside and outside the IERs, as depicted in Fig. 4. The modified interpolating field takes the form of Eq. (5). It is shown that high deformation and strain gradients can be presented in the sub-domain of individual elements. (2) The standard linear representation of the interpolating field is adopted in the direction perpendicular to the material property gradient (g-axis and f-axis). In 2-D or 3-D analysis, hexahedral or quadrilateral elements can be obtained from the simple product of an onedimensional (1-D) special element with log-linear mixed shape functions in the direction parallel to the material property gradient and one or two 1-D standard elements in the direction perpendicular to the material property gradient. ILEs for the left zones of the interfaces/joints (left ILEs) are



uðnÞ ¼

c1 þ c2 n

1 6 n 6 nm

c3 þ c4 n þ c5 ln jd1 n þ d2 j nm 6 n 6 1

0 E1 Þ ; d1 ¼ EðE1 ð1n mÞ

ð5-aÞ

E0 nm d2 ¼ EE11ð1n mÞ

ILEs for the right zones of the interfaces/joints (right ILEs) are

 uðnÞ ¼

c3 þ c4 n þ c5 ln jd1 n þ d2 j 1 6 n 6 nm c1 þ c2 n

2 E0 ; d1 ¼ E0Eð1þn mÞ

nm 6 n 6 1

ð5-bÞ

þE0 nm d2 ¼ EE02ð1þn mÞ

where uðnÞ is the displacement in the n-direction; d1 and d2 are constants determined by material property variations; and c1 ; . . . ; c5 are undetermined constants.

Fig. 3. Schematic diagrams of ILEs: (a) 3-D interfacial layer elements and (b) 2-D interfacial layer elements for left zones of interfaces/joints.

404

D. Luo et al. / Engineering Structures 128 (2016) 400–414

Fig. 4. The material property field, interpolating field and strain field in n-axis: (a) interfacial layer elements for left zones of interfaces/joints and (b) interfacial layer elements for right zones of interfaces/joints.

The interpolating field in Eq. (5) yields the following shape functions for the 1-D 3-node element in the direction parallel to the material property gradient (Eq. (6)). Left ILEs:

8 n n m 1 6 n 6 nm > 1þn > >   < m   E0  1  E0  nm N1 ðnÞ ¼ ð1nm Þ ln E1 ð1nm Þ ln E1 nþln jd1 nþd2 j >  i h nm 6 n 6 1 > > E  : 1 ln E0  ð1þnm Þ d1 þð1n mÞ 1 8 1 6 n 6 nm <0 d1 nmhd1 nþln jd1 nþd2 j i N2 ðnÞ ¼ nm 6 n 6 1 : ð1n Þ d þ 1 ln E0  m 1 ð1nm Þ E1 8 1þn 1 6 n 6 nm > > > 1þnm   h  i <  E0   E0  N3 ðnÞ ¼ d1 ð1þnm Þþln E1 þ d1 ð1þnm Þþln E1  n2 ln jd1 nþd2 j >  i h nm 6 n 6 1 > > E  : 1 ð1þnm Þð1nm Þ d1 þð1n ln E0  mÞ 1 ð6-aÞ

functions of the material parameters (E1 ; E2 ; E0 ) and the characteristic length scale of IERs (nm ). The values of the n coordinates of node 3 (nm ) can be set to any value from 1 to 1. Thus, the shape functions are adaptive to different thicknesses and different material property variations of IERs. The 2-D or 3-D shape functions of each element are derived by multiplying the 1-D shape functions belonging to the adjacent element edges of the respective node, as given by Eq. (7).

ui ¼

n X Nj dij ; Nðn; g; fÞ ¼ NðnÞ  NðgÞ  NðfÞ

ð7Þ

j¼1

where ui is the nodal displacement corresponding to node i and N represents the shape functions. Each shape function depends on the coefficients nm and E1 ; E2 ; E0 of the adjacent element edges in P the n-axis. Partition of unity (e.g., ni¼1 N i ¼ 1) and Kronecker delta property (e.g., N ij ¼ dij ) hold for ILEs. In addition, it is noted that the geometry of ILEs is approximated by the linear shape functions. The piecewise linear and standard linear shape functions are adopted for the n-axis and gðfÞaxis, respectively.

Right ILEs:

8   E E E > ln 2 n d 0 þd 0 nln jd1 nþd2 j > < E0  mh 1 E2 1 E2  i 1 6 n 6 nm E  E N1 ðnÞ ¼ 1 ð1þnm Þ d1 E0 þð1þn ln E2  Þ m > 2 0 > : 0 nm 6 n 6 1     8 E2  1 E2  1 > ln þ ln nln jd nþd > ð1þnm Þ E0  ð1þnm Þ E0  1 2j > <  i h 1 6 n 6 nm   E E 1 N2 ðnÞ ¼ ln 2 ð1nm Þ E0 d1 þð1þn Þ  E0  m 2 > > > : nnm nm 6 n 6 1 1n  i h  i 8 hm E2  E0 E  E0 >  E d1 ð1nm Þþln E   E d1 ð1nm Þþln E2  nþ2 ln jd1 nþd2 j > > 0 0 < 2  i h2 1 6 n 6 nm E  E 1 N3 ðnÞ ¼ ln E2  ð1nm Þð1þnm Þ E0 d1 þð1þn Þ m 2 0 > > > : 1n nm 6 n 6 1 1nm ð6-bÞ Unlike the standard shape functions expressed in terms of the local nondimensional variable n only, the special shape functions are also

2.3. Stiffness matrix of interfacial layer elements To account for complex and continuous material property variations at the element level, the core concept of the micro-element method is adopted, and an improvement is made in the formulation of the material property matrix. An ILE can be divided into a number of nonuniform sub-elements, as shown in Fig. 3b. The sub-elements have their own three degrees of freedom and can be defined as different types of elements. Eight-node hexahedron elements with standard linear shape functions and logarithmic shape functions are selected for the sub-elements inside and outside the IERs, respectively. To avoid increasing the computational elements and degrees of freedom, the degrees of freedom of all sub-elements’ nodes are translated into the degrees of freedom of ILEs by the compatibility conditions. The stiffness matrix of ILEs can be obtained by adding the transformed stiffness matrices of all sub-elements.

405

D. Luo et al. / Engineering Structures 128 (2016) 400–414

2.3.1. Stiffness matrix of sub-elements The nodal displacement vector of 3-D sub-elements is given by

v1

fde g ¼ ½ u1

   u8

w1

v8

T

ð8Þ

w8 

Then, the displacements of sub-elements can be written as

fue g ¼ ½Ne fde g

ð9Þ

For the sub-elements outside the IERs, the standard linear shape functions are

N0i ¼

1 ð1 þ ni nÞð1 þ gi gÞð1 þ fi fÞ ði ¼ 1; 2; . . . ; 8Þ 8

ð10-aÞ

For the sub-elements inside the IERs, the logarithmic shape functions are

0

 1n 0 1þn 0 1 di1 ð 2 xis þ 2 xie Þþdi2  C 0 B di1 x0is þdi2  N0i ¼ 14 @1  Að1 þ gi gÞð1 þ fi fÞ ði ¼ 1; 4; 5; 8Þ di1 xie þdi2  ln   di1 x0 þdi2 is 0  1n 0 1þn 0 1 di1 ð 2 xis þ 2 xie Þþdi2  ln  C 0 þd d x B  i1 0is i2  N0i ¼ 14 @ Að1 þ gi gÞð1 þ fi fÞ ði ¼ 2; 3; 6; 7Þ di1 xie þdi2  ln   d x0 þd ln 

i1 is

i2

ð10-bÞ where N0i represents the shape functions; ðni ; gi ; fi Þ denote the local coordinates of node i; x0is and x0ie are the x-coordinate values of the start and end nodes belonging to the adjacent element edges of node i, respectively; and di1 and di2 are the constants related to the adjacent element edges of node i and are determined according to Eq. (5). Strains are obtained from the displacements by differentiation as follows:

feg ¼ ½L½Ne fde g ¼ ½Be fde g

ð11Þ

where ½Be  is the strain-displacement matrix. Because the geometry of ILEs is approximated by linear shape functions, the Jacobian matrix ½J remains the same as the standard eight-node hexahedron elements. Then, the strain-stress relationships are given by

frg ¼ ½De feg ¼ ½De ½Be fde g

where ½D  is the constitutive matrix, which is a function of the spatial coordinates for the nonhomogeneous material. To show the continuity of actual material property variations at the element level, the determination of the material properties is improved to direct sample at the Gauss points of sub-elements.

De ¼ De ðx; y; zÞ

ð13Þ

The principle of virtual work yields the following element stiffness equations:

½K e fde g ¼ fF e g

½K  ¼

e T

Xe

Z

1

Z

Z

1

¼ 1

e

e

½B  ½D ½B dV 1

1

T

½Be  ½De ½Be jJjdndgdf

X

ð15Þ

1

Z fF e g ¼

fd g ¼

½   ueit

#T n X    ¼    Nj ðnt ; gt ; ft Þdij    ¼ ½Afdg T

j¼1

ði ¼ 1; 2; 3; j ¼ 1; 2; . . . ; 8; n ¼ 1; 2; . . . ; 12Þ

ð17Þ

where ½A is the transformation matrix from the nodal displacement vector of the sub-elements to that of the ILEs. Then, the strain-displacement matrix, stiffness matrix, and load vector of sub-elements with respect to the nodal displacement vector of the sub-elements can be transformed into the corresponding ones with respect to the nodal displacement vector of the ILEs.

e fue g ¼ ½Ne fde g ¼ ½Ne ½Afdg ¼ ½ Nfdg

ð18Þ

e fee g ¼ ½Be fde g ¼ ½Be ½Afdg ¼ ½ Bfdg

ð19Þ

½K   ¼

ZZZ X

e T ½De ½ Bd e X¼ ½ B

ZZZ

T

X

½AT ½Be  ½De ½Be ½AdX

¼ ½AT ½K e ½A fF  g ¼ ½AT ½F e 

ð20Þ ð21Þ

2.3.3. Stiffness matrix of interfacial layer elements After the transition analysis is completed for all sub-elements, the stiffness matrix and load vector of the ILEs can be obtained by adding the transformed stiffness matrices and load vectors of all sub-elements:

½K ¼

X  ½K 

ð22Þ

e

fFg ¼

X  fF g

ð23Þ

e

3. Numerical verification The 2-D and 3-D FEM codes for ILEs developed in this study are implemented in ABAQUS software as user subroutines. In this section, two numerical examples of bimaterial plates are presented to evaluate the validity and accuracy of the ILEs. The comparative analysis of the ILE method using a simpler mesh with analytical solutions or the thin-layer standard element method using a finer mesh is carried out. In addition, the merits of ILEs in mesh simplification is discussed.

ð14Þ

where ½K e  and fF e g are the element stiffness matrix and the load vector of the sub-elements, respectively.

ZZZ

"

e

ð12Þ

e

e

of freedom of ILEs by the compatibility conditions. The transition method is as follows: the local coordinates of sub-elements’ nodes (nt ; gt ; ft ) are substituted into the displacement mode of ILEs (Eq. (7)); then, the nodal displacement vector of the sub-elements fde g (Eq. (8)) can be given by that of the ILEs fdg (Eq. (4)).

Z NeT ðn; g; fÞpðn; g; fÞdX þ

NeT ðn; g; fÞtðn; g; fÞdS

ð16Þ

S

2.3.2. Transition of sub-element degrees of freedom A number of sub-elements are set up in ILEs, and the degrees of freedom of all sub-elements’ nodes are translated into the degrees

3.1. Numerical case 1: plate subjected to tension loads Fig. 5a illustrates a bimaterial plate with material variation in the x-direction subjected to uniform tension in the direction parallel to the material property gradient. The calculation parameters for case 1 are as follows: (1) The plate size is 500 mm  200 mm (L  W), and the thickness of the IERs is 48:8 mm. (2) The load rt is 1  106 Pa. (3) The elastic modulus varies along the x-direction as

8 EðxÞ ¼ E1 x1 6 x 6 x2 > > > > < EðxÞ ¼ E1 þ ðxx2 Þ ðE0  E1 Þ x2 6 x 6 x3 ðx3 x2 Þ 4Þ > EðxÞ ¼ E2 þ ðxðxx ðE0  E2 Þ x3 6 x 6 x4 > > 3 x4 Þ > : EðxÞ ¼ E2 x4 6 x 6 x5

ð24Þ

406

D. Luo et al. / Engineering Structures 128 (2016) 400–414

Interfacial layer element 1

(a)

Joint

Interfacial layer element 2

(b)

Fig. 5. The calculation condition and element mesh for case 1: (a) the calculation condition and (b) the element mesh using interfacial layer elements.

Fig. 6. Comparison of the displacement, strain and stress results obtained from numerical simulation and theoretical analysis: (a) ux at y ¼ 0:02 m; (b) rxx at y ¼ 0:011 m; and (c) exx at y ¼ 0:011 m.

407

D. Luo et al. / Engineering Structures 128 (2016) 400–414

where 10

10

9

E1 ¼ 2:5  10 Pa; E2 ¼ 3:0  10 Pa; E0 ¼ 3:0  10 Pa x1 ¼ 0:25 m; x2 ¼ 0:0244 m; x3 ¼ 0; x4 ¼ 0:0244 m; x5 ¼ 0:25 m In addition, Poisson’s ratio m is zero. In the numerical simulation using ILEs, the element mesh for case 1 consists of two ILEs placed along the interfaces/joints (Fig. 5b). An ILE is divided into 5  2 and 5  14 (row  column) sub-elements inside and outside the IERs, respectively. A 2  2 Gauss quadrature is employed for each sub-element. Fig. 6 shows the comparison of the displacement, strain, and stress results obtained from numerical simulation and theoretical analysis (Eq. (3)). It can be seen from Fig. 6 that the ILE numerical results exactly match the analytical solutions in terms of the distribution all over the structures as well as the detailed characteristics of the IERs. The relative error of ux at characteristic points A (0.0244, 0) and B (0.0244, 0) are 0.004% and 0.001%, respectively. It is shown that a high accuracy is achieved by using ILEs with a very coarse mesh. According to Kim et al. [28] and Santare et al. [29], for this case, the graded four-node isoparametric elements with linear shape functions provide fluctuating stress distribution in the direction parallel to the material property gradient and predict sharp stress jumps at the element borders; these results are not exactly consistent with the analytical solution. However, it can be seen from Fig. 6 that the stress results for the ILEs exactly match the analytical solution. This is because the exact displacement field and actual material property field are perfectly modeled by the ILEs. However, the linear interpolating field of graded finite elements cannot approximate the logarithmic behavior of the displacement solution. This indicates that the improvement in the shape functions in ILEs results in superior performance at the element level for some composite elements with standard shape functions (e.g., graded finite elements, micro-elements, and virtual laminated elements).

plane stress elements, white elements in Fig. 7b) are generated for the entire structure. A ILE is divided into 5  2 and 5  8 (row  column) sub-elements inside and outside the IERs, respectively. A 2  2 Gauss quadrature is employed for each sub-element. In the standard element model (Fig. 7c), a 16-layer mesh is generated in the thickness direction of the interfaces/joints, and 6200 standard elements (four-node plane stress elements, white elements in Fig. 7c) are generated for the entire structure. A 2  2 Gauss quadrature is employed for each standard element. Fig. 8 shows the comparison of the displacement results obtained by using different interface/joint modeling. The two sets of results obtained by using ILEs with a simpler mesh and using standard elements with a finer mesh are found to be in good agreement. According to the relative deviation defined by Eq. (25), the relative deviation of ux at characteristic points A (0.0075, 0) and B (0.0075, 0) (black points in Fig. 7a) are 0.08% and 0.09%, respectively. The relative deviation of uy at points A and B are 2.25% and 2.18%, respectively. The relative deviation of uy not only results from the different elements used for the interface/joint modeling but also results from the great difference in the standard element mesh outside the IERs (white zone in Fig. 7b). When an identical mesh is generated outside the IERs for the two numerical models, the relative deviation of uy at points A and B are reduced to 0.04% and 0.02%, respectively. This indicates that the ILE method is efficient and shows little relative deviation under complicated calculation conditions.

3.2. Numerical case 2: plate subjected to complicated loads Because the analytical solutions of nonhomogeneous plates subjected to complicated loads are not available, the thin-layer standard element method with a fine mesh, in which step-wise constant approximation is used to simulate continuous material property variations, is employed to provide a reference case. Then, the numerical results obtained using the ILEs with a simpler mesh are compared with those obtained using the standard elements with a finer mesh. Fig. 7a illustrates a plate with material variation in the x-direction subjected to complicated loads. The calculation parameters for case 2 are as follows: (1) The plate size is 500 mm  200 mm (L  W), and the thickness of the IERs is

(a) Interfacial layer element

Joint

Standard element

(b)

15 mm. (2) The loads are q1 ¼ 1  106 Pa, q2 ¼ 0:5  106 Pa,

Joint

F 1 ¼ 0:3  106 N, F 2 ¼ 0:3  106 N, and F 3 ¼ 0:1  106 N. (3) The variation in the elastic modulus along the x-direction can be expressed by Eq. (24), where

Standard element

E1 ¼ 2:5  1010 Pa; E2 ¼ 3:0  1010 Pa; E0 ¼ 3:0  109 Pa x1 ¼ 0:25 m; x2 ¼ 0:0075 m; x3 ¼ 0;

x4 ¼ 0:0075 m;

x5 ¼ 0:25 m In addition, Poisson’s ratio is constant (m ¼ 0:17) for the zone ½x1 ; x5 . The element mesh for case 2 using the ILEs and standard elements is shown in Fig. 7. In the ILE model (Fig. 7b), large-sized two-layer ILEs (gray elements in Fig. 7b) are placed along the interfaces/joints, and 20 ILEs and 100 standard elements (four-node

(c) Fig. 7. The calculation condition and element mesh for case 2: (a) the calculation condition; (b) the element mesh using interfacial layer elements with a simple mesh; and (c) the element mesh using standard elements with a fine mesh.

408

D. Luo et al. / Engineering Structures 128 (2016) 400–414

Fig. 8. Comparison of the displacement results obtained by using different interface/joint modeling: (a) ux at y ¼ 0 m and (b) uy at y ¼ 0 m.

Fig. 9. Comparison of the strain and stress results obtained by using different interface/joint modeling: (a) rxx at y ¼ 0:011 m; (b) ryy at y ¼ 0:011 m; (c) exx at y ¼ 0:011 m; and (d) eyy at y ¼ 0:011 m.

409

D. Luo et al. / Engineering Structures 128 (2016) 400–414

  uILE  uSE    100 e ¼  uSE 

ð25Þ

where uILE and uSE are the displacement results obtained by using the ILEs and standard elements, respectively. Fig. 9 shows the comparison of the strain and stress results obtained by using different interface/joint modeling. The strain and stress distributions obtained by using ILEs with a simpler mesh and using standard elements with a finer mesh are in good agreement. The high deformation gradients and stress concentrations in the IERs can be depicted by using the ILEs with a very coarse mesh,

and the strain and stress values obtained by using the ILEs at the middle plane of the interfaces/joints are close to those obtained by using the standard elements. 3.3. Discussion on merits of interface layer element in mesh simplification To quantitatively assess the effectiveness of ILEs in mesh simplification, several numerical cases with different material variations are analyzed using the ILEs and standard elements, and the element costs of the numerical models are compared.

Table 1 Comparison of element cost for the modeling of interfaces/joints by using interfacial layer elements and standard elements. Case no.

Calculation parameters

Layers of standard element in IERs

Number of standard elements in IERs

Number of ILE elements in IERs

Relative deviation of ux (%)

Relative deviation of uy (%)

Case 1

E0 = 3 GPa E1 = 25 GPa E2 = 30 GPa t = 0.0075 m

1 4 8 16

30 60 100 180

20

3.29 1.58 0.62 0.19

0.58 0.23 0.05 0.04

Case 2

E0 = 0.6 GPa E1 = 25 GPa E2 = 30 GPa t = 0.0075 m

1 4 8 16

30 60 100 180

20

9.17 6.51 5.70 3.17

1.79 1.23 1.05 0.53

Fig. 10. Numerical models of the Xiluodu arch dam: (a) dam-foundation system model; (b) the modeling of joints using interfacial layer elements; (c) the modeling of joints using thin-layer standard elements; and (d) monolithic structure without the joints.

410

D. Luo et al. / Engineering Structures 128 (2016) 400–414

A plate with the same basic geometry, boundary conditions, and material variation like numerical case 2 (Section 3.2) is simulated. The

plate

is

subjected

to

pressures

q1 ¼ 1  106 Pa

and

6

q2 ¼ 0:2  10 Pa on faces x ¼ 0:25 m and y ¼ 0:1 m, respectively. When the numerical cases are calculated using the standard elements (four-node plane stress elements, CPS4), the element layers in the IERs are gradually increased to obtain increasingly accurate results. Then, the results using the standard elements are compared with those using the large-sized two-layer ILEs, as listed in Table 1. Table 1 indicates that the displacements at the characteristic point Pð0:0075; 0:0Þ obtained using the standard elements become increasingly close to those obtained using the ILEs as the number of standard elements in the IERs increases. The low relative deviation of the two results can be obtained until the 16-layer standard elements are generated in the IERs. In other words, the results

Table 2 Physical-mechanical parameters of the concrete dam, rock masses and inter/intrastrata weak zones. Material

Density (t/m3)

Elastic modulus (GPa)

Poisson’s ratio

Dam concrete Outlets (reinforced concrete) Joints (middle plane) Left bank-1 (reinforced foundation) Left bank-2 (reinforced foundation) Right bank-1 (reinforced foundation) Right bank-2 (reinforced foundation) Grouted riverbed-1 Grouted riverbed -2 Grouted riverbed -3 Replaced foundation Dam toe concrete blocks Type II rock mass Type III1 rock mass Type III2 rock mass C2 weak zone LC3 weak zone RC3 weak zone C5 weak zone LC6 weak zone RC6 weak zone LC7 weak zone RC7 weak zone LC8 weak zone RC8 weak zone C9 weak zone P2BN weak zone

2.663 2.663 2.663 2.75 2.75 2.75 2.75 2.75 2.75 2.75 2.75 2.663 2.85 2.85 2.75 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4

46 46 15 25 25 25 25 17 17 17 38 46 21.45 15.275 7.15 0.5 0.8 0.8 0.8 0.8 0.9 1.7 1.3 0.8 0.85 0.5 1.5

0.167 0.167 0.167 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.25 0.167 0.2 0.25 0.28 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3

obtained by the abundant layers of the standard elements in the IER can be obtained by just two-layer ILEs along the interfaces/ joints. Therefore, it is concluded that ILEs show good performance in terms of decreasing the complexity of mesh refinement, and they can simplify the meshing procedure and reduce the computational elements. Furthermore, the use of ILEs would achieve its greatest potential for mesh simplification in 3-D cases. Based on the two numerical cases described above, it is inferred that the ILEs proposed in this paper are very effective for the deformation and stress analysis of structures with interfaces/joints. The use of ILEs substantially reduces the complexity of mesh refinement around interface/joint zones without compromising on the qualitative behavior of the deformation and stress field or on the result accuracy. 4. Engineering application 4.1. Presentation of the Xiluodu hydropower station The Xiluodu hydropower station is located on the Jinshajiang River, Leibo County, Sichuan Province in southwestern China. It is the third largest power station in the world and is capable of generating 13.86 million kW of power. The overall volume of the reservoir is as large as 12.67  109 m3. The principal structure consists of a double-curvature arch dam with a height of 285.5 m and a crest length of 603 m with spillways and underground powerhouses. 4.2. Numerical models and material parameters The Xiluodu dam is a double-curvature arch dam with 31 monoliths and 30 transverse joints. The thickness of the dam is 14 m at the crest and 60 m at the lowest foundation. 3-D numerical models have been established for the dam-foundation system, as shown in Fig. 10a. The transverse joints are modeled by largesized ILEs containing the joints (Fig. 10b) or by thin-layer standard elements (Fig. 10c), or the dam structure is modeled without the joints (Fig. 10d). In these models, the complex geological and fracture features in the bedrock are explicitly simulated, which include the different types of rock masses, poured concrete that replaced the excavated rock masses, and major faults. The mechanical parameters of the concrete dam, rock masses, and inter/intrastrata weak zones are listed in Table 2. Some major material divisions are illustrated in Fig. 11. According to the geometric features and material properties of the transverse joints, the thickness of

Dam joints

Dam concrete

Left bank 1

Right bank 1 Right bank 2

Left bank 2

Outlets Grouted riverbed-1

Dam toe concrete blocks Grouted riverbed -2 Grouted riverbed -3

Downstream face Fig. 11. Material divisions of the Xiluodu arch dam.

411

D. Luo et al. / Engineering Structures 128 (2016) 400–414

the transverse joints on each side is approximately once the maximum aggregate size of dam concrete and the elastic modulus of the middle plane is approximately 1/3th that of dam concrete in these calculations. In our studies, the effect of impoundment on the behavior of the dam-foundation system is mainly analyzed, so the main calculation conditions of the Xiluodu arch dam are as follows: (1) The pouring of dam concrete has been completed, and arch-closure grouting also has been finished. So the grouted joints are modeled in the dam numerical analysis. The surface roughness and shear keys of the construction joints are not considered in this studies. (2) During the impounding period of the Xiluodu project, the dam has gone through twice rapid impounding process from dead water level to normal water level. So the behavior of the damfoundation system is calculated at dead water level and normal water level. Moreover, the increments of numerical results at different water levels are compared with field monitoring results to validate the results. Important loads in the impounding period are considered in the calculation, e.g., the essential dead weight, water pressure, and the effect of the shrinkage of river valley and foundation settlement. These loads are specifically considered as follows: (1) The dead weight of whole dam is applied at once; (2) Water pressure at dead water level or normal water level is applied on the upstream surface of the dam; (3) The effect of the shrinkage of river valley is considered by applying the actual monitoring deformation on the two sides of dam foundation. (4) The effect of foundation settlement is considered by applying the pressure of reservoir water on the surface of dam base.

4.3. Results and analyses 4.3.1. Impounding effects on dam deformation Fig. 12 shows the dam deformation at the dead water level (540 m) and normal water level (600 m) obtained by the ILE models. The results indicate the following: (1) The dam incline upstream due to the dam self-weight at the dead water level. As the water level increases the dam continually deforms downstream. (2) When the water level is 600 m, the deformation of

Fig. 13. Comparison of the displacement increments at monolith #15 obtained from numerical simulation and field monitoring.

Y Displacement (unit: m)

Positive value: downstream

Minimum principal stress (unit: Pa)

Positive value: tension

Upstream face

(a) Y Displacement (unit: m)

Downstream face Positive value: downstream

(a) Minimum principal stress (unit: Pa)

Upstream face

(b) Fig. 12. The displacement along the river direction of the dam obtained by the interfacial layer element models: (a) at dead water level and (b) at normal water level.

Positive value: tension

Downstream face

(b) Fig. 14. The minimum principal stress of the dam obtained by the interfacial layer element models: (a) at dead water level and (b) at normal water level.

412

D. Luo et al. / Engineering Structures 128 (2016) 400–414

the crown cantilever is the most notable and the maximum downstream displacement appears at about 1/2 of arch height, and the displacements of the dam monoliths located at the two banks are basically similar. The comparison between the displacement increments at the crown cantilever (monolith #15) obtained from numerical simulation and field monitoring is illustrated in Fig. 13, when the water level rises from EL 540 to EL 600 m. The displacement increments obtained by the three numerical models have similar distribution, and the values are considerably close to each other. The displacement distributions of numerical results and monitoring results are similar, and the maximum displacement appears at the same height of the monolith. 4.3.2. Impounding effects on dam stress Fig. 14 shows the minimum principal stress of the dam at the dead water level and the normal water level obtained by the ILE models. During the impounding process, the dam is mainly compressed (negative value in Fig. 14), which indicates an enhancement of the dam thrust into the banks. Compressive stress in the arch direction increases while compressive stress in the cantilever direction decreases. The highest compressive stress appears on the downstream surface around the bottom of the dam monoliths located on the steep slope of the two banks at the normal water level. Based on the data of the strainmeters embedded at the dam site, the dam stress is calculated using a deformation method. Fig. 15 compares the stress of the characteristic points at EL 442.2 m inside the monoliths obtained from numerical simulation and field monitoring. It is found that the dam stress distributions obtained by the three numerical models are similar, and the normal stress obtained from the numerical simulation and field monitoring are in agreement. During the impounding process, the normal stress in the tangential directions of arch (rtt ) increases significantly and the normal stress in the vertical direction (rvv ) decreases, indicating an enhancement of the effect in the arch direction and a reduction of the effect in the cantilever direction. 4.3.3. Impounding effects on transverse joints The stress of the transverse joints between dam monoliths at the dead water level and at the normal water level obtained by the ILE models is shown in Fig. 16. From the figure, the following

Maximum principal stress (unit: Pa)

Positive value: tension

Upstream face

(a) Maximum principal stress (unit: Pa)

Positive value: tension

Upstream face

(b) Fig. 16. The maximum principal stress of transverse joints obtained by the interfacial layer element models: (a) at the dead water level and (b) at normal water level.

conclusions are drawn: (1) The transverse joints are mainly under compression (negative value in Fig. 16) over the entire dam at the dead water level, but the transverse joints located on the steep slope of the two banks are potentially under slight tension. (2) With the increasing of water level, the sealed transverse joints

Fig. 15. Comparison of the stress of characteristic points at EL 442.2 m inside monoliths obtained from numerical simulation and field monitoring: (a) normal stress in the tangential direction of arch (rtt ) and (b) normal stress in the vertical direction (rvv ).

D. Luo et al. / Engineering Structures 128 (2016) 400–414

413

Fig. 17. Comparison of the increments of aperture variation for transverse joints obtained from numerical simulation and field monitoring: (a) joint #15 and (b) joint #18.

close further. The joints are almost under compression over the entire dam at the normal water level, and the joints are usually not opened. Fig. 17 compares the increments of aperture variation for the transverse joints obtained from numerical simulation and field monitoring when the water level rises from EL 540 to EL 600 m. The deformation distributions of the joints obtained from field monitoring and the ILE model are similar, and the values are considerably close to each other. However, the results obtained from the thin-layer standard element model and structure model without joints are considerably smaller than those obtained from field monitoring and the ILE model. It is found that the ILE model has the ability to provide higher-precision results on the detailed characteristics of high deformation gradients across transverse joints. This is because the use of ILEs provides a good numerical approximation of the weakening of joints by considering the complicatedly and drastically varying displacement fields and material property fields in joint zones. All in all, the comparative analysis results from both the field monitoring and numerical modeling reveal that overall stability of the dam during the filling process meets the dam operational requirements. Meanwhile, the case study on the arch dam presented in this section further justifies the large potential of ILEs in practical engineering applications. The results of ILE model, thin-layer element model and structure model without joints show similar deformation and stress results over the entire structure, but the deformation of transverse joints obtained from the ILE model is in a better agreement with those obtained from field monitoring compared with the other two models. It is demonstrated that the ILE model is capable of providing higher-precision results on the detailed characteristics of high deformation gradients across transverse joints. 5. Conclusions The bonded zone of interfaces/joints is not a simple geometrical plane but a thin and weak interfacial layer with complex material property variations and severe deformation gradients. To model real geometric features and material properties of thin interface/ joint regions without special mesh refinement, an interfacial layer element is developed in this paper for finite element analysis of the structures with nonhomogeneous interfacial layers, like grouted transvers joints. The conclusions of this paper are as follows:

(1) An interfacial layer element can be used to simulate arbitrary material property variations and local high deformation gradients. The inhomogeneity is simulated by a new method combining micro-elements and graded finite elements. The local high deformation gradient across thin interface/joint regions is simulated by new piecewise shape functions. The thin interface/joint regions can be finely simulated in structural analysis by using interfacial layer elements without specially generating fine meshes around interfaces/joints. The interfacial layer elements are characterized by high computational precision and convenience for the modeling of complex structures with interfaces/ joints. (2) By using the interfacial layer elements, the filling process of Xiluodu super-high arch dam is simulated and analyzed. Based on the field monitoring results, it is found that three modeling methods for the grouted transverse joints show similar deformation and stress results over the entire structure, but the interfacial layer element model is capable of providing higher-precision results on the detailed characteristics of high deformation gradients across transverse joints compared with the thin-layer element model and structure model without joints. Hence, it can be concluded that interfacial layer elements have large potential in practical engineering applications. In addition, the comparative analysis results from both field monitoring and numerical modeling reveal that the running status of Xiluodu arch dam is accorded with the performance regulation under water loading and the integrity stability during the filling process meets the safety running requirement.

Acknowledgements This research work was supported by National Natural Science Foundation of China (Nos.: 51339003 and 51579134), National Basic Research Program of China (973 Program, Grant No. 2013CB035902). References [1] Zhang GX, Liu Y, Zheng CY, Feng F. Simulation of influence of multi-defects on long-term working performance of high arch dam. Sci China Technol Sci 2011;54:1–8.

414

D. Luo et al. / Engineering Structures 128 (2016) 400–414

[2] Luo DN, Lin P, Li QB, Zheng D, Liu HY. Effect of the impounding process on the overall stability of a high arch dam: a case study of the Xiluodu dam, China. Arab J Geosci 2015;8:9023–41. [3] Lin P, Ma TH, Liang ZZ, Tang CA, Wang RK. Failure and overall stability analysis on high arch dam based on DFPA code. Eng Fail Anal 2014;45:164–84. [4] Wang CL, Forth JP, Nikitas N, Sarhosis V. Retrofitting of masonry walls by using a mortar joint technique: experiments and numerical validation. Eng Struct 2016;117:58–70. [5] Niwa J, Fakhruddin, Matsumoto K, Sato Y, Yamada M, Yamauchi T. Experimental study on shear behavior of the interface between old and new deck slabs. Eng Struct 2016;126:278–91. [6] Lee CH, Chang KH, Van Do VN. Modeling the high cycle fatigue behavior of Tjoint fillet welds considering weld-induced residual stresses based on continuum damage mechanics. Eng Struct 2016;125:205–16. [7] O’dowd NP, Stout MG, Shih CF. Fracture toughness of alumina-niobium interfaces: experiments and analyses. Philos Mag A 1992;66:1037–64. [8] Rahle M, Evans AG, Ashby MF, Hirth JP. Metal-ceramic interfaces. In: Actascripta metallurgica proceeding series 4. New York: Pergamon Press; 1990. [9] Kang YL, Jia YQ, Qiu Y, Zhang JP. Experimental analysis of mechanical behavior in interfacial minor region of bonded bimaterial. J Exp Mech 1995;10:285–91 (in Chinese). [10] Delale F, Erdogan F. On the mechanical modeling of the interfacial region in bonded half-planes. J Appl Mech 1988;55:317–24. [11] Post D, Wood JD, Han B, Parks VJ, Gerstle FP. Thermal stresses in a bimaterial joint: an experimental analysis. J Appl Mech 1994;61:192–8. [12] Bennegadi ML, Hadjazi K, Sereir Z, Amziane S, Mahi BE. General cohesive zone model for prediction of interfacial stresses induced by intermediate flexural crack of FRP-plated RC beams. Eng Struct 2016;126:147–57. [13] Guin WE, Wang JL. Theoretical model of adhesively bonded single lap joints with functionally graded adherends. Eng Struct 2016;124:316–32. [14] Karabatakis DA, Hatzigogos TN. Analysis of creep behaviour using interface elements. Comput Geotech 2002;29:257–77. [15] Meo M, Thieulot E. Delamination modelling in a double cantilever beam. Compos Struct 2005;71:429–34.

[16] Ranzi G, Gara F, Ansourian P. General method of analysis for composite beams with longitudinal and transverse partial interaction. Compos Struct 2006;84:2373–84. [17] Goodman RE, Taylor RL, Brekke TL. A model for the mechanics of jointed rock. J Soil Mech Foundations Div 1968;94:637–59. [18] Arabshahi H, Lotfi V. Earthquake response of concrete gravity dams including dam–foundation interface nonlinearities. Eng Struct 2008;30:3065–73. [19] Desai CS, Zaman MM, Lightner JG, Siriwardane HJ. Thin - layer element for interfaces and joints. Int J Numer Anal Met 1984;8:19–43. [20] Sharma KG, Desai CS. Analysis and implementation of thin-layer element for interfaces and joints. J Eng Math 1992;118:2442–62. [21] Nam SH, Song HW, Byun KJ, Maekawa K. Seismic analysis of underground reinforced concrete structures considering elasto-plastic interface element with thickness. Eng Struct 2006;28:1122–31. [22] Manzoli OL, Gamino AL, Rodrigues EA, Claro GKS. Modeling of interfaces in two-dimensional problems using solid finite elements with high aspect ratio. Comput Struct 2012;94:70–82. [23] Kaliakin VN, Li J. Insight into deficiencies associated with commonly used zero-thickness interface elements. Comput Geotech 1995;17:225–52. [24] Zhang GA, Zhang JM. Monotonic and cyclic tests of interface between structure and gravelly soil. Soils Found 2006;46:505–18. [25] Pande GN, Sharma KG. On joint/interface elements and associated problems of numerical ill-conditioning. Int J Numer Anal Methods 1979;3:293–300. [26] Zuo Z, Hu Y, Li QB, Liu GW. An extended finite element method for pipeembedded plane thermal analysis. Finite Elem Anal Des 2015;102:52–64. [27] Cavalli F, Gastaldi L. Local enrichment of finite elements for interface problems. Comput Struct 2014;133:111–21. [28] Kim JH, Paulino GH. Isoparametric graded finite elements for nonhomogeneous isotropic and orthotropic materials. J Appl Mech 2002;69:502–14. [29] Santare MH, Lambros J. Use of graded finite elements to model the behavior of nonhomogeneous materials. J Appl Mech 2000;67:819–22. [30] Chen SH, Su PF, Shahrour I. Composite element algorithm for the thermal analysis of mass concrete: simulation of lift joint. Finite Elem Anal Des 2011;47:536–42.