ELASTIC/VISCOPLASTIC MODELS

FOR GEOMATERIALS AND POWDER-LIKE MATERIALS

WITH APPLICATIONS

By

JISHAN JIN

A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL

OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT

OF THE REQUIREMENTS FOR THE DEGREE OF

DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA

To My Family,

In Memory of My Father

ACKNOWLEDGEMENTS

I would like to express my deepest gratitude to my supervisory committee chairman,

Graduate Research Professor Nicolae D. Cristescu, for his guidance, tremendous support, and

encouragement during the study. I also sincerely appreciate Professor Cristescu's patience,

kindness, and discussions to his students.

I am deeply indebted to Professor Renwei Mei for helpful discussions on the

numerical analysis aspects during the work and for the critical reading of the early draft of

this dissertation. I am very grateful to Professor Frank C. Townsend for his guidance and

helpful discussions on the experimental aspects. I would like to thank Professors Bhavani V.

Sankar and Edward K. Walsh for their support and for serving on the committee.

Many thanks also go to my wife, Qianhong, for her understanding and full support,

and to my family members for their endless support in my life.

Finally, the financial support from the Engineering Research Center at University of

Florida, Particle Science and Technology, is gratefully acknowledged.

TABLE OF CONTENTS

ACKNOWLEDGEMENT ................ ............................ iii

ABSTRACT ............................................ ............ vii

CHAPTERS

1 INTRODUCTION ............................................... 1

1.1 Background .................................... .......... 1

1.1.1 Rock Salt ......................................... 2

1.1.2 Powder M materials .................................... 4

1.1.3 Models ............................................ 6

1.2 Literature Survey ............................................. 8

1.2.1 Mohr-Coulomb Model and Drucker-Prager Model ............ 8

1.2.2 Cap Model .............. ............................ 10

1.2.3 Cam-Clay M odel ...................................... 13

1.2.4 Lade Model ........................................ 16

1.2.5 Endochronic Model ................................... 17

1.2.6 Models for Compaction (Consolidation) of Powders .......... 18

1.2.7 The Elastic/Viscoplastic Models .......................... 21

1.2.8 Other M odels ........................................ 24

1.2.9 Remark on Recent Development of Plasticity-Based Models ... 25

1.2.10 Finite Element Method ........................... 25

1.3 Objective and Outline of the Present Work ......................... 27

2 THE ELASTIC/VISCOPLASTIC (CRISTESCU) THEORY ................ 31

2.1 The Elastic/Viscoplastic Constitutive Equation ...................... 31

2.2 Determination of the Elastic/Viscoplastic Constitutive Equation ........ 34

2.2.1 Elastic Parameters and Viscosity Coefficient ................ 34

2.2.2 Yield Function .................................. ..... 37

2.2.3 Viscoplastic Potential ................................ 38

2.3 Loading Unloading Conditions ............................ .. 39

3 AN ELASTIC/VISCOPLASTIC MODEL FOR TRANSIENT CREEP OF

ROCK SALT ................... ........... .............. 40

3.1 Determination of the ElasticViscoplastic Constitutive Model .......... 41

3.1.1 Determination of the Elastic Parameters .................... 41

3.1.2 Determination of the Compressibility/Dilatancy Boundary

and the Failure Condition ............................ 43

3.1.3 Determination of the Yield Function ...................... 45

3.1.4 Determination of the Viscoplastic Potential ................. 50

3.2 Comparison with the Data ................................... .. 58

3.3 The Finite Element Analysis .................................... 63

3.3.1 Formulation of the Elastic/iscoplastic Theory .............. 63

3.3.2 Example I: Axial Compression with Confining Pressure ....... 65

3.3.3 Example II. Stress Analysis for A Cylindrical Cavity in Rock

Salt ...................... .. .................... 67

3.4 Discussion and Conclusion ........................... ....... .. 72

4 TRIAXIAL EXPERIMENTAL RESULTS ON ALUMINA POWDERS ....... 74

4.1 Introduction ....................... .................. 75

4.2 Experimental Procedures ............................ ...... 78

4.2.1 Triaxial Equipment Setup ............................ 78

4.2.2 Specimen Preparation ................................ 83

4.2.3 Back Pressure Saturation ........................... 84

4.2.4 Hydrostatic Loading Test ............................ 85

4.2.5 Deviatoric Loading Test .............................. 86

4.3 Experimental Results .. ................... ................ 88

4.3.1 Hydrostatic Tests .......................... ... 88

4.3.2 Deviatoric Tests ...................................... 98

4.4 Discussion and Conclusion ................................ 111

5 THE ELASTIC/ISCOPLASTIC MODEL FOR ALUMINA POWDER A10.. 116

5.1 Determination of the ElasticViscoplastic Model ................... 116

5.1.1 Elastic Parameters .................................... 117

5.1.2 Compressibility/Dilatancy (C/D) Boundary and Failure

Condition ................... .... ............. 119

5.1.3 Yield Function .................................... 121

5.1.4 Viscoplastic Potential ............................... 126

5.2 Comparison with Experimental Data ............................. 135

5.2.1 Triaxial Tests .................. ..................... 135

5.2.2 Hydrostatic Tests ........................... ... 141

5.2.3 Rate Dependent Tests ............................ 144

5.3 Discussion and Conclusion .................................. 147

6 A NEW METHOD TO DETERMINE THE ELASTIC/VISCOPLASTIC

MODEL ............... ................................. 149

6.1 Introduction ............. ................. .............. 149

6.2 Elastic/Viscoplastic Constitutive Model .......................... 151

6.2.1 Elastic/Viscoplastic Constitutive Equations ................ 151

6.2.2 Loading-Unloading Conditions .......................... 153

6.3 Determination of the Elastic/Viscoplastic Model for Alumina Powder

A16-SG ............. .............................. 155

6.3.1 Elastic Parameters .................................. 155

6.3.2 Yield Surface ..................................... 155

6.3.3 Irreversible Strain Rate Orientation Tensor ................ 162

6.3.4 Viscosity Coefficient ............................... 167

6.3.5 Remark ........................................ 168

6.4 Validation of the Model ..................................... 169

6.4.1 Creep Type Formula with Stepwise Stress Variations ........ 169

6.4.2 Conventional Triaxial Tests ............................ 170

6.4.3 Constant Mean Stress Tests ............................ 173

6.4.4 Hydrostatic Creep Tests .............................. 175

6.5 Discussion and Conclusion .................................. 175

7 CONCLUSION AND FUTURE WORK .............................. 179

7.1 Conclusion ................... ........ ..................... 179

7.2 Future W ork ............. .... ....................... ..... 180

APPENDICES

A THE FORMULATION OF ELASTIC/VISCOPLASTICITY FOR THE

FINITE ELEMENT METHOD .......................................... 181

B TRIAXIAL EXPERIMENTAL RESULTS FOR ALUMINA POWDER A10

(DENSE) .................. ............... ............ 191

C TRIAXIAL EXPERIMENTAL RESULTS FOR ALUMINA POWDER A10

(LOOSE) ................... ............... ........ 228

D TRIAXIAL EXPERIMENTAL RESULTS FOR ALUMINA POWDER A16-

SG ......................................... ............ 248

REFERENCES ................................................. 270

BIOGRAPHICAL SKETCH .......... ......................... 281

Abstract of Dissertation Presented to the Graduate School

of the University of Florida in Partial Fulfillment of the

Requirements for the Degree of Doctor of Philosophy

ELASTIC/VISCOPLASTIC MODELS

FOR GEOMATERIALS AND POWDER-LIKE MATERIALS

WITH APPLICATIONS

By

JISHAN JIN

December 1996

Chairperson: Graduate Research Professor Nicolae D. Cristescu

Major Department: Aerospace Engineering, Mechanics and Engineering Science

An elastic/viscoplastic model with a nonassociated flow rule for transient creep of

rock salt is formulated based on a set of triaxial tests. A new procedure to determine the yield

surfaces and the potential surfaces is proposed. The asymptotic behavior and the physical

grounds of yield surfaces and of potential surfaces are incorporated into the model so that

these surfaces are smooth and the model better matches the data. This type of model with

both yield surfaces and potential surfaces is implemented for the first time into a finite

element program. It is shown that the prediction of the model matches quite well the

experimental data.

The mechanical behavior of alumina powders is investigated experimentally. Three

series of triaxial tests are performed. It is shown that the mechanical (compaction) behavior

of alumina powders is strongly dependent on the particle size and on the initial density. For

the powder A16-SG with smaller particle size (0.4-lpm), the volumetric strain always

exhibits compressibility, while for the powder A10 with larger particle size (40-200pm), the

volumetric strain exhibits both compressibility and dilatancy.

An elastic/viscoplastic model with a nonassociated flow rule to describe the nonlinear

time-dependent behavior of alumina powder A10 is formulated based on five conventional

triaxial tests. By using such type of tests, a new procedure to determine the yield surfaces and

the potential surfaces is proposed. The model is checked against the data. A good agreement

between the model prediction and the data is obtained.

A new methodology is proposed to determine the elastic/viscoplastic model for

compressible materials. In such models, the irreversible volumetric strain is chosen as a work

hardening parameter. A model to describe the behavior of alumina powder A16-SG is

formulated. The irreversible strain rate orientation tensor is determined. The model is easier

to handle and has fewer parameters than the previous ones. A very good agreement between

the data and the model prediction is obtained.

CHAPTER 1

INTRODUCTION

1.1 Background

The formulation of general constitutive equations of three-dimensional models for

time dependent materials has been studied since the pioneering works by Maxwell in 1868

and by Kelvin in 1875 (see Malvem [1969]). Many materials possess time dependent

properties: for instance, polymers exhibit viscoelastic properties; a rock flows in a long

period of time although it seems strong and hard; a metal creeps at high temperature, etc. A

great number of contributions have been made in this field (see Perzyna [1966]; Cristescu

[1967]; Christensen [1982]; Cristescu & Suliciu [1982]; Phillips [1986]; Cristescu [1989];

Chaboche [1989]; Lubliner [1990]). Recently three-dimensional modeling for rock salt and

powder materials has attracted the attention of many researchers, due to the demand of

various engineering industries (see Hardy & Langer [1984, 1988]; Chenot et al. [1990]; Jinka

& Lewis [1994]; Krishnaswami & Trasorras [1995]; Xu & McMeeking [1995]; Aubertin

[1996]; Tszeng & Wu [1996]). They wish to develop general constitutive equations or

modify the existing ones which, on the one hand, can be used to describe better the behavior

of materials, and, on the other hand, can be incorporated into the finite element method.

Further, the models will be used for optimal designs or/and safety evaluations.

1.1.1 Rock Salt

Rock salt exhibits time dependent properties as revealed by creep, relaxation and rate

dependent tests (Hansen & Carter [1984]; Handin et al. [1986]). It is found under natural

underground conditions. Rock salt consists mainly of halite with little anhydrides and clay,

and has small permeability. Due to the theological behavior and the physical properties of

rock salt, the large underground cavities in mass of rock salt are considered as ideal places

to store radioactive waste or hazardous chemical waste. Many salt cavities were excavated

around the world for this purpose (See Hugout [1988] in France; Eekelen [1988] in

Netherlands; Matalucci & Munson [1988] in USA; Kappei & Gessler [1988] in Germany).

The deformations of those cavities are continuously recorded. The stability and the safety of

those large underground cavities are being evaluated. Since the stability and the safety are

not precisely predicted nowadays, millions of extra dollars are spent on maintaining those

cavities in a workable state. Thus it is highly necessary to formulate a general constitutive

equation in order to predict the mechanical behavior of rock salt and to evaluate the safety

of cavities. Also a general constitutive equation would be beneficial, not only for the

evaluation of the stability and the safety of large underground cavities, but also for other rock

salt mining engineering problems such as the improvement of salt mining design procedures

(Rolnik [1988]).

Most of the research done so far was devoted to the long term behavior of rock salt,

i.e., stationary creep. A variety of models have been proposed (e.g., Langer [1988]; Chan et

al. [1994]; Aubertin et al. [1996]). Most of them can be used to predict the long term

behavior of rock salt by neglecting the deformation in the transient creep phase. However,

occasionally, a great discrepancy is found between the measured creep deformation of

underground cavities or openings and that predicted by the models (Aubertin et al. [1993]).

It is believed that this discrepancy is largely due to the neglect of the transient creep

deformation, because the stress around an opening (or cavity) is redistributed during the

excavation of the opening in the transient creep phase. The redistributed stresses will

influence the stationary creep of rock salt around the opening. Therefore it is crucial to

formulate correctly the constitutive equations for transient creep, if one wishes to describe

the behavior of rock salt in the rock surrounding the excavation during the stationary creep.

Several models for transient creep of rock salt have been formulated (e.g., Langer

[1984]; Munson & Dawson [1984]; Desai & Varadarajan [1987]; Aubertin et al. [1993]). In

most cases one starts from empirical formulae using a viscoelastic or viscoplastic model with

Drucker-Prager type of yield surface or internal state variable concepts. However, these

models either do not give good matching with the experimental data or could hardly be

incorporated into a finite element program. Thus these models should be judged from the

point of view of how well they match the experimental data and how easily they can be

incorporated into a finite element program. In the present thesis, we start from a set of true

triaxial tests on rock salt and follow a procedure to formulate a general three-dimensional

elastic/viscoplastic constitutive equation for the transient creep phase. No empirical formulae

are used since empirical formulae may not be adequate for other cases. Nonassociated flow

rule is explored and the viscoplastic potential is determined based on the experimental data.

The model not only matches better the experimental results, but can also be incorporated

easily into the finite element method. The model will be presented in Chapter 3.

1.1.2 Powder Materials

The mechanical behavior of powder materials has been studied for a long time, since

powder materials are used in many industrial applications such as metal powder metallurgy

(Johnson [1992a, 1992b]), ceramics injection molding (Mutsuddy & Ford [1995]) or food

products (Puri et al. [1995]). Powder materials such as metal powder or ceramics powder,

by some compaction or/and sintering processes, become useful products with very dense and

hard properties. However, in shape forming processes such as uniaxial compressing,

ceramics injection molding or pressure casting, the density distribution in the product is

nonuniform due to the friction between powder and the wall of a die. The nonuniform

density distribution would cause warping or other defects in the product in a later sintering

process. These warpings or defects will damage the products and unqualify the products

(Reed [1988]). Nowadays, most shape forming techniques are the ad hoc ones. Due to the

