%================================= % get global stiffness matrix %================================= function sK = getK(E,opt) switch opt case 1 % B / 2, length of cell = 1 % E = [E1111 E1122 sqrt(2)*E1112 % E2222 sqrt(2)*E2212 % sym: 2E1212] % syms d11 d12 d13 d22 d23 d33 % D = [d11; d12; d13; d22; d23; d33]; % d11 d12 d13 d22 d23 d33 keCoeff = [1/3 0 1/4 0 0 1/12; 0 1/4 1/6 0 1/6 1/16; -1/3 0 0 0 0 1/24; 0 1/4 -1/6 0 1/12 -1/16; -1/6 0 -1/4 0 0 -1/24; 0 -1/4 -1/12 0 -1/12 -1/16; 1/6 0 0 0 0 -1/12; 0 -1/4 1/12 0 -1/6 1/16; % 8 0 1/4 1/6 0 1/6 1/16; 0 0 0 1/3 1/4 1/12; 0 -1/4 -1/6 0 1/12 1/16; 0 0 0 1/6 0 -1/12; 0 -1/4 -1/12 0 -1/12 -1/16; 0 0 0 -1/6 -1/4 -1/24; 0 1/4 1/12 0 -1/6 -1/16; 0 0 0 -1/3 0 1/24; % 8 -1/3 0 0 0 0 1/24; 0 -1/4 -1/6 0 1/12 1/16; 1/3 0 -1/4 0 0 1/12; 0 -1/4 1/6 0 1/6 -1/16; 1/6 0 0 0 0 -1/12; 0 1/4 1/12 0 -1/6 -1/16; -1/6 0 1/4 0 0 -1/24; 0 1/4 -1/12 0 -1/12 1/16; % 8 0 1/4 -1/6 0 1/12 -1/16 0 0 0 1/6 0 -1/12; 0 -1/4 1/6 0 1/6 -1/16; 0 0 0 1/3 -1/4 1/12; 0 -1/4 1/12 0 -1/6 1/16; 0 0 0 -1/3 0 1/24; 0 1/4 -1/12 0 -1/12 1/16; 0 0 0 -1/6 1/4 -1/24; % 8 -1/6 0 -1/4 0 0 -1/24; 0 -1/4 -1/12 0 -1/12 -1/16; 1/6 0 0 0 0 -1/12; 0 -1/4 1/12 0 -1/6 1/16; 1/3 0 1/4 0 0 1/12; 0 1/4 1/6 0 1/6 1/16; -1/3 0 0 0 0 1/24; 0 1/4 -1/6 0 1/12 -1/16; % 8 0 -1/4 -1/12 0 -1/12 -1/16; 0 0 0 -1/6 -1/4 -1/24; 0 1/4 1/12 0 -1/6 -1/16; 0 0 0 -1/3 0 1/24; 0 1/4 1/6 0 1/6 1/16; 0 0 0 1/3 1/4 1/12; 0 -1/4 -1/6 0 1/12 1/16; 0 0 0 1/6 0 -1/12; % 8 1/6 0 0 0 0 -1/12; 0 1/4 1/12 0 -1/6 -1/16; -1/6 0 1/4 0 0 -1/24; 0 1/4 -1/12 0 -1/12 1/16; -1/3 0 0 0 0 1/24; 0 -1/4 -1/6 0 1/12 1/16; 1/3 0 -1/4 0 0 1/12; 0 -1/4 1/6 0 1/6 -1/16; % 8 0 -1/4 1/12 0 -1/6 1/16; 0 0 0 -1/3 0 1/24; 0 1/4 -1/12 0 -1/12 1/16; 0 0 0 -1/6 1/4 -1/24; 0 1/4 -1/6 0 1/12 -1/16; 0 0 0 1/6 0 -1/12; 0 -1/4 1/6 0 1/6 -1/16; 0 0 0 1/3 -1/4 1/12;]; K = keCoeff*E; sK = K(:); % aaa = reshape(KE,8,8); % % K = [ d11/3 + d13/4 + d33/12, d12/4 + d13/6 + d23/6 + d33/16, d33/24 - d11/3, d12/4 - d13/6 + d23/12 - d33/16, - d11/6 - d13/4 - d33/24, - d12/4 - d13/12 - d23/12 - d33/16, d11/6 - d33/12, d13/12 - d12/4 - d23/6 + d33/16; % d12/4 + d13/6 + d23/6 + d33/16, d22/3 + d23/4 + d33/12, d23/12 - d13/6 - d12/4 + d33/16, d22/6 - d33/12, - d12/4 - d13/12 - d23/12 - d33/16, - d22/6 - d23/4 - d33/24, d12/4 + d13/12 - d23/6 - d33/16, d33/24 - d22/3; % d33/24 - d11/3, d23/12 - d13/6 - d12/4 + d33/16, d11/3 - d13/4 + d33/12, d13/6 - d12/4 + d23/6 - d33/16, d11/6 - d33/12, d12/4 + d13/12 - d23/6 - d33/16, d13/4 - d11/6 - d33/24, d12/4 - d13/12 - d23/12 + d33/16; % d12/4 - d13/6 + d23/12 - d33/16, d22/6 - d33/12, d13/6 - d12/4 + d23/6 - d33/16, d22/3 - d23/4 + d33/12, d13/12 - d12/4 - d23/6 + d33/16, d33/24 - d22/3, d12/4 - d13/12 - d23/12 + d33/16, d23/4 - d22/6 - d33/24; % - d11/6 - d13/4 - d33/24, - d12/4 - d13/12 - d23/12 - d33/16, d11/6 - d33/12, d13/12 - d12/4 - d23/6 + d33/16, d11/3 + d13/4 + d33/12, d12/4 + d13/6 + d23/6 + d33/16, d33/24 - d11/3, d12/4 - d13/6 + d23/12 - d33/16; % - d12/4 - d13/12 - d23/12 - d33/16, - d22/6 - d23/4 - d33/24, d12/4 + d13/12 - d23/6 - d33/16, d33/24 - d22/3, d12/4 + d13/6 + d23/6 + d33/16, d22/3 + d23/4 + d33/12, d23/12 - d13/6 - d12/4 + d33/16, d22/6 - d33/12; % d11/6 - d33/12, d12/4 + d13/12 - d23/6 - d33/16, d13/4 - d11/6 - d33/24, d12/4 - d13/12 - d23/12 + d33/16, d33/24 - d11/3, d23/12 - d13/6 - d12/4 + d33/16, d11/3 - d13/4 + d33/12, d13/6 - d12/4 + d23/6 - d33/16; % d13/12 - d12/4 - d23/6 + d33/16, d33/24 - d22/3, d12/4 - d13/12 - d23/12 + d33/16, d23/4 - d22/6 - d33/24, d12/4 - d13/6 + d23/12 - d33/16, d22/6 - d33/12, d13/6 - d12/4 + d23/6 - d33/16, d22/3 - d23/4 + d33/12]; case 2 % B not / 2, length of cell = 1 % E = [E1111 E1122 E1112 % E2222 E2212 % sym: E1212] % d11 d12 d13 d22 d23 d33 keCoeff = [1/3, 0, 1/2, 0, 0, 1/3; 0, 1/4, 1/3, 0, 1/3, 1/4; -1/3, 0, 0, 0, 0, 1/6; 0, 1/4, -1/3, 0, 1/6, -1/4; -1/6, 0, -1/2, 0, 0, -1/6; 0, -1/4, -1/6, 0, -1/6, -1/4; 1/6, 0, 0, 0, 0, -1/3; 0, -1/4, 1/6, 0, -1/3, 1/4; % column 1 end 0, 1/4, 1/3, 0, 1/3, 1/4; 0, 0, 0, 1/3, 1/2, 1/3; 0, -1/4, -1/3, 0, 1/6, 1/4; 0, 0, 0, 1/6, 0, -1/3; 0, -1/4, -1/6, 0, -1/6, -1/4; 0, 0, 0, -1/6, -1/2, -1/6; 0, 1/4, 1/6, 0, -1/3, -1/4; 0, 0, 0, -1/3, 0, 1/6; % column 2 end -1/3, 0, 0, 0, 0, 1/6; 0, -1/4, -1/3, 0, 1/6, 1/4; 1/3, 0, -1/2, 0, 0, 1/3; 0, -1/4, 1/3, 0, 1/3, -1/4; 1/6, 0, 0, 0, 0, -1/3; 0, 1/4, 1/6, 0, -1/3, -1/4; -1/6, 0, 1/2, 0, 0, -1/6; 0, 1/4, -1/6, 0, -1/6, 1/4; % column 3 end 0, 1/4, -1/3, 0, 1/6, -1/4; 0, 0, 0, 1/6, 0, -1/3; 0, -1/4, 1/3, 0, 1/3, -1/4; 0, 0, 0, 1/3, -1/2, 1/3; 0, -1/4, 1/6, 0, -1/3, 1/4; 0, 0, 0, -1/3, 0, 1/6; 0, 1/4, -1/6, 0, -1/6, 1/4; 0, 0, 0, -1/6, 1/2, -1/6; % column 4 end -1/6, 0, -1/2, 0, 0, -1/6; 0, -1/4, -1/6, 0, -1/6, -1/4; 1/6, 0, 0, 0, 0, -1/3; 0, -1/4, 1/6, 0, -1/3, 1/4; 1/3, 0, 1/2, 0, 0, 1/3; 0, 1/4, 1/3, 0, 1/3, 1/4; -1/3 0, 0, 0, 0, 1/6; 0, 1/4, -1/3, 0, 1/6, -1/4; % column 5 end 0, -1/4, -1/6, 0, -1/6, -1/4; 0, 0, 0, -1/6, -1/2, -1/6; 0, 1/4, 1/6, 0, -1/3, -1/4; 0, 0, 0, -1/3, 0, 1/6; 0, 1/4, 1/3, 0, 1/3, 1/4; 0, 0, 0, 1/3, 1/2, 1/3; 0, -1/4, -1/3, 0, 1/6, 1/4; 0, 0, 0, 1/6, 0, -1/3; % column 6 end 1/6 0, 0, 0, 0, -1/3; 0, 1/4, 1/6, 0, -1/3, -1/4; -1/6 0, 1/2, 0, 0, -1/6; 0, 1/4, -1/6, 0, -1/6, 1/4; -1/3, 0, 0, 0, 0, 1/6; 0, -1/4, -1/3, 0, 1/6, 1/4; 1/3 0, -1/2, 0, 0, 1/3; 0, -1/4, 1/3, 0, 1/3, -1/4; % column 7 end 0, -1/4, 1/6, 0, -1/3, 1/4; 0, 0, 0, -1/3, 0, 1/6; 0, 1/4, -1/6, 0, -1/6, 1/4; 0, 0, 0, -1/6, 1/2, -1/6; 0, 1/4, -1/3, 0, 1/6, -1/4; 0, 0, 0, 1/6, 0, -1/3; 0, -1/4, 1/3, 0, 1/3, -1/4; 0, 0, 0, 1/3, -1/2, 1/3; % column 8 end ]; K = keCoeff*E; sK = K(:); % Ke % [ 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] case 3 % B/2, length of cell = 1/50 % d11 d12 d13 d22 d23 d33 keCoeff = [ 0.25 0 0.2 0 0 0.06251; 0 0.25 0.125 0 0.125 0.0625; -0.25 0 0 0 0 0.06249; 0 0.25 -0.125 0 0.125 -0.0625; -0.25 0 -1/4 0 0 -0.06249; 0 -1/4 -0.125 0 -0.125 -0.0625; 0.25 0 0 0 0 -0.06251; 0 -1/4 0.125 0 -0.125 0.0625; % 8 0 1/4 0.125 0 0.125 0.0625; 0 0 0 0.25 1/4 0.06251; 0 -1/4 -0.125 0 0.125 0.0625; 0 0 0 0.25 0 -0.06251; %0.25*d22 - 0.06251*d33 0 -1/4 -0.125 0 -0.125 -0.0625; 0 0 0 -0.25 -1/4 -0.06249; 0 1/4 0.125 0 -0.125 -0.0625; 0 0 0 -0.25 0 0.06249; % 8 -0.25 0 0 0 0 0.06249; 0 -1/4 -0.125 0 0.125 0.0625; 0.25 0 -1/4 0 0 0.06251; 0 -1/4 0.125 0 0.125 -0.0625; 0.25 0 0 0 0 -0.06251; 0 1/4 0.125 0 -0.125 -0.0625; -0.25 0 1/4 0 0 -0.06249; 0 1/4 -0.125 0 -0.125 0.0625; % 8 0 1/4 -0.125 0 0.125 -0.0625; 0 0 0 0.25 0 -0.06251; % 0.25*d22 - 0.06251*d33 0 -1/4 0.125 0 0.125 -0.0625; 0 0 0 0.25 -1/4 0.06251; 0 -1/4 0.125 0 -0.125 0.0625; 0 0 0 -0.25 0 0.06249; 0 1/4 -0.125 0 -0.125 0.0625; 0 0 0 -0.25 1/4 -0.06249; % 8 -0.25 0 -1/4 0 0 -0.06249; 0 -1/4 -0.125 0 -0.125 -0.0625; 0.25 0 0 0 0 -0.06251; 0 -1/4 0.125 0 -0.125 0.0625; 0.25 0 1/4 0 0 0.06251; 0 1/4 0.125 0 0.125 0.0625; -0.25 0 0 0 0 0.06249; 0 0.25 -0.125 0 0.125 -0.0625; % 8 % 0.25*d12 - 0.125*d13 + 0.125*d23 - 0.0625*d33 0 -1/4 -0.125 0 -0.125 -0.0625; 0 0 0 -0.25 -1/4 -0.06249; 0 1/4 0.125 0 -0.125 -0.0625; 0 0 0 -0.25 0 0.06249; 0 1/4 0.125 0 0.125 0.0625; 0 0 0 0.25 1/4 0.06251; 0 -1/4 -0.125 0 0.125 0.0625; 0 0 0 0.25 0 -0.06251; % 8 0.25 0 0 0 0 -0.06251; 0 1/4 0.125 0 -0.125 -0.0625; -0.25 0 1/4 0 0 -0.06249; 0 1/4 -0.125 0 -0.125 0.0625; -0.25 0 0 0 0 0.06249; 0 -1/4 -0.125 0 0.125 0.0625; 0.25 0 -1/4 0 0 0.06251; 0 -1/4 0.125 0 0.125 -0.0625; % 8 0 -1/4 0.125 0 -0.125 0.0625; 0 0 0 -0.25 0 0.06249; 0 1/4 -0.125 0 -0.125 0.0625; 0 0 0 -0.25 1/4 -0.06249; 0 1/4 -0.125 0 0.12 -0.0625;% 0 0 0 0.25 0 -0.06251; 0 -1/4 0.125 0 0.125 -0.0625; 0 0 0 0.25 -1/4 0.06251;]; K = keCoeff*E; sK = K(:); % d11 = D(1,1); % d12 = D(1,2); % d13 = D(1,3); % d22 = D(2,2); % d23 = D(2,3); % d33 = D(3,3); % Ke = ... % [ 0.25*d11 + 0.25*d13 + 0.06251*d33, 0.25*d12 + 0.125*d13 + 0.125*d23 + 0.0625*d33, 0.06249*d33 - 0.25*d11, 0.25*d12 - 0.125*d13 + 0.125*d23 - 0.0625*d33, - 0.25*d11 - 0.25*d13 - 0.06249*d33, - 0.25*d12 - 0.125*d13 - 0.125*d23 - 0.0625*d33, 0.25*d11 - 0.06251*d33, 0.125*d13 - 0.25*d12 - 0.125*d23 + 0.0625*d33 % 0.25*d12 + 0.125*d13 + 0.125*d23 + 0.0625*d33, 0.25*d22 + 0.25*d23 + 0.06251*d33, 0.125*d23 - 0.125*d13 - 0.25*d12 + 0.0625*d33, 0.25*d22 - 0.06251*d33, - 0.25*d12 - 0.125*d13 - 0.125*d23 - 0.0625*d33, - 0.25*d22 - 0.25*d23 - 0.06249*d33, 0.25*d12 + 0.125*d13 - 0.125*d23 - 0.0625*d33, 0.06249*d33 - 0.25*d22 % 0.06249*d33 - 0.25*d11, 0.125*d23 - 0.125*d13 - 0.25*d12 + 0.0625*d33, 0.25*d11 - 0.25*d13 + 0.06251*d33, 0.125*d13 - 0.25*d12 + 0.125*d23 - 0.0625*d33, 0.25*d11 - 0.06251*d33, 0.25*d12 + 0.125*d13 - 0.125*d23 - 0.0625*d33, 0.25*d13 - 0.25*d11 - 0.06249*d33, 0.25*d12 - 0.125*d13 - 0.125*d23 + 0.0625*d33 % 0.25*d12 - 0.125*d13 + 0.125*d23 - 0.0625*d33, 0.25*d22 - 0.06251*d33, 0.125*d13 - 0.25*d12 + 0.125*d23 - 0.0625*d33, 0.25*d22 - 0.25*d23 + 0.06251*d33, 0.125*d13 - 0.25*d12 - 0.125*d23 + 0.0625*d33, 0.06249*d33 - 0.25*d22, 0.25*d12 - 0.125*d13 - 0.125*d23 + 0.0625*d33, 0.25*d23 - 0.25*d22 - 0.06249*d33 % - 0.25*d11 - 0.25*d13 - 0.06249*d33, - 0.25*d12 - 0.125*d13 - 0.125*d23 - 0.0625*d33, 0.25*d11 - 0.06251*d33, 0.125*d13 - 0.25*d12 - 0.125*d23 + 0.0625*d33, 0.25*d11 + 0.25*d13 + 0.06251*d33, 0.25*d12 + 0.125*d13 + 0.125*d23 + 0.0625*d33, 0.06249*d33 - 0.25*d11, 0.25*d12 - 0.125*d13 + 0.125*d23 - 0.0625*d33 % - 0.25*d12 - 0.125*d13 - 0.125*d23 - 0.0625*d33, - 0.25*d22 - 0.25*d23 - 0.06249*d33, 0.25*d12 + 0.125*d13 - 0.125*d23 - 0.0625*d33, 0.06249*d33 - 0.25*d22, 0.25*d12 + 0.125*d13 + 0.125*d23 + 0.0625*d33, 0.25*d22 + 0.25*d23 + 0.06251*d33, 0.125*d23 - 0.125*d13 - 0.25*d12 + 0.0625*d33, 0.25*d22 - 0.06251*d33 % 0.25*d11 - 0.06251*d33, 0.25*d12 + 0.125*d13 - 0.125*d23 - 0.0625*d33, 0.25*d13 - 0.25*d11 - 0.06249*d33, 0.25*d12 - 0.125*d13 - 0.125*d23 + 0.0625*d33, 0.06249*d33 - 0.25*d11, 0.125*d23 - 0.125*d13 - 0.25*d12 + 0.0625*d33, 0.25*d11 - 0.25*d13 + 0.06251*d33, 0.125*d13 - 0.25*d12 + 0.125*d23 - 0.0625*d33 % 0.125*d13 - 0.25*d12 - 0.125*d23 + 0.0625*d33, 0.06249*d33 - 0.25*d22, 0.25*d12 - 0.125*d13 - 0.125*d23 + 0.0625*d33, 0.25*d23 - 0.25*d22 - 0.06249*d33, 0.25*d12 - 0.125*d13 + 0.125*d23 - 0.0625*d33, 0.25*d22 - 0.06251*d33, 0.125*d13 - 0.25*d12 + 0.125*d23 - 0.0625*d33, 0.25*d22 - 0.25*d23 + 0.06251*d33]; end