User:Eml4500.f08.echo.mott/hw6

=HW Report 6=

Electric Pylon
% Cross sectional area of elements A = [400]; % mm^2

% Elastic Modulus e = [200]*10^6; % kPa

% Density rho = [7.8]*1000; % g/mm^3

% Coordinates of the nodes coord = [1.5 0 0; 4.53 0 0; 1.88 1.18 0; 3.03 1.18 0; 4.17 1.18 0; 2.2 2.181 0; 3.03 2.181 0; 3.85 2.181 0; 2.53 3.181 0; 3.53 3.181 0; 2.34 3.681 0;3.03 3.681 0; 3.71 3.681 0; 2.16 4.181 0; 2.62 4.181 0; 3.44 4.181 0; 3.89 4.181 0; 2.36 4.5 0; 3.7 4.5 0; 1.8 5.172 0; 4.25 5.172 0; 0 6 0; 1.06 6 0; 1.5 6 0; 1.82 6 0; 2.22 6 0; 2.71 6 0; 3.34 6 0; 3.84 6 0; 4.23 6 0; 4.55 6 0; 5 6 0; 6.05 6 0; .8 6.23 0; 5.25 6.23 0; 1.24 6.355 0;  4.82 6.355 0; 1.5 6.43 0; 2 6.43 0; 2.48 6.43 0; 3.03 6.43 0; 3.58 6.43 0; 4.05 6.43 0; 4.55 6.43 0; 1.5 7.43 0; 4.55 7.43 0]*1000; %mm

% Connectivity Array conn = [ 1 3; 1 4; 4 2; 2 5; 3 4; 4 5; 3 6; 3 7; 5 7; 5 8; 6 7; 7 8; 6 9; 7 9; 7 10; 8 10; 9 10; 9 11; 9 12; 10 12; 10 13; 11 12; 12 13; 11 14; 11 15; 12 15; 12 16; 13 16; 13 17; 14 15; 15 16; 16 17; 14 20; 14 18; 15 18; 16 19; 17 19; 17 21; 18 20; 19 21; 20 24; 20 26; 21 29; 21 31; 22 23; 23 24; 24 25; 25 26; 26 27; 27 28; 28 29; 29 30; 30 31; 31 32; 32 33; 22 34; 23 34;  23 36; 24 36; 24 38; 25 38; 25 39; 26 39; 26 40; 27 40; 27 41; 28 41; 28 42; 29 42; 29 43; 30 43; 30 44; 31 44; 31 37; 32 37;  32 35; 33 35; 34 36; 35 37; 36 38; 37 44; 38 39; 39 40; 40 41; 41 42; 42 43; 43 44; 38 45; 39 45; 43 46; 44 46];

% Location Master Matrix lmm = [1 2 5 6; 1 2 7 8; 7 8 3 4; 3 4 9 10; 5 6 7 8; 7 8 9 10; 5 6 11 12; 5 6 13 14; 9 10 13 14; 9 10 15 16; 11 12 13 14; 13 14 15 16; 11 12 17 18; 13 14 17 18; 13 14 19 20; 15 16 19 20; 17 18 19 20; 17 18 21 22; 17 18 23 24; 19 20 23 24; 19 20 25 26;  21 22 23 24;   23 24 25 26; 21 22 27 28; 21 22 29 30; 23 24 29 30; 23 24 31 32; 25 26 31 32; 25 26 33 34; 27 28 29 30; 29 30 31 32; 31 32 33 34; 27  28 39 40; 27 28 35 36; 29 30 35 36; 31 32 37 38; 33 34 37 38; 33 34 41 42; 35 36 39 40; 37 38 41 42; 39 40 47 48; 39 40 51 52; 41 42  57 58; 41 42 61 62; 43 44 45 46; 45 46 47 48; 47 48 49 50; 49 50 51 52; 51 52 53 54; 53 54 55 56;  55 56 57 58; 57 58 59 60; 59 60 61  62; 61 62 63 64; 63 64 65 66; 43 44 67 68; 45 46 67 68; 45 46 71 72; 47 48 71 72; 47 48 75 76;  49 50 75 76; 49 50 77 78; 51 52 77 78;      51 52 79 80; 53 54 79 80; 53 54 81 82; 55 56 81 82; 55 56 83 84; 57 58 83 84; 57 58 85  86; 59 60 85 86; 59 60 87 88; 61 62 87 88; 61  62 73 74; 63 64 73 74; 63 64 69 70; 65 66 69 70; 67 68 71 72; 69 70 73 74; 71 72 75   76; 73 74 87 88; 75 76 77 78; 77 78 79 80; 79 80  81 82; 81 82 83 84; 83 84 85 86; 85 86 87 88; 75 76 89 90; 77 78 89 90; 85 86 91   92; 87 88 91 92];

% Generating the stiffness matrix K

K = zeros(92); M = zeros(92);

