exs_bar2¶

Purpose:

Analysis of a plane truss.

Description:

Consider a plane truss consisting of three bars with the properties \(E=200\) GPa, \(A_1=6.0 \times 10^{-4}\) m², \(A_2=3.0 \times 10^{-4}\) m² and \(A_3=10.0 \times 10^{-4}\) m², and loaded by a single force \(P=80\) kN. The corresponding finite element model consists of three elements and eight degrees of freedom.

Example:

The topology is defined by the matrix:

>> Edof=[1   1  2  5  6;
         2   5  6  7  8;
         3   3  4  5  6];

The stiffness matrix K and the load vector f, are defined by:

>> K=zeros(8);
   f=zeros(8,1); f(6)=-80e3;

The element property vectors ep1, ep2 and ep3 are defined by:

>> E=2.0e11;
>> A1=6.0e-4;     A2=3.0e-4;     A3=10.0e-4;
>> ep1=[E A1];    ep2=[E A2];    ep3=[E A3];

And the element coordinate vectors ex1, ex2, ex3, ey1, ey2 and ey3 by:

>> ex1=[0 1.6];    ex2=[1.6 1.6];    ex3=[0 1.6];
>> ey1=[0 0];      ey2=[0 1.2];      ey3=[1.2 0];

The element stiffness matrices Ke1, Ke2 and Ke3 are computed using bar2e:

>> Ke1=bar2e(ex1,ey1,ep1)

Ke1 =

  1.0e+007 *

    7.5000         0   -7.5000         0
         0         0         0         0
   -7.5000         0    7.5000         0
         0         0         0         0

>> Ke2=bar2e(ex2,ey2,ep2)

Ke2 =

  1.0e+007 *

         0         0         0         0
         0    5.0000         0   -5.0000
         0         0         0         0
         0   -5.0000         0    5.0000

>> Ke3=bar2e(ex3,ey3,ep3)

Ke3 =

  1.0e+007 *

    6.4000   -4.8000   -6.4000    4.8000
   -4.8000    3.6000    4.8000   -3.6000
   -6.4000    4.8000    6.4000   -4.8000
    4.8000   -3.6000   -4.8000    3.6000

Based on the topology information, the global stiffness matrix can be generated by assembling the element stiffness matrices:

>> K=assem(Edof(1,:),K,Ke1);
>> K=assem(Edof(2,:),K,Ke2);
>> K=assem(Edof(3,:),K,Ke3)

K =

  1.0e+008 *

  Columns 1 through 7

    0.7500         0         0         0   -0.7500         0         0
         0         0         0         0         0         0         0
         0         0    0.6400   -0.4800   -0.6400    0.4800         0
         0         0   -0.4800    0.3600    0.4800   -0.3600         0
   -0.7500         0   -0.6400    0.4800    1.3900   -0.4800         0
         0         0    0.4800   -0.3600   -0.4800    0.8600         0
         0         0         0         0         0         0         0
         0         0         0         0         0   -0.5000         0

  Column 8

         0
         0
         0
         0
         0
   -0.5000
         0
    0.5000

Considering the prescribed displacements in bc, the system of equations is solved using the function solveq, yielding displacements a and support forces r:

>> bc= [1 0;2 0;3 0;4 0;7 0;8 0];
>> [a,r]=solveq(K,f,bc)

a =

  1.0e-002 *

         0
         0
         0
         0
   -0.0398
   -0.1152
         0
         0

r =

  1.0e+004 *

    2.9845
         0
   -2.9845
    2.2383
    0.0000
    0.0000
         0
    5.7617

The vertical displacement at the point of loading is 1.15 mm. The section forces es1, es2 and es3 are calculated using bar2s from element displacements ed1, ed2 and ed3 obtained using extract:

>> ed1=extract_ed(Edof(1,:),a);
>> es1=bar2s(ex1,ey1,ep1,ed1)

es1 =

  1.0e+004 *

   -2.9845
   -2.9845

>> ed2=extract_ed(Edof(2,:),a);
>> es2=bar2s(ex2,ey2,ep2,ed2)

es2 =

  1.0e+004 *

    5.7617
    5.7617

>> ed3=extract_ed(Edof(3,:),a);
>> es3=bar2s(ex3,ey3,ep3,ed3)

es3 =

  1.0e+004 *

    3.7306
    3.7306

The normal forces are \(N_1=-29.84\) kN, \(N_2=57.62\) kN and \(N_3=37.31\) kN.

A displacement diagram is displayed using the function eldisp2 and normal force diagram using the function secforce2:

>> figure(1)
>> plotpar=[2 1 0];
>> eldraw2(ex1,ey1,plotpar);
>> eldraw2(ex2,ey2,plotpar);
>> eldraw2(ex3,ey3,plotpar);
>> sfac=scalfact2(ex1,ey1,ed1,0.1);
>> plotpar=[1 2 1];
>> eldisp2(ex1,ey1,ed1,plotpar,sfac);
>> eldisp2(ex2,ey2,ed2,plotpar,sfac);
>> eldisp2(ex3,ey3,ed3,plotpar,sfac);
>> axis([-0.4 2.0 -0.4 1.4]);
>> scalgraph2(sfac,[1e-3 0 -0.3]);
>> title('Displacements')

>> figure(2)
>> plotpar=[2 1];
>> sfac=scalfact2(ex1,ey1,N2(:,1),0.1);
>> secforce2(ex1,ey1,N1(:,1),plotpar,sfac);
>> secforce2(ex2,ey2,N2(:,1),plotpar,sfac);
>> secforce2(ex3,ey3,N3(:,1),plotpar,sfac);
>> axis([-0.4 2.0 -0.4 1.4]);
>> scalgraph2(sfac,[5e4 0 -0.3]);
>> title('Normal force')