The Mobilization Force in a Confined Granular Column
Dominik Schillinger Pre-diploma (Civil Engineering), University of Stuttgart, Germany, 2004
A Thesis Submitted in Partial Fulfillment of the Requirements for the Degree of Master of Science at the University of Connecticut 2006
ACKNOWLEDGEMENTS I would like to express my deepest and sincere gratitude to my major advisor Dr. Ramesh B. Malla for his guidance and encouragement during the last year. Working with Dr. Malla has been a very fruitful and enriching experience, from which I have benefited both academically and personally. His hard work in research, his dedication to teaching and his constant concern for his students contributed largely to a continuous and efficient advancement in my graduate work. I would also like to thank Dr. Malla for organizing financial support throughout the year. It was partially provided by Hamilton Sundstrand, Windsor Locks, CT, the NASA EPSCoR program (NASA Experimental Program to Stimulate Competitive Research), and the UConn School of Engineering, to whom all I would like to extend my sincere thanks. I would also like to thank my associate advisors Dr. Michael L. Accorsi and Dr. Jeong-Ho Kim for reviewing the manuscript of this thesis and providing valuable comments. Special appreciation goes to Johnathan Drasdis and Robert Wakefield, who helped in conducting experiments and gave feed-back and lots of new ideas, and to the Department of Civil & Environmental Engineering for its state-of-the-art research facilities and logistic support. Moreover, I wish to express my deepest gratitude to Dr. Renate Seitz from the Department of Higher Education of the State of Connecticut. My studies at the University of Connecticut could not have taken place without her commitment. Finally, I would like to thank my family and my friends for their love. Without their encouragement, academic success would have never been possible.
iii
TABLE OF CONTENTS
Title Page
i
Approval Page
ii
Acknowledgements
iii
Table of Contents
iv
List of Figures
ix
1. Introduction 1.1 Research Motivation
1
1.2 Objectives of the Work:
3
The “Zero-axial-stress” Problem and the Mobilization Problem 5
1.3 Scope and Structure of the Work
2. Multi-Disciplinary Background of the Physics of Granular Materials 2.1 Coulomb’s Law of Friction 2.1.1
Physical Background and Continuum Assumption
7
2.1.2
Coulomb’s Friction Coefficient
8
2.1.3
Contact Formulation between Grains and Tube
8
2.1.4
Mohr-Coulomb Criterion and Internal Friction Angle
9
2.2 Janssen’s Bin Analysis 2.2.1
Assumptions and Model Derivation
12
2.2.2
Typical Model Behaviour
14
iv
2.2.3
K-value Determination
15
2.3 Force Chains and Arching 2.3.1
The Cannon Ball Example
19
2.3.2
Force Chains
20
2.3.3
Arching
22
2.4 Geotechnical Classification Parameters 2.4.1
Grain-size Distribution
24
2.4.2
Weight-volume Relations
26
2.4.3
Unit Weight, Porosity and Void Ratio for Alumina Granulate
27
3. An Analytical Elastic Solution Approach by the Fourier-Bessel Potential 3.1 Fundamentals of Theory of Elasticity 3.1.1
Basic Assumptions of Linear Elasticity
29
3.1.2
Equilibrium and Kinematic Equations
30
3.1.3
Constitutive Relations, Compatibility and Boundary Conditions
32
3.1.4
Potential Functions
34
3.1.5
Axisymmetry and Cylindrical Coordinates
35
3.2 Love’s Potential Function for Axisymmetry 3.2.1
Derivation of the Stress Relations
38
3.2.2
Derivation of the Displacement Relations
43
3.3 Fourier-Bessel Series Approach 3.3.1
The Fourier-Bessel Potential for Z
45
3.3.2
Stress and Displacement Expressions
47
3.3.3
Functional Orthogonality
48
3.3.4
Fourier Series and Finite Fourier Transforms
49
v
3.4 Boundary Conditions 3.4.1
Basic Requirements
51
3.4.2
The Janssen Shear Stress Assumption
52
3.4.3
Application of Transform Methods
53
3.4.4
Superposition of Polynomial Solutions
56
3.4.5
Eigenfunction Analysis
57
3.5 Error Analysis and Improvement 3.5.1
Periodicity vs. Monotonic Decrease
60
3.5.2
Sampling and Truncation Error
63
3.5.3
The Janssen Error and the Improvement Algorithm
66
3.6 The “Zero-axial-stress” Problem in Continuum Modeling 3.6.1
The Simple Janssen Model
70
3.6.2
Continuum Theory vs. Discrete Particle Behaviour
71
3.6.3
The Remedy of a Small Constant Bottom Stress
72
4. Experimental Parameter Identification and Mobilization Testing 4.1 Coulomb’s Friction Coefficient 4.1.1
Basic Testing Concept
75
4.1.2
Critical Discussion
76
4.1.3
Results
77
4.2 Elastic Material Properties 4.2.1
Loading-Unloading Behaviour of Granular Materials
80
4.2.2
The Uniaxial Strain Test
81
4.2.3
The Triaxial Test Set-up
82
4.2.4
Triaxial Testing for Young’s Modulus
83
vi
4.3 Material Properties for the Cap Model 4.3.1
The Mohr-Coulomb Failure Envelope
86
4.3.2
Porosity Change and Plastic Volumetric Strains
87
4.3.3
The Triaxial Flowmeter Test
88
4.4 Mobilization Test 4.4.1
Experimental Set-up in “Reverse”
91
4.4.2
Results
93
5. The Elasto-plastic Drucker-Prager Cap Model 5.1 Elasto-plastic Modeling in Geomechanics 5.1.1
The Incremental Flow Theory of Plasticity
95
5.1.2
Yield Criterion, Flow Rule and Plastic Potential
96
5.1.3
Generalized Stress-Strain Relations
99
5.2 The Drucker-Prager Model 5.2.1
The Drucker-Prager Yield Criterion
101
5.2.2
Flow Theory Formulation for Drucker-Prager
103
5.2.3
The Radial Return Algorithm for FE Implementation
104
5.3 The End Cap Extension 5.3.1
The Concept of the End Cap
106
5.3.2
Cap Surface and Plastic Hardening Computations
107
5.3.3
A Modular Code Structure for FE Implementation
110
5.4 The Drucker-Prager Cap Model in Simple Tests 5.4.1
General Incremental Finite Element Procedure
113
5.4.2
A Simple One Element FE Model
114
vii
5.4.3
Isotropic Consolidation
116
5.4.4
Triaxial Shear Testing
119
6. Results and Discussion 6.1 Elastic Stresses and Displacements in a Confined Granular Column 6.1.1
Exact Elasticity Solution from the Fourier-Bessel Potential
121
6.1.2
Elasticity Solution from Finite Element Analysis
124
6.1.3
Result Plots
125
6.2 Elasticity vs. Elasto-plastic and Experimental Material Behaviour 6.2.1
The Uniaxial Strain Test
132
6.2.2
Results Plots
133
6.2.3
Discussion and Model Evaluation
138
6.3 The Mobilization Force in a Confined Granular Column 6.3.1
The Mobilization Force and Continuum Modeling
142
6.3.2
A Janssen and a FEM Approach
143
6.3.3
Discussion and Evaluation
148
7. Conclusions and Recommendations 7.1 Summary of the Work
151
7.2 Conclusions
152
7.3 Recommendations for Future Research
154
Appendix A: MATLAB code for the Fourier-Bessel Method
158
Appendix B: Fortran User subroutine UMAT for the Drucker-Prager cap model
161
Complete References & Bibliography
169
viii
LIST OF FIGURES
Fig. 1.1
A finite length of granular material stuck in a tube
1
Fig. 1.2
Finite length granular column and pushing top piston
2
Fig. 1.3
Schematic of typical water processing bed
2
Fig. 1.4
“Zero-axial-tress” problem: Formulation as an Elastic Boundary Value
4
Problem Fig. 2.1
Free body diagram: Single grain – Tube wall
7
Fig. 2.2
Microscopic asperities
8
Fig. 2.3
Rigid-plastic behaviour of an ideal Coulomb material
10
Fig. 2.4
The Mohr-Coulomb failure criterion
11
Fig. 2.5
Janssen’s cylinder
12
Fig. 2.6
Janssen’s differential slice
13
Fig. 2.7
Mohr’s circle for incipient failure
16
Fig. 2.8
Uniform compression
17
Fig. 2.9
Example of Janssen formulae
18
Fig. 2.10 A possible actual force scenario and corresponding contact points
19
Fig. 2.11 Force distribution for idealized balls and stresses for elastic model
19
Fig. 2.12 Force chain
20
Fig. 2.13 Stress pattern in a 2D granular material under compression after Dantu
21
(1957) Fig. 2.14 Photo-elastic experiment done by Travers et al. (1987)
21
Fig. 2.15 Arch structures are stable and endurable (Aqueduct of Maintenon)
22
Fig. 2.16 Stress pattern from a Contact Dynamics simulation by Radjai (1997)
22
ix
Fig. 2.17 Particle size classification according to the American Association of State
24
Highway and Transportation Officials (AASHTO) Fig. 2.18 Particle size distribution of alumina granulate (dots) and Ottawa sand
25
(triangles) as shown in (Anandakumar 2005) Fig. 2.19 Parameters for weight-volume relations and separation of the three phase
26
system Fig. 3.1
Stresses on a differential slice in plane strain
31
Fig. 3.2
Displacements and strains
31
Fig. 3.3
Stress and displacement boundary conditions
33
Fig. 3.4
Polar coordinates r, θ
35
Fig. 3.5
Stresses on a differential element in cylindrical coordinates
36
Fig. 3.6
“Zero-axial-tress”: Problem formulation
51
Fig. 3.7
Uniform compression without shear or friction at the wall
57
Fig. 3.8
First three Eigenfunctions for the Fourier- Bessel potential
59
Fig. 3.9
Exponential Janssen shear and its Fourier series representation
62
Fig. 3.10 Sampling error on z ∈ [12; L=20] for n=20
62
Fig. 3.11 Truncation error on z ∈ [8; 10] for n=20
65
Fig. 3.12 Coulomb’s coefficient μ close for n=50
65
Fig. 3.13 Simple algorithm for improving μ close
68
Fig. 3.14 K-functions after 1st and 2nd iterations
69
Fig. 3.15
μ close after 1st and 2nd iteration
69
Fig. 3.16 Janssen’s differential slice
70
Fig. 3.17 Arching and force chains in a cylinder
73
Fig. 3.18 Bottom stress for continuum approach
73
Fig. 4.1
Experimental set-up for friction coefficient test
75
Fig. 4.2
Sensor-measured axial stress induced by piston loadings of 0.25 – 0.5 –
78
1.0 – 1.5 N/mm² x
Fig. 4.3
Frictional shearing stress at which the cylinder moved downwards
79
Fig. 4.4
Friction coefficient
79
Fig. 4.5
Experimental set-up for the uniaxial strain test
81
Fig. 4.6
Principal stresses on the surface of the triaxial specimen
82
Fig. 4.6
Experimental set-up for the triaxial test
83
Fig. 4.8
Separation of the Triaxial into a Hydrostatic and a Uniaxial stress state
84
Fig. 4.9
Loading-unloading lines in the uniaxial strain test for max stress 0.1 N/mm²
85
determined from the polynomial regression functions
Fig. 4.10 Loading-unloading lines in a triaxial test for max piston stress 0.1 N/mm²
85
Fig. 4.11 Direct Reading Flowmeter
88
Fig. 4.12 Mohr-Coulomb failure criterion for a the cohesionless alumina granulate
90
with an internal friction angle of 37° Fig. 4.13 Hardening parameter
90
Fig. 4.14 Experimental set-up for the mobilization test
92
Fig. 4.15 Mobilizing stresses for three different bed lengths
94
Fig. 4.16 Load-displacement for bottom (Loading machine platen) and top (Dial
94
gage) of the granular column with 0.0045 N/mm² Top Weight Fig. 5.1
Common yield criteria in soil plasticity
98
Fig. 5.2
Matching of yield criteria in the deviatoric plane
102
Fig. 5.3
Drucker-Prager type of strain hardening cap model
108
Fig. 5.4
Yield surfaces for the Drucker-Prager model in the meridional plane
108
Fig. 5.5
Axisymmetric element and its position in the cylinder
115
Fig. 5.6
Single element model with boundary conditions as used in the following
115
Fig. 5.7
Stress path / cap movement in the meridional plane for hydrostatic stress
118
Fig. 5.8
Evolution of volumetric strain under loading–unloading isotropic
118
Consolidation Fig. 5.9
Stress path/cap movement in the meridional plane for triaxial shear testing
Fig. 5.10 Evolution of volumetric strain in triaxial shear testing up to failure xi
120 120
Fig. 6.1
Boundary functions from Fourier-Bessel analysis
122
Fig. 6.2
Deformed axisymmetric elastic FE mesh (middle part not displayed)
123
Fig. 6.3
Axisymmetric plane over which the stress and displacement results are
125
plotted as the third dimension Fig. 6.4
Radial stress σ z
126
Fig. 6.5
Circumferential stress σ z
127
Fig. 6.6
Axial stress σ z
128
Fig. 6.7
Shear stress τ r z
129
Fig. 6.8
Radial displacement u r
130
Fig. 6.9
Axial displacement u z
131
Fig. 6.10 Axisymmetric finite element model for the uniaxial strain test
133
Fig. 6.11 Radial stress σ r
134
Fig. 6.12 Axial stress σ z
135
Fig. 6.13 Shear stress τ r z
136
Fig. 6.14 Radial and axial displacements u r and u z
137
Fig. 6.15 Shear behaviour of elements close to the wall (stress path 1);
139
Fig. 6.16 Finite length granular column and pushing top piston
143
Fig. 6.17 FEM for determination of the mobilization force
145
Fig. 6.18 Mobilization stress at the loading end for 3 different bottom support
146
stresses (0.002, 0.006, 0.010 N/mm²) determined by FEM and Janssen Fig. 6.19 Frictional shear stresses along the wall for 3 different bottom support
146
stresses (0.002, 0.006, 0.010 N/mm²) determined by FEM and Janssen Fig. 6.20 Janssen, FEM and experimental results showing the linear dependence of
147
the mobilization stress at top from the bottom support stress Fig. 6.21 Janssen and experimental results showing the exponential dependence of the mobilization stress at top from the column length
xii
147
Chapter 1: Introduction
The case of a cylindrical container filled with granular material is a basic situation where the presence of the wall boundaries induces horizontal pressure and shearing forces caused by friction between wall and grains. It can be found in various engineering applications, e.g. in caissons construction and pile driving, in the design of silos and storage bins or in space life support systems, where granular material contained in a steel tube is compacted by a spring.
1.1 Research Motivation
The phenomenon commonly referred to as arching can make vertical pressure at a certain depth of the granular cylinder zero, which can lead to a column of granular material with free top and bottom surfaces. It is held in vertical equilibrium only by the friction forces at the tube wall (see Fig. 1.1). When a piston with gradually increased load is applied on top (see Fig. 1.2), the whole granular column starts to move downwards at a Fig. 1.1 A finite length of granular material stuck in a tube
critical piston force, as soon as the top load
1
cannot be equilibrated by the frictional wall shear. This critical force, which is referred to as the mobilization force, and the corresponding stress and deformation behaviour can be used beneficially in the above mentioned practical fields. We will give an example from the design of water-processing beds in
Fig. 1.2 Finite length granular column and pushing top piston
space life support systems (Malla, Gopal 2002a), (Malla, Gopal 2002b),
(Malla, Anandakumar 2004), Malla, Anandakumar 2006). Fig. 1.3 shows a typical water processing bed intended for purifying wastewater on board of space stations. The granular packed bed is enclosed in a steel tube and compressed by a spring. In order to work efficiently, the granular bed is required to be in a compressed state throughout its axial length.
Fig. 1.3 Schematic of Typical Water Processing Bed
2
The spring on top of the bed induces axial compression that is, however, reduced by wall friction with increasing depth. A situation as shown in Fig 1.1 may occur at a certain depth leaving the rest of the bed uncompressed, which leads to inefficient operation conditions or a malfunction of the bed in the worst case. If the friction mobilization force for the total length of the bed is known, a compression of the whole bed can be insured by applying this force at the top by the spring.
1.2 Objectives of the Work: The “Zero-axial-stress” Problem and the Mobilization Problem
In this work, we mainly deal with modelling the stress and deformation behaviour of a granular column, which is pushed from the top surface (see Fig. 1.2). We concentrate on two aspects of this situation: First, we want to know if it is possible to make axial stresses zero at a finite depth in a confined elastic medium. If yes, we can model the arching effect in the granular column directly by a linear elastic model and the well-known stress-strain relations from continuum mechanics. Therefore, the subsequent set of stress and displacement boundary conditions is applied to an isotropic elastic cylindrical body (see Fig. 1.4):
At z = 0 :
σ z (r ) = σ 0
(Loading stress from the piston)
At z = L :
σ z (r ) = 0
(Zero axial stress at a finite depth L)
At r = R :
ur ( z) = 0
(Confinement by the tube wall)
3
μ=
τ rz ( z ) σ r ( z)
(Coulomb’s friction)
In the following, we refer to this question as the “zero-axial-stress” problem. Second, we want to inquire how accurate linear elasticity can predict stresses and displacements in a confined granular column in comparison to a refined elasto-plastic Drucker-Prager cap model and to experimental results. Third, we want to know at what critical piston force the granular column will start to slide down in a rigid body motion and try to develop theoretical models to predict this mobilization force. Additionally, we are interested in the difference of linear elastic and elasto-plastic results for the mobi-
Fig. 1.4 “Zero-axial-tress” problem: Formulation as an Elastic Boundary Value Problem
lization force. For the description of stresses and displacements, we follow a continuum approach. It is assumed that length scales are sufficiently large, so that the discrete nature of matter can be neglected, and the material and its properties can be assumed to be continuously distributed. For modelling the mechanical stress-displacement behaviour of the confined granular material, we can hence use the well established methods of continuum mechanics. Note that in the scope of this work, the focus lies on constitutive modelling of the granular body rather than on a detailed tribological description of the wallgrain contact interface, which is therefore approximated by simple Coulomb’s friction.
4
For experimental verification, we use alumina granulate that is confined in a plastic tube made out of polymethylmethacrylate (PMMA). Note that in order not to destroy the particles during the loading tests, stresses are not to exceed 0.2 N / mm 2 anytime.
1.3 Structure and Scope of the Work
This thesis is organized in seven chapters. Chapter 1 gives explains the research motivation (section 1.1), defines the main problems that are examined (section 1.2), and give an outline of the structure of this work (section 1.3). Chapter 2 is primarily intended to provide preliminary information on Coulomb’s friction, arching and common material parameters for granular materials in general and alumina granulate in particular. Furthermore, we try to give a short introduction to the field of granular materials from the perspective of different disciplines that have been working on the topic, namely agricultural engineers (section 2.2), physicists (section 2.3) and geotechnical engineers (section 2.4). Chapter 3 is devoted to the “zero-axial-stress” problem. We find an analytical, exact solution by means of Fourier-Bessel analysis (sections 3.2 - 3.4) and prove that the “zero-axial-stress” problem with the set of boundary conditions stated above cannot be solved by elasticity (section 3.5). Therefore, we propose the remedy of a constant bottom stress for the representation of the arching effect in continuum models (section 3.6). Chapter 4 deals with all experiments carried out in the scope of this work. On the one hand, it discusses the determination of material parameters for the example case of alumina granulate confined in the plastic PMMA tube, i. e. the friction coefficient between
5
the plastic tube surface and grains (section 4.1) and material parameters required for the elastic and elasto-plastic constitutive relations (sections 4.2 and 4.3). On the other hand, we suggest a testing method to determine the friction mobilization force in a confined column of alumina granulate of different depths (section 4.4). Chapter 5 provides a brief introduction to elasto-plastic constitutive modelling (section 5.1). We focus on the Drucker-Prager cap model, its implementation in the finite element method (sections 5.2 and 5.3) and its modelling capabilities in terms of volumetric compaction and shear failure (section 5.4). The first part of the final Chapter 6 presents the exact, analytical stress and displacement fields for a linear elastic granular body confined in a rigid tube, that have been found by the Fourier-Bessel approach of chapter 3. The results are checked by a finite element analysis. The second part of chapter 6 deals with the evaluation of linear elasticity for modelling stresses and displacements in granular materials in comparison with elasto-plastic (Drucker-Prager cap model) and experimental results. In the third part of chapter 6, we attempt to predict the mobilization force in the alumina granulate confined in the PMMA tube using a simple elastic Janssen model and a computational elastoplastic FE model. Chapter 7 provides a short summary of all analytical, computational and experimental methods and techniques that are used in this work (section 7.1) and presents conclusions that we can draw from chapter 6 (section 7.2). We conclude by recommending a set of future studies (section 7.3).
6
Chapter 2: Multi-Disciplinary Background of the Physics of Granular Materials
2.1 Coulomb’s Law of Friction
The frictional transfer mechanism that brings shear stresses from the grains to the tube wall constitutes a very important aspect of both the “zero-axial-stress” problem and the mobilization problem. In this work, we will confine ourselves to the elementary concept of Coulomb’s friction.
2.1.1 Physical background and continuum assumption. Consider a single grain that touches the wall as shown in Fig. 2.1. Friction between these two bodies bases on the fact that their physical boundaries are rough from a microscopic point of view. Actual contacts are not continuous, but occur between asperities of the two contacting surfaces (see Fig. 2.2). It is assumed that length scales of the total contact area are sufficiently large. Contact forces that are physically present only at the contact points of the asperities can thus be averaged over the total contact area and therefore assumed to be continuously distributed. This averaged contact force can be split into a tangential part T and a normal part N .
7
Fig. 2.1 Free body diagram: Single grain – Tube wall
The tangential force T tends to produce a relative motion between the two contact surfaces. The asperities, however, offer resistance against that relative motion, because they are pushed against each other by the normal force N and get stuck - analogous to a locking mechanism.
Fig. 2.2 Microscopic asperities
2.1.2 Coulomb’s Friction Coefficient. The tangential force T has to reach a critical value to release the locking and make the asperities slide over each other. Coulomb (1776) postulated that this critical value
T crit can be found proportional to the normal part N in the form T crit = μ ⋅ N
(2.1)
where μ is called Coulomb’s coefficient of friction. It can be differentiated in a static coefficient μ s that is valid before any relative motion between the bodies occurs, i. e. the material is just about to fail, and a smaller dynamic coefficient μ d that is valid after the onset of sliding.
2.1.3 Contact Formulation between Grains and Tube. Frictional contact problems in analytical mechanics are very involved and solutions only possible for special cases (Johnson 1985). Therefore, our first concern in the description of the contact between granular cylinder and tube wall must be simplicity that keeps the cylinder problem accessible for analytical solution methods.
8
We therefore adopt the elementary concept of Coulomb’s friction over the whole contact surface and postulate
At r = R :
μ=
τ rz = constant σr
(2.2)
The friction induced shear τ r z always acts in the opposite direction to relative sliding. The friction coefficient μ = 0.36 for activated alumina on PMMA plastic is determined experimentally in section 4.1. Nevertheless, Eqn. (2.2) is a simplification and is inconsistent in several aspects. Experiments indicate that for activated alumina Eqn. (2.2) is not constant, but exponential (see section 4.1). Moreover, it only applies in the case of incipient failure when the two bodies are about to slide. This is satisfied when the whole cylinder is about to be mobilized. The static coefficient of friction is used all the time, although the granular cylinder shows considerable axial deformation before sliding, which implies that there must be small relative motion over the surface before mobilization.
2.1.4 Mohr-Coulomb Criterion and Internal Friction Angle. The Coulomb’s condition Eqn. (2.1) can also be used as a simple approach to describe the stability of a granular body itself. We define the material as a rigid-plastic, ideal Coulomb material where the Coulomb’s yield criterion with
τ ≤ μ ⋅σ + c
(2.3)
must be satisfied at any point of the body.
9
The parameter c is called the cohesion and can be intuitively interpreted as the material’s capacity of taking some small amount of tension (see Fig. 2.4). If the two sides of Eqn. (2.3) are equal, the material is about to slide (incipient failure); if Eqn. (2.3) cannot be satisfied, the material fails. This behaviour can be illustrated in a simple example: Consider a granular body subjected to a force as shown in Fig. 2.3. At some point, if P is gradually increased, Eqn. (2.3) cannot be satisfied anymore at all points. Eventually, the material will divide itself into two blocks which slide past Fig. 2.3 Rigid-plastic behaviour of an ideal Coulomb material
each other in a so-called slip-plane along those points.
We define the so-called angle of internal friction φ as
μ = tan(φ )
(2.4)
and find the Mohr-Coulomb criterion from Eqn. (2.3)
τ ≤ tan(φ ) ⋅ σ + c
(2.5)
Eqn. (2.5) can be interpreted graphically in a Mohr’s circle. For a material point not to be sliding, the Mohr’s circle must either stay below the internal yield locus IYL, which is determined by φ , or touch it (incipient failure). Thus Eqn. (2.5) imposes a limit on the magnitude of the shear stress that can occur in a granular material before failure.
10
Fig. 2.4 The Mohr-Coulomb failure criterion
11
2.2 Janssen’s Bin Analysis
Janssen (1895) proposed his well-known analysis to predict stresses developed by granular materials in silos. This simple solution is still widely used, i. e. in codes of practice for bunker design or as basis or checking for more refined models in research.
2.2.1 Assumptions and Model Derivation. Janssen’s starting point was the observation that grains confined in a silo had a marked proclivity to redirect vertically applied stresses towards the sides, where they are resisted by wall friction. He decided to treat the grains as a continuous medium with homogeneous density γ and focused on a differential slice taken out at some depth z from the cylinder in Fig. 2.5. Furthermore, Janssen based his analysis on two simplifying assumptions that contradict basic princi-
Fig. 2.5 Janssen’s cylinder
ples of continuum mechanics:
(1) The axial stress σ z is uniform across any horizontal section of the material.
(2) The axial and radial stresses σ z and σ r are principal stresses. Moreover, at any point of the granular body their ratio K is constant in the form
K=
σz = constant σr
(2.6)
12
With assumptions (1) and (2) in mind, we can perform a force balance on the differential slice as shown in Fig. 2.6
π R 2 σ z + π R 2 γ δ z = π R 2 (σ z + δσ z ) + π R δ z τ w
(2.7)
We can simplify Eqn. (2.7) by using Eqn. (2.6) and Coulomb’s friction
τ w = μ ⋅σ r
(2.8)
with μ being the Coulomb’s coefficient of friction at the wall, and find
∂σ z 2 μ K + σz =γ R ∂z
(2.9)
Fig. 2.6 Janssen’s differential slice
For this simple first order differential equation, we can write down the solution by inspection
σz =
γR ⎛ 2μ K + C exp⎜ − R 2μ K ⎝
⎞ z⎟ ⎠
(2.10)
13
where C is an arbitrary constant to be determined by the boundary condition for σ z at
z = 0 . We assume that the top surface is subjected to a uniform surcharge σ 0 .
σz =σ0
At z = 0 :
(2.11)
Substitution of Eqn. (2.11) into Eqn. (2.10) and application of Eqns. (2.6) and (2.8) bring about the following classical Janssen formulae
σz =
γR 2μ K
⎞ z⎟ ⎠
(2.12)
σr =
γR⎡ ⎛ 2μ K ⎞ ⎛ 2 μ K ⎞⎤ 1 − exp⎜ − z ⎟⎥ + K σ 0 exp⎜ − z⎟ ⎢ 2μ ⎣ R ⎠⎦ R ⎠ ⎝ ⎝
(2.13)
τw =
⎡ ⎛ 2μ K ⎢1 − exp⎜ − R ⎝ ⎣
⎞⎤ ⎛ 2μ K z ⎟⎥ + σ 0 exp⎜ − R ⎠⎦ ⎝
γR⎡
⎛ 2μ K ⎞ ⎛ 2μ K ⎞⎤ 1 − exp⎜ − z ⎟⎥ + μ K σ 0 exp⎜ − z⎟ ⎢ 2 ⎣ R ⎠⎦ R ⎠ ⎝ ⎝
(2.14)
2.2.2 Typical Model Behaviour. Eqn. (2.12) – (2.14) are monotonic decreasing and exponential. All three of them tend asymptotically to constant values at great depths, since
For z → ∞ :
⎛ 2μ K ⎞ z⎟ → 0 exp⎜ − R ⎠ ⎝
(2.15)
Accordingly, the asymptotic shear stress at the wall is
τw =
γR
(2.16)
2
14
Thus, at infinite depths the weight of the material is carried by the wall alone and stresses resulting from any top surcharge decay to zero. However, this means in turn that at finite depths stresses cannot be found zero in case of a non-zero surcharge, even in case of a weightless material. Assumptions (1) and (2) are highly arguable from the point of view of continuum mechanics. Equilibrium requires the existence of shear stresses on the horizontal section that have to be zero in the center of the slice because of symmetry. In turn, that indicates that σ z can neither be constant over the cross sections nor a principal stress. Moreover, the Janssen analysis only gives the shear stress at the wall, but cannot say anything about the shear stress distribution inside the cylinder. Nevertheless, recent experiments (Orvalez, Clément 2005) on constrained granular beds show that the Janssen model and its exponential behaviour reflect the average stress behaviour surprisingly well. Several researchers have dealt with the Janssen model. Terzaghi (1943) pointed out the theoretical inconsistency of (1) and (2), but found experimentally that “they are not so important that they cannot be used as the basis for rough estimates”. Walker (1966) attempted to correct the inconsistency of the principal stress assumption (2) by introducing a distribution factor. Cowin (1977) proved that the assumption (1) of constant
σ z is not necessary to obtain the classical Janssen formulae.
2.2.3 K-value determination. The choice of the K-value Eqn. (2.1) is decisive for the accuracy of the Janssen model. Since Janssen, who determined his K-value experimentally, several theoretical
15
methods have been introduced to compute the K-value (Brown, Richards 1966), (Drescher 1991), (Nedderman 1992). The two most important are
1. The Rankine states: It is assumed that the material is in a state of incipient failure everywhere throughout the cylinder. Thus the Mohr’s circle at any point must touch the internal yield locus φ .
Fig. 2.7 Mohr’s circle for incipient failure
For a cohesionless material, we can find from Fig. 2.7 the principal stresses σ 1 and
σ 3 and the corresponding K-values as σ 1 = p (1 + sin φ )
(2.16)
σ 3 = p (1 − sin φ )
(2.17)
16
K active =
σ 3 1 − sin φ = σ 1 1 + sin φ
K passive =
(2.18)
σ 1 1 + sin φ = σ 3 1 − sin φ
(2.19)
Eqns. (2.18) and (2.19) first developed by Rankine (1857) give the lower and upper limit ratio for principal stresses in cohesionless soil. K active is known as Rankine’s coefficient of active earth pressure and K passive as coefficient of passive earth pressure. Assumption (2) does not specify which of the two normal stresses the major principal stress is. However, experimental experience shows that it is reasonable to assume
σ 1 as the axial stress σ z and hence take K active Eqn. (2.18) as the actual K-value.
2. Elastic state of uniform compression Consider an elastic medium that is constrained as shown in Fig. 2.8 without wall friction and subjected to a surcharge together with the kinematic assumption of plane strain. In this elementary case of uniform compression, σ r , σ z and σ θ perpendicular to the plane are principal stresses and thus assumption (2) is true. With the stipulation that the material cannot bulge out in r and θ directions we find from the stress-strain relations Eqns. (3.21) – (3.24) derived in Fig. 2.8 Uniform compression
section 2.1
εr =
1 [σ r − ν (σ θ + σ z )] = 0 E
(2.20)
17
εθ =
1 [σ θ − ν (σ r + σ z )] = 0 E
(2.21)
and can express from these two relations σ r in terms of σ z and ν . Substitution in Eqn. (2.1) gives K in terms of Poisson’s ratio ν for the elastic state of uniform compression
K=
ν
(2.22)
1 −ν
In the elastic solution approach for the “zero-axial-stress” problem in chapter 3, we will use Eqn. (2.22), because we can derive K directly from the elastic environment. In Fig. 2.9, Janssen’s expressions for σ r , σ z and τ w are plotted for the case of the activated alumina granulate with material parameters [ν = 0.32 ; μ = 0.36 ; K = 0.4706 ] and confined in the plastic tube with R = 3.75 cm . We assume a surcharge of
σ 0 = 0.1 N / mm 2 . The material is idealized as weightless; therefore the asymptotic value is zero. Janssen analysis 0.1 Janssen axial stress Janssen radial stress Janssen shear stress
0.09 0.08
Stress [N/mm²]
0.07 0.06 0.05 0.04 0.03 0.02 0.01 0
0
5
10
15 20 25 Depth from loading end z [cm]
30
Fig. 2.9 Example of Janssen formulae
18
35
40
2.3 Force Chains and Arching
Although we have decided to treat the cylinder as a homogeneous material and apply the principles of continuum mechanics to find stresses and displacements, it is important that we always keep in mind the highly complex discrete nature of granular materials. Therefore, we may look briefly into the two key concepts of discrete force propagation: Force chains and arching.
2.3.1 The cannon ball example. To demonstrate the complexity of what might actually happen inside our cylinder, we consider the simple example of a cannon ball stack subjected to a force P (Duran 1999). Idealized balls are stacked in hyperstatic equilibrium, involving six points of contact for each
Fig. 2.10 A possible actual force scenario and corresponding contact points
ball. As a result, we can afford to randomly lose several contacts - because of tiny variations in shapes and diameters of real balls. This situation leads to an almost random spatial distribution of contact points (see Fig. 2.10). Given the large number of ways to
Fig. 2.11 Force distribution for idealized balls and stresses for elastic model
dispose of the nonessential contacts, 19
the force distribution shown is only one of many possibilities. The forces at the bottom of a cannon ball stack are thus not really predictable, since contacts are expected to differ totally from stack to stack. If the cannon balls are assumed perfectly smooth and absolutely identical, a triangular solution can be easily found. If we apply spatial homogenization – analogous to the continuum assumption for the “zero-axial-stress” problem, we find the parabolic elasticity solution which is similar to the one of the idealized balls (see Fig. 2.11). The cannon ball example indicates the significant simplification of a continuum approach in granular materials. Continuum stresses and displacements are incapable of predicting the actual behaviour of a material sample. They can be at best predictions that come close to averaged real-world stresses and displacements of a range of samples. Relatively good accordance with experimental results is expected in the case of idealized spheres.
2.3.2 Force chains. The disorder of contacts in the packing structure of the grains leads to the formation of force chains. Forces do not propagate continuously, but along discrete paths that are determined by particle contacts (see Fig. 2.12). The idea of force chains was discovered experimentally for the first time by Dantu (1957). He Fig. 2.12 Force chain
used glass cylinders stacked up in a transparent container with their axes all lined up perpendicu-
20
lar to the plane in Fig. 2.13. Upon compressing the cylinders with a piston, he lighted up the entire assembly from behind and used crossed polarizers on either side of the container. Due to the photo-elastic effect, the luminosity of the cylinders is proportional to the stress carried and a stress pattern can be observed (see also Figs. 2.13 and 2.14). From these pictures, we can clearly see the tendency of the force chains to redirect horizontal forces to the sides. This observance gave rise to the idea of arch effects inside granular materials.
Fig. 2.13 Stress pattern in a 2-D granular material under compression after Dantu (1957)
21
Fig. 2.14 Photo-elastic experiment done by Travers et al. (1987)
2.3.3 Arching. From the point of view of structural engineering, an arch is a special structure that is able to transfer loads to its outer extremities. It can be built of loose blocks of stone, because the arch works under normal compression without the need of developing large bending moments. By looking at Fig. 2.15, we can
Fig. 2.15 Arch structures are stable and endurable (Aqueduct of Maintenon)
imagine easily that force chains inside a granular body can work analogously to structural arches. In terms of the “zero-axial-stress” problem, this idea provides an intuitive explanation why the granular material can stay in the tube without being supported from the bottom. Force chains redirect horizontal forces to the tube walls where they are taken by wall fric-
Fig. 2.16 Stress pattern from a Contact Dynamics simulation by Radjai (1997)
tion. At some distance there is a final arch that carries all residual forces including the weight of the material. However, a general theory that is able to model the phenomena of force chains and arching reliably is still lacking, although various approaches have been proposed to es-
22
tablish a general theoretical framework. Besides different continuum approaches (Horne, Nedderman 1976), (Cantelaube, Goddard 1997), (Goddard 1998), (Savage 1998), (Claudin 1999), researchers mainly focus on statistical methods that are based on the framework of statistical mechanics (Liu, Nagel, Coppersmith et al. 1995), (Claudin, Bouchaud 1998) and computational methods that take into account the behaviour of any particle in a structure. Well-known examples for the latter are the Discrete Element Method (Cundall, Strack 1979), (Sakaguchi, Ozaki 1993) and the Contact Dynamics Method (Radjai 1997) (see Fig. 2.16).
23
2.4 Geotechnical Classification Parameters
Geotechnical engineers have gathered extensive experience in describing analytically and experimentally the behaviour of granular aggregates in the form of soils. We will review briefly geotechnical classifications in terms of grain size distribution and weight-volume relationships (Das 2000) that are used in chapter 4 and 5, and find the respective parameter data for alumina granulate.
2.4.1 Grain-Size Distribution Soil consists of particles whose size is usually spread over a wide range. Depending on the predominant particle size in the aggregate, soils are called gravel, sand, silt or clay. The corresponding soil-separate-size limits according to the American Association of State Highway and Transportation Officials (AASHTO) are shown in Fig. 2.17. The particle size distribution affects the mechanical properties such as elastic moduli or shear strength considerably.
Gravel 76.2 – 2.0
Sand 2.0 – 0.075
Silt 0.075 – 0.002
Clay < 0.002
Fig. 2.17 Particle size classification according to the American Association of State Highway and Transportation Officials (AASHTO)
The determination of the size range of particles is usually done by a sieve analysis. It consists of shaking the soil sample through a set of sieves that have progressively smaller openings. After the shaking period, the mass of soil retained in each sieve is measured. The results are presented on a semi-logarithmic plot known as particle- or grain-size distribution curve. 24
Anandakumar (2005) determined the particle size distribution for the original alumina granulate sample by a sieve analysis, whose distribution curve is shown in Fig. 2.18. According to this data, alumina granulate can be well compared to coarse-grained dry sand in terms of particle-size distribution. For the evaluation of material parameters and the mobilization testing in chapter 4, we take out all particles smaller than 1.7 mm by sieving the original alumina granulate used in (Anandakumar 2005) through a sieve sized 1.7 mm (sieve number 12 according to U.S. Standard Sieve Sizes). Thus, the grain size is practically restricted to a small range from 1.7 to 2 mm (see Fig. 2.18). Compared to the original sample, this gives a far more homogeneous aggregate, whose mechanical behaviour is expected to be fairly consistent in terms of experimental repeatability and small scattering of the results.
100
Percent Passing
80
Ottawa Sand
60
1.7 mm point
40
Alumina Granulate
20 0 0.01
0.1
1
Particle size (mm) Fig. 2.18 Particle size distribution of alumina granulate (dots) and Ottawa sand (triangles) as shown in (Anandakumar 2005)
25
10
2.4.2 Weight-volume Relations Soils are three-phase systems consisting of soil solids, water and air. Weightvolume relationships between them can be developed by separation of the phases (see Fig. 2.19).
Fig. 2.19 Parameters for weight-volume relations and separation of the three phase system
The total volume of a given soil sample can be expressed as
V = V s + Vv = V s + V w + V a
(2.23)
where the indices s, v, w, a stand for solids, voids, water and air, respectively. Since the alumina granulate is assumed to be completely dry, we neglect water and deal with an aggregate consisting only of solid grains and voids filled with air.
26
The volume relationships commonly used to classify dry soils in geotechnical engineering are void ratio and porosity. They are defined as follows:
Void Ratio e :
e=
Vv Vs
(2.24)
Porosity n :
n=
Vv V
(2.25)
The relationship between void ratio and porosity reads
e=
n 1− n
and
n=
e e +1
(2.26)
The (dry) unit weight is the weight of soil per unit volume and is defined as
Unit Weight γ :
γ =
Ws Ws = V V s + Vv
(2.27)
where Ws is the weight of the solid particles in a sample with total volume V .
2.4.3 Unit Weight and Void Ratio in Alumina Granulate To measure unit weight and void ratio in alumina granulate, we fill Ws = 100 g of material that was previously sieved according to section 2.4.1, in the PMMA tube (∅ 7.5 cm) through a funnel and homogenize the packing by slightly shaking the tube. We repeat this procedure 10 times and find the average depth of the granular column
l = 3,28 cm from which V = 145 cm 3 can be computed.
27
Anandakumar (2005) identified the particle density as γ p = 1,8 g / cm 2 from which we can compute Vs = 56 cm 3 for the case of Ws = 100 g .
Accordingly we find for the weight-volume relationships for our activated alumina:
Void Ratio e :
e = 1.56
(2.28)
Porosity n :
n = 0.61
(2.29)
Unit Weight γ :
γ = 0.69 g / cm 3
(2.30)
In the testing chapter 4, we will guarantee equivalent packing of the alumina granulate in all experiments by providing samples with consistent unit weight and void ratio Eqns. (2.28) – (2.30). This is accomplished by controlling the unit weight in the same way as described above.
28
Chapter 3: An Analytical Elastic Solution Approach by the Fourier-Bessel Potential
3.1 Fundamentals of Theory of Elasticity
The theory of elasticity provides a tool for describing stresses and deformations in engineering materials. Before dealing with the advanced theory of elastic potentials to derive an analytical solution to the “zero-axial-stress” problem we may begin by recalling the derivation of basic formulae for the simple two-dimensional case of plane strain. We will also introduce the concepts of axisymmetry and cylindrical coordinates.
3.1.1 Basic assumptions of linear elasticity It is of major importance that we keep in mind basic idealizations and assumptions of elasticity (Timoshenko, Goodier 1970) before going into the mathematical formulation. It is this simple set of assumptions that will help us later on to distinguish between elasticity and plasticity and to estimate the capabilities of an elastic model in describing stresses and displacements in a granular medium:
(1)
Elasticity is a continuum theory and does not consider discrete particle structures. It is assumed that the matter of an elastic body is homogeneous and continuously distributed over its volume. This implies that the smallest element cut from the body possesses the same specific physical properties as the body itself. Further-
29
more, these properties are assumed isotropic, which means they are the same in all directions.
(2)
All deformations produced by external forces disappear after removal of these forces. This implies that elastic bodies resume their initial form when unloaded.
(3)
Deformations are small and the corresponding displacements do not affect substantially the action of the external forces. Calculations can be based on initial dimensions of the body. Superposition can be applied.
(4)
The problem is restricted to elasto-statics.
For the following sections, we take the kinematic assumption of plane strain
( ε z = γ x z = γ y z = 0 ) . Details on the following can be found in (Fung 1965), (Timoshenko, Goodier 1970), (Saada 1974), (Cook, Young 1999).
3.1.2 Equilibrium and Kinematic Equations We consider 2-D stresses σ x , σ y , τ xy , that act on a differential plane element as shown in Fig. 3.1. For static equilibrium to prevail, forces acting on the differential element must sum to zero in both x and y directions.
∂ σ x ∂ τ xy + + Bx = 0 ∂x ∂y ∂ τ xy ∂x
+
∂σ y ∂y
(3.1)
+ By = 0
(3.2)
30
These are the differential equations of equilibrium for plane strain. A physically possible state of stress satisfies the differential equations of equilibrium at every point.
Fig. 3.1 Stresses on a differential slice in plane strain
Fig. 3.2 Displacements and strains
In the next step, we consider 2-D displacements with the components u ( x, y ) parallel to the x-axis and v( x, y ) parallel to the y-axis. As shown in Fig. 3.2, strains ε x , ε y and γ xy are expressed in terms of u and v for small deformations as follows
εx =
∂u ∂x
(3.3)
εy =
∂v ∂y
(3.4)
γ xy =
∂u ∂v + ∂y ∂x
(3.5)
These are the strain-displacement relations for plane strain.
31
3.1.3 Constitutive Relations, Compatibility and Boundary Conditions For an isotropic elastic material, the components of stress and the components of strain can be related by only two independent material parameters. The most common set of parameters are Young’s modulus E and Poisson’s ratio ν .
[
]
(3.6)
[
]
(3.7)
εx =
1 σ x − ν (σ y + σ z ) E
εy =
1 σ y − ν (σ x + σ z ) E
γ xy =
1 τ xy G
(3.8)
G is the shear modulus or modulus of rigidity and can be obtained by
G=
E 2 (1 + ν)
(3.9)
These are the stress-strain relations for plane strain. Keep in mind that
σ z = ν (σ x + σ y ) ≠ 0 .
The three strain components (3.3) – (3.5) are expressed by two functions u and v ; hence they cannot be taken arbitrarily, and an additional relation between the strains must exist. Differentiating (3.3) twice with respect to y , (3.4) twice with respect to x and (3.5) once with respect to x and y gives
2 2 ∂ 2 ε x ∂ ε y ∂ γ xy = + ∂x ∂ y ∂x2 ∂ y2
(3.10)
32
This is the compatibility equation for a plane problem. It ensures that displacements are continuous and no cracks or overlapping of material points appear. If we substitute Eqns. (3.6) – (3.8) in (3.10), we find
∂ 2 τ xy ∂2 ∂2 σ x − ν (σ y + σ z ) + 2 σ y − ν (σ x + σ z ) = 2 (1 + ν) ∂x ∂ y ∂ y2 ∂x
[
]
[
]
(3.11)
Substitution of σ z = ν (σ x + σ y ) and Eqns. (3.1) and (3.2) yields
⎛ ∂2 ∂2 ⎜⎜ 2 + 2 ∂y ⎝ ∂x
⎞ 1 ⎟⎟ (σ x + σ y ) = − (1 − ν) ⎠
⎛ ∂B x ∂B y ⎜⎜ + ∂y ⎝ ∂x
⎞ ⎟⎟ ⎠
(3.12)
This is the compatibility equation in plane strain expressed in stress components for the general case of body forces.
Stresses and displacements do not only have to satisfy equilibrium and compatibility, but also prescribed values at the boundary. Accordingly, boundary conditions in elastostatics appear as surface tractions or displacements (see Fig. 3.3). We may prescribe
Fig. 3.3 Stress and displacement boundary conditions
components of each vector on the same boundary region. However, the same component of the stress traction vector and the displacement vector cannot be specified at the same point on a boundary (Little 1973). We may keep the latter in mind for it becomes important in section 3.4. For further general information on boundary conditions we refer to (Fung 1965), (Little 1973) and (Saada 1974). 33
3.1.4 Potential Functions The equations of equilibrium (3.1) and (3.2) and the compatibility equation (3.12) together with the boundary conditions provide a system that is sufficient for the complete determination of the stress distribution in elementary 2-D problems. It can be solved by integration of Eqns. (3.1), (3.2) and (3.12) and evaluating integration constants with boundary conditions (Timoshenko, Goodier 1970). A more powerful method of solving this system of equations, that enables us to solve more complex questions such as the “zero-axial-stress” problem, is the introduction of a potential function φ . Assume that body forces are zero. It may be verified by substitution that Eqns. (3.1) and (3.2) are satisfied if we assume the following relations.
σx =
∂ 2φ ∂ y2
σy =
∂ 2φ ∂x2
τ xy = −
∂ 2φ ∂x ∂ y
(3.13)
The true solution of the problem is that, which satisfies also the compatibility equation (3.12). Substitution of Eqns. (3.13) into (3.12) yields
∂ 4φ ∂ 4φ ∂ 4φ + + 2 =0 ∂x4 ∂x2∂ y2 ∂ y4
(3.14)
Thus, the solution of a 2-D problem reduces to finding a solution of Eqn. (3.14) that satisfies the boundary conditions of the problem.
34
3.1.5 Axisymmetry and Cylindrical Coordinates The kinematic assumption of axisymmetry is based on rotational symmetry of both boundary conditions and geometry about an axis. In this case, all field equations in cylindrical coordinates r , θ , z are independent of θ . It is shown in (Fung 1965), (Timoshenko, Goodier 1970) that for axisymmetry uθ , γ rθ , γ θ z , τ rθ , τ θ z are zero.
For effective use in cylindrical problems, the relevant field equations (3.1) – (3.9) need to be transformed in cylindrical coordinates. We assume axisymmetric conditions.
Kinematic equations. Consider first the cross sectional plane as shown in Fig. 3.4 with displacement u r only. If u r is the radial displacement of the side ad of the element abcd, the radial displacement of
Fig. 3.4 Polar coordinates r, θ
the side bc is u r + (∂ u / ∂ r ) dr . Hence
εr =
∂ ur ∂r
(3.15)
Although uθ is zero, there is a non-zero strain component in tangential direction. With radial displacement u r , the new length of the arc ad is (u r + r ) dθ and the tangential strain is therefore
εθ =
ur r
(3.16)
35
In the rotational plane, there are both displacements u r , u z . Strains can be found analogous to the plane strain case.
εz =
∂ uz ∂z
(3.17)
γ rz =
∂ ur ∂ u z + ∂z ∂r
(3.18)
Consider a differential element in the rotational plane as shown in Fig. 3.5. Taking equilibrium in radial and axial directions and assuming that there are no body forces, we arrive at the following differential equations of equilibrium (Fung 1965) Timoshenko, Goodier 1970).
∂ σ rr ∂ τ rz σ rr − σ θθ + + =0 r ∂r ∂z
(3.19)
∂ τ rz ∂ σ zz τ rz + + =0 ∂r ∂z r
(3.20)
Stress-strain relations only require mutually perpendicular directions (Cook, Young 1999), which is satisfied by
r, θ , z . Fig. 3.5 Stresses on a differential element in cylindrical coordinates
Accordingly, they have the same form as before.
36
εr =
1 [σ r − ν (σ θ + σ z )] E
(3.21)
εθ =
1 [σ θ − ν (σ r + σ z )] E
(3.22)
εz =
1 [σ z − ν (σ r + σ θ )] E
(3.23)
γ rz =
2 (1 + ν) τ rz E
(3.24)
The same relations (3.15) – (3.24) can also be obtained from the equations in Cartesian coordinates (3.1) – (3.8) by means of coordinate transformation (Fung 1965) or be derived from the general case of curvilinear coordinates (Saada 1974).
37
3.2 Love’s Potential Function for Axisymmetry
The set of relations for plane problems Eqns. (3.13) are not valid for cylindrical coordinates. These equations have to be derived anew for the axisymmetric cylindrical case. In 1906, Love introduced a potential function Z for considering problems of solids of revolution under the influence of axisymmetric loadings. In the following, we present a brief overview on how the stress and displacement relations can be obtained. Details can be found in (Love 1944), (Timoshenko, Goodier 1970). The derivation is based on the strain-displacement (3.15) – (3.18), the equilibrium (3.19) – (3.20) and the stressstrain relations (3.21) – (3.24). In the 1930s, Galerkin and Papkovich developed three potential functions for the general problem of elasticity, commonly known as the Galerkin vector. The derivation of Love’s potential function for axisymmetry from this general theory can be found in (Fung 1965), (Little 1973).
3.2.1 Derivation of the Stress Relations In analogy to the plane strain case Eqn. (3.13), we assume a potential function φ and postulate that
τ rz = −
∂ 2φ ∂r ∂ z
(3.25)
Substituting (3.25) in Eqn. (3.20) and integrating with respect to z yields
∂ 2φ 1 ∂ φ σz = 2 + r ∂r ∂r
(3.26)
38
Note that adding an integration constant is not required in this case, since it can be readily included in the potential function. This will also hold for all other integrations throughout this chapter. From equations (3.15) and (3.16) we find that
εr =
∂ (r ε θ ) ∂r
(3.27)
With Eqns. (3.21) and (3.23) we can rewrite Eqn. (3.27)
σ r − ν (σ θ + σ z ) =
∂ [(σ θ − ν σ r − ν σ z ) r ] ∂r
(3.28)
Using the product rule for differentiation on the right side of Eqn. (3.28), we find
(1 + ν ) (σ r − σ θ ) =
∂ [σ θ − ν (σ r + σ z )] r ∂r
(3.29)
Now we introduce a subsidiary function ψ by the equation
σr =
∂ 2φ +ψ ∂z2
(3.30)
Rewriting Eqn. (3.29) by substitution of Eqns. (3.19) and (3.30) yields
(1 + ν )
∂ψ ∂ [σ θ − ν (σ r + σ z ) ] = 0 + ∂r ∂r
(3.31)
If we integrate Eqn. (3.31) with respect to r and substitute Eqns. (3.26) and (3.30), we get
39
⎛ ∂ 2φ ⎞ ⎛ ∂ 2φ 1 ∂ φ ⎞ ⎟ σ θ = (1 + ν ) ψ − ν ⎜⎜ 2 + ψ ⎟⎟ − ν ⎜⎜ 2 + r ∂ r ⎟⎠ ⎝ ∂z ⎠ ⎝ ∂r
(3.32)
Eqn. (3.32) can be expressed simply as
σ θ = ν ∇ 2φ − ψ
(3.33)
with the axisymmetric Laplacian operator in cylindrical coordinates.
∇2 =
1 ∂ ∂2 ∂2 + + ∂ r2 r ∂ r ∂ z2
(3.34)
The functions φ and ψ are not independent of each other. To find the relation between them, we use the compatibility condition. In general, this equation does not need to have the form of Eqn. (3.14), but can be any additional relation that ensures the compatibility of the three strain fields Eqns. (3.15), (3.17) and (3.18). A simple compatibility condition can be found as follows. We rewrite Eqn. (3.16) with Eqn. (3.22)
ur =
r [σ θ − ν (σ r + σ r )] E
(3.35)
Integration of Eqn. (3.31) and substitution into (3.35) yields
u r = −(1 + ν )
rψ E
(3.36)
Using Eqns. (3.18) and (3.24)
τ rz =
2 (1 + ν ) ⎛ ∂ur ∂u z ⎞ ⎜ ⎟ + E ⎜⎝ ∂ z ∂ r ⎟⎠
(3.37) 40
we can find a first expression for the derivative of u z with respect to r by inserting Eqns. (3.25) and (3.26) into (3.37).
∂u z 2 (1 + ν ) ∂ 2φ (1 + ν ) ∂ψ =− + r ∂r E ∂r ∂ z E ∂r
(3.38)
We obtain a second expression for the derivative of u z with respect to z by substitution of Eqns. (3.26), (3.30) and (3.33) into Eqn. (3.23)
∂u z (1 + ν ) = E ∂z
⎞ ⎛ ∂ 2φ 1 ∂φ ⎜⎜ 2 + − ν ∇ 2φ ⎟⎟ r ∂r ⎠ ⎝ ∂r
(3.39)
We can find the compatibility condition, which gives a relation between φ and ψ , by integrating Eqns. (3.38) and (3.39) with respect to r and z respectively, and equating the results.
(1 − ν )
∂ 2 ∂ 3φ ∂ 2ψ ∇ φ+ = r ∂r ∂r ∂ z 2 ∂z2
(3.40)
It turns out to be convenient later on if we introduce a new function λ at this point by means of the relation
rψ =
∂ φ ∂λ + ∂r ∂r
(3.41)
With Eqn. (3.41) we find from the compatibility condition (3.40) the simpler relation between φ and λ .
41
∂ 2λ ∂ 2 = (1 − ν ) ∇φ 2 ∂r ∂z
(3.42)
With this relation (3.42) at hand, we now go back to the stress components Eqns. (3.26), (3.30) and (3.33). We want to express these relations in terms of only one single potential function Z. Therefore, we try to replace all other subsidiary potentials by assuming the following new potential
Z=
∂ [φ + λ ] ∂z
(3.43)
For σ r , inserting Eqns. (3.41) and (3.42) into (3.30) yields
σr =
∂ ⎛ ∂ 2φ 1 ∂φ 1 ∂Z ⎞ ⎜⎜ 2 + ⎟ ++ ∂z ⎝ ∂z r ∂r r ∂ r ⎟⎠
σr =
∂ 2φ ∂ 2 Z ∂ 2 Z ⎞ ∂ ⎛ 2 ⎜⎜ ∇ Z − 2 − 2 − 2 ⎟⎟ ∂z ⎝ ∂r ∂z ∂r ⎠
∂ ⎡ 2 ∂2Z ⎤ σr = ⎢ν ∇ Z − 2 ⎥ ∂z ⎣ ∂z ⎦
(3.44)
In the same way we can find
σθ =
∂ ⎡ 2 1 ∂Z⎤ ⎢ν ∇ Z − ⎥ ∂z ⎣ r ∂r ⎦
(3.45)
σz =
∂2 Z ⎤ ∂ ⎡ 2 ( 2 ν ) Z − ∇ − ⎢ ⎥ ∂z⎣ ∂ z2 ⎦
(3.46)
Substitution of Eqns. (3.44) and (3.45) into Eqn. (3.19) brings about
42
τ rz =
∂2 Z ⎤ ∂ ⎡ 2 ( 1 ) Z − ν ∇ − ⎢ ⎥ ∂r ⎣ ∂ z2 ⎦
(3.47)
It can now be verified by substitution that relations (3.44) – (3.47) satisfy the equilibrium equations (3.19) and (3.20) only if
∇ 2∇ 2 Z = 0
(3.48)
Eqn. (3.48) is called a biharmonic equation, and its solution is called a biharmonic function. For the solution of our problem, we can assume any potential function Z that satisfies Eqn. (3.28). Equilibrium and compatibility conditions are automatically satisfied by Z inside the volume of the body. Z is generally referred to as Love’s potential function or Love’s strain function. A compendium of the history of the biharmonic problem together with various examples and solutions can be fond in (Meleshko 2003).
3.2.2 Derivation of the Displacement Relations We want to find the displacements u r , u z corresponding to the stresses expressed in Eqns. (3.44) – (3.47). Substitution of these relations into Eqn. (3.35) yields
ur =
r [σ θθ − ν (σ rr + σ zz )] = E
− 1 +ν ∂ 2 Z E ∂r ∂ z
(3.49)
For u z , we know ∂ u z / ∂ z from (1.17) and ∂ u z / ∂ r from (3.18). Substituting Eqn. (3.17) into (3.23) yields with Eqns. (3.44) - (3.46)
E
∂ uz ∂2 Z ⎤ ∂ ⎡ 2 2 2 ( 1 ) ( 1 ) Z = σ zz − ν (σ rr + σ θθ ) = − ν ∇ − + ν ⎢ ⎥ ∂z ∂z ⎣ ∂2 z2 ⎦
43
(3.50)
Integration with respect to z brings about
⎡ ∂2 Z ⎤ E u z = (1 + ν) ⎢2 (1 − ν 2 ) ∇ 2 Z − 2 ⎥ + f (r ) ∂ z⎦ ⎣
(3.51)
Substituting (3.47) and (3.49) into (3.18) yields
G
∂ uz ∂u 1 ∂2 Z ⎤ ∂ ⎡ 2 ( 1 ) ν Z = τ rz − G r = − ∇ − ⎢ ⎥ 2 ∂2 z ⎦ ∂z ∂z ∂r ⎣
(3.52)
If we recall that G = E / 2 (1 + ν ) , integration with respect to r brings about
⎡ ∂2 Z ⎤ E u z = (1 + ν) ⎢2 (1 − ν) ∇ 2 Z − 2 ⎥ + g ( z ) ∂ z⎦ ⎣
(3.53)
Since (3.51) and (3.53) must be equal, f (r ) and g ( z ) must be equal as well. Since
f (r ) is a function of r only and g ( z ) a function of z only, they must be equal to an arbitrary constant A. This constant corresponds to an axial rigid-body translation in z direction and will be discarded for now on the understanding that it can be restored whenever needed to meet boundary conditions.
f (r ) = g ( z ) = A
(3.54)
Accordingly, the displacement components can be derived from Z as follows:
2 G ur = −
∂2 Z ∂r ∂z
2 G u z = 2 (1 − ν) ∇ 2 Z −
(3.55)
∂2 Z ∂2 z
(3.56)
44
3.3 Fourier-Bessel Series Approach
We developed a set of equations (3.44) – (3.47) and (3.55) - (3.56) that can be applied to a potential function Z to find stresses and displacements in axisymmetric cases. Equilibrium and compatibility are automatically satisfied at any point inside the cylinder, if Z satisfies the biharmonic condition
∇ 2∇ 2 Z = 0
(3.57)
Thus, we reduced our problem to finding an appropriate biharmonic potential function Z that can be adjusted to a set of boundary conditions. We will now give a short overview of the derivation of the Fourier-Bessel potential for Z and then compute stress and displacement expressions for this elastic cylinder problem. Details on the derivation can be found in (Pickett 1944) and (Little 1973). Furthermore, we will discuss some mathematical background on Fourier transform methods that will be needed in chapter 3.4.
3.3.1 The Fourier-Bessel Potential for Z Little (1966, 1973) gives a general solution for the biharmonic equation (3.57) for zero body forces in the following form
Z = [ A J 0 (γ r ) + F Y0 γ z (γ r ) + ∂ / ∂ r {C J 0 (γ r ) + D Y0 (γ r )} (3.58)
+ Eγ z J 0 (γ r ) + F γ z Y0 (γ r )] e
±γ z
where J 0 and Y0 are the Bessel functions of order zero of the first and second kind, respectively. γ is called the separation constant and can be real or imaginary. 45
The potential function Z for the elastic cylinder problem is derived from (3.58) as follows. The coefficients of the Bessel functions of the second kind Y0 are chosen to be zero. If γ is imaginary, we obtain trigonometric behaviour in the z -direction and modified Bessel functions I 0 and I 1 in the r -direction. The former is easily demonstrated by application of the Euler formula
e i z + e − i z = cos z + sin z
(3.59)
to the factor e ± γ z of Eqn. (3.58). We drop the sine part and consider only the cosine part (see also section 3.5 “Aliasing and truncation error”). Thus
Z = cos(γ z ) [ A I o (γ r ) + Bγ r I 1 (γ r )]
(3.60)
where I n (γ r ) is the modified Bessel function, which is defined as
⎛ 1 ⎞ I n ( x) = exp⎜ − n π i ⎟ J n (i x) ⎝ 2 ⎠
(3.61)
Besides condition (3.57), Z must offer enough flexibility in terms of its constants
A, B , which have to be determined such that Z meets the boundary conditions. For the application of Fourier transforms, we need Z to be stated in the form of a series of sines and cosines. This kind of series is called Fourier series [see Eqn. (3.75)]. Since we adopted the framework of linear elasticity, we can choose n → ∞ different separation constants γ n for Eqn. (3.60) and superpose the n resulting equations to a final potential function Z in the form
∞
Z = ∑ cos(γ n z ) [An I o (γ n r ) + Bn γ n r I 1 (γ n r )] n =1
46
(3.62)
which is called the Fourier-Bessel potential in the following. It may be verified by substitution that the Fourier-Bessel potential Eqn. (3.62) satisfies the biharmonic condition (3.57). A rigorous and complete mathematical derivation of Eqn. (3.62) and a discussion of the special Bessel functions involved can be found in (Little 1966, 1973). The potential function Z in this particular form was first proposed by Pickett (1944).
3.3.2 Stress and Displacement Expressions Application of Eqns. (3.44) – (3.47) and (3.55) – (3.56) to the Fourier-Bessel potential (3.62) yields stress and displacement expression for the elastic cylinder
σ r (r , z ) =
∞
∑α
3 n
n =1
σ θ (r , z ) =
σ z (r , z ) =
∞
αn2
n =1
r
∑
⎧⎪ ⎡ I (α r ) ⎤ sin(α n z ) ⎨ An ⎢ I 0 (α n r ) − 1 n ⎥ + αnr ⎦ ⎪⎩ ⎣ Bn [(1 − 2ν ) I 0 (α n r ) + α n rI 1 (α n r )] }
(3.63)
sin(α n z ){An I 1 (α n r ) + Bn (1 − 2ν ) I 0 (α n r )}
(3.64)
∞
∑ −α
3 n
sin(α n z ) {An I 0 (α n r ) +
n =1
τ rz (r , z ) =
Bn [α n rI 1 (α n r ) + 2(2 − ν ) I 0 (α n r )] }
∞
∑α
3 n
cos(α n z ) {An I 1 (α n r ) +
n =1
u r (r , z ) =
(1 + ν ) E
Bn [α n rI 0 (α n r ) + 2(1 − ν ) I 1 (α n r )] }
∑ sin(α z )[A α ∞
n
n
2 n
I 1 (α n r ) + α n rBn I 0 (α n r )
n =1
47
2
]
(3.65)
(3.66)
(3.67)
u z (r , z ) =
∑ cos(α z ) [A α
(1 + ν ) E
∞
n
n
n =1
2 n
{
I 0 (α n r ) +
}]
(3.68)
Bn 4α n (1 − ν ) I 0 (α n r ) + α n rI 1 (α n r ) 2
2
3.3.3 Functional Orthogonality
r
r
Consider two orthogonal vectors u ⊥ v in 3-D space. Recall the inner product (or dot product) that yields
r r r r u , v = 0, if u ≠ v
(3.69)
r r r r u , v = 1, if u = v
(3.70)
The functional analogy to the inner product of two vectors is the inner product of two functions f and g , which is defined as
b
f , g ≡ ∫ f g dx
(3.71)
a
b
f , g = ∫ f g dx = 0
(3.72)
a
b
b
a
a
f , f = g , g = ∫ f f dx = ∫ g g dx = 1
(3.73)
Functions f , g have to be orthonormal, if Eqn. (3.73) is expected to yield 1. Analogously to vectors, we can normalize functions by dividing them by
b
N = ∫ f f dx
(3.74)
a
48
3.3.4 Fourier Series and Finite Fourier Transforms The basic idea of a Fourier series is to express an arbitrary periodic function f (z ) as infinite series of sines and cosines. The function has greatest period 2 L .
f ( z ) = b0
∞ 1 1 1 ⎛ iπ z ⎞ ∞ ⎛ iπ z ⎞ + ∑ bi cos ⎜ sin ⎜ ⎟ + ∑ ai ⎟ L L ⎝ L ⎠ 2 L i =1 ⎝ L ⎠ i =1
i = 1, 2, ... ∞
(3.75)
For simplicity, we introduce
∞
∞
i =1
i =1
f ( z ) = b0ψ 0C + ∑ biψ iC + ∑ aiψ iS
With
ψ 0C =
(3.76)
1 2L
(3.77)
ψ iC =
1 ⎛ iπ z ⎞ cos ⎜ ⎟ L ⎝ L ⎠
(3.78)
ψ iS =
1 ⎛ iπ z ⎞ sin ⎜ ⎟ L ⎝ L ⎠
(3.79)
Since Eqns. (3.77) – (3.79) constitute an infinite set of orthonormal functions, we can determine the unknown coefficients a i , b0 , bi as follows. Presuming that we know
f ( z ) , we find the coefficients ai by multiplying Eqn. (3.76) by Eqn. (3.79) and integrating over one period
L
∫
−L
L
∞
−L
i =1
L
∞
−L
i =1
L
f ( z ) ψ Sj dz = b0 ∫ ψ 0C ψ Sj + ∑ bi ∫ ψ iC ψ Sj + ∑ ai ∫ ψ iS ψ Sj
(3.80)
−L
The first two terms of the right side cancel out on account of the orthogonality condition (3.70) over the period 2 L . 49
1 ⎛ iπ z ⎞ 1 ⎛ iπ z ⎞ cos ⎜ sin ⎜ ⎟, ⎟ = L ⎝ L ⎠ L ⎝ L ⎠
1 1 ⎛ iπ z ⎞ , sin ⎜ ⎟ =0 2L L ⎝ L ⎠
(3.81)
For the third term of the right side, we find
L ⎧ 0, i ≠ j 1 ⎛ iπ z ⎞ 1 ⎛ iπ z ⎞ ψ Sj ψ iS dz = ⎨ sin ⎜ , sin = ⎟ ⎜ ⎟ ∫ L ⎝ L ⎠ ⎝ L ⎠ L ⎩ 1, i = j −L
(3.82)
From Eqn. (3.82), we find with Eqns. (3.72) and (3.73)
L
a1 ⋅ 0 + a 2 ⋅ 0 + ... + a j ⋅ 1 + ... =
∫ f ( z) ψ
S j
dz
(3.83)
−L
Within the power series only the jth term remains and we can determine constants
a1 , a 2 ... a j ... by repeating this procedure with ψ 1S ,ψ 2S ...ψ Sj ... . Similarly, constants b0 , bi can be found and we can summarize L
aj =
∫ f ( z) ψ
S j
dz
(3.84)
C 0
dz
(3.85)
C j
dz
(3.86)
−L L
b0 =
∫ f ( z) ψ
−L L
bj =
∫ f ( z) ψ
−L
This method for obtaining constants a i , b0 , bi is known as the finite Fourier transform.
50
3.4 Boundary conditions
We obtained from the Fourier-Bessel potential a set of equations (3.63) – (3.68) for stresses and displacements in the elastic cylinder. We will now establish appropriate boundary conditions on the cylindrical surface and then adjust those stress and displacement equations to these boundary functions. For this, finite Fourier transform methods are used.
3.4.1 Basic Requirements We may begin by recalling the “zero-axialstress” problem from the introduction: Is it possible in a laterally constrained elastic cylinder that is pushed from one end, to find axial stresses σ z = 0 at some finite depth, if we impose Coulomb’s friction over the radial surface (see Fig. 3.6)? This question comprises the complete set of boundary statements. Note that we don’t specify σ z = 0 explicitly, but find this condition implicitly – in case it is possible – by satisfying Eqns. (3.87) – (3.89) and computing Fig. 3.6 “Zero-axial-tress”: Problem formulation
σ z at various depths. At z = 0 :
σ z (r ) = σ 0
(3.87)
At r = R :
ur ( z) = 0
(3.88)
51
μ=
τ rz ( z ) σ r ( z)
(3.89)
At the top of the cylinder, Eqn. (3.87) represents the constant pushing stress σ 0 . At the radial boundary, Eqn. (3.88) constitutes the rigid constraint by the tube wall and Eqn. (3.87) represents Coulomb’s friction. However, Eqn. (3.89) is ill-posed (see section 3.1.3). This relation includes σ r and is therefore not allowed to be applied directly (Little 1973), since the displacement component in radial direction is already prescribed by Eqn. (3.67). Moreover, if we insert Eqns. (3.63) and (3.66) into Eqn. (3.89), we will find that the resulting relation is ineligible for the application of transform methods. Other valuable information can be found by looking at the Fourier-Bessel potential: Since there are only two sets of constants An , Bn , we can only satisfy two of Eqns. (3.87) – (3.89) by transform methods. Fourier transforms can be applied only for boundary functions f (z ) , we must hence choose Eqns. (3.88) and (3.89) to satisfy An and Bn . Fortunately, it will turn out that Eqn. (3.87) can be satisfied by superposition of a simple uniform compression state on the Fourier transform solution.
3.4.2 The Janssen Shear Stress Assumption As for Eqn. (3.89), only σ r is problematic, but using the shear stress τ rz as a boundary condition is still allowed. Recall the derivation of the Janssen equations (2.12) – (2.14) from section 2.2. We expect these relations to be a fair approximation of the exact solution that we are actually looking for. Janssen uses Coulomb’s condition at the wall of his differential slice to derive his stress relations. In turn, the Janssen shearing stress at the wall inherently satisfies Coulomb’s condition.
52
The basic idea is to specify the Janssen shear stress distribution as a boundary function over the radial surface. Having obtained a solution for the Bessel potential, we can substitute the resulting normal stresses σ rPotential together with the Janssen shear into the Coulomb’s condition and get τ rJanssen z
τ rJanssen z σ rPotential
= μ close ≈ μ
(3.90)
Since Janssen’s shear stress is supposed to be close to the exact solution, only a very small deviation of μ close from the actual μ is expected. Accordingly, we modify the set of boundary functions for the radial surface as
At r = R :
⎛ 2μ K ⎞ z⎟ R ⎠ ⎝
τ rz = μ K σ 0 exp⎜ −
(3.91)
ur ( z) = 0
(3.92)
3.4.3 Application of Transform Methods By means of the finite Fourier cosine transform, we can determine the constants
An , Bn in the Fourier-Bessel potential such that the shear stress expression Eqn. (3.76) adds up to Eqn. (3.92). We begin by choosing the separation constants α n . Recall from section 3.3 that the n cosine functions cos(α n z ) from the Fourier-Bessel potential must form an orthogonal set for the application of transform methods. In other words, we have to assume α n such that cos(α n z ) are orthogonal to each other. Let’s consider first the simpler case of cosine functions in the form of
53
f n ( z ) = cos(n z )
(3.93)
This set of functions is orthogonal on the interval [− 1 / 2 π ; 1 / 2 π ] such that Eqn. (3.82) together with Eqn. (3.78) are satisfied. We can find the following orthogonality condition (3.94) from substituting Eqn. (3.93) into these relations and integrating by parts twice. From Eqn. (3.94) the separation constants n can then be determined.
⎛ 1 ⎞ cos⎜ n π ⎟ = 0 ⎝ 2 ⎠
→
n = 1, 3, ... ∞
(3.94)
If we change the interval to [− L; L ] , it may be verified by substitution that Eqn. (3.95) is equivalent to Eqn. (3.94). Hence we find the separation constants α n as
cos(α n L) = 0
αn =
→
(2n − 1) π 2 L
(n = 1, 2, ..., ∞)
(3.95)
Next, we choose a relation between the constants An and Bn by looking at the boundary condition Eqn. (3.92). Consider the radial displacement expression
u r (r , z ) =
(1 + ν ) E
∑ sin(α z)[A α ∞
n
n
2 n
I 1 (α n r ) + α n rBn I 0 (α n r ) 2
]
(3.96)
n =1
For Eqn. (3.96) to be zero all over the radial cylinder surface [R ; z ], the bracket term has to be zero at r = R , since the sine term is not everywhere zero there. Accordingly, we find
An = − Bn
R I 0 (α n R) I 1 (α n R)
(3.97)
54
Substitution of the relation (3.97) into the stress expression (3.76) yields
τ rz (r , z ) =
∞
∑ n =1
α n 3 cos(α n z ) {− 5α n R I1 (α n r ) I 0 (α n R) +7 I 1 (α n r ) ⋅ Bn 5 I 1 (α n R ) I 1 (α n R) + 5 α n r I 1 (α n r ) I 0 (α n R) }
(3.98)
If we substitute the second boundary function into Eqn. (3.98), we find
τ rz ( R, z ) =
∞
∑ n =1
α n 3 cos(α n z ) {− 5α n R I1 (α n R) I 0 (α n R) + 7 I 1 (α n R) I 1 (α n R) Bn 5 I 1 (α n R) ⎛ 2μ K ⎞ z⎟ + 5 α n R I 1 (α n R ) I 0 (α n R )} = μ K exp⎜ − R ⎝ ⎠
(3.99)
Now we summarize in Eqn. (3.99) all functions other than the cosine to a scalar value
Gn Eqn. (3.100) and bring it to the right side by division. Denoting the Janssen shear stress Eqn. (3.91) by f (z ) , we get for each n
α n3 {− 5α n R I1 (α n R) I 0 (α n R) + 7 I 1 (α n R) I 1 (α n R) Gn = 5 I 1 (α n R ) + 5 α n R I 1 (α n R ) I 0 (α n R )}
(3.100)
1 f ( z ) = Bn cos(α n z ) Gn
(3.101)
Eqn. (3.101) has the form of Eqn. (3.76) with the slight difference that there are only cosine terms and no normalization factor
1 / L . Hence the constants can be obtained
by application of the finite Fourier cosine transform from Eqn. (3.86) as follows: First we normalize the left side of Eqn. (1.101) by
1 / L . Then we multiply both
sides by ψ nC = 1 / L cos(α n z ) and integrate with respect to z from [− L ; L ] . Considering Eqn. (3.83), we find the constants Bn consecutively for n = 1, 2, ..., ∞ by repeating this procedure for each n . The series is truncated after n = 50 . The corresponding con55
stants An can be readily obtained by substitution of the evaluated constants Bn in Eqn. (3.97). We can then insert cos(α n z ) into Eqns. (3.63) – (3.68) and superpose the first 50 series terms to the stress and displacement fields. However, this solution only satisfies the boundary conditions at the radial surface.
3.4.4 Superposition of Polynomial Solutions We are looking for a solution also satisfying Eqn. (3.97) at z = 0 . Recall the FourierBessel expression for axial stress
σ z (r , z ) =
∞
∑ −α
3 n
sin(α n z ) {An I 0 (α n r ) +
n =1
Bn [α n rI 1 (α n r ) + 2(2 − ν ) I 0 (α n r )] }
(3.102)
On account of the sine factor, the axial stress in Eqn. (3.102) must be zero all over the top surface z = 0 . Thus, if we can find a second solution that has uniform compression σ 0 at z = 0 and does not interfere with the radial boundary conditions already satisfied, we can superpose them to the final solution. Little (1973) gives a number of polynomial potentials that can be taken as a Love’s potential function Z. Consider the following polynomial potential.
Z = C ( z 3 − 3 2 z r 2 ) + D (r 2 + z 2 ) z
(3.103)
It may be verified by substitution that Eqn. (3.103) satisfies the biharmonic condition Eqn. (3.48). A discussion on the derivation of Eqn. (3.102) from Legendre polynomials can also be found in (Little 1973). With relations (3.44) - (3.48) and (3.55) – (3.56) we find stresses and displacements as 56
σ r = σ θ = 3 C + (10ν − 2) D
(3.104)
σ z = − 6 C + (14 − 10ν ) D z
(3.105)
τ rz = 0
(3.106)
2 G u r = − 6 C z + (14 − 20ν ) D z
(3.107)
2 G u z = 3 C r − 2D r
(3.108)
and assume the following two boundary conditions to determine C, D
At z = 0 :
σ z (r ) = σ 0
(3.109)
At r = R :
ur ( z) = 0
(3.110)
Fig. 3.7 Uniform compression without shear or friction at the wall
The solution corresponds to a constrained elastic cylinder in uniform compression σ 0 without friction at the wall (see Fig. 3.7). If we superpose this solution to the one obtained from the Fourier-Bessel potential, we find stress and displacement fields that satisfy all three initial boundary requirements Eqns. (3.87) – (3.89). The complete set of computational steps is implemented in a MATLAB routine that can be found in Appendix A.
3.4.5 Eigenfunction Analysis Conclusively, it should be mentioned that a more refined mathematical insight can be provided by the general theory of Eigenfunction analysis for boundary value prob-
57
lems. In view of the frequent appearance in all fields of applied mechanics, we may formulate the elastic cylinder problem briefly in terms of this important mathematical theory. In our cylinder, the coordinate system ( r , z ) is chosen such that the radial boundary of the elastic body corresponds to the coordinate surface ( R, z ) . In other words, the solution on the radial boundary surface involves only one coordinate variable z . Additionally, we can characterize the separation constants α n as a second variable. First, we determined the separation constants by satisfying orthogonality conditions at the boundary edges. Eqn. (3.95) is called characteristic equation and the resulting α n are the corresponding Eigenvalues. The boundary functions, the single series terms of Eqn. (3.98) containing these Eigenvalues, are called Eigenfunctions. It can be shown that a series of these Eigenfunctions, for all possible Eigenvalues at the edges, is able to represent any given boundary function in between. Since trigonometric functions are orthogonal, the constants of the Eigenfunctions can be easily obtained by transform methods as demonstrated. This procedure is therefore called Eigenfunction expansion. The first 3 Eigenfunctions of the Fourier-Bessel potential are shown in Fig. 3.8. One could imagine that a superposition of an infinite sum of these functions could add up to any specified boundary function. This concept is the functional equivalent to the wellknown concept of linear combination from vector algebra. It turns out that the elastic cylinder problem belongs to the general class of Sturm-Liouville problems. Further discussion on Sturm-Liouville problems and their solution by Eigenfunction expansion can be found in (Kreyszig 1999), (Little 1973).
58
-3
20
Eigenfunction expansion
x 10
1. Eigenfunction (alpha=1/2*pi/L) 2. Eigenfunction (alpha=3/2*pi/L) 3. Eigenfunction (alpha=5/2*pi/L) Janssen shear stress (boundary function)
Shear stress [N/mm²]
15
10
5
0
-5
0
5
10 Depth from loading end z [cm]
15
20
Fig. 3.8 First three Eigenfunctions for the Fourier- Bessel potential
59
3.5 Error Analysis and Improvement
Stress and displacement results on the basis of the Fourier-Bessel formulation problem can be obtained from the MATLAB routine in Appendix A and are displayed graphically in section 6.1. We adopt the following setting with material parameters for the case of alumina granulate from section 4.2:
-
Geometry:
Radius R = 3.75cm, Cylinder Length Lcyl = 15cm , Sampling interval [ 0; L = 20cm ]
-
Material:
-
Surcharge:
-
E = 15 ⋅ 10 6 N / mm 2 , ν = 0.32, μ = 0.36 ,weightless
σ 0 = 0.1 N / mm 2 ν Janssen shear ass.: K = = 0.4706 (see section 2.2) 1 −ν
Here, we first give some additional considerations on the Fourier-Bessel potential from a geometric point of view and assess the solution fields in terms of their computational accuracy. Finally, we will discuss the Janssen error resulting from the Janssen shear stress assumption and propose an algorithm that reduces this error source.
3.5.1 Periodicity vs. Monotonic Decrease The nature of the solution fields that are derived from the Fourier-Bessel potential is periodic in the z -direction, since their behaviour is determined by sines or cosines. Thus stresses and displacements are either anti-symmetric (odd functions) or symmetric (even functions) with regard to z = 0 . Orthogonality requires that Eigenvalues α n are chosen such that cos(α n L) = 0 in Eqn. (3.95). Accordingly, solutions fields with cosines
60
or sines are zero at
[ n L ; n = 1, 3, ... ]
or
[ n L ; n = 0, 2, 4, ... ] ,
respectively, and their
period is 2 L . However, the stress and displacement behaviour of the confined elastic medium is not periodic, but exponential and monotonic decreasing with respect to z . Despite the periodicity of the trigonometric series, it is possible to obtain monotonic behaviour within a period (see Fig. 3.9). Let’s look at an interval z ∈ [ 0; L cyl ] for a finite cylinder length. In order to get non-periodic solutions we choose this interval [ 0; L cyl ] to remain within the sampling interval [ 0; L ]. Consider the Janssen shear assumption Eqn. (3.111) and the shear stress derived from the Fourier-Bessel potential Eqn. (3.112) at r = R
⎛ 2μ K R ⎝
τ rz = μ K σ 0 exp⎜ − τ rz ( R, z ) =
∞
∑α n =1
3 n
⎞ z⎟ ⎠
(3.111)
cos(α n z ) {An I 1 (α n R) + Bn [α n R I 0 (α n R) + 2(1 − ν ) I 1 (α n R )] }
(3.112)
We determined the constants An , Bn by means of the finite Fourier transform such that Eqn. (3.112) represents Eqn. (3.111). Both functions are plotted in Fig. 3.9. We clearly see the periodicity of the series and the large deviation from the Janssen exponential function. The series only consists of cosine terms and thus represents an even function symmetric to z = 0 . Nevertheless, within the sampling interval [ 0; L ] the series is able to represent exponential behaviour very accurately.
61
Trigonometric vs. exponential behavior *10-2 Fourier series expansion at r=R Exponential Janssen shear
Shear stress [N/mm²]
2
1
0
-1
-2 -20
0
20 40 60 Depth from loading end z [cm]
80
100
Fig. 3.9 Exponential Janssen shear and its Fourier series representation
Sampling error *10-2 Fourier series expansion at r=R Exponential Janssen shear
Shear stress [N/mm²]
0.5
0.4
0.3
0.2
0.1
0 12
13
14
15 16 17 Depth from loading end z [cm]
Fig. 3.10 Sampling error on z ∈ [12; L=20] for n=20
62
18
19
20
3.5.2 Sampling and Truncation Error However, the series representation is not exact on the sampling interval either. Compare the exponential behaviour of Eqn. (3.111) to the series behaviour of Eqn. (3.112) on [ 0; L ]. Whereas the exponential shear function asymptotically tends towards zero and reaches it only in infinity, the cosine series must be zero already at the end of the sampling interval z = L on account of the orthogonality condition Eqn. (3.72). Therefore, the closer the z in the cosine term approaches L , the more difficult it becomes for the cosine series to model exponential behaviour appropriately (see Fig. 3.10). This error is known as the aliasing, sampling or discretization error (Jerri 1992). It occurs towards the end and beginning of each period and results from the fact that the series periodicity is forced upon a monotonic function. As a result, we have to choose the sampling interval [ 0; L ] larger than the length of the cylinder [ 0; L cyl ] to avoid major influences of the sampling error on L cyl . In Fig. 3.10, the error decays to a reasonably small amount at depths smaller than 15 cm (see Fig. 3.11), which is in accordance with our choice of L cyl = 15 cm . With the sampling error in mind, we can now easily explain the specific choice of only cosine functions for the Fourier-Bessel potential Eqn. (3.60). For modelling an exponential function that gives 1 at z = 0 and a much smaller number at z = L , the cosine function giving 1 at z = 0 and 0 at z = L is more appropriate than the sine function giving 0 at z = 0 and 1 at z = L . We could have added an additional sine term into the Fourier-Bessel potential Eqn. (3.4) obtaining a potential in the form
63
Z =
∞
∑ cos(γ
n
z ) [An I o (γ n r ) + Bn γ n r I 1 (γ n r )]
n =1 ∞
+ ∑ sin( β m z ) [C m I o ( β m r ) + Dm β m r I 1 ( β m r )]
(3.113)
m =1
Eqn. (3.113) is expected to give a much smaller error at the end and beginning of the period which is basically reduced to the well-known Gibbs phenomenon that results from the sudden jump at z = L (Sneddon 1972), (Jerri 1992). But it is far easier for us to deal with the sampling error here than to go through the additional computational expense of Eqn. (3.113). A second error comes from truncation of the infinite Fourier series Eqn. (3.62). We evaluate the constants in Eqn. (3.112) only for a finite n and compute the series as a finite sum of n terms. In Fig. 3.11 we can see the incomplete series expression oscillating around the true Janssen shear stress. The truncation error can readily be interpreted as the difference from the true solution that would have been cancelled out by the truncated terms. It can be easily reduced to any small difference by increasing the number of series terms. Thus, the accuracy of the series representation within the interval depends on the one hand on L and the corresponding sampling error, on the other hand on n and the corresponding truncation error. Additionally, it must be taken into account that these two errors influence each other. Increasing the sampling interval decreases the sampling error. But at the same time, we increase the truncation error, since for the same accuracy on a larger interval, more series terms are needed. Accordingly, we have to increase the number of series terms also, if we want to keep the truncation error constant.
64
Truncation error *10-2 Fourier series expansion at r=R Exponential Janssen shear
Shear stress [N/mm²]
0.8
0.76
0.72
0.68
8
8.5
9 Depth from loading end z [cm]
9.5
10
Fig. 3.11 Truncation error on z ∈ [8; 10] for n=20
Coulomb´s condition mu(close) Janssen/E-Theory FEM
Coulomb's coefficient mu [-]
mu exact = 0.36 0.36
0.355
0.35
0.345 0
5 10 Depth from loading end z [cm]
Fig. 3.12 Coulomb’s coefficient
μ close
for n=50
65
15
3.5.3 The Janssen Error and Improvement Algorithm As indicated in section 3.4, we expect a third error source in the form that the boundary condition Eqn. (3.89)
At r = R :
μ=
τ rz ( z ) σ r ( z)
(3.114)
is not exactly satisfied due to the Janssen shear stress assumption that is expected to differ slightly from the exact shear solution. Instead, we get a Coulomb’s coefficient μ close from the series expressions in the form
τ rzJanssen σ rFourier − Bessel
= μ close ≈ μ exact
(3.115)
The Coulomb’s coefficient μ close is plotted in Fig. 3.12. Besides the oscillations of the truncation error and their increase towards the end due to the sampling error, the Janssen error, which is the difference between μ close and μ exact , can clearly be seen. From section 2.3, we recall Janssen’s stress analysis Eqns. (2.12) - (2.14), which is based on equilibrium of a differential slice as shown in Fig. 2.6, an appropriate K-value
σ rJanssen ( z ) = K ⋅ σ zJanssen ( z )
(3.116)
and the assumption of Coulomb’s friction at the wall
τ rJanssen ( z ) = μ ⋅ σ rWall ( z ) = μ ⋅ K ⋅ σ zJanssen ( z ) z
(3.117)
We interpret the uniform axial stress σ zJanssen ( z ) as an average stress of the FourierBessel axial stress field σ zFourier − Bessel (r , z ) across each horizontal section. It can be easily
66
verified by equilibrium considerations for any depth of the Janssen and the Fourier-Bessel solutions that this assumption must be correct. We can get the exact shear stress distribution from the Janssen model, provided that we know the exact K-value
K exact ( z ) =
σ rexact ( R, z ) ave σ zexact ( z )
(3.118)
where σ rexact ( R, z ) represents the exact stress distribution over the radial surface,
ave σ zexact ( z ) the exact averaged axial stress over each cross section and K exact ( z ) is expected to be a function of z rather than a constant. We actually do not know any of the terms involved in Eqn. (3.117), but we can make an approximating guess from the Fourier-Bessel solutions. From Fig. 3.12, we know that the constant Coulomb’s friction coefficient μ is underestimated with increasing z , which means that in Eqn. (3.117) the Fourier-Bessel radial stress σ rFourier − Bessel ( R, z ) at the wall is too high. If we replace in Eqn. (3.116) σ rJanssen ( z ) by σ rFourier − Bessel ( R, z ) , we find a new expression
K
closer
σ rFourier − Bessel ( R, z ) ( z) = σ zJanssen ( z )
(3.119)
that is much closer to Eqn. (3.118) than the original constant K-value. The new K-function K closer ( z ) can be plugged into the Janssen shear formulation Eqn. (3.111) to find a better shear distribution τ rcloser (z ) at the radial surface. This shear z expression can be used as a boundary condition in the Fourier-Bessel model to find the corresponding radial stress field σ rFourier − Bessel ( R, z ) and we can thus return to Eqn. (3.119). 67
This algorithm may be repeated until a specified small error of tolerance is reached. A simple iteration scheme is presented in Fig. 3.13 and results for the first two iterations are shown in Figs. 3.14 and Fig. 3.15. Whereas the original Janssen error is around 3%, it can be reduced to < 0.5 % with only two iterations. With an increasing number of iterations, the Fourier-Bessel solution can be improved towards the exact solution for constant μ over the radial boundary.
⎡i = 1 : number of iterations ⎢ ⎢ i F −B ⎢ σ i +1 series K = irJ ; ⎢ σz ⎢ ⎢ ⎢ Exponential regression ⎢ ⎢ i +1 K function = a 0 + a1 ⋅ exp(− z ) ⎢ ⎢ + a 2 ⋅ z exp(− z ) + a3 ⋅ z 2 exp(− z ) ; ⎢ ⎢ ⎢ Update with i +1K function : ⎢ ⎢ Janssen axial stress i +1σ zJ ⎢ Fourier − Bessel analysis i +1σ rF − B ⎢ ⎢⎣
Fig. 3.13 Simple algorithm for improving
68
μ close
K-functions 0.495 After 1st Iteration After 2nd Iteration 0.49
K-value
0.485
0.48
0.475 original K-value = 0.4706 0.47 0
1
2
3 4 5 6 7 Depth from loading end z [cm]
8
9
10
Fig. 3.14 K-functions after 1st and 2nd iterations
Improved Coulomb´s coefficient mu(close) Original Fourier-Bessel After 1st Iteration After 2nd Iteration
Coulomb´s coefficient mu [-]
mu exact = 0.36 0.36
0.355
0.35
0.345 0
1
Fig. 3.15
2
μ close
3 4 5 6 7 Depth from loading end z [cm]
after 1st and 2nd iteration
69
8
9
10
3.6 The “Zero-axial-stress” Problem in Continuum Modelling
Completing chapter 3, we discuss the elastic results, which are graphically illustrated in section 5.1, in terms of the “zero-axial-stress” problem. First, we prove that this problem cannot be solved by elasticity. We then discuss problems of continuous modelling in discrete materials and draw conclusions for the further treatment of the “zeroaxial-stress” problem in continuum theory.
3.6.1 The Simple Janssen Model Recall the original motivation for the previous analyses from Fig. 3.6 in section 3.4: Is it possible in a laterally constrained elastic cylinder that is pushed from one end, to find axial stresses σ z = 0 at some finite depth, if we impose Coulomb’s friction over the radial surface? We can now prove that the answer must be no.
Fig. 3.16 Janssen’s differential slice
Let’s first look at Janssen’s differential slice in Fig 3.16. Assume that the stress
σ z + dσ z at the bottom of the slice is zero. Equilibrium in z -direction then yields σ z ⋅ π R 2 = τ r z ⋅ 2π R ⋅ dz
(3.120) 70
With Eqn. (3.117) we find
σ z = μ ⋅ K ⋅ σ z ⋅ 2 R ⋅ dz 1 = μ ⋅ K ⋅ 2 R ⋅ dz
(3.121)
Since in Eqn. (3.121) a constant factor multiplied by an infinitesimal length can never be 1 , vertical equilibrium will not be satisfied in Fig. (3.16), unless a stress
σ z + dσ z ≠ 0 exists at the bottom of the slice. In a more practical interpretation, we could say that the actual axial stress can be reduced only by a very small portion of itself in each slice. We may take as a simple visualization a box of flour where a fixed portion of the actual content is taken out successively. This can be repeated indefinitely, but the box will never be cleared completely. We proved that refined solution fields for stresses and displacements that were determined by analytical methods of elasticity theory support this prediction of the simple Janssen model. We can attempt any cylinder length with the help of the Fourier-Bessel potential by increasing sampling interval [ 0; L ] and truncation term n , but eventually we will always find a residual stress at the bottom z = Lcyl .
3.6.2 Continuum Theory vs. Discrete Particle Behaviour The statement that axial stresses cannot be zero at any finite length of the cylinder is not restricted to elasticity alone, but holds for any continuous medium. Elasto-plasticity as introduced in chapter 5 is unable as well to find zero axial stresses. Hence, the reason for the impossibility of that problem formulation must be more general and is due to the concept of continuum theory.
71
When we describe granular materials as a continuum, we disregard the discrete particle structure and assume that the overall, large-scale behaviour of the granular body can be approximated by continuous stress and displacement fields. These solutions have to satisfy not only equilibrium, but also the restrictive requirement of compatibility, i. e. no discontinuities whatsoever are allowed. We may recall from the introduction, that there is a real-world counterpart to our “zero-axial-stress” problem. Experiments prove that a finite granular column supported only by wall friction can stay inside a tube and axial stress at the bottom must hence be zero. As we have seen in section 1.3, physicists are presently working on the phenomenon of arching and discrete force propagation. Various experiments indicate that grains can locally develop discrete mechanical structures (“arches”) that are responsible for non-uniform stress propagation (“force chains”). Hence the phenomenon of arching in granular matter is mainly characterized by non-homogeneous, discontinuous material behaviour.
3.6.3 The Remedy of a Constant Bottom Stress At this point, we have to differentiate clearly between the discontinuous local behaviour of discrete particles that continuum theory fails to account for and the overall stress and displacement characteristics of the whole granular body that continuum theory is able to predict in average very well, as we will see from experiments later on.
72
P
Fig. 3.17 Arching and force chains in a cylinder
Fig. 3.18 Bottom stress for continuum approach
The peculiar characteristics of the granular material allow the formation of an arching structure. It discretizes the continuous flow of stresses and at the bottom a complete transference of axial forces to the wall becomes possible (see Fig. 3.17). If this discrete, local behaviour is homogenized by a continuum approach - just like at any other location inside the granular body - , it loses its most crucial effect of making the grains stay inside the tube. We have to conclude that a continuum model cancelling out any discrete material characteristics cannot be an appropriate arching model.
73
Note that we still refer to an idealized weightless material to simplify our case. The concept of discrete force propagation is to be applied only to the forces resulting from the top surcharge P. The weight of the material can never be distributed to the walls completely, since it must be present at every particle, or it must be taken by another mechanism (see introduction). We suggest a simple interpretation of the discrete arching behaviour of the bottom part into continuum theory. Consider the arch as shown in Fig. 3.18. The granular body is cut at some depth over the arch and forces are idealized as constant axial stress. The effect of the arching part on the rest of the granular body is modelled by a constant supporting axial stress σ res which is provided at the bottom of the cut part in upward direction to account for the arching effect in continuum theory. To find the maximum loadbearing capacity of a certain length in terms of P, we will assume in chapter 5 that we know the maximum stress that the arching structure can withstand, and apply the simple model by substituting this σ res into the continuum analysis. By this stipulation we acknowledge that continuum theory is unable to model the behaviour in the cut bottom part. On the other hand, we claim that the strength of the cut arch is the determining factor for the maximum top surcharge σ 0 , that the pushed granular cylinder will finally be able to take. We will see that the latter assumption is supported by the experimental part of chapter 4. With σ res at the bottom it is possible to model stress and displacement fields by continuum theories, such as elasticity (FourierBessel model in chapter 3) and plasticity (Drucker-Prager cap model in chapter 5). The final results in terms of the mobilization force and stress and displacement fields in an elastic column are given in chapter 6.
74
Chapter 4: Experimental Parameter Identification and Mobilization Testing
4.1 Coulomb’s Friction Coefficient μ
We have to determine the Coulomb’s coefficient for friction between the plastic PMMA tube and the alumina granulate. The Jenike shear cell is the conventional method for measuring friction coefficients experimentally (Jenike 1970), which, however, does not allow easily for testing with a cylindrical specimen. Apart from finding flat PMMA material with the same properties, we also cannot estimate a possible influence of the geometric form of the contact surfaces on μ .
4.1.1 Basic Testing Concept To avoid these problems we are measuring the friction coefficient in situ, i. e. we modify the loading conditions in the original tube such that μ can eventually be determined. The resulting experimental set-up is shown in Fig 4.1. The basic concept can be
Fig. 4.1 Experimental set-up for friction coefficient test
summarized as follows:
75
(1)
A 5 cm layer of material is subjected to a constant axial force F by the loading machine. The tube (∅ 7.5 cm) is supported by frictional contact between plastic and grains only.
(2)
The corresponding radial bulging force N in the grains is measured at the wall by force sensors (∅ 0.95 cm).
(3)
Weights are stacked on the PMMA tube and successively increased until the granular cylinder collapses and moves down at a critical weight force T .
(4)
With the assumption that over a relatively small depth the radial stresses computed from N as well as the frictional shear stresses computed from T are evenly distributed over the contact surfaces, we find the coefficient μ from Eqn. (2.2) as
μ=
ave τ r z ave σ r
=
T π (0.95 / 2 cm) 2 ⋅ N 5cm ⋅ 2π ⋅ 3.75cm
(4.1)
4.1.2 Critical Discussion The outlined experimental method has theoretical and practical shortcomings whose influence we will review briefly in the following:
(1)
We could observe a considerable variation of radial wall stresses over the depth during the experiment. This is due to the weight force on the plastic tube that is completely transferred into the granular layer via friction, where it gradually adds to 76
the axial force from the loading machine. Accordingly, the bulging radial force must be smaller at the top of the granular layer than at the bottom. We account for this phenomenon by placing sensors at different depths and taking their reading while increasing the weights on the tube.
(2)
The sensor results are extremely scattered (see Fig. 4.2), which is mainly due to the discrete force propagation in the grains (see section 2.4). The response depends on whether a main force chain is hitting the sensing area or not. The scattering might also be amplified by the frictional shear of the weight force which does not add up to the axial stress homogeneously, but in a discontinuous manner. We eliminate the scattering by means of a simple polynomial regression in the form
y = a0 + a1 x + a 2 x 2 .
(3)
Experimental experience strongly indicates inaccurate sensor measurements for part of the results. However, any attempt to correlate this to the scattering could not be verified. In the scope of this work, we therefore assume this error source to be small and the sensor readings to be fairly accurate. For issues of technical details, calibration and error estimation we refer to the literature (Gopal 2001), (Anandakumar 2005), (FlexiForce User Manual 2006).
4.1.3 Results We carry out the experimental procedure at 4 different axial stress levels repeating the testing 5 times at each of them. Radial forces are measured by 4 sensors at different depths (see Fig. 4.1). We first compute the polynomial regression function N as shown
77
in Fig. 4.2 for radial stress readings from the sensors, whose different locations are not considered independently. We then compute the polynomial regression function T for the frictional shear stress that is taken from the maximum weight forces on the tube as shown in Fig. 4.3. We can relate T and N in the form of Eqn. (4.1) and find a function for μ dependent on the axial stress as shown in Fig. 4.4. We take as an average the final constant μ = 0.36 .
Radial stress for different axial loadings Radial stress at the wall from sensor readings [N/mm²]
0.07
0.06
Polynomial regression Experimental values
0.05
0.04
0.03
0.02
0.01
0 0.02
0.04
0.06 0.08 0.1 0.12 0.14 Axial stress from the loading machine [N/mm²]
0.16
Fig. 4.2 Sensor-measured axial stress induced by piston loadings of 0.25 – 0.5 – 1.0 – 1.5 N/mm²; 4 sensors x 5 repetitions are 20 data points per axial stress level
78
Frictional shear stress vs. axial stress 0.025 Polynomial regression Experimental values
Frictional shear stress [N/mm²]
0.02
0.015
0.01
0.005
0 0.02
0.04
0.06 0.08 0.1 0.12 Axial stress from loading machine [N/mm²]
0.14
0.16
Fig. 4.3 Frictional shearing stress at which the cylinder moved downwards; the recorded weight force on the cylinder is evenly distributed over the contact surface
Friction coefficient vs. axial stress
Friction coefficient mu [-]
0.4
0.36
0.3
mu from experimental regressions assumption of mu averaged 0.02
0.04
0.06 0.08 0.1 0.12 Axial stress from loading machine [N/mm²]
0.14
0.16
Fig. 4.4 Friction coefficient μ determined from the polynomial regression functions
79
4.2 Elastic Material Properties
Purely elastic material behaviour can be described by two independent material constants (see section 4.1). We determine the constrained modulus M from a uniaxial strain test and Young’s modulus E from a triaxial test for the alumina granulate.
4.2.1 Loading-Unloading Behaviour of Granular Materials Before presenting the experimental results, we may begin by recalling some theoretical background of elasticity (see section 2.1). The basic assumption of elastic material behaviour states that the material resumes its original shape when unloaded, i. e. all deformations are recoverable. However, this formulation is not in accordance with real granular materials, whose deformations are to a considerable extent not recoverable. These plastic deformations do not vanish upon unloading, but stay permanently. Accordingly, we can separate the strains of each material element into an elastic and a plastic part. This is called the elastic-plastic split and reads
ε ij = ε ije + ε ijp
(4.2)
Eqn. (4.2) indicates that the behaviour during loading is not governed by elasticity alone, but by a modified set of constitutive relations that can account for plastic deformation. However, if we consider only the unloading process, the material behaviour must be purely elastic, since the plastic strains ε
p
stay constant, when the load is removed. In
chapter 4, we will show an example of a modified set of constitutive relations in the form of the Drucker-Prager cap model, where the loading phase is described by an elastoplastic material law and the unloading phase is purely elastic.
80
Experiments in soils mechanics show that reliable elastic material properties can be derived from the unloading data of standardized testing methods (Chen, Baladi 1985), (Bardet 1997).
4.2.2 The Uniaxial Strain Test An alumina sample of small depth (5 cm) confined in the PMMA tube is loaded axially (see Fig. 4.5). During both loading and unloading we record axial stress and top displacement, from which axial strain can be computed. The loading-unloading curve for the 0.1 N/mm² stress level is shown in
Fig. 4.5 Experimental set-up for the uniaxial strain test
Fig. 4.9. Imagine a small material element in Cartesian coordinates in the center of the cylinder. For the unloading case, we write the elastic Hooke’s law in the form of Eqns. (3.6) – (3.8). Since the stiff tube prevents the sample from expanding radially, we can assume all horizontal strain elements to be zero and the horizontal stresses to be equal to σ radial .
ε x = ε y = 0; ε z ≠ 0;
(4.3)
σ x = σ y = σ radial
(4.4)
Substitution of Eqns. (4.3) and (4.4) into Eqns. (3.6) – (3.8) yield
σz = M εz
(4.5)
81
M=
E (1 −ν ) (1 + ν )(1 − 2ν )
(4.6)
Eqn. (4.6) is called the constrained modulus M and is represented by the slope of the unloading line in Fig. 4.9. We compute the slopes for 10 different samples and find as an average
M = 21.5 N / mm 2 .
4.2.3 Triaxial Test Set-up The triaxial test set-up consists of a cylindrical granular specimen which is encased in a triaxial cell filled with water as shown in Fig 4.7. The grains are protected by a rubber membrane and end caps against water infiltration. The lower cap is equipped with a porous disk to allow air to flow out (undrained test). Lateral stress in the form of water pressure can be adjusted by air pressure connected to the cell (consolidated test). Water pressure is preferred to
Fig. 4.6 Principal stresses on the surface of the triaxial specimen
direct air pressure because it is inert and has fewer disturbances. Axial stress is applied by both water pressure and piston force on the upper cap. During the experiment, the triaxial cell is sealed against the atmosphere to hold the inner air and water pressure.
82
The piston force and the displacement of the upper cap, which both can be read from the loading machine, and the air pressure are recorded
and transformed into
axial stress and strain and lateral stress, respectively. Note that the specimen is in a principal stress state, because there are no shearing stresses over the boundaries (see Fig. 4.6). Within the scope of this work,
Fig. 4.6 Experimental set-up for the triaxial test
we confine ourselves to consolidated, undrained triaxial compression tests that we will use here and in the next section in three different variations for determining the Young’s modulus, the Mohr-Coulomb failure envelope and the volumetric plastic strains. Further information on different types of triaxial tests can be found in (Bardet 1997).
4.2.4 Triaxial Test for Young’s Modulus We will now determine Young’s modulus from the unloading data of a triaxial test. First, we set a constant air pressure on the triaxial cell that exerts the hydrostatic stress
σ radial = σ axial = σ hyd on the specimen. Then, the axial stress is gradually increased up to a maximum stress level σ axial = σ hyd + σ piston and slowly removed again.
83
We recall that the unloading phase is considered purely elastic. In a small deformation setting, superposition is allowed, i. e. we can split the stresses in arbitrary portions. Consider Fig. 4.8: If we subtract the hydrostatic stress state, we are left with a uniaxial stress state induced by the piston force. Since the hydrostatic stress is not changed during unloading, the strain change can be entirely assigned to the uniaxial stress change and hence we find the following relation for unloading
σ Piston = σ axial − σ hyd = E ε axial
(4.7)
The loading-unloading data for the 0.1 N/mm² stress level is shown in Fig. 4.10. We conclude that the unloading slope represent Young’s modulus E . We record σ Piston and
ε axial during unloading of the specimen, we can compute E
from Eqn. (4.7). The average
Young’s modulus of all 3 tests carried out reads
E = 15 N / mm 2 . Using Eqn. (4.6) and the results E and M we can also compute the Poisson’s ratio ν for alumina granulate, which is ν = 0.32 .
84
Fig. 4.8 Separation of the Triaxial into a Hydrostatic and a Uniaxial Stress State
Loading - Unloading for Uniaxial Strain 0.11 0.1
Axial Stress [N/mm²]
0.08
0.06
Loading
0.04
M
0.02
0
Unloading
0
0.005
0.01 0.015 Axial Strain [-]
0.02
Fig. 4.9 Loading-unloading lines in the uniaxial strain test for max stress 0.1 N/mm²
Loading - Unloading for Triaxial Test 0.11 0.1
Piston stress [N/mm²]
0.08
Loading
0.06
0.04 E
0.02
0
Unloading
0
0.005
0.01
0.015 0.02 Axial strain [-]
0.025
0.03
0.035
Fig. 4.10 Loading-unloading lines in a triaxial test for max piston stress 0.1 N/mm²; note that the piston stress is axial stress minus radial stress
85
4.3 Material Properties for the Cap Model
For the implementation of the elasto-plastic Drucker-Prager cap model that we will discuss in detail in chapter 5, we need a set of material constants that account for the plastic behaviour of the alumina granulate during loading. Those are the granular shear strength (Mohr-Coulomb criterion) and the permanent compaction (evolution of plastic volumetric strains).
4.3.1 The Mohr-Coulomb Failure Envelope The Mohr-Coulomb criterion can be determined by a series of triaxial failure tests: We apply a constant hydrostatic compression to the specimen inside the triaxial cell (see Fig. 4.6) by setting a constant pressure. Then, the piston force is increased until a shearing deformation in the form of Fig. 2.2 can be observed. We repeat this procedure for 4 different hydrostatic pressure levels recording the peak piston force. From section 2.1, we recall the concepts of the Mohr-Coulomb failure criterion and the internal friction angle, which indicate the internal shear strength of a granular body. From the recorded data, the lateral stress ad the peak axial stress can be computed. Since they are both principal stresses, we can draw them in the form of Mohr’s circles (see Fig. 4.12). The Mohr-Coulomb criterion can be determined by fitting a tangential curve to the 4 circles. Since this curve is almost linear within the stress range of interest (0 – 0.2 N/mm²), we choose a straight line, which corresponds to the internal yield locus from section 2.1 and can be described by the internal friction angle φ = 37° and the cohesion c = 0 (see Fig. 2.4 and Fig. 4.12). The curved failure envelope at higher stress ranges is a phenomenon, which is well-known from various soils (Bardet 1997).
86
4.3.2 Porosity Change and Plastic Volumetric Strains We introduce the deviatoric split of strain and stress tensors in the form
e p ε ij = ε vol + eij = (ε vol + ε vol ) δ ij + (eije + eijp )
(4.8)
σ ij = σ hyd δ ij + sij
(4.9)
with volumetric strains ε vol and the hydrostatic stresses σ hyd ; and the deviatoric strains
eij and deviatoric stresses s ij that are defined as follows:
ε vol = 1 / 3 ε kk
and
eij = ε ij − ε vol δ ij and
σ hyd = 1 / 3 σ kk
(4.10)
sij = σ ij − σ hyd δ ij
(4.11)
Consider a granular specimen subjected to a hydrostatic load in the triaxial test. Upon removal of the load, an unrecoverable plastic deformation can be observed (permanent compaction), which can be expressed as an increase in the specific density (decrease of void ratio) of the material. Since we limited the stress range to levels where the alumina granulate cannot yield or be otherwise destroyed (< 0.2 N/mm²) (see section 1.2), the density increase cannot be due to plastic deformation in the alumina particles themselves, but must be due to a change in the porosity n or void ratio e of the particle aggregate. We can imagine that the specific behaviour of the specimen under pure shear deformation is not characterized by porosity change, but by tangential sliding as shown in Fig. 1.2. Therefore, we assume a correlation between hydrostatic stress and change of porosity such that the corresponding volumetric plastic strains represent the amount of pore volume loss. The corresponding volumetric elastic strains can be interpreted as the elastic recoverable deformation of the individual particles. 87
4.3.3 The Triaxial Flowmeter Test To incorporate this plastic behaviour into the elastoplastic model in chapter 5, the relation between hydrostatic stresses and plastic volumetric strains representing the porosity change has to be determined experimentally for alumina granulate. Consider again the granular specimen under hydrostatic load inside the triaxial cell (see Fig. 4.6). Since the granular medium inside the rubber membrane is connected to the outer atmosphere via the porous disk and the open drainage at the bottom, air is flowing out as the specimen deforms under load. In fact, the flow rate of air can be directly correlated with the change of porosity dn as
p = dn = dε vol
1 ∂ Vair ⋅ V0 ∂σ hyd
(4.12)
Fig. 4.11 Direct Reading Flowmeter: 1 – 250 mL/min Flow Range +/- 5% Reading Accuracy +/- 1% Repeatability
where V 0 represents the initial volume of the specimen. p For determining dε vol , the triaxial cell is extended by a flowmeter that is connected to
the drainage at the bottom. This device measures the amount of air that is flowing out with respect to time. Therefore we need to modify Eqn. (4.12) by chain rule
p dε vol =
1 ∂ Vair ∂ t ⋅ V0 ∂ t ∂σ hyd
(4.13)
88
where ∂ t / ∂σ hyd is the inverse of the rate of hydrostatic pressure on the cell. To evaluate Eqn. (4.13) practically, we increase the pressure by a constant rate and videotape the flowmeter. This device shown in Fig. 4.11 indicates the current flow rate by the position of a ball moving vertically on a direct reading scale. Afterwards, we take discrete flow values for very small time steps from the videotaping, multiply these flow p (Simpson rule). The hydrovalues by the time step size and sum the results up to ε vol
static stress and the corresponding plastic volumetric strains, which are used as the hardening parameter in chapter 5, are shown in Fig. 4.13. For the implementation in the cap model, an exponential regression of these discrete results is required in the form p ε vol = W [1 − exp(− Dσ hyd )]
(4.14)
The parameters W = 0.28 and D = 7.25 are found by trial and error. The best fit curve for the alumina granulate is also shown in Fig. 4.13. For details on flowmeters, we refer to the technical reference website of the manufacturer (Omega Engineering, Inc. 2006).
89
Mohr-Coulomb Failure Envelope 0.5 0.4
Shear Stress [N/mm²]
0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 0
0.2
0.4
0.6 0.8 Principal Stress [N/mm²]
1
1.2
Fig. 4.12 Mohr-Coulomb failure criterion for a the cohesionless alumina granulate with an internal friction angle of 37°
Hardening Parameter - Plastic Volumetric Strains 0.45 Hydrostatic Test Data Exponential Regression
0.4
Hydrostatic Stress [N/mm²]
0.35 0.3 0.25 0.2 0.15 0.1 0.05 0
0
0.05
0.1 0.15 0.2 Plastic Volumetric Strain [-]
0.25
0.3
Fig. 4.13 Hardening parameter; note that the regression was chosen in order to best fit the behaviour in the valid stress range [0 – 0.2 N/mm²]
90
4.4 Mobilization Test
We suggest a simple testing method to acquire experimental data for the mobilization force of a finite length bed of alumina granulate confined in the PMMA tube. The experimental outcome will be compared in chapter 6 with results from the analytical elastic and computational elasto-plastic model.
4.4.1 Experimental Set-up We recall the original situation as described in section 1.1: A finite column of alumina granulate gets stuck in the rigid plastic tube, where it is supported only by the frictional shear stress over the wall contact surface. We want to find out how much force we have to apply at the top to mobilize the granular column, i. e. it is sliding down the wall. However, we demonstrated in section 3.6 that continuum models have difficulties in describing the granular behaviour at the bottom of the column, where stresses are zero in reality. We presented recent theories of discrete force propagation that give a flavour of the complex mechanisms leading to zero stresses at some finite depths of the material. Finally, we suggested a simple remedy that cuts away the lower part of the column and replaces it by a distributed force at the upper part (see Fig 3.18). To obtain experimental results that reflect this specific remedial situation, we modify the original situation as shown in Fig. 4.14: The granular column is not pushed from the top, but from the bottom (“upside down”) by the piston force F , while the top of the granular cylinder is held by the upper piston. The important advantage of this reverse testing is the opportunity of putting a weight force on top of the granular column that can
91
be controlled and adjusted easily. This weight force models the distributed load in the remedial situation in Fig. (3.18) accurately. Note that the upside down testing is possible, since the weight of the alumina granulate is very small and its influence can thus be neglected. The experimental procedure is simple. The piston force F is increased slowly with a corresponding stress rate of about 0.0005 N/mm² per second so that one complete
Fig. 4.14 Experimental set-up for the mobilization test
test may take up to 30 minutes. Thus, all dynamic effects can be minimized and F is assumed static. The bottom displacement is given by the loading machine and the top displacement is measured by a dial gage. The peak force before the granular column begins to slide freely is taken as the mobilization force. A series of similar tests, where glass spheres in a tube were pushed in reverse, have recently been carried out by (Kolb et al. 1999) and (Vanel, Clément 1999).
92
4.4.2 Results The experimental procedure was carried out for column lengths of 7.5 cm, 15 cm and 22.5 cm with different weight forces on top. Note that in case of no top weight there is still a small distributed load on top, which is due to the material weight of the top layer. The corresponding mobilization forces for each case are listed in Fig. 4.15. From analytical and computational models presented in section 6.3, we know that the mobilization stress should increase linearly with increasing weight force on top. Although the amount of scattering in this experiment is much smaller than in section 4.1, since discrete effects are more likely to average out over the bigger dimension, the linear behaviour can still not be clearly confirmed by the experimental results in Fig. 4.15. The corresponding top and bottom displacements are shown in Fig. 4.16. They show the typical behaviour that was already pointed out in (Gopal 2001) and (Anandakumar 2005). A large amount of scattering is due to the sensitivity of the material to the moisture content of the air (humidity) that influences the mechanical properties considerably. See section 5.3 for a more detailed discussion on that.
93
Fig. 4.15 Mobilizing stresses for three different bed lengths; note the considerable scattering; the anticipated linear tendency cannot be definitely confirmed
Top and Bottom Displacement Behaviour 0.25
Top Displacements
Bottom Displacements
Applied Stress [N/mm²]
0.2
Length 22.5cm 0.15
0.1 Length 15.0cm 0.05 Length 7.5cm 0
0
0.5
1
1.5 2 2.5 Displacement [mm]
3
3.5
4
Fig. 4.16 Load-displacement for bottom (Loading machine platen) and top (Dial gage) of the granular column in mobilizing tests with 0.0045 N/mm² Top Weight
94
Chapter 5: The Elasto-plastic Drucker-Prager Cap Model
5.1 Elasto-plastic Modelling in Geomechanics
Soil or geomechanics is a continuum approach to describe stresses and displacements in soils. This branch of solid mechanics uses the same framework of equilibrium, kinematic and compatibility equations (3.15) – (3.20), but tries to account for the complex behaviour of granular soils by refining and specialising the constitutive relations. The state of the art in geomechanics in terms of mathematical formulation, numerical solution procedures and experimental verification techniques is highly mature and has been successfully applied for decades in practical geotechnical engineering.
5.1.1 The Incremental Flow Theory of Plasticity Under working load conditions, soil may be idealized as elastic, i. e. the material returns to its original shape after load removal. In other words, there exists a one-to-one coordination between stress and strain increments. However, these formulations fail to identify plastic deformations when unloading occurs. The flow theory of plasticity represents an extension of elastic stress-strain relations into the plastic range. Its basic premise is the assumption that an elasto-plastic material can undergo elastic (recoverable) as well as plastic (permanent) deformations at each loading increment. Therefore, we define - in analogy to the elasto-plastic split Eqn. (5.2) - the total strain increment tensor as the sum of the elastic and plastic strain increment tensors as
95
dε ij = dε ije + dε ijp
(5.1)
The flow theory is furthermore based on three fundamental assumptions: (1) the existence of a yield surface (2) the formulation of an appropriate flow rule (3) the evolution of subsequent loading surfaces (hardening rule); note that we confine ourselves in sections 5.1 and 5.2 to the case of perfect plasticity, i. e. no hardening occurs.
5.1.2 Yield Criterion, Flow Rule and Plastic Potential The yield surface defines the point at which elastic relations cease to be valid, i. e. the onset of plastic behaviour. It can be described mathematically as a function of the stress tensor
f (σ ij ) = f c
(5.2)
where f c is a constant for a perfectly plastic material. Eqn. (5.2) is represented by a hypersurface in the 6-dimensional stress space, but can also be graphically illustrated in the 3-dimensional principal stress space. The best known yield criteria for soils out of a great variety are shown in Fig. 5.1. The general shape of a yield surface can be best described by its cross-sectional shapes in the deviatoric plane (left-hand side) and its meridians in the hydrostatic or meridional plane (right-hand side). The shear strength of soils increases with the hydrostatic stress, which is well reflected in coned shapes. Note that the sign convention in soil mechanics denotes compressive stresses as positive. 96
Elastic behaviour occurs, if after an increment of stress dσ ij , the new state of stress is in the elastic domain, i. e. lies below the yield surface, that is
f < fc
(5.3)
or the state of stress moves off the yield surface
f = fc
and
df =
∂f dσ ij < 0 ∂σ ij
(unloading)
(5.4)
Plastic deformation occurs, if after an increment of stress, the state of stress is in the plastic domain, i. e. lies on the yield surface, that is
f = fc
and
df =
∂f dσ ij = 0 ∂σ ij
(loading)
(5.5)
The state of stress can never be over the yield surface. This set of equations (5.3) – (5.5) is called the loading-unloading or consistency conditions and determines, if the material behaves elastically or plastically. The flow rule constitutes the relationship between the next increment of plastic strain dε ijp and the present state of stress σ ij for an element on the yield surface. It is defined by the plastic potential function g in the form
dε ijp = dλ
∂g ∂σ ij
(5.6)
where dλ is a positive scalar of proportionality depending on stress sate + load history.
97
Fig. 5.1 Common yield criteria in soil plasticity; note that we will focus on the Drucker-Prager criterion in the following
98
If the yield and potential surfaces coincide with each other ( f = g ), we call Eqn. (5.6) an associated flow rule. In this case, it gives the direction of the incremental plastic strain tensor normal to the plastic potential g . Drucker (1951) showed that this normality condition together with a convex shape of the plastic potential guarantees the uniqueness of the solution for a boundary value problem in perfect plasticity. For further details on flow theory in geomechanics (hardening rules, non-associated flow, dilatancy modelling), we refer to (Chen, Mizuno 1990), (Potts, Zdravkovich 1999) and (Dunne, Petrinic 2005).
5.1.3 Generalized Stress-Strain Relations We will now give a short overview on the derivation of the general stress-strain relations for an ideal plastic material. The main target is to find an expression for the yet unknown proportionality factor dλ . We assume a displacement-driven setting; i. e. dε ij is specified at the beginning of each incremental step.
dε ij = dε ije + dε ijp =
d I1 ∂f 1 δ ij + dsij + dλ 9K 2G ∂σ ij
(5.7)
We know from Eqns. (5.2) – (5.5) that
dλ = 0 , if f < f c or f = f c but df < 0 dλ > 0 , if f = f c and df =
(5.8)
∂f dσ ij = 0 ∂σ ij
(5.9)
We can now determine the form of dλ by combining Eqns. (5.7) and (5.9). Solving Eqn. (5.7) for dsij , we find the stress increment tensor
99
dσ ij = dsij + 1 / 3 d I1δ ij = 2G dε ij − 2G dλ
∂ f ⎛ 1 2G ⎞ +⎜ − ⎟dI1 δ ij ∂σ ij ⎝ 3 9 K ⎠
(5.10)
Substituting Eqn. (5.10) into (5.9) yields
2G
∂f ∂ f ∂ f ⎛ 1 2G ⎞ ∂f dε ij − 2G dλ +⎜ − δ ij = 0 ⎟ dI 1 ∂σ ij ∂σ ij ∂σ ij ⎝ 3 9 K ⎠ ∂σ ij
(5.11)
We can now eliminate dI 1 by using Eqn. (5.7) with i = j
⎛ ⎞ ∂f ⎟ dI 1 = 3K ⎜ dε kk − dλ δ ij ⎜ ⎟ σ ∂ ij ⎝ ⎠
(5.12)
Substitution of Eqn. (5.12) into Eqn. (5.11) gives the proportionality factor
dλ =
∂f ∂f 3K − 2G δ ij dε ij + dε kk ∂σ ij 6G ∂σ ij ⎞ ∂f ∂f 3K − 2G ⎛ ∂ f ⎜⎜ δ mn ⎟⎟ 2 + 6G ⎝ ∂σ mn ∂σ mn ∂σ mn ⎠
(5.13)
For a displacement-driven setting, the factor dλ can now be determined uniquely. Finally, we can determine the stress increment dσ ij by substituting dλ into Eqn. (5.10).
100
5.2 The Drucker-Prager Model
We focus on the Drucker-Prager failure criterion and its general formulation for the implementation as the constitutive relation into a finite element code.
5.2.1 The Drucker-Prager Yield Criterion The Drucker-Prager yield criterion attempts to approximate the well-known Coulomb criterion by a smooth function without corners and is expressed as a simple function of the first invariant of the stress tensor I 1 and the second invariant of the deviatoric stress tensor J 2 , together with two material constants α and k .
I 1 = σ ii = σ 11 + σ 22 + σ 33
[
]
J 2 = 1 / 6 (σ 11 − σ 22 ) 2 + (σ 22 − σ 33 ) 2 + (σ 33 − σ 11 ) 2 + σ 12 + σ 23 + σ 13 2
2
f ( I1 , J 2 ) = α I1 + J 2 = k
2
(5.14)
It may be verified by inspection that Eqn. (5.14) depicts a circular cone with symmetry about the hydrostatic axis as shown in Fig. 5.1. Since the Drucker-Prager criterion has been established as an approximation of the Coulomb criterion, it is common to determine the material constants α , k by matching two particular points with those of the Coulomb criterion. We thus get α , k in terms of the familiar Coulomb constants φ , c . In the principal stress space, the Drucker-Prager criterion can be matched with the apex of the Coulomb criterion as an outer cone as shown in Fig. 5.2.
101
Fig. 5.2 Matching of yield criteria in the deviatoric plane
Since a line element connecting the apex O and point A contains the same line for both criteria, the relation between the two sets of material constants can be found by comparing the invariant form of the Coulomb criterion for θ = 0 Eqn. (2.5) to Eqn. (5.15) (Zienkiewicz 1978), (Chen, Mizuno 1990).
2 sin φ 3 (3 − sin φ )
α=
k=
I1 + J 2 =
6 c cos φ
(5.15)
3 (3 − sin φ )
2 sin φ
(5.16)
3 (3 − sin φ ) 6 c cos φ
(5.17)
3 (3 − sin φ )
102
From the Mohr-Coulomb parameters φ = 37° and c = 0 that we found from experiments in section 3.3, we find α = 0.29 (angle of friction tan −1 (α ) = 16.16° ) and the cohesion k = 0 for the Drucker-Prager model.
5.2.2 Flow Theory Formulation for Drucker-Prager If we substitute the Drucker-Prager yield criterion Eqn. (5.14) into the generalized stress- strain relations given in sections 5.1.3, we find the constitutive relations of the Drucker-Prager material model that can be integrated in a FE code:
a) Stress-strain relations
For f = α I 1 +
dε ij =
J 2 − k < 0 , or f = α I 1 + J 2 − k = 0 and df =
⎫ ⎪ ⎬ elastic domain ⎪ ⎭
d I1 1 δ ij + ds ij 9K 2G
dσ ij = K dε kk δ ij + 2G deij
For f = α I 1 +
dε ij =
J 2 − k = 0 and df =
∂f dσ ij < 0 ∂σ ij (5.18)
∂f dσ ij = 0 ∂σ ij
⎡ sij d I1 1 δ ij + dsij + dλ ⎢− αδ ij + 9K 2G 2 J2 ⎢⎣
⎤ ⎥ ⎥⎦
⎡ G sij ⎤ dσ ij = K dε kk δ ij + 2G deij − dλ ⎢− 3Kαδ ij + ⎥ J 2 ⎦⎥ ⎣⎢
⎫ ⎪ ⎪⎪ ⎬ elasto-plastic domain (5.19) ⎪ ⎪ ⎪⎭
where dλ is given by
− 3Kα dε ij + dλ =
G J2
s mn demn (5.20)
9 Kα 2 + G
103
b) Stiffness formulation
For the direct use in a FE code, Eqn. (5.19) has to be reformulated in the form
ep dσ ij = C ijkl dε ij
(5.21)
whose matrix form in Voigt notation for the axisymmetric case reads
⎧ dσ r ⎫ ⎪ dσ ⎪ ⎪ θ⎪ ep ⎬ = C ⎨ ⎪ dσ z ⎪ ⎪⎩dσ r z ⎪⎭ 4 x1
[ ]
4 x3
⎧ dσ 11 ⎫ ⎧ dε r ⎫ ⎪dσ ⎪ ⎪ ⎪ ⎪ 22 ⎪ ep or ⎨ ⎨ dε z ⎬ ⎬ = C ⎪2 dγ ⎪ ⎪dσ 33 ⎪ r z ⎭ 3 x1 ⎩ ⎪⎩ dσ 13 ⎪⎭ 4 x1
[ ]
4 x3
⎧ dε 11 ⎫ ⎪ ⎪ ⎨ dε 33 ⎬ ⎪2 dγ ⎪ 13 ⎭ 3 x1 ⎩
(5.22)
By reducing Eqn. (5.19) axisymmetrically and substituting into Eqn. (5.21) we get the elasto-plastic tangent moduli after some algebraic manipulations (Chen, Mizuno 1990)
[C ] ep
⎡ K + 4 / 3G ⎢ K − 2 / 3G = ⎢ ⎢ K − 2 / 3G ⎢ 0 ⎣
⎡ H 11 2 0⎤ ⎢ K − 2 / 3G 0 ⎥⎥ 1 ⎢ H 22 H 11 − K + 4 / 3G 0 ⎥ H ⎢ H 33 H 11 ⎢ ⎥ 0 G⎦ ⎢⎣ H 13 H 11 K − 2 / 3G
H 11 H 33 H 22 H 33 H 33
2
H 13 H 33
H 11 H 13 ⎤ ⎥ H 22 H 13 ⎥ (5.23) H 33 H 13 ⎥ 2 ⎥ H 13 ⎥⎦
with the abbreviations H = 9 Kα 2 + G , H11 = 3Kα + G s11 / J 2 ,
H 22 = 3Kα + G s 22 / J 2 , H 33 = 3Kα + G s33 / J 2 and H 13 = G s13 / J 2 .
5.2.3 The Radial Return Algorithm for FE implementation At this point, we want to confine ourselves to a brief algorithm scheme for the radial return algorithm for an associated, perfectly plastic Drucker-Prager material. This special algorithm, which is implicit and unconditionally stable, has been developed to incorpo-
104
rate elasto-plasticity in finite element computations. Detailed information on the radial return algorithm can be found in (Simo, Hughes 1987), (Chen, Mizuno 1990), (Lubarda 2002), (Dunne, Petrinic 2005).
1) Input from FE code: σ ijn , dε ijn +1 2) Deviatoric split:
dε kkn +1 = trace(dε ijn +1 ) ;
deijn +1 = dε ijn +1 − 1 / 3 dε kkn +1
3) Compute elastic trial stresses: ∗
I 1e = I 1n + 3K dε kkn +1 ;
∗ e ij
s = sijn + 2G deijn +1 ;
4) Check trial yield criterion function: If
∗
f ( ∗I 1e , ∗ s ije ) ≤ 0 Then Elastic step: I 1n +1 = ∗I 1e ; s ijn +1 = ∗ s ije ; H ij = 0 ;
If
∗
f ( ∗I 1e , ∗ s ije ) > 0 Then Scale back stresses such that
∗
f ( ∗I 1e , ∗ s ije ) ≅ 0 (radial return method)
∗ e ⎡ J 2 − g ( ∗I 1e ) ⎤ ∂ g ⎥ ; = −3 ⎢ 2 ⎢⎣ 9 K (dg / dI 1 ) + G ⎥⎦ ∂ I 1 g ( I 1n +1 ) ∗ e I 1n +1 = ∗I 1e − 3K dε kkp n +1 ; sijn +1 = sij ; H ij g ( ∗I 1e )
Elasto-plastic step: dε kkp
n +1
5) Compute tangent moduli C ijkl [see Eqn. (4.20)] 6) Output to FE code: σ ijn +1 , Cijkl
105
5.3 The End Cap Extension
The Drucker-Prager model can accurately represent the shear strength of soil, but fails to account for its behaviour in case of similar principal stresses σ 1 ≈ σ 2 ≈ σ 3 . This can be easily demonstrated in the principal stress space (see Fig 5.1): If stresses stay close to the space diagonal, which represents the hydrostatic stress state, they will never encounter the Drucker-Pager yield cone and are thus modelled elastically even for very high compression. However, hydrostatic experiments show that soils exhibit considerable unrecoverable (plastic) deformation in such stress states, which is mainly due to the change of void ratio (porosity).
5.3.1 The Concept of the End Cap Drucker himself introduced the basic two innovations for the development of a model that can represent plastic deformations for stresses close to the hydrostatic diagonal (Drucker et al. 1957): (1)
A spherical end cap is fitted to the Drucker-Prager cone that serves as an additional yield surface. The cap is treated as an elasto-plastic hardening material, i. e. in a yielded state both the cone and the end cap expand along the space diagonal (see Fig. 5.3)
(2)
The specific porosity of each corresponding hydrostatic stress state is used as the hardening parameter to determine the successive loading surfaces (see also section 4.3 for details).
106
From these ideas, two specific types of isotropic strain-hardening models were developed that are both widely used in geotechnical engineering applications: On the one hand, there are the Cambridge models (Cam-clay models) that are mostly used for cohesive soils and whose principles are summarized in (Schofield, Wroth 1968). Further developments and FE implementations can be found in (Potts, Zdravkovich 1999). On the other hand, there are the generalized cap models that were originally developed for non-cohesive soils such as loose sands (DiMaggio, Sandler 1971). Since we concentrate on non-cohesive materials that show similar properties to sand, we will focus on the Drucker-Prager Cap Model, which belongs to the latter category. It bases on Drucker’s original idea to cap a cone-shaped yield criterion such as the Drucker-Prager model by a strain-hardening elliptical end sphere. In the following, we very briefly outline the mathematical formulation.
5.3.2 Cap Surface and Plastic Hardening Computations We describe the end cap in the meridional plane as a successive series of ellipses that are defined by a constant ratio R of the minor and major axes and the current intersection points L(ε kkp ) and X (ε kkp ) with the Drucker-Prager cone and the I 1 -axis, respectively (see Fig. 5.4). X (ε kkp ) can be expressed as
X (ε kkp ) = −
1 ⎛ ε kkp ln⎜1 − D ⎜⎝ W
⎞ ⎟⎟ ⎠
(5.24)
where D and W are material constants and ε kkp is the hardening parameter (see section 4.3 for introduction and experimental determination).
107
Fig. 5.3 Drucker-Prager type of strain hardening cap model
Fig. 5.4 Yield surfaces for the Drucker-Prager model in the meridional plane
108
The elliptic 3-D sphere function can hence be written as
f cap ( I 1 , J 2 , ε kkp ) =
[
1 ( X (ε kkp ) − L(ε kkp )) 2 − ( I 1 − L(ε kkp )) 2 R
]
1 2
− J2 = 0
(5.25)
According to Drucker’s stability postulate for hardening materials (Drucker et al. 1957), we assume an associated plastic potential for the cap surface in the form
f cap = g cap . The incremental hardening computation in the case of cap yielding can then be briefly outlined as follows:
(1)
We find from the incremental strains dε ij and the corresponding elastic trial stresses ∗ I 1e and ∗ sije that the cap yield condition is violated in the form
(
f cap ∗ I 1e , (2)
∗
J 2e , ε kkp
n
)> 0
and
∗
I 1e ≥ L(ε kkp n )
(5.26)
Analogous to the Drucker-Prager yield surface, the target is now to find the new stress state I 1n+1 , sijn+1 , the corresponding elasto-plastic tangent moduli Cijkl and additionally the new hardening parameter ε kkp A trial hardening parameter ε kkp
t
n +1
.
can be determined by an iterative scheme
given in (Chen, Baladi 1985). With ε kkp t , we find from the approximated yield condition Eqn. (5.27) and flow rule Eqn. (5.28) the relation (5.29)
f cap ( I 1t , J 2t , ε kkp t ) ≅ 0 dε
p ij
n +1
∂ g cap = dλ = dλ ∂σ ij
(5.27)
⎛ sijt ⎞ ⎜ − a δ ij ⎟ ⎜2 Jt ⎟ 2 ⎝ ⎠
109
(5.28)
[
∂ ( X (ε kkp t ) − L(ε kkp t )) 2 − ( I 1 − L(ε kkp t )) 2 with a = ∂ I 1t de
(3)
p ij
n +1
=s
dε kkp
n +1 ij
]
1 2
n +1
(5.29)
6a J 2t
As proven in (Chen, Baladi 1985), the final stress deviators and the elastic trial stress are related by
s ijn +1 + 2G deijp
n +1
= ∗ s ije
(5.30)
Substitution of Eqn. (5.29) into (5.30) and multiplication of each side by itself gives
J 2n +1 + G
(4)
dε kkp n +1 = 3a
∗
J 2e
(5.31)
This relation can be evaluated for dε kkp
n +1
and
J 2n +1 and hence I 1n+1 and sijn+1 by
means of a modified regula falsi method and the elasto-plastic tangent moduli can be calculated.
Details on this procedure are given in (Chen, Baladi 1985) and (Chen, Mizuno 1990).
5.3.3 A Modular Code Structure for FE Implementation According to the complex combination of yield surfaces, several elastic trial stress paths are possible and the final FE material code has to account for any of the following scenarios:
110
(1)
Elastic step: The trial stresses do not violate any of the yield surfaces and the element behaviour is truly elastic.
(2)
Elasto-plastic step: The trial stresses violate the Drucker-Prager yield criterion and the radial return algorithm from section 5.2.3 applies.
(3)
Hardening step: The cap surface is violated and the hardening computation as briefly outlined in section 5.3.2 applies.
(4)
Corner step: The final stress state is on the corner at which the yield surface intersects the cap. This will always occur, when the hardening surface is allowed to move back (dilatancy reduction). An appropriate treatment of this situation is given in (Chen, Mizuno 1990).
A basic algorithm scheme for the FE implementation of the Drucker-Prager cap model indicates the modular structure, which allows the splitting of the code in several independent routines that are later run within the main material routine.
1) Input from FE code: σ ijn , dε ijn +1 2) Deviatoric split:
dε kkn +1 = trace(dε ijn +1 ) ;
deijn +1 = dε ijn +1 − 1 / 3 dε kkn +1
3) Compute elastic trial stresses: ∗
I 1e = I 1n + 3K dε kkn +1 ;
∗ e ij
s = sijn + 2G deijn +1 ;
111
4) Check yield surfaces and cap position:
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩
[
]
⎡If ∗f yield ( ∗I 1e , ∗ sije ) ≤ 0 and ∗f cap ∗ I 1e , ∗ sije , X n (ε kkp ) ≤ 0 Then ⎢ Elastic step : I 1n +1 = ∗I 1e ; sijn +1 = ∗ sije ⎢⎣ ⎡If ∗f yield ( ∗I 1e , ∗ sije ) > 0 and ∗I 1e ≥ Ln (ε kkp ) Then ⎢ ∗ e ⎡ ⎢ J 2 − g ( ∗I 1e ) ⎤ ∂ g ⎥ Elasto- plastic step : dε kkp n +1 = −3 ⎢ ⎢ 2 ⎢⎣ 9 K (dg / dI 1 ) + G ⎥⎦ ∂ I 1 ⎢ ⎢ g ( I 1n +1 ) ∗ e ⎢ n +1 ∗ e p n +1 n +1 3 ; I I K d s sij ε = − = kk ij 1 1 ⎢ g ( ∗I 1e ) ⎢ ⎢ ⎢ ⎡If I 1n +1 < L(ε kkp n +1 ) Then ⎢ ⎢ ⎢ Corner step : see (Chen, Mizuno 1990) ⎣ ⎢ ⎣
[
]
⎡If ∗f cap ∗ I 1e , ∗ sije , X n (ε kkp ) > 0 and ∗I 1e < Ln (ε kkp ) Then ⎢ dε kkp n +1 ⎢ n +1 Hardening step : R egula falsi with J G + = 2 ⎢⎣ 3a
∗
J 2e
5) Compute tangent moduli C ijkl [see section 3.2.3 and (Chen, Mizuno 1990)] 6) Output to FE code: σ ijn +1 , C ijkl
This algorithm is written in a Fortran language for the integration as a material model in a FE code and can be found in Appendix B. It is a combination of code samples given by (Sandler et al. 1976), (Chen, Mizuno 1990) and (Dunne, Petrinic 2005).
112
5.4 The Drucker-Prager Cap Model in Simple Tests
The Drucker-Prager cap model is integrated as a user subroutine into the finite element software package ABAQUS. Before using the Drucker-Prager cap model in larger finite element structures, we study its behaviour for two basic loading-unloading scenarios in a single axisymmetric, linear element.
5.4.1 General Incremental Finite Element Procedure In general, the standard finite element procedure for static analysis transforms the non-linear equations of equilibrium for an elasto-plastic problem to the incremental form of linearized finite element equations (Bathe 1996).
[K ]n {dU }n+1 = {dR}n+1
(5.32)
From the total stiffness matrix [K ] and the (n+1)-th incremental load vector {dR}
n +1
nodal displacement increments {dU }
n +1
strain increments {dε }
n +1
can be computed. From the known {dU }
n +1
,
, the
can be calculated according to the element formulation at each
integration point. These basic finite element operations can be solved by any commercial FE package. Here, we will use the program ABAQUS (2005). Keep in mind that the coding in Appendix B needs input quantities in terms of the Cauchy stress tensor σ ij and the linearized strain tensor ε ij and is therefore restricted to small deformations. We find the corresponding stress increments {dσ }
n +1
ments {dε }
n +1
by putting the strain incre-
into the constitutive relations. In our case, these stress-strain relations are
represented by the Drucker-Prager cap model, whose output quantities are the updated material stress response
113
{σ }n+1 = {σ }n + {dσ }n+1
(5.33)
and the updated elasto-plastic tangent moduli that can be assembled to the total stiffness matrix [K ]
n +1
for the (n+2)-th increment.
In section 5.3.3 we briefly outlined the algorithm scheme of the Drucker-Prager cap model and implemented it into a full Fortran code, which is given in Appendix B. This modular code structure comprises the Drucker-Prager subroutine, which is based on the radial return algorithm given by (Dunne, Petrinic 2005), and the cap subroutine, which is given by (Chen, Baladi 1985) and (Chen, Mizuno 1990). It furthermore contains the commands that are necessary for ABAQUS to recognize this code as the user subroutine UMAT.
5.4.2 A Simple One Element FE Model In the following, we demonstrate the characteristic behaviour of the Drucker-Prager cap model for a single axisymmetric linear element as shown in Figs. 5.5 and 5.6 under hydrostatic stress and triaxial shearing. The software package ABAQUS - together with the material code in Appendix B that is integrated as the user subroutine UMAT - are used for all computations. Note that ABAQUS requires input and output quantities in axisymmetric Voigt notation Eqn. (5.22). We recall the material parameters for alumina granulate from sections 4.2, 4.3 and 5.2 and integrate them in UMAT.
114
Fig. 5.5 Axisymmetric element and its position in the cylinder
Fig. 5.6 Single element model with boundary conditions as used in the following
115
-
Elastic modulus E = 15 N / mm 2 and Poisson’s ratio ν = 0.32
-
Angle of friction α = 16.16° and cohesion k = 0
-
W = 0.28 and D = 7.25 for the hardening parameter curve
-
The initial cap position in the meridional plane is 0 and the cap eccentricity
R = 1.4 is assumed from the original cap model for sand (DiMaggio, Sandler 1971). We discussed the similarity between alumina and sand in section 2.4.
5.4.3 Isotropic Consolidation During isotropic consolidation (hydrostatic pressure only), the volumetric strain ε kk consisting of the elastic part ε kke and plastic part ε kkp is generated from the hardening cap portion of the model. The following relations hold:
Hydrostatic Stresses:
σ r = σ θ = σ z = 1 / 3 I1 = p
(5.34)
Volumetric Strains:
ε r = ε θ = ε z = 1 / 3 ε kk
(5.35)
Deviatoric Stresses:
s r = sθ = s z = s r z = 0
(5.36)
Deviatoric Strains:
er = eθ = e z = er z = 0
(5.37)
The relationship between the increment of hydrostatic pressure dp and the elastic volumetric strain increment dε kke can be expressed as usual
dp = K dε kke
(5.38)
On the other hand, the relation between the plastic volumetric strain increment dε kkp and
I 1 is expressed by Eqn. (5.24), where X (ε kkp ) is I 1 . Thus we can write 116
ε kkp = W exp( D I 1 ) − W
(5.39)
∂ ε kkp ∂ I 1 ∂I = DW exp( D I 1 ) ⋅ 1 ⋅ ∂p ∂ I1 ∂ p
(5.40)
1 dε kkp 3 D W exp( D I 1 )
(5.41)
dp =
Hence, we can write the relationship between the total volumetric strain increment dε kk and the hydrostatic pressure dp as
~ dp = K dε kkp
(5.42)
~
with the apparent bulk modulus K =
1 1 / K + 3 D W exp( D x)
(5.43)
The second term in the denominator of Eqn. (5.43) produces an apparent softening of the bulk modulus due to a plastic volumetric compaction (void ratio or porosity change). At high pressure levels, the softening term goes to zero, and the apparent bulk
~
modulus K approaches the elastic bulk modulus K . We subject a linear axisymmetric element to hydrostatic pressure and evaluate the stress-strain behaviour of the cap model at an integration point as shown in Figs. (5.7) and (5.8). The element is first isotropically consolidated with p = 0.1 N / mm 2 (point 1 to point 2), then completely unloaded (point 2 to point 3), and then reloaded to
p = 0.1 N / mm 2 again (point 3 to point 2). We see that the model dictates a purely elastic unloading-reloading behaviour.
117
Fig. 5.7 Stress path / cap movement in the meridional plane for hydrostatic stress
Stress-strain paths for isotropic consolidation
Slope = 3*K
Total volumetric strains eps kk [-]
0.15
Point 3
Point 2
0.1
Slope = 3*apparent K
0.05
Point 1 0
0
0.05
0.1 0.15 0.2 0.25 First stress invariant I1 = 3p [N/mm²]
0.3
Fig. 5.8 Evolution of volumetric strain under loading–unloading isotropic consolidation
118
5.4.4 Triaxial Shear Testing During triaxial testing, the sample is first isotropically consolidated by a hydrostatic pressure p = 0.05 N / mm 2 , which follows the path from point 1 to point 2 as shown in Fig 5.9. During the following shear phase, only the axial load is increased (the piston force - see section 4.3), which is the shear phase (point 2 to point 3). As soon as the stress path reaches the Drucker-Prager cone at p axial ≈ 0.07 N / mm 2 , the sample fails (see Fig. 5.10). With σ r , σ z and ε r , ε z being principal stresses and strains respectively, we can find
~
the apparent shear modulus G for the material under triaxial loading conditions from the ratio of stress difference increment to strain difference increment (Chen, Baladi 1985), (Bardet 1997) as follows
∂ (σ z − σ r ) ~ = 2G = 2G ∂ (ε z − ε r )
⎡ ∂ (ε zp − ε rp ) ⎤ ⎢1 − ⎥ ∂ (ε z − ε r ) ⎦ ⎣
(5.44)
The second term on the right-hand side generates a softening of the elastic shear modulus G due to plastic flow. Because of the elliptical cap shape, shear stress also contributes to compaction (see stress path in Fig 5.9).
In the chapter 6, we will use the elasto-plastic Drucker-Prager cap model, first, to compare the stress and displacement results of the refined constitutive relations to linear elasticity, second, to study to what extent elasto-plasticity improves stress and displacement solutions compared to experimental results of alumina granulate, and third, to find out what effect elasto-plasticity has on the prediction of the mobilization force.
119
Fig. 5.9 Stress path / cap movement in the meridional plane for triaxial shear testing
Stress path under triaxial shear testing 0.15
Total volumetric strains eps kk [-]
Point 3
Point 2
0.1
0.05
Point 1
0
0
Note: One element model fails after point 3 because of perfectly plastic material 0.05 0.1 0.15 First stress invariant I1 = 3p [N/mm²]
0.2
Fig. 5.10 Evolution of volumetric strain in triaxial shear testing up to failure
120
Chapter 6: Results and Discussion
6.1 Elastic Stresses and Displacements in a Confined Granular Column
Before moving on to the mobilization force for a granular column of finite depth, we present the exact elastic stress and displacement fields for a long granular column of alumina granulate using the analytical Fourier-Bessel method from chapter 3. The solution method holds for the general case of any elastic cylinder with Coulomb’s friction over its lateral surface and is not restricted to granular materials. We also present corresponding FE computations with linear elastic constitutive relations in order to back up the exact Fourier-Bessel solutions. Note that a checking for accuracy by the FE results is not appropriate, since the Fourier-Bessel solutions are analytically exact.
6.1.1 Exact Elasticity Solutions from the Fourier-Bessel Potential We compute stress and displacement fields for the upper 15 cm of long confined granular column. The results are graphically illustrated in section 6.1.3. For the input in the MATLAB routine in Appendix A, the following parameters are chosen:
-
Geometry:
Radius R = 3.75cm, Cylinder Length Lcyl = 15cm , Sampling interval [ 0; L = 20cm ]
121
-
Material:
E = 15 ⋅ 10 6 N / mm 2 , ν = 0.32, μ = 0.36 , (weightless alumina granulate)
-
Surcharge:
σ 0 = 0.1 N / mm 2
The question arises whether all boundary results apart from the specified boundary functions are maintainable from a physical point of view. From the plots we find the following:
τ rz ≠ 0
(6.1)
ur = 0
(6.2)
u z = f (r )
(6.3)
At r = R :
u z ( z) ≠ 0
(6.4)
At z = L cyl :
u z ( R) = A
(6.5)
At z = 0 :
Eqn. (6.1) is only surprising insofar that it is different from the well-known Janssen model, but necessary from the point of view of equilibrium (see Fig. 6.1). In Eqn. (6.2) it is assumed that the rigid platen prevents radial sliding of the elastic cylinder at any point of the top surface. It can be found as a specified boundary function
in
(Moghe,
Neff
1971)
and
(Ti-
moshenko, Goodier 1970).
Fig. 6.1 Boundary functions from Fourier-Bessel analysis
122
Eqn (6.3) must be considered as inconsistent. We find a considerable bending effect throughout the cylinder, which cannot be there at the top, since the rigid platen keeps the top surface flat. However, in finite element tests with constant axial top displacement, the bending effect started to re-appear right below the flat top. Therefore, we accept this inconsistency which is obviously limited to some of the upper portion of the cylinder. For a displacement-driven analysis, i. e. an axial displacement at the top is specified instead of the constant axial stress, the solution cannot be found with the method described in chapter 2, since the superposition in section 2.4.4 does not hold anymore. If there is a constant displacement at the top, the axial stress over the top surface is likely to be not constant. We recall from chapter 2.1 that in case of Coulomb’s friction particles at the cylinder surface are not allowed to move, because they are only about to slide. This condition is obviously not met in Eqn. (6.4), which shows clearly the limitation of Coulomb’s simple friction model. Experimental experience shows that an axial strain at the wall is absolutely reasonable. The material slides locally over the cylinder surface and after a short period of sliding, the deformed particles are im-
123
Fig. 6.2 Deformed axisymmetric elastic FE mesh (middle part not displayed)
mobile again. The axial displacement u z ( R, L cyl ) is not determinate by the boundary conditions Eqns. (3.87) – (3.89), which we know from Eqn. (3.54). In our case, we get some nonzero value A at (r , z ) = ( R, L cyl ) as stated in Eqn. (6.5). To make axial deformations comparable to the FE solution fields, we impose a rigid body motion on the column such that A = 0 .
5.1.2 Elasticity Solutions from Finite Element Analysis We want to check our analytical Fourier-Bessel results by comparison to corresponding FE solution fields. A long axisymmetric elastic column is therefore modelled with linear elastic constitutive relations in the FE software package ABAQUS (2005) as shown in Fig. 6.2. It can be described by the following parameters that are chosen in accordance with the Fourier-Bessel solutions from 6.1.1:
-
Geometry:
R = 3.75cm, L = 35 cm
-
Elements:
Axisymmetric, 4-node linear, solid mechanics
-
Spatial Resolution:
50 / 350 elements in radial / axial direction, respectively
-
Material:
Linear elasticity for small deformations ( E = 15 ⋅ 10 6 N / mm 2 , ν = 0.32 ), weightless (as obtained for alumina granulate)
-
Boundary Conditions:
see Fig. 6.2
-
Contact Constraint:
Tangential behaviour with friction ( μ = 0.36 ), Lagrangian multiplier (exact), master surface (undeformable wall) and slave surface (elastic elements), for details see (Be-
124
lytschko et al. 2000) and (Wriggers 2002)
The solution fields for the upper 15 cm of the FE model are graphically displayed in section 6.1.3. To make them comparable to the corresponding Fourier-Bessel solution, axial displacements are subjected to a rigid body motion such that they are zero at
(r , z ) = ( R, 15 cm) . The Fourier-Bessel and FE solution fields fit very well. Note that the finite element analysis predicts a singularity for all stress fields at the point where the top surface meets the wall [ (r , z ) = ( R, L = 0 cm) ]. This effect is minimized by the fine mesh to the area close to the corner, but can still be seen distinctly in the plots.
6.1.3 Result Plots The Fourier-Bessel and FEM stress and displacement fields are plotted as the third dimension over the axisymmetric plane (r , z ) (see Fig. 6.3).
Fig. 6.3 Axisymmetric plane over which the stress and displacement results are plotted as the third dimension
125
Fig. 6.4 Radial stress
126
σr
Fig. 6.5 Circumferential stress
127
σθ
Fig. 6.6 Axial stress
128
σz
Fig. 6.7 Shear stress
129
τ rz
Fig. 6.8 Radial displacement u r
130
Fig. 6.9 Axial displacement u z
131
6.2 Elasticity vs. Elasto-plastic and Experimental Material Behaviour
First, we compare elastic, elasto-plastic and experimental material behaviour for the simple case of a uniaxial strain test, whose principles have been discussed in section 4.3.2. From the results, we can evaluate the adequacy of elasticity for modelling stresses and displacements in alumina granulate.
6.2.1 The Uniaxial Strain Test The FEM model shown in Fig. 6.10 is run in ABAQUS with both linear elastic constitutive relations and the Drucker-Prager cap model from Appendix B, which is integrated by the user subroutine UMAT. The following model parameters are chosen:
-
Geometry:
R = 3.75cm, L = 5 cm
-
Elements:
Axisymmetric, 4-node linear, solid mechanics
-
Spatial Resolution:
20 / 25 elements in radial / axial direction, respectively
-
Material:
Linear elasticity ( E = 15 ⋅ 10 6 N / mm 2 , ν = 0.32 ) or Cap model with material parameters from 4.4
-
Boundary Conditions: see Fig. 6.10
-
Contact Constraint:
Tangential behaviour with friction ( μ = 0.36 ), Penalty method, master surface (undeformable wall) and slave surface (elastic elements), for details see (Belytschko et al. 2000) and (Wriggers 2002)
132
Since we deal only with FEA in the following, we can mend the bending inconsistency Eqn. (6.3) by exchanging the constant axial surcharge σ z at the top surface by a constant axial displacement u z = 0.5 mm . Under the assumption that the axial strain does not vary much along the small depth of 5 cm, this corresponds to an axial strain of around 1% , so that a small deformation Fig. 6.10 Axisymmetric finite element model for the uniaxial strain test
setting can still be tolerated.
6.2.2 Result Plots We plot stress and displacement results in Figs. 6.11 – 6.14 along the wall, the center of the cylinder as well as over cross sectional lines of different depths. From the experiments in section 4.3.2, the loading end displacement u zexp = 0.5 mm for a 5 cm depth of alumina granulate is known in relation to the piston force F . Under the assumption of constant stress at the loading end, we can compute σ zexp (r , 0) = const from F . Both experimental results u zexp and σ zexp are added for comparison in Figs. 6.12 and 6.14, respectively. Note that some stress results are multiplied by the factor 10 for better illustration as marked in the plots.
133
Radial stress over the top surface 0.16 Elastic Elastic-Plastic
0.14
Stress [N/mm²]
0.12 0.1 0.08 10 * enlarged
0.06 0.04 0.02 0
0
1.25
2.5
3.75
Radius r [cm]
Radial stress along the wall 0.16 Elastic Elastic-Plastic
0.14
Stress [N/mm²]
0.12 0.1 0.08 10 * enlarged
0.06 0.04 0.02 0
0
0.5
1
1.5 2 2.5 3 3.5 Depth from loading end z [cm]
Fig. 6.11 Radial stress
134
σr
4
4.5
5
Axial stress over the loading surface 0.4 Elastic Elasto-Plastic Experimental
Stress [N/mm²]
0.3
0.2
10 * enlarged 0.1 0.071
0
0
1.25
2.5
3.75
Radius r [cm]
Axial stress along the wall 0.4 Elastic Elastic-Plastic
0.35
Stress [N/mm²]
0.3 0.25 0.2 0.15 0.1
10 * enlarged
0.05 0
0
0.5
1
1.5 2 2.5 3 3.5 Depth from loading end z [cm]
Fig. 6.12 Axial stress
135
σz
4
4.5
5
Shear stress after third element layer 0,04 Elastic Elastic-Plastic
Stress [N/mm²]
0.03
0.02
10 * enlarged 0,01
0 -0.005
0
1.25
2.5
3.75
Radius r [cm]
Shear stress along the wall 0.04 Elastic Elastic-Plastic
Stress [N/mm²]
0.03
0.02
10 * enlarged
0.01
0
0
1
2 3 Depth from loading end z [cm]
Fig. 6.13 Shear stress
136
τ rz
4
5
Radial displacement halfways at L = 2.5 cm 0.03 Elastic Elastic-Plastic
Displacement [mm²]
0.02
0.01
0 0
1.25 2.5 Depth from loading end z [cm]
3.75
Axial displacement 0.6 Elastic at r = 0 Elasto-Plastic at r = 0 Elastic at r = R Elasto-Plastic at r = R
Experimental value at loading end = 0.5 mm
Displacement [mm]
0.5
0.4
0.3
0.2
0.1
0
0
1
2 3 Depth from loading end z [cm]
4
Fig. 6.14 Radial and axial displacements u r and u z
137
5
6.2.3 Discussion and Model Evaluation Fig. 6.12 clearly reveals that elasticity is not able to adequately model stresses and displacements at the same time in the confined alumina granulate. Whereas the displacements in Fig. 6.14 are somewhat reasonable (which is no surprise, since the matching loading end boundary condition u zexp is specified), the stresses overestimate the experimental value σ zexp by the factor 40. On the other hand, the Drucker-Prager cap model gives a set of far better stress and displacement results. In Fig. 6.12, its axial stress at the loading end lies within a few percent from the averaged σ zexp from the piston force. This can be explained by its capability to undergo deformation with much less stress increase because of hydrostatic softening Eqn. (5.43) and shear softening Eqn. (5.44). Accordingly, the stress-strain relations of the Drucker-Prager cap model can give a set of consistent stress and displacement solution fields, whereas elasticity can only approximate either displacements or stresses depending on whether the former or the latter is prescribed as a matching boundary condition.
The Drucker-Prager solution has some more interesting characteristics that, however, we cannot check with the alumina granulate, since we lack appropriate experimental data for back up: The stress path in the elements close to the wall must reach the Drucker-Prager yield surface. On the one hand, the shear stress produced by the wall friction is transformed into plastic shear deformation close to the wall (see. Fig. 6.13), which leads to much lower shear stresses further inside the cylinder than predicted by elasticity.
138
On the other hand, shear failure in the loading case (corner state) does not produce any volumetric strains (compaction) (see Fig 6.15). Recall in this context that strains depend on both the current state of stress and the loading history (stress path). Therefore, axial displacements are almost zero close to the wall, but large in the center of the cylinder, where shear stresses are zero (see Fig. 6.14). This means there is a considerable bending effect along the depth, although displacements are constant over the top and bottom cross sections. The frictional shear stress at the loading end is zero in both elastic and elastoplastic solutions (see Fig 6.13), which cannot be true, if we consider the exact elasticity solution in Fig. 6.7. This inconsistency is likely to originate in the finite element contact formulation, which goes beyond the scope of this work.
Fig. 6.15 Shear behaviour of elements close to the wall (stress path 1); Recall that principal strain and stress directions are assumed to be equivalent. The incremental plastic strain tensor on the yield surface has a component opposite to the compaction direction (volumetric expansion). In the unloading case, the cap can thus be drawn back by the reduction of plastic volumetric strains (dilatancy control). To which extent compaction and volumetric extension cancel depends primarily on the loading history. Stress path 1 has no compaction since the states of stress are always corner states that yield in pure shear. However, stress path 2 shows considerable compaction (see Fig. 4.8), although the current state of stress is equal to stress path 1.
139
Gudehus (1985) formulated four basic requirements from which the quality of a set of constitutive relations can be assessed: (1)
The mathematical formulation should result in a unique and stable stressstrain relationship.
(2)
This mathematical relationship should be defined only by a few parameters which can be determined from standard test data.
(3)
The mathematical model should encompass the well-known Mohr - Coulomb criterion as a special case.
(4)
The constitutive equations should reflect the key characteristics of experimental data.
On the one hand, the Drucker-Prager cap model satisfies all four requirements excellently: The cap algorithm in Appendix B is implicit and unconditionally stable (1), all important material parameters have been determined by standardized soils testing methods in sections 4.2 and 4.3 (2) and the Drucker-Prager yield cone has been matched to the Mohr-Coulomb criterion in section 5.2 (3). Most importantly, the cap model is able to model a variety of granular key characteristics (4), such as hydrostatic and shear softening, change of void ratio and unrecoverable deformation in the unloading case. On the other hand, elasticity does not have to be totally rejected. Linear elasticity can of course never fulfil (3), since there is no plasticity involved. However, (1) and (2) are always satisfied and requirement (4) depends on the specific engineering problem. In some applications, where only stresses are of interest and displacements are not considered, elasticity has been applied very successfully. The most well-known example
140
is the Janssen model from section 2.2, which is based on an equilibrium formulation only without considering any displacements or strains, and has been proven to predict fairly accurate stresses (no displacements!) in storage bins. If displacements are to be included, nonlinear elastic models have been applied successfully (Chen, Mizuno 1990). In the final section, we will inquire how accurately - compared to elasto-plastic and experimental results - an elastic model can predict the mobilization force for a confined granular column.
141
6.3 The Mobilization Force in a Confined Granular Column
We call the particular force, which sets a confined granular column in motion when applied on the top loading end, the mobilization force F mob . As explained in the introduction, we are looking for a method to approximate the mobilization force by a continuum theory. In section 3.6, the concept of a supporting bottom stress σ zs up has been introduced that represents the effect of arching. We compute the mobilization force by a simple Janssen model and compare the results to a computational elasto-plastic approach based on chapter 5 and the experimental results from section 4.4.
6.3.1 Mobilization Force and Continuum Modelling The mobilization force F mob depends on the frictional shear stress that is developed at the wall – grain interface. The larger the portion of force that can be transferred to the wall, the higher is the total mobilization force. In the elastic and elasto-plastic continuum models that we presented in chapter 3 and 5, the frictional shear stress depends primarily on two material characteristics:
(1)
The relation between axial and radial normal stress in the confined material, which is described by ν in linear elasticity and the K-value in the Janssen model.
(2)
The friction coefficient μ at the contact surface between grains and tube wall that we assumed to be constant (Coulomb friction).
142
We recall from section 3.6 the assumption
of
a
supporting
bottom
stress σ zs up that is provided to account for the arching effect at the bottom (see Figs. 3.18 and 6.16). To move the whole granular column, σ zs up has to be overcome by the bottom residual of the mobilization stress σ zload = F mob / A , which is ap-
Fig. 6.16 Finite length granular column and pushing top piston
plied on the loading end (see Fig. 6.16). σ zload has to be much larger than σ zs up , because it is constantly reduced by the counteracting wall friction, as it moves down the depth of the column. Therefore, the mobilization force is largely affected by the depth of the granular column, which can be backed up by experimental experience (see Fig. 4.15). We suggest the following simple test to approximate the mobilization stress σ zload : We choose a continuum model with zero displacement and frictional boundary conditions at r = R and the constant supporting bottom stress σ zs up at z = L (see Fig. 6.16). The constant surcharge σ z at z = 0 is stepwise increased, until equilibrium cannot be found anymore. From the last step, for which equilibrium could be satisfied, we find the corresponding surcharge σ zload and can compute the mobilization force F mob .
6.3.2 A Janssen and a FEM Approach This procedure can be evaluated either analytically by the Janssen model or computationally by the finite element method.
143
a) The analytical Janssen method Recall the Janssen formula Eqn. (2.12) from section 2.2. If we assume the material to be weightless as in the case of the alumina granulate we can write the expression for axial stress as
⎛ 2μ K ⎞ z⎟ R ⎠ ⎝
σ z = σ 0 exp⎜ −
(6.6)
We know that the axial stress σ z must be equal to σ zs up at z = L . If we insert this boundary condition into Eqn. (6.6) we find an expression for the mobilization surcharge
σ 0 = σ zload at z = 0 as
σ
load z
=σ
s up z
⎡ ⎛ 2 μ K ⎞⎤ ⎢exp⎜ − R L ⎟⎥ ⎠⎦ ⎣ ⎝
−1
(6.7)
For the case of a 15 cm column of alumina granulate confined in the plastic PMMA tube, we adopt the material parameters K = 0.47 and μ = 0.36 from section 4.5. The results are shown in Figs. 6.18 and 6.19.
b) The computational FEM method Consider an FEM model implemented in the software package ABAQUS (2005) as shown in Fig 6.17. The following model parameters are chosen:
-
Geometry:
R = 3.75cm, L = 15 cm
-
Elements:
Axisymmetric, 4-node linear, solid mechanics
-
Spatial Resolution:
75 / 25 elements in radial / axial direction, respectively 144
-
Material:
Cap model with material parameters from section 5.4
-
Boundary Conditions: see Fig. 6.17; the loading end is displacement driven
-
Contact Constraint:
Penalty Method ( μ = 0.36 )
Furthermore, we abandon the small deformation assumption, since strains will be much higher than the 1% margin from section 6.2. However, a large deformation setting is not compatible to the cap model formulation from chapter 5. Therefore, the elastic constitutive relations provided by the user subroutine UMAT are exchanged for the elasto-plastic DruckerPrager cap model that is provided by the ABAQUS material library. Apart from the nonassociated plastic potential, it works exactly the same way as our material code in Appendix B, but can additionally be applied for large deformations. The axial surcharge σ z at the loading end (which on our case is represented by a constant displacement u z ) is divided in several small loading steps. At a critical step, equilibrium cannot be found and the FE software gives out an error message.
145
Fig. 6.17 FEM for determination of the mobilization force
Axial stress over the top loading surface
0.05
FEM + Bottom support 0.010 FEM + Bottom support 0.006 FEM + Bottom support 0.002 Corresponding Janssen results
0.04
Janssen 0.0388 vs. averaged FEM 0.0380
Axial Stress [N/mm²]
0.06
0.03 Janssen 0.0233 vs. averaged FEM 0.0234 0.02
Janssen 0.0078 vs. averaged FEM 0.087
0.01
0
0
0.5
1
1.5 2 Radius r [cm]
2.5
3
3.5
Fig. 6.18 Mobilization stress at the loading end for 3 different bottom support stresses (0.002, 0.006, 0.010 N/mm²) determined by FEM and Janssen
-3
7
Shear stress along the wall
x 10
FEM Janssen 0.010 N/mm² 0.006 N/mm² 0.002 N/mm²
6
Shear Stress [N/mm²]
5
4
3
2
1
0
0
3
6 9 Depth from loading end z [cm]
12
15
Fig. 6.19 Frictional shear stresses along the wall for 3 different bottom support stresses (0.002, 0.006, 0.010 N/mm²) determined by FEM and Janssen
146
Janssen, FEM and Experimental Results 1
Top loading end stress [N/mm²]
FEM Janssen Experimental
0.1
0.05
0
0
0.002
0.004 0.006 0.008 Supporting bottom stress [N/mm²]
0.01
0.012
Fig. 6.20 Janssen, FEM and experimental results showing the linear dependence of the mobilization stress at top from the bottom support stress
Mobilization Stress for a Bottom Support of 0.005 N/mm² 0.25
Mobilization stress at loading end [N/mm²]
Janssen Experimental 0.2
0.15
0.1
0.05
0
7.5
15 Depth of the column [cm]
22.5
Fig. 6.21 Janssen and experimental results showing the exponential dependence of the mobilization stress at top from the column length
147
We decrease the top displacement increment by trial and error to accurately find the transition point and give out the corresponding mobilization stress σ zload for the last stable step. From there, the mobilization force F mob = σ zload A can be computed easily. Note that we still use ABAQUS standard (implicit FEM). The results are paired with the Janssen solution and the experimental results from section 4.4 in Figs. 6.18 - 6.21.
6.3.3 Discussion and Evaluation Despite the controversial results of elasticity modelling the uniaxial strain test (see 6.3.2), we find a surprisingly accurate accordance of the Janssen model with the elastoplastic FEM in the mobilization test (see Figs. 6.18 – 6.21). Since the frictional shear stresses of the Janssen model and FEM are very close, the mobilization stresses σ zload are consistent as well. Note again that the shear stress are zero at z = 0 in the FE model, which is due to the contact formulation (see 6.3.2 for a short discussion).
The predicted mobilization stresses are -
linearly dependent on the supporting bottom stress (see Fig. 6.20).
-
exponentially dependent on the length of the bed (Fig. 6.21).
This behaviour is at least qualitatively supported by the experimental results from the mobilization tests on alumina granulate (see section 4.4). However, the experimental values of the mobilization stresses are about 3 to 5 times higher than the elastic and elasto-plastic results.
148
This discrepancy can be due to several reasons:
(1)
The assumption of a small supporting bottom stress is not appropriate for modelling arching effects (see sections 2.3 and 3.6).
(2)
There is another mechanism involved that we have failed to discover so far.
(3)
The alumina granulate gives inconsistent results.
(4)
Coulomb’s friction model is inappropriate for modelling the frictional behaviour between the grains and the wall.
(5)
The mobilization testing method does not represent the conditions assumed in the modelling.
(6)
Continuum mechanics and its constitutive relations are unable to model granular material behaviour in limit states.
(7)
The material parameter tests in chapter 4 are flawed and give wrong material parameters.
The last three last arguments (5), (6) and (7) can be rejected. First, there are several independent research groups that published consistent results for the mobilization force tested equivalently to the method presented in section 4.4 (see references there). As for (6), we can refer to soil mechanics where continuum models have been used extensively. Their application is not restricted to static problems, but they can also account for limit state calculations such as the slope failure in embankments (Chen, Mizuno 1990), (Potts, Zdravkovich 1999). Therefore, there is no reason why the Drucker-Prager
149
cap model is to be blamed for such big deviations in our case. Furthermore, the consistent results of the cap model for the uniaxial stain test indicate appropriate material parameters for the alumina granulate (see Fig. 6.12). Hence, we have to concentrate on arguments (1) – (4). On the one hand, experimental experience shows that the alumina granulate indeed gives inconsistent results with variation up to 30%, if experimental procedures are repeated in the same day. This scattering is most likely due to the discrete nature of the material in mind, which affects the repeatability of tests as discussed in section 4.1. On the other hand, a scattering up to 150% has been observed in equivalent tests that were carried out on different days. Alumina granulate is extremely sensitive to water and it starts to react heavily if brought into contact with it. This gives rise to the suspicion that the moisture content (humidity) that varied considerably between the test dates affects the material parameters thus amplifying the scattering effect. Since all material parameters for the cap model were carried out in the winter, when humidity is stable and only varying slightly, we get consistent results in section 6.2. However, it is not reasonable to blame exclusively the material for the large deviation between theory and experiment. But concerning the adequacy of the simplification that is included in the assumption of the small supporting bottom stress for modelling arching effects, we can hardly more than guess at this point. Therefore, we propose a series of further investigations on arguments (1) – (4) in the next chapter.
150
Chapter 7: Conclusions and Recommendations
As presented in the first chapter, the two main objectives of this work are, first, to clarify the existence of the “zero-axial-stress” problem in a confined elastic medium and, second, to study the mobilization problem.
7.1 Summary of the Work
The analytical, computational and experimental methods and techniques that we used to study these two aspects are briefly summarized in the following:
•
For the “zero-axial-stress” problem, we chose an analytical Fourier-Bessel potential approach, which was purely elastic. The Coulomb’s friction condition at the radial boundary could be well approximated by the assumption of Janssen shear stress. Furthermore, we presented an improvement algorithm to iterate the solution fields towards the exact Coulomb’s condition.
•
We introduced a small bottom support stress to represent the arching effect at the bottom of the granular column for continuum modelling.
•
We briefly derived the elasto-plastic Drucker-Prager cap model that had been developed in geomechanics, and outlined its implementation into the finite element
151
method. We studied the characteristics of elasto-plastic material behaviour in a simple axisymmetric element.
•
For experimental verification, we chose the example of a column of alumina granulate confined in a plastic tube. The elastic and elasto-plastic material parameters were derived by standardized testing methods from soil mechanics. To determine the friction coefficient at the grain – tube wall interface, we proposed a special in situ friction testing procedure.
•
With the help of the bottom support stress representing the arching effect, the mobilization force could predicted analytically by the Janssen model and computationally by FEM. We also proposed an “upside down” testing method to verify the mobilization force experimentally.
7.2 Conclusions of the Research
We presented and discussed the results of the analytical, computational and experimental models in chapter 6. From there, we can draw the following conclusions:
•
The Fourier-Bessel method provided the exact stress and displacement solutions for a confined elastic cylinder that is infinitesimally long in z -direction. We found out that neither of the stress and displacement fields is constant over any cross sec-
152
tional area except for the top surface of the column where constant axial stress was specified as a boundary condition.
•
The “zero-axial-stress” problem is not possible in continuum mechanics. The phenomenon of arching in a finite confined granular column cannot be represented in a continuum model in the context of our set of boundary conditions.
•
Axial and radial stress fields and the shear stress distribution over the radial boundary that are predicted by the Janssen model are very close to the corresponding exact solution fields for a confined elastic medium with Coulomb’s friction over the radial boundary. We can assess the accuracy of the Janssen model from the point of view of an exact linear elasticity solution by comparing μ close given by the Janssen shear stress assumption to the actual constant Coulomb’s friction coefficient μ .
•
The exact K-value for elasticity is a function of z and can be determined by the improvement algorithm.
•
Linear elasticity is not appropriate for modelling both stress and corresponding displacement fields together, since it gives either appropriate stresses with corresponding displacements that are too small, or appropriate displacements with corresponding stresses that are too high. However, if only one kind of solution fields is taken into consideration, i. e. either only stresses or only displacements, linear elas-
153
ticity is able to approximate stresses or displacements in granular materials very well. As an example, we studied the Janssen model that considers stresses only.
•
Elasto-plasticity is an excellent tool for modelling stress fields together with the corresponding displacement fields, since it can account for plastic deformations also. Elasto-plastic results for the uniaxial strain test match very well with experimental verification.
•
The mobilization force could be predicted analytically (Janssen), computationally (FEM) and experimentally (“upside down” testing for confined alumina granulate), which yield qualitatively concordant results. On the one hand, the mobilization force is linearly dependent on the strength of the bottom arch (represented by the bottom support stress in continuum modelling and the top weight in the experimental “upside down” testing). On the other hand, the mobilization force is exponentially dependent on the length of the column.
•
Linear elasticity (analytical Janssen model) and elasto-plasticity (FEM with Drucker-Prager cap material) give very close results that differ less than 1% from each other. This is surprising, since we have seen that stress and displacement fields are very different from each other.
•
Alumina granulate is not the ideal material for verification purposes, since it shows a considerable interference by environmental conditions. For the same testing under the same conditions, the results differ about 150% when tested at different days. We
154
believe that this is caused by the sensitivity to water, so that the material properties of alumina granulate change with the moisture content of the air (humidity).
7.3 Recommendations for Future Work
Whereas the mathematical formulations of both the analytical elastic Fourier-Bessel model and the computational elasto-plastic Drucker-Prager model have been examined thoroughly in the scope of this work, the experimental verification of these theoretical methods is still uncompleted. We could successfully gather experimental data from alumina granulate samples to show that elasto-plasticity is far more appropriate for modelling stress and displacement fields than linear elasticity (uniaxial strain test). However, we could not obtain reliable data for the mobilization force in a confined column of alumina granulate. We recall from section 6.3.3 that the reasons for this could not be clarified definitely, because we could not estimate the following hypotheses at this point:
(1)
The assumption of a small supporting bottom stress is not appropriate for modelling arching effects (see sections 2.3 and 3.6).
(2)
There is another mechanism involved that we have failed to discover so far.
(3)
The example material (alumina granulate) gives inconsistent results.
(4)
Coulomb’s friction model is inappropriate for modelling the frictional behaviour between the grains and the wall.
155
In order to make appropriate statements on the influence of (1), (2), (3) and (4), we suggest the following future work to be undertaken:
•
Apart from continuum mechanics there are other theories that give a different approach to stresses and displacements in granular materials. Some of these theories such as statistical mechanics and computational discrete elements have been briefly described in section 2.3. Research groups in both fields are currently working on methods to describe arching effects. Their application to the mobilization problem will certainly give new insights in problems related to (1) and (2) and provide valuable data for comparison to continuum results obtained in this work.
•
We recommend the repetition of the material testing in chapter 5 for materials that are easy to handle, little affected by environmental influences and are known to be appropriate for homogenization. Thus a reliable estimation of argument (3) can be given. Suitable material examples are sieved Ottawa sand or, in the best case, glass spheres.
•
The adequacy and accuracy of the Coulomb’s friction model that we applied at the grain – tube wall interface should be checked. Some research on the frictional behaviour of grains on plastic surfaces has already been carried out in the field of tribology. Further experimental investigations on confined columns of various granular materials can also provide some insight on the real tribological behaviour. The results can be used to assess the influence of argument (4).
156
•
Additionally, other testing procedures for the mobilization force should be explored and set up to definitely exclude argument an influence of the “upside down” testing procedure as presented in section 4.4.
Unfortunately, the considerable discrepancy between real-world behaviour and theoretical predictions makes it impossible at this point to evaluate the method that we presented for the prediction of the mobilization force. With more experimental and theoretical experience and knowledge from these research suggestions, it will become also possible to give a detailed and sound evaluation on the adequacy and accuracy of the methods presented within this work.
157
Appendix A: MATLAB routine for the Fourier-Bessel model % Fourier-Bessel solution and polynomial superposition % Closed-form solution for the elastic cylinder problem % % % %
Dominik Schillinger Department of Civil & Environmental Engineering University of Connecticut May 9th, 2006
syms r z A B C D alpha % Geometry, Material and Load R=3.75; L=20; P=0.1; nu=0.32; N=50;
% % % % %
Radius Sampling interval (Depth of the cylinder +33%) Stress at z=0 (Top surcharge) Poisson’s ratio Number of Fourier series terms
% --------------------------------------------------------------------% Fourier-Bessel potential % --------------------------------------------------------------------% Janssen Shear Assumption mu = 0.36; % Coulomb’s coefficient K = nu/(1-nu); % K-value shear = -Kvalue*mu*P*exp(-2*mu*Kvalue/R*z);
% Janssen shear stress
% Orthogonality conditions for trigonometric functions for n=1:N a(n) = (2*n-1)/2*pi/L; end % Fourier-Bessel Potential (1st Love’s Strain Function) Z1 = cos(alpha*z)*[A*besseli(0,alpha*r)+B*alpha*r*besseli(1,alpha*r)]; % Taking the Laplacian in cylindrical coordinates del = diff(diff(Z1,r),r)+1/r*diff(Z1,r)+diff(diff(Z1,z),z);
% Getting the stresses srr1 sth1 szz1 srz1
= = = =
simplify(diff(nu*del-diff(diff(Z1,r),r),z)); simplify(diff(nu*del-1/r*diff(Z1,r),z)); simplify(diff((2-nu)*del-diff(diff(Z1,z),z),z)); simplify(diff((1-nu)*del-diff(diff(Z1,z),z),r));
158
% Getting the displacements urr1 = simplify(-diff(diff(Z1,r),z)); uzz1 = simplify(2*(1-nu)*del-diff(diff(Z1,z),z)); % Solving for A~B relation S = simplify(subs(solve(urr1/sin(alpha*z),A),r,R));
% Finite Fourier Cosine Transform for B for n=1:N G =double(simplify(subs(srz1/cos(a(n)*z), (next line) {A1,alpha,r},{S1,a(n),R})/B)); B1(n) = double(1/G * 2/L * int(shear*cos(a(n)*z),z,0,L)); end
% Determining the constant D for n=1:N S1(n) = double(subs(S,{B,alpha,r},{B1(n),a(n),R})); end
% Summing up the first 50 series terms % Radial stress Sigrr=0; for n=1:N Sigrr = Sigrr + simplify(subs(srr1,{A,B,alpha},{S1(n),B1(n),a(n)})); end % Circumferential stress Sigth=0; for n=1:N Sigth = Sigth + simplify(subs(sth1,{A,B,alpha},{S1(n),B1(n),a(n)})); end % Axial stress Sigzz=0; for n=1:N Sigzz = Sigzz + simplify(subs(szz1,{A,B,alpha},{S1(n),B1(n),a(n)})); end % Shear stress Sigrz=0; for n=1:N Sigrz = Sigrz + simplify(subs(srzn,{A,B,alpha},{S1(n),B1(n),a(n)})); end % Radial displacement Disrr=0; for n=1:N Disrr = Disrr + simplify(subs(urr1,{A,B,alpha},{S1(n),B1(n),a(n)}));
159
end % Axial displacement Diszz=0; for n=1:N Diszz = Diszz + simplify(subs(uzz1,{A,B,alpha},{S1(n),B1(n),a(n)})); end
% --------------------------------------------------------------------% Polynomial potential % ---------------------------------------------------------------------
% Love Strain Function (polynomials from Legendre's equation) Z2 = C*(z^3 - 3/2*z*r^2) + D*(r^2 + z^2)*z; % Taking the Laplacian del = diff(diff(Z2,r),r)+1/r*diff(Z2,r)+diff(diff(Z2,z),z);
% Getting the stresses szz2 = simplify(diff((2-nu)*del-diff(diff(Z2,z),z),z)) srr2 = simplify(diff(nu*del-diff(diff(Z2,r),r),z))
% Getting the displacements uzz2 = simplify(2*(1-nu)*del-diff(diff(Z2,z),z)) urr2 = simplify(-diff(diff(Z2,r),z))
% Evaluating constants with boundary conditions S2 = simplify(subs(solve(urr2,C),r,R)) S3 = simplify(solve(subs(szz2-P,{A,z},{S2,L}),D))
% Getting the polynomial stress and displacement fields srr sth szz uzz
% % % %
= = = =
double(subs(srr2,{C,D},{S2,S3})) double(subs(srr2,{C,D},{S2,S3})) double(subs(szz2,{C,D},{S2,S3})) simple(subs(uzz2,{C,D},{S2,S3}))
--------------------------------------------------------------------Superposition of Fourier-Bessel and polynomial solutions to the final stress and displacements fields ---------------------------------------------------------------------
Sigrr Sigzz Diszz Disth
= = = =
Sigrr-srr; Sigzz-szz; Diszz-uzz; Disth-sth;
160
Appendix B: Fortran user subroutine UMAT for the DruckerPrager cap model
C C C C C
*********************************************************************************** ** ABAQUS USER SUBROUTINE - UMAT DPCMDL (DRUCKER-PRAGER CAP MODEL) ** *********************************************************************************** ***********************************************************************************
1 2 3 4
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, RPL,DDSDDT,DRPLDE,DRPLDT, STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C INCLUDE 'ABA_PARAM.INC' C CHARACTER*80 CMNAME C C DOUBLE PRECISION STRESS(4),STATEV(2),DDSDDE(4,4),DSTRAN(4), DOUBLE PRECISION PROPS(6) C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C
NPROPS = INDEX OF MATERIAL CONSTANTS IN ABAQUS INPUT FILE (PROPS) --------------------------------------------------------------------NPROPS --------------------------------------------------------------------1 E (ELASTIC DOMAIN) 2 NU (ELASTIC DOMAIN) 3 ALPHA (COHESIONLESS DRUCKER-PRAGER, K=0) 4 AD (CAP) 5 AW (CAP) 6 R (CAP) ---------------------------------------------------------------------NSTATV = SOLUTION DEPENDENT STATE VARIABLES (STATEV) --------------------------------------------------------------------NSTATV --------------------------------------------------------------------1 EL CAP LOCATION (INTERSECTION CAP/FAILURE SURFACE) 2 MSTATE MATERIAL STATE = 0 ... ELASTIC STATE = 1 ... YIELDING STATE = 2 ... HARDENING STATE = 3 ... CORNER STATE ----------------------------------------------------------------------
CALL SUBROUTINE CAPMDL CALL CAPMDL(STRESS,DSTRAN,STATEV,DDSDDE,PROPS)
C RETURN END C C C C C C C C C
SUBROUTINE CAPMDL(STRESS,DSTRAN,STATEV,EPCM) *********************************************************************************** ** MAIN SUBROUTINE CAPMDL COMPUTES A CONSTITUTIVE MATRIX FOR CAP MODEL ** ** UNDER PLANE STRAIN CONDITION ** *********************************************************************************** ***********************************************************************************
161
INCLUDE 'ABA_PARAM.INC' C DOUBLE PRECISION DOUBLE PRECISION DOUBLE PRECISION DOUBLE PRECISION DOUBLE PRECISION DOUBLE PRECISION DOUBLE PRECISION INTEGER NIT
EV,DER,DEZ,DET,DERZ,PROPS(6) EL,E,NU,ALPHA,AD,AW,R,MSTATE STRESS(4),DSTRAN(4),EPCM(4,4),STATEV(2) XLL,BMOD,SMOD,THREEK,TWOG INCRNT,CONV,ELCTS,TSI1E,TSI1 TSCTS,ELL,CAPL,FL,FR,XR,XL,ELR SIGR,SIGZ,SIGT,SIGRZ,SLOPE
C EL=STATEV(1) E=PROPS(1) NU=PROPS(2) ALPHA=PROPS(3) AD=PROPS(4) AW=PROPS(5) R=PROPS(6) C INCRNT=10.D0 NIT=60 CONV=0.001 C C C
COMPUTE ELASTIC MODULI AND MATERIAL PROPERTIES BMOD =E/(3.D0*(1.D0-2.D0*NU)) SMOD =E/(2.D0*(1.D0+NU)) THREEK=3.D0*BMOD TWOG =2.D0*SMOD
C C C
CHANGE ENGINEERING SHEAR STRAIN TO TENSORIAL STRAIN DSTRAN(4)=0.5*DSTRAN(4)
C C C
COMPUTE SUBDIVIDED STRAIN INCREMENTS FOR GIVEN STRAIN INCREMENTS DSTRAN(1)=DSTRAN(1)/INCRNT DSTRAN(2)=DSTRAN(2)/INCRNT DSTRAN(3)=DSTRAN(3)/INCRNT DSTRAN(4)=DSTRAN(4)/INCRNT
C C C
DO LOOP FOR COMPUTATION OF STRESS STATE FOR NEW STRAIN STATE DO 170 INC=1,INCRNT
C C C
COMPUTE STRESS INVARIANTS AND DEVIATORIC STRESSES AT PREVIOUS STEP SI1I =STRESS(1)+STRESS(2)+STRESS(3) PRESS =SI1I/3.D0 DSIGR =STRESS(1)-PRESS DSIGZ =STRESS(2)-PRESS DSIGT =STRESS(3)-PRESS DSIGRZ=STRESS(4)
C C C
COMPUTE VOLUMETRIC STRAIN INCREMENT AND DEVIATORIC STRAIN INCREMENTS EV =DSTRAN(1)+DSTRAN(2)+DSTRAN(3) DER =DSTRAN(1)-EV/3.D0 DEZ =DSTRAN(2)-EV/3.D0 DET =DSTRAN(3)-EV/3.D0 DERZ=DSTRAN(4)
C C C
COMPUTE ACCUMULATED PLASTIC VOLUMETRIC STRAIN XL =XV(EL,R,ALPHA) EVPI =EVP(XL,AD,AW)
C C C
COMPUTE ELASTIC DEVIATORIC STRESSES AND STRESS INVARIANTS TDSIGR =DSIGR +TWOG*DER TDSIGZ =DSIGZ +TWOG*DEZ
162
1 C C C
TDSIGT =DSIGT +TWOG*DET TDSIGRZ=DSIGRZ+TWOG*DERZ TSI1 =SI1I+THREEK*EV TSRJ2 =SQRT(TDSIGR*TDSIGR+TDSIGT*TDSIGT+TDSIGR*TDSIGT+ TDSIGRZ*TDSIGRZ) SET INITIAL MSTATE AND SCALING FACTOR RATIO =1.D0 MSTATE=0
C C C
CHECK FOR HARDENING CODING CAPL=AMIN1(0.D0,EL) IF(TSI1.LT.CAPL) GO TO 70
C C C C C
****** FAILURE SURFACE CODING ******j CHECK FOR FAILURE SURFACE VMISES=SRJ2C(CAPL,XL,CAPL,R) FI1 =FF(TSI1,ALPHA) FD =TSRJ2-AMIN1(FI1,VMISES) IF(FD.LE.0) GO TO 160
C C C
FAILURE SURFACE CALCULATION DFDI1=0.D0 TSCTS=TSI1+CONV*TSRJ2 IF (FI1.LT.VMISES) DFDI1=(FI1-FF(TSCTS,ALPHA))/(CONV*TSRJ2) DEVP =3.D0*DFDI1*FD/(3.D0*THREEK*DFDI1**2+0.5*TWOG) TSI1E=TSI1 TSI1 =TSI1-THREEK*DEVP
C C C
SET MSTATE FOR YIELD SURFACE ZONE MSTATE=1
C C C
20
30
C C C 70 C C C 80
****** DILATANCY AND CORNER CODING ****** IF(EL.LT.0) GO TO 20 TSI1 = AMAX1(TSI1,CAPL) GO TO 30 ELCTS =EL+CONV*FF(EL,ALPHA) ELL =AMIN1(TSI1E,ELCTS) XLL =XV(ELL,R,ALPHA) DLDEVP=(ELL-EL)/((XLL-XL)*AD*AW*EXP(AD*XL)) ELINT =EL ELDD =ELINT+DEVP*DLDEVP EL =AMIN1(0.D0,ELDD) IF(TSI1.GE.EL) GO TO 30 TSI1 =(DLDEVP*TSI1E+THREEK*ELINT)/(DLDEVP+THREEK) EL =AMIN1(0.D0,TSI1) IF(TSI1.LE.0) MSTATE=3 FI1 =FF(TSI1,ALPHA) SLOPE =ALPHA IF(VMISES.LT.FI1) SLOPE=0.D0 RATIO =AMIN1(FI1,VMISES)/TSRJ2 GO TO 160 ****** HARDENING CODING FOR ELLIPTIC CAP SURFACE ****** IF(TSI1.LT.XL) GO TO 80 IF(TSRJ2.LE.SRJ2C(TSI1,XL,CAPL,R)) GO TO 160 INITIAL SET OF CAP LOCATION TSI1E =TSI1 TSRJ2E=TSRJ2 ELL =EL ELR =TSI1E
163
1
C C C
100
110
120
130 C C C
140 C C C
IF(TSI1E.LE.XL) FL=(EL-TSI1E)/(EL-XL) IF(TSI1E.GT.XL) FL=2.D0*TSRJ2E/ (TSRJ2E+SRJ2C(TSI1E,XL,CAPL,R))-1.D0 XR =AMIN1(0.D0,XV(ELR,R,ALPHA)) TSI1R =TSI1E - THREEK*(EVP(XR,AD,AW)-EVPI) FR =(XR-TSI1R)/(ELR-XR) COMP =CONV*FF((FL*XR-FR*XL)/(FL-FR),ALPHA) IF((ABS(TSI1)+TSRJ2).LT.COMP) GO TO 160 MSTATE=2 FOLD =0.D0 EMPLOY A MODIFIED REGULA FALSI METHOD TO DETERMINE CAP LOCATION
DO 130 IT=1,NIT EL =(FL*ELR-FR*ELL)/(FL-FR) XL =AMIN1(0.D0,XV(EL,R,ALPHA)) DEVP =EVP(XL,AD,AW)-EVPI TSI1 =TSI1E-THREEK*DEVP CAPL =AMIN1(0.D0,EL) IF(TSI1.LT.XL) FC=(EL-TSI1)/(EL-XL) IF(TSI1.GT.CAPL) FC=(XL-TSI1)/(CAPL-XL) IF(TSI1.LT.XL.OR.TSI1.GT.CAPL) GO TO 110 TSRJ2=SRJ2C(TSI1,XL,CAPL,R) DELI1=CONV*(XL-TSI1) DESP =0.D0 IF(DELI1.EQ.0) GO TO 100 DESP =(DEVP/6.D0)*R*SQRT((XL-CAPL)**2-(TSI1-CAPL)**2) 1 /(TSI1-CAPL) SRJ2T=TSRJ2+TWOG*DESP FC =(TSRJ2E-SRJ2T)/(TSRJ2E+SRJ2T) IF(ABS(TSRJ2E-SRJ2T).LE.COMP) GO TO 140 IF(FC.GT.0.AND.(TSI1-CAPL).GE.DELI1) GO TO 140 IF(FC.GT.0) GO TO 120 ELR=EL FL =FC IF(FOLD.LT.0) FL=0.5*FL GO TO 130 ELL=EL FL =FC IF(FOLD.GT.0.) FR=0.5*FR FOLD=FC IF THE MODIFIELD REGULA FALSI METHOD DOES NOT CONVERGE, SET A FINAL CAP TSI1=AMAX1(TSI1,XL) IF(TSI1.GT.AMIN1(0.D0,EL)) TSI1=CAPL TSRJ2=SRJ2C(TSI1,XL,CAPL,R) RATIO=0.D0 CHECK CORNER SRJ2T=FF(TSI1,ALPHA) ERR =TSRJ2/SRJ2T IF(ERR.GE.0.99 .AND.ERR.LE.1) MSTATE=3 IF(TSRJ2E.EQ.0.) GO TO 160 RATIO=TSRJ2/TSRJ2E
C C C C C 160
C C C
****** COMPUTE STRESS STATE CORRESPONDING TO NEW STRAIN STATE ****** SCALE BACK ELASTIC DEVIATORIC TRIAL STRESS DSIGR =RATIO*TDSIGR DSIGZ =RATIO*TDSIGZ DSIGT =RATIO*TDSIGT DSIGRZ=RATIO*TDSIGRZ COMPUTE FINAL STRESSES PRESS=TSI1/3.D0 SIGR =PRESS+DSIGR SIGZ =PRESS+DSIGZ
164
SIGT =PRESS+DSIGT SIGRZ=DSIGRZ C C C
COMPUTE ACCUMULATED PLASTIC VOLUMETRIC STRAIN XL =XV(EL,R,ALPHA) TEVP=EVP(XL,AD,AW)
C C C
INITIALIZE STRESS STATE AS PREVIOUS STRESS STATE STRESS(1)=SIGR STRESS(2)=SIGZ STRESS(3)=SIGT STRESS(4)=SIGRZ CONTINUE
170 C C C
UPDATE SOLUTION DEPENDENT VARIABLES STATEV(1)=EL STATEV(2)=MSTATE
C C C
CONSTRUCT A NEW CONSTITUTIVE MATRIX
1 1
C C C C C C C C
IF(MSTATE.EQ.0) CALL LINELS(BMOD,SMOD,EPCM) IF(MSTATE.EQ.1) CALL EPDP(SIGR,SIGZ,SIGT,SIGRZ,BMOD,SMOD, SLOPE,EPCM) IF(MSTATE.EQ.2.OR.MSTATE.EQ.3) CALL EPCAP(SIGR,SIGZ,SIGT,SIGRZ, EL,R,AD,AW,TEVP,BMOD,SMOD,ALPHA,EPCM) RETURN END
*********************************************************************************** ** UTILITY SUBROUTINES AND FUNCTIONS FOR MAIN SUBROUTINE CAPMDL ** *********************************************************************************** ***********************************************************************************
SUBROUTINE LINELS(BMOD,SMOD,EPCM) C *********************************************************************************** C ** UTILITY SUBROUTINE LINELS COMPUTES LINEAR ELASTIC CONSTITUTIVE MATRIX ** C ** (UMAT JACOBIAN) UNDER PLANE STRAIN CONDITION ** C *********************************************************************************** C INCLUDE 'ABA_PARAM.INC' C DOUBLE PRECISION EPCM(4,4),BMOD,SMOD C C COMPUTE COMPONENTS OF MATRIX EPCM C EPCM(1,1)=BMOD+4.D0/3.D0*SMOD EPCM(1,2)=BMOD-2.D0/3.D0*SMOD EPCM(1,3)=BMOD-2.D0/3.D0*SMOD EPCM(1,4)=0.D0 EPCM(2,1)=EPCM(1,2) EPCM(2,2)=EPCM(1,1) EPCM(2,3)=EPCM(1,3) EPCM(2,4)=0.D0 EPCM(3,1)=EPCM(1,2) EPCM(3,2)=EPCM(1,3) EPCM(3,3)=EPCM(1,1) EPCM(4,1)=0.D0 EPCM(4,2)=0.D0 EPCM(4,3)=0.D0 EPCM(4,4)=SMOD RETURN END C C C
165
C C C C C
SUBROUTINE EPDP(SIGR,SIGZ,SIGT,SIGRZ,BMOD,SMOD,ALPHA,EPCM) *********************************************************************************** ** UTILITY SUBROUTINE EPDP COMPUTES ELASTO-PLASTIC CONSTITUTIVE MATRIX FOR ** ** (UMAT JACOBIAN)THE DRUCKER-PRAGER MATERIAL UNDER PLANE STRAIN CONDITION ** *********************************************************************************** INCLUDE 'ABA_PARAM.INC'
C DOUBLE PRECISION EPCM(4,4),SI1,SJ2,SRJ2,H,HR,HZ,HT,HRZ DOUBLE PRECISION DSIGR,DSIGZ,DSIGT,DSIGRZ,BMOD,SMOD DOUBLE PRECISION SIGR,SIGZ,SIGT,SIGRZ,ALPHA C C C
COMPUTE STRESS INVARIANTS
1 C C C
SI1 =SIGR+SIGZ+SIGT SJ2 =((SIGR-SIGZ)**2+(SIGZ-SIGT)**2+(SIGT-SIGR)**2)/6.D0 +SIGRZ**2 SRJ2=SQRT(SJ2) COMPUTE DEVIATORIC STRESSES DSIGR =SIGR-SI1/3.D0 DSIGZ =SIGZ-SI1/3.D0 DSIGT =SIGT-SI1/3.D0 DSIGRZ=SIGRZ
C C C
****** COMPUTE H,HX,HY,HZ,HXY FOR DRUCKER-PRAGER MATERIAL ****** H =9.D0*BMOD*ALPHA**2+SMOD HR =3.D0*BMOD*ALPHA+SMOD/SRJ2*DSIGR HZ =3.D0*BMOD*ALPHA+SMOD/SRJ2*DSIGZ HT =3.D0*BMOD*ALPHA+SMOD/SRJ2*DSIGT HRZ=SMOD/SRJ2*DSIGRZ
C C C
COMPUTE COMPONENTS OF MATRIX EPCM EPCM(1,1)=BMOD+4.D0/3.D0*SMOD-HR*HR/H EPCM(1,2)=BMOD-2.D0/3.D0*SMOD-HZ*HZ/H EPCM(1,3)=BMOD-2.D0/3.D0*SMOD-HT*HT/H EPCM(1,4)=-HR*HRZ/H EPCM(2,1)=EPCM(1,2) EPCM(2,2)=BMOD+4.D0/3.D0*SMOD-HZ*HZ/H EPCM(2,3)=BMOD-2.D0/3.D0*SMOD-HT*HT/H EPCM(2,4)=-HZ*HRZ/H EPCM(3,1)=EPCM(1,3) EPCM(3,2)=EPCM(2,3) EPCM(3,3)=BMOD+4.D0/3.D0*SMOD-HT*HT/H EPCM(3,4)=-HT*HRZ/H EPCM(4,1)=EPCM(1,4) EPCM(4,2)=EPCM(2,4) EPCM(4,3)=EPCM(3,4) EPCM(4,4)=-HRZ*HRZ/H RETURN END
C C C SUBROUTINE EPCAP(SIGR,SIGZ,SIGT,SIGRZ,EL,R,AD,AW,TEVP,BMOD, 1 SMOD,ALPHA,EPCM) C *********************************************************************************** C ** UTILITY SUBROUTINE EPCAP COMPUTES ELASTO-PLASTIC CONSTITUTIVE MATRIX ** C ** (UMAT JACOBIAN) FOR THE CAP MATERIAL UNDER PLANE STRAIN CONDITION ** C *********************************************************************************** C INCLUDE 'ABA_PARAM.INC' C DOUBLE PRECISION EPCM(4,4),SJ2,H,HR,HZ,HT,HRZ DOUBLE PRECISION DSIGR,DSIGZ,DSIGT,DSIGRZ,BMOD,SMOD DOUBLE PRECISION SIGR,SIGZ,SIGT,SIGRZ,EL,TEVP,AW DOUBLE PRECISION CAPL,AH,BH,DFDX,R,AD,ALPHA
166
C C C
COMPUTE STRESS INVARIANTS AND DEVIATORIC STRESSES
1
C C C
PRESS =(SIGR+SIGZ+SIGT)/3. SJ2 =((SIGR-SIGZ)**2+(SIGZ-SIGT)**2+(SIGT-SIGR)**2)/6.D0 +SIGRZ**2 DSIGR =SIGR-PRESS DSIGZ =SIGZ-PRESS DSIGT =SIGT-PRESS DSIGRZ=SIGRZ COMPUTE STIFFNESS COEFFICIENTS OF ELLIPTIC CAP MODEL CAPL=AMIN1(0.,EL) AH =2.D0*(3.D0*PRESS-CAPL) BH =R*R DFDX=-2.D0*(XV(EL,R,ALPHA)-CAPL)
C C C
****** COMPUTE H,HX,HY,HZ,HXY FOR CAP MATERIAL ******
1
C C C
H =9.D0*BMOD*AH*AH+4.D0*BH*BH*SJ2*SMOD-3.D0*AH*DFDX /(AD*(AW+TEVP)) HR =3.D0*BMOD*AH+2.D0*BH*DSIGR*SMOD HZ =3.D0*BMOD*AH+2.D0*BH*DSIGZ*SMOD HT =3.D0*BMOD*AH+2.D0*BH*DSIGT*SMOD HRZ=2.D0*BH*DSIGRZ*SMOD COMPUTE COMPONENTS OF A CONSTITUTIVE MATRIX EPCM(1,1)=BMOD+4./3.*SMOD-HR*HR/H EPCM(1,2)=BMOD-2./3.*SMOD-HZ*HZ/H EPCM(1,3)=BMOD-2./3.*SMOD-HT*HT/H EPCM(1,4)=-HR*HRZ/H EPCM(2,1)=EPCM(1,2) EPCM(2,2)=BMOD+4./3.*SMOD-HZ*HZ/H EPCM(2,3)=BMOD-2./3.*SMOD-HT*HT/H EPCM(2,4)=-HZ*HRZ/H EPCM(3,1)=EPCM(1,3) EPCM(3,2)=EPCM(2,3) EPCM(3,3)=BMOD+4./3.*SMOD-HT*HT/H EPCM(3,4)=-HT*HRZ/H EPCM(4,1)=EPCM(1,4) EPCM(4,2)=EPCM(2,4) EPCM(4,3)=EPCM(3,4) EPCM(4,4)=-HRZ*HRZ/H RETURN END
C C C C *********************************************************************************** C ** UTILITY FUNCTIONS ** C *********************************************************************************** C FUNCTION XV(EL,R,ALPHA) C *********************************************************************************** C ** COMPUTES THE VALUE OF SI1 WHERE CAP SURFACE INTERSECTS WITH HYDROSTATIC AXIS ** C *********************************************************************************** DOUBLE PRECISION EL,R,ALPHA XV=EL-R*(-ALPHA*EL) RETURN END C FUNCTION FF(SI1,ALPHA) C *********************************************************************************** C ** COMPUTES THE VALUE OF SQRT(J2) CORRESPONDING TO SI1 ** C *********************************************************************************** DOUBLE PRECISION SI1,ALPHA FF=-ALPHA*SI1
167
RETURN END C FUNCTION EVP(XLL,AD,AW) C *********************************************************************************** C ** COMPUTES THE ACCUMULATED VOLUMETRIC PLASTIC STRAIN ** C *********************************************************************************** DOUBLE PRECISION XLL,AD,AW EVP=AW*(EXP(AD*XLL)-1.) RETURN END C FUNCTION SRJ2C(SI1,XL,CAPL,R) C *********************************************************************************** C ** COMPUTES THE VALUE OF SRJ2 TO SI1 ON THE ELLIPTIC CAP ** C *********************************************************************************** DOUBLE PRECISION SI1,XL,CAPL,R SRJ2C=SQRT(ABS(XL-CAPL)**2-(SI1-CAPL)**2)/R RETURN END
168
Complete References and Bibliography [1]
ABAQUS User Manual (2005).
[2]
Anandakumar, G. (2005), “Stresses and Deformations in Full-scale Granular Activated Alumina Packed Columns Used for Water Processing Applications in Space”, M.S. Thesis, The University of Connecticut, Storrs, CT.
[3]
Bardet, J.-P. (1997), Experimental Soil Mechanics, Prentice-Hall.
[4]
Bathe, K.-J. (1996), Finite Element Procedures in Engineering Analysis, Prentice-Hall.
[5]
Belytschko, T., Liu, W. K., Moran, B. (2000), Nonlinear Finite Elements for Continua and Structures, Wiley.
[6]
Brown, R. L., Richards, J. C. (1966), Principles of Powder Mechanics, Pergamon.
[7]
Cantelaube F., Goddard, J. D. (1997), “Elastoplastic Arching in a 2-D Granular Heap”, Proc. Powders & Grains ’97, Behringer, Jenkins (Eds.), Balkema, 211.
[8]
Chen, W. F., Baladi, G. Y. (1985), Soil Plasticity, Elsevier.
[9]
Chen, W. F., Mizuno, E. (1990), Nonlinear Analysis in Soil Mechanics, Elsevier.
[10] Claudin, Ph. (1999), “La Physique de Tas de Sable“, Annales de Physique 24, 1-207. [11] Cook, R. D., Young, W. C. (1999), Advanced Mechanics of Materials, Prentice-Hall. [12] Coulomb C. A. (1776), “Sur Une Application des Régles de Maximis & Minimis à Quelques Problèmes de Statique Relatifs à l'Architecture“, Mémoires de Mathématique de l'Académie Royale des Sciences, 343-382. [13] Cowin, S. C. (1977), “The Theory of Static Loads in Bins”, J Appl Mech 44, 409. [14] Cundall P. A., Strack O. D. L. (1979), “A Discrete Numerical Model for Granular Assemblies”, Géotechniqe 29, 47. [15] Dantu, P. (1957), “Contribution à l’Étude Mécanique et Géométrique des Milieux Pulvérulents“, Proc. 4th Int. Conf. on Soil Mech. and Foundation Eng., 133. [16] Das, B. M. (2000), Fundamentals of Geotechnical Engineering, Brooks/Cole. [17] DiMaggio, F. L., Sandler, I. S. (1971), “Material Models for Granular Soils”, J Eng Mech Div ASCE 97, 935-950. [18] Drescher, A. (1991), Analytical Methods in Bin Load Analysis, Elsevier.
169
[19] Drucker, D. C. (1951), “A More Fundamental Approach to Plastic Stress-Strain Relations”, Proc. First U.S. National Congress of Applied Mechanics, 487-491. [20] Drucker, D. C., Prager, W. (1952), “Soil Mechanics and Plastic Analysis or Limit Design”, Q Appl Math 10, 157-164. [21] Drucker, D. C., Gibson, R. E., Henkel, D. J. (1957), “Soil Mechanics and Work-hardening Theories of Plasticity”, Trans ASCE 122, 338-346. [22] Dunne, F., Petrinic, N. (2005), Introduction to Computational Plasticity, Oxford University Press. [23] Duran, J. (1999), Sands, Powders, and Grains: An introduction to the Physics of Granular Materials, Springer. [24] Fung, Y. C. (1965), Foundations of Solid Mechanics, Prentice-Hall. [25] Goddard, J. D. (1998), “Continuum Modeling in Granular Assemblies”, in Physics of Dry Granular Media, Herrmann, H. J. et al. (Eds.), NATO ASI 1997, 1-24. [26] Gopal, J. (2001), “Experimental Load-Deflection Studies of Water De-ionizing Beds”, M.S. Thesis, The University of Connecticut. [27] Gudehus, G. (1985), “Requirements of Constitutive Relations for Soils”, in IUTAM William Prager Symposium Vol. on Mechanics of Geomaterials, Bazánt, Z. P. (Ed.), Wiley. [28] Horne, R. M., Nedderman, R. M. (1976), “Analysis of the Stress Distribution in TwoDimensional Bins by the Method of Characteristics”, Powder Technology 14, 93. [29] Hill, R. (1950), The Mathematical Theory of Plasticity, Clarendon Press. [30] Janssen, H. A. (1895), “Versuche über Getreidedruck in Silozellen“, Zeitung des Vereins Deutscher Ingenieure 39, 1045-1049. [31] Jenike, A.W. (1970), “Storage and Flow of Solids”, Bull. No. 123, Engng. Exp. Station, Univ. Utah. [32] Jerri, A. J. (1992), Integral and Discrete Transforms with Application and Error Analysis, Dekker. [33] Johnson, K. L. (1985), Contact Mechanics, Cambridge University Press. [34] Kolb, E., Mazozi, T., Clément, E., Duran, J. (1999), “Force Fluctuations in a vertically pushed column”, Eur. Phys. J. B 8, 482-493.
170
[35] Kreyszig, E. (1999), Advanced Engineering Mathematics, Wiley. [36] Little, R. W. (1973), Elasticity, Prentice-Hall. [37] Little, R. W., Childs, S. B. (1966), “Elastostatic Boundary Problems in Solid Cylinders”, Quart Appl Math, 15-3, 261-274. [38] Liu C.-H., Nagel S. R., Coppersmith S. N. et al. (1995), “Force Fluctuations in Bead Packs”, Science 269, 513. [39] Love, A. E. H. (1944), A Treatise on the Mathematical Theory of Elasticity, Dover Publications. [40] Lubarda, V. A. (2002), Elastoplasticity Theory, CRC Press. [41] Malla, R. and Gopal, J., (2002a), "Experimental and Analytical Studies of Full-Scale Amberlite Water De-Ionizing Bed for Space Applications," Proc. Of ASCE Engineering Mechanics Division Conference 2002, ASCE, (CD ROM). [42] Malla, R. and Gopal, J., (2002b) “Load and Deflection Characteristics of Water De-Ionizing Medium for Space Applications,” Proc. of SPACE & ROBOTICS 2002, ASCE, p. 356. [43] Malla, R. and Anandakumar, A., (2004), "Force-Displacement Behavior of Activated Alumina Packed Beds for Water Processing Application in Space Life Support Systems," Earth & Space 2004, Malla, Maji (Eds.), ASCE, pp 413-420. [44] Malla, R., and Anandakumar, G., (2006), “Determination of Stress and Deformation Variation in a Cylindrical Bed of Granular Materials with Applications in Space”, Earth & Space 2006, Malla, Binienda, Maji (Eds.), ASCE (CD ROM). [45] Meleshko, V. V. (2003), “Selected Topics in the History of the Two-dimensional Biharmonic Problem”, Appl Mech Rev 56/1, 33-85. [46] Moghe, S. R., Neff, H. F. (1971), “Elastic Deformations of Constrained Cylinders”, J Appl Mech 38, 393-399. [47] Nedderman, R. M. (1992), Statics and Kinematics of Granular Materials, Cambridge University Press. [48] Omega® Engineering, Inc. (2006), “Technical Reference Website on Flowmeters”, (May 16, 2006).
171
[49] Orvalez, G., Clément, E. (2005), “Elastic Medium Confined in a Cylinder versus the Janssen Experiment”, Eur. Phys. J. E 16, 421-438. [50] Pickett, G. (1944), “Applications of the Fourier Method to the Solution of Certain Boundary Problems”, J Appl Mech 11, 176-182. [51] Prevost, J. H. (1987), “Modeling the Behavior of Geomaterials”, in Geotechnical Modeling and Applications, Sayed, M. (Ed.), Gulf Publishing. [52] Potts, D. M., Zdravkovic, L. (1999), Finite Element Analysis in Geotechnical Engineering, Thomas Telford. [53] Radjai F., Wolf D., Jean M., Roux S., Moreau J.-J. (1997), “Force Network in Dense Granular Media”, Proc. Powders & Grains ’97, Behringer, Jenkins (Eds.), Balkema, 211-214. [54] Rankine, W. J. M. (1857), “On the Stability of Loose Earth”, Phil. Trans. Royal Soc. 147. [55] Saada, A. S., (1974), Elasticity – Theory and Application, Pergamon Press. [56] Sandler, I. S., DiMaggio, F. L., Baladi, G. Y. (1976), “Generalized Cap Model for Geological Materials”, ASCE J Geotechnical Division 7, 683-699. [57] Sakaguchi H, Ozaki, E. (1993), “Analysis of the Formation of Arches Plugging the Flow of Granular Materials”, Proc. Powders & Grains ’93, Thornton (Ed.), Balkema. [58] Savage S.B. (1997), “Problems in the Statics and Dynamics of Granular Materials”, Proc. Powders & Grains ’97, Behringer, Jenkins (Eds.), Balkema, 185-194. [59] Savage, S. B. (1998), “Modeling and Boundary Value Problems”, in Physics of Dry Granular Media, Herrmann, H. J. et al. (Eds.), NATO ASI 1997, 25-94. [60] Schofield, A., Wroth, P. (1968), Critical State Soil Mechanics, McGraw-Hill. [61] Sneddon, I. N. (1972), The use of integral transforms, McGraw-Hill. [62] Simo, J. C., Hughes, T. J. R. (1987), “General Return Mapping Algorithms for RateIndependent Plasticity”, in Proc. 2nd Int. Conf. on Constitutive Laws for Engineering Materials Vol. 1, Desai, C. S., Krempl, E. (Eds.), Elsevier. [63] Terzaghi, K. (1943), Theoretical Soil Mechanics, Wiley. [64] Tekscan Inc. (2006), “FlexiForce® User Manual”, (June 12, 2006) [65] Timoshenko, S. P., Goodier, J. N., (1970), Theory of Elasticity, McGraw-Hill.
172
[66] Travers T., Ammi M., Bideau D., Gervois A., Messager J.-C., Troadec J.-P. (1987), “Uniaxial Compression of 2D Packings of Cylinders. Effects of Weak Disorder”, Europhys. Lett. 4, 329. [67] Vanel, L., Clément, E. (1999), “Pressure screening and fluctuations at the bottom of a granular column”, Eur. Phys. J. B 11, 525-533. [68] Walker, D. M. (1966), “An Approximate Theory for Pressures and Arching in Hoppers”, Chem. Eng. Sci. 21, 975-997. [69] Wriggers, P. (2002), Computational Contact Mechanics, Wiley. [70] Zienkiewicz, O. C. (1978), The Finite Element Method, McGraw-Hill.
173
