assemble global stiffness matrix with varying young's modulus

I am trying to assemble stifness matrix of two element beam whose elastic modus is different in each element. the elastic modulus can take any values as [1 1],[1 2],[2 1] and [2 2];

my code is

nel=2; % number of element
nnel=2; % node per element
ndof=2; % number of degree of freedom per elment
 nnode=(nnel-1)*nel+1;	% total number of nodes in system 
sdof=nnode*ndof;	% total system dofs  
v=[1 2];
N=2;
p=permn(v,N);
[l m]=size(p);
kk=zeros(sdof, sdof) ;
mm=zeros(sdof,sdof); 
force=zeros(sdof, 1); 
index=zeros(nnel*ndof,1); 
kk=zeros(sdof, sdof) ;
i=1:4;
v=[1 2];
N=2;
p=permn(v,N);
[l m]=size(p);
kk(:,:,i)=zeros(sdof, sdof,4);
 for l=1:l
for iel=1:nel
    
    edof = nnel*ndof; 
start = (iel-1)*(nnel-1)*ndof; 
for i=1:edof 
index(i)=start+i; 
end 
x=p(l,iel)
k(:,:,1)=ones(4, 4) ;
k(:,:,2)=2*ones(4, 4) ;
edof = length (index);
for i=1:edof 
ii=index(i) ;
for j=1:edof 
jj=index(j) ;
kk(ii,jj,l)=kk(ii,jj)+k(i,j,x) 

end 
end
end
 end
 %the values of kk(::,1) is coming correct but (k:,:,2) matrix is begining
 %with 2 but it should start with 1,similarlly (k:,:,3)is begining
 %with 3 but it should start with 2 and same proble is with (k:,:,4)
 

Answer 1

[Ayush]

Hi,

It seems like you are trying to assemble the stiffness matrix "kk" for a two-element beam system with varying elastic modulus. The modulus for each element can take the values defined by the permutations of [1 2] over "N"=2 elements, resulting in the combinations [1 1],[1 2],[2 1],[2 2].

The nested loop structure is the core part of your code that is responsible for populating the "kk" matrix with the appropriate values based on the permutations "p". The "kk" matrix is a 3-dimensional array, and thus, instead of "kk(ii,jj)", I tried using "kk(ii, jj,l)" this way, all three indices can access an element. This change gave the required output, as mentioned in the comments of the given code. Refer to the below-modified code for a better understanding:

nel=2; % number of element
nnel=2; % node per element
ndof=2; % number of degree of freedom per element
nnode=(nnel-1)*nel+1; % total number of nodes in system 
sdof=nnode*ndof; % total system dofs  
v=[1 2];
N=2;

% Generate permutations 
[V1, V2] = ndgrid(v, v);
p = [V2(:), V1(:)]; % Concatenate to get the permutations matrix

[l, m]=size(p);
i = 1:4;
kk(:,:,i)=zeros(sdof, sdof, 4); % Initialize stiffness matrix for each permutation
index=zeros(nnel*ndof, 1); % Initialize index matrix

for perm_index=1:l
    for iel=1:nel
        edof = nnel*ndof;
        start = (iel-1)*(nnel-1)*ndof;
        for i=1:edof
            index(i)=start+i;
        end
        x = p(perm_index,iel);
        k(:,:,1)=ones(4, 4);
        k(:,:,2)=2*ones(4, 4);
        edof = length (index);
        for i=1:edof
            ii=index(i);
            for j=1:edof
                jj=index(j);
                kk(ii,jj,perm_index)=kk(ii,jj,perm_index)+k(i,j,x); % Assemble the stiffness matrix
            end
        end
    end
end

disp(kk)

(:,:,1) =

     1     1     1     1     0     0
     1     1     1     1     0     0
     1     1     2     2     1     1
     1     1     2     2     1     1
     0     0     1     1     1     1
     0     0     1     1     1     1


(:,:,2) =

     1     1     1     1     0     0
     1     1     1     1     0     0
     1     1     3     3     2     2
     1     1     3     3     2     2
     0     0     2     2     2     2
     0     0     2     2     2     2


(:,:,3) =

     2     2     2     2     0     0
     2     2     2     2     0     0
     2     2     3     3     1     1
     2     2     3     3     1     1
     0     0     1     1     1     1
     0     0     1     1     1     1


(:,:,4) =

     2     2     2     2     0     0
     2     2     2     2     0     0
     2     2     4     4     2     2
     2     2     4     4     2     2
     0     0     2     2     2     2
     0     0     2     2     2     2

I have used "ndgrid" for generating permutations. For more information on "ndgrid" function refer the below documentation:

https://www.mathworks.com/help/matlab/ref/ndgrid.html

I hope this helps!