Getting error in simplified static analysis

Forum for OpenSees users to post questions, comments, etc. on the use of the OpenSees interpreter, OpenSees.exe

Moderators: silvia, selimgunay, Moderators

Post Reply
armkam
Posts: 30
Joined: Thu Oct 13, 2016 4:02 pm
Location: The University of Auckland

Getting error in simplified static analysis

Post by armkam » Sun Nov 27, 2016 4:28 pm

Hi,

I wanted to model a soil-structure interaction problem (Lshaped retaining wall) and I modeled a simplified one by just fixing far and base nodes of the soil and without modeling an interface between wall and soil at this stage.
I am sure that nodes numbers and coordinates and elements are correct since I have checked them using with these two commands "print nodes.out -node" and "print ele.out -ele"
I also changed my ndMaterial to an elastic one but nothing happened.
I also modeled a very simple plate on a 8-element soil and still I get the error (you can find the model code of this at the end).
However I still get these errors

WARNING DOF_GROUP::setID - invalid location 2 in ID of size 2

ProfileSPDLinDirectSolver::solve() - aii < 0 (i,aii): (73,0)
WARNING AcceleratedNEWTON::solveCurrentStep() -the LinearSysOfEqn failed in solve ()
Static Analysis::analyze() - the Algorithm failed at iteration : 0 with domain at load factor 0.1
OpenSees > analyze failed, returned: -3 error flag



Here are the codes for your reference:

*****************************************************************************************************************************************************************************
Lshaped cantilever wall
*****************************************************************************************************************************************************************************

wipe

#-------------------------------------------------------------------------------------------
# 1. DEFINE ANALYSIS PARAMETERS
#-------------------------------------------------------------------------------------------

#---GEOMETRY
# thickness of soil profile
set soilThickBot 20.0
# length of soil profile base
set soilLengthBot 40.0
# wall dimension
set wallDim 3.0
# wall thickness
set wallThick 0.5
# length of soil profile behind the wall
set soilLengthBehWall [expr ($soilLengthBot/2)-($wallDim/3)]

#---MATERIAL PROPERTIES
# soil mass density (Mg/m^3)
set rho 1.7
# soil shear velocity
set Vs 250.0
# soil shear modulus
set Gmax [expr $rho*$Vs*$Vs]
# poisson's ratio of soil
set nu 0.2
# soil elastic modulus
set Es [expr 2*$Gmax*(1+$nu)]
# soil bulk modulus
set bulk [expr $Es/(3*(1-2*$nu))]
# soil friction angle
set phi 35.0
# soil cohesion
set coh 0.0
# peak shear strain
set gammaPeak 0.1
# reference pressure
set refPress 80.0
# pressure dependency coefficient
set pressDependCoe 0.5
# phase transformation angle
set PTAng 27
# rate of shear-induced volume decrease (contraction) or pore pressure buildup
set contrac 0.07
# the rate of shear-induced volume increase (dilation)
set dilat1 0.4
set dilat2 2.0
# parameters controlling the mechanism of cyclic mobility (Kpa)
set liquefac1 10
set liquefac2 0.01
set liquefac3 1

# soil-wall friction angle
set phiFric [expr $phi/2]
set pi [expr acos(-1.0)]

# bedrock shear wave velocity
set rockVs 760
# bedrock mass density (Mg/m^3)
set rockDen 2.4

# retaining wall thickness
set Wthick 0.5
# spacing of walls
set dSpace 1.0
# wall area section
set Warea [expr $Wthick*$dSpace]
# inertia about local z-axis of concrete wall
set Izc [expr $dSpace*pow($Wthick,3.0)/3]

# compression strength of concrete wall (Mpa)
set fc -25.0
# elastic modulus of concrete wall (Kpa)
set Ec [expr 4700*pow(-$fc,0.5)*1000]
# max strain
set epsc0 -0.002
# crushing strength
set fpcu 0.0
# strain at cruching
set epsU -0.005
# steel yield strength (Mpa)
set Fy 300
# steel young modulus (Kpa)
set E0 200000000
# steel strain-hardening ratio
set bSteel 0.01

