User:Egm4507.s13.team4.colocar/FEA.s13.R2.3.cc

Problem Statement
KS 2008 p.103 sec.2.7 pb.21 write a matlab program to assemble the global stiffness matrix, compute the unknown disp dofs, the reaction forces, the member forces; compare the results with exact hand calculation Statics, Mechanics of Materials).

run CALFEM to verify the results, but now do this problem in a different way (as in commercial FE codes).

first, establish 2 arrays: (1) the global node coordinate array “coord”, and (2) the elem connectivity array “conn”. next, write a matlab code to loop over the number of element to call the function “bar2e” of CALFEM to compute the element stiffness matrices. with this method, you let your matlab code construct the arrays ex1, ey1, ex2, ey2, etc. for you, using the info from the arrays “coord” and “conn”.

MATLAB Code
The following is a MATLAB code that assembles the global stiffness matrix, computes the unknown displacement dofs, reaction forces, and member forces function FEA_R2p3 E = 10000; % Young's modulus A = pi*2^2/4; % Area % Lengths of each element L1 = 12; L2 = 12; L3 = 12; L4 = sqrt((12^2)+(9^2)); L5 = 9; L6 = sqrt((12^2)+(9^2)); L7 = 12; L8 = 9; L9 = sqrt((12^2)+(9^2)); % topology matrix % [element# node1x node1y node2x node2y] Edof = [1 1 2  3  4;2  3  4 5 6;3  5  6 7 8; 4 1 2  9 10;5  9 10 3 4;6  9 10 5 6;            7 9 10 11 12;8 11 12 5 6;9 11 12 7 8];    % element stiffness matrix in local coordinates for elements 1 through 9 esmloc1 = (E*A/L1)*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0]; esmloc2 = (E*A/L2)*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0]; esmloc3 = (E*A/L3)*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0]; esmloc4 = (E*A/L4)*[(L1/L4)^2 (L1/L4)*(L5/L4) -(L1/L4)^2 -(L1/L4)*(L5/L4); (L1/L4)*(L5/L4) (L5/L4)^2 -(L1/L4)*(L5/L4) -(L5/L4)^2; -(L1/L4)^2 -(L1/L4)*(L5/L4) (L1/L4)^2 (L1/L4)*(L5/L4); -(L1/L4)*(L5/L4) -(L5/L4)^2 (L1/L4)*(L5/L4) (L5/L4)^2]; % [cos^2 cossin -cos^2 -cossin;cossin sin^2 -cossin -sin^2; % -cos^2 -cossin cos^2 cossin;-cossin -sin^2 cossin sin^2] esmloc5 = (E*A/L5)*[0 0 0 0;0 1 0 -1;0 0 0 0;0 -1 0 1]; esmloc6 = (E*A/L6)*[(L1/L4)^2 -(L1/L4)*(L5/L4) -(L1/L4)^2 (L1/L4)*(L5/L4); -(L1/L4)*(L5/L4) (L5/L4)^2 (L1/L4)*(L5/L4) -(L5/L4)^2; -(L1/L4)^2 (L1/L4)*(L5/L4) (L1/L4)^2 -(L1/L4)*(L5/L4); (L1/L4)*(L5/L4) -(L5/L4)^2 -(L1/L4)*(L5/L4) (L5/L4)^2]; % [cos^2 -cossin -cos^2 cossin;-cossin sin^2 cossin -sin^2; % -cos^2 cossin cos^2 -cossin;cossin -sin^2 -cossin sin^2] esmloc7 = (E*A/L7)*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0]; esmloc8 = (E*A/L8)*[0 0 0 0;0 1 0 -1;0 0 0 0;0 -1 0 1]; esmloc9 = (E*A/L9)*[(L1/L4)^2 -(L1/L4)*(L5/L4) -(L1/L4)^2 (L1/L4)*(L5/L4); -(L1/L4)*(L5/L4) (L5/L4)^2 (L1/L4)*(L5/L4) -(L5/L4)^2; -(L1/L4)^2 (L1/L4)*(L5/L4) (L1/L4)^2 -(L1/L4)*(L5/L4); (L1/L4)*(L5/L4) -(L5/L4)^2 -(L1/L4)*(L5/L4) (L5/L4)^2]; % [cos^2 -cossin -cos^2 cossin;-cossin sin^2 cossin -sin^2; % -cos^2 cossin cos^2 -cossin;cossin -sin^2 -cossin sin^2] % element nodal coordinates % ex1 and ey1 are the x and y coordinates to the nodes of element 1 ex1 = [0 12]; ey1 = [0 0]; ex2 = [12 24]; ey2 = [0 0]; ex3 = [24 36]; ey3 = [0 0]; ex4 = [0 12]; ey4 = [0 9]; ex5 = [12 12]; ey5 = [9 0]; ex6 = [12 24]; ey6 = [9 0]; ex7 = [12 24]; ey7 = [9 9]; ex8 = [24 24]; ey8 = [9 0]; ex9 = [24 36]; ey9 = [9 0]; K = zeros(12); F = zeros(12,1); F(6)=-1200; F(11)=400; ep = [E A]; % assembly of the global stiffness matrix K = assem(Edof(1,:),K,esmloc1); K = assem(Edof(2,:),K,esmloc2); K = assem(Edof(3,:),K,esmloc3); K = assem(Edof(4,:),K,esmloc4); K = assem(Edof(5,:),K,esmloc5); K = assem(Edof(6,:),K,esmloc6); K = assem(Edof(7,:),K,esmloc7); K = assem(Edof(8,:),K,esmloc8); K = assem(Edof(9,:),K,esmloc9); % boundary conditions bc = [1 0;2 0;8 0]; % reaction force vector Q = solveq(K,F,bc); % finding element displacements ed1 = extract(Edof(1,:),Q); ed2 = extract(Edof(2,:),Q); ed3 = extract(Edof(3,:),Q); ed4 = extract(Edof(4,:),Q); ed5 = extract(Edof(5,:),Q); ed6 = extract(Edof(6,:),Q); ed7 = extract(Edof(7,:),Q); ed8 = extract(Edof(8,:),Q); ed9 = extract(Edof(9,:),Q); disp = [ed1;ed2;ed3;ed4;ed5;ed6;ed7;ed8;ed9]; % reaction forces at each element; member forces (stress) N1 = bar2s(ex1,ey1,ep,ed1); MF1 = N1/A; N2 = bar2s(ex2,ey2,ep,ed2); MF2 = N2/A; N3 = bar2s(ex3,ey3,ep,ed3); MF3 = N3/A; N4 = bar2s(ex4,ey4,ep,ed4); MF4 = N4/A; N5 = bar2s(ex5,ey5,ep,ed5); MF5 = N5/A; N6 = bar2s(ex6,ey6,ep,ed6); MF6 = N6/A; N7 = bar2s(ex7,ey7,ep,ed7); MF7 = N7/A; N8 = bar2s(ex8,ey8,ep,ed8); MF8 = N8/A; N9 = bar2s(ex9,ey9,ep,ed9); MF9 = N9/A; RF = [N1;N2;N3;N4;N5;N6;N7;N8;N9]; MF = [MF1;MF2;MF3;MF4;MF5;MF6;MF7;MF8;MF9];

