Efficient response surface method for practical geotechnical reliability analysis

Computers and Geotechnics 69 (2015) 496–505

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

Efficient response surface method for practical geotechnical reliability analysis J. Zhang a, H.Z. Chen a, H.W. Huang a, Z. Luo b,⇑ a Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education and Department of Geotechnical Engineering, Tongji University, 1239 Siping Road, Shanghai 200092, China b Department of Civil Engineering, University of Akron, Akron, OH 44325, USA

a r t i c l e

i n f o

Article history: Received 2 February 2015 Received in revised form 10 June 2015 Accepted 11 June 2015

Keywords: Reliability Response surface method Shallow foundation Deep mixed columns

a b s t r a c t Although numerical models have been widely used in the geotechnical profession, their applications for reliability analysis are still rather limited mainly because most geotechnical numerical programs lack the probabilistic function. In this study, an efficient response surface method is suggested for geotechnical reliability analysis using existing numerical programs. The developed approach can effectively avoid the occurrence of negative values for positive random parameters, and thus it solves the key limitation of the classical response surface method in geotechnical applications. To facilitate its application, a procedure is designed for automating reliability analysis with the commercially available program, FLAC3D. The developed method and procedure are demonstrated in detail using a slope example. The versatility of the suggested procedure is illustrated through the serviceability reliability analysis of a soft ground improved with deep mixing columns. For the reinforced ground studied in this paper, the reliability is mainly governed by the uncertainties in Young’s moduli of the soil and reinforced columns. The suggested method can be conveniently used for reliability-based design of deep mixing columns. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction To assess the effect of uncertainties on geotechnical predictions, reliability methods are increasingly being adopted in geotechnical engineering (e.g., Refs. [1,2]). Several geotechnical design codes have been developed based on the reliability theory (e.g., Refs. [3–5]). Probably because most geotechnical numerical programs do not have the probabilistic analysis function, the applications of reliability methods are still limited to problems with relatively simple limit state functions. Compared with the wide application of numerical models in geotechnical engineering, the lack of ability for reliability analysis based on numerical models greatly limits the value of reliability analysis in geotechnical engineering. How to implement reliability analysis utilising sophisticated geotechnical numerical models has been one of the key challenges for the application of reliability methods in geotechnical engineering. Efforts have been dedicated to the development of specific numerical programs capable of geotechnical reliability analysis (e.g., Refs. [6,7]). Such programs can profoundly facilitate the applications of reliability methods for complex geotechnical problems. ⇑ Corresponding author. E-mail addresses: [email protected] (J. Zhang), [email protected] (H.Z. Chen), [email protected] (H.W. Huang), [email protected] (Z. Luo). http://dx.doi.org/10.1016/j.compgeo.2015.06.010 0266-352X/Ó 2015 Elsevier Ltd. All rights reserved.

However, developing such numerical programs requires expertise in both geotechnical numerical analysis and geotechnical reliability. In this regard, response surface methods (RSM) as alternatives are adopted for complex geotechnical reliability problems, in which the reliability analysis is realised through the computationally efficient models that approximate the deterministic numerical solutions (e.g., Refs. [8–15]). Among the RSM available, the iterative method suggested by Bucher and Bourgund [16] based on the first-order reliability method (FORM) is shown to be quite efficient and thus commonly used in various fields (e.g., Refs. [17–22]). The method by Bucher and Bourgund [16] is termed the classical RSM in this paper. As will be shown subsequently in this study, the classical RSM may be subjected to serious convergence difficulties and may not be suitable for geotechnical applications. The objective of this paper is to suggest an efficient and robust RSM for geotechnical reliability analysis that can overcome the convergence difficulty of classical RSM. This method can automate reliability analysis using commercial geotechnical programs such as FLAC3D [23]. As such, the suggested method can profoundly extend the capability of the profession for reliability analysis of complex geotechnical problems. The structure of this paper is as follows: First, FORM is reviewed and the limitation of the classical RSM is discussed. Then, a modified RSM is introduced to deal with

J. Zhang et al. / Computers and Geotechnics 69 (2015) 496–505

the limitation of the classical RSM. Thereafter, a procedure for automating geotechnical reliability analysis using FLAC3D based on the suggested RSM is described. Finally, the applicability and effectiveness of the suggested method are illustrated through the serviceability reliability analysis of a soft ground reinforced with deep mixed (DM) columns.

approximated by Eq. (3), xD computed by Eq. (3) may not be close to the design point of the actual performance function. In such a case, the response surface needs to be updated with a new set of calibration points around a new centre point determined using the following equation:

xc ¼ lx  gðlx Þ 2. First-order reliability method (FORM) Let x denote the uncertain variables in the performance function g(x) with g(x) < 0 indicating failure. Various methods such as FORM (e.g., Ref. [24]) and Monte Carlo simulation (e.g., Ref. [25]) are available for the estimation of failure probability, defined as the probability of g(x) < 0. Although Monte Carlo simulation is one of the most flexible and versatile approaches, its application is often challenged by the considerable computational work involved. For comparison, as FORM is reasonably accurate and also computationally efficient, it is applicable to a large number of practical problems and widely used in practice [26]. In this paper, the focus is on geotechnical reliability analysis using solutions of numerical models, and FORM is employed to achieve computational efficiency. With FORM, the failure probability (pf) is calculated as follows (e.g., Ref. [24]):

pf ¼ 1  UðbÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b ¼ min yR1 yT

ð1Þ ð2Þ

gðxÞ¼0

where U = cumulative distribution function (CDF) of the standard normal variable, b = reliability index, y = reduced variables of x, and R = correlation matrix of the random variables. FORM can be efficiently executed in a spreadsheet when the limit state function has an explicit form [24]. When the explicit solution of a problem is not available and it has to be evaluated with a stand-alone numerical program, the implementation of FORM is challenging due to the coupling between the numerical program and FORM. To deal with such a challenge, Bucher and Bourgund [16] suggested a FORM-based RSM for conducting reliability analysis based on computationally expensive deterministic models, which is now widely employed in various fields. As mentioned previously, the method suggested by Bucher and Bourgund [16] is termed the classical RMS in this study.

