Finite elements for cohesive-frictional material

Finite Elements in Analysis and Design 47 (2011) 784–795

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

Finite elements for cohesive-frictional material G.B. Muravskii Faculty of Civil Engineering, Technion—Israel Institute of Technology, Technion, 32000 Haifa, Israel

a r t i c l e i n f o

a b s t r a c t

Article history: Received 7 November 2009 Received in revised form 22 November 2010 Accepted 25 February 2011 Available online 21 March 2011

New finite elements comprising parts (rods and spiral springs) with independent behaviors, which can be determined from the corresponding independent experiments – shear, deviatoric loading, hydrostatic loading – are suggested for isotropic material in 2D and 3D cases and applied for solving some problems for cohesive-frictional material. The material is assumed to be non-linear with the stress– strain relationship of hyperbolic type and failure conditions determined by Mohr–Coulomb law. The difficulty of the problems consists in the fact that the limit stresses themselves depend on unknown stress distribution, so a stress–strain equation cannot be written in the explicit form before solving the problem. This difficulty is avoided by iterations: a suitable initial stress distribution (e. g. from the corresponding elastic problem) is supposed and, knowing the material properties, the new stress distribution (appearing as a result of slow gradual load application) can be found which is taken as the initial for the next iteration, and so on. The following two problems are considered: failure of a horizontal layer under action of a uniform pressure applied to a rectangular area on the layer surface and the problem of slope stability. & 2011 Elsevier B.V. All rights reserved.

Keywords: Finite elements Coulomb–Mohr law Hyperbolic law Slope stability

1. Introduction An idea of application of finite elements built from rods and springs is rather attractive since we can deal with one-dimensional objects and therefore taking into account complicated properties of materials in a relatively simple way provided a sufficient number of simple parts is taken into consideration. This gives a possibility of satisfying all energy requirements related to the deformable material. Actually a treatment of elastic media as a set of particles, which interact as being connected to each other with elastic springs, was considered by Poisson [1] who found that the treatment leads to an elastic isotropic material with Poisson’s ratio n ¼1/4. In the book by Rzhanitsin [2] two cases of a plane structure of elastic rods are considered with square and triangle forms of a periodic part of the structure. It was found that the structures are equivalent to an elastic isotropic material with Poisson’s ratio n ¼1/3 for plane stress conditions. In the case of plane strain conditions, these structures lead to Poisson’s ratio equal to 1/4. Another modeling 2D and 3D continuum by a set of one-dimensional elements was suggested in works of Streeter and Wylie [3] and Papadakis [4]. The rod elements can transmit shear and pressure waves and at interior nodes of the latticework, specific transfer elements receive and transmit the waves. Note that this modeling implies zero value of Poisson’s ratio for the material under consideration. In the paper by Muravskii [5], one-

E-mail address: [email protected] 0168-874X/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2011.02.009

dimensional parts of the finite elements are applied in such a way that, for a linearly elastic isotropic body, arbitrary values of Poisson’s ratio can be approached; however, whole independence in the behavior of separate components of the finite element was not achieved. In the present paper, new finite elements are suggested for 2D and 3D problems. Stipulating isotropy and small strains, the properties of the component of the element can be determined from independent experiments—shear, deviatoric loading, hydrostatic loading. Note that the shear and deviatoric experiments are assumed to be performed in the presence of an additional confining pressure. In geomechanics, it is important to study a non-linear material with failure conditions determined by Mohr–Coulomb law. The difficulty of the appearing problem consists in the fact that the limit stresses depend themselves on the unknown stress distribution, so a stress–strain equation cannot be written in the explicit form before solving the problem. This difficulty is avoided by iterations: a suitable initial stress distribution is assumed, and knowing the supposed material properties, the new stress distribution (occurring as a result of slow gradual load application) is found which is taken as the initial for the next iteration, and so on. The stress–strain relationship for shear and deviatoric loading is taken as hyperbolic [6], and the relationship between the hydrostatic pressure and the relative volume variation is assumed as linearly elastic with the values G0 and n0 of shear modulus and Poisson’s ratio. Note that this assumption can be changed without significant difficulties. At each iteration, the corresponding static problem is treated dynamically provided that the load is applied

G.B. Muravskii / Finite Elements in Analysis and Design 47 (2011) 784–795

where p is the pressure (tension) and D is the relative volume change of the element. Although relation (1) is written as a function, it is assumed that this relationship can take into account the history of deformation (as in the case of the hysteretic behavior) and derivatives of D. In the case of a medium with pores containing an additional component, gas or water, Eq. (1) relates to the main material (skeleton), and pore pressure, pp, should be taken into account additionally, so the whole pressure is equal to p þpp. Introducing pore pressure into calculations leads to well-known conception of ‘effective’ values of stresses in the main material. We apply the following important assumption related to the rods entering the elements: they have zero stiffness for axial deforming and are able to transfer to the corresponding nodes of the element the cross forces appearing due to moments in the spiral springs.

slowly achieving the final values at the end of some assigned time interval (0, t0). These dynamic problems are solved with the help of the method of constant average acceleration at each time step, which itself requires (internal) iterations because of non-linearity of the problems. The suggested method works well under the condition that the stress–strain state of the considered system is not very close to the state of failure when a small variation in stresses can result in significant changes in limit stresses and therefore in the material behavior. The problems requiring the consideration of material failure can be treated using an extrapolation: for clarification of failure state, the stresses are used, which correspond to a state possibly close to the failure state but allowing yet the convergence of the above iteration process. This treatment is applied in an example relating to slope stability in which shear and deviatoric stresses obey the hyperbolic law with the Mohr–Coulomb failure condition, and the hydrostatic deformations are linearly elastic. The 3D finite elements are demonstrated in a problem with a horizontal layer under the action of a uniform normal pressure applied to a rectangular area on the layer surface.

2.1. Finite elements for plane problems Consider a volume of an isotropic deformable body of dimensions a  b  d (d is the measure in orthogonal direction to the plane of the element), which is modeled by the finite element of Fig. 1a. Suppose that in a shear deformation experiment with shear stresses txy on the faces of the element we obtain a relationship between shear stress t ¼ txy and shear strain g:

2. Finite elements suitable for description of behavior of materials with complicated properties Consider rectangular elements (for plane problems) and elements in the form of the rectangular parallelepiped (for 3D problems) consisting of rods and spiral springs (Fig. 1a and b). The adjacent diagonal rods and adjacent vertical and horizontal rods are connected with spiral springs, which generate a moment depending on variation of the angle between the rods. Such springs are assumed in all the faces of the 3D element of Fig. 1b (where only the springs for one of the parallelepiped faces are shown), so we have 30 springs for each 3D element. We will denote rectangle apexes for each face in the manner as shown in Fig. 1b using indices for identification of faces: 1, 4 for the faces parallel to the plane (x, y); 3, 6 for the faces parallel to the plane (y, z); and 2, 5 for the faces parallel to the plane (x, z). The apex denoted as A represents simultaneously the apexes A1, A2, A3 in the corresponding rectangles lying in the planes 1, 2, 3; the apex B serves as the apexes B1, B2, A6; the apex C serves as the apexes C1, B5, D6; and the apex D serves as the apexes D1, D3, A5. Analogously, E is for D2, A4, B3; F is for C2, B4, B6; G is for C4, C5, C6; and K is for D4, D5, C3. This notation is suitable for organization of calculations. Suppose that interior of the elements is filled with a deformable material of the kind of compressible fluid (possibly able to withstand also tension), which reacts only on volume variation: p ¼ pðDÞ

