Settlement analysis of piled-raft foundations by means of a multiphase model accounting for soil-pile interactions

Computers and Geotechnics 46 (2012) 26–38

Contents lists available at SciVerse ScienceDirect

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

Settlement analysis of piled-raft foundations by means of a multiphase model accounting for soil-pile interactions E. Bourgeois a, P. de Buhan b,⇑, G. Hassen b a b

Université Paris-Est, IFSTTAR-MACS, 58 bd Lefebvre, 75732 Paris Cedex 15, France Université Paris-Est, Laboratoire Navier, Ecole des Ponts ParisTech/IFSTTAR/CNRS, 6-8 Avenue Blaise Pascal, 77455 Marne-la-Vallée Cedex 2, France

a r t i c l e

i n f o

Article history: Received 17 February 2012 Received in revised form 18 May 2012 Accepted 28 May 2012 Available online 29 June 2012 Keywords: Piled raft foundations Settlement analysis Load bearing capacity Multiphase model Soil pile interactions Finite element Elastoplasticity

a b s t r a c t The behavior of vertically loaded raft foundations, strengthened by a large group of floating piles, is investigated by means of a multiphase approach, where special emphasis is put on the ground-pile interactions, which play an important role on the foundation performance. A new specific interaction law relating to the pile tip resistance is incorporated into the multiphase model, in addition to that already developed for pile shaft friction in a previous version. A numerical identification procedure is presented, which allows evaluation of the stiffness and yield strength parameters governing such interaction laws, as functions of the different geometric characteristics defining the pile layout (pile diameter and spacing), along with the material constitutive properties of the soil, modeled as a purely cohesive clay. The improved multiphase model, implemented in a finite element code, is used to analyze the response of a square raft lying upon a group of piles, varying the pile number and pile length. A relatively good agreement is obtained between the results of this approach, expressed in terms of load-settlement behavior as well as pile force distributions, and those derived from direct computationally intensive finite element calculations, where piles are regarded as individual elements embedded in the soil. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction The rational design of foundations reinforced by large groups of piles is an important issue in geotechnical engineering, both in terms of settlement reduction and increase of the load bearing capacity to be expected from such a reinforcement scheme. Various numerical methods have been proposed for addressing this kind of problem, ranging from the ‘‘hybrid model’’ [8,12,10], to variational approach-based methods [9,18,20] or to conventional, but more comprehensive and potentially accurate, finite element analyses [2,16,17,19,21,22,24]. Unfortunately, as a result of the very strong heterogeneity of the pile reinforced zone, the above mentioned finite element simulations are very difficult to implement, first as concerns data preparation and mesh generation, then in terms of oversized numerical problem and correlative computational costs, notably when three dimensional configurations as well as elastoplastic behavior have to be considered. In order to overcome such difficulties, a so called ‘‘multiphase model’’ of the pile reinforced ground has been proposed, aimed at predicting the global response of piled raft foundations subject to purely vertical [23] or combined loadings [13]. In an early version of this model, the composite reinforced ground was implicitly regarded as an anisotropic homogenized medium, ⇑ Corresponding author. E-mail address: [email protected] (P. de Buhan). 0266-352X/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compgeo.2012.05.015

with no specific attention paid to the interactions prevailing at the soil-pile interfaces allowing for possible failure and slippage. A first extended version of the model has been therefore developed, where the previous ‘‘perfect bonding’’ assumption should be abandoned and soil-inclusion interaction accounted for Ref. [3]. According to this generalized model, the reinforced ground is no more regarded as a single homogenized medium, but as the superposition of two mutually interacting continuous media, called ‘‘phases’’ (hence the name of ‘‘multiphase’’ model). A specific law accounting for soil pile interaction along the pile shaft has been formulated in the framework of the model, then implemented in the finite element package software CESAR-LCPC [15] and applied to the numerical simulation of the convergence of bolt-supported tunnels [5]. It has also been applied in a quite recent paper [4] to the simulation of piled raft foundations, providing more realistic predictions for the settlement behavior and pile force distributions than those obtained with the model based on the ‘‘perfect bonding’’ hypothesis. More specifically, comparisons with ‘‘direct’’ finite element calculations (where the piles and surrounding soil are discretized as separate volume elements) clearly showed that the perfect bonding model underestimated the foundation settlement under a given load level, while at the same time quite significantly overestimated the pile force distributions. The comparisons showed that the extended model significantly improved the results, but still remained insufficiently accurate, because it led to axial forces

27

E. Bourgeois et al. / Computers and Geotechnics 46 (2012) 26–38

vanishing at the pile tips, contrary to what is observed in ‘‘direct’’ simulations, where the axial force at the pile tips represents a significant proportion of the maximum force at the pile head. In order to overcome such a drawback, the present paper proposes a significant additional enhancement of the multiphase model, where a specific soil-pile interaction law, associated with the ‘‘pile tip resistance’’, is incorporated into the multiphase model, in addition to the already introduced ‘‘pile shaft friction’’ interaction law. This second kind of soil pile interaction behavior is numerically implemented in a finite element formulation of the multiphase model and then applied to the settlement analysis of a piled raft foundation, where the influence of the soil-pile tip interaction on the foundation response, as well as the pile force distributions, is clearly highlighted. Furthermore, a procedure for evaluating the soil-pile interaction parameters to be introduced in the multiphase calculations is outlined and discussed. 2. Multiphase model incorporating ground-pile interactions 2.1. Problem statement and general principle of the multiphase description The typical geotechnical problem to be dealt with is that of a vertical loading applied to a semi-infinite soil mass through a square surface footing of side B as sketched in Fig. 1a. In order to increase its bearing capacity, the soil (soft clay) has been previously reinforced by a group of vertical piles of length L placed just beneath the footing. Denoting by q the pile’s radius and by s the spacing between two adjacent piles assumed to be distributed in the groundmass following a regular pattern, a key parameter of such a reinforcement scheme is the reinforcement volume fraction defined as:

g ¼ pðq=sÞ2

ð1Þ

which generally proves to remain a small quantity. The multiphase approach applied to the design of such a piledraft foundation consists in replacing the composite pile-strengthened volume of soil of height L and side B, not by one single equivalent medium as in the classical homogenization method, but by two superposed interacting continua, called ‘‘phases’’. According to this model, a general presentation of which may be found for instance in [23], the soil is represented by the matrix phase, while the group of piles is represented by the reinforcement phase. More precisely, two coincident particles attached to each phase respectively, are located at any geometrical point of the reinforced zone; each particle is attributed its own kinematics, namely a displace-

