G2 Example 14

G2Example14.py

# ------------------------------------------------------------------------
# The following Python code is implemented by Professor Terje Haukaas at
# the University of British Columbia in Vancouver, Canada. It is made
# freely available online at terje.civil.ubc.ca together with notes,
# examples, and additional Python code. Please be cautious when using
# this code; it may contain bugs and comes without warranty of any kind.
# ------------------------------------------------------------------------

# Input [N, m, kg, sec]
elementType = 12
L = 15.0
nel = 5
E = 200e9
fy = 350e6
alpha = 0.02
rho = 7850.0
hw = 0.355
bf = 0.365
tf = 0.018
tw = 0.011
nf = 2
nw = 8
nsec = 5
maxIter = 100
tol = 1e-6
trackNodes = [nel+1]
trackDOFs = [1]

from GroundMotionHalfSineWave import *
gmdt = 0.005
T = 0.4
amplitudeINg = 20
numHalfSineWaves = 1
duration = 4
groundMotionFile = 'Pulse.txt'
createHalfSineWave(T, numHalfSineWaves, gmdt, amplitudeINg, groundMotionFile)
gmScaling = 1
dtAnalysis = gmdt
targetDamping = 0.05

A = tw * (hw - 2 * tf) + 2 * bf * tf
I = tw * (hw - 2 * tf) ** 3 / 12.0 + 2 * bf * tf * (0.5 * (hw - tf)) ** 2
M = A * L/nel * rho
print('\n'"Analytical first natural frequency: %.2frad" % (1.875 ** 2 * np.sqrt(E * I / (rho * A * L**4))))

# *************************** #
#          OpenSees           #
# *************************** #
import matplotlib.pyplot as plt
from sys import platform
if platform == "darwin":
    from openseespymac.opensees import *
elif platform == "win32":
    from openseespy.opensees import *
else:
    print("Cannot handle this type of operating system")
    import sys
    sys.exit()
model('basic', '-ndm', 2, '-ndf', 3)
for i in range(nel+1):
    node(i+1, 0.0, i*L/nel)
fix(1, 1, 1, 1)
uniaxialMaterial('Steel01', 1, fy, E, alpha)
section('WFSection2d', 1, 1, hw, tw, bf, tf, nw, nf)
beamIntegration('Legendre', 1, 1, nsec)
geomTransf('Linear', 1)
for i in range(nel):
    element('dispBeamColumn', i+1, i+1, i+2, 1, 1)
for i in range(nel-1):
    mass(i+2, M, 0.0, 0.0)
mass(nel+1, 0.5*M, 0.0, 0.0)
timeSeries('Path', 1, '-filePath', 'Pulse.txt', '-dt', dtAnalysis, '-factor', 9.81*gmScaling)
pattern('UniformExcitation', 1, 1, '-accel', 1)
eigenValues = eigen(2)
omega1 = np.sqrt(eigenValues[0])
omega2 = np.sqrt(eigenValues[1])
print('\n'"The first natural frequency from OpenSees is %.2f radians" % omega1)
factor = 2 * targetDamping / (omega1 + omega2)
cM = omega1 * omega2 * factor
cK = factor
rayleigh(cM, 0.0, cK, 0.0)
#modalDamping(targetDamping)
system('BandGeneral')
constraints('Plain')
test('NormDispIncr', 1.0e-12, 10)
algorithm('Newton')
numberer('RCM')
integrator('Newmark', 0.5, 0.25)
analysis('Transient')
time = [getTime()]
u = [0.0]
while getTime() < duration:
    ok = analyze(1, .02)
    time.append(getTime())
    u.append(nodeDisp(nel+1, 1))
plt.ion()
plt.figure()
plt.plot(time, u, 'k-')
plt.xlabel("Time [sec.]")
plt.ylabel("Displacement")
plt.title("OpenSees Response (max=%.3f, min=%.3f)" % (np.max(u), np.min(u)))
print('\n'"Click somewhere in the plot to continue...")
plt.waitforbuttonpress()

# *************************** #
#           G2                #
# *************************** #
from G2AnalysisNonlinearDynamic import *
from G2Model import *

NODES = []
for i in range(nel+1):
    NODES.append([0.0, i*L/nel])

CONSTRAINTS = [[1, 1, 1]]
for i in range(nel):
    CONSTRAINTS.append([0, 0, 0])

ELEMENTS = []
for i in range(nel):
    ELEMENTS.append([elementType, nsec, 0.0, i+1, i+2])

SECTIONS = []
for i in range(nel):
    SECTIONS.append(['WideFlange', hw, bf, tf, tw, nf, nw])

MATERIALS = []
for i in range(nel):
    MATERIALS.append(['Bilinear', E, fy, alpha])

LOADS = np.zeros((nel+1, 3))

MASS = [[0, 0, 0]]
for i in range(nel-1):
    MASS.append([M, 0, 0])
MASS.append([0.5*M, 0, 0])

a = [NODES, CONSTRAINTS, ELEMENTS, SECTIONS, MATERIALS, LOADS, MASS]
m = model(a)
dampingModel = ['Rayleigh', 'Initial', 'Given', cM, cK]

nonlinearDynamicAnalysis(m, dampingModel, groundMotionFile, gmScaling, gmdt, dtAnalysis, duration, trackNodes, trackDOFs, [])