A Modified B-Differentiable Equation Method for the Seismic Analysis of Arch Dams Considering the Initial Strength of Contraction Joints

Abstract

For arch dams, the joint surface has a certain bond strength after grouting the contraction joints, which can withstand the arch-wise tensile stress to a certain extent and influence the stress distribution of the dam blocks on both sides. The seismic response analysis of arch dams needs to consider the influence of the initial tensile and shear strengths generated by the contraction joint grouting. Thus, a modified B-differentiable equation method is proposed by introducing the initial tensile strength and the initial shear strength of contraction joints the into traditional B-differentiable equation method. In the proposed method, the shear strength varies with the contraction joint opening. The modified B-differentiable equation method can still be solved by the B-differentiable damped Newton method with theoretical guarantee of convergence. Then, a seismic calculation model for the dam–reservoir–foundation system is developed based on the modified B-differentiable equation method, the Westergaard additional mass method, and multiple transmission boundary conditions. The influence of the initial tensile and shear strength of the contraction joints on the arch dam seismic response is discussed. The results show that the initial tensile and shear strength of the contraction joints have an influence on the opening and distribution range of the contraction joints and the maximum values of the principal tensile and compressive stress in the dam body. The initial tensile and compressive strength of the contraction joints should be considered when carrying out seismic response analysis of arch dams.

1. Introduction

During the construction work of arch dams, in order to facilitate heat dissipation and reduce the shrinkage stress of the concrete to prevent the dam from cracking, the dam needs to be poured in sections with contraction joints between each section, after the temperature of the dam has cooled to a stable temperature, the dam is then pressure-grouted to make the dam an integral whole [1]. In order to improve the shear resistance of the seam surface, the keyway is usually provided in the seam surface. In practice, the contraction joints are grouted and provided with keyways, but their tensile and shear strength is also limited. Under seismic action, the contraction joints of arch dams may open, close, and misalign, thus causing redistribution of stress in the arch dam body, which has a considerable impact on the seismic analysis and safety evaluation of arch dams [2,3,4].

There are currently three main types of methods for modelling contraction joint contact problems in arch dams: the joint element method, the crack smeared method, and the contact force model. The joint element method is actually a way of simulating the opening, closing, and misalignment of contraction joints by the constitutive relationship of the joint elements, thus transforming the contact nonlinearity into a material nonlinearity problem [5]. The joint element method includes nonlinear spring element and interface element. The nonlinear spring element [6,7] consists of the corresponding nodes in the contraction joints, and it uses a translational spring and a rotational spring to represent the relationship between the pressure and bending moment acting at the contact pair and the relative displacement and relative rotation of the contact pair, respectively. The interface element [8,9,10,11,12,13] consists of two faces that overlap each other or have a small gap, each defined by nodes that are not in the same plane. The element can simulate the initial tensile strength of the contraction joints, the initial gap, and the nonlinear constitutive relationship between normal pressure and tangential friction, but mutual embedding may occur between the contacting bodies. Meanwhile, to ensure the convergence of the algorithm and the accuracy of the results, the spring value or the normal and tangential stiffness coefficients of the coupling units need to be selected reasonably [14]. The crack smeared method [15] simulates the opening, misalignment, and closing of cracks and contraction joints in the dam body by modifying the material properties at the integration points in the conventional elements, which can achieve relatively good approximation results under certain conditions [16].

The contact force model takes the contact force between the contraction joints as an unknown quantity and solves for the contact force that satisfies the constraint by taking the frictional contact condition as the constraint [17,18,19]. Zhao et al. [20] applied the finite element hybrid method for the numerical analysis of arch dam contraction joints. The method decomposes the force acting on the contact body into the external force and the contact force on the contact surface, treats the displacement of the contact body and the contact force on the contact interface as hybrid variables, and iterates the contact state and the contact force simultaneously. Du et al. [21] used the nonlinear dynamic contact constitutive model and Coulomb friction model to simulate the normal interactions and tangential interactions of contraction joints. Sun [22] adopted a surface-to-surface contact model to investigate the effectivity of the shape memory alloy restrainer bars to reduce the seismic response of arch dams. Mirzabozorg [23] studied the effects of wave passage and incoherency on the seismic response of arch dams by using a node-to-node contact model to simulate the contraction joints. Guo et al. [24,25] adopted the contact nonlinear model based on the contact force method to perform parallel analysis of the seismic response of arch dams. Xue et al. [20] expressed the frictional contact conditions between the contraction joints in the form of the B-differentiable equation [26,27] and performed a comparative analysis of the seismic response of an arch dam under the action of ground shaking with different canonical spectra, but the study did not consider the effect of the initial strength of the contraction joints on the dam seismic response. In general, the joint element method can be computationally expensive, especially when there is a large number of joints or when the joints have complex geometries. The crack smeared method has a lower computational cost compared to the joint element method for large-scale problems. The computational cost of the contact force model depends on the complexity of the contact algorithm and the number of contact pairs.

