Singularity problem

[yusf]

Hi i recently developed a code for my analysis of scissors mechanism (imagine a lift ) here is my code

%----------D matrix (includes connectivity information) ------------
%-----------------------------------------------------------------------------------------
D = [1 2 4 3 5 8 10 9 11 14 16 15 19 20 21 22 24 25 27 28 30 31 33 34 37 39 41 43 ;
3 4 6 7 9 10 12 13 15 16 17 18 22 21 23 26 28 27 29 32 34 33 35 36 38 40 42 44 ];
%-----------------------------------------------------------------------------------------

%--------- Transformation matrices--------------------------------------------------------
T11=[-1/sqrt(2) 0 1/sqrt(2);
0 1 0 ;
-1/sqrt(2) 0 -1/sqrt(2)];

T22=[1/sqrt(2) 0 1/sqrt(2);
0 1 0 ;
-1/sqrt(2) 0 1/sqrt(2)];

T33=[0 1 0;
-1 0 0;
0 0 1];

O=zeros(3,3);

T1=[T11 O O O;
O T11 O O;
O O T11 O;
O O O T11];

T2=[T22 O O O;
O T22 O O;
O O T22 O;
O O O T22];

T3=[T33 O O O;
O T33 O O;
O O T33 O;
O O O T33];
%----------------------------------------------------------------------------------------------
numel=28; % numel=number of element
nnode=44; % nnode=number of nodes
ndof=6; % ndof= number of dof of a node

K=zeros(nnode*ndof,nnode*ndof);

%L=zeros(1,numelb);EA=zeros(1,numelb);EIy=zeros(1,numelb);EIz=zeros(1,numelb);GIx=zeros(1,numelb)

for n=1:numel

%--------assigning lenghts,material properties etc..------------------------------
if n>=1 & n<=24

L=750;EA=1.792e8;EIy=1.423e11;EIz=4.601e10;GIx=5.461e10;

else

L=760;EA=9.048e7;EIy=1.484e10;EIz=1.484e10;GIx=1.187e10;

end

%------stiffness matrix(12*12)----------------------------------------------------

k=[EA/L, 0, 0, 0, 0, 0, -EA/L, 0, 0, 0, 0, 0;
0, 12*EIz/L^3, 0, 0, 0, 6*EIz/L^2, 0,-12*EIz/L^3,0, 0, 0, 6*EIz/L^2;
0, 0, 12*EIy/L^3,0,6*EIy/L^2,0, 0, 0,-12*EIy/L^3,0,-6*EIy/L^2,0;
0, 0, 0, GIx/L 0, 0, 0, 0, 0, -GIx/L 0, 0;
0, 0,6*EIy/L^2, 0, 4*EIy/L, 0, 0, 0, 6*EIy/L^2,0, 2*EIy/L, 0;
0, 6*EIz/L^2 0, 0, 0, 4*EIy/L, 0, -6*EIz/L^2 0, 0, 0, 2*EIz/L;
-EA/L 0, 0, 0, 0, 0, EA/L 0, 0, 0, 0, 0;
0, -12*EIz/L^3 0, 0, 0,-6*EIz/L^2 0, 12*EIz/L^3 0, 0, 0, -6*EIz/L^2;
0, 0, -12*EIy/L^3,0,6*EIy/L^2,0, 0, 0, 12*EIy/L^3,0, 6*EIy/L^2,0;
0, 0, 0, -GIx/L 0, 0, 0, 0, 0, GIx/L 0, 0;
0, 0,-6*EIy/L^2, 0, 2*EIy/L, 0, 0, 0, 6*EIy/L^2,0, 4*EIy/L, 0;
0, 6*EIz/L^2 0, 0, 0, 2*EIy/L, 0, -6*EIz/L^2 0, 0, 0, 4*EIz/L];


%----coordinate transformation------------------------------------------------------

if n==1|n==4|n==5|n==8|n==9|n==12|n==13|n==16|n==17|n==20|n==21|n==24

