基于Wolfe准则修正牛顿法的三维网格质量优化方法

本发明属于三维四面体delaunay网格生成的网格质量优化，涉及一种基于wolfe准则修正牛顿法的三维网格质量优化方法。

背景技术：

1、网格生成是数值计算的一个重要分支，特别是计算流体力学(computationalfluid dynamics，cfd)对网格质量的要求更高，作为cfd数值计算与模拟的第一步，也是最基础与最重要的一步，它的自动化程度与网格质量能直接影响cfd计算周期的总耗时、数值计算模拟的精度与效率。近年来，随着高性能计算机和cfd的飞速发展，cfd应用领域越来越广泛，问题也越来越复杂，现有的网格生成技术似乎已经不足以承担复杂外形的高质量数值模拟，因此具有高质量的网格生成技术仍然是目前的研究热点；在网格生成领域中，delaunay三角化法有强大完整的理论支撑，并且由于这种方法生成的网格在基础、效率、质量方面都比较优异，其已成为目前通用的全自动网格生成方法之一。

2、三维四面体有限元网格剖分较之二维三角形网格剖分，出现形状恶劣单元的几率要大得多，原因是四面体单元扭曲变形的形式比三角形单元要多。因此，三维四面体有限元网格优化是一个值得关注的问题。现有的网格优化算法可以分为三种类型：点删除点插入操作、网格拓扑变换、网格光顺。网格光顺一般又可以分为laplacian光顺、智能laplacian光顺和基于优化的光顺。

3、其中，基于优化的光顺利用了“最优化”理论，通过对目标函数求最优解的方法来确定每个网格节点的最优位置，因此在不考虑时空效率的情况下，理论上单独使用基于优化的光顺基本可以实现对初始网格最大限度的优化效果。但网格质量优化的整个过程一般会将各种算法结合在一起，在这过程中，基于优化的光顺也是最消耗计算时间和内存的步骤之一。

4、基于优化的光顺在网格优化中是指应用优化算法来改善和优化网格的拓扑结构和几何形状，以获得更平滑和自然的网格形态，其有几大核心问题，包括单元度量准则的选择，目标函数的建立，优化算法的实现等，其中优化算法的实现是一大难点，它能很大程度决定基于优化的光顺对网格的优化效率以及劣质网格的优化率，因此，一个好的优化算法是至关重要的。

5、一种现有的基于优化的光顺的优化算法(董亮,刘厚林,代翠,等.基于光顺的网格优化算法及其关键因素分析[j].江苏大学学报:自然科学版,2012,33(5):533-537.)，该算法利用数字1与半径比的差值作为目标函数，在每次迭代时会先求解目标函数在当前点位置下的梯度，然后与利用前后信息构造的多个矩阵相乘得到下降方向，进而计算点最优位置，达到优化效果。该方法的计算量比较大，每次迭代对上一次迭代结果的依赖性较大，导致最终优化效果不是很好。因此，需要构造一种能适应三维四面体网格特性的优化算法，在保证稳健性的情况下降低计算量的同时达到更高的优化率，最终实现高稳健性、高优化率的网格质量优化。

技术实现思路

1、针对上述存在问题或不足，为解决构造一种能适应三维四面体网格特性的优化算法，在保证稳健性的情况下降低计算量的同时达到更高的优化率，实现高稳健性、高效率的网格质量优化；本发明提供了一种基于wolfe准则修正牛顿法的三维网格质量优化方法。

2、一种基于wolfe准则修正牛顿法的三维网格质量优化方法，包括以下步骤：

3、步骤1、将目标飞行器进行初始有限元建模，进行网格剖分，建立三维四面体网格，得到网格信息；网格信息包括所有网格节点的坐标信息、模型表面三角形网格信息和模型内部四面体网格信息，其中模型内部四面体网格在后续的步骤中统称为网格单元。

4、所述三角形网格由3个不同的网格节点组成，四面体网格由4个不同的网格节点组成，每个网格节点拥有唯一且互不相同的编号。

5、步骤2、将步骤1得到的三维四面体网格标记为初始网格，并依据模型表面三角形网格信息将所有网格节点进行分类，其中参与构成模型表面三角形网格的网格节点称为边界网格节点，剩下的所有网格节点称为内部网格节点。

6、步骤3、确定网格单元质量度量准则，设定目标函数a(x)。

7、本发明使用的网格单元质量度量准则是半径比ρ，具体如下式

8、

9、其中v是四面体网格单元的体积，sm是四面体第m面的三角形面积，lj是四面体第j条边的边长，在该网格单元质量度量准则下，正四面体取得最大值1。

10、在进行网格质量优化前，会首先对每个待优化网格节点构造局部优化域，局部优化域的具体定义如下：

11、首先在当前网格中查找所有以待优化网格节点为顶点的网格单元，这些网格单元构成一个类似于球体的三维实体，构成该三维实体表面的所有网格面叫做待优化网格节点的“壳”，“壳”与“壳”内部空间一起称为待优化网格节点的局部优化域。

12、接下来定义网格质量函数：

13、

14、其中x是待优化网格节点坐标，fi(x)是局部优化域内第i个网格的网格质量函数。

15、则目标函数a(x)的定义如下式

16、

17、其中ei(x)是fi(x)的倒数，称为错误函数，n是局部优化域内的网格数量，由于最优化理论是求函数的最小值，所以错误函数的存在保证了由其定义的目标函数a(x)完美适应最优化理论的求值逻辑，使得当目标函数取最小值时的坐标x就是该局部优化域对应的待优化网格节点的最优坐标位置。