A dynamic calculation model of the dam–reservoir–foundation system considering the effect of the contraction joints will be developed based on a modified B-differentiable equation method, the Westergaard additional mass method, and multiple transmission boundary conditions to discuss the effect of the initial tensile and compressive strength of the contraction joints on the seismic response of the arch dam. The results of the study will provide a theoretical basis and suggestions for the seismic safety evaluation of arch dams.

2. The Modified B-Differentiable Equation Method

Under seismic action, the contraction joints in an arch dam may undergo opening, closing, and frictional misalignment, which can be solved as a contact nonlinear problem. Under the assumption that the seam surface is planar, the B-differentiable equation method writes the frictional contact conditions between the contraction joints in the form of a B-differentiable equation so that they are strictly satisfied [28]. The following is a brief introduction of the B-differentiable equation method [27] for the 3D elastic frictional contact problem, taking the contact of two objects (represented by body 1 and body 2, respectively) as an example.

Based on the assumption of small deformation and small strain, a point-to-point contact model is applied to the contact surface after the finite element discretization. That is, the nodes of two contact surfaces are in one-to-one correspondence along the normal direction of the contact surface, constituting multiple contact point pairs. The contact variables and contact conditions along the normal and tangential directions of the contact surface are often used to describe the contact problem. Therefore, a local coordinate system nab is defined at the contact point pair, where $\overrightarrow{n}$ denotes the normal direction of the contact surface, pointing from body 2 to body 1, and $\overrightarrow{a}$ and $\overrightarrow{b}$ denote the two tangential unit vectors at this contact point that are perpendicular to each other on the contact surface. $P_n$, $P_a$, and $P_b$ denote the contact forces of the contact point pairs in the local coordinate system. Assuming that the total number of contact point pairs in the contact system is NC, for the three-dimensional elastic friction contact problem, the contact conditions at each contact pair i can be expressed in the form of the following set of B-differentiable equations:

$${\mathrm{H}}_2^i\left(d{u}_c^i,d{P}^i\right)=\min \left\{r\Delta u_n^i,P_n^i\right\}=0$$

(1)

$${\mathrm{H}}_3^i\left(d{u}_c^i,d{P}^i\right)=P_a^i-\lambda P_a^i\left(r\right)=0$$

(2)

$${\mathrm{H}}_4^i\left(d{u}_c^i,d{P}^i\right)=P_b^i-\lambda P_b^i\left(r\right)=0$$

(3)

where

$$P_a^i(r)=P_a^i-r\Delta d{u}_a^i$$

(4)

$$P_b^i(r)=P_b^i-r\Delta d{u}_b^i$$

(5)

$$\lambda =\min \left\{{\displaystyle \frac{\mu \min \left\{P_n^i,0\right\}}{\sqrt{{\left(P_a^i\left(r\right)\right)}^2+{\left(P_b^i\left(r\right)\right)}^2}}},1\right\}$$

(6)

where ${\mathrm{H}}_2^i, {\mathrm{H}}_3^i, {\mathrm{H}}_4^i$ denote the functions of the normal contact condition and tangential contact condition. $d{\mathbf{u}}_{\mathrm{c}}^i, d{\mathbf{P}}^i$ denote the incremental displacement and incremental contact force of the ith contact point pair, respectively. μ is the friction coefficient; $\Delta u_n^i, \Delta d{u}_a^i, \Delta d{u}_b^i$ denote the incremental normal relative displacement and tangential relative displacement of the ith contact point pair, respectively. r takes a positive value, which is equivalent to a penalty factor, but unlike the penalty factor of the penalty function method. The contact condition can be satisfied accurately without taking a large value of r [21]. $P_n^i, P_a^i, P_b^i$ denote the normal contact force and two tangential contact forces in the ith contact point pair, respectively.

