Composite Plate: Bending Analysis With Matlab Code

For interior node (i,j):

%% Compute ABD Matrix A = zeros(3,3); B = zeros(3,3); D = zeros(3,3); for k = 1:num_plies theta_k = theta(k) * pi/180; m = cos(theta_k); n = sin(theta_k); % Transformation matrix T = [m^2, n^2, 2 m n; n^2, m^2, -2 m n; -m n, m n, m^2-n^2]; % Q_bar = T * Q * T_inv Q = [Q11, Q12, 0; Q12, Q22, 0; 0, 0, Q66]; Q_bar = T * Q * T'; % Integrate through thickness A = A + Q_bar * (z(k+1)-z(k)); B = B + Q_bar * 0.5 * (z(k+1)^2 - z(k)^2); D = D + Q_bar * (1/3) * (z(k+1)^3 - z(k)^3); end % For symmetric laminate, B should be zero (numerically small) B = zeros(3,3); % enforce symmetry Composite Plate Bending Analysis With Matlab Code

% Load (uniform pressure) F(n) = 1000; % Pa end end For interior node (i,j): %% Compute ABD Matrix

[ D_{11} \frac{\partial^4 w}{\partial x^4} + 2(D_{12}+2D_{66}) \frac{\partial^4 w}{\partial x^2 \partial y^2} + D_{22} \frac{\partial^4 w}{\partial y^4} = q(x,y) ] % Central difference coefficients c1 = D(1,1)/dx^4; c2

We assemble a sparse linear system ( [K] {w} = {f} ) and solve. Below is the complete code. It computes deflections, curvatures, and then stresses in each ply at Gauss points.

% Central difference coefficients c1 = D(1,1)/dx^4; c2 = (2*(D(1,2)+2 D(3,3)))/(dx^2 dy^2); c3 = D(2,2)/dy^4;

% Build coefficient matrix for D11 w,xxxx + 2(D12+2D66) w,xxyy + D22 w,yyyy = q N = Nx*Ny; K = sparse(N,N); F = zeros(N,1);