Problem with damping force

[amkopensees]

Hi,

I have a problem with the damping force I get from my analysis. I have a simple cantilever column which is a concrete column modeled with fiber section. The mass is applied to the top node.
My problem is when I check the damping forces at different nodes, other than the base node (fixed node at the bottom), there are very abrupt increases in the forces which does not make sense to me. This column is taken out of a larger bridge model that I am currently working on. I have the same problem in the bridge model as well. The following in the column model.
I would be thankful if you can help me with this.

wipe
model BasicBuilder -ndm 2 -ndf 3

set spe 12.2e3; # length of the end spans in mm
set spm 24.4e3; # length of the mid-span in mm
set deckWe 0.039; # unit weight if the end span deck in KN/mm
set deckWm 0.052; # unit weight if the mid span deck in KN/mm
set founh 1092.0; # height of foundation in mm
set benth 1220.0; # height of the multi-column beam bent in mm
set ncol 3.0; # number of columns in each bent
set colh 4.6e3; # height of the columns in mm
set coldiam 914.4; # column diameter in mm
set colrad [expr $coldiam/2.0]; # column radius in mm
set colcov 50.8; # concrete cover for columns in mm

set deckMe [expr $deckWe/9.98e3]; # mass density of the deck of the end span in KN.Sec^2/mm^2
set deckMm [expr $deckWm/9.98e3]; # mass density of the deck of the mid span in KN.Sec^2/mm^2

#================================== Material Properties ========================================#

set fy 0.414; # fy steel yield stress in KN/mm^2
set fc -0.0276; # fc concrete compressive strength in KN/mm^2

# ###### Unconfined Concrete Properties
set epsc0 -2.0e-3; # strain at f'c
set fcu 0.0; # crushing strength in KN/mm^2
set epsu -0.0055; # strain at crushing

# Reinforcing Steel Properties
set Es 200.0; # Es steel modulus of elasticity in KN/mm^2
set hr 0.018; # hr hardening ratio
set Rein_As 645.0; # #9 bar cross-sectional area in mm^2

# ###### Confined Concrete Properties
set rho_s 1.031e-3; # look at the equation (C.10)
set K_confined [expr 1.0+($rho_s*$fy)/(-1.0*$fc)];

set fcc [expr $K_confined*$fc]; # f'cc confined concrete strength in KN/mm^2
set epsc0c [expr $K_confined*$epsc0]; # strain at f'cc
set fcuc -8.0e-3; # crushing strength in KN/mm^2
set epsuc -0.05; # strain at crushing

#===================================== Column Section ==========================================#

# 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
}

# ====================================== Transformation =========================================#

set ColTrsTag 2;
geomTransf Linear $ColTrsTag;

#========================================= Nodes ==============================================#

set DeckLevel [expr $founh/2.0+$colh+$benth];

# Left column nodes
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 4000 0.0 $DeckLevel -mass [expr ($spm*$deckMm/(2.0*$ncol))+($spe*$deckMe/(2.0*$ncol))] [expr ($spm*$deckMm/(2.0*$ncol))+($spe*$deckMe/(2.0*$ncol))] 0.0;

#========================================== Supports ===========================================#
fix 4302 1 1 1;
#========================================== Elements ==========================================#

set colfound_rFlag 1;
set bearing_rFlag 1;
set abut_rFlag 1;
set pound_rFlag 0;
# Left bent column elements
# element forceBeamColumn $eleTag $iNode $jNode $numIntgrPts $secTag $transfTag <-mass $massDens> <-integration $intType>
element forceBeamColumn 4302 4302 4303 5 $ColSecTag $ColTrsTag
element forceBeamColumn 43031 4303 4304 5 $ColSecTag $ColTrsTag
element forceBeamColumn 43041 4304 4305 5 $ColSecTag $ColTrsTag
element forceBeamColumn 43051 4305 4306 5 $ColSecTag $ColTrsTag
element forceBeamColumn 43061 4306 4307 5 $ColSecTag $ColTrsTag
element forceBeamColumn 43071 4307 4308 5 $ColSecTag $ColTrsTag
element forceBeamColumn 43081 4308 4309 5 $ColSecTag $ColTrsTag
element forceBeamColumn 43091 4309 4310 5 $ColSecTag $ColTrsTag
element forceBeamColumn 4310 4310 4000 5 $ColSecTag $ColTrsTag

#=========================================== Damping ===========================================#

# ------------ 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 DampingRatio 0.05; # damping ratio
set xDamp $DampingRatio; # 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 2; # mode 2
set T {};
set pi 3.141593;
set lambdaN [eigen -fullGenLapack [expr $nEigenJ]]; # eigenvalue analysis for nEigenJ modes
foreach lam $lambdaN {
lappend T [expr (2*$pi)/sqrt($lam)];
}

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)]; # M-prop. damping; D = alphaM*M
set betaKcurr [expr $KcurrSwitch*2.*$xDamp/($omegaI)]; # current-K; +beatKcurr*KCurrent
set betaKcomm [expr $KcommSwitch*2.*$xDamp/($omegaI)]; # last-committed K; +betaKcomm*KlastCommitt
set betaKinit [expr $KinitSwitch*2.*$xDamp/($omegaI)]; # initial-K; +beatKinit*Kini
rayleigh $alphaM $betaKcurr $betaKinit $betaKcomm; # RAYLEIGH damping

#===================== AXIAL LOAD ========================
pattern Plain 1 "Constant" {
load 4000 0.0 [expr -1.0*(($spm*$deckWm/(2.0*$ncol))+($spe*$deckWe/(2.0*$ncol)))] 0.0
}

integrator LoadControl 0.1
system BandGeneral
set tol 0.0001
set iter 1000
test NormDispIncr $tol $iter
numberer RCM;
constraints Penalty 1.0e15 1.0e15
algorithm Newton
analysis Static
analyze 10
puts "right after gravity analysis"

recorder Node -file DampingForce.out -time -node 4000 4302 4303 4304 4305 4306 4307 4308 4309 4310 -dof 1 rayleighForces

loadConst -time 0.0
# ============================================================================
recorder display "OpenSees Real Time" 300 0 1000 1000 -wipe
prp 0 2575 1000;
vup 0 1 0;
vpn -1.5 1.5 1;
viewWindow -35000 35000 -35000 35000
# display [$modeshape (-1=1st, -2=2nd..)] [$magnifying factor for nodes] [$magnifying factor for displacement]
display 1 5 100
#===================================== Time History Analysis ==================================#
set nbrModes 1
set AccelSeriesTag 2;
set DTime 0.01
set NPTS 4096
set endTime [expr ($NPTS-1.0)*$DTime];
# timeSeries Path $tag -dt $dt -filePath $filePath <-factor $cFactor>
timeSeries Path $AccelSeriesTag -dt $DTime -filePath NIS000.AT2 -factor 9.98e3;
# pattern UniformExcitation $patternTag $dir -accel $tsTag
pattern UniformExcitation 2 1 -accel $AccelSeriesTag;

wipeAnalysis;
set tol 0.0001
set iter 1000
constraints Penalty 1.0e15 1.0e15
system BandGeneral
test NormDispIncr $tol $iter
numberer RCM
algorithm Newton
integrator Newmark 0.5 0.25;
analysis Transient
analyze $NPTS $DTime

wipe

[fmk]

try a Plain constraint handler

[amkopensees]

Dear fmk,

I tried Plain constraint and still I have the same problem!

[fmk]

sorry for delay in responding i have been on holiday with only occasonal internet access .. did you look at the nodal velocities to see if there are corresponding spikes?

[amkopensees]

Dear fmk,

Its absolutely fine. Hope you had a nice holiday.
I checked the velocities at different nodes. I think they look reasonable. For example, I looked at the base node which is fixed in all DOFs. It shows zero velocity (as it should be), but there are those spikes in the damping force.
The velocity in all nodes looks reasonable. But the damping force shows the spikes, other than the top node.
The only node whose damping force does not show those spikes is the top node of the column to which the concentrated mass is attached.