三维八节点有限元编程问题

11713264515

用C语言写了一个程序求解三维八节点有限元问题，但是求解的单元刚度矩阵与ANSYS中的单元刚度矩阵不同，附上求解单元刚度矩阵的程序。望各位仁兄赐教，不慎感激！
void elem_severity_3d(double (*Ke)[24],int elem_id)
{
        int m,n,p,i,j,k;
        double J[3]={-0.7745966692,0,0.7745966692};
        double H[3]={0.5555555556,0.8888888889,0.5555555556};
        double Jc[9]={0};//Jc为三维的雅可比矩阵
        //计算Ke
        for(i=0;i<24;i++)
        {
                for(j=0;j<24;j++)
                {
                        Ke[i][j]=0;
                }
        }
        for(m=0;m<3;m++)
        {
                for(n=0;n<3;n++)
                {
                        for(p=0;p<3;p++)
                        {
                        double D[6][6]={0},B[6][24]={0},Bt[24][6]={0},DB[6][24]={0},BDB[24][24]={0};
                        double x[8]={0},y[8]={0},z[8]={0};double X,Y,Z,Jh;
                        X=J[m];Y=J[n];Z=J[p];
                        //计算D
                        D[5][5]=D[4][4]=D[3][3]=E/(2*(1+mu));
                        D[1][0]=D[1][2]=D[2][0]=D[2][1]=D[0][2]=D[0][1]=E*mu/((1+mu)*(1-2*mu));
                        D[1][1]=D[2][2]=D[0][0]=D[1][0]+2*D[3][3];
                        for(i=0;i<8;i++)
                        {
                                x[i]=node[9*(elem[8*elem_id+i]-1)];
                        }
                        for(i=0;i<8;i++)
                        {
                                y[i]=node[9*(elem[8*elem_id+i]-1)+1];
                        }
                        for(i=0;i<8;i++)
                        {
                                z[i]=node[9*(elem[8*elem_id+i]-1)+2];
                        }
                        double Na1,Na2,Na3,Na4,Na5,Na6,Na7,Na8,Nb1,Nb2,Nb3,Nb4,Nb5,Nb6,Nb7,Nb8,Nc1,Nc2,Nc3,Nc4,Nc5,Nc6,Nc7,Nc8;//Na,Nb,Nc分别为形函数N对局部坐标的偏导数
                        double xa,xb,xc,ya,yb,yc,za,zb,zc;//xa,xb,xc,ya,yb,yc,za,zb,zc分别为总体坐标x,y,z对局部坐标的偏导数
                       
                        Na1=-0.125*(1-Y)*(1-Z);Na2=-Na1;Na3=0.125*(1+Y)*(1-Z);Na4=-Na3;Na5=-0.125*(1-Y)*(1+Z);Na6=-Na5;Na7=0.125*(1+Y)*(1+Z);Na8=-Na7;
                        Nb1=-0.125*(1-X)*(1-Z);Nb2=-0.125*(1+X)*(1-Z);Nb3=-Nb2;Nb4=-Nb1;Nb5=-0.125*(1-X)*(1+Z);Nb6=-0.125*(1+X)*(1+Z);Nb7=-Nb6;Nb8=-Nb5;
                        Nc1=-0.125*(1-X)*(1-Y);Nc2=-0.125*(1+X)*(1-Y);Nc3=-0.125*(1+X)*(1+Y);Nc4=-0.125*(1-X)*(1+Y);Nc5=-Nc1;Nc6=-Nc2;Nc7=-Nc3;Nc8=-Nc4;
                        xa=Na1*x[0]+Na2*x[1]+Na3*x[2]+Na4*x[3]+Na5*x[4]+Na6*x[5]+Na7*x[6]+Na8*x[7];
                        xb=Nb1*x[0]+Nb2*x[1]+Nb3*x[2]+Nb4*x[3]+Nb5*x[4]+Nb6*x[5]+Nb7*x[6]+Nb8*x[7];
                        xc=Nc1*x[0]+Nc2*x[1]+Nc3*x[2]+Nc4*x[3]+Nc5*x[4]+Nc6*x[5]+Nc7*x[6]+Nc8*x[7];
                        ya=Na1*y[0]+Na2*y[1]+Na3*y[2]+Na4*y[3]+Na5*y[4]+Na6*y[5]+Na7*y[6]+Na8*y[7];
                        yb=Nb1*y[0]+Nb2*y[1]+Nb3*y[2]+Nb4*y[3]+Nb5*y[4]+Nb6*y[5]+Nb7*y[6]+Nb8*y[7];
                        yc=Nc1*y[0]+Nc2*y[1]+Nc3*y[2]+Nc4*y[3]+Nc5*y[4]+Nc6*y[5]+Nc7*y[6]+Nc8*y[7];
                        za=Na1*z[0]+Na2*z[1]+Na3*z[2]+Na4*z[3]+Na5*z[4]+Na6*z[5]+Na7*z[6]+Na8*z[7];
                        zb=Nb1*z[0]+Nb2*z[1]+Nb3*z[2]+Nb4*z[3]+Nb5*z[4]+Nb6*z[5]+Nb7*z[6]+Nb8*z[7];
                        zc=Nc1*z[0]+Nc2*z[1]+Nc3*z[2]+Nc4*z[3]+Nc5*z[4]+Nc6*z[5]+Nc7*z[6]+Nc8*z[7];
                        double Nabc1[3]={Na1,Nb1,Nc1},Nabc2[3]={Na2,Nb2,Nc2},Nabc3[3]={Na3,Nb3,Nc3},Nabc4[3]={Na4,Nb4,Nc4},Nabc5[3]={Na5,Nb5,Nc5},Nabc6[3]={Na6,Nb6,Nc6},Nabc7[3]={Na7,Nb7,Nc7},Nabc8[3]={Na8,Nb8,Nc8};
                        double Nxyz1[3]={0},Nxyz2[3]={0},Nxyz3[3]={0},Nxyz4[3]={0},Nxyz5[3]={0},Nxyz6[3]={0},Nxyz7[3]={0},Nxyz8[3]={0};//Nxyz[0]、Nxyz[1]、Nxyz[2]分别为形函数N对总体坐标x、y、z的偏导数
                        Jc[0]=xa;Jc[1]=xb;Jc[2]=xc;Jc[3]=ya;Jc[4]=yb;Jc[5]=yc;Jc[6]=za;Jc[7]=zb;Jc[8]=zc;
                       
                        Jh=Jc[0]*Jc[4]*Jc[8]+Jc[3]*Jc[7]*Jc[2]+Jc[6]*Jc[1]*Jc[5]-Jc[2]*Jc[4]*Jc[6]-Jc[1]*Jc[3]*Jc[8]-Jc[0]*Jc[5]*Jc[7];//Jh为雅可比矩阵Jc的行列式
                        mx_inver(Jc,3);//对Jc矩阵求逆
                        for(i=0;i<3;i++)
                        {
                                for(j=0;j<3;j++)
                                {
                                    Nxyz1[i]+=Jc[i+3*j]*Nabc1[j];
                                        Nxyz2[i]+=Jc[i+3*j]*Nabc2[j];
                                        Nxyz3[i]+=Jc[i+3*j]*Nabc3[j];
                                        Nxyz4[i]+=Jc[i+3*j]*Nabc4[j];
                                        Nxyz5[i]+=Jc[i+3*j]*Nabc5[j];
                                        Nxyz6[i]+=Jc[i+3*j]*Nabc6[j];
                                        Nxyz7[i]+=Jc[i+3*j]*Nabc7[j];
                                        Nxyz8[i]+=Jc[i+3*j]*Nabc8[j];
                                }       
                        }
                        B[0][0]=Nxyz1[0];B[1][1]=Nxyz1[1];B[2][2]=Nxyz1[2];B[3][1]=B[5][2]=B[0][0];B[3][0]=B[4][2]=B[1][1];B[4][1]=B[5][0]=B[2][2];
                        B[0][3]=Nxyz2[0];B[1][4]=Nxyz2[1];B[2][5]=Nxyz2[2];B[3][4]=B[5][5]=B[0][3];B[3][3]=B[4][5]=B[1][4];B[4][4]=B[5][3]=B[2][5];
                        B[0][6]=Nxyz3[0];B[1][7]=Nxyz3[1];B[2][8]=Nxyz3[2];B[3][7]=B[5][8]=B[0][6];B[3][6]=B[4][8]=B[1][7];B[4][7]=B[5][6]=B[2][8];
                        B[0][9]=Nxyz4[0];B[1][10]=Nxyz4[1];B[2][11]=Nxyz4[2];B[3][10]=B[5][11]=B[0][9];B[3][9]=B[4][11]=B[1][10];B[4][10]=B[5][9]=B[2][11];
                        B[0][12]=Nxyz5[0];B[1][13]=Nxyz5[1];B[2][14]=Nxyz5[2];B[3][13]=B[5][14]=B[0][12];B[3][12]=B[4][14]=B[1][13];B[4][13]=B[5][12]=B[2][14];
                        B[0][15]=Nxyz6[0];B[1][16]=Nxyz6[1];B[2][17]=Nxyz6[2];B[3][16]=B[5][17]=B[0][15];B[3][15]=B[4][17]=B[1][16];B[4][16]=B[5][15]=B[2][17];
                        B[0][18]=Nxyz7[0];B[1][19]=Nxyz7[1];B[2][20]=Nxyz7[2];B[3][19]=B[5][20]=B[0][18];B[3][18]=B[4][20]=B[1][19];B[4][19]=B[5][18]=B[2][20];
                        B[0][21]=Nxyz8[0];B[1][22]=Nxyz8[1];B[2][23]=Nxyz8[2];B[3][22]=B[5][23]=B[0][21];B[3][21]=B[4][23]=B[1][22];B[4][22]=B[5][21]=B[2][23];
                for(i=0;i<6;i++)
                {
                        for(j=0;j<24;j++)
                       {
                              Bt[j][i]=B[i][j];
                       }
                }
                        //计算DB
               for(i=0;i<6;i++)
               {
                      for(j=0;j<24;j++)
                      {
                             for(k=0;k<6;k++)
                             {
                                    DB[i][j]+=D[i][k]*B[k][j];
                             }
                      }
               }
               //计算BDB
               for(i=0;i<24;i++)
               {
                      for(j=0;j<24;j++)
                      {
                             for(k=0;k<6;k++)
                             {
                                    BDB[i][j]+=Bt[i][k]*DB[k][j];
                             }
                      }
                }
                   for(i=0;i<24;i++)
                        {
                                for(j=0;j<24;j++)
                                {
                                        Ke[i][j]+=BDB[i][j]*Jh*H[m]*H[n]*H[p];
                                }
                        }
                        }       
                }
        }
       
        }
}

