Uncertainty modelling and sensitivity analysis of tunnel face stability

Structural Safety 32 (2010) 402–410

Contents lists available at ScienceDirect

Structural Safety journal homepage: www.elsevier.com/locate/strusafe

Uncertainty modelling and sensitivity analysis of tunnel face stability Wolfgang Fellin a, Julian King b, Ansgar Kirsch a, Michael Oberguggenberger c,* a

Unit of Geotechnical and Tunnel Engineering, University of Innsbruck, Innsbruck A-6020, Austria Breath Research Unit, Austrian Academy of Sciences, A-6850 Dornbirn, Austria c Unit of Engineering Mathematics, University of Innsbruck, Innsbruck A-6020, Austria b

a r t i c l e

i n f o

Article history: Available online 2 July 2010 Keywords: Tunnel face stability Model error Sensitivity analysis Monte Carlo simulation

a b s t r a c t This paper proposes an approach to the choice and evaluation of engineering models with the aid of a typical application in geotechnics. An important issue in the construction of shallow tunnels, especially in weak ground conditions, is the tunnel face stability. Various theoretical and numerical models for predicting the necessary support pressure have been put forth in the literature. In this paper, we combine laboratory experiments performed at the University of Innsbruck with current methods of uncertainty and sensitivity analysis for assessing adequacy, predictive power and robustness of the models. The major issues are the handling of the twofold uncertainty of test results and of model predictions as well as the decision about what are the influential input parameters. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction This article addresses the question of model choice and model adequacy in engineering design, especially in geotechnics. Experimental and mathematical methods will be combined to achieve this task. In fact, various types of simplifications and assumptions have to be introduced in geotechnical calculations. This can lead to different models for the same geotechnical problem. These models do not predict the same system behavior in general. The question arises how to assess adequacy, predictive power and robustness of the models. We set out to investigate this issue using laboratory data on the one hand and methods from uncertainty analysis on the other hand. Predictive power can be assessed by comparison of experimental results and theoretical prediction. Here the uncertainty lies in the experimental results, the input data of the models and the propagation of uncertainty to the theoretical output. Robustness and adequacy of the models can be best understood by means of sensitivity analysis [1–4]. When combining experimental data and theoretical models, sampling based sensitivity analysis – with its recently developed powerful statistical indicators – suggests itself as a suitable approach [5–10]. As a further important tool in the assessment of the joint uncertainty of the model parameters, we employ bootstrap resampling techniques [11–14]. The construction of shallow tunnels is an engineering challenge up to the present day. Tunnels with low cover are often headed using the shield technique. In this context the face stability is an * Corresponding author. Tel.: +43 512 507 6824; fax: +43 512 507 2941. E-mail address: [email protected] (M. Oberguggenberger). 0167-4730/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.strusafe.2010.06.001

important issue. In order to minimize settlements at the ground surface and to prevent failure of the soil ahead of the face, the tunnel face must be supported. It has been a long-standing topic of research how to predict the necessary support pressure for shield tunnelling. A variety of theoretical and numerical models for estimation of the minimum required support pressure have been proposed. The theoretical approaches can be subdivided into kinematic approaches with failure mechanisms (e.g. [15–27]) and static approaches with admissible stress fields (e.g. [23,28]). Some additional approaches are neither purely kinematic nor purely static [29,30]. We will use some of these models to exemplify the proposed strategy for assessing the predictive power of a geotechnical model. Experimental investigations of face stability range from experiments at single gravity, so-called 1g-model tests (e.g. [31–34]) to centrifuge tests at multiples of g (e.g. [35–40]). Large scale tests are rare (e.g. [41]). We use a series of 1g-model tests [33,42] for comparison with the prediction of the chosen theoretical models. In the theoretical models under scrutiny, the output parameter was the necessary support pressure ps. The input (soil) parameters possessing the largest degree of random variability were identified as the actual void ratio e and the loose and dense state void ratios el and ed, respectively. All these parameters were estimated in smallscale laboratory experiments. Another important model parameter is the friction angle u of the soil. This parameter is estimated by means of a linear model u  b0 + b1Id, with the relative density Id (which in turn is a function of e, el and ed). In order to assess the influence of the regression coefficients b0, b1 on the output ps, we needed to determine their statistical distribution. We achieved this by means of the so-called resampling technique, producing a large bootstrap sample of the experimental data and thereby simulating

403

W. Fellin et al. / Structural Safety 32 (2010) 402–410

