Geotech Book (2)

Analysis of Geotechnical Problems with Abaqus

2

Overview Day 1 • Lecture 1

Introduction

• Lecture 2

Physical Testing

• Lecture 3

Constitutive Models: Part 1

• Lecture 4

Constitutive Models: Part 2

• Workshop 1

Material Models for Geotechnical Applications

Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

3

Overview Day 2 • Lecture 5

Analysis of Porous Media

• Workshop 2 • Lecture 6

Pore Fluid Flow Analysis: Consolidation

Modeling Aspects

• Workshop 3

Pore Fluid Flow Analysis: Wicking


4

Overview Additional Material • Appendix 1

Stress Equilibrium and Fluid Continuity Equations

• Appendix 2

Bibliography of Geotechnical Example Problems


5

Legal Notices The Abaqus Software described in this documentation is available only under license from Dassault Systèmes and its subsidiary and may be used or reproduced only in accordance with the terms of such license. This documentation and the software described in this documentation are subject to change without prior notice. Dassault Systèmes and its subsidiaries shall not be responsible for the consequences of any errors or omissions that may appear in this documentation. No part of this documentation may be reproduced or distributed in any form without prior written permission of Dassault Systèmes or its subsidiary. © Dassault Systèmes, 2008. Printed in the United States of America Abaqus, the 3DS logo, SIMULIA and CATIA are trademarks or registered trademarks of Dassault Systèmes or its subsidiaries in the US and/or other countries. Other company, product, and service names may be trademarks or service marks of their respective owners. For additional information concerning trademarks, copyrights, and licenses, see the Legal Notices in the Abaqus Version 6.8 Release Notes and the notices at: http://www.simulia.com/products/products_legal.html.


6

Revision Status Lecture 1

8/08

Updated for 6.8

Lecture 2

8/08

Updated for 6.8

Lecture 3

8/08

Updated for 6.8

Lecture 4

8/08

Updated for 6.8

Lecture 5

8/08

Updated for 6.8

Lecture 6

8/08

Updated for 6.8

Appendix 1

8/08

Updated for 6.8

Appendix 2

8/08

Updated for 6.8

Workshop 1

8/08

Updated for 6.8

Workshop 2

8/08

Updated for 6.8

Workshop 3

8/08

Updated for 6.8


Introduction Lecture 1

© Dassault Systèmes, 2008

L1.2

Overview • Introduction • Overview of Geotechnical Applications • Classical and Modern Design Approaches • Some Cases for Numerical (FE) Analysis • Experimental Testing and Numerical Analysis • Requirements for Realistic Constitutive Theories


Introduction


L1.4

Introduction • In this lecture we present an overview of geotechnical applications and discuss the philosophy on which the usage of numerical (finite element) analysis for geotechnical problems is based. • Lecture 2 deals with experimental testing and how it relates to the calibration of constitutive models for geotechnical materials. • The different Abaqus constitutive models applicable to geotechnical materials are presented in Lectures 3 and 4. • Their usage, calibration, implementation, and limitations are discussed. • In Lecture 5 we outline the treatment of porous media in Abaqus and discuss the coupling between fluid flow and stress/deformation. • Several modeling issues relating to geotechnical situations are discussed in Lecture 6.


Overview of Geotechnical Applications


L1.6

Overview of Geotechnical Applications • Some typical geotechnical applications


L1.7

Overview of Geotechnical Applications • Main characteristics of soil materials • Soil is composed of solid grains and voids; the voids may be fully or partially filled with pore fluid. • When the soil is loaded, the pressure in the pore fluid gets modified. The pore fluid then flows and the soil deforms with time. • The deformations and pore fluid flow will be governed by the mechanical and permeability properties of the soil, the mechanical and pore pressure boundary conditions, contact conditions, presence of reinforcements, etc.


L1.8

Overview of Geotechnical Applications • The following need to be considered before modeling geotechnical problems: Is coupled pore fluid diffusion – stress analysis necessary?

In the coupled setting the pore fluid flow affects the deformations and the deformations affect the pore fluid flow In an uncoupled setting the pore fluid flow does not affect the deformations, and the problem can then be modeled as a deformation problem or just a pore fluid flow problem.

Is contact modeling important for the particular application?

In some cases detailed modeling of contact conditions is not necessary.

Which material model will be appropriate for representing the physical material behavior?

Select appropriate material model according to the physical material behavior.


L1.9

Overview of Geotechnical Applications • The following need to be considered before modeling geotechnical problems: Does the application contain reinforcement?

Use embedded elements, if necessary, to model reinforcements. For example, for modeling of ties, soil nails, etc.

Is modeling of excavation and construction important?

Use model change, if necessary, in setting up the analysis steps. Using model change one can activate or deactivate elements and contact pairs and can model construction and excavation.

How detailed is the information one has about the initial conditions in the model?

For ascertaining equilibrium at time=0, one needs to begin with a Geostatic step, wherein known initial conditions, boundary conditions, and external loads are applied.


L1.10

Overview of Geotechnical Applications • Dams

Saturation values inside a dam

• For the modeling of dams one needs to take into account: • Fully and partially saturated flow • Appropriate material models for the dam material • Mechanical behavior • Pore fluid flow behavior Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L1.11

Overview of Geotechnical Applications • The following features also need to be taken into account for the modeling of dams: • Sequential construction—model change • Seepage analysis—transient consolidation • Earthquake loading—dynamics


L1.12

Overview of Geotechnical Applications • Modeling of jointed rocks • Response depends on joint orientations and asperities of adjacent surfaces • Joints can open and close, which can produce induce anisotropy at macro scale • Need to take into account these features in a continuum setting


http://www.outreach.canterbury.ac.nz/resources /geology/glossary/joints.jpg inclined joint set

vertical joint set

L1.13

Overview of Geotechnical Applications • Modeling of pile foundations • Contact modeling is important • Modeling of friction between the pile and the soil • Cohesive elements with damage characteristics • Soil plasticity • Coupled consolidation analysis for modeling short and long term behavior Distribution of pore pressure immediately after a pile is driven


L1.14

Overview of Geotechnical Applications • Slope failure analysis • Unreinforced and reinforced slopes • Soil plasticity • Material failure • Embedded elements

Stability analysis using Abaqus/Explicit Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Stability analysis of unreinforced and reinforced slope at working load levels using Abaqus/Standard

L1.15

Overview of Geotechnical Applications • Retaining walls • Contact • Soil plasticity models • Coupled consolidation • Embedded elements http://www.versa-loc.biz/versatech/images/wallcomp_soils.jpg

• Landfills • Quasi-static analysis • Creep modeling • Geomembranes and geo-synthetics • Embedded elements • Contact

http://www.keepingavclean.com/images/landfill105.jpg


L1.16

Overview of Geotechnical Applications • Mining and tunneling • Excavations and construction

Shotcrete liners

Top heading

layer 4 layer 3

• Liners and ties • Material needs to be added or removed • Contact conditions need to be taken into account


Invert and Concrete floor

Gravity


layer 1

Gravity

• Transient material degradation needs to be represented in the model

layer 2

Bench

L1.17

Overview of Geotechnical Applications • Nuclear waste disposal • Sequentially coupled thermo-hydromechanical (THM) analysis • Ground water seepage modeling • Cross fluid flow, heat flow, and mechanical contact interactions between soil and non-soil regions • Examine worst-case scenarios, provide data on environmental impact Bentonite blocks

Independent thermal expansion properties for: • Soil matrix • Pore fluid

Steel liner Heater (diameter 0.9) Granite Granite Bentonite barrier

Concrete plug

Service zone, control and data acquisition system Principal access tunnel to KWO

Heaters

2.28

Void ratio-dependent permeability

4.54

1.00

4.54

17.4

4.34 2.7 70.4


(Dimensions in meters)

Courtesy: SIMULIA/Scandinavia


L1.18

Overview of Geotechnical Applications • Abaqus provides the following features for geotechnical applications: • Full coupling between stress and pore fluid diffusion • Anisotropic and void ratio-dependent permeability, velocity-dependent permeability • Wide range of soil plasticity models including dilation and compaction available in Abaqus/Standard and Abaqus/Explicit • Jointed material model for modeling jointed rocks • Saturated and unsaturated media and automatic computation of phreatic surface and capillary zone


L1.19

Overview of Geotechnical Applications • Abaqus features for geotechnical applications (contd.): • Contact between permeable/permeable and permeable/impermeable regions • Stress equilibrium as well as fluid continuity maintained • Embedded elements for modeling reinforcement • Elements can be added or deleted for modeling construction and excavation • Thermal expansion for pore fluid and solid matrix • Independent bulk moduli for soil grains, pore fluid, and elastic properties for soil matrix • Moisture and gel swelling


Classical and Modern Design Approaches


L1.21

Classical and Modern Design Approaches • In the classical approach two basic types of calculations are done: failure estimates and deformation estimates. • Failure estimates are based on rigid perfectly plastic stress-strain assumptions:

• Examples:


L1.22

Classical and Modern Design Approaches • Examples (cont.):

Bearing capacity of a foundation

Retaining wall stability

• The result is a factor of safety, which is evaluated based on experience (design code). Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L1.23

Classical and Modern Design Approaches • Deformation estimates assume linear elastic behavior with average elastic properties:

• Foundation settlement example: settlement

p is bearing pressure

⎛ 1 − v2 ⎞ w = pb ⎜ ⎟ f, ⎝ E ⎠

b is width of foundation E, v are average elastic properties f is shape factor (based on small scale tests)


L1.24

Classical and Modern Design Approaches • In the modern approach failure and deformation characteristics are obtained from the same analysis. The analysis requires: • a complete constitutive model, which can include the loading and unloading behavior, and

• the numerical solution of a boundary value problem. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L1.25

Classical and Modern Design Approaches

• Numerical (finite element) analysis can handle arbitrary geometries. • When finite elements are used, we have the option of allowing different shear strains along the failure surface, and thus the stress state along the failure surface can change from point to point. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Some Cases for Numerical (FE) Analysis


L1.27

Some Cases for Numerical (FE) Analysis • Slope stability and deformation:

• Classical limit failure calculations can predict ultimate stability of a slope. • Numerical finite element analysis is necessary for the calculation of deformations and detailed soil behavior.


L1.28

Some Cases for Numerical (FE) Analysis • Construction of an earth dam and subsequent filling of reservoir: • Detailed soil behavior plays an important role.

• Local soil failure may trigger overall collapse or hamper functionality. • The event sequence - construction, filling of reservoir, and long-term consolidation - must be considered using numerical modeling. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L1.29

Some Cases for Numerical (FE) Analysis • Tunneling: • The initial state of stress of the soil or rock mass is important.

• Stability depends on: • Initial state of stress • Sequence of excavation • Stabilizing aids such as liners and rock bolts Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Experimental Testing and Numerical Analysis


L1.31

Experimental Testing and Numerical Analysis Laboratory testing

Constitutive model

Calibration

Finite element model

Small or large scale testing

Comparison of model predictions with test results

Design prototype Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L1.32

Experimental Testing and Numerical Analysis • Experimental testing and constitutive model calibration: • Laboratory tests should represent the load or deformation conditions of the physical problem. • Choose a constitutive model that represents the major characteristics; minor features may be ignored. • Obtain test measurements suitable for the specific constitutive model. • Calibrate the model parameters using measured values of repeatable test results.


L1.33

Experimental Testing and Numerical Analysis

• Finite Element Analysis: • Numerical model should capture the important features of the physical problem without irrelevant detail. • Use of an appropriate constitutive model is critical, although simplifications are often justifiable. • Provides a tool to aid design, supplementing engineering judgment and experience.


Requirements for Realistic Constitutive Theories


L1.35

Requirements for Realistic Constitutive Theories • Realistic constitutive models should be: • able to represent the macro behavior as dictated by micromechanics • suitable for numerical implementation • able to represent material behavior in any relevant spatial situation: • one-dimensional • plane strain • axisymmetric • three-dimensional • able to reasonably extrapolate to conditions that cannot be reproduced with laboratory testing equipment


Physical Testing Lecture 2


L2.2

Overview • Physical Testing • Basic Experimental Observations • Testing Requirements and Calibration of Constitutive Models


Physical Testing


L2.4

Physical Testing • Geotechnical materials: • are generally voided and sensitive to volume changes. • Volume changes are closely tied to the magnitude of the hydrostatic pressure stress: • Hence, it is important to test the materials over the range of hydrostatic pressure of interest. • Standard tests: • Hydrostatic (or isotropic) compression tests • Oedometer (or uniaxial strain) tests • Triaxial compression and extension tests • Uniaxial compression tests (a special case of triaxial compression) • Shear tests


L2.5

Physical Testing • A practical constitutive model should require for calibration only the information generated by these standard tests. • Assumptions need to be made regarding tensile and true triaxial behavior as direct tensile tests and true triaxial tests are difficult to perform. • The diversity of geotechnical materials means that a wide range of behaviors is possible. • What follows are some very general observations for frictional materials.


Basic Experimental Observations


L2.7

Basic Experimental Observations • Hydrostatic or isotropic compression test:


L2.8

Basic Experimental Observations • Oedometer or uniaxial strain test:


L2.9

Basic Experimental Observations • Triaxial compression tests:


L2.10

Basic Experimental Observations • The “critical state” concept:

Void ratio e

• Casagrande defined “critical state” as the state (for monotonic loading) at which continued shear deformation can occur without further change in effective stress and volume (void ratio) of the material. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L2.11

Basic Experimental Observations • Cyclic tests: • Isotropic compression or uniaxial strain (oedometer) tests

• Triaxial or uniaxial compression tests provide the deviatoric behavior


L2.12

Basic Experimental Observations • Truly triaxial (cubical) tests:


L2.13

Basic Experimental Observations • Essential aspects of behavior of voided frictional materials: + + + + + +

Nonlinear stress-strain behavior Irreversible deformations Influence of hydrostatic pressure stress on “strength” Influence of hydrostatic pressure on stress-strain behavior Influence of intermediate principal stress on “strength” Shear stressing-dilatancy coupling


L2.14

Basic Experimental Observations • Essential aspects of behavior of voided frictional materials (contd.): + + + +

Influence of hydrostatic pressure stress on volume changes Hardening/softening related to volume changes Stress path dependency Effects of small stress reversals

– Effects of large stress reversals (hysteresis) – Degradation of elastic stiffness after large stress reversals + Included in Abaqus – Not included in Abaqus


Testing Requirements and Calibration of Constitutive Models


L2.16

Testing Requirements and Calibration of Constitutive Models • Basic requirements for laboratory testing of geotechnical materials: • The specimens must be tested under the assumption that they represent an average material behavior: • In the ground there will be some variation within the same soil or rock mass, and the spatial scale of such variations may be large compared to laboratory test specimens. • The tests must simulate the in situ conditions as closely as possible: • drainage conditions • density of the material • range of stresses • For deep mining cases these may be difficult to obtain.

• All stresses and strains must be measured throughout the stressstrain response: • to allow complete characterization of the constitutive behavior. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L2.17

Testing Requirements and Calibration of Constitutive Models • Laboratory testing should be guided by a previously proposed constitutive model: • The understanding of this model is necessary for the correct interpretation of laboratory tests. • The model parameters should be physical and measurable in practicable experiments.


L2.18

Testing Requirements and Calibration of Constitutive Models • Laboratory tests for calibrating relevant components of constitutive models: • Isotropic compression test and uniaxial strain (oedometer) test: • One test is required to calibrate hydrostatic behavior. • One unloading is necessary to calibrate the elastic part of this behavior.

Isotropic compression

• Triaxial compression tests: • Two (preferably more) tests are required to calibrate the shear behavior and its hydrostatic pressure dependence. • One unloading in each test is necessary to calibrate the elastic part of this behavior. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Triaxial compression

L2.19

Testing Requirements and Calibration of Constitutive Models • Triaxial extension tests: • The triaxial extension test is performed in a triaxial machine. • In the presence of confining lateral pressures the vertical stress is gradually decreased till failure.

• Two (preferably more) such tests are required to calibrate the intermediate principal stress dependence of the shear behavior. • Direct tension test: • The direct tension test or the uniaxial tension test is performed by applying uniaxial strain along a particular direction. • One test is required to calibrate the tensile behavior of cohesive materials like rocks or soils with cohesion.


L2.20

Testing Requirements and Calibration of Constitutive Models

• Truly triaxial or cubical tests: • Many tests are required to calibrate the behavior of the material when subjected to different stresses in all directions.


L2.21

Testing Requirements and Calibration of Constitutive Models • Shear box tests and indirect tension (Brazilian) test: • These are useful in calibrating the cohesive properties of materials.

Shear box test

Indirect tension test

• Multiple unloading-reloading cycles in any of above tests: • These are useful for calibrating the effects of large stress reversals such as hysteresis and elastic degradation. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Constitutive Models: Part 1 Lecture 3


L3.2

Overview • Stress Invariants and Spaces • Overview of Constitutive Models • Elasticity • Plastic Behavior of Soils • Mohr-Coulomb Model • Extended Drucker-Prager Models


Stress Invariants and Spaces


L3.4

Stress Invariants and Spaces • Stress in three dimensions: three direct and three shear components.

• Symmetry of stress tensor:

σ 12 = σ 21, σ 13 = σ 31, σ 23 = σ 32 . • Abaqus convention: tensile stress is positive. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.5

Stress Invariants and Spaces • Principal stresses: stresses normal to planes in which shear stresses are zero.


L3.6

Stress Invariants and Spaces • In two dimensions (Mohr’s circle):

σ 1,2 =

σ 11 + σ 22


⎛ σ − σ 22 ⎞ 2 ± ⎜ 11 ⎟ + τ 12 2 2 ⎝ ⎠ 2τ 12 tan 2θ = σ 11 − σ 22 2

L3.7

Stress Invariants and Spaces • Stress decomposition: deviatoric plus hydrostatic:

σ = S − pI . • Abaqus invariants: pressure stress, p = −

1 trace(σ ), 3

Mises equivalent stress, q =

3 ( S : S ), 2 1

⎛9 ⎞3 third invariant, r = ⎜ S ⋅ S : S ⎟ . 2 ⎝ ⎠


L3.8

Stress Invariants and Spaces • For all models except the Mohr-Coulomb, also define the deviatoric stress measure:

t=

q⎡ 1 ⎛ 1 ⎞⎛ r ⎞ ⎢1 + − ⎜ 1 − ⎟⎜ ⎟ 2 ⎢ K ⎝ K ⎠⎝ q ⎠ ⎣

3⎤

⎥, ⎥⎦

so that t = q/K in triaxial tension (r = q) and t = q in triaxial compression (r = − q). If K = 1, t = q. K is typically between 0.8 and 1.0. • t = constant is a “rounded” surface in the deviatoric plane.


L3.9

Stress Invariants and Spaces • Useful planes:


L3.10

Stress Invariants and Spaces • Meridional plane:


L3.11

Stress Invariants and Spaces • Deviatoric (or Π ) plane:


Overview of Constitutive Models


L3.13

Overview of Constitutive Models • Elasticity models: • Linear, isotropic • Porous, isotropic (nonlinear) • Damaged, orthotropic (nonlinear; used in concrete, jointed material) • Plasticity models: • Open surface, pressure independent (Mises) • Open surface, pressure dependent (Drucker-Prager, Mohr-Coulomb) • Closed surface (Cam-clay, Drucker-Prager with Cap) • Multisurface (jointed material) • Nested surfaces (bounding surface*)


L3.14

Overview of Constitutive Models • Other inelastic models: • Continuum damage theories* • Endochronic theories*

None of the available models (with the possible exception of the jointed material model) is capable of accurately handling large stress reversals such as those occurring during cyclic loading or severe dynamic events.

*Not available in Abaqus for geotechnical materials Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Elasticity


L3.16

Elasticity • Linear isotropic elasticity:

Input file usage: *ELASTIC E,ν


Young’s modulus and Poisson’s ratio

L3.17

Elasticity • Porous elasticity: • Nonlinear, isotropic • Pressure stress varies as an exponential function of volumetric strain:

(

( ))

⎡1 + e0 ⎤ el p = − ptel + p0 + ptel exp ⎢ 1 − exp ε vol ⎥ ⎣ κ ⎦ ⎛ p0 + ptel ⎞ κ el ln ⎜ or ⎟ = J − 1, (1 + e0 ) ⎝ p + ptel ⎠

(

)

where Jel − 1 is the nominal volumetric elastic strain.

J = dV/dV ° is the ratio of current volume to reference volume. el = ln(Jel) ε vol

and Jel = exp (ε vol). el

el is the logarithmic measure of volumetric elastic strain. ε vol

ptel is the elastic tensile stress limit, can be zero or non zero. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.18

Elasticity • Deviatoric behavior: • Constant shear modulus S = 2Ge , (deviatoric elastic stiffness independent of pressure stress) el

ˆ e , • Constant Poisson’s ratio dS = 2Gd (deviatoric stiffness dependent on the pressure stress) el

where the instantaneous shear modulus, Gˆ , is

3 (1 − 2v ) (1 + e0 ) el Gˆ = p + ptel exp ε vol . 2 (1 + v ) κ

(

) ( )

κ , ptel , G, v are material parameters, p0 is the

initial value of hydrostatic pressure stress, and e0 is the initial void ratio.


L3.19

Elasticity

Porous elasticity el el • ε vol has an arbitrary origin, defined so that p = p0 at ε vol = 0.


L3.20

Elasticity • Usage *POROUS ELASTIC, SHEAR=POISSON κ, ν, ptel *POROUS ELASTIC, SHEAR=G κ, G, ptel


L3.21

Elasticity • ∗INITIAL CONDITIONS, TYPE=RATIO is required to define the initial void ratio (porosity) of the material when porous elasticity is used. • User subroutine VOIDRI can be used to specify complex initial void ratio distributions.


Plastic Behavior of Soils


L3.23

Plastic Behavior of Soils • Plasticity models are characterized by the following: • Yield surface • It is a surface defining the yield criterion for the plasticity model. • Flow rule • The condition determining the plastic strain after yielding. • Plastic flow can be associated or nonassociated. • Associated flow means that the plastic strain direction is along the direction of yielding, or normal to the yield surface. • Nonassociated flow means that the plastic strain direction is not along the direction of yielding.


L3.24

Plastic Behavior of Soils • Dilation • Dilation during plastic flow means that the material changes its volume while undergoing yielding.

Volume increases

Yielding without dilation (ductile metals) Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Yielding with dilation (soils)

L3.25

Plastic Behavior of Soils • Yielding in soils is associated with internal friction between soil grains. • Internal friction between soil grains results from interference at the grain level. Direction of yielding Friction angle φ

u

Dilation angle ψ

v Volume increases

• Dilation angle ψ = arctan(v/u). • The plastic flow will be associative when ψ = φ. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.26

Plastic Behavior of Soils • Yielding in ductile metals is generally not accompanied by volume change. • The yielding is due to dislocation. • The flow is associative as there is no internal friction and the material does not dilate during yielding. • Both φ and ψ are zero. • Yield stress does not depend on pressure. Yield surface Shear stress

Pressure


L3.27

Plastic Behavior of Soils • Yielding in soils is generally accompanied by volume change. • Existence of internal friction (nonzero φ) • Existence of cohesion • similar to resistance against dislocation in ductile metals • The material may dilate during yielding in shear. • The flow can be associative (ψ = φ ), or nonassociative. • Yield stress depends on pressure.

Yield surface Shear stress

Pressure Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Mohr-Coulomb Model


L3.29

Mohr-Coulomb Model • The characteristics of the Mohr-Coulomb plasticity model in Abaqus: • The model is intended for granular materials like soils under monotonic loading. • It does not consider rate dependence. • The linear isotropic elastic response is followed by non-recoverable response idealized as being plastic. • The yield behavior depends on the hydrostatic pressure: • The material becomes stronger as the confining pressure increases. • The yield behavior may be influenced by the magnitude of the intermediate principal stress.


L3.30

Mohr-Coulomb Model

• The model includes isotropic hardening or softening. • The inelastic behavior is generally accompanied by volume change. • The flow rule may include: • inelastic dilation as well as • inelastic shearing • The plastic flow potential is smooth and nonassociated. • Material properties can be temperature dependent.


L3.31

Mohr-Coulomb Model • Description of the Mohr-Coulomb Model: • The model requires linear isotropic elasticity. • The Mohr-Coulomb yield function is given by:

F = Rmc q − p tan φ − c = 0, where

Rmc(Θ,φ) is a measure of the shape of the yield surface in the deviatoric

plane,

Rmc =

1 π⎞ 1 ⎛ π⎞ ⎛ sin ⎜ Θ + ⎟ + cos ⎜ Θ + ⎟ tan φ , 3 3 3⎠ 3 cos φ ⎝ ⎠ ⎝

φ is the slope of the Mohr-Coulomb yield surface in the Rmcq − p stress

plane, which is commonly referred to as the friction angle of the material,

0 ≤ φ < 90;


L3.32

Mohr-Coulomb Model c is the cohesion of the material; and Θ is the deviatoric polar angle defined as cos(3Θ) =

r3 . q3

• The Mohr-Coulomb model assumes that the hardening is defined in terms of the material’s cohesion, c. • The cohesion can be defined as a function of plastic strain, temperature, or field variables. • The hardening is isotropic.


L3.33

Mohr-Coulomb Model

(a)

(b)

Yield surface in the meridionial plane (a) and the deviatoric plane (b) Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.34

Mohr-Coulomb Model • The flow potential, G, is chosen as a hyperbolic function in the meridional stress plane and the smooth elliptic function proposed by Menétrey and Willam (1995) in the deviatoric stress plane:

G=

( ε c |0 tanψ )2 + ( Rmwq )2 − p tanψ .

(

)

• The initial cohesion of the material, c |0 = c ε pl = 0 ; the dilation angle, ψ ; and the meridional eccentricity, ε, control the shape of G in the meridional plane.


L3.35

Mohr-Coulomb Model

ε defines the rate at which G approaches the asymptote. • The flow potential tends to a straight line in the meridional stress plane as the meridional eccentricity tends to zero.

Mohr-Coulomb flow potential in the meridional plane


L3.36

Mohr-Coulomb Model • Rmw(Θ, e,φ) controls the shape of G in the deviatoric plane:

Rmw =

(

)

4 1 − e2 ( cos Θ ) + ( 2e − 1)

(

)

2

(

2 1 − e cos Θ + ( 2e − 1) 4 1 − e 2

2

2

) ( cos Θ )

2

⎛π ⎞ Rmc ⎜ ,φ ⎟ ⎝3 ⎠ + 5e − 4e 2

• The deviatoric eccentricity, e, describes the “out-of-roundedness” of the deviatoric section in terms of the ratio between the shear stress along the extension meridian (Θ = 0) and the shear stress along the compression meridian (Θ = π / 3). • The default value of the deviatoric eccentricity is calculated as e =

3 − sin φ 3 + sin φ

and allows the Abaqus Mohr-Coulomb model to match the behavior of the classical Mohr-Coulomb model in triaxial compression and tension. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.37

Mohr-Coulomb Model

• The deviatoric eccentricity may have the following range: ½ ≤ e ≤ 1.0. • The plastic flow in the deviatoric plane is always nonassociated.


L3.38

Mohr-Coulomb Model • Usage and calibration • Linear, isotropic elasticity must be used. • The ∗MOHR COULOMB option is used

e ε

to define ψ, φ, e, and ε. • The ECCENTRICITY parameter is used to define ε. The default value is 0.1. • The DEVIATORIC ECCENTRICITY parameter

φ

can be used to define e. The deviatoric eccentricity may have the following range: 1 < e ≤ 1.0.

2

• If e is defined directly, Abaqus matches the classical Mohr-Coulomb model only in triaxial compression. • The friction angle, φ, and the dilation angle, ψ, can be functions of temperature and field variables. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

ψ

L3.39

Mohr-Coulomb Model • The *MOHR-COULOMB option must always be accompanied by the *MOHR-COULOMB HARDENING material option: • It defines the evolution of the cohesion, c: • It can be made temperature and/or field dependent. • *EXPANSION can be used to introduce thermal volume change effects.


L3.40

Mohr-Coulomb Model • The plastic flow in the deviatoric plane is always nonassociated; • Needs the unsymmetric solver (∗STEP, UNSYMM=YES).


L3.41

Mohr-Coulomb Model • Calibration of the Mohr-Coulomb model: • Use the critical stress states from several different triaxial tests. • The critical stress states are plotted in the meridional plane to provide an estimate of: • friction angle, φ, and • cohesion, c.

φ Rmc q

c Pressure

p

• The dilation angle, ψ, is chosen so that the volume change during the plastic deformation matches that seen experimentally. • Hardening data is obtained from one of the triaxial tests. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.42

Mohr-Coulomb Model • The Abaqus Mohr-Coulomb model uses a smooth plastic flow potential: • It does not always provide the same plastic behavior as the classical (associated) Mohr-Coulomb model, which has a faceted flow potential. • With the default value of deviatoric eccentricity, e, Abaqus does match classical Mohr-Coulomb behavior under triaxial extension or compression. • Benchmark Problem 1.14.5, Finite deformation of an elastic-plastic granular material, shows how to match the Abaqus Mohr-Coulomb model to the classical Mohr-Coulomb model for plane strain deformation.


L3.43

Mohr-Coulomb Model • Example: Limit Load Analysis of a Strip Foundation • This example involves the limit load analysis of a strip foundation (Benchmark Problem 1.14.4). • A rigid footing 10 feet wide on a granular soil material of 12 feet depth is loaded and the average pressure vs. displacement of the footing is computed. • The model is assumed to be in a plane strain condition. Right half of footing

Soil

Symmetry plane

• CPE8R finite elements and CINPE5R infinite elements are used. • The footing is driven by prescribed displacement. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.44

Mohr-Coulomb Model • Boundary conditions: • Symmetry boundary conditions along the symmetry plane • Fixed at the bottom • Rough friction assumed between the footing and the soil • Soil nodes lying on the footing are constrained to move vertically in unison and with zero horizontal displacements using linear constraint equations. • Mohr-Coulomb model: • Elastic modulus 30,000 psi and Poisson’s ratio 0.3 • Friction angle 20o • Cohesion 10 psi • Dilation angles used: • 0o (nondilatant flow) • 20o (dilatant flow) Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.45

