Manzari Dafalias Material: Difference between revisions

From OpenSeesWiki
Jump to navigation Jump to search
(Created page with '{{CommandManualMenu}} This command is used to construct a multi-dimensional Manzari-Dafalias(2004) material. {| | style="background:yellow; color:black; width:800px" | '''nDma...')
 
No edit summary
 
(7 intermediate revisions by 2 users not shown)
Line 4: Line 4:


{|  
{|  
| style="background:yellow; color:black; width:800px" | '''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>'''
| style="background:yellow; color:black; width:800px" | '''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'''
|}   
|}   


Line 10: Line 10:
|  style="width:150px" | '''$matTag ''' || integer tag identifying material
|  style="width:150px" | '''$matTag ''' || integer tag identifying material
|-
|-
|  '''$G0 ''' || bulk modulus constant
|  '''$G0 ''' || shear modulus constant
|-
|-
|  '''$nu ''' || poisson ratio
|  '''$nu ''' || poisson ratio
Line 45: Line 45:
|-
|-
|  '''$Den ''' || mass density of the material  
|  '''$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)
|}
|}


Line 81: Line 71:
|-
|-
  |Elastoplastic: ||updateMaterialStage -material $matTag -stage 1
  |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>
|}
|}


Line 102: Line 76:
==Theory==
==Theory==


:::<math> p = \frac{1}{3} \mathrm{tr}(\mathbf{\sigma}) </math>
<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> \mathbf{s} = \mathrm{dev} (\mathbf{\sigma}) = \mathbf{\sigma} - \frac{1}{3} p \mathbf{1}
</math>
</math>


Line 187: Line 161:
==Example==
==Example==


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




Line 248: Line 222:


# Create material
# 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
#          ManzariDafalias  tag    G0  nu  e_init  Mc    c    lambda_c    e0    ksi  P_atm  m    h0  ch    nb  A0      nd  z_max  cz    Den   
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
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


# Create element
# Create element

Latest revision as of 02:22, 27 January 2018




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
$matTag integer tag identifying material
$G0 shear 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


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


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 SSPBrickUP 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  
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  

# 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