查看: 735|回复: 3|关注: 0

[已解决] 警告: 矩阵为奇异值、接近奇异值或缩放错误。结果可能不准确。RCOND = NaN。 > In Tri_main (line 179)

[复制链接]

新手

9 麦片

财富积分


050


1

主题

3

帖子

0

最佳答案
本帖最后由 2463476757 于 2020-6-28 10:34 编辑

   % 位移边界条件的处理---对角元素改1法
   Knew=K;
   Rnew=R;
   for i=1:length(d)
    if d(i)==0
        Knew(i,:)=0;
        Knew(:,i)=0;
        Rnew(i)=0;
        Knew(i,i)=1;
    end
   end

   %求解位移
179   d=inv(Knew)*Rnew;

%---------------------------------------------------------输出节点位移
disp('----------------------------------------节点位移');
警告: 矩阵为奇异值、接近奇异值或缩放错误。结果可能不准确。RCOND = NaN。
> In Tri_main (line 179)  d=inv(Knew)*Rnew;

                               
登录/注册后可看大图


Node.docx

16.99 KB, 下载次数: 0

回复主题 已获打赏: 0 积分

举报

论坛优秀回答者

9

主题

1630

帖子

342

最佳答案
  • 关注者: 81
发表于 2020-6-28 09:19:10 | 显示全部楼层
  1. d=inv(Knew)*Rnew;
复制代码

改成
  1. d=pinv(Knew)*Rnew;
复制代码

若不行请上传完整代码
回复此楼 已获打赏: 0 积分

举报

新手

9 麦片

财富积分


050


1

主题

3

帖子

0

最佳答案
 楼主| 发表于 2020-6-28 10:33:35 | 显示全部楼层
20141303 发表于 2020-6-28 09:19

改成
若不行请上传完整代码

E=200000;u=0.33;   %材料参数
t=1;%几何参数
Fax=1200;Fay=1000;Fbx=800;Fby=750;%载荷边界条件
Nodes=dlmread('Nodes2.txt');   %读取节点坐标
Elements=dlmread('Elements2.txt');    %读取单元编号及其相应节点


n=length(Nodes(:,1)); %节点个数
m=length(Elements(:,1));  %单元个数

%节点力

R=zeros(2*n,1);
R(17)=Fbx;
R(18)=Fby;
R(23)=Fax;
R(24)=Fay;

%节点位移
d=ones(2*n,1);
d(19:20)=0;
d(21:22)=0;



Num_n=n;
Num_e=m;


%-------------------------------------------------------------------------------------计算单元刚度及总体刚度
D=  E/(1-u^2)*[ 1  u  0;...
                           u 1   0;...
                           0  0 (1-u)/2];        %弹性矩阵
                       
K=zeros(2*Num_n,2*Num_n);      %构建总体2n*2n的空矩阵
for i=1:Num_n
%-----------------------------------------------------计算单元刚度
    ni=Elements(i,1);     %调取单元的i, j, m节点
    nj=Elements(i,2);
    nm=Elements(i,3);
   
    xi=Nodes(ni,1);  %i节点的坐标
    yi=Nodes(ni,2);
   
    xj=Nodes(nj,1);   % j节点的坐标
    yj=Nodes(nj,2);
   
    xm=Nodes(nm,1);  % m节点的坐标
    ym=Nodes(nm,2);
   
    ai=xj*ym-xm*yj;   % 分别计算ai,bi,ci   (3-5)
    bi=yj-ym;
    ci=-xj+xm;
   
    aj=xm*yi-xi*ym;
    bj=ym-yi;
    cj=-xm+xi;
   
    am=xi*yj-xj*yi;
    bm=yi-yj;
    cm=-xi+xj;

    A=1/2*det([1  xi     yi;...           % 计算面积   (3-6)
                        1  xj     yj;...
                        1  xm   ym]);

   
     
     B=1/(2*A)*[bi   0   bj  0   bm    0;...
                         0    ci   0   cj   0      cm;...
                         ci    bi  cj   bj   cm   bm];           %应变矩阵 (3-14)     
                     
  
     
     S=D*B;    %应力矩阵 (3-17)   
      

     SS{i}=S; % 将单元应力矩阵赋值给SS,用于后续单元应力计算
     
     Ke=B'*S*t*A;   % 刚度矩阵 (3-23)   
