计算流体力学—周正贵[(P111-P125)+P176]

1、7.3 网格生成图7.3（a）为直角坐标系下求解域局部网格，按照节点序号给出网格线对应的贴体坐标数值，和，。通常可简单地将节点序号作为在曲线坐标中的坐标值，比如。这样在计算域内即得到正方形网格，如图7.3（b）所示。如果在求解域内网格已分布完成，对应每一节点的x y坐标值已知，则坐标变换关系可用节点坐标值的差分表达式表示，比如点雅可比可表示为 图 7.3（a）求解域内局部网格 图7.3（b）计算域内局部网格 因此要得到坐标变换关系必须在求解域内分布网格线，即进行网格生成，而后由这些网格线构成曲线贴体坐标。合理的网格线分布不仅对计算精度有直接影响，甚至会影响计算过程的收敛性。因此，网格生成是流场

2、数值计算中一个比较关键的问题。通常要求生成的网格线要平滑，避免有局部太大的扭曲，不同簇网格线要尽可能正交。网格分布要与物理问题本身相匹配，也就是说疏密分布应与物理量变化率相适应。比如在求解黏性绕流时，在壁面附面层区流动参数沿壁面法向变化很剧烈，因此沿此方向网格线要密集分布。对于超音速流动问题在激波附近网格线要加密。网格生成主要有三种方法：代数方法、微分方程方法和保角变换方法。微分方程方法可以处理各种类型的不规则边界，具有较强的通用性，因而应用最为广泛；代数方法在工程实际中也有一定的应用；保角变换方法由于适用范围的局限性很大，目前实际应用较少。在此仅介绍代数方法和微分方程方法。7.3.1 代数生

3、成方法代数方法实际上是一种插值方法，下面以叶栅通道内H型网格生成来说明这一方法的基本思想。如图7.4，如果采用NS方程计算叶栅通道内流场，要求网格在叶片表面附面层加密。图7.4 叶栅绕流网格在叶片前后缘附近，由于流动参数变化较剧烈，因而此区域网格点分布也要适当加密。为此将流向网格线叶片通道内部分、切向网格线归为一类，此类网格线上节点分布要求二端密而中间稀。采用下列数学表达式对此类网格线上节点进行分布： （7.26） （7.27） （7.28） ， （7.29）分别为线段起始和终止坐标。为此线段上分布的节点数。分别为伸展和阻尼因子。当时网格节点在此线段上均布，当时网格节点在线段两端对称加密，越小

4、，越大，则二端节点密度越大，如图7.5所示。图7.5 网格节点分示意 将H型网格进口截面至叶栅通道进口和叶栅通道出口至H型网格出口截面流向网格线归为另一类，此类网格线上节点分布要求从一端向另一端节点分布由密而稀。对于此类网格线采用下述方法进行节点分布（以进口流向网格线为例）： （7.30） （7.31） （7.32） （7.33）分别为中弧线在前缘点坐标及导数值。为线段上节点数，为线段长度，为前缘处相邻节点间距离，为第i个节点坐标值。采用上式分布网格节点，可做到节点间距过渡光滑，在前缘点左右网格节点间距相等。7.3.2 微分方程生成方法微分方程生成方法根据所采用微分方程类型可分为双曲型方程、抛

5、物型方程以及椭圆型方程生成方法。其中椭圆型方程生成方法应用较为广泛，在此作详细阐述。如果一温度分布均匀的高温机翼置于静止流场中，在稳定状态下，四周液体由于热传导而呈不同的温度分布，如图7.6（a）画出温度分布等值线，可将这一簇温度等值线作为网格线（图中等线）。为求温度分布首先要给出控制方程。对于传热问题，有： 图7.6（a） 翼型绕流求解域网格示意 图7.6（b） 翼型绕流计算域网格示意当定常时，于是。在直角坐标系下展开，得或者写成因此可以采用上述椭圆型方程生成网格线。还可从液体运动角度分析。假设液体在由内边界HABCD和外边界GFE构成的“C”型流道内流动，液体不可压无粘无旋。这时有。如果采

