（飞行器设计专业论文）非定常euler方程解及在颤振中的应用.pdf

西北工业大学硕士学位论文 摘要 本文在深入研究非结构运动网格技术的基础上，结合基于a l e 描述的e u l e r 方程对二维和三维非定常流场进行了数值模拟；并以e u l e r 方程为主控方程，耦 合结构运动方程对二维翼型进行了颤振特性研究。 在经典线性弹簧方法的基础上，通过引入诸如抗拉、抗扭、抗压和边界修正 等方法对其进行了改进，由计算结果可知改进后的方法明显提高了非结构运动网 格的稳定性和适用性。 在流场求解技术方面，以非结构运动网格为基础，在时间推进上采用全隐式 的双时间推迸方法；在空间上采用格心格式的有限体积法对非定常e u l e r 方程进 行空间离散，用四步r u n g e - - k u t t a 方法作显式的拟时间推进，同时采用当地时 间步长、隐式残值光顺等加速收敛措施；求解了二维n a c a6 4 a 0 1 0 翼型的跨、 亚音速的非定常流场以及三维m 6 机翼的跨音速非定常流场。计算结果与文献及 实验结果符会的比较好，说明了本文所采取的流场求解技术是可靠的。 在气动弹性方面，在求解e u l e r 方程的基础上，耦合颤振运动方程，发展了 一种求解二维颤振问题的解算器。并以二维n a c a6 4 a 0 1 0 翼型为算例，求解了 在跨、亚音速下的颤振边晃，得到了其结构响应的时间历程。由计算结果和实验 值的对比可以看出本文所发展的解算器在进行二维颤振特性分析方面是可行的。 关键字：非结构动网格，e u l e r 方程，非定常，颤振 西北工业大学硕士学位论文 a b s t r a c t i nt h i s d i s s e r t a t i o n ，t h ee u l e re q u a t i o n s b a s e do nt h e a r b i t r a r y l a g r a n g i a n - e u l e r i a n ( a l e ) f o r m u l a t i o na r ea p p l i e dt os o l v e2 - da n d3 - du n s t e a d y f l o wf i e l dw i t hf u r t h e ri n v e s t i g a t i n gt h eu n s t r u c t u r e dd y n a m i cg r i dt e c h n i q u e t h e nt h e a e r o e l a s t i c i t yi si n v e s t i g a t e db yu s i n gt h ee u l e re q u a t i o n sa n ds t r u c t u r ee q u a t i o n s o nt h eb a s i so fc l a s s i cl i n e a rs p r i n ga n a l o g y , s e v e r a lm o d i f i e dm e t h o d sf o rt h e l i n e a rs p r i n gs t i f f n e s sc o e f f i c i e n ta r ei n t r o d u c e d t nt h e s em e t h o d st h es e g m e n t s ，f a c e s a n db o d i e so f t h ec e l la r ec o n s i d e r e d i nt h ec a s eo f n o tm u c hp e n a l t yi nc o m p u t i n g ，t h e a d a p t a t i o no f t h eu n s t r u c t u r e dd y n a m i cm e s h i sg r e a t e re n h a n c e d i nt h ep a p e r , 2 - da n d3 - du n s t e a d yf l o wf i e l d sa r ec o m p u t e d t h er e s u l t so f s o l v i n gt h eu n s t e a d yf l o wf i e l d st h a ta l ei na g r e e m e n tw i t ht h ee x p e r i m e n t ss h o wt h a t t h eu n s t e a d yf l o wf i e l ds o l v e ri sc o r r e c ta n dt h eu n s t r u c t u r e dd y n a m i cm e s ht e c h n i q u e o ft h ed i s s e r t a t i o ni sw e l lr o b u s t s ot h eu n s t e a d yf l o wf i e l ds o l v e rc a nb ea p p l i e dt o v a r i o u sk i n d so f p r o b l e m s 、丽t l lm o v i n ga n dd e f o r m i n gb o d i e s f o ra e r o e l a s t i c i t y , t h ef l u t t e rc h a r a c t e r so f2 - da i r f o i la f ef u r t h e ra n a l y z e db y c o u p l i n gt h ee u l e re q u a t i o n sa n ds t r u c t u r ee q u a t i o n s n l er e s u l t so ft h i sp a p e ra l e b a s i c a l l yi d e n t i c a lw i t l lt h ee x p e r i m e n t s s ot h ec a l c u l a t o ri sg o o df o ra n a l y z i n gf l u t t e r c h a r a c t e r s k e yw o r d s ：d y n a m i cu n s t r u c t u r e d 鲥d ，e u l e re q u a t i o n s , u n s t e a d y , f l u t t e r 1 1 1 西北工业大学 学位论文知识产权声明书 本人完全了解学校有关保护知识产权的规定即：研究生在校攻读学位期间论文上作的 知识产权单位属于西北工业大学。学校有权保留并向国家有关部或机构送交论文的复印僻 和电子版。本人允许论文被查阅和借阅。学校町以将本学位论文的全部或部分内容编入有关 数据库进行检索，可以采刚影印、缩印或扫描等复制手段保存和汇编本学位论文。同时本人 保证，毕业后结合学位论文研究课题再撰写的文章一律注明作者单位为两- i l ： ：业大学。 保密论文待解密后适用本声明。 ，学位论文作者签名：i 鲢登指导教师签名：乇7 瓮了芭、 五司年3 月上。日z - 。0 1 年0 月z 。日 西北工业大学 学位论文原创性声明 秉承学校严谨的学风和优良的科学道德，本人郑重声明：所呈交的学位论文，是本 人在导师的指导下进行研究1 作所取得的成果。尽我所知，除文中已经注明引用的内容 和致谢的地方外，本论文不包含任何其他个人或集体已经公开发表或撰写过的研究成果， 不包含本人或其他己申请学位或其他用途使用过的成果。对本文的研究做出重要贡献的 个人手u 集体，均己在文中以明确方式表明。 本人学位论文与资料若有不实，愿意承担一切相关的法律责任。 学位论文作者签名赵逸毽 j 7 年3 月2 - o h 西北工业大学硕士学位论文第一章绪论 1 1 引言 第一章绪论 流体力学是研究流体与物体之间有相对运动时流体运动的基本规律以及流体 与物体之间的作用力的一门自然科学。而计算流体力学( c f d ) 是现代流体力学 中的一个重要学科分支，它是计算数学、计算机科学和流体力学相结合的产物。 随着计算技术的迅猛发展，即高速计算机的出现和数值计算方法的不断发展，计 算流体力学技术上已经日趋成熟。现在已经能够通过求解e u l e r 或n s 方程来进 行流场的数值模拟，。这种数值模拟相比较风洞实验而言，有设计周期短、成本低、 通用性好等优点。由于这些优点的存在，使得计算流体力学已经作为一种基本的 研究方法和设计手段逐步进入并深入了航空、航天、船舶等领域【1 】【2 】，已经成为 人类航空航天科学技术的一个不可缺少的重要组成部分。 计算流体力学经历了数值求解拉普拉斯方程、小扰动速势方程、全速势方程、 欧拉方程和n a v i e r - s t o k e s 方程( 以下简称n s 方程) 等发展阶段。2 0 世纪8 0 年代以前，由于计算机水平的限制，计算流体力学的数值模拟主要以求解拉普拉 斯方程、小扰动速势方程、全速势方程为主，其中有代表性的是基于拉普拉斯方 程的面源法以及有限差分法求解小扰动速势方程和全速势方程。在随后的二十多 年中，由于计算机技术和数值计算方法的迅猛发展，计算流体力学在求解欧拉方 程和n s 方程以及数值模拟复杂流场方面取得了重大突破。在此期间，计算流 体力学数值模拟有了很多新方法：有限差分法、有限体积法、有限元法等。随着 计算对象复杂性与计算程度要求的提高，计算方法与计算格式有了许多创新，诸 如t v d 格式、e n o 格式、n n d 格式等高精度、高分辨率差分格式的提出。计 算流体力学对激波、漩涡等复杂问题的模拟能力也有了很大的提高。与此同时， 计算流体力学的网格生成技术也获得了飞速的发展。在结构网格方面出现了代数 生成网格法、解微分方程生成网格法、保角变换法等多种网格生成方法。网格类 型也由单一的c 型网格、o 型网格、h 型网格发展到嵌套网格和多块对接网格 等。非结构网格方面出现了阵面推进法、d e l a u n a y 方法、八叉树方法等多种非 结构网格生成方法。近年来，各种类型的结构、非结构杂交网格技术更是层出不 西北工业大学硕士学位论文第一章绪论 穷，使计算流体力学解决复杂流动和精细流动问题的能力有了很大的提高。 以计算流体力学( c f d ) 和计算结构力学( c s d ) 为基础，数值模拟流固耦 合问题并且进行流固祸合问题的研究是近年来发展的主要方向。根据学科间求解 问题的密切程度，提出了学科间紧耦合( t i g h tc o u p l i n g ) 和松耦合( l o o s ec o u p l i n g ) 的概念。所谓紧耦合，是指将流体和固体控制方程合写成统一的有限元形式，将 c f d 和c s d 程序合并成一个整体，在每个时间步同时求出所有变量；而松耦合 则是每个时间步先后分别计算c f d 和c s d 程序。松耦合方法易于实现并且能够 利用现有技术，应用较广，尤其是模拟已知物体运动状态的非定常气动力的情况。 l d h n e r 等人的其它研究工作还有二维和三维外挂物分离【3 l 【4 1 、三维爆炸波与掩体 的相互作用【5 1 、飞行员、座椅与载机的分离【6 1 、全机外形副油箱投放过程的模拟 7 1 等，到1 9 9 9 年已经发展到可以模拟数以百计的弹片飞散过程i 研的水平。h o s a n g a d i 、 l e e 、c a v a l l o 、d a s h 等人则在考虑了粘性、湍流、化学非平衡的计算方法中还包 含了网格的变形运动，成功求解了子弹药散开绕流i 9 , 埘，到2 0 0 0 年，增加局部重 构网格方法模拟弹头防护罩抛开的研究工作 1 1 , 1 2 也已完成。 1 2 非结构运动网格的发展现状 计算网格的生成技术一直是研究计算流体力学的重要前提条件，网格生成的 质量直接影响数值计算的结果。 由于非结构网格有利于处理复杂几何外形以及拥有随机的数据存储结构，使 得近年来非结构网格技术得到了顺速的发展和应用。目前，常用的非结构网格生 成方法有d e l a u n a y 方法【1 3 】【1 4 1 和阵面推进方法( a d v a n c i n g f r o n t m e t h o d ) 1 5 1 7 j 。 d e l a u n e y 方法的基本思想是：首先在计算空间按某种方法分靠一系列的网格 节点，然后按照一定的准则将这些点连接成网格单元。这些准则中最常用的是 d e l a u n c y 球面准则，对二维情况该准则要求每个网格单元的外接圆不包含其它节 点，对于三维情况要求四面体单元的外接球不包含其他网格节点。d e l a u n a y 方法 的主要优点是网格生成速度快、网格质量较高，其缺点是需要事先提供初始网格， 网格的疏密变化不易控制，对比较复杂的外形，网格生成时的壁面穿透问题也很 难避免。阵面推进方法( 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 方法而言，其网格生成的效率较低，网格质量稍差，但通过改进数 据结构和算法以及对生成的网格进行光顺和优化等后置处理，上述阵面推进方法 的缺点是可以克服的。 在非结构网格范畴内，有多种处理边界运动的方法。其中一种方法是网格固 定，物体边界运动时与网格单元形成切割。这种方法在每个时间步都必须处理边 界上的切割单元，并且要不停地加密和粗化网格以满足分布要求 8 1 。此方法多用 于非结构笛卡尔网格，还有一种方法是网格与运动边界固连，随边界一起作刚性 运动。刚性运动方法用于边界无变形和粘性绕流模拟中的高度的拉伸的边界层网 格。 在处理边界变形和较大的边界运动时，非结构动网格一般采用网格变形 ( d e f o r m i n g ) 与局部重构( l o c a lr e m e s h i n g ) 相结合的方法。如果动边界位移较 小，如静气动弹性问题，仅仅依靠网格的变形就能适应边界的运动。网格变形区 可以局限在运动边界周围一定范围内，以此提高计算效率，也可以是去除外边界 和固定边界外的整个流场。当处理位移较大的边界运动，如机弹分离问题时，可 将动边界周围一定范围设成网格变形区，该区网格随着边界的运动而变形，当出 现严重扭曲的网格单元时，重新生成变形区内的网格。流动参数通过插值运算从 旧网格映射到生成的新网格。 网格变形能力及变形后网格质量的好坏取决于控制网格变形的物理模型的构 建。常用的变形模型有弹簧近似( s p r i n ga n a l o g y ) 模型【1 8 , - 2 m 、弹性体( e l a s t i c i t y m e t h o d ) 模型 2 t 】和代数模型1 2 2 1 。弹簧近似方法最早被b a t i n a 用来求解翼型强迫振 动绕流。其基本思想是将网格单元的各条边看作弹簧，弹簧系数与网格边的长度 有关。当边界运动后，通过求解弹簧系统节点受力平衡问题确定网格点的新位置。 弹性体方法是将计算区域比作一个线性弹性体，通过求解弹性力学方程组来确定 网格点的位置。所谓代数模型，是指网格点的位移由动边界位移乘以一个系数得 到，该系数在动边界上取l ，在流场外边界取o ，其问按一定函数规律插值。这种 方法效率最高，但仅适于变形区形状非常简单的情况。弹簧近似方法由于没有考 西北工业大学硕士学位论文第一章绪论 虑剪切应力的作用，变形能力不如弹性体方法，但计算效率高，已成功用于气动 弹性分析、外挂物分离、网格自适应调整、气动设计等。弹性体方法考虑了所有 应力，变形能力强，但计算工作量大，效率低。b l o m 对标准弹簧近似方法进行了 改进，加入了扭转弹簧效应，使其变形能力大为提高。k e n n o n 等还提出了一种不 可压流模型实现网格点的运动。其基本思想是将网格点或单元看作不可压流体微 团，当物体边界运动时，使网格点绕运动边界“流动”，然后用d e l a u n a y 方法重 新连接网格点成为新的网格。这种方法实质上也是一种重构网格。由于拓扑结构 每个时间步都有变化，频繁的插值运算是不可避免的。 局部重构网格方法始见于f o r m a g g i a 2 3 】对二维外挂物投放的数值模拟中。此 外，l i i h n e r 及其研究小组发展了这一方法i 鲫】。他们的共同之处在于变形区即重 构区内的网格重新生成时采用标准的网格生成方法，如阵面推进法等。b a u m 指 出【5 j ，重构网格后的插值运算能导致数值耗散，重构越频繁，耗散损失就越大。 增大变形区的范围，使更多的网格单元分担边界的位移，可以降低局部重构操作 的频度，但同时也增加了每一步网格变形的计算工作量，这是一个优化问题【7 l 。 为了降低耗散损失，b a u m 及l u h n e r 等近年的研究工作中变形区范围逐渐增大， 1 9 9 4 年的文献1 5 】中由原来得两层单元扩大到动边界周围得十层单元，到1 9 9 7 年 的文献【7 l 再次增大到二十层单元，并发展了一种基于拉普拉斯算子的刚性模型控 制网格变形。但总的来说，这些研究工作仍然可以在网格变形模型方面作进一步 的改进。c a v a l t 曾提出一种分区概念【lo 】，变形区设成四边形并随运动物体刚性平 动，在分区界面上用层推进方法实现网格添加，运动物体姿态变化由变形区内网 格变形来适应。这种方案没有被真正应用。其后来的文献 n , 1 2 1 中应用了局部重构 方法，在变形区采用了弹性体模型，取得了很好的结果，但计算效率较低。 1 3 流场数值模拟的发展现状 流场数值模拟的实质就是通过差分或其他办法把流动的控制方程由非线性的 偏微分方程组转化为只包含简单加减乘除等运算的代数方程组，并且将这些代数 方程组和具体的流动边界条件结合起来，编制成程序后让计算机进行求解，从而 获得所需流场参数的过程。 最近二十年里，以e u l e r 方程和n s 方程为控制方程的流场数值模拟技术获 4 西北工业大学硕士学位论文第一章绪论 得了巨大的成功，取得了丰硕的应用成果。基于欧拉方程的流场数值模拟技术不 但能够模拟和捕捉激波，而且能够有效地模拟有旋流和自动捕获由压力差引起的 分离涡，自动形成分离线和尾迹涡面，揭示许多流动现象的本质特征。相对e u l e r 方程而言，基于n s 方程的数值模拟技术发展相对缓慢一些。由于计算机容量 和速度的限制，完全n s 方程的求解在短时期内还难以进入实用阶段。目前，n s 方程的求解还是建立在雷诺平均n s 方程的基础上，借助适当的湍流模型 来实现的。而湍流问题至今仍然是困扰整个流体力学界的一个难题，湍流的基本 机理至今还没有完全弄清，这就决定了目前的各种湍流模型都有各自的局限。目 前尚不存在一个通用的湍流模型，既能在相当广泛的流动情况中反映出比较正确 的物理特征，从而得到足够准确的计算结果，又能够为当前的计算资源所接受， 可以这么说，湍流模型问题在一定程度上制约了n s 方程数值模拟技术的发展。 当前数值模拟e u l e r 方程和n s 方程来求解非定常运动问题的时间推进方法 主要有三种；显式、隐式和双时间( d u a l - - t i m e ) 。由于显式方法受稳定性条件的 限制，时间步长只能取得非常小，并且不能使用定常计算的一些加速收敛措施， 如：当地时间步长，焓阻尼和多块网格技术等。因此采用显式方法求解非定常问 题需要太长的计算时间；隐式方法的时间步长不再受稳定性条件限制，时间步长 麓 可以取得较大，但求解庞大的矩阵方程对计算机的速度和内存等方面提出了严峻 的考验，特别是求解三维问题；j a m e s o n 提出的用于非定常计算的双时间方法在 计算效率方面综合了上述显式和隐式方法的优势。该方法是一种全隐式方法，物 理时间步长可以取得足够大，真实物理时间的收敛精度可达二阶以上，拟时间的 推进又可以采用显式的加速措施，此外编程方便，能进行向量运算和分区的并行 运算，也是该方法的另一个优势。 当前数值模拟e u l e r 方程和n - - s 方程的空间离散方法也主要有三种：有限差 分法、有限体积法 2 4 1 和有限元方法，在非结构网格中常用的是后面两种。有限体 积法的主要思想是把流体力学的积分型控制方程应用到每一个由网格剖分得到的 控制体微元上，从而把积分方程的求解转化为代数方程或常微分方程的求解，再 利用各种时间推进格式进行迭代，最后在每个控制体微元上得到流场的稳定解。 有限元方法的数学基础是变分原理和加权余量法，其主要思想是在每个网格单元 上选择若干合适的点作为解函数的插值节点，方程的近似解由各网格单元的近似 函数逼近，而单元中的近似函数又由单元基函数的线性组合得到。这两种方法各 5 西北工业大学硕士学位论文第一章绪论 有所长，有限元方法的数学概念较浓厚，对网格形状的适应性较好，有限体积法 的物理概念十分明确，方法比较简便，容易被接受和掌握。许多研究者已经证明 了分片线性插值近似的情况下的有限元方法和有限体积法是等价的口5 堋。本文采 用的是有限体积法。 1 4 飞行器气动弹性研究的现状和前景 一直以来，飞行器的气动弹性问题都是飞行器设计必须考虑的重要课题。从 固定机翼飞行器的问世的第一天起，人类就开始面临飞行器的气动弹性问题。1 9 0 3 年l a n g l e y 的单翼机在作有动力的首次飞行时，就因为典型的静气动弹性扭转发 散而折断了机翼。近1 0 0 年来，随着飞行器研制的不断进步，飞机的气动弹性问 题也日益成为影响飞行器飞行性能的难题，美国的f 1 8 一e f 战斗机的研制过程 中，机翼在一定迎角下飞行时就遇到静气动弹性的扭转发散的问题。我国在国产 飞机的自行研制和型号改进的时候，也不只一次的遇到过抖振、嗡鸣等气动弹性 问题。由于飞行器的气动弹性问题常常严重影响飞行器的飞行性能，甚至会造成 飞机结构的破坏，造成空难事故的发生，所以必须深入研究和解决。 由于飞行器气动弹性力学属于多学科交叉的学科，其研究的问题同时涉及到 流体力学、固体力学以及控制理论等多个学科的难点，因此飞行器的气动弹性一 直以来都是飞行器设计者必须面对的重要研究难题。在传统的数值计算方法( 包 括从国外进口的商业软件包) 中，使用的气动力模型都是以势流、线化方程为理 论基础的非定常升力面方法，虽然此技术途径的计算效率较高，但其有一定的缺 点： 1 非定常升力面理论是以零迎角为基础的数值方法，不能反映出迎角的影 响，只能适用于传统的附着流型，但许多飞行器已经广泛采用分离流型 ( 即飞行器的机翼在中等迎角、甚至5 、6 度迎角时，机翼上已出现前缘 分离涡) ，此时机翼上的气动力分布时有差别的，因此，线性非定常升力 面方法不再适用。 2 由于使用的基本方程是势流理论的线化方程，所以在亚音速、超音速阶 段较适用，在高亚音速和跨音速阶段不适用，近些年国内学者在跨音速 小扰动方程的前提下，对非定常升力面理论进行了改进，扩大了气动弹 6 西北工业大学硕士学位论文第一章绪论 性的计算范围，但严格来讲，其理论上的适用范围还很有限，它在激波 的模拟方面存在不足，另外，它无法模拟激波随结构变形所产生的运动 变化。 3 由于只是在频域中计算，所以只能计算颤振速度、颤振频率两个基本参 数，不能计算亚临界、超临界时结构弹性振动的过程，而人们有时很希 望知道亚临界、超临界时结构弹性振动的过程。 4 传统方法不能分辨出超l 临界机翼与普通机翼的非定常气动力之间的差 别。 5 不能考虑实际飞行中有迎角时的静载对弹性结构引起的静变形，而静变 形实际上又改变了机翼各部分的迎角分布，甚至静载引起的静变形会超 出弹性结构的线性范围。此外，传统的颤振计算方法难以加入结构非线 性的影响。 由以上可以看出，常规的、以线性小扰动为基本假设的气动弹性颤振分析理 论，由于无法正确地揭示非线性气动弹性现象发生的机理，因此，必须建立和发 展相应的处理非线性气动弹性问题的分析理论与技术。二十世纪九十年代中后期 以来，由于c f d 和c s d 的发展，特别使c f d 中基于e u l e r 方程和n s 方程的数 值模拟。c s d 中有限元计算方法的日益成熟，为飞行器气动弹性问题，特别是关 键的跨音速气动弹性问题的研究提供了条件 2 7 - 3 0 l 。由于气动弹性的数值计算是建 立在e u l e r 和n s 方程的基础上进行定常和非定常求解的，所以计算量非常的大。 直到二十世纪九十年代后期，由于计算机硬件的发展，多c p u 计算i 作站和大规 模并行机群的研制和应用才得以解决这问题。目前，国内外在静气动弹性、颤 振、抖振、嗡鸣、伺服气动弹性等气动弹性问题的数值计算已经有大量的文献 3 3 - 3 6 1 1 5 9 - 6 1 】发表。 1 5 本文工作重点 根据目前国际上气动弹性研究的发展趋势以及国内航空航天技术发展的需 要，本文以计算流体力学方法为基础，主要工作有以下几个方面： 1 ) 在非结构网格的基础上，深入的研究了非结构运动网格技术。通过对标 准线性弹簧及其改进方法在非结构运动网格中的应用对比，可以看出改 7 西北工业大学硕士学位论文第一章绪论 进后的方法明显提高了非结构运动网格的适应性和稳定性。 2 ) 在二维非结构运动网格技术的基础上，通过求解二维e u i e r 方程深入细 致的研究了谐和振动n a c a6 4 a 0 1 0 翼型的跨音速和亚音速非定常气动 力。 3 ) 在三维非结构运动网格生成技术的基础上，通过求解三维e u l e r 方程研 究了m 6 机翼柔性谐和振动的跨音速非定常气动力 4 ) 在求解a l e 形式的e u l e r 方程的基础上，耦合颤振运动方程，发展了求 解二维颤振问题的计算方法和程序，得到了结构响应的时间历程。发展 了一种颤振系统响应特性的分析方法。 3 西北工业大学硕士学位论文第二章非结构运动网格技术 第二章非结构运动网格技术 在非结构运动网格的生成方法中，无论是将网格固定，边界运动切割网格的 切割网格法还是将网格和边界固连，并一起运动的刚体运动法，都有其方法的缺 点和局限性。所以我们在处理这类问题时基本上都是采用网格变形和局部重构相 结合的方法。常用的网格变形方法有弹簧近似方法和弹性体方法。尽管弹性体方 法具有较强的变形能力，但由于其工作量大，而且本文处理的是小位移问题，所 以本文采用的是弹簧近似方法，这样在很大的程度上可以提高计算效率。由于标 准的弹簧近似方法没有考虑到剪切应力的作用，变形能力不强，故本文采用了改 进的标准弹簧近似方法，计算表明改进后的方法大幅度提高了网格的变形能力。 2 1 顶点弹簧方法 弹簧近似方法认为构成网格的每一条边都是具有一定倔强系数的弹簧，所以 在计算弹簧张力时就要定义弹簧的平衡长度。根据b l o m 的定义，弹簧近似方法 中计算平衡长度有两种方法，种是顶点弹簧( v e r t e xs p r i n g s ) 法，规定弹簧的 平衡长度为零；另一种是棱边弹簧( s e g m e n ts p r i n g s ) 法，规定弹簧的平衡长度 等于边界未运动和变形前边的原长。本文采用的是顶点弹簧法。 传统的顶点弹簧近似方法一般用于网格优化，倔强系数取为常数，网格点移 动后的新位置坐标相当于周围点的坐标的平均。本文由于要实现网格的变形运动， 故对这一传统方法进行了如下的改进。 对于顶点弹簧，节点f ，_ ，间的弹簧张力为： 元= k o e 一再) ( 2 1 ) 式中k u 为连接节点f ，的弹簧的倔强系数；只，弓分别为节点j ，的位置矢 量。则对任意的一个节点i 列它的合力方程为： 丘：兰巧e i ) ( 2 2 ) j - i 式中m 为计算区域内所有与节点i 相连的非结构网格的节点总数。不同于网格优 9 西北工业大学硕士学位论文 第二章非结构运动网格技术 化的情况，修改后的顶点弹簧方法令计算区域内所有网格节点的受力都为初始状 态即边界没有发生运动和变形时所受的合力。对计算区域内所有的m 个非结构网 格的节点应用( 2 2 ) 式可得到初始状态下的该线性系统的表达式： 一 k 尹= f ( 2 3 ) 式中： k = n 。 一k i ，k l ， k l j = i ； k 。 ； k _ 1 m 一蜀 k j - t m 一艺 j - l r = 一 ，f = 曩 ： 二 e ： 二 f 卅 式中为包括边界点的网格节点数。很显然( 2 3 ) 式是一个i t t 阶的线性方程组， 应用g a u s s - s e i d e l 迭代法求解得到其相应的分量迭代式为： 尹虬球夸川一黔 眨4 ， ( f = 1 ，2 ，m ；七= 0 ，1 ，2 ，) 由于该方程组的刚度矩阵k 正定对称，故( 2 4 ) 式关于任意初始向量均收敛。当 边界变形或运动到栉+ 1 时刻，利用( 2 4 ) 式经过数次迭代便可以达到需要的精度， 得到”+ 1 时刻的网格点的位置。 在与用中心有限体积法耦合求解流场时，需要用到网格单元的运动速度 t = 置7 + 只7 + 刁f ，可以由下列公式求出： v g = 3 ”一4 蛩町+ 曩” 2 f 一1 亡 2 i 己 | ，= l + o b 2 ) ( 2 5 ) 式中乏为网格形心点的位置矢量，霉为网格四个顶点的位置矢量，a t 为时间步长， 1 0 西北工业大学硕士学位论文第二章非结构运动网格技术 上标o + 1 ) 、o ) 和o 1 ) 分别表示不同的时间层。 2 2 弹簧倔强系数的选取 本节主要阐述了标准线性弹簧方法倔强系数的确定以及对标准线性弹簧改进 方法的引入。数值试验表明。改进后的弹簧方法在计算量增加不多的情况下，明 显的提高了该方法的变形能力和稳定性。 2 2 1 标准线性弹簧方法倔强系数的确定 标准线性弹簧方法大都选取网格边的边长为参照变量来确定弹簧的倔强系 数。根据大量的数值经验，一般按公式( 2 6 ) 计算，可以很好的保证二维或三维 网格的疏密特征。 k ：万1 ( 2 6 ) 式中f 。为网格边的边长。于是当z ，哼o 时，鬈。m o o 。这样便很好的保证了 任意的网格节点i ，_ 不会碰撞。 2 2 2 改进的标准线性弹簧方法倔强系数的确定 ( 2 6 ) 式只是考虑了弹簧在其直线方向的伸缩稳定性，却没有考虑到网格节 点穿透所在网格单元的面( 二维为边) 的情况。所以，当边界运动或变形比较大 时就很容易产生体积( 二维为面积) 为负的“无效”网格单元。为了防止产生这 种“无效”网格单元，我们引入了改进的标准线性弹簧方法【5 8 1 。这种方法主要有 以下改进： 第一，在标准的线性弹簧倔强系数公式中引入单元几何形状的变量并作以下 改进： 圹壶 q j 上式中为拥有网格边ij 的所有单元中体积最小的单元。当寸0 时， 西北工业大学硕士学位论文 第二章非结构运动网格技术 五乙。j m 。这样就保证了与网格边fj 相对的所有单元中体积最小的单元不会 成为体积为负的“无效”网格单元。 第二，为了防止网格运动时对四面体单元中的三角形单元的挤压，引入如下的 防挤压系数： k 2 赤 ( 2 9 ) 上式中吼为与网格边i ，所对的三角形内角。 第三，为了防止网格运动时对四面体单元的挤压，采用了文献中的简化方法， 引入防挤压系数k ，删，具体表达式为： 足删2 而1 2 1 0 ) 上式中吼为与网格边f ，所对的四面体单元所夹内角。 经过上述的改进和修正，最后得到总的弹簧系数为： 巧= 后d 0 乇。+ k 。+ k 删，自) ( 2 1 1 ) 上式中硒= 1 d 为倔强系数修正因子，d 计算区域内点到物面的最小距离。 2 3 非结构动网格的实现步骤及算例分析 上面所述的方法需要每一时间步都更新弹簧倔强系数，这样当网格作多步连 续运动时会大大增加计算量，甚至有可能使( 2 3 ) 式所示的线性弹簧系统变成非 线性系统，从而导致运动网格系统发散。所以本文在实际计算时，均将弹簧倔强 系数固定在初始状态下，数值计算表明该方法是可行的、可靠的。 本文算例中生成运动网格的算法流程图见图2 1 。图2 2 和图2 3 给出了 n a c a 0 0 1 2 翼型在每个时间步都更新弹簧倔强系数的情况下做俯仰运动的网格 图，由图可以看出改进的方法明显提高了动网格的适应性。图2 4 的( a ) 和( b ) 分别给 出的是应用标准线性弹簧方法和改进方法使n a c a 0 0 0 6 翼型偏转4 0 度时的状况， 由图可以看出改进后的方法明显提高了非结构运动网格的稳定性。 西北工业大学硕士学位论文 第二章非结构运动网格技术 2 4 小结 本章主要介绍标准的线性弹簧方法及其改进方法，并通过数值计算将两种方 法进行了比较，得到以下结论： i ) 改进的线性弹簧方法大幅度的提高了非结构运动网格的适应性。 2 ) 改进的线性弹簧方法保持了原有的标准线性弹簧方法计算效率高的特 点，并在保证网格稳定性的同时大大提高了网格的变形能力。 西北工业大学硕士学位论文 第二章非结构运动网格技术 2 5 附图 图2 1 非结构运动网格生成流程图 西= l g i 业大学硕士学位论文第二章非结构运动网格技术 图2 2 标准弹簧方法进行多步反复网格变化的情况 图2 3 改进弹簧方法进行多步反复网格变化的情况 西北工业大学硕士学位论文 第二章非结构运动网格技术 图2 4 ( a )n a c a 0 0 0 6 翼型在标准弹簧方法偏转4 0 。时的整体网格图以及 后缘“无效”网格放大状态 图2 a ( b ) n a c a 0 0 0 6 翼型在改进弹簧方法偏转4 0 。时的整体网格图以及 后缘“无效”网格放大状态 1 6 西北工业大学硕士学位论文第三章非定常e u l e r 方程的求解 第三章非定常e u i e l 方程的求解 本文采用j a m e s o n 提出的有限体积法对非定常e u l e r 方程进行求解。空间离 散上采用的是格心格式的有限体积法，这种离散方法的优点是算法比较容易实现， 数据存储简单，流场求解的分辨率较高等，缺点是离散精度和稳定性容易受网格 质量影响，而且很难用多重网格进行加速收敛；在时间推进上采用了j a m e s o n 提 出的用于非定常计算的双时间推进法，该方法有以下优点： 1 ) 能保证在真实物理时间的推进过程中流场计算的收敛精度达到使用者所 期望的高度，这对于模拟比较长的时间历程问题而言，是保证计算结果不漂移的 有力措施。 2 ) 它是一种全隐格式方法，物理时间步长理论上可以按使用者的愿望取足够 大，从而保证了计算方法具有一定的高效性，并且方法的时间步长不是受格式限 制而是受物理因素的影响。 3 ) 该计算方法又具有显式方法的优点，可以采用显式方法的加速收敛措施例 如当地时间步长、隐式残值光顺等。 3 1 基于a l e 描述的e u l e r 控制方程 描述流体运动有两大基本方法，、即拉格朗日( l a g r a n g i a n ) 方法和欧拉 ( e u l e r i a n ) 方法。拉格朗日方法是指网格随流体一起运动，网格点的速度与当 地流体微团的速度相同；而欧拉方法是指网格的空间位置固定，即网格点的速度 为零，流体微团穿过网格单元所构成的控制体。任意拉格朗日欧拉( a l e ) 方法 是这两种基本方法的统一，它允许网格以任意速度运动。当网格速度为零时则为 欧拉方法，当网格速度等于流体运动速度时则为拉格朗日方法。 本文应用的控制方程是a l e 描述的守恒型三维可压缩非定常积分形式的 e u l e r 方程，其直角坐标系下的积分形式为： 立s t 巧豆d v + 豇户( 跏五嬲= 。 1 7 ( 3 1 ) 西北工业大学硕士学位论文第三章非定常e u i e r 方程的求解 式中守恒变量项亘为： q = p p u p w e 对流通量矢量项户 ) 在x ，y 和z 轴的分量为； 曩 ) = 只 ) = 丘 ) ： p u p u u + p 钟u p w u ( e + p ) u + x ，p p v p u v p v v + p p w v ( p + p ) v + y p p w p u w p v w p w w + p ( 口+ p ) w + z ，p 逆变速度矢量圣= 听+ 巧+ 厩沿x ，y 和z 轴的分量u ，哪定义为； 此外对于理想气体压强p 为： ( 3 2 ) ( 3 3 ) ( 3 4 ) ( 3 5 ) ( 3 6 ) p = ( ，- 1 ) ( e 一妄p ( u 2 + i t 2 + w 2 ) ) ( 3 7 ) 以上各式中q 为控制体，a q 表示控制体单元的边界，d v 是体积微元，元是控制 体边界外法向单位向量，d s 是面积微元；p ，甜，v ，w ，e 分别表示流体的密度、 工，闸2 轴方向的流体速度和单位体积流体的总内能；t ，只和4 分别为网格沿 t b t 所刁 一 一 一 “ v 甜 = = = u 矿矿 ，j、【 西北工业大学硕士学位论文第三章非定常e u l e r 方程的求解 工，y 和z 轴方向的网格速度；对于理想气体比热比y = 1 4 。 3 2 无量纲化参数的引入 数值计算中，方程组中的控制参数，p ，甜，v ，x ，y ) 多采用无量纲化形式，选取 翼型弦长c ，自由来流的密度成，自由来流音速a 。为无量纲化的参考量，则无量 纲化变量如下： 孑：三 罗：上 三= 三 cc c 万：上 几 一 “ = 一 氏 芦2 石p 一 1 ， 钆 茁= 仁一u o ) 4 lg _ - 一虬) 3 3 空间离散 f ：鱼f c w w = ：一 有限体积法对e u l c r 方程进行空间离散的基本思想是：将积分形式的e u l e r 方程直接应用到每一个划分好的网格单元上，然后通过插值和差分的办法计算通 量项，将积分方程转化为以网格单元中心物理量为未知量的代数离散方程组进行 求解。本节的重点是在非结构网格下如何利用中心差分形式的格心有限体积法对 e u l e r 方程进行空间离散并添加消除振荡和抑制波动的人工粘性。 3 3 1 格心有限体积法 考虑某一四面体网格单元i 并将其简记为q ，定义磊为q ，内各物理量磊的 平均值，并且假设在网格中心位置磊= 磊，则方程( 3 1 ) 在q ，内可以作如下的 近似： 1 9 西北工业大学硕士学位论文 第三章非定常e u l e r 方程的求解 等毋d v + 豇施m 豳= 。 ( 3 8 ) 其中a 噙表示q 。的边界。上式在单元q ，上写成离散形式为： 烈dq j v 。) + e = 。 ( 3 9 ) 这里v 。为q ，的体积，e 为单元f 的通量项，计算用的表达式为： ( 3 1 0 ) 式中，瓦，和嬲。分别表示网格单元f 的第七个表面的外法线单位矢量和面积， 晟，( 亘) 表示网格单元珀q 第七个表面上户( 勐平均值。在数值计算中我们用网格各 表面面心处的户( 亘) 值来近似替代该表面上户( 豆) 的平均值。 由表达式( 3 1 0 ) 可知，通量项己的计算过程中需要用到网格各个表面面心 上的流动参数。在中心差分形式的格心有限体积法中，表面面心上的流动参数是 利用共面的两个相邻网格单元格心处的流动参数进行插值得到的。 3 3 2 人工粘性 使用中心格式的有限体积法，在用差分项代替微分项时，舍弃了高阶导数 项，导致某些高频误差分量在求解流场的过程中不衰减。如果在通量项中适当 地加入一些高阶导数项，就会避免这种振荡现象，使数值解完全地收敛到定常 状态。这些人为加入的高阶项，称为人工耗散。 e u l e r 方程的空间有限体积离散方程( 3 9 ) 在引入人工粘性后变为： 鲁( 蚕。v 。) 托_ = 0 其中丘为网格单元i 上的人工粘性耗散通量，与有限体积法的通量计算形式一致， 西，也是网格各单元面上的人工粘性通量之和。为了处理流场光滑区域和压力梯度 较大的区域里不同的振荡情况，需要引入不同形式的人工稻性项。在流场光滑区 i 嘏 t | 、jql r 。糊 l l c 西北工业大学硕士学位论文 第三章非定常e u l e r 方程的求解 域主要是引入四阶人工粘性，在压力梯度较大的区域引入二阶人工粘性。这样人 工粘性通量可以表达为： f i , = 犀2 + 西j 4 ( 3 1 2 ) 式中叫”，叫4 分别是网格单元f 的四个面上二阶和四阶人工粘性通量之和，为： 厨2 = 西 ( 3 1 3 ) 厨4 = 韶 ( 3 1 4 ) 这里下标壤表示该面由网格单元f 和网格单元七所共用，耀和耀的具体表达式 如下： 稚= 吼鼯( 磊一磊) 。 ( 3 1 5 ) 孑= 一口m 占( v2 亘t v 2 磊) m ( 3 1 6 ) 式中，v 2 磊= ( 色一磊) ，筇和为自适应开关系数，口耻为尺度因子二阶 和四阶人工粘性的比例分配是由f 、和盯“共同决定的。在激波附近二阶人 工粘性应该起主导作用，因此这里引入激波感受因子的概念，以实现人工粘性对 压力梯度的感应。其定义为： 驴剜 m 这样自适应开关系数和彰即由激波感受因子所决定 簖= 后2 v 雎 搿= m a x ( o ，k 4 一露) 以上式子中七( 2 ) 和七( 4