Three-Dimensional Numerical Analysis of Compound Lining in Complex Underground Surge-Shaft Structure

The mechanical behavior of lining structure of deep-embedded cylinder surge shaft with multifork tunnel is analyzed using threedimensional nonlinear FEM. With the elastic-plastic constitutive relations of rock mass imported and the implicit bolt element and distributed concrete cracking model adopted, a computing method of complex surge shaft is presented for the simulation of underground excavations and concrete lining cracks. In order to reflect the interaction and initial gap between rock mass and concrete lining, a three-dimensional nonlinear interface element is adopted, which can take into account both the normal and tangential characteristics. By an actual engineering computation, the distortion characteristics and stress distribution rules of the dimensional multifork surge-shaft lining structure under different behavior are revealed. The results verify the rationality and feasibility of this computation model and method and provide a new idea and reference for the complex surge-shaft design and construction.


Introduction
The surge shaft of hydropower station is a typical building in underground engineering.In recent years, with the construction of many large-scale hydropower stations, the inner diameter of surge shafts has become larger and larger, and the constructive type and geological condition have become more and more complex.According to NATM, the compound lining structure is adopted in many large-scale surge shafts, with anchoring and shotcreting support applied outside and integral lining built inside.As an important component of hydropower station structure, how to properly analyze the stability of surrounding rocks in the construction period and the structural security in the operation period is a critical technique problem.
With regard to the structural analysis theories of surge shaft, many scholars did a lot of research work [1][2][3][4].The two most popular methods at present are traditional structural mechanics method and FEM.As to the former, the effect of surrounding rocks on lining was reflected by considering the pressure and elastic resistance of rocks, without considering the deformation compatibility and internal force adjustment in the lining-rock interaction, and the results cannot properly reflect the real stress state of structure.The latter is capable of considering the characteristics of surrounding rocks, external loads, and edge-restraint condition, and therefore it has highly forecasting precision and is applicable to the surgeshaft characteristics of multivariate bodily form and complex strained condition.Papers [5][6][7][8][9][10] use FEM to analyze the stability of surrounding rocks and operating characteristics of lining in simple formal underground structure but did not involve the analysis of complex multifork surge-shaft structure.
Based on the elastic-plastic constitutive relations of rock mass, this paper simulates the progress of excavation and lining buildup in deeply buried circular restricted orifice surge shaft by FEM, simulates the interaction and initial gap between lining and surrounding rocks with the 3D eight-node nonlinear interface element, and further proposes a calculation method for the lining cracks and structure reinforcement.In this calculation model, the main influencing factors (such as the property of surrounding rocks, the excavating and supporting process of cavern, and the interaction and initial gap between surrounding rocks and lining) of lining stress characteristics could be effectively reflected.Finally, by an actual engineering computation, the structural stability, the distortion, cracking characteristics, and the stress distribution of complex multifork surge shaft lining under different condition are revealed, which verifies the rationality and feasibility of this computing method and provides some helpful ideas and references for the structure design and construction.

Numerical Analysis Method of Surge-Shaft
Lining Structure 2.1.Simulation of Excavation, Support, and Lining Buildup.
After the underground chambers are excavated, the stress variation in surrounding rocks is easy to cause cracks and damage.This paper describes the stress variation state of rock mass after excavation with 3D elastic-plastic damage constitutive relations [11]: where {  }  is the stress of rock mass after cracking and damage, {  } is the stress corresponding to elastic-plastic constitutive relation, {  } are three normal stresses, and  is the damage coefficient calculated as follows: where   , , and  are curve fitting constants in material damage experiment and   is the 2nd invariant of rock plastic strain tensor.
The excavation load of rock mass could be divided into two parts: elastic load {  } and plastic load {  }: where  is the structural elastic coefficient after rock excavation and the value is set to elastic coefficient of the element which firstly enters in plastic phase; that is, where   is elastic coefficients of all elements that enter in plastic state after excavation.For the elastic coefficient of a certain element, theyield function  of rock mass is calculated as follows: where { 0 } is the initial stress of rock mass before excavation and {Δ} is the stress incremental of rock mass after excavation.
The elastic load {  } is calculated once, and the plastic load {  } is calculated by steps; at every step the increment load is iteratively calculated as follows: where [  ] is elastic rigidity matrix, [  ] is damage rigidity matrix, and the relationship of rock damage stress matrix, elastic stress matrix, and plastic stress matrix is as follows: To simulate the anchor bolt and anchor cable support, the implicit anchor bolt and anchor cable element are adopted (for details see [12][13][14]).
The concrete lining structure is built after excavation and support of surge-shaft chamber, considering the influence of chamber deformation on lining structure; the plastic load {  } released from surge-shaft excavation could be divided into two parts: where  is the excavation release coefficient decided by release property of rock mass after excavation.Before the surge-shaft lining buildup (the lining elements are still air elements), the whole structure is calculated iteratively only with the { 1 }.After lining buildup (here lining elements are changed into concrete elements), the global stiffness is reformed, and then { 2 } is imposed on the whole structure for iterative computations.Finally, based on the water level diversification of surge shaft and load combinations in operating period, the deformation and stress characteristics of surge-shaft lining under different conditions are calculated; furthermore, the reinforcement ratio is figured out.

