![]() |
![]() |
Viscoplastic Finite Element Model for Expansive Soils by Assistant Professor, Institut für Statik, Technical University of Braunschweig, Germany |
ABSTRACT
Finite element implementation of the constitutive law, which illustrate viscoplastic or time-dependent behaviour of expansive soils is described along with verification of the model with respect to field behaviour of structures built in and on expansive soils. The model takes into account the current stress, water content and clay content as well as environmental factors. Procedures to derive model parameters are presented together with validations of the model. The model is verified with respect to the observed behaviour of expansive clay samples found in literature.
KEYWORDS: Finite element method, expansive soils, clay, swelling
INTRODUCTION
The finite element analysis of soil structures has been the object of many studies in the past and much progress has been made since the first publications appeared. The finite element method can help to provide a clearer understanding of soil behaviour, even in cases where soil data are limited, and can thus assist in making decisions on design. The recent increase in development and application of many computer-based techniques such as finite element, finite difference, and boundary element methods allow to apply complex constitutive models as those covering swelling. As the finite element method is well described and reported in literature, only a brief presentation of the used program is presented
OBJECTIVES
One of the aims of this study is to develop a finite element algorithm following the variational formulation of the elasto-viscoplastic problems of expansive soils. The purpose of the analysis is to evaluate the displacement pattern (heave and shrinkage values)within the expansive soil areas. This may help in the identification of the behaviour of the swelling soils and may produce data which can be used for the design of structures built on and in such type of soils. The time-dependent deformation and stress changes are not associated with pore-water migration only, but also with the swelling and viscous nature of the considered material. Here emphasis is put on the behaviour of the swelling soils rather than that of its effects on structures.
THE FINITE ELEMENT PROGRAM
A non-linear elasto-viscoplastic finite element computer program is developed in this study. This program has the capability to treat the problem as elastic, plastic, and viscoplastic each separately or in any required combination of these three material phenomenon. As the swelling behaviour of soils is a time dependent one, a viscoplastic theory is considered as it allows the modelling of time rate effects in the plastic deformation process. Consequently after initial yielding of the material the plastic flow, and the resulting stresses and strains, are time dependent. The program can apply isoparametric quadratic 8- as well as 9-Nodes elements. The tangential stiffness method and other two combined algorithm are used for the calculation of the stiffness matrix. All final results are achieved by applying the tangential stiffness method. Von Mises criterion, Drucker-Prager criterion and a modified one from Drucker-Prager discussed in the following sections are enhanced in this program. The program can generate any shape of FE-mesh.
CONSTITUTIVE LAWS FOR EXPANSIVE SOILS
One of the main objectives of this work is to produce a unified model for the swelling and shrinkage behaviour of expansive soil. Expansive soils are those clayey materials which exhibit significant volume changes caused by changes in the subsoil moisture. The expansive soil swells if its moisture content increases and it shrinks if its moisture content decreases. This phenomenon depends mainly on the mineralogical content of the clay soils. The montmorillonite clay minerals are considered as highly expansive and the most effective ones for swelling behaviour. This type of soil is found in many places all over the world especially in the arid and semi-arid regions. Many damages to structures constructed in and on expansive soils in different places are reported in literature [1]. The early papers report on structural damages caused by the soil movements and on considered measures for preventing or reducing these damages. Then studies followed dealing with the identification and classification of expansive soils and with the study of their properties and the associated soil surface heave. A historical collection of those studies with the relative discussion are available in literature [2], [3]. The proper constitutive models for swelling soils have to consider factors such as elastic, plastic, and viscous(time-dependent) deformations, volume change, and water diffusion in a unified manner.
The magnitudes of swell percent and swelling pressure are influenced by a large number of factors, which can broadly be divided into the following:
i. Compositional factors, which include the type and amount of clay mineral present in the soil as well as the pore fluid characteristics.
ii. Environmental factors, which include initial moisture content, initial density, initial degree of saturation, initial soil structure, stress history, availability and composition of ambient water and temperature.
iii. Procedural factors in laboratory testing e.g. size and shape of soil sample, degree of disturbance and testing procedure and techniques used.
In addition to the above factors, the field swelling behavior of expansive soils is largely governed by several other factors e.g.
iv. climate,
v. depth of active zone,
vi. location and thickness of the expansive soil layer,
vii. applied loads(weight of structure and soil overburden),
viii. vegetation,
ix. site topography and surface drainage, and
x. confinement.
Using the above general assumptions the general form of the constitutive equation considering the effect of current stress, water content and clay content on the performance of the rate of swelling strain has been introduced as follows
|
(1) |
where
is the strain rate
is the current stress
A, A´, B, a and b are constants
w is the water content as percentage.
The constant A in the above constitutive model represents the climatic site topography and confinement and is constant for a certain expansive soil zone. Also the constant b ¢ is the same per expansive soil zone and takes into account the consideration of the weight of structure and soil overburden.
The effects of the initial moisture content, initial degree of saturation and availability and composition of ambient water and temperature on the constitutive model are covered by the constants B, a, b. The constant A´ which is different from layer to layer under the same climatic factor considers the compositional and environmental factors discussed above.
k: constant value for each layer depends on the clay content and mineralogical constant.
The constitutive law incorporates the visco-plastic behaviour of such soil when its moisture content increases or decreases. This material model is calibrated to experimental tests found in literature. The model is validated in this study with respect to viscoplastic behaviour observed in laboratory tests available in literature.
IMPLEMENTATION OF THE MODEL INTO A FINITE ELEMENT PROGRAM
In this section a general description is presented of the basic expressions of the elasto-viscoplastic media of a solid continua and the detailed description of the developed constitutive law.
General statement of Constitutive Model
The total strain rate is divided quite generally into
elastic
viscoplastic
and initial
, such as that due to temperature[4]
|
(2) |
in which the overdot denotes the rate of change with respect to time.
Useful constitutive models for complex geological materials are obtained by using the incremental forms. Then the relation between increments of stress and strain can be expressed as follows:
|
(3) |
where D-1 is a symmetric elasticity matrix (compliance matrix).
The onset of viscoplastic behavior is governed by a scalar yield condition of the form
|
(4) |
in which y is the yield stress and k is a history dependent hardening (or softening) parameter. It is assumed that viscoplastic flow occurs for values of F>0 only. The general form of the viscoplastic flow rule given by Zienkiewicz and Cormeau [4] is used in this study.
|
(5) |
in which Q = Q(s evp, k) is a plastic potential and g is a fluidity parameter controlling the plastic flow rate which may depend on some state variables--such as time, total strain invariants, etc. The term F(F/Fo) is a positive monotonic increasing function for F>0 and the notation < > implies
|
(6) |
Fo denotes here any convenient reference value of F to render the expressions non-dimensional.
Equation (7) shows an analogy with the flow rule of conventional non-associated plasticity:
|
(7) |
where dl is a proportionality constant termed the plastic multiplier. If the associated plasticity situations, in which Fº Q, is considered, the general equation of viscoplastic flow rule reduces to
|
(8) |
The function F is taken to be equal to the viscoplastic equation (1) which depends on the current stress and water content.
|
(9) |
Since the effect of clay content in the constitutive equation is constant for each assumed homogeneous layer, it is considered to be included in the fluidity parameter g
|
(10) |
The constants A, A¢, B, b, b´ in the above equations depend on the compositional and environmental factors.
The term ¶ F/ s is the flow vector a and defined as follows
|
(11) |
Computation of the Strain and Stress Increments
The finite element program employed in this study uses three time stepping schemes,
The three methods predict approximately the same steady state values. The semi-implicit and explicit methods need smaller efforts of computation than the implicit method and producing rapid results for the steady state condition in respect of the transient stage. For this reason they are used for primary computation to test the suggested parameters as well as to have an idea over the steady state values and their time. The final computation for every problem in this study is done by using the implicit method.
In the following section the strain and stress increments are computed according to the three schemes individually.
Implicit method
The strain increment Dn, occurring in a time
interval
(12) A fully implicit (or backward difference) scheme is obtained by considering Q=1. The strain increment in this case is determined from the strain rate corresponding to the end of the time interval. Taylor series expansion is used to define the viscoplastic strain rate at time n+1, (13) where Dsn is the stress change occurring in the time interval
Dtn=tn+1 - tn and (14) The matrix H is the matrix whose eigenvalues determine the limiting time step length, Dt which can be employed in the explicit integration schemes and is depending on the stress level. Matrix Hn must be explicitly determined for the yield criterion assumed for the material behavior. From equations (8) and (14) (15) where the symbols á F and the superscript n are dropped for convenience. The function F in this study was taken to be dependent on the first invariant only, (16) where the flow vector is: (17)
(18) (19) For a three dimensional situation (20) or from equation (19) (21) (22) Substituting into equation (15) and restoring the symbols< > (23) The form of d f / dF depends on the explicit form of F discussed in the previous sections. For two dimensional situations (plane stress, plane strain and axial symmetry) the M matrix reduces to
(24) For plane stress and plane strain problems only the upper 3X3 partition is employed while for axisymmetric situations the complete matrix is utilized. Thus, equation (12) for a purely implicit time stepping scheme can be written as (25) (26) where (27) The incremental form of the stress strain relation presented in equation (3) is: The total strain increment in terms of the displacement increment can be written as (28) Since the associative case is considered in this study, the matrix Dn is a symmetric matrix. (29) Substituting for Devpn, from (25), then (27) becomes: (30) Semi implicit method The case Q = 1/2 results in the implicit trapezoidal scheme which is also known as the Crank-Nicolson rule in the context of linear equations. The semi implicit method follows the same derivation process of the equation of the fully implicit method described above, except that equation (26) becomes (31) Explicit method The Euler time integration scheme which is also referred to as fully explicit is obtained by considering Q = 0. The explicit method is also called forward difference method, since the strain increment is completely determined by conditions existing at time tn. For the solution of linear elastic problems by the explicit scheme (Q = 0), equation (29) simplifies considerably to (32) Equations of Equilibrium The incremental form of the equations of equilibrium, which should be satisfied at any instant of time tn, are (33) where the vector Dƒn represents the change in the vector of equivalent nodal loads due to applied surface tractions, body forces, thermal loads, etc, during the time interval Dtn. Since in this study the load increments are applied as discrete steps, the change of loads for all time steps other than the first step within an increment is equal to zero (34) Considering equations (25) and (32), the displacement increments occurring during time step
Dtn can be calculated from (35) where KTn is the tangential stiffness matrix of the following form (36) and DVn are termed the incremental pseudo-loads and are given by (37) Substituting the values of the displacement increments in (32) yields the stress increments Dsn. The total stresses and displacements are given by (38) From equations (27) and (28) (39) and then (40) The stress increment calculation based on the linearized form of the incremental equilibrium equations (33), the total stresses sn+1, obtained by accumulating all such stress increments are not strictly correct. They will not exactly satisfy the equations of equilibrium (41) The equilibrium corrections of this program are to evaluate sn+1 according to (32) and (33) and then compute the residual, or out-of-balance, forces Y as (42) The computed residual force is then added to the applied force increment at the next time step. Magnitude of Time Step for Numerical Stability The time step length adopted in this study is basically empirical. The value of the time step length can be constant through all time intervals or varied for each time interval. In the variable scheme the
magnitude of the time step is controlled by a factor τ which limits the maximum effective viscoplastic strain increment (43) The value of the time step length is taken as the least value of
Dtn satisfying the above equation at the Gaussian integration points. A variant on the above is to limit the time step length according to (44) in which (45) The value of the time increment parameter τ must be specified by the user. For
explicit time schemes accurate results have been obtained[4], [5] in the range 0.01<
τ <0.15. For implicit schemes, values of τ up to 10 have been found to be of stable
accuracy. A further limit is generally imposed relating the (46) change of time step between successive intervals i.e. where k is a specified constant taken to be normally 1.5. This value is suggested from
experience. Convergence Criteria It is essential to include a reliable convergence criteria in the program in order to identify correctly the increment in the viscoplastic strain. The attainment of steady state conditions can be monitored by accumulating some measures of the viscoplastic strain rate for all Gauss points in the structure. The degree of total viscoplastic flow at any point is best monitored by evaluating the total
effective viscoplastic strain rate at all Gauss points according to (47) Considering a global measure of convergence the steady state condition is assumed to be achieved when (48) in which the summation is taken over all Gauss points and TOLER is a convergence tolerance value given as input data. In this program the load is applied in discrete increments. An increment of load is applied and the time stepping process is followed until either steady state conditions are achieved, or a specified number of time steps is reached. Then a further increment is applied and the process repeated for all load increments.
as a fraction of the total effective strain
, so that
is the first total strain invariant and
is the first viscoplastic strain rate invariant. Thus Dtn can be formally written for this case as
Verification of the developed Constitutive Model
The constitutive model is verified in this work by comparing its performance with experimental data given in literature. Two comparison between the computed data and experimental ones by Madsen and by Mazurik and Komornik are discussed.
Mazurik and Komornik Experiment
The developed constitutive law is verified by comparing its performance with the experimental data given by Mazurik and Komornik[6]. The experimental data produced from many tests performed on almost identical clay samples(a=45.34 sq.cm, Gs=2.71, Wl=80%, Wp=24%), compacted into rings at a controlled initial moisture content of 22%, γd=1.480 ±2% gr/cu.cm, an initial height of 1.840 ±2% cm, and an initial degree of saturation S=72%.
where
a: area,
Gs: specific gravity
Wl: liquid limit
Wp: plastic limit
γd: dry unit weight
S: degree of saturation.
These tests were carried out on unsaturated clay samples subjected to inundation under various initial non-restrained constant loadings, with the intention of determining the log time-percentage swell and pressure-swell percentage curves.
The percentage swell at any time for different stress level can be determined by integrating the equation for the slope of percentage swelling with respect to time.
|
(49) |
The constants A, b and c for the above equation are computed by considering the data produced from the above experiment. The curves produced by the equation and the experimental data are plotted together in Figure 1.
It can be seen from Figure 1 that the calculated values of the percentage swell at different times are in a good agreement with the experimental ones.
Madsen Experiment
In this section the results of swell due to the expansive phenomena of the soil obtained from the finite element program are compared to the experimental ones from Madsen[7] on expansive soils. Table 1 shows the protocol of the swelling pressure experiments. In this experiment the sample of soil is loaded with a normal stress of value 5 kN/m².
![]() | |
Table 1. Protocol of the Madsen experiment | |
![]() | |
Drilling number: | Station 546 |
Orientation: | vertical downwards |
Depth of sample: | 2.65 m |
Time of sampling: | 25.10.1983 |
Sample location: | 18°c |
Sample location time | 25.10.1983-8.11.1983 |
Dimentions of Specimen: | d=70 mm, l=20 mm |
Density of the sample: | 2.45 g/cm³ |
Water content: | 3.3% |
Test climate: | 15°c, 80% Humidity. |
![]() |
Figure 1. Computed time - % swell curves of swelling clay under various constant loads in comparison with experimental ones.
Figure 2. Comparison of finite element results and experimental results by Madsen.
Figure 2 presents the experimental and the computed results of the considered sample of soil. The computed results seem to be in a good agreement with the experimental ones. The maximum swelling strain recorded by the experiment after 35 days is equal to 7.02% while the theoretical swell strain at the same time is equal to 7.33%. This small value of error (%4) verifies the ability of the program to f(6) predict the swelling effects.
CONCLUSION
A constitutive model is derived to represent the viscoplastic behavior of the expansive soil upon wetting and drying. The constitutive model is derived as a product function of three parts namely the stress, moisture content, and clay content.
The performance of this material model agrees with the experimental results.
The implementation of this constitutive model into the finite element algorithm is tested and verified by applying the model to various types of loading and boundary conditions. The verifications ensure its validity and capability of tackling different types of problems of swelling and contraction as well as the effects on the structures built in or on such types of soil.
REFERENCES
Jones, D. E., and Holtz, W. G.,: Expansive soils-the hidden disaster. Civil Eng., ASCE 43(8),pp 49-51,1973.
Rao, R.; Smart, P.,: Significance of Particle Size Distribution in Prediction of Swell Properties. Proc. of 4th Int. Con. on Expansive Soils, Vol. 1, pp. 96-105, 1980.
Snethen, D. R.; Huang, G.: Evaluation of Soil Suction-Heave Prediction Methods, Proc. of 7th Int. Conf. on Expansive soils, Dallas, Texas, Vol 1, pp. 12-17, 1992.
Zienkiewicz, O. C.; Cormeau, I. C.: Visco-plasticity and creep in elastic solids- a unified numerical solution approach", International Jour. of Numerical Methods in Engineering, 8, 821-845(1974)
Dinis, L. M. S.; Owen , D. R. J.: Elastic-viscoplastic analysis of plates by the finite element method, Computers & Structures, 8, 207-215 (1978).
Mazurik, A.; Komornik, A.: Interaction of Superstructure and Swelling Clay, Proceedings of the Third International Conference on Expansive Soils, Vol. I, Divisions 1-5, Haifa, July, 1973, pp 309-317
Madsen F. T.: Quelldruckmessung an Tongesteinen und Berechnung des Quelldruck nach der DLVO - Theorie, Zürich, 1976
![]() | |
© 2002 ejge |