function.h

来自「计算弹性模量的小程序」· C头文件 代码 · 共 162 行

void FN(double A,double B,double *N,int *Ele);
void FNA(double A,double B,double *NA,int *Ele);
void FNB(double A,double B,double *NB,int *Ele);
double  Inv_jaco(double A,double B,double *xo,double *yo,
                            double *Nx, double *Ny,int *Ele);
void eld(double E,double u, double (*D)[4],int Iopt);
void CacuEk(double E, double u,double Thick,double *xo,double *yo,
            double *ek,int *Ele,int Iopt,int nGauss);

void FN(double A,double B,double *N,int *Ele)
/* A、B—局部坐标?和?；
    N—长度为8的数组的指针，存放八个插值基函数值；
    Ele—长度为8的数组的指针，存放插值区域的节点
              构成；如果某元素为零，则意味该节点不存在。*/
{ if(Ele[4])  N[4]=0.5*(1-A*A)*(1-B);else N[4]=0.0;
  if(Ele[5])  N[5]=0.5*(1+A)*(1-B*B);else N[5]=0.0;
  if(Ele[6])  N[6]=0.5*(1-A*A)*(1+B);else N[6]=0.0;
  if(Ele[7])  N[7]=0.5*(1-A)*(1-B*B);else N[7]=0.0;
  N[0]=0.25*(1-A)*(1-B)-0.5*N[7]-0.5*N[4];
  N[1]=0.25*(1+A)*(1-B)-0.5*N[4]-0.5*N[5];
  N[2]=0.25*(1+A)*(1+B)-0.5*N[5]-0.5*N[6];
  N[3]=0.25*(1-A)*(1+B)-0.5*N[6]-0.5*N[7];
  return;
}

void FNA(double A,double B,double *NA,int *Ele)
/* A、B—局部坐标?和?；
    NA—长度为8的数组的指针，存放插值基函数偏导数值；
    Ele—长度为8的数组的指针，存放插值区域的节点
              构成；如果某元素为零，则意味该节点不存在。*/
{ if(Ele[4])  NA[4]=-A*(1-B);   else NA[4]=0.0;
  if(Ele[5])  NA[5]=0.5*(1-B*B);else NA[5]=0.0;
  if(Ele[6])  NA[6]=-A*(1+B);   else NA[6]=0.0;
  if(Ele[7])  NA[7]=-0.5*(1-B*B); else NA[7]=0.0;
  NA[0]=-0.25*(1-B)-0.5*NA[7]-0.5*NA[4];
  NA[1]=0.25*(1-B)-0.5*NA[4]-0.5*NA[5];
  NA[2]=0.25*(1+B)-0.5*NA[5]-0.5*NA[6];
  NA[3]=-0.25*(1+B)-0.5*NA[6]-0.5*NA[7];
  return;
}
void FNB(double A,double B,double *NB,int *Ele)
/* A、B—局部坐标?和?；
    NB—长度为8的数组的指针，存放插值基函数偏导数值；
    Ele—长度为8的数组的指针，存放插值区域的节点
              构成；如果某元素为零，则意味该节点不存在。*/
{ if(Ele[4])  NB[4]=-0.5*(1-A*A);else NB[4]=0.0;
  if(Ele[5])  NB[5]=-B*(1+A);    else NB[5]=0.0;
  if(Ele[6])  NB[6]=0.5*(1-A*A); else NB[6]=0.0;
  if(Ele[7])  NB[7]=-B*(1-A);    else NB[7]=0.0;
  NB[0]=-0.25*(1-A)-0.5*NB[7]-0.5*NB[4];
  NB[1]=-0.25*(1+A)-0.5*NB[4]-0.5*NB[5];
  NB[2]=0.25*(1+A)-0.5*NB[5]-0.5*NB[6];
  NB[3]=0.25*(1-A)-0.5*NB[6]-0.5*NB[7];
  return;
}
double  Inv_jaco(double A,double B,double *xo,double *yo,
                            double *Nx, double *Ny,int *Ele)
