3. Laminate theory and shell kinematics

Chapter 2 established how a single orthotropic ply behaves under plane stress. Real composite parts, however, are stacks of plies with different fiber orientations, and their response depends on how those plies are arranged through the thickness. This chapter develops the two kinematic frameworks that connect ply-level stiffness to laminate-level behavior — Classical Lamination Theory (CLT) and First-Order Shear Deformation Theory (FSDT) — and then shows how each is discretized in a finite element shell.

Coordinate Conventions and Laminate Bookkeeping

Before any stiffness matrix can be assembled, the reference frame must be fixed. We adopt the following convention throughout this course:

• The laminate mid-surface is the z = 0 plane; z is positive upward, measured from the mid-surface.
• Ply k occupies the interval [z_k−1, z_k], with z₀ = −h/2 and z_N = +h/2 for a symmetric laminate of total thickness h.
• Ply fiber angle θ is measured counter-clockwise from the global x-axis to the fiber direction.
• Membrane forces N are force per unit width; moments M are moment per unit width.

These definitions matter because sign errors in z or θ propagate directly into the coupling matrix B and produce non-physical bending–extension coupling. Every result in this chapter uses the convention above.

Classical Lamination Theory

CLT extends Kirchhoff plate assumptions to a layered medium. The three kinematic assumptions are: straight lines normal to the mid-surface remain straight and normal after deformation; the transverse normal strain ε_zz is neglected; and each ply is in a state of plane stress. The displacement field is therefore

u(x, y, z) = u0(x, y) − z · ∂w0/∂x
v(x, y, z) = v0(x, y) − z · ∂w0/∂y
w(x, y, z) = w0(x, y)

Differentiating gives the strain distribution ε = ε⁰ + z·κ, where ε⁰ is the mid-surface membrane strain and κ is the curvature. Substituting the ply constitutive law σ = Q̄_k·ε and integrating through the thickness yields the ABD relation:

[ N ]   [ A  B ] [ ε0 ]
[ M ] = [ B  D ] [ κ  ]

A = Σ Q̄_k (z_k − z_{k−1})
B = ½ Σ Q̄_k (z_k² − z_{k−1}²)
D = ⅓ Σ Q̄_k (z_k³ − z_{k−1}³)

The extensional stiffness A governs in-plane response, the bending stiffness D governs flexure, and B couples the two. For a symmetric laminate, B = 0 and membrane and bending behavior decouple — the single most useful design rule in laminate tailoring. A quasi-isotropic [0/±45/90]_s laminate, for example, has an A matrix identical to an isotropic plate of the same thickness, which is why it is often used as a baseline for comparison.

First-Order Shear Deformation Theory

CLT neglects transverse shear deformation, which is acceptable for thin laminates but increasingly inaccurate as the span-to-thickness ratio falls below roughly 20. FSDT (also called Mindlin–Reissner theory) relaxes the normality assumption: normals remain straight but not necessarily normal to the mid-surface. Two independent rotations φ_x and φ_y are introduced:

u = u0 + z·φx
v = v0 + z·φy
w = w0

γxz = φx + ∂w0/∂x
γyz = φy + ∂w0/∂y

The constitutive relation gains transverse shear terms:

[ N ]   [ A  B   0 ] [ ε0 ]
[ M ] = [ B  D   0 ] [ κ  ]
[ Q ]   [ 0  0   H ] [ γ  ]

H = k · Σ G_k (z_k − z_{k−1})

