Moment Curvature Analysis-Fail to converge

[mairead]

Hi there,

I am trying to do a moment curvature analysis of a RC circular cross section but it is failing to converge. I was wondering if you could tell me what I am doing wrong? I have copied my code below(It is based on the mom-curv example on the website).

Thanks


# Define model builder
# --------------------
model basic -ndm 2 -ndf 3

# Define materials for nonlinear columns
# ------------------------------------------
#Units N, m
set PI 3.14;

#UNCONFINED CONCRETE
set fco 30.e6; #Concrete Strength
set Eco [expr 5000000.*pow($fco,.5)]; #Elastic modulus in N/m^2

#STEEL REINFORCEMENT
set fy 460.e6; #Steel Strength
set fyh 460.e6; #Transverse Steel Strength
set Es 200.e9; #modulus of steel N/m^2
set esu .12; #Longitudinal steel maximum strain
set fsu [expr 1.5*$fy]; #Longitudinal steel ultimate stress N/m^2


#Calculate confined concrete properties using manders model
set Dh .016; #Diameter of transverse reinforcement in m
set s .125; #Spacing of transverse reinforcement in m

set Dcol 1.0; #Columm Diameter
set clb .035; #Concrete Cover
set Dcore [expr $Dcol-(2*($clb+$Dh))]; #Core Diameter
set Ash [expr ($PI*pow($Dh,2))/4]; #Transverse Bar Area
set Ros [expr (4*$Ash)/($Dcore*$s)];
set fl [expr (2*$fyh*$Ash)/($Dcore*$s)]; #Maximum effective lateral pressure in N/m^2
set fll [expr .95*$fl]; #Effective lateral confining stress in N/m^2

set fcc [expr $fco*((2.254*(pow((1+((7.94*$fll)/$fco)),.5)))-((2*$fll)/$fco)-1.254)]; #Confined concrete compressive strength
set ecc [expr .002*(1+(5*(($fcc/$fco)-1)))]; #Confined concrete strain
set Esec [expr $fcc/$ecc]; #Secant Modulus of elasticity
set r [expr $Eco/($Eco-$Esec)];
set ecu [expr .004+((1.4*$Ros*$fyh*$esu)/$fcc)]; #Ultimate compression strain in concrete
set x [expr $ecu/$ecc];
set fcu [expr ($fcc*$x*$r)/($r-1+pow($x,$r))]; #Ultimate compressive strength in concrete


#MATERIAL GENERATION

set IDConcCoreCol 1;
set IDConcCovCol 2;
set IDReinfCol 3;
#CONCRETE tag f'c ec0 f'cu ecu
#Core concrete (confined)
uniaxialMaterial Concrete01 $IDConcCoreCol [expr -$fcc] [expr -$ecc] [expr -$fcu] [expr -$ecu];

#Cover concrete (un-confined)
uniaxialMaterial Concrete01 $IDConcCovCol [expr -$fco] [expr -(2.0*$fco/$Eco)] -10. -0.00550;

#REINFORCING STEEL Tag fy E b
uniaxialMaterial Steel01 $IDReinfCol $fyh $Es 0.018;



# Define cross-section for nonlinear columns
# ------------------------------------------

set ColSecTag 1;

set numBars 18; # Number of Longitudinal Reinforcing Bars
set barArea .001256; # Area of 1 reinforcing bar
set theta [expr 360.0/$numBars]; # Determine angle increment between bars
set Rcol [expr $Dcol/2]; # Column Radius
set Rcore [expr $Dcore/2]; # Core Radius

set nfCoreR 8; # number of radial divisions in the core (number of "rings")
set nfCoreT 12; # number of theta divisions in the core (number of "wedges")
set nfCoverR 2; # number of radial divisions in the cover
set nfCoverT 12; # number of theta divisions in the cover

section fiberSec $ColSecTag {
# tag div raddiv Y-cen Z-cen int-rad out-rad start end
#CORE PATCH
patch circ $IDConcCoreCol $nfCoreT $nfCoreR 0 0 0.0 $Rcore 0. 360.;
#COVER PATCH
patch circ $IDConcCovCol $nfCoverT $nfCoverR 0 0 $Rcore $Rcol 0. 360.;
#REINFORCING LAYER
# tag numBars areaBar Y-cen Z-cen rad start end
layer circ $IDReinfCol $numBars $barArea 0 0 $Rcore $theta 360.;

}

# Estimate yield curvature
# (Assuming no axial load and only top and bottom steel)
set d [expr $Dcol-$clb] ;# d -- from cover to rebar
set epsy [expr $fy/$Es] ;# steel yield strain
set Ky [expr $epsy/(0.7*$d)]

# Print estimate to standard output
puts "Estimated yield curvature: $Ky"

# Set axial load
set P -1.8e6;

set mu 15; # Target ductility for analysis
set numIncr 100; # Number of analysis increments

# Call the section analysis procedure
source MomentCurvature.tcl
MomentCurvature 1 $P [expr $Ky*$mu] $numIncr

----------------------------------------------------------------------------------------------------------------------

proc MomentCurvature {secTag axialLoad maxK {numIncr 100} } {
# Define two nodes at (0,0)
node 1 0.0 0.0
node 2 0.0 0.0

# Fix all degrees of freedom except axial and bending
fix 1 1 1 1
fix 2 0 1 0

# Define element
# tag ndI ndJ secTag
element zeroLengthSection 1 1 2 $secTag

# Create recorder
recorder Node -file section$secTag.out -time -node 2 -dof 3 disp

# Define constant axial load
pattern Plain 1 "Constant" {
load 2 $axialLoad 0.0 0.0
}

# Define analysis parameters
integrator LoadControl 0.0
system SparseGeneral -piv; # Overkill, but may need the pivoting!
test NormUnbalance 1.0e-9 10
numberer Plain
constraints Plain
algorithm Newton
analysis Static

# Do one analysis for constant axial load
analyze 1

# Define reference moment
pattern Plain 2 "Linear" {
load 2 0.0 0.0 1.0
}

# Compute curvature increment
set dK [expr $maxK/$numIncr]

# Use displacement control at node 2 for section analysis
integrator DisplacementControl 2 3 $dK 1 $dK $dK

# Do the section analysis
analyze $numIncr

[vesna]

Try to do the analysis in the loop. Here is an example: http://opensees.berkeley.edu/wiki/index ... e_Analysis

[mairead]

Thank you!