Study on Influence of Numerical Simulation Accuracy on High Core Wall Rockfill Dam Deformation and Crack Analysis

Abstract

With the increasing height of rockfill dams, detailed and effective numerical simulation for the dam design and construction has become increasingly important. This paper develops a detailed and effective finite element simulation program for the construction and impoundment process of high large-scale earth rockfill dams and studies the influence of numerical simulation accuracy on dam deformation and crack analysis. By establishing the subdivision mode of conventional elements, the mesh generation of detailed finite element simulation is realized. Using the row compression storage method and the MPI parallel method, the calculation program of the total stiffness matrix is developed on the Linux operating system. On the basis of the detailed simulation program above, the filling and impoundment process of Pubugou high core wall rockfill dam is simulated, and it is found that the influence of simulation accuracy on the core wall zone is greater than that on the rockfill zone. Compared with the dam construction process, the detailed simulation of the impoundment process is more necessary. When simulating the impoundment process, there is a significant effect of simulation accuracy on the simulation result of the upstream side of the core wall and the upstream filtration zone. So, it is necessary to simulate the impoundment process of the earth rockfill dam in detail. Moreover, using the detailed method to simulate the impoundment process of earth rockfill dams can more accurately predict dam cracks.

1. Introduction

Rockfill dams are one of most preferred dam types in hydraulic engineering due to their use of local materials, good adaptation to poor geological conditions, simple construction sequence, short construction period and good seismic performance [1]. They are currently constructed with great frequency all over the world, especially in China. With the development of construction technology, the height of rockfill dams has significantly increased. Before the 21st century, the height of rockfill dams constructed in China was below 200 m. In the new century, China built many rockfill dams that are more than 200 m high, with some dams even reaching more than 300 m. Due to the great height and large storage capacity of the dams, the safety of high rockfill dams is a major concern that attracts the attention of engineers and scientific researchers [2,3]. High rockfill dams mainly contain core wall rockfill dams and concrete face rockfill dams. The structural behaviors of these two types of dams, such as stress and deformation, uneven settlement, and arch effects, are directly related to the safe operation of dams [4]. The currently used concept and method of high dam design are basically the same as those that have been used since the middle of the last century. It is difficult to adapt to the rapid development of high dam construction practices; thus far, there is no design specification for high dams above 200 m, which is a prominent issue in the field of dam safety [5,6].

Deformation control is a key problem in high earth rockfill dam design and construction [7]. For the concrete face rockfill dams, long-term time-dependent deformation control is one of the most important research objects [8,9]. It has been reported that some rockfill dams have continuous settlements decades after their construction [10]. Cethana Dam in Australia is still deformed after 10 years of operation [11]. Within one and a half years of completion, the settlement of downstream rockfill area in Aguamilpa Dam in Mexico increased by 250 mm, the deflection of the top panel increased by about 50~70 mm, and the deformation of the dam body is still in progress [12]. As the supporting structure of concrete face slab, large deformation will inevitably affect the face slab. During operation, when the lower part of the face slab loses the support of the rockfill, the face slab will be empty, and it is more likely to be crushed or collapse under the action of the reservoir water pressure [13,14]. Excessive creep deformation is one of the main reasons for the voids of Tianshengqiao concrete face rockfill dam in China [15]. For the core wall rockfill dams, the wetting deformation of upstream rockfill is considered to be the main cause of the collapse settlement and uncoordinated impoundment deformation of core wall rockfill dams, so the wetting deformation control is one of the key problems [16,17,18]. During the impoundment, wetting deformation causes the collapse settlement in the upstream shell of core wall rockfill dams, which may increase the settlement of the rockfill relative to the core; the discordant deformation of a core wall rockfill dam tends to cause cracks in the dam body and aggravates the core wall arching effect, which can cause the dam body to fail [19,20,21]. A significant number of engineering cases of differential settlement, crack and even the failure of rockfill dams and highway embankments caused by wetting-induced settlement have been reported. As is shown in

Table 1. The case of earth-rock dams showing cracks or even collapsing at first impoundment.

