参考2多维导热问题的数值解原理

第二章多维导热问题2.1二维非稳态导热全隐格式的通用离散方程三种二维坐标系中的网格系统见下图2-1。采用控制容积积分法导出的离散方程以二维直角坐标系圆柱轴对称坐标系图2-1三种坐标系极坐标系直角坐标系圆柱轴对称坐标系图2-1三种坐标系极坐标系直角坐标系下的为例，根据二维非稳态导热方程：2.1)取全隐格式，假设节点之间温度线性分布，界面上热流密度均匀分布。非稳态项积分：dxdydt=(pc)(T一T0)AxAydt pppswt扩散项积分：鬥J£広明）dxdydt+鬥埠広寻）dxdydtswswsw.T一T.T一T.T一T.T一T九—E P—九e (8x)e-—P Ww(5x)wAyAt+九nT—T.T—T—N P一九—P (Sx) s(Sx)AxAt源项积分：tneJJJSdxdydt源项积分：tneJJJSdxdydt=(Sc+ST)AxAyAtPPtsw上述结果整理成aT

P二aaT

P二aPET+aT+aT+aT+bEWWNNSS2.2)其中各系数为：Aya=E (5x)/九Aya=E (5x)/九ee=a+a+a+a+a0一SAxAySWNPP2.2)AyAx(5x)/九'ww@y)/九’nnAx(Sy)/九ss2.3a)2.3b)让AxA2.3b)At2.3c)b=SAxAy+a0T02.3c)cPP仍然需要记住，式（2.3a）表示的是各节点之间的热导（热阻的倒数），分子上的Ax、Ay代表的是各控制容积面上的面积；在二维问题中，AxAy的乘积是控制容积的体积。a0代表P的是控制容积的热惯性。由此可见，利用上述系数计算式的物理含义，很容易写出三维导热问题的离散化方程及它的系数。对于圆柱轴对称坐标和极坐标，同样可以利用系数的物理含义写出各系数计算式，离散方程与式（2.2）相同。不过要注意，在圆柱轴对称坐标中，选用一个弧度角的范围，极坐标取垂直于纸面一个单位长度（lm）。这样三种坐标系下的离散方程的系数可以表示为表2.1以便于编写统一的计算程序。二维导热问题中三种坐标系中系数的通用表达式 表2.1坐标系直角圆柱轴对称极坐标通用表达式东西坐标xx9X南北坐标yrrY半 径1rrR东西尺度系数11rSX东西节点间距5x5xr59(5X)(SX)南北节点间距5y5r5r5Y东西导热面积AyrArArR(AY)/SX南北导热面积AxrAxrA9R(AX)控制容积体积AxAyrArAxrA9ArR(AX)(AY)a(a)W EAyrArArR(AY)(SX)2(5x)/九e e(5x)/九e e(5x)/九e e(59)/九e ea(a)N SAxrArrA9R(AX)(5y)/九n n(5y)/九n nn(5r)/九n nn(5r)/九n na0p(pc)R(AX)(AY)/AtpbSR(AX)(AY)+aoT0c PPapa+a+a+a+a0—SR(AX)(AY)E W N S P P面得到的是计算域内内节点的离散化方程，对于边界节点，可以采用边界控制容积热平衡方程导出节点方程。

2.2边界节点方程第一类边界条件是给定边界上的温度值，所以求解区域是内接点方程组第二类边界条件给出的是边界上的热流密度，通常表示为这样的表达在求解时还不能直接引入到节点上，需要根据能量守恒方程变换为dT一dnB2.4)第三类边界条件为对流换热条件，已知参数为边界面上的对流换热系数和流体温度，表示为q=h(T-T)B Bf同样需要经过变换后才能进行计算，一般变换成-卫

dnB=h(T-2.2边界节点方程第一类边界条件是给定边界上的温度值，所以求解区域是内接点方程组第二类边界条件给出的是边界上的热流密度，通常表示为这样的表达在求解时还不能直接引入到节点上，需要根据能量守恒方程变换为dT一dnB2.4)第三类边界条件为对流换热条件，已知参数为边界面上的对流换热系数和流体温度，表示为q=h(T-T)B Bf同样需要经过变换后才能进行计算，一般变换成-卫

