Analytical solutions for rheological processes around bores and tunnels

Five analytical solutions are presented for processes of linear viscoelastic, homogeneous, and isotropic solids around freshly opened bores / tunnels in various initial stress fields. The solutions are obtained via a simple and direct realization of Volterra’s principle. This realization is based on an appropriate decomposition of the known solution of the corresponding elastic problem and leads to ordinary differential equations in the time variable in the viscoelastic case. Fairly rich temporal behaviors are revealed.

fitted based on experimental data, as is done, for example, in the Anelastic Strain Recovery method of determining 3D in situ stress [16][17][18], where all the Kluitenberg-Verhás coefficients have been found to be needed, both in the deviatoric part of stress and strain and in the spherical part.
Various other fields of engineering also find such rheological models relevant, for solid materials like plastics, asphalt, biomaterials, etc., Pressurizing of thick-walled plastic tubes (an example that also turns out to be related to the present paper, as revealed at the end of the Conclusions) is just one example among the many, and the damping and delaying properties of such materials are not necessarily disadvantageous but can also be technologically benefitted, as utilized, for instance, for absorbing vibration [19]. All these motivate the development of solution methods for rheology and, as part of it, for linear viscoelasticity.
A known general treatment of linear viscoelastic problems is based on integral equations [20][21][22][23]. Within this framework, Volterra's principle [24] has also been realized. Volterra's principle says, in brief, that the solution of a linear viscoelastic problem can be derived from the solution of the corresponding elastic one. This integral equation-based approach has been utilized to find various numerical solutions. Obtaining analytical solutions is far less straightforward since an operator inverse is to be calculated analytically.
In order to achieve analytical solutions, the strong heuristic power of Volterra's principle has recently been revealed to have simpler and more direct manifestation in certain examples [25]. In these cases, the viscoelastic solution has been obtained by replacing the elasticity coefficients of the corresponding elastic solution by timedependent functions and by determining these time dependencies from the emerging ordinary differential equations. In one of the treated examples, the result has proved to reproduce the known solution gained by direct means in [26] (while the other cases appear to be new results).
The realizations presented in [25] have, however, enabled to treat only a limited range of problems. Here, we present a considerable generalization that-while still not being as general as the integral equation-based formulation-is shown here to cover numerous practically interesting situations.
The approach introduced here can be summarized as follows: 1. we decompose the known stress and strain solutions of the elastic problem corresponding to the fully opened bore/tunnel according to different dependencies on the elastic coefficients, 1 2. we model the gradual opening of the bore/tunnel via time-dependent rescaling of the boundary conditions and rescale the components of the elastic solution accordingly, 3. we assume the solution of the viscoelastic problem as a time-dependent linear combination of the identified elastic components and solve the emerging set of linear ordinary differential equations, 4. the displacement field is derived from the strain solution via Cesàro's formula.
In Sect. 2 of the paper, the framework is specified, including the basic equations and the aspect of boundary conditions. Section 3 describes the method in general terms. Then, in Sect. 4, five applications are provided, with figures that demonstrate the characteristics of the obtained solutions. Section 5 closes the paper with conclusions and an outlook.