the joint distribution of b0 and b1. We believe that this is a novel method for obtaining joint distributions – including correlations – of geotechnical data. As a first application of the statistical data model, we could assess the ranges of the output parameter ps by means of the FirstOrder-Second-Moment-Method and compare them with the experimental results. The model with the best fit was then scrutinized further: we calculated the sensitivities of the output ps with respect to the five input parameters described above. Here we used Monte Carlo simulation based on the input distributions obtained before. Going beyond the rather crude picture obtained by scatterplots, we computed stronger statistical measures of sensitivity, such as partial correlation coefficients. These indicators are designed so as to remove hidden influences of co-variates. In addition, this method lends itself to a further application of resampling, allowing to determine the statistical significance of the resulting sensitivities. Further, these methods are applicable in numerical models as well – accordingly, we included a Finite Element calculation in our list of models. In short, the goal of the paper is to propose an approach to model choice and model assessment with the aid of a typical application in geotechnics. Experiments play a twofold role here. On the one hand, 1g-model tests are performed to investigate the behavior of a tunnel face close to failure. On the other hand, the outcome of these tests are contrasted with the predictions of theoretical models. These theoretical models contain material parameters that in turn are determined from (different) experiments. Thus both the outcome of the 1g-model tests and the predictions are uncertain. In the presence of this twofold uncertainty, the assessment of model quality requires sophisticated methods from data analysis and uncertainty analysis. Bootstrap resampling techniques are used to assess the statistical distributions of the input parameters, resulting in variability intervals for the model predictions that can be compared with confidence intervals of the test results. This enables a comparison of the range of the predicted output with the range of the measured output and thus allows to assess the model quality. We take one further step, once the fittest model has been chosen. This step is sensitivity analysis, determining a ranking of the input parameters according to their influence on the model output. The significance of the ranks is again assessed by bootstrap methods. Highly influential input parameters should be known more precisely than less influential parameters. This aids in deciding where to focus the effort in further experimental, in situ or laboratory investigations. In addition, if the influence of an input parameter is classified as non-zero, this supports the structure of the model that considers it as a factor to be accounted for. The paper is organized as follows. In Section 2, we briefly present the theoretical models under investigation. In Section 3, we describe the experimental set-up. Section 4 is devoted to uncertainty and sensitivity analysis. We formulate the statistical data models and explain how joint distributions were obtained by resampling. Then we do the FOSM calculation that leads to an overall assessment of the predictive power of the models. This section concludes with the sensitivity analysis based on Monte Carlo simulations and the statistical indicators mentioned above. In Section 5, we discuss the Finite Element model and the corresponding uncertainty/sensitivity analysis. The final section summarizes our conclusions. The methods of uncertainty/sensitivity analysis are based on our earlier paper [4]; for a general survey of sampling based sensitivity analysis we recommend [8].

for tunnels in soft ground. We use five models for our comparison, namely those of Horn, Kolymbas, Krause, Leca and Dormieux, Vermeer and Ruse. Horn was the first to present a kinematic mechanism with a sliding wedge for the given problem (Fig. 1). The silo theory [43] is used to calculate the vertical force acting on the top of the wedge. Force equilibrium yields the support force as function of the wedge geometry. The necessary support force is the maximum value of this function. We use the original version of the Horn model [15] in our study. Note that there are various variations of the original model [16–19,44,20,22,45,22]. They differ by assumptions about the lateral earth pressure coefficients used in the silo theory, the vertical distribution of vertical earth pressure, cohesional and frictional forces on the top of the wedge, the approximation of the non-rectangular tunnel cross section and the shape of the basal boundary of the sliding wedge. An overview of these variations can be found in [33]. Kolymbas derived a single equation for the support pressure, which follows from the equilibrium condition by assuming an admissible stress field above the tunnel and full mobilization of soil strength at the tunnel crown [29]. Krause considered a single hemispherical body of soil as threedimensional failure mechanism. Cohesion and frictional forces in the contact area build up a resistance against a rotational type of failure. By equilibrating the moments acting on the body, he derived an expression for the necessary support pressure ps [22]. Leca and Dormieux [23] presented a sliding wedge mechanism of one or two solid conical wedges with circular cross sections. Making use of the upper bound theorem, they calculated the support pressure as a function of the inclination angle a of the cone axis with respect to the horizontal. By maximizing the resulting pressure over a range of a, expressions for the necessary support pressure were obtained. Vermeer and Ruse presented a simple equation [46–48], which was developed using Finite Element calculations. The soil was modelled with a linear elastic, perfectly plastic constitutive model with a Mohr–Coulomb failure condition. Failure was triggered by a reduction of the pressure on the tunnel face. The principal input parameters of the above models are the unit weight c, the cohesion c and friction angle u of the soil, as well as the overburden C and the diameter D of the tunnel. The cohesion is set to zero in this investigation, as dry sand was used in the laboratory experiments [33]. The geometry parameters C and D can be determined with high accuracy, do not change throughout the test and are independent of the other parameters.

wid

th w

gth

l

len

cover C cuboid

D prismatic wedge diameter D

2. Theoretical models A lot of researchers have put forward theoretical models and empirical relations to predict the necessary support pressure ps

Fig. 1. Geometry of Horn’s rigid block mechanism.

404

W. Fellin et al. / Structural Safety 32 (2010) 402–410

The dry unit weight cd = cs/(1 + e) is a function of the unit weight of the solids cs and the actual void ratio e. The friction angle u depends generally on the stress level and the density of the soil. While the stress level in the small scale experiments does not change significantly with the depth we assume the friction angle to be a function of the relative density only u = u(Id). The relative density is

Id ¼

el  e el  ed

carriage

goniometer

ð1Þ

wherein the void ratios el (loose state) and ed (dense state) follow from standardized laboratory experiments. We applied a linear relationship u  b0 + b1Id. In summary, we use the d = 5 parameters e, el, ed, b0 and b1 as input for the theoretical models.

load cell knob Fig. 3. Carriage construction with load cell, goniometer and turning knob.

3. Experimental investigation In a research project at the Unit of Geotechnical and Tunnel Engineering 1g-model tests were developed to investigate the behavior of a tunnel face close to failure [33,42]. We note that scaling is always a problem in geotechnical experiments. For the face stability experiments the main issues are: similarity and deterministic size effect (grain size in relation to dimension of the experiment [49,50] and soil non-linearity [51]) and stochastic size effects (decreasing shear resistance with increasing sample size [52,53]). A detailed discussion of these issues is given in [33,42]. However, mechanically sound face stability models should perform well at all scales as long as appropriate input parameters are used. We put effort into the choice of appropriate parameter values for a stress level comparable to the one prevailing in the small scale experiment. Obviously, for large scale experiments other input parameters must be used. The experimental set-up we used is described in the following subsection. 3.1. Model box The model box (Fig. 2) consisted of a steel frame with inner dimensions 37.2  37.2  41.0 (width  depth  height in centimeter) and wooden walls. The tunnel was modelled by a 4 mm thick hollow aluminium cylinder with an inner diameter of Di = 10 cm, which reached approximately 7 cm into the soil domain. The face of the model tunnel was supported by an aluminium disk with a slightly smaller diameter than the inner diameter of the tunnel to exclude friction between disk and tunnel. The piston rod was supported by a linear roller bearing, embedded in the side wall. It allowed for a horizon-