785

t ¼ tðgÞ

ð2Þ

As above, this relationship can be understood more broadly than a simple function. Note that the shear deformation g is equal to the angle change j of the initially right apex angles of the rectangle in Fig. 1a, and the angle between diagonals does not change (in the linear approach, i.e. neglecting powers of j greater than one in the angle change). It should be noted that the relationship (2) can depend on normal stresses. When considering below the application of Mohr–Coulomb law along with the hyperbolic stress–strain relation we see how the normal stresses influence the relationship (2). Consider (Fig. 2) equilibrium of one of element’s nodes, say, C under action of the forces stemming from stresses t (horizontal tad/2 and vertical tbd/2) and from moments m (horizontal 2 m/b and vertical 2 m/a), which act oppositely to the forces from stresses. Equating the corresponding forces from stresses and from moments we come to the following relationship between moments m in the spiral springs located at the apexes of the element and the shear stress: mðjÞ ¼

tðjÞabd

ð3Þ

4

To obtain the constitutive equation for the moment M in the spring, which connects the diagonal rods, consider (Fig.3) the

ð1Þ

a

C

D

4 C mC1

D b

G

5

K

mD1

6 b

3 M1

y

F

A

y x

E

z

a

c

A

B

x

mA1

1

Fig. 1. Finite elements for 2D problems (a) and 3D problems (b).

2

mB1

B

786

G.B. Muravskii / Finite Elements in Analysis and Design 47 (2011) 784–795

Dm

τ

elements boundaries will be

m C

1 ðmB þ mD mA mC Þ ð8Þ abd Here the moments generated in the spiral springs in points A, B, C, and D (Fig. 1a) are taken as positive if the corresponding angle increases. The averaged normal stresses on boundaries of the element arise from the hydrostatic stress p and the deviatoric addition due to the action of the moment M. We obtain (positive stresses correspond to pressure):

txy ¼

b τ

A m

τ

τ a

m B

Fig. 2. Stresses and moments in spiral springs under shearing.

sx ¼ p2 q

D

b

q

C

M

q

A

B

Fig. 3. Stresses and moments in spiral spring under deviatoric loading.

deviatoric loading of the element with the pressure q in the vertical direction and the tension q in the horizontal direction; corresponding strains will be  e and e. Note that this loading can be meant as an additional loading to the previously existing hydrostatic pressure, thus the material under total load can be free of tension. For the spiral spring connecting the diagonal rods, the increase of the upper angle will be (in the linear approach) 4abe a2 þ b2

ð4Þ

Assume that the relationship (as above, in the wide sense) between q and e is known: q ¼ qðeÞ

ð5Þ

Consider equilibrium of one of element’s nodes, say, C under action of the forces stemming from stresses q (horizontal qbd/2 and vertical qad/2) and from the moment M, M/(a2 þb2)1/2, whose direction is orthogonal to the diagonal AC. The horizontal and vertical components of this force act oppositely to the corresponding forces representing stresses. The equilibrium equations lead to the following relationship between moment M in the spiral spring connecting the diagonal of the element and the deviatoric stress q: M¼

qða2 þ b2 Þd 2

 2  a þb a þb2 Mðjd Þ ¼ dq jd 2 4ab

M ða2 þ b2 Þd

ð10Þ

ð6Þ

A more accurate treatment of the normal stresses leads to linearly distributed supplements with zero values at the centrum of the corresponding side of the element contour, which results in the moments Mx and My of normal stresses (Fig.4). The equilibrium equations written for the nodes of the finite elements lead to the following expressions for the ‘‘bending’’ moments (in addition to (8)–(10)): Mx ¼  My ¼

ð7Þ

The expression in the brackets represents e through jd according to (4). In the considered case, the value p in the relationship (1) should be meant as the mean normal stress in the (x, y) plane. Eqs. (3) and (7) allow us to express forces, with which the element acts on its nodes, through the node displacements and develop rather simple algorithm for solving different problems of deformable body. If a problem is solved and for all elements the pressure p and moments M, mA, mB, mC, mD have been found, stresses on the faces of the elements can be determined by uniform distribution along the corresponding sides, and the node forces appearing from the moments in the spiral springs with the corresponding balancing by stresses. The shear stresses on the

ð11Þ ð12Þ

An important advantage of the suggested finite elements consists in the fact that their constituents (the springs in the vertexes of the rectangle, the spring connecting the diagonals, and the filling) have independent behaviors, which are determined by the shear stress–shear strain relationship (the springs in the vertexes), the deviatoric normal stress–deviatoric normal strain relationship (the spring connecting the diagonal rods), and the p–D relationship for the filling. Thus, we can consider the constituents as one-dimensional parts with known behavior. 2.1.1. Linearly elastic isotropic body Consider as an example a linearly elastic body with shear modulus G and Poisson’s ratio n. The obtained results can serve as an approbation of the transformations made. For modeling the body using finite elements of Fig. 1a in the plane strain conditions, we apply Eqs. (1), (4), and (7) in the form p¼

m¼

G G ðeAD þ eBC þ eAB þ eDC Þ D¼ 12v 2ð12vÞ

ð13Þ

Gj abd 4

M¼G

2

mA þ mB mD mC 2

mA þ mD mB mC 2

or 2

ð9Þ

sy ¼ p þ2

q

a

jd ¼

M ða2 þb2 Þd

ð14Þ

ða2 þ b2 Þ2 jd d 4ab

ð15Þ

where in Eq. (13) the mean stress in (x, y) plane, i.e. p¼ (sx þ sy)/2, is expressed through rod deformations. Note that in the

My

σy

Mx

Mx

σx

σx

My

σy

Fig. 4. Additional bending moments appearing at element faces.

G.B. Muravskii / Finite Elements in Analysis and Design 47 (2011) 784–795

considered case for deviatoric loading, Eq. (5) has the form q(e)¼ t(2e)¼2Ge. For plane stress conditions, Eq. (13) is changed to p¼

Gð1 þ nÞ ðeAD þ eBC þ eAB þ eDC Þ 2ð1vÞ

ð16Þ

