⭐ 欢迎来到虫虫下载站！ | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站
📤 上传资源

📄 jiekou.cpp

📁 c源程序及与ANSYS衔接(JIEKOU) 提供了一个接口程序及方法
💻 CPP
📖 第 1 页 / 共 2 页
字号:
🔍
12 

void Astk()											//组装总刚矩阵
{float tc,eo,po,th,c;
int tll,i,mm,idn,n,j;
eo=WK[0];po=WK[1];									//eo为弹性模量，po为泊松比
for(i=0;i<3;i++)
	{for(mm=0;mm<3;mm++)	D[i][mm]=0.0;}
if(PTYPE==2)										//平面应变问题对eo和po转化
	{eo=eo/(1.0-po*po);	po=po/(1.0-po);}
D[0][0]=eo/(1.0-po*po);  D[1][1]=D[0][0];	 D[0][1]=D[0][0]*po;
D[1][0]=D[0][1];		D[2][2]=eo/2/(1.+po); 		//计算弹性矩阵D[]
NJ2=NODE_SUM*2;
Oned();
idn=ID[NJ2-1];
if(!(stiff=(float *)malloc(idn*sizeof(float))))		//为总刚矩阵的存储开辟临时空间
{printf("out of memory.\n");
	exit(0); }
for(i=0;i<idn;i++)*(stiff+i)=0;
for(i=1;i<=ELE_SUM;i++)								//叠加形成总刚矩阵
	{Asek(i,3);
	for(mm=1;mm<=6;mm++)
		{tll=LL[mm-1];
		for(n=1;n<=6;n++)
			{tc=LL[n-1];
			if(tll>=tc)
				{j=ID[tll-1]-tll+tc;
				*(stiff+j-1)=*(stiff+j-1)+EK[mm-1][n-1];}
			}
		}
	}
}

void Astp()											//组装总载向量函数
{float pe[4],gm,th,b,eg,ae,rij,qi,qj,bb;
int i,k,j;
gm=WK[2];th=WK[3];
for(i=0;i<NJ2;i++)DISPLACE[i]=0.0;
if(load_num>0)										//将分布力转化到节点上
{for(k=1;k<=load_num;k++)
	{i=INFO_LOADLINE[k-1][0];
	j=INFO_LOADLINE[k-1][1];
	rij=sqrt((NODE_X[j-1]-NODE_X[i-1])*(NODE_X[j-1]-NODE_X[i-1])+(NODE_Y[j-1]-NODE_Y[i-1])*(NODE_Y[j-1]-NODE_Y[i-1]));
	qi=(2.0*INFO_LOADLINE[k-1][2]+INFO_LOADLINE[k-1][3])*rij/6.0;
	qj=(2.0*INFO_LOADLINE[k-1][3]+INFO_LOADLINE[k-1][2])*rij/6.0;
	bb=INFO_LOADLINE[k-1][4];
	pe[0]=qi*cos(bb*PI/180);pe[1]=qi*sin(bb*PI/180);
	pe[2]=qj*cos(bb*PI/180);pe[3]=qj*sin(bb*PI/180);
	DISPLACE[2*i-1-1]=DISPLACE[2*i-1-1]+pe[0];
	DISPLACE[2*i-1]=DISPLACE[2*i-1]+pe[1];
	DISPLACE[2*j-1-1]=DISPLACE[2*j-1-1]+pe[2];
	DISPLACE[2*j-1]=DISPLACE[2*j-1]+pe[3];}
}
if(LOAD_NODE>0)										//将集中力转化到节点上
{for(i=1;i<=LOAD_NODE;i++)
	{j=INFO_LOADNODE[i-1][0];b=INFO_LOADNODE[i-1][2];
	DISPLACE[2*j-1-1]=DISPLACE[2*j-1-1]+INFO_LOADNODE[i-1][1]*cos(b*PI/180);
	DISPLACE[2*j-1]=DISPLACE[2*j-1]+INFO_LOADNODE[i-1][1]*sin(b*PI/180);}
}
if(gm>0)											//将重力(体积力)转化到结点上
{for(i=1;i<=ELE_SUM;i++)
		{ae=Asek(i,1);
		eg=-gm*ae*th/3.0;
		for(k=0;k<=2;k++)
			{j=LL[2*k];
			DISPLACE[j-1]=DISPLACE[j-1]+eg;}		
		}
}
}

void Cons()											//约束处理函数
{int i,j,it,itd;
for(i=1;i<=bc_num;i++)
	{for(j=2;j<=3;j++)
		{if(BC_INFO[i-1][j-1]==1)
			{it=2*BC_INFO[i-1][0]+j-3;
			itd=ID[it-1];
			*(stiff+itd-1)=*(stiff+itd-1)*1.0e+20;}	//有约束的主对角元乘大数
		}
	}
}