Table 1 Properties of the applied Ottendorf–Okrilla sand [54]. Property

Sand S1

Mean grain size Coeff. of uniformity Grain shape Particle density

d50 U

qs

tal displacement of the piston into (and out of) the cylinder, which is denoted by s. The support force was measured by a load cell at the outside end of the piston rod (Fig. 3) while slowly retracting the piston into the cylinder. 3.2. Sand properties Commercially available Ottendorf–Okrilla sand with grain diameters between 0.1 and 2.0 mm was used. Some properties are listed in Table 1. 3.3. Preparation of the sandbox experiments The tests were performed with dry sand for various cover-todiameter ratios (C/D-ratios) and different initial relative densities Id. Dense samples were prepared by pouring sand into the box with a funnel. With a drop height of approximately 10 cm, layers of 5 cm thickness were installed with an approximately homogeneous

steel frame

model tunnel

model tunnel

0.58 mm 3.6 (sub)angular 2.635 g/cm3

piston

Fig. 2. Box for 1g-model tests.

405

W. Fellin et al. / Structural Safety 32 (2010) 402–410

density distribution. These layers were then compacted by means of a tamper until a high density was reached. Loose samples were prepared by carefully putting the sand into the box with a small shovel, trying to prevent any compacting action. The void ratio was determined from the total mass of the sand inside the box and the occupied volume. The resulting void ratios e and relative densities Id are therefore average values for the material in the box. Spatial variability of void ratio and therefore friction angle can have an impact on the failure mechanism and, hence, face pressure. However, this is currently beyond the ability of the experimental set-up. A first attempt of a kinematic model which considers spatial variability of the friction angle can be found in [27]. 3.4. Results of the sandbox experiments For each combination of C/D and relative density Id two separate tests were performed to check the reproducibility of the results, which were in total 31 tests. Load–displacement curves for loose specimen and different C/D-ratios can be seen in Fig. 4. The resulting force on the piston was normalized with the area of the piston A, the dry unit weight cd and the piston diameter D to F/(AcdD) = p/ (cdD). Fig. 4 shows that all curves (for different overburden) achieve a constant (residual) value after relative displacements of 2–3%. This value is interpreted as (dimensionless) necessary support pressure ND, because application of a lower pressure (in a pressure-controlled test) would lead to infinite displacements. As ND could not be determined precisely, the results must be understood as mean values with a precision of ±0.01. Table 2 compiles results of experiments with different C/D ratios. The results convey that the cover-to-diameter ratio does not significantly influence the results. This is supported by experimental evidence [23] and numerical results [46]. Arching has been ob-

3.5. Determination of the friction angle u In soil mechanics, it is an accepted fact that u depends on stress level and initial density of the soil. Standard soil mechanics tests are not suited for pressure levels below, say, 50 kPa. In the performed experiments stresses around 2 kPa prevailed. Therefore, a series of 75 sand chute experiments were performed with the applied sand (Fig. 5). The chute was slowly tilted until the material slid down. The sand chute has a width of 11.7 cm. This is large compared to the grain size of the sand (0.1–2 mm). Therefore, we do assume that the side walls do not affect the results. All experiments showed an initial slip surface inside the sand layer. So, the influence of the rough chute bottom should also be negligible. Therefore, we assumed that an infinite slope model is applicable and interpreted the tilt angle a at the initiation of the sand avalanche as the peak angle upeak of the sand at very low stress level. We used a simplified method for the estimation of the criticalstate friction angle. The angle of repose of a loose tip of dry soil subjected to toe excavation can serve as an approximation of the critical-state angle uc [55,56]. Herle [57] suggested to slowly pour a pile of sand from a funnel and measure the slope of the pile. This method is standard for the determination of uc for the hypoplastic constitutive model of sand [58] and is used in the following. From the average angle of repose we estimated uc = 32.5° (from 44 separate measurements, standard deviation 1.1°) for the applied sand, which fits into the typical range for quartz sand uc = 33° ± 1° [55]. 4. Uncertainty/sensitivity analysis 4.1. Input distributions

0.40 C/D=0.5, Id=0.34

0.35

A linear relationship between u and Id was suggested by [55,59– 61]. We thus fitted the linear model

C/D=0.5, I =0.25 d

C/D=1.0, Id=0.33

0.30

u ¼ b0 þ b1 Id þ e; Efeg ¼ 0; covfeg ¼ r2 I

C/D=1.0, Id=0.33

p/(γd D)

served in a Particle Image Velocimetry (PIV) analysis of the failure mechanism [33,42]. The arch redirects the weight above the sliding wedge onto the tunnel lining and the surrounding ground in such a way that the soil in front of the piston face builds up less pressure on the face.

ð2Þ

0.25

where E{e} is the expectation value of the error and cov{e} its covariance matrix. The regression model was fitted to the sample (ui, Id,i), i = 1, . . . , n of size n = 75 obtained by the sand chute and

0.20 0.15

ND

0.10 0.05 0.00 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

s/D (%) Fig. 4. Evolution of the normalized support pressure ND = p/(cdD) with the normalized piston advance s/D.