B

2.2. Statics of the pile-reinforced zone as a two-phase system The system is referred to an orthonormal frame Oxyz, with the origin O placed at the center of the contact area between the raft and the soil. According to the two-phase description of the piledraft foundation, the matrix phase is modelled as a classical continuum at any point of which the stress is described by a second order tensor rm (Fig. 2a), while the stress at each point of the reinforcement phase, which occupies the reinforced domain V, is defined by a uniaxial tensor nrex  ex, where ex is the downwards oriented unit vector parallel to the pile orientation (Fig. 2b), and the scalar nr represents the density of axial force in the piles per unit cross sectional area of reinforced ground, which may therefore be calculated as:

nr ¼ Nr =s2

ð2Þ

where Nr is the axial force in one individual pile. The equilibrium equations can thus be expressed for each phase separately as:

divrm þ cm ex þ I ¼ 0

ð3Þ

for the matrix phase and:

divðnr ex  ex Þ þ cr ex  I ¼ 0

ð4Þ

for the reinforcement phase. In these equations cm (respectively cr) denotes the specific weight of the matrix (resp. reinforcement) phase, and I (resp. –I) is the body volume force density modelling at any point the action of the reinforcement (resp. matrix) phase onto the matrix (resp. reinforcement) phase. Owing to the fact that the reinforcement volume fraction g is small (typically less than 5%), the reinforcement phase specific weight defined as:

B

y

L

y

L

ρ

x soft clay

ment vector nm for the matrix, and nr for the reinforcement phase (Fig. 1b). It should be noted that, in order to capture the shear and flexural behavior of the reinforcing piles, a more general multiphase model could be adopted, in which the kinematics of the reinforcement phase is characterized by a rotation xr in addition to the sole translation nr introduced above. Nevertheless, it has been shown in [13], that the shear and flexural contributions of the reinforcing piles play a negligible role in the case of vertically loaded foundations, as that considered in this paper. That is why the simplified multiphase model, in which such shear and flexural effects are disregarded, will be adopted from now on.

pile s

(a)

x

matrix

ξm

s reinforcem ent

(b)

Fig. 1. Piled raft under vertical loading: (a) initial problem and (b) multiphase description.

ξr

28

E. Bourgeois et al. / Computers and Geotechnics 46 (2012) 26–38

y

σm

γ m ex

L

− I ex

I ex

nr e x ⊗ e x

interaction forces

Σ

pe x

γ r ex ≅ 0 Σ

− pe x

x

x

(a)

(b)

Fig. 2. Statics of the piled-raft foundation modelled as a two-phase system (although geometrically superposed, both phases have been shifted apart for the sake of clarity).

cr ¼ gcp

ð5Þ

where cp is the specific weight of the pile constituent material (such as concrete), which is of the same order of magnitude as that of the soil (cs), can be neglected with respect to the matrix phase specific mass. Indeed:

g  1 ) cr ¼ gcp  cm ¼ ð1  gÞcs ffi cs

ð6Þ

Under such circumstances, the reinforcement phase equilibrium Eq. (4) simplifies to:

divðnr ex  ex Þ þ cr ex  I ¼ ð@nr =@xÞex  I ¼ 0 |{z}

ð7Þ

ffi0

which implies that the interaction body force density is reduced to one single component along the pile direction:

I ¼ Iex

with I ¼ @nr =@x

ð8Þ

2.3. Two kinds of soil-pile interactions Starting from the equilibrium Eq. (8), combined with definition (2), one obtains:

@nr 1 @Nr I¼ ¼ 2 s @x @x

@Nr ¼ 2pqs @x

where s–s is the shear stress exerted by the soil (resp. pile) onto the pile shaft (resp. soil), which yields:

I¼

2pqs s2

ð11Þ

The latter relationship allows interpretation of the body force volume density I, at the macroscopic scale of the multiphase model, as the equivalent shear stress distribution applied by the pile onto the surrounding soil, due to the side or shaft friction developed along the soil-pile interface. A second kind of soil-pile interaction, which will play an important role in the analysis, is that developed at the tips, and more specifically in the present case, at the lower tips of the piles, located at depth x = L, as shown in Fig. 3b. This interaction is modelled at the macroscopic scale of the multiphase model, by means of a surface density of forces p applied to the reinforcement phase along its bottom surface R (Fig. 2b). This implies the following boundary condition for the reinforcement phase:

p ¼ nr ðx ¼ LÞ ¼

ð9Þ

Referring to the equilibrium of an infinitesimal element of pile, as sketched in Fig. 3a, we may write:

ð10Þ

Nr ðx ¼ LÞ s2

ð12Þ

while generating at the same time in the matrix phase, as a direct consequence of the principle of mutual interactions, a discontinuity of the vertical stress component across R (Fig. 2a):

Nr s

ρ

s

s

dx

τ s

N r + dN r (a)

N r ( x = L) (b)

Fig. 3. Soil-pile interactions: (a) along the pile length; (b) at the bottom pile tips.

29

E. Bourgeois et al. / Computers and Geotechnics 46 (2012) 26–38

rmxx ðx ¼ Lþ Þ  rmxx ðx ¼ L Þ ¼ p

ð13Þ

2.4. Constitutive equations of the multiphase model The procedure for deriving the constitutive behavior of each individual phase, already explained for instance in [23], is quite straightforward. Assuming that the piles remain elastic, the constitutive equation of the reinforcement phase may simply be written:

nr ¼ aer

ð14Þ

where er ¼ @nrx =@x is the axial strain of the reinforcement phase, while a represents the axial stiffness of the piles per unit transverse area, which can be evaluated as the product of the reinforcement volume fraction by the Young’s modulus of the pile constituent material:

a¼

nr

e

r

¼

Nr Ep ðpq2 Þer ¼ ¼ gEp r 2 s e s2 er

ð15Þ

On the other hand, the matrix phase constitutive relations are simply identified with those of the soil, namely a soft clay, which will be modeled in the sequel as an elastic perfectly plastic purely cohesive material (Tresca’s yield condition). Furthermore, the two types of soil-pile interactions, previously introduced, are governed by specific constitutive laws which can be expressed as follows: The constitutive behavior of the first kind of interaction relating to the soil-pile ‘‘shaft’’ or ‘‘side’’ friction is formulated by means of a relationship linking the interaction volume density I to the relative axial displacement between the reinforcement and the matrix phases, defined as:

D¼

nrx



nm x

ð16Þ

Similarly, the second type of interaction, associated with the ‘‘pile tip resistance’’, will be expressed by a relation between the interaction force surface density p and the axial reinforcement/matrix relative displacement D(L) on R (x = L). Again, adopting an elastic perfectly plastic constitutive model, we thus obtain:

8 0 > < 0 if jpj 6 p p p p _ p ¼ c ðDðLÞ  D ðLÞÞ with D ¼ P 0 if p ¼ þp0 ; p_ ¼ 0 > : 6 0 if p ¼ p0 ; p_ ¼ 0

ð18Þ

Both ‘‘shaft friction’’ and ‘‘tip resistance’’ interaction constitutive laws are represented in Fig. 4 below, in the form of classical bilinear stress-stress diagrams. Evidently, these interaction constitutive laws are strongly similar to load-transfer curves (‘‘t–z’’ curves) classically introduced in the literature on pile foundations for designing the load bearing capacity of individual piles driven in a soil (see for instance among the most recent references: [1]. It should be emphasized that, apart from the fact that such load transfer curves refer to individual piles and not to a ‘‘reinforcement phase’’ perceived as a homogenized 3D continuum, an important difference between the above interaction laws and the usual load-transfer curves, is that the relevant kinematic variable associated with the pile-ground interaction forces is the relative displacement (D) between the matrix and reinforcement phases and not the absolute settlement (‘‘z’’) of the pile. This basically comes from the fact that in most design methods based upon such load-transfer curves, the sole reaction of the soil mass on the system of piles and raft is considered through such curves, without any explicit consideration of the behaviour of the soil mass as a loaded deformable continuum. 3. A numerical procedure for evaluating the soil-pile interaction multiphase parameters

In the context of an elastic perfectly plastic behavior, such a constitutive law will take the following form:

3.1. Outline of the procedure

8 0 > < 0 if jIj 6 I p I p I ¼ c ðD  D Þ with D_ ¼ P 0 if I ¼ þI0 ; I_ ¼ 0 > : 6 0 if I ¼ I0 ; I_ ¼ 0

It has been shown that the constitutive elastic parameter of the reinforcement phase, describing the overall behavior of a group of piles can be readily determined from Eq. (15), while the elastoplastic parameters of the matrix phase are identical to those of the soil. On the other hand, the evaluation of the stiffness (cI and cp) as well as the yield strength (I0 and p0) parameters governing the two kinds of interaction behavior between phases, is a little bit more complex. A possible way could be for instance to perform a backanalysis of static loading tests carried out on single piles, fitting the different parameters from the comparison between the experimental results and those obtained from a numerical simulation of this kind of test [11].

ð17Þ

where Dp is the plastic component of the relative displacement, while cI and I0 are coefficients characterizing the stiffness of the interaction and the threshold value of the interaction force density for which an irreversible (plastic) relative displacement Dp between the matrix and reinforcement phases occurs. The latter parameter can for instance be related to the maximum skin friction between the piles and the ground [4].

p

I

p0

I0 cI

Δp -I0 (a)

Δ

cp

Δ (L )

Δp ( L)

- p0 (b)

Fig. 4. Stress–strain diagrams for (a) ‘‘shaft friction’’ and (b) ‘‘pile tip’’ interaction laws.

30

E. Bourgeois et al. / Computers and Geotechnics 46 (2012) 26–38

The procedure adopted in this paper for evaluating such parameters is of the numerical type, quite similar to that developed in the context a purely elastic behavior for the design of piled embankments [7,14], where the pile tip interaction takes place at the pile head, and not at the lower end as in our situation. This procedure can be generally described as follows: The procedure is based on considering the simplified boundary value problem (or auxiliary problem) shown in Fig. 5, which is closely related to the configuration of pile loading tests. A soil layer of thickness L + h and infinite horizontal extent is placed on a rigid substratum. It is reinforced by a uniform distribution of vertical piles driven in the soil from its upper surface down to a depth L, an increasing vertical load Q being applied on top of each pile. Analyzing the response of such a reinforced soil layer in terms of load-settlement curves can be achieved either through a numerical simulation of the auxiliary problem displayed in Fig. 5b (see later on), or by solving the problem where the reinforced soil layer is modeled as a two-phase system involving the different interaction parameters to be identified, as shown in Fig. 5a.

D00 ¼

D l

ð22Þ

2

where l ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi aEm oed

cI ðaþEm oed Þ

may be written:

D ¼ A sinhðx=lÞ þ B coshðx=lÞ

ð23Þ

Constant A is determined from the boundary conditions imposed on the upper surface x = 0 (Fig. 5a):

rmxx ð0Þ ¼ Emoed m0 ð0Þ ¼ 0

nr ð0Þ ¼ ar 0 ð0Þ ¼ q and

ð24Þ

which gives:

q ql D0 ð0Þ ¼ r0 ð0Þ  m0 ð0Þ ¼  ) A ¼ 

a

a

ð25Þ

The second constant B is calculated as follows from the different equations governing the pile tip interaction taking place at x = L. As a result of Eqs. (12) and (13), we have:

rmxx ðx ¼ Lþ Þ rmxx ðx ¼ L Þ ¼ p |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl}

p ¼ nr ðx ¼ LÞ;

3.2. Solving the elastic multiphase boundary value problem

is a characteristic length. The solution of (22)

ð26Þ

q

Referring first to a purely elastic analysis of the multiphase problem, where m(x) and r(x) denote the (vertical) settlements of the matrix and reinforcement phases at depth x, respectively, the settlement of the lower un-reinforced zone can be easily calculated as follows. Due to the global vertical equilibrium, the vertical stress at the top of this zone (or at the base of the reinforced zone) is:

rmxx ðx ¼ Lþ Þ ¼ q ¼ Q =s2

ð19Þ

so that, making use of the matrix and reinforcement phases elastic laws, along with the pile tip interaction elastic constitutive behavior (18), one finally obtains:

 m 0   E m ðL Þ  cp ½r  mðLÞ ¼ q x ¼ L :  oed ar0 ðLÞ þ cp ½r  mðLÞ ¼ 0 and consequently:

and the settlement profile across this zone is simply:

L6xþh:

mðxÞ ¼

q Em oed

ðx  ðL þ hÞÞ

ð20Þ

m m where ðkm ; Gm Þ are the Lamé’s constants and Em oed ¼ k þ 2G represents the constrained modulus of the matrix phase. In the upper reinforced zone, modeled as a two-phase system, the combination of the equilibrium and constitutive equations pertaining to each phase, leads to the following differential system (see [7] or [14], for more details):

 m 00  Eoed m þ cI ðr  mÞ ¼ 0 |fflffl{zfflffl}  0 6 x 6 L :  D  ar 00  cI ðr  mÞ ¼ 0

ð21Þ

DðLÞ q ¼ m l Eoed

ð28Þ

where j is a non-dimensional factor defined as:

j ¼ lcp



a þ Emoed aEmoed

 ð29Þ

Introducing (23) and (25) into Eq. (28) yields the value of constant B and then:   ql ða=Em oed Þ=coshðL=lÞ þ 1 þ j tanhðL=lÞ DðxÞ ¼ coshðx=lÞ  sinhðx=lÞ a j þ tanhðL=lÞ

As a final result, the stress distribution in the reinforcement phase can be calculated from integrating the second equation of (21):

nr ðxÞ ¼ ar 0 ðxÞ ¼ nr ð0Þ þcI |fflffl{zfflffl} q

q = Q / s2

Q

n ( x = 0) = −q r

L

σ xxm ( x = 0) = 0 x=L

Z

x

Dð1Þd1

ð31Þ

0

that is from (30):

a Em oed nr ðxÞ ¼ q m þq a þ Eoed a þ Emoed   ða=Em oed Þ=coshðL=lÞ þ 1 þ j tanhðL=lÞ  sinhðx=lÞ  coshðx=lÞ j þ tanhðL=lÞ ð32Þ

h x = L+ h

(a)

D0 ðLÞ þ j

ð30Þ

2

(with ðÞ00 ¼ d ðÞ=dx2 ). This system of equations leads to the following second order differential equation in D = r  m:

ξ xm ( x = L + h ) = 0

ð27Þ

x

s (b)

Fig. 5. Simplified boundary value problem for the evaluation of soil-pile interaction parameters: (a) multiphase description; (b) representative volume.

The corresponding stress profiles have been drawn in Fig. 6 in a non-dimensional form, for typical values of the elastic parameters (a ¼ 10Em oed and l ¼ 0:3 L), while varying the value of the nondimensional parameter j from 0 (no pile tip reaction) to 10. As it is clearly apparent from this figure, the pile tip reaction represents a growing proportion of the total load applied on the pile head as the parameter j is increased.

31

E. Bourgeois et al. / Computers and Geotechnics 46 (2012) 26–38

0 x L : I ¼ divrm ¼ divðnr ex  ex Þ ¼ I0 ex

nr / q 0

-0.2

-0.4

-0.6

-0.8

-1

while the pile tip interaction force surface density can be deduced from either (12), (13):

pile head ( x = 0)

0

ð35Þ

0.1

þ  m r 0 p ¼ rm xx ðL Þ þ rxx ðL Þ ¼ n ðLÞ ¼ p

0.2

thereby proving that the generalized stress field is satisfying the different yield strength conditions (including the matrix phase) while remaining in equilibrium with the following load value:

0.3 0.4

q ¼ nr ð0Þ ¼ p0 þ I0 L

x / L 0.5

ð36Þ

which therefore represents a lower bound estimate for the limit load:

0.6

qþ p0 þ I0 L

0.7

10

0.8

κ =0

0.9

0.2

0.5

1

It can be easily shown from using the upper bound kinematic approach of yield design, that this value is actually equal to the ultimate bearing capacity:

2 pile tip ( x = L)

1

ð37Þ

Fig. 6. Elastic stress distributions along the reinforcing piles evaluated from the multiphase model with pile shaft and tip interactions.

qþ ¼ p0 þ I0 L

ð38Þ

4. Identification of the interaction parameters

3.3. Ultimate bearing capacity

4.1. Elastic parameters (cI, cp)

Our attention is now focused on evaluating from the same multiphase model, the maximum load q+ which can be applied on top of the piles on account of the pile-soil interaction yield strength parameters (I0, p0). This problem can be formulated in the framework of the yield design (or limit analysis) approach applied to multiphase systems, as explained in de Buhan and Hassen [6], where the soil reinforcement tip interaction was not considered. Thus, referring to the lower bound static method of yield design, the following ‘‘generalized stress field’’ (rm, nr, I, p) is considered in the pile-strengthened soil layer modeled as a two-phase system (Fig. 7): Reinforcement phase (0 6 x 6 L):

Generally speaking, both elastic interaction parameters can be written as functions of the pile and soil elastic coefficients, as well as geometrical characteristics such as the pile radius and spacing:

nr ðxÞ ¼ p0 þ I0 ðx  LÞ

ð33Þ

Matrix phase (0 6 x 6 L + h):

 0  I x1

for 0 6 L

rm ðxÞ ¼ 

ð34Þ

 ðI0 L þ p0 Þ1 for Lþ 6 x 6 L þ h

As a result of both equilibrium equations (3) and (4), the pile shaft interaction force volume density is:

cI ¼ cI ðEs ; ms ; Ep ; mp ; q; sÞ;

cp ¼ cp ðEs ; ms ; Ep ; mp ; q; sÞ

ð39Þ

It follows from simple dimensional analysis considerations, that the above relationships can be rewritten in a non-dimensional form as:

c I s2 ¼ fI Es



ms ;

 Ep p q ; m ; ; s Es

cp s ¼ fp Es



ms ;

 Ep p q ; m ; s Es

ð40Þ

that is on account of (1):

cI ¼

Es I g s2



ms ;

 Ep p s ;m ;g ; E

cp ¼

Es p g s



ms ;

 Ep p s ;m ;g E

ð41Þ

This means that, all the other parameters being kept constant (notably the reinforcement volume fraction g), the pile shaft interaction coefficient cI (resp. pile tip interaction coefficient cp) varies in direct proportion to 1/s2 (resp. 1/s). These relationships are confirmed by finite element simulations performed on the auxiliary problem sketched in Fig. 5b with the following characteristics:

Es ¼ 45 MPa; ms ¼ 0:3; Ep ¼ 20 GPa;mp ¼ 0:2; L ¼ 20 m; h ¼ 4 m

σ xxm

− I 0L

− ( p 0 + I 0 L)

ð42Þ

nr

where three configurations corresponding to different values of the spacing s and pile radius q have been considered:

8 q ¼ 0:15 m > < s ¼ 2 m; s ¼ 2:5 m; q ¼ 0:187 m > : s ¼ 3:33 m; q ¼ 0:25 m

L −((pp00++II00LL))

so that the reinforcement volume fraction remains constant equal to:

− p0

g ¼ pðq=sÞ2 ffi 1:77% L+h

ð43Þ

x

x Fig. 7. Distributions of vertical stress at failure in the matrix and reinforcement phases.

ð44Þ

As it will be seen later on, these configurations correspond for instance to a square raft of side B = 10 m lying upon a soil reinforced by a group n = 9, 16 or 25 piles, as is shown for n = 9 in Fig. 8, where the finite element mesh of one eighth of the cross section is also displayed.

32

E. Bourgeois et al. / Computers and Geotechnics 46 (2012) 26–38

0

-0.0005

0 2 4

soil

-0.001

-0.0015

settlement (m)

-0.002

-0.0025

f.e.m (n=9) f.e.m. (n=16) f.e.m. (n=25) mult. (n=9) mult. (n=16)

6

ρ = 0.25m

mult. (n=25)

pile 8

s / 2 = 1.67m

10 depth (m)

Fig. 8. Finite element discretization of the auxiliary problem.

In each of the three configurations, a uniform pressure of r r– 1 MPa has been imposed on the top section of the pile, so that the value of the compressive stress in the reinforcement phase is equal to:

nr ðx ¼ 0Þ ¼ q ¼ 

rðpq Þ s2

¼ gr ¼ 0:0177 MPa

ð45Þ

  n ¼ 9 : cI ¼ 7:98 MN=m4 ; cp ¼ 4:66 MN=m3    n ¼ 16 : cI ¼ 14:03 MN=m4 ; cp ¼ 6:32 MN=m3   n ¼ 25 : cI ¼ 21:92 MN=m4 ; cp ¼ 7:95 MN=m3

ð46Þ

It can be verified that, with a good accuracy, cI is proportional to n (and then to 1/s2) while cp is proportional to n1/2 (and thus to 1/ s):

cp ffi 1:68n1=2

ð47Þ

The same trend is observed from the identification of the interaction parameters on the basis of the settlement profiles in the reinforcement phase (Fig. 10). These parameters are reasonably well approximated by almost the same formulas as (47):

cI ffi 0:83n;

4.2. Yield strength parameters (I0, p0)

2

The stress distribution in the reinforcement is then calculated at different depths from the results of the finite element simulations, through the same equation as (45) (points in Fig. 9). An optimization procedure is then carried out selecting the couple of interaction stiffness parameters (cIcp) for which the numerical series of points are best fitted by the analytical expression (32) derived from the multiphase analysis (solid curves of Fig. 9). The following couples of interaction stiffness parameters are obtained from the above procedure:

cI ffi 0:88n;

Fig. 10. Identification of interaction stiffness parameters from settlement curves in the reinforcement phase.

cp ffi 1:55n1=2

The evaluation of I0 can be straightforwardly derived from Eq. (11) involving the shear stress s exerted by the soil onto the pile shaft. Since the latter is limited by the cohesion C of the surrounding soil, assumed to be perfectly bonded to the pile, one obtains on account of (1):

I0 ¼

pffiffiffiffiffiffiffi C 2pqC ¼ 2 pg / n1=2 s2 s

ð49Þ

implying that for g being kept constant, the value of I0 is proportional to 1/s or n1/2. The determination of the other yield strength parameter p0 is then simply deduced from Eq. (38), where q+ denotes the limit load of the auxiliary problem, determined from a finite element based elasto-plastic simulation conducted with the same numerical model as for the elastic simulation. Thus:

p0 ¼ qþ  I0 L

ð50Þ

For a purely cohesive soil with a cohesion of C = 30 kPa, 9 numerical simulations, varying the pile density (n = 9–16–25) as well as the pile length (L = 8 m–12 m–16 m), have been performed. The results are as follows: It should be noted that, unlike what has been observed for the interaction stiffness cp, the pile tip interaction strength p0 is a (roughly linear) increasing function of the pile length.

ð48Þ 4.3. Comments on the ‘‘perfect bonding’’ multiphase model

0

-0.005

0

-0.01

-0.015

-0.02 r

n (MPa )

f.e.m. (n=9) f.e.m.(n=16)

2

f.e.m. (n=25) mult. (n=9) mult. (n=16)

4

The simplified multiphase model based upon the ‘‘perfect bonding’’ assumption is a particular case of the above described model, where both phases have the same kinematics, or equivalently when both stiffness and yield strength parameters governing the matrix-reinforcement interaction laws take infinite values. This perfect bonding assumption requires that the two following conditions are satisfied.

mult. (n=25)

6

8

10

depth (m)

Fig. 9. Fitting the interaction stiffness parameters from the finite element simulations.

(a) The individual piles are perfectly adherent to the surrounding soil, which means that no slippage or displacement discontinuity is allowed at the pile-soil interface. (b) As shown by the previous identification procedure of the interaction parameters, the later are increasing functions of the pile density (that is the pile number n: see Eqs. (48) and (49)). In other words, the multiphase model with perfect bonding assumption is obtained as a limit situation when the pile density tends to infinity.

33

E. Bourgeois et al. / Computers and Geotechnics 46 (2012) 26–38

5. Application to the settlement analysis of piled-raft foundations 5.1. Problem statement As an application of the above developments, we are now focused on predicting the settlement of a square raft of side B = 10 m lying upon a homogeneous soft clay, with the elastic properties previously defined in Eq. (42) and a cohesion C = 30 kPa. The raft being subjected to an increasing vertical uniform pressure, the analysis of this problem can be restricted to one eighth of the structure as shown in Fig. 11 which provides a perspective view along with a top view of the calculation model (finite element code CESAR-LCPC), where the piles have been discretized as 3D volume elements. As an indication of the numerical problem size, the total number of nodes varies from 20,221 for a group of 900208 m-long piles to 28,900 for 25 16 m-long piles. Fig. 12 displays the finite element mesh of the same structure where the pile-reinforced ground volume placed beneath the raft has been modelled as a two-phase system. Since there is no need to refine the mesh in the reinforced zone as in the previous case, the number of nodes can be reduced to 7000 with, as a direct consequence, a strongly reduced computational time. Another decisive advantage of the multiphase model with respect to a direct simulation where the piles are discretized as individual elements, lies in the fact that the same mesh can be used whatever the number on