Dam Name	Height/m	Country	Cases
El Isiro Dam	30.0	Venezuela	It showed an approximately 90-m long and 40–60 cm wide crack on the downstream slope during impoundment [22].
Guanyinyan Dam	79.0	China	Six cracks appeared on the upstream side at the junction and four cracks appeared at the dam crest along the dam axis during impoundment [3].
Teton Dam	126.5	America	Wetting-induced collapse settlement has been cited as the key component collapse settlement of the Teton dam [23].
Infiernillo Dam	148.0	Mexico	It developed a sudden settlement in the upstream shell and crest. A crack formed along the dam axis due to settlement during impoundment [22].
Xiaolangdi Dam	154.0	China	A longitudinal crack appeared at the crest during first impoundment [18].
Ataturk Dam	169.0	Turkey	It developed serious settlement, and even longitudinal cracks at the dam crest and an upstream rockfill slide at the first impoundment [24,25].
Pubugou Dam	186.0	China	A crack 230 m long, 5 cm wide, and 1.0–2.5 m deep appeared at the crest during the first impoundment, 5.5–6.0 m from the axis [26,27].

, many earth-rock dams show cracks or even collapse at the first impoundment. At the same time, in the case of the dam break in Sanpan region, Laos, which was reported by media in July 2018, the main cause of the initial cracks in the dam was the wetting deformation [18].

Moreover, many cases have been reported in the literature which involve damage or distress related to a collapse in fills other than dams [28]. Thus, it is of great engineering significance to study the deformation control and safety evaluation of high rockfill dams.

Numerical analysis plays a crucial role in modern geotechnical engineering. After approximately five decades of development, numerical analysis in geotechnical engineering has reached a certain degree of maturity [29]. Various numerical methods are adopted in geotechnical modeling. Among them, grid-based methods such as the finite element method (FEM), finite difference method (FDM) and finite volume method (FVM) are well-established and successfully employed in various applications [30,31]. Therefore, numerical simulation is an efficient approach to deformation control analysis and the safety evaluation of high rockfill dams. Currently, the finite element method (FEM) is the primary method used for the numerical simulation of dams. The statistical data show that compared with the field monitoring data, the current phenomenon in the calculation of rockfill dams is that “the calculated deformation of low dam (dam height less than 100 m) is large, the calculated deformation of high dam (dam height more than 100 m) is small” [32]. Thus, it is very important to improve the numerical simulation accuracy of rockfill dams for dam construction and safety evaluation [33,34].

The simulation accuracy is influenced by several factors, such as the constitutive model of the soil and rockfill materials, the similarity between the working conditions assumed in the numerical simulation and those applied in practice, the calculation scale of the simulation, and the parameters used in the calculations [27,35].

From the perspective of the constitutive model, several models have been used to calculate the deformation of the earth-rock dam, including the Duncan–Chang model, the Cambridge model, the double yield surfaces model, the generalized plastic model and the KG model. Among them, the Duncan–Chang model is widely used in engineering and research because of its simplicity, clear parameter meaning, accuracy and relevance to the working conditions, etc. [4,36]. In terms of determining the parameters of the constitutive model, the parameters used in calculation are generally obtained by arranging the data of three triaxial tests of dam materials. However, the specimens in the indoor large axial test are generally cylindrical, and the maximum particle size allowed is generally less than 5 cm, which is much smaller than the maximum particle size of the on-site dam material that is about 100 cm. Although the scale transformation of the sample gradation is carried out in the experiment, the parameters of the constitutive model obtained still have great problems [37]. Many studies have attempted to use the back-analysis method to accurately obtain the parameters based on in situ monitoring data, e.g., the improved genetic algorithm [4], the intelligent evolutionary algorithm [38] and parallel mutation particle swarm optimization [39].

From the perspective of working condition simulation, the FEM element size is a crucial factor, the smaller the element size discretizes, the more detailed the simulation of construction and the impoundment process and the higher the accuracy that can be obtained. Due to huge bulk volume and extremely complex geometrical and geological features, it is forbiddingly difficult to perform a dam structural analysis with even moderate geometry fidelity in engineering practices [40,41]. Rockfill dams generally employ layered rolling construction: the thickness of each layer is 0.6~1.0 m, so the size of the element should not be greater than 1.0 m. However, in terms of the current calculation scale of rockfill dams, the calculation element number is generally less than 100 thousand, and the maximum element size up to ten meters, even tens of meters. The reasons for this are as follows: (a) mesh generation: rockfill dams have a complex topographic condition, dam material zoning and construction schemes, which make it impossible to realize automatic mesh generation; (b) calculation software: by the complexity of rockfill dam structure, dam materials’ stress–strain relationship, and construction and impoundment schemes, the method of calculating rockfill dams is usually by using software developed by researchers themselves, or using general FEM software after a rather lengthy secondary development, and there is a significant challenge for implementing large-scale computing [42].

