李佳龍，李 鋼，余丁浩

(大連理工大學海岸和近海工程國家重點實驗室，遼寧，大連 116024)

比例邊界有限元法(Scaled boundary finite element method，簡稱SBFEM)是由Wolf 和Song[1?2]提出的一種半解析的數(shù)值求解方法，最早用于解決土-結(jié)構(gòu)相互作用問題中的無限域動力剛度的計算[3]。多邊形比例邊界有限元法[4]結(jié)合了有限單元法與邊界元法的優(yōu)勢，具有較高的計算精度以及收斂性[5]，目前已運用到眾多工程問題的數(shù)值分析中[6?9]。多邊形比例邊界單元由于其網(wǎng)格的靈活性，在模擬裂紋擴展過程、處理局部網(wǎng)格重剖分等方面相較于有限元單元具有更明顯的優(yōu)勢。目前，比例邊界有限元法更多關(guān)注的是線彈性問題的求解，對于非線性問題的求解則研究較少。通過對多邊形比例邊界單元的彈性特性進行分析，可獲得與有限單元法等價的形函數(shù)以及應(yīng)變-位移矩陣，且對于線彈性問題與非線性問題均適用?；诖?，Ooi 等[10]開創(chuàng)了多邊形比例邊界單元的非線性分析，并將其運用于彈塑性帶裂紋問題的求解，其中單元的彈塑性本構(gòu)矩陣采用最小二乘法進行數(shù)值擬合，但該方法計算過程復雜且要求多邊形單元尺寸足夠小才能夠精確表達塑性區(qū)的變化。Chen等[11]提出在每個線單元覆蓋的扇形區(qū)域引入3 個數(shù)值積分點來進行多邊形比例邊界單元的非線性分析，具有較高的計算精度。無論是Ooi 等[10]亦或是Chen[11]等所提出的多邊形比例邊界有限元非線性分析過程，兩者的基本思路均是：先計算多邊形比例邊界單元的彈塑性剛度矩陣，然后再集成整個求解區(qū)域的總體剛度矩陣。然而，同傳統(tǒng)有限單元法一樣，均是對總體剛度矩陣進行分解求解從而獲得單元的位移響應(yīng)，但大規(guī)模剛度矩陣的分解依舊是限制非線性問題高效求解最主要的因素，因而有必要提出一種針對多邊形比例邊界單元的高效求解方法。

隔離非線性有限元法(Inelasticity-separated finite element method，簡稱IS-FEM)是由Li 和Yu[12]基于有限單元法的基本框架提出的一種結(jié)構(gòu)非線性分析新方法，采用Woodbury 公式[13]直接求解控制方程，實現(xiàn)了局部非線性問題的高效求解。該方法將單元應(yīng)變分解為線彈性應(yīng)變與非線性應(yīng)變兩部分，非線性應(yīng)變通過在單元中設(shè)置非線性應(yīng)變插值點以插值的形式來表示，從而將結(jié)構(gòu)的非線性部分隔離開來。對于非局部非線性問題，采用Woodbury 近似法[14?15]對控制方程進行求解，將非線性問題的迭代求解過程轉(zhuǎn)換為多次的常剛度矩陣回代以及稀疏矩陣與向量的乘積，對大規(guī)模非線性問題的求解具有明顯的效率優(yōu)勢。目前已成功運用于鋼筋混凝土框架結(jié)構(gòu)[16]、二維[17]以及三維[18]等問題的非線性分析，顯示出其較高的計算效率。

