温度场有限容积法程序入门之三：2D温度场显式迭代计算（暂不考虑潜热）

          我们首先介绍温度场的求解吧，假设边界条件和初始条件已经设定。在贴代码之前，我们先谈谈这个类需要什么属性和行为：节点数组用于存储计算变量、网格大小、维度定义、计算函数，也就这么多了。如何计算某节点的温度？计算其东南西北方位相接节点对该节点的穿导热之和即可，读者这里可以考虑一下如何添加源相和对流换热进去。

package Soong.Solver
{
public class TSolver
{
public  var Tlist:Vector.<Node>;

private var xGridNum:uint = 1;//Number of Grid Allocated in X Direction
private var yGridNum:uint = 1;//Number of Grid Allocated in X Direction

public var dx:Number = 1;//Grid Size in X Direction
public var dy:Number = 1;//Grid Size in Y Direction

public var Sx:Number = 0;//Area of Heat Interface in X Direction
public var Sy:Number = 0;//Area of Heat Interface in X Direction

public var cellVol:Number = 0;//Volume of Control Volume

public var Freezing:Boolean=false;//If Time to Freeze

public function TSolver(xGridNum:uint,yGridNum:uint,dx:Number,dy:Number)
{
this.xGridNum = xGridNum;
this.yGridNum = yGridNum;

this.dx = dx;
this.dy = dy;

this.Sx = dy * 1;
this.Sy = dx * 1;

this.cellVol = dx * dy * 1;
}

public function Step(timeStep:Number):void
{
var col:uint = 0;
var row:uint = 0;
var node:Node = null;

for (col = 2; col < xGridNum - 2; col++ )
{
for (row = 2; row < yGridNum-2; row++ )
{
node = Tlist[Index(col, row)] as Node;

CalTnext(timeStep,node,col,row);

node.T0=node.T;
}
}
}

public function CalTnext(timeStep:Number,node:Node,col:uint,row:uint):void
{
var conner:Boolean=false;
var node_Adj:Node = null;

var conductionHeat:Number = 0;

//For Node on/in Connor or Edge
var SxFix:Number=1;//Area Fix Factor For Non-Interior Region in X Direction
var SyFix:Number=1;//Area Fix Factor For Non-Interior Region in Y Direction

var VolFix:Number=1;//Volume Fix Factor For Non-Interior Region in Y Direction

if(((col==2)&&(row==2))||((col==2)&&(row==yGridNum-3))||((col==xGridNum-3)&&(row==2))||((col==xGridNum-3)&&(row==yGridNum-3)))
{
SxFix=1/2.0;
SyFix=1/2.0;

conner=true;
}

if((col==2)||(col==xGridNum-3))
{
VolFix/=2;

if(!conner)
{
SyFix=1/2.0;
}
}

if((row==2)||(row==yGridNum-3))
{
VolFix/=2;

if(!conner)
{
SxFix=1/2.0;
}
}

node_Adj = Tlist[Index(col+1, row)] as Node;
conductionHeat+=node.eHeatExchangeFactor*(node_Adj.T0-node.T0)*Sx*SxFix;

node_Adj = Tlist[Index(col-1, row)] as Node;
conductionHeat+=node.wHeatExchangeFactor*(node_Adj.T0-node.T0)*Sx*SxFix;

node_Adj = Tlist[Index(col, row+1)] as Node;
conductionHeat+=node.nHeatExchangeFactor*(node_Adj.T0-node.T0)*Sy*SyFix;

node_Adj = Tlist[Index(col, row - 1)] as Node;
conductionHeat+=node.sHeatExchangeFactor*(node_Adj.T0-node.T0)*Sy*SyFix;

var dT:Number = conductionHeat * timeStep;
dT /= cellVol * VolFix * node.Rho * node.Cp;

node.T = node.T0 + dT;
}

public function LatentHeatRelease(node:Node):void
{

}

//Apply the Boundary Condition
public function ApplyBC():void
{

}

private function Index(col:uint=0,row:uint=0):uint
{
return row * xGridNum + col;
}
}
}

        简单吧，需要注意的是不同位置的节点其传热面积以及控制体体积不尽相同，需要Fix一下。可以预见，如果将SxFix、SyFix，VolFix设置为Node类的成员变量，计算会更快。这里给出初步的计算结果（迭代100s的结果）。目前笔者没有贴出所有代码，这时按照笔者提供的程序是无法运行的，读者想想，还缺点什么？



      将其对称得到整个界面：



    有点样子了，这还不是最终的计算结果，凝固潜热还没有考虑进去，可以使用物理意义明确的温度回升法处理。另外，我们没有离散偏微分方程，但是我们的方法和离散偏微分方程殊途同归的。也许读者可以理解有限差分和有限容积的连续与区别了。