Third-order shear deformation theory with a laminate offset#

This page documents the displacement field of the models based on Reddy’s third-order shear deformation theory (TSDT), 'plate_tsdt_donnell', 'cylshell_tsdt_donnell' and 'cylshell_tsdt_sanders', when the laminate is offset from the reference surface, Shell.offset, and how panels with different offsets are connected, see panels.models and panels.multidomain.connections.kCsdt.

The coordinate \(z\) is measured from the reference surface, \(z = 0\), and the laminate of thickness \(h\) occupies \(d - h/2 \le z \le d + h/2\), with the offset \(d\), the convention of composites.laminated_plate(), whose stiffness matrices are the integrals through the thickness from \(z = d - h/2\).

Displacement field#

The in-plane displacements are sought as

(1)#\[u(z) = u + z \phi_x + p(z) \left(\phi_x + w_{,x}\right) \qquad v(z) = v + z \Phi_y + p(z) \left(\phi_y + w_{,y}\right)\]

with \(p(z)\) a cubic polynomial, \(\Phi_y = \phi_y + v/r\) with the Sanders-Koiter kinematics and \(\Phi_y = \phi_y\) otherwise. The conditions on \(p(z)\) are:

  • \(p(0) = 0\), such that \(u\) and \(v\) are the displacements of the reference surface, as for the classical (CLPT) and the first-order (FSDT) theories;

  • \(p'(d \pm h/2) = -1\), such that the transverse shear strain \(\gamma_{xz} = u_{,z} + w_{,x} = (\phi_x + w_{,x}) (1 + p'(z))\) vanishes at both faces.

These three conditions leave one free parameter in the cubic \(p(z)\). It does not change the family of displacement fields: \(\gamma_{xz}\) is a quadratic polynomial in \(z\) vanishing at both faces, therefore \(\gamma_{xz} = \hat{\gamma}_{xz} (1 - 4 (z - d)^2/h^2)\) for any choice, and the integration of \(u_{,z} = \gamma_{xz} - w_{,x}\) with \(u(0) = u\) gives

(2)#\[u(z) = u - z w_{,x} + \hat{\gamma}_{xz} \left[ z - \frac{4}{3 h^2} \left( (z - d)^3 + d^3 \right) \right]\]

which depends only on \(u\), \(w\) and the amplitude \(\hat{\gamma}_{xz}\) of the transverse shear strain. The free parameter only selects which rotation of the normal is the unknown \(\phi_x = u_{,z}(z_\phi) = -w_{,x} + \hat{\gamma}_{xz} (1 - 4 (z_\phi - d)^2/h^2)\), i.e. the coordinate \(z_\phi\) where it is measured:

  • at the reference surface, \(z_\phi = 0\): \(\hat{\gamma}_{xz} = (\phi_x + w_{,x})/(1 - 4 d^2/h^2)\), singular for \(d = \pm h/2\). When the reference surface is a face of the laminate, the transverse shear vanishes there and \(u_{,z}(0) = -w_{,x}\), which cannot be an independent variable; the determinant of the relation, proportional to \((2 d - h)(2 d + h)\), vanishes, and the formulation is ill-conditioned when the reference surface approaches a face, the common case of a skin and a stiffener connected at their common face;

  • at the mid-surface, \(z_\phi = d\): \(\hat{\gamma}_{xz} = \phi_x + w_{,x}\), regular for any offset;

  • the average through the thickness, \(\bar\phi_x = \frac{1}{h} \int u_{,z} dz = -w_{,x} + \frac{2}{3} \hat{\gamma}_{xz}\), also regular, but not Reddy’s variable for \(d = 0\).

The rotations at the mid-surface are adopted, which gives

(3)#\[p(z) = -c_1 \left[ (z - d)^3 + d^3 \right], \qquad c_1 = \frac{4}{3 h^2}\]

Reddy’s field for \(d = 0\). The conditions are solved with SymPy.

Strains, stiffness and mass matrices#

With \(\zeta = z - d\) measured from the mid-surface, the in-plane strains and the transverse shear strains are

(4)#\[\varepsilon = \left( \varepsilon^{(0)} + d \varepsilon^{(1)} + d^3 \varepsilon^{(3)} \right) + \zeta \varepsilon^{(1)} + \zeta^3 \varepsilon^{(3)} \qquad \gamma = \gamma^{(0)} + \zeta^2 \gamma^{(2)}\]

