基于复杂网格处理的高精度数值积分技术

林义淼

摘要：目前，对飞机表面进行网格划分并进行流体计算以优化整体结构的手段被广泛运用，然而由于外形复杂，采用笛卡尔网格划分时会出现带曲边或曲面的复杂网格，若直接进行数值积分则精度有限。该文提出了对复杂网格进行映射变换处理后再进行数值积分，并且数值实验表明采用该处理方法大大提高了精度。

关键词：网格划分；数值积分；复杂网格处理

中图分类号：TP393 文献标识码：A 文章编号：1009-3044（2015）29-0179-02

High-accuracy Numerical Integration Technology Based on Processing Complex Grids

LIN Yi-miao

（School of Mathematics and System Sciences， Beihang University， Beijing 100191， China）

Abstract： Currently， the method of aircraft surface meshing and computational fluid dynamics to optimize the overall structure is widely used， however， due to complex shape， there exist complex grids with curved edge or surface when using Descartes grid. The accuracy is limited if doing the numerical integration directly. This paper proposes that processing complex grids with mapping transformation before numerical integration， and the numerical experiment shows that using the method of processing complex grids can greatly improve the accuracy.

Key words： meshing； numerical integration； processing complex grids

随着数字化风洞的不断研究，使用计算机模拟飞机在风洞中的表现，大大降低实验成本来优化飞机结构，将飞机表面进行网格划分并且进行流体力学计算是在所难免的手段。

然而，由于飞机外形结构复杂，采用笛卡尔网格对飞机表面进行网格划分时，所出现的数值积分区域就不是标准的积分区域了。例如二维翼型网格划分后会出现带曲边边界的三角形和四边形网格，三维翼型网格划分后会出现带曲面的四面体网格。

目前，一般的数值积分技术中，对于带曲边边界的网格是直接线段化处理，然而这样做精度有限，而如果将网格加密处理，精度提高有限却会增加计算负担和降低计算效率，所以在本文将提出一种保留曲边边界进行处理并且提高了数值积分精度的方法。

1 复杂网格介绍

网格划分属于流体力学计算中的预处理部分，是将计算区域划分为较小的、不重叠的子域或单元（网格），从而得到有限单元及节点。而单元的数目决定了CFD的计算精度，网格的细密程度决定了数值解的精度和计算时间。

在采用笛卡尔网格进行划分时，由于飞机外形复杂，例如对飞机机翼表面进行网格划分时，会出现一边为曲边的四边形网格（或三角形网格）和两边为曲边的四边形网格（或三角形网格），而如今随着不断研究，飞机外形也是日趋复杂化，相应地对飞机表面进行网格划分所得到的网格也是日渐复杂。

当然，若是网格划分足够密，每个网格足够小，那么曲边在一定程度上也就相当于直边了，所以在目前的数值积分技术中，都是采用直接将曲边看成直线进行处理，但如此做法必然会影响到最后的数值积分精度，所以本文提出一种基于复杂网格处理的数值积分技术，来提高数值积分的精度。

2 复杂网格处理

现假设在二维平面[?2]上有一个带曲边边界的四边形网格，要求解函数[fx，y]在其积分区域[Ω]上的数值积分，一般情况下，直接将曲边看成直线边进行处理，利用上面介绍的数值积分方法（Simpson或Gauss方法等）求数值积分得到所要的结果，下面来介绍一下我们的处理方法。

首先建立一般化的直角坐标，坐标平面为[xOy]平面，四边形网格的四个顶点坐标所对应的向量依次为[γ1，γ2，γ3，γ4]，其中向量[γ2，γ4]所对应的顶点连接的正是网格的曲边。现有另一坐标系[O-ξη]，其上有一个单位为1的正方形区域[Ω]，接下来找到一个一一映射：

[φ：?2（ξ，η）??2（x，y）]

将区域[Ω]一一映射到区域[Ω]上。

很自然地，有如下映射关系：

[（0，1）?γ1，（0，0）?γ2，（1，1）?γ3，（1，0）?γ4]

那么就可以构造一个一一映射：

[γ（x，y）=γ1（1-ξ）η+γ2（1-ξ）（1-η）+γ3ξη+γ4ξ（1-η）]

则有：

[x=γ1（x）（1-ξ）η+γ2（x）（1-ξ）（1-η）+γ3（x）ξη+γ4（x）ξ（1-η）y=γ1（y）（1-ξ）η+γ2（y）（1-ξ）（1-η）+γ3（y）ξη+γ4（y）ξ（1-η）] （8）

