Este cuaderno realiza el análisis dinámico 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 dynamic 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.
pip install httpimport
pip install openseespy
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
import matplotlib
print(matplotlib.__version__)
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
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 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 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)
# 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')
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)
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)
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]))
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)
dt = 0.02
accX = np.genfromtxt('https://jpi-ingenieria.com/ASEP/1974oct_X.txt')
steps = len(accX)
tFinal = steps*dt # 32.
timeSeries('Path', 2, '-dt', dt, '-values', *accX, '-factor', 2.0*cm)
pattern('UniformExcitation', 2, 1, '-accel', 2)
# damping ratio
damp = 0.04
# # lower frequency
# omega1 = 2 * np.pi / Tmodes[0] # * 0.2
# # upper frequency
# omega2 = 2 * np.pi / Tmodes[3] # * 20
# a0 = 2*damp*omega1*omega2/(omega1 + omega2)
# a1 = 2*damp/(omega1 + omega2)
# print("a0,a1:",a0,a1)
# rayleigh(a0, a1, 0.0, 0.0)
power = 2*np.pi/Tmodes[0]
betaKcomm = 2 * (damp/power)
rayleigh(0.0, 0.0, 0.0, betaKcomm)
constraints('Transformation')
numberer('RCM')
system('UmfPack')
test('NormDispIncr', 1e-8, 100)
algorithm('RaphsonNewton')
integrator('Newmark', 0.5, 0.25)
analysis('Transient')
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'}
dtA = 0.01
N = int(steps*dt/dtA)
x = np.zeros((N+1,2))
y = np.zeros((N+1,2))
tCurrent = getTime()
k = 0
while tCurrent < tFinal:
ok = analyze(1, 0.01)
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, 0.01)
if ok==0: break
if ok==0: break
#
if ok!=0:
print('\n'+'='*80)
print('\n\n\nNo se encontró solución!!! Paso: %i, Tiempo: %.2fs\n\n'%(k,tCurrent))
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]
k = k + 1
tCurrent = getTime()
if k % int(N/50) == 0:
print("Paso: %i, Tiempo: %.2fs"%(k,tCurrent))
# opsv2.plot_defo(10)
# plt.show()
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()
Periodos de Vibración después del Análisis Sísmico
Modes after Earthquake Analysis
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]))