Coupled beam-boundary-element model (BE-BEM) for analysis of underground openings

Computers and Geotechnics 2 (1986) 239-256

couPu~ B m m - s o ~ , - m m ~ a ~

. o m ~ (sE-am) FOR ANALYSIS OF [~NDE]M3ROUND OPENINGS

R. P~ttler ILF Consulting Engineers Innsbruok Austria and G. A. Swoboda Institut f"dr Baustatik University of Innsbruck Austria

ABSTRACT

The coupled beam-boundary-element (BE-BEM) model D ,2] is an efficient computation model to determine the support behavior of the rook mess and the support measures in a realistic way without requiring excessive computation time. In order to avoid time-consuming FE-ccmputations the "coefficient of subgrade reaction method" (CSRM) was often applied for stability analyses in tunneling especially in GermAny [33. The difficulties in using the CSRM are well-known, a fact which often gave rise to discussions about the input data and the computation results obtained. In this paper the results of a BE-BEM computation will be compared with the approxi,mte solution obtained by the CSRM. It will be shown that the time-intensive FE-BEM model [4,5] can be replaced by the BE-BEM model as far as the support behavior of the rook mess is predominately elastic.

THEORY

Coupling of Boundary Elements (BE) with beam elements is done on the same basis as coupling of BE with Finite Elements (FE) [4~5]. For the BE region a fully populated, but symmetrical stiffness matrix was developed.

[[c1

Ecl]

239 Computers and Geotechnics 0266-352X/86/$03.50 © Elsevier Applied Science Publishers Ltd, England, 1986. Printed in Great Britain

(1)

240

[c] : [B] -1 [A]

[A]

(2)

and [B] are the matrices in the well-known BE equation

[A] Ial= EB] Ibl [M] is a banded

matrix and is obtained

(3)

by assembling

the single

element

matrices [M] e.

E

[.] : ~

(4)

[M]e

e= I

When coupling beam-boundary element, linear shape functions for the two-node boundary elements (4 degrees of freedom) are chosen (Table I).

TABLE I Shape

Functions of Beam and Boundary Element.

BEAM-ELEMENT

$HAPE~

NODE 1

NI, N.

p

L

NODE 2

BOUNDARY-ELEMENT NODE 1

NOOE 2

u:l

-I

CONFORM

nJrrn~ Nz,Ns 1-31~2.2~3

3~ 2- 2~"I

~

w=l NONCONFORM

~o=1

N3,NI 'L(1~-2¢ 2. I~3)

L(-¢2o1~ 3 )

241

NI=

N4= 1

I -~

-~

(5) N2=

N s --

Defining IN] e

[:1o.20][:oN20] [N]e

N.

0

N5

NI

0

c6

N2

it follows

[H]e

.-/[~]e,~ ["7 dx:

C7)

[N21 [I] .22 where [I] is the 2 x 2 unit matrix for the two-dimensional problem. Nij is obtained by integration:

Nij =

f N i Nj dx

(8)

Using equations (5) and (8) the submatrix [M]e becomes

1

0 1

sym.

1

o

I

~

0

(9)

1

0 I

A two-node beam element with 6 degrees of freedom is used. The "Integrated Beam Element" considers in the element stiffness matrix the nonlinearity of the concrete, reinforcement, cracking and interaction between concrete and rein-

242 forcement

(tension-stiffening).

Cracks

are

"smeared

out"

in a continuous

fashion along the element length so that cracking does not cause geometrical changes of the system.

The stiffness matrix of the beam elements can be written as follows:

[K]B~M : /[B] ~ [D]pl [B] dx

(I0~

The coefficients of [Dip I depend on the stress level of the cross-section and in this way on the material law ~,2]. Y

¥

-X

~X

FIGURE I. Two-node constraint element.

Coupling of the beam element with the BE region is done by means of twonode constraint elements D ~2,6] (Fig I). Besides the displacements additional constraints in the form of forces k,

and ~s

are included. These act between

the two nodes. The additional constraint equations can be expressed as:

[c] lal -Ibl:0 [c]

oonstrai~t coefficient

lal

displacement vector

(11)

~tri~

[1,2]

Ibl specifiedconstants [1,2] Coupling of equation (11) with the standard finite element displacement formulation

243

lal -Iql

-- o

(12)

leads to the following general equation in matrix formulation:

I:l-I:

= o

(13)

The nonlinear behavior of the constraint element is governed by the MohrCoulomb failure criterion.

[C] and

Ibl depend on the actual state of the

constraint element. For the three possible interface states (free, sllp, fix), [C] and Ibl are defined in references [1,2].

By coupling FE and BE, for finite elements and the boundary elements the same

shape functions are used.

When coupling boundary elements with beam

elements the first degree Hermitian polynomials are used as the shape function for the beam elements for rotation and displacement in cross direction to the local beam axis and only for the displacement in the direction of the local beam axis a linear shape function is used (Table I). In consequence geometric compatibility is given only in the nodes while

there

is a gap

in cross

direction of the beam axis. In order to be able to evaluate the accuracy of the computation result, it is necessary to know the size of this gap at any point (Fig 2).

FIGURE 2. Gap between beam and boundary element.

244 In references D,2] the formula of the gap between boundary element and beam element is developed.