for i = 1:91 l = lmm(i, :); c = conn(i, :); k = PlaneTrussElement(e, A, coord(c,:)); m = LumpMassMatrix(rho, A, coord(c,:)); K(l, l) = K(l, l) + k; M(l, l) = M(l, l) + m;

end

% Force matrix

R = zeros(92, 1); R(66,1) = -1000; [d, rf] = NodalSoln(K, R, [1,2,3,4], zeros(4,1)); truss_reax = [];

for i = 1:91 truss_reax = [truss_reax; PlaneTrussResults(e, A, coord(conn(i,:),:), d(lmm(i,:)))]; [Deps(i), sigma(i,1), force(i,1)] = PlaneTrussResults(e, A, coord(conn(i,:),:), d(lmm(i,:))); end d = d*1000;

% Plot the "superman" structure % model with 3-D beam elements

dof = 2;         %  dof per node:  1= axial disp x,  2= disp y

% obtain the coordinates of all nodes

n_node = 46;            % number of nodes n_elem = 91;            % number of elements total_dof = 2 * n_node; % total dof of system

position(:, 1) = [ 1.5; 0; 0]; position(:, 2) = [ 4.53; 0; 0];

position(:, 3) = [ 1.88; 1.18; 0]; position(:, 4) = [ 3.03; 1.18; 0]; position(:, 5) = [ 4.17; 1.18; 0];

position(:, 6) = [ 2.2; 2.181; 0]; position(:, 7) = [ 3.03; 2.181; 0]; position(:, 8) = [ 3.85; 2.181; 0];

position(:, 9) = [ 2.53; 3.181; 0]; position(:, 10) = [ 3.53; 3.181; 0];

position(:, 11) = [ 2.34; 3.681; 0]; position(:, 12) = [ 3.03; 3.681; 0]; position(:, 13) = [ 3.71; 3.681; 0];

position(:, 14) = [ 2.16; 4.181; 0]; position(:, 15) = [ 2.62; 4.181; 0]; position(:, 16) = [ 3.44; 4.181; 0]; position(:, 17) = [ 3.89; 4.181; 0];

position(:, 18) = [ 2.36; 4.5; 0]; position(:, 19) = [ 3.7; 4.5; 0];

position(:, 20) = [ 1.8; 5.172; 0]; position(:, 21) = [ 4.25; 5.172; 0];

position(:, 22) = [ 0; 6; 0]; position(:, 23) = [ 1.06; 6; 0]; position(:, 24) = [ 1.5; 6; 0]; position(:, 25) = [ 1.82; 6; 0]; position(:, 26) = [ 2.22; 6; 0]; position(:, 27) = [ 2.71; 6; 0]; position(:, 28) = [ 3.34; 6; 0]; position(:, 29) = [ 3.84; 6; 0]; position(:, 30) = [ 4.23; 6; 0]; position(:, 31) = [ 4.55; 6; 0]; position(:, 32) = [ 5; 6; 0]; position(:, 33) = [ 6.05; 6; 0];

position(:, 34) = [ .8; 6.23; 0]; position(:, 35) = [ 5.25; 6.23; 0];

position(:, 36) = [ 1.24; 6.355; 0]; position(:, 37) = [ 4.82; 6.355; 0];

position(:, 38) = [ 1.5; 6.43; 0]; position(:, 39) = [ 2; 6.43; 0]; position(:, 40) = [ 2.48; 6.43; 0]; position(:, 41) = [ 3.03; 6.43; 0]; position(:, 42) = [ 3.58; 6.43; 0]; position(:, 43) = [ 4.05; 6.43; 0]; position(:, 44) = [ 4.55; 6.43; 0];

position(:, 45) = [ 1.5; 7.43; 0]; position(:, 46) = [ 4.55; 7.43; 0]; % Deformed Truss Position dposition(:, 1) = [ 1.5; 0; 0]; dposition(:, 2) = [ 4.53; 0; 0];

dposition(:, 3) = [ 1.88 + d(5); 1.18 + d(6); 0]; dposition(:, 4) = [ 3.03 + d(7); 1.18 + d(8); 0]; dposition(:, 5) = [ 4.17 + d(9); 1.18 + d(10); 0];

dposition(:, 6) = [ 2.2 + d(11); 2.181 + d(12); 0]; dposition(:, 7) = [ 3.03 + d(13); 2.181 + d(14); 0]; dposition(:, 8) = [ 3.85 + d(15); 2.181 + d(16); 0];

dposition(:, 9) = [ 2.53 + d(17); 3.181 + d(18); 0]; dposition(:, 10) = [ 3.53 + d(19); 3.181 + d(20); 0];

dposition(:, 11) = [ 2.34 + d(21); 3.681 + d(22); 0]; dposition(:, 12) = [ 3.03 + d(23); 3.681 + d(24); 0]; dposition(:, 13) = [ 3.71 + d(25); 3.681 + d(26); 0];

