new IMKBilin stays elastic

[silviamazzoni]

@amaelkady can you, please, help me out with the new IMKBilin? I am trying to compare it to the old Bilin and HystereticSM, but I can't get it to yield. I am using a script from Laura Eads that I found in the wiki. I am pretty sure I should have the input right. so the script is long, but i minimized it to what is relevant. also, your documentation has "set c_A 1.00;" in the example, but I don't see it in the documentation, so I am concerned I am missing an input. here is the script:

[silviamazzoni]

`#!/usr/bin/env python
# coding: utf-8

# ## BUT I AM JUST TESTING MATERIALS HERE
# 
# # Concentrated-Plasticity Modeling
# # -- Concentrated-Plasticity Integration --
# ## Pushover Analysis of 2-Story Moment Frame (without panel zones)
# ### Silvia Mazzoni, May 2022
# #### silviamazzoni@yahoo.com
# 
# This model uses the a DIFFERENT methodology from the original example.
# Here we use the force beam-column element with the concentrated-plasticity integration scheme
# 

# In[1]:

# This example was taken from Laura Eads example in the OpenSees Wiki page:
# Pushover Analysis of 2-Story Moment Frame (without panel zones)
# https://opensees.berkeley.edu/wiki/index.php/Pushover_Analysis_of_2-Story_Moment_Frame
# I have converted the Tcl code to Python and have added the Concentrated-Plasticity Model I developed

## BUT I AM JUST TESTING MATERIALS HERE

# In[2]:

# import sys 
# OpenSeesPyPath = r'D:\Projects\OpenSees\Development\OpenSeesFork\Win64\bin' 
# sys.path.append(OpenSeesPyPath) 
# import opensees as ops 

import openseespy.opensees as ops
import matplotlib.pyplot as plt
import numpy as np
plt.rc('font',size=6)

# In[3]:

## modified in 2022 to Python by Silvia Mazzoni
# https://opensees.berkeley.edu/wiki/index.php/Pushover_and_Dynamic_Analyses_of_2-Story_Moment_Frame_with_Panel_Zones_and_RBS

# In[4]:

def defineStrainHistory(peaksArray,scaleFactor,nSteps,nCycles):
    strain = []
    for thisPeak in peaksArray:
        for i in range(nCycles):
            strain = np.append(strain,np.linspace(0,thisPeak*scaleFactor,nSteps))
            strain = np.append(strain,np.linspace(thisPeak*scaleFactor,-thisPeak*scaleFactor,nSteps))
            strain = np.append(strain,np.linspace(-thisPeak*scaleFactor,0,nSteps))

    return strain

peaksArray=np.array([1,2,3,4,5,6,7,8,9,10])
# peaksArray=[1,10]
scaleFactor = 0.1
# scaleFactor = 0.01
nSteps = 100
nCycles = 1
strain=defineStrainHistory(peaksArray,scaleFactor,nSteps,nCycles)
plt.plot(strain)

def testMaterial(materialTag,strainAmp,title):

    ops.testUniaxialMaterial(materialTag)
    outstrain = []
    outstress = []
    MUy = []
    for eps in strain:
        ops.setStrain(strainAmp*eps)
        outstrain.append(ops.getStrain())
        outstress.append(ops.getStress())
        tangent = ops.getTangent() # Not used

    plt.plot(outstrain,outstress)
    plt.title(title)
    plt.grid()
    plt.show()
    return outstrain,outstress

# In[5]:

############################################################################################
# rotSect2DModIKModel.tcl

#
# This routine creates a section with elastic axial and bilinear flexural response
# 
# Behavior follows: Bilinear Response based on Modified Ibarra Krawinkler Deterioration Model 
#
# Written by: Dimitrios G. Lignos, Ph.D.
#
# Variables
#   $secID =    Element Identification (integer) 
#   $E =        Young's modulus
#   $A =        Cross-sectional area
#   $K =        Initial stiffness after the modification for n (see Ibarra and Krawinkler, 2005)
#   $asPos =    Strain hardening ratio after n modification (see Ibarra and Krawinkler, 2005)
#   $asNeg =    Strain hardening ratio after n modification (see Ibarra and Krawinkler, 2005)
#   $MyPos =    Positive yield moment (with sign)
#   $MyNeg =    Negative yield moment (with sign)
#   $LS =       Basic strength deterioration parameter (see Lignos and Krawinkler, 2009)
#   $LK =       Unloading stiffness deterioration parameter (see Lignos and Krawinkler, 2009)
#   $LA =       Accelerated reloading stiffness deterioration parameter (see Lignos and Krawinkler, 2009)
#   $LD =       Post-capping strength deterioration parameter (see Lignos and Krawinkler, 2009)
#   $cS =       Exponent for basic strength deterioration
#   $cK =       Exponent for unloading stiffness deterioration
#   $cA =       Exponent for accelerated reloading stiffness deterioration
#   $cD =       Exponent for post-capping strength deterioration
#   $th_pP =    Plastic rotation capacity for positive loading direction
#   $th_pN =    Plastic rotation capacity for negative loading direction
#   $th_pcP =   Post-capping rotation capacity for positive loading direction
#   $th_pcN =   Post-capping rotation capacity for negative loading direction
#   $ResP =     Residual strength ratio for positive loading direction
#   $ResN =     Residual strength ratio for negative loading direction
#   $th_uP =    Ultimate rotation capacity for positive loading direction
#   $th_uN =    Ultimate rotation capacity for negative loading direction
#   $DP =       Rate of cyclic deterioration for positive loading direction
#   $DN =       Rate of cyclic deterioration for negative loading direction
#
# References:
#       Ibarra, L. F., and Krawinkler, H. (2005). “Global collapse of frame structures under seismic excitations,” Technical Report 152, The John A. Blume Earthquake Engineering Research Center, Department of Civil Engineering, Stanford University, Stanford, CA.
#       Ibarra, L. F., Medina, R. A., and Krawinkler, H. (2005). “Hysteretic models that incorporate strength and stiffness deterioration,” International Journal for Earthquake Engineering and Structural Dynamics, Vol. 34, No.12, pp. 1489-1511.
#       Lignos, D. G., and Krawinkler, H. (2010). “Deterioration Modeling of Steel Beams and Columns in Support to Collapse Prediction of Steel Moment Frames”, ASCE, Journal of Structural Engineering (under review).
#       Lignos, D. G., and Krawinkler, H. (2009). “Sidesway Collapse of Deteriorating Structural Systems under Seismic Excitations,” Technical Report 172, The John A. Blume Earthquake Engineering Research Center, Department of Civil Engineering, Stanford University, Stanford, CA.
#
############################################################################################
## modified in 2022 to Python by Silvia Mazzoni

def rotSect2DModIKModel(secID,E,A,K,asPos,asNeg,MyPos,MyNeg,LS,LK,LA,LD,cS,cK,cA,cD,th_pP,th_pN,th_pcP,th_pcN,ResP,ResN,th_uP,th_uN,DP,DN):
    print('---------------------------------------------------------')
    # Create the material and section
    ops.uniaxialMaterial('Bilin',secID,K,asPos,asNeg,MyPos,MyNeg,LS,LK,LA,LD,cS,cK,cA,cD,th_pP,th_pN,th_pcP,th_pcN,ResP,ResN,th_uP,th_uN,DP,DN);
#     ops.uniaxialMaterial('Elastic',10*secID+1,E*A);       # this is not used as a material, this is an axial-force-strain response
    ops.section('Aggregator',secID,secID,'Mz');     # only bending is plastic

    testMaterial(secID,0.1,'Bilin') # (MyPos/K+th_pP+th_pcP)/500

    # convert to hysteretic:
    M1 = MyPos
    r1 = M1/K
    M2 = M1 + th_pP*asPos*K
    r2 = r1 + th_pP
    k3 = (0-M2)/th_pcP
    M3 = ResP*M1
    r3 = r2+th_pcP*(M2-M3)/(M2)
    M4 = M3
    r4 = th_uP
    M5 = 0.01*M4
    r5 = 1.01*r4
    posEnv = [r1,M1,r2,M2,r3,M3,r4,M4,r5,M5]
    pinch = [.01,.99]
    ops.uniaxialMaterial('HystereticSM',secID+100,'-posEnv',*posEnv,'-pinch',*pinch,'-damage', 0, 0,'-degEnv',0,'-XYorder')
    testMaterial(secID+100,0.1,'HystereticSM')

    # convert to new IMKBilin
    dp_pos = th_pP
    dpc_pos = th_pcP
    du_pos = th_uP
    Fy_pos = MyPos
    FmaxFy_pos = MyPos + th_pP*asPos*K
    FresFy_pos = ResP*MyPos
    EnvPos = [dp_pos,dpc_pos,du_pos,Fy_pos,FmaxFy_pos,FresFy_pos]
    print('EnvPos',EnvPos)

    dp_neg = th_pN
    dpc_neg = th_pcN
    du_neg = th_uN
    Fy_neg = MyNeg
    FmaxFy_neg = MyNeg + th_pN*asNeg*K
    FresFy_neg = ResN*MyNeg
    EnvNeg = [dp_neg,dpc_neg,du_neg,Fy_neg,FmaxFy_neg,FresFy_neg]

    Lamdas = [.5,.5,.5]
    Cs = [1,1,1]
    Ds = [1,1]
    ops.uniaxialMaterial('IMKBilin',secID+10000,K,*EnvPos,*EnvNeg,*Lamdas,*Cs,*Ds)

    testMaterial(secID+10000,0.1,'IMKBilin')

# In[ ]:

# In[6]:

# #### LATERAL ANALYSIS

# # --------------------------------------------------------------------------------------------------
# # Example: 2-Story Steel Moment Frame with Concentrated Plasticity
# # Centerline Model with Concentrated Plastic Hinges at Beam-Column Joint
# # Created by:  Laura Eads, Stanford University, 2010
# # Units: kips, inches, seconds

###################################################################################################
#          Up & Source Definition                                     
###################################################################################################
ops.wipe();                         # clear memory of past model definitions
ops.model('BasicBuilder','-ndm',2,'-ndf',3);    # Define the model builder, ndm = #dimension, ndf = #dofs

# ###################################################################################################
# #          Define Building Geometry, Nodes, and Constraints                                             
# ###################################################################################################
# # define structure-geometry parameters
WBay=30.0*12.0;     # bay width in inches
HStory1=15.0*12.0;      # 1st story height in inches
HStoryTyp=12.0*12.0;        # story height of other stories in inches

# In[7]:

###################################################################################################
#          Define Section Properties and Elements                                                     
###################################################################################################
# define material properties
Es=29000.0;         # steel Young's modulus

# define column section W24x131 for Story 1 & 2
Acol_12=38.5;       # cross-sectional area
Icol_12=4020.0; # moment of inertia
Mycol_12=20350.0;   # yield moment

# define beam section W27x102 for Floor 2 & 3
Abeam_23=30.0;      # cross-sectional area (full section properties)
Ibeam_23=3620.0;    # moment of inertia  (full section properties)
Mybeam_23=10938.0;  # yield moment at plastic hinge location (i.e., My of RBS section, if used)
# note: In this example the hinges form right at the beam-column joint, so using an RBS doesn't make sense; 
#       however, it is done here simply for illustrative purposes.

# determine stiffness modifications to equate the stiffness of the spring-elastic element-spring subassembly to the stiffness of the actual frame member
# Reference:  Ibarra, L. F., and Krawinkler, H. (2005). "Global collapse of frame structures under seismic excitations," Technical Report 152,
#             The John A. Blume Earthquake Engineering Research Center, Department of Civil Engineering, Stanford University, Stanford, CA.
# calculate modified section properties to account for spring stiffness being in series with the elastic ops.element(stiffness
n=10.0;     # stiffness multiplier for rotational spring

# calculate modified moment of inertia for elastic elements
Icol_12mod=Icol_12*(n+1.0)/n;   # modified moment of inertia for columns in Story 1 & 2
Ibeam_23mod=Ibeam_23*(n+1.0)/n; # modified moment of inertia for beams in Floor 2 & 3
# calculate modified rotational stiffness for plastic hinge springs
Ks_col_1=   n*6.0*Es*Icol_12mod/HStory1;        # rotational stiffness of Story 1 column springs
Ks_col_2=   n*6.0*Es*Icol_12mod/HStoryTyp;  # rotational stiffness of Story 2 column springs
Ks_beam_23=n*6.0*Es*Ibeam_23mod/WBay;       # rotational stiffness of Floor 2 & 3 beam springs

# In[8]:

###################################################################################################
#          Define Rotational Springs for Plastic Hinges                                               
###################################################################################################
# define rotational spring properties and create spring elements using "rotSpring2DModIKModel" procedure
# rotSpring2DModIKModel(creates a uniaxial material spring with a bilinear response based on Modified Ibarra Krawinkler Deterioration Model
# references provided in rotSpring2DModIKModel.tcl
# input values for Story 1 column springs
McMy=1.05;          # ratio of capping moment to yield moment, Mc / My
LS=1000.0;          # basic strength deterioration (a very large # = no cyclic deterioration)
LK=1000.0;          # unloading stiffness deterioration (a very large # = no cyclic deterioration)
LA=1000.0;          # accelerated reloading stiffness deterioration (a very large # = no cyclic deterioration)
LD=1000.0;          # post-capping strength deterioration (a very large # = no deterioration)
cS=1.0;             # exponent for basic strength deterioration (c = 1.0 for no deterioration)
cK=1.0;             # exponent for unloading stiffness deterioration (c = 1.0 for no deterioration)
cA=1.0;             # exponent for accelerated reloading stiffness deterioration (c = 1.0 for no deterioration)
cD=1.0;             # exponent for post-capping strength deterioration (c = 1.0 for no deterioration)
th_pP=0.025;        # plastic rot capacity for pos loading
th_pN=0.025;        # plastic rot capacity for neg loading
th_pcP=0.3;         # post-capping rot capacity for pos loading
th_pcN=0.3;         # post-capping rot capacity for neg loading
ResP=0.4;           # residual strength ratio for pos loading
ResN=0.4;           # residual strength ratio for neg loading
th_uP=0.4;          # ultimate rot capacity for pos loading
th_uN=0.4;          # ultimate rot capacity for neg loading
DP=1.0;             # rate of cyclic deterioration for pos loading
DN=1.0;             # rate of cyclic deterioration for neg loading
a_mem=(n+1.0)*(Mycol_12*(McMy-1.0)) / (Ks_col_1*th_pP); # strain hardening ratio of spring
b=(a_mem)/(1.0+n*(1.0-a_mem));                          # modified strain hardening ratio of spring (Ibarra & Krawinkler 2005, note: Eqn B.5 is incorrect)

# define column springs

# Spring ID: "3xya", where 3 = col spring, x = Pier #, y = Story #, a = location in story
# "a" convention: 1 = bottom of story, 2 = top of story
# command: rotSpring2DModIKModel    id    ndR  ndC     K   asPos  asNeg  MyPos      MyNeg      LS    LK    LA    LD   cS   cK   cA   cD  th_p+   th_p-   th_pc+   th_pc-  Res+   Res-   th_u+   th_u-    D+     D-
# col springs @ bottom of Story 1 (at base)
rotSect2DModIKModel(3111,Es,Acol_12,Ks_col_1,b,b,Mycol_12,-Mycol_12,LS,LK,LA,LD,cS,cK,cA,cD,th_pP,th_pN,th_pcP,th_pcN,ResP,ResN,th_uP,th_uN,DP,DN);

# recompute strain hardening since Story 2 is not the same height as Story 1
a_mem=(n+1.0)*(Mycol_12*(McMy-1.0)) / (Ks_col_2*th_pP); # strain hardening ratio of spring
b=(a_mem)/(1.0+n*(1.0-a_mem));                          # modified strain hardening ratio of spring (Ibarra & Krawinkler 2005, note: there is mistake in Eqn B.5)
# col springs @ bottom of Story 2 (above Floor 2)
rotSect2DModIKModel(3121,Es,Acol_12,Ks_col_2,b,b,Mycol_12,-Mycol_12,LS,LK,LA,LD,cS,cK,cA,cD,th_pP,th_pN,th_pcP,th_pcN,ResP,ResN,th_uP,th_uN,DP,DN);

# define beam springs
# Spring ID: "4xya", where 4 = beam spring, x = Bay #, y = Floor #, a = location in bay
# "a" convention: 1 = left end, 2 = right end
# redefine the rotations since they are not the same
th_pP=0.02;
th_pN=0.02;
th_pcP=0.16;
th_pcN=0.16;
a_mem=(n+1.0)*(Mybeam_23*(McMy-1.0)) / (Ks_beam_23*th_pP);  # strain hardening ratio of spring
b=(a_mem)/(1.0+n*(1.0-a_mem));                              # modified strain hardening ratio of spring (Ibarra & Krawinkler 2005, note: there is mistake in Eqn B.5)
#beam springs at Floor 2
rotSect2DModIKModel(4121,Es,Abeam_23,Ks_beam_23,b,b,Mybeam_23,-Mybeam_23,LS,LK,LA,LD,cS,cK,cA,cD,th_pP,th_pN,th_pcP,th_pcN,ResP,ResN,th_uP,th_uN,DP,DN);

[silviamazzoni]

this is what i get for Bilin: HystereticSM: IMKBilin:

(they may not be all from the same case)

[amaelkady]

@silviamazzoni Hi silvia, I suspect your issue is related to the definition of the lamda parameter which may be different between the two models. For the IMKBilin documentation (https://outlook.office.com/mail/inbox/id/AAQkAGY4OWI2ODZjLWZiYTEtNGZmOS05ODdiLWU1YjdmNWYwOTdlOAAQAKzbKkuSQk%2FVm00WV%2Fvcq6Y%3D), it mentions that "Lamda is used to compute the reference energy based on the following equation Ref_Energy = Lamda * Fy" while perhaps in theis HystereticSM model, the Ref_Energy is taken directly as Lamda.

[silviamazzoni]

@amaelkady not sure about which link you wanted to paste. I am able to replicate the Bilin behavior well enough with HystereticSM. I am unable to do it with IMKBilin. why is it not yielding? also, what is the meaning of c_A in your table in the documentation, but not in the command? This is the documentation page I am using for IMKBilin: https://opensees.github.io/OpenSeesDocumentation/user/manual/material/uniaxialMaterials/IMKBilin.html

here are my inputs for the two materials: Bilin $matTag $K0 $as_Plus $as_Neg $My_Plus $My_Neg $Lamda_S $Lamda_C $Lamda_A $Lamda_K $c_S $c_C $c_A $c_K $theta_p_Plus $theta_p_Neg $theta_pc_Plus $theta_pc_Neg $Res_Pos $Res_Neg $theta_u_Plus $theta_u_Neg $D_Plus $D_Neg <$nFactor> Bilin 3111 42746000.0 0.00096 0.00096 20350.0 -20350.0 1000.0 1000.0 1000.0 1000.0 1.0 1.0 1.0 1.0 0.025 0.025 0.3 0.3 0.4 0.4 0.4 0.4 1.0 1.0

IMKBilin $matTag $Ke $dp_pos $dpc_pos $du_pos $Fy_pos $FmaxFy_pos $FresFy_pos $dp_neg $dpc_neg $du_neg $Fy_neg $FmaxFy_neg $FresFy_neg $Lamda_S $Lamda_C $Lamda_K $c_S $c_C $c_K $D_pos $D_neg IMKBilin 13111 42746000.0 0.025 0.3 0.4 20350.0 21377.28 8140.0 0.025 0.3 0.4 -20350.0 -19322.71888802286 -8140.0 0.5 0.5 0.5 1 1 1 1 1

[amaelkady]

@silviamazzoni there is no c_A in the IMKBilin definition or the parameter documentation table (see below). There is a difference in the lamda definition as I mentioned in my comment above which is also mentioned as note in the documentation (see highlighted text below). In the IMKBilin, "Lamda is used to compute the reference energy based on the following equation Ref_Energy = Lamda * Fy" while perhaps in the older models model, the Ref_Energy was taken directly as Lamda.

[silviamazzoni]

@amaelkady

1. as i had stated in the original post, the example on that page does have c_A:

but, the main question I asked is:

2. In the input I sent you I am using Lambda=0.5, as recommended in the example. but why does the model not replicate the original Bilin, what am I doing wrong? i like how you not define the maximum strength, and not the tangent, but I can't even get it to yield. 20350./42746000.0 = 0.000476, so the 3rd plot above. should show yielding... <<<<

[amaelkady]

@silviamazzoni 1. for c_A, this just on the parameters used in the IMK models (pinching and peakoriented). It is only in this example script (not the actual documentation of the inputs) becasue the same example is shown in all other models. anayway, this is the deteriroration exponent for accelerated stiffness degradation.

2. I understand the issue but I'm not able to check the code above because it's too long/complicated and I can't tell what parameters went into the model. If possible, can you create a simple script with just one zero length element.

[amaelkady]

@silviamazzoni actually, i think i know what the issue might be: in your definition you have -ve values. The IMKBilin model only takes positive values even for the negative backbone

[amaelkady]

@silviamazzoni also it seems you are not using proper values for FmaxFy_pos (generaly 1 to 1.5)

[silviamazzoni]

yes the code was there to show you how i computed the values in case you were curious, i gave you the actual input. good catch on the negative numbers. I will change them and test it. i would recommend that you make your code a bit more user-resistant and just convert them to positive if negative has no meaning.

[silviamazzoni]

@silviamazzoni also it seems you are not using proper values for FmaxFy_pos (generaly 1 to 1.5) for positive I have 21k/20k, which is 1.05. good catch on the negative one, as yes it goes down. i can fix that.

[mhscott]

yes the code was there to show you how i computed the values in case you were curious, i gave you the actual input. good catch on the negative numbers. I will change them and test it. i would recommend that you make your code a bit more user-resistant and just convert them to positive if negative has no meaning.

I agree. Just change the signs in the constructor. It's an easy modification @amaelkady

[silviamazzoni]

here is the new input, positive and negative look the same, i still don't get yielding: IMKBilin 13111 42746000.0 0.025 0.3 0.4 20350.0 21377.28111197714 8140.0 0.025 0.3 0.4 20350.0 21377.28111197714 8140.0 0.5 0.5 0.5 1 1 1 1 1

[silviamazzoni]

@amaelkady p.s. i do like your format better than the original mess, which is why i'm trying to get it to work.

[amaelkady]

@silviamazzoni thanks. Just to clarify, the model is working and is thoroughly validated. Your definition is still incorrect as the FmaxFy_pos is still incorrect. you should use 1.05 instead

[amaelkady]

@silviamazzoni and also FresFy_pos is a ratio not an absolute force value. this needds to be modified as well in your definition

[mhscott]

@amaelkady It would be helpful if this guidance was added to the documentation page

[amaelkady]

@mhscott it is already mentioned

[silviamazzoni]

oooo, i see, It's a ratio, I was wondering why it had Fy in the name. sorry, I was looking for the word ratio. Yes, i know the model is working no question of it, which is why I asked what I was doing wrong. I have adhd I can't read all the words in documentation. I count on the names being intuitive.

[silviamazzoni]

so this makes sense on the ratio sign always being positive! and why I don't see the yield since the post-yield slope is so high! :)

[mhscott]

@amaelkady I don't see it in the parameter list, but I'll take your word for it. Thanks.

[amaelkady]

the notes section under the table

[silviamazzoni]

okey, i got it to work. I had to use lambda = 0 to match the old stuff:

yes everything is there, I just don't have the neurological ability to read it.

thank you!

[mhscott]

@amaelkady Thanks, that's very helpful. I was also referring to the recommended ranges of parameter values, like for the Fy ratios and the c values and what not. If you can add those to the parameter list, that would be very helpful to users not familiar with all the intimate details of the models.