数值流形法独立覆盖区域的一种自动选取方法

，,b

(长江科学院 a. 非连续变形分析实验室；b. 材料与结构研究所，武汉 430010)

1 研究背景

数值流形法(简称流形法)[1]采用2套相互独立的网格，即反映数值解精度的数学网格及反映几何边界和材料分区的物理网格，将研究区域划分成有限个相互重叠的覆盖，在各个覆盖上独立定义局部覆盖函数。覆盖函数可以灵活选取，如用可反映物理场局部特征的解析函数。为了保证连续性，在覆盖的重叠区域引入权函数，并使权函数之和为1。

目前，流形法多采用有限元网格作为数学网格，用有限元的形函数作为权函数。局部覆盖函数定义在结点上(以下称为结点覆盖)，每个有限元正好是其各个结点上的局部覆盖的重叠区域，即这几个覆盖在有限元的区域内完全重叠。这种方法存在如下一些问题：在结构的一些关键部位，如果2套网格的边界不吻合，其计算精度不高；局部网格加密不方便等。为了克服上述缺点，文献[2]提出了部分重叠覆盖的流形法，即采用大小单元相间的矩形数学网格，以大单元作为独立覆盖区域(即覆盖中没有与其他覆盖相互重叠的区域)，或一个大单元的所有结点覆盖合并为一个覆盖，小单元作为各相邻覆盖的重叠区域，以保证协调性。在此基础上，文献[3]提出了简单易行的局部加密技术。其实，在覆盖函数为非局部解析函数的情况下，如一般的多项式覆盖函数，重叠区域的精度高于独立覆盖区域，把重叠区域仅限于小单元显得不经济。此外，在三角形和四面体数学网格中，不易自动生成用作重叠区域的小单元。

为了改进独立覆盖区域的选取，本文采用常规尺寸单元的数学网格，应用图的着色方法将数学网格的单元(简称数学单元)染上不同的颜色，同色单元互不关联。在此基础上，尽可能多地将互不关联的数学单元选作独立覆盖区域，不属于独立覆盖区域的结点用结点覆盖(即常规的完全重叠覆盖)定义。这种独立覆盖区域的选取方法适用于任何单元类型的数学网格，而且选取过程完全由程序自动完成。

覆盖函数形式为高阶多项式函数时，完全重叠覆盖存在线性相关性(即整体刚度矩阵存在零特征值)。本文通过算例说明，线性相关仅存在于数学网格边界上的结点覆盖，覆盖合并可以减少整体刚度矩阵零特征值的数目，甚至完全避免线性相关。

2 数学单元的多色划分与独立覆盖的选取

图的着色是指对图的元素(如点、线、面等)染上一些颜色，使得同色元素互不相邻。有限元网格单元多色划分就是网格中的任意2个相邻单元都染上不同的颜色，使得同色单元互不相邻，成为一独立集。具有公共点、边、面的单元都属于相邻单元，而具有公共边、面意味着具有公共点。所以，只需以具有公共点作为判断单元相邻的依据，建立单元间的关系：相邻的单元数目和单元编号。

图的着色问题是著名的NP完全问题(多项式复杂程度的非确定性问题)，不存在多项式时间的精确算法。关于顶点着色问题的近似算法有很多种，本文采用贪心法(Greedy algorithm) 。该法所用颜色数不一定是理论上的最小值，但编程简单，速度快，并可满足本文的需要。

对于单元着色问题而言，贪心法的策略是：选择一种颜色，从任一单元开始，依次考察未被着色的每个单元，如果一个单元没有这种颜色的相邻单元，则用这种颜色着色；如果未被着色的单元都不能以这种颜色着色，则选择一种新颜色，再尽可能多地给这些单元着色。如图1所示的矩形网格，与一个结点关联的单元数最多为4，至少需用4种颜色着色。

图1 数学单元多色划分

经多色划分后，数学单元被分为若干个独立集。在此基础上，以独立覆盖区域最多为原则，选择用于定义独立覆盖的单元。首先将单元数最多的独立集内的所有单元选作独立覆盖区域；其次，对于剩余单元，依次选择一个单元，如与当前已有独立覆盖区域不相邻，则也选作独立覆盖区域；最后，将不属于独立覆盖区域的结点作为结点覆盖。如图1所示，单元数最多的独立集是蓝色单元集和红色单元集。如果选蓝色单元作为独立覆盖区域，则其他单元均与蓝色单元相邻，故不能选作独立覆盖区域。网格顶边的6个结点不是蓝色单元的结点，故定义为结点覆盖。这样包括6个蓝色单元的独立覆盖和6个结点覆盖，如图2(a)所示。

这种独立覆盖的选取方法适用于任何单元类型的数学网格。图2(b)中所示的三角形网格单元，可用6种颜色着色，因而有6个蓝色单元的独立覆盖和右边界上的2个结点覆盖。

