✅RC section mesh#

Here we regard the RC section as a single-material section, ignoring the different material properties of steel bars and concrete.

[1]:
import numpy as np
import openseespy.opensees as ops
import opstool as opst

Generate geometric objects#

Functions: offset() and add_polygon() are used.

[2]:
# the points of the outer contour line, only the turning point of the line is needed, counterclockwise or clockwise.
outlines = [[0.5, 0], [7.5, 0], [8, 0.5], [8, 4.5],
            [7.5, 5], [0.5, 5], [0, 4.5], [0, 0.5]]
# cover thick
cover_d = 0.08
# Offset to get the inner boundary of the cover layer
coverlines = opst.offset(outlines, d=cover_d)

# Generate polygonal geometry object for cover layer
cover = opst.add_polygon(outlines, holes=[coverlines])

# Creating core with voids
holelines1 = [[1, 1], [3.5, 1], [3.5, 4], [1, 4]]
holelines2 = [[4.5, 1], [7, 1], [7, 4], [4.5, 4]]
core = opst.add_polygon(coverlines, holes=[holelines1, holelines2])

Generate mesh#

by class SecMesh.

[3]:
sec = opst.SecMesh(sec_name="My Fiber Section")
# Grouping, the dict key is the group name, which can be arbitrary.
sec.assign_group({"cover": cover, "core":core})
# Specify the mesh size
sec.assign_mesh_size(dict(cover=0.2, core=0.4))
# Specify the region color
sec.assign_group_color(dict(cover="gray", core="green"))
# Specify the material tag in the opensees, the material needs to be defined by you beforehand.
ops.uniaxialMaterial('Concrete01', 1, -30, -0.002, -15, -0.005)
ops.uniaxialMaterial('Concrete01', 2, -40, -0.006, -30, -0.015)
sec.assign_ops_matTag(dict(cover=1, core=2))
# mesh!
sec.mesh()

add rebars#

by class Rebars.

[4]:
ops.uniaxialMaterial('Steel01', 3, 200, 2.E5, 0.02)
# Instantiating the rebar class
rebars = opst.Rebars()

rebar_d_outer = 0.032   # dia of rebar
rebar_d_inner = 0.02
# Offset to obtain the rebars arranged along the contour line, Inward offset is positive
rebar_lines1 = opst.offset(outlines, d=cover_d + rebar_d_outer / 2)
# add the rebar line, gap is the spacing of the rebars, matTag is the opensees material tag predefined.
rebars.add_rebar_line(
    points=rebar_lines1, dia=rebar_d_outer, gap=0.15, color="red", matTag=3
)
# Offset to obtain the rebars arranged along the holes
rebar_lines2 = opst.offset(holelines1, d=-(cover_d + rebar_d_inner / 2))
rebars.add_rebar_line(
    points=rebar_lines2, dia=rebar_d_inner, gap=0.2, color="black", matTag=3
)
rebar_lines3 = opst.offset(holelines2, d=-(cover_d + rebar_d_inner / 2))
rebars.add_rebar_line(
    points=rebar_lines3, dia=rebar_d_inner, gap=0.2, color="black", matTag=3
)
# ---------------------------
# add rebars to the sec
sec.add_rebars(rebars)

Get the section properties#

section props dict by method get_sec_props(), including:

  • Cross-sectional area (A)

  • Shear area (Asy, Asz)

  • Elastic centroid (centroid)

  • Second moments of area about the centroidal axis (Iy, Iz, Iyz)

  • Elastic section moduli about the centroidal axis with respect to the top and bottom fibres (Wyt, Wyb, Wzt, Wzb)

  • Torsion constant (J)

  • Principal axis angle (phi)

  • Section mass (mass), only true if material density is defined, otherwise geometric area (mass density is 1)

  • ratio of reinforcement (rho_rebar)

If you need a faster calculation of section properties, you can use method get_frame_props(), except that the shear area cannot be calculated.

Tip

Since the finite element method is used to calculate the section properties by sectionproperties pacakge, the calculation may be a bit slow, which can be ignored if there is no need to calculate the section properties.

[5]:
sec_props = sec.get_sec_props(display_results=True, plot_centroids=False)
# or more faster
sec_props2 = sec.get_frame_props(display_results=True)
# for key, value in sec_props.items():
#     print(f"{key}: {value}")
                              Section Properties                               
