MVLEM issue with Dynamic Analysis

[zwaezi]

Dear OpenSEES Community,

I am reaching out to seek assistance with an issue I've encountered while conducting dynamic analysis using MVLEM (Multiple-Vertical-Line-Element-Model) element for Reinforced Concrete (RC) walls. During the dynamic analysis of a simple cantilever wall subjected to horizontal excitation, I have observed an unexpected discrepancy. Despite setting a very small NormDispIncr tolerance of 1e-12 and assuming a zero vertical mass to mitigate the influence of vertical inertial forces, the axial reaction at the base does not match the applied weight. Furthermore, I've noted that this element fails to converge when the NormUnbalance is used as test criteria for dynamic analysis except when tolerance is set to be a very high value such as 1e-2 or 1e-3.

In an effort to troubleshoot this issue, I've experimented with various parameters and configurations. Surprisingly, setting the Rayleigh damping coefficient to zero appears to resolve the problem, aligning the axial reaction at the base with the applied gravity load.

I have prepared a comparison plot illustrating the base vertical reaction versus the applied gravity load, which I have attached for reference.

Thank you for your attention to this matter. Your insights, suggestions, or guidance would be invaluable in resolving this issue.

Warm regards, Zak

[mhscott]

Hello @zwaezi, thank you for summarizing the issue. It's unclear if this is a modelling issue of but of some kind.

So, please post a minimal working example as a reply to this thread.

P.S. It's capitalized OpenSees, not OpenSEES

[zwaezi]

Thank you Dr @mhscott. I am using a cantilever wall element with just two nodes employing MVLEM. While I have assigned a near-zero mass to the vertical DOF of the upper node and applied earthquake excitation to the system, I am noticing that the vertical reaction at the base of the wall changes with time (as shown in the above figure). I am using Rayleigh K-dependent damping. However, when I set the Rayleigh damping to be zero this discrepancy between the support reaction and the applied gravity load is resolved. I am very grateful for your help.

Zak

[silviamazzoni]

Which stiffness are you using?

[mhscott]

Please post a simple script.

[zwaezi]

Which stiffness are you using?

Hello Dr. @silviamazzoni I am using Committed Stiffness.

[zwaezi]

Thank you Dr. @mhscott . Sorry for taking too long to prepare the script.

#    
# by Zakariya Waezi

#    units: N, mm, MPa
#from openseespy.opensees import *
import sys
import os
import shutil
sys.path.append('c:/programdata/anaconda3/envs/py38/lib/site-packages/')
from openseespy.opensees import *
import opsvis as opsv
#import openseespy.postprocessing.Get_Rendering as opsplt
import numpy as np
import matplotlib.pyplot as plt
from math import asin, sqrt, atan, sin, cos
import csv
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
import time

mypath= 'IDA_Nishida-Rectangular'
foldername=    mypath
#-----------------------------------------------
# Concrete Properties
#-----------------------------------------------    
g = 9810
pi=np.pi
# Values provided in N, mm
fc=30.4 # default value 36.2
fc_unconfined=fc
epc_unconfined=0.002
fcu_unconfined=6.1
ecu_unconfined=0.006
lambda_unconfined=0.1
ft_unconfined=0.36*np.sqrt(fc_unconfined)
Ec_unconfined=6900+3320*np.sqrt(fc_unconfined)
Et_unconfined=Ec_unconfined
G=0.4*Ec_unconfined
Ec=Ec_unconfined

fc_confined=34.9
epc_confined=0.0036
fcu_confined=7.2
ecu_confined=0.018
lambda_confined=0.1
ft_confined=0.36*np.sqrt(fc_confined)
Ec_confined=6900+3320*np.sqrt(fc_confined)
Et_confined=Ec_confined

fc_compr_strut=0.34*fc# default values 0.34*fc
Ec_diagonal=Ec_unconfined#6900+3320*np.sqrt(fc_compr_strut)
eps0_diagonal=2*fc_compr_strut/Ec_diagonal

#-----------------------------------------------
# Steel Properties
#-----------------------------------------------
fy=373.
fu=500.
fyt=340
Es=200000
Esh=6000.
Esh=Es*0.01   #default Es*0.01
eps_sh=0.01
eps_ult=0.15
lsr=30
beta_buckling=1.0
r_buckling=0.6
gamma_buckling=0.5
Cf_fatigure=0.26#0.26
alpha_fatigue=0.506
Cd_fatigue=0.26 #0.26   

#----------------------------------------------
# Column Properties 
#----------------------------------------------
L=3000.#.              # Total Length
bCol=800.               # Column width
hCol=450.               # Column Length
LpCol=416;              # Column Height

