賈 悅，ANITESCU Cosmin，李 春

基于改進的PHT-樣條自適應等幾何配點法

賈 悅1，ANITESCU Cosmin2，李 春1

(1. 西北工業(yè)大學力學與土木建筑學院，陜西 西安 710129；2.Institute of Structure Mechanics, Bauhaus University Weimar, Thuringia Weimar 99423)

將傳統等幾何配點法擴展至任意高階單元并且滿足自適應局部細分功能，提出一種基于改進的PHT樣條單元的自適應等幾何配點法。改進的PHT樣條單元依然具有傳統PHT樣條單元局部細分功能，但因為傳統PHT樣條函數在層級網格劃分后需要對部分基函數的定義域進行截斷處理，所以在層級細分過于頻繁區(qū)域，部分函數可能因為嚴重變形而影響計算穩(wěn)定性，而改進的PHT樣條函數無需截斷處理，定義域內基函數始終具有穩(wěn)定形態(tài)，這使得改進的PHT樣條單元更適合高階連續(xù)性計算及多層網格細分。該算法結合PHT樣條單元的特點，選取高斯點作為配置點。為了簡化邊界施加條件，采用了耦合線性方程組的方法，在問題域內采用高斯配點法，在問題域邊界采用傳統伽遼金方法，最終耦合2組線性方程組。本算法的局部細分準則基于復原解和復原解誤差。實例計算結果表明，基于改進的PHT樣條的自適應等幾何配點法可以擴展至任意高階單元計算，并滿足最佳收斂率，且與理論值吻合。

等幾何分析方法；配點法；PHT樣條函數；高斯配置點；高階單元；局部細分

等幾何分析(isogeometric analysis，IGA)方法[1]是一種高階連續(xù)有限元計算方法，該方法應用一種基函數同時進行幾何建模和結構分析，從而實現從幾何建模到結構分析的無縫連接。傳統IGA方法是基于非均勻有理B樣條函數(non-uniform rational B-spline，NURBS)[2]，其主要用于計算機輔助設計(computer aided design，CAD)，具有較強的幾何建模作用，因此IGA方法可以用最少的單元精確建模。較傳統有限元法(finite element method，FEM)，因無需幾何重新建模，避免了引入誤差，節(jié)省重新建模時間。另外，因為傳統有限元方法是基于拉格朗日多項式函數，因此該方法在高階連續(xù)性計算中存在穩(wěn)定性問題[3]，所以主流商業(yè)軟件也僅僅基于一階或二階基函數單元計算，但是工程或自然科學中有一些問題所涉及的變量本身是由高階導數定義的，這使得有限元方法很難近似這些變量。較傳統有限元方法，IGA法本質上是一種高階連續(xù)有限元計算方法，因為IGA法應用NURBS樣條函數作為其基函數，該函數具有非常好的高階連續(xù)特性，并且隨著單元階數的升高，計算結果更加穩(wěn)定，所以IGA法可以提供任意高階連續(xù)性計算的需要[4]，其中最顯著的應用是殼單元結構應用[5]。IGA法比較直接的均勻網格細分，其中主要類型包括h-細分、p-細分和k-細分[6-7]。為了改進傳統IGA法的局部細分功能，該方法與一些具有局部細分功能的樣條函數相結合，例如T-樣條[8-9]和PHT樣條[10]。對于高階連續(xù)性計算，IGA方法可以保證任意高階連續(xù)性計算的需要，并且隨著單元階數的升高，該方法計算依然穩(wěn)定[11-12]。目前，IGA法已成功應用于固體力學[6]、流體力學[13]、斷裂力學[14]、生物工程[15]領域中。