┏━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Symbol     Value                   Definition                             ┃
┡━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│ A          2.450E+01               Cross-sectional area                   │
│ Asy        1.400E+01               Shear area y-axis                      │
│ Asz        1.215E+01               Shear area z-axis                      │
│ centroid   (4.000E+00, 2.500E+00)  Elastic centroid                       │
│ Iy         6.935E+01               Moment of inertia y-axis               │
│ Iz         1.522E+02               Moment of inertia z-axis               │
│ Iyz        -3.126E-13              Product of inertia                     │
│ Wyt        2.774E+01               Section moduli of top fibres y-axis    │
│ Wyb        2.774E+01               Section moduli of bottom fibres y-axis │
│ Wzt        3.806E+01               Section moduli of top fibres z-axis    │
│ Wzb        3.806E+01               Section moduli of bottom fibres z-axis │
│ J          1.584E+02               Torsion constant                       │
│ phi        -9.000E+01              Principal axis angle                   │
│ mass       2.450E+01               Section mass                           │
│ rho_rebar  6.747E-03               Ratio of reinforcement                 │
└───────────┴────────────────────────┴────────────────────────────────────────┘
                           Frame Section Properties                            
┏━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Symbol     Value                   Definition                             ┃
┡━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│ A          2.450E+01               Cross-sectional area                   │
│ centroid   (4.000E+00, 2.500E+00)  Elastic centroid                       │
│ Iy         6.935E+01               Moment of inertia y-axis               │
│ Iz         1.522E+02               Moment of inertia z-axis               │
│ Iyz        -3.126E-13              Product of inertia                     │
│ Wyt        2.774E+01               Section moduli of top fibres y-axis    │
│ Wyb        2.774E+01               Section moduli of bottom fibres y-axis │
│ Wzt        3.806E+01               Section moduli of top fibres z-axis    │
│ Wzb        3.806E+01               Section moduli of bottom fibres z-axis │
│ J          1.584E+02               Torsion constant                       │
│ phi        -9.000E+01              Principal axis angle                   │
│ mass       2.450E+01               Section mass                           │
│ rho_rebar  6.747E-03               Ratio of reinforcement                 │
└───────────┴────────────────────────┴────────────────────────────────────────┘

centering or rotate the section#

[6]:
sec.centring()
#sec.rotate(90)

View the section mesh#

engine=’plotly’

[7]:
sec.view(fill=True, engine='plotly', save_html=None, on_notebook=True)

and engine=’matplotlib’

[8]:
sec.view(fill=True, engine='mpl')
../../_images/src_notebooks_mod_secmesh_22_0.png

plot and get stresses#

plot the stress by method get_stress().

This is just for demonstration, the stresses are only realistic if you specify the correct material for each geometry region.

[9]:
stresses = sec.get_stress(N=-1000, Vy=1000, Vz=1000,
                          Myy=1000, Mzz=1000, Mxx=1000,
                          plot_stress='xx', cmap="jet")
for stress in stresses:
    print('Region: {0}'.format(stress['Region']))
    print('List Size: {0}'.format(len(stress['sig_xx'])))
    print('Normal Stresses: {0}'.format(stress['sig_xx']))
    print('von Mises Stresses: {0}'.format(stress['sig_vm']))
../../_images/src_notebooks_mod_secmesh_25_0.png
Region: cover
List Size: 3080
Normal Stresses: [-53.87153561 -99.85483929 -95.9299888  ...   0.           0.
   0.        ]
von Mises Stresses: [ 55.07925233 108.01423572 105.22433508 ...   0.           0.
   0.        ]
Region: core
List Size: 3080
Normal Stresses: [  0.           0.           0.         ... -51.04937768 -42.12135335
 -42.93883398]
von Mises Stresses: [  0.           0.           0.         ... 185.49767954 176.84526638
 172.62404039]

Generate py or tcl file#

[10]:
G = 10000   # Shear modulus
sec.to_file("mysec.py", secTag=1, GJ=G * sec_props['J'])
# If you have calculated the section properties
# sec.to_file("mysec.tcl", secTag=1, GJ=G * sec_props['J'])

Generate openseespy cmds implicitly#

This command can be used after defining the OpenSees material. This defines the openseespy fiber command directly, without outputting the file.

G = 10000   # Shear modulus
J = sec_props['J']  # or other number if you don't care
sec.opspy_cmds(secTag=1, GJ=G * J)

It’s done, you don’t need to do anything more. Your section mesh has been delivered to the openseespy domain.

✅Composite Section Mesh#

Of course, we can also build composite material fiber sections, we need to use function add_material().

[11]:
import numpy as np
import opstool as opst

Specify the characteristics of each material

[12]:
Ec = 3.45E7
Es = 2.0E8
Nus = 0.3
Nuc = 0.2
pho_c = 2.55
pho_s = 7.86
steel_mat = opst.add_material(name='steel', elastic_modulus=Es, poissons_ratio=Nus, density=pho_s)
conc_mat = opst.add_material(name='conc', elastic_modulus=Ec, poissons_ratio=Nuc, density=pho_c)

Use predefined materials when generating geometric objects

