tet10 单元刚度矩阵的 FEM 实现：（提供的代码片段）

计算科学 有限元
2021-12-01 02:41:21

我正在尝试在 tet10 元素上实现 fem 代码。目前我不喜欢使用开源，因为我想对算法有基本的感觉。我密切关注 tet10 实施的讲义

http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch10.d/AFEM.Ch10.pdf

然而，当我在上面给出的示例四面体上测试我的代码时，我得到一个与文档中给出的非常不同的刚度矩阵。然而，我得到了我的刚度矩阵的特征值趋势的匹配。我在这里包含了一个测试代码来计算基本的刚度矩阵。

    #include<cstdlib>
    #include<iostream>
    #include<fstream>
    #include<ctime>
    #include<cmath>
    #include<string>
    #include<vector>
    #include<unistd.h>
    #include <iterator>
    #include <fstream>
    using namespace std;
    #include<Eigen/Core>
    #include<Eigen/Dense>
    #include<Eigen/Eigenvalues> 
    #include<Eigen/SVD>
    #include<Eigen/Sparse>
    #include<Eigen/SparseCholesky>
    #include<Eigen/SparseLU>
    using namespace Eigen;
    typedef Eigen::SparseMatrix<double> SpMat;
    struct NodeInfo
     {
      double x,y,z;
     };
    struct ElementInfo
    {
      int node[10];
     double volume, area[4],nx[4],ny[4],nz[4];
    };
   struct constants
   {
    double YoungsModulus, PoissonRatio, Density, gravity;
   };
   std::istream& operator>>(std::istream& is, ElementInfo& element)
   {
    is >> element.node[0] >> element.node[1] >> element.node[2] >>      element.node[3] 
   >> element.node[4] >> element.node[5] >> element.node[6] >> element.node[7]
   >> element.node[8] >> element.node[9] >> element.volume
   >> element.nx[0] >> element.ny[0] >>element.nz[0] >> element.area[0] 
   >> element.nx[1] >> element.ny[1] >>element.nz[1] >> element.area[1]
   >> element.nx[2] >> element.ny[2] >>element.nz[2] >> element.area[2]
   >> element.nx[3] >> element.ny[3] >>element.nz[3] >> element.area[3];
  return is;
  }
   std::istream& operator>>(std::istream& is, NodeInfo& node)
  {
    is >> node.x >>node.y >> node.z ;
    return is;
  }

  void ComputeElasticityMat(Eigen::MatrixXd& El, constants& M);
  void ComputeElementStiffness(ElementInfo& ithEle, vector<NodeInfo>& node,
  Eigen::MatrixXd& K, Eigen::MatrixXd& E, double (*eta)[4], double weight);

