Not seeing difference between -ele and -eleOnly using Region command and element mass density

[jfhuang]

Dear community,

I have created a 1story-1bay frame model and tried to assign rayleigh damping (mass and initial stiffness proportional) to the whole model using the region command. According to the manual there are several options available in the region command including the '-ele' and '-eleOnly' options. If I understand it well (please correct me if otherwise), using the '-eleOnly' option will exclude the nodal masses and hence only result in stiffness proportional damping while the '-ele' option will include the nodal mass proportional damping as intended. However, when testing these two options, I would get different results (nodal disp and accel time histories) when the mass is directly assigned to nodes, and would get the same results when I assign mass to the model using element mass density which I do not quite understand. Shouldn't it also be different when assigning mass to the model using element mass density? Or will the mass assigned by element mass density be internally considered irrespective of the '-ele' and '-eleOnly' options when forming the damping matrix? Any thoughts shall be greatly appreciated, thanks!

My script is as follows:

################################################################################
# build model
################################################################################
# units and consts
set pi 3.141593;
set in 1.;
set kip 1.;
set sec 1.;
set ft [expr 12.*$in];
set mm [expr $in/25.4];
set cm [expr 10.*$mm];
set meter [expr 100.*$cm];
set meter2 [expr $meter*$meter];
set N [expr $kip/4448.22];
set kN [expr 1000.*$N];
set kg [expr $N*$sec*$sec/$meter];
set Pa [expr $N/$meter/$meter];

# building model is built with 3 dimension, 6 dofs, with the dofs out-of-plane restrained for 2D behavior
model basic -ndm 3 -ndf 6 ;
set transfTag 1;
geomTransf Corotational $transfTag 0 1 0;

# variables
set Length [expr 30.0*$ft/$meter]; ## define bay width (unit: in) ##
set H [expr 13.0*$ft/$meter]; ## story height ##
set L [expr 30.0*$ft/$meter]; ## bay width ##
set N1 4; ## number of column elements at each story ####
set N2 [expr $N1]; ## number of beam elements at each bay ####
set L1 [expr $H/$N1]; ## define element length in column unit: in ###
set L2 [expr $L/$N2]; ## define element length in beam unit: in ###

set Nelecol [expr 1*$N1]; ## Total number of elements per column
set Nelebeam [expr 1*$N2]; ## Total number of elements per beam
set Nnodecol [expr 1*$N1+1]; ## Total number of nodes per column
set Nnodebeam [expr 1*($N2-1)]; ## Total number of internal nodes per beam
set Nnodecoltotal [expr 2*$Nnodecol]; ## Total column nodes for all columns.

set Nnode 3234; # Start node number for the frame
set Nele 2600; # Start element number for the frame
set x0 68.576; # offset in x direction
set y0 36.572; # offset in y direction
set z0 96; # offset in z direction

# nodes
for {set i 1} {$i<[expr 1*$N1+2]} {incr i 1} {
node [expr $Nnode+$i] [expr $Length*0.0+$x0] [expr 0.0+$y0] [expr $L1*($i-1)+$z0];
node [expr $Nnode+$i+1*$N1+1] [expr $Length*1.0+$x0] [expr 0.0+$y0] [expr $L1*($i-1)+$z0];
}

for {set i 1} {$i<[expr $N2]} {incr i 1} {
node [expr $Nnode+$i+(1*$N1+1)*2] [expr $Length*0.0+$L2*$i+$x0] [expr 0.0+$y0] [expr $H+$z0];
}

# node tags for column base
set nodeColumnBase {}
for {set i 1} {$i<3} {incr i 1} {
lappend nodeColumnBase [expr $Nnode+1+($i-1)*$Nnodecol]
}

# fix column base nodes for all dofs; fix all other building nodes in global uy, rx, rz for 2D system
for {set i 0} {$i < $Nnodecoltotal+1*$Nnodebeam} {incr i} {
set currNodeTag [expr $i+1+$Nnode]
if {[lsearch $nodeColumnBase $currNodeTag] != -1} {
fix $currNodeTag 1 1 1 1 1 1;
puts base=$currNodeTag
} else {
fix $currNodeTag 0 1 0 1 0 1
}
}


# material
uniaxialMaterial Steel01 1 [expr 65.0/$Pa] [expr 29000.0/$Pa] [expr 1.0/400.0];
uniaxialMaterial Steel01 2 [expr 50.0/$Pa] [expr 29000.0/$Pa] [expr 1.0/400.0];

