I am using a simple model with 1 force-based beam-column element with a fiber section to represent the configuration of the wall. Concrete02 and Steel02 material models are used to represent unconfined concrete and reinforcing steel respectively. 5 integration points are used along the height of the wall. Material regularization is accounted for in the material model parameters.
With these modeling assumptions, I am unable to capture the strength degradation at higher drifts. After initial rapid strength degradation post-peak, the wall model shows hardening behavior which is not expected. I am attaching my OpenSeePy code here, if you have any insights on what might be wrong with the model it'll be highly appreciated.
Code: Select all
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 2 13:04:47 2022
@author: pkakoty
"""
# -*- coding: utf-8 -*-
"""
Created on Sun Jan 30 16:49:42 2022
@author: pkakoty
"""
import openseespy.opensees as ops
import numpy as np
ops.wipe()
#########################
# Define Units
#########################
mm = 1.
m = 1000 * mm
kg = 1.
s = 1.
N = 1.
kN = N * 1000
Pa = N / m ** 2
kPa = Pa * 10 ** 3
MPa = Pa * 10 ** 6
GPa = Pa * 10 ** 9
ksi = MPa * 6.89476
psi = ksi * 10**(-3)
g = 9.8 * m * s**-2
########################
# Geometric definitions
#########################
story_height = 3.84 * m # Height of wall
no_story = 1
dbe_ele = 1 # six diplacement-based beam-column element per story
wall_length = 1920 * mm # Total length of the wall
wall_thickness = 200 * mm # Thickness of the wall
#########################
# Model Definitions
#########################
ops.model('basic','-ndm', 2)
#########################
# Define Nodes
#########################
print("Defining Nodes")
# Nodes for the wall system
node_data = np.zeros((no_story*dbe_ele + 1, 3))
node_height = 0
node_ind = 0
for i in range(no_story + 1):
if i < no_story:
for j in range(dbe_ele):
node_data[node_ind, 0] = 1000*(i+1) + (j+1)
node_data[node_ind, 1] = 0
node_data[node_ind, 2] = node_height
node_ind = node_ind + 1
node_height = node_height + story_height/(dbe_ele)
else:
node_data[node_ind, 0] = 1000*(i+1) + 1
node_data[node_ind, 1] = 0
node_data[node_ind, 2] = node_height
for i in range(node_data.shape[0]):
ops.node(int(node_data[i,0]), node_data[i,1], node_data[i,2])
ops.node(100, 0, 0)
ops.fix(int(node_data[0,0]),*[1, 1, 1]) # Wall node at the base
ops.fix(100,*[0, 1, 1]) # Shear spring node at the base
###########################################################################
############################# Define Materials ############################
###########################################################################
print("Defining Materials")
# Unconfined Concrete Properties
fpc = - 41.73 * MPa # concrete compressive strength at 28 days
Ecc = 23.3 * GPa
epsc0 = -0.00225 # concrete strain at maximum strength
fpcu = 0.2*fpc # Concrete crushing strength
epsU = -0.004289 # Concrete strain at crushing strength
lambda_ = 0.1 # ratio between unloading slope at $epscu and initial slope
ft = 2.73 * MPa # Tensile strength
Ets = 0.05*Ecc # Tenson softening stiffness
# ops.uniaxialMaterial('Concrete02', matTag, fpc, epsc0, fpcu, epsU, lambda, ft, Ets)
ops.uniaxialMaterial('Concrete02', 1, fpc, epsc0, fpcu, epsU, lambda_, ft, Ets)
# Steel Section
Fy = 296.7 * MPa
E0 = 200. * GPa
b = 0.0027
params = [12.95, 0.925, 0.15]
# uniaxialMaterial('Steel02', matTag, Fy, E0, b, *params, a1=a2*Fy/E0, a2=1.0, a3=a4*Fy/E0, a4=1.0, sigInit=0.0)
ops.uniaxialMaterial('Steel02', 3, Fy, E0, b, *params)
# Min Max Wrapper
maxStrain = 0.13
minStrain = epsU
# uniaxialMaterial('MinMax', matTag, otherTag, '-min', minStrain=1e-16, '-max', maxStrain=1e16)
ops.uniaxialMaterial('MinMax', 30, 3, '-max', maxStrain, '-min', minStrain)
# Shear Spring Section
G = 0.4 * Ecc # As per ACI 318 [Marafi et al. (2019)]
kG = G * (5/6) * wall_thickness # Marafi et al. (2019)
# uniaxialMaterial('Elastic', matTag, E, eta=0.0, Eneg=E)
ops.uniaxialMaterial('Elastic', 4, kG)
# Leaning column spring material
# uniaxialMaterial('Elastic', matTag, E, eta=0.0, Eneg=E)
ops.uniaxialMaterial('Elastic', 5, kG/1000000.)
# Leaning column and truss material
# uniaxialMaterial('Elastic', matTag, E, eta=0.0, Eneg=E)
ops.uniaxialMaterial('Elastic', 6, kG*1000000.)
###########################################################################
########################### Define Fiber Sections #########################
###########################################################################
# Section pre-processing
# Single layer reinforcement at the mid-layer of the wall
y_rebar_mid = np.array([-684, -456, -228, 0, 228, 456, 684]) * mm
y_rebar_end = np.array([-912, 912]) * mm
z_rebar = np.array([0]) * mm
Abar_mid = np.pi * (12 * mm/ 2.)**2
Abar_end = np.pi * (20 * mm/ 2.)**2
Nbar_mid = len(y_rebar_mid) * len(z_rebar)
Nbar_end = len(y_rebar_end) * len(z_rebar)
rebar_coord_mid = np.zeros((Nbar_mid,2))
rebar_coord_end = np.zeros((Nbar_end,2))
for ii, y_cord in enumerate(y_rebar_mid):
for jj, z_cord in enumerate(z_rebar):
rebar_coord_mid[ii*len(z_rebar)+jj, :] = [y_cord, z_cord]
for ii, y_cord in enumerate(y_rebar_end):
for jj, z_cord in enumerate(z_rebar):
rebar_coord_end[ii*len(z_rebar)+jj, :] = [y_cord, z_cord]
# Unconfined concrete in bottom boundary zone for RECT patch
unconfined_vertices = [-(wall_length/2), -(wall_thickness/2),
(wall_length/2), (wall_thickness/2)]
#########################
# Define Fiber Section
#########################
print("Defining Fiber Sections")
# section('Fiber', secTag, '-GJ', GJ)
ops.section('Fiber', 1)
# Define fiber sections for all the rebars
# fiber(yloc, zloc, A, matTag)
for ii in range(rebar_coord_mid.shape[0]):
ops.fiber(rebar_coord_mid[ii,0], rebar_coord_mid[ii,1], Abar_mid, 3)
for ii in range(rebar_coord_end.shape[0]):
ops.fiber(rebar_coord_end[ii,0], rebar_coord_end[ii,1], Abar_end, 3)
# Rectangular wall section: unconfined concrete
# patch('rect', matTag, numSubdivY, numSubdivZ, *crdsI, *crdsJ)
ops.patch('rect', 1, 10, 4, *unconfined_vertices)
# Define transform and integration
# geomTransf('Linear', transfTag, '-jntOffset', *dI, *dJ)
ops.geomTransf('Linear', 1)
# Define beam integration method and points
# beamIntegration('Lobatto', tag, secTag, N)
ops.beamIntegration('NewtonCotes', 1, 1, 5)
#########################
# Define Elements
#########################
# element('dispBeamColumn', eleTag, *eleNodes, transfTag, integrationTag, '-cMass', '-mass', mass=0.0)
# element('elasticBeamColumn', eleTag, *eleNodes, Area, E_mod, Iz, transfTag, <'-mass', mass>, <'-cMass'>, <'-release', releaseCode>)
# ele tags - 1-x
print("Defining Elements (Nonlinear Elements for cantilevered wall system)")
ele_counter = 0
for i in range(node_data.shape[0]-1):
# Elements in the cantilevered wall section
ele_nodes = [int(node_data[i,0]), int(node_data[i+1,0])]
# print(ele_nodes)
# element('dispBeamColumn', eleTag, *eleNodes, transfTag, integrationTag, '-cMass', '-mass', mass=0.0)
# ops.element('dispBeamColumn', i+1, *ele_nodes, 1, 1)
# # element('forceBeamColumn', eleTag, *eleNodes, transfTag, integrationTag, '-iter', maxIter=10, tol=1e-12, '-mass', mass=0.0)
ops.element('forceBeamColumn', i+1, *ele_nodes, 1, 1)
ele_counter += 1
ops.element('zeroLength', 1001, *[100,1001], '-mat', *[4, 5, 5], '-dir', *[1, 2, 3])
#########################
# Load definitions
#########################
axial_load = 0.035 * fpc * wall_length * wall_thickness
# timeSeries('Constant', tag, '-factor', factor=1.0)
ops.timeSeries("Constant", 1)
# pattern('Plain', patternTag, tsTag, '-fact', fact)
ops.pattern('Plain', 1, 1)
# load(nodeTag, *loadValues)
ops.load(2001, *[1e-10, axial_load, 1e-10])
# mass = 2400 * wall_length * wall_thickness * story_height
# ops.mass(int(2001), *[mass, 1e-10, 1e-10])
###########################################################################
########################### Define Gravity Analysis #######################
###########################################################################
ops.wipeAnalysis()
ops.system('BandSPD')
ops.constraints('Transformation')
ops.numberer('RCM')
ops.test('NormDispIncr', 1.0e-12, 10)
ops.algorithm('Newton')
# ops.integrator('LoadControl', 0.1, 4, 0.001, 1)
ops.integrator('LoadControl', 1)
ops.analysis('Static')
ops.analyze(1)
ops.loadConst('-time', 0.0)
###########################################################################
########################### Define Pushover Analysis #######################
###########################################################################
sdr_file = 'pushover.out'
ops.recorder('Node', '-file', sdr_file,'-time', '-node', 2001, '-dof', 1,2,3, 'disp')
# Bottom left corner
conc_file = 'concrete_fiber1.out'
ops.recorder('Element', '-file', conc_file,'-ele', 1, 'section', 1, 'fiber', -960, -100, 1, 'stressStrain')
stl_file = 'steel_fiber1.out'
ops.recorder('Element', '-file', stl_file,'-ele', 1, 'section', 1, 'fiber', -960, -100, 3, 'stressStrain')
# Top left corner
conc_file = 'concrete_fiber2.out'
ops.recorder('Element', '-file', conc_file,'-ele', 1, 'section', 1, 'fiber', - 960, 100, 1, 'stressStrain')
stl_file = 'steel_fiber2.out'
ops.recorder('Element', '-file', stl_file,'-ele', 1, 'section', 1, 'fiber', -960, 100, 3, 'stressStrain')
# Top right corner
conc_file = 'concrete_fiber3.out'
ops.recorder('Element', '-file', conc_file,'-ele', 1, 'section', 1, 'fiber', 960, 100, 1, 'stressStrain')
stl_file = 'steel_fiber3.out'
ops.recorder('Element', '-file', stl_file,'-ele', 1, 'section', 1, 'fiber', 960, 100, 3, 'stressStrain')
# Bottom right corner
conc_file = 'concrete_fiber4.out'
ops.recorder('Element', '-file', conc_file,'-ele', 1, 'section', 1, 'fiber', 960, -100, 1, 'stressStrain')
stl_file = 'steel_fiber4.out'
ops.recorder('Element', '-file', stl_file,'-ele', 1, 'section', 1, 'fiber', 960, -100, 3, 'stressStrain')
# iDmax = [0.0015, 0.002, 0.0025, 0.0035, 0.005, 0.0075, 0.01, 0.015, 0.02, 0.03, 0.04]
# sw_test = pd.read_csv("experimental_results.csv")
# iDstep = np.array(sw_test["Disp"])
# iDmax = [0.0015, 0.002, 0.0025, 0.0035, 0.005, 0.0075, 0.01, 0.015]
Dincr = 0.5 #0.2 * story_height
# Fact = story_height
# CycleType = 'Full'
# Ncycles = 2
ControlNode = 2001
ControlNodeDof = 1
du = 0.01 * mm
maxU = 150 * mm
# currentDisp = 0.
# ok = 0
force = [0]
disp = [0]
# ops.timeSeries("Constant", 2)
# ##########################################################################
# ########################### Monotonic Pushover ###########################
# ##########################################################################
# pushover_series = 100
# pushover_pattern = 100
# pushover_direction = 1
# stop analysis at 10% drift
topLateralDisp = 0.05 * no_story * story_height
# set displacement increment as 0.001 for now, should adjust based on a expected ductility of the building
initial_disp_increment= 0.001 * no_story * story_height
ops.timeSeries('Linear', 3)
ops.pattern('Plain', 2, 3)
# set reference force as 1N at top floor
# ref_force = 1
# ops.sp(ControlNode, ControlNodeDof, 1)
ops.load(ControlNode, 1, 0., 0.)
# ops.wipeAnalysis()
ops.wipeAnalysis()
ops.constraints('Transformation')
ops.numberer('RCM')
ops.test('NormDispIncr', 1.0e-6, 100)
ops.system(('BandGen'))
ops.algorithm('KrylovNewton')
ops.integrator('DisplacementControl', ControlNode, ControlNodeDof, Dincr/3)
ops.analysis('Static')
analysis_step = 0
# topLateralDisp = 1
# current setup only adjusts step size and convergence criteria, add more analysis loop options if necessary
while ops.nodeDisp(ControlNode, ControlNodeDof) < topLateralDisp:
analysis_result = ops.analyze(1)
if analysis_result != 0:
print("***Smaller step of 0.0001 at step", analysis_step)
ops.integrator('DisplacementControl', ControlNode, ControlNodeDof, Dincr/ 10)
analysis_result = ops.analyze(1)
ops.integrator('DisplacementControl', ControlNode, ControlNodeDof, Dincr)
if analysis_result != 0:
print("***Smaller step of 0.00001 at step",analysis_step)
ops.integrator('DisplacementControl', ControlNode, ControlNodeDof, Dincr/ 100)
analysis_result = ops.analyze(1)
ops.integrator('DisplacementControl', ControlNode, ControlNodeDof, Dincr)
if analysis_result != 0:
print("***Smaller step of 0.000001 at step", analysis_step)
ops.integrator('DisplacementControl', ControlNode, ControlNodeDof, initial_disp_increment/ 1000)
ops.test('NormDispIncr',1.0e-6,1000)
analysis_result = ops.analyze(1)
ops.integrator('DisplacementControl', ControlNode, ControlNodeDof, initial_disp_increment)
ops.test('NormDispIncr',1.0e-6,100)
if analysis_result != 0:
print("***Analysis failed at step", analysis_step)
break
if analysis_result == 0:
analysis_step += 1
disp.append(ops.nodeDisp(ControlNode, ControlNodeDof))
force.append(-1 * ops.eleForce(dbe_ele, ControlNodeDof)/1000)
print(analysis_step)
list_size = np.size(disp)
hysteresis = np.zeros((list_size, 2))
hysteresis[:, 0] = disp
hysteresis[:, 1] = force
np.savetxt('hysteresis.csv', hysteresis, delimiter = ",")