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

1、精品文档你我共享AAAAAA)dxdydt -yIL (x)eTP TW, LAyAt +ru.(x)wUTN -Tp r Tp _TNn（、X）n.X筑(x)n(2.2)第二章多维导热问题2.1二维非稳态导热全隐格式的通用离散方程三种二维坐标系中的网格系统见下图2-1。采用控制容积积分法导出的离散方程以二维图2-1三种坐标系直角坐标系下的为例，根据二维非稳态导热方程：(2.1)取全隐格式，假设节点之间温度线性分布，界面上热流密度均匀分布。非稳态项积分:n e t：盘I I I c dxdydt = （，c）P （TP -TP） ：x ：ys w t扩散项积分:t T n e -)dxdydt

2、(Xt sw 內源项积分:t 7 n eSdxdydt = (Sc SpTp) ：x y tt s w上述结果整理成:apTp = &ETE awTw aNTN asTs b其中各系数为:极坐标系精品文档你我共享AAAAAAaE(XM e aP = aEaw a”aW-C；X)w/wasaNaP0 _ Spxiy(2.2)0aP二-，aS(y)n / nCy)s / s(2.3a)(2.3b)b = Sc.x：y aP0Tp0（2.3c）仍然需要记住，式（2.3a）表示的是各节点之间的热导（热阻的倒数），分子上的ex、厶y代表的是各控制容积面上的面积；在二维问题中，：x ：y的乘积是控制容积的

3、体积。aP0代表的是控制容积的热惯性。由此可见，利用上述系数计算式的物理含义，很容易写出三维导热问题的离散化方程及它的系数。对于圆柱轴对称坐标和极坐标，同样可以利用系数的物理含义写出各系数计算式，离散方程与式（2.2）相同。不过要注意，在圆柱轴对称坐标中，选用一个弧度角的范围，极坐标取垂直于纸面一个单位长度（1m）。这样三种坐标系下的离散方程的系数可以表示为表2.1以便于编写统一的计算程序。二维导热问题中三种坐标系中系数的通用表达式表2.1坐标系直角圆柱轴对称极坐标通用表达式东西坐标xxeX南北坐标yrrY半径irrR东西尺度系数i1rSX东西节点间距5xrd6(SX)(SX)南北节点间距勿靳

4、dr6Y东西导热面积也yrArrR(AY)/SX南北导热面积zrAxreRgX)控制容积体积xAyrArAxr也印lrR(AX )(也Y)RQY)r冷AraW (aE)(&X)e/爲(叽/(胡)ee(SX)2(6x)eexrn也r侶aN ( as )()n 入(6r)n n(r)n 入RQX)()n/纭n0ap(Pc)pR(AX)QY)/t精品文档你我共享AAAAAAbScR3X)9Y)+ap0Tp0apaE +aW +aN +as +aP。_SPRQX)(AY)上面得到的是计算域内内节点的离散化方程，对于边界节点，可以采用边界控制容积热平衡方程导出节点方程。精品文档你我共享AAAAAAcn

5、B(2.4)整理后:0apTp = QETE aN TN QSTS b(2.6)其中:aE(X)e / eaNaSXap0y，taP=aEa” a(2.7)(2.8)2.2边界节点方程第一类边界条件是给定边界上的温度值，所以求解区域是内接点方程组。 第二类边界条件给出的是边界上的热流密度，通常表示为qB 二 f (x, y,z)这样的表达在求解时还不能直接引入到节点上，需要根据能量守恒方程变换为第三类边界条件为对流换热条件，已知参数为边界面上的对流换热系数和流体温度，表示为qB = h(TB - Tf )同样需要经过变换后才能进行计算，一般变换成FT九=h(Tb Tf)( 2.5)cn B容易

6、看出，第二类和第三类边界条件根据上述表达可以用统一的方式(边界上的热流密度)离散，参见图2-2。对于P节点，若采用显格式并考虑有内热源：00亠丄丄0qB yaE (TE- TP ) aN(TN- TP )aS(TS- TP )(SSpTp)：Xy =ap(Tp -Tp)b =Sc x y ap0Tp0 - qB y若采用隐格式离散方程为：apTp 二 aEE aNN a$Ts b方程中的系数与上相同。这样，对于第二类边界问题， 边界面上的温度被排除在外， 待计算完毕后通过插值方式 获得。对于第三类边界条件，容易看出，从流体到 P节点的传热热阻有两部分组成，半个 控制容积的导热热阻和边界面上的对