with the generalized strains of the TSDT, see panels.models.fsdt_tsdt_field. The strain energy is that of the laminate about its mid-surface, whose matrices \(A, B, D, E, F, H\) are the integrals of the stiffness of the plies multiplied by \(1, \zeta, \zeta^2, \zeta^3, \zeta^4, \zeta^6\), computed by composites with a zero offset, Shell._get_lam_mid(), for the generalized strains \(T e\), see Shell._get_lam_ABD(). In \(z\) the products of the strains would require the integrals of \(z^5\), which do not appear in \(\zeta\). Likewise, the mass matrix uses the moments of the density about the mid-surface, \(J_n = \int \rho \zeta^n d\zeta\), \(n = 0, 1, 2, 3, 4, 6\), see Shell.calc_kM().

Connections at the reference surface#

Panels with different offsets share the reference surface, where the connections of panels.multidomain.connections.kCsdt, panels.multidomain.connections.kCBFycte, panels.multidomain.connections.kCBFxcte and of the null-space method tie \(u, v, w\) and the rotations of the normal at \(z = 0\):

(5)#\[\theta_x = u_{,z}\big|_{z=0} = \phi_x - \frac{4 d^2}{h^2} \left(\phi_x + w_{,x}\right) \qquad \theta_y = v_{,z}\big|_{z=0} = \Phi_y - \frac{4 d^2}{h^2} \left(\phi_y + w_{,y}\right)\]

computed with the offset \(d\) and the thickness \(h\) of each panel, see kCsdt._rotation_x() and kCsdt._rotation_y(). For \(d = \pm h/2\) they become \(\theta_x = -w_{,x}\) and \(\theta_y = \Phi_y - \phi_y - w_{,y}\). The rotations \(\phi_x, \phi_y\) belong to the mid-surfaces, at the coordinates \(z = d\) of each panel, and are not continuous across the connection of panels with different offsets. For the FSDT the rotation is constant through the thickness, \(\theta = \phi\), and for the CLPT it is \(-w_{,x}\) and \(-w_{,y}\) for any offset: their connections are unchanged.

Prescribed rotations#

The rotations phix and phiy prescribed with Shell.add_point_pd(), Shell.add_distr_pd_fixed_x() and Shell.add_distr_pd_fixed_y(), imposed with pd_method='null-space', see Prescribed displacements imposed exactly, are the same rotations at the reference surface, \(\theta_x\) and \(\theta_y\) of Eq. (5), computed with the offset and the thickness of the panel, see nullspace.pd_constraint_matrix(). With an offset, the rotations at the mid-surface \(\phi_x\) and \(\phi_y\) differ from the prescribed values by the terms of the transverse shear strains, and for \(d = \pm h/2\) the prescribed rotations are \(-w_{,x}\) and \(\Phi_y - \phi_y - w_{,y}\). For the FSDT, and for the TSDT with \(d = 0\), the prescribed rotations are \(\phi_x\) and \(\Phi_y\), as before. The penalty method, pd_method='penalty', does not prescribe rotations.

@pytest.mark.parametrize('model', ['plate_tsdt_donnell',
                                   'cylshell_tsdt_sanders'])
@pytest.mark.parametrize('kind', ['point', 'line'])
@pytest.mark.parametrize('comp', [PHIX, PHIY])
@pytest.mark.parametrize('offset', [0.3*H, -0.4*H, 0.5*H])
def test_prescribed_rotation_with_offset(model, kind, comp, offset):
    r"""Pressure loaded panel with a prescribed rotation, at a point or
    along an edge: the rotation at the reference surface evaluated from the
    field is the prescribed value, not the rotation `\phi` at the
    mid-surface, which differs by `\alpha \gamma` with the transverse shear
    strain `\gamma \neq 0`"""
    p = _shell(model, offset)
    p.add_pressure_load(-2.e5)
    xs, ys, values = _prescribe(p, kind, comp, THETA)
    md = MultiDomain([p], conn=[], pd_method='null-space')
    c = md.linear_static()
    thx, thy, phix, phiy = _rotations(p, c, xs, ys)
    theta, phi = (thx, phix) if comp == PHIX else (thy, phiy)
    assert np.allclose(theta, values, rtol=0, atol=1e-8*THETA)
    # the shear strain is not zero at the prescription, phi != theta
    assert np.abs(phi - values).max() > 1e-3*THETA

Clamped edges#