Mohr-Coulomb Model • The Mohr-Coulomb material data is specified as follows:

*MATERIAL,NAME=A1

Elastic modulus and Poisson’s ratio

*ELASTIC 30000.,0.3

Friction angle and dilation angle

*MOHR COULOMB 20.,0. *MOHR COULOMB HARDENING

Cohesion

10.,

• A prescribed displacement of 5 inches is applied and the reaction is computed. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.46

Mohr-Coulomb Model • Results Terzaghi (175 psi)

Prandtl (143 psi)

Mohr Coulomb with no dilation Mohr Coulomb with dilation


L3.47

Mohr-Coulomb Model • Results • The Abaqus results lie close to the Terzaghi solution for strip foundations on granular soils. • The solutions by the Prandtl theory and the Terzaghi theory are from: • Chen, W. F., Limit Analysis and Soil Plasticity, Elsevier, Amsterdam, 1975. • The Prandtl solution considers the stresses due to the weight of the soil as small compared to the soil strength. • This theory was often used for predicting the behavior of metals indented by rigid materials. • The Terzaghi solution takes into account the soil and overburden weight and also the material cohesion.


Extended Drucker-Prager Models


L3.49

Extended Drucker-Prager Models • The characteristics of the Extended Drucker-Prager models: • These models are intended for monotonic loading. • For example, the limit load analysis of a soil foundation. • These models are the simplest available models for simulating frictional materials. • The elastic response in these models is followed by a non-recoverable response, which is idealized as being plastic. • The material in these models is initially isotropic. • The yield behavior of these models depends on the hydrostatic pressure. • The material becomes stronger as the confining pressure increases. • The yield behavior may be influenced by the magnitude of the intermediate principal stress.


L3.50

Extended Drucker-Prager Models • The characteristics of the Extended Drucker-Prager models (contd.): • These models differ in the manner in which the hydrostatic pressure dependence is introduced. • These models include isotropic hardening or softening. • The inelastic behavior in these models is generally accompanied by volume change: • The flow rule may include inelastic dilation as well as inelastic shearing. • These models can incorporate strain-rate dependent material properties. • The material properties in these models can be made temperature dependent. • Either linear elasticity or nonlinear porous elasticity can be used with these models.


L3.51

Extended Drucker-Prager Models • Three Extended Drucker-Prager models based on the yield surface shape in the meridional plane are available in Abaqus: • Linear • Hyperbolic • Exponent • The particular form of the material model can be chosen based on: • the kind of material being modeled • the available experimental data for calibration of the model parameters • the range of pressure stress values that the material is likely to see


L3.52

Extended Drucker-Prager Models • The Linear Drucker-Prager model • The yield surface of the linear model is written as

F = t − p tan β − d = 0.

• p is the equivalent pressure stress • t is the deviatoric stress measure • β is the friction angle

• d is cohesion and is related to the hardening input data • The model allows for separate dilatation and friction angles


L3.53

Extended Drucker-Prager Models • The model provides a non-circular yield surface in the deviatoric plane: • It allows for the matching of different yield values in tension and compression. • It provides flexibility in fitting experimental results. • However, the surface is too smooth to be a close approximation to the Mohr-Coulomb surface.

K is the ratio of the yield stress in triaxial tension to the yield stress in triaxial compression.


L3.54

Extended Drucker-Prager Models • We assume a (possibly) nonassociated flow rule, where the direction of the inelastic deformation vector is normal to a linear plastic potential, G:

d ε pl =

d ε pl ∂G , c ∂σ

where G = t − p tanψ ,

c is a constant that depends on the type of hardening data, dε

pl

dε

pl

= d ε11pl

= d ε11pl d γ pl d ε pl = 3

in uniaxial compression, in uniaxial tension, and in pure shear.

ψ is the dilation angle in the p–t plane. This flow rule definition precludes dilation angles ψ > 71.5° ( tan ψ > 3) , which is not likely to be a limitation for real materials.


L3.55

Extended Drucker-Prager Models • The flow is associated in the deviatoric plane but nonassociated in the p–t plane if ψ ≠ β. • For ψ = 0, the material is nondilatational. • If ψ = β, the model is fully associated.


L3.56

Extended Drucker-Prager Models • Hyperbolic model • The hyperbolic yield criterion is a continuous combination of the maximum tensile stress condition of Rankine (tensile cut-off) and the linear DruckerPrager condition at high confining stress. It is written as

F = l02 + q 2 − p tan β − d ′ = 0, where d' is the hardening parameter that is related to the hardening input data as

d ′ = l02 + σ c2 − d ′ = l02 + σ t2 +

σc 3

σt 3

tan β if hardening is defined by uniaxial compression, σc; tan β if hardening is defined by uniaxial tension, σt;

d ′ = l02 + d 2 if hardening is defined by shear (cohesion), d. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.57

Extended Drucker-Prager Models • l0 = d ′ 0 − pt 0 tan β determines how quickly the hyperbola approaches its asymptote (see sketch). • pt 0 is the initial hydrostatic tension strength of the material, d ′ 0 is the initial value of d', and β is the friction angle measured at high confining pressure. • The model treats β and l0 as constants during hardening.


L3.58


• The yield surface is a von Mises circle in the deviatoric stress plane. (The K parameter is not available for this model.)


L3.59

Extended Drucker-Prager Models • Exponent model • The general exponent form provides the most general yield criterion available in this class of models. The yield function is written as

F = aq b − p − pt = 0, where a and b are material parameters independent of plastic deformation and pt is the hardening parameter that represents the hydrostatic tension strength of the material and is related to the input data as

pt = aσ cb − pt = aσ tb +

σc 3

σt 3

if hardening is defined by uniaxial compression, σc ; if hardening is defined by uniaxial tension, σt ; and

pt = ad b if hardening is defined by shear (cohesion), d.


L3.60


• The yield surface is a von Mises circle in the deviatoric stress plane. (The parameter K is not available for this model.) • The material parameters a, b, and pt : • Can be specified directly, or • Abaqus will determine them from specified triaxial test data using a least squares fit.


L3.61

Extended Drucker-Prager Models • Flow in the hyperbolic and exponent models • Flow potential governing the plastic flow in the hyperbolic and general exponent models:

G=

(ε σ 0 tanψ )

2

+ q 2 − p tanψ ,

where,

ψ is the dilation angle in the meridional plane at high confining pressure, σ 0 is the initial yield stress, ε is the eccentricity defining the rate at which the function approaches its asymptote • The flow potential tends to a straight line as the eccentricity tends to zero.


L3.62

Extended Drucker-Prager Models • The function approaches the linear Drucker-Prager flow potential asymptotically at high confining pressure and intersects the hydrostatic pressure axis at 90°. • The potential is continuous and smooth. • It ensures that the flow direction is always defined uniquely.

• The potential is the von Mises circle in the deviatoric stress plane. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.63

Extended Drucker-Prager Models • Associated flow is obtained in the hyperbolic model if β = ψ and

ε=

l0 . σ 0 tanψ

• In the general exponent model the flow is always nonassociated in the meridional plane. • The default flow potential eccentricity is ε = 0.1. • It provides an almost constant dilatational angle over a wide range of confining pressure stress values. • Increasing the value of ε provides more roundedness to the flow potential. • The dilation angle increases smoothly as the confining pressure decreases.

ε < 0.1 may lead to convergence problems if the material is subjected to low confining pressures because of the very sharp curvature of the flow potential near its intersection with the p-axis.


L3.64

Extended Drucker-Prager Models • Usage • These models are invoked by the ∗DRUCKER PRAGER material option. • The SHEAR CRITERION parameter is set to LINEAR, HYPERBOLIC, or EXPONENT to define the yield surface shape. • The ∗DRUCKER PRAGER HARDENING option should also be used. • This option defines the evolution of the yield stress in uniaxial compression (TYPE=COMPRESSION), in uniaxial tension (TYPE=TENSION), or in pure shear (TYPE=SHEAR). Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.65

Extended Drucker-Prager Models • The yield function can be made rate dependent by using the ∗RATE DEPENDENT option or by specifying the yield stress as a function of the plastic strain rate. • Rate dependency is rarely used for geotechnical materials. • It may however be important for polymers.


L3.66

Extended Drucker-Prager Models • The elasticity is defined by: • linear elasticity, or • porous elasticity. • All material parameters in these models can be made temperature or field dependent. – Thermal expansion can be used to introduce thermal volume change effects. – ∗INITIAL CONDITIONS, TYPE=RATIO is required to define the initial void ratio of the material if porous elasticity is used.


L3.67

Extended Drucker-Prager Models • Analyses using a nonassociated flow version of the model may require the use of the unsymmetric solver because of the resulting unsymmetric plasticity equations. • If the default symmetric solver is used when the flow is nonassociated, Abaqus may not find a converged solution.


L3.68

Extended Drucker-Prager Models • Matching experimental data • The simplest modified Drucker-Prager model requires at least two experiments for calibration. • Simplest model: linear, rate independent, temperature independent, and yielding independent of the third stress invariant. • Most common experiments performed: • uniaxial compression (for cohesive materials) • triaxial compression or tension tests • shear tests for cohesive materials • The uniaxial compression test involves compressing the sample between two rigid platens. • Record the load and displacement in the direction of loading • Record the lateral displacements for measuring volume changes • Triaxial test data are required for a more accurate calibration. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.69

Extended Drucker-Prager Models • A triaxial machine enables application of a confining pressure and a differential stress. • Several tests covering the range of confining pressures of interest are usually performed. • The stress and strain in the direction of loading are recorded, together with the lateral strain, to enable calibration of volume changes.


L3.70

Extended Drucker-Prager Models • In a triaxial compression test the specimen is confined by pressure and an additional compression stress is superposed in one direction. • Thus, the principal stresses are all negative, with 0 ≥ σ1 = σ2 ≥ σ3.

Triaxial compression and tension tests

• The stress invariant values in triaxial compression are

1 ( 2σ 1 + σ 3 ), q = σ 1 − σ 3, r = −q, t = q. 3 • The triaxial results can, thus, be plotted in the q–p plane. p=−


L3.71

Extended Drucker-Prager Models • The stress state corresponding to some user-chosen critical level provides one data point for calibrating the yield surface material parameters. • Choose the stress at the onset of inelastic behavior or the ultimate yield stress. • Additional data points are obtained from triaxial tests at different levels of confinement. • These data points define the shape and position of the yield surface in the meridional plane.


L3.72

Extended Drucker-Prager Models • Defining the shape and position of the yield surface is adequate to define the model if it is to be used as a failure surface. • To incorporate isotropic hardening, one of the stress-strains curves from the triaxial tests can be used to define the hardening behavior. • The curve that represents hardening most accurately over a wide range of loading conditions should be chosen. • Unloading measurements in these tests are useful to calibrate the elasticity, particularly in cases where the initial elastic region is not well defined.


L3.73

Extended Drucker-Prager Models • Linear Drucker-Prager model • Fitting the best straight line through the results provides the friction angle β.


L3.74

Extended Drucker-Prager Models • Triaxial tension test data are also needed to define K: • Confine the specimen by pressure and then reduce the pressure in one direction. • In this case the principal stresses are 0 ≥ σ1 ≥ σ2 = σ3. • The stress invariants are now

p=−

q 1 (σ 1 + 2σ 3 ) , q = σ 1 − σ 3, r = q, t = . 3 K

K can be found by plotting q versus p. – K is the ratio of the values of q for triaxial tension and compression at the same value of p. • The dilation angle ψ is obtained from shear tests, and it must be chosen such that a reasonable match of the volume changes during yielding is obtained. • Generally, 0 ≤ ψ ≤ β. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.75

Extended Drucker-Prager Models • Hyperbolic model • Use the triaxial compression results at high confining pressures to obtain β and d' for the hyperbolic model. • Hydrostatic tension, pt , is also needed to complete the calibration. β q

d′

− d ′ tan β

− pt

Hyperbolic :

p F =

(d′

0

− pt

tan β ) + q − p tan β − d ′ = 0 2

0

2


L3.76

Extended Drucker-Prager Models • Exponent model • Abaqus provides a capability to determine the material parameters a, b, and pt required for the exponent model from triaxial data: • A “best fit” of the triaxial test data at different levels of confining stress is performed.


L3.77

Extended Drucker-Prager Models • The data points obtained from triaxial tests are specified using the ∗TRIAXIAL TEST DATA option. • The TEST DATA parameter is required on the ∗DRUCKER PRAGER option to use this feature. • The ∗TRIAXIAL TEST DATA option must be used with the ∗DRUCKER PRAGER option. • This capability allows for all three parameters, a, b, and pt , to be calibrated, or, if some of the parameters are known, to calibrate only the unknown parameters.


L3.78

Extended Drucker-Prager Models • Example: Limit Load Analysis of a Strip Foundation • This is a continuation of the limit load analysis of a strip foundation example previously analyzed with the Mohr-Coulomb material model. • Review of the model: • A rigid footing 10 feet wide on a granular soil material of 12 feet depth is loaded and the average pressure vs. displacement of the footing is computed. • The model is assumed to be in a plane strain condition, and uses CPE8R and CINPE5R elements. Right half of footing

Soil

Symmetry plane

• The footing is driven by prescribed displacement. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.79

Extended Drucker-Prager Models • Extended Drucker-Prager model: • Elastic modulus 30,000 psi and Poisson’s ratio 0.3 • Friction angle 30.64o for dilatant flow and 30.16 for nondilatant flow • Stress ratio 1.0 • Yield stress 20.2 psi for dilatant flow and 19.8 psi for nondilatant flow • Dilation angles used: 0o (nondilatant flow) and 30.16o (dilatant flow) • These material parameters have been selected such that they match with the Mohr-Coulomb material data under plane strain conditions. • See “Matching Mohr-Coulomb parameters to the Drucker-Prager model,” Section 18.3.1 of the Abaqus Analysis User’s Manual.


L3.80


• The Drucker-Prager material data is specified as follows:

*MATERIAL,NAME= A1


*ELASTIC 30000.,0.3 *DRUCKER PRAGER,SHEAR CRITERION=LINEAR Friction angle, stress ratio, and dilation angle 30.64,1.0,0. *DRUCKER PRAGER HARDENING Yield stress and plastic strain 20.2,0.

• A prescribed displacement of 5 inches is applied and the reaction is computed. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L3.81

Extended Drucker-Prager Models • Results Terzaghi (175 psi)

Prandtl (143 psi) Drucker Prager with no dilation Drucker Prager with dilation Mohr Coulomb with no dilation Mohr Coulomb with dilation


L3.82

Extended Drucker-Prager Models • Results • It is observed that the Drucker-Prager results for the dilatant and nondilatant cases are close to the respective Mohr-Coulomb solutions. • The nondilatant Drucker-Prager and Mohr-Coulomb models lead to a softer response and a lower limit load than the corresponding dilatant versions. • Presence of dilation results in stiffening of a volumetrically compressed material. • The limit load values obtained for the nondilatant Drucker-Prager and Mohr-Coulomb models lie between the Prandtl and Terzaghi solutions.


Constitutive Models: Part 2 Lecture 4


L4.2

Overview • Modified Drucker-Prager/Cap Model • Critical State (Clay) Plasticity Model • Jointed Material Model • Soil Plasticity Models - Summary • Comments on the Numerical Implementation


Modified Drucker-Prager/Cap Model


L4.4

Modified Drucker-Prager/Cap Model • The characteristics of the Modified Drucker-Prager/Cap model: • This model is intended to simulate the constitutive response of cohesive geological materials. • It adds a “cap” yield surface to the linear Drucker-Prager model • to bound the model in hydrostatic compression, and • to help control volume dilatancy when the material yields in shear. • The elastic response is followed by a non-recoverable response idealized as being plastic. • The material is initially isotropic.


L4.5


• The yield behavior depends on the hydrostatic pressure. • There are two distinct regions of behavior. • On the failure surface the material is perfectly plastic. • On the cap yield surface it hardens and also stiffens. The hardening/softening behavior is a function of the volumetric plastic strain.

• The yield behavior may be influenced by the magnitude of the intermediate principal stress.


L4.6

Modified Drucker-Prager/Cap Model • The inelastic behavior is generally accompanied by volume changes. • On the failure surface the material dilates. • On the cap surface it compacts. • At the intersection of these surfaces, the material can yield indefinitely at constant shear stress without changing volume. • Under large stress reversals the model provides reasonable material response on the cap region; • however, on the failure surface region the model is acceptable only for essentially monotonic loading. • The material properties can be temperature dependent.


L4.7

Modified Drucker-Prager/Cap Model • Description • The model can use linear elasticity or nonlinear porous elasticity. • The model uses two main yield surface segments. • a linearly pressure-dependent Drucker-Prager shear failure surface • a compression cap yield surface • The Drucker-Prager failure surface itself is perfectly plastic (no hardening), but plastic flow on this surface produces inelastic volume increase, which causes the cap to soften. • The Drucker-Prager failure surface is

Fs = t − p tan β − d = 0. • β is the angle of friction; d is the cohesion of the material. • t is the measure of the deviatoric stress, and it allows matching of different stress values in tension and compression in the deviatoric plane. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.8



L4.9

Modified Drucker-Prager/Cap Model • The cap yield surface is

Fc =

⎡

2

⎤ Rt ⎥ − R ( d + pa tan β ) = 0, ⎣ (1 + α − α cos β ) ⎦

( p − pa )2 + ⎢

where

R is a material parameter that controls the shape of the cap,

α is a small number, and

( )

pl pa ε vol is an evolution parameter that represents the volumetric

plastic strain driven hardening/softening.


L4.10


• The cap yield surface • is elliptical with constant eccentricity in the meridional ( p–t ) plane. • includes dependence on the third stress invariant in the deviatoric plane. • hardens or softens as a function of the volumetric plastic strain. • Volumetric plastic compaction (when yielding on the cap) causes hardening. • Volumetric plastic dilation (when yielding on the shear failure surface) causes softening.


L4.11

Modified Drucker-Prager/Cap Model • The hardening/softening law is a user-defined piecewise linear function relating the hydrostatic compression yield stress, pb, and the corresponding volumetric plastic strain,

(

)

pl pl pb = pb ε vol + ε vol . 0


L4.12

Modified Drucker-Prager/Cap Model • The volumetric plastic strain axis in the hardening curve has an arbitrary origin. pl • ε vol is the position on this axis corresponding to the initial state of 0 the material when the analysis begins, thus defining the position of

the cap ( pb) at the start of the analysis. • The evolution parameter pa is then given as

pa =

pb − Rd . (1 + R tan β )

• The parameter α is a small number (typically 0.01 to 0.05) used to define a transition yield surface

Ft =

( p − pa )

2

2

⎡ ⎛ ⎤ α ⎞ + ⎢t − ⎜ 1 − ⎟ ( d + pa tan β ) ⎥ − α ( d + pa tan β ) = 0, β cos ⎠ ⎣ ⎝ ⎦

so that the model provides a smooth intersection between the cap and the failure surfaces. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.13

Modified Drucker-Prager/Cap Model • A larger value of α can be used to construct more complex (curved) failure surfaces.

• A large value of α can provide a softening surface with a better fit of experimentally observed shear failure. • If α = 0 then there is no transition surface and no true softening behavior (i.e., softening while yielding on the same surface) in the model.


L4.14

Modified Drucker-Prager/Cap Model • The plastic flow is defined by a flow potential that is • associated in the deviatoric plane, • associated in the cap region in the meridional plane, and • nonassociated in the failure surface and transition regions in the meridional plane. The flow potential surface used in the meridional plane is made up of an elliptical portion in the cap region that is identical to the cap yield surface: 2

Gc =

( p − pa )


2

⎡ ⎤ Rt +⎢ ⎥ . ⎣ (1 + α − α cos β ) ⎦

L4.15

Modified Drucker-Prager/Cap Model • In the failure and transition regions another elliptical portion provides the nonassociated flow component in the model: 2

⎡ ⎤ t Gs = ⎡⎣( pa − p ) tan β ⎤⎦ + ⎢ ⎥ . ⎣ (1 + α − α cos β ) ⎦ 2

• The two elliptical portions form a continuous and smooth potential surface.


L4.16

Modified Drucker-Prager/Cap Model • On the failure and transition regions:

• The volumetric plastic strain rate is proportional to ( pa − p ) tan β . • The flow is purely deviatoric at the intersection of the DruckerPrager cone with the cap (i.e., the critical state condition). • The deviatoric plastic strain rate is proportional to t . • The flow is purely volumetric at the apex of the Drucker-Prager cone.

• Along the Drucker-Prager failure surface between these two extremes, the ratio between the volumetric and deviatoric plastic strain rates varies linearly.


L4.17


• If the initial stress is given such that the stress point lies outside the initially defined cap or transition yield surfaces, then Abaqus will try to adjust the initial position of the cap to make the stress point lie on the yield surface and a warning will be issued. • If the stress point lies outside the Drucker-Prager failure surface, then an error message will be issued and the analysis will be terminated.


L4.18

Modified Drucker-Prager/Cap Model • Usage • The Modified Cap model is invoked with the ∗CAP PLASTICITY material option. • It also allows for the yield and flow rule parameters to be made dependent on temperature and predefined field variables.


L4.19

Modified Drucker-Prager/Cap Model • The ∗CAP PLASTICITY option must always be accompanied by the ∗CAP HARDENING material option. • It defines the evolution of the yield stress in hydrostatic compression. • It can be temperature and field dependent. • The range of values for which pressure yield stress, pb , is defined should be sufficient to include all values of effective pressure stress that the material will be subjected to during the analysis.


L4.20

Modified Drucker-Prager/Cap Model • The elastic behavior is defined with linear elasticity or with porous elasticity. • The elasticity parameters can be temperature and field dependent. – Initial void ratio (porosity) is required if nonlinear elasticity is used. • User subroutine VOIDRI can be used to specify complex initial void ratio distributions. – Thermal expansion can be used to introduce thermal volume change effects.


L4.21

Modified Drucker-Prager/Cap Model • The nonassociated flow present in the failure surface of the model produces an unsymmetric material stiffness matrix. • The unsymmetric solver should be used if there is significant plastic flow due to shearing. • When Abaqus uses a symmetric solver with unsymmetric equations, the rate of convergence may be slow. • If the region of the model in which nonassociated flow occurs is confined, it may be possible to use the default symmetric solver and still get an acceptable rate of convergence.


L4.22

Modified Drucker-Prager/Cap Model • Calibration • At least three experiments are required to calibrate the simplest version of the Cap model: • A hydrostatic compression test • A uniaxial strain (oedometer) test is also acceptable. • Two uniaxial and/or triaxial compression tests • Use more than two tests for a more accurate calibration. • The hydrostatic compression test is performed by pressurizing the sample equally in all directions. • The applied pressure and volume change are recorded. • It provides the evolution of the hydrostatic compression yield stress, pl pb ε vol , required for the cap hardening option.

( )


L4.23

Modified Drucker-Prager/Cap Model • The uniaxial compression test involves compressing the sample between two rigid platens. • Record the load and displacement in the direction of loading. • The lateral displacements should also be recorded so that correct volume changes can be calibrated. • Triaxial compression experiments are performed using a standard triaxial machine where a fixed confining pressure is maintained while a differential stress is applied. • Perform several tests covering the range of confining pressures of interest. • Record the stress and strain in the direction of loading, together with the lateral strain, so that correct volume changes can be calibrated. • Unloading measurements in these tests are useful to calibrate the elasticity, particularly in cases where the initial elastic region is not well defined. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.24

Modified Drucker-Prager/Cap Model • To obtain the friction angle, β, and cohesion, d, which define the shear failure dependence on hydrostatic pressure, • Plot the failure stresses of any two uniaxial and/or triaxial compression experiments in the pressure stress (p) versus shear stress (q) space. • The slope of the straight line passing through the two points gives the angle β, and the intersection with the q-axis gives d.


L4.25

Modified Drucker-Prager/Cap Model • K controls the yield dependence on the third stress invariant, and its calibration needs experimental results from a true triaxial (cubical) test. • True triaxial tests are difficult to conduct. • Assume a value of K (the value of K is generally between 0.8 and 1.0) or ignore this effect.

• R represents the curvature of the cap part of the yield surface. • It can be calibrated from a number of triaxial tests at high confining pressures (in the cap region).

• R must be between 0 and 1.

Y

R=

X Y X


L4.26

Modified Drucker-Prager/Cap Model • Example: Limit Load Analysis of a Strip Foundation • This is a continuation of the limit load analysis of a strip foundation example previously analyzed using the Mohr-Coulomb and DruckerPrager models. • Review of the model: • A rigid footing 10 feet wide on a granular soil material of 12 feet depth is loaded by a prescribed displacement of 5 inches and the average pressure vs. displacement of the footing is computed. • The model is assumed to be in a plane strain condition and is discretized by CPE8R and CINPE5R elements. Right half of footing

Soil

Symmetry plane Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.27

Modified Drucker-Prager/Cap Model • Cap hardening is used for the Modified DruckerPrager/Cap model. • Adds a cap yield surface to the modified Drucker-Prager model. • The cap bounds the yield surface in hydrostatic compression, thus providing an inelastic hardening mechanism to represent plastic compaction. • The cap also helps control volume dilatancy when the material yields in shear by providing softening as a function of the inelastic volume increase created as the material yields on the Drucker-Prager shear failure and transition yield surfaces. • The model uses associated flow in the cap region and nonassociated flow in the shear failure and transition regions. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.28

Modified Drucker-Prager/Cap Model • The Drucker-Prager material data is specified as follows: *MATERIAL,NAME= A1 *ELASTIC


30000.,0.3 *CAP PLASTICITY 16.212,30.64,0.1,.00041,.01,1.0 *CAP HARDENING 2.15,0. 20.96,.0005

Cohesion, Friction angle, Cap eccentricity, Initial volumetric plastic strain, Transition surface radius parameter, Stress ratio

46.6,.001 79.67,.0015 126.28,.002 205.95,.0025 311.27,.0028 655.6,.00299


Table of pressure yield stress and volumetric inelastic strain

L4.29

Modified Drucker-Prager/Cap Model • The Drucker-Prager material data is specified in Abaqus/CAE as follows:

Cohesion, Friction angle, Cap eccentricity, Initial volumetric plastic strain, Transition surface radius parameter, Stress ratio

Table of pressure yield stress and volumetric inelastic strain


L4.30

Modified Drucker-Prager/Cap Model • Results Terzaghi (175 psi)

Prandtl (143 psi)

Drucker- Prager/Cap Drucker Prager with no dilation Drucker Prager with dilation Mohr Coulomb with no dilation Mohr Coulomb with dilation Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.31

Modified Drucker-Prager/Cap Model • Results • It is observed that the Modified Drucker-Prager/Cap model provides a response that is comparable to the corresponding Drucker-Prager nondilatant response. • This response is due to the addition of the cap and the nonassociated flow in the failure region, which combine to reduce dilation in the model and therefore approximate the Drucker-Prager nondilatant flow. • The results with the Modified Drucker-Prager/Cap model are also close to the results obtained from the nondilatant Mohr-Coulomb model. • The Modified Drucker-Prager/Cap model as well as the nondilatant DruckerPrager and Mohr-Coulomb models lead to a softer response and a lower limit load than the corresponding dilatant versions. • The limit load values obtained for the Modified Drucker-Prager/Cap model and the nondilatant Drucker-Prager and Mohr-Coulomb models lie between the Prandtl and Terzaghi solutions. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Critical State (Clay) Plasticity Model


L4.33

Critical State (Clay) Plasticity Model • The characteristics of the Critical State (Clay) Plasticity model: • This model is intended to simulate the constitutive behavior of cohesionless materials. • It is an extended version of the classical critical state theories originally developed at Cambridge from 1960–1970. • The elastic response is followed by a non-recoverable response idealized as being plastic. • The material is initially isotropic. • The yield behavior depends on the hydrostatic pressure. • The critical state line separates two distinct regions of behavior. • On the “dry” side of critical state the material softens. • On the “wet” side it hardens (and also stiffens). • The hardening/softening behavior is a function of the volumetric plastic strain.


L4.34

Critical State (Clay) Plasticity Model • The inelastic behavior is generally accompanied by volume changes. • On the “dry” side the material dilates. • On the “wet” side it compacts. • On the critical state line the material can yield indefinitely at constant shear stress without changing volume. • The yield behavior may be influenced by the magnitude of the intermediate principal stress. • The model assumes that the material is cohesionless. • Under large stress reversals the model provides a reasonable material response on the “wet” (cap) side of the critical state. • On the “dry” side the model is acceptable only for essentially monotonic loading. • The material properties can be temperature dependent.


L4.35


Response on the “dry” side of critical state Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.36


Response on the “wet” side of critical state


L4.37

Critical State (Clay) Plasticity Model • Description of the Critical State (Clay) Plasticity Model: • Either linear elasticity or nonlinear porous elasticity can be used with the Clay model. • The Clay yield surface is elliptical in the meridional plane and also includes dependence on the third stress invariant. • for p > a (the “wet” side of critical state),

f ( p , q, r ) =

2

2

1 ⎛p ⎞ ⎛ t ⎞ ⎜ − 1⎟ + ⎜ ⎟ −1 = 0 β 2 ⎝ a ⎠ ⎝ Ma ⎠

• for p ≤ a (the “dry” side of critical state), 2

2

⎛p ⎞ ⎛ t ⎞ f ( p , q , r ) = ⎜ − 1⎟ + ⎜ ⎟ −1 = 0 ⎝ a ⎠ ⎝ Ma ⎠


L4.38



L4.39

Critical State (Clay) Plasticity Model • β is a constant used for modifying the yield surface shape on the “wet” side of the critical state. •

β = 1 on the “dry” side of the critical state.

•

β < 1 in most cases on the “wet” side.

• The measure of deviatoric stress, t, allows matching of different stress values in tension and compression in the deviatoric plane.


L4.40

Critical State (Clay) Plasticity Model M is the slope of the critical state line in the p–t plane (the ratio of t to p at critical state).

Clay yield and critical state surfaces in principal stress space


L4.41

Critical State (Clay) Plasticity Model • Associated flow is used with the Clay model. • The size of the yield surface is defined by a. • The evolution of a characterizes the hardening or softening of the material.

( )

pl • Two approaches can be used to define the evolution of a ε vol .

