Gravity Analysis is failded

[yuqing]

Hi everyone
I am doing a pushover anaysis
As long as the gravity load is 0, everything will be OK.
However, when I apply the gravity load, the gravity analysis will not converge no matter how small it is !!!
Could someone give me some advice? THX!!


wipe all; # clear memory of past model definitions
model BasicBuilder -ndm 2 -ndf 3; # Define the model builder, ndm = #dimension, ndf = #dofs
source DisplayModel2D.tcl; # procedure for displaying a 2D perspective of model
source DisplayPlane.tcl; # procedure for displaying a plane in a model
source rotSpring2DModIKModel.tcl; # procedure for defining a rotational spring (zero-length element)

###################################################################################################
# Define Analysis Type
###################################################################################################
# Define type of analysis: "pushover" = pushover; "dynamic" = dynamic
set analysisType "pushover";

if {$analysisType == "pushover"} {
set dataDir Concentrated-Pushover-Output; # name of output folder
file mkdir $dataDir; # create output folder
}
if {$analysisType == "dynamic"} {
set dataDir Concentrated-Dynamic-Output; # name of output folder
file mkdir $dataDir; # create output folder
}


###################################################################################################
# Define Building Geometry, Nodes, and Constraints
###################################################################################################
# define structure-geometry parameters
set NStories 1; # number of stories
set NBays 1; # number of frame bays (excludes bay for P-delta column)
set WBay 5000 ; # bay width
set HStory1 3000; # 1st story height
set HBuilding 3000; # height of building

# calculate locations of beam/column joints:
set Pier1 0.0; # leftmost column line
set Pier2 [expr $Pier1 + $WBay];

set Floor1 0.0; # ground floor
set Floor2 [expr $Floor1 + $HStory1];


# calculate joint offset distance for beam plastic hinges
set phlat [expr 0.0]; # lateral dist from beam-col joint to loc of hinge on Floor 2


# calculate nodal masses -- lump floor masses at frame nodes
set g 9800; # acceleration due to gravity
set Floor2Weight 30161.42094 ; # weight of Floor 2 19994.29937
set NodalMass2 [expr (3.077696014)/ (4.0)]; # mass at each node on Floor 2
set Negligible 1e-9; # a very smnumber to avoid problems with zero

# define nodes and assign masses to beam-column intersections of frame
# command: node nodeID xcoord ycoord -mass mass_dof1 mass_dof2 mass_dof3
# nodeID convention: "xy" where x = Pier # and y = Floor #
node 11 $Pier1 $Floor1;
node 21 $Pier2 $Floor1;

node 1001 $Pier1 [expr $Floor1+1000] -mass [expr $NodalMass2 / (4.0)] $Negligible $Negligible;
node 1002 $Pier1 [expr $Floor1+2000] -mass [expr $NodalMass2 / (4.0)] $Negligible $Negligible;
node 1003 $Pier2 [expr $Floor1+1000] -mass [expr $NodalMass2 / (4.0)] $Negligible $Negligible;
node 1004 $Pier2 [expr $Floor1+2000] -mass [expr $NodalMass2 / (4.0)] $Negligible $Negligible;

node 12 $Pier1 $Floor2 -mass $NodalMass2 $Negligible $Negligible;
node 22 $Pier2 $Floor2 -mass $NodalMass2 $Negligible $Negligible;

node 2001 [expr $Pier1+1000] $Floor2 ;
node 2002 [expr $Pier1+2000] $Floor2 -mass [expr $NodalMass2 / (2.0)] $Negligible $Negligible;
node 2003 [expr $Pier1+3000] $Floor2 -mass [expr $NodalMass2 / (2.0)] $Negligible $Negligible;
node 2004 [expr $Pier1+4000] $Floor2 ;

# define extra nodes for plastic hinge rotational springs
# nodeID convention: "xya" where x = Pier #, y = Floor #, a = location relative to beam-column joint
# "a" convention: 2 = left; 3 = right;
# "a" convention: 6 = below; 7 = above;

# column hinges at bottom of Story 1 (base)
node 117 $Pier1 $Floor1;
node 217 $Pier2 $Floor1;

# column hinges at top of Story 1
node 126 $Pier1 $Floor2;
node 226 $Pier2 $Floor2;