For testing the new finite element consider a plane stress problem for a plate (deep beam) having length L¼2.4 m, height H¼1.6 m and width¼0.1 m (Fig. 5). The beam is loaded along its upper side by the vertical uniformly distributed load P0 ¼ 10 N/m, which is in equilibrium with the shear stresses applied to the left and right boundaries of the beam at x ¼ 7L/2:   3P L H2 y2 txy ¼ 7 03 ð17Þ 4 dH An exact solution of this problem is constructed in [7] by applying additional self-balanced normal stresses to the left and right boundaries:   6P 2y3 H2 y ð18Þ  sx ¼  03 10 3 dH The problem has been solved by finite element method using the new finite elements and also traditional rectangular finite elements, in which to linear distribution of displacement within elements the term with xy is added. The results for G ¼1000 Pa, n ¼ 1/3 are compared with exact solution from [7] in Table 1. In the first row are results for dividing the height H to 16 parts and the half of the length, L/2, to 12 parts (192 square elements for half of the beam), in the second row are results corresponding to 4  16 rectangular (a ¼3b) elements for the half of the beam. The vertical displacements at the points (0,0), (0.6,0) and stresses sx at the point (0.15, 0.75) in the center of corresponding elements (i.e. without consideration of the above-mentioned ‘‘bending’’ moments) are presented. The data for the traditional rectangular elements are given in brackets. For the considered finite elements, the accuracy is of the same order. The new elements lead to a more accurate solution as compared to traditional finite elements in the case of square elements whereas a better accuracy corresponds to the traditional finite elements in the case a ¼3b. It is advisable to use the new elements with a and b close to each other. Though, the above linearly elastic problem is considered here only for an illustration, the suggested new elements are y

P0 = 10N/m

787

dedicated mainly for studying materials with complicated behavior. 2.1.2. Non-linear material with limit stresses obeying Mohr– Coulomb law (plane strain conditions) Let the material for infinitesimal small strains be linearly elastic and isotropic with Poisson’s ratio n0 and shear modulus G0. Assume initially that for an element in equilibrium conditions, the stresses sx and sy (in the center of the element) and therefore the pressure p in the filling and the deviatoric normal stress q are known. From Fig. 6 the following maximum values for q and txy can be obtained qmax ¼ R ¼ c cosðfÞ þp sinðfÞ

tmax ¼

ð19Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2 q2

ð20Þ

where c and f are cohesive shear stress and angle of friction, respectively; for a material comprising water or gas in pores with an own pressure they are meant as effective values. It should be noted that acting in the frame of method of iteration mentioned above, the value of q in an element can be more than R for the taken stress distribution at an iteration start. In this case, the element will be in the state of failure, and in calculations one can accept q ¼0.999R with corresponding small value of tmax according to (20). This means a state very close to failure. In subsequent iterations the state of the element can be changed as a result of the change in stresses. Emphasize that the suggested method considers two causes of element failure connected with (19) and (20). Knowing the limiting values for q and txy the constitutive equations for deviatoric and shear stresses can be constructed. In the case of hyperbolic law, for example, we obtain: qðeÞ ¼ qmax

tðgÞ ¼ tmax

e

ð21Þ

qmax =ð2G0 Þ þ 9e9

g tmax =G0 þ9g9

ð22Þ

In this treatment, a hidden assumption is applied that material properties are fixed, being defined by the final stress state. If an external load, which corresponds to this state, will gradually increase from zero to its actual value then Eqs. (21) and (22) should lead to the stress state used for the material properties definition. These considerations form a basis for above-mentioned

H = 1.6m

τ xy

τ

σx O

x

τmax τxy

R

φ

L = 2.4m Fig. 5. Plane stress problem for deep beam.

c

p

σ p+q

Fig. 6. Determination of limit values for deviatoric stress q and shear stress txy.

Table 1 Comparison of results for suggested finite elements, traditional finite elements (in brackets) and exact solution. Elements for half of beam

192 (a ¼b) 64 (a¼ 3b) Exact solution

Stress sx (Pa)

Vertical displacements (m) at x ¼0, y ¼0

at x ¼0.6, y ¼0

at x¼ 0.15, y ¼0.75

 0.09619 ( 0.09611)  0.09339 ( 0.09427)  0.09640

 0.07035 ( 0.07029)  0.06844 ( 0.06891)  0.07052

168.75 (168.5) 160.77 (164.35) 168.81

788

G.B. Muravskii / Finite Elements in Analysis and Design 47 (2011) 784–795

method of iteration designated for determination of unknown stresses in the non-linear material, whose properties are influenced themselves by stresses. In the case of a material with hysteretic behavior [8,9], relationships (21) and (22) can serve for definition of the corresponding backbone curves, which allow to study vibrations following the static loading applied previously. When solving problems relating to soil mechanics initial shear modulus G0 should be taken in turn as a function of hydrostatic pressure initially unknown. The relation of the type pffiffiffi ð23Þ G0 ¼ a þ b p is applied where the value of p is assumed to be zero in the case of its negative values. Emphasize that for determination of the initial shear modulus G0 for the hyperbolic or another law we use the internal pressure p corresponding to the actual stress–strain condition (see above considerations relating to Eqs. (21) and (22)). As an example of Eq. (23), consider the following equation [10,11] suitable for many undisturbed cohesive soils, as well as sands: G0 ¼ 1230

ð2:973eÞ2 ðOKRÞK s0 1=2 1þe

ð24Þ

in which e¼void ratio, OKR¼overconsolidation ratio, s0 ¼mean principal effective stress, and both s0 and G0 are in pounds per square inch (psi). The value of K depends on the plasticity index, PI, of soil and lies between 0 and 0.5. Taking in (24) e¼1, K ¼0, s0 ¼ 2p(1þ n0)/3Ep and going from psi to Pa for G0 and s0 (or p) we obtain the following assessed value for the coefficient b E200 kPa, which will be applied in the examples. According to (24) the value of a equals zero. In the following calculations, we take for a a relatively small value (10 kPa) instead of zero. The uncertainty relating to the coefficients a and b does not play an important role, since noticeable changes in the values of a and b do not influence practically the results concerning stress distribution or failure conditions. Note that in the case of the plane strain problem for a horizontal linearly elastic isotropic layer the depth distribution of normal stresses is not influenced at all by stiffness of the layer and is determined only by Poisson’s ratio n: p¼

w0 h 2ð1nÞ

ð25aÞ

sy ¼ w0 h

ð25bÞ

nw h sx ¼ 0 1n

ð25cÞ

where w0 is the unit weight and h is depth of a current point. In the case of hyperbolic law, Eqs. (7), (21), (4) and (3), (22) allow us to write the following constitutive equations for M(jd) and m(j): Mðjd Þ ¼

mðjÞ ¼

a2 þb2 jd dqmax 2 2abqmax =ðða2 þ b2 ÞG0 Þ þ 9jd 9

abd j tmax 4 tmax =G0 þ9j9

ð26Þ

ð27Þ

For the third constituent of the finite element, i.e. the filling reacting only on volume change, we assume at this stage the linearly elastic behavior, which is determined by the parameters G0 and n0. According to (13) we obtain for plane strain problems p¼

G0 ðeAD þ eBC þ eAB þ eDC Þ 2ð12v0 Þ

Coulomb, on unknown final stresses. In the method developed below the iterations are suggested: a feasible distribution of stresses (for example as in Eq. (25) in the case of the self-weight action) is assumed as an initial approach; after solving the problem we improve the initial approach (using Eqs. (9) and (10), which determine the normal stresses). After several iterations, the initial and final stress distributions practically coincide, which shows that the required stress distribution has been found. Convergence of the method becomes worse when the sought-for final equilibrium state is near to failure; this is demonstrated below.