18、步骤4、基于步骤2中的初始网格，利用步骤3确定的网格单元质量度量准则计算所有网格单元质量。

19、步骤5、设定网格单元质量优化阈值为0.1-0.6之间，数值越大优化效果越好、时空消耗越高，依据步骤4中的网格单元质量计算结果，将低于该优化阈值的网格单元标记为劣质单元，并放入待优化网格单元容器φm。

20、步骤6、创建一个空容器，定义为待优化网格节点容器φv，遍历待优化网格单元容器φm，将每个网格单元上的网格节点放入待优化网格节点容器φv。

21、步骤7、依据网格节点的唯一编号信息，删除待优化网格节点容器φv中的重复网格节点，然后再删除所有边界网格节点，确保容器φv内只有内部网格节点且互不相同，这些内部网格节点称为待优化网格节点。

22、步骤8、逐个取出待优化网格节点容器φv中的待优化网格节点，对该待优化网格节点构造局部优化域，计算局部优化域的目标函数，并采用基于wolfe准则的修正牛顿法，求解局部优化域相应的目标函数最小值，从而对该待优化网格节坐标位置进行优化。

23、具体优化流程如下：

24、在每个待优化网格节点的最优值求解迭代开始前，设置迭代次数k＝0，每次迭代开始前的待优化网格节点坐标为xk，目标函数为a(xk)。

25、迭代开始后，首先计算待优化网格节点当前位置下的目标函数梯度与初始hessian矩阵每次求解与时，对于一阶微分，采用如下的精度为o(h4)的一阶导数中心差分公式

26、

27、其中h为自变量x的微分。

28、对于二阶微分，采用如下的精度为o(h2)的二阶导数中心差分公式

29、

30、对于具体的微分过程，一阶微分和二阶微分均采用具有自适应精度的函数数值微分方法，保证与求解的准确性与稳健性。

31、当不为正定矩阵时，确定一个修正矩阵ek对hessian矩阵进行修改，使得修改后的hessian矩阵bk正定，具体公式如下

32、

33、对于修正矩阵ek的选取，直接取ek＝τki，其中τk是条件数，i是单位矩阵，当τk充分大时，总可以保证修改后的hessian矩阵bk正定。在本发明中为了降低对hessian矩阵的修改导致的计算误差，对τk值的选取有着严格的规定，首先通过hessian矩阵的最小特征值来设定τk的初始值，将其带入公式(6)，若此时的矩阵bk正定，接受该τk值和矩阵bk，继续后续的计算；若此时的矩阵bk不为正定矩阵，则不断将τk加1并将加1后的结果赋值给τk，直到矩阵bk正定，接受此时的τk值和矩阵bk，继续后续的计算。

34、在和bk计算完成后，就可以求解修正的牛顿方程来得到修正的牛顿方向dk了，修正的牛顿方程定义如下

35、

36、牛顿方向dk就是目标函数a(xk)在当前位置下的一个较优下降方向，当待优化网格节点沿着dk方向移动时，目标函数a(xk)的值总是在不断减少。

37、因此，每次迭代后，进行待优化网格节点坐标位置的更新，如下式

38、xk+1＝xk+αkdk   (8)

39、其中αk是步长，不同于经典牛顿法对步长αk的取值恒为1，本流程中对步长αk的确定是基于wolfe准则的，会利用满足wolfe准则的线搜索算法来查找每次迭代的步长αk，相比于其他线搜索算法，该方法求得的步长αk更为接近于最优步长，其中wolfe准则的公式定义如下

40、

41、

42、满足上式的步长αk满足wolfe准则，其中c1,c2∈(0,1)为给定的常数，且c1<c2，步长αk查找完毕后就能通过式(8)来确定本次迭代后待优化网格节点坐标位置xk+1了。

43、迭代的终止条件比较简单，当如下不等式中任意一个不满足时，迭代终止

44、

45、其中αmin为最小步长，设定为1×10-6；kmax为最大迭代次数，设定为7；amin为每次迭代前后的目标函数差值的最小值，设定为1×10-4；在每次迭代后，若以上不等式有任意一个不满足时，说明网格节点当前坐标位置已经非常接近最优坐标位置，再进行迭代所付出的时空消耗与其对局部优化域的优化效果不成正比。

46、步骤9、对经过步骤8优化后的网格整体，更新所有网格单元质量，得到网格质量优化后的高质量网格并输出。

47、与现有的基于优化的光顺的优化算法相比，本发明的有益效果：

48、本发明提出的基于wolfe准则修正牛顿法的三维网格质量优化方法利用了半径比的倒数和作为目标函数，对整个局部优化域进行最优网格节点位置计算，保证了整个局部优化域的计算收敛性和约束完整性，而现有的基于优化的光顺的优化算法则利用数字1与半径比的最大差值作为目标函数，这种方法并不能保证计算结果的有效性，而且因为只计算单个网格的数据，大多数情况下还会导致整个局部优化域的网格质量降低。本发明提出的优化方法能适应三维四面体网格特性，并且由于在三维坐标下hessian矩阵为3×3矩阵，其求解相对比较简单，而现有的优化算法的每次迭代对上一次迭代结果的依赖性较大，计算量也非常大，导致最终优化效果不是很好。因此，本发明提出的优化方法能在保证稳健性的情况下降低计算量的同时达到更高的优化率，最终实现高稳健性、高优化率的网格质量优化。