G2 Example 7

G2Example7.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.
# ------------------------------------------------------------------------
import matplotlib.pyplot as plt

from G2AnalysisNonlinearStatic import *
from G2Model import *

# Specify element type
# 5 : linear elastic frame element
# 6 : linear elastic frame element with P-delta
# 7 : linear elastic with end releases
# 8 : two-component parallel model
# 9: one-component series model with hysteretic materials
# 10: rigid interior with hysteretic materials
# 12: displacement-based distributed plasticity frame element
# 13: force-based distributed plasticity frame element
elementType = 9
PdeltaOption = "BigSmall"

# Portal frame dimensions [N, m]
L = 4.0
H = 4.0

# Rectangular cross-section dimensions [N, m]
b = 0.1
h = 0.1

# Material properties [N, m]
E = 200e9
fy = 350e6
alpha = 0.05

# Hardening in springs
springAlpha = alpha * 2

# Horizontal point load in the upper left corner of the frame [N, m]
F = 150e3

# Derived cross-section parameters
A = b*h
I = b*h**3/12

# Select yield moment
factorMy = 1.7
My = fy*I/(h/2) * factorMy

# Newton-Raphson parameters
nsteps = 10
dt = 1.0/nsteps
KcalcFrequency = 1
maxIter = 100
tol = 1e-5

# Tracked displacement
trackNode = 2
trackDOF = 1

# Nodal coordinates
NODES = [[0.0, 0.0],
         [0.0, H],
         [L,   H],
         [L,   0.0]]

# Boundary conditions (0=free, 1=fixed, sets #DOFs per node)
CONSTRAINTS = [[1, 1, 1],
               [0, 0, 0],
               [0, 0, 0],
               [1, 1, 1]]

# Lateral load, and axial force for Element 6
P = 0.5 * (np.pi**2 * E*I)/((2*L)**2)
LOADS = [[0.0, 0.0, 0.0],
         [F,   -P,  0.0],
         [0.0, -P,  0.0],
         [0.0, 0.0, 0.0]]

# Empty mass matrix
MASS = np.zeros((4, 3))

# Elements, sections, and materials, depending on the selected element type
def elementSetup(elementType, PdeltaOption):

    if elementType == 5: # Linear elastic

        q = 0.0
        ELEMENTS = [[elementType, E, A, I, q, 1, 2],
                    [elementType, E, A, I, q, 2, 3],
                    [elementType, E, A, I, q, 3, 4]]
        SECTIONS = np.zeros(3)
        MATERIALS = np.zeros(3)

    elif elementType == 6: # Linear elastic with P-delta

        q = 0.0
        ELEMENTS = [[elementType, PdeltaOption, E, A, I, q, 1, 2],
                    [elementType, PdeltaOption, E, A, I, q, 2, 3],
                    [elementType, PdeltaOption, E, A, I, q, 3, 4]]
        SECTIONS = np.zeros(3)
        MATERIALS = np.zeros(3)

    elif elementType == 7: # End releases

        ELEMENTS = [[elementType, E, A, I, My, 1, 2],
                    [elementType, E, A, I, My, 2, 3],
                    [elementType, E, A, I, My,     4, 3]]
        SECTIONS = np.zeros(3)
        MATERIALS = np.zeros(3)

    elif elementType == 8: # Parallel system model

        ELEMENTS = [[elementType, E, A, I, My, springAlpha, 1, 2],
                    [elementType, E, A, I, My, springAlpha, 2, 3],
                    [elementType, E, A, I, My, springAlpha, 4, 3]]
        SECTIONS = np.zeros(3)
        MATERIALS = np.zeros(3)

    elif elementType == 9: # Series system model with uniaxial materials

        ELEMENTS = [[elementType, 1, 2],
                    [elementType, 2, 3],
                    [elementType, 4, 3]]

        # Linear elastic stiffness in the interior
        SECTIONS = [['Interior of element type 9', E*A, E*I],
                    ['Interior of element type 9', E*A, E*I],
                    ['Interior of element type 9', E*A, E*I]]

        # Moment-rotation relationship for hinges, with high initial stiffness
        factor = 1000
        MATERIALS = [['Bilinear', 5.6*E*I/L*factor, My, springAlpha/factor /(1-springAlpha)],
                     ['Bilinear', 5.6*E*I/L*factor, My, springAlpha/factor /(1-springAlpha)],
                     ['Bilinear', 5.6*E*I/L*factor, My, springAlpha/factor /(1-springAlpha)]]

    elif elementType == 10: # Rigid interior

        ELEMENTS = [[elementType, 1, 2],
                    [elementType, 2, 3],
                    [elementType, 4, 3]]

        # Linear elastic axial stiffness in the interior
        SECTIONS = [['Interior of element type 10', E*A],
                    ['Interior of element type 10', E*A],
                    ['Interior of element type 10', E*A]]

        # Moment-rotation relationship for hinges
        MATERIALS = [['Bilinear', 5.6*E*I/L, My, springAlpha],
                     ['Bilinear', 5.6*E*I/L, My, springAlpha],
                     ['Bilinear', 5.6*E*I/L, My, springAlpha]]

    elif elementType == 12 or elementType == 13: # Displacement-based or force-based

        nsec = 5
        nfib = 10
        q = 0.0
        ELEMENTS = [[elementType, nsec, q, 1, 2],
                    [elementType, nsec, q, 2, 3],
                    [elementType, nsec, q, 4, 3]]

        SECTIONS = [['Rectangular', h, b, nfib],
                    ['Rectangular', h, b, nfib],
                    ['Rectangular', h, b, nfib]]

        MATERIALS = [['Bilinear', E, fy, alpha],
                     ['Bilinear', E, fy, alpha],
                     ['Bilinear', E, fy, alpha]]

    else:
        print("Error: Cannot understand the element type specified in the input file")
        import sys
        sys.exit()

    return ELEMENTS, SECTIONS, MATERIALS