2.2. Finite elements for 3-D problems First, we consider shear deformations of the 3-dimentional element shown in Fig. 1b. The corresponding moments in spiral springs are denoted with subscripts 1 (4), 2 (5), and 3 (6) for planes parallel to (x, y) plane, (x, z) plane, and (y, z) plane, respectively (Fig. 1b). For all moments corresponding to rectangle apexes we obtain identical constitutive equations, thus the indices can be omitted mðjÞ ¼

8

ð29Þ

However the behavior t(j) can be different for the mentioned planes and actually we use in calculations three different equations for moments relating to the planes 1 (4), 2 (5), and 3 (6). This will be explained below when considering Mohr–Coulomb condition of failure. More complicated are equations relating to springs, which connect the diagonal rods. Consider deviatoric loading of the element of isotropic material with a stress q (tension) in the direction along x-axis and the stress q in the direction of y-axis; corresponding strains will be e,  e, 0 in directions of axes x, y, z, respectively. Only deforming of springs connecting diagonal roads occurs. Moments for parallel faces behave identical in the considered loading, so one can deal only with moments M1, M2, M3. For the finite element, the three equilibrium equations for the node C (Fig. 1b) have the form: M1 M2 qbc bþ 2 c¼ 4 a2 þ b2 a þ c2

ð30Þ

M1 M3 qac aþ 2 c¼ 4 a2 þ b2 b þ c2

ð31Þ

M2 M3 a¼ 2 b a2 þ c2 b þc2

ð32Þ

The last equation corresponds to the absence of external forces in the z-direction. As it is seen, this equation follows from (30) and (31), i.e. only two equations are independent. The angle change for the spiral springs generating moment M1 is expressed similarly to Eq. (4):

jd1 ¼

4abe a2 þb2

ð33Þ

For angle changes corresponding to moments M2 and M3, we have (using the absence of displacements in the z-direction):

ð28Þ

As has been mentioned above, the most difficulty arising when solving the considering problems consists in the fact that the behavior of the material during loading depends on the material strength, which in turn depends, according to the law of Mohr–

tðjÞabc

jd2 ¼

2ace a2 þc2

ð34Þ

jd3 ¼

2cbe c2 þb2

ð35Þ

G.B. Muravskii / Finite Elements in Analysis and Design 47 (2011) 784–795

789

Further, consider deviatoric loading parallel to the plane (y, z). The following equations analogous to (30)–(32) can be obtained: M2 M3 qab aþ 2 b¼ 4 a2 þc2 b þc2

ð36Þ

M1 M3 qac aþ 2 c¼ 4 a2 þb2 b þ c2

ð37Þ

M1 M2 b¼ 2 c a2 þb2 a þc2

ð38Þ

with the angle changes

jd3 ¼

4cbe c2 þ b2

ð39Þ

jd2 ¼

2ace c2 þ a2

ð40Þ

jd1 ¼

2abe a2 þ b2

ð41Þ

Fig. 7. Exact and approximate (dotted line) definition of moment in spring connecting diagonal roads.

Note that the angle argument for M2 is the same in Eqs. (30), (32) and (36), (38). Substituting M2 from (38) into Eq. (30) and taking account on values of angles results in  j  a2 þ b2 a2 þ b2 M1 ðjd1 Þ þ M1 d1 ¼ cq jd1 ð42Þ 2 4 4ab

j 4 jn, define an integer N rounding up the value Log2(j/jn). The value j0 ¼ j/2N lies in the interval 0r j r jn, and M(j0) can be

where in relationship q¼q(e) the expression for e through jd1 according to (33) is used. Analogously the constitutive equation for moment M3 is constructed (M2 from (32) is substituted into Eq. (36)):  j  c2 þ b2 c2 þ b2 aq M3 ðjd3 Þ þ M3 d3 ¼ jd3 ð43Þ 2 4 4cb

where M(jj  1) is known from the previous step. Obviously jN ¼ j, and the desired value M(j) has been found. For negative values of j, the antisymmetric continuation should be used. Let the function Q(j) have a limit value Ql when j increases. From Eq. (45) it is seen that the limit value of M(j) corresponding to increasing values of j equals Ql/2 and the initial derivative of the function M(j) equals

The corresponding result for the moment M2 can be deduced from an additional consideration of the deviatoric loading parallel to (x, z) plane:  j  c2 þ a2 c2 þ a2 M2 ðjd2 Þ þ M2 d2 ¼ bq jd2 ð44Þ 2 4 4ca Whereas in the case of the plane finite element, the moment in the spring connecting diagonal rods is expressed directly using the relationship q¼q(e) (see (7)), in the 3D-case, the corresponding equations become more difficult requiring the solution of an equation of the type (42)–(44). A suitable algorithm can be constructed as follows. Consider the general equation with initial values of the functions equal to zero: MðjÞ þMðj=2Þ ¼ Q ðjÞ

ð45Þ

For a sufficient small interval 0 r j r jn a suitable representation of the function M(j) can be taken in the form: MðjÞ ¼ S1 j þ S2 j2

ð46Þ

The requirement that the function on left-hand side in (45) equals Q(jn/2) and Q(jn) at j ¼ jn/2 and j ¼ jn, respectively, results in S1 ¼

8Q ðj=2Þ2Q ðjÞ , 3j

S2 ¼

16Q ðj=2Þ þ 8Q ðjÞ 5j2

ð47Þ

The value of jn is defined from the requirement that the maximum error in Eq. (45) after substitution of the function (46) is less than a chosen value. Let for example the function Q(j) be as follows: Q ð jÞ ¼

j 0:2 þ 9j9

ð48Þ

Taking jn ¼0.02 we have the error in Eq. (45) less than 0.00004 at the considered interval. In order to find a value of the function M(j) for an arbitrary magnitude of the argument

found directly by (46). Using Eq. (45) we perform N times the calculations:

jj ¼ 2jj1 , Mðjj Þ ¼ Q ðjj ÞMðjj1 Þ ðj ¼ 1,. . .,NÞ

M0 ð0Þ ¼ 2Q 0 ð0Þ=3

ð49Þ

ð50Þ

Constructing a function of the same type as Q(j) with this initial derivative and with limit value Ml ¼Ql/2 we obtain an approximate function for M(j): Ma ðjÞ ¼

1 Q ð4j=3Þ 2

ð51Þ

For the example (48), this approximate function will be Ma ðjÞ ¼

0:5j 0:15 þ 9j9

ð52Þ

In Fig. 7 this function is shown (dashed line) along with exact function defined according to the above described algorithm. Analogously approximate functions for M1(j), M2(j), and M3(j) are obtained using the known behavior of functions in right-hand side of Eqs. (42)–(44). The suggested approximation results in slightly decreased stiffness. Using the approximate treatment the subsequent correction of the function q(e) according to Eqs. (42)–(44) is required. Pursuant to (45) the function Q(j) after the correction will be Qa ðjÞ ¼ Ma ðjÞ þMa ðj=2Þ

ð53Þ

Since for an isotropic material the behaviors of q(e) and t(j) should be similar, for the sake of compatibility, it is advisable to make the corresponding correction also for the function t(j) (see (29)):

ta ðjÞ ¼ Ta ðjÞ þ Ta ðj=2Þ

