Moment curvature

[yfeng]

Hi all,
I am doing moment curvature analysis. The analysis can not converge when I use stress unit MPa, force unit N and length unit mm( N/mm^2=Mpa). But I change the unit to kN/m^2, kN and m, the analysis worked well. There must be something I miss when changing the units. Can anybody find out where is the problem. Thank you so much.

Here is my code.

# SET UP ----------------------------------------------------------------------------
wipe; # clear memory of all past model definitions
model BasicBuilder -ndm 2 -ndf 3; # Define the model builder, ndm=#dimension, ndf=#dofs
set dataDir DataOriModel_MomentCurvature2; # set up name of data directory (can remove this)
file mkdir $dataDir; # create data directory
set GMdir "../GMfiles/"; # ground-motion file directory


source RCcircSection2.tcl; # procedure for definining RC fiber section

# Define SECTIONS -------------------------------------------------------------
set SectionType FiberSection; # options: Elastic FiberSection

# define section tags:
set ColSecTag 1

# Definition of parameters for fiber RC section

# MATERIAL parameters

# confined and unconfined CONCRETE-----------------------------------------------------------------
# nominal concrete compressive strength

set fc -28; # CONCRETE Compressive Strength, MPa (+Tension, -Compression)
set Ec [expr 5000*sqrt(-$fc)]; # Concrete Elastic Modulus,MPa


# confined concrete
set Kfcb 1.46; # ratio of confined to unconfined concrete strength
set Kresb 0.6; # ratio of residual/ultimate to maximum stress

set fc1C [expr $Kfcb*$fc]; # CONFINED concrete (mander model), maximum stress;MPa
# CONFINED concrete (mander model), maximum stress;MPa

set eps0 -0.002; # UNCONFINED concrete strain

set fac [expr $fc1C/$fc];


set eps1C [expr $eps0*(1+5*($fac-1))]; # strain at confined concrete strength

set fc2C [expr $Kresb*$fc1C]; # crushing strength;MPa

set eps2C -0.05;

# unconfined concrete---------------------------------------------------------------------------------

set fc1U $fc; # UNCONFINED concrete (todeschini parabolic model), maximum stress
set eps1U $eps0; # strain at maximum strength of unconfined concrete
set fc2U -0.001; # crushing strength;MPa
set eps2U -0.0064; # strain at ultimate stress


set IDconcCore 1;
set IDconcCover 2;
set IDconcCorec 3;

uniaxialMaterial Concrete01 $IDconcCore $fc1C $eps1C $fc2C $eps2C ; # Core concrete (confined)

uniaxialMaterial Concrete01 $IDconcCover $fc1U $eps1U $fc2U $eps2U; # Cover concrete (unconfined)

# REINFORCING STEEL parameters--------------------------------------------------------------------------
#
set Fy 450; # STEEL yield stress;1000kN/m^2=MPa
set Es 200000; # modulus of steel;MPa
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
set Fu 675;
set esh 0.008;
set Esh 8000;
set eult 0.12;

set IDreinf 4

# constants alpha = 0.506, Cf = 0.26, Cd = 0.45

set alpha 0.506;
set cf 0.26;
set cd 0.45;

uniaxialMaterial ReinforcingSteel $IDreinf $Fy $Fu $Es $Esh $esh $eult ;

# section GEOMETRY -------------------------------------------------------------
set pi 3.1415;
set DSec 609.6; # Column Diameter
set Across [expr ($pi/4)*pow($DSec,2)];
set coverSec 28.6; # Column cover to reinforcing steel NA.
set numBarsSec 16; # number of longitudinal-reinforcement bars in column
set dbar 19.05 ; # bar diameter
set barAreaSec [expr ($pi/4)*pow($dbar,2)]; # area of longitudinal-reinforcement bars
set ColSecTag 1; # set tag for symmetric section





# Generate a circular reinforced concrete section
# with one layer of steel evenly distributed around the perimeter and a confined core.
# confined core.
#
#
# Notes
# The center of the reinforcing bars are placed at the inner radius
# The core concrete ends at the inner radius (same as reinforcing bars)
# The reinforcing bars are all the same size
# The center of the section is at (0,0) in the local axis system
# Zero degrees is along section y-axis
#
set ri 0.0; # inner radius of the section, only for hollow sections
set ro [expr $DSec/2]; # overall (outer) radius of the section
set nfCoreR 8; # number of radial divisions in the core (number of "rings")
set nfCoreT 24; # number of theta divisions in the core (number of "wedges")
set nfCoverR 2; # number of radial divisions in the cover
set nfCoverT 24; # number of theta divisions in the cover



