CE 890 Introduction to Matlab
Matlab Project Assembly of Global Stiffness Matrix
Submitted to Professor Dr. Phanikumar S. Mantha Civil & Environmental Engineering
Submitted By Aqeel Ahmed PID 36846644
Introduction In finite element method (& structural analysis approach), a structure is modeled as an assembly of elements or components with various forms of connection between them. A continuous discrete system is modeled with a finite number of elements interconnected at finite number of nodes. The behaviour of individual elements is characterised by the element's stiffness or flexibility relation, which altogether leads to the system's stiffness or flexibility relation. To establish the element's stiffness or flexibility relation, further leading to the global stiffness/flexibility matrix, MATLAB programming can be effectively used. In this project paper, stiffness matrix has been obtained using different approaches for spring elements and then extended to bar and beam elements. A general code has also been included that is capable of reading from any text file the connectivity matrix and compute the global stiffness matrix. Also, the knowledge of cells in Matlab has been included in the codes which necessarily eased the work. All codes have been developed for a defined problem in hand and results compared to solutions for verification.
Assembling the Global Stiffness Matrix for Spring Elements To develop the stiffness matrix, we take an example of two springs connected together and a force P equal to 15 kN is applied to it. The spring constants are k1 = 100 kN and k2 = 200 kN. The layout is as follows:
Figure 1: Spring System for Two Elements
Solution 1.
Approach to Solution
a. Step 1. It involves discretization of the problem. The domain consists of two springs/elements and connected at nodes.
b. Step 2. Elements need to have connectivity as follows:
Element Number
Node i
Node j
1
1
2
2
2
3
c. Step 3. (Element Stiffness Matrix ). To formulate the stiffness matrix for each spring, we have the stiffness’s of each spring. ( k1 = 100 kN, k2 = 200 kN). Calling the function SpringElementStiffness will give us the 2x2 stiffness matrix for each spring. The details are: MATLAB Code function y = SpringElementStiffness(k)
% This Function claculates the element stiffness matrix for springs with spring stiffness as k. It returns 2x2 stiffness matirx
y = [k -k; -k k];
MATLAB Output
>> k1 = SpringElementStiffness(100)
k1 = 100 -100 -100 100 >> k2 = SpringElementStiffness(200)
k2 = 200 -200 -200 200
d. Step 4 (Assembling the Global Stiffness Matrix for the System) . The system has three nodes; therefore the global stiffness matrix will be 3x3 matrix. To obtain the K matrix, first we setup the zero matrix of size 3x3 and then call the Matlab function “SpringAssemble” to obtain the matrix. The details are: MATLAB Code function y = SpringAssemble(K,k,i,j) % This function will assemble the element stiffness matrix k at node i(left node) and j (right hand node) into global stiffness matrix K K(i,i)=K(i,i)+k(1,1); K(i,j)=K(i,j)+k(1,2); K(j,i)=K(j,i)+k(2,1); K(j,j)=K(j,j)+k(2,2); y = K;
MATLAB Output
>> K = zeros(3,3) K=
0
0
0
0
0
0
0
0
0
>> K = SpringAssemble(K,k1,1,2)
K=
100
-100
0
-100
100
0
0
0
0
>> K = SpringAssemble(K,k2,2,3)
K=
100
-100
0
-100
300
-200
0
-200
200
The same approach is tested for a six spring system having different connectivity of nodes. The details are:
Figure 2: Six-element Spring System
Solution Step 1.
1.
The domain consists of six elements and five nodes. The
connectivity will be:
2.
Step 2.
Element Number
Node i
Node j
1
1
3
2
3
4
3
3
5
4
3
5
5
3
4
6
4
2
Each element has 2x2 stiffness matrix and since there are five
nodes, therefore, K (global) size will be 5x5. Each element stiffness matrix will be obtained by plugging in the ‘k’ (spring constant in kN) of respective spring. The out put is as follows:
>> k1= SpringElementStiffness(100)
k1 = 100 -100 -100 100
>> k2= SpringElementStiffness(200)
k2 = 200 -200 -200 200
>> k3= SpringElementStiffness(300)
k3 = 300 -300 -300 300
>> k4= SpringElementStiffness(400)
k4 = 400 -400 -400 400
>> k5= SpringElementStiffness(500)
k5 = 500 -500 -500 500
>> k6 = SpringElementStiffness(600)
k6 = 600 -600 -600 600
>> K = zeros(5,5)
K=
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
>> K = SpringAssemble(K,k1,1,3)
K=
100
0 -100
0
0
0
0
0
0 100
0
0
0
0
0
0
0
0
0
0
0
0
0 -100
0
>> K = SpringAssemble(K,k2,3,4)
K=
100
0
-100
0
0
0
0
0
0
0
0 300
-200
0
0
0 -200
200
0
0
0
0
0
-100
0
>> K = SpringAssemble(K,k3,3,5)
K=
100
0
-100
0
0
0
-100
0
600
0 0
0 0
-200 -300
0
0
-200
200
0
0
0
-300
0
300
>> K = SpringAssemble(K,k4,3,5)
K=
100
0
-100
0
0
0
0
0
-100
0
1000
-200
-700
0
0
-200
200
0
0
0
-700
0
700
0
0
>> K = SpringAssemble(K,k5,3,4)
K=
100
0
-100
0
0
0
0
0
0
0
-100
0
1500
-700
-700
0
0
-700
700
0
0
0
-700
0
700
>> K = SpringAssemble(K,k6,4,2)
K=
100
0
-100
0
0
0
600
0
-600
0
-100
0
1500
-700
-700
0
-600
-700
1300
0
0
-700
0
0 700
Stiffness Matrix for Bar Element Dealing with bar elements involves 2 degree of freedom (dof) per node (similar to springs). The problem in hand is to obtain the global stiffness matrix of 4 bars connected with node connectivity as shown:
Figure 3: Bar Elements with Node Numbering
Solution Approach - 1
The connectivity will be read through a text file and used in the main program to obtain the global stiffness matrix. For the problem EA/L is assumed to be constant. The connectivity is read from the text file (Node1.txt) and can be varied for any number of elements. The code is as follows:
MATLAB Code clc, clear all elcon = load('Node1.txt') % To read the file regarding the connectivity of the elements [row, col] = size(elcon)
% Arranging the data in matrix form
% Creating the Stiffness Matrix of Zeros Stiffness = zeros(row + 1) number of elements
% The size of K(global) is one plus the
%********************************************** % Defining the Element Stiffness matrix % ********************************************* a = [1 -1; -1 1]
% Assuming EA/L is constant
% ********************************************* % Assembly of Stiffness Matrix %**********************************************
for i=1:row m = elcon(i,2); n = elcon(i,3); Stiffness(m,m) = Stiffness(m,n) = Stiffness(n,m) = Stiffness(n,n) = end Stiffness
MATLAB Output
elcon =
1
1
2
2
2
3
3
3
4
4
4
5
row = 4
col = 3
Stiffness =
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
a=
1
-1
-1
1
Stiffness(m,m) Stiffness(m,n) Stiffness(n,m) Stiffness(n,n)
+ + + +
a(1,1); a(1,2); a(2,1); a(2,2);
Stiffness =
1
-1
0
0
0
-1
2
-1
0
0
0
-1
2
-1
0
0
0
-1
2
-1
0
0
0
-1
1
Approach – 2
The same problem has been addressed by writing the code in a very generalized form. This code requires the input of number of elements and length (L) and computes the global stiffness matrix.
MATLAB Code clc,clear all %********************************************** % Input Data %********************************************** numel = 4
% The number of elements
numnodes = numel + 1 neq = numnodes
% Total number of nodes
connection = [1:numel; 2:numel+1]' % Take care of any number of elements % Location of nodes L = 1 x = [0:numel]'/numel*L K = zeros(neq,neq) % The Assembly of the Global Stiffness Matrix for nel = 1:numel n1 = connection(nel,1); n2 = connection(nel,2); x1 = x(n1); x2 = x(n2); ke = [1 -1;-1 1]; % Assembly of element matrix into Global K Matrix K([n1,n2],[n1,n2])=K([n1,n2],[n1,n2])+ke; end
K
MATLAB Output
numel = 4
numnodes = 5
neq = 5
connection = 1
2
2
3
3
4
4
5
L= 1
x= 0 0.2500 0.5000 0.7500 1.0000
K= 1
-1
0
0
0
-1
2
-1
0
0
0
-1
2
-1
0
0
0
-1
2
-1
0
0
0
-1
1
Approach 3
Another approach to obtain the stiffness matrix is using the cell array. The same has been done using following MATLAB Code MATLAB Code clc, clear all a = [1 -1;-1 1] % Input b1 = [1 b2 = [2 b3 = [3 b4 = [4
the connectivity of the nodes of elements 2] 3] 4] 5]
% Assigning the connectivity to cell b = {b1,b2,b3,b4} K = zeros(5,5) for i = 1:4 for m = 1:2 for n = 1:2 K(b{i}(1,m),b{i}(1,n))=K(b{i}(1,m),b{i}(1,n)) + a(m,n) end end end K
Stiffness Matrix for Beams
The methodology can be developed for the beam elements using 2 degree of freedom per node. The element stiffness matrix will become 4x4 and accordingly the global stiffness matrix dimensions will change. Consider a beam discretized into 3 elements (4 nodes per element) as shown below:
1
3
2
4
3
2
1
4
5
7
6
8
Figure 4: Beam dicretized (4 nodes)
The global stiffness matrix will be 8x8. The MATLAB code to assemble it using arbitrary element stiffness matrix (4x4) is as follows: MATLAB Code clc, clear all numel = 3 nnodes = numel+1 dof = {[1 2 3 4],[3 4 5 6],[5 6 7 8]} K = zeros(nnodes*2) k = {rand(4), rand(4), rand(4)} % Assembling the Global Stiffness Matrix for i = 1:numel for m = 1:4 for n = 1:4 K(dof{i}(1,m),dof{i}(1,n))= K(dof{i}(1,m),dof{i}(1,n))+k{i}(m,n) end end end K
Conclusion The global stiffness matrix can be assembled using different techniques as described above. The approach to address the problem has been improved as understanding of the MATLAB functions and code writing improved. These can be easily extended to account for the matrix multiplication to get nodal degree of freedom and nodal forces/reactions as required.
