1、1第二章 多维导热问题2.1 二维非稳态导热全隐格式的通用离散方程三种二维坐标系中的网格系统见下图 2-1。采用控制容积积分法导出的离散方程以二维直角坐标系下的为例，根据二维非稳态导热方程：SyTxtTc)()(（2.1）取全隐格式，假设节点之间温度线性分布，界面上热流密度均匀分布。非稳态项积分： yxTcdxyttTcPnsewtt )() 0扩散项积分： ttnsewttnsew dxytyxy)()( txTTtTT nNPsnNWPeEe )()()()( 源项积分： tyxSSdxytPcttnsew )(上述结果整理成：bTaTaTaSNWEP （2.2）其中各系数为：yxPSNW
2、EP 0（2.2）， ， ， eExya/)(wxya/)(nNya/)(sSyxa/)(直角坐标系 圆柱轴对称坐标系 极坐标系图 2-1 三种坐标系2（2.3a）yxtcaPP) (0（2.3b）0PcTaSb（2.3c）仍然需要记住，式（2.3a）表示的是各节点之间的热导（热阻的倒数） ，分子上的 、x代表的是各控制容积面上的面积；在二维问题中， 的乘积是控制容积的体积。y yx代表的是控制容积的热惯性。由此可见，利用上述系数计算式的物理含义，很容易写0Pa出三维导热问题的离散化方程及它的系数。对于圆柱轴对称坐标和极坐标，同样可以利用系数的物理含义写出各系数计算式，离散方程与式（2.2）相
3、同。不过要注意，在圆柱轴对称坐标中，选用一个弧度角的范围，极坐标取垂直于纸面一个单位长度（1m ） 。这样三种坐标系下的离散方程的系数可以表示为表 2.1 以便于编写统一的计算程序。二维导热问题中三种坐标系中系数的通用表达式 表 2.1坐标系 直角 圆柱轴对称 极坐标 通用表达式东西坐标 xxX南北坐标 yrrY半 径 1 R东西尺度系数 1 1 S东西节点间距 xxr)(南北节点间距 yrY东西导热面积 rSXR/)(南北导热面积 xxr控制容积体积 yr)(Y( )WaEex/)(ePxr/)(e/)( exSXR/)(2( )NSny/)( nr/)(nr/)(ny/)(0Pa tYXR
4、cPb 0)(PcTaS3Pa )(0YXRSaaPSNWE 上面得到的是计算域内内节点的离散化方程，对于边界节点，可以采用边界控制容积热平衡方程导出节点方程。2.2 边界节点方程第一类边界条件是给定边界上的温度值，所以求解区域是内接点方程组。第二类边界条件给出的是边界上的热流密度，通常表示为),(zyxfqB这样的表达在求解时还不能直接引入到节点上，需要根据能量守恒方程变换为BnT（2.4） 第三类边界条件为对流换热条件，已知参数为边界面上的对流换热系数和流体温度，表示为 )(fBThq同样需要经过变换后才能进行计算，一般变换成)(fBn（2.5）容易看出，第二类和第三类边界条件根据上述表达
5、可以用统一的方式（边界上的热流密度）离散，参见图 2-2。对于 P 节点，若采用显格式并考虑有内热源： )()()()()( 00000 PPCPSPNEB TayxTSTaTaTayq 整理后：bSNEP000（2.6）其中：， ，eExya/)(nNyxa/)(， ，sSy/)(tcPP 0，yxSaaSNEP 0qTyxbBPc 图 2-2 边界节点的离散图4（2.7）若采用隐格式离散方程为：bTaTaSNEP（2.8）方程中的系数与上相同。这样，对于第二类边界问题，边界面上的温度被排除在外，待计算完毕后通过插值方式获得。对于第三类边界条件，容易看出，从流体到 P 节点的传热热阻有两部分
6、组成，半个控制容积的导热热阻和边界面上的对流换热热阻，即边界上的热流为：BwPfBxhTq/)(/1（2.9）同样可以将边界面上的温度排除在外，最后才插值计算获得。上述处理结果，使得内节点和靠近边界的节点的代数方程取得了相同的形式，只不过靠近边界的节点方程相应有一个系数为 0。2.3 代数方程的求解方法求解线性方程组的两类方法是直接求解和迭代求解，直接求解是通过一次计算来获得代数方程的精确解，但计算工作量特别大；迭代计算是将计算分成许多轮次，每次计算量减少，只要迭代方式组织合理，可以获得比直接解法更好的经济性，在计算流体力学和传热学中经常采用这种方式，尤其在节点数很大时，即使收敛慢的迭代方法也
7、可能比消元法更加有效。在迭代计算中有两个问题，一是迭代的收敛性问题；其次是如何加快迭代速度问题。一般对如导热这一类问题，迭代收敛条件为1Pna（2.10）由于采用有限容积法生成的离散方程，这一条件上述条件一定是满足的。常用的迭代方式中，Jacobi 法的收敛速度最慢，Gauss-Seidel 迭代比较快。交替方向线迭代方法（ADI）是最有利的迭代方法。2.4 边界上不规则区域的处理方法常见的处理方法有阶梯形边界逼近真实边界、坐标变换法、区域扩充法及边界节点单独建立方程法等。在采用商用软件时这些方法通常不需要我们专门去考虑。例题 1：假设图 2-3 所示的矩形截面肋片，肋片根部为温度 T0，肋的
8、上下两侧及端部为对流换热条件。试研究在不同的 Bi（=h / ）数下肋片中的温度分布，并比较数值计算所得出的导热量与按一维假定得出的导热量的区别。图 2-3 例题 1 附图5解：该问题要计算肋片的导热量，实际是肋片表面与周围流体之间的对流换热量，可以在肋片根部取得导热量的数据，所以先要计算出肋片内的温度场。计算条件按照 Bi 数确定，分成 2、1、0.1 和 0.05 等 4 组。假定 L/(2 )=4，并取肋片材料导热系数为 16.27W/(mK)，比热容为 502.48J/(kgK)，密度为 8030kg/m3。L=80mm。具体参数见右表。左边界给定温度 373K，周围流体温度 293K
9、。利用FLUENT 求解。1. 利用 GAMBIT 建立计算几何模型方法见附录。肋片长度为 8 个单位，厚度为 2个单位（建模过程中根据肋片导热上下的对称情况，只画出一半，即 1 个厚度单位） 。2边界条件，对流面根据上述列出的 Bi 数确定对流换热系数。序号 Bi 数 h/W/(m2K) Q/W1 2 3254 21536.42 1 1627 21168.983 0.1 162.7 2402.514 0.05 81.35 2273.7373. 计算结果分析：如果将该问题作为一维模型计算，当端面绝热时，端面温度为305.4K，数值计算得到的端面温度上下各为 306.4K 和 305.8K，平均
10、约为 306.1K，两者相对差别为 0.26%；导热量计算一维模型的为 821.4W，数值计算的为 2400.2W，两者相对差别为 2.6%。所以在 Bi 数为 0.1 时，数值计算结果表明此条件下肋片可以简化为一维模型。Bi=2.0 Bi=1.0图 2-4 不同 Bi 数下肋片内部温度分布Bi=0.1 Bi=0.05Bi=2.0 Bi=1.0图 2-4 中心对称面上温度分布Bi=0.1 Bi=0.05图 2-4 中心对称面上温度分布图 2-5 不同 Bi 数下上下两面温度分布图8从计算结果分析，在 Bi 数小于 .1 条件下，一维模型是对实际的一个近似，但肋片内的温度分布的二维属性仍然存在，
11、即便在 Bi 数为 0.05 的条件下，肋片的上表面和中心对称面上的温度还是有微小的差异。例题 2：一钢锭，大小为 0.50.71.0m，初始温度均匀为 20，其导热系数为40.5W/(mK)，密度 7900kg/m3，比热容 710.05J/(kgK)，试确定将其置入 1200的加热炉中 4 小时后的最低温度和最高温度。假设炉内烟气与钢锭之间的换热系数为 348 W/(m2K)。解：这是三维非稳态导热问题，导热体外表是对流换热条件。可以先估计，最高温度在顶角上，最低温度在中心点上角（见图2-6 中顶角点和中心点） 。计算区域可以按对称性选取为原物体的 1/8，见图2-6 。正面、右面和顶面为
12、对流边界面，背面、左面及底面为对称面（绝热面） 。1利用 GAMBIT 建模：1）确定导热体大小：启动 GAMBIT， R ，弹出对话框如图 2-7 所示。在宽度（Width）栏、深度（Depth）栏、高度（Height）栏中分别输入 0.125、0.175和 0.25；点击 Apply、Close，点击 Fit to Window 按钮查看图形。2）划分网格：依次点击 按钮，打开图2-8 所示对话框。点击 Volumes 右侧的黄色框后，用 Shift+鼠标左键点击工作显示区内几何图形的边线后，在 Spacing 项选择 Interval size，输入 0.005；点击 Apply、Clo
13、se ，所划出的网格为 253550。3）设置边界类型：为方便，先将网格显示关闭，点击，在 Specify Display Attributes 中点击 Mesh 项中的 Off，点击确定，点击 Close；依次点击 ，在 Specify Boundary Types 对话框中加入六个面的边界面属性，分别在Name 项中给定边界面名称、Type、选定边界、点击 Apply。右边界 Right，属性 Wall，在工作区中 Shift+鼠标左键点右边界线变为红色，点击 Apply；各边界面属性分别如下表：名称 Right Left Top Bottom Front Back图 2-6 导热体 图
14、2-7 六面体设置对话框图 2-8 网格设置对话框9Type Wall Symmetry Wall Symmetry Wall Symmetry最后点击 Close，存上 Mesh 文件，退出 GAMBIT。三维建模的步骤比二维要少，原因是建模的开始是以面开始的，而二维要将线转化为面。2启用 FLUENT：按条件分别确定材料、边界、初始等具体数据，下面是时间步长1 秒所获得的数据：10要注意：111）该问题可以通过精确解获得数据，顶角点温度 1470K，中心点温度 1431K；2）因计算时间长，故开始计算时可以考虑两个问题：一是时间步长多少为好？二是空间步长是多少为好？习题：1 将例题 1 实践一次，列出计算报告。2 将例题 2 实践一次，列出计算报告。