本文將隔離非線性有限元法的思想用于比例邊界有限元非線性分析，提出一種高效的隔離非線性比例邊界有限元法(Inelasticity-separated scaled boundary finite element method，簡稱IS-SBFEM)。主要思想是：將多邊形比例邊界單元所有扇形區(qū)域獨立考慮，在每個扇形單元內(nèi)設(shè)置若干個非線性應(yīng)變插值點，從而建立相互獨立的非線性應(yīng)變場。每個扇形區(qū)域的徑向自然坐標以及環(huán)向自然坐標分別為ξ=[0,1]、η=[?1,1]，因此本文針對多邊形比例邊界單元的扇形區(qū)域提出一種新的數(shù)值積分方案，即：采用4 點高斯積分方案對線性的扇形子單元進行數(shù)值積分。對于更高階的比例邊界單元，選取更多的數(shù)值積分點進行數(shù)值積分即可獲得足夠的精度。該數(shù)值求解方案對于線彈性問題以及非線性問題均適用，且精度保持不變。由于隔離非線性比例邊界單元相較于隔離非線性有限元單元引入了更多的非線性應(yīng)變插值點，即便是局部非線性問題其舒爾補矩陣維數(shù)也往往較大，因此建議直接采用Woodbury 近似法對控制方程進行求解，但該方法對較小規(guī)模的非線性問題計算效率偏低。鑒于此，對于較小規(guī)模的非線性問題，建議采用多邊形比例邊界有限元法直接求解；對于大規(guī)模的非線性問題，采用隔離非線性比例邊界有限元法求解具有更高的計算效率。基于MATLAB 平臺分別編制了比例邊界有限元法以及隔離非線性比例邊界有限元法兩套非線性分析程序，數(shù)值算例驗證了本文所提出的隔離非線性比例邊界有限元法的正確性與高效性。

1 SBFEM 基本理論

1.1 多邊形比例邊界有限元

在求解任意的平面問題時，可采用滿足比例中心要求的多邊形比例邊界單元進行離散。圖1為一典型的多邊形比例邊界單元，連接比例中心O與外邊界線單元端部結(jié)點可將其劃分為多個扇形子單元。邊界線單元可采用一維線單元進行離散，其自然坐標為?1≤η≤1，因而邊界線單元上的任意一點坐標(xb(η),yb(η))可表示為：

式中：xb=[x1,x2, ···,xm]T與yb=[y1,y2, ···,ym]T為該邊界線單元上m個結(jié)點的坐標；N(η)為一維Gauss-Lobatto-Lagrange 形函數(shù)。

圖1 多邊形比例邊界單元Fig. 1 Polygon scaled boundary element

在比例中心沿邊界的徑向方向引入自然坐標ξ，其坐標范圍為0≤ξ≤1，其中：ξ=0 為比例中心位置，ξ=1 為邊界位置。因而該求解域內(nèi)任意點坐標(x(ξ,η),y(ξ,η))可表示為：

比例邊界有限元法假定單元位移在環(huán)向方向為數(shù)值解，而在徑向方向則為一解析解。將其在環(huán)向與徑向兩個方向上解耦，可表示為：

式中：Nu(η)為邊界線單元的插值形函數(shù)；u(ξ)為徑向方向的位移函數(shù)，u(ξ)的求解可參考比例邊界有限元法相關(guān)文獻[1 ? 2，10 ? 11]，此處不再贅述。通過對多邊形比例邊界單元的彈性特性進行分析，可獲得每個扇形區(qū)與有限元單元等價的單元形函數(shù)Φ(ξ,η)以及應(yīng)變-位移矩陣B(ξ,η)。兩者僅與單元本身形狀相關(guān)而與材料的性質(zhì)無關(guān)，對線彈性問題以及非線性問題均適用，分別為：

式中：ψu為位移模態(tài)對應(yīng)的轉(zhuǎn)換矩陣；Sn是由單元的Hamilton 矩陣的負特征值組成的對角矩陣；B1(η)與B2(η)為單元由笛卡爾坐標到局部坐標的轉(zhuǎn)換矩陣。此時，單元的應(yīng)變場ε(ξ,η)為：

式中，ub為多邊形比例邊界單元的結(jié)點位移。

1.2 非線性比例邊界有限元

比例邊界有限元法作為一種基于位移的有限元求解方法，采用最小位能原理建立其非線性求解的控制方程。對于滿足平衡條件的變形體結(jié)構(gòu)有：

式中：ft與fb分別為單元邊界載荷與體積載荷；σ(ξ,η)為單元上一個荷載步收斂后的應(yīng)力；δε 為虛位移場δub對應(yīng)的虛應(yīng)變場，且有δε=B(ξ,η)δub。單元的應(yīng)力增量Δσ(ξ,η)可通過單元的彈塑性矩陣Dep以及結(jié)點位移增量Δub表示為：

