李子陽,李 季,向 衍

(1.南京水利科學(xué)研究院,江蘇南京 210029;2.水利部大壩安全管理中心,江蘇南京 210029;3.黃河上游水電開發(fā)有限責(zé)任公司,青海西寧 810008)

在水利工程中,裂縫是影響建筑物安全與正常運(yùn)行的不利因素之一。但由于裂縫物理特性的非線性以及尖端處應(yīng)力應(yīng)變的奇異性,使得有限元模擬計(jì)算困難,往往難以得出準(zhǔn)確結(jié)果。目前進(jìn)行裂縫穩(wěn)定性分析主要有4種判據(jù):在工程斷裂安全分析方面,由于縫端應(yīng)力強(qiáng)度因子KⅠ可以采用位移直接法計(jì)算,分析手續(xù)簡單,應(yīng)用最為廣泛[1-4]。但KⅠ是基于線彈性的公式,對(duì)彈塑性材料存在塑性區(qū)修正的問題,為滿足精度要求,分析中對(duì)縫端單元尺寸及選用外推節(jié)點(diǎn)有很多限制[5-6]。即便如此,采用靠近縫端的節(jié)點(diǎn)位移和采用遠(yuǎn)離縫端的節(jié)點(diǎn)位移進(jìn)行外推得到的KⅠ仍常存在較大差異,難以進(jìn)行規(guī)范化比較。而J積分由于是從能量角度對(duì)含裂縫物體的應(yīng)力應(yīng)變場進(jìn)行分析,可以避開裂縫尖端應(yīng)力應(yīng)變求解的困難,已成為研究彈塑性斷裂力學(xué)的主流。但J積分的研究分析多限于平面問題,三維J積分的理論雖然已經(jīng)比較成熟,但由于還存在積分路線確定及J積分有限元計(jì)算等定量實(shí)現(xiàn)上的技術(shù)問題,在水利工程中的應(yīng)用還未能開展。據(jù)此,本文就三維J積分并入有限元計(jì)算程序中的計(jì)算處理方法進(jìn)行討論,以提供一些三維J積分應(yīng)用的實(shí)現(xiàn)技術(shù)。

1 三維J積分定義

為了避開求解裂縫前緣塑性區(qū)應(yīng)力應(yīng)變場時(shí)數(shù)學(xué)上的困難,Rice[7]從能量守恒出發(fā),利用格林公式分析了平面裂縫尖端區(qū)域應(yīng)力應(yīng)變場的強(qiáng)度,首次將J積分應(yīng)用于斷裂力學(xué)。在實(shí)踐中,為解決J積分對(duì)處理空間含裂紋問題存在的困難,陳國棟等[8]利用勢能原理和格林定理導(dǎo)出了三維J積分表達(dá)式,并用高斯公式證明在滿足不計(jì)體力、小應(yīng)變以及單調(diào)加載等條件下,三維J積分值與積分曲面無關(guān)。

三維J積分可定義為以下曲面的積分:)式中:Σ為裂縫下表面到上表面任一包圍縫端的積分曲面;W為曲面上任一點(diǎn)的形變能密度函數(shù),分別為應(yīng)力和應(yīng)變),若采用線彈性本構(gòu)關(guān)系,則;n為 Σ上的向外單位法線;T為作用在Σ上沿外法線方向的張力矢量;b為縫端位移方向(表明其擴(kuò)展方向)的單位法線矢量;u為Σ上的位移矢量。

為解決分析及評(píng)判的一致性問題,宮本博[9]進(jìn)一步提出了單位長度J積分的計(jì)算形式(取裂縫擴(kuò)展方向?yàn)閤方向,如圖1所示):式中:B為計(jì)算裂縫長度;ui為i方向上的位移;i=1,2,3。

圖1 三維J積分有限元積分曲面

2 三維J積分有限元計(jì)算的程序?qū)崿F(xiàn)

J積分并入有限元計(jì)算主要存在積分路徑的搜索及依據(jù)有限元結(jié)果的積分實(shí)現(xiàn)等技術(shù)問題,對(duì)此作以下討論。

