{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Moment Curvature Example\n",
"This example is based on https://opensees.berkeley.edu/wiki/index.php/Moment_Curvature_Example and https://openseespydoc.readthedocs.io/en/latest/src/MomentCurvature.html but using the Reinforcing Steel Material (A706 Grade 60) develop by Moehle https://opensees.berkeley.edu/wiki/index.php/Reinforcing_Steel_Material\n",
"\n",
"For your own moment-curvature analyses take into account that the cover should be outsude of the diameter of the rebar including the transverse reinforcement. In this way the example in the OpenSees wiki doesn't represent a realistic scenario.\n",
"\n",
"The properties of the material A706 grade 60 are the following (Final Report - Performance Characterization of Beams with High-Strength Reinforcement https://pdfs.semanticscholar.org/d74a/710b727f26a06c72ac8729d904fba95fb75e.pdf):\n",
"\n",
"Steel: Uniaxial material ReinforcingSteel (Kunnath et al. 2009) based on Chang and Mander (1994)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"| **Model** |**Variable**|**Description**|**Value**|\n",
"|:--:|---:|-----:|:----:|\n",
"|**ReinforcingSteel**| | | |\n",
"| |$f_y$|Yield stress (ksi)|$64.5$|\n",
"| |$f_u$|Ultimate stress (ksi)|$95.5$|\n",
"| |$E_s$|Young’s modulus (ksi)|$29000$|\n",
"| |$E_{sh}$|Tangent stiffness at initiation of strain hardening (ksi)|$950$|\n",
"| |$\\varepsilon_{sh}$|Strain at initiation of strain hardening|$0.0055$|\n",
"| |$\\varepsilon_{su}$|Strain at ultimate stress|$0.15$|\n",
"| |$l_{RS},$ $\\beta,$ $r,$ $\\gamma,$|Parameters for buckling model based on Gomes and Appleton (1997)|$5, 1.0, 0.75, 0.0$|\n",
"| |$a_1$, $limit$|Parameters for controlling isotropic hardening|$4.3, 0.01$|\n",
"| |$R1,$ $R2,$ $R3$|Parameters for controlling transition from elastic to plastic branches|$0.333, 20, 6$|"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from openseespy.opensees import *\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 1: ReinforcingSteel\n",
"\n",
"Using the properties of the [table](#Moment-Curvature-Example) I'll try to replicate the section of the [OpenSees example](#Moment-Curvature-Example)."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Start MomentCurvature.py example\n",
"Estimated yield curvature: 0.0001318007662835249\n",
"Estimated ultimate curvature: 0.0033609195402298894\n",
"==========================\n"
]
}
],
"source": [
"def MomentCurvature(secTag, axialLoad, maxK, numIncr=100):\n",
" \n",
" # Define two nodes at (0,0)\n",
" node(1, 0.0, 0.0)\n",
" node(2, 0.0, 0.0)\n",
"\n",
" # Fix all degrees of freedom except axial and bending\n",
" fix(1, 1, 1, 1)\n",
" fix(2, 0, 1, 0)\n",
" \n",
" # Define element\n",
" # tag ndI ndJ secTag\n",
" element('zeroLengthSection', 1, 1, 2, secTag)\n",
"\n",
" # Define constant axial load\n",
" timeSeries('Constant', 1)\n",
" pattern('Plain', 1, 1)\n",
" load(2, axialLoad, 0.0, 0.0)\n",
"\n",
" # Define analysis parameters\n",
" integrator('LoadControl', 0.0)\n",
" system('SparseGeneral', '-piv')\n",
" test('NormUnbalance', 1e-9, 500)\n",
" numberer('Plain')\n",
" constraints('Plain')\n",
" #algorithm('Newton')\n",
" algorithm('ModifiedNewton', 'False', 'False')\n",
" analysis('Static')\n",
"\n",
" # Do one analysis for constant axial load\n",
" analyze(1)\n",
"\n",
" # Define reference moment\n",
" timeSeries('Linear', 2)\n",
" pattern('Plain',2, 2)\n",
" load(2, 0.0, 0.0, 1.0)\n",
"\n",
" # Compute curvature increment\n",
" dK = maxK / numIncr\n",
"\n",
" # Use displacement control at node 2 for section analysis\n",
" integrator('DisplacementControl', 2,3,dK,1,dK,dK)\n",
" \n",
" # Create recorder\n",
" #recorder Node -file section$secTag.out -time -node 2 -dof 3 disp\n",
" import os\n",
" recorder('Node', '-file', 'node2.txt','-time', '-closeOnWrite', '-node', *[2],'-dof', *[3], 'disp')\n",
" recorder('Element', '-file', 'node13.txt', '-closeOnWrite', '-ele', *[1],*['section',1,str(10.5),str(0),\"1\",\"stressStrain\"])\n",
" recorder('Element', '-file', 'node14.txt', '-closeOnWrite', '-ele', *[1],*['section',1,str(10.5),str(0),\"2\",\"stressStrain\"])\n",
" # Do the section analysis\n",
" analyze(numIncr)\n",
" \n",
"wipe()\n",
"print(\"Start MomentCurvature.py example\")\n",
"\n",
"# Define model builder\n",
"# --------------------\n",
"model('basic','-ndm',2,'-ndf',3)\n",
"\n",
"# Define materials for nonlinear columns\n",
"# ------------------------------------------\n",
"# CONCRETE tag f'c ec0 f'cu ecu\n",
"# Core concrete (confined)\n",
"uniaxialMaterial('Concrete01',1, -6.0, -0.0033, -5.0, -0.014)\n",
"\n",
"# Cover concrete (unconfined)\n",
"uniaxialMaterial('Concrete01',2, -5.0, -0.002, 0.0, -0.006)\n",
"\n",
"# STEEL\n",
"# Reinforcing steel \n",
"fy = 64.5; # Yield stress\n",
"fu=95.5;\n",
"Es = 29000.0; # Young's modulus\n",
"Esh=950;\n",
"esh=0.0055;\n",
"eult=0.15;\n",
"lsr=2;\n",
"beta=1;\n",
"r=1;\n",
"gamma=0;\n",
"a1=4.3;\n",
"limit=0.01\n",
"R1=0.333;\n",
"R2=20;\n",
"R3=6\n",
"\n",
"# tag fy E0 b\n",
"#uniaxialMaterial('Steel01', 3, fy, Es, 1/100)\n",
"uniaxialMaterial('ReinforcingSteel', 3, fy, fu, Es, Esh, esh, eult, '-GABuck', lsr, beta, r, gamma, '-IsoHard', a1, limit, '-MPCurveParams', R1, R2, R3)\n",
"#uniaxialMaterial('ReinforcingSteel', 3, fy, fu, Es, Esh, esh, eult)\n",
"#uniaxialMaterial('ReinforcingSteel', 3, fy, fu, Es, Esh, esh, eult, '-IsoHard', a1, limit, '-MPCurveParams', R1, R2, R3)\n",
"#uniaxialMaterial('Steel02', 3, fy, Es, 0.01,*[15,0.925,0.15])\n",
"\n",
"# Define cross-section for nonlinear columns\n",
"# ------------------------------------------\n",
"\n",
"# set some paramaters\n",
"colWidth = 15\n",
"colDepth = 24\n",
"\n",
"cover = 1.5\n",
"As = 0.60; # area of no. 7 bars\n",
"\n",
"# some variables derived from the parameters\n",
"y1 = colDepth/2.0\n",
"z1 = colWidth/2.0\n",
"\n",
"section('Fiber', 1)\n",
"\n",
"# Create the concrete core fibers\n",
"patch('rect',1,20*2,2*2 ,cover-y1, cover-z1, y1-cover, z1-cover)\n",
"\n",
"# Create the concrete cover fibers (top, bottom, left, right)\n",
"patch('rect',2,20*2,2*2 ,-y1, z1-cover, y1, z1)\n",
"patch('rect',2,20*2,2*2 ,-y1, -z1, y1, cover-z1)\n",
"patch('rect',2,4*2,2*2 ,-y1, cover-z1, cover-y1, z1-cover)\n",
"patch('rect',2,4*2,2*2 ,y1-cover, cover-z1, y1, z1-cover)\n",
"\n",
"# Create the reinforcing fibers (left, middle, right)\n",
"layer('straight', 3, 3, As, y1-cover, z1-cover, y1-cover, cover-z1)\n",
"layer('straight', 3, 2, As, 0.0 , z1-cover, 0.0 , cover-z1)\n",
"layer('straight', 3, 3, As, cover-y1, z1-cover, cover-y1, cover-z1)\n",
"\n",
"# Estimate yield curvature\n",
"# (Assuming no axial load and only top and bottom steel)\n",
"# d -- from cover to rebar\n",
"d = colDepth-cover\n",
"# steel yield strain\n",
"epsy = fy/Es\n",
"Ky = epsy/(0.75*d)\n",
"\n",
"# Print estimate to standard output\n",
"print(\"Estimated yield curvature: \", Ky)\n",
"\n",
"# Set axial load \n",
"P = -180\n",
"\n",
"# Target ductility for analysis\n",
"mu = 25.5\n",
"\n",
"# Number of analysis increments\n",
"numIncr = 100\n",
"\n",
"# Call the section analysis procedure\n",
"MomentCurvature(1, P, Ky*mu, numIncr)\n",
"\n",
"results = open('results.out','a+')\n",
"\n",
"u = nodeDisp(2,3)\n",
"print(\"Estimated ultimate curvature: \", u)\n",
"#if abs(u-0.00190476190476190541)<1e-12:\n",
"# results.write('PASSED : MomentCurvature.py\\n');\n",
"# print(\"Passed!\")\n",
"#else:\n",
"# results.write('FAILED : MomentCurvature.py\\n');\n",
"# print(\"Failed!\")\n",
"\n",
"results.close()\n",
"\n",
"print(\"==========================\")"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\r\n",
"\r\n",
"\r\n",
"\r\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\r\n",
"\r\n",
"\r\n",
"\r\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"from IPython.display import set_matplotlib_formats\n",
"set_matplotlib_formats('svg')\n",
"import matplotlib.font_manager as font_manager\n",
"font = font_manager.FontProperties(family='Times New Roman',\n",
" size=10)\n",
"\n",
"with open('node2.txt', 'r') as f:\n",
" lines = f.readlines()\n",
" y11 = [float(line.split()[0]) for line in lines]\n",
" x11 = [float(line.split()[1]) for line in lines]\n",
"plt.plot(x11 ,y11,color='red',label='ReinforcingSteel')\n",
"plt.title('Moment Curvature ReinforcingSteel', fontsize=14,fontfamily='Times New Roman',weight='bold')\n",
"plt.xlabel('Curvature (1/in)', fontsize=12,fontfamily='Times New Roman',weight='bold')\n",
"plt.ylabel('Moment (kip-in)', fontsize=12,fontfamily='Times New Roman',weight='bold')\n",
"plt.xticks(rotation=70)\n",
"plt.grid(True)\n",
"plt.show()\n",
"\n",
"with open('node13.txt', 'r') as f:\n",
" lines = f.readlines()\n",
" y13 = [float(line.split()[0]) for line in lines]\n",
" x13 = [float(line.split()[1]) for line in lines]\n",
"plt.plot(x13 ,y13,color='blue',label='Confined')\n",
"\n",
"with open('node14.txt', 'r') as f:\n",
" lines = f.readlines()\n",
" y14 = [float(line.split()[0]) for line in lines]\n",
" x14 = [float(line.split()[1]) for line in lines]\n",
"plt.plot(x14 ,y14,color='green',label='Unconfined')\n",
"\n",
"plt.xlabel('Strain (in/in)', fontsize=12,fontfamily='Times New Roman',weight='bold')\n",
"plt.ylabel('Stress (ksi)', fontsize=12,fontfamily='Times New Roman',weight='bold')\n",
"plt.title('Concrete Material Stress-Strain', fontsize=14,fontfamily='Times New Roman',weight='bold')\n",
"legend = plt.legend(loc='best', shadow=False, fontsize='8',ncol=3,prop=font)\n",
"\n",
"plt.xticks(rotation=70)\n",
"plt.grid(True)\n",
"plt.show()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now using the same method but using a bilinear steel material the same method as XTRACT"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 2: Steel02\n",
"\n",
"Now we will use a bilinear steel material similar in how XTRACT does it"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Start MomentCurvature.py example\n",
"Estimated yield curvature: 0.0001318007662835249\n",
"Estimated ultimate curvature: 0.0033609195402298894\n",
"==========================\n"
]
}
],
"source": [
"wipe()\n",
"print(\"Start MomentCurvature.py example\")\n",
"\n",
"# Define model builder\n",
"# --------------------\n",
"model('basic','-ndm',2,'-ndf',3)\n",
"\n",
"# Define materials for nonlinear columns\n",
"# ------------------------------------------\n",
"# CONCRETE tag f'c ec0 f'cu ecu\n",
"# Core concrete (confined)\n",
"uniaxialMaterial('Concrete01',1, -6.0, -0.0033, -5.0, -0.014)\n",
"\n",
"# Cover concrete (unconfined)\n",
"uniaxialMaterial('Concrete01',2, -5.0, -0.002, 0.0, -0.006)\n",
"\n",
"# tag fy E0 b\n",
"#uniaxialMaterial('Steel01', 3, fy, Es, 1/100)\n",
"#uniaxialMaterial('ReinforcingSteel', 3, fy, fu, Es, Esh, esh, eult, '-GABuck', lsr, beta, r, gamma, '-IsoHard', a1, limit, '-MPCurveParams', R1, R2, R3)\n",
"#uniaxialMaterial('ReinforcingSteel', 3, fy, fu, Es, Esh, esh, eult)\n",
"#uniaxialMaterial('ReinforcingSteel', 3, fy, fu, Es, Esh, esh, eult, '-IsoHard', a1, limit, '-MPCurveParams', R1, R2, R3)\n",
"uniaxialMaterial('Steel02', 3, fy, Es, 0.01,*[15,0.925,0.15])\n",
"\n",
"# Define cross-section for nonlinear columns\n",
"# ------------------------------------------\n",
"\n",
"section('Fiber', 1)\n",
"\n",
"# Create the concrete core fibers\n",
"patch('rect',1,20*2,2*2 ,cover-y1, cover-z1, y1-cover, z1-cover)\n",
"\n",
"# Create the concrete cover fibers (top, bottom, left, right)\n",
"patch('rect',2,20*2,2*2 ,-y1, z1-cover, y1, z1)\n",
"patch('rect',2,20*2,2*2 ,-y1, -z1, y1, cover-z1)\n",
"patch('rect',2,4*2,2*2 ,-y1, cover-z1, cover-y1, z1-cover)\n",
"patch('rect',2,4*2,2*2 ,y1-cover, cover-z1, y1, z1-cover)\n",
"\n",
"# Create the reinforcing fibers (left, middle, right)\n",
"layer('straight', 3, 3, As, y1-cover, z1-cover, y1-cover, cover-z1)\n",
"layer('straight', 3, 2, As, 0.0 , z1-cover, 0.0 , cover-z1)\n",
"layer('straight', 3, 3, As, cover-y1, z1-cover, cover-y1, cover-z1)\n",
"\n",
"# Estimate yield curvature\n",
"# (Assuming no axial load and only top and bottom steel)\n",
"# d -- from cover to rebar\n",
"d = colDepth-cover\n",
"# steel yield strain\n",
"epsy = fy/Es\n",
"Ky = epsy/(0.75*d)\n",
"\n",
"# Print estimate to standard output\n",
"print(\"Estimated yield curvature: \", Ky)\n",
"\n",
"\n",
"# Number of analysis increments\n",
"numIncr = 100\n",
"\n",
"# Call the section analysis procedure\n",
"MomentCurvature(1, P, Ky*mu, numIncr)\n",
"\n",
"results = open('results.out','a+')\n",
"\n",
"u = nodeDisp(2,3)\n",
"print(\"Estimated ultimate curvature: \", u)\n",
"#if abs(u-0.00190476190476190541)<1e-12:\n",
"# results.write('PASSED : MomentCurvature.py\\n');\n",
"# print(\"Passed!\")\n",
"#else:\n",
"# results.write('FAILED : MomentCurvature.py\\n');\n",
"# print(\"Failed!\")\n",
"\n",
"results.close()\n",
"\n",
"print(\"==========================\")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\r\n",
"\r\n",
"\r\n",
"\r\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\r\n",
"\r\n",
"\r\n",
"\r\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"from IPython.display import set_matplotlib_formats\n",
"set_matplotlib_formats('svg')\n",
"with open('node2.txt', 'r') as f:\n",
" lines = f.readlines()\n",
" y22 = [float(line.split()[0]) for line in lines]\n",
" x22 = [float(line.split()[1]) for line in lines]\n",
"plt.plot(x22 ,y22,color='blue',label='Steel02')\n",
"plt.title('Moment Curvature Steel02', fontsize=14,fontfamily='Times New Roman',weight='bold')\n",
"plt.xlabel('Curvature (1/in)', fontsize=12,fontfamily='Times New Roman',weight='bold')\n",
"plt.ylabel('Moment (kip-in)', fontsize=12,fontfamily='Times New Roman',weight='bold')\n",
"plt.xticks(rotation=70)\n",
"plt.grid(True)\n",
"plt.show()\n",
"\n",
"plt.show()\n",
"\n",
"with open('node13.txt', 'r') as f:\n",
" lines = f.readlines()\n",
" y13 = [float(line.split()[0]) for line in lines]\n",
" x13 = [float(line.split()[1]) for line in lines]\n",
"plt.plot(x13 ,y13,color='blue',label='Confined')\n",
"\n",
"with open('node14.txt', 'r') as f:\n",
" lines = f.readlines()\n",
" y14 = [float(line.split()[0]) for line in lines]\n",
" x14 = [float(line.split()[1]) for line in lines]\n",
"plt.plot(x14 ,y14,color='green',label='Unconfined')\n",
"\n",
"plt.xlabel('Strain (in/in)', fontsize=12,fontfamily='Times New Roman',weight='bold')\n",
"plt.ylabel('Stress (ksi)', fontsize=12,fontfamily='Times New Roman',weight='bold')\n",
"plt.title('Concrete Material Stress-Strain', fontsize=14,fontfamily='Times New Roman',weight='bold')\n",
"legend = plt.legend(loc='best', shadow=False, fontsize='8',ncol=3,prop=font)\n",
"\n",
"plt.xticks(rotation=70)\n",
"plt.grid(True)\n",
"plt.show()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## XTRACT Comparison\n",
"\n",
"Now we will compare results"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\r\n",
"\r\n",
"\r\n",
"\r\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"with open('pm.txt', 'r') as f:\n",
" lines = f.readlines()\n",
" y33 = [float(line.split()[0]) for line in lines]\n",
" x33 = [float(line.split()[1]) for line in lines]\n",
"\n",
"\n",
"x=[x11,x22,x33]\n",
"y=[y11,y22,y33]\n",
"\n",
"plt.plot(x[0] ,y[0],color='red',label='ReinforcingSteel')\n",
"plt.plot(x[1] ,y[1],color='blue',label='Steel02')\n",
"plt.plot(x[2] ,y[2],linestyle='dashed',color='green',label='XTRACT')\n",
"plt.title('Moment Curvature XTRACT Comparison', fontsize=14,fontfamily='Times New Roman',weight='bold')\n",
"plt.xlabel('Curvature (1/in)', fontsize=12,fontfamily='Times New Roman',weight='bold')\n",
"plt.ylabel('Moment (kip-in)', fontsize=12,fontfamily='Times New Roman',weight='bold')\n",
"plt.xticks(rotation=70)\n",
"plt.grid(True)\n",
"legend = plt.legend(loc='best', shadow=False, fontsize='8',ncol=3,prop=font)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Conclusions\n",
"\n",
" * The Steel02 steel material can replicate the results of XTRACT accurately because both use the same material stress-strain curve.\n",
" * OpenSees it is faster than XTRACT but the lack of graphical interface it is a big disavantage, also Opensees needs to be checked for maximum concrete strain otherwise it can give unaccurate results. This check needs to be programed manually.\n",
" * OpenSees is a good alternative as long as the results are validated using XTRACT and maximum strains checked."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}