With this, the output is as follows: GlobalStiffnessMatrix = K    dispDofs = disp reactionForces = RF   memberForces = MF

GlobalStiffnessMatrix = 1.0e+003 * 3.9584   1.0053   -2.6180         0         0         0         0         0   -1.3404   -1.0053         0         0    1.0053    0.7540         0         0         0         0         0         0   -1.0053   -0.7540         0         0   -2.6180         0    5.2360         0   -2.6180         0         0         0         0         0         0         0         0         0         0    3.4907         0         0         0         0         0   -3.4907         0         0         0         0   -2.6180         0    6.5764   -1.0053   -2.6180         0   -1.3404    1.0053         0         0         0         0         0         0   -1.0053    4.2446         0         0    1.0053   -0.7540         0   -3.4907         0         0         0         0   -2.6180         0    3.9584   -1.0053         0         0   -1.3404    1.0053         0         0         0         0         0         0   -1.0053    0.7540         0         0    1.0053   -0.7540   -1.3404   -1.0053         0         0   -1.3404    1.0053         0         0    5.2988         0   -2.6180         0   -1.0053   -0.7540         0   -3.4907    1.0053   -0.7540         0         0         0    4.9986         0         0         0         0         0         0         0         0   -1.3404    1.0053   -2.6180         0    3.9584   -1.0053         0         0         0         0         0   -3.4907    1.0053   -0.7540         0         0   -1.0053    4.2446

