Axial load in forceBeamColumn

[xantis85]

Hi,
I've a problem with the axial load in a forceBeamColumn. I'm using this element to model the nonlinear beam-column in series with the limit state material, but the applied axial load varies during the analysis although the element has not yet failed axially. I'm using the displacement control mode for the pushover, if I use the load control, the axial load remains constant but the displacement history is not what I had imposed. If I delete the limit state material from the script the problem remains. There is something wrong in the integrators or algorithm? Thank you in advance.

# Tags.tcl


set coreTag 1 ;# core concrete
set coverTag 2 ;# cover concrete
set steelTag 3 ;# steel
set shearTag 4 ;# shear limit state material
set momTag 5 ;# mom-curv hysteretic model
set axialTag 6 ;# elastic axial force-strain model
set momDegTag 7 ;# degrading moment-curv hysteretic model
set axialFailTag 8 ;# axial limit state material
set centerSlipTag 9 ;# elastic slip spring
set rigidMatTag 10
set softMatTag 11
# Section tags
set flexSec 1
set shearSec 2
set flexTopSec 3
set flexBotSec 4
set axialSec 5
set shearAxialSec 6
set axialOnlySec 7
set flexShearSec 8
set shearAxialOnlySec 9
# Limit Curve tags
set shearCurveTag 1
set axialCurveTag 2
# Element tags
set bcTag 99

set iterazione 10;



# LibUnits.tcl


set in 1.; # define basic units -- output units
set kip 1.; # define basic units -- output units
set sec 1.; # define basic units -- output units
set LunitTXT "inch"; # define basic-unit text for output
set FunitTXT "kip"; # define basic-unit text for output
set TunitTXT "sec"; # define basic-unit text for output
set ft [expr 12.*$in]; # define engineering units
set ksi [expr $kip/pow($in,2)];
set psi [expr $ksi/1000.];
set lbf [expr $psi*$in*$in]; # pounds force
#puts "lbf is $lbf"
set pcf [expr $lbf/pow($ft,3)]; # pounds per cubic foot
#puts "pcf is $pcf"
set psf [expr $lbf/pow($ft,3)]; # pounds per square foot
set in2 [expr $in*$in]; # inch^2
set in4 [expr $in*$in*$in*$in]; # inch^4
set cm [expr $in/2.54]; # centimeter, needed for displacement input in MultipleSupport excitation
set PI [expr 2*asin(1.0)]; # define constants
set g [expr 32.2*$ft/pow($sec,2)]; # gravitational acceleration
set Ubig 1.e10; # a really large number
set Usmall [expr 1/$Ubig]; # a really small number


#MODEL



set HCol 14; # rectangular-Column width
set BCol 20;
set LCol1 151; # column height for first floor senza considerare il joint
set nSteps 6000
set dlamda 0.1
###########################
# Build model
###########################
model BasicBuilder -ndm 2 -ndf 3
set dataDir PushoverCiclica
file mkdir $dataDir
################################
# Define nodal mesh and B.C.s
################################
set L $LCol1
# tag X Y
node 2 0.0 0.0
node 3 0.0 $L
node 4 0.0 $L
node 5 0.0 $L

# tag DX DY RZ
fix 2 1 1 1;#messo adesso
fix 4 0 0 1
fix 5 1 1 1
##############################
# Create column section
##############################



# General Material parameters


set G $Ubig; # make stiff shear modulus
set J 1.0; # torsional section stiffness (G makes GJ large)
set GJ [expr $G*$J];