將式(5)、式(9)代入式(8)所示的虛功方程，即可得到多邊形比例邊界單元非線性求解的控制方程，即：

式中，kep、Δf分別為多邊形比例邊界單元的彈塑性剛度矩陣以及外荷載增量。

多邊形比例邊界單元的剛度矩陣kep與外荷載Δf的具體計算步驟為：首先，通過式(5)、式(6)獲得多邊形單元每一個扇形區(qū)的形函數(shù)Φ(ξ,η)以及應(yīng)變-位移矩陣B(ξ,η)，然后，通過式(11)、式(12)計算該扇形區(qū)對單元剛度矩陣以及外荷載的貢獻，最后，再集成該多邊形單元的剛度矩陣以及外荷載向量[11]。通過對所有的多邊形單元進行計算并按自由度組裝獲得計算域的總體剛度矩陣以及外荷載向量，即可得到整個計算域的控制方程：

式中：Kep即為整個計算域的彈塑性剛度矩陣；ΔX與ΔF分別為總的結(jié)點位移增量與外荷載向量；Ne為整個計算域多邊形單元個數(shù)。采用上述方法，可將半解析的比例邊界有限元法轉(zhuǎn)換成與有限單元法一樣的純數(shù)值求解方法，從而實現(xiàn)了多邊形比例邊界單元的非線性分析。

對于式(13)中的剛度矩陣Kep同樣具有稀疏性，可采用LDLT分解法直接進行求解。然而，對于大規(guī)模問題的求解其計算效率偏低，因而有必要提出一種針對多邊形比例邊界單元的高效求解方法。

2 隔離非線性比例邊界有限元法

2.1 應(yīng)變分解與非線性應(yīng)變插值

隔離非線性有限元法(IS-FEM)[12]作為一種高效的材料非線性分析方法，已在二維、三維等有限元問題中展現(xiàn)其較高的計算效率，在此將其推廣至多邊形比例邊界單元的非線性分析，提出一種高效的隔離非線性比例邊界有限元法(IS-SBFEM)。對于一個多邊形比例邊界單元，采用增量形式，每個扇形區(qū)域內(nèi)任意點處的應(yīng)變增量Δε(ξ,η)可通過該扇形區(qū)的應(yīng)變矩陣B(ξ,η)以及單元邊界結(jié)點位移增量Δub表示為：

基于彈塑性力學的基本理論，單元內(nèi)任意點處的應(yīng)變增量Δε(ξ,η)均可分解為線彈性的應(yīng)變增

2.2 控制方程

比例邊界有限元法作為一種基于位移的有限元分析方法，采用虛功原理可建立隔離非線性比例邊界有限元法的控制方程。

式中，δ(Δε(ξ,η))為虛結(jié)點位移δ(Δub)對應(yīng)的虛應(yīng)變，即：δ(Δε(ξ,η))=B(ξ,η)δ(Δub)。其中，虛應(yīng)變δ(Δε(ξ,η))同樣可以分解為虛線彈性應(yīng)變δ(Δε′(ξ,η))以及虛非線性應(yīng)變δ(Δε″(ξ,η))兩部分，即δ(Δε(ξ,η)) =δ(Δε′(ξ,η)) + δ(Δε″(ξ,η))，且同樣滿足式(15)所示的非線性應(yīng)變插值關(guān)系：δ(Δε″(ξ,η))=Csδ( ?ε′s′)。將式(18)、式(19)兩式所示的應(yīng)力增量關(guān)系代入式(20)，即可得到隔離非線性比例邊界單元的控制方程：

2.3 多邊形單元數(shù)值積分

一個多邊形單元可拆分為ns個扇形區(qū)域，因此，式(22)所示的多邊形比例邊界單元的剛度矩陣ke的計算過程為：先單獨對每個扇形區(qū)域Ωs進行積分運算，然后再對所有的扇形區(qū)域進行集成。每個扇形子單元積分區(qū)域Ωs的自然坐標均為ξ=[0,1]與η=[?1,1]，且有dΩ=hξ|J(η)|dηdξ。因此，初始彈性剛度矩陣ke可表示為：

式中：ke,s為第s個扇形區(qū)域?qū)φ麄€多邊形單元剛度的貢獻；h為單元厚度；|J(η)|為邊界線單元的雅克比矩陣。在扇形區(qū)域內(nèi)引入g個高斯積分點，矩陣ke,s可采用高斯積分方案進行數(shù)值積分，有：

式中，ωi為積分權(quán)重。當多邊形比例邊界單元的階數(shù)較低時，可采用表1 所示4 點高斯積分方案進行計算，圖2 所示為多邊形比例邊界單元的數(shù)值積分點布置圖，其中線單元上的高斯積分點僅用來計算多邊形單元的形函數(shù)Φ(ξ,η)以及應(yīng)變-位移矩陣B(ξ,η)。對于高階的多邊形比例邊界單元，可設(shè)置更多的高斯積分點。

表1 扇形單元積分點坐標Table 1 Integral point coordinates of sector element

圖2 多邊形單元高斯點布置方案Fig. 2 Gaussian points layout scheme of polygon element

該積分方案對式(11)的計算同樣適用，僅需將式(26)中的彈性矩陣De更換為彈塑性矩陣Dep，即可得到多邊形比例邊界單元的彈塑性剛度矩陣kep。多邊形比例邊界單元的彈塑性矩陣Dep的計算可直接調(diào)用有限元程序中的彈塑性矩陣計算模塊，從而完成多邊形比例邊界單元的非線性分析。

式(23)與式(24)所示的單元矩陣k′以及k″的計算同式(22)一致，均是先通過高斯積分方案計算每個扇形區(qū)域的值，然后再進行集成。在對第s個扇形區(qū)進行數(shù)值積分時，由于假定每個扇形區(qū)域非線性應(yīng)變場相互獨立，則此時僅有該扇形區(qū)插值矩陣Cs存在，其余扇形區(qū)域的插值矩陣均為零矩陣。因而單元矩陣k′與k″是以扇形區(qū)為單元的分塊矩陣以及分塊對角矩陣。

對于較低階的比例邊界單元，每個扇形區(qū)可采用圖2 所示的4 個高斯積分點作為單元的非線性應(yīng)變插值點。由于插值矩陣Cs為Kronecker delta 函數(shù)，因此式(30)、式(31)中的矩陣分別是以插值點為單位的分塊矩陣以及分塊對角矩陣。

其中，子矩陣分別為：

3 控制方程求解

對于一個ns邊形的隔離非線性比例邊界單元，其內(nèi)部共有4ns個非線性應(yīng)變插值點，因而其舒爾補矩陣KS規(guī)模較大，直接求解的計算量可能比傳統(tǒng)的剛度矩陣直接分解求解的計算量還大。在此，可采用Woodbury 近似法[14?15](Woodbury approximation approach)對式(39)進行近似求解，整個求解過程僅為多次的常剛度矩陣回代以及稀疏矩陣與向量的乘積。該方法對于大規(guī)模非線性問題的求解具有明顯優(yōu)勢，且規(guī)模越大其計算效率相對于直接的剛度矩陣分解求解算法效率更高。文獻[17]已從理論上說明了對于平面有限元問題，當結(jié)點自由度數(shù)目超過2000 時即可體現(xiàn)出近似隔離非線性有限元法的高效性。

專職輔導員在專業(yè)知識和能力方面的缺乏，難以對學生做出專業(yè)的引導和相關(guān)職業(yè)的規(guī)劃。培養(yǎng)專業(yè)輔導員既彌補了此方面的不足，又對專業(yè)課教師參與學生思政工作起到了拋磚引玉的作用。

基于本文所提出的隔離非線性比例邊界有限元模型，編制了相應(yīng)的計算程序，整體分析流程如圖3 所示。

圖3 隔離非線性比例邊界有限元法分析流程Fig. 3 Analysis flow of IS-SBFEM

4 數(shù)值算例

4.1 懸臂梁靜力分析

圖4 為一懸臂梁模型，其幾何尺寸如圖所示，厚度為0.1 m。材料參數(shù)為：彈性模量E=2.0×105MPa，泊松比為ν=0.2。采用多邊形網(wǎng)格以及線性四邊形有限元單元進行網(wǎng)格剖分，其中多邊形網(wǎng)格共計109 個單元以及264 個結(jié)點，如圖5所示。有限元網(wǎng)格采用三種網(wǎng)格尺寸，分別為80 個單元與105 個結(jié)點、320 個單元與369 個結(jié)點以及1280 個單元與1377 個結(jié)點，分別如圖6(a)、圖6(b)以及圖6(c)所示。

