why i get this ERROR???

Question

0 个投票

V

2 个评论
显示无隐藏无

ammar ahmed 2020-12-8

clc;
clear all;
E=200e9;
p=0.3
d=1;
t=1; 
L=5;
endload=-1000
nmeshx=32;
nmeshy=8;
%........%
Xspacing=(L/nmeshx)/2;
yspacing=(d/nmeshy)/2;
nel=nmeshx*nmeshy;
nnode=(nmeshx*2+1)*(nmeshy+1)+nmeshy*(nmeshx+1);
%.........coordinates 
nodecoord=zeros(nnode,2);
c=0;
for i=1:(nmeshy*2+1)
    if mod(i,2)==1
        for j=1:(2*nmeshx+1)
            c=c+1;
      xcoord=Xspacing*(j-1);
            ycoord=yspacing*(i-1);
            nodecoord(c,:) = [xcoord ycoord];
        end
    else
    for j=1:(nmeshx+1)
        c=c+1;
        xcoord=(2*Xspacing)*(j-1);
        ycoord=yspacing*(i-1);
        nodecoord(c,:)=[xcoord ycoord];
      end
    end
end
%..........conenectivy
conn=zeros(nel,4);
rdiff=nmeshx*2+nmeshx+2;
rdiffsmall=nmeshx*2+1;
c=0;
for j=0:nmeshy-1
    for i=0:nmeshx-1
        c=c+1;
        c1=(2*i+1)+rdiff*j;
        c2=c1+2;
        c4=c1+rdiff;
        c3=c4+2;
        conn(c,:)=[c1;c2;c3;c4];
    end
end
%..........load matrix
Load = zeros(nmeshy*2+1,3);
laod(1,:)=[rdiffsmall 0 endload/(nmeshy*2+1)];
for i=1:nmeshy*2
    Load(i+1,:)=[Load(i,1)+round((nmeshx+1)*(abs(cos((i+1)*pi()/2)))+rdiffsmall*(abs(cos(i*pi()/2 )))) 0 endload/(nmeshy*2+1)];
end
nLoad=size(Load,1);
%..................support
s=zeros(nmeshy*2+1,3);
s(1,:)=[1 1 1];
for i=1:nmeshy*2                 %s=[nodelid ux uy ] 1..fixed 0...free
    s(i+1,:)=[Load(i,1)+ 1 1 1];
end
numsupport=size(s,1);
c=0;
for i=1:nel
     nodecoordelem=zeros(4,2);
     for j=1:4
   nodecoordelem(j,:)=nodecoord(conn(i,j),:);
     end
     c=c+1;
     nodecoordelemlist(c).nodecoordelem=nodecoordelem;
end
EM=E/(1-p^2).*[1 p 0;p 1 0; 0 0 (1-p)/2];
%shape functions
syms ks ta
    n1 = 0.25*(1-ks)*(1-ta);
    n2 = 0.25*(1+ks)*(1-ta);
    n3 = 0.25*(1+ks)*(1+ta);
    n4 = 0.25*(1-ks)*(1+ta);
N8=[n1;n2;n3;n4];
dnks = [-0.25*(1-ta) 0.25*1*(1-ta)  0.25*1*(1+ta)  0.25*-1*(1+ta)];
dnta = [-0.25*(1-ks) 0.25*-1*(1+ks) 0.25*1*(1+ks) 0.25*(1-ks)];
D8=[dnks 
    dnta];
%element stiffness matrix
for g=1:nel
jac=D8*nodecoordelemlist(g).nodecoordelem; %jacobian matrix
capgamma = inv(jac);      %inverse jacobian matrix
B=sym(zeros(3,8)); 
for k=0:3
m=2*k+1;
n=m+1;
B(1,m) = capgamma(1,1)*dnks(k+1)+capgamma(1,2)*dnta(k+1);
B(2,n) = capgamma(2,1)*dnks(k+1)+capgamma(2,2)*dnta(k+1);
B(3,m) = B(2,n);
B(3,n) = B(1,m);
end
kelemlist(g).kelem = int(int(transpose(B)*EM*B*t*det(jac),ks,-1,1),ta,-1,1);
                                 
