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

(1.大連理工大學(xué)工程力學(xué)系/工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，遼寧，大連 116024；2.大連理工大學(xué)汽車工程學(xué)院，遼寧，大連 116024；3.大連交通大學(xué)土木工程學(xué)院，遼寧，大連 116028)

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

對時空耦合的粘彈性問題的數(shù)值求解，需同時考慮時間和空間的離散計算及計算精度。基于分段積分/差分的時域粘彈性問題數(shù)值求解方法[5－8]，當(dāng)積分/差分階數(shù)較低，計算精度受時間步長大小影響較大，而恰當(dāng)?shù)牟介L大小往往是預(yù)先難以確定的?；诜e分變換與對應(yīng)性原理的數(shù)值方法[9－12]，通常需要進(jìn)行數(shù)值的積分逆變換，其間可能產(chǎn)生的計算誤差將影響到時域的計算結(jié)果。

另一方面，可在有限元方法[13]、無網(wǎng)格方法[14]、邊界元方法[15]等求解粘彈性問題的基礎(chǔ)上，發(fā)展新的數(shù)值方法，以進(jìn)一步提高粘彈性問題的空間數(shù)值求解精度。

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

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

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

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

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

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

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

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

1 遞推控制方程

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

邊界條件為：

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

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

導(dǎo)數(shù)的轉(zhuǎn)換關(guān)系可為：

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

2 遞推本構(gòu)方程

考慮圖1所示的三參數(shù)固體粘彈性模型，其微分本構(gòu)關(guān)系為：

其中：

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

對于平面應(yīng)力問題：

對于平面應(yīng)變問題：

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

其中：

對于平面應(yīng)力問題：

其中：

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

在第一時段：

在第(k+1)時段：

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

3 遞推IG-SBFEM方程

3.1 NURBS基函數(shù)

NURBS是對B樣條進(jìn)行有理化構(gòu)造后得到的有理樣條函數(shù)。在一維參數(shù)坐標(biāo)空間中選定一組節(jié)點(diǎn)，記作設(shè)n為基函數(shù)個數(shù)，p為基函數(shù)最高次數(shù)。

一維B樣條基函數(shù)的遞推關(guān)系定義為：

當(dāng)p＝0時，

當(dāng)p≥1時，

通過引入權(quán)參數(shù)對 B樣條進(jìn)行有理化，得到NURBS基函數(shù)：

式中，為對應(yīng)的權(quán)參數(shù)。

利用下式，任意形狀曲線均可由NURBS基函數(shù)和控制點(diǎn)坐標(biāo)精確表示：

3.2 遞推IG-SBFEM求解格式

采用如圖2所示的比例邊界坐標(biāo)系，其中徑向坐標(biāo)?比例中心處?＝0，在邊界處?＝1。環(huán)向坐標(biāo)表示曲線到起始點(diǎn)的距離，比例邊界坐標(biāo)系與笛卡爾坐標(biāo)系的轉(zhuǎn)換關(guān)系如下：

圖2 比例邊界坐標(biāo)系Fig.2 Scaled boundary coordinate system

對于一個SBFEM建模的子域邊界，在環(huán)向采用NURBS離散，節(jié)點(diǎn)向量為任一點(diǎn)坐標(biāo)表示為[22]：

在比例坐標(biāo)系統(tǒng)中：

對位移采用NURBS基函數(shù)離散，得到：

比例邊界坐標(biāo)系中的應(yīng)變算子為[31]：

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

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

其中：

其中：

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

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

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

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

因此，邊界節(jié)點(diǎn)位移向量的m階求解方程為：

在第k個時段起點(diǎn)，

在各時間內(nèi)的收斂準(zhǔn)則為：

式中：β是規(guī)定的誤差限；表示L2-范數(shù)。

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

對式(57)的非負(fù)基函數(shù)項(xiàng)利用 NURBS基函數(shù)的單位分解性，將求解域上的邊界約束條件轉(zhuǎn)化到控制點(diǎn)上：

與上式相對應(yīng)的遞推邊界條件為：

4 數(shù)值算例

4.1 半無限域蠕變問題

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

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

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

圖4 應(yīng)力誤差范數(shù)隨自由度的變化曲線Fig.4 The variation of stress error norm with DOF

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

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

其中：

圖5 A點(diǎn)蠕變曲線Fig.5 Creep curve of point A

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

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

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

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

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

圖7通過A點(diǎn)的位移解，描述了誤差限β對展開階次的影響。計算結(jié)果顯示，與相比，的誤差限需要更多的展開階數(shù)，可見所提方法可根據(jù)不同的精度要求自適應(yīng)地調(diào)整展開階數(shù)。

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

4.2 帶裂紋粘彈性平板的應(yīng)力松弛問題

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

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

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

表2 不同時間步長下A點(diǎn)IG-SBFEM應(yīng)力解的比較Table 2 Comparion for IG-SBFEM displacement solutions at point A with different time steps

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

圖10給出了點(diǎn)A、B和C的應(yīng)力松弛曲線，計算結(jié)果顯示，IG-SBFEM的數(shù)值解與參考解符合良好。

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

圖10 應(yīng)力隨時間變化曲線Fig.10 The variation of stress with time

4.3 粘彈性巖體中的襯砌結(jié)構(gòu)蠕變分析

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

式中：上標(biāo)C表示混凝土襯砌；R表示巖體。

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

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

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

圖13 有限元網(wǎng)格Fig.13 FEM mesh

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

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

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

4.4 計算結(jié)果小結(jié)

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

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

圖14 A點(diǎn)和B點(diǎn)的蠕變曲線Fig.14 The creep curves for points A and B

5 結(jié)論

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