CLT Theory and MATLAB

A formal teaching page on laminated structures, laminate stiffness, ply stress, strength assessment, and code interpretation

Teaching code package. The MATLAB scripts used in this page are available below.

Open MATLAB code page Download MATLAB code package

Contents

1. Basic definitions 2. Why CLT is needed 3. Material data and engineering meaning 4. Lamina constitutive law 5. Coordinate transformation and rotated ply stiffness 6. Laminate A, B, and D matrices 7. Load resultants, mid-plane strain, and curvature 8. Ply stress recovery 9. Strength and first-ply failure 10. Worked laminate example 11. MATLAB implementation 12. Example numerical results and how to read them 13. Recommended plots and outputs

1. Basic definitions

A ply is one thin layer of composite material. When that ply is treated as a mechanically homogeneous orthotropic layer, it is also called a lamina. A laminate is a bonded stack of plies.

The purpose of CLT is to connect ply-level mechanics to laminate-level response.

2. Why CLT is needed

For an isotropic plate, global stiffness is described by only a few constants. For a laminate, that is no longer enough. The response depends on the ply orientation set \(\{\theta_k\}\), the thickness distribution \(\{t_k\}\), and the lamina orthotropic properties \((E_1,E_2,G_{12},\nu_{12})\).

CLT is useful because it answers the first structural questions quickly:

3. Material data and engineering meaning

The basic lamina elastic constants are:

The strength parameters are:

Do not mix them. The elastic constants determine deformation and stress distribution. The strength parameters determine whether the predicted stress state is safe or unsafe.

Material-axis meaning.

  • Direction 1 is the fiber direction, usually the stiff and strong direction.
  • Direction 2 is the in-plane transverse direction, usually much softer and weaker.

This directional contrast is exactly why a composite lamina cannot be treated like an isotropic sheet.

4. Lamina constitutive law

Conditions under which the zero terms are valid.

The reduced stiffness matrix written as

\[ \mathbf{Q}= \begin{bmatrix} Q_{11} & Q_{12} & 0\\ Q_{12} & Q_{22} & 0\\ 0 & 0 & Q_{66} \end{bmatrix} \]

is not the most general 2D constitutive form. It is valid only under the following conditions:

  1. The material is treated as an orthotropic lamina.
  2. The constitutive law is written in the lamina material principal axes \((1,2)\).
  3. The in-plane response is written in the plane-stress reduced form used in Classical Lamination Theory.
  4. The ply is described in its own local material frame, not yet in a rotated laminate frame.

Under these conditions, the normal–shear coupling terms vanish, so \(Q_{16}=Q_{26}=0\). In a more general coordinate system, these terms do not have to be zero.

At the lamina level, the first task is to write the constitutive relation of a single orthotropic ply in its own material axes \((1,2)\). These axes are attached to the lamina itself. Direction \(1\) is the fiber direction, and direction \(2\) is the in-plane transverse direction.

In the material axes \((1,2)\), define the in-plane stress and strain vectors as

\[ \boldsymbol{\sigma}_{12}= \begin{bmatrix} \sigma_1\\ \sigma_2\\ \tau_{12} \end{bmatrix}, \qquad \boldsymbol{\varepsilon}_{12}= \begin{bmatrix} \varepsilon_1\\ \varepsilon_2\\ \gamma_{12} \end{bmatrix}. \]

Here each symbol has a specific physical meaning:

Important convention. In CLT, \(\gamma_{12}\) is the engineering shear strain, not the tensor shear strain.

Therefore the shear stiffness appears directly as \(Q_{66}=G_{12}\).

Under plane stress,

\[ \boldsymbol{\sigma}_{12}=\mathbf{Q}\boldsymbol{\varepsilon}_{12}, \qquad \mathbf{Q}= \begin{bmatrix} Q_{11} & Q_{12} & 0\\ Q_{12} & Q_{22} & 0\\ 0 & 0 & Q_{66} \end{bmatrix}. \]

The zero terms in \(\mathbf{Q}\) must be interpreted carefully. They do not mean that shear is absent. They mean only that, in the material principal axes of an orthotropic lamina, the normal response and the in-plane shear response are uncoupled.

Therefore the zeros in the reduced stiffness matrix reflect a special coordinate choice, not the disappearance of shear physics.

\[ \sigma_3=\tau_{13}=\tau_{23}=0. \]

The plane-stress assumption means that the lamina is thin enough that out-of-plane stress components are neglected. This is the standard reduction used in Classical Lamination Theory.

What \(\mathbf{Q}\) represents.

\(\mathbf{Q}\) is the reduced stiffness matrix of one lamina in its own material frame. “Reduced” means reduced from a full three-dimensional orthotropic constitutive law to the two-dimensional plane-stress form used in CLT.

The reciprocal Poisson ratio is

\[ \nu_{21}=\nu_{12}\frac{E_2}{E_1}. \]

This relation follows from reciprocity and constitutive symmetry in linear elasticity. It is required for a physically consistent orthotropic lamina model.

The reduced stiffness terms are

\[ Q_{11}=\frac{E_1}{1-\nu_{12}\nu_{21}},\qquad Q_{22}=\frac{E_2}{1-\nu_{12}\nu_{21}},\qquad Q_{12}=\frac{\nu_{12}E_2}{1-\nu_{12}\nu_{21}},\qquad Q_{66}=G_{12}. \]