dispDofs = 0        0    0.3056   -1.4992    0.3056   -1.4992    0.6112   -2.1836    0.6112   -2.1836    1.0695         0         0         0    0.8260   -1.4992    0.8260   -1.4992    0.3056   -1.4992    0.8260   -1.4992    0.6112   -2.1836    0.8260   -1.4992    0.5204   -1.9258    0.5204   -1.9258    0.6112   -2.1836    0.5204   -1.9258    1.0695         0

reactionForces = 1.0e+003 * 0.8000   0.8000    1.2000   -0.5000         0    0.5000   -0.8000    0.9000   -1.5000

memberForces = 254.6479 254.6479  381.9719 -159.1549         0  159.1549  254.6479  286.4789 -477.4648

It can be seen here that the problem done by hand corresponds to the force output of the MATLAB program.

CALFEM MATLAB
% CALFEM % global coordinate matrix Coord = [0 0;12 0;24 0;36 0;12 9;24 9]; % element connectivity array conn = Edof; % global nodal dof matrix Dof = [1 2;3 4;5 6;7 8;9 10;11 12]; % contructing element arrays [ex,ey]=coordxtr(conn,Coord,Dof,2); K=zeros(12); EA=[E A]; % Young's modulus, Area % compute element stiffness matrix Ke and global stiffness matrix for i=1:9 Ke = bar2e(ex(i,:),ey(i,:),EA); % element stiffness matrix K = assem(conn(i,:),K,Ke);     % global stiffness matrix end % element displacement vector ed=extract(conn,Q); % reaction forces for i=1:9 N(i)=bar2s(ex(i,:),ey(i,:),ep,ed(i,:)); end % member forces MF1C = N(1)/A; % stress in AC    MF2C = N(2)/A;  % stress in CE    MF3C = N(3)/A;  % stress in EF    MF4C = N(4)/A;  % stress in AB    MF5C = 0;  % stress in BC    MF6C = N(6)/A;  % stress in BE    MF7C = N(7)/A;  % stress in BD    MF8C = N(8)/A;  % stress in DE    MF9C = N(9)/A;  % stress in DF    MFC = [MF1C;MF2C;MF3C;MF4C;MF5C;MF6C;MF7C;MF8C;MF9C];

The output to the CALFEM code is as follows: CGlobalStiffnessMatrix = K    CdispDofs = ed    CreactionForces = N'    CmemberForces = MFC

CGlobalStiffnessMatrix = 1.0e+003 * 3.9584   1.0053   -2.6180         0         0         0         0         0   -1.3404   -1.0053         0         0    1.0053    0.7540         0         0         0         0         0         0   -1.0053   -0.7540         0         0   -2.6180         0    5.2360         0   -2.6180         0         0         0         0         0         0         0         0         0         0    3.4907         0         0         0         0         0   -3.4907         0         0         0         0   -2.6180         0    6.5764   -1.0053   -2.6180         0   -1.3404    1.0053         0         0         0         0         0         0   -1.0053    4.2446         0         0    1.0053   -0.7540         0   -3.4907         0         0         0         0   -2.6180         0    3.9584   -1.0053         0         0   -1.3404    1.0053         0         0         0         0         0         0   -1.0053    0.7540         0         0    1.0053   -0.7540   -1.3404   -1.0053         0         0   -1.3404    1.0053         0         0    5.2988         0   -2.6180         0   -1.0053   -0.7540         0   -3.4907    1.0053   -0.7540         0         0         0    4.9986         0         0         0         0         0         0         0         0   -1.3404    1.0053   -2.6180         0    3.9584   -1.0053         0         0         0         0         0   -3.4907    1.0053   -0.7540         0         0   -1.0053    4.2446

CdispDofs = 0        0    0.3056   -1.4992    0.3056   -1.4992    0.6112   -2.1836    0.6112   -2.1836    1.0695         0         0         0    0.8260   -1.4992    0.8260   -1.4992    0.3056   -1.4992    0.8260   -1.4992    0.6112   -2.1836    0.8260   -1.4992    0.5204   -1.9258    0.5204   -1.9258    0.6112   -2.1836    0.5204   -1.9258    1.0695         0

CreactionForces = 1.0e+003 * 0.8000   0.8000    1.2000   -0.5000         0    0.5000   -0.8000    0.9000   -1.5000

CmemberForces = 254.6479 254.6479  381.9719 -159.1549         0  159.1549  254.6479  286.4789 -477.4648

The results of all methods verifies their accuracy.