w(X)BE_BEM = ( ~ - 3 ~2 + 2 ~3)wI + L ( ~ - 2 ~2 + ~3)

+(-,~ +3,~2_2

~I

(14)

~3)w2 + L(- ~2 + ~3) ~2

By integration over the length of the beam and by summation over the number of elements E the area of gap is defined as follows: E ABE_BEM :

Numerical integration is done

E e:1

V / W ( x ) 2BE_BEMdX '

using

the GauB-Legendre.

(15)

The

parameter of

accuracy related to the length of the beam (Ie) is

a =

ABE-BEM E :C e=l

Value

a

(16)

ie

(equation 16) divided by the total related displacements in cross

direction of beam resp. boundary elements (WBE resp. WBE M) leads to the values

aB~ and aBEM"

a

E ~ e=1

aBE=E

ie J/

2 ' w(X)BE

(17) dx

e=1

aBE and aBEM show the area of gap expressed as a percentage of the total displacement.

245 E a

C

ie

e=1 aBEM=

E

j/

2 1 w(X)Bg M

(18) dx

e=1 In order to be able to evaluate the parameter of accuracy due to the nonconform

shape

functions

a

lined

tunnel

cross-section

was

studied

for

different loads (Fig 3) as a function of the element size.

LOAD CASES 1

2

3

4

FIGURE 3. Assessment of accuracy: Load configuration.

TABLE 2 Assessment of Accuracy: Input Data for Parameter Study.

Shotcrete Rock Mass

E

~

d

MN/m 2

-

m

1000 1

0.17 0.20

0.05 -

Figures 4, 5 and 6 show the characteristic values for different element sizes (q, 8, 16 and 32 subdivisions) according to equation (16), equation (17) and equation (18).

246 m 2 (xlooo)

FIGURE 4.

Assessment of accuracy: Total area of gap related to the length of beam (a), equation (16). aKu(X)

~2---£~ NUMBER OF ELEMENTS

FIGURE 5.

Assessment of accuracy: Total area of gap related to the boundary element displacement (aBEM) , equation (17). a H [X)

"°I 20, i

NUMBI~q OF ELEMENTS

FIGURE 6.

Assessment

of accuracy:

Total

area

of gap

element displacement (aBE) , equation (18).

related

to the beam

247 By dividing the cross-section into 16 elements the error is, due to the nonconforming shape functions, less than 5 % independent of the load configuration. As expected, the results for load figuration 4 (internal pressure) are, independent of the number of elements, of utmost accuracy compared with the other load cases. Obviously the Eap is the smaller, the smaller the bending moments are. AnaloEous results are obtained if, instead of the area of gap, the absolute and relative

characteristic values of displacement w(x) are con-

sidered. The "apparently" high accuracy of load case I with 4 subdivisions is caused by the stiffness ratio of rock mass and beams as no bending moments occur.

The

efficiency

of

the coupled

BE-BEM model is shown by the results

obtained from a comparison of load case 4 (with different subdivisions) with the exact solution [7] in Fig 7. For 8 subdivisions the error for the radial displacements of the BE-BEM solution is 8 . 1 % , for 16 subdivisions only 2.2 % and for 32 subdivisions 0.5 %.

