Dafalias Manzari Model

[avytin]

I am performing a plane strain shearing of one brick (stdBrick and Brick8N) element using Dafalias Manzari model. The problem is that the shear induced volumetric strains that I see are fully recoverable at every cycle. Is there any a mistake in my .tcl file or is there a specific flag that I have to activate to get plasticity?

Thanks in advance!!

[avytin]

model BasicBuilder -ndm 3 -ndf 3


node 1 0.000 0.000 0.000
node 2 1.000 0.000 0.000
node 3 0.000 1.000 0.000
node 4 1.000 1.000 0.000
node 5 0.000 0.000 1.000
node 6 1.000 0.000 1.000
node 7 0.000 1.000 1.000
node 8 1.000 1.000 1.000

# Define the three dimensional material
## Dafalias-Manzari (2004), with fabric dilation
set rho 1.5
set e0 0.74
set G0 125.0
set v 0.05
set Pat 100.0
set kc 0.8
set M_cal 1.25
set c 0.712
set lambda_c 0.019
set xi 0.7
set er 0.934
set m_low 0.01
set h0 7.05
set ch 0.968
set nb 1.1
set A0 0.704
set nd 3.5
set z_max 4.0
set cz 600.0
set zT "0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0"
set p 0.0
set initS "$p 0.0 0.0 0.0 $p 0.0 0.0 0.0 $p"
set mc "$rho $e0 $G0 $v $Pat $kc $M_cal $c $lambda_c $xi $er \
$m_low $h0 $ch $nb $A0 $nd $z_max $cz"
set it "$zT $zT"
set MP "MaterialConstant 19 $mc InternalTensor 2 $it"
set EL "Dafalias-Manzari 3 4 5 6 2 $initS"
set YF "Dafalias-Manzari 0 12 2 1"
set PF "Dafalias-Manzari 0 2 0 11 0 9 0 10 0 5 0 12 0 7 0 8 0 16 0 17 2 1 2 2"
set te1 "Dafalias-Manzari 2 11 9 10 5 12 7 8 15 13 14 3 1 2"
set te2 "Dafalias-Manzari-fabric 12 19 18 1 2"
set TE "$te1 $te2"
nDMaterial NewTemplate3Dep 1 -MaterialParameter $MP \
-ElasticState $EL \
-YieldFunction $YF \
-PlasticFlow $PF \
-TensorEvolution $TE

# End of material definition

# Define Elements
set gravX 0.0
set gravY 0.
set press 0.
set bulk_f 2.2e9
set n_por 0.333
set bulk [expr $bulk_f/$n_por]

# element stdBrick 1 5 7 8 6 1 3 4 2 1 0.0 0.0 0.0 0.00
element Brick8N 1 1 3 4 2 5 7 8 6 1 0.0 0.0 0.0 0.00

fix 1 1 1 1
fix 2 1 1 1
fix 3 0 0 1
fix 4 0 0 1
fix 5 1 1 1
fix 6 1 1 1
fix 7 0 0 1
fix 8 0 0 1

equalDOF 3 4 1
equalDOF 7 8 1

# Depending on the excitation one should change the Timeseriesini factor
set Timeseriesini "Constant -factor 50."

# pattern UniformExcitation 1 1 -accel $Timeseries;
pattern Plain 1 $Timeseriesini {
load 4 0. -.5 0.
load 3 0. -.5 0.
load 7 0. -.5 0.
load 8 0. -.5 0.
}

# # Set material to elastic for gravity loading
# updateMaterialStage -material 1 -stage 1

# GRAVITY APPLICATION (elastic behavior)
set gamma 1.5

# create the SOE, ConstraintHandler, Integrator, Algorithm and Numberer
integrator Newmark $gamma [expr pow($gamma+0.5, 2)/4] \
0.00 0.0 0.002 0.0
test NormDispIncr 1.0e-5 5 1;
constraints Transformation
algorithm Newton
numberer RCM
system ProfileSPD
analysis Transient
analyze 3 5.e5
puts "End of Elastic Phase of Gravity Application"

# # Set material to elasto-plastic for the rest of the loading
# updateMaterialStage -material 1 -stage 1
analyze 5 5.e5
puts "End of Gravity Application. Starting Dynamic Excitation..."

