ganjian.cpp

有限元编程问题

#include<stdio.h>
#include<math.h>


#define ne 3                                                   //单元数
#define nj 4                                                   //节点数
#define nz 6                                                   //支撑数
#define npj 0                                                  //节点载荷数
#define npf 1                                                  //非节点载荷数
#define nj3 12                                                 //节点位移总数
#define dd 6                                                   //半带宽
#define e0 2.1E8                                              //弹性模量
#define a0 0.008                                               //截面积
#define i0 1.22E-4                                             //单元惯性距 
#define pi 3.141592654                                


int jm[ne+1][3]={{0,0,0},{0,1,2},{0,2,3},{0,4,3}};             /*gghjghg*/             
double gc[ne+1]={0.0,1.0,2.0,1.0};                             
double gj[ne+1]={0.0,90.0,0.0,90.0};
double mj[ne+1]={0.0,a0,a0,a0};
double gx[ne+1]={0.0,i0,i0,i0};
int zc[nz+1]={0,1,2,3,10,11,12};
double pj[npj+1][3]={{0.0,0.0,0.0}};
double pf[npf+1][5]={{0,0,0,0,0},{0,-20,1.0,2.0,2.0}};
double kz[nj3+1][dd+1],p[nj3+1];
double pe[7],f[7],f0[7],t[7][7];
double ke[7][7],kd[7][7];


//**kz[][]—整体刚度矩阵
//**ke[][]—整体坐标下的单元刚度矩阵
//**kd[][]—局部坐标下的单位刚度矩阵
//**t[][]—坐标变换矩阵

//**这是函数声明
void jdugd(int);
void zb(int);
void gdnl(int);
void dugd(int);