int main()
 {
char filename[] = "elements.dat";
char filename1[] ="nodes.dat";
int nelem, nnodes, TotalDOF;
int ndf=3;//number of degrees of freedom at each node
int npe=10;//number of points per element

// Gauss Integration points
double alfa=(5.0+3.0*sqrt(5.0))/20;
double beta =(5.0-sqrt(5.0))/20;
double weight=1.0/4.0;
double GaussPoints[4][4]={{alfa,beta,beta,beta},{beta,alfa,beta,beta},  {beta,beta,alfa,beta},{beta,beta,beta,alfa}};

// To define the element properties
constants MaterialProperties={.YoungsModulus=480.0,.PoissonRatio=1.0/3.0, .Density=2830, .gravity=9.81}; 
// Read the data file for the elements and the nodes
std::vector<ElementInfo> Elem;
std::vector<NodeInfo> Node;
std::ifstream ifs(filename), ifs1(filename1);
ifs>>nelem;
cout<<"The total number of elements is="<<nelem<<endl;
if (ifs) {
    std::copy(std::istream_iterator<ElementInfo>(ifs), 
            std::istream_iterator<ElementInfo>(),
            std::back_inserter(Elem));
}
else {
    std::cerr << "Couldn't open " << filename << " for reading\n";
}
cout<<"Element file read sucessfuly  "<<endl;

  if (ifs1) {
    std::copy(std::istream_iterator<NodeInfo>(ifs1), 
            std::istream_iterator<NodeInfo>(),
            std::back_inserter(Node));
}
else {
    std::cerr << "Couldn't open " << filename1 << " for reading\n";
}
    cout<<"Node file read sucessfuly  "<<endl;

  cout<<"The total number of Nodes is="<<Node.size() <<endl;
  nnodes=Node.size();

  //To form the Elasticity Matrix
  TotalDOF=nnodes*ndf;
  cout<<"The Total number of Degrees of freedom=:\n"<<TotalDOF<<endl;


   MatrixXd EelemMat(6,6), KelemMat(30,30),  KGl(TotalDOF, TotalDOF);
   VectorXd BodyForce(30), BFGl(TotalDOF), SurfaceForce(30);
   EelemMat.setZero(6,6);

   ComputeElasticityMat(EelemMat,MaterialProperties);
   cout<<"Elasticity Matrix=:\n"<<EelemMat<<endl;

   //To find the Stiffness Matrix
   if(nelem!=Elem.size()){
   cout<<"Something wrong in element file...please check!"<<endl;
   exit(0);
  }

 for(int i=0;i<nelem;i++)
   {
    KelemMat.setZero(30,30);
       ComputeElementStiffness(Elem[i],Node,KelemMat,
       EelemMat,GaussPoints,weight);

      cout<<"Stiffness Matrix="<<KelemMat<<endl;
   }
   JacobiSVD<MatrixXd> svd(KelemMat, ComputeThinU | ComputeThinV);
   return 0;    
 }   

 void ComputeElasticityMat(Eigen::MatrixXd& El, constants& M)
  {
   double K=M.YoungsModulus/((1.0+M.PoissonRatio)*  (1.0-2.0*M.PoissonRatio));
    El(0,0)=1.0-M.PoissonRatio;
    El(0,1)=M.PoissonRatio;
    El(0,2)=M.PoissonRatio;

    El(1,0)=M.PoissonRatio;
    El(1,1)=1.0-M.PoissonRatio;
    El(1,2)=M.PoissonRatio;

    El(2,0)=M.PoissonRatio;
    El(2,1)=M.PoissonRatio;
    El(2,2)=1.0-M.PoissonRatio;

    El(3,3)=0.5-M.PoissonRatio;
    El(4,4)=El(3,3);
    El(5,5)=El(3,3);

    El=El*K;

    return;
  }


  void ComputeElementStiffness(ElementInfo& ithEle, vector<NodeInfo>& node, Eigen::MatrixXd& K,
  Eigen::MatrixXd& E, double (*eta)[4], double weight)
 {
  double jx1,jx2,jx3,jx4,jy1,jy2,jy3,jy4,jz1,jz2,jz3,jz4,
       a1,a2,a3,a4,b1,b2,b3,b4,c1,c2,c3,c4,Jdet;
  VectorXd Nfx(10),Nfy(10),Nfz(10),xi(4);
  MatrixXd B(6,30), Bt; 
  MatrixXd J(4,4), P, Iaug(4,3),Jinv(4,4);
  double v01,v02,v03,v04,V;
 for(int intPoints=0;intPoints<4;intPoints++)
  {
    for(int j=0;j<4;j++)
       xi(j) = eta[intPoints][j];


       jx1=4.0*(node[ithEle.node[0]].x*(xi(0)-0.25)+node[ithEle.node[4]].x*xi(1)+node[ithEle.node[6]].x*xi(2)+node[ithEle.node[7]].x*xi(3));
       jy1=4.0*(node[ithEle.node[0]].y*(xi(0)-0.25)+node[ithEle.node[4]].y*xi(1)+node[ithEle.node[6]].y*xi(2)+node[ithEle.node[7]].y*xi(3));
       jz1=4.0*(node[ithEle.node[0]].z*(xi(0)-0.25)+node[ithEle.node[4]].z*xi(1)+node[ithEle.node[6]].z*xi(2)+node[ithEle.node[7]].z*xi(3));

       jx2=4.0*(node[ithEle.node[4]].x*xi(0)+node[ithEle.node[1]].x*(xi(1)-0.25)+node[ithEle.node[5]].x*xi(2)+node[ithEle.node[8]].x*xi(3));
       jy2=4.0*(node[ithEle.node[4]].y*xi(0)+node[ithEle.node[1]].y*(xi(1)-0.25)+node[ithEle.node[5]].y*xi(2)+node[ithEle.node[8]].y*xi(3));
       jz2=4.0*(node[ithEle.node[4]].z*xi(0)+node[ithEle.node[1]].z*(xi(1)-0.25)+node[ithEle.node[5]].z*xi(2)+node[ithEle.node[8]].z*xi(3));

       jx3=4.0*(node[ithEle.node[6]].x*xi(0)+node[ithEle.node[5]].x*xi(1)+node[ithEle.node[2]].x*(xi(2)-0.25)+node[ithEle.node[9]].x*xi(3));
       jy3=4.0*(node[ithEle.node[6]].y*xi(0)+node[ithEle.node[5]].y*xi(1)+node[ithEle.node[2]].y*(xi(2)-0.25)+node[ithEle.node[9]].y*xi(3));
       jz3=4.0*(node[ithEle.node[6]].z*xi(0)+node[ithEle.node[5]].z*xi(1)+node[ithEle.node[2]].z*(xi(2)-0.25)+node[ithEle.node[9]].z*xi(3));

       jx4=4.0*(node[ithEle.node[7]].x*xi(0)+node[ithEle.node[8]].x*xi(1)+node[ithEle.node[9]].x*xi(2)+node[ithEle.node[3]].x*(xi(3)-0.25));
       jy4=4.0*(node[ithEle.node[7]].y*xi(0)+node[ithEle.node[8]].y*xi(1)+node[ithEle.node[9]].y*xi(2)+node[ithEle.node[3]].y*(xi(3)-0.25));
       jz4=4.0*(node[ithEle.node[7]].z*xi(0)+node[ithEle.node[8]].z*xi(1)+node[ithEle.node[9]].z*xi(2)+node[ithEle.node[3]].z*(xi(3)-0.25));

       J.row(0)<<1,1,1,1;
       J.row(1)<<jx1,jx2,jx3,jx4;
       J.row(2)<<jy1,jy2,jy3,jy4;
       J.row(3)<<jz1,jz2,jz3,jz4;

       Jdet=J.determinant();
       Jinv=J.inverse();

       Iaug<<0,0,0,
             1,0,0,
             0,1,0,
             0,0,1;

       P=Jinv*Iaug;

       a1=P(0,0);
       a2=P(1,0);
       a3=P(2,0);
       a4=P(3,0);

       b1=P(0,1);
       b2=P(1,1);
       b3=P(2,1);
       b4=P(3,1);

       c1=P(0,2);
       c2=P(1,2);
       c3=P(2,2);
       c4=P(3,2);

       Nfx << (4.0*xi(0)-1)*a1, (4.0*xi(1)-1)*a2, (4.0*xi(2)-1)*a3, (4.0*xi(3)-1)*a4,
              4.0*(a1*xi(1)+a2*xi(0)), 4.0*(a2*xi(2)+a3*xi(1)), 4.0*(a1*xi(2)+a3*xi(0)),
              4.0*(a1*xi(3)+a4*xi(0)), 4.0*(a2*xi(3)+a4*xi(1)), 4.0*(a3*xi(3)+a4*xi(2));

       Nfy << (4.0*xi(0)-1)*b1, (4.0*xi(1)-1)*b2, (4.0*xi(2)-1)*b3, (4.0*xi(3)-1)*b4,
              4.0*(b1*xi(1)+b2*xi(0)), 4.0*(b2*xi(2)+b3*xi(1)), 4.0*(b1*xi(2)+b3*xi(0)),
              4.0*(b1*xi(3)+b4*xi(0)), 4.0*(b2*xi(3)+b4*xi(1)), 4.0*(b3*xi(3)+b4*xi(2));

       Nfz << (4.0*xi(0)-1)*c1, (4.0*xi(1)-1)*c2, (4.0*xi(2)-1)*c3, (4.0*xi(3)-1)*c4,
              4.0*(c1*xi(1)+c2*xi(0)), 4.0*(c2*xi(2)+c3*xi(1)), 4.0*(c1*xi(2)+c3*xi(0)),
              4.0*(c1*xi(3)+c4*xi(0)), 4.0*(c2*xi(3)+c4*xi(1)), 4.0*(c3*xi(3)+c4*xi(2));          

       B.row(0) << Nfx(0),0,0, Nfx(1),0,0, Nfx(2),0,0, Nfx(3),0,0, Nfx(4),0,0, Nfx(5),0,0, Nfx(6),0,0, Nfx(7),0,0, Nfx(8),0,0, Nfx(9),0,0;
       B.row(1) << 0,Nfy(0),0, 0,Nfy(1),0, 0,Nfy(2),0, 0,Nfy(3),0, 0,Nfy(4),0, 0,Nfy(5),0, 0,Nfy(6),0, 0,Nfy(7),0, 0,Nfy(8),0, 0,Nfy(9),0;
       B.row(2) << 0,0,Nfz(0), 0,0,Nfz(1), 0,0,Nfz(2), 0,0,Nfz(3), 0,0,Nfz(4), 0,0,Nfz(5), 0,0,Nfz(6), 0,0,Nfz(7), 0,0,Nfz(8), 0,0,Nfz(9);
       B.row(3) << Nfy(0),Nfx(0),0,  Nfy(1),Nfx(1),0, Nfy(2),Nfx(2),0, Nfy(3),Nfx(3),0, Nfy(4),Nfx(4),0, Nfy(5),Nfx(5),0, Nfy(6),Nfx(6),0,
                   Nfy(7),Nfx(7),0,Nfy(8),Nfx(8),0,  Nfy(9),Nfx(9),0;
       B.row(4) << 0,Nfz(0),Nfy(0), 0,Nfz(1),Nfy(1), 0,Nfz(2),Nfy(2), 0,Nfz(3),Nfy(3), 0,Nfz(4),Nfy(4), 0,Nfz(5),Nfy(5), 0,Nfz(6),Nfy(6),
                   0,Nfz(7),Nfy(7), 0,Nfz(8),Nfy(8), 0,Nfz(9),Nfy(9);
       B.row(5) << Nfz(0),0,Nfx(0), Nfz(1),0,Nfx(1), Nfz(2),0,Nfx(2), Nfz(3),0,Nfx(3), Nfz(4),0,Nfx(4), Nfz(5),0,Nfx(5), Nfz(6),0,Nfx(6),
                   Nfz(7),0,Nfx(7), Nfz(8),0,Nfx(8), Nfz(9),0,Nfx(9);

       Bt = B.transpose();
       K += weight*(Bt*E*B*Jdet/6.0);
       cout<<"Jdet="<<Jdet<<endl;
 }
   return;
}

