11.1 结构网格

第 11 章 网格生成原理

在数值求解控制方程之前，必须离散所有边界表面，并在流动区域内部生成体网格。如第3章开头的3.1节所述，基本上可选择：

• 结构网格；

• 非结构网格。

图11.1和11.2给出民用飞机结构网格的示例 $ [1, 2] $；作为对比，图11.3和11.4展示类似构型的非结构表面网格与体网格 $ [3-6] $。

结构网格和非结构网格各有优点与不足，3.1节已有介绍。然而，无论采用哪种网格，从 CAD（计算机辅助设计）系统导入网格生成程序的数据质量，往往才是瓶颈。解析表面定义通常通过标准格式传递，例如 IGES（初始图形交换规范 [7,8]）、STEP（产品模型数据交换标准 [9–11]）或 DXF（图形交换格式 [12]）。另一种常用方法是通过 STL（立体光刻）格式 [13–15]，直接传递三角形表面网格；所有现代 CAD 程序均支持这一格式。直接使用 CAD 原生数据则相对少见 [16]。经验表明，传递过程可能降低数据精度。此外，CAD 系统自身的表面表示也常不够精确，通常导致相邻表面片之间出现缝隙、重叠或不连续。必须先消除这些缺陷，才能离散表面或利用表面生成体网格。这一过程称为“CAD 修复”[17,18]。

下面介绍生成结构网格（11.1节）和非结构网格（11.2节）的基本方法。受篇幅限制，这里仅作简要描述。关于表面建模，以及结构、非结构表面网格与体网格生成的深入讨论，见文献 [19,20]；网格生成方法的发展综述见文献 [21]。

11.1 结构网格

结构网格的显著特征是：物理空间中的网格点，可唯一映射到三个连续整数 i、j、k 的组合，每个整数对应一个坐标方向。这些整数组合定义了所谓的计算空间，见图3.2。计算空间中的坐标 $ \xi $、$ \eta $、$ \zeta $ 与 i、j、k 的关系为

$$ \begin{aligned}\xi&=i/i_{\max},\quad i=0,1,2,\cdots,i_{\max},\\\eta&=j/j_{\max},\quad j=0,1,2,\cdots,j_{\max},\\\zeta&=k/k_{\max},\ k=0,1,2,\cdots,k_{\max}.\end{aligned}\tag{11.1} $$

这种映射意味着 $ 0 \leq \xi \leq 1 $、$ 0 \leq \eta \leq 1 $ 和 $ 0 \leq \zeta \leq 1 $。连接相邻网格点，可在计算空间中形成立方体，在物理空间中形成六面体；二维情形则形成四边形。结构网格生成系统采用四边形离散流动区域的边界表面，称为表面网格，并以六面体填充内部。计算域内部的网格称为体网格。

生成结构网格时，首先沿边界曲线，即表面片的边界，布置网格点。通常在曲率较大的区域布置更密集的节点。利用边界曲线上的点分布，可生成表面网格；再根据包围物理区域的表面网格，构造体网格。因此，共同的问题是：如何依据已知的边界离散，生成区域内部的网格。可采用两种不同方法：

• 代数网格生成；

• 基于偏微分方程（PDE）的网格生成。

应用偏微分方程时，需要一个有效的初始表面网格或体网格，通常由代数方法生成。常用的偏微分方程有两类：

• 椭圆型方程；

• 双曲型方程。

代数网格生成（11.1.2节）通过直接的函数关系，描述计算空间与物理空间之间的坐标变换。最常用的代数技术是超限插值（Transfinite Interpolation，TFI）。给定全部边界上的点分布后，通过插值生成物理区域内部的网格点。某些 TFI 形式还允许控制边界处的夹角与网格间距。

基于椭圆型偏微分方程的方法（11.1.3节）最为常用。用户可指定网格线与边界的夹角，控制表面附近的网格间距及其增长比。椭圆型网格生成还能保证整个区域内的网格平滑，因此可生成高质量、边界正交的网格。其不足是，与代数或双曲型方法相比，计算时间长得多，而且还会遇到数值困难。

