Zero length element in a multistory concrete frame

[ele_cam]

Hi all,

I'm trying to model a concrete multistory frame by using zero length elements for the plastic hinges at the ends of both columns and beams. I don't know if this is correct since the structure apparently has some lability (neither the eigen value analysis is performed correctly). Is this because three zero length elements converge at the same node i.e. the intersection between the columns and the beam?
If I fix by "equalDOF" the translations and the rotations of the zero length elements then I get the correct results as if the frame was elastic, otherwise the pushover analysis is not performed... does anybody can help me?

Thank you!
I attach the source code in the following.

###################################################################################################
wipe all; # clear memory of past model definitions
model BasicBuilder -ndm 2 -ndf 3; # Define the model builder, ndm = #dimension, ndf = #dofs
source LibUnits_Ele.tcl; # list of units of measure
source DisplayModel2D.tcl; # procedure for displaying a 2D perspective of model
source DisplayPlane.tcl; # procedure for displaying a plane in a model

###################################################################################################
# Define type of analysis: "pushover" = pushover
set analysisType "pushover";
if {$analysisType == "pushover"} {
set dataDir OutputPO; # name of output folder
file mkdir $dataDir; # create output folder
}

###################################################################################################
# define structure-geometry parameters
set NStories 3; # number of stories
set NBays 1; # number of frame bays (excludes bay for P-delta column)
set WBay [expr 8]; # bay width
set HStory [expr 4]; # story height
set HBuilding [expr $HStory*$NStories]; # height of building

# calculate locations of beam/column joints:
set Pier1 0.0; # leftmost column line
set Pier2 [expr $Pier1 + $WBay];
set Floor1 0.0; # ground floor
set Floor2 [expr $Floor1 + $HStory];
set Floor3 [expr $Floor2 + $HStory];
set Floor4 [expr $Floor3 + $HStory];

# calculate nodal masses -- lump floor masses at frame nodes
set Floor2Weight 500.0; # weight of Floor 2 in kN
set Floor3Weight 500.0; # weight of Floor 3 in kN
set Floor4Weight 500.0; # weight of Floor 4 in kN
set WBuilding [expr $Floor2Weight + $Floor3Weight + $Floor4Weight]; # total building weight
set NodalMass2 [expr ($Floor2Weight/$g) / (2.0)]; # mass at each node on Floor 2
set NodalMass3 [expr ($Floor3Weight/$g) / (2.0)]; # mass at each node on Floor 3
set NodalMass4 [expr ($Floor4Weight/$g) / (2.0)]; # mass at each node on Floor 4

# define nodes and assign masses to beam-column intersections of frame
# command: node nodeID xcoord ycoord -mass mass_dof1 mass_dof2 mass_dof3
# nodeID convention: "xy" where x = Pier # and y = Floor #
node 11 $Pier1 $Floor1;
node 21 $Pier2 $Floor1;
node 12 $Pier1 $Floor2 -mass $NodalMass2 $Usmall $Usmall;
node 22 $Pier2 $Floor2 -mass $NodalMass2 $Usmall $Usmall;
node 13 $Pier1 $Floor3 -mass $NodalMass3 $Usmall $Usmall;
node 23 $Pier2 $Floor3 -mass $NodalMass3 $Usmall $Usmall;
node 14 $Pier1 $Floor4 -mass $NodalMass4 $Usmall $Usmall;
node 24 $Pier2 $Floor4 -mass $NodalMass4 $Usmall $Usmall;
node 117 $Pier1 $Floor1;
node 217 $Pier2 $Floor1;
node 126 $Pier1 $Floor2;
node 122 $Pier1 $Floor2;
node 127 $Pier1 $Floor2;
node 226 $Pier2 $Floor2;
node 223 $Pier2 $Floor2;
node 227 $Pier2 $Floor2;
node 136 $Pier1 $Floor3;
node 132 $Pier1 $Floor3;
node 137 $Pier1 $Floor3;
node 236 $Pier2 $Floor3;
node 233 $Pier2 $Floor3;
node 237 $Pier2 $Floor3;
node 146 $Pier1 $Floor4;
node 142 $Pier1 $Floor4;
node 246 $Pier2 $Floor4;
node 243 $Pier2 $Floor4;

