text2.cpp

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

//假设一根长为L的圆台杆，小端半径为r，大端半径为R。小端温度为T0，大端温度为TL，源项为线型分布，即S=Sc+SpT.
#include <iostream.h>
#include <string.h>
#include <math.h>
#include <ostream.h>
#include <fstream.h>
#define N 6                                                       //边界个数
#define PI 3.14159
#define L 10                                                       //杆的长度
#define R 2
#define r 1
#define M 5                                                        //t最大值除以平均时间变量dt,要求能够整除
#define E 100
#define dt 2
#define rou 0.01
#define C 1
void main()
{
	float CVarea(float);                                           //求解界面面积的函数
	float CVface_x[N];                                             //界面的坐标值
	float point_x[N+1];                                            //结点的坐标值
 	float Gam[N];                                                  //界面处的扩散导热系数
	float Coe_a[N-1];                                                //计算系数a
	float Coe_ap[N];
	float T[N+1];                                                  //所求时刻结点处的温度值
	float T00[N+1]={5,8.06029,13.4066,17.9386,21.5924,24.2598,25};                                                //上一时刻的结点温度值
	float b[N-1],Sc,Sp,se,sw,ap0,error;                              //各个参数
	float u[N],l[N],y[N];
	float T1[M]={25,25,25,25,25},T2[M]={5,5,5,5,5};
	float tt;
	int i,j,num;
	//划分控制容积边界
	cout<<"Please input CV_face Coordinate values:"<<endl;
	for(i=0;i<N;i++)
		cin>>CVface_x[i];                                          //输入边界处的坐标值
	//划分结点
	point_x[0]=0;                                                  //结点初始值为0
	for(i=0;i<N;i++)
		point_x[i+1]=CVface_x[i]+(CVface_x[i+1]-CVface_x[i])/2;
	point_x[N]=L;                                                  //结点末点坐标值为L
	//输入每个结点处的物体的扩散传导系数
	cout<<"输入每个结点处的物体的扩散传导系数:"<<endl;
	for(i=0;i<N;i++)
		cin>>Gam[i];
	num=0;
	for(j=1;j<=M;)
	{
	//计算系数a及计算b
		Sc=1;
		Sp=0;
		Coe_ap[0]=0;
		ap0=rou*C;
		for(i=0;i<N;i++)
			{
			se=2*CVarea(CVface_x[i]);
			sw=2*CVarea(CVface_x[i+1]);
			b[i]=(se*sw)/(se+sw)*(CVface_x[i+1]-CVface_x[i])*(Sc+ap0*T00[i+1]);             //调和平均值
			Coe_a[i]=(Gam[i]*se/2)/(point_x[i+1]-point_x[i]);
			}
		for(i=1;i<N;i++)
			{
			se=2*CVarea(CVface_x[i-1]);
			sw=2*CVarea(CVface_x[i]);
			Coe_ap[i]=(Coe_a[i-1]+Coe_a[i])+(ap0-Sp)*(se*sw)/(se+sw)*(CVface_x[i]-CVface_x[i-1]);
			}
	//追赶法解方程组
		b[0]=b[0]+Coe_a[0]*T1[j-1];
		b[N-2]=b[N-2]+Coe_a[N-1]*T2[j-1];
		u[0]=-Coe_ap[1];
		for(i=1;i<(N-1);i++)
			{
			l[i]=Coe_a[i]/u[i-1];
			u[i]=-Coe_ap[i+1]-l[i]*Coe_a[i];
			}
		y[0]=-b[0];
		for(i=1;i<(N-1);i++)
			y[i]=-b[i]-l[i]*y[i-1];
		T[N-1]=y[N-2]/u[N-2];
		for(i=N-3;i>=0;i--)
			T[i+1]=(y[i]-Coe_a[i+1]*T[i+2])/u[i];
		T[0]=T1[j-1];
		T[N]=T2[j-1];
		error=0;
		for(i=0;i<=N;i++)
			error=error+(T00[i]-T[i])*(T00[i]-T[i]);
		error=sqrt(error);
		if(error<=E)
		{
			num=0;
			for(i=0;i<=N;i++)
				T00[i]=T[i];
			tt=j*dt;
			cout<<"The results is:"<<tt<<endl;
			for(i=0;i<=N;i++)
				cout<<T[i]<<endl;
			ofstream ft("text2.txt");
			ft<<"结点坐标    结点温度"<<endl;
			for(i=0;i<=N;i++)
				{
				ft<<point_x[i]<<"           ";
				ft<<T[i]<<endl;
				}
			ft.close();
			j++;
		}
		else
		{
			for(i=0;i<=N;i++)
				T00[i]=T[i];
			num++;
			if(num>10)
			{
				cout<<"错误!"<<endl;
				j=M+1;
			}
		}
	}
}
float CVarea(float dx)
{
	float dr,sf;
	dr=dx/L*(R-r)+r;
	sf=dr*dr*PI;
	return sf;
}