问题的元素文件（elements.dat）是

    1
    0     1     2    3    4    5    6    7    8    9      
    24.0  -0.6295797707 -0.7709500359   -0.0962567111   0.0018562616
   -0.5930703716    0.2172077372    -0.7752988670   0.0051374961
    0.9488273288    -0.1721170560   0.2647686142    0.0037463536
    0.1956924749    0.2842456298    0.9385674601    0.0033773339

节点文件 (nodes.dat)

 2.0000000000   3.0000000000    4.0000000000    
 6.0000000000   3.0000000000    2.0000000000    
 2.0000000000   5.0000000000    1.0000000000    
 4.0000000000   3.0000000000    6.0000000000    
 4.0000000000   3.0000000000    3.0000000000    
 4.0000000000   4.0000000000    1.5000000000    
 2.0000000000   4.0000000000    2.5000000000    
 3.0000000000   3.0000000000    5.0000000000    
 5.0000000000   3.0000000000    4.0000000000    
 3.0000000000   4.0000000000    3.5000000000

图 10.8 中给出的正确刚度矩阵仅适用于前三个值

   KelemMat=[447, 324,72.....]

而我得到

   KelemMat=[ 70.617   134.209  -108.193 .... ]

还有一件事是在高斯点上得到一个负行列式是可以的……我确实得到了一个