Their physical meaning is:

Material constants vs. constitutive coefficients.

\(E_1\), \(E_2\), \(G_{12}\), and \(\nu_{12}\) are the engineering input data. \(Q_{11}\), \(Q_{22}\), \(Q_{12}\), and \(Q_{66}\) are the constitutive coefficients derived from them. In CLT calculations, the stiffness matrix \(\mathbf{Q}\) is what is actually used.

Numerical and physical illustration: one lamina, one material, four strain states.

The MATLAB program below considers a single lamina with one fixed set of material properties \((E_1,E_2,\nu_{12},G_{12})\). The material is the same in all four cases. Only the imposed in-plane strain state is changed.

Therefore the four plots do not represent four different materials. They represent the same orthotropic lamina subjected to four different strain inputs:

  1. longitudinal tension: \(\varepsilon_1 \neq 0\)
  2. transverse tension: \(\varepsilon_2 \neq 0\)
  3. pure shear: \(\gamma_{12} \neq 0\)
  4. biaxial tension: \(\varepsilon_1 \neq 0\) and \(\varepsilon_2 \neq 0\)

These different deformations can be predicted and calculated directly from the constitutive discussion in this section. Once the strain vector \(\boldsymbol{\varepsilon}_{12}\) is prescribed, the stress vector \(\boldsymbol{\sigma}_{12}\) follows from

\[ \boldsymbol{\sigma}_{12}=\mathbf{Q}\boldsymbol{\varepsilon}_{12}. \]

The same constitutive law also explains why each imposed strain pattern produces a different geometric deformation of the same square lamina.

4.1 MATLAB illustration for one lamina under four prescribed strain states

% Lamina Analysis: Material Properties and Deformation Visualization
clear; clc; clf;

% ==========================================================
% PART 1: Define Lamina Material Properties
% ==========================================================
% Example: Carbon/Epoxy T300 lamina properties
E1  = 181e9;    % Young's modulus in fiber direction (Pa)
E2  = 10.3e9;   % Young's modulus in transverse direction (Pa)
v12 = 0.28;     % Major Poisson's ratio
G12 = 7.17e9;   % In-plane shear modulus (Pa)

% Calculate minor Poisson's ratio using reciprocity
v21 = v12 * E2 / E1;

% Compute reduced stiffness matrix [Q] components
denom = 1 - v12 * v21;
Q11 = E1 / denom;
Q12 = v12 * E2 / denom;
Q22 = E2 / denom;
Q66 = G12;

% ==========================================================
% PART 2: Set Plot Parameters
% ==========================================================
ax_fs    = 14;   % Axis label font size
title_fs = 16;   % Subplot title font size
label_fs = 20;   % Right-side strain label font size
text_fs  = 11;   % General annotation font size

% Original unit square in the 1-2 material coordinate system
x_orig = [0 1 1 0 0];
y_orig = [0 0 1 1 0];

% Define four prescribed strain cases: [e1; e2; g12]
cases = {
    [0.2; 0;   0  ],   % Case 1: Longitudinal tension
    [0;   0.2; 0  ],   % Case 2: Transverse tension
    [0;   0;   0.4],   % Case 3: Pure shear
    [0.2; 0.2; 0  ]    % Case 4: Biaxial tension
};

titles = {
    'Case 1: Longitudinal Tension', ...
    'Case 2: Transverse Tension', ...
    'Case 3: Pure Shear', ...
    'Case 4: Biaxial Tension'
};

figure('Color', 'w', 'Name', 'Lamina Analysis - Final Version');

