% Analytical solution (Navier, simply supported, symmetric laminate) % For square plate a=b, load p=p0, D11, D22, D12, D66, D16=D26=0 if symmetric balanced % Assume D16=0, D26=0 for cross-ply [0/90]s D11 = D(1,1); D12 = D(1,2); D22 = D(2,2); D66 = D(3,3); w_analytical = 0; for m = 1:2:25 for n = 1:2:25 denom = pi^4 * ( D11*(m/a)^4 + 2*(D12+2 D66) (m/a)^2*(n/b)^2 + D22*(n/b)^4 ); w_analytical = w_analytical + (16 p0/(m n pi^2)) / denom; end end w_analytical = w_analytical * sin(m pi/2) sin(n pi/2); % at center
Mxx ; Myy ; Mxy = [D] * κxx ; κyy ; κxy We use a 4-node rectangular element (size 2a×2b in local coordinates). Each node has 3 DOF: w, θx = ∂w/∂y, θy = -∂w/∂x. 2.1 Shape Functions (non-conforming but widely used) The deflection w is approximated by a 12-term polynomial: Composite Plate Bending Analysis With Matlab Code
% Dummy B (3x12) - replace with actual derivatives in real code B = zeros(3,12); % B matrix structure: row1: d2w/dx2, row2: d2w/dy2, row3: 2*d2w/dxdy % For actual implementation, please refer to standard FEA textbooks. % Assemble into global matrix dof_map = zeros(1,12);
% Assemble into global matrix dof_map = zeros(1,12); for inode = 1:4 global_node = nodes(inode); dof_map(3*(inode-1)+1) = 3*(global_node-1) + 1; % w dof_map(3*(inode-1)+2) = 3*(global_node-1) + 2; % theta_x dof_map(3*(inode-1)+3) = 3*(global_node-1) + 3; % theta_y end K_global(dof_map, dof_map) = K_global(dof_map, dof_map) + Ke; F_global(dof_map) = F_global(dof_map) + Fe; end dof_map(3*(inode-1)+1) = 3*(global_node-1) + 1
% Local coordinates x = xi * a_elem; y = eta * b_elem;
% Find center deflection center_x = floor(nx/2)+1; center_y = floor(ny/2)+1; w_center_FEM = W(center_x, center_y);