# constrain beam-column joints in a floor to have the same lateral displacement using the "equalDOF" command
# command: equalDOF $MasterNodeID $SlaveNodeID $dof1 $dof2...
set dof1 1; # constrain movement in dof 1 (x-direction)
equalDOF 12 22 $dof1; # Floor 2: Pier 1 to Pier 2
equalDOF 13 23 $dof1; # Floor 3: Pier 1 to Pier 2
equalDOF 14 24 $dof1; # Floor 3: Pier 1 to Pier 2

# assign boundary condidtions
# command: fix nodeID dxFixity dyFixity rzFixity
# fixity values: 1 = constrained; 0 = unconstrained
# fix the base of the building;
fix 11 1 1 1;
fix 21 1 1 1;

###################################################################################################
# define material properties
set fy [expr 430.*$MPa]; # steel strength at yielding
set Es [expr 206000.*$MPa]; # steel Young's modulus
set Esh 0; # steel Young's modulus at hardening
set esh [expr $fy/$Es]; # strain corresponding to initial strain hardening
set esu 0.06; # strain at peak stress
set fc [expr -35.*$MPa]; # concrete compressive strength
set ec -0.002; # concrete strain at maximum strength
set ecu -0.0035; # concrete strain at crushing strength
set Ec [expr 9500*pow(abs($fc/$MPa),(1./3.))*$MPa]; # concrete Young's modulus
set fct [expr 0.25*pow(abs($fc/$MPa),(2./3.))*$MPa]; # concrete maximum tensile strength
set et [expr abs($ecu/2.)]; # concrete ultimate tensile strain
set EIbeam [expr 0.5*$Ec*0.5*(0.8**3)/12]; # initial tangent of the bending moment behavior of beams

# define elastic section properties

# define column section 70x70
set HCol 0.7
set BCol 0.7
set HCol2 [expr $HCol/2]; # Mid Column Depth
set BCol2 [expr $BCol/2]; # Mid Column Width
set cc 0.04; # concrete cover

set Acol [expr $HCol*$BCol]; # cross-sectional area
set Icol [expr (0.7**4)/12]; # moment of inertia

# define beam section 50x80
set Abeam 0.4; # cross-sectional area
set Ibeam [expr 0.5*(0.8**3)/12]; # moment of inertia

# set up geometric transformations of element
set PDeltaTransf 1;
geomTransf PDelta $PDeltaTransf; # PDelta transformation

# define column and beam elements (elastic part)
element elasticBeamColumn 111 117 126 $Acol $Ec $Icol $PDeltaTransf;
element elasticBeamColumn 121 217 226 $Acol $Ec $Icol $PDeltaTransf;
element elasticBeamColumn 112 127 136 $Acol $Ec $Icol $PDeltaTransf;
element elasticBeamColumn 122 227 236 $Acol $Ec $Icol $PDeltaTransf;
element elasticBeamColumn 113 137 146 $Acol $Ec $Icol $PDeltaTransf;
element elasticBeamColumn 123 237 246 $Acol $Ec $Icol $PDeltaTransf;

element elasticBeamColumn 212 122 223 $Abeam $Ec $Ibeam $PDeltaTransf;
element elasticBeamColumn 213 132 233 $Abeam $Ec $Ibeam $PDeltaTransf;
element elasticBeamColumn 214 142 243 $Abeam $Ec $Ibeam $PDeltaTransf;

# ROTATIONAL SPRINGS

# BEAMS rotational springs
set matTagB1 21
set matTagB2 22
set matTagB3 23
set secTagB1 221
set secTagB2 222
set secTagB3 223