图2 独立覆盖(蓝色单元)与结点覆盖(蓝色结点)

在独立覆盖的选取过程中，还可以人为事先设定某些独立覆盖区域或重叠区域。

3 部分重叠覆盖流形法的实现技术

在完全重叠覆盖流形法中，数学网格的每个结点上都有独立的局部覆盖函数。在部分重叠覆盖流形法中，局部覆盖函数定义在选作独立覆盖的数学单元或结点覆盖上。如图2所示，不同的蓝色单元，其局部覆盖函数不同，但同一蓝色单元的所有结点，其局部覆盖函数都相同。其他单元为重叠区域，其结点上的局部覆盖函数为所在蓝色单元的局部覆盖函数或顶边结点的局部覆盖函数。所以，只需增加结点自由度之间的关联信息，在完全重叠覆盖流形法程序的基础上，就可以实现部分重叠覆盖，实现过程简单方便。

部分重叠覆盖流形法的局部覆盖数目少于完全重叠覆盖流形法。如图2(a)所示，共有30个结点，但只需12个局部覆盖。然而，独立覆盖单元的权函数恒等于1，覆盖函数为多项式函数时，其精度不如重叠区域，所以覆盖函数的阶数要比完全重叠覆盖情况高1阶，部分重叠覆盖流形法的精度才能达到完全重叠覆盖流形法的精度。这时重叠区域具有更高的精度，将其选在精度要求高的部位是一种合适的选择。整体刚度矩阵中非零元素的个数是影响方程组求解计算工作量的因素之一。由于覆盖函数阶数不尽相同，2种流形法的整体刚度矩阵中非零元素个数的多寡，视具体情况而定。

刚度矩阵的性态是影响方程组迭代求解收敛速度的重要因素。目前，流形法的覆盖函数多为多项式函数，权函数(有限元的形函数)也为多项式函数。与其他单位分解方法一样，在高阶多项式覆盖函数情况，完全重叠覆盖流形法的广义自由度存在线性相关[4-6]，整体刚度矩阵是奇异的。从一维问题来看，在完全重叠覆盖流形法中，某结点的覆盖函数影响区域是该结点两侧的2个单元，该结点对应的权函数在2个单元中的分布形状是个三角形；在部分重叠覆盖流形法中，某独立覆盖单元的覆盖函数影响区域包括该单元及其两侧的2个单元，对应的权函数在3个单元中的分布形状是个梯形，干扰了线性相关性，有利于改善刚度矩阵的性态，其改善作用将在本文第4章中通过数值试验评价。

4 算 例

算例1为8 m×1 m悬臂梁，弹性模量为30 GPa，泊松比为0.167。左端固定，右端垂直向下作用力为100 kN。如图3(a)所示，均匀划分8个数学单元，物理网格与数学网格吻合。这种网格只需2种颜色着色，每种颜色4个单元。2种不同的独立覆盖区域布置分别见图3(b)和(c)，前者左端单元为独立覆盖区域，后者左端单元为重叠覆盖区域(即2个结点覆盖与其右边的一个独立覆盖构成重叠区域)。

图3 8 m×1 m悬臂梁网格与覆盖

平面应力问题的弹性力学解为，A点挠度=6.89 mm，B点拉应力=4.80 MPa。流形法的计算结果和整体刚度矩阵的零特征值个数见表1。从表1可见，完全重叠覆盖流形法采用1阶多项式覆盖函数，获得与弹性力学解答相当一致的计算结果；部分重叠覆盖流形法采用比完全重叠覆盖高一阶的覆盖函数，即2阶覆盖函数，左端单元为重叠覆盖区域情况的计算结果也与弹性力学解答基本一致，而左端单元为独立覆盖区域情况的计算结果与弹性力学解答相差较大，说明重要区域为重叠覆盖区域，有利于提高计算精度。

表1 算例1计算结果与刚度矩阵零特征值个数

从表1可见，对于本例，不同覆盖情况，整体刚度矩阵都存在零特征值，但是部分重叠覆盖情况零特征值个数均少于完全重叠覆盖情况。存在零特征值，表明刚度矩阵为半正定阵，用直接法求解方程组时，需采用高斯消去法或配以摄动法求解[4-5]。

笔者采用预处理共轭梯度法(PCG)求解这类方程，对不同类型、规模的算例都可顺利获得正常解。采用多项式覆盖函数的改进形式[7]，可以有效提高迭代收敛速度。

算例2如图4所示，为9 m×1 m悬臂梁，均匀划分9个数学单元，其他条件同算例1。算例3如图5(a)所示，为平面问题，划分9×9个数学单元，物理网格与数学网格吻合，底部全约束。

图4 9 m×1 m悬臂梁覆盖