根据数学分析中积分变换的内容有如下关系：

[Ωf（x，y）dxdy=Ωf（ξ，η）J（φ）dξdη]

其中[Ω]表示原来曲边边界网格的区域，[Ω]表示单位为1的标准正方形区域，而[J（φ）]为映射[φ]的[Jacobi]矩阵，

[J=?x?ξ?x?η?y?ξ?y?η]

我们还可进行如下处理：

根据一一映射关系，可以假设[?2（x，y）]中曲边的曲线函数为[Γ（ξ）（0≤ξ≤1）]，那么在曲边的两个端点处，就有如下关系：

[γ2=Γ（0），γ4=Γ（1）]

则在式（8）中就有

[x=γ1（x）（1-ξ）η+Γ（0）（x）（1-ξ）（1-η）+γ3（x）ξη+Γ（1）（x）ξ（1-η）]

其中

[Γ（0）（x）（1-ξ）（1-η）+Γ（1）（x）ξ（1-η）=（1-η）[Γ（0）（x）（1-ξ）+Γ（1）（x）ξ]]

通过对[ξ]的取值讨论，容易得到上式右边等于[（1-η）Γ（ξ）（x）]，于是式（8）就可以写成：

[x=γ1（x）（1-ξ）η+γ3（x）ξη+Γ（ξ）（x）（1-η）y=γ1（y）（1-ξ）η+γ3（y）ξη+Γ（ξ）（y）（1-η）]

对带曲边的网格区域进行如上处理后再利用已有的数值积分技术进行数值积分即可得到精度提高的数值结果了。

3 数值实验和结果

为了便于进行实验，我们取了一个相对简单的算例，并且能得到精确结果的例子，最后能够进行多方对比来验证我们的方法。

现假设有一个带曲边的四边形网格，并且曲边函数为[f（x）=ex]，另外三边分别为[x=0，y=0，x=1]，区域表示为[Ω]。

现在要求解[ΩF（x，y）dxdy]的数值积分结果，其中[F（x，y）=x2+y2]，我们都知道该结果完全可以精确求解出来：

[I=ΩF（x，y）dxdy=010exx2+y2dydx=19e3+e-199]

数值结果为[I≈2.8388970421]。

而在现有的数值积分技术中，都是忽略了曲边，直接连接了曲边的两点，而后利用现有的数值积分方法进行数值计算。我们选取了比较常用的Simpson方法和Gauss方法进行数值积分，得到结果如下：

[IS=2.7954089752，IG=2.8133023936.]

下面我们先用本文第三部分介绍的方法对曲边网格进行处理后有：

[φ：?2??2：（ξ，η）?（ξ，ηeξ）]

则[Jacobi]矩阵为

[J（φ）=?x?ξ?x?η?y?ξ?y?η=10ηeξeξ=eξ]

则所求的数值积分变成了

[I=ΩF（x，y）dxdy=Ω（ξ2+η2e2ξ）eξdξdη=Ωξ2eξ+η2e3ξdξdη]

其中[Ω]表示标准正面形区域。

之后我们依然采用现有的Simpson方法和Gauss方法，并且取同样的节点进行数值积分计算，得到的结果为：

[I′S=2.8384170008，I′G=2.8388960361.]

我们可以对以上数据进行简单的处理和误差分析，得到结果如表1所示。

其中[α]表示对于网格曲边直接连接曲边两点做线段化处理，而[β]表示对网格曲边进行映射变换处理，处理方法见本文第三部分。

从上表可以很容易看出，对曲边边界进行映射变换后再进行数值积分的话，精确度有显著提高，尽管这只是一个简单算例，但也说明了我们对于带曲边的网格进行适当的处理再数值积分的话，的确可以提高精度，相信这种处理方法和思想在高精度数值积分技术中会发挥其应有的作用。

参考文献：

[1] 徐利治，周蕴时.高维数值积分[M]. 科学出版社，1980.

[2] 徐利治，周蕴时，何天晓.高维数值积分选讲[M]. 安徽教育出版社，1985.

[3] Philip J.Davis， Philip Rabinowitz著，冯振兴，伍富良译.数值积分法[M].高等教育出版社，1986.

[4] Curtis F.Gerald， Patrick O.Wheatley.应用数值分析[M]. 吕淑娟，译.机械工业出版，2006.