• The first approach: • For some materials, over the range of confining pressures, ( p ) , of interest, it is observed that during plastic deformation,

de = −λ d ( ln p ) , where λ is the logarithmic hardening constant and e is the void ratio.


L4.42

Critical State (Clay) Plasticity Model ⎡

• Integrating this equation we obtain, a = a0 exp ⎢(1 + e0 )

⎣

1 − J pl ⎤ ⎥, λ − κ J pl ⎦

where e0 is the initial void ratio, κ is the porous elasticity volumetric constant (the logarithmic bulk modulus), and a0 defines the position of a at the start of the analysis—the initial overconsolidation of the material. • The value of a0 can be given directly or may be computed as

1 ⎛ e − e − κ ln p0 ⎞ a0 = exp ⎜ 1 0 ⎟, λ −κ 2 ⎝ ⎠ where p0 is the initial value of the pressure stress and e1 is the intercept of the virgin consolidation line with the void ratio axis in a plot of void ratio versus pressure stress. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.43



L4.44


( )

pl • The second approach to define the evolution of a ε vol :

• The evolution of the yield surface can be defined as a piecewise linear function relating yield stress in hydrostatic compression, pc, pl and volumetric plastic strain, ε vol :

( )

pl pc = pc ε vol .

The evolution parameter is then given by:

a=


pc . (1 + β )

L4.45

Critical State (Clay) Plasticity Model • The volumetric plastic strain axis has an arbitrary origin. pl • ε vol is the position on this axis corresponding to the initial state of 0 the material.

• It defines the initial hydrostatic pressure, pc 0 . • It defines the initial yield surface size, a0 .

• Data must be provided over a wide enough range of values of pc to cover all situations that will arise in the application.


L4.46

Critical State (Clay) Plasticity Model • Abaqus checks that the yield surface is not violated by: • the user-specified material constants, and • the initial effective stress conditions defined at each point in the material. • At any material point where the yield function is violated: • a warning message is issued, and

• a0 is adjusted so that the yield function is satisfied exactly. • Thus, the initial stress state lies on the yield surface.


L4.47



L4.48

Critical State (Clay) Plasticity Model • Usage • The Clay model is invoked with the ∗CLAY PLASTICITY option. • It defines the yield. • It defines the hardening if logarithmic hardening is chosen. • The ∗CLAY HARDENING option can be used to define piecewise linear hardening. • These parameters can be entered as temperature or field dependent.


L4.49

Critical State (Clay) Plasticity Model • The elasticity is defined as: • linear elastic or • porous elastic. • If clay plasticity with logarithmic hardening is used, then porous elasticity must be used. • In the porous elasticity case el zero tensile strength pt = 0 is normally used since the material is assumed to be cohesionless.

(

)


L4.50

Critical State (Clay) Plasticity Model • *EXPANSION can be used to introduce thermal volume change effects. • *INITIAL CONDITIONS, TYPE=RATIO is required to define the initial void ratio of the material. • User subroutine VOIDRI can be used to specify complex initial void ratio distributions. • *INITIAL CONDITIONS, TYPE=STRESS is required to define the initial effective stress state in the material. • User subroutine SIGINI can be used to specify complex initial stress distributions. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.51

Critical State (Clay) Plasticity Model • Calibration • At least two experiments are required to calibrate the simplest version of the Clay model: • a hydrostatic compression test • A uniaxial strain (oedometer) test is also acceptable. • and a triaxial compression test • More than one triaxial test is useful for a more accurate calibration. • The hydrostatic compression test is performed by pressurizing the sample equally in all directions. • The applied pressure and the volume change are recorded.


L4.52

Critical State (Clay) Plasticity Model • Triaxial compression experiments are performed using a standard triaxial machine. • A fixed confining pressure is maintained while a differential stress is applied. • Several tests covering the range of confining pressures of interest are usually performed. • Again, the stress and strain in the direction of loading are recorded together with the lateral strain so that correct volume changes can be calibrated.


L4.53

Critical State (Clay) Plasticity Model • Unloading measurements in these tests are useful to calibrate the elasticity: • in cases where the initial elastic region is not well defined, and • to help determine whether a constant shear modulus or a constant Poisson’s ratio should be used for porous elasticity, and to determine their values.


L4.54

Critical State (Clay) Plasticity Model • The onset of yielding in the hydrostatic compression test provides the initial position of the yield surface, a0 . • The logarithmic bulk moduli, κ and λ, are determined from the hydrostatic compression experimental data by plotting the logarithm of pressure versus void ratio. • The void ratio, e, is related to the measured volume change as

J = exp ( ε vol ) =

(1 + e ) . (1 + e0 )

• The slope of the line obtained for the elastic regime is −κ, and the slope in the inelastic range is −λ. • For a valid model λ > κ.


L4.55

Critical State (Clay) Plasticity Model • The triaxial compression tests allow the calibration of parameters M and β.

M is the ratio of the shear stress, q, to the pressure stress, p, at critical state.

Critical state line q=Mp

It can be obtained from the stress values when the material has become perfectly plastic (critical state).

β (0 < β < 1) represents the curvature of the cap part of the yield surface. It can be calibrated from a number of triaxial tests at high confining pressures (on the “wet” side of critical state). Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.56


• Calibration of K, which controls the yield dependence on the third stress invariant, needs experimental results from a true triaxial (cubical) test. • These results are generally not available. • Assume a value of K (the value of K is generally between 0.8 and 1.0) or ignore this effect.


Jointed Material Model


L4.58

Jointed Material Model • The Jointed Material model is intended to provide a simple, continuum model for materials containing a high density of parallel joint surfaces in different orientations. • The spacing of the joints of a particular orientation is assumed to be sufficiently close compared to the characteristic dimensions of the model. • The joints can be smeared into a continuum of slip systems. • This model can be applied to the modeling of geotechnical problems where the medium of interest is composed of significantly faulted rock.


L4.59

Jointed Material Model • The basic characteristics of the Jointed Material model are: • The elastic response is followed by a non-recoverable response idealized as being plastic. • The model provides for up to three joint systems. • The joint systems may exhibit frictional sliding and may also open and close. • The material response is orthotropic when any joint system is open.

• The model includes a bulk failure mechanism based on a DruckerPrager model. • The inelastic sliding mechanisms on the joints and bulk material may be purely frictional or include dilation. • The model provides a reasonable material response under large stress reversals (including joint opening/closing and cyclic shear). • The material properties can be temperature dependent. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.60


2 1 3


L4.61

Jointed Material Model • Description • Consider a joint a oriented by the normal to the joint surface na.

taα, α = 1, 2, are two unit, orthogonal vectors in the joint surface. • The local stresses are the pressure stress and the shear stresses across the joint: p = −n ⋅ σ ⋅ n , a

a

a

2

τ aα = n a ⋅ σ ⋅ t aα .

1 3

• We define the shear stress magnitude as τ a = τ aατ aα . • The local strains are the normal strain across the joint

ε an = n a ⋅ ε ⋅ n a • and the engineering shear strain in the α-direction in the joint surface

γ aα = n a ⋅ ε ⋅ t aα + t aα ⋅ ε ⋅ n a .


L4.62

Jointed Material Model el pl • We use a linear strain rate decomposition: d ε = d ε + d ε .

• If several systems are active (we designate an active system by i, where i = b indicates the bulk material system and i = a is a joint system a):

d ε pl =

∑ dε

pl i .

i

• When all joints at a point are closed, the elastic behavior of the material is isotropic and linear. • We use a stress-based joint opening criterion, whereas joint closing is monitored based on strain. • Joint system a opens when the estimated pressure stress across the joint (normal to the joint surface) is no longer positive:

pa ≤ 0. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.63

Jointed Material Model • In this case the material has no elastic stiffness with respect to direct strain across the joint system. • It may or may not have stiffness with respect to shearing associated with this direction. • Thus, open joints create anisotropic elastic response at a point. el el • The joint system remains open so long as ε an ( ps ) ≤ ε an ,

whereε an is the component of direct elastic strain across the joint and ε an( ps ) is the component of direct elastic strain across the joint calculated in plane stress as v el ε an ( ps ) = − E (σ a1 + σ a 2 ) , el

el

where E is the Young’s modulus of the material, v is the Poisson’s ratio, and

σ aα = taα ⋅ σ ⋅ taα are the direct stresses in the plane of the joint. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.64

Jointed Material Model • The failure surface for sliding on joint a is

f a = τ a − pa tan β a − d a = 0, where β is the friction angle and da is the cohesion for system a. • fa < 0 implies joint a does not slip. • fa = 0 implies joint a slips: • The inelastic strain on the system is

⎛ ⎞ τ d ε apl = d ε apl ⎜ sinψ a n a n a + aα cosψ a ( n a t aα + t aα n a ) ⎟ , τa ⎝ ⎠ pl d ε where is the magnitude of the inelastic strain rate. a Also: ψa is the dilation angle for this joint system.

ψa = 0 provides pure shear flow on the joint. ψa > 0 causes dilation of the joint as it slips.


L4.65

Jointed Material Model • The sliding of the different joint systems at a point is independent. • Sliding on one system does not change the failure criterion or the dilation angle for any other joint system at the same point. • Joint system model:


L4.66

Jointed Material Model • In addition to the joint systems, the model includes a bulk material failure mechanism based on the linear Drucker-Prager failure criterion:

q − p tan β b − db = 0, where βb is the friction angle for the bulk material, and db is the cohesion for the bulk material. If this failure criterion is reached, the bulk inelastic flow is defined by:

∂gb , 1 1 − tanψ b ∂σ 3 gb = q − p tanψ b ,

d ε bpl = d ε bpl

1

where ψb is the dilation angle for the bulk material. • This bulk failure system is independent of the joint systems (bulk inelastic flow does not change the behavior of any joint system). Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.67

Jointed Material Model • Bulk material model:


L4.68

Jointed Material Model • Usage • The Jointed Material model is invoked with the ∗JOINTED MATERIAL material option. • The yield and flow parameters can be temperature and field dependent. • This option must be repeated for each existing system (bulk material and up to three joints); it may appear four times. • Use the JOINT DIRECTION parameter for referring to the orientation definition corresponding to the joint. • The ∗ORIENTATION option must then be used to define the joint orientation in the original configuration. • Stress and strain components will still be output in global directions unless the ∗ORIENTATION option is also used on the section definition option associated with the material definition. • When defining the bulk material, omit the JOINT DIRECTION parameter. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.69

Jointed Material Model • The ∗JOINTED MATERIAL option may appear a fifth time with the SHEAR RETENTION parameter to define a nonzero shear modulus for open joints. • The elasticity must be defined with the ∗ELASTIC, TYPE=ISOTROPIC material option since we assume that the material is linear elastic and isotropic when all joints are closed. • The material cannot be elastically incompressible (Poisson’s ratio must be less than 0.5). – ∗EXPANSION can be used to introduce thermal volume change effects. • Analyses using a nonassociated flow version of the model may require the use of the unsymmetric solver because of the resulting nonsymmetry of the plasticity equations. • The Jointed Material model is not supported by Abaqus/CAE. • Use the Keywords Editor to add the jointed material options to an isotropic elastic model. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.70

Jointed Material Model • Calibration • At least two experiments are required to calibrate the behavior of each of the existing inelastic mechanisms. • For the bulk material calibration the following experiments can be used: • Uniaxial compression test (for cohesive materials), and, • Triaxial compression tests. • The uniaxial compression test involves compressing the sample between two rigid platens. • The load and displacement in the direction of loading are recorded. • The lateral displacements are also recorded so that correct volume changes can be calibrated. • The triaxial compression tests involve the use of a standard triaxial machine where a fixed confining pressure is maintained while a differential stress is applied. • Stress and strain in the direction of loading are recorded, together with the lateral strain, so that correct volume changes can be calibrated. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.71

Jointed Material Model • For the joint systems the most common experiments are shear tests. • A normal pressure is applied across the joint and then the joint is sheared. • Several tests covering the range of pressures of interest are performed. • Stresses and strains in the normal and shear directions are recorded. • Unloading measurements in all of the above tests are useful. • These help in the calibration of elasticity, particularly in cases where the initial elastic region is not well-defined.


L4.72

Jointed Material Model • The angle of friction, β, and the cohesion, d, defining the failure stress dependence on pressure are calculated by plotting the failure stresses of any two experiments in the pressure stress versus shear stress space. • The slope of the straight line passing through the two points gives β, and the intersection with the shear stress axis gives d.


L4.73

Jointed Material Model • If more than two experiments are available, a best fit straight line over the range of interest of pressure stress must be used to calculate the slope of the failure surface. • The dilation angle, ψ , (0 < ψ < β ), measured from a shear test, must be chosen such that a reasonable match of the volume changes during yielding is obtained. • The volume changes are calculated from the strains in all directions.

• The calibration of the temperature-dependent model requires repetition of the experiments at different temperatures over the range of interest. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.74

Jointed Material Model • Example: Jointed rock slope stability • We examine the stability of the excavation of part of a jointed rock mass, leaving a sloped embankment (Example Problem 1.1.6). • This example illustrates the use of the jointed material model. • This problem has been studied previously by the following: • Barton (1971) and Hoek (1970), who used limit equilibrium methods and • Zienkiewicz and Pande (1977), who used a finite element model.


L4.75

Jointed Material Model • As in most geotechnical problems, the analysis begins from a nonzero state of stress. • The active “loading” in this case consists of removal of material to represent the excavation.

Excavation of part of a jointed rock mass

• We examine the effect of joint cohesion on slope collapse through a sequence of solutions with different values of joint cohesion, with all other parameters kept fixed.


L4.76

Jointed Material Model • The rock mass contains two sets of planes of weakness: one vertical set of joints and one set of inclined joints.

JOINT1

JOINT2

Excavation geometry

Model analyzed with cohesion (d) values up to 120 kPa Model analyzed with both nonassociated flow (shown) and associated flow (ψ =β = 45º). Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.77

Jointed Material Model A

B

Deformed shape: d=30, nonassociated flow Deformation scale factor = 3,000 Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.78

Jointed Material Model • The variation of horizontal displacements as cohesion is reduced suggest that the slope collapses if the cohesion is less than 24 kPa for the case of associated flow or less than 26 kPa for the case of non-dilatant flow.

Horizontal displacement Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.79

Jointed Material Model • The cohesion values for collapse of the slope (24 kPa for associated flow or 26 kPa for non-dilatant flow) compare well with the value calculated by Barton (26 kPa) • Barton used a planar failure assumption in his limit equilibrium calculations. • Barton's calculations also include “tension cracking” (akin to joint opening with no tension strength) as we do. • Hoek calculated a cohesion value of 24 kPa for collapse of the slope. • Although he also makes the planar failure assumption, he does not include tension cracking. • This may be why Hoek’s calculated value is lower than Barton's. • Zienkiewicz and Pande calculated the cohesion value necessary for collapse as 23 kPa for associated flow and 25 kPa for non-dilatant flow. • They assumed the joints have a tension strength of 1/10 of the cohesion. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.80


Joint set 1 (vertical joints)

Joint set 2 (inclined joints)

plastic strain (d=30, nonassociated flow) Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Soil Plasticity Models - Summary


L4.82

Soil Plasticity Models - Summary • Soil plasticity models in Abaqus Model

Features

Abaqus/Standard

Abaqus/Explicit

Mohr Coulomb

Open surface model, suitable for sand as well as clays

Yes

No

Extended DruckerPrager models

Open surface model, suitable for sand as well as clays

Yes

Yes

Modified DruckerPrager/Cap model

Closed surface model with a cap, exhibits hardening due to pressure

Yes

Yes

Critical state (clay) plasticity model

Based on the critical state concept, suitable for sands or materials without cohesion

Yes

No

Jointed material model

Three independent material behaviors along the planes of joints plus one bulk material behavior

Yes

No


Comments on the Numerical Implementation


L4.84

Comments on the Numerical Implementation • For all the material models presented in this set of notes, the constitutive behavior is given in rate form. • The inelastic strain is defined only as a strain rate. • This strain rate must be integrated over each finite time increment. • We use backward Euler integration:

∆ε pl =

d ε pl dt

∆t, t +∆t

where ∆t is the time increment and t is the time at the beginning of the increment. • The backward Euler method is chosen because it is unconditionally stable. In addition, it provides acceptable accuracy, especially when the strain increment is large compared to the strain to cause yield. • The backward Euler method is more traditionally called “radial return” because of its simple geometric interpretation for the basic case of a Mises yield and associated flow model. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L4.85

Comments on the Numerical Implementation • In some simple cases, such as a perfectly plastic Drucker-Prager model, the backward difference equations can be solved in closed form. • However, in most cases, Newton’s method is used for the numerical solution of the integrated plasticity equations. • For some material models additional logic is used to handle all the possible combinations of behavior available within the model. • This is required for models that involve opening/closing conditions, such as the Jointed Material model. • This is required for models that involve more than one yield surface, such as the Modified Drucker-Prager/Cap model.


L4.86

Comments on the Numerical Implementation • Occasionally Abaqus fails to find a solution to the integrated constitutive equations. • It then gives a warning message: “PLASTICITY ALGORITHM FAILS TO CONVERGE AT n INTEGRATION POINTS.” • If automatic time incrementation has been chosen (which is always recommended), Abaqus tries again with a smaller time increment and the problem is, thereby, resolved. • If necessary, detailed information about where the algorithms are failing can be obtained by setting ∗PRINT, PLASTICITY=YES.


L4.87

Comments on the Numerical Implementation • Abaqus provides consistent material Jacobians (tangent stiffnesses) for the global equilibrium iterations so that quadratic convergence can be achieved.


Analysis of Porous Media Lecture 5


L5.2

Overview • Basic Assumptions and Effective Stress • Stress Equilibrium and Flow Continuity • Types of Analyses and Usage • Saturated Example Problems • Partially Saturated Example Problems


L5.3

Overview • Abaqus has capabilities for modeling single phase flow through porous media. • The following cases can be analyzed: • Fully saturated flow (encountered in many geotechnical applications) • Partially saturated flow (encountered in irrigation problems and hydrology problems) • A combination of the above two (calculation of phreatic surfaces) • Fluid gravity effects may or may not be considered. • Total pore pressure analysis includes fluid weight. • This analysis is required when the loading provided by fluid gravity is large. • It is also required when the “wicking” phenomenon (transient capillary suction of liquid into a dry body) is of interest. • Excess pore pressure analysis disregards fluid weight. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.4

Overview • Two other effects, gel swelling and moisture swelling, may be included in partially saturated cases. • These swelling effects are usually associated with moisture absorption into polymeric systems (such as paper towels) rather than with geotechnical systems and are hence not discussed in these notes.


Basic Assumptions and Effective Stress


L5.6

Basic Assumptions and Effective Stress • We model a deforming porous medium using the conventional approach. • The medium is considered to be a multi-phase material. • An effective stress principle is used to describe the material behavior. • An elementary volume, dV, is made up of a volume of grains of solid material, dVg, and a volume of voids, dVv, which is either fully or partly saturated with a volume of wetting fluid, dVw.


L5.7

Basic Assumptions and Effective Stress • The porosity of the medium, n, is defined as the ratio of the volume of voids to the total volume:

n=

dVv . dV

(

)

• Abaqus generally uses void ratio, e = dVv dVg , instead of porosity. • Conversion relationships are

e=

n e 1 , n= , 1− n = . 1− n 1+ e 1+ e


L5.8

Basic Assumptions and Effective Stress • Saturation, s, is defined as the ratio of wetting fluid volume to void volume:

s=

dVw . dVv

• For a fully saturated medium s = 1, while for a completely dry medium s = 0.


L5.9

Basic Assumptions and Effective Stress • The total stress acting at a point, σ, is assumed to be made up of: • An average pressure stress in the wetting fluid, uw, called the “pore pressure” times a factor, χ, that depends on the saturation • An “effective stress” in the material skeleton,σ :

σ = σ + χ uwI. • In general, χ = χ(s) can be measured experimentally. • Typical experimental data are shown at right.


L5.10

Basic Assumptions and Effective Stress • Since such data are difficult to obtain, Abaqus assumes χ = s.

σ = σ + suwI. • We assume the following: • The constitutive response of the porous medium consists of simple bulk elasticity relationships for the fluid and for the solid grains. • A constitutive theory for the material skeleton is chosen wherein σ is defined as a function of the strain history and temperature of the material.

σ =σ

(strain history, temperature, state variables).


L5.11

Basic Assumptions and Effective Stress • Any of the constitutive models in Abaqus can be used to model the material skeleton of voided materials. • The strain rate decomposition is then given by

(

)

d ε = d ε gvol + d ε wvol I + d ε el + d ε pl , where d ε gvol , d ε wvol are the volume strain rates in the solid grains and fluid and d ε el , d ε pl are the elastic and plastic strain rates in the material skeleton.


Stress Equilibrium and Flow Continuity


L5.13

Stress Equilibrium and Flow Continuity • Stress equilibrium for the solid phase of the material is expressed by writing the principle of virtual work for the volume under consideration in its current configuration at time t:

∫ (σ − χ u I ) : δε dV = ∫ t ⋅ δ vdS + ∫ f ⋅ δ vdV + ∫ sn ρ g ⋅ δ vdV , w

V

where

w

S

V

V

def

δε = sym ( ∂δ v ∂x ) is the virtual rate of deformation, σ is the true (Cauchy) effective stress,

δ v is a virtual velocity field, t are surface tractions per unit area, f are body forces (excluding fluid weight) per unit volume,

ρw is the density of the fluid, and g is the gravitational acceleration (assumed constant and in a fixed direction). Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.14

Stress Equilibrium and Flow Continuity • This equation is then discretized using a Lagrangian formulation for the solid phase, with displacements as the nodal variables. • The porous medium is thus modeled by attaching the finite element mesh to the solid phase. • Fluid may flow through this mesh.


L5.15

Stress Equilibrium and Flow Continuity • A continuity equation is therefore required for the fluid, equating the rate of increase in fluid volume stored at a point to the rate of volume of fluid flowing into the point within the time increment:

⎞ ρ d ⎛ ρw ⎜ ⎟ = − w0 sn n ⋅ v w dS , sn dV 0 ⎟ dt ⎜ ρ w ρ S w ⎝V ⎠

∫

where

∫

• vw is the average velocity of the fluid relative to the solid phase (the seepage velocity) and

• n is the outward normal to S. • This equation has been normalized by ρ w , the reference density of the fluid. 0

• The continuity equation is integrated in time using the backward Euler approximation and discretized with finite elements using pore pressure as the variable. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.16

Stress Equilibrium and Flow Continuity • The pore fluid flow behavior is assumed to be governed either by Darcy's law or by Forchheimer's law. • Darcy's law is generally applicable to low fluid flow velocities. • Forchheimer's law is used for higher flow velocities. • Darcy's law may be thought of as a linearized version of Forchheimer's law.


L5.17

Stress Equilibrium and Flow Continuity • Forchheimer's law describes pore fluid flow as:

vw = −

1

(

sng ρ w 1 + β v w ⋅ v w

)

⎛ ∂u ⎞ kˆ ⋅ ⎜ w − ρ wg ⎟ , ∂ x ⎝ ⎠

where

• g is the magnitude of the gravitational acceleration, • kˆ ( s, e ) is the permeability of the medium (possibly anisotropic) with units of length/time, and

• β ( e ) is a “velocity coefficient”. (Some texts use a different definition of permeability. This is discussed later.)


L5.18

Stress Equilibrium and Flow Continuity • Darcy's law is obtained by setting β = 0. • We see that, as the fluid velocity tends to zero, Forchheimer's law approaches Darcy’s law. • The permeability depends on the saturation of the fluid and on the porosity of the medium. We assume these dependencies are separable, so that

kˆ = ks k, where

• ks(s) gives the saturation dependency, with ks(1) = 1.0, and k(e) is the fully saturated permeability. • For isotropic materials k = kI.


L5.19

Stress Equilibrium and Flow Continuity • Experimental observation often suggests that, in steady flow through a partially saturated medium, the permeability varies with s3.

• We, therefore, take ks = s3 by default. • Different forms of behavior for ks(s) can be defined by using the ∗PERMEABILITY option. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.20

Stress Equilibrium and Flow Continuity • Because uw measures pressure in the wetting fluid, the medium is fully saturated for uw > 0. • Negative values of uw represent capillary effects in the medium. • For uw < 0 it is known that, at a given value of capillary pressure, − uw, the saturation lies within certain limits.


L5.21

Stress Equilibrium and Flow Continuity • These limits are defined by using the ∗SORPTION option. • Typical forms are shown below:


L5.22

Stress Equilibrium and Flow Continuity a e • We write these limits as s ≤ s ≤ s , where

• sa(uw) is the limit at which absorption will occur (so that s > 0 ), and • se(uw) is the limit at which exsorption will occur; thus, s < 0. • We assume that these relationships are uniquely invertible and can also be written as • uwa ( s ) during absorption, and • uwe ( s ) during exsorption. • We also assume that some wetting fluid will always be present in the medium: s > 0. • The transition between absorption and exsorption, and vice-versa, takes place along “scanning” lines, which are approximated by a single value of duw/ds at all saturation levels.


L5.23

Stress Equilibrium and Flow Continuity • Saturation is treated as a state variable that may have to change in certain situations. • Such a situation arises when the wetting liquid pressure is outside the range for which its value is admissible according to the actual data.


L5.24

Stress Equilibrium and Flow Continuity • Note on permeability units: In Abaqus we define permeability in the flow constitutive equation

vw = −

(

1

sng ρ w 1 + β v w ⋅ v w

)

⎛ ∂u ⎞ kˆ ⋅ ⎜ w − ρ wg ⎟ ⎝ ∂x ⎠

as kˆ , with units of length/time. • In this equation

ρw is the mass density of the fluid, • s is the saturation, • n is the porosity,

β is the velocity coefficient, and • g is the gravitational acceleration. • It is then clear that both sides of the equation have units of length/time.


L5.25

Stress Equilibrium and Flow Continuity • However, one other definition of permeability (K) is often used:

K=

µ ˆ k,

g ρw

where µ is the fluid viscosity in poise units (mass/time-length). • In this context, permeability K has length squared units (or Darcy) and what we refer to in Abaqus as permeability, k ,ˆ is called the hydraulic conductivity.


L5.26

Stress Equilibrium and Flow Continuity • In the coupled problem the stress equilibrium and fluid continuity equations must be solved simultaneously. • In the general nonlinear case we use a Newton scheme to solve the equations. • The Newton equations for various cases of the formulation (transient or steady-state flow, etc.) are described in Appendix A.


Types of Analyses and Usage


L5.28

Types of Analyses and Usage • The analysis of flow through porous media in Abaqus is available for plane strain, axisymmetric, and three-dimensional problems. • Special coupled displacement/pore pressure elements must be used. • These elements have a linear distribution of pore pressure. • These elements use either a first-order or a second-order distribution of displacement. • The coupled stress/flow problems are solved by Abaqus using the ∗SOILS procedure: • By default, Abaqus will solve the steadystate problem (∗SOILS, STEADY STATE). • The transient problem is invoked with ∗SOILS, CONSOLIDATION. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.29

Types of Analyses and Usage • The steady-state problem assumes that the time scale is long enough to ensure that there is no transient effect in the pore fluid diffusion part of the problem. • The time scale chosen is then only relevant to any possible rate effects in the constitutive model used for the material skeleton. • Mechanical loads and boundary conditions can be changed gradually over a step to accommodate nonlinearities in the response.


L5.30

Types of Analyses and Usage • Uncoupled (purely diffusion) pore pressure elements are not available in Abaqus. • Such elements are useful in cases where only the pore fluid flow part of the problem is of interest. • In Abaqus we can use coupled elements and constrain all the displacement degrees of freedom. • Coupled deformation/pore pressure infinite elements are not available in Abaqus. • Dynamic coupled stress/fluid flow analysis cannot be performed since we assume no inertia in the existing coupled analysis capability. • The inertia is important in cases such as the behavior of a dam subjected to earthquake loading. • Liquefaction effects may arise in dynamic situations.


L5.31

Types of Analyses and Usage • Three-way coupled stress/fluid flow/temperature analysis cannot be performed. • This could be important in applications such as those encountered in nuclear waste disposal and oil reservoir simulation. • The capability for fluid flow through porous media assumes single-phase fluid flow. • Multi-phase fluid flow is important in cases such as oil reservoir simulation. • The vapor phase is ignored in the partially saturated flow formulation. • There are situations for which this assumption is not adequate.


L5.32

Types of Analyses and Usage • The transient problem includes the time integration of the diffusion effects; therefore, the choice of time increment is important. • The integration procedure used in Abaqus introduces a relationship between the minimum usable time increment and the element size. • This minimum is a requirement only for second-order elements, but it is recommended for all diffusion elements.


L5.33

Types of Analyses and Usage • A simple guideline that can be used for fully saturated flow is:

ρw g ⎛

E ∆t > ⎜1 − 6 Ek ⎜⎝ K g

2

⎞ 2 ⎟⎟ ( ∆l ) , ⎠

where:

∆t is the time increment, E is the Young's modulus of the material skeleton, k is the permeability of the saturated medium (in units of length/time), Kg is the bulk modulus of the solid grains, and ∆l is a typical element dimension.


L5.34

Types of Analyses and Usage • A simple guideline that can be used for partially saturated flow is

∆t >

ρ w gn0 ds 6k s k duw

( ∆l )2 ,

where

n0 is the initial porosity of the material, ks is the permeability-saturation relationship, and

ds/duw is the rate of change of saturation with respect to pore pressure as defined in the ∗SORPTION material option.


L5.35

Types of Analyses and Usage • Spurious oscillations may appear in the solution if time increments smaller than this value are used. • If the problem requires analysis with smaller time increments than the above relationship allows, a finer mesh is required. • Generally there is no upper limit on the time increment size except accuracy, since the integration procedure is unconditionally stable unless nonlinearities cause numerical problems.


L5.36

Types of Analyses and Usage • The accuracy of the time integration is controlled by the tolerance UTOL (maximum allowable change in pore pressure during the increment). • UTOL is also used to drive the automatic incrementation procedure for *SOILS, CONSOLIDATION analysis. • Transient analysis may be terminated by • completing a specified time period; • or it may be continued until steadystate conditions are reached, steady state being defined by all fluid pressures changing at less than a user-defined rate.


