Working now, it gives 0.249 seconds
, I do not know if it is ok, I just made it run.
What I did to your code was:
1. changed the builder to "model BasicBuilder -ndm 2 -ndf 3"
2. Based on the node coordinates I assumed [3] -- zero length -- [1] -- FBE -- [2], so I changed your code accordingly
3. The steel02 definition lacked R0 R1 and R2, which are not optional, so I uncommented them and add them to the definition
4. At "section fiberSec $ColSecTag" you had not only fibre commands but "set" commands inside the { } as well. I took them out of the {}
5. The definition of the "layer circ" had at the end this "$theta 360", I think that $theta was wrong. I deleted it as I assumed the reinforcement was evenly distributed around the section
6. When you define the BS section the concrete material should have an Elastic Perfectly Plastic behaviour or define it in a way that do not loose strength under large strains. This means that the Concrete04 that you are using for the fibre section of the pier cannot be used in the definition of the Bond_SP section. I used Concrete01 (because it was simpler) to define 2 new materials with EPP behaviour. They are based on the material properties used for the pier fibre section.
7. The shear behaviour of the zeroLengthSection should be stiff, you had assigned G·A, so I changed it to 1e8 to make it stiff. You can include G·A in the FBE element (with an aggregator) to include the shear deformations happening in the pier. Note: If using DBEs aggregator will not work, so the shear defs should be included with an aggregator to the zeroLengthSection element, using G·A/H. The use of the constants assumes linear elastic shear response.
8. As mentioned in 2. I re-defined the FBE from node 1 to 2.
9. For the ZeroLengthSection from 3 to 1, I changed the vector you had, I used "-orient 0 1 0 1 0 0", so local x points in Global Y and local y points in Global X. This way it matches with the FBE as well. Please note that in 2D the the local Y should be kept in the XY plane, yours was pointing in global Z (0 0 -1).
The error you had pointed out a singular matrix, so some of the elements in your model had to be defined in such a way that does not make sense to OpenSees. This, I think, it was the Vy of the BS section pointing in global Z, which is a dimension that does not exist in a 2D domain.
The revised code below:
# Performance Modeling Strategies for Modern Reinforced Concrete Bridge Columns
# date: 16/12/2018
# Units in kip/in
# column.tcl
# This model is a 2D model
# SET UP ----------------------------------------------------------------------------
model BasicBuilder -ndm 2 -ndf 3; # 2 spacial dimensions, 3 DOF's per nodey
set dataDir DataBS; # set up name for data directory
file mkdir $dataDir/; # create data directory
set GMdir "../GMfiles"; # ground-motion file directory
# -----------------------------------------------------
# define section geometry
set Dcol 24.0; # Column Diameter (in)
set LCol 96.0; # column length (in) # Specimen Lehman et al. N°415. L/D=4.0
set P 147.0 ; # weight (kips)
set Weight [expr $P]
set g 386.22; #[inch/s2]
# calculated parameters
set Mass [expr $P/$g]; # nodal mass
# nodal coordinates:
node 3 0. 0.;
node 1 0. 0.; # node#, X, Y
node 2 0. $LCol
# Single point constraints -- Boundary Conditions
fix 3 1 1 1;
#fix 3 1 0 0;
#equalDOF $rNodeTag $cNodeTag $dof1 $dof2
#equalDOF 1 3 1;
# nodal masses:
mass 2 $Mass 1e-9 0.;
# we need to set up parameters that are particular to the model.
set IDctrlNode 2; # node where displacement is read for displacement control
set IDctrlDOF 1; # degree of freedom of displacement read for displacement control
set iSupportNode "3"; # define support node, if needed.
# Define Column geometry
#-------------------------------------------------------------------
set coverSec 0.874; # Column cover to reinforcing steel (in)
set numBarsSec 22; # number of longitudinal-reinforcement bars in the column section
set dlSec 0.626; # [in] Steel diameter T15.9
set PI 3.14;
set barAreaSec [expr $PI*pow($dlSec,2)/4]; #area of longitudinal-reinforcement bars # (in^2)
# MATERIAL parameters
#-------------------------------------------------------------------
set IDconcC 1; # material ID tag -- confined core concrete
set IDconcU 2; # material ID tag -- unconfined cover concrete
set IDreinf 3; # material ID tag -- reinforcement
# Define uniaxial materials
#-------------------------------------------------------------------
# Mander concrete model for confined and unconfined concrete
# characteristics of materials test of Lehman et al. (PEER 98-01)
# confined concrete (core)
#-------------------------
set fc 4.496; # [Ksi] CONCRETE Compressive Strength
set ec 0.002294; # concrete strain at maximum strength ??????
set ecu 0.008101; # [Ksi] concrete strain at crushing strength ????
set Ec 3822.0; # [Ksi]Concrete Elastic Modulus
set fct 0.5046; #[Ksi] maximum tensile strength of concrete
set et 0.000132; # ultimate tensile strain of concrete
#set beta 0.1; #value defining the exponential curve parameter to define the residual stress
#uniaxialMaterial Concrete04 $matTag $fc $ec $ecu $Ec <$fct $et> <$beta>;
#uniaxialMaterial Concrete04 $IDconcC $fc $ec $ecu $Ec $fct $et;
#uniaxialMaterial Concrete04 $IDconcC [expr -$fc] [expr -$ec] [expr -$ecu] $Ec $fct $et;
uniaxialMaterial Concrete04 $IDconcC -$fc -$ec -$ecu $Ec $fct $et;
set IDconcC2 10
uniaxialMaterial Concrete01 $IDconcC2 -$fc -$ec [expr -$fc*0.99] -$ecu
# unconfined concrete (cover)
#----------------------------
set fcU 4.496; # [Ksi] CONCRETE Compressive Strength
set ecU 0.002; # concrete strain at maximum strength (#I did not find the experimentation values)
set ecuU 0.0035; # [Ksi] concrete strain at crushing strength (#I did not find the experimentation values)
set EcU 3822.0; # [Ksi]Concrete Elastic Modulus
set fctU 0.5046; #[Ksi] maximum tensile strength of concrete
set etU 0.000132; # ultimate tensile strain of concrete
#set betaU 0.1; #value defining the exponential curve parameter to define the residual stress
# uniaxialMaterial Concrete04 $matTag $fc $ec $ecu $Ec <$fct $et> <$beta>;
# niaxialMaterial Concrete04 $IDconcU [expr -$fcU] [expr -$ecU] [expr -$ecuU] $EcU $fctU $etU;
uniaxialMaterial Concrete04 $IDconcU -$fcU -$ecU -$ecuU $EcU $fctU $etU;
set IDconcU2 11
uniaxialMaterial Concrete01 $IDconcU2 -$fcU -$ecU [expr -$fcU*0.99] -$ecuU
# ----------- Reinforcing Steel --------
set Es 29000.0 ; # modulus of steel (Ksi)??????
set Fy 67.0; # [Ksi] yield stress Lehman et al.2000 (specimen 415.)
set Fu 91.37; #[Ksi] ultimate stress Lehman et al.2000
set B 0.01; # strain-hardening ratio (Lehman et al. 2000)
set R0 18.5; # control the transition from elastic to plastic branches
set cR1 0.925; # control the transition from elastic to plastic branches
set cR2 0.15; # control the transition from elastic to plastic branches
#set a1 0.04;
#set a2 1.0;
#set a3 0.04;
#set a4 1.0;
#uniaxialMaterial Steel02 $IDreinf $Fy $Es $B $R0 $cR1 $cR2 $a1 $a2 $a3 $a4;
uniaxialMaterial Steel02 $IDreinf $Fy $Es $B $R0 $cR1 $cR2;
puts "materials defined"
#===============================================================Bond-Slip=================================================================
set alphaBS 0.4 ; #alphaBS is a parameter used in the local bond-slip relation and can be taken as 0.4
set Pow [expr (1.0/$alphaBS)]
#Eq for a model with units (ksi, kips, inches)
#set Sy [expr (((((($dlSec*$Fy*1000.0)/(4000.0*(sqrt($fc*-1000.0))))*((2.0*$alphaBS)+1.0))*pow($Pow,1))*0.1)+0.013)]
set Sy 0.08337;
set Su [expr $Sy*35.0]
set b 0.5
set R 0.7
set bondslipMat 5
#uniaxialMaterial Bond_SP01 $matTag $Fy $Sy $Fu $Su $b $R
uniaxialMaterial Bond_SP01 $bondslipMat $Fy $Sy $Fu $Su $b $R
puts "materials defined "
#----------------------------------------------------------------
# Generate a circular reinforced concrete section
# with one layer of steel evenly distributed around the perimeter and a confined core.
# confined core.
# by: Michael H. Scott, 2003
# Notes
# The center of the reinforcing bars are placed at the inner radius
# The core concrete ends at the inner radius (same as reinforcing bars)
# The reinforcing bars are all the same size
# The center of the section is at (0,0) in the local axis system
# Zero degrees is along section y-axis
#Following recommendations using uniformly distibutions ratdial disretization shema
set ri 0.0; # inner radius of the section, only for hollow sections
set ro [expr $Dcol/2]; # overall (outer) radius of the section
set nfCoreR 10; # number of radial divisions in the core (number of "rings")
set nfCoreT 20; # number of theta divisions in the core (number of "wedges")
set nfCoverR 1; # number of radial divisions in the cover
set nfCoverT 20; # number of theta divisions in the cover
# Define ELEMENTS & SECTIONS
#----------------------------
set ColSecTag 1; # assign a tag number to the column section
# Define the fiber section
#-------------------------
set rc [expr $ro-$coverSec]; # Core radius
set bcent 6.39 ; # expr[$coversec+$Dsec/2]
set rb [expr $ro-$bcent];
set theta [expr 360.0/$numBarsSec]; # Determine angle increment between bars
section fiberSec $ColSecTag {
#patch circ $matTag $numSubdivCirc $numSubdivRad $yCenter $zCenter $intRad $extRad <$startAng $endAng>
patch circ $IDconcC $nfCoreT $nfCoreR 0 0 $ri $rc 0 360; # Define the core patch
patch circ $IDconcU $nfCoverT $nfCoverR 0 0 $rc $ro 0 360; # Define the cover patch
#layer circ $matTag $numFiber $areaFiber $yCenter $zCenter $radius <$startAng $endAng>
layer circ $IDreinf $numBarsSec $barAreaSec 0 0 $rc; # Define the reinforcing layer
}
puts "section geometry defined of Column"
# Define the fiber section of Bond Slip
#---------------------------------------
set SecTagBS 3;
section fiberSec $SecTagBS {; # Define the fiber section
patch circ $IDconcC2 $nfCoreT $nfCoreR 0 0 $ri $rc 0 360; # Define the core patch
patch circ $IDconcU2 $nfCoverT $nfCoverR 0 0 $rc $ro 0 360; # Define the cover patch
layer circ $bondslipMat $numBarsSec $barAreaSec 0 0 $rc ; # Define the reinforcing layer
}
puts "section geometry defined of Bond Slip"
#Aggregator
#------------------------
set shearsec 4;
set TagsecAgg 6;
set AreaSec [expr $PI*pow($Dcol,2)/4]; #gross cross-sectional area (in^2)
set G [expr ($Ec)/(2.4)];
puts "G=$G"
set GA [expr ($G*$AreaSec)];
puts "GA=$GA"
#Assign shear to the model of flexure
#------------------------------------
uniaxialMaterial Elastic $shearsec 1e8
puts "shear defined"
#section Aggregator $TagsecAgg $shearsec Vy -section $ColSecTag
section Aggregator $TagsecAgg $shearsec Vy -section $SecTagBS
puts "section aggregator defined"
# define geometric transformation
#--------------------------------
#performs a linear geometric transformation of beam stiffness and resisting force
#from the basic system to the global-coordinate system
set ColTransfTag 1; # associate a tag to column transformation
set ColTransfType Linear ; # options, Linear PDelta Corotational
geomTransf $ColTransfType $ColTransfTag ;
# element connectivity:
set numIntgrPts 5; # number of integration points for force-based element
#element nonlinearBeamColumn $eleTag $iNode $jNode $numIntgrPts $secTag $transfTag
#element nonlinearBeamColumn 1 2 3 $numIntgrPts $TagsecAgg $ColTransfTag
element nonlinearBeamColumn 1 1 2 $numIntgrPts $ColSecTag $ColTransfTag
puts "element nonlinear defined"
######ZeroElement of Bond-Slip
#element zeroLengthSection $EleTag $ibNode $jbNode $SecTagBS
element zeroLengthSection 2 3 1 $TagsecAgg -orient 0 1 0 1 0 0
puts "zeroLength element defined"
puts "model etablished"
#-------------------------------------------------------------------------------
# Column- Pile element connectivity:
# Create recorder
recorder Node -file $dataDir/DFree.out -time -node 2 -dof 1 2 3 disp; # displacements of free nodes
recorder Node -file $dataDir/RBase.out -time -node 3 -dof 1 2 3 reaction; # support reaction ????? reaction node3
#recorder Element <-file $fileName> <-time> <-ele ($ele1 $ele2 ...)> <-eleRange $startEle $endEle> <-region $regTag> <-ele all> ($arg1 $arg2 ...)
#$secNum refers to the integration point whose data is to be output (In this case we want 3 Output= section 3)
recorder Element -file $dataDir/secstrSTELneg.out -time -ele 1 section 3 fiber -$ro 0. $IDreinf stressStrain #??? A revoir data ro
recorder Element -file $dataDir/secstrSTEELpos.out -time -ele 1 section 3 fiber $ro 0. $IDreinf stressStrain
recorder Element -file $dataDir/secstrCONC.out -time -ele 1 section 3 fiber $rc 0. $IDconcC stressStrain
#recorder plot $dataDir/topdispac.out Node2_Ydisp 10 10 300 300 -columns 2 1
# Static analysis
#-------------------------------------------------------------
#define GRAVITY
pattern Plain 1 Linear {
load 2 0 -$P 0
}
# Gravity-analysis parameters -- load-controlled static analysis
set Tol 1.0e-5; # convergence tolerance for test
constraints Plain; # how it handles boundary conditions
numberer Plain; # renumber dof's to minimize band-width (optimization), if you want to
system BandGeneral; # how to store and solve the system of equations in the analysis
test NormDispIncr $Tol 100 ; # determine if convergence has been achieved at the end of an iteration step
algorithm Newton; # use Newton's solution algorithm: updates tangent stiffness at every iteration
set NstepGravity 1; # apply gravity in 10 steps
set DGravity [expr 1./$NstepGravity]; # first load increment;
integrator LoadControl $DGravity; # determine the next time step for an analysis
analysis Static; # define type of analysis static or transient
analyze $NstepGravity; # apply gravity
#-----maintain constant gravity loads and reset time to zero------
loadConst -time 0.0
print -node 2
puts "Model Built"
#EIGEN MODE
#-----------------------------
set omega2 [eigen 1];
set omega1 [expr pow($omega2 ,0.5)];
set period1 [expr 2*3.1415/$omega1];
puts " OpenSees Calculated Period of First mode= $period1"