Help ! Where is error?

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

Moderators: silvia, selimgunay, Moderators

Post Reply
zhangchuanchao
Posts: 29
Joined: Wed Oct 15, 2008 2:01 am
Location: china

Help ! Where is error?

Post by zhangchuanchao »

hello everyone, I am analysising a 2-run single-layer structure, however the error is always exist . I can not find out , please help me find it Thank you!

As follows

wipe; # clear memory of all past model definitions
model BasicBuilder -ndm 2 -ndf 3; # Define the model builder, ndm=#dimension, ndf=#dofs
set dataDir Data; # set up name of data directory (you can remove this)
file mkdir $dataDir; # create data directory
set GMdir "../GMfiles/"; # ground-motion file directory
source LibUnits.tcl; #
# define GEOMETRY -------------------------------------------------------------
# define structure-geometry paramters
set LCol [expr 12*$ft]; # column height
set LBeam1 [expr 20*$ft]; # beam length
set LBeam2 [expr 22*$ft];

# calculate locations of beam/column intersections:
set X1 0.;
set X2 [expr $X1 + $LBeam1];
set X3 [expr $X2 + $LBeam2];
set Y1 0.;
set Y2 [expr $Y1 + $LCol];

# define nodal coordinates
node 11 $X1 $Y1
node 12 $X2 $Y1
node 13 $X3 $Y1
node 21 $X1 $Y2
node 22 $X2 $Y2
node 23 $X3 $Y2


# Set up parameters that are particular to the model for displacement control
set IDctrlNode 21; # node where displacement is read for displacement control
set IDctrlDOF 1; # degree of freedom of displacement read for displacement control
set NStory 1; # number of stories above ground level
set NBay 2; # number of bays
set LBuilding $Y2; # total building height

# BOUNDARY CONDITIONS
fix 11 1 1 1
fix 12 1 1 1
fix 13 1 1 1


# Structural-Steel W-section properties
# material properties:
set Es [expr 29028*$ksi]; # Steel Young's Modulus set nu 0.3;
set Gs [expr $Es/2./[expr 1+$nu]]; # Torsional stiffness Modulus


# column sections: A1 60x56 A2 60x66
set AgCol1 [expr 520.8*pow($in,2)]; # cross-sectional area
set IzCol1 [expr 21096*pow($in,4)]; # moment of Inertia
set AgCol2 [expr 613.8*pow($in,2)];
set IzCol2 [expr 34536*pow($in,4)];
# beam sections: W30x56.5
set AgBeam [expr 262.7*pow($in,2)]; # cross-sectional area
set IzBeam [expr 10832*pow($in,4)]; # moment of Inertia

set ColSecTag1 1
set ColSecTag2 2
set BeamSecTag 3
section Elastic $ColSecTag1 $Es $AgCol1 $IzCol1
section Elastic $ColSecTag2 $Es $AgCol2 $IzCol2
section Elastic $BeamSecTag $Es $AgBeam $IzBeam

# define ELEMENTS
# set up geometric transformations of element
# separate columns and beams, in case of P-Delta analysis for columns
set IDColTransf1 1; # all columns
set IDColTransf2 2;
set IDBeamTransf 3; # all beams
set ColTransfType1 Linear ; # options, Linear PDelta Corotational
set ColTransfType2 Linear ;
geomTransf $ColTransfType1 $IDColTransf1 ; # only columns can have PDelta effects (gravity effects)
geomTransf $ColTransfType2 $IDColTransf2 ;
geomTransf Linear $IDBeamTransf


# Define Beam-Column Elements
set np 5; # number of Gauss integration points for nonlinear curvature distribution-- np=2 for linear distribution ok
# columns
element nonlinearBeamColumn 111 11 21 $np $ColSecTag1 $IDColTransf1; # level 1-2
element nonlinearBeamColumn 112 12 22 $np $ColSecTag2 $IDColTransf2
element nonlinearBeamColumn 113 13 23 $np $ColSecTag1 $IDColTransf1
# beams
element nonlinearBeamColumn 221 21 22 $np $BeamSecTag $IDBeamTransf; # level 2
element nonlinearBeamColumn 222 22 23 $np $BeamSecTag $IDBeamTransf;