In order to simulate the construction and impoundment process of rockfill dams in detail, the authors have proposed a detailed simulation method for dam simulation and implemented it through a self-developed program. Moreover, the detailed simulation method is used to study the influence of the numerical simulation accuracy on the simulation results of high core wall rockfill dam in this paper. The method mainly contains a mesh refining method for initial FEM elements and some improvement on the original FEM calculation method for rockfill dams. In the first part of this paper, the detailed simulation method including the mesh refining method and numerical simulation improvement is introduced. Then, two simulation models representing different simulation precision of the Pubugou high core rockfill dam are established, and the two models are used to simulate the dam construction and impoundment process with different precision. On the basis of studying the simulation results of different precision, the paper finally points out the influence of numerical simulation precision on the simulation results of dam construction and the impoundment process.

2. Detailed Simulation Method

The authors used the mesh generation method to establish the detailed model of earth-rock dam by determining the subdivision rules of common elements and programming methods. The total stiffness calculation and equation solution are carried out by means of parallel computing, and the detailed numerical simulation of earth rockfill dam with more than one million elements is realized. This section introduces the mesh generation process, the parallel computing of total stiffness matrix and the simulation process of earth-rock dam construction and the impoundment process.

2.1. Mesh Generation

The model of a 3-D earth rockfill dam is mainly composed of cube elements, as shown in Figure 1a. However, in the dam slope, material partition boundary and the connection between the dam body and valley, other types of elements are often added to meet the need for structural transitions. The surface of solid elements in the finite element method numerical simulation of the earth rockfill dam is usually quadrilateral or triangle. According to the surface-to-surface, surface-to-line and surface-to-point transition forms, we can know all the element types of the 3-D earth rockfill dam model: as shown in Figure 1a, the cube element formed by the transition from quadrilateral to quadrilateral surface; as shown in Figure 1b, the heptahedron element generated by a quadrilateral-to-triangular transition; as shown in Figure 1c, the triangular prism element generated by the quadrilateral-to-straight line transition; as shown in Figure 1d, the pyramid element generated by the concentration of quadrilaterals-to-a point; as shown in Figure 1e, the triangular prism element generated by triangular face to line transition (the same as Figure 1c); as shown in Figure 1f, the pyramid element generated by triangle-to-straight line transition (the same as Figure 1d); and as shown in Figure 1g, the tetrahedral element generated by the triangle concentrated to a point. When the element is refined, the segmentation method of the contact surface or line between the elements is fixed, and the refined model will not appear at the suspension point.

There are two types of contact surfaces: quadrilateral and triangle. According to the idea of dividing each line into two segments, the quadrilateral surface adds a face center node besides the midpoint of the line, and no other nodes are added in the triangular surface. As shown in Figure 2a, each quadrilateral is divided into four quadrilaterals, and as shown in Figure 2b, each triangular surface is divided into four triangles.

In addition, when subdividing the solid element in Figure 1, the corresponding internal nodes of the element need to be added to realize the secondary subdivision of the mesh. In cube subdivision, a body center node must be added to connect to the center node of the six faces. Then, each cube is divided into eight cubes, as shown in Figure 3a. When the heptahedron is subdivided, a node similar to the center of the body is added to connect the center of the three quadrilateral surfaces and the midpoint of each side of the triangular surface that is not collinear with the quadrilateral surface. The element can be divided into eight parts, including four cubes, three heptahedrons and a tetrahedron. However, one surface of the three cubes needs to be segmented by diagonal lines, so the three cubes are divided into two prisms, and the remaining tetrahedron is divided into a heptahedron and a tetrahedron. As shown in Figure 3b, a heptahedron is divided into four heptahedrons, six prisms and two tetrahedrons. When the triangular prism is subdivided, there is no need to add internal nodes, directly according to the surface segmentation method, such as in Figure 3c, where a triangular prism is subdivided into eight triangular prisms. Pyramid elements also need not increase internal nodes, as shown in Figure 3d; after subdivision, they are divided into six pyramid elements and four tetrahedral elements. Tetrahedron elements also do not need to increase internal nodes. After cutting four vertices, an octahedron is formed. The octahedron can be divided into four tetrahedrons by properly connecting the two vertices corresponding to the octahedron, as shown in Figure 3e, where each tetrahedron is divided into eight tetrahedrons.

