基于改進(jìn)的PHT-樣條自適應(yīng)等幾何配點(diǎn)法

賈 悅，ANITESCU Cosmin，李 春

基于改進(jìn)的PHT-樣條自適應(yīng)等幾何配點(diǎn)法

賈 悅1，ANITESCU Cosmin2，李 春1

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

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

等幾何分析方法；配點(diǎn)法；PHT樣條函數(shù)；高斯配置點(diǎn)；高階單元；局部細(xì)分

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

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

圖1 等幾何配點(diǎn)法流程圖

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

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

1 自適應(yīng)等幾何配點(diǎn)法

1.1 等幾何配點(diǎn)法

等幾何配點(diǎn)法需要定義2個(gè)空間，分別是參數(shù)空間和物理空間(圖2)。2個(gè)空間可通過關(guān)系式

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

與等幾何伽遼金法相比，等幾何配點(diǎn)法直接計(jì)算偏微分方程強(qiáng)形式，邊界值問題強(qiáng)形式為

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

在定義配置線性方程組對應(yīng)的系數(shù)矩陣時(shí)，配置點(diǎn)的個(gè)數(shù)決定線性方程組系數(shù)矩陣的行數(shù)，而基函數(shù)的個(gè)數(shù)決定其列數(shù)。因此，為了保證生成的系數(shù)矩陣是有效方陣，需要保證配置點(diǎn)的個(gè)數(shù)等于基函數(shù)的個(gè)數(shù)。本文采用高斯點(diǎn)作為配置點(diǎn)，隨著單元階數(shù)的增加，每個(gè)單元內(nèi)配置點(diǎn)的個(gè)數(shù)也隨之增加(圖6和圖7)，但對于非均勻網(wǎng)格劃分，會(huì)引入T節(jié)點(diǎn)，并在T節(jié)點(diǎn)周圍增加單元個(gè)數(shù)，但不增加基函數(shù)。若保持T節(jié)點(diǎn)周圍單元內(nèi)配置點(diǎn)的個(gè)數(shù)不變，配置點(diǎn)的個(gè)數(shù)就需大于基函數(shù)的個(gè)數(shù)，最終生成的配置線性方程組的系數(shù)矩陣是超定矩陣，于是增加了求解的難度。因此需要對T節(jié)點(diǎn)周圍的配置點(diǎn)做刪減和重置，如圖6和圖7，從而保證基函數(shù)的個(gè)數(shù)等于高斯配置點(diǎn)的個(gè)數(shù)。

圖3 高階單元內(nèi)基函數(shù)分類

圖4 傳統(tǒng)PHT樣條函數(shù)(三階單元)((a)第1層網(wǎng)格細(xì)分；(b)截?cái)喽x域后的基函數(shù)；(c)第2層網(wǎng)格細(xì)分；(d)截?cái)喽x域后的基函數(shù)；(e)第3層網(wǎng)格細(xì)分；(f)截?cái)喽x域后的基函數(shù)；(g～i)基函數(shù)正面圖)

圖5 傳統(tǒng)PHT樣條函數(shù)作用域與改進(jìn)后PHT樣條函數(shù)作用域(三階單元) ((a)傳統(tǒng)PHT樣條函數(shù)作用域；(b)改進(jìn)后PHT樣條函數(shù)作用域；(c)傳統(tǒng)PHT樣條函數(shù)；(d)改進(jìn)后PHT樣條函數(shù))

圖6 層級劃分網(wǎng)格(三階單元)((a)第1層；(b)第2層；(c)第3層；(d)第4層)

圖7 層級劃分網(wǎng)格(四階單元) ((a)第1層；(b)第2層；(c)第3層；(d)第4層)

1.2 耦合線性方程組

等幾何伽遼金法會(huì)自動(dòng)生成對稱系數(shù)矩陣，而等幾何配點(diǎn)法對應(yīng)的系數(shù)矩陣是非對稱的，因?yàn)橄禂?shù)矩陣的行數(shù)是由配置點(diǎn)的個(gè)數(shù)決定的，系數(shù)矩陣的列數(shù)是由基函數(shù)空間中基函數(shù)的個(gè)數(shù)決定的。因此，只有當(dāng)配置點(diǎn)的個(gè)數(shù)等于基函數(shù)的個(gè)數(shù)時(shí)，才能保證系數(shù)線性方程組對應(yīng)的系數(shù)矩陣是有效方陣。如圖8所示，該問題域由2片結(jié)構(gòu)組成，每片包括4個(gè)單元，每個(gè)單元應(yīng)用四階樣條基函數(shù)，由于片內(nèi)要求1連續(xù)，片間單元要求0連續(xù)，所以共有120個(gè)基函數(shù)(圖8中黑色數(shù)字表示域內(nèi)基函數(shù)，紅色數(shù)字表示邊界基函數(shù))，需要注意的是片間兩側(cè)的基函數(shù)也視為邊界基函數(shù)。因?yàn)樵搯栴}應(yīng)用四階單元，所以每個(gè)單元內(nèi)有6個(gè)高斯配置點(diǎn)(圖8中實(shí)心點(diǎn))，整個(gè)問題域內(nèi)共有72個(gè)高斯配置點(diǎn)，加上邊界上的48個(gè)基函數(shù)(圖8中紅色數(shù)字)，可以定義120個(gè)線性方程。因此，這樣耦合等幾何伽遼金法和等幾何配點(diǎn)法，保證了生成的線性方程組對應(yīng)的系數(shù)矩陣是有效方陣。

