jiekou.cpp

c源程序及与ANSYS衔接(JIEKOU) 提供了一个接口程序及方法

#define ANSYS
#include"string.h"
#include"stdio.h"
#include"conio.h"
#include"stdlib.h"
#include"process.h"
#include"malloc.h"
#include"math.h"
#define PI 3.1415926

float NODE_X[500],NODE_Y[500],INFO_LOADLINE[30][5],INFO_LOADNODE[30][3],WK[4],*stiff;
int ELE_NUM[1000][3],BC_INFO[50][3],NODE_SUM,ELE_SUM,BC_LINE,BC_NODE,LOAD_LINE,LOAD_NODE,PTYPE;
int bc_num,load_num,NJ2,SH[8],LL[6],ID[1000];
float DISPLACE[1000],D[3][3],EK[6][6],S[3][6];
FILE *fp,*fp1;

char Input(char *);							//读入相关参数
float Asek(int n,int iask);					//求解单元刚度矩阵
void Astk(),Astp(),Cons(),Oned(),Solve();
void Result();								//输出相关结果

void main()
{char instr[30];							//存放输入和输出的文件名
int i;
printf("\nplease input the datafile name:");
scanf("%s",instr);							//输入原始数据文件名
i=Input(instr);								//读入数据，返回信息变量
if(i==1){printf("Data file input OK!Any key continue\n");getch();}	
if(i==2){printf("Data file open error!\n");exit(0);}		
Astk();										//组装总刚矩阵
Astp();										//组装总载矩阵
Cons();										//约束处理
Solve();									//求解方程组
printf("\nCalculating finish.Anykey to continue...\n");
getch();
Result();									//给出结果并输出
printf("\nAnykey to QUIT!\n");
getch();
}