The classical RSM approximates the performance function g(x) by a second-order polynomial function:

gðxÞ  b0 þ

k k X X b i xi þ bkþi x2i i¼1

lx  xD gðlx Þ  gðxD Þ

ð3Þ

i¼1

where xi = the ith element of x, k = dimension of x, and bi (i = 0, 1, . . ., 2k) = unknown deterministic coefficients. To determine the (2k + 1) unknown coefficients, the performance function can be first evaluated around a centre point xc = {xc1, xc2, . . . , xck} and other 2k points around xc: {xc1 ± mrx1, xc2, . . . , xck}, {xc1, xc2 ± mrx2, . . . , xck}, . . . , and {xc1, xc2, . . . , xck ± mrxk}, where m is a parameter determining the relative distance of the calibration points and rxi = standard deviation of xi. Equating the values of the performance function with those calculated using Eq. (3) at the prescribed (2k + 1) calibration points, the unknown coefficients can then be solved. After the second-order polynomial function is established, the reliability index can be calculated with Eq. (3) instead of the numerical model. Let xD denote the design point found based on Eq. (3) using FORM. When the performance function is not well

ð4Þ

where lx = mean of x. Through the updating given by Eq. (4), the sampling points are expected to move closer to the limit state function g(x) = 0, indicating a more accurate search for the design point. With the new response surface (Eq. (3) with updated coefficients), the reliability index can be recalculated. This updating procedure is iterated until the resulting FORM reliability index does not change within a tolerable error eb such as eb = 0.01. The above RMS has extensive applications in many fields (e.g., Refs. [17–22]). When determining the calibration points in the classical RSM, however, negative values may occur for positive random variables [21]. Such a phenomenon often occurs when the coefficient of variation (COV) of the random variables is large or when the failure probability is small. In geotechnical engineering, many random variables, such as the cohesion and the friction angle, are non-negative engineering properties. In addition, the uncertainties involved in geotechnical engineering are often higher than those in other fields [27]. Hence, the occurrence of negative values for positive random variables when generating the calibrating points using the classical RSM is quite common in geotechnical engineering, which has previously been noticed by Mollon et al. [21]. Such physically impermissible values result in the iteration termination of numerical programs before the reliability index converges. Therefore, the classical RSM may not be suitable for geotechnical reliability analysis. Although the step size during the iteration process can be adjusted manually to avoid negative values for positive random variables [21], such a procedure has to be realised in a trial-and-error manner and may not always be effective. 4. Modified response surface method Let y denote the reduced variables of x, and let x = T(y) denote the transformation relationship between reduced variables and original variables. Substituting this relationship into g(x), the performance function can be expressed in terms of y as follows:

gðxÞ ¼ g½TðyÞ ¼ GðyÞ

3. Classical RSM and its limitation

497

ð5Þ

where G(y) is the performance function in the reduced space. A set of transformation equations for different types of random variables is documented in the study by Low and Tang [24]. For instance, if xi is a lognormal random variable with a mean of lxi and a COV of dxi, it can be related to its reduced variable yi as follows:

xi ¼ Tðyi Þ ¼ expðki þ ni yi Þ

ð6Þ

where ki and ni are the mean and standard deviation of ln xi, respectively, which can be calculated with the following equations:

ki ¼ ln lxi  0:5n2i qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi ni ¼ ln 1 þ d2xi

ð7Þ ð8Þ

As shown in Eq. (6), xi is always non-negative as assumed in the lognormal distribution regardless of the value of yi. This procedure is also suitable for other types of non-negative random variables. Thus, if the calibration points are determined in the reduced space instead of the original space, the physically impermissible sampling points can be effectively avoided. With this realisation, we suggest that one can approximate the performance function in the reduced space as follows:

498

J. Zhang et al. / Computers and Geotechnics 69 (2015) 496–505

GðyÞ  b0 þ

k 2k X X bi yi þ bi y2i i¼1

ð9Þ

i¼kþ1

where bi (i = 0, 1, 2, . . . , 2k + 1) = coefficients to be determined. In Eq. (9), the coefficients can be determined based on the response of G(y) at the following 2n + 1 points: yc = {yc1, yc2, . . . , ycn}, {yc1 ± mry1, yc2, . . . , ycn}, {yc1, yc2 ± mry2, . . . , ycn}, . . . , and {yc1, yc2, . . . , ycr ± mryn}, where ryi = standard deviation of yi. As the standard deviation of a reduced variable is 1 (i.e., ryi = 1), the 2k calibration points around yc are indeed {yc1 ± m, yc2, . . . , ycn}, {yc1, yc2 ± m, . . . , ycn}, . . . , and {yc1, yc2, . . . , ycn ± m}. The reliability analysis can then be conducted using the approximated performance function Eq. (9). Similar to the classical RSM, if the performance function is not well approximated by a second-order polynomial function, the reliability index obtained based on Eq. (9) may not be accurate. In such a case, the response surface can be updated based on a new set of calibration points with the centre determined as follows:

yc ¼ ly  Gðly Þ

ly  yD Gðly Þ  GðyD Þ

ð10Þ

Note that y is a vector of reduced variables, so the elements of

ly are all zero. Eq. (10) can be further written as follows:

yc ¼

yD Gðly Þ Gðly Þ  GðyD Þ

ð11Þ

Eq. (9) can then be calibrated again using responses of G(y) evaluated at the new calibration points, and the reliability analysis can be conducted based on the updated response surface. Such a process is iterated until the obtained reliability index does not change within a tolerable error eb. This modified RSM needs repetitive execution of the numerical program. When stand-alone deterministic programs are used to solve the geotechnical model, one should first use the deterministic geotechnical program to evaluate the responses of the numerical model at the calibration points, and then fit the response surface and perform reliability analysis with another software such as Excel or Matlab. The above process shall be continued until convergence is achieved. As repetitive switching between the reliability analysis software and the deterministic geotechnical program manually is quite cumbersome and may be subjected to human error, a procedure that can automate the above process is highly desirable. Motivated by the work of Youssef Abdel Massih and Soubra [20], a procedure will be introduced later in this study to implement the modified RSM automatically with FLAC3D. It should be noted that, although FLAC3D is used in this paper, this procedure may also be extended for reliability analysis with other deterministic programs.

