基于改进的PHT-样条自适应等几何配点法

贾 悦，ANITESCU Cosmin，李 春

基于改进的PHT-样条自适应等几何配点法

贾 悦1，ANITESCU Cosmin2，李 春1

(1. 西北工业大学力学与土木建筑学院，陕西 西安 710129；2.Institute of Structure Mechanics, Bauhaus University Weimar, Thuringia Weimar 99423)

将传统等几何配点法扩展至任意高阶单元并且满足自适应局部细分功能，提出一种基于改进的PHT样条单元的自适应等几何配点法。改进的PHT样条单元依然具有传统PHT样条单元局部细分功能，但因为传统PHT样条函数在层级网格划分后需要对部分基函数的定义域进行截断处理，所以在层级细分过于频繁区域，部分函数可能因为严重变形而影响计算稳定性，而改进的PHT样条函数无需截断处理，定义域内基函数始终具有稳定形态，这使得改进的PHT样条单元更适合高阶连续性计算及多层网格细分。该算法结合PHT样条单元的特点，选取高斯点作为配置点。为了简化边界施加条件，采用了耦合线性方程组的方法，在问题域内采用高斯配点法，在问题域边界采用传统伽辽金方法，最终耦合2组线性方程组。本算法的局部细分准则基于复原解和复原解误差。实例计算结果表明，基于改进的PHT样条的自适应等几何配点法可以扩展至任意高阶单元计算，并满足最佳收敛率，且与理论值吻合。

等几何分析方法；配点法；PHT样条函数；高斯配置点；高阶单元；局部细分

等几何分析(isogeometric analysis，IGA)方法[1]是一种高阶连续有限元计算方法，该方法应用一种基函数同时进行几何建模和结构分析，从而实现从几何建模到结构分析的无缝连接。传统IGA方法是基于非均匀有理B样条函数(non-uniform rational B-spline，NURBS)[2]，其主要用于计算机辅助设计(computer aided design，CAD)，具有较强的几何建模作用，因此IGA方法可以用最少的单元精确建模。较传统有限元法(finite element method，FEM)，因无需几何重新建模，避免了引入误差，节省重新建模时间。另外，因为传统有限元方法是基于拉格朗日多项式函数，因此该方法在高阶连续性计算中存在稳定性问题[3]，所以主流商业软件也仅仅基于一阶或二阶基函数单元计算，但是工程或自然科学中有一些问题所涉及的变量本身是由高阶导数定义的，这使得有限元方法很难近似这些变量。较传统有限元方法，IGA法本质上是一种高阶连续有限元计算方法，因为IGA法应用NURBS样条函数作为其基函数，该函数具有非常好的高阶连续特性，并且随着单元阶数的升高，计算结果更加稳定，所以IGA法可以提供任意高阶连续性计算的需要[4]，其中最显著的应用是壳单元结构应用[5]。IGA法比较直接的均匀网格细分，其中主要类型包括h-细分、p-细分和k-细分[6-7]。为了改进传统IGA法的局部细分功能，该方法与一些具有局部细分功能的样条函数相结合，例如T-样条[8-9]和PHT样条[10]。对于高阶连续性计算，IGA方法可以保证任意高阶连续性计算的需要，并且随着单元阶数的升高，该方法计算依然稳定[11-12]。目前，IGA法已成功应用于固体力学[6]、流体力学[13]、断裂力学[14]、生物工程[15]领域中。

IGA方法主要分为：等几何伽辽金法[1]和等几何配点法[16-17]2类，并用于计算偏微分方程边界值问题。其中，传统等几何伽辽金方法通过高斯散度定理和积分可加性将偏微分方程的强形式弱化，并进行离散化处理，从而得到线性方程组，最终求解未知量。传统等几何伽辽金法具有广泛的适用性和计算精度。较传统等几何伽辽金法，等几何配点法直接在问题域内定义配置点，然后将配置点带入偏微分方程强形式，离散化后即可得到配置线性方程组。因此该方法更加容易实现，并且可以降低计算成本和节省存储空间，其具有重要的工程应用价值。但是等几何配点法依然有2大关键性问题亟待解决，即高阶单元内配置点的选取和边界条件施加稳定性问题。图1是等几何配点法的工作流程图。

图1 等几何配点法流程图