dposition(:, 14) = [ 2.16 + d(27); 4.181 + d(28); 0]; dposition(:, 15) = [ 2.62 + d(29); 4.181 + d(30); 0]; dposition(:, 16) = [ 3.44 + d(31); 4.181 + d(32); 0]; dposition(:, 17) = [ 3.89 + d(33); 4.181 + d(34); 0];

dposition(:, 18) = [ 2.36 + d(35); 4.5 + d(36); 0]; dposition(:, 19) = [ 3.7 + d(37); 4.5 + d(38); 0];

dposition(:, 20) = [ 1.8 + d(39); 5.172 + d(40); 0]; dposition(:, 21) = [ 4.25 + d(41); 5.172 + d(42); 0];

dposition(:, 22) = [ 0 + d(43); 6 + d(44); 0]; dposition(:, 23) = [ 1.06 + d(45); 6 + d(46); 0]; dposition(:, 24) = [ 1.5 + d(47); 6 + d(48); 0]; dposition(:, 25) = [ 1.82 + d(49); 6 + d(50); 0]; dposition(:, 26) = [ 2.22 + d(51); 6 + d(52); 0]; dposition(:, 27) = [ 2.71 + d(53); 6 + d(54); 0]; dposition(:, 28) = [ 3.34 + d(55); 6 + d(56); 0]; dposition(:, 29) = [ 3.84 + d(57); 6 + d(58); 0]; dposition(:, 30) = [ 4.23 + d(59); 6 + d(60); 0]; dposition(:, 31) = [ 4.55 + d(61); 6 + d(62); 0]; dposition(:, 32) = [ 5 + d(63); 6 + d(64); 0]; dposition(:, 33) = [ 6.05 + d(65); 6 + d(66); 0];

dposition(:, 34) = [ .8 + d(67); 6.23 + d(68); 0]; dposition(:, 35) = [ 5.25 + d(69); 6.23 + d(70); 0];

dposition(:, 36) = [ 1.24 + d(71); 6.355 + d(72); 0]; dposition(:, 37) = [ 4.82 + d(73); 6.355 + d(74); 0];

dposition(:, 38) = [ 1.5 + d(75); 6.43 + d(76); 0]; dposition(:, 39) = [ 2 + d(77); 6.43 + d(78); 0]; dposition(:, 40) = [ 2.48 + d(79); 6.43 + d(80); 0]; dposition(:, 41) = [ 3.03 + d(81); 6.43 + d(82); 0]; dposition(:, 42) = [ 3.58 + d(83); 6.43 + d(84); 0]; dposition(:, 43) = [ 4.05 + d(85); 6.43 + d(86); 0]; dposition(:, 44) = [ 4.55 + d(87); 6.43 + d(88); 0];

dposition(:, 45) = [ 1.5 + d(89); 7.43 + d(90); 0]; dposition(:, 46) = [ 4.55 + d(91); 7.43 + d(92); 0];

% print the node coord. for i = 1 : n_node; x(i) = position(1,i)*8.3; y(i) = position(2,i)*8.3; z(i) = position(3,i)*8.3; end

for i = 1 : n_node; dx(i) = dposition(1,i)*8.3; dy(i) = dposition(2,i)*8.3; dz(i) = dposition(3,i)*8.3; end

node_connect(1, 1) = 1;  % element 1 node_connect(2, 1) = 3;

node_connect(1, 2) = 1;  % element 2 node_connect(2, 2) = 4;

node_connect(1, 3) = 3;  % element 5 node_connect(2, 3) = 4;

node_connect(1, 4) = 4;  % element 3 node_connect(2, 4) = 2;

node_connect(1, 5) = 4;  % element 6 node_connect(2, 5) = 5;

node_connect(1, 6) = 2;  % element 4 node_connect(2, 6) = 5;

node_connect(1, 7) = 3;  % element 8 node_connect(2, 7) = 7;

node_connect(1, 8) = 5;  % element 9 node_connect(2, 8) = 7;

node_connect(1, 9) = 3;  % element 7 node_connect(2, 9) = 6;

node_connect(1, 10) = 5;  % element 10 node_connect(2, 10) = 8;

node_connect(1, 11) = 6;  % element 11 node_connect(2, 11) = 7;

node_connect(1, 12) = 7;  % element 12 node_connect(2, 12) = 8;

node_connect(1, 13) = 6;  % element 13 node_connect(2, 13) = 9;

node_connect(1, 14) = 7;  % element 14 node_connect(2, 14) = 9;

node_connect(1, 15) = 7;  % element 15 node_connect(2, 15) = 10;

node_connect(1, 16) = 8;  % element 16 node_connect(2, 16) = 10;

node_connect(1, 17) = 9;  % element 17 node_connect(2, 17) = 10;

node_connect(1, 18) = 9;  % element 18 node_connect(2, 18) = 11;

node_connect(1, 19) = 9;  % element 19 node_connect(2, 19) = 12;

node_connect(1, 20) = 10;  % element 20 node_connect(2, 20) = 12;

node_connect(1, 21) = 10;  % element 21 node_connect(2, 21) = 13;