lack of the correct prediction of density distribution in a mold, many products are not

qualified, and millions of dollars are wasted each year (Johnson [1992b]). Thus it is crucial

to understand the mechanical behavior of powder materials during these processes and it is

urgent for designers to simulate these shape forming processes and to predict the density

distribution in a mold.

In order to predict the density distribution of powder in a mold, and to simulate the

shape forming process, the three-dimensional stress strain relationships, i.e., the constitutive

equations for powders have to be correctly and precisely formulated. Triaxial tests are often

used to formulate constitutive equations (Desai & Siriwardane [1984]; Cristescu [1989]).

Unfortunately only very limited triaxial tests have been reported (Koerner & Quirus [1971];

5

Shima & Mimura [1986]; Brown & Abou-Chedid [1994]; Gurson & Yuan [1995]) for

different applications. Without complete experimental data, it is difficult to formulate

constitutive equations correctly. In order to understand the compaction behavior of powder

materials, to formulate a general three-dimensional model, or/and to examine the existing

models, additional triaxial tests are needed.

Alumina powder (AzO03) has been widely used in ceramic engineering for various

applications due to the specific properties of alumina such as high electrical resistivity, high

erosion resistance, high melting temperature and high abrasiveness etc. For instance, alumina

powder can be used in the manufacture of porcelain, alumina laboratory ware, wear-resistant

parts, sandblast nozzles, medical components, abrasives and refractories etc. Hundreds of

tons of alumina powder and alumina-based articles are produced each year (Richerson

[1982]). However, very few experiments have been performed for alumina powders

regarding to the mechanical (compaction) behavior in three-dimensional stress states. Only

some tests in uniaxial conditions have been reported for ceramics powders (e.g., Stanley-

Wood [1988]; Gethin et al. [1994]; Chen et al. [1994]). Unfortunately, uniaxial tests reveal

mainly one-dimensional relationship between pressure and volume changes. Such tests are

unable to provide the necessary data to formulate three-dimensional constitutive equations.

On the contrary, triaxial tests reveal the deformation characteristics in three dimensions

under both hydrostatic pressure and deviatoric loading conditions. Such tests can furnish the

necessary data for the formulation of the general constitutive equations. In the present thesis,

a series of triaxial tests on alumina powders have been performed. Many interesting results

have been found. The experimental results will be presented in Chapter 4.

1.1.3 Models

There are some geomaterials-type models to describe the three-dimensional behavior

of frictional materials such as the Mohr-Coulomb, Drucker-Prager, Cam-Clay, Cap and Lade

models etc. These models are developed based on mainly the experiments on sand and clay.

In some models, the yield surfaces are a prior assumed and the associated flow rule is used.

In some other models, the non-associated flow rule is used and the irreversible strain

potential is assumed (see the following sections on the review of models). Most models are

time independent ones. It is difficult to use these model to describe the behavior of rock salt

(Langer [1988]). These models have a common disadvantage: they may not predict the

volumetric dilatancy behavior correctly even though some models predict the volumetric

compressibility pretty well (e.g., Desai & Siriwardane [1984] for artificial soil). For example,

the Cap model was modified and many additional constants were introduced in the model

by Chen and Baladi [1985]. The model still does not predict the volumetric dilatancy

precisely although the prediction matches well the other data. The Drucker-Prager model

overpredicts the dilatancy behavior. The Cam-clay model does not predict the dilatancy

behavior until the materials reach the failure stress and become softening. One reason for

these weak aspects is that the assumed yield surface is appropriate for some materials, but

it may not be good for other materials. Another reason is the skepticism of the normality in

the associated flow rule. Many experimental results do show that the behaviors of

geomaterials do not follow the normality conditions (Maier & Hueckel [1979]; Kim & Lade

[1988]; Jin et al. [1991]; Anandarajah et al. [1995]). The non-associated flow rule has been

explored by several researchers (Lade & Kim [1989a, 1989b]; Desai & Hashmi [1989]).

7

However, since in their models the potential surfaces are assumed in a priori form for some

materials, it is difficult to use these models for other type of materials. It is advisable to

examine the models carefully before using them.

There are also several models for the compaction of powders. Usually, the elliptical

yield curves or other type of curve in the (1,, Aj) plane are used, where I, and J' are the

first invariant of stress tensor and the second invariant of deviatoric stress tensor. These

models are mainly determined from the micromechanics analysis for porous metal materials.

The associated flow rule is usually used in these models. Nowadays, no formulation exists

regarding to the non-associated plasticity formulation for powder materials even though

some powder materials do have non-associated flow properties (Brown & Abou-Chedid

[1994]; Gurson & Yaun [1995]). It is important and necessary to explore the non-associated

flow rule. It is urgent and valuable to develop a new methodology to formulate non-

associated constitutive equations based on experimental results.

In section 1.2, we review most models which are often referred in the literature. Both

time dependent models and time independent models are reviewed. The advantages and

disadvantages of these models are highlighted. In Chapter 5 the elastic/viscoplastic model

is used to formulate the constitutive equation for alumina powder. The viscoplastic potential

is determined based on the experimental data presented in Chapter 4. Non-associated flow

rule is used in the model. In Chapter 6, the elastic/viscoplastic model is developed for

compressible materials. The orientation of irreversible strain rate tensor is determined based

on the experimental results on A16-SG alumina powder. A good agreement between the

model prediction and experimental data is obtained.

1.2 Literature Survey

In terms of the yielding or plastic flow behavior, materials can generally be classified

in two types. One type of material is the so-called nonfrictional material. For this kind of

material, such as metals, the yielding or plastic flow behavior is independent on the first

invariant of stress tensor. Several yield surfaces can be used to describe their yielding

behavior, as for instance, the Mises or Tresca yield conditions (see Malvern [1969]). An

associated flow rule is usually assumed for such type of models. The other type of material

is called frictional material. For such type of material, the yielding of plastic flow is

dependent on the first invariant of stress tensor. In other words, if there was friction within

the mass of materials, the frictional resistance will be proportional to the normal force. Sand,

clay, rock, concrete, powder-like materials etc. are falling into this category. The Mises yield

criterion is not any more suitable to describe yield behavior of this type of materials. Other

types of models should be used for such materials. In the following sections, we will review

several models for frictional materials and discuss the advantages and disadvantages of these

models. These models include Mohr-Coulomb model, Drucker-Prager model, Cap model,

Cam-Clay model and Lade model etc. More advanced models are also discussed.

1.2.1 Mohr-Coulomb Model and Drucker-Prager Model

According to the Mohr-Coulomb model, the shear strength (or stress at failure)

increases with increasing normal stress on the failure plane. The failure criteria can be

written as

T=c+otan#( (1.1)

where T is the shear stress on the failure plane, c the cohesion of the material, o the normal