# rezero time
wipeAnalysis
setTime 0.0
set omega 1.
set pi 3.141593

# Calculate the period from the timeseries
set T [expr 2*$pi/$omega]
set Timeseries "Sine 0.0 50. $T -shift 0.0 -factor 20."
pattern Plain 2 $Timeseries {
load 4 1. 0. 0.
load 3 1. 0. 0.
load 7 1. 0. 0.
load 8 1. 0. 0.
}

set gamma 0.6
test NormDispIncr 1.0e-4 200 1;
integrator Newmark $gamma [expr pow($gamma+0.5, 2)/4] 0.00 0.0 0.002 0.0
constraints Transformation
algorithm KrylovNewton
numberer RCM
system ProfileSPD

# create the analysis object
analysis Transient

# Recorders
recorder Node -file output_disp.txt -time -dof 1 2 3 disp
recorder Element -file stress_1.txt -time -eleRange 1 1 material 1 stress
recorder Element -file stress_2.txt -time -eleRange 1 1 material 2 stress
recorder Element -file stress_3.txt -time -eleRange 1 1 material 3 stress
recorder Element -file stress_4.txt -time -eleRange 1 1 material 4 stress
recorder Element -file strain_1.txt -time -eleRange 1 1 material 1 strain
recorder Element -file strain_2.txt -time -eleRange 1 1 material 2 strain
recorder Element -file strain_3.txt -time -eleRange 1 1 material 3 strain
recorder Element -file strain_4.txt -time -eleRange 1 1 material 4 strain

# perform analysis
set startT [clock seconds]
analyze 1500 0.01
set endT [clock seconds]
puts "Execution time: [expr $endT-$startT] seconds."

[rjaeger]

I believe the Dafalias-Manzari model uses non-associative plasticity, which means that the material tangent and stiffness matrices would be asymmetric. So the use of a ProfileSPD system might be invalid. Try switching to UmfPack.

Robbie Jaeger

[aneeman]

[quote="rjaeger"]I believe the Dafalias-Manzari model uses non-associative plasticity, which means that the material tangent and stiffness matrices would be asymmetric. So the use of a ProfileSPD system might be invalid. Try switching to UmfPack.

Robbie Jaeger[/quote]

here's a script that does 500 steps isotropic compression followed by 1000 steps triaxial shearing. The material stiffness definately loses stiffness and symmetry (jn fact some of the time its bifurcating).

Hope this helps,
alisa

#[jeremic,zhyang]@ucdavis.edu
# tcl example file using Template3Dep material
# Inelastic static analysis -- Triaxial test
#
#two load stages:
# stage 1: isotropic compression
# stage 2: triaxial shearing
#May 7, 2001

source Units.tcl

# ################################
# create the modelbuilder
# #################################

model BasicBuilder -ndm 3 -ndf 3


set g 9.81
#set ld -8.0
# ################################
# build the model
# #################################

node 1 1.0 1.0 0.0
node 2 0.0 1.0 0.0
node 3 0.0 0.0 0.0
node 4 1.0 0.0 0.0
node 5 1.0 1.0 1.0
node 6 0.0 1.0 1.0
node 7 0.0 0.0 1.0
node 8 1.0 0.0 1.0

#triaxial test boundary
fix 1 0 1 1
fix 2 1 1 1
fix 3 1 0 1
fix 4 0 0 1

set mat1 1

set rho [expr 1478.0*$kg/pow($m, 3)]; #unit weight

set rho1 $rho

# ##DM04############################
set e0 0.65
set G0 125.0
set v 0.05
set Pat 100.0
set kc 0.8
set M_cal 1.25
set c 0.712
set lambda_c 0.019
set xi 0.7
set er 0.934
set m_low 0.01
set h0 7.05
set ch 0.968
set nb 1.1
set A0 0.704
set nd 3.5
set z_max 4.0
set cz 600.0
set zZ "0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0"
set p -1.0
set initS "$p 0.0 0.0 0.0 $p 0.0 0.0 0.0 $p"
set zS "0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0"