# Define GRAVITY LOADS, weight and masses
# calculate dead load of frame, assume this to be an internal frame (do LL in a similar manner)
# calculate distributed weight along the beam length
set GammaConcrete [expr 400*$pcf]; # Reinforced-Concrete floor slabs
set Tslab [expr 6*$in]; # 6-inch slab
set Lslab1 [expr 2*$LBeam1/2]; # assume slab extends a distance of $LBeam1/2 in/out of plane
set Lslab2 [expr 2*$LBeam2/2];
set Qslab1 [expr $GammaConcrete*$Tslab*$Lslab1];
set Qslab2 [expr $GammaConcrete*$Tslab*$Lslab2];
set QBeam1 [expr 231*$lbf/$ft]; # W-section weight per length
set QBeam2 [expr 231*$lbf/$ft];
set QdlBeam1 [expr $Qslab1 + $QBeam1]; # dead load distributed along beam.
set QdlBeam2 [expr $Qslab2 + $QBeam2];
set QdlCol1 [expr 617*$lbf/$ft]; # W-section weight per length
set QdlCol2 [expr 718*$lbf/$ft];
set WeightCol1 [expr $QdlCol1*$LCol]; # total Column weight
set WeightCol2 [expr $QdlCol2*$LCol];
set WeightBeam1 [expr $QdlBeam1*$LBeam1]; # total Beam weight
set WeightBeam2 [expr $QdlBeam2*$LBeam2];

# assign masses to the nodes that the columns are connected to
# each connection takes the mass of 1/2 of each element framing into it (mass=weight/$g)
mass 21 [expr ($WeightCol1/2 + $WeightBeam1/2)/$g] 0. 0.; # level 2
mass 22 [expr ($WeightCol2/2 + $WeightBeam1/2 +$WeightBeam2/2)/$g] 0. 0.;
mass 23 [expr ($WeightCol1/2 + $WeightBeam2/2)/$g] 0. 0.;

# calculate total Floor Mass
set WeightFloor2 [expr $WeightCol1*2/2+$WeightCol2*1/2+$WeightBeam1+$WeightBeam2]; # level 2 weight
set WeightTotal $WeightFloor2; # total frame weight
set MassFloor2 [expr $WeightFloor2/$g];
set MassTotal $MassFloor2; # total frame mass

# LATERAL-LOAD distribution for static pushover analysis
# calculate distribution of lateral load based on mass/weight distributions along building height
# Fj = WjHj/sum(WiHi) * Weight at each floor j
set sumWiHi [expr $WeightFloor2*$Y2]; # denominator
set Fj2 [expr $WeightFloor2*$Y2/$sumWiHi*$WeightTotal]; # total for floor 2
set Fi2 [expr $Fj2/4]; # per node on floor 2
set iFi "$Fi2 "; # vectorize

# Define RECORDERS ------------------------------------------------------------
recorder Node -file $dataDir/DFree.out -time -node 21 -dof 1 2 3 disp; # displacements of free node
recorder Node -file $dataDir/DFree.out -time -node 22 -dof 1 2 3 disp;
recorder Node -file $dataDir/DFree.out -time -node 23 -dof 1 2 3 disp;
recorder Node -file $dataDir/DBase.out -time -node 11 12 13 -dof 1 2 3 disp; # displacements of support nodes
recorder Node -file $dataDir/DBase.out -time -node 21 22 23 -dof 1 2 3 disp;
recorder Node -file $dataDir/RBase.out -time -node 11 12 13 -dof 1 2 3 reaction; # support reaction
recorder Drift -file $dataDir/DrNode.out -time -iNode 21 -jNode 11 -dof 1 -perpDirn 2; # lateral drift
recorder Element -file $dataDir/Fel1.out -time -ele 111 localForce; # element forces in local coordinates
recorder Element -file $dataDir/Fel1.out -time -ele 121 localForce;
recorder Element -file $dataDir/Fel1.out -time -ele 131 localForce;
recorder Element -file $dataDir/Fel1.out -time -ele 221 localForce;
recorder Element -file $dataDir/Fel1.out -time -ele 222 localForce;
recorder Element -file $dataDir/ForceEle1sec1.out -time -ele 111 section 1 force; # section forces, axial and moment, node i
recorder Element -file $dataDir/DefoEle1sec1.out -time -ele 111 section 1 deformation; # section deformations, axial and curvature, node i
recorder Element -file $dataDir/ForceEle1sec$np.out -time -ele 111 section $np force; # section forces, axial and moment, node j
recorder Element -file $dataDir/DefoEle1sec$np.out -time -ele 111 section $np deformation; # section deformations, axial and curvature, node j


