User:Eml4507.s13.team7.herran

Given
Consider a plane truss consisting of tree bars with the properties E = 200 GPa, A1 = 6.0 · 10-4 m2, A2 = 3.0 · 10-4 m2 and A3 = 10.0 · 10-4 m2, and loaded by a single force P = 80 kN. The corresponding finite element model consists of three elements and eight degrees of freedom.



Find
Analyze the 2-D truss system utilizing the CALFEM toolbox in MATLAB.

Solution
The finite element model (shown below) matching the truss is made up of twelve degrees of freedom and ten elements.



The topology matrix, Edof, is calculated to be

Edof =

[ 1 1 2 5 6;

2 5 6 7 8;

3 3 4 5 6];

K, a global stiffness matrix, and f, a load vector, are defined. The load P is separated into x and y components then inserted in proper location in the load vector f.

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];

The element coordinate vectors ex1, ex2, ex3, ey1, ey2 and ey3 are defined 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 function bar2e from the CALFEM toolbox is used to compute the element matrices, Ke. The global stiffness matrix is then created using the function assem, integrating these element matrices.

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

The bar2e function takes the element node coordinates, young’s modulus, and cross section area and creates a 4x4 stiffness matrix for a 2-D bar element. This is done using the following script.

E=ep(1);

A=ep(2);

b =[ ex(2) -ex(1); ey(2) -ey(1)];

L=sqrt(b'*b);

Kle =E*A/L*[ 1 -1; -1 1];

n=b'/L;

G =[ n zeros(size(n)); zeros(size(n))  n];

Ke=G'*Kle*G;

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 *

0.7500      0        0          0        -0.7500        0          0       0

0        0        0          0           0           0          0       0

0        0      0.6400    -0.4800     -0.6400      0.4800       0       0

0        0     -0.4800     0.3600      0.4800     -0.3600       0       0

-0.7500     0     -0.6400     0.4800      1.3900     -0.4800       0       0

0        0      0.4800    -0.3600     -0.4800      0.8600       0    -0.5000

0        0        0         0            0           0          0       0

0        0        0         0            0        -0.5000       0     0.5000

The assem function uses this script to assemble element matrices Ke and fe into the global stiffness matrix K and global force vector f according to the topology matrix edof. The script takes in the edof, K, Ke, f, and fe and puts out the K and f matrix and vector, respectively.

[ nie,n]=size(edof);

t=edof(:,2:n);

for i = 1:nie

K(t(i,:),t(i,:)) = K(t(i,:),t(i,:))+Ke;

if nargin==5

f(t(i,:))=f(t(i,:))+fe;

end

end

Using the following loop, the element matrices are computed and assembled.

for i=1:10

Ke=bar2e(Ex(i,:),Ey(i,:),ep);

K=assem(Edof(i,:),K,Ke);

end;

The displacements a and the support forces r are computed by solving the system of equations considering the prescribed displacements in bc. This is done using the solveq function.

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 solveq function uses this script to solve the system of equations and puts the results in a global displacement matrix and reaction force matrix.

if nargin==2 ;

d=K\f ;

elseif nargin==3;

[ nd,nd]=size(K);

fdof =[1:nd]';

d=zeros(size(fdof));

Q=zeros(size(fdof));

pdof=bc(:,1);

dp=bc(:,2);

fdof(pdof )=[];

s=K(fdof,fdof)\(f(fdof)-K(fdof,pdof)*dp);

d(pdof)=dp;

d(fdof)=s;

end

Q=K*d-f;

The results of a show that the displacement at the point of loading is 1.15 mm in the y-direction. The section forces es1, es2 and es3 are calculated using bar2s from element displacements ed1, ed2 and ed3, which are obtained using extract. The element displacements are used to evaluate normal forces. The extract function obtains the global displacements from the Edof and a matrices, and the normal forces are evaluated using the bar2s function.

ed1=extract(Edof(1,:),a);

N1=bar2s(ex1,ey1,ep1,ed1)

N1 =

-2.9845e+004

ed2=extract(Edof(2,:),a);

N2=bar2s(ex2,ey2,ep2,ed2)

N2 =

5.7617e+004

ed3=extract(Edof(3,:),a);

N3=bar2s(ex3,ey3,ep3,ed3)

N3 =

3.7306e+004

Therefore, the normal forces are N1 = −29.84 kN, N2 = 57.62 kN and N3 = 37.31 kN.

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.