REALI和HUGHES[18]在2015年指出，基于Greville abscissae配点法在奇数阶单元计算中收敛率低于理论值。2021年WANG等[19]基于Greville abscissae配点法实现奇数阶单元最佳收敛率的计算，该工作采用ANITESCU等[20]提出的超收敛点理论，将伽辽金计算结果二阶导数误差值最小的点作为超收敛点。当然，超收敛点只需在参数单元中计算一次，然后通过单元映射即可得到均匀网格内的超收敛点。但该方法仅限于匀匀网格划分，对于更一般的情况，例如非均匀网格单元中的超收敛点依然是一个开放性的难题[18]。针对为等几何配点法边界条件施加稳定性问题，可以寻找预调节算子[21]，DE LORENZIS等[22]通过加强Neumann边界条件方法改进等几何配点法施加的边界条件。ZHU等[23]提出了基于B++样条的配点法。另外，混合伽辽金-配点法[24-25]的提出为改善等几何配点法的稳定性提供了新思路。

本文提出一种基于改进的PHT样条函数的自适应等几何配点方法。改进的PHT样条函数[24]继承传统PHT样条函数的局部细分功能，依然通过十字插入实现网格细分。但是，较传统PHT样条函数，改进后的PHT样条基函数始终具有稳定形态。另外，为了保证基函数与配置点个数的相同，采用高斯点作为配置点。在处理边界条件问题时，为了降低传统配点法施加边界条件复杂度且提高计算稳定性，本研究采用一种耦合的算法。在问题域内应用等几何配点法定义线性方程组，在边界处应用等几何伽辽金法定义线性方程组，最后将2组耦合为一个整体，并求解未知量。本算法的局部细分准则基于复原解误差，取误差较小的点作为差值点，本文选取的差值点为超收敛点。应用高阶多项式函数在超收敛点处插值拟合，以得到较近似解及更为精确的复原解，用复原解和近似解之间的差值作为复原解误差，然后定义一个阈值，将误差大于该阈值的单元进行细分，最终实现局部细分功能。

1 自适应等几何配点法

1.1 等几何配点法

等几何配点法需要定义2个空间，分别是参数空间和物理空间(图2)。2个空间可通过关系式

进行相互转换。其中，为基函数，在本文中即为PHT样条函数；为基函数对应的控制点，用于参数空间到物理空间的映射变化。参数空间的形状均为矩形(2D)，物理空间即为任意实际形状。

与等几何伽辽金法相比，等几何配点法直接计算偏微分方程强形式，边界值问题强形式为

传统PHT样条函数继承了B样条函数的一些优点，如：非负性、局部定义、归一性。但文献[26]发现，在一些过度网格细分时，传统PHT样条函数会出现衰变现象。如一些函数的定义域会迅速趋于零，或出现严重变形(图4，函数的最大值在逐渐变小)，该现象被称为基函数的衰变现象，在实际应用中会影响计算稳定性。这种衰变现象是由截断定义域操作引起的，所以为了防止衰变现象，最直接的办法是避免截断定义域操作。因此，文献[26]基于局部张量积定义一种新的PHT样条函数，其无需截断函数定义域，在任何细分情况避免出现基函数衰变现象。定义新的PHT样条函数，首先要找到与基函数节点相关联的定义域网格，然后在该定义域网格上定义4个相关的张量积B样条方程作为基函数，其始终被定义在矩形网格定义域之内(如图5(b))，所以无需截断其定义域。

在定义配置线性方程组对应的系数矩阵时，配置点的个数决定线性方程组系数矩阵的行数，而基函数的个数决定其列数。因此，为了保证生成的系数矩阵是有效方阵，需要保证配置点的个数等于基函数的个数。本文采用高斯点作为配置点，随着单元阶数的增加，每个单元内配置点的个数也随之增加(图6和图7)，但对于非均匀网格划分，会引入T节点，并在T节点周围增加单元个数，但不增加基函数。若保持T节点周围单元内配置点的个数不变，配置点的个数就需大于基函数的个数，最终生成的配置线性方程组的系数矩阵是超定矩阵，于是增加了求解的难度。因此需要对T节点周围的配置点做删减和重置，如图6和图7，从而保证基函数的个数等于高斯配置点的个数。

图3 高阶单元内基函数分类

图4 传统PHT样条函数(三阶单元)((a)第1层网格细分；(b)截断定义域后的基函数；(c)第2层网格细分；(d)截断定义域后的基函数；(e)第3层网格细分；(f)截断定义域后的基函数；(g～i)基函数正面图)