% ==========================================================
% PART 3: Loop Through Cases and Plot Deformation
% ==========================================================
for i = 1:4
    subplot(2, 2, i);
    hold on;
    
    % Current prescribed strain components
    e1  = cases{i}(1);
    e2  = cases{i}(2);
    g12 = cases{i}(3);
    
    % ------------------------------------------------------
    % Effective strains used only for visualization
    % Poisson contraction is included in uniaxial cases
    % ------------------------------------------------------
    e1_plot = e1;
    e2_plot = e2;
    g12_plot = g12;
    
    if i == 1
        % Case 1: Longitudinal tension causes transverse contraction
        e2_plot = -v12 * e1;
    elseif i == 2
        % Case 2: Transverse tension causes longitudinal contraction
        e1_plot = -v21 * e2;
    end
    
    % ------------------------------------------------------
    % Stress calculation in principal material coordinates
    % sigma = Q * epsilon
    % Note: Q16 = Q26 = 0 in the principal material axes
    % ------------------------------------------------------
    s1    = Q11 * e1 + Q12 * e2;
    s2    = Q12 * e1 + Q22 * e2;
    tau12 = Q66 * g12;
    
    % ------------------------------------------------------
    % Deformation mapping
    % x' = (1 + e1)x + gamma12*y
    % y' = (1 + e2)y
    % ------------------------------------------------------
    x_new = (1 + e1_plot) * x_orig + g12_plot * y_orig;
    y_new = (1 + e2_plot) * y_orig;
    
    % Plot original and deformed shapes
    plot(x_orig, y_orig, '--', 'LineWidth', 1.5, 'Color', [0.6 0.6 0.6]);
    plot(x_new, y_new, 'r-s', 'LineWidth', 2.5, 'MarkerFaceColor', 'r', 'MarkerSize', 7);
    
    % Plot fiber direction arrow
    quiver(0.1, 0.5, 0.45, 0, 0, 'b', 'LineWidth', 2.2, 'MaxHeadSize', 0.5);
    text(0.14, 0.58, 'Fiber', 'Color', 'b', 'FontWeight', 'bold', 'FontSize', 12);
    
    % ------------------------------------------------------
    % Add loading/boundary-condition annotations
    % ------------------------------------------------------
    switch i
        case 1
            % Longitudinal tension: pull left and right
            quiver(-0.08, 0.5, -0.18, 0, 0, 'Color', [0.85 0 0], 'LineWidth', 2, 'MaxHeadSize', 0.8);
            quiver(1.22, 0.5,  0.18, 0, 0, 'Color', [0.85 0 0], 'LineWidth', 2, 'MaxHeadSize', 0.8);
            text(-0.18, 0.63, '\sigma_1', 'Color', [0.85 0 0], 'FontSize', text_fs, 'FontWeight', 'bold');
            text(1.24, 0.63, '\sigma_1', 'Color', [0.85 0 0], 'FontSize', text_fs, 'FontWeight', 'bold');
            
        case 2
            % Transverse tension: pull top and bottom
            quiver(0.5, -0.08, 0, -0.18, 0, 'Color', [0.85 0 0], 'LineWidth', 2, 'MaxHeadSize', 0.8);
            quiver(0.5, 1.22, 0,  0.18, 0, 'Color', [0.85 0 0], 'LineWidth', 2, 'MaxHeadSize', 0.8);
            text(0.58, -0.18, '\sigma_2', 'Color', [0.85 0 0], 'FontSize', text_fs, 'FontWeight', 'bold');
            text(0.58,  1.33, '\sigma_2', 'Color', [0.85 0 0], 'FontSize', text_fs, 'FontWeight', 'bold');
            
        case 3
            % Pure shear: bottom fixed, top edge sheared to the right
            quiver(0.35, 1.10, 0.35, 0, 0, 'Color', [0.85 0 0], 'LineWidth', 2, 'MaxHeadSize', 0.8);
            text(0.74, 1.16, '\tau_{12}', 'Color', [0.85 0 0], 'FontSize', text_fs, 'FontWeight', 'bold');
            text(0.18, -0.18, 'Bottom edge fixed', 'Color', [0.2 0.2 0.2], 'FontSize', text_fs, 'FontWeight', 'bold');
            
        case 4
            % Biaxial tension: pull in both 1- and 2-directions
            quiver(-0.08, 0.6, -0.18, 0, 0, 'Color', [0.85 0 0], 'LineWidth', 2, 'MaxHeadSize', 0.8);
            quiver(1.22, 0.6,  0.18, 0, 0, 'Color', [0.85 0 0], 'LineWidth', 2, 'MaxHeadSize', 0.8);
            quiver(0.6, -0.08, 0, -0.18, 0, 'Color', [0.85 0 0], 'LineWidth', 2, 'MaxHeadSize', 0.8);
            quiver(0.6, 1.22, 0,  0.18, 0, 'Color', [0.85 0 0], 'LineWidth', 2, 'MaxHeadSize', 0.8);
            text(-0.18, 0.73, '\sigma_1', 'Color', [0.85 0 0], 'FontSize', text_fs, 'FontWeight', 'bold');
            text(1.24, 0.73, '\sigma_1', 'Color', [0.85 0 0], 'FontSize', text_fs, 'FontWeight', 'bold');
            text(0.68, -0.18, '\sigma_2', 'Color', [0.85 0 0], 'FontSize', text_fs, 'FontWeight', 'bold');
            text(0.68,  1.33, '\sigma_2', 'Color', [0.85 0 0], 'FontSize', text_fs, 'FontWeight', 'bold');
    end
    
    % Axes formatting
    axis([-0.25 2.5 -0.25 1.65]);
    axis equal;
    grid on;
    box on;
    title(titles{i}, 'FontSize', title_fs, 'FontWeight', 'bold');
    xlabel('Direction 1 (Fiber)', 'FontSize', ax_fs);
    ylabel('Direction 2 (Transverse)', 'FontSize', ax_fs);
    
    % ------------------------------------------------------
    % Display prescribed strain labels on the right side
    % ------------------------------------------------------
    label_x = 1.72;
    
    if e1 > 0
        text(label_x, 1.15, ['\epsilon_1 = ' num2str(e1)], ...
            'Color', 'r', 'FontSize', label_fs, 'FontWeight', 'bold');
    end
    
    if e2 > 0
        if e1 > 0
            y_pos = 0.85;
        else
            y_pos = 1.15;
        end
        text(label_x, y_pos, ['\epsilon_2 = ' num2str(e2)], ...
            'Color', 'r', 'FontSize', label_fs, 'FontWeight', 'bold');
    end
    
    if g12 > 0
        text(label_x, 0.60, ['\gamma_{12} = ' num2str(g12)], ...
            'Color', 'r', 'FontSize', label_fs, 'FontWeight', 'bold');
    end
