求解粘彈性問題的時(shí)域自適應(yīng)等幾何比例邊界有限元法

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

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

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

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

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

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

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

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

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

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

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

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

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

1 遞推控制方程

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

邊界條件為：

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

式中，tk_1和Tk分別表示第k段時(shí)間間隔的起始點(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

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

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

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

其中：

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

其中：

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

在第一時(shí)段：

在第(k+1)時(shí)段：

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

3 遞推IG-SBFEM方程

3.1 NURBS基函數(shù)

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

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

當(dāng)p＝0時(shí)，

當(dāng)p≥1時(shí)，

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

式中，為對(duì)應(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

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

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

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

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

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

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

其中：

其中：

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

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

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

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

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

在第k個(gè)時(shí)段起點(diǎn)，

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

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

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

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

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

4 數(shù)值算例

4.1 半無限域蠕變問題

考慮一個(gè)受單向拉伸的帶圓孔的粘彈性無限域，受x方向的體力f= 100 N/m2，利用對(duì)稱性取如圖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]。

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

圖3 帶有1/4圓孔的半無限長(zhǎng)粘彈性體區(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)行計(jì)算，參考解由二次SBFEM得到的收斂解提供[32]。IG-SBFEM采用5個(gè)控制點(diǎn)，常規(guī)線性和二次SBFEM使用5個(gè)節(jié)點(diǎn)。計(jì)算結(jié)果顯示，在相同節(jié)點(diǎn)情況下，IG-SBFEM計(jì)算得到的蠕變曲線更加接近參考解。

采用相對(duì)偏差范數(shù)評(píng)估計(jì)算精度[33]：

其中：

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

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

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

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

表1給出了點(diǎn)A在不同時(shí)刻的位移解，在相同控制點(diǎn)/節(jié)點(diǎn)情況下，IG-SBFEM計(jì)算結(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)的位移解，描述了誤差限β對(duì)展開階次的影響。計(jì)算結(jié)果顯示，與相比，的誤差限需要更多的展開階數(shù)，可見所提方法可根據(jù)不同的精度要求自適應(yīng)地調(diào)整展開階數(shù)。

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

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

考慮一個(gè)裂紋的粘彈性平板，如圖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)所示)在不同的時(shí)間步長(zhǎng)下σy解的相對(duì)偏差比較。隨著時(shí)間步長(zhǎng)從0.005 s變化到0.01 s，不同時(shí)刻的應(yīng)力解相對(duì)偏差相差在1.6%以內(nèi)。

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

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

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

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

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

圖10 應(yīng)力隨時(shí)間變化曲線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ù)對(duì)稱性，選取如圖12所示的1/4結(jié)構(gòu)進(jìn)行分析，IG-SBFEM數(shù)值計(jì)算模型如圖12所示，其中控制點(diǎn)數(shù)量為 20個(gè)，在此基礎(chǔ)上進(jìn)行節(jié)點(diǎn)插入加密可得到48和97個(gè)控制點(diǎn)的二次NURBS離散形式[34]。參考解由ABAQUS的收斂解提供，計(jì)算模型如圖13所示，共使用49670個(gè)節(jié)點(diǎn)。

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

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

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

表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解與參考解的比較，計(jì)算結(jié)果顯示，所提方法與參考解符合良好。

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

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

2)圖6(a)和圖6(b)和表1的計(jì)算結(jié)果及結(jié)果表明，當(dāng)未知量數(shù)目相同時(shí)，與常規(guī)線性和二次SBFEM 相比，IG-SBFEM 具有更高的計(jì)算精度；3)由圖7和表2的計(jì)算結(jié)果可見，所提算法可對(duì)不同的時(shí)間步長(zhǎng)，通過自適應(yīng)計(jì)算，保持穩(wěn)定的計(jì)算精度，并通過展開階數(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ù)、比例邊界有限元、時(shí)域分段自適應(yīng)算法的優(yōu)點(diǎn)，提出了一種新的求解粘彈性問題的時(shí)域自適應(yīng)IG-SBFEM方法。所提方法將時(shí)空耦合的粘彈性性問題解耦為一系列遞推形式的空間彈性問題，并建立了基于IG-SBFEM的遞推求解方程。IG-SBFEM具有半解析性，可更準(zhǔn)確地描述問題的幾何模型，可將SBFEM與CAD模型無縫融合，并適于求解無限域和應(yīng)力奇異性相關(guān)的粘彈性問題。所提算法時(shí)域計(jì)算精度高，對(duì)不同的時(shí)間步長(zhǎng)可保持穩(wěn)定的計(jì)算精度。數(shù)值算例的結(jié)果表明，所提方法具有良好的計(jì)算精度和收斂性。