Follower (hydrostatic) pressure#

A pressure added with follower=True in Shell.add_pressure_load() stays normal to the deformed mid-surface and acts on its deformed area. In the undeformed state it equals the dead load returned by Shell.calc_fext(). The configuration-dependent part of its force vector enters Shell.calc_fint(), and its load stiffness Shell.calc_kCfollower(), which is unsymmetric in general, enters Shell.calc_kT(), both at the load factor inc. The theory, the truncation of the area vector and the conditions for a symmetric load stiffness are derived in Follower pressure loads: theory.

Linear buckling of rings and long cylinders#

The linear buckling problem of a hydrostatic pressure is lb(kC0, kG + kCfollower), where structsolve.lb() selects the solvers of unsymmetric matrices when needed. A long cylinder buckles as a ring, modelled with a quarter of the circumference and symmetry planes:

def ring_shell(model, R, stack, plyts, laminaprops, m=4, n=15, **kwargs):
    r"""Quarter ring in generalized plane strain, see the module notes"""
    s = Shell(model=model, **kwargs)
    s.r = R
    s.a = 0.5*R
    s.b = np.pi/2*R
    s.stack = list(stack)
    s.plyts = list(plyts)
    s.laminaprops = list(laminaprops)
    s.m = m
    s.n = n
    for e in ('x1', 'x2', 'y1', 'y2'):
        for f in 'uvw':
            setattr(s, e + f, 1.)
            setattr(s, e + f + 'r', 1.)
    # symmetry planes at the ends: generalized plane strain
    s.x1u = s.x2u = 0.
    s.x1vr = s.x2vr = 0.
    s.x1wr = s.x2wr = 0.
    # symmetry planes at the straight edges, modes cos(n y/r), n even
    s.y1v = s.y2v = 0.
    s.y1wr = s.y2wr = 0.
    # rotations: phiy free at the ends, phix zero there (symmetry), phiy
    # zero at the straight edges (it varies as sin(n y/r))
    s.x1phiy = s.x2phiy = 1.
    s.x1phix = s.x2phix = 0.
    s.y1phix = s.y2phix = 1.
    s.y1phiy = s.y2phiy = 0.
    s._rebuild()
    return s

def ring_buckling_pressure(s, load='follower', zp=0., return_mode=False):
    r"""Critical external pressure of the ring ``s`` by linear buckling, as a
    positive value, the load multiplier of a unit external pressure"""
    s.clear_loads()
    s.add_pressure_load(-1., cte=False, follower=(load == 'follower'), zp=zp)
    k0 = s.calc_kC()
    c0 = solve(k0, s.calc_fext(), silent=True)
    K = s.calc_kG(c=c0)
    if load == 'follower':
        K = K + s.calc_kCfollower()
    elif load == 'central':
        K = K + central_load_stiffness(s, -1.*(1 + zp/s.r))
    eigvals, eigvecs = lb(k0, K, silent=True, num_eigvalues=4)
    if return_mode:
        return eigvals[0], eigvecs[:, 0]
    return eigvals[0]

With the Sanders kinematics the critical pressure of the mode n = 2 is 3 D/r**3 for the follower pressure and 4 D/r**3 for the dead load:

@pytest.mark.parametrize('load', ['follower', 'dead', 'central'])
@pytest.mark.parametrize('model', ['cylshell_clpt_sanders',
                                   'cylshell_clpt_donnell',
                                   'cylshell_fsdt_sanders',
                                   'cylshell_fsdt_donnell'])