IgCol=1/12*bCol*hCol**3      # Column's Uncracked Moment of Inertia
I_effCol=0.4*IgCol        #Column's Effective Moment of Inertia
E_rigidCol=Es*100      # Column's Rigid modulus of elasticity for Relative member of the plastic hinge
AgCol=bCol*hCol        #Column's Area

AsTen=17.*pi/4*10.**2.
AsComp=17.*pi/4*10.**2.
AsMid=2.*pi/4*10.**2.
As_list=[AsMid for i in range(9)]
As_list[0]=AsTen
As_list[8]=AsComp

AsShear=2.*pi/4*6.**2.
db_ten=10.0 
db_comp=10.0
Shear_Spacing=75.

cover_net_Col=40.
dpCol=40+db_ten/2.
dCol=hCol-dpCol
AshearCol=bCol*dCol

Ec_Column=Ec_unconfined
G_Column=G

Axial_Compressive_stress=-1.0
Axial_Load=Axial_Compressive_stress*AgCol # Axial Load
Weight=abs(Axial_Load)
Weight_vert=Weight

Mass =  Weight/g
Mass_vert = Weight_vert/g
Rotational_Mass=1e-12

Conf_Concrete_mat_index=1#10#default value =10
Unconf_Concrete_mat_index=2
Diagonal_Concrete_mat_index=3#514#511#previous value=4----411 chang and mander---510&511 ConcreteMO
Steel_mat_index_col=4
ShearSteelMat=5#555#557#556

#Remove existing model
wipe()

if not os.path.exists(mypath):
    os.makedirs(mypath)

#    Create ModelBuilder (with two-dimensions and 3 DOF/node)
model('Basic', '-ndm', 2,'-ndf', 3)
logFile(f'{foldername}/LogFileNishida_Cantilever.txt', '-noEcho')

#======================================================================
# Define materials for nonlinear columns
# ------------------------------------------
# CONCRETE                   tag  f'c    ec0    f'cu   ecu
# ------------------------------------------
# unconfined concrete

uniaxialMaterial('Concrete02', 1, -fc_confined, -epc_confined,  -fcu_confined,  -ecu_confined, lambda_confined, ft_confined, Et_confined) #

uniaxialMaterial('Concrete02', 2, -fc_unconfined, -epc_unconfined,  -fcu_unconfined,  -ecu_unconfined, lambda_unconfined, ft_unconfined, Et_unconfined) #

uniaxialMaterial('ReinforcingSteel', 4, fy, fu, Es, Esh, eps_sh, eps_ult,'-IsoHard', 4.3, 0.4, '-GABuck', lsr, beta_buckling, r_buckling,gamma_buckling, '-CMFatigue', Cf_fatigure, alpha_fatigue, Cd_fatigue)  

uniaxialMaterial("Steel02",    5,    fyt,    Es,  0.004    ,15.,    0.925,    0.15)

#======================================================================
# ---------------
# Node definition
# ---------------    
node(    1,  *[ 0.  ,    0.0 ]) # Support
node(    2,  *[ 0.  ,    0.0 ]) # Left of the Plastic Hinge Model
node(    3,  *[ 0.  ,    LpCol]) # Right of the Plastic Hinge Model
node(    4,  *[ 0.  ,    L   ]) # Mass concentration point 
#======================================================================
# ------------------
# Element definition
# ------------------  
    #            tag
geomTransf('Corotational',    1)
geomTransf('PDelta',    2)
    #                            Tag    NodeI NodeJ         A            E            Iz            transfTag   
# Elastic Beam 
element('ElasticTimoshenkoBeam', 1000, *[3, 4], Ec_Column, G, AgCol,   I_effCol, AshearCol, 2   ) # Elastic Element  
#====================================================================  

Num_Concrete_Fibers=len(As_list)#2*(Num_upperhalf_nodes-1)           # Number of the concrete fibers
conf_thicknessCol=(hCol-2*dpCol)/(Num_Concrete_Fibers-1)
width_list=[0.0 for i in range(Num_Concrete_Fibers)]
matConcreteTags=[0 for i in range(Num_Concrete_Fibers)]
matSteelTags=[Steel_mat_index_col for i in range(Num_Concrete_Fibers)]

for i in range(Num_Concrete_Fibers):
    if i==0 or i==Num_Concrete_Fibers-1:
          width_list[i]=dpCol+0.5*conf_thicknessCol
          matConcreteTags[i]=Unconf_Concrete_mat_index
    else:
          width_list[i]=conf_thicknessCol
          matConcreteTags[i]=Conf_Concrete_mat_index

thickness_list=[bCol for i in range(Num_Concrete_Fibers)]
    #---------------------------------------------------------------
rhos_=list(np.array(As_list)/(np.array(width_list)*np.array(thickness_list)))