Table 2 Results for the necessary dimensionless support pressure. Test no. 1–4 5–13 14–17 18–23 24–27 28–31

C/D 0.25 0.50 0.75 1.00 1.50 2.00

Material

ND

S1 S1 S1 S1 S1 S1

0.11 0.11 0.12 0.12 0.12 0.13

α

Fig. 5. Sand chute experiments: tilt angle a at the onset of the sand sliding.

W. Fellin et al. / Structural Safety 32 (2010) 402–410

 T   di ¼ ðb0 ; b1 Þi  b R1 ðb0 ; b1 Þi  b b

ð3Þ

which is well known to follow a v2(2)-distribution given normality [66]. In our case this can be consolidated either quantitatively (passing a KS-test) or graphically via quantile–quantile-plots: if di is a sample of a v2(2)-distribution, then the sorted values d[i] are   estimators of the quantiles v2 ð2Þ1 i0:5 , i.e. the corresponding B points should lie on a straight line (see Fig. 7, lower panel). From the viewpoint of stochastic simulation we hence are able to realize ðb0 ; b1 Þ  Nðb; Rb Þ by employing the usual direct scheme

42 40

ϕ [°]

38 36 34

9

β

1

8 7 6 32

33

33.5

20

p−value KS: 0.568 10

0

0

5

10

15

20

d[i] Fig. 7. Empirical joint distribution of (b0, b1) and multivariate normality test.

[67], i.e. generation of an uncorrelated bivariate Gaussian, multiplication by the Cholesky factor of Rb and addition of b. 4.2. FOSM How do different models propagate the uncertainties of the input parameters? To answer this question, we employed the classical First-Order-Second-Moment-Method, thereby obtaining an approximation for mean and variance of the predicted dimensionless support pressure. More specifically, let f be a sufficiently smooth model function depending on random variables Xi (i = 1, . . . , d) with corresponding nominal (mean) values li and standard deviations ri. By linearizing f around l = (l1, l2, . . . , ld), i.e.

f  f ðlÞ þ

 d X  @f    X i  li @X i l i¼1

ð4Þ

and exploiting the properties of mean and variance, we arrive at lf  f(l) and 2 f

r 

d X i¼1

 !2  d X i1 X @f  @f @f  2 r þ 2 rij ; i @X i l @X i @X j l i¼1 j¼1

ð5Þ

where rij denotes possible covariances between the input variables. Here, the required partial derivatives were determined by numerically differentiating the model function (finite difference approximation). For the given task, the three void ratios were assigned independent, univariate normal distributions. As mentioned above, the void ratio e was determined with the total mass of sand and its occupied volume in the sandbox to be e = 0.65. The mean values and standard deviations for the regression parameters b0 and b1 were extracted from the estimated bivariate distribution. The mean values and standard deviations of the limit void ratios el and ed are based on 20 separate laboratory experiments for each limit. We assume that the variation of e is similar to that of el. All values are summarized in Table 3. As b0 and b1 are correlated, also the covariance rb0 b1 and correlation coefficient qb0 b1 were Table 3 Statistical parameters for the five input parameters.

32 30

32.5

β0 2

angle of repose experiments outlined above. Here, it was assumed that the critical void ratio ecrit at a very low stress level (as in our small scale tests) can be approximated by the loose state void ratio el [57,58,62–65]. For e < ecrit (Id > 0) the friction angle u is a peak friction angle upeak; for e = ecrit (Id = 0) peak and critical friction angle are assumed to be equal, thus upeak = ucrit = u. In Fig. 6, the measured friction angles are plotted against the relative density Id. The aforementioned scheme leads to the estimators b0 = 32.72, b1 = 7.26. Testing for the significance of regression by means of the standard F-statistic leads to a clear rejection of the null hypothesis H0:b1 = 0. Taking into account the high coefficient of determination (R2 = 0.8) we hence have strong indications that the linear model is appropriate. We turn to assessing the distribution of the regression parameters b0, b1. Regressing on the data obtained from the sand chute experiments supplies us with a single realization of b0, b1. To obtain a sample of the values (b0, b1) and thus have a possibility to assess the bivariate distribution of (b0, b1) we employed a resampling technique called bootstrapping cases [13]. For background on resampling in general we refer to [14]. More specifically, B = 1000 bootstrap samples of size n = 75 each are generated by the following procedure: from the (n  2) data matrix (ui, Id,i) we randomly draw n-times, so that each row has equal probability of being drawn. The results are combined to produce a bootstrap (n  2) matrix (note that some rows of the bootstrap matrix may be repetitions of rows of the original data matrix). This is repeated B = 1000 times. We perform a linear regression on each of these matrices and obtain a population of n = 1000 pairs of regression coefficients (b0, b1). This population is an approximation to the empirical joint distribution of (b0, b1) and depicted in the upper panel of Fig. 7. From the given scatterplot as well as from the form of the corresponding marginal distributions we expect a bivariate Gaussian distribution ðb0 ; b1 Þ  Nðb; Rb Þ, where b and Rb are the sample mean, respectively the sample covariance of (b0, b1)i, i = 1, . . . , B. This assumption might be tested by examining the set

quantiles χ

406

0

0.2

0.4

0.6

0.8

I

d

Fig. 6. Sand chute results and linear regression.

1

Parameter

Mean li

Std. dev. ri

Distribution type

ed el e b0 b1

0.48 0.71 0.65 32.7 7.28

0.03 0.03 0.05 0.16 0.4

Normal Normal Normal Normal Normal

407

W. Fellin et al. / Structural Safety 32 (2010) 402–410

0.20

0.16

0.16

0.14

0.14

0.12