双曲型偏微分方程也可用于网格生成（11.1.4节）。这种技术沿某一个计算坐标方向，在两个表面之间逐步推进生成体网格。推进过程是显式的：依据已知的表面点分布，生成新的一层。双曲型网格生成的一个固有限制是，无法完全控制网格外边界的形状；但在外流等情况下，这未必构成实际问题。该方法可在整个区域内获得近似正交网格，而且计算成本低。

11.1.1 C 型、H 型与 O 型网格拓扑

在生成网格之前，必须先考虑其拓扑，即确定所需网格块的数量，以及各块之间的排列方式。对于复杂几何，这项工作可能耗时数周至数月。对于每个网格块，必须将计算域中的边界或部分边界，指定给物理空间中的某个边界，例如固壁或远场。物理空间中的网格形状很大程度上取决于这种对应关系。实践中形成了三种标准单块网格拓扑，称为 C 型、H 型和 O 型，因为在平面视图中，网格线类似相应的大写字母。三维网格可用两种拓扑的组合描述。例如，机翼周围网格通常在流向采用 C 型，见图11.2，在展向采用 O 型或 H 型；前者称为 C–O 型网格。下面详细讨论这三种拓扑。

C 型网格拓扑

C 型拓扑中，一族网格线包围气动外形，并在存在尾迹时同时构成尾迹区。图11.5给出示意图。可见，网格线 $ \eta = const. $ 从远场（$ \xi = 0 $）出发，沿尾迹经过后缘节点 b，顺时针绕过物体，最后重新延伸至远场（$ \xi = 1 $）。另一族网格线（$ \xi = const. $）从物体和尾迹沿法向伸出。网格线 $ \eta = 0 $ 的 a–b 线段构成坐标割线，即物理空间的 a–b 线段在计算空间中映射为两个线段，分别为 $ a \leq \xi \leq b $ 与 $ b' \leq \xi \leq a' $。因此，割线上侧（$ b' - a' $）与下侧 a–b 的节点在计算机内存中分别存储。相应边界条件见8.7节，网格示例见图11.6。

H 型网格拓扑

H 型网格拓扑常用于叶轮机械的叶片通道。图11.7展示这种拓扑。可见，气动外形的表面由两条不同网格线，即 $ \eta = 0 $ 和 $ \eta = 1 $ 描述。与 C 型网格不同，其中一族网格线（$ \eta = const. $）紧随流线，入口位于 $ \xi = 0 $，出口位于 $ \xi = 1 $。

乍看之下，似乎没有明显的坐标割线。然而，在叶轮机械中，线段 a–b 与 e–f 互为周期边界，三维情形为旋转周期；c–d 与 g–h 也如此。这类边界条件见8.8节。图11.8给出涡轮叶片之间非正交 H 型网格的示例。

O 型网格拓扑

由图11.9的 O 型拓扑可见，一族网格线 $ (\eta = const.) $ 在气动外形周围形成闭合曲线。另一族网格线 $ (\xi=\mathrm{const.}) $ 则沿径向分布于物体与外边界之间，此处外边界为远场。完整的边界线 $ \eta=0 $ 描述物体的轮廓，从 a 到 $ a' $。在计算空间中，坐标割线由边界 $ \xi=0 $（节点 $ a-c $）与 $ \xi=1 $（节点 $ a'-c' $）定义。图11.10展示一种标准 O 型网格，用于模拟翼型无黏绕流。O 型拓扑的不足是，在尖锐后缘处网格质量较差。

11.1.2 代数网格生成

依据指定的边界点分布生成表面网格或体网格时，最常用的代数技术是 TFI。Gordon 和 Hall 于1973年首次描述这一方法 $ [22] $。TFI 在计算空间各坐标方向采用一维单变量插值，其一般形式为

$$ \begin{gathered}\vec{U}=\sum_{i=1}^{L}\sum_{n=0}^{P}\alpha_{i}^{n}(\xi)\frac{\partial^{n}\vec{r}(\xi_{i},\eta,\zeta)}{\partial\xi^{n}},\\\vec{V}=\sum_{j=1}^{M}\sum_{m=0}^{Q}\beta_{j}^{m}(\eta)\frac{\partial^{m}\vec{r}(\xi,\eta_{j},\zeta)}{\partial\eta^{m}},\\\vec{W}=\sum_{k=1}^{N}\sum_{l=0}^{R}\gamma_{k}^{l}(\zeta)\frac{\partial^{l}\vec{r}(\xi,\eta,\zeta_{k})}{\partial\zeta^{l}}.\end{gathered}\tag{11.2} $$

