% Build finite difference matrix N_total = nx * ny; A_mat = sparse(N_total, N_total); F = zeros(N_total,1);
[ 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) ] We discretize the plate domain ( 0 \le x \le a ), ( 0 \le y \le b ) into ( n_x \times n_y ) grid points. 2.1 Central Difference Approximations At interior node ( (i,j) ): Composite Plate Bending Analysis With Matlab Code
for i = 2:nx-1 for j = 2:ny-1 idx = node(i,j); % Finite difference coefficients F(idx) = q0; % uniform pressure % Build finite difference matrix N_total = nx
[ \frac{\partial^4 w}{\partial x^4} \approx \frac{w_{i-2,j} - 4 w_{i-1,j} + 6 w_{i,j} - 4 w_{i+1,j} + w_{i+2,j}}{h_x^4} ] [ \frac{\partial^4 w}{\partial y^4} \approx \frac{w_{i,j-2} - 4 w_{i,j-1} + 6 w_{i,j} - 4 w_{i,j+1} + w_{i,j+2}}{h_y^4} ] [ \frac{\partial^4 w}{\partial x^2 \partial y^2} \approx \frac{w_{i+1,j+1} - 2w_{i+1,j} + w_{i+1,j-1} - 2w_{i,j+1} + 4w_{i,j} - 2w_{i,j-1} + w_{i-1,j+1} - 2w_{i-1,j} + w_{i-1,j-1}}{h_x^2 h_y^2} ] [ w = 0, \quad M_x = -D_{11} \frac{\partial^2 w}{\partial x^2} - D_{12} \frac{\partial^2 w}{\partial y^2} = 0 \quad \text{at } x=0,a ] [ w = 0, \quad M_y = -D_{12} \frac{\partial^2 w}{\partial x^2} - D_{22} \frac{\partial^2 w}{\partial y^2} = 0 \quad \text{at } y=0,b ] A_mat = sparse(N_total
[w, x, y] = CompositePlateBending(a, b, layup, thicknesses, q0, nx, ny);
[ D_{11} \frac{\partial^4 w}{\partial x^4} + 4 D_{16} \frac{\partial^4 w}{\partial x^3 \partial y} + 2(D_{12} + 2 D_{66}) \frac{\partial^4 w}{\partial x^2 \partial y^2} + 4 D_{26} \frac{\partial^4 w}{\partial x \partial y^3} + D_{22} \frac{\partial^4 w}{\partial y^4} = q(x,y) ]
% w_xxxx term if i-2 >= 1, A_mat(idx, node(i-2,j)) = A_mat(idx, node(i-2,j)) + Dxx/dx^4; end A_mat(idx, node(i-1,j)) = A_mat(idx, node(i-1,j)) -4*Dxx/dx^4; A_mat(idx, node(i,j)) = A_mat(idx, node(i,j)) +6*Dxx/dx^4; A_mat(idx, node(i+1,j)) = A_mat(idx, node(i+1,j)) -4*Dxx/dx^4; if i+2 <= nx, A_mat(idx, node(i+2,j)) = A_mat(idx, node(i+2,j)) + Dxx/dx^4; end % w_yyyy term if j-2 >= 1, A_mat(idx, node(i,j-2)) = A_mat(idx, node(i,j-2)) + Dyy/dy^4; end A_mat(idx, node(i,j-1)) = A_mat(idx, node(i,j-1)) -4*Dyy/dy^4; A_mat(idx, node(i,j)) = A_mat(idx, node(i,j)) +6*Dyy/dy^4; A_mat(idx, node(i,j+1)) = A_mat(idx, node(i,j+1)) -4*Dyy/dy^4; if j+2 <= ny, A_mat(idx, node(i,j+2)) = A_mat(idx, node(i,j+2)) + Dyy/dy^4; end % w_xxyy term coef = 2*Dxy/(dx^2 * dy^2); if i-1>=1 && j-1>=1, A_mat(idx, node(i-1,j-1)) = A_mat(idx, node(i-1,j-1)) + coef; end if i-1>=1, A_mat(idx, node(i-1,j)) = A_mat(idx, node(i-1,j)) -2*coef; end if i-1>=1 && j+1<=ny, A_mat(idx, node(i-1,j+1)) = A_mat(idx, node(i-1,j+1)) + coef; end if j-1>=1, A_mat(idx, node(i,j-1)) = A_mat(idx, node(i,j-1)) -2*coef; end A_mat(idx, idx) = A_mat(idx, idx) +4*coef; if j+1<=ny, A_mat(idx, node(i,j+1)) = A_mat(idx, node(i,j+1)) -2*coef; end if i+1<=nx && j-1>=1, A_mat(idx, node(i+1,j-1)) = A_mat(idx, node(i+1,j-1)) + coef; end if i+1<=nx, A_mat(idx, node(i+1,j)) = A_mat(idx, node(i+1,j)) -2*coef; end if i+1<=nx && j+1<=ny, A_mat(idx, node(i+1,j+1)) = A_mat(idx, node(i+1,j+1)) + coef; end end end