T-Beam Moment Curvature Analysis

Forum for OpenSees users to post questions, comments, etc. on the use of the OpenSees interpreter, OpenSees.exe

Moderators: silvia, selimgunay, Moderators

Post Reply
stuarts
Posts: 6
Joined: Sat Sep 27, 2008 8:00 pm
Location: UW

T-Beam Moment Curvature Analysis

Post by stuarts »

I am trying to run a Moment Curvature analysis on a reinforced concrete T-Beam section but when I run my code it returns what seems to be a convergence error. I have been trying to figure it out but have thus far been unsuccessful. Below is the code that I have been using:

Code: Select all

puts " -------------------------2D Model------------------------"

puts " -------------------------------------------------------------- "

puts "                           RC T-Beam Section"

puts " -------------------------------------------------------------- "

# Define model builder
# --------------------
model basic -ndm 2 -ndf 3
set secTag 1

# define UNITS ----------------------------------------------------------------------------
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
set pcf [expr $lbf/pow($ft,3)];		# pounds per cubic foot
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

# set up parameters for column section and element definition
#--------------------------------------------------
set IDcore 1; 		# ID tag for core concrete
set IDcover 2; 		# ID tag for cover concrete
set IDsteel 3; 		# ID tag for steel

# Define materials for nonlinear columns
# ------------------------------------------

# Confined concrete: (CORE CONCRETE)
set fc [expr -4*$ksi]; 			# CONCRETE Compressive Strength, ksi (+Tension, -Compression)
set Ec [expr 57*$ksi*sqrt(-$fc/$psi)]; 	# Concrete Elastic Modulus
set fc1C [expr 1.2649*$fc]; 		# CONFINED concrete (mander model), maximum stress
set eps1C [expr 2.*$fc1C/$Ec]; 		# strain at maximum stress
set fc2C $fc; 			# ultimate stress
set eps2C [expr 5*$eps1C]; 		# strain at ultimate stress

uniaxialMaterial Concrete01 $IDcore $fc1C $eps1C $fc2C $eps2C;

# Unconfined concrete: (COVER CONCRETE)
set fc1U $fc; 			# UNCONFINED concrete (todeschini parabolic model), maximum stress
set eps1U -0.003; 			# strain at maximum stress
set fc2U [expr 0.1*$fc]; 		# ultimate stress
set eps2U -0.006; 			# strain at ultimate stress

uniaxialMaterial Concrete01 $IDcover $fc1U $eps1U $fc2U $eps2U; 

# STEEL
# Reinforcing steel 
set fy 60.0;      		# Yield stress
set E 29000.0;    		# Young's modulus
#                                tag  fy E0    b
uniaxialMaterial Steel01  3  $fy $E 0.01

# Define cross-section for T-Beam section
# ------------------------------------------

# set some paramaters
set bf [expr 24*$in];		# flange width (For a rect. beam this is the width or 0)
set d [expr 18*$in];		# nominal depth
set bw [expr 6*$in];		# web thickness (For a rect beam this is the width b)
set tf [expr 4*$in];		# flange thickness (For a rectangular beam set this equal to 0)

set nfdw 10;		# number of fibers along web depth 
set nftw 10;		# number of fibers along web thickness
set nfbf 10;		# number of fibers along flange width (you want this many in a bi-directional loading)
set nftf 10;		# number of fibers along flange thickness

set cover  [expr 1.5*$in];
set As [expr 1.0*$in*$in];    	 # area of no. 9 bars
puts "As = $As"

# some variables derived from the parameters

set y1 $tf	
set y2 [expr $d-$y1]		 
set z1 [expr $bf/2]
set z2 [expr $bw/2]

# Discretize the section with layers and patches
section Fiber $secTag {
     	
	# Define Web Section
	#             nfIJ  nfJK       yI    zI    yJ      zJ  yK  zK  yL   zL
     	patch quad  $IDcore $nftw $nfdw  -$y2 $z2 -$y2 -$z2 0 -$z2 0 $z2
	# Define Flange Section
	#               nfIJ  nfJK    yI   zI  yJ   zJ    yK    zK     yL    zL
     	patch quad  $IDcore $nfbf $nftf       0 $z2 0 -$z1 $y1 -$z1 $y1 $z1

	# Define Reinforcing Steel
	# 	[y start]  [z start] [y end]  [z end]
	#layer straight $IDsteel 3 $As [expr -$y2+$cover] [expr $z2-$cover] [expr -$y2t+$cover] [expr -$z2+$cover]
  }

