榮吉利，熊麗園，劉鋮，辛鵬飛,2，劉志超

(1.北京理工大學(xué) 宇航學(xué)院，北京 100081；2.北京空間飛行器總體設(shè)計部，北京 100094；3.空間物理重點實驗室，北京 100076)

為了實現(xiàn)計算機輔助設(shè)計(CAD)與計算機輔助工程(CAE)的融合，Hughes等[1]把非均勻有理B樣條(NURBS)基函數(shù)引入到等參有限元分析中，建立了等幾何分析方法(IGA).與傳統(tǒng)有限元[2-3]拉格朗日多項式插值方法不同，等幾何分析方法將CAD中樣條基函數(shù)作為插值函數(shù)直接用于有限元分析中，能夠精確描述復(fù)雜幾何構(gòu)型，消除了幾何誤差.

在等幾何分析中，B樣條和NURBS是兩類常用的樣條基函數(shù).對于2維或3維問題，多變量B樣條與NURBS通常采用張量積形式構(gòu)造，其網(wǎng)格局部加密將會導(dǎo)致全局網(wǎng)格的變化，無法實現(xiàn)局部細(xì)化，極大地降低了等幾何分析方法的計算效率.而為實現(xiàn)網(wǎng)格局部細(xì)化及自適應(yīng)網(wǎng)格，多種可局部加密的樣條函數(shù)應(yīng)運而生，例如，T樣條、LR樣條(Local Refined-Splines)以及層次B樣條(Hierarchical B-Splines，簡稱HB樣條).在CAD領(lǐng)域，Sederber[4]首次提出了T樣條概念，允許網(wǎng)格T節(jié)點的存在，實現(xiàn)了網(wǎng)格的局部細(xì)化.Bazilevst[5]將T樣條引入等幾何分析，研究表明T樣條基函數(shù)可能存在線性相關(guān)性，將導(dǎo)致系統(tǒng)剛度矩陣嚴(yán)重數(shù)值病態(tài).同時，Dokken等[6]推廣NURBS，提出了LR樣條，它通過在張量積形式的B樣條網(wǎng)格中插入節(jié)點線段(一整條節(jié)點線中的一段)，實現(xiàn)網(wǎng)格局部細(xì)化.但LR樣條基函數(shù)仍可能存在線性相關(guān)性，并且該樣條基函數(shù)構(gòu)造較為復(fù)雜.此外，F(xiàn)orsey與Bartels[7]提出了層次B樣條(HB樣條)，HB樣條基函數(shù)由一系列嵌套的張量積樣條基函數(shù)構(gòu)成，豐富了樣條模型的局部特性.Vuong[8]給出了一種構(gòu)造HB樣條空間的一組線性無關(guān)的基函數(shù)的方法，并將其應(yīng)用在自適應(yīng)等幾何分析中，實現(xiàn)了L型區(qū)域靜態(tài)熱傳導(dǎo)問題的網(wǎng)格局部細(xì)化.Henning等[9]將HB樣條應(yīng)用于裂紋擴展問題的等幾何分析中，在裂紋擴展過程中構(gòu)造了基于HB樣條的自適應(yīng)網(wǎng)格，結(jié)果表明該方法能夠有效地提高計算效率.HB樣條具有嵌套性，它的基函數(shù)是由多個細(xì)分層的基函數(shù)組成，這使得HB樣條在自適應(yīng)等幾何分析中表現(xiàn)出巨大的優(yōu)勢.此外，HB樣條的構(gòu)造簡單，這些優(yōu)點使得HB樣條成為等幾何分析的研究熱點之一.上述研究大多僅涉及結(jié)構(gòu)靜力學(xué)問題，而尚無關(guān)于等幾何分析處理結(jié)構(gòu)動力學(xué)問題中局部細(xì)化算法研究.

在等幾何分析中，為了實現(xiàn)IGA單元形函數(shù)的統(tǒng)一化，Borden等[10]首次提出了NURBS的貝塞爾提取方法(Bézier extraction)，采用描述貝賽爾曲線的伯恩斯坦基函數(shù)表示NURBS基函數(shù)，使得每個NURBS單元的形函數(shù)都相同，由此可以將等幾何分析程序與傳統(tǒng)有限元程序融合.之后，Scott等[11]將貝塞爾提取應(yīng)用于基于T樣條的等幾何分析中.D’Angella等[12]將貝塞爾提取用于HB樣條中，提出了多層貝塞爾提取方法(multi-level Bézier extraction)，并證明該方法適用于非線性問題和并行計算.

為此，在上述研究基礎(chǔ)上，探討了大變形平面結(jié)構(gòu)動力學(xué)問題中的網(wǎng)格自適應(yīng)算法，給出了基于HB樣條的網(wǎng)格局部加密算法，并分析了HB樣條基函數(shù)的構(gòu)造以及多層貝塞爾提取方法.

1 HB樣條和多層貝塞爾提取算法

1.1 貝塞爾曲線、B樣條和NURBS

一條p階貝塞爾曲線是p+1個伯恩斯坦基函數(shù)的線性組合，其表達(dá)式為

(1)

式中：Pi為貝塞爾控制點；ξ∈[-1,1]為貝塞爾單元的參數(shù)坐標(biāo).伯恩斯坦基函數(shù)Bi,p可表示為

(1-ξ)p-(i-1)(1+ξ)(i-1),

(2)

其中，1 ≤i≤p+1.從伯恩斯坦基的表達(dá)式可以看出，給定基函數(shù)的階數(shù)p,p+1個伯恩斯坦基函數(shù)可唯一給出，即每個單元上的基函數(shù)都相同，與有限元方法中的單元形函數(shù)具有類似的性質(zhì).

此外，一個單變量p階B樣條基函數(shù)定義在一個遞增的節(jié)點矢量Ξ上，

Ξ=[ξ1ξ2…ξn+p+1],

(3)

式中：ξ為量綱一參數(shù)坐標(biāo)；n為基函數(shù)個數(shù).節(jié)點矢量組成的區(qū)域稱為參數(shù)域，參數(shù)域含有的節(jié)點矢量的個數(shù)決定了幾何模型的維度.基于Cox-de Boor遞推公式，一個p階B樣條基函數(shù)可表示為兩個p-1階B樣條基函數(shù)的線性組合為

(4)

其中，規(guī)定0/0=0.第i個基函數(shù)Ni,p與ξi,ξi+1,,ξi+p+1共p+2個節(jié)點矢量元素有關(guān)，這些元素組成的區(qū)間[ξi,ξi+p+1]稱為基函數(shù)Ni,p的支持域，記任意基函數(shù)f的支持域為suppf.B樣條基函數(shù)連續(xù)性由節(jié)點矢量元素重復(fù)度決定，對于重復(fù)度為k的節(jié)點，基函數(shù)Cp-k連續(xù)，而在節(jié)點之間C∞連續(xù).一般而言，對于開放式節(jié)點矢量，第1個結(jié)點和最后一個結(jié)點的重復(fù)度為k=p+1.

在B樣條基函數(shù)基礎(chǔ)上，引入權(quán)因子可推廣至NURBS樣條基函數(shù)，NURBS基函數(shù)除具備B樣條良好性質(zhì)外，還能精確描述各類二次曲線與曲面，例如圓、雙曲線、橢圓面等.一個p階NURBS樣條基函數(shù)為

(5)

上述單變量樣條函數(shù)可方便推廣至多維情況，以2維問題為例，雙變量張量積形式的NURBS基函數(shù)為

(6)

式中：wi,j為權(quán)因子；Ni,p(ξ)和Mi,q(η)分別為p階和q階B樣條基函數(shù).由此，NURBS曲面中任意一點的位置坐標(biāo)矢量可表示為

(7)

式中：{Pi,j}(i=1,2,,n，j=1,2,,m)為曲面控制點，Ξ=[ξiξ2…ξn+p+1]和Η=[η1η2…ηm+q+1]分別為ξ和η方向的節(jié)點向量.

在等幾何分析中，對于1維情況，將節(jié)點矢量中兩個相鄰的不等節(jié)點區(qū)間e=[ξi,ξi+1],ξi≠ξi+1視為一個單元.由于參數(shù)域與幾何模型存在映射關(guān)系，這些單元把實際物理模型劃分為多個一一對應(yīng)的子區(qū)域.圖1為一個定義在節(jié)點矢量Ξ1=[0 0 0 1 2 2 3 4 4 4]的單變量B樣條基函數(shù)分布圖，該節(jié)點矢量包含4個單元，單元e2=[1,2]用點劃線表示.圖2表示定義在相同節(jié)點矢量Ξ2×Ξ2上的雙變量基函數(shù)分布圖，Ξ2=[0 0 0 1 2 3 4 4 4].

1.2 貝塞爾提取與節(jié)點插入算法

在傳統(tǒng)有限元方法中，一類單元一般具有相同的形函數(shù).然而，在等幾何分析中，采用的B樣條基函數(shù)定義在全局的節(jié)點矢量上，單元的基函數(shù)通常不同.為了實現(xiàn)等幾何分析程序與現(xiàn)有的有限元程序相融合，貝塞爾提取技術(shù)被廣泛應(yīng)用于各類樣條基函數(shù)，如B樣條，NURBS以及T樣條.貝塞爾提取的核心思想是采用節(jié)點插入算法，實現(xiàn)伯恩斯坦基函數(shù)來表示全局定義下的B樣條基函數(shù).

(8)

其中，

(9)

為了保持節(jié)點插入中幾何構(gòu)型不變，應(yīng)滿足

(CTP)TB(ξ)=PTCB(ξ).

(10)

因此，B樣條基函數(shù)與伯恩斯坦基函數(shù)之間的關(guān)系為

N(ξ)=CB(ξ),

(11)

其中，線性算子C稱為全局貝塞爾提取算子[11]，它能將任意階連續(xù)的B樣條基函數(shù)表示成分段C0連續(xù)的伯恩斯坦基函數(shù)的線性組合.單元貝塞爾提取算子Ce可從全局貝塞爾提取算子C獲得.由此，單元B樣條基函數(shù)可表示為

(12)

(13)

其中，兩個矩陣的克羅內(nèi)克積A?B的表達(dá)式為

(14)

上述方法可直接推廣至NURBS基函數(shù)貝塞爾提取中，其基函數(shù)可表示為

(15)

式中：we為單元的權(quán)系數(shù)矢量；We為單元權(quán)系數(shù)的對角矩陣.

1.3 HB樣條構(gòu)造

為解決張量積B樣條局部細(xì)化問題，本小節(jié)給出HB樣條基函數(shù)構(gòu)造方式.

HB樣條基函數(shù)包含有各個細(xì)分層的B樣條基函數(shù)，它是一組函數(shù)N=∪lNl作為分析和幾何描述的基函數(shù).HB樣條基函數(shù)的構(gòu)造從一個定義在區(qū)域Ω0的NURBS基函數(shù)N0上展開的.HB樣條的基函數(shù)K通過下面的遞推公式給出.

① 初始：K0={τ∈N0:suppτ≠?}.

③K=KN-1.

1.4 多層貝塞爾提取

為了實現(xiàn)HB樣條貝塞爾提取，首先定義相鄰兩層間基函數(shù)關(guān)系.基于節(jié)點插入算法，第l層的B樣條基Nl都可以表示成細(xì)化的第l+1層的B樣條基Nl+1的線性組合為

Nl=Ml,l+1Nl+1,

(16)

式中M為細(xì)分算子[13-14].因此，相鄰連續(xù)層之間的控制點的關(guān)系為

(17)

上述細(xì)化算法可保持幾何構(gòu)型不變.兩個不連續(xù)層I,J(J>I)之間的細(xì)分算子可由相鄰細(xì)分算子連乘得到

(18)

(19)

通過細(xì)分算子建立各個細(xì)分層基函數(shù)之間的關(guān)系.

其次，對于HB樣條基函數(shù)的貝塞爾提取算子，由于HB樣條的基函數(shù)包含有各個不同細(xì)分層的基函數(shù)，每個HB樣條的單元的形函數(shù)包括多層B樣條基函數(shù).將單元ε′ei(第l層的第ei個單元)的形函數(shù)統(tǒng)一用該單元所在層l的B樣條基函數(shù)Nl表示.在第l層進(jìn)行貝塞爾提取，從而得到多層貝塞爾提取.下面通過一個例子，如圖5所示，具體說明多層貝塞爾提取算子構(gòu)造過程.

(20)

式中MLε3稱為多層細(xì)分算子，它是通過組合多個細(xì)分算子的行向量來推導(dǎo)得到.例如，單元的多層細(xì)分算子是通過組合R0,2，R1,2和R2,2中相應(yīng)的行向量得到.通過多層細(xì)分算子，則層次B樣條的基函數(shù)可以表示成單一層標(biāo)準(zhǔn)B樣條的基函數(shù).將第2層的B樣條基函數(shù)通過貝塞爾提取，表示成伯恩斯坦基函數(shù)的組合，

(21)

式中Cε5為單元的貝塞爾提取算子.最后，單元ε5的形函數(shù)就表示為伯恩斯坦基函數(shù)的線性組合，

(22)

式中MLε5Cε5稱為多層貝塞爾提取算子[12].

由此，對于HB樣條基函數(shù)，采用伯恩斯坦基函數(shù)來表示，實現(xiàn)單元形函數(shù)統(tǒng)一化，有效地簡化等幾何分析程序，更加方便地與傳統(tǒng)有限元程序相融合.

2 大變形平面結(jié)構(gòu)動力學(xué)建模方法

2.1 大變形平面等幾何有限單元

基于等幾何分析方法，采用HB樣條基函數(shù)描述彈性體構(gòu)型，并通過連續(xù)介質(zhì)力理論描述彈性體的變形，從而建立平面板結(jié)構(gòu)的動力學(xué)模型.平面單元在發(fā)生大變形后，其變形后的構(gòu)型為當(dāng)前構(gòu)型，初始構(gòu)型為參考構(gòu)型，單元的參數(shù)域單元可視為母單元構(gòu)型，它們之間的關(guān)系如圖6所示.

考慮任意平面板結(jié)構(gòu)，其初始幾何構(gòu)型通過NURBS曲面描述.以NURBS曲面片中的一個單元[ξi,ξi+1]×[ηj,ηj+1]為研究對象.在發(fā)生大轉(zhuǎn)動大變形后其上任意一點P的位置坐標(biāo)可表示為

r(ξ,η,t)=r0(ξ,η)+u(ξ,η,t)=S(ξ,η)e(t),

(23)

式中：u(ξ,η,t)為點P的位移，S(ξ,η)為NURBS基函數(shù)在定義域[ξi,ξi+1]×[ηj,ηj+1]內(nèi)的非零元素，也就是單元的形函數(shù)，e(t)為NURBS單元控制點，稱為單元廣義坐標(biāo).為了提高計算效率，對公式(23)中的單元廣義坐標(biāo)進(jìn)行重排，點P的位置坐標(biāo)可重新表示為[15]

(24)

? ?

(25)

基于連續(xù)介質(zhì)力學(xué)理論，平面等幾何單元的變形梯度張量F可表示為

(26)

(27)

其中，i=1,2,,2×[(p+1)×(q+1)]，m=1,2,a=1,2,n=1,2，單元彈性力可表示為應(yīng)變能對廣義坐標(biāo)的1次偏導(dǎo)數(shù)，

(28)

式中D為4階彈性張量.將式(27)帶入式(28)，則彈性力的分量可進(jìn)一步表示為

