Shell and beam properties#

Shell properties and laminates#

The shell elements read their constitutive matrices and mass integrals from a ShellProp, passed to every update_* method that needs it. The stiffnesses of the first-order shear deformation theory are

\[\begin{split}\begin{Bmatrix} N \\ M \end{Bmatrix} = \begin{bmatrix} A & B \\ B & D \end{bmatrix} \begin{Bmatrix} \varepsilon^{(0)} \\ \varepsilon^{(1)} \end{Bmatrix}, \qquad \begin{Bmatrix} Q_y \\ Q_x \end{Bmatrix} = \boldsymbol{A}_{ts} \begin{Bmatrix} \gamma_{yz} \\ \gamma_{xz} \end{Bmatrix}\end{split}\]

with \(A_{ij}\), \(B_{ij}\), \(D_{ij}\) the thickness integrals of the ply stiffnesses times \(1\), \(z\) and \(z^2\), and \(\boldsymbol{A}_{ts}\) = [[A44, A45], [A45, A55]], see ShellProp. A property is created by:

The ply stiffnesses are those of plane stress, see Lamina.rebuild(). Plane strain and thermal expansion are not modelled: the thermal expansion coefficients and the reference temperature of MatLamina are not used.

The ply angles are measured from the material direction of each element, the vector xmati, xmatj, xmatk of update_rotation_matrix projected onto the element, or from the element \(x\) axis when it is not given, see Conventions. The elements rotate the matrices from the material direction to the element axes with m11, m12, m21, m22; ShellProp.calc_constitutive_element() gives the rotated matrices.

Offset of the reference surface#

The nodes of the shell elements are on the reference surface, and \(z\) is measured from it. With offset \(= e\), the laminate occupies \(-h/2 + e \le z \le h/2 + e\), i.e. its mid-surface is at the distance \(e\) from the nodes, along the element \(z\) axis. For a single isotropic layer the integrals give

\[\boldsymbol{B} = e \boldsymbol{A}, \qquad \boldsymbol{D} = \boldsymbol{D}_{e=0} + e^2 \boldsymbol{A}, \qquad \int \rho z \, dz = \rho h e, \qquad \int \rho z^2 \, dz = \rho \left(\frac{h^3}{12} + h e^2\right)\]

such that an offset couples the membrane and bending responses. It models e.g. a skin whose nodes are on its inner surface, shared with the beams of a stiffener. The mass integrals intrho, intrhoz and intrhoz2 enter the mass matrices of the shells.

Transverse shear correction#

The stiffnesses A44, A45 and A55 of a ShellProp are the ones used by the elements, with the shear correction already applied, and the elements apply no other correction factor. The correction is selected with shear_correction, see ShellProp.calc_transverse_shear_stiffness():

  • 'rohwer', the default: the equilibrium approach of Rohwer (1988), giving the full matrix \(\boldsymbol{A}_{ts}\), re-evaluated in the element axes by ShellProp.calc_Ats_element() when the material direction differs from the element \(x\) axis;

  • 'vlachoutsis': the two factors of Vlachoutsis (1992), doi:10.1002/nme.1620330712, identical to 'whitney' (Whitney, 1973) and, for laminates symmetric about the mid-surface, to 'chow' (Chow, 1971);

  • 'birman_bert': the average shear strain factors of Birman and Bert (2002);

  • 'thickness_shear': the factors that match the first thickness-shear frequency of the laminate, which require the ply densities;

  • 'constant': \(k = 5/6\);

  • None: no correction, \(\boldsymbol{A}_{ts} = \sum_k \boldsymbol{C}_s^{(k)} h_k\).

For an isotropic plate 'rohwer', 'vlachoutsis', 'whitney', 'chow' and 'constant' give \(5/6 \, G h\), 'birman_bert' gives \(G h\) and 'thickness_shear' gives \(\pi^2/12 \, G h\). Except for 'rohwer', the elements rotate \(\boldsymbol{A}_{ts}\) to the element axes as a tensor. None of the corrections depends on offset. The uncorrected stiffnesses Abar44, Abar45, Abar55 and the constant-stress ones Abarbar44, Abarbar45, Abarbar55 are kept for comparison, and scf_k13, scf_k23 are the resulting ratios, informative only. A ShellProp created from lamination parameters has no plies and no correction; A44, A45, A55 can then be set directly.

