Este cuaderno obtiene los modos de vibración de una viga modelada con 4 tipo de elementos: Elementos Viga Euler-Bernoulli, Elementos Viga Timoshenko, Elementos cuadrados 2D y Elementos cúbicos 3D. Este código usa las librerías de OpenSeesPy. Puede ejecutar este cuaderno a partir de este enlace de Colab.
Cuaderno realizado como material del curso "Programando el Método de Elmentos Finitos".
This notebook obtains mode shapes of a beam using four types of Elements: Euler-Bernoulli Beam Elements, Timoshenko Beam Elements, 2D Quad Elements, and 3D Brick Elements. This code import libraries from OpenSeesPy. You can run this notebook in the next colab link.
Importing Libraries¶
In [1]:
pip install openseespy
In [2]:
pip install matplotlib==3.3.3
Verificando la versión de Matplotlib, debe ser 3.3.3
Checking the version of Matplotlib, must be 3.3.3
In [3]:
import matplotlib
print(matplotlib.__version__)
In [4]:
from openseespy.opensees import *
import openseespy.postprocessing.ops_vis as opsv
import openseespy.postprocessing.Get_Rendering as opsplt
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Definimos unidades
Defining Units
In [5]:
# Unidades Base
m = 1
kg = 1
s = 1
# Otras Unidades
cm = 0.01*m
N = kg*m/s**2
kgf = 9.81*N
tonf = 1000*kgf
Pa = N/m**2
MPa = 10**6*Pa
# Constantes Físicas
g = 9.80665*m/s**2
Definimos propiedades y geometría
Defining properties and geometry
In [6]:
# Properties
E = 2.1*10**5*MPa
u = 0.25
G = 0.5*E/(1+u)
# Beam
a = 5*cm
A = a**2
Avy, Avz = 5/6*A, 5/6*A
Iz = a**4/12
Iy = a**4/12
β= 1/3-0.21*1.*(1-(1.)**4/12)
Jxx = β*a**4
ρ = 7850*kg/m**3
# Geometry
Lx = 20*cm # Change this length for short or long beams.
# Modes
Nmodes = 5
Euler-Bernoulli Beam Elements¶
In [7]:
wipe()
model('basic', '-ndm', 3, '-ndf', 6)
Ne = 10
# Creamos los nodos para Ne Elementos
# (Creating nodes for Ne Elements)
for i in range(Ne+1):
node(i+1, *[i*Lx/Ne,0.,0.])
# Restricciones
fix(1,*[1,1,1,0,0,0])
fix(Ne+1,*[0,1,1,0,0,0])
fixY(0.,*[0,1,0,0,0,0])
geomTransf('Linear', 1, *[1, -1, 0])
# Creamos Ne Elementos Viga Euler-Bernoulli
# (Creating Ne Euler-Bernoulli Elements)
for i in range(Ne):
element('elasticBeamColumn', i+1, i+1, i+2, A, E, G, Jxx, Iy, Iz, 1,'-mass', ρ*A)
In [8]:
opsv.plot_model(fig_wi_he=(40., 40.),az_el=(-90,0))
Obtenemos las frecuencias de los primeros Nmodes
Obtaining the first Nmodes frequencies
In [9]:
vals = eigen('-fullGenLapack',Nmodes)
Fmodes = np.zeros(len(vals))
for i in range(Nmodes):
Fmodes[i] = vals[i]**0.5/(2*np.pi)
print("F[%i]: %.5f" % (i+1, Fmodes[i]))
In [10]:
for i in range(1,Nmodes+1):
opsv.plot_mode_shape(i, 1)
plt.show()
Timoshenko Beam Elements¶
In [11]:
wipe()
model('basic', '-ndm', 3, '-ndf', 6)
# Creamos los nodos para Ne Elementos
# (Creating nodes for Ne Elements)
Ne = 20
for i in range(Ne+1):
node(i+1, *[i*Lx/Ne,0.,0.])
# Restricciones
fix(1,*[1,1,1,0,0,0])
fix(Ne+1,*[0,1,1,0,0,0])
fixY(0.,*[0,1,0,0,0,0])
geomTransf('Linear', 1, *[1, -1, 0])
# Creamos Ne Elementos viga Timoshenko
# (Creating Ne Timoshenko Elements)
for i in range(Ne):
element('ElasticTimoshenkoBeam', i+1, i+1, i+2, E, G, A, Jxx, Iy, Iz, Avy, Avz, 1,'-mass', ρ*A)
In [12]:
opsv.plot_model(fig_wi_he=(40., 40.),az_el=(-90,0))
Obtenemos las frecuencias de los primeros Nmodes
Obtaining the first Nmodes frequencies
In [13]:
vals = eigen('-fullGenLapack',Nmodes)
Fmodes = np.zeros(len(vals))
for i in range(Nmodes):
Fmodes[i] = vals[i]**0.5/(2*np.pi)
print("F[%i]: %.5f" % (i+1, Fmodes[i]))
In [14]:
for i in range(1,Nmodes+1):
opsplt.plot_modeshape(i, 0.02)
2D Quad Elements¶
In [15]:
wipe()
model('basic', '-ndm', 2, '-ndf', 2)
# Definimos un material Elástico e Isotrópico
nDMaterial('ElasticIsotropic', 1, E, u, ρ)
# Configuramos el tamaño de malla
# Setting the mesh size
ms = a/8 # try: ms = a/2,a/4,a/6,a/8
nx = int(Lx/ms)
ny = int(a/ms)
# Creamos los Nodos
# Creating Nodes
ni = 1
for i in range(nx+1):
for j in range(ny+1):
node(ni,*[ms*i,ms*j])
ni = ni + 1
# Creamos los Elementos
# Creating Elements
ne = 1
for i in range(nx):
for j in range(ny):
element('quad', ne, *[i*(ny+1)+j+1,i*(ny+1)+j+2,(i+1)*(ny+1)+j+2,(i+1)*(ny+1)+j+1], a,'PlaneStrain', 1)
ne = ne + 1
# Asignamos restricciones en (0,a/2) y (Lx,a/2)
# Fixing at (0,a/2) and (Lx,a/2)
for i in range(1,ni):
if round(nodeCoord(i)[0],5)== 0.0 and round(nodeCoord(i)[1],5) == a/2:
fix(i,*[1,1])
if round(nodeCoord(i)[0],5) == Lx and round(nodeCoord(i)[1],5) == a/2:
fix(i,*[0,1])
In [16]:
opsv.plot_model(node_labels=0, element_labels=0,fig_wi_he=(50,8))
Obtenemos las frecuencias de los primeros Nmodes
Obtaining the first Nmodes frequencies
In [17]:
vals = eigen('-fullGenLapack',Nmodes)#,
Fmodes = np.zeros(len(vals))
for i in range(Nmodes):
Fmodes[i] = vals[i]**0.5/(2*np.pi)
print("F[%i]: %.5f" % (i+1, Fmodes[i]))
In [18]:
for i in range(Nmodes):
opsv.plot_mode_shape(i+1,0.01)
plt.show()
3D Brick Elements¶
In [19]:
wipe()
model('basic', '-ndm', 3, '-ndf', 3)
# Definimos un material Elástico e Isotrópico
nDMaterial('ElasticIsotropic', 1, E, u, ρ)
# Configuramos el tamaño de malla
# Setting the mesh size
ms = a/4 # try: ms = a/2,a/4,a/6,a/8
nx, ny, nz = int(Lx/ms), int(a/ms), int(a/ms)
# Creamos los Nodos
# Creating Nodes
ni = 1
for i in range(nx+1):
for j in range(ny+1):
for k in range(nz+1):
node(ni,i*ms,j*ms,k*ms)
fix(ni,*[0,1,0,0,0,0])
ni = ni + 1
# Creamos los Eleemntos
# Creating Elements
ne = 1
for k in range(nx):
for j in range(ny):
for i in range(nz):
nodos = [(nz+1)*j+i+1+k*(ny+1)*(nz+1), (nz+1)*j+i+2+k*(ny+1)*(nz+1), (nz+1)*(j+1)+i+2+k*(ny+1)*(nz+1), (nz+1)*(j+1)+i+1+k*(ny+1)*(nz+1),
(nz+1)*j+i+1+(k+1)*(ny+1)*(nz+1), (nz+1)*j+i+2+(k+1)*(ny+1)*(nz+1), (nz+1)*(j+1)+i+2+(k+1)*(ny+1)*(nz+1),(nz+1)*(j+1)+i+1+(k+1)*(ny+1)*(nz+1)]
element('stdBrick', ne, *nodos, 1, 0, 0, 0.0)
ne = ne + 1
# Asignamos restricciones en (0,a/2) y (Lx,a/2)
# Fixing at (0,a/2) and (Lx,a/2)
for i in range(1,ni):
if nodeCoord(i)[0] == 0 and nodeCoord(i)[2] == a/2:
fix(i,*[1,0,1,0,0,0])
if nodeCoord(i)[0] == Lx and nodeCoord(i)[2] == a/2:
fix(i,*[0,0,1,0,0,0])
In [20]:
opsv.plot_model(node_labels=0, element_labels=0,fig_wi_he=(40.0, 40.0))
Obtenemos las frecuencias de los primeros Nmodes
Obtaining the first Nmodes frequencies
In [21]:
vals = eigen('-fullGenLapack',Nmodes)#'-fullGenLapack',
Fmodes = np.zeros(len(vals))
for i in range(Nmodes):
Fmodes[i] = vals[i]**0.5/(2*np.pi)
print("F[%i]: %.5f" % (i+1, Fmodes[i]))
In [22]:
for i in range(Nmodes):
opsv.plot_mode_shape(i+1,0.05)
plt.show()
