2D Structural Linear Link : Truss 2D problems

When we consider the 1D bars we study previously, but now living in a 2D (or 3D), we can assemble them in a structure named TRUSS.

The shape of the final structure can be any shape connected in any stable configuration, trusses typically comprise five or more triangular units constructed with straight members whose ends are connected at joints referred to as nodes. Tipically these structures are bridges, roofs or electrical pylons (see the image below).

(image from Wikipedia)

Example of structural problem 2D: Loaded Bridge

Consider the structure of a bridge (see figure) where the bars are made on iron with

(Giga Pascals) and a section area

Compute the displacements when vertical loads are applied (we don't consider the pillars of the bridge in this problem)

Geometry

Define coordinates and connectivity of the nodes. Here the length units heve been taken in

and this affects the units for the Young modulus (

)

% coordinates for the top nodes are computed from the angle (3600*sin(60�))

high = 3600*sin(60*pi/180);

nodes = [

0 0;

3600 0;

7200 0;

10800 0;

1800 high;

5400 high;

9000 high;

];

elem = [

1 2;

2 3;

3 4;

1 5;

5 2;

5 6;

2 6;

6 3;

6 7;

3 7;

7 4;

];

numNod = size(nodes,1);

numElem = size(elem,1);

ndim = size(nodes,2); %nodes are defined in the plane (2 coordinates)

plotElements(nodes,elem,1);

axis equal

Real constants: Materials and section Area

We have to state the materials involved (Young Modulus) and the cross sectional area (A). We assume all bars are equal in this case, same material and same area. Units must be coherent, usually expressed in International System, but here we use

% Young modulus
Y = 200e+3; %(Conversion: 1GP--->1000N/mm^2)
E = Y*ones(1, numElem); %all elements are made with the same material
% Section Area
secArea = 3250; %here in mm^2
A = secArea*ones(1,numElem);  %all elements have the same cross sectional area

Assembly

K = zeros(ndim*numNod); %initialize the global Stiff Matrix
for e = 1: numElem 
       Ke = planeLinkStiffMatrix(nodes,elem,e,E,A);
       rows = [elem(e,1)*ndim-1; elem(e,1)*ndim; elem(e,2)*ndim-1; elem(e,2)*ndim];
       colums = rows; %notice that here there are 4 rows/columns involved: (x,y) for each node	
       K(rows,colums) = K(rows,colums)+Ke; %assembly
end
Kini = K;

Loads

% The Forces and Loads (natural BC) applied to the structure.
Q = zeros(ndim*numNod,1); %initialization
nodf = 1; %node where the force is applied
Q(nodf*ndim-1) = 0;     % x force component
Q(nodf*ndim) = -280000; % y force component
nodf = 2; %node where the force is applied
Q(nodf*ndim-1) = 0;
Q(nodf*ndim) = -210000;
nodf = 3; %node where the force is applied
Q(nodf*ndim-1) = 0;
Q(nodf*ndim) = -280000;
nodf = 4; %node where the force is applied
Q(nodf*ndim-1) = 0;
Q(nodf*ndim) = -360000;

Boundary Conditions (fix displacements)

% fix displacements (essential BC)
fixedNodes = [];
nNod = 1;
fixedNodes = [fixedNodes, nNod*ndim-1]; %(u1_x=0)
fixedNodes = [fixedNodes, nNod*ndim];   %(u1_y=0)
nNod = 4;
%fixedNodes = [fixedNodes, nNod*ndim-1]; %(u_x=0)
fixedNodes = [fixedNodes, nNod*ndim];    %(u_y=0)
freeNodes = setdiff(1:ndim*numNod,fixedNodes); %Complementary of fixed nodes
% modify the linear system
% Q=Q-K(:,fixedNodes)*u(fixedNodes);   %not needed because u(fixedNodes) = 0.

Solve the System

Km = K(freeNodes,freeNodes);
Qm = Q(freeNodes);
format short e; %just to a better view of the numbers 
um=Km\Qm;
u(freeNodes,1)=um; %this way u is a column vector
u(fixedNodes,1)=0;
displacements=[u(1:ndim:end), u(2:ndim:end)]
displacements = 7×2
            0            0
   7.4611e-01  -6.5764e+00
   2.3130e+00  -6.9928e+00
   3.1337e+00            0
   3.0839e+00  -3.5036e+00
   1.5917e+00  -7.2369e+00
  -4.9741e-02  -3.7333e+00

Post Process

% Compute Reaction Forces (only acting on fixed nodes)

rForces=Kini*u-Q; % K has not been modified

ReacForces=[rForces(1:ndim:end), rForces(2:ndim:end)]

ReacForces = 7×2

  -4.0745e-10   5.1333e+05
  -1.1642e-10  -1.1642e-10
   1.1642e-10  -5.2387e-10
            0   6.1667e+05
   1.1642e-10  -3.4925e-10
  -1.5280e-10  -2.3283e-10
   1.1096e-10            0

% Show original structure and displaced one

figure()

esc=8; %escale factor to magnify the visual displacements

plotDeformedTruss(nodes,elem,u,esc);

% find max displacement and the associated nodes

[val, ind] = max(abs(displacements)) %both in x and y coordinates together

val = 1×2

   3.1337e+00   7.2369e+00

ind = 1×2

     4     6

Exercise 1:

Suppose that we double the section area of the elements on the bottom. Compute the maximum displacement in y in this new design with same BC.

Sol= 6.3323e+00

Exercise 2:

Compute the displacements of the point E in the Howe roof structure (see figure) when using the same bars as in the initial bridge example.

Sol: ux = 8.4615e-04m, uy = -2.8943e-03m

Exercise 3:

Build the 3D version of these truss problems. First you have to modify the planeLinkStiffMatrix.m function (below) to a 3D version named spatialLinkStiffMatrix.m with the same input and output name of variables. For that, you have to use in the appropriate lines the following code:

    Le = sqrt(x21*x21+y21*y21+z21*z21); 
    coef = ((E(e)*A(e))/Le^3); 
    M = [x21*x21, x21*y21, x21*z21;  
       y21*x21, y21*y21, y21*z21; 
       z21*x21, z21*y21, z21*z21]; 
    Ke = [ M, -M;
          -M, M]; 
    Ke = coef*Ke;  
    

After that you have to adapt the used bar example on this practice to the 3D case, modifying all the needed code lines in order to solve the example shown on the figure using values L=1000 mm, F=200 N for the 3D bars with material constants E=2000N/mm^2 and section area A=200mm^2.

Solution: Displacements for node 4 are (0, -9.1733e+00, -3.0925e+01)

function Ke = planeLinkStiffMatrix(nod,elem,e,E,A)
%
% Stiffness Matrix for a linear bar element living in 2D
% 
% (c)Numerical Factory
%
x1=nod(elem(e,1),1);
y1=nod(elem(e,1),2);
x2=nod(elem(e,2),1);
y2=nod(elem(e,2),2);
 
x21=x2-x1;
y21=y2-y1;
Le=sqrt(x21*x21+y21*y21);
coef=((E(e)*A(e))/Le^3);
M=[x21*x21, x21*y21;
   x21*y21, y21*y21];    
Ke=[M,-M;
    -M,M];
Ke=coef*Ke;
end