# Nonlinear seismic response of a MRF¶

#### March 2020, By Amir Hossein Namadchi¶

This is an OpenSeesPy simulation of a moment resisting frame subjected to seismic excitation. The model was introduced by C. Kolay & J. M. Ricles in their paper entitled Assessment of explicit and semi-explicit classes of model-based algorithms for direct integration in structural dynamics. The beams and columns of the MRF are modeled using dispBeamColumn fiber elements. The gravity load resisting system associated with the MRF is modeled using a lean-on column composed of linear elastic beam-column elements with 2nd Order $P-\Delta$ effects [1].

In [1]:
import numpy as np
import openseespy.opensees as ops
import matplotlib.pyplot as plt
import eSEESminiPy
%matplotlib notebook


## Units¶

In [2]:
## Units
m = 1.0               # Meters
KN = 1.0              # KiloNewtons
sec = 1.0             # Seconds
inch = 0.0254*m       # inches

kg = KN*(sec**2)/m    # mass unit (derived)
g = 9.81*(m/sec**2)   # gravitational constant


## Earthquake record¶

This will load 1994 Northridge earthquake ground motion record (Canyon Country - W Lost Cany station) downloaded for the PEER website. Then, the record is scaled by a factor of 3 as follows (scaling could also be done when defining timeSeries):

In [3]:
dt = 0.01*sec
northridge = np.column_stack((np.arange(0,len(northridge)*dt, dt),
northridge*3*g))


## Model Definition¶

### Geometry¶

Node coordinates and element connectivity are defined.

In [4]:
ops.wipe()
ops.model('basic','-ndm',2,'-ndf',3)

## Main Nodes
# Node Coordinates Matrix (size : nn x 2)
node_coords = np.array([[0,0],[6,0],
[0,1],[6,1],
[0,2],[6,2],
[0,2.5],[6,2.5],
[0,3],[0.5,3],[1,3],[3,3],[5,3],[5.5,3],[6,3],
[0,3.5],[6,3.5],
[0,4],[6,4],
[0,5],[6,5],
[0,5.5],[6,5.5],
[0,6],[0.5,6],[1,6],[3,6],[5,6],[5.5,6],[6,6]
], dtype = np.float64)*m
## Main Elements
# Element Connectivity Matrix (size: nel x 2)
connectivity = [[1, 3], [3, 5], [5, 7], [7, 9],
[9, 10], [10, 11], [11, 12], [12, 13],
[13, 14], [14, 15], [15, 8], [8, 6],
[6, 4], [4, 2], [9, 16], [16, 18],
[18, 20], [20, 22], [22, 24], [24, 25],
[25, 26], [26, 27], [27, 28], [28, 29],
[29, 30], [30, 23], [23, 21], [21, 19],
[19, 17], [17, 15]]

# Get Number of elements
nel = len(connectivity)

# Distinguish beams and columns by their element tag ID
all_the_beams = list(range(5, 10+1)) + list(range(20, 25+1))
all_the_cols = list(np.setdiff1d(np.arange(1, nel+1),
all_the_beams))


### Sections & Material¶

sections are defined in dict which is quite self-explanatory.

In [5]:
# Main Beams and Columns
sections = {'W24x55':{'d':23.57*inch, 'tw':0.395*inch,
'bf':7.005*inch, 'tf':0.505*inch,
'A':16.2*(inch**2),
'I1':1350*(inch**4), 'I2':29.1*(inch**4)},

'W14x120':{'d':14.48*inch, 'tw':0.590*inch,
'bf':14.670*inch, 'tf':0.940*inch,
'A':35.3*(inch**2),
'I1':1380*(inch**4), 'I2':495*(inch**4)}
}

# Leaning columns section properties
leaning_col = {'A':(9.76e-2)*m**2,
'I1':(7.125e-4)*m**4}

# Material properties
F_y = 345000*(KN/m**2)     # yield strength
E_0 = 2e8*(KN/m**2)        # initial elastic tangent
eta = 0.01                 # strain-hardening ratio
rho = 7850*(kg/m**3)       # mass density


In [6]:
# Nodal loads and Masses
lumped_mass = 50.97*kg     # seismic floor mass
P_1 = 500*KN               # Nodal loads

# Adding nodes to the domain
## Main Nodes
[ops.node(n+1,*node_coords[n])
for n in range(len(node_coords))];
## Fictitious Nodes (Leaning columns)
ops.node(100,*[7.0, 0.0])                                # @ Base
ops.node(101,*[7.0, 3.0],
'-mass', *[lumped_mass, 0.00001, 0.00001])      # @ Story 1
ops.node(102,*[7.0, 6.0],
'-mass', *[lumped_mass, 0.00001, 0.00001])      # @ Story 2 (roof)

# Material
# -> uniaxial bilinear steel material with kinematic hardening
ops.uniaxialMaterial('Steel01', 1,
F_y, E_0, eta)

## Beams
ops.section('WFSection2d', 1, 1,
sections['W24x55']['d'],
sections['W24x55']['tw'],
sections['W24x55']['bf'],
sections['W24x55']['tf'], 10, 3)
## Columns
ops.section('WFSection2d', 2, 1,
sections['W14x120']['d'],
sections['W14x120']['tw'],
sections['W14x120']['bf'],
sections['W14x120']['tf'], 10, 3)