0.12

0.1

0.1

N

D

0.15

0.10

0.05

0.4

0.45

0.5

0.55

0.65

e 0

Horn

Leca Krause

Ruse

0.7

0.75

0.8

el

d

0.16

0.16

0.14

0.14

0.12

0.12

0.1

0.1

Kolymbas

Fig. 8. Variation of predicted ND for distributed input data; experimental results shaded in grey.

important; from the empirical joint distribution they were quantified as rb0 b1 ¼ 0:0388 and qb0 b1 ¼ 0:6, respectively. Fig. 8 illustrates the result of the FOSM calculations for the different models. The error bars indicate approximate 95% confidence intervals based on the calculated model variance. Apart from the different predictions of the mean lND (central points in Fig. 8), it becomes obvious that the different models process the input uncertainties differently: e.g. the distribution of ND for the Krause model is much narrower than that for the Horn model, although they predict roughly the same nominal value lND . It becomes also clear from Fig. 8 that some of the models fail to predict the obtained experimental results on a 95% confidence level. The best agreement is achieved by the Vermeer–Ruse model, which was investigated further by a Monte Carlo study. Due to the non-linearity of the face pressure models, the best model for our small scale tests may not be the best for large scale tests. However, the decision process would be as outlined here if large scale tests were used as bases for any model comparison.

0.5

0.6

0.7

0.8

32.2

32.7

33.2

β0

e 0.16 0.14 0.12 0.1

6.5

7

7.5

8

β

1

4.3. Monte Carlo study To determine the influence of each parameter on the necessary dimensionless support pressure ND predicted by the respective model, we conducted a Monte Carlo study of size N = 1000, with the common input samples generated according to the postulated probability distributions. Care was taken to ensure that the physical constraints ed,i 6 el,i and ed,i 6 ei hold for each sample. Correlation control was applied in order to keep the sample correlation matrix associated with (ed, el, e) close to diagonal [68]. Generally, the output statistics obtained for each model closely resemble those calculated via FOSM, thus reconfirming our previous ranking. Consequently, results will be presented for the model of Vermeer–Ruse only. We shall stress the fact, however, that all schemes presented below are equally applicable for all models discussed and give very similar qualitative results. A first scan of the influence of the individual input parameters on the output can be done by means of scatterplots (Fig. 9). A decidedly positive correlation between e and ND is visible, a negative correlation between el and ND appears to be present, but is less obvious. However, such plots might be misleading as they cannot take into account hidden interactions between the input parameters under scrutiny. 4.4. Partial correlation coefficients We therefore turn to a method which intends to remove the influence of co-variates on the correlation between a given input

Fig. 9. Scatterplots of five input parameters ed, el, e, b0, b1 versus output ND (Vermeer–Ruse).

variable Xi and the output variable Y. This method is based on the partial correlation coefficient (PCC). In our case, the model output is Y = ND, the five input variables X1, . . . , X5 are ed, el, e, b0 and b1. In order not to overload the notation, we use the generic X and Y for random variables (rather than the specific denotations ed, etc.). The partial correlation between two random variables Xi and Y given a set of co-variates Xni = {X1, . . . , Xi1, Xi+1, . . . , Xd} is defined as the correlation between the two residuals sX i X ni and sYX ni obtained by regressing Xi on Xni and Y on Xni, respectively. More precisely, one first constructs the two regression models

b i ¼ n0 þ X

X

nj X j ;

b ¼g þ Y 0

j–i

X

gj X j ;

ð6Þ

j–i

obtaining the residuals

b i; sX i X ni ¼ X i  X

b: sYX ni ¼ Y  Y

ð7Þ

Since sX i X ni and sYX ni are those parts of Xi and Y that remain after subtraction of the best linear estimates in terms of Xni, the partial correlation coefficient

qXi ;YXni ¼ qðsXi Xni ; sYXni Þ

ð8Þ

quantifies the linear relationship between Xi and Y after removal of any part of the variation due to the linear influence of Xni. As compared to the usual correlation coefficient, the partial correlation

408

W. Fellin et al. / Structural Safety 32 (2010) 402–410

1

0.5

0

−0.5

−1

e

d

e

l

e

β

0

β

1

Fig. 11. Finite Element mesh with 1438 elements (piston not shown).

Fig. 10. PCCs and corresponding 95% bootstrap confidence intervals of five input parameters ed, el, e, b0, b1 with output ND (Vermeer–Ruse).

1

coefficient is a much more discriminating, if not to say radical, measure of influence. In fact, if the input–output map is a truly linear function, the PCC of an input variable that enters with a non-zero coefficient is equal to one. On the other hand, if two input variables are multiples of each other, their PCCs with the output are zero. In reality, input–output maps are not ideally linear functions and so the effect is somewhat moderated. Still the PCCs are an accentuating measure of influence. For further background on PCCs and computational issues see [8,69]. Calculation of the PCCs is based on the original Monte Carlo sample and the associated support pressures. Again, we employed bootstrapping to gain further statistical information on this sensitivity index (Fig. 10). As only support pressure values already available are used in this calculation, empirical confidence intervals for the PCCs can be obtained with a minimum of additional computational cost. Accordingly, all PCCs test to be non-zero and hence can be regarded as assertive. As above, el and e are recognized as the most influential factors. On could suppose that the smallness of the PCCs of b0 and b1 (though they are still significantly different from zero) could originate from the fact that they enter as correlated variables in the model. However, in this case the scatterplots (Fig. 9) confirm the weakness of the correlation with the output ND.

0.5

0

−0.5

−1

d