Simulation of Interface between Lining and Surrounding
Rocks.The restraining effect of surrounding rocks on lining, which has great influence on the internal force of lining, could be reflected by setting interface elements [15] between lining and surrounding rocks.Based on the regular thin finite element, this paper introduces an ameliorative 3D eight-node nonlinear interface element (Figure 1) to avoid the drawbacks that only consider the gap area deformation and arbitrarily set the values of normal rigidity in Goodman interface element [16].In constitutive relations only, with nonlinear characteristics in normal and tangential direction considered and nonlinear elastic constitutive relation imported, the interface deformation patterns, such as felt, slippage and puff, are taken into account, which can make the calculation results more actual.Supposing the strain of interface element distributes evenly along thickness, the increment form of constitutive relation could be expressed as where   is the normal of interface and [  ] is the elastic stress matrix of interface element expressed as where ,  is modulus of filling material or contact material between surrounding rocks and lining, and  is Poisson ratio.
(1) To describe the normal stress deformation relation of interface, the hyperbolic model of rock joint normal deformation [17,18] is introduced: where      0 is the initial normal rigidity of interface, its value could be set to  1 ,  is normal relative deformation, and  is thickness of interface.As  is very small,      = /.In this paper, negative value is given to the compressive stress and closing displacement of interface, with positive value given to the tensile stress and opening displacement.
(2) To describe the tangential stress-strains relation, the - hyperbolic model [19] proposed by Clough and Duncan is adopted: where  is the shear stress,  is the tangential displacement,      and      are the tangential rigidity,  is cohesion,  is frictional angle,   is the volume weight of water,   is atmospheric pressure,  and  are two parameters, and   is the destruction ratio (  < 1).
Due to the basic principle of FEM, by the transformation relations between local coordinate system and integral coordinate system, the interface element stiffness matrix in integral coordinate system could be defined as where [  ] is transformation matrix between local coordinate system and integral coordinate system.
The following strength discriminate criterions are adopted.(1) If the normal stress is greater than zero or tensile strength, the interface opens.(2) If the normal stress is less than zero or tensile strength, the Mohr-Coulomb criterion is used to discriminate whether the interface slips.During the iteration, the overloading stress occurring at the interface opening and the shear stress Δ =  − (  + ) which exceeds the shear strength at the interface slipping will be transformed into the nodal loads and distributed to the surrounding elements.If the interface element crazes, the shear modulus and normal rigidity will be diminished.

Simulation of Initial Gap.
Due to the quality and method of construction, the concrete shrinks when the temperature drops, the initial gaps often exist between lining and surrounding rocks, especially at the crown of chamber.Supposing the concrete lining bears all internal water pressure before the gap closes, only if the lining touches surrounding rocks after the gap closes, they both share the internal water pressure.So the initial gap can be simulated by the following methods.
(1) Under internal water operating condition, the interface elements set outside the lining of chamber crown is defined as initial gap elements with thickness of  (Figure 2), and their parameters are given as air elements before the gap closes.In the classification calculation process, at each step of load iteration after the nodal displacement calculation, the following progress will be cycled to all gap elements: first, it calculates the radial relative displacement Δ() of the inner face and outer face, where  is the serial number of elements.If Δ() ≥ , it is regarded that the initial gap is compacted and the lining touches surrounding rocks, so the initial gap elements will be changed into the compacted interface elements, and their stiffness matrix will be calculated by the method hereinbefore.
If () ≤ , no operation is required.Lastly, the loop will be closed after the judgment of all joint elements and enter the next iteration.
(2) Under external water operating condition, supposing that the external water pressure acts on the outer face of lining in the form of surface force, and the crown, sidewall, and base plate of lining are in good agglomeration, so the impact of initial gap is not taken into account.

