By Vat Sun and Chandara Ven
Development of truss Equation with finite element method
1. Local and global coordinates
In this development, we are examining a bar that is positioned at an angle of 𝜃 in relation to the
global 𝑥 axis. The bar's local axis, denoted by 𝑥 ′ , extends from node 1 to node 2 in the direction of
the bar, as illustrated in the Fig. 1.
𝑦 𝑇 𝑥′
2
𝑢2′ , 𝑓′2𝑥
𝐿
𝑦′
1
𝑇 𝑢1′ , 𝑓′1𝑥
𝜃 𝑥
Fig. 1 Tensile force on bar element in local and global coordinates.
Hook’s law
𝜎𝑥 = 𝐸𝜀𝑥 (1)
𝑑𝑢 𝑢2′ − 𝑢1′ (2)
𝜀𝑥 = =
𝑑𝑥 𝐿
The element of stiffness matrix is derived.
𝑢2′ − 𝑢1′ (3)
𝑇 = 𝐴𝜎𝑥 = 𝐴𝐸
𝐿
′
𝑢2′ − 𝑢1′ (4)
𝑓2𝑥 = 𝑇 = 𝐴𝐸
𝐿
In matrix form
𝑓′ 𝐴𝐸 1 −1 𝑢1′ (5)
{ 1𝑥
′ } = [ ]{ }
𝑓1𝑥 𝐿 −1 1 𝑢2′
{𝑓′} = [𝑘′]{𝑑′} (6)
𝐴𝐸 1 −1 (7)
[𝑘′] = [ ]
𝐿 −1 1
1
By Vat Sun and Chandara Ven
2. Coordinate of displacement in rotation axes with angle θ
In numerous situations, it is useful to utilize both local (𝑥 ′ , 𝑦 ′ ) and global (𝑥, 𝑦) coordinates. The local
coordinates are always selected to represent the individual element usefully, while the global
coordinates are chosen to be useful for the entire structure.
𝑦′ 𝑦
𝑣
𝑥′
𝑣′
𝐴⃗ 𝑢′
𝜃 𝑥
𝑢
Fig. 2 Coordinate in rotational axes.
The coordinates of vector 𝐴⃗ can be expressed in both local (𝑥 ′ , 𝑦 ′ ) and global (𝑥, 𝑦) coordinates. The
relationship between these two coordinate systems can be established by a transformation matrix.
𝑢′ = 𝑢 cos 𝜃 + 𝑣 sin 𝜃 (8)
′
𝑣 = −𝑢 sin 𝜃 + 𝑣 cos 𝜃 (9)
In matrix form
cos 𝜃 sin 𝜃 𝑢
{𝑢′} = [ ]{ }
𝑣′ − sin 𝜃 cos 𝜃 𝑣
𝑐 𝑠 𝑢
{𝑢′} = [ ]{ }
𝑣′ −𝑠 𝑐 𝑣
Where 𝑐 = cos 𝜃 , 𝑠 = sin 𝜃
3. Global stiffness matrix in two-dimensional coordinates
From eq. (8):
𝑢1′ = 𝑢1 cos 𝜃 + 𝑣1 sin 𝜃
𝑢2′ = 𝑢2 cos 𝜃 + 𝑣2 sin 𝜃
In matrix form
𝑢1
𝑢′ 𝑐 𝑠 0 0 𝑣1
{ 1′ } = [ ]{ }
𝑢2 0 0 𝑐 𝑠 𝑢2
𝑣2
𝑐 𝑠 0 0
𝑀=[ ]
0 0 𝑐 𝑠
2
By Vat Sun and Chandara Ven
{𝑑′} = [𝑀]{𝑑} (10)
Similarly
𝑓1𝑥
𝑓′ 𝑐 𝑠 0 0 𝑓1𝑦
{ 1𝑥
′ } =[ ]
𝑓2𝑥 0 0 𝑐 𝑠 𝑓2𝑥
{𝑓2𝑦 }
{𝑓′} = [𝑀]{𝑓} (11)
𝑢1 𝑓1𝑥
𝑣 𝑓1𝑦
{𝑑} = {𝑢1 } , {𝑓} =
2 𝑓2𝑥
𝑣2 {𝑓2𝑦 }
From eq. (7)
{𝑓′} = [𝑘′]{𝑑′}
[𝑀]{𝑓} = [𝑘′][𝑀]{𝑑} (12)
By
𝑐 0
𝑐 𝑠 0 0 𝑠 0 2 2
0 ] = [1 0
𝑀 ∙ 𝑀′ = [ ]∙[ ] = [𝑐 + 𝑠 ]=𝐼=1
0 0 𝑐 𝑠 0 𝑐 0 𝑐 + 𝑠2
2 0 1
0 𝑠
From eq. (12) above
[𝑀]{𝑓} = 𝐼[𝑘′][𝑀]{𝑑}
[𝑀]{𝑓} = 𝑀 ∙ 𝑀′ [𝑘′][𝑀]{𝑑}
Then we got
{𝑓} = 𝑴′ [𝒌′][𝑴]{𝑑} (13)
We get the global stiff ness matrix for an element:
𝑐 0
𝒌 = 𝑀′ [𝑘′][𝑀] = [
𝑠 0] ∙ 𝐴𝐸 [ 1 −1] ∙ [ 𝑐 𝑠 0 0
]
0 𝑐 𝐿 −1 1 0 0 𝑐 𝑠
0 𝑠
𝑐 −𝑐 𝑐2 𝑐𝑠 −𝑐 2 −𝑐𝑠
𝐴𝐸 𝑠 −𝑠] ∙ [ 𝑐 𝑠 0 0 𝐴𝐸 𝑐𝑠 𝑠2 −𝑐𝑠 −𝑠 2 ]
𝑘 = 𝑀′ [𝑘′][𝑀] = [−𝑐 𝑐 ]= [
𝐿 0 0 𝑐 𝑠 𝐿 −𝑐 2 −𝑐𝑠 𝑐2 𝑐𝑠
−𝑠 𝑠 −𝑐𝑠 −𝑠 2 𝑐𝑠 𝑠2
then,
𝑐2 𝑐𝑠 −𝑐 2 −𝑐𝑠 (14)
𝐴𝐸 𝑐𝑠 𝑠2 −𝑐𝑠 −𝑠 2 ]
𝑘= [
𝐿 −𝑐 2 −𝑐𝑠 𝑐2 𝑐𝑠
−𝑐𝑠 −𝑠 2 𝑐𝑠 𝑠2
3
By Vat Sun and Chandara Ven
The eq. (14) denotes the explicit stiffness matrix for a bar positioned at any angle within the x-y
coordinate (global coordinate).
The total stiffness matrix obtained from summation of the stiffness matrix from each element in
truss.
𝑁 (15)
∑[𝑘 𝑒 ] = [𝐾]
𝑒=1
The global nodal force matrix
𝑁 (16)
∑{𝑓 𝑒 } = {𝐹}
𝑒=1
For global nodal displacements {𝐷}
{𝐹} = [𝐾]{𝐷} (17)
4. Degree of freedom in two-dimensional coordinates
In one element there are 2 joint nodes, each node has two degrees of freedom which can move on x
and y directions. The degree of each node can determine with formular:
Horizontal degree (x-direction) 𝑑𝑜𝑓 = 2 × 𝑛𝑜𝑑𝑒 − 1 (18)
Vertical degree (y-direction) 𝑑𝑜𝑓 = 2 × 𝑛𝑜𝑑𝑒 (19)
where 𝑛𝑜𝑑𝑒 is the node number.
The displacement happens only at the nodes without support.
5. Calculate stress of element
Stress of element is axial force per cross-sectional area:
′
𝜎 = 𝐸𝜀 = 𝑓2𝑥 /𝐴
From eq. (3)
𝑢2′ − 𝑢1′ 𝐸 𝑢′
𝜎=𝐸 = [−1 1] { 1′ }
𝐿 𝐿 𝑢2
From eq. (10)
𝑢1
𝐸 𝑐 𝑠 0 0 𝑣1
𝜎 = [−1 1] [ ]{ }
𝐿 0 0 𝑐 𝑠 𝑢2
𝑣2
𝑢1 (20)
𝐸 𝑣
𝜎 = [−𝑐 −𝑠 𝑐 𝑠] {𝑢1 }
𝐿 2
𝑣2
4
By Vat Sun and Chandara Ven
6. Calculate the bond reactions
The reaction of supporters is given by substitute {𝑑} into the original finite element equation:
{𝑅} = [𝐾]{𝐷} − {𝐹} (21)
We only need to use those rows of the structural stiffness matrix that correspond to the fixed
supports. At these supports, we are not supplying an external force so {𝐹} = 0. Our equation
becomes:
{𝑅} = [𝐾]{𝐷} (22)
7. Truss example
𝑅4𝑦
𝐹3𝑦 = −25 𝑘𝑁
𝑅4𝑥 4 4
3
3
3m 2
1 1 2
𝑅1𝑥
𝐹2𝑥 = 20 𝑘𝑁
4m
𝑅1𝑦 𝑅2𝑦
Fig. 3 Truss example
In the truss structure, we want to determine some variables unknown such nodal displacements {𝑑},
reaction of support {𝑅}, normal force in bars 𝑁, and normal stress in bar 𝜎.
7.1. Define the variable input
8 6
4 4 5
7 3
3
2
1 4
1 1 2
3
2
Fig. 4 Degrees of freedom
5
By Vat Sun and Chandara Ven
Table 1 Coordinates of nodes
Node x-coordinate (m) y-coordinate (m)
1 0 0
2 4 0
3 4 3
4 0 3
Each element has two nodes.
Table 2 The element and connection nodes
Element From node To Node E (GPa) (mm)
1 1 2 210 200
2 2 3 210 200
3 1 3 210 200
4 3 4 210 200
MATLAB code
% input: node x(m) y(m)
Nodes=[1 0 0
2 4 0
3 4 3
4 0 3];
% node1 node2 E(GPa) A(mm^2)
Elements=[ 1 2 210 200
2 3 210 200
1 3 210 200
3 4 210 200];
NUMND=size(Nodes,1); % Number of nodes
NUMBAR=size(Elements,1); % Number of elements
NUMDEG=2*NUMND; % Number of degrees of freedom
7.2. Apply forces
MATLAB code
Force=zeros(NUMDEG,1); %{F}
Force(3)=20; % kN
Force(6)=-25; % kN
7.3. Support
The degrees of freedom (dof) for constrained displacements are numbered as 1, 2, 4, 7, and 8.
6
By Vat Sun and Chandara Ven
MATLAB code
Support=[1 2 4 7 8];
7.4. Structural stiffness matrix
Element 1
1 2 3 4
10500 0 −10500 0 1
𝑘1 = [ 0 0
0
0 0] 2
0 3
−10500 10500
0 0 0 0 4
or in structural matrix
1 2 3 4 5 6 7 8
10500 0 −10500 0 0 0 0 0 1
0 0 0 0 0 0 0 0 2
−10500 0 10500 0 0 0 0 0 3
𝑘1 = 0 0 0 0 0 0 0 0 4
0 0 0 0 0 0 0 0 5
0 0 0 0 0 0 0 0 6
0 0 0 0 0 0 0 0 7
[ 0 0 0 0 0 0 0 0] 8
Element 2
3 4 5 6
0 0 0 −14000 3
𝑘2 = [0 14000 0 0 ]4
0 0 0 0 5
0 −14000 0 14000 6
Element 3
1 2 5 6
5376 4032 −5376 −4032 1
𝑘3 = [ 4032 3024 −4032 −5376] 2
−5376 −4032 5376 4032 5
−4032 −5376 4032 5376 6
Element 4
5 6 7 8
10500 0 −10500 0 5
𝑘4 = [ 0 0 0 0] 6
−10500 0 10500 0 7
0 0 0 0 8
The result of structural stiffness matrix
4
∑[𝑘 𝑒 ] = [𝐾]
𝑒=1
7
By Vat Sun and Chandara Ven
1 2 3 4 5 6 7 8
15876 4032 −10500 0 −5376 −4032 0 0 1
4032 3024 0 0 −4032 −3024 0 0 2
−10500 0 10500 0 0 0 0 0 3
[𝐾] = 0 0 0 14000 0 −14000 0 0 4
−5376 −4032 0 0 15876 4032 −10500 0 5
−4032 −3024 0 −14000 4032 17024 0 0 6
0 0 0 0 −10500 0 10500 0 7
[ 0 0 0 0 0 0 0 0] 8
MATLAB code
% Calculate two dimensional stiffness matrix [K]
K=zeros(NUMDEG);
for n=1:NUMBAR
E=Elements(n,3);A=Elements(n,4);
NA=Elements(n,1);NB=Elements(n,2);
x1=Nodes(NA,2);
y1=Nodes(NA,3);
x2=Nodes(NB,2);
y2=Nodes(NB,3);
Le=sqrt((x2-x1)^2+(y2-y1)^2);
c=(x2-x1)/Le;
s=(y2-y1)/Le;
Ke=E*A/Le*[c^2 c*s -c^2 -c*s
c*s s^2 -c*s -s^2
-c^2 -c*s c^2 c*s
-c*s -s^2 c*s s^2];
ae=[2*NA-1 2*NA 2*NB-1 2*NB];
K(ae,ae)=K(ae,ae)+Ke;
end
{𝐹} = [𝐾]{𝑑}
In the given equation, the structural or global stiffness matrix is represented by [𝐾], displacement of
each node is denoted by {𝑑}, and the external force matrix is represented by {𝐹}.
𝑑1
15876 4032 −10500 0 −5376 −4032 0 0 0
4032 3024 0 0 −4032 −3024 0 0 𝑑2 0
−10500 0 10500 0 0 0 0 0 𝑑3 20
0 0 0 14000 0 −14000 0 0 𝑑4 = 0
−5376 −4032 0 0 15876 4032 −10500 0 𝑑5 0
−4032 −3024 0 −14000 4032 17024 0 0 𝑑6 −25
0 0 0 0 −10500 0 10500 0 𝑑7 0
[ 0 0 0 0 0 0 0 0] {𝑑 } { 0 }
8
Boundary conditions exist at the fixed supports, where we assume that these joints are immovable
in the restricted direction. Accordingly, we exclude them from our matrix. The degrees of freedom
(dof) for constrained displacements are numbered as 1, 2, 4, 7, and 8. The rows and columns to be
removed in equation above are indicated by the lines.
8
By Vat Sun and Chandara Ven
𝑑1
15876 4032 −10500 0 −5376 −4032 0 0 0
4032 3024 0 0 −4032 −3024 0 0 𝑑2 0
−10500 0 10500 0 0 0 0 0 𝑑3 20
0 0 0 14000 0 −14000 0 0 𝑑4 = 0
−5376 −4032 0 0 15876 4032 −10500 0 𝑑5 0
−4032 −3024 0 −14000 4032 17024 0 0 𝑑6 −25
0 0 0 0 −10500 0 10500 0 𝑑7 0
[ 0 0 0 0 0 0 0 0] {𝑑 } { 0 }
8
Resulting
10500 0 0 𝑑3 20
[ 0 15876 4032 ] {𝑑5 } = { 0 }
0 4032 17024 𝑑6 −25
Then
𝑑3 0.001905
{𝑑5 } = { 0.000397 }
𝑑6 −0.00156
MATLAB code
% Calcule nodal displacement
ad=setdiff(1:NUMDEG,Support);
d=K(ad,ad)\Force(ad);
D=zeros(NUMDEG,1);
D(ad)=d;
7.5. Stress computation
From eq. (20)
𝑢1
𝐸 𝑣
𝜎 = [−𝑐 −𝑠 𝑐 𝑠] {𝑢1 }
𝐿 2
𝑣2
MATLAB code
% Stress computation
Sigma=zeros(NUMBAR,1);
for n=1:NUMBAR
E=Elements(n,3);
NA=Elements(n,1);NB=Elements(n,2);
x1=Nodes(Elements(n,1),2);
y1=Nodes(Elements(n,1),3);
x2=Nodes(Elements(n,2),2);
y2=Nodes(Elements(n,2),3);
Le=sqrt((x2-x1)^2+(y2-y1)^2);
c=(x2-x1)/Le;
s=(y2-y1)/Le;
ae=[2*NA-1 2*NA 2*NB-1 2*NB];
9
By Vat Sun and Chandara Ven
Sigma(n)=E/Le*[-c -s c s]*D(ae);
end
7.6. Calculate normal force in the element
MATLAB code
% Calculates normal forces in bars N
N=zeros(NUMND,1);
for n=1:NUMBAR
A=Elements(n,4);
N(n)=Sigma(n)*A;
end
7.7. Calculate bond reaction
From eq. (22)
{𝑅} = [𝐾]{𝐷}
0
𝑅1 0
15876 4032 −10500 0 −5376 −4032 0 0
R2 0.001905
4032 3024 0 0 −4032 −3024 0 0
𝑅4 = 0 0 0 14000 0 −14000 0 0 0
0.000397
𝑅7 0 0 0 0 −10500 0 10500 0
{𝑅8 } [ 0 0 0 0 0 0 0 0] −0.00156
0
{ 0 }
Then,
𝑅1 −15.8333
R2 3.1250
𝑅4 = 21.8750
𝑅7 −4.1667
{𝑅8 } { 0 }
MATLAB code
% Calculate bond reaction
R=K(Support,:)*D;
10
