Simulation of fracture propagation in fibre-reinforced concrete using FDEM: an application to tunnel linings

The application of the combined finite-discrete element method (FDEM) to simulate fracture propagation in fibre-reinforced-concrete (FRC)-lined tunnels has been investigated. This constitutes the first attempt of using FDEM for the simulation of fracture in FRC structures. The mathematical implementations of the new FDEM joint-element constitutive model are first introduced, and the numerical model is then validated comparing the results for plain and FRC beams with three-point bending experimental data. The code has also been applied to two practical tunnel design case studies, showing different behaviours depending on the type of concrete and shape of tunnel section. The FDEM simulations of the linings are also compared with results from a finite element code that is commonly used in the engineering design practise. These results show the capabilities of FDEM for better understanding of the fracture mechanics and crack propagation in FRC tunnels. A methodology for directly inferring the numerical parameters from three-point bending tests is also illustrated. The results of this research can be applied to any FRC structure.


Introduction
Further developments in fibre-reinforced concrete (FRC) technology, with the utilisation of a variety of fibre materials (e.g. steel, plastic, carbon, etc.) and geometries, have increased the interest in FRC structures. FRC is now widely used for tunnel linings. Following the collapse of the Heathrow Express Rail Link tunnel in 1994 (and 39 other major fibre-reinforced sprayed tunnel collapses worldwide), a review by the Health and Safety Executive [1] concluded that some safety-critical aspects of fibre-reinforced tunnel design and construction were poorly understood. Moreover, the complex interaction of fibres and concrete for the cracking performance of FRC for its structural design still need to be better understood and quantified. This opens the door to B A. Farsi ado.farsi@imperial.ac.uk 1 Applied Modelling and Computation Group, Department of Earth Science and Engineering , Imperial College London, South Kensington Campus, London SW7 2AZ, UK 2 Bedi Consulting Ltd, London, UK 3 London Underground Limited, London, UK the application of computational methods for simulating the mechanical behaviour and fracturing of FRC structures.
Interest in simulating fracture propagation extends across a variety of scientific and engineering fields, such as structural analysis, material design, nuclear waste disposal risk assessment, oil and gas reservoir engineering and subsurface ore mining [2]. Numerical simulations are performed in order to predict the formation and behaviour of these fracture systems, due to the geometric and physical complexity inherent in fracture phenomena.
Two main modelling approaches can be identified in the literature for fracture analysis: discrete crack and smeared crack models, also known as geometric/non-geometric, or grid/subgrid methods. They were introduced in the late 1960s by Ngo and Scordelis and Rashid in application to concrete structural analysis [3]. The smeared crack model is based on the assumption that in concrete, due to its heterogeneity and the presence of reinforcement, many small cracks nucleate which only in a later stage of the loading process link up to form one or more dominant cracks. Since each individual crack is not numerically resolved, the smeared crack model captures the deterioration process through a constitutive relation, thus smearing out the cracks over the continuum. In this kind of analysis, cracks are represented as an isotropic or anisotropic damage concentration band within a mesh element from which fracture geometry can be inferred, instead of being explicitly defined [4]. Similarly in phase-field approaches, fracture surfaces are approximated by the definition of a field that smears the boundary of the cracks over a small region [5]. One of the major advantages of phase-field approaches is that the evolution of fractures follows from the solution of a coupled system of partial differential equations. This drastically reduces the complexity of the solver that is required to calculate the simulated problems. In contrast, the discrete crack model represents cracks discretely and aims to simulate the initiation and propagation of dominant cracks. This kind of analysis can be performed with different approaches, e.g. boundary element-based methods [6], peridynamics [7,8], finite element simulations [9,10] extended finite element method [11], discrete element method (DEM) [12][13][14][15][16] or combined finitediscrete element method (FDEM) [17]. An important aspect to underline is that the great majority of geometric methods do not maintain a representation of the fractures separate from the mesh and rely on mesh editing techniques, such as in-situ insertion of new crack nodes, edges and faces, to capture mesh growth.
There are also several approaches presented in the literature aimed at describing solid fragmentation using DEM. In [14], numerical simulations of uniaxial compression and cutting processes of rock have been presented employing packed discrete elements bonded together to represent the bulk rock material. Similar examples are presented to describe tunnelling, flexural and Brazilian tests in [16]. The standard approach consists of defining a packed structure of particles (generally spheres) with a determined particle size distribution and then defining laws for the contact and interactions between particles on the basis of parameters such as penalty numbers, stiffness of the bonding and to obtain forces that, once they are applied to the elements, define the movement of the simulated bodies with Newton's laws of motion.
Algorithms for FDEM simulations started to be proposed from the 90s. Extensive developments and applications of the FDEM method have been carried out after the release of the open source Y-code in [17], and different versions have been released, including the code developed from the collaboration between Queen Mary University and Los Alamos National Laboratory [18,19], the Y-Geo and Y-GUI software that have been developed by the Geomechanics Group led by Giovanni Grasselli at Toronto University [20,21], and Virtual Geoscience Simulation Tools (VGeST) released by the Applied Modelling and Computation Group (AMCG) at Imperial College, London [22,23]. Recently, the AMCG has upgraded and renamed VGeST as 'Solidity'. A commercial FDEM code developed by Geomechanica (www.geomechanica.com) has also been released in Canada, although its application has been limited to modelling geomaterials. While the first Y-code employed finite strain elasticity coupled with a smeared crack model to capture deformation, rotation, contact interaction and fragmentation, the AMCG has greatly improved the code, implementing a range of constitutive models in 3D [24,25], thermal coupling [26], parallelisation and a faster contact detection algorithm [27].
Several models for the simulation of fibre-reinforced concrete have been proposed [28][29][30]. Most of the techniques that have been presented, instead of inserting discrete fractures, represent the concrete failure in a continuous fashion with plastic deformations. A separate numerical representation of the fibres and concrete matrix has been presented in [31], where the spatial location and orientation of the reinforcements have to be known or assumed a priori. Another approach that has been proven to be effective in simulating fibres in cracked concrete is to smear their effects in the traction versus crack opening law, which governs bridging stress across the crack faces of plain concrete [31]. In [32], a cohesive fracture model for the simulation of FRC has been presented. One of the limitations of this approach is that the locations of the fractures must be known a priori in order to insert cohesive elements in the right locations in the model. The advantage of using FDEM is that fractures can initiate anywhere in the simulated system due to stress criteria, without the need to know the failure mechanism of the simulated structure a priori, (i.e. in FDEM, there is no requirement to seed pre-existing flaws or crack tips within regions of the material and thereby constrain crack propagation in the model). Furthermore, the multibody framework of FDEM allows for the fragmented material interactions resulting from cracking to contribute to postpeak strength behaviour, for example when lining failures include compressive shear fracture behaviour and frictional sliding.
FDEM simulations of geological and concrete structures have attracted more and more attention in recent years. This work constitutes the first attempt to extend the FDEM framework to fibre-reinforced concrete structures. This will allow designers to better assess the serviceability and ultimate limit state behaviour of FRC tunnel linings, and it will contribute to the development of detailed design guidance for FRC structures. The mechanical behaviour of FRC is different from the brittle materials that were previously simulated with FDEM codes (mainly rock and concrete). For this reason, a novel joint-element constitutive model that introduces the softening behaviour induced by the presence of the fibres in the concrete matrix was developed. This constitutive model is presented in Sect. 2, and it contains several aspects of novelty from previously proposed cohesive models for FRC, e.g. in [32]. Instead of employing just multilinear functions, this constitutive model utilises a more sophisticated stress-opening curve for the concrete matrix (i.e. the early post-peak behaviour), which is based on heuristic curves parameterised with actual tests on concrete. In Sect. 3, this novel joint-element constitutive model is validated through direct comparison with three-point bending experimental results. In Sect. 4, the model is applied to two different tunnel lining designs, both a circular and a typical SCL lining profile used at Heathrow in London Clay, both loaded up to the ultimate limit state. The results are then compared to those from standard elasto-plastic finite element models that are commonly used in the engineering design practise but lack discrete fracture representation. In Sect. 5, an extension of the proposed joint-element constitutive model to represent more complex post-peak behaviours (e.g. hardening) is introduced. A methodology for directly inferring the numerical parameters from three-point bending tests is also illustrated. Section 6 contains a detailed discussion of the results that are shown in this work. The key findings are then summarised in Sect. 7.