set AgCol [expr $HCol*$BCol]; # rectangular-Column cross-sectional area
set AeCol [expr ($HCol-2.75)*$BCol];
set IzCol [expr 0.5*1./12*$BCol*pow($HCol,3)]; # about-local-z Rect-Column gross moment of inertial
# -----------------------------------# confined and unconfined CONCRETE#--------------------------------------------------
# confined concrete
set Kfc 1.3; # ratio of confined to unconfined concrete strength
set Kres 0.2; # ratio of residual/ultimate to maximum stress
set lambda 0.1;
uniaxialMaterial Elastic 1111 9.9e9;
set FC 5;
# nominal concrete compressive strength
set fc [expr -$FC*$ksi];
set Ec [expr 57*$ksi*sqrt(-$fc/$psi)];
set fc1C [expr $Kfc*$fc]; # CONFINED concrete (mander model), maximum stress
set eps1C [expr 2.*$fc1C/$Ec]; # strain at maximum stress
set fc2C [expr $Kres*$fc1C]; # ultimate stress
set eps2C [expr 20*$eps1C]; # strain at ultimate stress
# ratio between unloading slope at $eps2 and initial slope $Ec
# unconfined concrete
set fc1U $fc; # UNCONFINED concrete (todeschini parabolic model), maximum stress
set eps1U -0.003; # strain at maximum strength of unconfined concrete
set fc2U [expr $Kres*$fc1U]; # ultimate stress
set eps2U -0.01; # strain at ultimate stress
# tensile-strength properties
set ftC [expr -0.14*$fc1C]; # tensile strength +tension
set ftU [expr -0.14*$fc1U]; # tensile strength +tension
set Ets [expr $ftU/0.002]; # tension softening stiffness
set KAxial [expr $AgCol*$Ec/$L]
uniaxialMaterial Elastic 110 $KAxial +0.000000E+000
uniaxialMaterial ElasticPP 120 $KAxial 0.001

set IDconcCore 2567
set IDconcCover 2568

uniaxialMaterial Concrete02 $IDconcCore $fc1C $eps1C $fc2C $eps2C $lambda $ftC $Ets; # Core concrete (confined)
uniaxialMaterial Concrete02 $IDconcCover $fc1U $eps1U $fc2U $eps2U $lambda $ftU $Ets; # Cover concrete (unconfined)

uniaxialMaterial ElasticPP 5555 1700 0.001

#set Fy [expr 66.8*$ksi]; # STEEL yield stress
set Fy40 [expr 40.*$ksi]; # yield stress for 40 ksi steel (beam reinforcement)
set Fy60 [expr 60.*$ksi]; # yield stress for 60 ksi steel (column reinforcement)

set IDsteel40 40; # ID tag for 40 ksi steel (beams & slabs)
set IDsteel60 60; # ID tag for 60 ksi steel (columns)


set Es [expr 29000.*$ksi]; # 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 Steel02 $IDsteel40 $Fy40 $Es $Bs $R0 $cR1 $cR2;
uniaxialMaterial Steel02 $IDsteel60 $Fy60 $Es $Bs $R0 $cR1 $cR2

set steelID $IDsteel60;

#SECTION FIBER

set col_cover [expr 2.75*$in]; # rectangular-RC-Column cover
set nfCoreY 20; # number of fibers in the core patch in the y direction
set nfCoreZ 20; # number of fibers in the core patch in the z direction
set nfCoverY 20; # number of fibers in the cover patches with long sides in the y direction
set nfCoverZ 20; # number of fibers in the cover patches with long sides in the z direction

set coverH 2.75;
set coverB 2.75;
set numBarsIntTot 0;
set barAreaInt 1;
set barAreaBot 1;
set barAreaTop 1;
set numBarsTop 4;
set numBarsBot 4;
set HSec $HCol
set BSec $BCol

set coverY [expr $HSec/2.0]; # The distance from the section z-axis to the edge of the cover concrete -- outer edge of cover concrete
set coverZ [expr $BSec/2.0]; # The distance from the section y-axis to the edge of the cover concrete -- outer edge of cover concrete
set coreY [expr $coverY-$coverH]; # The distance from the section z-axis to the edge of the core concrete -- edge of the core concrete/inner edge of cover concrete
set coreZ [expr $coverZ-$coverB]; # The distance from the section y-axis to the edge of the core concrete -- edge of the core concrete/inner edge of cover concrete
set numBarsInt [expr $numBarsIntTot/2]; # number of intermediate bars per side

