Source code

[amia]

Dear Frank

I was wondering if you could please make the source code for Pinching Limit State Material available online?
I have looked at the database available at the link below but I think the Pinching Limit State Material is not made available yet?
http://opensees.berkeley.edu/WebSVN/lis ... a8c827eb5a

I would like to have a look at the source code because I am not sure how the rotation limit is being applied in the material model (through the associated shear Limit Curve). From my understanding, the rotation limit should be applied to the nodes which are selected to monitor the rotation/drift. So the rotation should be: {(displacement of node J) - (displacement of node I)} / perpendicular distance between node J and I.
However, I have noticed that this is not the case when I am using Pinching Limit State Material, and the rotation limit appears to be reached earlier than the defined limit. Below is a simple code I have written for a cantilever column using a single elastic element for the column and a zero length spring at the bottom defined by the Pinching Limit State Material.

Thank you in advance
Anita
________________________________________________________________

# Units: N, mm, MPa, s

## Create ModelBuilder space
wipe all;
model BasicBuilder -ndm 2 -ndf 3; # 2D with 3 degress of freedom

# ----- Elastic column sections ------------------------------------------------ #
set colDepth 457.0
set colWidth 457.0
set Ec [expr 5000*sqrt(21.1)]
set Ac [expr $colWidth*$colDepth]
set Izc [expr 0.2*$colWidth*pow($colDepth,3)/12]

# ---- Nodes ------------------------------------------------------------------- #
set Lc 1473.0

node 1 0.0 0.0
node 11 0.0 0.0
node 2 0.0 $Lc

# ---- Elastic element connectivity -------------------------------------------- #
set colTag 1
set transTag 1

geomTransf PDelta $transTag
element elasticBeamColumn $colTag 11 2 $Ac $Ec $Izc $transTag

# ---- Shear Spring ------------------------------------------------------------ #
#limitCurve RotationShearCurve $crvTag $eleTag $ndI $ndJ $rotAxis $Vn $Vr $Kdeg $rotLim
set shearCurveTag 1
set eleTag $colTag
set ndIb 11
set ndJb 2
set rotAxis 3
set Vn 10.0e5
set Vr 2000.0
set Kdeg [expr -4762.0]
set rotLim [expr 0.02]
limitCurve RotationShearCurve $shearCurveTag $eleTag $ndIb $ndJb $rotAxis $Vn $Vr $Kdeg $rotLim

#uniaxialMaterial PinchingLimitStateMaterial $matTag $nodeT $nodeB $driftAxis $Kelas $crvTyp $crvTag $eleTag $b $d $h $a $st $As $Acc $ld $db $rhot $f'c $fy $fyt
set shearMatTag 100
set nodeT 11
set nodeB 2
set driftAxis 1
set Kelas -3
set crvTyp 2
#set crvTag $shearCurveTag
set eleTag $colTag
set b 457.0
set d [expr 457.0-65.13]
set h 457.0
set a [expr $Lc]
set st 350.0
set As [expr 8*645.0]
set Acc [expr 317.2*317.2]
set ld 1000.0
set db 28.65
set rhot 0.0020
set fc 21.1
set fy 438.8
set fyt 476.0

# Calibrated Model
uniaxialMaterial PinchingLimitStateMaterial $shearMatTag $nodeT $nodeB $driftAxis $Kelas $crvTyp $shearCurveTag $eleTag $b $d $h $a $st $As $Acc $ld $db $rhot [expr -$fc] $fy $fyt

# zero-length for shear response
set softMaterial 999
set rigidMaterial 998
uniaxialMaterial Elastic $softMaterial 0.9
uniaxialMaterial Elastic $rigidMaterial 0.9e15

set shearEleTag 2
set SoftEleTag 3

element zeroLength $shearEleTag 1 11 -mat $shearMatTag $rigidMaterial $rigidMaterial -dir 1 2 3
element zeroLength $SoftEleTag 1 11 -mat $softMaterial -dir 1 ; # required to get gravity analysis to run

# ---- Node restraints ------------------------------------------------------------ #
fix 1 1 1 1

# ---- Gravity Analysis ----------------------------------------------------------- #
file mkdir Results/Gravity
recorder Node -file Results/Gravity/node1Reactions.out -load -node 1 -dof 1 2 3 reaction
puts " ----------------------------------------------- "
puts "Running gravity analysis..."