char Input(char *str)
{int i,j,nd,k1,k2,j1,n,sum,ss=0,sp=0;			//sp标识当前已读入的信息数据个数个数
float x1,y1,x2,y2,nr,dx1,dy1,dk,dx,dyx,q1,q2,aL;
float a1,a2,a3,a4,w2[20];
float *w1;										//指向读入的信息数据的临时存储空间
if((fp=fopen(str,"r+"))==NULL)return 2;			//打开数据文件
if(ferror(fp))return 2;
fscanf(fp,"%d %d %d %d %d %d %d",&PTYPE,&NODE_SUM,&ELE_SUM,&BC_LINE,&BC_NODE,&LOAD_LINE,&LOAD_NODE);
//读入控制数据，PTYPE为1表示平面应力问题，为2表示平面应变问题
#ifdef ANSYS									//条件编译，若定义了“ANSYS”则执行下面的语句
sum=(3*NODE_SUM)+(14*ELE_SUM)+(5*BC_LINE)+(3*BC_NODE)+(7*LOAD_LINE)+(3*LOAD_NODE)+4;
//计算需要读入的信息数据个数
if(!(w1=(float *)malloc(sum*sizeof(float))))		//开辟临时空间给信息数据，首地址为w1
	{printf("out of memory.\n");	exit(0); }		//若空间不够则报警
if((fp1=fopen("node_ansys.in","r+"))==NULL)return 2;
for(ss=0;ss<3*NODE_SUM;ss++)
fscanf(fp1,"%f",w1+ss);								//从w1地址开始存入节点信息数据
fclose(fp1);
if((fp1=fopen("element_ansys.in","r+"))==NULL)return 2;
for(ss=3*NODE_SUM;ss<(3*NODE_SUM+14*ELE_SUM);ss++)
fscanf(fp1,"%f",w1+ss);								//继续存入单元信息数据
fclose(fp1);
for(ss=(3*NODE_SUM+14*ELE_SUM);ss<sum;ss++)
fscanf(fp,"%f",w1+ss);								//继续存入其他信息数据
#else												//条件编译，若没有定义“ANSYS”则执行下面的语句
sum=(3*NODE_SUM)+(4*ELE_SUM)+(5*BC_LINE)+(3*BC_NODE)+(7*LOAD_LINE)+(3*LOAD_NODE)+4;
if(!(w1=(float *)malloc(sum*sizeof(float))))
	{printf("out of memory.\n");	exit(0); }
if((fp1=fopen("node.in","r+"))==NULL)return 2;
for(ss=0;ss<3*NODE_SUM;ss++)
fscanf(fp1,"%f",w1+ss);
fclose(fp1);
if((fp1=fopen("element.in","r+"))==NULL)return 2;
for(ss=3*NODE_SUM;ss<(3*NODE_SUM+4*ELE_SUM);ss++)
fscanf(fp1,"%f",w1+ss);
fclose(fp1);
for(ss=(3*NODE_SUM+4*ELE_SUM);ss<sum;ss++)
fscanf(fp,"%f",w1+ss);
#endif											//条件编译结束
fclose(fp);
if(NODE_SUM>0)									//NODE_SUM表示节点总数
{for(i=1;i<=NODE_SUM;i++)
	{nd=*(w1+sp+3*i-2-1);						//nd读入节点编号
	NODE_X[nd-1]=*(w1+sp+3*i-1-1);				//NODE_X[nd-1]存储该节点横坐标
	NODE_Y[nd-1]=*(w1+sp+3*i-0-1);	}			//NODE_Y[nd-1]存储该节点纵坐标
sp=sp+3*NODE_SUM;}
if(ELE_SUM>0)									//ELE_SUM表示单元总数
#ifdef ANSYS									//条件编译，若定义了“ANSYS”则执行下面的语句
{for(i=1;i<=ELE_SUM;i++)
	{nd=*(w1+sp+14*i-1);						//nd读入单元编号
	ELE_NUM[nd-1][0]=*(w1+sp+14*i-13-1);
	ELE_NUM[nd-1][1]=*(w1+sp+14*i-12-1);
	ELE_NUM[nd-1][2]=*(w1+sp+14*i-11-1);	}	//ELE_NUM[nd-1]存储该单元三节点编码
sp=sp+14*ELE_SUM;}
#else										//条件编译，若没有定义“ANSYS”则执行下面的语句
{for(i=1;i<=ELE_SUM;i++)
	{nd=*(w1+sp+4*i-3-1);					//nd读入单元编号
	ELE_NUM[nd-1][0]=*(w1+sp+4*i-2-1);
	ELE_NUM[nd-1][1]=*(w1+sp+4*i-1-1);
	ELE_NUM[nd-1][2]=*(w1+sp+4*i-1); }		//ELE_NUM[nd-1]存储该单元三节点编码
sp=sp+4*ELE_SUM;}
#endif										//条件编译结束
bc_num=0;
if(BC_LINE>0)								//BC_LINE表示有约束的直线边数
{for(i=1;i<=BC_LINE;i++)
	{x1=*(w1+sp+5*i-5);y1=*(w1+sp+5*i-4);
	x2=*(w1+sp+5*i-3);y2=*(w1+sp+5*i-2);	//存储该边的起终点坐标
	nr=*(w1+sp+5*i-1);
	//nr读入该边的约束信息，-1则x方向有约束，0则x,y方向均有约束，1则y方向有约束
	k1=0;k2=0;
	if(nr<0.1) k1=1;						//x方向有约束则令k1为1
	if(nr>-0.1) k2=1;						//y方向有约束则令k2为1
	dx=x2-x1;
	//下面分斜率是否存在两种情况分别寻找在直线上的节点，并把边的约束信息转化到节点上
	if(fabs(dx)<0.0001)					//斜率不存在
		{for(j=1;j<=NODE_SUM;j++)
			{dx1=NODE_X[j-1]-x1;
			if(fabs(dx1)<0.0001)			//顺次判断各个节点是否在该直线上
				{bc_num++;					//存储数组元素加1
				BC_INFO[bc_num-1][0]=j;
				BC_INFO[bc_num-1][1]=k1;
				BC_INFO[bc_num-1][2]=k2; }
	//BC_INFO[]存储在该直线上的节点的约束信息，依次为节点编号和x,y方向的约束情况
			}
		}
	else							//斜率存在
		{dk=(y2-y1)/dx;					//求斜率
		for(j=1;j<=NODE_SUM;j++)
			{dy1=NODE_Y[j-1]-y1-dk*(NODE_X[j-1]-x1);
			if(fabs(dy1)<0.0001)			//顺次判断各个节点是否在该直线上
				{bc_num++;
				BC_INFO[bc_num-1][0]=j;
				BC_INFO[bc_num-1][1]=k1;
				BC_INFO[bc_num-1][2]=k2;}
			}
		}
	}
sp=sp+5*BC_LINE;}
if(BC_NODE>0)							//BC_NODE表示有约束的节点数
{for(i=1;i<=BC_NODE;i++)
	{bc_num++;
	BC_INFO[bc_num-1][0]=*(w1+sp+3*i-2-1);
	BC_INFO[bc_num-1][1]=*(w1+sp+3*i-1-1);
	BC_INFO[bc_num-1][2]=*(w1+sp+3*i-1);}
//BC_INFO[]存储节点约束信息，依次为节点编号和x,y方向的约束情况
	sp=sp+3*BC_NODE;}
load_num=0;
if(LOAD_LINE>0)							//LOAD_LINE表示有分布载荷的直线边数
{for(i=1;i<=LOAD_LINE;i++)
	{j1=0;
	x1=*(w1+sp+7*i-6-1);y1=*(w1+sp+7*i-5-1);
	x2=*(w1+sp+7*i-4-1);y2=*(w1+sp+7*i-3-1);		//读入该边的约束始终点坐标
	q1=*(w1+sp+7*i-2-1);q2=*(w1+sp+7*i-1-1);		//读入始终点的载荷集度
	aL=*(w1+sp+7*i-1);			//读入从x轴按右手法则转至该载荷外作用线的角度
	dx=x2-x1;		
	if(fabs(dx)<0.0001)						//斜率不存在
		{for(j=1;j<=NODE_SUM;j++)
			{dx1=NODE_X[j-1]-x1;
			if(fabs(dx1)<0.0001)			//顺次判断各节点是否在该直线上
				{j1++;		
				w2[j1-1]=j;	}				//w2[]存储在该直线上的节点编号
			}
		}
	else									//斜率存在
		{dk=(y2-y1)/dx;						//求斜率
		for(j=1;j<=NODE_SUM;j++)
			{dy1=NODE_Y[j-1]-y1-dk*(NODE_X[j-1]-x1);
			if(fabs(dy1)<0.0001)		//顺次判断各节点是否在该直线上
				{j1++;					//j1最后得到在该直线上的节点数
				w2[j1-1]=j;}			//w2[]存储在该直线上的节点编号
			}
		}
	a1=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
	a2=(q2-q1)/a1;					//求出载荷梯度
	//下面一个循环将对直线边的分布载荷转化到各个单元边上
	for(n=1;n<=j1-1;n++)
		{load_num++;
		k1=w2[n-1];k2=w2[n];				//单元边始终节点编码
		a3=sqrt((NODE_X[k1-1]-x1)*(NODE_X[k1-1]-x1)+(NODE_Y[k1-1]-y1)*(NODE_Y[k1-1]-y1));
		a4=sqrt((NODE_X[k2-1]-x1)*(NODE_X[k2-1]-x1)+(NODE_Y[k2-1]-y1)*(NODE_Y[k2-1]-y1));
		//a3,a4分别保存始终点到约束起点的距离
		INFO_LOADLINE[load_num-1][0]=k1; 
		INFO_LOADLINE[load_num-1][1]=k2;				//第load_num条单元边的始终节点编码
		INFO_LOADLINE[load_num-1][2]=q1+a3*a2;
		INFO_LOADLINE[load_num-1][3]=q1+a4*a2;			//第load_num条单元边的始终载荷
		INFO_LOADLINE[load_num-1][4]=aL;}				//第load_num条单元边的载荷方向
	}
sp=sp+7*LOAD_LINE;}
if(LOAD_NODE>0)											//LOAD_NODE表示集中力作用节点数
{for(i=1;i<=LOAD_NODE;i++)
	{INFO_LOADNODE[i-1][0]=*(w1+sp+3*i-2-1);
	INFO_LOADNODE[i-1][1]=*(w1+sp+3*i-1-1);
	INFO_LOADNODE[i-1][2]=*(w1+sp+3*i-1);	}
	//INFO_LOADNODE[]存储集中力作用点信息，依次为节点坐标，力的大小，作用方向
	sp=sp+3*LOAD_NODE;}
for(i=0;i<4;i++)WK[i]=*(w1+sp+i);		//WK[]依次读入材料的弹性模量,泊松比,容重和板厚
free(w1);								//释放存储信息数据的临时空间
return 1;}