6、用流函数等值线构成网格线，而势函数等值线构成网格线，则有。在直角坐标系下展开，得 （7.34）边界条件设定为 （7.35a） （7.35b）图7.6(a)和图7.6(b)反映了求解与计算域网格对应关系。利用边界条件（7.35），通过数值求解式（7.34）即可得到坐标变换数值关系。但采用式（7.34）在x y坐标系下计算，求解域边界不规则，因而需要将式（7.34）式进行变换，转化成在计算域内计算。利用式（7.19），（7.34）可转化成 （7.36）其中由于，故方程保持椭圆性。如果求解内边界HAB、CD、GFE、GH、DE4个边界形状确定，即每一个边界节点(x，y)坐标已知，则在计算域矩形边界上

7、每一边界节点x y值都为已知。采用式（7.36）变换是均匀的或者称为自然的，也就是说求解域内网格线分布仅与边界上的节点分布有关。事实上，如果加上一些辅助项还能使坐标拉伸或压缩，进而对网格线分布密度进行调整。下面就来说明这一点。考察方程它的解是，这就建立了坐标变换关系。显然这是一个线性变换，坐标上的等间距节点分布变换到x坐标上也是等间距的。如果方程改成它的解是。为了方便起见，设端点条件为。则有：在坐标上进行等间距节点分布，对应x坐标上已不是等节点间距了。当时，处的节点加密；当时则变稀疏，如图7.7。可见，改变c值可改变网格节点的疏密分布。如果希望在确定的所对应的附近改变节点分布的疏密程度，可以采

8、用下列方程：图7.7 网格节点分布规律 （7.37）其中： 这时如果式（7.37）中则对应的处附近网格节点分布变稀疏；如果则变密集。如果有n个节点附近网格疏密程度要调整，则可采用下列关系式： （7.38）当比较大时，由于指数函数衰减很快，各点间影响很小。将上述方法推广到二维问题中，得到： （7.39）其中： （7.40a） （7.40b）源项实现对线分布控制，其中（7.40a）中第1项实现对线与相邻网格线间距的整体控制图7.8（a）；第2项实现对线关于控制，如图7.8（b）。源项实现对线分布控制与源项实现对线分布控制相类似。图7.8 网格分布控制实际求解是在计算域内，因而要采用式（7.39）的

9、逆变换方程： （7.41）7.3.3 壁面处网格正交性分析 在任意线坐标系下，二维非定常可压缩流NS方程可写成： （7.42）将任一流动参数f的x，y偏导数表示成任意曲线坐标关系： （7.43）其中，J为坐标变换雅可比。如图7.9，考察壁面节点。采用一阶精度向前差分，采用二阶精度中心差分，即图7.9 近壁区网格对于叶栅通道流，如果采用H型网格计算，等线（叶栅切向）通常与y轴平行，这时。可推导出的离散误差表达式如下：由于，所以有： （7.44a）又有：。这时方程7.44（a）变为 （7.44b）同理可得： （7.44c）方程7.44（a）、（b）、（c）表明，流动参数y方向偏导数离散误差与网格线

10、夹角无关，而x方向偏导数离散误差与成反比。由于壁面附近流动参数变化剧烈，其二阶导数值相应较大，因此提高网格在壁面处正交性，可减小离散误差，提高计算精度。下例为平板紊流附层流场计算，考察网格正交性对计算结果影响。平板长度为0.06m，来流。计算网格由平行于平板和斜交于平板两簇网格线构成，斜交于平板网格线与平板交角分别为，与平板相邻网格线距平板约为2.0m。取距离平板前缘处垂直截面上计算结果进行比较。图7.10（a）、（b）分别为附面层内速度整体与近壁区局部分布比较。表1为壁面剪切应力，其中经验值为采用经验公式计算值，其中DE表示与经验值相对差。由图7.10（a）、（b）和表1知，当时，附面层内速