#---WAVELENGTH PARAMETERS
# highest frequency desired to be well resolved (Hz)
set fMax 100.0
# wavelength of highest resolved frequency
set wave [expr $Vs/$fMax]
# number of elements per one wavelengths
set nWaveEle 10

#-------------------------------------------------------------------------------------------
# 2. DEFINE MESH GEOMETRY
#-------------------------------------------------------------------------------------------

# trial vertical element size
set hTrial [expr $wave/$nWaveEle]
# trial number of elements
set nTrialVer [expr $soilThickBot/$hTrial]
# round up if not an integer number of elements
if { $nTrialVer > [expr floor($soilThickBot/$hTrial)] } {
set [expr int(floor($soilThickBot/$hTrial)+1)]
} else {
set numEleY [expr int($nTrialVer)]
}
puts "vertical number of elements: $numEleY"

# vertical size of elements
set sizeEleY [expr $soilThickBot/$numEleY]
puts "vertical size of elements: $sizeEleY"


# trial horizental element size
set lTrial [expr $wave/$nWaveEle]
# trial number of elements
set nTrialHor [expr $soilLengthBot/$lTrial]
# round up if not an integer number of elements
if { $nTrialHor > [expr floor($soilLengthBot/$lTrial)] } {
set numEleX [expr int(floor($soilLengthBot/$lTrial)+1)]
} else {
set numEleX [expr int($nTrialHor)]
}
puts "horizental number size of elements: $numEleX"

# horizental size of elements
set sizeEleX [expr $soilLengthBot/$numEleX]
puts "horizental size of elements: $sizeEleX"



# trial vertical element size behind the wall
set hWallTrial [expr $wave/$nWaveEle]
# trial number of elements
set nTrialWallVer [expr $wallDim/$hWallTrial]
# round up if not an integer number of elements
if { $nTrialWallVer > [expr floor($wallDim/$hWallTrial)] } {
set numEleWallY [expr int(floor($wallDim/$hWallTrial)+1)]
} else {
set numEleWallY [expr int($nTrialWallVer)]
}
puts "vertical number of elements behind the wall: $numEleWallY"

# final vertical size of elements behind the wall
set sizeEleWallY [expr $wallDim/$numEleWallY]
puts "vertical size of elements behind the wall: $sizeEleWallY"



# trial horizental element size behind the wall
set lWalltrial [expr $wave/$nWaveEle]
# trial number of elements
set nWallTrialHor [expr $soilLengthBehWall/$lWalltrial]
# round up if not an integer number of elements
if { $nWallTrialHor > [expr floor($soilLengthBehWall/$lWalltrial)] } {
set numEleWallX [expr int(floor($soilLengthBehWall/$lWalltrial)+1)]
} else {
set numEleWallX [expr int($nWallTrialHor)]
}
puts "horizental number size of elements behind the wall: $numEleWallX"

# horizental size of elements behind the wall
set sizeEleWallX [expr $soilLengthBehWall/$numEleWallX]
puts "horizental size of elements behind the wall: $sizeEleWallX"


# number of nodes underlying the wall
set NumNodeUndWall [expr ($numEleX+1)*($numEleY+1)]
puts "node numbers underlying the soil as should be: $NumNodeUndWall"

# number of nodes behind the wall
set NumNodeBehWall [expr ($numEleWallX+1)*($numEleWallY+1-1)]
puts "total node numbers as should be: [expr $NumNodeBehWall+$NumNodeUndWall]"


# trial horizental element size of the wall
# round up if not an integer number of elements
set numWallHor [expr $wallDim/$sizeEleWallX]
if { $numWallHor > [expr floor($wallDim/$sizeEleWallX)] } {
set numWallHor [expr int(floor($wallDim/$sizeEleWallX)+1)]
} else {
set numWallHor [expr int($numWallHor)]
}
puts "horizental number of elements of the wall: $numWallHor"

