Influence of the cracking consideration in the displacements of a shallow tunnel lined in steel-reinforced concrete: numerical model

Alegre, RS, Brasil Abstract: The study of shallow tunnels introduces different aspects when compared to deep tunnels analysis. Among these characteristics can be emphasized the non-uniform shape of the deformed cross-section and the impossibility of some simplifications such as the consideration of homogeneous stresses around the excavation. Since the tunnel lining is subject to combined bending and compression, shallow tunnels are more susceptible to the development of final tensile stresses (not founded in deep tunnels) and the consequent concrete cracking. This paper presents a numerical simulation in finite elements with the ANSYS software. The purpose of the study is to analyze the cracking influence in the displacements of shallow tunnels, treating the concrete behavior by three different models: elastic, viscoelastic and viscoelastic method, concrete, viscoelasticity. Resumo: O estudo de túneis superficiais apresenta peculiaridades se comparado à análise de túneis profundos. Dentre estes aspectos, pode-se enfatizar o formato não uniforme da seção deformada e a impossibilidade de adotar algumas simplificações, como a consideração de tensões uniformes ao redor da escavação. Uma vez que o revestimento do túnel está sujeito à flexo-compressão, túneis superficiais são mais suscetíveis ao desenvolvimento de tensões finais de tração (não encontradas em túneis profundos) e a consequente fissuração do concreto. Este artigo apresenta uma simulação numérica pelo método dos elementos finitos, realizada com o software ANSYS. O estudo tem o objetivo de analisar a influência da fissuração nos deslocamentos de túneis superficiais, ao tratar o comportamento do concreto por meio de três diferentes modelos: elástico, viscoelástico e com Palavras-chave: túneis superficiais, fissuração, método dos elementos This paper presents a finite element model elaborated with ANSYS software. The tunnel excavation and lining placement processes are simulated with the tool of activation and deactivation of elements, step by step. Soil mass behavior is represented by an elastoplastic model with the Mohr-Coulomb plasticity criterion. Three different models represent the concrete behavior: linear elastic, viscoelastic and viscoelastic considering the cracking effects. Since these two last-mentioned models are not previously available in the software chosen, ANSYS's User Programmable Features (UPF) tool is used, based on the studies of Quevedo and Schmitz Additionally,


INTRODUCTION
The strains resulting from an excavation process reduce the existing stresses in the soil mass and induce forces in the lining that correspond to a part of this initial geostatic pressures [1], [2]. So, the design of underground facilities requires the correct evaluation of both the soil mass deformations and the lining stress levels. The magnitude of these solicitations depends on several factors, such as the excavation method, the lining type, the material characteristics and, mainly, the interaction of the lining and the soil mass.
In tunnel engineering, there is a distinction between treating a tunnel as a shallow or as a deep one, since they present different characteristics. Stress and strain fields are the main aspects that distinguish then. Dealing with shallow tunnels, located at a small depth compared with their cross-section size, the complexity of the problem increases, since this type of tunnel is influenced by the free ground surface above it. Shallow tunnels present, then, an ovalized cross-sectional deformed shape, different from the uniform deformed shape usually seen in deep tunnels [3].
Since shallow tunnel linings are subjected to lower compression levels than deep tunnels, it is important to consider the actual coefficient of earth pressure at rest ( K ), which modifies the value of the horizontal pressures in relation to the vertical pressures, considering a relation value different from one. On the other hand, treating the horizontal pressure acting on the soil as equal to the vertical pressure (simplification often adopted in deep tunnels analysis) may have a significant impact on the results. The actual K coefficient causes some cross-sectional regions less compressed than others. As the lining is subject to bending-compressive stresses, a shallow tunnel becomes more susceptible to the development of tensile stresses and consequent concrete cracking, which may reduce its stiffness and induce higher strains.
This paper presents a finite element model elaborated with ANSYS software. The tunnel excavation and lining placement processes are simulated with the tool of activation and deactivation of elements, step by step. Soil mass behavior is represented by an elastoplastic model with the Mohr-Coulomb plasticity criterion. Three different models represent the concrete behavior: linear elastic, viscoelastic and viscoelastic considering the cracking effects. Since these two last-mentioned models are not previously available in the software chosen, ANSYS's User Programmable Features (UPF) tool is used, based on the studies of Quevedo [4] and Schmitz [5]. Additionally, the steel reinforcement is represented by an embedded reinforcement element, with perfect elastoplastic behavior.
The purpose of the performed analysis is, thereby, to determine the difference between the consideration of concrete behavior with the help of the three mentioned models. In these analyses, tunnel depth is also varied to establish its influence. Furthermore, the impact of other involved parameters (unlined length ( 0 d ), concrete strength ( ck f ) and soil mass Young modulus ( m E )) is studied by considering a given depth value and the concrete behavior as viscoelastic with cracking.