# beam hinges at Floor 2
node 122 [expr $Pier1 + $phlat] $Floor2;
node 223 [expr $Pier2 - $phlat] $Floor2;



# constrain beam-column joints in a floor to have the same lateral displacement using the "equalDOF" command
# command: equalDOF $MasterNodeID $SlaveNodeID $dof1 $dof2...
#set dof1 1; # constrain movement in dof 1 (x-direction)
#equalDOF 12 22 $dof1; # Floor 2: Pier 1 to Pier 2
set dof1 1; # constrain movement in dof 1 (x-direction)
equalDOF 12 22 $dof1; # Floor 2: Pier 1 to Pier 2

# assign boundary condidtions
# command: fix nodeID dxFixity dyFixity rzFixity
# fixity values: 1 = constrained; 0 = unconstrained
# fix the base of the building; pin P-delta column at base
fix 11 1 1 1;
fix 21 1 1 1;


###################################################################################################
# Define Section Properties and Elements
###################################################################################################
# define material properties
set Es 2.01e5; # steel Young's modulus

# define column section W14x257 for Story 1
set Acol_ex 48774.096; # cross-sectional area
set Icol_ex 1409904188; # moment of inertia
set Mycol_ex 1694384582; # yield moment

# define beam section 改W30x116 （850x250x16x14）for Floor 2
set Abeam_2 19387.43958; # cross-sectional area
set Ibeam_2 2019046661; # moment of inertia
set Mybeam_2 1187674506; # yield moment at plastic hinge location



# determine stiffness modifications to equate the stiffness of the spring-elastic element-spring subassembly to the stiffness of the actual frame member

set n 1000; # stiffness multiplier for rotational spring

# calculate modified moment of inertia for elastic elements
set Icol_exmod [expr $Icol_ex*($n+1.0)/$n]; # modified moment of inertia for columns exterior
set Ibeam_2mod [expr $Ibeam_2*($n+1.0)/$n]; # modified moment of inertia for beams in Floor 2


# calculate modified rotational stiffness for plastic hinge springs
set Ks_col_ex [expr $n*6.0*$Es*$Icol_exmod/$HStory1]; # rotational stiffness of Story 1 column springs
set Ks_beam_2 [expr $n*6.0*$Es*$Ibeam_2mod/$WBay]; # rotational stiffness of Floor 2 beam springs


# set up geometric transformations of element
set PDeltaTransf 1;
geomTransf PDelta $PDeltaTransf; # PDelta transformation

# define elastic column elements using "element" command

# eleID convention: "1xy" where 1 = col, x = Pier #, y = Story #
# Columns Story 1
element elasticBeamColumn 111 117 1001 $Acol_ex $Es $Icol_exmod $PDeltaTransf; # Pier 1
element elasticBeamColumn 112 1001 1002 $Acol_ex $Es $Icol_exmod $PDeltaTransf; # Pier 1
element elasticBeamColumn 113 1002 126 $Acol_ex $Es $Icol_exmod $PDeltaTransf; # Pier 1

element elasticBeamColumn 121 217 1003 $Acol_ex $Es $Icol_exmod $PDeltaTransf; # Pier 2
element elasticBeamColumn 122 1003 1004 $Acol_ex $Es $Icol_exmod $PDeltaTransf; # Pier 2
element elasticBeamColumn 123 1004 226 $Acol_ex $Es $Icol_exmod $PDeltaTransf; # Pier 2


# define elastic beam elements
# eleID convention: "2xy" where 2 = beam, x = Bay #, y = Floor #
# Beams Story 2
element elasticBeamColumn 212 122 2001 $Abeam_2 $Es $Ibeam_2mod $PDeltaTransf;
element elasticBeamColumn 213 2001 2002 $Abeam_2 $Es $Ibeam_2mod $PDeltaTransf;
element elasticBeamColumn 214 2002 2003 $Abeam_2 $Es $Ibeam_2mod $PDeltaTransf;
element elasticBeamColumn 215 2003 2004 $Abeam_2 $Es $Ibeam_2mod $PDeltaTransf;
element elasticBeamColumn 216 2004 223 $Abeam_2 $Es $Ibeam_2mod $PDeltaTransf;


# display the model with the node numbers
DisplayModel2D NodeNumbers

###################################################################################################
# Define Rotational Springs for Plastic Hinges
###################################################################################################
# define rotational spring properties and create spring elements using "rotSpring2DModIKModel" procedure
# input values for column springs
set McMy 1.17; # ratio of capping moment to yield moment, Mc / My
set LS 1000.0; # basic strength deterioration (a very large # = no cyclic deterioration)
set LK 1000.0; # unloading stiffness deterioration (a very large # = no cyclic deterioration)
set LA 1000.0; # accelerated reloading stiffness deterioration (a very large # = no cyclic deterioration)
set LD 1000.0; # post-capping strength deterioration (a very large # = no deterioration)
set cS 1.0; # exponent for basic strength deterioration (c = 1.0 for no deterioration)
set cK 1.0; # exponent for unloading stiffness deterioration (c = 1.0 for no deterioration)
set cA 1.0; # exponent for accelerated reloading stiffness deterioration (c = 1.0 for no deterioration)
set cD 1.0; # exponent for post-capping strength deterioration (c = 1.0 for no deterioration)
set th_pP 0.0015; # plastic rot capacity for pos loading
set th_pN 0.0015; # plastic rot capacity for neg loading
set th_pcP 10; # post-capping rot capacity for pos loading
set th_pcN 10; # post-capping rot capacity for neg loading
set ResP 0.4; # residual strength ratio for pos loading
set ResN 0.4; # residual strength ratio for neg loading
set th_uP 10; # ultimate rot capacity for pos loading
set th_uN 10; # ultimate rot capacity for neg loading
set DP 1.0; # rate of cyclic deterioration for pos loading
set DN 1.0; # rate of cyclic deterioration for neg loading
set a_mem [expr ($n+1.0)*($Mycol_ex*($McMy-1.0)) / ($Ks_col_ex*$th_pP)]; # strain hardening ratio of spring
set b [expr ($a_mem)/(1.0+$n*(1.0-$a_mem))]; # modified strain hardening ratio of spring (Ibarra & Krawinkler 2005, note: Eqn B.5 is incorrect)

# define column springs
# col springs @ bottom of Story 1 (at base)
rotSpring2DModIKModel 3111 11 117 $Ks_col_ex $b $b $Mycol_ex [expr -$Mycol_ex] $LS $LK $LA $LD $cS $cK $cA $cD $th_pP $th_pN $th_pcP $th_pcN $ResP $ResN $th_uP $th_uN $DP $DN;
rotSpring2DModIKModel 3211 21 217 $Ks_col_ex $b $b $Mycol_ex [expr -$Mycol_ex] $LS $LK $LA $LD $cS $cK $cA $cD $th_pP $th_pN $th_pcP $th_pcN $ResP $ResN $th_uP $th_uN $DP $DN;

#col springs @ top of Story 1 (below Floor 2)
rotSpring2DModIKModel 3112 12 126 $Ks_col_ex $b $b $Mycol_ex [expr -$Mycol_ex] $LS $LK $LA $LD $cS $cK $cA $cD $th_pP $th_pN $th_pcP $th_pcN $ResP $ResN $th_uP $th_uN $DP $DN;
rotSpring2DModIKModel 3212 22 226 $Ks_col_ex $b $b $Mycol_ex [expr -$Mycol_ex] $LS $LK $LA $LD $cS $cK $cA $cD $th_pP $th_pN $th_pcP $th_pcN $ResP $ResN $th_uP $th_uN $DP $DN;



# create region for frame column springs
# command: region $regionID -ele $ele_1_ID $ele_2_ID...
region 1 -ele 3111 3211 3112 3212 ;



# define beam springs
# Spring ID: "4xya", where 4 = beam spring, x = Bay #, y = Floor #, a = location in bay
# "a" convention: 1 = left end, 2 = right end

#beam springs at Floor 2
# redefine the rotations since they are not the same
set McMy 1.17; # ratio of capping moment to yield moment, Mc / My
set LS 1.43; # basic strength deterioration (a very large # = no cyclic deterioration)
set LK 1.43; # unloading stiffness deterioration (a very large # = no cyclic deterioration)
set LA 1.43; # accelerated reloading stiffness deterioration (a very large # = no cyclic deterioration)
set LD 1.43; # post-capping strength deterioration (a very large # = no deterioration)
set cS 1.0; # exponent for basic strength deterioration (c = 1.0 for no deterioration)
set cK 1.0; # exponent for unloading stiffness deterioration (c = 1.0 for no deterioration)
set cA 1.0; # exponent for accelerated reloading stiffness deterioration (c = 1.0 for no deterioration)
set cD 1.0; # exponent for post-capping strength deterioration (c = 1.0 for no deterioration)
set th_pP 0.011;
set th_pN 0.011;
set th_pcP 0.155;
set th_pcN 0.155;
set th_uP 0.2; # ultimate rot capacity for pos loading
set th_uN 0.2; # ultimate rot capacity for neg loading
set a_mem [expr ($n+1.0)*($Mybeam_2*($McMy-1.0)) / ($Ks_beam_2*$th_pP)]; # strain hardening ratio of spring
set b [expr ($a_mem)/(1.0+$n*(1.0-$a_mem))]; # modified strain hardening ratio of spring (Ibarra & Krawinkler 2005, note: there is mistake in Eqn B.5)

rotSpring2DModIKModel 4121 12 122 $Ks_beam_2 $b $b $Mybeam_2 [expr -$Mybeam_2] $LS $LK $LA $LD $cS $cK $cA $cD $th_pP $th_pN $th_pcP $th_pcN $ResP $ResN $th_uP $th_uN $DP $DN;
rotSpring2DModIKModel 4122 22 223 $Ks_beam_2 $b $b $Mybeam_2 [expr -$Mybeam_2] $LS $LK $LA $LD $cS $cK $cA $cD $th_pP $th_pN $th_pcP $th_pcN $ResP $ResN $th_uP $th_uN $DP $DN;


# create region for beam springs Floor2
region 2 -ele 4121 4122 ;



############################################################################
# Eigenvalue Analysis
############################################################################
set pi [expr 2.0*asin(1.0)]; # Definition of pi
set nEigenI 1; # mode i = 1
set nEigenJ 2; # mode j = 2
# set nEigenK 3; # mode K = 3
set lambdaN [eigen [expr $nEigenJ]]; # eigenvalue analysis for nEigenJ modes
set lambdaI [lindex $lambdaN [expr 0]]; # eigenvalue mode i = 1
set lambdaJ [lindex $lambdaN [expr $nEigenJ-1]]; # eigenvalue mode j = 2
# set lambdaK [lindex $lambdaN [expr $nEigenK-1]]; # eigenvalue mode j = 3
set w1 [expr pow($lambdaI,0.5)]; # w1 (1st mode circular frequency)
set w2 [expr pow($lambdaJ,0.5)]; # w2 (2nd mode circular frequency)
# set w3 [expr pow($lambdaK,0.5)]; # w3 (3rd mode circular frequency)
set T1 [expr 2.0*$pi/$w1]; # 1st mode period of the structure
set T2 [expr 2.0*$pi/$w2]; # 2nd mode period of the structure
# set T3 [expr 2.0*$pi/$w3]; # 3rd mode period of the structure
puts "T1 = $T1 s"; # display the first mode period in the command window
puts "T2 = $T2 s"; # display the second mode period in the command window
# puts "T3 = $T3 s"; # display the third mode period in the command window

############################################################################

puts "Model Built"

#############################################################################
#Gravity Loads & Gravity Analysis
############################################################################
# define gravity
pattern Plain 1 Linear {
load 12 0 -10000 0
load 22 0 -10000 0
}


# Gravity-analysis: load-controlled static analysis
# define parameters of gravity
constraints Plain
numberer Plain
system BandGeneral
test NormDispIncr 10e-3 20 2
algorithm ModifiedNewton
integrator LoadControl 0.1
analysis Static
analyze 10

# maintain constant gravity loads and reset time to zero
loadConst -time 0.0
puts "Analysis of gravity is completed..."

############################################################################
# Recorders
############################################################################
# record drift histories
# drift recorder command: recorder Drift -file $filename -time -iNode $NodeI_ID -jNode $NodeJ_ID -dof $dof -perpDirn $Record.drift.perpendicular.to.this.direction
recorder Drift -file $dataDir/Drift-Story1.out -time -iNode 11 -jNode 12 -dof 1 -perpDirn 2;


# record floor displacements
recorder Node -file $dataDir/Disp.out -time -node 12 -dof 1 disp;

# record base shear reactions
recorder Node -file $dataDir/Vbase.out -node 117 217 -dof 1 reaction;

# record story 1 column forces in global coordinates
recorder Element -file $dataDir/Fcol111.out -time -ele 111 force;
recorder Element -file $dataDir/Fcol121.out -time -ele 121 force;


# record response history of all frame column springs (one file for moment, one for rotation)
recorder Element -file $dataDir/Col-Mom-Hist.out -time -region 1 force;
recorder Element -file $dataDir/Col-Rot-Hist.out -time -region 1 deformation;

# record response history of all frame beam springs (one file for moment, one for rotation)
recorder Element -file $dataDir/Beam-Mom-Hist.out -time -region 2 force;
recorder Element -file $dataDir/Beam-Rot-Hist.out -time -region 2 deformation;


#######################################################################################
# #
# Analysis Section #
# #
#######################################################################################

############################################################################
# Pushover Analysis #
############################################################################
if {$analysisType == "pushover"} {
puts "Running Pushover..."
# assign lateral loads and create load pattern: use ASCE 7-10 distribution 1.4.3 Fx=0.01Wx
set lat2 200; # force on each frame node in Floor 2


pattern Plain 200 Linear {
load 12 $lat2 0.0 0.0;
load 22 $lat2 0.0 0.0;

}

#display deformed shape:
set ViewScale 5;
DisplayModel2D DeformedShape $ViewScale ; # display deformed shape, the scaling factor needs to be adjusted for each model

#displacement parameters
set IDctrlNode 12; # node where disp is read for disp control
set IDctrlDOF 1; # degree of freedom read for disp control (1 = x displacement)
set Dmax [expr 0.1*$HBuilding]; # maximum displacement of pushover: 10% roof drift
set Dincr [expr 0.5]; # displacement increment

#analysis commands
constraints Plain; # how it handles boundary conditions
numberer RCM; # renumber dof's to minimize band-width (optimization)
system BandGeneral; # how to store and solve the system of equations in the analysis (large model: try UmfPack)
test NormUnbalance 1.0e-1 400 1; # type of convergence criteria with tolerance, max iterations
algorithm Newton; # use Newton's solution algorithm: updates tangent stiffness at every iteration
integrator DisplacementControl $IDctrlNode $IDctrlDOF $Dincr; # use displacement-controlled analysis
analysis Static; # define type of analysis: static for pushover
set Nsteps [expr int($Dmax/$Dincr)];# number of pushover analysis steps
set ok [analyze $Nsteps]; # this will return zero if no convergence problems were encountered
puts "Pushover complete"; # display this message in the command window
}

[fmk]

the problem lies with the parameters passed to the Bilin cmmand in the rotSpring2DModIKModel.tcl file. if i change the metaerial to elastic your script works. so obviously somethin is wrong with the parameters you are passing.

[yuqing]

THX！ However, my friends use this code and solve the problem, could u see why it works?
# define gravity
pattern Plain 1 Linear {
load 12 0 -1000 0
load 22 0 -1000 0
}


# Gravity-analysis: load-controlled static analysis
# define parameters of gravity
set Tol 1.0e-4; # convergence tolerance for test
constraints Transformation
numberer Plain
system BandGeneral
set NstepGravity 10; # apply gravity in 6 steps
set DGravity [expr 1./$NstepGravity]; # first load increment;
integrator LoadControl $DGravity; # determine the next time step for an analysis
puts "applying gravity ..."
for {set i 0} {$i < $NstepGravity} {incr i 1} {
test NormDispIncr $Tol 6; # determine if convergence has been achieved at the end of an iteration step
algorithm Newton; # use Newton's solution algorithm: updates tangent stiffness at every iteration
analysis Static; # define type of analysis static or transient
set ok [analyze 1];
if {$ok != 0} {
puts "trying Newton with Initial Tangent ..."
test NormDispIncr $Tol 10000 0; # determine if convergence has been achieved at the end of an iteration step
algorithm Newton -initial; # use Newton's solution algorithm: updates tangent stiffness at every iteration
analysis Static; # define type of analysis static or transient
set ok [analyze 1];
}
if {$ok != 0} {
# remove recorders;
return -1
}
};
loadConst -time 0.0