2.1 積分回路的拓?fù)渌阉?/h3>

雖然理論上J積分對(duì)回路的選擇沒有要求,但在實(shí)際的有限元模擬計(jì)算中考慮單元的離散性及計(jì)算的易實(shí)現(xiàn)性,實(shí)際的回路選擇都是以已有的線面為參考。積分路徑可以直接定義,但需要一定的計(jì)算經(jīng)驗(yàn),顯然不利于程序自動(dòng)化的實(shí)現(xiàn)。為此,根據(jù)有限元網(wǎng)格中點(diǎn)面之間連接的拓?fù)潢P(guān)系,采用拓?fù)渌阉髯詣?dòng)定義積分路線。

考慮單元的劃分及計(jì)算方便,三維J積分曲面一般選擇成由2個(gè)側(cè)面S2和S3與1個(gè)沿裂縫長度方向的主面S1組成(如圖1)。由于兩側(cè)面積分對(duì)J積分值貢獻(xiàn)很小,基本不超過1%[9],故忽略不計(jì),只對(duì)主積分面 S1進(jìn)行積分路徑尋找及曲面積分,以簡化J積分的計(jì)算。

將積分曲面尋找轉(zhuǎn)化成面上各積分回路節(jié)點(diǎn)尋找(然后對(duì)應(yīng)節(jié)點(diǎn)相連構(gòu)成面信息)的平面拓?fù)渌阉鲉栴},則按照拓?fù)渌阉鞯母拍?積分回路的定義就是根據(jù)縫端區(qū)域單元或節(jié)點(diǎn)間的拓?fù)潢P(guān)系,首先拓?fù)浒l(fā)現(xiàn)路徑所在的區(qū)域,然后通過建立節(jié)點(diǎn)間的鄰接矩陣拓?fù)渌阉骷s束化外回路。積分路徑拓?fù)渌阉鞯膶?shí)現(xiàn)步驟包括收集拓?fù)湫畔?、研究回路區(qū)域、生成節(jié)點(diǎn)鄰接矩陣和回路的自動(dòng)搜索。

2.1.1收集拓?fù)湫畔?/p>

根據(jù)有限元程序獲得的縫端區(qū)域單元及節(jié)點(diǎn)信息,收集并設(shè)定如下所需拓?fù)湫畔?①縫端節(jié)點(diǎn)及各積分回路起始、結(jié)束節(jié)點(diǎn)編號(hào);②裂縫周圍積分區(qū)域的單元面拓?fù)湫畔?包括單元面編號(hào)及單元面節(jié)點(diǎn)組成等,并以單元面信息為基礎(chǔ),利用對(duì)象選擇集函數(shù)得到邊線段信息,由此確定面中節(jié)點(diǎn)的拓?fù)溥B接關(guān)系。

2.1.2確定回路區(qū)域

尋找包含上一積分回路節(jié)點(diǎn)(縫端節(jié)點(diǎn)作為初始積分回路)且不為上一回路區(qū)域(積分回路所包圍的區(qū)域)的所有單元面,刪除面法向與縫長度方向一致的單元面,由此得到第i個(gè)回路區(qū)間,將該區(qū)間所有節(jié)點(diǎn)中上一回路區(qū)域的節(jié)點(diǎn)刪去即得到構(gòu)成該積分回路的所有節(jié)點(diǎn)。

2.1.3生成節(jié)點(diǎn)鄰接矩陣

在拓?fù)鋱D中,點(diǎn)與點(diǎn)之間的連接關(guān)系可以用鄰接矩陣來描述。即把N個(gè)點(diǎn)的拓?fù)鋱D表示為N×N個(gè)元素的對(duì)稱矩陣,行與列的元素都代表點(diǎn)的編號(hào),如果拓?fù)鋱D的i點(diǎn)與j點(diǎn)之間有邊連接的話,其鄰接矩陣A中aij元素的值取為1;否則為0。

