平面四节点等参单元有限元程序
#include
#include
#include
#include
#include
float **float_two_array_malloc(int m,int n) //实型二维动态数组分配
{
float **a;
int i,j;
a=(float **)malloc(m*sizeof(float *));
for (i=0;i
a[i]=(float *)malloc(n*sizeof(float));
for (j=0; j
a[i][j]=0;
}
}
return a;
}
float *float_one_array_malloc(int n) //实型一维动态数组分配
{
float *a;
int j;
a=(float *)malloc(n*sizeof(float));
for (j=0; j
a[j]=0;
}
return a;
}
int **int_two_array_malloc(int m,int n) //整型二维动态数组分配
{
int **a;
int i,j;
a=(int **)malloc(m*sizeof(int *));
for (i=0;i
a[i]=(int *)malloc(n*sizeof(int));
for (j=0; j
a[i][j]=0;
}
}
return a;
}
int *int_one_array_malloc(int n) //整型一维动态数组分配
{
int *a;
int j;
a=(int *)malloc(n*sizeof(int));
for (j=0; j
a[j]=0;
}
return a;
}
//全局变量定义
int np,ne,nr1,nr2,nd,nf,n,ndf,ld,nm,nu1,nu2;
float *x,*y,*p,**pp,u1,u2; //x,y,me,nrr,p为输入信息
float coords[2][4],det,fun[4],pn[2][4],xjac[2][2],rjac[2][2],dnx[2][4]; //求单刚定义的变量
int *LD,**me,*nrr,*nu,*IS,nn; //nn为变带宽一维组中元素个数
//IS表示单元节点位移在结构位移阵列中对应的序号
float **EK,**s,*A,**ns,**mat; //EK为总体坐标系下单元的刚度矩阵
//A存储结构刚度矩阵的一维变带宽数组
//s存放各单元应力中心点的应力
//ns存放各节点的应力
void input() //输入函数定义
{
int i,j;
FILE *fp1;
fp1=fopen("input.dat","r"); //要输入的内容放在input.dat中
if(fp1==NULL)
{
printf("cann't open the file !\n");
exit(1);
}
//控制信息读入
fscanf(fp1,"%d",&np);
fscanf(fp1,"%d",&ne);
fscanf(fp1,"%d",&nr1);
fscanf(fp1,"%d",&nr2);
fscanf(fp1,"%d",&nd);
fscanf(fp1,"%d",&nf);
fscanf(fp1,"%d",&ld);
fscanf(fp1,"%d",&nm);
fscanf(fp1,"%d",&nu1);
fscanf(fp1,"%f",&u1);
fscanf(fp1,"%d",&nu2);
fscanf(fp1,"%f",&u2);
n=nf*np; ndf=nf*nd;
// printf("%d\t%f\t%d\t%f\n",nu1,u1,nu2,u2);
//动态数组生成
x=float_one_array_malloc(np);
y=float_one_array_malloc(np);
nrr=int_one_array_malloc(nr1+nr2); //位移约束为零的节点，nr1--x方向约束，nr2---y方向约束
nu=int_one_array_malloc(nu1+nu2); //位移约束非零的节点
pp=float_two_array_malloc(ld,2); //非零载荷 （集中力）,pp[ld][0]对应的自由度序号，pp[ld][1]集中力大小
p=float_one_array_malloc(n); //载荷矩阵
me=int_two_array_malloc(nd,ne); //单元节点的总体编号
mat=float_two_array_malloc(4,nm);//单元材料数组，先定义为各向同性材料
//原始数据读入
for(i
=0;i
fscanf(fp1,"%f",&x[i]);
fscanf(fp1,"%f",&y[i]);
}
for(i=0;i
for(j=0;j
fscanf(fp1,"%d",&me[j][i]);
me[j][i]=me[j][i]-1; //abaqus模型节点编号从1开始，C语言中从0开始
}
}
for(i=0;i
fscanf(fp1,"%d",&nrr[i]);//x方向约束的节点号
nrr[i]=(nrr[i]-1)*2; //转化为自由度序号，abaqus模型节点编号从1开始，C语言中从0开始
}
for(i=0;i
fscanf(fp1,"%d",&nrr[nr1+i]);//y方向约束的节点号，abaqus模型节点编号从1开始，C语言中从0开始
nrr[nr1+i]=(nrr[nr1+i]-1)*2+1;
}
for(i=0;i
for(i=0;i
for(i=0;i
j=int(pp[i][0]);
p[j]=pp[i][1];
}
for(i=0;i
for(j=0;j<4;j++)
{
fscanf(fp1,"%f",&mat[j][i]);
}
mat[0][i]=mat[0][i]-1;mat[1][i]=mat[1][i]-1; //abaqus模型节点编号从1开始，C语言中从0开始
}
for(i=0;i
fscanf(fp1,"%d",&nu[i]);//x方向约束的节点号
nu[i]=(nu[i]-1)*2; //转化为自由度序号，abaqus模型节点编号从1开始，C语言中从0开始
}
for(i=0;i
fscanf(fp1,"%d",&nu[nu1+i]);//y方向约束的节点号，abaqus模型节点编号从1开始，C语言中从0开始
nu[nu1+i]=(nu[nu1+i]-1)*2+1;
}
fclose(fp1);
printf("读入信息完成......\n");
}
void output() //结果输出函数
{ int i,j;
FILE *fp2;
fp2=fopen("output.dat","w"); //运算结果放在output.dat中
if(fp2==NULL)
{
printf("cann't open the file !\n");
exit(1);
}
fprintf(fp2,"位移为：\n"); //输出位移
for(i=0;i
fprintf(fp2,"\n各单元应力：Sx,\tSy,\tSxy,\tS1,\tS2,\tSmises\n"); //输出单元应力
for(i=0;i
fprintf(fp2,"%d\t",i+1);
for(j=0;j<6;j++)
fprintf(fp2,"%.4f\t",s[j][i]);
fprintf(fp2,"\n");
}
fclose(fp2);
printf("写入output.dat完成......\n");
}
void tecplot() //结果生成tecplot文件
{ int i,j;
float bl=1.0; //位移显示比例
FILE *fp3;
fp3=fopen("tecplot.plt","w"); //运算结果放在tecplot.plt中
if(fp3==NULL)
{
printf("cann't open the file !\n");
exit(1);
}
fprintf(fp3," TITLE = ' test '\n");
fprintf(fp3," VARIABLES=x,y,Ux,Uy,Sx,Sy,Sxy,Smax,Smin,Smises\n");
fprintf(fp3," ZONE N=%d,E=%d,F=FEPOINT,ET=quadrilateral\n",np,ne);
for(j=0;j
fprintf(fp3,"%.4f\t%.4f\t%.4f\t%.4f\t",x[j]+bl*p[j*2],y[j]+bl*p[j*2+1],p[j*2],p[j*2+1]);
for(i=0;i<6;i++)
{
fprintf(fp3,"%.4f\t",ns[i][j]);
}
fprintf(fp3,"\n");
}
for(i=0;i
for(j=0;j
fprintf(fp3,"\n");
}
fclose(fp3);
printf("写入tecplot.plt完成......\n");
}
void fld() //生成LD数组，存放组对角线元素在A中的位置
{
int k,i,j,j1,ig,nw;
LD=int_one_array_malloc(n);
LD[0]=0;
for(k=0;k
ig=np;
for(i=0;i
for(j=0;j
if(me[j][i]==k)
{ for(j1=0;j1
}
}
for(i=0;i
j=nf*k+i;
nw=nf*(k-ig)+i+1;
if(j!=0) LD[j]=LD[j-1]+nw;
}
}
nn=LD[n-1]+1; //nn为变带宽一维组中元素个数,LD下标从零开始
printf("生成LD数组完成......\n");
}
void fis(int ie) //形成单元ie的节点自由度号IS数组
{
int i,j,ii;
//IS=int_one_array_malloc(ndf);
for(i=0;i
for(j=0;j
IS[ii]=nf*me[i][ie]+j;
}
}
}
void fdnx(int ie,float gs,float gr)//-------------计算形函数对整体坐标的导数dnx[2][4]
{
int m,u,i,v;
float g1,g2,xi[4],eta[4],rec;
xi[0]=-1;xi[1]=1;xi[2]=1;xi[3]=-1;
eta[0]=-1;eta[1]=-1;eta[2]=1;eta[3]=1;
//--------读入单元节点坐标
for (i=0;i<4;i++)
{
coords[0][i]=x[me[i][ie]];
coords[1][i]=y[me[i][ie]];
}
//---------计算形函数和|J|
for(m=0;m<4;m++)
{
g1=(1+xi[m]*gr);
g2=(1+eta[m]*gs);
fun[m]=0.25*g1*g2;
pn[0][m]=0.25*xi[m]*g2;
pn[1][m]=0.25*eta[m]*g1;
}
for(u=0;u<2;u++)
{
for(v=0;v<2;v++)
{
det=0.0;
for(m=0;m<4;m++)
{
det=det+pn[u][m]*coords[v][m];
}
xjac[u][v]=det;
}
}
det=xjac[0][0]*xjac[1][1]-xjac[0][1]*xjac[1][0];
//-------------雅格比矩阵求逆----------------------
rec=1.0/det;
rjac[0][0]=rec*xjac[1][1];
rjac[1][1]=rec*xjac[0][0];
rjac[0][1]=-rec*xjac[0][1];
rjac[1][0]=-rec*xjac[1][0];
//------------形函数对整体坐标的导数-----------
for(u=0;u<4;u++)
{
for(v=0;v<2;v++)
{
dnx[v][u]=0.0;
for(m=0;m<2;m++)
{
dnx[v][u]=dnx[v][u]+rjac[v][m]*pn[m][u];
}
}
}
}
void fek(int ie) //形成总体坐标系下四节点平面等参单元ie的刚度矩阵EK（8，8）
{
int i,j,ii,jj,i1,i2,j1,j2,k,l;
float D[3][3],E,MU,d1,d2,d3,dxx,dxy,dyx,dyy,gs,gsh,gr,grh,rstg[3],h[3],dnix,dniy,dnjx,dnjy;
rstg[0]=-0.77459666924;rstg[1]=0;rstg[2]=-rstg[0];
h[0]=0.555555555556;h[1]=0.888888888889;h[2]=h[0];
//-----形成D矩阵-------
for(i=0;i
if((ie<=mat[1][i])&&(ie>=mat[0][i]))
{
E=mat[2][i];
MU=mat[3][i];
}
}
D[0][0]=E/(1.0-MU*MU);
D[1][1]=D[0][0];
D[0][1]=D[0][0]*MU;
D[1][0]=D[0][0]*MU;
D[0][2]=0;
D[1][2]=0;
D[2][0]=0;
D[2][1]=0;
D[2][2]=E/(2*(1+MU));
d1=D[0][0];
d2=D[0][1];
d3=D[2][2];
//---------单元刚度矩阵清零--------------
for(i=0;i
}
//---------------------------------
for(i=0;i
//--------------------形成I1，I2，J1，J2---------
ii=2*i
;
i1=ii;
i2=ii+1;
for(j=0;j
jj=2*j;
j1=jj;
j2=jj+1;
//---------------dxx,dxy,dyx,dyy-----------
dxx=0;
dxy=0;
dyx=0;
dyy=0;
for(k=0;k<3;k++)
{
gs=rstg[k]; //高斯点的值
gsh=h[k]; //加权系数
for(l=0;l<3;l++)
{
gr=rstg[l];
grh=h[l];
fdnx(ie,gs,gr);
dnix=dnx[0][i];
dniy=dnx[1][i];
dnjx=dnx[0][j];
dnjy=dnx[1][j];
dxx=dxx+dnix*dnjx*det*gsh*grh;
dxy=dxy+dnix*dnjy*det*gsh*grh;
dyx=dyx+dniy*dnjx*det*gsh*grh;
dyy=dyy+dniy*dnjy*det*gsh*grh;
}
}
//--------------分块形成单元刚度矩阵----------
EK[i1][j1]=dxx*d1+dyy*d3;
EK[i2][j2]=dyy*d1+dxx*d3;
EK[i1][j2]=dxy*d2+dyx*d3;
EK[i2][j1]=dyx*d2+dxy*d3;
}
}
}
void fasum(float*A) //将单元刚度矩阵叠加到结构刚度矩阵中并存到A（nn）中
{ int i,j,ni,nj,ij;
for(i=0;i
for(j=0;j
if(nj<=ni)
{ij=LD[ni]-(ni-nj);
A[ij]=A[ij]+EK[i][j];
}
}
}
}
void fr(float*A)//约束处理子程序
{
int i,j,t;
for(i=0;i
j=nrr[i];
t=LD[j];
A[t]=(float)1.0E30;
p[j]=0.0;
}
for(i=0;i
j=nu[i];
t=LD[j];
A[t]=(float)1.0E30;
p[j]=(float)1.0E30*u1;
}
for(i=nu1;i
j=nu[i];
t=LD[j];
A[t]=(float)1.0E30;
p[j]=(float)1.0E30*u2;
}
printf("约束处理完成......\n");
}
void cholesky(float*A) //用cholesky方法求解线性方程组
{
int i, j, k,ij,mi;
float** z; //增广矩阵
z = float_two_array_malloc(n, n+1);
/* 将原始数据导入z中 */
for(i=0; i
for(j=0; j<=i; j++)
{
if(i==0) mi=0;
else mi=i-(LD[i]-LD[i-1])+1;//第i行第一个非零元素的列号
if(j
{ ij=LD[i]-(i-j);
z[i][j]=A[ij];
}
}
}
for(i=0;i
for(j=i+1;j
}
for(i=0;i
float temp;
/* 分解过程 */
for(k=0; k
temp = 0;
for(i=0; i
temp = temp + z[i][k] * z[i][k];
}
z[k][k] = (float)sqrt(z[k][k] - temp);
for(i=k+1; i
temp = 0;
for(j=0; j
temp = temp + z[j][k] * z[j][i];
}
z[k][i] = (z[k][i] - temp) / z[k][k];
}
}
/* 回代过程 */
z[n-1][n] = z[n-1][n] / z[n-1][n-1];
for(k=n-2; k>=0; k--)
{
temp = 0;
for(i=k+1; i
temp = temp + z[k][i] * z[i][n];
}
z[k][n] = (z[k][n] - temp) / z[k][k];
}
/*将位移阵列存入p中*/
for(i=0; i
/*释放z数组*/
for (i=0;i
delete[] z[i];
}
delete[] z;
printf("求解线性方程组完成......\n");
}
void stress(int ie)// 求单元中心点的应力
{ int i,j,ii,j1,j2;
float d1,d2,d3,D[3][3],B[3][8],r[8],sig[3],E,MU,ss,rr,bi,ci,h1,h2,sx,sy,sxy;
//-----形成D矩阵-------
for(i=0;i
if((ie<=mat[1][i])&&(ie>=mat[0][i]))
{
E=mat[2][i];
MU=mat[3][i];
}
}
D[0][0]=E/(1.0-MU*MU);
D[1][1]=D[0][0];
D[0][1]=D[0][0]*MU;
D[1][0]=D[0][0]*MU;
D[0][2]=0;
D[1][2]=0;
D[2][0]=0;
D[2][1]=0;
D[2][2]=E/(2*(1+MU));
d1=D[0][0];
d2=D[0][1];
d3=D[2][2];
ss=0.0;
rr=0.0;
fdnx(ie,ss,rr);
//------形成B矩阵--------
for(i=0;i<4;i++)
{
ii=2*i;
j1=ii;
j2=ii+1;
bi=dnx[0][i];
ci=dnx[1][i];
B[0][j1]=bi;
B[1][j1]=0.0;
B[2][j1]=ci;
B[0][j2]=0.0;
B[1][j2]=ci;
B[2][j2]=bi;
}
//-----找出该单元四个节点的八个位移，写到r[8]中----
for(i=0;i<4;i++)
{
r[i*2]=p[me[i][ie]*2];
r[i*2+1]=p[me[i][ie]*2+1];
}
//-------计算应变sig[3]--------
for(i=0;i<3;i++)
{ sig[i]=0.0;}
for(i=0;i<8;i++)
{
for(j=0;j<3;j++)
{
sig[j]=sig[j]+B[j][i]*r[i];
}
}
//-------计算三个应力分量，写到s[0][ie]——s[2][ie]中-----
s[0][ie]=d1*sig[0]+d2*sig[1];
s[1][ie]=d2*sig[0]+d1*sig[1];
s[2][ie]=d3*sig[2];
//---------计算主应力、mises应力------写到s[3][ie]——s[5][ie]中----
sx=s[0][ie];
sy=s[1][ie];
sxy=s[2][ie];
h1=sx+sy;
h2=sqrt((sx-sy)*(sx-sy)+0.25*sxy*sxy);
s[3][ie]=(h1+h2)/2;
s[4][ie]=(h1-h2)/2;
s[5][ie]=sqrt(((s[3][ie]-s[4][ie])*(s[3][ie]-s[4][ie])+s[3][ie]*s[3][ie]+s[4][ie]*s[4][ie])/2);
}
void stress_average()// 将单元中心点的应力平均到节点上
{
int i,j,k;
for(i=0;i
for(j=0;j<4;j++)
{
for(k=0;k<6;k++)
{
ns[k][me[j][i]]=ns[k][me[j][i]]+s[k][i];
}
ns[6][me[j][i]]=ns[6][me[j][i]]+1;
}
}
for(i=0;i
for(j=0;j<6;j++)
{
ns[j][i]=ns[j][i]/ns[6][i];
}
}
printf("求节点应力完成......\n");
}
void main()
{
int i,j,k;
float *A;
input();
fld();
A=float_one_array_malloc(nn);
for(i=0;i
s=float_two_array_malloc(6,ne);
ns=float_two_array_malloc(7,np);
IS=int_one_array_malloc(ndf);
for(k=0;k
fis(k);
fek(k);
fasum(A);
}
printf("生成总体刚度矩阵完成......\n");
fr(A);
cholesky(A);
for(j=0;j
stress(j);
}
printf("求单元应力完成......\n");
stress_average();
output();
tecplot();
//测试input函数代码
// printf("%d,%d,%d,%d,%d,%d,%d,%d\n",np,ne,nr1,nr2,nd,nf,ld,nm);
// for(i=0;i
i++)
// printf("%.1f ",x[i]);
// printf("\n");
// for(i=0;i
// printf("\n");
// for(i=0;i
// printf("\n");
// for(i=0;i
// printf("\n");
// for(i=0;i
// printf("\n");
// }
//测试fld函数代码
//fld();
//for(i=0;i
//printf("\n");
//printf("%d\n",nn);
//测试fis函数代码
//fis(0);
//for(i=0;i
//printf("\n");
//测试fek函数代码
//for(i=0;i
// printf("\n"); printf("\n");
// }
//显示数组A
//for(i=0;i
//printf("\n");
//显示位移ui,vi
//printf("各点位移为：\n");
//for(j=0;j
//显示各单元应力分量
// printf("\n\n各单元应力分量为：\n");
//for(i=0;i
// printf("%d\t",i+1);
// for(j=0;j<6;j++)
// printf("%.4f\t",s[j][i]);
// printf("\n");
//}
//显示各节点应力分量
// printf("\n\n各节点应力分量为：\n");
//for(i=0;i
// printf("%d\t",i+1);
// for(j=0;j<6;j++)
// printf("%.4f\t",ns[j][i]);
// printf("\n");
//}
}
一下是输入文件内容
Input.dat中的内容，即输入的原始数据：
9 4 3 1 4 2 3 1 3 1.0 0 0.0//依次为np,ne,nr1,nr2,nd,nf,ld,nm,nu1,u1,nu2,u2
-20 15 //各节点的x,y坐标
-20 0
-20 -15
0 15
0 0
0 -15
20 15
20 0
20 -15
1 2 5 4 //各单元节点的整体编号
2 3 6 5
4 5 8 7
5 6 9 8
1 2 3 //x方向约束为零的位移对应的节点编号
3 //y方向约束为零的位移对应的节点编号
12 10 //集中力载荷
14 40
16 10
1 4 200 0.3 //材料信息——起、止单元号，E，MU
7 8 9 //位移非零的节点号