piles or their geometrical characteristics (spacing, diameter). The multiphase calculations have been performed for the three reinforcement configurations corresponding to the different values of radius and spacing given in (43) and for the same mechanical properties (42) of piles, supposed to remain elastic. For these three configurations, the interaction stiffness and strength parameters have been already identified and are reported in Eq. (46) and Table 1. 5.2. Load-settlement curves Fig. 13 represents different load-settlement curves obtained from the f.e.m. multiphase simulation, where the length L and Table 1 identification of yield strength interaction parameters as functions of pile length and spacing. Pile number

I0 (kPa/m)

p0 (kPa)

n=9

4.24

n = 16

5.65

n = 25

7.06

L = 8 m: 3.53 L = 12 m: 4.24 L = 16 m: 5.48 L = 8 m: 5.12 L = 12 m: 6.54 L = 16 m: 7.77 L = 8 m: 6.54 L = 12 m: 8.30 L = 16 m: 9.89

10 m raft

L = 16 m

24 m piles

B/2 = 5 m Fig. 11. Finite element mesh of the piled raft foundation (n = 25 16 m-long piles).

two - phase zone

Fig. 12. Finite element mesh of the piled raft foundation modelled as a two-phase system.

34

E. Bourgeois et al. / Computers and Geotechnics 46 (2012) 26–38

0

10

20

30

40

total load (MN) -20 -50 -80

In order to be more specific, computed settlements associated with a load level equal to 20 MN have been reported in Fig. 14. While this load value is approximately equal to the ultimate bearing capacity of the un-reinforced foundation, and would thus induce a failure of the latter, the beneficial role of both pile length

perfect bonding

Settlement (mm) for 20 MN

180

L =8m

-110

160

unreinforced

-140

140 120

25

-170

L=8m

100

0

10

L=12m

80

-200 settlement (mm) n = 9 16

L=16m

60

20

30

40

total load (MN)

40 20 0

-20 -50

0

-140

25

-170

0

10

0

5

15

20

n=9

25

30

total load (MN) 8m/M

-100

16

12 m / M

20

30

40

total load (MN)

-150

16 m / M 8m/D

-200

n =16

12m / D

-250

16m / D

-50

-300

unreinforced

perfect bonding

settlement (mm)

-110 -140

10

-50

-20

-80

n→ ∞

n=25

perfect bonding

L = 12m

-200 settlement (mm)

n=16

Fig. 14. Evaluations of the foundation settlement for a prescribed load level (20 MN).

unreinforced -80 -110

n=9

0

L = 16m

-170 -200 settlement (mm)

25

0

5

Fig. 13. Load settlement curves for different pile lengths and densities (multiphase simulations).

20

25

30

total load (MN) 8m/M

-100

12 m / M 16 m / M

-150

8m/D 12m / D

-200

the number n of piles under the raft have been varied. The curves corresponding to the un-reinforced foundation as well as, for each pile length, to the case of perfect bonding (where all elastic and strength interaction parameters are set equal to infinity), have also been drawn in the same figures. Several comments ought to be made upon examining these results. First, as could be expected, for a given number of piles, the settlement is a decreasing function of the pile length. Second, for a fixed value of the pile length, the number of piles (or more precisely the pile density, since it should be kept in mind that the reinforcement volume fraction and thus the elastic parameter a are kept constant), has a quite significant influence on the settlement reduction. It thus appears that reinforcing the foundation by a large number of small diameter piles is more efficient in terms of settlement reduction than using a reinforcement scheme involving relatively few large diameter piles. The limit theoretical configuration would correspond to an infinite number of piles (n ? 1) of infinitely small diameter, thus recovering the perfect bonding case where the interaction parameters are assigned very large values (upper curves in Fig. 13).

15

-50

n=9 16

10

n =9

16 m / D -250

settlement (mm) 0

0

5

10

15

-150 -200 -250

25

30

total load (MN)

-50 -100

20

8m/M 12 m / M 16 m / M 8m/D

n =25

12 m /D 16 m / D

settlement (mm) Fig. 15. Comparison of load-settlement curves computed from the direct (D) and multiphase (M) approaches.

35

E. Bourgeois et al. / Computers and Geotechnics 46 (2012) 26–38

dation, than, the direct simulation. One should however remain cautious about the latter comment, since the results of the direct numerical simulations appear to be rather sensitive to the finite element discretization, notably in the vicinity of the pile tips, making their numerical accuracy somewhat questionable.

and density is clearly apparent from this Figure. Indeed, the predicted value of the settlement, equal to 165 mm for 3  3 8 m-long piles drops to 50 mm for 5  5 16 m-long piles, that is a reduction factor larger than 3. 5.3. Multiphase vs. direct approach

5.4. Pile force distributions For comparison and validation purposes, the same load-settlement curves have been drawn, derived from ‘‘direct’’ numerical computations using finite element models such as that sketched in Fig. 11, where the mesh has been refined around each individual pile in order to capture the pile-ground interactions with as much accuracy as possible. A fairly good agreement may be observed from Fig. 15 between the two numerical predictions, with in most cases a small relative difference in terms of settlement, never exceeding 20% in the worst cases. It should be noted that the agreement is better for short piles than for long piles. In all the situations considered in the analysis, the multiphase model leads to slightly lower values of the settlement, and thus to a somewhat stiffer global behaviour of the foun-

n=9, L=12m

n=9, L=8m 0

-0.2

-0.4

-0.6

0

0

-0.2

axial force (MN)

4

n=9, L=16m

-0.6

4 6

8

8 center (D)

10

center (M)

12

corner (D)

14

corner (M)

depth (m)

-0.1

-0.2

axial force (MN)

4

center (D)

10

center (M)

12

corner (D)

14

corner (M) depth (m)

12 14 depth (m)

16

0

-0.2

-0.1

center (D)

center (M)

12

center (M)

corner (D)

14

corner (D)

corner (M)

16 depth (m)

corner (M)

n=16, L=16m

-0.4

axial force (MN)

2 4