%-----------------------------------------------------计算单元刚度结束

%-----------------------------------------------------计算总体刚度  P34


   
   Ki=zeros(2*Num_n,2*Num_n);  %将单元刚度6*6矩阵扩展为2n*2n的空矩阵
   
    Ki(2*ni-1:2*ni,2*ni-1:2*ni)=Ke(1:2,1:2);  %按照自由度将单元刚度矩阵Ke填写至扩展后的刚度矩阵Ki
    Ki(2*ni-1:2*ni,2*nj-1:2*nj)=Ke(1:2,3:4);
    Ki(2*ni-1:2*ni,2*nm-1:2*nm)=Ke(1:2,5:6);
   
    Ki(2*nj-1:2*nj,2*ni-1:2*ni)=Ke(3:4,1:2);
    Ki(2*nj-1:2*nj,2*nj-1:2*nj)=Ke(3:4,3:4);
    Ki(2*nj-1:2*nj,2*nm-1:2*nm)=Ke(3:4,5:6);
   
    Ki(2*nm-1:2*nm,2*ni-1:2*ni)=Ke(5:6,1:2);
    Ki(2*nm-1:2*nm,2*nj-1:2*nj)=Ke(5:6,3:4);
    Ki(2*nm-1:2*nm,2*nm-1:2*nm)=Ke(5:6,5:6);
   
    K=K+Ki; %将所有单元刚度矩阵相加得总刚
   
%     disp(['单元号',num2str(i)]);
%     bijm=[bi,bj,bm];
%     cijm=[ci,cj,cm];
   
    disp(['单元号',num2str(i)]);
    Ke;
   
end
disp(['总刚']);
K;


   % 位移边界条件的处理---对角元素改1法
   Knew=K;
   Rnew=R;
   for i=1:length(d)
    if d(i)==0
        Knew(i,:)=0;
        Knew(:,i)=0;
        Rnew(i)=0;
        Knew(i,i)=1;
    end
   end

   %求解位移
d=pinv(Knew)*Rnew;

%---------------------------------------------------------输出节点位移
disp('----------------------------------------节点位移');
for i=1:Num_n
    disp(['节点号',num2str(i),'  x方向位移:',num2str(d(2*i-1,1)),';  y方向位移:',num2str(d(2*i,1))]);
end
uv(:,1)=d(1:2:end);  uv(:,2)=d(2:2:end);
dsum=sqrt(uv(:,1).^2+uv(:,2).^2);
dmax=max(dsum);
index=find(dsum==dmax);
disp(['最大位移为:',num2str(dmax),';  最大位移节点为:',num2str(index)]);

%--------------------------------------------------------------------------------------------------------------------画变形后的图
nnew=Nodes+max(abs(Nodes(:,1)))/max(abs(d))/15*uv;
plot(nnew(:,1),nnew(:,2),'blao')

for i=1:Num_e
    ele=Elements(i,:);  
    ele(:,end+1)=ele(:,1);
    xx=nnew(ele(1,:),1);
    yy=nnew(ele(1,:),2);
    plot(xx,yy,'bla-.')
end
%--------------------------------------------------------------------------------------------------------------------画变形后的图

% 求解单元应力
for i=1:Num_e
    ni=Elements(i,1); %调取单元的i, j, m节点
    nj=Elements(i,2);
    nm=Elements(i,3);
    delta=[d(2*ni-1),d(2*ni),d(2*nj-1),d(2*nj),d(2*nm-1),d(2*nm)]';  % 单元节点位移
    S=SS{i};  % 单元应力矩阵
    sigma=S*delta;   % 单元应力
    s1(i,:)=sigma';  % 储存第i个单元应力