end

% Optional overall title
sgtitle('In-Plane Deformation Modes of an Orthotropic Lamina', ...
        'FontSize', 18, 'FontWeight', 'bold');
One lamina under four prescribed in-plane strain states
One lamina, one material, four prescribed in-plane strain states. The same orthotropic single-ply lamina is shown under longitudinal tension, transverse tension, pure shear, and biaxial tension. The different deformed shapes arise from different prescribed strain vectors \(\boldsymbol{\varepsilon}_{12}\), while the material itself remains the same. Section 4 provides the constitutive relation used to predict and calculate the corresponding stresses.

How to read the four panels.

  • Case 1: \(\varepsilon_1 \neq 0\), so the lamina stretches mainly in the fiber direction.
  • Case 2: \(\varepsilon_2 \neq 0\), so the lamina stretches mainly transverse to the fiber direction.
  • Case 3: \(\gamma_{12} \neq 0\), so the square distorts into a parallelogram under pure in-plane shear.
  • Case 4: both \(\varepsilon_1\) and \(\varepsilon_2\) are nonzero, so the lamina expands biaxially.

This is exactly the point of the lamina constitutive law: the same material responds differently depending on the imposed strain state, and the corresponding stress state can be computed from \(\mathbf{Q}\).

What this section accomplishes.

Up to this point, we only know the constitutive behavior of one lamina in its own material axes. A real laminate, however, is built in the global laminate axes \((x,y)\), and different plies are rotated by different angles. Therefore the next necessary step is to express each ply stiffness in the laminate coordinate system.

5. Coordinate transformation and rotated ply stiffness

A laminate is defined in the global axes \((x,y)\), not in the material axes of every ply. Therefore the ply stiffness must be transformed into the laminate frame.

\[ \boldsymbol{\sigma}_{xy}=\bar{\mathbf{Q}}\boldsymbol{\varepsilon}_{xy}, \]

where \(\bar{\mathbf{Q}}\) is the transformed reduced stiffness matrix of a ply oriented at angle \(\theta\). This is the step that makes the stacking sequence mechanically meaningful.

In other words, the zero coupling terms in the material-axis matrix \(\mathbf{Q}\) do not remain zero in general; they remain zero only for special alignment cases such as \(0^\circ\) and \(90^\circ\) plies.

In the laminate axes, define

\[ \boldsymbol{\sigma}_{xy}= \begin{bmatrix} \sigma_x\\ \sigma_y\\ \tau_{xy} \end{bmatrix}, \qquad \boldsymbol{\varepsilon}_{xy}= \begin{bmatrix} \varepsilon_x\\ \varepsilon_y\\ \gamma_{xy} \end{bmatrix}. \]

In general, the transformed reduced stiffness has the form

\[ \bar{\mathbf{Q}}= \begin{bmatrix} \bar{Q}_{11} & \bar{Q}_{12} & \bar{Q}_{16}\\ \bar{Q}_{12} & \bar{Q}_{22} & \bar{Q}_{26}\\ \bar{Q}_{16} & \bar{Q}_{26} & \bar{Q}_{66} \end{bmatrix}. \]

This matrix is the general plane-stress constitutive form of a rotated orthotropic lamina written in the laminate axes. The matrix form itself is valid for any ply orientation \(\theta\), but the values of its entries depend on both the lamina properties and the orientation angle.

To compute \(\bar{\mathbf{Q}}\), we start from the reduced stiffness matrix \(\mathbf{Q}\) in the material axes \((1,2)\) and transform it into the laminate axes \((x,y)\). Define

\[ m=\cos\theta,\qquad n=\sin\theta. \]

The transformed reduced stiffness components are

\[ \bar{Q}_{11}=Q_{11}m^4+2\left(Q_{12}+2Q_{66}\right)m^2n^2+Q_{22}n^4, \] \[ \bar{Q}_{22}=Q_{11}n^4+2\left(Q_{12}+2Q_{66}\right)m^2n^2+Q_{22}m^4, \] \[ \bar{Q}_{12}=\left(Q_{11}+Q_{22}-4Q_{66}\right)m^2n^2+Q_{12}\left(m^4+n^4\right), \] \[ \bar{Q}_{16}=\left(Q_{11}-Q_{12}-2Q_{66}\right)m^3n-\left(Q_{22}-Q_{12}-2Q_{66}\right)mn^3, \] \[ \bar{Q}_{26}=\left(Q_{11}-Q_{12}-2Q_{66}\right)mn^3-\left(Q_{22}-Q_{12}-2Q_{66}\right)m^3n, \] \[ \bar{Q}_{66}=\left(Q_{11}+Q_{22}-2Q_{12}-2Q_{66}\right)m^2n^2+Q_{66}\left(m^4+n^4\right). \]

How the transformation should be read.

The material-axis matrix \(\mathbf{Q}\) is first built from the lamina engineering constants \(E_1\), \(E_2\), \(G_{12}\), and \(\nu_{12}\). The trigonometric factors \(m=\cos\theta\) and \(n=\sin\theta\) then rotate that local stiffness into the laminate axes. In other words, the laminate does not change the material; it changes the coordinate frame in which the same ply stiffness is written.

