Linear buckling analysis

The linear buckling eigenvalue problem \(([K_C] + \lambda [K_G])\{c\} = \{0\}\) is solved with structsolve.lb(), returning the load multipliers \(\lambda\).

Constant pre-buckling stress state

The pre-buckling stress state can be defined by the constant stress resultants Nxx, Nyy and Nxy of the Shell object, used by Shell.calc_kG() to calculate the geometric stiffness matrix analytically. The example below verifies the critical loads of laminated plates and shells, simply supported on all edges or with one free edge, under axial or transverse compression:

def test_shell_lb():
    for model in ['plate_clpt_donnell',
                  'cylshell_clpt_donnell',
                  'cylshell_clpt_sanders',
                  ]:
        # ssss
        s = Shell()
        s.m = 12
        s.n = 12
        s.stack = [0, 90, -45, +45, +45, -45, 90, 0]
        s.plyt = 0.125e-3
        s.laminaprop = (142.5e9, 8.7e9, 0.28, 5.1e9, 5.1e9, 5.1e9)
        s.model = model
        s.a = 2.
        s.b = 1.
        s.r = 1.e8
        s.Nxx = -1
        eigvals, eigvecs = lb(s.calc_kC(), s.calc_kG(), silent=True)
        assert np.isclose(eigvals[0], 157.5369886, atol=0.1, rtol=0)

        s.Nxx = 0
        s.Nyy = -1
        eigvals, eigvecs = lb(s.calc_kC(), s.calc_kG(), silent=True)
        assert np.isclose(eigvals[0], 58.1488566, atol=0.1, rtol=0)

        # ssfs
        s = Shell()
        s.y2u = 1
        s.y2v = 1
        s.y2w = 1
        s.m = 12
        s.n = 13
        s.stack = [0, 90, -45, +45]
        s.plyt = 0.125e-3
        s.laminaprop = (142.5e9, 8.7e9, 0.28, 5.1e9, 5.1e9, 5.1e9)
        s.model = model
        s.a = 1.
        s.b = 0.5
        s.r = 1.e8
        s.Nxx = -1
        eigvals, eigvecs = lb(s.calc_kC(), s.calc_kG(), silent=True)
        assert np.isclose(eigvals[0], 15.8106238, atol=0.1, rtol=0)

        s.x2u = 1
        s.x2v = 1
        s.x2w = 1
        s.y2u = 0
        s.y2v = 0
        s.y2w = 0
        s.Nxx = 0
        s.Nyy = -1
        eigvals, eigvecs = lb(s.calc_kC(), s.calc_kG(), silent=True)
        assert np.isclose(eigvals[0], 13.9105, atol=0.1, rtol=0)
        plot_shell(s, eigvecs[:, 0], vec='w')

Isotropic materials are defined with laminaprop = (E, nu), and combined load cases are defined by more than one stress resultant, e.g. compression and shear:

def test_lb_isotropic():
    #NOTE ssss boundary conditions by default
    s = Shell()
    s.m = 12
    s.n = 12
    s.stack = [0, 90, -45, +45, +45, -45, 90, 0]
    thickness = 1e-3
    E = 70e9
    nu = 0.33
    s.a = 2.
    s.b = 1.
    s.r = 2.

    # compression
    # it is possible to apply any static load if necessary
    s.Nxx = -1
    # shear
    s.Nxy = -1

    s.plyt = thickness
    s.laminaprop = (E, nu)
    s.model = 'cylshell_clpt_donnell'
    # radius

    eigvals, eigvecs = lb(s.calc_kC(), s.calc_kG(), silent=True)
    plot_shell(s, eigvecs[:, 0], vec='w', filename='example_isotropic.png')

Pre-buckling stress state from a static analysis

The geometric stiffness matrix can also be calculated from the Ritz constants of a linear static solution, passed as c to Shell.calc_kG(), which then integrates the pre-buckling stress field numerically:

def test_panel_fkG_num():
    for model in ['plate_clpt_donnell',
                  'cylshell_clpt_donnell',
                  'cylshell_clpt_sanders',
                  'plate_fsdt_donnell',
                  'plate_tsdt_donnell']:
        print('Checking fkG_num for model {0}'.format(model))
        # ssss
        p = Shell()
        p.a = 8.
        p.b = 4.
        p.r = 1.e8
        p.stack = [0, 90, 90, 0, -45, +45]
        p.plyt = 1e-3*0.125
        p.laminaprop = (142.5e9, 8.7e9, 0.28, 5.1e9, 5.1e9, 5.1e9)
        p.model = model

        Nxx = -1.

        p.m = 8
        p.n = 9

        p.x1ur = 1
        p.x2u = 1
        p.x2ur = 1

        p.y1u = 1
        p.y1ur = 1
        p.y2u = 1
        p.y2ur = 1

        p.x1v = 0
        p.x1vr = 0
        p.x2v = 1
        p.x2vr = 1
        p.y1v = 1
        p.y1vr = 1
        p.y2v = 1
        p.y2vr = 1

        # from constant Nxx
        p.Nxx = Nxx
        k0 = p.calc_kC()
        kG0 = p.calc_kG()
        eigvals, eigvecs = lb(k0, kG0, silent=True)
        assert np.isclose(eigvals[0], 4.47698, atol=0.01, rtol=0)

        # from pre-buckling static solution
        p.Nxx = 0.
        p.add_distr_load_fixed_x(p.a, lambda y: Nxx, None, None)
        fext = p.calc_fext()
        c0 = solve(k0, fext)
        kG = p.calc_kG(c=c0)
        eigvals, eigvecs = lb(k0, kG, silent=True)
        assert np.isclose(eigvals[0], 4.39212, atol=0.01, rtol=0)