end

disp('----------------------------------------单元应力');
for i=1:Num_e
    disp(['单元号',num2str(i),' sigmax:',num2str(s1(i,1)),';   sigmay:',num2str(s1(i,2)),';  tao:',num2str(s1(i,3))]);
end

% 绘制应力云图
% for i=1:Num_e
%    ijm=Elements(i,:);
%     fill(Nodes(ijm,1),Nodes(ijm,2),s1(i,1)*ones(3,1))
% end
% shading interp;
% colorbar;

%绕节点平均法
for i=1:Num_n
    s_node(i,1:3)=0; %存储第i个节点的总应力
    n=0;  %计算节点周围单元个数
    for j=1:Num_e  %遍历所有单元
        index=find(Elements(j,:)==i);
        if isempty(index)
        else
            n=n+1; %计算节点周围单元个数
            s_node(i,1:3)=s1(j,1:3)+s_node(i,1:3);
        end
    end
    s_node(i,1:3)=s_node(i,1:3)/n;  %求平均
end

%---------------------------------------------------------输出绕节点平均后的节点应力
disp('----------------------------------------绕节点平均后的节点应力');
for i=1:Num_n
    disp(['节点号',num2str(i),' sigmax:',num2str(s_node(i,1)),';   sigmay:',num2str(s_node(i,2)),';  tao:',num2str(s_node(i,3))]);
end

% 节点平均后的应力云图
for i=1:Num_e
    ijm=Elements(i,:);  
    fill(Nodes(ijm,1),Nodes(ijm,2),s_node(ijm,1))
end
shading interp;
colorbar;     
错误使用 svd
SVD 的输入不能包含 NaN 或 Inf。

出错 pinv (line 18)
[U,S,V] = svd(A,'econ');

出错 Tri_main (line 179)
d=pinv(Knew)*Rnew;
回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

9

主题

1630

帖子

342

最佳答案
  • 关注者: 81
