why my stiffness matrix is singular even after applying correct boundary conditions

19 views (last 30 days)
abhishek tiwari
abhishek tiwari on 28 Jun 2014
Answered: John D'Errico on 29 Jun 2014
Hi all, I have written a program for solving a 20-noded brick element by fEM. I have checked the determinants of the jacobians; all are positive and non zero. I have checked almost everything I could but still when I apply BC my stiffness matrix is singular. Please please please help. I am sending the codes with coordinates, connectivity and integration points.
xx = [-55
-55
-55
-55
-35
-35
-35
-35
-35
-35
-35
-35
-55
-55
-55
-55
-45
-45
-45
-45
];
yy = [10
25
10
25
10
25
10
25
10
17.5
25
17.5
17.5
25
17.5
10
25
10
25
10
];
zz = [5
5
0
0
5
5
0
0
2.5
0
2.5
5
5
2.5
0
2.5
5
5
0
0
];
% total_nodes = size(xx,2);
% % MCO = [-1 -1 -1; 1 -1 -1; 1 1 -1;-1 1 -1;-1 -1 1;1 -1 1; 1 1 1; -1 1 1;...
% 0 -1 -1;1 0 -1; 0 1 -1;-1 0 -1;-1 -1 0; 1 -1 0; 1 1 0; -1 1 0;...
% 0 -1 1; 1 0 1; 0 1 1; -1 0 1];
conn = [1 5 6 8 7 1 2 4 3 12 11 10 9 13 14 15 16 18 17 19 20];
%mapped coordinates of elemen nodes for C3D20 elements MCO = maped coord
% MCO = [-1 -1 -1; 1 -1 -1; 1 1 -1;-1 1 -1;-1 -1 1;1 -1 1; 1 1 1; -1 1 1;...
% 0 -1 -1;1 0 -1; 0 1 -1;-1 0 -1;-1 -1 0; 1 -1 0; 1 1 0; -1 1 0;...
% 0 -1 1; 1 0 1; 0 1 1; -1 0 1];
syms eps tau eta;
total_dof = 3*(size(xx,1));
phi = sparse(60,60);
Kg = sparse(total_dof,total_dof);
K_local = sparse(60,60);
total_no_element = 1;
Dl = sparse(3,20);
Dm = sparse(3,20);
X = sparse(60,2);
P = sparse(60,1);
b = zeros(6,3);
B = sparse(6,60);
x = sparse(3,20);
g = sparse(20,1);
h = sparse(20,1);
y = sparse(20,1);
z = sparse(20,1);
dof = zeros(3,1,20,total_no_element);
%-----------Material property matrix--------------------------------------
prompt = 'What is the Elastic modulus? ';
result = input(prompt);
E = result;
prompt = 'What is the poissons ratio? ';
result = input(prompt);
v = result;
f = 1-v;
g1 = (1-2*v)/2;
M = (E/((1+v)*(1-2*v)))*[f v v 0 0 0;v f v 0 0 0; v v f 0 0 0;...
0 0 0 g1 0 0;0 0 0 0 g1 0; 0 0 0 0 0 g1];
N(1).n = (eps/8 - 1/8)*(eta - 1)*(tau - 1)*(eps + eta + tau + 2);
N(9).n = 0.25*(1-eps^2)*(1-eta)*(1-tau);
N(2).n = (eps/8 + 1/8)*(eta - 1)*(tau - 1)*(eps - eta - tau - 2);
N(10).n = 0.25*(1+eps)*(1-eta^2)*(1-tau);
N(3).n = (eps/8 + 1/8)*(eta + 1)*( 1-tau )*(eps + eta - tau - 2);
N(11).n = 0.25*(1-eps^2)*(1+eta)*(1-tau);
N(4).n = -(eps/8 - 1/8)*(eta + 1)*(-tau + 1)*(-eps + eta - tau - 2);
N(12).n = 0.25*(1-eps)*(1-eta^2)*(1-tau);
N(17).n = 0.25*(1-eps)*(1-eta)*(1-tau^2);
N(18).n = 0.25*(1+eps)*(1-eta)*(1-tau^2);
N(19).n = 0.25*(1+eps)*(1+eta)*(1-tau^2);
N(20).n = 0.25*(1-eps)*(1+eta)*(1-tau^2);
N(5).n = (eps/8 - 1/8)*(eta - 1)*(tau + 1)*(-eps - eta + tau - 2);
N(13).n = 0.25*(1-eps^2)*(1-eta)*(1+tau);
N(14).n = 0.25*(1+eps)*(1-eta^2)*(1+tau);
N(15).n = 0.25*(1-eps^2)*(1+eta)*(1+tau);
N(16).n = 0.25*(1-eps)*(1-eta^2)*(1+tau);
N(6).n = (1/8)*(1+eps)*(1-eta)*(1+tau)*(-2+eps-eta+tau);
N(7).n = (1/8)*(1+eps)*(1+eta)*(1+tau)*(-2+eps+eta+tau);
N(8).n = (1/8)*(1-eps)*(1+eta)*(1+tau)*(-2-eps+eta+tau);
for i = 1:20;
J(i).eps = diff(N(i).n,eps);
J(i).eta = diff(N(i).n,eta);
J(i).tau = diff(N(i).n,tau);
end
Dl = [J(1).eps J(2).eps J(3).eps J(4).eps J(5).eps J(6).eps J(7).eps ...
J(8).eps J(9).eps J(10).eps J(11).eps J(12).eps J(13).eps J(14).eps...
J(15).eps J(16).eps J(17).eps J(18).eps J(19).eps J(20).eps;J(1).eta...
J(2).eta J(3).eta J(4).eta J(5).eta J(6).eta J(7).eta ...
J(8).eta J(9).eta J(10).eta J(11).eta J(12).eta J(13).eta J(14).eta...
J(15).eta J(16).eta J(17).eta J(18).eta J(19).eta J(20).eta;J(1).tau...
J(2).tau J(3).tau J(4).tau J(5).tau J(6).tau J(7).tau ...
J(8).tau J(9).tau J(10).tau J(11).tau J(12).tau J(13).tau J(14).tau...
J(15).tau J(16).tau J(17).tau J(18).tau J(19).tau J(20).tau];
coord = [double(-0.774596669241483) 0 double(0.774596669241483)];
weights = [5/9 8/9 5/9];
IP = zeros(27,3);
IP(1,:) = [coord(1) coord(1) coord(1)];
%6 points
IP(2,:) = [coord(2) coord(1) coord(1)];
IP(3,:) = [coord(3) coord(1) coord(1)];
IP(4,:) = [coord(1) coord(2) coord(1)];
IP(5,:) = [coord(2) coord(2) coord(1)];
IP(6,:) = [coord(3) coord(2) coord(1)];
IP(7,:) = [coord(1) coord(3) coord(1)];
%12 points
IP(8,:) = [coord(2) coord(3) coord(1)];
IP(9,:) = [coord(3) coord(3) coord(1)];
IP(10,:) = [coord(1) coord(1) coord(2)];
IP(11,:) = [coord(2) coord(1) coord(2)];
IP(12,:) = [coord(3) coord(1) coord(2)];
IP(13,:) = [coord(1) coord(2) coord(2)];
IP(14,:) = [coord(2) coord(2) coord(2)];
IP(15,:) = [coord(3) coord(2) coord(2)];
IP(16,:) = [coord(1) coord(3) coord(2)];
IP(17,:) = [coord(2) coord(3) coord(2)];
IP(18,:) = [coord(3) coord(3) coord(2)];
IP(19,:) = [coord(1) coord(1) coord(3)];
%8 points
IP(20,:) = [coord(2) coord(1) coord(3)];
IP(21,:) = [coord(3) coord(1) coord(3)];
IP(22,:) = [coord(1) coord(2) coord(3)];
IP(23,:) = [coord(2) coord(2) coord(3)];
IP(24,:) = [coord(3) coord(2) coord(3)];
IP(25,:) = [coord(1) coord(3) coord(3)];
IP(26,:) = [coord(2) coord(3) coord(3)];
IP(27,:) = [coord(3) coord(3) coord(3)];
for element_number = 1:total_no_element;
%connectivity vector 'g' have no of nodes for each element
%making a local dof vector which has global dofs assigned for each element
% for element_number = 1:total_no_element;
g = conn(element_number,2:21);
for i = 1:20;
h(i) = xx(g(i));
y(i) = yy(g(i));
z(i) = zz(g(i));
a = (g(i)*3) ;
dof(:,1,i,element_number) = [a-2 a-1 a];
m(i).m = [a-2 a-1 a];
end
P = [m(1).m m(2).m m(3).m m(4).m m(5).m m(6).m m(7).m m(8).m m(9).m...
m(10).m m(11).m m(12).m m(13).m m(14).m m(15).m m(16).m m(17).m...
m(18).m m(19).m m(20).m]';
d(element_number).dof = P;
Cn(element_number).c = [h';y';z'];
%Generating D matrix for 27 Gauss Points
for GP = 1:27;
% Dl = subs(Dl,{eps,eta,tau},{IP(GP,1),IP(GP,2),IP(GP,3)});
Dm = subs(Dl,{eps,eta,tau},{(IP(GP,1)),(IP(GP,2)),(IP(GP,3))});
Jac(element_number).IP(GP).j = Dm*Cn(element_number).c';
Dg(element_number).IP(GP).d = inv(Jac(element_number).IP(GP).j)*Dm;
x = Dg(element_number).IP(GP).d;
b = [x(1,1) 0 0;0 x(2,1) 0; 0 0 x(3,1);x(2,1) x(1,1) 0;0 x(3,1) x(2,1);x(3,1) 0 x(1,1)];
B = [b b b b b b b b b b b b b b b b b b b b];
K(GP).k = double(B'*M*B*(det(Jac(element_number).IP(GP).j)));
end
phi = double(((weights(1))^3)*K(1).k + (weights(1)^2)*((weights(2))*K(2).k) ...
+ (weights(1)^2)*((weights(3))*K(3).k) + ...
(weights(1)^2)*((weights(2))*K(4).k) + ...
(weights(2)^2)*((weights(1))*K(5).k) + ...
(weights(1))*(weights(2))*(weights(3))*K(6).k +...
(weights(1)^2)*(weights(3))*K(7).k + ...
(weights(1))*(weights(2))*(weights(3))*K(8).k +...
(weights(3)^2)*(weights(1))*K(9).k + ...
(weights(1)^2)*(weights(2))*K(10).k + ...
(weights(2)^2)*(weights(1))*K(11).k + ...
(weights(1))*(weights(2))*(weights(3))*K(12).k +...
(weights(2)^2)*(weights(1))*K(13).k + ...
((weights(2))^3)*K(14).k + ...
(weights(2)^2)*(weights(3))*K(15).k + ...
(weights(1))*(weights(2))*(weights(3))*K(16).k +...
(weights(2)^2)*(weights(3))*K(17).k + ...
(weights(3)^2)*(weights(2))*K(18).k + ...
(weights(1)^2)*(weights(3))*K(19).k + ...
(weights(1))*(weights(2))*(weights(3))*K(20).k +...
(weights(3)^2)*(weights(1))*K(21).k + ...
(weights(1))*(weights(2))*(weights(3))*K(22).k +...
(weights(2)^2)*(weights(3))*K(23).k + ...
(weights(3)^2)*(weights(2))*K(24).k + ...
(weights(3)^2)*(weights(1))*K(25).k + ...
(weights(3)^2)*(weights(2))*K(26).k + ...
((weights(3))^3)*K(27).k);
K_local = double(phi);
%%%%%Assembly still has a problem
X(:,2) = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60];
X(:,1) = d(element_number).dof;
for u = 1:60;
for t = 1:60;
if Kg(X(u,1),X(t,1)) ==0;
Kg(X(u,1),X(t,1))= K_local(u,t);
else
Kg(X(u,1),X(t,1))= Kg(X(u,1),X(t,1))+ K_local(u,t);
end
end
end
end

Sign in to answer this question.

Answers (1)

John D'Errico
John D'Errico on 29 Jun 2014
It is virtually impossible to know what you have done wrong, without spending many hours tracing through your code, then deciphering what the boundary conditions were! I'd need to figure out exactly what you were doing. Did you do it all right? Did you make a silly mistake? Just a typo?
The first question might be to ask what your BC were. For example, suppose the BC are entirely derivative based? The your stiffness matrix MUST be singular, since then there is an arbitrary displacement factor that will not have been specified. You could freely translate the whole thing with no penalty.
But I really don't know (unless my guess is correct), and without far more clarity on your part, nobody else will probably make the effort.

Categories

Find more on Execution Speed in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!