Get Quad4 Stiffness 2D

getQuad4Stiffness2D.py

# ------------------------------------------------------------------------
# The following function is implemented in Python 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 any form of warranty.
# ------------------------------------------------------------------------

def getQuad4Stiffness2D(E, nu, t, ps, x, y):

    # Set material stiffness
    D = np.zeros((3, 3))
    if ps == 1:  # Plane stress

        factor = E / (1.0-nu * nu)

        D[0, 0] = factor
        D[1, 0] = nu * factor
        D[2, 0] = 0.0
        D[0, 1] = nu * factor
        D[0, 2] = 0.0
        D[1, 1] = factor
        D[1, 2] = 0.0
        D[2, 1] = 0.0
        D[2, 2] = 0.5 * factor * (1.0 - nu)

    else:  # Plane strain

        factor = E / ((1.0+nu) * (1.0-2.0 * nu))

        D[0, 0] = factor * (1.0-nu)
        D[1, 0] = nu * factor
        D[2, 0] = 0.0
        D[0, 1] = nu * factor
        D[0, 2] = 0.0
        D[1, 1] = factor * (1.0 - nu)
        D[1, 2] = 0.0
        D[2, 1] = 0.0
        D[2, 2] = 0.5 * factor * (1.0 - 2.0 * nu)

    # Set quadrature rule
    order = 2
    integrationPoints = np.array([-0.577350269189626, 0.577350269189626])
    integrationWeights = np.array([1, 1])

    # Initialize loop variables
    N = np.zeros(4)
    dNdxi = np.zeros(4)
    dNdeta = np.zeros(4)
    dNdx = np.zeros(4)
    dNdy = np.zeros(4)
    iJac = np.zeros((2, 2))
    B = np.zeros((3, 8))
    Kg = np.zeros((8, 8))

    # Loop over integration points
    for k in range(order):
        for l in range(order):

            # Re-initialize the Jacobian matrix at each integration points
            jacobianMatrix = np.zeros((2, 2))

            # Get location of integeration points
            xi = integrationPoints[k]
            eta = integrationPoints[l]

            # Get shape functions, derivative and Jacobian at this quadrature point
            factor = 0.25

            N[0] = factor * (1.0-xi) * (1.0-eta)
            N[1] = factor * (1.0+xi) * (1.0-eta)
            N[2] = factor * (1.0+xi) * (1.0+eta)
            N[3] = factor * (1.0-xi) * (1.0+eta)

            dNdxi[0] = factor * -1.0 * (1.0-eta)
            dNdxi[1] = factor * 1.0 * (1.0-eta)
            dNdxi[2] = factor * 1.0 * (1.0+eta)
            dNdxi[3] = factor * -1.0 * (1.0+eta)

            dNdeta[0] = factor * (1.0-xi) * -1.0
            dNdeta[1] = factor * (1.0+xi) * -1.0
            dNdeta[2] = factor * (1.0+xi) * 1.0
            dNdeta[3] = factor * (1.0-xi) * 1.0

            for m in range(4):
                jacobianMatrix[0, 0] = jacobianMatrix[0, 0] + dNdxi[m] * x[m]
                jacobianMatrix[1, 0] = jacobianMatrix[1, 0] + dNdxi[m] * y[m]
                jacobianMatrix[0, 1] = jacobianMatrix[0, 1] + dNdeta[m] * x[m]
                jacobianMatrix[1, 1] = jacobianMatrix[1, 1] + dNdeta[m] * y[m]

            # Determine the Jacobian determinant
            jacDet = jacobianMatrix[0, 0] * jacobianMatrix[1, 1] - jacobianMatrix[0, 1] * jacobianMatrix[1, 0]

            # Determine the inverse of the Jacobian matrix
            iJac[0, 0] = 1.0 / jacDet * jacobianMatrix[1, 1]
            iJac[1, 1] = 1.0 / jacDet * jacobianMatrix[0, 0]
            iJac[0, 1] = -1.0 / jacDet * jacobianMatrix[0, 1]
            iJac[1, 0] = -1.0 / jacDet * jacobianMatrix[1, 0]

            # Determine components of the B-matrix
            for m in range(4):
                dNdx[m] = dNdxi[m] * iJac[0, 0] + dNdeta[m] * iJac[1, 0]
                dNdy[m] = dNdxi[m] * iJac[0, 1] + dNdeta[m] * iJac[1, 1]

            # Set the strain-displacement matrix B
            for m in range(4):
                B[0][(2 * m)] = dNdx[m]
                B[1][(2 * m)] = 0.0
                B[2][(2 * m)] = dNdy[m]

                B[0][(2 * m + 1)] = 0.0
                B[1][(2 * m + 1)] = dNdy[m]
                B[2][(2 * m + 1)] = dNdx[m]

            # Pre-compute Jacobian times quadrature weight (volume change)
            dvol = jacDet * t * integrationWeights[k] * integrationWeights[l]

            # Compute 8x8 stiffness matrix as B ^ T * D * B * dvol
            Kg += np.multiply((np.transpose(B).dot(D)).dot(B), dvol)

    return Kg