effective stress on the failure plane and (P the angle of internal friction.

However, in the Mohr-Coulomb model, the intermediate stress is not taken into

account for the failure criteria. A generalization to account for the effects of all principal

stress was suggested by Drucker and Prager [1952] by using the invariants of stress tensor.

This generalized criterion can be written as

f= 2-aI,-K (1.2)

where a and K are positive material parameters. I, is the first invariant of stress tensor and

J2 is the second invariant of the deviatoric stress tensor. The associated flow (normality) rule

is assumed in this model. However, if the associated flow rule is assumed, the model usually

over-predicts the dilatant volumetric strain of the material comparing to the volumetric strain

obtained in experiments. Some problems can be solved without using the normality rule. If

the irreversible compressibility is taken into account, a cap is usually added at the open end

of the Drucker-Prager failure surface (see Chen & Baladi [1985]). If the Mohr-Coulomb

failure coincides with the Drucker-Prager failure condition at triaxial compression, the above

coefficients satisfy the following relations

2sin() 6ccos(F

V (3-sin) (3-sin() (1.3)

1.2.2 Cap Model

Based on experimental results, DiMaggio and Sandier [1971] and Sandier et al.

[1976] proposed a cap model. The yield surface for this model consists of two yield segments

(Fig. 1.1). One is a fixed yield surface and the other one is a yield cap surface. The fixed

yield segment can be written as

fA ')=0 (1.4)

There are several possible functions for f,. The expression for f, adopted by DiMaggio and

Sandler [1971] is given by

f, = --ye -'-a=0 (1.5)

and Desai and Siriwardane [1984] use

f,=c -ye -'-I,-a=0 (1.6)

for an artificial soil, where y, p, 0 and a are material parameters. The yield surfaces for yield

caps can be expressed as

f2(i, ,K)=0 (1.7)

where K defines the deformation history. DiMaggio and Sandler [1971] and Sandier et al.

[1976] have used an elliptic cap to represent the yield surface for the cohesionless material

they have considered

f =R 2J2+1-L)2-R2b =0 (1.8)

where Rb=(X-L) (Fig. 1.1), R is the ratio of the major to the minor axis of the ellipse, X the

value of I, at the intersection of the cap with the I -axis, L the value of J, at the center of the

ellipse, and b the value of (J2)n0 when I, =L. Because the two yield surfaces are intersecting

at I,=L, the X-L is related with f, through

X-L=R(ye -'+a+0I) (1.9)

The value of X depends on the plastic volumetric strain ep and is assumed as

X -- Inl--+Z, (1.10)

D WI

where D, Z, and W are material parameters to be determined. Usually the associated flow

rule is assumed for this model.

In the model, the cap surface may shrink or expand dependent on the increase or the

decrease of the irreversible volumetric strain. When the loading stress is on the cap surface

and the stress increment points outside, work-hardening takes place and the irreversible

volumetric strain increases so that the cap surface expands. When the loading stress is on the

failure surface, the volumetric irreversible strain decreases so that the cap surface shrinks.

Thus, this model can describe a softening behavior. However, for some frictional materials,

no softening behavior is present in the region considered. Sander et al. [1976] modified eqn

L X

z

Fig. 1.1 Yield surfaces for cap model.

(1.10) to prevent this softening behavior. Later, a great lot of progress has been made by

many researchers (see Desai & Siriwardane [1984]; Chen & Baladi [1985]). The model has

been used to describe the behavior of sand, artificial soil, concrete, rock and powder-like

materials etc.

The original cap model is further extended to account for rate effect by means of a

viscoplastic theory proposed by Perzyna [1966]. A lot of research have been done

(Zienkiewicz & Cormeau [1974]; Hughes & Taylor [1978]; Katona [1984]; Simo et al.

[1988]; Hoftstetter et al. [1993]). The model is nowadays extended to account for tensile

stresses and becomes more efficient when incorporated into the finite element method.

However, there are some limitations in this model, for instance, dilatancy is only

obtained when the stress state reaches failure surface. For some materials, this assumption

13

is not correct. Materials such as rock (Brace et al. [1966]; Cristescu [1989]), sand (Kim &

Lade [1988]) and some powders (Jin & Cristescu [1996b]; Jin et al. [1996d]), exhibit

dilatancy far away from the failure surface. Even though the model was modified by Chen

and Baladi [1985], it still would not describe correctly volumetric dilatancy. In other words,

this model could not give the reasonable prediction of volumetric strain at least for some

frictional materials. Another limitation is the assumption in eqn (1.10). It may not describe

correctly the work hardening behavior for some materials.

1.2.3 Cam-Clay Model

The Cam-Clay model was developed by a research group in Cambridge lead by

Roscoe (see Desai & Siriwardane [1984]; Chen & Baladi [1985]). The model is based on the

concept of Critical Void Ratio (or Critical Density). They found that the yielding of loose or

dense soil continues under both drained and undrained conditions until the material reaches

a critical void ratio. After the critical void ratio reaches, the volumetric strain does not

change any more. For normally consolidated soil, the stress state at this critical void ratio

reaches failure, while for over-consolidated soil, the stress state at the critical void ratio

reaches the residual stress after passing through failure and softening. After the critical void

ratio, the void ratio remains constant during subsequent deformations, that is, the material

will reach a state in which the arrangement of the particles is such that no volume change

takes place during shearing. This particular void ratio is called the critical void ratio. This

can be considered as the critical state of the material.

14

This idea can be expressed in the (p, q, e) space, where p is mean stress, q=o0-

o3=(3J2). and e is void ratio as shown in Fig. 1.2. The State Boundary (or Roscoe Surface)

can be used to describe the behavior of normally consolidated clay while the Hvorslev

Surface is used to describe the behavior of over consolidated clay (see Desai & Siriwardane

[1984]).

State boundary

or Roscoe

surface \

Fig.1.2 Critical State Line on the (p,q,e) space.

Fig. 1.3 shows a projection of the critical state line (CSL) on the p-q space together

with projections of typical section of the state boundary surface. The CSL is usually a straight

line passing through the origin. The projections of the state boundary surfaces are represented

by continuous curves and referred to as yield surfaces or yield caps. The yield surfaces can

15

be spherical, elliptical or of other shapes. In this model, the associated flow rule (normality)

is assumed.

If the yield surface is expressed as

-pMlu ( (1.11)

the model is called Cam-Clay model, which was first derived by Roscoe in 1958 (see Desai

& Siriwardane [1984]). If the yield surface is an ellipse

M2p2-M2pop+q2=0 (1.12)

the model is called Modified Cam-Clay model, which was first developed by Roscoe in 1968

(see Desai & Siriwardane [1984]), where M is the slope of the CSL and po is a work-

hardening parameter on the p axis. The above yield surface passes through (po,0) in the (p,q)

space. po can be expressed as a function of irreversible volumetric strain which is obtained

from hydrostatic tests only. It can be shown that the normal to the above yield surface at a

critical point is perpendicular to p axis, i.e., no volume changes at this point.

This model was used by many researchers (see Desai & Siriwardane [1984]; Chen

& Baladi [1985]; Britto & Gunn [1987]; Tripodi et al. [1994]) for normal consolidated clay,

over-consolidated clay, and other materials. This model was incorporated into several

commercial finite element programs (e.g., ABAQUS [1993]). However, this model has some

limitations. For some materials, the critical void ratio does not exist. For some materials, the

critical state line does not coincide with the compressibility/dilatancy boundary. Therefore,

it could not predict volumetric strain correctly.

yield cap

Fig. 1.3 Yield locus in q-p space, projection of Fig. 1.2 on q-p space.

1.2.4 Lade Model

In this model, Lade and Duncan [1975] first proposed the following function for the

failure surface, based on experimental observation for sand

1-K 113=0 (1.13)

where I, and I, are the first and third invariants of stress tensor, Ki is constant depending on

the density of sand. The expansion of the yield surfaces is defined by function f given by

If

f-- (1.14)

I3

17

whose values vary with loading and can reach the value K, at failure. In the model, the plastic

potential function is assumed as

Q=I-K213 (1.15)

where for a given f, K2 is a constant. The theory allows for a non-associated flow rule.

Subsequently, the theory was modified to include curved surfaces (Lade [1977]).

Later the theory was modified again (Kim & Lade [1988]; Lade & Kim [1989a]; Lade & Kim

[1989b]; Lade [1989]). The new model employs a single, isotropic yield surface. The yield

surfaces, expressed in terms of three invariants of stress tensor, describe the contours of total

plastic work. The non-associated flow rule is used in the model. The potential function is

also assumed as a function of three invariants of stress tensor. The elastic Young's modulus

is expressed as a function of the confining pressure. This model has been used to describe

the behavior of sand, clay and concrete etc. This model was also examined by Reddy and

Saxena [1992] for cemented and uncemented sand. They obtained the model from one

hydrostatic test and three conventional triaxial compression tests, then checked the model

with different loading paths. They concluded that the predicted results match reasonably the

experimental data except for the volumetric strain at the conventional triaxial tests.

1.2.5 Endochronic Model

The endochronic theory is derived from the laws of thermodynamics, i.e., the

conservation of energy and dissipation (Clausius-Duhem) inequality. It is based on the idea

of internal variable development. The theory was first proposed by Valanis ([197 la, 1971b,

18

1975]) by choosing the internal variable to be an intrinsic time scale. It was shown by

Valanis that mechanical properties such as hystersis and hardening in metals can be predicted

accurately by this model. Subsequently, Bazant [1977, 1978] has applied the model for soil

and concrete. Reddy and Saxena [1992] have used this model for cemented sand and

obtained the model with one hydrostatic compression test and one conventional triaxial

compression test. One function in the model was modified. The model was checked against

different loading paths. A good agreement between the tests and the prediction was obtained.

However, the endochronic theory is still in the development stage (see Yeh et al [1994]; Wu

& Ho [1995]) and very few cases of implementation of FEM are available. It was shown by

Sandler [1978] that in the present form, the endochronic model can cause difficulties in

numerical implementation, particularly in accounting for unloading, which may be remedied

by introducing internal barriers. In that case, the theory may exhibit some features similar as

the yield loading conditions from the plasticity theory, which the endochronic theory would

initially like to avoid.

1.2.6 Models for Compaction (Consolidation) of Powders

Based on the micromechanics analysis for porous metal, one kind of models for

porous metal or metal powder was first proposed by Green [1972] and Oyan et al. [1973].

The yield surface ( in this kind model can be expressed as a function of the first invariant

of stress tensor, I,, and the second invariant of deviatoric stress tensor, J2,,

0=A(T))J2+B(Ti)II( -(i)o=0( ,

(1.16)

19

where ri is relative density, A(Tr), B(rl) and 8(rl) are functions of relative density and oy is

yield stress of powder materials at full density. Later, Shima and Oyane [1976], Corapcioglu

and Uz [1978], Doraivelu et al. [1984] and Shima and Mimura [1986] also obtained this kind

of yield surface with different functions of A(T), B(r)) and 8(rl) based on the different

assumptions. In the work of Oyan et al. [1973], for instance, A, B and 6 are expressed as

A= 3; B 1

1+ 2 (1.17)

An associated flow rule is usually assumed in this kind of models.

The yield surface (1.16) is an elliptical surface in the principal stress space with the

major axis coinciding with the principal stress. In (1,J2) space, the yield surfaces are a family

of elliptical curves. Along each elliptical curve shown in Fig. 1.4 the relative density is

constant. The yield surfaces have symmetric properties for the tension and compression cases

and approach the von Mises yield surfaces when the relative density tends towards one, i.e.,

full density. The deformation history was incorporated into the model through the use of

relative density as a state variable. The dependence of deformation history can be generalized

(Brown & Weber [1988]; Krishnaswami & Trasorras [1995]; Tszeng & Wu [1996]) to

include a second state variable representing a scalar measure of the average plastic strain

experienced by the powder particles, thereby allowing the yield stress to change with

inelastic deformation. These models can be improved. However, several parameters are

involved, which become difficult to be determined from experimental data.

20

The model represented by eqn (1.16) was extended to viscoplasticity and thermo-

viscoplasticity (Abouaf et al. [1988]; Chenot et al. [1990]). The finite element methods were

used for the prediction of relative density variations in experimental powder compacts

(Nohara et al. [1988]; Jinka & Lewis [1994]; Krishnaswami & Trasorras [1995]; Oliver et

al. [1996]; Xu & McMeeking [1996]).

Fig.l.4 Typical yield surfaces for the model (1.16).

Another kind of micromechanics model was proposed by Gurson [1977]. For the case

of a spherical cavity within a perfect plastic matrix, the yield function can be expressed as

.2fos(1 1 1 ,

<)=3 +2fcosh I -1-f2=O (1.18)

y li 2a

21

where oy is the yield stress of the matrix and f is the void volume fraction defined by void

volume over total volume of porous metal material. This model is suitable to less porous

metal materials only. Other micromechanics models (Torquato [1991]; Fleck et al. [1992a,

1992b, 1995]; Nemart-Nasser & Hori [1993]) are also available. However, all these

micromechanics models are based on many idealized assumptions. Due to very few

experimental data available, it is difficult to validate the models for the complicated stress

states (Abou-Chedid & Brown [1992]; Brown & Abou-Chedid [1994]; Gurson & Yuan

[1995]).

1.2.7 The ElasticViscoplastic Models

The first viscoplastic model was proposed by Bingham in 1922 (see Malvem [1969]).

He considered the case of a simple shear in the x-direction and supposed that no motion takes

place until the stress T reaches a critical value (see Fig. 1.5), after which the magnitude of rate

deformation D is proportional to the amount by which T exceeds k. It can be shown

schematically by Fig.1.5. The viscoplastic equation for simple shear is assumed as

D 0 for T

2rlD= (1- )T for T>k (1.19)

T

where Tr is a viscosity coefficient.

T T

k

Fig.1.5 Bingham viscoplastic model.

Hohenemser and Prager in 1932 (see Malvem [1969]) assumed incompressibility and

generalized the equation to

f 0 for F<0

2lDq= Fo' for F>0 (1.20)

where oa is stress deviator tensor and F is measurement of overstress given by

k

F=l- ,

S(1.21)

where J2 is the second invariant of deviatoric stress. Subsequently, a lot of research has been

done in the viscoplasticity theory (Perzyna [1966]; Malvem [1969]; see comprehensive

review by Cristescu [1982]). One theory was proposed by Perzyna [1966]. The formulation

can be written as

S 1 1-2v .. F

e= --o +---- j+y*o (1.22)*

'i 2G ii E Go (1.22)

23

where G, E are elastic shear modulus and Young's modulus respectively, .ijis rate of stress

tensor, 6 is mean stress rate.

F- i-1 (1.23)

k

was suggested by Perzyna [1966]. In this theory, the associated flow rule has been used. F

should satisfy the convex conditions. Drucker's stable material postulate was used to define

F. Due to considerable research works (see Cristescu & Suliciu [1982]; Lubliner [1990]), the

elastic/viscoplastic theory can be reduced to the elasto-plastic theory if time effect is

negligible while the elasto-plastic theory can be generalized to elastic/viscoplastic theory if

time effect is taken into account (see Simo et al. [1988]).

Recently, one kind of elastic/viscoplastic model was proposed by Cristescu [1989,

1991, 1994]. The theory is based on the elastic/viscoplastic theory and the characteristics of

geomaterials (mainly, rock and sand). The basic formulation can be written as

1 1 AI / W(t)\9F(o)

e = +( I )61+k W(t J(o (1.24)

2G 3K 2G \H(o)l ao(

where K, G are elastic moduli, H is yield function, F is viscoplastic potential, W(t) is

irreversible work per unit volume, kT is a viscosity coefficient. This formulation is a

generalized form of the Bingham material. However, Cristescu [1991] gave a procedure to

determine the yield surfaces and potential surfaces based on a set of hydrostatic tests and

conventional triaxial compression tests. There is no a prior assumption regarding the yield

function and viscoplastic potential. This model has been used by many researchers (Florea

24

[1994a,b]; Cristescu et al. [1994]; Dahou et al. [1995]). The procedure to determine the yield

function and the viscoplastic potential was improved by Jin and Cristescu [1996a] and Jin

et al. [1996c] so that the model could better match the experimental results and be easily

incorporated into a finite element method. Based on the conventional triaxial compression

tests on alumina powder (Jin & Cristescu [1996b]; Jin et al. [1996d]) which will be presented

in detail in the present thesis, the improved procedure was used for the modeling of

consolidation of alumina powder A10 (Jin et al. [1996a, 1996b]). Another procedure to

determine the orientation of the irreversible volumetric strain rate was given by Cazacu et

al. [1996] and Cristescu et al. [1996]. They used a discontinuity function to describe the

effect of the compressibility and dilatancy boundary. A methodology to determine the

elastic/viscoplastic model in terms of the irreversible volumetric strain as a working

parameter is developed in this thesis. The model prediction matches well the experimental

data. This methodology is more attractive for the powder materials since the work hardening

parameter volumetricc strain) is directly associated with volumetric reduction, and the

procedure to determine the model is much simpler.

1.2.8 Other Models

There exist also many other models. Green elastic model (existence of a strain energy

function) and Cauchy elastic model (stress is assumed as a function of strain) can be used to

describe nonlinear elastic behavior (see Desai & Siriwardane [1984]). Hypoelasticity (the

increment of stress is expressed as a function of stress and increment of strain) originally

proposed by Truesdell [1955, 1966] could serve as a general nonlinear model (also see Green

25

[1956]). However these models also have some limitations (see comments by Chen & Baladi

[1985]).

There are also many publications regarding to the homogenization approach for

modeling (Christensen [1991]; Nemat-Nasser [1993]; Huang et al. [1994]). However, due

to complicated modeling and assumptions, this approach does not yet mature for practical

problems.

1.2.9 Remark on Recent Development of Plasticity-Based Models

For the plasticity-based models, the isotropic hardening and the kinematic hardening

rules were developed and used. Several theories were developed for work hardening. There

are several models for the kinematic work hardening (see comprehensive review by Desai

& Siriwardane [1984]; Chen & Baladi [1985]; Drucker [1988]). Sometimes, the isotropic

hardening rule and the kinematic hardening rules are incorporated into the above models

presented in the previous sections. However, due to the difficulties to determine the

parameters in the work hardening rule from experiments, the application of these hardening

rules is limited.

1.2.10 Finite Element Method

During the past several decades, the finite element method has rapidly become a very

popular technique for the numerical solution of complex problems in engineering with the

help of large-digit computers (see Bathe [1982]; Zienkiewicz & Taylor [1989]; Hughes

[1989]). Applications range from the stress analysis of structures to the solution of acoustical,

26

neutron physics and fluid dynamics problems. Indeed, the finite element method has been

developed so rapidly that the finite element process is established as a general numerical

method for the solution of partial differential equation systems subject to known boundary

and/or initial conditions.

The success of the finite element method is based largely on the basic finite element

procedures used: (i) the formulation of the problem in variational or weighted residual form

(weak form), (ii) the finite element discretization of this formulation, and (iii) the effective

solution of the resulting finite element equations. These basic steps are the same whichever

problem is considered. These basic steps provide a general framework and a quite natural

approach to engineering analysis in conjunction with the use of the large digit computers.

The basic concepts of finite element method and basic formulation for problems in

the isoparametric finite element representation can be referred in the books by Zienkiewicz

and Taylor [ 1989] and Hughes [1989]. The basic idea and procedures are easy to extend. The

basic formulation can serve for different problems, such as elastic-plasticity and

elastic/viscoplasticity problems.

In Appendix 1, we will give brief review of the discrete formulation of the

elastic/viscoplasticity theory for the finite element method. These formulations have been

implemented into a finite element method program. The program has been used for the

analysis of several problems in the present thesis. The main results will be presented in

Chapter 3.

1.3 Objective and Outline of the Present Work

The accurate constitutive formulation of nonlinear time dependent behavior of

materials is extremely important both for a theoretical analysis and for the finite element

simulation. With the existing models presented in the previous sections, it is difficult to

describe accurately three-dimensional nonlinear time dependent behavior, especially for the

nonlinear volumetric behavior from compressibility to dilatancy. Without an accurate

constitutive equation, there is no way to simulate engineering problems efficiently and

accurately when the nonlinear time dependent materials are involved. Thus, today it is a

challenge to develop a methodology to formulate nonlinear time dependent behavior of

materials accurately.

Cristescu [1989, 1991, 1994], as a pioneer, was trying to formulate nonlinear time

dependent behavior with the elastic/viscoplastic theory based on a set of triaxial tests. In his

formulation, the dilatancy and compressibility of volumetric strain are highly considered. The

yield surfaces and potential surfaces are determined from tests without any a priori

assumption. His formulation was used to describe rock salt and sand. However, his model

should be checked with more experimental data and examined from the point of view of the

application of the finite element method. In the present work, Cristescu's approach is

employed to establish the constitutive equations for rock salt and powder materials. We

found that some steps in his approach of the determination of the elastic/viscoplastic model

should be modified or changed in order to match better experimental data and to be easily

incorporated into the finite element method.

28

The main objective of the present work is to investigate the behavior of powder

materials experimentally by triaxial tests and to formulate three-dimensional constitutive

models for powder materials and for rock salt based on the elastic/viscoplastic theory and the

experimental data. Emphasis is placed on the improvement of the procedure (Cristescu

[1991,1994]) to determine the elastic/viscoplastic models. In the determination of the

models, the asymptotic behavior and the physical requirements of yield surfaces and potential

surfaces are incorporated into these procedures so that the yield surfaces and the viscoplastic

potentials are kept smooth. These models can better match the data and can easily be

incorporated into the finite element program.

The present work is associated with three aspects: (1) triaxial tests, (2) constitutive

modeling, and (3) finite element analysis. First, a constitutive equation for rock salt is

formulated based on the true triaxial experimental data obtained by Hunsche (see Cristescu

& Hunsche [1992]). Besides the properties considered in the Cristescu's approach, we have

considered the asymptotic behavior of yield surfaces and of the potential surfaces in the

procedure to determine the model. A new method to determine the viscoplastic potential

surface is proposed. Non-associated flow rule is used. The model matches very well the data.

Also, the model is incorporated into a finite element program. Several engineering problems

are analyzed with the program. The finite element analysis gives a good agreement with

experimental data. Such type of the elastic/viscoplastic model (both yield surfaces and

potential surfaces) is used for the first time in the finite element programs.

Secondly we have performed a series of triaxial tests on alumina powders using a

triaxial apparatus. The influences of the particle size and of the initial density on the

29

mechanical behavior of alumina powders are investigated experimentally. It is shown that

the volumetric strain behavior is strongly influenced by the particle size of powders. For the

alumina powder A16-SG of smaller particle size (0.4-1 pm), the volumetric strain always

exhibits compressibility in the range of applied stress state, while for the alumina powder

A 0 with larger particle size (40-200pm), the volumetric strain exhibits both compressibility

and dilatancy. The influence of the initial density on the behavior of alumina powder A10

has also been investigated. The elastic parameters of these alumina powders have been

measured according to the loading-reloading process. The set of triaxial tests furnishes the

necessary data for the formulation of three-dimensional constitutive equations.

Thirdly, a three-dimensional elastic/viscoplastic constitutive model for alumina

powder A10 is formulated based on the experimental data (a set of conventional triaxial

tests). The yield surfaces and the potential surfaces are determined. A new procedure to

determine the model is proposed based on such type of conventional triaxial compression

tests. The model is checked carefully against the experimental data with different loading

paths. A good agreement between the data and the model prediction is obtained.

Finally, another methodology to formulate the behavior of the materials which are

compressible only is proposed. The irreversible volumetric strain is taken as a work

hardening parameter instead of the irreversible stress work used in the previous method. The

methodology is used to formulate the behavior of alumina powder A16-SG. The irreversible

strain rate orientation tensor is used and determined in this model. It is found that this model

is much easier to handle and has fewer parameters involved than the previous ones. A

reasonable matching between the experimental data and the model prediction is obtained.

30

The organization of the dissertation is as follows: in the first Chapter, the models

existing in the literature are reviewed. In Chapter 2, the Cristescu's approach is outlined. A

new constitutive equation for transient creep of rock salt is proposed in Chapter 3. In the

fourth Chapter, three sets of triaxial tests on alumina powders (two sets for A10, one set for

A16-SG) have been performed and presented. In Chapter 5, an elastic/viscoplastic

constitutive equation for alumina powder A10 is formulated based on the hydrostatic tests

and the deviatoric tests presented in Chapter 4. In Chapter 6, a new methodology to

formulate the behavior of the materials which are compressible only has been proposed.

Finally, some conclusions and suggestions for the future work are given in Chapter 7.

CHAPTER 2

THE ELASTIC/VISCOPLASTIC (CRISTESCU) THEORY

2.1 The Elastic/iscoplastic Constitutive Equation

Many materials such as rock salt and powder materials can be considered as

homogeneous and isotropic ones. These materials have no preference directions. In many

cases, the deformations and rotations of particles are small. In the framework of the

elastic/viscoplastic theory, it is assumed that (1) the materials considered are homogenous

and isotropic; (2) the deformations and rotations of particles are small.

Based on the above assumptions, the total strain rate tensor e can be obtained by

adding the elastic strain rate tensor c and the irreversible strain rate tensor e :

"E "1

e=e +e (2.1)

For the elastic response of material, let us consider the fact that both longitudinal and

transverse extended body seismic waves can propagate through most materials (Cristescu

[1993a]). The fact suggests that most materials exhibit an elastic "instantaneous response."

The elastic parameters can be obtained from the speeds of longitudinal and transverse

extended body seismic waves. In the framework of the elastic/viscoplastic theory, the elastic

32

parameters can also be obtained by unloading process suggested by Cristescu [1989] for time

dependent materials. Thus, the instantaneous response can be expressed by a rate type elastic

relation

E 1 1

e- +(-. )61 (2.2)

2G 3K 2G "

where 6 is stress rate tensor, 6 mean stress rate, 1 unit tensor and G and K are elastic shear

and bulk moduli respectively.

For the irreversible part, it is assumed that the material follows a viscoplastic type

behavior (Cristescu [1991, 1994]), i.e.,

-k. 1 W(t) OF(o)

H(o)l o (2.3a)

where kT is viscosity coefficient, W(t) is total irreversible stress work per unit volume at time

t, F(o) is the viscoplastic potential and H(o) is the yield function with

H(o(t))=W(t) (2.3b)

the relaxation boundary or the equation of the stabilization boundary (where e =0, 6=0).

The bracket o denotes the positive part of a function, i.e.,

= (A +41)=A =A if A (2.4)

2 0 if A 0. (2.4)

If the yield function coincides with the viscoplastic potential, the flow rule is associated,

otherwise, the non-associated flow rule is assumed. For some materials, the viscoplastic

33

potential may not exist. In that case, the irreversible strain rate tensor is assumed to be

proportional to a tensor N (Cazacu et al. [1996]; Cristescu et al. [1996]):

e' k (l- w ) (2.5)

where N(o) is a tensor governing the orientation of irreversible strain rate tensor.

From eqn (2.1-2.3) and (2.5), the elastic/viscoplastic constitutive equation will be

written as

e 1 (- 1 )61+k 1 W(t) F().6a)

2G 3K 2G H(o) o (aa

or

= +( 1 1 +k I W(t) ) (2.6b)

2G 3K 2G H(o)b

In general, the yield function H and the potential function F (or N) are all dependent on stress

tensor. However, if H and F (or N) are assumed to be dependent on the first stress invariant

1, and on the second invariant J2 of deviatoric stress tensor only, disregarding the influence

of the third invariant of stress tensor, the whole constitutive equation can be determined from

a couple of triaxial tests (Cristescu [1991, 1994]). That is, H and F (or N) are assumed to be

dependent on mean stress

1 1

o=- -(1o+02+0 3) (2.7)

3 3

and octahedral shear stress r

t= -2 (2.8)

or equivalent stress

D=- ^3 (2.9)

2.2 Determination of the Elastic/Viscoplastic Constitutive Equation

In the general constitutive equation (2.6a) or (2.6b), there are three parameters, G, K

and kT, and two functions, H and F (or N) to be determined. Cristescu [1991, 1994] provided

a primary procedure to determine them from a couple of triaxial tests for rock salt and

dry/wet sand. The procedure is shown by scheme in Fig. 2.1. The details in the procedure

may be referred in Chapter 3 and 5.

2.2.1 Elastic Parameters and Viscosity Coefficient

Elastic parameters G and K are involved in all steps of the model determination. Thus

the elastic parameters must first be determined. Generally, there are two methods used to

determine elastic parameters, i.e., dynamic and static methods. In the former case, the seismic

wave velocities propagating through materials are measured, and from them are calculated

the elastic constants. In the latter case, the elastic parameters are determined using an

35

unloading procedure (Cristescu [1989]) for time dependent materials. The main idea of the

unloading procedure is to separate time effect from the pure elastic response during

unloading. For example, in stress-controlled tests, we keep stress constant at a chosen value

and allow the specimen to have a sufficient time to creep until the rates of strains are small

enough so that no significant interference between creep and elastic unloading would take

place during the unloading performed in a comparatively short time period. Elastic

parameters can afterwards be determined from the slope of the first portion of the unloading

curves (usually one-third or one-fourth of total axial stress). This procedure also can be used

for strain-controlled tests (see Chapter 4).

The viscosity coefficient can be determined from creep tests such as hydrostatic creep

tests after the yield surface and the potential surface are determined. If, for instance, a

hydrostatic creep test is performed, we integrate eqn (2.3) and obtain the following formula

for the viscosity coefficient

W (t) 'f H(o) 1

H(o)) l Fo(tf-ti) (2.10)

ao

where ti is the "initial" moment of the creep test, tf is the "final" moment, o(t,)=o(t)=

const., and H is yield function and F is potential surface. If we use a nonassociated

constitutive equation, then kT could be determined simultaneously with F. However, creep

tests can still be used for the adjustment of viscosity coefficient.

Functions & parameters

needed in the model Procedure used

Elastic parameters from both hydrostatic and deviatoric tests

G, K through unloading-reloading

Failure condition from deviatr t

--| from deviatoric tests }

H, from hydrostatic tests

Yield surface H- from deviatoqric ests

H

SAssumption: H=H,+Hd

SHydrostatic potential from hydrostatic tests .

G(o,T) function based on E, in deviatoric tests

Potential surface F

k F

F fkT F = fG(o,t)ds+ g(r)

g(r) based on deviatoric tests

Viscositycoefficient rom creep tests

Fig. 2.1 Schematic of the procedure to determine the elastic/viscoplastic model.

2.2.2 Yield Function

The yield function is defined as the irreversible stress work per unit volume on the

stabilization boundary or relaxation boundary according to the elastic/viscoplasticity theory

(see Cristescu [1967, 1989]; Cristescu & Suliciu [1982]; Lubliner [1990]), i.e.,

H(o(t))=W(t) (2.11)

Thus, the yield function can be obtained by evaluating the irreversible stress work per unit

volume along the relaxation boundary. W(t) can be computed as

W(t) =fo(t ):'(t)dt (2.12)

0

where (:) means double contraction on tensors. W(t) is used as an internal state variable or

as a work-hardening parameter. The constitutive equation (2.6) may be used even if the

material is assumed to have a zero initial yield stress.

The yield surface H is assumed as the sum of two parts, hydrostatic one H and

deviatoric one Hd, according to the loading paths used in the triaxial tests, i.e., H=H,+Hd.

H, is determined from the hydrostatic loading tests, while Hdis determined from deviatoric

tests by the regression of the irreversible stress work along the relaxation boundaries as

indicated in Fig.2.1. For different materials, the yield function might not be the same. For the

details of the determination of yield function, one can refer to Chapter 3 and 5 or refer

Cristescu [1991, 1994] for specific materials.

38

It is noted that most models in the literature (see Chapter 1) end at this step (see Fig.

2.1) and associated flow rule is assumed. However, most friction materials are not satisfied

with the associated flow rule. Non-associated flow rule should be used. The viscoplastic

potential should be explored.

2.2.3 Viscoplastic Potential

The viscoplastic potential F or the tensor N in eqn (2.6) governs the orientation of

irreversible strain rate tensor. It can also be determined from experimental data. Several steps

are required to determine the potential as shown in Fig. 2.1. According to the procedure,

some properties of the material such as a failure condition or/and compressibility/dilatancy

boundary are incorporated into the potential. The potential might not be the same for

different materials. In Chapter 3, 5 and 6, the procedure will be given in details for specific

materials. However, many steps have been changed from the original steps (Cristescu [1991,

1994]). The potentials obtained satisfy additional conditions which are required by physical

observation and by the finite element method.

If the yield function H(o) coincides with the viscoplastic potential, i.e., H(o) = F(o),

the flow rule is associated, while if H(o) F(o), the constitutive equation is called non-

associated. Both the associated and nonassociated flow rule can be used for modeling.

However, in general, the associated flow rule could not give good results, while the

nonassociated flow rule can.

2.3 Loading-Unloading Conditions

Let us assume that at an initial stress state (the so-called "primary" stress) oa=o(t0)

is an equilibrium stress state, i.e.H(o(to))=WP, with WP being the value of W for the

primary stress state. A stress variation from o(t) to o(t) with t>to is called unloading if

H(o(t))

loading, instantaneous elastic response given in eqn (2.2) is assumed. A stress variation

from (t) to o(t) (*o(t)) with t>to is called loading if H(o(t))>W(to) with one of the

three possible subcases defined by the following inequalities satisfied at (t):

dF 8F

:1>0 or >0 compressibility

do ao

OF 1F

-:1=0 or a=0 compress./dilatancy boundary (2.13a,b,c)

-:1<0 or -<0 dilatancy

do ao

The three inequalities correspond to the inequalities to be satisfied by the rate of irreversible

volumetric strain e, involved in eqn (2.3). Also, these inequalities define the region of

compressibility, compressibility/dilatancy boundary, and the region of dilatancy. These

concepts are very important not only in the theory but also in engineering applications

(Cristescu [1993a]).

CHAPTER 3

AN ELASTICNISCOPLASTIC MODEL FOR TRANSIENT CREEP

OF ROCK SALT

In the formulation of transient creep of rock salt, we start from a complete and

accurate set of true triaxial data on Gorleben rock salt obtained by Hunsche (see Cristescu

& Hunsche [1992]). The set of data consists of six true triaxial tests under the conditions of

different constant mean stress, i.e., o=8, 14, 25, 30, 35 and 40 MPa. For the formulation of

hydrostatic constitutive equation, one can refer to Cristescu and Hunsche [1992] and

Cristescu [1993a]. For the formulation of the mechanical behavior of rock salt in transient

creep, we follow Cristescu's approach (Chapter 2 and Cristescu [1989, 1991, 1993a, 1993b,

1994]). In this approach both compressibility and dilatancy properties of rock salt will be

considered, which is an advantage of the model over other models (See Chapter 1). The

concept of compressibility/dilatancy (C/D) boundary is introduced, thus one could well

determine with such a constitutive equation, where around an opening, for instance, the rock

becomes dilatant and where compressible. The procedure in this approach allows the

determination of a nonassociated elastic/viscoplastic constitutive equation able to describe

creep, relaxation, dilatancy and/or compressibility during creep, work-hardening and failure.

In this chapter a new elastic/viscoplastic constitutive equation for transient creep is

proposed. Compared with the early models (Cristescu [1993a, 1994]) the determination of

41

the yield function and of the viscoplastic potential is improved: the yield function has a

singularity at failure and the viscoplastic potential surfaces are requested to satisfy additional

conditions required to connect hydrostatic tests with deviatoric ones. A new procedure to

determine the viscoplastic potential is proposed. The viscoplastic potential is formulated in

simpler analytical expression than that in the previous papers. The present model on the one

side will be in a better agreement with the data and on the other side could be easily

incorporated into a FEM program. An explicit integration on viscoplastic strain components

is used in the constitutive equation. The results obtained by FEM are in good agreement with

the experimental data. The stress distribution around a vertical cylindrical cavern is analyzed.

Some interesting results are found.

3.1 Determination of the Elastic/Viscoplastic Constitutive Model

The procedure to determine the elastic/viscoplastic model has the following steps:

(1) determination of the elastic parameters, (2) determination of the C/D boundary and the

failure condition, (3) determination of the yield function, and (4) determination of the

viscoplastic potential and viscosity coefficient. Each step is directly related to experimental

data. We will discuss in details the above steps in the following sections.

3.1.1 Determination of the Elastic Parameters

The elastic parameters are determined by using the unloading processes in a stress

controlled apparatus mentioned in Chapter 2. For a stress controlled test, stress is kept

42

constant at a chosen level. The specimen is allowed to have a sufficient time to creep until

the rates of strains are small enough so that no significant interference between creep and

unloading phenomena would take place during the unloading performed in a comparatively

short time period as schematically shown in Fig. 3.1. Elastic parameters can afterwards be

determined from the slope of the first portion of the unloading curves (Fig. 3.1). For rock

salt, an average value of bulk modulus K=21.7 GPa and of shear modulus G=l 1.8 GPa have

been measured by Hunsche (see Cristescu & Hunsche [1992]). The values obtained by this

method are very close to those obtained by dynamic method (Cristescu [1989]). The elastic

parameters are calculated with the formula (3.1)

1 1

a, =arctg --- ; a=arctg

1G 9K6G

3G 9K 9K 6G

E2

CX

0 = constant

loading

Transverse strains Axial strain

Fig. 3.1 Procedure to determine the elastic parameters in

unloading processes following short creep periods

3.1.2 Determination of the Compressibilitv/Dilatancv Boundary and the Failure Condition

The compressibility/dilatancy (C/D) boundary is defined by eqn (2.13b). The

boundary has specific meanings in engineering applications (Cristescu [1989, 1993a]). Below

the boundary shown in Fig. 3.2 there is the compressibility region. In the compressibility

region, the rate of irreversible volumetric strain is positive (compression strain is

conventionally taken as positive), which means that the absolute value of irreversible volume

is reduced. The decrease of irreversible volumetric strain may close microcracks and pores

which may be present in geomaterials. In the dilatancy region shown in Fig.3.2, the rate of

irreversible volumetric strain is negative, which means that the absolute value of irreversible

volume is increasing. The increase of irreversible volume is directly related to the opening

of microcracks and of pores, or of any kinds of damage which may be present in

geomaterials. Sometimes in engineering applications, dilatancy is to be avoided. The C/D

boundary can be obtained straightforward from experimental data. In classical triaxial tests

where hydrostatic loading is followed by deviatoric loading under constant confining

pressure, the curve ao-e: is obtained in deviatoric tests. Here o and el are the stress and

volumetric strain during deviatoric loading with respect to the reference configuration at the

end of hydrostatic tests superscriptt R means "relative", i.e., with respect to the mentioned

reference configuration). The C/D boundary is determined as the stress loci where the slopes

(or tangents) of the stress-volumetric strain curves are equal to the elastic ones (i.e., e' =0).

In true triaxial tests where the hydrostatic loading is followed by a deviatoric loading under

44

constant mean stress, the C/D boundary is obtained directly from T eR there where the

slope of this curve is vertical.

For the rock salt from true triaxial tests, Cristescu and Hunsche [1992] have obtained

the C/D boundary as follows

o 02

X(o,r )=--i+ fi + f2- 0 (3.2)

O, j f a

where f,=-0.017, f,=0.9 and o=1 MPa(see Fig. 3.2 where the full line is just this

boundary).

Fig. 3.2 shows by diamonds the locus of maximum values of octahedral shear stresses

obtained in several tests using the data by Hunsche (see Cristescu & Hunsche [1992]); they

correspond to the mean stress values: o = 14, 20, 25, 30, 35 and 40 MPa, respectively. The

following equation can be used to approximate the locus of the maximum stresses

tf(o)=,-Yex2 Y3 (3.3)

where T, (o) is octahedral shear stress at failure, o is mean stress and y1=38.0 (MPa),

Y2=34.9(MPa), and y3=0.04 obtained through curve fitting. The stress states at failure and

the fitting curve are shown in Fig. 3.2.

35

(MPa) 30

Failure condition

25

dilatancy region

20 E <0

15 / =0

C/D boundary

10

5 ev >0

compressibility region

0 10 20 30 40 50 60

a (MPa)

Fig. 3.2 Domains of compressibility, dilatancy, C/D boundary and failure condition

3.1.3 Determination of the Yield Function

The yield function is defined as the stabilization boundary or relaxation boundary

according to the viscoplasticity theory (see Cristescu [1967, 1989]; Cristescu & Suliciu

[1982]; Lubliner [1990]). Due to the available short term experimental data on rock salt, it

is assumed that the stress strain curves are "almost" reaching the relaxation boundary for

transient creep. Thus the yield function can be obtained by evaluating the irreversible stress

work from the experimental data using eqn (2.12). W(t) and H(o,t) can be separated into two

parts according to the two loading paths used in true triaxial tests. In the first stage, the

hydrostatic loading, the stress components are kept equal (o0=o2=03), and increased

according to a certain conventionally chosen time interval t e [0, th). The second stage, the

46

deviatoric loading, takes place in the time interval t e [th, t] with o = constant, and it is only

the octahedral shear stress which is increased. Thus the complete expression of the yield

function is the sum of the yield functions obtained in the two stages (Cristescu [1994]):

H(o,T): =H(o)+Hdo,1) (3.4)

as is the total irreversible stress work:

W(t)=W,(t)+W(t) (3.5)

with the relations

Hh(o)=W,(t); Hd(o,T)= Wdt) (3.6)

where the subscripts h and d indicate hydrostatic and deviatoric stages respectively. Wh and

Wd can be computed from eqn (2.12) using the corresponding integral intervals.

From the calculation of irreversible stress work in hydrostatic tests and curve fitting,

Cristescu and Hunsche [1992] obtained the following expression for H,(a)

hoisin Wo--+(p +hi if as o(

H,(o): = Ga ) (3.7)

h+hI if a >0o

where ho0. 116 MPa, h=-0.103 MPa, o=2.880, (p=-62.60 and 0o=53 MPa.

Wd

(MPa)

0.1

S= 20 (MPa)

0.01

0.001 .I

0 5 10 15 20 25

T (MPa)

Fig. 3.3 Typical variation of irreversible stress work Wd with T for o =20 MPa,

test results points, a straight solid line suggesting a model behavior.

Using the experimental data obtained in the deviatoric loading stage, we calculated

the irreversible stress work using eqn (2.12). Fig. 3.3 shows a typical computed result for the

mean stress 20 MPa. The deviatoric irreversible stress work can be approximated by two

terms, for instance; one is linear in the exponential of T, the other has a singularity at failure.

Thus, the formula for the deviatoric part of the yield function becomes,

Hjo,t)=a(a) exp b(o)- + k(o)-exp 8( -1) (3.8)

Va.8

where

a(o) =a+ac-+..a ; b(o)=P, exp -P,- +

j r 0 [ *)(3.9)

k(o)=Kl a ; T/)=YY2 exp -Yz- .

There are two terms involved in the approximation of Hd(o,r) as shown in eqn (3.8).

First, from the log(Wd) ~ T plot of the experimental data (see Fig. 3.3), it was found that

besides a small region near failure Hd(o,t) can well be approximated by an exponential

function, i.e., Hd(o,t) a exp(bc). However, for a given o, near the maximum octahedral

shear stress t~m it was found that WfT) increases much faster than described by an

