domain error: argument not in valid range

[linzhidan]

dear vesna
i just created a 2D column for rubber bearing.
i just modify the TransfType from "Linear" to "PDelta" in "geomTransf $TransfType $TransfTag", when i run my input, I encounter such an error.

domain error: argument not in valid range
while executing
"expr pow($lambda,0.5)"
invoked from within
"set omega [expr pow($lambda,0.5)]"
(file "E:\OPENSEESWORK\RCcolumn_PDelte_question.tcl" line 169)

i would like to know why the "Linear" is ok, the "PDelta" is wrong.
looking forward to your reply! thanks!

my input:
# --------------------------------------------------------------------------------------------------
# column 2D for rubber bearing
# all units are in kg, m, s, N
# forceBeamColumn ELEMENT with multinode lumped mass
# consider Linear or PDelter
# ^Y
# |
# 2 ---
# | |
# | |
# | |
# (1) Lcol=0.164 m
# | |
# | |
# | |
# =1= ---- -------->X
#

# SET UP ----------------------------------------------------------------------------
wipe; # clear opensees model
model basic -ndm 2 -ndf 3; # 2 dimensions, 3 dof per node

# Create nodes
# ------------
# Set parameters for overall model geometry
set Lcol 0.164; # column length
set PI [expr 2*asin(1.0)]; # define constant
# Create nodes--------nodal coordinates
node 1 0 0; # base node
node 2 0 $Lcol; # top node

# Boundary Conditions
fix 1 1 1 1; # fix supports at base of column
fix 2 0 0 1; # fix ratation at top of column

# Define materials for nonlinear columns
# -------------------------------------------------------------------------------------------------------
# physical dimension parameter
set intRad 0.015; # isolater bearing inner radius
set extRad 0.250; # isolater bearing external radius
set nr 20; # number of rubber layer
set ns 19; # number of steel plate
set tr 0.0048; # thickness of the rubber layer
set Tr [expr $tr*$nr]; # total thickness of the rubber layer
set ts 0.002; # thickness of steel plate
set tc 0.015; # thickness of closure steel plate
set Aiso [expr $PI*(pow($extRad,2)-pow($intRad,2))]; # the cross-sectional area of isolater bearing
set Iiso [expr $PI*(pow($extRad,4)-pow($intRad,4))*1./4];# inertia moment of isolater bearing

# material parameter
set G 0.392; # standard shear modulus of the rubber material
set Ev 1960e+6; # constrained volume elastic modulus of the rubber material, unit: N/m^2
set S1 [expr 2.*($extRad-$intRad)/(4*$tr)];# first shape factor
set k [expr 0.97939+0.17734*$G-1.4516*pow($G,2)+0.86783*pow($G,3)];# Correction factor
set G2 0.392e+6; # standard shear modulus of the rubber material,unit: N/m^2
set E0 [expr 3*$G2]; # standard elastic modulus of rubber material
set Ec [expr $E0*(1+2*$k*pow($S1,2))]; # compression modulus of rubber bearing
set Eb [expr $E0*(1+2/3*$k*pow($S1,2))];# equibalent flexural elastic modulus of rubber bearing
set Ecv [expr $Ev*$Ec/($Ev+$Ec)]; # correction compress elastic modulus of rubber bearing
set Ebv [expr $Ev*$Eb/($Ev+$Eb)]; # correction flexural elastic modulus of rubber bearing
set rhos [expr 7850]; # density of steel, unit: kg/pow($m,3)
set rhor [expr 1150]; # density of rubber,
set rhoiso [expr ($rhor*$Aiso*$tr*$nr+$rhos*$Aiso*$ts*$ns+2.*$rhos*$Aiso*$tc)/($Aiso*$Lcol)];# bearing equivalent density

# calculated stiffness parameters
set EIiso [expr $Ebv*$Iiso*$Lcol/$Tr]; # EI, for moment-curvature relationship
set EAiso [expr $Ecv*$Aiso]; # EA, for axial-force-strain relationship
set GAiso [expr $G2*$Aiso*$Lcol/$Tr]; # GA, for shear-force-strain relationship

# nodal masses
set PA 10e+6; # vertical area load at top node, unit: N/m^2
set P [expr $PA*$Aiso]; # axial force
set mt [expr $P/9.8];

set massID 2
mass $massID [expr $rhoiso*$Lcol*$Aiso+$mt] 1.e-9 0.;

# Define ELEMENTS & SECTIONS
# define sections of column
set colMatTagFlex 11; # assign a tag number to the concrete column flexural behavior
set colMatTagAxial 12; # assign a tag number to the concrete column axial behavior
set colMatTagShear 13; # assign a tag number to the concrete column shear behavior
set colSecTag 1; # assign a tag number to the concrete column section tag
uniaxialMaterial Elastic $colMatTagFlex $EIiso; # linear behavior for flexure
uniaxialMaterial Elastic $colMatTagAxial $EAiso; # this is not used as a material, this is an axial-force-strain response
uniaxialMaterial Elastic $colMatTagShear $GAiso; # this is not used as a material, this is an shear-force-strain response
section Aggregator $colSecTag $colMatTagAxial P $colMatTagFlex Mz $colMatTagShear Vy; # combine axial and flexural and shear behavior into one section

# define geometric transformation
set TransfTag 1; # associate a tag to column transformation
set TransfType PDelta; # options, Linear PDelta Corotational
geomTransf $TransfType $TransfTag ;
set numIntgrPts 5; # element connectivity

# define element
set eleID 1;
set nodeI 1;
set nodeJ 2;
element forceBeamColumn $eleID $nodeI $nodeJ $numIntgrPts $colSecTag $TransfTag;


# Define gravity loads
pattern Plain 1 Linear {
set loadID 2;
load $loadID 0. -[expr $P+$rhoiso*$Lcol*$Aiso*9.8] 0.; # node FX FY MZ
}

# ------------------------------
# End of model generation
# ------------------------------



# ------------------------------
# Start of analysis generation
# ------------------------------
system BandGeneral; # Create the system of equation, a sparse solver with partial pivoting
constraints Plain; # Create the constraint handler, the transformation method
numberer Plain; # Create the DOF numberer, the reverse Cuthill-McKee algorithm
test NormDispIncr 1.0e-12 10; # Create the convergence test, the norm of the residual with a tolerance of
algorithm Newton; # Create the solution algorithm, a Newton-Raphson algorithm
integrator LoadControl 0.1; # Create the integration scheme, the LoadControl scheme using steps of 0.1
analysis Static; # Create the analysis object

# ------------------------------
# End of analysis generation
# ------------------------------


# ------------------------------
# Finally perform the analysis
# ------------------------------
analyze 10; # perform the gravity load analysis, requires 10 steps to reach the load level
set nEigen 1; # show mode shape for nEigen eigenvalue
set lambdaN [eigen -generalized -fullGenLapack $nEigen]; # perform eigenvalue analysis for ModeShape -generalized -fullGenLapack
set lambda [lindex $lambdaN [expr $nEigen-1]];
set omega [expr pow($lambda,0.5)]
set Tperiod [expr 2*$PI/$omega]; # period (sec.)

[vesna]

Print to terminal the eigenvalues for the two cases of geometric transformation you are considering to make sure the values look ok. From the error message it looks like there is a problem with eigenvalues.

[linzhidan]

dear vesna,
thank you for your reply!
i print the two cases' eigenvalues.
for Linear geometric transformation , $lambdaN =3.817169e+000,
for PDelta, $lambdaN = -5.593893e+001,
it is the problem you have decided.

my question is : i just modify the Linear to PDelta, why it like is?
waiting for your reply,thanks!

[linzhidan]

dear vesna,
my element is this:
# ^Y
# |
# 2 ---
# | |
# | |
# | |
# (1) Lcol=0.164m
# | |
# | |
# | |
# =1= ---- -------->X

[vesna]

PDelta reduces the stiffness of the system, which in nonlinear range of behavior can lead to system's negative stiffness. Negative eigenvalue implies negative stiffness of the system.

[linzhidan]

dear vesna
thank you for your reply.
i would ask another question follow this problem.
it is :
when i do the dynamic ananysis with PDelta, it will put in damping, but the $lambdaI is also negative.
how can i do for this?


# RAYLEIGH damping parameters,
# D=$alphaM*M + $betaKcurr*Kcurrent + $betaKcomm*KlastCommit + $beatKinit*$Kinitial
set xDamp 0.05; # damping ratio of concrete
set MpropSwitch 1.0;
set KcurrSwitch 0.0;
set KcommSwitch 1.0;
set KinitSwitch 0.0;
set nEigenI 1; # mode 1
set nEigenJ 1; # 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)];
puts "$omegaI"
set omegaJ [expr pow($lambdaJ,0.5)];
puts "$omegaJ"
set alphaM [expr $MpropSwitch*$xDamp*(2*$omegaI*$omegaJ)/($omegaI+$omegaJ)]; # M-prop. damping; D = alphaM*M
puts "$alphaM"
set betaKcurr [expr $KcurrSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # current-K; +beatKcurr*KCurrent
puts "$betaKcurr"
set betaKcomm [expr $KcommSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # last-committed K; +betaKcomm*KlastCommitt
puts "$betaKcomm"
set betaKinit [expr $KinitSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # initial-K; +beatKinit*Kini
puts "$betaKinit"
# define damping
rayleigh $alphaM $betaKcurr $betaKinit $betaKcomm; # RAYLEIGH damping

[vesna]

Your axial load is probably vary big. Plot shear force-displacement relationship to see when the stiffness degradation gets initiated. It can help you understand the behavior of your model better.

[linzhidan]

dear vesna,
thank you for your attention.
my model simulates the rubber bearing using columnelement with section Aggregator which combining axial and flexural and shear behavior into one section.
axial stiffness is much more bigger than the shear stiffness in the rubber bearing.

so i will try the pushover analysis on rubber bearing with PDelta using this modal parameters.
Thank you much, and see you tomorrow!