# Boundary Conditions
## Fixing the Base Nodes
[ops.fix(n, 1, 1, 0)
for n in [1, 2, 100]];
## Rigid floor diaphragm
ops.equalDOF(12, 101, 1)
ops.equalDOF(27, 102, 1)

# Transformations & Integration
## Transformation
ops.geomTransf('Linear', 1)    # For Beams
ops.geomTransf('PDelta', 2)    # For leaning Columns
## Integration scheme
ops.beamIntegration('Lobatto', 1, 1, 5)   # For Beams
ops.beamIntegration('Lobatto', 2, 2, 5)   # For Columns

## Beams
[ops.element('dispBeamColumn',
e, *connectivity[e-1], 1, 1,
'-cMass', rho*sections['W24x55']['A'])
for e in all_the_beams];
## Columns
## -> OpenseesPy cannot handle numpy int types
## -> so I had to convert them to primitive python int type
[ops.element('dispBeamColumn',
e, *connectivity[e-1], 1, 2,
'-cMass', rho*sections['W14x120']['A'])
for e in list(map(int, all_the_cols))];

## Leaning Columns
ops.element('elasticBeamColumn', nel+1, *[100, 101],
leaning_col['A'], E_0, leaning_col['I1'], 2)
ops.element('elasticBeamColumn', nel+2, *[101, 102],
leaning_col['A'], E_0, leaning_col['I1'], 2)


### Draw Model¶

The model can now be drawn using eSEESminiPy:

In [7]:
eSEESminiPy.drawModel()


### Damping Model¶

The model assumes 2% damping for the first and second modes of the system according to rayleigh's damping model. for two modes, the damping coefficients can be obtained by:

$$\left( \begin{array}{l} {\alpha _0}\\ {\alpha _1} \end{array} \right) = 2\frac{{{\omega _m}{\omega _n}}}{{\omega _n^2 - \omega _m^2}}\left[ {\begin{array}{*{20}{c}} {{\omega _n}}&{ - {\omega _m}}\\ { - 1/{\omega _n}}&{1/{\omega _m}} \end{array}} \right]\left( \begin{array}{l} {\zeta _m}\\ {\zeta _n} \end{array} \right)\$$

So, we need to perform an eigen analysis to obtain first two natural frequencies.

In [7]:
# Building Rayleigh damping model
omega = np.sqrt(ops.eigen('-fullGenLapack', 2))
print('Two first periods are:', 2*np.pi/omega)

a_m, b_k = 2*((omega[0]*omega[1])/(omega[1]**2-omega[0]**2))*(
np.array([[omega[1],-omega[0]],
[-1/omega[1],1/omega[0]]])@np.array([0.02,0.02]))
## Rayleigh damping based on initial stiffness
ops.rayleigh(a_m, 0, b_k, 0)

Two first periods are: [0.61688251 0.11741945]


## Analysis¶

### Gravity Analysis¶

In [8]:
# Time Series
ops.timeSeries('Linear', 1)                    # For gravitional loads

ops.pattern('Plain', 1, 1)

# Settings
ops.constraints('Transformation')
ops.numberer('RCM')
ops.system('ProfileSPD')
ops.test('NormUnbalance', 0.000001, 100)
ops.algorithm('Newton')
ops.analysis('Static')

# Perform static analysis
ops.analyze(10)

Out[8]:
0

### Time History Analysis¶

In [9]:
# Set time to zero
ops.wipeAnalysis()

# Time Series
ops.timeSeries('Path', 2, '-dt', dt,           # For EQ
'-values', *northridge[:,1],
'-time', *northridge[:,0])

ops.pattern('UniformExcitation', 2, 1, '-accel', 2)

# Settings
ops.constraints('Plain')
ops.numberer('RCM')
ops.system('ProfileSPD')
ops.test('NormUnbalance', 0.0000001, 100)
ops.algorithm('Newton')
ops.integrator('Newmark', 0.5, 0.25)
ops.analysis('Transient')

# Record some responses to plot
time_lst =[]     # list to hold time stations for plotting
d_lst = []      # list to hold roof displacments
for i in range(len(northridge)):
ops.analyze(1, dt)
time_lst.append(ops.getTime())
d_lst.append(ops.nodeDisp(27,1))


## Visualization¶

Time history of the horizontal displacement of the roof is plotted here

In [10]:
plt.figure(figsize=(12,4))
plt.plot(time_lst, np.array(d_lst), color = '#d62d20', linewidth=1.75)
plt.ylabel('Horizontal Displacement (m)', {'fontname':'Cambria', 'fontstyle':'italic','size':14})
plt.xlabel('Time (sec)', {'fontname':'Cambria', 'fontstyle':'italic','size':14})
plt.grid()
plt.yticks(fontname = 'Cambria', fontsize = 14)
plt.xticks(fontname = 'Cambria', fontsize = 14);


### References¶

• Kolay, C. and Ricles, J.M., 2016. Assessment of explicit and semi‐explicit classes of model‐based algorithms for direct integration in structural dynamics. International Journal for Numerical Methods in Engineering, 107(1), pp.49-73.