导热问题有限元刚度矩阵三维立方体8节点

导热问题有限元刚度矩阵三维立方体8节点

三维立方体8节点等参有限元的热传导刚度矩阵计算很复杂，这里直接给了计算结果

采用求导积分方法的matlab程序如下，DH是xyz三个方向的导热系数[kx,ky,kz]

function KE = elementMatVec3DNew(DH)
a=0.5; b=0.5;c=0.5;
k = [DH(1) 0 0;0 DH(2) 0;0 0 DH(3)];
u1=-a;v1=-b;w1=-c;
u2=a;v2=-b;w2=-c;
u3=a;v3=b;w3=-c;
u4=-a;v4=b;w4=-c;
u5=-a;v5=-b;w5=c;
u6=a;v6=-b;w6=c;
u7=a;v7=b;w7=c;
u8=-a;v8=b;w8=c;

                          
syms x y z;
N1=(1/8)*(1+(u1./a)*(x./a))*(1+(v1./b)*(y./b))*(1+(w1./c)*(z./c));
N2=(1/8)*(1+(u2/a)*(x./a))*(1+(v2/b)*(y./b))*(1+(w2/c)*(z./c));
N3=(1/8)*(1+(u3/a)*(x./a))*(1+(v3/b)*(y./b))*(1+(w3/c)*(z./c));
N4=(1/8)*(1+(u4/a)*(x./a))*(1+(v4/b)*(y./b))*(1+(w4/c)*(z./c));
N5=(1/8)*(1+(u5/a)*(x./a))*(1+(v5/b)*(y./b))*(1+(w5/c)*(z./c));
N6=(1/8)*(1+(u6/a)*(x./a))*(1+(v6/b)*(y./b))*(1+(w6/c)*(z./c));
N7=(1/8)*(1+(u7/a)*(x./a))*(1+(v7/b)*(y./b))*(1+(w7/c)*(z./c));
N8=(1/8)*(1+(u8/a)*(x./a))*(1+(v8/b)*(y./b))*(1+(w8/c)*(z./c));

B=[diff(N1,x) diff(N2,x) diff(N3,x)  diff(N4,x) diff(N5,x) diff(N6,x)  diff(N7,x) diff(N8,x);
   diff(N1,y) diff(N2,y) diff(N3,y)  diff(N4,y) diff(N5,y) diff(N6,y)  diff(N7,y) diff(N8,y);
   diff(N1,z) diff(N2,z) diff(N3,z)  diff(N4,z) diff(N5,z) diff(N6,z)  diff(N7,z) diff(N8,z)];
KE=double(int(int(int((B'*k*B),x,-a,a), y,-b,b),z,-c,c));
end

对于节点距离为1的标准格式

计算过程如下

function KE = elementMatVec3DNew2(DH)

KE1=[4 -4 -2 2 2 -2 -1 1;
    -4 4 2 -2 -2 2 1 -1;
    -2 2 4 -4 -1 1 2 -2;
    2 -2 -4 4 1 -1 -2 2;
    2 -2 -1 1 4 -4 -2 2;
    -2 2 1 -1 -4 4 2 -2;
    -1 1 2 -2 -2 2 4 -4;
    1 -1 -2 2 2 -2 -4 4];
KE2=[4 2 -2 -4 2 1 -1 -2;
    2 4 -4 -2 1 2 -2 -1;
    -2 -4 4 2 -1 -2 2 1;
    -4 -2 2 4 -2 -1 1 2;
    2 1 -1 -2 4 2 -2 -4;
    1 2 -2 -1 2 4 -4 -2;
    -1 -2 2 1 -2 -4 4 2;
    -2 -1 1 2 -4 -2 2 4;];
KE3=[4 2 1 2 -4 -2 -1 -2;
    2 4 2 1 -2 -4 -2 -1;
    1 2 4 2 -1 -2 -4 -2;
    2 1 2 4 -2 -1 -2 -4;
    -4 -2 -1 -2 4 2 1 2;
    -2 -4 -2 -1 2 4 2 1;
    -1 -2 -4 -2 1 2 4 2;
    -2 -1 -2 -4 2 1 2 4;];
KE = (KE1*DH(1)+KE2*DH(2)+KE3*DH(3))/36;
end

---

转账请注明来源，欢迎评论区讨论