exponential function. It looks as if Wd(t) blows up near the maximum octahedral shear

stress. To capture this seemingly singular behavior of H,(o,T) near failure in a simple

manner, a first order algebraic singularity is introduced, HI(o,r) kl (tf-t), in which k and

tf are determined from the fit of the experimental data. Because this function decays too

slowly away from Tf, which would offset the close fit already obtained for the exponential

growth part, an exponential decay factor exp[8(r/Tr 1)] is introduced to reduce the error in

joining these two asymptotes. Thus the second part of Hd(o,t) becomes k/(Tz-)exp[8(t/r-

1)]. The factor 8 is introduced by error and trial method to fit better the curves. In general,

T" obtained from the fit of the data is a little larger than the actual maximum octahedral shear

stress in each of the experiments and the difference between T, and T are less than 1%.

Hence practically T, can be used, in lieu of max, to represent the octahedral shear stress at

failure as shown in Fig. 3.2.

12

(MPa) 40

10 calculated values from curve fitting

++++ calculated value from experiments

8

6 30

20 25

4-

14

2 a=8MPa

0 I

0 5 10 15 20 25 30 35

T (MPa)

Fig. 3.4 Fitting curves of irreversible work. Crosses experimental results,

solid lines prediction of eqn (3.8)

All the coefficients a, b, k and T, depend on o. After analyzing the data for each o,