(29)

2.2 大變形結(jié)構(gòu)動力學(xué)方程時間積分算法

對于動力學(xué)系統(tǒng)而言，動力學(xué)方程可由第一類拉格朗日方程得出：

(30)

式中：M為系統(tǒng)質(zhì)量陣；q為系統(tǒng)廣義坐標(biāo)；F為系統(tǒng)彈性力；Q為系統(tǒng)廣義外力；Φ和Φq分別為系統(tǒng)約束方程以及對廣義坐標(biāo)的雅克比矩陣；λ為拉格朗日乘子.

動力學(xué)方程的求解采用廣義-α方法[16].動力學(xué)方程進(jìn)行時間離散后，需在在每個時間步長內(nèi)通過Newton-Raphson迭代求解方程組：

(31)

式中Δq和Δλ是廣義坐標(biāo)和拉格朗日乘子的增量，

(32)

2.3 SPR誤差估計方法

在等幾何分析中，將HB樣條基函數(shù)作為幾何離散和分析求解的基函數(shù)，可實現(xiàn)網(wǎng)格的局部細(xì)化.在實際計算中，對于一個求解模型，需要確定具體的細(xì)分區(qū)域.例如，對于具有孔結(jié)構(gòu)的模型，其孔周應(yīng)力集中，位移梯度變化較大，需要在孔周區(qū)域建立多個細(xì)分層.對于復(fù)雜結(jié)構(gòu)，加密區(qū)域無法事先獲得，而為了實現(xiàn)自適應(yīng)細(xì)化，誤差估計算法是其中的一個重要環(huán)節(jié).

等幾何分析的能量范數(shù)誤差定義為

(33)

Δi/Δ>const,

(34)

其中，參數(shù)const是一個常量，則該單元的誤差較大，需要進(jìn)一步細(xì)化.

通常，大變形結(jié)構(gòu)動力學(xué)問題無法獲取解析解，可采用超收斂修復(fù)方法(superconvergence patch recovery，SPR)修復(fù)應(yīng)力，可得到應(yīng)力的參考解.SPR將高斯點上的應(yīng)力值映射至控制點上，可得到整個單元修復(fù)應(yīng)力場，從而可獲得整體修復(fù)應(yīng)力場.將該修復(fù)的應(yīng)力場當(dāng)作參考解，從而可獲得系統(tǒng)的應(yīng)力誤差或能量誤差.

通常情況下，數(shù)值求解之后，僅可獲得高斯點處的應(yīng)力值，即應(yīng)力場的強形式為

(35)

將式(35)進(jìn)行變分，同時在單元內(nèi)或整個求解域積分，可得到對應(yīng)的弱形式為

(36)

離散后進(jìn)行高斯積分后可表示為

∑wNTNTnode=∑wNTTGauss,

(37)

因此，節(jié)點的應(yīng)力值可表示為

Tnode=Ψ/PGauss,

(38)

式中：Ψ=∑wNTN，PGauss=∑wNT，N為單元的形函數(shù)；w為高斯積分系數(shù).

求解域的應(yīng)力場可通過B樣條基函數(shù)插值得到：T(ξ,η)=S(ξ,η)Tnode，將該修復(fù)應(yīng)力場當(dāng)作視為應(yīng)力參考解，可獲得系統(tǒng)的能量范數(shù)誤差.

3 數(shù)值算例

3.1 圓孔線彈性問題

第一個算例是一個無限大平板的孔周應(yīng)力集中問題.本文研究了在平面應(yīng)力假設(shè)下的線彈性問題：

divσ(u)=f,

(39)

由于對稱性，計算域選取1/4區(qū)域，如圖7所示，該算例的仿真參數(shù)來源于文獻(xiàn)[8]中的第4.2小節(jié).該靜力學(xué)問題的解析解為

(40)