For reasons of comparison, we subjected a Finite Element model of the sandbox experiments to the same uncertainty and sensitivity analysis. As the FOSM cannot yield any information about the distribution of the predicted ND values, we turned to Monte Carlo simulation performing N = 1000 Finite Element calculations with the same input sample as in the preceding section. For the FE calculations, described in the following, the linear elastic, perfectly plastic Mohr–Coulomb model was used. This model is frequently applied to geotechnical boundary value problems. The elastoplastic Mohr–Coulomb model is implemented in the Finite Element code ABAQUS, which was used for the simulations. It requires five input parameters:  Young’s modulus E and Poisson’s ratio m describe the material behavior in the elastic domain.  The friction angle u, cohesion c and dilation angle w govern the plastic behavior of the material.

e

l

e

β

0

β

1

Fig. 12. PCCs and corresponding 95% bootstrap confidence intervals of five input parameters ed, el, e, b0, b1 with output ND (FE model).

estimated at E = 825 kPa, which accounts for the low stress levels prevailing in the sand specimen. The cohesion c was neglected for dry sand, the dilation angle assumed to be w = 4°. Poisson’s ratio was derived from u with the help of Jaky’s empirical formula, K0  1  sin u, and an expression for K0 from elasticity theory

K0 ¼ 5. Finite Element model

e

m 1m

,

m

1  sin u : 2  sin u

In the model (Fig. 11) soil and piston were modelled as independent parts, which interacted via frictionless contact elements. This way the laboratory conditions were best reproduced. The tunnel lining was not modelled explicitly; instead the nodes on the perimeter of the tunnel were fixed in all directions, thus representing a rigid lining. During the numerical simulation the piston was retracted into the tunnel and the resultant force on the piston was logged. The ND values of the Finite Element calculations were derived from the support force after 0.5 mm piston displacement, when all the curves had reached a well-defined plateau. As before, we calculated the partial correlation coefficients between input and the necessary dimensionless support pressure estimated by the FE model (Fig. 12). As can readily be seen, the FE model and the Vermeer–Ruse model exhibit very similar sensitivities with respect to the uncertain soil parameters. 6. Conclusions

The self-weight of the soil cd and the friction angle u were calculated from the input sample exactly in the same way as described in the previous section. The Young modulus was

In this paper, we demonstrated how statistical estimators in conjunction with laboratory experiments can be used to assess

W. Fellin et al. / Structural Safety 32 (2010) 402–410

the predictive power and robustness of models in geotechnical engineering. We performed a case study with a number of theoretical and numerical models for the required support pressure for shield tunnelling. Ultimately, this may serve as a basis for determining the safety of the tunnel face during construction and for gaining deeper understanding of the relevant mechanisms leading to failure. In particular, we showed that novel statistical methods constitute a powerful tool for identifying the most influential model parameters, and we showed that a comparison of models is possible in terms of analyzing the uncertainty of the model answer and contrasting it with test results, which also contain uncertainties. It turned out that bootstrapping is a suitable method for producing samples of multi-dimensional random variables and for assessing the significance of the computed sensitivities. This is particularly convenient in a Finite Element analysis, because it does not require additional computational cost beyond the original Monte Carlo simulation. We remark that the sensitivities obtained in the Vermeer–Ruse model and the FE-model turned out to be quite similar. This may be explained by the fact that the Vermeer–Ruse model was originally calibrated on a similar FE-model. Acknowledgements The first two authors (W. F., A. K.) gratefully acknowledge support by the Tyrolean Science Foundation, the third author (J. K.) is recipient of a DOC-fellowship of the Austrian Academy of Sciences at the Breath Research Unit. References [1] Bucher C. Robustness analysis in structural optimization. Struct Infrastruct Eng 2009;5(4):287–93. [2] Konietzky H, Schlegel R. Using sensitivity analysis and optimization for civil engineering and geotechnical applications. In: Proceedings optimization and stochastic days 4.0, Weimar; 2007. . [3] Oberguggenberger M, Fellin W. Reliability bounds through random sets: nonparametric methods and geotechnical applications. Comput Struct 2008;86:1093–101. [4] Oberguggenberger M, King J, Schmelzer B. Classical and imprecise probability methods for sensitivity analysis in engineering: a case study. Int J Approx Reason 2009;50:680–93. [5] Archer G, Saltelli A, Sobol I. Sensitivity measures, ANOVA-like techniques and the use of the bootstrap. J Stat Comput Simulat 1997;58:99–120. [6] Ferson S, Tucker W. Sensitivity analysis using probability bounding. Reliab Eng Syst Safety 2006;91:1435–42. [7] Hall J. Uncertainty-based sensitivity indices for imprecise probability distributions. Reliab Eng Syst Safety 2006;91:1443–51. [8] Helton JC, Johnson JD, Sallaberry CJ, Storlie CB. Survey of sampling-based methods for uncertainty and sensitivity analysis. Reliab Eng Syst Safety 2006;91:1120–75. [9] Kreinovich V, Beck J, Ferregut C, Sanchez A, Keller G, Averill M, et al. MonteCarlo-type techniques for processing interval uncertainty, and their potential engineering applications. Reliab Comput 2007;13:25–69. [10] Sobol I. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math Comput Simulat 2001;55:271–80. [11] Good P. Resampling methods: a practical guide to data analysis. 2nd ed. Boston (MA): Birkhäuser; 2001. [12] Lunneborg C. Data analysis by resampling: concepts and applications. Pacific Grove, CA: Duxbury Press/Thomson Learning; 2000. [13] Montgomery D, Peck E, Vining G. Introduction to linear regression analysis. New York: Wiley; 2001. [14] Shao J, Tu D-S. The Jackknife and the Bootstrap. New York: Springer-Verlag; 1995. [15] Horn M. Horizontaler Erddruck auf senkrechte Abschlussflächen von Tunneln, Landeskonferenz der ungarischen Tiefbauindustrie (Deutsche Überarbeitung STUVA), Düsseldorf; 1961. [16] Anagnostou G, Kovári K, Ein Beitrag zur Statik der Ortsbrust beim Hydroschildvortrieb. In: Symposium ’92, Probleme bei maschinellen Tunnelvortrieben?, Gerätehersteller und Anwender berichten.; 1992. [17] Anagnostou G, Kovári K. Face stability conditions with earth-pressurebalanced shields. Tunn Undergr Space Technol 1996;11(2):165–73. [18] Girmscheid G. Tunnelbohrmaschinen-Vortriebsmethoden und Logistik. In: Betonkalender. Berlin: Ernst & Sohn; 2005. p. 119–256 [chapter 1.3].