node_connect(1, 22) = 11;  % element 22 node_connect(2, 22) = 12;

node_connect(1, 23) = 12;  % element 23 node_connect(2, 23) = 13;

node_connect(1, 24) = 11;  % element 24 node_connect(2, 24) = 14;

node_connect(1, 25) = 11;  % element 25 node_connect(2, 25) = 15;

node_connect(1, 26) = 12;  % element 26 node_connect(2, 26) = 15;

node_connect(1, 27) = 12;  % element 27 node_connect(2, 27) = 16;

node_connect(1, 28) = 13;  % element 28 node_connect(2, 28) = 16;

node_connect(1, 29) = 13;  % element 29 node_connect(2, 29) = 17;

node_connect(1, 30) = 14;  % element 30 node_connect(2, 30) = 15;

node_connect(1, 31) = 15;  % element 31 node_connect(2, 31) = 16;

node_connect(1, 32) = 16;  % element 32 node_connect(2, 32) = 17;

node_connect(1, 33) = 14;  % element 33 node_connect(2, 33) = 20;

node_connect(1, 34) = 14;  % element 34 node_connect(2, 34) = 18;

node_connect(1, 35) = 15;  % element 35 node_connect(2, 35) = 18;

node_connect(1, 36) = 16;  % element 36 node_connect(2, 36) = 19;

node_connect(1, 37) = 17;  % element 37 node_connect(2, 37) = 19;

node_connect(1, 38) = 17;  % element 38 node_connect(2, 38) = 21;

node_connect(1, 39) = 18;  % element 39 node_connect(2, 39) = 20;

node_connect(1, 40) = 19;  % element 40 node_connect(2, 40) = 21;

node_connect(1, 41) = 20;  % element 41 node_connect(2, 41) = 24;

node_connect(1, 42) = 20;  % element 42 node_connect(2, 42) = 26;

node_connect(1, 43) = 21;  % element 43 node_connect(2, 43) = 29;

node_connect(1, 44) = 21;  % element 44 node_connect(2, 44) = 31;

node_connect(1, 45) = 22;  % element 45 node_connect(2, 45) = 23;

node_connect(1, 46) = 23;  % element 46 node_connect(2, 46) = 24;

node_connect(1, 47) = 24;  % element 47 node_connect(2, 47) = 25;

node_connect(1, 48) = 25;  % element 48 node_connect(2, 48) = 26;

node_connect(1, 49) = 26;  % element 49 node_connect(2, 49) = 27;

node_connect(1, 50) = 27;  % element 50 node_connect(2, 50) = 28;

node_connect(1, 51) = 28;  % element 51 node_connect(2, 51) = 29;

node_connect(1, 52) = 29;  % element 52 node_connect(2, 52) = 30;

node_connect(1, 53) = 30;  % element 53 node_connect(2, 53) = 31;

node_connect(1, 54) = 31;  % element 54 node_connect(2, 54) = 32;

node_connect(1, 55) = 32;  % element 55 node_connect(2, 55) = 33;

node_connect(1, 56) = 22;  % element 56 node_connect(2, 56) = 34;

node_connect(1, 57) = 23;  % element 57 node_connect(2, 57) = 34;

node_connect(1, 58) = 23;  % element 58 node_connect(2, 58) = 36;

node_connect(1, 59) = 24;  % element 59 node_connect(2, 59) = 36;

node_connect(1, 60) = 24;  % element 60 node_connect(2, 60) = 38;

node_connect(1, 61) = 25;  % element 61 node_connect(2, 61) = 38;

node_connect(1, 62) = 25;  % element 62 node_connect(2, 62) = 39;

node_connect(1, 63) = 26;  % element 63 node_connect(2, 63) = 39;

node_connect(1, 64) = 26;  % element 64 node_connect(2, 64) = 40;

node_connect(1, 65) = 27;  % element 65 node_connect(2, 65) = 40;

node_connect(1, 66) = 27;  % element 66 node_connect(2, 66) = 41;

node_connect(1, 67) = 28;  % element 67 node_connect(2, 67) = 41;

node_connect(1, 68) = 28;  % element 68 node_connect(2, 68) = 42;

node_connect(1, 69) = 29;  % element 69 node_connect(2, 69) = 42;

node_connect(1, 70) = 29;  % element 70 node_connect(2, 70) = 43;

node_connect(1, 71) = 30;  % element 71 node_connect(2, 71) = 43;

node_connect(1, 72) = 30;  % element 72 node_connect(2, 72) = 44;

node_connect(1, 73) = 31;  % element 73 node_connect(2, 73) = 44;

node_connect(1, 74) = 31;  % element 74 node_connect(2, 74) = 37;

node_connect(1, 75) = 32;  % element 75 node_connect(2, 75) = 37;

node_connect(1, 76) = 32;  % element 76 node_connect(2, 76) = 35;

node_connect(1, 77) = 33;  % element 77 node_connect(2, 77) = 35;