ð54Þ

where (see (51)) Ta ðjÞ ¼

1 tð4j=3Þ 2

ð55Þ

790

G.B. Muravskii / Finite Elements in Analysis and Design 47 (2011) 784–795

From Fig. 7 it is clear that deviations of the corrected functions q(e) and t(j) from their initially assessed behaviors will be acceptable in engineering applications. Consider pressure (tension) p in the filling, which reacts on volume change only. The relative volume change D can be written in the form (for small deformations): 1 D ¼ ðeA1B1 þ eD1C1 þ eA4B4 þ eD4C4 þ eA1D1 þ eB1C1 4 þ eA4D4 þ eB4C4 þ eA1A4 þ eB1B4 þ eC1C4 þ eD1D4 Þ ð56Þ The notation of apexes of the 3D element corresponds to the explanation given before Eq. (1). Whereas in the plane case p was meant as average normal stress in the (x, y) plane, (sx þ sy)/2, now p is the whole average normal stress. In the case of the linear elasticity: p¼

2Gð1 þ nÞ D 3ð12nÞ

z O

x M 2m

y

2m 1m Fig. 8. Prism under action of self-weight; at upper face, uniform distributed tension equilibrates the weight.

ð57Þ The equation for moments in spiral springs at the apexes of the finite element has the form corresponding to Eq. (29):

2.2.1. Determination of stresses As in the plane case (see (8) and (9)) stresses relating to boundary faces of the 3D element are determined using averaging node forces generating by spiral springs over the corresponding areas. Remember that moments acting in planes parallel to (x, y) plane, (x, z) plane, and (y, z) plane are designated with subscripts 1 (4), 2 (5), and 3 (6), respectively (see Fig. 1b). Eq. (8) is generalized as follows:

txy ¼

1 ðmB1 þ mD1 mA1 mC1 þ mB4 þ mD4 mA4 mC4 Þ abc

ð58Þ

tyz ¼

1 ðmA3 mD3 þ mC3 mB3 þ mA6 mD6 þmC6 mB6 Þ abc

ð59Þ

txz ¼

1 ðmA2 þ mC2 mB2 mD2 þ mA5 þmC5 mB5 mD5 Þ abc

ð60Þ

where indices determine both the point and plane. The normal stresses are determined by the pressure (tension) in the filling and the moments in spiral springs, which connect diagonal rods:

sx ¼ p

2ðM1 þM4 Þ 2ðM2 þ M5 Þ  2 ða2 þ b2 Þc ða þ c2 Þb

ð61Þ

sy ¼ p þ

2ðM1 þ M4 Þ 2ðM3 þ M6 Þ þ ða2 þ b2 Þc ðb2 þ c2 Þa

ð62Þ

sz ¼ p

2ðM3 þ M6 Þ 2ðM2 þ M5 Þ þ ðc2 þb2 Þa ða2 þ c2 Þb

ð63Þ

The stresses given in Eqs. (58)–(63) should be related to the center of the element. 2.2.2. Linearly elastic isotropic body Consider an isotropic linearly elastic material characterized with shear modulus G and Poisson’s ratio n. The relationship between q and e results in q(e)¼ t(2e)¼2Ge and we obtain from Eqs. (42)–(44): 3 ða2 þ b2 Þ2 M1 ¼ G cjd1 2 8ab

ð64Þ

3 ða2 þ c2 Þ2 M2 ¼ G bjd2 2 8ac

ð65Þ

3 ðc2 þ b2 Þ2 M3 ¼ G ajd3 2 8cb

ð66Þ

m¼

Gabc j 8

ð67Þ

In addition to Eqs. (64)–(67), Eqs. (57) and (56) should be applied in the considered linearly elastic case. The obtained equations can be used for determination of linear relationships connecting the apexes displacements and forces acting upon the apexes from rods and the filling. The corresponding matrix has the order 24 and can be constructed without difficulties. However, determination of the stiffness matrix can be avoided. When solving dynamic problems or static problems with the dynamic treatment we apply iterations on each time step, and the corresponding forces can be calculated directly through displacements of nodes of the element. The stiffness matrix is needless in such a treatment. For testing 3D elements in a relatively simple static problem, consider the following example. A prism (Fig. 8) is under selfweight action equilibrated by the tension uniformly distributed over the upper face of the prism; the center point of the upper face is fixed. This problem has a simple solution given in [12], which results in the following expressions for horizontal (along axis x) and vertical displacements at the upper face: u¼

nrgxL , 2ð1 þ nÞG

v¼

nrgðx2 þ z2 Þ 4ð1 þ nÞG

ð68Þ

where L is size of the prism in vertical direction. In calculations the quarter of the prism corresponding to positive coordinates x and z is considered. Using elements with a¼0.125 m, b¼0.5 m, c¼0.25 m we have 4  4  4 ¼64 elements for the quarter. For G¼ 100,000 Pa, n ¼1/3, r ¼5000 kg/m3, g¼9.8 m/s2 displacements at the point M(0.5, 0, 0) are u¼  0.06110, v¼0.007540 whereas the corresponding exact values are u¼  0.06125, v¼0.007656. Taking a¼0.0625 m, b¼0.25 m, c ¼0.125 m, i.e. 8  8  8¼512 elements, the values u ¼  0.06120, v¼0.007617 are obtained.

2.2.3. Non-linear material with limit stresses obeying Coulomb– Mohr law In contrast to the plane case, the conditions of failure of the type of (19) and (20) should be written now for the three planes: (x, y), (y, z), and (x, z) using mean stresses (sx þ sy)/2, (sy þ sz)/2, (sx þ sz)/2 instead of p in Eq. (19) for determining limit deviatoric stresses for these planes Rxy, Ryz, Rxz. For determining limit values of the shear stresses txy, tyz, txz one should use deviatoric stresses qxy, qyz, qxz instead of q in Eq. (20). The following 6 equations

G.B. Muravskii / Finite Elements in Analysis and Design 47 (2011) 784–795

determine the failure conditions for the 3D finite element: Rxy ¼ c cosðfÞ þ Ryz ¼ c cosðfÞ þ Rxz ¼ c cosðfÞ þ

sx þ sy 2

sy þ sz 2

sx þ sz 2

sinðfÞ

ð69Þ

sinðfÞ

ð70Þ

sinðfÞ

ð71Þ

txymax ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2xy q2xy

ð72Þ

tyzmax ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2yz q2yz

ð73Þ

txzmax ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2xz q2xz

ð74Þ

791

Using Eqs. (78)–(83) instead of the exact hyperbolic relations for shear and deviatoric stresses we consider functions close to hyperbolic functions requiring limiting values and initial derivatives. This treatment allows us to diminish computational time, which increases when solving exact Eqs. (42)–(44). The closeness of results for the approximate and exact moment determination is shown in an example below. For the pressure p in the filling, Eq. (57) can be used with n0 and G0 instead of n and G, respectively. Initial shear modulus G0 can be determined in the case of soils using Eq. (23) with values of a and b given after Eq. (24).

3. Examples

where  sx þ sy   qxy ¼ sy   2  sy þ sz   qyz ¼ sy   2  sx þ sz   qxz ¼ sx   2

