Development of TD-BEM Formulation for Dynamic Analysis for Twin-Parallel Circular Tunnels in an Elastic Semi-Infinite Medium

Weidong Lei,Hai Zhou,＊,Hongjun Li and Rui Chen

1School of Civil and Environmental Engineering,Harbin Institute of Technology,Shenzhen,518055,China

2Urban and Rural Construction Institute,Hebei Agricultural University,Baoding,071001,China

ABSTRACT In order to simulate the propagation process of subway vibration of parallel tunnels in semi-infinite rocks or soils,time domain boundary element method (TD-BEM) formulation for analyzing the dynamic response of twin-parallel circular tunnels in an elastic semi-infinite medium is developed in this paper.The time domain boundary integral equations of displacement and stress for the elastodynamic problem are presented based on Betti’s reciprocal work theorem,ignoring contributions from initial conditions and body forces.In the process of establishing time domain boundary integral equations,some virtual boundaries are constructed between finite boundaries and the free boundary to form a boundary to refer to the time domain boundary integral equations for a single circular tunnel under dynamic loads.The numerical treatment and solving process of time domain boundary integral equations are given in detail,including temporal discretization,spatial discretization and the assembly of the influencing coefficients.In the process of the numerical implementation,infinite boundary elements are incorporated in time domain boundary element method formulation to satisfy stress free conditions on the ground surface,in addition,to reduce the discretization of the boundary of the ground surface.The applicability and efficiency of the presented time domain boundary element formulation are verified by a deliberately designed example.

KEYWORDS Time-domain boundary element method;twin-parallel tunnels;infinite boundary elements;elastic semiinfinite medium

1 Introduction

Due to the demanding development and the limited space in urban areas,complex underground structures including the parallel subway tunnels,spaced at a limited distance,are being constructed.Although the subway brings convenience,a moving vehicle on uneven rails produces a dynamic force between wheels and rails and causes vibration,which has potential negative effects on sub-infrastructure beneath the rails,twin-parallel tunnels and even the structures on the ground [1,2].

Without loss of generality,taking the unballasted track in a subway project as an example,the whole propagation process of subway vibration could be roughly divided into two phases:First,the coupled vibration between rails and the wheels of the vehicle acts on the track slab,sequentially passes through the CA mortar layer,the railway ballast,the lining of the tunnel and the outer protective layer of the lining.And after that,the vibration propagates in the semi-infinite surrounding rocks or soils.Since the layout of the parallel subway tunnels is generally in the form of long cylinders in semi-infinite rocks or soils,the subway vibration analysis system could be illustrated by a plane model,as shown in Fig.1.

Figure 1:2-D vibration analysis model

In order to fully understand the propagation law of subway vibration,great efforts have been devoted to developing numerical models for the dynamic track-tunnel-soil system,among which finite element method (FEM) is widely used.However,FEM requires an artificial boundary when dealing with infinite or semi-infinite problems,so it is not the ideal numerical tool.Since boundary element method (BEM) is good at simulating the radiational wave propagation in semi-infinite soils,the FEM-BEM coupled model would be most attractive [3-5].Usually,FEM should be employed to dynamically simulate the track-tunnel system,where vibration forces in the boundary of the tunnel due to the passage of a train could be determined.Then,BEM should be introduced to model the dynamic response of the surrounding soil.

In comparison with FEM,BEM has irreplaceable advantages in dealing with infinite or semi-infinite problems,automatically satisfying Sommerfeld’s radiation conditions at the infinity for infinite and semi-infinite domains,restricting the discretization only on the boundary of the computational domain for elastic analysis [6-10].Besides,reducing the dimensionality for the problem and possessing good computational accuracy are its other advantages [11-13].

BEM for elastodynamic problems could be classified according to the nature of the adopted fundamental solutions [14].The time domain boundary element method (TD-BEM) employs time dependent fundamental solutions [15],Laplace transform method adopts the Laplace transform of the fundamental solutions [16] and the dual reciprocity method (DR-BEM) uses the fundamental solutions of elastostatics [17].It is obvious that TD-BEM provides a direct way of obtaining the time history of the response for transient problems without any domain transforms [18].The use of time domain fundamental solutions turns TD-BEM formulation very elegant and attractive from the mathematical point of view [14,15,18].In addition,the fulfillment of the radiation condition by time domain fundamental solutions turns the formulation suitable for performing infinite domain analysis,since there are no reflected waves from the infinity [19,20].

