Assembling the Global Stiffness Matrix for Spring Elements

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

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

> 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

1.Step 1.The domain consists of six elements and five nodes. The connectivity will be:

Element Number / Node i / Node j
1 / 1 / 3
2 / 3 / 4
3 / 3 / 5
4 / 3 / 5
5 / 3 / 4
6 / 4 / 2

2.Step 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

> K = SpringAssemble(K,k1,1,3)

K =

100 0 -100 0 0

0 0 0 0 0

-100 0 100 0 0

0 0 0 0 0

> K = SpringAssemble(K,k2,3,4)

K =

100 0 -100 0 0

0 0 0 0 0

-100 0 300 -200 0

0 0 -200 200 0

0 0 0 0 0

> K = SpringAssemble(K,k3,3,5)

K =

100 0 -100 0 0

0 0 0 0 0

-100 0 600 -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 0 0

-100 0 1000 -200 -700

0 0 -200 200 0

0 0 -700 0 700

> 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 0 -700 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) % The size of K(global) is one plus the number of elements

%**********************************************

% 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,m) + a(1,1);

Stiffness(m,n) = Stiffness(m,n) + a(1,2);

Stiffness(n,m) = Stiffness(n,m) + a(2,1);

Stiffness(n,n) = Stiffness(n,n) + a(2,2);

end

Stiffness

MATLAB Output

elcon =

1 1 2

2 2 3

3 3 4

4 4 5

row =

col =

Stiffness =

0 0 0 0 0

a =

1 -1

-1 1

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 % Total number of nodes

neq = numnodes

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

MATLAB Output

numel =

numnodes =

neq =

connection =

1 2

2 3

3 4

4 5

L =

x =

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 the connectivity of the nodes of elements

b1 = [1 2]

b2 = [2 3]

b3 = [3 4]

b4 = [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

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:

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

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.