Composite Plate Bending Analysis With Matlab Code [hot]

Unlike black-box commercial FEA software (ANSYS, Abaqus), the MATLAB code lets you see every matrix: the ABD matrix, the element stiffness matrix, the shear correction factor, and the assembly process. You truly learn why a composite plate bends differently from an isotropic one.

Dij=13∑k=1N(Q̄ij)k(zk3−zk−13)cap D sub i j end-sub equals one-third sum from k equals 1 to cap N of open paren cap Q bar sub i j end-sub close paren sub k open paren z sub k cubed minus z sub k minus 1 end-sub cubed close paren Transverse Shear Stiffness Matrix (A_s) For FSDT formulations, an additional matrix Ascap A sub s accounts for transverse shear stiffness:

This article explores the theoretical foundation of composite plate bending and provides a complete Matlab script to calculate deflections and stresses. Theoretical Framework Classical Laminated Plate Theory (CLPT) Composite Plate Bending Analysis With Matlab Code

𝜕Nx𝜕x+𝜕Nxy𝜕y=0the fraction with numerator partial cap N sub x and denominator partial x end-fraction plus the fraction with numerator partial cap N sub x y end-sub and denominator partial y end-fraction equals 0

[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; % 2x2 points and weights 1/3

% Material properties (T300/5208 carbon/epoxy) E1 = 181e9; % longitudinal modulus [Pa] E2 = 10.3e9; % transverse modulus [Pa] nu12 = 0.28; % major Poisson's ratio G12 = 7.17e9; % in-plane shear modulus [Pa]

%% COMPOSITE PLATE BENDING ANALYSIS USING FSDT & FEM clc; clear; close all; %% 1. GEOMETRY & MESH DEFINITION L_x = 1.0; % Plate length along X-axis (m) L_y = 1.0; % Plate width along Y-axis (m) Nx_elem = 10; % Number of elements along X Ny_elem = 10; % Number of elements along Y % Material Properties (e.g., Graphite-Epoxy) E1 = 175e9; % Longitudinal Elastic Modulus (Pa) E2 = 7e9; % Transverse Elastic Modulus (Pa) G12 = 3.5e9; % In-plane Shear Modulus (Pa) G13 = 3.5e9; % Transverse Shear Modulus (Pa) G23 = 1.4e9; % Transverse Shear Modulus (Pa) nu12 = 0.25; % Major Poisson's ratio nu21 = (E2/E1)*nu12; % Stacking Sequence and Thicknesses % [0/90/90/0] Symmetric Cross-ply Laminate theta = [0, 90, 90, 0]; % Ply angles in degrees ply_t = [0.005, 0.005, 0.005, 0.005]; % Individual ply thicknesses (m) total_h = sum(ply_t); % Total laminate thickness num_plies = length(theta); % Loading q0 = -10000; % Uniformly distributed transverse load (N/m^2) %% 2. CALCULATE LAMINATE ABD & As MATRICES [A, B, D, As] = compute_ABD_matrices(E1, E2, G12, G13, G23, nu12, nu21, theta, ply_t); %% 3. FINITE ELEMENT MESH GENERATION [node_coords, elem_nodes] = generate_mesh(L_x, L_y, Nx_elem, Ny_elem); num_nodes = size(node_coords, 1); num_elements = size(elem_nodes, 1); dofs_per_node = 5; % [u, v, w, rot_x, rot_y] total_dofs = num_nodes * dofs_per_node; %% 4. GLOBAL STIFFNESS MATRIX AND LOAD VECTOR ASSEMBLY K_global = zeros(total_dofs, total_dofs); F_global = zeros(total_dofs, 1); % Gauss Quadrature points and weights for integration % Bending (Full Integration: 2x2) gauss_pt = [-1/sqrt(3), 1/sqrt(3)]; gauss_wt = [1.0, 1.0]; % Shear (Reduced Integration: 1x1 to avoid shear locking) gauss_pt_r = 0.0; gauss_wt_r = 2.0; for e = 1:num_elements nodes = elem_nodes(e, :); X = node_coords(nodes, 1); Y = node_coords(nodes, 2); K_elem = zeros(20, 20); F_elem = zeros(20, 1); % Element sizing (for rectangular elements) a_elem = (X(2) - X(1)) / 2; b_elem = (Y(4) - Y(1)) / 2; Jacobian_det = a_elem * b_elem; % --- Bending and Extension Component (2x2 Integration) --- for i = 1:2 for j = 1:2 xi = gauss_pt(i); eta = gauss_pt(j); [~, B_b] = get_shape_functions(xi, eta, a_elem, b_elem); % Constitutive stiffness mapping ABD_matrix = [A, B; B, D]; K_elem = K_elem + (B_b' * ABD_matrix * B_b) * Jacobian_det * gauss_wt(i) * gauss_wt(j); end end % --- Transverse Shear Component (1x1 Reduced Integration) --- for i = 1:1 for j = 1:1 xi = gauss_pt_r(i); eta = gauss_pt_r(j); [N, ~] = get_shape_functions(xi, eta, a_elem, b_elem); [~, B_s] = get_shear_B_matrix(xi, eta, a_elem, b_elem); K_elem = K_elem + (B_s' * As * B_s) * Jacobian_det * gauss_wt_r(i) * gauss_wt_r(j); % Distribute uniform force to element transverse DOFs (w parameter) for n = 1:4 F_elem((n-1)*5 + 3) = F_elem((n-1)*5 + 3) + N(n) * q0 * Jacobian_det * gauss_wt_r(i) * gauss_wt_r(j); end end end % Assembly into Global System global_dofs = []; for n = 1:4 global_dofs = [global_dofs, (nodes(n)-1)*5 + (1:5)]; end K_global(global_dofs, global_dofs) = K_global(global_dofs, global_dofs) + K_elem; F_global(global_dofs) = F_global(global_dofs) + F_elem; end %% 5. BOUNDARY CONDITIONS (Simply Supported On All Edges) % SSSS-1 condition: w = 0, rot_x = 0 or rot_y = 0 along edges bc_nodes = []; for i = 1:num_nodes x = node_coords(i, 1); y = node_coords(i, 2); if x == 0 || x == L_x || y == 0 || y == L_y % Constrain displacement w (DOF 3) bc_nodes = [bc_nodes, (i-1)*5 + 3]; if x == 0 || x == L_x bc_nodes = [bc_nodes, (i-1)*5 + 5]; % Constrain rot_y end if y == 0 || y == L_y bc_nodes = [bc_nodes, (i-1)*5 + 4]; % Constrain rot_x end end end bc_nodes = unique(bc_nodes); % Apply constraints by partitioning all_dofs = 1:total_dofs; free_dofs = setdiff(all_dofs, bc_nodes); %% 6. SYSTEM SOLUTION U_free = K_global(free_dofs, free_dofs) \ F_global(free_dofs); U_global = zeros(total_dofs, 1); U_global(free_dofs) = U_free; %% 7. POST-PROCESSING & VISUALIZATION w_deflections = U_global(3:5:end); fprintf('Maximum Mid-span Deflection: %e meters\n', min(w_deflections)); % Visualizing Deflection Field X_grid = reshape(node_coords(:,1), [Nx_elem+1, Ny_elem+1]); Y_grid = reshape(node_coords(:,2), [Nx_elem+1, Ny_elem+1]); W_grid = reshape(w_deflections, [Nx_elem+1, Ny_elem+1]); figure; surf(X_grid, Y_grid, W_grid, 'EdgeColor', 'interp'); colorbar; title('Laminated Composite Plate Displacement Contour (w)'); xlabel('X Dimension (m)'); ylabel('Y Dimension (m)'); zlabel('Deflection (m)'); shading interp; view(3); %% --- HELPER FUNCTIONS --- function [A, B, D, As] = compute_ABD_matrices(E1, E2, G12, G13, G23, nu12, nu21, theta, ply_t) A = zeros(3,3); B = zeros(3,3); D = zeros(3,3); As = zeros(2,2); kappa = 5/6; % Shear correction factor % Reduced stiffness components in principal material coordinates Q11 = E1 / (1 - nu12*nu21); Q12 = (nu21*E1) / (1 - nu12*nu21); Q22 = E2 / (1 - nu12*nu21); Q66 = G12; Q44 = G23; Q55 = G13; Q = [Q11, Q12, 0; Q12, Q22, 0; 0, 0, Q66]; Q_shear = [Q44, 0; 0, Q55]; % Coordinates of ply interfaces relative to geometric mid-plane h_total = sum(ply_t); z = zeros(length(ply_t)+1, 1); z(1) = -h_total / 2; for k = 1:length(ply_t) z(k+1) = z(k) + ply_t(k); end % Accumulate matrices layer by layer for k = 1:length(theta) a = deg2rad(theta(k)); m = cos(a); n = sin(a); % Transformation matrix for in-plane components T = [m^2, n^2, 2*m*n; n^2, m^2, -2*m*n; -m*n, m*n, m^2-n^2]; R = [1 0 0; 0 1 0; 0 0 2]; % Reissner profile strain correction Q_bar = T \ Q * R * T / R; % Transformation for transverse shear components T_s = [m, -n; n, m]; Q_bar_shear = T_s \ Q_shear * T_s; % Integration over ply bounds 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); As = As + kappa * Q_bar_shear * (z(k+1) - z(k)); end end function [node_coords, elem_nodes] = generate_mesh(Lx, Ly, Nx, Ny) [X, Y] = meshgrid(linspace(0, Lx, Nx+1), linspace(0, Ly, Ny+1)); node_coords = [X(:), Y(:)]; elem_nodes = []; for j = 1:Ny for i = 1:Nx n1 = (j-1)*(Nx+1) + i; n2 = n1 + 1; n4 = j*(Nx+1) + i; n3 = n4 + 1; elem_nodes = [elem_nodes; n1, n2, n3, n4]; end end end function [N, B_b] = get_shape_functions(xi, eta, a, b) % Shape functions for 4-node rectangular bilinear element N = 0.25 * [(1-xi)*(1-eta); (1+xi)*(1-eta); (1+xi)*(1+eta); (1-xi)*(1+eta)]; % Derivatives with respect to natural coordinates xi and eta dN_dxi = 0.25 * [-(1-eta); (1-eta); (1+eta); -(1+eta)]; dN_deta = 0.25 * [-(1-xi); -(1+xi); (1+xi); (1-xi)]; % Derivatives with respect to physical coordinates x and y dN_dx = dN_dxi / a; dN_dy = dN_deta / b; B_b = zeros(6, 20); for i = 1:4 col = (i-1)*5 + 1; % Strain-displacement configuration matching A, B, D ordering B_b(1, col) = dN_dx(i); B_b(2, col+1) = dN_dy(i); B_b(3, col) = dN_dy(i); B_b(3, col+1) = dN_dx(i); B_b(4, col+3) = dN_dx(i); B_b(5, col+4) = dN_dy(i); B_b(6, col+3) = dN_dy(i); B_b(6, col+4) = dN_dx(i); end end function [N, B_s] = get_shear_B_matrix(xi, eta, a, b) N = 0.25 * [(1-xi)*(1-eta); (1+xi)*(1-eta); (1+xi)*(1+eta); (1-xi)*(1+eta)]; dN_dxi = 0.25 * [-(1-eta); (1-eta); (1+eta); -(1+eta)]; dN_deta = 0.25 * [-(1-xi); -(1+xi); (1+xi); (1-xi)]; dN_dx = dN_dxi / a; dN_dy = dN_deta / b; B_s = zeros(2, 20); for i = 1:4 col = (i-1)*5 + 1; % Transverse shear strains: gamma_yz, gamma_xz B_s(1, col+2) = dN_dy(i); B_s(1, col+4) = N(i); B_s(2, col+2) = dN_dx(i); B_s(2, col+3) = N(i); end end Use code with caution. Interpretation of Analysis and Results -1/3] * sqrt(3)

The following MATLAB script automates the calculation of the ABDcap A cap B cap D

function Ke = element_stiffness(xy, ABD, As) % xy: 4x2 matrix of element nodal coordinates (x,y) % ABD: 6x6 matrix [A B; B D] % As: 2x2 shear stiffness matrix % Returns 20x20 element stiffness matrix

Modify the theta array (e.g., [45, -45, -45, 45] ) to study angle-ply laminates. Evaluating Stresses: Once the deflection

% Gauss points for 2x2 (full) and 1x1 (reduced) gauss2 = [-1/3, 1/3; % 2x2 points and weights 1/3, 1/3; 1/3, -1/3; -1/3, -1/3] * sqrt(3); w2 = [1,1,1,1];