確定回路節(jié)點(diǎn)集的鄰接矩陣步驟如下:①鄰接矩陣A的行、列數(shù)均為節(jié)點(diǎn)數(shù)m;②把矩陣A的所有元素值aij置為零;③對(duì)照回路區(qū)間線段表,若節(jié)點(diǎn) ai與aj相連的話,則在第 i行第j列賦值為1,依次填充矩陣。由此確定的鄰接矩陣為對(duì)稱矩陣。鄰接矩陣建立了節(jié)點(diǎn)間的拓?fù)潢P(guān)系,為實(shí)現(xiàn)回路的拓?fù)溥B接搜索提供了數(shù)據(jù)結(jié)構(gòu)上的支持。

2.1.4回路的自動(dòng)搜索

利用鄰接矩陣自動(dòng)搜索回路。設(shè)回路節(jié)點(diǎn)為m個(gè),為討論方便,用ni,k表示某一行的列號(hào),表示ni,k的鄰接標(biāo)號(hào),即

整個(gè)搜索過程由第1行開始,依次判斷鄰接矩陣大于0的元素作為起點(diǎn)。為討論任意性,設(shè)[i1,i2]處的值 ai1,i2大于0,記錄 i1為第 1點(diǎn),i2為第2點(diǎn),然后轉(zhuǎn)到[i1,i2]的對(duì)稱元素[i2,i1];在第 i2行搜索ai2,i1的鄰接標(biāo)號(hào)及再鄰接標(biāo)號(hào),直至找到ai2,i3＞0,把列號(hào) i3記錄為第3點(diǎn),然后再在i3行搜索第4點(diǎn)……這樣經(jīng)過若干次搜索之后,搜索到所給的結(jié)束點(diǎn)時(shí)即可結(jié)束此次搜索。搜索的具體順序可以表示成式(4)的形式。

由上述步驟,使用網(wǎng)格拓?fù)渌阉骶涂梢宰詣?dòng)定義多個(gè)不斷增加尺寸的回路區(qū)域,而這些區(qū)域的邊界節(jié)點(diǎn)拓?fù)溥B接即成為積分回路。其中,第1回路區(qū)域由所有連接到裂紋尖端單元的節(jié)點(diǎn)構(gòu)成,下一區(qū)域在前一區(qū)域的基礎(chǔ)上加上與前一區(qū)域節(jié)點(diǎn)相連的所有單元節(jié)點(diǎn)。對(duì)于每個(gè)縫端節(jié)點(diǎn)分別確定積分區(qū)域,各區(qū)域呈輻射狀增長,對(duì)于規(guī)則網(wǎng)格,相當(dāng)于創(chuàng)建了一個(gè)垂直于裂紋前緣的節(jié)點(diǎn)盤。完成所有縫端節(jié)點(diǎn)積分路徑的尋找,即可根據(jù)節(jié)點(diǎn)對(duì)應(yīng)關(guān)系確定積分曲面。為了搜索及計(jì)算方便,一般要求裂縫前沿網(wǎng)格應(yīng)由退化的規(guī)則六面體單元組成,由此三維裂縫拓?fù)渌阉鞔_定的積分回路如圖2所示。

圖2 三維裂縫拓?fù)渌阉鱆積分路線

2.2 積分的計(jì)算

在式(2)中,令

在有限元中,積分曲面 Σ被分割為一系列的單元面。如果可以求出 Σ上每個(gè)單元面的JW與JT,那么將所有的單元面結(jié)果相加,即可求得 Σ上的J。據(jù)此,將J的計(jì)算考慮如下:對(duì)于單個(gè)單元面,采用線彈性本構(gòu)關(guān)系,則式中:ne為該單元面外法線方向矢量。

對(duì)此,如果能求得單元面上的應(yīng)力、應(yīng)變函數(shù)及面邊線的空間方程,則可進(jìn)行精確計(jì)算。但有限元分析往往只給出節(jié)點(diǎn)的計(jì)算結(jié)果,函數(shù)形式還需要利用插值手段構(gòu)建;另一方面需要積分的單元面較多而且形狀可能很不規(guī)則,積分邊界的確定比較麻煩。

