（流体力学专业论文）二维非结构网格Euler方程的LUSGS算法.pdf

摘要 本文对将二维非结构三角形网格重新排序之后、用l u s g s ( l o w e r u p p e r s y m m e t r i cg a u s s s e i d e l ) 隐式方法计算e u l e r 方程的方法进行了尝试。在l u 近似因 式分解的过程中，本文采用了一种新的推导方式。陔方法与传统方法在结构网格中 结果一致，但能更广泛地用于重排后的非结构网格中。文中提出了一种适用于格心 格式计算的网格重排方法，且经重排的网格能达到矢量化计算的目的。三角形网格 重排后出现了“不平衡”问题( 即某一网格的邻居中编号比该网格大的网格数与编 号比该网格小的网格数为1 ：2 或者2 ：1 ) ，本文对此情况下的计算可行性及计算效果 进行了探讨。，皋文对不同翼型在不同流动条件下的二维e u l e r 方程进行了计算，并 和四步r u n g e k u t t a 显式方法进行了结果比较，验证了此方法的良好的计算、收敛 效果；同时通过物面附近网格加密后的实例计算，发现了该方法在粘性计算方面将 有广泛用途。该方法克服了以往的隐式方法大量耗费内存的弱点，达到了计算费时 短和占用内存少的统一。 关键词：网格重新排序，l u s g s 隐式方法，笤笼运亟亘垂丕蕉垂蠹e 五五瑶蠢澎7 赠三堡翼等翥嚣委夏k _ a 籼博褥渤窆； f 司榕连叙j 一，1 抖y 莉黝； a b s t r a c t am e t h o df o r2 de u l e r e q u a t i o nc o m p u t a t i o n w h i c h u s e s i m p l i c i t l u - s g s a p p r o x i m a t ef a c t o r i z a t i o ns c h e m ea f t e rr e o r d e r i n go f 2 - du n s t r u c t u r e dt r i a n g l eg n d s ，i s d i s c u s s e d an e w w a y f o rl u - s g s a p p r o x i m a t ef a c t o r i z a t i o nd e d u c i n gi sp r o p o s e d t h e n e ww a yc a ng e tt h es a m er e s u l ta st h eo r i g i n a lw a yo ns t r u c t u r e dg r i d s ，b u tc a nb e w i d e - s c a l e a p p l i e d o nu n s t r u c t u r e d g r i d s w h i l et h e o r i g i n a l o n ec a l l t t h e g r i d s r e o r d e r i n gb r o u g h to u th e r e i sf o rc e l l - c e n t e r e df i n i t e - v o l u m es c h e m ea n di tl c a d st o v e c t o r i a lc a l c u l a t i o n t h er e o r d e r i n gf o r2 - dt r i a n g l eg r i d sc a u s e si m b a l a n c eo f 酣d s ( i t m e a n st h a tf o rg r i di ，t h er a t i oo ft h es u mo fi t sn e i g h b o r sw h o s en u m b e re x c e e di ta n d t h a tl e s st h a ni ti s1 ：2o r2 ：1 1t h ep o s s i b i l i t ya n dt h ee m c i e n c yf o rc o m p u t a t i o nu n d e r t h i sc o n d i t i o na r ed i s c u s s e d t h em e t h o di sa p p l i e dt os e v e r a l2 df l o wp r o b l e m sa n dt h e r e s u l t sa r ec o m p a r e dw i t ht h o s ef r o mt h e4 - s t e pe x p l i c i tr u n g e k u t t as c h e m e t h e c o m p a r i s o nd e m o n s t r a t e s t h em e t h o d se f f i c i e n c ya n d g o o dc o n v e r g e n c e t h ep o t e n t i a l i t y o ft 1 1 i sm e t h o di nv i s c o u sc a l c u l a t i o ni sp r o p o s e d u n l i k em a n yo t h e rt r a d i t i o n a li m p l i c i t m e t h o d s i tr e q u i r e sm u c h1 e s sm e m o r ya n da l s oc o m p u t e sf a s t k e yw o r d s ：g i d sr e o r d e r i n g ，i m p l i c i t l u s g ss c h e m e ，l u s g s a p p r o x i m a t e f a c t o r i z a t i o n d e d u c i n g ，c e l l - c e n t e r e ds c h e m e ，i m b a l a n c e o f g r i d s ，4 - s t e pe x p l i c i t r u n g e - k u t t as c h e m e ，p o t e n t i a l i t yi nv i s c o u sc a l c u l a t i o n 2 塑塞堕奎堕丕盔兰堡圭兰垡笙塞一 第章绪论 “网格生成”和“数值计算”是支撑起计算流体力学( c f d ) 的两块具有举足轻重作 用的巨大磐石。“网格生成”是“源”，没有了它，计算流体力学便成了无源之水；而网 格生成技术的长江后浪推前浪，为c f d 的发展起到了关键而显赫的作用。随着飞行器外 形从简单而日趋复杂，随着人们在该领域的研究的步步深入，网格生成的方法也随之在 调整、完善和更新。追溯“网格生成”方法历经的足迹，我们可以看到它发展的过程： 2 0 世纪7 0 年代，t h o m p s o n 等提出的采用求解椭圆型方程方法生成单域的贴体网格为网 格生成技术开创了先河。紧随其后，s t e g e r 等又提出采用求解双曲型方程方法生成贴体 网格。这样的网格能较准确地满足边界条件，适合于较为简单的外形。当外形变得逐渐 复杂时，生成单域贴体网格就有些勉为其难了，即使牵强生成了，网格质量也无法得到 保证，会影响到在整个流场中数值求解的效果这在8 0 年代曾经圜扰了人们一段时 间。但很快人们就研究出了应对的策略，即采用分区网格：根据外形的特点将总体流场 拆分成若于个子区域，在各个区域内分别建立网格并在交界处形成交界面网格。分区网 格又分为各子域无重叠部分的“分区对接网格”和有重叠部分的“分区重叠网格”。尽箭 分区网格已是一种成熟的构造复杂外形网格的方法，但旦要对外形做局部修改或改变 其构形时，重新分区和构造网格就要耗费不少的时间和精力了。于是在9 0 年代时，人们 的注意力逐步转向了非结构网格与笛卡尔自适应网格。与结构网格相比，使用非结构网 格能大大减少网格的生成总数，同时非结构网格在对复杂外形的适应性方面也具有相当 大的优势；而笛卡尔自适应网格的优势则体现在网格生成过程简单、省时，且网格容易 加密，有利于提高计算精度。这两种方法的迅速发展使具有复杂外形的网格生成技术呈 现了一片欣欣向荣的局面。“网格生成是c f d 作为空气动力学工程应用的有效工具所面 临的关键技术之一” 1 ，沿着网格生成技术的足迹一路走来，我们看到了它在不断的变 化中“汰弱留强”、推陈出新的适应过程。 在为具有复杂外形的几何体“度身定做”结构或非结构网格的这些方面我们已作出 了很多的努力和尝试，得到了不少成果。而在如何进行数值计算以及如何加速计算的收 敛性方面，我们也付出了不懈的努力，取得的发展也是有目共睹。早在计算机诞牛之前， 就有为数不少的“先行者”们在计算方法上研究出了许多有意义的成果 2 。这其中很多 人的名字是我们耳熟能详的，他们创立的一些理论和方法也为我们熟知和相传：在2 0 世纪1 0 年代，r i c h a r d s o n 提出了解l a p l a c e 方程的点迭代格式及松弛技术，随后又有 l i e b m a n n ，s o u t h w e l l 等改进了松弛法，使之被广泛应用于求解流体和固体力学问题： l e w y 、f r i e d r i c h s 和c o u r a n t 提出了特征线方法，给出了著名的稳定性判别条件c f l 条件。这些学者的研究成果和其他一些数学家研究的偏微分方程数学理论相结合，构成 了有限差分的数学理论基础。在此之后，v o nn e u m a n n ，r i c h t m y e r ，l a x 等研究并建立 三丝韭堕塑旦整呈坐竺查篓塑兰竺：! 曼! 墅鲨一 了非线性双曲型方程守恒律的数值方法理论，尤其以其弱解理论为含有激波及其他间断 的气体流动的数值模拟打下了理论基础；v o nn e u m a n n 还提出了一种时间推进问题数值 方法线性稳定性的分析方法，l a x 和k r e i s s 等也给出了非定常偏微分方程差分逼近的稳 定性理论；v o nn e u m a n n 和r i c h t m y e r 明确提出了采用人工粘性的方法作为带激波计算 的流动，这也成为了“激波捕捉”法的起源ff r a n h a l 提出了解l a p l a c e 方程的逐次超 松弛( s o r ) 法，大大提高了收敛速度；而p e a c e m a n 和r a c h f o r d ，d o u g l a s 等一起采用 交替变换迭代方向和不受限制的时间步长，把无条件稳定性和计算简便性结合了起来， 形成了隐式( a d i ) 格式。 计算机的问世和其技术的迅速发展也为计算流体力学起到了推波助澜的作用：6 0 年代面元法的应用；7 0 年代的无粘非线性方程数值求解的突破；以及7 0 年代以来在采 用时间相关方法求解可压缩e u l e r 和n - s 方程技术之模拟复杂流场等研究领域的突破性 进展。经过长时间的努力之后，c f d 工作者们发展了相当数量的高精度、高分辨率的差 分格式，如总差分减小( t v d ) 格式，基本无振荡( e n o ) 格式，无波动、无自由参数的 耗散( n n d ) 格式，耗散比拟方法，界值为限( b ) 格式及矢通量分裂和通量差分裂格 式等。而现在，人们正在致力于更高精度的计算格式和方法等方面的研究。 在进行数值计算时，如何加快计算的收敛一直是我们关注的问题。由于结构网格的 有序性，很多行之有效的加速收敛的方法已经很好地运用于结构网格上，如：近似因式 分解迭代算法( 简称a f 算法) ，网格逐次加密法，多重网格法等。当非结构网格的运用 越来越广泛时，我们非常希望它也能“分享”结构网格的这些成果，从而可以节约不少 “白手起家”、重新探求适应予非结构网格的加速收敛方法的时间。这样一个美好的心愿 曾经处处碰壁，因为我们通过不断的研究而得的、能加速在结构网格上的计算收敛的方 法往往依赖于网格节点的有序性，这对舍去了网格节点结构性限制的非结构网格就不适 用了。这曾经让非结构网格加速收敛的步伐一度停滞不前，也带来了不少问题：如非结 构网格的计算不得不是显式的，运用于非结构网格的隐式计算代价昂贵且不能得出令人 满意的结果( 与结构网格带给我们的相比) 等。目前，已有一些隐式的方法运用在的非 结构网格上来加速收敛c 3 、 4 ，如带有预处理的刷r e s 法和c g s 法。但是，这些隐式 方法却由于其非结构网格的任意性而大大增加了记录每个单元信息的存储量。虽然种 无需建立显式矩阵的n e w t o n - k r y l o v 5 方法也引起了人们的注意，但是在它的预处理过 程中仍需要计算大量的矩阵。 在非结构网格中采用了隐式方法却付出了增加不少的内存存储量的“代价”，加速 收敛和减少存储量似乎成了矛盾体。所以，当在结构网格计算中行之有效的l u - s g s 近似 因式分解法被成功“挪用”到了非结构网格中 6 j 时，让我们为之振奋因为，而且通 过理论与实例我们已经验证，这样做不仅不会占用额外的存储空间，而且它的“性能” 也能和它在结构网格中的运用相媲美。 值得我们注意的是，l u - s g s 方法的运用很关键的一点就是：要能在有序的网格节点 南京航空航天大学硕士学位论文 上“扫过”从编号最小的到编号最大的，或反向。对于这，结构网格有得天独厚的 优势：它的各节点是有序的；而非结构网格则不然，在其他方面给我们带来了很多便利 的各节点的随机性，在这儿就有些一筹奠展了。不过“车到山前必有路”，本文就要介绍 一种“网格重排”的方法，重排之后的各网格由无序变为有序，l u s g s 方法便又有了用 武之地。网格经重排之后，显示了很好的计算收敛性。在参考 6 中提及了这种网格重排 的方式，但它是为非结构网格的格点格式设计的；本文中提及的这种重排方法却能用于 非结构网格的格心格式，并已在二维非结构三角形网格上作了尝试。 在网格重排的过程中，为了避免过多地出现某一个网格周围的“邻居”编号都比其 大或都比其小的情况，网格重排的起点不再任意选取，而是从物面边界的整个一圈开始， 逐层向外推进。采用网格重排的初衷是为了使网格变得有序，能采用l u - s g s 法进行计算。 网格重排后的关键是要使一个网格的周围邻居中，编号既有比该网格大的，又有比该网 格小的，而且在邻居中，这两种网格的数量越接近则越具有平衡性。但由于现在网格重 排是用于二维非结构三角形网格的，所以对每个网格周围的邻居进行编号时，在绝大多 数情况下，必定出现的是编号比该网格大的网格数量和编号比该网格小的网格数量之比 为i ：2 或者2 ：1 ，这样就打破了两种情况下的网格数应接近相等的“平衡性”。网格重排 后出现了这种不平衡的情况，用l u - s g s 法进行计算的效果如何? 本文作者进行了尝试。 结果证明在不平衡情况之下，用l u - s g s 法进行计算仍能取得较好的计算效果，而且其所 占内存少、计算费时少等优势已经初见端倪。因此推测当网格重排后的“平衡问题”能 更好地处理之后如在二维情况下对非三角形网格进行重捧、采用格点格式对网格点 进行重排，或对三维非结构网格进行重排该方法在计算收敛速度快方面的优势将有 更突出的体现。 本文在推导l u 近似因式分解时还采用了一种不同于传统方法的新的推导方法，文中 将新的推导过程和传统的推导过程进行了对比。我们可以看到，在新的推导过程中，虽 然站在了和传统方法相同的推导起点上，可是采用的却是与传统方法不同的思考角度， 而最终结果是殊途而同归。通过对两种推导方式的比较我们可以看到，传统的l u 近似因 式分解的推导是基于结构网格的，无法推广到非结构网格；而新的推导方式却是基于非 结构网格的，当它从非结构网格“退化”到结构网格时，与传统方法的推导结果完全一 致。因而新的推导方式具有更好的普适性。 本文将“网格重排后使用l u - s g s 算法”的方法运用在了二维e u l e r 方程计算中，并 进行了丰富的实例计算。如0 0 1 2 翼型在亚音速、超音速情况下的计算，两段翼型在亚音 速情况下的计算，上下交错两段翼型在亚音速下的计算，单段翼型物面网格加密后( 粘 性网格) 在亚音速情况下的计算等等。计算结果和显式四步“r u n g e k u t t a 法 7 的计算 结果进行了比较，发现这种“网格重排后使用l u - s g s 算法”的方法在各种情况的计算中 都体现了较好的解的精确性，并在绝大多数情况下都要比r u n g e k u t t a 法得收敛速度要 明显地快。值得一提的是，在单段翼型物面网格加密后( 粘性网格) 的计算结果对比中， 二维非结构网格e u l e r 方程的l u s o s 算法 我们可以看到运用这种方法能较快地收敛，而“r u n g e k u t t a ”法几乎无法收敛，这就体 现了该方法在粘性计算中具有的巨大潜力。 第二章网格的生成方法 网格生成技术是计算流体力学发展的一个重要分支。为了较为容易地满足边界条 件，传统的方法是形成贴体( 单域) 计算网格；而随着外形复杂程度的提高，要形成单 域的贴体网格变得较为困难了，因此人们又发展了新的分区结构网格和非结构网格的方 法。下面就生成结构网格的单域( 贴体) 网格的代数生成法、分区对接网格的生成法以 及非结构网格的d e l a u n a y 生成法和推进阵面法 2 向读者做一个介绍。 2 1 结构网格的生成方法 2 1 1 单域结构( 贴体) 网格的生成方法代数生成方法 通常来说，二维流场的基本做法是将原来的物理流动平面( x ，y ) 变换到计算平面 ( ，n ) ，并将原来的求解方程变换为( e ，q ) 平面上的方程，并在( ；，n ) 面中 的均匀网格中求解此方程，为此必须知道相应的导数关系t ，i ，n ，n ，。网格生成采 用某种方法实现这种变换关系，并找出相应的导数关系。单域结构( 贴体) 网格的生成 方法大致可以分为三类：代数生成方法，椭圆型微分方程生成方法，双曲型微分方程生 成方法。以下要介绍的是二维情况下的代数生成方法，其方法可以推广至二维情况。 利用已知的边界值采用中间插值方式来生成网格系的方法称为代数生成方法。给出 一个包含某些待定系数的特殊曲线坐标函数f ( 掌，7 7 )( 即插值函数) ，根据拟合边界上或 某些中间点上给出的笛卡尔坐标值，以及可能还有的导数值要求来确定这些待定系数。 如此一来，在曲线坐标的常值处即可用此插值函数来确定笛卡尔坐标值，从而确定整个 网格系统。由此看来，插值函数是代数方法的核心。一维插值函数可以写为： 上， f ( 善) = ：妒( 习东 ( 1 ) = i 其中帆函数是满足 畎分 的多项式，称为型函数( “b l e n d i n g ”f u n c t i o n ) 。若采用l a g r a n g e 多项式，则 ( 2 ) 堕塞塾窒塾蒌查兰堡主鲎垡笙苎。一 州争：卉鍪洲圳 ( 3 ) 仇( 书= 兀 兰；( ，h ) ( 3 ) 采用不同形式的插值函数，即可形成不同的代数方法。多维时需要多方向的插值，此处 仅介绍“超限插值”( t r a n s f i n i t ei n t e r p o l a t i o n ) 的思想 8 ，即二维时在每一个方向 单独采用多项式插值： 芦( 孝，7 ) = 纯( 号) i ( 孝。，柙) 肿】 1 ( 4 ) 村 一 i ( 乎，叩) = 矿。( 导) 尹( 善，卵。) 两式相加，则右端项为： 蜃( 翻) ：兰吼( 雩) 尹( 朋+ 兰( 导) 珊，巩) ( 5 ) 在e = 0 的边界处， 蜃( o ，叩) = 讯叩) + ( 争酏巩) ( 6 ) 这表明简单地在各方向单独采用多项式插值不能在边界处拟合给定值。为此，插值 函数应取为： i ( 孝，可) ：窆吼( 争i ( 知即) + 兰y 。( 旦) 尹皤，刁。) 一芝兰纯( 萼) 矿。( 号) f ( 乞，) ( 7 )i ( 孝，可) = 吼( 争i ( 孝。，即) + y 。( 号) 尹皤，刁。) 一妒。( 号) 矿。( 号) f ( 乞，) ( 7 ) n - i - l j 抽ia t - i o 这样就能在各边界处都拟合给定值。上式右端最后一项称为插值的张量一乘积形式。上 式也可以改写成 7(4，rd：m7(4 r d 妒，( 旦) 尹( 孝，) + 兰吼( 萼) 湫厶 叩) 一兰( 号) i ( 掌，r ，) j ( 8 ) ，= 妒，( 旦) 尹( 孝，) + 吼( 号) f i ( 厶，7 ) 一( 号) i ( 掌，) j ( 8 ) 具体计算时： 丘( 掌，7 ) = y 。( 卫) i ( 孝，) ；n 方向的插值 ( 9 ) 元( 乎，叮) = 仉( 争( 尹( 六，智) 一丘( 磊，玎) 】；方向的插值 ( 1 0 ) f ( 孝，i ) = e ( 孝，r ) + e ( 掌，r ) ( 1 1 ) 考虑到只给定两个边界面的函数较难控制所需要的网格的分布，且在壁面附近网格线的 正交性对流场的计算结果有较大的影响，因此可以采用广义的超限插值 9 ，即不仅给出 函数本身，还给出外法线方向的导数，从而有效地控制在壁面附近插值函数的特点。 三丝j ! 堕塑塑竖里型! ! 互堡箜! 望：! ! 望竺堕一 2 1 2 分区对接网格 对于复杂多部件的飞行器外形，要生成单域( 贴体) 网格就有相当的困难了。因此 目前采用分区网格及分区计算技术，即根据外形的特点将总流场分为若干个了域，对每 个子域分别建立网格并在其中对流动主控方程求解，各子域的解在相邻子域边界处通 过耦合条件来实现光滑。区域分解的基本原则是：尽量使每个子域的边界简单以便于网 格的建立；各子域大小也尽量相同以实现计算负载的平衡。 实际飞行器的外形通常是以离散数据给出的，常常不能满足c f d 计算的要求，因此 在开始生成网格之前，必须进行几何处理，所定义的外形曲线要与原始数据所描述的吻 合一致。然后在几何处理的基础上生成计算网格。图l 表示了网格生成的过程。 分区网格生成的步骤通常可归结为； 也引至h 崦 画 图 l 1 根据外形和流动特点分区，并确定每区中的网格拓扑； 2 生成表面网格在上述几何处理基础上，按要求的网格疏密分布生成各部件表面 的网格。现取一个机翼机身尾翼的全机外形为例，简单介绍生成表面网格的过程 1 0 。 若机身表面形状可以解析表达( 如旋成体机身) ，则在求出机身与机翼、平尾、腹鳍等部 件的交线后，机身及各部件表面网格的划分都较容易。但实际飞行器机身外形都以离散 的数据方式给定，并且只给出若干个横截面的外形数据，多个截面数据的分布还不完全 对应，这样，机身外形的模拟及其表面网格的划分需要依靠插值来完成。飞行器表面网 格的划分方法还取决于总体网格的拓扑结构。对于沿机身轴线方向取一系列横截面的方 式来划分拓扑时，表面网格的生成按下列步骤进行； ( 1 ) 重新整理机身外形的原始数据；用以弧长为自变量的线性插值对每个机身截面 的原始数据重新进行分布，可以等弧长分布、也可变弧长分布，重新分布的节点数至少 南京航空航天大学硕士学位论文 要大于原始数据中截面分布节点最大数，以保证重新分布的数据所描述的截面外形与原 始外形吻合一致。重新分布形成的新数据在各个截面的节点数是相同的，新数据是机身 表面网格生成的基础数据。 ( 2 ) 求出机翼、平尾、立尾、及腹鳍等部件的表面函数及与机身的交线函数。 ( 3 ) 按网格数及疏密分布的要求沿机身轴向取一系列横截面，用基础数据沿x 方向 作样条函数插值，获得这些截面上与基础数据节点数相同的外形曲线( y ，z ) 坐标的节 点值，以此节点值为基础，按照截面上的网格数及网格点疏密要求，在横截面上用以弧 长为自变量的线性插值求得各网格点的坐标( y ，z ) 值，即所需各截面上机身表面的网 格，或整个机身表面的网格分布。 以弧长为自变量的线性插值原理为： 暾采截回外彤曲线有n 个数话节点( ”，z ，) ，i = 1 ，2 ，no 足义s l 2 0 ，依此求得母一节 点处的弧长坐标 s ，= s 。+ 厄了石再i i 了，i = 2 , 3 ， ( 1 2 ) 对于任弧长值s s ，s ， ，当s 晦，。 时，对应s 的y ，z 值为： p ，+ 丽y i + i y 1 ( n s ) 、 ( 1 3 ) h + 嚣_ 弘踟 ( 4 ) 其余部件先在弦向生成网格，然后在展向插值获得表面的网格点。 3 生成交界面网格空间流场分区后，相邻区之间的公共交界面一般是一个空 间曲面，它在空间的位置、走向及其上的网格点分布极大地影响者以它为边界的两相邻 空间网格的生成过程和结果的质量。在一空间曲面上生成网格或在空间生成一个网格曲 面，原则上讲，可使用的方法是很多的。这里介绍一种矢性三次多项式插值的方法 9 ， 通过在两条边界线上的对应节点之间生成的两格线来构造交界面网格曲面。设a ，1 3 分别 是这两条边界线上的一对对应点，从a 到b 之间的矢性三次多项式插值函数为 芦( 善) = 磊+ 西善+ 磊善2 + 毛毒3 ( 1 4 ) 其中尹是节点的位置矢量，氐，蟊，昂2 ，厅3 是待定参数。假设要在a ，b 之间形成条n 个节 点的网格线，则可取专= ( i i ) ，( 一i ) i = 1 , 2 “3 ，n 。 三丝! ! 丝塑旦鱼里! ! ! ! 查堡盟! 旦：! 鱼! 墨鲨 一 一一 根据i ( o ) = ，i ( 1 ) = 名，i ( o ) = l ，i ( 1 ) = l 可以确定待定参数 口02 o l 2 r 。 ( 15 ) 五：= 3 ( 毛一) 一( 矗+ 2 ) 瓦= 毛+ t 一2 ( 瑶一) 其中无，不代表了从a 点到b 点的网格线在a ，b 点的切向矢量，其模代表了网格在a ，b 点附近的密度。通过调整t ，矗的大小和方向可控制网格在a ，1 3 点的密度和走向，牛成 符合要求的网格。在两条边界线的每一对对应点之间生成一条这样的网格线，这些网格 线的集合便在这两条边界线之间构成了一个网格曲面，这就是所要求的相邻公共交界面 的网格曲面。 当交界面正好是一个平面时，其上的网格可用求解二维椭圆型方程的方法方便地生 成。 4 空间网格的生成当表面和交界面上的网格生成后，各区的边界即已确定，各 区内空间网格原则上可以用代数方法或求解椭圆型方程的方法生成。若用椭圆型方程方 法，可先求解l a p l a c e 方程生成初场网格，再用修正源项方法调整网格线与边界的夹角 及与边界的间距，以改善空间网格的质量。 2 2 非结构网格的生成方法 与结构网格相比，非结构网格舍去了网格节点的结构性限制，易于控制网格单元的 大小、形状和网格点的位置，因此在处理复杂外形问题时有更大的灵活性。此外，对于 结构网格，在计算平面内网格线和平面都应该保持连续，并且正交于物体边界和相邻的 网格线和面；非结构网格没有这方面的约束，这就消除了网格生成过程中的一个主要的 “绊脚石”。而且，非结构网格中个点周围的点数和单元都是不固定的，这样就可以便 于自适应计算，合理分布网格的疏密，提高计算精度。随着需要处理的气动外形的曰渐 复杂，集众多优点于一身的非结构网格越来越“得人心”，而非结构网格技术也得到了迅 速的发展。 目前用得最多的非结构网格单元形状是三角形( 二维) 和四面体( 二维) 。将边界上 与流场中规定的节点连成互不重叠又覆盏全流场的许多三角形或四面体的过程，称为二 角形化。在现有的许多生成非结构网格的方法中，最常用的三角形化的方法有两种，一 是d e l a u n a y 三角化法，另一种是推进阵面法。因为本文中介绍的网格重排是基于非结构 网格之上的，所以下面将介绍一下这两种非结构网格生成方法。 南京航空航天大学颂：t 学位论文 2 2 1 d e l a u n a y 三角化方法 2 。2 1 1 基本原理 平面上的节点可按照参考 2 生成。给出节点之后，下一步就要进行将平面二角化的 过程。 平面化的依据是d i r i c h l e t 提出的一种利用已知点集将平面划分成【一、多边形的理 论。这种理论的基本思想是：假设平面上给出一组点 p 。) ，则 p 。) 中的每一个点( 不失一 般性，取其中的一个点p ；) 都有属于自己的一个区域s ，使得区域s 内的任意一点与该 ( 点) 的距离都比与中其他点的距离近。这种划分方法把平面划分成了一系列不重叠的 凸多边形，叫做v o r o n o i 区域，覆盖整个控制体，如在图2 中，九个点组成的点集 a ) 按照d i r e c h l i t 理论将平面划分成了若干个凸多边形( 其中有的凸多边形顶点在无穷远 处) 。以点5 为例，点5 所拥有的凸多边形v 。v ：，v 。v 。v 6 中的每一点距离点5 都比其余的 八个点近。实现过程难度不大：将点5 和点集( a 中的其它八点连线，做每条线段的垂 直平分线，则这些垂直平分线就会 6 图2 二维非结构网格e u l e r 方程的l u - s o s 算法 自动构成一个凸多边形。凸多边形的每一条边都对应着点集 a 中的两个点，如构成i “i 多 边形v 。v ，v ，v 6v 8 中的边v 3v 4 就对应着点对( 4 ，5 ) ，而边uv e 则对应的是( 5 ，7 ) 。将 所有的点对连线，则整个平面就被三角化了。对给定的这九个点组成的点集来说，这种 三角化是唯一确定的，而且深沉的三角形的形状都比较接近正三角形这样的划分对 流场的离散和求解非常有利，也是d e l a u n a y 方法的优点之一。 我们可以给出一张表格来更清晰地说明格点之间的数据关系。我们先作如下定义： 三角形的顶点被称为“形成点”，“形成点”包括下面将提到的“初始化点”和点集( a j 中的点；定义三角形的垂心( 即凸多边形的顶点) 为“顶点”，“顶点”集记为 v t ) 。 表l 中给出了图2 所对应的数据关系。每个凸多边形的顶点周围都有三个点集 中的点，而每一个顶点周围又都有三个相邻的项点。边界上的顶点对应的相邻顶点在无 穷远处，在表格中记为一1 。这样的数据关系将在下面的三角化的过程中起到重要的作用。 顶点形成点相邻顶点 v 1 24v ，一1一l 也23 5u一1 v 3254v iv 2v 4 v 4457v ；kv ； v s4 76v v ，一l v 658 7 v iv r 一1 v ，679v 。一11 v 8385v 8v 2 1 表1 2 2 1 2 d e l a u n a y 三角化进程 d e l a u n a y 三角化由下面各步骤组成： ( 一) 初始化 设给定的点集为( p 。) 。首先定义一个包含 p ) 的凸多边形外壳，此外可进行初始 三角化，例如可以定义一个四边形，建立起这四个点所组成的v o r o n o i 图结构。图3 中给出了这一典型的初始化图，表2 给出了相应的数据结构。顶点l ，2 ，3 ，4 并未 给出严格的定义，因为他们位于凸多边形外壳的外侧，因此不拥有三个形成点，也 没有三个相邻的顶点。在表格中。我们用一1 来表示缺少的形成点和顶点。 ( 二) 引入新点 凸壳内引入一点pe 只) 。新近引入的这一点将破坏原来的三角化结构，将一些 三角形删除，并在其周围建立新的三角化结构。 f 0 南京航空航天人学硕士学位论文 顶点形成点相邻顶点 v ll24v ：jv bv i v 2234v lv jv 1 v 3l41v 一ll v 443一lv 2一l一1 v 532一lv 21一l v 6211v 一l一1 表2 ( 三) 确定将要被删除的顶点号 每一个顶点都有三个形成点，该顶点又恰好是三个形成点所构成的二角形的外 接圆圆心。如果新近引入的一点p ( 见图3 ) 恰好落入某个顶点所属的外界圆，则该 顶点将被删除。在本例中，新加入点p 落入了三角形i - 2 - 4 ( 顶点vi ) 和2 3 4 ( 顶 点u ) 的外接圆中，因此顶点v 。、v 。将被删除( 见图4 ) 。 ( 四) 确定每一个被删除顶点的形成点 这些形成点将和新引进的点成为“相邻点”，把它们分别与新引进的点连线，构 成新的三角形和新的顶点( 也是每个三角形的外心) 。在此例中，点i ，2 ，3 ，4 都 是符合条件的形成点。 ( 五) 确定邻接顶点 确定与被删除的顶点相邻而自己又未被删除的顶点。这些顶点在构造新的顶点结 构中将起重要作用。在本例中，v 。，v 。v 。，v 。符合上述条件。 v 图3 图4 鬻 二维非结构网格e u j c r 方程的l u s g s 算法 ( 六) 确定新顶点的形成点 确定新的v o r o n o i 顶点的形成点。新的顶点的形成点必然是包含新近引入点和 新近引入点的另外两个邻居，而且这两个点必须构成相邻三角形的一条边。在此例 中，新顶点v 。的三个形成点是新近引入点p 和p 的两个相邻点3 ，4 。其他顶点的形 成点在表3 中。 ( 七) 确定新顶点的相邻顶点 根据第六步，所有新顶点的形成点都已被找到，对每一个新顶点，搜索在步骤 五中所确定的每个相邻点的三个形成点以及其他新顶点的三个形成点，如果其中的 某一个顶点的三个形成点中有两个与该新顶点中的两个形成点重合，则这一个顶点 是该新顶点的相邻顶点，在本例中，新顶点v 。的三个相邻顶点是旧顶点v 。，以及新顶 点v v 7 。 ( 八) 重新整理v o r o n o i 图的数据结构 将v o r o n o i 图的数据结构进行调整，改写被删除的顶点。在本例中v 。、v ：是被 删除的顶点，为了有效利用空间，我们将这些被删除的顶点改写成新的顶点，并给 出相应的数据结构。在本例中分别用v ，、v ：来取代砜、v ，这样就完成了按理工作， 形成了新的三角化结构。整理后的数据结构由表3 给出。 ( 九) 引入新的点 重复第( 二) 步第( 八) 步的过程，不断引入新的点，直至 p 。 中的点全部参 加到平面的划分中去，这样也就完成了d e l a u n a y 三角化的进程。 、 以上的步骤相对比较繁琐，工作量会随着f p 。) 中点的个数的增加而呈平方级增 大。因此，我们可以采用一种树状搜索的方法，将有效地加速其进程。只要能给出 生成足够数量节点的方法，并使生成过程系统化，则d e l a u n a y 二角化方法可具有相 当的自动化程度。 顶点形成点相邻顶点 v il2pv 2v 6v 7 v 2l p4v -v v t v 3141v 2一ll v 44 3一l v 。 一11 v 532一lv ，一1一l v 6211v 一l一1 v ，3p2 v rv iv 5 v 8 4p3 v 2v ，v 1 表3 2 2 1 3 退化的处理 通常情况下，任意一组点在平面上的d e l a u n a y 三角化是唯一的，但在某些情况 1 2 南京航空航天大学硕 学位论文 下，二角化却并不唯一，这种情况被称作退化。这一点从v o r o n o i 多边形的定义即 可指导。在二位问题中，可能发生的退化有以下几种： ( 1 ) 两个点重合 ( 2 ) 应该构成三角形的三个顶点共线 ( 3 ) 四个或四个以上的点共圆 对于这二种情况，只要在计算中稍加处理即可避免。 2 2 1 4 网格的优化 按以上方法所生成的三角形，其形状尽管相对于给定的节点是优化的，但它必然要 受到节点位置的影响，所以不可能所有网格的性状都尽如人意。如果适当地调整网格点 的位置，网格质量将进一步提高，这称为网格的光顺。一种常用的网格光顺方法是 l a p l a c i a n 光顺方法。这种光顺技术是通过将节点向这个节点周围的三角形所构成多边 形的形心移动来实现。如果p ( x ，y ，) 为一个内部节点，n ( p ) 为与p 。相连的节点总数，则 光顺迭代技术可表示如下： 世 x ，= 2 j x t u ( t ) 怠) ( 1 6 ) ( p ) 、 y 。= y t n 【p 1 ) 光顺的方法可以使网格更均匀，三角形更接近于正三角形。均匀的网格对于保证计算的 精度和收敛性都有着极其重要的作用。图5 给出了优化前后网格局部图。 v 、， j 一 q 一入万 卜朋 优化前 图5 优化前后的网格 优化后 三堡韭堕塑旦堡! ! ! ! ! 查堡塑! 旦：! 璺! 茎鲨 一 2 2 1 5d e l a u n a y 方法的发展 构造d e l a u n a y v o r o n o i 图形的方法有b o w y e r 算法 1 i 和w a t s o n 算法 1 2 等。在网 格生成的应用中较广泛的是采用b o w y e r 算法。d e l a u n a y v o r o n o i 图形只解决了在给定 点集情况下如何连成一组三角形的问题，从网格生成角度看还需首先产牛点集。早期应 用d e l a u n a y 方法的工作使用结构网格技术产生点集，然后用d e l a u n a y 方法连接以形成 整个全机外形的非结构网格。作为网格生成的应用，d e l a u n a y 方法后来发展的重点在于 研究如何自动生成任意外形的点集。一种常用的方法是将计算域分成一些较大的了域， 在对各子域分别按修正的四叉树( q u a d t r e e ) a 叉树( o c t r e e ) 分解技术剖分，直到最 大单元的尺度与边界上离散的尺度相当为止。近来的发展为推进阵面的思想产生点集， 再用d e l a u n a y 方法连接，或用迭代点插入和局部重连等方法连接。重连是根据满足所要 求的质量准则优化的反复连接过程，因而大大改进了所生成网格的质量。 2 2 2 推进阵面法 2 2 2 1 布置控制点和生成背景网格 推进阵面法是网格和节点同时产生的非结构网格生成方法 1 3 。根据网格密度的需 要，在平面上布置一些控制点并给每一控制点定义一个尺度，连接这些控制点即可将平 面划分成大块的三角形背景网格。每一三角形背景网格中任一点的尺度都可以根据该二 角形三个顶点的尺度插值得到，因此相当于布置了一个遍及整个平面的网格尺度函数。 推进阵面法生成非结构网格的质量对背景网格的依赖性很大，如何生成较好的背景网格 对此法是否成功极为重要。 背景网格一般有结构和非结构两种。传统的非结构背景网格一般由人工构造，对复 杂三维问题这是相当复杂且费时的工作，故不易适应复杂边界的情况；而且在非结构背 景网格上插值需要复杂的查询过程，长时间的查询与插值也必然影响推进效率。结构背 景网格通常通过求解带源项的热传导方程来光滑分布控制参数 1 4 ，自动化程度较高， 但需求解非线性偏微分方程。文献 1 5 提出的采用d e l a u n a y 方法生成背景网格，可实现 背景网格生成的自动化。其主要步骤如下： ( 1 ) 沿流场边界布置若干节点控制点。非结构网格的优点之一就是可以根据要求 在任意的位置上进行加密。这种加密是通过在相应的位置上布置控制点和给控制点一个 大小合适的尺度来实现的。网格要求加密的地方尺度要小些，反之尺度要大些。控制点 上尺度的大小将与该点周围所生成的三角形的边长大致相当。 ( 2 ) 用d e l a u n a y 准则将控制点连成三角形。控制点布置好以后，可以根据d e l a u n a y 南京航空航天大学硕i ：学位论文 三角方法对平面进行三角化，所生成的三角形网格即背景网格。 ( 3 ) 在三角形顶点上随意定义网格尺度，但为实现物体边界与远场边界网格尺度的非 线性变化，在流场中嵌入一个制定网格尺度的椭圆人为边界以更好的控制物面边界附近 网格的质量。 ( 4 ) 生成背景网格的目的在于将离散点上的尺度插值到整个平面上。每一个背景网格 内的每一个网格内的尺度可以由它的三个顶点插值得到。因此，在以后的网格生成过程 中，当一个新的面积单元和一个新的节点要产生时，只需根据基边中点处的尺度即可确 定面积单元的大小和节点的位置。 2 2 2 2 定义边界和阵面初始化 在对控制体进行三角化以前，先要定义边界( 如图6 ) 。控制体的外边界上的边要按 逆时针方向组织，相邻的两点连接，形成一个闭合的有向环。内边界也是一个一个的有 相闭环，有多少个物体，就有多少个闭环。每一个内边界的边都要按顺时针方向排列。 这些边将被作为保留边，在最终的三角化结构中，这些边都将成为三角形的边。 图6 三角化的最初阵面一般就是所有的内边界的总和。阵面由两个整型向量组成。有多 少条活动边，阵面向量就有多少个元素。在一个向量里，储存活动边的一个点，在另一 个向量里储存活动边的另一个点。因为内边界是按照顺时针方向组织的，所以将要被二 角化的控制体始终在活动阵面的左侧。边界始终是不变的，而阵亟却要随着三角化的推 进而不断更新。随着三角形的生成，新生成的三角形的两条新边将成为活动边，丽二角 形的老边将从活动边向量中去掉。直到遇到外边界，阵面上活动边的个数减少到零，二 = 维非结构网格e u l e r 方程的l u - s g s 算法 角化过程便结束。 2 2 2 3 推进生成三角形 初始化阵面形成以后，就可以根据有关规则推进生成格点和三角形了 1 6 。在牛成 三角形的进程中，遵循以下步骤： ( 1 ) 在阵面上选择作为基边的活动边。一般来说，基边就取最小边( 经大小排序后 的最后一边就是最小边) 。在图7 中，我们选取a b 边作为基边。 ( 2 ) 假设a 和b 是选取的基边的两个端点，而边a b 的中心为m ( 如图7 所示) 。按 照一定的方式对背景网格进行搜索，直到找出m 点落在了哪个背景网格里面。由前 面已经知道，背景网