element('MVLEM',1, 0.0, *[2,3] , Num_Concrete_Fibers, 0.5,\
        '-thick', *thickness_list,\
        '-width', *width_list,\
        '-rho', *rhos_, \
        '-matConcrete', *matConcreteTags, \
        '-matSteel', *matSteelTags, \
        '-matShear', ShearSteelMat) 
#====================================================================

fix(2,   *[1,1,1])

fix(1, 1, 1, 1)
print('Constraints have been built!') 
#======================================================================       
#------------
# Apply Mass
#------------

mass(4, Mass,  Mass, 1e-12)
mass(3, 1e-12, 1e-12, 1e-12)
#mass(3,0.1*Mass,0.1*Mass,1e-12)

recorder('Node', '-file', f'{mypath}/DFree.out','-time', '-node', 4, '-dof', 1,2,3, 'disp')
recorder('Node', '-file', f'{mypath}/RBase.out','-time', '-node', 2, '-dof', 1,2,3, 'reaction')
recorder('Node', '-file', f'{mypath}/AccFree.out','-time', '-node', 4, '-dof', 1,2,3, 'accel')

print('====================================================================')

# ---------------------------------
# Apply axial load and set constant
# ---------------------------------

    #          type     tag
timeSeries("Linear",  1)
    #          type  node    time series tag
pattern("Plain",    1,                1)
    #          node  load
load(4, *[0.0 , Axial_Load, 0.0])

Tol = 1e-12 # convergence tolerance for test

NstepGravity = 10

DGravity = 1/NstepGravity

integrator('LoadControl', DGravity) # determine the next time step for an analysis

numberer('Plain') # renumber dof's to minimize band-width (optimization), if you want to

system('BandGeneral')#'BandGeneral') # how to store and solve the system of equations in the analysis

constraints('Plain') # how it handles boundary conditions

test('NormDispIncr', Tol, 6) # determine if convergence has been achieved at the end of an iteration step

algorithm('Newton') # use Newton's solution algorithm: updates tangent stiffness at every iteration

analysis('Static') # define type of analysis static or transient

analyze(NstepGravity) # apply gravity

loadConst('-time', 0.0) #maintain constant gravity loads and reset time to zero

# calculate eigenvalues & print results

numEigen = 2

eigenValues = eigen('-fullGenLapack',numEigen)

PI = 2 * asin(1.0)

Period = np.zeros((numEigen,1))

Omegas=np.sqrt(eigenValues)

# ---------------------------------
# Apply damping
# ---------------------------------
T1=2*pi/Omegas[0]           # 1st Period

xDamp = 0.02                # 2% damping ratio

betaKcomm = 2 * xDamp/(Omegas[0]+Omegas[1])

alphaM =  2 * (xDamp*Omegas[0]*Omegas[1])/(Omegas[0]+Omegas[1])             # M-prop. damping; D = alphaM*M 

betaKcurr =  0.0        # K-proportional damping;      +beatKcurr*KCurrent

betaKinit = 0.0 # initial-stiffness proportional damping      +beatKinit*Kini

rayleigh(alphaM,betaKcurr, betaKinit, betaKcomm) # RAYLEIGH damping

AmplificationFactor=2.0

GMfact = g* AmplificationFactor         # data in input file is in g Units -- ACCELERATION TH

maxNumIter = 10

IDloadTag = 400         # load tag

dt = 0.01           # time step for input ground motion 

GMfile='Nishida_Kobe-NS-000.AT2'

GMdirection =1

EQ=np.loadtxt(GMfile)

timeSeries('Path', 2, '-dt', dt, '-filePath', GMfile, '-factor', GMfact)

pattern('UniformExcitation', IDloadTag, GMdirection, '-accel', 2)

wipeAnalysis()

Tol = 1e-12 # 1e-8 convergence tolerance for test

constraints('Transformation')

numberer('Plain')

system('FullGeneral')

test('EnergyIncr', Tol, maxNumIter)

algorithm('ModifiedNewton')

# ---------------------------------
# Main Analysis Part
# ---------------------------------
NewmarkGamma = 0.5

NewmarkBeta = 0.25

integrator('Newmark', NewmarkGamma, NewmarkBeta)

#integrator('CentralDifference')

analysis('Transient')

constraints('Transformation')

numberer('Plain')

system('FullGeneral')

test('EnergyIncr', Tol, maxNumIter)

algorithm('ModifiedNewton')

# set some variables
DtAnalysis=0.01
TmaxAnalysis = len(EQ)*dt
tCurrent = getTime()
ok = 0

time = [tCurrent]
u3 = [0.0]

u3 = [0.0]

tests = {1:'NormDispIncr'}#{1:'NormDispIncr', 2: 'RelativeEnergyIncr', 4: 'RelativeNormUnbalance',5: 'RelativeNormDispIncr', 6: 'NormUnbalance'}
algorithms = {1:'KrylovNewton'}#, 2: 'SecantNewton' , 4: 'RaphsonNewton',5: 'PeriodicNewton', 6: 'BFGS', 7: 'Broyden', 8: 'NewtonLineSearch'}