5. Reliability analysis of a shallow foundation with explicit performance function Consider a shallow strip foundation subjected to a loading of q = 100 kPa as shown in Fig. 1. Its bearing capacity qu can be calculated using the following equation:

qu ¼ 0:5c1 BN c þ cN c þ c2 Df Nq

where B is the width of the foundation; Df is the depth of the foundation base; c1 is the unit weight of the soil under foundation base; c2 is the unit weight of the soil above the foundation base; c is the cohesion of the soil; and Nc, Nq, and Nc are bearing-capacity factors. Sieffert and Bay-Gress [28] surveyed and compared the bearing capacity factors used in European countries. In this study, the bearing capacity factors suggested by Cherubini [29] are used:

Nc ¼ 1:8ðNq  1Þ tan u p u þ Nq ¼ tan2 ep tan u 4 2 Nc ¼ ðNq  1Þ cot u

ð13Þ ð14Þ ð15Þ

where u is the friction angle of the soil. In this example, c1 = c2 = 18 kN/m3, B = 1.5 m, and Df = 0 m. In this example, c and u are modelled as random variables. Thus, x = {c, u}. The mean and standard deviation of c are 8 and 4 kPa, respectively. The mean and standard deviation of u are 30° and 6°, respectively. Both c and u are assumed to be log-normally distributed and statistically independent. Based on FORM, the reliability index of this shallow foundation is 2.439 (pf = 7.4  103), and the design point is {3.907 kPa, 19.503°} in the original space and {1.281, 2.075} in the reduced space. Based on Monte Carlo simulation with one million samples, the failure probability of the shallow foundation is 5.7  103. The failure probability calculated based on FORM is not exactly the same as that calculated based on Monte Carlo simulation due to the non-linearity in the performance function. The reliability index of this foundation is first analysed with the classical RSM with m = 1, 2, and 3, respectively, which are step sizes commonly used in the literature (e.g., Refs. [16–18]). The values of xc and xD found during the iteration process are shown in Table 1. In the case of m = 1, the value of xci  mrxi of c becomes negative in the second iteration. A similar phenomenon can also be observed when m = 2 and m = 3 are adopted. Thus, the occurrence of a negative value for a positive random variable is observed in all the cases of m = 1, 2, and 3. It was suggested (e.g., Ref. [17]) that the classical RSM may be improved by initially using a large step size and then reducing the step size in the subsequent iterations. To study the effect of such a scheme, the classical RSM is also implemented by adopting m = 2 in the first iteration and m = 1 in the subsequent iterations, and the results are shown in Table 1.

Table 1 Intermediate results from the classical RSM (shallow foundation example, eb = 0.01). Iter. no.

xc

m = 1a

1 2

8.000 30.000 3.275 23.447 2.012 2.635 22.559

m = 2a

1 2

8.000 30.000 6.315 25.548 0.760 5.173 22.529

m = 3a

1

8.000 30.000

m = 2 in the first iteration, and 1 2 m = 1 in subsequent 3 iterations

Fig. 1. A shallow strip foundation subjected to uniform loading.

ð12Þ

c

xD

u

c

b

u

8.000 30.000 6.315 25.548 0.760 5.173 22.529 4.046 19.955 2.302 3.970 19.763

a In the cases of m = 1, m = 2, and m = 3, negative values are generated for c when determining calibration points with xci  mrxi in the second iteration, the second iteration, and the first iteration, respectively.

499

J. Zhang et al. / Computers and Geotechnics 69 (2015) 496–505

In such a case, the value of xci  mrxi of c becomes negative in the third iteration. It seems difficult, if not impossible, to find a universal step size changing scheme that can effectively avoid the occurrence of negative values for positive random variables. In the above analysis, negative values occur when xci is less than mrxi. To avoid the occurrence of such negative values, one may impose a lower bound value for xi if xci is less than mrxi. For instance, if xci  mrxi < 0, a lower bound value of 1 can be adopted for xi and the corresponding calibration point can be changed from {xc1, xc2, . . . , xci + mrxi, xci  mrxi, . . . , xck} to {xc1, xc2, . . . , xci + mrxi, 1, . . . , xck}. The reliability analysis results based on such a scheme are shown in Table 2, where a lower bound value of 1 kPa is adopted for c. We can see that after imposing a lower bound value, convergence can be achieved in all the cases of m = 1, 2, and 3, and the computed reliability indexes are 2.431, 2.463 and 2.542, respectively. The dependence of the converged reliability index on the value of m has also been previously observed in Guan and Melchers [18]. Among the converged reliability indexes obtained, the one calculated based on m = 1 is closest to the actual FORM reliability index. This is probably because a response surface based on closer calibration points around the design point can provide more accurate information about the performance function around the design point, and hence it is more effective for failure probability calculation. To study the effect of the changing step size during iteration, the above RSM is also implemented with m = 2 or m = 3 in the first iteration and m = 1 in the subsequent iterations, and the results are summarised in Table 2. When such step size changing schemes are adopted, the reliability indexes obtained are closer to the actual FORM reliability index. Thus, while changing the step size during iteration may not be able to avoid the occurrence of negative values for positive random variables, it can improve the accuracy of the calculated reliability index. Finally, the above reliability problem is also studied with the method suggested in this paper, and the results obtained based on different m values are summarised in Table 3. Convergent results can be found in all the cases of m = 1, 2, and 3, and the

Table 2 Intermediate results obtained from the classical RSM improved by imposing lower bound values for non-negative random variables (shallow foundation example, eb = 0.01). Iter. no.

xc

m=1

1 2 3

8.000 30.000 3.275 23.447 2.012 2.635 22.559 3.697 19.857 2.427 3.687 19.832 3.987 19.438 2.431

m=2

1 2 3 4 5

8.000 5.173 4.564 4.404 4.372

30.000 22.529 19.686 19.039 18.930

6.315 4.759 4.437 4.377 4.367

25.548 20.271 19.138 18.945 18.916

0.760 2.069 2.395 2.454 2.463

m=3

