Este cuaderno realiza el análisis estático no lineal de una estructura de concreto armado compuesta por pórticos y muros de concreto armado. 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 "Análisis Sísmico de Edificaciones en OpenSees".

This notebook perform the nonlinear static analysis of a RC building with frames and walls. This code import libraries from OpenSeesPy. You can run this notebook in the next colab link.

Importamos funciones y librerías¶

Importing functions and libraries¶

pip install httpimport

Requirement already satisfied: httpimport in /usr/local/lib/python3.7/dist-packages (0.5.18)

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.3.0.1)

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: pillow>=6.2.0 in /usr/local/lib/python3.7/dist-packages (from matplotlib==3.3.3) (7.1.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: python-dateutil>=2.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib==3.3.3) (2.8.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: cycler>=0.10 in /usr/local/lib/python3.7/dist-packages (from matplotlib==3.3.3) (0.11.0)
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

import matplotlib
print(matplotlib.__version__)

3.3.3

Importamos librerías modificadas ligeramente de Ops_vis. Además, importamos variables relacionadas a las unidades y parámetros de los materiales.
Importing Libraries slightly modified from Ops_vis. Likewise, units and material parameters are imported.

from openseespy.opensees import *
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import httpimport 
with httpimport.remote_repo(['ops_vis2','Get_Rendering', 'unidades', 'parametros'],'https://jpi-ingenieria.com/ASEP'): 
    import ops_vis2 as opsv2
    import Get_Rendering as opsplt
    from unidades import *
    from parametros import *

%matplotlib inline

Inicializamos el modelo¶

Starting the model¶

Se define que el modelo será un modelo en 3 dimensiones y que considerarán 6 GDL por cada nodo, además se define que se considerarán diafragmas rígidos.

wipe()
model('basic', '-ndm', 3, '-ndf', 6)

Definimos los materiales¶

Material properties¶

Consideramos un modelo Elástico para el concreto que trabajará a cortante y modelos no lineales para el concreto confinado, el concreto no confinado y el acero de refuerzo longitudinal

Programa para calibrar:http://www.jdcui.com/?p=2571

# Concreto Cortante
uniaxialMaterial('Elastic', 3, G)  
# Concreto confinado          tag  f'c   ec0     f'cu   ecu
uniaxialMaterial('Concrete02', 4, fpc1, epsc01, fpcu1, epsU1, lambda1, ft1, Ets)
# Concreto no confinado
uniaxialMaterial('Concrete02', 5, fpc2, epsc02, fpcu2, epsU2, lambda1, ft2, Ets)
# Acero de refuerzo       tag  fy  E0  b
uniaxialMaterial('Steel02', 6, fy, Es, 0.01, 18, 0.925, 0.15)

Definimos los tipos de secciones¶

Creating RC sections¶

Definimos las secciones en base al modelo de fibras, considerando el concreto confinado para la zona ubicada dentro de los estribos de confinamiento y el concreto no confinado para la parte restante.

l1, l2 = 0.4, 0.4

# Sección tipo fibra para Columna
cuant, nc_bars = 0.01, 8
As = l1*l2*cuant/nc_bars
Ac, ρlc, Izc, Iyc, k, Jc = prop_col(l1, l2) # Función que devuelve propiedades de la columna
y1 = l1/ 2.0
z1 = l2/ 2.0

n_y = int((y1)/cover*2)
n_z = int((z1)/cover*2)

fib_sec_c = [['section', 'Fiber', 1, '-GJ', G*Jc],
              ['patch', 'rect', 4, n_y, n_z, cover - y1, cover - z1, y1 - cover, z1 - cover], # Patch 1
              ['patch', 'rect', 5, n_y+2, 1, -y1, z1 - cover, y1, z1], # Patch 2
              ['patch', 'rect', 5, n_y+2, 1, -y1, -z1, y1, cover - z1], # Patch 3
              ['patch', 'rect', 5, 1, n_z, -y1, cover - z1, cover - y1, z1 - cover], # Patch 4
              ['patch', 'rect', 5, 1, n_z, y1 -cover, cover - z1, y1, z1 - cover], # Patch 5
              ['layer', 'straight', 6, 3, As, y1 - cover,
              z1 - cover, y1 - cover, cover - z1],
              ['layer', 'straight', 6, 2, As, 0.0, z1 - cover, 0.0, cover - z1],
              ['layer', 'straight', 6, 3, As, cover - y1, z1 - cover, cover - y1, cover - z1]]

# Ploteo de sección
matcolor = ['r', 'lightgrey', 'gold', 'r', 'lightgrey', 'gold']
opsv2.plot_fiber_section(fib_sec_c, matcolor=matcolor)
plt.axis('equal')
plt.show()
plt.close()

# Crear sección tipo fibra en OpenSees
for li in fib_sec_c:
    
    if li[0] == 'section':
        # section('Fiber',1,'-GJ', G*Jc)
        eval('%s("%s",%s,"%s",%s)' % tuple(li))
    else:
        eval('%s("%s",%s,%s,%s,%s,%s,%s,%s)' % tuple(li))

b,h = 0.2, 0.4

# Sección tipo fibra para Viga
cuant, nc_bars = 0.01, 8
As = b*h*cuant/nc_bars
Av, ρlv, Izv, Iyv, k, Jv = prop_vig(b, h) # Función que devuelve propiedades de la viga
y1 = h / 2.0
z1 = b / 2.0

n_y = int((y1)/cover*2)  # Número de divisiones en la dirección Y
n_z = int((z1)/cover*2)  # Número de divisiones en la dirección Z

fib_sec_v = [['section', 'Fiber', 2, '-GJ', G*Jv],
              ['patch', 'rect', 4, n_y, n_z, cover - y1, cover - z1, y1 - cover, z1 - cover],
              ['patch', 'rect', 5, n_y+2, 1, -y1, z1 - cover, y1, z1],
              ['patch', 'rect', 5, n_y+2, 1, -y1, -z1, y1, cover - z1],
              ['patch', 'rect', 5, 1, n_z, -y1, cover - z1, cover - y1, z1 - cover],
              ['patch', 'rect', 5, 1, n_z, y1 -cover, cover - z1, y1, z1 - cover],
              ['layer', 'straight', 6, 3, As, y1 - cover,
              z1 - cover, y1 - cover, cover - z1],
              ['layer', 'straight', 6, 2, As, 0.0, z1 - cover, 0.0, cover - z1],
              ['layer', 'straight', 6, 3, As, cover - y1, z1 - cover, cover - y1, cover - z1]]

# Ploteo de sección
matcolor = ['r', 'lightgrey', 'gold', 'r', 'lightgrey', 'gold']
opsv2.plot_fiber_section(fib_sec_v, matcolor=matcolor)
plt.axis('equal')
plt.show()
plt.close()

# Crear sección tipo fibra en OpenSees
for li in fib_sec_v:
    if li[0] == 'section':
        eval('%s("%s",%s,"%s",%s)' % tuple(li))
    else:
        eval('%s("%s",%s,%s,%s,%s,%s,%s,%s)' % tuple(li))

Las funciones mx y my calculan los parámetros necesarios para definir los elementos MVLEM_3D.

Los parámetros de entrada son:
- L: Longitud del muro
- t: Espesor del muro
Los parámetros de salidad son:
- muf : Cantidad de elementos verticales
- tt : Espesor de cada elemento vertical
- ww : Ancho de cada elemento vertical
- pp : Cuantía de cada elemento vertical
- conc : Etiqueta del material concreto de cada elemento vertical
- acero : Etiqueta del material acero de cada elemento vertical

def funtion_mx(L, t):
    Lmx = L*m   # Longitud del muro
    t = t*m  # Espesor del muro
    ancho = 50*cm  # Discretización del muro
    mufx = int(round(Lmx/ancho))
    ttx = np.zeros(mufx) 
    ttx[:] = t # Arreglo de espesores
    wwx = np.zeros(mufx)
    wwx[:] = Lmx/(mufx) # Arreglo de anchos
    ρρx = np.zeros(mufx)
    ρρx[:] = 0.0064  # Cuantía vertical en muros
    ρρx[0], ρρx[-1] = 0.01, 0.01  # Cuantía vertical en núcleo
    concx = np.zeros(mufx)
    concx[:] = 5 # Concreto sin confinar
    concx[0], concx[-1] = 4, 4 # Concreto confinado
    acerox = np.zeros(mufx) 
    acerox[:] = 6 # Modelo No Lineal del acero
    return mufx, ttx, wwx, ρρx, concx, acerox

Se crean los nodos con la información importada

# Creamos los nodos
Lx = 1.0
node(1, *[0.,0.,0.])
node(2, *[3.,0.,0.])
node(3, *[3.+Lx,0.,0.])
node(4, *[0.,0.,3.])
node(5, *[3.,0.,3.])
node(6, *[3.+Lx,0.,3.])
node(7, *[0.,0.,6.])
node(8, *[3.,0.,6.])
node(9, *[3.+Lx,0.,6.])

Aplicamos las restricciones¶

Constraint assignment¶

Aplicamos las restricciones de los 6 GDL en la base

fixY(0.0, *[0, 1, 0, 0, 0, 0], '-tol', 1e-6)
fixZ(0.0, *[1, 1, 1, 1, 1, 1], '-tol', 1e-6)

Establecemos la transformación geométrica para los ejes locales

# Establecemos transformación geométrica
geomTransf('PDelta', int(1), *[1, 0, 0]) # Columna
geomTransf('Linear', int(2), *[1, -1, 0]) # Vigas

Se define el método de integración y la cantidad de puntos para cada sección (vigas y columnas)

# Lobatto integration
ni = 5  # Número de puntos de integración a lo largo del elemento
# beamIntegration('Lobatto', tag, secTag, N)
beamIntegration('Lobatto', 1, 1, ni)
    
# beamIntegration('Lobatto', tag, secTag, N)
beamIntegration('Lobatto', 2, 2, ni)

Definimos los elementos con sus respectivas propiedades¶

Creation of RC Elements¶

# Creamos los elementos

Ac = l1*l2
element('forceBeamColumn', 1, 1, 4, 1, 1, '-mass', 2400*Ac*m**2 )
element('forceBeamColumn', 2, 4, 7, 1, 1, '-mass', 2400*Ac*m**2 )

Av = b*h
element('forceBeamColumn', 3, 4, 5, 2, 2, '-mass', 2400*Av*m**2)
element('forceBeamColumn', 4, 7, 8, 2, 2, '-mass', 2400*Av*m**2)

muf, tt, ww, ρρ, conc, acero = funtion_mx(Lx, 0.2)
element('MVLEM_3D', 5, 3, 2, 5, 6, int(muf), '-thick', *tt[:], '-width', *ww[:],'-rho', *ρρ[:], 
        '-matConcrete', *conc[:], '-matSteel', *acero[:], '-matShear', int(3), '-Poisson', 0.2, '-Density', ρ)
element('MVLEM_3D', 6, 6, 5, 8, 9, int(muf), '-thick', *tt[:], '-width', *ww[:],'-rho', *ρρ[:], 
        '-matConcrete', *conc[:], '-matSteel', *acero[:], '-matShear', int(3), '-Poisson', 0.2, '-Density', ρ)

opsplt.plot_model('nodes','elements')

No Model_ODB specified, trying to get data from the active model.
3D model

(<Figure size 432x288 with 1 Axes>, <Axes3DSubplot:xlabel='X', ylabel='Y'>)

Independientemente de las cargas que se tendrán que asignar a la estructura para realizar el análisis por gravedad, tenemos que asignar la masa que tendrá la estructura para poder realizar el análisis modal de la estructura.

Mass assignment and show modes

# ASIGNACIÓN DE MASAS Y MODOS DE VIBRACIÓN
aPlanta = (3.+Lx)*5.
Carga = wTotal*aPlanta*m**2
mass(4, Carga/3, Carga/3, 0.0)
mass(5, Carga/3, Carga/3, 0.0)
mass(6, Carga/3, Carga/3, 0.0)
mass(7, Carga/3, Carga/3, 0.0)
mass(8, Carga/3, Carga/3, 0.0)
mass(9, Carga/3, Carga/3, 0.0)

Nmodes = 4

vals = eigen(Nmodes)
Tmodes = np.zeros(len(vals))
for i in range(Nmodes):
    Tmodes[i] = 2*np.pi/vals[i]**0.5
    print("T[%i]: %.5f" % (i+1, Tmodes[i]))

opsplt.plot_modeshape(1, 100)

T[1]: 0.22656
T[2]: 0.04200
T[3]: 0.01596
T[4]: 0.01552
No Model_ODB specified to plot modeshapes
3D model

(<Figure size 432x288 with 1 Axes>, <Axes3DSubplot:xlabel='X', ylabel='Y'>)

Análisis por cargas de gravedad¶

Gravity Analysis¶

Asignación de cargas¶

Load Assignment¶

En el siguiente fragmento de código se crean los patrones de carga, además se asigna el peso propio de los elementos vigas y columnas a través de cargas distribuidas y el peso propio de los muros en forma puntual en cada vértice de este (25% del peso del muro en cada vértice) en la dirección de la gravedad.

timeSeries('Linear', 1)  # tag
pattern('Plain', 1, 1)  # tag, timeSeries_tag

# Peso propio de los elementos
eleLoad('-ele', 1, '-type','-beamUniform', 0.0, 0.0, -ρ*Ac*g)
eleLoad('-ele', 2, '-type','-beamUniform', 0.0, 0.0, -ρ*Ac*g)

eleLoad('-ele', 3, '-type','-beamUniform', -ρ*Av*g, 0.0, 0.0)
eleLoad('-ele', 4, '-type','-beamUniform', -ρ*Av*g, 0.0, 0.0)

dz = 3.
t = 0.2
load(3, 0., 0., -0.25*dz*Lx*t*ρ*g, 0., 0., 0.)
load(2, 0., 0., -0.25*dz*Lx*t*ρ*g, 0., 0., 0.)
load(5, 0., 0., -0.25*dz*Lx*t*ρ*g, 0., 0., 0.)
load(6, 0., 0., -0.25*dz*Lx*t*ρ*g, 0., 0., 0.)

load(6, 0., 0., -0.25*dz*Lx*t*ρ*g, 0., 0., 0.)
load(5, 0., 0., -0.25*dz*Lx*t*ρ*g, 0., 0., 0.)
load(8, 0., 0., -0.25*dz*Lx*t*ρ*g, 0., 0., 0.)
load(9, 0., 0., -0.25*dz*Lx*t*ρ*g, 0., 0., 0.)

Adicionalmente se tiene que agregar el peso propio de la tabiquería, aligerado, acabados además de la Carga Viva (dependiente del tipo de estructura) directamente a los nodos.

# Carga viva y muerta sobre los nodos

load(4, 0.0, 0.0, -Carga*g/3, 0.0, 0.0, 0.0)
load(5, 0.0, 0.0, -Carga*g/3, 0.0, 0.0, 0.0)
load(6, 0.0, 0.0, -Carga*g/3, 0.0, 0.0, 0.0)
load(7, 0.0, 0.0, -Carga*g/3, 0.0, 0.0, 0.0)
load(8, 0.0, 0.0, -Carga*g/3, 0.0, 0.0, 0.0)
load(9, 0.0, 0.0, -Carga*g/3, 0.0, 0.0, 0.0)

Se establecen los parámetros necesarios para realizar el análisis por cargas de gravedad.

# Análisis de cargas de gravedad
system('UmfPack')
constraints('Transformation')
numberer('RCM')
test('NormDispIncr', 1.0e-12, 10, 3)
algorithm('Newton')
integrator('LoadControl', 0.1)
analysis('Static')
analyze(10)

0

Ploteamos la estructura deformada por cargas de gravedad

opsv2.plot_defo(1000)
plt.show()

vals = eigen(Nmodes)
Tmodes = np.zeros(len(vals))
for i in range(Nmodes):
    Tmodes[i] = 2*np.pi/vals[i]**0.5
    print("T[%i]: %.5f" % (i+1, Tmodes[i]))

T[1]: 0.22756
T[2]: 0.04217
T[3]: 0.01597
T[4]: 0.01553

Análisis Estático No Lineal (Pushover)¶

Nonlinear Static analysis¶

Enlistamos algunos nodos de la base¶

Listing some nodes and elements¶

base  = [1,2,3]
piso1 = [4,5,6]
col1 = [1]
wall1 = [5]
col2 = [2]
wall2 = [6]

Inicializamos el análisis¶

wipeAnalysis()
loadConst('-time', 0.0)

timeSeries('Linear', 2)
pattern('Plain', 2, 2)

print(nodeEigenvector(5, 1)[0])
print(nodeEigenvector(8, 1)[0])
r = nodeEigenvector(5, 1)[0]/nodeEigenvector(8, 1)[0]
print(r)

FF = 10
# load(1, FF*tonf, 0., 0., 0., 0., 0.)
# load(2, FF*tonf, 0., 0., 0., 0., 0.)
# load(3, FF*tonf, 0., 0., 0., 0., 0.)
load(4, FF*tonf*r, 0., 0., 0., 0., 0.)
load(5, FF*tonf*r, 0., 0., 0., 0., 0.)
load(6, FF*tonf*r, 0., 0., 0., 0., 0.)
load(7, FF*tonf, 0., 0., 0., 0., 0.)
load(8, FF*tonf, 0., 0., 0., 0., 0.)
load(9, FF*tonf, 0., 0., 0., 0., 0.)

# eleLoad('-ele', 1, '-type', '-beamPoint', 0., FF*tonf/10, 0.5, 0.0)
# eleLoad('-ele', 2, '-type', '-beamPoint', 0., FF*tonf/10, 0.5, 0.0)

steps = 120
system("UmfPack")
constraints("Transformation")
numberer("RCM")
integrator("DisplacementControl", 8, 1, 0.001)
# integrator('LoadControl', 1./steps)
# integrator('ArcLength',0.0006,0.01)
test('NormDispIncr', 1.0e-8, 100)
algorithm("RaphsonNewton")
analysis("Static")

############################## PUSHOVER ##################################
x = np.zeros((steps+1,2))
y = np.zeros((steps+1,2))

tests = {1:'NormDispIncr', 2: 'RelativeEnergyIncr', 4: 'RelativeNormUnbalance',5: 'RelativeNormDispIncr', 6: 'NormUnbalance'}
algts = {1:'KrylovNewton', 2: 'SecantNewton' , 4: 'RaphsonNewton',5: 'PeriodicNewton', 6: 'BFGS', 7: 'Broyden', 8: 'NewtonLineSearch'}

for k in range(1,steps+1):
  ok = analyze(1)
  if ok!=0:
    for i in tests:
      test(tests[i],1.0e-8, 1000)
      for j in algts:
        if j < 4:
          algorithm(algts[j],'-initial')
        else:
          algorithm(algts[j])
        print('\n\nIntentando con:',tests[i],algts[j])
        ok = analyze(1)
        if ok==0: break
      if ok==0: break
  #
  if ok!=0:
    print('\n'+'='*80)
    print('\n\n\nNo se encontró solución!!! Paso: %i\n\n'%(k))
    break

  V1 = 0
  for bi in col1:
    V1 = V1 - eleForce(bi)[0]
  for bi in wall1:
    V1 = V1 - eleForce(bi)[0]
    V1 = V1 - eleForce(bi)[6]

  V2 = 0
  for bi in col2:
    V2 = V2 - eleForce(bi)[0]
  for bi in wall2:
    V2 = V2 - eleForce(bi)[0]
    V2 = V2 - eleForce(bi)[6]

  x[k]=[nodeDisp(5, 1), nodeDisp(8, 1)-nodeDisp(5, 1)]
  y[k]=[V1, V2]

  if k % int(steps/5) == 0:
    print("Paso:",k)
    opsv2.plot_defo(10)
    plt.show()

0.002898889483759376
0.007838857491450418
0.369810203453896
Paso: 24

Paso: 48

Paso: 72

Paso: 96

WARNING: CTestNormDispIncr::test() - failed to converge 
after: 100 iterations  current Norm: 1.46386e-05 (max: 1e-08, Norm deltaR: 3259.88)
AcceleratedNewton::solveCurrentStep() -The ConvergenceTest object failed in test()
StaticAnalysis::analyze() - the Algorithm failed at iteration: 0 with domain at load factor 0.41541
OpenSees > analyze failed, returned: -3 error flag

Intentando con: NormDispIncr KrylovNewton
Paso: 120

Ploteamos Cortantes de entrepiso y Distorsiones¶

Plotting storey shear vs drift¶

plt.figure(figsize=(15,15))
plt.plot(1000*x[:k,:]/dz, y[:k,:]/kN,':')
plt.axis([-25,25,-170,170])
plt.xlabel('Distorsión (‰)')
plt.ylabel('Cortante de Entrepiso (kN)')
plt.show()