# fiber sections
set d1 [expr 29.0 /$meter ];
set d11 [expr 12.73 /$meter ];
set tw1 [expr 0.98 /$meter ];
set bf1 [expr 14.3 /$meter ];
set GJ1 [expr 716540.0/$N/$meter2];

section Fiber 1 -GJ $GJ1 {
patch rect 1 1 11 [expr -$d1/2.0] [expr -$bf1/2.0] [expr -$d11] [expr $bf1/2.0];
patch rect 1 10 1 [expr -$d11] [expr -$tw1/2.0] [expr $d11] [expr $tw1/2.0];
patch rect 1 1 11 [expr $d11] [expr -$bf1/2.0] [expr $d1/2.0] [expr $bf1/2.0];
}

set d2 [expr 29.6 /$meter ];
set d21 [expr 12.71 /$meter ];
set tw2 [expr 1.16 /$meter ];
set bf2 [expr 14.4 /$meter ];
set GJ2 [expr 1171600.0/$N/$meter2];

section Fiber 2 -GJ $GJ2 {
patch rect 1 1 11 [expr -$d2/2.0] [expr -$bf2/2.0] [expr -$d21] [expr $bf2/2.0];
patch rect 1 10 1 [expr -$d21] [expr -$tw2/2.0] [expr $d21] [expr $tw2/2.0];
patch rect 1 1 11 [expr $d21] [expr -$bf2/2.0] [expr $d2/2.0] [expr $bf2/2.0];
}

set d3 [expr 30.7 /$meter ];
set d31 [expr 14.16 /$meter ];
set tw3 [expr 0.71 /$meter ];
set bf3 [expr 15.0 /$meter ];
set GJ3 [expr 243600.0/$N/$meter2];

section Fiber 3 -GJ $GJ3 {
patch rect 2 1 11 [expr -$d3/2.0] [expr -$bf3/2.0] [expr -$d31] [expr $bf3/2.0];
patch rect 2 10 1 [expr -$d31] [expr -$tw3/2.0] [expr $d31] [expr $tw3/2.0];
patch rect 2 1 11 [expr $d31] [expr -$bf3/2.0] [expr $d3/2.0] [expr $bf3/2.0];
}

# elements
set numIntgrPts 3;
set intType Lobatto;
for {set i 1} {$i<[expr 1*$N1+1]} {incr i 1} {
element dispBeamColumn [expr $Nele+$i] [expr $Nnode+$i] [expr $Nnode+$i+1] $numIntgrPts 1 1 -integration $intType;
element dispBeamColumn [expr $Nele+$i+1*$N1] [expr $Nnode+$i+1*$N1+1] [expr $Nnode+$i+1*$N1+1+1] $numIntgrPts 2 1 -integration $intType;
}

set a [expr 1*$N1*2]; ### total number of elements for all columns ######
set b [expr 1*$N1+1]; ### total number of nodes for each column ########
set c [expr 2*(1*$N1+1)]; ### total number of nodes for all column ########
set d [expr 1*$N2]; ### total number of elements for beam at the same level ########
set e [expr 1*($N2-1)]; ### total number of internal nodes for beam at the same level ########

set m1 [expr 100*0.001955839/$kg*$meter]; ### line mass of beam at second floor ###

element dispBeamColumn [expr $Nele+$a+1] [expr $Nnode+1+$N1] [expr $Nnode+$c+1] $numIntgrPts 3 1 -mass $m1 -integration $intType;
element dispBeamColumn [expr $Nele+$a+2] [expr $Nnode+$c+1] [expr $Nnode+$c+2] $numIntgrPts 3 1 -mass $m1 -integration $intType;
element dispBeamColumn [expr $Nele+$a+3] [expr $Nnode+$c+2] [expr $Nnode+$c+3] $numIntgrPts 3 1 -mass $m1 -integration $intType;
element dispBeamColumn [expr $Nele+$a+4] [expr $Nnode+$c+3] [expr $Nnode+$c] $numIntgrPts 3 1 -mass $m1 -integration $intType;

# uncomment to test when assigning nodal mass directly
# mass 3236 500000 500000 500000 0 0 0
# mass 3237 500000 500000 500000 0 0 0
# mass 3238 500000 500000 500000 0 0 0
# mass 3241 500000 500000 500000 0 0 0
# mass 3242 500000 500000 500000 0 0 0
# mass 3243 500000 500000 500000 0 0 0