TD-BEM formulation for elastodynamics was originally proposed for the first time in two journal papers by Mansur et al.[21,22] and after that in Mansur’s Ph.D thesis [23].Subsequently,Carrer et al.[24,25] found the procedure to work out the integral equation for stresses at internal points,where kernels were calculated by adopting the concept of the finite part of an integral,rather than following the kernel regularization procedure presented by Mansur [21-23].Based on the above researches,extensive studies on applying and improving TD-BEM formulation for elastodynamic problems were conducted.Zhang [26] applied TD-BEM for transient elastodynamic crack analysis,and he and his colleagues [27] developed the numerical scheme in orthotropic media.Alielahi et al.[28] used TD-BEM to study the seismic ground response in the vicinity of an embedded circular cavity.Lei et al.[29] presented the analytical treatment to the singular integral terms in boundary integral equations,which improved the computation accuracy of TDBEM for 2-D elastodynamic problems.Several other studies [30-33] were devoted to improving computational stability.

The above mentioned studies are focused mainly on infinite domain problems,where only the discretization on the boundary is required and external boundaries are not required to be artificially truncated.However,many practical geotechnical engineering problems require models associated with the semi-infinite domain.The stress free conditions on the ground surface should be taken into consideration when dealing with semi-infinite problems by TD-BEM [34].Three approaches were proposed for the analysis of the problem with the semi-infinite medium.The first approach is to find semi-infinite fundamental solutions [35,36],where the stress free conditions are automatically satisfied and the discretization of the boundary of the ground surface is not required.However,semi-infinite fundamental solutions are obtained with the assistance of the method of source image,which means that the smooth and regular ground surface is a must.The second approach is to artificially truncate the boundary of the ground surface to the bounded region of the manageable size [37].Therefore,it increases the degree of freedom and the amount of calculation to obtain acceptable results.The other approach is to incorporate the infinite boundary element into the boundary element analysis [38,39],making full use of the simplicity of full plane fundamental solutions and decreasing the discretization of the free surface.Therefore,the last method is the most effective,compared with the other two methods.

The idea of the infinite boundary element and the reciprocal decaying shape function were first proposed by Watson [40],and after that,Zhang et al.[41] incorporated the infinite boundary element in BEM formulation and applied the formulation to the static problem of a dam foundation.Davies et al.[42] developed the semi analytical approach based on the circular region of exclusion in the far field and described the order adaptive criteria for the numerical integrations.The algorithm was further improved in the subsequent paper [43] by means of the analytical integration of the strongly singular traction kernel.All the aforementioned works are valuable efforts to try to incorporate the infinite boundary element in BEM formulation for elastostatic problems.Unfortunately,to the authors’ best knowledge,the TD-BEM formulation,which incorporates the infinite boundary element to deal with dynamic problems,has not been reported in current references.

By removing the FEM sub-model in the ideal FEM-BEM coupled model in Fig.1 for the propagation of the vehicle-rail coupled vibration in the semi-infinite soils or rocks,the remaining model becomes the BEM sub-model,as shown in Fig.2.Therefore,the development of TD-BEM formulation for dynamic analysis of twin tunnels under dynamic loads in an elastic semi-infinite medium is one of the key points in the ideal FEM-BEM coupled model for the dynamic track-tunnel-soil system.

This paper presents the semi-infinite TD-BEM formulation for analyzing the dynamic response of twin-paralleled circular tunnels under dynamic loads.Time domain boundary integral equations for this problem are presented,and then,the numerical implementation by introducing the infinite boundary element is presented.Finally,the presented TD-BEM formulation for twin-parallel tunnels in an elastic semi-infinite medium,in which the infinite boundary element and time domain fundamental solutions are incorporated,is verified by a deliberately designed example.

Figure 2:TD-BEM model in the vibration analysis system

2 Boundary Integral Equations in an Elastic Semi-Infinite Medium

2.1 Displacement Boundary Integral Equation

The boundaries for the semi-infinite domain problem include finite boundariesΓ1andΓ2,the boundary of ground surfaceΓ3and nominal boundaries at the infinityΓ∞,as shown in Fig.3.By referring to [35],the displacement boundary integral equation for the elastic,isotropic and homogeneous semi-infinite domain containing twin-parallel circular tunnels based on Betti’s reciprocal work theorem,ignoring contributions from initial conditions and body forces,can be expressed as:

wherecij(P)is a parameter that depends on the position of the source pointP;uk(P,t)andpk(P,t)are displacement and traction components in thekdirection,respectively;u∗ik(P,τ;Q,t)andp∗ik(P,τ;Q,t)stand for time domain displacement and traction fundamental solutions respectively,representing the displacement and traction at the field pointQin thekdirection attinstant due to a unit point force applied at the source pointPin theidirection at any time instantτ.

Figure 3:Schematic diagram of the semi-infinite domain with twin-parallel tunnels

It is noted that the boundaries AB,BA,CD,and DC are virtual boundaries,and their relations are shown in Eqs.(2) and (3).The infinite terms in Eq.(1) can be simplified in Eqs.(4)and (5).

Thus,the displacement boundary integral equation can be re-formulated,as:

where Γ=Γ1∪Γ2∪Γ3.

The displacement and traction fundamental solutions in Eq.(6) can be formulated by Eqs.(7)and (8),according to [24,29]:

The corresponding integrals of fundamental solutions are expressed by Eqs.(9) and (10):

In Eqs.(9) and (10),ρis the density,andris the distance between the field pointQand the source pointP.csis the shear wave velocity orSwave velocity,andcdis the compression wave velocity orPwave velocity.These parameters are expressed as follows:

where

The tensorsEik,FikandJik(in Eq.(9)) andAik,BikandDik(in Eq.(10)),which are related only to space,are defined by:

where

The parameters that are relevant to both space and time in the Eqs.(9) and (10) are expressed,as:

It should be noted that the subscriptwcould be substituted withsordin the Eqs.(26)-(28),representing the contribution of the shear wave (Swave) or the compression wave(Pwave),respectively.

In the following context,all the finite parts of the singular integrals are analogously calculated.

2.2 Stress Integral Equation

The stress integral equation for internal points is obtained by combining Hooke’s law with the derivatives of Eq.(6),which is expressed,as:

where In references [24,25],the stresses at surface nodes are computed using simple finite difference formulae,which employ the displacements and tractions computed by the TD-BEM formulation.The stress integral equation at internal points is expressed,as:

whereDpwandSuwrepresent the contribution of shear and compression waves to the integrals related to the tractions and displacements,respectively.According to references [24,25,29],the terms are expressed in Eqs.(34)-(37):

The tensorsEijk,Fijk,Jijk,,and(in Eqs.(34) and (35)) andAijk,Bijk,Dijk,,and(in Eqs.(36) and (37)),which are related to space,are expressed,as:

By recalling the expressions,Eqs.(11)-(30),these tensors can be calculated,and the boundary integral equation in terms of the internal stress Eq.(31) is determined.

3 Numerical Implementation

3.1 Temporal Discretization

The displacement is assumed to have the linear variation in time,whereas the traction is assumed to be constant.For the given time interval [tm−1,tm],the variables of the displacement and the traction can be written:

where the time interpolation functions are defined,as:

3.2 Spatial Discretization

Fig.4 shows the boundary discretization and the distribution of boundary elements for twin-parallel tunnels in a semi-infinite domain,including infinite boundary elements.The finite boundaries,Γ1andΓ2,are linear boundary elements.Quadratic boundary elements are adopted to describe infinite boundary elements in spatial discretization.

Figure 4:Schematic diagram of boundary discretization

From the expressions of fundamental solutions,Eqs.(7) and (8),it can be seen thatu∗ikandp∗ikare the same orders as the functions of 1/rand 1/r2,respectively.Whenrapproaches∞,meaning that the field pointQtends to the infinity,u∗ikandp∗ikare the infinite small functions with the same orders as the functions of 1/rand 1/r2,which could be mathematically expressed,as:

The displacement solution and traction solution can be obtained by the superpositions of displacement and traction fundamental solutions respectively,so decay rates ofuk(Q,τ)andpk(Q,τ)also possess the same orders as the functions of 1/rand 1/r2,respectively,when the field pointQtends to the infinity,as:

Therefore,besides the satisfaction of the convergence criterion,the interpolation mode of the infinite boundary element needs to reflect the above mentioned decaying characteristics at the infinity.The requirements about the shape functions of infinite boundary elements are

(i) Delta function property,

(iii) When the field pointQis at the infinity,Qk(x)→∞,uk(x)→0,pk(x)→0

There are left and right forms for the infinite boundary elements,as shown in Fig.5.The shape functions of the displacement and traction for infinite boundary elements are respectively expressed in Eqs.(56) and (57).

Figure 5:Infinite boundary elements models (a) left (b) right

It is obvious that the above shape functions possess delta function property and completeness.Moreover,whenξapproaches−1 in the left infinite boundary element,meaning that the field pointQtends to the infinity,the displacement solution and the traction solution have the following relations:

For the right infinite boundary element,these relations are also applicable whenξapproaches 1.It has been proved that the above shape functions meet the above requirements for infinite boundary elements.

It is also noted that the stress on the ground surface is free,i.e.,pk(Q,τ)is null.One has,

The linear interpolation functions of the displacement and the traction in space can be expressed as:

Therefore,the variables of the displacement and the traction after discretization in both time and space can be rewritten as:

whereNqrepresents the number of nodes in the boundary element.

So far,according to the decaying characteristics of fundamental solutions from the finite boundary to the infinity,the shape functions for infinite boundary elements are constructed to mathematically describe the relationships between the infinite boundary and the infinity for the field variables.Moreover,the shape functions satisfy the convergence criterion.

3.3 Discretization of Boundary Integral Equations

The numerical solutions to Eqs.(6) and (31) are obtained after discretization both in space and time.The boundaries are discretized intoNboundary elements,and the time interval of the analysis [0,t] is divided intoMtime intervals,each with the same length Δt.The terms of the displacement and the traction in boundary integral equations are respectively expressed as follows:

In Eq.(63),the coefficientsandrespectively describe the influences of the traction and displacement at nodeaof elementeat instantmΔton the displacement of the source point.In Eq.(64),the coefficientsanddescribe the influences of the traction and displacement at nodeaof elementeat instantmΔton the stress of the source point.

The values of influencing coefficients are calculated by using the effective method in [29],where the analytical treatments of the singular integral terms were presented,and to which the interested readers could refer for details.

3.4 Matrix form of Boundary Integral Equations

The influencing coefficients at nodeaof elementeat instantmΔtare assembled based on one-to-one correspondence in time and space.After that,boundary integral equations are expressed in matrix forms,and the programs of the TD-BEM formulation could be coded.

The influencing coefficients of the displacement at nodeaof elementeat instantmΔtare obtained by adding the influencing coefficients caused by the (m+1)th time step and those caused by themth time step,as shown in Eqs.(65) and (66)

where the superscripts1and2respectively refer to the forward temporal node in the (m+1)th time step and the backward temporal node in the (m)th time step.

To further simplify,the influencing coefficients of linear spatial elements are also assembled.The influencing coefficients of elementeare obtained by adding the influencing coefficients at the first node of elementeand those at the second node of element (e−1).The specific assemblage equations are expressed as follows:

whereNei(i=1,2,3) represents the total number of elements for the corresponding boundariesΓ1,Γ2andΓ3.

When the source pointPtraverses all boundary nodes,the matrix forms of boundary integral equations after the assemblage can be expressed as follows:

where

In Eq.(71),uMandpMdenote the displacement and traction vectors at instantMΔtrespectively,and matrixCis constituted by thecik(P)coefficient.HandGare boundary influence matrices corresponding to the coefficientsandrespectively.In Eq.(72),σMis the vector representing the stress components at instantMΔt,and matricesSandDare correspondingly constructed by coefficientsandrespectively.

4 Numerical Verification and Discussion

This section is to verify the applicability and the accuracy of the presented TD-BEM formulation for twin-parallel tunnels in an elastic semi-infinite medium,where the infinite boundary elements are incorporated and the time domain fundamental solutions are adopted.

For the purpose of verifying the algorithm,the example with analytical solutions should be the best,because of the uniqueness and the objectivity of the analytical solution.Unfortunately or fortunately,although the analytical solution to the response of the twin-parallel tunnels under dynamic load in an elastic semi-infinite medium is not available,the analytical solution to the response of the single tunnel under dynamic load in an elastic infinite medium is available.Therefore,the numerical example for twin-parallel tunnels is deliberately designed to utilize the analytical solution to the response of the single tunnel under dynamic load in an elastic infinite medium to verify the presented TD-BEM formulation for twin-parallel tunnels under dynamic load in an elastic semi-infinite medium.As shown in Fig.6a,two parallel circular subway tunnels(r0=3 m) are deliberately assumed to be buried in enough deep depth (h=30 m) beneath the ground surface and are spaced at enough faraway distance (d=120 m).Monitoring point A is positioned at a nearer distance (6 m) to the center of the left circular tunnel than to the right.The twin-parallel tunnels are assumed to be subjected to the same Heaviside uniform loadsp=200 MPa,shown in Fig.6b,suddenly at the same time.The mechanical parameters of the surrounding soil are:the elastic modulusE=8×109Pa,Poisson ratioν=0.25 and the densityρ=2350 kg/m3.

