ultimate curvature

[sisa]

At first happy Easter to everybody!
I have a problem with the result of Moment-Curvature Analysis. If I consider a RC section (for simplicity all the section is of confined concrete) the result of ultimate curvature it’s very different to the result I get with another program (Response 2000). In fact with Opensess the ultimate curvature (that I want to be at ultimate compressive strain for concrete that is in this case 0,012) is very larger than one of Response2000 (6.8e-5 rad/mm in Opensees and 3.49e-5 rad/mm in Response2000) while the yield curvature (that correspond to reinforcement strain of 0,002) are very similar. I don’t understand what is the problem because the section and the materials would be the same. In order to read the ultimate curvauture I read in the recorder of stressstrain fiber (Cls.out) the time when strain reaches the 0,012 and in the Moment-Curvature recorder (MC.out) the corrisponding curvature.Is it the right procedure?
I hope someone can help me.
Thank you
Annalisa
These are the script (one for the section and one for Moment-Curvature analysis)

section
#MATERIALE
set IDarmature 2
set IDclsConf 3
set Ec [expr 26500]
set fcc -37.09
set epsc1c -0.0052
set fcuc -33.28
set epscuc -0.012

set fy 450
set fyd [expr $fy/1.15]
set Es [expr 2.1e5]
set epsyd 0.002
set b 0.005
uniaxialMaterial Concrete04 $IDclsConf $fcc $epsc1c $epscuc $Ec;
uniaxialMaterial Steel01 $IDarmature $fy $Es $b

#GEOMETRIA
set H 300
set B 300
set cf 30
set numArm 4
set diamArm 20
set areaArm [expr 3.14*$diamArm**2/4]
set SezTag 1

#FIBRE
# y orizzontale z verticale
set fibreyCore 10
set fibrezCore 10
set fibreyCover 10
set fibrezCover 10
set coverZ [expr $H/2]
set coverY [expr $B/2]
set coreY [expr $coverY-$cf]
set coreZ [expr $coverZ-$cf]

section Fiber $SezTag {
patch quadr $IDclsConf $fibrezCore $fibreyCore -$coverY -$coverZ $coverY -$coverZ $coverY $coverZ -$coverY $coverZ;
layer straight $IDarmature $numArm $areaArm -$coreY $coreZ -$coreY -$coreZ;# ;#armatura superiore
layer straight $IDarmature $numArm $areaArm $coreY $coreZ $coreY -$coreZ;#; #armatura inferiore
}

source MomentCurvature2D.tcl

# set AXIAL LOAD --------------------------------------------------------
set P 0; # + Tension, - Compression

# set maximum Curvature:
set Ku [expr 7e-4];#curvatura ultima giusta
set numIncr 100; # Number of analysis increments to maximum curvature (default=100)
# Call the section analysis procedure
MomentCurvature2D $SezTag $P $Ku $numIncr $IDclsConf