RCcircSection $ColSecTag $ri $ro $coverSec $IDconcCore $IDconcCover $IDreinf $numBarsSec $barAreaSec $nfCoreR $nfCoreT $nfCoverR $nfCoverT

set rc [expr $ro-$coverSec];

set IDColTransf 1; # all columns
set ColTransfType PDelta ; # options: Linear, PDelta, Corotational
geomTransf $ColTransfType $IDColTransf; # only columns can have PDelta effects (gravity effects)

#Estimate yield curvature
#(Assuming no aixal load)
set epsy [expr $Fy/$Es];
set Ky 0.0083;

puts "Estimated yield strain: $epsy"

puts "Estimated yield curvature: $Ky"

#set axial load

set P -805000; # force in N
set numIncr 100; #Number of analysis increments

#Call the section analysis procedure
source MomentCurvature.tcl
MomentCurvature 1 $P 0.25 $numIncr

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

[silvia]

i recommend you actually use the units as variables as i do in the examples, this will help you find the inconsistency.
fore example, the 5000 in front of the square root for E should have a unit factor....

[Macedo]

Hi,
I had the same problem a few weeks ago. So I decide to make a file defining the units. I use the following units definitions and the models works fine:

# Units_(SI).tcl

puts "Units Definition Started"
set m 1.; # Define Basic Unit - Length
set kg 1.; # Define Basic Unit - Mass
set sec 1.; # Define Basic Unit - Time
set rad 1.; # Define Basic Unit - Radians

set N [expr $kg*$m/pow($sec,2)]; # Define Engineering Units - Newton
set kN [expr $N*1000]; # Define Engineering Units - Kilo Newton
set Pa [expr $N/pow($m,2)]; # Define Engineering Units - Pascals
set MPa [expr $Pa*1000000]; # Define Engineering Units - Mega Pascal
set cm [expr $m/100]; # Define Engineering Units - Centimeter
set mm [expr $m/1000]; # Define Engineering Units - Milimeter

#Constants
set g [expr 9.81*$m/pow($sec,2)]; # Gravitational Acceleration
set PI [expr 2*asin(1.0)]; # PI constant 3.14159

I hope it helps...

[silvia]

exactly.
but you have to remember to ALWAYS use the unit variables throughout your input! it makes it soooo much better!!!

[silvia]

Macedo,
be careful here because you have not defined g in the conversion from mass to force -- this is tricky with the ambiguous kg....

[yfeng]

[quote="silvia"]i recommend you actually use the units as variables as i do in the examples, this will help you find the inconsistency.
fore example, the 5000 in front of the square root for E should have a unit factor....[/quote]


Thanks. I will do that.

[yfeng]

[quote="Macedo"]Hi,
I had the same problem a few weeks ago. So I decide to make a file defining the units. I use the following units definitions and the models works fine:

# Units_(SI).tcl

puts "Units Definition Started"
set m 1.; # Define Basic Unit - Length
set kg 1.; # Define Basic Unit - Mass
set sec 1.; # Define Basic Unit - Time
set rad 1.; # Define Basic Unit - Radians

set N [expr $kg*$m/pow($sec,2)]; # Define Engineering Units - Newton
set kN [expr $N*1000]; # Define Engineering Units - Kilo Newton
set Pa [expr $N/pow($m,2)]; # Define Engineering Units - Pascals
set MPa [expr $Pa*1000000]; # Define Engineering Units - Mega Pascal
set cm [expr $m/100]; # Define Engineering Units - Centimeter
set mm [expr $m/1000]; # Define Engineering Units - Milimeter

#Constants
set g [expr 9.81*$m/pow($sec,2)]; # Gravitational Acceleration
set PI [expr 2*asin(1.0)]; # PI constant 3.14159

I hope it helps...[/quote]

It helps. Thank you.

[andrewfirst]

I just encountered the same problem.
I modified the tcl file by using the unit system of SI: kg,m,sec,N...
The nonconvergence problem is solved.
But what is the reason? Numerical precision?

[silvia]

likely, you had the wrong scaling factor for masses and accelerations.
the units file keeps you honest and correct.

[andrewfirst]

The former units is : Ton, mm, sec, N, .... The result does not converge.
Then the unit is changed into SI: kg, m, sec, N... The result converge.
I think both two unit system is right. The problem may lie in the program itsself.

[silvia]

the program has no units!
you have some sort of conversion issue. also, you can only define 3 of the units, as only 3 are independent.