Figure 6:Sketch of the example (a) geometry of tunnels and distribution of loads and(b) load definition

In this scenario,theoretically,the response of monitoring point A could be classified into the early phase and the late phase.In the early phase,the wave initiated from the left tunnel has arrived at and passed by point A,while the reflected wave from the ground,due to the wave from the left tunnel,has not yet arrived,let alone the more faraway wave from the right tunnel.In the late phase,the wave at monitoring point A is subsequently disturbed by the reflected wave from the ground and by the propagating wave from the right tunnel.Therefore,the response of monitoring point A for the problem of twin-parallel tunnels in half-plane in the early phase should be the same as the response of the problem of a single circular tunnel in an infinite medium.According to the analytical wave velocity (v=2021 m/s) in the surrounding ideal elastic soil,the early phase is from 0.000 s to 0.029 s,and the late phase is from the instant of the occurrence of the disturbance of the reflected wave from the ground at the instant of 0.029 s to the instant of the analysis,during which the additional disturbance of the wave from the right tunnel joins at the instant of 0.055 s.Therefore,the response of monitoring point A in the early phase for the problem of twin-parallel circular tunnels in the semi-infinite medium from the presented TD-BEM could be compared with the analytical response from the method of characteristics for the problem of the single circular cavity under dynamic load in an infinite medium given by Chou et al.[45] to verify the presented TD-BEM formulation.It is noted that the analytical solution of the single circular cavity in an infinite medium in [45] has been used to verify the TD-BEM formulation in several references [24,29],to which the interested readers could refer for more information and mutual comparison.

Figure 7:Comparison between results at point A (2ro,0).(a) Radial stress σr.(b) Circumferential stress σt.(c) Velocity vr

Figs.7a-7c show the comparison between the calculated results from the presented TDBEM for the problem of twin-parallel circular tunnels in an elastic semi-infinite medium and from the method of characteristics for the problem of single circular cavity in an infinite medium in terms of radial stress (σr),circumferential stress (σt) and velocity (vr) respectively at monitoring point A.

From Figs.7a-7c,it can be seen that,in the early phase (0.000-0.029 s),good agreement between calculated results from the presented TD-BEM and those from the method of characteristics is obtained.This observation quantitively invalidates the applicability and accuracy for the presented TD-BEM formulation for twin-parallel circular tunnels in an elastic semi-infinite medium in the early phase.It can also be seen from the figures that,in the late phase (0.029-0.14 s),the calculated results from the presented TD-BEM fluctuate around the analytical results from the method of the characteristics,with dramatic differences in general in the whole late phase and with sudden changes at the occurrence of the disturbance of the reflected wave from the ground surface and at the occurrence of the additional disturbance of the wave from the right tunnel,indicating the effects of the reflected wave and the interaction of the waves.Meanwhile,good agreement between the calculated instants of the occurrences of the disturbances from the presented TD-BEM and the corresponding analytical instants are obtained.These two observations qualitatively invalidate the applicability and accuracy for the presented TD-BEM formulation for twin-parallel circular tunnels in an elastic semi-infinite medium in the late phase.

5 Conclusions

This study is concerned with developing TD-BEM formulation for dynamic analysis for twinparallel circular tunnels in a semi-infinite medium,with the following conclusions:

• Time domain fundamental solutions are adopted in the boundary integral equations in TDBEM formulation for twin-parallel tunnels under dynamic loads in an elastic semi-infinite medium.The time domain fundamental solutions make it possible to directly and effectively obtain the time history of the response in the presented TD-BEM formulation.

• In the process of the numerical implementation,infinite boundary elements are incorporated in TD-BEM formulation to satisfy the stress free conditions on the ground surface,moreover,to reduce the discretization of the boundary of the ground surface.

• The applicability and the accuracy of the TD-BEM formulation are verified by comparing the calculated results from the presented TD-BEM with those from the method of characteristics in the deliberately designed example.

• The presented TD-BEM formulation could be a step forward in the ideal FEM-BEM coupled model for twin-parallel tunnels under dynamic loads in an elastic semi-infinite medium and could be applicable for blasting projects with twin-parallel blasting holes.

Funding Statement:The authors would like to acknowledge the financial support from the research Grants,Nos.2019YFC1511105 and 2019YFC1511104 by National Key R&D Program of China,No.51778193 provided by National Natural Science Foundation of China,and No.QN2020135 provided by Hebei Education Department.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.