6 center (M)

12 14 depth (m)

-0.4

-0.6 axial force (MN)

2 4 6

center (D)

8

10

center (M)

10

center (D)

12

center (M)

14

corner (D)

12

corner (D)

14

corner (M) depth (m)

16

corner (M)

16 depth (m)

0

-0.1

-0.2

n=25, L=16m

-0.3

-0.4

2 4

axial force (MN)

corner (D)

12

corner (M)

14 16 depth (m)

-0.1

-0.2

-0.3

-0.4 axial force (MN)

2 4 6 8

8 10

0 0

6 center (D)

10

-0.2

8

-0.3 axial force (MN)

8

0

-0.6

0

4

axial force (MN)

4

n=25, L=12m

-0.2

0 2

2

0

n=25, L=8m 0

-0.8

10

6

6

-0.6

center (D)

0

2

-0.4

0

8

10

-0.4

0

16

-0.2

n=16, L=12m

-0.3

8

0

6

n=16, L=8m 0

-0.8 axial force (MN)

2

6

16

-0.4

0

2

16

A quite significant output of the numerical simulations concerns the distribution of forces along the piles. The charts presented in Figs. 16 and 17 refer to such predicted distributions for a pile placed at the centre and at the corner of the raft, respectively, in the different configurations obtained from varying the pile length and density. Fig. 16 displays such results obtained from the direct and multiphase numerical simulations for a load level equal to 40% of the ultimate bearing capacity. Starting from the pile tips, where the axial force is significantly different from zero, the axial force is increasing when moving upwards, first in a linear way in the bottom part, then more progressively in the upper part up to the pile

center (D)

10

center (M)

12

corner (D)

14

corner (M)

16

center (D) center (M) corner (D) depth (m)

corner (M)

Fig. 16. Pile force distributions at 40% of the ultimate load bearing capacity (in the legend (M) refers to the multiphase approach and (D) to the standard finite element simulations).

36

E. Bourgeois et al. / Computers and Geotechnics 46 (2012) 26–38

n=16, L=8m 0

n=16, L=12m

-0.2

0

-0.4

0

axial force (MN)

2 4

8

10

10

12

12

16

axial force (MN)

-1

axial force (MN)

5

10

M

14

M

D

depth (m)

depth (m) n=16, L=12m

n=16, L=8m 0

-0.2

-0.4

0

0

axial force (MN)

2 4

-0.2

n=16, L=16m

-0.4

-0.6

0

-0.5

-1

0

0

axial force (MN)

2 4

axial force (MN)

5

6

6 8

8

10

10

12

12

M

14

D

16

10

M

14

M

D

D

15

depth (m)

16

depth (m)

depth (m)

n=25, L=8m 0

D

15

16

depth (m)

-0.5

0

4

8

D

0

-0.6

2

6

M

n=16, L=16m

-0.4

0

6

14

-0.2

-0.1

n=25, L=12m

-0.2

-0.3

0

axial force (MN)

2

0

-0.2

-0.4

0

n=25, L=16m

axial force

2

0

-0.6

6

6

6

8

8

10

10

M

D

14

16

16

depth (m)

8 10

12

14

-0.6

axial force (MN)

2 4

12

-0.4

0

4

4

-0.2

M

D

12

M

D

14

depth (m)

depth (m)

16

Fig. 17. Pile force distributions at 90% of the ultimate load bearing capacity ((M): multiphase approach; (D): standard finite element simulations).

-0.5

-1

n=25, L=16m

n=25, L=12m

n=25, L=8m 0 0

-1.5

0

-0.5

axial force (MN)

4 6

2 4 6

axial force (MN)

2 4 6

8

8

8

10

10

10

12

12

general model

12 14

perfect bonding

16

depth (m)

14 16

depth (m)

-0.5

-1

-1.5

0

0

2

0

-1

general model perfect bonding

14

axial force (MN) general model perfect bonding

16

depth (m)

Fig. 18. Comparison of pile force distributions predicted by the general multiphase model and those based on the perfect bonding assumption.

E. Bourgeois et al. / Computers and Geotechnics 46 (2012) 26–38

cap. This kind of profile can be attributed to the fact that the pile tip resistance is fully mobilized (p = p0) as well as plasticity of the shaft interaction in the lower zone (I = I0), hence the linear variation of the axial force in this region (see equilibrium equation (9)). As can be seen from these charts, the axial force distributions obtained from the direct and multiphase simulations, respectively, appear to be relatively close to each other. They are almost coincident in the lower part for, notably for a small pile density (n = 9), while the multiphase simulation produces slightly larger values when approaching the pile cap, which is consistent with the above mentioned global stiffer response predicted by the multiphase calculations (see Section 5.3.). As regards the influence of the position of the pile under the raft, it is clearly apparent from Fig. 15 that the axial force distributions in the lower part are the same for a pile located at the centre or at a corner of the raft, whereas the central piles undergo smaller axial forces than the corner piles in the upper part. The discrepancy between the two profiles becomes all the more important as the pile length and density increase. These observations are confirmed by Fig. 17 which displays the same kind of results for a load level equal to 90% of the ultimate bearing capacity. Indeed, the pile force distributions are the same, whatever the location of the pile. They correspond to a linear variation of the form given by Eq. (33) where both the pile tip and pile shaft resistance along the entire pile length, are mobilized. Meanwhile, the direct and multiphase simulations produce almost coincident predictions, with the possible exception of loosely packed short piles (n = 9, L = 8–12 m). This is again to be related with a much better agreement between the corresponding load-settlement curves when approaching the ultimate bearing capacity of the pile raft foundation. So as to further illustrate the crucial role played by the interaction parameters in the design of piled raft foundations, Fig. 18 shows the completely different pile force distributions predicted by the same multiphase model where perfect bonding is assumed by assigning very large values to both stiffness and strength interaction parameters. Such distributions are clearly overestimated from making such an assumption, notably at the pile tip level, resulting in a much stiffer global response of the foundation settlement and unrealistic value of the ultimate bearing capacity (Fig. 13).