k=T1.'*k*T1;

elseif n==2|n==3|n==6|n==7|n==10|n==11|n==14|n==15|n==18|n==19|n==22|n==23
k=T2.'*k*T2;

else

k=T3.'*k*T3;

end
%--------assembly--------------------------------------------------------------------

a= D(1,n); % a= first node of nth element
b= D(2,n); % b= second node of nth element
for i=1:ndof
for j=1:ndof % ndof= number of dof of a node
K((a-1)*ndof+i,(a-1)*ndof+j)=K((a-1)*ndof+i,(a-1)*ndof+j)+k(i,j);
K((b-1)*ndof+i,(b-1)*ndof+j)=K((b-1)*ndof+i,(b-1)*ndof+j)+k(ndof+i,ndof+j);
K((a-1)*ndof+i,(b-1)*ndof+j)=K((a-1)*ndof+i,(b-1)*ndof+j)+k(i,ndof+j);
K((b-1)*ndof+i,(a-1)*ndof+j)=K((b-1)*ndof+i,(a-1)*ndof+j)+k(ndof+i,j);
end
end
end

%---------force vector-----------------------------------------------------------------1e8

syms dd
F(1,1)=dd;
F(2:264,1)=0;
syms q
F(99,1)=-498198*q;F(105,1)=-498198*q; F(207,1)=-498198*q; F(213,1)=-498198*q;
F(1,1)=0;
%----------------------------imposing boundary conditions------------------------------


K(1,1)=1e8;K(2,2)=1e8;K(3,3)=1e8;K(4,4)=1e8;K(6,6)=1e8;K(7,7)=1e8;K(8,8)=1e8;K(9,9)=1e8;K(10,10)=1e8;K(12,12)=1e8;
K(109,109)=1e8;K(110,110)=1e8;K(111,111)=1e8;K(112,112)=1e8;K(114,114)=1e8;K(115,115)=1e8;K(116,116)=1e8;
K(117,117)=1e8;K(118,118)=1e8;K(120,120)=1e8;

...........

Well my problem is that i couldnt take inverse of Assembly stiffnes matrix K which is K 264x264 matrix.. what would be the reason for that ..

i used space frame elements to model beams in 3D
can any one give me idea about it..how can i solve this..

 

[johnhors]

If you can't invert your stiffness matrix then it is singular, this is caused by insufficient restraints being applied to the model. In the last section of your "code" you appear to be replacing coefficients of the stiffness matrix with the value 1e8, if the intention is to apply a spring to earth at these locations then this value should be added to the existing coefficient value and not replacing it.

 

[yusf]

hi;

does it make any difference adding or completely replacing it from K..it is juct a number i am adding or replacing..i assigned that large value for that i need to constrain it (fix at some corresponding dof)

 

[xerf]

Do you get an error ?

Some remarks:

If you use the penalty approach at the global K level, in order to impose essential boundary conditions, you must also modify accordingly the global force vector.

If you want to enforce DOF_n=A then
replace:

K_nn = G*K_nn
F_n = G*K_nn*A

where G is the penalty number

Also a penalty number of 1e8 is not sufficiently bigconsidering the usual K_ij you listed are of the same order of magnitude.

A good penalty number should be at least 8-10 orders of magnitude higher than K_ij .

I was thinking that solving systems by using the inverse is not the best approach ( at least for the reason it is about 4 times computationally expensive than elimination).

 

[johnhors]

"does it make any difference adding or completely replacing it from K" - with the values you are using it most certainly does, in comparison with your existing coefficient values you are not assigning a large value, as xerf pointed out. Unless you carry out the method outlined by xerf exactly as stated, you MUST ADD spring rate values to your leading diagonal terms, otherwise by replacing terms you have probably got an ill-conditioned nonsense matrix.

Are you really doing an inverse operation? I had assumed that you were actually doing a Gaussian Elimination.

 

[dfish67]

Why not fix those DOFs? IE zero out the rows and columns.

Damon