//**主程序开始
void main()
{
 int i,j,k,e,dh,h,ii,jj,hz,al,bl,m,l,dl,zl,z,j0;
 double cl,wy[7];
 int im,in,jn;

//***********************************************
//<功能：形成矩阵P>
//***********************************************

 if(npj>0)
 {
	 for(i=1;i<=npj;i++)
	 {                                                       //把节点载荷送入P
		 j=pj[i][2];                         
         p[j]=pj[i][1];
	 }
 }
 if(npf>0)
 {
	 for(i=1;i<=npf;i++)
	 {                                                      //求固端反力F0
		 hz=i;
		 gdnl(hz);
		 e=(int)pf[hz][3];
		 zb(e);                                             //求单元号码
		 for(j=1;j<=6;j++)                                  //求坐标变换矩阵T
		 {
			 pe[j]=0.0;
			 for(k=1;k<=6;k++)                              //求等效节点载荷
			 { 
				 pe[j]=pe[j]-t[k][j]*f0[k];
			 }
		 }
		 al=jm[e][1];
		 bl=jm[e][2];
		 p[3*al-2]=p[3*al-2]+pe[1];                         //将等效节点载荷送到P中
		 p[3*al-1]=p[3*al-1]+pe[2];
		 p[3*al]=p[3*al]+pe[3];
		 p[3*bl-2]=p[3*bl-2]+pe[4];
		 p[3*bl-1]=p[3*bl-1]+pe[5];
		 p[3*bl]=p[3*bl]+pe[6];
	 }
 }


 //*************************************************
 //<功能：生成整体刚度矩阵kz[][]>
 for(e=1;e<=ne;e++)                                        //按单元循环
 { 
  	 dugd(e);                                              //求整体坐标系中的单元刚度矩阵ke
	 for(i=1;i<=2;i++)                                     //对行码循环
	 {
		 for(ii=1;ii<=3;ii++)
		 { 
			 h=3*(i-1)+ii;                                 //元素在ke中的行码
			 dh=3*(jm[e][i]-1)+ii;                         //该元素在KZ中的行码
			 for(j=1;j<=2;j++)
			 {
				 for(jj=1;jj<=3;jj++)                      //对列码循环
				 {
					 l=3*(j-1)+jj;                         //元素在ke中的列码
					 zl=3*(jm[e][j]-1)+jj;                 //该元素在KZ中的行码
					 dl=zl-dh+1;                           //该元素在KZ*中的行码
					 if(dl>0)
						 kz[dh][dl]=kz[dh][dl]+ke[h][l];   //刚度集成
				 }
			 }
		 }
	 }
 }

//**引入边界条件**
  for(i=1;i<=nz;i++)                                       //按支撑循环
 {
	 z=zc[i];                                              //支撑对应的位移数
	 kz[z][l]=1.0;                                         //第一列置1
	 for(j=2;j<=dd;j++)
	 {
		 kz[z][j]=0.0;                                     //行置0
	 }
	 if((z!=1))
	 {
		 if(z>dd)
			 j0=dd;
		 else if(z<=dd)
			 j0=z;                                         //列（45度斜线）置0
		 for(j=2;j<=j0;j++)
			 kz[z-j+1][j]=0.0;
	 }
	 p[z]=0.0;                                             //P置0
 }




 for(k=1;k<=nj3-1;k++)
 {
	 if(nj3>k+dd-1)                                        //求最大行码
		 im=k+dd-1;
	 else if(nj3<=k+dd-1)
		 im=nj3;
	 in=k+1;
	 for(i=in;i<=im;i++)
	 {
		 l=i-k+1;
		 cl=kz[k][l]/kz[k][1];                             //修改KZ
		 jn=dd-l+1;
		 for(j=1;j<=jn;j++)
		 {
			 m=j+i-k;
		     kz[i][j]=kz[i][j]-cl*kz[k][m];
		 }
         p[i]=p[i]-cl*p[k];                               //修改P
	 }
 }




 p[nj3]=p[nj3]/kz[nj3][1];                                //求最后一个位移分量
 for(i=nj3-1;i>=1;i--)
 {
	 if(dd>nj3-i+1)
		 j0=nj3-i+1;
	 else j0=dd;                                          //求最大列码j0
	 for(j=2;j<=j0;j++)
	 {
		 h=j+i-1;
		 p[i]=p[i]-kz[i][j]*p[h];
	 }
	 p[i]=p[i]/kz[i][1];                                  //求其他位移分量
 }
	 printf("\n");
	 printf("_____________________________________________________________\n");
	 printf("NJ            U                 V                CETA       \n");     //输出位移
	 for(i=1;i<=nj;i++)
	 {
		 printf(" %-5d %14.11f     %14.11f     %14.11f\n",i,p[3*i-2],p[3*i-1],p[3*i]);
	 }
	 printf("_____________________________________________________________\n");
	 //*根据E的值输出相应E单元的N,Q,M(A,B)的结果**
	 printf("E             N                  Q                 M     \n");
	 //*计算轴力N,剪力Q,弯矩M*
	 for(e=1;e<=ne;e++)                                  //按单元循环
	 {
		 jdugd(e);                                       //求局部单元刚度矩阵kd
		 zb(e);                                          //求坐标变换矩阵T
		 for(i=1;i<=2;i++)
		 {
			 for(ii=1;ii<=3;ii++)
			 {
				 h=3*(i-1)+ii;
				 dh=3*(jm[e][i]-1)+ii;                   //给出整体坐标下单元节点位移
				 wy[h]=p[dh];
			 }
		 }
		 for(i=1;i<=6;i++)
		 {
			 f[i]=0.0;
			 for(j=1;j<=6;j++)
			 {
				 for(k=1;k<=6;k++)                      //求由节点位移引起的单元节点力
				 {
					 f[i]=f[i]+kd[i][j]*t[j][k]*wy[k];
				 }
			 }
		 }
		 if(npf>0)
		 {
			 for(i=1;i<=npf;i++)                        //按非节点载荷数循环
				 if(pf[i][3]==e)                        //找到荷载所在的单元
				 {
					 hz=i;
					 gdnl(hz);                          //求固端反力
					 for(j=1;j<=6;j++)                  //将固端反力累加
					 {
						 f[j]=f[j]+f0[j];
					 }
				 }
		 }
		 printf("%-3d(A)  %9.5f           %9.5f         %9.5f\n",e,f[1],f[2],f[3]);         //输出单元A（i）端内力
		 printf("   (B)  %9.5f           %9.5f         %9.5f\n",f[4],f[5],f[6]);            //输出单元B（i）端内力
	 }
      return;
}
//**主程序结束**

//************************************************************
//<功能：将非节点载荷下的杆端力计算出来存入f0[]>
//************************************************************