This is why \(\bar{Q}_{16}\) and \(\bar{Q}_{26}\) appear only after rotation. They are not new material constants. They are transformed stiffness terms produced by expressing the orthotropic ply in a rotated global frame.

When this matrix form is valid.

  1. The ply is modeled as an orthotropic lamina.
  2. The constitutive law is written in the plane-stress reduced form of CLT.
  3. The stresses and strains are expressed in the laminate axes \((x,y)\), not in the local material axes \((1,2)\).
  4. The ply may be rotated by any angle \(\theta\) relative to the laminate axes.

Therefore, the 3×3 structure of \(\bar{\mathbf{Q}}\) is the correct general form for a rotated ply in CLT. What changes from case to case is not the matrix size or structure, but whether the coupling terms \(\bar{Q}_{16}\) and \(\bar{Q}_{26}\) vanish or remain nonzero.

Key physical point.

In the material axes, the orthotropic lamina had no normal–shear coupling terms in the constitutive matrix. After rotation, such coupling generally appears through \(\bar{Q}_{16}\) and \(\bar{Q}_{26}\). This is one of the main reasons orientation matters so strongly in laminate mechanics.

The transformed matrix \(\bar{\mathbf{Q}}\) is what is actually assembled through the thickness to form the laminate stiffness matrices \(\mathbf{A}\), \(\mathbf{B}\), and \(\mathbf{D}\). Therefore every ply contributes to the global laminate response through its own rotated stiffness, not through the unrotated \(\mathbf{Q}\) matrix.

How the general form should be interpreted.

The matrix

\[ \bar{\mathbf{Q}}= \begin{bmatrix} \bar{Q}_{11} & \bar{Q}_{12} & \bar{Q}_{16}\\ \bar{Q}_{12} & \bar{Q}_{22} & \bar{Q}_{26}\\ \bar{Q}_{16} & \bar{Q}_{26} & \bar{Q}_{66} \end{bmatrix} \]

should be understood as a single general constitutive template for all rotated plies under plane stress. It does not imply that \(\bar{Q}_{16}\) and \(\bar{Q}_{26}\) are always nonzero. Instead, those two terms are angle-dependent and identify whether normal and shear response are coupled in the laminate frame.

Special and general cases.

Case A: \(0^\circ\) ply.

If the ply is aligned with the laminate axes, then the laminate axes and the material principal axes coincide. In that case, \(\bar{\mathbf{Q}}=\mathbf{Q}\), and the coupling terms remain zero: \(\bar{Q}_{16}=\bar{Q}_{26}=0\).

Case B: \(90^\circ\) ply.

A \(90^\circ\) ply is still aligned with a material symmetry direction, even though the fiber and transverse directions are exchanged relative to the laminate axes. Because the laminate axes are still along principal material directions, the normal–shear coupling terms again vanish: \(\bar{Q}_{16}=\bar{Q}_{26}=0\).

Case C: off-axis ply such as \(45^\circ\).

For an off-axis ply, the laminate axes no longer coincide with the material principal axes. Then normal and shear response generally become coupled in the transformed constitutive law, and \(\bar{Q}_{16}\) and \(\bar{Q}_{26}\) are generally nonzero. This is the usual case for rotated plies in a laminate.

Most important takeaway.

The local material-axis matrix \(\mathbf{Q}\) and the global laminate-axis matrix \(\bar{\mathbf{Q}}\) should never be mixed. The matrix \(\mathbf{Q}\) is valid in the local principal material axes of a ply. The matrix \(\bar{\mathbf{Q}}\) is the rotated form used when that ply is embedded in a laminate and described in global coordinates.

For \(0^\circ\) and \(90^\circ\) plies, the global matrix still has the same 3×3 form, but \(\bar{Q}_{16}=\bar{Q}_{26}=0\). For a general off-axis ply, those coupling terms are generally nonzero.

6. Laminate A, B, and D matrices

Let \(z_0,z_1,\ldots,z_n\) denote the through-thickness interface coordinates. Then the laminate stiffness matrices are

\[ \mathbf{A}=\sum_{k=1}^{n}\bar{\mathbf{Q}}^{(k)}(z_k-z_{k-1}), \] \[ \mathbf{B}=\frac{1}{2}\sum_{k=1}^{n}\bar{\mathbf{Q}}^{(k)}(z_k^2-z_{k-1}^2), \] \[ \mathbf{D}=\frac{1}{3}\sum_{k=1}^{n}\bar{\mathbf{Q}}^{(k)}(z_k^3-z_{k-1}^3). \]

If the laminate is symmetric about the mid-plane, \(\mathbf{B}=0\).

Physical interpretation of the \(\mathbf{A}\), \(\mathbf{B}\), and \(\mathbf{D}\) matrices.

This set of equations is the core of Classical Lamination Theory (CLT). It explains how the stiffness of the full laminate is built by accumulating the stiffness contribution of each individual ply through the thickness direction. In essence, the laminate stiffness is obtained by a discrete summation of ply stiffness contributions over the thickness coordinate \(z\).

1. \(\mathbf{A}\) matrix: extensional stiffness