图5 传统PHT样条函数作用域与改进后PHT样条函数作用域(三阶单元) ((a)传统PHT样条函数作用域；(b)改进后PHT样条函数作用域；(c)传统PHT样条函数；(d)改进后PHT样条函数)

图6 层级划分网格(三阶单元)((a)第1层；(b)第2层；(c)第3层；(d)第4层)

图7 层级划分网格(四阶单元) ((a)第1层；(b)第2层；(c)第3层；(d)第4层)

1.2 耦合线性方程组

等几何伽辽金法会自动生成对称系数矩阵，而等几何配点法对应的系数矩阵是非对称的，因为系数矩阵的行数是由配置点的个数决定的，系数矩阵的列数是由基函数空间中基函数的个数决定的。因此，只有当配置点的个数等于基函数的个数时，才能保证系数线性方程组对应的系数矩阵是有效方阵。如图8所示，该问题域由2片结构组成，每片包括4个单元，每个单元应用四阶样条基函数，由于片内要求1连续，片间单元要求0连续，所以共有120个基函数(图8中黑色数字表示域内基函数，红色数字表示边界基函数)，需要注意的是片间两侧的基函数也视为边界基函数。因为该问题应用四阶单元，所以每个单元内有6个高斯配置点(图8中实心点)，整个问题域内共有72个高斯配置点，加上边界上的48个基函数(图8中红色数字)，可以定义120个线性方程。因此，这样耦合等几何伽辽金法和等几何配点法，保证了生成的线性方程组对应的系数矩阵是有效方阵。

图8 基函数分布和配置点(四阶单元)

1.3 基于复原解的局部细分准则

自适应等几何配点法的局部细分准则基于复原解和复原解误差，在本研究工作中，3种误差类型分别是实际误差、复原解误差及估计误差。实际误差是解析解与近似解之间的差值；复原解误差是解析解与复原解之间的差值；估计误差是近似解与复原解之间的差值。其中，因为估计误差避免了解析解，所以其具有重要的应用价值。本文计算复原解的过程主要分为2步：①在数值解的基础上选取插值点，这些点可以是超收敛点；②在一些数学假设的前提下，应用高阶多项式函数在超收敛点处做插值计算，从而得到更为精确的近似解，即复原解。在本工作中，对于超收敛点的坐标，只需要在参考单元中计算一次，然后将这些点通过坐标变换映射到任意均匀单元，该计算过程可参考文献[20]。因为复原解是近似解的一个有效上界，得到复原解后，将其作为解析解，从而定义估计误差。给定一个阈值，标记出误差值大于该阈值的单元并将其细分。

2 举例——2D平板孔洞模型

首先，通过经典的平板孔洞模型测试基于改进后的PHT样条函数的等几何配置法。极坐标下的精确解为

图9比较了自适应网格细分和均匀网格细分，图中红色和蓝色网格分别表示连接的2片结构。图10计算了应力及该方向上应力误差。图11比较了等几何伽辽金法和高斯等几何配点法的计算结果，随着阶数的升高，收敛率也随之增加，并且各阶单元均满足最优收敛率，与理论值吻合。

图9 比较自适应网格划分与均匀网格划分((a，c，e，g)自适应网格细分；(b，d，f，h)均匀网格细分)

图10 网格细分、应力图及其应力误差((a～c)局部网格细分；(d～f)应力图；(g～i)应力误差)

图11 高斯配点法误差分析，三阶、四阶、五阶及六阶单元收敛率((a～d)等几何伽辽金法收敛率；(e～h)高斯配点法收敛率)

3 结束语