Plain concrete
The Y-code, which was presented in [17], has provided the fundamental theoretical aspects of the contact and cohesive crack model employed in this work. The transition from a continuous domain to a discontinuous domain is carried out through fracture initiation and propagation processes. The progressive fracture mechanisms implemented in the code for the description of plain concrete is idealised with a distributed micro-crack zone, a bridging zone and a traction free macrocrack zone. This is based on the assumption that the stressstrain curved consists of a hardening branch (before the peak) and a strain-softening part (where the stress decreases with the strain increasing), as illustrated in Fig. 1.   Figure 2 shows the strain-softening relation that has been implemented in the code through the constitutive law of the joint-elements in terms of stress and displacements. The strain-softening part is also defined in terms of stress and displacements in order to avoid the ill-posedness of the problem generated by a stress-strain definition. The area under the graph in Fig. 2 is the energy release rate (G c ), i.e. the energy dissipated in order to extend the surface of the crack. G c is also equal to twice the surface energy, which quantifies the disruption of bonds that occur when a surface is created. The relationship between stresses and displacements is modelled through a single-crack model: when the size of separation is zero, the bonding stress is equal to the tensile strength ( f t ), implying that the separation begins only after reaching this stress value equal to f t . Once the separation starts to increase, there will be a decrease in bonding stress. When it reaches a limit value of separation (δ p ), the bonding stress tends to zero. In the actual implementation of this model, the separation of adjacent element edges is assumed in advance by introducing joint-elements and describing the topology of adjacent elements with different nodes. As no two elements share any nodes, the continuity between elements is enforced through the penalty function method. Before the bonding stress reaches the tensile strength, its value is given by Eq. (1), where δ t is the separation corresponding to when the bonding stress is equal to the tensile strength (δ t = 2h f t / p), h is the size of that particular finite element and p is the penalty parameter.
function representing an approximation of the experimental stress-displacement curves and the parameters a = 0.63, b = 1.8 and c = 6 are obtained from the interpolation of experimental stress-displacement curves. These parameters only define the shape of the softening curve, which is then stretched depending on the material properties of the material (i.e. the energy release rate and the tensile strength). However, the heuristic curve that was implemented in the first Y-code was derived from direct tension experiments on concrete samples [33]. There are other materials that have been tested to obtain a more representative softening curve shape, such as granite, e.g. see [19]. The variable D is given by Eq. (3), and the complete relationship for the normal bonding stress as a function of separation can be written as shown in Eq. (4).