where k is the shear correction factor (commonly 5/6 for a homogeneous section, or computed via a shear-energy equivalence for layered sections) and G_k is the ply transverse shear stiffness. The shear correction factor is a known limitation of FSDT: it is a single scalar that cannot capture the true parabolic shear stress distribution through the thickness. Higher-order theories (e.g., Reddy's third-order theory) remove this approximation at the cost of additional degrees of freedom.

Finite Element Implementation

Both theories map naturally onto shell elements. The mid-surface is discretized with 2D elements, and the through-thickness integration is handled either by the element's own integration scheme or by explicit ply-by-ply integration.

CLT shell elements use a C¹-continuous interpolation for w (because curvature involves second derivatives) or, more commonly in commercial codes, a discrete Kirchhoff formulation that enforces zero transverse shear at selected points. The element stiffness is assembled as:

K_e = ∫ B^T · [A B; B D] · B dA

where B is the strain–displacement matrix relating nodal degrees of freedom to ε⁰ and κ. The ABD matrix is evaluated once per element (or per integration point if properties vary) and passed to the quadrature routine.

FSDT shell elements are C⁰-continuous and require only first derivatives of the displacement field. They carry five degrees of freedom per node (u₀, v₀, w₀, φ_x, φ_y) and are prone to shear locking when fully integrated with low-order interpolation. Standard remedies are selective reduced integration (under-integrate the shear term) or assumed natural strain (ANS) interpolation for the transverse shear components. In practice, a four-node FSDT shell with ANS shear interpolation and 2×2 in-plane Gauss integration is the workhorse element for laminated composite analysis.

Through-thickness integration deserves explicit attention. For a laminate with N plies, the ABD matrices are computed by summing over plies, not by a fixed Gauss rule through the thickness. If the element uses a single integration point through the thickness, the ply-wise variation of Q̄ is lost and the bending stiffness will be wrong. The correct approach is to evaluate the ABD sum analytically (as in the formulas above) or to use at least one integration point per ply.

Choosing Between CLT and FSDT

The decision is driven by the span-to-thickness ratio and the quantity of interest. For a thin panel (h/L < 0.05) with in-plane loading dominant, CLT is accurate and computationally cheaper. For a thick section, a sandwich panel with a soft core, or any case where transverse shear stress or delamination-driven failure is being assessed, FSDT is required. A useful rule of thumb: if the first natural frequency predicted by CLT and FSDT differ by more than 5%, the structure is thick enough that FSDT should be used.

Worked Example: [0/90]_s Laminate Stiffness

Consider a carbon/epoxy laminate with ply properties E₁ = 135 GPa, E₂ = 10 GPa, G₁₂ = 5 GPa, ν₁₂ = 0.3, and ply thickness t = 0.125 mm. The layup is [0/90]_s, giving four plies and total thickness h = 0.5 mm. The ply z-coordinates are: z₀ = −0.25, z₁ = −0.125, z₂ = 0, z₃ = +0.125, z₄ = +0.25 mm.

Computing the reduced stiffness Q̄ for each ply orientation and summing through the thickness gives the following approximate values (in SI units, N/m and N·m):

These values are illustrative and rounded; a production analysis would compute them with full precision from the ply properties. The key observation is that B = 0, confirming the symmetric layup decouples membrane and bending response. If the layup were changed to [0/90] (unsymmetric), B would become non-zero and the panel would exhibit bending under pure in-plane load — a behavior that is often undesirable and sometimes overlooked in preliminary design.

Implementation Checklist for FE Modeling

1. Define the mid-surface geometry and mesh with shell elements (quad or tri) sized to capture the expected deformation wavelength.
2. Assign ply stacking sequence, fiber angles, and ply thicknesses using the convention from the start of this chapter.
3. Select the element formulation: CLT (thin) or FSDT (thick/shear-sensitive). Confirm the number of DOF per node matches the theory.
4. Verify that through-thickness integration captures each ply — at minimum one integration point per ply, or use the analytical ABD sum.
5. Apply boundary conditions consistent with the shell theory: for FSDT, constrain φ_x and φ_y appropriately at clamped edges; for CLT, constrain the slope ∂w/∂n.
6. Run a convergence study on mesh density and, for FSDT, on the shear correction factor if it is user-specified.
7. Post-process ply-by-ply stresses using the recovered strain field and the ply constitutive law, not the laminate-averaged stress.

The next chapter builds on this foundation by introducing failure criteria that operate on the ply-level stress state recovered from these shell analyses.