clear all % lx = 1; ly = 1; syms x y lx ly N(1) = 1/4*(1-x/lx)*(1-y/ly); N(2) = 1/4*(1+x/lx)*(1-y/ly); N(3) = 1/4*(1+x/lx)*(1+y/ly); N(4) = 1/4*(1+x/lx)*(1+y/ly); Be= [diff(N(1),x) 0 diff(N(2),x) 0 diff(N(3),x) 0 diff(N(4),x) 0; 0 diff(N(1),y) 0 diff(N(2),y) 0 diff(N(3),y) 0 diff(N(4),y); diff(N(1),y) diff(N(1),x) diff(N(2),y) diff(N(2),x) diff(N(3),y) diff(N(3),x) diff(N(4),y) diff(N(4),x)]; % n5 = 1/2 * (1 - (x/lx)^2) * (1 - y/ly); % n6 = 1/2 * (1 + (x/lx)) * (1 - (y/ly)^2); % n7 = 1/2 * (1 - (x/lx)^2) * (1 + (y/ly)); % n8 = 1/2 * (1 - (x/lx)) * (1 - (y/ly)^2); % n1 = N(1) - 1/2*n5 - 1/2*n8; % n2 = N(2) - 1/2*n5 - 1/2*n6; % n3 = N(3) - 1/2*n6 - 1/2*n7; % n4 = N(4) - 1/2*n7 - 1/2*n8; % % Be= [diff(n1,x), 0 , diff(n2,x), 0 , diff(n3,x), 0 , diff(n4,x), 0 , diff(n5,x), 0 , diff(n6,x), 0 , diff(n7,x), 0 , diff(n8,x), 0; % % 0 , diff(n1,y), 0 , diff(n2,y), 0 , diff(n3,y), 0 , diff(n4,y), 0 , diff(n5,y), 0 , diff(n6,y), 0 , diff(n7,y), 0 , diff(n8,y); % % diff(n1,y), diff(n1,x), diff(n2,y), diff(n2,x), diff(n3,y), diff(n3,x), diff(n4,y), diff(n4,x), diff(n5,y), diff(n5,x), diff(n6,y), diff(n6,x), diff(n7,y), diff(n7,x), diff(n8,y), diff(n8,x)]; % Be= [diff(n1,x), 0 , diff(n5,x), 0 , diff(n2,x), 0 , diff(n6,x), 0 , diff(n3,x), 0 , diff(n7,x), 0 , diff(n4,x), 0 , diff(n8,x), 0; % 0 , diff(n1,y), 0 , diff(n5,y), 0 , diff(n2,y), 0 , diff(n6,y), 0 , diff(n3,y), 0 , diff(n7,y), 0 , diff(n4,y), 0 , diff(n8,y); % diff(n1,y), diff(n1,x), diff(n5,y), diff(n5,x), diff(n2,y), diff(n2,x), diff(n6,y), diff(n6,x), diff(n3,y), diff(n3,x), diff(n7,y), diff(n7,x), diff(n4,y), diff(n4,x), diff(n8,y), diff(n8,x)]; % syms E mu % syms d11 d12 d13 d22 d23 d33 % D = [d11 d12 d13 % d12 d22 d23 % d13 d23 d33]; E=1; mu=0.3; nu = 0.3; D=(E/(1-mu^2))*[1 mu 0 mu 1 0 0 0 (1-mu)/2]; % E = 0.0387 % syms u1 u2 u3 u4 u5 u6 u7 u8 % u = [u1 u2 u3 u4 u5 u6 u7 u8] % lx = 1; ly = 1; syms x y lx ly N(1)=(1-x/lx)*(1-y/ly); N(2)=x/lx*(1-y/ly); N(3)=x/lx*y/ly; N(4)=(1-x/lx)*y/ly; Be= [diff(N(1),x) 0 diff(N(2),x) 0 diff(N(3),x) 0 diff(N(4),x) 0; 0 diff(N(1),y) 0 diff(N(2),y) 0 diff(N(3),y) 0 diff(N(4),y); diff(N(1),y) diff(N(1),x) diff(N(2),y) diff(N(2),x) diff(N(3),y) diff(N(3),x) diff(N(4),y) diff(N(4),x)]; % K=Be.'*D*Be; % Ke=sym(zeros(8,8)); K=Be; He = sym(zeros(3,8)); for i=1:3 for j=1:8 He(i,j) = int(int(K(i,j),'x',0.0,1),'y',0.0,1); % changing % Ke(i,j)=int(int(K(i,j),'x',[0,lx]),'y',[0,ly]); end end He vpa(Ke,4) % B divide 2 Se = ... [- d11/2 - d13/4, - d12/2 - d13/4, d11/2 - d13/4, d13/4 - d12/2, d11/2 + d13/4, d12/2 + d13/4, d13/4 - d11/2, d12/2 - d13/4 - d12/2 - d23/4, - d22/2 - d23/4, d12/2 - d23/4, d23/4 - d22/2, d12/2 + d23/4, d22/2 + d23/4, d23/4 - d12/2, d22/2 - d23/4 - d13/2 - d33/4, - d23/2 - d33/4, d13/2 - d33/4, d33/4 - d23/2, d13/2 + d33/4, d23/2 + d33/4, d33/4 - d13/2, d23/2 - d33/4;]; % a = 1; b = 1; % temp = E/(a*b*(1-nu^2))*[1/3*(b^2+(1-nu)/2*a^2), a*b/8*(1+nu), -1/3*(b^2-(1-nu)/4*a^2), -a*b/8*(1-3*nu), -1/6*(b^2+(1-nu)/2*a^2), -a*b/8*(1+nu), 1/6*(b^2-(1-nu)*a^2), a*b/8*(1-3*nu); % a*b/8*(1+nu), 1/3*(a^2+(1-nu)/2*b^2), a*b/8*(1-3*nu), -1/6*(a^2-(1-nu)*b^2), -a*b/8*(1+nu), -1/6*(a^2+(1-nu)/2*b^2), -a*b/8*(1-3*nu), 1/3*(-a^2+(1-nu)/4*b^2); % -1/3*(b^2-(1-nu)/4*a^2), a*b/8*(1-3*nu), 1/3*(b^2+(1-nu)/2*a^2), -a*b/8*(1+nu), 1/6*(b^2-(1-nu)*a^2), -a*b/8*(1-3*nu), -1/6*(b^2+(1-nu)/2*a^2), a*b/8*(1+nu); % -a*b/8*(1-3*nu), -1/6*(a^2-(1-nu)*b^2), -a*b/8*(1+nu), 1/3*(a^2+(1-nu)/2*b^2), a*b/8*(1-3*nu), 1/3*(-a^2+(1-nu)/4*b^2), a*b/8*(1+nu), -1/6*(a^2+(1-nu)/2*b^2); % -1/6*(b^2+(1-nu)/2*a^2), -a*b/8*(1+nu), 1/6*(b^2-(1-nu)*a^2), a*b/8*(1-3*nu), 1/3*(b^2+(1-nu)/2*a^2), a*b/8*(1+nu), 1/3*(-b^2+(1-nu)/4*a^2), -a*b/8*(1-3*nu); % -a*b/8*(1+nu), -1/6*(a^2+(1-nu)/2*b^2), -a*b/8*(1-3*nu), 1/3*(-a^2+(1-nu)/4*b^2), a*b/8*(1+nu), 1/3*(a^2+(1-nu)/2*b^2), a*b/8*(1-3*nu), 1/6*(a^2-(1-nu)*b^2); % 1/6*(b^2-(1-nu)*a^2), -a*b/8*(1-3*nu), -1/6*(b^2+(1-nu)/2*a^2), a*b/8*(1+nu), 1/3*(-b^2+(1-nu)/4*a^2), a*b/8*(1-3*nu), 1/3*(b^2+(1-nu)/2*a^2), -a*b/8*(1+nu); % a*b/8*(1-3*nu), 1/3*(-a^2+(1-nu)/4*b^2), a*b/8*(1+nu), -1/6*(a^2+(1-nu)/2*b^2), -a*b/8*(1-3*nu), 1/6*(a^2-(1-nu)*b^2), -a*b/8*(1+nu), 1/3*(a^2+(1-nu)/2*b^2);] % K = vpa(Ke,4); % B divide 2 K = ... [ 0.3333*d11 + 0.25*d13 + 0.08333*d33, 0.25*d12 + 0.1667*d13 + 0.1667*d23 + 0.0625*d33, 0.04167*d33 - 0.3333*d11, 0.25*d12 - 0.1667*d13 + 0.08333*d23 - 0.0625*d33, - 0.1667*d11 - 0.25*d13 - 0.04167*d33, - 0.25*d12 - 0.08333*d13 - 0.08333*d23 - 0.0625*d33, 0.1667*d11 - 0.08333*d33, 0.08333*d13 - 0.25*d12 - 0.1667*d23 + 0.0625*d33; 0.25*d12 + 0.1667*d13 + 0.1667*d23 + 0.0625*d33, 0.3333*d22 + 0.25*d23 + 0.08333*d33, 0.08333*d23 - 0.1667*d13 - 0.25*d12 + 0.0625*d33, 0.1667*d22 - 0.08333*d33, - 0.25*d12 - 0.08333*d13 - 0.08333*d23 - 0.0625*d33, - 0.1667*d22 - 0.25*d23 - 0.04167*d33, 0.25*d12 + 0.08333*d13 - 0.1667*d23 - 0.0625*d33, 0.04167*d33 - 0.3333*d22; 0.04167*d33 - 0.3333*d11, 0.08333*d23 - 0.1667*d13 - 0.25*d12 + 0.0625*d33, 0.3333*d11 - 0.25*d13 + 0.08333*d33, 0.1667*d13 - 0.25*d12 + 0.1667*d23 - 0.0625*d33, 0.1667*d11 - 0.08333*d33, 0.25*d12 + 0.08333*d13 - 0.1667*d23 - 0.0625*d33, 0.25*d13 - 0.1667*d11 - 0.04167*d33, 0.25*d12 - 0.08333*d13 - 0.08333*d23 + 0.0625*d33; 0.25*d12 - 0.1667*d13 + 0.08333*d23 - 0.0625*d33, 0.1667*d22 - 0.08333*d33, 0.1667*d13 - 0.25*d12 + 0.1667*d23 - 0.0625*d33, 0.3333*d22 - 0.25*d23 + 0.08333*d33, 0.08333*d13 - 0.25*d12 - 0.1667*d23 + 0.0625*d33, 0.04167*d33 - 0.3333*d22, 0.25*d12 - 0.08333*d13 - 0.08333*d23 + 0.0625*d33, 0.25*d23 - 0.1667*d22 - 0.04167*d33; - 0.1667*d11 - 0.25*d13 - 0.04167*d33, - 0.25*d12 - 0.08333*d13 - 0.08333*d23 - 0.0625*d33, 0.1667*d11 - 0.08333*d33, 0.08333*d13 - 0.25*d12 - 0.1667*d23 + 0.0625*d33, 0.3333*d11 + 0.25*d13 + 0.08333*d33, 0.25*d12 + 0.1667*d13 + 0.1667*d23 + 0.0625*d33, 0.04167*d33 - 0.3333*d11, 0.25*d12 - 0.1667*d13 + 0.08333*d23 - 0.0625*d33; -0.25*d12 - 0.08333*d13 - 0.08333*d23 - 0.0625*d33, - 0.1667*d22 - 0.25*d23 - 0.04167*d33, 0.25*d12 + 0.08333*d13 - 0.1667*d23 - 0.0625*d33, 0.04167*d33 - 0.3333*d22, 0.25*d12 + 0.1667*d13 + 0.1667*d23 + 0.0625*d33, 0.3333*d22 + 0.25*d23 + 0.08333*d33, 0.08333*d23 - 0.1667*d13 - 0.25*d12 + 0.0625*d33, 0.1667*d22 - 0.08333*d33; 0.1667*d11 - 0.08333*d33, 0.25*d12 + 0.08333*d13 - 0.1667*d23 - 0.0625*d33, 0.25*d13 - 0.1667*d11 - 0.04167*d33, 0.25*d12 - 0.08333*d13 - 0.08333*d23 + 0.0625*d33, 0.04167*d33 - 0.3333*d11, 0.08333*d23 - 0.1667*d13 - 0.25*d12 + 0.0625*d33, 0.3333*d11 - 0.25*d13 + 0.08333*d33, 0.1667*d13 - 0.25*d12 + 0.1667*d23 - 0.0625*d33; 0.08333*d13 - 0.25*d12 - 0.1667*d23 + 0.0625*d33, 0.04167*d33 - 0.3333*d22, 0.25*d12 - 0.08333*d13 - 0.08333*d23 + 0.0625*d33, 0.25*d23 - 0.1667*d22 - 0.04167*d33, 0.25*d12 - 0.1667*d13 + 0.08333*d23 - 0.0625*d33, 0.1667*d22 - 0.08333*d33, 0.1667*d13 - 0.25*d12 + 0.1667*d23 - 0.0625*d33, 0.3333*d22 - 0.25*d23 + 0.08333*d33;]; % B not divide 2 K = ... [ d11/3 + d13/2 + d33/3, d12/4 + d13/3 + d23/3 + d33/4, d33/6 - d11/3, d12/4 - d13/3 + d23/6 - d33/4, - d11/6 - d13/2 - d33/6, - d12/4 - d13/6 - d23/6 - d33/4, d11/6 - d33/3, d13/6 - d12/4 - d23/3 + d33/4; d12/4 + d13/3 + d23/3 + d33/4, d22/3 + d23/2 + d33/3, d23/6 - d13/3 - d12/4 + d33/4, d22/6 - d33/3, - d12/4 - d13/6 - d23/6 - d33/4, - d22/6 - d23/2 - d33/6, d12/4 + d13/6 - d23/3 - d33/4, d33/6 - d22/3; d33/6 - d11/3, d23/6 - d13/3 - d12/4 + d33/4, d11/3 - d13/2 + d33/3, d13/3 - d12/4 + d23/3 - d33/4, d11/6 - d33/3, d12/4 + d13/6 - d23/3 - d33/4, d13/2 - d11/6 - d33/6, d12/4 - d13/6 - d23/6 + d33/4; d12/4 - d13/3 + d23/6 - d33/4, d22/6 - d33/3, d13/3 - d12/4 + d23/3 - d33/4, d22/3 - d23/2 + d33/3, d13/6 - d12/4 - d23/3 + d33/4, d33/6 - d22/3, d12/4 - d13/6 - d23/6 + d33/4, d23/2 - d22/6 - d33/6; - d11/6 - d13/2 - d33/6, - d12/4 - d13/6 - d23/6 - d33/4, d11/6 - d33/3, d13/6 - d12/4 - d23/3 + d33/4, d11/3 + d13/2 + d33/3, d12/4 + d13/3 + d23/3 + d33/4, d33/6 - d11/3, d12/4 - d13/3 + d23/6 - d33/4; - d12/4 - d13/6 - d23/6 - d33/4, - d22/6 - d23/2 - d33/6, d12/4 + d13/6 - d23/3 - d33/4, d33/6 - d22/3, d12/4 + d13/3 + d23/3 + d33/4, d22/3 + d23/2 + d33/3, d23/6 - d13/3 - d12/4 + d33/4, d22/6 - d33/3; d11/6 - d33/3, d12/4 + d13/6 - d23/3 - d33/4, d13/2 - d11/6 - d33/6, d12/4 - d13/6 - d23/6 + d33/4, d33/6 - d11/3, d23/6 - d13/3 - d12/4 + d33/4, d11/3 - d13/2 + d33/3, d13/3 - d12/4 + d23/3 - d33/4; d13/6 - d12/4 - d23/3 + d33/4, d33/6 - d22/3, d12/4 - d13/6 - d23/6 + d33/4, d23/2 - d22/6 - d33/6, d12/4 - d13/3 + d23/6 - d33/4, d22/6 - d33/3, d13/3 - d12/4 + d23/3 - d33/4, d22/3 - d23/2 + d33/3]; % Ke = [ (E*(mu - 9))/(24*(mu^2 - 1)), -(E*(7*mu + 1))/(32*(mu^2 - 1)), (E*(mu + 15))/(48*(mu^2 - 1)), -(E*(9*mu - 1))/(32*(mu^2 - 1)), -(E*(mu - 9))/(48*(mu^2 - 1)), (E*(7*mu + 1))/(32*(mu^2 - 1)), -(E*(mu + 3))/(24*(mu^2 - 1)), (E*(9*mu - 1))/(32*(mu^2 - 1)) % -(E*(7*mu + 1))/(32*(mu^2 - 1)), (E*(mu - 9))/(24*(mu^2 - 1)), (E*(9*mu - 1))/(32*(mu^2 - 1)), -(E*(mu + 3))/(24*(mu^2 - 1)), (E*(7*mu + 1))/(32*(mu^2 - 1)), -(E*(mu - 9))/(48*(mu^2 - 1)), -(E*(9*mu - 1))/(32*(mu^2 - 1)), (E*(mu + 15))/(48*(mu^2 - 1)) % (E*(mu + 15))/(48*(mu^2 - 1)), (E*(9*mu - 1))/(32*(mu^2 - 1)), (E*(mu - 9))/(24*(mu^2 - 1)), (E*(7*mu + 1))/(32*(mu^2 - 1)), -(E*(mu + 3))/(24*(mu^2 - 1)), -(E*(9*mu - 1))/(32*(mu^2 - 1)), -(E*(mu - 9))/(48*(mu^2 - 1)), -(E*(7*mu + 1))/(32*(mu^2 - 1)) % -(E*(9*mu - 1))/(32*(mu^2 - 1)), -(E*(mu + 3))/(24*(mu^2 - 1)), (E*(7*mu + 1))/(32*(mu^2 - 1)), (E*(mu - 9))/(24*(mu^2 - 1)), (E*(9*mu - 1))/(32*(mu^2 - 1)), (E*(mu + 15))/(48*(mu^2 - 1)), -(E*(7*mu + 1))/(32*(mu^2 - 1)), -(E*(mu - 9))/(48*(mu^2 - 1)) % -(E*(mu - 9))/(48*(mu^2 - 1)), (E*(7*mu + 1))/(32*(mu^2 - 1)), -(E*(mu + 3))/(24*(mu^2 - 1)), (E*(9*mu - 1))/(32*(mu^2 - 1)), (E*(mu - 9))/(24*(mu^2 - 1)), -(E*(7*mu + 1))/(32*(mu^2 - 1)), (E*(mu + 15))/(48*(mu^2 - 1)), -(E*(9*mu - 1))/(32*(mu^2 - 1)) % (E*(7*mu + 1))/(32*(mu^2 - 1)), -(E*(mu - 9))/(48*(mu^2 - 1)), -(E*(9*mu - 1))/(32*(mu^2 - 1)), (E*(mu + 15))/(48*(mu^2 - 1)), -(E*(7*mu + 1))/(32*(mu^2 - 1)), (E*(mu - 9))/(24*(mu^2 - 1)), (E*(9*mu - 1))/(32*(mu^2 - 1)), -(E*(mu + 3))/(24*(mu^2 - 1)) % -(E*(mu + 3))/(24*(mu^2 - 1)), -(E*(9*mu - 1))/(32*(mu^2 - 1)), -(E*(mu - 9))/(48*(mu^2 - 1)), -(E*(7*mu + 1))/(32*(mu^2 - 1)), (E*(mu + 15))/(48*(mu^2 - 1)), (E*(9*mu - 1))/(32*(mu^2 - 1)), (E*(mu - 9))/(24*(mu^2 - 1)), (E*(7*mu + 1))/(32*(mu^2 - 1)) % (E*(9*mu - 1))/(32*(mu^2 - 1)), (E*(mu + 15))/(48*(mu^2 - 1)), -(E*(7*mu + 1))/(32*(mu^2 - 1)), -(E*(mu - 9))/(48*(mu^2 - 1)), -(E*(9*mu - 1))/(32*(mu^2 - 1)), -(E*(mu + 3))/(24*(mu^2 - 1)), (E*(7*mu + 1))/(32*(mu^2 - 1)), (E*(mu - 9))/(24*(mu^2 - 1))]; % Ke % syms E mu % ke with B is original (not divide 2) % ke = [(E*(mu - 3)) / (6 * (mu*mu - 1)), -E / (8 * (mu - 1)), (E*(mu + 3)) / (12 * (mu*mu - 1)), -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)), -(E*(mu - 3)) / (12 * (mu*mu - 1)), E / (8 * (mu - 1)), -(E*mu) / (6 * (mu*mu - 1)), (E*(3 * mu - 1)) / (8 * (mu*mu - 1)) % -E / (8 * (mu - 1)), (E*(mu - 3)) / (6 * (mu*mu - 1)), (E*(3 * mu - 1)) / (8 * (mu*mu - 1)), -(E*mu) / (6 * (mu*mu - 1)), E / (8 * (mu - 1)), -(E*(mu - 3)) / (12 * (mu*mu - 1)), -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)), (E*(mu + 3)) / (12 * (mu*mu - 1)) % (E*(mu + 3)) / (12 * (mu*mu - 1)), (E*(3 * mu - 1)) / (8 * (mu*mu - 1)), (E*(mu - 3)) / (6 * (mu*mu - 1)), E / (8 * (mu - 1)), -(E*mu) / (6 * (mu*mu - 1)), -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)), -(E*(mu - 3)) / (12 * (mu*mu - 1)), -E / (8 * (mu - 1)) % -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)), -(E*mu) / (6 * (mu*mu - 1)), E / (8 * (mu - 1)), (E*(mu - 3)) / (6 * (mu*mu - 1)), (E*(3 * mu - 1)) / (8 * (mu*mu - 1)), (E*(mu + 3)) / (12 * (mu*mu - 1)), -E / (8 * (mu - 1)), -(E*(mu - 3)) / (12 * (mu*mu - 1)) % -(E*(mu - 3)) / (12 * (mu*mu - 1)), E / (8 * (mu - 1)), -(E*mu) / (6 * (mu*mu - 1)), (E*(3 * mu - 1)) / (8 * (mu*mu - 1)), (E*(mu - 3)) / (6 * (mu*mu - 1)), -E / (8 * (mu - 1)), (E*(mu + 3)) / (12 * (mu*mu - 1)), -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)) % E / (8 * (mu - 1)), -(E*(mu - 3)) / (12 * (mu*mu - 1)), -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)), (E*(mu + 3)) / (12 * (mu*mu - 1)), -E / (8 * (mu - 1)), (E*(mu - 3)) / (6 * (mu*mu - 1)), (E*(3 * mu - 1)) / (8 * (mu*mu - 1)), -(E*mu) / (6 * (mu*mu - 1)) % -(E*mu) / (6 * (mu*mu - 1)), -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)), -(E*(mu - 3)) / (12 * (mu*mu - 1)), -E / (8 * (mu - 1)), (E*(mu + 3)) / (12 * (mu*mu - 1)), (E*(3 * mu - 1)) / (8 * (mu*mu - 1)), (E*(mu - 3)) / (6 * (mu*mu - 1)), E / (8 * (mu - 1)) % (E*(3 * mu - 1)) / (8 * (mu*mu - 1)), (E*(mu + 3)) / (12 * (mu*mu - 1)), -E / (8 * (mu - 1)), -(E*(mu - 3)) / (12 * (mu*mu - 1)), -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)), -(E*mu) / (6 * (mu*mu - 1)), E / (8 * (mu - 1)), (E*(mu - 3)) / (6 * (mu*mu - 1)) ]