pattern Plain 1 Linear {
load 2 0 -667000.0 0
}

test NormDispIncr 1e-12 150
constraints Penalty 1.0e16 1.0e16 ;
numberer RCM
algorithm Newton
system BandGeneral
set NstepGravity 100;
set DGravity [expr 1.0/$NstepGravity]
integrator LoadControl $DGravity
analysis Static
set ok [analyze $NstepGravity] ; # this will return zero if no convergence problems were encountered
if {$ok == 0} {
puts "Gravity analysis completed SUCCESSFULLY" ;
} else {
puts "Gravity analysis FAILED";
}

remove recorders ; # needed if not recorders will continue recording other type of analyses
loadConst -time 0.0


# ---- Pushover Analysis ----------------------------------------------------------- #
file mkdir Results

recorder Node -file Results/DispNode2.out -load -node 2 -dof 1 disp
recorder Node -file Results/DispNode11.out -load -node 11 -dof 1 disp
recorder Node -file Results/DispNode1.out -load -node 1 -dof 1 disp
recorder Node -file Results/BaseShear.out -load -node 1 -dof 1 2 3 reaction
recorder Element -file Results/Ele1Globalforce.out -ele 1 -dof 1 globalForce
recorder Element -file Results/EleHingeLocalforce.out -ele $colTag -dof 1 2 3 localForce

recorder Element -file Results/shearEleTag.out -ele $shearEleTag -load -dof 1 2 3 localForce
puts " ----------------------------------------------- "
puts "Running pushover analysis..."

set nodeTag 2
set Dmax [expr 150.0]
set Dincr 0.1
set Nsteps [expr int($Dmax/$Dincr)]
set ok 0
set n 1

pattern Plain 2 "Linear -factor 1" {
sp $nodeTag 1 1.0
}

integrator LoadControl $Dincr 1 $Dincr $Dincr
constraints Penalty 1.0e16 1.0e16 ;
numberer RCM ;
system UmfPack ;
test EnergyIncr 1.0e-10 1000 0
algorithm KrylovNewton ;
analysis Static ;

set ok [analyze $Nsteps]

if {$ok != 0} {
puts "PO ANALYSIS FAILED"
} else {
puts "PO SUCCESSFULL"
}

[fmk]

pensees.berkeley.edu/WebSVN/filedetails.php?repname=OpenSees&path=%2Ftrunk%2FSRC%2Fmaterial%2Funiaxial%2FlimitState%2FPinchingLimitStateMaterial.cpp&

[Jayanthan]

Hi Anita,

This is Jayanthan from UBC, Canada. I am also working on shear critical column. As you used LeBorgne's shear critical model, I too use the same. Could you please tell me how you calculated the stiffness value for rotational spring in soft material and rigid material tag?

Thanks & Regards
Jayanthan
m.jayanth1992@gmail.com

[amia]

Hi Jayanthan

My apologies for the late reply.
The stiffness for soft and rigid material are arbitrary to allow to define the material behaviour. For example the soft material is basically a material with zero stiffness, but to avoid numerical issues I've used 0.9. The rigid material is set to a very high number, usually this has to be 10 or 100 higher than the elastic stiffness of the beams or columns in your model.
I hope that helps!

May I also ask you how you are finding LeBorgne's shear critical model? I can get some models to work correctly, however, I seem to have a lot of convergence problems when the rotation limit is reached.. Once the limit is reached the response is unable to follow the defined slope for the shear spring.

[Jayanthan]

Hi Anita,