mm

Q~C T

A::.|,/,

,~=a.T'/. A:I.I'I. ~o..44

LEXACT 8OLUTION REF.I 7|

D

NUMBEROFELEMENTS FIGURE 7.

Assessment of accuracy: Comparison of results of load case 4 with the exact solution.

~ARISON

OF I B - B E M W I T H

CSRM

Using the CSRM (coefficient of subgrade reaction method) the support of the

shotcrete

against

the

rock

is

determined

by

a

supplementary

matrix

248 [K] B in addition to the stiffness m t r i x

[K] of the beam. For a bedding in

optional direction to the beam axis, the fully populated supplementary matrix [K]B is formulated in reference [2].

By means of the following example of a stability analysis for the heading with an invert arch installed immediately after excavation the calculation is compared with a coupled beam-boundary element

system using the CSRM.

The

assumed load configuration is verified by the measuring computations which have been performed during tunnel construction in the course of an express railway

line

in Germany [8,9].

In Fig 8 the

load

configuration

and

the

dimensions of the system are shown.

L

FIGURE 8.

WlIIII

~' : 0 . 0 2 0 M N / m !

°

Load configuration and dimensions of system.

The beam elements simulating the shotcrete lining consider the nonlinear material behavior of reinforced concrete [2]. Here the tension stiffening is considered by modifying the stress-strain curve of concrete according to ~0]. In the computation the reinforcement is taken into account according to its location using the bilinear stress-strain curve according to DIN 1045. Two layers

of

(3.77 cm2/m,

welded

steel

fabric

mats

Q 377

BSt

500 M

are

installed

a • = 500 N/mm2) whereas the concrete cover is 6 cm (Fig 9). e

249

~

118111~mll

O" [ ml/m=] CONCRETE

1=

FIGURE 9.

3.5

2.

EEL

d-

= ~ tl..]

Stress-strain curve of concrete and steel.

(Assumptlc~s A, B , C)

The rock mass supports the outer lining in radial direction. Tangential bedding is not considered modulus

of subgrade

according

reaction

to

[3].

is calculated.

From the soll Since

parameters

the modulus

the

of subgrade

reaction is not a material parameter and, apart from the structural geometry, is dependent on the load configuration, an exact determination is possible only for trivial cases. In literature approximations are given for tunnel crosssections.

different assumptions were studied [2]:

Three

A:

The radial modulus of subgrade reaction is defined as follows: ES K

=

-

(20)

-

R

m

with

ES=

E .

I -

V

I-

P-2P

2

1

Rm = ( A ) ~

A

... ..... average radius

.................

area

of

excavation

250 SYSTEM , i 1¢m : 1.0m

M-h N-h

M-bl

I

,y-b ~

M,-!

COEFFI~ENT OF SUBGRADE

.sA~,o. x [m~,]

FIGURE 10. Moduli

I

of subgrade

J

.=, ~ "~' i /

reaction

/ A

21.7

C

4.,

s

\ \

=1.T

/

~.1

for assumptions

\

A, B,

C. Notations

Table 3 and Table 4.

B!

In the area of high curvature the radial modulus of subgrade reaction is doubled compared with the modulus determined by equation (20).

C:

The

radial

modulus

of

subgrade

reaction

is

determined

according

to

equation (20), R m is replaced by the actual radius at any point.

The different distributions of the moduli of subgrade reaction (K) are shown in Fig 10 using E = 100 MN/m 2 and

B-IEH

(&ssumptioo

v = 0.20 for calculating K.

]))

For the BE-BEM calculation the stiffness matrix of the two-dimensional (plain-straln)

fullspace is used. The modulus of elasticity is chosen with

E = 100 MM/m 2, the Poisson's ratio with

Y = 0.20, analogous to assumptions A,

B~ C. In this calculation the boundary elements are only used for simulating the support of the lining and not for simulating the rock mass and its bearing capacity as done in the following chapter. The interaction between shotcrete and rock mass is simulated by constraint elements with nonlinear behavior [2]. According to the statements given above tangential forces are not transferred.

The