MECHANICAL BEHAVIOR OF TUNNELS
In the study of mechanical and structural behavior of a tunnel, it is particularly important to consider in the analysis the constructive phases, related to the soil mass excavation. During the excavation process, the region near the tunnel advancing cross-section (called excavation face) is subjected to higher stresses and strains gradients. Hence, this part of the tunnel consists of an influence zone that starts from a cross-section located inside the unexcavated soil mass, where the radial displacements are small, to another section located in the excavated part, where the displacements reach maximum values [6]- [8].
The radial displacements ( r u ) define the so-called tunnel convergence ( i U ), shown in Equation 1. This parameter is described by the relation between these displacements in each section and the tunnel radius ( e R ).

( )
(1) Figure 1 shows the convergence curve of a circular tunnel, relative to the excavation sections. From it, we can observe that the displacement values are null in a region of the unexcavated soil mass far from the excavation face (undisturbed). They present a maximum gradient in the region around the excavation face and show the greatest magnitudes in an excavated region far from the front (in which the tunnel reached the equilibrium).
Numerical methods based on the finite element method (FEM) are an important tool to determine soil mass displacements and to design the tunnel lining, since they present several advantages over analytical and empirical methods. The numerical methods allow the use of nonlinear and distinct constitutive models for the soil and the lining, the modeling of complex geometries and the consideration of boundary and loading conditions like those found in the field. Additionally, numerical simulations can reproduce the excavation process and the lining installation by the activation and deactivation method, as done in the studies of Bernaud [9], Fiore et al. [10] and Quevedo [4].
In this model, the excavation and lining placement sequences are simulated by changing the stiffness value of the affected elements (soil or concrete) at each excavation step (face advance). To represent the soil mass removal, the stiffness of the excavated elements is reduced. In order to simulate the lining placement (at a specified distance 0 d from the advancing face), the mechanical characteristics of the respective elements, which were previously related to the soil, are replaced to those of the concrete lining [11]. The typical mesh employed in the models which are used in the construction stages analysis is shown in Figure 2 (illustrative).

Tunnel lining
The lining is the structural system installed to provide the necessary support to the tunnel. Thus, it confers stiffness and assists the soil mass in the stabilization of strains during the construction and the different structural stages. It acts also with other functions such as limiting water infiltration and representing the basis for the final inner tunnel surface [12]. In general, this is a hyperstatic problem, since the forces acting in the lining, the stresses resulting from these actions, and the resulting displacements are interdependent and related to the soil-lining interaction.
Several factors are associated with the stresses developed in a tunnel lining as, for example, the lining stiffness. For a better understanding of stiffness influence, Peck [13] suggests that two main situations should be assumed, both for a circular tunnel with lining. If the tunnel lining is perfectly flexible but able to withstand the radial compressive pressures, tangential and shear stresses would not appear, neither bending moments. On the other hand, assuming this lining to be perfectly rigid, the pressures would cause bending moments due to the lining resistance to the imposed forces. In practice, however, the actual lining stiffness is acting between the two proposed conditions. In this way, the tunnel equilibrium cannot be reached only by changing its cross-section size, and this distortion may induce the development of bending moments.
The different structural behavior of a tunnel lining can be explained as a case of bending-compressive stresses. For deep tunnels, the resulting main stresses are purely compressive and more uniform than those found in shallow tunnels. As a result, the bending causes essentially areas with a higher concentration of stresses (usually compressive ones). On the other hand, in shallow tunnels, due to the smaller magnitude of the compressive forces, this bending can cause final tensile stresses.

