yt.c

平板焊接过程的二维温度场的数值模拟程序。我在课程设计中做的作业.采用显示差分法

//编写：杨通2008-6-13
//调试：杨通2008-6-13
//二维平板焊接过程FDM
//计算结果保存在temp.txt下
/*'调用c语言的执行文件进行计算，一边计算一边保存，（c语言显示执行进度，写中间文件process，vb读取。）然后退出c执行
'执行程序首先进行网格实际剖分，xy方向节点坐标数组》建立m*n的两个数组保存温度（+温度初始化）》每个节点的参数，或者用节点
'类型索引。》分类计算，最外层进行时间步长迭代+写process+写温度场。*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define x_nod 50
#define y_nod 50
#define dx  (dx_min+dx_inc*abs(inter-j))
main()
{
float temp1[x_nod][y_nod],temp2[x_nod][y_nod];//主场
float max_t,inc_t,max_inct,min_inct,t,dy,dx_min,dx_inc,dz,t_last=0;
float t_pro,t_dec;
float ro,cp,namenda;//物性参数
float q_in;//焊枪热流密度
float fac1,fac2;
double Q1=0,Q2=0;/////必须是double
int i,j,inter,wherein=1;//下标
float init_in,init_out;
float saokongbi;
FILE  *fp;
///////////////////////初始化
t=0;
fp=fopen("temp.txt","w");
t_pro=3;t_dec=0.3;
inter=floor(y_nod/2);
init_in=20;init_out=20;
saokongbi=0.3;
fac1=50;fac2=10000;q_in=2000;
ro=500;cp=20;namenda=40;
max_t=50;inc_t=0.05;max_inct=0.1;min_inct=0.01;dy=0.1;dx_min=0.05;dx_inc=0.005;dz=0.3;
for(i=0;i<x_nod;i++)  for(j=0;j<y_nod;j++)  {temp1[i][j]=init_in;temp2[i][j]=init_in;}
///////////////////////

while (t<=max_t)
{t=t+inc_t;
for (i=1;i<x_nod-1;i++)
for (j=1;j<y_nod-1;j++)
{
if (j==inter&&i==wherein)//焊点
{
Q1=0;
Q1+=inc_t*namenda*(temp1[i][j-1]+temp1[i][j+1]-2*temp1[i][j])*dy*dz/dx;
Q1+=inc_t*namenda*(temp1[i-1][j]+temp1[i+1][j]-2*temp1[i][j])*dx*dz/dy;
Q1+=inc_t*q_in*fac2*dx*dy;//上
temp2[i][j]=temp1[i][j]+2*Q1/dx/dy/dz/ro/cp;
}
else
{
Q1=0;
Q1=Q1+inc_t*namenda*(temp1[i][j-1]+temp1[i][j+1]-2*temp1[i][j])*dy*dz/dx;
Q1+=inc_t*namenda*(temp1[i-1][j]+temp1[i+1][j]-2*temp1[i][j])*dx*dz/dy;
Q1+=inc_t*(init_out-temp1[i][j])*fac1*dx*dy;//上
temp2[i][j]=temp1[i][j]+2*Q1/dx/dy/dz/ro/cp;
}
}//end for 2
for(j=1;j<y_nod;j++)
{
Q1=0;Q2=0;
Q1+=inc_t*namenda*(temp1[0][j-1]+temp1[0][j+1]-2*temp1[0][j])*dy*dz/dx;
Q1+=inc_t*namenda*(temp1[1][j]-temp1[0][j])*dx*dz/dy;
Q1+=inc_t*(init_out-temp1[0][j])*fac1*dx*dz;
Q1+=inc_t*(init_out-temp1[0][j])*fac1*dx*dy;//上
temp2[0][j]=temp1[0][j]+2*Q1/dx/dy/dz/ro/cp;/////////////////
Q2+=inc_t*namenda*(temp1[x_nod-1][j-1]+temp1[x_nod-1][j+1]-2*temp1[x_nod-1][j])*dz*dy/dx;
Q2+=inc_t*namenda*(temp1[x_nod-2][j]-temp1[x_nod-1][j])*dz*dx/dy;
Q2+=inc_t*(init_out-temp1[x_nod-1][j])*fac1*dx*dz;
Q2+=inc_t*(init_out-temp1[x_nod-1][j])*fac1*dx*dy;//上
temp2[x_nod-1][j]=temp1[x_nod-1][j]+2*Q2/dx/dy/dz/ro/cp;
}
for(i=1;i<x_nod;i++)
{
Q1=0;Q2=0;
j=0;
Q1+=inc_t*namenda*(temp1[i-1][0]+temp1[i+1][0]-2*temp1[i][0])*dx*dz/dy;
Q1+=inc_t*namenda*(temp1[i][1]-temp1[i][0])*dy*dz/dx;
Q1+=inc_t*(init_out-temp1[i][0])*fac1*dy*dz;
Q1+=inc_t*(init_out-temp1[i][0])*fac1*dx*dy;//上
temp2[i][0]=temp1[i][0]+2*Q1/dx/dy/dz/ro/cp;/////////////////
j=y_nod-1;
Q2+=inc_t*namenda*(temp1[i-1][y_nod-1]+temp1[i+1][y_nod-1]-2*temp1[i][y_nod-1])*dx*dz/dy;
Q2+=inc_t*namenda*(temp1[i][y_nod-2]-temp1[i][y_nod-1])*dy*dz/dx;
Q2+=inc_t*(init_out-temp1[i][y_nod-1])*fac1*dy*dz;
Q2+=inc_t*(init_out-temp1[i][y_nod-1])*fac1*dx*dy;//上
temp2[i][y_nod-1]=temp1[i][y_nod-1]+2*Q2/dx/dy/dz/ro/cp;/////////////////
}
//4point
Q1=0;
i=0;j=0;
Q1+=inc_t*namenda*(temp1[i][j+1]-temp1[i][j])*dy*dz/dx;
Q1+=inc_t*namenda*(temp1[i+1][j]-temp1[i][j])*dx*dz/dy;
Q1+=inc_t*(init_out-temp1[i][j])*fac1*dx*dz;
Q1+=inc_t*(init_out-temp1[i][j])*fac1*dy*dz;
Q1+=inc_t*(init_out-temp1[i][i])*fac1*dx*dy;//上
temp2[i][j]=temp1[i][j]+2*Q1/dx/dy/dz/ro/cp;
Q1=0;
i=x_nod-1;j=0;
Q1+=inc_t*namenda*(temp1[i-1][j]-temp1[i][j])*dx*dz/dy;
Q1+=inc_t*namenda*(temp1[i][j+1]-temp1[i][j])*dy*dz/dx;
Q1+=inc_t*(init_out-temp1[i][j])*fac1*dx*dz;
Q1+=inc_t*(init_out-temp1[i][j])*fac1*dy*dz;
Q1+=inc_t*(init_out-temp1[i][i])*fac1*dx*dy;//上
temp2[i][j]=temp1[i][j]+2*Q1/dx/dy/dz/ro/cp;
Q1=0;
i=0;j=y_nod-1;
Q1+=inc_t*namenda*(temp1[i][j-1]-temp1[i][j])*dy*dz/dx;
Q1+=inc_t*namenda*(temp1[i+1][j]-temp1[i][j])*dx*dz/dy;
Q1+=inc_t*(init_out-temp1[i][j])*fac1*dx*dz;
Q1+=inc_t*(init_out-temp1[i][j])*fac1*dy*dz;
Q1+=inc_t*(init_out-temp1[i][i])*fac1*dx*dy;//上
temp2[i][j]=temp1[i][j]+2*Q1/dx/dy/dz/ro/cp;
Q1=0;
i=x_nod-1;j=y_nod-1;
Q1+=inc_t*namenda*(temp1[i][j-1]-temp1[i][j])*dy*dz/dx;
Q1+=inc_t*namenda*(temp1[i-1][j]-temp1[i][j])*dx*dz/dy;
Q1+=inc_t*(init_out-temp1[i][j])*fac1*dx*dz;
Q1+=inc_t*(init_out-temp1[i][j])*fac1*dy*dz;
Q1+=inc_t*(init_out-temp1[i][i])*fac1*dx*dy;//上
temp2[i][j]=temp1[i][j]+2*Q1/dx/dy/dz/ro/cp;
//判断时间步

/*if((t-t_last)>0.5*max_t)
{Q1=0;Q2=0;
for(i=0;i<x_nod;i++)  for(j=0;j<y_nod;j++)  {Q1+=temp1[i][j]/x_nod;Q2+=temp2[i][j]/x_nod;}
if ((abs(Q1/y_nod-Q2/y_nod)>(0.5*Q1/y_nod))&&((t_dec*inc_t)>=min_inct)) inc_t=t_dec*inc_t;//时间自适应
if ((abs(Q1/y_nod-Q2/y_nod)<(0.001*Q1/y_nod))&&((t_pro*inc_t)<=max_inct)) inc_t=t_pro*inc_t;//时间自适应
t_last=t;
}*/

for(i=0;i<x_nod;i++)  for(j=0;j<y_nod;j++)  temp1[i][j]=temp2[i][j];

//打印
printf("%3.2f\n",t);
//fprintf(fp,"%3.2f\n",t);
for(i=0;i<x_nod;i++)  
{for(j=0;j<y_nod;j++)
{
//printf("%3.2f\t",temp1[i][j]);
fprintf(fp,"%3.2f\t",temp1[i][j]);
} 
printf("%3.2f\n",temp1[i][j]);
fprintf(fp,"\n",temp1[i][j]);
}
//计算位置更新
wherein=(int)(t/max_t/saokongbi*y_nod);
}//end while
fclose(fp);
}//end main