為了避開積分運(yùn)算并獲得較高的精度,仍然采用有限元離散的思想,即將單元面各自邊的中點(diǎn)連接拆分成4個(gè)區(qū)域(圖3),對(duì)每個(gè)區(qū)域繼續(xù)劃分還可分解成16個(gè)、64個(gè)等。對(duì)于每個(gè)小的區(qū)域,可近似認(rèn)為應(yīng)力、應(yīng)變值保持不變(取四節(jié)點(diǎn)平均值)。如此有

式中:σij(k),εij(k)分別為第 k個(gè)區(qū)域的應(yīng)力、應(yīng)變值;a(k)為第k個(gè)區(qū)域的面積。

圖3 積分剖分處理示意圖

對(duì)于JT,設(shè)n1,n2,n3為面元da外法線n的方向余弦,由力的平衡條件有

可求得

同樣地,把每個(gè)單元面劃分成多個(gè)積分區(qū)域,有

由此得

3 實(shí)例分析

某重力拱壩壩高178 m,水庫正常高水位為2600m。壩址區(qū)屬大陸性氣候,全年冷期長,氣溫變差大。自運(yùn)行以來,在壩體下游面出現(xiàn)130多條裂縫,且以水平裂縫居多,其中30多條水平裂縫貫穿壩段下游面。裂縫成因分析表明裂縫主要由低溫拉應(yīng)力引起,為Ⅰ型裂縫,比較危險(xiǎn)。為此,選取拱冠梁2560.00m高程處水平裂縫作為典型裂縫,應(yīng)用J積分法對(duì)其進(jìn)行穩(wěn)定性分析,并與應(yīng)力強(qiáng)度因子法所得結(jié)果進(jìn)行對(duì)比。由于沒有成熟的模擬裂縫自動(dòng)擴(kuò)展的有限元應(yīng)用模型,對(duì)裂縫開裂深度進(jìn)行假定(分別假定為1.0m,2.0m,3.0m,4.0m等),分別分析其穩(wěn)定性。

對(duì)大壩建立三維有限元分析整體模型,見圖4。其中水平裂縫的縫端局部細(xì)化如圖5所示(開裂深度為2.0m,其余不同開裂深度的模型類似)。裂縫區(qū)(填充區(qū)所包含區(qū)域)采用退化六面體單元,縫端加密單元尺寸0.125m左右,沿裂縫開裂方向單元尺寸定義為1m,簡化了模擬不同開裂深度時(shí)模型的重新劃分(只需將裂縫區(qū)沿開裂方向平移),裂縫區(qū)與整體采用Glub技術(shù)接觸連接。計(jì)算時(shí)采用反演和設(shè)計(jì)參數(shù),考慮材料非線性,采用Drucker-Prager屈服準(zhǔn)則。溫度場計(jì)算以水溫、氣溫和壩內(nèi)溫度計(jì)測值近似作為邊界約束條件,封拱溫度場取為7℃。由于裂縫主要是由低溫引起,為Ⅰ類拉開型,裂縫及橫縫均采用Contact單元模擬其接觸及受拉張開性態(tài)[10]。

模型計(jì)算中選取典型荷載工況為:水壓、自重、低溫等單項(xiàng)荷載工況以及水壓+自重+低溫組合工況。經(jīng)對(duì)各種工況下總體有限元模型分析,結(jié)果表明水平裂縫在“2590m水位+低溫+壩體自重”工況下梁向拉應(yīng)力最大,對(duì)穩(wěn)定最不利,因此,分析此最不利荷載工況下裂縫的穩(wěn)定性。

裂縫計(jì)算所得不同開裂深度時(shí)的JⅠ積分值與縫端應(yīng)力強(qiáng)度因子KⅠ外推值列于表1。由表1可以看出:在小范圍屈服及合理選取外推節(jié)點(diǎn)計(jì)算KⅠ的情況下,JⅠ與KⅠ具有類似的計(jì)算結(jié)果;當(dāng)水平裂縫開裂深度為4.0m時(shí),KⅠ＜KⅠC,且JⅠ＜JⅠC,說明此時(shí)裂縫達(dá)到穩(wěn)定狀態(tài)。

表1 縫端應(yīng)力強(qiáng)度因子KⅠ及 JⅠ積分值