ð75Þ ð76Þ ð77Þ

The considerations given after Eq. (20) relating to the case q4R hold also here for values Rxy, Ryz, Rxz and qxy, qyz, qxz in which corresponding changes are made if necessary. Assuming, for example, the hyperbolic law with the given limiting values, we can write 6 equations describing the behavior of shear and deviatoric stresses and afterwards to construct the constitutive equations for the moments m1(m4), m2(m5), m3(m6), M1(M4), M2(M5), and M3(M6) using Eq. (29) (three times) and Eqs. (42)–(44). Let us apply the above-mentioned approximate treatment. In the case of the hyperbolic law, according to (51), (42)–(44) the following equations for moments M1, M2 and M3 hold: M1 ðjd1 Þ ¼

a2 þ b2 l1 jd1 , cRxy 8 Rxy =ð2G0 Þ þ l1 9jd1 9

l1 ¼

a2 þ b2 3ab

ð78Þ

M2 ðjd2 Þ ¼

a2 þ c2 l2 jd2 bRxz , 8 Rxz =ð2G0 Þ þ l2 9jd2 9

l2 ¼

a2 þ c2 3ac

ð79Þ

M3 ðjd3 Þ ¼

b2 þc2 l3 jd3 , aRyz 8 Ryz =ð2G0 Þ þ l3 9jd3 9

l3 ¼

b2 þ c2 3bc

ð80Þ

The exact equations for m1, m2, m3 can be obtained using (29) and hyperbolic laws for txy, tyz, txz with the corresponding limit values (72)–(74). As mentioned above, when applying approximate equations for moments M1, M2, and M3 (i.e. correcting the behavior of deviatoric stresses), it is appropriate to make the corresponding changes in shear stresses. Using Eqs. (54) and (55) results in ! abc 1 1 m1 ðjÞ ¼ þ txymax j 16 3txymax =ð4G0 Þ þ 9j9 3txymax =ð2G0 Þ þ 9j9

