Manzari Dafalias Material

• Command_Manual
  • Tcl Commands
  • Modeling_Commands
    • model
    • uniaxialMaterial
    • ndMaterial
    • frictionModel
    • section
    • geometricTransf
    • element
    • node
    • sp commands
      • fix
      • fixX
      • fixY
      • fixZ
    • mp commands
      • equalDOF
      • rigidDiaphragm
      • rigidLink
    • timeSeries
    • pattern
    • mass
    • block commands
      • block2D
      • block3D
    • region
    • rayleigh
  • Analysis Commands
    • constraints Command
    • numberer Command
    • system Command
    • test Command
    • algorithm Command
    • integrator Command
    • analysis Command
    • eigen Command
    • analyze Command
  • Output Commands
    • recorder
    • record
    • print
    • printA
    • logFile
    • RealTime_Output_Commands
  • Misc Commands
    • eleResponse Command
    • getTime Command
    • getEleTags Command
    • getNodeTags
    • loadConst Command
    • nodeCoord Command
    • nodeDisp Command
    • nodeVel Command
    • nodeAccel Command
    • nodeEigenvector Command
    • nodeBounds
    • print Command
    • remove Command
    • reset Command
    • testIter
    • testNorms
    • wipe Command
    • exit Command
  • DataBase Commands
    • Template:DataBaseCommandsMenu

---

This command is used to construct a multi-dimensional Manzari-Dafalias(2004) material.

nDmaterial ManzariDafalias $matTag $G0 $nu $e_init $Mc $c $lambda_c $e0 $ksi $P_atm $m $h0 $ch $nb $A0 $nd $z_max $cz $Den <$intScheme $TanType $JacoType $TolF $TolR>

$matTag	integer tag identifying material
$G0	bulk modulus constant
$nu	poisson ratio
$e_init	initial void ratio
$Mc	critical state stress ratio
$c	ratio of critical state stress ratio in extension and compression
$lambda_c	critical state line constant
$e0	critical void ratio at p = 0
$ksi	critical state line constant
$P_atm	atmospheric pressure
$m	yield surface constant (radius of yield surface in stress ratio space)
$h0	constant parameter
$ch	constant parameter
$nb	bounding surface parameter, $nb ≥ 0
$A0	dilatancy parameter
$nd	dilatancy surface parameter $nd ≥ 0
$z_max	fabric-dilatancy tensor parameter
$cz	fabric-dilatancy tensor parameter
$Den	mass density of the material
<$intScheme>	plastic integration scheme, 0: forward explicit, 1: backward implicit, 2: improved backward implicit (default = 2)
<$TanType>	material modulus matrix, 0: elastic stiffness, 1: continuum elastoplastic stiffness, 2: consistent elastoplastic stiffness (default = 2)
<$JacoType>	jacobian matrix used for newton iterations, 0: finite difference jacobian, 1: analytical jacobian (default = 1)
<$TolF>	drift from yield surface tolerance (default = 1.0e-7)
<$TolR>	newton iterations residuals tolerance (default = 1.0e-7)

The material formulations for the Manzari-Dafalias object are "ThreeDimensional" and "PlaneStrain"

---

Code Developed by: Alborz Ghofrani, Pedro Arduino, U Washington

---

Notes

• Valid Element Recorder queries are
  • stress, strain
  • alpha (or backstressratio) for <math>\mathbf{\alpha}</math>
  • fabric for <math>\mathbf{z}</math>
  • alpha_in (or alphain) for <math>\mathbf{\alpha_{in}}</math>

e.g.
 recorder Element -eleRange 1 $numElem -time -file stress.out  stress

• Elastic or Elastoplastic response could be enforced by

Elastic:	updateMaterialStage -material $matTag -stage 0
Elastoplastic:	updateMaterialStage -material $matTag -stage 1

• Integration scheme ($intScheme):

0:	Forward Euler integration scheme with no error control on the yield surface drift
1:	Backward Euler integration scheme (This is the same as Closest Point Projection Method), Newton-Raphson algorithm is used to solve the implicit equation. May not converge due to sensitivity of Newton-Raphson's method to the initial guess.
2:	The same as 'option 1' but in case of failure in Newton-Raphson, it uses different strategies to insure convergence.

• Jacobian ($jacoType):

0:	Uses a finite difference jacobian <math> \mathbf{J}(x) = \frac{\mathbf{R}(x+h) - \mathbf{R}(x)}{h} </math>
1:	Uses analytical jacobian <math> \mathbf{J}(x) = \frac{\partial{\mathbf{R}}}{\partial{x}} </math>

Theory

<math> p = \frac{1}{3} \mathrm{tr}(\mathbf{\sigma}) </math>

<math> \mathbf{s} = \mathrm{dev} (\mathbf{\sigma}) = \mathbf{\sigma} - \frac{1}{3} p \mathbf{1}

</math>

Elasticity

Elastic moduli are considered to be functions of p and current void ratio:

<math> G = G_0 p_{atm}\frac{\left(2.97-e\right)^2}{1+e}\left(\frac{p}{p_{atm}}\right)^{1/2}

</math>

<math> K = \frac{2(1+\nu)}{3(1-2\nu)} G </math>

The elastic stress-strain relationship is:

<math> d\mathbf{e}^\mathrm{e} = \frac{d\mathbf{s}}{2G}

</math>

<math> d\varepsilon^\mathrm{e}_v = \frac{dp}{K}

</math>

Critical State Line

A power relationship is assumed for the critical state line:

<math> e_c = e_0 - \lambda_c\left(\frac{p_c}{p_{atm}}\right)^\xi </math>

where <math>e_0</math> is the void ratio at <math> p_c = 0 </math> and <math> \lambda_c </math> and <math> \xi </math> constants.

Yield Surface

Yield surface is a stress-ratio dependent surface in this model and is defined as

<math> \left\| \mathbf{s} - p \mathbf{\alpha} \right\| - \sqrt\frac{2}{3}pm = 0 </math>

with <math> \mathbf{\alpha} </math> being the deviatoric back stress-ratio.

Plastic Strain Increment

The increment of the plastic strain tensor is given by

<math> d\mathbf{\varepsilon}^p = \langle L \rangle \mathbf{R}

</math>

where