set Mbeam1 954.4; # resistant bending moment of beams 1st story - High Ductility
set Mbeam2 834.0; # resistant bending moment of beams 2nd story - High Ductility
set Mbeam3 386.8; # resistant bending moment of beams 3rd story - High Ductility
set Lpbeam [expr 0.1*$WBay]; # length of the plastic hinge
# uniaxialMaterial Steel01 $matTag $Fy $E0 $b <$a1 $a2 $a3 $a4>
uniaxialMaterial Steel01 $matTagB1 $Mbeam1 [expr $Mbeam1/(0.00535*$Lpbeam)] $Usmall
uniaxialMaterial Steel01 $matTagB2 $Mbeam2 [expr $Mbeam2/(0.00535*$Lpbeam)] $Usmall
uniaxialMaterial Steel01 $matTagB3 $Mbeam3 [expr $Mbeam3/(0.00535*$Lpbeam)] $Usmall

# BEAMS sections
# section Aggregator $secTag $matTag1 $dof1 $matTag2 $dof2 ....... <-section $sectionTag>
section Aggregator $secTagB1 $matTagB1 Mz; # bending behavior about z axis of the beam 1st story
section Aggregator $secTagB2 $matTagB2 Mz; # bending behavior about z axis of the beam 2nd story
section Aggregator $secTagB3 $matTagB3 Mz; # bending behavior about z axis of the beam 3rd story

# COLUMNS rotational springs
set matTagC1 11
set matTagC2 12
set matTagC3 13
set secTagC1 111
set secTagC2 112
set secTagC3 113

set Mcolumn1 630.; # resistant bending moment of columns 1st story
set Mcolumn2 560.; # resistant bending moment of columns 2nd story
set Mcolumn3 485.; # resistant bending moment of columns 3rd story
set Kc1 0.005; # yield curvature of columns 1st story
set Kc2 0.0053; # yield curvature of columns 2nd story
set Kc3 0.0051; # yield curvature of columns 3rd story
set Lpcol [expr 0.1*$HStory]; # length of the plastic hinge
# uniaxialMaterial Steel01 $matTag $Fy $E0 $b <$a1 $a2 $a3 $a4>
uniaxialMaterial Steel01 $matTagC1 $Mcolumn1 [expr $Mcolumn1/($Kc1*$Lpcol)] $Usmall
uniaxialMaterial Steel01 $matTagC2 $Mcolumn2 [expr $Mcolumn2/($Kc2*$Lpcol)] $Usmall
uniaxialMaterial Steel01 $matTagC3 $Mcolumn3 [expr $Mcolumn3/($Kc3*$Lpcol)] $Usmall

# BEAMS sections
# section Aggregator $secTag $matTag1 $dof1 $matTag2 $dof2 ....... <-section $sectionTag>
section Aggregator $secTagC1 $matTagC1 Mz; # bending behavior about z axis of the columns 1st story
section Aggregator $secTagC2 $matTagC2 Mz; # bending behavior about z axis of the columns 1st story
section Aggregator $secTagC3 $matTagC3 Mz; # bending behavior about z axis of the columns 1st story

# define plastic hinge elements

element zeroLength 3111 11 117 -mat $matTagC1 -dir 3;
element zeroLength 3211 21 217 -mat $matTagC1 -dir 3;
element zeroLength 3112 126 12 -mat $matTagC1 -dir 3;
element zeroLength 3212 226 22 -mat $matTagC1 -dir 3;
element zeroLength 4121 12 122 -mat $matTagB1 -dir 3; # beam 1st
element zeroLength 4122 223 22 -mat $matTagB1 -dir 3; # beam 1st
element zeroLength 3121 12 127 -mat $matTagC2 -dir 3;
element zeroLength 3221 22 227 -mat $matTagC2 -dir 3;
element zeroLength 3122 136 13 -mat $matTagC2 -dir 3;
element zeroLength 3222 236 23 -mat $matTagC2 -dir 3;
element zeroLength 4131 13 132 -mat $matTagB2 -dir 3; # beam 2nd
element zeroLength 4132 233 23 -mat $matTagB2 -dir 3; # beam 2nd
element zeroLength 3131 13 137 -mat $matTagC3 -dir 3;
element zeroLength 3231 23 237 -mat $matTagC3 -dir 3;
element zeroLength 3132 146 14 -mat $matTagC3 -dir 3;
element zeroLength 3232 246 24 -mat $matTagC3 -dir 3;
element zeroLength 4141 14 142 -mat $matTagB3 -dir 3; # beam 3rd
element zeroLength 4142 243 24 -mat $matTagB3 -dir 3; # beam 3rd