# trial vertical element size of the wall
# round up if not an integer number of elements
set numWallVer [expr $wallDim/$sizeEleWallY]
if { $numWallVer > [expr floor($wallDim/$sizeEleWallY)] } {
set numWallVer [expr int(floor($wallDim/$sizeEleWallY)+1)]
} else {
set numWallVer [expr int($numWallHor)]
}
puts "vertical number of elements of the wall: $numWallVer"

#-------------------------------------------------------------------------------------------
# 3. DEFINE NODES FOR SOIL ELEMENTS
#-------------------------------------------------------------------------------------------

model BasicBuilder -ndm 2 -ndf 2


# loop over nodes underlying the wall
set h 0.0
for {set k 1} {$k <= [expr $numEleY+1]} {incr k 1} {
for {set j 1} {$j <= [expr $numEleX+1]} {incr j 1} {
node [expr int($h+$j)] [expr ($j-1)*$sizeEleX] [expr ($k-1)*$sizeEleY]
}

set h [expr $h+$j-1]
}
puts "Finish creating soil nodes underlying the wall..."


# loop over nodes behind the wall
set hh $NumNodeUndWall
set xBegin [expr $soilLengthBot-$soilLengthBehWall]
set yBegin $soilThickBot

for {set kk 1} {$kk <= [expr $numEleWallY]} {incr kk 1} {
for {set jj 1} {$jj <= [expr $numEleWallX+1]} {incr jj 1} {
if { $jj == [expr 1] & $hh == $NumNodeUndWall } {
puts "Are node numbers underlying the soil ($NumNodeUndWall) ok? [expr int($hh+$jj)-1]"
}
node [expr int($hh+$jj)] [expr $xBegin+(($jj-1)*$sizeEleX)] [expr $yBegin+(($kk)*$sizeEleY)]
}

set hh [expr $hh+$jj-1]
}
puts "Are total node numbers soil ([expr $NumNodeBehWall+$NumNodeUndWall]) ok? [expr int($hh)]"
puts "Finish creating soil nodes underlying the wall..."

set nodeBehWallStart [expr int(($numEleX+1)*(($numEleY-1)+1)+(($soilLengthBot-$soilLengthBehWall)/$sizeEleX)+1)]


#######################################################################################

#-------------------------------------------------------------------------------------------
# 4. DEFINE DASHPOT NODES
#-------------------------------------------------------------------------------------------



# FIX SOIL FAR ENDS UNDER THE WALL
for {set j [expr $numEleY+1]} {$j <= [expr $numEleY+1]} {incr j 1} {
fix [expr int(1+($j-1)*($numEleX+1))] 1 1 1
fix [expr int($j*($numEleX+1))] 1 1 1
}


# FIX SOIL FAR END BEHIND THE WALL
for {set j 1} {$j <= [expr $numEleWallY]} {incr j 1} {
fix [expr int(($numEleX+1)*($numEleY+1)+($j)*($numEleWallX+1))] 1 1 1
}


# FIX SOIL BASE NODES
for {set j 1} {$j <= [expr $numEleX+1]} {incr j 1} {
fix $j 1 1 1
}

puts "Fix Far Ends of Soil"


# SOIL MATERIAL
#nDMaterial PressureDependMultiYield $tag $nd $rho $refShearModul $refBulkModul $frictionAng $peakShearStra $refPress $pressDependCoe $PTAng $contrac $dilat1 $dilat2 $liquefac1 $liquefac2 $liquefac3 <$noYieldSurf=20 <$r1 $Gs1 ?> $e=0.6 $cs1=0.9 $cs2=0.02 $cs3=0.7 $pa=101>
nDMaterial PressureDependMultiYield 100 2 $rho $Gmax $bulk $phi $gammaPeak $refPress $pressDependCoe $PTAng $contrac $dilat1 $dilat2 $liquefac1 $liquefac2 $liquefac3
#nDMaterial ElasticIsotropic $matTag $E $v <$rho>
#nDMaterial ElasticIsotropic 1 $Es $nu $rho
nDMaterial InitialStateAnalysisWrapper 1 100 2