7、流换热热阻，即边界上的热流为：图精品文档你我共享AAAAAAqB =1/h(X)w/ B(2.9)同样可以将边界面上的温度排除在外，最后才插值计算获得。上述处理结果，使得内节点和靠近边界的节点的代数方程取得了相同的形式，只不过靠近边界的节点方程相应有一个系数为0。2.3代数方程的求解方法求解线性方程组的两类方法是直接求解和迭代求解，直接求解是通过一次计算来获得代数方程的精确解，但计算工作量特别大；迭代计算是将计算分成许多轮次，每次计算量减少,只要迭代方式组织合理，可以获得比直接解法更好的经济性，在计算流体力学和传热学中经常采用这种方式，尤其在节点数很大时，即使收敛慢的迭代方法也可能比消元法更加

8、有效。在迭代计算中有两个问题，一是迭代的收敛性问题；其次是如何加快迭代速度问题。 般对如导热这一类问题，迭代收敛条件为 |anap-1(2.10)由于采用有限容积法生成的离散方程，这一条件上述条件一定是满足的。常用的迭代方式中，Jacobi法的收敛速度最慢，Gauss-Seidel迭代比较快。交替方向线 迭代方法（ADI ）是最有利的迭代方法。2.4边界上不规则区域的处理方法常见的处理方法有阶梯形边界逼近真实边界、坐标变换法、区域扩充法及边界节点单独建立方程法等。在采用商用软件时这些方法通常不需要我们专门去考虑。例题1假设图2-3所示的矩形截面肋片， 肋片根部为 温度T。，肋的上下两侧及端部为

9、对流换热条件。试研 究在不同的Bi （= h S/入）数下肋片中的温度分布，并 比较数值计算所得出的导热量与按一维假定得出的导 热量的区别。解：该问题要计算肋片的导热量，实际是肋片表面与 周围流体之间的对流换热量，可以在肋片根部取得导r.-一图2-3例题1附图热量的数据，所以先要计算出肋片内的温度场。计算条件按照Bi数确定，分成 2、1、0.13502.48J/(kg K)；密度为 8030kg/m。序号Bi数h/W/(m 2 K)Q/WL=80mm。具体参数见右表。左边界给1232542 X 1536.4定温度373K，周围流体温度293K。禾U2116272X 1168.98用FLUENT

10、求解。30.1162.72 X 402.511.利用GAMBIT建立计算几何模40.0581.352 X 273.73和0.05等4组。假定 L/（2 S ）=4，并取肋片材料导热系数为16.27W/（mK），比热容为型方法见附录1。肋片长度为8个单位，厚度为2个单位（建模过程中根据肋片导热上下的对称情况，只画出一半，即1个厚度单位）。2 边界条件，对流面根据上述列出的Bi数确定对流换热系数。3计算结果分析：如果将该问题作为一维模型计算，当端面绝热时，端面温度为305.4K ,精品文档你我共享AAAAAA精品文档你我共享AAAAAAidAamHtrfI4F*09 -U*町w/曲tt鼻K3-J?

11、4K3-JtaR?3金tt3 lbg1. fFtn瞰9-OhK*D?MMBi=0.1从计算结果分析，在数值计算得到的端面温度上下各为306.4K和305.8K，平均约为306.1K，两者相对差别为0.26%;导热量计算一维模型的为821.4W，数值计算的为2X 400.2W，两者相对差别为2.6%。所以在Bi数为0.1时，数值计算结果表明此条件下肋片可以简化为一维模型。Bi=0.1图2-5不同Bi数下上下两面温度分布图Bi数小于0 .1条件下，一维模型是对实际的一个近似，但肋片内3 蓟33+lr4323SBi=0.05图2-4不同Bi数下肋片内部温度分布PMitian mJPcsitioft

