#### Abstract

The contraction joints of arch dams with and without shear keys are simplified to be with no-slip condition and with relative sliding condition, respectively. Based on the Lagrange multiplier method, a contact model considering the manner of independent cantilever dead load type 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. Different from the conventional Gauss iterative method, the strategy of the alternating iterative solution of normal force and tangential force is employed. The parallelization based on overlapping domain decomposition method (ODDM) and explicit message passing using distributed memory parallel computers is employed to improve the computational efficiency. 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. The result shows that the values of maximum principal tensile stress under relative sliding condition are significantly greater than those under no-slip condition.

#### 1. 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. The nonlinear model was established based on the framework of plasticity, and fewer degrees of freedom (DOFs) were adopted for analysis in [1–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. The result demonstrates that the extent of joint opening can be limited. The 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. The tangential sliding is restricted by introducing large tangential stiffness at the joints in [5–7]. Jiang et al. [8] adopted soft contact model to consider joint opening. The nonlinear tangential spring is introduced to simulate the tangential movement of an arch-gravity dam during earthquake based on the Abaqus software. The 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. There is no introduction of the penalty function in [11–13].

The 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. The 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. Thirdly, 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.

#### 2. 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. The 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. The equations considering the manner of the independent cantilever dead load type for arch dams including Lagrange multiplier can be written as

The detailed information of the above matrix and vector can refer to [13].

The contact force equation can be obtained from (1):where and are the corresponding vectors in the local coordinate system; and ; is the transformation matrix in different coordinate systems.

The 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.

The 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 , then the joint is in the opening state, and the normal contact forces vanish. The correction is applied as and the iterative error directly can be got by .

If , then the joint is in the closed state, and no correction is required. The iterative error directly can be got by .

Secondly, the tangential contact forces are obtained as follows:

There is no limit in the tangential force, and no correction is needed. The iterative error can be got directly by .(b)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. The solution of each iterative step is as follows.

Firstly, the normal contact force is obtained as follows:

If , then the joint is in the opening state, and the normal and tangential contact forces vanish. The correction is applied as , , and and the iterative error directly can be got by . Then go to next iterative step.

If , then the contact surface is in the closed state, and no correction is required. The iterative error directly can be got by .

Secondly, the tangential contact forces are obtained as follows:

Then, .

If , then the joint is in the stick state, and no correction is required. The iterative error can be got directly by .

If , then the joint is in the sliding state and the scaled tangential forces are defined as and . The iterative error can be got directly by , 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. The model is special to the nonlinear analysis of arch dams considering the manner of the independent cantilever dead load type.

##### 2.2. Mass Foundation Model with Viscoelastic Artificial Boundary

###### 2.2.1. Viscoelastic Artificial Boundary

The viscoelastic artificial boundary is adopted to consider the effect of radiation damping of far-filed foundation [27–29], as shown in Figure 1.

###### 2.2.2. Seismic Input Mechanism

As shown in 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. Then, the time histories of the equivalent force on the boundary nodes can be computed using the equation represented in matrix form aswhere is force vector, is free field stress tensor computed from free field displacement, and is outer normal vector of the boundary.

is the damping coefficient matrix added to the node of boundary, is the spring coefficient matrix added to the node of boundary, is the free field displacement vector, , and the free field velocity vector, .

#### 3. Parallelization Strategy

##### 3.1. Domain Decomposition Method (DDM)

The basic idea of domain decomposition method (DDM) is to divide a complex computing system into several subsystems by using the strategy of partitioning. The solution of the original system is transformed into the solution of the subsystems, and data exchange is achieved among subsystems through message passing.

The DDM consists of overlapping and nonoverlapping domain decomposition method (ODDM and NODDM). The former is based on Schwarz alternation method, and the overlapping region is used to pass message among different partitions. The 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. The 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. The basic principle of Schwarz alternating method is illustrated with a simple static plane stress problem. The 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.

The procedure for the alternative solution method is as follows (Figure 4).