1个回答

这是相同的matlab片段。

clc;
node = [2 3 4
        6 3 2 
        2 5 1
        4 3 6
        4 3 3
        4 4 1.5
        2 4 2.5
        3 3 5
        5 3 4
        3 4 3.5]; 
Y=480; %Youngs Modulus
nu=1/3; % Poisson ratio
alfa=(5.0+3.0*sqrt(5.0))/20; %Gauss Points
beta =(5.0-sqrt(5.0))/20;    %Gauss Points
weight=0.25; %weights for integration points
GaussPoints=[alfa beta beta beta;beta alfa beta beta;beta beta alfa  beta;beta beta beta alfa];
%Elasticity Matrix
E=[1-nu nu nu 0 0 0;nu 1-nu nu 0 0 0;nu nu 1-nu 0 0 0;0 0 0 0.5-nu 0 0;0 0 0 0 0.5-nu 0;0 0 0 0 0 0.5-nu];
E=Y/((1+nu)*(1-2*nu))*E;
%Stiffness Matrix
K=zeros(30,30);

for j=1:4
xi=GaussPoints(:,j);
jx1=4.0*(node(1,1)*(xi(1)-0.25)+node(5,1)*xi(2)+node(7,1)*xi(3)+node(8,1)*xi(4));
jy1=4.0*(node(1,2)*(xi(1)-0.25)+node(5,2)*xi(2)+node(7,2)*xi(3)+node(8,2)*xi(4));
jz1=4.0*(node(1,3)*(xi(1)-0.25)+node(5,3)*xi(2)+node(7,3)*xi(3)+node(8,3)*xi(4));