本文提出了一种基于改进的PHT样条函数的自适应等几何配点法，在定义传统PHT样条函数时，为了满足基函数归一性，对T节点周围的基函数做截断定义域处理。但是对于一些网格过度细分的情况，一些基函数可能会出现严重变形而影响计算稳定性，所以改进后的PHT样条函数避免截断定义域，在任何细分情况始终保持钟形稳定形状。本文将改进的PHT样条函数作为自适应等几何配点法的基函数。为了保证基函数的个数与配置点的个数相等，该算法采用高斯点作为配置点，在局部细分过程中，由于T节点不会增加基函数，所以需要对T节点周围单元内的配置点做删减和重置处理，从而保证在局部细分过程中依然满足配置点的个数等于基函数的个数，最终得到的线性方程是有效方阵。另外，该算法采用一种耦合的方法施加边界条件，在问题域内应用等几何配点定义线性方程组，在问题域边界应用等几何伽辽金法定义线性方程组，最后耦合2组线性方程组为一个整体，求解未知量。而局部细分准则基于复原解和复原解误差，复原解作为近似解的一个有效上界，可以当作解析解，用于计算复原解误差，然后给定阈值，选定误差值大于该阈值的单元进行网格细分。该算法已扩展至四阶、五阶和六阶高阶单元计算，其计算结果可满足最佳收敛率计算，与理论值吻合。

[1] HUGHES T J R, COTTRELL J A, BAZILEVS Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[J]. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39-41): 4135-4195.

[2] PIEGL L, TILLER W. The NURBS book[M]. Heidelberg: Springer, 1997: 117-122.

[3] COTTRELL J A, HUGHES T J R, BAZILEVS Y. Isogeometric analysis[M]. Chichester: John Wiley & Sons, Ltd, 2009: 31-33.

[4] LIPTON S, EVANS J A, BAZILEVS Y, et al. Robustness of isogeometric structural discretizations under severe mesh distortion[J]. Computer Methods in Applied Mechanics and Engineering, 2010, 199(5-8): 357-373.

[5] CASQUERO H, WEI X D, TOSHNIWAL D, et al. Seamless integration of design and Kirchhoff-Love shell analysis using analysis-suitable unstructured T-splines[J]. Computer Methods in Applied Mechanics and Engineering, 2020, 360: 112765.

[6] COTTRELL J A, REALI A, BAZILEVS Y, et al. Isogeometric analysis of structural vibrations[J]. Computer Methods in Applied Mechanics and Engineering, 2006, 195(41-43): 5257-5296.

[7] 徐岗, 王毅刚, 胡维华. 等几何分析中的r-p型细化方法[J]. 计算机辅助设计与图形学学报, 2011, 23(12): 2019-2024.

XU G, WANG Y G, HU W H. R-p-refinement in isogeometric analysis[J]. Journal of Computer-Aided Design & Computer Graphics, 2011, 23(12): 2019-2024 (in Chinese).

[8] SCOTT M A, LI X, SEDERBERG T W, et al. Local refinement of analysis-suitable T-splines[J]. Computer Methods in Applied Mechanics and Engineering, 2012, 213-216: 206-222.

[9] LI X, SCOTT M A. Analysis-suitable T-splines: Characterization, refineability, and approximation[J]. Mathematical Models and Methods in Applied Sciences, 2014, 24(6): 1141-1164.

[10] DENG J S, CHEN F L, LI X, et al. Polynomial splines over hierarchical T-meshes[J]. Graphical Models, 2008, 70(4): 76-86.

[11] VERHOOSEL C V, SCOTT M A, HUGHES T J R, et al. An isogeometric analysis approach to gradient damage models[J]. International Journal for Numerical Methods in Engineering, 2011, 86(1): 115-134.

[12] 徐岗, 李新, 黄章进, 等. 面向等几何分析的几何计算[J]. 计算机辅助设计与图形学学报, 2015, 27(4): 570-581.

XU G, LI X, HUANG Z J, et al. Geometric computing for isogeometric analysis[J]. Journal of Computer-Aided Design & Computer Graphics, 2015, 27(4): 570-581 (in Chinese).

[13] BAZILEVS Y, CALO V M, HUGHES T J R, et al. Isogeometric fluid-structure interaction: theory, algorithms, and computations[J]. Computational Mechanics, 2008, 43(1): 3-37.

[14] DE LUYCKER E, BENSON D J, BELYTSCHKO T, et al. X-FEM in isogeometric analysis for linear fracture mechanics[J]. International Journal for Numerical Methods in Engineering, 2011, 87(6): 541-565.

[15] KAMENSKY D, HSU M C, SCHILLINGER D, et al. An immersogeometric variational framework for fluid-structure interaction: application to bioprosthetic heart valves[J]. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 1005-1053.

[16] AURICCHIO F, BEIRÃO DA VEIGA L, HUGHES T J R, et al. Isogeometric collocation methods[J]. Mathematical Models and Methods in Applied Sciences, 2010, 20(11): 2075-2107.

