Auiui = (26*h)/35 Aduidui = (2*h^3)/105 Auidui = 0 Aui1ui = (9*h)/70 Adui1dui = -h^3/140 Adui1ui = (13*h^2)/420 Aui1dui = -(13*h^2)/420 Buiui = 12/(5*h) Bduidui = (4*h)/15 Buidui = 0 Bui1ui = -6/(5*h) Bdui1dui = -h/30 Bdui1ui = -1/10 Bui1dui = 1/10 +- -+ | #1 - #2 + 1 | | | | 2 3 | | 2 conj(x) conj(x) | | conj(x) - ---------- + -------- | | conj(h) 2 | | conj(h) | | | | #2 - #1 | | | | 3 2 | | conj(x) conj(x) | | -------- - -------- | | 2 conj(h) | | conj(h) | +- -+ where 3 2 conj(x) #1 == ---------- 3 conj(h) 2 3 conj(x) #2 == ---------- 2 conj(h) A = [ 1200000/h^3, 600000/h^2, -1200000/h^3, 600000/h^2] [ 600000/h^2, 400000/h, -600000/h^2, 200000/h] [ -1200000/h^3, -600000/h^2, 1200000/h^3, -600000/h^2] [ 600000/h^2, 200000/h, -600000/h^2, 400000/h] b = [ -h/2, -h^2/12, -h/2, h^2/12] >>
%% Meng Yee CHUAH % Week 6 - Finite Elements %% 9.1 (c) close all; clear all; clc; syms x h; cubicMatrix = [1 0 0 0; 0 1 0 0; 1 h h^2 h^3; 0 1 2*h 3*h^2]; u = [1 x x^2 x^3]'; phi = (cubicMatrix\eye(4))'*u; dphi = diff(phi,x); Auiui = int(phi(1)^2,x,0,h)+int(phi(3)^2,x,0,h) Aduidui = int(phi(2)^2,x,0,h)+int(phi(4)^2,x,0,h) Auidui = int(phi(1)*phi(2),x,0,h)+int(phi(3)*phi(4),x,0,h) Aui1ui = int(phi(1)*phi(3),x,0,h) Adui1dui = int(phi(2)*phi(4),x,0,h) Adui1ui = int(phi(2)*phi(3),x,0,h) Aui1dui = int(phi(1)*phi(4),x,0,h) Buiui = int(dphi(1)^2,x,0,h)+int(dphi(3)^2,x,0,h) Bduidui = int(dphi(2)^2,x,0,h)+int(dphi(4)^2,x,0,h) Buidui = int(dphi(1)*dphi(2),x,0,h)+int(dphi(3)*dphi(4),x,0,h) Bui1ui = int(dphi(1)*dphi(3),x,0,h) Bdui1dui = int(dphi(2)*dphi(4),x,0,h) Bdui1ui = int(dphi(2)*dphi(3),x,0,h) Bui1dui = int(dphi(1)*dphi(4),x,0,h) pretty(phi) %% 9.2 Beam bending E = 1e4; % Young's modulus I = 10; % Moment of inertia nElements = 10; % Number of elements L = 100; % Length of beam nNodes = 2*(nElements+1); % Distributed load at each element % f = 0; f = -1; % Global load, F = [force, moment, force, moment, ...] F = zeros(nNodes,1); % F(end-1) = -1e2; % Tip force d2phi = diff(phi,x,2); for j = 1:4 for i = 1:4 A(i,j) = int(E*I*d2phi(i)*d2phi(j),x,0,h); end b(j) = int(phi(j)*f,x,0,h); end A % Element stiffness matrix b % Element force vector % Elements of length h = 1 A = subs(A,h,L/nElements); b = subs(b',h,L/nElements); K = zeros(nNodes); r = zeros(nNodes,1); for i = 1:2:nNodes-2 % Form global stiffness matrix K K(i:i+3,i:i+3) = K(i:i+3,i:i+3)+A; % Form global force vector r r(i:i+3) = r(i:i+3)+b; end r = r+F; % Find displacement vector u = K\r % For boundary condition u0 = 0 and udot0 = 0 (Cantilevered beam) u = K(3:end,3:end)\r(3:end); x = 0:L/nElements:L; % Exact solution for tip force y = 1e2/(6*E*I)*(x.^3-3*L*x.^2); % Exact solution for distributed force y2 = -1/(24*E*I)*(x.^4-4*L*x.^3+6*(L^2)*x.^2); figure('Units','inches','Position',[0 0 6 4.5],'PaperPosition',[0 0 6 4.5]); plot([0; u(1:2:end)]); hold on; plot(y2,'-r'); xlim([0 nElements]); legend('Finite element solution','True solution','Location','Best'); xlabel('Number of elements'); ylabel('Vertical displacement'); title('Beam bending under distributed load'); saveas(gcf,sprintf('9-2 Beam bending under distributed load 2.png')); saveas(gcf,sprintf('9-2 Beam bending under distributed load 2.fig'));