有限单元平面节点计算TC上机编程.docx
- 文档编号:9987582
- 上传时间:2023-02-07
- 格式:DOCX
- 页数:23
- 大小:94.50KB
有限单元平面节点计算TC上机编程.docx
《有限单元平面节点计算TC上机编程.docx》由会员分享,可在线阅读,更多相关《有限单元平面节点计算TC上机编程.docx(23页珍藏版)》请在冰豆网上搜索。
有限单元平面节点计算TC上机编程
#include
#include
#include
#include
voidswap(float*a,float*b);
float*alloc1float(intn);
voidfree1float(float*p);
int*alloc1lint(intn);
voidfree1int(int*p);
float**alloc2float(intn1,intn2);
voidfree2float(float**p);
int**alloc2int(intn1,intn2);
voidfree2int(int**p);
voidmatrix_mul(float**a,float**b,float**c,intn,intk,intm);
voidmatrix_tran(float**a,float**b,intn,intm);
voidprint_m(FILE*fp,float**a,introw,intcol,char*str=NULL);
intselect_main_factor(float**a,float*b,intn,intm);
intequation_solve_gauss_main(float**a,intn,float*b);
voidplane_D(floatE,floatmu,float**D,intflag);
voidPlane3Node_B(floatxi,floatyi,floatxj,floatyj,floatxm,floatym,float**B);
voidPlane3Node_S(float**D,float**B,float**S);
voidPlane3Node_ke(floatxi,floatyi,floatxj,floatyj,floatxm,floatym,floatE0,floatmu0,floath,float**ke,intflag);
voidPlaneT3_load(intdir,intI,intJ,intM,intn1,floatfv1,intn2,floatfv2,floath,floatlen,float*fe);
voidMakeAK(float**ke,int*me,intn,float**AK);
voidMakeAF(float*fe,int*me,intn,float*AF);
voidrestrict_condition_AKF(float**AK,float*AF,floatn,intr,floatur);
voidmain(intargc,char*argv[])
{
inti,j,m,itmp;
intnumN;
intnumE;
intnumMP;
intnumFix;
intnumLoadN;
intnumLoadE;
intflag;
float*E,*mu,*h;
float*X,*Y;
int*eMP,*I,*J,*M;
int*noNF;
int*dirNF;
float*valNF;
int*noFix;
float*valFix;
int*noEF;
int*dirEF;
int**ndEF;
float**valEF;
floatxi,yi,xj,yj,xm,ym;
int*me;
float**ke,*fe;
float**AK,*AF;
intmaxD;
FILE*fpin,*fpout;
charfilein[100],fileout[100];
printf("请输入参数文件名:
***.txt\n");
scanf("%s",filein);
printf("请输入结果输出文件名:
***.txt\n");
scanf("%s",fileout);
fpin=fopen(filein,"r");
if(fpin==NULL){
printf("不能打开文件\n");
exit(0);
}
fpout=fopen(fileout,"w");
fscanf(fpin,"%d%d%d%d%d%d%d\n",&numN,&numE,&numMP,&numLoadN,&numLoadE,&numFix,&flag);
fprintf(fpout,"结点数:
%d单元数:
%d材料数:
%d结点荷载数:
%d单元荷载数:
%d位移约束数:
%d\n",numN,numE,numMP,numLoadN,numLoadE,numFix);
if(flag>0)fprintf(fpout,"平面应力问题\n");
elsefprintf(fpout,"平面应变问题\n");
maxD=2*numN;
E=(float*)malloc(numMP*sizeof(float));
mu=(float*)malloc(numMP*sizeof(float));
h=(float*)malloc(numMP*sizeof(float));
Y=(float*)malloc(numN*sizeof(float));
X=(float*)malloc(numN*sizeof(float));
Y=(float*)malloc(numN*sizeof(float));
eMP=(int*)malloc(numE*sizeof(int));
I=(int*)malloc(numE*sizeof(int));
J=(int*)malloc(numE*sizeof(int));
M=(int*)malloc(numE*sizeof(int));
noNF=(int*)malloc(numLoadN*sizeof(int));
dirNF=(int*)malloc(numLoadN*sizeof(int));
valNF=(float*)malloc(numLoadN*sizeof(float));
noFix=(int*)malloc(numFix*sizeof(int));
valFix=(float*)malloc(numFix*sizeof(float));
noEF=(int*)malloc(numLoadE*sizeof(int));
dirEF=(int*)malloc(numLoadE*sizeof(int));
ndEF=alloc2int(numLoadE,2);
valEF=alloc2float(numLoadE,2);
me=(int*)malloc(6*sizeof(int));
ke=alloc2float(6,6);
AK=alloc2float(maxD,maxD);
AF=(float*)malloc(maxD*sizeof(float));
for(i=0;i AF[i]=0.0; for(j=0;j } for(i=0;i fscanf(fpin,"%d%f%f%f\n",&itmp,&E[i],&mu[i],&h[i]); } for(i=0;i fscanf(fpin,"%d%f%f\n",&itmp,&X[i],&Y[i]); } for(i=0;i fscanf(fpin,"%d%d%d%d%d\n",&itmp,&eMP[i],&I[i],&J[i],&M[i]); fprintf(fpout,"%d%d%d%d\n",itmp,eMP[i],I[i],J[i],M[i]); } for(i=0;i fscanf(fpin,"%d%d%f\n",&noNF[i],&dirNF[i],&valNF[i]); } for(i=0;i fscanf(fpin,"%d%d%d%f%d%f\n",&noEF[i],&dirEF[i],&ndEF[i][0],&valEF[i][0],&ndEF[i][1],&valEF[i][1]); fprintf(fpout,"%d%d%d%f%d%f\n\n",noEF[i],dirEF[i],ndEF[i][0],valEF[i][0],ndEF[i][1],valEF[i][1]); } for(i=0;i fscanf(fpin,"%d%f\n",&noFix[i],&valFix[i]); } for(i=0;i xi=X[I[i]-1];yi=Y[I[i]-1]; xj=X[J[i]-1];yj=Y[J[i]-1]; xm=X[M[i]-1];ym=Y[M[i]-1]; Plane3Node_ke(xi,yi,xj,yj,xm,ym,E[eMP[i]-1],mu[eMP[i]-1],h[eMP[i]-1],ke,flag); print_m(fpout,ke,6,6,"ke: "); me[0]=2*I[i]-1;me[1]=2*I[i]; me[2]=2*J[i]-1;me[3]=2*J[i]; me[4]=2*M[i]-1;me[5]=2*M[i]; MakeAK(ke,me,6,AK); } print_m(fpout,AK,maxD,maxD,"AK: "); for(i=0;i itmp=2*(noNF[i]-1)+dirNF[i]-1; AF[itmp]+=valNF[i]; } floatlen,*FE; FE=(float*)malloc(6*sizeof(float)); for(i=0;i xi=X[ndEF[i][0]-1];yi=Y[ndEF[i][0]-1]; xj=X[ndEF[i][1]-1];yj=Y[ndEF[i][1]-1]; len=sqrt((xj-xi)*(xj-xi)+(yj-yi)*(yj-yi)); PlaneT3_load(dirEF[i],I[noEF[i]-1],J[noEF[i]-1],M[noEF[i]-1],ndEF[i][0],valEF[i][0],ndEF[i][1],valEF[i][1],h[eMP[noEF[i]-1]-1],len,FE); fprintf(fpout,"单元%d等效结点荷载: ",noEF[i]); for(j=0;j<6;j++)fprintf(fpout,"%8g",FE[j]); fprintf(fpout,"\n"); me[0]=2*I[noEF[i]-1]-1;me[1]=2*I[noEF[i]-1]; me[2]=2*J[noEF[i]-1]-1;me[3]=2*J[noEF[i]-1]; me[4]=2*M[noEF[i]-1]-1;me[5]=2*M[noEF[i]-1]; MakeAF(FE,me,6,AF); } for(i=0;i restrict_condition_AKF(AK,AF,maxD,noFix[i],valFix[i]); } equation_solve_gauss_main(AK,maxD,AF); fprintf(fpout,"----------结点位移----------\n"); fprintf(fpout,"结点号Displacement_XDisplacement_Y\n"); for(i=0;i fprintf(fpout,"%8d%15g%15g\n",i+1,AF[2*i],AF[2*i+1]); } fprintf(fpout,"\n"); float**D,**S,**B; D=alloc2float(3,3); B=alloc2float(3,6); S=alloc2float(3,6); float*sigma; float*u; u=(float*)malloc(6*sizeof(float)); sigma=(float*)malloc(6*sizeof(float)); fprintf(fpout,"-----------单元应力-----------\n"); fprintf(fpout,"单元号Sigma_XSigma_YTua_XY\n"); for(i=0;i xi=X[I[i]-1];xj=X[J[i]-1];xm=X[M[i]-1]; yi=Y[I[i]-1];yj=Y[J[i]-1];ym=Y[M[i]-1]; me[0]=2*I[i]-1;me[1]=2*I[i]; me[2]=2*J[i]-1;me[3]=2*J[i]; me[4]=2*M[i]-1;me[5]=2*M[i]; for(j=0;j<6;j++) u[j]=AF[me[j]-1]; plane_D(E[eMP[i]-1],mu[eMP[i]-1],D,flag); plane_D(E[eMP[i]-1],mu[eMP[i]-1],D,flag); Plane3Node_B(xi,yi,xj,yj,xm,ym,B); Plane3Node_S(D,B,S); for(j=0;j<3;j++){ sigma[j]=0.0; for(m=0;m<6;m++)sigma[j]+=S[j][m]*u[m]; } fprintf(fpout,"%8d",i+1); for(j=0;j<3;j++){ fprintf(fpout,"%12g",sigma[j]); } fprintf(fpout,"\n"); } fprintf(fpout,"\n"); printf("请打开输出文件查看结果\n"); } voidswap(float*a,float*b) { floattmp; tmp=*a; *a=*b; *b=tmp; } float*alloc1float(intn) { float*a; a=(float*)malloc(n*sizeof(float)); returna; } voidfree1float(float*p) { free(p); } int*alloc1lint(intn) { int*a; a=(int*)malloc(n*sizeof(int)); returna; } voidfree1int(int*p) { free(p); } float**alloc2float(intn1,intn2) { inti; intsize; void**a; if(n2<=0)n2=n1; size=sizeof(float); a=(void**)malloc(n1*sizeof(void*)); a[0]=(void*)malloc(n1*n2*size); for(i=0;i a[i]=(char*)a[0]+size*n2*i; return(float**)a; } voidfree2float(float**p) { free(p[0]); free(p); } int**alloc2int(intn1,intn2) { inti; intsize; void**a; if(n2<=0)n2=n1; size=sizeof(int); a=(void**)malloc(n1*sizeof(void*)); a[0]=(void*)malloc(n1*n2*size); for(i=0;i a[i]=(char*)a[0]+size*n2*i; return(int**)a; } voidfree2int(int**p) { free(p[0]); free(p); } voidmatrix_mul(float**a,float**b,float**c,intn,intk,intm) { inti,j,l; for(i=0;i for(j=0;j c[i][j]=0.0; for(l=0;l c[i][j]+=a[i][l]*b[l][j]; } } } } voidmatrix_tran(float**a,float**b,intn,intm) { inti,j; for(i=0;i for(j=0;j b[j][i]=a[i][j]; } } } voidprint_m(FILE*fp,float**a,introw,intcol,char*str) { inti,j; if(str==NULL) fprintf(fp,"矩阵元素: \n"); else fprintf(fp,"%s\n",str); for(i=0;i for(j=0;j fprintf(fp,"%10g",a[i][j]); } fprintf(fp,"\n"); } fprintf(fp,"\n"); } intequation_solve_gauss_main(float**a,intn,float*b) { inti,j,k; floattmp,res,fmax,lik; for(k=0;k res=select_main_factor(a,b,n,k); if(res==-1)returnres; for(i=k+1;i lik=a[i][k]/a[k][k]; a[i][k]=0.0; for(j=k+1;j a[i][j]=a[i][j]-lik*a[k][j]; } b[i]=b[i]-lik*b[k]; } } b[n-1]=b[n-1]/a[n-1][n-1]; for(i=n-2;i>=0;i--){ tmp=0.0; for(j=i+1;j b[i]=(b[i]-tmp)/a[i][i]; } return0; } intselect_main_factor(float**a,float*b,intn,intm) { floatfmax,tmp; intis,i,j; fmax=-1.0; for(i=m;i tmp=fabs(a[i][m]); if(fmax fmax=tmp; is=i; } } if(fmax<=0.0){ printf("Nosoluctionforthequation! ! ! \n"); return-1; } if(is! =m){ for(j=m;j swap(&a[is][j],&a[m][j]); } swap(&b[is],&b[m]); } } voidplane_D(floatE,floatmu,float**D,intflag) { inti,j; floatcoef; if(flag<=0){ E=E/(1-mu*mu); mu=mu/(1-mu); } coef=E/(1-mu*mu); D[0][0]=coef;D[0][1]=coef*mu;D[0][2]=0.0; D[1][0]=D[0][1];D[1][1]=coef;D[1][2]=0.0; D[2][0]=0.0;D[2][1]=0.0;D[2][2]=coef*(1-mu)/2.0; } voidPlane3Node_B(floatxi,floatyi,floatxj,floatyj,floatxm,floatym,float**B) { inti,j; floatai,aj,am,bi,bj,bm,ci,cj,cm; floatA2; ai=xj*ym-xm*yj;aj=xm*yi-xi*ym;am=xi*yj-xj*yi; bi=yj-ym;bj=ym-yi;bm=yi-yj; ci=-xj+xm;cj=-xm+xi;cm=-xi+xj; A2=ai+aj+am; for(i=0;i<3;i++){ for(j=0;j<6;j++)B[i][j]=0.0; } B[0][0]=bi/A2;B[0][2]=bj/A2;B[0][4]=bm/A2; B[1][1]=ci/A2;B[1][3]=cj/A2;B[1][5]=cm/A2; B[2][0]=ci/A2;B[2][2]=cj/A2;B[2][4]=cm/A2; B[2][1]=bi/A2;B[2][3]=bj/A2;B[2][5]=bm/A2; } voidPlane3Node_S(float**D,float**B,float**S) { matrix_mul(D,B,S,3,3,6); } voidPlane3Node_ke(floatxi,floatyi,floatxj,floatyj,floatxm,floatym,floatE,floatmu,floath,float**ke,intflag) { inti,j; floattmp,mu2,bi,bj,bm,ci,cj,cm,A; bi=yj-ym;bj=ym-
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限 单元 平面 节点 计算 TC 上机 编程