UTOL

L5.37

Types of Analyses and Usage • The porous medium coupled analysis capability can provide solutions either in terms of total or of “excess” pore fluid pressure. • The difference between total and excess pressure is relevant only for cases in which gravitational loading is important. • Total pressure solutions are provided when the GRAV distributed load type is used to define the gravity load on the model. • In total pore pressure problems the *DENSITY material option must be used to specify the density of the dry material only. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.38

Types of Analyses and Usage • The gravity contribution from the pore fluid is defined through the SPECIFIC (weight) parameter on the ∗PERMEABILITY option, together with the direction of the gravity vector specified in the ∗DLOAD option with load type GRAV. • Excess pressure solutions are provided in all other cases (for example, when gravity loading is defined with a body force –i.e., ∗DLOAD types BX, BY, or BZ).


L5.39

Types of Analyses and Usage • The ∗PERMEABILITY option is used to specify the permeability of the saturated medium, which can be isotropic or anisotropic and a function of the void ratio. • The ∗PERMEABILITY option can be repeated with TYPE=VELOCITY to invoke Forchheimer's flow law instead of the default Darcy's law. • In partially saturated cases the ∗PERMEABILITY option can also be repeated with TYPE=SATURATION to specify the dependence of permeability on saturation, ks(s). • By default, ks = s3 . Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.40

Types of Analyses and Usage • In partially saturated cases the *SORPTION option is used to specify the dependence of negative (partially saturated) pore pressure on saturation. • When it is used with TYPE=ABSORPTION (default), it defines the absorption curve. • When it is used with TYPE=EXSORPTION, it defines the exsorption curve (by default, the exsorption curve is the same as absorption). • Analytical (logarithmic) or tabular input data are permitted for these options. • In partially saturated cases the *SORPTION, TYPE=SCANNING option is used to define the scanning slope between absorption and exsorption. • Abaqus will generate this slope automatically if the user does not specify it. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.41

Types of Analyses and Usage • The *POROUS BULK MODULI material option can be used to specify the bulk modulus of the solid grains and fluid if the user wants to include the compressibility of these components of the porous medium. • *EXPANSION can be used to introduce thermal volume change effects for the solid grains and the pore fluid.


L5.42

Types of Analyses and Usage • *INITIAL CONDITIONS, TYPE=RATIO is required to define the initial void ratio (porosity) of the medium. • User subroutine VOIDRI can be used to specify complex initial void ratio distributions. • *INITIAL CONDITIONS, TYPE=STRESS can be used to define the initial effective stress state in the material. • User subroutine SIGINI can be used to specify complex initial effective stress distributions.


L5.43

Types of Analyses and Usage • *INITIAL CONDITIONS, TYPE=PORE PRESSURE can be used to define the initial pore fluid pressure in the medium. • The initial pore pressures can be defined as a linear function of elevation in the model or as a constant value. • Abaqus assumes that the vertical (elevation) direction is the 3-direction in three-dimensional models and is the 2-direction in two-dimensional (or axisymmetric) models. • User subroutine UPOREP is also available for complicated cases. • The specification of initial conditions in coupled stress partially saturated flow problems with gravity is not always trivial. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.44

Types of Analyses and Usage • The ∗DFLOW option allows the outward normal flow velocity, vn, to be prescribed across a surface. • Complex dependencies of vn on time and position can be coded in user subroutine DFLOW. • The ∗FLOW option defines the outward flow velocity as

(

)

vn = k s uw − uw∞ , ∞

where ks anduw are known values. • Again, complex conditions can be coded in user subroutine FLOW. • The ∗DFLOW and ∗FLOW options are not supported by Abaqus/CAE. • Use the Keywords Editor to include these options.


L5.45

Types of Analyses and Usage • If large-deformation analysis is required because of the presence of large strains, the NLGEOM parameter can be included on the ∗STEP option.


L5.46

Types of Analyses and Usage • The unsymmetric solver is used automatically if the user requests steadystate analysis or any kind of partially saturated flow analysis. • Abaqus automatically uses the unsymmetric solver when fluid gravity effects are included in a step. • For other simulations such as those involving nonlinear geometric effects, nonlinear permeability, or nonassociated plastic flow, • using the unsymmetric solver (∗STEP, UNSYMM=YES) may improve the rate of convergence. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.47

Types of Analyses and Usage • Two other saturation-dependent effects can be included in partially saturated flow analysis. • The ∗MOISTURE SWELLING material option can be used to define saturation driven volumetric swelling (or shrinkage as negative swelling) of the solid skeleton: • In this option, the reversible swelling strain is defined as a function of saturation. • Anisotropic swelling can be specified by using the ∗RATIOS option,

ε iims = rii

(

( ))

1 ms ε ( s ) − ε ms s I . 3


L5.48

Types of Analyses and Usage • The ∗GEL material option can be used to define the growth of gel particles that swell and trap fluid. • The growth of the gel particles depends on: • the saturation of the wetting fluid, • the size of the gel particles, and • their number per unit of volume of porous material.


Saturated Example Problems


L5.50

Saturated Example Problems • Example: Finite-strain consolidation of a two-dimensional solid • This example involves the large-scale consolidation of a two-dimensional solid (Benchmark Problem 1.14.3). • The modeled strip of soil is fully saturated and the analysis is transient. B

Boundary conditions: – Free drainage across top surface – Other surfaces impermeable and smooth

H

H

Model of soil strip Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Mesh: – 35, quadratic, reduced integration, plane strain elements with pore pressure (CPE8RP)

Loading: – Pressure = 3.4457 MPa – The loaded portion of the strip is 1/5 of the width

L5.51

Saturated Example Problems • Nonlinearities caused by the large geometry changes are considered, as well as the effects of the change in the void ratio on the permeability of the material.

*MATERIAL, NAME=Soil *ELASTIC 1000. *PERMEABILITY 2.E-6, 1. 2.E-5, 1.5


L5.52

Saturated Example Problems • The load is applied in two equal time increments during the first transient soils consolidation step, and it is kept constant thereafter. • Practical consolidation analyses require solutions across several orders of magnitude of time, and the automatic time incrementation scheme is designed to generate cost-effective solutions for such cases. • The algorithm is based on the user-supplied tolerance on the pore pressure change permitted in any increment, UTOL. • Abaqus uses this value in the following manner: • If the maximum change in pore pressure at any node is greater than UTOL, the increment is repeated with a proportionally reduced time increment. • If the maximum change in pore pressure at any node is consistently less than UTOL, the time increment size is proportionally increased.


L5.53

Saturated Example Problems • In this analysis the maximum pore pressure change per increment (UTOL) is set to 15 psi. • This is about 3% of the maximum pore pressure in the model following application of the load. • With this value the first time increment is 7.2 seconds and the final time increment is 5,437 seconds. • This is quite typical of diffusion processes: at early times the time rates of pore pressure are significant and at later times these time rates are very low. • The analysis is performed both with and without finite-strain effects. • The soil permeability varies with void ratio in the model that includes finite-strain effects. • Permeability is constant in the small-strain analysis.


L5.54

Saturated Example Problems • Results of the finite-strain analysis:

A

Point B

B

Point A

Final deformed shape


Pore pressure histories

L5.55

Saturated Example Problems • The finite-strain and small-strain analyses predict large differences in the final consolidation: Finite strain

• The small-strain result shows about 40% more deformation than the finite-strain case. • This is consistent with results from the one-dimensional Terzaghi consolidation solutions.

Small strain

• Clearly, in cases where settlement magnitudes are significant, finite-strain effects are important.

Deformation histories


L5.56

Saturated Example Problems • Example: Steady-state analysis of a dam • This problem involves the analysis of a concrete dam on a rock foundation. • The problem description provided by Dr. G. Pande of Swansea University formed the basis of the analysis. Mesh: – 110, quadratic, reduced integration, plane strain elements with pore pressure (CPE8RP)

Dam

30 m

Boundary conditions: – Exterior edges initially impervious – Symmetry at the foundation sides and bottom

Foundation

30 m

30 m

21 m

30 m

Concrete dam on rock foundation Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.57

Saturated Example Problems • We assume the rock behaves as a Drucker-Prager material with nonassociated flow and orthotropic permeability. Drucker-Prager constants: *MATERIAL, NAME=ROCK Cohesion, d = 0.1 MPa *DENSITY Friction angle, β = 40° 2400. Dilation angle, ψ = 20° *ELASTIC 30000.E6, .2 *DRUCKER PRAGER 40., 1., 20. *DRUCKER PRAGER HARDENING .139E6 *PERMEABILITY, SPECIFIC=9810., TYPE=ORTHO .0002, .00001, .0002

Vertical permeability, kv = 0.00001 m/sec Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.58

Saturated Example Problems • The dam is composed of concrete. The concrete model cracks in tension at a stress of 0.15 MPa. This is a very low tensile strength, perhaps less than one-tenth of the real tensile strength of the material. *MATERIAL, NAME=CONCRETE *DENSITY 2400. *ELASTIC 20000.E6, .25 *PERMEABILITY, SPECIFIC=9810. .00001 *CONCRETE Isotropic permeability 10.E6 k = 0.00001 m/sec 20.E6, .001 *FAILURE RATIOS 1.16, .0075 *TENSION STIFFENING 1., 0. 0., .5E-3 *SHEAR RETENTION 1., 1000. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.59

Saturated Example Problems • The analysis is done in three stages. • The first step establishes geostatic equilibrium where the geostatic stresses balance the gravity loads in the rock. • The vertical stress in the rock foundation is zero at the surface and increases linearly to approximately −0.7 MPa at the impervious bottom boundary. • The horizontal stresses are 0.6 times the vertical stress.

g

• The excess pore pressure in the rock is zero. Vertical stress

• This state of stress corresponds to the undeformed configuration of the rock mass. • This type of analysis step is discussed in Lecture 6. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.60

Saturated Example Problems • The second step of the analysis includes the following additional loads: • Gravity loads due to the concrete dam construction • Pressure loads on the upstream side of the dam and foundation representing the filling of the reservoir. • In this step we assume that no pore pressures develop, because we are concerned with short-term behavior and the modeled materials have very low permeabilities.

Deformation scale factor: 1.3E+4

g g

Short-term deformed shape


L5.61

Saturated Example Problems • In the final stage of analysis we perform a steady-state pore fluid diffusion/stress analysis to calculate the state of the structure 25 years after filling the reservoir. • Applied pore pressure boundary conditions represent the upstream and downstream conditions for the seepage part of the problem. • We also include a phreatic surface in the concrete dam wall. • The position of this zero pore pressure surface is assumed, the analysis is performed, and the validity of the assumption is checked. • An iterative procedure can then be used to correct the location of the phreatic surface at steady state. • Although not done in this case, Abaqus is capable of calculating the phreatic surface automatically.


L5.62

Saturated Example Problems • Results

Nonzero excess pore pressures corresponding to the hydrostatic pressure caused by the head of water in the full reservoir

Zero excess pore pressure phreatic surface

Zero flux boundary Pore pressure boundary conditions and long-term pore pressure contours Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.63

