Direct Simple Shear Test Simulation

[mirkomorosi]

Dear all,
I'm trying to simulate a Direct Simple Shear test simulation with opensees, enforcing drained condition to an isotropically consolidated soil sample.
Quad finite element has been used to implement this analysis, and PDMY02 material defines the soil constitutive model properties. The simulation is split into two stage: isotropic consolidation (gravity analysis) and dynamic shearing.
During the isotropic consolidation the two upper nodes of the Quad finite element were free to move, whereas the bottom nodes were fixed. Since a sufficiently high value of the Bulk Modulus has been assigned (B=1000*G), even though a 2D finite element was implemented, at the end of the gravity analysis (downward vertical forces at the two upper nodes) the isotropic stress stage has been reached (s’11=s’22=s’33=p’=300kPa).

After this stage dynamic loading starts (shearing, strain controlled, sinusoidal horizontal displacement at the upper nodes). Several loops have been implemented, characterized by different maximum horizontal peak amplitude (2000 steps for each cycle, with two cycle for each different amplitude). At the beginning of the dynamic loading the vertical strains do not practically change with respect to the final stage of isotropic consolidation, and normal stresses kept the isotropic state.
After the smaller peak displacement loops the output starts to be instable, meaning that around 20600 step of the analysis (beginning of the gamma,max=1*10-3 loop) normal stresses start peaking at much higher levels and afterwards dropping at zero, with simultaneous similar oscillation of the vertical upper node strains (node 3 and 4). At the positive compressive peak stress correspond the apparent shear failure of the soil sample (the ratio between deviatoris shear stress and peak shear strength eta,r goes up to one) even though the soil sample stress and strain state characterizing the previous step would not justify this failure.
During this onset of unstable phase normal stresses and strains keeps oscillating up to a point where normal stress drop to zero.

Since similar simple shear tests have been already modeled with other software and due to available experimental test results it seems that this instability captured by opensees dynamic analysis would not be correctly reflecting the real test behavior.
In order to deeper understand the hypothesis and possible mistakes that I made during this implementation I have two question to ask.

QUESTION 1.
The first one deals with the way constraints and restraints have been imposed. In order to avoid vertical strain oscillation during the dynamic stage I’ve tried to fix the vertical settlement get after the gravity load application throughout the dynamic stage, simply implemented as follows:
(…)
# Adjust some fixities for shear loading
fix 3 0 1
fix 4 0 1
equalDOF 3 4 1 2
(…)
This syntax let the vertical strains of the upper nodes (3 and 4) drop almost to zero, with simultaneous normal stresses drop to zero level. Thus vanishing the isotropic initial state.
To this aim I’ve also tried to apply simultaneous horizontal displacement sinusoidal pattern and zero vertical settlement during the displacement controlled analysis, but without getting the expected results (== zero vertical settlement during dynamic analysis).

Could you please tell me what I did wrong? And what should be the correct implementation to fix certain DoFs during the following dynamic stage of the analysis?

QUESTION 2.
Since I have assigned an extremely high value to the Bulk Modulus, I want to check which is the effect of updating this modulus to its correct value before running the dynamic analysis. The syntax used for this implementation is the following:
(…)
set ctr 10000
parameter $ctr element 1 material $matTag refB
updateParameter $ctr 279.1e3
(…)
Where refB refers to the reference Bulk Modulus, input parameter of PDMY02 material constitutive model:
(…)
nDMaterial PressureDependMultiYield02 $matTag 2 $massDen $refG $refB $frinctionAng \
$peakShearStrain $refPress $pressDependCoe $phaseTransAng \
$contractionParam1 $contractionParam3 $dilationParam1 $dilationParam3 -19\
1.00e-6 0.998 1.78e-6 0.997 3.16e-6 0.994 5.62e-6 0.990 1.00e-5 0.982\
1.78e-5 0.970 3.16e-5 0.948 5.62e-5 0.914 1.00e-4 0.862 1.78e-4 0.789\
3.16e-4 0.698 5.62e-4 0.593 1.00e-3 0.484 1.78e-3 0.380 3.16e-3 0.286\
5.62e-3 0.207 1.00e-2 0.142 1.78e-2 0.092 3.16e-2 0.057\
$contractionParam2 $dilationParam2 $liqParam1 $liqParam2 \
$void $cs1 $cs2 $cs3 $pa $c;
}
(…)

