How to generate an nxn matrix using n number of 2x2 matrices?

Hi,

I am having trouble generating a code in the exisiting loop to create a global stiffness matrix from each element stiffness matrix. I am trying to program matlab to do something similar to this.

I have written the code and loop to generate the connectivity matrix and each element stiffness matrix. I am hoping to use the same loop to generate the global stiffness matrix.

e = %connectivity table

     1     1     2
     2     2     3
     3     3     4
     4     4     5
     5     5     6
     6     6     7

%element stiffness matricies
ke(:,:,1) =

    1.5708   -1.5708
   -1.5708    1.5708

ke(:,:,2) =

    0.8836   -0.8836
   -0.8836    0.8836

ke(:,:,3) =

    0.6981   -0.6981
   -0.6981    0.6981

ke(:,:,4) =

    0.6136   -0.6136
   -0.6136    0.6136

ke(:,:,5) =

    0.5655   -0.5655
   -0.5655    0.5655

ke(:,:,6) =

    0.5345   -0.5345
   -0.5345    0.5345
>> 

Can anyone help me?

Thank you!

%Element connectivity Table
n=6; %creates 6 elements 
e=[]; %creates empty connectivity matrix 
k=[1 -1;-1 1]; %standard format of K
ke=[]; %creates empty element stiffness matrix
Ae=[]; %creates empty element area amatrix
ks=zeros(n,n);
E=1; 
L=2;
R1=1;
R2=0.5;
i=1;
for m=1:n; %goes through each element
    i=i;
    j=i+1;
    x=L/m;%thickness of each slice of the beam
    A1=pi()/L^2*(R2*L+(R1-R2)*x)^2;%area of each slice of the beam 
    Ae(:,m)=[A1];%stores each area in a matrix
    e(m,:)= [m i j];%stores results in matrix 
    ke(:,:,m)=E*Ae(:,m)/L*k; %stores each element stiffness matrix
end

[Jan]

What exactly is your question? You ask for creating a nxn matrix using n 2x2 matrices. In the picture there is an example to create a 3c3 matrix from 2 2x2 matrices.

" I am hoping to use the same loop to generate the global stiffness matrix." - If you show us the code of this loop, we might be able to suggest a mdoficiation.

[Jan]

Please use the tools to format your code. See the bar on top of the field for writing the message. Instead of posting the code twice, it would be useful to mention, what exactly your question is.

[Guillaume]

I'm also not sure what the question is. Furthermore, as far as I can tell you explain how to form a matrix for two springs in series but not when they're in parallel. I don't know anything about stiffness matrices, but I would assume the result is not the same for parallel springs.

I can see several possible question with what you have described. Given a connectivity table, the first task would probably be to find which springs are in parallel and which are in series. Once that is done you would then have to assemble the matrices according to the relevant formula. It's not clear which part of this you need help with.

[Monique Embury]

In my problem, the springs are in series.

The image posted is an example of what I am trying to do. In the example, 2 2x2 matrices are combined to make a 3x3 matrix.

I am trying to combine 6 2x2 matrices to create a 6x6 matrix in a similar fashion. I am not sure how to program this. That is what I am seeking help on.

Answer 1

[Guillaume]

"I am trying to combine 6 2x2 matrices to create a 6x6 matrix in a similar fashion"

According to your algorithm the stiffness matrix will be 7x7 not 6x6. It is trivial to achieve

% inputs:
% ke: a 2x2xN matrix

stiffness_matrix = zeros(size(ke, 3) + 1);  %create a (N+1)x(N+1) matrix of 0s
for idx = 1:size(ke, 3)
    stiffness_matrix(idx:idx+1, idx:idx+1) = stiffness_matrix(idx:idx+1, idx:idx+1) + ke(:, :, idx);
end

[Monique Embury]

Thank you very much!

[Monique Embury]

I now have another problem silimar to the other. My element matrices are now 4x4 with 11 of them. I would like to add all the element matricies to make a global one. The problem this time, the elements and nodes are not in numerical order. Can you help with the KG matrix?

E=1.9E6; % lb/in^2 modulus of elasticity
A2=21.5; % in^2 area of cross section
EA2=E*A2;
node=[1 0 0; 2 8.7*12 0; 3 6*12 2.2*12; 4 12*12 4.4*12; 5 15.3*12 0; 6 18*12 2.2*12; 7 24*12 0]%node, x, y 
e2=[1 1 3; 2 1 2; 3 2 3; 4 3 4; 5 2 4; 6 2 5; 7 4 5; 8 4 6; 9 5 6; 10 5 7; 11 6 7] %element, node1, node2
d=11; %11 elements 
KE=zeros(4,4); %empty matrix to store element stiffness maxtrix
for c=1:d;
    I=e2(c,2);
    J=e2(c,3);
    xI=node(I,2);
    yI=node(I,3);
    xJ=node(J,2);
    yJ=node(J,3);
    L=sqrt((xJ-xI)^2+(yJ-yI)^2);
    l=(xJ-xI)/L;
    m=(yJ-yI)/L;
    KE(:,:,c)=EA2/L*[l^2 l*m -l^2 -l*m; l*m m^2 -l*m -m^2; -l^2 -l*m l^2 l*m;  -l*m -m^2 l*m m^2]
    for n=1:2; %number of nodes in each element 
        for o=1:2;
            I=e2(c,n);
            J=e2(c,o);
            for p=1:2; %DOF for each element 
                for t=1:2;
                    KG((I-1)*2+p, (J-1)*2+t)=KG((I-1)*2+p, (J-1)*2+t)+KE((n-1)*2+p, (o-1)*2+t);
                end 
            end
        end 
    end 