/* A、B—局部坐标?和?；
    xo—长度为8的数组的指针，存放插值节点总体x坐标；
    yo—长度为8的数组的指针，存放插值节点总体y坐标；
    Nx—长度为8的数组的指针，存放插值基函数对x偏导数；
    Ny—长度为8的数组的指针，存放插值基函数对y偏导数；
    Ele—长度为8的数组的指针，存放插值区域的节点
              构成；如果某元素为零，则意味该节点不存在。
*/
{  int i,j;  double detj,temp,NA[8],NB[8],jaco[2][2];
   for(i=0;i<2;i++) for(j=0;j<2;j++) jaco[i][j]=0.0;
   FNA(A,B,NA,Ele);FNB(A,B,NB,Ele);
   for(i=0;i<8;i++) if(Ele[i])
   { jaco[0][0]+=NA[i]*xo[i];jaco[0][1]+=NA[i]*yo[i];
     jaco[1][0]+=NB[i]*xo[i];jaco[1][1]+=NB[i]*yo[i];
    }
  detj=jaco[0][0]*jaco[1][1]-jaco[0][1]*jaco[1][0];
  /*if(detj<1.0e-16)
  { ofstream Out("file.err"); Out<<"detj<=0\n";Out.close();exit(1); }*/
  temp=jaco[0][0]/detj;  jaco[0][0]=jaco[1][1]/detj;  jaco[1][1]=temp;
  jaco[0][1]=-jaco[0][1]/detj;  jaco[1][0]=-jaco[1][0]/detj;
  for(i=0;i<8;i++) if(Ele[i])
  { Nx[i]=jaco[0][0]*NA[i]+jaco[0][1]*NB[i];
    Ny[i]=jaco[1][0]*NA[i]+jaco[1][1]*NB[i];
  }
  return(detj);
}

void eld(double E,double u, double (*D)[4],int Iopt)
/* 参数说明：
   E—杨氏模量；
   u—泊松比；
   D—4*4矩阵的指针，存储本构矩阵
      Iopt=0,平面应力；=1,平面应变；=2,轴对称
*/
{int i,j; double t,a,b;
  for(i=0;i<4;i++) for(j=0;j<4;j++) D[i][j]=0.0;  
  if(Iopt)                                               //处理平面应变或轴对称情况
  { t=E*(1-u)/(1+u)/(1-2*u);a=t*u/(1-u); b=0.5*t*(1-2*u)/(1-u);
     D[0][0]=t; D[1][1]=t; D[2][2]=b; D[3][3]=t; 
     D[0][1]=a; D[1][0]=a; D[0][3]=a; D[3][0]=a; D[1][3]=a; D[3][1]=a;
   }
  else                                                        //处理平面应力情况
  { t=E/(1-u*u);  D[0][0]=t; D[1][1]=t; D[2][2]= 0.5*t*(1.0-u);
    D[0][1]=t*u; D[1][0]=t*u ;
  }
}

void CacuEk(double E, double u,double Thick,double *xo,double *yo,
            double *ek,int *Ele,int Iopt,int nGauss)
/* 参数说明：
   E—杨氏模量；
   u—泊松比；
   Thick—单元厚度，仅用于平面应力情况下；
   xo—长度为8的数组指针，存放单元节点的x坐标
   yo—长度为8的数组指针，存放单元节点的y坐标
   ek—长度为256的数组指针，存放按行排列的单元刚度矩阵
   Ele—长度为8的数组指针，存放单元节点编号
   Iopt—单元选项，Iopt=0，平面应力；Iopt=1，平面应变，Iopt=2，轴对称
   nGauss—积分点个数，nGauss=2，采用2*2的高斯积分，否则采用3*3的高斯积分。
*/
{ 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(Iopt==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.5773502692;AB[2]=-AB[0];
	  H[0]=1.0;H[2]=1.0;
	  ij=2;
  }
  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矩阵清零
  eld(E,u,D,Iopt);                                  //计算D矩阵
  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(Iopt==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(Iopt==2) B[3][2*i]=N[i]/rx;
		  }
	  }//计算B矩阵
      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*Thick;
  }                        //完成单元刚度矩阵的高斯积分
}