Error in Soil and Wall Interface

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

Error in Soil and Wall Interface

Post by armkam »

Hi all,
I have modeled a Lshaped retaining wall and dry sand and wanted to define interface by BeamContact2D but I've gotten these kind errors that you can find an example following:

"element id: 30001
Connected external nodes: 20002 20003 12953 30002 BeamContact2D"
Node 20002 (17.785,20.125) and 20003 (18.125,20.125) are nodes of the wall (3dof), 12953 (18,20) is the soil node (2dof) which is located right under and between of nodes 20002 and 20003 and 30002 is lagrangian node with (0,0) coordinates (2dof).

I am not sure what this error means.
if it is applicable the model is copied below:



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
#-------------------------------------------------------------------------------------------

node 20000 0.0 0.0
node 20001 0.0 0.0

puts "Finished creating dashpot nodes..."

#################################################################
#-------------------------------------------------------------------------------------------
# 6. DEFINE BOUNDARY CONDITIONS AND EQUAL DOF
#-------------------------------------------------------------------------------------------




# define fixity of base nodes
for {set j 1} {$j <= [expr $numEleX+1]} {incr j 1} {
fix $j 0 1
}


# define fixity of dashpot nodes
fix 20000 1 1
fix 20001 0 1

# define equal DOF for simple shear deformation of soil elements
for {set k [expr $numEleX+1]} {$k <= [expr $NumNodeUndWall-$numEleX]} {incr k [expr $numEleX]} {
for {set h $k} {$h <= [expr $k+$numEleX-1]} {incr h 1} {
equalDOF $k [expr $h+1] 1 2
}
}

for {set kk [expr $NumNodeUndWall+1]} {$kk <= [expr $NumNodeBehWall+$NumNodeUndWall-$numEleWallX]} {incr kk [expr $numEleWallX+1]} {
for {set hh $kk} {$hh <= [expr $kk+$numEleWallX-1]} {incr hh 1} {
equalDOF $kk [expr $hh+1] 1 2
}
}



# define equal DOF for dashpot and base soil nodes
for {set j 2} {$j <= [expr $numEleX+1]} {incr j 1} {
equalDOF 1 [expr int($j)] 1
}
equalDOF 1 20001 1

puts "Finished creating all boundary conditions and equalDOF..."

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

# 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 1 2 $rho $Gmax $bulk $phi $gammaPeak $refPress $pressDependCoe $PTAng $contrac $dilat1 $dilat2 $liquefac1 $liquefac2 $liquefac3



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

# 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 "Toal 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..."

# lagrangian nodes

for {set ff 30002} {$ff <= [expr 30002+$numWallHor+$numWallVer-1]} {incr ff 1} {
node $ff 20.875 20.125
puts "ff $ff"
}

puts "ff out loop $ff"

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

model BasicBuilder -ndm 2 -ndf 3


# horizental part of Wall
for {set tt 20002} {$tt <= [expr 20002+$numWallHor]} {incr tt 1} {
node $tt [expr $soilLengthBot-$soilLengthBehWall-$wallDim+($tt-20002)*$sizeEleWallX-$sizeEleWallX/2] [expr $soilThickBot+$sizeEleWallY/2]
puts " tt $tt"
}
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-$sizeEleWallX/2] [expr $soilThickBot+(($pp+1)-$tt)*$sizeEleWallY+$sizeEleWallY/2]
puts " pp $pp"
}
puts "pp out loop $pp"


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


# Wall Material

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"


# interface material
nDMaterial ContactMaterial2D 2 [expr tan($phiFric*$pi/10)] 1000.0 0.0 0.0



puts "Finished creating all materials..."

# set gap and force tolerances for beam contact elements
set gTol 0
set fTol 0

# interface elements
for {set tt 1} {$tt <= $numWallHor} {incr tt 1} {
set nConSoil [expr int(($numEleX+1)*(($numEleY-1)+1)+(($soilLengthBot-$soilLengthBehWall-$wallDim)/$sizeEleX)+$tt)]
set nINode [expr 20001+$tt]
set nJNode [expr 20001+($tt+1)]
set nLag [expr 30001+$tt]
element BeamContact2D [expr int(30000+$tt)] $nINode $nJNode $nConSoil $nLag 2 $sizeEleX $gTol $fTol
puts "tt $tt"
}
puts "tt out loop $tt"


#for {set pp $tt} {$pp <= [expr $numWallHor+$numWallVer]} {incr pp 1} {
set nConSoil [expr int(($numEleX+1)*($numEleY+1)+($pp-$tt)*($numEleWallX+1)+1)]
set nINode [expr 20001+$pp]
set nJNode [expr 20001+($pp+1)]
set nLag [expr 30001+$pp]
element BeamContact2D [expr int(30000+$pp)] $nINode $nJNode $nConSoil $nLag 2 $sizeEleX $gTol $fTol
}


#-------------------------------------------------------------------------------------------
# 9. DEFINE MATERIAL AND ELEMENTS FOR VISCOUS DAMPERS
#-------------------------------------------------------------------------------------------


# dashpot coefficient
set mC [expr $sizeEleX*$rockDen*$rockVs]

# material
uniaxialMaterial Viscous 4000 $mC 1

# elements
element zeroLength 50000 20000 20001 -mat 4000 -dir 1

puts "Finished creating dashpot material and element..."

#-------------------------------------------------------------------------------------------
# 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.25
analysis Static
analyze 4

puts "Finished with elastic gravity analysis..."
YuriWong
Posts: 2
Joined: Tue Sep 27, 2016 6:26 pm
Location: University of Auckland

Re: Error in Soil and Wall Interface

Post by YuriWong »

Hi armkam,

Arman? Hahaha.

I think that "BeamContact2D, element id: 30001 Connected external nodes: 20002 20003 12953 30002" is just the output of "print ele.out -ele". The error is the "WARNING" part following this output.

However, I am not sure what is wrong with your codes. You may need to check the boundary conditions and add a gravity analysis before the "STATIC LOADING", or try a transient analysis.

Good luck.

Yuri
Post Reply