Hello all,
I am trying to model VELACS 1, it is one of my first attempts to use OpenSees so I have a little bit of problems.
I have modeled the geometry of the problem, and in gravity loading (soil self weight) I have faced a problem and get errors while running the code.
I have distributed the soil weight forces to the nodes.
I hope somebody can help me with it,
here's the code:
# Units: m, sec, ton, kN
wipe all
set startT [clock seconds]
################################### Soil Properties ##########################################
set Gs 2.68 ;# Soil Grain mass density (Gs)
set nou 0.05 ;# Poisson's Ratio
set rhof 1.0 ;# Fluid mass density
set Phi 31.4 ;# Friction angle at peak shear strength, in degree
set Dr 40.0 ;# Relative Density
set e0 [expr 0.894-(0.894-0.516)*$Dr/100]
set poro0 [expr $e0/(1+$e0)]
set rhos [expr ((1.0-$poro0)*$Gs+$poro0)*$rhof] ;# Saturated Soil mass density
# Permeability (m/s)
set kx 2.1e-5
set ky 2.1e-5
set kz 2.1e-5
set gravX 0.0000
set gravY 0.0000
set gravZ -9.8100
################################### Define Number of Analysis steps ###########################
set NumIncr1 50 ; # Number of Steps in Static Analysis
################################### Define Soil Nodes ##########################################
model basic -ndm 2 -ndf 3
set xincr 2.54
set yincr 2.00
for {set i 1} {$i<=[expr int(22.86/$xincr)+1]} {incr i 1} {
for {set j 1} {$j<=[expr int(10.00/$yincr)+1]} {incr j 1} {
set xdim [expr ($i-1)*$xincr ]
set ydim [expr ($j-1)*$yincr ]
set NodeNum [expr $j+ (int(10.00/$yincr)+1)*($i-1)]
node $NodeNum $xdim $ydim
puts "$NodeNum $xdim $ydim"
}
}
########################## Define Constitutive Model for soil (Dafalias-Manzari) ###############
set Gr 7.5e4
set Br 2.0e5
set PeakShearStrain 0.1
set Pr 80
set pressDependCoe 0.5
set PTAng 27
set contrac 0.07
set dilat1 0.4
set dilat2 2
set liquefac1 10
set liquefac2 0.01
set liquefac3 1
nDMaterial PressureDependMultiYield 1 2 $rhos $Gr $Br $Phi $PeakShearStrain $Pr $pressDependCoe $PTAng $contrac $dilat1 $dilat2 $liquefac1 $liquefac2 $liquefac3
################################### Define Soil Elements #######################################
set m 1
for {set i 1} {$i <= [expr int(22.86/$xincr)]} {incr i 1} {
for {set j 1} {$j <= [expr int(10.00/$yincr)]} {incr j 1} {
set n1 [expr 6*$i-(6-$j)]
set n2 [expr 6*$i-(6-$j)+1]
set n3 [expr 6*$i+$j]
set n4 [expr 6*$i+$j+1]
set EleNum [expr ($m)]
element quadUP $EleNum $n1 $n2 $n3 $n4 1.0 1 [expr 2.2e6/$poro0] $rhof $kx $ky 0 -9.81 0
set m [expr ($m+1)]
puts "$EleNum $n1 $n2 $n3 $n4 $j"
}
}
### update material stage to zero for elastic behavior
updateMaterialStage -material 1 -stage 0
################################### Boundary Conditions #######################################
### Lateral Boundaries
fixX 0.0 1 0 0 -tol 1e-10
fixY 0.0 1 1 0 -tol 1e-10
fixX 22.86 1 0 0 -tol 1e-10
for {set m 1} {$m<=5} {incr m 1} {
equalDOF $m [expr $m+1] 2
equalDOF [expr $m+54] [expr $m+55] 2
}
### Surface Pore Pressure
for {set m 1} {$m<=10} {incr m 1} {
fix [expr 6*$m] 0 0 1
}
########################### Applying Gravity Loads to Nodes ###########################
set m 1
for {set i 1} {$i <= [expr int(22.86/$xincr)]} {incr i 1} {
for {set j 1} {$j <= [expr int(10.00/$yincr)]} {incr j 1} {
set n1 [expr 6*$i-(6-$j)]
set n2 [expr 6*$i-(6-$j)+1]
set n3 [expr 6*$i+$j]
set n4 [expr 6*$i+$j+1]
pattern Plain $m Linear {
load $n1 0 [expr $xincr*$yincr*9.81*($rhos-$rhof)/4] 0
load $n2 0 [expr $xincr*$yincr*9.81*($rhos-$rhof)/4] 0
load $n3 0 [expr $xincr*$yincr*9.81*($rhos-$rhof)/4] 0
load $n4 0 [expr $xincr*$yincr*9.81*($rhos-$rhof)/4] 0
set m [expr ($m+1)]
}
}
}
################################# Gravity Analysis #################################
system UmfPack
constraints Transformation
test NormDispIncr 0.01 50 1
algorithm Newton
numberer RCM
integrator LoadControl [expr 1.0/10]
analysis Static
analyze 10
updateMaterialStage -material 1 -stage 1
analyze 40
loadConst -time 0.0
VELACS model-1 Modelling
Moderator: Moderators