Dear Opensees Users
Below is the model I developed to analyze the nonlinear behavior of pier type structure and soil. I am utilizing the Contact Interface Model and an RC column for this simulation. The model works as long as the seismic excitation keeps the structure with in elastic range (I have tested with 0.1g level for the same exitation) but the moment the it exceeds the moment capacity (40Mpa) I get convergence problems. I am a relatively new user to opensees and I do not know any experienced users to turn for advise. I would be quite grateful if someone can direct me to get this model to run properly. I have not included the takatori excitation to this post (it's really big) but should you need it I would be happy to send it.
Thank you very much for your time!
################## FUKAE TAKATORI CIM ROCKING #######################
###### UNITS USED: NEWTON, METER, KG & SEC
wipe
#------build a 2D model with 3 DOF-----#
model BasicBuilder -ndm 2 -ndf 3
#----structure properties-----------#
set Efoot 7.0e10
set Afoot 49
set Ifoot 200
set Hfoot 1
set Hfootcg 0.75
set Hcol 12
set Mdeck 12.0e5
#----define parameters for macro element use---#
set Rv 0.1
set deltaL 0.01
set L 7
set B $L
set Vtot [expr $Mdeck*9.81]
set Df 1
#----soil parameters unit_weight angle_of_friction & cohesion----#
set gamma 16000
set phi 40
set c 0
set pi [expr acos(-1)]
source FScalc.tcl
set Vult [expr $Vtot*$FS]
#------CIM parameter calcuation-------#
set nu 0.3
set N160 25
source embeddedstiffness.tcl
#----define nodes------#
node 1 0 0
node 2 0 0
node 3 0 $Hfoot
node 4 0 $Hfootcg
node 5 0 [expr $Hfoot+$Hcol]
#---geometry of footing/ column elements----#
geomTransf Linear 1
#------foundation/soil element-------#
#section SFS2d matID FS Vult L Kv Kh Rv deltaL
section soilFootingSection2d 1 $FS $Vult $L $Kv $Kh $Rv $deltaL
#element ZLS eleID ndi ndj matID
element zeroLengthSection 1 1 2 1 -orient 0 -1 0 1 0 0
#------foundation element------#
#footing element eleID ndi ndj A E I Transtag
element elasticBeamColumn 2 2 3 $Afoot $Efoot $Ifoot 1
#--------------define material for nonlinear column-----------------#
#confined concrete mat.Tag fc1C eps1C fc2C eps2C lambda ftC Ets
uniaxialMaterial Concrete02 1 -56e6 -0.003 -44.8e6 -0.02 0.1 5.6e6 3.2e9
#unconfined concrete mat.Tag fc1U eps1U fc2U eps2U lambda ftU Ets
uniaxialMaterial Concrete02 2 -44.8e6 -0.003 -4.5e6 -0.006 0.1 4.5e6 3.2e9
#uniaxialMaterial Hysteretic $IDHysteretic $Fy $epsY $Fy1 $epsY1 $Fu $epsU -$Fy -$epsY -$Fy1 -$epsY1 -$Fu -$epsU $pinchX $pinchY $damage1 $damage2 $betaMUsteel
uniaxialMaterial Hysteretic 3 620e6 0.003 690e6 0.04 655e6 0.25 -620e6 -0.003 -690e6 -0.04 -620e6 -0.25 1 1 0 0.2 0.5
#--------------define column cross-section dimensions----------------#
#column diameter
set D 3
#number of reinforments
set numBar 100
#cover concrete
set Cover 0.0381
#area of steel for each longitudinal #10rebars
set As 0.000819
#-------creating fiber sections for core/ cover concrete and reinforcements-------#
set ri 0
set r0 [expr $D/2]
set nfCoreR 8
set nfCoreT 8
set nfCoverR 4
set nfCoverT 8
set rc [expr $r0-$Cover]
section Fiber 2 {
#create core fiber mat.Tag numSubDivCirc numSubDivRad yCenter zCenter intRad extRad startAng endAng
patch circ 1 $nfCoreT $nfCoreR 0 0 $ri $rc 0 360
#create cover fiber mat.Tag numSubDivCirc numSubDivRad yCenter zCenter intRad extRad startAng endAng
patch circ 2 $nfCoverT $nfCoverR 0 0 $rc $r0 0 360
#create reinforcement fibers mat.Tag $numFiber $areaFiber $yCenter $zCenter $radius <$startAng $endAng>
set theta [expr 360.0/$numBar]
layer circ 3 $numBar $As 0 0 $rc $theta 360
}
#-------------------define column elements----------------------#
#number of integration points along length of element
set npbase 2
set npcolumn 10
set eleType forceBeamColumn
#create column using beam-column elements tag ndI ndJ nsecs secID transferTag
element $eleType 3 3 4 $npbase 2 1
element $eleType 4 4 5 $npcolumn 2 1
#-----------boundary condition------------------------#
#fix base node in all directions nodetag x y rotation
fix 1 1 1 1
#-------------gravity analysis------------------------#
#define gravity loads - in 10 increments - deck mass 1200Mg
pattern Plain 1 Linear {
load 5 0 -1177200 0
}
#define analysis - gravity load
test NormDispIncr 1e-8 10 1
algorithm Newton
system UmfPack
constraints Plain
numberer Plain
analysis Static
#-------define recorders--------#
set name "node"
for {set n 1} {$n <= 5} {incr n 1} {
set fileName [join [list $name $n] {} ]
recorder Node -file $fileName -node $n -dof 1 2 3 disp
}
set name "element"
for {set n 1} {$n <= 4} {incr n 1} {
set fileName [join [list $name $n] {} ]
recorder Element -file $fileName -time -ele $n force
}
#apply gravity loads in 10 steps
analyze 10
#setting time back to 0 before dynamic analysis
loadConst -time 0.0
#indicate factor of safety
puts "factor of safety $FS"
puts "-----Gravity Analysis Done-----"
wipeAnalysis
#----------stuperstructure node--------------#
#define mass at node 5 in direction x y rotation
mass 5 $Mdeck $Mdeck 1e-3
#-----------dynamic analysis------------#
test EnergyIncr 1e-5 100 1
algorithm Newton
system UmfPack
constraints Plain
numberer RCM
integrator Newmark 0.6 0.32
analysis VariableTransient
#-----------apply damping---------#
set xDamp 0.06
set lambda [eigen 1]
set omega [expr pow($lambda,0.5)]
set alphaM 0.05
set betaK 0
set betaKcomm [expr 2*$xDamp/$omega]
set betaKinit 0
rayleigh $alphaM $betaK $betaKinit $betaKcomm
#-----------define ground motion----------#
set dT 0.0001
set dTmin [expr $dT/1000]
set dTmax $dT
# acceleration time history is read from an external file
set Series "Path -filePath Takatorinew_input.txt -dt $dT -factor 4.9"
# acceleration is applied at the fixed node in horizontal direction (1)
pattern UniformExcitation 2 1 -accel $Series
#apply shake
set steps 409500
set itr 100
for {set i 1} {$i < $steps} {incr i 1} {
test EnergyIncr 1e-4 $itr 0
set ok [analyze 1 $dT $dTmin $dTmax $itr]
if {$ok != 0} {
test EnergyIncr 1e-3 $itr 0
set ok [analyze 1 $dT $dTmin $dTmax $itr]
}
if {$ok != 0} {
test EnergyIncr 1e-2 $itr 0
set ok [analyze 1 $dT $dTmin $dTmax $itr]
}
if {$ok != 0} {
test EnergyIncr 1e-1 $itr 0
set ok [analyze 1 $dT $dTmin $dTmax $itr]
}
}
# print out final node and element outputs on screen
for {set n 1} {$n <= 5} {incr n 1} {
print node $n
}
print ele
puts "done"
#end
wipe
############ FOOTING STIFFNESS CALCULATION SUBROUTINE ###########
########## UNITS USED: NEWTON, METER, KG & SEC ##################
######### FOOTING EMBEDMENT Df ############
#---------shear modulus (FEMA 356)----------#
set sv [expr $Vtot/($B*$L)]
set Go [expr 20000*pow($N160,0.33)*sqrt($sv)]
set G_overGo 1.0
set G [expr $Go*$G_overGo]
set G [expr $G*47.88]
puts "G $G"
#----------stiffness at surface----------#
set Kxsurf [expr $G*$B/(2-$nu)*(3.4*pow(($L/$B),0.65)+1.2)]
set Kzsurf [expr $G*$B/(1-$nu)*(1.55*pow(($L/$B),0.75)+0.8 )]
set d $Hfoot
#embedment depth center
set h [expr ($Df-$d)+$d/2]
set BetaX [expr (1+0.21*sqrt($Df/$B))*(1+1.6*pow(($h*$d*($B+$L)/($B*$L*$L)),0.4))]
set BetaZ [expr (1+($Df/(21*$B))*(2+2.6*$B/$L))*(1+0.32*pow(($d*($B+$L)/($B*$L)),(2/3)))]
#-------------embedded stiffness-----------#
set Kh [expr $BetaX*$Kxsurf]
set Kv [expr $BetaZ*$Kzsurf]
puts "Kh $Kh"
puts "Kv $Kv"
########### FS CALCULATION SUBROUTINE ##############
######## UNITS USED: NEWTON, METER, KG & SEC #######
#----bearing capacity modification factors-----#
set Phi [expr $phi*$pi/180]
set Nq [expr pow(tan($phi/2+$pi/4),2)*exp($pi*tan($Phi))]
set Nc [expr ($Nq-1)/tan($Phi)]
set Ng [expr 2*($Nq+1)*tan($Phi)]
set Fqs [expr 1+($B/$L)*tan($Phi)]
set Fcs [expr 1+($B/$L)*($Nq/$Nc)]
set Fgs [expr 1-0.4*($B/$L)]
set Fcd [expr 1+0.4*($Df/$B)]
set Fqd [expr 1+2*tan($Phi)*($Df/$B)*pow((1-sin($Phi)),2)]
set Fgd 1
#-----bearing capacity calculation-----#
set qu [expr $c*$Nc*$Fcs*$Fcd+$gamma*$Df*$Nq*$Fqs*$Fqd+0.5*$gamma*$B*$Ng*$Fgs*$Fgd]
set q [expr $Vtot/($B*$L)]
set FS [expr $qu/$q]
Need Help with Convergence Analysis
Moderators: silvia, selimgunay, Moderators