################################################################################

# soil elements
set wgtX 0.0
set wgtY [expr -9.81*$rho]

for {set k 1} {$k <= $numEleY} {incr k 1} {
for {set j 1} {$j <= $numEleX} {incr j 1} {

set nI [expr $j+($k-1)*($numEleX+1)]
set nJ [expr ($j+1)+($k-1)*($numEleX+1)]
set nK [expr ($j+1)+($k*($numEleX+1))]
set nL [expr $j+($k*($numEleX+1))]

element quad [expr int($j+($k-1)*$numEleX)] $nI $nJ $nK $nL 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY
}
}


for {set kk 1} {$kk <= $numEleWallY} {incr kk 1} {
for {set jj 1} {$jj <= $numEleWallX} {incr jj 1} {

set nI [expr $nodeBehWallStart+($jj-1)+($kk-1)*($numEleWallX+1)]
set nJ [expr $nodeBehWallStart+$jj+($kk-1)*($numEleWallX+1)]
set nK [expr $nodeBehWallStart+$jj+($kk*($numEleWallX+1))]
set nL [expr $nodeBehWallStart+($jj-1)+($kk*($numEleWallX+1))]

element quad [expr int(($numEleX*$numEleY)+$jj+($kk-1)*$numEleWallX)] $nI $nJ $nK $nL 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY
}
}
puts "Total number of elements is ([expr ($numEleX*$numEleY)+($numEleWallX*$numEleWallY)]) is the loop ok? [expr int(($numEleX*$numEleY)+$jj-1+($kk-2)*$numEleWallX)]"
puts "Finished creating all soil elements..."


#-------------------------------------------------------------------------------------------
# 5. DEFINE WALL AND LAGRANGIAN NODES
#-------------------------------------------------------------------------------------------

model BasicBuilder -ndm 2 -ndf 3



for {set tt 20002} {$tt <= [expr 20002+$numWallHor]} {incr tt 1} {
node $tt [expr $soilLengthBot-$soilLengthBehWall-$wallDim+($tt-20002)*$sizeEleWallX] $soilThickBot
}
puts " tt out loop $tt"


# vertical part of Wall
for {set pp $tt} {$pp <= [expr $tt+$numWallVer-1]} {incr pp 1} {
node $pp [expr $soilLengthBot-$soilLengthBehWall] [expr $soilThickBot+(($pp+1)-$tt)*$sizeEleWallY]
puts " pp $pp"
}
puts "pp out loop $pp"

#-------------------------------------------------------------------------------------------
# 7. DEFINE MATERIALS
#-------------------------------------------------------------------------------------------


# Wall Material
#uniaxialMaterial Concrete01 2 $fc $epsc0 $fpcu $epsU
#uniaxialMaterial Steel01 3 $Fy $E0 $bSteel
section Elastic 1 $E0 $Warea $Izc



#-------------------------------------------------------------------------------------------
# 8. DEFINE Elements
#-------------------------------------------------------------------------------------------



# wall elements

set numIntPts 3
set transTag 1
# geometric transformation
geomTransf Linear $transTag
for {set tt 20002} {$tt <= [expr 20002+$numWallHor-1]} {incr tt 1} {
element dispBeamColumn $tt $tt [expr int($tt+1)] $numIntPts 1 $transTag
puts "tt $tt"
}
puts "tt out loop $tt"



for {set pp $tt} {$pp <= [expr $tt+$numWallVer-1]} {incr pp 1} {
element dispBeamColumn $pp $pp [expr int($pp+1)] $numIntPts 1 $transTag
puts " pp $pp"
}
puts "pp out loop $pp"



#-------------------------------------------------------------------------------------------
# 10. CREATE GRAVITY RECORDERS
#-------------------------------------------------------------------------------------------


