The Effect of Shear Sliding of Vertical Contraction Joints on Seismic Response of High Arch Dams with Fine Finite Element Model

State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, Beijing 100048, China China Institute of Water Resources and Hydropower Research, Beijing 100048, China China (ree Gorges Projects Development Co., Ltd., Beijing 100048, China Construction and Administration Bureau of South-to-North Water Diversion Middle Route Project, Beijing 100038, China North China University of Water Resources and Electric Power, Zhengzhou 450046, China


Introduction
Due to the existence of the thermal stress of mass concrete, arch dams are divided into several dam sections by contraction joints with or without shear keys. Under strong earthquake, the joint opening and sliding may happen. Since the shear keys are used to increase shear resistance, the effect of shear keys should be considered in analysis. Extensive research about the nonlinear analysis of joints has been conducted. Lau et al. [1] proposed a nonlinear model based on plasticity to consider joint opening and sliding based on ADAP88. Arabshahi and Lotfi [2] employed interface elements based on plasticity to consider joint opening and sliding under earthquake. Ahmadi et al. [3] put forward a joint element model considering joint opening and sliding similar to nonassociated plasticity theory. Omid and Lotfi [4] adopted interface elements to simulate the joint nonlinearity, and the constant tangential penalty value is adopted to model the effect of shear keys. e nonlinear model was established based on the framework of plasticity, and fewer degrees of freedom (DOFs) were adopted for analysis in [1][2][3][4]. Fenves et al. [5] presented a joint element to consider the joint opening, which indicated that the high arch tensile stress of arch dams would be released during earthquake. Zhang et al. [6] considered the effect of joint opening and joint reinforcement based on ADAP88. e result demonstrates that the extent of joint opening can be limited. e penalty functions are introduced to guarantee the displacement constraint condition of the interface faces which may have an impact on the convergence of numerical calculation, and the efficiency and stability of nonlinear dynamic analysis need to be improved [1]. Pan et al. [7] used contact surface model to consider joint opening based on the Abaqus software. e tangential sliding is restricted by introducing large tangential stiffness at the joints in [5][6][7]. Jiang et al. [8] adopted soft contact model to consider joint opening. e nonlinear tangential spring is introduced to simulate the tangential movement of an arch-gravity dam during earthquake based on the Abaqus software. e contact surface model used exponential relationship between force and deflection to prevent normal penetration, and the selection of nonlinear tangential stiffness value is complicated and difficult. Du and Jiang [9] considered the shear key's real arrangement to model its nonlinear behavior. Wang et al. [10] considered joint grouting materials with hard and soft contact model of Abaqus software, and shear keys are not considered. Lin and Hu [11] used the nonsmooth equation method to explore the influence of joint opening, tangential sliding, shear key, and reinforcement. Fan et al. [12] combined partitioned finite element with interface element to simulate the seismic failure of dams and seismic sliding of a rock wedge. Guo et al. [13] used a contact force model to consider joint opening. ere is no introduction of the penalty function in [11][12][13].
e study on the parallel solution of dynamic finite element analysis has been widely explored. Belytschko et al. [14] used a Single Instruction Multiple Data (SIMD) to implement the parallel approach on Connection Machine computers. With the advent of message passing libraries, general purpose codes using Single Program Multiple Data (SPMD) approach such as ParaDyn (Hoover et al. [15]) and PRONTO3D (Plimpton et al. [16]) were emerging. Danielson and Namburu [17] presented a code for use on scalable computers using SPMD approach. Parallelization of explicit dynamic finite element was discussed (Fahmy and Namini [18]; Malone and Johnson [19,20]; Namburu et al. [21]; Clinckemaillie et al. [22]) and demonstrated as an enabling technology with explicit message passing (Krysl and Bittnar [23]). Recently, parallel computation has made some progress in the seismic analysis of arch dams. Chen et al. [24] completed the seismic analysis of an arch dam including the opening of contraction joints and dam stress distribution using 6 nodes (12CPUS) with the parallel program generated by the software PFEPG [25]. Zhong et al. [26] studied the seismic damage progress of an arch dam without considering the effect of contraction joints using 4 nodes with the parallel program PDPAD on the basis of PFEPG.
In this paper, the joints with and without shear keys are simplified to be with no-slip condition and with relative sliding condition, respectively. e contact model adopted is discussed with no-slip condition and with relative sliding condition for contraction joints of arch dams considering the manner of independent cantilever dead load type firstly. Secondly, the massed foundation model with viscoelastic artificial boundary and nonuniform ground motion input is explained. irdly, the parallelization strategy employed to increase the computational efficiency is described and the performance of parallel computation is tested. Finally, an existing high arch dam with fine finite element model is analyzed to investigate the effect of shear sliding of vertical joints on seismic response of the arch dam.

Analysis Model
2.1. Contact Model. In this paper, the manner of the independent cantilever dead load type is adopted instead of staged construction and sequence of grouting during construction process of arch dams for simplification. e dead load is assumed to be applied to individual cantilevers without considering transferring force between cantilevers. After joint grouting, water load and other loads such as seasonal temperature load are applied. e equations considering the manner of the independent cantilever dead load type for arch dams including Lagrange multiplier can be written as (1) e detailed information of the above matrix and vector can refer to [13]. e contact force equation can be obtained from (1): where λ l and γ l are the corresponding vectors in the local coordinate system; ; T is the transformation matrix in different coordinate systems. e normal and tangential contact force can be solved iteratively by (2) with constitutive law. In this paper, the strategy of the alternating iterative solution of normal force and tangential force is employed.
e contraction joint has no tensile strength and can only be capable of transferring compressive force in the normal direction. For no-slip condition with shear keys, unlimited force is assumed to be transferred in the tangential direction. For relative sliding condition without shear keys, Coulomb model is suitable for tangential movement.
(a) For no-slip condition with shear keys, the solution of each iterative step is as follows.
Firstly, the normal contact force is obtained as follows: If λ k+1 n ≤ 0, then the joint is in the opening state, and the normal contact forces vanish. e correction is applied as λ k+1 n � 0 and the iterative error directly can be got by err � err + (λ k+1 n − λ k n ) 2 . If λ k+1 n > 0, then the joint is in the closed state, and no correction is required. e iterative error directly can be got by err � err + (λ k+1 n − λ k n ) 2 .

Advances in Civil Engineering
Secondly, the tangential contact forces are obtained as follows: ere is no limit in the tangential force, and no correction is needed. e iterative error can be got directly by For relative sliding condition without shear keys, Coulomb friction criterion is introduced to the iteration process to model sliding of contraction joints without shear keys. e solution of each iterative step is as follows.
Firstly, the normal contact force is obtained as follows: If λ k+1 n ≤ 0, then the joint is in the opening state, and the normal and tangential contact forces vanish. e correction is applied as λ k+1 n � 0, λ k+1 s1 � 0, and λ k+1 s2 � 0 and the iterative error directly can be got by en go to next iterative step.
If λ k+1 n > 0, then the contact surface is in the closed state, and no correction is required. e iterative error directly can be got by err � err + (λ k+1 n − λ k n ) 2 . Secondly, the tangential contact forces are obtained as follows: s ≤ μλ k+1 n , then the joint is in the stick state, and no correction is required. e iterative error can be got directly by err � err n , then the joint is in the sliding state and the scaled tangential forces are defined as λ k+1 s1 � λ k+1 s1 (μλ k+1 n /λ k+1 s ) and λ k+1 s2 � λ k+1 s2 (μλ k+1 n /λ k+1 s ). e iterative error can be got directly by err � err where μ is the friction coefficient. After the iteration is completed, the displacement can be obtained through substituting the contact force into the dynamic equation. Obviously, there is no need to introduce any penalty value for no-slip and relative sliding conditions during the solution process, which avoids the potential stability and uncertainty of numerical calculation caused by parameter selection. e model is special to the nonlinear analysis of arch dams considering the manner of the independent cantilever dead load type.

Mass Foundation Model with Viscoelastic
Artificial Boundary 2.2.1. Viscoelastic Artificial Boundary. e viscoelastic artificial boundary is adopted to consider the effect of radiation damping of far-filed foundation [27][28][29], as shown in Figure 1. Figure 2, for each node of the bottom and side, based on the assumption of one-dimensional wave theory, the free field displacement and velocity time histories are the superposition of the incident wave from bottom to top and the reflected wave from top back to bottom. en, the time histories of the equivalent force on the boundary nodes can be computed using the equation represented in matrix form as

Seismic Input Mechanism. As shown in
where F B is force vector, σ is free field stress tensor computed from free field displacement, and n is outer normal vector of the boundary. C B is the damping coefficient matrix added to the node of boundary, K B is the spring coefficient matrix added to the node of boundary, u is the free field displacement vector, u � u v w T , and _ u the free field velocity vector,

Domain Decomposition Method (DDM).
e basic idea of domain decomposition method (DDM) is to divide a complex computing system into several subsystems by using the strategy of partitioning. e solution of the original system is transformed into the solution of the subsystems, and data exchange is achieved among subsystems through message passing. e DDM consists of overlapping and nonoverlapping domain decomposition method (ODDM and NODDM). e former is based on Schwarz alternation method, and the overlapping region is used to pass message among different partitions. e latter is based on substructure method, and the interface of each partition is used to pass message.
In this paper, the parallel programming based on ODDM is presented. e solution in a whole complex region can be divided into that in several overlapping simple partitions based on the idea of breaking up the whole into parts through Schwarz alternating method, which provides a mathematical basis for distributed parallel computing. e basic principle of Schwarz alternating method is illustrated with a simple static plane stress problem. e plane stress problem in an ABCD region is decomposed into that in two overlapping regions including ABGH and CDFE, and EFHG is the overlapping region shown in Figure 3.
e procedure for the alternative solution method is as follows (Figure 4).  [30], if two nodes of a contact pair are in different subregions, nonoverlapping domain decomposition will appear. To ensure that two nodes of a contact pair are in the same subregion, the above node partitioning based on the Metis program needs to be reprocessed. Two nodes of each contact pair are forced to place in the same subregion by looping over all contact pairs, and the new node partition is obtained. en the internal and overlapping elements of each partition are obtained and the nodes are classified according to the attribute. e specific process is as follows ( Figure 5).

ODDM for Explicit Solution of Wave
As shown in Figure 6, the elements of each partition are divided into internal elements and overlapping elements. Figure 7 shows that the classification of nodes of each partition. e internal nodes and boundary nodes are the nodes to be solved for each partition. e boundary nodes and external nodes are the nodes to exchange information with other partitions. e boundary nodes of this partition are the external nodes of other partitions. e external nodes provide boundary conditions for this partition.

Parallel Computing Framework.
e distributed parallel computing is employed, which is suitable for large-scale cluster. It adopts master-slave programming mode and consists of one main process and several slave processes. e master process is a control program, which does not participate in calculation. It is responsible for sending data to the slave processes and receiving and sorting out the data from the slave processes. e message passing occurs among the master processes and the slave processes, and no message passing occurs among the slave processes. e slave processes are only responsible for the calculation of the corresponding subregions. e message passing is realized by blocking communication between computation and communication as shown in Figure 8. Figure 9 shows the box model and domain decomposition into 95 partitions. e total number of freedom of the model is 10,744,731. e length of input ground motion is 10 s. e performance of parallel computation is tested.

Parallel Performance Test.
From Table 1, running time is from 1313.1 h (54d17.1 h) with 1 slave processor to 40.5 h with 95 slave processors, which shows that the parallel computation software has good speedup. With the increase of processors, the running time is gradually decreasing and speedup is gradually increasing, which shows the parallel computation program has good scalability. With the increase of processors, the computing of each processor is decreasing, and the data communication between processors is increasing. Although  Go through all computing nodes in subregion to find out the nodes sharing the same elements with the computing nodes while they are not currently part of the partition. They are called external nodes, and these elements are overlapping region elements.
Go through all the overlapping area elements of each partition. The nodes in all overlapping area elements are boundary nodes except external nodes.
Go through all the elements of each partition. If the node of the element is composed of internal nodes and boundary nodes, the element is the internal element of the partition. Firstly, initialize the displacement boundary GH in ABGH region, and then the displacement solution u 1 of ABGH region is obtained.
Make EF displacement boundary g v = u 1 in CDFE region, and then the displacement solution v 1 in CDFE region is obtained.
Make GH displacement boundary g u = v 1 in ABGH region, and then the displacement solution in ABGH region is obtained.
Does u n v n converge to the same solution? Advances in Civil Engineering the benefit of parallel computation is increasing, the parallel efficiency is decreasing.

Model Construction.
e dam-foundation finite element model includes 1182559 nodes, 1078026 elements, and 3547677 DOFs, respectively. e maximum dam height is 285.5 m and the crest elevation is 610 m. Figure 10 shows the dam-foundation system and Figure 11 shows the dam model. Figure 12 shows the arrangement of 30 vertical contraction joints which divide the dam into 31 sections (sequent numbering from left to right).
In the present study, the vertical contraction joint is the only source of nonlinearity in the model of the system. e hydrodynamic effect is considered by the modified Westergaard added mass model [31] according to Chinese code.
Vertical contraction joints are as follows: for without shear keys, the frictional coefficient f � 0.8.
e Rayleigh type damping with damping ratio of 5% is selected for analysis.

Load.
e load condition can be seen in Table 2. e seismic wave time histories are shown in Figure 13.

Calculation Scheme
(i) No-slip condition of vertical contraction joints with shear keys: only opening and closing at the joints are considered (ii) Relative sliding condition of vertical contraction joints without shear keys: opening and closing and sliding at the joints may happen

Discussion of Results
(1) Figures 14 and 15 show a comparison of the joint opening under two conditions. Note that the overall distribution of the joint opening is similar in consideration of the two conditions due to the whole load combination including static and seismic loads, although there are some slight differences. Under noslip condition of vertical contraction joints with shear keys, the maximum joint opening value is 12.6 mm which locates at the downstream side of the dam crest for #16 joint. e corresponding time history is shown in Figure 16. In contrast, under relative sliding condition of vertical contraction joints without shear keys, the maximum joint opening value is 13.6 mm which locates at the downstream side of the dam crest for #18 joint. e corresponding time history shown in Figures 17 and  18 shows joint slippage for relative sliding condition of vertical contraction joints without shear keys. e slippage reaches 23.1 mm which locates at the upstream side of the dam crest for #12 joint.   Advances in Civil Engineering Note that the maximum principal tensile stresses around the middle and upper elevations of the upstream surface are dominated by vertical stresses under two kinds of conditions. e values of maximum principal tensile stress and     (2) Figures 19-26 show envelopes of maximum principal tensile stress distribution and vertical stress distribution under two conditions. e results can be obtained as follows.
Note that the maximum principal tensile stress around the middle and upper elevations of the downstream surface is dominated by vertical stress under relative sliding condition, while that under no-slip condition is not dominated by vertical stress. Under relative sliding condition, the maximum principal tensile stress and vertical stress at these    Figure 19: Envelope of maximum principal tensile stress distribution of upstream surface of the dam for no-slip condition [13]. Advances in Civil Engineering 9

Conclusions
(1) Based on the Lagrange multiplier method, a contact model considering the independent cantilever dead load with no-slip condition and relative sliding condition is proposed to model the nonlinearities of vertical contraction joins, which is special to the nonlinear analysis of arch dams considering the manner of dead load type. ere is no need to introduce any penalty value for no-slip and relative sliding conditions, which avoids the potential stability and uncertainty of numerical calculation caused by parameter selection. Different from the conventional Gauss iterative method, the strategy of the alternating iterative solution of normal force and tangential force is employed.        (2) e parallelization based on ODDM and explicit message passing using distributed memory parallel computers is employed to improve the computational efficiency. (3) e effect of shear sliding of vertical joints on seismic response of a high arch dam in China is investigated by nonlinear comprehensive analysis model of fine finite element. e analysis results show that the shear sliding has slight effect on the joint opening characteristics and significantly changes the stress distribution. (4) In summary, the values of maximum principal tensile stress under relative sliding condition are significantly greater than those under no-slip condition. For the former, when the contraction joint opening occurs, the tangential movement is free, which aggravates the vibration of the corresponding dam section in the stream direction and the increase of vertical stress. For the latter, the tangential movement between dam sections is restricted, and the integrity is strengthened, which limits the vibration of the dam in the stream direction. In fact, for high arch dams in China, spherical crown-shaped shear keys are mostly adopted, whose mechanism is more complicated for seismic analysis. e two extreme conditions presented in this paper maybe provide the corresponding envelopes.

Data Availability
All data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.