NUMERICAL MODELING OF COMPACTION
OF PARTICULATE SYSTEMS
By
WEIJIAN WANG
A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY
OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR
THE DEGREE OF MASTER OF SCIENCE
UNIVERSITY OF FLORIDA
1999
ACKNOWLEDGMENTS
I would like to express my deepest gratitude to Dr. Nicolae D. Cristescu, my
advisor, for his guidance, tremendous help, and encouragement for the whole two years
passed. Dr. Bhavani V. Sankar, Dr. Peter G. Ifju and Dr. Edward K. Walsh are so kind to
serve on the committee. I am so grateful to their support and help. I also appreciate Dr.
Oana Cazacu for her guidance and discussions with me last year. Many other thanks
should go to my friends. I would have given up without their encouragement.
Finally, the financial support of the Engineering Research Center (ERC) at
University of Florida is gratefully acknowledged.
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS ...................................................... ii
ABSTRACT ...................................................................................... v
CHAPTERS
1 IN T R O D U C T IO N ................................................................................ 1
2 CLASSICAL MODELS FOR POWDER MATERIALS ............................ 4
2 .1 In tro d u ctio n ........................................................................... 4
2.2 Critical State Model .......... ..................................... 5
2.3 DruckerPrager/Cap M odel .................................. ............. 7
2 .4 F riction ................................................. ........... .. . ........ 10
2.5 Compaction Cases Study ................................ ......... 11
2.5.1 B background ......................... ......................... ......... 11
2.5.2 Computational Model ................................. ......... 12
2.5.3 FEM Model and Mesh .................... 14
2.5.4 Material Parameters ............................................. 14
2.6 R esult and A analysis .................................. .............. .............. 18
2 .7 C onclu sion .................................................. ........... ......... 27
3 A TIME DEPENDENT MODEL FOR ALUMINA POWDER CRISTESCU
M ODEL ........................................... ......... 29
3.1 Introduction .......... ...... .......................... ........... 29
3.2 Description of Cristescu Model .................................... .......... 30
3.2.1 General Form .......... .......................... ............. 30
3.2.2 M material P aram eters ...................................................... 31
3.2.3 Determination of the Yield Function H(c) ....................... 33
3.2.4 Determination of the Viscoplastic Strain Rate
Orientation Tensor N(c) ......................... 34
3.3 Simulation Algorithm for General Loading Test ........................... 34
3.3.1 Description ..................................... ............. ......... 34
3.3.2 Result and A analysis ........................................ ......... 37
3.4 Transient Creep Computation ................................................ 39
3.4.1 Description ..................................... ............. ......... 39
3.4.2 Result and Analysis ...................................................... 41
3.5 Analytical Solution for Transient Creep ........................................ 51
3.5.1 D description .......................................................... 51
3.5.2 Result and Analysis ..................................................... 52
3 .6 C o n clu sio n ........................... ........ ..... ............. ............. 54
4 REALIZATION OF CRISTESCU MODEL WITH FEM ............................ 55
4 .1 G e n e ra l ........... ... ..... ............................................. .......... 5 5
4.2 Algorithm for Stress Control Test ............................................ 55
4.3 Algorithm for Strain Control Test ....................................... 57
4.4 Algorithm for ABAQUS ........................ .......... 58
REFERENCES ........................................ ........... ......................... 62
BIO GRA PH ICAL SK ETCH ........................................................................ 64
Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Science
NUMERICAL MODELING OF COMPACTION OF PARTICULATE SYSTEMS
By
Weijian Wang
August 1999
Chairman: Graduate Research Professor Nicolae D. Cristescu
Major Department: Aerospace Engineering, Mechanics, and Engineering Science
The goal of this dissertation is to combine some constitutive models for particulate
systems with finite element methods(FEM) in order to understand the mechanical behavior
of cohesive powders during the compaction process and make some proper prediction for
manufacture control and tools design. One classical constitutive model the Drucker
Prager model is implemented with advanced FEM software ABAQUS. Also, a more
complicated model the Cristescu model will be introduced. This model will
hopefully be used with ABAQUS. Some algorithms will be introduced in this thesis.
CHAPTER 1
INTRODUCTION
Powder compaction is important in many aspects of industry, including the flow
and storage of agricultural products, the forming of ceramic components, powder
metallurgy parts, and pharmaceutical tablets, road and dam construction, subway and
tunnel excavation, mine and oilwell locating, and rock slides and avalanches. For
example, it is estimated that more than 80% of all medication doses are administered as
tablets, that is, in unit dosage forms prepared by compacting powders in dies.
The object of powder pressing is to quickly and reproducibly form a defectfree,
homogeneously dense powder to net shape. The pressing pressure should be high enough
to impart sufficient strength for subsequent handling but low enough to avoid excessive
wear on the press and tooling. Pharmaceutical tablets must have sufficient strength to
survive packaging and transportation but must be weak enough to disintegrate after
administration. Ceramic compacts may require higher strength to withstand the high
stresses during green compact machining.
In recent decades, a substantial amount of research and development has been
conducted in different areas and significant technological advances have been made.
Scientists have been trying to find appropriate models that can describe the characteristics
of different kinds of powders and predict the behaviors of the powder precisely during the
manufacturing process for the industry. It is of great significance because this will bring
more economic effects to the industry and engineers can control the particle compaction.
In modeling of powder compaction, constitutive relation and numerical methods
are developed for simulating the response of powder during compaction. A good model
should be able to predict the state of stress, strain and the displacements everywhere in a
powder sample subjected to some external and internal forces. Traditional soil mechanics
and rock mechanics concepts and theories are introduced and adopted into the constitutive
models for the powders. At the same time, finite element methods (FEM) have become
the most important numerical methods because they can accurately predict the density
variation in a compaction procedure and can provide information of defect detection.
More and more scientists are engaged in the work of combining classical or promising
constitutive models with FEM and compare with the results from the experiments and
improve them. This shows a great potential in the research work of powder materials.
Many theoretical models have been issued for understanding the constitutive
relations of powders. Among them the CamClay model21 and the DruckerPrager
model13 are undoubtedly the most popular and most widely used. They will be introduced
briefly in chapter two and some computation examples will be introduced at the same
time. These examples combine the DruckerPrager model with some advanced FEM
software (ABAQUS)14. They not only show interesting results but also tell us how
significant it is to implement the constitutive model into the modem numerical methods.
The advantages of the CamClay model and the DruckerPrager model are their easyto
use property and comparatively accurate agreement with experimental results under
various conditions. Unavoidably, the disadvantages limit the use of these models when
more experimental conditions are to be considered. For example, both of these two
models do not consider the time effects while it is very important in the manufacturing
process. This requires some more sophisticated models. Cristescu model is one of these.
Dr. Nicolae D. Cristescu has been working in the area of geomaterials and powder
materials for dozens of years.4,610,1618 His models for rock materials have been widely
used in scientific research and industry. One of his models for ceramic powder is
introduced in chapter three. Efforts will be made to implement this model with ABAQUS
software since ABAQUS provides some useful interfaces for users. Although no
satisfactory result can be presented in this dissertation, I would still like to introduce and
discuss the methods and algorithms briefly in chapter four.
CHAPTER 2
CLASSICAL MODELS FOR POWDER MATERIALS
2.1 Introduction
The powder behaves quite differently than the same material in solid phase and in a
more complex way, but under essentially monotonic loading conditions rather simple
constitutive models provide useful design information. The body composed of a large
amount of particles is thought of as a continuous media as are other solid materials, so
that the constitutive model of similar form to that used for other solid materials is also
used for powder materials.
The strain rate is decomposed into two parts, sE is the elastic strain rate and sF is
the irreversible (plastic) strain rate.
S= sE sI (2.1)
The elastic part is defined by a differential Hooke's Law, which is
F 1 1 (
sE [ +( ) l (2.2)
_J 2G 3K 2G 1(
in which K and G are elastic moduli.
For metals, sometimes the elastic part is enough for many applications. Some solid
materials accept a linear stressstrain relationship under a considerably large load. While
for powders the situation is quite different. The yield stress of the powder is often much
lower than that of metals. The postyield behavior cannot be described by the elastic
theory, and the constitutive equation for powder is of great concern for the industry. One
comparison between the implementation of elastic theory and plastic theory will be shown
in this chapter.
There are many kinds of plastic models that can be used to describe the plastic
behavior of the powder. We cannot use the perfect plasticity models because the powder
exhibits work hardening. Fortunately ABAQUS provides several options for nonmetal
materials. The plastic part is defined by the flow rule:
k g (2.3)
where g is a "plastic potential" and X is the flow parameter. (Mises, 1928)
2.2 Critical state model
The critical state model (Camclay model) was developed by Roscoe and Schofield
et al., 1968 and Parry, 1972, at Cambridge University. This model is used often for
modeling noncohesive materials. It contains mainly two parts: a yield surface, and a
critical state surface that is the locus of effective stress.
The equation of the yield surface is
f(p,q) (P 1)2 +( q 2 1= 0 (2.4)
p a Ma
The equation of the critical state line is
q=Mp (2.5)
p and q are the mean stress and the equivalent stress respectively. They are defined
1
P = (, + CY 2 +3)
S:S (2.6)
q= 2:S
S is the stress deviator tensor.
S = T +pI (2.7)
T is the stress tensor.
M and P are userspecified material parameters and a is the evolution parameter.
Figure 21 shows how the yield surface of the Camclay model looks like.
Critical State Line
" Yield surface
Figure 21. Typical yield surface for Camclay model
Camclay model is often used for soil and sand materials. It doesn't consider the
cohesion of the material. Unfortunately for most of the metallic powders and
geomaterials, cohesion cannot be neglected. Thus another model is introduced here, it is
the modified DruckerPrager/cap model.
2.3 DruckerPrager/Cap model
A typical yield surface of this model is shown in Figure 22. This surface is shown in
the qp space in which p and q are the mean stress and the equivalent stress and S is the
stress deviator tensor. Their definitions are the same as that defined in (2.6) and (2.7).
Cap, F,
R(d+P tanp3)
Figure 22. Typical yield surface for DruckerPrager/Cap model
Similar to the Camclay model, the yield surface in the modified DruckerPrager/cap
model is expressed by two segments: a shear failure surface and a cap that intersects the
equivalent pressure axis. There is also a transition region between these segments, which
provides a smooth connection. The cap surface serves two main purposes. It bounds the
" 
yield surface in pure hydrostatic compression. It also controls the volume dilatancy when
the material yields in shear.
The equation of the shear failure surface is written as,
q ptan3 d = 0 (2.8)
where P3 is the internal friction angle and d is the cohesion of the material.
The equation of the shear failure surface is written as,
/ ,2
(ppa)2 + R(d+pa tan3)= 0 (2.9)
l^ cosp]
where R and ca are material parameters that control the shape of the cap. Pa is an
evolution parameter and defined as
Pb Rd
Pa = b R (2.10)
1+Rtan p
The parameter ca is a small number, typically 0.01 ~ 0.05. It is used to provide a
smooth connection between the shear failure surface and the cap surface.
The cap surface has an elliptical shape with a constant eccentricity in the meridian q
p plane.
The function involved in the transition segment is
(Ppa)2+ q(1 )(d+patanp) u(d+patanp) =0 (2.11)
1 cosp
Plastic flow is nonassociated to the shear failure surface and the transition segment,
but is associated to the cap surface. As shown in Figure 23, the flow potential surface is
made up of two parts 
shear surface respectively.
The functions of Go
G, = (p Pa)2 +
Gs = [(pPa)tanf
Sllmilar
ellipses
Go and Gs. They are associated with the cap surface and the
and Gs are written as,
Rq )2
l+(a
cos P3 (2.12)
2 +( q )2
l+aC
cosB
G, (Shear failure)
_
 \ ..   \. I ap) (l4Uaa Seci
R(d*p.ian~'l
Figure 23. Flow surface for DruckerPrager/Cap model
While in the manufacturing of powder materials, the variation of the density is of
great concern and additional derivation about the density variation is given here.
We use the formula
v V (2.13)
Vo Vo
no" '
in which ev represents the volumetric strain, Vo and Vi are the initial and final
volumes respectively.
Then,
m m
8 = PI o = PO 1 (2.14)
m p1
Po
So,
p1 =Po (2.15)
1+Sv
Equation 2.16 will be used to compute the density variation during the compaction
described in the following examples.
2.4. Friction
Usually friction exists between the powder sample and the container wall. Thus
friction cannot be neglected in the computational model. The standard friction model is
the classical Coulomb friction, which is accepted and used in ABAQUS. [11][12] This model
states that at the contact surfaces the two bodies do not slide over each other as long as
the shear stress magnitude is less than P3, the friction coefficient times the pressure stress
between them. For alumina powder we assume P3 is a constant. As a matter of fact, we
do not consider the deformation of the container wall. So the container wall is always
modeled as a rigid surface in a FEM model by ABAQUS.
2.5 Compaction Cases study
2.5.1 Background
In the industry of manufacture or transport and handling, the shapes for products are
versatile and the containers for storing powders have various shapes. See Figure 24 and
Figure 25.
a. Solid cylinder
b. Cylindrical pipe
c. Tshape d. Cupshape
Figure 24. Different shape powder products
r to I
a. Pyramidshape container b. Coneshape container
Figure 25. Different shape containers
2.5.2 Computational model
In this chapter, two cases are used as typical samples for computation. See Figure
26 and Figure 27. In case I, the powder is surrounded by the die with a cylindrical
shape. Correspondingly, a cylindrical shaped powder sample forms after the compaction.
There are two punches acting on the top and bottom of the powder sample. One can
assert a downward load on the top punch or an upward load on the bottom punch.
The second case is a little more complicated. Different from the last case, there are
two punches at the bottom of the sample. Thus two forces can be loaded from the
bottom. Finally a cup shaped sample forms.
In the uniaxial test, the lateral surfaces of the sample are regarded as fixed because
of the rigid die wall of the container. In the computation, we have 2 =3 = 0 (in
Cartesian coordinates) or sr = e = 0 (for cylindrical coordinates).
13
P2
Punch
Powder
Punch
Pi
Figure 26. Case I cylindrical sample apparatus for compaction
Pi
7 Punch
Powder
Die Die
PPunch
SPunch
P3 P3
Figure 27. Case II Cup shape sample apparatus for compaction
2.5.3 FEM model and mesh
In both cases 8node rectangular element is used. Due to the symmetry, only half of
the cross section is considered. Figure28, Figure 29 and Figure 210 describe the model
size, boundary conditions and load conditions.
In the first case, there are altogether 49 elements. The left side is the axis of the
sample. Because of the symmetry, only the horizontal displacement is prescribed. The
bottom is assumed to be fixed. The load acts on the top. Two kinds of loads are
considered in this case. One is a uniform displacement load, that is, the top surface will
keep parallel to the horizontal plane. The other one is the uniform distributed load.
2.5.4 Material Parameters
Table 21 lists the parameters used in the computation. The values of these
parameters are obtained from the experiments.
Bulk modulus, K 6777.35MPa
Shear modulus, G 3531.86MPa
Elastic modulus, E 9027.43MPa
Poisson's ratio, v 0.278
Internal friction angle, f3 16.50
Cohesion, d 5.5Mpa
Cap shape parameter, R 0.558
Transition surface parameter, ac 0.03
Table 21 Table of Material parameters
1224
1624
,or
2C24
" 1424 1824 2224
1222 1622 2C22
1422 1822 2222
#1^, 1220 1620 2C20
14:0 1820 2220
89) r (B1@ (B (B@ B
1018
, 1218 1618 2018
1418 18 8 2218
116 1216 1616 2 16
 1416 18 6 2216
S114 1614 2(14
1414 1814 2:14
#^, 1212 1612 2 12
1412 18 2 2212
1010l 1410 1810 2210
0.01 m 
Figure 29. FEM mesh and boundary conditions for cylindrical shape sample
0.01 m
2424
2422
2420
2418
2416
2414
2412
0.385 m
2
410
@ 0 @ @@@@
0 @ @@@@
@ @ @ @ @@(@
@ @ @ @ @@@@
1 3 5 7 9 11
13 
15 
19 
8611
8411
8211
8011
 H 8013
 H 8015
4811
8017
 8021
Figure 210. FEM mesh for cup shape sample
801 803 805
807
1 18 48 1 1 88 1 1
i  8019
I~p1' 8 8j.1~


1
21 8021
Figure 211. Boundary conditions for cup shape sample
801I
8811,
The cap hardening curve shown in Figure 211 is also from experimental data. It
determines the corresponding parameter option in the ABAQUS input file.
0 05 0 1 0 15 02 0 25 03 0 35
plastic volumetric strain
Figure 211. Hardening curve
2.6 Result and analysis
2.6.1 Case 1.1
Cylindrical shape sample, distributated load acting on the top
This is a load control case. Note that the top surface is not flat any more after the
compaction. Due to the friction existing between the sample and the wall, the
displacement at the edge is smaller than at the center. Thus the initial horizontal lines form
into the curves with the concave downward shape, as shown in Figure 212.
Figure 213, Figure 214 and Figure 215 show the distribution of the axial stress,
the volumetric strain and the density after the compaction respectively. All of them
distribute stratiformly because of the existence of friction. If there is no friction, their
distribution should be uniform over the crosssection.
19
The volumetric strain and the density actually show the same result because they
differ only by a constant (see equation 2.11). Density is of great concern in the industry
because people want to know how the powder mass is distributed after the compaction,
and this distribution determines the quality and the characteristic of the products.
LFigure 212. Deformation after loading
Figure 212. Deformation after loading
T. 49L01
 97L01
6 45F1>o
.5 DIXot
S 42k01
4, DOz~O1
A 311,01
3 671.Ot
3 35LOt
 2. BJE+o1
3 1z.Oi
 9 OIO1
I31 O
Figure 213. Axial stress distribution after loading
I
LI
DfBRK1 f lLVE
 3. 3; i01
 7. 1zO
3. 11101
B3 2I01
2 Mlr01
2 r5101
J 96L01i
1. 93iOi
] TO 01
Mu,
Figure 214. Volumetric strain distribution after loading
i. AE, QD
91. BE OD
. ,IE UD
2 2DE D0
S24E, 00D
,2 DE 3 OD
3E tE UD
imur
Figure 215. Density distribution after loading
2.6.2 Case 1.2 Cylindrical shape sample, uniform displacement on the top
Figure 216, Figure 217, Figure 218 and Figure 219 show the results under the
uniform displacement load. They are very similar to the results under the uniform
distributed load. The volumetric strain and the density distribute more evenly than in the
former case.
L1
Figure 216. Deformation after loading
~a VALVE
.4 *~5*qI
4 ~4E.!l
4 ~1L.!1
~. OJt' II
.k ....I
: T~t'II
4 4E'*1
i h.t
1 1~54*i
L1
Figure 217. Axial stress distribution after loading
TraHl rJitli
.9 aoE1
s i .
5, OT 1 I 1
4m
i *;.i I
t j11 .11
Figure 218. Volumetric strain distribution after loading
ti M00
tZ, iiEtDO
+2, 222E4O
+2.3 +00
+2. 56O +00
*2. 56i+SO
+2. 67 +00
+2. 1t OQ
3, DE+tDO0
t6, 5EtaO
Figure 219. Density distribution after loading
The difference between the elastic model and plastic model is apparent, as shown
in Figure 220, especially for powder materials. Thus, only the elastic model is not
enough for the study of the powders and geomaterials.
60
50
40
30
20
10
0
0.00
0.20 0.25 0.30 0.35
Figure 220. Comparison of the results from the elastic and the plastic model
Friction does not play a big role in this case. In the elastic stage, there is almost no
effect from friction. However, only the basic Coulomb friction is implemented in this
chapter. If another more complicated friction model is used, we may see how it works in
more detail.
80
without friction
with friction
0.05 0.10 0.15
axial strain
0.20 0.25 0.30 0.35
Figure 221. Comparison of the cases with and without friction
0.05 0.10 0.15
Strain
70
60
50
40
J 30
" 20
10
t<
0 
0.00
Figure 222 and Figure 223 show the density variation of such a point lying at the
center of the half cross section under two different load conditions. They exhibit the same
tendency except at the two ends, that is, at the beginning and at the end of loading that
could cause some unstable results.
2.50
2.00
r
E 1.50
1.00
0.50
0.00
0.000 0.001 0.002 0.003
displacement (m)
0.004 0.005 0.006
Figure 222. Under displacement control, the density variation at point 1818
2.50
2.00
o
1.50
0.50
0.00
0 10 20 30
stress (MPa)
40 50 60
Figure 223. Under stress control, the density variation at point 1818
Figure 224 shows the general trend of the stressstrain curve. It conforms to the
tendency exhibited in the hardening curve of this material. (See Figure 211)
60
50
40
30
20
10
0
0.00 0.05 0.10 0.15 0.20 0.25 0.30
axis strain
Figure 224. Stressstrain curve at point 1818 under displacement control
2.6.3 Case 2 Cup shape sample, uniform displacement on the top
It is available to handle those samples with a little more sophisticated shapes using
DruckerPrager/Cap model. Figure 225, Figure 226 and Figure 227 show the
distribution of the axial stress, the volumetric strain and the density respectively. Tran et
al did the computations for the models with cylinder shape and Tshape.161 This is
interesting and challenging because either the industrial products or containers have
versatile shapes, sometimes some of them are even weired. Sophisticated shape brings
difficulty to the modeling work. One needs to consider the complicated boundary
conditions and some negative effects caused by the complicated geometry such as stress
concentration. In the example above, we should be aware that around the interior corner
of the cup singularities may exist. We need to use more advanced theory if we want to
explain some uncommon pheonomena. Here we consider only the horizontal and vertical
section. Along the vertical section the density distributes in layers while along the
horizontal section the density is almost homogeneous everywhere. This is reasonable and
was proved by the experiments.[61
4 96tOoQ
4. 02 E+00
). 08E+00
2 .L4+DD
J. 2.E+00
L
+1. L E+00
Figure 225. Axial stress distribution after loading
1" 3LV
Figure 226. Volumetric distribution after loading
U VA.h2 VALUB EI
iti SsBOs0
t],66EtO
t1], U EtDO
4], I Et00
4 1, l 6 E40i0
.a, 78s0Q
1. 00 E+00
J] .2Et0 
L
Figure 227. Density distribution after loading
2.7 Conclusion
DruckerPrager/cap model can be used to describe the material behaviors during
loading and unloading. It can be used for the models with different shapes shown in
Figure 24 and Figure 25 or even more sophisticated geometry shapes. Some
experimental data show excellent agreement with the computational results obtained using
this model. However, the plots of the contours of the density show that the distribution of
the density is not uniform. Only the density values in the middle are close to the
experimental data. At the same time the DruckerPrager/cap model is relatively simple, it
assumes that the material parameters, such as the elastic modulus and friction angle, are
constant. As a matter of fact they are not invariable throughout the loading and unloading
procedure. And it doesn't take into account the time effects. We know that for powder
28
materials the time effect is very important. Also, it is believed that it is not appropriate for
some kinds of powder like consolidated clays and loose sands. Therefore more advanced
models are expected to be used to simulate more complicated experimental conditions and
processes.
CHAPTER 3
A TIME DEPENDENT MODEL FOR ALUMINA POWDER CRISTESCU MODEL
3.1 Introduction
For the geomaterials and metallic powders, it is necessary to incorporate the
irreversible volumetric changes and time effects into the constitutive equation. The
uniaxial compression tests are revealing that the behavior of the powder material is always
nonlinear. The volume is first compressible and afterwards dilatant. Just before failure the
volume increases very much. Most of the volume change is irreversible and should be
considered into the constitutive model. Since before failure the volume is dilatant, it is a
property that may be used to predict failure. Time effect is another important factor. It is
apparent in a uniaxial test that the faster the loading rate is, the higher the stressstrain
curve is. That is, for the same deformation different stress levels can be reached.
During the decades of 1980s and 1990s, Dr.N. Cristescu and his colleagues studied
extensively various geomaterials.4'610'1618 They found, first, that the irreversible
volumetric changes must be introduced into the constitutive equations. Secondly, in a
uniaxial compression test the volumetric deformation is nonlinear. The volume shows
compressibility at the beginning and dilatancy afterwards. Most of the volume change is
irreversible. The strain rate or loading history has a significant influence on dilatancy,
failure and in a lesser degree on compressibility. Secondly, time effects are also very
important. When loading rate is increased the whole stressstrain curve raises starting
from the lowervalue stress and strain. Thirdly, they introduced a very important concept
 the compressibility/dilatancy boundary. Most powder materials are compressible for
small octahedral shear stress, but dilatant for higher values. This boundary indicates
where the material is compressible and where it is dilatant.7'9
3.2 Description of Cristescu Model
To completely describe the constitutive relations for an elastoplastic material, four
different assumptions are required:
1. A yield function, which generalizes the concept of the yield stress to 2D and 3D
stress status.
2. A relationship between the directions of the principal plastic strain increments and
the principal stresses.
3. A flow rule that specifies the relative magnitudes of the incremental plastic strains
when the material is yielding.
4. A hardening law, a relationship between the equivalent stress and the equivalent
plastic strain or the stress works that is done on the material when it is yielding.
Let us see how Cristescu model works.
3.2.1 General Form
The general form of Cristescu model9 is
1 1 1 W(t) F (3.1)
+( )ol+k 1 (3.1)
2G 3K 2G H(cy) /Bc
in which,
s is the total strain rate.
G and K are elastic moduli.
k is the viscosity parameter.
=1/2(A+A).
H(c) is the yield function.
F(Cy,T) is the flow potential function. For alumina powder, F can be replaced with a tensor
evaluated function N(c) which defines the orientation of the viscoplastic strain rate.4
3.2.2 Material Parameters
In this model, the material parameters such as elastic modulus are not constant any
more. Some fitting formulas are obtained from the hydrostatic and triaxial experimental
data.
K(cy) = K paa exp b b
L Pab (3.2)
E(c) = E pa Pexp(dG)
where K", E" are the constant limiting values that K and E approach as the pressure
increases. c, P3, b, and d are dimensionless material constants.
Figure 31 and Figure 32 show how K and E vary with the mean stress. K and E
are both asymptotically approaching to some constant values when the stress level
increases.
11000
10000 
9000 
8000
0 7000 
6000
0 5000
= 4000
3000 
2000 
1000 
0
0 20 40 60 80 100
mean stress(MPa)
Figure 31. Bulk modulus variation
800
700
600
2 500 
400
S300
200
100
0
0 100 200 300 400 500 600
mean stress(MPa)
Figure 32 Elastic modulus variation
3.2.3 Determination of the Yield Function H(y)
For the isotropic materials the yield function H(y) depends on stress invariants only.
It can be determined from the data obtained in drained triaxial compression tests in which
the first stage is hydrostatic and the second is deviatoric. Thus we assume that H(c) has
two parts:
H(c) = Hh (G) + Hd (Y, T) (3.3)
In the hydrostatic stage, Hd(o, c) = 0.
Hh ()= ah ()2 +bh () (3.4)
Pa Pa
CY
Hd (,T) = a( +a2)m()9 +B( )2 (3.5)
Pa Pz Pa
where ah, bh, a1, a2, m and B are fitting constants.
3.2.4 Determination of the Viscoplastic Strain Rate Orientation Tensor N(G)
N(c) is a tensor function to describe the orientation of the irreversible strain
increments.4 For isotropic material, N(c) is assumed of the form
N(y) = N1I+N2 (3.6)
where cy' is the stress deviator.
Thus,
; =k I (NI+N2 ) (3.7)
H )T
kNi and kN2 can be calculated by the following formulas,
kN, = q(C) + sgn X(o, Tc)Y(c, c)
Pa
cy a (3.8)
kN2 =m(_)14 m2( ) +n2(J)2exp n2( +r
Pa Pa Pa Pa
in which
p(c) = ak C2 + bkC
X(c,'c) = 0.555 c (3.9)
5
T T
S2Y(oT)= z y )y exp(zb, )
WQ(,'c)=Zc 2) 2c 
Pa Pa Pa
where ak, bk, Zc, ya, yc, Zb are material constants.
3.3 Simulation Algorithm for General Loading Test
3.3.1 Description
In a general loading test, one can get the stressstrain curves from different test
conditions such as loading time, loading speed, ambient temperature. The influences of
these conditions are important for evaluating the mechanical behaviors of the material.
Here, the test procedure simulated with the algorithms introduced in this chapter always
starts from a hydrostatic loading. For example, if a hydrostatic loading is carried out up to
C71 = 02 = C3 = 294kPa, afterwards only the axial load is increased while the other two
stress components are kept constant. There are three kinds of algorithms for three
different types of test, the universal continuous loading test, and the creep test. There are
two kinds of creep tests, one considers the transient stage only and the test stops as the
stabilization is reached. In another algorithm a fixed time length is considered at each load
step. In the general test, we don't care about the particular load path. After reaching the
yielding, the irreversible strain rate is
S(1 W)(kNI + kN2 ) (3.10)
H T
Next, we make the inner product with stress tensor a on both sides. Recall the
definition of W, W(t) = f : s'dr On the left side,
'to
W = (1 W)(kN I + kN2 ) (3.11)
H T
in which,
l l = 3C (3.12)
= j (cj C61) = cj 3C2 (3.13)
then,
dW = [3kNkNN +C( 3Z2 3)]dt (3.14)
H
Integrate it in the time interval [tn, tn+i],
W'n+ ft [3kN, N2 ( 3C2)]dt (3.15)
1
H
H ln(1 ) Wn+1 = tn [3kNC + (kN2 3C2)]dt (3.16)
Let FI = 3kNY + kN2 ( o 3C2)
1 Wn+
HW,
SW.
exp( 1n FIdt)
H tn
(3.17)
(3.18)
Wn,+ = H+, [1(1 ) exp( Fdt)]
n+l H t
This is the iteration scheme for calculating W at every step. Numerically, the
integration in the formula above can be calculated using the trapezoidal rule.
jtn1 F~dt
(3.19)
 (FI + FIn +)
2
Table 21 lists the parameter values used in the computation.
Table 21 Table of material parameters
Confining pressure 294 kPa
Axial stress increment 20 kPa
Time increment 0.01 second
ah 5.833x106
bh 1. Ox 106
a, 0.27
a2 37
m 8.5
B 4.6x105
ze 1.127x1015
ya 0.842
ye 56.06
Zb 0.01
ak 8.0x1011
bk 1.641 x 106
ml 2.38x1033
m2 0.038
ni 4.0x107
n2 0.005
r 8.0x 104
P3 6.95x105
d 0.002
E __7.0x 105
3.3.2 Result and Analysis
Confining pressure is an important factor in the triaxial test. As shown in Figure 3
3, if the confining pressure is higher, then higher stress level is needed to get to the same
deformation. At the same time, more volumetric change can be obtained under the higher
confining pressure.
392kPa 700 (kPa)
392kPa
200
100
0
0.01 0.008 0.006 0.004 0.002 0 0.002 0.004
Figure 33. Stressstrain curves for different confining pressures in a general test
Figure 34 and Figure 35 show the influences of some factors, load increment and
time increment. Actually, these two figures express the same thing, the role of the loading
rate Ao/At in the general test. For any time dependent constitutive model like this, we
should show clearly how the time effects are described. In the triaxial test the loading rate
is often adjusted to different value levels as required in order to get different results.
Ac=20kPa Ac=30kPa
o3=294kPa
200 400 600
S(kPa)
Figure 34 Irreversible stress work for different load increments, same time increments
5
At=0.01s At=0.001s
a3=294kPa
200
t (kPa)
Figure 35.Irreversible stress work for different time increments, same load increments
2500
2000
1500
1000
500
0
Viscosity is another important factor. Figure36 exhibits its influence on stress work
W. If the viscosity is higher, which means the powder is more viscous, more difficult to
"flow", then more stress work should be done under the same stress level.
100 200 300 400
(kPa)
Figure 36.Irreversible stress work under different values of k
3.4 Transient Creep Computation
3.4.1 Description
Creep test is one of the most important tests to study the constitutive relation of
the materials, especially for the time dependent models. In the creep test, the load (stress)
level is kept constant at each step while the deformation still continues. There may be
many loading steps in a creep test. The mechanism of creep is discussed in detail in the
books by Boresi, Malvern and Cristescu.3'10'19
In this section, only the transient stage is considered. The calculation stops when
the deformation becomes stable. When the deformation approaches stabilization, the
irreversible strain rate approaches to zero. This rule is used as a condition to know about
the current state of the material. In order to compute the irreversible strain rate, the
irreversible stress work is required. The difficulty in the algorithm for the transient creep
is the iteration scheme the irreversible strain rate and the irreversible stress work.
A brief computing flow chart is introduced in this section. See Figure 37.
In a deviatoric test, the confining pressure keeps constant in the whole test
procedure while the hydrostatic stress and the deviatoric stress keep changing. We need to
know the influence from the hydrostatic stress and the deviatoric stress separately. Thus,
the irreversible stress work can be decomposed into two parts, the volumetric part and the
deviatoric part.
W(t) = J; dT
0
= t (76 1+ )(1' 6 + E, )d' (3.20)
to j d
Let o (3.21)
Wd =J t dT
Then,
W = Wv+ Wd (3.22)
41
Calculate the stress at the
current step
Calculate a, a', and c for the
current step
Calculate material parameters
for the current step
E(c), K(c), G(c)
Calculate H and W for the
current step
W N
(s =( Nextl)( ad tep
SCalculate totalversible tra instrain
4Yes current step
w = ft a:FrNext load stepd
Figure 37 The algorithm for the transient stage in the creep test
3.4.2 Result and analysis
In this test, after each load increment, the stress will be kept constant until the
material reaches the stabilization. The time to get to the stabilization is different at
different stress levels. Finally, if we connect the successive steps, we can not get a linear
response, as Figure 38 shown.
500
450 
400 03 = 294 kPa
350 Acy = 20 kPa
300 
(kPa)
250 
200 
150 
100 
50
0 10 20 30 40 50 60 70 80 90
Time(s
Figure 38 Octahedral shear stress in triaxial test
Figure 39 is the typical stressstrain curve that can be seen in many experimental
results. Axial strain si and lateral strain 82 expand in different directions. Volumetric
strain Ev is also presented in this figure. It shows clearly the transition between the
compressibility and dilatancy.
Figure 310, Figure 311 and Figure 312 show how the irreversible stress work
and its components as they vary with the increasing load. Although the total irreversible
stress work increases monotonically as the load increases, its volumetric part shows a
complete different tendency. This illustrates the concept of compressibility and dilatancy.
When the volume of the sample begins to increase, the sign of the volumetric irreversible
stress work becomes negative.
than that of the deviatoric part.
500 
450
400
350
300
250
r (kPa)
200 
150
100
50
0
0.1
However, the magnitude of volumetric part is much less
0.05 0 0.05
Strain
0.1 0.15 0.2 0.25
Figure 39 Stressstrain curves
4
2
0
2
Wv (kPa) 4 
6
8
10 
0
100 200 300
r (kPa)
400 500
Figure 310 Volumetric irreversible stress work as a function of c
44
150
100 
Wd (kPa)
50
0
0  
0 100 200 300 400 500
c (kPa)
Figure 311 Deviatoric irreversible stress work as a function of c
140
120 
100 
80 
W (kPa)
60 
40 
20
0
0 100 200 300 400 500
c (kPa)
Figure 312 Total irreversible stress work as a function of c
In Figure 313, we observe clearly the transition from the compressibility to
dilatancy. Also, the confining pressure is one factor that can influence the result. The
higher the confining pressure is, the higher the load is needed to get to the same
deformation. At the compressibility/dilatancy transition points that locate on the rightmost
of the curves respectively, the volumetric change is larger if the confining pressure is
higher.
Figure 314 shows the deviatoric irreversible stress work for the different confining
pressures. The deviatoric part dominates the magnitude of the irreversible stress work.
As already mentioned, the higher the confining pressure is, the higher load level is needed
to get the same stress work. Moreover, under the same stress level, the results of the
deviatoric part of the irreversible stress work are much different for different confining
pressures.
392kPa 600(kPa)
294kPa 500
400
196kPa
200
100
0.01 0.008 0.006 0.004 0.002 0 0.002 0.004 0.006 0.008
Figure 313 Influence of confining pressure on the volumetric strain
In Figure 315, we can clearly observe that the volume firstly decreases, which
means the compressibility. Afterwards it increases, which means dilatancy.
301
25
20
15
10
5
0
98kPa
196kPa
294kPa 392kPa
100 200 300 400
S(kPa)
Figure 314 Variation of the irreversible stress work as a function of C for different
confining pressure
x 103
0 10 20 30 40 50 60 70 80 90
Time(s)
Figure 315 Volumetric strain during the loading procedure
Figure 316 shows the results under different load increments. The confining
pressure is kept constant (294kPa) while the axial stress increment takes three different
values from low to high, lOkPa, 20kPa, and 30kPa. If the load increment is lower, the
value of ev is larger at the C/D transition point, which means higher compressibility is
obtained. At the same time, the time to get to the C/D boundary is longer if the load
increment is lower. The constitutive relation here is only related to the load increment,
instead of the loading rate. To make this point more straightforward, let us see what will
happen if the load changes dramatically, see Figure 317.
X 103
4
3.5 
3
2.5
ev
2
1.5
1
0.5
0
20 40 60 80
Time (s)
100 120
Figure 316 Volumetric strain for different loading increment and same confining pressure
The curves in Figure 317 exhibit the same tendency as that shown in Figure 316.
When the load increment is larger, the number of the step is fewer.
X 103
3 i
h
2 CD e f
1 B d
A c E
a
1 failure
2
33 = 294 kPa
4
5
failur
6
0 5 10 15 20 25
Time (s)
Figure 317 Volumetric strain for larger loading increment, same confining pressure
The viscousity parameter k may vary due to some reasons(temperature, etc). As
shown in Figure 318, if k is smaller, the material is more flowable and longer time is
needed to get to the stabilization.
Figure 319 shows the Compressibility/Dilatancy boundary, failure surface, the
curve OH/Wc = 0 and several yield surfaces.
The equation of the Compressibility/Dilatancy boundary is
0.555 c = 0 (3.23)
The equation of the failure surface is
ar + b, c = 0 (3.24)
For this kind of alumina powder, ar = 0.679, br = 9.308.
OH
From = 0, we can derive out
T
2ah bh am 2 m F2 90
( ) + ( + a2) () = 0 (3.25)
Pa Pa Pa Pa Pa Pa
As a matter of fact, this curve passes in the 2D plane through the points of the
maxima on the yield surfaces H = constant.
X 103
3.5 
3
k=10.0 k=.0
2.5
Sv 2
1.5 
1
failurefailure
0.5 
0
0 10 20 30 40 50 60 70 80
Time (s)
Figure 318. Volumetric strain curve for different value of k
The points shown in Figure 320 correspond to those in Figure 317. It shows the
load path in the whole test procedure, from where the deviatoric test begins, how it passes
through different yield surfaces, the C/D boundary and reach the failure surface. The
small letters correspond to the case with the axial stress increment of Aci=98kPa while the
capital letters correspond to the case with Acy=196kPa. We can see that these two stress
states coincide every second point. Connecting these points we can get a line which has a
slope of 1.4142. It is just the equation of the deviatoric loading path, which is
S = const (3.26)
2
Expand (3.26),
( 2= ( +2 2)2 2(Y 2) =2 (3.27)
Equation (3.27) means the confining pressure doesn't change along this loading
path, that is, in the deviatoric test, although cy and c keep increasing. It is correct for
Karman test in which 02 = 3.
800 C/D bo ndary
OH/y0=0
600 Failure
400
H=8.03kPa H=19.22kP
200
=1.88kPa
0
0 500 1000 1500 2000
a (kPa)
Figure 319 Compressibility/Dilatancy boundary, failure surface 8H/8y = 0, and several
yield surfaces
800
Failure C/D boundary
600
k F
& 400 1
200 e0C Arctan(1.4142
c, B
b
0 A
0 500 1000 1500 2000
cy (kPa)
Figure 320 Deviatoric test procedure in Tc plane
3.5 Analytical Solution for the Creep Test with Fixed Time Interval
3.5.1 Description
In this type of creep test, we assume that at time to a load starts. The stresses
increase quite fast from the previous stress level and afterwards they are kept constant for
a fixed time interval. More loading steps like this repeat with the increase of the stress
level if required. A new formula can be obtained.6
1 W =(t) 1 ( exp (kN 1 kN2 ).(to t) (3.28)
H(y,T) [ H(cy, ) 1 1 TI
where o(t) = o(to) = constant and W(to) is the initial value of W at time to. If the
initial state is on a yield surface then W(to) = H(co). The formula for calculating strain at
arbitrary time is
I L W(to kN + kN y(t0
s(t)= sE + exp (kN,1+kN, (tot)
l(kN + kN2 ). L
H T
3.29)
with the initial conditions
= 0, W(t0) Wo
I 3K 2G 2G
This time the time interval during which the load at each step is constant. It can be
set up at different magnitudes.
3.5.2 Result and Analysis
Figure 321 exhibits the same story as Figure 317 does. For bigger value of k, the
time is enough for the material to reach the stabilization while it is not enough for smaller
k. Thus we see two different curves, one with the apparent creep stage while the other
only with the transient stage.
Figure 322 shows the results for different step time intervals. If the time interval
is relatively small, we can get the transient creep only.
X 103 ..
4
3
2 0k=.1
1 
Sv 0
1
At = 10 sec
2 
3
failur
4
0 10 20 30 40 50 60 70 80
time
Figure 321 Volumtric strain change for different k values
x 103
4
3 At=3sec
At=10sec
2
k=1.0
0
Sv
0 10 20
40 50 60 70 80 90
time
Figure 322 Volumetric strain change for different step time intervals
54
3.6 Conclusion
Cristescu model is in reasonable qualitative agreement with the experiment
results,4'6'10 although a quantitative comparison is not shown here. It describes the time
effects in the compaction and dilatant process. A more complete discussion could be done
if it could be combined with FEM methods. This will be introduced in next chapter.
CHAPTER 4
REALIZATION OF CRISTESCU MODEL WITH FEM
4.1 General
Cristescu model is a timedependent semilinear model. Usually, the variable
stiffness method is the most acceptable method to solve the nonlinear problem.24 In this
method an effective iteration scheme needs to be erected to describe the constitutive
relation. Any computation process should proceed in an incremental manner considering
suitable small time intervals. In each time interval the stress level should be known so that
the increments of strain can be found.
Owen and Hinton present the formulation of elastic/viscoplasticity theory in
discrete time form in the book. They provide the detailed algorithm and an example
program for solving some simple problems.20 Jin presented his algorithm in his dissertation
for alumina powder combining with a large amount of experimental work in his research.17
These are all good reference for designing the FEM algorithm for powder materials.
4.2 Algorithm for Stress Control Test
In stress control test, the loading rate or the stress increment is known. The
unknown is the strain rate or the strain increment.
Let us consider an arbitrary load step and define it as the n+lth step between the
state n and state n+1. According to the basic assumption,
A+= A n+ As n+(4.1)
The elastic part can be calculated using Hooke's Law,
Asn+ (D n+) 1A~n+1 (4.2)
in which Dn+1 is the elastic matrix at the n+lth step. In Cristescu model, the
elastic moduli are not constant. D is the function of 7.
Implement the implicit Euler method to get the irreversible strain increment,
AsI = At (l 0)sn+s0n+ (4.3)
0 is a constant and 0 < 0 < 1. To get the stable result, 0 should satisfy 0 < 0 < 1.
Plug (4.2) and (4.3) into (4.1),
n+1 Atn+i = D 1A + At(1 0)sn+s (4.4)
Equation (4.4) forms the iteration scheme of two successive steps.
Define two state variables p and q.
I
qn = D'Ac + At(16)sn+sn (4.5)
I
pn = n At0sn (4.6)
Then equation (4.4) changes into
Pn+l = qn (4.7)
If the state variables at the nth step are given, that is, tn,sn,Cn,sn,Sn,pn are
known.
Conduct another loop from the nth step to the n+lth step. Now consider the ith step
in this loop.
q = p, + P 6s, (4.8)
From equation (4.6),
S= 1 (4.9)
So, equation (4.8) changes into,
6s, = qn p, (4.10)
Update
1+1 = + 6 ,
(4.11)
P +1 = 81 AtOg1+1
The loops stops when 6Ei < tolerance and then continue the n+2th step.
4.3 Algorithm for Strain Control Test
In strain control test, the strain rate or the strain increment is known. The
unknown is the loading rate or the stress increment.
Let us consider an arbitrary load step and define it as the n+lth step between the
state n and state n+1. From equations (4.1), (4.2), (4.3)
D lo +1 + AtOn+ = A At(1 0)+D l c (4.12)
Equation (4.12) forms the iteration scheme of two successive steps.
Define two state variables p and q.
I
q. = As At(1 0)n + D 'A (4.13)
I
p. = D c +AtO n (4.14)
Then equation (4.12) changes into
Pn.+ = qn (4.15)
If the state variables at the nth step are given, that is, tn,sn,CY,sn,Sn,pn are
known.
Conduct another loop from the nth step to the n+lth step. Now consider the ith step
in this loop.
qn = p, + O o, (4.16)
From equation (4.14),
p D 1 + AtO (4.17)
So, equation (4.16) changes into,
P' = q n P, (4.18)
Update
(4.19)
Pl+1 = ,1+1 AtOsg+,
The loops stops when 60i < tolerance and then continue the n+2th step.
4.4 Algorithm for ABAQUS
Since ABAQUS is hopefully to be implemented with Cristescu model, we need to
pay more attention to the characteristics provided by ABAQUS. Besides those internal
constitutive models, ABAQUS provides user subroutine as the interface with the users.
User subroutine is a powerful and flexible tool for analysis because users can define the
constitutive relations by themselves. UMAT is one of these subroutines that could be
used. First, let us look at the interface defined in this subroutine.
*USER SUBROUTINE
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,
2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,
3 NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,
4 PNEWDT,CELENT,DFGRDO,DFGRD1,NOEL,NPT,KSLAY,
5 KSPT,KSTEP,KINC)
Figure 41. The interface of UMAT subroutine
In Figure 41, we notice that there are many parameters defined for this subroutine.
Among them, the most important parameter is DDSDDE that is elasticity (Jacobian)
matrix of the constitutive model. DDSDDT is defined in an incremental form,
OAA/0As. It has to be updated every iteration step. Figure 42 is my planning algorithm.
This method actually uses the predictorcorrector method. In detail, first we need the
elasticity matrix from the last step to calculate the new stress level, that is, to get a pseudo
new stress. Then, calculate values the modulus and construct the new elasticity matrix.
Recalculate the new stress with this new elasticity matrix. Next, consider the constitutive
60
equation to find the strain rate. Implementing the classical numerical methods such as the
trapezoidal method, we can get the strain tensor in the current step.
The algorithms introduced in Section 4.2 and Section 4.3 cannot be used with
ABAQUS. Another algorithm is designed, as shown in Figure 42. Basically, it uses the
predictorcorrector method to calculate the elasticity matrix.
No results from any of these algorithms can be presented in this dissertation.
Figure 42. Flow chart of the algorithm for Cristescu model combining with ABAQUS
REFERENCES
[1] Aydin I., Briscoe, B.J., Sanliturk, K.Y., "The internal form of compacted ceramic
components: a comparison of a finite element modeling with experiment," Technology and
Medicine, Powder Technology, 1996, Vol.89.
[2] Bickford W.B., Finite Element Method (2nd Edition), Richard D.Irwin, Inc., Boston,
1990 and 1994.
[3] Boresi, A.P., Schmidt, R.J., Sidebottom, O.M., Advanced Mechanics of Materials (5th
Edition), John Wiley & Sons, Inc, New York, 1988.
[4] Cazacu, 0., Jin J., Cristescu, N. D., "A New Constitutive Model for Alumina Powder
Compaction," KONA No. 15, 1997.
[5] Chotourou, H., Gakwaya, A., Guillot, M., and Hrairi, M., "Implementing a cap
material model for the simulation of metal powder compaction," AMDVol.216, Net
Shape Processing of Powder Materials, ASME 1995.
[6] Cristescu, N.D., "A procedure for determining the constitutive equations for materials
exhibiting both timedependent and timeindependent plasticity," Int. J. Solids Structures,
1972, Vol.8.
[7] Cristescu, N.D., "A procedure to determine nonassociated constitutive equations for
geomaterials," Int. J. Plasticity, 1994, Vol.10.
[8] Cristescu, N.D., "Plasticity of Porous and Particulate Materials," Journal of
Engineering Materials and Technology, April 1996, Vol.118.
[9] Cristescu, N.D., Cazacu, 0., and Jin, J., "Constitutive equation for compaction of
ceramic powders," IUTAM Symposium on Mechanics of Granular and Porous Materials,
University of Cambridge, July 1517, 1996.
[10] Cristescu, N.D., Hunsche, U., Time Effects in Rock Mechanics, John Wiley & Sons,
Chichester, 1998.
[ 11] Desai, C.S., Siriwardane, H.J., Constitutive Laws for Engineering Materials i ilh
Emphasis on Geologic Materials, PrenticeHall, Inc., Englewood Cliffs, New Jersey,
1983.
[12] DiMaggio, F.L., Sandler, I.S., "Material model for granular soils," Journal of the
Engineering Mechanics Division, 1971, Vol.97.
[13] Drucker, D.C., and Prager, W., "Soil mechanics and plastic analysis or limit design,"
Quarterly of Applied Mechanics, Vol.10, 1952.
[14] Hibbitt, Karlsson and Sorensen, Inc., ABAQUS Theory Manual Version 5.4, 1994.
[15] Hibbitt, Karlsson and Sorensen, Inc., ABAQUS User's Manual Version 5.4, 1994.
[16] Jin, J., "Elastic/Viscoplastic Models for Geomaterials and Powderlike Materials with
Applications," Ph.D. Dissertation, University of Florida, 1996.
[17] Jin, J., Cazacu, 0., and Cristescu, N.D., "An elastic/viscoplastic consolidation model
for alumina powder," Proc. 5th World Congress of Chemical Engineering, Vol.5, July
1996, San Diego.
[18] Jin, J., Cristescu, N.D., "Experimental study on the consolidation of alumina
powder," Preprint ASME Mechanics and Materials Conference. The Johns Hopkins
University, 1996.
[19] Malvern, L.E., Introduction to the Mechanics of a Continuous Medium, Prentice
Hall, Inc., Englewood Cliffs, New Jersey, 1969.
[20] Owen, D.R.J. and Hinton, E., Finite Elements in Plasticity: Theory and Practice,
Pineridge Press, Swansea, U.K., 1980.
[21] Schofield, A., Wroth, P., Critical State Soil Mechanics, McGrawHill, London,
1968.
[22] Tran, D.V., Lewis, R.W., Gethin, D.T., and Ariffin, A.K., "Numerical Modeling of
Powder Compaction Processes: Displacement Based Finite Element Method," Powder
Metallurgy, 1993, Vol.36, No.4.
[23] Watson, S.A., Adams, M.J., Rough, S.L., Briscoe, B.J., and Papathanasiou, T., "The
modeling of the influence of wall friction on the characteristics of pressed ceramic parts,"
IUTAM Symposium on Mechanics of Granular and Porous Materials, Kluwer Academic
Publishers, Netherlands, 1997.
[24] Zienkiewicz, O.C., The Finite Element Method in Engineering Science (The second
edition, expanded and revised, edition of The Finite Element Method in Structural and
Continuum Mechanics), McGrawHill, London, 1971.
BIOGRAPHICAL SKETCH
Weijian Wang was born on May 10, 1971, in Hebei Province, People's Republic of
China. He received his B.E. degree in the Department of Engineering Mechanics,
Tsinghua University in July 1993. Afterwards he worked as an assistant engineer in China
Aerospace Agency. In May 1997, he was admitted into the Department of Aerospace
Engineering, Mechanics, and Engineering Science, University of Florida. Two years later,
in August 1999, he received his M.S. degree from University of Florida.