set coreID 2567
set coverID 2568

# Define the fiber section
section fiberSec $flexSec {
# Define the core patch
patch quadr $coreID $nfCoreZ $nfCoreY -$coreY $coreZ -$coreY -$coreZ $coreY -$coreZ $coreY $coreZ

# Define the four cover patches
patch quadr $coverID 2 $nfCoverY -$coverY $coverZ -$coreY $coreZ $coreY $coreZ $coverY $coverZ
patch quadr $coverID 2 $nfCoverY -$coreY -$coreZ -$coverY -$coverZ $coverY -$coverZ $coreY -$coreZ
patch quadr $coverID $nfCoverZ 2 -$coverY $coverZ -$coverY -$coverZ -$coreY -$coreZ -$coreY $coreZ
patch quadr $coverID $nfCoverZ 2 $coreY $coreZ $coreY -$coreZ $coverY -$coverZ $coverY $coverZ

# define reinforcing layers
layer straight $steelID $numBarsInt $barAreaInt -$coreY $coreZ $coreY $coreZ; # intermediate skin reinf. (+z parallele ad y sulla faccia z+)
layer straight $steelID $numBarsInt $barAreaInt -$coreY -$coreZ $coreY -$coreZ; # intermediate skin reinf. -z
layer straight $steelID $numBarsTop $barAreaTop $coreY $coreZ $coreY -$coreZ; # top layer reinfocement
layer straight $steelID $numBarsBot $barAreaBot -$coreY $coreZ -$coreY -$coreZ; # bottom layer reinforcement

}; # end of fibersection definition



#LIMIT STATE MATERIALS


# Create axial failure spring
#source CenterColAxialSpring.tcl
limitCurve Axial $axialCurveTag $bcTag 21.27 24.7 3 2 2 2 4 1 2 0 0
uniaxialMaterial LimitState $axialFailTag\
25 0.000029 50 0.0000579 75 0.0000869 -25 -0.000029 -50 -0.0000579 -75 -0.0000869 0.5 0.4 0 0 0 $axialCurveTag 1

# Create shear failure spring
#source CenterColShearSpring.tcl
limitCurve Shear $shearCurveTag $bcTag\
0.001571429 5000 20 14 11.25 21.27 24.7 3 2 0 2 4 1 2 0
uniaxialMaterial LimitState $shearTag\
25 0.014706 30 0.017647 45 0.026471 -25 -0.0147059 -30 -0.017647 -45 -0.026471 0.5 0.4 0 0 0 $shearCurveTag 2 0

section Aggregator $shearAxialOnlySec $shearTag Vy $axialFailTag P;



#################################
# Define the beam-column element
#################################
geomTransf Linear 1
set nint 5; #che cos'é?
element forceBeamColumn $bcTag 2 3 $nint $flexSec 1 -iter 5 1e-15
####################################
# Define the zero-length end springs
####################################
uniaxialMaterial Elastic $rigidMatTag 9.9e9
element zeroLength 3 3 4 -mat $shearTag $axialFailTag $rigidMatTag -dir 1 2 6
#soft spring to take axial load after failure
uniaxialMaterial Elastic $softMatTag 99.9
element zeroLength 4 4 5 -mat $softMatTag -dir 2
#####################
# Recorders
#####################
recorder Node -file PushoverCiclica/node4DispX.out -time -node 4 -dof 1 disp
recorder Node -file PushoverCiclica/nodesReactY.out -time -node 2 -dof 2 reaction
recorder Element -file PushoverCiclica/Column.out -time -ele $bcTag force; #mio $bcTag