1 2 3 4 5 6 7

8.000 5.050 5.614 5.162 4.959 4.896 4.878

30.000 23.965 20.292 18.858 18.413 18.293 18.262

7.072 6.048 5.300 4.998 4.906 4.881 4.875

28.102 22.059 19.402 18.561 18.331 18.271 18.257

0.232 1.496 2.196 2.446 2.519 2.537 2.542

m = 2 in the first iteration, and 1 2 m = 1 in subsequent 3 iterations 4

8.000 5.173 3.970 4.007

30.000 22.529 19.763 19.396

6.315 4.046 4.012 4.024

25.548 19.955 19.409 19.357

0.760 2.302 2.431 2.440

m = 3 in the first iteration, and 1 2 m = 1 in subsequent 3 iterations 4 5

8.000 5.050 3.822 3.983 4.020

30.000 23.965 20.074 19.442 19.362

7.072 3.933 3.991 4.021 4.026

28.102 20.336 19.462 19.365 19.351

0.232 2.254 2.425 2.439 2.440

c

xD

u

c

b

u

Table 3 Intermediate results from the suggested RSM (shallow foundation example, eb = 0.01). Iter. no.

yc

m=1

1 2 3 4

0.000 0.000 2.156 1.025 2.483 1.180 1.586 1.921 1.582 1.916 1.278 2.072 1.279 2.073 1.246 2.095

2.387 2.492 2.434 2.438

m=2

1 2 3 4 5 6

0.000 0.290 0.865 1.070 1.118 1.126

0.000 1.286 1.939 2.133 2.168 2.173

0.112 0.780 1.052 1.115 1.126 1.128

0.496 1.749 2.098 2.163 2.173 2.174

0.508 1.915 2.347 2.433 2.447 2.449

m=3

1 2 3 4 5 6 7

0.000 0.041 0.362 0.700 0.864 0.921 0.938

0.000 1.125 1.815 2.140 2.256 2.292 2.302

0.003 0.271 0.648 0.846 0.916 0.937 0.943

0.088 1.362 1.983 2.209 2.278 2.298 2.304

0.088 1.389 2.087 2.365 2.455 2.482 2.489

m = 2 in the first iteration, and m = 1 in subsequent iteration

1 2 3 4 5

0.000 0.290 1.200 1.246 1.243

0.000 1.286 1.983 2.088 2.098

0.112 1.162 1.244 1.243 1.242

0.496 1.920 2.085 2.098 2.099

0.508 2.244 2.428 2.439 2.439

m = 3 in the first iteration, and m = 1 in subsequent iterations

1 2 3 4 5

0.000 0.041 1.182 1.247 1.244

0.000 1.125 1.950 2.085 2.098

0.003 1.131 1.244 1.243 1.243

0.088 1.867 2.080 2.098 2.099

0.088 2.183 2.424 2.438 2.439

c

yD

u

c

b

u

corresponding computed reliability indexes are 2.438, 2.449, and 2.489, respectively. Similar to the observation in Table 2, the calculated reliability index depends on the value of m used, and the reliability index obtained based on m = 1 is closest to the actual FORM reliability index. However, when the value of m is the same, the reliability index obtained based on the method suggested in this paper is closer to the actual FORM reliability index, probably because a response surface calibrated with non-uniformly distributed calibration points is less accurate for locating the design point. To study the effect of changing step size during iteration, the results of using m = 2 or m = 3 in the first iteration and m = 1 in the subsequent iterations are also shown in Table 3. Compared with the result obtained based on using a constant step size of m = 2 or m = 3, the reliability indexes obtained based on changing step size scheme are also closer to the actual FORM reliability index. It should be noted that although the classical RSM can be improved by imposing a lower bound value for xi when xci is less than mrxi, such a scheme cannot eliminate the possibility of generating negative values for positive random variables when determining xc with Eq. (4), which will also terminate the iteration. This will be illustrated in a numerical example studied later in this paper.

6. Automatic reliability analysis based on FLAC3D 6.1. General procedure Fig. 2 shows the flowchart of the procedure for automating reliability analysis with Matlab and FLAC3D. In this procedure, Matlab is used for reliability analysis and FLAC3D is used for assessing the response of the performance function, respectively. The FLAC3D model is called from Matlab using the routine !f3300.exe. The data exchange between FLAC3D and Matlab is realised through text

500

J. Zhang et al. / Computers and Geotechnics 69 (2015) 496–505

Fig. 4. Finite difference mesh of the slope example.