11、度分布及壁面剪切应力与经验公式符合很好。比较图7.10（a）和（b），角对近壁处速度分布影响更大。近壁处速度分布决定壁面剪切应力，由图7.10（b）和表1，要准确预测壁面阻力，不应小于。 图7.10（a）附面层内速度分布 图7.10（b）非常靠近壁面处速度分布表1 壁面剪切应力实验9080706050403048.9249.4449.5251.6058.7666.5678.2095.61DE(%)0.01.11.25.520.136.059.595.47.3.4 自适应网格简介在进行网格生成时要求在流动参数变化较剧烈的区域分布密集，比如激波附近、尾流区、附面层区等，因此在生成网格前要先估计在流

12、场中哪些区域流动参数变化较快。对于黏性绕流可以认定在物面附近的附面层区流动参数沿物面法向变化较快；但对于有些流动问题，比如超音速流，流场中有无激波、激波位置如何 往往很难预估。所以，在计算过程中，网格的分布最好能根据计算出的流动参数空间变化情况不断调整，这样就产生了自适应网格（Adaptive Meches）。调整的方法一般有两种：根据计算所得流动参数值，确定新的网格分布，用插值方法计算出新网格节点上的参数值；用动网格方法，网格坐标与流动控制方程联立求解。下面就一维问题阐述自适应网格的基本思想。引入权函数，使 （7.45）上式离散化的形式为 （7.46）其中。可见值取大，网格密度增加，变小；值

13、取小，网格变稀。引入计算域坐标，使，为常数，通常取网格点数，这时等于1。于是，式（7.46）变为 （7.47）上式中常数若为，x区间为(0，L)，对应的区间为(1，N)。则对应的方程为 对上式积分得 于是有： 因此， （7.48）以上积分中以为自变量，如果以x为自变量也可得类似的积分关系。由于 积分得 和 于是有： （7.49） （7.50a） （7.50b）以上给出的式（7.48）和式（7.50）即可用于生成一维自适应网格。它们之间不同的是式（7.48）中权函数；而在式（7.50）中。网格分布和流动控制方程求解可以分别进行，即每求若干步控制方程的解（这里假设问题是非定常的或迭代求定常解的问题

14、）之后，再根据生成网格的方程生成新网格。现在的问题是如何取值。最简单的方法是取 （7.51）为速度的一阶导数，于是根据式（7.47），有：这样在速度梯度较大的区域网格会自动加密。上式进一步可变成 （7.52）这就是说每一网格间距对应于相同的速度增量，即在坐标上是按等取网格点，如图7.11（a）所示。这种方法的=缺点是当时（在流动均匀的区域），网格间距。为此可将权重因子作一改变，取相等相等 （7.53）相等（a） （b） （c）图7.11 权函数与网格节点分布这时有 即 为微段弧长，于是式（7.47）变为 或 （7.54）也就是说网格间距由曲线等弧长条件确定，如图7.11（b）所示。权函数更一般

15、的取法是 （7.55）于是在关系曲线上：故 代入式（7.47），得或 （7.56）即在坐标上按等弧长取网格，如图7.11（c）所示。当时原来斜率大的地方变得更陡，网格加密程度更大。以上作法也有一些缺点，在有些问题中处曲率很大（即很大），如果希望在这些区域网格加密度，权函数可用下式： （7.57）其中k为曲线的曲率，为兼顾斜率和曲率的影响，权函数可以用 （7.58）或 （7.59）以上是一维自适应网格生成的基本方法，这种方法同样可推广到二维和三维自适应网格生成。7.3.5 计算网格生成实例卡门翼型绕流计算网格 采用椭圆型微分方程生成求解域中的网格，要将在计算域内的椭圆型方程（7.41）差分离散。如图7.12，令即（坐标值采用网格节点序号），得图7.12 局部网格示意图 （7.60a） 7.60b）其中： 方程（7.60a、b）可采用点迭代和线迭代方法进行求解，在此采用黎曼点迭代，得 （7.61a） （7.61b）以下采用上述差分方程生成网格，其求解域流场外边界为一圆，内边界为卡门翼型，如图7.13（a），其型面坐标数据见附录V卡门翼型x y坐标值。 图7.13（a）初始网格 图7.13（b） 生成的网格图7.13（c） 计算域网格 采用数值方法生成网格首先要给定边界节点x y坐标，并且还要给出假定初始时刻内部节点的x