温度场有限容积法程序入门之五：展望及问题

  首先总结一下以上[传输过程数值模拟可视化编程开发]程序的特点：

1，边界上的节点当作内部节点处理，所有节点公用一个迭代公式；

2，可以设置多种材质，多种材料时，注意及时更新每个节点的物性参数，如密度，比热，热导，一般使用调和平均，参见文献[1]。

      存在的问题：

1 ，笔者花了大约一天冒着猝死的危险写完的，时间仓促，通用性不够；

2，当前程序最致命的缺点在于似乎只能计算矩形区域和矩形网格，对于这些问题，可以参考直接差分法（DFDM，但是直接差分在网格扭曲度较大的情况下误差较大），可以有效处理不规则图形。

      其次，程序的拓展：

1， 3D拓展：

         节点需要除了东南西北增加上下两个热导，没做过，但是原理是一样的。

2，对流传热：

          节点热量来源除了扩散传导，源相外还有对流传热部分，如何计算对流引起的热量传递呢？其实传输现象的微分方程中对流项是最难以计算的，表面上看是一个一阶导数，实际上它还是个一阶导数，只不过需要注意传导热和对流热的量级关系，只有他俩量级相当，使用中心差分才能得到有意义的解。

 

          这里讲一个小故事：狼羊同饮，羊位于下游，狼意欲吃羊，曰：“你他妈污染老子的水源，拿命来”，羊说“小弟我在下游怎么会污染大哥的水源”。如果我是狼而且懂的Pe准数（读者快查Pe准数的意义）的话，也不会浪费口舌去吃。在羊的思维里，下游是不会污染上游的，智商低下，果然可以被任何肉丝动物吃。事实上，当水速很低，下游会通过扩散污染上游的，所以嘛，狼吃羊就是应该的，不懂还死搅蛮缠。于是人们就发明了一阶迎风、混合、二阶迎风、指数、QUICK等格式去离散对流项，以考虑上游对当前节点传输行为的重要性。话说笔者年青时小解顺风冒三丈，可见流动的上游确实对传输现象有不可忽略的影响，如今老了：“廉颇老矣，尚能尿远？”。

        曾经一位老师大吹自己年轻时流体流动传输编程多么强悍，后来我一个师弟问他什么是“一阶迎风格式”，他说不知，我笑了。

       扯远了，编程时如果考虑对流换热，请查阅文献[1]，里面有离散格式，在本文后面附上关于1D对流-扩散方程迭代计算学习笔记摘录一二，就不在这里班门弄斧了。

3，非线性材料：

         当物性参数随温度或位置等参数变化时，如何计算？留给读者思考吧，笔者也不擅长。有人说将物性参数线性化，其实未必完全妥当，因为有时就是非线性的，线性化了，反而误差大了。提示一下，更新物性参数，多次迭代，直到两次计算结果接近到一定程度，此时，该时刻下方程的解收敛了，以该收敛解为初值可以迭代计算下一个时间步长的解了。

 

4，非均匀网格：

         这个也留给读者吧，需要对过程的物理意义十分明晰。如果对整个物理方程的建立都不清楚，肯定算不对的。本文引入热导的概念，使得计算在一定程度上得到简化。

         

       到此，贴出来的代码可以独立运行了。

  附1： 一维扩散-对流方程迭代格式学习笔记（摘录自参考文献1）：

    比如2D对流-扩散方程（pp.177）如下：

      对其在控制体体积内积分（pp.177）：

     整理为代数迭代方式，只考虑1D情形（pp.178）：

    其中，不同坐标系下系数计算不同（pp.179）：

    上述表达式中有函数A（），其计算根据不同离散格式而不同，具体如下（pp.151）：

需要特别注意的是：与积累项和源项不同的是，扩散项和对流项是从控制体边界进入到控制体的，所以Ae，Aw的计算使用到的物性参数及速度不是主节点的物性参数和主节点的速度（使用同位网格的话），而是需要插值计算得到的边界上的值，具体加权平均计算公式要根据具体物理意义而定。   

 

    对于已知气体压强分布的一维N-S方程，就不用SIMPLE算法迭代计算压强了，如果使用同位网格，N-S方程解等同于对流-扩散方程的求解。附带介绍一个自己开发的程序Demo：

    附2，FAQ：