409

[19] Mayer P-M, Hartwig U, Schwab C. Standsicherheitsuntersuchungen der Ortsbrust mittels Bruchkörpermodell und FEM. Bautechnik 2003;80:452–67. [20] Holzhäuser J. Problematik der Standsicherheit der Ortsbrust beim TBMVortrieb im Betriebszustand Druckluftstützung. In: Beiträge anlässlich des 50. Geburtstages von Herrn Professor Dr.-Ing. Rolf Katzenbach, no. 52 in Mitteilungen des Institutes und der Versuchsanstalt für Geotechnik, Darmstadt University of Technology; 2000. p. 49–62. [21] Vermeer P, Ruse N, Dong Z, Härle D. Ortsbruststabilität von Tunnelbauwerken am Beispiel des Rennsteig-Tunnels. In: Tagungsband 2. TAE Kolloquium ‘‘Bauen in Boden und Fels”; 2000. p. 195–202. [22] Krause T. Schildvortrieb mit flüssigkeits-und erdgestützter Ortsbrust, Mitteilung des Instituts für Grundbau und Bodenmechanik, Technische Universität Braunschweig, no. 24; 1987. [23] Leca E, Dormieux L. Upper and lower bound solutions for the face stability of shallow circular tunnels in frictional material. Géotechnique 1990;40(4):581–606. [24] Soubra A-H. Three-dimensional face stability analysis of shallow circular tunnels. In: International conference on geotechnical and geological engineering, Melbourne, Australia; November 19–24, 2000. p. 1–6. [25] Soubra A-H. Kinematical approach to the face stability analysis of shallow circular tunnels. In: 8th International symposium on plasticity, British Columbia, Canada; 2000. p. 443–5. [26] Soubra A-H, Dias D, Emeriault F, Kastner R. Three-dimensional face stability analysis of circular tunnels by a kinematical approach. In: GeoCongress 2008: characterization, monitoring, and modeling of geosystems (GSP 179), New Orleans; 2008. p. 894–901. [27] Mollon G, Phoon K, Dias D, Soubra A-H. A new 2D failure mechanism for face stability analysis of a pressurized tunnel in spatially variable sands. In: GeoFlorida 2010: advances in analysis, modeling & design (GSP 199), ASCE, Reston; 2010. p. 2052–61. [28] Atkinson JH, Potts DM. Stability of a shallow circular tunnel in cohesionless soil. Géotechnique 1977;27(2):203–15. [29] Kolymbas D. Tunnelling and tunnel mechanics. Berlin: Springer-Verlag; 2005. [30] Balthaus H. Standsicherheit der flüssigkeitsgestützten Ortsbrust bei schildvorgetriebenen Tunneln. In: Scheer J, Ahrens H, Bargstädt HJ, editors. Festschrift Heinz Duddeck. Institut für Statik der Technischen Universität Braunschweig. Berlin: Springer; 1988. p. 477–92. [31] Takano D, Otani J, Fukushige S, Natagani H. Investigation of interaction behavior between soil and face bolts using X-ray CT. In: Desrues J, Viggiani G, Bèsuelle P, editors. Advances in X-ray tomography for geomaterials. London: ISTE Ltd.; 2006. p. 389–95. [32] Sterpi D, Cividini A. A physical and numerical investigation on the stability of shallow tunnels in strain softening media. Rock Mech Rock Eng 2004;37(4):277–98. [33] Kirsch A. On the face stability of shallow tunnels in sand. In: Kolymbas D, editor. Advances in geotechnical engineering and tunneling. Berlin: Logos Verlag; 2009. [34] Subrin D, Branque D, Berthoz N, Wong H. Kinematic 3D approaches to evaluate TBM face stability: comparison with experimental laboratory observations. In: EURO:TUN 2009 – computational methods in tunnelling; 2009. p. 801–7. [35] Chambon P, Corté J-F, Garnier J, König D. Face stability of shallow tunnels in granular soils. In: Ko H-Y, McLean F, editors. Centrifuge 91. Rotterdam: Balkema; 1991. p. 99–105. [36] Chambon P, Corté J-F. Shallow tunnels in cohesionless soil: stability of tunnel face. ASCE J Geotech Eng 1994;120(7):1148–65. [37] Al Hallak R, Garnier J, Léca E. Experimental study of the stability of a tunnel face reinforced by bolts. In: Kusakabe O, Fujita K, Miyazaki Y, editors. Geotechnical aspects of underground construction in soft ground. Rotterdam: Balkema; 2000. p. 65–8. [38] Plekkenpol JW, van der Schrier JS, Hergarden HJ. Shield tunnelling in saturated sand – face support pressure and soil deformations. In: Bezuijen A, van Lottum H, editors. Tunnelling: a decade of progress, geodelft 1995– 2005. London: Taylor & Francis; 2006. [39] Kamata H, Mashimo H. Centrifuge model test of tunnel face reinforcement by bolting. Tunnell Undergr Space Technol 2003;18:205–12. [40] Kimura T, Mair RJ. Centrifugal testing of model tunnels in soft soil. In: Proceedings of the 10th international conference on soil mechanics and foundation engineering, Stockholm, vol. 1; 1981. p. 319–22. [41] Renpeng C, Jun L, Xuecheng B, Yunmin C. Large scale experimental investigation on face stability of shallow tunnels in dry cohesionless soil. In: EURO:TUN 2009 – computational methods in tunnelling; 2009. p. 809–16. [42] Kirsch A. Experimental investigation of the face stability of shallow tunnels in sand. Acta Geotech 2010;5(1):43–62. [43] Janssen HA. Versuche über Getreidedruck in Silozellen. Zeitschrift des Vereins Deutscher Ingenieure 1895;39(35):1045–9. [44] Kirsch A, Kolymbas D. Theoretische Untersuchung zur Ortsbruststabilität. Bautechnik 2005;82(7):449–56. [45] Mohkam M, Wong YW. Three dimensional stability analysis of the tunnel face under fluid pressure. In: Swoboda G. (Ed.), Proceedings of the 6th international conference on numerical methods in geomechanics, Innsbruck, vol. 4; 1988. p. 2271–8. [46] Ruse NM. Räumliche Betrachtung der Standsicherheit der Ortsbrust beim Tunnelvortrieb, Mitteilungen des Instituts für Geotechnik, Universität Stuttgart, no. 51; 2004. [47] Vermeer P, Ruse N. Die Stabilität der Tunnelortsbrust in homogenem Baugrund. Geotechnik 2001;24(3):186–93.