void Result()										//结果输出函数
{float w[3],w1[6],st[500][3],th,pyll,ryll,s1,s2,thet,stx,sty,stu;
int i,l,j,k1;
fp=fopen("node_displace.dat","w");						//打开结果输出文件
fprintf(fp,"\n");								//将单元原始信息和节点位移输出
fprintf(fp,"---PLANE PROBLEM ANALYSIS BY TRIANGULAR ELEMENT---\n");
fprintf(fp,"NUMBER OF NODES=%-9d NUMBER OF ELEMENTS=%-3d",NODE_SUM,ELE_SUM);
if(PTYPE==1)fprintf(fp,"\nPLANE STRESS PROBLEM");
else if(PTYPE==2)fprintf(fp,"\nPLANE STRAIN PROMLEM");
fprintf(fp,"PTYPE=%d",PTYPE);
fprintf(fp,"\nELASTICITY MODULUS=%-10.2f POISSON MODULUS=%-8.4f",WK[0],WK[1]);
fprintf(fp,"\nSPECIFIC GRAVITY=%-8.4f PLATE THICKNESS=%-8.4f",WK[2],WK[3]);
fprintf(fp,"\n\nNODAL COORDINATES&DISPLAYMENTS:");
fprintf(fp,"\nNODE:	 X	    Y	        U	         V");
for(i=1;i<=NODE_SUM;i++)
	{fprintf(fp,"\n%3d",i);
	fprintf(fp,"%8.3f     %8.3f     ",NODE_X[i-1] ,NODE_Y[i-1]);
	fprintf(fp,"%8.3e     %8.3e     ",DISPLACE[2*i-1-1] ,DISPLACE[2*i-1]);
}
fclose(fp);										//关闭结果输出文件
fp=fopen("ele_stress.dat","w");						//打开单元应力数据文件
fprintf(fp,"单元 	  x应力	  	y应力       方向\n");
for(l=1;l<=ELE_SUM;l++)
	{Asek(l,2);
	for(i=1;i<=6;i++)w1[i-1]=DISPLACE[LL[i-1]-1];	//计算应力
	for(i=0;i<3;i++)
		{w[i]=0.0;
		for(k1=0;k1<6;k1++)w[i]=w[i]+S[i][k1]*w1[k1];
		}
	st[l-1][0]=w[0];
	st[l-1][1]=w[1];
	st[l-1][2]=w[2];
	fprintf(fp,"%3d    %f     %f      %f\n",l,w[0],w[1],w[2]);}	//输出单元应力信息
fclose(fp);								//关闭单元应力数据文件
fp1=fopen("FOR_POST.DAT","w");
fprintf(fp1,"%9.4f%9.4f\n",(float)NODE_SUM,(float)ELE_SUM);
fp=fopen("node_stress.dat","w");				//打开节点应力数据文件
fprintf(fp,"节点       x应力         y应力      剪应力\n");
for(l=1;l<=NODE_SUM;l++)
	{k1=0;stx=0.0;sty=0.0;stu=0.0;
	for(j=1;j<=ELE_SUM;j++)
		{if((ELE_NUM[j-1][0]==l)||(ELE_NUM[j-1][1]==l)||(ELE_NUM[j-1][2]==l))
			{stx=stx+st[j-1][0];
			sty=sty+st[j-1][1];
			stu=stu+st[j-1][2];
			k1++;}
		}
	stx=stx/k1;sty=sty/k1;stu=stu/k1;
	fprintf(fp,"%3d   %6.4e    %6.4e   %6f\n",l,stx,sty,stu);	//输出节点应力
	fprintf(fp1,"%9.4f%9.4f%9.4f%9.4f",NODE_X[l-1],NODE_Y[l-1],DISPLACE[2*l-1-1],DISPLACE[2*l-1]);
	fprintf(fp1,"%9.4f%9.4f%9.4f\n",stx,sty,stu);}
fclose(fp);							//关闭节点应力数据文件
for(i=1;i<=ELE_SUM;i++)
	{fprintf(fp1,"%9.4f%9.4f%9.4f%9.4f",(float)ELE_NUM[i-1][0],(float)ELE_NUM[i-1][1],(float)ELE_NUM[i-1][2],(float)ELE_NUM[i-1][2]);
     fprintf(fp1,"%9.4f%9.4f%9.4f\n",st[i-1][0],st[i-1][1],st[i-1][2]);}
fclose(fp1);
fp=fopen("USER_POST.log","w");			//创立与ANSYS后处理的接口文件
fprintf(fp,"/PREP7\n");
fprintf(fp,"ET,1,PLANE42\n");
fprintf(fp,"DOF,rotx\n");
fprintf(fp,"DOF,roty\n");
fprintf(fp,"DOF,rotz\n");			//定义节点自由度
fprintf(fp,"*dim,info,,2\n");
fprintf(fp,"*vread,info(1),FOR_POST,dat\n");
fprintf(fp,"(2f9.4)\n");
fprintf(fp,"*dim,nd_info,,info(1),7\n");
fprintf(fp,"*vread,nd_info(1,1),FOR_POST,dat,,JIK,7,info(1),,1\n");
fprintf(fp,"(7f9.4)\n");
fprintf(fp,"*dim,ele_info,,info(2),7\n");
fprintf(fp,"*vread,ele_info(1,1),FOR_POST,dat,,JIK,7,info(2),,info(1)+1\n");
fprintf(fp,"(7f9.4)\n");
fprintf(fp,"*do,i,1,info(1)\n");
fprintf(fp,"N,i,nd_info(i,1),nd_info(i,2)\n");
fprintf(fp,"*enddo\n");
fprintf(fp,"*do,i,1,info(2)\n");
fprintf(fp,"E,ele_info(i,1),ele_info(i,2),ele_info(i,3),ele_info(i,4)\n");
fprintf(fp,"*enddo\n");
fprintf(fp,"/post1\n");
fprintf(fp,"*do,i,1,info(1)\n");
fprintf(fp,"DNSOL,i,u,x,nd_info(i,3)\n");
fprintf(fp,"DNSOL,i,u,y,nd_info(i,4)\n");
fprintf(fp,"DNSOL,i,rot,x,nd_info(i,5)\n");
fprintf(fp,"DNSOL,i,rot,y,nd_info(i,6)\n");
fprintf(fp,"DNSOL,i,rot,z,nd_info(i,7)\n");
fprintf(fp,"*enddo\n");				//把数组存放的结果赋给节点
fprintf(fp,"*do,i,1,info(2)\n");
fprintf(fp,"*do,j,1,3\n");
fprintf(fp,"num=ele_info(i,j)\n");
fprintf(fp,"desol,i,num,s,x,ele_info(i,5)\n");
fprintf(fp,"desol,i,num,s,y,ele_info(i,6)\n");
fprintf(fp,"desol,i,num,s,z,ele_info(i,7)\n");
fprintf(fp,"*enddo\n");
fprintf(fp,"*enddo");				//把数组存放的结果赋给单元
fclose(fp);}						//关闭接口文件

