University of Florida/Eml4507/s13.team3.GuzyR6

CALFEM Verification edit

First begin by constructing the Edof matrix. Column 1 corresponds to the element number; columns 2 and 3 correspond to the global node that corresponds to the first and second local nodes respectively.

 >> Edof=[1 1 2 3 4
          2 1 2 5 6
          3 3 4 7 8
          4 1 2 5 6
          5 3 4 9 10
          6 3 4 5 6
          7 3 4 11 12
          8 1 2 5 6
          9 1 2 9 10
          10 9 10 11 12
          11 13 14 15 16
          12 5 6 3 4
          13 5 6 9 10
          14 17 18 7 8
          15 9 10 5 6  
          16 19 20 5 6
          17 17 18 7 8
          18 9 10 3 4
          19 17 18 7 8
          20 3 4 15 16
          21 15 16 3 4
          22 17 18 7 8
          23 9 10 3 4
          24 17 18 7 8
          25 1 2 11 12];

Construct an empty stiffness matrix K that has side lengths equal to the number of degrees of freedom times the number of global nodes. For this problem, there are 10 global nodes and motion has 3 degrees of freedom.

 >> K=zeros(30,30);

Constructing the force matrix and constructing the displacement matrix Q. Global nodes and known displacements as shown.

The function bar2e generates the element stiffness matrices for a bar with 3 degrees of freedom.

