求解粘彈性問題的時域自適應等幾何比例邊界有限元法

何宜謙，王霄騰，祝雪峰，楊海天，薛齊文

(1.大連理工大學工程力學系/工業裝備結構分析國家重點實驗室，遼寧，大連 116024；2.大連理工大學汽車工程學院，遼寧，大連 116024；3.大連交通大學土木工程學院，遼寧，大連 116028)

許多材料具有粘彈性性質，如混凝土、高聚物、瀝青、陶瓷、生物組織等，粘彈性問題的求解具有重要的工程背景和理論探討價值[1－4]。由于時間相關的本構關系，加之復雜的邊界形狀/邊界條件，粘彈性問題通常難以解析求解，因此，發展準確有效的數值方法十分必要。

對時空耦合的粘彈性問題的數值求解，需同時考慮時間和空間的離散計算及計算精度?；诜侄畏e分/差分的時域粘彈性問題數值求解方法[5－8]，當積分/差分階數較低，計算精度受時間步長大小影響較大，而恰當的步長大小往往是預先難以確定的?；诜e分變換與對應性原理的數值方法[9－12]，通常需要進行數值的積分逆變換，其間可能產生的計算誤差將影響到時域的計算結果。

另一方面，可在有限元方法[13]、無網格方法[14]、邊界元方法[15]等求解粘彈性問題的基礎上，發展新的數值方法，以進一步提高粘彈性問題的空間數值求解精度。

等幾何分析[16](isogeometric analysis, IGA)基于有限元分析方法的等參單元思想，將計算機輔助幾何設計(CAGD)中用于表達幾何模型的非均勻有理B樣條(NURBS)的基函數作為場函數的形函數。由于NURBS基函數具有高階連續性，可精確描述各種幾何，因此基于NURBS的IGA具有良好的計算精度和收斂性[17]，并實現了CAD與CAE的無縫結合。

將等幾何分析技術用于粘彈性問題的求解，有望利用等幾何分析的以上優點，以進一步提高粘彈性問題的空間數值求解精度。但目前有關IGA在粘彈性數值求解中應用的文獻報道鮮有涉及。Erfan等[18]結合 IGA有限元和時域差分方法建立了一個粘彈性振動問題的數值計算模型，但限于一維問題且時域采用一階差分格式。

比例邊界有限元法(scaled boundary finite element method, SBFEM)是一種半解析的邊界類的數值方法[19]，適于處理無限域和應力奇異性問題，并已在粘彈性問題的數值求解中得到應用[20]。

林皋、宋崇民、劉俊、王文淵等將IGA和SBFEM相結合，建立了IG-SBFEM模型，并應用于靜力學、動力學、電磁場、結構地基相互作用、熱傳導、液體晃蕩等問題的求解[21－28]。但目前似尚未見到IG-SBFEM求解粘彈性問題的直接相關報道。

時域分段自適應算法是一種求解時空耦合問題的時域高精度數值算法，已與有限元、邊界元、無網格等空間數值方法相結合，用于瞬態熱傳導、動力學和粘彈性等問題的數值求解[29－30]。該算法的主要特點是可將一個時空耦合轉化為一系列具有遞推格式的空間問題，可在每個時間段內通過自適應計算，保證對不同時間步長保持穩定的計算精度。

基于以上考慮，本文將時域自適應算法和等幾何比例邊界元法(IG-SBFEM)相結合，提出了一種新的求解粘彈性問題的數值方法。利用時域分段自適應算法將粘彈性問題解耦，建立了遞推格式的IG-SBFEM計算模型，并采用自適應技術遞推求解，并將求解結果與解析解和參考解進行了比對。

本文算法在時域通過分段時域自適應計算，保證不同時間步長的計算精度；在空間具有徑向解析性，且可對幾何模型進行更精確描述。

論文章節安排如下：第1部分和第2部分分別推導了基于時域分段展開的粘彈性遞推控制方程和本構方程，第3部分推導了IG-SBFEM的求解格式，第4部分給出了三個算例驗證，第5部分為結論。

1 遞推控制方程

粘彈性問題的控制方程為：

邊界條件為：

將時域劃分為若干離散的時間段，在各離散時段內，將各變量展開為：

式中，tk_1和Tk分別表示第k段時間間隔的起始點和大小。

導數的轉換關系可為：

將式(5)~式(12)代入式(1)~式(4)中，原初邊值問題式(1)~式(4)轉化為以下具有遞推形式的邊值問題：

2 遞推本構方程

考慮圖1所示的三參數固體粘彈性模型，其微分本構關系為：

其中：

圖1 三參數固體模型Fig.1 Three-parameter solid model

對于平面應力問題：

對于平面應變問題：

將式(5)~式(6)和式(13)代入式(18)，并比較等式兩邊同次冪數的系數：

其中：

對于平面應力問題：

其中：

對于平面應變問題，需要將E1、E2和v分別替換成

在第一時段：

在第(k+1)時段：

其中，下標(k+1)和k分別表示第(k+1)和第k次時間間隔。

3 遞推IG-SBFEM方程

3.1 NURBS基函數

NURBS是對B樣條進行有理化構造后得到的有理樣條函數。在一維參數坐標空間中選定一組節點，記作設n為基函數個數，p為基函數最高次數。

一維B樣條基函數的遞推關系定義為：

當p＝0時，

當p≥1時，

通過引入權參數對 B樣條進行有理化，得到NURBS基函數：

式中，為對應的權參數。

利用下式，任意形狀曲線均可由NURBS基函數和控制點坐標精確表示：

3.2 遞推IG-SBFEM求解格式

采用如圖2所示的比例邊界坐標系，其中徑向坐標?比例中心處?＝0，在邊界處?＝1。環向坐標表示曲線到起始點的距離，比例邊界坐標系與笛卡爾坐標系的轉換關系如下：

圖2 比例邊界坐標系Fig.2 Scaled boundary coordinate system

對于一個SBFEM建模的子域邊界，在環向采用NURBS離散，節點向量為任一點坐標表示為[22]：

在比例坐標系統中：

對位移采用NURBS基函數離散，得到：

比例邊界坐標系中的應變算子為[31]：

在無體力的情況下，由虛功原理可得：

將式(24)代入式(41)可得：

其中：

其中：

式(42)~式(43)的解為[31]：

將式(50)代入式(42)、式(43)可得：

求解以上的特征值方程，可得到求解域的剛度矩陣為：

式中：[φ1]由n個特征向量組成；[Q1]由特征向量{q}組成[31]。

因此，邊界節點位移向量的m階求解方程為：

在第k個時段起點，

在各時間內的收斂準則為：

式中：β是規定的誤差限；表示L2-范數。

由于 NURBS邊界單元的控制點不具有插值性，本文采用一種利用NURBS基函數的單位分解性施加本質邊界條件的方法[23]，對于常位移邊界條件，將與邊界相關的控制點分成兩組，在上恒為零的和不恒為零的有：

對式(57)的非負基函數項利用 NURBS基函數的單位分解性，將求解域上的邊界約束條件轉化到控制點上：

與上式相對應的遞推邊界條件為：

4 數值算例

4.1 半無限域蠕變問題

考慮一個受單向拉伸的帶圓孔的粘彈性無限域，受x方向的體力f= 100 N/m2，利用對稱性取如圖3所示的1/4區域分析，位移邊界條件如圖3所示。粘彈性材料本構采用三參數固體模型，并采用巖體材料參數E1＝1 960N/m2，E2＝9 800 N/m2，?1＝5 2083(Pa·d)/m2[20]。

考慮到本文算法將時間相關的粘彈性問題轉化為一系列彈性問題求解如式(53)所示，因此，本文算法的計算精度依賴于彈性問題的IG-SBFEM計算精度。以文獻[32]中式(44)的能量誤差范數?E為指標，將IG-SBFEM計算結果與文獻[32]的計算結果相比較，計算結果如圖4所示，對于本算例，在相同的自由度情況下，IG-SBFEM比線性和二次的SBFEM具有更高的計算精度。

圖3 帶有1/4圓孔的半無限長粘彈性體區域Fig.3 A semi-unbounded viscoelastic plate with 1/4 circular hole

圖4 應力誤差范數隨自由度的變化曲線Fig.4 The variation of stress error norm with DOF

圖5給出了點A的蠕變曲線，分別采用線性和二次SBFEM和IG-SBFEM進行計算，參考解由二次SBFEM得到的收斂解提供[32]。IG-SBFEM采用5個控制點，常規線性和二次SBFEM使用5個節點。計算結果顯示，在相同節點情況下，IG-SBFEM計算得到的蠕變曲線更加接近參考解。

采用相對偏差范數評估計算精度[33]：

其中：

圖5 A點蠕變曲線Fig.5 Creep curve of point A

圖6(a)和圖6(b)通過A點在不同時刻的位移解，比較了SBFEM與IG-SBFEM在不同自由度下的考慮1＜?＜3區域內的偏差范數曲線，采用平均斜率計算收斂階次mc。計算結果顯示，IG-SBFEM具有更高的計算精度和更快的收斂性。

圖6(c)比較了SBFEM與IG-SBFEM在三種自由度下的計算時間，結果顯示，在相同自由度下，IG-SBFEM和SBFME的計算耗時相差不大。

圖6 偏差范數和計算時間隨自由度的變化曲線Fig.6 The variation curve of the error norm and computation time with DOF

表1給出了點A在不同時刻的位移解，在相同控制點/節點情況下，IG-SBFEM計算結果的最大偏差僅為0.5068%，而二次SBFEM和線性SBFEM的最大偏差分別為1.6312%和5.5051%。

表1 位移數值解與參考解的比較Table 1 Comparison of numerical solutions and reference solution for displacement

圖7通過A點的位移解，描述了誤差限β對展開階次的影響。計算結果顯示，與相比，的誤差限需要更多的展開階數，可見所提方法可根據不同的精度要求自適應地調整展開階數。

圖7 展開階次隨時間的變化Fig.7 The variation of expansion order with time

4.2 帶裂紋粘彈性平板的應力松弛問題

考慮一個裂紋的粘彈性平板，如圖8(a)所示，h＝b＝1 m ，c＝0.15m ，平板上側施加常位移邊界條件u0＝5×10-4m ，粘彈性材料本構采用三參數固體模型，并采用混凝土材料參數E1＝1 9600 N/m2，E2＝1 9600N/m2，?1＝6 35420(Pa·d)/m2[17]。IG-SBFEM數值模型如圖8(b)所示。采用ABAQUS軟件的收斂解作為參考解，節點總數為 12343，單元分布如圖8(c)所示。

表2給出了點A(如圖8(a)所示)在不同的時間步長下σy解的相對偏差比較。隨著時間步長從0.005 s變化到0.01 s，不同時刻的應力解相對偏差相差在1.6%以內。

圖8 帶裂紋平板Fig.8 A plate with a crack

表2 不同時間步長下A點IG-SBFEM應力解的比較Table 2 Comparion for IG-SBFEM displacement solutions at point A with different time steps

圖9給出了本文算法和時域非自適應算法的比較，其中非自適應計算由ABAQUS軟件提供。計算結果表明，當時間步長由0.01增加到0.1時，非自適應算法的計算結果明顯偏離于參考解，而本文所提的自適應算法的計算結果受時間步長變化的影響很小。

圖10給出了點A、B和C的應力松弛曲線，計算結果顯示，IG-SBFEM的數值解與參考解符合良好。

圖9 本文算法和時域非自適應算法的比較Fig.9 The comparison of the proposed model with non adaptive algorithm in time domain

圖10 應力隨時間變化曲線Fig.10 The variation of stress with time

4.3 粘彈性巖體中的襯砌結構蠕變分析

考慮如圖11所示的在無限域粘彈性巖體中的混凝土襯砌，襯砌內部受到作用于內壁的徑向均布荷載P=1 N/m2，幾何參數a=1.2 m，b=1.0 m，混凝土襯砌和巖體材料的三參數固體粘彈性參數[17]：

式中：上標C表示混凝土襯砌；R表示巖體。

圖11 無限域粘彈性巖體中的混凝土襯砌Fig.11 Circular concrete lining in unbounded viscoelastic rock

根據對稱性，選取如圖12所示的1/4結構進行分析，IG-SBFEM數值計算模型如圖12所示，其中控制點數量為 20個，在此基礎上進行節點插入加密可得到48和97個控制點的二次NURBS離散形式[34]。參考解由ABAQUS的收斂解提供，計算模型如圖13所示，共使用49670個節點。

圖12 IG-SBFEM數值計算模型(20個控制點)Fig.12 IG-SBFEM numerical model(20 control points)

圖13 有限元網格Fig.13 FEM mesh

表3給出了不同控制點數目對點A在x方向位移的影響，隨著控制點數目的增加，數值解的相對偏差逐漸減小。

表3 不同控制點條件下的A點x軸方向位移解Table 3 Displacement solution in x-direction at point A with the different control points

圖14給出了A點和B兩點蠕變曲線IG-SBFEM解與參考解的比較，計算結果顯示，所提方法與參考解符合良好。

4.4 計算結果小結

1)圖4、圖10、圖14以及表1、表2的計算結果及比較表明，所提算法可對粘彈性問題進行準確有效的數值求解；

2)圖6(a)和圖6(b)和表1的計算結果及結果表明，當未知量數目相同時，與常規線性和二次SBFEM 相比，IG-SBFEM 具有更高的計算精度；3)由圖7和表2的計算結果可見，所提算法可對不同的時間步長，通過自適應計算，保持穩定的計算精度，并通過展開階數的調整滿足不同大小誤差限要求。

圖14 A點和B點的蠕變曲線Fig.14 The creep curves for points A and B

5 結論

本文的主要貢獻是集成等幾何分析技術、比例邊界有限元、時域分段自適應算法的優點，提出了一種新的求解粘彈性問題的時域自適應IG-SBFEM方法。所提方法將時空耦合的粘彈性性問題解耦為一系列遞推形式的空間彈性問題，并建立了基于IG-SBFEM的遞推求解方程。IG-SBFEM具有半解析性，可更準確地描述問題的幾何模型，可將SBFEM與CAD模型無縫融合，并適于求解無限域和應力奇異性相關的粘彈性問題。所提算法時域計算精度高，對不同的時間步長可保持穩定的計算精度。數值算例的結果表明，所提方法具有良好的計算精度和收斂性。