PDMY03 elementdriver: Difference between revisions

From OpenSeesWiki
Jump to navigation Jump to search
(Created page with '{{CommandManualMenu}} {| | style="color:black; font-size: 150%; " | <center>'''Solid-fluid fully coupled (u-p) plane-strain 9-4 noded element: saturated soil element with press...')
 
 
(7 intermediate revisions by the same user not shown)
Line 2: Line 2:


{|  
{|  
| style="color:black; font-size: 150%; " | <center>'''Solid-fluid fully coupled (u-p) plane-strain 9-4 noded element: saturated soil element with pressure dependent material, subjected to 1D sinusoidal base shaking'''</center>
| style="color:black; font-size: 150%; " | <center>'''Element Drivers for PressureDependMultiYield03 Material'''</center>
|}
|}
   
   


== Input File ==
== DSS_main.tcl file ==
<syntaxhighlight lang="tcl">
<syntaxhighlight lang="tcl">
# Written by Jinchi Lu and Zhaohui Yang (May 2004)
###########################################################################
# Model with 9-4-quad-UP element and PDMY03 material       #
# By: Arash Khosravifar  January 2009   #
# Last edit:   August 2018   #
###########################################################################


wipe
wipe;
set matOpt  1          ;# 1 = pressure depend;
wipe all;
                        ;# 2 = pressure independ;
set startT [clock seconds];
set fmass  1            ;# fluid mass density
set smass  2.0          ;# saturated soil mass density
set G    9.0e4   
set B    2.2e5   
set bulk  2.2e6                        ;#fluid-solid combined bulk modulus
set vperm  5.e-6        ;#vertical permeability (m/s)
set hperm  [expr $vperm]       ;#horizontal permeability (m/s)


set accGravity  9.81    ;#acceleration of gravity
######################### MAIN USER INPUTS #################################
set vperm  [expr $vperm/$accGravity/$fmass]  ;#actual value used in computation
set matTag 2; # 1 for (N1)60=35, 2 for (N1)60=25, 3 for (N1)60=15 and 4 for (N1)60=5
set hperm  [expr $hperm/$accGravity/$fmass]  ;#actual value used in computation
set vertPress [expr 1.0*100]; # kPa vertical confining pressure
set loadBias 0.0                        ;# Static shear load, in percentage
set loadbias 0.3; # Alpha = tau_xy/s'vo
                                        ;# of gravity load (=sin(inclination angle))
set Analysis_case "undrained_cyclic"
# Options:
# undrained_cyclic
# undrained_monotonic: define target max shear strain
# drained_cyclic: applies 2 cycles for each strain (9 shear strain ranging from 0.0003% to 3%) for the purpose of G/Gmax and Damping curves
# drained_monotonic: define target max shear strain
######################### END OF MAIN USER INPUTS ###########################


set accMul  2.          ;# acc. multiplier
# Load material properties
set period  1.0        ;# Period for applied Sine wave
source DSS_material_properties.tcl
set deltaT  0.01      ;# time step for analysis
set numSteps 2500      ;# number of time steps
set gamma    0.6        ;# Newmark integration parameter


set massProportionalDamping  0. ;
# Some variables for the ELEMENT
set InitStiffnessProportionalDamping 0.002;
set fluidDen 1.0; # Fluid mass density
set waterbulk 2.2e6; # kPa
set kdrain 1.e1; # permeability for drained loading
set kundrain 1.e-8; # permeability for undrained loading
set gravY 0.0;
set gravX 0.0;


#############################################################
# Some loading related variables
# BUILD MODEL
if {$Analysis_case == "undrained_cyclic"} {
set target_shear_strain [expr 3.0/100.]
set csr 0.45;
set period      1.;
set deltaT 0.005;
set kPerm $kundrain
}
if {$Analysis_case == "undrained_monotonic"} {
set target_shear_strain [expr 3.0/100.]
set period      1.;
set deltaT 0.005;
set numSteps 200;
set kPerm $kundrain
# also change penalty from 1e15 and 1e5 for sand to clay
}
if {$Analysis_case == "drained_monotonic"} {
# also change the damping to 0.005
set target_shear_strain [expr 1.0/100.]; #unit in strain - it will be multiplied by ySize later. Note: DON'T forget the ".0" in 1.0
set period      1.;
set deltaT 0.005;
set numSteps 3000;
 
set kPerm $kdrain
set phaseTransAng $frinctionAng
set contraction_a 0.0;
set contraction_b 0.0;
set contraction_c 0.0;
set contraction_d 0.0;
set contraction_e 0.0;
set dilation_a 0.0;
set dilation_b 0.0;
set dilation_c 0.0;
set liqParam1 0.0;
set liqParam2 0.0;
}
if {$Analysis_case == "drained_cyclic"} {
set period      10.;
set deltaT 0.005;
set numSteps 4000;
set kPerm $kdrain
set phaseTransAng $frinctionAng
set contraction_a 0.0;
set contraction_b 0.0;
set contraction_c 0.0;
set contraction_d 0.0;
set contraction_e 0.0;
set dilation_a 0.0;
set dilation_b 0.0;
set dilation_c 0.0;
set liqParam1 0.0;
set liqParam2 0.0;
}
 
# Define Rayleigh Damping
set rayleigh_damping1 0.001;
set rayleigh_damping2 0.001;
set rayleigh_w1 [expr 2*3.14/$period];
set rayleigh_w2 [expr 2*3.14/($period/2.0)];
set Kinitproprayleigh [expr 2*$rayleigh_w1*$rayleigh_w2/($rayleigh_w1*$rayleigh_w1-$rayleigh_w2*$rayleigh_w2)*(-1*$rayleigh_damping1/$rayleigh_w1+1*$rayleigh_damping2/$rayleigh_w2)]; #This is a1 [expr 2*$dampratio/(2*3.1416/0.02)];
set Mproprayleigh [expr 2*$rayleigh_w1*$rayleigh_w2/($rayleigh_w1*$rayleigh_w1-$rayleigh_w2*$rayleigh_w2)*($rayleigh_damping1*$rayleigh_w1-$rayleigh_damping2*$rayleigh_w2)]; #This is a0 [expr 2*$dampratio*(2*3.1416/2.0)];


#create the ModelBuilder
####################################################################
set xsize 0.1;
set ysize 0.1;
set thickness 0.1;
 
# define the 3DOF nodes
model basic -ndm 2 -ndf 3
model basic -ndm 2 -ndf 3
node 1 0 0
node 1   0.0 0.0  
node 2 2.5 0
node 2   $xsize 0.0  
node 3 2.5 2
node 3   $xsize $ysize
node 4 0 2
node 4   0.0 $ysize


fix 1 1 1 0  
fix 1 1 1 0
fix 2 1 1 0
fix 2 1 1 0
fix 3 0 0 1
fix 3 0 0 1
fix 4 0 0 1
fix 4 0 0 1
equalDOF 3 4 1 2
equalDOF 1 2  3
equalDOF 3 4 1 2


# define the 2DOF nodes
model basic -ndm 2 -ndf 2
model basic -ndm 2 -ndf 2
node 5 1.25 0.
node 5   [expr $xsize/2] 0.0
node 6 2.5 1
node 6   $xsize [expr $ysize/2]
node 7 1.25 2
node 7   [expr $xsize/2] $ysize
node 8 0 1
node 8   0.0 [expr $ysize/2]
node 9 1.25 1
node 9   [expr $xsize/2] [expr $ysize/2]


fix 5 1 1
fix 5 1 1
equalDOF 3 7 1 2
equalDOF 3 7 1 2
equalDOF 6 8 1 2
equalDOF 6 8 1 2
equalDOF 6 9 1 2
# equalDOF 6 9 1 2


set gravY [expr -$accGravity]  ;#calc. gravity
set gravX [expr -$gravY*$loadBias]


# define material and properties
################# define material ##################################
switch $matOpt {
nDMaterial PressureDependMultiYield03 $matTag 2 $massDen $refG $refB $frinctionAng \
   1 {
$peakShearStrain $refPress $pressDependCoe $phaseTransAng \
    nDMaterial PressureDependMultiYield02 1 2 1.8 $G $B 32 .1 80 0.5\
0 $contraction_a $contraction_b $contraction_c $contraction_d $contraction_e \
                                          26. 0.067 0.23 0.06 0.27
$dilation_a $dilation_b $dilation_c \
  }
$noYieldSurf  $liqParam1 $liqParam2 \
  2 {
$pa $S0;
    nDMaterial PressureIndependMultiYield 2 2 1.8 4.e4 2.e5 40 .1
  }
####################################################################
# define the element                    thick      maTag  bulk        fmass      hPerm   vPerm
element 9_4_QuadUP  1 1 2 3 4 5 6 7 8 9  $thickness $matTag    [expr $waterbulk/0.4]  $fluidDen  $kPerm $kPerm $gravX $gravY;
 
####################################################################
# Recorders
file mkdir output
eval "recorder Node -file output/disp1_$Analysis_case.out -time -node 1 2 3 4 5 6 7 8 9 -dof 1 disp";
eval "recorder Node -file output/disp2_$Analysis_case.out -time -node 1 2 3 4 5 6 7 8 9 -dof 2 disp";
eval "recorder Node -file output/pwp_$Analysis_case.out -time -node 1 2 3 4 -dof 3 vel";
eval "recorder Element -file output/stress_$Analysis_case.out -time -ele 1 material 9 stress";  # s11 s22 s33 s12 eta (EFFECTIVE stresses)
eval "recorder Element -file output/strain_$Analysis_case.out -time -ele 1 material 9 strain";  # e11 e22 gama12 (it is not e33 nor e12)
eval "recorder Element -file output/backbone_$Analysis_case.out -ele 1 material 9 backbone 100 1600";
####################################################################
# GRAVITY APPLICATION (elastic behavior)
model basic -ndm 2 -ndf 3
pattern Plain 1 Constant {
load 3 [expr $loadbias*$vertPress*$xsize*$thickness*0.25] [expr -$vertPress*$xsize*$thickness*0.25] 0
load 4 [expr $loadbias*$vertPress*$xsize*$thickness*0.25] [expr -$vertPress*$xsize*$thickness*0.25] 0
}
model basic -ndm 2 -ndf 2
pattern Plain 2 Constant {
load 7 [expr $loadbias*$vertPress*$xsize*$thickness*0.5] [expr -$vertPress*$xsize*$thickness*0.5]
}
}


              #   ele#                    thick maTag    bulk  mDensity  perm1  perm2  gravity     
numberer RCM
element 9_4_QuadUP  1     1 2 3 4 5 6 7 8 9  1.0   1       $bulk $fmass  $hperm  $vperm $gravX $gravY   
system ProfileSPD
test NormDispIncr 1.e-5 50 2
constraints Plain; #Penalty 1.e18 1.e18
rayleigh [expr 1.0*$Mproprayleigh] 0.0 $Kinitproprayleigh 0.0
set gama 1.5;
set betta [expr pow($gama+0.5, 2)/4]; #[expr 1./4.]; # 0.25 for const accel, 1/6 for linear accel (in this case deltaT<period/pi)
integrator Newmark $gama $betta;
algorithm KrylovNewton
analysis VariableTransient


#set material to elastic for gravity loading
updateMaterialStage -material 1 -stage 0
updateMaterialStage -material $matOpt -stage 0 
updateMaterialStage -material 2 -stage 0
updateMaterialStage -material 3 -stage 0
updateMaterialStage -material 4 -stage 0
for {set i 1} {$i <= 100} {incr i 1} {
analyze 1 1000
set tCurrent [getTime]
puts "time = $tCurrent sec"
}
for {set i 1} {$i <= 100} {incr i 1} {
analyze 1 1000
set tCurrent [getTime]
puts "time = $tCurrent sec"
}
puts "Done Elastic"
updateMaterialStage -material 1 -stage 1
updateMaterialStage -material 2 -stage 1
updateMaterialStage -material 3 -stage 1
updateMaterialStage -material 4 -stage 1
for {set i 1} {$i <= 100} {incr i 1} {
analyze 1 100
set tCurrent [getTime]
puts "time = $tCurrent sec"
}


#recorder for nodal variables along the vertical center line.
for {set i 1} {$i <= 200} {incr i 1} {
set SnodeList {}
analyze 1 1000
for {set i 0} {$i < 9} {incr i 1} {
set tCurrent [getTime]
  lappend SnodeList [expr $i+1]
puts "time = $tCurrent sec"
}
}
puts "Done Plastic"
puts "Gravity Done!"


set FnodeList {}
####################################################################
for {set i 0} {$i < 4} {incr i 1} {
# Adjust some fixities for shear loading
  lappend FnodeList [expr $i+1]
 
# remove surface pwp fixity
model basic -ndm 2 -ndf 3
remove sp 3 3
remove sp 4 3
equalDOF 3 4  3
 
####################################################################
loadConst -time 0.0
wipeAnalysis
 
####################################################################
# NOW APPLY LOADING SEQUENCE AND ANALYZE (plastic)
 
# Load control cyclic
if {$Analysis_case == "undrained_cyclic"} {
 
# Define analysis parameters
numberer RCM
system ProfileSPD
test NormDispIncr 1.e-5 50 0
constraints Plain
rayleigh $Mproprayleigh  0.0 $Kinitproprayleigh 0.0
set gama 1.5; #0.5; #gama=alfa in other texts
set betta [expr pow($gama+0.5, 2)/4]; # 0.25 for const accel, 1/6 for linear accel (in this case deltaT<period/pi)
integrator Newmark $gama $betta;
algorithm KrylovNewton
analysis VariableTransient
model basic -ndm 2 -ndf 3
pattern Plain 4 "Sine 0 [expr $period*1000] $period" {
load 3 [expr $csr*$vertPress*$xsize*$thickness*0.25] 0 0
load 4 [expr $csr*$vertPress*$xsize*$thickness*0.25] 0 0
}
model basic -ndm 2 -ndf 2
pattern Plain 5 "Sine 0 [expr $period*1000] $period" {
load 7 [expr $csr*$vertPress*$xsize*$thickness*0.5] 0
}
puts "$Analysis_case"
while {[expr abs([nodeDisp 7 1]/$ysize)]<$target_shear_strain} {
# while {[getTime]<0.25}
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
puts "time = [getTime] sec"
}
}
}


#############################################################
# Disp control monotonic
# GRAVITY APPLICATION (elastic behavior)
if {$Analysis_case == "undrained_monotonic"} {
# Define analysis parameters
numberer RCM
system BandGeneral; #ProfileSPD
test NormDispIncr 1.e-5 50 0; #1.e-5 50 0
constraints Penalty 1.e10 1.e10; #for clay use 1.e5, for sand use 1.e15
rayleigh $Mproprayleigh  0.0 $Kinitproprayleigh 0.0
set gama 1.5; #0.5; #gama=alfa in other texts
set betta [expr pow($gama+0.5, 2)/4]; # 0.25 for const accel, 1/6 for linear accel (in this case deltaT<period/pi)
integrator Newmark $gama $betta;
algorithm KrylovNewton
analysis VariableTransient
model basic -ndm 2 -ndf 2
pattern Plain 14 Linear {
sp 3 1 [expr $ysize*$target_shear_strain]
sp 4 1 [expr $ysize*$target_shear_strain]
sp 7 1 [expr $ysize*$target_shear_strain]
}
puts "$Analysis_case"
for {set i 1} {$i <= [expr $numSteps]} {incr i 1} {
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
set tCurrent [getTime]
puts "time = $tCurrent sec"
}
}


# create the SOE, ConstraintHandler, Integrator, Algorithm and Numberer
# Disp cotrol monotonic
numberer RCM
if {$Analysis_case == "drained_monotonic"} {
system ProfileSPD
test NormDispIncr 1.0e-8 30 0
# Define analysis parameters
algorithm KrylovNewton
numberer RCM
constraints Penalty 1.e18 1.e18
system BandGeneral; #ProfileSPD
set nw 1.5
test NormDispIncr 1.e-5 50 0
set nw2 [expr pow($nw+0.5, 2)/4]
constraints Penalty 1.e18 1.e18
integrator Newmark $nw $nw2
rayleigh $Mproprayleigh  0.0 $Kinitproprayleigh 0.0
analysis Transient
set gama 1.5; #0.5; #gama=alfa in other texts
set betta [expr pow($gama+0.5, 2)/4]; # 0.25 for const accel, 1/6 for linear accel (in this case deltaT<period/pi)
integrator Newmark $gama $betta;
algorithm KrylovNewton
analysis VariableTransient
model basic -ndm 2 -ndf 3
pattern Plain 15 Linear {
sp 3 1 [expr $ysize*$target_shear_strain]
sp 4 1 [expr $ysize*$target_shear_strain]
sp 7 1 [expr $ysize*$target_shear_strain]
}
puts "$Analysis_case"
for {set i 1} {$i <= [expr $numSteps]} {incr i 1} {
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
set tCurrent [getTime]
puts "time = $tCurrent sec"
}
}


analyze 10 5e3
# Disp control cyclic for G/Gmax and Xi curves
if {$Analysis_case == "drained_cyclic"} {
# Define analysis parameters
numberer RCM
system BandGeneral; #ProfileSPD
test NormDispIncr 1.e-5 50 0
constraints Penalty 1.e18 1.e18
rayleigh $Mproprayleigh  0.0 $Kinitproprayleigh 0.0
set gama 1.5; #0.5; #gama=alfa in other texts
set betta [expr pow($gama+0.5, 2)/4]; # 0.25 for const accel, 1/6 for linear accel (in this case deltaT<period/pi)
integrator Newmark $gama $betta;
algorithm KrylovNewton
analysis VariableTransient
model basic -ndm 2 -ndf 2
pattern Plain 4 "Sine 0 1000 $period" {
sp 3 1 [expr $ysize*0.0003/100.]
}
puts "$Analysis_case"
for {set i 1} {$i <= $numSteps} {incr i 1} {
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
set tCurrent [getTime]
puts "time = $tCurrent sec"
}
remove loadPattern 4
loadConst;


updateMaterialStage -material $matOpt -stage 1
pattern Plain 5 "Sine 0 1000 $period" {
sp 3 1 [expr $ysize*0.001/100.]
}
for {set i 1} {$i <= $numSteps} {incr i 1} {
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
set tCurrent [getTime]
puts "time = $tCurrent sec"
}
remove loadPattern 5
loadConst


analyze 100 1.e0
pattern Plain 6 "Sine 0 1000 $period" {
sp 3 1 [expr $ysize*0.003*0.01]
}
for {set i 1} {$i <= $numSteps} {incr i 1} {
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
set tCurrent [getTime]
puts "time = $tCurrent sec"
}
remove loadPattern 6
loadConst


# rezero time
pattern Plain 7 "Sine 0 1000 $period" {
wipeAnalysis
sp 3 1 [expr $ysize*0.01*0.01]
setTime 0.0
}
#############################################################
for {set i 1} {$i <= $numSteps} {incr i 1} {
# NOW APPLY LOADING SEQUENCE AND ANALYZE (plastic)
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
set tCurrent [getTime]
puts "time = $tCurrent sec"
}
remove loadPattern 7
loadConst


#  base input motion                   
pattern Plain 8 "Sine 0 1000 $period" {
pattern UniformExcitation  1  1    -accel "Sine 0. 10. $period -factor $accMul"
sp 3 1 [expr $ysize*0.03*0.01]
}
for {set i 1} {$i <= $numSteps} {incr i 1} {
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
set tCurrent [getTime]
puts "time = $tCurrent sec"
}
remove loadPattern 8
loadConst


eval "recorder Node -file disp  -time -node $SnodeList -dof 1 2 -dT $deltaT disp"
pattern Plain 9 "Sine 0 1000 $period" {
eval "recorder Node -file  pwp -time -node $FnodeList -dof 3 -dT $deltaT vel"
sp 3 1 [expr $ysize*0.1*0.01]
eval "recorder Node -file acc  -time -node $SnodeList -dof 1 2 -dT $deltaT accel"
}
recorder Element -ele 1  -time -file stress1  -dT $deltaT material 1 stress
for {set i 1} {$i <= $numSteps} {incr i 1} {
recorder Element -ele 1 -time -file strain1  -dT $deltaT material 1 strain
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
recorder Element -ele 1  -time -file stress5  -dT $deltaT material 5 stress
set tCurrent [getTime]
recorder Element -ele 1  -time -file strain5  -dT $deltaT material 5 strain
puts "time = $tCurrent sec"
recorder Element -ele 1  -time -file stress9  -dT $deltaT material 9 stress
}
recorder Element -ele 1  -time -file strain9  -dT $deltaT material 9 strain
remove loadPattern 9
loadConst


pattern Plain 10 "Sine 0 1000 $period" {
sp 3 1 [expr $ysize*0.3*0.01]
}
for {set i 1} {$i <= $numSteps} {incr i 1} {
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
set tCurrent [getTime]
puts "time = $tCurrent sec"
}
remove loadPattern 10
loadConst


constraints Penalty 1.e18 1.e18
pattern Plain 11 "Sine 0 1000 $period" {
test NormDispIncr 1.e-4 25 0
sp 3 1 [expr $ysize*1.0*0.01]
numberer RCM
}
algorithm KrylovNewton
for {set i 1} {$i <= $numSteps} {incr i 1} {
system ProfileSPD
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
integrator Newmark $gamma  [expr pow($gamma+0.5, 2)/4]
set tCurrent [getTime]
rayleigh $massProportionalDamping 0.0 $InitStiffnessProportionalDamping 0.0
puts "time = $tCurrent sec"
analysis VariableTransient
}
remove loadPattern 11
loadConst


set startT [clock seconds]
pattern Plain 12 "Sine 0 1000 $period" {
analyze $numSteps $deltaT [expr $deltaT/100] $deltaT 15
sp 3 1 [expr $ysize*3.0*0.01]
}
for {set i 1} {$i <= [expr $numSteps*1.25]} {incr i 1} {
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
set tCurrent [getTime]
puts "time = $tCurrent sec"
}
}
####################################################################
set endT [clock seconds]
set endT [clock seconds]
puts "Execution time: [expr $endT-$startT] seconds."
puts "Execution time: [expr $endT-$startT] seconds."
puts "Finished with cyclic loading"
wipe all;
wipe;  # flush ouput stream
</syntaxhighlight>
== DSS_material_properties.tcl file ==
<syntaxhighlight lang="tcl">
#############################################################################
# PDMY03 input parameters materials     #
# By: Arash Khosravifar  January 2009 #
# Last edit:   August 2018 #
#############################################################################
if {$matTag == 1} {
# defined variables for Sand (N1)60=35 "nonliq" (Mat 1)
set massDen 2.06;  # (ton/m3)
set refG 111.9e3;# (kPa)
set refB 298.3e3;# (kPa)
set frinctionAng 42.2; # (degree)
set peakShearStrain 0.1;
set refPress 100.;  # (kPa)
set pressDependCoe 0.5;
set phaseTransAng 37.2; # (degree)
set contraction_a 0.001;  # Contraction rate.
set contraction_b 0.0;  # fabric damage
set contraction_c 0.8;  # k_sigma effect
set contraction_d 2.2; # ~CRR*[3/(1+2k0)]/2
set contraction_e 0.0;
set dilation_a 0.6;
set dilation_b 3.0;
set dilation_c -0.5; # k_sigma effect
set liqParam1 1.0;
set liqParam2 0.0;
set noYieldSurf 20;
# set void 0.55;
# set cs1 0.9;
# set cs2 0.02;
# set cs3 0.0;
set pa 100;  # (kPa)
set S0 1.73;  # (kPa)
}
if {$matTag == 2} {
# defined variables for Sand (N1)60=25 (Mat 2)
set massDen 2.03;  #(ton/m3)
set refG 94.6e3; # (kPa)
set refB 252.6e3;# (kPa)
set frinctionAng 35.8;
set peakShearStrain 0.1;
set refPress 100; # (kPa)
set pressDependCoe 0.5;
set phaseTransAng 30.8;
set contraction_a 0.005;
set contraction_b 1.0;
set contraction_c 0.6; 
set contraction_d 4.6; # ~CRR*[3/(1+2k0)]/2
set contraction_e -1.0;
set dilation_a 0.45;
set dilation_b 3.0;
set dilation_c -0.4;
set liqParam1 1.0;
set liqParam2 0.0;
set noYieldSurf 20;
# set void 0.60;
# set cs1 0.9;
# set cs2 0.02;
# set cs3 0.0;
set pa 100;  # (kPa)
set S0 1.73; # (kPa)
}
if {$matTag == 3} {
# defined variables for Sand (N1)60=15 (Mat 3)
set massDen 1.99; #(ton/m3)
set refG 73.7e3; # (kPa)
set refB 196.8e3;# (kPa)
set frinctionAng 30.3;
set peakShearStrain 0.1;
set refPress 100; # (kPa)
set pressDependCoe 0.5;
set phaseTransAng 25.3;
set contraction_a 0.012;
set contraction_b 3.0;
set contraction_c 0.4;
set contraction_d 9.0; # ~CRR*[3/(1+2k0)]/2
set contraction_e 0.0;
set dilation_a 0.3;
set dilation_b 3.0;
set dilation_c -0.3;
set liqParam1 1.0;
set liqParam2 0.0;
set noYieldSurf 20;
# set void 0.67;
# set cs1 0.9;
# set cs2 0.02;
# set cs3 0.0;
set pa 100;  # (kPa)
set S0 1.73; # (kPa)
}


wipe  #flush ouput stream
if {$matTag == 4} {
</syntaxhighlight>
# defined variables for Sand (N1)60=5 (Mat 4)
set massDen 1.94; #(ton/m3)
set refG 46.9e3; # (kPa)
set refB 125.1e3;# (kPa)
set frinctionAng 25.4;
set peakShearStrain 0.1;
set refPress 100; # (kPa)
set pressDependCoe 0.5;
set phaseTransAng 20.0;


set contraction_a 0.03;
set contraction_b 5.0;
set contraction_c 0.2;
set contraction_d 16.0; # ~CRR*[3/(1+2k0)]/2
set contraction_e 2.0;
set dilation_a 0.15;
set dilation_b 3.0;
set dilation_c -0.2;


set liqParam1 1.0;
set liqParam2 0.0;
set noYieldSurf 20;
# set void 0.76;
# set cs1 0.9;
# set cs2 0.02;
# set cs3 0.0;
set pa 100;  # (kPa)
set S0 1.73; # (kPa)
}
</syntaxhighlight>


== MATLAB Plotting File ==
== MATLAB Plotting File ==
<syntaxhighlight lang="matlab">
<syntaxhighlight lang="matlab">
clear all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Post processor for 9-4quadUP element                  %
% By: Arash Khosravifar    January 2009                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
clc
clear all
close all
 
% Type of analysis. Options:
% undrained_cyclic
% undrained_monotonic
% drained_cyclic
% drained_monotonic
Analysis_case = 'undrained_cyclic'; 
 
% Loading files
input_disp1      = load(['output\disp1_',Analysis_case,'.out']);
input_disp2      = load(['output\disp2_',Analysis_case,'.out']);
input_stress    = load(['output\stress_',Analysis_case,'.out']);
input_strain    = load(['output\strain_',Analysis_case,'.out']);
input_pwp      = load(['output\pwp_',Analysis_case,'.out']);
input_backbone  = load(['output\backbone_',Analysis_case,'.out']);
 
 
% Displacements
Disp_1 = input_disp1(:,2:end);
Disp_2 = input_disp2(:,2:end);
% Stresses  compression:+  tension:-
s_11  = -input_stress(:,2);
s_22  = -input_stress(:,3);
s_33  = -input_stress(:,4);
s_12  = input_stress(:,5);
eta  = input_stress(:,6); % defined in the manual - not used


a1=load('acc');
tau_oct  = 1/3*sqrt((s_11-s_22).^2 + (s_11-s_33).^2 + (s_22-s_33).^2 + 6*s_12.^2); % Octahedral shear stress
d1=load('disp');
q    = 3/sqrt(2)*tau_oct; %1/3*sign(s_12).*sqrt((s_11-s_22).^2 + (s_11-s_33).^2 + (s_22-s_33).^2 + 6*s_12.^2);
p1=load('pwp');
pwp  = (input_pwp(:,2)+input_pwp(:,3)+input_pwp(:,4)+input_pwp(:,5))/4;
s1=load('stress1');
p    = (s_11 + s_22 + s_33)/3;
e1=load('strain1');
p_tot = p + pwp;
s5=load('stress5');
time  = input_pwp(:,1);
e5=load('strain5');
s9=load('stress9');
e9=load('strain9');


fs=[0.5, 0.2, 4, 6];
%Note: The stresses are already effective. no need to subtract
fs2=[0.5, 0.2, 4, 3];
% Strains  Upward (dilation):-  Downward (contraction):+
accMul = 2;
e_11  = -input_strain(:,2);
e_22  = -input_strain(:,3);
gama_12 = input_strain(:,4);
e_33  = zeros(size(e_11));      % plain strain
e_v  = e_11 + e_22 + e_33;     % e_33 is zero
gama_oct = 2/3*sqrt((e_11-e_22).^2+(e_22-e_33).^2+(e_11-e_33).^2+6*(gama_12/2).^2);  


%integration point 1 p-q
% CSL line
po=(s1(:,2)+s1(:,3)+s1(:,4))/3;
cs1 = 0.9;
for i=1:size(s1,1)
cs2 = 0.1;
    qo(i)=(s1(i,2)-s1(i,3))^2 + (s1(i,3)-s1(i,4))^2 +(s1(i,2)-s1(i,4))^2 + 6.0* s1(i,5)^2;
cs3 = 1.0;
  qo(i)=sign(s1(i,5))*1/3.0*qo(i)^0.5;
void_i = 0.8;    % initial void ratio
end
void = void_i - e_22*(1+void_i);


figure(1); close 1; figure(1);
% Number of steps in "consolidation phase"
%integration point 1 stress-strain
consol_step = 500;
subplot(2,1,1), plot(e1(:,4),s1(:,5),'r');
total_step  = size(gama_12,1);
title ('shear stress \tau_x_y VS. shear strain \epsilon_x_y at integration point 1');
xLabel('Shear strain \epsilon_x_y');
yLabel('Shear stress \tau_x_y (kPa)');
subplot(2,1,2), plot(-po,qo,'r');
title ('confinement p VS. deviatoric stress q at integration point 1');
xLabel('confinement p (kPa)');
yLabel('q (kPa)');
set(gcf,'paperposition',fs);
saveas(gcf,'SS_PQ_p1','jpg');


% Ru
ru_1 = pwp./max(p,0.001);
ru_2 = pwp/p(consol_step);
ru_3 = pwp/s_22(consol_step);


%integration point 5 p-q
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
po=(s5(:,2)+s5(:,3)+s5(:,4))/3;
if strcmp(Analysis_case,'drained_cyclic') == 1
for i=1:size(s5,1)
     period      = 10;
     qo(i)=(s5(i,2)-s5(i,3))^2 + (s5(i,3)-s5(i,4))^2 +(s5(i,2)-s5(i,4))^2 + 6.0* s5(i,5)^2;
    cycle      = time/period;
  qo(i)=sign(s5(i,5))*1/3.0*qo(i)^0.5;
end


figure(5); close 5; figure(5);
    % G/Gmax
%integration point 5 stress-strain
    NoStepsCycle  = 2000;
subplot(2,1,1), plot(e5(:,4),s5(:,5),'r');
    NoCycleEach  = 2;
title ('shear stress \tau_x_y VS. shear strain \epsilon_x_y at integration point 5');
    corr_steps = [consol_step+NoStepsCycle/4:NoStepsCycle*NoCycleEach : consol_step+9*NoStepsCycle*NoCycleEach];
xLabel('Shear strain \epsilon_x_y');
    G_sec = s_12(corr_steps)./gama_12(corr_steps);     %Secant stiffness
yLabel('Shear stress \tau_x_y (kPa)');
    G_max = s_12(consol_step+1)/gama_12(consol_step+1);             %Secant stiffness at the first load incr.
subplot(2,1,2), plot(-po,qo,'r');
    G_max_oct = tau_oct(consol_step+1)/gama_oct(consol_step+1);    
title ('confinement p VS. deviatoric stress q at integration point 5');
    % Note that the Gmax,r assigned for the material by the user is not the
xLabel('confinement p (kPa)');
    % conventional Gmax (s_12/gama_12), but it is the OCTAHEDRAL one
yLabel('q (kPa)');
    % (tau_oct/gama_oct)
set(gcf,'paperposition',fs);
saveas(gcf,'SS_PQ_p5','jpg');


%integration point 9 p-q
    % EPRI 1993
po=(s9(:,2)+s9(:,3)+s9(:,4))/3;
    GratioEPRI6 =  [1,0.98,0.919,0.751,0.519,0.271,0.121,0.045,0.024];  % EPRI1993 for 0-6m corresponding to strain=[0.000003,0.00001,0.00003,0.0001,0.0003,0.001,0.003,0.01,0.03]
for i=1:size(s1,1)
    GratioEPRI76 = [1,0.999,0.982,0.9,0.752,0.502,0.282,0.124,0.07];   % EPRI1993 for 36-76m
     qo(i)=(s9(i,2)-s9(i,3))^2 + (s9(i,3)-s9(i,4))^2 +(s9(i,2)-s9(i,4))^2 + 6.0* s9(i,5)^2;
    gama_EPRI  = [0.000003,0.00001,0.00003,0.0001,0.0003,0.001,0.003,0.01,0.03];
  qo(i)=sign(s9(i,5))*1/3.0*qo(i)^0.5;
     gama_EPRI_oct = 2/3*sqrt((e_11(500)-e_22(500)).^2+(e_22(500)-e_33(500)).^2+(e_11(500)-e_33(500)).^2+6*(gama_EPRI/2).^2);
end


figure(6); close 6; figure(6);
    % Use backbone record
%integration point 9 stress-strain
    backbone = input_backbone(1,:)';
subplot(2,1,1), plot(e9(:,4),s9(:,5),'r');
    gama_backbone_100_oct = input_backbone(5:4:end);
title ('shear stress \tau_x_y VS. shear strain \epsilon_x_y at integration point 9');
    G_backbone_100 = input_backbone(6:4:end);
xLabel('Shear strain \epsilon_x_y');
    gama_backbone_1600_oct = input_backbone(7:4:end);
yLabel('Shear stress \tau_x_y (kPa)');
    G_backbone_1600 = input_backbone(8:4:end);
subplot(2,1,2), plot(-po,qo,'r');
title ('confinement p VS. deviatoric stress q at integration point 9');
xLabel('confinement p (kPa)');
yLabel('q (kPa)');
set(gcf,'paperposition',fs);
saveas(gcf,'SS_PQ_p9','jpg');


figure(2); close 2; figure(2);
    % Xi
%node 3 displacement relative to node 1
    for ii = 1:length(corr_steps)
plot(d1(:,1),d1(:,6));
        hysteresis = trapz(gama_12(corr_steps(ii):corr_steps(ii)+NoStepsCycle),s_12(corr_steps(ii):corr_steps(ii)+NoStepsCycle));
title ('Lateral displacement at element top');
        triangle  = 1/2*gama_12(corr_steps(ii))*s_12(corr_steps(ii));
xLabel('Time (s)');
        Xi(ii)    = hysteresis/(4*pi*triangle)*100;   %(%)
yLabel('Displacement (m)');
    end
set(gcf,'paperposition',fs2);
    XiEPRI6 =  [0.791,1.27,2.28,4.26,7.75,14.74,22.03,27.12,30.21];    % (%) EPRI1993 for 0-6m corresponding to strain=[0.000003,0.00001,0.00003,0.0001,0.0003,0.001,0.003,0.01,0.03]
saveas(gcf,'Disp','jpg');
    XiEPRI76 = [0.591,0.88,1.271,2.46,4.551,8.34,14.13,21.22,25.31];    % (%) EPRI1993 for 36-76m
   
% Plots
    figure
    semilogx(gama_12(corr_steps),G_sec/G_max,'-ob','LineWidth',2);
    hold on; grid on;
    semilogx(gama_EPRI,GratioEPRI6,'--r','LineWidth',2);
    semilogx(gama_EPRI,GratioEPRI76,'--r','LineWidth',2);
    xlabel('\gamma_{12}');
    ylabel('G_{sec}/G_{max}');
    legend('PDMY03 - \sigma_{vo}=100 KPa','EPRI1993 0-6m','EPRI1993 36-76m')


s=accMul*sin(0:pi/50:20*pi);
    % Xi
s=[s';zeros(3000,1)];
    figure
s1=interp1(0:0.01:40,s,a1(:,1));
    semilogx(gama_12(corr_steps),Xi,'-ob','LineWidth',2);
    hold on; grid on;
    semilogx(gama_12(corr_steps),XiEPRI6,'--r','LineWidth',2);
    semilogx(gama_12(corr_steps),XiEPRI76,'--r','LineWidth',2);
    xlabel('\gamma_{12}');
    ylabel('Equivalent damping ratio (%)');
    legend('PDMY03 - \sigma_{vo}=100 KPa','EPRI1993 0-6m','EPRI1993 36-76m')
   
    figure
    subplot(2,1,1);
    plot(gama_12,s_12/s_22(consol_step),'blue');
    xlabel('\gamma_{12}');
    ylabel('CSR = (\sigma_{12}/\sigma''_{vc})');
    title(['\sigma''_{vc} = ',num2str(s_22(consol_step)),' kPa']);


figure(3); close 3; figure(3);
    subplot(2,1,2);
%node acceleration
    plot(cycle(consol_step+1:end),gama_12(consol_step+1:end),'blue');
a = plot(a1(:,1),s1+a1(:,6),'r');
    xlabel('No. of cycles');
title ('Lateral acceleration at element top');
    ylabel('\gamma_{12}');
xLabel('Time (s)');
       
yLabel('Acceleration (m/s^2)');
 
set(gcf,'paperposition',fs2);
end;   % (end of drained_cyclic)
saveas(gcf,'Acc','jpg');
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


figure(4); close 4; figure(4);
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=plot(p1(:,1),p1(:,2));
if strcmp(Analysis_case,'undrained_cyclic') == 1
title ('Pore pressure at base');
    period      = 1;
xLabel('Time (s)');
    cycle      = time/period;
yLabel('Pore pressure (kPa)');
set(gcf,'paperposition',fs2);
saveas(gcf,'EPWP','jpg');


</syntaxhighlight>
    % find "Liquefaction Triggering" time
    gama_consol = gama_12(consol_step);
 
    T_1  = find(abs(gama_12-gama_consol)>0.01,1);
    if any(T_1) == 0
        T_1 = consol_step+1;
    end
    N_1 = cycle(T_1)
 
    T_3 = find(abs(gama_12-gama_consol)>0.03,1);
    if any(T_3) == 0
        T_3 = consol_step+1;
    end
    N_3 = cycle(T_3)
 
    % T_r = find(ru_3(consol_step:total_step)>0.98,1);
    T_r = find(ru_3>0.98,1);
    if any(T_r) == 0
        T_r = consol_step+1;
    end
    N_r = cycle(T_r)
 
    % Plotting
    figure
    subplot(2,2,1);
    plot(gama_12,s_12/s_22(consol_step));grid on; hold on;xlabel('\gamma_{12}');ylabel('CSR = (\sigma_{12}/\sigma''_{vc})')
    plot(gama_12(T_1),s_12(T_1)/s_22(consol_step),'g.','MarkerSize',15);
    plot(gama_12(T_3),s_12(T_3)/s_22(consol_step),'r.','MarkerSize',15);
    plot(gama_12(T_r),s_12(T_r)/s_22(consol_step),'m.','MarkerSize',15);
    title(['CSR = ',num2str(max(abs(s_12))/s_22(consol_step),2),'  and  \sigma''_{vc} = ',num2str(s_22(consol_step)),' kPa']);
%    xlim([-0.1,0.1]);
 
    subplot(2,2,2);plot(s_22/s_22(consol_step),s_12/s_22(consol_step));grid on; hold on;xlabel('vertical effective stress ratio (\sigma''_{v}/\sigma''_{vc})');ylabel('CSR = (\sigma_{12}/\sigma''_{vc})');
    plot(s_22(T_1)/s_22(consol_step),s_12(T_1)/s_22(consol_step),'g.','MarkerSize',15);
    plot(s_22(T_3)/s_22(consol_step),s_12(T_3)/s_22(consol_step),'r.','MarkerSize',15);
    plot(s_22(T_r)/s_22(consol_step),s_12(T_r)/s_22(consol_step),'m.','MarkerSize',15);
    xlim([min(0,min(s_22/s_22(consol_step))),max(s_22)/s_22(consol_step)]);
    title(['CSR = ',num2str(max(abs(s_12))/s_22(consol_step),2),'  and  \sigma''_{vc} = ',num2str(s_22(consol_step)),' kPa']);
 
    subplot(2,2,3);plot(cycle(consol_step+1:end),pwp(consol_step+1:end));grid on;ylabel('PWP (kPa)');xlabel('No. of cycles');
    hold on;
    plot(cycle(T_1),pwp(T_1),'g.','MarkerSize',15);
    plot(cycle(T_3),pwp(T_3),'r.','MarkerSize',15);
    plot(cycle(T_r),pwp(T_r),'m.','MarkerSize',15);
    ylim([-10,100]);
 
    subplot(2,2,4);plot(cycle(consol_step+1:end),gama_12(consol_step+1:end));grid on;xlabel('No. of cycles');ylabel('\gamma_{12}');
 
end;    % (end of undrained_cyclic)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(Analysis_case,'drained_monotonic') == 1
 
    G_max = s_12(consol_step+1)/gama_12(consol_step+1);
    G_max_oct = tau_oct(consol_step+1)/gama_oct(consol_step+1);     
   
   
    figure
    subplot(1,2,1);hold on; grid on;
    plot(gama_12(consol_step+1:end),s_12(consol_step+1:end));
    plot(gama_oct(consol_step+1:end),tau_oct(consol_step+1:end),'r');
    legend('\sigma_{12} vs \gamma_{12}','\tau_{oct} vs \gamma_{oct}',4);
    xlabel('\gamma_{12} and \gamma_{oct}');ylabel('\sigma_{12} and \tau_{oct}');
    ylim([0,60]);
    title(['\sigma''_{vc} (kPa) = ',num2str(s_22(consol_step))]);


    subplot(1,2,2);grid on;hold on
    plot(p(consol_step+1:end),tau_oct(consol_step+1:end),'r');
    plot(s_22(consol_step+1:end),s_12(consol_step+1:end),'b')
    ylim([0,60]);
 
    slope_oct = tau_oct(end)/p(end);
    Phi_oct_back_calculated = asin(3*slope_oct/(2*sqrt(2)+slope_oct))/3.1416*180;
    slope_conventional = s_12(end)/s_22(end);
    Phi_conventional_back_calculated = atan(slope_conventional)/3.1416*180;


== Displacement Output File ==
    plot([0,max(p)*1.2],[0,(max(p)*1.2)*slope_oct],'--r')
[[File:PD02_Ex23Disp.jpg]]
    plot([0,max(s_22)*1.2],[0,(max(s_22)*1.2)*slope_conventional],'--b')
    xlabel('\sigma''_{v} and p'' (kPa)');ylabel('\sigma_{12} and \tau_{oct}');
    xlim([0,max(p)*1.2]);
    legend(['\phi_{oct} = ',num2str(Phi_oct_back_calculated)],['\phi_{DSS} = ',num2str(Phi_conventional_back_calculated)],2)


== Stress-Strain Output File (integration point 1) ==
end;    % (end of drained_monotonic)
[[File:PD02_Ex23SS_PQ1.jpg]]
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


== Stress-Strain Output File (integration point 5) ==
[[File:PD02_Ex23SS_PQ5.jpg]]


== Stress-Strain Output File (integration point 9) ==
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[[File:PD02_Ex23SS_PQ9.jpg]]
if strcmp(Analysis_case,'undrained_monotonic') == 1
    figure
    subplot(2,2,1);plot(gama_12(consol_step+1:end),s_12(consol_step+1:end));
    grid on;
    xlabel('\gamma_{12}');
    ylabel('\sigma_{12} (kPa)');
   
    subplot(2,2,2);plot(s_22(consol_step+1:end),s_12(consol_step+1:end));
    grid on; hold on;
    xlim([0,max(s_22)])
    plot([0,s_22(end)],[0,s_12(end)],'r')
    xlabel('\sigma''_{v} (kPa)');
    ylabel('\sigma_{12} (kPa)');
    slope_conventional = s_12(end)/s_22(end);
    Phi_conventional_back_calculated = atan(slope_conventional)/3.1416*180;
    legend(['\phi_{DSS} = ',num2str(Phi_conventional_back_calculated)],2)


== Excess Pore Pressure Output File ==
    subplot(2,2,3);plot(gama_12,pwp);grid on;
[[File:PD02_Ex23EPP.jpg]]
    xlabel('\gamma_{12}')
    ylabel('pwp (kPa)')
    subplot(2,2,4)
    plot(gama_12,e_v-e_v(consol_step)); grid on;
    set(gca,'YDir','rev');
    xlabel('\gamma_{12}')
    ylabel('Vomuletric strain \epsilon_v - \epsilon_o');


== Acceleration Output File ==
end;    % (end of undrained_monotonic')
[[File:PD02_Ex23Accel.jpg]]
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
</syntaxhighlight>




Latest revision as of 19:50, 20 December 2019



Element Drivers for PressureDependMultiYield03 Material


DSS_main.tcl file

###########################################################################
# Model with 9-4-quad-UP element and PDMY03 material				      #
# By: Arash Khosravifar   January 2009									  #
# Last edit:			  August 2018									  #
###########################################################################

wipe;
wipe all;
set startT [clock seconds];

######################### MAIN USER INPUTS #################################
set matTag				2;				 # 1 for (N1)60=35, 2 for (N1)60=25, 3 for (N1)60=15 and 4 for (N1)60=5
set vertPress 			[expr 1.0*100];  # kPa vertical confining pressure
set loadbias			0.3; 			 # Alpha = tau_xy/s'vo
set Analysis_case "undrained_cyclic"
# Options:
# undrained_cyclic
# undrained_monotonic: define target max shear strain
# drained_cyclic: applies 2 cycles for each strain (9 shear strain ranging from 0.0003% to 3%) for the purpose of G/Gmax and Damping curves
# drained_monotonic: define target max shear strain
######################### END OF MAIN USER INPUTS ###########################

# Load material properties
source DSS_material_properties.tcl

# Some variables for the ELEMENT
set fluidDen 			1.0;		# Fluid mass density
set waterbulk 			2.2e6;		# kPa
set kdrain 				1.e1;		# permeability for drained loading
set kundrain 			1.e-8;		# permeability for undrained loading
set gravY 				0.0;
set gravX 				0.0;

# Some loading related variables
if {$Analysis_case == "undrained_cyclic"} {
	set target_shear_strain [expr 3.0/100.]
	set csr 				0.45;
	set period      		1.;
	set deltaT 				0.005;
	set kPerm				$kundrain
}
if {$Analysis_case == "undrained_monotonic"} {
	set target_shear_strain [expr 3.0/100.]
	set period      		1.;
	set deltaT 				0.005;
	set numSteps 			200;
	set kPerm				$kundrain
	# also change penalty from 1e15 and 1e5 for sand to clay
}
if {$Analysis_case == "drained_monotonic"} {
	# also change the damping to 0.005
	set target_shear_strain [expr 1.0/100.];	#unit in strain - it will be multiplied by ySize later. Note: DON'T forget the ".0" in 1.0
	set period      		1.;
	set deltaT 				0.005;
	set numSteps 			3000;

	set kPerm				$kdrain
	set phaseTransAng $frinctionAng
	set contraction_a 	0.0;
	set contraction_b 	0.0;
	set contraction_c 	0.0;
	set contraction_d 	0.0;
	set contraction_e 	0.0;
	set dilation_a 		0.0;
	set dilation_b 		0.0;
	set dilation_c 		0.0;
	set liqParam1 			0.0;
	set liqParam2 			0.0;
}
if {$Analysis_case == "drained_cyclic"} {
	set period      		10.;
	set deltaT 				0.005;
	set numSteps 			4000;
	
	set kPerm				$kdrain
	set phaseTransAng $frinctionAng
	set contraction_a 	0.0;
	set contraction_b 	0.0;
	set contraction_c 	0.0;
	set contraction_d 	0.0;
	set contraction_e 	0.0;
	set dilation_a 		0.0;
	set dilation_b 		0.0;
	set dilation_c 		0.0;
	set liqParam1 			0.0;
	set liqParam2 			0.0;
}

# Define Rayleigh Damping
set rayleigh_damping1	0.001;
set rayleigh_damping2	0.001;
set rayleigh_w1			[expr 2*3.14/$period];
set rayleigh_w2			[expr 2*3.14/($period/2.0)];
set Kinitproprayleigh	[expr 2*$rayleigh_w1*$rayleigh_w2/($rayleigh_w1*$rayleigh_w1-$rayleigh_w2*$rayleigh_w2)*(-1*$rayleigh_damping1/$rayleigh_w1+1*$rayleigh_damping2/$rayleigh_w2)]; #This is a1 [expr 2*$dampratio/(2*3.1416/0.02)];
set Mproprayleigh		[expr 2*$rayleigh_w1*$rayleigh_w2/($rayleigh_w1*$rayleigh_w1-$rayleigh_w2*$rayleigh_w2)*($rayleigh_damping1*$rayleigh_w1-$rayleigh_damping2*$rayleigh_w2)]; #This is a0 [expr 2*$dampratio*(2*3.1416/2.0)];

####################################################################
set xsize 0.1;
set ysize 0.1;
set thickness 0.1;

# define the 3DOF nodes
model basic -ndm 2 -ndf 3
node 1   0.0 0.0 
node 2   $xsize 0.0 
node 3   $xsize $ysize 
node 4   0.0 $ysize

fix 1  1 1 0
fix 2  1 1 0
fix 3  0 0 1
fix 4  0 0 1
equalDOF 1 2  3
equalDOF 3 4  1 2

# define the 2DOF nodes
model basic -ndm 2 -ndf 2
node 5   [expr $xsize/2] 0.0
node 6   $xsize [expr $ysize/2]
node 7   [expr $xsize/2] $ysize
node 8   0.0 [expr $ysize/2]
node 9   [expr $xsize/2] [expr $ysize/2]

fix 5 1 1
equalDOF 3 7  1 2
equalDOF 6 8  1 2
# equalDOF 6 9  1 2


################# define material ##################################
 nDMaterial PressureDependMultiYield03 $matTag 2 $massDen $refG $refB $frinctionAng \
	 $peakShearStrain $refPress $pressDependCoe $phaseTransAng \
	 0 $contraction_a $contraction_b $contraction_c $contraction_d $contraction_e \
	 $dilation_a $dilation_b $dilation_c \
	 $noYieldSurf  $liqParam1 $liqParam2 \
	 $pa $S0;	
	
	
####################################################################
# define the element                     thick       maTag  		bulk         fmass      hPerm   vPerm
element 9_4_QuadUP  1  1 2 3 4 5 6 7 8 9  $thickness $matTag    [expr $waterbulk/0.4]  $fluidDen  $kPerm $kPerm $gravX $gravY;

####################################################################
# Recorders
file mkdir output
eval "recorder Node -file output/disp1_$Analysis_case.out -time -node 1 2 3 4 5 6 7 8 9 -dof 1 disp";
eval "recorder Node -file output/disp2_$Analysis_case.out -time -node 1 2 3 4 5 6 7 8 9 -dof 2 disp";
eval "recorder Node -file output/pwp_$Analysis_case.out -time -node 1 2 3 4 -dof 3 vel";
eval "recorder Element -file output/stress_$Analysis_case.out -time -ele 1 material 9 stress";  # s11 s22 s33 s12 eta	(EFFECTIVE stresses)
eval "recorder Element -file output/strain_$Analysis_case.out -time -ele 1 material 9 strain";  # e11 e22 gama12 (it is not e33 nor e12)
eval "recorder Element -file output/backbone_$Analysis_case.out -ele 1 material 9 backbone 100 1600";
####################################################################
# GRAVITY APPLICATION (elastic behavior)
model basic -ndm 2 -ndf 3
pattern Plain 1 Constant {
	load 3 [expr $loadbias*$vertPress*$xsize*$thickness*0.25] [expr -$vertPress*$xsize*$thickness*0.25] 0
	load 4 [expr $loadbias*$vertPress*$xsize*$thickness*0.25] [expr -$vertPress*$xsize*$thickness*0.25] 0
}
model basic -ndm 2 -ndf 2
pattern Plain 2 Constant {
	load 7 [expr $loadbias*$vertPress*$xsize*$thickness*0.5] [expr -$vertPress*$xsize*$thickness*0.5]
}

numberer RCM
system ProfileSPD
test NormDispIncr 1.e-5 50 2
constraints Plain; #Penalty 1.e18 1.e18
rayleigh [expr 1.0*$Mproprayleigh] 0.0 $Kinitproprayleigh 0.0
set gama 1.5;
set betta [expr pow($gama+0.5, 2)/4]; #[expr 1./4.]; # 0.25 for const accel, 1/6 for linear accel (in this case deltaT<period/pi)
integrator Newmark $gama $betta;
algorithm KrylovNewton
analysis VariableTransient

updateMaterialStage -material 1 -stage 0
updateMaterialStage -material 2 -stage 0
updateMaterialStage -material 3 -stage 0
updateMaterialStage -material 4 -stage 0
for {set i 1} {$i <= 100} {incr i 1} {
analyze 1 1000
set tCurrent [getTime]
puts "time = $tCurrent sec"
}
for {set i 1} {$i <= 100} {incr i 1} {
analyze 1 1000
set tCurrent [getTime]
puts "time = $tCurrent sec"
}
puts "Done Elastic"
updateMaterialStage -material 1 -stage 1
updateMaterialStage -material 2 -stage 1
updateMaterialStage -material 3 -stage 1
updateMaterialStage -material 4 -stage 1
for {set i 1} {$i <= 100} {incr i 1} {
analyze 1 100
set tCurrent [getTime]
puts "time = $tCurrent sec"
}

for {set i 1} {$i <= 200} {incr i 1} {
analyze 1 1000
set tCurrent [getTime]
puts "time = $tCurrent sec"
}
puts "Done Plastic"
puts "Gravity Done!"

####################################################################
# Adjust some fixities for shear loading

# remove surface pwp fixity
model basic -ndm 2 -ndf 3
remove sp 3 3
remove sp 4 3
equalDOF 3 4  3

####################################################################
loadConst -time 0.0
wipeAnalysis

####################################################################
# NOW APPLY LOADING SEQUENCE AND ANALYZE (plastic)

# Load control cyclic
if {$Analysis_case == "undrained_cyclic"} {

	# Define analysis parameters
	numberer RCM
	system ProfileSPD
	test NormDispIncr 1.e-5 50 0
	constraints Plain
	rayleigh $Mproprayleigh  0.0 $Kinitproprayleigh 0.0
	set gama 1.5; #0.5;	#gama=alfa in other texts
	set betta [expr pow($gama+0.5, 2)/4]; # 0.25 for const accel, 1/6 for linear accel (in this case deltaT<period/pi)
	integrator Newmark $gama $betta;
	algorithm KrylovNewton
	analysis VariableTransient
	
	model basic -ndm 2 -ndf 3
	pattern Plain 4 "Sine 0 [expr $period*1000] $period" {
		load 3 [expr $csr*$vertPress*$xsize*$thickness*0.25] 0 0
		load 4 [expr $csr*$vertPress*$xsize*$thickness*0.25] 0 0
	}
	model basic -ndm 2 -ndf 2
	pattern Plain 5 "Sine 0 [expr $period*1000] $period" {
		load 7 [expr $csr*$vertPress*$xsize*$thickness*0.5] 0
	}
	puts "$Analysis_case"
	while {[expr abs([nodeDisp 7 1]/$ysize)]<$target_shear_strain} {
	# while {[getTime]<0.25}
	analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
	puts "time = [getTime] sec"
	}
}

# Disp control monotonic
if {$Analysis_case == "undrained_monotonic"} {
	
	# Define analysis parameters
	numberer RCM
	system BandGeneral; #ProfileSPD
	test NormDispIncr 1.e-5 50 0; #1.e-5 50 0
	constraints Penalty 1.e10 1.e10; #for clay use 1.e5, for sand use 1.e15
	rayleigh $Mproprayleigh  0.0 $Kinitproprayleigh 0.0
	set gama 1.5; #0.5;	#gama=alfa in other texts
	set betta [expr pow($gama+0.5, 2)/4]; # 0.25 for const accel, 1/6 for linear accel (in this case deltaT<period/pi)
	integrator Newmark $gama $betta;
	algorithm KrylovNewton
	analysis VariableTransient
	
	model basic -ndm 2 -ndf 2
	pattern Plain 14 Linear {
		sp 3 1 [expr $ysize*$target_shear_strain]
		sp 4 1 [expr $ysize*$target_shear_strain]
		sp 7 1 [expr $ysize*$target_shear_strain]
	}
	puts "$Analysis_case"
	for {set i 1} {$i <= [expr $numSteps]} {incr i 1} {
	analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
	set tCurrent [getTime]
	puts "time = $tCurrent sec"
	}
}

# Disp cotrol monotonic
if {$Analysis_case == "drained_monotonic"} {
	
	# Define analysis parameters
	numberer RCM
	system BandGeneral; #ProfileSPD
	test NormDispIncr 1.e-5 50 0
	constraints Penalty 1.e18 1.e18
	rayleigh $Mproprayleigh  0.0 $Kinitproprayleigh 0.0
	set gama 1.5; #0.5;	#gama=alfa in other texts
	set betta [expr pow($gama+0.5, 2)/4]; # 0.25 for const accel, 1/6 for linear accel (in this case deltaT<period/pi)
	integrator Newmark $gama $betta;
	algorithm KrylovNewton
	analysis VariableTransient
	
	model basic -ndm 2 -ndf 3
	pattern Plain 15 Linear {
		sp 3 1 [expr $ysize*$target_shear_strain]
		sp 4 1 [expr $ysize*$target_shear_strain]
		sp 7 1 [expr $ysize*$target_shear_strain]
	}
	puts "$Analysis_case"
	for {set i 1} {$i <= [expr $numSteps]} {incr i 1} {
	analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
	set tCurrent [getTime]
	puts "time = $tCurrent sec"
	}
}

# Disp control cyclic for G/Gmax and Xi curves
if {$Analysis_case == "drained_cyclic"} {
	
	# Define analysis parameters
	numberer RCM
	system BandGeneral; #ProfileSPD
	test NormDispIncr 1.e-5 50 0
	constraints Penalty 1.e18 1.e18
	rayleigh $Mproprayleigh  0.0 $Kinitproprayleigh 0.0
	set gama 1.5; #0.5;	#gama=alfa in other texts
	set betta [expr pow($gama+0.5, 2)/4]; # 0.25 for const accel, 1/6 for linear accel (in this case deltaT<period/pi)
	integrator Newmark $gama $betta;
	algorithm KrylovNewton
	analysis VariableTransient
		
	model basic -ndm 2 -ndf 2
	pattern Plain 4 "Sine 0 1000 $period" {
		sp 3 1 [expr $ysize*0.0003/100.]
	}
	puts "$Analysis_case"
	for {set i 1} {$i <= $numSteps} {incr i 1} {
	analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
	set tCurrent [getTime]
	puts "time = $tCurrent sec"
	}
	remove loadPattern 4
	loadConst;	

	pattern Plain 5 "Sine 0 1000 $period" {
		sp 3 1 [expr $ysize*0.001/100.]
	}
	for {set i 1} {$i <= $numSteps} {incr i 1} {
	analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
	set tCurrent [getTime]
	puts "time = $tCurrent sec"
	}
	remove loadPattern 5
	loadConst

	pattern Plain 6 "Sine 0 1000 $period" {
		sp 3 1 [expr $ysize*0.003*0.01]
	}
	for {set i 1} {$i <= $numSteps} {incr i 1} {
	analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
	set tCurrent [getTime]
	puts "time = $tCurrent sec"
	}
	remove loadPattern 6
	loadConst

	pattern Plain 7 "Sine 0 1000 $period" {
		sp 3 1 [expr $ysize*0.01*0.01]
	}
	for {set i 1} {$i <= $numSteps} {incr i 1} {
	analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
	set tCurrent [getTime]
	puts "time = $tCurrent sec"
	}
	remove loadPattern 7
	loadConst

	pattern Plain 8 "Sine 0 1000 $period" {
		sp 3 1 [expr $ysize*0.03*0.01]
	}
	for {set i 1} {$i <= $numSteps} {incr i 1} {
	analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
	set tCurrent [getTime]
	puts "time = $tCurrent sec"
	}
	remove loadPattern 8
	loadConst

	pattern Plain 9 "Sine 0 1000 $period" {
		sp 3 1 [expr $ysize*0.1*0.01]
	}
	for {set i 1} {$i <= $numSteps} {incr i 1} {
	analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
	set tCurrent [getTime]
	puts "time = $tCurrent sec"
	}
	remove loadPattern 9
	loadConst

	pattern Plain 10 "Sine 0 1000 $period" {
		sp 3 1 [expr $ysize*0.3*0.01]
	}
	for {set i 1} {$i <= $numSteps} {incr i 1} {
	analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
	set tCurrent [getTime]
	puts "time = $tCurrent sec"
	}
	remove loadPattern 10
	loadConst

	pattern Plain 11 "Sine 0 1000 $period" {
		sp 3 1 [expr $ysize*1.0*0.01]
	}
	for {set i 1} {$i <= $numSteps} {incr i 1} {
	analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
	set tCurrent [getTime]
	puts "time = $tCurrent sec"
	}
	remove loadPattern 11
	loadConst

	pattern Plain 12 "Sine 0 1000 $period" {
		sp 3 1 [expr $ysize*3.0*0.01]
	}
	for {set i 1} {$i <= [expr $numSteps*1.25]} {incr i 1} {
	analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
	set tCurrent [getTime]
	puts "time = $tCurrent sec"
	}
}
####################################################################
set endT [clock seconds]
puts "Execution time: [expr $endT-$startT] seconds."
puts "Finished with cyclic loading"
wipe all;
wipe;  # flush ouput stream

DSS_material_properties.tcl file

#############################################################################
# PDMY03 input parameters materials 									    #
# By: Arash Khosravifar   January 2009										#
# Last edit:			  August 2018										#
#############################################################################

if {$matTag == 1} {
# defined variables for Sand (N1)60=35 "nonliq" (Mat 1)
set massDen 			2.06;   # (ton/m3)
set refG 				111.9e3;# (kPa)
set refB 				298.3e3;# (kPa)
set frinctionAng 		42.2;	# (degree)
set peakShearStrain 	0.1;
set refPress 			100.;   # (kPa)
set pressDependCoe 		0.5;
set phaseTransAng 		37.2;	# (degree)

set contraction_a 		0.001;  # Contraction rate. 
set contraction_b 		0.0;  	# fabric damage
set contraction_c 		0.8;  	# k_sigma effect
set contraction_d 		2.2; 	# ~CRR*[3/(1+2k0)]/2
set contraction_e 		0.0; 	
set dilation_a 			0.6;	 
set dilation_b 			3.0;	
set dilation_c 			-0.5;	# k_sigma effect

set liqParam1 			1.0;
set liqParam2 			0.0;

set noYieldSurf 		20;
# set void 				0.55;
# set cs1 				0.9;
# set cs2 				0.02;
# set cs3 				0.0;
set pa					100;   # (kPa)
set S0					1.73;  # (kPa)
}

if {$matTag == 2} {
# defined variables for Sand (N1)60=25 (Mat 2)
set massDen 			2.03;   #(ton/m3)
set refG 				94.6e3; # (kPa)
set refB 				252.6e3;# (kPa)
set frinctionAng 		35.8;
set peakShearStrain 	0.1;
set refPress 			100; 	# (kPa)
set pressDependCoe 		0.5;
set phaseTransAng 		30.8;

set contraction_a 		0.005;
set contraction_b 		1.0;
set contraction_c 		0.6;  
set contraction_d 		4.6; 	# ~CRR*[3/(1+2k0)]/2
set contraction_e 		-1.0;
set dilation_a 			0.45;
set dilation_b 			3.0;
set dilation_c 			-0.4; 

set liqParam1 			1.0;
set liqParam2 			0.0;

set noYieldSurf 		20;
# set void 				0.60;
# set cs1 				0.9;
# set cs2 				0.02;
# set cs3 				0.0;
set pa					100;  	# (kPa)
set S0					1.73; 	# (kPa)
}

if {$matTag == 3} {
# defined variables for Sand (N1)60=15 (Mat 3)
set massDen 			1.99; 	#(ton/m3)
set refG 				73.7e3; # (kPa)
set refB 				196.8e3;# (kPa)
set frinctionAng 		30.3;
set peakShearStrain 	0.1;
set refPress 			100; 	# (kPa)
set pressDependCoe 		0.5;
set phaseTransAng 		25.3;

set contraction_a 		0.012;
set contraction_b 		3.0;
set contraction_c 		0.4; 
set contraction_d 		9.0; 	# ~CRR*[3/(1+2k0)]/2
set contraction_e 		0.0; 
set dilation_a 			0.3; 
set dilation_b 			3.0;
set dilation_c 			-0.3;

set liqParam1 			1.0;
set liqParam2 			0.0;

set noYieldSurf 		20;
# set void 				0.67;
# set cs1 				0.9;
# set cs2 				0.02;
# set cs3 				0.0;
set pa					100;  	# (kPa)
set S0					1.73; 	# (kPa)
}

if {$matTag == 4} {
# defined variables for Sand (N1)60=5 (Mat 4)
set massDen 			1.94; 	#(ton/m3)
set refG 				46.9e3; # (kPa)
set refB 				125.1e3;# (kPa)
set frinctionAng 		25.4;
set peakShearStrain 	0.1;
set refPress 			100; 	# (kPa)
set pressDependCoe 		0.5;
set phaseTransAng 		20.0;

set contraction_a 	 	0.03;
set contraction_b 	 	5.0;
set contraction_c 	 	0.2; 
set contraction_d 	 	16.0; 	# ~CRR*[3/(1+2k0)]/2
set contraction_e 	 	2.0; 
set dilation_a 		 	0.15;
set dilation_b 		 	3.0;
set dilation_c 		 	-0.2;

set liqParam1 			1.0;
set liqParam2 			0.0;

set noYieldSurf 		20;
# set void 				0.76;
# set cs1 				0.9;
# set cs2 				0.02;
# set cs3 				0.0;
set pa					100;  	# (kPa)
set S0					1.73; 	# (kPa)
}

MATLAB Plotting File

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Post processor for 9-4quadUP element                  %
% By: Arash Khosravifar    January 2009                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc
clear all
close all

% Type of analysis. Options:
% undrained_cyclic
% undrained_monotonic
% drained_cyclic
% drained_monotonic
Analysis_case = 'undrained_cyclic';   

% Loading files
input_disp1      = load(['output\disp1_',Analysis_case,'.out']);
input_disp2      = load(['output\disp2_',Analysis_case,'.out']);
input_stress    = load(['output\stress_',Analysis_case,'.out']);
input_strain    = load(['output\strain_',Analysis_case,'.out']);
input_pwp       = load(['output\pwp_',Analysis_case,'.out']);
input_backbone  = load(['output\backbone_',Analysis_case,'.out']);


% Displacements
Disp_1 = input_disp1(:,2:end);
Disp_2 = input_disp2(:,2:end);
% Stresses  compression:+   tension:-
s_11  = -input_stress(:,2);
s_22  = -input_stress(:,3);
s_33  = -input_stress(:,4);
s_12  = input_stress(:,5);
eta   = input_stress(:,6);  % defined in the manual - not used

tau_oct  = 1/3*sqrt((s_11-s_22).^2 + (s_11-s_33).^2 + (s_22-s_33).^2 + 6*s_12.^2); % Octahedral shear stress
q     = 3/sqrt(2)*tau_oct; %1/3*sign(s_12).*sqrt((s_11-s_22).^2 + (s_11-s_33).^2 + (s_22-s_33).^2 + 6*s_12.^2);
pwp   = (input_pwp(:,2)+input_pwp(:,3)+input_pwp(:,4)+input_pwp(:,5))/4;
p     = (s_11 + s_22 + s_33)/3;
p_tot = p + pwp;
time  = input_pwp(:,1);

%Note: The stresses are already effective. no need to subtract
% Strains   Upward (dilation):-   Downward (contraction):+
e_11  = -input_strain(:,2);
e_22  = -input_strain(:,3);
gama_12 = input_strain(:,4);
e_33  = zeros(size(e_11));       % plain strain
e_v   = e_11 + e_22 + e_33;      % e_33 is zero
gama_oct = 2/3*sqrt((e_11-e_22).^2+(e_22-e_33).^2+(e_11-e_33).^2+6*(gama_12/2).^2); 

% CSL line
cs1 = 0.9;
cs2 = 0.1;
cs3 = 1.0;
void_i = 0.8;     % initial void ratio
void = void_i - e_22*(1+void_i);

% Number of steps in "consolidation phase"
consol_step = 500;
total_step  = size(gama_12,1);

% Ru
ru_1 = pwp./max(p,0.001);
ru_2 = pwp/p(consol_step);
ru_3 = pwp/s_22(consol_step);

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(Analysis_case,'drained_cyclic') == 1
    period      = 10;
    cycle       = time/period;

    % G/Gmax
    NoStepsCycle  = 2000;
    NoCycleEach   = 2;
    corr_steps = [consol_step+NoStepsCycle/4:NoStepsCycle*NoCycleEach : consol_step+9*NoStepsCycle*NoCycleEach];
    G_sec = s_12(corr_steps)./gama_12(corr_steps);      %Secant stiffness
    G_max = s_12(consol_step+1)/gama_12(consol_step+1);              %Secant stiffness at the first load incr.
    G_max_oct = tau_oct(consol_step+1)/gama_oct(consol_step+1);      
    % Note that the Gmax,r assigned for the material by the user is not the
    % conventional Gmax (s_12/gama_12), but it is the OCTAHEDRAL one
    % (tau_oct/gama_oct)

    % EPRI 1993
    GratioEPRI6 =  [1,0.98,0.919,0.751,0.519,0.271,0.121,0.045,0.024];  % EPRI1993 for 0-6m corresponding to strain=[0.000003,0.00001,0.00003,0.0001,0.0003,0.001,0.003,0.01,0.03]
    GratioEPRI76 = [1,0.999,0.982,0.9,0.752,0.502,0.282,0.124,0.07];    % EPRI1993 for 36-76m
    gama_EPRI  =  [0.000003,0.00001,0.00003,0.0001,0.0003,0.001,0.003,0.01,0.03];
    gama_EPRI_oct = 2/3*sqrt((e_11(500)-e_22(500)).^2+(e_22(500)-e_33(500)).^2+(e_11(500)-e_33(500)).^2+6*(gama_EPRI/2).^2);

    % Use backbone record
    backbone = input_backbone(1,:)';
    gama_backbone_100_oct = input_backbone(5:4:end);
    G_backbone_100 = input_backbone(6:4:end);
    gama_backbone_1600_oct = input_backbone(7:4:end);
    G_backbone_1600 = input_backbone(8:4:end);

    % Xi
    for ii = 1:length(corr_steps)
        hysteresis = trapz(gama_12(corr_steps(ii):corr_steps(ii)+NoStepsCycle),s_12(corr_steps(ii):corr_steps(ii)+NoStepsCycle));
        triangle   = 1/2*gama_12(corr_steps(ii))*s_12(corr_steps(ii));
        Xi(ii)     = hysteresis/(4*pi*triangle)*100;    %(%)
    end
    XiEPRI6 =  [0.791,1.27,2.28,4.26,7.75,14.74,22.03,27.12,30.21];     % (%) EPRI1993 for 0-6m corresponding to strain=[0.000003,0.00001,0.00003,0.0001,0.0003,0.001,0.003,0.01,0.03]
    XiEPRI76 = [0.591,0.88,1.271,2.46,4.551,8.34,14.13,21.22,25.31];    % (%) EPRI1993 for 36-76m
    
	% Plots
    figure
    semilogx(gama_12(corr_steps),G_sec/G_max,'-ob','LineWidth',2);
    hold on; grid on;
    semilogx(gama_EPRI,GratioEPRI6,'--r','LineWidth',2);
    semilogx(gama_EPRI,GratioEPRI76,'--r','LineWidth',2);
    xlabel('\gamma_{12}');
    ylabel('G_{sec}/G_{max}');
    legend('PDMY03 - \sigma_{vo}=100 KPa','EPRI1993 0-6m','EPRI1993 36-76m')

    % Xi
    figure
    semilogx(gama_12(corr_steps),Xi,'-ob','LineWidth',2);
    hold on; grid on;
    semilogx(gama_12(corr_steps),XiEPRI6,'--r','LineWidth',2);
    semilogx(gama_12(corr_steps),XiEPRI76,'--r','LineWidth',2);
    xlabel('\gamma_{12}');
    ylabel('Equivalent damping ratio (%)');
    legend('PDMY03 - \sigma_{vo}=100 KPa','EPRI1993 0-6m','EPRI1993 36-76m')
    
    figure
    subplot(2,1,1);
    plot(gama_12,s_12/s_22(consol_step),'blue');
    xlabel('\gamma_{12}');
    ylabel('CSR = (\sigma_{12}/\sigma''_{vc})');
    title(['\sigma''_{vc} = ',num2str(s_22(consol_step)),' kPa']);

    subplot(2,1,2);
    plot(cycle(consol_step+1:end),gama_12(consol_step+1:end),'blue');
    xlabel('No. of cycles');
    ylabel('\gamma_{12}');
        
   
end;    % (end of drained_cyclic)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(Analysis_case,'undrained_cyclic') == 1
    period      = 1;
    cycle       = time/period;

    % find "Liquefaction Triggering" time
    gama_consol = gama_12(consol_step);

    T_1  = find(abs(gama_12-gama_consol)>0.01,1);
    if any(T_1) == 0
        T_1 = consol_step+1;
    end
    N_1 = cycle(T_1)

    T_3 = find(abs(gama_12-gama_consol)>0.03,1);
    if any(T_3) == 0
        T_3 = consol_step+1;
    end
    N_3 = cycle(T_3)

    % T_r = find(ru_3(consol_step:total_step)>0.98,1);
    T_r = find(ru_3>0.98,1);
    if any(T_r) == 0
        T_r = consol_step+1;
    end
    N_r = cycle(T_r)

    % Plotting
    figure
    subplot(2,2,1);
    plot(gama_12,s_12/s_22(consol_step));grid on; hold on;xlabel('\gamma_{12}');ylabel('CSR = (\sigma_{12}/\sigma''_{vc})')
    plot(gama_12(T_1),s_12(T_1)/s_22(consol_step),'g.','MarkerSize',15);
    plot(gama_12(T_3),s_12(T_3)/s_22(consol_step),'r.','MarkerSize',15);
    plot(gama_12(T_r),s_12(T_r)/s_22(consol_step),'m.','MarkerSize',15);
    title(['CSR = ',num2str(max(abs(s_12))/s_22(consol_step),2),'   and   \sigma''_{vc} = ',num2str(s_22(consol_step)),' kPa']);
%     xlim([-0.1,0.1]);

    subplot(2,2,2);plot(s_22/s_22(consol_step),s_12/s_22(consol_step));grid on; hold on;xlabel('vertical effective stress ratio (\sigma''_{v}/\sigma''_{vc})');ylabel('CSR = (\sigma_{12}/\sigma''_{vc})');
    plot(s_22(T_1)/s_22(consol_step),s_12(T_1)/s_22(consol_step),'g.','MarkerSize',15);
    plot(s_22(T_3)/s_22(consol_step),s_12(T_3)/s_22(consol_step),'r.','MarkerSize',15);
    plot(s_22(T_r)/s_22(consol_step),s_12(T_r)/s_22(consol_step),'m.','MarkerSize',15);
    xlim([min(0,min(s_22/s_22(consol_step))),max(s_22)/s_22(consol_step)]);
    title(['CSR = ',num2str(max(abs(s_12))/s_22(consol_step),2),'   and   \sigma''_{vc} = ',num2str(s_22(consol_step)),' kPa']);

    subplot(2,2,3);plot(cycle(consol_step+1:end),pwp(consol_step+1:end));grid on;ylabel('PWP (kPa)');xlabel('No. of cycles');
    hold on;
    plot(cycle(T_1),pwp(T_1),'g.','MarkerSize',15);
    plot(cycle(T_3),pwp(T_3),'r.','MarkerSize',15);
    plot(cycle(T_r),pwp(T_r),'m.','MarkerSize',15);
    ylim([-10,100]);

    subplot(2,2,4);plot(cycle(consol_step+1:end),gama_12(consol_step+1:end));grid on;xlabel('No. of cycles');ylabel('\gamma_{12}');

end;    % (end of undrained_cyclic)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(Analysis_case,'drained_monotonic') == 1

    G_max = s_12(consol_step+1)/gama_12(consol_step+1);
    G_max_oct = tau_oct(consol_step+1)/gama_oct(consol_step+1);      
    
    
    figure
    subplot(1,2,1);hold on; grid on;
    plot(gama_12(consol_step+1:end),s_12(consol_step+1:end));
    plot(gama_oct(consol_step+1:end),tau_oct(consol_step+1:end),'r');
    legend('\sigma_{12} vs \gamma_{12}','\tau_{oct} vs \gamma_{oct}',4);
    xlabel('\gamma_{12} and \gamma_{oct}');ylabel('\sigma_{12} and \tau_{oct}');
    ylim([0,60]);
    title(['\sigma''_{vc} (kPa) = ',num2str(s_22(consol_step))]);

    subplot(1,2,2);grid on;hold on
    plot(p(consol_step+1:end),tau_oct(consol_step+1:end),'r');
    plot(s_22(consol_step+1:end),s_12(consol_step+1:end),'b')
    ylim([0,60]);
   
    slope_oct = tau_oct(end)/p(end);
    Phi_oct_back_calculated = asin(3*slope_oct/(2*sqrt(2)+slope_oct))/3.1416*180;
    slope_conventional = s_12(end)/s_22(end);
    Phi_conventional_back_calculated = atan(slope_conventional)/3.1416*180;

    plot([0,max(p)*1.2],[0,(max(p)*1.2)*slope_oct],'--r')
    plot([0,max(s_22)*1.2],[0,(max(s_22)*1.2)*slope_conventional],'--b')
    xlabel('\sigma''_{v} and p'' (kPa)');ylabel('\sigma_{12} and \tau_{oct}');
    xlim([0,max(p)*1.2]);
    legend(['\phi_{oct} = ',num2str(Phi_oct_back_calculated)],['\phi_{DSS} = ',num2str(Phi_conventional_back_calculated)],2)

end;    % (end of drained_monotonic)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(Analysis_case,'undrained_monotonic') == 1
    figure
    subplot(2,2,1);plot(gama_12(consol_step+1:end),s_12(consol_step+1:end));
    grid on;
    xlabel('\gamma_{12}');
    ylabel('\sigma_{12} (kPa)');
    
    subplot(2,2,2);plot(s_22(consol_step+1:end),s_12(consol_step+1:end));
    grid on; hold on;
    xlim([0,max(s_22)])
    plot([0,s_22(end)],[0,s_12(end)],'r')
    xlabel('\sigma''_{v} (kPa)');
    ylabel('\sigma_{12} (kPa)');
    slope_conventional = s_12(end)/s_22(end);
    Phi_conventional_back_calculated = atan(slope_conventional)/3.1416*180;
    legend(['\phi_{DSS} = ',num2str(Phi_conventional_back_calculated)],2)

    subplot(2,2,3);plot(gama_12,pwp);grid on;
    xlabel('\gamma_{12}')
    ylabel('pwp (kPa)')
    subplot(2,2,4)
    plot(gama_12,e_v-e_v(consol_step)); grid on;
    set(gca,'YDir','rev');
    xlabel('\gamma_{12}')
    ylabel('Vomuletric strain \epsilon_v - \epsilon_o');

end;    % (end of undrained_monotonic')
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Return to:



Retrieved from "https://opensees.ist.berkeley.edu/wiki/index.php?title=PDMY03_elementdriver&oldid=12753"