圖6與圖7分別列出了假設(shè)開裂深度為2.0m時(shí)采用距縫端不同距離節(jié)點(diǎn)有限元結(jié)果計(jì)算出的KⅠ值及不同積分路徑的JⅠ積分值,可以看出:KⅠ有很強(qiáng)的距離依賴性,選取不同的節(jié)點(diǎn)值外推,其結(jié)果可能差別較大;而JⅠ積分值的路徑依賴性很小,基本一致。故JⅠ積分作為裂縫穩(wěn)定性與否的判斷依據(jù)更有可比性及參考性。

圖4 大壩三維有限元模型

圖5 裂縫細(xì)化有限元分析模型

圖6 距縫端不同距離應(yīng)力強(qiáng)度因子計(jì)算值

圖7 不同積分路徑JⅠ積分計(jì)算值

4 小 結(jié)

鑒于J積分從能量角度對(duì)含裂縫物體的應(yīng)力應(yīng)變場進(jìn)行分析,可以避開裂紋尖端應(yīng)力應(yīng)變求解的困難,將三維J積分計(jì)算并入有限元計(jì)算程序中,并對(duì)計(jì)算實(shí)現(xiàn)的技術(shù)方法進(jìn)行了討論,包括通過建立節(jié)點(diǎn)間的鄰接矩陣實(shí)現(xiàn)積分回路拓?fù)渌阉饕约胺e分計(jì)算的離散方法等。實(shí)例分析表明其與合理選取外推節(jié)點(diǎn)計(jì)算的KⅠ具有類似的計(jì)算結(jié)果。但J積分對(duì)單元尺寸大小要求不高,路徑依賴性也很小,作為裂縫穩(wěn)定性與否的判斷依據(jù)更有可比性及參考性,可有效提高對(duì)裂縫的分析水平。

[1]黃耀英,吳中如,顧沖時(shí),等.縫端應(yīng)力強(qiáng)度因子對(duì)網(wǎng)格尺寸的敏感性分析[J].巖石力學(xué)與工程學(xué)報(bào),2007,26(S2):41-47.

[2]李先明,秦忠國.閘墩混凝土表面溫度裂縫穩(wěn)定性斷裂力學(xué)分析[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版,2007,35(4):422-424.

[3]李子陽,王建,金永強(qiáng).大壩水平裂縫穩(wěn)定性分析[J].水力發(fā)電,2007,33(7):36-41.

[4]李同春,劉曉青,呂泰仁,等.柘溪水電站多聯(lián)體壩裂縫穩(wěn)定三維隨機(jī)分析[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版,1997,25(增刊):86-89.

[5]朱伯芳.有限單元法原理與應(yīng)用[M].北京:中國水利水電出版社,1998.

[6]于驍中,張彥秋,曹建國,等.混凝土復(fù)合型(Ⅰ、Ⅱ型)裂紋斷裂準(zhǔn)則的計(jì)算和試驗(yàn)研究[J].水利學(xué)報(bào),1982(6):27-37.

[7]RICE J R.A path independent integral and the approximate analysis of strain concentration by notches and crack[J].Appl Mech,1968,35:379-386.

[8]陳國棟,呂運(yùn)冰,汪冬生.空間體J積分表達(dá)式及其守恒性[J].武漢理工大學(xué)學(xué)報(bào):交通科學(xué)與工程版,2007,31(1):130-132.

[9]宮本博.彈塑性斷裂力學(xué)[M].楊秉憲,王幼復(fù),編譯.太原:山西人民出版社,1983.

[10]李子陽,谷艷昌,張磊.基于Marc的高拱壩蓄水期變形分析[J].水利水電科技進(jìn)展,2008,28(4):15-19.

[11]沈長松,陸紹俊,林益才.混凝土重力拱壩下游面裂縫斷裂穩(wěn)定性初析[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版,1995,23(1):22-29.

[12]KANNINEN M V F,POPELAR C H.高等斷裂力學(xué)[M].洪其麟,鄭光華,鄭祺選,等,譯.北京:北京航空學(xué)院出版社,1987