2.2. Parallel Computation of Total Stiffness Matrix

In the finite element simulation of the earth-rock dam, the total stiffness matrix is a typical symmetric sparse matrix, whose storage methods include the diagonal storage method, variable bandwidth storage method, coordinate storage method and its improvement method, Ellpack–Itpack storage method, row compression storage method (CSR), Sherman’s storage method, super matrix storage method, etc. [43]. However, the computing time and storage space restrict the large-scale numerical simulation in the engineering field. The design of the large-scale stiffness matrix storage method and its parallel solution technology can penetrate this bottleneck. Only storing non-zero elements in the overall stiffness matrix is the most storage-saving space method, and the corresponding solver should be designed according to the corresponding storage format [44].

In this paper, the row compression storage method is used, which requires three one-dimensional arrays $\mathbf{A}$, $\mathbf{I}\mathbf{A}$ and $\mathbf{J}\mathbf{A}$ to determine the position and value of non-zero elements in the sparse matrix. As shown in Figure 4, $\mathbf{A}$ stores the values of non-zero elements; $\mathbf{J}\mathbf{A}$ stores the column number of elements in $\mathbf{A}$; and $\mathbf{I}\mathbf{A}$ stores the serial number of the first non-zero element in each row in $\mathbf{A}$. For the symmetric matrix, only the upper triangular part is stored, and the first non-zero element in each row is its diagonal element.

In the process of earth-rock dam construction simulation, with the increase in construction elevation, the stiffness matrix of the element is constantly changing. Each construction load step needs to calculate the stiffness matrix of all the current completed construction elements and use it to integrate the total stiffness matrix. If each load step calculates three arrays of $\mathbf{A}$, $\mathbf{I}\mathbf{A}$ and $\mathbf{J}\mathbf{A}$ above, the time consumed in determining the total stiffness matrix and solving the equations is basically in the same order of magnitude [42]. It is lucky to discover that in the calculation of the total stiffness matrix, the stiffness indicator matrix of $\mathbf{I}\mathbf{A}$ and $\mathbf{J}\mathbf{A}$ does not change with the change of the element stiffness matrix, only related to the dam mesh model and constraints. So, the matrix of $\mathbf{I}\mathbf{A}$ and $\mathbf{J}\mathbf{A}$ can be determined before the stress–strain calculation. $\mathbf{I}\mathbf{A}$ and $\mathbf{J}\mathbf{A}$ are determined to find the position of non-zero elements in the total stiffness matrix. The method for determining a row (column) non-zero element in the total stiffness matrix is: (1) determine the nodes corresponding to this row (column); (2) find the nodes in the same structure element as this node, record and sort; (3) the non-zero elements in the row (column) of the total stiffness matrix are other non-constraint directions of this node itself and all non-constraint directions corresponding to other nodes in the common structure element of this node.

The array $\mathbf{I}\mathbf{A}$ and $\mathbf{J}\mathbf{A}$ can be determined by the above method. However, a large number of elements need to be looped when a node looks for the nodes that are associated with it. In large-scale numerical simulation, such a nested loop is very time-consuming, and in this search process, no matter the node cycle or the element cycle, there is no inevitable dependence, so we can use parallel means to improve its efficiency.

Based on the above characteristics, this paper uses the MPI (Message Passing Interface) parallel method on the Linux system to determine the non-zero element indicator matrix IA and JA in the total stiffness matrix. The specific process is shown in Figure 5. Moreover, the solution of the equations is realized by calling the parallel solver after the calculation of the total stiffness matrix and load vector has been completed.

2.3. Simulation of Construction and Impoundment Process

The numerical simulation of earth-rock dam needs to be combined with the actual situation to simulate the construction and impoundment process. During the construction period, the dam produces corresponding deformation under its own gravity. During the impoundment period, the dam deforms under the action of the water load, buoyancy and wetting deformation of the upstream rockfill materials. In addition, the construction and water storage period of high earth-rock dam is long, and can reach several years or even more than a decade. Therefore, it is necessary to consider the creep deformation during the construction and impoundment process. Figure 6 shows the numerical simulation process of an earth-rock dam’s construction and impoundment process.