float Asek(int n,int iask)				//计算单元刚度矩阵，n为单元编号
{float bi,ci,cm,bm,cj,bj,b[3][6],bb[6][3];
int i,j,m,k1=0,k2=0,k3,k4;
float th,ae,a2,at;
float x1,x2,x3,y1,y2,y3;
i=ELE_NUM[n-1][0];j=ELE_NUM[n-1][1];m=ELE_NUM[n-1][2];		//节点单元码i,j,m
cm=NODE_X[j-1]-NODE_X[i-1];bm=NODE_Y[i-1]-NODE_Y[j-1];
cj=NODE_X[i-1]-NODE_X[m-1];bj=NODE_Y[m-1]-NODE_Y[i-1];
ae=(bj*cm-bm*cj)/2.0;										//ae得到单元面积
th=WK[3];													//th得到板厚
if(iask>1)
	{for(k1=0;k1<3;k1++)
		{for(k2=0;k2<6;k2++)
			{b[k1][k2]=0.0;bb[k2][k1]=0.0;}
		}
	b[0][0]=-bj-bm;	b[2][1]=b[0][0];  b[0][2]=bj; 	 b[2][3]=bj;
	b[0][4]=bm;  	b[2][5]=bm;	 b[1][1]=-cj-cm;  b[2][0]=b[1][1];
	b[1][3]=cj;		b[2][2]=cj;		 b[1][5]=cm;	 b[2][4]=cm;
	a2=0.5/ae;
	for(k3=0;k3<=2;k3++)
		{for(k4=0;k4<=5;k4++)
		b[k3][k4]=a2*b[k3][k4];}							//生成几何矩阵b[3][6]
	for(k3=0;k3<=2;k3++)
		{for(k4=0;k4<=5;k4++)
			{S[k3][k4]=0.0;
			for(k1=0;k1<=2;k1++)
			S[k3][k4]=S[k3][k4]+D[k3][k1]*b[k1][k4];}
		}
	}														//生成应力矩阵s[3][6]
if(iask>2)
	{for(k1=0;k1<3;k1++)
		{for(k2=0;k2<6;k2++)
		bb[k2][k1]=b[k1][k2];}			//生成几何矩阵的转置矩阵bb[k2][k1]=b[k1][k2]
	for(k3=0;k3<6;k3++)
		{for(k4=0;k4<6;k4++)
			{EK[k3][k4]=0.0;
			for(k1=0;k1<3;k1++)
			EK[k3][k4]=EK[k3][k4]+bb[k3][k1]*S[k1][k4];}
		}
	at=ae*th;
	for(k1=0;k1<6;k1++)
		{for(k4=0;k4<6;k4++)
		EK[k1][k4]=at*EK[k1][k4];	}
	}									//生成单元刚度矩阵EK[6][6]
LL[0]=i+i-1;LL[1]=i+i;LL[2]=j+j-1;LL[3]=j+j;LL[4]=m+m-1;LL[5]=m+m; 	//三节点总自由度码
return ae;															//返回单元面积
}