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
Requirement already satisfied: openseespy in /usr/local/lib/python3.7/dist-packages (3.3.0.1.1)
Requirement already satisfied: openseespylinux>=3.3.0.1 in /usr/local/lib/python3.7/dist-packages (from openseespy) (3.4.0.1)
In [2]:
pip install matplotlib==3.3.3
Requirement already satisfied: matplotlib==3.3.3 in /usr/local/lib/python3.7/dist-packages (3.3.3)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib==3.3.3) (1.3.2)
Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.3 in /usr/local/lib/python3.7/dist-packages (from matplotlib==3.3.3) (3.0.6)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.7/dist-packages (from matplotlib==3.3.3) (7.1.2)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.7/dist-packages (from matplotlib==3.3.3) (0.11.0)
Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib==3.3.3) (2.8.2)
Requirement already satisfied: numpy>=1.15 in /usr/local/lib/python3.7/dist-packages (from matplotlib==3.3.3) (1.19.5)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.7/dist-packages (from python-dateutil>=2.1->matplotlib==3.3.3) (1.15.0)

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__)
3.3.3
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]))
F[1]: 2931.64297
F[2]: 6458.59793
F[3]: 11725.26061
F[4]: 19216.76167
F[5]: 26367.22179
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]))
F[1]: 2728.42051
F[2]: 6463.58112
F[3]: 9207.75283
F[4]: 16982.11690
F[5]: 19350.89322
In [14]:
for i in range(1,Nmodes+1):
  opsplt.plot_modeshape(i, 0.02)
No Model_ODB specified to plot modeshapes
3D model
No Model_ODB specified to plot modeshapes
3D model
No Model_ODB specified to plot modeshapes
3D model
No Model_ODB specified to plot modeshapes
3D model
No Model_ODB specified to plot modeshapes
3D model

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]))
F[1]: 2720.36201
F[2]: 4925.33962
F[3]: 8693.47779
F[4]: 15419.23896
F[5]: 15901.05496
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))