equalDOF 11 117 1 2 ;
equalDOF 21 217 1 2 ;
equalDOF 12 126 1 2 ;
equalDOF 22 226 1 2 ;
equalDOF 12 122 1 2 ; # beam 1st
equalDOF 22 223 1 2 ; # beam 1st
equalDOF 12 127 1 2 ;
equalDOF 22 227 1 2 ;
equalDOF 13 136 1 2 ;
equalDOF 23 236 1 2 ;
equalDOF 13 132 1 2 ; # beam 2nd
equalDOF 23 233 1 2 ; # beam 2nd
equalDOF 13 137 1 2 ;
equalDOF 23 237 1 2 ;
equalDOF 14 146 1 2 ;
equalDOF 24 246 1 2 ;
equalDOF 14 142 1 2 ; # beam 3rd
equalDOF 24 243 1 2 ; # beam 3rd

# display the model with the node numbers
DisplayModel2D NodeNumbers

############################################################################
set pi [expr 2.0*asin(1.0)]; # Definition of pi
set nEigenI 1; # mode i = 1
set nEigenJ 2; # mode j = 2
set nEigenK 3; # mode k = 3
set lambdaN [eigen [expr $nEigenK]]; # eigenvalue analysis for nEigenK modes
set lambdaI [lindex $lambdaN [expr 0]]; # eigenvalue mode i = 1
set lambdaJ [lindex $lambdaN [expr $nEigenK-2]]; # eigenvalue mode j = 2
set lambdaK [lindex $lambdaN [expr $nEigenK-1]]; # eigenvalue mode j = 2
set w1 [expr pow($lambdaI,0.5)]; # w1 (1st mode circular frequency)
set w2 [expr pow($lambdaJ,0.5)]; # w2 (2nd mode circular frequency)
set w3 [expr pow($lambdaK,0.5)]; # w3 (3rd mode circular frequency)
set T1 [expr 2.0*$pi/$w1]; # 1st mode period of the structure
set T2 [expr 2.0*$pi/$w2]; # 2nd mode period of the structure
set T3 [expr 2.0*$pi/$w3]; # 3rd mode period of the structure
puts "T1 = $T1 s"; # display the first mode period in the command window
puts "T2 = $T2 s"; # display the second mode period in the command window
puts "T3 = $T3 s"; # display the second mode period in the command window

############################################################################
# apply gravity loads
#command: pattern PatternType $PatternID TimeSeriesType
pattern Plain 101 Constant {

# point loads on leaning column nodes
# command: load node Fx Fy Mz
set P_P [expr -250]; # Peso proprio del solaio ripartito su ciascuna colonna
load 12 0.0 $P_P 0.0; # Floor 2
load 22 0.0 $P_P 0.0; # Floor 2
load 13 0.0 $P_P 0.0; # Floor 3
load 23 0.0 $P_P 0.0; # Floor 3
load 14 0.0 $P_P 0.0; # Floor 4
load 24 0.0 $P_P 0.0; # Floor 4
}

# Gravity-analysis: load-controlled static analysis
set Tol 1.0e-6; # convergence tolerance for test
constraints Plain; # how it handles boundary conditions
numberer RCM; # renumber dof's to minimize band-width (optimization)
system BandGeneral; # how to store and solve the system of equations in the analysis (large model: try UmfPack)
test NormDispIncr $Tol 6; # 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 10; # apply gravity in 10 steps
set DGravity [expr 1.0/$NstepGravity]; # 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
puts "Model Built"

