（流体力学专业论文）非结构化网格上非稳态导热问题的求解与应用.pdf

碗士孝愤谚盘非鲢镪化嘲格上非毪卷译执村题的求蛘与鱼曝 摘要 在非结构耽踺格孛，鼗谴求解糕量扩敷方程( 警热方程) 是当今数僮计箨鳇基硪惩 题。由于嚣绩掏他网措对；蕊赠区域鼠育良好的适斑谯，眈及容易实现嘲穰鑫动生成 等优点，在数值计算预域得到了广泛的成用。 本文主器研究非结构化网格中非稳奎导热问题的求解过程，其中包括空闻和时间 物理量熬嵩散方法，三类边器条停熬处理，菇爱不嗣漂瑗豹德醒等。嚣结毒毒耽列辂鹣 生成采用商用网格生成软件导热控制方程中扩散项的离散采用基元有限容积法，二 阶迎风格式 在时问上的离敖，二维闷鼹分别采用了照式格式，c f a n k n i c o l n 烙式 和垒隐格式，三维问题只袋用了全隐格式；并且对格戏靛稳定性进行了奶步讨论。 蘑f o r t r a n 语言缟笃模块纯计篝程序，以茈增强冀通用住。为了验证此求解遵 控的正确性，定整性和通用性，对非稳卷热传导问题本文给出了两类算例，总共十一 个问题。第一类六个闻题同精确解进行了比较，第二爽五个问题同结构化阿格中的数 缓解或文敷审的给的结象避哿魄较，镄舞了比较瑾您戆结鬃。这燕耀为工程计算窝 设计提供有利的数值计算耩j 芋。 关毽词：嚣结构往潮播非稳态等热显式格式c r a n k - n i c o l s o n 捂式 企隐格式 硕士学伊论文 非结构化嘲棒上非稳卷导执问题的求解与戍用 a b s t r a c t i ti so n eo ft h eb a s i cp r o b l e m st os t u d yt h eh e a tc o n d u c t i o ne q u a t i o n ( t h ed i f f u s i o n e q u a t i o n ) u s i n gn u m e r i c a lm e t h o do nu n s t r u c t u r e d 鲥d s i th a sb e e nw i d e l yu s e di na l o to f a r e a s , b e c a u s eo ft h et w om a i na d v a n t a g e so ft h eu n s t r u c t u r e d 鲥d a ，v i z ，t h ee a s eo f 扩破 g e n e r a t i o nf o rc o m p l e x ，r e a l i s t i cg e o m e t r i e s ，a n dt h ea b i l i t yt od y n a m i c a l l ya d a p tt h eg r i d t ol o c a lf e a t u r e so f i n t e r e s l t h ep r o c e d u r eo fs o l v i n gu n s t e a d yh e a tc o n d u c t i o np r o b l e m so i lu n s t r u c t u r e dg r i d a w a si n v e s t i g a t e di nt h i st h e s i s i ti n c l u d e dt h ed i s c r e t i z a t i o nm e t h o do fg o v e r n i n ge q u a t i o n ， t h et r e a t m e n to fb o u n d a r yc o n d i t i o n sa n ds o u r c et e r m s c o m m e r c i a ls o f t w a r ew a su s e dt o g e n e r a t eu n s t r u c t u r e dm e s h e s ，t h ed i f f u s i o nt e r mw a sd i s e r e t i z e db yc e l l - b a s e df m i t a v o l u m em e t h o d , a n dt h es e c o n do r d e ra c c u r a c y t h eu n s t e a d yt e r mw a sd i s c r e t i z e db y e x p l i c i ts c h e m e ，c r a n k n i c o l s o ns c h e m ea n di m p l i c i ts c h e m ei nt w o d i m e n s i o np r o b l e m s a n di m p l i c i ts c h e m ew a su s e di nt h r e ed i m e n s i o np r o b l e m s t h ef e a t u r eo ft h es c h e m e s w a ss t u d i e d c o m p u t a t i o n a lp r o g r a mw a sm a d eu s i n gv i s u a lf o r t r a n ，a n dt h ec o d e sw e r e d i v i d e di n t od i f f e r e n tm o d u l e sa c c o r d i n gt od i f f c r e n tp u r p o s e s t h ee l e v e ne x a m p l e so f t w o g r o u p sw e r ep r e s e n t e dt ov a l i d a t et h ec o r r e c t i o n , i n t e g r i t ya n du n i v e r s a l i t yo fm ep r o g r a m t h en u m e r i c a ls o l u t i o n so f s i xe x a m p l e si nt h ef i r s tg r o u pw e r ec o m p a r e dw i t ht h ea n a l y t i c s o l u t i o n s ，a n dt h e r ea r ef i v ee x a m p l e si nt h es e c o n dg r o u p s ，a n dt h er e s u l t sw e r ea n a l y z e d a n dc o m p a r e dw i t ht h en u m e r i c a ls o l u t i o n sc o m p u t e do ns t r u c t u r e d 鲥d so rt h er e s u l t s s h o w ni nt h er e f e r e n c e a l lo ft h e s ep r o v et h a tt h em e t h o da n dt h ep r o g r a ma r ec o r r e c ta n d p e r f e c t i tw i l ls u p p l yu s e f i dn u m e r i c a lc o d e sa sr e f e r e n c ef o re n g i n e e r i n ga n dd e s i g n i n g k e y w o r d s ：u n s t r u c t u r e dg r i d ，u n s t e a d yh e a tt r a n s f e r , e x p l i c i ts c h e m e ， c r a n k n i c o l s o ns c h e m e ，i m p l i c i ts c h e m e 声明 奉带使论文慧我在导师的指导下取褥的研究成聚尽我所辩，在 本学位论文中，除了加以标注和致谢的部分外，不包岔其他人已经发 表或公布过豹研究成果，也不露会我为获褥饪何教育飘构豹学位或学 历而使用过的材料。与我一同工作的同事对本学位论文做出的贡献均 已在论文中捧了明确戆滢明。 磅究生笈名：壹童莲2 0 0 6 年g 月够医 学位论文使用授权声明 南京理工大学有权保存本学位论文的电子和纸质文档，可以借阅 或上潮公布本学位逾文的全部或帮分内容，可以向奄关部门或艉构送 交并授权其保存、借阅或上嘲公布本学位论文的全黼戏部分内容。对 于保密论文，按缣密的有关规定秘程序处理。 研究生签名：j 蠢缝2 0 0 6 年f 月妫目 南京理工大学硕士学位论文 非结构化网格e 稳态导执问题的求解与廊用 1 绪论 1 1 课题应用背景 自然界和工程上很多导热过程是非稳态的，即温度场是随时间而变化的。例如， 室外空气温度和太阳辐射的周期性变化所引起的房屋维护结构( 墙壁、屋顶等) 温度 场随时间的变化；采暖设备间歇供暖时引起墙内温度随时间的变化，锅炉炉墙和冰箱 外壁中的热传递，地下热网的热损失、核反应堆内热源的放热、电子元器件的散热等 等，这些都是非稳态导热过程。 许多工程实际问题需要确定物体内部的温度场随时间的变化，或确定其内部温度 达到某一限定值所需的时间。例如：机器在启动、停机及变动工况时，急剧的温度变 化会使部件因热应力而破坏，因此需要确定物体内部的瞬时温度场。钢制工件的热处 理是一个典型的非稳态导热过程，掌握工件中温度变化的速率是控制工件热处理质量 的重要因素。再例如，金属在加热炉内加热时，需要确定它在加热炉内停留的时间， 以保证达到规定的中心温度。 由此可见非稳态导热是个实际意义重大的课题。非稳态导热的特点是温度分布随 空间坐标和时问坐标变化，实际应用中我们需要找出不同时刻的温度分布情况，以此 确定热流密度随时间空间的变化等。此外，工程流动与换热过程中的不少现象，其控 制方程与导热方程属于同一类型，例如二维势位流动，常物性流体在平直管道内充分 发展的对流换热、质扩散过程及某些通过多介质的流动和轴承的润滑流动等。导热问 题的计算方法同样适用于上述问题的求解。 1 2 非稳态导热问题的研究方法 工程上计算非稳态导热过程中温度分布问题的方法一般有三种：一是分析法；二 是实验方法；三是数值解法。分析法即按f o u r i e r 的热传导方程求解，能够获得所研 究问题的精确解，可以为实验和数值计算提供比较依据。但是，该方法局限性很大， 即使是求解简单的导热问题也是相当复杂和繁难的，对于复杂几何形状的物体和非线 性边界条件下的导热问题，应用分析解法几乎是不可能的；实验方法是传热学的基本 求解方法，通过实验可以得出相对准确的结果，但是此种方法的适应性不好，针对某 个特定的问题总是要进行新的实验，实验研究耗费的人力、物力和财力巨大，并且实 验周期长；数值解法在很大程度上弥补了分析解法和实验方法的缺点，它适应性强， 特别是对于复杂问题更显其优越性“一。 1 南京理工大学硕士学位论文 非结构化埘格上非稳志导执蜘题的求解与应用 数值模拟与实验方法相比有以下优势： ( 1 ) 缩短设计和研发的准备阶段。数值模拟对不同的结构、试验条件在计算机上 运行只需要较短的时间； ( 2 ) 可以很容易得模拟实验中难以做到的甚至不可能做到的条件； ( 3 ) 数值模拟通过对偏微分方程的简化求解，可以给出温度场的详细信息； ( 4 ) 更能节约成本和降低能量消耗。其数值模拟的计算成本随着计算机硬件的发 展而迅速下降。 随着计算机科学的飞速发展，对数值计算方法的研究有了长足得进展，并得到了 日益广泛的应用。结构化网格中导热问题的数值计算方法已比较成熟，非结构化网格 导热问题也有了长足的进展。数值解法的实质在于使一个连续过程离散化，用一系列 代数方程式代替微分方程式，通过对这一系列代数方程式的四则运算来获得温度场的 近似数值解，可以使许多复杂的导热问题得到铰满意的数值解。 目前，有限容积法是传热数值计算中应用较为广泛的离散化方法，它以控制容积 中心节点作为控制容积的代表，并着眼于控制容积的积分熟平衡。采用结构化网格有 限差分方法和有限容积方法，数值模拟各种导热问题，已经非常成熟，但是，在应用 于复杂区域的求解上，结构化网格有很大的局限性。非结构化网格在有限元法中已经 被采用了多年，然而，它真正应用于计算流体力学( c f d ) 领域中则是在二十世纪8 0 年代，最初是由b a k e r ，w e a t h e r il i 和l o h n e r 将此引入有限容积法的。 除了对复杂区域的良好适应性这一显著优点之外，应用非结构化网格进行计算求 解还具有结构化网格无法企及的其他的优点： ( 1 ) 组合和分解简单，当需要添加和删除某些区域网格时，只需在原来的网格系 统上添加或删除局部网格单元，而无须在整个区域上重新生成网格； ( 2 ) 利用非结构化网格求解时。坐标系和维数的变化对解算器的影响很小。因为 在非结构网格上偏微分方程的离散不依赖于坐标系和维数，只与网格单元自身有关。 自二十世纪8 0 年代以后的2 0 多年中，由于对不规则区域的良好适应性，及其容 易实现网格自动化生成等特点，非结构化网格在计算流体力学、计算传热学等很多领 域中得到了迅速的发展n ，。 1 3 国内外的发展状况及课题的提出 目前，数值计算方法在计算流体力学、计算传热学领域中的应用已经相当成熟， 二十世纪八十年代美国学者s u h a sv p a t a n k a r 教授发展了一套比较成熟的基于结构 化网格的导热问题的计算程序。随着数值计算方法的迅速发展，国外已经开发了许 多基于非结构网格应用于此领域的大型商业通用软件，这些软件的开发和应用，对实 2 南京珲t 大学碗十学付论文非结构化刚格卜非稳态导执问题的求解与府用 际工程问题的研究、设计、优化、管理以及评价都是重要的工具和手段，起到十分有 益的作用。应用于计算流动和传热的比较著名的几个大型商业软件有：f l u e n t 、 p h o e n i c s 、s t a r - c d 和c f x 等都采用了非结构化网格有限容积法。 国内在这方面的研究也取得了大量的成果和进展。西安交大的陶文铨教授、清华 大学的周力行教授等都在基本理论和实际问题的应用上取得了许多重大成果”；此 外，又如华中理工大学的王昌凌、王弘将实体造型理论与导热理论有机结合，实现了 导热三维温度场的动态可视化，为热传导问题的数值仿真提供了有效途径”；华南理 工大学的司光树等”1 以一维非稳态导热为例，研究了步长和差分格式对数值计算的影 响；甄永杰用三角形网格以显式方法对肋片和加热元件等的非稳态导热过程进行了数 值模拟旧”1 ；中南大学土木工程学院对铜筋混凝土结合梁试件的降温过程进行了数值 计算“”；南京工业大学对加气混凝土墙体非稳态传热的数值模拟“”与实验结果吻合较 好等等。这些研究主要涉及到建筑物维护结构“1 、土壤源热泵“”、半导体电极“”、 热电堆例、高炉冷却壁”】、隔热层材料m 矧等等。 但是，据笔者对所掌握文献的调查和研究，国内研究者应用数值计算方法解决工 程实际中的各类导热问题时，通常采用的还是传统的结构化网格有限差分方法和有限 容积方法，有限元方法也有应用嘲，而对非结构化网格应用于非稳态导热问题温度场 求解方面的研究和探讨还相当少，而且尚没有适用于该领域的通用性较强的数值求解 程序，这是与非结构化网格生成技术比较复杂是有很大关系的，非结构化网格中由于 一个节点与其邻点的关系不是固定不变的，因此这种联结信息必须对每一个节点显式 的确定下来并加存储，而在结构化网格中这种联结关系是固有地隐含着的。因此，非 结构化网格对网格信息的存储量较结构化网格大得多，生成技术非常复杂。鉴于这一 点，本文作者试图在非结构化网格在非稳态导热问题中的应用做一些研究和探讨，编 写通用性强的计算程序，为工程应用提供参考。 1 4 本文的主要研究内容 本课题主要研究数值方法求解非稳态导热的温度场分布问题，将非结构化网格有 限容积法应用于求解非稳态温度场分布问题的各种情况，探讨不同时问离散格式的应 用，开发应用于导热问题的基于非结构化网格有限容积法的计算程序。 非结构化网格的生成采用商用网格生成软件，对网格信息进行处理，编写出适用 性较强的通用计算机程序。并选取若干算例与精确解或可信的数值解进行比较，以验 证计算程序的正确性和通用性，为工程计算和设计提供具有指导意义的实用数值计算 方法和通用程序。其具体内容如下： 一、非结构网格有限容积法在二维非稳态温度场数值求解中的应用，包括： 1 有限容积法建立二维温度场的离散方程； 南京理工大学硕上学位论文非结构化嘲格上非稳态导热问题的求解与应用 2 编写显式，c r a n k - n i c o l s o n 格式和全隐格式三种格式的计算程序： 3 选取算例，采用三种格式计算得到的数值解与精确解比较，分析误差产生 的原因，讨论格式的稳定性。 4 应用于实际问题，进行计算和分析。 二、非结构网格有限容积法在三维非稳态温度场数值求解中的应用，包括： 1 有限容积法建立三维温度场的离散方程； 2 编写计算程序； 3 选取算例，将数值解与精确解傲比较，非结构化网格法和结构网格法求解 非稳态导热问题进行比较，对结果进行分析讨论。 4 应用于实际问题，进行计算和分析。 三、对计算程序进行规范和完善，增强通用性。 4 南京理工大学硕士学位论文非结构化州格e 非稳态导捣问题的求解与应甩 2 网格生成技术 对于传热问题进行数值计算的第一步是生成网格，即要对空间上连续的计算区域 进行剖分，把它划分成许多个子区域，并确定每个子区域中的节点。从总体上来说， 数值计算中采用的网格可以大致分为结构化网格和非结构化网格两大类。本章将就这 两种不同的网格及其生成方法进行简要介绍。 2 1 结构化网格生成技术 结构化网格足数值计算中使用的最传统的网格。一般数值计算中正交与非正交曲 线坐标系中生成的网格都是结构化网格，其特点是每一节点与其相邻节点之问的联结 关系固定不变且隐含在所生成的网格中，因而我们不必专门设置数据去确认节点与邻 点之间的这种联系。 生成适体啦标的方法原则上都是一唑特定的变换，即把物坪空间上的一些不规则 区域变换成为计算空间上的规则区域。目前结构化网格的生成技术已相当成熟，大致 可以分为以下几种：正交曲线坐标系中的常规网格、适体坐标法、对角直角坐标法。 其中，适体坐标法又可分为保角变换法、代数法( 包括边界规范化方法，多面法和无 限插值法) 、微分方程法( 包括椭圆型方程法、双曲型方程法和抛物型方程法) 等。 关于结构化网格各种生成方法的原理和特点，在此不作详细介绍。 2 2 非结构化网格生成技术 非结构网格的概念基于如下假设：四面体是三维空间最简单的形状，任何空间区 域都可以被以四面体为单元的网格所填满，对于平面网格而言，这就是三角形单元。 其生成方法是利用三角形( 二维) 或四面体( 三维) 在定义复杂外形时的灵活性，以 d e l a u n a y 法或前沿推进法为基础。全部采用三角形( 四面体) 来填充二维( 三维) 空间它消除了结构网格中节点的结构性限制，节点和单元的分布可控性好，因而能 较好地处理边界，适用于模拟真实复杂外型。非结构网格生成方法在其生成过程中采 用一定的准则进行优化判断，因而能生成高质量的网格，很容易控制网格的大小和节 点的密度，它采用随机的数据结构有利于进行网格自适应。一旦在边界上指定网格的 分布，在边界之间可以自动生成网格，无需分块或用户的干预，而且不需要在子域之 间传递信息。 非结构化网格由于其对不规则区域的特别适应性，自2 0 世纪8 0 年代以来得到迅 速的发展，在这种网格中单元与节点的编号无固定规则可遵循，因而除了每一个单元 南京理工大学硕士学伊论文非结构化嘲格上非稳态导热问题的求解与应用 及其节点的几何信息必须加于存储外，与该单元相邻的那些单元的编号等也必须作为 联接关系的信息存储起来，使非结构化网格的存储信息量比较大。所谓“非结构化”， 就是在这种网格系统中节点的编号命名并无一定规则，甚至是完全随意的，而且每一 个节点的邻点个数也不一定是固定不变的。 总之，与结构化网格中节点排列有序、每个节点与邻点的关系固定不变的严密结 构相比，非结构化网格表现出一种不规则、无固定结构的特点。但正是这种特点使该 结构对不规则区域具有十分灵活的适应能力。非结构化网格中由于一个节点与其邻点 的关系不是固定不变的，因此这种联结信息必须对每一个节点显式地确定下来并加以 存储。然而，与块结构化网格相比，在非结构化网格中这种联结关系是全域一致的， 这就是说，只要设计一种数据结构就可以用来描述整个计算区域中每一个节点的这种 关系，而在块结构化网格中则必须区分块重叠区以及内部区域，因为位于块重叠区内 各计算节点的数据结构是与块内部节点的数据结构完全不同的。非结构化网格的这一 特点对于刚格的自动生成、自适应处珲及平行计算的实施带来了不少方便。从刚格本 身及与此相关的算法来看，结构化网格与偏微分方程的场理论更密切，而非结构化网 格则更接近于计算几何学及计算机科学的原则。 生成非结构化网格的方法可以大致分为三类： ( 1 ) 前沿推进法( a d v a n c i n gf r o n tm e t h o d ) ； ( 2 ) d e l a u n a y 三角形化方法( d e l a u n a yt r i a n g u l a t i o n ) ： ( 3 ) 其它方法。 就方法本身而言，前沿推进法更具有启发性，而d e l a u n a y 三角形化方法则严格 基于计算几何的基础上。其它的方法包括：集上述两法于一体的一些方法，以及在构 思上完全与上述方法不同的所谓四叉树法或者八叉树法等。 生成非结构化网格的前沿推进法是1 9 8 5 年l 0 提出的。在前沿推进法中，从边 界上的网格点所形成的一系列线段( 前沿) 出发，逐一与区域内部的点形成三角形， 所形成的三角形的新的边加入到前沿的行列，而生成该三角形的出发的边则从前沿行 列消去，如此不断向区域内部推进，直到前沿的行列为空，所生成的三角形覆盖全域 为止。前沿推进法的不同实施方法的主要区别在于： ( 1 ) 如何形成边界上的前沿； ( 2 ) 如何向区域内部加点； ( 3 ) 如何推进前沿使与内点连成三角形。 生成非结构化网格的另一种重要方法是d e l a u n a y 三角形化方法汹1 。该方法是一 种将平面上一组已给定的点联接成三角形的方法，这些所联接而成的三角形具有以下 特点：( 1 ) 所形成的三角形互不重叠；( 2 ) 所形成的三角形可以覆盖整个计算区域；( 3 ) 每一个点均不位于不包含该点的三角形的外接圆内。用d e l a u n a y 三角形化方法所联 6 南京理工大学硕士学竹论丈 接成的各三角形中的最小角对给定的这组点是各联接方案中的最大者，称之为最大最 小角准则，对给定的这组点，用该方法生成的三角形的边长的均匀性是最好的。这种 方法基于几何学的严密规则之上，在2 0 世纪的8 0 9 0 年代得到了迅速的发展潞”， 目前已成为生成非结构化网格的一种重要方法。有关非结构化网格生成的知识参见专 门的网格生成技术专著，在此不作为本文的重点进行讨论。 目前，很多通用的商用c f d 计算软件的前处理软件，如t e c p l o t 、g a m b i t 、i c e m c f d 等，已经可以在很短的时问内方便迅速地生成质量相当高的非结构化网格，大大地为 数值计算提供了便利条件。 由于非结构化网格的生成技术十分复杂，工作量相当繁重，为了节约时间，加快 课题进度，本文采用的非结构化网格使用商用网格生成软件生成。 7 碗士学付论文 非结构化阿格上非稳态导执问题的求解与应用 3 非稳态导热方程的数值解法 在非结构化网格上，控制方程的离散不受坐标系和计算区域形状的影响，离散方 法一般采用有限容积法。有限容积法( f i n i t ev o l u m em e t h o d ) 是一种在数值传热学中 被广泛使用的数值计算方法。它的基本原理是将所要求解的区域划分成一系列控制容 积，每个控制容积都由一个节点做代表。通过将守恒型的控制方程对控制容积做积分 来导出离散方程，用有限容积法导出的离散方程可以保证具有守恒特点，而且离散方 程的物理意义明确。 本章将具体介绍二维情况下采用基元有限容积法离散热传导( 扩散) 方程的情况， 并直接推广应用于三维的情况。 3 1 非结构网格的几何描述m 1 如图3 1 1 所示是典型的非正交非结构网格几何关系示意图。 n 图3 1 1 非正交非结构网格几何关系示意图 两单元中心的连线为f 方向，单元表面的切向为昭方向。根据本单元( p ) 与周 围单元( e ) 的关系，在三角形单元上有3 组f ，_ 1 7 。从图中可见 铲蛩y ；5 警 ， 铲等胪等 ( 3 1 2 ) ! ! 兰兰兰 ! ! 竺兰竺竺! ! 苎堡查! 竺望嬖竺苎竺! 竺! 蟛：瓜历瓦；了 ( 3 1 3 ) ，7 = 瓜i 了瓦f 万 ( 3 1 4 ) 乓= x j + y j = x , i 。i + y ，j 一 ( 3 1 5 ) 以= 抄6 一只) a y = - - ( x b x o ) ( 3 1 6 ) a 。= a 3 + a 0 ( 3 1 7 ) 玎= a j 。毡y ，一碱= 等 ( 3 1 - 8 ) 3 2 非稳态导热控制方程在非结构化网格上的离散9 蚓 在直角坐标系中将稳态的导热( 扩散) 方程写为， d i v ( f g r a d 妒) + s = 0 ( 3 2 1 a ) 加入非稳态项即可得到非稳态的导热控制方程， 伊p 等= d i v ( r g r a d 妒) 钳 ( 3 2 i b ) 方程可简写为， p c p 警- d 心 ( 3 2 2 ) 在扩散方程中莎表示标量物理量，r 为扩散系数，d 表示扩散项，对于导热问题，矿 即为温度nf 为导热系数a ，s 为源项，p 为密度，c ，为比定压热容对于固体和不可 压缩的流体其比定容热容f 。= q = c 。啊以写成式( 3 2 3 ) 的形式 s = s 。+ s 。矿 非稳态项的积分， 一f i , p c ，矽y 鸭等叫) 上标护表示上个时层的计算结果。 在控制容积p 内对扩散项进行积分， ( 3 2 3 ) ( 3 2 4 ) 9 硕士学付论文 非结构化嗍格上非稳态导趣问题的求解与应用 盯( d + s ) d 矿= d ，j + d ，+ ( s 。+ s ，九) 4 ( 3 2 5 a ) p 。，l 。 ( d + s ) d y = 口。( 谚一啡) + d ，j + ( 足+ s ，砟) 4 ( 3 2 s b ) b ：导学生 ( 3 2 5 。) 幽i a 。a “ 其中，下标p 表示计算单元的中心点，下标j 表示与之相邻单元的中心点，柏表示与 计算单元相邻的单元的数目，圪表示控制单元的体积。q 是总扩散项。它可以表 示为基本扩散项d p 。和二次扩散项d ，。之和( 下标p 和j 分别表示p r i m a r y 和 s e c o n d a r y ) ，见式( 3 2 6 ) 。 d = f v a = d p j + d ，j ( 3 2 6 ) 三者的关系如图3 2 1 。引入二次扩散项之后离散方程相当于二阶精度。对正交网格 二次扩散项为零，对于四边形、六面体等非结构化网格，计算二次扩散项均比较容易， 但对于三角形，四面体或其它多面体部比较难。对控制单元p 主扩散项和二次扩散项 可以写成， 。圹r l 必d a i 蕊a i o 4 1 ( 3 2 7 a ) 驴c 掣筹碱 慨z ，n ， 这里下标p 代表计算单元的中心点，下标e 代表相邻单元的中心点，a “= 砖，色，= 艺。 图3 2 1 扩散项示意图图3 2 2 单元矢量示意图 1 0 硕士学付论文 非结构化嗍格上非稳态导执问题的求解与应用 式( 3 2 7 ) 的另一种表达方式为， _ - + d p j = l ( 亨妒) 删乞孚坐 a i * e s j ( 3 2 8 a ) 驴暑卜矗西。射， z 鼬， d s 。为两单元中心点p 和e 之间的距离。为方便起见，可将式( 3 2 5 ) 写成 引入加权因子0 ，1 ，式( 3 2 1 b ) 的积分式可化为， ( 3 2 9 ) ( 呜) ，鲁( 九一形) = ，x r p ( 妒) + ( 1 一) r ；( 妒) ( 3 2 1 0 a ) 将式( 3 2 9 ) 代入式( 3 2 1 0 a ) 得式( 3 2 i o b ) i ，善州刚，等呐卜i = i 凸 i 兰e 阮+ ( 1 一门彤一( 1 一门晖】+ l 厂兰+ ( 1 一，) 量d ol ( 3 2 1 0 b ) j - ili = i i = l j + 晦”一) s f + ( 1 一脾o ，o 】咿v 嘛) ，譬 将( 3 2 1 0 b ) 进行整理化简为。 曲 口，九= 口。 矾+ ( 1 一厂) ( 一妒；) 】+ b j 。i 铲耻岳嚣 珊j ( 3 2 1 1 a ) ( 3 2 1 t b ) 6 = ，+ ( 1 一) d ，o j + j s c a g , + ( 1 一f ) s 。a v p 9 ”删0 0 0 ，等露 3 。 + ( 1 一门彤n r ，+ ( 肛，) ，二露 a ，= 差占。+ 学一万，4 。 ( 3 2 1 l d ) 0 ) p 矿 p s+ c s ( +d + 、， p 矿 一 谚 ， 8 = )矿( p r 硕士学侍论交 非结构化嘲格t 非稳态导热问题的求解与应用 以上( 3 2 1 l a - 3 2 1 l d ) 四式为非稳态导热方程的两层格式( 既离散方程中仅出现相 邻两时间层上的值) 的一种通用形式，式中上标0 均表示上一时层的值。当f = 0 时， 为显式格式；f = l 时，为全隐格式；f = i 2 时，为c r a n k - n i c o l s o n 格式啪1 。此 外将方程中的非稳态项去掉，令= l 即得到稳态扩散方程的离散形式 月6 绵幺= q 谚+ 6 ( 3 2 1 2 a ) i = l 铲耻岳嚣 础i ： 6 = d j j + & f i 柚 口，= b ，一s p 4 ( 3 2 1 2 d ) 需要指出的是，坐标系是针对结构化网格而言，对非结构化网格，控制方程和离 散方稃都是在笛卡几直角坐标系下的，它适用丁：任何形状的几何体，以上式( 3 2 1 1 ) 与式( 3 2 1 2 ) 可以直接推广应用于三维的情况。 3 3 初始条件和边界条件 初始条件是研究物理现象在过程开始时刻求解变量的空间分布，对于非稳态热传 导问题，初始条件即求解区域内初始的温度场分布情况。其初始温度场可能是均匀的 常数，也可能是空间坐标的函数。 边界条件是在求解区域的边界上所求解的变量或其一阶导数随地点及时间的变 化规律。对于导热问题其边界条件的表达方式可以分为三类： ( a ) 第一类边界条件 已知任何物体边界上面的温度值，即，7 1 。= 易，式中占表示边界面，是温度 在边界面厅上的给定值。对于非稳态导热过程，乃可以不随时间发生变化，即 易= c o n a t ，若边界面上温度随时间而变化，还应给出乃= 厂( f ) 的函数关系。 ( b ) 第二类边界条件 已知任何时刻物体边界面上的热流密度值。可表示为， 鼋l 。2 q 。或一乳= 警 ( 3 3 1 ) 式中口是给定的通过边界面占的热流密度，对于非稳态导热过程。函可以是常数，若 硕士学付论文 非结构化网格上非稳态导执问题的求解与应用 边界面上的热流密度足随时间变化的，还要给出q w = ，( r ) 的函数关系若某一个边 界面占是绝热的，根据傅立叶定律，该边界上温度变化率数值为零，即 一乳= 。 s z ， 如图3 3 1 所示肋片肋基处的边界条件就是x = o 的界面处热流密度值恒定为零。 = ：0 。：一 婴i d ji 。o 图3 3 1 肋片的第二、三类边界条件 ( c ) 第三类边界条件 已知边界面周围流体温度乃和边界面与流体之问的表面传热系数h 。根据牛顿冷 却定律，物体边界面占与流体间的对流换热量可以写为， q = 。一弓) ( 3 3 3 ) 于是，第三类边界条件可以表示为， 一a 矧。= 一。1 。一乃) ( 3 。t ) 对于非稳态导热过程，h 和弓可以不随时间变化，也可以是时间的函数，这时还要给 出它们和时间之间的函数关系。 3 4 离散方程中相关系数的计算方法”蝴 3 4 1 f d s , 的计算方法 如图3 4 1 1 所示，以一维情况为例，介绍控制单元界面扩散系数f 计算方法， 可以假设f 是通过将单元p 和单元e 的厂值进行线性插值得到，即 e = 石r p + ( 1 一f a k ( 3 4 1 ) 其中， 正= 等 。 ( 3 - 4 2 ) 对于导热方程离散变量为温度r ，扩散系数即为导热系数a 。 硕士学付论文 非结掏化嘲椿上非稳态肆担纠题的求解与应用 界面i 上的热流密度为 p 图3 4 i 1 控制单元p 和e 之间的距离 吼= n ( 丁t p - t d 引用热阻的概念，热流密度可表示成 耳一 一a x i l 堕 砟。以 其中足为有效热阻，由式( 3 4 3 ) 和( 3 4 4 ) 得 去鲁+ 钉出l 疋j 其中f ，出代表尸和f 之间的有效热阻，式( 3 4 5 ) 可写为 三：生生 缸纽t r e 缸t r p 则界面上的当量导热系数可按式( 3 4 6 ) 计算。 对丁= 二维的情况计算方法如下： 由式( 3 4 5 ) 可以类推给出二维界面热导率的表达式： 老= 睁钉凼，l r pkj ( 3 4 3 ) ( 3 4 4 ) ( 3 4 5 ) ( 3 4 6 ) 两边同乘以击f 得 c 咄( 矗 = ( 1 + 七 ( 器 ( 3 4 8 ) 由图3 4 1 2 中所示的几何关系，式( 3 4 1 8 ) 中七可以表示为 后：拿：了d z ( 3 4 9 ) 凼l4 一。 两相邻单元之间的界面f ( 见图3 4 1 3 ) 的方程可表示为 1 4 硕士学付论文 非结构化嘲将上非稳态导执问题的求解与应用 a c + b y + c z + d = 0 ( 3 4 1 0 a ) 图3 4 1 2 控制单元p 和e 之日j 的距离示意图图3 4 1 3 单元p 和e 之间的界面f 对于二维的情况界面i 的方程为 a ) c + b y + c = 0 式( 3 4 1 l a ) 中 p 而p y - i ：o i x - - 工2y - y 2 l 一= y 2 一y b = 工i x 2 c = 工2 y l 一工l y 2 从点p ( x o ，只) 到界面i 的距离为 拈爷首 ( 3 4 1 1 a ) ( 3 4 i i b ) ( 3 4 1 2 a ) ( 3 4 1 2 b ) ( 3 4 1 2 c ) ( 3 4 1 3 ) 将( 3 4 1 3 ) 式代入( 3 4 9 ) 式得 拈嬲( 3 4 1 4 )l 瓜，+ 缈p + q 则界面i 上的离散系数按式( 3 4 8 ) 和( 3 4 1 4 ) 来计算 3 4 2 ( v ) 。j 的计算方法 _ 对于区域内部的单元和边界上的单元，方程( 3 2 8 ) 中( v 矿) 一的计算方法是 不同的，此处我们分两部分来阐述。 1 5 竺主兰竺笙兰 ! ! 竺竺些竺兰兰苎堡查! 兰竺璺竺苎竺：! 堕：生 ( a ) 内部单元p ( b ) 边界早兀p 图3 4 2 1 内部单元和边界单元 ( a ) 内部单元 一 如图( 3 4 2 1 a ) 所示的内部单元，界面j 上( v 。的计算式为， ( 寻妒) 删：业唑 ( 3 4 1 5 ) 其中( 寻矿) p 和( 亨妒) 为单元的梯度。对单元只单元的梯度可由离散定理得 c p = 击善如( 3 4 1 f i ) 其中为单元尸的面数。无是面七上妒白勺平均值，可取谚，和破。的算术平均值，对 于图3 4 1 3 中所示的单元p 的界面i 处的矿值由式( 3 4 1 7 ) 计算 蛾p ：加+ ( 寻矿) ， p 而 ( 3 4 1 7 ) 且要满足 m i n ( 妒p ，妒e ) 矿1 psm a x ( 加，妒) ( 3 4 1 8 ) 其中( 寻妒) 。为重构梯度，d ；。是连接点p 与界面f 中点的向量，如图3 4 2 2 所示。 向量d ；，只与单元的几何形状有关。引入限制系数( v e n k a t a k r i s h n a n ，1 9 9 3 ) 啪1 ， 单元p 中面i 的重构梯度按式( 3 4 1 9 ) 计算 ( 寻妒) ，= 口p 孑妒) p ( 3 4 1 9 ) 将式( 3 4 1 9 ) 代入式( 3 4 1 7 ) 得 咖- p ：加+ 口p 【寻妒) p ；i ( 3 4 2 0 ) 式( 3 4 2 0 ) 中希妒) ，由式( 3 4 1 f i ) 得到。令坼= j 可得到界面上的驴值 1 6 硕士学付论文 非结构化网格上非稳态导执问题的求解与应用 晚爿= 华 但是，这样并不一定能使谚，满足式( 3 4 1 8 ) 需要引入系数q ，对于界面i 有 其中， 对单元p 有 嘶2 商i l ( 。 m 辽t ， 破p 加 蛾， 加 识。p = 咖p 秽= m a x ( o r ，如) ，秽= m i n ( o e ，九) 盯p =r a i n ( a l l 口i ) ( 3 4 2 2 ) 图3 4 2 2 向量d r l ( b ) 边界单元 图( 3 4 2 1 b ) 给出了典型的边界单元，界面l 与物理边界口重合。对于这种情 况有， ( 亨妒) 删：( 亨妒) p ( 3 4 2 3 ) 对于单元的温度梯度( v 矿) p 可用式( 3 4 2 4 ) 计算， 嘞一嘞嘞一嘞 坚缸生嘶 硕士学付论文非结掏化埘捧上非稳态导执纠题的求解与虚用 c 而p = 击砉缸 界面上的平均矿值由式( 3 4 2 5 ) 计算， 荔= 加 磊= 学 面= 华 ( 3 4 2 4 ) 界面f 上妒值的计算在前一节已经介绍过，不再赘述。式( 3 4 2 5 a ) 中边界上 的靠假定是已知的。对于第一类边界条件，边界上的九值给出之后，无须再添加其 它条件。对于第二，第三类边界条件，边界上的热流密度或边界上的热流密度和矿之 间的关系是已知的，九的计算方法将在下面介绍。 3 4 3 边界上温度值九的确定 对于第一类边界条件，边界上的温度值是已知的，而对于第二类和第三类边界条 件：( a ) 给出热流密度，如式( 3 4 2 6 ) ；( b ) 给出热流密度和矿值的关系( 式3 4 2 7 ) ， 边界e 的温度值是未知的要通过计算得出。以下介绍这两种情况下边界上温度值的计 算方法( 图3 4 2 1 b 给出了其中一个边与物理边界重合的边界单元p ) 。 2 9 胂 钆= q 。+ q 。九 式( 3 4 2 6 ) 是式( 3 4 2 7 ) 的一种特殊情况， 的能量平衡关系如下 q 口= d b ( 3 4 2 6 ) ( 3 4 2 7 ) 即= g g f 。，吼= 0 。物理边界上 g 。l j 。l = 。，。+ 。；，口 瓴垤九朴加，麓山 ( 3 4 2 8 ) ( 3 4 2 9 ) ( 3 4 3 0 ) 1 8 硕士学竹论文非结构化阿格上非稳态导趣叫题的求解与应用 其中 九=b 加也川。刚 口一g i - - 4 ) 。 i ( 3 4 3 1 ) 口：孚当坐 ( 3 4 3 2 ) 出口j 。白。 “一 得到边界上的温度值九之后，边界上的热通量即可由式( 3 4 3 3 ) 计算 铲南一以 t s 。， 3 5 稳定性条件讨论1 1 5 1 呻1 式( 3 2 1 1 ) 中，当取= o 时，得到非结构化网格上非稳态导热方程的显式格式 表达式为， 嘶= q 簖一晖) + 6 ( 3 5 1 a ) 铲弘岳等 慨5 埘 6 = 兰或+ 幢+ 鄙睇+ ( p c ，) ，形f 迭名( 3 5 1 。) 口。：! 坐丝坠( 3 5 1 d )口o = t j a 1 0 ， 由式( 3 5 1 a ) 到式( 3 5 1 d ) 四式可得， 办；鍪n b 羞。竺：罡m ：竺：竺：竺竺竺 s z ，。 啦0p p a i 1 9 硕士学付论史非结构化网格上非稳态导执问题的求解与应用 卟善丢等( p c p ) p a 。i a 妒 萋象象+ 萋n b 0 弛叫删 。5 _ 当内热源不随时间变化时，式( 3 5 3 ) 中右边第二项中秽的系数是不小于零的，第 一项中蝶的系数必须满足， 卜器筹犯驴。 慨， 也即， 善暑著7 慨例 c s s s ， 式( 3 5 5 ) 即为采用显式格式时，为加快计算进程而调整时问步长f 和网格单 元尺寸时应满足的关系。因为当该式小于零时，织的系数为负数。这样不同时刻的 计算值会出现波动，导致出现违反热力学第二定律的结论“”。因此在应用显式格式 进行计算时