Simple shear test Model

[ZNizamani]

Hi everyone,

I am trying to model DSS cyclic behavior on sand using load controlled method. So far I have good results which by the way resemble with un-drained conditions but when I do drained analysis I am having some errors. I think I am not using proper drained loading condition If something else is wrong you are more welcome to discuss.

Please, Could any one help me to solve this issue?

Kind regards,
Zubair


Code

wipe
wipe all
set startT [clock seconds]

set massDen 2.23; #(ton/m3)
set refG 114684.3356 ; # (kPa)
set refB 299078.7577; # (kPa)
set frinctionAng 41;
set peakShearStrain 0.1;
set refPress 92.6; # (kPa)
set pressDependCoe 0.5;
set phaseTransAng 20;

set contractionParam1 0.035;
set contractionParam2 5.20;
set contractionParam3 0.2;

set dilationParam1 0.045;
set dilationParam2 1.50;
set dilationParam3 -0.170;

set liqParam1 0.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 c 0.1; # (kPa)

set vertPress [expr 1.0*92.6]; # kPa vertical confining pressure
set loadbias 0.0; # Alpha = tau_xy/s'vo

set matTag 1;

#Analysis Selection
set Analysis_case "drained_cyclic"

# 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-3; # 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 5.0/100.]
set csr 0.25;
set period 1;
set deltaT 0.005;
set kPerm $kundrain
}

if {$Analysis_case == "drained_cyclic"} {
set target_shear_strain [expr 5.0/100.]
set csr 0.1;
set period 10.;
set deltaT 0.005;
set numSteps 4000;

set kPerm $kdrain

set frinctionAng 41;
set peakShearStrain 0.1;
set refPress 92.6; # (kPa)
set pressDependCoe 0.5;
set phaseTransAng 20;

set kPerm $kdrain
set contractionParam1 0.0;
set contractionParam2 0.0;
set contractionParam3 0.0;
set dilationParam1 0.0;
set dilationParam2 0.0;
set dilationParam3 0.0;
set liqParam1 0.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 c 0.1; # (kPa)
}
# Define Rayleigh Damping
set rayleigh_damping2 0.001;
set rayleigh_damping1 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

################# define material ##################################

nDMaterial PressureDependMultiYield02 $matTag 2 $massDen $refG $refB $frinctionAng \
$peakShearStrain $refPress $pressDependCoe $phaseTransAng \
$contractionParam1 $contractionParam3 $dilationParam1 $dilationParam3 \
$noYieldSurf $contractionParam2 $dilationParam2 $liqParam1 $liqParam2 \
$void $cs1 $cs2 $cs3 $pa $c;

################# define Element ##################################
# 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 g12 (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 $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
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
updateMaterials -material 2 refB [expr $refG*2/3.];

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 && [getTime]<200} {
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
puts "time = [getTime] sec"
}
}

# Disp control cyclic
if {$Analysis_case == "drained_cyclic"} {


# Define analysis parameters
numberer RCM
system BandGeneral
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 6 "Sine 0 [expr $period*1000] $period" {
sp 3 1 [expr $csr*$vertPress*$xsize*$thickness*0.25]
sp 4 1 [expr $csr*$vertPress*$xsize*$thickness*0.25]
}
model basic -ndm 2 -ndf 2
pattern Plain 7 "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 && [getTime]<200} {
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
puts "time = [getTime] sec"
}
}


####################################################################
set endT [clock seconds]
puts "Execution time: [expr $endT-$startT] seconds."
puts "Finished with cyclic loading"
wipe all
wipe; # flush output stream

[skamalzare]

Hi,

The third DOF with quadUP elements is the PWP. For your drained analysis, you need to assign zero PWP to all your nodes (e.g. fix #NodeNo 1 1 1).

In the section: "# Adjust some fixities for shear loading" -- You remove the PWP constraints to apply shear loading. You probably shouldn't do that.

Bests,
Soheil

[ZNizamani]

Hi Soheil,
Thank you for your kind response,

I changed fixity to drained condition using same command as you said
fix 1 1 1 1
fix 2 1 1 1
fix 3 0 0 1
fix 4 0 0 1
and I used drained syntax as follows

# Load controlled
if {$Analysis_case == "drained_cyclic"} {


# Define analysis parameters
numberer RCM
system BandGeneral
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 6 "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 7 "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 && [getTime]<300} {
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
puts "time = [getTime] sec"
}
}
I get the response of 0 shear strain and change in vertical stress which is in actual it should not. My shear stress-strain loop doesn't look any right. Would you please look at the code again? I will be thankful.

Best regards,
Zubair

[skamalzare]

Zubair,

The fact that you are getting almost zero change in vertical stress is good! It is actually what is expected. In DSS tests, we keep the vertical total stress constant and only apply shear stress to the sample. As your test is drained no PWP will be generated. Therefore, effective stresses are the same as total stresses and they should not change during the test.

I ran your code with the modifications you suggested. I am monitoring +-0.02% shear strain with 0 change in vertical stresses. The code looks good to me.

In section "Adjust some fixities for shear loading", you are removing 3rd DOF (PWP) constrains from your top nodes (i.e. remove sp 3 3). This has minor effects on your results, but in general, it is not correct for drained conditions. It resembles closing drainage valves on top of your sample (needed for undrained simulations but not in drained conditions).