3. Detailed Simulation of Pubugou Dam

In order to verify the feasibility of the proposed method above, it is necessary to carry out the detailed numerical simulation of the actual high core wall rockfill dam project. In this paper, the influence of numerical simulation accuracy on high core wall rockfill dam deformation and crack analysis is studied. Therefore, the high core rockfill dam with cracks in the process of impoundment—Pubugou Dam is selected.

3.1. Project Description

The Pubugou Hydropower Station is located on the middle reaches of the Dadu River. It is an extra-large hydro-junction project with power generation as the main purpose and the comprehensive benefits of flood control and sediment interception. The normal reservoir water level is at an elevation of 850.0 m, and the dead water level at an elevation of 790.0 m. The total storage capacity is 5.39 billion m³. The total installed capacity of the power station is 3600 MW, and the average annual power generation is 14.79 billion degrees.

The dam is a gravel soil core wall rockfill dam with a crest elevation of 856.0 m, a maximum dam height of 186.0 m, a crest width of 14.0 m and a crest length of about 540.5 m. The top elevation of the core wall is 854.0 m, and the width is 4.0 m. The slope of both the upstream and downstream dam slopes is 1:0.25, and the bottom elevation is 670.0 m. The section consists of five zones, namely gravel soil core wall, filter layer, transition layer, main rockfill area and secondary rockfill area within the downstream dam body. There are two layers of filter layers on both sides of the upstream and downstream reaches of the core wall, with the upstream width of 4.0 m and the downstream width of 6.0 m. There is a transition layer between the rockfill of dam shell and the filter layer, and the slope of the contact surface between the transition layer and the rockfill of the dam shell is 1:0.4. The typical section of the dam is shown in Figure 7 [45].

The construction of the power station began on 30 March 2004, the closure began in late November 2005, the construction of the upstream and downstream rockfill began in March 2006, the absolute section construction began in April 2007 and the dam was filled to the top on 20 September 2009. From 1 November 2009 to 13 December 2009, the first impoundment stored the water level to 790.0 m of elevation. From 8 May 2010 to 13 October 2010, the second impoundment increase the water level to 850.0 m of elevation. However, when the reservoir was impounded to about 842.2 m of elevation on the morning of 26 August 2010, the longitudinal cracks appeared at the dam crest. The cracks are located on the downstream side of the dam crest. The distance from the dam axis is about 5.5 to 6.0 m, and almost parallel to the dam axis. The crack length is about 230 m, and the maximum crack width is about 5 cm. Through the pit inspection, it is found that the crack depth is about 1.0~2.5 m [46,47].

3.2. Model and Simulation

3.2.1. Finite Element Models

According to the actual situation of cracks during the impoundment of earth-rock dams, tensile stress may occur in the dam crest, the upstream dam slope and the downstream dam slope near the dam crest—that is, cracks may occur in these zones. Therefore, the denser the mesh in these parts, the higher the accuracy of the calculation results. In this part, the initial mesh generation of Pubugou dam is carried out based on the traditional mesh generation technology. As shown in Figure 8a, the finite element model has 15,967 elements and 16,682 nodes. Then, by using the mesh subdivision method in Section 2.1, the initial mesh is subdivided twice, and the detailed finite element model is obtained. As shown in Figure 8b, there are 1,032,907 elements and 981,194 nodes.

3.2.2. Construction and Impoundment Process Simulation

The actual construction and impoundment process of the dam and its simulation are shown in Figure 9. The solid line is the actual construction and impoundment process, and the dashed line is the simulation process. In the process of dam impoundment, the water pressure is directly applied to the upstream side of the core wall, and the seepage in the dam is not considered. However, the buoyancy and wetting deformation of the upstream rockfill materials in the water level change area and the creep deformation during the construction and impoundment process are considered.

In order to accurately simulate the dam construction and operation, each construction and impoundment load step simulates one layer of the mesh at most. The initial finite element model used 35 construction load steps and 25 impoundment load steps to simulate the dam construction and impoundment process. Considering the creep deformation of dam materials in the process of construction and impoundment, 25 creep load steps were inserted into the construction load step according to the actual construction progress. Combined with the process of water level change, 15 creep load steps were inserted into the impoundment load step. After the impoundment, 11 creep load steps were added to simulate the date to 31 December 2012.