Fibre-reinforced concrete
The fracture mechanism involved in the fibre-reinforced spayed concrete is affected by the presence of discrete fibres in the concrete matrix. Although the presence of fibres increases the fracture process zone size, neither the tensile strength nor the early post-peak behaviour is assumed to be influenced by the presence of fibres at low volume fractions. Figure 3 shows the modified stress-opening relation that has been implemented in the code to account for the concrete cracking, fibre de-bonding and fibre pull-out. The concrete matrix fracture mechanisms dominate the post-peak behaviour for small crack opening. For this reason, the stress, after reaching the value of tensile strength ( f c t ), follows the relation for the plain concrete with the previous exponential softening curve defined by Eq. (2) up to a crack opening of δ f p , defined in Eq. (5). The value of crack opening δ f p is reached when the fibres in the concrete matrix are activated. At the beginning of the fibre pull-out phase, the bonding stress between two fracture walls is defined by the post-crack residual strength f f t . The stress then decreases linearly with the crack opening δ, up to a value of f f , corresponding to the crack opening limit δ f . The complete relationship for the normal bonding stress as a function of separation is shown in Eq. (6).
The constitutive laws for the joint-elements that have been defined in this section can be represented as the interposition of normal and shear springs between the joint nodes of the elements in contact, as shown in Fig. 4. The grave normal springs have a nonlinear stress-displacement law given by Eq. (4), whereas the shear springs have analogous laws representing shear failures. The normal and shear springs between the joint nodes of the elements in contact are removed once the separation reaches the value δ f , meaning that the fracture has propagated through the edge.

