[ W_mn = \fracQ_mn \pi^4 \left[ D_11\left(\fracma\right)^4 + 2(D_12+2D_66)\left(\fracma\right)^2\left(\fracnb\right)^2 + D_22\left(\fracnb\right)^4 + 4 D_16\left(\fracma\right)^3\left(\fracnb\right) + 4 D_26\left(\fracma\right)\left(\fracnb\right)^3 \right] ]
end
function [N, dN_dxi, detJ] = shape_functions(xy, xi, eta) % xy: nodal coordinates, 4x2 N = zeros(4,1); dN_dxi = zeros(2,4); % dN/dxi and dN/deta for i = 1:4 xi_i = [-1, 1, 1, -1]; eta_i = [-1, -1, 1, 1]; N(i) = 0.25 * (1 + xi_i(i)*xi) * (1 + eta_i(i)*eta); dN_dxi(1,i) = 0.25 * xi_i(i) * (1 + eta_i(i)*eta); dN_dxi(2,i) = 0.25 * eta_i(i) * (1 + xi_i(i)*xi); end % Jacobian J = zeros(2,2); for i = 1:4 J(1,1) = J(1,1) + dN_dxi(1,i) * xy(i,1); J(1,2) = J(1,2) + dN_dxi(1,i) * xy(i,2); J(2,1) = J(2,1) + dN_dxi(2,i) * xy(i,1); J(2,2) = J(2,2) + dN_dxi(2,i) * xy(i,2); end detJ = det(J); % Derivatives w.r.t. physical coordinates dN_dxi = J \ dN_dxi; end
A = A + Q_bar * (z_top - z_bot); B = B + 0.5 * Q_bar * (z_top^2 - z_bot^2); D = D + (1/3) * Q_bar * (z_top^3 - z_bot^3); As = As + Q_s_bar * (z_top - z_bot);
w(x,y)=∑m=1∞∑n=1∞Wmnsin(mπxa)sin(nπyb)w open paren x comma y close paren equals sum from m equals 1 to infinity of sum from n equals 1 to infinity of cap W sub m n end-sub sine open paren the fraction with numerator m pi x and denominator a end-fraction close paren sine open paren the fraction with numerator n pi y and denominator b end-fraction close paren
The double-loop summation converges rapidly to find the maximum deflection point located precisely at the center of the plate ( Visualization Output Composite Plate Bending Analysis With Matlab Code
For a uniformly distributed transverse load (q) (force per area), the equivalent nodal forces for (w) are:
% ========================================================================= % COMPOSITE PLATE BENDING ANALYSIS USING NAVIER'S SOLUTION (CLPT) % ========================================================================= clear; clc; close all; %% 1. Geometry and Material Properties a = 0.5; % Length of plate (m) b = 0.5; % Width of plate (m) q0 = -10000; % Uniform distributed load (Pa, negative downwards) % Material: Graphite/Epoxy E1 = 138e9; % Longitudinal Elastic Modulus (Pa) E2 = 8.9e9; % Transverse Elastic Modulus (Pa) G12 = 7.1e9; % Shear Modulus (Pa) nu12 = 0.3; % Major Poisson's ratio nu21 = (E2/E1)*nu12; % Minor Poisson's ratio %% 2. Laminate Stacking Sequence % Symmetric cross-ply laminate [0 / 90 / 90 / 0] theta = [0, 90, 90, 0]; ply_thickness = 0.00125 * ones(1, length(theta)); % 1.25mm per ply total_thickness = sum(ply_thickness); % Calculate ply interface z-coordinates relative to the mid-plane n_plies = length(theta); z = zeros(1, n_plies + 1); z(1) = -total_thickness / 2; for k = 1:n_plies z(k+1) = z(k) + ply_thickness(k); end %% 3. Calculate ABD Matrices A = zeros(3,3); B = zeros(3,3); D = zeros(3,3); % Reduced Stiffness Matrix [Q] in material coordinates 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]; for k = 1:n_plies rad = deg2rad(theta(k)); m = cos(rad); n = sin(rad); % Transformation Matrix [T] T = [m^2, n^2, 2*m*n; n^2, m^2, -2*m*n; -m*n, m*n, m^2-n^2]; % Reuter's Matrix [R] R = [1 0 0; 0 1 0; 0 0 2]; % Transformed Reduced Stiffness Matrix [Q_bar] Q_bar = inv(T) * Q * R * T * inv(R); % Integrate across thickness to populate A, B, D matrices A = A + Q_bar * (z(k+1) - z(k)); B = B + Q_bar * (z(k+1)^2 - z(k)^2) / 2; D = D + Q_bar * (z(k+1)^3 - z(k)^3) / 3; end disp('--- Calculated Bending Stiffness Matrix [D] (N-m) ---'); disp(D); %% 4. Navier's Solution Implementation Nx = 50; Ny = 50; % Grid resolution for plotting x = linspace(0, a, Nx); y = linspace(0, b, Ny); [X, Y] = meshgrid(x, y); W = zeros(Nx, Ny); % Convergence limits for Fourier expansion m_max = 29; n_max = 29; for m = 1:2:m_max for n = 1:2:n_max % Fourier load amplitude coefficient Qmn = (16 * q0) / (pi^2 * m * n); % Denominator formula for symmetric cross-ply bending denom = pi^4 * (D(1,1)*(m/a)^4 + 2*(D(1,2) + 2*D(3,3))*(m/a)^2*(n/b)^2 + D(2,2)*(n/b)^4); Wmn = Qmn / denom; % Accumulate spatial components W = W + Wmn * sin(m*pi*X/a) .* sin(n*pi*Y/b); end end %% 5. Results Presentation and Plotting max_deflection = min(W(:)); % Deflection is down (negative) fprintf('Maximum Center Deflection: %.4f mm\n', abs(max_deflection) * 1000); figure('Color', [1 1 1]); surf(X, Y, W * 1000, 'EdgeColor', 'none'); colormap('jet'); colorbar; title('Displacement Field of Simply Supported Composite Plate'); xlabel('X-Axis / Length (m)'); ylabel('Y-Axis / Width (m)'); zlabel('Deflection w (mm)'); view(-37.5, 30); grid on; Use code with caution. 4. Engineering Analysis & Interpretation
% Loop over odd m and n (others zero for uniform load) for m = 1:2:Mmax for n = 1:2:Nmax % Determine Qmn based on load type if strcmp(load_type, 'sinusoidal') if m == 1 && n == 1 Qmn(m,n) = q0; else Qmn(m,n) = 0; end elseif strcmp(load_type, 'uniform') Qmn(m,n) = 16 * q0 / (pi^2 * m * n); else error('Load type not recognized. Use ''sinusoidal'' or ''uniform''.'); end
end
% Extract D matrix components for specially orthotropic plate D11 = D(1,1); D12 = D(1,2); D22 = D(2,2); D66 = D(3,3); D16 = D(1,3); D26 = D(2,3); [ W_mn = \fracQ_mn \pi^4 \left[ D_11\left(\fracma\right)^4 +
matrix calculates directly to 0. This guarantees that pure bending moments will not generate unwanted twisting or stretching deformations.
Matrix (Coupling Stiffness): Relates in-plane forces to curvatures (and moments to strains). For symmetric laminates, Matrix (Bending Stiffness): Relates moments to curvatures. 1.2. Governing Equations for Bending
% Transformation for shear: Q_s_bar = T_s * Q_s * T_s' T_s = [c, s; -s, c]; Q_s_bar = T_s * Q_s * T_s';
Because the boundaries are simply supported and the pressure load is completely uniform, the peak transverse deflection happens exactly at the geometric center of the plate:
[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; Laminate Stacking Sequence % Symmetric cross-ply laminate [0
(Coupling Stiffness Matrix): Relates in-plane forces to bending curvatures, and bending moments to in-plane strains. For symmetric laminates,
% Expand to 20x20 (u,v,w,θx,θy per node) % Here we assemble directly into 5 DOF format % For simplicity, we use block matrices % Actual implementation would map correctly % We'll assemble Ke as 5x5 blocks per node
Straight lines normal to the mid-surface remain normal to the mid-surface after deformation. The thickness of the plate does not change. 1.1. Laminate Constitutive Relations The relationship between forces , mid-plane strains , and curvatures is given by the ABD matrix:
% Pre-allocate A = zeros(3,3); B = zeros(3,3); D = zeros(3,3);