圖4 懸臂梁模型Fig. 4 Cantilever beam model

圖5 多邊形網(wǎng)格Fig. 5 Polygon mesh

圖6 有限元網(wǎng)格Fig. 6 Finite element mesh

4.1.1 彈性分析

在自由端結(jié)點A處施加豎向荷載F=100 kN，并假定網(wǎng)格最密的有限元模型c 為參考解。表2給出了結(jié)點A在不同網(wǎng)格模型下橫向和豎向位移，可以看出采用本文提出的數(shù)值積分方案所得的計算結(jié)果與半解析的SBFEM 結(jié)果基本一致，且與有限元模型b 結(jié)果也基本一致。

表2 A 點處的位移Table 2 Displacement of node A

4.1.2 彈塑性分析

假定該懸臂梁為雙線性隨動強化彈塑性模型，其初始屈服應(yīng)力σy=360 MPa，切線模量Eh=0.002E。自由端結(jié)點A處施加往復豎向位移u，每個荷載步附加位移約束增量為d0=0.2 mm，共計2400 個荷載步[19]。采用SBFEM 以及IS-SBFEM兩種求解方案，其中隔離非線性方案采用直接分解舒爾補矩陣的精確求解以及近似求解兩種方法，共計四種求解方案。

通過彈性分析可知，有限元模型b 的剛度與多邊形模型基本一致，因而選擇模型b 作為對比模型來驗證非線性分析的精度。圖7 為兩種網(wǎng)格模型結(jié)點A處的豎向荷載F與結(jié)點位移u的變化曲線，可以看出四種求解方案結(jié)果基本吻合，多邊形單元的三種求解方案之間相對誤差基本可以忽略。多邊形模型與有限元模型計算結(jié)果一致，表明本文針對多邊形比例邊界單元所提出的數(shù)值積分方案以及IS-SBFEM 的正確性。圖8 給出了IS-SBFEM 與SBFEM 兩種求解方案的計算自由度曲線，比例邊界有限元模型的剛度矩陣維數(shù)為510，而隔離非線性比例邊界有限元模型的計算自由度隨著非線性程度的增加與減少呈現(xiàn)高低變化。

圖7 豎向荷載-位移曲線Fig. 7 Curves of the vertical load vs. displacement

圖8 計算自由度曲線Fig. 8 Calculate DOFs curves

表3 給出了多邊形網(wǎng)格三種求解方案控制方程的具體求解時間，可以看出SBFEM 求解效率最高，近似IS-SBFEM 法次之，精確IS-SBFEM 最低。這是因為當問題規(guī)模較小時，剛度矩陣的分解與回代耗時兩者差別不大，且僅需對剛度矩陣進行一次分解與回代即可得到當前迭代步的計算結(jié)果。近似的IS-SBFEM 由于要進行多次的常剛度矩陣回代以及稀疏矩陣與向量的乘積，因而其計算效率相對于直接的剛度矩陣分解法反而更低。采用精確的IS-SBFEM 進行求解時，由于舒爾補矩陣為高階滿陣且遠大于剛度矩陣維數(shù)，因而其求解效率最低。

表3 控制方程計算時間Table 3 Calculate time of the governing equation

4.2 Koyna 重力壩動力非線性分析

印度Koyna 重力壩作為少數(shù)幾個在強震中破壞且較有完整記錄的重力壩之一，長期以來均被科研工作者奉為經(jīng)典算例[20]。該壩長850 m，高103 m，符合平面應(yīng)變模型幾何特征，其具體幾何尺寸如圖9(a)所示。根據(jù)實際損傷以及相關(guān)文獻可知，壩踵以及坡折面處均發(fā)生破壞，因而對這兩處的網(wǎng)格進行局部加密，多邊形有限元模型如圖9(b)所示。其中壩踵以及坡折面處的網(wǎng)格細化圖分別如圖9(c)以及圖9(d)所示，整個模型共計36904 個多邊形單元以及18753 個結(jié)點。

