Linear static analysis of beams and frames

An L-shaped frame of BeamC elements, clamped at the base of the column and loaded at the free end of the beam, in its plane and out of its plane. The orientation vector of update_rotation_matrix puts the element \(z\) axis of both members along the global \(Z\) axis, such that the in-plane bending is governed by Izz and the out-of-plane bending by Iyy, see Conventions. The displacements of the free end, with the bending, shear, axial and torsional deformations, are exact with one element per member, and the internal forces at the base are the support reactions:

def frame(num_elements, H=3., L=2.):
    r"""Nodes, elements, property and ``KC0`` of the L-shaped frame, with
    ``num_elements`` elements per member"""
    E = 200.e9
    nu = 0.3
    b = 0.05 # width along the element z axis
    h = 0.10 # height along the element y axis
    prop = BeamProp()
    prop.E = E
    # NOTE no shear correction, A is the shear area
    prop.G = E/(2*(1 + nu))
    prop.A = b*h
    prop.Izz = b*h**3/12
    prop.Iyy = h*b**3/12
    prop.J = 1.5e-6

    # NOTE column from (0, 0, 0) to (0, H, 0), then beam to (L, H, 0)
    s = np.linspace(0, 1, num_elements + 1)
    column = np.column_stack((0*s, H*s, 0*s))
    beam = np.column_stack((L*s[1:], H + 0*s[1:], 0*s[1:]))
    ncoords = np.vstack((column, beam))
    x = ncoords.flatten()
    num_nodes = ncoords.shape[0]
    N = DOF*num_nodes

    data = BeamCData()
    probe = BeamCProbe()
    num_beams = num_nodes - 1
    KC0r = np.zeros(data.KC0_SPARSE_SIZE*num_beams, dtype=INT)
    KC0c = np.zeros(data.KC0_SPARSE_SIZE*num_beams, dtype=INT)
    KC0v = np.zeros(data.KC0_SPARSE_SIZE*num_beams, dtype=DOUBLE)
    elements = []
    init_k_KC0 = 0
    for pos1 in range(num_beams):
        pos2 = pos1 + 1
        elem = BeamC(probe)
        elem.n1 = pos1 + 1
        elem.n2 = pos2 + 1
        elem.c1 = DOF*pos1
        elem.c2 = DOF*pos2
        elem.init_k_KC0 = init_k_KC0
        # NOTE the element z axis is Z for both members: in-plane bending
        #      with Izz, out-of-plane bending with Iyy
        if pos1 < num_elements:
            elem.update_rotation_matrix(-1., 0., 0., x) # column
        else:
            elem.update_rotation_matrix(0., 1., 0., x) # beam
        elem.update_probe_xe(x)
        elem.update_KC0(KC0r, KC0c, KC0v, prop)
        elements.append(elem)
        init_k_KC0 += data.KC0_SPARSE_SIZE
    KC0 = coo_matrix((KC0v, (KC0r, KC0c)), shape=(N, N)).tocsc()
    return elements, prop, x, N, KC0

@pytest.mark.parametrize('num_elements', [1, 4])
def test_l_frame_tip_displacements(num_elements):
    H = 3.
    L = 2.
    P = 1.e3
    elements, prop, x, N, KC0 = frame(num_elements, H=H, L=L)
    E, G, A = prop.E, prop.G, prop.A
    Izz, Iyy, J = prop.Izz, prop.Iyy, prop.J

    # NOTE clamped base, the first node
    bk = np.zeros(N, dtype=bool)
    bk[:DOF] = True
    bu = ~bk
    KC0uu = KC0[bu, :][:, bu]
    tip = N - DOF

    # in-plane load at the free end
    f = np.zeros(N)
    f[tip + 1] = -P
    u = np.zeros(N)
    u[bu] = solve(KC0uu, f[bu])
    v_ref = -P*(L**3/(3*E*Izz) + L/(G*A) + L**2*H/(E*Izz) + H/(E*A))
    u_ref = P*L*H**2/(2*E*Izz)
    rz_ref = -P*(L*H/(E*Izz) + L**2/(2*E*Izz))
    print('v', u[tip + 1], v_ref, 'u', u[tip + 0], u_ref)
    assert np.isclose(u[tip + 1], v_ref, rtol=1e-8)
    assert np.isclose(u[tip + 0], u_ref, rtol=1e-8)
    assert np.isclose(u[tip + 5], rz_ref, rtol=1e-8)

    # NOTE the internal forces at the base are the support reactions,
    #      balancing the load and its moment P*L about Z
    fint = np.zeros(N)
    for elem in elements:
        elem.update_probe_ue(u)
        elem.update_fint(fint, prop)
    assert np.allclose(fint[:DOF], [0, P, 0, 0, 0, P*L], atol=1e-6*P*L)
    assert np.allclose(fint[bu], f[bu], atol=1e-6*P)

    # out-of-plane load at the free end
    f = np.zeros(N)
    f[tip + 2] = P
    u = np.zeros(N)
    u[bu] = solve(KC0uu, f[bu])
    w_ref = P*((L**3 + H**3)/(3*E*Iyy) + (L + H)/(G*A) + L**2*H/(G*J))
    print('w', u[tip + 2], w_ref)
    assert np.isclose(u[tip + 2], w_ref, rtol=1e-8)

The tip deflection of a cantilever with shear deformation, Eq. (1.2) of Bathe (2014), Finite Element Procedures, 2nd ed., Sec. 1.2, p. 5, exact with one BeamC element, whereas BeamLR converges with mesh refinement:

def cantilever_tip(name, prop, L, num_elements):
    """Tip deflections along y and z under unit tip forces"""
    KC0, _ = assemble(name, prop, L, num_elements)
    KC0 = KC0.toarray()
    N = KC0.shape[0]
    bu = np.ones(N, dtype=bool)
    bu[:DOF] = False
    tip = []
    for direction in [1, 2]:
        f = np.zeros(N)
        f[N - DOF + direction] = 1.
        u = np.zeros(N)
        u[bu] = np.linalg.solve(KC0[np.ix_(bu, bu)], f[bu])
        tip.append(u[N - DOF + direction])
    return np.array(tip)

@pytest.mark.parametrize('name, num_elements, rtol', [
    ('BeamC', 1, 1e-10),
    ('BeamC', 3, 1e-10),
    # NOTE BeamLR converges quadratically; error 0.02 % with 32 elements
    ('BeamLR', 100, 1e-4),
    ])
def test_bathe_cantilever_shear_deflection(name, num_elements, rtol):
    L = 28. # cm, L + r_N in Fig. 1.2(b) of Bathe (2014)
    W = 1000. # N
    prop = bathe_prop()
    delta = W*cantilever_tip(name, prop, L, num_elements)
    # Eq. (1.2) of Bathe (2014)
    delta_ref = W*L**3/(3*prop.E*prop.Izz) + W*L/(prop.A*prop.G)
    print(name, num_elements, delta, delta_ref)
    assert np.allclose(delta, delta_ref, rtol=rtol)
    # NOTE published value, given with two significant digits
    assert np.allclose(delta, 0.053, atol=0.0005)