Calculation Method of Concrete Lining Cracking.
The concrete material is apt to crack under tension that causes a sudden stress variation and stiffness reduction [20].In the iteration process, if the strain of lining elements   >   , where   is the limit tension strain of concrete, the lining element is considered as cracked.In this case, the stress Δ =   −   which exceeds the concrete ultimate tensile strength   will be changed into the nodal loads and applied to the structure to enter the next iterative calculation.After concrete cracking, the elasticity modulus that is normal to the fracture face falls off, and the stress-strains relation changes; at this moment the crack elements are treated as anisotropic material, and the expression  = 1 −   is used to describe its damage breakage properties, where   is the damage variable of concrete material calculated as follows [21,22]: where  is the calculation value of element principal tensile strain,   is the concrete ultimate tensile strain, and   and   are two curve fitting constants; for normal concrete, 0.7 <   < 1, 10 4 <   < 10 5 .Then the stress matrix after concrete cracking is transformed into the matrix in global coordinate system, and the total stiffness matrix is reset to enter the next load iteration: where [𝑅] is the transfer matrix between the three principal directions of elements and the global coordinate.
Based on the element strain value   and structure reinforcement condition calculated by the nonlinear iteration above, the Maximal probability crack width of lining elements can be estimated as follows [23]: where  = 2,  is the protection cover of reinforcement,  is spacing of bars,  =   /  is the elasticity modulus ratio of reinforcement to concrete,  is the reinforcement ratio, and  is the damage coefficient.
If the reinforcement is calculated by the limitative crack width, the concrete after cracking is considered as inactive; namely, the damage coefficient  = 1; then  = 0, and the reinforcement ratio of lining can be calculated as

Computation Model and Material Parameter. Xiaowan
Hydropower Station is a state large project in China, which adopts the underground power house structure, and the installed capacity is 4200 MW (6 × 700 MW).The tailrace surge chamber uses the form of double cylindrical restrictedorifice surge shaft, three gensets share one shaft, the height of surge shaft is 90.0 m, and the diameter is 32.0 m; the lower part uses the form of complicated multifork, the tailwater tunnel is in horseshoe shape, and the diameter is 18 m.The combined form of shotcrete-bolt support and concrete lining is used in the well bore and multifork section.The chamber is deeply embedded in the lightly weathered granite gneiss, and the burial depth of number 1 surge shaft is nearly 500 m.The initial stress field of chamber excavation is simulated by 3D feedback stress field, and the three principal stresses in the center of surge shaft are −37.5 Mpa, −27.53 Mpa, and −11.49Mpa.
Taking number 1 surge shaft and the surrounding rocks as research objects, this paper conducts computational analysis to the load-bearing characteristic of concrete lining and impedance board structure in different operating conditions.The 3D FEM mesh contains 28240 twenty-nodded isoparametric elements, which includes the tailrace surge chamber, bifurcated pipe section, impedance board, and four faults; the whole FEM mesh is shown in Figure 3.The model ranges along -, -, -axes are 131.04,168.18, and 135.00 m, respectively.The thickness of lining at side wall, multifork, and impedance board is about 0.8 m, 1.0 m, and 1.5 m; the FEM mesh of lining structure and impedance board are shown in Figure 4.The C20 concrete is adopted in lining structure, interface elements are set between surrounding rocks and lining, and the physical mechanical parameters of materials are listed in Table 1.

