G2 Example 16

G2Example16.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.
# ------------------------------------------------------------------------

from G2AnalysisLinearDynamic import *
from G2Model import *
from GroundMotionHalfSineWave import *

#              |
#              | P
#              |
#              V
#        ----> * -----> F
#        ----> |
#        ----> |
#        ----> |
#        ----> |
#      q ----> | L
#        ----> |
#        ----> |
#        ----> |
#        ----> |
#        ----> |
#            -----

# Input [N, m, kg, sec]
L = 15.0           # Total length of cantilever
elementType = 5    # Linear frame element
nel = 5            # Number of elements along cantilever
E = 200e9          # Modulus of elasticity
rho = 7850.0       # Mass density
hw = 0.355         # Web height
bf = 0.365         # Flange width
tf = 0.018         # Flange thickness
tw = 0.011         # Web thickness
trackNode = nel+1  # Node to be plotted
trackDOF = 1       # DOF to be plotted

# Ground motion
dt = 0.02
duration = 1.0
groundMotionFile = 'Pulse.txt'
groundMotionPeriod = 0.5
numGroundMotionHalfSineWaves = 1
groundMotionAmplitude = 1
createHalfSineWave(groundMotionPeriod, numGroundMotionHalfSineWaves, dt, groundMotionAmplitude, groundMotionFile)

# Damping
dampingModel = ['Rayleigh', 'Initial', 'UseEigenvalues', 1, 2, 0.05]

# Area, moment of inertia, nodal mass
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

# Calculate analytical eigenvalue(s)
print('\n'"Analytical first natural frequency: %.2frad" % (1.875 ** 2 * np.sqrt(E * I / (rho * A * L**4))))

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

# Boundary conditions (0=free, 1=fixed, sets #DOFs per node)
CONSTRAINTS = [[1, 1, 1]]
for i in range(nel):
    CONSTRAINTS.append([0, 0, 0])

# Element connectivity and type
ELEMENTS = []
for i in range(nel):
    ELEMENTS.append([elementType, E, A, I, 0.0, i+1, i+2])

# Empty arrays
SECTIONS = np.zeros(nel)
MATERIALS = np.zeros(nel)

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

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

# Check Young's modulus DDM
a = [NODES, CONSTRAINTS, ELEMENTS, SECTIONS, MATERIALS, LOADS, MASS]
m = model(a)
DDMparameters = [['Element', 'E', range(1, nel+1)]]
u, dudx = linearDynamicAnalysis(m, dampingModel, groundMotionFile, dt, duration, trackNode, trackDOF, DDMparameters)
perturbationFraction = 1e-8
m.setParameters([['Element', 'E', range(1, nel+1), E*(1.0+perturbationFraction)]])
responsePert, blank = linearDynamicAnalysis(m, dampingModel, groundMotionFile, dt, duration, trackNode, trackDOF, [])
fdmSensitivity = 1.0/(E*perturbationFraction) * np.subtract(responsePert, u)
print('\n'"FDM-DDM difference:", np.max(np.abs(np.subtract(fdmSensitivity, dudx[:,0]))))
plt.ion()
plt.figure()
plt.autoscale(True)
plt.title("Checking Young's Modulus Response Sensitivity")
numTimePoints = len(dudx[:,0])
t = np.linspace(0, dt * numTimePoints, numTimePoints)
plt.plot(t, fdmSensitivity, 'ro-', label='FDM')
plt.plot(t, dudx[:,0], 'k-', label='DDM')
plt.xlabel("Time [sec]")
plt.ylabel("Derivative")
plt.legend(loc='upper right')
print('\n'"Click somewhere in the plot to continue...")
plt.waitforbuttonpress()