end

Thank you!

[Guillaume]

Note: If your original query has been answered, you ought to accept the answer.

"The problem this time, the elements and nodes are not in numerical order."

Not only that, but you now have springs in parallel. In fact, they form a mesh. So, you need to explain how they should be combined. If you plot your connectivity table:

e2=[1 1 3; 2 1 2; 3 2 3; 4 3 4; 5 2 4; 6 2 5; 7 4 5; 8 4 6; 9 5 6; 10 5 7; 11 6 7] %element, node1, node2
g = graph(e2(:, 2), e2(:, 3));
plot(g)

You get

[Monique Embury]

Apologies, I usually am not one to post on forms and don't really know the courtesies or expectations.

That plot would be correct. This is the given image. We are now working with truss elements instead of springs, but the same methods can be applied.

[Guillaume]

" but the same methods can be applied."

As I said earlier, I know nothing about stiffness matrix but I assume that you can no longer sum the matrices along the diagonal. So, how does one construct a stiffness matrix in the general case of a mesh?

[Monique Embury]

They can be added along the diagonal again. Using the element number and corresponding nodes (node 1=i and node 2=j). The element matrix should be added to column and rows: 2*i - 1, 2*i, 2*j - 1, 2*j.

[Guillaume]

Ah ok, then it's easy.

stiffness_matrix = zeros(2*max([e2(:, 2); e2(:, 3)]));
for idx = 1 : size(e2, 1)
    indices = 2 * e2(idx, 2:3) - [1;0];  %requires R2016b or later creates a 2x2 matrix of indices [2i-1, 2j-1; 2i, 2j]
    stiffness_matrix(indices, indices) = stiffness_matrix(indices, indices) + KE(:, :, idx);
end

[Monique Embury]

With this code:

E=1.9E6; % lb/in^2 modulus of elasticity
A2=21.5; % in^2 area of cross section
EA2=E*A2;
node=[1 0 0; 2 8.7*12 0; 3 6*12 2.2*12; 4 12*12 4.4*12; 5 15.3*12 0; 6 18*12 2.2*12; 7 24*12 0]%node, x, y 
e2=[1 1 3; 2 1 2; 3 2 3; 4 3 4; 5 2 4; 6 2 5; 7 4 5; 8 4 6; 9 5 6; 10 5 7; 11 6 7] %element, node1, node2
d=11; %11 elements 
KE=zeros(4,4); %empty matrix to store element stiffness maxtrix
KG=[];
for c=1:d;
    I=e2(c,2);
    J=e2(c,3);
    xI=node(I,2);
    yI=node(I,3);
    xJ=node(J,2);
    yJ=node(J,3);
    L=sqrt((xJ-xI)^2+(yJ-yI)^2);
    l=(xJ-xI)/L;
    m=(yJ-yI)/L;
    KE(:,:,c)=EA2/L*[l^2 l*m -l^2 -l*m; l*m m^2 -l*m -m^2; -l^2 -l*m l^2 l*m;  -l*m -m^2 l*m m^2]
    stiffness_matrix = zeros(2*max([e2(:, 2); e2(:, 3)]));
        for idx = 1 : size(e2, 1)
            indices = 2 * e2(idx, 2:3) - [1;0];  %creates a 2x2 matrix of indices [2i-1, 2j-1; 2i, 2j]
            stiffness_matrix(indices, indices) = stiffness_matrix(indices, indices) + KE(:, :, idx);
        end
end

I get this error:

Index in position 3 exceeds array bounds (must not exceed 1).

Error in MECH532_Project1_Embury (line 98)
            stiffness_matrix(indices, indices) = stiffness_matrix(indices, indices) + KE(:, :, idx);

Any Ideas on how to fix it?

Thank you!

[Guillaume]

The construction of the stiffness matrix can only happen after KE has been fully constructed, not during the construction, so move that code after the c loop.

[Monique Embury]

You are amazing and thank you so very much for your help!

Answer 2

[jahanzaib ahmad]

u just want to add stiffness matrix of each element to get n x n stiffness matrix

423601-assemble-global-stiffness-matrix

[jahanzaib ahmad]

else share code and explain bit more about problem please