# define GRAVITY -------------------------------------------------------------
# GRAVITY LOADS # define gravity load applied to beams and columns -- eleLoad applies loads in local coordinate axis
pattern Plain 101 Linear {
eleLoad -ele 221 -type -beamUniform -$QdlBeam1; # beams level 2 (in -ydirection)
eleLoad -ele 222 -type -beamUniform -$QdlBeam2;
eleLoad -ele 111 113 -type -beamUniform 0 -$QdlCol1; # columns level 1-2 (in -xdirection)
eleLoad -ele 112 -type -beamUniform 0 -$QdlCol2;

}
# Gravity-analysis parameters -- load-controlled static analysis
set Tol 1.0e-8; # convergence tolerance for test
variable constraintsTypeGravity Plain; # default;
if { [info exists RigidDiaphragm] == 1} {
if {$RigidDiaphragm=="ON"} {
variable constraintsTypeGravity Lagrange; # large model: try Transformation
}; # if rigid diaphragm is on
}; # if rigid diaphragm exists
constraints $constraintsTypeGravity ; # how it handles boundary conditions
numberer RCM; # renumber dof's to minimize band-width (optimization), if you want to
system BandGeneral ; # how to store and solve the system of equations in the analysis (large model: try UmfPack)
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
set NstepGravity 10; # apply gravity in 10 steps
set DGravity [expr 1./$NstepGravity]; # first load increment;
integrator LoadControl $DGravity; # determine the next time step for an analysis
analysis Static; # define type of analysis static or transient
analyze $NstepGravity; # apply gravity

# ------------------------------------------------- maintain constant gravity loads and reset time to zero
loadConst -time 0.0

puts "Model Built"


# Dynamic Earthquake Analysis
# 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 up analysis parameters
source LibAnalysisDynamicParameters.tcl; # constraintsHandler,DOFnumberer,system-ofequations,convergenceTest,solutionAlgorithm,integrator

# ------------ 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 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]"

The following are wrong tips

Model Built
BandArpackSolver::Error with _saupd info = -9999
Could not build an Arnoldi factorization.IPARAM(5) returns the size of the curre
nt Arnoldi
factorization. The user is advised to check thatenough workspace and array stora
ge has been allocated.
Warning FrequencyAlgo::solveCurrentStep() - the EigenSOE failed in solve().
EigenAnalysis::analyze() - algorithm failed
syntax error in expression "pow(, 0.5)": commas can only separate function argum
ents
while executing
"expr pow($lambdaI, 0.5)"
invoked from within
"set omegaI [expr pow($lambdaI, 0.5)]"
silvia
Posts: 3909
Joined: Tue Jan 11, 2005 7:44 am
Location: Degenkolb Engineers
Contact:

Post by silvia »

the eigensolver is having a hard time, meaning your model is not fully constrained.
Silvia Mazzoni, PhD
Structural Consultant
Degenkolb Engineers
235 Montgomery Street, Suite 500
San Francisco, CA. 94104
silvia
Posts: 3909
Joined: Tue Jan 11, 2005 7:44 am
Location: Degenkolb Engineers
Contact:

Post by silvia »

make sure that the following commands are not giving zero input:
mass 21 [expr ($WeightCol1/2 + $WeightBeam1/2)/$g] 0. 0.; # level 2
mass 22 [expr ($WeightCol2/2 + $WeightBeam1/2 +$WeightBeam2/2)/$g] 0. 0.;
mass 23 [expr ($WeightCol1/2 + $WeightBeam2/2)/$g] 0. 0.;


try putting 2. instead of 2
Silvia Mazzoni, PhD
Structural Consultant
Degenkolb Engineers
235 Montgomery Street, Suite 500
San Francisco, CA. 94104
zhangchuanchao
Posts: 29
Joined: Wed Oct 15, 2008 2:01 am
Location: china

help again

Post by zhangchuanchao »

hello,silvia
sorry ,I can not understand you idea---try putting 2. instead of 2
where?and how to change ---the eigensolver is having a hard time, meaning your model is not fully constrained.
and I still can not find the error
can you point out to me in detail?
thank you!
silvia
Posts: 3909
Joined: Tue Jan 11, 2005 7:44 am
Location: Degenkolb Engineers
Contact:

Post by silvia »

likely, your stiffness is singular, your model is not properly constrained.
Silvia Mazzoni, PhD
Structural Consultant
Degenkolb Engineers
235 Montgomery Street, Suite 500
San Francisco, CA. 94104
Post Reply