2D Structural Linear Link Problem: Truss 2D problem

Example of structural problem 2D: Loaded Bridge

Consider the structure of a bridge (see figure) where the bars are made on iron with E=200GPa and a section area A=3250mm^2.

Compute the displacements when vertical loads are applied.

Geometry

Define coordinates and connectivity of the nodes.

nodes = [

0 0;

3600 0;

7200 0;

10800 0;

1800 3118;

5400 3118;

9000 3118;

];

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)

plotLinkNodElem(nodes,elem);

Real constants: Materials and section Area

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

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

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;	
       K(rows,colums)=K(rows,colums)+Ke; %assembly
end

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 (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, this is only valid if BC = 0.
Km=K(freeNodes,freeNodes);
Qm=Q(freeNodes);

Solve the System

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.4604e-01  -6.5759e+00
   2.3127e+00  -6.9923e+00
   3.1334e+00            0
   3.0836e+00  -3.5033e+00
   1.5916e+00  -7.2363e+00
  -4.9736e-02  -3.7330e+00

Post Process

% Reaction Forces

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

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

ReacForces = 7×2

   4.1681e-10   5.1333e+05
   1.2509e-12  -2.9104e-10
  -2.1745e-10   1.1642e-10
  -7.3590e-11   6.1667e+05
  -2.2008e-10   1.2808e-10
   1.5142e-10  -1.1864e-10
  -3.7664e-11  -1.2645e-10

% Show original structure and displaced one

figure()

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

plotLinkNodElemDespl(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.1334e+00   7.2363e+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 this new design with same BC.

Sol= 6.3319e+00 mm

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

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