Example Shell Element simply supported Beam

[Prafullamalla]

Code: Select all

# -*- coding: utf-8 -*-
"""
Created on Sat Dec 17 10:01:11 2022

@author: praf_mall@hotmail.com
"""

import openseespy.opensees as ops
import numpy as np
import matplotlib.pyplot as plt

ops.wipe()
ops.model('basic', '-ndm', 2, '-ndf', 2)


L=2000;h=250;b=250

nEleX=20;nEleY=4
delX=L/nEleX;delY=h/nEleY
EndNode=(nEleX+1)*(nEleY+1)

################################################
## Create Nodes for Beam
j = 0
while j < nEleY+1:
    i=0
    while i < nEleX+1:
        nodeTag= j*(nEleX+1)+i+1
        Xcoord=(i*delX)
        Ycoord=(j*delY)
        ops.node(nodeTag,Xcoord,Ycoord)
        #print(nodeTag,Xcoord,Ycoord )
        i = i+1
    j = j+1

# --------------------------------------------------------
# Define shell elements
# --------------------------------------------------------
# Create quad elements - command:
E=36000
matTag=1    
# element quad eleID node1 node2 node3 node4 thick type matID
ops.nDMaterial ('ElasticIsotropic', matTag, E, 0.2)
j= 0;
while j < nEleY:   
    i = 0
    while i < nEleX:
        EleTag= j*nEleX+i+1
        nodeI= j*(nEleX+1)+i+1
        nodeJ= j*(nEleX+1)+i+2
        nodeK=(j+1)*(nEleX+1)+i+2
        nodeL=(j+1)*(nEleX+1)+i+1
        ops.element('quad',EleTag,nodeI,nodeJ,nodeK,nodeL,b,'PlaneStrain', matTag)
       # print(EleTag,nodeI,nodeJ,nodeK,nodeL) 
        i=i+1       
    j = j+1
    
ops.fix(1,1,1)
ops.fix(nEleX+1,0,1)

import opsvis as opsv
# Display Model
opsv.plot_model() 

controlNode=int((nEleX+1)*(nEleY)+nEleX/2+1)
controlDOF=2
Py=-1000
# Create a Plain load pattern with a Linear TimeSeries
ops.timeSeries ('Linear', 1)
ops.pattern    ('Plain', 1, 1)

ops.load(controlNode, 0.0, Py)

Tol = 1.0e-6
ops.system      ('BandGeneral')
ops.constraints ('Plain')
ops.numberer    ('RCM')
ops.algorithm   ('Newton')
ops.integrator  ('LoadControl', 1)
ops.analysis    ('Static')
ops.analyze     (1)   

print('Deflection calculated',ops.nodeDisp(controlNode))
I=(b*h**3)/12
Deflection=Py*L**3/(48*E*I)
print('Deflection',Deflection)
ops.reactions()
print('Reaction 1',ops.nodeReaction(1,2))
print('Reaction 2',ops.nodeReaction(nEleX+1,2))
print(ops.nodeDisp(1,2))
print(ops.nodeDisp(nEleX+1,2))

opsv.plot_loads_2d()
# - plot deformation
opsv.plot_defo(unDefoFlag=1)
plt.axis('equal')

# get values at OpenSees nodes
sig_out = opsv.sig_out_per_node()

# !!! select from sig_out: e.g. vmises
#j, jstr = 0, 'sxx'
j, jstr = 1, 'syy'
#j, jstr = 2, 'sxy'
#j, jstr = 3, 'vmis'
# j, jstr = 4, 's1'
# j, jstr = 5, 's2'
#j, jstr = 6, 'alfa'

nds_val = sig_out[:, j]

plt.figure()
opsv.plot_stress_2d(nds_val)
plt.show()