Moment Curvature analysis
proc MomentCurvature2D { SezTag axialLoad maxK numIncr mat} {

# Define two nodes at (0,0)
node 1001 0.0 0.0
node 1002 0.0 0.0

# Fix all degrees of freedom except axial and bending
fix 1001 1 1 1
fix 1002 0 1 0;#lascio libera la deformazione assiale (per consentire la deformazione del materiale) e la rotazione

# Define element
# tag ndI ndJ secTag
element zeroLengthSection 2001 1001 1002 $SezTag

# Create recorder
recorder Node -file SezioneCA/Momento.out -time -node 1003 -dof 3 reaction; # stampa il valore del momento
recorder Node -file SezioneCA/Curvatura.out -node 1002 -dof 3 disp; #stampa il valore della curvatura
recorder Node -file SezioneCA/MC.out -time -node 1002 -dof 3 disp; #stampa momento e curvatura in 2 colonne
recorder Element -file SezioneCA/armaturaDef.out -ele 2001 section fiber -120 120 strain
recorder Element -file SezioneCA/armaturaTens.out -time -ele 2001 section fiber -120 120 stress
recorder Element -file SezioneCA/Cls.out -time -ele 2001 section fiber 150 12 stressStrain

# Define constant axial load
pattern Plain 3001 "Constant" {
load 1002 $axialLoad 0.0 0.0;# controllare che il carico assiale sia verticale
}

# Define analysis parameters
integrator LoadControl 0 1 0 0
system SparseGeneral -piv
test EnergyIncr 1.0e-9 10
numberer Plain
constraints Plain
algorithm Newton
analysis Static

# Do one analysis for constant axial load
analyze 1

loadConst -time 0.0

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

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

# Use displacement control at node 1002 for section analysis, dof 3
integrator DisplacementControl 1002 3 $dK 1 $dK $dK

# Do the section analysis
set ok [analyze $numIncr]
if {$ok == 0} {puts "convergenza raggiunta al primo tentativo"}

# ----------------------------------------------if convergence failure-------------------------
set IDctrlNode 1002
set IDctrlDOF 3
set Dmax $maxK
set Dincr $dK
set TolStatic 1.e-9;
set testTypeStatic EnergyIncr
set maxNumIterStatic 6
set algorithmTypeStatic Newton
if {$ok != 0} {
# if analysis fails, we try some other stuff, performance is slower inside this loop
set Dstep 0.0;
set ok 0
while {$Dstep <= 1.0 && $ok == 0} {
set controlDisp [nodeDisp $IDctrlNode $IDctrlDOF ]
set Dstep [expr $controlDisp/$Dmax]
set ok [analyze 1]; # this will return zero if no convergence problems were encountered
if {$ok != 0} {; # reduce step size if still fails to converge
set Nk 4; # reduce step size
set DincrReduced [expr $Dincr/$Nk];
integrator DisplacementControl $IDctrlNode $IDctrlDOF $DincrReduced
for {set ik 1} {$ik <=$Nk} {incr ik 1} {
set ok [analyze 1]; # this will return zero if no convergence problems were encountered
if {$ok != 0} {
# if analysis fails, we try some other stuff
# performance is slower inside this loop global maxNumIterStatic; # max no. of iterations performed before "failure to converge" is ret'd
puts "Trying Newton with Initial Tangent .."
test NormDispIncr $TolStatic 2000 0
algorithm Newton -initial
set ok [analyze 1]
test $testTypeStatic $TolStatic $maxNumIterStatic 0
algorithm $algorithmTypeStatic
}
if {$ok != 0} {
puts "Trying Broyden .."
algorithm Broyden 8
set ok [analyze 1 ]
algorithm $algorithmTypeStatic
}
if {$ok != 0} {
puts "Trying NewtonWithLineSearch .."
algorithm NewtonLineSearch 0.8
set ok [analyze 1]
algorithm $algorithmTypeStatic
}
if {$ok != 0} {; # stop if still fails to converge
puts "PROBLEM $IDctrlNode $IDctrlDOF [nodeDisp $IDctrlNode $IDctrlDOF]"
return -1
}; # end if
}; # end for
integrator DisplacementControl $IDctrlNode $IDctrlDOF $Dincr; # bring back to original increment
}; # end if
}; # end while loop
}; # end if ok !0
# -----------------------------------------------------------------------------------------------------


if {$ok != 0 } {
puts "PROBLEM"
} else {
puts "DONE"
}

}

[andrett7]

Hi Annalisa,

about the difference I don't know really because I don't use Response 2000. May be near the ultimate value the behavior of the materials is sligthly different.
About the procedure I think is correct.

Buona Pasqua anche a te.

Andrea

[sisa]

Grazie Andrea per il sostegno!
Can someone check this program in order to see if there's something wrong? I really don't know why the result are so different from those I get whit manual calculation or with other program...I am considering a RC section with 2 layer each of 4 bar with diameter 20mm. The problem is that concrete reach ultimate stress (0.012) with a moment-curvature very larger than those I need with manual calculation...
Thank you
Annalisa

[sisa]