<math> \mathbf{R} = \mathbf{R'} + \frac{1}{3} D \mathbf{1} </math>

therefore

<math> d\mathbf{e}^p = \langle L \rangle \mathbf{R'}</math> and <math> d\varepsilon^p_v = \langle L \rangle D </math>

The hardening modulus in this model is defined as

<math> K_p = \frac{2}{3} p h (\mathbf{\alpha}^b_{\theta} - \mathbf{\alpha}): \mathbf{n}

</math> where <math>\mathbf{n}</math> is the deviatoric part of the gradient to yield surface.

<math> \mathbf{\alpha}^b_{\theta} = \sqrt{\frac{2}{3}} \left[g(\theta,c) M_c exp(-n^b\Psi) - m\right] \mathbf{n}

</math>, <math>\Psi</math> being the state parameter.

the hardening parameter <math> h </math> is defined as

<math> h = \frac{b_0}{(\mathbf{\alpha}-\mathbf{\alpha_{in}}):\mathbf{n}}

</math>, <math>\mathbf{\alpha_{in}}</math> is the value of <math>\mathbf{\alpha}</math> at initiation of loading cycle.

<math>b_0 = G_0 h_0 (1-c_h e) \left(\frac{p}{p_{atm}}\right)^{-1/2} </math>

Also the dilation parameters are defined as

<math> D = A_d (\mathbf{\alpha}^d_{\theta}-\mathbf{\alpha}) : \mathbf{n}

</math>

<math> \mathbf{\alpha}^d_{\theta} = \sqrt{\frac{2}{3}} \left[g(\theta,c) M_c exp(n^d\Psi) - m\right] \mathbf{n}

</math>

<math> A_d = A_0 (1+\langle \mathbf{z : n}\rangle) </math>, where <math> \mathbf{z} </math> is the fabric tensor.

The evolution of fabric and the back stress-ratio tensors are defined as

<math> d\mathbf{z} = - c_z \langle -d\varepsilon^p_v \rangle (z_{max}\mathbf{n}+\mathbf{z}) </math>

<math> d\mathbf{\alpha} = \langle L \rangle (2/3) h (\mathbf{\alpha}^b_{\theta} - \mathbf{\alpha})

</math>

Example

This example, provides an undrained confined triaxial compression test using one 8-node SSPBrichUP element and ManzariDafalias material model.

# HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH #
# 3D Undrained Conventional Triaxial Compression Test Using One Element #
# University of Washington, Department of Civil and Environmental Eng   #
# Geotechnical Eng Group, A. Ghofrani, P. Arduino - Dec 2013            #
# Basic units are m, Ton(metric), s										#
# HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH #

wipe

# ------------------------ #
# Test Specific parameters #
# ------------------------ #
# Confinement Stress
set pConf -300.0
# Deviatoric strain
set devDisp -0.3
# Permeablity
set perm 1.0e-10
# Initial void ratio
set vR 0.8

# Rayleigh damping parameter
set damp   0.1
set omega1 0.0157
set omega2 64.123
set a1 [expr 2.0*$damp/($omega1+$omega2)]
set a0 [expr $a1*$omega1*$omega2]

# HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
# HHHHHHHHHHHHHHHHHHHHHHHHHHHCreate ModelHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
# HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH

# Create a 3D model with 4 Degrees of Freedom
model BasicBuilder -ndm 3 -ndf 4

# Create nodes
node 1	1.0	0.0	0.0
node 2	1.0	1.0	0.0
node 3 	0.0	1.0	0.0	
node 4	0.0	0.0	0.0
node 5	1.0	0.0	1.0
node 6 	1.0	1.0	1.0
node 7 	0.0	1.0	1.0
node 8 	0.0	0.0	1.0

# Create Fixities
fix 1 	0 1 1 1
fix 2 	0 0 1 1
fix 3	1 0 1 1
fix 4 	1 1 1 1
fix 5	0 1 0 1
fix 6 	0 0 0 1
fix 7	1 0 0 1
fix 8 	1 1 0 1


# Create material
#          ManzariDafalias  tag    G0   nu   e_init   Mc    c    lambda_c    e0    ksi   P_atm   m    h0   ch    nb  A0      nd   z_max   cz    Den   intScheme  TanType  JacoType   TolF     TolR
nDMaterial ManzariDafalias   1    125  0.05   $vR    1.25  0.712   0.019    0.934  0.7    100   0.01 7.05 0.968 1.1 0.704    3.5    4     600  1.42      2          2         1     1.0e-10 1.0e-10

# Create element
#       SSPbrickUP  tag    i j k l m n p q  matTag  fBulk  fDen    k1    k2   k3   void   alpha    <b1 b2 b3>
element SSPbrickUP   1     1 2 3 4 5 6 7 8    1     2.2e6   1.0  $perm $perm $perm  $vR   1.5e-9 

# Create recorders
recorder Node    -file disp.out   -time -nodeRange 1 8 -dof 1 2 3 disp
recorder Node    -file press.out  -time -nodeRange 1 8 -dof 4     vel
recorder Element -file stress.out -time stress
recorder Element -file strain.out -time strain
recorder Element -file alpha.out  -time alpha
recorder Element -file fabric.out -time fabric


# Create analysis
constraints Penalty 1.0e18 1.0e18
test        NormDispIncr 1.0e-5 20 1
algorithm   Newton
numberer    RCM
system      BandGeneral
integrator  Newmark 0.5 0.25
rayleigh    $a0 0. $a1 0.0
analysis    Transient

# Apply confinement pressure
set pNode [expr $pConf / 4.0]
pattern Plain 1 {Series -time {0 10000 1e10} -values {0 1 1} -factor 1} {
    load 1  $pNode  0.0    0.0    0.0
    load 2  $pNode  $pNode 0.0    0.0
    load 3  0.0     $pNode 0.0    0.0
    load 4  0.0     0.0    0.0    0.0
    load 5  $pNode  0.0    $pNode 0.0
    load 6  $pNode  $pNode $pNode 0.0
    load 7  0.0     $pNode $pNode 0.0
    load 8  0.0     0.0    $pNode 0.0
}
analyze 100 100

# Let the model rest and waves damp out
analyze 50  100

# Close drainage valves
for {set x 1} {$x<9} {incr x} {
   remove sp $x 4
}
analyze 50 100

# Read vertical displacement of top plane
set vertDisp [nodeDisp 5 3]
# Apply deviatoric strain
set lValues [list 1 [expr 1+$devDisp/$vertDisp] [expr 1+$devDisp/$vertDisp]]
set ts "{Series -time {20000 1020000 10020000} -values {$lValues} -factor 1}"

# loading object deviator stress
eval "pattern Plain 2 $ts { 
	sp 5  3	$vertDisp
	sp 6  3	$vertDisp
	sp 7  3 $vertDisp
	sp 8  3 $vertDisp
}"

# Set number and length of (pseudo)time steps
set dT      100
set numStep 10000

# Analyze and use substepping if needed
set remStep $numStep
set success 0
proc subStepAnalyze {dT subStep} {
	if {$subStep > 10} {
		return -10
	}
	for {set i 1} {$i < 3} {incr i} {
		puts "Try dT = $dT"
		set success [analyze 1 $dT]
		if {$success != 0} {
			set success [subStepAnalyze [expr $dT/2.0] [expr $subStep+1]]
			if {$success == -10} {
				puts "Did not converge."
				return success
			}
		} else {
			if {$i==1} {
				puts "Substep $subStep : Left side converged with dT = $dT"
			} else {
				puts "Substep $subStep : Right side converged with dT = $dT"
			}
		}
	}
	return success
}

puts "Start analysis"
set startT [clock seconds]

while {$success != -10} {
	set subStep 0
	set success [analyze $remStep  $dT]
	if {$success == 0} {
		puts "Analysis Finished"
		break
	} else {
		set curTime  [getTime]
		puts "Analysis failed at $curTime . Try substepping."
		set success  [subStepAnalyze [expr $dT/2.0] [incr subStep]]
        set curStep  [expr int(($curTime-20000)/$dT + 1)]
        set remStep  [expr int($numStep-$curStep)]
		puts "Current step: $curStep , Remaining steps: $remStep"
	}
}
set endT [clock seconds]
puts "loading analysis execution time: [expr $endT-$startT] seconds."

wipe

References

Dafalias YF, Manzari MT. "Simple plasticity sand model accounting for fabric change effects". Journal of Engineering Mechanics 2004