NUMERICAL MODELING OF COMPACTION
OF PARTICULATE SYSTEMS
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
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
ACKNOWLEDGMENTS ...................................................... ii
ABSTRACT ...................................................................................... v
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 Drucker-Prager/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
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.
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 oil-well 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 defect-free,
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 Cam-Clay model21 and the Drucker-Prager
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 Drucker-Prager 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 Cam-Clay model and the Drucker-Prager model are their easy-to-
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,6-10,16-18 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.
CLASSICAL MODELS FOR POWDER MATERIALS
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 stress-strain 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 post-yield 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 non-metal
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 (Cam-clay model) was developed by Roscoe and Schofield
et al., 1968 and Parry, 1972, at Cambridge University. This model is used often for
modeling non-cohesive 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
p and q are the mean stress and the equivalent stress respectively. They are defined
P =- (, + CY 2 +3)
S is the stress deviator tensor.
S = T +pI (2.7)
T is the stress tensor.
M and P are user-specified material parameters and a is the evolution parameter.
Figure 2-1 shows how the yield surface of the Cam-clay model looks like.
Critical State Line
" Yield surface
Figure 2-1. Typical yield surface for Cam-clay model
Cam-clay 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 Drucker-Prager/cap model.
2.3 Drucker-Prager/Cap model
A typical yield surface of this model is shown in Figure 2-2. This surface is shown in
the q-p 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).
Figure 2-2. Typical yield surface for Drucker-Prager/Cap model
Similar to the Cam-clay model, the yield surface in the modified Drucker-Prager/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,
(p-pa)2 + -R(d+pa tan3)= 0 (2.9)
where R and ca are material parameters that control the shape of the cap. Pa is an
evolution parameter and defined as
Pa = b R (2.10)
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-
The function involved in the transition segment is
(P-pa)2+ q-(1- )(d+patanp) -u(d+patanp) =0 (2.11)
Plastic flow is non-associated to the shear failure surface and the transition segment,
but is associated to the cap surface. As shown in Figure 2-3, the flow potential surface is
made up of two parts -
shear surface respectively.
The functions of Go
G, = (p -Pa)2 +
Gs = [(p-Pa)tanf
Go and Gs. They are associated with the cap surface and the
and Gs are written as,
cos P3 (2.12)
2 +( q )2
G, (Shear failure)
- \ ..- --- - \. I ap) (l4Ua-a Seci
Figure 2-3. Flow surface for Drucker-Prager/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)
in which ev represents the volumetric strain, Vo and Vi are the initial and final
8 = PI o = PO 1 (2.14)
p1 =Po (2.15)
Equation 2.16 will be used to compute the density variation during the compaction
described in the following examples.
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.  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
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 2-4 and
a. Solid cylinder
b. Cylindrical pipe
c. T-shape d. Cup-shape
Figure 2-4. Different shape powder products
r to I
a. Pyramid-shape container b. Cone-shape container
Figure 2-5. Different shape containers
2.5.2 Computational model
In this chapter, two cases are used as typical samples for computation. See Figure
2-6 and Figure 2-7. 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).
Figure 2-6. Case I- cylindrical sample apparatus for compaction
Figure 2-7. Case II Cup shape sample apparatus for compaction
2.5.3 FEM model and mesh
In both cases 8-node rectangular element is used. Due to the symmetry, only half of
the cross section is considered. Figure2-8, Figure 2-9 and Figure 2-10 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 2-1 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 2-1 Table of Material parameters
-"- 1424 1824 2224
1222 1622 2C22
1422 1822 2222
#1^, 1220 1620 2C20
14:0 1820 2220
89) r (B1@ (B (B@ B
, 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 2-9. FEM mesh and boundary conditions for cylindrical shape sample
@ 0 @ @@@@
0 @ @@@@
@ @ @ @ @@(@
@ @ @ @ @@@@
1 3 5 7 9 11
- -H- 8013
- -H- 8015
Figure 2-10. FEM mesh for cup shape sample
801 803 805
1 18 48 1 1 88 1 1
i | 8019
I~p1' 8 8j.1~
Figure 2-11. Boundary conditions for cup shape sample
The cap hardening curve shown in Figure 2-11 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 2-11. 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 2-12.
Figure 2-13, Figure 2-14 and Figure 2-15 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 cross-section.
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 2-12. Deformation after loading
Figure 2-12. Deformation after loading
- 2. BJE+o1
- 9 OIO1
Figure 2-13. Axial stress distribution after loading
DfBRK1 f lLVE
-- 3. 3; i-01
-- 7. 1z-O
--] TO -01
Figure 2-14. Volumetric strain distribution after loading
-i. AE, QD
91. BE OD
. ,IE UD
2 2DE D0
,2 DE 3 OD
3E tE UD
Figure 2-15. Density distribution after loading
2.6.2 Case 1.2 Cylindrical shape sample, uniform displacement on the top
Figure 2-16, Figure 2-17, Figure 2-18 and Figure 2-19 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
Figure 2-16. Deformation after loading
~. OJt' II
Figure 2-17. Axial stress distribution after loading
s i- .
-5, OT 1- I 1
-i *;|.i -I
-t j11 .11
Figure 2-18. Volumetric strain distribution after loading
+2. 56O +00
+2. 67 +00
+2. 1t OQ
Figure 2-19. Density distribution after loading
The difference between the elastic model and plastic model is apparent, as shown
in Figure 2-20, especially for powder materials. Thus, only the elastic model is not
enough for the study of the powders and geomaterials.
0.20 0.25 0.30 0.35
Figure 2-20. 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
0.05 0.10 0.15
0.20 0.25 0.30 0.35
Figure 2-21. Comparison of the cases with and without friction
0.05 0.10 0.15
Figure 2-22 and Figure 2-23 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.
0.000 0.001 0.002 0.003
0.004 0.005 0.006
Figure 2-22. Under displacement control, the density variation at point 1818
0 10 20 30
40 50 60
Figure 2-23. Under stress control, the density variation at point 1818
Figure 2-24 shows the general trend of the stress-strain curve. It conforms to the
tendency exhibited in the hardening curve of this material. (See Figure 2-11)
0.00 0.05 0.10 0.15 0.20 0.25 0.30
Figure 2-24. Stress-strain 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
Drucker-Prager/Cap model. Figure 2-25, Figure 2-26 and Figure 2-27 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 T-shape.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. 02 E+00
+1. L E+00
Figure 2-25. Axial stress distribution after loading
Figure 2-26. Volumetric distribution after loading
U VA.h2 VALUB EI
t1], U EtDO
4], I Et00
4 1, l 6 E4-0i0
1. 00 E+00
J] .2Et0- -
Figure 2-27. Density distribution after loading
Drucker-Prager/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 2-4 and Figure 2-5 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 Drucker-Prager/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
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
A TIME DEPENDENT MODEL FOR ALUMINA POWDER CRISTESCU MODEL
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 stress-strain
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'6-10'16-18 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 non-linear. 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 stress-strain curve raises starting
from the lower-value 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 elasto-plastic material, four
different assumptions are required:
1. A yield function, which generalizes the concept of the yield stress to 2-D and 3-D
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
s is the total strain rate.
G and K are elastic moduli.
k is the viscosity parameter.
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
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 3-1 and Figure 3-2 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
0- 7000 -
0 20 40 60 80 100
Figure 3-1. Bulk modulus variation
2 500 -
0 100 200 300 400 500 600
Figure 3-2 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
H(c) = Hh (G) + Hd (Y, T) (3.3)
In the hydrostatic stage, Hd(o, c) = 0.
Hh ()= ah ()2 +bh () (3.4)
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.
; =k I- (NI+N2 -) (3.7)
kNi and kN2 can be calculated by the following formulas,
kN, = q(C) + sgn X(o, Tc)Y(c, -c)
cy a-- (3.8)
kN2 =m(_)14 m2(- ) +n2(J)2exp n2(- +r
Pa Pa Pa Pa
p(c) = ak C2 + bkC
X(c,'c) = 0.555 -c (3.9)
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
In a general loading test, one can get the stress-strain 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)
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,
W = (1 W)(kN I + kN2 ) (3.11)
l l = 3C (3.12)
= j (cj C61) = cj 3C2 (3.13)
dW = [3kNkNN +C( -3Z2 3)]dt (3.14)
Integrate it in the time interval [tn, tn+i],
W'n+ ft [3kN, N2 ( 3C2)]dt (3.15)
H ln(1- ) Wn+1 = tn [3kNC + (kN2 3C2)]dt (3.16)
Let FI = 3kNY + kN2 ( o 3C2)
exp(- 1n FIdt)
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.
- (FI + FIn +)
Table 2-1 lists the parameter values used in the computation.
Table 2-1 Table of material parameters
Confining pressure 294 kPa
Axial stress increment 20 kPa
Time increment 0.01 second
bh 1. Ox 10-6
bk 1.641 x 10-6
r 8.0x 10-4
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
392kPa 700 (kPa)
-0.01 -0.008 -0.006 -0.004 -0.002 0 0.002 0.004
Figure 3-3. Stress-strain curves for different confining pressures in a general test
Figure 3-4 and Figure 3-5 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.
200 400 600
Figure 3-4 Irreversible stress work for different load increments, same time increments
Figure 3-5.Irreversible stress work for different time increments, same load increments
Viscosity is another important factor. Figure3-6 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
Figure 3-6.Irreversible stress work under different values of k
3.4 Transient Creep Computation
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 3-7.
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
W(t) = J; dT
= t (76 1+ )(1' 6 + E, )d' (3.20)
to j d
Let o (3.21)
Wd =J t dT
W = Wv+ Wd (3.22)
Calculate the stress at the
Calculate a, a', and c for the
Calculate material parameters
for the current step
E(c), K(c), G(c)
Calculate H and W for the
(s =( Nextl)( ad tep
SCalculate totalversible tra instrain
4Yes current step
w- = ft a:FrNext load stepd
Figure 3-7 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 3-8 shown.
400- 03 = 294 kPa
350- Acy = 20 kPa
0 10 20 30 40 50 60 70 80 90
Figure 3-8 Octahedral shear stress in triaxial test
Figure 3-9 is the typical stress-strain 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 3-10, Figure 3-11 and Figure 3-12 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.
However, the magnitude of volumetric part is much less
-0.05 0 0.05
0.1 0.15 0.2 0.25
Figure 3-9 Stress-strain curves
Wv (kPa) -4 -
100 200 300
Figure 3-10 Volumetric irreversible stress work as a function of c
0 ---------- ------------
0 100 200 300 400 500
Figure 3-11 Deviatoric irreversible stress work as a function of c
0 100 200 300 400 500
Figure 3-12 Total irreversible stress work as a function of c
In Figure 3-13, 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
Figure 3-14 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
-0.01 -0.008 -0.006 -0.004 -0.002 0 0.002 0.004 0.006 0.008
Figure 3-13 Influence of confining pressure on the volumetric strain
In Figure 3-15, we can clearly observe that the volume firstly decreases, which
means the compressibility. Afterwards it increases, which means dilatancy.
100 200 300 400
Figure 3-14 Variation of the irreversible stress work as a function of C for different
0 10 20 30 40 50 60 70 80 90
Figure 3-15 Volumetric strain during the loading procedure
Figure 3-16 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 3-17.
20 40 60 80
Figure 3-16 Volumetric strain for different loading increment and same confining pressure
The curves in Figure 3-17 exhibit the same tendency as that shown in Figure 3-16.
When the load increment is larger, the number of the step is fewer.
2- CD e f
1 B d
A c E
-33 = 294 kPa
0 5 10 15 20 25
Figure 3-17 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 3-18, if k is smaller, the material is more flowable and longer time is
needed to get to the stabilization.
Figure 3-19 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.
From = 0, we can derive out
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 2-D plane through the points of the
maxima on the yield surfaces H = constant.
0 10 20 30 40 50 60 70 80
Figure 3-18. Volumetric strain curve for different value of k
The points shown in Figure 3-20 correspond to those in Figure 3-17. 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= ( +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
0 500 1000 1500 2000
Figure 3-19 Compressibility/Dilatancy boundary, failure surface 8H/8y = 0, and several
Failure C/D boundary
& 400 1
200- e0C Arctan(1.4142
0 500 1000 1500 2000
Figure 3-20 Deviatoric test procedure in T-c- plane
3.5 Analytical Solution for the Creep Test with Fixed Time Interval
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, (to-t)
l(kN + kN2 ). L
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 3-21 exhibits the same story as Figure 3-17 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 3-22 shows the results for different step time intervals. If the time interval
is relatively small, we can get the transient creep only.
X 10-3 ..
At = 10 sec
0 10 20 30 40 50 60 70 80
Figure 3-21 Volumtric strain change for different k values
0 10 20
40 50 60 70 80 90
Figure 3-22 Volumetric strain change for different step time intervals
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.
REALIZATION OF CRISTESCU MODEL WITH FEM
Cristescu model is a time-dependent semi-linear model. Usually, the variable
stiffness method is the most acceptable method to solve the non-linear 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.
qn = D-'Ac + At(1-6)sn+sn (4.5)
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
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)
1+1 = + -6 ,
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.
q. = As At(1 0)n + D 'A (4.13)
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
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)
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.
Figure 4-1. The interface of UMAT subroutine
In Figure 4-1, 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 4-2 is my planning algorithm.
This method actually uses the predictor-corrector 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
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 4-2. Basically, it uses the
predictor-corrector method to calculate the elasticity matrix.
No results from any of these algorithms can be presented in this dissertation.
Figure 4-2. Flow chart of the algorithm for Cristescu model combining with ABAQUS
 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.
 Bickford W.B., Finite Element Method (2nd Edition), Richard D.Irwin, Inc., Boston,
1990 and 1994.
 Boresi, A.P., Schmidt, R.J., Sidebottom, O.M., Advanced Mechanics of Materials (5th
Edition), John Wiley & Sons, Inc, New York, 1988.
 Cazacu, 0., Jin J., Cristescu, N. D., "A New Constitutive Model for Alumina Powder
Compaction," KONA No. 15, 1997.
 Chotourou, H., Gakwaya, A., Guillot, M., and Hrairi, M., "Implementing a cap
material model for the simulation of metal powder compaction," AMD-Vol.216, Net
Shape Processing of Powder Materials, ASME 1995.
 Cristescu, N.D., "A procedure for determining the constitutive equations for materials
exhibiting both time-dependent and time-independent plasticity," Int. J. Solids Structures,
 Cristescu, N.D., "A procedure to determine nonassociated constitutive equations for
geomaterials," Int. J. Plasticity, 1994, Vol.10.
 Cristescu, N.D., "Plasticity of Porous and Particulate Materials," Journal of
Engineering Materials and Technology, April 1996, Vol.118.
 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 15-17, 1996.
 Cristescu, N.D., Hunsche, U., Time Effects in Rock Mechanics, John Wiley & Sons,
[ 11] Desai, C.S., Siriwardane, H.J., Constitutive Laws for Engineering Materials i ilh
Emphasis on Geologic Materials, Prentice-Hall, Inc., Englewood Cliffs, New Jersey,
 DiMaggio, F.L., Sandler, I.S., "Material model for granular soils," Journal of the
Engineering Mechanics Division, 1971, Vol.97.
 Drucker, D.C., and Prager, W., "Soil mechanics and plastic analysis or limit design,"
Quarterly of Applied Mechanics, Vol.10, 1952.
 Hibbitt, Karlsson and Sorensen, Inc., ABAQUS Theory Manual Version 5.4, 1994.
 Hibbitt, Karlsson and Sorensen, Inc., ABAQUS User's Manual Version 5.4, 1994.
 Jin, J., "Elastic/Viscoplastic Models for Geomaterials and Powder-like Materials with
Applications," Ph.D. Dissertation, University of Florida, 1996.
 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.
 Jin, J., Cristescu, N.D., "Experimental study on the consolidation of alumina
powder," Preprint ASME Mechanics and Materials Conference. The Johns Hopkins
 Malvern, L.E., Introduction to the Mechanics of a Continuous Medium, Prentice-
Hall, Inc., Englewood Cliffs, New Jersey, 1969.
 Owen, D.R.J. and Hinton, E., Finite Elements in Plasticity: Theory and Practice,
Pineridge Press, Swansea, U.K., 1980.
 Schofield, A., Wroth, P., Critical State Soil Mechanics, McGraw-Hill, London,
 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.
 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.
 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), McGraw-Hill, London, 1971.
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.