采用NURBS樣條基函數(shù)作為幾何離散的基函數(shù)，能夠精確描述圓孔構(gòu)型.層次細(xì)化是一個保證幾何不變的過程，通過簡單的節(jié)點插入算法就可以計算出NURBS樣條的新的權(quán)系數(shù)和控制點.初始構(gòu)型采用2階NURBS樣條描述，初始廣義坐標(biāo)個數(shù)是28，網(wǎng)格數(shù)是8，最終通過3次網(wǎng)格局部層次細(xì)化產(chǎn)生了一個4層自適應(yīng)網(wǎng)格，如圖8所示.

3.2 平面板的結(jié)構(gòu)動力學(xué)問題

本算例研究一大變形平面板的結(jié)構(gòu)動力學(xué)問題.一個矩形板在重力作用下振動，板的左端固支，如圖10所示，板的長度寬度以及厚度分別是2，1和1 m，密度是7 810 kg/m3，彈性模量10 MPa，泊松比0.3.仿真計算時間是2 s，時間積分步長是10-4s.

平面板結(jié)構(gòu)的初始構(gòu)型采用2階NURBS樣條描述，初始網(wǎng)格數(shù)是32，初始廣義坐標(biāo)個數(shù)是60，在進(jìn)行網(wǎng)格自適應(yīng)細(xì)化的過程中，采用HB樣條的網(wǎng)格細(xì)分方法，通過SPR方法修復(fù)應(yīng)力，得到能量范數(shù)誤差來指導(dǎo)自適應(yīng)過程.在平面板的振動過程中，在誤差較大的區(qū)域進(jìn)行網(wǎng)格細(xì)化，共進(jìn)行了3次網(wǎng)格細(xì)化，最后形成了4層HB樣條細(xì)分網(wǎng)格，最終網(wǎng)格數(shù)是368，初始廣義坐標(biāo)個數(shù)是420.圖11是該結(jié)構(gòu)不同細(xì)化程度的網(wǎng)格圖，1/8左邊界部分最終網(wǎng)格最密集.參考模型也采用NURBS樣條來建模計算求解，其網(wǎng)格數(shù)和初始廣義坐標(biāo)個數(shù)分別為368、450，它與求解模型中最密網(wǎng)格數(shù)和初始廣義坐標(biāo)個數(shù)基本一致.圖12給出了整個計算過程中系統(tǒng)的能量范數(shù)誤差，相同初始廣義坐標(biāo)個數(shù)時，均勻網(wǎng)格的參考模型的誤差峰值為12.33，而采用HB樣條建模計算的誤差峰值為9.20，減小了25.4%.圖13表示1/8左邊界區(qū)域能量范數(shù)誤差，該部分是誤差較大區(qū)域，從圖中可看出，誤差峰值從均勻網(wǎng)格的9.32減小到6.59，降低了35.5%，平均誤差從4.37減小到3.10，降低了29.0%.由此說明基于HB樣條能夠有效地提高大變形平面結(jié)構(gòu)動力學(xué)問題的計算效率與精度.

4 結(jié) 論

基于等幾何分析方法，研究了大變形結(jié)構(gòu)動力學(xué)問題的網(wǎng)格局部加密算法.采用HB樣條和多層貝塞爾提取方法分別對圓孔線彈性問題和平面大變形板結(jié)構(gòu)動力學(xué)問題進(jìn)行建模.分析表明：

① 采用多層貝塞爾提取技術(shù)，將HB樣條的單元形函數(shù)統(tǒng)一，簡化了等幾何分析程序，實現(xiàn)了等幾何分析程序與有限元程序的融合；

② 采用超收斂修復(fù)方法(SPR)得到的修復(fù)應(yīng)力解可獲得動力學(xué)系統(tǒng)的能量范數(shù)誤差，且基于HB樣條的自適應(yīng)網(wǎng)格技術(shù)獲得的誤差相比均勻網(wǎng)格誤差明顯降低，表明該方法能夠有效地提高大變形平面結(jié)構(gòu)動力學(xué)問題的計算效率與精度.