node_connect(1, 78) = 34;  % element 78 node_connect(2, 78) = 36;

node_connect(1, 79) = 35;  % element 79 node_connect(2, 79) = 37;

node_connect(1, 80) = 36;  % element 80 node_connect(2, 80) = 38;

node_connect(1, 81) = 37;  % element 81 node_connect(2, 81) = 44;

node_connect(1, 82) = 38;  % element 82 node_connect(2, 82) = 39;

node_connect(1, 83) = 39;  % element 83 node_connect(2, 83) = 40;

node_connect(1, 84) = 40;  % element 84 node_connect(2, 84) = 41;

node_connect(1, 85) = 41;  % element 85 node_connect(2, 85) = 42;

node_connect(1, 86) = 42;  % element 86 node_connect(2, 86) = 43;

node_connect(1, 87) = 43;  % element 87 node_connect(2, 87) = 44;

node_connect(1, 88) = 38;  % element 88 node_connect(2, 88) = 45;

node_connect(1, 89) = 39;  % element 89 node_connect(2, 89) = 45;

node_connect(1, 90) = 43;  % element 90 node_connect(2, 90) = 46;

node_connect(1, 91) = 44;  % element 91 node_connect(2, 91) = 46;

% connect all beam elements by connectivity array for i = 1 : n_elem; node_1 = node_connect(1,i); node_2 = node_connect(2,i); %dnode_1 = dnode_connect(1,i); %dnode_2 = dnode_connect(2,i); xx = [x(node_1),x(node_2)]; dxx = [dx(node_1),dx(node_2)]; yy = [y(node_1),y(node_2)]; dyy = [dy(node_1),dy(node_2)]; zz = [z(node_1),z(node_2)]; dzz = [dz(node_1),dz(node_2)]; plot3(xx,yy,zz,'r--') hold on     plot3(dxx,dyy,zz,'-') hold on  end

title('Electric Pylon') xlabel('x') ylabel('y') zlabel('z')

view([0 0 1])     % xy plane view % view([0 1 0])    % xz plane view %view([1 0 0])    % yz plane view

% produce postscript file for later printing % print -dps electric_pylon.ps  %print -

Results

[Max_ax_Stress, Loc] = max(Deps)

% Finding the Maximum Tensile stress and where it is located Max_ax_Stress =

4.5256e-008

Loc =

81 % This means that the maximum tensile stress is located in the 81st dof which correlates to element number 66.

[Min_ax_Stress, Loc] = min(Deps)

% Finding the Maximum Compressive Stress and where it is located.

Min_ax_Stress =

-4.3478e-008

Loc =

55

% This means that the maximum compressive stress is located in the 81st dof which correlates to element number 67.



In order to run this coding, the following three codes used in previous homeworks were used, i.e. PlaneTrussElements, PlaneTrussResults, and NodalSln, and they were all taken from the website for the text Fundamental Finite Element Analysis and Applications by M. Asghar Bhatti.

In addition to these codes, the following m.file had to be created in order to assemble the lump mass matrix for the Electric Pylon system:

This code created the following mass matrix for the pylon:

Monday, November 3rd
Initial conditions can be found when $$t = 0$$.

Notation:

1. Prescribe $$u(x, t=0) = \bar{u} (x)$$


 * Where $$\bar{u} (x)$$ is a known function of displacement

2. Also let $$\frac{\partial u}{\partial t} (x, t=0) = \dot{u} (x, t=0) = \bar{v} (x)$$


 * Where $$\bar{v}(x)$$ is a known function of velocity


 * Now proceed with the Finite Element Method to solve Partial Differential Equations.

First, look at the principle of virtual work (continuous) for the dynamics of an elastic bar:


 * $$ {(1)} \qquad { \frac{\partial}{\partial x} ( E A \frac{\partial u}{\partial x}) + f = m  \ddot u } $$

The discrete equation of motion for the elastic bar can be written as


 * $$ {(2)} \qquad {\underline{M} \underline{\ddot d} + \underline{K}\ \underline{d} = \underline{F}} $$

By inspection and comparison of equations (1) & (2) it is clear that


 * $${\frac{\partial}{\partial x} [ E A \frac{\partial u}{\partial x}]} \equiv {-\underline{K}\ \underline{d}}$$


 * $${f} \equiv {\underline{F}}$$


 * $$ {m} \equiv {\underline{M}\ \underline}$$

Equations (1) & (2) are applicable to a multiple degree of freedom system.

In order to derive (2) from (1), look at a single dof system - a spring/mass combination:


 * [[Image:SpringMassSystem.jpg|270px]]

For the system shown above, $$ m \ddot d + k d = F$$

And,


 * $$(3) \qquad \int_{0}^{L} W(x) \Big\{ \frac{\partial}{\partial x} [ E A \frac{\partial u}{\partial x}] +f - m  \ddot u  \Big\} dx = 0 \qquad for \ all \ W(x)$$

