Dear All,
I have a 2D cantilever reinforced concrete column modeled in OpenSees. The Column section is molded with fiber section. Also, the first element at the bottom of the column is a zero length element with lines springs to simulate the effect of foundation. I want to perform an eigen value analysis. The column as concentrated mass at the top node. I considered the mass only acting in horizontal degree of freedom. Therefore, the only solver that can solve eigen value problem is "fullGenLapack". This one returns the fist eigen value with corresponding natural frequency equal to about 50 rad/sec. With some Test I figured out that this cannot be the correct natural frequency. I tried another option. I considered the mass in vertical DOF as well like:
mass $NodeNumber $ConsColMass $ConsColMass 0.0;
and tried both "fullGenLapack" solver and the default solver. fullGenLapack returns the same frequency as before, but the default solver returns the frequency about 17 rad/sec which is more close the what I expect (basically I tried to excite the column under a sine excitation with different loading frequencies to see in which loading frequency resonance happens).
Can anyone tell me why these two give me two different answers? And also why the fullGenLapack is returning the wrong frequency?
Eigen Value Analysis
Moderators: silvia, selimgunay, Moderators
-
amkopensees
- Posts: 24
- Joined: Mon Jun 18, 2012 9:29 am
- Location: Lehigh University
Re: Eigen Value Analysis
it should not happen. did you check the eigenvalues to see which mode shapes you are getting.
here is a 2d RC concrete column example that demonstrates correct behaviour (NOTE: to get the default working even though you only have 1 mass just use a very small value for some other dof):
set solverTypes {-genBandArpack -fullGenLapack}
foreach solverType $solverTypes {
foreach massDirn {1 2} {
wipe
model basic -ndm 2 -ndf 3
set width 360
set height 144
node 1 0.0 0.0
node 3 0.0 $height
fix 1 1 1 1
uniaxialMaterial Concrete01 1 -6.0 -0.004 -5.0 -0.014
uniaxialMaterial Concrete01 2 -5.0 -0.002 0.0 -0.006
set fy 60.0; # Yield stress
set E 30000.0; # Young's modulus
# tag fy E0 b
uniaxialMaterial Steel01 3 $fy $E 0.01
#uniaxialMaterial Steel01 1 $fy $E 0.01
#uniaxialMaterial Steel01 2 $fy $E 0.01
# Define cross-section for nonlinear columns
# ------------------------------------------
# set some paramaters
set colWidth 15
set colDepth 24
set cover 1.5
set As 0.60; # area of no. 7 bars
# some variables derived from the parameters
set y1 [expr $colDepth/2.0]
set z1 [expr $colWidth/2.0]
section Fiber 1 {
# Create the concrete core fibers
patch rect 1 10 1 [expr $cover-$y1] [expr $cover-$z1] [expr $y1-$cover] [expr $z1-$cover]
# Create the concrete cover fibers (top, bottom, left, right)
patch rect 2 10 1 [expr -$y1] [expr $z1-$cover] $y1 $z1
patch rect 2 10 1 [expr -$y1] [expr -$z1] $y1 [expr $cover-$z1]
patch rect 2 2 1 [expr -$y1] [expr $cover-$z1] [expr $cover-$y1] [expr $z1-$cover]
patch rect 2 2 1 [expr $y1-$cover] [expr $cover-$z1] $y1 [expr $z1-$cover]
# Create the reinforcing fibers (left, middle, right)
layer straight 3 3 $As [expr $y1-$cover] [expr $z1-$cover] [expr $y1-$cover] [expr $cover-$z1]
layer straight 3 2 $As 0.0 [expr $z1-$cover] 0.0 [expr $cover-$z1]
layer straight 3 3 $As [expr $cover-$y1] [expr $z1-$cover] [expr $cover-$y1] [expr $cover-$z1]
}
# Define column elements
# ----------------------
# Geometry of column elements
# tag
geomTransf Linear 1
set np 5
set eleType forceBeamColumn
element $eleType 1 1 3 $np 1 1
# Define nodal mass in terms of axial load on columns
set P 10.0
set g 386.4
set m [expr $P/$g]; # expr command to evaluate an expression
if {$massDirn == 2} {
mass 3 $m $m 0
} else {
mass 3 $m [expr $m/1e12] 0
}
puts "#dirn mass acting: $massDirn solverType: $solverType eigenvalue: [eigen $solverType 1]"
}
}
when you run it you should see something like:
#dirn mass acting: 1 solverType: -genBandArpack eigenvalue: 3.079542e+03
#dirn mass acting: 2 solverType: -genBandArpack eigenvalue: 3.079542e+03
#dirn mass acting: 1 solverType: -fullGenLapack eigenvalue: 3.079542e+03
#dirn mass acting: 2 solverType: -fullGenLapack eigenvalue: 3.079542e+03
here is a 2d RC concrete column example that demonstrates correct behaviour (NOTE: to get the default working even though you only have 1 mass just use a very small value for some other dof):
set solverTypes {-genBandArpack -fullGenLapack}
foreach solverType $solverTypes {
foreach massDirn {1 2} {
wipe
model basic -ndm 2 -ndf 3
set width 360
set height 144
node 1 0.0 0.0
node 3 0.0 $height
fix 1 1 1 1
uniaxialMaterial Concrete01 1 -6.0 -0.004 -5.0 -0.014
uniaxialMaterial Concrete01 2 -5.0 -0.002 0.0 -0.006
set fy 60.0; # Yield stress
set E 30000.0; # Young's modulus
# tag fy E0 b
uniaxialMaterial Steel01 3 $fy $E 0.01
#uniaxialMaterial Steel01 1 $fy $E 0.01
#uniaxialMaterial Steel01 2 $fy $E 0.01
# Define cross-section for nonlinear columns
# ------------------------------------------
# set some paramaters
set colWidth 15
set colDepth 24
set cover 1.5
set As 0.60; # area of no. 7 bars
# some variables derived from the parameters
set y1 [expr $colDepth/2.0]
set z1 [expr $colWidth/2.0]
section Fiber 1 {
# Create the concrete core fibers
patch rect 1 10 1 [expr $cover-$y1] [expr $cover-$z1] [expr $y1-$cover] [expr $z1-$cover]
# Create the concrete cover fibers (top, bottom, left, right)
patch rect 2 10 1 [expr -$y1] [expr $z1-$cover] $y1 $z1
patch rect 2 10 1 [expr -$y1] [expr -$z1] $y1 [expr $cover-$z1]
patch rect 2 2 1 [expr -$y1] [expr $cover-$z1] [expr $cover-$y1] [expr $z1-$cover]
patch rect 2 2 1 [expr $y1-$cover] [expr $cover-$z1] $y1 [expr $z1-$cover]
# Create the reinforcing fibers (left, middle, right)
layer straight 3 3 $As [expr $y1-$cover] [expr $z1-$cover] [expr $y1-$cover] [expr $cover-$z1]
layer straight 3 2 $As 0.0 [expr $z1-$cover] 0.0 [expr $cover-$z1]
layer straight 3 3 $As [expr $cover-$y1] [expr $z1-$cover] [expr $cover-$y1] [expr $cover-$z1]
}
# Define column elements
# ----------------------
# Geometry of column elements
# tag
geomTransf Linear 1
set np 5
set eleType forceBeamColumn
element $eleType 1 1 3 $np 1 1
# Define nodal mass in terms of axial load on columns
set P 10.0
set g 386.4
set m [expr $P/$g]; # expr command to evaluate an expression
if {$massDirn == 2} {
mass 3 $m $m 0
} else {
mass 3 $m [expr $m/1e12] 0
}
puts "#dirn mass acting: $massDirn solverType: $solverType eigenvalue: [eigen $solverType 1]"
}
}
when you run it you should see something like:
#dirn mass acting: 1 solverType: -genBandArpack eigenvalue: 3.079542e+03
#dirn mass acting: 2 solverType: -genBandArpack eigenvalue: 3.079542e+03
#dirn mass acting: 1 solverType: -fullGenLapack eigenvalue: 3.079542e+03
#dirn mass acting: 2 solverType: -fullGenLapack eigenvalue: 3.079542e+03
-
amkopensees
- Posts: 24
- Joined: Mon Jun 18, 2012 9:29 am
- Location: Lehigh University
Re: Eigen Value Analysis
Dear fmk,
I think you forgot to put the solver type right at the end of the code. I mean you have:
puts "#dirn mass acting: $massDirn solverType: $solverType eigenvalue: [eigen 1]"
So, the solver is always the default. I think it should be as follows if you want to try it with both solvers:
puts "#dirn mass acting: $massDirn solverType: $solverType eigenvalue: [eigen $solverType 1]"
With the model you posted its ok. But for the case of my column I tried and again I got completely different answers from these two solvers. My model is as follows. I appreciate it if you take a look at it and tell me what is wrong.
set solverTypes {-genBandArpack -fullGenLapack}
foreach solverType $solverTypes {
foreach massDirn {1 2} {
wipe
model BasicBuilder -ndm 2 -ndf 3
set spe 12.2e3;
set spm 24.4e3;
set colh 4.6e3;
set coldiam 914.4;
set colrad [expr $coldiam/2.0];
set colcov 50.8;
set founh 1092;
set kpileh 8.8207;
set kpilev 155.7795;
set npiles 8.0;
set founkh [expr $npiles*$kpileh];
set founkrh [expr $kpilev*6.0*0.76e3*0.76e3];
set deckWe 0.04279;
set deckMe [expr $deckWe/9.98e3];
set deckWm 0.05604;
set deckMm [expr $deckWm/9.98e3];
set fc -0.04179;
set epsc0 -2.0e-3;
set fcu 0.0;
set epsu -0.0055;
set fy 0.43751;
set Es 200.0;
set hr 0.018;
set Rein_As 645.0;
set rho_s 1.031e-3;
set K_confined [expr 1.0+($rho_s*$fy)/(-1.0*$fc)];
set fcc [expr $K_confined*$fc];
set epsc0c [expr $K_confined*$epsc0];
set fcuc -8.0e-3;
set epsuc -0.05;
set ConsColMass [expr $deckMe*$spe/6.0+$deckMm*$spm/6.0]; # consentrated nodal mass at the top
set DampingRatio 0.06649;
set AxialLoadNode 4210;
set AxialLoad [expr -1.0*($deckWe*$spe/6.0+$deckWm*$spm/6.0)]
# Unconfined Concrete
set ConcTag 700
# uniaxialMaterial Concrete01 $matTag $fpc $epsc0 $fpcu $epsU
uniaxialMaterial Concrete01 $ConcTag $fc $epsc0 $fcu $epsu
# Confined Concrete
set ConfConcTag 800
# uniaxialMaterial Concrete01 $matTag $fpc $epsc0 $fpcu $epsU
uniaxialMaterial Concrete01 $ConfConcTag $fcc $epsc0c $fcuc $epsuc
# Reinforceing Steel
set ReinTag 900
# uniaxialMaterial Steel01 $matTag $Fy $E0 $b <$a1 $a2 $a3 $a4>
uniaxialMaterial Steel01 $ReinTag $fy $Es $hr
# Defining fiber cross-section
set ColSecTag 100
section Fiber $ColSecTag {
# Confined concrete layers in the core
set ConfConc_numSubdivCirc 12
set ConfConc_numSubdivRad 8
set ConfConc_yCenter 0.0
set ConfConc_zCenter 0.0
set ConfConc_intRad 0.0
set ConfConc_extRad [expr $colrad-$colcov]
set ConfConc_startAng 0.0
set ConfConc_endAng 360.0
# patch circ $matTag $numSubdivCirc $numSubdivRad $yCenter $zCenter $intRad $extRad <$startAng $endAng>
patch circ $ConfConcTag $ConfConc_numSubdivCirc $ConfConc_numSubdivRad $ConfConc_yCenter $ConfConc_zCenter $ConfConc_intRad $ConfConc_extRad $ConfConc_startAng $ConfConc_endAng
# Unconfined concrete layers in cover
set Conc_numSubdivCirc 12
set Conc_numSubdivRad 2
set Conc_yCenter 0.0
set Conc_zCenter 0.0
set Conc_intRad [expr $colrad-$colcov]
set Conc_extRad $colrad
set Conc_startAng 0.0
set Conc_endAng 360.0
# patch circ $matTag $numSubdivCirc $numSubdivRad $yCenter $zCenter $intRad $extRad <$startAng $endAng>
patch circ $ConcTag $Conc_numSubdivCirc $Conc_numSubdivRad $Conc_yCenter $Conc_zCenter $Conc_intRad $Conc_extRad $Conc_startAng $Conc_endAng
# Reinforcement leyers
set St_numFiber 12
set St_areaFiber $Rein_As
set St_yCenter 0.0
set St_zCenter 0.0
set St_radius [expr $colrad-$colcov]
set St_startAng 0.0
set St_endAng 330.0
# layer circ $matTag $numFiber $areaFiber $yCenter $zCenter $radius <$startAng $endAng>
layer circ $ReinTag $St_numFiber $St_areaFiber $St_yCenter $St_zCenter $St_radius $St_startAng $St_endAng
}
set large 1.0e10;
#======================================== Foundation ============================================#
set HFoundTag 110
set RFoundTag 120
uniaxialMaterial Elastic $HFoundTag $founkh
uniaxialMaterial Elastic $RFoundTag $founkrh
set ColTrsTag 2;
geomTransf Linear $ColTrsTag;
node 4301 0.0 0.0;
node 4302 0.0 0.0;
node 4303 0.0 546.0;
node 4304 0.0 852.66;
node 4305 0.0 1159.32;
node 4306 0.0 1465.98;
node 4307 0.0 1772.64;
node 4308 0.0 2079.3;
node 4309 0.0 3612.66;
node 4310 0.0 5146.0;
node 4210 0.0 5756.0;
fix 4301 1 1 1;
fix 4302 0 1 0;
#========================================= Nodal Mass ==========================================#
puts "Mass= $ConsColMass"
set MassNode 4210
if {$massDirn == 2} {
mass $MassNode $ConsColMass $ConsColMass 0.0;
} else {
mass $MassNode $ConsColMass [expr $ConsColMass/1e12] 0.0;
}
# element zeroLength $eleTag $iNode $jNode -mat $matTag1 $matTag2 -dir $dir1 $dir2 ...<-doRayleigh $rFlag>
element zeroLength 4301 4301 4302 -mat $HFoundTag $RFoundTag -dir 1 6;
# element elasticBeamColumn $eleTag $iNode $jNode $A $E $Iz $transfTag <-mass $massDens>
element elasticBeamColumn 4302 4302 4303 $large $large $large $ColTrsTag
# element dispBeamColumn $eleTag $iNode $jNode $numIntgrPts $secTag $transfTag <-mass $massDens> <-integration $intType>
element dispBeamColumn 4303 4303 4304 5 $ColSecTag $ColTrsTag
element dispBeamColumn 4304 4304 4305 5 $ColSecTag $ColTrsTag
element dispBeamColumn 4305 4305 4306 5 $ColSecTag $ColTrsTag
element dispBeamColumn 4306 4306 4307 5 $ColSecTag $ColTrsTag
element dispBeamColumn 4307 4307 4308 5 $ColSecTag $ColTrsTag
element dispBeamColumn 4308 4308 4309 5 $ColSecTag $ColTrsTag
element dispBeamColumn 4309 4309 4310 5 $ColSecTag $ColTrsTag
element elasticBeamColumn 4310 4310 4210 $large $large $large $ColTrsTag
puts "#dirn mass acting: $massDirn solverType: $solverType eigenvalue: [eigen $solverType 1]"
}
}
I think you forgot to put the solver type right at the end of the code. I mean you have:
puts "#dirn mass acting: $massDirn solverType: $solverType eigenvalue: [eigen 1]"
So, the solver is always the default. I think it should be as follows if you want to try it with both solvers:
puts "#dirn mass acting: $massDirn solverType: $solverType eigenvalue: [eigen $solverType 1]"
With the model you posted its ok. But for the case of my column I tried and again I got completely different answers from these two solvers. My model is as follows. I appreciate it if you take a look at it and tell me what is wrong.
set solverTypes {-genBandArpack -fullGenLapack}
foreach solverType $solverTypes {
foreach massDirn {1 2} {
wipe
model BasicBuilder -ndm 2 -ndf 3
set spe 12.2e3;
set spm 24.4e3;
set colh 4.6e3;
set coldiam 914.4;
set colrad [expr $coldiam/2.0];
set colcov 50.8;
set founh 1092;
set kpileh 8.8207;
set kpilev 155.7795;
set npiles 8.0;
set founkh [expr $npiles*$kpileh];
set founkrh [expr $kpilev*6.0*0.76e3*0.76e3];
set deckWe 0.04279;
set deckMe [expr $deckWe/9.98e3];
set deckWm 0.05604;
set deckMm [expr $deckWm/9.98e3];
set fc -0.04179;
set epsc0 -2.0e-3;
set fcu 0.0;
set epsu -0.0055;
set fy 0.43751;
set Es 200.0;
set hr 0.018;
set Rein_As 645.0;
set rho_s 1.031e-3;
set K_confined [expr 1.0+($rho_s*$fy)/(-1.0*$fc)];
set fcc [expr $K_confined*$fc];
set epsc0c [expr $K_confined*$epsc0];
set fcuc -8.0e-3;
set epsuc -0.05;
set ConsColMass [expr $deckMe*$spe/6.0+$deckMm*$spm/6.0]; # consentrated nodal mass at the top
set DampingRatio 0.06649;
set AxialLoadNode 4210;
set AxialLoad [expr -1.0*($deckWe*$spe/6.0+$deckWm*$spm/6.0)]
# Unconfined Concrete
set ConcTag 700
# uniaxialMaterial Concrete01 $matTag $fpc $epsc0 $fpcu $epsU
uniaxialMaterial Concrete01 $ConcTag $fc $epsc0 $fcu $epsu
# Confined Concrete
set ConfConcTag 800
# uniaxialMaterial Concrete01 $matTag $fpc $epsc0 $fpcu $epsU
uniaxialMaterial Concrete01 $ConfConcTag $fcc $epsc0c $fcuc $epsuc
# Reinforceing Steel
set ReinTag 900
# uniaxialMaterial Steel01 $matTag $Fy $E0 $b <$a1 $a2 $a3 $a4>
uniaxialMaterial Steel01 $ReinTag $fy $Es $hr
# Defining fiber cross-section
set ColSecTag 100
section Fiber $ColSecTag {
# Confined concrete layers in the core
set ConfConc_numSubdivCirc 12
set ConfConc_numSubdivRad 8
set ConfConc_yCenter 0.0
set ConfConc_zCenter 0.0
set ConfConc_intRad 0.0
set ConfConc_extRad [expr $colrad-$colcov]
set ConfConc_startAng 0.0
set ConfConc_endAng 360.0
# patch circ $matTag $numSubdivCirc $numSubdivRad $yCenter $zCenter $intRad $extRad <$startAng $endAng>
patch circ $ConfConcTag $ConfConc_numSubdivCirc $ConfConc_numSubdivRad $ConfConc_yCenter $ConfConc_zCenter $ConfConc_intRad $ConfConc_extRad $ConfConc_startAng $ConfConc_endAng
# Unconfined concrete layers in cover
set Conc_numSubdivCirc 12
set Conc_numSubdivRad 2
set Conc_yCenter 0.0
set Conc_zCenter 0.0
set Conc_intRad [expr $colrad-$colcov]
set Conc_extRad $colrad
set Conc_startAng 0.0
set Conc_endAng 360.0
# patch circ $matTag $numSubdivCirc $numSubdivRad $yCenter $zCenter $intRad $extRad <$startAng $endAng>
patch circ $ConcTag $Conc_numSubdivCirc $Conc_numSubdivRad $Conc_yCenter $Conc_zCenter $Conc_intRad $Conc_extRad $Conc_startAng $Conc_endAng
# Reinforcement leyers
set St_numFiber 12
set St_areaFiber $Rein_As
set St_yCenter 0.0
set St_zCenter 0.0
set St_radius [expr $colrad-$colcov]
set St_startAng 0.0
set St_endAng 330.0
# layer circ $matTag $numFiber $areaFiber $yCenter $zCenter $radius <$startAng $endAng>
layer circ $ReinTag $St_numFiber $St_areaFiber $St_yCenter $St_zCenter $St_radius $St_startAng $St_endAng
}
set large 1.0e10;
#======================================== Foundation ============================================#
set HFoundTag 110
set RFoundTag 120
uniaxialMaterial Elastic $HFoundTag $founkh
uniaxialMaterial Elastic $RFoundTag $founkrh
set ColTrsTag 2;
geomTransf Linear $ColTrsTag;
node 4301 0.0 0.0;
node 4302 0.0 0.0;
node 4303 0.0 546.0;
node 4304 0.0 852.66;
node 4305 0.0 1159.32;
node 4306 0.0 1465.98;
node 4307 0.0 1772.64;
node 4308 0.0 2079.3;
node 4309 0.0 3612.66;
node 4310 0.0 5146.0;
node 4210 0.0 5756.0;
fix 4301 1 1 1;
fix 4302 0 1 0;
#========================================= Nodal Mass ==========================================#
puts "Mass= $ConsColMass"
set MassNode 4210
if {$massDirn == 2} {
mass $MassNode $ConsColMass $ConsColMass 0.0;
} else {
mass $MassNode $ConsColMass [expr $ConsColMass/1e12] 0.0;
}
# element zeroLength $eleTag $iNode $jNode -mat $matTag1 $matTag2 -dir $dir1 $dir2 ...<-doRayleigh $rFlag>
element zeroLength 4301 4301 4302 -mat $HFoundTag $RFoundTag -dir 1 6;
# element elasticBeamColumn $eleTag $iNode $jNode $A $E $Iz $transfTag <-mass $massDens>
element elasticBeamColumn 4302 4302 4303 $large $large $large $ColTrsTag
# element dispBeamColumn $eleTag $iNode $jNode $numIntgrPts $secTag $transfTag <-mass $massDens> <-integration $intType>
element dispBeamColumn 4303 4303 4304 5 $ColSecTag $ColTrsTag
element dispBeamColumn 4304 4304 4305 5 $ColSecTag $ColTrsTag
element dispBeamColumn 4305 4305 4306 5 $ColSecTag $ColTrsTag
element dispBeamColumn 4306 4306 4307 5 $ColSecTag $ColTrsTag
element dispBeamColumn 4307 4307 4308 5 $ColSecTag $ColTrsTag
element dispBeamColumn 4308 4308 4309 5 $ColSecTag $ColTrsTag
element dispBeamColumn 4309 4309 4310 5 $ColSecTag $ColTrsTag
element elasticBeamColumn 4310 4310 4210 $large $large $large $ColTrsTag
puts "#dirn mass acting: $massDirn solverType: $solverType eigenvalue: [eigen $solverType 1]"
}
}
Re: Eigen Value Analysis
your problem is numerical and related to the variable large. large cannot be too large otherwise relative to the other terms as it effects the results due to conditioning.
e.g. for the elastic element the axial stiffness term is: A*E/L and 1e10*1e10/L is too large!
e.g. for the elastic element the axial stiffness term is: A*E/L and 1e10*1e10/L is too large!
-
amkopensees
- Posts: 24
- Joined: Mon Jun 18, 2012 9:29 am
- Location: Lehigh University
Re: Eigen Value Analysis
Thanks a lot.
-
newarrived
- Posts: 14
- Joined: Sat Aug 16, 2014 11:14 am
- Location: Iran, Shiraz
Re: Eigen Value Analysis
hello, when I use [eigen 5] in 5-story 3D concrete frame included fiber section, it's modes are same as SAP result, but when i use [eigen 1] i get this error: "ArpackSolver::Maximum number of iteration reached." and about [eigen 9] or greater number i get this error:
"ArpackSolver::Error with _saupd info = -9999
Could not build an Arnoldi factorization.IPARAM(5) the size of the current Arnoldi factorization: is 15factorization. The user is advised to check thatenough workspace and array storage has been allocated. "
in order to use [eigen -fullGenLapack] i get completely different results.
what is the problem?
"ArpackSolver::Error with _saupd info = -9999
Could not build an Arnoldi factorization.IPARAM(5) the size of the current Arnoldi factorization: is 15factorization. The user is advised to check thatenough workspace and array storage has been allocated. "
in order to use [eigen -fullGenLapack] i get completely different results.
what is the problem?
Re: Eigen Value Analysis
do you have duplicate eigenvalues in the results when you ask for 5?
-
newarrived
- Posts: 14
- Joined: Sat Aug 16, 2014 11:14 am
- Location: Iran, Shiraz
Re: Eigen Value Analysis
yes, i got the same eigenvalues for 2 first modes, and so on for 4th and 5th modes
Re: Eigen Value Analysis
Arpack is an iteractive solver .. it can't determine the 1'st eigenvalue as it keeps iterating between 1 and 2 as they are the same .. hence eigen 1 will not work, eigen 2 should.