the dependencies of a, b, k and T, on o are found out as shown in eqn (3.9). The following

constants are obtained:

a,=1.88x10-4(MPa), a2=1.13x10-4(MPa), 3=0.016x10- 4(MPa);

P,=0.566, P,=0.106, P3=0.217; (3.10)

K,=0.00171(MPa), K2=1.356 .

In the process of curve fitting, it was found that the data corresponding to 0=8 MPa

is not consistent with the trend established for o>8 MPa. Since during this test the mean

stress o is not really constant. For this reason, these data are further omitted. The data

corresponding to o=35 MPa were not used for the determination of a, b, k, and ,. Hence the

50

experimental data at o=35 MPa are used as an independent check to test the accuracy of the

Hd(o,r). Fig. 3.4 shows very close matching for o=35 MPa as it is for the data obtained with

other values of o.

3.1.4 Determination of the Viscoplastic Potential

After the determination of the yield function, we can use these results to get the

viscoplastic potential which governs the orientation of viscoplastic strain rate. The potential

has to satisfy several constrains from the elastic/viscoplasticity theory of

compressible/dilatant materials (Cristescu [1993a]). Let us shortly remind them. From eqn

(2.3) and (2.13), we have:

OF E

G(o,t): kT v (1

Oo I VW(t) \ (3.11)

H(o,t)

For the hydrostatic loading, G(o,z) satisfies the following condition

G(o,t)I,= -kT--F =4(o) (3.12)

a IT1=0

where 4(o) is determined from hydrostatic tests, while for general loading paths from eqn

(2.13), G (T,o) is required to satisfy the following conditions:

G(o,t) >0 X(o,t)>0 compressibilityy)

G(o,T) =0 X(o,r)=0 (C/D boundary) (3.13)

G(o,t) <0 X(o,T)<0 dilatancyy)

51

where X(o,t)=0 is C/D boundary. The rate of irreversible volumetric strain will increase in

compressibility domain, decrease in dilatancy domain and be equal to zero at the boundary

between compressibility and dilatancy regions. Also the rate of irreversible volumetric strain

may reach negative infinite values when stress state is approaching the failure condition

denoted as Y/o,r)-=O. The simplest form which has the above properties is

aF X(o,T)(o)

kT a( Y(or) (3.14)

a Y/or,t)

where the function P is related to hydrostatic behavior by

X(o,0) t(o)

Y0o,0) 3.15)

Cristescu [1991, 1994] has employed this form to get the viscoplastic potential. Although

this form satisfies exactly eqn (3.13), it may not match well enough the data within the

compressibility and dilatancy regions and may cause a discontinuity of e' at T-r. In order

to match better the data within the mentioned regions and ensure the continuity of e', we

propose a procedure to determine the irreversible volumetric strain rate directly from

experiments, i.e., to determine the function G(o,t) by fitting the experimental data with the

requirement that G has to satisfy the conditions (3.13).

Let us start from eqn (3.12). G(o,T) may be divided into two parts according to two

physical meanings: compressibility and dilatancy, denoted by G, and G, respectively. Each

part of G(o,T) is considered to be an asymptotic approximation and to be a dominant part in

52

the corresponding region. In order to satisfy the condition G(o,t) = 0 at the C/D boundary,

we take G,(o,t) to be strictly equal to zero and G, to be approximately equal to zero on this

boundary in the sense that G, is dominant in the compressibility region as compared with G,

From eqn (2.6), aF/Ta should be zero along and in the immediate neighborhood of =-O in

order to obtain a correct strain variation along a loading path very close to the hydrostatic

one, i.e., an additional very small deviatoric stress was superposed. This condition is imposed

in order to avoid the presence of a vertex of the surface at the hydrostatic axis and to keep

the continuity of e when passing from hydrostatic to deviatoric loading. Since aG/at is also

related (see the following sections) to aF(o,t)/at, at this stage we take aG/at=0 at T=O. Thus

we take G,/ T-=0 and aG2 / 1 =0 along =-0. Depending on the evaluation of the right-hand

side term of eqn (3.11) using the data in the compressibility region, we will be able to choose

an appropriate form for G, in order to match the data and satisfy the above two conditions.

However, since in the compressibility region the experimental data available have not had

enough accurate digits, the results obtained from the right-hand side of eqn (3.12) show a

large scatter in that region. For this reason we use the simplest form for G,, i.e., a second

order polynomial in t as

Gi(o, )= Q(o) -e(o),, (3.16)

where as before a.= 1 MPa. From the requirement 9G,(o,t)/r-=0 at T-O, follows that a

possible first order term in eqn (3.16) is zero.

53

From hydrostatic tests Cristescu and Hunsche [1992] have obtained for 4(o) involved

in eqn (3.15) and (3.16):

