Composite materials, particularly fiber-reinforced laminates, are extensively used in aerospace, automotive, and civil engineering due to their high strength-to-weight and stiffness-to-weight ratios. However, their anisotropic and layered nature makes bending analysis more complex than isotropic plates.
%% 4. Element Stiffness Matrix (4-node quadrilateral, 5 DOF/node) K_global = sparse(nDofs, nDofs); F_global = zeros(nDofs, 1);
For an orthotropic lamina (e.g., fiber-reinforced composite) oriented at an angle θ, the stress-strain relationship in the lamina's principal coordinates (1,2) is:
D11𝜕4w𝜕x4+2(D12+2D66)𝜕4w𝜕x2𝜕y2+D22𝜕4w𝜕y4=q(x,y)cap D sub 11 partial to the fourth power w over partial x to the fourth power end-fraction plus 2 open paren cap D sub 12 plus 2 cap D sub 66 close paren the fraction with numerator partial to the fourth power w and denominator partial x squared partial y squared end-fraction plus cap D sub 22 partial to the fourth power w over partial y to the fourth power end-fraction equals q open paren x comma y close paren Navier's Solution Method Composite Plate Bending Analysis With Matlab Code
%% Navier analytical solution (simply supported, symmetric laminate -> B=0) % Number of terms m_max = 51; n_max = 51; w_max = 0; x_center = a/2; y_center = b/2; for m = 1:2:m_max for n = 1:2:n_max % Fourier coefficient for uniform load qmn = (16 q0)/(m n pi^2); % Denominator term1 = D(1,1) (m pi/a)^4; term2 = 2 (D(1,2)+2 D(3,3)) (m pi/a)^2 (n pi/b)^2; term3 = D(2,2) (n pi/b)^4; term4 = 4 D(1,3) (m pi/a)^3*(n pi/b); term5 = 4 D(2,3) (m pi/a) (n pi/b)^3; denom = term1 + term2 + term3 + term4 + term5; Wmn = qmn / denom; w_contrib = Wmn * sin(m pi x_center/a) * sin(n pi y_center/b); w_max = w_max + w_contrib; end end fprintf('Navier solution: Maximum deflection = %.6f mm\n', w_max*1000);
fprintf('\nMaximum deflection: %.4f mm\n', max(W)*1000); fprintf('Minimum deflection: %.4f mm\n', min(W)*1000);
Define supports (e.g., simply supported or clamped). For symmetric laminates
Running the code for a 0.5m×0.5m, 5mm thick cross-ply [0/90]s plate under 1 kPa pressure yields:
clear; clc; close all;
% ========================================================================= % Composite Plate Bending Analysis using CLPT and Navier's Method % ========================================================================= clear; clc; close all; %% 1. Input Parameters % Plate Geometry a = 0.5; % Length along x-axis (m) b = 0.5; % Width along y-axis (m) q0 = 10e3; % Uniformly distributed load (Pa) % Material Properties (Carbon/Epoxy) E1 = 150e9; % Longitudinal Modulus (Pa) E2 = 10e9; % Transverse Modulus (Pa) G12 = 5e9; % Shear Modulus (Pa) nu12 = 0.3; % Poisson's ratio nu21 = (E2/E1)*nu12; % Layup Definition % Symmetric Cross-ply [0 / 90 / 90 / 0] theta =; ply_thickness = 0.00125 * ones(1, length(theta)); % Thickness of each ply (m) h = sum(ply_thickness); %% 2. Calculate Laminate Stiffness (ABD Matrices) % Reduced stiffness matrix in material coordinates (Q) Q11 = E1 / (1 - nu12*nu21); Q12 = (nu12*E2) / (1 - nu12*nu21); Q22 = E2 / (1 - nu12*nu21); Q66 = G12; Q = [Q11, Q12, 0; Q12, Q22, 0; 0, 0, Q66]; % Initialize ABD Matrices A = zeros(3,3); B = zeros(3,3); D = zeros(3,3); % Z-coordinates for each ply boundary z = zeros(1, length(theta) + 1); z(1) = -h/2; for k = 1:length(theta) z(k+1) = z(k) + ply_thickness(k); end % Assemble matrices for k = 1:length(theta) % Angle transformation t = deg2rad(theta(k)); m = cos(t); n = sin(t); % Transformed reduced stiffness matrix (Qbar) T = [m^2, n^2, 2*m*n; n^2, m^2, -2*m*n; -m*n, m*n, m^2-n^2]; Qbar = T \ Q * (T.^-1).'; % Simplified transformation representation % Alternative explicit Qbar calculation for accuracy Qbar(1,1) = Q11*m^4 + 2*(Q12 + 2*Q66)*m^2*n^2 + Q22*n^4; Qbar(1,2) = (Q11 + Q22 - 4*Q66)*m^2*n^2 + Q12*(m^4 + n^4); Qbar(2,2) = Q11*n^4 + 2*(Q12 + 2*Q66)*m^2*n^2 + Q22*m^4; Qbar(1,3) = (Q11 - Q12 - 2*Q66)*m^3*n + (Q12 - Q22 + 2*Q66)*m*n^3; Qbar(2,3) = (Q11 - Q12 - 2*Q66)*m*n^3 + (Q12 - Q22 + 2*Q66)*m^3*n; Qbar(3,3) = (Q11 + Q22 - 2*Q12 - 2*Q66)*m^2*n^2 + Q66*(m^4 + n^4); Qbar(2,1) = Qbar(1,2); Qbar(3,1) = Qbar(1,3); Qbar(3,2) = Qbar(2,3); % Integrate across thickness A = A + Qbar * (z(k+1) - z(k)); B = B + 0.5 * Qbar * (z(k+1)^2 - z(k)^2); D = D + (1/3) * Qbar * (z(k+1)^3 - z(k)^3); end %% 3. Navier Solution for Deflection Nx = 50; Ny = 50; x = linspace(0, a, Nx); y = linspace(0, b, Ny); [X, Y] = meshgrid(x, y); w = zeros(size(X)); M_terms = 49; % Number of Fourier terms (odd numbers) N_terms = 49; for m = 1:2:M_terms for n = 1:2:N_terms % Load coefficient Qmn = (16 * q0) / (pi^2 * m * n); % Denominator factor from structural stiffness alpha = m * pi / a; beta = n * pi / b; denom = D(1,1)*alpha^4 + 2*(D(1,2) + 2*D(3,3))*alpha^2*beta^2 + D(2,2)*beta^4; % Deflection coefficient Wmn = Qmn / denom; % Accumulate deflection field w = w + Wmn * sin(alpha * X) .* sin(beta * Y); end end %% 4. Plot Results figure; surf(X, Y, w * 1e3, 'EdgeColor', 'none'); colormap jet; colorbar; title('Transverse Deflection of Composite Plate'); xlabel('X-axis (m)'); ylabel('Y-axis (m)'); zlabel('Deflection (mm)'); view(45, 30); Use code with caution. Result Interpretation Nwx = zeros(1
% Shape functions for w (Hermitian-type, non-conforming) % We use standard Kirchhoff plate element (Zienkiewicz's non-conforming) % Define basis functions: Nw = zeros(1,4); Nwx = zeros(1,4); % dNw/dx Nwy = zeros(1,4); % dNw/dy
We use bilinear shape functions for w and rotations derived from the Kirchhoff constraint. A practical alternative is the discrete Kirchhoff quadrilateral (DKQ) element, but for simplicity we adopt the conforming rectangular element with 12 DOFs.
Matrix (Coupling Stiffness): Links in-plane forces to bending curvatures. For symmetric laminates,