12、(m)精品文档你我共享AAAAAA的温度分布的二维属性仍然存在，即便在 面上的温度还是有微小的差异。Bi数为0.05的条件下，肋片的上表面和中心对称例题2: 一钢锭，大小为0.5X 0.7X 1.0m ,初始温度均匀为20C,其导热系数为40.5W/（m -K）,密度7900kg/m3，比热容710.05J/（kg K），试确定将其置入1200 C的加热炉中4小时后的最低温度和最高温度。假设炉内烟气与钢锭之间的换热系数为348 W/（m 2 K）。解：这是三维非稳态导热问 题，导热体外表是对流换热 条件。可以先估计，最高温 度在顶角上，最低温度在中 心点上角（见图2-6中顶角 点和中心点）。计

13、算区域可 以按对称性选取为原物体 的1/8,见图2-6。正面、 右面和顶面为对流边界面， 背面、左面及底面为对称面（绝热面）。1.利用 GAMBIT 建 模：1）确定导热体大小：启动 如图2-7所示。GAMBIT ,卬一匸戸）Brick，弹出对话框0.25;点击 Apply、Close，点击 Fit to Window形。2）戈U分网格：依次点击Mesh VolumesVolumes兀回簟按钮，打开图2-8所示对话框。点击 Volumes右侧的黄色框后，用 Shift+鼠标 左键点击工作显示区内几何图形的边线后，在Spaci ng项选择Interval size，输入0.005;点击Apply

14、 Close,所划出的网格为25X 35X 50。3）设置边界类型：为方便，先将网格显示关闭,点击口在Specify Display Attributes中点击 Mesh项中的 Off,点击确定， 点击Close;依次点击丨辛 I ，在Specify Boundary Types对话框中加入六个面的边界面属性，分别在Name项中给定边界面名称、Type、选定边界、点击 Apply。右边界Right，属性Wall，在工作区中Shift+鼠标左键点右边界线变为红色，点击Spacing：F ftpply DefaultJo.00intervalOptions:F MeshRemove old me：

15、FseniOO ir?w（-）rnJ Ignore size func | RMET | Cose I名称RightLeftTopBottomFrontBackTypeWallSymmetryWallSymmetryWallSymmetry在宽度（Width ）栏、深度（Depth）栏、高度（Height）栏中分别输入 0.125、0.175和三维建模的步骤比二维要少， 原因是建模的开始是以面开始的，而二维要将线转化为面。图2-6导热体图2-7六面体设置对话框Scheme: Appl/ Defaull图2-8网格设置对话框Apply ;各边界面属性分别如下表:最后点击 Close，存上 Mes

16、h文件，退出 GAMBIT。精品文档你我共享AAAAAA1 BOO .0000Average of Surface Vertex Values (k)Time StepConvergsince history of Static Temperature on mnerpoint (Time=1 -4400a+D4)Jun 21,2007FLUENT 6.S (3d, pbns, lam, unsteady)2启用FLUENT :按条件分别确定材料、边界、初始等具体数据，下面是时间步长 秒所获得的数据：Average of Surface Vertex Values (k)Time StepCo

17、nverg 的饰 history of Static Tempei：ratur on toprighl ( Tiirii= L4400e-H4fJun 21,2007FLUENT 6,3 (3d. pbns. lam, unsteady)iiKjnitor-21400.000012QO.OQOO1QOD.OOOOBOO .0000600.0000400.00002QD.OOOD精品文档你我共享AAAAAAAverage of Surface Vertex Values (k)Conv&rghistory of Static Temperatur on topighi (Time =3.0DOe

18、+03*Jun 21,2007 FLUENT 63 (3d. pbns, lamr unsteady!iiKjnitor-?1100.0000Average of SurfaceVertexValues/JUn 21,2007FLUENT 6,3 (3d, pbns, lam, unsteady)Time StepConvrgince history of Static Temperature on mnerpoint (Timj&=3.6000e+D3 JVTTT、 - 、；要注意：1） 该问题可以通过精确解获得数据，顶角点温度1470K，中心点温度1431K;2）因计算时间长，故开始计算时可以考虑两个问题：一是时间步长多少为好？二是空 间步长是多少为好？习题：1000.MM900.0000soo.oooo700.0000600.00005QD.OOOO40