Q1：为什么用热导而不是热阻的概念？

A：其一通用性；其二笔者印象里，乘法的CPU时间小于除法；也许也有人会继续问，你算热导不是也用了除法吗？免不了用除法，在本算例下，只在初始化时用除法算了热导一批，而迭代过程，很多循环中用了乘法，利大于弊。

 

Q2：为什么网格四周添加2层多余的网格，是为了浪费内存吗？

A：比如对端点处某量进行边界条件处理，一阶导数、二阶导数等，对流项中高阶离散格式（比如QUICK格式）的计算会用到。

 

Q3：为什么没有计算迭代的稳定性？

A：楼主很懒惰，爱收敛不收敛，收敛不了把时间步长弄小。我计算每一节点的最大时间步长，浪费多少CPU时间啊。

 

 Q4：为什么不用隐式迭代方法，据说这种方法绝对稳定？

A：因为显示迭代格式物理意义清晰，隐式迭代格式精度与解析解存在一定误差，一般用C-N格式。

这里给出一段显示迭代过程可能会用到的代码，切忌不要草草照搬代码，我的代码main函数中矩阵系数A中首尾项在程序中无用，用于求解上对焦矩阵线性方程组，而且只是一维的；二维的是5对角矩阵；三维的是7对角矩阵。

// TDMA.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"

#define REAL float

#define idxA(i,j) (i*2+j-2)
#define idxL(i,j) (i+j-2)
#define idxU(i,j) (i-1)
#define idxb(i) (i-1)
#define idxz(i) (i-1)

void TDMASolver(int dim,REAL* A,REAL*b,REAL* root);

void main()
{
	REAL A[12]={0,1,2,
			  2,3,4,
			  4,5,6,
			  6,7,8};

	REAL b[4]={3,9,15,13};

	REAL root[4];

	TDMASolver(4,A,b,root);

	printf("Hello World!\n");

	
}

void TDMASolver(int dim,REAL* A,REAL* b,REAL* root)
{
	int i=0;

	REAL* L=new REAL[dim<<1-1];
	REAL* U=new REAL[dim-1];
	REAL* z=new REAL[dim];

	// Step 1
	L[idxL(1,1)]=A[idxA(1,1)];
	U[idxU(1,2)]=A[idxA(1,2)]/L[idxL(1,1)];
	z[idxz(1)]=b[idxb(1)]/L[idxL(1,1)];
	
	//Step 2
	for (i=2;i<dim;i++)
	{
		L[idxL(i,i-1)]=A[idxA(i,i-1)];
		L[idxL(i,i)]=A[idxA(i,i)]-L[idxL(i,i-1)]*U[idxU(i-1,i)];
		U[idxU(i,i+1)]=A[idxA(i,i+1)]/L[idxL(i,i)];
		z[idxz(i)]=(b[idxb(i)]-L[idxL(i,i-1)]*z[idxz(i-1)])/L[idxL(i,i)];
	}
	
	//Step 3, Now i=n
	L[idxL(i,i-1)]=A[idxA(i,i-1)];
	L[idxL(i,i)]=A[idxA(i,i)]-L[idxL(i,i-1)]*U[idxU(i-1,i)];
	z[idxz(i)]=(b[idxb(i)]-L[idxL(i,i-1)]*z[idxb(i-1)])/L[idxL(i,i)];

	//Step 4
	root[dim-1]=z[dim-1];
	
	//Step 5
	for (i=dim-1;i>0;i--)
	{
		root[i-1]=z[idxz(i)]-U[idxU(i,i+1)]*root[i];
	}

	delete[] L;
	delete[] U;
	delete[] z;
}

 

Q5：用AS3写性能低下？

A：必须的，Adobe太懒惰了，停滞不前，曾经答应使AS3性能达到C的80%，后来取消了，昨天Unity还取消支持Flash了，笔者也很郁闷。不过如果非要使用Flash以便于网络展示，可以使用其FlasCC编译技术；如果希望使用其低学习成本的3D技术且主要为桌面部署，还是在应用程序中嵌入Flash控件吧。最后还是希望大家使用原生语言开发本地（Native）程序，而不是脚本语言。

参考文献：

1，陶文铨：数值传热学(第2版)