jx2=4.0*(node(5,1)*xi(1)+node(2,1)*(xi(2)-0.25)+node(6,1)*xi(3)+node(9,1)*xi(4));
jy2=4.0*(node(5,2)*xi(1)+node(2,2)*(xi(2)-0.25)+node(6,2)*xi(3)+node(9,2)*xi(4));
jz2=4.0*(node(5,3)*xi(1)+node(2,3)*(xi(2)-0.25)+node(6,3)*xi(3)+node(9,3)*xi(4));

jx3=4.0*(node(7,1)*xi(1)+node(6,1)*xi(2)+node(3,1)*(xi(3)-0.25)+node(10,1)*xi(4));
jy3=4.0*(node(7,2)*xi(1)+node(6,2)*xi(2)+node(3,2)*(xi(3)-0.25)+node(10,2)*xi(4));
jz3=4.0*(node(7,3)*xi(1)+node(6,3)*xi(2)+node(3,3)*(xi(3)-0.25)+node(10,3)*xi(4));

jx4=4.0*(node(8,1)*xi(1)+node(9,1)*xi(2)+node(10,1)*xi(3)+node(4,1)*(xi(4)-0.25));
jy4=4.0*(node(8,2)*xi(1)+node(9,2)*xi(2)+node(10,2)*xi(3)+node(4,2)*(xi(4)-0.25));
jz4=4.0*(node(8,3)*xi(1)+node(9,3)*xi(2)+node(10,3)*xi(3)+node(4,3)*(xi(4)-0.25));


J=[1 1 1 1;jx1 jx2 jx3 jx4;jy1 jy2 jy3 jy4;jz1 jz2 jz3 jz4];

Jdet=det(J);
Jinv=inv(J);
Iaug=[0 0 0;1 0 0;0 1 0;0 0 1];
P=Jinv*Iaug;

       a1=P(1,1);
       a2=P(2,1);
       a3=P(3,1);
       a4=P(4,1);

       b1=P(1,2);
       b2=P(2,2);
       b3=P(3,2);
       b4=P(4,2);

       c1=P(1,3);
       c2=P(2,3);
       c3=P(3,3);
       c4=P(4,3);

