Brainery Byte
R210706
Silvia Mazzoni, 2021
[email protected]
NOTES:
Process Outline:
When you are done, you can print the notebook. You can also download the notebook to your local computer.
Because we are working in Binder, and Binder sessions are meant to be ephemeral, it is not possible for you to save any changes you make to your Jupyter Notebook. If you do make changes or notes, you will need to download the notebook to your own computer by clicking File > Download as > Notebook (.ipynb). The only way you will be able to run these is if you have the appropriate software to run Jupyter Notebooks in Python and pip install OpenSeesPy and eSEESminiPy in your Python configuration. You may view my videos on how to install Anaconda, Jupyter Notebooks and OpenSeesPy (https://www.youtube.com/c/silviasbrainery).
This Code has been developed by Silvia Mazzoni. Please acknowledge this in your scripts, when applicable
# Initalize Python Libraries
import openseespy.opensees as ops
import eSEESminiPy
import matplotlib.pyplot as plt
import math
%%javascript
// expand output box so graphs are viewed in notebook properly
IPython.OutputArea.auto_scroll_threshold = 1000;
# define UNITS ----------------------------------------------------------------------------
inch = 1. # define basic units -- output units
kip = 1. # define basic units -- output units
sec = 1. # define basic units -- output units
LunitTXT = 'inch' # define basic-unit text for output
FunitTXT = 'kip' # define basic-unit text for output
TunitTXT = 'sec' # define basic-unit text for output
ft = 12.*inch # define engineering units
ksi = kip/math.pow(inch,2)
psi = ksi/1000.
lbf = psi*inch*inch # pounds force
pcf = lbf/math.pow(ft,3) # pounds per cubic foot
psf = lbf/math.pow(ft,3) # pounds per square foot
inch2 = inch*inch # inch^2
inch4 = inch*inch*inch*inch # inch^4
cm = inch/2.54 # centimeter, needed for displacement input in MultipleSupport excitation
PI = 2*math.asin(1.0) # define constants
g = 32.2*ft/math.pow(sec,2) # gravitational acceleration
Ubig = 1.e10 # a really large number
Usmall = 1/Ubig # a really small number
inchInv = 1/inch;
# Initialize Arrays and variables
global matTag
matTag = 0
secTag = 0
SectionTagMap = {}
MaterialTagMap = {}
allSectionData = {}
WsectionData = {}
MomCurvData = {}
MatStressStrainData = {}
SectionLabelList = []
plt.close('all')
# Define bending directions for moment-curvature analysis
# you may remove either y or z if needed.
BendingDirectionList = ['z','y']
# UTILITIES
# Define some local handy procs
def addMaterial(MaterialType,InputArray):
global matTag
matTag += 1
ops.uniaxialMaterial(MaterialType,matTag,*InputArray)
return matTag
### Plot Moment-Curvature Responses
def plotMomCurv(testSectionLabel,MomCurvData):
thisFig = {}
axModel = {}
thisFigMomCurv = plt.figure('Moment Curvature ' + testSectionLabel,figsize=(6,3), dpi=200, facecolor='w', edgecolor='k' )
iAx = 0
for BendingDirection in MomCurvData[testSectionLabel].keys():
iAx += 1
axMomCurv = thisFigMomCurv.add_subplot(1,len(MomCurvData[testSectionLabel].keys()),iAx)
line, = axMomCurv.plot(MomCurvData[testSectionLabel][BendingDirection]['Curvature'], MomCurvData[testSectionLabel][BendingDirection]['Moment'])
eSEESminiPy.formatAx(axMomCurv,'Moment-Curvature ' + testSectionLabel+ ' ' + BendingDirection,'Curvature ' + BendingDirection,'Moment ' + BendingDirection,6,6)
for thisFiberLabel in MomCurvData[testSectionLabel][BendingDirection]['StrainData'][SectionTagMap[testSectionLabel]].keys():
if iAx == 1:
thisFig[thisFiberLabel] = plt.figure('Stress Strain ' + thisFiberLabel,figsize=(4,1), dpi=200, facecolor='w', edgecolor='k' )
axModel[thisFiberLabel] = thisFig[thisFiberLabel].add_subplot(1,len(MomCurvData[testSectionLabel].keys()),iAx)
line, = axModel[thisFiberLabel].plot(MomCurvData[testSectionLabel][BendingDirection]['StrainData'][SectionTagMap[testSectionLabel]][thisFiberLabel],MomCurvData[testSectionLabel][BendingDirection]['StressData'][SectionTagMap[testSectionLabel]][thisFiberLabel], 'k',linewidth = 1)
eSEESminiPy.formatAx(axModel[thisFiberLabel],'Stress Strain ' + thisFiberLabel + " " + BendingDirection,'Strain ' + BendingDirection,'Stress ' + BendingDirection,6,6)
# Plot P-M Interaction Data
def plotPMInterax(testSectionLabel,PMinteraxData):
FiberFig = {}
FiberAx = {}
thisFigPMI = plt.figure('P-M Interax ' + testSectionLabel,figsize=(6,3), dpi=200, facecolor='w', edgecolor='k' )
thisFigMomCurv = plt.figure('Moment Curvature ' + testSectionLabel,figsize=(6,3), dpi=200, facecolor='w', edgecolor='k' )
iAx = 0
for BendingDirection in BendingDirectionList:
iAx += 1
axPMI = thisFigPMI.add_subplot(1,2,iAx)
axMomCurv = thisFigMomCurv.add_subplot(1,2,iAx)
for CurvatureSign in PMinteraxData[testSectionLabel][BendingDirection].keys():
for thisAxial in PMinteraxData[testSectionLabel][BendingDirection][CurvatureSign].keys():
thisMomCurvData = PMinteraxData[testSectionLabel][BendingDirection][CurvatureSign][thisAxial]
line, = axPMI.plot(thisMomCurvData['Moment'], thisMomCurvData['AxialForce'])
axPMI.plot(max(thisMomCurvData['Moment'], key=abs), thisAxial,'ro',markerSize = 2,linewidth = 1)
line, = axMomCurv.plot(thisMomCurvData['Curvature'], thisMomCurvData['Moment'])
StrainData = thisMomCurvData['StrainData'][SectionTagMap[testSectionLabel]]
StressData = thisMomCurvData['StressData'][SectionTagMap[testSectionLabel]]
for thisFiberLabel,thisFiberStrain in StrainData.items():
if not thisFiberLabel in FiberFig.keys():
FiberFig[thisFiberLabel]= plt.figure('Fiber SS ' + testSectionLabel + thisFiberLabel,figsize=(4,1), dpi=200, facecolor='w', edgecolor='k' )
if not thisFiberLabel+ ' ' + BendingDirection in FiberAx.keys():
FiberAx[thisFiberLabel+ ' ' + BendingDirection] = FiberFig[thisFiberLabel].add_subplot(1,2,iAx)
thisFiberStress = StressData[thisFiberLabel]
line, = FiberAx[thisFiberLabel+ ' ' + BendingDirection].plot(thisFiberStrain,thisFiberStress,linewidth = 1)
eSEESminiPy.formatAx(axPMI,'P-M Interax ' + testSectionLabel + ' ' + BendingDirection,'Bending Moment ' + BendingDirection,'Axial Force',6,6)
eSEESminiPy.formatAx(axMomCurv,'P-M Interax MomentCurv ' + testSectionLabel + ' ' + BendingDirection ,'Curvature ' + BendingDirection,'Bending Moment ' + BendingDirection,6,6)
eSEESminiPy.reverseYaxis(axPMI) # plot Axial Force Compression up
for thisAxLabel,thisValue in FiberAx.items():
eSEESminiPy.formatAx(thisValue,thisAxLabel,'Strain','Stress',6,6)
# Define Normalized Loading Cycles for Material and Section Testing.
# These cycles are normalized to have a maximum value of 1.0
# You can later change the amplitude factor depending on the application
peaksArray = [0.1,0.2,0.3,0.5,0.7,0.85,1.0]
nCycles = 3
nSteps = 200
StrainSeries = eSEESminiPy.defineStrainHistory(peaksArray,nCycles,nSteps)
eSEESminiPy.plotYoneList(StrainSeries,title='Strain History',Ylabel='Strain')
# for concrete you want more loading in compression for testing
peaksArray2 = [0.1,0,-0.2,0,0.15,0,-0.3,0,0.3,0,-0.6,0,-0.8,0,-1.0,0.0]
nCycles2 = 0; # do a single segment from one value to the next
nSteps = 100
StrainSeries2 = eSEESminiPy.defineStrainHistory(peaksArray2,nCycles2,nSteps)
eSEESminiPy.plotYoneList(StrainSeries2,title='Strain History 2',Ylabel='Strain')
# Start-up OpenSees
ops.wipe()
Dim = '3D'
ops.model('basic','-ndm',3,'-ndf',6)
Fiber Section, Circular RC Section
################################################################################################
### REINFORCED-CONCRETE SECTION -- Circular
################################################################################################
MainSectionLabel = 'RCsectionCirc'
thisSectionType = '_FiberOnly'
thisSectionLabel = MainSectionLabel + thisSectionType
# initialize Section Data
thisSectionData = {}
thisMaterialTagMap = {}
## MATERIALS
# ##### Reinforcing steel:
# Define Steel Material
Fy=60*ksi
YoungsModulus = 2.900E+04*ksi
Bsh = 0.01
# Backbone points for Hysteretic Material
f1 = Fy
eps1 = Fy/YoungsModulus
eps2 = 20*eps1
f2 = f1+Bsh*(eps2-eps1)*(YoungsModulus)
f3 = 0.9*f2
eps3 = 50*eps1
PositiveEnvelopeData = [f1,eps1,f2,eps2,f3,eps3]
NegativeEnvelopeData = [-f1,-eps1,-f2,-eps2,-f3,-eps3]
# hyteretic-material properties
[pinchX,pinchY] = [1,1]
[damage1,damage2] = [0.0,0.0]
beta = 0.5
hystereticData = [pinchX,pinchY,damage1,damage2,beta]
thisMaterialTagMap['Steel'] = addMaterial('Hysteretic',[*PositiveEnvelopeData,*NegativeEnvelopeData,*hystereticData])
# ##### Concrete:
# nominal concrete compressive strength
fc = -4.0*ksi # CONCRETE Compressive Strength, ksi (+Tension, -Compression)
Ec = 57*ksi*math.sqrt(-fc/psi) # Concrete Elastic Modulus
nu = 0.2
ShearModulus = Ec/2./(1+nu) # Torsional stiffness Modulus
# confined concrete
Kfc = 1.3 # ratio of confined to unconfined concrete strength
Kres = 0.2 # ratio of residual/ultimate to maximum stress
fc1C = Kfc*fc # CONFINED concrete (mander model), maximum stress
eps1C = 2.*fc1C/Ec # strain at maximum stress
fc2C = Kres*fc1C # ultimate stress
eps2C = 20*eps1C # strain at ultimate stress
Lambda = 0.1 # ratio between unloading slope at eps2 and initial slope Ec
# unconfined concrete
fc1U = fc # UNCONFINED concrete (todeschini parabolic model), maximum stress
eps1U = -0.003 # strain at maximum strength of unconfined concrete
fc2U = Kres*fc1U # ultimate stress
eps2U = -0.01 # strain at ultimate stress
# tensile-strength properties
ftC = -0.14*fc1C # tensile strength+tension
ftU = -0.14*fc1U # tensile strength+tension
Ets = ftU/0.002 # tension softening stiffness
# Define Materials:
thisMaterialTagMap['UnconfinedConcrete'] = addMaterial('Concrete02',[fc1U,eps1U,fc2U,eps2U,Lambda,ftU,Ets])
thisMaterialTagMap['ConfinedConcrete'] = addMaterial('Concrete02',[fc1C,eps1C,fc2C,eps2C,Lambda,ftC,Ets])
thisMaterialTagMap['ConfinedConcrete2'] = addMaterial('Concrete02',[1.05*fc1C,eps1C,1.05*fc2C,eps2C,Lambda,ftC,Ets])
# Add material maps to all
MaterialTagMap[thisSectionLabel] = thisMaterialTagMap
# Test Materials
thisMatStressStrainData = {}
# you may have some convergence problems here because I can't add a stabilizing elastic element in the procedure.
# use symmetric loading for Steel
ampFactor = 0.1
thisMatStressStrainData['Steel'] = eSEESminiPy.getMaterialStressStrain(Dim,thisMaterialTagMap['Steel'],StrainSeries,ampFactor)
# use more-compression loading for concrete
ampFactor = 0.003
thisMatStressStrainData['UnconfinedConcrete'] = eSEESminiPy.getMaterialStressStrain(Dim,thisMaterialTagMap['UnconfinedConcrete'],StrainSeries2,ampFactor)
thisMatStressStrainData['ConfinedConcrete'] = eSEESminiPy.getMaterialStressStrain(Dim,thisMaterialTagMap['ConfinedConcrete'],StrainSeries2,ampFactor)
thisMatStressStrainData['ConfinedConcrete2'] = eSEESminiPy.getMaterialStressStrain(Dim,thisMaterialTagMap['ConfinedConcrete2'],StrainSeries2,ampFactor)
for thisMaterialLabel in thisMatStressStrainData.keys():
eSEESminiPy.plotXYoneList(thisMatStressStrainData[thisMaterialLabel]['strain'],thisMatStressStrainData[thisMaterialLabel]['stress'],thisMaterialLabel)
# Add material response to all
MatStressStrainData[thisSectionLabel] = thisMatStressStrainData
# FIBER SECTION properties
# assign a section Tag for OpenSees and save it in the data map
secTag += 1
SectionTagMap[thisSectionLabel] = secTag
thisSectionData['SectionTag'] = secTag
# Column section geometry:
DCol = 28*inch
DColInner = 6*inch
thisSectionData['Ag'] = math.pi*DCol*DCol/4 - math.pi*DColInner*DColInner/4
cover = 2*inch
layerspacing = 2.0*inch
numBarsOuter = 22 # number of longitudinal-reinforcement bars
numBarsInner = 12
barAreaOuter = 1.*inch*inch # longitudinal-reinforcement bar area
barAreaInner = 1.*inch*inch
# axial strength of column, Compression and tension
thisSectionData['PmaxC'] = thisSectionData['Ag']*fc - (numBarsOuter*barAreaOuter + numBarsInner*barAreaInner)*Fy
thisSectionData['PmaxT'] = 0.99*(numBarsOuter*barAreaOuter + numBarsInner*barAreaInner)*Fy
nfCoreR = 8 # number of radial divisions in the core (number of "rings")
nfCoreT = 8 # number of theta divisions in the core (number of "wedges")
nfCoreMidR = 4 # number of radial divisions in the cover
nfCoreMidT = 16 # number of theta divisions in the cover
nfCoverR = 4 # number of radial divisions in the cover
nfCoverT = 16 # number of theta divisions in the cover
# Torsional Stiffness
GJ = 1e9
thisSectionData['GJ'] = GJ
[coreID,coverID,steelID] = [thisMaterialTagMap['ConfinedConcrete'],thisMaterialTagMap['UnconfinedConcrete'],thisMaterialTagMap['Steel']]
Rinner = DColInner/2
Router = DCol/2
Rcore = Router - cover
Rcore2 = Router - cover - layerspacing
# Define the fiber section
FiberSectionCommandData = {}
FiberSectionCommandData['sectionCall'] = ['section','Fiber',secTag,'-GJ',GJ]
# Define the core patch
FiberSectionCommandData['CorePatch'] = ['patch','circ',coreID,nfCoreT,nfCoreR,0,0,Rinner,Rcore2,0,360]
# Define the cover patch
FiberSectionCommandData['MidPatch'] = ['patch','circ',coreID,nfCoreMidT,nfCoreMidR,0,0,Rcore2,Rcore,0,360]
FiberSectionCommandData['CoverPatch'] = ['patch','circ',coverID,nfCoverT,nfCoverR,0,0,Rcore,Router,0,360]
# Define reinforcing layers
thetaOuter = 360.0/numBarsOuter # Determine angle increment between bars for the start of the layer, otherwise you'd be double counting
thetaInner = 360.0/numBarsInner # Determine angle increment between bars for the start of the layer, otherwise you'd be double counting
FiberSectionCommandData['OuterReinf'] = ['layer','circ',steelID,numBarsOuter,barAreaOuter,0,0,Rcore,0,360-thetaOuter]
FiberSectionCommandData['InnerReinf'] = ['layer','circ',steelID,numBarsInner,barAreaInner,0,0,Rcore2,0,360-thetaInner]
# Add recorder fibers (zero area) and define strain limit states
strainLimitsSteel = {}
strainLimitsSteel['SteelYield'] = [-1.1*Fy/YoungsModulus,1.1*Fy/YoungsModulus]
strainLimitsSteel['LS'] = [-4*Fy/YoungsModulus,4*Fy/YoungsModulus]
strainLimitsSteel['CP'] = [-8*Fy/YoungsModulus,8*Fy/YoungsModulus]
strainLimitsConcrete = {}
strainLimitsConcrete['LS'] = [-2*0.004,999]
strainLimitsConcrete['CP'] = [-4*0.004,999]
strainLimitsConcrete['ConcreteCrushing'] = [-0.004,999]
FiberSectionCommandData['BottomSteel'] = ['fiber',-Rcore,0,0,steelID,strainLimitsSteel]
FiberSectionCommandData['TopSteel'] = ['fiber',Rcore,0,0,steelID,strainLimitsSteel]
FiberSectionCommandData['RightSteel'] = ['fiber',0,Rcore,0,steelID,strainLimitsSteel]
FiberSectionCommandData['LeftSteel'] = ['fiber',0,-Rcore,0,steelID,strainLimitsSteel]
FiberSectionCommandData['BottomCore'] = ['fiber',-Rcore,0,0,coreID,strainLimitsConcrete]
FiberSectionCommandData['TopCore'] = ['fiber',Rcore,0,0,coreID,strainLimitsConcrete]
FiberSectionCommandData['RightCore'] = ['fiber',0,Rcore,0,coreID,strainLimitsConcrete]
FiberSectionCommandData['LeftCore'] = ['fiber',0,-Rcore,0,coreID,strainLimitsConcrete]
FiberSectionCommandData['Origin'] = ['fiber',Rcore,0,0,coreID]
# assign all these data to section
thisSectionData['FiberSectionCommandData'] = FiberSectionCommandData
# draw Fiber Section:
eSEESminiPy.drawFiberSection(FiberSectionCommandData,thisSectionLabel)
# buildSection:
ops.section(*FiberSectionCommandData['sectionCall'][1:])
for thisFiberLabel,thisFiberCommandData in FiberSectionCommandData.items():
fiberType = thisFiberCommandData[0]
if fiberType == 'patch':
ops.patch(*thisFiberCommandData[1:12+1])
elif fiberType == 'layer':
if thisFiberCommandData[3]>0:
ops.layer(*thisFiberCommandData[1:8+1])
elif fiberType == 'fiber':
ops.fiber(*thisFiberCommandData[1:4+1])
# add to all
allSectionData[thisSectionLabel] = thisSectionData
## RUN Moment-Curvature Analysis
testSectionLabel = MainSectionLabel + thisSectionType
MomCurvData[testSectionLabel] = {}
# The InputSection is a single argument to the Moment-Curvature Analysis
# If you want strain-recorder data, include the FiberSectionDommandData into a list:
FiberSectionCommandData = allSectionData[testSectionLabel]['FiberSectionCommandData']
InputSection = [SectionTagMap[testSectionLabel],FiberSectionCommandData]
### MOMENT-CURVATURE ANALYSIS
axialLoad = 180
ampFactor = 0.01
thisLimitState = '' # no limit state as we are doing a cyclic analysis
for BendingDirection in BendingDirectionList:
MomCurvData[testSectionLabel][BendingDirection] = eSEESminiPy.runMomentCurvature(InputSection,axialLoad,StrainSeries,5,Dim+BendingDirection,ampFactor,thisLimitState)
### Plot Moment-Curvature Responses
plotMomCurv(testSectionLabel,MomCurvData)