Elasticity, viscoelasticity, and the class of problems to treat
We consider purely mechanical problems of homogeneous and isotropic continuous media, and our aim is to determine the displacement field u, the strain field ε and the stress field σ (where both ε and σ are symmetric tensors). We work in the force equilibrial approximation, i.e., when acceleration is neglected: 2 σ · ← ∇ = − g, (1) where is mass density and g is volumetric force density (assumed to be time independent), and ← ∇ and → ∇ are the nabla operators acting to the left and to the right, respectively (reflecting proper tensorial order 3 ).
Concerning ε, we stay in the small-strain approximation, which then imposes the geometric compatibility equation in the form According to mathematics, to a symmetric tensor field ε with property (2), there exists a vector field u Cauchy -called hereafter Cauchy vector potential-from which ε can be obtained as The Cauchy vector potential is not unique for a given ε, and all Cauchy vector potentials can be derived from the strain field according to Cesàro's formula (see, e.g., [27]), with A 1,3 denoting antisymmetrization in the first and third tensorial indices, where each of the position vector r 0 , the path of integration, the vector function u 0 (t), and the antisymmetric tensor function (t) is arbitrary. The displacement field u is one of these Cauchy vector potentials so when we wish to reconstruct u from the strain field then we need to fix these uncertainties (this rigid-body motion freedom) using symmetry arguments and other physically plausible considerations. In case of linear elasticity (for a homogeneous and isotropic medium), connection between stress and strain is provided by Hooke's law, in the deviatoric-spherical separation, where σ s = 1 3 (tr σ ) 1 denotes the spherical part-which is proportional to the identity tensor 1-, while σ d = σ −σ s is the deviatoric (traceless) part; furthermore, E d = 2G is the deviatoric elasticity coefficient and E s = 3K is the spherical one (G being the shear modulus and K the bulk modulus).
For linear viscoelastic models of solids, one can generalize Hooke's law by replacing the elasticity coefficients with polynomials of the time derivative operator. Namely, where the stress-related operators S d , S s , and the strain related ones E d , E s are Each time derivative term brings in a time scale (via its coefficient). 3 Correspondingly, in indexed notation with Cartesian tensorial components, (1) reads ∂ j σ i j = − g i , and (3) is For an illustration of what viscoelastic material models are covered by (6)- (8), see Appendix A, where it is shown that models given in a Prony series form can all be re-expressed as (6)- (8).
In our applications, we concentrate on the Kluitenberg-Verhás model family [15,28] hereafter, overdot abbreviates partial time derivative. (Note that the small-strain assumption allows to approximate the substantial time derivative with the partial time derivative.) In the case of elasticity, equations (1), (2), and (5) form the system of equations to be solved. Together with appropriate boundary conditions imposed at the boundaries of the spatial domain considered, the solution exists and is unique; however, to obtain this solution is not necessarily simple, since (5) poses separate conditions for the deviatoric and spherical parts, while (1), (2), and the boundary conditions prescribe the requirements for the sum of the deviatoric and spherical parts.
When one deals with the above-described viscoelastic generalization of the problem then, in addition to (1), (2), (6) and the boundary conditions (which may be time dependent in general), initial conditions are also required to ensure uniqueness of the solution, since the constitutive equations contain time derivatives [see (7) and (8)]. In the viscoelastic case, all fields are functions of both time and space, and all equations and boundary conditions have to be satisfied for all time instants, which raises an even more complicated task than for the elastic counterpart. Hence, inspired by Volterra's above-mentioned idea, we wish to utilize the fact that the elastic solution satisfies all the spatial requirements so, by relying on the structure of the elastic solution, we reduce the problem to solving a set of temporal differential equations.
As the first step, we separate the effect of the force density by subtracting some such time-independent fieldsσ andε-henceforth: primary fields-that satisfy the equations [note that, for time-independent stress and strain fields, (6) gets simplified to Hooke's law with E d = E d 0 , E s = E s 0 , cf. (7)- (8)]. Then the difference fields-henceforth: complementary fields- satisfy the homogeneous equationŝ Naturally, in general, this transformation modifies the initial and boundary conditions, which are to be taken into account during the calculations. If the spatial domain filled with the medium has more than one boundary (more than one connected boundary surface) then the problem can be divided into subproblems in which only one boundary condition is nonzero (the boundary condition is nonzero only on one connected boundary surface). Henceforth, we always analyze one such subproblem. Thanks to linearity of all equations involved, the sum of such subsolutions yields solution for a whole problem.
In order to provide models for tunnels, we consider problems with boundary conditions that are specified for stress rather than displacement. More closely, the normal component of stress is to be prescribed as the function of time. This time dependence will be allowed with the limitation that it must mean a space-independent rescaling λ(t) of the boundary condition-like gradual loading/unloading of a surface where loading may be space dependent but the ratio of normal stress values at two different boundary points is time independent. In other words, the boundary condition must realize a homogeneous amplification/tuning along the boundary. The time-dependent multiplier λ(t) can be quite arbitrary, the only restriction being that it be sufficiently many times differentiable. For gradual switching on, like when modeling drilling, it can be chosen as Corresponding to such a time-dependent homogeneous rescaling of the boundary condition, the solution of the elastic problem also gets rescaled-space independently rescaled-by the factor λ(t) .
Another remarkable property of the elastic solution (at any fixed t) is that, based on dimensional reasoning, dependence of the stress solution on E d and E s must be such that stress only depends on the dimensionless ratio In other words, it depends only on the Poisson's ratio ν, which is related to η as This property is, naturally, directly visible in the examples considered below. Accordingly, it proves beneficial to use, instead of strain, a multiple of it that has the dimension of stress. This can be simply achieved by called hereafter stress-dimensioned strain. (Note that, for solids, E s is always positive.) Correspondingly, The compatibility equation (15) remains in the same form forζ : Withζ , Hooke's law is simplified tô Then, it is apparent that, in a solution of an elastic problem with stress boundary condition,ζ also depends on the elasticity coefficients through η only [see (18)]. It is to be emphasized that, while customarily one only focuses on the space dependence of an elastic solution, for the solution method described here, dependence on the elasticity coefficients E d , E s -more closely, on η-will also be of central importance.
In parallel, the viscoelastic problem characterized by (16) imposes between stress and stress-dimensioned strain, where [cf. (7)- (8)] (as has been mentioned earlier, E s 0 of the viscoelastic case plays the role of E s of the Hooke case).

The method of four elastic spatial pattern sets
For a fixed η, the elastic solution satisfies the spatial requirements (14), (22), and the boundary condition (at a given instant). This holds for any allowable η so the idea is to seek for the viscoelastic solution as a linear combination of elastic solutions with various η's (the sum of their weights being 1, to fulfill the boundary condition). If the initial conditions can also be given as such a linear combination-and this is the case whenever the initial state of the continuum is a static, therefore, elastic one 4 -and if the time-dependent weights respect the temporal condition (6) then we have obtained the solution of the viscoelastic problem. Although continuously many allowable ηs exist, the corresponding elastic solutions may not be linearly independent. In order to explore the range of freedom we have, we decompose the elastic solution (for the boundary condition at some given instant): where I , J K, and L are some integers-assumed finite throughout this paper-, a i (η)'s are linearly independent coefficient functions, α i (r)'s are linearly independent space dependencies-called hereafter 'spatial patterns'-, b j (η)'s are also a set of linearly independent coefficient functions [not necessarily the same ones as the a i (η)'s], β j (r)'s are also linearly independent spatial patterns [not necessarily the same ones as the α i (r)'s], and the c k (η)'s, d l (η)'s, γ k (r)'s, δ l (r)'s are analogous function sets. These decompositions are different from expansions (Fourier etc.) where the expansion happens with respect some standard predefined function set. Here, the spatial patterns are specific to the given problem. Nevertheless, in practice, the decompositions may be straightforward to perform (as is the case in each of the examples considered below, as we will see): one needs to identify linearly independent η dependencies, and the corresponding r dependencies-the spatial patterns-follow uniquely.
From (27), the expansionŝ follow due to Hooke's law (23). Further relationships among the expansions are less direct, e.g., in all the examples presented in Sect. 4, the spherical parts of the α i (r)'s are already linearly dependent and, hence, are ineligible as the independent spatial patterns δ l (r) inζ s el andσ s el . Nevertheless, in all the examples, the decompositions can be readily performed. In parallel, the general conditions K ≤ I, L ≤ J, I ≤ K + L , J ≤ K + L are also easy to find: the first two follow from that taking the deviatoric (or spherical) part of linearly independent tensor-valued functions may decrease but cannot increase the number of linearly independent tensor-valued functions (accidental coalescences may happen), while the two remaining inequalities are due to that adding K terms and L terms may produce at most K + L linearly independent terms.
When the boundary condition is multiplied by the time-dependent factor λ(t) , the elastic solution gets rescaled accordingly: Next, the solution of the viscoelastic problem σ rheol (t, r),ζ rheol (t, r) is sought for in a form where the ηdependent factors of the spatial patterns are replaced by time-dependent factors. Now we see that, in spite of the continuously many possible η's, the range of freedom is only finite dimensional: over the I spatial patterns α i (r) themselves is that the former ones satisfy the spatial condition (14) and, analogously, thê ζ el (η n , r)'s respect (22). Accordingly, our ansatẑ for the viscoelastic problem obeys both (14) and (22). The boundary condition imposes and the only remaining requirement on the unknown coefficient functions ϕ m (t) , ψ n (t) is (24). For discussing (24), we need the relationship between the deviatoric part of the expansion ofσ el (η, r) given in (26) and that ofσ d el (η, r) seen in (28): the α d i (r)'s may not be linearly independent but the γ k (r)'s are, hence, any α d i (r) must be expressible as a linear combination of the γ k (r)'s. Analogous statements apply to the β d j (r)'s, α s i (r)'s, and β s j (r)'s so altogether hold with appropriate matrices A ik , B jk , C il , and D jl . Substituting these into the deviatoric and spherical parts of (30) yieldŝ The viscoelastic conditions (24) generate the system of equations Since both the γ (r)'s and δ(r)'s are linearly independent sets, the equality of the corresponding coefficients follows: This is the set of temporal ordinary differential equations one has to solve for the ϕ m (t) 's and ψ n (t) 's. Altogether, we have K + L + 1 equations for the I + J unknowns (the +1 is (32)). In all the examples considered in Sect. 4, K + L + 1 = I + J so the viscoelastic problem can be solved uniquely.

Examples
Five concrete examples follow, with figures that demonstrate the typical behavior revealed. We use for the function λ(t) that switches the boundary condition on [cf. (17)]. Three practically interesting different cases are considered: • when the switch-on λ(t) is very slow-i.e., the switch-on time scale t 2 −t 1 is small-compared to the viscoelastic time scales, • when the switch-on time scale is comparable to the viscoelastic time scales, • when the switch-on λ(t) is very fast compared to the viscoelastic time scales.
In the figures below, the functions d denote the time-dependent factors of the spatial patterns γ k (r) and δ l (r), read off from (39)-(40). As already indicated, while there is arbitrariness in the choice of the η k 's, the spatial patterns themselves have been fixed uniquely at the level of the elastic solution so the coefficients of the patterns do not contain any arbitrariness. The displacement field-understood, naturally, with respect to the primary initial state of the continuum-is, via (4), where the rigid-body like displacement-and-rotation freedom is fixed by prescribing u → 0 for asymptotically distant points in the first two examples and, in the further three examples, by the heuristically analogous condition that the highest points of the ground level have no vertical displacement. For the plots, two concrete viscoelastic models are considered. Since, according to experience, viscoelastic material behavior is more prevalent in the deviatoric part, we take Hooke's elasticity model for the spherical part. For the deviatoric part, our two choices are the Kelvin model and the Kluitenberg-Verhás model. Remarkably, already the deviatoric Kelvin model leads to a Poynting-Thomson-Zener behavior in uniaxial loadings, exhibiting creep and stress relaxation.
In case of the Kelvin model in the deviatoric part and Hooke model spherically, the viscoelastic time scale of the model is Using this as the unit of time, in the plotted examples, slow switch-on of the boundary condition is represented by t 2 − t 1 = 5, the medium case by t 2 − t 1 = 1, and fast switch-on by t 2 − t 1 = 0.1. The plots are calculated with η = E d /E s = E d 0 /E s = 0.4 [to which the Poisson's ratio ν = 0.25 corresponds, cf. (19)].
When taking the Kluitenberg-Verhás model in the deviatoric part, for simplicity, we choose zero for the coefficient ofσ d -in this case the index of inertia (see [15]) is necessarily positive, so (damped) viscoelastic oscillation will be present-: Notably, models with positive index of inertia may be relevant for materials with remarkable micro-or mesoscopic inertia (similarly to micropolar continua).
Since the bore/tunnel reaches its fully open final state within finite time, the viscoelastic solution is expected to asymptotically tend to the elastic counterpart. For the coefficient functions (44)-(45) (the time-dependent factors of γ k (r) and δ l (r), respectively), this tells, in view of (28) resp. (27) In each of the examples below, we have performed this consistency check and found agreement. The details are presented for the first example below.

Cylindrical tunnel opened in infinite, homogeneous, and isotropic stress field
In an infinite, homogeneous, and isotropic stress fieldσ , we open an infinite cylindrical tunnel with radius R (see Fig. 1). In cylindrical coordinates, the boundary conditions specifying the elastic solution for the completely open bore are which are rewritten for the complementary field aŝ The solution of the elastic problem, for this completely opened bore, is well known (see, e.g., [4, p. 155] or [5,6]): Comparison with the general formulae (26) and (34) finds Then, substituting into (41), (42), and (32), one obtains the following system of differential equations for the viscoelastic problem: two equations for the two unknowns ϕ(t), ψ(t).
For this example, the outcome [i.e., (57)-(58)] can be checked against an already known result. Namely, the same problem has already been solved in [26], there via a different approach. As one can check, the present outcome is in agreement with the result found in [26] (i.e., Eq. (56) of [26]): the obtained differential equations are the same, with the same initial conditions. 5 Figure 1 illustrates the solution.
The consistency check (49) says for this example (and there are no nonzero coefficient functions s l , s l ). Reassuringly, Fig. 1, which has been prepared with η = 0.4, 1 η = 2.5, is consistent with these asymptotic values 1 and 2.5, respectively. As a further illustration, Appendix 1 presents, for this example, a comparison of the analytically determined coefficient functions to numerically calculated ones (and also finds good agreement).

Cylindrical tunnel opened in an infinite and homogeneous but anisotropic stress field
Next, we generalize the previous example by allowing the initial stress field to be anisotropic.
Hence, the expansion ofσ el (η, r) in (26) holds with Taking the deviatoric part yields the expansion ofσ d el (η, r) with Since α s 1 (r) = α s 2 (r), the expansion ofσ s el (η, r) only contains one term: The matrices in the linear relationships (34) can be identified as All these together specify the two deviatoric viscoelastic equations and the spherical equation These are supplemented by the requirement from the boundary condition, See Figs. 2 and 3 for typical solutions.

Cylindrical tunnel opened in a homogeneous medium under the action of gravity
The remaining three examples discuss the viscoelastic process of a semi-infinite gravitating domain (also called "dead load") weakened by a bore/tunnel, for three lateral pressure coefficients. The elastic solution of the problem can be found in [1][2][3]. 7 The primary stress field, in an appropriate Cartesian coordinate system (see the outline in  Fig. 2 Outline and typical results for a cylindrical bore (tunnel) opened in an infinite and homogeneous but anisotropic stress field in a Kelvin-Hooke material. Upper right: (Enlarged) time evolution tendency of the displacement field, for a fast opening (blue: original boundaries; green, yellow, orange, and red: later snapshots of the boundaries). Such distortions are not uncommon for tunnels, see [10][11][12][13]. Second and third rows: Time evolution of the coefficient functions for slow (left column), medium-speed (middle column), and fast (right column) openings (indices: d1 black, d2 green, s1 red). In addition to the anisotropic distortion of the contours, note the nontrivial time dependencies depicted by the greenish lines in the middle and right columns: they start with opposite sign with respect to the later and asymptotically taken sign. This leads to temporary counter-intuitive motion components. (Color figure online) Fig. 4), can be written as where γ = g describes the homogeneous force density, d is the depth of the center of the bore from the surface, and parameter k is the lateral pressure coefficient. [1] gives the solution of the problem for three different values of k: • When one considers a hydrostatic pressure distribution for the primary field then k = 1.
• When the dilatation of the medium is laterally inhibited then the strain componentsε x x andε zz are zeros so one can derive k = ν 1−ν = 1−η 1+2η . • When the dilatation of the medium is laterally free then the stress componentsσ x x andσ zz are zeros so k = 0. Transforming (73) to a cylindrical coordinate system yields the stress components The boundary conditions are prescribed for the contour of the cylinder and for the plane surface-the horizontal boundary-after the drilling; on these boundaries, the normal component of stress is zero.
Mindlin gives the solution in form of an infinite series in the bipolar coordinate system that suits to both the cylinder and the plane [1]. If the ratio of the depth d of the center of the bore and the radius R of the cylinder satisfies d/R > 1.5-large-depth approximation-then it suffices to take the leading order term from the infinite series.
This first term is transformed to cylindrical coordinate system in [2,3]: The entries of the complementary field (having a plane strain situation) arê The three above-mentioned special cases k = 1, k = ν 1−ν = 1−η 1+2η , k = 0 are discussed in the following three subsections.

Gravitating solid, hydrostatic initial stress state (k = 1)
In this case, the complementary field can be written aŝ with so there are two expansion coefficients, The deviatoric part ofσ el (η, r) delivers the expansion ofσ d el (η, r), with In parallel, the spherical parts merge again, with Correspondingly,ζ el (η, r) iŝ telling Reading off the matrices yields The resulting deviatoric viscoelastic equations are while the spherical equation is Finally, the boundary condition imposes Typical outcomes are visible in Fig. 4. The corresponding matrices are The resulting viscoelastic equations read while from the boundary condition, we have Figures 5 and 6 illustrate the solutions.

Gravitating solid, free lateral deformations (k = 0)
In this case, the complementary stress field consists of three independent spatial patterns again: with ). Upper right: (Enlarged) time evolution tendency of the displacement field, for a fast opening (blue: original boundaries; green, yellow, orange, and red: later snapshots of the boundaries). Second and third rows: Time evolution of the coefficient functions for slow (left column), medium-speed (middle column), and fast (right column) openings (indices: d1 black, d2 green, d3 orange, s1 red, s2 blue). The more elaborate initial condition yields a richer structure of spatial patterns, some of which start with opposite sign compared to the asymptotic equilibrial one. (Color figure online) Taking the deviatoric part yields the expansion ofσ d el (η, r) with In parallel, the spherical part gives one less spatial pattern again, as follows from α s 3 (r) = 2 α s 1 (r) − α s 2 (r) : The resulting expansion ofζ el (η, r) has The matrices prove to be With these, the deviatoric viscoelastic equations are and the spherical ones are The boundary condition requires to fulfill Solutions are displayed in Fig. 7.

Conclusions
Viscoelastic properties of solids can lead, around a freshly opened bore/tunnel, not only to delayed and damped response but-even in the simplest 'deviatoric Kelvin-spherical Hooke' case-to direction dependent and temporarily counter-intuitive motion, as the above figures illustrate. Analytical solutions provide means to easily explore such and other related features, qualitatively and, depending on situation, quantitatively as well. When some of the viscoelastic coefficients are unknown or only vaguely known, the explicit formulae allow for fast evaluation for many trial coefficient values. When a numerical-primarily, finite-element-solution is also available, it can be validated using the analytical solution. This is actually one of our plans for the future: to realize the problems presented here numerically as well, and after satisfactory agreement, the numerical version is safe to generalize to settings where the geometry or some other aspect make the analytical approach unfeasible. Here, as a first step, we have provided numericalfinite-difference-comparison to the analytical result obtained for one of the example problems (cf. Sect. 4.1 and Appendix 1). For problems with less geometric symmetry, a finite-element approach appears more feasible. This, however, is not straightforward as finite-element softwares prefer to consider viscoelastic problems with boundary conditions on displacement (see, e.g., the manual [30] of the software Abaqus, which realizes approaches taken from [31,32], among others) while our analytical results are for boundary conditions on stress. Another nontrivial issue regarding comparison of analytical and numerical results is numerical error: simulating complex geometries for many time steps at small numerical error (to make the comparison to the analytical prediction meaningful) in a stress-oriented-rather than displacement-oriented-finite-element realization requires very many discrete degrees of freedom, rendering the simulation severely resource-intensive. Some suitable efficient approach is necessary, and this is the direction we plan to pursue in future.

Fig. 7
Outline and typical results for a cylindrical bore (tunnel) opened in a homogeneous Kelvin-Hooke medium under gravity, with free lateral deformations in the primary field (k = 0). Upper right: (Enlarged) time evolution tendency of the displacement field, for a fast opening (blue: original boundaries; green, yellow, orange, and red: later snapshots of the boundaries). Second and third rows: Time evolution of the coefficient functions for slow (left column), medium-speed (middle column), and fast (right column) openings (indices: d1 black, d2 green, d3 orange, s1 red, s2 blue). A remarkable difference from the previous two examples is that, here, the contour of the tunnel expands in the horizontal direction. (Color figure online) In parallel, further investigation would be beneficial for the analytical solution method (presented here) itself, exploring the possibilities of generalization, for example, in the direction of infinitely many expansion terms (see e.g., [1,9]) with a hierarchy so any finite truncation of the expansion would lead to informative approximation of the full result. Here, we have initiated this truncation methodology by creating the solution of the gravity-related problems at the level of the first term of such an expansion. Such a hierarchic approach is not obvious to realize in the general integral equation formalism.
At last, it is worth repeating that the approach and solutions presented here (and the solutions for thick-walled tubes and spherical tanks and hollows in [25] and in the preprint [33]) find applications in other fields of engineering and technology as well, for describing plastics, asphalt, biomaterials, composites,and structured materials. etc..

Funding Information
Open access funding provided by Budapest University of Technology and Economics.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Fig. 8
Finite-difference numerical results compared to the analytical results, for the cylindrical tunnel opened in infinite, homogeneous, and isotropic stress field. Up to numerical error, a good agreement can be observed