Nfx=[(4.0*xi(1)-1)*a1 (4.0*xi(2)-1)*a2  (4.0*xi(3)-1)*a3 (4.0*xi(4)-1)*a4 ...
   4.0*(a1*xi(2)+a2*xi(1)) 4.0*(a2*xi(3)+a3*xi(2)) 4.0*(a1*xi(3)+a3*xi(1)) ...
   4.0*(a1*xi(4)+a4*xi(1)) 4.0*(a2*xi(4)+a4*xi(2)) 4.0*(a3*xi(4)+a4*xi(3))];

Nfy=[(4.0*xi(1)-1)*b1 (4.0*xi(2)-1)*b2 (4.0*xi(3)-1)*b3 (4.0*xi(4)-1)*b4 ...
  4.0*(b1*xi(2)+b2*xi(1)) 4.0*(b2*xi(3)+b3*xi(2)) 4.0*(b1*xi(3)+b3*xi(1)) ...
  4.0*(b1*xi(4)+b4*xi(1)) 4.0*(b2*xi(4)+b4*xi(2)) 4.0*(b3*xi(4)+b4*xi(3))];

Nfz=[(4.0*xi(1)-1)*c1 (4.0*xi(2)-1)*c2 (4.0*xi(3)-1)*c3 (4.0*xi(4)-1)*c4 ...
  4.0*(c1*xi(2)+c2*xi(1)) 4.0*(c2*xi(3)+c3*xi(2)) 4.0*(c1*xi(3)+c3*xi(1)) ...
  4.0*(c1*xi(4)+c4*xi(1)) 4.0*(c2*xi(4)+c4*xi(2)) 4.0*(c3*xi(4)+c4*xi(3))];

B=[Nfx(1) 0 0 Nfx(2) 0 0 Nfx(3) 0 0 Nfx(4) 0 0 Nfx(5) 0 0 Nfx(6) 0 0 Nfx(7) 0 0 Nfx(8) 0 0 Nfx(9) 0 0 Nfx(10) 0 0;
  0 Nfy(1) 0 0 Nfy(2) 0 0 Nfy(3) 0 0 Nfy(4) 0 0 Nfy(5) 0 0 Nfy(6) 0 0 Nfy(7) 0 0 Nfy(8) 0 0 Nfy(9) 0 0 Nfy(10) 0;
  0 0 Nfz(1) 0 0 Nfz(2) 0 0 Nfz(3) 0 0 Nfz(4) 0 0 Nfz(5) 0 0 Nfz(6) 0 0 Nfz(7) 0 0 Nfz(8) 0 0 Nfz(9) 0 0 Nfz(10);
  Nfy(1) Nfx(1) 0 Nfy(2) Nfx(2) 0 Nfy(3) Nfx(3) 0 Nfy(4) Nfx(4) 0 Nfy(5) Nfx(5) 0 Nfy(6) Nfx(6) 0 Nfy(7) Nfx(7) 0 Nfy(8) Nfx(8) 0 Nfy(9) Nfx(9) 0 Nfy(10) Nfx(10) 0;
  0 Nfz(1) Nfy(1) 0 Nfz(2) Nfy(2) 0 Nfz(3) Nfy(3) 0 Nfz(4) Nfy(4) 0 Nfz(5) Nfy(5) 0 Nfz(6) Nfy(6) 0 Nfz(7) Nfy(7) 0 Nfz(8) Nfy(8) 0 Nfz(9) Nfy(9) 0 Nfz(10) Nfy(10);
  Nfz(1) 0 Nfx(1) Nfz(2) 0 Nfx(2) Nfz(3) 0 Nfx(3) Nfz(4) 0 Nfx(4) Nfz(5) 0 Nfx(5) Nfz(6) 0 Nfx(6) Nfz(7) 0 Nfx(7) Nfz(8) 0 Nfx(8) Nfz(9) 0 Nfx(9) Nfz(10) 0 Nfx(10)];  

K=K+weight*(B'*E*B*Jdet/6.0);
Jdet
end
其它你可能感兴趣的问题
上一篇如何使用 Ansys 设计建模器（或任何其他 3D CAD 软件）构建椭圆体 下一篇自动分类一系列 2D 轮廓或 3D 闭合曲面