Numerical model
The implemented joint-element constitutive model has been validated by comparing the simulated response with the experimental data for three bending tests on plain and FRC beams reported in [32]. The geometry of the tested samples is shown in Fig. 5a.
The experiments have been modelled with 2D FDEM simulations. The boundary conditions and the 2D triangular mesh are illustrated in Fig. 5b. The simulation is performed with a plane-strain assumption. The top roller is constrained with constant velocity. To reduce the run time of the numerical simulation, the velocity of the constraint is set to 0.01 m/s. Although this loading rate is higher than the one in the laboratory experiment, it induces a quasi-static response in the bar as there is no significant difference between the force applied by the punch and the one applied by the two constraints. To further reduce the calculation time, when the simulation starts, the roller is in contact with the specimen. Since fracture initiates and propagates above the notch, the specimen is discretised with an unstructured coarse mesh at the two edges. Moreover, a fine mesh is employed in the centre of the specimen to better represent both the de-bonding stress during the opening of the crack and the fracture path along the element boundaries. The total number of elements employed in the model is around 2200. The material prop-erties used to describe the three-point bending apparatus are E s =210 GPa, ν s =0.3 and ρ s =7850 kg/m 2 , where E s is the Young's modulus, ν s is the Poisson's ratio and ρ s is the density. The material properties used for the specimens are E c = 26.5 GPa, ν c =0.19, ρ c = 2333 kg/m 3 , f t =3 MPa and G I = 90 J/m 2 . The frictional interaction between the rollers and the samples is modelled using a Coulomb coefficient of friction equal to 0.01. The parameters that describe the FRC post-peak behaviour due to its fibres are f t f = f f = 0.5 MPa and δ f = 3.5 mm. The materiel properties used in the numerical simulations are summarised in Table 1.

Results
Results from the FDEM simulations on plain and fibrereinforced concrete and the corresponding experimental load-crack mouth opening displacement (CMOD) curves of three bending tests on the two types of concrete, are shown in Fig. 6. The numerical results for the plain concrete specimens are compared with the corresponding experimental curves in Fig. 6a. The value of the peak load is correctly represented in the simulated experiment. The numerical model slightly overestimates the load in the initial part of the post-peak load-CMOD curve and underestimates the load after a CMOD of about 0.2 mm. This difference is due to the values of the parameters a, b and c that have been chosen to represent the decay of the concrete stiffness during crack propagation.
Overall, the behaviour of the actual plain concrete specimens is well represented in the corresponding numerical simulation, with an ultimate CMOD of 0.8-0.9 mm. Figure 6b shows the results for the FRC specimens. Here, the numerical results match very well the experimental data, not only in terms of peak load, but also for the residual strength during fibre de-bonding and fibre pull-out. Figure 7 shows the difference in the horizontal stress filed for the two types of concrete for an equivalent fracture width. In the plain concrete beam, the two fracture walls are unloaded, as shown in Fig. 7a. Differently, in the FRC beam shown in Fig. 7b, even though the fracture has almost split the specimen in two halves, the fracture walls are still joined by some steel fibres and the cracked section is under horizontal stress, in order to support the applied load that is transmitted by the punch.

