Composite Plate Bending Analysis With Matlab Code May 2026

If the fibers are oriented at an angle $\theta$ relative to the plate axis ($x-y$), we transform the stiffness matrix: $$ [\barQ] = [T]^-1 [Q] [T]^-T $$ (This is handled via transformation matrices involving $\sin\theta$ and $\cos\theta$).

%% Composite Plate Bending Analysis using FEM (FSDT)
% Analysis of a simply supported rectangular composite plate
clc; clear; close all;
%% 1. Input Parameters
% Geometry
Lx = 1.0;       % Length in x (m)
Ly = 1.0;       % Length in y (m)
thick = 0.01;   % Total thickness (m)
% Mesh Density
nx = 20;        % Elements in x
ny = 20;        % Elements in y
% Material Properties (Orthotropic - Typical Carbon/Epoxy)
E1 = 140e9;     % Longitudinal Modulus (Pa)
E2 = 10e9;      % Transverse Modulus (Pa)
G12 = 5e9;      % In-plane Shear Modulus (Pa)
G23 = 3.5e9;    % Out-of-plane Shear Modulus (Pa) (approx)
nu12 = 0.3;     % Major Poisson's ratio
% Layup Definition
% Format: [Angle (deg), Thickness (m)]
% Symmetric 4-layer layup [0/90]_s
layup = [
0,   thick/4;
90,  thick/4;
90,  thick/4;
0,   thick/4
];
% Loading
q = 1000; % Uniform transverse load (N/m^2)
%% 2. Calculate ABD Matrix
% Uses Classical Lamination Theory (CLT)
ABD = calculate_ABD(layup, E1, E2, G12, nu12, G23);
A = ABD.A;
B = ABD.B;
D = ABD.D;
Hs = ABD.Hs; % Shear stiffness matrix
% Extract constitutive matrices for FEM
% For pure bending analysis, we focus on D and H.
Db = D;
Ds = Hs;
%% 3. Mesh Generation
[node, element] = create_mesh(Lx, Ly, nx, ny);
n_node = size(node, 1);
n_elem = size(element, 1);
n_dof = 5 * n_node; % 5 DOFs per node (u,v,w,thx,thy)
% Plot Mesh
figure(1);
patch('Faces', element, 'Vertices', node, 'FaceColor', 'c', 'EdgeColor', 'k');
axis equal; title('Mesh'); xlabel('x'); ylabel('y');
%% 4. Global Stiffness Matrix Assembly
K = sparse(n_dof, n_dof);
gp = [-1/sqrt(3), 1/sqrt(3)]; % Gaussian points (2x2 integration)
fprintf('Assembling Stiffness Matrix...\n');
for e = 1:n_elem
% Get node IDs and coordinates
sctr = element(e, :);
coords = node(sctr, :);
% DOF mapping for this element
sctrB = zeros(1, 20); % 4 nodes * 5 DOFs
for i = 1:4
    sctrB((i-1)*5+

Introduction

Composite plates are widely used in various engineering applications, such as aerospace, automotive, and civil engineering, due to their high strength-to-weight ratio, stiffness, and resistance to corrosion. However, analyzing the bending behavior of composite plates can be challenging due to their complex material properties and laminate configurations.

Classical Laminate Theory (CLT)

The Classical Laminate Theory (CLT) is a widely used method for analyzing the bending behavior of composite plates. The CLT assumes that the plate is thin, and the deformations are small. The theory is based on the following assumptions:

The CLT provides a set of equations that relate the mid-plane strains and curvatures to the applied loads. The equations are:

$$\beginbmatrix \epsilon_x \ \epsilon_y \ \gamma_xy \endbmatrix = \beginbmatrix \epsilon_x^0 \ \epsilon_y^0 \ \gamma_xy^0 \endbmatrix + z \beginbmatrix \kappa_x \ \kappa_y \ \kappa_xy \endbmatrix$$

$$\beginbmatrix \sigma_x \ \sigma_y \ \tau_xy \endbmatrix = \beginbmatrix Q_11 & Q_12 & Q_16 \ Q_12 & Q_22 & Q_26 \ Q_16 & Q_26 & Q_66 \endbmatrix \beginbmatrix \epsilon_x \ \epsilon_y \ \gamma_xy \endbmatrix$$

$$\beginbmatrix N_x \ N_y \ N_xy \endbmatrix = \int_-h/2^h/2 \beginbmatrix \sigma_x \ \sigma_y \ \tau_xy \endbmatrix dz$$ Composite Plate Bending Analysis With Matlab Code

$$\beginbmatrix M_x \ M_y \ M_xy \endbmatrix = \int_-h/2^h/2 \beginbmatrix \sigma_x \ \sigma_y \ \tau_xy \endbmatrix z dz$$

MATLAB Code

Here is a sample MATLAB code to perform composite plate bending analysis using the CLT:

function [displacements, stresses, strains] = composite_plate_bending_analysis(E1, E2, nu12, G12, t, Lx, Ly, q)
% Material properties
Q11 = E1 / (1 - nu12^2);
Q22 = E2 / (1 - nu12^2);
Q12 = nu12 * E2 / (1 - nu12^2);
Q16 = 0;
Q26 = 0;
Q66 = G12;
% Laminate properties
h = sum(t);
z = [-h/2; h/2];
% Applied load
q = 1;  % uniform load
% Mid-plane strains and curvatures
ex0 = 0;
ey0 = 0;
gxy0 = 0;
kx = -q / (24 * D);
ky = -q / (24 * D);
kxy = 0;
% Bending stiffness matrix
D = zeros(3,3);
for i = 1:length(t)
    zi = z(i);
    D = D + (1/3) * (zi^3 - zi^2 * t(i)) * [Q11, Q12, Q16; Q12, Q22, Q26; Q16, Q26, Q66];
end
% Displacements
w = (q / (24 * D)) * (Lx^4 + Ly^4);
% Stresses and strains
stresses = zeros(3,1);
strains = zeros(3,1);
for i = 1:length(t)
    zi = z(i);
    stresses = [Q11, Q12, Q16; Q12, Q22, Q26; Q16, Q26, Q66] * (ex0 + zi * kx);
    strains = [ex0 + zi * kx; ey0 + zi * ky; gxy0 + zi * kxy];
end
end

Example Usage

To use the code, simply call the function with the required input arguments:

E1 = 100e9;  % Young's modulus in the 1-direction (Pa)
E2 = 50e9;   % Young's modulus in the 2-direction (Pa)
nu12 = 0.3;  % Poisson's ratio
G12 = 20e9;  % Shear modulus (Pa)
t = [0.001; 0.002];  % layer thicknesses (m)
Lx = 1;      % plate length in the x-direction (m)
Ly = 1;      % plate length in the y-direction (m)
q = 1000;    % uniform load (Pa)
[displacements, stresses, strains] = composite_plate_bending_analysis(E1, E2, nu12, G12, t, Lx, Ly, q);

This code provides a basic framework for analyzing the bending behavior of composite plates using the Classical Laminate Theory. However, please note that this is a simplified example and real-world applications may require more complex analyses, such as considering non-uniform loads, boundary conditions, and material nonlinearity.

Introduction

Composite plates are widely used in various engineering applications, such as aerospace, automotive, and civil engineering, due to their high strength-to-weight ratio and stiffness. However, analyzing the bending behavior of composite plates can be complex due to their anisotropic material properties. This guide provides an overview of composite plate bending analysis using MATLAB code.

Theoretical Background

The bending behavior of composite plates can be analyzed using the following theories:

Governing Equations

The governing equations for composite plate bending analysis using FSDT are:

$$\beginbmatrix \frac\partial^2 M_x\partial x^2 + 2\frac\partial^2 M_xy\partial x \partial y + \frac\partial^2 M_y\partial y^2 = q \ \frac\partial^2 M_x\partial x \partial y + \frac\partial^2 M_xy\partial x^2 + \frac\partial^2 M_y\partial y^2 = 0 \endbmatrix$$

$$\beginbmatrix M_x \ M_y \ M_xy \endbmatrix = \beginbmatrix D_11 & D_12 & D_16 \ D_12 & D_22 & D_26 \ D_16 & D_26 & D_66 \endbmatrix \beginbmatrix \kappa_x \ \kappa_y \ \kappa_xy \endbmatrix$$

where $M_x$, $M_y$, and $M_xy$ are the bending and twisting moments, $q$ is the transverse load, $D_ij$ are the flexural stiffnesses, and $\kappa_x$, $\kappa_y$, and $\kappa_xy$ are the curvatures.

MATLAB Code

The following MATLAB code performs a bending analysis of a composite plate using FSDT:

% Define plate properties
a = 10;  % plate length (m)
b = 10;  % plate width (m)
h = 0.1;  % plate thickness (m)
E1 = 100e9;  % Young's modulus in x-direction (Pa)
E2 = 50e9;  % Young's modulus in y-direction (Pa)
G12 = 20e9;  % shear modulus (Pa)
nu12 = 0.3;  % Poisson's ratio
q = 1000;  % transverse load (Pa)
% Define material stiffness matrix
Q11 = E1 / (1 - nu12^2);
Q22 = E2 / (1 - nu12^2);
Q12 = nu12 * Q11;
Q66 = G12;
Q16 = 0;
Q26 = 0;
% Define flexural stiffness matrix
D11 = (1/3) * (Q11 * h^3);
D22 = (1/3) * (Q22 * h^3);
D12 = (1/3) * (Q12 * h^3);
D66 = (1/3) * (Q66 * h^3);
D16 = (1/3) * (Q16 * h^3);
D26 = (1/3) * (Q26 * h^3);
% Assemble global stiffness matrix
K = [D11, D12, D16; D12, D22, D26; D16, D26, D66];
% Solve for deflection and rotation
w = q / (D11 * (1 - nu12^2));
theta_x = - (D12 / D11) * w;
theta_y = - (D26 / D22) * w;
% Display results
fprintf('Deflection: %.2f mm\n', w * 1000);
fprintf('Rotation (x): %.2f degrees\n', theta_x * 180 / pi);
fprintf('Rotation (y): %.2f degrees\n', theta_y * 180 / pi);

This code defines the plate properties, material stiffness matrix, and flexural stiffness matrix. It then assembles the global stiffness matrix and solves for the deflection and rotation of the plate under a transverse load.

Example Use Case

Consider a square composite plate with a length and width of 10 m and a thickness of 0.1 m. The plate is made of a carbon/epoxy material with the following properties:

The plate is subjected to a transverse load of 1000 Pa. Using the MATLAB code above, the deflection and rotation of the plate can be calculated. If the fibers are oriented at an angle

Limitations and Future Work

This guide provides a basic introduction to composite plate bending analysis using MATLAB code. However, there are several limitations and areas for future work:

Analyzing composite plate bending in MATLAB typically involves implementing Classical Laminated Plate Theory (CLPT) or First-order Shear Deformation Theory (FSDT) to calculate structural responses like deflection and stress distributions. Key Analytical Concepts

ABD Matrix: The core of composite analysis, where A represents extensional stiffness, B represents coupling stiffness (essential for unsymmetric layups), and D represents bending stiffness. Theories used: CLPT: Best for thin plates ( ) where shear deformation is negligible.

FSDT (Mindlin-Reissner): Accounts for shear deformation, making it necessary for thicker plates.

Quasi-3D Theory: Provides higher accuracy for transverse displacement by accounting for variation through the laminate thickness. Implementation in MATLAB A typical script for bending analysis follows these steps:

Composite plate bending analysis evaluates how laminated structures—made of layers with varying fiber orientations—deform under transverse loads. Unlike isotropic materials, these plates exhibit directional mechanical properties (anisotropy), requiring specialized theories like Classical Laminated Plate Theory (CLPT) for thin plates or First-order Shear Deformation Theory (FSDT) for thicker ones. 1. Calculate Laminate Stiffness (ABD Matrix)

The first step is determining the ABD matrix, which relates mid-plane strains and curvatures to applied resultants (forces and moments). Determine Reduced Stiffness ( ): Calculate the matrix for each layer using its material properties ( Transform Stiffness (

[Q]̄modified open bracket cap Q close bracket with bar above ): Use a transformation matrix based on each ply's orientation angle ( ) to convert local stiffness to global coordinates. Assemble Matrices: Integrate the

[Q]̄modified open bracket cap Q close bracket with bar above matrices through the thickness to find:

(Extensional Stiffness): Relates in-plane forces to strains.

(Coupling Stiffness): Relates in-plane forces to curvatures (zero for symmetric laminates). (Bending Stiffness): Relates moments to curvatures. 2. Formulate Governing Equations

Using the determined stiffness, the relationship between resultants and deformations is expressed as:

[NM]=[ABBD][ϵ0κ]the 2 by 1 column matrix; cap N, cap M end-matrix; equals the 2 by 2 matrix; Row 1: cap A, cap B; Row 2: cap B, cap D end-matrix; the 2 by 1 column matrix; epsilon to the 0 power, kappa end-matrix; For a simple bending analysis (where ), you solve for mid-plane curvatures ( ) to find the plate's deflection. 3. MATLAB Code Implementation

The following code snippet demonstrates the calculation of the ABD matrix for a laminated composite, a critical precursor to bending analysis.

% Material Properties (Graphite/Epoxy Example) E1 = 181e9; E2 = 10.3e9; nu12 = 0.28; G12 = 7.17e9; nu21 = nu12 * E2 / E1; % Layup Definition [theta1, theta2, ...] and thickness angles = [0, 45, -45, 90]; % degrees tk = 0.005; % thickness of each ply (m) n = length(angles); h = n * tk; z = -h/2 : tk : h/2; % coordinate of each interface % Initialize ABD matrices A = zeros(3,3); B = zeros(3,3); D = zeros(3,3); % Reduced stiffness matrix [Q] in principal material directions Q = [E1/(1-nu12*nu21), nu12*E2/(1-nu12*nu21), 0; nu12*E2/(1-nu12*nu21), E2/(1-nu12*nu21), 0; 0, 0, G12]; for k = 1:n theta = deg2rad(angles(k)); m = cos(theta); n_s = sin(theta); % Transformation Matrix [T] T = [m^2, n_s^2, 2*m*n_s; n_s^2, m^2, -2*m*n_s; -m*n_s, m*n_s, m^2-n_s^2]; % Reuter's Matrix [R] R = [1 0 0; 0 1 0; 0 0 2]; % Transformed Stiffness Matrix [Q_bar] Q_bar = inv(T) * Q * R * T * inv(R); % Accumulate into ABD A = A + Q_bar * (z(k+1) - z(k)); B = B + 0.5 * Q_bar * (z(k+1)^2 - z(k)^2); D = D + (1/3) * Q_bar * (z(k+1)^3 - z(k)^3); end disp('Bending Stiffness Matrix [D]:'); disp(D); Use code with caution. Copied to clipboard 4. Displacement and Stress Analysis

Once the stiffness is known, you can solve for displacement (

) using numerical methods like the Finite Element Method (FEM).

Mesh Generation: Divide the plate into elements (e.g., 4-node Q4 elements).

Global Assembly: Assemble the global stiffness matrix from element matrices derived via FSDT or CLPT.

Apply Boundary Conditions: Define supports (e.g., simply supported or clamped). Introduction Composite plates are widely used in various

Solve & Post-process: Calculate deflections and then retrieve global/local stresses for each layer to check for failure (using criteria like Tsai-Hill).

For more complex geometries or non-linear effects, practitioners often transition from custom MATLAB scripts to specialized software like ABAQUS or ANSYS. Composite Plate Bending Analysis With Matlab Code

Bending analysis of composite plates typically uses Classical Laminate Plate Theory (CLPT) for thin plates or First-order Shear Deformation Theory (FSDT)

for thicker structures. The process involves calculating the laminate stiffness (the ABD matrix), solving for mid-plane deformations, and then determining layer-by-layer stresses. ScienceDirect.com 1. Define Material and Layer Properties

First, define the properties of each lamina (layer), including Young's moduli ( ), shear modulus ( cap G sub 12 ), and Poisson's ratio ( ). For each layer , specify the thickness and the fiber orientation angle theta sub k 2. Calculate the Reduced Stiffness Matrix ( The reduced stiffness matrix for an orthotropic lamina in its principal directions is:

open bracket cap Q close bracket equals the 3 by 3 matrix; Row 1: cap Q sub 11, cap Q sub 12, 0; Row 2: cap Q sub 12, cap Q sub 22, 0; Row 3: 0, 0, cap Q sub 66 end-matrix; SCIRP Open Access 3. Transform Stiffness to Global Coordinates ( Each layer's stiffness must be transformed into the global

coordinate system of the plate using the transformation matrix based on the fiber angle SCIRP Open Access 4. Compute the ABD Matrix The ABD matrix relates the applied forces and moments to the mid-plane strains and curvatures SCIRP Open Access

the 2 by 1 column matrix; Row 1: the set cap N end-set, Row 2: the set cap M end-set end-matrix; equals the 2 by 2 matrix; Row 1: Column 1: open bracket cap A close bracket, Column 2: open bracket cap B close bracket; Row 2: Column 1: open bracket cap B close bracket, Column 2: open bracket cap D close bracket end-matrix; the 2 by 1 column matrix; Row 1: the set epsilon to the 0 power end-set, Row 2: the set kappa end-set end-matrix; (Extensional Stiffness): (Coupling Stiffness): (Bending Stiffness): is the vertical position of the k raised to the t h power layer relative to the mid-plane. SCIRP Open Access 5. Solve for Deformations and Stresses

Given the applied moments (e.g., from a distributed load), solve the simultaneous equations to find

. Use these to find the global and local stresses in each layer to perform failure analysis. Harvard University Matlab Code Implementation This script calculates the ABD matrix for a symmetric % Composite Plate Bending Analysis - ABD Matrix Calculation % Material properties (Pa) v21 = v12 * E2 / E1; angles = [ % Layer orientations (degrees) % Thickness of each layer (m) n = length(angles); h_total = n * t; z = -h_total/ : t : h_total/ % Z-coordinates of layer interfaces % Reduced stiffness matrix [Q] Q11 = E1/( -v12*v21); Q22 = E2/( -v12*v21); Q12 = v12*E2/( -v12*v21); Q66 = G12; Q = [Q11 Q12 Q66];

A = zeros( ); B = zeros( ); D = zeros(

:n theta = deg2rad(angles(k)); c = cos(theta); s = sin(theta); % Transformation matrix [T] *s*c; -s*c s*c c^ ]; R = [ % Reuter's matrix Qbar = inv(T) * Q * R * T * inv(R); % Accumulate A, B, D matrices A = A + Qbar * (z(k+ ) - z(k)); B = B + * Qbar * (z(k+ ); D = D + ( ) * Qbar * (z(k+ 'Bending Stiffness Matrix [D]:' ); disp(D); Use code with caution. Copied to clipboard

For advanced FEM-based plate bending (including meshing and displacement plots), you can use the Partial Differential Equation Toolbox Layup Analysis Tool MATLAB Central deflection

for a specific boundary condition like a simply supported plate? Structural Analysis Using Finite Element Method in MATLAB

Learn how to perform structural analysis using the finite element method with Partial Differential Equation ToolboxTM in MATLAB®.

Flexural analysis of laminated composite plates - ScienceDirect

| Issue | What to Check in the Code | | :--- | :--- | | Shear Locking for Thin Plates | If using FSDT (Mindlin) with linear shape functions, the code may be overly stiff for thin plates. Look for selective reduced integration (using int points differently for shear vs. bending). | | Classical Theory (CLPT) Overly Stiff | Code using CLPT ignores transverse shear deformation. It will be accurate for very thin plates (span/thickness > 50) but will under-predict deflection for moderately thick composite plates. | | Boundary Conditions | Many student codes only handle simply supported (SS) or fully clamped (CC). Be wary if you need free edges or symmetry conditions. | | Stress Recovery | The best codes output stress per layer (top, middle, bottom). Weak codes only output global moments. Ensure the code you review includes Q_bar back-transformation to get stresses in material coordinates. | | Convergence | A good code will have a convergence study. A bad one assumes one mesh works for all. |

The laminate stiffness is represented by three matrices:

For a laminate with N layers, the bending stiffness matrix D (3×3) is defined as:

D_ij = (1/3) * Σ_k=1^N (Q_ij)_k * (z_k^3 - z_k-1^3)

Where ( Q_ij ) are transformed reduced stiffnesses of the k-th layer at angle θ.

The moment-curvature relation:

 Mxx ; Myy ; Mxy  = [D] *  κxx ; κyy ; κxy 

The code checks if the laminate is symmetric ($[B]$ is effectively zero).