void gdnl(int hz)
{
	int ind,e;
	double g,c,l0,d;


	g=pf[hz][1];                                       //载荷值
	c=pf[hz][2];                                       //载荷位置
	e=(int)pf[hz][3];                                  //作用单元
	ind=(int)pf[hz][4];                                //载荷类型
	l0=gc[e];                                          //杆长
	d=l0-c;
	if(ind==1)
	{
		f0[1]=0.0;
		f0[2]=-(g*c*(2-2*c*c/(l0*l0)+(c*c*c)/(l0*l0*l0)))/2;         //均布载荷的固端反力
		f0[3]=-(g*c*c)*(6-8*c/l0+3*c*c/(l0*l0))/12;
		f0[4]=0.0;
		f0[5]=-g*c-f0[2];
		f0[6]=(g*c*c*c)*(4-3*c/l0)/(12*l0);
	}
	else
	{
		if(ind==2)                                    //横向集中力的固端反力
		{
			f0[1]=0.0;
			f0[2]=(-(g*d*d)*(l0+2*c))/(l0*l0*l0);
			f0[3]=-(g*c*d*d)/(l0*l0);
			f0[4]=0.0;
			f0[5]=(-g*c*c*(l0+2*d))/(l0*l0*l0);
			f0[6]=(g*c*c*d)/(l0*l0);
		}
		else
		{
			f0[1]=-(g*d/l0);                         //纵向集中力的固端反力
			f0[2]=0.0;
			f0[3]=0.0;
			f0[4]=-g*c/l0;
			f0[5]=0.0;
			f0[6]=0.0;
		}
	}
}

//************************************************************
//<功能：构成坐标变换矩阵>
//************************************************************
void zb(int e)
{
	double ceta,co,si;
	int i,j;
	ceta=(gj[e]*pi)/180;                             //角度变弧度
	co=cos(ceta);
	si=sin(ceta);
	t[1][1]=co;                                      //计算T右上角元素
	t[1][2]=si;
	t[2][1]=-si;
	t[2][2]=co;
	t[3][3]=1.0;
	for(i=1;i<=3;i++)
	{
		for(j=1;j<=3;j++)                            //计算T的左下角元素
		{
			t[i+3][j+3]=t[i][j];
		}
	}
}



//*****************************************************
//<计算局部坐标下单元刚度矩阵kd[][]>
//*****************************************************
void jdugd(int e)
{   
	double A0,l0,j0;
	int i;
	int j;
    
    
    A0=mj[e];                                       //面积
	l0=gc[e];                                       //杆长
	j0=gx[e];                                       //惯性钜

    
	for(i=0;i<=6;i++)
		for(j=0;j<=6;j++)                           //kd清0
		kd[i][j]=0.0;

    kd[1][1]=e0*A0/l0;
	kd[2][2]=12*e0*j0/pow(l0,3);
	kd[3][2]=6*e0*j0/pow(l0,2);
	kd[3][3]=4*e0*j0/l0;
	kd[4][1]=-kd[1][1];
	kd[4][4]=kd[1][1];
	kd[5][2]=-kd[2][2];                             //计算kd左下角各元素
	kd[5][3]=-kd[3][2];
	kd[5][5]=kd[2][2];
	kd[6][2]=kd[3][2];
	kd[6][3]=2*e0*j0/l0;
	kd[6][5]=-kd[3][2];
	kd[6][6]=kd[3][3];

    for(i=1;i<=6;i++)
		for(j=1;j<=i;j++)                          //将kd左下角元素按对称原则送到右下角
		kd[j][i]=kd[i][j];
}


//**********************************************************
//<计算整体坐标下单元刚度矩阵ke[][]>
//**********************************************************
void dugd(int e)
{
	int i,k,j,m;
	jdugd(e);                                     //计算局部单元刚度矩阵kd
	zb(e);                                        //计算坐标变换矩阵T
	for(i=1;i<=6;i++)
	{
		for(j=1;j<=6;j++)
		{
			ke[i][j]=0.0;
			for(k=1;k<=6;k++)
			{
				for(m=1;m<=6;m++)
				{
					ke[i][j]=ke[i][j]+t[k][i]*kd[k][m]*t[m][j];     //计算刚度坐标内单元刚度矩阵ke
				}
			}
		}
	}
}


//**程序结束**