However, when using the detailed finite element model simulation, each construction and impoundment load step is divided into four load steps to ensure that each construction and impoundment load step simulates a layer of the mesh at most. Thus, there are 140 construction load steps and 100 impoundment load steps. To ensure the comparability of the simulation results, the insertion time of the creep load step is the same as that of the initial model, with 51 creep load steps.

Duncan and Chang’s EB constitutive model is used in this paper to simulate the instantaneous deformation under dead load. There are seven parameters, i.e., $c$, $\varphi$ (or ${\varphi }_0$, $\Delta \varphi$), $R_f$, $K$, $n$, $K_b$ and $m$, which can be evaluated by using a group of conventional triaxial tests. The model parameters are shown in

Table 2. The EB model parameters of the main materials in Pubugou Dam.

Dam Material	$\mathit{\gamma }$	$\mathit{K}$	$\mathit{n}$	${\mathit{R}}_{\mathit{f}}$	${\mathit{K}}_{\mathit{b}}$	$\mathit{m}$	$\mathit{\varphi }(\circ )$	$\Delta \mathit{\varphi }(\circ )$	$\mathbf{c}$
Core wall	2.36	550	0.42	0.76	240	0.29	35.0	0.0	0.12
Filter	2.03	790	0.59	0.81	400	0.30	35.5	0.0	14.1
Transition	2.15	986	0.36	0.74	550	0.32	38.8	0.0	11.5
Major rockfill	2.10	1000	0.52	0.68	420	0.34	54.0	10.0	0.0
Minor rockfill	2.10	800	0.50	0.70	318	0.30	51.0	10.0	0.0

[48].

The seven parameters creep model is used to simulate the creep deformation of the dam, and the Merchant equation is used to describe the creep curve in the model,

$$\mathsf{\epsilon }\left(\mathrm{t}\right)={\mathsf{\epsilon }}_{\mathrm{i}}+{\mathsf{\epsilon }}_{\mathrm{f}}\left(1-{\mathrm{e}}^{-\mathsf{\omega }\mathrm{t}}\right)$$

(1)

where $\mathsf{\epsilon }\left(\mathrm{t}\right)$ is the creep strain developed at time t; ${\mathsf{\epsilon }}_{\mathrm{i}}$ and ${\mathsf{\epsilon }}_{\mathrm{f}}$ are the initial and permanent creep strains, respectively; $\mathrm{e}$ is the natural index; and $\omega$ is a parameter representing the initial relative deformation rate (or creep strain during the first day). Fang [49] gave the improved calculation formulas of permanent volumetric creep strain ${\mathsf{\epsilon }}_{\mathrm{vf}}$ and permanent shear creep strain ${\gamma }_{\mathrm{f}}$, as shown in Equation (2),

