基于三維J積分的裂縫穩定性仿真計算

李子陽,李 季,向 衍

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

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

1 三維J積分定義

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

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

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

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

2 三維J積分有限元計算的程序實現

J積分并入有限元計算主要存在積分路徑的搜索及依據有限元結果的積分實現等技術問題,對此作以下討論。

2.1 積分回路的拓撲搜索

雖然理論上J積分對回路的選擇沒有要求,但在實際的有限元模擬計算中考慮單元的離散性及計算的易實現性,實際的回路選擇都是以已有的線面為參考。積分路徑可以直接定義,但需要一定的計算經驗,顯然不利于程序自動化的實現。為此,根據有限元網格中點面之間連接的拓撲關系,采用拓撲搜索自動定義積分路線。

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

將積分曲面尋找轉化成面上各積分回路節點尋找(然后對應節點相連構成面信息)的平面拓撲搜索問題,則按照拓撲搜索的概念,積分回路的定義就是根據縫端區域單元或節點間的拓撲關系,首先拓撲發現路徑所在的區域,然后通過建立節點間的鄰接矩陣拓撲搜索約束化外回路。積分路徑拓撲搜索的實現步驟包括收集拓撲信息、研究回路區域、生成節點鄰接矩陣和回路的自動搜索。

2.1.1收集拓撲信息

根據有限元程序獲得的縫端區域單元及節點信息,收集并設定如下所需拓撲信息:①縫端節點及各積分回路起始、結束節點編號;②裂縫周圍積分區域的單元面拓撲信息,包括單元面編號及單元面節點組成等,并以單元面信息為基礎,利用對象選擇集函數得到邊線段信息,由此確定面中節點的拓撲連接關系。

2.1.2確定回路區域

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

2.1.3生成節點鄰接矩陣

在拓撲圖中,點與點之間的連接關系可以用鄰接矩陣來描述。即把N個點的拓撲圖表示為N×N個元素的對稱矩陣,行與列的元素都代表點的編號,如果拓撲圖的i點與j點之間有邊連接的話,其鄰接矩陣A中aij元素的值取為1;否則為0。

確定回路節點集的鄰接矩陣步驟如下:①鄰接矩陣A的行、列數均為節點數m;②把矩陣A的所有元素值aij置為零;③對照回路區間線段表,若節點 ai與aj相連的話,則在第 i行第j列賦值為1,依次填充矩陣。由此確定的鄰接矩陣為對稱矩陣。鄰接矩陣建立了節點間的拓撲關系,為實現回路的拓撲連接搜索提供了數據結構上的支持。

2.1.4回路的自動搜索

利用鄰接矩陣自動搜索回路。設回路節點為m個,為討論方便,用ni,k表示某一行的列號,表示ni,k的鄰接標號,即

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

由上述步驟,使用網格拓撲搜索就可以自動定義多個不斷增加尺寸的回路區域,而這些區域的邊界節點拓撲連接即成為積分回路。其中,第1回路區域由所有連接到裂紋尖端單元的節點構成,下一區域在前一區域的基礎上加上與前一區域節點相連的所有單元節點。對于每個縫端節點分別確定積分區域,各區域呈輻射狀增長,對于規則網格,相當于創建了一個垂直于裂紋前緣的節點盤。完成所有縫端節點積分路徑的尋找,即可根據節點對應關系確定積分曲面。為了搜索及計算方便,一般要求裂縫前沿網格應由退化的規則六面體單元組成,由此三維裂縫拓撲搜索確定的積分回路如圖2所示。

圖2 三維裂縫拓撲搜索J積分路線

2.2 積分的計算

在式(2)中,令

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

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

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

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

圖3 積分剖分處理示意圖

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

可求得

同樣地,把每個單元面劃分成多個積分區域,有

由此得

3 實例分析

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

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

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

裂縫計算所得不同開裂深度時的JⅠ積分值與縫端應力強度因子KⅠ外推值列于表1。由表1可以看出:在小范圍屈服及合理選取外推節點計算KⅠ的情況下,JⅠ與KⅠ具有類似的計算結果;當水平裂縫開裂深度為4.0m時,KⅠ＜KⅠC,且JⅠ＜JⅠC,說明此時裂縫達到穩定狀態。

表1 縫端應力強度因子KⅠ及 JⅠ積分值

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

圖4 大壩三維有限元模型

圖5 裂縫細化有限元分析模型

圖6 距縫端不同距離應力強度因子計算值

圖7 不同積分路徑JⅠ積分計算值

4 小 結

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

[1]黃耀英,吳中如,顧沖時,等.縫端應力強度因子對網格尺寸的敏感性分析[J].巖石力學與工程學報,2007,26(S2):41-47.

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

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

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

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

[6]于驍中,張彥秋,曹建國,等.混凝土復合型(Ⅰ、Ⅱ型)裂紋斷裂準則的計算和試驗研究[J].水利學報,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]陳國棟,呂運冰,汪冬生.空間體J積分表達式及其守恒性[J].武漢理工大學學報:交通科學與工程版,2007,31(1):130-132.

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

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

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

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