发表于 2020-6-28 11:42:07 | 显示全部楼层 |此回复为最佳答案
前面部分修改如下试试
  1. E=200000;u=0.33;   %材料参数
  2. t=1;%几何参数
  3. Fax=1200;Fay=1000;Fbx=800;Fby=750;%载荷边界条件
  4. Nodes=dlmread('Nodes2.txt');   %读取节点坐标
  5. Elements=dlmread('Elements2.txt');    %读取单元编号及其相应节点


  6. n=length(Nodes(:,1)); %节点个数
  7. m=length(Elements(:,1));  %单元个数

  8. %节点力

  9. R=zeros(2*n,1);
  10. R(17)=Fbx;
  11. R(18)=Fby;
  12. R(23)=Fax;
  13. R(24)=Fay;

  14. %节点位移
  15. d=ones(2*n,1);
  16. d(19:20)=0;
  17. d(21:22)=0;



  18. Num_n=n;
  19. Num_e=m;


  20. %-------------------------------------------------------------------------------------计算单元刚度及总体刚度
  21. D=  E/(1-u^2)*[ 1  u  0;...
  22.                            u 1   0;...
  23.                            0  0 (1-u)/2];        %弹性矩阵
  24.                        
  25. K=zeros(2*Num_n,2*Num_n);      %构建总体2n*2n的空矩阵
  26. for i=1:Num_n
  27. %-----------------------------------------------------计算单元刚度
  28.     ni=Elements(i,1);     %调取单元的i, j, m节点
  29.     nj=Elements(i,2);
  30.     nm=Elements(i,3);
  31.    
  32.     xi=Nodes(ni,1);  %i节点的坐标
  33.     yi=Nodes(ni,2);
  34.    
  35.     xj=Nodes(nj,1);   % j节点的坐标
  36.     yj=Nodes(nj,2);
  37.    
  38.     xm=Nodes(nm,1);  % m节点的坐标
  39.     ym=Nodes(nm,2);
  40.    
  41.     ai=xj*ym-xm*yj;   % 分别计算ai,bi,ci   (3-5)
  42.     bi=yj-ym;
  43.     ci=-xj+xm;
  44.    
  45.     aj=xm*yi-xi*ym;
  46.     bj=ym-yi;
  47.     cj=-xm+xi;
  48.    
  49.     am=xi*yj-xj*yi;
  50.     bm=yi-yj;
  51.     cm=-xi+xj;

  52.     A=1/2*det([1  xi     yi;...           % 计算面积   (3-6)
  53.                         1  xj     yj;...
  54.                         1  xm   ym]);

  55.    
  56.      
  57.      B=1/(2*A+eps)*[bi   0   bj  0   bm    0;...
  58.                          0    ci   0   cj   0      cm;...
  59.                          ci    bi  cj   bj   cm   bm];           %应变矩阵 (3-14)     
  60.                      
  61.   
  62.      
  63.      S=D*B;    %应力矩阵 (3-17)   
  64.       

  65.      SS{i}=S; % 将单元应力矩阵赋值给SS,用于后续单元应力计算
  66.      
  67.      Ke=B'*S*t*A;   % 刚度矩阵 (3-23)   
  68. %-----------------------------------------------------计算单元刚度结束

  69. %-----------------------------------------------------计算总体刚度  P34


  70.    
  71.    Ki=zeros(2*Num_n,2*Num_n);  %将单元刚度6*6矩阵扩展为2n*2n的空矩阵
  72.    
  73.     Ki(2*ni-1:2*ni,2*ni-1:2*ni)=Ke(1:2,1:2);  %按照自由度将单元刚度矩阵Ke填写至扩展后的刚度矩阵Ki
  74.     Ki(2*ni-1:2*ni,2*nj-1:2*nj)=Ke(1:2,3:4);
  75.     Ki(2*ni-1:2*ni,2*nm-1:2*nm)=Ke(1:2,5:6);
  76.    
  77.     Ki(2*nj-1:2*nj,2*ni-1:2*ni)=Ke(3:4,1:2);
  78.     Ki(2*nj-1:2*nj,2*nj-1:2*nj)=Ke(3:4,3:4);
  79.     Ki(2*nj-1:2*nj,2*nm-1:2*nm)=Ke(3:4,5:6);
  80.    
  81.     Ki(2*nm-1:2*nm,2*ni-1:2*ni)=Ke(5:6,1:2);
  82.     Ki(2*nm-1:2*nm,2*nj-1:2*nj)=Ke(5:6,3:4);
  83.     Ki(2*nm-1:2*nm,2*nm-1:2*nm)=Ke(5:6,5:6);
  84.    
  85.     K=K+Ki; %将所有单元刚度矩阵相加得总刚
  86.    
  87.     K(1,1)
  88.     i
  89. %     disp(['单元号',num2str(i)]);
  90. %     bijm=[bi,bj,bm];
  91. %     cijm=[ci,cj,cm];
  92.    
  93.     disp(['单元号',num2str(i)]);
  94.     Ke;
  95.    
  96. end
  97. disp(['总刚']);
  98. K;


  99.    % 位移边界条件的处理---对角元素改1法
  100.    Knew=K;
  101.    Rnew=R;
  102.    for i=1:length(d)
  103.     if d(i)==0
  104.         Knew(i,:)=0;
  105.         Knew(:,i)=0;
  106.         Rnew(i)=0;
  107.         Knew(i,i)=1;
  108.     end
  109.    end

  110.    %求解位移
  111. d=pinv(Knew)*Rnew;
复制代码
回复此楼 已获打赏: 0 积分

举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /4 下一条

快速回复 返回顶部 返回列表