有限单元平面节点计算TC上机编程Word文档下载推荐.docx
- 文档编号:5789492
- 上传时间:2023-05-05
- 格式:DOCX
- 页数:23
- 大小:94.50KB
有限单元平面节点计算TC上机编程Word文档下载推荐.docx
《有限单元平面节点计算TC上机编程Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《有限单元平面节点计算TC上机编程Word文档下载推荐.docx(23页珍藏版)》请在冰点文库上搜索。
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);
请输入结果输出文件名:
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));
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<
maxD;
i++){
AF[i]=0.0;
for(j=0;
j<
j++)AK[i][j]=0.0;
numMP;
i++){
fscanf(fpin,"
%d%f%f%f\n"
itmp,&
E[i],&
mu[i],&
h[i]);
numN;
fscanf(fpin,"
%d%f%f\n"
X[i],&
Y[i]);
numE;
i++){
%d%d%d%d%d\n"
eMP[i],&
I[i],&
J[i],&
M[i]);
%d%d%d%d\n"
itmp,eMP[i],I[i],J[i],M[i]);
}
numLoadN;
%d%d%f\n"
noNF[i],&
dirNF[i],&
valNF[i]);
numLoadE;
%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]);
}
numFix;
%d%f\n"
noFix[i],&
valFix[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:
itmp=2*(noNF[i]-1)+dirNF[i]-1;
AF[itmp]+=valNF[i];
floatlen,*FE;
FE=(float*)malloc(6*sizeof(float));
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;
6;
j++)fprintf(fpout,"
%8g"
FE[j]);
\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;
restrict_condition_AKF(AK,AF,maxD,noFix[i],valFix[i]);
equation_solve_gauss_main(AK,maxD,AF);
fprintf(fpout,"
----------结点位移----------\n"
结点号Displacement_XDisplacement_Y\n"
%8d%15g%15g\n"
i+1,AF[2*i],AF[2*i+1]);
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));
-----------单元应力-----------\n"
单元号Sigma_XSigma_YTua_XY\n"
xj=X[J[i]-1];
xm=X[M[i]-1];
yi=Y[I[i]-1];
me[0]=2*I[i]-1;
for(j=0;
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);
3;
j++){
sigma[j]=0.0;
for(m=0;
m<
m++)sigma[j]+=S[j][m]*u[m];
}
%8d"
i+1);
fprintf(fpout,"
%12g"
sigma[j]);
请打开输出文件查看结果\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));
voidfree1int(int*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);
n1;
i++)
a[i]=(char*)a[0]+size*n2*i;
return(float**)a;
voidfree2float(float**p)
free(p[0]);
int**alloc2int(intn1,intn2)
size=sizeof(int);
return(int**)a;
voidfree2int(int**p)
voidmatrix_mul(float**a,float**b,float**c,intn,intk,intm)
inti,j,l;
n;
m;
c[i][j]=0.0;
for(l=0;
l<
k;
l++){
c[i][j]+=a[i][l]*b[l][j];
voidmatrix_tran(float**a,float**b,intn,intm)
inti,j;
b[j][i]=a[i][j];
voidprint_m(FILE*fp,float**a,introw,intcol,char*str)
if(str==NULL)
fprintf(fp,"
矩阵元素:
else
fprintf(fp,"
%s\n"
str);
row;
col;
fprintf(fp,"
%10g"
a[i][j]);
intequation_solve_gauss_main(float**a,intn,float*b)
inti,j,k;
floattmp,res,fmax,lik;
for(k=0;
k<
n-1;
k++){
res=select_main_factor(a,b,n,k);
if(res==-1)returnres;
for(i=k+1;
lik=a[i][k]/a[k][k];
a[i][k]=0.0;
for(j=k+1;
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++)tmp=tmp+a[i][j]*b[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;
tmp=fabs(a[i][m]);
if(fmax<
tmp){
fmax=tmp;
is=i;
if(fmax<
=0.0){
Nosoluctionforthequation!
!
return-1;
if(is!
=m){
for(j=m;
swap(&
a[is][j],&
a[m][j]);
swap(&
b[is],&
b[m]);
voidplane_D(floatE,floatmu,float**D,intflag)
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)
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;
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)
{
floattmp,mu2,bi,bj,bm,ci,cj,cm,A;
bj=ym-
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限 单元 平面 节点 计算 TC 上机 编程