o (-q2 /

()= a if o o00 (3.17)

0 if 0>00

where q,=1.5 x 108 (s l), q2=o, = 53 MPa.

From the above conditions, G,(o,T)lj=0= (o) and from the requirement G,(o,T)=0

along X(o,t)=0, we can determine Q[ and Q2 from eqn (3.17) (note f2 / f, = -q2 / .) as

,(o)=q o-q2 Qq2 () q1

o o)( f2 (3.18)

In the dilatancy region, the function G2 mentioned above should be dominant in that

region, but must have much smaller values than G, in the compressibility region in order to

obtain an asymptotic approximation of the conditions (3.13), Higher order polynomials in

z seem to be good expressions for the approximation of G, in the dilatancy region. Based on

the data available and the approximated evaluation of the right-hand side of eqn (3.11) (see

Cristescu [1991]), we can choose the 5th and 9th order powers to fit the data (Fig. 3.5),

G2(Fo,)=A(O)(i J +B(o)( -J (3.19)

G. a

54

A(o) and B(o) can be determined using the data from true triaxial experiments performed at

several mean stresses. When analyzing the data and trying to choose proper polynomial

functions, we have found that at the beginning of the dilatancy region, by combining G,(o,t)

with a fifth order power function we get a better fitting of the experimental data. Thus, the

A coefficients are determined from several tests performed with various constant values of

a. By combining a ninth order power function with G,(o,T) and the above fifth order

function, we were able to obtain a better fitting of the curves in the dilatancy region at higher

value of T. Thus, the B coefficients are determined. After obtaining the values of A and B

from the tests performed under different mean stresses, we use the least-square method to get

the expression of the functions A(o) and B(o),

A(o)a B(o)=a9 (,) (3.20)

where a5=-5.73x 10-7 (s'), b5= -2.12, a,=-1.63x 108(s-'), b9=-4.77.

In the neighborhood of failure, the dilatancy behavior may become singular.

However, the singularity is very difficult to make precise from the data shown in Fig. 3.5.

Thus, for the expression of G we use both G, describing compressibility and G2 describing

dilatancy. From eqns (3.16) and (3.19), it is easy to show that in the compressibility region

G, is a dominant function as compared to G2. One can use a parameter

p=Ic G2 da/fc Glda (where domain C is the compressibility region) to decide which

term is dominant. The value p=0.11 of this parameter indicates that G, is dominant in the

compressible region. The function G is approximated by

(3.21)

G(o,T)=G,+G2=Qi,(o)+Q2() + A(o) -I '+B(o) -

Then, we integrate eqn (3.11) with respect to a to get

k7F(o,t)=fG(T,o)do+g(T)=F,(o,T) +g(r)

where

F,(o,r):=fG(o,t)do .

(3.22)

(3.23)

0.005

G(o,T)

(s-1) 0

-0.005

-0.01

-0.015

a= 14MPa \

-0.02

20

25 o 30

-0 .025 i I I I

0 5 10 15 20 25 30 35

T (MPa)

Fig. 3.5. Data for various confining pressure shown used to determine the function G(o,t),

symbols experimental data, solid lines prediction of eqn (3.21).

56

In order to obtain g(t), we differentiate eqn (3.22) with respect to T and combine the

irreversible strain components from eqn (2.3)

FW() F 8F aT

S I W) ) (3.24)

SH(o) ~ a o o (3.

Finally we get the formula for g '(s):

1I 1 / 2 *.. ._ .2 *,'i .2 OF

r)=(eI-E2) +(E), (3-e) -" (.5F,

W(t) r2 3 (3.25)

H(o)

In eqn (3.25), all terms are known besides the strain rates involved in the square root

which are obtained from the experimental data (Cristescu [1991]). For the determination of

g'(T) we plot Fig. 3.6 using the data for several constant mean stresses. It was found that g'(r)

is independent on mean stress as predicted by eqn (3.25), which implies that the viscoplastic

potential exists. Thus the procedure to find the potential is logical.

The solid line curve in Fig. 3.6 is an approximation for g'(T), which can be described

by a polynomial function:

g'(C +. (3.26)

where g=0.0013(MPa s '), g,=6.5x106(MPas -), g,=2.8xlO-3(MPa s -).Thusg(r)

can be written as

0.3

g'(C)

(sl) o MPa

cr MPaa

8

0.2 14

20

30

40

0.1

0 5 10 15 20 25 30 35

t (MPa)

Fig. 3.6 Data for various confining pressure shown used to determine g'(t);

symbols test results, solid line prediction of (3.26).

1 aT I3 a 9 9

) g, -- + & i 4 (3.27)

It is noted that in eqn (3.26), we start with the first order power. Since aF/8l should

be equal to zero at =-0, we let g'(t)=0 at T--0 in order to obtain the continuity of the

derivatives of the viscoplastic potential.

Therefore the viscoplastic potential for Gorleben rock salt is completely determined

based on the experimental data and can be written as follows

S 014 2q1 o 2 2 12 2 o

kfl(o,t)=o 4 -- i

4 ,J 3 oI a 2 ,) f2t o,) (3.28)

Using the expression (3.28), we get the shapes of the viscoplastic potential surfaces

while by using eqn (3.4) for H(o,T) we can plot the shapes of the yield surfaces. Fig. 3.7

shows the two sets of curves. It is clear that these two families of curves are quite distinct

both in the compressibility and in the dilatancy regions. Let us note also that the C/D

boundary is quite distinct from the curve aH/Oo=O. Thus, an associated flow rule cannot

predict correctly the C/D boundary. Let us remind that the yield function and the viscoplastic

potential are determined from the data following two distinct procedures. It is found that

FsH. Thus, an associated flow rule is not suitable for this rock salt. Non-associated

constitutive equation should be used with the potential furnished by eqn (3.28).

3.2 Comparison with the Data

One simple way to test a model is to compare the experimental data with the model

prediction. For the model developed in the previous sections, we use three different criteria

to check it. First, this can be done by trying to reproduce with the model the kind of tests

which have been used to formulate the model. Secondly, as an independent check we try to

40

(MPa) 35

awao= 0 ^_ r- 2.0

30 .o05--

5 fail --5 5.02-

25

20

15-

0.1 0.2 CD

5

0 10 20 30 40 50 60

o (MPa)

Fig. 3.7 Shape of yield surfaces (thin solid lines), potential surfaces (thick solid lines),

C/D boundary (dash-dot line) and lWH/o=0 (dashed line).

match a set of data which have not been used to formulate the model. Third, we incorporate

the model into a finite element program and try to describe triaxial tests. In order to

reproduce with the model the data which have been used to formulate the model, let us

consider constitutive equation (2.6a). It will be assumed that stresses are increased by small

successive steps, according to the same law as in the experiments done to establish the

model, and with the same global loading rate as in the experiments (Cristescu [1989, 1991,

1994]). At each very small loading step, the stresses are assumed to increase instantly at time

to, say. Afterwards a creep follows under constant stresses in the time interval t-t,. We can

integrate eqn (2.3) by multiplying it with o to get the following equation for the irreversible

stress work

W(t)=H(o(t))-(H(o)-W(to))exp (to t) (3.29)

The total strains as a result of one stress step can be written as follows

e(t) =E) H 0 1 -exp 8oaF (tot) (3.30)

--o---:F L HO J

H oF

with the initial conditions

e' =0

t=to : e )1 1 1 (3.31)

S -3K 2G 2G

where o(to) is the initial relative stress at each stress step.

The above constitutive equations can be used to reproduce the experimental results.

The viscosity coefficient kT can be determined from creep tests and used as a parameter

which takes into account the speed of the tests as mentioned in the paper by Cristescu [1991,

1994]. The constitutive equations are used to predict the stress-strain curves for the cases

o=14, 0=20, o =40 and 0=35 MPa. The prediction of the model matches well with the

experimental data as shown in Figs. 3.8a-3.8d where circles correspond to the experimental

61

results, and the solid lines to the model prediction. As an independent check, for the

experimental data at 0=35 MPa which were not used to formulate the constitutive equation,

there is also a very close agreement between experimental data and the predicted results (Fig.

3.8d).

20-

S 18

(MPa)

16

14

12

10

8

6

4

2

0

-0.06

-0.04 -0.02

0 0.02 0.04 0.06

Fig. 3.8a Experimental stress-strain curves (circles) compared with

predicted results (solid lines) for a = 14 MPa.

25

(MPa) E, Ev

20

15

S= 20 MPa

10

5

0 I I

-0.04 -0.02 0 0.02 0.04 0.06

Fig. 3.8b Same as in Fig. 3.8a. but for o=20 MPa

35

E2 EF

(MPa)

25

20

15 H= 40 MPa

10

5

-0.1 -0.05 0 0.05 0.1 0.15 0.2

Fig. 3.8c Same as in Fig. 3.8a but for o=40 MPa (note: the data e2 are slightly

different from those for E3)

30

T : E2 Ev 000

(MPa) 25

20

Y= 35 MPa

15

10

5

-0.1 -0.05 0 0.05 0.1 0.15 0.2

Fig. 3.8d Same as in Fig. 3.8c but for o=35 MPa as an independent check

3.3 The Finite Element Analysis

3.3.1 Formulation of the ElasticViscoplasticity Theory

The formulation of elastic/viscoplasticity theory in discrete time form is presented

in the book by Owen and Hinton [1980], with a classic approach by truncated Taylor series

for the rate of viscoplastic strain (Zienkiewicz & Cormeau [1974]; Hughes & Taylor [1978];

Owen & Hinton [1980]; Marques & Owen [1983]; Szabo [1990]; Desai et al. [1995]). There

are many finite element programs with source code available, such as VISCO (Owen &

Hinton [1980]) programs. These codes are usually for education and research purpose, which

are explained in detail. There are also many commercial finite element programs without

64

source code in public, such as ABAQUS [1993]. Several standard open connections in these

commercial programs can be used for users to write their own subroutines, but the open

connections are limited. In order to verify the model developed above, we use VISCO

program (Owen & Hinton [1980]) as our basic program. Several subroutines were

additionally developed for the implementation of the model. A brief description of the finite

element algorithm for elastic/viscoplastic model is given in appendix A.

During the computation, the process of time step marching is stopped at certain time

or when the stresses at all Gauss points satisfy H=W(t) and continues in the next loading

step. The time step At., at step (n+l) was selected subject to the following empirical criteria

A -fn -; Ati lKAta (3.32)

where I, and 1 are second invariants of the viscoplastic strain and strain-rate tensors

respectively; 0 and K are specified constants to control time step. The first criterion selects

a variable time-step size such that the maximum effective viscoplastic strain increment

occurring during next time step is a fraction of the total effective strain accumulated before.

At,,- is evaluated at each Gaussian integration point and the least value is taken for

computation. The second criterion imposes a restriction on the variable step size between

successive intervals calculated by the first criterion to prevent oscillations in the solutions

as steady-state conditions are approaching.

In the following examples, only explicit integration scheme for viscoplastic strains