Thanks for replying my post. I am also facing the same issue. Initially, I tried LeBorgne's model for sezen column with one shear spring at the bottom of the column that seems to work fine, where I used calibrated input but I am not getting the perfect degrading slope for shear spring. Then, I introduced shear spring at top and bottom, it is not working properly.
Another question, I need to ask, Do we need to use zero length section for defining bar slip fiber section (as given the LeBorgne's paper)?
If you wish, I can also email my opensees code with you.

Thanks & Regards
Jayanthan
m.jayanth1992@gmail.com

[amia]

Hi Jayanthan

With bar-slip I think you can also model it with a spring rather than zero length section, but I think the zero length section is meant to be a more accurate method for modelling bar-slip..
As I've mentioned I have problems getting the shear spring to work consistently as well, especially for columns with higher axial load..And i have simplified my models to cantilevers for the time being. I know that LeBorgne discusses in his thesis that sometimes having two shear springs may cause numerical issues for a double curvature column if they have the same rotational limit.. Have you tried changing the rotation limit for one end just slightly?

And thank you, it may help if I have a look at your code as well!
If you would like to you may email it to me (amia@student.unimelb.edu.au) or upload it here!

Thank you
Anita

[Jayanthan]

Hi Anita,

I would like to post my script where I placed shear spring at bottom and modeled zerolength barslip as fiber section, but when I place the shear spring at top instead of placing it at bottom, the script is not working properly....

I also sent you an email regarding this long back,

# UBCO
# Sezen (2002)_Specimen1_Flexure Shear Critical Column
# unit kip,inches
# SET UP ------------------------------------------------

# Tags.tcl
set coreTag 1 ;# core concrete
set coverTag 2 ;# cover concrete
set steelTag 3 ;# steel

set coreTagc 6 ;# core concrete
set coverTagc 7 ;# cover concrete
set steelTagc 8 ;# steel

set rigidMatTag 12;
set flexSec 1;
set flexSecc 2;

#--------------------------------------------------------
# Build model
#--------------------------------------------------------
puts "Starting model"
wipe;
source DisplayModel2D.tcl; #procedure for displaying a 2D perspective of model
source DisplayPlane_IR.tcl; #procedure for displaying a plane in a model
set dataDir Sezen(2002)_Specimen1_Result
file mkdir $dataDir;
model BasicBuilder -ndm 2 -ndf 3;

#--------------------------------------------------------
# Define nodal mesh and B.C.s
#--------------------------------------------------------
set L 116;

# tag X Y
node 1 0.0 0.0
node 2 0.0 0.0
node 3 0.0 18
#node 4 0.0 98
#node 5 0.0 $L
node 6 0.0 $L

# tag DX DY RZ
fix 1 1 1 1

#--------------------------------------------------------
# Create column section
#--------------------------------------------------------
# CenterColSecFiber.tcl

set h 18
set b 18

# Set parameters for fiber section
set cover 2;
set Y1 [expr $h/2.0]
set Z1 [expr $b/2.0]

#-----------------------------------------------------------------------------------------
# material for column element
#-----------------------------------------------------------------------------------------
# nominal concrete compressive strength
set fc [expr -3.059]; # CONCRETE Compressive Strength, ksi (+Tension, -Compression)
set Ec [expr 57*sqrt(3059)]; # Concrete Elastic Modulus
# confined concrete
set Kfc 1.1; # ratio of confined to unconfined concrete strength
set fc1C [expr $Kfc*$fc]; # CONFINED concrete (mander model), maximum stress
set eps1C [expr 0.002*(1+5*($Kfc-1))]; #strain at maximum stress
set fc2C [expr 0.9*$fc1C]; # ultimate stress
set eps2C [expr 0.009];# strain at ultimate stress
#set eps2C [expr 3*$eps1C]; # strain at ultimate stress
# unconfined concrete
set fc1U $fc; # UNCONFINED concrete (todeschini parabolic model), maximum stress
set eps1U -0.002; # strain at maximum strength of unconfined concrete
set fc2U [expr ($fc1U/3)]; # ultimate stress
set eps2U -0.006; # strain at ultimate stress
set lambda 3; # ratio between unloading slope at $eps2 and initial slope $Ec
# tensile-strength properties
#set ftC [expr 7.5*sqrt(3.670)]; # tensile strength +tension
#set ftU [expr 7.5*sqrt(3.059)]; # tensile strength +tension
#set Etc [expr (2450/5)]; # tension softening stiffness
#set Etu [expr (3060/5)];
#set Etcb [expr (150/5)]
#set Etub [expr (190/5)]
# -----------
set Fy [expr 63.5]; # STEEL yield stress
set Es [expr 29000]; # modulus of steel
set Bs 0.01; # strain-hardening ratio
set R0 18; # 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
#uniaxialMaterial Concrete02 $coreTag $fc1C $eps1C $fc2C $eps2C $lambda $ftC $Etc; # build core concrete (confined)
#uniaxialMaterial Concrete02 $coverTag $fc1U $eps1U $fc2U $eps2U $lambda $ftU $Etu; # build cover concrete (unconfined)

uniaxialMaterial Concrete01 $coreTag $fc1C $eps1C $fc2C $eps2C; # build core concrete (confined)
uniaxialMaterial Concrete01 $coverTag $fc1U $eps1U $fc2U $eps2U; # build cover concrete (unconfined)
uniaxialMaterial Steel02 $steelTag $Fy $Es $Bs $R0 $cR1 $cR2; # build reinforcement material

#-----------------------------------------------------------------------------------------
#Define the fiber section

set barAreaCol 0.99;
section Fiber $flexSec {

# Create the concrete core fibers
patch rect $coreTag 10 1 [expr $cover-$Y1] [expr $cover-$Z1] [expr $Y1-$cover] [expr $Z1-$cover]; # Define the concrete core patch

# Create the concrete cover fibers (top, bottom, left, right)
patch rect $coverTag 10 1 [expr -$Y1] [expr $Z1-$cover] $Y1 $Z1
patch rect $coverTag 10 1 [expr -$Y1] [expr -$Z1] $Y1 [expr $cover-$Z1]
patch rect $coverTag 2 1 [expr -$Y1] [expr $cover-$Z1] [expr $cover-$Y1] [expr $Z1-$cover]
patch rect $coverTag 2 1 [expr $Y1-$cover] [expr $cover-$Z1] $Y1 [expr $Z1-$cover]

# create the reinforcing fibers
layer straight $steelTag 3 $barAreaCol [expr $Y1-$cover] [expr $Z1-$cover] [expr $Y1-$cover] [expr $cover-$Z1]
layer straight $steelTag 2 $barAreaCol 0.0 [expr $Z1-$cover] 0.0 [expr $cover-$Z1]
layer straight $steelTag 3 $barAreaCol [expr $cover-$Y1] [expr $Z1-$cover] [expr $cover-$Y1] [expr $cover-$Z1]
}

#zerolength fiber section material properties
## Core concrete (confined)
uniaxialMaterial Concrete01 $coreTagc $fc1C -0.0500 $fc2C -0.1500 ; # build core concrete (confined)
uniaxialMaterial Concrete01 $coverTagc $fc1U -0.0333 $fc2U -0.1000; # build cover concrete (unconfined)
uniaxialMaterial Steel02 $steelTagc 80 1800 0.01 $R0 $cR1 $cR2; # build reinforcement material

section Fiber $flexSecc {

# Create the concrete core fibers (Confined)
patch rect $coreTagc 10 1 [expr $cover-$Y1] [expr $cover-$Z1] [expr $Y1-$cover] [expr $Z1-$cover]

# Create the concrete cover fibers (top, bottom, left, right)
patch rect $coverTagc 10 1 [expr -$Y1] [expr $Z1-$cover] $Y1 $Z1
patch rect $coverTagc 10 1 [expr -$Y1] [expr -$Z1] $Y1 [expr $cover-$Z1]
patch rect $coverTagc 2 1 [expr -$Y1] [expr $cover-$Z1] [expr $cover-$Y1] [expr $Z1-$cover]
patch rect $coverTagc 2 1 [expr $Y1-$cover] [expr $cover-$Z1] $Y1 [expr $Z1-$cover]

# Create the reinforcing fibers ( top, middle, bottom)
# explanations for the reinforcement placement
# symmetric section
# y
# ^
# |
# --------------------- -- --
# | o o o | | -- cover
# | | |
# | | |
# z <--- | o + o | H
# | | |
# | | |
# | o o o | | -- cover
# --------------------- -- --
# |-------- B --------|
layer straight $steelTagc 3 $barAreaCol [expr $Y1-$cover] [expr $Z1-$cover] [expr $Y1-$cover] [expr $cover-$Z1]
layer straight $steelTagc 2 $barAreaCol 0.0 [expr $Z1-$cover] 0.0 [expr $cover-$Z1]
layer straight $steelTagc 3 $barAreaCol [expr $cover-$Y1] [expr $Z1-$cover] [expr $cover-$Y1] [expr $cover-$Z1]
}

#--------------------------------------------------------
# Define the beam-column element
#--------------------------------------------------------
geomTransf PDelta 1


#Distributed Plasticity Elements
set nint 5
element zeroLengthSection 2 1 2 $flexSecc -orient 0 1 0 1 0 0;
element nonlinearBeamColumn 11 2 3 $nint $flexSec 1 -iter 6 1e-12
element nonlinearBeamColumn 12 3 6 $nint $flexSec 1 -iter 6 1e-12
#element forceBeamColumn 13 4 5 $nint $flexSec 1 -iter 5 1e-12

#Calibrated Input for shear spring--------------------------------
#limitCurve RotationShearCurve $crvTag $eleTag $ndI $ndJ $rotAxis $Vn $Vr $Kdeg $defType $b $d $h $L $st $As $Acc $ld $db $rhot $f'c $fy $fyt $delta
limitCurve RotationShearCurve 2 11 1 3 3 -1 -1 0 1 18 15.44 18 116 12 7.99 256 56 1.128 0.0015 -3.059 63.5 69.02 0.0

#uniaxialMaterial PinchingLimitStateMaterial $matTag $nodeT $nodeB $driftAxis $Kelas $crvTyp $crvTag $eleTag $b $d $h $a $st $As $Acc $ld $db $rhot $f'c $fy $fyt
uniaxialMaterial PinchingLimitStateMaterial 9 6 2 1 -3 2 2 11 18 15.44 18 58 12 7.99 256 56 1.128 0.0015 -3.059 63.5 69.02
#limitCurve RotationShearCurve 3 13 6 4 3 0 -1 0 1 18 15.44 18 116 12 7.99 256 56 1.128 0.0017 -3.059 63.5 69.02 0.0
#uniaxialMaterial PinchingLimitStateMaterial 10 5 1 1 -3 2 3 13 18 15.44 18 58 12 7.99 256 56 1.128 0.0017 -3.059 63.5 69.02

#--------------------------------------------------------
# Define the zero-length end springs
#--------------------------------------------------------

#bottom of column slip spring

element zeroLength 1 1 2 -mat 9 -dir 1

#--------------------------------------------------------
# Define GRAVITY Load
#--------------------------------------------------------
# Define Axial Load
set P [expr 150]

timeSeries Linear 1
pattern Plain 1 1 {
load 6 0 [expr -$P] 0
}

# Gravity-analysis parameters -- load-controlled static analysis
set Tol 1.0e-6; # 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 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./$NstepGravity]; # first load increment;
integrator LoadControl $DGravity; # determine the next time step for an analysis
analysis Static; # define type of analysis static or transient
puts "Check"
analyze $NstepGravity; # apply gravity
# ------------------------------------------------- maintain constant gravity loads and reset time to zero
puts "Check"
loadConst -time 0.0

puts "Model Built"
# display deformed shape:
set ViewScale 5;
DisplayModel2D DeformedShape $ViewScale ; # display deformed shape, the scaling factor needs to be adjusted for each model

#--------------------------------------------------------
# Recorders
#--------------------------------------------------------

# record displacements for node with applied displacement
recorder Node -file $dataDir/TopDisp.out -time -node 6 -dof 1 disp
recorder Node -file $dataDir/RX.out -time -node 1 -dof 1 reaction

source Cyclic.tcl


Regards,
Jayanthan

[amia]

Hi Jayanthan

I have not received an email from you, have you sent it to amia@student.unimelb.edu.au ?

Your script looks correct to me. The only thing that I would suggest is to define the zero length element behaviour in the other two directions as well (i.e Y and rotation), assign a rigid material to these directions.. I think this might help to get your model to work when you place the spring at the top of the column. Also if you are modelling a double curvature column, then you should fix the top node in the rotational direction as well.

I have actually stopped using the pinching limit state material as I could not overcome the convergence issues! And at this stage I have moved onto a lumped plasticity approach..
But, please let me know how you go with your modelling as I am still very interested!

All the best
Anita

[Jayanthan]

Hi Anita,

Thank you for your reply.. I held up with other work so I didn't see the forum.. sorry for that.. I sent an email to you, please check your spam folder as well.

Regards,
Jayanthan