3.1. 3D problem for horizontal layer under action of normal load uniformly distributed over square area on layer surface For simplification of the problem, we will neglect the own weight of the layer and, additionally, assume that the shear modulus G0 ¼10 MPa is constant independent of stresses in the layer; the pore pressure is not taken into account. Let n0 ¼0.3, c¼10,000 Pa, f ¼201, the value of given vertical pressure p0 ¼120 kPa, thickness of the layer H¼11 m, and size of loaded square at the surface B¼H (Fig. 9). By virtue of symmetry, only the quarter of the system can be considered (x Z0, z Z0, 0ry r11 m) with the corresponding conditions imposed on displacements at the planes x¼0 and z ¼0. Instead of infinite dimensions in x and z directions we take finite lengths 66 m. The bottom face of the layer is fixed, at the vertical faces x¼33 m the horizontal displacements ux ¼ 0 and at the vertical faces z ¼33 m the horizontal displacements uz ¼0. We use cubic finite elements with a¼b ¼c¼1 m (11  33  33 ¼11,979 elements for the quarter), the border of the loaded area at the layer surface lies on the lines x¼5.5 m and z¼5.5 m at the surface, i.e. between element borders. Such a choice of dimensions of elements and loaded area allows diminishing the influence of discontinuity in the applied loading. Appearing static problems are solved dynamically using slow growth of the applied load from 0 to p0 according to: ( 1cosðpt=t Þ 0 p0 ðt rt0 Þ 2 ð84Þ pðtÞ ¼ ðt 4t0 Þ p0 where t is time, t0 ¼time of load growth (t0 ¼5 s in calculations). For chosen values a, b, c, G0, n0, r ¼2000 kg/m3 (density), the time step Dt ¼0.001 s ensures the required accuracy. The mass of an element is transferred to element’s nodes. For excluding concomitant vibrations, external dissipative forces proportional to velocities (with the corresponding coefficient 500 N s/m) are applied in nodes of elements. On a current time step, knowing initial displacements, velocities and accelerations of nodes of the

p0

z

x B

ð81Þ m2 ðjÞ ¼

! abc 1 1 þ txzmax j 16 3txzmax =ð4G0 Þ þ 9j9 3txzmax =ð2G0 Þ þ 9j9

y H

ð82Þ abc 1 1 þ tyzmax j m3 ðjÞ ¼ 16 3tyzmax =ð4G0 Þ þ9j9 3tyzmax =ð2G0 Þ þ9j9

!

ð83Þ

Fig. 9. Layer loaded by surface uniform pressure applied to square area.

792

G.B. Muravskii / Finite Elements in Analysis and Design 47 (2011) 784–795

0.3

Table 2 Vertical displacements of center of loaded area for different values of strength diminishing coefficient fi after 4 and 5 iterations.

0.2

After 4 iterations After 5 iterations

stress/p0

σy σx τxy

fi ¼1.3

1.35

1.4

1.45

0.20277 0.20383

0.25507 0.26348

0.31609 0.33518

0.50796 0.50992

8 0.1

6 0 0.2

0.4

0.6

0.8

1

1/Uy

0

depth/H

4

Fig. 10. Normalized stresses under corner of loaded area; symbols correspond to Ref. [14].

2

system, the corresponding values at the step end are determined by (internal) iterations using the method of mean accelerations. Firstly, the end accelerations are taken equal to the accelerations at the start step, and the end displacements and velocities are found using the accelerations. This allows us to determine the forces with which elements act on their nodes and from Newton’s equations of motion to find new values for the end accelerations. After two or three repetitions, the end accelerations, velocities and displacements are determined with acceptable accuracy, and one can deal with the next time step. Numeration of elements is performed using indices i, j, k which increase in x, y, z directions from 1 to 33, 11, 33, respectively, analogously the nodes are numerated with changing lowest index to 0. This numeration is used when determining forces acting on elements nodes (reactions of elements) and nodes accelerations at the end of each time step in the process of internal iterations, which realize the method of constant average accelerations. Firstly, the linearly elastic case is addressed, for which the external iterations usually associated with initially unknown material properties (depending on unknown stresses) are not needed. In Fig. 10 normalized stresses under the corner of the loaded area are presented for n0 ¼0.3 being related to the centers of the corresponding cubic elements. Symbols in the figure are for results from [13,14]. The authors of [14] accepted, apparently, zero magnitude for surface values of txy, however, because of discontinuity in the applied loading, the point limit normalized value for this stress equals 0.159155 (for the given conditions), which can be found from considering this problem for the halfspace instead of the layer. Paradoxically, the exact values of shear stresses under the corner of loaded area increase while approaching the surface. Results corresponding to the elastic solution are used further for constructing solutions for non-linear (hyperbolic) material with Coulomb–Mohr failure conditions, namely the stresses in the elastic solutions serve as an initial approach for above described iterations. Numerous calculations show that instead of simple iterations, in which an end result ri after i iterations is taken as an initial result ri þ 1,0 for the (i þ1)th iteration, the following scheme is more successful: ri þ 1,0 ¼ lri,0 þ ð1lÞri

ð85Þ

with l ¼0.3 (for i¼0 we take ri,0 ¼ri). In the case of exact determination of moments in springs connecting the diagonal rods, the value of jn entering Eq. (47) is taken equal to 0.000001.

0 1.1

1.2

1.3

1.4

1.5

1.6

f Fig. 11. FOS determination using inverse displacements.

Below we demonstrate how the new finite elements allow solving problems associated with material failure. Introduce the factor of safety FOS in the usual manner, i.e. as the coefficient of strength parameters decrease leading to system failure. Consider a series of values of c and f: cj ¼

  c tanðfÞ , fj ¼ arctan f f

ð86Þ

where f¼fj ¼1.1, 1.2, 1.3, 1.35, 1.4, 1.45. For these values of strength parameters, the corresponding solution is obtained by the method of iterations (5 iterations are used) starting from the elastic solutions. In Table 2, vertical displacements Uy of the center of the loaded area are given for fourth and fifth iterations correspondingly to the taken values of fj. Increase in the displacements and worsening of convergence indicate the nearing to the failure. In Fig. 11 the values inverse to displacements (taken after 5 iterations) are given versus f ¼fj (black squares); also the results (black circles) are shown, which correspond to the above described approximate determination of moments in spiral springs (Eqs. (78)–(83)). The results are close to a straight line whose continuation until f-axis delivers an assessed value of FOS¼1.549. As an alternative to this method of FOS determination, which appears to be very simple here due to the closeness of the corresponding plot to the straight line (which not always takes place), consider the following methodology. For each performed solution corresponding to fj, we fix normal stresses (i.e. strength properties of the material are prescribed and the iterations are abandoned), and determine for which value of f¼fn the failure takes place. Doing so we consider the vertical velocity V of the point in the center of the loaded area and determine by repeated calculations the value of f, which leads to unbounded growth of the velocity; in these calculations the forces of resistance proportional to velocity mentioned after Eq. (84) are not introduced into the dynamic treatment of the static problems. For clarification, address Fig. 12 where the stresses obtained for

G.B. Muravskii / Finite Elements in Analysis and Design 47 (2011) 784–795

793

30m

1.6 A

f = 1.49 10m

20m

1.2

f = f* = 1.484

V

10m

f = 1.48

0.8

60m

0.4

0

Fig. 14. Scheme of embankment.

0

4

8

12

Table 3 Displacements of apex (point A) of embankment different values of strength diminishing coefficient fi after 4 and 5 iterations.

t (sec) Fig. 12. On FOS determination using velocity behavior when fixing initial stress distribution for fj ¼ 1.35.

1.6

1.5

f*

1.4

1.3

f* = f

1.2

1.1 1.1

1.2

1.3

1.4

1.5

1.6

f Fig. 13. FOS determination using fixed initial stresses taken for different f ¼fj and evaluating failure values f n.

fj ¼1.35 after 5 iterations are used for determination of strength parameters. For 3 values of f the velocity V as a function of time is presented; the desired value of f n equals 1.484. The value of FOS can be found considering functions f n ¼f n(f), where f means the coefficient of strength decrease for which stress distribution is taken when determining f n, and linear function f n ¼f (Fig. 13). The closer is the taken coefficient fj to the FOS, the closer are f n and fj to each other. Therefore, the point of intersection of the corresponding plots determines the FOS. Points f n obtained for f ¼1.3, 1.35, 1.4, 1.45 allow a good approximation with straight line which leads to FOS¼1.548 which is very close to the result found from the displacements consideration. Exact and approximate methods of determination of the moments in the springs lead practically to the same results relating to the FOS (Fig. 11). 3.2. Application of new finite elements to slope stability investigation For this important problem of geomechanics, various different methods have been developed based on the concept of limit equilibrium or on finite element technique. As the objective of

After 4 iterations After 5 iterations

fi ¼ 0.86

0.87

0.88

0.89

0.025696 0.026055

0.026792 0.027392

0.028308 0.029323

0.030925 0.033298

this paper is only to illustrate how new finite elements can be successful for solving some problems in geomechanics, we do not present an analysis of numerous works in the area of slope stability confining ourselves to few references ([15–19]). Consider an embankment of height H¼10 m, sloping at angle 451 with friction angle f ¼201, unit weight w0 ¼20 kN/m3, cohesion c¼12.38 kPa. According to the solution of Chen [19] the value of FOS equals 1.0 without taking into account the pore pressure. This example has been studied also in [18] where a suggestion is made to adjust Poisson’s ratio for avoiding a spurious failure, which is connected with assumption that material is elastic–perfectly plastic obeying Mohr–Coulomb’s yield criterion. However, doing so we actually consider another material instead of the one given in the problem. Therefore it is desirable to avoid the assumption of elastic–perfect plasticity. In calculations we consider the system shown in Fig. 14; the bottom plane is fixed and at vertical boundaries, horizontal displacements equal zero. We apply square elements with a¼b ¼0.5 m; the inclined line becomes stepped. The same methodology as in the previous subsection is applied here for assessment of the factor of safety. The initial shear modulus is taken in the form of (23) with a ¼10 kPa and b ¼ 200 kPa. Going from the elastic solution we construct solutions for the following values of f¼fj ¼0.86, 0.865, 0.87, 0.875, 0.88, 0.885, 0.89, 0.895 in Eq. (86). Note that the elastic solution itself requires the iterations since the initial shear modulus depends on unknown stresses; at this stage the stresses before iterations are taken according to Eq. (25) using depths of current points and n ¼ n0 ¼0.35 (initial Poisson’s ratio). In dynamic calculations we apply Eq. (84) with t0 ¼5 s changing pressure p to unit weight w; time step Dt ¼0.0005 s; scheme (85) is used for iterations. In Table 3, full displacements U of the point A (Fig. 14) are given for fourth and fifth iterations for the taken values of fj. Unproportional growth of displacements and differences in results for fourth and fifth iterations (worsening of convergence) with increase of fj take place. It is interesting to compare stresses at the beginning and the end of the fifth iteration, which could estimate accuracy of results. Such a comparison is made in Fig. 15 (for fj ¼0.86) and Fig. 16 (for fj ¼0.89) for stresses in elements lying in several horizontal rows (counting from bottom). For fj ¼0.89 the noticeable difference is seen for k¼21 (the first row in the part with slope) for values of i for which the elements are close to failure. In Fig. 17 the values inverse to displacements (taken for fifth iteration) are given

794

G.B. Muravskii / Finite Elements in Analysis and Design 47 (2011) 784–795

300000

200000 k = 15

fj = 0.86

160000

k = 15

250000

21

200000

21 σx

σy

120000 σx

80000

σy

28

σx0

28 σy0

100000

35

40000

0

150000

35

0

40

80

50000

0

120

0

40

80

i

120

i

Fig. 15. Stresses sx (a) and sy (b) in elements for fj ¼ 0.86 after 5 iterations; i is number in row from left, k is number of row from bottom; dotted lines are for initial values.

200000

300000 k = 15

fj = 0.89

160000

k = 15

250000

21

200000

21 σx

σy

120000 150000 28 80000

28 100000 35

40000

35

50000

0 0

40

i

80

120

0

0

40

i

80

120

Fig. 16. Stresses sx (a) and sy (b) in elements for fj ¼ 0.89 after 5 iterations; i is number in row from left, k is number of row from bottom; dotted lines are for initial values.

versus f¼fj (polygonal line). An assessment for FOS can be made continuing until f-axis the straight sections of the polygonal line; this treatment should deliver overestimated values of FOS (judging from the character of the obtained line). The following interval for FOS can be assessed:

show that the value of FOS is very close to 0.9. The value FOS ¼1 obtained in [19] and based on the method of limit equilibrium with assumed kinematics of failure is overestimated.

0:895 o FOS o0:922

4. Conclusions

ð87Þ

The alternative method presented in the previous subsection leads in the considered case to more narrow bounds. Taking stress distribution (after 5 iterations) for each fj we determine the corresponding coefficient f¼f n leading to failure. One should find the point of intersection of polygonal line f n ¼f n(f) (Fig. 18) and the straight line f n ¼ f . Assessed overestimated values of FOS found for different sections of the polygonal line according to Fig. 18 are as follows: 0.898, 0.904, 0.907, 0.912, 0.915, 0.918, and 0.924. Accuracy of these results are higher, the less are fj, which determine the corresponding section of the polygonal line. Note that the section between points for fj ¼0.86 and fj ¼0.865 leads to an upper bound for FOS (0.924) sufficiently close to more conservative (however less accurate) value (0.898) corresponding to the section with fj ¼0.89 and fj ¼0.895. All these considerations

An important characteristic of the suggested finite elements is the independent behavior of their components, which allow us to obtain simple constitutive equations for the components and to study effectively the behavior of materials with complicated properties, particularly, a non-linear material obeying the Mohr–Coulomb failure condition. Interestingly, the condition leads to the separate strength values relating to the two kinds of springs entering the finite elements. An additional point to emphasize is the mechanical simplicity of the elements. As a rule, the considered problems are characterized by impossibility to write from the outset constitutive equations in an explicit form, since the limit stresses depend on the unknown actual stress distribution. For avoiding this difficulty, the method is applied, which consists in the iterative improvement of initially prescribed

G.B. Muravskii / Finite Elements in Analysis and Design 47 (2011) 784–795

40

1/U

30

20

10

0 0.84

0.88

0.92 f

0.96

1

Fig. 17. FOS determination based on consideration of displacements after 5 iterations.

that application of traditional finite elements seems to be unable to overcome the 3D problem in which failure conditions, defining limiting stresses, go along with a non-linear behavior such as obeying the hyperbolic law. Using the methodology described above we have to determine the pressure in the filling and moments in the spiral springs through the displacements of element’s nodes and afterwards to find forces acting on the nodes (at each step of the gradual loading); all are realized straightforwardly. Trying to solve the problem with the help of conventional continuum elements we have the possibility to find directly (knowing nod displacements) in an internal point of an element all the shear stresses and the deviatoric stresses for coordinate planes xy, xz, yz: (sx  sy)/2, (sx– sz)/2, (sy–sz)/2 using, say, hyperbolic law for these components. However, the deviatoric stresses found by this means usually prove to be incompatible (the second of them does not equal the sum of the first and third). Thus we have no possibility to determine the node forces. Apparently the new suggested elements have the properties which are not inherent with the conventional finite elements, whose applications for solution of non-linear problems where limiting stress conditions define the non-linear behavior are absent in publications. References

1

0.96

f*

f* = f*( f )

f* = f

0.92

0.88 0.86

795

0.88

0.9

f

0.92

0.94

Fig. 18. FOS determination using fixed initial stresses taken for different f ¼fj and evaluating failure values f n.

stresses, until stresses at the start and end of an iteration are close enough to each others. The convergence of the method is worsening with nearing to the state of failure, for which estimation an extrapolation is performed going from the results calculated for parameters still ensuring the convergence. Testing the finite elements in linearly elastic problems shows a good accuracy, especially in the case of cubic (square) or close to cubic (square) form of the elements. 3D problems require (in the nonlinear case) a specific algorithm for determination of moments in springs connecting diagonal elements. An alternative simple approximate method seems to be enough accurate for engineering applications. On two examples important for geomechanics the possibility is demonstrated for an adequate assessment of failure conditions of the considered systems. It should be noted

[1] A.E.H. Love, A Treatise on the Mathematical Theory of Elasticity, Dover, New York, 1944. [2] A.R. Rzhanitsin, Structural Mechanics, Vysshaia Shkola, Moscow, 1982 (in Russian). [3] V.L. Streeter, E.B. Wylie, Two- and three-dimensional transients, Journal of Basic Engineering. Transactions of ASME 90D (1968) 501–510. [4] C.N. Papadakis, Soil transients by characteristics method. Dissertation, University of Michigan, 1973. [5] G.B. Muravskii, 2D and 3D seismic response analysis using experimental data on cyclic behaviour of soils, in: Proceedings of the 4th International Conference on Cyclic Behaviour of Soils and Liquefaction Phenomena, 2004, pp. 572–600. [6] J.M. Duncan, C.Y. Chang, Nonlinear analysis of stress and strain in soils, Journal of the Soil Mechanics and Foundations Division, ASCE 96 (1970) 1629–1653. [7] S.P. Timoshenko, J. Gudier, Theory of Elasticity, McGraw-Hill, New York, 1970. [8] G. Muravskii, On description of hysteretic behaviour of materials, International Journal of Solids and Structures 42 (2005) 2625–2644. [9] G. Muravskii, Application of hysteresis functions in vibration problems, Journal of Sound and Vibration 319 (2009) 476–490. [10] B.O. Hardin, V.P. Drnevich, Shear modulus and damping in soils: measurement and parameter effects, Journal of the Soil Mechanics and Foundations Division, ASCE 98 (1972) 603–624. [11] B.O. Hardin, V.P. Drnevich, Shear modulus and damping in soils: design equations and curves, Journal of the Soil Mechanics and Foundations Division, ASCE 98 (1972) 667–692. [12] S.P. Timoshenko, Theory of Elasticity, McGraw-Hill, New York, 1970. [13] H.G. Poulos, E.H. Davis, Elastic solutions for soil and rock mechanics, John Wiley & Sons, New York, 1974. [14] D.M. Milovic, J.P. Tournier, Stresses and displacements due to rectangular load on a layer of finite thickness, Soil and Foundations 11 (1971) 1–27. [15] J.M. Duncan, State of the art: limit equilibrium and finite-element analysis of slopes, Journal of Geotechnical Engineering 122 (1996) 577–596. [16] R. Baker, Determination of the critical slip surface in slope stability computations, International Journal for Numerical and Analytical Methods in Geomechanics 4 (1980) 333–359. [17] M.M. Farias, D.J. Naylor, Safety analysis using finite elements, Computers and Geotechnics 22 (1998) 165–181. [18] H. Zheng, D.F. Liu, C.G. Li, Slope stability analysis based on elasto-plastic finite element method, International Journal for Numerical Methods in Engineering 64 (2005) 1871–1888. [19] W.F. Chen, Limit Analysis and Soil Plasticity, Elsevier, Amsterdam, 1975.