is employed. For nonassociated flow rule, the stiffness matrix [Kt"] is symmetric and does

65

not vary. The explicit scheme is simpler and easy to implement into the finite element

program. However, the explicit scheme may not be stable when time steps become large. The

restrictions expressed by eqn (3.32) have to be imposed on the scheme so that the explicit

scheme is kept stable. The instability problem can be overcome by implicit schemes.

3.3.2 Example I: Axial Compression with Confining Pressure

In true triaxial tests for rock salt (see Cristescu & Hunsche [1992]), the specimen is

cubic. Stresses in two horizonal directions remain equal during the tests. Thus this kind of

tests will be approximated by an axial symmetrical test where the friction between the

cylindrical specimen and the piston is neglected. Due to the axial symmetry of the cylindrical

specimen, we can use a quarter of the specimen only in order to create a finite element mesh.

6 Serendipity quadratic elements and 29 nodes are used as shown in Fig. 3.9. Each element

has eight nodes. The boundary conditions are also shown in Fig. 3.9.

We will use the same loading procedure as used in the experiments. Hydrostatic

loading is applied first and afterwards the deviatoric loading is performed under constant

mean stress. Let us consider the test for the case o-40 MPa. The total hydrostatic loading is

divided into 10 steps. The computation remains stable for all the steps. The total deviatoric

load is also divided into 10 steps with unequal values. In each loading step, the computation

is carried on until the condition H(o(t))=W(t) in every Gaussian integration point or a

convergence condition is satisfied. As a convergence condition we use here the ratio of the

second invariants of irreversible strain rate evaluated at the initial and the current moment.

The tolerance chosen is 0.1%. In the first eight steps shown in Fig. 3.10, the scheme is stable

66

and convergent. However, when we apply the ninth stress step, the scheme becomes

unstable and blows up. In this computation, Q-0.0009 and = 1.1 were used. The results are

shown in Fig. 3.10. The results obtained with the finite element program give a reasonable

matching with the experimental results although the scheme in the finite element program

blows up near failure. The instability problem in finite element programs can be overcome

by implicit schemes (see Hughes & Taylor [1978]; Simo & Govindjee [1991]).

"o---

Fig. 3.9 Mesh of a quarter of cylindrical specimen used in the finite element analysis.

3.3.3 Example II: Stress Analysis for A Cylindrical Cavity in Rock Salt

Let us consider now a cylindrical cavity such as a borehole or a shaft with the same

initial or primary stress in all horizonal directions (oa=o h). For simplicity, primary

hydrostatic stress is assumed, i.e., the primary stress in the vertical direction o, is taken to

be the same with that in the horizonal direction. At excavation, the diameter a of the cavity

is assumed to be 1 meter. We make some idealized assumptions: the problem is a plane strain

one with primary stress x =hy =av= 10 MPa in the first example and 30 MPa in the second

example, and that the primary stress is instantaneously released at the cavity surface due to

excavation at time to. Afterwards the rock salt around the cavity creeps according to the

above constitutive equations. We use a quarter profile of the cavity for our computation. The

corresponding boundary conditions are shown in Fig. 3.11. Here the cylindrical coordinate

r is used for the distance from the axis of the cavity to some point inside the rock mass (e.g.,

r=a represents the surface of the wall). Two meshes are used for computation. An area with

7 m in length and 7 m in width is first considered with 253 nodes and 72 elements as shown

in Fig. 3.11. Each element has 8 nodes. In the second variant, we enlarge the area to 20 m in

length and 20 m in width and refine the mesh to a total of 333 nodes and 96 elements. The

computed results obtained with the two meshes do not differ significantly. Figs. 3.12-3.16

show the results obtained with the finer mesh.

In the computation, we choose for the restriction parameters in eqn (3.32): =-0.01

and K=1.3. The computation is stable and convergent. Figs. 12-15 show the variations of the

stresses and octahedral shear stress with the distance from the surface of the wall for the two

cases. A comparison with the elastic solution of the ultimate elastic/viscoplastic one (i.e.,

35

T 30 ev El

(MPa) E2

25

20

15 / 40 (MPa)

10

5

-0.1 -0.05 0 0.05 0.1 0.15

E

Fig. 3.10 Comparison between experimental data (solid line) with FEM results (triangles).

computation was carried out up to stabilization) shows a significant change in stress

distribution. The circumferential stress oe decreases near the wall of the cavity and then

slightly increases at farther distances (Fig. 3.12 and 3.14). The vertical stress oz and the radial

stress o, decrease with respect to the elastic solution. The octahedral shear stress t is smaller

near the surface of the wall, but greater at farther distance (Fig. 3.15 and 3.15). At greater

distances from the surface of the wall, the difference between the two solutions becomes

negligible and the stresses approach asymptotically the primary values. In the case of initial

stress 10 MPa, all regions are compressible. In the case of initial stress 30 MPa, there is a

dilatancy region expanding up to about 0.4 m from the surface of the wall. The other regions

around the opening are compressible. The displacements inside the rock vary in time. Fig.

69

3.16 shows the variation of the displacement at the surface of the wall, starting from the

elastic "instantaneous" response up to stabilization of transient creep. The upper curve

corresponds to the case of initial stress 30 MPa and lower curve to the case of initial stress

10 MPa. A dimensionless quantity kTxt is used for the "time" parameter, where k, is the

viscosity coefficient and t is time. If kT is determined from creep tests (Cristescu 1989), the

time t can also be made precise. Here the intention was to give an example only, obtained

with some simplifying assumptions in order to illustrate the use of the model in a mining

problem. More elaborate engineering problems could also be analyzed.

Fig. 3.11 The finite element mesh used for a vertical cavity.

(MPa)

(MPa)

1 2 3 4 5 6

r/a

Fig. 3.12 Stress distribution with distance from the surface of the cavity wall for the

primary hydrostatic stress 10 MPa (thin lines-instantaneous elastic solution; thick

lines-ultimate elastic/viscoplastic solution obtained when (t-to)xk.=6.6).

8

(MPa) 7

6

5

4

3

2

0

1 3 5 7 9 11 13 15 17

Fig. 3.13 Variation of the octahedral shear stress with distance from surface of the cavity

for the primary hydrostatic stress 10 MPa (thin line instantaneous elastic

solution; thick line ultimate elastic/viscoplastic solution obtained when (t-to)xkT=6.6).

71

60

a

(MPa)

50

40

30

20

10

1 4 7 r/a 10

Fig. 3.14 Same as in Fig. 3.15 but for the primary hydrostatic stress 30 MPa,

(ultimate elastic/viscoplastic solution obtained when (t-to)xkT=5.6).

25

(MPa)

20

15

10

5-

0

1 4 7 10 13 16 a

Fig. 3.15 Same as in Fig. 3.16 but for the primary hydrostatic stress 30 MPa,

(ultimate elastic/viscoplastic solution obtained when (t-to)xkr=5.6).

0.004

u/a

0.0035

0.003

0.0025

0.002

0.0015

0.001

0.0005

0 ....

0 1 2 3 4 5 6 7

kT x (t-to)

Fig. 3.16 Variation of relative radial displacement (u/a) of the surface of the wall with

dimensionless "time" kTx(t-to) (upper curve, for primary stress 30 MPa;

lower curve for primary stress 10 MPa).

3.4 Discussion and Conclusion

In present model transient creep is taken into account only. Since this model is a

phenomenological one based on short term tests, it may underestimate the strain magnitudes.

Actually, stationary creep should also be taken into account when one wants to predict stress

or displacement distributions around some underground openings. Stationary creep may also

contribute to the stress redistribution around an underground opening due to the nonlinear

behavior of rock salt. The finite element program for both transient and stationary creep

should be used in this case. The possible influence of kinematic hardening will also be

considered in the future papers, if enough experimental data would be available. Temperature

73

and humidity are other important factors which influence the behavior of rock salt. We have

not considered them in the present model.

In this formulation, non-associated flow rule was used. It was obtained from

experimental data. Many models adopted nonassociated flow rule (Desai et al. [1987]; Kim

& Lade [19881; Lade & Kim [1989a, 1989b]; Klisinski et al. [1992]) although non-associated

flow rule has been found to cause some problems (Sandier & Rubin [1987]).

Using the triaxial experimental data on Gorleben rock salt, a new nonassociated

elastic/viscoplastic model for transient creep is proposed. This elastic/viscoplastic model can

be used to describe creep, relaxation, dilatancy and/or compressibility and work hardening.

Singularity and asymptotic properties of the yield function and of the viscoplastic potential

which has to satisfy several restrictions are considered. A new procedure to determine the

viscoplastic potential is developed. This procedure allows to better fit data and keep the

viscoplastic potential smooth. The examples given show that the model predicts quite well

the experimental data. Since both the yield function and viscoplastic potential are smooth

functions, the present model can be incorporated quite easily into a finite element program.

An explicit integration scheme with very small time steps for viscoplastic strain components

is used. Several examples show that the model is easy to use in the finite element program.

Instability problems during our computation have not been found except the unstable

computation caused by explicit scheme for viscoplastic strains although non-associated flow

rule has been found to cause some problems (Sandler & Rubin [1987]). Implicit integration

schemes may be better to overcome unstable computation problems.

CHAPTER 4

TRIAXIAL EXPERIMENTAL RESULTS ON ALUMINA POWDERS

In this chapter, the mechanical behavior of alumina powder has been investigated

experimentally under hydrostatic and deviatoric loading conditions in a triaxial apparatus.

Two types of alumina powder, as received A10 and A16-SG, were tested. The specimens

were fully saturated with water. The elastic parameters of alumina powders were measured

based on the loading-reloading processes following a creep or relaxation period. The

experimental results reveal that the mechanical behavior of alumina powder is strongly

dependent on the particle size and on the initial density. The deformation of alumina powder

is time dependent and stress history dependent. The application of deviatoric stress can

produce additional consolidation or dilatancy. The compressibility/dilatancy boundaries and

the failure conditions have been obtained for the tested powders.

Also, the experiments have been performed according to the requirements presented

in the previous chapters to formulate a three dimensional constitutive equation. Thus the

experimental results provide a set of data for the formulation of three-dimensional, either

viscoelastic or viscoplastic, models.

In section 4.2, the basic testing procedure is presented and discussed, which consists

of four sequential steps: (1) specimen preparation, (2) back pressure saturation, (3)

hydrostatic consolidation or hydrostatic loading, and (4) shearing (deviatoric loading). In

75

section 4.3, the experimental results under hydrostatic loading conditions are presented. The

results under deviatoric loading conditions are presented and discussed in section 4.4.

4.1 Introduction

Forming processes are crucial in the ceramics industry today. Among various forming

processes, casting is commonly used for the forming of products. One of the casting

processes is pressure casting. This process has been investigated as a forming technique for

porcelain, complex refractory shapes or hard ferrite magnets, etc. The casting time is usually

controlled by regulation of the external pressure (Reed [1988]). In pressure casting, the mold

serves as a filter. The externally applied pressure is usually increased up to 1.5 MPa. Water

in the cast is squeezed out. This process reduces the water content significantly so that the

drying shrinkage can be reduced. However, pressure casting of slurries or pastes in a mold

creates density gradients due to the cast geometry or the friction between the materials and

the wall of the mold. In other words, the density distribution is not uniform in the cast. This

non-uniform density distribution will cause problems such as warping in the later sintering

process. In order to understand the process and avoid non-uniform density distribution in a

mold, the consolidation behavior, or stress strain relationships, of slurries or pastes must be

investigated.

Slurry or paste of alumina powder fully saturated with water is often used for pressure

casting. Therefore, the consolidation behavior of alumina paste plays an important role in

pressure casting. Whether or not we can get good products from this kind of casting depends

76

on how well the consolidation behavior of powder or powder paste is understood. To

understand the consolidation behavior of alumina paste at different stress states, we need to

know the three dimensional stress strain relationships, i.e., constitutive equations, of these

powders. To formulate constitutive equations, triaxial tests are often used.

Triaxial tests have been successfully used in civil engineering for geomaterials such

as sand, clay, and rock to obtain general three dimensional constitutive equations (Bishop &

Henkel [1962]; Cristescu [1989]). However, powders such as alumina powder or alumina

paste are different from sand or clay. Usually powder possesses more pure constituents and

surface chemistry of particles is also involved (Reed [1988]). Hence the mechanical

properties of powders cannot be obtained directly from those of sand, clay, or rock. Several

triaxial tests for metal powders or dry ceramic powders at high pressure were reported

(Shima & Mimura [1986]; Brown & Abou-Chedid [1994]; Gurson & Yuan [1995]).

Unfortunately, very few triaxial tests have been performed for ceramic paste such as alumina

powder (Stanley-Wood [1988]). Triaxial tests are quite different from uniaxial tests. In a

uniaxial confined test, powder is placed into a rigid cylindrical mold or die and pressure is

applied along axial direction using a piston. Such tests reveal mainly one dimensional

relationship between pressure and volume changes (Gethin et al. [1994]; Chen et al. [1994];

Shapiro [1995]). However, due to the friction between the powder and the wall of a mold,

there exist very complicated stress states in the mold and the final density distribution is

nonuniform. Consequently, uniaxial tests cannot provide good basis for the formulation of

three-dimensional relationship between stress and strain. On the contrary, triaxial tests reveal

the deformability characteristics not only under hydrostatic pressures, but also under

deviatoric loading conditions (e.g., shearing under different constant confining pressures or

under constant mean stress conditions). The stress-strain relationships obtained in deviatoric

tests reveal the three-dimensional consolidation behavior. These tests furnish necessary data

for the formulation of the three dimensional constitutive equations (Cristescu [1994]; Gurson

& Yuan [1995]).

In this chapter, triaxial tests are performed to investigate the volume change behavior

of alumina powder under different stress states. The effects of particle size, initial density,

and time are studied. A series of triaxial tests on three types of alumina powders provided

by ALCOA (Aluminum Company Of American) are performed. The first is as-received A10

alumina powder with an initial volume fraction of 0.36; this type will be labeled as "dense

A10". The second one is also A10 alumina powder, but with a smaller prepared initial

volume fraction of 0.334, labeled as "loose A10". The third one is as-received A16-SG

alumina powder with a prepared initial volume fraction of 0.41. The specifications of these

materials are shown in Table 4.1. The elastic parameters were measured for all three types

of alumina powders based on the method proposed by Cristescu [1989] for time-dependent

materials. The triaxial test results show that the consolidation behavior of alumina powder

is strongly dependent on particle size, initial density, and time. In all cases a higher pressure

results in more consolidation. However, deviatoric stress can lead to either consolidation or

dilatancy. For the stress ranges applied to the A16-SG alumina powder, no dilatancy was

observed, while dilatancy occurred for both A10 alumina powders. The

compressibility/dilatancy boundary and the failure condition have also been measured. These

data can be used to formulate a general visco-elastic or elastic/viscoplastic model to describe

78

the consolidation behavior of these alumina powders. In the present research, a total of thirty-

three triaxial tests have been performed. These tests are listed in Table 4.2 and Table 4.3.

Most of the data are given in Appendices B, C and D.

Table 4.1 Specification of Alumina Powders

type of prepared initial density' Particle relative void ratio'

alumina (g/cm3) size(pm) density"

A10 (dense) 1.42 (two layer+vacuum) 40-200 0.36 1.78

A10 (loose) 1.32 (pluviation+vacuum) 40-200 0.334 1.99

A16-SG 1.62 (tap+vacuum) 0.5-1.0 0.41 1.44

SThe density after saturation.

"Theoretical Density of Alumina=3.95 g/cm3 (see Richerson [1982])

+ void ratio = V/V,,

where Vv is volume of void, Vs is volume of solid material in the specimen.

4.2 Experimental Procedures

4.2.1 Triaxial Equipment Setup

Fig. 4.1 shows the schematic diagram of triaxial equipment used (Bishop & Henkel

[1962]).

Table 4.2 List of triaxial tests for alumina powder (dense and loose A10)

No. date type initial confining "B" loading Number Number failure C/D comments

1995 of test density pressure value rate unloading unloading stress stress

-1996 (g/cm3) (kg/cm') (s") on hydro. on deviat. (kg/cm2) (kg/cm2)

1 2/14-17 CT 1.466 3=(7-4) 5.6e-5 no 1 9.28 5.2-5.6 Dry, try test

2 3/7-8 CT 1.388 3=(6-3) 0.93 5.6e-5 2 1 7.02 5-5.5 not max stress

3 3/13-16 CT 1.424 3=(6-3) 0.95 5.6e-5 1(8steps) 1 8.6 4.5-5

4 4/11-13 CT 1.424 3=(7-4) 0.97 l.le-5 4(7steps) 3 8.99 5-5.2

5 4/17-18 CT 1.424 2=(6-4) 0.96 1.1e-5 4(7steps) 6 6.02 2.75-3 unload too much

6' 4/19-20 CT 1.423 2=(6-4) 0.94 1.le-5 3(4steps) 5 5.99 3-3.2

7' 4/25-26 CT 1.423 1=(5-4) 0.97 1.le-5 1 2 3.11 1-1.2

8 5/1-2 CT 1.421 4=(7-3) 0.97 l.le-5 7(8steps) 5 10.93 8.3-8.5 not max stress

9" 5/3-8 CT 1.42 4=(5.5-1.5) 0.94 1.le-5 3 6 11.38 7.8-8 wrong on hydro

10' 5/10-11 CT 1.428 5=(6.5-1.5) 0.945 l.le-5 5 6 14.41 10.5-10.8

11' 6/27-29 CT 1.415 5=(6.5-1.5) 0.94 .le-5 5 5 13.79 10.6-10.9

12 7/3-5 CT 1.425 5=(6.5-1.5) 0.94 1.le-5 5 0 14.28 10.6-10.9

13' 7/18-19 CT 1.422 3=(4.5-1.5) 0.95 1.le-5 5 3 8.76 5.0-5.2

14' 8/2-3 CT 1.415 4=(5.5-1.5) 0.95 1.le-5 4 0 11.42 7.0-7.2

15" 8/8-9 CT 1.325 4=(5.5-1.5) 0.95 3.4e-5 4 2 9.58 9.3 pluviation test

Table 4.2 Continued

No. date type initial confining "B" loading Number Number failure C/D comments

1995 of test density pressure value rate unloading unloading stress stress

-1996 (g/cm3) (kg/cm2) (s") on hydro on deviat (kg/cm2) (kg/cm')

16" 8/16-17 CT 1.338 5=(6.5-1.5) 0.96 3.4e-5 5 3 12.23 12.1 pluviation test

17" 10/25-26 CT 1.316 3=(4.5-1.5) 0.97 3.4e-5 3 3 7.04 6.85 pluviation test

18" 10/30-11/1 CT 1.325 1=(2.5-1.5) 0.97 3.4e-5 1 1 2.46 2.1 pluviation test

19" 11/7-8 CT 1.313 2=(3.5-1.5) 0.96 3.4e-5 2 2 4.75 4.3-4.4 pluviation test

20 11/14-15 CT 1.412 3.5=(5-1.5) 0.95 l.le-5 4 2 9.64 6-6.1 chamber leak

21 12/12-13 CT 1.44 3.5=(5-1.5) 0.93 l.le-5 N/A N/A too large density

22" 12/19-20 CT 1.421 3.5=(5-1.5) 0.96 l.le-5 4 2 10.17 6-6.1

23' 12/27-28 RL 1.425 4=(5.5-1.5) 0.96 relax (2 steps) (8 steps) 10.8(min) 7.0-7.2 12.6 (max stress)

24 1/3-4 CT 1.411 4=(5.5-1.5) 0.97 5.4e-5 4 0 11.34 7.3-7.5

25' 1/8-9 CT 1.42 4=(5.5-1.5) 0.96 3.3e-4 4 0 11.77 8-8.1 very small dilat.

26' 1/11-12 CP 1.419 4=(5.5-1.5) 0.96 (1 step) 0 5.98 1.6 no effect on fail

Note: CT=conventional triaxial test (Kgnrmfn test); RL=relaxation test; CP=constant mean stress test

Confining pressure = (chamber pressure back pressure)

* The tests (total 11) have been listed in Appendix B

SThe tests (total 5) have been listed in Appendix C

Table 4.3 List of triaxial tests for Alumina powder A16-SG

No. date type initial confining-p "B" loading Number Number failure C/D comments

1996 of test density (cham.-back) value rate unloading unloading stress stress

(g/cm3) (kg/cm2) (s ") on hydro on deviat (kg/cm2) (kg/cm')

1' 1/17-18 CT 1.632 4=(5.5-1.5) 0.96 3.4e-5 4 I 11.88 N/A compress.-only

2* 2/16-18 CT 1.604 2=(3.5-1.5) 0.95 3.4e-5 2 2 5.65 N/A compress.-only

3* 3/6-7 CT 1.635 3=(4.5-1.5) 0.96 3.4e-5 3 3 8.79 N/A compress.-only

4' 3/11-12 CT 1.628 5=(6.5-1.5) 0.95 3.4e-5 5 2 14.83 N/A compress.-only

5' 3/13-14 CT 1.623 1=(2.5-1.5) 0.95 3.4e-5 1 1 2.97 N/A compress.-only

6* 3/26-27 CP 1.634 3=(4.5-1.5) 0.95 3.4e-5 (1 step) 0 4.34 N/A compress.-only

7 4/8-10 CP 1.646 4=(5.5-1.5) 0.98 3.4e-5 (1 step) 0 5.66 N/A compress.-only

Note: CT=conventional triaxial test (KirmAn test); CP=constant mean stress test

Confining pressure = (chamber pressure back pressure)

* the tests (total 6) have been listed in Appendix D

vacuum/backpressure

burette

air pressure dial

(Chamber pressure)

displacement

- measurement

Fig. 4.1. Schematic diagram of triaxial equipment setup.

The equipment consists of primarily a chamber, framework, panel, and base motor.

There are several valves to control the flow of water and/or air. Two inlets/outlets at the

bottom and the top of the specimen are connected to a vacuum or/and a burette on the panel.

The following quantities have been measured: axial force by load cell, axial displacement

of specimen by digital gauge, chamber pressure and backpressure by pressure transducers

respectively as shown in Fig. 4.1. Chamber pressure can provide an all-around pressure on

the top and the lateral surface of specimen. Back pressure acting as pore pressure can provide

the pressure of fluid in the specimen.

4.2.2 Specimen Preparation

The specimen is prepared inside a split mold. The split mold consists of two half

cylindrical shells which are connected with a vacuum pump. First, two clamps are put to hold

the mold. A membrane of about 0.12 mm thickness is put inside of the mold and the ends of

the membrane are tightened on the mold by several O-rings. A porous stone is placed at the

bottom of the mold. The mold is then put on the pedestal which connects all necessary fluid

inputs/outputs. A vacuum is applied to pull the membrane to the mold sides. There are two

methods to prepare the specimen in order to obtain a desired initial density. One is the so-

called layer method. The other one is called pluviation method. In the layer method, the

specimen is prepared in the successive layers. For each layer, a constant density is maintained

by certain height with a given amount of powder material. A small vibration is applied by

taping the side of the mold for densification. The maximum variation in relative density

between specimens can be kept within 0.5%. In the pluviation method, powder is poured into

the mold slowly without tapping the mold until the mold is full with powder. After the mold

is filled with powder, the top porous stone and cap are placed on the specimen. The vacuum

pump is disconnected from the mold and then applied to the top and bottom inlets of the

specimen. The mold is then removed. The specimen is supported by about 100 kPa vacuum

pressure. The average diameter Do of the specimen is obtained by measuring the diameter on

the top, middle, and bottom. The average initial height of the specimen, H0, can also be

obtained by measuring heights at four different points. By knowing the theoretical material

density y, of powder, weight of specimen W (powder is assumed dry), and the volume of the

specimen V, the volume fraction of the powder can be computed:

Y W

relative density=volume fraction (4.1)

Yt Vy, (4)

Afterwards the prepared specimen is put into the chamber.

A cylindrical specimen, about 70 mm in diameter and 145 mm in height, is used for

all tests. It is assumed that the specimens deform uniformly and the end effect is neglected.

Since the particle size is less than 0.2 mm, the membrane penetration effect is also neglected

(Frydman [1973]).

4.2.3 Back Pressure Saturation

The next step is to saturate the specimen. First, water is deaired in the burette for

about 20 minutes. Then water is allowed to flow into the bottom of specimen due to vacuum

pressure existing in the specimen. An alternating sequence of vacuuming the top of specimen

and allowing water in from the bottom is followed until most of the trapped air is removed.

A small confining pressure (03=30 kPa) is maintained during saturation to support the

specimen with a certain back pressure in the specimen. A "B" check is performed to verify

the saturation. Here, "B" is defined by B:=Au/Ao3, in which Ao, is the increment of chamber

pressure and Au is the increment of the pore water pressure in the specimen with undrained

conditions. During the tests, the volume of specimen changes. Saturation is essential as the

volume change is measured via the water volume entering/exiting the sample. Usually the

"B" value should be higher than 90% when the volumetric strain is to be measured under

drained conditions. In order to achieve a high accuracy of the measurement of volumetric

strain, the "B" value is required to be as large as possible. In our tests, municipal tap water

85

is used. In order to get higher B values, the prepared specimen is subjected to a backpressure

of 147 kPa over night. The B valve measured the next day is higher than or at least equal to

93% in all tests. See Table 4.2.

The change of the height of specimen, AH, during saturation is measured by attaching

the dial gage to the piston rod before saturation begins, and noting zero. Subsequently, the

small confining pressure (30 kPa) applied during saturation and the back pressure saturation

will cause volume reduction of the specimen. Since AD (the change of the diameter of

specimen) cannot be measured while AH can, it is assumed that AH / H0= -AD / D (axial

and radial strains are equal). The corrected section area A, after saturation is thus obtained,

Ho*,+2AH

A,=A o =A[1+2e] (4.2)

Ho

where e is axial strain after saturation.

4.2.4 Hydrostatic Loading Test

We begin with an incremental hydrostatic loading test (also called consolidation test

or isostatic test). Hydrostatic stress is computed as the chamber pressure minus the

backpressure. At each loading step, the stress is kept constant for several minutes and the

sample is allowed to shrink by creep until the rate of the volumetric strain is very small.

Usually it takes ten minutes. Afterwards the specimen is unloaded with small stress

increments and the stress is kept constant for 5 minutes. Then the specimen is reloaded to

reach the previous stress level. In this way accurate values of the elastic bulk modulus can

86

be calculated (Cristescu [1989]). From such tests, the volumetric change of specimen due to

both hydrostatic stress and, for a given stress level, time can be obtained. Since the

hydrostatic pressure causes specimen compression, the section area of the specimen is

corrected according to the following formula

H(V -Ac)

where: Vo is initial volume of specimen, AVc is volume change of the specimen during

consolidation at each step, AV,,=3VoAH / I H is the volume change of the specimen during

saturation, and AH is the change in height of the specimen at each loading step plus the

change during saturation.

4.2.5 Deviatoric Loading Test

The second stage of the triaxial test is the deviatoric loading test. There are two ways

to perform deviatoric tests. One is called "conventional triaxial test". In the deviatoric stage

of conventional triaxial tests, the confining pressure (lateral stress) is kept constant and only

the axial stress is increased. The confining pressure is the difference between chamber

pressure and back-pressure. The deviatoric stress is defined as the difference of stress

between the axial stress and the lateral stress. The other is called "constant mean stress test"

or "constant-p test". In the deviatoric stage of constant mean stress tests, the mean stress is

kept constant and only the deviatoric stress is increased.

87

Loading was performed at a constant displacement rate and the following quantities

in all deviatoric tests were measured: axial force, chamber pressure, back pressure, the length

change of specimen, and the volume change of specimen. The axial stress, confining

pressure, axial strain, and lateral strain are computed from the measured quantities. During

the test, the section area of cylindrical specimen changes. The following correction is applied

during tests:

V-AV I -E

A H l-A e ,l (4.4)

where VcH and Acare volume, height, and section area of specimen, respectively after the

hydrostatic test, AV is the volume change, AH is the change of the height during deviatoric

test, Ae, = AV/Vc, and e =AH /H. The deviatoric stress is defined as

P

O d (4.5)

where P is the axial load. Eqn (4.4) is also suitable for area correction for other type of tests

such as constant mean stress test or relaxation test.

4.3 Experimental Results

4.3.1 Hydrostatic Tests

Three series of hydrostatic tests (first stage of a triaxial test) on alumina powders

provided by ALCOA have been performed (see Table 4.2 and Appendices B, C and D). The

results of these hydrostatic tests are presented and discussed below.

The hydrostatic tests were performed with an unloading-reloading process in order

to obtain the elastic parameters. Fig. 4.2a shows a typical relationship between volumetric

strain and time in a hydrostatic test for dense A-10 alumina powder at different stress levels,

while Figs. 4.2b and 4.2c show similar results for loose A10 and A16-SG, respectively. At

each plateau the number shown marks the pressure applied. For each plateau shown, the

successive points correspond to 1.0, 2.0, 5.0, and 10.0 minutes (a total of about 10 minutes)

after each reloading. Afterwards an unloading increment of 20 kPa pressure is performed and

the remaining stress is kept constant for 5 minutes, followed by a reloading with the same

20 kPa, and then the final stress same as the previous one is kept constant for 5 minutes. The

cycle of unloading-reloading curve allows us to calculate the elastic bulk modulus. It was

found that A16-SG exhibits much more volumetric creep deformation than A10's for the

same testing time interval and same pressure, i.e., the A16-SG alumina is much more

compressible than the A 10's. After each loading, the volumetric strain continues to vary in

time (creep) as shown in Figs 4.2a-4.2c. The volumetric strain for A16-SG exceeds 6%,

while the volumetric strain for dense A10 and for loose A10 is only about 0.6 % and 1.1 %

under the same conditions (pressure 490 kPa), although A16-SG is initially more dense than

89

A 10's. This difference in compressibility is mainly due to the difference in particle size of

the powders. The particle size of A10 alumina is in the range of 40-200pm, while the

particle size of A16-SG is in the range of 0.4-lim. The particle shape of both alumina

powders is irregular. However, due to the different particle size, the surface area of particle

in these two powders is also different. The surface chemical properties of particle are

certainly involved, which will affect mechanical properties. This aspect should be

investigated separately.

As already mentioned, the unloading processes are used to obtain the elastic

parameters. The unloading processes also allow us to separate the time effect from elastic

properties. After each loading, the specimen is allowed to creep under constant pressure for

10 minutes. When the rate of volumetric strain becomes very small, the unloading is

performed. During the short time necessary to perform the unloading, the contribution of the

deformation by creep is negligible since then the material has reached nearly a quasi-static

state (no much strain can be obtained with the same stress afterwards as shown in Fig.4.3).

Thus, the influence of the time effect on the unloading process is negligible. The elastic bulk

modulus can be obtained from these unloading processes, as shown in Fig. 4.3, as the slope

of the unloading curve. The bulk modulus obtained by following this procedure is very close

to that obtained by using dynamic method for other materials such as rocks (Cristescu

[1989]). It is expected that for alumina powder the value obtained by the unloading processes

is also close to that obtained by dynamic method.

0.006

Sf 490

392

0.004

294

0.002 = 196 (kPa) dense AI0

98

0 I I

0 0.5 1 1.5 2 2.5

t(h)

Fig. 4.2a Volumetric strain versus time for dense A10 (points-reading data).

0.012

490

0.01

392

0.008

0.006

0.006 a =294(kPa)

0.004 196 loose A10

0.002 98

0 .... i, I. ..1 .liiiI

0 0.5 1 1.5 2 2.5 3

t (h)

Fig. 4.2b Same as 4.2a but for loose A10 (points-reading data).

0.07

0.06

490

0.05

392

0.04 -

294

0.03

0.02 = 196 (kPa)

0.01

0

98

0 0.5 1 1.5 2 2.5 3 3.5

t (h)

Fig. 4.2c Volumetric strain versus time for A16-SG (points-reading data).

It was found that elastic bulk moduli change smoothly with mean stress as shown in

Fig. 4.3 and Fig. 4.4. The elastic bulk moduli may also depend on other material properties

such as the density (Gurson & Yuan [1995]). As a first approximation, it is assumed that the

elastic bulk modulus depends on mean stress only. There are several possible functions

which could be used to fit the data if the behaviors of the elastic moduli at very high pressure

or/and very low pressure are also taken into account (Jin et al. [1996a]). Here, a linear

approximation is used to fit these data. It gives a reasonable fitting of the data in the pressure

interval shown in Fig. 4.4a and 4.4b. The bulk modulus K can be expressed as

K(o)= ko+ kjo (4.6)

where ko and k, are material constants. For comparison, the values of the coefficients k, and

92

k, for different types of alumina powders are listed in Table 4.4. It is seen that the denser

material has larger k,. However, the linear relationship (4.6) is not expected to be valid at

high pressure since K is expected to reach a limiting value when all pores are closed at very

high pressures.

500

(kPa)

400

300

200

dense A10

100

0

0 0.002 0.004 0.006

Fig. 4.3a Hydrostatic stress versus volumetric strain in hydrostatic

loading for dense A10 (points-recorded data).