# damping
set MpropSwitch 1.0;
set KcurrSwitch 0.0;
set KcommSwitch 0.0;
set KinitSwitch 1.0;


set Lambda [eigen 4]
set T1 [expr 4*asin(1.0)/pow([lindex $Lambda 0],0.5)];
set T2 [expr 4*asin(1.0)/pow([lindex $Lambda 2],0.5)];
puts "T1= $T1"
puts "T2= $T2"

set w1 [expr 2*$pi/$T1];
set w2 [expr 2*$pi/$T2];
set xDamp5 0.03; # damping ratio
set alphaM5 [expr $MpropSwitch*$xDamp5*(2*$w1*$w2)/($w1+$w2)]; # M-prop. damping; D = alphaM*M
set betaKcurr5 [expr $KcurrSwitch*2.0*$xDamp5/($w1+$w2)]; # current-K; +beatKcurr*KCurrent
set betaKcomm5 [expr $KcommSwitch*2.0*$xDamp5/($w1+$w2)]; # last-committed K; +betaKcomm*KlastCommitt
set betaKinit5 [expr $KinitSwitch*2.0*$xDamp5/($w1+$w2)]; # initial-K; +beatKinit*Kini

set region5 []
for {set i 1} {$i < 2*$N1+$N2+1} {incr i} {
lappend region5 [expr $i+$Nele]
}
# assign damping to all building elements
eval "region 5 -eleOnly [concat $region5] -rayleigh $alphaM5 $betaKcurr5 $betaKinit5 $betaKcomm5;"
# eval "region 5 -ele [concat $region5] -rayleigh $alphaM5 $betaKcurr5 $betaKinit5 $betaKcomm5;"

################################################################################
# recorder
################################################################################
set iNode {}
for {set i 0} {$i <= 3} {incr i} {
lappend iNode [expr $Nnode+$i+1]
}

eval "region 7 -nodeOnly $iNode"
file mkdir osRawResult
recorder Node -xml osRawResult/floorDisp.xml -time -region 7 -dof 1 disp
recorder Node -xml osRawResult/floorAcce.xml -time -region 7 -dof 1 accel


################################################################################
# analysis
################################################################################
set dt 2.256940091e-03
set numSteps 7000

# create timeSeries and load pattern
set motionFolder "."
# timeSeries Path 1 -dt $dt -filePath center_x.txt -factor 1.0;
timeSeries Path 1 -dt $dt -filePath $motionFolder/14365_acce_x.txt -factor 1.0;
# pattern UniformExcitation $patternTag $dir -accel $tsTag <-vel0 $vel0> <-fact $cFactor>
pattern UniformExcitation 2 1 -accel 1;

system ProfileSPD
constraints Plain
numberer RCM
test NormDispIncr 1.0e-7 2000 0
algorithm Newton
integrator Newmark 0.5 0.25
analysis Transient

for {set i 1} {$i <= $numSteps} {incr i} {
puts "Step $i / $numSteps"
analyze 1 $dt
}

[selimgunay]

Are you using consistent mass matrix for the distributed mass case? That could be the reason of the differences. How do the periods compare in the two cases?

[jfhuang]

Hi selimgunay,

Thank you very much for your reply. I have only used dispBeamColumn in the model and haven't used the -cMass option so I think the model has lumped mass. Both of the two cases (assign masses directly to the nodes or assign distributed masses) have the same periods (e.g. T1= 0.2303070539808556 s) while when I use the -cMass option the model has a different T1 of 0.2303453273130979 s, which sort of confirms the lumped mass assumption.

I have simplified the frame model to only have 4 nodes and 3 elements, and the mass was only applied to element 2 which is the beam of the portal frame. It would be great if you could test it out. Thanks!

Simplified script is as follows:

Code: Select all

################################################################################
# build model
################################################################################
model basic -ndm 3 -ndf 6 ;
geomTransf Corotational 1 0 1 0;

# nodes
node 1 0. 0. 0.
node 2 0. 0. 10. 
node 3 2. 0. 10.
node 4 2. 0. 0.

fix 1 1 1 1 1 1 1
fix 4 1 1 1 1 1 1

# material
uniaxialMaterial Steel01 1 [expr 65.0] [expr 29000.0] [expr 1.0/400.0];