set P -79.0
pattern Plain 1 Constant {
load 4 0.0 $P 0.0
}
set loading cyclic
set control displ;
#######################
# Apply gravity loads #
#######################
initialize
system BandGeneral
constraints Penalty 1.0e12 1.0e12
numberer RCM
test NormDispIncr 1.0e-9 50 0
algorithm ModifiedNewton
integrator LoadControl 0 1 0 0
analysis Static
analyze 10

loadConst -time 0.0

#################################################
# Apply transverse displacements or timehistory #
#################################################

if {$iterazione==10} {
set nSteps_0 [list 788 1576 2000 2000 3000 3000 1 1 1]}

set count 30001;
set stepp -1;

foreach du {0.005 -0.005 0.005 -0.005 +0.005 -0.005 +0.005 -0.005 0.005} {
incr stepp 1;
incr count 1;

set nSteps [lindex $nSteps_0 $stepp]
pattern Plain $count "Linear -factor $du" {
sp 4 1 -1.0}
###############################################
# Set Analysis options for transverse loading #
###############################################
wipeAnalysis

variable numbererType Plain ;#RCM
numberer $numbererType

variable systemType BandGeneral;
system $systemType;

if {$control == "load"} {
integrator LoadControl $dlamda 1 $dlamda $dlamda
constraints Penalty 1.0e14 1.0e14

} elseif {$control == "displ"} {
integrator DisplacementControl 4 1 $du 1 $du $du
constraints Penalty 1.0e14 1.0e14

} elseif {$control == "newmark"} {

integrator Newmark 0.5 0.25 ;# no need to consider damping here as considered in Rayleigh command
constraints Transformation

} else {
puts stderr "Invalid control option: $control"
}

logFile logfile.txt

variable Tol 1.e-8;
variable maxNumIter 80;
variable printFlag 0;
variable testType NormDispIncr;
test $testType $Tol $maxNumIter $printFlag;

variable algorithmType Newton;
algorithm $algorithmType;
analysis Static
#########################################
# Analyze model with transverse loading
#########################################
set ok 0 ;# ok=0 for converged analysis step
set n 1
while {$n < $nSteps && $ok == 0} {
set ok [analyze 1 [expr 0.01*$dlamda] ] ;# analyse 1 step with original settings for algorithm and test

if {$ok != 0} {
#blah
test NormDispIncr 1.0e-10 50000 5 ;# increase max number of iterations
algorithm ModifiedNewton -initial ;# use initial elastic stiffness for NewtonRaphson
set ok [analyze 1 [expr 0.01*$dlamda] ] ;# analyse 1 step with new settings
if {$ok != 0} {
# #blah
test NormDispIncr 1.0e-6 30000 1 ;# increase max number of iterations
algorithm KrylovNewton -initial ;# use initial elastic stiffness for NewtonRaphson
set ok [analyze 1 [expr 0.01*$dlamda] ] ;# analyse 1 step with new settings
}

if {$ok != 0} {
puts "Trying Broyden .."
algorithm Broyden 15
set ok [analyze 1 [expr 0.01*$dlamda]]
}

if {$ok != 0} {
puts "Trying NewtonWithLineSearch .."
algorithm NewtonLineSearch .8
set ok [analyze 1 [expr 0.01*$dlamda]]
}

algorithm Newton ;# restore algorithm to Newton with current stiffness
test NormDispIncr 1.0e-10 80 0 ;# restore max number of iterations
}
incr n
}


if {$ok != 0} {

puts "ANALYSIS FAILED"
} else {

puts "SUCCESSFULL ANALYSIS"
}
}

[fmk]

you need to modify the script so i can see what it is you expect versus what OpenSees is giving you so i can figure out where you are messing up .. even then you might need to come up with a simpler script .. look at it from our point of view, we don't want to try and figure out what it is you are doing and then figure out where you messed up in a script containing over 300 lines of code!!

eg. instead of puts "ANALYSIS FAILED"
puts "RESULTS NOT RIGHT BECAUSE $A is not equal $B"