2-story MDOF analysis

[shiro]

I'm trying a two-story MDOF analysis.

I created the code as follows and calculated it, but I got the following error.

If you understand the cause, I would appreciate it if you could let me know.



WARNING: analysis .. existing StaticAnalysis exists => wipeAnalysis not invoked, problems may arise
WARNING BandGenLinLapackSolver::solve() -factorization failed, matrix singular U(i,i) = 0, i= 2
WARNING BandGenLinLapackSolver::solve() -factorization failed, matrix singular U(i,i) = 0, i= 2
WARNING BandGenLinLapackSolver::solve() -factorization failed, matrix singular U(i,i) = 0, i= 2
WARNING BandGenLinLapackSolver::solve() -factorization failed, matrix singular U(i,i) = 0, i= 2
WARNING BandGenLinLapackSolver::solve() -factorization failed, matrix singular U(i,i) = 0, i= 2
WARNING BandGenLinLapackSolver::solve() -factorization failed, matrix singular U(i,i) = 0, i= 2
ArpackSolver::Error with _saupd info = -9999
Could not build an Arnoldi factorization.IPARAM(5) the size of the current Arnoldi factorization: is 2factorization. The user is advised to check thatenough workspace and array storage has been allocated.
WARNING DirectIntegrationAnalysis::eigen() - EigenSOE failed in solve()
missing function argument at _@_
in expression "pow(_@_,0.5)"



# ------------------------------
# Start of model generation
# ------------------------------
# Create ModelBuilder (with two-dimensions and 3 DOF/node)
model BasicBuilder -ndm 2 -ndf 3

source DisplayModel2D.tcl; # procedure for displaying a 2D perspective of model
source DisplayPlane.tcl; # procedure for displaying a plane in a model

source units.tcl;

set GMdir "./GMfiles/"; # ground-motion file directory

#####################################################################
# Define Analysis Type #
######################################################################

# Define type of analysis: "pushover" = pushover; "dynamic" = dynamic

set analysisType "dynamic";

#set analysisType "cyclic";

#set analysisType "pushover";

#set analysisType "static";