The \(\mathbf{A}\) matrix represents the in-plane or extensional stiffness of the laminate. Its meaning is the most direct: the resistance of the laminate to in-plane stretching depends on how stiff each ply is and how much thickness that ply occupies.

Physically, every ply contributes to extensional stiffness regardless of whether it is located near the top, bottom, or middle of the laminate. As long as the ply exists, it adds to \(\mathbf{A}\).

\[ \mathbf{A}=\int \bar{\mathbf{Q}}\,dz \;\approx\; \sum_k \left(\bar{\mathbf{Q}}^{(k)}\Delta z_k\right). \]

2. \(\mathbf{B}\) matrix: coupling stiffness

The \(\mathbf{B}\) matrix governs membrane–bending coupling. This is the matrix that determines whether in-plane loading can induce bending or whether bending can induce in-plane strain. Its physical meaning is tied directly to the symmetry of stiffness about the laminate mid-plane \((z=0)\).

In practical terms, the contribution of a ply to \(\mathbf{B}\) depends not only on its stiffness, but also on its distance from the mid-plane. If the laminate is perfectly symmetric, the positive and negative \(z\)-side contributions cancel, giving \(\mathbf{B}=0\). If the laminate is not symmetric, the cancellation is incomplete, and \(\mathbf{B}\neq 0\), which leads to coupling and possible warpage.

For example, if the upper half and lower half of the laminate do not mirror each other in orientation and stiffness, their weighted contributions about the mid-plane will not balance. That is the origin of extension–bending coupling.

\[ \mathbf{B}=\int \bar{\mathbf{Q}}\,z\,dz \;\approx\; \sum_k \left(\bar{\mathbf{Q}}^{(k)}\Delta z_k\,\bar z_k\right), \]

where \(\bar z_k\) is the distance from the mid-plane to the center of the \(k\)-th ply.

3. \(\mathbf{D}\) matrix: bending stiffness

The \(\mathbf{D}\) matrix measures the resistance of the laminate to bending. The key idea is the same as in elementary beam theory: material located farther from the neutral axis contributes much more strongly to bending resistance than material placed near the middle.

This is closely related to the parallel-axis principle. The dominant term scales with the square of the distance from the mid-plane, which means that plies near the outer surfaces are much more effective in increasing bending stiffness. This is why placing stiff plies near the top and bottom surfaces is highly effective when bending resistance is important.

In the discrete form, one contribution comes from the ply position relative to the mid-plane, and another smaller contribution comes from the bending inertia of the ply about its own centroid.

\[ \mathbf{D}=\int \bar{\mathbf{Q}}\,z^2\,dz \;\approx\; \sum_k \left[ \bar{\mathbf{Q}}^{(k)} \left( \Delta z_k\,\bar z_k^2+\frac{\Delta z_k^3}{12} \right) \right]. \]

Summary of the accumulation logic

  • \(\mathbf{A}\): depends mainly on how much stiffness material is present through the thickness.
  • \(\mathbf{B}\): depends on whether the stiffness distribution is balanced or unbalanced about the mid-plane.
  • \(\mathbf{D}\): depends strongly on whether stiff material is placed near the middle or near the outer surfaces.

In that sense, the ABD formulation acts like a mechanical bookkeeping system for the laminate: \(\mathbf{A}\) tracks extensional stiffness, \(\mathbf{B}\) tracks coupling caused by asymmetry, and \(\mathbf{D}\) tracks bending resistance controlled by through-thickness stiffness placement.

The thickness coordinates \(z_k\) are measured through the laminate thickness, usually from the laminate mid-plane. Therefore \(z=0\) is the mid-surface, negative \(z\) lies below it, and positive \(z\) lies above it.

Units matter.

  • \(\bar{\mathbf{Q}}\): units Pa = N/m\(^2\)
  • \(\mathbf{A}\): units N/m
  • \(\mathbf{B}\): units N
  • \(\mathbf{D}\): units N\(\cdot\)m

7. Load resultants, mid-plane strain, and curvature

Once the laminate stiffness matrices \(\mathbf{A}\), \(\mathbf{B}\), and \(\mathbf{D}\) are known, the next step in Classical Lamination Theory is to connect the applied laminate loads to the global deformation of the plate. This is done through the ABD constitutive system:

\[ \begin{bmatrix} \mathbf{N}\\ \mathbf{M} \end{bmatrix} = \begin{bmatrix} \mathbf{A} & \mathbf{B}\\ \mathbf{B} & \mathbf{D} \end{bmatrix} \begin{bmatrix} \boldsymbol{\varepsilon}_0\\ \boldsymbol{\kappa} \end{bmatrix}. \]

This is the main laminate-level equilibrium relation solved in CLT. It tells us how in-plane force resultants and bending/twisting moment resultants generate mid-plane strains and curvatures.

Define

\[ \mathbf{N}= \begin{bmatrix} N_x\\ N_y\\ N_{xy} \end{bmatrix}, \qquad \mathbf{M}= \begin{bmatrix} M_x\\ M_y\\ M_{xy} \end{bmatrix}, \] \[ \boldsymbol{\varepsilon}_0= \begin{bmatrix} \varepsilon_x^0\\ \varepsilon_y^0\\ \gamma_{xy}^0 \end{bmatrix}, \qquad \boldsymbol{\kappa}= \begin{bmatrix} \kappa_x\\ \kappa_y\\ \kappa_{xy} \end{bmatrix}. \]