式（11.2）中，$ \vec{U} $、$ \vec{V} $、$ \vec{W} $ 分别表示沿 $ \xi- $、$ \eta- $、$ \zeta $ 方向的单变量插值函数；$ \alpha_{i}^{n}(\xi) $、$ \beta_{j}^{m}(\eta) $、$ \gamma_{k}^{l}(\zeta) $ 为混合函数，$ \vec{r} $ 表示物理空间中网格点的位置。计算插值函数时，必须指定位置 $ \vec{r} $、导数 $ \partial^{n}\vec{r}/\partial\xi^{n} $ 等。由于边界曲线或表面已经离散，可将这些值代入式（11.2）。已知 $ \vec{U} $、$ \vec{V} $、$ \vec{W} $ 后，利用插值函数的 Boolean 和即可生成区域内部的网格，即

$$ \vec{r}=\vec{U}\oplus\vec{V}\oplus\vec{W}=\vec{U}+\vec{V}+\vec{W}-\vec{U}\vec{V}-\vec{U}\vec{W}-\vec{V}\vec{W}+\vec{U}\vec{V}\vec{W}.\tag{11.3} $$

式（11.3）保证二维情况下与全部四条边界曲线相符，三维情况下与全部六个边界面相符。式（11.3）的张量积按以下方式计算：

$$ \begin{aligned}\vec{U}\vec{V}&=\sum_{i=1}^{L}\sum_{j=1}^{M}\sum_{m=0}^{Q}\sum_{n=0}^{P}\alpha_{i}^{n}\beta_{j}^{m}\frac{\partial^{mn}\vec{r}(\xi_{i},\eta_{j},\zeta)}{\partial\eta^{m}\partial\xi^{n}},\\\vec{U}\vec{W}&=\sum_{i=1}^{L}\sum_{k=1}^{N}\sum_{l=0}^{R}\sum_{n=0}^{P}\alpha_{i}^{n}\gamma_{k}^{l}\frac{\partial^{ln}\vec{r}(\xi_{i},\eta,\zeta_{k})}{\partial\zeta^{l}\partial\xi^{n}},\\\vec{V}\vec{W}&=\sum_{j=1}^{M}\sum_{k=1}^{N}\sum_{l=0}^{R}\sum_{m=0}^{Q}\beta_{j}^{m}\gamma_{k}^{l}\frac{\partial^{lm}\vec{r}(\xi,\eta_{j},\zeta_{k})}{\partial\zeta^{l}\partial\eta^{m}},\\\vec{U}\vec{V}\vec{W}&=\sum_{i=1}^{L}\sum_{j=1}^{M}\sum_{k=1}^{N}\sum_{l=0}^{R}\sum_{m=0}^{Q}\sum_{n=0}^{P}\alpha_{i}^{n}\beta_{j}^{m}\gamma_{k}^{l}\frac{\partial^{lmn}\vec{r}(\xi_{i},\eta_{j},\zeta_{k})}{\partial\zeta^{l}\partial\eta^{m}\partial\xi^{n}}.\end{aligned}\tag{11.4} $$

关于网格生成中 Boolean 算子及张量积的更多细节，见文献 $ [23, 24] $。

可采用多种插值函数，例如线性、Lagrange、Hermite 和样条函数。最常用的是线性 TFI，即在式（11.2）和（11.4）中令 L=M=N=2、P=Q=R=0。这样，仅依据六个边界表面上的点分布即可生成体网格，其中令 $ \xi_{1} = 0 $、$ \xi_{2} = 1 $，对 $ \eta_{j} $ 和 $ \zeta_{k} $ 也作类似设置。线性 TFI 的混合函数为 [25]

$$ \begin{aligned}\alpha_{1}^{0}(\xi)&=1-\xi,\quad\alpha_{2}^{0}(\xi)=\xi,\\\beta_{1}^{0}(\eta)&=1-\eta,\quad\beta_{2}^{0}(\eta)=\eta,\\\gamma_{1}^{0}(\zeta)&=1-\zeta,\quad\gamma_{2}^{0}(\zeta)=\zeta.\end{aligned}\tag{11.5} $$

线性 TFI 的计算效率很高。图11.8给出采用线性 TFI 代数方法生成的网格示例。

采用三次混合函数的 Hermite TFI，还允许指定边界处网格线的斜率。令式（11.2）和（11.4）中 L=M=N=2、P=Q=R=1，即得到 Hermite TFI。也可在不同计算坐标方向混合使用线性插值与 Hermite 插值。Hermite TFI 及其他 TFI 扩展，例如网格间距控制，见文献 [20，第3章]。

11.1.3 椭圆型网格生成

基于椭圆型偏微分方程的网格生成方法，能产生单元尺寸及网格线斜率平滑变化的网格。它还允许控制边界附近的正交性与间距，这对模拟黏性流动尤为重要。Thompson 等人于1974年首次将椭圆型偏微分方程用于网格生成 $ [26] $。二维数值实现的详细描述例如见文献 $ [27, 28] $。

在三维情况下，对于网格点未知的笛卡尔坐标 $ \vec{r}=[x,y,z]^{T} $，Poisson 方程组可写为

$$ \begin{aligned}\alpha_{11}\frac{\partial^{2}\vec{r}}{\partial\xi^{2}}+\alpha_{22}\frac{\partial^{2}\vec{r}}{\partial\eta^{2}}+\alpha_{33}\frac{\partial^{2}\vec{r}}{\partial\zeta^{2}}+2\left(\alpha_{12}\frac{\partial^{2}\vec{r}}{\partial\xi\partial\eta}\right.\left.+\alpha_{13}\frac{\partial^{2}\vec{r}}{\partial\xi\partial\zeta}+\alpha_{23}\frac{\partial^{2}\vec{r}}{\partial\eta\partial\zeta}\right)\\&=-\frac{1}{J^{2}}\left(P\frac{\partial\vec{r}}{\partial\xi}+Q\frac{\partial\vec{r}}{\partial\eta}+R\frac{\partial\vec{r}}{\partial\zeta}\right),\end{aligned}\tag{11.6} $$

其中，P、Q、R 为控制函数。度量系数 $ \alpha $ 由以下关系给出：

$$ \begin{gathered}\alpha_{11}=\left(\frac{\partial\vec{r}}{\partial\eta}\cdot\frac{\partial\vec{r}}{\partial\eta}\right)\left(\frac{\partial\vec{r}}{\partial\zeta}\cdot\frac{\partial\vec{r}}{\partial\zeta}\right)-\left(\frac{\partial\vec{r}}{\partial\eta}\cdot\frac{\partial\vec{r}}{\partial\zeta}\right)^{2},\\\alpha_{22}=\left(\frac{\partial\vec{r}}{\partial\xi}\cdot\frac{\partial\vec{r}}{\partial\xi}\right)\left(\frac{\partial\vec{r}}{\partial\zeta}\cdot\frac{\partial\vec{r}}{\partial\zeta}\right)-\left(\frac{\partial\vec{r}}{\partial\xi}\cdot\frac{\partial\vec{r}}{\partial\zeta}\right)^{2},\\\alpha_{33}=\left(\frac{\partial\vec{r}}{\partial\xi}\cdot\frac{\partial\vec{r}}{\partial\xi}\right)\left(\frac{\partial\vec{r}}{\partial\eta}\cdot\frac{\partial\vec{r}}{\partial\eta}\right)-\left(\frac{\partial\vec{r}}{\partial\xi}\cdot\frac{\partial\vec{r}}{\partial\eta}\right)^{2},\\\alpha_{12}=\left(\frac{\partial\vec{r}}{\partial\xi}\cdot\frac{\partial\vec{r}}{\partial\zeta}\right)\left(\frac{\partial\vec{r}}{\partial\eta}\cdot\frac{\partial\vec{r}}{\partial\zeta}\right)-\left(\frac{\partial\vec{r}}{\partial\xi}\cdot\frac{\partial\vec{r}}{\partial\eta}\right)\left(\frac{\partial\vec{r}}{\partial\zeta}\cdot\frac{\partial\vec{r}}{\partial\zeta}\right),\\\alpha_{13}=\left(\frac{\partial\vec{r}}{\partial\xi}\cdot\frac{\partial\vec{r}}{\partial\eta}\right)\left(\frac{\partial\vec{r}}{\partial\eta}\cdot\frac{\partial\vec{r}}{\partial\zeta}\right)-\left(\frac{\partial\vec{r}}{\partial\xi}\cdot\frac{\partial\vec{r}}{\partial\zeta}\right)\left(\frac{\partial\vec{r}}{\partial\eta}\cdot\frac{\partial\vec{r}}{\partial\eta}\right),\\\alpha_{23}=\left(\frac{\partial\vec{r}}{\partial\xi}\cdot\frac{\partial\vec{r}}{\partial\zeta}\right)\left(\frac{\partial\vec{r}}{\partial\xi}\cdot\frac{\partial\vec{r}}{\partial\eta}\right)-\left(\frac{\partial\vec{r}}{\partial\eta}\cdot\frac{\partial\vec{r}}{\partial\zeta}\right)\left(\frac{\partial\vec{r}}{\partial\xi}\cdot\frac{\partial\vec{r}}{\partial\xi}\right).\end{gathered}\tag{11.7} $$

式（11.7）中，括号内的各项表示标量积。坐标变换 Jacobian 行列式的倒数 $ (J^{-1}) $ 按式（A.19）计算。令 P=Q=R=0 时，式（11.6）退化为 Laplace 方程。可利用其固有的光顺性质，改善代数生成网格的质量，但无法控制内部区域的点分布。

忽略式（11.6）中沿 $ \zeta $ 方向的导数，即可得到生成二维网格或表面网格的椭圆型方程：

$$ \begin{gathered}\alpha_{11}\frac{\partial^{2}x}{\partial\xi^{2}}-2\alpha_{12}\frac{\partial^{2}x}{\partial\xi\partial\eta}+\alpha_{22}\frac{\partial^{2}x}{\partial\eta^{2}}=-\frac{1}{J^{2}}\left(P\frac{\partial x}{\partial\xi}+Q\frac{\partial x}{\partial\eta}\right),\\\alpha_{11}\frac{\partial^{2} y}{\partial\xi^{2}}-2\alpha_{12}\frac{\partial^{2} y}{\partial\xi\partial\eta}+\alpha_{22}\frac{\partial^{2} y}{\partial\eta^{2}}=-\frac{1}{J^{2}}\left(P\frac{\partial y}{\partial\xi}+Q\frac{\partial y}{\partial\eta}\right).\end{gathered}\tag{11.8} $$

式（11.8）的度量系数为

$$ \begin{gathered}\alpha_{11}=\left(\frac{\partial x}{\partial\eta}\right)^{2}+\left(\frac{\partial y}{\partial\eta}\right)^{2},\\\alpha_{12}=\frac{\partial x}{\partial\xi}\frac{\partial x}{\partial\eta}+\frac{\partial y}{\partial\xi}\frac{\partial y}{\partial\eta},\\\alpha_{22}=\left(\frac{\partial x}{\partial\xi}\right)^{2}+\left(\frac{\partial y}{\partial\xi}\right)^{2},\end{gathered}\tag{11.9} $$

坐标变换 Jacobian 行列式的倒数 $ J^{-1} $ 则由式（A.25）定义。

要生成边界正交网格，必须在求解式（11.6）或（11.8）时指定合适的边界条件，基本上可采用 Neumann 条件或 Dirichlet 条件。Neumann 边界条件允许直接指定网格线与边界的交角，此时不需要控制函数 $ (P = Q = R = 0) $。但是，无法控制边界上的点分布及间距。实际上，边界节点会自动重新分布，以满足给定的网格线倾斜角。因此，该方法不适用于黏性壁面边界。椭圆型网格生成中 Neumann 条件的更多细节，例如见文献 [20，第6章]。

若边界点位置必须保持固定，则采用 Dirichlet 边界条件，并利用控制函数实现所需的交角与间距。图11.11展示二维情况下控制函数对网格的影响。可见，由于边界节点固定，式（11.8）中 P 为负时，网格线 $ \xi = const. $ 向左旋转；Q<0 时，则使网格线 $ \eta = const. $ 更靠近边界。首先，根据指定的倾斜角与间距同实际值之间的差，在相应边界上计算控制函数值；再将 P、Q、R 的边界值插值到区域内部，求解方程组（11.6）或（11.8）。重复这一过程，直至达到所需网格性质。基于 Dirichlet 条件的二维方法由 Sorenson [27] 及 Thomas 和 Middlecoff [29] 发展，后来由 Sorenson [30] 和 Thompson [31] 推广到三维。

