Possible numeric error using Two Node Link element or any other object in the model

[wless]

Hi:

We are working on a research that consists in evaluating the effect of live load on nonlinear response history analysis of low-rise buildings. As part of the research program, we are executing a parametric study of "idealized structures". These structures are based on Two-Node Link Elements with a trilinear plasticity model using Multilinear Material. These elements represent a lateral force-resisting system of "real" steel buildings.

The following picture shows a geometric configuration of this type of structures. The two-story building has two links acting on each orthogonal direction and has a master point that contraints these links.

https://ibb.co/cg0FxHS

The links have trilinear constitutive model like show below:

https://ibb.co/XFMy344

In a calibration phase, we are comparing a response on multiple platforms (SAP2000, MATLAB and OpenSees). In this process, we have seen OpenSees has possibly a numerical problem because somecases its response doesn't coincide with others (MATLAB and SAP2000).

The following figure shows a response of hysteresis curve for the three different softwares:

https://ibb.co/d57MKzj

It can be seen that response of OpenSees is different with others. No matter if the response no coincide it shouldn't be out of backbone curve.

Here it's a model code on OpenSees:
_____________________________________________________________________________________________________________

# Regular idealized structure, R = 6
# -------------------------------------------------------------------

# Units: kN, m, s

# -----------------------------------
# Modeling
# -----------------------------------

# Remove an existing model
wipe

# Create ModelBuilder object
# del BasicBuilder -ndm $ndm
model BasicBuilder -ndm 3

# Variables
set Lx 46.8; # Structure length in x direction - m
set Lz 46.8; # Structure length in z direction - m
set hz 3.9; # Story height - m
set e1_1 0.00851750055232987; # Shear deformation at yield limit state for first story - m
set e2_1 0.0298112519331545; # Shear deformatión at ultimate limit state for first story - m
set e3_1 2.98112519331545; # Shear deformatión at rupture limit state for first story - m
set s1_1 10923.0950780772; # Shear force at yield limit state for first story - kN
set s2_1 21846.1901561544; # Shear force at ultimate limit state for first story - kN
set s3_1 21846.1901561544; # Shear force at rupture limit state for first story - kN
set e1_2 0.00530821365771109; # Shear deformation at yield limit state for second story - m
set e2_2 0.0185787478019888; # Shear deformatión at ultimate limit state for second story - m
set e3_2 1.85787478019888; # Shear deformatión at rupture limit state for second story - m
set s1_2 6807.41047466863; # Shear force at yield limit state for second story - kN
set s2_2 13614.8209493373; # Shear force at rupture limit state for second story - kN
set s3_2 13614.8209493373; # Shear force at rupture limit state for second story - kN
set mx 2233.42323831278; # Mass in x direction - Mg
set mz 2233.42323831278; # Mass in z direction - Mg
set mt 815288.818913696; # Rotation mass about y direction - Mg-rad
set psi 0.03; # damping coefficient
set dt 0.01; # Time step - s
set nSign 2221; # Number of time steps
set g 9.80665; # Gravity accelearion - m/s^2
set alpha 0.0; # "Coefficient of irregularity
set beta 0.0; # "Coefficient of irregularity
set gamma 0.0; # "Coefficient of irregularity
set delta 0.0; # "Coefficient of irregularity

# Create nodes
# de $nodeTag (ndm $coords)
node 1 [expr $Lx/2] 0 [expr $alpha*$Lz]
node 2 [expr $Lx/2] 0 [expr $Lz - $beta*$Lz]
node 3 [expr $gamma*$Lx] 0 [expr $Lz/2]
node 4 [expr $Lx - $delta*$Lx] 0 [expr $Lz/2]
node 5 [expr $Lx/2] $hz [expr $alpha*$Lz]
node 6 [expr $Lx/2] $hz [expr $Lz - $beta*$Lz]
node 7 [expr $gamma*$Lx] $hz [expr $Lz/2]
node 8 [expr $Lx - $delta*$Lx] $hz [expr $Lz/2]
node 9 [expr $Lx/2] $hz [expr $Lz/2]
node 10 [expr $Lx/2] [expr 2*$hz] [expr $alpha*$Lz]
node 11 [expr $Lx/2] [expr 2*$hz] [expr $Lz - $beta*$Lz]
node 12 [expr $gamma*$Lx] [expr 2*$hz] [expr $Lz/2]
node 13 [expr $Lx - $delta*$Lx] [expr 2*$hz] [expr $Lz/2]
node 14 [expr $Lx/2] [expr 2*$hz] [expr $Lz/2]

# Assign masses
# ss $nodeTag (ndf $massValues)
mass 9 $mx 0 $mz 0 $mt 0
mass 14 $mx 0 $mz 0 $mt 0

# Define trilinear plasticity model
# iaxialMaterial MultiLinear $matTag $u1 $f1 $u2 $f2 $u3 $f3
uniaxialMaterial MultiLinear 3 $e1_1 $s1_1 $e2_1 $s2_1 $e3_1 $s3_1
uniaxialMaterial MultiLinear 4 $e1_2 $s1_2 $e2_2 $s2_2 $e3_2 $s3_2

# Define elements
# ement twoNodeLink $eleTag $iNode $jNode -mat $matTags -dir $dirs -orient $y1 $y2 $y3 -doRayleigh
element twoNodeLink 1 1 5 -mat 3 -dir 2 -orient -1 0 0 -doRayleigh
element twoNodeLink 2 2 6 -mat 3 -dir 2 -orient -1 0 0 -doRayleigh
element twoNodeLink 3 3 7 -mat 3 -dir 2 -orient 0 0 1 -doRayleigh
element twoNodeLink 4 4 8 -mat 3 -dir 2 -orient 0 0 1 -doRayleigh
element twoNodeLink 5 5 10 -mat 4 -dir 2 -orient -1 0 0 -doRayleigh
element twoNodeLink 6 6 11 -mat 4 -dir 2 -orient -1 0 0 -doRayleigh
element twoNodeLink 7 7 12 -mat 4 -dir 2 -orient 0 0 1 -doRayleigh
element twoNodeLink 8 8 13 -mat 4 -dir 2 -orient 0 0 1 -doRayleigh

# Set fixes
# x $nodeTag (ndf $constrValues)
fix 1 1 1 1 1 1 1
fix 2 1 1 1 1 1 1
fix 3 1 1 1 1 1 1
fix 4 1 1 1 1 1 1
fix 5 0 1 0 1 0 1
fix 6 0 1 0 1 0 1
fix 7 0 1 0 1 0 1
fix 8 0 1 0 1 0 1
fix 9 0 1 0 1 0 1
fix 10 0 1 0 1 0 1
fix 11 0 1 0 1 0 1
fix 12 0 1 0 1 0 1
fix 13 0 1 0 1 0 1
fix 14 0 1 0 1 0 1

# Define diaphrams
# gidDiaphragm $perpDirn $masterNodeTag $slaveNodeTags
rigidDiaphragm 2 9 5 6 7 8
rigidDiaphragm 2 14 10 11 12 13

# Do eigen values (Peridos)
set lambda [eigen -fullGenLapack 6];

# Perform a frecuencies and periods of the structure
set omega {}
set f {}
set T {}
set pi [expr acos(-1.0)]

foreach lam $lambda {
lappend omega [expr sqrt($lam)]
lappend f [expr sqrt($lam)/(2*$pi)]
lappend T [expr (2*$pi)/sqrt($lam)]
}

puts "Periods are $T"

# Write a file with fundamental periods
set period "periodos.out"
set Periods [open $period "w"]
foreach t $T {
puts $Periods "$t"
}
close $Periods

# Rayleigh damping
set w1 [expr [lindex $omega 0]/1.5]
set w2 [lindex $omega 1]
set a0 [expr $psi*2*$w1*$w2/($w1 + $w2)]
set a1 [expr $psi*2/($w1 + $w2)]

# yleigh $alphaM $betaK $betaKinit $betaKcomm
rayleigh $a0 0.0 $a1 0.0

# Restart every previous analysis
reset

# Define acceleration records
# meSeries Path $tag -dt $dt -filePath $filePath -factor $cFactor
timeSeries Path 1 -dt $dt -filePath GM_01_1.txt -factor [expr 1*$g]
timeSeries Path 2 -dt $dt -filePath GM_01_2.txt -factor [expr 1*$g]

# Set patterns with previous-defined records
# ttern UniformExcitation $patternTag $dir -accel $tsTag
pattern UniformExcitation 1 1 -accel 1
pattern UniformExcitation 2 3 -accel 2

# ------------------------------------
# Analysis
# ------------------------------------

# gorithm algorithmType?
algorithm Newton
# tegrator Newmark $gamma $beta
integrator Newmark 0.5 0.25
# nstraints constraintType?
constraints Transformation
# mberer numbererType?
numberer RCM
# stem systemType?
system FullGeneral
# st NormDispIncr $tol $iter <$pFlag>
test NormDispIncr 1E-5 150
# Crear objeto analysis
analysis Transient

# -----------------------------------
# Recording
# -----------------------------------

# Nodal displacement and accelerations of CM nodes
# corder Node -file $fileName -precision $nSD -time -node $node -dof $dofs $respType
recorder Node -file displ.out -precision 10 -time -node 9 -dof 1 3 disp
recorder Node -file accel.out -precision 10 -time -node 9 -dof 1 3 accel

# Drift of CM nodes
# corder Drift -file $fileName -precision $nSD -time -iNode $node -jNode $node -dof $dof -perpDirn $perpDirn
recorder Drift -file driftX.out -precision 10 -time -iNode 3 -jNode 9 -dof 1 -perpDirn 2
recorder Drift -file driftY.out -precision 10 -time -iNode 3 -jNode 9 -dof 3 -perpDirn 2

# Hysteresis of elements (kN - m)
# corder Element -file $fileName -precision $nSD -time -ele $ele $arg1
recorder Element -file basicForce.out -precision 10 -time -ele 1 basicForce
recorder Element -file basicDisp.out -precision 10 -time -ele 1 basicDisplacement

# -------------------------------
# Execute analysis
# -------------------------------

# Ejecutar el análisis
analyze $nSign $dt

_____________________________________________________________________________________________________________


Thank you for your attention

[wless]

We have to mention there are some building configurations that show good results. Here, we upload a figure that demostrates an equal results for SAP2000, MATLAB and OpenSees. For this reason, we consider OpenSees has some type of numerical error:

https://ibb.co/RgFXZ07