# Create the structural model based on the above input
ELEMENTS, SECTIONS, MATERIALS = elementSetup(elementType, PdeltaOption)
a = [NODES, CONSTRAINTS, ELEMENTS, SECTIONS, MATERIALS, LOADS, MASS]
m = model(a)

# Run the analysis
loadFactor, u, dudx = nonlinearStaticAnalysis(m, nsteps, dt, maxIter, KcalcFrequency, tol, trackNode, trackDOF, [], True)

# Compare elements, if requested
compareElements = True
if compareElements:
    plt.ion()
    plt.figure()
    plt.title('Different Response for Different Element Types')
    plt.grid(True)
    plt.autoscale(True)
    plt.ylabel('Load Factor')
    plt.xlabel('Displacement')
    colors = ['black','blue','green','red','cyan','magenta','gray','yellow', 'orange', 'lime']
    colorCounter = 0
    for elementType in [5, 6, 7, 8, 9, 10, 12, 13]:

        print('\n'"Now analyzing Element Type %i" % elementType)

        if elementType == 6:

            ELEMENTS, SECTIONS, MATERIALS = elementSetup(elementType, PdeltaOption)
            a = [NODES, CONSTRAINTS, ELEMENTS, SECTIONS, MATERIALS, LOADS, MASS]
            m = model(a)
            loadFactor, u, dudx = nonlinearStaticAnalysis(m, nsteps, dt, maxIter, KcalcFrequency, tol, trackNode, trackDOF, [], False)
            plt.plot(u, loadFactor, '-', color=colors[colorCounter], linewidth=1.0, label=('Element ' + str(elementType) + PdeltaOption))
            colorCounter += 1

            PdeltaOption = "Big"
            ELEMENTS, SECTIONS, MATERIALS = elementSetup(elementType, PdeltaOption)
            a = [NODES, CONSTRAINTS, ELEMENTS, SECTIONS, MATERIALS, LOADS, MASS]
            m = model(a)
            loadFactor, u, dudx = nonlinearStaticAnalysis(m, nsteps, dt, maxIter, KcalcFrequency, tol, trackNode, trackDOF, [], False)
            plt.plot(u, loadFactor, '-', color=colors[colorCounter], linewidth=1.0, label=('Element ' + str(elementType) + PdeltaOption))
            colorCounter += 1

        else:

            ELEMENTS, SECTIONS, MATERIALS = elementSetup(elementType, PdeltaOption)
            a = [NODES, CONSTRAINTS, ELEMENTS, SECTIONS, MATERIALS, LOADS, MASS]
            m = model(a)
            loadFactor, u, dudx = nonlinearStaticAnalysis(m, nsteps, dt, maxIter, KcalcFrequency, tol, trackNode, trackDOF, [], False)
            if elementType == 7:
                plt.plot(u, loadFactor, 'o-', color=colors[colorCounter], linewidth=1.0, label=('Element ' + str(elementType)))
                colorCounter += 1
            else:
                plt.plot(u, loadFactor, '-', color=colors[colorCounter], linewidth=1.0, label=('Element ' + str(elementType)))
                colorCounter += 1

    plt.legend()
    print('\n'"Click somewhere in the plot to continue...")
    plt.waitforbuttonpress()