def test_isotropic_ring(model, load):
    R = 1.
    h = 0.01
    s = ring_shell(model, R, [0.], [h], [ISO], fsdt_shear_correction=5/6)
    A22 = s.ABD[1, 1]
    D22 = s.ABD[4, 4]
    S = s.ABD[6, 6] if 'fsdt' in model else None
    kin = 'sanders' if 'sanders' in model else 'donnell'
    ref = -ring_closed_form(kin, load, 2, A22, D22, R, S)
    p_cr = ring_buckling_pressure(s, load)
    # the discretization error of 15 terms along y is below 1e-8
    assert np.isclose(p_cr, ref, rtol=1.e-6), (p_cr, ref)
    inextensional = {('sanders', 'follower'): 3., ('sanders', 'dead'): 4.,
                     ('sanders', 'central'): 4.5,
                     ('donnell', 'follower'): 16/5, ('donnell', 'dead'): 4.,
                     ('donnell', 'central'): 64/15}
    # r/h = 100: extensibility and transverse shear of order (h/r)**2
    assert np.isclose(p_cr*R**3/D22, inextensional[(kin, load)], rtol=1e-3)

The long sandwich cylinders of Han, Kardomateas and Simitses [han2004SandwichCylinder] are compared with their shell formulas and, by running the module, with their elasticity solution, with the pressure on the outer face given by zp = h/2:

@pytest.mark.parametrize('face', ['boron', 'graphite', 'kevlar'])
def test_han2004_shell_formulas(face):
    r"""The FSDT ring reproduces the shell formulas of Table 2"""
    for i, ratio in enumerate(RATIOS):
        elast, core_only, gbar, s8r, c3d20 = HAN2004[face][i]
        for kind, ref in (('core', core_only), ('Gbar', gbar)):
            S = han_shear_stiffness(face, kind)
            s = sandwich_ring(face, ratio, shear=S)
            p_cr = ring_buckling_pressure(s)
            # the closed form of the same discrete theory
            closed = -ring_closed_form('sanders', 'follower', 2, s.ABD[1, 1],
                                       s.ABD[4, 4], s.r, S)
            assert np.isclose(p_cr, closed, rtol=1e-6)
            # the shell formula of Han et al. uses (EI)eq with Ef, panels the
            # plane-strain D22, E/(1 - nu12 nu21) is 0.5% larger for boron
            # and less for the other faces; the extensibility of the ring is
            # also neglected by the formula
            assert np.isclose(p_cr, ref, rtol=1e-2), (face, ratio, kind,
                                                     p_cr, ref)

Non-linear analyses#

The callables of the Shell are passed directly to structsolve.Analysis, which passes its load factor to the callables that accept inc, and uses calc_fext(inc=1., c=c) as the load vector of the current configuration in the arc-length methods:

@pytest.mark.parametrize('model', ['plate_clpt_donnell',
                                   'cylshell_clpt_sanders',
                                   'cylshell_fsdt_sanders'])
def test_newton_raphson_converges_quadratically(model):
    s = follower_shell(model)
    rec = Recorder(s)
    an = Analysis(s.calc_fext, rec.calc_fint, s.calc_kC, s.calc_kG)
    rec.an = an
    an.initialInc = 0.25
    an.relTOL = 1.e-12
    an.maxNumIter = 20
    an.static(NLgeom=True, silent=True)
    assert np.isclose(an.increments[-1], 1.)
    c = an.cs[-1]
    # deflections well above the thickness, the non-linear terms matter
    h = len(s.stack)*s.plyt
    wmax = np.abs(s.uvw(c, gridx=9, gridy=9)[1]['w']).max()
    assert wmax > 2*h, wmax
    R = s.calc_fext() - s.calc_fint(c)
    assert np.linalg.norm(R) <= 1.e-12*np.linalg.norm(s.calc_fext())
    for step in range(len(an.increments)):
        errors = [e for inc, e, st in rec.calls if st == step]
        checked, orders = check_quadratic(errors)
        assert checked >= 1, (step, errors)
    # the follower terms matter: the dead load gives another solution
    s_dead = make_shell(model)
    for load in s.pressure_loads_inc:
        s_dead.pressure_loads_inc.append(load[:5])
    an_dead = Analysis(s_dead.calc_fext, s_dead.calc_fint, s_dead.calc_kC,
                       s_dead.calc_kG)
    an_dead.initialInc = 0.25
    an_dead.static(NLgeom=True, silent=True)
    assert np.isclose(an_dead.increments[-1], 1.)
    assert (np.linalg.norm(an_dead.cs[-1] - c)
            > 1.e-4*np.linalg.norm(c))

