Difference between revisions of "PDMY03 elementdriver"
From OpenSeesWiki
Jump to navigationJump to search| Line 441: | Line 441: | ||
wipe; # flush ouput stream | wipe; # flush ouput stream | ||
</syntaxhighlight> | </syntaxhighlight> | ||
| − | |||
| − | |||
| Line 847: | Line 845: | ||
end; % (end of undrained_monotonic') | end; % (end of undrained_monotonic') | ||
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | ||
| − | |||
| − | |||
</syntaxhighlight> | </syntaxhighlight> | ||
Revision as of 23:48, 10 June 2019
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 4; # 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.0; # 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.085;
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 \
$contraction_a $contraction_c $dilation_a $dilation_c \
$noYieldSurf $contraction_b $dilation_b $liqParam1 $liqParam2 \
$void $cs1 $cs2 $cs3 $pa $S0 0 1. 0. 0. 0. $contraction_d 3. $contraction_e 2.;
####################################################################
# 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:
- NDMaterial Command
- UC San Diego soil models (Linear/Nonlinear, dry/drained/undrained soil response under general 2D/3D static/cyclic loading conditions (please visit UCSD for examples)
- UC San Diego Saturated Undrained soil
- Element Command
- UC San Diego u-p element (saturated soil)
- Related References