############################################################################
if {$analysisType == "pushover"} {
puts "Running Pushover..."
# assign lateral loads and create load pattern:
set lat 100; # lateral force at each floor
pattern Plain 200 Linear {
load 12 $lat 0.0 0.0;
load 13 [expr 2*$lat] 0.0 0.0;
load 14 [expr 3*$lat] 0.0 0.0;
}

# display deformed shape:
set ViewScale 0.5;
DisplayModel2D DeformedShape $ViewScale ; # display deformed shape, the scaling factor needs to be adjusted for each model

# displacement parameters
set IDctrlNode 13; # node where disp is read for disp control
set IDctrlDOF 1; # degree of freedom read for disp control (1 = x displacement)
set Dmax [expr 0.05*$HBuilding]; # maximum displacement of pushover: 10% roof drift
set Dincr [expr 0.0005]; # displacement increment

# analysis commands
constraints Plain; # how it handles boundary conditions
numberer RCM; # renumber dof's to minimize band-width (optimization)
system BandGeneral; # how to store and solve the system of equations in the analysis (large model: try UmfPack)
test NormUnbalance 1.0e-6 400; # tolerance, max iterations
algorithm Newton; # use Newton's solution algorithm: updates tangent stiffness at every iteration
integrator DisplacementControl $IDctrlNode $IDctrlDOF $Dincr; # use displacement-controlled analysis
# integrator LoadControl $lambda <$numIter $minLambda $maxLambda> -- qualche problema sul controllo in forza
# set lambda 0.1;
# integrator LoadControl $lambda;
analysis Static; # define type of analysis: static for pushover
set Nsteps [expr int($Dmax/$Dincr)];# number of pushover analysis steps
# set ok [analyze $Nsteps]; # this will return zero if no convergence problems were encountered
# if I wanna try other algorithms before convergence not reached...
set ok 0;
set Incr 0;

while {$ok == 0 && $Incr<$Nsteps} {

set ok [analyze 1]

# if the analysis fails try another method and than stop it
if {$ok != 0} {
# if the analysis fails try initial tangent iteration
puts "regular newton failed .. lets try an initial stiffness for this step"
test NormDispIncr 1.0e-12 1000
algorithm ModifiedNewton -initial
set ok [analyze 1]
if {$ok == 0} {
puts "that worked .. back to regular newton"
} elseif {$ok != 0} {
puts "Convergence not reached!"
break
}
test NormUnbalance 1.0e-6 400
algorithm Newton
# continue

}

set Incr [expr $Incr+1]

}
puts "Pushover complete"; # display this message in the command window
}

[vesna]

I do not understand your question? Try to describe it better.

[ele_cam]

Dear Vesna,

I try to explain the problem better. I modeled a three story concrete frame by using elasticBeamColumn elements for beams and columns. Then I used zeroLength elements in order to model the non linear behavior of the plastic hinges that can develop at the top/bottoms of the columns and at the ends of the beams. In order to define these elements I duplicated the nodes at the intersection of columns and beams.

If I fix all the degrees of freedom of these nodes by master-slave relationships (equalDOF command) then I get the results of the pushover as if the frame was elastic, so I assume that the model is correct. Then if I release the rotation of these nodes I don't get any result.

So the question is: am I doing right by inserting the zeroLength elements?
Thank you very much!

[gerber]

I suggest you to use Beam With Hinges Element for concentrated plasticity.

[ele_cam]

Hi gerber,

I already modeled the frame with Beam With Hinges Elements and it worked. But I would like to try another kind of modeling where the flexural behavior of elements is not influenced by the changing in the axial force loading the elements during the analysis. So I tried the zero length elements where I can define a priori the bending capacity -- just to make a comparison!

[vesna]

Hi ele_cam,

Here is a good example for you to look at: http://opensees.berkeley.edu/wiki/index ... ment_Frame

I believe it will help you find the errors in your model.