G2 Element 7

G2Element7.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.
# ------------------------------------------------------------------------
# Element 7: Perfectly plastic hinges form at element ends

import numpy as np

class element7():

    # -------------------------------------------------
    # Constructor
    # -------------------------------------------------
    def __init__(self, E, A, I, My, elno):

        self.type = "element7"

        # Store input
        self.no    = elno
        self.E     = E
        self.My    = My
        self.A     = A
        self.I     = I

        # No trial and committed variables are present because the element formulation is
        # in terms of total displacement; path-dependent behavior is not represented

    # -------------------------------------------------
    # Initialize
    # -------------------------------------------------
    def initialize(self, xyz):

        # Geometry restricted to 1,2 plane (x,y)
        dx = xyz[1, :] - xyz[0, :]

        L = np.sqrt(dx.dot(dx))

        # Check length
        if L < 1e-10:
            print('\n'"ERROR: Element number", self.no, "has zero length")
            import sys
            sys.exit()

    # -------------------------------------------------
    # State determination
    # -------------------------------------------------
    def state(self, xyz, ugMatrix, theLambda):

        # Extract *TOTAL* displacements (cannot deal with path-dependent behaviour)
        ug = ugMatrix[:, 0]

        # Geometry
        dx = xyz[1, :] - xyz[0, :]
        L = np.sqrt(dx.dot(dx))
        dx = dx / L

        # Form kinematic matrix
        Tbg = np.array([[-dx[1] / L, dx[0] / L, 1.0, dx[1] / L, -dx[0] / L, 0.0],
                        [-dx[1] / L, dx[0] / L, 0.0, dx[1] / L, -dx[0] / L, 1.0],
                        [-dx[0], -dx[1], 0.0, dx[0], dx[1], 0.0]])

        # Determine basic deformations from global deformations
        ub = Tbg.dot(ug)

        # Calculate EA/L and EI/L
        EAL = self.E * self.A / L
        EIL = self.E * self.I / L

        # Basic stiffness matrix, for different scenarios
        KbElastic = np.array([[4*EIL, 2*EIL, 0.0],
                              [2*EIL, 4*EIL, 0.0],
                              [0.0,   0.0,   EAL]])
        KbYieldEnd1 = np.array([[0.0, 0.0, 0.0],
                                [0.0,  3*EIL, 0.0],
                                [0.0,  0.0,   EAL]])
        KbYieldEnd2 = np.array([[3*EIL, 0.0, 0.0],
                                [0.0,    0.0, 0.0],
                                [0.0,    0.0, EAL]])
        KbYieldBoth = np.array([[0.0, 0.0, 0.0],
                                [0.0, 0.0, 0.0],
                                [0.0, 0.0, EAL]])

        # Calculate end moments for the elastsic scenario
        Fb = np.dot(KbElastic, ub)
        M1 = Fb[0]
        M2 = Fb[1]

        # Check the largest moment against the yield moment
        if np.abs(M1) > np.abs(M2) and np.abs(M1) > self.My:

            # Event factor for first yielding
            eta1 = self.My/np.abs(M1)

            # Element forces if only one yielding takes place
            FbFirstSlope = np.dot(KbElastic, ub) * eta1
            FbSecondSlope = np.dot(KbYieldEnd1, ub) * (1.0-eta1)
            Fb = FbFirstSlope + FbSecondSlope

            # Second-slope stiffness
            Kb = KbYieldEnd1

            # Moment at the other end
            M2 = Fb[1]

            # Check if that moment causes yield
            if np.abs(M2) > self.My:

                # Second event factor
                eta2 = (self.My-np.abs(FbFirstSlope[1]))/(np.abs(M2)-np.abs(FbFirstSlope[1]))

                # Element forces with scaled-back second-slope contribution
                FbThirdSlope = (1-eta1)*(1.0-eta2)*np.dot(KbYieldBoth, ub)
                Fb = FbFirstSlope + FbSecondSlope*eta2 + FbThirdSlope

                # Third-slope stiffness, with non-zero axial contribution
                Kb = KbYieldBoth

        elif np.abs(M2) > np.abs(M1) and np.abs(M2) > self.My:

            # Event factor for first yielding
            eta2 = self.My/np.abs(M2)

            # Element forces if only one yielding takes place
            FbFirstSlope = np.dot(KbElastic, ub) * eta2
            FbSecondSlope = np.dot(KbYieldEnd2, ub) * (1.0-eta2)
            Fb = FbFirstSlope + FbSecondSlope

            # Second-slope stiffness
            Kb = KbYieldEnd2

            # Moment at the other end
            M1 = Fb[0]

            # Check if that moment causes yield
            if np.abs(M1) > self.My:

                # Second event factor
                eta1 = (self.My-np.abs(FbFirstSlope[0]))/(np.abs(M1)-np.abs(FbFirstSlope[0]))

                # Element forces with scaled-back second-slope contribution
                FbThirdSlope = (1-eta2)*(1.0-eta1)*np.dot(KbYieldBoth, ub)
                Fb = FbFirstSlope + FbSecondSlope*eta1 + FbThirdSlope

                # Third-slope stiffness, with non-zero axial contribution
                Kb = KbYieldBoth
        else:
            Fb = np.dot(KbElastic, ub)
            Kb = KbElastic

        # Transform from basic to global
        Fg = (np.transpose(Tbg)).dot(Fb)
        Kg = (np.transpose(Tbg).dot(Kb)).dot(Tbg)

        return Fg, Kg

    # -------------------------------------------------
    # Commit
    # -------------------------------------------------
    def commit(self, xyz, u):

        # No action here because this element has a "total displacement" material model
        return

    # -------------------------------------------------
    # Print response
    # -------------------------------------------------
    def getResponse(self, xyz):

        print("Element %d (%s): No information" % (self.no, self.type))