{'fraction_non_alphanumeric': 0.0236540202049852, 'fraction_numerical': 0.006702380078765197, 'mean_word_length': 4.548588490770901, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient and realistic 3-D Boundary Element simulations of underground construction using isogeometric analysis

The paper outlines some recent developments of the boundary element method (BEM) that makes it more user friendly and suitable for a realistic simulation in geomechanics, especially for underground excavations and tunnelling. The innovations refer to the introduction of isogeometric concepts, elasto-plastic analysis and the simulation of ground support. The introduction of isogeometric concepts for the description of the excavation boundaries results in less user and analysis effort, since complex geometries can be modelled with few parameters and degrees of freedom. No mesh generation is necessary. In order to consider heterogeneous and inelastic ground conditions a domain discretisation is used but the use of isogeometric methods results in a minimum of user effort. Efficient methods are also presented for modelling rock bolts. A comparison of results of test examples with other numerical methods and analytical solutions confirm the efficiency and accuracy of the proposed implementation. A practical example with a complex geometry is presented.


Introduction
Since the publication of the first paper on the topic [14], isogeometric analysis has gained increased popularity. The majority of applications have been with the Finite Element method (FEM) and much less with the BEM. However, the advantage of the BEM, that requires only the discretisation of the boundary, makes it an ideal companion to Computer Aided Design (CAD). First applications of the isogeometric BEM (IGABEM) were published in elasticity in 2-D [20,21] and in 3-D [19]. Other applications followed (for example see [2], [11]). In [16] the concept of a geometry independent field approximation, which involved a decoupling of the geometry definition and the approximation of the unknown, was first introduced and was later adopted by others [3]. The seamless integration of BEM and CAD was discussed in [17]. In a recent book published on the isogeometric BEM [8] it was shown how geometrical information can be taken directly from CAD data and that efficient and accurate simulations with very few unknowns can be obtained.
One fact that has hampered the widespread use of the BEM is that fundamental solutions, on which the method is based, exist only for elastic material properties and homogeneous domains. Fundamental solutions can be obtained for anisotropic materials, but they are very complicated [18]. To overcome the fact that the original BEM can only deal with homogeneous and elastic domains, several workarounds were introduced. Among them we mention the introduction of boundary element regions, to consider a piecewise heterogeneous domain [6] and the coupling of the BEM with the Finite Element Method (FEM) where the FEM regions model non-linear behaviour [25]. The topic of elastic inclusions was recently dealt with in [23]. The concept of including non-linear effects by adding a volume integral was first introduced by Brebbia [9] and Banerjee [5] and later expanded in [12]. Various ways of avoiding the generation of a volume mesh were presented (see for example [24]), but their application is limited to finite domain problems. All of the mentioned solutions involve the introduction of errors or increase the discretisation effort.
In this paper we show several innovations that make the BEM suitable for realistic simulations in underground construction. This means that the BEM simulation can consider ground support, heterogenous ground conditions and non-linear material behaviour. We first introduce the theoretical background of the BEM with volume effects. Then the evaluation of the arising boundary and volume integrals is discussed in some detail. This includes the description of the excavation geometry with NURBS patches and the definition of subdomains, where material properties differ from the ones used to compute the fundamental solutions or which behave in an inelastic way. Rock bolts are also modelled as (linear) subdomains and their analytical integration allows many of them to be used in a simulation with a small increase in the numerical effort.
Two test examples are included, which test the accuracy and efficiency of the simulation of ground support and elasto-plastic material behaviour. Finally a practical example with some complexity is presented.
It is emphasised that the simulation approach presented here does not involve any mesh generation.
Instead, geometries are defined by NURBS patches using data generated by CAD programs or input data, in a similar data format, that are user generated.

Theory
In the following we will use the word inclusion to specify part of the analysis domain that is not modelled by boundary elements, i.e. parts of the domain that has different material properties or behaves inelastically. This also applies to the ground support.
As will be explained, we use the concept of initial stresses inside inclusions to consider those volume effects, which are not considered by the boundary discretisation.
In the following we will use matrix algebra and it is therefore necessary to convert the stress and strain tensors σ ij , ij to pseudo-vectors σ, using Voigt notation: The initial stress vector σ 0 , due the fact that a point inside the inclusion has properties that are different to the ones used for computing the fundamental solutions, is given by: where is the total strain, D is the constitutive matrix for computing the fundamental solutions and D incl is the corresponding matrix for the inclusion point.
Remark: While the matrix D is restricted to an isotropic elasticity matrix, the constitutive matrix D incl can be quite general, ranging from the sparsely populated isotropic elasticity matrix to a fully populated anisotropic elasticity matrix. Indeed, later on we will use an elasto-plastic constitutive matrix for D incl .
In the following we first establish the governing integral equations and then discuss in detail how the arising volume and surface integrals are evaluated.

Governing integral equations
Consider a domain Ω with a boundary Γ, containing a subdomain Ω 0 where initial stresses σ 0 (x) are present. We apply the theorem by Betti and the collocation method to arrive at the governing integral equations. This means that we set the work done on the boundary Γ by tractions T times displacements u equal to the work done by displacements U times tractions t. We assume T and U to be fundamental solutions of the governing differential equation atx due to a source atx n and u, t to be boundary values.
If initial stresses are present, additional work is done in the domain Ω 0 by the initial stresses σ 0 (x) times the fundamental solution for strains E(x n ,x). The integral equation can be written as (see [10], [1], [4]): wherex n are the coordinates of the collocation point n. To be able to solve the integral equations they have to be regularised. The regularised integral equations are written as : where A n = 0 for finite domain problems and A n = I for infinite domain problems. The derivation of Eq. (4) and the fundamental solutions U und T are presented in [8]. The fundamental solution E is given by: where r = |x −x|, r ,i = r i r and δ ij is the Kronecker Delta. The constants are: C = 1 16πG(1−ν ), C 3 = 1 − 2ν and C 4 = 3 where G is the shear modulus and ν the Poisson's ratio.
The tensor E ijk is converted to a matrix E:

Discretisation of integral equations
To be able to solve Equations (4) we have to discretise them. This involves 2 steps: • The subdivision of the boundary domain into patches and the volume domain into inclusions • The approximation of the unknown boundary values and the approximation of initial stresses.
This will be discussed in the subsequent sections separately for the boundary and volume integrals.

Discretisation of the boundary integrals
For the numerical solution of the boundary integral equations the integrals are expressed as sum of integrals over patches: where e specifies the patch number and E is the total number of patches. In the following the geometry of patches is specified using NURBS basis functions. The advantage of this is that some geometrical shapes such as cylinder and spheres can be described exactly with few parameters. For further information on NURBS and how (7) is obtained the reader is referred to [8].
There are 3 types of patches that are useful for geomechanics modelling: finite, infinite and trimmed patches. In addition we introduce a special patch with a cutout that can be used for modelling intersections.

Geometry definition of finite patches
In Figure 1 we show an example of a finite patch. The mapping from the local ξ(ξ, η) to the global x coordinate system is given by where R i (ξ, η) are NURBS basis functions and the control points (coordinates x i ) are numbered consecutively, first in the ξ-and then in the η-direction.
The vectors tangential to the surface are given by and the unit vector normal is The Jacobian is The direction of the "outward normal" depends on how the control points are numbered.

Geometry definition of infinite patches
Here we introduce a patch definition that is useful for the simulation in geomechanics where one sometimes has to consider a surface that tends to infinity [7]. In this case we define an infinite patch as shown in Figure 2. The mapping for a patch that extends to infinity in the η-direction is given by and the special infinite basis functions are The vectors in the tangential directions are given by The unit vector normal is computed as for the finite patch. It is noted that the Jacobian J tends to infinity as η tends to 1.

Trimmed patches
Patches can be trimmed using trimming curves, resulting in more complex geometries. The trimming curves are defined in patch coordinates ξ, η. A trimmed patch is shown in Figure 3. More information about trimming can be found in [7].

Special patches
Special patches can be used to model intersections of curved surfaces with flat surfaces with a minimum of effort. A special patch is shown in Figure 4.
Its geometry is defined by: The superscript II refers to the inner and I to the outer bounding curve.
The vectors in ξ and η directions are given by: where The outward normal and the Jacobian are computed the same way as for finite patches.

Defining geometry with NURBS
CAD programs use NURBS to describe geometrical shapes (cylinder, spheres or general smooth shapes). If two shapes intersect trimming is applied. Since CAD programs are mainly designed for visualisation, the resulting intersection geometry may have small gaps. The interested reader may consult [8], where a whole chapter is devoted to the topic of how data from CAD can be used for a BEM simulation.
However, instead of asking the CAD program to compute the intersection geometry and then extract the necessary data, which is complicated, one may write a MATLAB function that computes the intersection geometry without gaps. For the practical example we have developed a function that does this.
Here we show on an example how easy it is to define complex geometrical shapes with NURBS. We start with the definition of a circular excavation and of 2 cross-passages. Only a few lines of input data (shown in Figure 5 on the left and middle) are required. The two geometries are then intersected resulting in the geometry definition on the right of Figure 5 consisting of a trimmed patch, normal patches, infinite patches and special patches. Note that no mesh generation is necessary. The dotted lines indicate integration regions (see the section on integration below).

Approximation of boundary values
To be able to solve the patch integrals in equation (7), the boundary values must be approximated.
For the approximation we also use NURBS basis functions. For normal patches the unknown boundary values are approximated byû whereR k (ξ, η) are NURBS basis functions (the hat indicating that they may be different to the ones used for describing the geometry) andû e k ,t e k are parameter values. It should be noted that in contrast to Lagrange polynomials, parameter values do not in general represent real values.
For infinite patches we have 2 choices for the displacements: • Plane strain: Displacements are constant to infinity: • Decay: Displacements decay to zero as infinity is approached (24) where K ∞ is the number of parameters andû e∞ k are the parameter values on the finite boundary.
Our refinement philosophy is to take the NURBS functions that define the geometry of the problem and refine them as necessary using knot insertion and order elevation.

Known values are defined byū
whereR k are basis functions, which may be different from the ones defining the geometry and the unknown values.
Inserting the approximations into the patch integrals allows the boundary parameters can be taken outside: where the hat and overbar has been omitted, because this depends if the values are known or unknown.
This requires only the integration of fundamental solutions times the basis functions, which will be discussed later.

Discretisation of volume integral.
The volume integral is solved numerically by dividing the volume into inclusions, defining each one geometrically. The integral is replaced by a sum of integrations over inclusions: where N i is the number of inclusions and Ω ni specifies the inclusion domain. For the numerical treatment an approximation of the initial stress σ 0 is assumed inside the inclusion and the integrals are solved numerically or analytically as will be explained later.
Two types of inclusions are considered: General inclusion that represent a volume of material that has different elastic properties or behaves inelastically and linear inclusions to represent rock bolts. It should be noted that a C 0 continuity of displacements exists at the interface between the inclusion and the domain.

Geometrical definition of general inclusion
General inclusions are defined by bounding NURBS surfaces. We establish a local coordinate system s = (s, t, r) T = [0, 1] 3 as shown in Figure 6 and map from local s coordinates to global x coordinates. The global coordinates of a point x with the local coordinates s are given by The superscript I relates to the bottom (red) surface and II to the top (green) bounding surface and x I k , x II k are control point coordinates. K I and K II represent the number of control points, R I k (s, t) and R II k (s, t) are NURBS basis functions. Note that there is a one to one mapping between the local surface coordinates ξ, η and the local coordinates s, t.
The derivatives are given by where for example: The Jacobi matrix of this mapping is and the Jacobian is J = |J|.

Geometry definition of linear inclusion
This type of inclusions is used to model cables and rock bolts. Here we assume that the geometry is defined by a linear NURBS curve and that the bar has a circular cross-section with radius R over which the stress and strain are assumed constant. The assumption is that the area of the cross-section of the inclusion is significantly smaller than that of the medium it is embedded in, allowing simplifications to be introduced for the integration. We establish a local coordinate system s = [0, 1] as shown on the right in Figure 7. The global coordinates of a point x with the local coordinate s are given by where K is the number of control points, R k (s) are NURBS basis functions and x k are control point coordinates. We also define a local coordinate system whereby the z axis is along the bar, specified by unit vector v z . The vector along the bar is given by The Jacobian is The unit vector in z direction is given by

Approximation of initial stress
For the numerical integration, introduced below, we need the values of the initial stress at Gauss points. It is inefficient to compute σ 0 at every Gauss point, whose location changes according to the location ofx n . Instead we compute the initial stress at a fixed grid of points inside the inclusion. The value of initial stress at a point with the local coordinates s = (s, t, r) T for general inclusions and s = s for linear inclusions can be obtained by interpolation between grid points: where σ 0l is the initial stress vector at grid point l with the local coordinate s l . L is the total number of inclusion points and M σ l (s) are linear or constant basis functions, which will be shown later.

Numerical integration of boundary integrals
The boundary integrals to be solved are: They are evaluated numerically, using Gauss Quadrature The integration scheme now depends on the location of the collocation point. If it is outside the patch we use regular integration otherwise we have to use singular integration.
Initially we divide the patch into integration regions depending on the following: • The location of the collocation points. They should be on integration region boundaries.
• The aspect ratios of each integration region. It should be moderate and this is particularly important for singular integration.

Regular integration
For regular integration we have to consider that the value of integrand tends to infinity as the collocation point is approached. To maintain an adequate precision of integration is crucial to the quality of the results. Therefore the number of Gauss points has to be increased near the collocation point. There is no analytical formula to determine the number of Gauss points required for a certain precision, but estimates have been worked out in [8]. The number depends on the size of the integration region and the proximity of the collocation point. The best strategy is to limit the number of Gauss points available and to subdivide the integration region into subregions. A Quadtree method that increases the number of Gauss points near the collocation point is most efficient.
Gauss Quadrature requires limits which range from -1 to +1. Therefore we introduce new local The transformation to the patch coordinate where ξ s × η s is the size of the subregion and ξ s1 , η s1 are the starting coordinates.
The integration can be written as:

Singular integration
If the collocation point is part of the subregion then the integral involving U is weakly singular. The integral is solved by subdividing the integration region into triangular subregions with the collocation point at the apex ( Figure 8). This means that the Jacobian tends to zero as the collocation point is approached. The singular integration can now be written as where N t is the number of triangles. There are now two Jacobians involved, one for the transformation from the patch coordinates to triangular coordinates (J ξ,n t (ξ i ,η j )) which tends to zero and one for the transformation from patch to global coordinates (J(ξ i ,η)).

Integration of integrals over Ω 0
For the inclusion integrals we use numerical integration for general inclusions and analytical integration for linear inclusions. For the numerical integration we subdivide the inclusion into integration regions with the same criteria as used for the boundary integration and apply Gauss quadrature. When point n is part of the integration region we have to invoke singular integration, if it is not a regular one.

General inclusions
The integral to be solved for each inclusion is : Introducing the interpolation of initial stresses (38) we have: where

Regular integration
To maintain adequate precision of integration we subdivide the integration region into subregions depending on the size of the integration region and the proximity of point n. For sub-region n s the transformation from the inclusion (s) coordinates to the coordinates used for Gauss integrationξ = (ξ,η,ζ) T = [−1, 1] 3 is given by r = ∆r n 2 (1 +ζ) + r n s where ∆s n × ∆t n × ∆r n denotes the size of the sub-region and s n , t n , r n are the edge coordinates. The Jacobian of this transformation is J n ξ = 1 8 ∆s n ∆t n ∆r n . We can write: where J(s) is the Jacobian of the mapping between s and x coordinate systems.
Applying Gauss integration we have: where N s is the number of subregions and G s , G t and G r are the number of Gauss points (which depends on the proximity of n) andξ g s ,η g t ,ζ g r the Gauss point coordinates in s, t and r directions, respectively.
W g s W g t W g r are Gauss weights.

Singular integration
If the integration region includes the point x n , then the integrand tends to infinity as the point is approached. To deal with the integration involving the weakly singular Kernel we compute the Gauss points in a local coordinate system, where the Jacobian tends to zero as the singularity point is approached.
Singular integration of general inclusions is discussed in detail in [8].

Linear inclusion, reinforcement bar
For the linear inclusions we can apply analytical integration. We model the bar as a cylindrical region with radius R and length H and assume the initial stress to be in the local z direction and to vary linearly along the bar. Since the initial stresses are computed from the strains, this means that the displacements along the bar can have a quadratic variation. We consider two types of integration: one where pointx is outside the inclusion (regular integration) and one where it is not (singular integration).

Analytical computation of regular integral
Since we assume that the cross-sectional area is significantly smaller than the surrounding medium we can assume that E is constant over the cross-section. In addition we note that the result will multiply with the initial stresses in local directions ( {σ 0 }). The integral to be solved is therefore: where the prime indicates that the result is computed in the local x , y , z coordinate system (Fig. 9). r c is the distance between the source point and a point on the axis of the inclusion .
The integral to be solved is: We choose the local axes such thatx = 0 as follows. The vector pointing in the x direction is given by: and the one in y direction is: where the capital letter indicates that the vector is not normalised.
If pointx is along the axis of the bar this computation does not work and then we assume where v y is a vector in global y-direction.
For the computation of the fundamental solution we have: and The integral to be solved is: where the linear interpolation functions are given by: The analytical solution in Voigt notation is provided in Appendix A.
Since the result of the multiplication with σ 0 has to be in global coordinates a transformation to the global system is necessary: where T is the transformation matrix given by:

Analytical computation of singular integral
Referring to Fig. 10 we subdivide the bolt into 2 subregions and obtain the following integrals in polar coordinates: if the collocation point is on the bottom of the cylinder. In such a way we have that: The terms of E = E 1 + E 2 , in Voigt notation, different from zero are provided in Appendix B.
As before a transformation to the global system is necessary:

System of equations
The discretised integral equations can be written in matrix form as: where L is an assembled left hand side, x is the vector of unknowns and r is the assembled right hand side involving known boundary values (for details of derivation see [8]). [B 0 ] is a matrix where the rows refer to the collocation points n and which multiplies with vector {σ 0 } that gathers all initial stress components at grid points inside the inclusions.

Computation of values at grid points inside the inclusion
To compute the initial stresses at grid points we need to compute the strains at these points. Even though it is possible to compute the strains directly using integral equations this is rather cumbersome as it involves the integration of strongly singular functions and involves complicated methods to isolate the singularity ( see for example [12] ). To avoid this we compute the displacements first and then compute the strains using a method that is used in the Finite Element Method.

Computation of displacements
The displacement vector u at a grid point x inside the inclusion is given by: After the solution the values u(x) and t(x) are known, so the integrals can be evaluated numerically.
where R u i (ξ, η) are the NURBS basis functions used for approximating the displacements in patch e, that contains the point x k and ξ k , η k are the local coordinates of the point. The matrix [Â] and the vector {c} have to be modified for these grid points, whereas [B 0 ] will contain zero rows in this case.

Computation of strains, general inclusions
To compute the strains we interpolate the displacements between grid points and obtain for the displacement at a point with the local coordinate s: where u n is the displacement vector at grid point n and N is the number of grid points. The interpolation functions M n (s, r, t) are obtained by considering grid coordinates s, t, r.
Rewriting (72) for a point with the local coordinates s, t, r in terms of local interpolation functions we have: where N is the total number of grid points and The derivatives of the displacements are given by Unfortunately we can not use NURBS for the interpolation functions as they are based on parameter values instead of real values. Therefore they can not be used to interpolate the real displacement values at internal points. We use Lagrange polynomials instead.
The strains are given by: The strains at grid point k can be written in matrix notation: and The global derivatives of M n are given by: where J is the Jacobian matrix Eq. (33). For a linear inclusion we compute the strain in local directions as is shown later.
Gathering all strain vectors at grid points in { } we can write: After substitution of {u} we obtain: The initial stresses are computed by where [D − D incl ] is a matrix containing D − D incl as sub-matrices on the diagonal.

Computation of strain for linear inclusions
For linear inclusions it is convenient to work with the strain in local coordinates. If we assume the bolt to be fully grouted, i.e. no slip is allowed between the bolt and the domain it is embedded in and that the Poisson's ratio of the bolt has no effect, the only strain that has to be considered is the one along the bar 1 : where J is the Jacobian and v z is a unit vector in z direction.

Eq. (91) now becomes
The local initial stress vector is given by: where E and E incl is the Young' modulus of the domain and the inclusion respectively and

Solution procedure
A solution that already includes the effect of inclusions that have different elastic properties is possible by combining equation (68) with (94). Eq. (94) can be written in the following form: where: Eq. (101) along with Eq. (68) form the following linear system of equations: that can be solved in terms of boundary unknowns and internal strains.
It is also possible to obtain a system of equations that only multiplies with the boundary unknown: where [L] and {r} are modified left and right hand sides.
We rewrite the strain vector as: where Eq. (105) can be inserted in Eq. (68) in order to obtain: and, hence, the following system of equations can be obtained The matrices in Eq. (104) are defined by:

Elasto-plasticity, Newton-Raphson method
If the stress at an inclusion point exceeds the elastic limit incremental/iterative elasto-plastic procedures, well known in the FEM, are applied. A detailed discussion of these methods is beyond the scope of this paper. A very good description can be found in [22]. The increment of stress in the plastic regime is: where * is the in-elastic strain increment (i.e. the one that occurs after the stress has reached a state where F(σ)=0) . The in-elastic ( * ) and elastic ( e ) strain increments can be computed by where is the total plastic strain increment and F new is the value of F at the end of the increment, F old is the value at the beginning.
The initial stress due to plasticity is given by: If the increment in total strain has occurred while traversing the yield surface then: The total initial stress increment, including plastic effects, is given by where If the strain increment is totally elastic (f=0) we have: If the strain increment is totally plastic (f=1) we have: To start the simulation we set f in Eq. (117) equal to zero i.e. D = D e,incl and obtain the first result: With this result we compute the value of the yield function at internal points and compute a load factor λ that reduces the load to the one where first yield occurred. We then reduce the results to first yield: We also adjust the stresses at internal points to the new load level. For the first load step (n=1) the matrix D is updated with f = 1 for the point where the first yield has occurred.
We apply the rest of the loading in n steps. The increment in load is: At each increment we solve: For the second and subsequent load steps we update the matrix D with f according to the current state of stress. Standard return alogrithms can be applied to ensure that the stresses stay on the yield surface.
Remark: Note that incremental/iterative procedures need only be applied for elasto-plastic behaviour.
When inclusions are defined that have different elastic behaviour the solution is obtained without iteration.

Test Examples
We test the implementation on an example of a circular, infinitely long, tunnel in an infinite domain subjected to a virgin stress (see Figure 11). Since the aim of the examples is to test the accuracy of the implementation and not to perform a realistic simulation, we use non-dimensional parameters. For the domain we assume E=1 and ν = 0 and for the radius of the tunnel, R=1. The results of the IGABEM simulation are compared with a FEM analysis using the software PLAXIS and an analytical calculation.

Elastic analysis with no ground support
In order to establish a comparable regime between the IGABEM and PLAXIS simulations we first conduct an elastic simulation without ground support and with a virgin stress σ x = 0, σ y = 0, σ z = −1.

Discretisation with PLAXIS
The infinite domain is approximated by providing an artificial boundary at some distance from the tunnel. Since this distance affects the results we examine the error introduced by varying the distances from the tunnel centre to the artificial boundary to 2.5, 5 and 10 times the tunnel diameter (2.5D, x y z Figure 11: Geometry of test example. 5D, 10D). The boundary conditions at the outer boundary of the mesh (shown in Figure 12) are that displacements normal to the faces of the cuboid are set to zero. The fact that the tunnel is infinitely long can be modelled by extending the mesh along the tunnel axis to 10m (=5D). The mesh consists of quadratic 10-noded tetrahedral elements. For load case 0 the virgin stresses are assigned to all elements of the mesh and the excavation is then simulated by deactivating the elements inside the tunnel. This will give the same result as a much simpler plane strain simulation but the same mesh will be used for the second test example where a 3-D analysis is necessary.

Discretisation with IGABEM
Only the boundary of the tunnel is discretised as the infinite domain is explicitly considered by the The discretisation is shown in Figure 13. For the approximation of the displacements the same basis functions as for the description of the geometry are used, resulting in the collocation points shown. The discretisation has 48 degrees of freedom.

Comparison of results
To investigate the effect of the artificial boundary on the results of the FEM analysis we compare the values of vertical displacements along a vertical line above the tunnel with the exact result (Kirsch solution [15]) and the IGABEM result. It can be seen in Figure 14 that the FEM results converge to the exact results, while the IGABEM result already is in excellent agreement. The computing time for the PLAXIS simulation was 3 minutes 53 seconds and for the IGABEM simulation 8 seconds.
x y z Figure 15: Geometry of the elastic simulation with bolts Figure 16: Bolts simulated as embedded beams with PLAXIS

Elastic simulation with rock bolts
This example is designed to test the implementation of rock bolts. Three rock bolts are installed at the top of the tunnel (Figure 15). The diameter of the bolts is 0.05 and the elastic modulus is twice the value of the domain (i.e. E bolt = 2). It is noted that the consideration of the rock bolts renders the analysis three dimensional. The aim is to study the local influence of the rock bolts, assuming that further away plane strain conditions prevail.

Discretisation with PLAXIS
In PLAXIS rock bolts are simulated as embedded beams. The geometry of the bolt can be arbitrarily located in the finite element mesh, thus nodes of the bolt do not need to coincide with nodes of the finite element mesh. As the beam cross section is small, compared with its length, the bending stiffness is almost zero, therefore it can be assumed that the beam acts as a bolt. In Figure 16 the three bolts are shown embedded in the surrounding FEM mesh.
x y z Figure 17: Discretisation of the problem with IGABEM. Shown are the bolts in blue with the associated control points as hollow squares. The collocation points obtained after knot insertion are shown as red filled squares.

Discretsation with IGABEM
The bolts are defined as linear inclusions. In order to get results comparable to PLAXIS the continuity of displacements is changed to C 0 at the point where the rock bolts meet the boundary patch. This is done be inserting knots into the basis functions that approximate the unknown. This increases the number of collocation points and the degrees of freedom. The IGABEM discretisation is shown in Figure 17 and has 108 degrees of freedom.

Comparison of results
We compare the displacements along the rock bolts in Figure 18. Good agreement can be observed.
The computation time of the PLAXIS simulation was 6 minutes and 12 seconds and for the IGABEM simulation 24 seconds.
x y z Figure 19: Volume discretisation for the plasticity example. Shown are the two (green and red) surfaces defining the inclusion and the control points as hollow squares. Also shown are the internal points as stars and the subdivision into integration regions for the volume integration. These regions are automatically determined depending on the location of internal points and on the aspect ratio.

Elasto-plastic simulation without rock bolts
For this test we can compare with an analytical solution in plane strain due to Duncan-Fama (see for example [13]). For a hydrostatic virgin stress of p 0 (compression positive) and a Mohr-Coulomb yield condition with a friction angle φ and cohesion c, the solution for the elasto-plastic radial displacement of the tunnel wall is given by: where the extent of the plastic zone is given by: and

Discretisation with IGABEM
For this simulation the excavation is surrounded with a general inclusion extending one radius from the excavation surface as shown in Figure 19. The properties assigned to the inclusion are E = 1, ν = 0, φ = 10 • , c = 0.5. A hydrostatic compressive virgin stress p 0 = 1 was applied.

Comparison of results
For the above input data the analytical solution for the extent of the plastic zone was computed as 1.3, which means it is within the inclusion. The theoretical solution for the radial displacement was 1.262 which compares well with the IGABEM solution of 1.269. Convergence to 1% of residual was achieved after 6 iterations.

Practical example
The  The aim of the simulation is to investigate if a pre-installation of ground support would lead to a safer and more economic excavation. The required excavation stages are shown on a 60 m section of the cavern in Figure 20. In the first stage a tunnel system consisting of circular tunnels (by raise boring) and rectangular horizontal tunnels (by blasting) is excavated. Next cables are installed from the circular tunnels of such length that they reach to the surface of the cavern to be excavated. Finally the cavern is excavated. The cross-sectional diameter of the cables is 25 mm with the modulus of elasticity of E=210 GPa.

The simulation model
For a preliminary analysis we analyse a 60m section of the cavern At the edges of the section we use infinite plane strain boundary elements. Nine cables per circular excavation are considered. We start with the definition of the excavation surfaces using the method outlined in section 3.1.5. For the description of the smooth shape of the cavern walls 4 control points and a basis function of order 3 has been used for each half. For the variation of the unknowns we use the same basis functions as for the description of the geometry except that for the bottom surface the basis function was elevated by one order from linear to quadratic across the cavern. Next the inclusions are specified. This relates the cables and to the part of the domain where it is assumed that non-linear behaviour is taking place. In this study we concentrate on the rock mass behaviour between the tunnels and the excavation surface, so this part is selected. The definition of the inclusions is shown in Figure 22 and the final model in Figure 23.

Preliminary results
The novel approach to simulation will be used for a detailed study to determine if pre-installation of ground support can make the excavation of large underground caverns more feasible, safe and economic.
Here only a preliminary result is presented in Figure 24. For this case the convergence to 1% residual was achieved in 10 iterations. It is clear that more cables are required to provide adequate pre-support but it should be noted that the number of cables can be increased substantially without significantly increasing the simulation effort. This is because cables are generated automatically and because the numerical effort only increases slightly due to the fact that analytical integration is used. It should be noted that the number of degrees of freedom is not increased by increasing the number of cables.

Summary and Conclusions
A new approach to simulation for underground excavations using isogeometric methods and NURBS has been presented, which makes it extremely user friendly and efficient. No mesh generation is required and geometric data can be specified in a CAD data format. Since CAD software also uses NURBS and trimmed NURBS surfaces the connection to CAD is natural. Parts of the domain that have different material properties or behave in an inelastic way can be considered but a definition of inclusions is necessary. This means that for most practical applications only near field effects are considered in the simulation, i.e inclusions have limited extent and are located near the region of interest. If this restriction can be accepted, then the software can do everything a FEM simulation can do, but with a drastic reduction of user effort and execution time.
One of the main contributions of the paper is the presentation of analytically integrated rock bolt inclusions. Numerical integration is the most compute intensive aspect of any BEM simulation. In case of large patches the number of Gauss points can be high. This is because of the nature of the integrals, which exhibit singular behaviour, i.e. where the value of the integrand increases rapidly as the source point is approached. For the numerical volume integration, where the integration is in 3 directions the computational effort is particularly high. However, it is noted that the matrices that require volume integration only depend on geometry and can be precomputed for a particular model. Also, the computations are able to exploit to a very high degree the use of multiple processor hardware.
At the time of writing the simulation model has been written in MATLAB, had a first application to a real problem and will be used for a detailed study to simulate the effect of pre-support for large caverns.
Although concrete arches and shotcrete can be modelled with the general inclusion approach presented here, it becomes cumbersome, when the thickness is small. Therefore, special shotcrete inclusions with a semi-analytical integration are being developed next. It should be pointed out that it is also possible to develop inclusion geometries, that are not restricted to the cuboid shape presented here (in the sense that only two surfaces can be of general shape).
Finally, it is hoped that this contribution will make the BEM more attractive for simulations of underground excavations.

Acknowledgements
The authors would like to acknowledge the significant contribution of Scott Sloan to the topic of limit analysis. The first author had the pleasure of spending some time working at the centre of excellence at the University of Newcastle directed by him and was impressed by his leadership and quest for innovation.