The transverse shear stresses through the thickness are recovered with ShellProp.calc_transverse_shear_stress(), from the shear forces and the distribution of Rohwer (1988), or with ShellProp.calc_equilibrium_transverse_shear(), from the derivatives of the generalized strains of a solution, whose strain energy is given by ShellProp.calc_aposteriori_energy().

Idealized laminates#

After ShellProp.calc_constitutive_matrix(), the methods ShellProp.make_balanced(), ShellProp.make_orthotropic(), ShellProp.make_symmetric() and ShellProp.make_smeared() set to zero the stiffnesses that vanish for a balanced, orthotropic or symmetric laminate, or replace \(\boldsymbol{D}\) by \(h^2 \boldsymbol{A}/12\), which requires offset = 0. make_balanced_LP(), make_orthotropic_LP() and make_symmetric_LP() do the same with the lamination parameters, such that shellprop_from_LaminationParameters() gives the same stiffnesses:

def test_make_orthotropic_LP_transverse_shear():
    # G13 != G23, such that xiAts2 gives Abar45
    laminaprop = (142.5e9, 8.7e9, 0.28, 5.1e9, 5.1e9, 3.0e9)
    prop = laminated_plate([30, -45, 60, 60, -45, 30], plyt=0.125e-3,
                           laminaprop=laminaprop)
    matlamina = prop.plies[0].matlamina
    lp = prop.calc_lamination_parameters()
    assert abs(lp.xiAts2) > 0
    prop.make_orthotropic()
    make_orthotropic_LP(lp)
    assert lp.xiAts2 == 0
    prop_2 = shellprop_from_LaminationParameters(prop.h, matlamina, lp)
    assert np.allclose(prop_2.A, prop.A)
    assert np.allclose(prop_2.B, prop.B)
    assert np.allclose(prop_2.D, prop.D)
    assert np.allclose(prop_2.Abar_ts, prop.Abar_ts)

Beam properties#

BeamProp holds the cross-section integrals, in the element axes \(y\) and \(z\), of BeamC, BeamLR and Truss: A, Iyy \(= \int z^2 dA\), Izz \(= \int y^2 dA\), Iyz, the torsion constant J, Ay and Az, which are null when the beam axis passes through the centroid, the moduli E and G, and the mass integrals intrho, intrhoy, intrhoz, intrhoy2, intrhoz2 and intrhoyz. Bending in the element \(xy\) plane, i.e. the displacement \(v\), is governed by Izz and by the transverse shear modulus G12, bending in the \(xz\) plane by Iyy and G13, and the torsion by the torsional stiffness GJ. These three attributes follow G until they are set, G12 = G13 = G and GJ = G*J, such that a property given with G and J only keeps its meaning:

def test_defaults_follow_G():
    prop = BeamProp()
    assert prop.G12 == prop.G13 == prop.GJ == 0
    prop.G = 3.
    prop.J = 5.
    assert prop.G12 == 3. and prop.G13 == 3. and prop.GJ == 15.
    prop.G12 = 1.
    prop.G13 = 2.
    prop.GJ = 7.
    prop.G = 4.
    prop.J = 6.
    assert prop.G12 == 1. and prop.G13 == 2. and prop.GJ == 7.
    for prop2 in [pickle.loads(pickle.dumps(prop)), copy.deepcopy(prop)]:
        assert prop2.G12 == 1. and prop2.G13 == 2. and prop2.GJ == 7.
        assert prop2.G == 4. and prop2.J == 6.
    del prop.G12
    del prop.G13
    del prop.GJ
    assert prop.G12 == 4. and prop.G13 == 4. and prop.GJ == 24.

No correction is applied by the elements: the shear correction factors of the two planes, e.g. \(5/6\) for a rectangular section, are applied by the user to G12 and G13. Applied to G, as in some of the examples, they also scale GJ, unless GJ is set.

Examples#

The constitutive matrices of an isotropic plate, and of the same plate with an offset:

def test_isotropic_plate():
    prop = helper_isotropic_plate()
    A = np.array([[9629991.31944444, 2696397.56944444,       0.   ],
                  [2696397.56944444, 9629991.31944444,       0.   ],
                  [      0.        ,       0.        , 3466796.875]])
    D = np.array([[0.01253905, 0.00351093, 0.        ],
                  [0.00351093, 0.01253905, 0.        ],
                  [0.        , 0.        , 0.00451406]])
    Abar_ts = np.array([[3466796.875,       0.   ],
                  [      0.   , 3466796.875]])
    assert np.allclose(prop.A, A)
    assert np.allclose(prop.B, 0)
    assert np.allclose(prop.D, D)
    assert np.allclose(prop.Abar_ts, Abar_ts)
    assert np.allclose(prop.Ats, 5/6*Abar_ts)

def test_offset_isotropic_plate():
    r"""With the offset `e`, the plate occupies `-h/2 + e \le z \le h/2 + e`
    and the thickness integrals give `B = e A`, `D = D_0 + e^2 A`,
    ``intrhoz`` `= \rho h e` and ``intrhoz2`` `= \rho (h^3/12 + h e^2)`; the
    transverse shear stiffness does not change"""
    E = 71e9
    nu = 0.28
    rho = 2700.
    h = 0.002
    e = 0.6*h
    prop0 = isotropic_plate(thickness=h, E=E, nu=nu, rho=rho)
    prop = isotropic_plate(thickness=h, E=E, nu=nu, rho=rho, offset=e)
    assert np.allclose(prop.A, prop0.A)
    assert np.allclose(prop.B, e*prop0.A, rtol=1e-12)
    assert np.allclose(prop.D, prop0.D + e**2*prop0.A, rtol=1e-12)
    assert np.allclose(prop.Ats, prop0.Ats, rtol=1e-12)
    assert np.isclose(prop.intrho, rho*h)
    assert np.isclose(prop.intrhoz, rho*h*e)
    assert np.isclose(prop.intrhoz2, rho*(h**3/12 + h*e**2))

The shear correction of an isotropic plate, and the stiffnesses of a cross-ply laminate given by Rohwer (1988):

def test_isotropic_five_sixths():
    E = 71e9
    nu = 0.33
    h = 0.002
    G = E/(2*(1 + nu))
    for offset in [0., 0.37*h]:
        for mode in ['rohwer', 'vlachoutsis', 'constant']:
            prop = isotropic_plate(thickness=h, E=E, nu=nu, offset=offset,
                                   shear_correction=mode)
            assert np.allclose(prop.Ats, 5/6*G*h*np.eye(2), rtol=1e-12)
            assert np.isclose(prop.scf_k13, 5/6)
            assert np.isclose(prop.scf_k23, 5/6)
            assert np.allclose(prop.Abar_ts, G*h*np.eye(2))
            assert np.allclose(prop.Abarbar_ts, G*h*np.eye(2))
        prop = isotropic_plate(thickness=h, E=E, nu=nu, offset=offset,
                               shear_correction=None)
        assert np.allclose(prop.Ats, G*h*np.eye(2))
        assert np.isclose(prop.scf_k13, 1.)

def test_rohwer_cross_ply_reference_values():
    # NOTE Rohwer (1988), and Table plate-shear-rotation of the theory
    prop = laminated_plate(stack=CROSS_PLY, plyt=PLYT, laminaprop=CFRP)
    assert np.isclose(prop.Abar44, 3.450, rtol=1e-4)
    assert np.isclose(prop.Abar55, 3.450, rtol=1e-4)
    assert np.isclose(prop.Abarbar44, 3.067, rtol=1e-3)
    assert np.isclose(prop.Abarbar55, 3.067, rtol=1e-3)
    assert np.isclose(prop.A44, 2.5252, rtol=1e-4)
    assert np.isclose(prop.A55, 2.3131, rtol=1e-4)
    assert np.isclose(prop.A45, 0., atol=1e-12)
    assert np.isclose(prop.scf_k13, prop.A55/prop.Abar55)
    assert np.isclose(prop.scf_k23, prop.A44/prop.Abar44)

Rohwer, K., Improved transverse shear stiffness for layered finite elements, DFVLR-FB 88-32, 1988.