hope this could help.

Best Regards,
Soheil

[ZNizamani]

Dear Soheil,

It feels really great hearing from you. I may ask you some questions about this Model.
-- I ran this code as you said about drainage pattern like, remove sp 3 for un-drained condition and for drained condition let it be as it is like;
fix 1 1 1 0
fix 2 1 1 0
fix 3 0 0 1
fix 4 0 0 1

My model runs fine, but the problem I get is very small shear strain as you said -+ 0.02%. If I change CSR from 0.25 to 0.35 or even 0.4 it doesn't affect the results and if I change density it also doesn't affect the results. In the results, model doesn't generate Pore-water pressure and no change in vertical stress which is good but still there is NO volumetric strain is produced. In drained condition there must be volumetric strain. Am I making some syntax errors or drained test would be performed some other way? Please help me

warm regards,
Zubair

[skamalzare]

Zubair,

I am not sure if you understood my last message correctly. Why should we block drainage path from top of our sample if it is a drained test?

In your drained simulations you need to assign zero PWP to all your nodes (i.e. fix 3rd DOF).
In your undrained simulations, you need to let the nodes have PWP.

Therefore, your boundary conditions need to be changed for drained simulations.

You are not getting large shear strains, because your material is very strong. Look at OpenSeesWiki to have an idea about max and min parameters you might want to use. Generally speaking friction angle 41 is too much, dialation3 cannot be negative and ...

Bests,
Soheil

[sadeghihabib]

Dear All
I am modeling 3D Consolidated undrained cyclic simple shear test the same as the experimental tests conducted by Ishihara-Youshimine 1992

The Script below is my script
the problem is that after applying cyclic loading no porewter pressure generates in the single SSPbrickUP element
Could anybody help me with my script

wipe

# ------------------------ #
# Test Specific parameters #
# ------------------------ #
# Confinement Stress

set Atop [expr 1.0*1.0]
set Aside [expr 1.0*1.0]

set Pl 200.0 ; #kPa
set Fl [expr $Pl*$Aside/4] ; #Lateral element surface unit force(kN) (positive)
set nFl [expr -1*$Fl] ; #Lateral element surface unit force(kN) (negative)

set Pv 200.0 ; #kPa
set Fv [expr $Pv*$Atop/4] ; #vertical element surface unit force(kPa) (from bottom to top dir) (positive)
set nFv [expr -1*$Fv] ; #vertical element surface unit force(kPa) (from top to bottom dir) (negative)


# Permeablity
set perm 1.0e-10
# Initial void ratio
set vR 0.8

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

# HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
# HHHHHHHHHHHHHHHHHHHHHHHHHHHCreate ModelHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
# HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH

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

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

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


equalDOF 5 6 1 2 3 4
equalDOF 5 8 1 2 3 4
equalDOF 8 7 1 2 3 4



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

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

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


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

# Apply confinement pressure

pattern Plain 1 {Series -time {0 10000 1e10} -values {0 1 1} -factor 1} {
load 1 $nFl 0.0 0.0 0.0
load 2 $nFl $nFl 0.0 0.0
load 3 0.0 $nFl 0.0 0.0
load 4 0.0 0.0 0.0 0.0
load 5 $nFl $Fl $nFv 0.0
load 6 $nFl $nFl $nFv 0.0
load 7 $Fl $nFl $nFv 0.0
load 8 $Fl $Fl $nFv 0.0
}

analyze 100 100

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

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

# Apply cyclic shear loading


set PCSRx 0.159
set MaxTowx [expr $PCSRx*$Pv]
set MaxShearForcex [expr $MaxTowx*$Atop/4]
set periodx 4.0
timeSeries Trig 20 0.0 20.0 2.0 -factor $MaxShearForcex

pattern Plain 2 20 {

load 5 1.0 0.0 0.0 0.0
load 6 1.0 0.0 0.0 0.0
load 7 1.0 0.0 0.0 0.0
load 8 1.0 0.0 0.0 0.0

}

set PCSRy 0.159
set MaxTowy [expr $PCSRy*$Pv]
set MaxShearForcey [expr $MaxTowy*$Atop/4]
set periody 4.0
timeSeries Trig 30 0.0 20.0 2.0 -factor $MaxShearForcey

pattern Plain 3 30 {

load 5 0.0 1.0 0.0 0.0
load 6 0.0 1.0 0.0 0.0
load 7 0.0 1.0 0.0 0.0
load 8 0.0 1.0 0.0 0.0

}


# Set number and length of (pseudo)time steps
set dtDyn 0.00100
set analT 20.0000
set incnum [expr floor($analT/$dtDyn)]

puts "Start analysis"
set startT [clock seconds]

for {set numIncr 1} {$numIncr <= $incnum } {incr numIncr 1} {
puts "##### Seismic excitation step : $numIncr #####";
set OK [analyze 1 $dtDyn [expr $dtDyn/100.0] $dtDyn 20]
if {$OK != 0} {
error " ### analysis Faild ### "
}
}


set endT [clock seconds]
puts "loading analysis execution time: [expr $endT-$startT] seconds."

wipe