pasuka

1、计算结果不同很正常啊！一模一样，分毫不差才邪门
2、只用3X3X3高斯积分？依据或出处在哪里？

Iorifo

ANSYS默认会是非协调或者B-bar单元，另外你用的是3*3*3的积分点，如果element formulation完全一致，在编译器，操作系统完全相同的情况下，结果可以丝毫不差。

骞耕

ANSYS提供了多种单元技术，如果选择的是完全积分的话，将使用B-bar方法求解，而楼主的程序中似乎没有B-bar修正。顺便在这里问一下，采用完全积分时，ABAQUS和ANSYS的高阶六面体单元的结果也是不一样的，应该是ANSYS采用了特殊技巧，向大家请教一下有没有知道这种技巧的。

pasuka

骞耕 发表于 2015-10-8 13:14
ANSYS提供了多种单元技术，如果选择的是完全积分的话，将使用B-bar方法求解，而楼主的程序中似乎没有B-bar ...

ANSYS的理论手册已经给出B-Bar单元的参考文献，再顺藤摸瓜网上找找不就行了，这个是从网上扒来的matlab代码，对照Hughes的论文或专著花几天功夫看懂应该不难

骞耕

pasuka 发表于 2015-10-8 15:10
ANSYS的理论手册已经给出B-Bar单元的参考文献，再顺藤摸瓜网上找找不就行了，这个是从网上扒来的matlab ...

请问您使用FEM_Bar.m程序计算过高阶六面体单元的案例吗？ 跟ANSYS能对的上吗？

skyismine

引入威尔逊非协调项再试试

pasuka

骞耕 发表于 2015-10-10 10:07
请问您使用FEM_Bar.m程序计算过高阶六面体单元的案例吗？ 跟ANSYS能对的上吗？ ...

可以对上，原理其实不复杂