if {$analysisType == "static"} {
set dataDir Concentrated-Static-Output; # name of output folder
file mkdir $dataDir; # create output folder
}
if {$analysisType == "pushover"} {
set dataDir Concentrated-Pushover-Output; # name of output folder
file mkdir $dataDir; # create output folder
}
if {$analysisType == "cyclic"} {
set dataDir Concentrated-Cyclic-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 geometry for model
# -------------------------
set g [expr 32.174*12.0]
#set P 9.0
set P 1.0

set mass [expr $P/$g]

node 1 0.0 0.0
node 2 0.0 0.0 -mass 0.0 $mass 0.0
node 3 0.0 0.0 -mass 0.0 $mass 0.0

# set the boundary conditions
# tag DX DY RZ
fix 1 1 1 1

# Define material models

# Physical and mechanical characteristics of SWP :

## Materials
set F0 [expr 8.0*0.225] ;# kip
set FI [expr 1.20*0.225] ;# kip
set DU [expr 15.0/25.4] ;# in
set S0 [expr 5.0*0.225*25.4] ;# kip/in
set R1 0.058
set R2 -0.050
set R3 1.00
set R4 0.020
set alph 0.60
set bet 1.10

# material ID

set matID 1

# add the material to domain through the use of a procedure

uniaxialMaterial SAWS $matID $F0 $FI $DU $S0 $R1 $R2 $R3 $R4 $alph $bet

uniaxialMaterial Elastic 2 2.03e+11

# Define geometric transformations
# --------------------------------
# geomTransf type tag
geomTransf Linear 1

# Define elements
# ---------------

element zeroLength 1 1 2 -mat 1 -dir 1
element zeroLength 2 2 3 -mat 1 -dir 1


# ------------------------------
# Start of analysis generation
# ------------------------------
# create the system of equation
system BandGeneral
# create the DOF numberer
numberer Plain
# create the constraint handler
constraints Plain
# create the convergence test
test NormDispIncr 1.0e-12 10
# create the integration scheme
integrator LoadControl 0.1
# create the solution algorithm
algorithm Newton
# create the analysis object
analysis Static
# ------------------------------
# End of analysis generation
# ------------------------------

# ------------------------------
# Perform the gravity analysis
# ------------------------------
# perform the gravity load analysis, requires 10 steps to reach the load level
#analyze 10
#puts "\nGravity load analysis completed";

# set the gravity loads to be constant & reset the time in the domain
loadConst -time 0.0
remove recorders

# ------------------------------
# Start of recorder generation
# ------------------------------
# create a Recorder object for the nodal displacements at node 2
recorder Node -file Node_Dsp.out -time -node 2 3 -dof 1 disp

recorder Element -file Elmt_Force.out -time -ele 1 2 force

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


############################################################################
# Pushover Analysis #
############################################################################
if {$analysisType == "pushover"} {
puts "Running Pushover..."
# assign lateral loads and create load pattern: use ASCE 7-10 distribution
# set lat2 16.255; # force on each frame node in Floor 2
# set lat3 31.636; # force on each frame node in Floor 3
# pattern Plain 200 Linear {
# load 12 $lat2 0.0 0.0;
# load 22 $lat2 0.0 0.0;
# load 13 $lat3 0.0 0.0;
# load 23 $lat3 0.0 0.0;
# }

# set lat5 [expr (1.0 * 0.22481)]; # force on each frame node in Floor 5 P=1kN
set lat5 1.0; # force on each frame node in Floor 5 P=1kip
pattern Plain 200 Linear {
load 15 $lat5 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 13; # node where disp is read for disp control

set IDctrlNode 15; # 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.01]; # 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-6 400; # 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
}

######################################################################
# Time History/Dynamic Analysis #
######################################################################

if {$analysisType == "dynamic"} {
puts "Running dynamic analysis..."
# display deformed shape:
# set ViewScale 2; # amplify display of deformed shape
# DisplayModel2D DeformedShape $ViewScale;


# source in procedures
source ReadSMDfile.tcl; # procedure for reading GM file and converting it to proper format


# Uniform Earthquake ground motion (uniform acceleration input at all support nodes)
set GMdirection 1; # ground-motion direction
set GMfile "H-e12140" ; # ground-motion filenames
set GMfact 1.5; # ground-motion scaling factor

# set up ground-motion-analysis parameters
set DtAnalysis [expr 0.01*$sec]; # time-step Dt for lateral analysis
set TmaxAnalysis [expr 10. *$sec]; # maximum duration of ground-motion analysis -- should be 50*$sec
#set TmaxAnalysis [expr 20. *$sec]; # maximum duration of ground-motion analysis -- should be 50*$sec

# ----------- set up analysis parameters
source LibAnalysisDynamicParameters.tcl; # constraintsHandler,DOFnumberer,system-ofequations,convergenceTest,solutionAlgorithm,integrator


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


# ------------ define & apply damping
# RAYLEIGH damping parameters, Where to put M/K-prop damping, switches (http://opensees.berkeley.edu/OpenSees/m ... l/1099.htm)
# D=$alphaM*M + $betaKcurr*Kcurrent + $betaKcomm*KlastCommit + $beatKinit*$Kinitial
set xDamp 0.02; # damping ratio
set MpropSwitch 1.0;
set KcurrSwitch 0.0;
set KcommSwitch 1.0;
set KinitSwitch 0.0;
set nEigenI 1; # mode 1
#set nEigenJ 3; # mode 3

set nEigenJ 2; # mode 2

set lambdaN [eigen [expr $nEigenJ]]; # eigenvalue analysis for nEigenJ modes
set lambdaI [lindex $lambdaN [expr $nEigenI-1]]; # eigenvalue mode i
set lambdaJ [lindex $lambdaN [expr $nEigenJ-1]]; # eigenvalue mode j
set omegaI [expr pow($lambdaI,0.5)];
set omegaJ [expr pow($lambdaJ,0.5)];
set alphaM [expr $MpropSwitch*$xDamp*(2*$omegaI*$omegaJ)/($omegaI+$omegaJ)]; # M-prop. damping; D = alphaM*M
set betaKcurr [expr $KcurrSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # current-K; +beatKcurr*KCurrent
set betaKcomm [expr $KcommSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # last-committed K; +betaKcomm*KlastCommitt
set betaKinit [expr $KinitSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # initial-K; +beatKinit*Kini
rayleigh $alphaM $betaKcurr $betaKinit $betaKcomm; # RAYLEIGH damping

# --------------------------------- perform Dynamic Ground-Motion Analysis
# the following commands are unique to the Uniform Earthquake excitation
set IDloadTag 400; # for uniformSupport excitation
# Uniform EXCITATION: acceleration input
set inFile $GMdir$GMfile.at2
set outFile $GMdir$GMfile.g3; # set variable holding new filename (PEER files have .at2/dt2 extension)
ReadSMDFile $inFile $outFile dt; # call procedure to convert the ground-motion file
set GMfatt [expr $g*$GMfact]; # data in input file is in g Unifts -- ACCELERATION TH
set AccelSeries "Series -dt $dt -filePath $outFile -factor $GMfatt"; # time series information
pattern UniformExcitation $IDloadTag $GMdirection -accel $AccelSeries ; # create Unifform excitation

set Nsteps [expr int($TmaxAnalysis/$DtAnalysis)];
set ok [analyze $Nsteps $DtAnalysis]; # actually perform analysis; returns ok=0 if analysis was successful

if {$ok != 0} { ; # analysis was not successful.
# --------------------------------------------------------------------------------------------------
# change some analysis parameters to achieve convergence
# performance is slower inside this loop
# Time-controlled analysis
set ok 0;
set controlTime [getTime];
while {$controlTime < $TmaxAnalysis && $ok == 0} {
set controlTime [getTime]
set ok [analyze 1 $DtAnalysis]
if {$ok != 0} {
puts "Trying Newton with Initial Tangent .."
test NormDispIncr $Tol 1000 0
algorithm Newton -initial
set ok [analyze 1 $DtAnalysis]
test $testTypeDynamic $TolDynamic $maxNumIterDynamic 0
algorithm $algorithmTypeDynamic
}
if {$ok != 0} {
puts "Trying Broyden .."
algorithm Broyden 8
set ok [analyze 1 $DtAnalysis]
algorithm $algorithmTypeDynamic
}
if {$ok != 0} {
puts "Trying NewtonWithLineSearch .."
algorithm NewtonLineSearch .8
set ok [analyze 1 $DtAnalysis]
algorithm $algorithmTypeDynamic
}
}
}; # end if ok !0


puts "Ground Motion Done. End Time: [getTime]"

};

[mhscott]

You have 6 degrees of freedom (two nodes not fixed) but the zero length elements only act in two DOFs. The stiffness matrix is singular (zeros on the diagonal). Add more fix commands to the vertical and rotational DOFs of nodes 2 and 3.

[selimgunay]

You should add

fix 2 0 1 1
fix 3 0 1 1

[shiro]

Dear mhscott and selimgunay

When add more fix commands to the vertical and rotational DOFs of nodes 2 and 3 was set as follows and calculated,
the following error occurred.

Thank you for your advice.

# set the boundary conditions
# tag DX DY RZ
fix 1 1 1 1
fix 2 0 1 1
fix 3 0 1 1

WARNING: analysis .. existing StaticAnalysis exists => wipeAnalysis not invoked, problems may arise
ArpackSolver::Error with _saupd info = -3
NCV must be greater than NEV and less than or equal to N.
WARNING DirectIntegrationAnalysis::eigen() - EigenSOE failed in solve()
missing function argument at _@_
in expression "pow(_@_,0.5)"

[selimgunay]

You can try

set lambdaN [eigen -fullGenLapack [expr $nEigenJ]]

[selimgunay]

If that does not work, you can try adding negligible masses instead of the zero masses

node 2 0.0 0.0 -mass 0.0 $mass 0.0
node 3 0.0 0.0 -mass 0.0 $mass 0.0

[shiro]

Dear selimgunay

I changed lambdaN as follows and succeeded.

When masses were changed, the X direction displacement became "nan" and an error occurred.

Thank you for the appropriate advice.

[shiro]

One last question.

Assume that the ground is node1, the height of the first floor is node2, and the height of the second floor is node3.

The element numbers were set in order with the zeroLength element,
but is it calculated automatically considering the floor in the MDOF stiffness matrix?

element zeroLength 1 1 2 -mat 1 -dir 1
element zeroLength 2 2 3 -mat 1 -dir 1


Also, the displacement and stress were recorded as follows.
Is the response displacement of node2 the absolute displacement,
and is the displacement of node1 minus the relative displacement of the second floor?

The numbers for Force and localForce were the same.

Thank you for your reply.


recorder Node -file Node_Disp.out -time -node 2 3 -dof 1 disp

recorder Element -file Elmt_Force.out -time -ele 1 2 localForce

recorder Element -file Elmt_localForce.out -time -ele 1 2 localForce

[mhscott]

If I understand your first question correctly, the stiffness matrix is based on equation numbers, not the element numbers.

The displacements with UniformExcitation are relative to the ground.

[selimgunay]

recorder Node -file Node_Disp.out -time -node 2 3 -dof 1 disp: Here the displacements of nodes 2 and 3 are recorded with respect to the ground.

If you want to record the story drift ratio, which is [(displacement of node 3) minus (displacement of node 2)] /length, you can do it either a) with postprocessing or b) use a drift recorder: https://opensees.berkeley.edu/wiki/inde ... t_Recorder

[shiro]

Dear mhscott and selimgunay

thank you for your answer.

Does the meaning of following the equation number mean that the following 2DOFs are arranged in
the following stiffness matrix and mass matrix in the shear mass system?

I would appreciate if you could answer.


node 1 0.0 0.0
node 2 0.0 0.0 -mass $mass1 0.0 0.0
node 3 0.0 0.0 -mass $mass2 0.0 0.0

fix 1 1 1 1
fix 2 0 1 1
fix 3 0 1 1

uniaxialMaterial SAWS 1 $F0 $FI $DU $S0 $R1 $R2 $R3 $R4 $alph $bet
uniaxialMaterial SAWS 2 $F0 $FI $DU $S0 $R1 $R2 $R3 $R4 $alph $bet

element zeroLength 1 1 2 -mat 1 -dir 1
element zeroLength 2 2 3 -mat 2 -dir 1


M = | m1 0 | = | mass1(node2) 0 |
| 0 m2 | | 0 mass2(node3) |

K = | k2+k1 -k2 | = | k(element zeroLength 2 -mat 2)+k(element zeroLength 1 -mat 1) -k(element zeroLength 2 -mat 2) |
| -k2 k2 | | -k(element zeroLength 2 -mat 2) k(element zeroLength 2 -mat 2) |

[selimgunay]

Yes that is correct.

[shiro]

For small houses, I thought that it might be more convenient to analyze with MDOF rather than 2D or 3D analysis.

Thank you for all the advice.