李 俊,汪基偉,冷 飛

(河海大學(xué)土木與交通學(xué)院,江蘇 南京 210098)

擴(kuò)展有限元法及其在鋼筋混凝土結(jié)構(gòu)中的應(yīng)用綜述

李 俊,汪基偉,冷 飛

(河海大學(xué)土木與交通學(xué)院,江蘇 南京 210098)

簡(jiǎn)述擴(kuò)展有限元法的基本思想和應(yīng)用價(jià)值,以鋼筋混凝土的裂紋問(wèn)題(強(qiáng)不連續(xù)問(wèn)題)為例對(duì)擴(kuò)展有限元法中描述裂紋的水平集法、不連續(xù)位移場(chǎng)富集策略、積分方案進(jìn)行了介紹,并給出相關(guān)領(lǐng)域近年來(lái)的研究成果。概述擴(kuò)展有限元法在斷裂力學(xué)領(lǐng)域中的應(yīng)用動(dòng)態(tài)。在此基礎(chǔ)上,分別從混凝土開裂和鋼筋與混凝土兩種材料的黏結(jié)滑移兩個(gè)方面介紹擴(kuò)展有限元法在鋼筋混凝土結(jié)構(gòu)開裂分析中的研究與應(yīng)用進(jìn)展,指出為將該方法應(yīng)用于水工鋼筋混凝土結(jié)構(gòu)等規(guī)模巨大、配筋和受力復(fù)雜的實(shí)際工程中,需對(duì)單元模型、裂縫處理方式以及裂縫處鋼筋應(yīng)力等問(wèn)題進(jìn)行深入研究。

擴(kuò)展有限元法;鋼筋混凝土結(jié)構(gòu);斷裂力學(xué);綜述

20世紀(jì)末期,Mo?s等[1]在標(biāo)準(zhǔn)有限元法的基礎(chǔ)上,結(jié)合單位分解法、水平集法和斷裂力學(xué)研究成果,提出了擴(kuò)展有限元法(extended finite element method, XFEM),為不連續(xù)問(wèn)題提供了新的解決方法。該方法通過(guò)在標(biāo)準(zhǔn)位移形函數(shù)中增加具有不連續(xù)性質(zhì)的富集函數(shù),在單元內(nèi)部描述位移間斷(強(qiáng)不連續(xù))、應(yīng)變間斷(弱不連續(xù)),使得有限元網(wǎng)格獨(dú)立于裂紋幾何路徑,在追蹤裂紋擴(kuò)展時(shí)亦無(wú)需重剖網(wǎng)格,提高了裂紋問(wèn)題的求解效率。

在水利水電工程中,鋼筋混凝土結(jié)構(gòu)應(yīng)用廣泛,其限裂分析對(duì)結(jié)構(gòu)的正常使用性能及耐久性有重要意義,重要的大體積非桿系混凝土結(jié)構(gòu)的限裂分析通常需要借助非線性有限元進(jìn)行數(shù)值模擬[2]。傳統(tǒng)鋼筋混凝土有限元中,彌散裂縫模型應(yīng)用廣泛[3],該模型不需隨開裂區(qū)擴(kuò)展而更新網(wǎng)格,適用于大型工程,但裂縫分布及裂縫寬度有一定網(wǎng)格依賴性,且無(wú)法準(zhǔn)確模擬裂縫處的鋼筋應(yīng)力。事實(shí)上,鋼筋混凝土結(jié)構(gòu)開裂分析時(shí),不僅要考慮混凝土基體的位移不連續(xù),還要考慮鋼筋的加強(qiáng)作用及其與混凝土之間的黏結(jié)滑移,較為復(fù)雜。將擴(kuò)展有限元法應(yīng)用于鋼筋混凝土結(jié)構(gòu)開裂分析,能調(diào)和傳統(tǒng)有限元計(jì)算效率和計(jì)算精度之間的矛盾,有較高的理論價(jià)值和廣闊的應(yīng)用前景。

1 擴(kuò)展有限元法基本原理

1.1 水平集法

水平集法(level set method, LSM)可以從幾何上描述靜態(tài)或動(dòng)態(tài)裂紋,以及孔洞、夾雜等不連續(xù)界面,它將界面的變化表示為比界面高一維(時(shí)間t)的水平集曲線[4],設(shè)不連續(xù)界面為Γ(t),任意考察點(diǎn)x∈Γ(t)滿足

(1)

式中:φ(x,t)為考察點(diǎn)水平集函數(shù)。

對(duì)于裂紋問(wèn)題,擴(kuò)展有限元中常用符號(hào)距離函數(shù)作為其水平集函數(shù),設(shè)xΓ為不連續(xù)界面上距x最近的點(diǎn),有

(2)

Sukumar等[5]介紹了如何在擴(kuò)展有限元中采用水平集法描述孔洞、夾雜等弱不連續(xù)界面。Ventura等[6]采用水平集法描述裂紋擴(kuò)展過(guò)程,提出一種改進(jìn)的向量水平集公式,使得在描述裂紋擴(kuò)展的同時(shí)不改變已有裂紋的水平集值,并且裂紋沿原裂紋面法向擴(kuò)展時(shí),新裂紋水平集值也不會(huì)改變。水平集法可以精確刻畫不連續(xù)面的幾何性質(zhì)及其動(dòng)態(tài)變化,易于從平面問(wèn)題推廣到三維問(wèn)題[7]。

1.2 裂紋問(wèn)題的擴(kuò)展有限元富集策略

Melenk等[8]研究了有限元形函數(shù)的單位分解法,單位分解法的基本思想是任意函數(shù)F(x)都可以在計(jì)算域內(nèi)表示為如下形式:

(3)

其中NI(x)滿足單位分解,即

(4)

式中:NI(x)為標(biāo)準(zhǔn)有限元形函數(shù);I為求解域節(jié)點(diǎn)集;f(x)為F(x)在計(jì)算域內(nèi)的局部近似函數(shù)。

單位分解法為改進(jìn)標(biāo)準(zhǔn)有限元形函數(shù)提供了理論依據(jù)?；趩挝环纸夥?可以根據(jù)不連續(xù)問(wèn)題的特點(diǎn),在標(biāo)準(zhǔn)有限元位移近似場(chǎng)的基礎(chǔ)上增加富集項(xiàng),用來(lái)模擬單元內(nèi)部存在的強(qiáng)、弱不連續(xù)性。

擴(kuò)展有限元中,改進(jìn)后的位移場(chǎng)uh(x)由標(biāo)準(zhǔn)有限元部分和富集項(xiàng)疊加得到[9]:

(5)

式中:usta(x)為標(biāo)準(zhǔn)有限元位移;uenr(x)為富集位移;ui,sta為標(biāo)準(zhǔn)有限元位移未知量;Ni(x)為標(biāo)準(zhǔn)有限元形函數(shù);Ni*(x)為單位分解函數(shù);Ψ(x)為富集函數(shù);ai為附加自由度;I*為富集節(jié)點(diǎn)集。Ni*(x)可以和Ni(x)相同,也可不同,為了編制程序方便,可取為相同。

不連續(xù)問(wèn)題包括強(qiáng)不連續(xù)問(wèn)題(如裂紋問(wèn)題中的位移不連續(xù))、弱不連續(xù)問(wèn)題(如夾雜問(wèn)題和復(fù)合材料問(wèn)題中的應(yīng)變不連續(xù))等,不同類型的問(wèn)題需要不同的富集位移構(gòu)建方法。以下以二維裂紋問(wèn)題為例,對(duì)擴(kuò)展有限元的位移富集思路進(jìn)行說(shuō)明。圖1為包含任意形狀裂紋的二維規(guī)則化網(wǎng)格,裂紋所在單元的相關(guān)節(jié)點(diǎn)為富集節(jié)點(diǎn),裂紋影響單元分為裂紋貫穿單元(圖1中2～6號(hào)單元)和裂尖單元(圖1中1、7號(hào)單元)兩部分,相應(yīng)富集節(jié)點(diǎn)亦分為兩部分，即裂尖單元富集節(jié)點(diǎn)和裂紋貫穿單元富集節(jié)點(diǎn),若某節(jié)點(diǎn)同時(shí)屬于裂紋貫穿單元和裂尖單元,為保證裂尖位移場(chǎng)計(jì)算精度,令其優(yōu)先屬于裂尖富集節(jié)點(diǎn)[1]。富集單元和標(biāo)準(zhǔn)單元之間由混合單元(圖1中陰影區(qū)單元)連接,混合單元中同時(shí)包含富集節(jié)點(diǎn)和標(biāo)準(zhǔn)節(jié)點(diǎn)(不參與富集的節(jié)點(diǎn))。

圖1 包含任意形狀裂紋的二維網(wǎng)格

對(duì)裂紋貫穿單元,采用階躍函數(shù)構(gòu)建富集函數(shù),以反映裂紋兩側(cè)的位移不連續(xù);對(duì)裂尖單元,采用斷裂力學(xué)中裂尖位移場(chǎng)解析解的各項(xiàng)作為基函數(shù)構(gòu)建富集函數(shù)[10]。參照式(5),二維裂紋問(wèn)題的位移模式可統(tǒng)一寫為

(7)

(8)

式中:NΓ為裂紋貫穿單元富集節(jié)點(diǎn)集;NΛ為裂尖單元富集節(jié)點(diǎn)集;aj和bw,k分別為上述兩類富集節(jié)點(diǎn)的附加自由度;Nj(x)、Nk(x)為單位分解函數(shù);H(x)為階躍函數(shù);Fw(x)為裂尖函數(shù);r、θ為裂尖極坐標(biāo)系所定義的位置參數(shù)。

Belytschko等[11]對(duì)裂紋貫穿單元富集函數(shù)做出修正,在斷裂問(wèn)題中應(yīng)用較多:

(9)

將Ψi(x)代替式(5)中Ψ(x),并令x=xj(xj為混合單元節(jié)點(diǎn)坐標(biāo)),有uh(xj)=uj,sta,式中uj,sta為j節(jié)點(diǎn)的標(biāo)準(zhǔn)有限元位移未知量?？梢钥吹?改進(jìn)后得到的節(jié)點(diǎn)位移為標(biāo)準(zhǔn)有限元形函數(shù)的系數(shù),這一點(diǎn)和標(biāo)準(zhǔn)有限元一致,便于位移后處理,其原因在于對(duì)富集函數(shù)進(jìn)行修正后,混合單元內(nèi)位移函數(shù)富集項(xiàng)取值為零,即混合單元中富集節(jié)點(diǎn)的附加自由度不對(duì)混合單元產(chǎn)生影響,只對(duì)裂紋貫穿單元內(nèi)部產(chǎn)生影響。對(duì)于裂尖置于單元內(nèi)部的情況,可在富集函數(shù)中引入斜坡函數(shù)R(x)[12],以提高混合單元收斂速度,此時(shí)用式(10)中的Ψi(x)代替式(5)中的Ψ(x):

(10)

1.3 擴(kuò)展有限元的積分方案

擴(kuò)展有限元中形函數(shù)具有單位分解屬性,使得剛度矩陣保持了對(duì)稱、稀疏的特性,這一點(diǎn)和傳統(tǒng)有限元一致。但對(duì)于不連續(xù)單元(如裂尖單元和裂紋貫穿單元),在獲得剛度矩陣時(shí),其積分不能依照傳統(tǒng)連續(xù)單元進(jìn)行。不連續(xù)單元的積分方案中,分區(qū)域積分較為常用[13]。

對(duì)于裂紋貫穿單元,由于裂紋兩側(cè)子區(qū)域的位移場(chǎng)是連續(xù)的,所以可以在子區(qū)域分別積分,然后將子區(qū)域的積分疊加得到單元?jiǎng)偠染仃?。?duì)于裂尖單元,當(dāng)以裂紋尖端為頂點(diǎn)分成若干小三角形進(jìn)行分區(qū)域積分時(shí),裂尖處會(huì)出現(xiàn)奇異性,Laborde等[14]介紹了通過(guò)坐標(biāo)變換消除奇異性的方法。

由于分區(qū)域積分方案計(jì)算量較大,也有學(xué)者提出采用連續(xù)積分方案作為替代。Ventura[15]針對(duì)簡(jiǎn)單富集函數(shù)提出一種可將分區(qū)域積分轉(zhuǎn)化為連續(xù)高斯積分的方法,但當(dāng)單元形態(tài)不是平行四邊形時(shí)有明顯誤差。Song等[16]對(duì)擴(kuò)展有限元中位移基函數(shù)和自由度進(jìn)行重置,采用疊合單元和虛擬節(jié)點(diǎn)表示裂縫,在此基礎(chǔ)上提出簡(jiǎn)易積分方案,該方案每個(gè)單元設(shè)置一個(gè)積分點(diǎn),為保證矩陣非奇異,需采取控制措施避免零能模式。喻葭臨等[17]結(jié)合文獻(xiàn)[16]的思路和傳統(tǒng)雙線性四邊形單元的四節(jié)點(diǎn)積分方案,提出同時(shí)適用于含裂縫單元與標(biāo)準(zhǔn)單元(不含裂縫)的簡(jiǎn)易積分方案,避免了分區(qū)域積分的復(fù)雜計(jì)算,亦無(wú)需對(duì)零能模式采取控制措施,有較高的實(shí)用價(jià)值。

2 擴(kuò)展有限元法應(yīng)用概況

擴(kuò)展有限元允許在單元內(nèi)部描述裂縫,裂縫擴(kuò)展時(shí)無(wú)需重剖網(wǎng)格,同時(shí)保留了傳統(tǒng)有限元求解非線性問(wèn)題的優(yōu)勢(shì)[18],被廣泛應(yīng)用于斷裂力學(xué)領(lǐng)域。

Sukumar等[19-20]針對(duì)準(zhǔn)靜態(tài)裂紋擴(kuò)展問(wèn)題對(duì)擴(kuò)展有限元數(shù)值實(shí)現(xiàn)步驟進(jìn)行說(shuō)明。董玉文等[21]采用擴(kuò)展有限元法計(jì)算應(yīng)力強(qiáng)度因子,指出其精度和采用J積分計(jì)算得到的應(yīng)力強(qiáng)度因子精度相當(dāng)。除靜態(tài)和準(zhǔn)靜態(tài)裂紋問(wèn)題外,擴(kuò)展有限元法亦被應(yīng)用于動(dòng)態(tài)裂紋問(wèn)題中[22-23]。Belytschko等[22]采用漸進(jìn)的動(dòng)態(tài)裂紋尖端解作為基函數(shù),把I型和II型應(yīng)力強(qiáng)度因子作為附加自由度,使用顯示積分方案求解彈性動(dòng)態(tài)裂紋擴(kuò)展問(wèn)題。Menouillard等[23]對(duì)顯示積分系統(tǒng)的穩(wěn)定性進(jìn)行研究,并給出質(zhì)量矩陣的對(duì)角化方法。擴(kuò)展有限元的優(yōu)勢(shì)之一在于不預(yù)設(shè)開裂路徑,盡可能地還原裂縫的真實(shí)形態(tài),這一點(diǎn)在分支裂紋、交叉裂紋[24-25]的模擬中得到展現(xiàn)。

真實(shí)裂紋之間可能存在摩擦接觸,在水工結(jié)構(gòu)中,裂紋中還可能存在滲透壓力作用,是否考慮裂縫內(nèi)的作用對(duì)裂縫預(yù)測(cè)有較大影響。Dolbow等[26]采用擴(kuò)展有限元法研究裂紋面間的摩擦接觸,采用非線性本構(gòu)關(guān)系來(lái)描述裂紋面的接觸,并采用LATIN迭代法求解;余天堂[27]采用線性互補(bǔ)法求解裂紋面的非線性接觸,能夠避免復(fù)雜的迭代求解過(guò)程;李建波等[28-29]采用虛功原理推導(dǎo)了考慮單元內(nèi)部裂紋面上分布荷載及縫內(nèi)粘連的擴(kuò)展有限元基本公式,提出和標(biāo)準(zhǔn)有限元協(xié)調(diào)一致的擴(kuò)展有限元?jiǎng)偠染仃囆纬赡Ｊ健?/p>

紙產(chǎn)品主要通過(guò)造紙機(jī)來(lái)生產(chǎn)制造，生產(chǎn)過(guò)程包括漿料準(zhǔn)備、成形、壓榨和干燥等主要工序。干燥因其體積最大、固定資產(chǎn)投資最大、能源消耗最高，被認(rèn)為是最關(guān)鍵的工序。因此，本研究選擇紙張干燥過(guò)程作為切入點(diǎn)，總結(jié)紙張干燥過(guò)程建模與智能模擬技術(shù)的研究進(jìn)展，為建立紙張干燥過(guò)程系統(tǒng)模型，實(shí)現(xiàn)智能模擬紙張干燥生產(chǎn)過(guò)程并最終助力傳統(tǒng)造紙工業(yè)轉(zhuǎn)型升級(jí)、實(shí)現(xiàn)智能制造，積累技術(shù)力量。

三維斷裂問(wèn)題一直是研究的熱點(diǎn)和難點(diǎn),Sukumar等[30]首次利用擴(kuò)展有限元法研究了三維裂紋擴(kuò)展問(wèn)題,其裂尖富集函數(shù)建立在與裂紋尖端垂直的平面內(nèi),和平面問(wèn)題類似,其函數(shù)形式仍采用極坐標(biāo)表示。三維裂紋擴(kuò)展問(wèn)題的難點(diǎn)在于裂紋面的連續(xù)性與光滑性較難模擬,Duan等[31]利用單元水平集描述三維裂縫,采用最小二乘法改進(jìn)了原裂縫面方向和新預(yù)測(cè)開裂方向的一致性,從而改善了三維裂縫面的光滑性。

除斷裂力學(xué)領(lǐng)域之外,擴(kuò)展有限元在復(fù)合材料界面失效(同時(shí)包含強(qiáng)、弱不連續(xù))[32]以及變形局部化(如剪切帶、損傷過(guò)程區(qū))[16,33]等領(lǐng)域也有應(yīng)用。值得注意的是,擴(kuò)展有限元研究集中于線彈性斷裂力學(xué)領(lǐng)域,當(dāng)材料屬性發(fā)生變化時(shí),其斷裂性能也會(huì)發(fā)生變化,相關(guān)參數(shù)如斷裂能和富集函數(shù)構(gòu)造形式也會(huì)發(fā)生變化,應(yīng)根據(jù)所求問(wèn)題具體分析。

3 擴(kuò)展有限元法在鋼筋混凝土結(jié)構(gòu)開裂分析中的應(yīng)用

鋼筋混凝土結(jié)構(gòu)開裂分析時(shí)存在兩種不連續(xù)問(wèn)題:一是混凝土開裂導(dǎo)致的位移不連續(xù);二是鋼筋與混凝土兩種材料的黏結(jié)滑移。擴(kuò)展有限元法在上述兩種不連續(xù)模擬中均有應(yīng)用,現(xiàn)簡(jiǎn)述如下。

3.1 在混凝土開裂分析中的應(yīng)用

混凝土斷裂力學(xué)認(rèn)為在裂縫尖端區(qū)域存在斷裂過(guò)程區(qū),斷裂過(guò)程區(qū)是微細(xì)裂縫發(fā)生、發(fā)展并逐步轉(zhuǎn)化為宏觀裂紋的區(qū)域。黏聚裂紋模型(又稱凝聚力模型)[34]采用宏觀裂紋尖端的黏聚區(qū)表示斷裂過(guò)程區(qū),黏聚區(qū)的傳力機(jī)理通過(guò)應(yīng)力-張開位移曲線進(jìn)行描述,該曲線與坐標(biāo)軸所圍面積代表相應(yīng)材料的斷裂能。

黏聚裂紋模型可以較為真實(shí)地模擬裂尖受力狀態(tài),但需預(yù)先設(shè)定起裂位置,杜效鵲等[39]結(jié)合擴(kuò)展有限元和黏聚裂紋模型,研究了預(yù)制縫重力壩模型的斷裂特性,得出與試驗(yàn)結(jié)果一致的荷載響應(yīng)曲線和裂縫擴(kuò)展路徑;張曉東等[40]利用黏聚裂紋模型研究了帶初始邊緣裂紋的混凝土板在單向拉伸作用下的裂紋擴(kuò)展過(guò)程,效果良好。然而實(shí)際工程中并不是總能預(yù)先獲得開裂位置,而且一般有多條裂紋發(fā)生,此時(shí)過(guò)分關(guān)注裂尖穩(wěn)定是不必要的。陳勝宏等[41]采用擴(kuò)展有限元法對(duì)小灣拱壩壩踵開裂進(jìn)行分析,采用最大拉應(yīng)力準(zhǔn)則判斷裂縫是否擴(kuò)展,并且認(rèn)為裂縫發(fā)生后即貫穿整個(gè)單元,該處理方法便于處理大型復(fù)雜結(jié)構(gòu),但計(jì)算結(jié)果表明裂縫擴(kuò)展范圍與網(wǎng)格尺寸有一定關(guān)系。

在細(xì)觀層次上,混凝土被認(rèn)為是由骨料、硬化砂漿和二者之間的過(guò)渡區(qū)組成的復(fù)合材料,其中骨料和硬化砂漿交界處易產(chǎn)生微細(xì)裂紋,微細(xì)裂紋的發(fā)展和貫通導(dǎo)致宏觀裂縫的出現(xiàn)。采用擴(kuò)展有限元法模擬混凝土細(xì)觀斷裂過(guò)程需要解決以下問(wèn)題:①不同材料間的弱不連續(xù)模擬,即復(fù)合材料界面(過(guò)渡區(qū))的模擬;②微細(xì)裂紋產(chǎn)生、發(fā)展導(dǎo)致的強(qiáng)不連續(xù)模擬;③復(fù)合材料界面附近斷裂參數(shù)的確定。江守燕等[42]采用擴(kuò)展有限元法對(duì)包含圓形隨機(jī)骨料的多夾雜問(wèn)題進(jìn)行弱不連續(xù)分析,用水平集法描述材料界面,采用標(biāo)準(zhǔn)網(wǎng)格進(jìn)行離散,減輕了網(wǎng)格剖分的負(fù)擔(dān),但該研究假定材料為線彈性,不涉及微細(xì)裂紋的產(chǎn)生和擴(kuò)展。杜修力等[43]采用擴(kuò)展有限元法模擬了混凝土微細(xì)裂紋的產(chǎn)生和發(fā)展過(guò)程,將骨料周圍2mm厚度范圍內(nèi)設(shè)置為過(guò)渡區(qū),給出骨料、過(guò)渡區(qū)、砂漿的力學(xué)參數(shù)(開裂前假定三者都是線彈性材料,開裂后采用線性的應(yīng)力-裂紋寬度關(guān)系表征混凝土軟化特性)。由于假定裂紋產(chǎn)生后隨即貫穿整個(gè)單元,為保證收斂性和求解精度,該研究網(wǎng)格剖分較為精細(xì),對(duì)過(guò)渡區(qū)需單獨(dú)考慮,計(jì)算成本較大。于紅軍[44]研究了含復(fù)雜界面的非均勻材料的斷裂參數(shù)確定方法,采用相互作用積分計(jì)算位于界面附近或界面上的裂尖應(yīng)力強(qiáng)度因子,驗(yàn)證了相互作用積分的穩(wěn)定性和區(qū)域無(wú)關(guān)性。

上述研究將擴(kuò)展有限元法應(yīng)用到素混凝土的開裂模擬中,可為擴(kuò)展有限元法應(yīng)用于鋼筋混凝土結(jié)構(gòu)提供參考,但由于沒(méi)有考慮鋼筋在限制混凝土開裂中起到的作用,因而無(wú)法直接應(yīng)用于鋼筋混凝土結(jié)構(gòu)開裂分析。

3.2 在鋼筋與混凝土界面模擬中的應(yīng)用

不同于素混凝土結(jié)構(gòu),鋼筋混凝土結(jié)構(gòu)需要考慮鋼筋的加強(qiáng)作用以及鋼筋與混凝土兩種材料的相互作用,其中鋼筋與混凝土的黏結(jié)滑移是關(guān)鍵問(wèn)題,是否考慮黏結(jié)滑移對(duì)裂縫分布及裂縫寬度有明顯影響。茹忠亮等[45]用擴(kuò)展有限元法對(duì)預(yù)設(shè)初始裂縫的鋼筋混凝土梁復(fù)合斷裂過(guò)程進(jìn)行模擬,假定鋼筋與混凝土之間沒(méi)有滑移,得到的極限承載力較試驗(yàn)值偏大,鋼筋布置區(qū)的裂縫分布和試驗(yàn)情況也有較大差別。傳統(tǒng)鋼筋混凝土有限元模型依據(jù)鋼筋模擬方法的不同分為3種:分離式、組合式、整體式,其中分離式和組合式模型可以模擬鋼筋與混凝土的黏結(jié)滑移,在裂縫分析中應(yīng)用較多[46],也可在擴(kuò)展有限元分析中借鑒采用。

Simone等[47]把擴(kuò)展有限元的思想引入鋼筋與混凝土界面的模擬,利用有限元形函數(shù)的單位分解性質(zhì),對(duì)鋼筋與混凝土界面上的節(jié)點(diǎn)預(yù)設(shè)附加自由度以考慮二者的相對(duì)滑移。Deb等[48]在研究灌漿錨桿時(shí),定義含錨桿的實(shí)體單元為富集單元,富集單元有附加自由度,用于確定錨桿的位移、應(yīng)力,富集單元的剛度矩陣由巖體、砂漿、錨桿共同組成。Simone等[47]的研究基于分離式模型,Deb等[48]的研究基于組合式模型,但二者在分析中沒(méi)有考慮混凝土開裂或巖體開裂導(dǎo)致的不連續(xù)位移。Ibrahimbegovic等[49]基于平面組合式模型,在包含鋼筋的混凝土單元位移場(chǎng)中加入黏結(jié)滑移附加項(xiàng),基體開裂則采用內(nèi)嵌不連續(xù)模型進(jìn)行模擬,然后把位移分解為不計(jì)滑移時(shí)鋼筋混凝土的位移與黏結(jié)滑移引起的位移之和,該研究?jī)?yōu)點(diǎn)在于同時(shí)考慮了基體開裂和黏結(jié)滑移兩種不連續(xù)。何學(xué)[50]建立了分離式鋼筋混凝土擴(kuò)展有限元模型,采用四節(jié)點(diǎn)擴(kuò)展有限元等參單元模擬混凝土基體,采用桿單元模擬鋼筋,并對(duì)傳統(tǒng)四邊形黏結(jié)單元進(jìn)行改進(jìn),與鋼筋單元相連的部分保持位移連續(xù),與混凝土相連的部分添加由裂縫產(chǎn)生的附加位移,使之可以考慮裂縫引起的混凝土位移間斷。文獻(xiàn)[49-50]保留了傳統(tǒng)鋼筋混凝土有限元的架構(gòu),對(duì)混凝土、鋼筋及二者的黏結(jié)滑移分別考慮,適用于不同配筋及復(fù)雜受力狀況,為擴(kuò)展有限元法應(yīng)用于實(shí)際鋼筋混凝土工程作出了有益探索。

非桿系鋼筋混凝土結(jié)構(gòu)常通過(guò)控制鋼筋應(yīng)力間接控制裂縫寬度[2],因而準(zhǔn)確得到裂縫處的鋼筋應(yīng)力有重要意義。Contrafatto等[51]研究了規(guī)則配筋軸心受拉構(gòu)件的鋼筋應(yīng)力分布,用變分原理推導(dǎo)了考慮基體開裂、黏結(jié)滑移的基本方程,探討了混凝土形函數(shù)和鋼筋形函數(shù)的改進(jìn)方法,從而能夠描述由裂縫、滑移引起的局部大梯度應(yīng)力(應(yīng)變)場(chǎng),該研究對(duì)準(zhǔn)確求得復(fù)雜配筋及復(fù)雜受力狀態(tài)下裂縫處鋼筋應(yīng)力有參考意義。

另外,通用有限元軟件ABAQUS在原有的斷裂力學(xué)分析功能基礎(chǔ)上增加了擴(kuò)展有限元的分析模塊,許多學(xué)者基于此進(jìn)行了素混凝土構(gòu)件和鋼筋混凝土構(gòu)件的斷裂分析?？梢猿浞掷猛ㄓ糜邢拊绦虻姆蔷€性求解功能[52],也可通過(guò)加入彈簧單元的方法考慮鋼筋與混凝土的黏結(jié)滑移[53],目前該軟件的擴(kuò)展有限元分析功能仍處在不斷發(fā)展和豐富中。

3.3 需要繼續(xù)研究的問(wèn)題

由于鋼筋混凝土結(jié)構(gòu)開裂機(jī)理的復(fù)雜性,多數(shù)研究側(cè)重于兩種不連續(xù)問(wèn)題中的一個(gè)方面,少數(shù)研究對(duì)兩種不連續(xù)問(wèn)題進(jìn)行統(tǒng)一考慮[49-51,53]。擴(kuò)展有限元的思想可以和傳統(tǒng)鋼筋混凝土有限元結(jié)合(文獻(xiàn)[50]稱之為鋼筋混凝土擴(kuò)展有限元),為同時(shí)解決鋼筋混凝土結(jié)構(gòu)中的兩類不連續(xù)問(wèn)題提供思路。隨著結(jié)構(gòu)設(shè)計(jì)中正常使用性能及耐久性的重要性逐漸增大,需要更高效的裂縫分布計(jì)算方法,既能保證一定的計(jì)算精度,其計(jì)算成本又不是太高,擴(kuò)展有限元在該領(lǐng)域有廣闊的應(yīng)用前景。但是,目前鋼筋混凝土擴(kuò)展有限元還處于研究階段,尚有許多問(wèn)題需要深入研究:

a. 構(gòu)建合適的單元模型,能夠模擬鋼筋與混凝土的黏結(jié)滑移,同時(shí)滿足實(shí)際工程中復(fù)雜配筋的需要。傳統(tǒng)鋼筋混凝土有限元中,組合式模型可以考慮任意方向、任意數(shù)量的配筋,通過(guò)界面單元模擬鋼筋與混凝土的黏結(jié)滑移,可借鑒文獻(xiàn)[50]的方法對(duì)傳統(tǒng)界面單元進(jìn)行改進(jìn),在界面單元與混凝土相連處增加附加位移項(xiàng),使之能夠考慮混凝土開裂帶來(lái)的影響。

b. 構(gòu)建合適的裂縫模型,既能滿足裂縫分布計(jì)算精度要求,又不過(guò)分復(fù)雜,符合工程需求。斷裂力學(xué)裂縫模型需要較多參數(shù),在分析、設(shè)計(jì)時(shí)存在較多不確定性,而且鋼筋混凝土結(jié)構(gòu)中由于鋼筋參與工作,混凝土裂縫的開展受到限制,裂縫不會(huì)因裂尖處應(yīng)力集中導(dǎo)致不穩(wěn)定擴(kuò)展,過(guò)分關(guān)注裂尖應(yīng)力狀態(tài)是不必要的,因而傳統(tǒng)的強(qiáng)度理論(如最大拉應(yīng)力準(zhǔn)則)在裂縫分析、限裂設(shè)計(jì)中仍有應(yīng)用價(jià)值。此外,通常大體積混凝土結(jié)構(gòu)裂縫開展深度較大,在較小的網(wǎng)格尺寸下,仍可假定裂縫產(chǎn)生后隨即貫穿整個(gè)單元。

c. 多裂紋問(wèn)題及復(fù)合型開裂問(wèn)題。當(dāng)采用斷裂力學(xué)裂縫模型時(shí),單一裂紋的I型開裂已有較多研究,但多裂紋、交叉裂紋以及復(fù)合型開裂問(wèn)題研究較少。多裂紋以及復(fù)合型開裂問(wèn)題需要建立相應(yīng)的位移模式、裂紋擴(kuò)展準(zhǔn)則和裂紋匯合準(zhǔn)則,隨著裂紋數(shù)目增多,裂紋間相互作用更趨復(fù)雜。鋼筋混凝土結(jié)構(gòu)一般會(huì)產(chǎn)生較多裂縫,因而要反映其真實(shí)開裂情況仍有一定難度。如采用傳統(tǒng)強(qiáng)度理論作為開裂準(zhǔn)則,可簡(jiǎn)化分析,但仍需對(duì)位移模式以及裂縫間作用進(jìn)行研究。

d. 基體開裂情況下鋼筋與混凝土界面滑移的高精度模擬以及裂縫處鋼筋應(yīng)力的高精度算法。裂縫處鋼筋應(yīng)力和裂縫寬度具有正相關(guān)性,同時(shí)鋼筋應(yīng)力也是結(jié)構(gòu)(構(gòu)件)承載力的重要指標(biāo)。鋼筋一般可采用桿單元進(jìn)行模擬,桿單元的平均應(yīng)力即代表該單元內(nèi)鋼筋應(yīng)力,而事實(shí)上由于黏結(jié)切應(yīng)力的存在,單元內(nèi)鋼筋應(yīng)力分布是不均勻的,當(dāng)鋼筋橫穿裂縫時(shí),裂縫處鋼筋應(yīng)力有突躍現(xiàn)象?？煽紤]利用擴(kuò)展有限元位移富集思想構(gòu)建能夠描述裂縫處局部變形集中的鋼筋形函數(shù),非均勻有理B樣條曲線有良好的局部修改能力,且有良好的可導(dǎo)性,可用來(lái)準(zhǔn)確模擬裂縫處鋼筋應(yīng)力突躍。

4 結(jié) 論

a. 擴(kuò)展有限元法充分利用有限元形函數(shù)的單位分解屬性,使得不連續(xù)界面的描述獨(dú)立于網(wǎng)格剖分,在追蹤不連續(xù)界面擴(kuò)展時(shí)亦無(wú)需更新網(wǎng)格,同時(shí)由于其保留了標(biāo)準(zhǔn)有限元求解非線性問(wèn)題的優(yōu)勢(shì),因而在多種材料的強(qiáng)、弱不連續(xù)問(wèn)題研究中應(yīng)用廣泛。

b. 將擴(kuò)展有限元法應(yīng)用于鋼筋混凝土結(jié)構(gòu)開裂分析,有理論價(jià)值和應(yīng)用前景。同時(shí)考慮混凝土開裂和黏結(jié)滑移兩種不連續(xù)問(wèn)題的研究較少,鋼筋混凝土擴(kuò)展有限元可為此提供新的思路,但該方法尚處于研究階段,其理論體系有待于進(jìn)一步完善。

[ 1 ] MOЁS N,DOLBOW J,BELYTSCHKO T.A finite element method for crack growth without remeshing[J].International Journal for Numerical Methods in Engineering,1999,46(1): 131-150.

[ 2 ] SL 91—2008 水工混凝土結(jié)構(gòu)設(shè)計(jì)規(guī)范[S].

[ 3 ] 江見鯨,陸新征.混凝土結(jié)構(gòu)有限元分析[M].2版.北京: 清華大學(xué)出版社,2013.

[ 4 ] OSHER S,SETHIAN J A.Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[J].Journal of Computational Fhysics,1988,79(1): 12-49.

[ 5 ] SUKUMAR N,CHOPP D L,MOЁS N,et al.Modeling holes and inclusions by level sets in the extended finite-element method[J].Computer Methods in Applied Mechanics and Engineering,2001,190(46): 6183-6200.

[ 6 ] VENTURA G,BUDYN E,BELYTSCHKO T.Vector level sets for description of propagating cracks in finite elements [J].International Journal for Numerical Methods in Engineering,2003,58: 1571-1592.

[ 7 ] STOLARSKA M,CHOPP D L,MOЁS N,et al.Modelling crack growth by level sets in the extended finite element method[J].International Journal for Numerical Methods in Engineering,2001,51(8): 943-960.

[ 8 ] MELENK J M,BABUKA I.The partition of unity finite element method: basic theory and applications[J].Computer Methods in Applied Mechanics and Engineering,1996,139(1/2/3/4): 289-314.

[ 9 ] CHESSA J,WANG H,BELYTSCHKO T.On the construction of blending elements for local partition of unity enriched finite elements[J].International Journal for Numerical Methods in Engineering,2003,57(7): 1015-1038.

[10] BELYTSCHKO T,BLACK T.Elastic crack growth in finite elements with minimal remeshing[J].International Journal for Numerical Methods in Engineering,1999,45(5): 601-620.

[11] BELYTSCHKO T,CHEN H,XU J,et al.Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment [J].International Journal for Numerical Methods in Engineering,2003,58(12): 1873-1905.

[12] FRIES T P.A corrected XFEM approximation without problems in blending elements[J].International Journal for Numerical Methods in Engineering,2008,75(5): 503-532.

[13] 余天堂.擴(kuò)展有限元法的數(shù)值方面[J].巖土力學(xué),2007,28(增刊1): 305-310.(YU Tiantang.Numerical aspects of the extended finite element method[J].Rock and Solid Mechanics,2007,28 (Sup1): 305-310.(in Chinese))

[14] LABORDE P,POMMIER J,RENARD Y,et al.High-order extended finite element method for cracked domains[J].International Journal for Numerical Methods in Engineering,2005,64(3): 354-381.

[15] VENTURA G.On the elimination of quadrature subcells for description of propagating cracks in finite elements [J].International Journal for Numerical Methods in Engineering,2006,66: 761-795.

[16] SONG J H,AREIAS P,BELYTSCHKO T.A method for dynamic crack and shear band propagation with phantom nodes[J].International Journal for Numerical Methods in Engineering,2006,67(6): 868-893.

[17] 喻葭臨,張其光,于玉貞,等.擴(kuò)展有限元中非連續(xù)區(qū)域的一種積分方案[J].清華大學(xué)學(xué)報(bào):自然科學(xué)版,2009 49(3): 351-354.(YU Jialin,ZHANG Qiguang,YU Yuzhen,et al.An integration scheme for discontinuities in the extended finite element method[J].Journal of Tsinghua University: Science and Technology,2009,49(3): 351-354.(in Chinese))

[18] KARIHALOO B L,XIAO Q Z.Modelling of stationary and growing cracks in FE framework without remeshing: a state-of-the-art review[J].Computers & Structures,2003,81(3): 119-129.

[19] SUKUMAR N,PRVOST J H.Modeling quasi-static crack growth with the extended finite element method Part I: Computer implementation[J].International Journal of Solids and Structures,2003,40(26): 7513-7537.

[20] HUANG R,SUKUMAR N,PRVOST J H.Modeling quasi-static crack growth with the extended finite element method:part Ⅱ numerical applications[J].International Journal of Solids and Structures,2003,40(26): 7539-7552.

[21] 董玉文,余天堂,任青文.直接計(jì)算應(yīng)力強(qiáng)度因子的擴(kuò)展有限元法[J].計(jì)算力學(xué)學(xué)報(bào),2008,25(1): 72-77.(DONG Yuwen,YU Tiantang,REN Qingwen.Extended finite element method for direct evaluation of strength intensity factors[J].Chinese Journal of Computational Mechanics,2008,25(1): 72-77.(in Chinese))

[22] BELYTSCHKO T,CHEN H.Singular enrichment finite element method for elastodynamic crack propagation[J].International Journal of Computational Methods,2004,1(1): 1-15.

[23] MENOUILLARD T,RTHORJ,COMBESCURE A,et al.Efficient explicit time stepping for the eXtended Finite Element Method (X-FEM)[J].International Journal for Numerical Methods in Engineering,2006,68(9): 911-939.

[24] ZI G,SONG J H,BUDYN E,et al.A method for growing multiple cracks without remeshing and its application to fatigue crack growth [J].Modeling and Simulations for Material Science and Engineering,2004,12: 901-915.

[25] 石路楊,余天堂.多裂紋擴(kuò)展的擴(kuò)展有限元法分析[J].巖土力學(xué),2014,35(1): 263-272.(SHI Luyang,YU Tiantang.Analysis of multiple crack growth using extended finite element method[J].Rock and Soil Mechanics,2014,35(1): 263-272.(in Chinese))

[26] DOLBOW J,MOЁS N,BELYTSCHKO T.An extended finite element method for modeling crack growth with frictional contact[J].Computer Methods in Applied Mechanics and Engineering,2001,190(51): 6825-6846.

[27] 余天堂.摩擦接觸裂紋問(wèn)題的擴(kuò)展有限元法[J].工程力學(xué),2010 (4): 84-89.(YU Tiantang.An extended finite element method for modeling crack problems with frictional contact[J].Engineering Mechanics,2010 (4): 84-89.(in Chinese))

[28] 李建波,陳健云,林皋.非網(wǎng)格重剖分模擬宏觀裂紋體的擴(kuò)展有限單元法(1:基礎(chǔ)理論)[J].計(jì)算力學(xué)學(xué)報(bào),2006,23(2): 207-213.(LI Jianbo,CHEN Jianyun,LIN Gao.Extended finite element method for modeling cracks without remeshing (I: basic theory)[J].Chinese Journal of Computational Mechanics,2006,23(2): 207-213.(in Chinese))

[29] 李建波,陳健云,林皋.非網(wǎng)格重剖分模擬宏觀裂紋體的擴(kuò)展有限單元法(2: 數(shù)值實(shí)現(xiàn))[J].計(jì)算力學(xué)學(xué)報(bào),2006,23(3): 317-323.(LI Jianbo,CHEN Jianyun,LIN Gao.Extended finite element method for modeling cracks without remeshing (2: numerical implementation)[J].Chinese Journal of Computational Mechanics,2006,23(3): 317-323.(in Chinese))

[30] SUKUMAR N,MOЁS N,MORAN B,et al.Extended finite element method for three-dimensional crack modelling[J].International Journal for Numerical Methods in Engineering,2000,48(11): 1549-1570.

[31] DUAN Q,SONG J H,MENOUILLARD T,et al.Element-local level set method for three-dimensional dynamic crack growth[J].International Journal for Numerical Methods in Engineering,2009,80(12): 1520-1543.

[32] HETTICH T,RAMM E.Interface material failure modeled by the extended finite-element method and level sets [J].Computer Methods in Applied Mechanics and Engineering,2006,195: 4753-4767.

[34] ELICES M,GUINEA G V,GOMEZ J,et al.The cohesive zone model: advantages,limitations and challenges[J].Engineering Fracture Mechanics,2002,69(2): 137-163.

[35] MOЁS N,BELYTSCHKO T.Extended finite element method for cohesive crack growth[J].Engineering Fracture Mechanics,2002,69(7): 813-833.

[36] MARIANI S,PEREGO U.Extended finite element method for quasi-brittle fracture [J].International Journal for Numerical Methods in Engineering,2003,58: 103-126.

[37] UNGER J F,ECKARDT S,K?NKE C.Modelling of cohesive crack growth in concrete structures with the extended finite element method[J].Computer Methods in Applied Mechanics and Engineering,2007,196(41): 4087-4100.

[38] 方修君,金峰,王進(jìn)廷.用擴(kuò)展有限元方法模擬混凝土的復(fù)合型開裂過(guò)程[J].工程力學(xué),2007,24(增刊1): 46-52.(FANG Xiujun,JIN Feng,WANG Jinting.Simulation of mixed-mode fracture of concrete using extended finite element method[J].Engineering Mechanics,2007,24(Sup1): 46-52.(in Chinese))

[39] 杜效鵲,段云嶺,王光綸.重力壩斷裂數(shù)值分析研究[J].水利學(xué)報(bào),2005,36(9): 1-10.(DU Xiaohu,DUAN Yunling,WANG Guanglun.Numerical analysis of fracture in gravity dam[J].Journal of Hydraulic Engineering,2005,36(9): 1-10.(in Chinese))

[40] 張曉東,丁勇,任旭春.混凝土裂紋擴(kuò)展過(guò)程模擬的擴(kuò)展有限元法研究[J].工程力學(xué),2013,30(7): 14-21,27.(ZHANG Xiaodong,DING Yong,REN Xuchun.Simulation of the concrete crack propagation process with the extended finite element method[J].Engineering Mechanics,2013,30(7): 14-21,27.(in Chinese))

[41] 陳勝宏,汪衛(wèi)明,徐明毅,等.小灣高拱壩壩踵開裂的有限單元法分析[J].水利學(xué)報(bào),2003,34(1): 66-71.(CHEN Shenghong,WANG Weiming,XU Mingyi,et al.Finite element analysis of the crack propagation in high arch dam heel of Xiaowan project[J].Journal of Hydraulic Engineering,2003,34(1): 66-71.(in Chinese))

[42] 江守燕,杜成斌.弱不連續(xù)問(wèn)題擴(kuò)展有限元法的數(shù)值精度研究[J].力學(xué)學(xué)報(bào),2012,44(6): 1005-1015.(JIANG Shouyan,DU Chengbin.Study on numerical precision of extended finite element methods for modeling weak discontinuities[J].Chinese Journal of Theoretical and Applied Mechanics,2012,44(6): 1005-1015.(in Chinese))

[43] 杜修力,金瀏,黃景琦.基于擴(kuò)展有限元法的混凝土細(xì)觀斷裂破壞過(guò)程模擬[J].計(jì)算力學(xué)學(xué)報(bào),2012,29(6): 940-947.(DU Xiuli,JIN Liu,HUANG Jingqi.Simulation of meso-fracture process of concrete using the extended finite element method[J].Chinese Journal of Computational Mechanics,2012,29(6): 940-947.(in Chinese))

[44] 于紅軍.含復(fù)雜界面非均勻材料斷裂力學(xué)研究[D].哈爾濱:哈爾濱工業(yè)大學(xué),2010.

[45] 茹忠亮,申崴,趙洪波.基于擴(kuò)展有限元法的鋼筋混凝土梁復(fù)合斷裂過(guò)程模擬研究[J].工程力學(xué),2013,30(5): 215-220.(RU Zhongliang,SHEN Wei,ZHAO Hongbo.Simulation of mixed-mode fracture process of reinforced concrete beam based on extended finite element method[J].Engineering Mechanics,2013,30(5): 215-220.(in Chinese))

[46] 汪基偉,張雄文,林新志.考慮黏結(jié)滑移的平面組合式單元模型研究與應(yīng)用[J].工程力學(xué),2008,25(1): 97-102.(WANG Jiwei,ZHANG Xiongwen,LIN Xinzhi.Research and applications of a plane embedded combined element model considering bond and slip[J].Engineering Mechanics,2008,25(1): 97-102.(in Chinese))

[47] SIMONE A,WELLS G N,SLUYS L J.A novel technique for modelling interfaces in reinforced brittle materials[C]//DEBORST R,MAZARS J,PIJAUDIERCABOT G.4th International Conference on Fracture Mechanics of Concrete Structures.Leiden: A A Balkema Publishers,2001: 841-846.

[48] DEB D,DAS K C.Enriched finite element procedures for analyzing decoupled bolts installed in rock mass[J].International Journal for Numerical and Analytical Methods in Geomechanics,2011,35(15): 1636-1655.

[49] IBRAHIMBEGOVIC A,BOULKERTOUS A,DAVENNE L,et al.Modelling of reinforced-concrete structures providing crack-spacing based on X-FEM,ED-FEM and novel operator split solution procedure[J].International Journal for Numerical Methods in Engineering,2010,83(4): 452-481.

[50] 何學(xué).鋼筋混凝土擴(kuò)展有限元單元模型研究[D].南京: 河海大學(xué),2013.

[51] CONTRAFATTO L,CUOMO M,FAZIO F.An enriched finite element for crack opening and rebar slip in reinforced concrete members[J].International Journal of Fracture,2012,178(1/2): 33-50.

[52] 方修君,金峰.基于ABAQUS平臺(tái)的擴(kuò)展有限元法[J].工程力學(xué),2007,24(7): 6-10.(FANG Xiujun,JIN Feng.Extended finite element method based on abaqus[J].Engineering Mechanics,2007,24(7): 6-10.(in Chinese))

[53] 楊濤,鄒道勤.基于XFEM的鋼筋混凝土梁開裂數(shù)值模擬[J].浙江大學(xué)學(xué)報(bào):工學(xué)版,2013,47(3): 495-501.(YANG Tao,ZOU Daoqin.Numerical simulation of crack growth of reinforced concrete beam based on XFEM[J].Journal of Zhejiang University: Engineering Science,2013,47(3): 495-501.(in Chinese))

Review of extended finite element method and its application in reinforced concrete structures//

LI Jun, WANG Jiwei, LENG Fei

(CollegeofCivilandTransportationEngineering,HohaiUniversity,Nanjing210098,China)

The basic theory and application value of extended finite element method(XFEM) are introduced briefly. The level set method for describing crack, the enrichment strategy of discontinuous displacement field, and the integration scheme of XFEM are presented with an example of crack problem(strong discontinuity) of reinforced concrete(RC), and recent achievements in these fields are exposed. Thereby, a summary of the situation of applications of XFEM in fracture machanics is given. On the basis of the foregoing, the research and application development of XFEM for the analysis of reinforced concrete structure crack are presented from two aspects: concrete cracking, and bond-slip between reinforcement and concrete. It is pointed out that element models, crack treatment, steel stress at crack section and other issues have to be further investigated in order to meet the requirements of characteristics of practical projects(e.g., hydraulic reinforced concrete structures), such as large scale, complex reinforcement, and multi-axial force state.

extended finite element method; reinforced concrete structures; fracture machanics; review

國(guó)家大壩安全工程技術(shù)研究中心開放基金(NDSKFJJ1101)

李俊(1990—),男,江蘇沛縣人,碩士研究生,主要從事混凝土結(jié)構(gòu)研究。E-mail:lijunhhucn@163.com

10.3880/j.issn.1006-7647.2015.03.020

TV313

A

1006-7647(2015)03-0106-08

2014-02-07 編輯：駱超)