Rock Mech. Rock Engng. (2008) 41 (4), 601–628 DOI 10.1007/s00603-006-0120-9 Printed in The Netherlands
Formulation of a Three-dimensional Numerical Manifold Method with Tetrahedron and Hexahedron Elements By
Y. M. Cheng and Y. H. Zhang Department of Civil and Structural Engineering, Hong Kong Polytechnic University, Hong Kong, P.R. China Received December 22, 2005; accepted July 31, 2006 Published online March 13, 2007 # Springer-Verlag 2007
Summary Numerical Manifold Method (NMM) is an extension of Discontinuous Deformation Analysis which contains and combines FEM and joint=block oriented DDA in a unified form. At present, application of NMM is mainly limited to two-dimensional analysis which is not realistic for many underground excavation works. In the present paper, three-dimensional NMM based on tetrahedron and hexahedron elements is proposed. The corresponding three-dimensional cover contact detection algorithms which have some fundamental properties different from the corresponding two-dimensional analysis are also developed. Based on the present formulations, engineering problems with arbitrarily oriented discontinuities can be modeled and realistic modeling of construction works becomes possible. Keywords: Manifold method, discontinuous deformation analysis, joints, contacts.
1. Introduction DDA is originated by Shi (1985) and is a powerful numerical technique for discontinuous deformation analysis of blocky systems. With the developments of 2D DDA by Shi (1993), Cheng (1998, 2000, 2002) and others, DDA is now a powerful tool for problems in discontinuous media and has already been used in the solution of some geotechical problems. Numerical Manifold Method (NMM) is proposed by Shi (1991) as an extension of DDA. It is a very flexible numerical method which contains and combines finite element method (FEM) and blocky systems considered in DDA in a unified form. This method combines the advantages of DDA and FEM together and is applicable to a much wider range of problems.
602
Y. M. Cheng and Y. H. Zhang
In NMM, continuous and discontinuous cover functions are used to solve problems in continuous and discontinuous domains. This method can also be reduced to the classical finite element method or DDA but can be used for more general conditions that are found in many practical geotechnical problems. At present, nearly all NMM cover functions are adopted directly from finite element shape functions due to the mature development of various techniques in FEM. Simple cover functions are usually adopted for NMM to avoid complicated formulation and the use of higher order cover functions are not commonly used. The authors (2002) have however shown that the use of the more advanced twodimensional Wilson non-conforming elements to NMM can be useful in obtaining more accurate results for simple covers. The choice of cover function is hence also important for a proper NMM analysis, particularly when the size of each block is significant. Practical geotechnical problems are generally three-dimensional in nature. The geometry of the construction works, the properties of the various rock materials and the discontinuities all possess strong three-dimensional characteristics. At present, 3DDEM has been developed by Cundall (1985) in a DEM code 3DEC while there are some simple applications of 3D-DDA by Wang (1996), Jiao (1998), Jiang (2000) and the authors (2002). NMM is currently mainly limited to two-dimensional analysis and there is little research work in the 3D development although it is very important for engineering problems. The difficulties for three-dimensional NMM are: 1. definition of cover functions; 2. three-dimensional contact detection; 3. three-dimensional contact stress-strain relation. In this paper, three-dimensional formulations of NMM based on tetrahedron and hexahedron elements are proposed. For 2D NMM, the choice of the cover weight function is practically arbitrary as 2D geometry calculation is relatively simple. For 3D NMM, a difficult problem is the suitable choice of cover weight function. For ease of programming and speed of computation, volume coordinates are suggested as the cover weight functions for 3D tetrahedron element, while numerical integration similar to that used in finite element analysis is proposed for hexahedron element. All the sub-matrices for the equilibrium equations of 3D NMM are derived based on energy consideration. For 2-D contacts, the directions of shear and the friction forces are easily defined along the contact edge. For the corresponding 3-D contacts, the directions of the contact forces=stresses cannot be determined in a simple way at the start of a solution step. To deal with this problem, two iterative methods for the detection of contact directions are developed for tetrahedron and hexahedron elements. The present formulations and the developed computer codes are verified by two simple examples and applications to more complicated coal extraction problems are also discussed. 2. Tetrahedron and Hexahedron Elements Tetrahedron and hexahedron elements as shown in Fig. 1a and b are the common three-dimensional elements used in finite element analysis. These elements possess the advantage of simplicity in formulation and convenience in mesh generation (possibly
Formulation of a Three-dimensional Numerical Manifold Method
603
Fig. 1. a) Tetrahedron element, b) mapping from ð; ; Þ local coordinates to ðx; y; zÞ global coordinates
the most important reason for its popularity in application). Combining the use of these two elements, many three-dimensional problems can be modeled with ease. For tetrahedron elements, volume coordinates are commonly used in finite element analysis and also adopted in three-dimensional NMM. The volume coordinates Ni of an arbitrary point Pðx; y; zÞ in the element as shown in Fig. 1a are given by: N1 ¼
volðP234Þ ; volð1234Þ
N2 ¼
volðP341Þ ; volð1234Þ
N3 ¼
volðP412Þ ; volð1234Þ
N4 ¼
volðP123Þ ð1aÞ volð1234Þ
For tetrahedron elements, the ‘‘stiffness’’ matrices are formulated directly in terms of the vertices coordinates, in a way which is similar to the original DDA or NMM formulation by Shi (1991). For hexahedron elements, the matrices are formulated using numerical integration instead of the coordinates of vertices as this will help to simplify the lengthy algebra. For hexahedron, the C8 shape functions of 3D brick elements from finite element method are adopted as the cover weight functions. The shape functions of any point Pð; ; Þ in a C8 isoparametric element are 1 Ni ¼ ð1 þ i Þð1 þ i Þð1 þ i Þ 8 Eqs. (1a) and (1b) fulfill the requirement of X Ni ¼ 1
ð1bÞ
ð2Þ
It is easy to prove that Ni ¼ 1; 0 4 Ni 4 1;
when P ¼ i
when P 2 volð1; 2; 3; 4Þ
or
volð1; 2; 3; 4; 5; 6; 7; 8Þ
ð3Þ
604
Y. M. Cheng and Y. H. Zhang
The cover weight function is usually set to be equal to the interpolation function Ni . The displacements of any point Pðx; y; zÞ in a tetrahedron or hexahedron element are given by: uðx; y; zÞ ¼
4;8 X
Ni ðx; y; zÞui
i¼1
ðx; y; zÞ ¼
4;8 X
Ni ðx; y; zÞi
ð4Þ
i¼1
wðx; y; zÞ ¼
4;8 X
Ni ðx; y; zÞwi
i¼1
Equation (4) can be rewritten as Eqs. (5a) and (5b) for tetrahedron and hexahedron elements as 9 8 8 9 Deð1Þ > > > =
eð2Þ ð5aÞ ¼ fTgfDg ¼ ð½Teð1Þ ½Teð2Þ ½Teð3Þ ½Teð4Þ Þ Deð3Þ > : ; > > > w ; : Deð4Þ 9 8 ½Deð1Þ > > > > > > > ½Deð2Þ > > > > > > > > > > ½Deð3Þ > 8 9 > > > > > > > > > = < ½Deð4Þ = ¼ fTgfDg ¼ ð½Teð1Þ ½Teð2Þ ½Teð3Þ ½Teð4Þ ½Teð5Þ ½Teð6Þ ½Teð7Þ ½Teð8Þ Þ > > > ; : > > ½Deð5Þ > > > w > > > > > > ½D eð6Þ > > > > > > > > > > ½D > > eð7Þ > > ; : ½Deð8Þ ð5bÞ where,
2
NeðiÞ ½TeðiÞ ¼ 4 0 0
0 NeðiÞ 0
3 0 0 5; NeðiÞ
8 9 < ueðiÞ = fDeðiÞ g ¼ eðiÞ ; : ; weðiÞ
ð6Þ
For tetrahedron NeðiÞ ¼ fi1 þ fi2 x þ fi3 y þ fi4 z
ð7aÞ
1 NeðiÞ ¼ ð1 þ i Þð1 þ i Þð1 þ i Þ 8
ð7bÞ
For hexahedron
The coefficients fij for tetrahedron are expressed in terms of the coordinates of vertices explicitly and are given in the appendix. In Eq. (7b) for hexahedron ð; ; Þ, are the
Formulation of a Three-dimensional Numerical Manifold Method
605
local coordinates corresponding to the global coordinates ðx; y; zÞ which can be calculated using an iterative method as suggested in a later section of the present paper. 3. Equilibrium Equations for 3-Dimensional NMM Based on Tetrahedron/Hexahedron Elements The matrices of equilibrium equations for 3-dimensional NMM can be considered as an extension of the corresponding equations in 2D NMM analysis. From (4), we can obtain the strain in a tetrahedron or hexahedron element as: 9 8 @u > > > > > > > > > @x > > > > > > > @ > > 8 9 > > > > > > " x > > > > @y > > > > 9 8 > > > > > > > > > > > Deð1Þ > " y > > > > > @w > > > > > > > > > > = = > < = <" > eð2Þ @z z ¼ ½Be fDe g ¼ ð ½Beð1Þ ½Beð2Þ ½Beð3Þ ½Beð4Þ Þ ¼ f"e g ¼ @u @ > > > > > D > > > > > > > þ > > > xy > > eð3Þ > > > > ; > : > > > > @y @x > > > Deð4Þ yz > > > > > > > > > : ; > @ @w > > > > zx > > þ > > > > @z @y > > > > > > > > @w @u > > > > ; : þ > @x @z ð8aÞ 9 8 @u > > > > > > > > > @x > > > > > > > @ > > 8 9 > > > > > > " x > 9 8 > > > @y > > > > > > > > Deð1Þ > > > > > > > > > > > > " y > > > > > @w > > > > > > > > > > = = > = < <" > < Deð2Þ > @z z ¼ ½Be fDe g ¼ ð ½Beð1Þ ½Beð2Þ ½Beð8Þ Þ ¼ f"e g ¼ . @u @ > > > > .. > xy > > > > > > > > > > > > þ > > > > > > > > > > > > > @y @x > > ; > > : > > > > > yz > D > > eð8Þ > : ; > @ @w > > > > zx > > þ > > > @z @y > > > > > > > > > > > @w @u > ; : þ > @x @z ð8bÞ where
2
fi2 60 6 60 ½BeðiÞ ¼ 6 6 fi3 6 40 fi4
0 fi3 0 fi2 fi4 0
3 0 07 7 fi4 7 7 for tetrahedron 07 7 fi3 5 fi2
ð9aÞ
606
Y. M. Cheng and Y. H. Zhang
2
@Ni 6 @x 6 6 6 0 6 6 6 6 0 6 ½BeðiÞ ¼ 6 @N 6 i 6 6 @y 6 6 0 6 6 4 @N i @z
0 @Ni @y 0 @Ni @x @Ni @z 0
3 0 7 7 7 0 7 7 @Ni 7 7 7 @z 7 for hexahedron 7 0 7 7 7 @Ni 7 7 7 @y 7 @Ni 5 @x
ð9bÞ
Equation (9b) is obtained from (9c) according to an iterative inverse mapping method which will be discussed in the following. 9 8 9 8 @Ni > @Ni > > > > > > > > > > @ > > @x > > > > > > > > = < @N > < @N = i i ð9cÞ ¼ J 1 > > @y > @ > > > > > > > > > > > > @Ni > > @Ni > > > > ; > : ; : @z @ The stress-strain relation fg ¼ ½Ef"g can be expressed as 3 2 1 0 0 0 8 9 8 9 6 x > "x > 1 0 0 0 7 7> 6 > > > > 7> 6 > > > > > > > > > 7 > > 6 " 1 0 0 0 y> y > > > > > > > 7> 6 > > > > = < = < 7 6 1 2 " E z z 7 6 0 0 0 7 0 0 ¼ 6 2 > 7> ð1 þ Þð1 2Þ 6 xy > xy > > > > > > > > 7> 6 > > > > 1 2 > > > > 7 6 > > > > 0 0 0 0 0 > > > > yz yz 7 6 > > > 2 ; ; : > : 7 6 5 4 zx zx 1 2 0 0 0 0 0 2 ð10Þ where [E], E and are the elasticity matrix, Young’s modulus and Poisson’s ratio respectively.
3.1 The Elastic Stiffness Matrix The elastic strain energy e induced by elastic stresses in an element e is ððð 1 f"gT fgdxdydz e ¼ 2
ð11Þ
Formulation of a Three-dimensional Numerical Manifold Method
Substituting (8) into (11), we get ððð 1 e ¼ fDe gT ½Be T ½E½Be fDe gdxdydz 2 ððð 1 ¼ fDe gT fDe g ½Be T ½E½Be dxdydz for tetrahedron 2 ð1 ð1 ð1 1 ½Be T ½E½Be jJjddd for hexahedron ¼ fDe gT fDe g 2 1 1 1
607
ð12aÞ ð12bÞ
The element stiffness matrix for tetrahedron or hexahedron element is then given as: Ve ½BeðrÞ ½E½BeðsÞ ! ½KeðrÞeðsÞ ; r; s ¼ 14 for tetrahedron and 18 for hexahedron ð13Þ V for tetrahedron is given in the appendix while V for hexahedron is evaluated by numerical integration.
3.2 The Initial Stress Matrix Initial stress in element e can be expressed as fe 0 g ¼ f x 0
y 0
z 0
xy 0
yz 0
zx 0 gT
The energy induced by the initial stress is given by ððð ¼ f"e gT fe 0 gdxdydz ððð fDe gT ½Be T fe 0 gdxdydz ¼ ððð T ¼ fDe g ½Be T dxdydzfe 0 g for tetrahedron ððð ½Be T jJjdddfe 0 g for hexahedron ¼ fDe gT
ð14Þ
ð15aÞ ð15bÞ
The initial stress sub-matrix for tetrahedron can then be obtained explicitly as follows: 8 09 x > 9> 8 > > T > > ½Beð1Þ > > > 0> > > > > > = => < y0 > < ½B T > z T eð2Þ 0 for tetrahedron ð16Þ Ve ½Be fe g ¼ Ve xy 0 > > > ½Beð3Þ T > > > > > > > > > ;> yz 0 > : > > ½Beð4Þ T > ; : 0> zx From (16), we have Ve ½BeðrÞ T fe 0 g ! fFeðrÞ g;
r ¼ 1; 2; 3; 4
ð17aÞ
608
Y. M. Cheng and Y. H. Zhang
For hexahedron, numerical integration is required to form the initial-stress sub-matrix which is given by ð1 ð1 ð1 ½BeðrÞ T jJjdddf0e g ! fFeðrÞ g; r ¼ 18 for hexahedron ð17bÞ 1 1 1
3.3 The Point Loading Matrix For a point force FðFx ; Fy ; Fz Þ acting on the point ðX0 ; Y0 ; Z0 Þ in an element e, the potential energy induced by F is F ¼ Fx u Fy Fz w 8 9 8 9 > > = = < Fx > < Fx > ð18Þ T T Fy ¼ ð u w Þ Fy ¼ fDe g ½Te ðx0 ; y0 ; z0 Þ > > > > : ; : ; Fz Fz The corresponding point loading matrix for tetrahedron and hexahedron elements is given by 8 9 < Fx = ½TeðrÞ ðx0 ; y0 ; z0 ÞT Fy ! FeðrÞ ; r ¼ 14 or 18 ð19Þ : ; Fz 3.4 The Body Force Matrix If there is a body force GðGx ; Gy ; Gz Þ acting on element e, the corresponding potential energy is ððð G ¼ ðGx u þ Gy þ Gz wÞdxdydz 8 9 > ððð = < Gx > ð u w Þ Gy dxdydz ¼ > ; : > Gz 8 9 > ððð = < Gx > T T ¼ ½De for tetrahedron ð20aÞ ½Te ðx; y; zÞ dxdydz Gy > ; : > Gz 8 9 > ð1 ð1 ð1 = < Gx > T T for hexahedron ¼ fDe g ½ Teð1Þ Teð2Þ Teð8Þ jJjddd Gy > 1 1 1 ; : > Gz ð20bÞ The body force matrix for tetrahedron is given as 9 8 3 2 T > > ½T ðx; y; zÞ 8 9 8 9 > > eð1Þ > > > > 7> Gx > 6 ððð > = = 6ððð > < ½T ðx; y; zÞT > < Gx > 7< = eð2Þ 7 6 dxdydz7 Gy ½Te ðx; y; zÞT dxdydz Gy ¼ 6 T> > > > > 7> 6 ½T ðx; y; zÞ > > : ; 4 eð3Þ > > 5: Gz ; Gz > > > ; : ½T ðx; y; zÞT > eð4Þ
ð21aÞ
Formulation of a Three-dimensional Numerical Manifold Method
609
For hexahedron, the body force matrix is determined by numerical integration as 8 9 2 3 ð1 ð1 ð1 0 0 VeðiÞ < Gx = T VeðiÞ 0 5 ½TeðiÞ jJjddd Gy ¼ 4 0 ð21bÞ : ; 1 1 1 0 0 V Gz eðiÞ For Eq. (21a), the integration is evaluated as ððð ½TeðiÞ T dxdydz 2 3 0 0 ððð fi1 þ fi2 x þ fi3 y þ fi4 z 6 7 ¼ 0 fi1 þ fi2 x þ fi3 y þ fi4 z 0 4 5dxdydx 0 0 fi1 þ fi2 x þ fi3 y þ fi4 z 2 3 0 VeðiÞ 0 6 7 ð22Þ ¼ 4 0 VeðiÞ 0 5 0 0 VeðiÞ where VeðiÞ ¼ fi1 V e þ fi2 Vxe þ fi3 Vye þ fi4 Vze VeðiÞ ¼
ð1 ð1 ð1
NeðiÞ jJjddd
for tetrahedron for hexahedron
ð23aÞ ð23bÞ
1 1 1
The detailed expressions V e ; Vxe ; Vye ; Vze for tetrahedron are given in the appendix and will not be repeated here. The body force matrix is then given by 2 38 9 VeðrÞ 0 0 < Gx = 4 0 VeðrÞ 0 5 Gy ! FeðrÞ ; r ¼ 14 or 18 ð24Þ : ; 0 0 VeðrÞ Gz 3.5 The Inertia Force Matrix Let ð u w Þ be the time-dependent displacements of any point in element e, the force of inertia in unit volume is 8 9 8 9 uðtÞ = < fx = @2 < @ 2 fDe ðtÞg fy ¼ 2 ðtÞ ¼ ½Te ð25Þ : ; ; @t2 @t : wðtÞ fz where is the mass of unit volume. The potential energy of the inertial force in an element is 8 9 > ððð = < fx > i ¼ ð u w Þ fy dxdydz > ; : > ð26Þ fz ððð 2 2 ¼ fDe gT ½Te T ½Te dxdydz fDe g fVe ð0Þg 2
610
Y. M. Cheng and Y. H. Zhang
where is the time step, fVe ð0Þg is the speed at initial time step. From (26), we have ððð 2 T ½T ½T dxdydz ! ½KeðrÞeðsÞ ; r; s ¼ 14 or 18 ð27Þ eðrÞ eðsÞ 2 ððð 2 ½TeðrÞ T ½TeðsÞ dxdydz fVeðsÞ ð0Þg ! fFeðrÞ g; r; s ¼ 14 or 18 ð28Þ where 8 9 u ð0Þ = @ < eðsÞ eðsÞ ð0Þ fVeðsÞ ð0Þg ¼ @t : w ð0Þ ;
ð29Þ
eðsÞ
The matrix integration can be performed as 2 3 0 0 ððð ððð NeðrÞ NeðsÞ 6 7 0 NeðrÞ NeðsÞ 0 ½TeðrÞ T ½TeðsÞ dxdydz ¼ 4 5dxdydz 0 0 NeðrÞ NeðsÞ 0 1 0 0 teðrÞeðsÞ B teðrÞeðsÞ 0 C ¼@ 0 A 0 0 teðrÞeðsÞ
ð30Þ
For tetrahedron, teðrÞeðsÞ is given by ððð NeðrÞ NeðsÞ dxdydz teðrÞeðsÞ ¼ ððð ¼ ðfr1 þ fr2 x þ fr3 y þ fr4 zÞðfs1 þ fs2 x þ fs3 y þ fs4 zÞdxdydz
ð31aÞ
¼ fr1 fs1 V e þ ðfr1 fs2 þ fs1 fr2 ÞVxe þ ðfr1 fs3 þ fs1 fr3 ÞVye e e þ ðfr1 fs4 þ fs1 fr4 ÞVze þ fr2 fs2 Vxx þ fr3 fs3 Vyy þ fr4 fs4 Vzze e e e þ ðfr2 fs3 þ fs2 fr3 ÞVxy þ ðfr2 fs4 þ fs2 fr4 ÞVxz þ ðfr3 fs4 þ fs3 fr4 ÞVyz
For hexahedron, teðrÞeðsÞ is given by numerical integration as ð1 ð1 ð1 NeðrÞ NeðsÞ jJjddd teðrÞeðsÞ ¼
ð31bÞ
1 1 1 e e e e ; Vyze ; Vzx ; Vxx ; Vyy ; Vzze for tetrahedron are given in the The details of V e ; Vxe ; Vye ; Vze ; Vxy appendix. The inertia force matrix is hence given by 2 3 t 0 0 2 4 eðrÞeðsÞ 0 teðrÞeðsÞ 0 5 ! ½KeðrÞeðsÞ ; r; s ¼ 14 or 18 ð32Þ 2 0 0 t eðrÞeðsÞ
2 t 2 4 eðrÞeðsÞ 0 0
0 teðrÞeðsÞ 0
0 0 teðrÞeðsÞ
3 5fVeðsÞ ð0Þg ! fFeðrÞ g;
r; s ¼ 14 or
18 ð33Þ
Formulation of a Three-dimensional Numerical Manifold Method
611
3.6 The Fixed Point Matrix To simulate a fixed point ðx0 ; y0 ; z0 Þ within the discontinuous medium, a spring with a large spring stiffness p is usually adopted. The potential energy of the constraint spring is given by p f ¼ ðuðx0 ; y0 ; z0 Þ2 þ ðx0 ; y0 ; z0 Þ2 þ wðx0 ; y0 ; z0 Þ2 Þ 2 p ¼ fDe gT ½Te ðx0 ; y0 ; z0 ÞT ½Te ðx0 ; y0 ; z0 ÞfDe g ð34Þ 2 The stiffness matrix for this spring in tetrahedron is given by p½Te ðx0 ; y0 ; z0 ÞT ½Te ðx0 ; y0 ; z0 Þ 9 8 T > > ðx ; y ; z Þ ½T 0 0 0 > > eð1Þ > > > > > = < ½T ðx0 ; y0 ; z0 ÞT > eð2Þ ¼p > ½Teð3Þ ðx0 ; y0 ; z0 ÞT > > > > > > > > ; : ½T ðx ; y ; z ÞT > eð4Þ
0
0
0
ð½Teð1Þ ðx0 ; y0 ; z0 Þ ½Teð2Þ ðx0 ; y0 ; z0 Þ ½Teð3Þ ðx0 ; y0 ; z0 Þ ½Teð4Þ ðx0 ; y0 ; z0 ÞÞ ð35Þ 0
1 pNeðrÞ ðx0 ; y0 ; z0 ÞNeðsÞ ðx0 ; y0 ; z0 Þ@ 0 0
0 1 0
1 0 0 A ! ½KeðrÞeðsÞ ; 1
r; s ¼ 1; 2; 3; 4
ð36aÞ
For hexahedron element, the isoparametric coordinates ð0 ; 0 ; 0 Þ corresponding to ðx0 ; y0 ; z0 Þ are calculated by an iterative method as discussed in a later section of the paper. The stiffness matrix for this spring is given by 0 1 1 0 0 pNeðrÞ ð0 ;0 ;0 ÞNeðsÞ ð0 ;0 ;0 Þ@ 0 1 0 A ! ½KeðrÞeðsÞ ; r;s ¼ 1;2;3;4;5;6;7;8 ð36bÞ 0 0 1 4. 3D Contact Detection 3D cover contact detection in NMM is much more complicated than the corresponding 2D NMM problem in the geometry evaluation. For 3D NMM analysis, existing 3D contact detection schemes cannot be used directly because the contact stiffness matrix is not clearly defined before each step. The 3D contact relation is not well defined unless the displacement within each step is defined. To deal with this difficulty, the authors propose a simple iterative technique. For the tetrahedron as shown in Fig. 2a, suppose P1 is the vertex of element i, P2 P3 P4 is the contact face of element j, ðxk ; yk ; zk Þ and ðuk ; k ; wk Þ are the coordinates and displacements of point Pk ðk ¼ 1; 2; 3; 4Þ, the normal distance from point P1 to face P2 P3 P4 is 1 x 1 y 1 z1 1 1 x2 y2 z2 ð37Þ dn ¼ 2S 1 x3 y3 z3 1 x 4 y 4 z4
612
Y. M. Cheng and Y. H. Zhang
Fig. 2. a) Shear and normal displacements for tetrahedron, b) the contact relationship of point P1and face P2P3P5P4 for hexahedron
where S is the area of triangle P2 P3 P4 , and 1 x2 1 S ¼ 1 x3 2 1 x4 After a movement of ðuk ; k ; wk Þ, 1 1 1 dn 0 ¼ 2S 1 1
y2 y3 y4
P1 moves to P0 and the normal distance is x1 þ u1 y1 þ 1 z1 þ w1 x2 þ u2 y2 þ 2 z2 þ w2 x3 þ u3 y3 þ 3 z3 þ w3 x4 þ u4 y4 þ 4 z4 þ w4
ð38Þ
ð39Þ
If dn 0 < 0, the point P1 will penetrate the face of P2 P3 P4 , which means a penetration takes place. As shown in Fig. 2b, the contact face for hexahedron has four vertexes: P2 ; P3 ; P4 ; P5 . In fact, one face can be defined by three vertexes while the fourth vertex can be considered as redundant in the contact detection for hexahedron element. We can take any three vertexes out of the four vertexes to form a triangular face and to calculate the contact sub-matrices. For example, we can take P2 ; P3 ; P4 in consideration and Eqs. (37)–(39) can be used to determine the contact relation. For each face of a hexahedron as shown in Fig. 3, two triangles can be formed and this will be checked before the above procedure is applied. This simple technique can greatly simplify the various relations associated with contact and will be adopted in the following sections. If a contact is modeled as a spring along the incremental displacement directions, numerical computation will be simpler but the constitutive behaviour of the joint
Fig. 3. Shear displacement and normal displacement for hexahedron
Formulation of a Three-dimensional Numerical Manifold Method
613
cannot be considered as the joint is not a one-dimensional spring. In the present paper, the spring at the contact is decomposed into normal and shear springs and the Mohr– Coulomb law of friction is applied which is an extension of the approach proposed by the authors for DDA formulation (1998). The present approach can be extended to more complicated stress-strain relations at contacts easily and is considered to be crucial in extending NMM to more general conditions. Shi (1991) has adopted a onedimensional spring to represent the contact relation which is actually far from realistic and is not able to reflect the constitutive relations of the joints.
4.1 The Normal Contact Matrix If the second order terms in Eq. (39) are neglected, which is acceptable if the time step is small, we can get 8 9 8 9 u1 > > > = = < < u2 > 1 1 0 dn dn þ ð L11 L12 L13 Þ 1 þ ð L21 L22 L23 Þ 2 > > 2S ; 2S ; : > : > w1 w2 8 9 8 9 ð40Þ u3 > u4 > > > = = < < 1 1 þ ð L31 L32 L33 Þ 3 þ ð L41 L42 L43 Þ 4 > > > 2S 2S ; : ; : > w3 w4 where 1 L11 ¼ 1 1 1 L22 ¼ 1 1 1 L33 ¼ 1 1
1 x2 y2 z2 y3 z3 ; L12 ¼ 1 x3 1 x4 y 4 z4 1 x1 x1 z1 x3 z3 ; L23 ¼ 1 x3 1 x4 x 4 z4 1 y1 x1 y1 x2 y2 ; L41 ¼ 1 y2 1 y3 x4 y4
1 z2 z3 ; L13 ¼ 1 1 z4 1 y1 y3 ; L31 ¼ 1 1 y4 1 z1 z2 ; L42 ¼ 1 1 z3
1 y1 x2 y2 x3 y3 ; L21 ¼ 1 y3 1 y4 x4 y4 1 x1 y1 z1 y2 z2 ; L32 ¼ 1 x2 1 x4 y 4 z4 1 x1 x1 z1 x2 z2 ; L43 ¼ 1 x2 1 x3 x 3 z3
z1 z3 ; z4 z1 z2 ; z4 y1 y2 y3 ð41Þ
Set 8 9 < L11 = 1 fHg ¼ ½Ti ðx1 ; y1 ; z1 ÞT L12 : ; 2S L13
for tetrahedron
ð42aÞ
8 9 < L11 = 1 T fHg ¼ ½Ti ð1 ; 1 ; 1 Þ for hexahedron L : 12 ; 2S L13
ð42bÞ
614
Y. M. Cheng and Y. H. Zhang
8 9 8 9 8 9 L21 > L31 > > > > = = = < < 1 1 1 T T T fGg ¼ ½Tj ðx2 ;y2 ;z2 Þ L22 þ ½Tj ðx3 ;y3 ;z3 Þ L32 þ ½Tj ðx4 ;y4 ;z4 Þ L42 > > > 2S ; 2S ; 2S ; : > : > : > L23 L33 L43 ð43aÞ 8 9 8 9 8 9 L21 > L31 > > > > = = = < < 1 1 1 T T T fGg ¼ ½Tj ð2 ;2 ;2 Þ L22 þ ½Tj ð3 ;3 ;3 Þ L32 þ ½Tj ð4 ;4 ;4 Þ L42 > > > 2S ; 2S ; 2S ; : > : > : > L23 L33 L43 ð43bÞ Then we have dn 0 ¼ dn þ fHgT fDi g þ fGgT fDj g
ð44Þ
The potential energy of the normal spring is 1 s ¼ kn dn 02 2 p ¼ ðdn þ fHgT fDi g þ fGgT fDj gÞ2 2
ð45Þ
where kn is the normal spring constant. So we have kn fHiðrÞ gfHiðsÞ gT ! ½KiðrÞiðsÞ ;
r; s ¼ 14
or
18
ð46Þ
kn fHiðrÞ gfGjðsÞ gT ! ½KiðrÞjðsÞ ;
r; s ¼ 14
or
18
ð47Þ
kn fGjðrÞ gfHiðsÞ gT ! ½KjðrÞiðsÞ ;
r; s ¼ 14
or
18
ð48Þ
kn fGjðrÞ gfGjðsÞ gT ! ½KjðrÞjðsÞ ;
r; s ¼ 14
or
18
ð49Þ
kn dn fHiðrÞ g ! ½FiðrÞ ;
r ¼ 14
or
18
ð50Þ
kn dn fGjðrÞ g ! ½FjðrÞ ;
r ¼ 14
or
18
ð51Þ
4.2 The Shear Contact Matrix With reference to Fig. 2a and b, the spatial shear direction on the face P2 P3 P4 of points P1 and P0 cannot be determined easily, which is very different from the corresponding 2D problems where the shear direction is simply along the contact edge. Furthermore, the calculation of shear displacement is also very different from the corresponding 2D NMM. In the present study, vector theory is used for such a purpose. As shown in Fig. 3, suppose P0 ðx0 ; y0 ; z0 Þ is a point on face P2 P3 P4 and is the contact point for vertex P1 . After one time step, P0 and P1 move to P0 0 and P1 0 . The normal distance from P1 0 to P2 P3 P4 is dn 0 and the shear displacement on the face
Formulation of a Three-dimensional Numerical Manifold Method
615
P2 P3 P4 of point P1 and P0 is d 0 . Then we have ! d 02 ¼ ðjp0 0 p1 0 jÞ2 dn 02 ¼ ðx1 þ u1 x0 u0 Þ2 þ ðy1 þ 1 y0 0 Þ2 þ ðz1 þ w1 z0 w0 Þ2 dn 02 ¼ ðx1 x0 Þ2 þ ðy1 y0 Þ2 þ ðz1 z0 Þ2 þ 2ðx1 x0 Þðu1 u0 Þ þ 2ðy1 y0 Þð1 0 Þ þ 2ðz1 z0 Þðw1 w0 Þ þ ðu1 u0 Þ2 þ ð1 0 Þ2 þ ðw1 w0 Þ2 dn 02 ð52Þ After elimination of second and higher order terms, we get 0 1T 8 9 0 18 9 x1 x0 > x1 x0 > = = < u1 > < u0 > B C B C d 02 S1 þ 2@ y1 y0 A 1 2@ y1 y0 A 0 dn 02 > > ; ; : > : > z1 z0 w1 z1 z0 w0
ð53Þ
¼ S1 þ 2fDi gT fH 0 g þ 2fDj gT fG0 g ðfHgT fDi g þ fGgT fDj g þ dn Þ2 where S1 ¼ ðx1 x0 Þ2 þ ðy1 y0 Þ2 þ ðz1 z0 Þ2 8 9 < x1 x0 = fH 0 g ¼ ½Ti ðx1 ; y1 ; z1 ÞT y1 y0 : ; z1 z0 8 9 < x1 x0 = T fG0 g ¼ ½Tj ðx0 ; y0 ; z0 Þ y y0 : 1 ; z1 z0
ð54Þ ð55Þ
ð56Þ
The potential energy of the shear spring is 1 s ¼ ks d 02 2 p P ¼ ðS1 þ 2fDi gT fH 0 g þ 2fDj gT fG0 gÞ ðdn þ fHgT fDi g þ fGgT fDj gÞ2 ð57Þ 2 2 where ks is the shear spring. So the shear spring stiffness is given by ks fHiðrÞ gfHiðsÞ gT ! ½KiðrÞiðsÞ ;
r; s ¼ 14
or
18
ð58Þ
ks fHiðrÞ gfGjðsÞ gT ! ½KiðrÞjðsÞ ;
r; s ¼ 14
or
18
ð59Þ
ks fGjðrÞ gfHiðsÞ gT ! ½KjðrÞiðsÞ ;
r; s ¼ 14
or
18
ð60Þ
ks fGjðrÞ gfGjðsÞ gT ! ½KjðrÞjðsÞ ;
r; s ¼ 14
or
18
ð61Þ
ks dn fHiðrÞ g pfHiðrÞ 0 g ! ½FiðrÞ ;
r ¼ 14
or
18
ð62Þ
ks dn fGjðrÞ g pfGjðrÞ 0 g ! ½FjðrÞ ;
r ¼ 14
or
18
ð63Þ
616
Y. M. Cheng and Y. H. Zhang
Fig. 4. Direction of friction force and friction displacement
4.3 The Friction Force Matrix Coulomb’s law is used to model the stress-strain relation along contacts. The shear contact is considered to yield if F N tanð Þ þ C
ð64Þ
where N is the normal force on the contact surface, and C are the friction angle and cohesion force over the area of the contact respectively. The friction force is given by the contact spring as F ¼ pd tanð Þ
ð65Þ
where p is the stiffness of normal spring, d is the normal penetration distance, tanð Þ is the friction coefficient. 3D contact detection is very complicated for general problems and the authors (2002) have proposed a relatively simple but fast method for general conditions. The penetration edge method by Cheng (2002, 2005) is robust in the analysis but is applied after the first iteration. Before the first iteration and hence any contact detection, the direction of friction force is estimated in the following way. As shown in Fig. 4, suppose the contact point P1 ðx1 ; y1 ; z1 Þ moves to P01 ðx1 þ u; y1 þ v; z1 þ wÞ, where ~ n is the normal direction of the face P2 P3 P4 . The shear vector on the face P2 P3 P4 from point P1 to P01 is ! ! nÞ~ n ð66Þ p1 p1 0 ðp1 p1 0 ~ The equation of face P2 P3 P4 is x y x2 y2 x3 y3 x4 y4
z z2 z3 z4
1 1 ¼0 1 1
ð67Þ
which can be expressed as ax þ by þ cz ¼ d where y2 z2 1 a ¼ y 3 z3 1 ; y 4 z4 1
x 2 z2 1 b ¼ x3 z3 1 ; x 4 z4 1
x2 y2 1 c ¼ x3 y3 1 ; x4 y4 1
ð68Þ x2 y2 z2 d ¼ x3 y3 z3 x4 y4 z4 ð69Þ
Formulation of a Three-dimensional Numerical Manifold Method
617
It is easy to get a b c ~ n ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi~i þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi~j þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ k 2 2 2 2 2 2 2 a þb þc a þb þc a þ b2 þ c 2 k ¼ a1~i þ b1~j þ c1~
ð70Þ
! P1 P0 1 ¼ u1~i þ v1~j þ w1~ k
ð71Þ
! ! n ~ n ¼ ðu1 a21 u1 a1 b1 v1 a1 c1 w1 Þ~i P1 P1 0 P1 P 1 0 ~ þ ðv1 a1 b1 u1 b21 v1 b1 c1 w1 Þ~j þ ðw1 a1 c1 u1 b1 c1 v1 c21 w1 Þ~ k 9T 8 9 8 9T 8 9 8 9T 8 9 8 2 > = > = > = > = > = > = < 1 a1 > < u1 > < a1 b1 > < u1 > < a1 c1 > < u1 > 2 ~ ~ k ð72Þ v 1 i þ 1 b1 v1 j þ b1 c1 v1 ~ ¼ a1 b1 > > > > > > ; > ; > ; > ; > ; : : > : : > : 2; : 1 c1 w1 b1 c1 w1 w1 a1 c1 Then we can get the displacement df and potential energy f qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! ! nÞ~ n ¼ ðu21 þ v21 þ w21 Þ ða1 u1 þ b1 v1 þ c1 w1 Þ2 df ¼ P1 P1 0 ðP1 P1 0 ~ f ¼ Fdf
ð73Þ ð74Þ
We cannot obtain the expression and value of @f =@fDi g directly because f contains the radical signs of the unknowns ðu1 v1 w1 Þ. For the corresponding 2D condition, @f =@fDi g can be determined explicitly and an iterative method is not required. It is very difficult to pre-define the direction of the friction force and to establish the corresponding shear friction matrix in general 3-D problems. A simple iterative method which is adequate if the time step is not large is developed below. Assume dx ; dy ; dz are the friction displacement vector df at the direction of x, y, z, we have qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð75Þ df ¼ dx2 þ dy2 þ dz2 When Eq. (64) is fulfilled, the shear constraint is relieved. In the initial step, assume friction force is equal to zero, after solving the global equilibrium equations which do not contain the friction force matrix, we can get the displacements ðu1 ; v1 ; w1 Þ. The iterative procedures are: Step I: Substituting ðu1 ; v1 ; w1 Þ into Eq. (72), we can get the friction displacement vector ! ! P1 P 0 P1 P1 0 ~ n ~ n ¼ dx~i þ dy~j þ dz~ k ð76Þ and the direction of the friction force is 1 ðdx dy dz Þ df
ð77Þ
618
Y. M. Cheng and Y. H. Zhang
So the potential energy of friction force F at the point P1 in element e is 8 9 > = < dx > F f ¼ Fdf 0 ¼ ðu1 v1 w1 Þ dy > df ; : > dz 8 9 > = < dx > F T T ¼ fDi g ½Ti ðx1 ; y1 ; z1 Þ dy ¼ FfDi gT fHg > df ; : > dz where
9 8 Hið1Þ > > > =
ð78Þ
8 9 < dx = 1 ið2Þ ¼ ½Ti ðx1 ; y1 ; z1 ÞT dy fHg ¼ Hið3Þ > > : ; df > > dz ; : Hið4Þ
ð79Þ
FfHiðrÞ g ! fFiðrÞ g;
ð80Þ
So we have r ¼ 14
or
18
Similar equation can be established for point P0 in element j, so we have FfGjðrÞ g ! fFjðrÞ g;
r ¼ 14
or
18
ð81Þ
where 9 8 Gjð1Þ > > > =
8 9 < dx = 1 jð2Þ ¼ ½Tj ðx0 ; y0 ; z0 ÞT dy fGg ¼ Gjð3Þ > > : ; df > > dz ; : Gjð4Þ 9 8 Gjð1Þ > > > > > =
8 9 < dx = 1 ¼ ½Tj ð0 ; 0 ; 0 ÞT dy fGg ¼ .. > . > : ; df > > dz > > ; : Gjð8Þ
for tetrahedron
ð82aÞ
for hexahedron
ð82bÞ
jð2Þ
Step II: Solving the new global equilibrium equations with Eqs. (80) and (81), we can get the new displacements ðu1 0 ; v1 0 ; w1 0 Þ of point P1 . Check the following constraint rules: 0 u1 u1 100% < 1% u 1 0 v1 v1 100% < 1% ð83Þ v 1 0 w1 w1 100% < 1% w 1 If Eq. (83) is not fulfilled, go to step I, set ðu1 ; v1 ; w1 Þ ¼ ðu1 0 ; v1 0 ; w1 0 Þ, begin a new cycle. If Eq. (83) is fulfilled, stop this cycle.
Formulation of a Three-dimensional Numerical Manifold Method
619
5. Mapping between (x; y; z) and (n; g; f) for Hexahedron The mapping from isoparametric coordinates to global coordinates for hexahedron is given by 8 X Ni ð; ; Þxi x¼ i¼1
y¼
8 X
Ni ð; ; Þyi
ð84Þ
i¼1
z¼
8 X
Ni ð; ; Þzi
i¼1
From Eq. (84), it is very easy to calculate the corresponding ðx; y; zÞ of any point ð; ; Þ. On the contrary, it is difficult to express ð; ; Þ directly by using ðx; y; zÞ. In fact, we need to use an iterative method to calculate the corresponding ð; ; Þ of any point ðx; y; zÞ. The three dimensional Jacobian transformation matrix is given by 3 2 @x @y @z 6 @ @ @ 7 7 6 6 @x @y @z 7 7 ð85Þ J¼6 6 @ @ @ 7 7 6 4 @x @y @z 5 @ @ @ so we have
8 9 8 9 < dx = < d = dy ¼ J d : ; : ; dz d
Equation (86) can be rewritten as the following vector expression: ! ! dx ¼ J d
ð86Þ
ð87Þ
Since ! ! d ¼ J 1 dx
ð88Þ ! ¼ ð; ; Þ;~ x ¼ ðx; y; zÞ; ! x0 ¼ ðx0 ; y0 ; z0 Þ, Set 0 ¼ ð0 ; 0 ; 0 Þ ¼ 0; ~ ! ! d ¼ ðd; d; dÞ; dx ¼ ðdx; dy; dzÞ: ! Step I: Substituting 0 into the following equations step by step, we can get a new value of ~ : ! x0 ¼
8 X
! Ni 0 xi
i¼1
! dx ¼ ~ x! x0 ! 1 ! d ¼ J dx ! ! ~ ¼ 0 þ d
ð89Þ
620
Y. M. Cheng and Y. H. Zhang
! Step II: (A) If the value of ~ 0 is greater than a small positive number which is ! preset by the user (for example: 0.0001, as adopted in present analysis), set 0 ¼ ~ and go to step I to continue the iteration process. (B) If the difference is smaller than the preset value, the iterative process will stop and the value of ~ corresponding to ~ x is obtained. It should be noted that this inverse mapping step is not required in tetrahedron cover functions as the explicit equation forms can be determined without the use of local coordinates. 6. Examples to Verify the Proposed Formulation In general, it is not easy to verify a discontinuous code as a closed- form solution is usually not available. In this section, the proposed method and the corresponding computer code are verified by considering first a continuous problem, where an analytical solution is available. For a simple plane strain 2D shallow foundation problem in an elastic medium, the authors (1998) have obtained the settlement from 2D DDA with an accuracy which is compared with both the analytical and UDEC solutions. For this analysis, triangle elements are for both 2D DDA code and UDEC while hexahedron elements with a narrow width are adopted in the 3D NMM analysis and the results are shown in Table 1. The normal and shear stiffness of the joint are : kn ¼ 1 108 Pa and ks ¼ 3 107 Pa; the Young’s modulus and Poisson’s ratio of the elastic medium are E ¼ 8.84 107 Pa and 0.25. The width of the uniformly applied loading is 2 m. For the second verification example, a simple discontinuous problem as shown in Fig. 5 is considered. The volume and weight of the sliding wedge are 20.833 m3 and Table 1. Settlement at the centre of a 2D shallow foundation, width of loading is 2 m Stress (Pa)
Analytical (mm)
2D DDA
2D UDEC
3D NMM
10000 20000 30000 40000 50000 60000
0.265 0.530 0.795 1.060 1.325 1.590
0.272 0.550 0.809 1.072 1.333 1.615
0.237 0.475 0.760 1.050 1.339 1.637
0.280 0.554 0.811 1.075 1.331 1.602
Fig. 5. Simple sliding wedge problem
Formulation of a Three-dimensional Numerical Manifold Method
621
416.67 kN, respectively. With reference to the material parameters of the joint as shown in Fig. 5, the factor of safety computed by the limit equilibrium method is found to be 0.726. With the adoption of a single tetrahedron element and the application of gravity load to the tetrahedron, the factor of safety is also found to be 0.726. The sliding angle of this block is located at 45 to the x-axis which is obviously the correct answer. 7. Numerical Examples to more Complicated Problems To illustrate the applicability of 3D NMM to more complicated problems, three numerical examples are reported in the following. In example 1, which is a simple block sliding problem, the physical mesh is shown in Fig. 6a, while the mathematical covers of the manifold are shown in Fig. 6b, the solution of the system is controlled by the actual mesh as shown in Fig. 6c. The bottom blocks are fixed so that the inclined top surface of the bottom block can simulate the sliding surface. The weights of the blocks are the driving forces in the sliding
Fig. 6. a) Physical mesh of slope, b) mathematical covers of manifold method, c) manifold elements by mathematical cover and physical mesh
Fig. 7. Sliding state after 40 time steps and principal stresses
Fig. 8. Sliding state after 60 time steps and principal stresses
622
Y. M. Cheng and Y. H. Zhang
Fig. 9. 3D-veiw of computing grid (each block is 10 m in width)
Fig. 10. Cutaway view of computing grid
process. The parameters for analysis are: Young’s modulus of block E ¼ 1 1010 Pa, Poisson’s ratio ¼ 0.34, friction angle ’ ¼ 25 , cohesion c ¼ 1 105 Pa, density ¼ 2 104 kg=m3, t ¼ 0.02s, kn ¼ 3 109 N=m2, ks ¼ 1 109 N=m2. Figures 7 and 8 give the results after 40 time steps and 60 time steps. The critical principal stresses and the displacement modes as shown in Figs. 7 and 8 appear to be reasonable. Example 2 is a coal extraction problem and the blocks are shown in Figs. 9 and 10. The material properties and parameters for the calculation are E ¼ 12 109 Pa, ¼ 0.34, ’ ¼ 25 , c ¼ 1 105 Pa, ¼ 2 104 kg=m3, t ¼ 0.02s, kn ¼ 3 109 N=m2, ks ¼ 1 109 N=m2. As shown in Figs. 9 and 10, the complete computing region is 200 m 60 m 60 m, while the region of coal seam is 100 m 60 m 10 m and is represented by many small blocks (the thickness of each seam is 2.5 m). The bottom coal seam was excavated in six steps with the length of each step being 10 m, which are the typical coal extraction sequences in China. The bottom of the domain is fixed in vertical direction while the four vertical sides of the domain are fixed against lateral displacements in x or y directions but not in the vertical direction. The results after 4 steps and 10 steps of excavation are shown in Figs. 11–14. In Figs. 11 and 14, the principal stresses are shown only for those blocks close to the excavation domain for clarity. In Fig. 11, the coal at the right hand side is extracted and there are significant stress changes and rotation of principal stress around the excavation which are clearly demonstrated. The vertical stresses at the blocks above and below the excavation domain are reduced to a very small value and this is a reasonable
Formulation of a Three-dimensional Numerical Manifold Method
Fig. 11. 3D-view and principal stresses after 4 step excavation
Fig. 12. Cutaway view after 4 step excavation
Fig. 13. 3D-view after 10 step excavation
Fig. 14. Cutaway view and principal stresses after 10 step excavation
623
624
Y. M. Cheng and Y. H. Zhang
Fig. 15. a) 3D view of computing grid, b) 3D view of coal mining collapse and principal stresses after 200 time step
prediction. In Fig. 14, the whole extraction process is completed and the top coal layer has dropped while upheaval is observed at the bottom of the excavation. This result is also in agreement with the expected trend of displacements during coal extraction. Example 3, as shown in Fig. 15a, is a coal mining problem where the joints are spaced irregularly and contact detection is more difficult than in the previous 2 examples. The parameters used for analysis are E ¼ 12 109 Pa, ¼ 0.34, ’ ¼ 20 , c ¼ 6 104 Pa, ¼ 2 104 kg=m3, t ¼ 0.03s, kn ¼ 3 1010 N=m2, ks ¼ 1 1010 N=m2. The applications of the boundary conditions are similar to the previous example. Blocks 1 and 2 are removed and the system starts to fail at about 40 time steps. The displacements and stresses at 200 time step as shown in Fig. 15b clearly demonstrate the failure mechanism which also appears to be reasonable. 8. Discussion and Conclusion For NMM to be useful in the solution of rock engineering problems, the extension of NMM to 3D analysis must be considered. In this paper, the theories of NMM based on
Formulation of a Three-dimensional Numerical Manifold Method
625
tetrahedron and hexahedron elements are developed. While the stiffness matrices for tetrahedron are formed explicitly by the vertices coordinates similar to that by Shi (1991) for 2D NMM, numerical integration of the matrices is necessary for hexahedron element to avoid lengthy equations. Combining these two kinds of cover functions, most of the practical engineering problems involving discontinuity can be modeled without difficulty. If the mathematical and physical covers are chosen to be the same, the present formulations or the corresponding 2D formulations can be viewed as practically a discrete finite element formulation. In this paper, all the sub-matrices and force vectors associated with tetrahedron or hexahedron elements are derived and these equations will be useful for numerical implementation of 3D NMM. The numerical examples discussed have demonstrated the applicability of 3D NMM and the results obtained are reasonable. Even though tremendous computer resources are required for 3D NMM, with the advance in computer technology, the authors expect that 3D NMM can become a useful technique in the future. In the present model, the contact spring is decomposed into normal and shear springs so that any constitutive model for the discontinuity can be adopted. A new problem is however introduced by such a decomposition which is not present in the original simple contact spring model by Shi (1991). Following the decomposition of contact spring into normal and shear directions, there is a major difficulty in determining the shear direction and friction direction for the 3-D contact problem which is not present in the corresponding 2D analysis. Because of this difficulty, the contact stiffness matrices cannot be determined directly. The authors proposed the use of vector theory and an iterative method which can overcome this difficulty. The simple contact spring model by Shi (1991) is not a realistic model of the joint behaviour but the limitation is overcome by the present proposal. A major outstanding problem in NMM will now be the constitutive relation of the contact joint which is largely controlled by the joint constitutive behaviour. The present formulations allow the incorporation of any complicated joint constitutive model as the yield criterion or the stress=strain relation can be easily modified. The proposed tetrahedron and hexahedron elements have been verified to be correct for simple problems where exact solutions are available. The authors have also applied these elements to more complicated problems without difficulty and the results obtained are shown to be reasonable. The formulation for tetrahedron and hexahedron covers can actually be combined in the same problem if necessary. The preliminary study by the authors has shown that such a combination is possible, but not recommended, because the computer coding of such an approach can be very complicated. The use of either cover function will be sufficient for most engineering applications. Based on the formulation as introduced in this paper, the use of higher order elements should be relatively simple to implement. References Cheng, Y. M. (1998): Advancement and improvements in discontinuous deformation analysis. Comput. Geotech. 22(2), 153–163. Cheng, Y. M., Zhang, Y. H. (2000): Rigid body rotation and internal block discretization in DDA analysis. Int. J. Numer. Anal. Met. Geomech. 24, 567–578. Cheng, Y. M., Zhang, Y. H. (2002): Coupling of FEM and DDA methods. Int. J. Geomech. 2(4), 503–517.
626
Y. M. Cheng and Y. H. Zhang
Cheng, Y. M., Zhang, Y. H., Chen, W. S. (2002): Wilson non-conforming element in numerical manifold method. Commun. Numer. Meth. Engng. 18, 877–884. Cheng, Y. M., Chen, W. S., Ger, X. R. (2002): Procedure to the detect of three-dimensional blocks using penetration edges method. Geotechnical Special Publication no. 117, ASCE, 79–85. Cheng, Y. M., Chen, W. S., Zhang, Y. H. (2006) New approach to determine three dimensional contacts. In blocks system – penetration edges method. Int. J. Geomech. ASCE (accepted for publication) Cundall, P. A., Hart, R. D. (1985): Development of generalized 2-D and 3-D distinct element programs for modeling jointed rock. Itasca Consulting Group. Misc. Paper SL-85-1. U.S. Army Corps of Engineering. Jiang, Q. H. (2000): Research on three dimensional discontinuous deformation analysis method, PhD Dissertation. Wuhan Institute of Rock & Soil Mechanics. Jiao, Y. Y. (1998): Three dimensional DDA and Its Application. PhD Dissertation. Wuhan Institute of Rock & Soil Mechanics. Shi, G. H. (1991): Manifold method of material analysis. Trans. 9th Army Conf. On Applied Mathematics and computing, Minneapolis, Minnesota, 57–76. Shi, G. H. (1993): Block system modeling by discontinuous deformation analysis. Comput. Mech. Publications. Shi, G. H. (1996): Manifold method. Proc. of 1st Int. Forum on DDA and Simulations of Discontinuous Media. Berkeley, California, USA, 52–204. Shi, G. H. (1997): Numerical manifold method. Proc. of 2nd International Conference on Analysis of Discontinuous Deformation, 1–35. Wang, R. L., Chen, N. M., Liu, B. S. (1996): Theory analysis of 3D-DDA. Chin. J. Rock Mech. Engng. 15, 219–224.
Appendix: Equations for Tetrahedron 1 x1 y 1 z1 1 1 x2 y2 z2 ¼ V 1 þ V2 þ V 3 þ V4 V ¼ 6 1 x3 y3 z3 1 x 4 y 4 z4 1 x1 y 1 z1 1 x 2 y 2 z2 ¼ 6V J ¼ 1 x 3 y 3 z3 1 x 4 y 4 z4
ð90Þ
ð91Þ
1 V1 ¼ ðx2 y3 z4 þ x3 y4 z2 þ x4 y2 z3 x4 y3 z2 x2 y4 z3 x3 y2 z4 Þ 6
ð92Þ
1 V2 ¼ ðx1 y4 z3 þ x3 y1 z4 þ x4 y3 z1 x1 y3 z4 x3 y4 z1 x4 y1 z3 Þ 6
ð93Þ
1 V3 ¼ ðx1 y2 z4 þ x2 y4 z1 þ x4 y1 z2 x1 y4 z2 x3 y1 z4 x4 y2 z1 Þ 6
ð94Þ
1 V4 ¼ ðx1 y2 z3 þ x2 y3 z1 þ x3 y1 z2 x1 y3 z2 x2 y1 z3 x3 y2 z1 Þ 6
ð95Þ
Formulation of a Three-dimensional Numerical Manifold Method
f11 ¼
V1 V
627
ð96Þ
1 f12 ¼ ðy2 z4 y3 z4 þ y4 z3 y2 z3 þ y3 z2 y4 z2 Þ J
ð97Þ
1 f13 ¼ ðz2 x4 z3 x4 þ z4 x3 z2 x3 þ z3 x2 z4 x2 Þ J
ð98Þ
1 f14 ¼ ðx2 y4 x3 y4 þ x4 y3 x2 y3 þ x3 y2 x4 y2 Þ J
ð99Þ
f21 ¼
V2 V
ð100Þ
1 f22 ¼ ðy1 z3 y4 z3 þ y3 z4 y1 z4 þ y4 z1 y3 z1 Þ J
ð101Þ
1 f23 ¼ ðz1 x3 z4 x3 þ z3 x4 z1 x4 þ z4 x1 z3 x1 Þ J
ð102Þ
1 f24 ¼ ðx1 y3 x4 y3 þ x3 y4 x1 y4 þ x4 y1 x3 y1 Þ J
ð103Þ
f31 ¼
V3 V
ð104Þ
1 f32 ¼ ðy1 z4 y2 z4 þ y4 z2 y1 z2 þ y2 z1 y4 z1 Þ J
ð105Þ
1 f33 ¼ ðz1 x4 z2 x4 þ z4 x2 z1 x2 þ z2 x1 z4 x1 Þ J
ð106Þ
1 f34 ¼ ðx1 y4 x2 y4 þ x4 y2 x1 y2 þ x2 y1 x4 y1 Þ J
ð107Þ
f41 ¼
V4 V
ð108Þ
1 f42 ¼ ðy1 z2 y3 z2 þ y2 z3 y1 z3 þ y3 z1 y2 z1 Þ J
ð109Þ
1 f43 ¼ ðz1 x2 z3 x2 þ z2 x3 z1 x3 þ z3 x1 z2 x1 Þ J
ð110Þ
1 f44 ¼ ðx1 y2 x3 y2 þ x2 y3 x1 y3 þ x3 y1 x2 y1 Þ J ð 1 Ve ¼ 1Dðx; y; zÞ ¼ J 6 P1 P2 P3 P4 ð 1 Vxe ¼ xDðx; y; zÞ ¼ J ðx1 þ x2 þ x3 þ x4 Þ 24 P1 P2 P3 P4
ð111Þ ð112Þ ð113Þ
628 Y. M. Cheng and Y. H. Zhang: Formulation of a 3-D Numerical Manifold Method
Vye ¼ Vze ¼ e Vxx
ð
1 ðy1 þ y2 þ y3 þ y4 Þ 24
ð114Þ
zDðx; y; zÞ ¼ J
1 ðz1 þ z2 þ z3 þ z4 Þ 24
ð115Þ
P1 P2 P3 P4
ð P1 P2 P3 P4
ð
¼
yDðx; y; zÞ ¼ J
x2 Dðx; y; zÞ ¼ J
P1 P2 P3 P4
1 ð2x1 x1 þ x1 x2 þ x1 x3 þ x1 x4 120
þ x2 x1 þ 2x2 x2 þ x2 x3 þ x2 x4 þ x3 x1 þ x3 x2 þ 2x3 x3 þ x3 x4 þ x4 x1 þ x4 x2 þ x4 x3 þ 2x4 x4 Þ e Vyy ¼
ð
y2 Dðx; y; zÞ ¼ J
P1 P2 P3 P4
1 ð2y1 y1 þ y1 y2 þ y1 y3 þ y1 y4 þ y2 y1 120
þ 2y2 y2 þ y2 y3 þ y2 y4 þ y3 y1 þ y3 y2 þ 2y3 y3 þ y3 y4 þ y4 y1 þ y4 y2 þ y4 y3 þ 2y4 y4 Þ Vzze ¼
ð
z2 Dðx; y; zÞ ¼ J P1 P2 P3 P4
ð
¼
xyDðx; y; zÞ ¼ J P1 P2 P3 P4
ð
xzDðx; y; zÞ ¼ J
P1 P2 P3 P4
¼
ð P1 P2 P3 P4
yzDðx; y; zÞ ¼ J
ð119Þ
1 ð2x1 z1 þ x1 z2 þ x1 z3 þ x1 z4 120
þ x2 z1 þ 2x2 z2 þ x2 z3 þ x2 z4 þ x3 z1 þ x3 z2 þ 2x3 z3 þ x3 z4 þ x4 z1 þ x4 z2 þ x4 z3 þ 2x4 z4 Þ Vyze
ð118Þ
1 ð2x1 y1 þ x1 y2 þ x1 y3 þ x1 y4 120
þ x2 y1 þ 2x2 y2 þ x2 y3 þ x2 y4 þ x3 y1 þ x3 y2 þ 2x3 y3 þ x3 y4 þ x4 y1 þ x4 y2 þ x4 y3 þ 2x4 y4 Þ Vxze ¼
ð117Þ
1 ð2z1 z1 þ z1 z2 þ z1 z3 þ z1 z4 120
þ z2 z1 þ 2z2 z2 þ z2 z3 þ z2 z4 þ z3 z1 þ z3 z2 þ 2z3 z3 þ z3 z4 þ z4 z1 þ z4 z2 þ z4 z3 þ 2z4 z4 Þ e Vxy
ð116Þ
ð120Þ
1 ð2y1 z1 þ y1 z2 þ y1 z3 þ y1 z4 120
þ y2 z1 þ 2y2 z2 þ y2 z3 þ y2 z4 þ y3 z1 þ y3 z2 þ 2y3 z3 þ y3 z4 þ y4 z1 þ y4 z2 þ y4 z3 þ 2y4 z4 Þ
ð121Þ
Author’s address: Y. M. Cheng, Department of Civil and Structural Engineering, Hong Kong Polytechnic University, Room TV707, Hung Hom Kowloon, Hong Kong, P.R. China; e-mail: [email protected]
