（流体力学专业论文）时—空守恒元解元方法数值研究.pdf

两北工业大学硕士学位论文摘要 摘要 时一空守恒元解元方法( 简称c e s e ) 是求解守恒律方程的一种全新的数值方 法。本论文对时一空守恒元解元方法进行了一系列的研究，建立了一套基于非结 构网格的c e s e 方法数值模拟程序，求解了二维e u l e r 和n a v i e r - s t o k e s 方程，并 将其应用于求解亚跨超音速流场和低速不可压缩粘性流场的计算中。 本文首先对c e s e 数值方法进行了介绍，包括守恒元c e 和解元s e 的构造， 守恒变量及其空间导数的计算，粘性项的处理以及边界的处理等，并且应用 f o r t r a n 编写了程序，进行了算例验证，计算了亚音速剪切流和超音速前台阶 流动。 本论文对方腔流动，超音速激波边界层相互作用和圆柱绕流进行了数值模 拟，证明本文的方法具有计算较大速度范围流动的能力。 最后，本文对c e s e 方法进行了并行程序设计的初步研究。采用改进的 c 彤s e 方法和并行分区技术，针对非结构网格，发展了一套求解二维e u l c r 方程 的并行程序。对n a c a 0 0 1 2 翼型的亚跨音速流动和多段翼型复杂绕流进行了分 区并行求解，多区并行计算的结果与实验的结果吻合较好，表明了本文研究方法 的正确性和有效性。同时文中的分区方法能实现各处理器之间的负载均衡，有效 的节省计算时间。 关键词：时空守恒元解元方法，非结构网格，n a v i e r - s t o k e s 方程，数值模拟， 并行计算 西北土业人学硕士学位论文 a b s t r a c t a b s t r a c t t h es p a c e t l m ec o n s e r v a t i o ne l e m e n ta n ds o l u t i o ne l e m e n t ( c e s e ) m e t h o di s an e wn u m e r i c a lm e t h o df o r t h ec o n s e r v a t i o nl a w s i nt h i st h e s i s ，t h ec e s em e t h o di s s t u d i e d o nt h eb a s e s ，t h ec e s ep r o g r a mu s i n gu n s t r u c t u r e d 鲥di sd e v e l o p e dt o s o l v et h ee u l e ra n dn 撕* s t o k e se q u a t i o n s ，a n dt h es u b s o n i c ，t r a n s o n i c ，s u p e r s o n i c f l o wa n dl o ws p e e di n c o m p r e s s i b l ef l o wa r ec o m p u t e d f i r s t l y , t h ec e s en u m e r i c a lm e t h o d sa r ei n t r o d u c e d ，i n c l u d i n gt h ec o n s t r u c t i o n o fc o n s e r v a t i o ne l e m e n t s ( c e ) a n ds o l u t i o ne l e m e n t s ( s e ) ，t h ec o m p u t a t i o no ft h e c o n s e r v a t i o nv a r i a b l e sa n dt h es p a c ed e r i v a t i v e s ，t h et r e a t m e n to ft h ev i s c o u st e r m s a n dt h et r e a t m e n to ft h eb o u n d a r yc o n d i t i o n s b a s e do nt h ea b o v em e t h o d , t h e f o r t r a n p r o g r a mi sc o m p l e t e d t h e 跚b s o m cs h e a rf l o wa n dt h ef r o n ts t e pf l o wa l e u s e dt ov a l i d a t et h ep r o g r a m s e c o n d l y , i nt h i st h e s i s ，s i m u l a t i o n so fad r i v e nc a v i t yf l o w , as h o c k b o u n d a r y l a y e ri n t e r a c t i o np r o b l e ma n df l o wo v e rc y l i n d e r s a l ec a r r i e do u t t h ep r e s e n t a p p r o a c hc a nc a l c u l a t ef l o w sa tw i d es p e e dr a n g e f i m a l y , t h ep a r a l l e lp r o g r a mb a s e do nt h ec e s em e t h o dh a sb e e ni n v e s t i g a t e d t h ei m p r o v e dc e s em e t h o da n das u i t a b l ep a r a l l e ld o m a i nd e c o m p o s i t i o ns t r a t e g y h a sb e e nu s e dt od e v e l o pt h ep a r a l l e lp r o g r a mf o rs o l v i n gt h ee u l e re q u a t i o n s ，b a s e d o nt h eu n s t r u c t u r e dg r i d s o m es t a n d a r dt e s tc a s e st h a tc o n t a i n e dn a c a 0 0 1 2a i r f o i l i ns u b s o n i ca n dt r a n s o n i cf l o wa n dm u l t i - e l e m e n ta i r f o i li nl o ws p e e df l o wa r e c o m p u t e d c o m p a r e dw i l l le x p e r i m e n t a ld a t a , t h ee f f e c t i v e n e s sa n dr e l i a b i l i t yo ft h e m e t h o da l ev e r i f i e db yn u m e r i c a lr e s u l t s a tt h es a m et i m e ，t h ep a r a l l e ld o m a i n d e c o m p o s i t i o nm e t h o dh e r ee , a na c h i e v el o a db a l a n c eb e t w e e nc o m p u t i n gp r o c e s s o r s a n ds a v ec o m p u t a t i o n a lt i m ee f f e c t i v e l y k e yw o r d s ：c e s e ，u n s t r u c t u r e dg r i d n a v i e r - s t o k e se q u a t i o n s ，n u m e r i c a l s i m u l a t i o n , p a r a l l e lc o m p u t a t i o n s 两北工业大学硕十学位论文 符号表 q q 姨 q f e e l e p e _ f 只 e e a ，b 户 甜 1 ， p r m a r e e e r ， 符号表 守恒变量 守恒变量x 方向的导数 守恒变量y 方向的导数 守恒变量时问方向的导数 x 方向的无粘通量 无粘通量e 的x 方向的导数 无粘通量e 的y 方向的导数 无粘通量e 的时间方向的导数 y 方向的无粘通量 无粘通量f 的x 方向的导数 无粘通量f 的y 方向的导数 无粘通量f 的时间方向的导数 流通量矢量五和户的雅克比矩阵 密度 x 方向的速度 y 方向的速度 压强 静温 动力学粘性系数 马赫数 雷诺数 x 方向的粘性通量 y 方向的粘性通量 气体常数，对于空气且= 2 8 7 0 5 , i 嘛k ) 气体比热比，对于空气，= 1 4 i i i 西北j ：业大学硕十学侍论文 符号表 p c v c p p r 单位质量流体所具有的内能 等容比热，对于空气c = 7 1 6 5 j ( k g 足) 等压比热，对于空气c 。= 1 0 0 3 5 j ( k g k ) 普朗特数，对于空气p r = 0 7 2 时间步长 i v 西北工业大学学位论文知识产权声明书及原创性声明 学位论文知识产权声明书 本人完全了解学校有关保护知识产权的规定，即：研究生在校攻读学位期间论文工作 的知识产权单位属于西北工业大学。学校有权保留并向国家有关部门或机构送交论文的复 印件和电子版。本人允许论文被查阅和借阅。学校可以将本学位论文的全部或部分内容编 入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存和汇编本学位论文。 同时本人保证，毕业后结合学位论文研究课题再撰写的文章一律注明作者单位为西北工业 大学。 保密论文待解密后适用本声明。 学位论文作者签名：走整堇 7 0 0 7 m 弓月? 日 指导教师签名 哆年j 月7 日 西北工业大学 学位论文原创性声明 秉承学校严谨的学风和优良的科学道德，本人郑重声明：所呈交的学位论文，是本 人在导师的指导下进行研究工作所取得的成果。尽我所知，除文中已经注明引用的内容 和致谢的地方外，本论文不包含任何其他个人或集体己经公开发表或撰写过的研究成 果，不包含本人或其他己申请学位或其他用途使用过的成果。对本文的研究做出重要贡 献的个人和集体，均已在文中以明确方式表明。 本人学位论文与资料若有不实愿意承担一切相关的法律责任。 学位论文作者签名：走缝堇 2 0 0 7 年；月9 日 西北工业人学硕七学位论文第一章绪论 第一章绪论 1 1 计算流体力学的发展现状 计算流体力学，英文名称为c o m p u t a t i o n a lf l u i dd y n a m i c s ，简称c f d 。它利 用计算机和数值方法求解满足定解条件的流体动力学方程以获得流动规律和解 决流动问题的专门学科。作为一门多学科交叉融合形成的新兴学科，它是计算机 科学、流体力学和计算数学相结合的产物。二十多年来，随着计算机性能的飞速 提高以及数值计算方法的不断发展，计算流体力学技术上已经日趋成熟。与理论 流体力学和实验流体力学相比，计算流体力学是一门年轻的学科，但它有很强的 生命力，发展速度很快。目前，计算流体力学作为一种基本的研究方法和设计手 段已经逐步进入了航空、航天、船舶、气象、能源、建筑等许多应用研究领域。 目前仍有许多人们感兴趣的复杂流动问题难以解决或是结果不令人满意。造 成这一困难的原因除了受计算机条件的限制和计算的理论模型需要完善外，数值 方法的好坏是一个至关重要的因素。就传统的数值方法而言，尚存在一些不足之 处：( 1 ) 大多数传统数值方法都是针对微分方程来进行数值离散的，且一般都把 时间和空间方向分开来处理，这就很难保证本质上时间与空间统一的物理守恒 律。( 2 ) 除紧致、超紧致格式以外，一般数值方法的精度不是特别高，对于一般 流动的计算还可以，难以满足一些复杂流动计算要求。( 3 ) 传统计算格式的分辨 率一般都不是很高，而一般的高分辨率格式通常都需要特殊的限制技术( 如通量 限制和梯度限制等) 或者局部网格加密技术等。而且构造复杂，通用性不好，不 便于推广到高阶格式和多维情况。发展高精度的数值方法是未来几十年计算流体 力学发展的一个重要方向。近几十年来，人们提出了不少好的数值方法和格式， 然而这些方法都各有其自身的特点和适用范围。本文采用的时一空守恒元解元方 法( c e s e ) 就是一种全新的数值方法，该方法无论是从概念上还是从方法上都与 传统的数值方法( 有限元法、有限体积法，有限差分法等) 有所不同：首先，它 从守恒型积分方程出发，严格遵守物理意义上的守恒律，通过设立守恒元和解元， 使格式局部和全局都严格保证其物理意义上的守恒律。其次，在计算中它把时间 l 西北工业人学硕七学位论文第一章绪论 和空问同等对待，从本质上保证了时间与空间统一的物理守恒律。第三，它除利 用了简单的泰勒展开式外，无需任何其它的数值逼近技术，也不需要任何单调性 限制或特征技术，因此该方法具有格式构造简单、计算方便和对激波间断分辨率 高等优点。最后，该方法可直接推广到多维情形，无需采用维数分解或交替方向 技术。更具吸引力的是，该方法不仅可用来求解连续流动问题，而且可用来求解 含有激波等的不连续流动问题，其高分辨率的数值结果甚至比目i 狰广泛流行的某 些高分辨率格式的计算结果还好。正是由于上述这些独特的优点，该方法正越来 越受到人们的广泛关注。 1 2c e s e 方法的基本原理 时一空守恒元解元方法是s c c h a n g 在文献【1 】中提出的求解双曲型守恒 律方程的一种崭新的数值方法。经过十几年的发展，该方法已经开始用于一些复 杂流动的计算，而且它还被美国n a s a 列为第二代c f d 程序中的主要算法之一。 c e s e 方法是基于守恒律的积分形式构造差分格式，真正意义上的守恒律是 指质量、动量、能量在空间和时间上流通量的守恒，守恒律的数学表现形式为一 组积分形式的守恒方程。在假定流场的物理解是光滑的前提下，可以将守恒律的 积分形式变换为微分形式。但由于在数值求解中不可能无限多地增加所用的网格 点数，在边界层等物理解急剧变化的区域内这种物理解光滑的假定就很难被数值 解反映出来。而且这种光滑的假定实际上在流场中有激波存在时已经不能成立 了，在数值计算中所采用的网格宽度的尺度范围内，激波是一个不连续的间断面。 我们可以看到，从微分形式的守恒律出发构造的数值格式，如果不能严格保证流 通量的守恒性，必然存在着很大的缺陷。而c e s e 方法从积分形式的守恒律出 发来构造格式避免了上述现象的发生。 c e s e 方法，即时一空守恒元解元方法，首先构造了解元( s o l u t i o n d e m e n t ， 简称为s e ) 和守恒元( c o n s e r v a t i o ne l e m e n t ，简称为c e ) 。解元是指某个网格点 附近的影响区域，在这个区域内的参变量和流通量向量分布可以表示为该网格点 上的参变量和流通量向量的一种简单的光滑函数。如可以用该网格点上的参变量 和流通量向量进行t a y l o r 展开得到。以下用s e 表示解元。守恒元是指满足守恒 律的区域，在这个区域内质量、动量、能量均满足守恒律，其概念类似于有限体 2 两北工业大学硕士学位论文 第一章绪论 积法中的控制体，不过在这里是整个时空域上的控制体。以下用c e 表示守恒元。 各个网格点的解元要求不能有交叉的区域，否则在交叉区域内物理量会发生 间断。然后将整个计算域划分成互不重叠的守恒元，并保证流通量在由任意一个 守恒元组成的区域内保持守恒。守恒元要求充满所要求解的时空域，这样就可以 保证从局部( 单个守恒元) 到整体( 整个计算域) 均满足守恒律。 在c e s e 方法里，时间和空间被完全统一起来同等对待。在传统的计算格 式中，在时间推进时，空间和时间通常是分别对待的，时间项通常作为一个独立 的未知变量出现，而在空间方向上将整个方程进行离散。这样，计算中的控制元 仅仅是空间计算域上的子域，在其上仅进行了空间离散。而实际上流体的质量、 动量、能量的流通量的守恒性在本质上是在空间一时间上的性质，这不能不说是 传统数值方法的一个不足之处。而在c e s e 方法里，将时间和空间紧密地结合 起来同等对待，该方法构造的解元和守恒元都是在空间一时间计算域上的子域。 在c e s e 方法里，将流场变量及其对空间的偏导数都作为独立变量同时进 行求解。在有限差分格式中，在网格点上的导数项是利用网格点上的气动参数通 过有限差分近似来得到的。当物理解在网格间距内变化不是很快时，这种近似的 精度是很高的，尤其在高精度格式中更是如此。但是在边界层内或激波附近等物 理解变化梯度较大的区域内，为了达到一定的精度，必须使用非常细小的网格， 这将导致计算量的大量增加。c e s e 方法在解元内将数值解表示成某类基函数的 展开形式，同时将展开式的系数看成是独立的待求解的未知数。这就避免了传统 数值方法的计算量增加的问题。c e s e 方法在解元中采用简单的t a y l o r 展开形式 表示数值解，并将t a y l o r 展开式的系数看成是数值解导数项的近似值。将导数项 作为未知量，通过求解过程直接得到。 在传统数值方法中，有限体积法是唯一试图强制流通量守恒的数值方法。但 在这一方法中，在两个相邻的控制单元交界面上，流通量是通过这两个控制单元 的网格点上的值经过插值得到的，所以计算结果依赖于插值模型的构造，并且不 可避免的带有一侧效应。c e s e 方法通过构造解元和守恒元解决了这个问题，即 在守恒元的交界面上，流通量不是通过插值得到，而是通过格式本身直接求解得 出。在传统数值方法中，有限差分法、有限元法以及谱方法都是求解微分形式的 守恒律。虽然在有限元法中所求解的是积分形式的守恒方程，但是在作了某些光 滑性假定后，实际上所求解的积分形式的守恒方程等价于微分形式的守恒方程。 3 西北1 ：业人学硕士学位论文 第一章绪论 c e s e 方法是个两步格式，随着时间的推进，网格点在空间方向上交错排列： 相应地，解元和守恒元也在空间方向上交错排列。这种处理方法使得在求解某未 知变量时，在任一个空间坐标方向上只用到与其相邻的两个解元，用传统的数值 方法的提法就是三点格式。在c e s e 方法中，格式精度的提高不受格式中所有 网格点数的影响，各种精度的格式都是三点格式。c e s e 方法通过提高解元内近 似解的t a y l o r 展开式的阶数来提高格式的精度。所以称c e s e 方法是一个两步 三点的通用格式。这种简单的三点格式模型与物理上流体流动的性质是一致的， 即除了体积力以外，流场某处流体的动力学特性只受与其相邻的流体的影响。 c e s e 方法的三点格式模型使得赋初值和边界处理简单、灵活、准确。在c e s e 方法中参变量的导数项在边界上也需要给出边界条件，这也给边界条件处理带来 了灵活性，并提高了边界处理的精度。 在传统的有限差分型高精度格式中，格式精度的提高的措施一般是增加离散 方程中的网格点数，如二阶t v d 或e n o 格式在每一空间方向上用到五个网格点。 五点格式在边界处理时会降低精度，除非在边界上使用与格式本身不一致的特殊 处理方法，而这又增加了计算、编程的复杂性。c e s e 方法随着t a y l o r 展开式阶 数的提高，所需求解的方程数目相应增加。如对于一维问题，一阶格式需求解两 个方程，而二阶格式则需要求解三个方程。 在传统的高精度格式中，在控制单元的交界面上，流通量一般是通过以特征 值为基础的插值方法得到的。而高精度格式一般是五点格式，所以要涉及到大量 的矩阵计算，使运行速度大大降低。在插值运算中一般又要用到特殊的处理技术， 如限定因子等，增加了格式的复杂性。c e s e 方法通过求解过程直接求得守恒元 表面上的流通量，避免了上述困难。 基于以上特点，c e s e 方法具有以下优点： 1 结构简单，局部性好，便于边界处理； 2 方法简单，通用性好，计算量小，便于推广到多维情形； 3 计算结果精度高，激波分辨率高； 4 除了利用简单的泰勒公式展开外，无需任何其它的数值逼近技术，也不 需要任何单调性限制或特征技术； 正是c e s e 方法具有以上独特的优点，它一经提出便引起了人们的广泛关 注，并且已经逐步开始用于一些复杂流动问题的计算。 4 两北f 业人学硕士学位论文第一章绪论 1 3 研究背景和研究意义 时一空守恒元解元方法是一种全新的高精度的数值研究方法，对翼型的高精 度设计有利，比与之相当的有限体积法的精度高，激波捕捉能力强，可以在一个 网格内精确捕捉激波以及反射波。求解相同的问题所需的网格数少，计算量小。 可以模拟一维、二维、三维定常和非定常流动，计算域宽，马赫数可以从0 0 0 2 8 模拟到1 0 。流动现象模型可以包括运动和反应激波，声波，涡脱落，激波和边 界层相互干扰作用，爆轰波，空化现象等。特别地，时一空守恒元解元方法有同 时求解强激波和小扰动( 声波) 的独特能力，使得它成为声压场计算的一种独特 的方法。二阶精度的时一空守恒元解元方法能够精确的计算航空声学问题，而其 他常见的方法的二阶格式不足以求解航空声学问题。在捕捉激波时必须加数值耗 散，而这样又抹煞了小扰动。所以一种求解方法如果能够同时处理强激波和小扰 动，必须克服上述困难，而时一空守恒元解元方法就具有这个能力。 在航空声学领域，喷流噪音是一个具有挑战性的话题。时一空守恒元解元方 法正好满足这方面的需求，它具有低耗散误差，准确捕捉激波等特点。而且无反 射边界条件的应用不会引入额外的误差，对产生噪音的旋涡可以精确模拟。而传 统的噪音计算方法，主要缺点是计算的精度很难达到噪音计算的要求，并且计算 量大，对旋涡噪音的模拟难以精确。 另外，多相流的模拟，超空化武器的优化设计是现在数值模拟的又一难题。 时一空守恒元解元方法可以模拟多相流，它引入了混合密度函数，发展了一种新 的空化模型，可以模拟翼型从空化到超空化的一系列变化过程。时一空守恒元解 元方法还可以模拟脉冲爆震发动机内的爆轰燃烧过程。 在国外，时一空守恒元解元方法已经应用在计算流体力学，计算声学和计算 电磁学等各个领域。基于时一空守恒元解元方法的计算机程序也发展到一维、二 维和三维流动。已经得到的数值结果包括：各种激波管问题，二维和三维e u l e r 和n a v i e r - s t o k e s 方程的计算，z n d 模型爆震波，内破裂和爆炸问题，前台阶激 波问题，喷管冲击波，各种声波，激波和声波的相互作用，二维翼型的超空化模 拟问题等等。在n a s a ( 美国国家航空航天局) 有一个专门的小组进行c e s e 方法研究。文献【1 _ 4 】是s c c h a n g 对c e s e 方法的格式最初提出和研究阶段。 在文献【5 】中，针对一维模型方程构造了a 一格式和a f 格式，并推广到一维 5 西北f ：业人学硕士学位论文第一章绪论 e u l e r 方程和一维n a v i e r - s t o k e s 方程。文献对一维格式进行了稳定性分析，并利 用一维e u l e r 方程的格式对七个激波管问题进行了计算，计算结果表明，c e s e 方法的一维格式精度很高，所捕获的激波仅仅占l 2 个网格，对接触间断和膨 胀波的数值模拟也相当准确。x yw a n g 在文献【6 】中发展了c e s e 方法的二 维格式。描述了二维格式的构造方法，构造了模型方程的各种格式，然后将其推 广应用到二维e u l e r 方程，并进行了稳定性分析。同时针对二维对流一一扩散方 程构造了时空守恒格式。并对c e s e 方法的二维e u l e r 方程格式在直角坐标系以 及曲线坐标系下的构造方法作了详细描述，同时附有大量算例。文献 7 】首次提 出了c e s e 方法的非结构网格的构造格式，并且将其推广应用到三维e u l e r 方程。 后面的研究工作多针对非结构网格开展的。文献 8 ，9 】将c e s e 方法应用到声学 领域，进行了一维喷管的声波传播的预测，并且模拟了声波和激波的相互作用， 计算结果和理论结果吻合很好，验证了c e s e 方法有求解航空声学问题的能力， 能同时捕捉声波和激波。文献【1 0 】描述了二维n a v i e r - s t o k e s 方程的非结构网格的 构造格式，粘性边界的处理，对激波与边界层的相互作用进行了计算，与实验结 果吻合很好。文献【1 1 】对超音速喷流噪音进行了初步探讨，文献 1 2 】阐述了三维 n a v i e r - s t o k e s 方程的格式。a h i m a n s u 等在文献 1 3 】中阐述了c e s e 的并行分块 计算。z e n g - c h a nz h a n g ，h a oh e 等在文献 1 4 ，1 5 q b 应用c e s e 方法对一维、二 维、三维爆轰进行了研究和计算，这为我们今后的计算提供了参考。 c e s e 方法在国内研究相对较少，主要集中在激波捕捉和爆轰计算。在1 9 9 7 ( 文献 1 6 1 9 1 ) 年国内出现了c e s e 方法的应用，但是仅仅到二维e u l e r 格式。 近几年( 文献 2 0 , 2 1 , 2 2 ) 张德良等对混合气体燃烧转爆轰( d d t ) 过程用c e s e 方法进行了数值模拟。现在国内公开发表的文章显示，主要研究了激波管问题， 二维欧拉方程，激波反射问题，二维n a v i e r - s t o k e s 方程的矩形边界网格的求解， 二维非定常爆轰计算等。 1 4 并行计算技术的发展状况 计算流体力学是并行计算的主要应用领域之一，随着计算流体力学的发展， 人们要解决更加复杂的问题，对计算机的计算能力要求越来越高，这也促进了并 行计算的发展。未来一段时期，并行计算将是计算流体力学发展的一大趋势。 6 西北t 业大学硕士学位论文第一章绪论 并行计算的含义就是在并行计算机或分布式计算机等高性能计算系统上所 做的计算，它的物质基础是高性能并行计算机，当然也包括分布式网络计算机。 早期的并行化技术研究主要集中在向量机的向量化技术，侧重于开发程序的细粒 度并行性，目前这种技术已经发展得相当成熟。为了进一步提高并行计算的计算 效率，从八十年代中期，人们开始侧重于研究和开发粗、中粒度的并行性。到九 十年代，并行计算已经成为国内外的研究热点。近几年，国内的研究所和大学也 逐渐发展起来一批高性能计算中心。 并行计算需要网格并行软件的支持。目前最具代表性的并行软件是并行虚拟 机p v m d l l ( p m a l l e lv k t u 且lm a c l l i ) 和消息传递接口m p i d 2 ( m e s s a g ep a s s i n g i n t e r f a c e ) ，这两种并行软件均支持异构的网络并行系统，它们具有形式简单，通 用性强，发展较成熟，源代码开放的特点。 m p i 是目前最为重要的并行编程工具，它具有可移植性好，功能强大，效率 高等多种优点，有多种不同的免费、高效实用的实现版本，几乎所有的并行计算 机厂商都提供支持。本文就是采用m p i 进行数据传递，完成并行计算。 在国外，有关c e s e 格式的并行计算已经有研究。文献【1 3 1 应用并行分区技 术，实现了程序的并行化，并且对并行效率进行了研究和分析。国内对c e s e 方法研究相对较少，尚没有c e s e 格式并行方面成果介绍。 1 5 本文的主要工作 本文的目的是研究高精度格式c e s e 格式，对于c e s e 格式在二维e u l c r ， n a v i e r - s t o k e s 方程中的应用以及c e s e 并行计算进行初步研究和探索，研究的 主要内容包括： 1 利用f o r t r a n 程序设计语言进行程序设计，并编写串行程序，应用 c e s e 格式实现了二维e u l e r 方程的求解。 2 应用c e s e 格式实现了二维n a v i e r - s t o k e s 方程的计算求解；并且进行 了算例验证。 3 基于非结构网格，采用改进的c e s e 格式和并行分区技术，发展了求解 二维e u l e r 方程的并行程序；并对并行加速比和并行效率进行分析。 7 西北t 业人学硕士学位论文 第_ 二章c e s e 数值方法简介 第二章c e s e 数值方法简介 2 1 一维e u l e r 方程的时一空守恒格式 2 1 1 格式构造 我们考虑一维守恒型e u l e r 方程： 害+ 罢：of ：删(2-1) 瓦+ 素_ o 2 删 其中口为常数，即对流方程的对流速度。令 = x ，而= t 为二维e u c l i d e a n 空 间易中的坐标，如图2 - 1 所示。 _ 定义：h = ( 删，“) ，则 图2 - 1 二维e u c l i d e a n 空间 v z ：口塑+ o u ：0 缸a 在易空间中运用高斯散度定理，可得( 2 1 ) 的积分形式的守恒律： 辱d ；= ( v ；沙= 0 ( 2 - 2 ) s ( r ) s 8 西北工业大学硕士学位论文第二章c e s e 数值方法简介 其中：( 1 ) h = ( a l ，”) 为e 2 空间的流通量向量，( 2 ) s ( v ) 为e 2 空间中任一 区域v 的边界，( 3 ) d s = d r y 刀为指向s ( y ) 外侧的法向微元表面的面积，d t r 为 s ( v ) 上的微元表面的面积( - - 维时退化为弧长) ，n 为指向s ( 矿) 外侧的单位法向 量。而d j 即为向量h 通过面积元d j 离开区域v 的流通量，我们可以看到在方 程( 2 2 ) 中，空间和时间被统一起来同等对待，这是c e s e 方法与大多数传统方 法的不同之处，也是c e s e 方法的基础。 下面介绍守恒元和解元的划分方法。如何设计守恒元( c e ) 和解元( s e ) 是c e s e 方法的关键之一。不同的设计方法将导致不同的计算格式。所谓的守 恒元是时一空选定的某一区域，在其中要求积分方程( 2 2 ) 成立。为此可把整个求 解区域按一定要求划分成许多个相互搭接而不重叠的守恒元，并要求积分方程 ( 2 - 2 ) 在其上满足，由此导出的若干个离散方程用来求解网格点上的流动变量及其 空间导数。这样做可保证积分方程( 2 - 2 ) 在局部( 即各个守恒元上) 和全局( 即所 有守恒元之和) 均得到满足。所谓解元则是指在每一个网格点上附近取定的一个 微小区域，在其中假定流动变量足够光滑并可用某一简单函数来近似。在原始的 c e s e 方法中，在每一个网格点上要求设置守恒元的个数与待求变量的个数相 等，以便得到相应数量的离散方程来求解这些变量。因此，在一维情形下，如果 每一个网格点上有两个独立的待求变l u 。和“。，那么c e s e 方法中在每一网格 点上需要设置两个守恒元。在这里，作者采用经张增产，沈孟育改进过的c e s e 方法【堋。该方法无论在一维，二维还是在三维情形，在每一个网格点上都只需设 置一个守恒元，由此得到一个主要的离散方程，而采用离散变量在相邻时间半层 解元的公共点上连续的要求( 或简单的采用差商逼近微商的方法) 来提供足够数 量的补充关系式。 具体讲，对于一维情形，在x t 平面上采用如图2 2 的交替网格，用q 代表易 中的网格点( j ，n ) 集合，其中j 和n 分别代表空间方向和时间方向的标号，n = 0 ， j 1 ，吾；对于每一个n ，i 玟j - - - n _ + j 1 ，n 3 2 ，- - 。( ，玎) q ，设有一个解元 s e ( j ，n ) 和一个守恒元c e ( j ，n ) 与之对应。这里，相应于网格点a ( j ，n ) 的解元取为 图2 3 中的四边形b d f g 所围区域，守恒元c e ( j ，n ) 取图2 3 中矩形b c e f 所围 区域。 9 两北工业人学硕十学位论文第二章c e s e 数值方法简介 j - 1 j j + l 图2 - 2 一维网格单元设置 g n n + l n n 一1 图2 3 解元s e ( j ，n ) 和守恒元c e ( j ，n ) 的划分 根据上面关于解元的假定，( x , t ) s e ( j ，疗) ，我们可采用不同形式的函数 ( x ，t ；j ，行) 和厂 ，t ；j ，n ) 来逼近方程中u ( x ，f ) 和f ( x ，r ) 。这里采用简单的一阶 t a y l o r 展开式来逼近u ( x ，f ) 和f ( x ，f ) ，即： 甜( 工，t ；j ，甩) = 甜；+ ( “，) ：( x 一) + ( 工，) ：o t 4 ) ( 2 3 ) 厂( x ，t ；y ，刀) = + c 疋) n ，( x x ) + ( z ) ：( f t ”) ( 2 - 4 ) 相应地可设 i l ( 五f ；，r 1 ) = ( f ( f ；，力，( 毛t , j ，竹) )( 2 - 5 ) 于是方程( 2 2 ) 便可以用下式来逼近 c 【h d s = 0 ，u ，行) q ( 2 6 ) 占( c f ( ，川口 一 现将4 = u ( x ，t ；j ，刀) 幕l f = f + “f ；工珂) ( 其中甜( x ，f ；_ ，n ) 和厂 ，t ；j ，珂) 分别 f l q ( 2 3 ) 和( 2 4 ) 决定) 代入式( 2 - 1 ) 可得 ( u t ) ：= 一( 正) ： ( 2 。7 ) 1 0 西北【业大学硕十学位论文 第二章c e s e 数值方法简介 由式( 2 7 ) 及简单分析不难看出，在每个网格点上要求的独立变量只有和 村，o 将式( 2 - 3 卜( 2 5 ) 代入式( 2 6 ) ，整理可得： 俨( 吉弘爿静嘶n - ：2 + s ，- t 圹2 踹纠 ( z - 8 ) 其中： 掣= 等沁牡石a t n + 篆 公式( 2 8 ) 提供了一个主要的离散方程，还需得到求”，的补充关系式。 为了得到求虬的补充关系式，要求“在两个不同时间半层的相邻解元的公 共点b 和f ( 如图2 3 ) 上连续( 或相容) ，分别可得： 彬一( 致) ：三缸= 一n ，- 1 1 ：2 + “) ：= ：；三f “；+ 饥) 磅缸= 哳n - 1 ：2 地崩n 1 1 ，：2 j 1 f 由此可得 ，) ：= ；) ：+ ：) ：】( 2 1 0 ) 其中： = 丢阮旷“力 ( 2 1 1 ) ( 二) 知：= “= 篮+ 等( ) 删n - 1 1 ：2 ( 2 - 1 2 ) 可见，离散方程式( 2 8 ) 和( 2 1 0 ) 就是我们得到的一维e u l e r 方程的时空守恒 格式。 对于不光滑或变化比较剧烈的流场计算，应对式( 2 8 ) 作适当的修正。我们可 以采用加权平均来替代其右端的简单的算术平均。 ，) ：= ( ；) ；， ；) ；，口)( 2 一x 3 ) 式中w 是一个限制函数，定义如下： ( 1 ) 矽( 0 ，o ，口) = 0 ( 2 1 4 ) ( z ) 矽( t ，口) = 紫当k 1 4 + i t l 4 。时 ( 2 一s ) 西北丁业人学硕十学位论文第二章c e s e 数值方法简介 这里，口是一个可调参数，一般取1 或2 即可。 自此，我们完成了一维c e s e 方法用于求解e u l e r 方程的格式的推导，对于 它的稳定性条件这里要求c f l l 即可，可参考文献【3 2 】。 2 2 二维e u l e r 方程的时一空守恒格式 将上面的一维e u l e r 方程的时一空守恒格式推广到二维情形 6 1 ，x yw a n g 在文献【6 】中发展了基于非结构三角形网格的c e s e 方法的二维格式，本节把 参考文献 2 5 】附录d 中c e s e 基本求解格式与参考文献【2 6 】中改进的二维c e s e 格式综合起来进行应用。 与结构化网格相比，复杂外形的非结构网格生成相对容易，而且，非结构网 格的存储方式，很适合c e s e 方法的求解格式，能体现c e s e 方法的优越性。 本节介绍基于非结构三角形网格的c e s e 格式。 2 2 1 二维非定常流动的e u l e r 方程的守恒格式 二维e u l e r 万程的微分形式为： 孥+ 譬+ 鲁- 0 ，：l 2 ，3 ，4 ) ( 2 - 1 6 ) 8 t瓠踟 ： 、 一。 + “ 其中，q l = p ，幺= p u ，q 3 = p v , 幺= p ( y 一1 ) + p ( u 2 + 1 ，2 ) 2 ， 局= q 2 ，如= 幺q q l ， 易= 一1 娩+ 【( 3 一y 炀一p l 妨j 2 q l ， 目= 鹏q 4 q l 一一1 埏妨+ 研j 2 骈， 互= q ，e = q q 3 q 1 ， 巧= p l 娩+ 【( 3 一y 妨一p l 砭扩2 9 ， e = 鹏q 4 q l 一( f - 1 ) q 3 眩+ gj 2 骈 遵照文献【5 1 的做法，g ，弘f ) 代表e u c l i d e a n 空间巨的坐标。应用g a u s s 散度定理 西北r 业大学硕士学位论文 第二章c e s e 数值方法简介 - 3 - 知，微分方程( 2 一1 6 ) 可写成如下积分型守恒律方程式 巩d s = 0 ，肌= 1 ，2 , 3 ，4 ( 2 1 7 ) 式中：s ( v ) 是e 3 空间中任一个区域的边界5 以= 娩，瓦，巴) ，m = 1 , 2 ，3 ，4( 2 - 1 8 ) 2 2 2 非结构网格c e s e 的构造 图2 - 4 c e s e 非结构网格 1 3 西北1 = 业人学硕士学位论文第二章c e s e 数值方法简介 ( b ) 图2 - 5 ( a ) 守恒元c e ( j ，n ) ， 解元s f a i 舢 现在我们在非结构三角形网格上应用c e s e 方法求解二维e u l e r 方程。图 2 _ 4 是一个任意三角形a v i v 2 v 3 和它的三个相邻的三角形。在图2 - 4 中点c ，c l ， c 2 ，c 3 分别是中心三角形和三个相邻三角形的质心。c a - - 个相邻三角形的质心 和中心三角形的顶点组成了一个六边形( c i v 3 c 2 v , c 3 v 2 ) 。中一b - - 角形的解点就是 六边形的质心s 点。每一个三角形都对应单一的解点，在图2 - 4 中用交叉点来表 示，如s ，s l ，s 2 ，s 3 。我们就是要求解解点处的守恒变量和它的空间导数。通 常情况下，解点s 与质心c 是不重合的，除非是均匀网格。六边形c l v 3 c 2 v i c 3 v 2 由三个四边形( c l v 3 c v 2 ，c 2 v l c v 3 ，c 3 v 2 c v l ) 组成。这三个四边形的质心分别是 e l ，e 2 ，e 3 ，如图2 _ 4 所示。 关于网格点j 在第n 时间层上的守恒元和解元示意图如图2 - 5 所示。上标( ) 和。分别表示n 和n + l 时间层。守恒元c e 0 ，n ) 定义为六边形柱面 c i v 3 c 2 v 1 c 3 v 2 c 1 v 3 c 2 v l c 3 v 2 。它可以分成三个子守恒元c e s ：c e r 0 ，n ) ，r = l ， 2 ，3 。每一个子守恒元是e 3 空间的一个四边形柱面。解元s e ( j ，n ) 是由三个竖直 t 4 西北t 业大学硕士学位论文第二章c e s e 数值方法简介 面和一个水平面组成的，如图2 - 5 0 0 ) 所示。 2 2 3 计算级 对于任一点瓯y ，0 s e 0 ，n ) ，通过下面的一阶1 a y l o r 展开用q ( x ，y ，n ；j ，n ) ， e y n ；j ，n ) 和f y ，n ；j ，n ) 来近似表示q y ，0 ，e ( 写y 0 ，f 伍y t ) ： q + k ) ，f ；，甩) = 研+ ( q z g 一_ ) + b ) ：( ) ，一乃) + ( q z ( f 一，) ( 2 - 1 9 ) e + g ，y ，t ；j , n ) = f + 佤r g 一_ ) + b z 一乃) + z ( f f ”) ( 2 - 2 0 ) f g ，y , t ；j , n ) = 曰+ 仁) ：g 一_ ) + k z 一乃) + 化片( f r ”) ( 2 - 2 1 ) 其中( z 代表坐标为b ，y j ，t “) 的解点o ，刀) 的相关量。彤和f 是关于研的 函数，佤z ，包) ；， ) ：，以z ，阮z ， e 是g ，缸z ，b z 和 z 的函数。 f = 砑，眩= a ，n 缸) ：，慨z = 协z ，佤z = 娩e ( 2 - 2 2 ) f = 彤研，以) ：= 彤缸z ，k z = 彤b z ，化) ：- b ， 娩)