I'm controlling the same script and I recognize that the stressStrain of concrete and of reinforcement bar don't respect the hypothesis of plane section. Infact, given a Curvature and a strain of bar, if I calculate the position of neutral axis and the strain of concrete (with a linear strain profile) they are different from those recorded by OpenSess. Is It possible? Is It right or there is a problem in the program?
please help me because I really dont' found where is the problem....

[vesna]

I looked at your code and it looks good to me. The way to find ultimate curvature is to look at moment-curvature response. After you see sudden reduction of moment capacity, you can consider that to be the ultimate curvature.

Ultimate curvature is sensitive to the definition of your materials in the post-yielding range. When you say your ultimate value does not match the manually calculated value, what do you mean by that? What do you mean by manual calculation?

The fiber section is developed under assumption of Bernoulli's theorem. So the section stays plain after deformation.

One more thing. I noticed that you define dimensions of your elements as integers. It is always better to define it as a floating point values. Also, it is better to have denominators as a floating point values.

[sisa]

Thank you Vesna for your answer.
My work is to compare result of Opensees with results I optain with my simplified manual calculation (I calculate Moment and Curvature in corrispondance to yeld of bar and to ultimate state of section,that I consider be in corrispondance to ultimate compression strain of concrete). So when I read the recorder of Opensees and I read for which Moment-Curvature concrete reaches its ultimate compression strain, the two result are very different. And the other problem is that if I take the recorders of reinforcement bar's strain (εs) and corresponding Curvature (C) I can calculate the position of neutral axis (x=(C*d-εs)/C, where d=H-cover) and so the strain of concrete εc (εc = C*x ), but this results are very different from the recorder of concrete strain made by Opensees (in file Cls.out). For this reason I don't know if section stays plain.
For the fact of floating point values do you intend I have to write in this way?
set H 300.0
set B 300.0
set cf 30.0
thank you
Annalisa

[andrett7]

Hi Annalisa,

about floating point values yes, it you should use 300.0 and so on, otherwise there is a huge approximation.

For example:

expr (1/2) = 0
expr(1./2)=0.5

Bye

[vesna]

You can never match manually calculated strain values and strains calculated using very sophisticated models as fiber sections is.

[mairead]

Hi everyone,

I too am trying to find the ultimate curvature. However my the moment does not seem to reduce at any point, it just seems to keep increasing with curvature? I'm wondering is there a problem with my code? Maybe the way I define my materials? If you could take a look I would really appreciate it. Thanks.

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


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

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

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


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

set clb .035;
set Dcol 1.0; #Columm Diameter
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


#CONCRETE tag f'c ec0 f'cu ecu
#Core concrete (confined)
uniaxialMaterial Concrete01 1 -$fcc -$ecc -$fcu -$ecu;

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

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



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


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 Fiber 1 {
# tag div raddiv Y-cen Z-cen int-rad out-rad start end
#CORE PATCH
patch circ 1 $nfCoreT $nfCoreR 0.0 0.0 0.0 $Rcore 0. 360.;
#COVER PATCH
patch circ 2 $nfCoverT $nfCoverR 0.0 0.0 $Rcore $Rcol 0. 360.;
#REINFORCING LAYER
# tag numBars areaBar Y-cen Z-cen rad start end
layer circ 3 $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 -1800000;

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

# Call the section analysis procedure
source MomentCurvatureDeterministic.tcl
MomentCurvatureDeterministic 1 $P [expr $Ky*$mu] $numIncr
puts "analysis done";

#===========================================================
# Moment Curvature Procedure
#===========================================================
proc MomentCurvatureDeterministic {secTag axialLoad maxK {numIncr 100}} {
set name MomentCurvatureDeterministic;
set name1 Deteriorated
file mkdir $name/$name1; # create data directory

# 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 [concat $name/$name1/MomCurve.out] -time -node 2 -dof 3 disp
recorder Element -file [concat $name/$name1/MomStrain.out] -time -ele 1 section fiber -0.5 0.0 stressStrain

# 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.e-8 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