[17] AURICCHIO F, BEIRÃO DA VEIGA L, HUGHES T J R, et al. Isogeometric collocation for elastostatics and explicit dynamics[J]. Computer Methods in Applied Mechanics and Engineering, 2012, 249-252: 2-14.

[18] REALI A, HUGHES T J R. An introduction to isogeometric collocation methods[M]//Isogeometric Methods for Numerical Simulation. Vienna: Springer Vienna, 2015: 173-204.

[19] WANG D D, QI D L, LI X W. Superconvergent isogeometric collocation method with Greville points[J]. Computer Methods in Applied Mechanics and Engineering, 2021, 377: 113689.

[20] ANITESCU C, JIA Y, ZHANG Y J, et al. An isogeometric collocation method using superconvergent points[J]. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 1073-1097.

[21] BEIRÃO DA VEIGA L, CHO D, PAVARINO L F, et al. Overlapping Schwarz preconditioners for isogeometric collocation methods[J]. Computer Methods in Applied Mechanics and Engineering, 2014, 278: 239-253.

[22] DE LORENZIS L, EVANS J A, HUGHES T J R, et al. Isogeometric collocation: Neumann boundary conditions and contact[J]. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 21-54.

[23] ZHU X F, HU P, MA Z D. B++ splines with applications to isogeometric analysis[J]. Computer Methods in Applied Mechanics and Engineering, 2016, 311: 503-536.

[24] KLINKEL S, CHEN L, DORNISCH W. A NURBS based hybrid collocation-Galerkin method for the analysis of boundary represented solids[J]. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 689-711.

[25] SCHILLINGER D, BORDEN M J, STOLARSKI H K. Isogeometric collocation for phase-field fracture models[J]. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 583-610.

[26] KANG H M, XU J L, CHEN F L, et al. A new basis for PHT-splines[J]. Graphical Models, 2015, 82: 149-159.

[27] DE BOOR C, SWARTZ B. Collocation at Gaussian points[J]. SIAM Journal on Numerical Analysis, 1973, 10(4): 582-606.

An adaptive isogeometric collocation method with improved PHT-splines

JIA Yue1, ANITESCU Cosmin2, LI Chun1

(1. School of Mechanics, Civil Engineering and Architecture, Northwestern Polytechnical University, Xi’an Shaanxi 710129, China; 2. Institute of Structure Mechanics, Bauhaus University Weimar, Weimar Thuringia 99423, Germany)

The Gaussian isogeometric analysis (IGA) collocation method was extended to arbitrary higher order polynomial degrees. The current IGA collocation method applied a new hierarchical basis over T-meshes (PHT-splines), which took advantage of the tensor product structure to prevent the decay phenomenon from happening in the original PHT basis. The improved method collocated at Gaussian points as the superconvergent points for the new PHT elements. Based on the new PHT basis, the current collocation method can be extended to arbitrary higher order approximation. In order to simplify the collocation boundary condition, a hybrid method was adopted to impose the boundary condition, using the Galerkin method for the boundary part and combing with the collocation solving system. The local refinement strategy was driven by a recovery-based error estimator that invoked computing an improved approximation without knowledge of the exact solution. The proposed collocation method can obtain the optimal convergent rates, compared with the IGA Galerkin method.

isogeometric analysis; collocation method; PHT-splines; Gaussian collocation points; higher order elements; adaptive refinement

15 June，2021；

TP 391

10.11996/JG.j.2095-302X.2022010110

A

2095-302X(2022)01-0110-08

2021-06-15；

2021-07-04

4 July，2021

国家自然科学基金项目(11902263)；陕西省自然科学基金项目(2019JQ-623)

National Natural Science Foundation of China (11902263);Shaanxi Province Science Foundation (2019JQ-623)

贾 悦(1986–)，男，讲师，博士。主要研究方向为等几何分析方法、等几何配点法、图像配准技术。E-mail：yuejia@nwpu.edu.cn

JIA Yue (1986–), lecturer, Ph.D. His main research interests cover isogeometric analysis, isogeometric collocation method, image registration. E-mail：yuejia@nwpu.edu.cn

李 春(1979–)，男，教授，博士。主要研究方向为微纳米物理力学。E-mail：lichun@nwpu.edu.cn

LI Chun (1979–), professor, Ph.D. His main research interest covers micro-nano physical mechanics. E-mail：lichun@nwpu.edu.cn