Shallow tunnel conditions
The distribution of stresses in the soil mass, the magnitude order of the displacements on the ground surface and the degree of displacements symmetry above and below the tunnel cross-section present some aspects that distinguish the tunnels between shallow or deep. In practice, it is possible to differentiate these two tunnel types through the relationship between their depth ( H ) and their diameter ( D ) and, although there is no accordance about the limit value, the adoption of the relation / H D 10 < (as in Benamar [14] and Ferrão [15] studies) seems appropriate to consider the tunnel as a shallow one.
The stresses field developed around a shallow tunnel is not purely radial like in deep tunnels, due to the influence of the proximity to the free ground surface. Consequently, the cross-sectional deformation of a shallow tunnel is not merely uniform and is not related only to the cavity volume change. Actually, in accordance with Pinto and Whittle [3], two other components define the final deformed state: distortion (or ovalization) and vertical translation. The difference in the deformed shape from shallow to deep tunnels is shown in Figure 3. Moreover, regarding shallow tunnel analysis, it is most adequate to consider as acting in the soil mass the actual variable pressures than the geostatic-hydrostatic ones, taking into account the value of the coefficient of earth pressure at rest ( K ), which is actually different from the unit. This coefficient modifies the value of the horizontal pressures in relation to the vertical ones, causing non-uniform displacements and stresses, which vary along the tunnel cross-section (causing zones of stress concentration and relaxation) and bending moments in the tunnel lining [16]. As the soil weight is smaller at lower depths, the compression level at the lining is, consequently, lower. Then, the anisotropy induced by K can cause the development of tensile stresses.
Jensen [17] carried out analysis concerning circular shallow tunnels lined in concrete (soil mass and lining with linear elastic behavior), in which the influence of both depth and earth pressure coefficient in the appearance of tension stresses were studied. The obtained results indicated that tunnels with a smaller / H D relation present more areas of tensile stresses compared to deeper tunnels. The same conclusion was obtained on the variation of the K coefficient: the smaller and distant from the unit, more tensile areas tended to form.
In shallow tunnels, the change in the stress state caused by the excavation also induces another displacement type, called surface settlement. Among the methods of analysis and evaluation of surface settlements, is cited the empirical method given by Peck [13], which indicates that the shape of these displacements resembles a reversed Gaussian curve. Thus, the maximum displacement is precisely on the axis of the tunnel cross-section and the magnitude is greater the closer the tunnel is to the surface.

CONSTITUTIVE MODELS OF MATERIALS
This item aims to describe the constitutive models chosen to represent the materials that compose the tunnels: the soil mass and the steel-reinforced concrete lining.

Soil mass
To represent the soil mass behavior, it is employed the constitutive model of Mohr-Coulomb, often adopted as a resistance criterion in the geotechnical engineering. A version of this model is provided by ANSYS [18] and can be applied by the insertion of the soil parameters as the friction angle ( ϕ ), the cohesion ( C ) and the dilatancy angle (ψ ) of the material.