Numerical model
The new joint-element constitutive model has been applied to simulate the response of two FRC tunnels. The geometry and mesh of the simulated circular tunnel are shown in Fig. 8a, whereas the typical NATM tunnel shape, which is Fig. 5 a Geometry of three-point bending specimens that have been tested in [32]. b Mesh and boundary conditions of the same beam that have simulated using FDEM that of the Heathrow Trial Tunnel (HTT) [34], is illustrated in Fig. 8b. The tunnels have been modelled with 2D FDEM simulations, assuming plane-strain conditions. The two tunnels are subjected to increasing in situ stresses: the vertical stress is ramped from 100 kPa to a value of 1 MPa. The horizontal stress is ramped from 40 to 400 kPa. Although the FDEM code is capable of simulating the interactions between the soil and the tunnel linings, the effects of the in situ soil stresses and subsequent passive reaction have been directly applied to the linings as an equivalent surface pressure. This has been done to reduce the run time of the simulations. A fine mesh is employed for the linings to better represent the de-bonding stress during the opening of the crack. The refinement of the triangular mesh employed for the tunnels is illustrated in Fig. 8c. The total number of elements employed in the model is around 36,100. The material properties used to describe the soil (London Clay) are E s =125 MPa, ν s =0.2 and ρ s =2000 kg/m 2 , where E s is the Young's modulus, ν s is the Poisson's ratio and ρ s is the density. The material properties used for the linings have been reported in the previous section. Figure 9 shows the horizontal stress field and fracture surfaces for the numerical simulations of the plain concrete and FRC circular tunnel under ramping in situ stresses. Fractures appear with an in situ vertical stress between 450 and 600 kPa in the inner arch at the top and bottom of the linings. Figure 9c and e shows the fractured sections of the plain concrete lining. After failure, the stress around the cracks drops to zero, meaning that portion of the tunnel is not supporting the in situ stresses. Differently, the fibres in the fractured sections of the FRC lining are still transmitting some stress in order to equilibrate the in situ stresses, as shown in Fig. 9d and f. The horizontal stress at the bottom and top inner arch of the tunnel Fig. 6 Experimental results (dashed black) for three bending tests on a plain and b FRC beams [32]. Numerical response (red) of the simulated a plain and b FRC beams with the joint-element constitutive model implemented in FDEM   Fig. 7 Horizontal stress [Pa] in the simulated a plain and b FRC beams with the joint-element constitutive model implemented in FDEM  Fig. 8 Geometries and finite element mesh of a the circular and b complex-shaped tunnel that have been analysed using FDEM. c Detail of the mesh discretisation used for the complex-shaped tunnel is plotted against the applied in situ vertical stress in Fig. 10. The numerical results show that cracking is more localised for the plain concrete lining (with a crack spacing of about 180-200 mm) compared to the FRC (with a crack spacing of about 120-160 mm). The stress in these most vulnerable sections of the lining increases with the applied in situ stresses, up to the value of tensile strength (3 MPa). At that point, the stress carried by the plain concrete drops to zero, whereas reaches the value of the residual strength (0.5 MPa) in the FRC. Figure 11 shows the results of a similar analysis undertaken with finite elements (FE). This figure shows the areas where the model predicts tensile stresses in the tunnel lining. The finite element model implements the nonlinear stress strain curve of FRC shown in Fig. 12.

Results of circular tunnel shape
Comparing the FE analysis results with the FDEM analysis, it is apparent that the areas and magnitudes of stress predicted by both are comparable. However, the FE model is not able to capture the crack initiation and propagation. Thus, crack widths, spacings and locations are not able to be captured with this 'smeared' approach to modelling cracking. Figure 13 shows the horizontal stress field and fracture surfaces for the numerical simulations of the plain concrete and FRC HTT shape under ramping in situ stresses.

