单元刚度矩阵和节点编号顺序的关系？

xiejiawei110

在做 曾攀的《有限元分析基础教程》一书中“基于3节点三角形单元的矩形薄板分析”分析中发现按照书中给的单元刚度程序：

function k=Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)
%该函数计算单元的刚度矩阵
%输入弹性模量E，泊松比NU，厚度t
%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym
%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)
%输出单元刚度矩阵k(6X6)
A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj))/2;
betai = yj-ym;
betaj = ym-yi;
betam = yi-yj;
gammai = xm-xj;
gammaj = xi-xm;
gammam = xj-xi;
B = [betai 0 betaj 0 betam 0 ;
0 gammai 0 gammaj 0 gammam ;
gammai betai gammaj betaj gammam betam]/(2*A);
if ID == 1
D = (E/(1-NU*NU))*[1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2];
elseif ID == 2
D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2];
end
k= t*A*B'*D*B;
end

按照节点2,3,4的顺序调取得到单刚

E=1E7;
NU=1/3;
t=0.1;
ID=1;
KK=zeros(8,8);

k1=Triangle2D3Node_Stiffness(E,NU,t,2,0,0,1,0,0,ID)
得到
k1 =

   1.0e+06 *

    0.2813         0         0    0.1875   -0.2813   -0.1875
         0    0.0938    0.1875         0   -0.1875   -0.0938
         0    0.1875    0.3750         0   -0.3750   -0.1875
    0.1875         0         0    1.1250   -0.1875   -1.1250
   -0.2813   -0.1875   -0.3750   -0.1875    0.6563    0.3750
   -0.1875   -0.0938   -0.1875   -1.1250    0.3750    1.2188

和按照节点4,3,2的顺序得到的单刚的符号是完全相反的
按照节点4,3,2的顺序调取得到单刚
k1 =


   1.0e+06 *


   -0.6563   -0.3750    0.3750    0.1875    0.2813    0.1875
   -0.3750   -1.2188    0.1875    1.1250    0.1875    0.0938
    0.3750    0.1875   -0.3750         0         0   -0.1875
    0.1875    1.1250         0   -1.1250   -0.1875         0
    0.2813    0.1875         0   -0.1875   -0.2813         0
    0.1875    0.0938   -0.1875         0         0   -0.0938

单刚的符号不同，以至于组装后总刚矩阵也不一样，代入边界条件后求得的应力和节点位移等都不一样。
百思不得其解，同一个单元的单元刚度矩阵难道还和输入节点的顺序有关吗？真心请教各位,谢谢！

penultimate

是呀，你算出来的有向面积都是负数的，一定要是逆时针方向输入节点

xiejiawei110

penultimate 发表于 2017-10-19 15:08
是呀，你算出来的有向面积都是负数的，一定要是逆时针方向输入节点

还真有这个讲究呀，我学的曾攀那本书里面完全没有提到，怪不得呢，谢谢啦！