I’ve followed the example detailed in this web-page http://opensees.berkeley.edu/wiki/index ... oil_Column
as detailed in “Material and Elements definition” Section. I’ve also try to follow an alternative procedure (the syntax is in this case described in this link: http://opensees.berkeley.edu/OpenSees/m ... l/1560.htm) but neither the former or the latter method seem able to correctly implement the parameter updating.

Could you tell me how to do this? I mean, how to update the reference bulk modulus of PDMY02 model?

In the following my input file is attached:

wipe
wipe all
set startT [clock seconds]

######################### MAIN USER INPUTS #################################
set matTag 2; #
set matType "PDMY02"; #
set vertPress [expr 1.0*300]; # kPa vertical confining pressure
set loadbias 0.0; # Alpha = tau_xy/s'vo
set Analysis_case "drained_cyclic"
# drained_cyclic: applies 2 cycles for each strain (9 shear strain)
######################### END OF MAIN USER INPUTS ###########################

# material properties
set massDen 1.895; # (ton/m3)
set refG 128831; # (kPa)
#set refB 279.1e3; # (kPa)
set refB 128.8e6; # (kPa)
set frinctionAng 34; # (degree)
set peakShearStrain 0.0316;
set refPress 300.; # (kPa)
set pressDependCoe 0.5;
set phaseTransAng 27; # (degree)

# 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-5; # permeability for undrained loading
set gravY 0.0;
set gravX 0.0;
set press 0.0;

set period 10.;
set deltaT 0.005;
set numSteps 4000;

set kPerm $kdrain
set phaseTransAng $frinctionAng
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;

#---RAYLEIGH DAMPING PARAMETERS
set pi 3.141592654
# damping ratio
set damp 0.002
# lower frequency
set omega1 0
# upper frequency
set omega2 [expr 2*$pi*1]
# damping coefficients
set a0 0
set a1 [expr 2*$damp/($omega1 + $omega2)]
puts "damping coefficients: a_0 = $a0; a_1 = $a1"

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

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

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

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

if {$matType == "PDMY02"} {
# SAND (PDMY02)
# G/Gmax == SUMDES fitting hr=0.3 (p'=300kPa, Dense Sand)

# Darendeli 2001 Clean Sands & YEE et al. method to tau,ff=80kPa (tau,ff=90kPa;gamma,1=1*e-01%)
nDMaterial PressureDependMultiYield02 $matTag 2 $massDen $refG $refB $frinctionAng \
$peakShearStrain $refPress $pressDependCoe $phaseTransAng \
$contractionParam1 $contractionParam3 $dilationParam1 $dilationParam3 -19\
1.00e-6 0.998 1.78e-6 0.997 3.16e-6 0.994 5.62e-6 0.990 1.00e-5 0.982\
1.78e-5 0.970 3.16e-5 0.948 5.62e-5 0.914 1.00e-4 0.862 1.78e-4 0.789\
3.16e-4 0.698 5.62e-4 0.593 1.00e-3 0.484 1.78e-3 0.380 3.16e-3 0.286\
5.62e-3 0.207 1.00e-2 0.142 1.78e-2 0.092 3.16e-2 0.057\
$contractionParam2 $dilationParam2 $liqParam1 $liqParam2 \
$void $cs1 $cs2 $cs3 $pa $c;
}

####################################################################
element quad 1 1 2 3 4 $thickness "PlaneStrain" $matTag $press 0.0 [expr $massDen*$gravX] [expr $massDen*$gravY]
####################################################################
# Recorders
file mkdir output
eval "recorder Node -file outputDR/disp1_$Analysis_case.out -time -node 1 2 3 4 -dof 1 disp";
eval "recorder Node -file outputDR/disp2_$Analysis_case.out -time -node 1 2 3 4 -dof 2 disp";
eval "recorder Node -file outputDR/pwp_$Analysis_case.out -time -node 1 2 3 4 -dof 3 vel";
eval "recorder Element -file outputDR/stress_$Analysis_case.out -time -ele 1 material 1 stress"; # s11 s22 s33 s12 eta (EFFECTIVE stresses)
eval "recorder Element -file outputDR/strain_$Analysis_case.out -time -ele 1 material 1 strain"; # e11 e22 g12 (it is not e33 nor e12)
eval "recorder Element -file outputDR/backbone_$Analysis_case.out -ele 1 material 1 backbone 300";
####################################################################
# GRAVITY APPLICATION (elastic behavior)
model basic -ndm 2 -ndf 2
pattern Plain 1 Constant {
load 3 [expr $loadbias*$vertPress*$xsize*$thickness*0.5] [expr -$vertPress*$xsize*$thickness*0.5]
load 4 [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 $a0 0.0 $a1 0.0
#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
updateMaterialStage -material 2 -stage 0
updateMaterialStage -material 3 -stage 0
updateMaterialStage -material 4 -stage 0
updateMaterialStage -material 5 -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
updateMaterialStage -material 5 -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
#fix 3 0 1
#fix 4 0 1
#equalDOF 3 4 1 2
####################################################################
loadConst -time 0.0
wipeAnalysis

#updateMaterials (-material?!?!) $matTag refB 279.1e1
set ctr 10000
parameter $ctr element 1 material $matTag refB
updateParameter $ctr 279.1e3
#updateParameter -material $matTag -refB 279.1e3
####################################################################
# NOW APPLY LOADING SEQUENCE AND ANALYZE (plastic)

# Load control cyclic
# Disp control cyclic for G/Gmax and Xi curves

# Define analysis parameters
numberer RCM
system BandGeneral; #ProfileSPD
test NormDispIncr 1.e-5 50 0
constraints Penalty 1.e18 1.e18
rayleigh $a0 0.0 $a1 0.0
#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.0001/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 4
loadConst
#I can actually comment this out. It doesn't make any change.

pattern Plain 5 "Sine 0 1000 $period" {
sp 3 1 [expr $ysize*0.00316*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 5
loadConst

pattern Plain 6 "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 6
loadConst

pattern Plain 7 "Sine 0 1000 $period" {
sp 3 1 [expr $ysize*0.0316*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.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 8
loadConst

pattern Plain 9 "Sine 0 1000 $period" {
sp 3 1 [expr $ysize*0.316*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 2000 $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 10
loadConst

pattern Plain 11 "Sine 0 1000 $period" {
sp 3 1 [expr $ysize*1.78*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.16*0.01]
}
puts "$Analysis_case"
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 output stream

[skamalzare]

Hi,

That's my question too! Were you able to figure out how we can update RefBulkModulus?

Thanks,
Soheil

[Rain]

Hi,

I have some ideas about your second question. The parameter-updating command have already been modified. It can only update the permeability parameters of some solid-fluid fully coupled elements, like the quadUP element in your first example. It's nothing to do with the reference bulk modulus of PDMY02 material. However, you can use "updateMaterials" command to update the modulus in materials including PDMY and PIMY. I don't know whether it can be useful to PDMY02 material, but you can have a try. The command is as follows:

updateMaterials -material $tag shearModulus $newVal
updateMaterials -material $tag bulkModulus $newVal

[skamalzare]

Thank you Rain, I used your suggestions and it worked.

Bests,
SK

[mirkomorosi]

Thank you Rain,
I had already realized how to update soil bulk modulus within my script, implemented in the same way you suggested me.

"...
# soil bulk modulus
set bulkDRY_UPD [expr $EDRY_UPD/(3*(1-2*$nuUPD))]
updateMaterials -material 1 bulkModulus $bulkDRY_UPD;
..."

And I have also realized that this command works fine with PDMY02 material.

Best Regards,
Mirko