Nonlinear secant stiffness matrix

[ag110]

Dear friends
I need to obtain nonlinear secant stiffness matrix for a simple model during pushover analysis at each increment. I want to know if it is possible in opensees software? If yes, how?
If a subroutine is needed, does anyone has this subroutine?

[fmk]

unless this is a single degree of freedom model it is not possible to do in OpenSees

[els83]

I have a problem with Beamcolumn element. I modeled a simple 7DOF system with BeamColumn elements. I restricted all DOFs except the DOF 1. I expected to have the stiffness matrix of a shear building and calculated the eigen frequencies. they are not what I get from opensees as modal frequencies. But when I model the same system with axial elements with the same stiffness the modal frequencies are the same.

What is the problem, is OPENSEES modeling the the system with BeamColumn element more flexible (modal frequencies are less than what is expected).

[vesna]

I'm sorry but I do not understand your problem. Would you please explain it with more details?

[mnazari]

How can we obtain secant stiffness for a SDOF model? Since in my model, the natural frequency achieved by "Eigen 1" command is based on the tangent stiffness, not the secant. My SDOF model is an elasticBeamColumn element jointed to a zero length element with elastic multilinear behavior and I need the value of secant stiffness at each step to calculate the damping coefficient (K commit) based on that.

[vesna]

OpenSees does not output secant stiffness. You will have to record force and deformation for each time step to calculate Ksec yourself.

K commit is not secant but tangent stiffness.

[mnazari]

Thanks Vesna, I know that Kcommit is tangent stiffness, but I have a 2D rotational spring at the base of my structure and its behavior is dictated by an elastic Multilinear material. By doing free vibration analysis I have understood that plugging damping coefficient based on secant stiffness at each level as K commit, calculated damping is closer to the assigned value. But I cannot repeat it for an earthquake loading!! Do you have specific recommendation for damping coefficient in this case?

[vesna]

Some researches suggest using K tangent for the systems that experience nonlinear deformations.

[mnazari]

Hi,
Defining this multilinear elastic material as my rotational spring (mentioned in my previous post), I do Eigen value analysis to calculate frequency of the system in order to plug it in Rayleigh damping command. But achieved frequency is always the same, no matter how big is my loading, and it is based on initial stiffness defined by that multilinear link. It cannot be adjusted to the tangent stiffness at each level of loading! How can I do that?
Thanks

[vesna]

Can you post your model so I can see when do you do eigen analysis.

[mnazari]

Sure.. Here it is. Thanks


# SET UP ----------------------------------------------------------------------------
wipe; # clear opensees model
model basic -ndm 2 -ndf 3; # 2 dimensions, 3 dof per node
file mkdir Data; # create data directory

# source "Input.tcl"

# define GEOMETRY -------------------------------------------------------------
# nodal coordinates:
node 1 0.0 0; # node#, X Y
node 2 0.0 0;
node 3 0.0 164;
node 33 56.5 164;
node 44 146.07 164;
node 4 203.07 164;
node 45 203.07 119.25;
node 54 203.07 43.75;
node 5 203.07 13.3;
node 6 128.74 13.3;
node 66 128.74 13.3;
node 7 277.4 13.3;
node 77 277.4 13.3;
#Supports of the Mass Rig Frame
node 8 128.74 -86.34;
node 9 277.4 -86.34;


uniaxialMaterial ElasticMultiLinear 1 -strain -0.03 -0.02 -0.01 -0.0075 -0.0050 -0.0025 -0.0015 -0.0010 -8.55E-05 0 8.55E-05 0.0010 0.0015 0.0025 0.0050 0.0075 0.01 0.02 0.03 -stress -5638.16 -5433.4 -4653.24 -4438.2 -4233.57 -3916.43 -3811.13 -3767.87 -1418.47 0 1418.47 3767.87 3811.13 3916.43 4233.57 4438.2 4653.24 5433.4 5638.16
uniaxialMaterial Elastic 2 29000


element zeroLength 1 1 2 -mat 1 -dir 6 -doRayleigh 1

# Single point constraints -- Boundary Conditions
fix 1 1 1 1 ; # node DX DY RZ
fix 8 1 1 0 ; # node DX DY RZ
fix 9 1 1 0 ; # node DX DY RZ

equalDOF 1 2 1 2 ;

# Mass Rig Frame
equalDOF 6 66 1 2;
equalDOF 7 77 1 2;

equalDOF 5 54 1;