Going from equation (1) to (3) is trivial but going from (3) to (1) is not.

(3) Can be re-written as, $$ \int_{0}^{L} W(x)\ G(x) dx = 0 \qquad for \ all \ W(x)$$

Since (3) is valid for all values of W, then let $$W(x) = G(x)$$

Now (3) becomes:

$$ \int_{0}^{L} G^2(x) dx = 0 $$

$$ G^2 \ge 0 $$ and so $$ G(x) = 0 $$

Wednesday, November 5th
Integration by Parts:

$$r^{ }_{ }(x), s(x)$$

$$(rs)' = r's^{ }_{ } + rs'$$

$$r'=\frac{dr}{dx^{ }_{ }}$$, $$s'=\frac{ds}{dx^{ }_{ }}$$

$$\int (rs)' = \int r's + \int rs'$$, $$\int (rs)' = rs$$

so, $$ rs = \int r's + \int rs'$$

rearranging: $$ \int r's = rs - \int rs'$$

Recalling the Continuous PVW method: 1st term: $$r(x)=(EA)\frac{\partial u}{\partial x}$$,  $$ s(x)=w(x)^{ }_{ }$$

Using Integration by parts:

$$s^{ }_{ }=w(x)$$, $$r^{ }_{ }=(EA)\frac{\partial u}{\partial x}$$

$$\begin{align} \int_{x=0}^{x=L}w(x)\frac{\partial}{\partial x} \Bigg[(EA)\frac{\partial u}{\partial x}\Bigg]dx &=\Bigg[w(EA)\frac{\partial u}{\partial x}\Bigg]^{x=L}_{x=0} -\int^{L}_{0}\frac{\partial w}{\partial x}(EA)\frac{\partial u}{\partial x}dx \\ &=w(L)(EA)(L)\frac{\partial u}{\partial x}(L,t)-w(0)(EA)(0)\frac{\partial u}{\partial x}(0,t) -\int^{L}_{0}\frac{\partial w}{\partial x}(EA)\frac{\partial u}{\partial x}dx \\ \end{align} $$

Model Problem:



Set two boundary conditions:

At x=0, we choose w(x) such that w(0)=0 (i.e. kinematically admissible)

We now have the unknown reaction $$N(0,t)=(EA)(0)\frac{\partial u}{\partial x}(0,t)$$

Using continuous PVW, we have $$w(L)$$ and $$F(t)$$:

$$-\int_{0}^{L}\frac{dw}{dx}(EA)\frac{\partial u}{\partial x}dx+\int_{0}^{L}w(x)[f-m\ddot{u}]dx=0$$ for all $$w^{}_{}(x)$$ such that $$w^{}_{}(x)=0$$.

This leads to the following "weak form" for continuous PVW: $$\int_{0}^{L}w(m\ddot{u})dx+\int_{0}^{L}\frac{dw}{dx}(EA)\frac{\partial u}{\partial x}dx=w(L)F(t)+\int_{0}^{L}wfdx$$ for all $$w^{}_{}(x)$$ such that $$w^{}_{}(0)=0$$

Motivation: Discrete PVW applied to Equation on pg. 10-1

$$\mathbf{w}_{6x1}\cdot([\ ]_{6x2} \begin{Bmatrix} d_3 \\  d_4 \end{Bmatrix} -\mathbf{F}_{6x1})=\mathbf{0}_{1x1}$$

$$\mathbf{F}^T= \begin{bmatrix} F_1 & F_2 & F_3 & F_4 & F_5 & F_6 \end{bmatrix}$$

$$F_3$$ and $$F_4$$ are known forces, but $$F_1$$, $$F_2$$, $$F_5$$, and $$F_6$$ are unknown reactions.

Since w can be selected arbitrarily, we select w such that $$w_1=w_2=w_5=w_6=0$$.

To eliminate equations involving unknown reactions, we eliminate rows 1,2,5, and 6 to obtain the following equation:

$$\mathbf{\bar{K}}_{2x2}\mathbf{\bar{d}}_{2x1}=\mathbf{\bar{F}}_{2x1}$$ or $$\mathbf{\bar{K}}_{2x2}\begin{Bmatrix}d_3\\d_4\end{Bmatrix}=\begin{Bmatrix}F_3\\F_4\end{Bmatrix}$$

It should be noted that $$\mathbf{\bar{w}}\cdot(\mathbf{\bar{K}}\mathbf{\bar{d}}-\mathbf{\bar{F}})=0$$ for all $$\mathbf{\bar{w}}$$ where $$\mathbf{\bar{w}}=\begin{Bmatrix}w_3\\w_4\end{Bmatrix}$$

Friday, November 7th




Assume displacement u(x) for $$x_i^{ } \le x \le x_{i+1}$$ (i.e. $$x \epsilon \left [ x_i, x_{i+1} \right ]$$)

Motivation for linear interpolation of the u(x) 2-bar truss is shown in the figure below.