6. Conclusion A quite significant improvement of the multiphase model of reinforced soils, as specifically dedicated to the design of piled raft foundations, has been achieved through the formulation, and subsequent development, of a specific interaction law pertaining to the ‘‘pile tip resistance’’. The approach has been successfully validated from comparing the results of the multiphase simulations, incorporating such an innovative feature, to those obtained from standard, comprehensive, but highly computational time consuming, finite element calculations. Indeed, provided that the interaction stiffness and yield strength parameters have been properly identified, the multiphase-based calculation model provides reliable predictions both in terms of global settlement and pile force distributions, in a much easier and quicker way than the direct, but fairly cumbersome, finite element approach. This calculation model and the computational tool which has been derived from it, pave the way for the engineering design of piled raft foundations on a rational basis. Meanwhile, the present paper has established some important guidelines concerning the decisive question of how to determine

37

the multiphase interaction parameters from the pile reinforcement material and geometric characteristics. Thus, making use of simple dimensional analysis arguments, it has been shown in particular that the ‘‘pile shaft’’ and ‘‘pile tip’’ interaction elastic coefficients can be easily related to the pile density, so that the identification procedure of those parameters can be greatly simplified. Likewise, as regards the ‘‘pile shaft’’ interaction strength, a quite simple relationship can be established for a purely cohesive soil between the value of this parameter and the pile diameter as well as surrounding soil undrained shear strength. On the other hand, the situation seems to be mode complex concerning the ‘‘pile tip’’ interaction strength parameter, where no clear relationship with the pile density has been revealed so far, while some geometric characteristics such as the pile length seem to be involved in such a strength characteristic. As a means to clarify the latter point, a careful systematic parametric study should be undertaken in much the same way as it has been done for instance, in the context of elasticity, for piled embankments [7,14], where interaction at the ‘‘pile head’’, instead of the ‘‘pile tip’’, was considered. The objective of such an analysis, which could then be extended to frictional soils, would be to produce simple formulas expressing the value of the interaction strength parameters as a function of geometric as well as strength characteristics of the soil and soil-pile interface. References [1] Ashour M, Norris GM, Elfass S, Al-Hamdan AZ. Mobilized side and tip resistances of piles in clay. Comput Geotech 2010;37(7–8):858–66. [2] Baker CN, Azam T, Joseph LM. Settlement analysis for the 450 m KLCC towers. vol. 40. ASCE special Technical Publications; 1994. p. 1650–71. [3] Bennis M, de Buhan P. A multiphase constitutive model of reinforced soils acounting for soil-inclusion interaction behavior. Math Comp Modell 2003;37:469–76. [4] Bourgeois E, Hassen G, de Buhan P. Finite element simulations of the behavior of piled-raft foundations using a multiphase model. Int J Numer Anal Meth Geomech 2012. http://dx.doi.org/10.1002/nag.2077. [5] de Buhan P, Bourgeois E, Hassen G. Numerical simulation of bolt supported tunnels by means of a multiphase model conceived as an improved homogenization procedure. Int J Numer Anal Meth Geomech 2008;32:1597–615. [6] de Buhan P, Hassen. Macroscopic yield strength of reinforced soils : from homogenization theory to a multiphase approach. C R Mécanique 2010;338:132–8. [7] Cartiaux F-B, Gellee A, de Buhan P, Hassen G. Modélisation multiphasique appliquée au calcul d’ouvrages en sols renforcés par inclusions rigides. Revue Française de Géotechnique 2007;118:43–52. [8] Chow YK. Analysis of vertically loaded pile groups. Int J Numer Anal Meth Geomech 1986;10:59–72. [9] Chow YK, Yong KY, Shen WY. Analysis of piled raft foundations using a variational approach. Int J Geomech ASCE 2001;1(2):129–47. [10] Clancy P, Randolph MF. An approximate analysis procedure for piled raft foundations. Int J Numer Anal Meth Geomech 1993;17:849–69. [11] Comodros EM, Papadopoulou MC, Rentzeperis IK. Pile foundation analysis and design using expremimental data and 3-D numerical analysis. Comput Geotech 2009;36:819–36. [12] Griffiths DV, Clancy P, Randolph MF. Piled raft foundation analysis by finite elements. In: Geomech, Beer, Booker, Carter (Eds.), Proc comp meth adv. Rotterdam: Balkema; 1991. p. 1153–57. [13] Hassen G, de Buhan P. Elastoplastic multiphase model for simulating the response of piled raft foundations subject to combined loadings. Int J Numer Anal Meth Geomech 2006;30:843–64. [14] Hassen G, Dias D, de Buhan P. Multiphase constitutive model for the design of piled-embankments: comparison with three-dimensional numerical simulations. Int J Geomech 2009;9(6):258–66. [15] Humbert P, Dubouchet A, Fezans G, Remaud D. CESAR-LCPC: a computation software package dedicated to civil engineering uses. Bull Laboratoires des Ponts et Chaussées 2005;256–257:7–37. [16] Katzenbach R, Arslan U, Gutwald J, Holzhaüser J, Quick H. Soil structure interaction of the 300 m high commerzbank tower in Frankfurt am Main – measurements and technical studies. In: 14th Int conf soils mechanics and foundations engineering; 1997. vol. 2. p. 1081–84. [17] Lee J, Kim Y, Jeong S. Three-dimensional analysis of bearing behaviour of piled raft on soft clay. Comput Geotech 2010;37(1–2):103–14. [18] Liang F-Y, Chen L-Z. A modified variational approach for the analysis of piled raft foundations. Mech Res Commun 2004;31:593–604.

38

E. Bourgeois et al. / Computers and Geotechnics 46 (2012) 26–38

[19] Reul O. Numerical study of the bearing behaviour of piled rafts. Int J Geomech, ASCE 2004;4(2):59–68. [20] Shen WY, Teh CI. A variational solution for downdrag force analysis of pile groups. Int J Geomech, ASCE 2002;2(1):75–91. [21] Small JC, Zhang HH. Behaviour of piled raft foundations under lateral and vertical loading. Int J Geomech, ASCE 2002;2(1):29–45. [22] Smith, Wang. Analysis of piled rafts. Int J Numer Anal Methods Geomech 1998;22:777–90.

[23] Sudret B, de Buhan P. Multiphase model for inclusion-reinforced geostructures, Application to rock-bolted tunnels and piled raft foundations. Int J Numer Anal Methods Geomech 2001;25:155–82. [24] Vetter K. Untersuhungen zum Traverhalten der Kombinierten Pfahlplattengründung des Messeturms in Frankfurt am Main auf des Basis von Messungen und Numerischen Computersimulationen. TU Darmstadt – Diplomarbeit;1998. p. 190.