Lining: steel-reinforced concrete
This research aimed to contemplate two behaviors presented by the concrete: one model represents the timedependent strains and the other one represents the low resistance to tensile stresses. Thus, to implement a viscoelastic model with the consideration of cracking effects, it was necessary to use the ANSYS customization tool, the User Programmable Features -UPF, through the UserMat subroutine, specific to customize material behavior. For this, the procedures presented by Quevedo [4], Schmitz [5], Lazzari [19] and Quevedo et al. [20] were followed. Further details on formulations and solution features can be found in the studies cited above.
The time-dependent behavior of the concrete can be explained in terms of the creep and shrinkage phenomena. The first one concerns the continuous and gradual strains increase under constant stresses and can be separated, as in the case of the Solidification Theory proposed by Bazant and Prasannan [21], in a portion dependent on the loading age and another part dependent on the concrete age, the same distinction made by the formulation given by the Comité Euro-International du Béton (CEB-FIP MC90) [22]. The second phenomenon, called shrinkage, refers to the material volume reduction given by the gradual water loss remaining in the capillary vessels inside the concrete, which it has not been completely used in the cement hydration reactions.
The Comité Euro-International du Béton [22] evaluates the total strain in the age t of a concrete element uniaxially loaded at age 0 t with uniform stress ( ) in which ( ) The strain depending on the stress can be expressed by Equation 3: , The It is worth noting that, in accordance with Quevedo et al. [20], the use of the Comité Euro-International du Béton [22] formulation is due to the fact that the creep model fits into the Bazant and Prasannan [21] Solidification Theory. This facilitates the numerical solution of the concrete time-dependent behavior considering aging as an isolated factor (related only to the solidified concrete volume over time).
In turn, the concrete characteristic of not resisting tensile stresses well, since its tensile strength is much lower than its compressive strength, can result in cracking even at low stress levels (to which shallow tunnels are generally subjected, for example). To verify if cracking occurs, the stress level of the integration points is evaluated by two procedures. It is verified whether the integration point stresses reached the yield surface, which in this research is defined by Ottosen [23]. Then, it is also evaluated if cracking or crushing failure occurs (the cracking failure occurs if the first principal stress of the sample point is equal or greater than the half of the concrete average tension strength).
The cracking consideration (for the points which have reached the mentioned criteria) is evaluated through a model of distributed cracks. So, only material properties are modified, without being necessary to change the finite element mesh. In this case, as in Schmitz [5] and Lazzari [19] studies, the crack is considered to be formed in a perpendicular plane to the main stress direction. Furthermore, the longitudinal and transverse elasticity modulus are reduced, in addition to neglecting the Poisson effect.
Dealing with reinforced concrete, when cracking occurs, the concrete between cracks continues to resist certain tensile stresses. The adherence effects between the concrete and the steel bars contribute to the total structure stiffness. This phenomenon is known as tension stiffening effect. One of the ways to represent this behavior in the model is to modify the stress-strain curve of the concrete, which becomes linear with softening, as shown in Figure 4. A constitutive relation, proposed by Martinelli [24], describes this behavior and is given by Equation 6: .
in which c σ = concrete stress; ci E = tangent elasticity modulus; t ε = tensile nominal strain in the cracked zone; c ε = concrete strain; t σ = tensile stress in the cracked zone; and cTU ε = tensile strain limit which defines the end of the softening portion.
In addition to the tension stiffening effect consideration, a gradual reduction of shear stresses in the crack plane is also evaluated, as suggested by Hinton [25]. To represent this phenomenon in an approximate way, the material transverse elasticity modulus is multiplied by a reducing factor.
The described formulations considering both the viscoelastic behavior and the cracking effects in the concrete are then added to the UserMat subroutine, which is called by ANSYS iterative procedure for each integration point. The main program computes the total stresses, deformations, and increments of total deformations in the current time increment and the subroutine returns the update stresses according to the implemented instructions [20]. Shortly, within the subroutine, after stress determination considering the time-dependent effects, the cracking checking is done. If it is reached, the mentioned modifications are considered, and the stresses are updated again. If it is not reached, the subroutine saves the results and the program goes on to the next point of the iteration procedure. Finally, the behavior of the steel reinforcement bars of the lining, which resist only axial forces, can be described by a uniaxial model. Thus, the steel is represented by a bilinear stress-strain diagram and has a perfect elastoplastic behavior, although a small hardening factor is applied to avoid numerical errors, as recommended by Schmitz [5].

NUMERICAL ANALYSIS
This item presents details about the numerical modeling, describing the finite elements chosen, the mesh adopted and other characteristics of the model. Afterward, the results of the performed simulations are shown.