void Oned()							//一维存储函数
{int npoe[300][8],mini[300],kk[300],i,j,l,m,n,kki,kll,imi,k;
for(i=0;i<300;i++)
	{kk[i]=0;mini[i]=0;}
for(i=0;i<8;i++)
	{for(j=0;j<300;j++)npoe[j][i]=0;}
for(i=1;i<=ELE_SUM;i++)
	{for(j=1;j<=3;j++)
		{l=ELE_NUM[i-1][j-1];
		kk[l-1]=kk[l-1]+1;			//节点相关单元数
		n=kk[l-1];
		npoe[l-1][n-1]=i;}
	}
for(i=1;i<=NODE_SUM;i++)
	{mini[i-1]=i;
	kki=kk[i-1];
	for(j=1;j<=kki;j++)
		{k=npoe[i-1][j-1];			//i节点相关单元码
		for(l=1;l<=3;l++)
			{kll=ELE_NUM[k-1][l-1];
			if(kll<mini[i-1]) mini[i-1]=kll;}	//确定i行最小节点码
		}
	imi=i-mini[i-1];				//节点半带宽（不包括主对角子块）
	for(m=1;m<=2;m++)
		{n=2*(i-1)+m;
		ID[n-1]=2*imi+m;}			//每行半带宽（包括主对角元）
	}
for(i=2;i<=NJ2;i++)
ID[i-1]=ID[i-1]+ID[i-2];}			//主对角元在一维存储中的位置

void Solve()						//求解线性方程组函数
{int i,j,mi,mj,mij,ik,ii,kk,jj,ij,jk,im1,jm1,k,io,jo,l;
for(i=1;i<=NJ2;i++)
	{io=ID[i-1]-i;
	if(i!=1)
		{mi=ID[i-1-1]-io+1;
		for(j=mi;j<=i;j++)
			{jo=ID[j-1]-j;
			mj=1;
			if(j>1)mj=ID[j-1-1]-jo+1;
			mij=mi;
			if(mj>mi)mij=mj;
			ij=io+j;
			jm1=j-1;
			for(k=mij;k<=jm1;k++)
				{if(mij<=jm1)
					{ik=io+k;kk=ID[k-1];jk=jo+k;
					*(stiff+ij-1)=*(stiff+ij-1)-*(stiff+ik-1)*(*(stiff+kk-1))*(*(stiff+jk-1));}
				}
			if(j==i)continue;
			jj=ID[j-1];
			*(stiff+ij-1)=(*(stiff+ij-1))/(*(stiff+jj-1));
			DISPLACE[i-1]=DISPLACE[i-1]-(*(stiff+ij-1))*(*(stiff+jj-1))*DISPLACE[j-1];}
		}
	ii=io+i;
	DISPLACE[i-1]=DISPLACE[i-1]/(*(stiff+ii-1));}
for(l=2;l<=NJ2;l++)
	{i=NJ2-l+2;
	io=ID[i-1]-i;
	mi=ID[i-1-1]-io+1;
	im1=i-1;
	for(j=mi;j<=im1;j++)
		{if(mi<=im1)
			{ij=io+j;
			DISPLACE[j-1]=DISPLACE[j-1]-(*(stiff+ij-1))*DISPLACE[i-1];}
		}
	}
for(i=0;i<NJ2;i++)
	{if(fabs(DISPLACE[i])<1.0e-9) DISPLACE[i]=0.0;}
free(stiff);}					//求解完毕，释放为总刚矩阵开辟的临时空间
12

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -