- Email: [email protected]
Numerical simulation of penetration in sand based on FEM
Computers and Geotechnics 9 (1990) 73-86
NUMERICAL SIMULATION OF PENETRATION IN SAND BASED ON FEM
Z. Sikora Institut for Soil and Rock Mechanics Technical University Karlsruhe West Germany and G. Gudehus Institut for Soil and Rock Mechanics Technical University Karlsruhe West Germany
ABSTRACT This article presents the application of the Finite Element Method (FEM) for estimating the penetration resistance and the stress and strain state in the soil around a penetrating rod. Numerical solutions are obtained assuming Lagrangian formulation with finite deformations. The Jauman terms are incorporated in t'he analysis. A constitutive law by Kolymbas is used. A condition representing skin friction has been defined for the surface of the penetrometer, and subsequently the unknown stress state on its surface is determined with an iterative method. The method is closely connected with the localization of deformation. Some numerical problems due to time-integration and convergence criteria for incrementally nonlinear constitutive laws are also discussed. In numerical computations, the friction coefficient parameter is varied. The influence of this parameter on the penetration resistance is discussed. INTRODUCTION The penetration test is widely used to investigate the mechanical properties of soil layers in the field. In this test, a probe with a conical tip is driven vertically into the soil, while the cone resistance and the total force are measured continuously. The cone resistance usually increases with depth. The measured cone resistance and shaft friction 73 Computers and Geotechnics0266-352X/90/S03.50 © 1990 Elsevier Science Publishers
Ltd, England. Printed in Great Britain
74 are used to estimate soil parameters, which are used in engineering practice. Much research has been done on the behaviour of soils during penetration. The information obtained from laboratory tests is limited because very little is known about the stress-strain state in the soil during penetration. Since stresses and strains are not uniformly distributed in this case, no techniques were available to estimate the stress distribution in the interior of the sample from the boundary conditions. The boundary condition for the surface of the penetrometer should be stated as an inequality describing the friction phenomenon. All techniques using interface-elements to describe the interaction between probe and soil are considered as subjective. Furthermore, it is dif~cult to define a constitutive relation for the soils; therefore introducing another constitutive law for the interface makes this problem even more obscure. The Coulomb friction condition can be used to describe an appropriate boundary condition. This paper presents new possibility of well-posedness of some boundary problems, where on a part of the boundary, a friction condition is valid. THE BOUNDARY PROBLEM FOR STATIC PENETRATION In the subsequent analysis, we consider a boundary value problem modelling a penetration process of a non-deformable penetrating rod into sand. Figure 1 shows the boundary conditions and a finite element mesh for the boundary problem which will be discussed here. Unless stated otherwise, upper case bold face letters are used to denote second order tensors while lower case bold face letters are used to denote vectors in R ~. To simplify the considerations, we have supposed that the penetrating rod is placed at a depth equal to the length of the penetrometer skin surface. An initial stress To is assumed to be known in the zone around the penetrometer. The lower part of the penetrating rod is taken as a hemisphere instead of a cone, which simplifies the calculations from a numerical point of view. Penetration experiments, carried out by Allersma [I], with different shapes of penetrometers showed that the shape of the penetrometer's tip has no major influence on the stress distribution in the granular material. The boundary condition on the surface of the penetrometer is defined as Coulomb friction, i.e. the magnitude of the tangential stress component must not be greater than the normal stress component multiplied b.v a friction coefficient/~. The friction coefficient/~ is assumed to be constant. The stress tensor T at any boundary point on the shaft of the penetrometer during the penetration process must fulfil the following condition on the surface of the penetrometer,
75
i" Boundary c o n d i t i o n s : AB :r.. = r . , = 0 B C :v~ = 0
and
Tr. = 0
CD:v,
--0
and
7",~ = 0
DE:v,
-- 0
and
7~,, = 0
EA : i f
I n T T t l < tJnTTn
v, = 0 else
and
then
v, = vm,
InTTtl =
pnTTn
b
-
¢
Figure 1: Boundary conditions and FE mesh for penetration problem.
InTTtl K pnTTn,
(1)
where vectors n and t are external normal and tangential vectors of the surface of the penetrometer respectively. T is the Cauchy stress tensor. The incremental equilibrlum will be described by the virtual work principle where time deviratives of Cauchy stress tensor, stretching, rate of surface traction and body force are components of the incremental work. The initial stress state To in the subsoil is the initial condition for the initial-boundary-value problem, and has to fulfil the condition (1). To determine the deformation of a body from the system of differential equations governing its mechanical deformations, we need side conditions such as boundary conditions. The remaining boundary conditions are those of rate of place vi(=,,,) (so-called Dirichlet boundary condition) and/or rate of traction (Neumann boundary condition), which are illustrated in Figure 1. To be able to solve the boundary problem with boundary condition in an inequality form, as defined in (1), a different numerical approach, in comparison to classical equilibrium problem with known boundary conditions, is needed. The main numerical problems due to this new approach will be presented in the sequel. The stress distribution is not known exactly on the boundary and is restricted by the friction condition. The Coulomb friction is defined with the assistance of the total stress, and not of its
76 increment. CONSTITUTIVE LAW The subsoil is assumed to be dry. The density of the sand is implemented in the parameters of the constitutive law. The tensor form of the constitutive law, which was proposed by Kolymbas [2] is as follows
"r = H ( T , D ) =
-• (TD+DT)+C2tr(TD)I+C3trTv/~+C4t--~ T2
(2)
is the rate of the corotational Cauchy stress tensor and Cl ..... C4 are known material constants evaluated on the basis of triaxial tests. 1 is the unit tensor. T and D denote Cauchy stress tensor and stretching. The material time derivative "~ of the Cauchy stress tensor T is assumed to be the Jaumann's type increment T minus T W - W T . i.e.
"~ = "[ + W T
- TW
(3)
where W denotes the spin tensor. There is a simple calibration method based on experimental data from a triaxial apparatus [2]. The results of this calibration can be used as input data to another numerical calibration method based on a random process within a nonlinear optimization program for obtaining a better approximation of the experimental data. The calibration process can be defined in general as a nonlinear optimization problem. The calibration is a very difficult problem from the mathematical viewpoint, specially for constitutive laws of the rate type. In this paper, the above mentioned calibration methods will be not presented because of the nature of the topic and size of this paper. As an example for this boundary value problem, we calibrated the four constitutive parameters C1,...,C4 to the triaxial experimental data with Oostershelde sand from Holland. The constitutive law was calibrated to the simple compression triaxial test with a confining pressure equal to 100 kPa. The obtained constitutive parameters are as follows
C, = -165.108,
C2 = 9.966,
C3 = 164.061,
C4 = -375.351.
(4)
The results of element test calculation are compared with triaxial compression test data as can be seen in Figure 2.
77
.8
..
.E O I~
.4
5.
/,
f
/o
5
•
s /
4-
4
Z
h
.<
3~ t~
4
t.J
4
W N
N
¢
Z C3
4
0.
n¢
-!
2
4 AXIAL
6
B
STRAIN
10
12
Ig)
Figure 2: Element test calculation in comparison to triaxial compression test M A I N P R O B L E M S OF N U M E R I C A L C A L C U L A T I O N S The presented constitutive law (2), which belongs to a large class of rate type constitutive laws, is very convenient from the numerical point of view. There is not, for instance, any artificial partition of deformation into elastic and plastic parts as with elasto-plastlc constitutive laws. The utilized constitutive law is a non-linear tensor function with resrect to the stretching D, i.e. the FEM equilibrium condition is a non-linear system of equations on the load increment level. In a Lagrangian incremental approach for large displacements we express the equilibrium of the body at a load increment number "i + 1" using the principle of virtual work as follows
K~,(x', Au'+~)~u~ ~ = gR~ ÷*
(5)
where K denotes a stiffness matrix (operator in the sense of FEM), x i is the actual configuration after load increment number ~i', ~ u I+I is the increment of displacement for the load increment A R I+1 within load increment number "i + 1". It is thought here that the problem becomes more complex as a non-linear system of equations must be solved for each load increment. The first problem consists of the choice of a particular numerical method for solving this non-linear system of equations. It is
78 clear that the amount of calculation for a non-linear system, for instance (5), is greater than in relation to linear ones (e.g. for load increment within an elasto-plastic constitutive relation). For the boundary problem of this paper a modified Newton-Raphson method was sucessfully used. Futhermore, the evaluation of the size of the convergence region of this non-linear problem is a very difficult problem. For the penetration problem it was necessary to introduce into the computer program a numerical method for enlarging the convergence region. Another fundamental difficulty in the general description of equations of equilibrium is that the configuration of the body for each load increment level "i + 1" is unknown. In our analysis, we consider the motion of the body in a fixed (stationary) Cartesian coordinate system where the configuration after each load increment is known, i.e.
x i+l = x i + A u i+1.
(6)
This problem with assumption of large deformation is treated as geometrically nonlinear, i.e. after Lagrange updating we obtain
j~+t
= Kt,,,(x,+,, & u , + t ) A u ~ , _ Kt.~(x', A u ' + ' ) A u ~ ' # 0.
(7)
where ji+l denotes a vector of correction forces which can be used as a load vector for the correction of solution accuracy in FE equilibrium equation as follows
K t , n ( x i , z ~ u i + l ) A u ~ 1 = AR~ +' - d}
and
do = O.
(s)
The presented method can be called a non-iterative correction method for geometric nonlinearity. As mentioned earlier, the constitutive relation is a non-linear tensor function of stretching, i.e. the boundary problem belongs to the group of physically non-linear problems. Introducing the Coulomb friction condition as boundary condition to the differential equations of equilibrium, which will be presented in the next section in detail, makes this problem aditionally non-linear, i.e. even in the case of a linear constitutive law and without assumption of large deformations, the boundary problem is formed in an inequality state (in the sense of boundary condition (1)) and remains therefore non-linear. Another problem of numerical solution concerns the time integration method. The solution of FE equilibrium state gives us an answer about distribution of the Cauchy stress
79 material derivative. For the penetration a second-order integration method (presentation of this method will be given in next section) is developed, which was very efficient specially for the problem considered. The last problem, which we will mention about numerical solution of penetration with utilizing the 'direct' boundary condition as friction condition, concerns of FE discretization. It is clear that because of large deformations around of the penetrating rod the mesh must be assumed finer as in the farther environment of the penetrometer, but very important for the convergence conditions of the non-linear equation system is that the discretizing nodes of penetrometer structure, which are also boundary points with friction condition, are placed in equal distances between them. A number of different calculations showed that a sufficient mesh width is about ten times the penetrometer radius, and the depth of the mesh can be equal to four times the penetrometer's length. FRICTIONAL BOUNDARY CONDITION A question arises: how does one take boundary friction conditions into account if instead of being defined 'explicitly' they are defined through their dependence on appropriate components of stress vector? As mentioned in the last section, the boundary problem with friction condition as boundary condition is a non-linear problem, and therefore the solving method of this problem is iterative. We are supposing the so-called kinematic load method, i.e. on the boundary upon which friction conditions are given are assumed as initial approximation of the friction condition the Dirichlet boundary conditions on the shaft of penetrometer. The penetrating rod is assumed as non-deformable, and the boundary conditions in the nodal points on the surface of penetrometer are as follows
v,=0
and
v,=v0.
(9)
The boundary conditions on the penetrometer shaft and on the other boundaries are fixed, and from the FE equilibrium condition (8) we can obtain velocity and stress rate fields. The current stress state and displacement field are obtained after time integration. The normal and tangential component of stress in nodal points on the surface of penetrometer caa now be calculated, and the condition (1) of friction can be checked (violated or not). In the case that the condition (1) is violated it means that the geometrical position of the point is wrong and must be corrected. Our task is to find such a configuration that in all boundary nodes the friction condition will be fulfilled. The
80
/':= t
=
(ex~a, rJ/,m)
~o := We = d
(,..,,,,,)
-
(0,,,,,)
=
~=,.',-,'+,,.o)
"_.,]_ __
Figure 3: Sliding conditions on penetrometer surface nodal point with violation of friction condition wiU be called sliding point or point in a sliding state. For the correct determination of geometrical position of sliding point we use, as additionaly equation of the FE non-linear equation system, the Coulomb friction condition in a form of equality. In Figure 3 we can see a detail of geometrical position of nodal point A0 which is assumed to be in a sliding state. The geometrical position of point A0 will be corrected to position A1 where the condition (1) is fulfilled with the equality symbol (Figure 3). We assume that the posible direction of nodal point motion on the shaft of penetrometer is tangential to the penetrometer surface. The geometrical relations are easily understood from the Figure 3. Thus the direction is defined in advance by the geometry of the penetrating rod structure. Boundary nodes themselves can only move along the surface of the penetrometer in accordance with the tangential vektor t at the skin surface. The proper correction vector wo, which starts at A0 and ends at the point A! for a particular nodal point, may be described as w0 = c . t with-c as a scalar. This scalar c is a further unknown in the node, which together with the unknown components of the rate vector Pl = (vrs, v,l) fulfils the following geometric equations
81 v,,
-
c. cosa
= 0
(10)
vz, + c - s i n a = v0.
(11)
The Coulomb friction equality is presented by the non-linear equation in relation to the stretching tensor, which via FEM reveals a non-linear vector rate function. Let us suppose, to simplify matters, that the new stress state, Tt, in the nodal point is defined a~
T,
= To + &ArT
(12)
where & is a time integration parameter, and At is the time step. The stress tensor TI may fulfil the following Coulomb friction equality, i.e.
InT(To + a A t T ) t ] = pnT(To + 6 A t T ) n
(13)
or in another form
- 1 7" n r H ( T 0 , D ) r = n r ' l ' r = ~-~-~n T0r where
;
r = t +
#sign(ro)n
(14)
r0 = nrTot.
The right hand side of the equation (14) depends on a known initial stress tensor To and time increment At. The left hand side of the equation (14) is the already mentioned non-linear function in relation to the rate field and becomes attached to the non-linear equation system resulting from the equilibrium conditions. For an arbitrary load increment we expect some problems with convergence because introducing the Coulomb friction equation into the global system of equations makes the condition number of the tangent stiffness matrix worse, and for this case, some methods for enlarging the convergence region must be applied. TIME INTEGRATION
METHOD
The method of time integration plays an important role in rate type constitutive laws using FEM. The solution using the appropriate load increment is obtained by integration of the load rate in relation to the time parameter. This time parameter in
82 relation to the time independent behaviour of the tested material is not perceived in its physical sense, but rather as a load parameter. Here we will present only the main idea of the integration method which lies in the class of methods with a global truncation error of 0(At3). The exact developement of this time integration method due to a rate type constitutive law with FEM will be presented in another paper [4]. The integration process takes place between to and to + 2xt, to obtain the new stress T" at time to + At, i.e.
t0+At
T" = T o +
/
H(T(r),D(r))dr.
(is)
t0
It is not possible, or practically impossible, to integrate analytically the equation (15). In applying (15) a second order method was introduced which turned out to be most effective for this particular problem. Based on FEM, a stress rate field T(o°) is chosen for the actual geometry at the time to with appropriate boundary conditions with an initial stress of To. For the stress rate field defined in such a way, a new stress state is indicated according to the following equation,
(:6)
TCa,) = To + ,~t'r~0 °).
For the stress state of T (A0 defined in such a way, but still with unchanged geometry and boundary conditions, the calculations of load increment are repeated to obtain a new stress rate TI(T(At}). On the basis of these two stress rate fields, it is possible to obtain the stress increment A T {A0 in accordance with the following equation
T" = To + A t ( ~ ( 0 °) + ~'r1(T(~'t))) = To + A T (A0,
(17)
wherein ~, ~, ~ are time integration parameters. The method described by (16) and (17) belongs to the class of second order Runge-Kutta methods with global truncation error of 0(AtS). For the error analysis, this method can be redefined to one free parameter ~, where -~ can be chosen, to optimize any desired feature of process. The global error becomes a minimum when the parameters are
1
3
83
NUMERICAL
STABILITY
To check the stability of the presented time integration method, calculations are repeated twice with a doubly reduced integration step, T at, t o + ~.L
T"=To+
f
tO+A t
H(T(r),D(r))dr+
f
H(T(r),D(r))d~'=To+AT
(2~).
(19)
to +
to
In this way, we are able to obtain a new stress state, T °'. Calculations are repeated until the following stability criteria are fulfilled
T`3 - 7, 3" 8(T0)=maxl~l ,: x / t r ( T 0~)
;
8=
8(To) 1 - 2.--~
• '
~>1.
(20)
The use of an estimate (20) defining the approximation error 6 with the method of order ~. The knowledge of this error enables an automatic choice of time increment according to the following procedure. In order to automatize the calculation on load increment level, two error limits el and e: (et < e2) must be introduced. If after calculations of the load increment - ~i < el then in the next step, the integration step At should be doubled, -
el < 6 < e2 then the time increment does not change,
-
6 > e2 then the calculation has to be repeated with a time step twice as small. In our analysis the following error limits were applied: el = 0.001, e2 = 0.01, ~ = 2. A NUMERICAL
EXAMPLE
Two examples ilustrating the presented numerical methods are shown below. The numerical tests of penetration are performed at vertical stress level tr~ = lOOkPa and at porosity n = 0.433
( n , ~ = 0.47 and n , ~ = 0.36), which is close to the critical density
of this sand. In the Figure 4 in the left diagram we present calculated penetrating resistance F, as a function of vertical displacement. Two different friction coefficients/J were supposed: for the curve No. 1 with pl = tan(10*) and the curve No. 2 with p2 = tan(27*). The comparison with experimental data from Smith [3] (right diagram of Figure 4) is shown in the left diagram of Figure 4 too. In Figure 5, the distribution of components of normal and tangential stress on the surface of penetrometer with friction ~ = tan(10*) is presented, s is a normalized current
84
0
5
I
20
15
....... fhlr~Jitl
0
iiii ........ Jll,l
-e
Ii"kl 11 I llt~ill i I I I I~kJ l
I~l.t
r i 501"1
t.o
c ...
I
,
I
-~0
"--]gxperiment: |b~ Sml~s
"60
i~'/r'[~'
r
~]
'\ x
-40
"q
I i I i I i 1~
;',,
--,,,
-30
~
~-U,'~,"l,i;Jiq
5
%
-20
"k,l l IN! ~N,
I
I II~I~I
~--'~ Compl
4
-t0
-2 I IiNlhlli'~r-k21 t ! l J l ~lllNJtl I Ill "s i[i ~II
I
\\
j I
Figure 4: Experimental and calculated penetration resistance.
"'; ~
"; - -
Incr. Load No. 1
Incr. Load No. 4
,./,, = ~.I[%]
3O
I
_.
• ,
Iro = 0.4[%]
J
10
I
L
'
. 2
i
i ' ; .-T
1
~
I
, I
I
:
I
i •
I
~
I
• •
!
i
I
I
Incr. Load No. 8
}0
• . I,•
J=0.0
I
i
4
= o.sl%l
I
•
[
I
"; ~
i
I
i
(
i
I
,
[
,
!
,
It
[
Incr. Load No. 18
13,
L I 2
I
I
1,0
I
J
2
I
.
i .4
J
t
1
1
.
I
,
t
i
Figure 5: Stress distribution during penetrating process for p = tan(10")
l •
J
I
i I
,
0
85 parameter of the skin surface: s = 0 at the lowest point of the penetrometer and s = 1 at the highest point of the penetrometr surface. From this figure, we can see the change of stress distribution and progressive violation of friction condition on the shaft of the penetrometer during the penetration process. In the case of axial symmetry, the lower point of the penetrometer can not be in a sliding state as the tangential stress at this point is equal to zero. The FEM calculations, suppossing large displacements, cannot be made without the possibility of changing the topology of mesh, as the deformation of boundary elements around the penetrometer foot for displacement greater than two times the penetrometer's diameter have an important influence on FEM errors, which is based on our numerical experience. CONCLUSIONS In the presented analysis we demonstrated the ngmerical posibility of directly describing a friction phenomenon on the boundary of the considered sand volume. This new method, which allows the utilization of the Coulomb friction condition as a boundary condition together with the principle of virtual work, shows that methods with interface elements to describe friction conditions on the boundary are not needed. The Coulomb friction condition was used only as an example. Instead of the friction condition, other relations as function either of the stress tensor or of its material derivative can be implemented in a FEM program. The penetration problem is non-linear on the load increment level and, this is why some methods of convergence checking must be introduced into the numerical analysis. In the case when divergence occurs a method for enlarging convergence region must be used. ACKNOWLEGMENTS The first author is grateful to Prof. Dr. H. Winter and Dr. habil. D. Kolymbas, IBF Technical University of Karlsruhe, for some extremely valuable discussions. The work reported in this paper was supported by the Volkswagenwerk Foundation under Grant Az: 1/63 374. We gratefully acknowledge this support. REFERENCES 1. Allersma, H.G.B., Optical analysis of stress and strain in photoelastic particle assemblies, Delft University of Technology, Delft 1987, p. 125 2. Kolymbas, D., Eine konstitutive Theorie ffir B6den und andere kSrnige Stoffe,
86 in german, Publication of Institut for Soil Mechanics and Rock Mechanics, Technical University Karlsruhe, Karlsruhe 1988 3. Smits, F.P., Cone penetration tests in dry sand, Proc. of the Second European Symposium on Penetration Testing, Amsterdam, May 1982, pp.877-881 4. Sikora, Z., On integration methods for rate independent constitutive laws, SIAM ,Journal Num. Anal. (to appear)
NUMERICAL SIMULATION OF PENETRATION IN SAND BASED ON FEM
Z. Sikora Institut for Soil and Rock Mechanics Technical University Karlsruhe West Germany and G. Gudehus Institut for Soil and Rock Mechanics Technical University Karlsruhe West Germany
ABSTRACT This article presents the application of the Finite Element Method (FEM) for estimating the penetration resistance and the stress and strain state in the soil around a penetrating rod. Numerical solutions are obtained assuming Lagrangian formulation with finite deformations. The Jauman terms are incorporated in t'he analysis. A constitutive law by Kolymbas is used. A condition representing skin friction has been defined for the surface of the penetrometer, and subsequently the unknown stress state on its surface is determined with an iterative method. The method is closely connected with the localization of deformation. Some numerical problems due to time-integration and convergence criteria for incrementally nonlinear constitutive laws are also discussed. In numerical computations, the friction coefficient parameter is varied. The influence of this parameter on the penetration resistance is discussed. INTRODUCTION The penetration test is widely used to investigate the mechanical properties of soil layers in the field. In this test, a probe with a conical tip is driven vertically into the soil, while the cone resistance and the total force are measured continuously. The cone resistance usually increases with depth. The measured cone resistance and shaft friction 73 Computers and Geotechnics0266-352X/90/S03.50 © 1990 Elsevier Science Publishers
Ltd, England. Printed in Great Britain
74 are used to estimate soil parameters, which are used in engineering practice. Much research has been done on the behaviour of soils during penetration. The information obtained from laboratory tests is limited because very little is known about the stress-strain state in the soil during penetration. Since stresses and strains are not uniformly distributed in this case, no techniques were available to estimate the stress distribution in the interior of the sample from the boundary conditions. The boundary condition for the surface of the penetrometer should be stated as an inequality describing the friction phenomenon. All techniques using interface-elements to describe the interaction between probe and soil are considered as subjective. Furthermore, it is dif~cult to define a constitutive relation for the soils; therefore introducing another constitutive law for the interface makes this problem even more obscure. The Coulomb friction condition can be used to describe an appropriate boundary condition. This paper presents new possibility of well-posedness of some boundary problems, where on a part of the boundary, a friction condition is valid. THE BOUNDARY PROBLEM FOR STATIC PENETRATION In the subsequent analysis, we consider a boundary value problem modelling a penetration process of a non-deformable penetrating rod into sand. Figure 1 shows the boundary conditions and a finite element mesh for the boundary problem which will be discussed here. Unless stated otherwise, upper case bold face letters are used to denote second order tensors while lower case bold face letters are used to denote vectors in R ~. To simplify the considerations, we have supposed that the penetrating rod is placed at a depth equal to the length of the penetrometer skin surface. An initial stress To is assumed to be known in the zone around the penetrometer. The lower part of the penetrating rod is taken as a hemisphere instead of a cone, which simplifies the calculations from a numerical point of view. Penetration experiments, carried out by Allersma [I], with different shapes of penetrometers showed that the shape of the penetrometer's tip has no major influence on the stress distribution in the granular material. The boundary condition on the surface of the penetrometer is defined as Coulomb friction, i.e. the magnitude of the tangential stress component must not be greater than the normal stress component multiplied b.v a friction coefficient/~. The friction coefficient/~ is assumed to be constant. The stress tensor T at any boundary point on the shaft of the penetrometer during the penetration process must fulfil the following condition on the surface of the penetrometer,
75
i" Boundary c o n d i t i o n s : AB :r.. = r . , = 0 B C :v~ = 0
and
Tr. = 0
CD:v,
--0
and
7",~ = 0
DE:v,
-- 0
and
7~,, = 0
EA : i f
I n T T t l < tJnTTn
v, = 0 else
and
then
v, = vm,
InTTtl =
pnTTn
b
-
¢
Figure 1: Boundary conditions and FE mesh for penetration problem.
InTTtl K pnTTn,
(1)
where vectors n and t are external normal and tangential vectors of the surface of the penetrometer respectively. T is the Cauchy stress tensor. The incremental equilibrlum will be described by the virtual work principle where time deviratives of Cauchy stress tensor, stretching, rate of surface traction and body force are components of the incremental work. The initial stress state To in the subsoil is the initial condition for the initial-boundary-value problem, and has to fulfil the condition (1). To determine the deformation of a body from the system of differential equations governing its mechanical deformations, we need side conditions such as boundary conditions. The remaining boundary conditions are those of rate of place vi(=,,,) (so-called Dirichlet boundary condition) and/or rate of traction (Neumann boundary condition), which are illustrated in Figure 1. To be able to solve the boundary problem with boundary condition in an inequality form, as defined in (1), a different numerical approach, in comparison to classical equilibrium problem with known boundary conditions, is needed. The main numerical problems due to this new approach will be presented in the sequel. The stress distribution is not known exactly on the boundary and is restricted by the friction condition. The Coulomb friction is defined with the assistance of the total stress, and not of its
76 increment. CONSTITUTIVE LAW The subsoil is assumed to be dry. The density of the sand is implemented in the parameters of the constitutive law. The tensor form of the constitutive law, which was proposed by Kolymbas [2] is as follows
"r = H ( T , D ) =
-• (TD+DT)+C2tr(TD)I+C3trTv/~+C4t--~ T2
(2)
is the rate of the corotational Cauchy stress tensor and Cl ..... C4 are known material constants evaluated on the basis of triaxial tests. 1 is the unit tensor. T and D denote Cauchy stress tensor and stretching. The material time derivative "~ of the Cauchy stress tensor T is assumed to be the Jaumann's type increment T minus T W - W T . i.e.
"~ = "[ + W T
- TW
(3)
where W denotes the spin tensor. There is a simple calibration method based on experimental data from a triaxial apparatus [2]. The results of this calibration can be used as input data to another numerical calibration method based on a random process within a nonlinear optimization program for obtaining a better approximation of the experimental data. The calibration process can be defined in general as a nonlinear optimization problem. The calibration is a very difficult problem from the mathematical viewpoint, specially for constitutive laws of the rate type. In this paper, the above mentioned calibration methods will be not presented because of the nature of the topic and size of this paper. As an example for this boundary value problem, we calibrated the four constitutive parameters C1,...,C4 to the triaxial experimental data with Oostershelde sand from Holland. The constitutive law was calibrated to the simple compression triaxial test with a confining pressure equal to 100 kPa. The obtained constitutive parameters are as follows
C, = -165.108,
C2 = 9.966,
C3 = 164.061,
C4 = -375.351.
(4)
The results of element test calculation are compared with triaxial compression test data as can be seen in Figure 2.
77
.8
..
.E O I~
.4
5.
/,
f
/o
5
•
s /
4-
4
Z
h
.<
3~ t~
4
t.J
4
W N
N
¢
Z C3
4
0.
n¢
-!
2
4 AXIAL
6
B
STRAIN
10
12
Ig)
Figure 2: Element test calculation in comparison to triaxial compression test M A I N P R O B L E M S OF N U M E R I C A L C A L C U L A T I O N S The presented constitutive law (2), which belongs to a large class of rate type constitutive laws, is very convenient from the numerical point of view. There is not, for instance, any artificial partition of deformation into elastic and plastic parts as with elasto-plastlc constitutive laws. The utilized constitutive law is a non-linear tensor function with resrect to the stretching D, i.e. the FEM equilibrium condition is a non-linear system of equations on the load increment level. In a Lagrangian incremental approach for large displacements we express the equilibrium of the body at a load increment number "i + 1" using the principle of virtual work as follows
K~,(x', Au'+~)~u~ ~ = gR~ ÷*
(5)
where K denotes a stiffness matrix (operator in the sense of FEM), x i is the actual configuration after load increment number ~i', ~ u I+I is the increment of displacement for the load increment A R I+1 within load increment number "i + 1". It is thought here that the problem becomes more complex as a non-linear system of equations must be solved for each load increment. The first problem consists of the choice of a particular numerical method for solving this non-linear system of equations. It is
78 clear that the amount of calculation for a non-linear system, for instance (5), is greater than in relation to linear ones (e.g. for load increment within an elasto-plastic constitutive relation). For the boundary problem of this paper a modified Newton-Raphson method was sucessfully used. Futhermore, the evaluation of the size of the convergence region of this non-linear problem is a very difficult problem. For the penetration problem it was necessary to introduce into the computer program a numerical method for enlarging the convergence region. Another fundamental difficulty in the general description of equations of equilibrium is that the configuration of the body for each load increment level "i + 1" is unknown. In our analysis, we consider the motion of the body in a fixed (stationary) Cartesian coordinate system where the configuration after each load increment is known, i.e.
x i+l = x i + A u i+1.
(6)
This problem with assumption of large deformation is treated as geometrically nonlinear, i.e. after Lagrange updating we obtain
j~+t
= Kt,,,(x,+,, & u , + t ) A u ~ , _ Kt.~(x', A u ' + ' ) A u ~ ' # 0.
(7)
where ji+l denotes a vector of correction forces which can be used as a load vector for the correction of solution accuracy in FE equilibrium equation as follows
K t , n ( x i , z ~ u i + l ) A u ~ 1 = AR~ +' - d}
and
do = O.
(s)
The presented method can be called a non-iterative correction method for geometric nonlinearity. As mentioned earlier, the constitutive relation is a non-linear tensor function of stretching, i.e. the boundary problem belongs to the group of physically non-linear problems. Introducing the Coulomb friction condition as boundary condition to the differential equations of equilibrium, which will be presented in the next section in detail, makes this problem aditionally non-linear, i.e. even in the case of a linear constitutive law and without assumption of large deformations, the boundary problem is formed in an inequality state (in the sense of boundary condition (1)) and remains therefore non-linear. Another problem of numerical solution concerns the time integration method. The solution of FE equilibrium state gives us an answer about distribution of the Cauchy stress
79 material derivative. For the penetration a second-order integration method (presentation of this method will be given in next section) is developed, which was very efficient specially for the problem considered. The last problem, which we will mention about numerical solution of penetration with utilizing the 'direct' boundary condition as friction condition, concerns of FE discretization. It is clear that because of large deformations around of the penetrating rod the mesh must be assumed finer as in the farther environment of the penetrometer, but very important for the convergence conditions of the non-linear equation system is that the discretizing nodes of penetrometer structure, which are also boundary points with friction condition, are placed in equal distances between them. A number of different calculations showed that a sufficient mesh width is about ten times the penetrometer radius, and the depth of the mesh can be equal to four times the penetrometer's length. FRICTIONAL BOUNDARY CONDITION A question arises: how does one take boundary friction conditions into account if instead of being defined 'explicitly' they are defined through their dependence on appropriate components of stress vector? As mentioned in the last section, the boundary problem with friction condition as boundary condition is a non-linear problem, and therefore the solving method of this problem is iterative. We are supposing the so-called kinematic load method, i.e. on the boundary upon which friction conditions are given are assumed as initial approximation of the friction condition the Dirichlet boundary conditions on the shaft of penetrometer. The penetrating rod is assumed as non-deformable, and the boundary conditions in the nodal points on the surface of penetrometer are as follows
v,=0
and
v,=v0.
(9)
The boundary conditions on the penetrometer shaft and on the other boundaries are fixed, and from the FE equilibrium condition (8) we can obtain velocity and stress rate fields. The current stress state and displacement field are obtained after time integration. The normal and tangential component of stress in nodal points on the surface of penetrometer caa now be calculated, and the condition (1) of friction can be checked (violated or not). In the case that the condition (1) is violated it means that the geometrical position of the point is wrong and must be corrected. Our task is to find such a configuration that in all boundary nodes the friction condition will be fulfilled. The
80
/':= t
=
(ex~a, rJ/,m)
~o := We = d
(,..,,,,,)
-
(0,,,,,)
=
~=,.',-,'+,,.o)
"_.,]_ __
Figure 3: Sliding conditions on penetrometer surface nodal point with violation of friction condition wiU be called sliding point or point in a sliding state. For the correct determination of geometrical position of sliding point we use, as additionaly equation of the FE non-linear equation system, the Coulomb friction condition in a form of equality. In Figure 3 we can see a detail of geometrical position of nodal point A0 which is assumed to be in a sliding state. The geometrical position of point A0 will be corrected to position A1 where the condition (1) is fulfilled with the equality symbol (Figure 3). We assume that the posible direction of nodal point motion on the shaft of penetrometer is tangential to the penetrometer surface. The geometrical relations are easily understood from the Figure 3. Thus the direction is defined in advance by the geometry of the penetrating rod structure. Boundary nodes themselves can only move along the surface of the penetrometer in accordance with the tangential vektor t at the skin surface. The proper correction vector wo, which starts at A0 and ends at the point A! for a particular nodal point, may be described as w0 = c . t with-c as a scalar. This scalar c is a further unknown in the node, which together with the unknown components of the rate vector Pl = (vrs, v,l) fulfils the following geometric equations
81 v,,
-
c. cosa
= 0
(10)
vz, + c - s i n a = v0.
(11)
The Coulomb friction equality is presented by the non-linear equation in relation to the stretching tensor, which via FEM reveals a non-linear vector rate function. Let us suppose, to simplify matters, that the new stress state, Tt, in the nodal point is defined a~
T,
= To + &ArT
(12)
where & is a time integration parameter, and At is the time step. The stress tensor TI may fulfil the following Coulomb friction equality, i.e.
InT(To + a A t T ) t ] = pnT(To + 6 A t T ) n
(13)
or in another form
- 1 7" n r H ( T 0 , D ) r = n r ' l ' r = ~-~-~n T0r where
;
r = t +
#sign(ro)n
(14)
r0 = nrTot.
The right hand side of the equation (14) depends on a known initial stress tensor To and time increment At. The left hand side of the equation (14) is the already mentioned non-linear function in relation to the rate field and becomes attached to the non-linear equation system resulting from the equilibrium conditions. For an arbitrary load increment we expect some problems with convergence because introducing the Coulomb friction equation into the global system of equations makes the condition number of the tangent stiffness matrix worse, and for this case, some methods for enlarging the convergence region must be applied. TIME INTEGRATION
METHOD
The method of time integration plays an important role in rate type constitutive laws using FEM. The solution using the appropriate load increment is obtained by integration of the load rate in relation to the time parameter. This time parameter in
82 relation to the time independent behaviour of the tested material is not perceived in its physical sense, but rather as a load parameter. Here we will present only the main idea of the integration method which lies in the class of methods with a global truncation error of 0(At3). The exact developement of this time integration method due to a rate type constitutive law with FEM will be presented in another paper [4]. The integration process takes place between to and to + 2xt, to obtain the new stress T" at time to + At, i.e.
t0+At
T" = T o +
/
H(T(r),D(r))dr.
(is)
t0
It is not possible, or practically impossible, to integrate analytically the equation (15). In applying (15) a second order method was introduced which turned out to be most effective for this particular problem. Based on FEM, a stress rate field T(o°) is chosen for the actual geometry at the time to with appropriate boundary conditions with an initial stress of To. For the stress rate field defined in such a way, a new stress state is indicated according to the following equation,
(:6)
TCa,) = To + ,~t'r~0 °).
For the stress state of T (A0 defined in such a way, but still with unchanged geometry and boundary conditions, the calculations of load increment are repeated to obtain a new stress rate TI(T(At}). On the basis of these two stress rate fields, it is possible to obtain the stress increment A T {A0 in accordance with the following equation
T" = To + A t ( ~ ( 0 °) + ~'r1(T(~'t))) = To + A T (A0,
(17)
wherein ~, ~, ~ are time integration parameters. The method described by (16) and (17) belongs to the class of second order Runge-Kutta methods with global truncation error of 0(AtS). For the error analysis, this method can be redefined to one free parameter ~, where -~ can be chosen, to optimize any desired feature of process. The global error becomes a minimum when the parameters are
1
3
83
NUMERICAL
STABILITY
To check the stability of the presented time integration method, calculations are repeated twice with a doubly reduced integration step, T at, t o + ~.L
T"=To+
f
tO+A t
H(T(r),D(r))dr+
f
H(T(r),D(r))d~'=To+AT
(2~).
(19)
to +
to
In this way, we are able to obtain a new stress state, T °'. Calculations are repeated until the following stability criteria are fulfilled
T`3 - 7, 3" 8(T0)=maxl~l ,: x / t r ( T 0~)
;
8=
8(To) 1 - 2.--~
• '
~>1.
(20)
The use of an estimate (20) defining the approximation error 6 with the method of order ~. The knowledge of this error enables an automatic choice of time increment according to the following procedure. In order to automatize the calculation on load increment level, two error limits el and e: (et < e2) must be introduced. If after calculations of the load increment - ~i < el then in the next step, the integration step At should be doubled, -
el < 6 < e2 then the time increment does not change,
-
6 > e2 then the calculation has to be repeated with a time step twice as small. In our analysis the following error limits were applied: el = 0.001, e2 = 0.01, ~ = 2. A NUMERICAL
EXAMPLE
Two examples ilustrating the presented numerical methods are shown below. The numerical tests of penetration are performed at vertical stress level tr~ = lOOkPa and at porosity n = 0.433
( n , ~ = 0.47 and n , ~ = 0.36), which is close to the critical density
of this sand. In the Figure 4 in the left diagram we present calculated penetrating resistance F, as a function of vertical displacement. Two different friction coefficients/J were supposed: for the curve No. 1 with pl = tan(10*) and the curve No. 2 with p2 = tan(27*). The comparison with experimental data from Smith [3] (right diagram of Figure 4) is shown in the left diagram of Figure 4 too. In Figure 5, the distribution of components of normal and tangential stress on the surface of penetrometer with friction ~ = tan(10*) is presented, s is a normalized current
84
0
5
I
20
15
....... fhlr~Jitl
0
iiii ........ Jll,l
-e
Ii"kl 11 I llt~ill i I I I I~kJ l
I~l.t
r i 501"1
t.o
c ...
I
,
I
-~0
"--]gxperiment: |b~ Sml~s
"60
i~'/r'[~'
r
~]
'\ x
-40
"q
I i I i I i 1~
;',,
--,,,
-30
~
~-U,'~,"l,i;Jiq
5
%
-20
"k,l l IN! ~N,
I
I II~I~I
~--'~ Compl
4
-t0
-2 I IiNlhlli'~r-k21 t ! l J l ~lllNJtl I Ill "s i[i ~II
I
\\
j I
Figure 4: Experimental and calculated penetration resistance.
"'; ~
"; - -
Incr. Load No. 1
Incr. Load No. 4
,./,, = ~.I[%]
3O
I
_.
• ,
Iro = 0.4[%]
J
10
I
L
'
. 2
i
i ' ; .-T
1
~
I
, I
I
:
I
i •
I
~
I
• •
!
i
I
I
Incr. Load No. 8
}0
• . I,•
J=0.0
I
i
4
= o.sl%l
I
•
[
I
"; ~
i
I
i
(
i
I
,
[
,
!
,
It
[
Incr. Load No. 18
13,
L I 2
I
I
1,0
I
J
2
I
.
i .4
J
t
1
1
.
I
,
t
i
Figure 5: Stress distribution during penetrating process for p = tan(10")
l •
J
I
i I
,
0
85 parameter of the skin surface: s = 0 at the lowest point of the penetrometer and s = 1 at the highest point of the penetrometr surface. From this figure, we can see the change of stress distribution and progressive violation of friction condition on the shaft of the penetrometer during the penetration process. In the case of axial symmetry, the lower point of the penetrometer can not be in a sliding state as the tangential stress at this point is equal to zero. The FEM calculations, suppossing large displacements, cannot be made without the possibility of changing the topology of mesh, as the deformation of boundary elements around the penetrometer foot for displacement greater than two times the penetrometer's diameter have an important influence on FEM errors, which is based on our numerical experience. CONCLUSIONS In the presented analysis we demonstrated the ngmerical posibility of directly describing a friction phenomenon on the boundary of the considered sand volume. This new method, which allows the utilization of the Coulomb friction condition as a boundary condition together with the principle of virtual work, shows that methods with interface elements to describe friction conditions on the boundary are not needed. The Coulomb friction condition was used only as an example. Instead of the friction condition, other relations as function either of the stress tensor or of its material derivative can be implemented in a FEM program. The penetration problem is non-linear on the load increment level and, this is why some methods of convergence checking must be introduced into the numerical analysis. In the case when divergence occurs a method for enlarging convergence region must be used. ACKNOWLEGMENTS The first author is grateful to Prof. Dr. H. Winter and Dr. habil. D. Kolymbas, IBF Technical University of Karlsruhe, for some extremely valuable discussions. The work reported in this paper was supported by the Volkswagenwerk Foundation under Grant Az: 1/63 374. We gratefully acknowledge this support. REFERENCES 1. Allersma, H.G.B., Optical analysis of stress and strain in photoelastic particle assemblies, Delft University of Technology, Delft 1987, p. 125 2. Kolymbas, D., Eine konstitutive Theorie ffir B6den und andere kSrnige Stoffe,
86 in german, Publication of Institut for Soil Mechanics and Rock Mechanics, Technical University Karlsruhe, Karlsruhe 1988 3. Smits, F.P., Cone penetration tests in dry sand, Proc. of the Second European Symposium on Penetration Testing, Amsterdam, May 1982, pp.877-881 4. Sikora, Z., On integration methods for rate independent constitutive laws, SIAM ,Journal Num. Anal. (to appear)