Why this section is useful.

This section is where the theory becomes structurally predictive. Sections 4–6 define ply stiffness, rotated ply stiffness, and laminate stiffness. Section 7 uses those quantities to answer the practical question: under a given mechanical loading, how will the laminate deform?

In other words, this section gives the direct bridge from loading to global laminate response. Once \(\boldsymbol{\varepsilon}_0\) and \(\boldsymbol{\kappa}\) are known, all ply strains and stresses can be recovered later through the thickness.

What this section can predict.

  • Whether the laminate mainly stretches, bends, twists, or does a combination of them
  • The magnitude of mid-plane in-plane deformation under membrane loading
  • The laminate curvature under bending or twisting moments
  • Whether an in-plane load can induce bending due to \(\mathbf{B}\neq 0\)
  • Whether a moment can induce in-plane strain due to membrane–bending coupling
  • Whether an asymmetric laminate is likely to warp after loading or due to internal mismatch effects

This makes Section 7 essential for both deformation prediction and later stress recovery.

Physical interpretation of coupling.

If \(\mathbf{B}=0\), stretching and bending are decoupled. Then an in-plane load produces mainly mid-plane strain, while a pure moment produces mainly curvature. If \(\mathbf{B}\neq 0\), they interact. In that case, a purely in-plane load can induce bending, and a pure moment can induce mid-plane strain.

This is exactly why laminate symmetry matters so much in design. A symmetric laminate is much easier to control, while an asymmetric laminate can deform in less intuitive ways.

Comparison of warpage response in asymmetric and symmetric balanced laminates
Comparison of asymmetric and symmetric balanced laminates under coupling-sensitive conditions. The left panel shows an asymmetric balanced laminate with a nonzero coupling term \(B_{16}\), which produces out-of-plane deformation. The right panel shows a symmetric balanced laminate for which the coupling term vanishes, leading to a flat response. This illustrates why the \(\mathbf{B}\) matrix is so important in laminate deformation prediction.

Why the figure matters.

The figure provides a direct geometric interpretation of what the ABD system predicts. In the asymmetric case, the nonzero coupling term means that membrane-type loading is no longer isolated from bending-type response, so the plate can warp out of plane. In the symmetric case, the coupling term vanishes, and the same type of unintended out-of-plane deformation is suppressed.

This is a very useful teaching example because students often understand the algebra of \(\mathbf{B}\) only after they see its geometric consequence: nonzero coupling can produce warpage.

Typical engineering applications that need this section.

  • Aerospace composite panels: prediction of extension, bending, and twist under combined service loads
  • 3DIC and electronic substrates: warpage control caused by asymmetric layer stacking and stiffness mismatch
  • Printed circuit boards and package laminates: curvature and flatness prediction during assembly and service
  • Wind turbine blades and sporting structures: tailoring bending–twist coupling for performance
  • Composite automotive panels: stiffness design and deformation control under distributed loading
  • Satellite and space structures: dimensional stability and coupling suppression in thin laminated components

In all of these applications, engineers do not only want the laminate to be strong. They also want it to deform in a controlled and predictable way. That is exactly what Section 7 helps quantify.

Most important takeaway.

Section 7 is the point where laminate mechanics becomes directly usable for design. It predicts how a laminate responds globally to applied force and moment resultants, and it identifies whether symmetry, asymmetry, or coupling will lead to desirable behavior or unwanted warpage.

8. Ply stress recovery

The in-plane strain at thickness coordinate \(z\) is

\[ \boldsymbol{\varepsilon}_{xy}(z)=\boldsymbol{\varepsilon}_0+z\boldsymbol{\kappa}. \]

The global ply stress is

\[ \boldsymbol{\sigma}_{xy}^{(k)}(z)=\bar{\mathbf{Q}}^{(k)}\boldsymbol{\varepsilon}_{xy}(z). \]

For failure criteria, this stress must often be transformed back into the material frame to obtain \(\sigma_1\), \(\sigma_2\), and \(\tau_{12}\).

Practical point. Stress is often evaluated at the top and bottom of each ply, because bending-driven peak values often occur at these surfaces.

9. Strength and first-ply failure

A simple failure indicator is the Tsai–Hill failure index

\[ FI= \left(\frac{\sigma_1}{X}\right)^2 -\frac{\sigma_1\sigma_2}{X^2} +\left(\frac{\sigma_2}{Y}\right)^2 +\left(\frac{\tau_{12}}{S}\right)^2. \]

Here \(X\) and \(Y\) are selected from tensile or compressive strength values according to the sign of \(\sigma_1\) and \(\sigma_2\).

Interpretation.

  • \(FI<1\): below the criterion surface
  • \(FI=1\): reaches the criterion surface
  • \(FI>1\): predicted failure by the criterion

10. Worked laminate example

Use representative carbon/epoxy lamina data:

QuantityValue
\(E_1\)135 GPa
\(E_2\)10 GPa
\(G_{12}\)5 GPa
\(\nu_{12}\)0.30
\(X_t\)1500 MPa
\(X_c\)1200 MPa
\(Y_t\)40 MPa
\(Y_c\)150 MPa
\(S\)70 MPa

