李佳龍，李 鋼，余丁浩

(大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧，大連 116024)

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

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

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

1 SBFEM 基本理論

1.1 多邊形比例邊界有限元

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

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

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

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

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

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

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

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

1.2 非線性比例邊界有限元

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

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

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

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

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

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

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

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

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

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

基于彈塑性力學(xué)的基本理論，單元內(nèi)任意點(diǎn)處的應(yīng)變?cè)隽喀う?ξ,η)均可分解為線彈性的應(yīng)變?cè)?/p>

2.2 控制方程

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

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

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

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

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

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

表1 扇形單元積分點(diǎn)坐標(biāo)Table 1 Integral point coordinates of sector element

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

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

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

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

其中，子矩陣分別為：

3 控制方程求解

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

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

基于本文所提出的隔離非線性比例邊界有限元模型，編制了相應(yīng)的計(jì)算程序，整體分析流程如圖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)格以及線性四邊形有限元單元進(jìn)行網(wǎng)格剖分，其中多邊形網(wǎng)格共計(jì)109 個(gè)單元以及264 個(gè)結(jié)點(diǎn)，如圖5所示。有限元網(wǎng)格采用三種網(wǎng)格尺寸，分別為80 個(gè)單元與105 個(gè)結(jié)點(diǎn)、320 個(gè)單元與369 個(gè)結(jié)點(diǎn)以及1280 個(gè)單元與1377 個(gè)結(jié)點(diǎn)，分別如圖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é)點(diǎn)A處施加豎向荷載F=100 kN，并假定網(wǎng)格最密的有限元模型c 為參考解。表2給出了結(jié)點(diǎn)A在不同網(wǎng)格模型下橫向和豎向位移，可以看出采用本文提出的數(shù)值積分方案所得的計(jì)算結(jié)果與半解析的SBFEM 結(jié)果基本一致，且與有限元模型b 結(jié)果也基本一致。

表2 A 點(diǎn)處的位移Table 2 Displacement of node A

4.1.2 彈塑性分析

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

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

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

圖8 計(jì)算自由度曲線Fig. 8 Calculate DOFs curves

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

表3 控制方程計(jì)算時(shí)間Table 3 Calculate time of the governing equation

4.2 Koyna 重力壩動(dòng)力非線性分析

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

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

圖10 Koyna 地震動(dòng)記錄Fig. 10 Seismic acceleration records of Koyna

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

圖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 給出了該重力壩模型計(jì)算自由度曲線，雖然部分荷載步的計(jì)算自由度大于剛度矩陣維數(shù)，但大多數(shù)荷載步的計(jì)算自由度卻很小甚至模型仍處于彈性狀態(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 計(jì)算自由度曲線Fig. 17 Calculate DOFs curves

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

表4 控制方程計(jì)算時(shí)間Table 4 Calculate time of governing equation

5 結(jié)論

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

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

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

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