算例2 和算例3的共同点是，在数学网格没有结点覆盖。计算结果表明，它们的整体刚度矩阵都不存在零特征值，即都是正定矩阵。文献[2]采用大、小单元相间的矩形数学网格，以大单元作为独立覆盖区域，小单元作为重叠覆盖区域，也不存在结点覆盖，可以类推，其整体刚度矩阵也是正定矩阵。

广义有限元法可以说是物理网格与数学网格吻合的完全重叠覆盖流形法。文献[4]将矩形区域划分为n×n(n=1,2,4,8,16,32,64,128)个单元，采用不同阶数P(1,2,3,4,5)的简谐函数作为局部覆盖函数，计算得到的广义有限元法整体刚度矩阵的零特征值个数与P有关，并且在P相同情况下，零特征值的个数随n线性增长。笔者采用多项式覆盖函数，将矩形区域划分为n×m个数学单元，无论是边界上存在结点覆盖的部分重叠覆盖，还是边界点和内部点都为结点覆盖的完全重叠覆盖，零特征值的个数都与边界上结点覆盖的个数呈线性增长。n×n网格的边界结点个数为4n，所以文献[4]所获得的零特征值个数随n线性增长的规律，实际上也是零特征值个数随边界结点个数呈线性增长。因此，笔者猜测，高阶多项式覆盖函数的线性相关性仅存在于边界结点，在边界上布置独立覆盖单元即可避免线性相关。为此，删去图5(a)网格内部的独立覆盖单元，如图5(b)所示，结果表明，其刚度矩阵的确不存在零特征值，即为正定矩阵。

图5 9×9网格2种不同覆盖布置

5 结 语

(1) 用图的着色方法选取独立覆盖区域，方便并可用程序自动实现。

(2) 部分重叠覆盖流形法可避免线性相关性或减少刚度矩阵的零特征值个数。

(3) 线性相关与边界上存在结点覆盖(即常规的完全重叠覆盖)有关。边界上无结点覆盖时，可避免线性相关，降低线性相关性应从边界点着手。直接对边界上的结点覆盖进行适当合并，或许也可以避免线性相关。

(4) 对于本文涉及的半正定线性方程组，基于优化思想的PCG法可以获得正常结果。

(5) 上述结论可供其他单位分解类算法，如广义有限元法参考。

参考文献：

[1] 石根华.数值流形方法与非连续变形分析[M]. 裴觉民译.北京: 清华大学出版社, 1997. (SHI Gen-hua. Numerical Manifold Method and Discontinuous Deformation Analysis [M].Translated by PEI Jue-ming. Beijing: Tsinghua University Press, 1997. (in Chinese))

[2] 祁勇峰，苏海东，崔建华. 部分重叠覆盖的数值流形方法初步研究[J]. 长江科学院院报，2013,30(1): 65-70. (QI Yong-feng, SU Hai-dong, CUI Jian-hua. Preliminary Study on Numerical Manifold Method with Partially Overlapping Covers[J]. Journal of Yangtze River Scientific Research Institute, 2013,30(1): 65-70.(in Chinese))

[3] 苏海东，祁勇峰. 部分重叠覆盖流形法的覆盖加密方法[J]. 长江科学院院报，2013,30(7): 95-100. (SU Hai-dong, QI Yong-feng. Cover Refinement for Numerical Manifold Method with Partially Overlapping Covers [J]. Journal of Yangtze River Scientific Research Institute, 2013,30(7): 95-100.(in Chinese))

[4] STROUBOULIS T，BABUSKKA I，COPPS K．The Design and Analysis of the Generalized Finite Element Method[J]. Computer Methods in Applied Mechanics and Engineering，2000，181(1-3)：43-69.

[5] 李录贤，刘书静，张慧华,等. 广义有限元方法研究进展[J]. 应用力学学报，2009，26(1): 96-108. (LI Lu-xian, LIU Shu-jing, ZHANG Hui-hua,etal. Research Advances in Generalized Finite Element Method [J]. Chinese Journal of Applied Mechanics, 2009，26(1): 96-108. (in Chinese))

[6] 刘 欣，朱德懋. 基于单位分解的新型有限元方法研究[J]. 计算力学学报，2000，17(4)：422-427.(LIU Xin, ZHU De-mao. Study on Partition of Unit FEM[J]. Chinese Journal of Computational Mechanics, 2000，17(4)：422-427.(in Chinese))

[7] 林绍忠，祁勇峰，苏海东. 数值流形法中覆盖函数的改进形式及其应用[J]. 长江科学院院报,2006，23(6):55-58.(LIN Shao-zhong, QI Yong-feng, SU Hai-dong. Improved Local Cover Function of Numerical Manifold Method and its Application[J]. Journal of Yangtze River Scientific Research Institute, 2006，23(6):55-58. (in Chinese))