Calculation Condition.
The procedure of construction including excavating, supporting and building lining is simulated.Two load conditions which are considered in FEM calculations are listed as follows.(1) Normal operating condition: The main load is internal water pressure, which is calculated by the top surge water level of 1020.5 m; the acting head on impedance board is calculated by 7.15 m; the structure dead weight is included; from more secure angle, the external water pressure and rock pressure are not considered in calculation.
(2) Overhauling condition: The main load is external water pressure, which is calculated by the underground water level of 1033 m multiplied by reduction coefficient of 0.85; the dead weight is included; in consideration of water ramollescence, partial pressure of rocks is included; no water pressure is considered on the impedance board.

Comparison of Calculation Results under Different Combinations of Materials.
The excavation of surge chamber is simulated by excavating and supporting by stages; then lining is built and surge shaft is filled with water.In normal operating condition, two methods are used for calculation and comparison: (1) materials combination 1, considering the lining as linear elastic material without cracking and the surrounding rocks as elastic-plastic material, and (2) materials combination 2, calculating lining as nonlinear material with cracking and surrounding rocks as elasticplastic material.The calculation results are showed in Figures 5 and 6, with positive value given to tensile axial force and inflexed bending moment.
By comparison of the results under internal water operating condition, the values of lining stress and internal force in material combination 1 are greater than the values in material combination 2. At the cracked place (such as the bottom of surge shaft, the intersection of chamber, the base angle of bifurcation, and top of impedance board) calculated in material combination 2, the principal stresses results of two methods differ between 16% and 27%.The results indicate that in the case of considering initial gap, while lining is regarded as linear elastic material in material combination 1, the stiffness of lining is greater than the stiffness of surrounding rocks.It is difficult to make contact and the lining will bear more load proportion after contacting, so the calculated stress values of lining and internal force are comparatively large.But in material combination 2, lining is calculated as nonlinear material with cracking.Because of the degradation of the section stiffness after lining cracking, the calculated stress and internal force are smaller than the results in material combination 1 and closer to the actual value.Overall, in the case of considering the surrounding rocks and lining as a whole and the nonlinear cracking character of lining, the calculated results are more effective to reflect the combining bearing properties of them and the actual state of stress.

Comparison of Calculation Results under Different Conditions.
In the case of materials combination 2, internal water condition (normal operating condition) and external water condition (overhauling condition) are considered.The calculation results are showed in Figures 7 and 8.

Stress Characteristics of Surge-Shaft Lining and
Impedance Board Structure.It can be seen that, under internal water operating condition, the stress of surge shaft lining augments gradually from top to bottom.It is chiefly because the internal water pressure increases from top to bottom of the surge shaft.The first principal stress is tangential tensile stress, with the maximal value of 2.13 MPa, which occurs at the junction of surge shaft and draft tubes (Figure 9).The third principal stress is radial compressive stress, with the maximal value of −1.07 MPa occurring at the lower part of surge shaft.The stress distributes equably in different layer of shaft.But at the variable cross-section of surge shaft and near the bifurcated pipes, affected by the junctions of chambers, the stress concentration phenomenon is obvious.The concrete showed the risk of cracking.The distribution rule of section axial force and bending moment is similar to the rule of stress; the maximum also occurs at the   mid-lower part.The tensile stress at the crown and base angle of bifurcated pipes is comparatively great, which partially exceeds the tensile strength of concrete.The concrete showed the risk of cracking.The value of compressive stress is generally small, and stress concentration only occurs at the sidewall near the junctions of chambers.In general, the value of the first and third stress in most areas is not large.The structure stability of the surge shaft is guaranteed.
Affected by the hydrodynamic pressure at the bottom of impedance board, the tensile stress at bottom is greater than the one at top; the maximal tensile stress reaches 2.56 MPa, situated at the downstream junction of impedance board and surge-shaft sidewall.Stress concentration appears at the orifice of impedance board, and the maximal tensile stress reaches 2.2 MPa, situated at the downstream edge of orifice.Under external water pressure of overhauling condition, the upside of surge shaft is mainly under tension, and the underside is mainly under compression, and the values of stress, internal force, and displacement are small (Figure 10).However, the tensile stress and compressive stress at the base angle of bifurcated pipes and near junctions of chambers are comparatively great.The impedance board is mainly in compression.Tensile stress occurs diminutively.The values of stress and distortion in whole impedance board are not big.This indicates that the external water pressure does not take controlling effect in the reinforcement of well bore and impedance board but only in the reinforcement of base angle of bifurcated pipes and junctions of chambers.