##### 3.2. ODDM for Explicit Solution of Wave Equation

For the explicit wave equation in ODDM’s solution, no iteration is needed and the result of the next time step is directly solved based on the last time step owing to the decoupling of equation. The procedure is as follows:(1)ABGH region: based on the displacement and velocity in the ABGH region and the displacement boundary condition on the GH boundary at the last time step, the displacement and velocity in the region at the next time step are solved(2)CDFE region: based on the displacement and velocity in the CDFE region and the displacement boundary condition on the EF boundary at the last time step, the displacement and velocity in the region at the next time step are solved(3)EF boundary: the new displacement boundary condition on the EF boundary is obtained based on the new displacement in the ABGH region(4)GH boundary: the new displacement boundary condition on the GH boundary is obtained based on the new displacement in the CDFE region(5)Then go to next time step

##### 3.3. Domain Decomposition Based on Overlapping Elements including Contact Boundary

After partitioning the finite element nodes based on the Metis program [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. Then the internal and overlapping elements of each partition are obtained and the nodes are classified according to the attribute. The specific process is as follows (Figure 5).

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. The internal nodes and boundary nodes are the nodes to be solved for each partition. The boundary nodes and external nodes are the nodes to exchange information with other partitions. The boundary nodes of this partition are the external nodes of other partitions. The external nodes provide boundary conditions for this partition.

##### 3.4. Parallel Computing Framework

The 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. The 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. The message passing occurs among the master processes and the slave processes, and no message passing occurs among the slave processes. The slave processes are only responsible for the calculation of the corresponding subregions. The message passing is realized by blocking communication between computation and communication as shown in Figure 8.

##### 3.5. Parallel Performance Test

Figure 9 shows the box model and domain decomposition into 95 partitions. The total number of freedom of the model is 10,744,731. The length of input ground motion is 10 s. The performance of parallel computation is tested.

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 the benefit of parallel computation is increasing, the parallel efficiency is decreasing.

#### 4. Dynamic Analysis

##### 4.1. Model Construction

The dam-foundation finite element model includes 1182559 nodes, 1078026 elements, and 3547677 DOFs, respectively. The 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. The hydrodynamic effect is considered by the modified Westergaard added mass model [31] according to Chinese code.

##### 4.2. Material Parameters

The concrete material parameters are as follows: density ; static elastic modulus ; dynamic elastic modulus ; Poisson’s ratio .

The foundation rock material parameters are as follows: density ; deformation modulus and Poisson’s ratio are given according to different elevations [13].

Vertical contraction joints are as follows: for without shear keys, the frictional coefficient .

The Rayleigh type damping with damping ratio of 5% is selected for analysis.

##### 4.3. Load

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

**(a)**

**(b)**

**(c)**

##### 4.4. 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

##### 4.5. 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 no-slip 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. The 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. The corresponding time history shown in Figures 17 and 18 shows joint slippage for relative sliding condition of vertical contraction joints without shear keys. The slippage reaches 23.1 mm which locates at the upstream side of the dam crest for #12 joint.

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. The values of maximum principal tensile stress and vertical stress under relative sliding condition are significantly greater than those under no-slip condition. Under no-slip condition, the maximum principal tensile stress and vertical stress at these locations are 2.7 MPa and 2.3 MPa, respectively, while those under relative sliding condition are 6.5 MPa and 5.6 MPa.

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 locations are 4.2 MPa and 3.8 MPa, respectively. From Table 3, it is concluded that under no-slip condition the maximum principal tensile stress at these locations of the dam is 3.35 MPa, which is dominated by shear stress sxz of 3.16 MPa.

#### 5. 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. There 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)The parallelization based on ODDM and explicit message passing using distributed memory parallel computers is employed to improve the computational efficiency.(3)The 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. The 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. The 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

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This paper was supported by the National Key Research and Development Program of China (no. 2017YFC0404903), National Natural Science Foundation of China (no. 51709283), and China Three Gorges Corporation Research Project (contract no. XLD/2115).