椭圆型方程（11.6）或（11.8）通常采用二阶中心有限差分离散，所得线性代数方程组可用标准方法求解，例如采用多重网格加速的 Gauss–Seidel 松弛格式。控制函数值在外迭代中更新。为提高过程的稳健性，宜先用 Laplace 方程 $ (P = Q = R = 0) $ 进行若干次迭代，光顺通常由代数方法生成的初始网格 [32]。

图11.12展示椭圆型方程的作用。首先，利用上一小节的 TFI 方法生成翼型周围网格，结果见图11.12顶部；随后用 Laplace 方程光顺代数网格，即采用 P=Q=0 的式（11.8）。可见，Laplace 光顺改善了翼型前缘附近的网格，但不能控制表面法向间距。只有借助控制函数 $ P(\xi, \eta) $ 和 $ Q(\xi, \eta) $ 才能实现这一点，如图11.12底部所示。控制函数还使网格线与翼型表面正交。更多椭圆型网格生成示例见图11.1、11.2和11.6。

11.1.4 双曲型网格生成

当外边界形状无需精确控制时，双曲型偏微分方程适用于网格生成。首先指定初始点分布，然后沿法向，从上一层网格点推进给定距离，生成新的一层。将 双曲型偏微分方程用于网格生成，最早由 Starius $ [33] $ 及 Steger 和 Chaussee $ [34] $ 提出。较新的方法讨论见文献 $ [20, Chapter 5] $，各种扩展与改进见文献 $ [35–42] $。

生成体网格的双曲型方程为

$$ \begin{gathered}\frac{\partial\vec{r}}{\partial\xi}\cdot\frac{\partial\vec{r}}{\partial\zeta}=0,\\\frac{\partial\vec{r}}{\partial\eta}\cdot\frac{\partial\vec{r}}{\partial\zeta}=0,\\\frac{\partial\vec{r}}{\partial\zeta}\cdot\left(\frac{\partial\vec{r}}{\partial\xi}\times\frac{\partial\vec{r}}{\partial\eta}\right)=\Omega,\end{gathered}\tag{11.10} $$

其中，$ \vec{r}=[x, y,z]^{T} $ 表示网格点的笛卡尔坐标，$ \xi,\eta,\zeta $ 表示计算坐标，见式（11.1）；$ \Omega $ 为用户指定的单元体积。式（11.10）假设表面 $ \zeta= $ 为常数，代表初始状态。前两个关系表示正交条件，即 $ \xi= $ 为常数和 $ \eta= $ 为常数的坐标面，与 $ \zeta= $ 为常数的坐标面正交。最后一个关系保证单元体积等于 $ \Omega $，因此可控制各网格层之间的间距。

双曲型方程的二维形式可写为

$$ \begin{aligned}&\frac{\partial x}{\partial\xi}\frac{\partial x}{\partial\eta}+\frac{\partial y}{\partial\xi}\frac{\partial y}{\partial\eta}=0,\\&\frac{\partial x}{\partial\xi}\frac{\partial y}{\partial\eta}-\frac{\partial y}{\partial\xi}\frac{\partial x}{\partial\eta}=\Omega,\end{aligned}\tag{11.11} $$

其中，$ \Omega $ 现在表示指定的单元面积；初始点分布给定于曲线 $ \eta=0 $。

网格生成双曲型方程（11.10）或（11.11），在三维情况下沿 $ \xi $ 和 $ \eta $ 坐标、二维情况下沿 $ \xi $ 坐标，采用二阶中心差分离散。式（11.10）沿 $ \zeta $ 方向、式（11.11）沿 $ \eta $ 方向的推进，采用一阶隐式格式，即 ADI 格式，见6.2.3节。这样，推进步长只需根据所需网格间距选择。隐式算子用标准三对角求解器求逆。为稳定推进过程，必须在离散方程中加入人工耗散项，即二阶差分。与上一小节的椭圆型系统相比，双曲型方程（11.10）或（11.11）求解更快，通常也更容易。