dnB=h(T-T)Bf2.5)容易看出，第二类和第三类边界条件根据上述表达可以用统一的方式(边界上的热流密度)离散，参见图2-2。对于P节点，若采用显格式并考虑有内热源：qAy+a(T0-T)+a(T0-T)+a(T0E P NN P SS-T)+(S+ST)AxAy=a0(T-T0)PPPPP整理后：aTP=aT0+aT0+aT0+bSSEENN2.6)其中：Ay(5x)/九’eea=NAx@y)/九’nnAx(Sy)/九'ssa0P=屮AxAy,AtC•p•EJ=a+a+aEN+a0-SAxAy,PPTOC\o"1-5"\h\z\o"CurrentDocument"b=SAxAy+a0T0-qAy (2.7)c PP B若采用隐格式离散方程为：\o"CurrentDocument"aT=aT+aT+aT+b (2.8)PPEENNSS方程中的系数与上相同。这样，对于第二类边界问题，边界面上的温度被排除在外，待计算完毕后通过插值方式获得。对于第三类边界条件，容易看出，从流体到P节点的传热热阻有两部分组成，半个控制容积的导热热阻和边界面上的对流换热热阻，即边界上的热流为：q二 f p (29)^b1/h+(8x)/九wB同样可以将边界面上的温度排除在外，最后才插值计算获得。上述处理结果，使得内节点和靠近边界的节点的代数方程取得了相同的形式，只不过靠近边界的节点方程相应有一个系数为0。2.3代数方程的求解方法求解线性方程组的两类方法是直接求解和迭代求解，直接求解是通过一次计算来获得代数方程的精确解，但计算工作量特别大；迭代计算是将计算分成许多轮次，每次计算量减少，只要迭代方式组织合理，可以获得比直接解法更好的经济性，在计算流体力学和传热学中经常采用这种方式，尤其在节点数很大时，即使收敛慢的迭代方法也可能比消元法更加有效在迭代计算中有两个问题，一是迭代的收敛性问题；其次是如何加快迭代速度问题。一般对如导热这一类问题，迭代收敛条件为工a 亠<1 (2.10)aP由于采用有限容积法生成的离散方程，这一条件上述条件一定是满足的。常用的迭代方式中，Jacobi法的收敛速度最慢，Gauss-Seidel迭代比较快。交替方向线迭代方法(ADI)是最有利的迭代方法。2.4边界上不规则区域的处理方法常见的处理方法有阶梯形边界逼近真实边界、坐标变换法、区域扩充法及边界节点单独建立方程法等。在采用商用软件时这些方法通常不需要我们专门去考虑。例题1：假设图2-3所示的矩形截面肋片，肋片根部为 h丁TOC\o"1-5"\h\z温度T0,肋的上下两侧及端部为对流换热条件。试研 1—究在不同的Bi(=h6/A)数下肋片中的温度分布，并T°_ 角!比较数值计算所得出的导热量与按一维假定得出的导 热量的区别。 *L”解：该问题要计算肋片的导热量，实际是肋片表面与周围流体之间的对流换热量，可以在肋片根部取得导 图2-3例题1附图热量的数据，所以先要计算出肋片内的温度场。计算条件按照Bi数确定，分成2、1、0.1和0.05等4组。假定L/(26)=4，并取肋片材料导热系数为16.27W/(mK)，比热容为502.48J/(kg・K)，密度为8030kg/m3。序号Bi数h/W/(m2・K)Q/WL=80mm。具体参数见右表。左边界给1232542X1536.4定温度373K，周围流体温度293K。利2116272X1168.98用FLUENT求解。30.1162.72X402.511.利用GAMBIT建立计算几何模40.0581.352X273.73型方法见附录1。肋片长度为8个单位，厚度为2个单位(建模过程中根据肋片导热上下的对称情况，只画出一半，即1个厚度单位)。边界条件，对流面根据上述列出的Bi数确定对流换热系数。计算结果分析:如果将该问题作为一维模型计算，当端面绝热时，端面温度为305.4K,数值计算得到的端面温度上下各为306.4K和305.8K,平均约为306.1K,两者相对差别为0.26%；导热量计算一维模型的为821.4W,数值计算的为2X400.2W,两者相对差别为2.6%。所以在Bi数为0.1时,数值计算结果表明此条件下肋片可以简化为一维模型。3.739+02369AO23.659*02数值计算得到的端面温度上下各为306.4K和305.8K,平均约为306.1K,两者相对差别为0.26%；导热量计算一维模型的为821.4W,数值计算的为2X400.2W,两者相对差别为2.6%。所以在Bi数为0.1时,数值计算结果表明此条件下肋片可以简化为一维模型。3.739+02369AO23.659*023.619+023.53e+013.4叶023.455+023419*02S.379+M3.339+023.19e+l>2325e+U23.215+023.179+023i3e+i>23C®e*02305102301e*022.970+022.939+02Y\373X02369e+02365e-*023.61S4023.578*023.53e-t02349®*023456-^023.41e*02332*023.33X023.刘屮02325e*023.210+023.1血Q2343*023.090*023.05e-tOJ301e-*02297。*02295e-*02Bi=2.0Bi=1.03.73®"023.70e-f-02364e+02362+02358e*02350"02347ZQ2344e+02341e-*023.35®+02332e+02329a02326e*02323e+02320^023.178+023.55e+023.52e+02367e-*02314e+02Bi=0.1 Bi=0.05图2-4不同Bi数下肋片内部温度分布Bi=2.0Bi=1.0Bi=2.0Bi=1.0Bi=0.1图2-5不同Bi数下上下两面温度分布图从计算结果分析，在Bi数小于0.1条件下，一维模型是对实际的一个近似，但肋片内