ok=0
while tCurrent<TmaxAnalysis and ok==0:
    algorithm('ModifiedNewton')
    test('EnergyIncr', Tol, maxNumIter)
    ok = analyze(1, DtAnalysis) 
    tCurrent = getTime()
    if ok != 0:

        for jj in range(0,4):
            if ok==0:
                    break

            for kk in range(0,4):
                if ok==0:
                    break

                for i in tests:
                    if ok==0:
                        break

                    for j in algorithms:
                        if ok==0:
                            break

                        if j < 4:
                            algorithm(algorithms[j], '-initial')

                        else:
                            algorithm(algorithms[j])

                        test(tests[i], Tol*10**kk, 1000)
                        ok = analyze(10**(jj), DtAnalysis*10**(-jj))                            

                        if ok == 0:
                            tCurrent = getTime()
                            #print('============================================')
                            print(f'Analysis succeeded at t={getTime():.4f} with {algorithms[j]} {tests[i]} and Tol={Tol*10**kk:.2E}')
                            #print('============================================')

                            break

                        else:
                            #print(f'Analysis failed at time {getTime():.4f} with {tests[i]}, {algorithms[j]}')
                            continue

    else:
        print(f'Analysis succeeded at time {getTime():.4f}')

u2 = nodeDisp(4, 1)
print("u2 = ", u2)

wipe()

Base_reaction=np.loadtxt(f'{mypath}/RBase.out', delimiter=" ")
free_node_disp=np.loadtxt(f'{mypath}/DFree.out', delimiter=" ")
free_node_acc=np.loadtxt(f'{mypath}/AccFree.out', delimiter=" ")

Time=free_node_disp[:,0]
t0=np.where(Time>=1.0)[0][0]+1
Time=Time[t0::]
AppliedDisplacement=free_node_disp[t0::,1]#[readDispdata[i][1] for i in range(len(readDispdata))]

plt.figure(figsize=(8, 6))
plt.plot(Time,AppliedDisplacement,c='b',label=f'MVLEM')
plt.ylim(-200,200)
plt.legend(loc='upper right', borderaxespad=0.)
plt.ylabel('Relative Displacement (mm)')
plt.xlabel('Time (sec)')
plt.xlim(0,TmaxAnalysis)
plt.grid(which='major', axis='both')
plt.show()

plt.figure(figsize=(8, 6))
plt.plot(Time, Base_reaction[t0::,2]/(Weight),c='b',label='Base Reaction')
plt.plot(Time, (1+free_node_acc[t0::,2]/g),c='r',label='External Forces inclduing Inertia')   
plt.legend(loc='upper right', borderaxespad=0.)
plt.ylabel('Axial Force/Applied Force')
plt.xlabel('Time (sec)')
plt.xlim(0,TmaxAnalysis)
plt.grid(which='major', axis='both')
plt.show()  

plt.figure(figsize=(8, 6))
plt.plot(Time, Base_reaction[t0::,2]/(Weight+Weight*free_node_acc[t0::,2]/g),c='b',label='Base Reaction')    
plt.legend(loc='upper right', borderaxespad=0.)
plt.ylabel('Axial Force/Applied Force')
plt.xlabel('Time (sec)')
plt.xlim(0,TmaxAnalysis)
plt.grid(which='major', axis='both')
plt.show()  

[mhscott]

@kkolozvari would you be able to take a look at this?

[kkolozvari]

I am not sure if the problem is with the element or with OpenSees settings. I dont use OpenSeesPy, so I cant say if there is something there, but I never had this type of problem with the element. Maybe try using disp beam column (or similar) and see if get the same problem?

[zwaezi]

Thanks Dr. @kkolozvari and Dr. @mhscott . Interestingly, for this analysis when I don't consider the Rayleigh damping at all, the following results will be found meaning that the problem is resolved! Can it be caused by the "spurious damping forces" applied to zero-mass DOfs mentioned in https://onlinelibrary.wiley.com/doi/10.1002/eqe.2622 ?

[kkolozvari]

I don't think this is the issue. I used this element with Rayleigh damping many times and didn't have issues.

Do you really need corrotational transformation for mvlem when you're using a leaning PDelta column? You said earlier you are "assuming a zero vertical mass". If so, why are you applying vertical mass at node 4? Why are you applying mass to the P-Delta column instead of wall element? I see a few "funny" things in the model just by a quick glance. I suggest you do a detailed checking of the model.

Good luck.

[fmckenna]

any updates?

[kkolozvari]

I was not able to replicate this issue. Not sure if this is a problem with this particular model, or something else...

[mhscott]

@zwaezi Any updates?