Distortion Cracking Characteristics and Reinforcement
Demands.Under internal water pressure, the displacement of upside surge-shaft lining is small and distributes from 0.18 mm to 0.56 mm, the displacement of downside lining is greater and distributes from 0.34 mm to 0.71 mm, and the maximal displacement occurs at the orifice of impedance board and reaches 1.49 mm.Seen from the displacement  vector (Figure 11), the displacement of impedance board under hydrodynamic pressure is greater than other parts, and the lining distorts downwards and outwards under effect of deadweight and internal water pressure.It means the displacement is less affected by the constraint of surrounding rocks to sidewall.The crack width and reinforcement ratio of lining in internal and external water conditions are shown in Table 2. Due to the initial gap elements set outside the crown lining, before lining contacts with surrounding rocks, the outside surface of lining is equivalent to free face; at this moment, the concrete is easy to crack under internal water pressure, and the calculated cracking width of draft tube lining and crown lining of tailrace tunnel is around 0.07 mm (Figure 12), but the cracking width of spandrel lining is smaller.At the upside of surge-shaft lining, the cracking width is small, but at the downside it is greater and distributes from 0.03 mm to 0.09 mm.The maximal cracking width is 0.15 mm and occurs at the downriver junction of impedance board and surge-shaft sidewall.Under external water effect, the structure is mainly in compression with rare cracks appearance.The calculated maximal reinforcement ratio under variable  operating conditions is 1.85%, also located in the downriver junction of impedance board and surge-shaft sidewall.

Conclusions
By the 3D nonlinear FEM calculation, the deformation property and stress distribution rules of lining structure in underground surge shaft and bifurcated pipe are analyzed, and the main conclusions are drawn as follows.
(1) By introducing the elastic-plastic damage constitutive relations of rock mass and the implicit bolt and anchor cable element, this paper properly simulates the process of excavation, support, and lining buildup of the surge-shaft structure; the excavation and support effect of chamber could be effectively reflected.Comparison is made with the calculated results of considering concrete lining as linear elastic material and nonlinear material with cracking, which indicates that the latter is more efficient to reflect the real stress condition of lining, and the calculating method of concrete lining structure with cracking is shown to have practical significance.
(2) The 3D interface element including the interface normal and tangential nonlinear property is clear and definite in concept, which could effectively reflect the nonlinear relation of interface material and has good practicability.Its application in project case shows that the combined bearing and interaction between lining structure and surrounding rocks could be reasonably simulated.
(3) Simulation and calculated results of initial gap between lining structure and surrounding rocks show that the initial gap has great influence on the combined bearing of lining structure and surrounding rocks, the crack distribution, and width of lining structure, so good grouting treatment is necessary.Besides, the simulation method is explained to be effective to reflect the real structure stress condition and make the calculation results more practical.
(4) The calculation and analysis show that the internal water pressure in operating period is the primary controlling condition, and the maximal stress occurs at the orifice of impedance board and the joint between impedance board and downstream sidewall, where it should be focused on in the reinforcement designs; the base angle of bifurcated pipe and the junctions of chambers under external water pressure should also be taken into account.Under design loading, the stress of concrete lining basically meets the material resistance requirements, the maximal crack width is in allowable range, and the deformation of the lining and impedance board basically meets the structure service demands, which shows that adopting concrete lining structure in multifork cylinder surge shaft is absolutely feasible.

Figure 4 :
Figure 4: FEM mesh of surge-shaft lining and impedance board.

Figure 5 :
Figure 5: The first principal stress contrast figure.

Figure 6 :
Figure 6: The third principal stress contrast figure.

Figure 7 :
Figure 7: The first principal stress contrast figure.

Figure 8 :
Figure 8: The third principal stress contrast figure.

Figure 9 :Figure 10 :
Figure 9: The first principal stress distribution under internal water pressure.

Figure 12 :
Figure 12: Crack width distribution of typical sections under internal water pressure.

Table 1 :
Mechanical parameters of materials.

Table 2 :
Crack width and reinforcement ratio of lining in different conditions.