global xmean xsd xr xc yc yd Gm j=1; xmean = [9.8 10]; xsd = [3 2]; xr=[1,0;0,1]; yc= [0 0]; xc=getx(yc); dlmwrite('x.txt',xc','delimiter','\t','precision','%.2e','coffset',6); !f3300.exe FLACModel.dat Fsm=dlmread('FLACRes.txt') Gm=Fsm-1 delete('FLACRes.txt'); Fig. 2. Flow chart for implementing the suggested response surface method.

files. The procedure can be described as follows. First, the calibration points are determined in Matlab and saved into a text file. The data in this text file are entered into FLAC3D to obtain the responses of the performance function, which are also saved in a text file. The results from FLAC3D are then used as input for Matlab to build the response surface and calculate the reliability index. With the design point obtained in this iteration, a new set of calibration points can be obtained to update the response surface. If the discrepancy between the reliability indexes obtained in two successive iterations is not within a tolerable error, the response surface shall be further updated until the convergent reliability index is reached.

6.2. Illustrative example Herein, the aforementioned procedure is illustrated using a slope example adapted from Yamagami and Ueta [30], as shown in Fig. 3. A more complex geotechnical problem will be demonstrated in the subsequent portion in this paper to show the versatility and usefulness of this method. As Young’s modulus and Poisson’s ratio do not have an obvious effect on the factor of safety (FOS) of a slope [31], only the cohesion and the friction angle are modelled as random variables. Let c and u denote the cohesion

Young’s modulus: 0.6 MPa Poisson’s ratio: 0.25 Unit weight: 17.64 kN/m3 μc = 9.8 kPa, σc = 3.0 kPa, μϕ = 10o, σϕ = 2o

Fig. 3. Cross-section of the slope example analysed with FLAC3D.

m=2; Betaerr=10; while Betaerr>0.01 if j==1 [beta(j),yc,yd]=Getyd(yc,m); YD(j,:)=yd; YC(j,:)=yc; else [beta(j),yc,yd]=Getyd(yc,m); YD(j,:)=yd; YC(j,:)=yc; Betaerr=abs(beta(j)-beta(j-1)); end j=j+1; end Fig. 5. Matlab codes for the slope example.

and the friction angle of the soil, which are assumed to follow the lognormal distribution. Thus, x = {c, u}. The mean and standard deviation of c are 9.8 and 3.0 kPa, respectively. The mean and COV of u are 10° and 2°, respectively. In this study, the FOS of the slope is calculated using the shear strength reduction method implemented in FLAC3D [23], in which the FOS is defined as the number by which the shear strength parameters must be factored down to bring the geotechnical system to failure (e.g., Refs. [31,32]). The finite difference mesh generated is shown in Fig. 4. The FOS of the slope evaluated using the mean of the random variables is 1.37. For comparison, the FOS of the slope calculated using the shear strength reduction method implemented in program Slide is 1.33 [33]. The Matlab codes of this example are shown in Fig. 5. The Matlab function dlmwrite.m is used to write data into the file x.txt such that FLAC can read. Another Matlab function dlmread.m is used to read data from the file FLACRes.txt, which stores the FOS calculated by FLAC. The function !f3300.exe is used to call the FLAC codes stored in FLACModel.data as shown in Fig. 6 for computing the FOS of the slope based on soil parameters stored in x.txt. Other important functions in the Matlab codes are RSM.m, which calibrates the response surface given by Eq. (9)

J. Zhang et al. / Computers and Geotechnics 69 (2015) 496–505

501

Fig. 6. FISH codes for the slope example.

and also evaluates the reliability index based on FORM, and Getx.m, which converts variables in the reduced space into variables in the original space based on Eq. (6). The Matlab codes as shown in Fig. 4 are stored in the file main.m. The FLAC codes as in Fig. 6 consist of three parts, including reading data from Matlab through the text file x.txt, saving results from FLAC simulations to the text file FLACRes.txt, and evaluating the FOS of the foundation via the shear strength reduction method, respectively. The FLAC codes are compiled in the file FLACModel.data. The Matlab codes main.m and FLAC codes FLACModel.data are saved in the same folder. The main.m in Matlab can then be invoked to calculate the reliability index automatically. As an illustration, Table 4 shows the values of yc, yD, and the reliability index during the iteration process when m = 2 and eb = 0.01 are employed. It is observed that the reliability index converges to 1.480 in three iterations and the corresponding design point is {1.272, 0.756} in the reduced space or {6.404 kPa, 8.442°} in the original space.

7. Reliability of ground reinforced with DM columns 7.1. Problem description In the following study, the developed method is applied to the reliability analysis of a soft ground reinforced with DM columns. This example is adapted from Wu [34], and it is shown in Fig. 7. The foundation soils consist of 12.1 m of clay and 5.5 m of silty clay. The ground is subjected to a surface surcharge of 80 kPa with a width of 13.6 m. To reduce the settlement of the ground, DM columns with diameter of 0.9 m and length of 14 m are installed in an isolated pattern (in reference to Fig. 7). For simplicity, the circular isolated columns are modelled by square columns with an equivalent area in the numerical model as suggested by Huang and Han [35]. The equivalent width of the square columns is 0.8 m. The centre-to-centre distance between adjacent DM columns is 1.6 m. The allowable settlement of the foundation is 10 cm in this study [36]. As the settlement of the soil is most sensitive to the properties of the clay and the DM columns, the properties of the clay

502

J. Zhang et al. / Computers and Geotechnics 69 (2015) 496–505

Table 4 Intermediate results from the suggested RSM (slope example, m = 2, eb = 0.01). Iter. no.

yc

yD

c

u

c

u

1 2 3

1.265 1.293 1.288

0.814 0.763 0.766

1.189 1.277 1.272

0.765 0.753 0.756

b

1.414 1.483 1.480

(a)

(b) Fig. 8. FLAC3D mesh in the DM ground example: (a) elevation view and (b) plan view.

consisting of 6808 elements and 9000 nodes. Due to the symmetry around the centreline and the repetition of the span along the longitudinal direction, only a three-dimensional slice is modelled herein. The DM columns, the clay and the silty clay are all modelled as linearly elastic–perfectly plastic materials with the Mohr– Coulomb failure criterion. The boundary conditions adopted are explicitly shown in Fig. 7. The sensitivity analysis shows the current model size and the mesh produces the converged solutions. 7.3. Limit state function

Fig. 7. Layout of the DM columns: (a) cross-section and (b) plan view.

and the DM columns are modelled as random variables. Let Es, cs, and us denote Young’s modulus, the cohesion, and the friction angle of the clay, respectively. Let Ep, cp, and up denote Young’s modulus, the cohesion, and the friction angle of the DM columns, respectively. In this example, x = {Es, cs, us, Ep, cp, up}. All these random variables are assumed to follow the log-normal distribution. The statistics of random variables for the clay and DM columns are summarised in Table 5. The properties of the silty clay are treated as deterministic values as shown in Table 5. The problem of interest is to evaluate the serviceability reliability of this reinforced foundation. 7.2. FLAC3D model The mesh of the finite difference model is shown in Fig. 8. The model is 35 m in width, 17.6 m in depth, and 1.6 m in thickness,

After the model reaches the equilibrium with gravitational loading, there are two approaches to evaluate whether the settlement of the ground is less than the allowable value. In the first approach, one can directly apply the load of 80 kPa to the foundation and then evaluate the resulting displacement. If the calculated displacement is larger than 10 cm, it indicates serviceability failure. However, when the bearing capacity of the ground is less than or close to 80 kPa, the FLAC3D model may not converge or it may take a long time to converge. Alternatively, one may also apply equal vertical velocities to the nodes at the top of the model that are subjected to surface loading to obtain the load–settlement curve of the footing. Then the allowable load corresponding to the settlement of 10 cm based on the load–settlement curve can be assessed, as shown in Fig. 9. If the allowable load is <80 kPa, it indicates serviceability failure. In this velocity control method, the vertical settlement at the ground surface is uniform. Such a method is most appropriate when the overlaying footing is rigid and when the differential settlement is small. The velocity control method is employed in this study. Let q10(x) denote the allowable load corresponding to a settlement of 10 cm calculated by FLAC3D with input parameters x. Let qa denote the load applied to the footing (qa = 80 kPa in this example). The limit state function can then be written as follows:

gðxÞ ¼ q10 ðxÞ  qa

ð16Þ

Table 5 Material properties in the deep mixed (DM) ground example. Young’s modulus (MPa)

Cohesion (kPa)

Friction angle (°)

Poisson’s ratio

Unit weight (kN/m3)

Clay

Mean COV

4.68 0.3

16 0.15

7.6 0.25

0.33 –

18 –

DM columns

Mean COV

120 0.4

200 0.3

35 0.2

0.20 –

20 –

12.8

29

12.5

0.25

20

Silty clay

503

J. Zhang et al. / Computers and Geotechnics 69 (2015) 496–505

Fig. 9. Determination of the allowable bearing capacity based on the load– displacement curve.

instance, if m = 2 is used in the first iteration and m = 1 is used in the subsequent iterations, the value of xci  mrxi of Ep becomes negative in the second iteration. It seems very challenging to solve the current problem with the classical RSM. As mentioned previously, when xci is less than mrxi, one can impose a lower bound value for xi to avoid the occurrence of negative values. Table 7 shows the intermediate results during iteration when a lower bound value of 1 MPa is imposed for Ep. In the case of m = 1, a negative value is generated for Ep when determining xc with Eq. (4). A Similar phenomenon is also observed in the cases of m = 2 and m = 3. Thus, while the above scheme of imposing a lower bound value for xi can avoid generating negative values caused by xci  mrxi, it cannot eliminate the chance of generating negative values when determining xc with Eq. (4), which is a key step to bring the calibration points closer to the limit state function. In such a case, the iteration will also be terminated before convergence is achieved.

7.4. Reliability analysis with the classical RSM

7.5. Reliability analysis with the suggested method

The reliability index of the reinforced ground is first analysed with the classical RSM (in reference to Eq. (3)). Table 6 shows the value of xc and xD obtained during the iteration process when m = 1, 2, and 3 are adopted, respectively. In the cases of m = 1 and m = 2, negative values are generated for Ep when determining the calibration points with xci  mrxi during the second iteration, and FLAC3D cannot proceed with such inputs. Thus, no convergent reliability index can be reached when m is set to be 1 or 2. In the case of m = 3, a negative value is generated for Ep with xci  mrxi in the first iteration. Thus, m = 3 is not suitable for reliability analysis either. One may alter the step size m during the iteration process to avoid the occurrence of physically impermissible values. Although many step size adjustment schemes were attempted in this study, no convergent reliability indexes were reached. For

Using the RSM suggested in this paper, the reliability index of the reinforced ground is first solved with m = 1, and the results are shown in Table 8. The suggested method converges in three iterations, yielding a reliability index of 3.445. To study the effect of step size, the reliability of the reinforced ground is also solved with m = 2 and m = 3, respectively. Convergent reliability indexes can also be found for these two m values, as shown in Table 8. The reliability index calculated based on m = 2 is very close to that calculated based on m = 1. The reliability index calculated based on m = 3 is slightly larger than those obtained based on m = 1 and 2. As noticed previously, adjusting the step size during iteration helps improve the accuracy of the suggested RSM. To study the effect of changing step size during iteration on the suggested RSM in this example, the reliability index of the DM ground is also calculated

Table 6 Intermediate results from the classical RSM (DM ground example, eb = 0.01). Iter. no.

xc

xD

b

Es (MPa)

cs (kPa)

us (°)

Ep (MPa)

cp (kPa)

up (°)

Es (MPa)

cs (kPa)

us (°)

Ep (MPa)

cp (kPa)

up (°)

m = 1a

1 2

4.680 3.060

16.000 15.462

7.600 7.060

120.000 36.508

200.000 190.091

35.000 34.397

2.854

15.394

6.991

25.935

188.836

34.321

4.094

m = 2a

1 2

4.680 3.346

16.000 15.579

7.600 7.124

120.000 28.543

200.000 140.929

35.000 34.181

3.457

15.614

7.164

36.136

145.834

34.249

3.195

m = 3a

1

4.680

16.000

7.600

120.000

200.000

35.000

a

In the cases of m = 1, m = 2, and m = 3, negative values are generated for Ep when determining calibration points with xci  mrxi in the second iteration, the second iteration, and the first iteration, respectively.

Table 7 Intermediate results from the classical RSM improved by imposing lower bound values for non-negative random variables (DM ground example, eb = 0.01). Iter. no.

xc

xD

Es (MPa)

cs (kPa)

us (°)

b

Ep (MPa)

cp (kPa)

up (°)

Es (MPa)

cs (kPa)

us (°)

Ep (MPa)

cp (kPa)

up (°)

m = 1a

1 2 3 4

4.680 3.060 3.316 5.570

16.000 15.462 15.637 16.592

7.600 7.060 7.208 8.338

120.000 36.508 31.279 77.109

200.000 190.091 191.529 225.980

35.000 34.397 34.318 37.094

2.854 3.320 4.391

15.394 15.638 15.808

6.991 7.210 7.360

25.935 31.587 183.993

188.836 191.558 191.565

34.321 34.320 34.320

4.094 3.430 1.304

m = 2a

1 2 3

4.680 3.346 5.662

16.000 15.579 16.642

7.600 7.124 8.407

120.000 28.543 75.793

200.000 140.929 228.709

35.000 34.181 37.276

3.457 4.387

15.614 15.808

7.164 7.359

36.136 178.471

145.834 191.426

34.249 34.320

3.195 1.225

m = 3a

1 2 3

4.680 3.700 8.598

16.000 15.589 18.945

7.600 7.088 11.348

120.000 19.373 68.456

200.000 105.172 327.747

35.000 33.550 45.754

3.976 4.432

15.705 15.814

7.233 7.363

47.803 131.910

131.964 191.927

33.960 34.320

2.572 0.440

a In the cases of m = 1, m = 2, and m = 3, negative values are generated for Ep when determining calibration points with Eq. (4) in the fourth iteration, the third iteration, and the third iteration, respectively.

504

J. Zhang et al. / Computers and Geotechnics 69 (2015) 496–505

Table 8 Intermediate results from the suggested RSM (DM ground example, eb = 0.01). Iter. no.

yc

m=1

1 2 3

0 0 0 0 0 0.704 0.068 0.081 3.339 0.008 1.054 0.115 0.124 3.277 1.23E04

m=2

1 2 3

0 0 0 0 0 0 0.655 0.064 0.092 3.164 0.140 7.46E05 3.236 0.691 0.067 0.097 3.337 0.147 7.87E05 1.049 0.113 0.152 3.274 1.55E04 4.11E06 3.443 1.051 0.113 0.152 3.278 1.55E04 4.11E06 1.040 0.096 0.132 3.280 4.25E06 4.34E06 3.445

m=3

1 2 3

0 0 0 0 0 0.692 0.073 0.160 3.366 0.445 1.037 0.114 0.215 3.283 0.001

0 0.016 5.30E05

0 0 0 0 0 0.692 0.073 0.160 3.366 0.445 1.066 0.113 0.122 3.277 1.50E02

0 0.676 0.071 0.156 3.285 0.435 1.60E02 3.386 1.60E02 1.065 0.113 0.122 3.274 1.50E02 5.32E04 3.447 5.33E04 1.049 0.099 0.111 3.280 1.97E04 1.40E05 3.447

m = 3 in the first iteration, 1 2 and m = 1 in 3 subsequent iterations

Es

yD cs

us

Ep

cp

using m = 3 in the first iteration and m = 1 in the subsequent iterations, and the results are shown in Table 8. The obtained reliability index in such a case is 3.447, which is closer to the reliability indexes calculated based on constant step sizes of m = 1 and m = 2. Overall, compared with the shallow foundation example, the reliability analysis results in this example is less sensitive to the value of m. 7.6. Implication on design of DM ground The values of the design points in the reduced space can be used to quantify the contribution of each random variable to the overall uncertainty in the performance prediction [24]. Let y1, y2, y3, y4, y5, and y6 denote the random variables of Es, cs, us, Ep, cp, and up in the reduced space, respectively. The more the value of yi deviates from zero at the design point, the more the corresponding random variable affects the reliability of the problem [24]. In the case of m = 1, the values of y1, y2, y3, y4, y5, and y6 at the design point are 1.045, 0.099, 0.111, 3.279, 0.000, and 0.000, respectively. Thus, for the current reliability problem, Ec and Ep have the most significant contribution to the serviceability of the foundation, whereas the uncertainties in cs and cu are less important. The uncertainties in the shear strength parameters of the DM columns do not affect the reliability index of the reinforced ground, mainly because within the possible range of the random variables the DM columns are still in the elastic stage. In this regard, to reduce the uncertainties in the serviceability of this foundation, more efforts should be exerted on reducing the uncertainties in Ec and Ep. A practical engineering problem is to evaluate the allowable load of the DM reinforced ground. The reliability analysis

Fig. 10. Relationship between reliability index and the applied load.

up

Es

b cs

us

Ep

cp

up

0 0.669 0.065 0.077 3.170 0.007 1.45E05 3.241 1.53E05 1.054 0.115 0.124 3.278 1.23E04 8.12E06 3.447 8.12E06 1.045 0.099 0.111 3.279 5.41E06 6.93E06 3.445

0.676 0.071 0.156 3.285 0.435 0.016 3.386 1.036 0.114 0.215 3.282 0.001 5.29E05 3.450 1.026 0.096 0.190 3.290 1.23E05 6.37E06 3.453

described previously can be used to determine the allowable load applied to the foundation given a target reliability level. For the purpose of demonstration, the reliability index of the foundation is recalculated assuming q = 60 kPa and q = 100 kPa, respectively. The resulting relationship between the applied load and the reliability index of the foundation serviceability is shown in Fig. 10. As expected, the reliability index decreases as the applied load increases. If a target reliability index of 3.0 is chosen, the allowable load is about 90 kPa for this case study.

8. Summary and conclusions Current geotechnical reliability analyses mainly focus on problems with relatively simple limit state functions, and most numerical programs do not have the reliability analysis capability. In this paper, we suggest a modified RSM such that geotechnical reliability analysis can be performed conveniently using stand-alone deterministic numerical programs. The research work and findings presented in this paper can be summarised as follows: (1) A modified RSM constructed in the reduced space of random variables is suggested for geotechnical reliability analysis based on deterministic numerical programs. The suggested method can effectively avoid the negative values for positive random variables during the iteration procedure, and thus it solves the major limitation of the classical RSM in geotechnical applications. (2) A procedure is developed for automating geotechnical reliability analysis using the modified response surface method and FLAC3D. The detailed procedure is illustrated with a simple slope stability example. (3) The versatility of the procedure is illustrated through the reliability analysis of a foundation on soft ground improved by DM columns. It is shown that the classical RSM cannot yield a convergent reliability index. Instead, the modified RSM can calculate the reliability index of this problem efficiently. (4) For the serviceability reliability analysis of the DM ground, the convergence of the modified RSM is insensitive to the step size and the updating scheme. In addition, the reliability index is mainly affected by the uncertainties in Young’s modulus of the soil and Young’s modulus of the DM columns. The reliability-based design of DM foundation can be realised through determining the maximum allowable load given a target reliability level.

J. Zhang et al. / Computers and Geotechnics 69 (2015) 496–505

Acknowledgments This research was substantially supported by the Natural Science Foundation of China (41372275), the National 973 Basic Research Program of China (2014CB049100), the Shanghai Rising-star Program (15QA1403800), and the Fundamental Research Funds for the Central Universities. Dr. Abdul-Hamid Soubra from the University of Nantes kindly provided the authors the files for implementing the classical response surface method with FLAC3D, which is very helpful for developing the procedure for reliability analysis using the modified response surface method in this study. The help from Dr. Soubra is greatly appreciated. References [1] Christian JT. Geotechnical engineering reliability: how well do we know what we are doing? J Geotech Geoenviron Eng 2004;130(10):985–1003. [2] Phoon KK, Ching J. Risk and reliability in geotechnical engineering. New York: CRC Press; 2015. [3] Honjo Y, Suzuki M, Shirato M, Fukui J. Determination of partial factors for a vertically loaded pile based on reliability analysis. Soils Found 2002;42(5):91–109. [4] Paikowsky SG, Birgisson B, McVay M, Nguyen T, Kuo C, Baecher G, et al. Load and resistance factor design (LRFD) for deep foundations. NCHRP Rep. No. 507. Washington, DC: Transportation Research Board, National Research Council; 2001. [5] American Association of State Highway and Transportation Officials (AASHTO). LRFD bridge design specifications. 3rd ed. Washington, DC: AASHTO; 2004. [6] Griffiths DV, Huang J, Fenton GA. Influence of spatial variability on slope reliability using 2-D random fields. J Geotech Geoenviron Eng 2009;135(10):1367–78. [7] Huang J, Griffiths DV, Fenton GA. System reliability of slopes by RFEM. Soils Found 2010;50(3):343–53. [8] Wong FS. Slope reliability and response surface method. J Geotech Eng 1985;111(1):32–53. [9] Bauer J, Puła W. Reliability with respect to settlement limit-states of shallow foundations on linearly-deformable subsoil. Comput Geotech 2000;26(3– 4):281–308. [10] Goh ATC, Kulhawy FH. Neural network approach to model the limit state surface for reliability analysis. Can Geotech J 2003;40(6):1235–44. [11] Zhang LL, Zhang J, Zhang LM, Tang WH. Back analysis of slope failure with Markov chain Monte Carlo simulation. Comput Geotech 2010;37(7– 8):905–12. [12] Zhang J, Zhang L, Tang W. Slope reliability analysis considering site-specific performance information. J Geotech Geoenviron Eng 2011;137(3):227–38. [13] Zhang J, Zhang LM, Tang WH. Kriging numerical models for geotechnical reliability analysis. Soils Found 2011;51(6):1169–77. [14] Cho SE. Probabilistic stability analyses of slopes using the ANN-based response surface. Comput Geotech 2009;36(5):787–97.

505

[15] Li DQ, Jiang SH, Cao ZJ, Zhou W, Zhou CB, Zhang LM. A multiple responsesurface method for slope reliability analysis considering spatial variability of soil properties. Eng Geol 2015;187:60–72. [16] Bucher CG, Bourgund U. A fast and efficient response surface approach for structural reliability problems. Struct Saf 1990;7:57–66. [17] Rajashekhar M, Ellingwood B. Reliability of reinforced-concrete cylindrical shells. J Struct Eng 1995;121(2):336–47. [18] Guan XL, Melchers RE. Effect of response surface parameter variation on structural reliability estimates. Struct Saf 2001;23:429–44. [19] Chryssanthopoulos MK. Probabilistic buckling analysis of plates and shells. Thin-Wall Struct 1998;30(1–4):135–57. [20] Youssef Abdel Massih D, Soubra A. Reliability-based analysis of strip footings using response surface methodology. Int J Geomech 2008;8(2):134–43. [21] Mollon G, Dias D, Soubra A. Probabilistic analysis of circular tunnels in homogeneous soil using response surface methodology. J Geotech Geoenviron Eng 2009;135(9):1314–25. [22] Ji J, Low B. Stratified response surfaces for system probabilistic evaluation of slopes. J Geotech Geoenviron Eng 2012;138(11):1398–406. [23] FLAC3D. Fast Lagrangian analysis of continua in 3 dimensions. Minneapolis: Itasca Consulting Group Inc; 2005. [24] Low BK, Tang WH. Efficient spreadsheet algorithm for first-order reliability method. J Eng Mech 2007;133(12):1378–87. [25] EI-Ramly H, Morgenstern NR, Cruden DM. Lodalen slide: a probabilistic assessment. Can Geotech J 2006;43(9):956–68. [26] Rackwitz R. Reliability analysis-a review and some perspectives. Struct Saf 2001;23:365–95. [27] Phoon KK, Kulhwawy FH. Characterization of geotechnical variability. Can Geotech J 1999;36(4):612–24. [28] Sieffert JG, Bay-Gress Ch. Comparison of European bearing capacity calculation methods for shallow foundations. P I Civil Eng-Geotech 2000;143(2):65–74. [29] Cherubini C. Reliability evaluation of shallow foundation bearing capacity on c0 , u0 soils. Can Geotech J 2000;37:264–9. [30] Yamagami T, Ueta Y. Search noncircular slip surfaces by the Morgenstern-Price method. In: Proceedings of the 6th international conference on numerical and analytical methods in geomechanics. Rotterdam: Balkema Publications; 1988. p. 1335–40. [31] Griffiths DV, Lane PA. Slope stability analysis by finite elements. Geotechnique 1999;49(3):387–403. [32] Dawson E, Motamed F, Nesarajah S, Roth W. Geotechnical stability analysis by shear strength reduction. In: Griffiths DV, Fenton GA, Martin TA, editors. Slope stability 2000, proceedings of sessions of geo-Denver 2000. Reston, VA: ASCE; 2000. p. 99–113. [33] Rocscience. Slope stability verification manual Part I. Toronto: Rocscience Inc; 2006. [34] Wu L. Study on the mechanics properties of the strengthen system and the design calculation of high-speed railway pile-raft composite foundation. Master Thesis. Changsha, China: Central South University; 2012 [in Chinese]. [35] Huang J, Han J. 3D coupled mechanical and hydraulic modeling of a geosynthetic-reinforced deep mixed column-supported embankment. Geotext Geomembr 2009;27(4):272–80. [36] Ministry of Transportation (MOT). Technical specifications for design and construction of highway embankment on soft ground (JTJ 017-96). Beijing: China Communications Press; 1996 [in Chinese].