410

W. Fellin et al. / Structural Safety 32 (2010) 402–410

[48] Vermeer P, Ruse N, Marcher T. Tunnel heading stability in drained ground. Felsbau 2002;20(6):8–18. [49] Technical Commitee 2 of ISSMGE – Physical Modelling in Geotechnics, catalogue of scaling laws and similitude questions in centrifuge modelling; 2007. [50] Chambon P, Couillaud A, Munch P, Schürmann A, König D. Stabilité du front de taille d’un tunnel: Étude de l’effet d’échelle. In: Geo 95; 1995. p. 3 [cited in [49]]. [51] Muir Wood D. Geotechnical modelling. London: Spon Press; 2004. [52] Tatsuoka F, Goto S, Tanaka T, Tani K, Kimura Y. Particle size effects on bearing capacity of footing on granular material. In: Asaoka A, Adachi T, editors. Deformation and progressive failure in geomechanics. New York: Pergamon; 1997. p. 133–8. [53] Tejchman J. FE-simulations of a direct wall shear box test. Soils Found 2004;44(4):67–81. [54] Laudahn A. An approach to 1g modelling in geotechnical engineering with soiltron. Advances in geotechnical engineering and tunnelling, vol. 11. Berlin: Logos; 2004. [55] Bolton M. The strength and dilatancy of sands. Géotechnique 1986;36(1):65–78. [56] Cornforth DH. Prediction of drained strength of sands from relative density measurements. In: Evaluation of relative density and its role in geotechnical projects involving cohesionless soils. In: ASTM special technical publication, vol. 523. American Society for Testing and Materials, Philadelphia; 1973. p. 281–303. [57] Herle I. Hypoplastizität und Granulometrie einfacher Korngerüste, Veröffentlichungen des Institutes für Bodenmechanik und Felsmechanik der Universität Fridericana Karlsruhe, no. 142; 1997. [58] Herle I, Gudehus G. Determination of parameters of a hypoplastic constitutive model from properties of grain assemblies. Mech Cohes–Frict Mater 1999;4:461–86.

[59] Fukushima S, Tatsuoka F. Strength and deformation characteristics of saturated sand at extremely low pressures. Soils Found 1984;24(4):30–48. [60] De Beer E. Influence of the mean normal stress on the shearing strength of sand. In: Proceedings of the 6th international conference soil mechanics and foundation engineering, Montreal, vol. 1; 1965. p. 165–9. [61] Wu W. Hypoplastizität als mathematisches Modell zum mechanischen Verhalten granularer Stoffe, Veröffentlichungen des Institutes für Bodenmechanik und Felsmechanik der Universität Fridericana Karlsruhe,no. 129; 1992. [62] Alarcon-Guzman A, Leonards G, Chameau J. Undrained monotonic and cyclic strength of sands. J Geotech Eng 1988;144(10):1089–109. [63] Ishihara K. Liquefaction and flow failure during earthquake. Géotechnique 1993;43(3):351–415. [64] Mohamad R, Dobry R. Undrained monotonic and cyclic triaxial strength of sand. J Geotech Eng 1986;112(10):941–58. [65] Verdugo R, Ishihara K. The steady state of sandy soils. Soils Found 1996;36(2):81–91. [66] Anderson TW. An introduction to multivariate statistical analysis. New York: Wiley; 1958. [67] Devroye L. Non-uniform random variate generation. Berlin: Springer-Verlag; 1986. [68] Iman R, Conover W. A distribution-free approach to inducing rank correlation among input variables. Commun Stat-Simul Comput 1982;11:311–34. [69] Iman R, Shortencarier M, Johnson J, A FORTRAN 77 program and user’s guide for the calculation of partial correlation and standardized regression coefficients. SANDIA Report SAND85-0044, Tech. rep., Sandia National Laboratories, Albuquerque; 1985.