The deformed shape is a straight line, i.e. there was an implicit assumption of linear interpolation of displacement between 2 nodes.

Consider the case where there are only axial displacements (i.e. zero transverse displacement) Express u(x) in terms of di=u(xi) & di+1=u(xi+1) as a linear function in x (i.e., linear interpolation)

$$u(x)=N_i(x)d_i+N_{i+1}(x)d_{i+1}$$

where Ni(x) & Ni+1(x) are linear functions in x

HW

Ni(x)=? (in book)

END HW

$$N_{i+1}(x)=\frac{x-x_i}{x_{i+1}-x_i}$$





Monday, November 10th
 Continuous PVW to Discrete PVW  (continued)

Lagrangian interpolation

Motivation for form of $$N_i(x)$$ and $$N_{i+1}$$:

1. $$N_i(x)$$ and $$N_{i+1}(x)$$ are linear (straight lines), thus any linear combination of $$N_i$$ and $$N_{i+1}$$ is also linear, and in particular, the expression for $$u(x)$$.
 * $$N_i^{}(x)=\alpha_i+\beta_i x$$, where $$\alpha _i$$ and $$\beta_i$$ are numbers
 * $$N_{i+1}^{}=\alpha_{i+1}+\beta_{i+1} x$$, where $$\alpha_{i+1}$$ and $$\beta_{i+1}$$ are numbers

The linear combination of $$N_i$$ and $$N_{i+1}$$, $$ N_i^{} d_i + N_{i+1}d_{i+1} = (\alpha_i + \beta_i x)d_i + (\alpha_{i+1} + \beta_{i+1}x)d_{i+1} = (\alpha_i d_i + \alpha_{i+1}d_{i+1}) + (\beta_i d_i + \beta_{i+1} d_{i+1})x $$ is clearly a linear function of $$x$$

2. Recall the equation for $$u(x)$$ (interpolation of $$u(x)$$):
 * $$u(x_i)=N_i^{}(x_i)d_i+N_{i+1}^{}(x_{i})d_{i+1}$$ where $$N_i^{}(x_i)=1$$ and $$N_{i+1}^{}(x_{i})=0$$, so $$u^{}_{}(x_i)=d_i$$

Wednesday, November 12th
FEM via PVW (continued)

HW: Show that $$u^{ }_{}(x_{i+1})=d_{i+1}$$

$$u(x_{i+1})=N_i^{}(x_{i+1})d_i+N_{i+1}(x_i)d_{i+1}$$

In this case, $$N_i^{}(x_{i+1})=0$$ and $$N_{i+1}^{}(x_{i+1})=1$$

Substituting these values into the equation, we get $$u(x_{i+1}^{})=(0)\cdot d_i+(1)\cdot d_{i+1}$$

Or $$u(x_{i+1}^{})=d_{i+1}$$

From page 31-3: Interpolate for u(x) Apply same interpolation for u(x), i.e.,

$$w(x)=N_i^{ }(x)w_i+N_{i+1}(x)w_{i+1}$$

Element Stiffness Matrix for element i

$$\beta=\int_{x_i}^{x_{i+1}} \left [ N_i^'w_i + N_{i+1}^'w_{i+1} \right ] (EA) \left [ N_i^'d_i+N_{i+1}^'d_{i+1} \right ]\, dx$$

Where $$w'(x)=N_i^'w_i + N_{i+1}^'w_{i+1}$$

and $$u'(x)=N_i^'d_i+N_{i+1}^'d_{i+1}$$

and $$N_i^':=\frac{dN_i(x)}{dx}$$

Likewise for $$N_{i+1}^'$$

Note $$u(x)=\begin{bmatrix} N_i^'(x) & N_{i+1}(x) \end{bmatrix}_{1x2} \begin{Bmatrix} d_i \\ d_{i+1} \end{Bmatrix}_{2x1}$$

$$\frac{du(x)}{dx}=\begin{bmatrix} N_i^'(x) & N_{i+1}^'(x) \end{bmatrix}_{1x2} \begin{Bmatrix} d_i \\ d_{i+1}\end{Bmatrix}_{2x1}$$

Similarly: $$w(x)=\mathbf{N}(x) \begin{Bmatrix} w_i \\ w_{i+1} \end{Bmatrix}$$

$$\frac{dW(x)}{dx}=B(x) \begin{Bmatrix} w_i \\ w_{i+1} \end{Bmatrix}$$

Recall the elemental degrees of freedom:



$$\begin{Bmatrix} d_i \\ d_{i+1} \end{Bmatrix}=\begin{Bmatrix} d_1^{(i)} \\ d_2^{(i)} \end{Bmatrix} = \mathbf{d}^{(i)}$$

$$\begin{Bmatrix} w_i \\ w_{i+1} \end{Bmatrix}=\begin{Bmatrix} w_1^{(i)} \\ w_2^{(i)} \end{Bmatrix} = \mathbf{w}^{(i)}$$