Laminate:

\[ [0/45/-45/90]_s,\qquad t_{ply}=0.125~\text{mm} \]

Applied load:

\[ \mathbf{N}= \begin{bmatrix} 1000\\ 0\\ 0 \end{bmatrix}\text{N/m},\qquad \mathbf{M}= \begin{bmatrix} 0\\ 0\\ 0 \end{bmatrix}. \]

Why this example is useful.

Because the laminate is symmetric, a correct implementation should produce \(\mathbf{B}\approx 0\), and under pure membrane loading the curvatures should be negligible.

11. MATLAB implementation

11.1 Main script

E1  = 135e9;  E2  = 10e9;  G12 = 5e9;  v12 = 0.30;
Xt = 1500e6;  Xc = 1200e6; Yt = 40e6;  Yc = 150e6; S = 70e6;

angles = [0 45 -45 90 90 -45 45 0];
tply   = 0.125e-3;
NM     = [1000;0;0;0;0;0];

[ABD, z, Qbar_all] = build_laminate_ABD(E1,E2,G12,v12,angles,tply);

ek    = ABD \ NM;
eps0  = ek(1:3);
kappa = ek(4:6);

11.2 Function: build laminate ABD

function [ABD, z, Qbar_all] = build_laminate_ABD(E1,E2,G12,v12,angles,tply)
v21 = v12 * E2 / E1;
den = 1 - v12*v21;

Q11 = E1 / den;   Q22 = E2 / den;
Q12 = v12 * E2 / den;   Q66 = G12;

nply = length(angles);
h = nply * tply;
z = linspace(-h/2, h/2, nply+1);

A = zeros(3,3); B = zeros(3,3); D = zeros(3,3);
Qbar_all = zeros(3,3,nply);

for k = 1:nply
    theta = angles(k);
    m = cosd(theta); n = sind(theta);

    Qbar11 = Q11*m^4 + 2*(Q12+2*Q66)*m^2*n^2 + Q22*n^4;
    Qbar22 = Q11*n^4 + 2*(Q12+2*Q66)*m^2*n^2 + Q22*m^4;
    Qbar12 = (Q11+Q22-4*Q66)*m^2*n^2 + Q12*(m^4+n^4);
    Qbar16 = (Q11-Q12-2*Q66)*m^3*n - (Q22-Q12-2*Q66)*m*n^3;
    Qbar26 = (Q11-Q12-2*Q66)*m*n^3 - (Q22-Q12-2*Q66)*m^3*n;
    Qbar66 = (Q11+Q22-2*Q12-2*Q66)*m^2*n^2 + Q66*(m^4+n^4);

    Qbar = [Qbar11 Qbar12 Qbar16;
            Qbar12 Qbar22 Qbar26;
            Qbar16 Qbar26 Qbar66];

    Qbar_all(:,:,k) = Qbar;

    zk1 = z(k); zk2 = z(k+1);
    A = A + Qbar * (zk2-zk1);
    B = B + 0.5 * Qbar * (zk2^2-zk1^2);
    D = D + (1/3) * Qbar * (zk2^3-zk1^3);
end

ABD = [A B; B D];
end

11.3 Function: transform global stress to ply axes

function sigma_12 = stress_xy_to_12(sigma_xy, theta)
m = cosd(theta); n = sind(theta);
sx = sigma_xy(1); sy = sigma_xy(2); txy = sigma_xy(3);

s1  = m^2*sx + n^2*sy + 2*m*n*txy;
s2  = n^2*sx + m^2*sy - 2*m*n*txy;
t12 = -m*n*sx + m*n*sy + (m^2-n^2)*txy;

sigma_12 = [s1; s2; t12];
end

12. Example numerical results and how to read them

For the symmetric laminate \([0/45/-45/90]_s\) under pure in-plane loading, the first structural checks are:

  1. Is \(\mathbf{B}\) numerically close to zero?
  2. Are the curvatures \(\boldsymbol{\kappa}\) negligible under pure membrane loading?
  3. Which plies carry the highest longitudinal stress \(\sigma_1\)?
  4. Which plies carry the highest transverse stress \(\sigma_2\) or shear stress \(\tau_{12}\)?
  5. Which ply gives the largest failure index \(FI\)?

Typical interpretation. In a symmetric laminate under pure in-plane loading, one expects \(\mathbf{B}\approx 0\), and therefore the induced curvature should also be very small. If the computed \(\boldsymbol{\kappa}\) is unexpectedly large, the first things to check are the stacking sequence, the sign convention, and the ABD assembly.

The stress distribution should also be interpreted physically:

This is exactly why the stress must be transformed back to the ply axes before any strength statement is made.

13. Recommended plots and outputs

A good teaching script should not end with raw console output alone. The following outputs are particularly useful:

Recommended plot logic:

  1. Plot ply number on the horizontal axis.
  2. Plot \(\sigma_1\), \(\sigma_2\), or \(\tau_{12}\) on the vertical axis.
  3. Use a separate plot for the failure index \(FI\).
  4. Mark the critical ply directly in the figure caption or title.

These plots make the code pedagogically useful, because they connect laminate theory to interpretable structural behavior.

Back to module home