[13]:
outlines = [[0, 0], [2, 0], [2, 2], [0, 2]]
coverlines = opst.offset(outlines, d=0.05)
cover = opst.add_polygon(outlines, holes=[coverlines], material=conc_mat)
bonelines = [[0.5, 0.5], [1.5, 0.5], [1.5, 0.7], [1.1, 0.7], [1.1, 1.3], [1.5, 1.3], [1.5, 1.5],
             [0.5, 1.5], [0.5, 1.3], [0.9, 1.3], [0.9, 0.7], [0.5, 0.7], [0.5, 0.5]]
core = opst.add_polygon(coverlines, holes=[bonelines], material=conc_mat)
bone = opst.add_polygon(bonelines, material=steel_mat)

mesh

[14]:
sec = opst.SecMesh()
sec.assign_group(dict(cover=cover, core=core, bone=bone))
sec.assign_mesh_size(dict(cover=0.02, core=0.05, bone=0.02))
# sec.assign_group_color(dict(cover="gray", core="#b84592", bone='#ffc168'))
sec.assign_ops_matTag(dict(cover=1, core=2, bone=4))
sec.mesh()

add rebars

[15]:
# add rebars
rebars = opst.Rebars()
rebar_lines1 = opst.offset(outlines, d=0.05 + 0.032 / 2)
rebars.add_rebar_line(
    points=rebar_lines1, dia=0.032, gap=0.1, color="red", matTag=3
)
# add to the sec
sec.add_rebars(rebars)

Important

Use the modulus of elasticity and shear to convert the geometric properties of the composite section into the equivalent properties of the reference material. Here, steel bar.

[16]:
Gs = Es / (2 * (1 + Nus))  # shear modulu of steel
sec_props = sec.get_sec_props(Eref=Es, Gref=Gs, display_results=True)
                              Section Properties                               
┏━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Symbol     Value                   Definition                             ┃
┡━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│ A          1.120E+00               Cross-sectional area                   │
│ Asy        9.263E-01               Shear area y-axis                      │
│ Asz        7.878E-01               Shear area z-axis                      │
│ centroid   (1.000E+00, 1.000E+00)  Elastic centroid                       │
│ Iy         2.870E-01               Moment of inertia y-axis               │
│ Iz         2.579E-01               Moment of inertia z-axis               │
│ Iyz        -1.043E-15              Product of inertia                     │
│ Wyt        2.870E-01               Section moduli of top fibres y-axis    │
│ Wyb        2.870E-01               Section moduli of bottom fibres y-axis │
│ Wzt        2.579E-01               Section moduli of top fibres z-axis    │
│ Wzb        2.579E-01               Section moduli of bottom fibres z-axis │
│ J          4.477E-01               Torsion constant                       │
│ phi        0.000E+00               Principal axis angle                   │
│ mass       1.296E+01               Section mass                           │
│ rho_rebar  5.241E-02               Ratio of reinforcement                 │
└───────────┴────────────────────────┴────────────────────────────────────────┘
[17]:
sec.centring()
sec.view(fill=True, engine='plotly', save_html=None, on_notebook=True)

plot the stress by method get_stress().

The stresses are realistic because you specify the correct material for each geometry region.

[18]:
stresses = sec.get_stress(N=-1000, Vy=1000, Vz=1000,
                          Myy=1000, Mzz=1000, Mxx=1000,
                          plot_stress='xyz', cmap="coolwarm")
for stress in stresses:
    print('Region: {0}'.format(stress['Region']))
    print('List Size: {0}'.format(len(stress['sig_xyz'])))
    print('Shear Stresses: {0}'.format(stress['sig_xyz']))
    print('von Mises Stresses: {0}'.format(stress['sig_vm']))
../../_images/src_notebooks_mod_secmesh_47_0.png
Region: cover
List Size: 1525
Shear Stresses: [31.29209065 40.719758   29.9125975  ...  0.          0.
  0.        ]
von Mises Stresses: [ 101.74438855 1425.5058128   227.81625064 ...    0.            0.
    0.        ]
Region: core
List Size: 1525
Shear Stresses: [  0.           0.           0.         ...  82.65107602 426.59627777
 475.53108712]
von Mises Stresses: [  0.           0.           0.         ... 383.25058152 757.80592711
 860.3477186 ]
Region: bone
List Size: 1525
Shear Stresses: [   0.            0.            0.         ...  135.03001902 2010.6504183
    0.        ]
von Mises Stresses: [   0.            0.            0.         ... 2074.15602317 3616.59594814
    0.        ]

output the file

[19]:
sec.to_file("mysec.py", secTag=1, GJ=Gs * sec_props['J'])
# or
# sec.opspy_cmds(secTag=1, GJ=Gs * sec_props['J'])

✅Variable section mesh#

[20]:
import numpy as np
import opstool as opst

Function: var_line_string() is used to generate geometric line points that vary linearly or parabolicly, and then generate a cross-sectional mesh.

First determine the geometric data at both ends.

[21]:
# I end
outlines1 = [[0.5, 0], [7.5, 0], [8, 0.5], [8, 4.5],
            [7.5, 5], [0.5, 5], [0, 4.5], [0, 0.5]]
cover_d = 0.08
coverlines1 = opst.offset(outlines1, d=cover_d)
holelines1i = [[1, 1], [3.5, 1], [3.5, 4], [1, 4]]
holelines2i = [[4.5, 1], [7, 1], [7, 4], [4.5, 4]]
# J end
outlines2 = [[0.5, 0], [7.5, 0], [8, 0.5], [8, 2.5],
             [7.5, 3], [0.5, 3], [0, 2.5], [0, 0.5]]
cover_d = 0.05
coverlines2 = opst.offset(outlines2, d=cover_d)
holelines1j = [[1, 1], [3.5, 1], [3.5, 2], [1, 2]]
holelines2j = [[4.5, 1], [7, 1], [7, 2], [4.5, 2]]

To define a path of a section normal, it is recommended to determine a beam element every two coordinate points, where n_sec can be the number of fiber section integration points for each beam element.

[22]:
path = [(0, 0, 0), (0, 0, 16), (0, 0, 24), (0, 0, 40)]   # 3 beam elements

Get geometric data at each integration point.

[23]:
locs = opst.get_lobatto_loc(n_sec=5)
[24]:

outlines = opst.var_line_string(pointsi=outlines1, pointsj=outlines2, path=path, n_sec=5, loc_sec=locs, closure=True, y_degree=1, y_sym_plane="j-0", z_degree=2, z_sym_plane="j-0") coverlines = opst.var_line_string(pointsi=coverlines1, pointsj=coverlines2, path=path, n_sec=5, loc_sec=locs, closure=True, y_degree=1, y_sym_plane="j-0", z_degree=2, z_sym_plane="j-0") holelines1 = opst.var_line_string(pointsi=holelines1i, pointsj=holelines1j, path=path, n_sec=5, loc_sec=locs, closure=True, y_degree=1, y_sym_plane="j-0", z_degree=2, z_sym_plane="j-0") holelines2 = opst.var_line_string(pointsi=holelines2i, pointsj=holelines2j, path=path, n_sec=5, loc_sec=locs, closure=True, y_degree=1, y_sym_plane="j-0", z_degree=2, z_sym_plane="j-0")

Generate the section mesh at each integration point.

[25]:
sec_meshes = []
for i in range(len(outlines)):
    cover = opst.add_polygon(outlines[i], holes=[coverlines[i]])
    core = opst.add_polygon(coverlines[i], holes=[holelines1[i], holelines2[i]])
    sec = opst.SecMesh(sec_name=f"Section{i+1}")
    sec.assign_group(dict(cover=cover, core=core))
    sec.assign_mesh_size(dict(cover=0.2, core=0.4))
    sec.assign_group_color(dict(cover="gray", core="green"))
    sec.assign_ops_matTag(dict(cover=1, core=2))
    # mesh!
    sec.mesh()
    rebars = opst.Rebars()
    rebar_d_outer = 0.06   # dia of rebar
    rebar_d_inner = 0.06
    rebar_lines1 = opst.offset(coverlines[i], d=rebar_d_outer / 2)
    rebars.add_rebar_line(points=rebar_lines1, dia=rebar_d_outer,
                          gap=0.15, color="red", matTag=3)
    rebar_lines2 = opst.offset(holelines1[i], d=-(cover_d + rebar_d_inner / 2))
    rebars.add_rebar_line(
        points=rebar_lines2, dia=rebar_d_inner, gap=0.2, color="black", matTag=3
    )
    rebar_lines3 = opst.offset(holelines2[i], d=-(cover_d + rebar_d_inner / 2))
    rebars.add_rebar_line(
        points=rebar_lines3, dia=rebar_d_inner, gap=0.2, color="black", matTag=3
    )
    sec.add_rebars(rebars)
    sec.centring()  # must here
    sec_meshes.append(sec)

Function: vis_var_sec() is used to visualize varied section meshes.

[26]:
opst.vis_var_sec(sec_meshes=sec_meshes, n_sec=5, loc_sec=locs,
                 path=path, on_notebook=False)
../../_images/var_sec_mesh.png

You can also generate openseespy commands directly.

secTags = []
for i, sec in enumerate(sec_meshes):
    sec.opspy_cmds(secTag=i+1, GJ=1000000)
    secTags.append(i+1)

Using the variable section integration scheme in openseespy:

n_sec = 5
n_ele = int(len(sec_meshes) / n_sec)
for i in range(n_ele):
    tag = i + 1
    sec_tags = secTags[n_sec*i : n_sec*(i+1)]
    ops.beamIntegration('Lobatto', tag, n_sec, *sec_tags)