Unable to capture strength degradation at large drift levels

Forum for asking and answering questions related to use of the OpenSeesPy module

Moderators: silvia, selimgunay, Moderators

Post Reply
pkakoty
Posts: 2
Joined: Mon Jul 18, 2022 4:07 pm

Unable to capture strength degradation at large drift levels

Post by pkakoty » Tue Aug 02, 2022 12:38 pm

I have been trying to calibrate material properties to develop models for non-ductile RC shear wall buildings using experimental data from a series of tests. These are non-ductile RC walls, with 200 mm thickness and a single layer of reinforcement. As we can expect, the walls exhibit rapid strength degradation and consequently fail at reasonable drifts (2-3%).

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 = ",")

chenwangNigel
Posts: 5
Joined: Wed Jun 20, 2018 6:31 am
Location: McGill University

Re: Unable to capture strength degradation at large drift levels

Post by chenwangNigel » Tue Sep 20, 2022 7:17 pm

Have you tried to use more elements and finer fiber mesh?

By the way, what do you mean by "material regularization"? Could you please give more information on it? Thank you.

Post Reply