Equation (1) represents the normal contact condition, including the conditions that the normal direction cannot be embedded in each other and the normal contact force is pressure, assuming that a positive value of the contact force indicates pressure and a negative value indicates tension; a positive value of the normal relative displacement indicates that the contact point pair is open and a negative value indicates that embedding occurs. Equations (2) and (3) represent the tangential contact condition, assuming that the tangential motion obeys the Coulomb friction contact condition, that is, when sliding, the magnitude of the tangential contact force is equal to the normal contact force multiplied by the friction coefficient, the direction of sliding and tangential contact force direction in a straight line and in the opposite direction. In the absence of relative sliding, the absolute value of the tangential contact force is less than the normal contact force multiplied by the friction coefficient $\mu P_n^i$. The functions ${\mathrm{H}}_2^i, {\mathrm{H}}_3^i, {\mathrm{H}}_4^i$ in Equations (1)–(3) are not differentiable, which is caused by operator min.

For the actual arch dam project, the joint surface also has a certain bond strength after grouting the contraction joints, which can withstand the arch-wise tensile stress to a certain extent, and influence the stress distribution of the dam blocks on both sides. On the other hand, the application of trapezoidal or spherical keyways at the contraction joints can provide a certain shear strength to restrain the tangential movement between the dam blocks and improve the integrality of the dam. For trapezoidal or spherical keyway, as long as the opening degree is less than the height of the keyway when the contraction joint is opened, although relative sliding may occur along the joint surface on both sides of the dam block, the keyway can still exert shear resistance. The shear strength is not a constant, but it is related to the opening degree of the contraction joints. In order to consider that it has a certain initial tensile strength in the normal direction, and the strength and deformation characteristics that the restriction of the keyway on the tangential relative movement of the dam blocks on both sides of the contraction joints varies with the opening, the B-differentiable contact equations for solving frictional contact problems are modified. That is, the initial tensile strength and the shear strength varying with the opening of the contraction joints are introduced in Equations (1) and (6), respectively:

$${\mathrm{H}}_2^i\left(d{u}_c^i,d{P}^i\right)=\min \left\{r\Delta u_n^i,P_n^i+{\widehat{P}}_n^i\right\}=0$$

(7)

$$\lambda =\min \left\{{\displaystyle \frac{\mu \min \left\{P_n^i,0\right\}+c_T\left(\Delta u_n^i\right)\times {\widehat{F}}_t^i}{\sqrt{{\left(P_a^i\left(r\right)\right)}^2+{\left(P_b^i\left(r\right)\right)}^2}}},1\right\}$$

(8)

where ${\widehat{P}}_n^i$ represents the normal initial tensile nodal force caused by the bond between the cement paste and the dam concrete after grouting the contraction joints at the ith contact point pair, which is the normal cohesive force, and its value is equal to the product of the tensile strength of the mortar and the area assumed by the current contact nodes, which becomes zero once the contraction joint is pulled apart and the normal contact condition is restored to (1). $c_T\left(\Delta u_n^i\right)\times {\widehat{F}}_t^i$ represents the tangential binding force, which is tangential cohesion, caused by the bonding of the keyway and the cement slurry to the concrete of the dam to the tangential movement of the dam blocks on both sides of the contraction joint, related to the opening of the contraction joint. ${\widehat{F}}_t^i$ is the tangential shear nodal force generated by the keyway when the contraction joint is closed, and its value is equal to the product of the shear strength of the concrete and the area assumed by the current contact nodes. $c_T\left(\Delta u_n^i\right)$ represents the tangential restraint reduction factor of the keyway and is related to the contraction joint opening as defined below:

$$c_T\left(\Delta u_n^i\right)=\left\{\begin{array}{c}1\\ 1-\Delta u_n^i/\mathrm{Hk}\\ 0\end{array}\right.\begin{array}{c}\Delta u_n^i=0\\ 0<\Delta u_n^i\le \mathrm{Hk}\\ \Delta u_n^i>\mathrm{Hk}\end{array}$$

(9)

where Hk is the height of keyway; when the contraction joint is fully closed, the value of $c_T$ is 1, when the contraction joint opening exceeds the height of keyway, the value is 0. When the contraction joint opening is between 0 and Hk, it is assumed that the coefficient changes linearly with the contraction joint opening. Other physical quantities in Equations (7)–(9) have the same meaning as in Equations (1)–(6).

The modified B-differentiable equation group can still be solved by the B-differentiable damped Newton method [28,29].

3. The FE Model for the Arch Dam–Reservoir–Foundation System

3.1. Computational Method

(a)

Dam–Infinite Foundation Interaction

For the calculation model of the dam–infinite foundation interactions, the radiation damping effect of infinite foundation must be considered. In this paper, the artificial boundary condition expressed by the multiple transmission formula proposed by Liao [30] is adopted to simulate dam–foundation interactions. The boundary condition is based on a plane wave assumption, and is applicable to the condition that the far field is homogeneous and the foundation is infinite. It is space-time decoupling, which means that the displacement and velocity of the artificial boundary point are only related to a few adjacent nodes along the normal direction of the boundary and the adjacent time at the current moment, so the calculation amount is greatly reduced.

The artificial boundary conditions are imposed at the boundaries of the computation model, with the purpose of eliminating or reducing the impact of the reflected waves at artificial boundaries. It can not only simulate the propagation process of seismic waves towards infinity, but also, the seismic input mechanism adapted to the artificial boundaries can reflect the non-uniform distribution of seismic waves along the valley.

(b)

Dam–Reservoir Water Interaction

The dynamic water added mass method is applied for reflecting the effect on dam–reservoir water dynamic interactions. The dynamic water added mass is calculated following the Westergaard formula [31].

$$m_w\left(h\right)={\displaystyle \frac{7}{8}}\rho \sqrt{H_{0}h}$$

(10)

where $m_w\left(h\right)$ represents the additional mass at the node with depth h, H₀ represents the depth of reservoir water, and ρ indicates the density of water.

(c)

Method for Solving the Dynamic Equilibrium Equations of Dam–Reservoir–Foundation System

At the moment t + dt, the dynamic equilibrium equation of the arch dam–reservoir–foundation system is:

$$\mathbf{M}{\ddot{\mathbf{u}}}^{t+dt}+C{\dot{u}}^{t+dt}+K{u}^{t+dt}=F^{t+dt}+P_c^{t+dt}$$

(11)

where K, C, and M are the stiffness, damping, and mass matrices of the structure, respectively. $F^{t+dt}$ and $P_c^{t+dt}$ represent the external load on the structure and the contact force on the contact surface at the moment t + dt, respectively. $u^{t+dt}$, ${\dot{u}}^{t+dt}$, and ${\ddot{u}}^{t+dt}$ represent the displacement, velocity, and acceleration at the moment t + dt, respectively.

An implicit–explicit integration method is used to discretize the dynamic equilibrium equations in the time domain. In the implicit–explicit integration method, the elements and nodes are divided into two kinds: implicit and explicit. The Newmark integration method and forecast-correction integration method are used for implicit and explicit integration, respectively. In the forecast-correction explicit integration method, the displacement and velocity at the moment t + dt can be assumed as the following form:

$$\begin{array}{l}{\mathbf{u}}^{t+dt}={\tilde{\mathbf{u}}}^{t+dt}+d{t}^2\cdot \eta \cdot {\ddot{\mathbf{u}}}^{t+dt}\hfill \\ {\dot{\mathbf{u}}}^{t+dt}={\widetilde{\dot{\mathbf{u}}}}^{t+dt}+dt\cdot \gamma \cdot {\ddot{\mathbf{u}}}^{t+dt}\hfill \\ {\tilde{\mathbf{u}}}^{t+dt}={\mathbf{u}}^t+dt\cdot {\dot{\mathbf{u}}}^t+d{t}^2\cdot (1-2\eta ){\ddot{\mathbf{u}}}^t/2\hfill \\ {\widetilde{\dot{\mathbf{u}}}}^{t+dt}={\dot{\mathbf{u}}}^t+dt\cdot (1-\gamma ){\ddot{\mathbf{u}}}^t\hfill \end{array}$$

(12)

where ${\tilde{u}}^{t+dt}$ and ${\tilde{\dot{u}}}^{t+dt}$ represent the predicted displacement and velocity at the moment t + dt, respectively. $\gamma$ and $\eta$ represent the integration constant. ${\ddot{u}}^{t+dt}$ is an unknown variable, and ${\tilde{u}}^{t+dt}$ and ${\tilde{\dot{u}}}^{t+dt}$ can be obtained from the known displacement, velocity, and acceleration at moment t.

According to the implicit and explicit degree of freedom, the dynamic equilibrium equation can be expressed as:

$$\left[\begin{array}{cc}{M}_i& \\ & M_e\end{array}\right]\left[\begin{array}{c}{\ddot{u}}_i^{t+dt}\\ {\ddot{u}}_e^{t+dt}\end{array}\right]+\left[\begin{array}{cc}{C}_{ii}& C_{ie}\\ C_{ei}& C_{ee}\end{array}\right]\left[\begin{array}{c}{\dot{u}}_i^{t+dt}\\ {\tilde{\dot{u}}}_e^{t+dt}\end{array}\right]+\left[\begin{array}{cc}{K}_{ii}& K_{ie}\\ K_{ei}& K_{ee}\end{array}\right]\left[\begin{array}{c}{u}_i^{t+dt}\\ {\tilde{u}}_e^{t+dt}\end{array}\right]=\left[\begin{array}{c}{F}_i^{t+dt}+P_c^{t+dt}\\ F_e^{t+dt}\end{array}\right]$$

(13)

Equation (13) can be rewritten as two sets of equilibrium equations:

$$M_i{\ddot{u}}_i^{t+dt}+C_{ii}{\dot{u}}_i^{t+dt}+K_{ii}{u}_i^{t+dt}=F_i^{t+dt}+P_c^{t+dt}-C_{ie}{\tilde{\dot{u}}}_e^{t+dt}-K_{ie}{\tilde{u}}_e^{t+dt}$$

(14)

$$M_e{\ddot{u}}_e^{t+dt}=F_e^{t+dt}-C_{ee}{\tilde{\dot{u}}}_e^{t+dt}-K_{ee}{\tilde{u}}_e^{t+dt}-C_{ei}{\dot{u}}_i^{t+dt}-K_{ei}{u}_i^{t+dt}$$

(15)

Then, the solution process of the dynamic equilibrium equations is converted to solve for Equations (14) and (15), alternately.

Using the Newmark integration method, one can discretize Equation (14) in the time domain:

$${\overline{K}}_{ii}{u}_i^{t+dt}={\overline{F}}_i^{t+dt}+P_c^{t+dt}-C_{ie}{\tilde{\dot{u}}}_e^{t+dt}-K_{ie}{\tilde{u}}_e^{t+dt}$$

(16)

with

$$\begin{array}{l}{\overline{K}}_{ii}=K_{ii}+\left(\gamma /{\displaystyle \raisebox{1ex}{$dt$}\!\left/ \!\raisebox{-1ex}{$\eta $}\right.}\right)C_{ii}+{\left(d{t}^2\eta \right)}^{-1}{M}_i\\ {\overline{F}}_i^{t+dt}=F_i^{t+dt}-F_i^t+\left(0.5d{t}^2{\eta }^{-1}{\ddot{u}}_i^t+d{t}^{-1}{\eta }^{-1}{\dot{u}}_i^t\right)M_i-\\ \left(\left(1-0.5\gamma {\eta }^{-1}\right)dt{\ddot{u}}_i^t-\gamma {\eta }^{-1}{\dot{u}}_i^t\right)C_{ii}\end{array}$$

(17)

By combining Equation (16) with Equations (2), (3), (7) and (8), the dynamic equilibrium equations at the moment t + dt can be written as B-differentiable equations. Then, the unknown variables $u_i^{t+dt}$ and $P_c^{t+dt}$ can be obtained by the B-differentiable damped Newton method [27,28]. Then, the unknown variables ${\ddot{u}}_e^{t+dt}$ can be obtained by Equation (15).

3.2. The FE Mesh of the Arch Dam–Foundation

A high arch dam located in northwest China is analyzed as an example [26]. The dam height is 240 m, the dam crest elevation is 990 m, and the valley bottom elevation is 750 m. The scope of the finite element model of the arch dam–foundation system is to extend by one time the height of the dam in the horizontal, vertical, and downstream directions, respectively. The finite element mesh dispersion of the arch dam–foundation system and the dam body are shown in Figure 1, and Figure 2 is the location diagram of the contraction joints of the arch dam. In order to accurately simulate the propagation of elastic waves, the size of the finite element mesh of the foundation is not more than 20 m in the vertical direction. The finite element model of the dam–foundation system uses a three-dimensional eight-node solid element, which is divided into three layers of grids along the thickness direction of the dam body, and a total of 34 contraction joints are simulated. The dam foundation system finite element model has 162,428 nodes and 148,695 elements, including 4518 total dam elements, 8596 total nodes, and 2288 contact node pairs.

In the calculation and analysis, it is assumed that the dam concrete and bedrock are isotropic linear materials, and the values of their physical and mechanical parameters are shown in

Table 1. Material parameters.

Material	Density (kg/m3)	Elastic Modulus (GPa)	Poisson’s Ratio	Coefficient of Linear Expansion (/°C)	Thermal Conductivity (W/m·°C)	Specific Heat Capacity (kJ/(kg·°C))
Concrete	2400	21	0.167	1.0 × 10−5	3.0	970
Bedrock I	2755	13.5	0.24	1.0 × 10−5	2.67	840
Bedrock II	2700	10.0	0.26	1.0 × 10−5	2.67	840

. Bedrock I is the red area in Figure 1, and bedrock II is the purple area in Figure 1.

First, the static analysis of the dam–foundation system is carried out considering the hydrostatic pressure, dam weight, sediment pressure, and temperature. The water level in front of the dam is 150 m and the water level behind the dam is 41.5 m. The depth, floating bulk density, and internal friction angle of the silty sand are 57.50 m, 9.0 KN/m³, and 12.0°, respectively. The sediment pressure is applied to the upstream face of the dam in accordance with the hydrostatic pressure [32]. The mass of the foundation is neglected and the stiffness of the foundation is considered. On the basis of the static analysis of the dam–foundation system, the dynamic analysis is performed by applying seismic loads. In the dynamic analysis, the influence of sediment pressure on the dam–foundation system is ignored, Rayleigh damping is applied to reflect the damping characteristics of the dam–foundation system, and the Rayleigh damping coefficient is taken as 5%. The mass and stiffness of the foundation are both considered. Assuming that the seismic wave is incident vertically, the seismic action effect will consider the combined effect of both horizontal earthquake (downstream and cross-river) and vertical earthquake. The normalized acceleration time course of the input ground motion is shown in Figure 3, and the peak ground motion acceleration is 0.357 g. The seismic exciting records are applied on the front, back, left, right, and bottom boundaries of the arch dam–foundation system.

3.3. The Parameters of Contraction Joints

In practical engineering, the normal tensile strength and tangential shear strength should be related to factors such as contraction joint opening before grouting and grouting construction quality, for which information is still lacking. In order to compare the effects of contraction joint normal tensile strength and contraction joint tangential shear strength reflecting the action of the keyway on the stress and deformation of the dam body, the normal tensile strength and tangential shear strength of the contraction joint are set to different strengths, as shown in

Table 2. Calculation model of contraction joints.

Model Number	Normal Tensile Strength (MPa)	Tangential Shear Strength (MPa)
CM1	1.0	3.0
CM2	0.5	1.0
CM3	0.0	1.0
CM4	0.0	0.0

, when the deformation of the contraction joint is simulated by adopting the modified B-differentiable equation group method.

In Table 2, for contact model CM1, the contraction joint normal tensile capability and the keyway shear capability are stronger; for contact model CM2, the contraction joint normal tensile capability and the keyway shear capability are slightly reduced. The normal tensile strength of the contraction joint is not considered in contact models CM3 and CM4. In contact model CM3, although the contraction joint normal tensile strength is not considered, the shear strength of the keyway is 1 MPa, so after the contraction joint is opened, as long as the contraction joint opening is smaller than the height of the keyway, the keyway can still provide some shear capacity, only that this shear capacity decreases with the increase of the contraction joint opening, refer to Equation (9). In contact model CM4, the contraction joint normal tensile strength and the shear strength of the keyway are not considered, and the dam blocks on both sides can move freely along the tangential direction when the contraction joint is opened. Notably, the normal tensile strength is taken with reference to the tensile strength of the cement paste and the tensile strength of the bond between the cement paste and the dam concrete, and the tangential shear strength is taken mainly based on the concrete shear strength.

4. Numerical Results

4.1. Contraction Joint Aperture

The initial tensile strength is set to 1 MPa when the contraction joint model CM1 is used, indicating good construction quality and strong bond strength after contraction joint grouting. As seen in Figure 4 of the contraction joint opening distribution, the opening of the contraction joints and the opening range showed more obvious spacing phenomena, and the opening of the contraction joints near the arch crown beam is larger. The main reason for these phenomena is that during the dam movement, the values of the tensile stress applied to each contraction joint vary in magnitude, and the tensile stress at some contraction joints may firstly exceed the initial tensile strength, and the contraction joint opens. According to the assumptions of the contraction joint calculation model, once the contraction joint is opened, the location no longer bears the tensile stress during the subsequent dynamic deformation of the dam, which will produce the following two phenomena. On the one hand, the tensile stress acting on these contraction joints is redistributed under the seismic load, making the tensile stress in the rest of the contraction joints increase and possibly exceed the initial tensile strength, making the contraction joint opening range further expand; on the other hand, the release of tensile stress will also make the tensile stress in the adjacent contraction joints always less than the tensile strength, and the opening degree is small or even impossible, thus making it difficult to further expand the opening range of these contraction joints. When the initial tensile strength of the contractions joint gradually decreases, that is, using the contraction joint models CM2, CM3, and CM4, it can be seen from Figure 5, Figure 6 and Figure 7 that as the initial tensile strength decreases, in the top part of the dam, the range of opening along the cross-river to the contraction joint increases. When the initial tensile strength is 0 (contraction joint models CM3 and CM4), all contraction joints are opened, and the opening degree of other contraction joints tends to be uniform, except for the contraction joints near the two dam heads.

(1)

Contraction joint model CM1.

Figure 4. CM1 distribution of maximum opening of contraction joints on upstream and downstream surfaces.

Figure 4. CM1 distribution of maximum opening of contraction joints on upstream and downstream surfaces.

(2)

Contraction joint model CM2.

Figure 5. CM2 distribution of maximum opening of contraction joints on upstream and downstream surfaces.

Figure 5. CM2 distribution of maximum opening of contraction joints on upstream and downstream surfaces.

(3)

Contraction joint model CM3.

Figure 6. CM3 distribution of maximum opening of contraction joints on upstream and downstream surfaces.

Figure 6. CM3 distribution of maximum opening of contraction joints on upstream and downstream surfaces.

(4)

Contraction joint model CM4.

Figure 7. CM4 distribution of maximum opening of contraction joints on upstream and downstream surfaces.

Figure 7. CM4 distribution of maximum opening of contraction joints on upstream and downstream surfaces.

In addition, when the contraction joint model CM1 is adopted, the contraction joint opening near the arch crown beam of the dam body is significantly higher than that at other parts. This is because, under seismic action, when the arch dam moves upstream along the river, a large arch tensile stress will be generated in the upper and middle elevations of the dam body, near the arch crown beam, which may exceed the initial tensile strength of the contraction joints in this area, causing the contraction joints to open, and leading to a reduction in the stiffness of the dam body, and in the subsequent deformation, it is easy to further expand the opening range. When the contraction joint models CM2, CM3, and CM4 are adopted, the initial tensile strength gradually decreases. Then, each contraction joint can only bear a small tensile stress or even no tensile stress, which makes each contraction joint easier to open, the opening degree tends to be uniform, and the maximum value has a relatively large reduction.

4.2. Stress Distribution of the Dam Body

It can be seen from Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13 that the distribution of tensile and compressive principal stress on the dam surface obtained by adopting the contraction joint models CM1, CM2, and CM3 is basically the same. As shown in

Table 3. The maximum tensile and compressive principal stress (MPa).

	MPS	TS-UF	TS-DF	CS-UF	CS-DF
Contact Model
CM1		2.74	3.82	7.40	5.89
CM2		2.61	3.20	7.30	6.49
CM3		2.59	3.09	7.30	6.50
CM4		2.53	3.83	7.66	7.45

Notes: MPS means maximum principal stress; TS means tensile stress; CS means compressive stress; UF means upstream surface; DF means downstream surface.

, the maximum values of tensile and compressive principal stress decrease with the decrease of the initial tensile strength and shear strength of the contraction joint. As shown in Figure 14 and Figure 15, compared with the first three models, the stress response obtained by adopting model CM4 shows obvious discontinuity in stress distribution and slight increase in stress value. This is because if the initial tensile strength and initial shear strength of contraction joints are not taken into account, relative movement between dam blocks occurs easily under the action of seismic reciprocating loads, resulting in the contraction joints opening more easily, which has a certain impact on the stress distribution of the dam.

(1)

Contact model CM1.

Figure 8. CM1 distribution of maximum tensile principal stress P1 of upstream and downstream dam surfaces (Pa).

Figure 8. CM1 distribution of maximum tensile principal stress P1 of upstream and downstream dam surfaces (Pa).

Figure 9. CM1 distribution of maximum compressive principal stress P3 of upstream and downstream dam surfaces (Pa).

Figure 9. CM1 distribution of maximum compressive principal stress P3 of upstream and downstream dam surfaces (Pa).

(2)

Contact model CM2.

Figure 10. CM2 distribution of maximum tensile principal stress P1 of upstream and downstream dam surfaces (Pa).

Figure 10. CM2 distribution of maximum tensile principal stress P1 of upstream and downstream dam surfaces (Pa).

Figure 11. CM2 distribution of maximum compressive principal stress P3 of upstream and downstream dam surfaces (Pa).

Figure 11. CM2 distribution of maximum compressive principal stress P3 of upstream and downstream dam surfaces (Pa).

(3)

Contact model CM3.

Figure 12. CM3 distribution of maximum tensile principal stress P1 of upstream and downstream dam surfaces (Pa).

Figure 12. CM3 distribution of maximum tensile principal stress P1 of upstream and downstream dam surfaces (Pa).

Figure 13. CM3 distribution of maximum compressive principal stress P3 of upstream and downstream dam surfaces (Pa).

Figure 13. CM3 distribution of maximum compressive principal stress P3 of upstream and downstream dam surfaces (Pa).

(4)

Contact model CM4.

Figure 14. CM4 distribution of maximum tensile principal stress P1 of upstream and downstream dam surfaces (Pa).

Figure 14. CM4 distribution of maximum tensile principal stress P1 of upstream and downstream dam surfaces (Pa).

Figure 15. CM4 distribution of maximum compressive principal stress P3 of upstream and downstream dam surfaces (Pa).

Figure 15. CM4 distribution of maximum compressive principal stress P3 of upstream and downstream dam surfaces (Pa).

5. Conclusions

In this paper, a modified B-differentiable equation method was proposed to simulate the initial tensile strength and shear strength of the contraction joints of an arch dam. Based on the modified B-differentiable equation method, Westergaard additional mass method, and multiple transmission boundary conditions, a seismic calculation model of a dam–reservoir–foundation system considering the influence of contraction joints is established. The influence of initial tensile and compressive strength of contraction joints on the seismic response of arch dams is discussed. When the initial tensile strength and the initial shear strength of the contraction joint are strong, the opening degree and opening range of the contraction joints appear as obvious spacing phenomenon. As the initial tensile strength and shear strength of the contraction joints gradually decreases, each contraction joint can only bear a small tensile stress or even no tensile stress, which makes each contraction joint easier to open, the opening degree tends to be uniform, and the maximum value has a relatively large reduction. The maximum values of tensile and compressive principal stress of arch dams decrease with the decrease of initial tensile strength and initial shear strength of contraction joints. The reduction rates of maximum values of tensile and compressive principal stress are 32.2% and 12.2%, respectively. When the initial tensile strength and initial shear strength of the contraction joints are not considered, the maximum values of tensile and compressive principal stress are 3.83 MPa and 7.66 MPa, respectively, which is bigger than those considering the initial tensile strength and initial shear strength. Therefore, in order to realistically simulate the seismic responses of the arch dam, the initial tensile strength and shear strength of contraction joints should be considered in the numerical analysis. It is recommended to carry out initial tensile strength and shear strength experiments of contraction joints of arch dams so as to obtain more accurate initial tensile strength and shear strength values.