kelemlist(g).B = B;
end
%global stiffness matrix
kk = zeros(nnode*2,nnode*2);
for k=1:nel
 kelem=kelemlist(k).kelem;
for j=1:4
 for i=1:4
    n=conn(k,i);
    m=conn(k,j);
kk(2*n-1,2*m-1)=kk(2*n-1,2*m-1)+kelem(2*i-1,2*j-1);
kk(2*n-1,2*m)=kk(2*n-1,2*m)+kelem(2*i-1,2*j);
kk(2*n,2*m-1)=kk(2*n,2*m-1)+kelem(2*i,2*j-1);
kk(2*n,2*m)=kk(2*n,2*m)+kelem(2*i,2*j);
end
end
end
%apply the boundry condition 
kkWosup=kk;
for i=1:numsupport
    n=s(i,1);
    if(s(i,2)==1)
        kk(2*n-1,:)=0;
        kk(:,2*n-1)=0;
        kk(2*n-1,2*n-1)=1;
    end
    if(s(i,3)==1)
        kk(2*n,:)=0;
        kk(:,2*n)=0;
        kk(2*n,2*n)=1;
        end
end
    %load vector
    Loadvector = zeros(nnode*2,1);
for i=1:nLoad
n =Load(i,1);
Loadvector(2*n-1)=Load(i,2);
Loadvector(2*n) = Load(i,3);
end
%displacement computation
u = kk\Loadvector;
for i=1:nnode
unodelist(i).node = [u(2*i-1);u(2*i)];
end
for i=1:nel
temp = [];
for j=1:4
temp = [temp ; unodelist(conn(i,j)).node];
end
uelemlist(i).elem = temp;
end
beamtheory = (endload*L^3)/(3*E*(t*d^3/12)+(endload*L))/((E/2*(1-p))*(5/6)*d*t);
Stress=beamtheory/(2*Xspacing*yspacing*-1e-2);
%stress computation
scp = [-1 -1;1 -1;1 1;-1 1];
    for i=1:nel
uelem = uelemlist(i).elem';
temp = [];
felemlist(i).felem = double(kelemlist(i).kelem*uelem');
sigma = EM*kelemlist(i).B*uelem';
temp = [];
for j=1:4
temp = [temp double(subs(subs(sigma,ks,scp(j,1)),ta,scp(j,2)))];
end
sigmaelemlist(i).sigma = temp;
 end
fkk = kkWosup*u;

ammar ahmed 2020-12-8

the code is above pleace help me as fast as you can i was like spending to much time trying to figur out the the problem

请先登录，再进行评论。

请先登录，再回答此问题。

Follow Question

Answer 1

Star Strider 2020-12-8

0 个投票

I suspect it is because ‘(2*n-1)’ is less than or equal to 0.

I have no idea what ‘n’ is, so it could also be the problem if it is not an integer.

Since apparently ‘Load’ is not throwing an error, I suspect it is a function you (or someone else) wrote. The MATLAB function has a lower-case ‘L’, and since MATLAB is case-sensitive, the difference is not the problem.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Walter Roberson 2020-12-8

编辑：Walter Roberson 2020-12-8

在 MATLAB Online 中打开

As outside observers, we have no reason to expect that the kk matrix you are creating will be non-singular.

Typically when a stiffness matrix is singular, the implication is that there are too many degrees of freedom, and that additional entries need to be added to pin down some part that is too free to move.

>> size(kk)
ans =
        1698        1698
>> rank(kk)
ans =
   591

Not even close :( You are missing pinning down a lot of values !

Star Strider 2020-12-8

Steven — I didn’t see ‘n’ in the screencap. Thank you for discovering it. (My Answer remains valid.)

Walter — Thank you!

请先登录，再进行评论。

why i get this ERROR???

2 个评论
显示无隐藏无

回答（1 个）

7 个评论
显示 5更早的评论隐藏 5更早的评论

类别

标签

Community Treasure Hunt

why i get this ERROR???

2 个评论 显示 无 隐藏 无

回答（1 个）

7 个评论 显示 5更早的评论 隐藏 5更早的评论

类别

标签

另请参阅

Community Treasure Hunt

2 个评论
显示无隐藏无

7 个评论
显示 5更早的评论隐藏 5更早的评论