MVLEM model - Issues with applying gravity load

[LiamPledger]

Hi All,

I've developed a simple wall structure using the MVLEM model in OpenSeesPy. I get a period of 0.378 sec before applying gravity load. (This matches well with the SAP2000 model and hand calcs done in MathCAD ~ 0.38 s). However, when I apply gravity load using the Load command and then get the period after gravity load analysis, the period is 0.93 sec. I know there is clearly something drastically wrong however I can't work out what exactly. The mode shape does not change much at all even though period increases. I have checked units and that is not the problem. I have developed an RC frame model in opensees previously and period changes approx 1% when applying gravity load.

Even changing the direction of the load to +ve (upwards) which should reduce the period by increasing stiffness leads to the same period of 0.93 sec. I have also tried changing the transformation from Corotational to Linear and this doesn't do anything.

Any help would be greatly appreciated,
Cheers,
Liam


Code for reference :

import openseespy.opensees as ops
import numpy as np
import openseespy.postprocessing.Get_Rendering as opsplt

ops.wipe()
ops.model('basic', '-ndm', 2, '-ndf', 3)
ops.geomTransf('Corotational', 1) # Uses the natural deformation scheme, accounting for PDelta and rigid body motion - typically the most accurate

# all units in kN, m, and sec

# Wall dimensions
t = 0.3 # m
lw = 6.0 # m

# concrete elastic modulus
Ec = 25* 10**6 # kN / m^2

# define structure-geometry parameters
NStories = 5 # number of stories
NBays = 5 # number of frame bays (excludes bay for P-delta column)
WBay = 6.0 # bay width in meters

# Define element nodes
Wall_nodes = [[101, 102], [102, 103], [103, 104], [104, 105], [105, 106]]

# Define nodes
ntag = [101, 102, 103, 104, 105, 106]
cor = np.array([[0, 0, 0, 0, 0, 0], [0, 3.6, 7.2, 10.8, 14.4, 18.0]])

# Nodes
for i in range(0, len(ntag)):
ops.node(int(ntag), cor[0, i], cor[1, i])

# Define material properties for ConcreteCM confined
matTag_concrete = 1 # confined concrete
fpcc = 53.37*10**3 # Peak compressive strength kN/m2
epcc = 0.0106 # Strain @ peak compressive strength
# fpcc and epcc based on Mander J.B., Priestley M.J.N., and Park R. (1988). “Theoretical Stress-Strain Model for Confined Concrete”, ASCE Journal of Structural Engineering, V. 114, No. 8, pp. 1804-1826.

rc = fpcc/1000/6.68 - 1.85 # Ratio of unloading slope to initial slope
# rc based on : Tsai W.T. (1988), “Uniaxial Compressional Stress-Strain Relation of Concrete”, ASCE Journal of Structural Engineering, V. 114, No. 9, pp. 2133-2136.

xcrn = 0.002 # Compressive strain at which reloading begins
ft = 2 # Tensile strength
et = 0.0001 # Ultimate tensile strain
rt = 1.2 # Tension softening stiffness ratio
xcrp = 10000 # Tensile strain at which tension stiffening begins
mon = 0 # Flag to turn monotonic behavior on (1) or off (0)
gap_close = 0 # Gap closing parameter


# Define material properties
ops.uniaxialMaterial('ConcreteCM', matTag_concrete, fpcc, epcc, Ec, rc, xcrn, ft, et, rt, xcrp, '-GapClose', gap_close)


# Define material properties for SteelMPF
matTag_steel = 2
fyp = 540*10**3 # Yield strength in tension
fyn = -540*10**3 # Yield strength in compression
E0_steel = 200*10**6 # Initial elastic modulus of steel
bp = 0.02 # Strain-hardening ratio in tension
bn = -0.02 # Strain-hardening ratio in compression
a1 = 0.0
a2 = 1.0
a3 = 0.0
a4 = 1.0

R0=20; cR1=0.925; cR2=0.15
params=[R0,cR1,cR2]

ops.uniaxialMaterial('SteelMPF', matTag_steel, fyp, fyn, E0_steel, bp, bn, *params, a1, a2, a3, a4)

# SHEAR ........................................................
# uniaxialMaterial Elastic $matTag $E <$eta> <$Eneg>
# NOTE: large shear stiffness assigned since only flexural response
Ag = t * lw # Gross area of the wall cross section
mu = 0.1
G=Ec/(2*(1+mu)) # Shear Modulus
shear_matTag = 3

# Build shear material
ops.uniaxialMaterial('Elastic', shear_matTag, G*Ag)

# Define section properties
m = 5
thick = [t] * m # Thickness of each element
concrete_list = [matTag_concrete] * m
steel_list = [matTag_steel] * m
widths = [0.9, 1.4, 1.4, 1.4, 0.9] # Width of each element in mm
rho = [0.03, 0.01, 0.01, 0.01, 0.03] # Reinforcement ratio of each element

density = 2350*9.81/1000 # density of the wall in units kN/m^3
density = 0
c = 0.4 # location of the center of rotation


# Define MVLEM elements
for i in range(NStories):
eleTag = i + 1
ops.element('MVLEM',int(eleTag), float(density), *Wall_nodes, int(m), float(c), '-thick', *thick, '-width', *widths, '-rho', *rho, '-matConcrete', *concrete_list,'-matSteel', *steel_list, '-matShear', shear_matTag)

ops.fix(101, 1, 1, 1) # fixed


" Define mass distribution of nodes "
mass_frame = 150 # in Tonnes

for i in range(0, len(ntag)):
"mass command is used to set the mass at each node"
ops.mass(int(ntag), mass_frame, mass_frame, 0)


"Applying vertical load to frame to determine the building period"

"Create the nodal load - command: load nodeID xForce yForce"

# Gravity load
opsplt.createODB("Nonlin_RCWall", "Gravity", Nmodes=3)

opsplt.plot_modeshape(1, 1, Model = "Nonlin_RCWall")


ops.timeSeries('Linear', 1) # applies the load in a linear manner (not all at once)
ops.pattern("Plain", 1, 1) # create a plain load pattern - similar to ELF

for i in range(0, len(ntag)):
"mass command is used to set the mass at each node"
ops.load(int(ntag), 0, -mass_frame*9.81, 0)


"(mode shape number, scale, ModelName - Displays the model saved in a database named"

# create DOF number
ops.numberer("RCM")
# create SOE
ops.system('BandGeneral')
# create constraint handler
ops.constraints("Transformation")
# create number of steps
nsteps=25
# create integrator
ops.integrator('LoadControl', 1/nsteps)
# create algorithm
ops.test('NormUnbalance', 1e-5, 200, 0) # convergence scheme applied to the model
ops.algorithm("Newton") # solution scheme - Newton-Raphson method of solving non-linear equations
# create analysis object
ops.analysis("Static") # analysis type - ie, static, transient etc..
# perform the analysis
ops.recorder('Node', '-file', "results/modal/eigen.out",'-closeOnWrite','-dof',1,2,3,'eigen')
ops.analyze(nsteps)


opsplt.createODB("Nonlin_RCWall", "Gravity", Nmodes=3)
opsplt.plot_modeshape(1, 1)
ops.eigen('-genBandArpack', 1)
ops.record
ops.loadConst('-time', 0.0)
ops.wipeAnalysis()