I understand that you are trying to find element stiffness matrix for beam element using MATLAB and you find difficulty in assembling the global stiffness matrix.
Below are some of the equations and MATLAB codes for finding out the Element stiffness matrix “K” and Global Stiffness matrix “KG”. Hope they’ll be of some help to reach you to your desired solution.
The element stiffness matrix ‘K’ for beam element is defined as:
where,
E = interval modulus of elasticity,
Ia = moment of inertia,
l = length.
It is of order 4x4 due to the four degrees of freedom, two at each nodes.
The following is the MATLAB code of a function which is used to find out Beam element stiffness matrix “K”.
function y = BeamElementStiffness(a, b, alpha, Ia, l) %% Here, a and b are the left and the right bounds of interval modulus of elasticity and alpha belongs to [0, 1].
E=a+(b-a)*alpha;
y=((E*Ia)/l^3)*[12 6*l -12 6*l;6*l 4*l^2 -6*l 2*l^2; -12 -6*l 12 -6*l;6*l 2*l^2 -6*l 4*l^2];
end
The following MATLAB function is used to assemble the elemental stiffness matrix K of the beam element whose nodes are i (left-end node) and j (right-end node) into the global stiffness matrix "KG". It gives 2nx2n global stiffness matrix "KG".
function y = BeamElementAssemble (KG, K, i, j)
KG(2*i-1,2*i-1)=KG(2*i-1,2*i-1)+K(1,1);
KG(2*i-1,2*i)=KG(2*i-1,2*i)+K(1,2);
KG(2*i-1,2*j-1)=KG(2*i-1,2*j-1)+K(1,3);
KG(2*i-1,2*j)=KG(2*i-1,2*j)+K(1,4);
KG(2*i,2*i-1)=KG(2*i,2*i-1)+K(2,1);
KG(2*i,2*i)=KG(2*i,2*i)+K(2,2);
KG(2*i,2*j-1)=KG(2*i,2*j-1)+K(2,3);
KG(2*i,2*j)=KG(2*i,2*j)+K(2,4);
KG(2*j-1,2*i-1)=KG(2*j-1,2*i-1)+K(3,1);
KG(2*j-1,2*i)=KG(2*j-1,2*i)+K(3,2);
KG(2*j-1,2*j-1)=KG(2*j-1,2*j-1)+K(3,3);
KG(2*j-1,2*j)=KG(2*j-1,2*j)+K(3,4);
KG(2*j,2*i-1)=KG(2*j,2*i-1)+K(4,1);
KG(2*j,2*i)=KG(2*j,2*i)+K(4,2);
KG(2*j,2*j-1)=KG(2*j,2*j-1)+K(4,3);
KG(2*j,2*j)=KG(2*j,2*j)+K(4,4);
y=KG;
end