# fiber sections
set d1  [expr 29.0 ];
set d11 [expr 12.73];    
set tw1 [expr 0.98 ];
set bf1 [expr 14.3 ];
set GJ1 [expr 716540.0];

section Fiber 1 -GJ $GJ1 {
patch rect 1 1 11 [expr -$d1/2.0] [expr -$bf1/2.0]  [expr -$d11]   [expr $bf1/2.0];
patch rect 1 10 1 [expr -$d11]    [expr -$tw1/2.0]  [expr $d11]    [expr $tw1/2.0];
patch rect 1 1 11 [expr $d11]     [expr -$bf1/2.0]  [expr $d1/2.0] [expr $bf1/2.0];
}

# elements
set numIntgrPts 3;
set intType Lobatto;
set m1 [expr 100]; ### line mass of beam at second floor ###

element dispBeamColumn   1         1         2        $numIntgrPts      1      1             -integration    $intType;
element dispBeamColumn   3         3         4        $numIntgrPts      1      1             -integration    $intType;

# element dispBeamColumn   2         2         3        $numIntgrPts      1      1             -integration    $intType;
element dispBeamColumn   2         2         3        $numIntgrPts      1      1  -mass $m1  -integration    $intType;
# element dispBeamColumn   2         2         3        $numIntgrPts      1      1  -mass $m1  -cMass -integration    $intType;

# uncomment to test when assigning nodal mass directly
# set nodeMass [expr $m1*2./2.]
# mass 2 $nodeMass $nodeMass $nodeMass 0 0 0
# mass 3 $nodeMass $nodeMass $nodeMass 0 0 0

# damping
set MpropSwitch 1.0;
set KcurrSwitch 0.0;
set KcommSwitch 0.0;
set KinitSwitch 1.0;

set Lambda [eigen 2]
set T1 [expr 4*asin(1.0)/pow([lindex $Lambda 0],0.5)];
set T2 [expr 4*asin(1.0)/pow([lindex $Lambda 1],0.5)];
puts "T1= $T1"
puts "T2= $T2"

set pi 3.1415926
set w1 [expr 2*$pi/$T1];
set w2 [expr 2*$pi/$T2];
set xDamp5 0.03;					# damping ratio
set alphaM5 [expr $MpropSwitch*$xDamp5*(2*$w1*$w2)/($w1+$w2)];	# M-prop. damping; D = alphaM*M
set betaKcurr5 [expr $KcurrSwitch*2.0*$xDamp5/($w1+$w2)];         		# current-K;      +beatKcurr*KCurrent
set betaKcomm5 [expr $KcommSwitch*2.0*$xDamp5/($w1+$w2)];   		      # last-committed K;   +betaKcomm*KlastCommitt
set betaKinit5 [expr $KinitSwitch*2.0*$xDamp5/($w1+$w2)];         	      # initial-K;     +beatKinit*Kini

# assign damping to all building elements
set region1 [list 1 2 3];
# eval "region 1 -eleOnly [concat $region1] -rayleigh $alphaM5 $betaKcurr5 $betaKinit5 $betaKcomm5;"
eval "region 1 -ele [concat $region1] -rayleigh $alphaM5 $betaKcurr5 $betaKinit5 $betaKcomm5;"

################################################################################
# recorder
################################################################################
set iNode [list 2 3]

eval "region 2 -nodeOnly $iNode"
file mkdir osRawResult
recorder Node -xml osRawResult/floorDisp.xml -time -region 2 -dof 1 disp
recorder Node -xml osRawResult/floorAcce.xml -time -region 2 -dof 1 accel

################################################################################
# analysis
################################################################################
set dt 2.256940091e-03
set numSteps 7000

# create timeSeries and load pattern
set motionFolder "."
# timeSeries Path 1 -dt $dt -filePath center_x.txt -factor  1.0;
timeSeries Path 1 -dt $dt -filePath $motionFolder/14365_acce_x.txt -factor  1.0;
# pattern UniformExcitation $patternTag $dir -accel $tsTag <-vel0 $vel0> <-fact $cFactor>
pattern UniformExcitation 2  1 -accel 1;

system ProfileSPD
constraints Plain
numberer RCM
test NormDispIncr 1.0e-7 2000  0
algorithm Newton
integrator Newmark 0.5 0.25
analysis Transient

for {set i 1} {$i <= $numSteps} {incr i} {
	puts "Step $i / $numSteps"
	analyze 1 $dt
}