computation

results

numerically in Table 3.

are

graphically

shown

in Flgs

11 and

12 and

251

. . . . . . . . . . .

/

FORCE ~

/

1

SYSTEM FORCE

FIGURE 11. Comparison of results: Moments and normal forces.

L

.....

C

FIGURE 12. Comparison of results: Displacements.

TABLE 3 Comparison of CSRM and BE-BEM: Results.

A

vy-heading vy-invert vx-bench vx-bench M-heading N-heading M-invert N-invert M-bench N-bench

93 4 25 44 0.254 -1.208 -0.012 1.558 -O. 379 -I .749

CSRM B

64 5 15 30 0.218 -1.274 0.020 -1.696 -0.228 1.799

C

73

3O 12 37 0.222 -1.261 -0.101 -I .634 -0. 125 -1.732

BE-BEM D

91 16 22

45 0.212 -I .247 0.031 1.603 -0.213 -I .725

mm mm mm mm

e~Tm/m HN/m MNm/m

HN/m t~aqm/m MN/m

252 The BE-BEM model provides the exact solution D as here the soil is described as continuum and not simulated by independent springs as in the CSRH. As shown in Table 3 and in Figs 11 and 12, assumption B often used in practice shows, all together,

the best agreement

for this geometry.

As expected,

assumption

A

taking a general average radius as the basis for the calculation of the modulus of subgrade reaction shows little conformity with the exact solution due to the special geometry. Only the normal forces and the displacements in the heading area show realistic actual

results.

The

radii at the appropriate

same is valid

for assumption C where

point are used to determine

the

the modulus of

subgrade reaction.

COMPARISON ~

BE-BEM WITH FE-B~q

For a shallow tunnel (Fig 13) the FE-BEM calculation was compared with the BE-BEM in reference

Eli. Apart from an excessive saving in computer time and

computation time there was considerable conformity in the results obtained by BE-BEM and FE-BEM. Both, the normal forces and the displacements of the shotcrete lining were studied for two different load cases (Fig 14).

-20

-15

-10

I

I

I

-5 I

0 I

5 t

10 I

15 I

i

20 I