# Estimate yield curvature
# (Assuming no axial load and only top and bottom steel)
set d [expr $d-$cover];	# d -- from cover to rebar
set epsy [expr $fy/$E];	# 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 [expr 0*$kip];		# + Tension, - Compression

# set maximum Curvature:
set Ku [expr 0.01/$in];
set numIncr 100;		# Number of analysis increments to maximum curvature (default=100)
puts "numIncr: $numIncr"

# ANALYSIS of section --------------------------------------------------

	# 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

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

	# Create recorder
	recorder Node -file data/Mphi.out -time -node 1002 -dof 3 disp;	# output moment (col 1) & curvature (col 2)
	
	# Define constant axial load
	pattern Plain 3001 "Constant" {
		load 1002 $P 0.0 0.0
	}

	# Define analysis parameters
	integrator LoadControl 0 1 0 0
	system SparseGeneral -piv;	# Overkill, but may need the pivoting!
	test EnergyIncr  1.0e-2 100
	numberer Plain
	constraints Plain
	algorithm Newton
	analysis Static

	# Do one analysis for constant axial load
	analyze 1

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

	# Compute curvature increment
	set dK [expr $Ku/$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 convergence failure-------------------------
	set IDctrlNode 1002
	set IDctrlDOF 3
	set Dmax $Ku
	set Dincr $dK
	set TolStatic 1.e-2;
	set testTypeStatic EnergyIncr  
	set maxNumIterStatic 6
	set algorithmTypeStatic Newton
	set fmt1 "%s Pushover analysis: CtrlNode %.3i, dof %.1i, Curv=%.4f /%s";	# format for screen/file output of DONE/PROBLEM 

analysis
	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 [format $fmt1 "PROBLEM" $IDctrlNode $IDctrlDOF [nodeDisp $IDctrlNode 

$IDctrlDOF] $LunitTXT]
						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
	# -----------------------------------------------------------------------------------------------------
	global LunitTXT;							# load time-unit text
	if {  [info exists LunitTXT] != 1} {set LunitTXT "Length"};			# set blank if it has not been defined previously.
	
	if {$ok != 0 } {
		puts [format $fmt1 "PROBLEM" $IDctrlNode $IDctrlDOF [nodeDisp $IDctrlNode $IDctrlDOF] $LunitTXT]
	} else {
		puts [format $fmt1 "DONE"  $IDctrlNode $IDctrlDOF [nodeDisp $IDctrlNode $IDctrlDOF] $LunitTXT]
	}

puts "Analysis Finished"
puts " -------------------------------------------------------------- "
Thanks for the help.
zvidrih
Posts: 39
Joined: Wed Apr 30, 2008 1:55 am
Location: Ljubljana, Slovenia

error

Post by zvidrih »

You have an

Code: Select all

"analysis" 
command in line 185 (just before if {$ok != 0} {...) but this could be just the because of the limited page width at the forum.

And you have also commented the reinforcement layer (BTW - "y2t" is not defined)

Code: Select all

 #layer straight $IDsteel 3 $As [expr -$y2+$cover] [expr $z2-$cover] [expr -$y2t+$cover] [expr -$z2+$cover] 
If that is intentional, then your task of performing the moment-curvature analysis for unreinforced concrete section (concrete 01 has no tension) without any axial load is pretty strange.

Hope it helps
Zlatko Vidrih
Institute of Structural Engineering, earthquake Engineering and Construction IT
Faculty of Civil and Geodetic Engineering
University of Ljubljana, Slovenia
stuarts
Posts: 6
Joined: Sat Sep 27, 2008 8:00 pm
Location: UW

Post by stuarts »

Thanks for the help. I fixed some of the things you pointed out then it ran perfectly. I guess I had been looking at the code for too long and couldn't see the obvious. Cheers.
Post Reply