$$\beta=\int_{x_i}^{x_{i+1}}(\mathbf{Bw}^{(i)})_{1x1}(EA)_{1x1}(\mathbf{Bd}^{(i)})_{1x1}\,dx=\mathbf{w}^{(i)}\cdot (\mathbf{k}^{(i)} \mathbf{d}^{(i)})$$

$$\beta=\int_{x_i}^{x_{i+1}}(EA)(\mathbf{Bw}^{(i)}) \cdot (\mathbf{Bd}^{(i)})\, dx$$

$$(\mathbf{Bw}^{(i)})^T (\mathbf{Bd}^{(i)})$$

Where $$(\mathbf{Bw}^{(i)})^T=\mathbf{w}^{(i)T} \mathbf{B}^T=\mathbf{w}^{(i)} \cdot \mathbf{B}^T$$

$$\beta=\mathbf{w}^{(i)}\cdot \int \mathbf{B}^T (EA) \mathbf{B}$$

$$\mathbf{k}_{2x2}^{(i)}=\int_{x_i}^{x_{i+1}} \mathbf{B}^T(x)_{2x1}(EA)_{1x1}\mathbf{B}(x)_{1x2}\,dx$$

HW $$B(x)=\begin{bmatrix} \mathbf{HW} & \frac{1}{L^{(i)}} \end{bmatrix}$$

L(i)=xi+1-xi (length of element i)

HW6: Consider EA=const $$\mathbf{k}^{(i)}=\frac{EA}{L^{(i)}}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}$$

Transfer of variable coordinates from x to $$\tilde x$$

$$\tilde x := x-x_i$$

$$d \tilde x = dx$$

$$\mathbf{k}^{(i)}=\int_{\tilde x=0}^{\tilde x=L^{(i)}}\mathbf{B}^T(\tilde x)(EA)(\tilde x)\mathbf{B}(\tilde x)\,d \tilde x$$

HW6: Find expression for $$\mathbf{k}^{(i)}$$ using the above forumula

HW6



$$A(\tilde{x})=N_1^{(i)}(\tilde{x})A_1+N_2^{(i)}(\tilde{x})A_2$$

$$E(\tilde{x})=N_1^{(i)}(\tilde{x})E_1+N_2^{(i)}(\tilde{x})E_2$$

Find $$\mathbf{k}^{(i)}$$

Friday, November 14th
HW6

Using the figures from 31-4: $$N_1^{(i)}(\tilde{x})=HW6$$

$$N_2^{(i)}(\tilde{x})=\frac{\tilde{x}}{L^{(i)}}=\begin{cases} \mbox{0 at }\tilde{x}=0\\ \mbox{1 at }\tilde{x}=L^{(i)} \end{cases}$$

Shape functions $$N_1^{(i)}, N_2^{(i)}$$ (basis)

HW6: On page 159 of the book, set $$E_1=E_2=E$$

Let $$A(\tilde{x})$$ be linear as on p33-5. Obtain $$\mathbf{k}^{(1)}$$ from the previous problem and compare to the expression given in book.

$$\frac{E}{L^{(i)}}\frac{(A_1+A_2)}{2}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}=\mathbf{k}^{(i)}$$



Next, compare the general $$\mathbf{k}^{(i)}$$ on p33-5 to the stiffness matrix obtained by using $$\frac{(A_1+A_2)}{2}$$ and $$\frac{(E_1+E_2)}{2}$$

Note $$E_1 \ne E_2$$

$$\frac{(E_1+E_2)(A_1+A_2)}{4L^{(i)}}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}=\mathbf{k}_{average}^{(i)}$$

Find $$\mathbf{k}^{(i)}-\mathbf{k}_{average}^{(i)}$$

Remember: Recall the mean value theorem (MVT) & its relation to the centroid:



$$\int_{x=a}^{x=b} f(x)\,dx=f(\bar x)[b-a]$$ for $$\bar x \epsilon [a,b]$$

$$ \mathbf{a} \le \bar x \le \mathbf{b}$$

$$\int_A x \,dA = \bar x \int_A dA = \bar x A$$

$$\int_{x=a}^{x=b} f(x)g(x)\,dx=f(\bar x)g(\bar x)[b-a]$$

However, in general: $$f(\bar x) \ne \frac{1}{b-a}\int_a^b f(x)\,dx$$ Where the right hand side of the equation is the average value of f

Similarly: $$g(\bar x) \ne \frac{1}{b-a}\int_a^b g(x)\,dx$$ Where the right hand side of the equation is the average value of g

Modify the 2-bar truss code to accomodate the general $$\mathbf{k}^{(i)}$$ as on p33-5



Contributing Team Members
Chris Mott 19:51, 21 November 2008 (UTC)

Bradley LaCroix 20:16, 21 November 2008 (UTC)

Brandon Sell 20:25, 21 November 2008 (UTC)

Chris Sell 20:40, 21 November 2008 (UTC)

Jason Bruce

Charles Marshall 21:12, 21 November 2008 (UTC)