Soil t k., 0.50 ;, • 0.33 E~ 7MN/ma 7,0.020 MN/m) Soil 2 ko: O.S0 v : 0.33 EL= 140 M N / m 2 7 = 0.022 MN/m ] F_u:2 8 O M N / m z (" = 35" c = 0.010 MN/m z

-10

-is

BOUNDARY

FEM

FIGURE 13. Geological situation.

-

BEM

Model

Elements

253 In view concrete

of

tensile

the

results,

stresses,

minimum

it was

moments

not

which

necessary

do not bring

to consider

about

the

any

nonlinear

behavior of the concrete in this case. Constraint elements simulate the interaction

between

shotcrete

c = 0.100 MN/m 2.

These

and

soil:

values

angle

correspond

of to

friction

the

shear

~ = 35 ° , parameters

cohesion of

soil

layer 2. The density of the soll was 7 = 0 . 0 2 1 M N / m 3 (mean value from layers I and 2), Poisson's

ratio

y = 0.33, coefficient

~=O.~h

FIGURE

14. Load configuration

In addition

of lateral pressure k 0 = 0.50.

I~ = O . ~ h

A and B

to the calculation

1 .

as discussed

in [I], the same system was

computed assuming a loading modulus of E 1 = 140 MN/m 2 and an unloading

modulus

of E u = 280 MN/m 2. The results (Fig 15, Table 4) showed that, as expected, there is bad conformity of the BE-BEM and the FE-BEM calculation in this case. The

distribution

different. varying

of the normal

Consequently

it

forces

is not

the modulus of elasticity

and

the displacements

possible

to obtain

are completely

better

conformity

by

of soil (e.g. average value of loading and

unloading modulus).

Compared

with

the FE-BEM

it is shown

that,

behavior of the soil can be sufficiently

described

elasticity,

fully

the

FE-calculation

coupled beam-boundary

can

be

element calculation.

in cases where

the support

by means of the theory of

replaced

by

the

time-saving

This is not the case if the support

behavior of the soil is characterized by anisotropy, very different loading and unloading

moduli or extensive plastic zones.

254

TABLE 4 Comparison of Results: Displacements and Normal Forces

Results

vy-heading vy-lnvert vx-bench vy-bench N-heading N-invert N-bench

BE-BEM

FE-BEM 140

(mm) (mm) (mm) (mm) (MN) (MR) (MR)

7.7 -7.2 2.5 2.1 -0.513 -0.317 -0.887

A

B

7.3 -7.4 2.2 1.8 -0.513

7.9 -8.0 2.0 3.1 -0.q86 - 0 . 320 -0.886

-0.465 -0.895

~

-

BE-

% = FE-BEM BE-BEM A B 5 2 14 0 47

2 10 20 47 5 1

1

0

11

FE-BEM 140/280

6.2 -3.6 1.0 2.1 -0. 396 -0.248 -0.814

BEM/A

.....

e~-.E./e

. . . . .

FE - BE/140

. . . . . . . . . . . PE- e E / . o / 2 o o

FIGURE 15. Comparison of results: Displacements and normal forces.

S~g46RY

The coupled beam-boundary element system (BE-BEM) is an efficient compu-

tation model combining the advantages of little computation time, little preparation time for input data and exact determination of the soll behavior. In contradistinction to the CSRM~ the BE-BEM will lead to exact results even with nearby the same computation time.

255 ACKBOgLE~S

The development of the Coupled Beam-Boundary Element Model was sponsored by ILF Consulting Munich.

Engineers

The

solution

(Ingenieurgemeinschaft of

the

Boundary

L~sser-Feizlmayr),

Stiffness

Matrix

was

Innsbruckfinancially

supported by the ~sterreichische Nationalbank (PrJ. Nos. 2147 and 2529).

REFEBENCES

I.

PSttler, R., Swoboda, G. A. and Beer G., Nonlinear coupled beam-boundary element (BE-BEM) analysis for tunnels on microcomputers. In: Microcomputers in Engineering: Development and Application of Software, Ed. by B. A. Schrefler and R. W. Lewis, Pineridge Press, Swansea U.K. (1986) 161-176.

2.

P~ttler, R., Gekoppeltes Stabwerk-Boundary-Element-Modell zur Berechnung von Tunnelbauten unter BerUcksichtigung des nichtlinearen Materialverhaltens von Stahlbeton. Dissertation, Universit~t Innsbruck (1986).

3.

Duddeck, H., St~ding, A. und Schrewe, F., Zu den Standsicherheitsuntersuchungen fur die Tunnel der Neubaustrecke der Deutschen Bundesbahn. Felsbau 2 (1984) 143-151.

4.

Beer, G. and Swoboda, G. A., On the Efficient Analysis of Shallow Tunnels. Computers and Geotechnics I (1985) 15-31.

5.

Beer, G., Implementation of Combined Boundary Element-Finite Element Analysis with Application in Geomechanics. Developments in Boundary Element Methods 4 (1986) 191-226.

6.

Katona, M. G., A Simple Contact-Friction Interface Element with Applications to Buried Culverts. Int. Journal for Numerical and Analytical Methods in Geomechanics 7 (1984)

371-384. 7.

Windels, R., Kreisring im elastischen Kontinuum. Der Bauin~enieur 42 (1967) 429-439.

8.

John, M., P6ttler, R. und Heissel, G., Vergleich verschiedener Rechenmethoden ergebnissen. Felsbau 3 (1985) 21-26.

untereinander

und mit den MeE-

256 9.

Niedermeyer, S., Rahn, W., Reik, E. und Stoll, R. D., Geologlsche Gegebenheiten im Unteren und Mittleren Buntsandstein und ihre geotechnische Bedeutung fur Baumagnahmen im Fels. Felsbau I (1983) 54-67.

10.

Bergan, P. G. and Holand, I., Nonlinear Finite Element Analysis of Concrete Structures. Computer Methods in Applied Mechanics and Engineerin~ 17/18 (1979) 443467.

11.

Duddeck, H., Empfehlungen zur Berechnung yon Tunneln im Lockergestein. Die Bauteehnik 57 (1980) 349-356.

Received 29 May 1986; accepted 30 July 1986