圖8 基函數(shù)分布和配置點(diǎn)(四階單元)

1.3 基于復(fù)原解的局部細(xì)分準(zhǔn)則

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

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

首先，通過經(jīng)典的平板孔洞模型測試基于改進(jìn)后的PHT樣條函數(shù)的等幾何配置法。極坐標(biāo)下的精確解為

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

圖9 比較自適應(yīng)網(wǎng)格劃分與均勻網(wǎng)格劃分((a，c，e，g)自適應(yīng)網(wǎng)格細(xì)分；(b，d，f，h)均勻網(wǎng)格細(xì)分)

圖10 網(wǎng)格細(xì)分、應(yīng)力圖及其應(yīng)力誤差((a～c)局部網(wǎng)格細(xì)分；(d～f)應(yīng)力圖；(g～i)應(yīng)力誤差)

圖11 高斯配點(diǎn)法誤差分析，三階、四階、五階及六階單元收斂率((a～d)等幾何伽遼金法收斂率；(e～h)高斯配點(diǎn)法收斂率)

3 結(jié)束語

本文提出了一種基于改進(jìn)的PHT樣條函數(shù)的自適應(yīng)等幾何配點(diǎn)法，在定義傳統(tǒng)PHT樣條函數(shù)時(shí)，為了滿足基函數(shù)歸一性，對T節(jié)點(diǎn)周圍的基函數(shù)做截?cái)喽x域處理。但是對于一些網(wǎng)格過度細(xì)分的情況，一些基函數(shù)可能會(huì)出現(xiàn)嚴(yán)重變形而影響計(jì)算穩(wěn)定性，所以改進(jìn)后的PHT樣條函數(shù)避免截?cái)喽x域，在任何細(xì)分情況始終保持鐘形穩(wěn)定形狀。本文將改進(jìn)的PHT樣條函數(shù)作為自適應(yīng)等幾何配點(diǎn)法的基函數(shù)。為了保證基函數(shù)的個(gè)數(shù)與配置點(diǎn)的個(gè)數(shù)相等，該算法采用高斯點(diǎn)作為配置點(diǎn)，在局部細(xì)分過程中，由于T節(jié)點(diǎn)不會(huì)增加基函數(shù)，所以需要對T節(jié)點(diǎn)周圍單元內(nèi)的配置點(diǎn)做刪減和重置處理，從而保證在局部細(xì)分過程中依然滿足配置點(diǎn)的個(gè)數(shù)等于基函數(shù)的個(gè)數(shù)，最終得到的線性方程是有效方陣。另外，該算法采用一種耦合的方法施加邊界條件，在問題域內(nèi)應(yīng)用等幾何配點(diǎn)定義線性方程組，在問題域邊界應(yīng)用等幾何伽遼金法定義線性方程組，最后耦合2組線性方程組為一個(gè)整體，求解未知量。而局部細(xì)分準(zhǔn)則基于復(fù)原解和復(fù)原解誤差，復(fù)原解作為近似解的一個(gè)有效上界，可以當(dāng)作解析解，用于計(jì)算復(fù)原解誤差，然后給定閾值，選定誤差值大于該閾值的單元進(jìn)行網(wǎng)格細(xì)分。該算法已擴(kuò)展至四階、五階和六階高階單元計(jì)算，其計(jì)算結(jié)果可滿足最佳收斂率計(jì)算，與理論值吻合。

[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型細(xì)化方法[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 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ìn), 等. 面向等幾何分析的幾何計(jì)算[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 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

國家自然科學(xué)基金項(xiàng)目(11902263)；陜西省自然科學(xué)基金項(xiàng)目(2019JQ-623)

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

賈 悅(1986–)，男，講師，博士。主要研究方向?yàn)榈葞缀畏治龇椒ā⒌葞缀闻潼c(diǎn)法、圖像配準(zhǔn)技術(shù)。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–)，男，教授，博士。主要研究方向?yàn)槲⒓{米物理力學(xué)。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