Local Stiffness Matrix Negative

4 views (last 30 days)
Ramses Young
Ramses Young on 10 Apr 2022
Answered: Alan Stevens on 11 Apr 2022
Hey Everyone,
I am trying to construct a small 5 element, 2 noded finite element model of a tapered beam using matlab.
I currently have an acceptable answer, but I was wondering if there was any way that the clever people of the forum may have a more elegant solution.
Currently I am creating a global stiffness matrix using a connectivity matrix, describing the nodal connections. Since it is a beam, the nodes just line up with eachother simply.
In my nested loop I currently have the code at the bottom as a way of looping and checking if it is the corner of the sub-matrix 2x2. Do you all know if there is a better way of doing this so i dont have to use a if statement? It seems excessive.
B=[1 2;2 3;3 4;4 5;5 6]; %connectivity matrix
E=10e6 %youngs modulus, psi
L=24 %length of beam, inches
kg=zeros(6,6)
A=[11449/15625; 7921/15625; 5041/15625; 2809/15625; 49/625]; %area of cross section at each node
for e=1:5 %number of elements
for i=1:length(B(1,:)) %number of nodes per element
for j=1:length(B(1,:)) %number of nodes per element
if i==2 & j==2
k=A(e)*E/(L/5);
elseif i==2 & j==1
k=-A(e)*E/(L/5);
elseif i==1 & j==2
k=-A(e)*E/(L/5) ;
else
k=A(e)*E/(L/5);
end
kg(B(e,i),B(e,j))=kg(B(e,i),B(e,j))+k;
end
end
end

Answers (1)

Alan Stevens
Alan Stevens on 11 Apr 2022
How about replacing
if i==2 & j==2
k=A(e)*E/(L/5);
elseif i==2 & j==1
k=-A(e)*E/(L/5);
elseif i==1 & j==2
k=-A(e)*E/(L/5) ;
else
k=A(e)*E/(L/5);
end
by
k=A(e)*E/(L/5)*((i==j) - (i~=j));

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!