加热炉板坯温度场数值模拟（二维，第三类边界条件）

加热炉板坯温度场数值模拟（二维，第三类边界条件）

现场热像仪只能监测板坯的表面温度，不能监测板坯内部温度。利用数值模拟机理建模模拟板坯加热过程，通过建模，计算出板坯的温度分布，对工艺改造提出建设性建议。

图与思路（摘自计算传热学大叔）

相关参数（摘自《步进式加热炉板坯温度场数值模拟》）

#include "stdafx.h"
//x方向控制体的个数
#define N 11
//y方向控制体的个数
#define M 5

//时间间隔
#define dtime 1
//最大时间步
#define maxtstep  int(21.43 * 60 + 0.5)

#include <iostream>
#include <iomanip>
#include <fstream>
using namespace std;

int _tmain(int argc, _TCHAR *argv[])
{
	ofstream file;
	file.open("Temperature.txt", ios::out | ios::trunc);
	file << fixed;
	double length,height; 
	double Tleft,Tright,Tbottom,Tup; 
	double Tmiddle; 
	double den, c, k;//物质密度,比热,导热系数
	double s;//热源
	double h;//对流换热系数
	double dx,dy;
	//左右边界2个体积为零的控制体和中间N个控制体
	//上下边界2个体积为零的控制体和中间M个控制体
	double Tinitial[N+2][M+2], T[N+2][M+2];
	double ae0[N+2][M+2], aw0[N+2][M+2], an0[N+2][M+2], as0[N+2][M+2], ap0[N+2][M+2], ap1[N+2][M+2], b[N+2][M+2];//系数
	int i,j;
	//加热炉内板坯尺寸15*1.123*0.135
	//将钢坯沿加热炉纵向取截面，即只考虑钢坯温度随宽度和厚度方向的变化
	length = 1.123; 
	height = 0.135;
	Tleft = 1150; 
	Tright = 1150;
	Tbottom = 1150;
	Tup = 1150;
	Tmiddle = 930;
	den = 7840;
	c = 465;
	k = 28.5;
	s = 0;
	h = 375; //换热系数
	dx = length / N; 
	dy = height / M;
	//中间控制体
	for (i = 2; i < N;i++) for (j = 2; j < M;j++)
	{
		ae0[i][j] = k * dy / dx; //右侧网格
		aw0[i][j] = k * dy / dx; //左侧网格
		an0[i][j] = k * dx / dy; //上侧网格
		as0[i][j] = k * dx / dy; //下侧网格
		ap0[i][j] = den * c * dx * dy / dtime - ae0[i][j] - aw0[i][j] - an0[i][j] - as0[i][j];
		ap1[i][j] = ap0[i][j] + ae0[i][j] + aw0[i][j] + an0[i][j] + as0[i][j];
		b[i][j] = s * dx * dy;
	}
	//边界控制体
	i = 1;//左侧网格
	for (j = 2; j < M;j++)
	{
		ae0[i][j] = k * dy / dx; //
		aw0[i][j] = 1/(1/(k * dy / (dx / 2))+1/h); 
		an0[i][j] = k * dx / dy;
		as0[i][j] = k * dx / dy;
		ap0[i][j] = den * c * dx * dy / dtime - ae0[i][j] - aw0[i][j] - an0[i][j] - as0[i][j];
		ap1[i][j] = ap0[i][j] + ae0[i][j] + aw0[i][j] + an0[i][j] + as0[i][j];
		b[i][j] = s * dx * dy;
	}
	i = N;//右侧网格
	for (j = 2; j < M;j++)
	{
		ae0[i][j] = 1/(1/(k * dy / (dx / 2))+1/h); 
		aw0[i][j] = k * dy / dx; 
		an0[i][j] = k * dx / dy;
		as0[i][j] = k * dx / dy;
		ap0[i][j] = den * c * dx * dy / dtime - ae0[i][j] - aw0[i][j] - an0[i][j] - as0[i][j];
		ap1[i][j] = ap0[i][j] + ae0[i][j] + aw0[i][j] + an0[i][j] + as0[i][j];
		b[i][j] = s * dx * dy;
	}
	j = 1;//下侧网格
	for (i = 2; i < N;i++)
	{
		ae0[i][j] = k * dy / dx; 
		aw0[i][j] = k * dy / dx; 
		an0[i][j] = k * dx / dy;
		as0[i][j] = 1/(1/(k * dx / (dy / 2))+1/h);
		ap0[i][j] = den * c * dx * dy / dtime - ae0[i][j] - aw0[i][j] - an0[i][j] - as0[i][j];
		ap1[i][j] = ap0[i][j] + ae0[i][j] + aw0[i][j] + an0[i][j] + as0[i][j];
		b[i][j] = s * dx * dy;
	}
	j = M;//上侧网格
	for (i = 2; i < N;i++)
	{
		ae0[i][j] = k * dy / dx; 
		aw0[i][j] = k * dy / dx; 
		an0[i][j] = 1/(1/(k * dx / (dy / 2))+1/h);
		as0[i][j] = k * dx / dy;
		ap0[i][j] = den * c * dx * dy / dtime - ae0[i][j] - aw0[i][j] - an0[i][j] - as0[i][j];
		ap1[i][j] = ap0[i][j] + ae0[i][j] + aw0[i][j] + an0[i][j] + as0[i][j];
		b[i][j] = s * dx * dy;
	}
	i = 1,j = 1;//左下
	{
		ae0[i][j] = k * dy / dx; 
		aw0[i][j] = 1/(1/(k * dy / (dx / 2))+1/h); 
		an0[i][j] = k * dx / dy;
		as0[i][j] = 1/(1/(k * dx / (dy / 2))+1/h);
		ap0[i][j] = den * c * dx * dy / dtime - ae0[i][j] - aw0[i][j] - an0[i][j] - as0[i][j];
		ap1[i][j] = ap0[i][j] + ae0[i][j] + aw0[i][j] + an0[i][j] + as0[i][j];
		b[i][j] = s * dx * dy;
	}
	i = 1,j = M;//右下
	{
		ae0[i][j] = k * dy / dx; 
		aw0[i][j] = 1/(1/(k * dy / (dx / 2))+1/h); 
		an0[i][j] = 1/(1/(k * dx / (dy / 2))+1/h);
		as0[i][j] = k * dx / dy;
		ap0[i][j] = den * c * dx * dy / dtime - ae0[i][j] - aw0[i][j] - an0[i][j] - as0[i][j];
		ap1[i][j] = ap0[i][j] + ae0[i][j] + aw0[i][j] + an0[i][j] + as0[i][j];
		b[i][j] = s * dx * dy;
	}
	i = N,j = 1;//
	{
		ae0[i][j] = 1/(1/(k * dy / (dx / 2))+1/h); 
		aw0[i][j] = k * dy / dx; 
		an0[i][j] = k * dx / dy;
		as0[i][j] = 1/(1/(k * dx / (dy / 2))+1/h);
		ap0[i][j] = den * c * dx * dy / dtime - ae0[i][j] - aw0[i][j] - an0[i][j] - as0[i][j];
		ap1[i][j] = ap0[i][j] + ae0[i][j] + aw0[i][j] + an0[i][j] + as0[i][j];
		b[i][j] = s * dx * dy;
	}
	i = N,j = M;
	{
		ae0[i][j] = 1/(1/(k * dy / (dx / 2))+1/h); 
		aw0[i][j] = k * dy / dx; 
		an0[i][j] = 1/(1/(k * dx / (dy / 2))+1/h);
		as0[i][j] = k * dx / dy;
		ap0[i][j] = den * c * dx * dy / dtime - ae0[i][j] - aw0[i][j] - an0[i][j] - as0[i][j];
		ap1[i][j] = ap0[i][j] + ae0[i][j] + aw0[i][j] + an0[i][j] + as0[i][j];
		b[i][j] = s * dx * dy;
	}
	//初始条件
	for (i = 1; i <= N; i++) for (j = 1; j <= M; j++)
	{
		Tinitial[i][j] = Tmiddle;
	}
	//边界条件
	//左侧虚拟控制体
	i = 0; 
	for (j = 1; j <= M; j++) Tinitial[i][j] = Tleft;
	//右侧虚拟控制体
	i = N+1; 
	for (j = 1; j <= M; j++) Tinitial[i][j] = Tright;
	//下侧虚拟控制体
	j = 0; 
	for (i = 1; i <= N; i++) Tinitial[i][j] = Tbottom;
	//上侧虚拟控制体
	j = M+1; 
	for (i = 1; i <= N; i++) Tinitial[i][j] = Tup;
	for (int tstep = 1; tstep <= maxtstep; tstep++)
	{
		//cout << tstep << '\t';
		for (j = 1; j <= M; j++) for (i = 1; i <= N; i++)  
		{
			T[i][j] = (ae0[i][j] * Tinitial[i + 1][j] + aw0[i][j] * Tinitial[i - 1][j] + 
				       an0[i][j] * Tinitial[i][j + 1] + as0[i][j] * Tinitial[i][j - 1] +
				       ap0[i][j] * Tinitial[i][j] + b[i][j]) / ap1[i][j];
			//cout << setprecision(5) << T[i][j] << '\t';
			cout << setprecision(5) << i<< '\t';
			cout << setprecision(5) << j<< '\t';
		}
		//cout << '\n';
		for (j = 1; j <= M; j++) for (i = 1; i <= N; i++)  
		{
			file << setprecision(5) << T[i][j] << '\t';
		}
		file << '\n';
		for (j = 1; j <= M; j++) for (i = 1; i <= N; i++)  
		{
			Tinitial[i][j] = T[i][j];
		}
	}
	return 0;
}

第三类边界条件运行结果

如果 1/h = 0；第三类边界条件变成了第一类边界条件。
在机理建模的过程中，参考了“B站计算传热学大叔”和“步进式加热炉板坯温度场数值模拟”的内容。
这篇博文的正确性需要考证，希望大家在这里指正。