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

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

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

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

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

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

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

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

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

(1)

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

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

(2)

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

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

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

(3)

式中：ξ為量綱一參數(shù)坐標；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é)點矢量元素有關，這些元素組成的區(qū)間[ξi,ξi+p+1]稱為基函數(shù)Ni,p的支持域，記任意基函數(shù)f的支持域為suppf.B樣條基函數(shù)連續(xù)性由節(jié)點矢量元素重復度決定，對于重復度為k的節(jié)點，基函數(shù)Cp-k連續(xù)，而在節(jié)點之間C∞連續(xù).一般而言，對于開放式節(jié)點矢量，第1個結點和最后一個結點的重復度為k=p+1.

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

(5)

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

(6)

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

(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ù)域與幾何模型存在映射關系，這些單元把實際物理模型劃分為多個一一對應的子區(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ù)，如B樣條，NURBS以及T樣條.貝塞爾提取的核心思想是采用節(jié)點插入算法，實現(xiàn)伯恩斯坦基函數(shù)來表示全局定義下的B樣條基函數(shù).

(8)

其中，

(9)

為了保持節(jié)點插入中幾何構型不變，應滿足

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

(10)

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

N(ξ)=CB(ξ),

(11)

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

(12)

(13)

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

(14)

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

(15)

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

1.3 HB樣條構造

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

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

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

③K=KN-1.

1.4 多層貝塞爾提取

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

Nl=Ml,l+1Nl+1,

(16)

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

(17)

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

(18)

(19)

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

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

(20)

式中MLε3稱為多層細分算子，它是通過組合多個細分算子的行向量來推導得到.例如，單元的多層細分算子是通過組合R0,2，R1,2和R2,2中相應的行向量得到.通過多層細分算子，則層次B樣條的基函數(shù)可以表示成單一層標準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 大變形平面結構動力學建模方法

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

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

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

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單元控制點，稱為單元廣義坐標.為了提高計算效率，對公式(23)中的單元廣義坐標進行重排，點P的位置坐標可重新表示為[15]

(24)

? ?

(25)

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

(26)

(27)

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

(28)

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

(29)

2.2 大變形結構動力學方程時間積分算法

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

(30)

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

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

(31)

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

(32)

2.3 SPR誤差估計方法

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

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

(33)

Δi/Δ>const,

(34)

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

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

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

(35)

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

(36)

離散后進行高斯積分后可表示為

∑wNTNTnode=∑wNTTGauss,

(37)

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

Tnode=Ψ/PGauss,

(38)

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

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

3 數(shù)值算例

3.1 圓孔線彈性問題

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

divσ(u)=f,

(39)

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

(40)

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

3.2 平面板的結構動力學問題

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

平面板結構的初始構型采用2階NURBS樣條描述，初始網(wǎng)格數(shù)是32，初始廣義坐標個數(shù)是60，在進行網(wǎng)格自適應細化的過程中，采用HB樣條的網(wǎng)格細分方法，通過SPR方法修復應力，得到能量范數(shù)誤差來指導自適應過程.在平面板的振動過程中，在誤差較大的區(qū)域進行網(wǎng)格細化，共進行了3次網(wǎng)格細化，最后形成了4層HB樣條細分網(wǎng)格，最終網(wǎng)格數(shù)是368，初始廣義坐標個數(shù)是420.圖11是該結構不同細化程度的網(wǎng)格圖，1/8左邊界部分最終網(wǎng)格最密集.參考模型也采用NURBS樣條來建模計算求解，其網(wǎng)格數(shù)和初始廣義坐標個數(shù)分別為368、450，它與求解模型中最密網(wǎng)格數(shù)和初始廣義坐標個數(shù)基本一致.圖12給出了整個計算過程中系統(tǒng)的能量范數(shù)誤差，相同初始廣義坐標個數(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樣條能夠有效地提高大變形平面結構動力學問題的計算效率與精度.

4 結 論

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

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

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