Results of Heathrow trial tunnel shape
Fractures appear with an in situ vertical stress around 400 kPa. The plain concrete lining in Fig. 13a shows a single crack in the centre of the inner bottom arch of the tunnel. High tensile stresses also occur on the extrados of the lining at the 'knees'. The FRC lining in Fig. 13b shows a more distributed crack pattern due to stress redistributions and post-crack ductility provided by the fibres. In this case, the fibres have been activated to keep the structure together. There is a certain degree of asymmetry in the results due to the fact that the computational mesh is not symmetric with respect to the vertical axis. Figure 14 shows the results of a similar analysis undertaken with finite elements (FE), similar to the circular tunnel. As with the circular shape, the FE analysis predicts similar stresses to the FDEM analysis, and again the FE model is not able to capture the crack initiation and propagation. Moreover, the zone of tensile stresses on the extrados at the knees The horizontal stress at the bottom (black) and top (red) inner arch of the tunnel is plotted against the applied in situ vertical stress of the lining in the FE model is slightly above the locations identified in the FDEM simulation. This is likely to be due to the stress redistribution around discrete cracks that is modelled with FDEM code but 'smeared' by the FE code.

Multilinear joint-element constitutive model
In the previous sections, the residual strength of fibres has been implemented in the FDEM code with a linear degra-dation of the joint-element stiffness. This allowed the model to correctly simulate the experimental results presented in [32]. In some cases, the post-peak behaviour of FRC cannot be described with a linear curve, e.g. when it shows some forms of hardening. With a similar approach to the one presented in the previous sections, more complex stress-opening relations can be described with a multilinear degradation of the joint-element stiffness. There is merit in extending the current implementation of the strain-softening joint element constitutive model to a multilinear function, which would allow generalisation of the FDEM approach. The parameters would then be calibrated to any measured FRC response. A proposed approach is described herein. Equation (7) describes a modified multilinear version of the joint-element constitutive model, as shown in Fig. 15a. The f j and δ j parameters can either be directly calculated from uniaxial tensile tests or inferred from load-CMOD curves of bending tests, like the one shown in Fig. 15b. The residual flexural tensile strength f j at any CMOD can be approximated from the Euler-Bernoulli beam theory with Eq. (9), in accordance with EN 14651 [35], where F j are the recorded values of load, b is the width, l is the span length and h sp is the distance between the tip of the notch and the top of the specimen.