## Specify some constants
set grvt 10
set Gr [expr -$grvt];
set poro [expr $e0/(1.0+$e0)]
set alpha 1.0
set rho_s 2.8
set rho_f 1.0
set bulk_s 1.0e12
set bulk_f 2.2e6
set kx [expr 5.0e-4/$grvt/$rho_f]
set ky [expr 5.0e-4/$grvt/$rho_f]
set kz [expr 5.0e-4/$grvt/$rho_f]
set rho [expr (1.0-$poro)*$rho_s + $poro*$rho_f ]

##-------1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
set mc "$rho $e0 $G0 $v $Pat $kc $M_cal $c $lambda_c $xi $er $m_low $h0 $ch $nb $A0 $nd $z_max $cz"
set it "$zS $zZ"
set mp "MaterialConstant 19 $mc InternalTensor 2 $it"
set el "Dafalias-Manzari 3 4 5 6 2 $initS"
set yf "Dafalias-Manzari 0 12 2 1"
set pf "Dafalias-Manzari 0 2 0 11 0 9 0 10 0 5 0 12 0 7 0 8 0 16 0 17 2 1 2 2"

set te1 "Dafalias-Manzari 2 11 9 10 5 12 7 8 15 13 14 3 1 2"
# originally, not enough parameters for openseesPDD
#set te2 "Dafalias-Manzari-fabric 12 19 18 1 2"

#alisa - using thesis values.. first 8 mimic te1
set te2 "Dafalias-Manzari-fabric 2 11 9 10 5 12 7 8 16 17 19 18 1 2"

set te "$te1 $te2"

nDMaterial NewTemplate3Dep 1 -MaterialParameter $mp \
-ElasticState $el \
-YieldFunction $yf \
-PlasticFlow $pf \
-TensorEvolution $te



#_______________tag_____8 nodes______matID__bforce1_bforce2_bforce3_rho
element Brick8N 1 5 6 7 8 1 2 3 4 1 0.0 0.0 -9.81 1.8

#stage 1
#===========================================================
#isotropic load, values changed from 5.0
set p 5.0
set np -5.0
pattern Plain 2 "Linear" {
load 1 $np 0 0 -pattern 2
load 3 0 $p 0 -pattern 2
load 4 $np $p 0 -pattern 2
load 5 $np $np $np -pattern 2
load 6 $p $np $np -pattern 2
load 7 $p $p $np -pattern 2
load 8 $np $p $np -pattern 2
}

#Set up recorder
recorder Node -file node_iso.out -time -node 5 -dof 3 disp
#recorder plot node_iso.out "disp load" 10 10 300 300 -columns 2 1


# alisa chaned NN1 from 1 to 1000
set lf1 0.01
set NN1 500
system UmfPack
constraints Penalty 1e8 1e8
test NormDispIncr 1.0e-5 30 1
integrator LoadControl $lf1 1 $lf1 $lf1
algorithm Newton
numberer RCM
analysis Static


analyze $NN1



#stage 2
#===========================================================
#Axial loading
wipeAnalysis
#loadConst

equalDOF 5 6 3 3
equalDOF 5 7 3 3
equalDOF 5 8 3 3

remove recorders
recorder domain
#set previous load constant
loadConst time 0

# Before shifting to disp. control, apply reference load on the top
pattern Plain 3 "Linear" {
load 5 0.0 0.0 $np
load 6 0.0 0.0 $np
load 7 0.0 0.0 $np
load 8 0.0 0.0 $np}

#Set up recorder for axial loading stage
recorder Node -file node_z.out -time -node 5 -dof 3 disp
#recorder plot node_z.out "disp load vertical" 10 400 300 300 -columns 1 2

recorder Element -ele 1 -file element.out force
#recorder plot element.out "nodal force displ horizontal" 10 800 300 300 -columns 1 3
#alisa changed NN2 from 10 to 1000
set NN2 1000
set ndz -0.0001

system UmfPack
constraints Penalty 1e12 1e12
test NormDispIncr 1.0e-08 30 1
integrator DisplacementControl 5 3 $ndz 1 $ndz $ndz
algorithm Newton
numberer RCM
analysis Static

analyze $NN2

[avytin]

interestingly the script you sent me with the different te2 variable does not work (DM04_z:invalid input). I've tried even the last version of the newtemplate3dep models, by downloading it through cvs.