$$\{\begin{array}{c}{\mathsf{\epsilon }}_{\mathrm{vf}}={\mathrm{b}}_1{\left(\frac{{\mathsf{\sigma }}_3}{{\mathrm{P}}_{\mathrm{a}}}\right)}^{{\mathrm{m}}_1}+{\mathrm{c}}_1{\left(\frac{\mathrm{q}}{{\mathrm{P}}_{\mathrm{a}}}\right)}^{{\mathrm{m}}_2}\\ {\gamma }_{\mathrm{f}}={\mathrm{d}}_1{\left(\frac{\mathrm{S}}{1-\mathrm{S}}\right)}^{{\mathrm{m}}_3}\end{array}$$

(2)

where ${\mathrm{b}}_1$, ${\mathrm{c}}_1$, ${\mathrm{d}}_1$, ${\mathrm{m}}_1$, ${\mathrm{m}}_2$ and ${\mathrm{m}}_3$ are parameters.

Table 3. Pubugou Dam materials’ creep model parameters.

Dam Material	${\mathrm{b}}_1(\%)$	${\mathrm{c}}_1(\%)$	${\mathrm{d}}_1(\%)$	${\mathit{m}}_1$	${\mathit{m}}_2$	${\mathit{m}}_3$	$\mathit{\alpha }$
Core wall	0.1075	0.385	0.717	0.936	0.679	0.518	0.00300
Major rockfill	0.0575	0.140	0.454	0.383	0.365	0.482	0.00649
Minor rockfill	0.0975	0.160	0.612	0.797	0.455	0.542	0.00600

shows the creep model parameters which were obtained by the Nanjing Institute of Water Resources and the Yangtze River Academy of Sciences, according to the indoor creep test [48].

The wetting deformation model using Shen’s improved model [50] is denoted in Equation (3).

$$\{\begin{array}{c}\Delta {\epsilon }_v^w=C_w{\left(\frac{{\sigma }_3}{P_a}\right)}^{n_w}\\ \Delta {\gamma }^w=D_w\frac{S_L}{1-S_L}\end{array}$$

(3)

Table 4. Pubugou Dam materials wetting deformation model parameters.

Parameters	${\mathrm{C}}_{\mathrm{w}}(\%)$		${\mathrm{D}}_{\mathrm{w}}(\%)$
Value	0.0547	1.367	0.265

shows the wetting model parameters with reference to the simulation of Pubugou Dam wetting deformation by Guo et al. [27].

3.3. Simulation Results

3.3.1. Completed Construction Period

Figure 10 shows the contour map of the displacement, stress and stress level of the dam’s 240.0 m section for the completed construction. The red line is the simulation result of the initial model, and the green line is the simulation result of the detailed model. It can be seen from the figure that the distribution and magnitude of displacement, stress and stress level in the rockfill zone by the two models are almost the same. There are differences in the stress and stress level in the core wall zone, especially in the lower part of the core wall.

Therefore, the accuracy of the construction process simulation of the core rockfill dam has more influence on the simulation results of the core wall zone than the rockfill zone. The displacement, stress and stress level of the rockfill zone and the upper part of the core wall zone are almost not affected by the simulation accuracy.

3.3.2. Impoundment Period

According to the in situ data, longitudinal cracks appear at the crest of the dam when the water level reaches above 842.2 m of elevation. Therefore, this section analyzes the displacement, stress and stress level when the water level is at 844.0 m of elevation. Figure 11 shows the contour map of the displacement, stress and stress level of the dam’s 240.0 m section when the water level reaches 844.0 m of elevation. The red line is the simulation result of the initial model, and the green line is the simulation result of the detailed model. It can be seen from the figures that the horizontal displacements simulated by the two methods are quite different in the upper middle rockfill zone, the upstream side of the core wall and the filter zone, and there is no obvious difference in the other zones. There is no significant difference in settlement. The stress distribution is also quite different in the upstream side of the core wall and the filter zone, and there is no significant difference in the other regions. Using fine simulation, the tensile stress zone appears in the downstream side of the dam crest and dam slope. The stress levels obtained by the two models are quite different in the whole section.

In order to explore the influence of detailed simulation on the impoundment process, the displacement increments caused by impoundment are analysed. Figure 12 shows the contour map of displacement increment of the dam’s 240.0 m section caused by impoundment when the water level is at 844.0 m of elevation. The red line is the simulation result of the initial model, and the green line is the simulation result of the detailed model. It can be seen from the figures that the horizontal displacement increment caused by the impoundment simulated by the two methods is quite different in the middle of the rockfill area in the upstream of the dam and the middle of the core wall. The horizontal displacement and settlement increment caused by impoundment simulated by the detailed model are larger than that simulated by the initial model.

3.3.3. Comparison of Simulation Results between the Two Models

Table 5. The eigenvalue and rates change of the Pubugou dam simulated displacement.

Period	Completed Construction Displacement (cm)			Impoundment Increment Displacement (cm)
Direction	Upstream	Downstream	Settlement	Upstream	Downstream	Settlement
Initial model	22.65	39.85	278.62	5.18	55.20	62.68
Detailed model	22.73	40.10	281.04	8.10	59.38	63.32
Change rate (%)	0.35	0.63	0.87	56.37	7.57	1.02

shows the eigenvalue and its change rate of the completed construction displacement and impoundment displacement increment of Pubugou Dam, which was obtained by the two models’ simulation above. It can be seen that all the eigenvalues of the detailed model simulation are greater than those of the initial model simulation. Compared with the results of the initial model simulation, the change rates of the detailed model simulation results during completed construction are all within 1 %, and the simulation detailed degree has little effect on the displacement simulation results in the dam construction process. However, the impoundment displacement increment of the detailed model simulation is obviously larger than that of the initial model simulation. All the change rates of the horizontal displacement are larger than 5%; the change rate of horizontal displacement in the upstream direction is even 56.37%. Therefore, it is necessary to simulate the impoundment process of the earth rockfill dam in detail.

To summarize the above, the simulation detailed degree of the core wall rockfill dam impoundment process has an impact on the displacement and stress simulation results of the core wall area and the rockfill area, especially on the upstream side of the core wall and the upstream filter areas. Compared with the construction process, the fine simulation of the water storage process is more necessary. A detailed simulation of the impoundment process of the core wall rockfill dam can more effectively predict the impoundment deformation of the dam, and the insufficient detailed degree of the simulation will underestimate the impact of impoundment deformation on the dam’s safety. The influence of the simulation detailed degree should be considered in the study of the arching effect and hydraulic fracturing.

3.3.4. Dam Cracks Analyzation

The Pubugou core wall rockfill dam shows cracks at first impoundment, so it is necessary to judge the cracks based on the numerical simulation results. The deformation inclination method is widely used in the analysis and judgment of cracks in earth-rock dams. The method is defined as the ratio of settlement difference between two points in a horizontal plane to horizontal distance, i.e.,

$${\gamma }_S=\frac{\Delta S}{\Delta L}$$

(4)

where ${\gamma }_S$ is the settlement inclination, $\Delta S$ is the settlement difference of two points at the same elevation, and $\Delta L$ is the horizontal distance between these two points. The greater the inclination, the more likely it is to produce cracks. Previous engineering experience shows that in the earth dam, if the inclination is less than 1.0%, no cracks will occur.

Figure 13 shows the two models’ simulated deformation inclination of the dam surface corresponding to the water level at the elevation of 844.0 m. It can be seen from the diagram that the deformation inclination value reaches the maximum value on the dam crest. The maximum deformation inclination of initial model simulation is about 0.7%, which is less than the critical value of 1.0%. If the deformation inclination method is used to analyze the crack at the dam crest, it is judged that there is no crack, which is inconsistent with the actual situation. However, the maximum deformation inclination at the dam crest obtained by the detailed model simulation reaches the critical value of 1.0 %, so cracking may occur, which is consistent with reality. Therefore, the use of the detailed method to simulate the impoundment process of earth rockfill dams can more accurately predict dam cracks.

4. Conclusions

This paper focuses on solving the problems of mesh generation and calculation in the large-scale finite element simulation of earth rockfill dams and develops a detailed finite element simulation method for the construction and impoundment process of high earth rockfill dams. The programming language is used to develop the corresponding modules in the numerical simulation, such as mesh generation, stiffness matrix calculation and equation solving.

For the mesh generation, the subdivision mode of conventional elements is established by specifying the subdivision rules of connection parts between elements, and the mesh generation of detailed finite element simulation is realized. For the stiffness matrix calculation, the row compression storage method is used to store the total stiffness matrix, and the MPI parallel method is used to develop the calculation program of the indicator matrix in the total stiffness matrix on the Linux operating system. To research the influence of the simulation detailed degree on the high rockfill dam numerical simulation results, the construction and impoundment process of Pubugou high core wall rockfill dam was simulated using the initial and detailed simulation method. The results show that when simulating the dam’s construction process, the influence of simulation detailed degree on the core wall zone is greater than that on the rockfill zone, and when simulating the impoundment process, there is a significant effect of simulation detailed degree on the upstream side of the core wall and the upstream filtration zone. The main reason is that the nonlinearity of the core wall materials is stronger than that of the rockfill area, and the stress condition is more complex. Compared with the construction process, the fine simulation of the water storage process is more necessary. The simulation detailed degree has little effect on the displacement simulation results in the dam construction process, but it is necessary to simulate the impoundment process of the earth rockfill dam in detail. The influence of the simulation detailed degree should be considered in the study of the arching effect and hydraulic fracturing. Moreover, using the detailed method to simulate the impoundment process of earth rockfill dams can more accurately predict dam cracks.

Due to the strong nonlinearity of rockfill dam material and the complex simulation process, there are few studies on the influence of the simulation accuracy on the deformation and crack analysis of the dam. This paper accounts for the deficiency in this aspect. In this article, it is found that the core wall area rather than the rockfill area is more sensitive to the simulation accuracy, and the water storage process is more sensitive than the simulation’s construction process detailed degree. Therefore, it is suggested that the scholars should attempt to refine the simulation of the core wall as much as possible in the numerical simulation of core wall rockfill dams. The dam construction process simulation can be simple, but the water storage process must be more detailed, and detailed simulation can be used to analyze dam cracks more accurately.