# nodal masses ( 60kip+ weight of wall)
mass 45 0.05176 1.e-9 0.; # node#, Mx My Mz, Mass=Weight/g.
mass 54 0.05176 1.e-9 0.; # node#, Mx My Mz, Mass=Weight/g.
mass 5 0.05176 1.e-9 0.; # node#, Mx My Mz, Mass=Weight/g.
mass 3 0.00673 1.e-9 0.; # node#, Mx My Mz, Mass=Weight/g.


# Define ELEMENTS -------------------------------------------------------------
# define geometric transformation: performs a linear geometric transformation of beam stiffness and resisting force from the basic system to the global-coordinate system
geomTransf Linear 1; # associate a tag to transformation

element elasticBeamColumn 2 2 3 375 5219.2 175781 1; # element elasticBeamColumn $eleTag $iNode $jNode $A $E $Iz $transfTag
rigidLink beam 3 33;
element truss 3 33 44 8.405 2 -doRayleigh 1
rigidLink beam 4 44;


#Stiff Section and Material
element elasticBeamColumn 4 4 45 7000 9.0E08 29.461 1; # element elasticBeamColumn $eleTag $iNode $jNode $A $E $Iz $transfTag
element elasticBeamColumn 5 45 54 7000 9.0E08 29.461 1; # element elasticBeamColumn $eleTag $iNode $jNode $A $E $Iz $transfTag
element elasticBeamColumn 6 54 5 7000 9.0E08 29.461 1; # element elasticBeamColumn $eleTag $iNode $jNode $A $E $Iz $transfTag


# IPE14 section for columns and beams; change to real sections when you have the correct size.
element elasticBeamColumn 7 5 66 2.542 29000 12.9976 1; # element elasticBeamColumn $eleTag $iNode $jNode $A $E $Iz $transfTag
element elasticBeamColumn 8 5 77 2.542 29000 12.9976 1; # element elasticBeamColumn $eleTag $iNode $jNode $A $E $Iz $transfTag
element elasticBeamColumn 9 6 8 2.542 29000 12.9976 1; # element elasticBeamColumn $eleTag $iNode $jNode $A $E $Iz $transfTag
element elasticBeamColumn 10 7 9 2.542 29000 12.9976 1; # element elasticBeamColumn $eleTag $iNode $jNode $A $E $Iz $transfTag


pattern MultipleSupport 1 {
set GMfile "DispTable4.acc"
set DispSeries "Series -dt 0.003906 -filePath $GMfile -factor -1"
groundMotion 1 Plain -disp $DispSeries
imposedMotion 1 1 1
}

set lambda [eigen 1]
# Frequency Calculation
set omega [expr pow($lambda,0.5)]
set pi 3.141592654
# Period Calculation
puts "period:"
puts "[expr (2*$pi)/pow($lambda,0.5)]"
puts "frequency(rad/s):"
puts "[expr pow($lambda,0.5)]"
# Damping
rayleigh 0. 0. 0. [expr 2*0.02/$omega] ; # set damping based on first eigen mode


# To remove an sp_Constraint
remove sp 1 1


# create the analysis
wipeAnalysis; # clear previously-define analysis parameters
constraints Penalty 1.0E15 1.0E15;
# ===================================================================
set numbererTypeDynamic Plain
numberer $numbererTypeDynamic
# ===================================================================
set systemTypeDynamic BandGeneral
system $systemTypeDynamic
# ===================================================================
set TolDynamic 1.e-8;
set maxNumIterDynamic 1000;
set printFlagDynamic 0;
set testTypeDynamic EnergyIncr
test $testTypeDynamic $TolDynamic $maxNumIterDynamic $printFlagDynamic;
# ===================================================================
set algorithmTypeDynamic KrylovNewton
algorithm $algorithmTypeDynamic;
# ===================================================================
# HHT instead of Newmark, avoids wierd increasing aacceleration at top of SIN5 model.
set integratorTypeDynamic HHT
integrator $integratorTypeDynamic 0.5
# ===================================================================
set analysisTypeDynamic Transient
analysis $analysisTypeDynamic


set DtAnalysis 0.003906;
set TmaxAnalysis 52.46410751;
set Nsteps [expr int($TmaxAnalysis/$DtAnalysis)];
set ok [analyze $Nsteps $DtAnalysis];
# ===================================================================
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 1.0e-2 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]"

[vesna]

You do the eigen analysis and define the damping before any other type of analysis. The way you have defined it, the damping will be proportional to the tangent stiffness of the system. Is that what you wanted?

[mnazari]

Yes, Exactly! but it should change due to the definition of my rotational spring during earthquake motion.

[vesna]

The way you defined it the stiffness will be updated based on the properties of your rotational spring. However, the coefficient $betaKcomm will not change during the analysis.

[mnazari]

So, is there anyway I can change damping coefficient based on current stiffness?