The Vector containing the Young's modulus and cross sectional area are found and and the x and y coordinates for each element are shown below:

 >> ex1=[0 0];
 >> ey1=[0 3];
 >> ez1=[0 0];
 >> ex2=[0 3];
 >> ey2=[0 1];
 >> ez2=[3 1];
 >> ex3=[0 8];
 >> ey3=[6 3;
 >> ez3=[.5 1];
 >> ex4=[6 0];
 >> ey4=[7 3];
 >> ez4=[.9 .4];
 >> ex5=[8 3];
 >> ey5=[4 1];
 >> ez5=[0 4];
 >> ex6=[4 8];
 >> ey6=[3 3;
 >> ez6=[.5 4];
 >> ex7=[4 0];
 >> ey7=[5 3];
 >> ez7=[.3 1];
 >> ex8=[6 3];
 >> ey8=[0 1];
 >> ez8=[5 .5];
 >> ex9=[9 8];
 >> ey9=[9 3];
 >> ez9=[4 .5];
 >> ex10=[0 0];
 >> ey10=[6 3];
 >> ez10=[0 8];
 >> ex11=[6 3];
 >> ey11=[8 1];
 >> ez11=[8 1];
 >> ex12=[9 8];
 >> ey12=[8 9];
 >> ez12=[8 1];
 >> ex13=[0 0];
 >> ey13=[5 3];
 >> ez13=[05 1];
 >> ex14=[9 3];
 >> ey14=[8 1];
 >> ez14=[4 5];
 >> ex15=[0 9];
 >> ey15=[8 3;
 >> ez15=[8 4];
 >> ex16=[6 0];
 >> ey16=[8 3];
 >> ez16=[4 4];
 >> ex17=[0 8];
 >> ey17=[0 1];
 >> ez17=[.5 1];
 >> ex18=[6 8];
 >> ey18=[7 3];
 >> ez18=[0 .3];
 >> ex19=[6 8];
 >> ey19=[8 3];
 >> ez19=[5 8];
 >> ex20=[3 3];
 >> ey20=[5 1];
 >> ez20=[0 0];
 >> ex21=[0 8];
 >> ey21=[0 3;
 >> ez21=[0 0];
 >> ex22=[8 4];
 >> ey22=[0 3];
 >> ez22=[0 0];
 >> ex23=[5 3];
 >> ey23=[8 1];
 >> ez23=[0 4];
 >> ex24=[0 8];
 >> ey24=[8 8];
 >> ez24=[0 1];
 >> ex25=[0 0];
 >> ey25=[0 8]; 
 >> ez25=[8 8];

The commands:

 >> Ke1=bar2e(ex1, ey1, ez1, ep);
 >> Ke2=bar2e(ex2, ey2, ez2 ep);
 >> Ke3=bar2e(ex3, ey3, ez3 ep);
 >> Ke4=bar2e(ex4, ey4, ez4, ep);
 >> Ke5=bar2e(ex5, ey5, ez5, ep);
 >> Ke6=bar2e(ex6, ey6, ez6, ep);  
 >> Ke7=bar2e(ex7, ey7, ez7, ep);
 >> Ke8=bar2e(ex8, ey8, ez8, ep);
 >> Ke9=bar2e(ex9, ey9, ez9, ep);
 >> Ke10=bar2e(ex10, ey10, ez10, ep);
 >> Ke11=bar2e(ex11, ey11, ez11, ep);
 >> Ke12=bar2e(ex12, ey12, ez12, ep);
 >> Ke13=bar2e(ex13, ey13, ez13, ep);
 >> Ke14=bar2e(ex14, ey14, ez14 ep);
 >> Ke15=bar2e(ex15, ey15, ez15 ep);
 >> Ke16=bar2e(ex16, ey16, ez16, ep);
 >> Ke17=bar2e(ex17, ey17, ez17, ep);
 >> Ke18=bar2e(ex18, ey18, ez18, ep);  
 >> Ke19=bar2e(ex19, ey19, ez19, ep);
 >> Ke20=bar2e(ex20, ey20, ez20, ep);
 >> Ke21=bar2e(ex21, ey21, ez21, ep);
 >> Ke22=bar2e(ex22, ey22, ez22, ep);
 >> Ke23=bar2e(ex23, ey23, ez23, ep);
 >> Ke24=bar2e(ex24, ey24, ez24, ep);
 >> Ke25=bar2e(ex25, ey25, ez25, ep);

This yields 25 stiffness matrices. The global stiffness matrix can be assembled with the element stiffness matrices. This results in a stiffness matrix that is in accord with that computed in Matlab in the previous section. The function solveq takes the stiffness matrix, known force matrix, and displacement matrix and returns the complete force and displacement matrices.

 >> [q r]=solveq(K,F,Q);

The displacement and force matrices are in accord with those computed in Matlab in the previous section. The function extract is used to evaluate the displacements of the local nodes for each beam, which is needed to determine the force in each beam.

 >> 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);
 >> ed10=extract(Edof(10,:),q);

This results in the displacements. The function bar2s can be used to determine the force in each beam. Its input parameters are the x and y coordinates for the beam (ex1 and ey1), the physical properties of the beam (Young's modulus and cross sectional area - ep) and the displacement (ed).

 >> es1=bar2s(ex1, ey1,ez1 ep,ed1)
 >> es1=bar2s(ex2, ey2,ez2 ep,ed2)
 >> es1=bar2s(ex3, ey3,ez3 ep,ed3)
 >> es1=bar2s(ex1, ey1,ez4 ep,ed1)
 >> es1=bar2s(ex2, ey2,ez5 ep,ed2)
 >> es1=bar2s(ex3, ey3,ez6 ep,ed3)
 >> es1=bar2s(ex1, ey1,ez7 ep,ed1)
 >> es1=bar2s(ex2, ey2,ez8 ep,ed2)
 >> es1=bar2s(ex3, ey3,ez9 ep,ed3)
 >> es1=bar2s(ex3, ey3,ez10 ep,ed3)

This results in the following spring forces.

   F_bar =
    1.0e+04 *
   Columns 1 through 9
     0.0000    3.5997    3.5997   -3.5997   -3.5997    5.5981   -5.5981    5.5981   -5.5981
    -0.0000   -3.5997   -3.5997    3.5997    3.5997   -5.5981    5.5981   -5.5981    5.5981
   Columns 10 through 18
    -0.0000   -0.0000   -0.9053    0.9053    1.8114   -1.8114    1.8114   -1.8114    3.4748
     0.0000    0.0000    0.9053   -0.9053   -1.8114    1.8114   -1.8114    1.8114   -3.4748
   Columns 19 through 25
     3.4748   -3.4748   -3.4748   -6.7822    6.7822    6.7822   -6.7822
    -3.4748    3.4748    3.4748    6.7822   -6.7822   -6.7822    6.7822

This is in accord with the forces computed in Matlab in the previous section. The axial stresses are given by dividing the axial forces by the area.

stress = 
    1.0e+04 *
   Columns 1 through 9
    -0.0000   -1.1458   -1.1458    1.1458    1.1458   -1.7819    1.7819   -1.7819    1.7819
   Columns 10 through 18
     0.0000    0.0000    0.2882   -0.2882   -0.5766    0.5766   -0.5766    0.5766   -1.1061
   Columns 19 through 25
    -1.1061    1.1061    1.1061    2.1588   -2.1588   -2.1588    2.1588