Linear static analysis and kinetic criterion#

A geometrically linear structure under a follower load solves the unsymmetric system (k0 + kCfollower) c = fext, which structsolve.Analysis does with static(NLgeom=False), since Shell.calc_fext() accepts the configuration c:

@pytest.mark.parametrize('model', ['plate_clpt_donnell', 'cylshell_fsdt_sanders'])
def test_linear_static_with_follower_load(model):
    r"""``structsolve.Analysis.static(NLgeom=False)`` evaluates the follower
    load in the current configuration of the geometrically linear structure,
    which for the first-order area vector is ``(k0 + kCfollower) c = fext``,
    a non-symmetric linear system"""
    from structsolve import static
    s = make_shell(model)
    s.add_pressure_load(-2.e5, follower=True, cte=False)
    s.add_pressure_load(lambda x, y: 1.e5*(1 + x/s.a), x1=0.1, follower=True,
                        cte=False)
    an = Analysis(s.calc_fext, s.calc_fint, s.calc_kC, s.calc_kG)
    an.relTOL = 1.e-12
    increments, cs = an.static(NLgeom=False, silent=True)
    k0 = s.calc_kC()
    kCf = s.calc_kCfollower()
    fext = s.calc_fext()
    ref = static(k0 + kCf, fext, silent=True)[1][0]
    assert np.allclose(cs[0], ref, rtol=1e-8, atol=1e-10*np.abs(ref).max())
    # the follower terms matter
    c_dead = solve(k0, fext, silent=True)
    assert np.linalg.norm(c_dead - ref) > 1e-4*np.linalg.norm(ref)

Non-conservative systems may lose stability by flutter, which only the kinetic criterion detects: structsolve.freq() solves the unsymmetric eigenvalue problem of the tangent stiffness matrix with the load stiffness, whose eigenvalues become complex at a flutter load. The cantilever cylinder of [schweizerhof1984Pressure] diverges instead, its lowest frequency vanishing at the critical load of the static criterion:

def test_kinetic_criterion_cantilever():
    r"""Kinetic criterion with :func:`structsolve.freq`, the general solver of
    unsymmetric matrices: the cantilever B1 of Schweizerhof and Ramm (1984),
    shell 1, with the unsymmetric follower load stiffness, loses stability by
    divergence (a hybrid system, Section 6.4 of the paper), not by flutter:
    the lowest natural frequency squared stays real and vanishes at the
    critical load of the static criterion"""
    from structsolve import freq
    R = 100.
    s = schweizerhof_sector(6*R, R, 1., 3, 'B1')
    s.rho = 7.8e-9
    s._rebuild()
    s.add_pressure_load(-1., cte=False, follower=True)
    k0 = s.calc_kC()
    c0 = solve(k0, s.calc_fext(), silent=True)
    KL = s.calc_kG(c=c0) + s.calc_kCfollower()
    lam_cr = np.real(lb(k0, KL, silent=True, num_eigvalues=4)[0][0])
    kM = s.calc_kM()
    w2 = []
    for frac in (0., 0.5, 0.9, 1.):
        lambda2 = freq(k0 + frac*lam_cr*KL, kM, silent=True,
                       sparse_solver=False, sort=False)[0]
        finite = np.isfinite(lambda2)
        w2.append(-lambda2[finite][np.argmin(np.abs(lambda2[finite]))])
    w2 = np.array(w2)
    assert np.all(np.abs(w2.imag) <= 1e-8*np.abs(w2[0]))
    assert w2[0].real > w2[1].real > w2[2].real > 0
    assert abs(w2[3]) < 1e-6*abs(w2[0])

Validation notebooks#

The comparisons with the literature are reproduced, per reference, by the notebooks of the folder notebooks, which use the models and the data of tests/tests_shell/test_follower_pressure_validation.py: