1D Finite Element: A load column

Example of structural problem 1D for a load column

$$-\frac{d}{dx}\bigg(E(x)A(x)\frac{du(x)}{dx}\bigg) = 0$

Considering constant coefficients of the model equation with $$a_1=EA, a_0=0, f=0$. For each element we get a constant stiffness matrix:

$$ K_e = \frac{EA}{h} \left ( \begin{tabular}{cc} 1 & -1 \\ -1 & 1 \end{tabular} \right )$

We will build the global system $$[K] u = Q$ , according to the element decomposition and Boundary Conditions

Contents

Geometry

The column is divided in four equal elements of length h.

numDiv=4; %just to force numDiv to be integer
L=18; %total length
h=L/numDiv;
nodes=(0:h:L)'; %colum vector
for i=1:numDiv
       elem(i,:)=[i, i+1];
end
numNod = size(nodes,1);
numElem = size(elem,1);

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).

E = 2e+11; %equivalent to 2·10^11N
A = 250e-4; %expressed in m^2
Ke=(E*A)/h*[1,-1;-1,1]; %local stiff matrix

Assembly

K=zeros(numNod); %initialize the global Stiff Matrix
for i=1: numElem
       rows=[elem(i,1); elem(i,2)];
       colums= rows;
       K(rows,colums)=K(rows,colums)+Ke; %assembly
end

Loads

The Forces and Loads applied to the structure.

P = 11e+4; %load value
Q = zeros(numNod,1);
Q(2:4) = -2*P;
Q(numNod) = -3e+5;

Boundary Conditions

fixedNodes=1; %Fixed Nodes (global num.)
u(fixedNodes)=0; %BC on the fixed nodes
freeNodes = setdiff(1:numNod,fixedNodes); %Complementary of fixed nodes
% modify the linear system (in this case is not necessary because BC=0)
Q=Q-K(:,fixedNodes)*u(fixedNodes);
Qm=Q(freeNodes);
Km=K(freeNodes,freeNodes);
% solve the System
format short e; %just to a better view of the numbers
um=Km\Qm; %displacement of the free nodes.
u(freeNodes)=um;
u

u =

            0  -8.6400e-04  -1.5300e-03  -1.9980e-03  -2.2680e-03

Post Process

Reaction forces, forces and tensions can be computed from the results. According to Hook's Law the force in each element is proportional to its displacement. Tension is defined as Force/Area.

reactForces = K*u' - Q; %notice that reaction forces appears on fixed nodes
% show results
reactForces'

ans =

   9.6000e+05  -2.9104e-11  -1.4552e-10   1.7462e-10   2.9104e-10


for e=1:numElem
    elongation(e)=u(elem(e,2))-u(elem(e,1)); %change in length for element i
    force(e)=E*A*elongation(e)/h; %assume constant length and Area
    stress(e)=force(e)/A;
end
disp('elongation,    Force     , Stress   for each element');
[elongation' force' stress']

elongation,    Force     , Stress   for each element

ans =

  -8.6400e-04  -9.6000e+05  -3.8400e+07
  -6.6600e-04  -7.4000e+05  -2.9600e+07
  -4.6800e-04  -5.2000e+05  -2.0800e+07
  -2.7000e-04  -3.0000e+05  -1.2000e+07

Exercise 1:

Modify the initial example to consider element length of h=2.25 (numDiv=8) and preserve the loads at the same high where they are placed. Compute the new displacements.

Sol = 0, -4.3200e-04, -8.6400e-04, -1.1970e-03, -1.5300e-03, -1.7640e-03, -1.9980e-03, -2.1330e-03, -2.2680e-03

Exercise 2:

Modify the initial example in order to change the material of element 2 for another one with $$E=5\cdot10^9 N/m^2 \quad and \quad A=125cm^2$.

Compute the new displacements

(Hint: consider E and A as vectors and compute Ke inside the for loop).

E=2e+11*ones(numElem,1);
E(2)=5.0e+9;
A=250e-4*ones(numElem,1);
A(2)=125e-4;
%Assembly
K=zeros(numNod); %initialize the global Stiff Matrix
for i=1:numElem
   Ke=(E(i)*A(i))/h*[1,-1;-1,1]; %local stiff matrix
     .......

Sol= 0, -8.6400e-04, -5.4144e-02, -5.4612e-02, -5.4882e-02

(c)Numerical Factory 2019


Published with MATLAB® R2018a