A clamped edge keeps the section planar and unrotated, \(u(z) = 0\) for every \(z\). Eq. (2) is a cubic in \(z\) with the three constants \(u\), \(w_{,x}\) and \(\hat{\gamma}_{xz} = \phi_x + w_{,x}\), all zero: \(u = \phi_x = w_{,x} = 0\), the flags x1u = x1phix = x1wr = 0 on \(x = 0\), and then \(\theta_x = 0\) for any offset. The rotation at the reference surface alone, Eq. (5), does not clamp the edge. The face clamp, \(u = 0\) at \(z = d \mp h/2\), gives instead \(\phi_x = w_{,x}/2\) and \(u = -(d/2)(1 - 4 d^2/h^2) w_{,x}\), i.e. \(u(d) = 0\) at the mid-surface and \(\hat{\gamma}_{xz} = 3 w_{,x}/2\), leaving the section free to warp, see Face-clamped edges. Both clamps are properties of the physical section and do not depend on the offset:

def test_face_clamp_does_not_depend_on_the_reference_surface():
    r"""The faces are physical surfaces: the face-clamped homogeneous plate
    gives the same deflection for the reference surface at the
    mid-surface, `d = 0`, and at `d = 0.3 h` or `d = h/2`, as the planar
    clamp of the flags"""
    model = 'plate_tsdt_donnell'
    xs, ys = np.meshgrid(np.linspace(0, A, 5), np.linspace(0, B, 5))
    res = {}
    for offset in (0., 0.3*H, 0.5*H):
        s = plate(model, offset)
        for x in (0., A):
            face_clamp(s, x, H)
        _, c = static(s)
        t = plate(model, offset)
        clamp_by_flags(t, 'x1')
        clamp_by_flags(t, 'x2')
        _, ct = static(t)
        res[offset] = (fields(s, c, xs.ravel(), ys.ravel())['w'],
                       fields(t, ct, xs.ravel(), ys.ravel())['w'])
    for offset in (0.3*H, 0.5*H):
        for i in (0, 1):
            assert np.allclose(res[offset][i], res[0.][i], rtol=1e-7, atol=0)
    assert res[0.][1].max() < 0.99*res[0.][0].max()

Example#

Two plates with the offsets \(d_1 = 0\) and \(d_2\), connected along \(x = x_s\) with the null-space method and loaded by a pressure, have continuous rotations at the reference surface, whereas the rotations at the mid-surfaces jump across the connection:

@pytest.mark.parametrize('d2', [0.3*H, 0.5*H, -0.5*H])
def test_reference_rotation_continuous_with_different_offsets(d2):
    r"""With the null-space method the rotations at the common reference
    surface are continuous across the edge, not those at the mid-surfaces,
    which differ by `-\alpha \gamma_{xz}` for `d_1 \neq d_2`"""
    md, p1, p2 = _two_plates(0., d2, conn_method='null-space')
    p1.add_pressure_load(-1.e4)
    p2.add_pressure_load(-1.e4)
    K = md.reduce(md.calc_kC())
    c = md.expand(solve(K, md.reduce(md.calc_fext()), silent=True))
    ys = np.linspace(0.05, 0.95, 7)
    th1x, th1y, phi1x, _ = _rotations(p1, c, np.full_like(ys, p1.a), ys)
    th2x, th2y, phi2x, _ = _rotations(p2, c, np.zeros_like(ys), ys)
    scale = np.abs(th1x).max()
    assert scale > 0
    assert np.allclose(th1x, th2x, rtol=0, atol=1e-9*scale)
    assert np.allclose(th1y, th2y, rtol=0, atol=1e-9*scale)
    # the rotations at the mid-surfaces are not continuous
    assert np.abs(phi1x - phi2x).max() > 1e-3*scale

and two plates with the same offset reproduce the single plate:

@pytest.mark.parametrize('conn_method', ['penalty', 'null-space'])
def test_same_offset_split_plate_equals_single_plate(conn_method):
    r"""Two panels with the same offset reproduce the single panel"""
    d = 0.3*H
    md, p1, p2 = _two_plates(d, d, conn_method=conn_method, m=14)
    if conn_method == 'null-space':
        K, M = md.reduce(md.calc_kC()), md.reduce(md.calc_kM(silent=True))
    else:
        K, M = md.calc_kC(), md.calc_kM(silent=True)
    omegas = _frequencies(K, M, 6)
    single = _plate(1., 1., d, 'SSSS', m=18, n=18)
    ref = _frequencies(single.calc_kC(), single.calc_kM(), 6)
    assert np.allclose(omegas, ref, rtol=1e-4)

For a free homogeneous plate, the choice of the reference surface does not change the structure, and the natural frequencies do not depend on the offset:

def test_free_plate_frequencies_do_not_depend_on_the_reference_surface():
    r"""A homogeneous free plate is the same structure for any choice of the
    reference surface, whose Ritz spaces coincide for free edges"""
    f0 = _free_plate_freqs(0.)
    for offset in (0.01, -0.013, 0.025):
        assert np.allclose(_free_plate_freqs(offset), f0, rtol=1e-7)