I have build a N-Story,N-Bay MRF based on the example in openseeswiki ,Pushover and Dynamic Analyses of 2-Story Moment Frame with Panel Zones and RBS ,but,it didn't work very well.Although I have examine the model many times,I
can't find where my mistakes are.It really makes me very depressed
.I would be very grateful if you can get everyone’s help.Here I will attach the relevant code.#2D model of N-stories 4-Bays SMRFs
#Unitis m,N,seconds
###################################################################################################
# Set Up & Source Definition
###################################################################################################
wipe all;
model basic -ndm 2 -ndf 3;
source LibUnits.tcl;
source DisplayModel2D.tcl;
source DisplayPlane.tcl;
source Steel2d.tcl;
source rotPanelZone2D.tcl;
source elemPanelZone2D.tcl;
####################################################################################################
#Difine analysis type
####################################################################################################
set analysistype "dynamics";
set Datadir Nsout;
file mkdir $Datadir;
####################################################################################################
# Define Building Geometry, Nodes, and Constraints
####################################################################################################
set Nstories 3; # number of stories
set NBays 1; # number of bays
set WBay [expr 9.14*$m]; # width 0f bay
set HNstoryTyp [expr 3.96*$m];#Nstory height of other stories
set Hbuilding [expr $Nstories*$HNstoryTyp]; #height of building;
# display data of the buidling
puts "Nstories=$Nstories "
puts "NBays=$NBays";
puts "WBay=$WBay";
puts "HNstoryTyp = $HNstoryTyp m";
puts "Hbuilding = $Hbuilding m"
# define the cordinate at the base
# the base
for {set Pier 1} {$Pier <= [expr $NBays + 2]} {incr Pier 1 } {
set NodeID [expr 100*$Pier + 1];
set X [expr ($Pier - 1)*$WBay];
node $NodeID $X 0;
if {$Pier != [expr $NBays + 2]} {
node [expr $NodeID*10 + 8] $X 0
}
}
# difine the singular restraints;
for {set Pier 1} {$Pier <= [expr $NBays + 2]} {incr Pier 1} {
set NodeID [expr 100*$Pier + 1];
if {$Pier != [expr $NBays + 2]} {
fix $NodeID 1 1 1;
}
if {$Pier == [expr $NBays + 2]} {
fix $NodeID 1 1 0;
}
}
# define the weight and mass
set WeightTyp [expr 8965*$kN]
set WeightRoof [expr 7372*$kN];
# define the distributitive mass
set MTyp [expr $WeightTyp/$g/5];
set MRoof [expr $WeightRoof/$g/5];
set Negligible 1e-9;
# define some parameters
set N [expr $Nstories+1];
# choose the W14X2333 for the column
# choose the W24X62 for the girder
# define the steel
set Fy [expr 345*$Mpa];
set Es [expr 2e5*$Mpa];
set n 10;
# define the section
# for W14X233 column
set d14233 [expr 407*$mm]; # depth
set bf14233 [expr 404*$mm]; # flange width
set tf14233 [expr 44*$mm]; # flange thickness
set tw14233 [expr 27*$mm]; # web thickness
set A14233 [expr 44193.46*pow($mm,2)]; # area of the section
set Ix14233 [expr 1.2499e+05*pow($cm,4)]; # moment of inertia of the section
set ry14233 [expr 86.85*$mm]; # Week axis
set Sx14233 [expr 7.1396e+03*pow($cm,3)]; # Section Modulus
set I14233mod [expr $Ix14233*($n+1.0)/$n]; # modified moment of inertia
# for W24X62 beam
set d24062 [expr 603*$mm]; # depth
set bf24062 [expr 178.82*$mm]; # flange width
set tf24062 [expr 14.99*$mm]; # flange thickness
set tw24062 [expr 10.92*$mm]; # web thickness
set A24062 [expr 117.86*pow($cm,2)]; # area of the section
set Ix24062 [expr 64813.35*pow($cm,4)]; # moment of inertia of the section
set ry24062 [expr 28.58*$mm]; # Week axis
set Sx24062 [expr 2.4726e+03*pow($cm,3)]; # Section Modulus
set I24062mod [expr $Ix24062*($n+1.0)/$n]; # modified moment of inertia
# difine the parameters for panel zone
set phlat [expr $d14233/2];
set phvert [expr $d24062/2];
# define the other nodes including panel zone
for {set Nstory 2} {$Nstory <= $N} {incr Nstory 1} {
set Y [expr ($Nstory - 1)*$HNstoryTyp];
for {set Pier 1} {$Pier <= [expr $NBays +1]} {incr Pier 1} {
set X [expr ($Pier - 1)*$WBay];
set NodeID [expr {$Pier*100 + $Nstory}];
if {$Nstory != $N && $Pier != 1 && $Pier != 5} {
node [expr 10*$NodeID + 1] [expr $X + $phlat] $Y;
node [expr 10*$NodeID + 3] [expr $X - $phlat] $Y;
node [expr 10*$NodeID + 5] $X [expr $Y - $phvert];
node [expr 10*$NodeID + 6] $X [expr $Y - $phvert];
node [expr 10*$NodeID + 7] $X [expr $Y + $phvert];
node [expr 10*$NodeID + 8] $X [expr $Y + $phvert];
}
if {$Nstory != $N && $Pier == 1} {
node [expr 10*$NodeID + 1] [expr $X + $phlat] $Y;
node [expr 10*$NodeID + 5] $X [expr $Y - $phvert];
node [expr 10*$NodeID + 6] $X [expr $Y - $phvert];
node [expr 10*$NodeID + 7] $X [expr $Y + $phvert];
node [expr 10*$NodeID + 8] $X [expr $Y + $phvert];
}
if {$Nstory != $N && $Pier == 5} {
node [expr 10*$NodeID + 3] [expr $X - $phlat] $Y;
node [expr 10*$NodeID + 5] $X [expr $Y - $phvert];
node [expr 10*$NodeID + 6] $X [expr $Y - $phvert];
node [expr 10*$NodeID + 7] $X [expr $Y + $phvert];
node [expr 10*$NodeID + 8] $X [expr $Y + $phvert];
}
if {$Nstory == $N && $Pier != 1 && $Pier != 5 } {
node [expr 10*$NodeID + 1] [expr $X + $phlat] $Y;
node [expr 10*$NodeID + 3] [expr $X - $phlat] $Y;
node [expr 10*$NodeID + 5] $X [expr $Y - $phvert];
node [expr 10*$NodeID + 6] $X [expr $Y - $phvert];
}
if {$Nstory == $N && $Pier == 1} {
node [expr 10*$NodeID + 1] [expr $X + $phlat] $Y;
node [expr 10*$NodeID + 5] $X [expr $Y - $phvert];
node [expr 10*$NodeID + 6] $X [expr $Y - $phvert];
}
if {$Nstory == $N && $Pier == 5} {
node [expr 10*$NodeID + 3] [expr $X - $phlat] $Y;
node [expr 10*$NodeID + 5] $X [expr $Y - $phvert];
node [expr 10*$NodeID + 6] $X [expr $Y - $phvert];
}
# define the nodes for panel zone
if {$Nstory != $N } {
node [expr $NodeID*100+1] [expr {$X - $phlat}] [expr {$Y + $phvert}];
node [expr $NodeID*100+2] [expr {$X - $phlat}] [expr {$Y + $phvert}];
node [expr $NodeID*100+3] [expr {$X + $phlat}] [expr {$Y + $phvert}];
node [expr $NodeID*100+4] [expr {$X + $phlat}] [expr {$Y + $phvert}];
node [expr $NodeID*100+5] [expr {$X + $phlat}] $Y -mass $MTyp $Negligible $Negligible;
node [expr $NodeID*100+6] [expr {$X + $phlat}] [expr {$Y - $phvert}];
node [expr $NodeID*100+7] [expr {$X + $phlat}] [expr {$Y - $phvert}];
node [expr $NodeID*100+8] [expr {$X - $phlat}] [expr {$Y - $phvert}];
node [expr $NodeID*100+9] [expr {$X - $phlat}] [expr {$Y - $phvert}];
node [expr $NodeID*100+10] [expr {$X - $phlat}] $Y;
}
if {$Nstory == $N} {
node [expr $NodeID*100+1] [expr {$X - $phlat}] [expr {$Y + $phvert}];
node [expr $NodeID*100+2] [expr {$X - $phlat}] [expr {$Y + $phvert}];
node [expr $NodeID*100+3] [expr {$X + $phlat}] [expr {$Y + $phvert}];
node [expr $NodeID*100+4] [expr {$X + $phlat}] [expr {$Y + $phvert}];
node [expr $NodeID*100+5] [expr {$X + $phlat}] $Y -mass $MRoof $Negligible $Negligible;
node [expr $NodeID*100+6] [expr {$X + $phlat}] [expr {$Y - $phvert}];
node [expr $NodeID*100+7] [expr {$X + $phlat}] [expr {$Y - $phvert}];
node [expr $NodeID*100+8] [expr {$X - $phlat}] [expr {$Y - $phvert}];
node [expr $NodeID*100+9] [expr {$X - $phlat}] [expr {$Y - $phvert}];
node [expr $NodeID*100+10] [expr {$X - $phlat}] $Y;
node [expr $NodeID*10+7] $X [expr {$Y + $phvert}];
}
}
}
# define the nodes for PDelta column;
for {set Nstory 2} {$Nstory <= $N} {incr Nstory 1} {
set NodeID [expr ($NBays + 2)*100+$Nstory];
set X [expr ($NBays + 1)*$WBay];
set Y [expr ($Nstory - 1)*$HNstoryTyp];
node $NodeID $X $Y;
if {$Nstory <$N} {
node [expr 10*$NodeID + 5] $X $Y;
node [expr 10*$NodeID + 8] $X $Y;
}
if {$Nstory == $N} {
node [expr 10*$NodeID + 5] $X $Y;
}
}
# define the mutiple restraints
for {set Nstory 2} {$Nstory <= $N} {incr Nstory 1} {
for {set Pier 2} {$Pier <= [expr $NBays + 2]} {incr Pier 1} {
set ConID [expr (100*1 + $Nstory)*100 + 5];
set dof1 1;
if {$Pier != [expr $NBays + 2]} {
equalDOF $ConID [expr (100*$Pier + $Nstory)*100 + 5] $dof1;
}
if {$Pier == [expr $NBays + 2 ]} {
equalDOF $ConID [expr 100*$Pier + $Nstory] $dof1;
}
}
}
####################################################################################################
# Define elements for main structure
####################################################################################################
# define the geomTransformation for elements
set ColTransf 1;
set BeamTransf 2;
set PDeltaTransf 3;
geomTransf PDelta $ColTransf;
geomTransf Linear $BeamTransf;
geomTransf PDelta $PDeltaTransf;
# elements for column
for {set Nstory 1} {$Nstory <= [expr $N - 1]} {incr Nstory 1} {
for {set Pier 1} {$Pier <= [expr $NBays + 1]} {incr Pier 1} {
set eleID [expr 1*1000 + 100*$Pier + $Nstory]
set NU [expr ($Pier*100 + $Nstory)*10 + 8];
set ND [expr ($Pier*100 + $Nstory + 1)*10 + 5]
element elasticBeamColumn $eleID $NU $ND $A14233 $Es $I14233mod $ColTransf;
}
}
# Vis
# elements for beam
for {set Nstory 2} {$Nstory <= $N} {incr Nstory 1} {
for {set Pier 1} {$Pier <= $NBays} {incr Pier 1} {
set eleID [expr 2*1000+$Pier*100+($Nstory)];
set NL [expr (100*$Pier + $Nstory)*10 + 1];
set NR [expr (100*($Pier + 1) + $Nstory)*10 + 3];
element elasticBeamColumn $eleID $NL $NR $A24062 $Es $I24062mod $BeamTransf;
}
}
# define the pdelta column and rigid connections
set TrussMatTag 600;
set Arigid [expr 5*pow($m,2)];
set Irigid [expr 100*pow($m,4)];
uniaxialMaterial Elastic $TrussMatTag $Es;
for {set Nstory 1} {$Nstory <= [expr $N - 1]} {incr Nstory 1} {
set eleID [expr 7*1000+($NBays + 2)*100+$Nstory];
if {$Nstory != 1} {
set NU [expr (($NBays + 2)*100 + $Nstory + 1)*10 + 5];
set ND [expr (($NBays + 2)*100 + $Nstory )*10 + 8];
element elasticBeamColumn $eleID $NU $ND $Arigid $Es $Irigid $PDeltaTransf;
}
if {$Nstory == 1} {
set ND [expr ($NBays + 2)*100 + $Nstory];
set NU [expr (($NBays + 2)*100 + $Nstory + 1)*10 + 5];
element elasticBeamColumn $eleID $NU $ND $Arigid $Es $Irigid $PDeltaTransf;
}
}
for {set Nstory 2} {$Nstory <= $N} {incr Nstory 1} {
set eleID [expr 6*1000+($NBays + 1)*100+$Nstory];
set NL [expr (($NBays + 1)*100 + $Nstory)*100 + 5];
set NR [expr ($NBays + 2)*100 +$Nstory];
element truss $eleID $NL $NR $Arigid $TrussMatTag;
}
# for panel zone
set Apz [expr 5*pow($m,2)];
set Ipz [expr 100*pow($m,4)];
for {set Nstory 2} {$Nstory <= $N} {incr Nstory 1} {
for {set Pier 1} {$Pier <= [expr $NBays + 1]} {incr Pier 1} {
set eleID [expr (5*100000+$Pier*100+$Nstory)*10+1];
set DisN [expr ($Pier*100 + $Nstory)*100 + 1];
elemPanelZone2D $eleID $DisN $Es $Apz $Ipz $PDeltaTransf;
}
}
###################################################################################################
# Define Rotational Springs for Plastic Hinges, Panel Zones, and Leaning Columns
###################################################################################################
# for column
for {set Nstory 1} {$Nstory <= [expr $N - 1]} {incr Nstory 1} {
for {set Pier 1} {$Pier <= [expr $NBays + 1]} {incr Pier 1} {
#define the rotation degratation model(IMK)
set SprIDB [expr (3*1000 + 100*$Pier + $Nstory)*10 + 1];
set SprIDT [expr (3*1000 + 100*$Pier + $Nstory)*10 + 2];
SteelWSectionMR $SprIDB $Es $Fy $I14233mod $Sx14233 $HNstoryTyp $HNstoryTyp $d14233 $tw14233 $bf14233 $tf14233 $HNstoryTyp $ry14233 "other-than-RBS" 0 "-non";
SteelWSectionMR $SprIDT $Es $Fy $I14233mod $Sx14233 $HNstoryTyp $HNstoryTyp $d14233 $tw14233 $bf14233 $tf14233 $HNstoryTyp $ry14233 "other-than-RBS" 0 "-non";
if {$Nstory == 1} {
element zeroLength $SprIDB [expr (100*$Pier + $Nstory)*10 + 8] [expr 100*$Pier +$Nstory] -mat $SprIDB -dir 6
element zeroLength $SprIDT [expr (100*$Pier + $Nstory + 1)*10 + 5] [expr (100*$Pier + $Nstory + 1)*10 + 6] -mat $SprIDT -dir 6
equalDOF [expr 100*$Pier +$Nstory] [expr (100*$Pier + $Nstory)*10 + 8] 1 2;
equalDOF [expr (100*$Pier + $Nstory + 1)*10 + 6] [expr (100*$Pier + $Nstory + 1)*10 + 5] 1 2;
}
if {$Nstory != 1} {
element zeroLength $SprIDB [expr (100*$Pier + $Nstory)*10 + 8] [expr (100*$Pier + $Nstory)*10 + 7] -mat $SprIDB -dir 6;
element zeroLength $SprIDT [expr (100*$Pier + $Nstory + 1)*10 + 5] [expr (100*$Pier + $Nstory + 1)*10 + 6] -mat $SprIDT -dir 6;
equalDOF [expr (100*$Pier + $Nstory)*10 + 7] [expr (100*$Pier + $Nstory)*10 + 8] 1 2;
equalDOF [expr (100*$Pier + $Nstory + 1)*10 + 6] [expr (100*$Pier + $Nstory + 1)*10 + 5] 1 2;
}
}
}
# for beam
for {set Nstory 2} {$Nstory <= $N} {incr Nstory 1} {
for {set Pier 1} {$Pier <= $NBays} {incr Pier 1} {
#define the rotation degratation model(IMK)
set SL [expr (4*1000 + 100*$Pier + $Nstory)*10 + 1];
set SR [expr (4*1000 + 100*$Pier + $Nstory)*10 + 2];
SteelWSectionMR $SL $Es $Fy $I24062mod $Sx24062 $WBay $WBay $d24062 $tw24062 $bf24062 $tf24062 $WBay $ry24062 "other-than-RBS" 0 "-non";
SteelWSectionMR $SR $Es $Fy $I24062mod $Sx24062 $WBay $WBay $d24062 $tw24062 $bf24062 $tf24062 $WBay $ry24062 "other-than-RBS" 0 "-non";
element zeroLength $SL [expr (100*$Pier + $Nstory)*100 +5] [expr (100*$Pier + $Nstory)*10 + 1] -mat $SL -dir 6;
element zeroLength $SR [expr (100*($Pier + 1) + $Nstory)*100 +10] [expr (100*($Pier+ 1) + $Nstory)*10 + 3] -mat $SL -dir 6;
equalDOF [expr (100*$Pier + $Nstory)*100 +5] [expr (100*$Pier + $Nstory)*10 + 1] 1 2;
equalDOF [expr (100*($Pier + 1) + $Nstory)*100 +10] [expr (100*($Pier+ 1) + $Nstory)*10 + 3] 1 2;
}
}
# for PDelta
source rotLeaningCol.tcl
for {set Nstory 1} {$Nstory <= [expr $N - 1]} {incr Nstory 1} {
set SD [expr (5*1000 + ($NBays + 2)*100 + $Nstory)*10 + 1];
set SU [expr (5*1000 + ($NBays + 2)*100 + $Nstory)*10 + 2];
set ND [expr 100*($NBays + 2) + $Nstory];
set NU [expr 100*($NBays + 2) + $Nstory + 1];
if {$Nstory == 1} {
rotLeaningCol $SU $NU [expr $NU*10 + 5];
}
if {$Nstory != 1} {
rotLeaningCol $SU $NU [expr $NU*10 + 5];
rotLeaningCol $SD $ND [expr $ND*10 + 8];
}
}
# define the rotation for the panel zone
set Ry 1.2; # expected yield strength multiplier
set as_PZ 0.03; # strain hardening of panel zones
for {set Nstroy 2} {$Nstory <= $N} {incr Nstory 1} {
for {set Pier 1} {$Pier <= [expr $NBays + 1]} {incr Pier 1} {
set eleID [expr (4*1000 + 100*$Pier + $Nstory)*100];
rotPanelZone2D $eleID [expr (100*$Pier + $Nstory)*100 + 3] [expr (100*$Pier + $Nstory)*100 + 4] $Es $Fy $d14233 $bf14233 $tf14233 $tw14233 $d24062 $Ry $as_PZ;
}
}
DisplayModel2D NodeNumbers;
set ViewScale 5;
DisplayModel2D DeformedShape $ViewScale;
set Lambda [eigen -fullGenLapack 1];
set omega [expr pow($Lambda,0.5)];
set T1 [expr 2*$pi/$omega];
puts "T1 = $T1 s "
############################################################################
# Gravity Loads & Gravity Analysis
############################################################################
set P_PDTyp [expr -(5767*$kN+0.25*31988*$kN)/2];
set P_PDRoof [expr -(5767*$kN+0.25*1605*$kN)/2];
set P_FTyp [expr (-$WeightTyp-$P_PDTyp)/5];
set P_FRoof [expr (-$WeightRoof-$P_PDRoof)/5];
pattern Plain 101 Constant {
load 20205 0 -$WeightRoof 0 (this is simplified for check)
}
# Gravity-analysis: load-controlled static analysis
set Tol 1.0e-6;
constraints Plain;
numberer RCM;
system BandGeneral;
test NormDispIncr $Tol 6;
algorithm Newton;
set NstepGravity 10;
set DGravity [expr 1.0/$NstepGravity];
integrator LoadControl $DGravity;
analysis Static;
analyze $NstepGravity;
# maintain constant gravity loads and reset time to zero
loadConst -time 0.0
puts "Model Built"
And the error attached here:
T1 = 0.0034291009149344702 s
WARNING BandGenLinLapackSolver::solve() -factorization failed, matrix singular U(i,i) = 0, i= 15
WARNING NewtonRaphson::solveCurrentStep() -the LinearSysOfEqn failed in solve()
StaticAnalysis::analyze() - the Algorithm failed at iteration: 0 with domain at load factor 0.1
I konw that means my model is unstable,and maybe wrong somewhere,but ,I can't find where.
I have tried use the same method to build a 2X1 MRF,not using a loop,it can work very well,i.e the method is true and some definitions are also correct.But in the N-Story ,N-Bay,model,it can't work very well.This really bothers me.Sincelely,I will be very grateful to everyone who helped me.