《网格生成与坐标变换.ppt》由会员分享，可在线阅读，更多相关《网格生成与坐标变换.ppt（92页珍藏版）》请在三一办公上搜索。
1、计算流体动力学Computational fluid mechanics,机械与动力工程学院凌祥,2,第五章 网格生成与坐标变换,方程的一般变换度量和雅可比行列式再论适合CFD使用的控制方程注释拉伸(压缩)网格贴体坐标系:椭圆型网格生成自适应网格网格生成的进展有限体积网格生成的进展,3,方程的一般变换,考虑二维非定常流场,其自变量为x,y和t.我们要将物理平面中的自变量(x,y,t)变换成计算平面中的一组新自变量(,),用方程,目前这个变换是以一般形式写出来的.对于实际应用必须用具体的解析关系给出.,(5-1a),(5-1b),(5-1c),4,方程的一般变换,求导的链式法则,有,在上述表达式
2、中添加下标是为了强调求偏导数的过程中哪些变量保持不变,(5-2),5,方程的一般变换,在连续性方程,动量方程和能量方程中,未知函数是以导数形式出现,并且要将方程从(x,y,t)空间变换到(,)空间.根据求导的链式法则,有,(5-3),(5-4),(5-5),6,方程的一般变换,考虑第二章推导的粘性流动控制.可以看到方程的粘性项中包含二阶导数,这些二阶导数导数也要变换.在式(5-3)中,令,则有,(5-6),7,方程的一般变换,B和C是一个”复合导数”,它们包含了关于(x,y,t)中的某一个自变量和(,)中另一个自变量的导数,因此还需要进一步推导.由(5-2)规则得,(5-7),(5-8),8,
3、方程的一般变换,用式(5-7)和式(5-8)代替式(5-6)中B和C所代表的项,并重新排列各项的顺序,得到,式(5-9)中利用关于和的一阶偏导数,和二阶偏导数和混合导数乘以不同的度量给出了关于X的二阶偏导数.,(5-9),9,方程的一般变换,接下来,推导关于y的二阶偏导数.令,从式(5-2),有,其中,和,(5-10),10,方程的一般变换,利用式(5-4),有,(5-11),(5-12),式(5-11)和式(5-12)代回式(5-10),并重新排列各项的顺序,得到,(5-13),11,方程的一般变换,式(5-13)利用关于和的一阶偏导数,二阶偏导数和混合偏导数乘以不同的度量,给出了关于y的二
4、阶偏导数.再接下来,推导二阶混合偏导数,即,其中出现了与式(5-6)中B和C含义相同的项.,(5-14),12,方程的一般变换,用式(5-7)和式(5-8)代替式(5-14)中的这两项,并重新排列各项的顺序,得到,(5-15),式(5-15)利用关于和的一阶偏导数,二阶偏导数和混合偏导数乘以不同的度量,给出了 关于x和y的二阶混合偏导数.,13,方程的一般变换,例5-1 拉普拉斯方程,将方程(5-16)从(x,y)变换到(,),其中=(x,y),=(x,y),(5-16),14,方程的一般变换,利用式(5-9)和式(5-13),得,15,方程的一般变换,将上式合并同类项,最终得到,(5-17)
5、,考察方程式(5-16)和式(5-17).前者是物理平面(x,y)上的拉普拉斯方程,后者是计算平面(,)上经过变换的拉普拉斯方程.,16,度量和雅可比行列式,在式(5-3)到式(5-15)各式中,涉及网格几何性质的项,如,称为度量.在许多应用中,使用逆变换式(5-1ac)的逆变换更方便.,(5-18a),(5-18b),(5-18c),在式(5-18ac)中,和是自变量.,17,度量和雅可比行列式,为了用逆变换式计算方程中的度量,需要建立度量 等与逆度量 等的关系式.考虑流动控制方程中的一个未知函数,例如速度的x分量u.令,由变换式(5-18a)式(5-18b),即 和,u的全微分为,(5-1
6、9),18,度量和雅可比行列式,由方程(5-19),有,和,(5-20),(5-21),19,度量和雅可比行列式,式(5-20)和式(5-21)可以看作是两个未知数 和 的方程组,用克莱姆法则.从方程组式(5-20)和式(5-21)中解出,得到,(5-22a),20,度量和雅可比行列式,在方程(5-22)中,分母上的行列式称为雅可比行列式,记作,(5-22b),因此,将式(5-22)分子上的行列式展开,可以写成下面的形式,21,度量和雅可比行列式,(5-23a),回到方程组式(5-20)和式(5-21),解出,或,(5-23b),22,度量和雅可比行列式,为了能够进行一般性的讨论,将方程式(5
7、-23a)和式(5-23b)写成更一般的形式,和,(5-24a),(5-24b),23,度量和雅可比行列式,考虑二维的直接变换,(5-25a),(5-25b),对比式(5-25ab)与式(5-1ac),会发现从前面讨论中省略了=t.因为这里讨论只对空间度量感兴趣,所以忽略了时间的变换.,24,度量和雅可比行列式,根据全微分表达式,从式(5-25a)和式(5-25b)可以得到,或者用矩阵形式表示,(5-26a),(5-26b),(5-27),25,度量和雅可比行列式）,现在考虑逆变换,进行全微分,得,或者矩阵形式表示,(5-28a),(5-28b),(5-29a),(5-29b),(5-30),
8、26,度量和雅可比行列式）,从式(5-30)解出右边的列向量,也就是用式中2x2系数矩阵的逆矩阵去乘,得到,比较式(5-27)和式(5-31),得到,(5-31),(5-32),27,度量和雅可比行列式,根据计算逆矩阵的标准法则,式(5-32)可写为,(5-33),28,度量和雅可比行列式,考虑式(5-33)分母上的行列式.因为行列式转置后,其直不变,我们有,(5-34),29,度量和雅可比行列式,由式(5-22a)给出的J的定义可以看出,式(5-34)中的行列式正好就是变换的雅可比行列式J.将式(5-34)代入到式(5-33),得到,(5-35),30,度量和雅可比行列式,比较等式(5-35
9、)中两个矩阵的对应元素,就得到用逆度量表示的直接度量的关系式,即,(5-36a),(5-36b),(5-36c),(5-36d),右则公式给出直接度量与逆度量之间的关系式.,31,度量和雅可比行列式,把式(5-36a)和式(5-36b)代入到式(5-3)中,得到,这两个等式与逆度量给出的变换式(5-24a)和式(5-24b)完全相同,表明结果与前面的分析是一致的.,32,再论适合CFD使用的控制方程,在2.10节,方程(2-93)给出了流动控制方程的强守恒形式.对于空间二维的非定常流,如果没有源项,则方程化为,(5-37),简明起见,只考虑空间二维x,y的情形,而不是三维x,y,z的情形.但以
10、下的分析可直接推广到三维.,33,再论适合CFD使用的控制方程,问题:在计算平面(,)中,方程(5-37)还能写成守恒形式么?也就是说,变换后的方程还能写成,(5-38),这样的形式么,?,34,再论适合CFD使用的控制方程,首先,按照导数的变换式(5-3)和(5-4),对方程(5-37)中的变量进行变换,得,用式(5-22a)定义的雅可比行列式J乘以方程(5-39)得,(5-39),(5-40),35,再论适合CFD使用的控制方程,暂时放下方程(5-40),先考虑将 的导数展开,即,整理后,得到,(5-41),(5-42),36,再论适合CFD使用的控制方程,同样,考虑 对的导数,整理得,(
11、5-43),用同样的方法展开 和 并整理,得,(5-44),(5-45),37,再论适合CFD使用的控制方程,将式(5-42)式(5-45)代回到方程(5-40)并合并同类项,得,(5-46),38,再论适合CFD使用的控制方程,但发现,式(5-46)最后两项中用括号括起的部分等于零,其实,将式(5-36ad)代入到这些项中,就会得到,39,再论适合CFD使用的控制方程,于是,方程(5-46)可以写成,其中,(5-48a),(5-48b),(5-48c),(5-47),40,再论适合CFD使用的控制方程,如果将式(5-36ad)代入 和 的表达式(5-48b)和式(5-48c),可以得到,(5
12、-49a),(5-49b),其中 和 是用逆度量表示的.,41,注释,对于实际的问题和实际的几何形状,情况往往是:要么因为流动问题自身的特性(例如,流过平板的粘性流,壁面附近需要密集分布大量的网格点)要么因为边界的形状(例如,需要建立贴体曲线坐标系的弯曲物面)需要通过网格变换物理平面中的非均匀网格变换成计算平面中的均匀往格.有限体积法不需要这样变换,它能够直接处理物理平面中的非均匀网格.,42,拉伸(压缩)网格,例5-2 考虑图5-4所示的物理平面和计算平面.假设研究的是流过平板的粘性流.在平板的表面附近,速度迅速变化,如物理平面左边画出的速度剖面所示.,43,拉伸(压缩)网格,为了计算这种流
13、动在平板附近的细节,在y方向上需要使用细的网格,而在远离物面的地方,网格可以粗一些,如图5-4a.在计算平面上应建立均匀网格,如图5-4b,考察可以发现,在物理平面中,网格被”拉伸”了.,44,拉伸(压缩)网格,这种被”拉伸”的网格,可以用一个简单的解析变换就能完成,(5-50a),(5-50b),(5-51a),(5-51b),逆变换是,45,拉伸(压缩)网格,在物理平面中,x始终是同一个值.在计算平面中,也始终不变,所以在x方向上,网格并没有被拉伸.水平线就不是这样了,计算平面中的水平线是均匀分布的,始终不变,然而在物理平面中,相应的y值发生了变化,对求导,或者,用有限增量代替 和,近似地
14、得到,(5-52),46,拉伸(压缩)网格,例5-3 考察流动控制方程在物理平面和计算平面发生了哪些变化.为简单起见,假设是定常流,用连续性方程来进行说明.取定常流的连续性方程(2-25),在笛卡儿坐标系下,为,(5-53),这是物理平面的连续性方程.,47,拉伸(压缩)网格,利用导数的变换式(5-3)和式(5-4),可将这个方程变换到计算平面,变换后的形式为,(5-54),方程(5-54)中的度量可从直接变换式(5-50a)和式(5-50b)中求得,即,(5-55),48,拉伸(压缩)网格,将式(5-55)中的这些度量代入式(5-54),得到,(5-56),由式(5-50b),.因此,方程(
15、5-56)变为,或者,(5-57),49,拉伸(压缩)网格,例5-4 为了做进一步的论述,我们重作上面的推导过程,但这次是从逆变换式(5-51a)和式(5-51b)入手.回到方程(5-53),利用导数的逆变换式(5-24a)和式(5-24b),有,(5-58),50,拉伸(压缩)网格,方程(5-58)中的逆度量可从逆变换式(5-51a)和式(5-51b)中求得,结果是,(5-59),将式(5-59)代回式(5-58),得到,(5-60),51,拉伸(压缩)网格,例5-5 现在考虑更复杂的网格拉伸变换,研究绕钝体底部的超声速粘性流.此时,网格在x方向和y方向都要进行拉伸.图5-5给出了物理平面和
16、计算平面.根据变换,(5-61),52,拉伸(压缩)网格,流动,物理平面,计算平面,图5-5 均匀网格与压缩网格的比较,53,拉伸(压缩)网格,进行x方向(流向)的拉伸,其中,(5-62),(5-63),式(5-61)中,是物理平面发生最大聚集的点在计算平面中相应的位置.是常数,控制向 点聚集的程度,越大,聚集得就越密集.,54,拉伸(压缩)网格,网格沿y方向的横向拉伸,需要把物理平面分成:底部后面和底部的上方上下两部分.利用变换,(5-64),进行y方向的拉伸,其中,而 和 都是常数,对于刚才分出的两个部分,这些常数取值是不同的,55,贴体坐标系:椭圆型网格生成,考虑流过图5-6a所示的扩张
17、管道的流动.曲线de是管道上壁,直线f g是中心线.对于这种流动,物理平面中用矩形网格是不行的,在图5-6b中画曲线网格,使边界de和中心线f g都成为坐标线,曲线网格完全贴着边界.,图5-6 简单的贴体坐标系,56,贴体坐标系:椭圆型网格生成,曲线网格(图5-6a)变换成矩形网格(图5-6b)用下述方法:令 是图5-6a中上表面de的纵坐标则下述变换将生成(,)平面内的矩形网格,其中,(5-65),(5-66),57,贴体坐标系:椭圆型网格生成,例如,考虑物理平面的d点,.将这点坐标代入式(5-66),得,在计算平面中d点位于 上,考虑物理平面c点,此时,在计算平面中,c点也位于 上,与计算
18、平面中点d的坐标值相同.,58,贴体坐标系:椭圆型网格生成,更复杂的例子,如图5-7所示它是围绕机翼的流动设计的.,图5-7 椭圆方法生成的贴体网格,59,贴体坐标系:椭圆型网格生成,问题:那种变换能将图5-7a的曲线网格变换成图5-7b计算平面中的均匀网格?方法:把物理平面的网格用笛卡尔坐标表示,沿着内边界1,物理坐标已知,同样外边界2的物理坐标已知,这就暗示存在一个边值问题,求解椭圆型偏微分方程需要沿着包围区域的边界在每一点上给定边界条件.,60,贴体坐标系:椭圆型网格生成,最简单的椭圆方程是L a p l a c e方程,(5-67),(5-68),61,贴体坐标系:椭圆型网格生成,在方
19、程（5-67)和(5-68)中,和是因变量,x和y是自变量.将两组变量对调一下,写出逆方程,使x和y变成因变量,结果是,(5-69),(5-70),62,贴体坐标系:椭圆型网格生成,其中,63,贴体坐标系:椭圆型网格生成,图5-8画出了计算平面,再次强调沿着所有四条边界线1,2,3和4,(x,y)的值都已知,这是求解椭圆偏微分方程适定问题的基础.,图5-8计算平面(标出了边界条件并画出了一个内点),64,贴体坐标系:椭圆型网格生成,流动控制方程中的度量可以用有限差分计算,在(不论是物理平面还是计算平面)标记为(i,j)的网格点处,可以写出度量计算公式为,这样计算出来的度量值将被代入到变换后的流
20、动控制方程中.,65,贴体坐标系:椭圆型网格生成,实用的绕机翼型网格见图5-9,它就是用上述椭圆型网格生成的.,图5-9 围绕M i l e y翼型的网格,66,贴体坐标系:椭圆型网格生成,图5-10给出了翼型附近网格的细节.,图5-10 图5-9中贴体网格的一小部分,显示了翼型附近网格分布的细节,67,自适应网格,应该将大量的,密集的网格点分布在流场变量存在大的梯度的那部分流动区域内,从而改进CFD计算的数值精度.这样做不仅是因为密网格能够减小截断误差,更是因为要捕捉流动特性,梯度大的地方就需要更多的网格点.图5-11给出的流过平板的粘性流,可用来定性地说明这个问题.,68,自适应网格,图5
21、-11 在边界层内需要集中大量的网格点,69,自适应网格,自适应网格是能够自动向流场中大梯度区域聚集的网格,它利用求解的流场特征确定网格点在物理平面中的位置.自适应网格的优点是:当网格数量固定时,可以提高计算精度给定精度时,可以减少的网格点来达到这一精度.,70,自适应网格,自适应网格的一个简单的例子是C o r d a求解绕后台阶的粘性超声流时所用的网格,其中的坐标变换为,(5-71),(5-72),式中,g是流场中的原始变量,如p,或 T.,71,自适应网格,图5-12自适应网格原理示意图,72,自适应网格,(图5-12a)在物理平面中,t时刻N和N+1点的位置用黑色的圆点表示,两点之间的
22、x方向距离,其中上标t时刻.t时刻N点的x坐标,依赖于1,2点之间的,2,3点之间的,也就是说,(5-73),73,t+t时刻N和N+1点的新位置在图(5-12a)中用十字形记号标出,图中同时还标出了 x的新值.t+t时刻N点的x坐标为,(5-74),自适应网格,74,自适应网格,可以用N点和N+1点的x坐标的相应改变量除以时间增量 t来表示时间度量,即,(5-75a),式中,和 分别由(5-74)和式(5-73)给出.,(5-75b),75,自适应网格,其中,和 由类似式(5-73)和式(5-74)的表达式给出,即,76,自适应网格,回到一般的逆变换式(5-18ac).先考察式(5-18a)
23、,即,(5-18a),写出全微分,有,(5-76),方程中x的变化量d x由,和的变化量d,d 和d 表示.,77,如果这些变换是关于时间的,保持x和y不变,方程(5-76)可以写成,自适应网格,(5-77),在方程(5-77中),由于x始终保持不变,所以 恒等于零,78,不失一般性地,设式(5-18c)中的t=t()的形式给出,则在式(5-77)中 代入这些值,式(5-77)变为,自适应网格,(5-78),79,自适应网格,现在考虑方程(5-18b),即,(5-18b),于是,(5-79),从这个结果可有,(5-80),或,(5-81),80,自适应网格,式(5-78)和式(5-81)中都有
24、度量 和.将两式联立,用克莱姆法则从中解出.,(5-82),81,自适应网格,注意=t,并且分母就是雅可比行列式J,式(5-82)边为(省略下标),(5-83),用同样的方法,从式(5-78)和式(5-81)中解出,得,(5-84),82,自适应网格,出现在(5-83)和式(5-84)以及雅可比行列式J中的空间度量,和,则可以用中心差分来代替.例如,式中,i=N,j=M,83,自适应网格,图5-13给出了绕后台阶超声速流自适应网格的例子,流向从左至右.,图5-13 后台阶问题的自适应网格,84,网格生成的进展,图5-14用椭圆型网格生成,结合自适应网格,生成的三维贴体网格,可以看到机身表面,对
25、称平面和机翼平面内的网格线.对于现代CFD求解已经在使用由两个或多个独立的网格块组合而成的网格,网格块之间由交界面分开.也就是说,网格由两个或多个网格块组成,其中每一个网格块都独立与其他网格块.这些不同的网格块覆盖着流场的不同区域,因此被称为分块网格.,85,网格生成的进展,图5-14 覆盖F-20飞机外形的椭圆型自适应网格,86,网格生成的进展,图5-15 覆盖F-16飞机的分块网格,87,有限体积网格生成的进展,首先,将物理平面中的非均匀网格看成是有限体积法的网格单元.有限体积法不需要像有限差分法那样利用均匀的矩形网格计算,在物理平面内的非均匀网格上就可以直接进行有限体积计算.其次,有限体积法根本不需要结构网格,它可以适用于任意形状的网格单元,于是产生了非结构网格.,88,有限体积网格生成的进展,图5-16 环绕多元翼型的非结构网格,89,有限体积网格生成的进展,图5-17 压缩拐角上的非结构网格,90,有限体积网格生成的进展,正当非结构网格边得流行起来的时候,在完全相反方向上也取得了新的进展,那就是使用具有规整的笛卡尔网格.,图5-18 物面附近的笛卡尔网格,91,有限体积网格生成的进展,图5-19 计算多元翼型亚声速绕流的笛卡尔网格,92,有限体积网格生成的进展,图5-20 计算双椭球高超声速绕流的笛卡尔网格,
