cacuek.h

这是等参单元的有限元程序

 #ifndef __CACUEK_
 #define __CACUEK_
 #include "shape2d.h"
 #include "eld.h"
  void CacuEK(double E,double u,double tt,double *xo,double *yo,double *ek,
            int *Ele,int IK,int nGauss);
 //计算单元刚度矩阵
  void CacuEK(double E,double u,double tt,double *xo,double *yo,double *ek,
            int *Ele,int IK,int nGauss)
/*参考说明:
  E--杨式模量;
  u--泊松比;
  tt--单元厚度,仅用与平面应力的情况下;
  xo--长度为8的数组指针,存放单元节点的x坐标;
  yo--长度为8的数组指针,存放单元节点的y坐标;
  ek--长度为16*16的数组指针,存放按行排列的单元刚度矩阵;
  Ele--长度为8的数组指针,存放单元节点编号;
  IK--=0,平面应力;=1,平面应变;=2,轴对称;
  nGauss--nGauss!=2;采用3*3的高斯积分,否则采用2*2的高斯积分;
 */
  {
	 int i,j,m,n,ii,jj,ij,LL;
	 double detj,rx,N[8],Nx[8],Ny[8],B[4][16],D[4][4],AB[3],H[3];
	 if(IK==2)
	 LL=4;
	 else LL=3;                            //设定高斯积分点坐标及权系数
	 if(nGauss!=2)
	  {
	    AB[0]=-0.7745966692; AB[1]=0; AB[2]=-AB[0]; H[0]=5.0/9.0; H[1]=8.0/9.0; H[2]=H[0]; ij=1;
      }
	  else
	  {
	    AB[0]=-0.5773502962; AB[2]=-AB[0]; H[0]=1.0;H[2]=1.0;ij=2;
	  }
   eld(E,u,D,IK);                    //计算D矩阵
	  for(i=0;i<256;i++)
	      ek[i]=0.0;
	  for(i=0;i<4;i++)
	     for(j=0;j<16;j++)
		    B[i][j]=0.0;                  //B清零矩阵
	  for(ii=0;ii<3;ii=ii+ij)
	     for(jj=0;jj<3;jj=jj+ij)
    {
	  for(i=0;i<8;i++)
	    {
		  N[i]=0.0;Nx[i]=0.0;Ny[i]=0.0;   //清零
		}
		detj=Inv_jaco(AB[ii],AB[jj],xo,yo,Nx,Ny,Ele)*H[ii]*H[jj];
	if(IK==2)                           //轴对称时计算半径
	   {
	     FN(AB[ii],AB[jj],N,Ele);
		 rx=0.0;
		   for(i=0;i<8;i++)
		     if(Ele[i])
			   rx+=xo[i]*N[i];
	     detj=detj*2.0*3.1415926*rx;
		}
	  for(i=0;i<8;i++)
	   if(Ele[i])
	  {
	    B[0][2*i]=Nx[i];B[1][2*i+1]=Ny[i];B[2][2*i]=Ny[i];B[2][2*i+1]=Nx[i];
		  if(IK==2)
		    B[3][2*i]=N[i]/rx;
	  }
	  for(i=0;i<16;i++)
	    for(j=0;j<16;j++)
		  for(m=0;m<LL;m++)
		    for(n=0;n<LL;n++)
			  ek[16*i+j]+=B[m][i]*D[m][n]*B[n][j]*detj*tt;    //完成刚度矩阵的高斯积分
	  }
   }

#endif