的温度分布的二维属性仍然存在，即便在Bi数为0.05的条件下，肋片的上表面和中心对称面上的温度还是有微小的差异。例题2：一钢锭，大小为0.5X0.7X1.0m,初始温度均匀为20°C,其导热系数为40.5W/（mK）,密度7900kg/m3，比热容710.05J/（kg・K），试确定将其置入1200C的加热炉中4小时后的最低温度和最高温度。假设炉内烟气与钢锭之间的换热系数为348W/（m2・K）。解：这是三维非稳态导热问题，导热体外表是对流换热条件。可以先估计，最高温度在顶角上，最低温度在中心点上角（见图2-6中顶角点和中心点）。计算区域可以按对称性选取为原物体的1/8，见图2-6。正面、右面和顶面为对流边界面，背面、左面及底面为对称面（绝热面）。图2-6导热体图2-7六面体设置对话框1）确定导热体大小如图2-7所示。在宽度（Width）栏、启动GAMBIT,同卩Brick，弹出对话框深度（Depth）栏、高度（Height）栏中分别输入0.125、0.175和0.25图2-6导热体图2-7六面体设置对话框1）确定导热体大小如图2-7所示。在宽度（Width）栏、启动GAMBIT,同卩Brick，弹出对话框深度（Depth）栏、高度（Height）栏中分别输入0.125、0.175和0.25；点击Apply、Close，形。点击FittoWindowI按钮查看图，打开图2）划分网格：依次点击2-8所示对话框。点击Volumes右侧的黄色框后，用Shift+鼠标左键点击工作显示区内几何图形的边线后，在Spacing项选择Intervalsize，输入0.005；点击Apply、Close，所划出的网格为25X35X5003）设置边界类型：为方便，先将网格显示关闭，点击在SpecifyDisplayAttributes中点击Mesh项中的Off,点击确定,点击Close；依次点击对话框中加入六个面的边界面属性，分别在Name项中给定边界面名称、Type、选定边界、点击Apply。右边界Right，属性Wall，在工作区中Shift+鼠标左键点右边界线变为红色，点击Apply；各边界面属性分别如下表：，在SpecifyBoundaryTypes图2-8网格设置对话框名称RightLeftTopBottomFrontBackTypeWallSymmetryWallSymmetryWallSymmetry最后点击Close，存上Mesh文件，退出GAMBIT。三维建模的步骤比二维要少，原因是建模的开始是以面开始的，而二维要将线转化为面。2•启用FLUENT:按条件分别确定材料、边界、初始等具体数据，下面是时间步长1秒所获得的数据：AverageofSurfaceVertexValues(k)TimeStepConvergencehistoryofStaticTemperatureontopright(Time=1,4400e+04)Jun21,2007FLUENT6.3(3d,pbns,lam,unsteady)1600.00001400.00001200.00001000.0000Jun21