圖9 幾何模型與多邊形網(wǎng)格Fig. 9 Geometry model and polygon mesh

圖10 Koyna 地震動記錄Fig. 10 Seismic acceleration records of Koyna

圖11～圖16 分別為壩頂橫向與豎向位移、速度以及加速度時程對比曲線。其中有限元模型采用線性四邊形等參單元進行網(wǎng)格劃分，整個計算區(qū)域共剖分103081 個結(jié)點以及102400 個單元?？梢钥闯?，F(xiàn)EM、SBFEM 以及IS-SBFEM 三種算法的計算結(jié)果完全吻合，再次驗證了本文所提出的隔離非線性比例邊界有限元法的正確性。

圖11 壩頂橫向位移反應(yīng)Fig. 11 Horizontal displacement responses of dam crest

圖12 壩頂豎向位移反應(yīng)Fig. 12 Vertical displacement responses of dam crest

圖13 壩頂橫向速度反應(yīng)Fig. 13 Horizontal velocity responses of dam crest

圖14 壩頂豎向速度反應(yīng)Fig. 14 Vertical velocity responses of dam crest

圖17 給出了該重力壩模型計算自由度曲線，雖然部分荷載步的計算自由度大于剛度矩陣維數(shù)，但大多數(shù)荷載步的計算自由度卻很小甚至模型仍處于彈性狀態(tài)。

圖15 壩頂橫向加速度反應(yīng)Fig. 15 Horizontal acceleration responses of dam crest

圖16 壩頂橫向加速度反應(yīng)Fig. 16 Vertical acceleration responses of dam crest

圖17 計算自由度曲線Fig. 17 Calculate DOFs curves

表4 給出了兩種比例邊界有限元非線性分析方法控制方程的求解時間，可以看出SBFEM 的計算時間遠大于近似IS-SBFEM 的計算時間，約為其7 倍。該多邊形模型其剛度矩陣維數(shù)為75406，其一次分解與回代的計算時間分別約為1.8 s 與0.03 s。由于近似IS-SBFEM 的一次迭代過程僅為多次常剛度矩陣回代以及稀疏矩陣與向量的乘積，對于大規(guī)模非線性問題，采用IS-SBFEM 計算優(yōu)勢明顯。由于存在高階的舒爾補矩陣，精確的IS-SBFEM 不再適用?？梢?，對于大規(guī)模問題，IS-SBFEM 的計算效率遠遠大于SBFEM，此時ISSBFEM 應(yīng)作為非線性分析的首選，既可在保證非線性計算精度的同時也具有極高的計算效率。

表4 控制方程計算時間Table 4 Calculate time of governing equation

5 結(jié)論

基于多邊形比例邊界有限元法基本理論，將隔離非線性思想用于多邊形比例邊界有限元非線性分析。在每個邊界線單元覆蓋的扇形區(qū)域內(nèi)設(shè)置若干個非線性應(yīng)變插值點，并以插值的形式建立其非線性應(yīng)變場，從而構(gòu)造了隔離非線性多邊形比例邊界單元。采用虛功原理建立了隔離非線性比例邊界單元的控制方程，單元的彈塑性矩陣可直接調(diào)用有限元程序的本構(gòu)模塊獲得，從而實現(xiàn)了多邊形比例邊界有限元高效非線性分析。得到如下結(jié)論：

(1)每個扇形區(qū)域其徑向自然坐標以及環(huán)向自然坐標分別為ξ=[0,1]以及η=[?1,1]，可采用標準高斯積分方案進行數(shù)值積分，且對于高階多邊形比例邊界單元也同樣適用。

(2)多邊形單元中引入較多的非線性應(yīng)變插值點使得舒爾補矩陣維數(shù)較大，采用Woodbury 近似法聯(lián)合算法對控制方程進行求解，可以保證隔離非線性比例邊界有限元法的計算精度與高效性。

(3)對于大規(guī)模問題，建議采用隔離非線性比例邊界有限元法進行分析，以便獲得更高的計算效率。將該方法進行推廣，對工程結(jié)構(gòu)的非線性分析具有重要意義。