## record stress and strain at each gauss point in the soil elements
#recorder Element -file Gstress1.out -time -eleRange 1 $numEleY material 1 stress
#recorder Element -file Gstress2.out -time -eleRange 1 $numEleY material 2 stress
#recorder Element -file Gstress3.out -time -eleRange 1 $numEleY material 3 stress
#recorder Element -file Gstress4.out -time -eleRange 1 $numEleY material 4 stress
#
#recorder Element -file Gstrain1.out -time -eleRange 1 $numEleY material 1 strain
#recorder Element -file Gstrain2.out -time -eleRange 1 $numEleY material 2 strain
#recorder Element -file Gstrain3.out -time -eleRange 1 $numEleY material 3 strain
#recorder Element -file Gstrain4.out -time -eleRange 1 $numEleY material 4 strain

puts "Finished creating gravity recorders..."

print nodes.out -node
print ele.out -ele

#-----------------------------------------------------------------------------------------------------------
# 11. APPLY STATIC LOADING
#-----------------------------------------------------------------------------------------------------------
# load weight of 1m of wall (KN)

pattern Plain 1 Constant {

set Wload [expr -2*1*$wallThick*$wallDim*2400*9.81/1000]
set nLoad [expr 20001+int(2*$wallDim/(3*$sizeEleX))]
load $nLoad 0 $Wload 0

set WloadX -100
set nloadX [expr 20001+int($wallDim/$sizeEleX)+int($wallDim/$sizeEleY)-2]
load $nloadX $Wload 0 0

}


#-----------------------------------------------------------------------------------------------------------
# 12. STATIC ANALYSIS
#-----------------------------------------------------------------------------------------------------------


constraints Transformation
test NormDispIncr 1e-4 60 1
algorithm KrylovNewton
numberer RCM
system ProfileSPD
integrator LoadControl 0.1
analysis Static
analyze 10

puts "Finished with elastic gravity analysis..."














*****************************************************************************************************************************************************************************
Simple Plate on Elastic Soil
*****************************************************************************************************************************************************************************

wipe
model BasicBuilder -ndm 2 -ndf 2

node 1 0 0
node 2 0.5 0
node 3 1.0 0
node 4 1.5 0
node 5 2.0 0
node 6 0.0 0.5
node 7 0.5 0.5
node 8 1.0 0.5
node 9 1.5 0.5
node 10 2.0 0.5
node 11 0.0 1.0
node 12 0.5 1.0
node 13 1.0 1.0
node 14 1.5 1.0
node 15 2.0 1.0


fix 1 1 1 1
fix 2 1 1 1
fix 3 1 1 1
fix 4 1 1 1
fix 5 1 1 1
fix 6 1 1 1
fix 10 1 1 1
fix 11 1 1 1
fix 15 1 1 1


nDMaterial ElasticIsotropic 1 177083 0.2 1.7

set wgtX 0.0
set wgtY [expr -9.81*1.7]
element quad 1 1 2 7 6 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY
element quad 2 2 3 8 7 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY
element quad 3 3 4 9 8 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY
element quad 4 4 5 10 9 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY
element quad 5 6 7 12 11 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY
element quad 6 7 8 13 12 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY
element quad 7 8 9 14 13 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY
element quad 8 9 10 15 14 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY


model BasicBuilder -ndm 2 -ndf 3

node 1012 0.5 1.0
node 1013 1.0 1.0
node 1014 1.5 1.0


section Elastic 1 200000000 0.5 0.0417
geomTransf Linear 1

element dispBeamColumn 106 1012 1013 3 1 1
element dispBeamColumn 107 1013 1014 3 1 1

print nodes.out -node
print ele.out -ele

pattern Plain 1 Constant {
set Wload [expr -2*1*0.5*1.0*2400*9.81/1000]
load 1013 0 $Wload 0
set WloadX -100
load 1012 $WloadX 0 0
}

constraints Transformation
test NormDispIncr 1e-4 60 1
algorithm KrylovNewton
numberer RCM
system ProfileSPD
integrator LoadControl 0.1
analysis Static
analyze 10

Post Reply