IGA方法主要分為：等幾何伽遼金法[1]和等幾何配點法[16-17]2類，并用于計算偏微分方程邊界值問題。其中，傳統等幾何伽遼金方法通過高斯散度定理和積分可加性將偏微分方程的強形式弱化，并進行離散化處理，從而得到線性方程組，最終求解未知量。傳統等幾何伽遼金法具有廣泛的適用性和計算精度。較傳統等幾何伽遼金法，等幾何配點法直接在問題域內定義配置點，然后將配置點帶入偏微分方程強形式，離散化后即可得到配置線性方程組。因此該方法更加容易實現，并且可以降低計算成本和節(jié)省存儲空間，其具有重要的工程應用價值。但是等幾何配點法依然有2大關鍵性問題亟待解決，即高階單元內配置點的選取和邊界條件施加穩(wěn)定性問題。圖1是等幾何配點法的工作流程圖。

圖1 等幾何配點法流程圖

REALI和HUGHES[18]在2015年指出，基于Greville abscissae配點法在奇數階單元計算中收斂率低于理論值。2021年WANG等[19]基于Greville abscissae配點法實現奇數階單元最佳收斂率的計算，該工作采用ANITESCU等[20]提出的超收斂點理論，將伽遼金計算結果二階導數誤差值最小的點作為超收斂點。當然，超收斂點只需在參數單元中計算一次，然后通過單元映射即可得到均勻網格內的超收斂點。但該方法僅限于勻勻網格劃分，對于更一般的情況，例如非均勻網格單元中的超收斂點依然是一個開放性的難題[18]。針對為等幾何配點法邊界條件施加穩(wěn)定性問題，可以尋找預調節(jié)算子[21]，DE LORENZIS等[22]通過加強Neumann邊界條件方法改進等幾何配點法施加的邊界條件。ZHU等[23]提出了基于B++樣條的配點法。另外，混合伽遼金-配點法[24-25]的提出為改善等幾何配點法的穩(wěn)定性提供了新思路。

本文提出一種基于改進的PHT樣條函數的自適應等幾何配點方法。改進的PHT樣條函數[24]繼承傳統PHT樣條函數的局部細分功能，依然通過十字插入實現網格細分。但是，較傳統PHT樣條函數，改進后的PHT樣條基函數始終具有穩(wěn)定形態(tài)。另外，為了保證基函數與配置點個數的相同，采用高斯點作為配置點。在處理邊界條件問題時，為了降低傳統配點法施加邊界條件復雜度且提高計算穩(wěn)定性，本研究采用一種耦合的算法。在問題域內應用等幾何配點法定義線性方程組，在邊界處應用等幾何伽遼金法定義線性方程組，最后將2組耦合為一個整體，并求解未知量。本算法的局部細分準則基于復原解誤差，取誤差較小的點作為差值點，本文選取的差值點為超收斂點。應用高階多項式函數在超收斂點處插值擬合，以得到較近似解及更為精確的復原解，用復原解和近似解之間的差值作為復原解誤差，然后定義一個閾值，將誤差大于該閾值的單元進行細分，最終實現局部細分功能。

1 自適應等幾何配點法

1.1 等幾何配點法

等幾何配點法需要定義2個空間，分別是參數空間和物理空間(圖2)。2個空間可通過關系式

進行相互轉換。其中，為基函數，在本文中即為PHT樣條函數；為基函數對應的控制點，用于參數空間到物理空間的映射變化。參數空間的形狀均為矩形(2D)，物理空間即為任意實際形狀。

與等幾何伽遼金法相比，等幾何配點法直接計算偏微分方程強形式，邊界值問題強形式為

傳統PHT樣條函數繼承了B樣條函數的一些優(yōu)點，如：非負性、局部定義、歸一性。但文獻[26]發(fā)現，在一些過度網格細分時，傳統PHT樣條函數會出現衰變現象。如一些函數的定義域會迅速趨于零，或出現嚴重變形(圖4，函數的最大值在逐漸變小)，該現象被稱為基函數的衰變現象，在實際應用中會影響計算穩(wěn)定性。這種衰變現象是由截斷定義域操作引起的，所以為了防止衰變現象，最直接的辦法是避免截斷定義域操作。因此，文獻[26]基于局部張量積定義一種新的PHT樣條函數，其無需截斷函數定義域，在任何細分情況避免出現基函數衰變現象。定義新的PHT樣條函數，首先要找到與基函數節(jié)點相關聯的定義域網格，然后在該定義域網格上定義4個相關的張量積B樣條方程作為基函數，其始終被定義在矩形網格定義域之內(如圖5(b))，所以無需截斷其定義域。