Discussion
In the bonding stress model that has been described above, the stress and strain fields close to the crack tip are influenced by the magnitude and distribution of the bonding stress close to the crack tip. In particular, the stress field is influenced by the mesh topology close to the crack tip. For this reason, in order to obtain accurate numerical representations of the stress relaxation during fracture propagation, a sufficiently fine mesh needs to be used. Since each joint-element transmits a constant value of bonding stress, at least four elements should be used to discretise the plastic zone during the crack opening. The size of the plastic zone, which is the portion of material that is deforming under stresses that have reached the value of the tensile strength, is defined by the material properties of the simulated structure (Young's modulus, fracture toughness and tensile strength) and is normally of the order of few millimetres for standard concretes. For many engineering applications, when a system of dozens of metres is simulated, this constraint on the mesh size might be impractical as it would extend too much the runtime of the simulations. For this reason, a coarser mesh is often employed when simulating big systems, reducing the accuracy of the calculations.
The mesh size is not the only constraint on the numerical discretisation of the simulated body. The fracture path is constrained to follow the element edges, where the joint-elements are located. For this reason, since a crack cannot propagate through finite elements but only through their edges, a structured mesh would artificially constrain the cracks to follow some preferential paths. The fracture paths are generally not known a priori, and they should be only governed by the direction of maximum tensile stress at the crack tip. The fracture grows at a crack tip when the stress field reaches the value of tensile strength in the orthogonal direction to the closest joint-element and the kinematics of the structure activated by the external loads needs to allow an opening displacement in that same direction. Unstructured meshes can minimise the mesh dependency on the final crack paths, but the fact that the joint-elements are not necessarily aligned to the direction of maximum stress at the crack tip can cause an artificial overestimation of the structural strength and toughness at the mesoscale.
The code has an explicit solver; in other words, it discretises the continuous time in time steps, and then, it calculates the state of a system at a later time step on the basis of the state of the system at the current time step. This implies that, even if only the status of the system at an exact time is required (e.g. when the tunnel is loaded up to the ultimate limit state), the code needs to calculate the output for every prior time step until it reaches the required one. To give an example: if you want to calculate the state of a modelled system after three seconds from the initial conditions, the total run time is equal to the run time of a single iteration multiplied by the number of time steps before the required one; in this case, three seconds divided by the length of the time step used for the time integration in the simulation (which is normally in the order of the nanoseconds). At the beginning of the simulation, the whole model is at rest, with zero stress everywhere. This means that the in situ stresses need to be applied during the simulation slowly enough to avoid dynamic effects in the simulated system. In other words, to assure the stability of the simulation, the code needs to calculate the whole system for many time steps. In some cases, this makes the run time of a simulation unfeasible.
In spite of the modelling limitations that have been illustrated in this section, the FDEM code is capable of capturing the failure mechanisms of FRC structures with a sufficient level of accuracy, as demonstrated in the previous sections. On the one hand, this alone constitutes a considerable improvement from the potentially over conservative continuum methods, such as FEM. These methods do not explicitly capture fracture initiation and propagation and consequently cannot describe with adequate accuracy the stress redistribution in the structure during failure. On the other hand, the discrete representation of fractured portions in the simulated structures allows a quantitative estimation of the damage under a defined loading condition. Not only is it possible to assess the overall performance of the simulated structures, but it is also possible to extract information such as crack widths and crack spacing. This capability of the FDEM code can also be useful during tunnel design in order to identify the most vulnerable sections of the lining. Moreover, parametric studies can be done to take into account uncertainties during the tunnel construction. Parameters such as the variation in the thickness of the lining and actual density of fibres in the mix can be changed in the simulation to study the serviceability and ultimate limit state behaviour of FRC tunnel linings under different scenarios.

Conclusions
In conclusion, the work that has been presented is the first attempt to use a joint-element constitutive model in FDEM for the simulation of fracture in FRC structures. The new FDEM joint-element constitutive model that has been developed provides a promising prediction of peak and post-crack behaviour of a published data set [32] for plain concrete and fibre-reinforced concrete three-point bending test. The proposed constitutive model is sufficiently versatile to be adapted to any fibre response based on its three-point bending test results.
The FDEM code has been applied to a circular and typical SCL lining profile in London Clay loaded up to the ultimate limit state. The numerical results show that cracking is more localised for a plain concrete lining compared to the FRC. Some of the limitations of the code have been presented and discussed. Although the FDEM code is capable of simulating the interactions between the soil and the tunnel linings, in the proposed examples, the effects of the in situ stresses have been directly applied to the linings as a surface pressure in order to reduce the run time of the simulations.
The simulations that have been shown in this work were limited to FRC tunnels, but the model can be applied to any FRC structure. Continuum methods, such as FEM, do not explicitly capture fracture initiation and propagation and consequently cannot describe with adequate accuracy the stress redistribution in the structure during failure. For this reason, they generally give over conservative results.
The proposed method allows the discrete representation of fractured portions in the simulated structures that enable a quantitative estimation of the damage under defined loading condition. For instance, it is possible to extract information such as crack widths and crack spacing. The novel techniques employed in this study will help to increase knowledge of the serviceability and ultimate limit state behaviour of FRC tunnel linings, and so significantly contribute to the development of detailed design guidance for FRC structures.
The FDEM code can be used during tunnel design in order to identify the most vulnerable sections of the lining and therefore to target improved design efficiency. Moreover, parameters such as the variation in the thickness of the lining and actual density of fibres in the mix can be changed in the simulation to study the serviceability and ultimate limit state behaviour of FRC tunnel linings under different scenarios. More research needs to be undertaken to reduce the run time of the simulations and model explicitly the interaction between the soil and the FRC linings.