ANSYS numerical modeling
The solid finite element adopted in the model to represent both soil and concrete is the SOLID185, which has eight nodes with three degrees of freedom per node (translation in X, Y and Z). To represent the steel bars, the element chosen is the REINF264, an embedded reinforcement element which has only axial stiffness and could be placed in any orientation within the base element. In the modeling, the rebars (two layers) are positioned in the middle of the base element, both on longitudinal and circumferential directions.
The numerical simulation of the tunnel construction process is done by 37 excavation steps, with the length, each one, of one-third of the external tunnel radius, so that the final model has an excavated part (more refined mesh) and an unexcavated part. The model has 15678 finite elements. The excavation process, simulated by the element's activation and deactivation tool, is done with the ANSYS Birth and Death commands. First, it is generated a duplicate mesh in the lining region; then, the elements that represent the concrete are soon disabled. Finally, step by step, the soil mass elements are deactivated, as those of the lining are reactivated.
It is important to emphasize that, as a simplification, the speed and the length of the excavation are equal in each step, disregarding possible time intervals without excavation (that occurs in practice) and other variations. The soil mass, in turn, is treated as a single material, homogeneous and isotropic, without the presence of heterogeneous layers. , with γ the soil specific weight, H the tunnel height, and K the earth pressure coefficient) are also indicated. Moreover, it is applied a symmetry condition in the ZY plane. Consequently, only half of the model is simulated.
In addition to the external loads, it is necessary to apply an initial stress condition with the same values of v P and h P , as the soil elements weight. This procedure is done to assign a null initial strain state, with no deformations before the excavation. At the same time, the stress state will be the geostatic one.