在定義配置線性方程組對應的系數矩陣時，配置點的個數決定線性方程組系數矩陣的行數，而基函數的個數決定其列數。因此，為了保證生成的系數矩陣是有效方陣，需要保證配置點的個數等于基函數的個數。本文采用高斯點作為配置點，隨著單元階數的增加，每個單元內配置點的個數也隨之增加(圖6和圖7)，但對于非均勻網格劃分，會引入T節(jié)點，并在T節(jié)點周圍增加單元個數，但不增加基函數。若保持T節(jié)點周圍單元內配置點的個數不變，配置點的個數就需大于基函數的個數，最終生成的配置線性方程組的系數矩陣是超定矩陣，于是增加了求解的難度。因此需要對T節(jié)點周圍的配置點做刪減和重置，如圖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連續(xù)，片間單元要求0連續(xù)，所以共有120個基函數(圖8中黑色數字表示域內基函數，紅色數字表示邊界基函數)，需要注意的是片間兩側的基函數也視為邊界基函數。因為該問題應用四階單元，所以每個單元內有6個高斯配置點(圖8中實心點)，整個問題域內共有72個高斯配置點，加上邊界上的48個基函數(圖8中紅色數字)，可以定義120個線性方程。因此，這樣耦合等幾何伽遼金法和等幾何配點法，保證了生成的線性方程組對應的系數矩陣是有效方陣。

圖8 基函數分布和配置點(四階單元)

1.3 基于復原解的局部細分準則

自適應等幾何配點法的局部細分準則基于復原解和復原解誤差，在本研究工作中，3種誤差類型分別是實際誤差、復原解誤差及估計誤差。實際誤差是解析解與近似解之間的差值；復原解誤差是解析解與復原解之間的差值；估計誤差是近似解與復原解之間的差值。其中，因為估計誤差避免了解析解，所以其具有重要的應用價值。本文計算復原解的過程主要分為2步：①在數值解的基礎上選取插值點，這些點可以是超收斂點；②在一些數學假設的前提下，應用高階多項式函數在超收斂點處做插值計算，從而得到更為精確的近似解，即復原解。在本工作中，對于超收斂點的坐標，只需要在參考單元中計算一次，然后將這些點通過坐標變換映射到任意均勻單元，該計算過程可參考文獻[20]。因為復原解是近似解的一個有效上界，得到復原解后，將其作為解析解，從而定義估計誤差。給定一個閾值，標記出誤差值大于該閾值的單元并將其細分。

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

首先，通過經典的平板孔洞模型測試基于改進后的PHT樣條函數的等幾何配置法。極坐標下的精確解為

圖9比較了自適應網格細分和均勻網格細分，圖中紅色和藍色網格分別表示連接的2片結構。圖10計算了應力及該方向上應力誤差。圖11比較了等幾何伽遼金法和高斯等幾何配點法的計算結果，隨著階數的升高，收斂率也隨之增加，并且各階單元均滿足最優(yōu)收斂率，與理論值吻合。

圖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節(jié)點周圍的基函數做截斷定義域處理。但是對于一些網格過度細分的情況，一些基函數可能會出現嚴重變形而影響計算穩(wěn)定性，所以改進后的PHT樣條函數避免截斷定義域，在任何細分情況始終保持鐘形穩(wěn)定形狀。本文將改進的PHT樣條函數作為自適應等幾何配點法的基函數。為了保證基函數的個數與配置點的個數相等，該算法采用高斯點作為配置點，在局部細分過程中，由于T節(jié)點不會增加基函數，所以需要對T節(jié)點周圍單元內的配置點做刪減和重置處理，從而保證在局部細分過程中依然滿足配置點的個數等于基函數的個數，最終得到的線性方程是有效方陣。另外，該算法采用一種耦合的方法施加邊界條件，在問題域內應用等幾何配點定義線性方程組，在問題域邊界應用等幾何伽遼金法定義線性方程組，最后耦合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