二维稳定传热.cpp

数值传热学中的二维稳定传热的例子

#include <iostream.h>
#include <string.h>
#include <ostream.h>
#include <fstream.h>
#define N 8                                  //x轴方向
#define M 8                                  //y轴方向
main()
{
	float GAM(float,float);              //导热系数随坐标的变化值 
	float Sp(float,float);               //不同坐标下源项的值
	float Sc(float,float);               //不同坐标下源项的值
	float max(float,float);
	int i,j,k;
	float X_CV[N],X_point[N+1],Y_CV[M],Y_point[M+1],aw[N][M],as[N][M],ae[N][M],an[N][M],ap[N][M],b[N][M];
	float dx,dy,lx=8,ly=8,u,e,esp=10;
	float midx[N+1],midy[M+1],T[N+1][M+1];
	dx=lx/N;                              //取相同的微元
	dy=ly/M;
	X_CV[0]=0;                            //控制容积法划分网格
	for(i=1;i<N;i++)
	{
		X_CV[i]=dx*i;
		X_point[i]=0.5*(X_CV[i-1]+X_CV[i]);         //结点划分
	}
	X_point[N]=N*dx;
        Y_CV[0]=0;
	for(i=1;i<M;i++)
	{
		Y_CV[i]=dx*i;
		Y_point[i]=0.5*(Y_CV[i-1]+Y_CV[i]);
	}
	Y_point[M]=M*dx;
	midx[0]=midx[N-1]=dx/2;
	for(i=1;i<(N-1);i++)
		midx[i]=dx;
	midy[0]=midy[M-1]=dy/2;
	for(i=1;i<(M-1);i++)
		midy[i]=dy;
	for(j=1;j<(M-1);j++)                   //求四个系数
	{
		for(i=0;i<N;i++)
		{
			aw[i][j]=dy*GAM(X_point[i],Y_point[j])/midx[i];
			ae[i+1][j]=dy*GAM(X_point[i+1],Y_point[j])/midx[i+1];
		}
	}
	for(i=1;i<(N-1);i++)
	{
		for(j=0;j<M;j++)
		{
			as[i][j]=dx*GAM(X_point[i],Y_point[j])/midy[j];
			an[i][j+1]=dx*GAM(X_point[i],Y_point[j+1])/midy[j+1];
		}
	}
	for(j=1;j<M;j++)                      //确定ap和b
	{
		for(i=1;i<N;i++)
		{
			ap[i][j]=aw[i-1][j]+ae[i+1][j]+as[i][j-1]+an[i][j+1]-dx*dy*Sp(X_point[i],Y_point[j]);
			b[i][j]=Sc(X_point[i],Y_point[j])*dx*dy;
		}
	}
	for(i=1;i<N;i++)                      //边界条件
	{
		T[i][0]=350;
		T[i][M]=480;
	}
	for(i=1;i<M;i++)
	{
		T[0][i]=450;
		T[N][i]=300;
	}
	for(j=1;j<M;j++)                     //赋迭代初值
	{
		for(i=1;i<N;i++)
			T[i][j]=400;
	}
	k=0;
	while(k<100&&esp>0.08)               //高斯迭代
	{
		esp=0.08;
		for(j=1;j<M;j++)
		{
			for(i=1;i<N;i++)
			{
				u=T[i][j];
				T[i][j]=(aw[i-1][j]*T[i-1][j]+ae[i+1][j]*T[i+1][j]+an[i][j+1]*T[i][j+1]+as[i][j-1]*T[i][j-1]+b[i][j])/ap[i][j];
				e=T[i][j]-u;
				esp=max(e,esp);
			}
		}
			k=k+1;
	}
	if(k==100)                           //屏幕输出
		cout<<"迭代法不收敛!"<<endl;
	else
	{
		T[0][0]=T[0][M]=T[N][0]=T[N][M]=0;
		cout<<"---------------------------计算结果为---------------------------"<<endl;
		ofstream ft("text1.txt");                                                             //结果输入进文件
		ft<<"---------------------------计算结果为---------------------------"<<endl;
		for(j=0;j<=M;j++)
		{
			for(i=0;i<=N;i++)
			{
				cout.width(10);
				cout<<T[i][j];
				ft.width(10);
				ft<<T[i][j];
			}
			cout<<endl;
			ft<<endl;
		}
		ft.close();
	}
}

float GAM(float x,float y)
{
	float u;
	u=0*x+0*y+1;
	return u;
}
float Sp(float x,float y)
{
	float u;
	u=0*x+0*y;
	return u;
}
float Sc(float x,float y)
{
	float u;
	u=0*x+0*y;
	return u;
}
float max(float x,float y)
{
	if(x>y)
		return x;
	else 
		return y;
}