Results
Some simulations are carried out to compare the results obtained by altering the behavior models of the lining material. Thus, it is expected to analyze the influence of considering or not the cracking of the concrete rather than contemplating only an elastic or viscoelastic behavior. The objective is to obtain the displacements in a region far from the excavation face, in three points of the tunnel cross-section: at 90º, 180º and 270º, for three distinct lining constitutive models.
For the soil mass, a plastic behavior is admitted by assuming the Mohr-Coulomb criterion, as already explained in item 3.1. For the lining, it is expected to evaluate the evolution of the tunnel convergence as the model changes between elastic, viscoelastic, and viscoelastic with cracking (all these with steel reinforcement). The tunnel depth will also be varied so that its influence on the results can be analyzed.
In order to intensify the solicitations undergone by the lining, which can contribute to the cracking occurrence in the concrete, the distance of the lining placement to the excavation front (or unlined length) is considered null in the examples. This means that the lining is placed immediately after the soil mass excavation.
The required parameters for the analysis of a circular tunnel are shown in Table 1. It is worth to mention that the examples consider a fictitious tunnel, with reduced dimensions. Thus, the results represent only an assessment of the procedures adopted and an indication of behavior.  Tables 2, 3 e 4, where the convergence, in %, is given in three points of the tunnel cross-section, for each concrete behavior model. grow from the elastic model to the viscoelastic one, and also increase when considering cracking. In the example with / , H D 5 = however, the difference only occurred from the elastic model to the viscoelastic, without considerable changes in the displacements when considering cracking. This occurs due to the lower pressure levels at which the tunnel is subjected in this smaller depth ( H 10m = ), not causing sufficient tensile stresses for the concrete to crack. In the same way, it can be seen that this difference considering cracking is greater when the tunnel is at a higher depth and subjected to higher pressures (tunnel case with H 20m This consideration is important, since tunnels with larger diameters than the studied ( D 2m = ) may be located at higher depths than those analyzed and still be considered as a shallow tunnel.
The relative difference, from the elastic model to the viscoelastic model and from the viscoelastic model to the viscoelastic model with cracking, for each relation example in the three analyzed positions, is shown in Table 5. Another example of the influence of the concrete behavior model can be observed in the surface settlement results. Figure 6 shows the transverse settlement basin for the tunnel with / .

H D 8
= Since this tunnel has not a considerable ground cover above it, the magnitude of the settlement is small in the three models of concrete behavior. However, it is still possible to visualize the difference in the results according to the model (in terms of maximum settlement: . %

07
from the elastic to the viscoelastic model and . % 7 07 from the viscoelastic to the viscoelastic with cracking). Regarding cracking effects, it is possible to contemplate the concrete behavior in tension, defined previously by Equation 6 and Figure 4. Before cracking, stresses increase linearly with strains. Then, once cracked the concrete, these stresses decrease for greater strains. Hence, the tendency for the cracked model is lower stresses for the lining points with greater strains, which cracks first than others. Figures 7, 8 and 9 illustrate the circumferential stresses ( θθ σ ) in the tunnel lining, from the elastic, viscoelastic and viscoelastic with cracking models, respectively. The stresses in this direction have the highest magnitudes. From the figures, the aspects highlighted above can be observed. In elastic and viscoelastic models, the stress distribution is the same, but the magnitudes change. Regarding tension, the maximum value decrease from .
3 87MPa to . 2 72MPa . On the other hand, from viscoelastic to viscoelastic with cracking model, in addition to the maximum value going from .
2 72MPa to . 1 54MPa , the location of this maximum value also changes. This is because in the previous location of the maximum tensile stress (exactly at 90°), the concrete has already cracked, and its tensile behavior is in the descending part of the curve of Figure 4.   For better understanding the changes, Figure 10 schematically illustrates the tensile concrete behavior for a general section. First, the tensile part behaves linearly, with the maximum tensile stress at the bottom of the section. Upon reaching the concrete tensile strength and, therefore, the deformation referring to the beginning of crack formation, this stress decreases. The decline remains until it reaches the maximum deformation to consider the collaboration of concrete between cracks, equal to . 0 001 . Figure 10. Tensile stress variation in a general concrete section.
In order to elucidate the difference of behavior for the tunnel with / H D 5 = in relation to the others, the stresses ( θθ σ ) according to the model are shown also for this example of tunnel (Figures 11, 12 and 13). From the results, it can be verified the change only from the elastic to viscoelastic model, with the stresses decreasing. Considering cracking in concrete behavior does not cause variation in stresses (there is only a residual difference). Since in this depth the tunnel is submitted to low levels of stresses, the lining does not crack.   A simplified parametric analysis is also performed with the aim to study the influence of some parameters in the final convergence magnitude when considering the cracking in the lining model. The properties of the tunnel and the materials are the same employed in the previous simulations and shown in    As exposed by the results shown in the figures above, the convergence is greater for larger unlined lengths. This is related to the contribution of the lining to the tunnel support: the longer it takes to place the lining, the more the tunnel deforms. On the other hand, this greater strain of the soil mass implies smaller pressures acting on the lining, which attenuates the appearance of cracking.

Em (MPa)
Furthermore, the convergence values decrease when both concrete and soil mass stiffness increase. When the material is stiffer, it can support more pressure and, as a result, fewer displacements are found in the cross-section. Cracking also presents consequences: the smaller the ck f , for example, the concrete tensile strength is also smaller, and, accordingly, this concrete is more propitious to crack, causing greater displacements. However, the influence of these parameters is correlated, since the relative stiffness between the soil mass and the lining also affects the stress distribution, so that the greater the concrete stiffness in comparison to the ground, the higher tensile stresses may form.

CONCLUSIONS
There are several peculiarities of tunnels closest to the ground surface. Among these characteristics, is the great unevenness of the stress field around the excavation and the possibility of the appearance of tensile stresses in the lining, that could induce the concrete to crack. So, the present research deals with the structural behavior of the reinforced concrete lining of shallow tunnels, with emphasis on the cracking analysis and the impact obtained on the displacement results when considering its effects.
From the results obtained, it was possible to observe that the consideration of cracking affects the behavior of reinforced concrete tunnel linings, causing larger displacements in the tunnel cross-section than those calculated considering the concrete with elastic or viscoelastic behavior. For the example with greater depth ( / H D 10 = ) and, therefore, submitted to larger loads, this difference reached . % 16 3 in one of the cross-section points, analyzing the final convergence at a point distant from the excavation face. The maximum surface settlement value was also modified by changing the model: . % 7 07 greater when considering cracking. Changes in stress behavior were also demonstrated. The tunnel with / H D 8 = shown a decrease in the maximum value of stresses by varying the concrete model. In addition, when considering cracking, the maximum stress location also changes, indicating lower stresses for the lining points with greater strains, which cracks first than others. For the tunnel with / H D 5 = , the stresses are the same regarding or not the cracking in the model, since the loads in this depth are not enough to the concrete to crack.
By setting other parameters and changing the unlined length ( 0 d ), the characteristic compressive strength of concrete ( ck f ) and the Young modulus of the soil mass ( m E ), convergence results were higher for larger unlined lengths and smaller values of ck f and m E . It is worth noting that these influences are interdependent because, while they interfere directly in the magnitude of the displacements, they also interfere in the crack formation. Cracking, in turn, causes greater strains and, consequently, greater displacements.