Saturated Example Problems • Results (cont'd) • It is observed from the pore fluid effective velocity (FLVEL) plot that the pore water flows from the upstream to the downstream of the dam. • The magnitude of FLVEL is greatest just under the concrete dam because the pressure gradient is greatest in this region of the bedrock. Pore fluid effective velocity vector plot


Partially Saturated Example Problems


L5.65

Partially Saturated Example Problems • Example: Demand Wettability of a porous medium (uncoupled) • This example illustrates the Abaqus capability to solve uncoupled partially saturated porous fluid flow problems (Benchmark Problem 1.8.1). • We consider a “constrained demand wettability” test. • The demand wettability test is a common way of measuring the absorption properties of porous materials. • In such a test fluid is made available to the material at a certain location and the material is allowed to absorb as much fluid as it can. • We consider a square specimen of material and allow it to absorb fluid at its center. • A quarter mesh of reduced-integration elements (CPE8RP) is used. • In order to perform the uncoupled fluid flow analysis, all displacement degrees of freedom are fixed.

point of absorption

Square specimen


L5.66

Partially Saturated Example Problems • We investigate two material models. • One material contains a large number of gel particles that entrap fluid and, as a result, enhance the fluid retention capability of the material: *MATERIAL,NAME=CORE *PERMEABILITY,SPECIFIC=10000. 3.7E-4, *POROUS BULK MODULI ,100000000. Gel *ELASTIC definition 100000.0, 0.0 *GEL 0.0005, 0.0015, 1.E8, 500. *SORPTION ...

• The other material does not contain gel. • We also study the cyclic wetting behavior in the case of the sample containing gel particles. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.67

Partially Saturated Example Problems • Both material models have the following absorption and exsorption behavior: *SORPTION -100000.,.04 -10000.,.05 -4500.,.1 -3500.,.18 -2000.,.45 -1000.,.91 0.,1. *SORPTION, TYPE=EXSORPTION -100000.,.09 -10000.,.1 -8000.,.11 -6000.,.18 -4500.,.33 -3000.,.79 -2000.,.91 0.,1.

Absorption/exsportion curves Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.68

Partially Saturated Example Problems • The “loading” consists of prescribing a zero pore pressure (corresponding to full saturation) at the center of the sample (point A). • This pore pressure is held fixed for 600 seconds to model the fluid acquisition process. • The analysis is performed with the transient soils consolidation procedure using automatic time incrementation. • The tolerance on the pore pressure change permitted (UTOL), is set to a large value since we expect the nonlinearity of the material to restrict the size of the time increments during the transient stages of the analysis. A Pore pressure at 18.88 seconds Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.69

Partially Saturated Example Problems • Since the volume occupied by the sample is fixed (all displacements have been constrained) the volume of fluid absorbed is the same in the cases of the sample with and without gel particles.

Fluid volume absorbed vs. time for samples with and without gel Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.70


F E D C B A

Pore pressure history for samples with and without gel Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Saturation history for samples with and without gel

L5.71

Partially Saturated Example Problems • The difference between the samples with and without the gel appears in the proportions of the volume of the sample that will be occupied by free fluid and fluid trapped in the gel particles.

Void ratio history

Gel volume ratio history for sample with gel


L5.72


Volume of different phases (without gel)

Volume of different phases (with gel) Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.73

Partially Saturated Example Problems • The cyclic wetting analysis of the sample with gel particles contains three steps. • The initial wetting occurs during the first 600 seconds.

F E

D

• This is the same as the original analysis. C

• Zero pore pressure is prescribed at the center of the sample, point A.

B A

• Draining occurs for the next 600 seconds • The draining is modeled by prescribing a pore pressure of −10000 at point A. • This pore pressure corresponds to a saturation of 10%, which is the least saturation the sample can have once it has been wetted. • Finally, the rewetting process occurs over a time period of 800 seconds. • Once again zero pressure is prescribed at point A, corresponding to full saturation at the center of the sample. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.74

Partially Saturated Example Problems • As expected, pore pressure and saturation increase and decrease as the volume of fluid in the sample increases and decreases,

Fluid volume absorbed vs. time Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Pore pressure history

L5.75

Partially Saturated Example Problems • In the rewetting stage, full saturation is achieved more quickly than in the first wetting stage since the sample starts off with a higher saturation and has less capacity to absorb fluid.

Saturation history in cyclic wettability test Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.76

Partially Saturated Example Problems • During the draining stage the gel particles stop growing when the saturation of the surrounding free fluid falls below the value required to keep the gel growing. • At the same time the void ratio to remains constant.

Gel volume ratio history Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Void ratio history

L5.77

Partially Saturated Example Problems • As in any transient problem the element size and the time increment size are related, to the extent that time increments smaller than a certain size give no useful information. • In this example the minimum time increment size is about 0.1 seconds, so an initial time increment of 1 second was used. • If time increments smaller than the critical value are used, spurious oscillations may appear in the solution. Pore pressure history when initial time increment is 0.09 seconds Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.78

Partially Saturated Example Problems • Example: Calculation of phreatic surface in an earth dam (steady-state) • In this example fluid fluid flow occurs in a gravity field and only part of the porous medium is fully saturated, so the location of the phreatic surface is a part of the analysis solution. • Such problems are common in hydrology. • In this example the location of the phreatic surface in an earth dam is identified (Example Problem 9.1.2). • Another example is the well draw-down problem, where the phreatic surface of an aquifer must be located, based on pumping rates at particular well locations. • The basic approach takes advantage of the Abaqus capability to perform partially and fully saturated analysis: • The phreatic surface is located at the boundary of the fully saturated part of the model. • This approach has the advantage that the capillary zone, just above the phreatic surface, is also identified. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.79

Partially Saturated Example Problems • In a typical dam phreatic surface problem the upstream face of the dam (surface S1) is exposed to water in the reservoir behind the dam. • Since Abaqus uses a total pore pressure formulation, the pore pressure on this face must be prescribed to be uw = ( H1 − z ) g ρ w .

Impermeable material

Typical dam phreatic surface problem

• Likewise, on the downstream face of the dam (surface S2 ), uw = ( H 2 − z ) g ρ w . Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.80

Partially Saturated Example Problems • The bottom of the dam (surface S3) is assumed to rest on an impermeable foundation. • This is the natural boundary condition in the pore fluid flow formulation, so no further specification is needed on this surface.


Typical dam phreatic surface problem


L5.81

Partially Saturated Example Problems • The phreatic surface in the dam, S4, is found as the locus of points at which the pore fluid pressure, uw, is zero. • Above this surface the pore fluid pressure is negative, representing capillary tension causing the fluid to rise against the gravitational force and thus creating a capillary zone. • The saturation associated with particular values of capillary pressure is given by the absorption/exsorption curves.


Typical earth dam Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.82

Partially Saturated Example Problems • A special boundary condition is needed if the phreatic surface reaches an open, freely draining surface, as indicated on surface S5. • In such a case the pore fluid can drain freely down the face of the dam, so uw = 0 at all points on this surface below its intersection with the phreatic surface. • Above this point uw< 0, with its particular value depending on the solution.


Typical earth dam Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.83

Partially Saturated Example Problems • The particular earth dam considered here is filled to two-thirds of its height, and only a part of its base is impermeable. • We consider fluid flow only: deformation is ignored. • An analytical solution is available for comparison (Harr, 1962).

H1 =

Earth dam and analytical phreatic surface Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.84

Partially Saturated Example Problems • This example is specifically chosen to illustrate the use of the Abaqus drainage-only flow boundary condition.

uw = ( H1 − z ) g ρw

Drainage-only flow boundary Impervious

Finite element mesh

• On drainage-only edges the pore fluid pressure, uw, is constrained by a penalty method to be less than or equal to zero, thus enforcing the proper drainage-only behavior. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.85

Partially Saturated Example Problems • The weight of the water is applied with a gravity load. • A steady-state soils analysis is performed in five increments to allow Abaqus to resolve the high degree of nonlinearity. • Examining the steady-state contours of pore pressure, we see that the upper-right part of the dam shows negative pore pressures, indicating that it is partially saturated or dry.

Steady-state pore pressure Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.86

Partially Saturated Example Problems • By reducing the number of pore pressure contours to two, one positive and the other negative, the phreatic surface is shown clearly. • This surface compares well with the analytical phreatic surface calculated by Harr (1962).

Phreatic surface Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L5.87

Partially Saturated Example Problems • The contours of saturation show a fully saturated region under the phreatic surface and decreasing saturation in and above the phreatic zone.

Saturation


Modeling Aspects Lecture 6


L6.2

Overview • Element Technology • Geostatic States of Stress • Infinite Domains • Pore Fluid Surface Interactions • Element Addition and Removal • Material Wear/Ablation through Adaptive Meshing • Hydraulic Fracture • Reinforced Soil Slopes • Other Applications


Element Technology


L6.4

Element Technology • The Abaqus geotechnical constitutive models are available for plane strain, generalized plane strain, axisymmetric, axisymmetric-asymmetric, and three dimensional problems. • All Drucker-Prager models are also available in plane stress, except for the linear Drucker-Prager model with creep. • The analysis of flow through porous media in Abaqus is available for plane strain, axisymmetric, axisymmetric-asymmetric, and threedimensional problems. • Special coupled displacement/pore pressure elements must be used: • These elements have a linear distribution of pore pressure and either a linear or quadratic distribution of displacement. • The modified tetrahedral element C3D10MP(H) is particularly well suited for meshing general, complex structures in three dimensions. • The element works well in contact interactions. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.5

Element Technology • Cylindrical (CCL...) elements are available for modeling structures that are initially circular but are subjected to general, nonaxisymmetric loading. • An example is the analysis of a pile foundation where the pile is subjected to axial, horizontal, and moment loading. • These elements permit a coarse yet accurate discretization of a structure. • They can be used in contact calculations using the standard surface-based contact modeling approach. • They provide an attractive alternative to axisymmetric-asymmetric (CAXA...) family of elements. • These elements are not available for the analysis of flow through porous media. • These elements are not currently supported by Abaqus/CAE.

Concrete Slump Test Benchmark 1.1.10


Geostatic States of Stress


L6.7

Geostatic States of Stress • In most geotechnical problems a nonzero state of stress exists in the medium. • This typically consists of a vertical stress increasing linearly with depth, equilibrated by the weight of the material, and horizontal stresses caused by tectonic effects. • The active loading is applied on this initial stress state. • Active loading could be the load on a foundation or the removal of material during an excavation. • Except for purely linear analyses, the response of the system will be different for different initial stress states. • This illustrates a point of nonlinear analysis: • The response of a system to external loading depends on the state of the system when that loading sequence begins (and, by extension, to the sequence of loading). • The linear analysis concept of superposing load cases does not apply. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.8

Geostatic States of Stress • Abaqus provides the geostatic procedure to allow the user to establish the initial stress state. • The user will normally specify the initial effective stresses as model data and, in the first step of analysis, apply the body (gravity) loads corresponding to the weight of the material. • The *INITIAL CONDITIONS, TYPE=STRESS, GEOSTATIC option allows for the definition of the initial stresses due to the geostatic loading. • This option is not currently supported by Abaqus/CAE. • Use the Keywords Editor to include it in your model. • The *GEOSTATIC procedure option is used for this first analysis step in which the geostatic loading is applied. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.9

Geostatic States of Stress • Ideally, the geostatic loads and initial stresses should exactly equilibrate and produce zero deformations. • However, in complex problems it may be difficult to specify initial stresses and loads that exactly equilibrate. • The geostatic procedure is used to reestablish initial equilibrium if the loads and initial stresses specified are not in equilibrium. • It will produce deformations while doing this. • If the deformations produced are significant compared to the deformations caused by subsequent loading, the definition of the initial state should be reexamined.


L6.10

Geostatic States of Stress • In a coupled deformation/flow analysis, the geostatic procedure is equivalent to the steady-state soils analysis procedure. • In these problems it is important to establish initial stress equilibrium as well as steady-state flow conditions. • In fully or partially saturated flow problems, the initial void ratio, as well as the initial pore pressure and the initial effective stress, must be defined. • The initial conditions discussion that follows is based on the total pore pressure formulation. • The magnitude and direction of the gravitational loading are defined by using the GRAV ∗DLOAD option. • In the following discussion the z-axis points vertically upward.


L6.11

Geostatic States of Stress • Assume that, in the geostatic state, the pore fluid is in hydrostatic equilibrium so that

duw = −γ w , dz where γ w is the specific weight of the pore fluid. • If we also take γ w to be independent of z (which is usually the case, since the fluid is almost incompressible), this equation can be integrated:

(

)

u = γ w zw0 − z , 0

where zw is the height of the phreatic surface, at which u = 0 and above which u < 0 and the pore fluid is only partially saturated.


L6.12

Geostatic States of Stress • We usually assume that there are no significant shear stresses τxz, τyz. Then vertical equilibrium gives

dσ zz = ρ g + sn0γ w , dz where ρ is the dry density of the porous solid material, g is the gravitational acceleration, n0 is the initial porosity, and s is the saturation (0.0 ≤ s ≤ 1.0). • The ∗INITIAL CONDITIONS, TYPE=STRESS, GEOSTATIC option defines the initial value of the effective stress, σ , as

σ = σ + suI. • Combining this definition with the equilibrium statement in the z-direction and hydrostatic equilibrium in the pore fluid gives

dσ zz ds 0 ⎛ ⎞ zw − z ⎟ , = ρ g − γ w ⎜ s 1 − n0 − dz dz ⎝ ⎠ using the assumption that γ w is independent of z.

(


)

(

)

L6.13

Geostatic States of Stress • In many cases s is constant. • For example, in fully saturated flow s = 1.0 everywhere. • If we further assume that the initial porosity, n0 , and the dry density of the porous medium, ρ, are also constant, the previous equation is readily integrated to give

(

)

(

)(

)

σ zz = ρ g z − z 0 − γ w s 1 − n0 z − z w0 , where z0 is the position of the surface of the medium. • In more complicated cases where s, n0, and/or ρ vary with height, the equation must be integrated in the vertical direction to define the initial values of σ zz ( z ) .


L6.14

Geostatic States of Stress • In partially saturated cases the initial pore pressure and saturation values must lie on or between the absorption and exsorption curves. • In many geotechnical applications there is also horizontal stress. • If the pore fluid is under hydrostatic equilibrium and τ xz = τ yz = 0, equilibrium in the horizontal directions requires that the horizontal components of effective stress do not vary with horizontal position: i.e., σ h ( z ) only, where σ h is any horizontal component of effective stress.


L6.15

Geostatic States of Stress • The horizontal stress is typically assumed to be a fraction of the vertical stress: • Those fractions are defined in the x- and y-directions with the ∗INITIAL CONDITIONS, TYPE=STRESS, GEOSTATIC option. • If the horizontal stress is nonzero, the boundary conditions on any nonhorizontal edges of the finite element model must be fixed in the horizontal direction or infinite elements must be used so that horizontal equilibrium is maintained.


L6.16

Geostatic States of Stress • Partially Saturated Example: Desaturation of Soil Column (Transient) • This example demonstrates the Abaqus capability to solve coupled fluid flow problems in partially saturated porous media where the effects of gravity are important (Benchmark Problem 1.8.4). • We compare Abaqus results with the experimental work of Liakopoulos (1965). • The Liakopoulos experiment consists of the drainage of water from a vertical column of sand. • A column is filled with sand and instrumented to measure the moisture pressure through its height.

Soil column mesh


L6.17

Geostatic States of Stress • Prior to the start of the experiment the column of sand is saturated: • Water is added continually at the top and allowed to drain freely at the bottom. • The flow is regulated until zero pore pressure readings are obtained throughout the column. • At this point flow is stopped and the experiment starts: • The top of the column is made impermeable and the water is allowed to drain out of the column under gravity. • Pore pressures are measured. • We investigate two cases: 1. An uncoupled flow problem, in which the column is not allowed to deform. 2. A coupled problem, in which vertical deformation of the sand is considered.


L6.18

Geostatic States of Stress • The sand column material definition: ∗ELASTIC 1.3E3, 0. *DENSITY 1.5, *POROUS BULK MODULI , 2.E6 *PERMEABILITY, SPECIFIC=10. 4.5E-6, *PERMEABILITY, TYPE=SATURATION 0.666666,.85 1.,1. *SORPTION -100.,.8 -10.,.85 ... *SORPTION,TYPE=EXSORPTION -100.,.8 -10.,.85 ... Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

-

-

-

Absorption/exsorption curve

L6.19

Geostatic States of Stress • The soil column is meshed with pore pressure plane strain elements (CPE8RP). • These elements are included in a rigid body for the rigid column analysis. • Initial conditions are used to define the initial fully saturated state, the initial pore pressure of zero, and the initial void ratio of 0.4235.

• The deformable model has an additional initial condition to define the initial geostatic stress state. Initial vertical stress (deformable model) Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.20

Geostatic States of Stress • Gravity loading is applied to the column at the start of the analysis history. • In the case of the deforming column, an initial geostatic analysis step is performed to establish the initial equilibrium state. • The initial conditions exactly balance the weight of the fluid and dry material so that no deformation takes place. • Zero pore pressure boundary conditions enforce the initial steady state of fluid flow. *STEP,INC=1 *GEOSTATIC 1.E-6,1.E-6 *DLOAD All_Elements,GRAV,10.,0.,-1.,0. *BOUNDARY POR_Corner_Nodes,8,,0.


L6.21

Geostatic States of Stress • The fluid is allowed to drain through the bottom of the column—by prescribing zero pore pressures at these nodes— during a soils consolidation step. • The fluid will drain until the pressure gradient is equal to the weight of the fluid, at which time equilibrium is established. • The transient analysis is performed using automatic time incrementation. • The pore pressure tolerance that controls the automatic incrementation (UTOL) is set to a large value since we expect the nonlinearity Deformableof the material to restrict the time increment size.

Rigid Deformable

Rigid Final pore pressure


L6.22

Geostatic States of Stress

Deformable column

Rigid column

Pore pressure profiles Numerical vs. experimental results Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.23

Geostatic States of Stress • The pore pressure results of the coupled analysis are closer to the experiment than those of the uncoupled analysis. • The uncoupled analysis overestimates the pore pressures in the early stages of the transient. • This suggests that the coupled analysis is a better approximation of reality (as expected). • At steady state, the pore pressure gradient equals the fluid weight density (as required by Darcy's law).


L6.24

Geostatic States of Stress • As the transient analysis progresses, the material deformation slows and therefore the rigid column assumption becomes closer to reality. • As steady state is approached, both numerical solutions are in good agreement with the experiment. A

F

B C D

E F

E

D C B A

Displacement histories (deformable column) Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.25

Geostatic States of Stress • The time histories of the volume of fluid lost through the bottom of the column are shown for both the deformable and rigid columns. • As expected, more fluid is lost in the deforming column case.

Fluid volume lost (deformable and rigid columns) Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Infinite Domains


L6.27

Infinite Domains • Infinite elements (CIN…) are used in conjunction with the standard elements in problems involving infinite or very large domains. • A family of axisymmetric, planar, and three-dimensional infinite elements is available.

Axisymmetric solid infinite elements

Plane stress and plane strain solid infinite elements


L6.28

Infinite Domains

Three-dimensional solid infinite elements Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.29

Infinite Domains • Standard finite elements are used to model the region of interest, with the infinite elements modeling the far-field region.

• The solution in the far field is assumed to be linear, so only linear behavior is provided in the infinite elements. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.30

Infinite Domains • It is important to make an appropriate choice of the position of the nodes in the infinite direction with respect to the origin (“pole”) of the far-field solution. • For example, the solution for a point load applied to the boundary of a halfspace has its pole at the point of application of the load. • The second node along each infinite element edge pointing in the infinite direction must be positioned so that it is twice as far from the pole as the node on the same edge at the boundary between the finite and the infinite elements.


Pole

L6.31

Infinite Domains • Be careful when specifying the second nodes in the infinite direction so that the element edges in the infinite direction do not cross over.

• Infinite elements are not currently supported by Abaqus/CAE. • If you want to assign these element types to a model, you must use a text editor to add them to the input file generated in the Job module.


L6.32

Infinite Domains • The ∗NCOPY, POLE option provides a convenient way of defining the second nodes in the infinite direction.

Pole node

old set ∗NCOPY, POLE option


new set

L6.33

Infinite Domains • The static behavior of the infinite elements is based on modeling the basic solution variable, u (in stress analysis u is a displacement component), with respect to spatial distance r measured from a “pole” of the solution so that: u → 0 as r → ∞ , and u → ∞ as r → 0. • The interpolation provides terms of order 1 r and 1 r 2 . The far-field behavior of many common cases, such as a point load on a half-space, is thereby included.

Pole

u→∞ r →0

r

u→0 r →∞


L6.34

Infinite Domains • In plane stress and plane strain problems in which the loading is not selfequilibrating, the far-field displacement solution is typically of the form u = ln ( r ) . • This implies the displacement approaches infinity as r → ∞. • Infinite elements can still be used for such cases, provided the displacement results are treated as having an arbitrary reference value. • Thus, strain, stress, and relative displacements within the finite element part of the model will converge to unique values as the model is refined. • The total displacements will depend on the size of the region modeled with finite elements. • If the loading is self-equilibrating, the total displacements will converge on a unique solution.


L6.35

Infinite Domains • In direct integration dynamic response analysis and in direct-solution steady-state harmonic response analysis, infinite elements provide “quiet” boundaries to the finite element model. • This means that they maintain the static force that was present at the start of the dynamic response analysis on the finite/infinite boundary. • As a consequence the far-field nodes in the infinite elements will not displace during the dynamic response (there is no dynamic response within the infinite elements). • Since infinite elements hold the static stress on the boundary constant but do not provide any stiffness, some rigid body motion of the region modeled will generally occur. • This effect is usually small. • The infinite elements will provide additional normal and shear tractions on the boundary, proportional to the normal and shear components of the velocity of the boundary. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.36

Infinite Domains • Damping is applied at the boundary between finite and infinite elements to eliminate the reflection of wave energy back into the finite element mesh when plane waves cross the plane boundary: Finite elements

Infinite elements

Plane waves do not reflect at the boundary

• Infinite elements do not provide any contribution to eigenmode-based analysis procedures.


L6.37

Infinite Domains • For dynamic analyses the mesh should be defined such that transmission of energy from the finite elements to the infinite elements is optimized. • Infinite elements are based on eliminating energy transmission for plane waves crossing a parallel plane boundary. • Therefore, the ability of the elements to transmit energy out of the finite element mesh without trapping or reflecting it is optimized by making the finite/infinite boundary: • as close as possible to being orthogonal to the direction from which the waves will impinge on this boundary • and far enough from the detailed part of the mesh to be considered relatively plane. • Close to a free surface where Rayleigh waves may be important, infinite elements are most effective if they are orthogonal to the free surface. • Similarly, close to a material interface where Love waves may be important, infinite elements are most effective if they are orthogonal to the interface surface. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.38

Infinite Domains • Example: Wave propagation in an infinite medium. • This example shows the effectiveness of the infinite element quiet boundary formulation in a wave propagation problem (Benchmark Problem 2.2.1). • The analysis is of a infinite half-space subjected to vertical pulse distributed load.

Elastic material without damping: E =73 GPa v = 0.33 ρ = 2842 kg/m3


Load Amplitude

L6.39

Infinite Domains • The results using three different plane strain meshes will be compared: C

C

C

B

B

B 2 mm

2 mm

A

A 2 mm

small finite element mesh

A 6 mm

small finite element mesh including infinite element quiet boundaries extended finite element mesh

• The finite element meshes are assumed to have free boundaries at the far field and will reflect the propagating waves. • Results obtained using the small mesh without the infinite element quiet boundaries are given to show how the solution is affected by the reflection of the propagating waves. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.40

Infinite Domains • Based on the material properties, the speed of propagation of push waves in the material is approximately 6169 m/s and the speed of propagation of shear waves is approximately 3107 m/s. • The predominant push waves should reach the boundary of the small mesh in 0.324 µ s and reach the boundary of the extended mesh in 0.97 µ s. The analyses are run for 1.5 µ s, so the waves are allowed to reflect into the meshes. • The following results show that: • Initially, all three meshes give very similar results. • The small finite element mesh response is very different as soon as the waves have had time to reflect.


L6.41

Infinite Domains

C B A

Point A

Point B

Extended Mesh

Extended Mesh

Point C

Quiet Boundary

Extended Mesh

Quiet Boundary

Quiet Boundary

Small mesh Small mesh Small mesh

Vertical displacement (m) at points A, B, and C. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.42

Infinite Domains

Velocity at 0.32 µ s Abaqus/Standard models Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Pore Fluid Surface Interactions


L6.44

Pore Fluid Surface Interactions • The standard contact options offered in Abaqus for stress analysis can be used for geotechnical applications. • Contact theory and usage are discussed further in the "Contact in Abaqus/Standard" lecture notes. • During an analysis of porous media, contact interactions may involve fluid flowing across contact surfaces. • Two configurations are possible: • Pore pressure continuum elements on either side of the interface. • In this case pore pressure continuity between material on opposite sides of an interface must be maintained. • Pore pressure elements on one side and “regular” elements on the other. • This models fluid interaction with an impermeable surface.


L6.45

Pore Fluid Surface Interactions • The flow pattern is shown below:

• It is assumed that no fluid flow takes place tangential to the contact surface. • In steady-state analysis this implies that fluid that flows out of one side flows into the other side. • In transient analysis the flow into the interface is balanced with the rate of separation of the sides of the interface.


L6.46

Pore Fluid Surface Interactions • Pore pressure is continuous across the interface, whether the element is open or closed. • The contact condition is based on effective stress. • Hence, contact points may open as the pore pressure increases. • The total stress is constant and the effective stress decreases due to the pore pressure increase. • In a consolidation analysis, fluid volume between the surfaces is taken into consideration by balancing the flow from each surface. • The fluid in the interface is assumed to be incompressible.


L6.47

Pore Fluid Surface Interactions • Example: Contact between two permeable rings. • Two axisymmetric, linear elastic rings with permeability of 1.0 × 10−4 m/s are in contact.

Pressure: 100 Pa

• The transient analysis is performed in two steps: Step 1: Apply pressure and pore fluid velocities. Step 2: Slide the outer ring in the vertical direction.

Pore fluid velocity: 1×10−4 m/s

• Pore pressure continuity is maintained across the contact interface.

Pore fluid velocity: 1×10−4 m/s

Pore pressure at the end of Step 1.


L6.48

Pore Fluid Surface Interactions • Example: Pile Setup • Pile set-up is the time-dependent increase in pile pull-out resistance of driven piles. • Occurs as the remolded soil surrounding the pile consolidates. • Pore pressure in the soil surrounding the pile gets increased as the pile is driven. • This pore pressure dissipates with time and the effective stress increases. • Increase in the radial effective stress increases the frictional resistance against pile pull-out. • Reliable estimates of pile set-up may help in reducing design pile sections.


L6.49

Pore Fluid Surface Interactions • Pile Set-up analysis: • Pore pressure in the soil surrounding the pile decreases with time.

Pore pressure after 10 seconds

Pore pressure after 6.9 days


L6.50

Pore Fluid Surface Interactions • Pile Set-up analysis: • Radial effective stress in the soil surrounding the pile increases with time.

Radial effective stress after 10 seconds Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Radial effective stress after 6.9 days

L6.51

Pore Fluid Surface Interactions

Radial compressive effective stress increases with time

Frictional pull-out resistance of pile increases with time

Pull-out resistance after 10 seconds

Pull-out resistance after 6.9 days


Element Addition and Removal


L6.53

Element Addition and Removal • Many geotechnical applications involve a sequence of steps; in each of these, some material is removed or added. • Examples: excavation, embankment building, and tunneling. • During the analysis some part of the material mass is removed. • Liners or retaining walls may be inserted during the process.

Excavation of a trench with lining.

• Thus, geotechnical problems offer an interesting perspective on the need for generality in creating and using a finite element model: • The model itself, and not just its response, changes with time. • Parts of the original model disappear, while other components that were not originally present are added. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.54

Element Addition and Removal • The ∗MODEL CHANGE option allows for the removal and addition of elements during the analysis history. • This option is not currently supported by Abaqus/CAE. • Use the Keywords Editor to include it in your model. • To ensure that element removal has a smooth effect on the solution, elements are removed gradually during the removal step. • During the removal step, forces on neighboring elements due to stresses in the removed elements are ramped down to zero. • Since the influence of the removed elements is not completely absent until the end of the step, the removal should be done in a step were all other loads and boundary conditions remain unchanged.


L6.55

Element Addition and Removal • Elements can be reactivated either with or without strain. • Use ∗MODEL CHANGE, ADD=STRAIN FREE to reactivate elements without strain. • For example, when simulating the addition of a new, strain-free layer to a strained construction. • Upon reactivation element variables (stress, strain, plastic strain, etc.) are reset to zero. • Use ∗MODEL CHANGE, ADD=WITH STRAIN to reactivate elements with strain. • For example, when simulating the refueling of a nuclear reactor, where the new fuel assembly must conform to the distortion of its old neighbors. • Upon reactivation the strain is ramped from zero up to the required value; hence, reactivation with strain should be done in its own step.


L6.56

Element Addition and Removal • Loads that are applied in an area where elements are removed or reactivated require special consideration. • Any distributed loads, fluxes, flows, and foundations specified for inactive elements are also inactive. A record of these loads is still kept. • Continuation of loads across steps is not affected by removal. • On reactivation, these loads are still present, unless they are explicitly removed. • Concentrated loads or fluxes are not inactivated when surrounding elements are removed. • You must remove any concentrated loads or fluxes that are carried solely by the elements being removed. • Failure to remove such loads will result in a solver problem (a force is applied to a degree of freedom with zero stiffness).


L6.57

Element Addition and Removal • It may be necessary to reintroduce removed elements with some change to the nodal degrees of freedom. • Nodal variables (displacement, rotation, pore pressures, temperatures, etc., depending on the element type) are not changed by the ∗MODEL CHANGE option. • However, you can reset these variables using a boundary condition while the elements are inactive. • For example, to reintroduce removed elements with a different displacement: 1. Reset the nodal displacements of the removed elements with boundary conditions in an intermediate step. 2. Remove the boundary conditions when the elements are reactivated.


L6.58

Element Addition and Removal • Excavation and Building Analysis Example: Tunneling problem • This analysis models the sequential excavation of a tunnel designed according to the principles of the New Austrian Tunneling Method (NATM). • The excavation is done in steps to sequentially remove the top heading, bench, and invert parts of the tunnel. • Shotcrete is used as the initial liner according to the NATM method. • After excavation, some concrete is cast at the bottom of the tunnel. • This tunneling example was provided by Sauer Corporation.


L6.59

Element Addition and Removal • The plane strain model of the tunnel cross-section to be excavated in layered rock is shown. • A coarse mesh is used for illustration purposes:

Top heading

Shotcrete liners layer 4 layer 3

• 164 CPE4 rock and concrete elements

layer 2

Bench

layer 1

• 15 B21 lining elements • Geostatic initial stress defined for each rock layer:

Invert and Concrete floor

Finite element model Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.60

Element Addition and Removal • The concrete and the shotcrete are modeled as linear elastic. • The rock is modeled with the modified Drucker-Prager model. • The elastic properties are made dependent on a predefined field variable, so that the stiffness can be degraded during the excavation process. *MATERIAL, NAME=MBENCH *ELASTIC, TYPE=ISOTROPIC 2414., .27, 0. 2414., .27, 1. Bench material 965.6, .27, 2. Elastic dependence on 965.6, .27, 8. *DRUCKER PRAGER temperature field 54.82, 1, 54.82 *DRUCKER PRAGER HARDENING .464 *DENSITY .0247


L6.61

Element Addition and Removal • Step 1: Remove the lining elements and the concrete bottom of the tunnel. • These elements must be included in the original mesh, but they are not in place at the start of the operation.

Gravity Top heading Bench Invert

• Apply the gravity loads that equilibrate the virgin stress state in the rock—this should cause no deformation.

Concrete floor and liners removed


L6.62

Element Addition and Removal

• Step 2: Degrade the stiffness of the top heading of the tunnel. *FIELD ALNOD, 1.

Top heading degraded

Gravity

Bench Invert

Deformed model shape after Step 2 scale factor 300


L6.63

Element Addition and Removal • Step 3: Activate (strain free) the shotcrete lining elements for removal of top heading. Top heading lining

Gravity

Bench

• Step 4: Excavate top heading of the tunnel by removing the corresponding elements.

Invert

Deformed model shape after Step 4 scale factor 300


L6.64

Element Addition and Removal • Steps 5–10: Repeat Steps 2-4 twice for the excavation of the bench and invert parts of the tunnel and the addition of their respective lining.

Deformed shape scale factor 300

Displacement after bench excavation/lining Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Displacement after invert excavation/lining

L6.65

Element Addition and Removal • Step 11: Concrete is cast in place at the bottom of the tunnel.

Deformed shape scale factor 300

Concrete floor

Final displacement after concrete casting Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.66

Element Addition and Removal • Analysis of an earthen dam subjected to rapid drawdown and earthquake excitation

• Length of base = 245.2 m • Height of dam = 55 m from ground Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.67

Element Addition and Removal • Stages in the construction of the dam:

First stage

Second stage

Third stage

Final stage


L6.68

Element Addition and Removal • Pore pressure distribution during the reservoir fill-up stages:

First stage of the fill-up

Second stage of the fill-up

Third stage of the fill-up

After steady-state is reached


L6.69

Element Addition and Removal • Pore pressure distribution after rapid drawdown:


L6.70

Element Addition and Removal • Response to an earthquake: • Implicit dynamic analysis performed using *DYNAMIC • Inertial effects due to pore fluid modeled using lumped mass user elements placed at the nodes • Mass values computed using user subroutine UVARM • Potential for liquefaction assessed using • Cyclic stress ratio (CSR) and Cyclic resistance ratio (CRR)

Red contours indicate higher likelihood of liquefaction Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Material Wear/Ablation through Adaptive Meshing


L6.72

Material Wear/Ablation through Adaptive Meshing • The material wear/ablation feature in Abaqus enables modeling of material wear/erosion on the surface of a body. • Surface erosion is implemented using the adaptive meshing technique, enabling wear to extend many element lengths. • The wear surface defines the external boundary of an adaptive mesh domain.

Fluid velocity dependent wear of a well bore


L6.73

Material Wear/Ablation through Adaptive Meshing • The adaptive meshing is performed as a post-converged-increment operation.

Lagrangian increment Converged

Output

Mesh smoothing Advection


L6.74

Material Wear/Ablation through Adaptive Meshing • Adaptive mesh constraints define wear extent or velocity. • General descriptions of wear criteria possible through the user subroutine UMESHMOTION. • UMESHMOTION provides access to solution variables: • Nodal • Material • Contact • A surface-normal local coordinate system is also provided. • Total volume lost due to wear/erosion is available through the output variable, VOLC.


L6.75

Material Wear/Ablation through Adaptive Meshing • Usage • Adaptive mesh options *Adaptive mesh, elset= *Adaptive mesh constraint, type=[velocity|displacement], User *Adaptive mesh controls


L6.76

Material Wear/Ablation through Adaptive Meshing • Usage (cont'd) • Through adaptive mesh constraints the user describes a spatial mesh constraint velocity or displacement. • Adaptive mesh constraints are not supported in Abaqus/CAE, but can be defined using the Keywords Editor.


L6.77

Material Wear/Ablation through Adaptive Meshing • User subroutine interface: subroutine umeshmotion(uref,ulocal,node,nndof,lnodetype,alocal, * ndim,time,dtime,pnewdt,kstep,kinc,kmeshsweep,jmatyp,jgvblock,lsmooth) c include 'aba_param.inc' c dimension ulocal(ndim), jelemlist(*) dimension alocal(ndim,*), time(2) dimension jmatyp(*), jgvblock(*) user coding to define ULOCAL and, optionally, PNEWDT return end

• Key subroutine arguments: • Local surface directions are provided. • Current mesh motion must be defined in a surface local coordinate system. • The ratio of the suggested new time increment to the time increment currently being used can be updated. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.78

Material Wear/Ablation through Adaptive Meshing • The following utility routines are also provided for access to local solution results: GETVRN, e.g. for local pore pressure GETVRMAVGATNODE, e.g. for local PEEQ or FLVEL GETNDTOELMCONN, for access to elements adjacent to nodes


L6.79

Material Wear/Ablation through Adaptive Meshing • Example: Erosion of rock in a wellbore • Abaqus Example Problem 1.1.21: Erosion of material (sand production) from oil bore hole perforation tunnel • Erosion occurs when oil is extracted under a sufficiently high pressure gradient. • Steps performed: • Geostatic • Model change removal of well bore and casing • Apply pore pressure boundary conditions at the nodes on the bore well surface

Perforation tunnel

• Reduce pore pressure to the extraction pore pressure level and perform a steady state analysis

Bore hole

• Consolidation with material erosion for a time period of four days Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.80

Material Wear/Ablation through Adaptive Meshing • Results: • It is seen that the largest amount of material is eroded near the junction of the bore hole and the perforation tunnel. • Further away from the bore hole boundary the amount of erosion progressively decreases. • This behavior is expected because • there are high strains near the junction of the bore hole and the perforation tunnel, and • the erosion criterion is active only for values of the equivalent plastic strain above a threshold value. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

Equivalent plastic strain at the end of four days

L6.81

Material Wear/Ablation through Adaptive Meshing • Results: • Total volume of sand produced (eroded) in cubic inches as a function of analysis time in days:


Hydraulic Fracture


L6.83

Hydraulic Fracture • Hydraulic fractures are created by pumping fluids into oil wells at very high pressures. • Surface area exposed to the hydrocarbon bearing rocks gets increased. • Conductive pathways get created for enhanced hydrocarbon flow to the well bore. • Oil-well productivity increases. • Production lifetime of reservoirs gets extended. Hydraulic fracture in a hydrocarbon-bearing strata


L6.84

Hydraulic Fracture • Special-purpose cohesive elements with pore pressure degrees of freedom: • COH2D4P

6-node two-dimensional pore pressure cohesive element

• COHAX4P

6-node axisymmetric pore pressure cohesive element

• COH3D6P

9-node three-dimensional pore pressure cohesive element

• COH3D8P

12-node three-dimensional pore pressure cohesive element

• Fluid pressure on the cohesive element contributes to its mechanical behavior. • Fluid flow continuity within the fracture gap and through the interface is maintained. • An additional layer of flow resistance on the cohesive element surface can be defined.


L6.85

Hydraulic Fracture • Tangential flow within the fracture gap can be specified as: • Newtonian flow • Power law model • Normal flow across the fracture gap can be specified. • It represents resistance due to caking or fouling effects.


L6.86

Hydraulic Fracture • Tangential flow is absent by default. • Newtonian tangential flow: • Volume flow rate density vector is expressed as: q = −kt∇p, where kt is the tangential permeability and ∇p is the pressure gradient along the element. • Tangential permeability is given by Reynold’s equation: kt = d3/12µ, where µ is the fluid viscosity and d is the cohesive element opening. • An upper limit on the value of kt can also be specified. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.87

Hydraulic Fracture • Power law tangential flow: • Constitutive relation is defined as:τ = K γα

τ is the shear stress,

where:

γ is the shear strain rate, K is the fluid consistency, and

α is the power law coefficient. • Tangential volume flow rate density is defined as: 1

⎛ 2α ⎞⎛ 1 ⎞ α ⎛ d ⎞ q = −⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎝ 1 + 2α ⎠⎝ K ⎠ ⎝ 2 ⎠

1+ 2α

α

∇p

1−α

α

∇p

where d is the cohesive element opening. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.88

Hydraulic Fracture • Normal flow across gap surfaces: • Occurs when fluid leak-off coefficient is defined for the pore fluid. • The coefficient defines a pressure-flow relationship between the cohesive element’s middle nodes and their adjacent surface nodes. • The coefficient can be interpreted as the permeability of a finite layer of material on the cohesive element surfaces. • Normal flows are defined as qt = ct (pi− pt) and qb = cb (pi− pb), where qt and qb are the flow rates into the top and bottom surfaces, pb is the bottom surface pore pressure, pi is the mid-face pore pressure, and pt is the top surface pore pressure.


L6.89

Hydraulic Fracture • Usage: Keywords interface • Tangential flow: *GAP FLOW, TYPE=NEWTONIAN *GAP FLOW, TYPE=POWER LAW

• Normal flow: *FLUID LEAKOFF *FLUID LEAKOFF, USER

• Requires user subroutine UFLUIDLEAKOFF • Initially open elements are defined by: *INITIAL CONDITIONS, TYPE=INITIAL GAP


L6.90

Hydraulic Fracture • Usage: Abaqus/CAE interface


L6.91

Hydraulic Fracture • Example • Abaqus Example Problem 9.1.5 • Model dimensions: • 400 m diameter • 50 m thickness • Target formation is sandwiched between two shale layers. • Rock is modeled by C3D8RP elements. • Fracture is modeled by COH3D8P pore pressure cohesive elements. • Bore hole casing is modeled by M3D4 elements. • Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.92

Hydraulic Fracture • Constitutive properties: • A linear Drucker-Prager model with hardening is used for the rock. • A quadratic traction-interaction failure criterion is used for modeling damage initiation in the cohesive elements. • A mixed-mode energy-based damage evolution law is used for modeling damage propagation. • Gap flow is specified as Newtonian with a viscosity of 1 × 10–6 kPa s (1 centipoise), close to that of water. • Fluid leak-off is specified as 5.879 × 10–10 m3/(kPa s) for the early stages. In the final stage, when the polymer is dissolved, the fluid leak-off coefficient is increased to 1 × 10–3 m3/(kPa s). • This step-dependent fluid leak-off coefficient is set in user subroutine UFLUIDLEAKOFF.


L6.93

Hydraulic Fracture • Initial pore pressure, stress, and void ratio distributions are specified. • The analysis consists of 4 steps: • Geostatic step • Consolidation step wherein fluid is injected into the well under pressure • The cohesive elements open up due to the fluid pressure and create hydraulic fractures. • Consolidation step to allow for the built-up pore pressure to dissipate • The fractures are prevented from closing by additional boundary conditions in order to simulate effects due to proppant material. • Consolidation step with a drawdown pressure applied to the well bore nodes • Simulates production from the well Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.94

Hydraulic Fracture • Results:

Fracture geometry following the injection stage (deformations 500x)


History of the fracture opening profile

L6.95

Hydraulic Fracture • Results (cont'd):

History of the fracture pore pressure profile

Well yield before and after the hydraulic fracture process


Reinforced Soil Slopes


L6.97

Reinforced Soil Slopes • The assessment of the strength of soil slopes and investigating the means for increasing their safety against failure are crucial in construction projects involving large soil masses. • Slope stability analyses have traditionally been performed using an analytical limit state approach. • The soil is assumed to fail along a curved surface • At each point along this curved surface the soil is assumed to have fully mobilized its maximum strength in shear • The likelihood of failure of the slope is expressed as a factor of safety for that particular slope • However, any presence of reinforcement or local heterogeneity necessitates the use of numerical techniques such as FEA. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.98

Reinforced Soil Slopes • Finite element analysis: • Advantages: • Suitable for inhomogeneous and reinforced soils • Can include a wide range of plasticity material models • Difficulties: • Convergence problems in an implicit code • A limit state is difficult to reach • Outstanding issues: • Strength reduction technique • What indicates failure?


L6.99

Reinforced Soil Slopes • Key Abaqus features suitable for modeling reinforced soil slopes: • Supports a variety of constitutive models for soils, including Mohr-Coulomb, Drucker-Prager, Cam-Clay, and Cap plasticity models • Can model contact and frictional interaction between soil and structures • Embedded elements can be used to model reinforcement structures such as nails and ties • Ability to model transient distribution of pore fluid pressure in soils in Abaqus/Standard • Ability to model consolidation, or the deformation in soils associated with the redistribution of pore water and pressure, in Abaqus/Standard • In Abaqus/Explicit, ability to model the damage and removal of failed elements, as well as the movement of dislodged regions of soil relative to the intact regions


L6.100

Reinforced Soil Slopes • Soil slope example: • Begin by excavating from a horizontal ground level • The soil in the colored regions is removed in sequential stages • Final slope is inclined at 55 degrees • Vertical height of 20 m • At the height of 20m, the inclination reduces to 10 degrees. • The problem is defined as a two dimensional plane strain analysis. • Linear Drucker-Prager model • Unreinforced slope is just stable based on the limit state approach


L6.101


Nail configurations used in the model Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.102


Configurations after successive excavations Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.103


Contour Plots of Equivalent Plastic Strain (PEEQ) Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

L6.104

Reinforced Soil Slopes • Using a Python script a local measure of “safety” was computed based on how far away the current state of stress is from the yield stress at the current pressure:

Mises stress

a DP Yield Factor = b

Pressure


a a+b

L6.105

Reinforced Soil Slopes Unreinforced case:

Contour Plot of PEEQ

Contour Plot of DP Yield Factor


L6.106

Reinforced Soil Slopes Reinforced case – 3 nails:




L6.107

Reinforced Soil Slopes Reinforced case – 7 nails:




L6.108

Reinforced Soil Slopes • Failure modeling using Abaqus/Explicit • Material damage • Automatic deletion of damaged elements • New contact surfaces created along the failure surface • General contact used to model the sliding of the dislodged soil mass


L6.109


This analysis is described in detail in the following Technical Brief: “Analysis of Reinforced and Un-reinforced Soil Slopes using Abaqus”


Other Applications


L6.111

Other Applications • Other applications where Abaqus geotechnical features can be used: • Powder metallurgy • Porous elasticity • Soil plasticity models with hardening, in addition to the porous metal plasticity model • Soft tissues and tissue-implant interaction • Fluid flow and resulting modifications in tissue properties • Coupled consolidation and contact analysis • Absorption of fluids by polymeric materials • Moisture swelling and gel swelling • Coupled consolidation analysis


Stress Equilibrium and Fluid Continuity Equations Appendix 1


A1.2

• To study the different couplings and nonlinearities of the coupled problem, we can write the resulting Newton equations at a node as

where dc represents the vector of displacement corrections, Fr are the force residuals conjugate to the displacements, uc is the pore pressure correction, and ∆Vr is the residual change in fluid volume over the time increment conjugate to the pore pressure. • Special cases of the fully and partly saturated problems are presented in the following sections.


A1.3

• Fully saturated fluid flow • Kdd, Kdu, Kud, Kuu have the following components:

K dd = K s (σ , d ) + L ( d , u ) + K gd ( d )

∫

K s = β : D : β dV V

∫

L = − β : I u I : β dV V

∫

K gd = − N : g f1 I : β dV V

K du = B ( d ) + K sg (σ , d )

∫

B = − β : I dV V


A1.4

∫ 3K 1

K sg =

V

β : D : I dV g

T K ud = BT ( d ) + K sg (σ , d ) + ∆tLc ( d , e ) + ∆tkec (σ , d , e )

∫

BT = − I : β dV V

∫ 3K 1

T = K sg

V

Lc = −

∂δ u

∫ ∂x

V

kec =

I : D : β dV g

∫

V

⎛ ∂u ⎞ ⋅ k * ⋅ ⎜ − ρ w g ⎟ I : β dV ⎝ ∂x ⎠

∂δ u dk * ⎛ ∂u ⎞ 1 ⋅ ⋅ ⎜ − ρw g ⎟ ∂x de ⎝ ∂x ⎠ (1 − n )2


⎡ I : D 1 − n0 − ⎢ JK g ⎢⎣ 3K g

⎤ I ⎥ : β dV ⎥⎦

A1.5 ∗ Kuu = ∆tk ( d , e ) + K sg (σ ) + K g∗ ( d ) + K w∗ ( d ) + ∆t ke (σ , d , e )

k =−

∂δ u

∫ ∂x

⋅ k∗ ⋅

V ∗ K sg

∫ 9K 1

=

V

K g∗ = −

∫

V

K w∗ =

∫

V

ke =

∫

V

2 g

∂du dV ∂x

I : D : I dV

1 − n0 dV J Kg

1 − n0 1 − dV J Kw Kw

∂δ u dk ∗ ⎛ ∂u ⎞ 1 ⋅ ⋅ ⎜ − ρw g ⎟ ∂x de ⎝ ∂x ⎠ (1 − n )2

⎡ I : D : I 1 − n0 ⎤ − ⎢ ⎥ dV 2 J K g ⎦⎥ ⎣⎢ 9 K g

where β is the strain-displacement matrix,D is the constitutive matrix, N is the interpolator, Kg and Kw are the solid grains and fluid bulk moduli. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

A1.6

• Special cases • To facilitate the understanding of the coupled equations, let us consider some special cases of the transient problem without fluid gravity effects: • Linear material, small strain, incompressible grains and fluid, constant permeability:

where Ks is the usual stress stiffness; B, BT are the stress/pore pressure coupling terms; and ∆t k is the porous medium permeability term. The resulting system of equations is linear and symmetric. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

A1.7

• Nonlinear material, large strain, incompressible grains and fluid, nonlinear permeability:

where K s (σ , d ) is the stress stiffness term with material and geometric nonlinearities; L (d, u) is a large volume change term; B (d) , BT (d) are the stress/pore pressure coupling terms with geometric nonlinearity; ∆t Lc (d, e) is a large-strain coupling term with nonlinear permeability; and ∆t k (d, e) is the permeability term with geometric and permeability nonlinearities. The problem now becomes nonlinear and unsymmetric. The loss of symmetry is due to the inclusion of finite strains. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

A1.8

• Linear material, small strain, compressible grains and fluid, constant permeability:

where Ks is the usual stress stiffness; B, BT are the stress/pore pressure coupling terms; Ksg, KTsg, K*sg, K*g are grain compressibility terms; K*w is a fluid compressibility term; and ∆t k is the porous medium permeability term. The resulting system of equations is linear and symmetric.


A1.9

• Nonlinear material, large strain, compressible grains and fluid, nonlinear permeability:

where ∆t kec (σ , d , e ) , ∆t ke (σ , d , e ) are additional nonlinear grain compressibility terms with large-strain effects. This most general problem is nonlinear and unsymmetric. The loss of symmetry is now due to two reasons: finite strains and nonlinear permeability.


A1.10

• Consider now the steady-state problem without fluid gravity effects. In its simplest case (linear material, small strain, incompressible grains and fluid, constant permeability), it is only one way coupled:

• The equations are clearly unsymmetric. • The steady-state problem becomes fully coupled if large strains and/or nonlinear permeability are considered. However, the equations retain their unsymmetric characteristic for these cases. • In Abaqus the steady-state equations are always solved directly in the coupled form for all cases. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

A1.11

• The porous media coupled analysis capability can provide solutions either in terms of total or of excess pore pressure. The excess pore pressure at a point is the pore pressure in excess of the hydrostatic pressure required to support the weight of pore fluid above the elevation of the material point. • [Total pore pressure solutions are provided when the GRAV distributed load type is used to define the gravity load on the model. Excess pore pressure solutions are provided in all other cases (for example, when gravity loading is defined with distributed load types BX, BY, or BZ)]. • One important aspect arising from the inclusion of pore fluid gravity (in total pore pressure analyses) is that it generates a Newton method Jacobian term that is always unsymmetric:

∫

K gd = − N : g f1 I : β dV .( K dd ) V


A1.12

• Partially saturated fluid flow

• Kdd, Kdu, Kud, Kuu have the following components: K dd = K s (σ , d ) + L ( d , u ) + K gd ( d )

∫

K s = β : D : β dV V

∫

L = − β : I su I : β dV V

∫

K gd = − N : g sf1 I : β dV V

K du = B ( d ) + Bs ( d , u ) + K sg (σ , d ) + K sgs (σ , d , u ) + K gu ( d )

∫

B = − s β : I dV V


A1.13

∫

Bs = − u V

ds β : I dV du

∫ 3K s

K sg =

V

K sgs =

β : D : I dV g

∫ 3K

u ds β : D : I dV g du

V

∫

K gu = − N : g f 2 V

ds dV du

T K ud = BT ( d ) + K sg (σ , d ) + ∆tLc ( d , e ) + ∆tkec (σ , d , e )

∫

BT = − s I : β dV V


A1.14

∫ 3K s

T K sg =

V

∫

Lc = − s V

kec =

∫

V

I : D : β dV g

∂δ u ∗ ⎛ ∂u ⎞ ⋅ k ⋅ ⎜ − ρ w g ⎟ I : β dV ∂x ⎝ ∂x ⎠

∂δ u dk ∗ ⎛ ∂u ⎞ s ⋅ ⋅ ⎜ − ρ wg ⎟ ∂x de ⎝ ∂x ⎠ (1 − n )2

⎡ I : D 1 − n0 ⎤ − I ⎥ : β dV ⎢ ⎣⎢ 3K g JK g ⎥⎦

∗ Kuu = ∆tk ( d , e ) + K sg (σ ) + K g∗ ( d ) + K w∗ ( d ) + K s∗ ( d ) + ∆t ke (σ , d , e ) + ∆t kes ( d , e )

∫

k =− s V

∗ K sg =

∂δ u ∗ ∂du dV ⋅k ⋅ ∂x ∂x

∫ 9K s

V

2 g

I : D : I dV


A1.15

∫

K g∗ = − s V

1 − n0 dV J Kg

⎛ 1 − n0 1 ⎞ K w∗ = s ⎜ − ⎟ dV J Kw Kw ⎠ V ⎝

∫

K s∗ = −

∫

V

ke =

∫

V

ds ⎛ 1 − n0 ⎞ ⎜1 − ⎟ dV du ⎝ J ⎠

∂δ u dk * ⎛ ∂u ⎞ s ⋅ ⋅ ⎜ − ρ wg ⎟ ∂x de ⎝ ∂x ⎠ (1 − n )2

kes = −

⎡ I : D : I 1 − n0 ⎤ − ⎢ ⎥ dV 2 J K g ⎥⎦ ⎢⎣ 9 K g

s ds ∂δ u ∗ ⎛ ∂u ⎞ ⋅ k ⋅ ⎜ − ρ w g ⎟ dV x x ∂ ∂ du ⎝ ⎠ s

∫k

V

where β is the strain-displacement matrix, D is the constitutive matrix, N is the interpolator, Kg and Kw are the solid grains and fluid bulk moduli. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

A1.16

• Special cases • To facilitate the understanding of the coupled equations, let us consider some special cases of the transient problem without fluid gravity effects: • Linear material, small strain, incompressible grains and fluid, constant permeability:

where Ks is the usual stress stiffness; B, BT are stress/pore pressure coupling terms; Bs is a partially saturated coupling term; ∆t k is the permeability term; and ∆t kes, Ks* are partially saturated permeability terms. • The resulting system of equations is unsymmetric. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

A1.17

• Nonlinear material, large strain, incompressible grains and fluid, nonlinear permeability:

where K s (σ , d ) is the stress stiffness term with material and geometric nonlinearities; L(d,u) is a large volume change term; B(d), BT(d) are coupling terms with geometric nonlinearity; Bs(d,u) is a partially saturated coupling term with geometric nonlinearity; ∆t Lc (d, e) is a large-strain coupling term with nonlinear permeability; ∆t k (d, e) is the permeability term with geometric and permeability nonlinearities; and ∆t kes (d, e) , Ks* (d) are partially saturated permeability terms with geometric and permeability nonlinearities. • The problem remains unsymmetric. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

A1.18

• Linear material, small strain, compressible grains and fluid, constant permeability:

where Ks is the usual stress stiffness; B, BT are stress/pore pressure coupling terms; Bs is a partially saturated coupling term; Ksg, KTsg, Ksgs, K*sg, K*g are grain compressibility terms; K*w is a fluid compressibility term; ∆t k is the porous medium permeability term; and ∆t kes , Ks* are partially saturated permeability terms. • The resulting system of equations is unsymmetric. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

A1.19

• Nonlinear material, large strain, compressible grains and fluid, nonlinear permeability:

where ∆t kec (σ , d , e ) , ∆t ke (σ , d , e ) are additional nonlinear grain compressibility terms with large-strain effects. • This most general problem is unsymmetric. • In conclusion, the transient partially saturated flow problem is always unsymmetric. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

A1.20

• Consider now the steady-state problem without fluid gravity effects. In its simplest case (linear material, small strain, incompressible grains and fluid, constant permeability) it is only one-way coupled:

• The equations are unsymmetric. • The steady-state problem becomes fully coupled if large strains and/or nonlinear permeability are considered. The equations retain their unsymmetric characteristic for these cases. • In Abaqus the steady-state equations are always solved directly in the coupled form for all cases. Analysis of Geotechnical Problems with Abaqus © Dassault Systèmes, 2008

A1.21

• The porous media coupled analysis capability can provide solutions either in terms of total or of excess pore pressure. • [Total pore pressure solutions are provided when the GRAV distributed load type is used to define the gravity load on the model. Excess pore pressure solutions are provided in all other cases (for example, when gravity loading is defined with distributed load types BX, BY, or BZ)]. • One important aspect arising from the inclusion of pore fluid gravity (in total pore pressure analyses) is that it generates Newton method Jacobian terms that are always unsymmetric:

∫

K gd = − N : g sf1 I : β dV in K dd V

∫

K gu = − N : g f 2 V


ds dV in K du du

Bibliography of Geotechnical Example Problems Appendix 2


A2.2

References The following is a list of Abaqus Example and Benchmark Problems that show the use of capabilities for geotechnical modeling: • Example problems: 1.1.6:

Jointed rock slope stability

1.1.11:

Stress-free element reactivation

1.1.22:

Erosion of material (sand production) from oil bore hole perforation tunnel

9.1.1:

Plane strain consolidation

9.1.2:

Calculation of phreatic surface in an earth dam

9.1.3:

Axisymmetric simulation of an oil well

9.1.4

Analysis of a pipeline buried in soil

9.1.5

Hydraulically induced fracture in a well bore


A2.3

References • Benchmark problems: 1.1.10:

Concrete slump test

1.8.1:

Partially saturated flow in porous media

1.8.2:

Demand wettability of a porous medium: coupled analysis

1.8.3:

Wicking in a partially saturated porous medium

1.8.4:

Desaturation in a column of porous material

1.14.1:

The Terzaghi consolidation problem

1.14.2:

Consolidation of triaxial test specimen

1.14.3:

Finite-strain consolidation of a two-dimensional solid

1.14.4:

Limit load calculations with granular materials

1.14.5:

Finite deformation of an elastic-plastic granular material


A2.4

References 2.2.1:

Wave propagation in an infinite medium

2.2.2:

Infinite elements: the Boussinesq and Flamant problems

2.2.3:

Infinite elements: circular load on half-space

2.2.4:

Spherical cavity in an infinite medium

3.2.4:

Triaxial tests on a saturated clay

3.2.5:

Uniaxial tests on jointed material


Workshop Preliminaries Setting up the workshop directories and files If you are taking a public seminar, the steps in the following section have already been done for you: skip to Basic Operating System Commands, (p. WP.2). If everyone in your group is familiar with the operating system, skip directly to the workshops. The workshop files are included on the Abaqus release CD. If you have problems finding the files or setting up the directories, ask your systems manager for help. Note for systems managers: If you are setting up these directories and files for someone else, please make sure that there are appropriate privileges on the directories and files so that the user can write to the files and create new files in the directories. Workshop file setup (Note: UNIX is case-sensitive. Therefore, lowercase and uppercase letters must be typed as they are shown or listed.) 1. Find out where the Abaqus release is installed by typing UNIX and Windows NT: abqxxx whereami where abqxxx is the name of the Abaqus execution procedure on your system. It can be defined to have a different name. For example, the command for the 6.8–1 version might be aliased to abq681. This command will give the full path to the directory where Abaqus is installed, referred to here as abaqus_dir. 2. Extract all the workshop files from the course tar file by typing UNIX: abqxxx perl abaqus_dir/samples/course_setup.pl Windows NT: abqxxx perl abaqus_dir\samples\course_setup.pl Note that if you have Perl and the compilers already installed on your machine, you may simply type: UNIX: abaqus_dir/samples/course_setup.pl Windows NT: abaqus_dir\samples\course_setup.pl 3. The script will install the files into the current working directory. You will be asked to verify this and to choose which files you wish to install. Choose “y” for the appropriate lecture series when prompted. Once you have selected the lecture series, type “q” to skip the remaining lectures and to proceed with the installation of the chosen workshops.


Preliminaries for Abaqus Workshops

WP.2

Basic operating system commands (You can skip this section and go directly to the workshops if everyone in your group is familiar with the operating system.) Note: The following commands are limited to those necessary for doing the workshop exercises. Working with directories 1. Start in the current working directory. List the directory contents by typing UNIX: ls Windows NT:

dir

Both subdirectories and files will be listed. On some systems the file type (directory, executable, etc.) will be indicated by a symbol. 2. Change directories to a workshop subdirectory by typing Both UNIX and Windows NT: cd dir_name 3. To list with a long format showing sizes, dates, and file, type UNIX: ls -l Windows NT:

dir

4. Return to your home directory: UNIX: cd Windows NT: cd home-dir List the directory contents to verify that you are back in your home directory. 5. Change to the workshop subdirectory again. 6. The * is a wildcard character and can be used to do a partial listing. For example, list only Abaqus input files by typing UNIX: ls *.inp Windows NT: dir *.inp Working with files Use one of these files, filename.inp, to perform the following tasks: 1. Copy filename.inp to a file with the name newcopy.inp by typing UNIX:

cp filename.inp newcopy.inp

Windows NT: copy filename.inp newcopy.inp 2. Rename (or move) this new file to newname.inp by typing UNIX: mv newcopy.inp newname.inp Windows NT: rename newcopy.inp newname.inp (Be careful when using cp and mv since UNIX will overwrite existing files without warning.)



WP.3

3. Delete this file by typing UNIX: rm newname.inp Windows NT: erase newname.inp 4. View the contents of the files filename.inp by typing UNIX:

more filename.inp

Windows NT:

type filename.inp | more

This step will scroll through the file one page at a time. Now you are ready to start the workshops.



Workshop 1 Material Models for Geotechnical Applications Interactive Version Note: This workshop provides instructions on using geotechnical material models in terms of the Abaqus GUI interface. If you wish to use the Abaqus Keywords interface instead, please see the “Keywords” version of these instructions. Please complete either the Keywords or Interactive version of this workshop.

Introduction When you complete this workshop you will be able to: • • •

Create stress state diagrams. Define a number of geotechnical material models. Use Abaqus/CAE to build simple models and view and interpret material behavior.

Figure W1–1 Configuration of the model cube.



W1.2

In this workshop you will become familiar with some of the geotechnical material models in Abaqus. You will study the behavior of the 2x2x2 m cube shown in Figure W1-1. This cube will be subjected to uniaxial tension, uniaxial compression, pure shear, hydrostatic tension, and hydrostatic compression. Different material models will be assigned to the cube and the response of the cube for these material models will be simulated.

Preliminaries 1. Enter the working directory for this workshop: ../geotechnical/interactive/geomaterials

2. Run the script ws_geotech_geomaterials.py using the following command: abaqus cae startup=ws_geotech_geomaterials.py

The above command creates an Abaqus/CAE database named cube.cae in the current directory. This database contains a complete model named cube. The cube material is linear elastic with a Young’s modulus of 70 GPa and a Poisson’s ratio of 0.3. The faces of the cube lying on the 1-2, 2-3, and 3-1 planes are constrained against displacements normal to these respective planes. Prescribed displacements of 0.01 along the outward normals to the other three faces are specified so that the cube is subjected to a state of hydrostatic tension. This cube is meshed with one reduced integrated element as shown in Figure W1-1. For each of the problem cases in this workshop you will copy the original cube model and modify the applied boundary conditions and material data.

Case 1: Linear Elasticity You will begin this workshop by running the model created by the setup script and studying the behavior of the cube under hydrostatic tension. You will then create new models by copying the original model and modifying the boundary conditions to subject the cube to (a) uniaxial tension, (b) uniaxial compression, (c) hydrostatic compression, and (d) pure shear. After the boundary conditions are modified for a model, you will run the analysis, check for any modeling errors, and then view the results. 1. In an Abaqus/CAE session, open the database named cube.cae, which was created by running the setup script. 2. In the Model Tree, click mouse button 3 on the job named Case_1_hydro_TN and select Submit in the menu that appears to submit the job for analysis. Tip: If the Model Tree is not visible, make sure there is a check mark next to Show Model Tree in the View menu. If the Model Tree is still not visible, drag the cursor from the left side of the Abaqus/CAE window to expand the Model Tree. Once the Model Tree is visible, expand the Jobs container by clicking the “+” symbol next to Jobs.



W1.3

3. While the job is running, click mouse button 3 on the job Case_1_hydro_TN in the Model Tree and select Monitor from the menu that appears to view the progress of the analysis. Check for errors and warning messages under the Errors and Warnings tabs in the job monitor dialog box. 4. After the job has completed, click mouse button 3 on the job Case_1_hydro_TN in the Model Tree and select Results from the menu that appears. 5. Study the analysis results in the Visualization module. Look at the deformed shape and plot contours of the nonzero components of stress (S11, S22, and S33). 6. In the Model Tree, click mouse button 3 on the model cube and select Copy Model from the menu that appears. Name the new model Case_1_uni_TN. 7. Modify the boundary conditions to subject the cube to uniaxial tension in the 1direction. a. In the Model Tree, click mouse button 3 on the BCs container for the model Case_1_uni_TN and select Manager from the menu that appears. b. In Step-1, select BC-5, which is the prescribed displacement of 0.01 in the 2-direction, and click Delete. This will remove this displacement boundary condition. c. Similarly delete BC-6, which refers to the prescribed boundary displacement in the 3-direction. Removing these boundary conditions will allow free expansion or contraction of the cube along the 2- and 3-directions. Do not delete BC-4, which will apply the tension in the 1-direction. This model is now ready for execution. 8. In the Model Tree, double-click on Jobs to create a new job. 9. Name the new job Case_1_uni_TN and select Case_1_uni_TN as the model name. Click Continue. In the Edit Job dialog box, click OK to accept the default job settings. 10. Submit the job Case_1_uni_TN for analysis and view the results 11. Similarly, create models with names Case_1_uni_CMP for uniaxial compression, Case_1_hydro_CMP for hydrostatic compression, and Case_1_shear for pure shear state of deformation. Modify these models as follows: • Apply a displacement of -0.01 m in BC-4 and delete BC-5 and BC-6 to obtain uniaxial compression. • Apply displacements of -0.01 m for BC-4, BC-5, and BC-6 in order to obtain hydrostatic compression. • A pure shear state of deformation can be obtained by removing BC-2, BC5, and BC-6, modifying BC-4 so that now it specifies a displacement of 0.01 m along the 2-direction, and modifying BC-1 to prescribe zero displacement along directions 1 and 2. 12. Create and submit jobs for these models and visualize the results.



W1.4

Plot the locus of stress points for uniaxial tension, uniaxial compression, hydrostatic tension, hydrostatic compression, and pure shear on the meridional plane (pressure p vs. Mises stress q). Note: Detailed instructions are included in the workshop answers. Question W1–2: Derive the relationship between the pressure and Mises stress for uniaxial tension. Compare with the plot in the meridional plane. Question W1–1:

Case 2: Drucker-Prager Plasticity In Case 1 the material of the cube was linear elastic. In Case 2 you will add DruckerPrager plasticity and investigate the cube’s response to uniaxial tension and compression by performing the following steps: 1. Copy the model cube to a new model named Case_2_uni_CMP. 2. In the Model Tree branch for the model Case_2_uni_CMP, click mouse button 3 on the material Elastic and select Rename from the menu that appears. Rename this material to Drucker-Prager. 3. Double-click the renamed material Drucker-Prager to edit the material. 4. In the Edit Material dialog box select Mechanical→Plasticity→DruckerPrager. In the respective data columns specify the angle of friction as 30.0, flow stress ratio as 1.0 and dilation angle as 30.0. The units of the angles are in degrees. 5. Under the Drucker-Prager options of the material editor, click Suboptions and select Drucker-Prager Hardening from the menu that appears. Specify the yield stress as 30.0e6 Pa at an absolute plastic strain of 0.0. This data is for perfect plasticity as the yield stress is not dependent on the plastic strain. Select OK in the Drucker-Prager Hardening dialog box and the Edit Material dialog box. 6. In the Model Tree, double-click Section-1 to edit the section definition. Abaqus warns you that the section’s material no longer exists. Click Dismiss to close the warning. In the section editor, click OK to define the section with the material Drucker-Prager. 7. Delete boundary conditions BC-4, BC-5, and BC-6. 8. In the Model Tree branch for the model Case_2_uni_CMP, double-click Loads to create a new load. Select Pressure as the load type and click Continue. 9. Select the surface parallel to the face lying on the 2-3 plane as the surface for applying the pressure load (as shown in Figure W1-2) and click on Done in the



W1.5

prompt area. Specify the magnitude of the pressure load as 25.0e6 Pa and click on OK. This pressure load will subject the cube to uniaxial compression.

Figure W1-2: Selection of the surface lying parallel to the face on the 2-3 plane. 10. In the Model Tree, double-click Steps to create a new general static step named Step-2. In the step editor click the Incrementation tab and change the initial increment size to 0.025. 11. In the Model Tree, click mouse button 3 on the Loads container and select Manager from the menu that appears. In Step-2 increase the magnitude of the pressure load to 40.0e6 Pa. 12. In the Model Tree, double-click Jobs to create a new analysis job. Name the job Case_2_uni_CMP. Submit this job. 13. The job will not complete due to some errors. After the execution is over, click mouse button 3 on the job and select Results from the menu that appears. You should be able to visualize results from Step-1 and Step-2 in the Visualization module. You will find that Step-1 completed successfully, but Step-2 did not run to completion. For Step-1, sketch the yield surface on the meridional plane and plot the stress history along the loading direction. Confirm that the stress is in equilibrium with the external load. Question W1–4: Why does the job not complete Step-2? Question W1–3:

14. Create a new model named Case_2_uni_CMP_TN by copying the model Case_2_uni_CMP. 15. Delete the pressure load defined for the previous compression test. In order to investigate the performance of the Drucker-Prager material model in compression as well as tension, you will now perform a uniaxial compression test followed by a uniaxial tension test. Use two analysis steps for this purpose.



W1.6

16. In Step-1 apply uniaxial compression by prescribing a displacement of -0.01 m along the 1-direction on the face of the cube that is parallel to that lying on the 2-3 plane (as shown in Figure W1-2). In Step-2 change the boundary conditions on this face so that this face now has a prescribed displacement of 0.01 m, which will subject the cube to uniaxial tension. 17. Create a job named Case_2_uni_CMP_TN. Submit this job. Question W1–5:

Plot the history of direct stress (S11), the Mises stress, and the stress locus on the meridional plane. Does the material yield at the correct yield stress in tension and compression?

18. Plot the equivalent deviatoric plastic strain as a function of the volumetric plastic strain. a. Switch to the Visualization module and open the output database file Case_2_uni_CMP_TN.odb with the read-only option toggled off in the Open Database dialog box. b. Run the Python script ws_geotech_dev_tensor.py. This script computes the deviatoric and volumetric components of tensor output variables. c. In the Get Input dialog box, specify PE as the output variable. d. The script results are written to the output database as field data. To save the field data as X-Y data do the following: i. In the Results Tree, double-click XYData. ii. In the Create XY Data dialog box, select ODB field output as the data source. Click Continue. iii. In the Variables tabbed page of the XY Data from ODB Field Output dialog box, choose Integration Point as the variable position. Toggle on EQ_DEV_PE: Deviatoric, which is the equivalent deviatoric plastic strain. Toggle on VOL_PE: Volumetric, which is the volumetric plastic strain. iv. In the Elements/Nodes tabbed page, click Edit Selection, and select the element displayed in the viewport. Click Done in the prompt area. v. In the XY Data from ODB Field Output dialog box, click Save. e. In the Results Tree, double-click XYData. f. In the Create XY Data dialog box, choose Operate on XY data as the source and click Continue. g. In the Operate on XY Data dialog box, select combine(X,X) from the list of operators. Select EQ_DEV_PE and click Add to Expression. Similarly, add VOL_PE to the expression. h. Click Plot Expression.



W1.7

Question W1–6:

Compare these plastic strain results to a Drucker-Prager material with a zero dilation angle. You will need to rerun the analysis after modifying the value of the dilation angle. What do you observe?

Case 3: Drucker-Prager with Cap Plasticity In Case 1 the material of the cube was linear elastic and in Case 2 the material was Drucker-Prager. In Case 3 you will modify the material of the model to Drucker-Prager with Cap Plasticity and investigate the response of the model to uniaxial compression. 1. Create a new model by copying the model cube. Name the new model Case_3. 2. Rename material Elastic to Cap. 3. In the Model Tree, double-click Cap to edit the material. In the Edit Material dialog box, select Mechanical→Plasticity→Cap Plasticity. In the respective data columns specify the material cohesion as 2.4226e7 Pa, angle of friction as 30.0, cap eccentricity as 0.1, initial yield surface position as 0.002, transition surface radius as 0.0, and flow stress ratio as 1.0. 4. Select Suboptions→Cap Hardening and specify the yield stress vs. volumetric plastic strain according to the data in Table W1-1. Click OK in the Suboptions Editor and the Edit Material dialog boxes. pb (Pa) 5.0e6 10.0e6 30.0e6 60.0e6 100.0e6 200.0e6 500.0e6

εvol 0.000 0.001 0.002 0.003 0.004 0.005 0.006

Table W1–1: Cap Hardening Data. 5. Edit Section-1 so that it refers to the material Cap. Question W1–7:


Sketch the initial yield surface and the flow potential on the meridional plane.


W1.8

Next you will modify the model to perform a uniaxial compression test by applying a displacement of -0.05 m along direction 1 to the face of the cube that is parallel to the one lying on the 2-3 plane. 5. Delete boundary conditions BC-5, and BC-6. Edit BC-4 to specify the prescribed displacement of -0.05 m along the 1-direction. 6. Create a job named Case_3 and submit it for analysis. After the job has completed view the results. Question W1–8: Plot the equivalent deviatoric plastic strain as a function of the

volumetric plastic strain. (Use the Python script ws_geotech_dev_tensor.py as described previously.) Compare the result with the behavior obtained with the Drucker-Prager model in Question W1–6 and explain the behavior of the Cap model. Question W1–9: Sketch the updated flow potential on the meridional plane. Question W1–10: Plot the stress component S11 vs. the strain component E11. Explain the response.

Case 4: Clay Plasticity In Case 4 you will use the Clay material model. The Clay model in Abaqus does not have cohesion, and hence the model will first be subjected to an initial stress state of hydrostatic compression. This will occur in Step 1 of the analysis. Subsequently, in Step 2 of the analysis, the response of the model to uniaxial compression will be investigated. Further, in Step 3, the compression will be replaced by tension. The compression and tension in Steps 2 and 3 will be applied through the use of appropriate prescribed displacement boundary conditions. 1. Create a new model by copying the model cube. Name the new model Case_4. 2. Rename the material Elastic to Clay. 3. In the Model Tree, double-click Clay to edit the material. In the Edit Material dialog box, select Mechanical→Plasticity→Clay Plasticity. From the list of Hardening options, choose Tabular. In the respective data columns specify the stress ratio as 0.577 (which is tan 30 degrees), initial volumetric plastic strain as 0.002, wet yield surface size as 1.0, and flow stress ratio as 1.0. 4. Select Suboptions→Clay Hardening and specify the yield stress vs. volumetric plastic strain according to the data in Table W1-1, and click OK in the Suboptions Editor and the Edit Material dialog boxes. 5. Edit Section-1 so that it refers to the material Clay.



W1.9

6. Delete BC-4, BC-5, and BC-6. Apply a pressure load of magnitude 30.0e6 Pa in Step-1 to the three faces of the cube that do not have any displacement boundary conditions applied. 7. Apply an initial hydrostatic compressive stress of 30.0e6 Pa that balances the applied pressure load as follows: a. From the main menu bar, select Model→Edit Keywords→Case_4. b. In the Edit Keywords, Model: Case_4 dialog box each keyword is displayed in its own block. The ∗INITIAL CONDITIONS option must be placed in the model data portion of the input file (before the ∗STEP option). Select the text block just before the *STEP option, and click Add After to add an empty text block. c. In the new text block, type the following: *Initial Conditions, type=stress Set-1, -30.e6, -30.e6, -30.e6, 0.0, 0.0, 0.0

d. Click OK. 7. Create a general static step named Step-2. In the step editor, select the Incrementation tab. Change the initial increment size to 0.025. Similarly, create Step-3 with an initial increment size of 0.025. 8. In Step-2, apply uniaxial compression relative to the initial state by prescribing a uniform displacement of −0.05 for the displacement along 1-direction on the face of the cube that is parallel to that lying on the 2-3 plane. In Step-3, change the boundary conditions on this face so that this face now has a prescribed displacement of 0.05, which will subject the cube to uniaxial tension. 9. Create a job named Case_4 and submit it for analysis. After the job has completed view the results. Question W1–11: Sketch the initial yield surface on the meridional plane. Question W1–12: Plot the direct stress component S11 vs. the analysis time.

Compare the behavior in tension and compression with an equivalent Cap model.

Note: A script that creates complete models described in these instructions is available for your convenience. Run this script if you encounter difficulties following the instructions outlined here or if you wish to check your work. The script is named ws_geotech_geomaterials_answer.py and is available using the Abaqus fetch utility.



Answers 1 Material Models for Geotechnical Applications Interactive Version

Question W1–1:

Answer:

Plot the locus of stress points for uniaxial tension, uniaxial compression, hydrostatic tension, hydrostatic compression, and pure shear on the meridional plane (pressure p vs. Mises stress q). A plot of pressure p vs. Mises stress, q, can be obtained for a representative problem as follows: 1. In the Results Tree, expand the History Output container underneath the output database named Case_2_uni_CMP_TN.odb. 2. From the list of available history data, click mouse button 3 on the Mises stress (SINV: MISES) variable. From the menu that appears, select Save As and name the X–Y data Mises. 3. Similarly, select the pressure (SINV: PRESS) variable. Save the X–Y data as Press. 4. In the Results Tree, double-click XYData. 5. In the Create XY Data dialog box, choose Operate on XY data as the source and click Continue. 6. In the Operate on XY Data dialog box, select combine(X,X) from the list of operators. Select Press and click Add to Expression. Repeat this for Mises. The final expression is: combine ("Press", "Mises")

7. Click Save As to save the Mises vs. Pressure data for this model. Repeat this procedure for all five stress states. Plot all five Mises vs. Pressure curves in one X–Y plot (expand the XYData container, select the curves, and click mouse button 3; from the menu that appears, select Plot). Customize the plot (double-click plot components in the viewport or use Options→XY Options → …).



WA1.2

Shear

UC

UT

HC

HT

Question W1–2:

Derive the relationship between the pressure and Mises stress for uniaxial tension. Compare with the plot in the meridional plane.

Answer:

For uniaxial tension we have σ11 = σt and all other stress components as zero. Pressure p = −(σ11 + σ22 + σ33)/3 = −σt /3 Deviatoric stress component

S11

= σt + p = 2σt /3

S22

= −p

S33 = −p All other deviatoric stress components are zero. Mises stress

= {(3/2) S:S }0.5 = σt

Thus for uniaxial tension, Mises stress is related to the pressure by a factor of −3. The plot in the meridional plane agrees with this relationship.



WA1.3

Question W1–3:

For Step-1, sketch the yield surface on the meridional plane and plot the stress history along the loading direction. Confirm that the stress is in equilibrium with the external load.

Answer:

The Mises stress, q, has the same magnitude as the direct stress σ11, but it is always positive, i.e. q = |σ11|. If we apply an external pressure of 25 MPa on a face along the global 1-direction, we find that σ11 = −25 MPa. As this state of stress is uniaxial compression, we have Mises, q, as 25 MPa, and pressure, p, as 25/3, that is, 8.333 MPa. The yield point in compression, σc, is 30 MPa (same as the input value). The yield point in tension, σt, has a lower value than the yield point in compression. The yield stress in tension, σt , is related to the yield stress in compression, σc , as follows: We have the yield criterion as: F = q – p tan β − d = 0 For uniaxial tension, q = σt and pressure p = −σt /3. For uniaxial compression, q = σc and pressure p = σc /3. Substituting these values in the expression for F and eliminating d, we get,

σt = σc {(1−(tan β)/3)/(1+(tan β)/3)}. Therefore, for σc = 30 MPa and β = 30 degrees, we get,

σt = 20.32 MPa. The yield surface and the stress history for Step 1 can be depicted as follows:

Yield surface Mises q

(10, 30) (25/3 , 25)

(−20.32/3, 20.32)

Stress history

(0, 0)


Pressure p


WA1.4

Question W1–4: Why does the job not complete Step-2? Answer: The analysis fails with a numerical singularity, because the

applied load is larger than the limit (or failure) yield stress in compression. Question W1–5:

Answer:

Plot the history of direct stress, (S11), the Mises stress, (q), and the stress locus on the meridional plane. Does the material yield at the correct yield stress in tension and compression? Yes, the material does yield at the expected stress values in compression and tension. The yield stress in compression is 30 MPa and the yield stress in tension is 20.32 MPa, as discussed in the answer to Question W1–3. History plots of direct stress, (S11), the Mises stress, (q), and the stress locus on the meridional plane follow:

Mises

σt = 20.32 MPa

|σc| = 30.0 MPa

S11



WA1.5

30.0 MPa

Compression 20.32 MPa

Tension

In this plot the unloading curve appears to deviate from the loading curve only because of the coarseness of the incrementation.

Question W1–6:

Compare these plastic strain results to a Drucker Prager material with a zero dilation angle. You will need to rerun the analysis after modifying the value of the dilation angle. What do you observe?

Answer:

The volume strain is zero when the dilation angle is zero; i.e. the material is incompressible. Hybrid elements may be required to prevent locking. The slope of the line in a plot of equivalent deviatoric plastic strain vs volumetric plastic strain is the tangent of the dilation angle. The Drucker Prager model predicts continuing dilation, regardless of the mode of deformation. In practice, materials show a limited amount of dilation.



WA1.6

Dilation angle = 30

Dilation angle = 0

Question W1–7:


Answer:

Yield surface

Flow potential

Mises q (0, 24.226)

β = 30o Pressure p


(30, 0)


WA1.7

Question W1–8: Plot the equivalent deviatoric plastic strain as a function of the

Answer:

volumetric plastic strain. (Use the Python script ws_geotech_dev_tensor.py as described previously.) Compare the result with the behavior obtained with the Drucker Prager model in Question W1–6 and explain the behavior of the Cap model. The Cap model exhibits a limit on the amount of dilation. The Drucker-Prager model exhibits no dilation limit.

Question W1–9: Sketch the updated flow potential on the meridional plane. Answer:

Flow potential

Mises q (0, 24.226)

Yield surface (10, 30)

β = 30o (13, 0)


Pressure p


WA1.8

Question W1–10: Plot the stress component S11 vs. the strain component E11.

Explain the response. Yielding on the shear surface results in perfect plastic response; i.e. the shear surface is a failure surface.

Answer:

Question W1–11: Sketch the initial yield surface on the meridional plane. Answer:

Critical state line

Mises q

Yield surface Slope 3 Tension

Compression

Μ = 30o Pressure p


(30, 0)


WA1.9

Question W1–12: Plot the direct stress component S11 vs. the analysis time.

Answer:

Compare the behavior in tension and compression with an equivalent Cap model. In the Clay model, the yielding results in softening on the dry side, while hardening takes place on the wet side. The failure state occurs when the stress reaches the critical state line. The Clay model exhibits softening during tensile failure and hardening in compressive failure. The Cap model, however, would predict perfect plastic failure on the failure surface. If the shear failure surface of the Cap model matches with the critical state line of this Clay model, then the Cap model would predict failure without any softening at a relative uniaxial tensile stress of 14.53 MPa (the tension limit in the figure above, that is, at S11= −15.47 MPa). In the compressive case, it will fail without any hardening at a relative uniaxial compressive stress of 21.45 MPa (the compression limit in the figure above, that is, at S11= −51.45 MPa).

Step-1 Hydrostatic compression © Dassault Systèmes, 2008

Step-2 Uniaxial compression

Step-3 Uniaxial tension


Workshop 2 Pore Fluid Flow Analysis: Consolidation Interactive Version Note: This workshop provides instructions on using the pore fluid flow modeling feature in Abaqus in terms of the Abaqus GUI interface. If you wish to use the Abaqus Keywords interface instead, please see the “Keywords” version of these instructions. Please complete either the Keywords or Interactive version of this workshop.


Model pore fluid flow in fully saturated problems. Specify initial conditions that equilibriate loading. Use Abaqus/CAE to build simple models, and view and interpret the pore fluid flow phenomenon.

In this workshop you will perform a consolidation analysis on a fully saturated model. This problem illustrates how consolidation analyses can be performed in Abaqus by using the available coupled stress diffusion analysis cabability. The one-dimensional Terzaghi consolidation problem will be modeled using two dimensional plane strain elements. In this problem a body of soil is confined by impermeable, smooth, rigid walls on all but the top surface where perfect drainage is possible. At the start of the analysis, a load is applied suddenly on the top surface of the body. We wish to predict the response of the soil as a function of time after the load application. Gravity is neglected. This problem is similar to “The Terzaghi Consolidation Problem,” Section 1.14.1 of the Abaqus Benchmarks Manual. The model is 50 m wide and 100 m high. The soil is considered to be elastic with a modulus of 1.0e8 Pa, Poisson’s ratio of 0.3, and permeability of 0.002 m/s. The initial spatially constant void ratio value is 1.1. The soil is considered to be fully saturated. The



W2.2

model is meshed with 10 elements through its height and 1 element across its width. The geometrical configuration of the model is shown in Figure W2-1.

Figure W2-1: Geometrical configuration of the consolidation specimen.

Preliminaries 1. Enter the working directory for this workshop: ../geotechnical/interactive/consolidate

2. Run the script ws_geotech_consolidate.py using the following command: abaqus cae startup=ws_geotech_consolidate.py

The above command creates an Abaqus/CAE database named consolidate.cae in the current directory. The model geometry, material properties, and displacement boundary conditions have already been created in the database in the model named consolidate. The model consolidate does not have any pore pressure boundary conditions specified on the exterior edges of the soil. This is the natural zero pore fluid flux boundary condition, which represents an undrained condition. The application of the external load will result in a build up of pore pressure in the model. You will simulate this phenomenon in the first step of the consolidation analysis, wherein you will apply a pressure load on the top surface of the model and not allow for any pore fluid drainage across the whole model boundary. In the second step of the analysis you will maintain the pressure load from the first step and allow for pore fluid drainage across the top surface of the model.



W2.3

Initial Void Ratio In this section you will specify an initial void ratio of 1.1 in the model as a pre-defined initial field. You will use the keywords editor to define the initial void ratio for the entire model (Set-1) using the *INITIAL CONDITIONS option This specification can be achieved in CAE as follows: 1. In an Abaqus/CAE session, open the database named consolidate.cae, which was created by running the setup script ws_geotech_consolidate.py. 2. From the main menu bar, select Model→Edit Keywords→consolidate. 3. In the keywords editor that appears, add a new option block just before the comment **STEP: Step-1. 4. In the new text block, type the following: *Initial conditions, type=ratio Set-1, 1.1, 0.0, 1.1, 100.0

5. Click OK to close the keywords editor.

Pressure loading The first analysis step was created by the setup script. In this section you will complete the step definition by defining the pressure load. 1. In the Model Tree, double-click Loads to create a new load. 2. Select Step-1 as the load definition step and Pressure as the load type. 3. Select the top surface in the model as the surface to which the load will be applied, and click Done in the prompt area. 4. In the Edit Load dialog box, specify a load magnitude of 100.0 Pa and click OK. Question W2–1:


What will be the state of effective stress in the model when a pressure load of 100.0 is applied on the top surface of the model and pore fluid drainage is prevented from taking place?


W2.4

Consolidation In this section you will create the second step with a total step time of 10.0. In this step you will allow for pore fluid drainage across the top boundary by specifying appropriate pore pressure boundary conditions on the top surface. You also need to specify an appropriate value for UTOL, which is a parameter that governs the accuracy of the time integration for the transient consolidation analysis. Question W2–2:

Determine an appropriate value for UTOL.

The element size and the time increment size are related to the extent that time increments smaller than a certain size give no useful information. This coupling of the spatial and temporal approximations is most obvious at the start of diffusion problems, immediately after prescribed changes in the boundary values. Question W2–3:

Using the time-space criterion described in Lecture 5, determine the initial time increment size to avoid overshoot in the initial solution.

The values of UTOL and the initial time increment size will be specified as data for the second step of the analysis and the analysis. 1. In the Model Tree, double-click Steps to create a new step. In the Create Step dialog box that appears, accept Soils as the procedure and click Continue. 2. In the Step editor, specify a time period of 10.0 and toggle off the Include creep/swelling/viscoelastic behavior option. 3. Select the Incrementation tab and specify the initial time increment size at the appropriate location. Specify the value of the Max. pore pressure change per increment (UTOL). Click OK. 4. In the Model Tree, double-click BCs to create a new boundary condition. 5. In the Create Boundary dialog box, select Step-2 as the boundary condition definition step, Other as the Category, and Pore pressure as the boundary condition type. 6. Select the top surface of the model as the region for the new boundary condition, and click Done in the prompt area. 7. In the Edit Boundary Condition dialog box specify a pore pressure of 0.0. This specification will lead to a constant zero pore pressure on the top surface and allow for pore fluid drainage across this boundary. 8. Click OK in the Edit Boundary Condition dialog box.



W2.5

Visualizing the results The model definition is now complete. In this section you will run the analysis and study the results. 1. In the Model Tree, click mouse button 3 on the job consolidate and select Submit from the menu that appears. 2. After the job completes, click mouse button 3 on the job consolidate and select Results from the menu that appears. Abaqus switches to the Visualization module and opens the output database file consolidate.odb. Plot the pore pressure and vertical stress along the vertical edge at different times during the analysis. Question W2–5: Plot the pore pressure and vertical stress as a function of time at the bottom of the specimen. Question W2–4:

Note: A script that creates the complete model described in these instructions is available for your convenience. Run this script if you encounter difficulties following the instructions outlined here or if you wish to check your work. The script is named ws_geotech_consolidate_answer.py and is available using the Abaqus fetch utility.



Answers 2 Pore Fluid Flow Analysis: Consolidation Interactive Version

Question W2–1:

Answer:

Question W2–2: Answer:

What will be the state of effective stress in the model when a pressure load of 100.0 Pa is applied on the top surface of the model and pore fluid drainage is prevented from taking place? A stress state with a uniform pore pressure equal to the external load throughout the body with no effective stress carried by the soil skeleton is obtained. The fluid carries the instantaneous load.

Determine an appropriate value for UTOL. UTOL which is one order of magnitude less than the applied pressure should give moderate accuracy.

Question W2–3: Using the time-space criterion described in Lecture 5,

Answer:

determine the initial time increment size to avoid overshoot in the initial solution. The initial time increment can be computed as follows: ∆t ≥ {γw/(6Ek)}(∆h)2 We have ∆h = 10 m, E = 108 Pa, k = 2×10-4 m/s, N/m3. Hence,

γw =1.0

∆t ≥ 0.00083 s.



Plot the pore pressure and vertical stress along the vertical edge at different times during the analysis. A path needs to be created in the Visualization module so that results along this path can be plotted. The path creation and


WA2.2

plotting of X-Y results along this path can be accomplished as follows: 1. In the Results Tree, double-click Paths. 2. In the Create Path dialog box, select Node list and click Continue. 3. In the Edit Node List Path dialog box, click Add After. 4. In the viewport, select all consecutive nodes on one vertical edge from the bottom to the top. Click Done in the prompt area. 5. The selected nodes are listed in the Edit Node List Path dialog box. Click OK to complete the path definition. 6. In the Results Tree, double-click XYData. 7. In the Create XY Data dialog box, select Path and click Continue. 8. In the XY Data from Path dialog box, select the frame and the field output variable by clicking Step/Frame and Field Output respectively. Plot and save stress (S22) and pore pressure (POR) X-Y data for a few frames. 9. Expand the XYData container. 10. Select the saved data names and click mouse button 3; from the menu that appears, select Plot. A combined X-Y plot for the selected curves will be displayed. The plots indicate that the pore pressure decreases and the values of S22 increase over time. The pore pressure values are almost zero after about t=1.25.



WA2.3

t = 0.001 t = 0.022

t = 0.251

t = 1.248

t = 0.001 t = 0.022

t = 0.251

t = 1.248



WA2.4


Plot the pore pressure and vertical stress as a function of time at the bottom of the specimen. The pore pressure and stress results were written to the output database as history data during the analysis. This data is accessible from the History Output dialog box (Result→ History Output). To plot the stress magnitude you must save the history data and operate on it (double-click XYData, and select Operate on XY data as the data source). The values of POR decrease over time as the pore water permeates away from the top surface. The values of S22 simultaneously increase in magnitude. Hence, as the consolidation progresses, the load initially carried by the pore fluid is transferred to the solid grains.

S22 magnitude

POR



Workshop 3 Pore Fluid Flow Analysis: Wicking Interactive Version Note: This workshop provides instructions on using the pore fluid flow modeling feature in Abaqus in terms of the Abaqus GUI interface. If you wish to use the Abaqus Keywords interface instead, please see the “Keywords” version of these instructions. Please complete either the Keywords or Interactive version of this workshop.


Model pore fluid flow in partially saturated problems. Specify initial conditions that equilibriate loading. Use Abaqus/CAE to build simple models, and view and interpret the pore fluid flow phenomenon.

In this workshop you will perform a pore fluid absorption analysis on a partially saturated specimen. This problem illustrates how pore fluid absorption analyses can be performed in Abaqus using the available coupled stress diffusion analysis cabability. We consider a one-dimensional wicking test where the absorption of fluid takes place against the gravity load caused by the weight of the fluid. In this test the fluid is made available at the bottom of a partially saturated material column, and the material absorbs as much fluid as the weight of the rising fluid permits. The column of material is kinematically constrained in the horizontal direction so that all material deformations will be in the vertical direction. In this sense the problem is one-dimensional. The geometrical configuration of the model is shown in Figure W3–1.



W3.2

Figure W3–1 Geometrical configuration of the wicking model.

Preliminaries 1. Enter the working directory for this workshop: ../geotechnical/interactive/wicking

2. Run the script ws_geotech_wicking.py using the following command: abaqus cae startup=ws_geotech_wicking.py

The above command creates an Abaqus/CAE database named wicking.cae in the current directory. The model geometry, material properties, and loading history have already been created in the database in the model named wicking. The units used for the material data in the database are in terms of metric tonnes, meters (m), and seconds (s), and therefore the force will be in kiloNewtons (kN) and stress in kiloPascals (kPa). The column of material is 1.0 m high and 0.1 m wide. It is meshed with 10 plane strain elements with displacement and pore pressure degrees of freedom at the nodes (CPE4P). For the mechanical properties we assume that the material is elastic



W3.3

with a Young's modulus of 50000 Pa and a Poisson's ratio of 0.0. The dry mass density of the material is 100 kg/m3. The permeability of the fully saturated material is 3.7 × 10−4 m/sec. We assume that the permeability varies as a cubic function of saturation, which is the default setting for partially saturated permeability. The specific weight of the fluid (water) is 104 N/m3 and the bulk modulus is 2 GPa. The capillary action in the porous medium is defined by absorption/exsorption curves shown in Figure W3–2. These curves give the (negative) pore pressure versus saturation relationships for the absorption and exsorption behaviors. The transition between absorption and exsorption and vice-versa takes place along a scanning slope which in this example is set by default to be 1.05 times the largest slope of any branch in the absorption/exsorption curves.

Figure W3–2 Absorption and exsorption curves for the model material.



W3.4

Initial conditions The initial condition for saturation is 5% and the initial void ratio is 5.0 throughout the column. The initial conditions for the pore pressure must have a gradient that is equal to the specific weight of the fluid so that, according to Darcy's law, there is no initial flow. For this purpose we assume the initial pore pressures vary linearly from −12000 Pa at the bottom of the column to −22000 Pa at the top of the column. These initial conditions satisfy the pore pressure/saturation relationship in the sense that they are between the absorption and exsorption curves. Question W3–1:

Determine the vertical coordinate of the phreatic surface.

The initial conditions for the effective stress are calculated from the densities of the dry material and the fluid, the initial saturation, the void ratio, and the initial pore pressures using equilibrium considerations and the effective stress principle. It is important to specify the correct initial conditions for this type of problem; otherwise, the system may be so far out of equilibrium initially that it may fail to start because converged solutions cannot be found. Question W3–2:

Determine the initial conditions for the effective stress.

The initial conditions for the pore pressure, effective stress, saturation, and void ratio must be specified as input data for the analysis. Perform the following procedures to incorporate these data into the model database: 1. In an Abaqus/CAE session, open the database named wicking.cae, which was created by running the setup script ws_geotech_wicking.py. 2. From the main menu bar, select Model→Edit Keywords→wicking. 3. In the keywords editor that appears, add a new option block just before the comment **STEP: Step-1. 4. In the new option block, type the following with the appropriate values for the quantities in italics (Recall that the model is defined in terms of metric tonnes, meters, and seconds): *INITIAL CONDITIONS,TYPE=SATURATION Set-1,.05 *INITIAL CONDITIONS,TYPE=PORE PRESSURE Set-1,

pore press. at elevation 1.0,


1.0,


0.0


W3.5 *INITIAL CONDITIONS,TYPE=RATIO Set-1,5. *INITIAL CONDITIONS,TYPE=STRESS,GEOSTATIC Set-1,

stress at elevation 1.0,

1.0,


0.0, 0., 0.

5. Click OK to close the keywords editor.

Analysis and results The analysis will be performed in two steps. The first step is a geostatic analysis step with a short time duration. This geostatic step ensures that the model is in equilibrium at the initiation of step 2, which models the transient wicking phenomenon. For step 2, the convergence criterion for the ratio of the largest solution correction to the largest corresponding incremental solution value has been increased from the default value of 0.01 to a value of 1.0. This increase helps in obtaining convergence in the early stages of the analysis wherein the rate of change of saturation and pore pressure is large at the bottom of the wicking column. Steps 1 and 2 are already included in the model wicking. In this section you will create and submit an analysis job. Once the job is complete, you will study the results 1. In the Model Tree, click mouse button 3 on the job wicking and select Submit from the menu that appears. 2. After the job completes, study the results. Question W3–3:

Plot the vertical stress distribution, the pore pressure, and saturation along the height of the sample at different times.

Question W3–4:

Plot the time history of the volume of the fluid absorbed by the column.

Note: A script that creates the complete model described in these instructions is available for your convenience. Run this script if you encounter difficulties following the instructions outlined here or if you wish to check your work. The script is named ws_geotech_wicking_answer.py and is available using the Abaqus fetch utility.



Answers 3 Pore Fluid Flow Analysis: Wicking Interactive Version

Question W3–1: Determine the vertical coordinate of the phreatic surface. Answer:

The phreatic surface lies at a vertical coordinate of −1.2 m. This value can be found by linearly extrapolating the pore pressure distribution to obtain the vertical coordinate at which the pore pressure is zero. The known pore pressure values are −22000 Pa, and −12000 Pa, at vertical coordinates of 1.0 m, and 0.0 m, respectively.

Question W3–2: Determine the initial conditions for the effective stress. Answer:

The initial void ratio is spatially constant and is 5.0. The initial saturation is spatially constant and is 0.05. Denoting γd as the dry weight density of the soil, γw as the weight density of the pore fluid, e as the void ratio, n as the porosity, s as the saturation, z as the vertical coordinate, σz as the total stress, and σ z as the effective stress, we have: (dσz/dz) = γd + snγw, which implies,

σz = (γd + snγw)(z−z0). For γd = 103, s = 0.05, n = e/(1+e) = 5/6, and z0 = 1, σz = 1416.7 (z−1). Effective stress σ z = σz + suw At z = 0.0, pore pressure uw = −12.0 kPa At z = 1.0, pore pressure uw = −22.0 kPa Therefore, At z = 0.0, effective stress σ z= −2.0167 kPa At z = 1.0, effective stress σ z= −1.1000 kPa



WA3.2


Plot the vertical stress distribution, the pore pressure, and saturation along the height of the sample at different times. At steady state the pore pressure gradient must equal the weight of the fluid and hence the pore pressure varies linearly with height and saturation. Thus, points close to the bottom of the column are fully saturated, while those at the top are still at 5% saturation. In the following plots the stress and pore pressure values are in kiloPascals:

t=0

t = 15 t = 171 t = 2000 t = 69526 t = 200000

Distance from the bottom surface (m)



WA3.3

t = 200000 t = 15 t = 69526 t = 171

t=0

t = 2000


t = 15 t = 171 t = 2000 t = 69526 t = 200000

t=0




WA3.4

Question W3–4: Plot the time history of the volume of the fluid absorbed by the

column. Answer:

The volume of the fluid absorbed by the column will be given by the sum of the values of the output variable RVT (reactive fluid total volume) at all boundary nodes where the pore pressure is prescribed. The pore pressure is prescribed at nodes on the bottom of the column. The following plot, showing the total RVT value in cubic meters, indicates that the rate of increase in RVT decreases with time and the total volume of the absorbed fluid approaches a limiting value.



Workshop 1 Material Models for Geotechnical Applications Keywords Version Note: This workshop provides instructions on using geotechnical material models in terms of the Abaqus Keywords interface. If you wish to use the Abaqus GUI interface instead, please see the “Interactive” version of these instructions. Please complete either the Keywords or Interactive version of this workshop.


Create stress state diagrams. Define a number of geotechnical material models. Use Abaqus/Viewer to view and interpret material behavior.

Figure W1–1 Single element cube.



W1.2

In this workshop you will become familiar with some of the geotechnical material models in Abaqus. You will study the behavior of the 2x2x2 m cube shown in Figure W1–1. This cube will be subjected to uniaxial tension, uniaxial compression, pure shear, hydrostatic tension, and hydrostatic compression. Different material models will be assigned to the cube and the response of the cube for these material models will be simulated.

Preliminaries Enter the working directory for this workshop: ../geotechnical/keywords/geomaterials

This directory contains the input file geo_elastic_hydro_t.inp. The file geo_elastic_hydro_t.inp contains the input data to analyze the response of a cube subjected to hydrostatic tension. The cube is meshed with one reduced-integration solid element (C3D8R) as shown in Figure W1–1. The cube material is linear elastic with a Young’s modulus of 70 GPa and a Poisson’s ratio of 0.3. Node sets named x0, y0, z0, x1, y1, and z1 for nodes lying on the six faces of the cube have been defined in this input deck to help in the specification of the boundary conditions. These sets are shown in Figure W1–2. Two surfaces have also been defined to help in the specification of the loading conditions. The faces of the cube lying on the 2-3, 3-1, and 1-2 planes are constrained against displacements normal to these respective planes using boundary conditions applied to the node sets x0, y0, and z0, respectively. Prescribed displacements of 0.01 m along the outward normals of the other three faces are specified so that the cube is subjected to a state of hydrostatic tension. These prescribed displacements are applied to the nodes in sets x1, y1, and z1. For each of the problem cases in this workshop you will copy the original input file and modify the material data and loading conditions.

x1 x0

z1 z0

y1 y0 Figure W1–2 Node sets.



W1.3

Case 1: Linear Elasticity You will begin this workshop by executing the geo_elastic_hydro_t.inp input file and studying the behavior of the cube under hydrostatic tension. You will then create new input files by copying the original input file and modifying the boundary conditions to subject the cube to (a) uniaxial tension, (b) uniaxial compression, (c) hydrostatic compression, and (d) pure shear. After the boundary conditions are modified in an input file, you will run the analysis, check for any modeling errors, and then view the results. 1. At the command prompt, enter the following command to submit the hydrostatic tension analysis: abaqus job=geo_elastic_hydro_t

2. While the job is running, you can monitor the execution of the job by reading the message (.msg) file, the status (.sta) file, or the log (.log) file. 3. After the job has completed, enter the following command to visualize the results in Abaqus/Viewer: abaqus viewer odb=geo_elastic_hydro_t

4. Look at the deformed shape and plot contours of the nonzero components of stress (S11, S22, and S33). 5. Copy the input file geo_elastic_hydro_t.inp to geo_elastic_uni_t.inp. 6. In a text editor, modify the boundary conditions in the input file geo_elastic_uni_t.inp to subject the cube to uniaxial tension in the 1direction. a. Open geo_elastic_uni_t.inp in a text file editor and review its contents. b. In Step-1, remove the displacement boundary conditions applied on node sets y1 and z1. Removing these boundary conditions will allow free expansion or contraction of the cube along the 2- and 3-directions. Do not delete the boundary condition on node set x1, which will apply the tension in the 1-direction. This input file is now ready for execution. 7. Save the model and use the following command to submit it for analysis: abaqus job=geo_elastic_uni_t

8. View the results after the job has completed. 9. Similarly, create input files named geo_elastic_uni_c.inp for uniaxial compression, geo_elastic_hydro_c.inp for hydrostatic compression, and geo_elastic_shear.inp for pure shear state of deformation. Modify these input files as follows: • Modify the input file geo_elastic_uni_c.inp so that a prescribed displacement of −0.01 m is now specified on node set x1, and delete the boundary conditions on node sets y1 and z1 in order to obtain uniaxial compression. © Dassault Systèmes, 2008


W1.4

• Modify the input file geo_elastic_hydro_c.inp, by changing the prescribed displacements at node sets x1, y1, and z1 to −0.01 m in order to obtain hydrostatic compression. • In the input file geo_elastic_shear.inp, a pure shear state of deformation can be obtained by removing the boundary conditions on node sets y0, y1, and z1, modifying the boundary condition on node set x1 so that now it specifies a displacement of 0.01 m along the 2-direction, and modifying the boundary condition on node set x0 to prescribe zero displacement along directions 1 and 2. 10. Execute the jobs geo_elastic_uni_c, geo_elastic_hydro_c, and geo_elastic_shear and visualize the results. Plot the locus of stress points for uniaxial tension, uniaxial compression, hydrostatic tension, hydrostatic compression, and pure shear on the meridional plane (pressure p vs. Mises stress q). Note: Detailed instructions are included in the workshop answers. Question W1–2: Derive the relationship between the pressure and Mises stress for uniaxial tension. Compare with the plot in the meridional plane. Question W1–1:

Case 2: Drucker-Prager Plasticity In Case 1 the material of the cube was linear elastic. In Case 2 you will add DruckerPrager plasticity and investigate the cube’s response to uniaxial tension and compression by performing the following steps: 1. Copy the input file geo_elastic_hydro_t.inp to a new file named geo_drucker_uni_c.inp. 2. In the new input file, rename the material ELASTIC to DRUCKER-PRAGER. Tip: You need to make this change on the *MATERIAL option and the *SOLID SECTION option which refers to the renamed material. 3. Add the following lines after the *MATERIAL option to include Drucker-Prager plasticity in the material model. In this plasticity definition the friction angle is 30 degrees, the ratio of flow stress in triaxial tension to flow stress in triaxial compression is 1.0, the dilation angle is 30 degrees, and the yield stress is 30 MPa at a plastic strain of 0.0. *DRUCKER PRAGER 30., 1., 30. *DRUCKER PRAGER HARDENING 3E+07, 0.



W1.5

4. Delete the prescribed boundary conditions on node sets x1, y1, and z1 applied in Step-1. 5. Replace the step boundary conditions with a 25 MPa uniform pressure load on the face of the cube that is parallel to that lying on the 2-3 plane (i.e., on surface surf_x1 shown in Figure W1–3): *DSLOAD surf_x1, P, 2.5E+07

surf_x1

Figure W1–3 Surface lying parallel to the face on the 2-3 plane. 6. After the *END STEP option, add the following lines to create a new static step named Step-2, wherein the applied pressure load is increased to 40 MPa. The initial time increment size for this step is 0.025. *STEP, name=Step-2 *STATIC 0.025, 1., 1E-05, 1. *DSLOAD surf_x1, P, 4E+07 *ENDSTEP

7. Save the file and execute job geo_drucker_uni_c. 8. The job will not complete due to some errors. The message (.msg) file, the status (.sta) file, and the log (.log) file will indicate that the job did not complete. You will find that Step-1 completed successfully, but Step-2 did not run to completion. Open the output database (.odb) file in Abaqus/Viewer in order to visualize the results from Step-1 and Step-2. For Step-1, sketch the yield surface on the meridional plane and plot the stress history along the loading direction. Confirm that the stress is in equilibrium with the external load. Question W1–4: Why does the job not complete Step-2? Question W1–3:



W1.6

In order to investigate the performance of the Drucker-Prager material model in compression as well as tension, you will now perform a uniaxial compression test followed by a uniaxial tension test. Use two analysis steps for this purpose. 9. Create a new input file named geo_drucker_uni_ct.inp by copying the input file geo_drucker_uni_c.inp. 10. Open the input file geo_drucker_uni_ct.inp in a text editor and delete the pressure load defined in both analysis steps. 11. In Step-1 apply uniaxial compression by prescribing a displacement of −0.01 m along the 1-direction on the face of the cube that is parallel to that lying on the 23 plane. Node set x1 contains these nodes. In Step-2 change the boundary conditions on this node set so that this face has a prescribed displacement of 0.01 m, which will subject the cube to uniaxial tension. 12. Execute job geo_drucker_uni_ct and visualize the results. Question W1–5:

Plot the history of direct stress (S11), the Mises stress, and the stress locus on the meridional plane. Does the material yield at the correct yield stress in tension and compression?

13. Plot the equivalent deviatoric plastic strain as a function of the volumetric plastic strain. a. In Abaqus/Viewer, open the output database file geo_drucker_uni_ct.odb with the read-only option toggled off in the Open Database dialog box. b. Run the Python script ws_geotech_dev_tensor.py. This script computes the deviatoric and volumetric components of tensor output variables. c. In the Get Input dialog box, specify PE as the output variable. d. The script results are written to the output database as field data. To save the field data as X-Y data do the following: i. In the Results Tree, double-click XYData. ii. In the Create XY Data dialog box, select ODB field output as the data source. Click Continue. iii. In the Variables tabbed page of the XY Data from ODB Field Output dialog box, choose Integration Point as the variable position. Toggle on EQ_DEV_PE: Deviatoric, which is the equivalent deviatoric plastic strain. Toggle on VOL_PE: Volumetric, which is the volumetric plastic strain. iv. In the Elements/Nodes tabbed page, click Edit Selection, and select the element displayed in the viewport. Click Done in the prompt area. v. In the XY Data from ODB Field Output dialog box, click Save. e. In the Results Tree, double-click XYData.



W1.7

f. In the Create XY Data dialog box, choose Operate on XY data as the source and click Continue. g. In the Operate on XY Data dialog box, select combine(X,X) from the list of operators. Select EQ_DEV_PE and click Add to Expression. Similarly, add VOL_PE to the expression. h. Click Plot Expression. Question W1–6:


Case 3: Drucker-Prager with Cap Plasticity In Case 1 the material of the cube was linear elastic and in Case 2 the material also included Drucker-Prager plasticity. In Case 3 you will modify the cube material to Drucker-Prager with Cap Plasticity and investigate the response of the model to uniaxial compression. 1. Create a new input file named geo_cap.inp by copying input file geo_elastic_hydro_t.inp. 2. In the new input file, rename the material ELASTIC to CAP. 3. Add the following lines after the *MATERIAL option in order to define a Cap plasticity model in which the material cohesion is 2.4226e7 Pa, material angle of friction is 30 degrees, cap eccentricity parameter is 0.1, initial yield surface position is 0.002, transition surface radius is 0.0, and flow stress ratio is 1.0: *CAP PLASTICITY 2.4226E+07, 30., 0.1, 0.002, 0., 1.

4. Yield stress vs. volumetric plastic strain data for this material model is included in Table W1–1. This cap hardening data needs to be specified in the input file. pb (Pa) εvol 5.0e6 0.000 10.0e6 0.001 30.0e6 0.002 60.0e6 0.003 100.0e6 0.004 200.0e6 0.005 500.0e6 0.006 Table W1–1 Hardening Data. Specify this data by adding the following lines after the *CAP PLASTICITY option:



W1.8

*CAP HARDENING 5E+06, 0. 1E+07, 0.001 3E+07, 0.002 6E+07, 0.003 1E+08, 0.004 2E+08, 0.005 5E+08, 0.006 Question W1–7:


Next you will modify the input deck in order to perform a uniaxial compression test by applying a displacement of −0.05 m along direction 1 to the face of the cube that is parallel to the one lying on the 2-3 plane. 5. Delete the boundary conditions on node sets y1 and z1. Edit the boundary condition on node set x1 in order to specify a prescribed displacement of −0.05 m along the 1-direction. 6. Execute job geo_cap. After the job has completed, view the results. Question W1–8: Plot the equivalent deviatoric plastic strain as a function of the

volumetric plastic strain. (Use the Python script ws_geotech_dev_tensor.py as described previously.) Compare the result with the behavior obtained with the Drucker-Prager model in Question W1–6 and explain the behavior of the Cap model. Question W1–9: Sketch the updated flow potential on the meridional plane. Question W1–10: Plot the stress component S11 vs. the strain component E11. Explain the response.

Case 4: Clay Plasticity In Case 4 you will use the Clay material model. The Clay model in Abaqus does not have cohesion, and hence the model will first be subjected to an initial stress state of hydrostatic compression. This will occur in Step 1 of the analysis. Subsequently, in Step 2 of the analysis, the response of the model to superimposed uniaxial compression will be investigated. Further, in Step 3, the compression will be replaced by tension. The compression and tension in Steps 2 and 3 will be applied through the use of appropriate prescribed displacement boundary conditions. 1. Create a new input file named geo_clay by copying the file geo_elastic_hydro_c.inp. 2. In the new input file, rename the material ELASTIC to CLAY.



W1.9

3. Add the following lines after the *MATERIAL option in order to include clay plasticity with a stress ratio of 0.577 (which is tan 30 degrees), an initial volumetric plastic strain of 0.002, a wet yield surface size of 1.0, and a flow stress ratio of 1.0: *CLAY PLASTICITY, HARDENING=TABULAR , 0.577, 0.002, 1., 1.

4. The hardening data for the clay plasticity model must be specified in terms of the hydrostatic pressure stress at yield as a function of the volumetric plastic strain. For this model use the hardening data provided in Table W1-1. Add the following lines after the *CLAY PLASTICITY option in order to specify this hardening data: *CLAY HARDENING 5E+06, 0. 1E+07, 0.001 3E+07, 0.002 6E+07, 0.003 1E+08, 0.004 2E+08, 0.005 5E+08, 0.006

The first step of this analysis involves balancing an applied hydrostatic pressure load with the internal forces arising from the initial hydrostatic stresses, thus ensuring equilibrium. The computed displacements in this first step will be small if the initial conditions are appropriate for the applied loading conditions. 5. In Step-1 delete the boundary conditions on node sets x1, y1, and z1 and replace them with a pressure load of 30.0e6 Pa. Apply the pressure load to the faces included in the surface surf_x1y1z1 shown in Figure W1–4. This pressure load will subject the cube to a state of hydrostatic compression. Replace the Step-1 boundary conditions with the following option: *DSLOAD surf_x1y1z1, P, 3E+07 surf_x1y1z1

Figure W1–4 Surface surf_x1y1z1.



W1.10

6. Define an initial hydrostatic compressive stress of 30.0e6 Pa that balances the applied pressure load so that the cube is close to equilibrium at the start of the analysis. The initial stress can be applied to the element set all. Add the following lines before the *STEP option to specify the initial hydrostatic stress: *INITIAL CONDITIONS, TYPE=STRESS all, -30.E6, -30.E6, -30.E6, 0., 0., 0.

In the second step you will subject this hydrostatically compressed cube to a superimposed uniaxial compression and in the third step you will replace the superimposed uniaxial compression with a superimposed uniaxial tension. The uniaxial compression and uniaxial tension are applied using appropriate displacement boundary conditions. 7. After the *END STEP option, add the following lines to create a new static step named Step-2, wherein the superimposed compression is achieved by applying prescribed displacements to the nodes in node set x1: *STEP, NAME=Step-2 *STATIC 0.025, 1., 1E-05, 1. *BOUNDARY x1, 1, 1, -0.05 *ENDSTEP

8. After the *END STEP option for Step-2, add the following lines to create a new static step named Step-3, wherein the superimposed tension is achieved by applying prescribed displacements to the nodes in node set x1: *STEP, NAME=Step-3 *STATIC 0.025, 1., 1E-05, 1. *BOUNDARY x1, 1, 1, 0.05 *ENDSTEP

9. Save the file and execute job geo_clay. After the job has completed view the results. Question W1–11: Sketch the initial yield surface on the meridional plane. Question W1–12: Plot the direct stress component S11 vs. the analysis time.

Compare the behavior in tension and compression with an equivalent Cap model.



W1.11

Note: Complete input files are available for your convenience. You may consult these files if you encounter difficulties following the instructions outlined here or if you wish to check your work. The input files are named geo_elastic_uni_t.inp geo_elastic_uni_c.inp geo_elastic_hydro_c.inp geo_elastic_shear.inp geo_drucker_uni_c.inp geo_drucker_uni_ct.inp geo_cap.inp geo_clay.inp and are available using the Abaqus fetch utility.



Answers 1 Material Models for Geotechnical Applications Keywords Version

Question W1–1:

Answer:

Plot the locus of stress points for uniaxial tension, uniaxial compression, hydrostatic tension, hydrostatic compression, and pure shear on the meridional plane (pressure p vs. Mises stress q). A plot of pressure p vs. Mises stress, q, can be obtained for a representative problem as follows: 1. In the Results Tree, expand the History Output container underneath the output database named geo_drucker_uni_ct.odb.odb.

2. From the list of available history data, click mouse button 3 on the Mises stress (SINV: MISES) variable. From the menu that appears, select Save As and name the X–Y data Mises. 3. Similarly, select the pressure (SINV: PRESS) variable. Save the X–Y data as Press. 4. In the Results Tree, double-click XYData. 5. In the Create XY Data dialog box, choose Operate on XY data as the source and click Continue. 6. In the Operate on XY Data dialog box, select combine(X,X) from the list of operators. Select Press and click Add to Expression. Repeat this for Mises. The final expression is: combine ("Press", "Mises")

7. Click Save As to save the Mises vs. Pressure data for this model. Repeat this procedure for all five stress states. Plot all five Mises vs. Pressure curves in one X–Y plot (expand the XYData container, select the curves, and click mouse button 3; from the menu that appears, select Plot). Customize the plot (double-click plot components in the viewport or use Options→XY Options → …).



WA1.2

Shear

UC

UT

HC

HT

Question W1–2:

Derive the relationship between the pressure and Mises stress for uniaxial tension. Compare with the plot in the meridional plane.

Answer:

For uniaxial tension we have σ11 = σt and all other stress components as zero. Pressure p = −(σ11 + σ22 + σ33)/3 = −σt /3 Deviatoric stress component

S11

= σt + p = 2σt /3

S22

= −p

S33 = −p All other deviatoric stress components are zero. Mises stress

= {(3/2) S:S }0.5 = σt

Thus for uniaxial tension, Mises stress is related to the pressure by a factor of −3. The plot in the meridional plane agrees with this relationship.



WA1.3

Question W1–3:

For Step-1, sketch the yield surface on the meridional plane and plot the stress history along the loading direction. Confirm that the stress is in equilibrium with the external load.

Answer:

The Mises stress, q, has the same magnitude as the direct stress σ11, but it is always positive, i.e. q = |σ11|. If we apply an external pressure of 25 MPa on a face along the global 1-direction, we find that σ11 = −25 MPa. As this state of stress is uniaxial compression, we have Mises, q, as 25 MPa, and pressure, p, as 25/3, that is, 8.333 MPa. The yield point in compression, σc, is 30 MPa (same as the input value). The yield point in tension, σt, has a lower value than the yield point in compression. The yield stress in tension, σt , is related to the yield stress in compression, σc , as follows: We have the yield criterion as: F = q – p tan β − d = 0 For uniaxial tension, q = σt and pressure p = −σt /3. For uniaxial compression, q = σc and pressure p = σc /3. Substituting these values in the expression for F and eliminating d, we get,

σt = σc {(1−(tan β)/3)/(1+(tan β)/3)}. Therefore, for σc = 30 MPa and β = 30 degrees, we get,

σt = 20.32 MPa. The yield surface and the stress history for Step 1 can be depicted as follows:

Yield surface Mises q

(10, 30) (25/3 , 25)

(−20.32/3, 20.32)

Stress history

(0, 0)


Pressure p


WA1.4

Question W1–4: Why does the job not complete Step-2? Answer: The analysis fails with a numerical singularity, because the

applied load is larger than the limit (or failure) yield stress in compression. Question W1–5:

Answer:

Plot the history of direct stress, (S11), the Mises stress, (q), and the stress locus on the meridional plane. Does the material yield at the correct yield stress in tension and compression? Yes, the material does yield at the expected stress values in compression and tension. The yield stress in compression is 30 MPa and the yield stress in tension is 20.32 MPa, as discussed in the answer to Question W1–3. History plots of direct stress, (S11), the Mises stress, (q), and the stress locus on the meridional plane follow:

Mises

σt = 20.32 MPa

|σc| = 30.0 MPa

S11



WA1.5

30.0 MPa

Compression 20.32 MPa

Tension

In this plot the unloading curve appears to deviate from the loading curve only because of the coarseness of the incrementation.

Question W1–6:


Answer:

The volume strain is zero when the dilation angle is zero; i.e. the material is incompressible. Hybrid elements may be required to prevent locking. The slope of the line in a plot of equivalent deviatoric plastic strain vs volumetric plastic strain is the tangent of the dilation angle. The Drucker-Prager model predicts continuing dilation, regardless of the mode of deformation. In practice, materials show a limited amount of dilation.



WA1.6

Dilation angle = 30

Dilation angle = 0

Question W1–7:


Answer:

Yield surface

Flow potential

Mises q (0, 24.226)

β = 30o Pressure p


(30, 0)


WA1.7

Question W1–8: Plot the equivalent deviatoric plastic strain as a function of the

Answer:

volumetric plastic strain. (Use the Python script ws_geotech_dev_tensor.py as described previously.) Compare the result with the behavior obtained with the Drucker-Prager model in Question W1–6 and explain the behavior of the Cap model. The Cap model exhibits a limit on the amount of dilation. The Drucker-Prager model exhibits no dilation limit.

Question W1–9: Sketch the updated flow potential on the meridional plane. Answer:

Flow potential

Mises q (0, 24.226)

Yield surface (10, 30)

β = 30o (13, 0)


Pressure p


WA1.8

Question W1–10: Plot the stress component S11 vs. the strain component E11.

Explain the response. Yielding on the shear surface results in perfect plastic response; i.e. the shear surface is a failure surface.

Answer:

Question W1–11: Sketch the initial yield surface on the meridional plane. Answer:

Critical state line

Mises q

Yield surface Slope 3 Tension

Compression

Μ = 30o Pressure p


(30, 0)


WA1.9

Question W1–12: Plot the direct stress component S11 vs. the analysis time.

Answer:

Compare the behavior in tension and compression with an equivalent Cap model. In the Clay model, the yielding results in softening on the dry side, while hardening takes place on the wet side. The failure state occurs when the stress reaches the critical state line. The Clay model exhibits softening during tensile failure and hardening in compressive failure. The Cap model, however, would predict perfect plastic failure on the failure surface. If the shear failure surface of the Cap model matches with the critical state line of this Clay model, then the Cap model would predict failure without any softening at a relative uniaxial tensile stress of 14.53 MPa (the tension limit in the figure above, that is, at S11= −15.47 MPa). In the compressive case, it will fail without any hardening at a relative uniaxial compressive stress of 21.45 MPa (the compression limit in the figure above, that is, at S11= −51.45 MPa).

Step-1 Hydrostatic compression © Dassault Systèmes, 2008

Step-2 Uniaxial compression

Step-3 Uniaxial tension


Workshop 2 Pore Fluid Flow Analysis: Consolidation Keywords Version Note: This workshop provides instructions on using the pore fluid flow modeling feature in Abaqus in terms of the Abaqus Keywords interface. If you wish to use the Abaqus GUI interface instead, please see the “Interactive” version of these instructions. Please complete either the Keywords or Interactive version of this workshop.


Model pore fluid flow in fully saturated problems. Specify initial conditions that equilibriate loading. Use Abaqus/Viewer to view and interpret the pore fluid flow phenomenon.

In this workshop you will perform a consolidation analysis on a fully saturated model. This problem illustrates how consolidation analyses can be performed in Abaqus by using the available coupled stress diffusion analysis cabability. The one-dimensional Terzaghi consolidation problem will be modeled using two dimensional plane strain elements. In this problem a body of soil is confined by impermeable, smooth, rigid walls on all but the top surface where perfect drainage is possible. At the start of the analysis, a load is applied suddenly on the top surface of the body. We wish to predict the response of the soil as a function of time after the load application. Gravity is neglected. This problem is similar to “The Terzaghi Consolidation Problem,” Section 1.14.1 of the Abaqus Benchmarks Manual. The model is 50 m wide and 100 m high. The soil is considered to be elastic with a modulus of 1.0e8 Pa, Poisson’s ratio of 0.3, and permeability of 0.002 m/s. The initial spatially constant void ratio value is 1.1. The soil is considered to be fully saturated. The model is meshed with 10 elements through its height and 1 element across its width. The meshed model is shown in Figure W2–1.



W2.2

Figure W2–1 Consolidation specimen mesh.

Preliminaries 1. Enter the working directory for this workshop: ../geotechnical/keywords/consolidate

This directory contains the input file geo_consolidate.inp. 2. Open geo_consolidate.inp in a text editor. The input file geo_consolidate.inp already contains the model mesh, material properties, and displacement boundary conditions. The model does not have any pore pressure boundary conditions specified on the exterior edges of the soil. This is the natural zero pore fluid flux boundary condition, which represents an undrained condition. The application of an external load will result in a build up of pore pressure in the model. You will simulate this phenomenon in the first step of the consolidation analysis, wherein you will apply a pressure load on the top surface of the model and not allow for any pore fluid drainage across the whole model boundary. In the second step of the analysis you will maintain the pressure load from the first step and allow for pore fluid drainage across the top surface of the model.



W2.3

Initial Void Ratio In this section you will specify an initial void ratio of 1.1 as a pre-defined initial field. Use the *INITIAL CONDITIONS option to specify the initial void ratio for element set all, which encompasses all the elements in the model. Add the following lines before the *STEP option in order to specify the initial void ratio: *Initial conditions, type=ratio all, 1.1, 0.0, 1.1, 1.0

Pressure loading The first analysis step already exists in the input file. You must complete the step definition by defining the pressure load. A surface named surfload has already been created in the input file; this surface covers the top boundary of the consolidation specimen. To apply a pressure load of 100.0 Pa to the surface surfload, add the following lines before the output requests in the first analysis step: *Dsload surfload, P, 100.0 Question W2–1:

What will be the state of effective stress in the model when a pressure load of 100.0 is applied on the top surface of the model and pore fluid drainage is prevented from taking place?

Consolidation In this section you will define the second step with a total step time of 10.0. In this step you will allow for pore fluid drainage across the top boundary by specifying appropriate pore pressure boundary conditions on the top surface. Use a prescribed pore pressure value of 0.0 for this purpose. A node set named porebc, which contains both the nodes on the top surface, has been created for the application of the pore pressure boundary conditions. You also need to specify an appropriate value for UTOL, which is a parameter that governs the accuracy of the time integration for the transient consolidation analysis. Question W2–2:


Determine an appropriate value for UTOL.


W2.4

The element size and the time increment size are related to the extent that time increments smaller than a certain size give no useful information. This coupling of the spatial and temporal approximations is most obvious at the start of diffusion problems, immediately after prescribed changes in the boundary values. Question W2–3:

Using the time-space criterion described in Lecture 5, determine the initial time increment size to avoid overshoot in the initial solution.

The values of UTOL and the initial time increment size will be specified as data for the second step of the analysis. Add the following lines after the *END STEP option of the first step and use the appropriate values for utol and the initial time increment: *Step, name=Step-2 *Soils, consolidation, utol=utol Initial time increment, 10., 1e-05, 10., *Boundary porebc, 8, 8, 0.0 *End Step

Save the input file geo_consolidate.inp and submit it for anlaysis.

Visualizing the results After the job completes, view the results in Abaqus/Viewer. Plot the pore pressure and vertical stress along the vertical edge at different times during the analysis. Question W2–5: Plot the pore pressure and vertical stress as a function of time at the bottom of the specimen. Question W2–4:

Note: The complete input file is available for your convenience. You may consult this file if you encounter difficulties following the instructions outlined here or if you wish to check your work. The input file is named geo_consolidate_complete.inp and is available in the working directory for this workshop.



Answers 2 Pore Fluid Flow Analysis: Consolidation Keywords Version

Question W2–1:

Answer:


What will be the state of effective stress in the model when a pressure load of 100.0 Pa is applied on the top surface of the model and pore fluid drainage is prevented from taking place? A stress state with a uniform pore pressure equal to the external load throughout the body with no effective stress carried by the soil skeleton is obtained. The fluid carries the instantaneous load.

Determine an appropriate value for UTOL. UTOL which is one order of magnitude less than the applied pressure should give moderate accuracy.

Question W2–3: Using the time-space criterion described in Lecture 5,

Answer:

determine the initial time increment size to avoid overshoot in the initial solution. The initial time increment can be computed as follows: ∆t ≥ {γw/(6Ek)}(∆h)2 We have ∆h = 10 m, E = 108 Pa, k = 2×10-4 m/s, N/m3. Hence,

γw =1.0

∆t ≥ 0.00083 s.



WA2.2


Plot the pore pressure and vertical stress along the vertical edge at different times during the analysis. A path needs to be created in Abaqus/Viewer so that results along this path can be plotted. The path creation and plotting of X-Y results along this path can be accomplished as follows: 1. In the Results Tree, double-click Paths. 2. In the Create Path dialog box, select Node list and click Continue. 3. In the Edit Node List Path dialog box, click Add After. 4. In the viewport, select all consecutive nodes on one vertical edge from the bottom to the top. Click Done in the prompt area. 5. The selected nodes are listed in the Edit Node List Path dialog box. Click OK to complete the path definition. 6. In the Results Tree, double-click XYData. 7. In the Create XY Data dialog box, select Path and click Continue. 8. In the XY Data from Path dialog box, select the frame and the field output variable by clicking Step/Frame and Field Output respectively. Plot and save stress (S22) and pore pressure (POR) X-Y data for a few frames. 9. Expand the XYData container. 10. Select the saved data names and click mouse button 3; from the menu that appears, select Plot. A combined X-Y plot for the selected curves will be displayed. The plots indicate that the pore pressure decreases and the values of S22 increase over time. The pore pressure values are almost zero after about t=1.25.



WA2.3

t = 0.001 t = 0.022

t = 0.251

t = 1.248

t = 0.001 t = 0.022

t = 0.251

t = 1.248



WA2.4


Plot the pore pressure and vertical stress as a function of time at the bottom of the specimen. The pore pressure and stress results were written to the output database as history data during the analysis. This data is accessible from the History Output dialog box (Result→ History Output). To plot the stress magnitude you must save the history data and operate on it (double-click XYData, and select Operate on XY data as the data source). The values of POR decrease over time as the pore water permeates away from the top surface. The values of S22 simultaneously increase in magnitude. Hence, as the consolidation progresses, the load initially carried by the pore fluid is transferred to the solid grains.

S22 magnitude

POR



Workshop 3 Pore Fluid Flow Analysis: Wicking Keywords Version Note: This workshop provides instructions on using the pore fluid flow modeling feature in Abaqus in terms of the Abaqus Keywords interface. If you wish to use the Abaqus GUI interface instead, please see the “GUI” version of these instructions. Please complete either the Keywords or Interactive version of this workshop.


Model pore fluid flow in partially saturated problems. Specify initial conditions that equilibriate loading. Use Abaqus/Viewer to view and interpret the pore fluid flow phenomenon.

In this workshop you will perform a pore fluid absorption analysis on a partially saturated specimen. This problem illustrates how pore fluid absorption analyses can be performed in Abaqus using the available coupled stress diffusion analysis cabability. We consider a one-dimensional wicking test where the absorption of fluid takes place against the gravity load caused by the weight of the fluid. In this test the fluid is made available at the bottom of a partially saturated material column, and the material absorbs as much fluid as the weight of the rising fluid permits. The column of material is kinematically constrained in the horizontal direction so that all material deformations will be in the vertical direction. In this sense the problem is one-dimensional. The model mesh is shown in Figure W3–1.



W3.2

Figure W3–1 Wicking model mesh.

Preliminaries 1. Enter the working directory for this workshop: ../geotechnical/keywords/wicking

The directory contains the input file geo_wicking.inp. 2. Open geo_wicking.inp in a text editor. The input file geo_wicking.inp already contains the model mesh, material properties, and loading history. The units used for the material data in the database are in terms of metric tonnes, meters (m), and seconds (s), and therefore the force will be in kiloNewtons (kN) and stress in kiloPascals (kPa). The column of material is 1.0 m high and 0.1 m wide. It is meshed with 10 plane strain elements with displacement and pore pressure degrees of freedom at the nodes (CPE4P). For the mechanical properties we assume that the material is elastic with a Young's modulus of 50000 Pa and a Poisson's ratio of 0.0. The dry mass density of the material is 100 kg/m3.



W3.3

The permeability of the fully saturated material is 3.7 × 10−4 m/sec. We assume that the permeability varies as a cubic function of saturation, which is the default setting for partially saturated permeability. The specific weight of the fluid (water) is 104 N/m3 and the bulk modulus is 2 GPa. The capillary action in the porous medium is defined by absorption/exsorption curves shown in Figure W3–2. These curves give the (negative) pore pressure versus saturation relationships for the absorption and exsorption behaviors. The transition between absorption and exsorption and vice-versa takes place along a scanning slope which in this example is set by default to be 1.05 times the largest slope of any branch in the absorption/exsorption curves.

Figure W3–2 Absorption and exsorption curves for the model material.



W3.4

Initial conditions The initial condition for saturation is 5% and the initial void ratio is 5.0 throughout the column. The initial conditions for the pore pressure must have a gradient that is equal to the specific weight of the fluid so that, according to Darcy's law, there is no initial flow. For this purpose we assume the initial pore pressures vary linearly from −12000 Pa at the bottom of the column to −22000 Pa at the top of the column. These initial conditions satisfy the pore pressure/saturation relationship in the sense that they are between the absorption and exsorption curves. Question W3–1:

Determine the vertical coordinate of the phreatic surface.

The initial conditions for the effective stress are calculated from the densities of the dry material and the fluid, the initial saturation, the void ratio, and the initial pore pressures using equilibrium considerations and the effective stress principle. It is important to specify the correct initial conditions for this type of problem; otherwise, the system may be so far out of equilibrium initially that it may fail to start because converged solutions cannot be found. Question W3–2:

Determine the initial conditions for the effective stress.

The initial conditions for the pore pressure, effective stress, saturation, and void ratio must be specified as input data for the analysis. An element set named all and a node set named all have been created encompassing all the elements and nodes in the model. Add the following lines before the *STEP option in order to specify the initial conditions (You will need to type in the appropriate values for the quantities in italics. Recall that the model is defined in terms of metric tonnes, meters, and seconds.): *INITIAL CONDITIONS,TYPE=SATURATION all,.05 *INITIAL CONDITIONS,TYPE=PORE PRESSURE all,


1.0,


0.0

*INITIAL CONDITIONS,TYPE=RATIO all,5. *INITIAL CONDITIONS,TYPE=STRESS,GEOSTATIC all,



1.0,


0.0, 0., 0.


W3.5

Analysis and results The input file is now ready for execution. The analysis will be performed in two steps. The first step is a geostatic analysis step with a short time duration. This geostatic step ensures that the model is in equilibrium at the initiation of step 2, which models the transient wicking phenomenon. For step 2, the convergence criterion for the ratio of the largest solution correction to the largest corresponding incremental solution value has been increased from the default value of 0.01 to a value of 1.0. This increase helps in obtaining convergence in the early stages of the analysis wherein the rate of change of saturation and pore pressure is large at the bottom of the wicking column. Steps 1 and 2 are already included in the input file geo_wicking.inp. Save the input file geo_wicking.inp and submit it for analysis. After the job completes, view the results in Abaqus/Viewer. Question W3–3:

Plot the vertical stress distribution, the pore pressure, and saturation along the height of the sample at different times.

Question W3–4:

Plot the time history of the volume of the fluid absorbed by the column.

Note: The complete input file is available for your convenience. You may consult this file if you encounter difficulties following the instructions outlined here or if you wish to check your work. The input file is named geo_wicking_complete.inp and is available in the working directory for this workshop.



Answers 3 Pore Fluid Flow Analysis: Wicking Keywords Version

Question W3–1: Determine the vertical coordinate of the phreatic surface. Answer:

The phreatic surface lies at a vertical coordinate of −1.2 m. This value can be found by linearly extrapolating the pore pressure distribution to obtain the vertical coordinate at which the pore pressure is zero. The known pore pressure values are −22000 Pa, and −12000 Pa, at vertical coordinates of 1.0 m, and 0.0 m, respectively.

Question W3–2: Determine the initial conditions for the effective stress. Answer:

The initial void ratio is spatially constant and is 5.0. The initial saturation is spatially constant and is 0.05. Denoting γd as the dry weight density of the soil, γw as the weight density of the pore fluid, e as the void ratio, n as the porosity, s as the saturation, z as the vertical coordinate, σz as the total stress, and σ z as the effective stress, we have: (dσz/dz) = γd + snγw, which implies,

σz = (γd + snγw)(z−z0). For γd = 103, s = 0.05, n = e/(1+e) = 5/6, and z0 = 1, σz = 1416.7 (z−1). Effective stress σ z = σz + suw At z = 0.0, pore pressure uw = −12.0 kPa At z = 1.0, pore pressure uw = −22.0 kPa Therefore, At z = 0.0, effective stress σ z= −2.0167 kPa At z = 1.0, effective stress σ z= −1.1000 kPa



WA3.2


Plot the vertical stress distribution, the pore pressure, and saturation along the height of the sample at different times. At steady state the pore pressure gradient must equal the weight of the fluid and hence the pore pressure varies linearly with height and saturation. Thus, points close to the bottom of the column are fully saturated, while those at the top are still at 5% saturation. In the following plots the stress and pore pressure values are in kiloPascals:

t=0

t = 15 t = 171 t = 2000 t = 69526 t = 200000




WA3.3

t = 200000 t = 15 t = 69526 t=0

t = 171 t = 2000


t = 15 t = 171 t = 2000 t = 69526 t = 200000

t=0




WA3.4

Question W3–4: Plot the time history of the volume of the fluid absorbed by the

column. Answer:

The volume of the fluid absorbed by the column will be given by the sum of the values of the output variable RVT (reactive fluid total volume) at all boundary nodes where the pore pressure is prescribed. The pore pressure is prescribed at nodes on the bottom of the column. The following plot, showing the total RVT value in cubic meters, indicates that the rate of increase in RVT decreases with time and the total volume of the absorbed fluid approaches a limiting value.



Geotech Book (2)

Recommend Documents