水工混凝土帶縫結(jié)構(gòu)溫度場問題的擴展有限元法

林學(xué)軍，陳 安，霍中艷，陳旭東

(1.浙江省臺州市港航管理局，臺州 317700；2.浙江海洋大學(xué)，舟山 316000；3.鄭州大學(xué)，鄭州 450001)

水工混凝土帶縫結(jié)構(gòu)溫度場問題的擴展有限元法

林學(xué)軍1，陳 安1，霍中艷2，陳旭東3

(1.浙江省臺州市港航管理局，臺州 317700；2.浙江海洋大學(xué)，舟山 316000；3.鄭州大學(xué)，鄭州 450001)

長期處在溫度荷載等作用下運行的水工結(jié)構(gòu)或多或少都會產(chǎn)生裂縫，有必要對其溫度場進行計算分析。鑒于擴展有限元法求解不連續(xù)問題的獨特優(yōu)勢，文章用其求解溫度場；首先分析給出了裂縫結(jié)構(gòu)溫度場的逼近形式，隨后推導(dǎo)了用擴展有限元法求解不連續(xù)溫度場的離散方程，其中裂縫采用等溫形式進行模擬，即溫度連續(xù)但其一階導(dǎo)數(shù)不連續(xù)，繼而建立起溫度場的擴展有限元法。并以帶縫十字形板的算例驗證了擴展有限元法求解溫度—應(yīng)力問題的正確性；最后，用該方法分析了某面板堆石壩面板裂縫在冬季降溫期的擴展情況，其結(jié)果與實際監(jiān)測成果較一致。

水工混凝土；帶縫結(jié)構(gòu)；溫度場；擴展有限元

從20世紀(jì)60年代起，混凝土結(jié)構(gòu)在施工澆筑過程中由于混凝土水化熱而產(chǎn)生的溫度場和溫度裂縫、在運營階段由于周圍環(huán)境及太陽輻射等作用產(chǎn)生溫度荷載而導(dǎo)致的開裂，一直倍受世界工程界學(xué)者的關(guān)注。大量的工程實踐和理論研究都表明，當(dāng)溫度場發(fā)生變化時，會導(dǎo)致混凝土體內(nèi)熱應(yīng)力的產(chǎn)生，溫度荷載的變化對大體積混凝土結(jié)構(gòu)的應(yīng)力狀態(tài)有顯著的影響，其導(dǎo)致的溫度應(yīng)力值可能超過其余外荷載引起的應(yīng)力值，甚至比其他各種荷載產(chǎn)生應(yīng)力的總和還要大[1]，長期處在溫度荷載作用下運行的水工混凝土結(jié)構(gòu)或多或少存在著裂縫，而其裂縫的存在不僅破壞了結(jié)構(gòu)的連續(xù)性，還擾亂了熱流的傳播，從而加劇裂縫端部的應(yīng)力集中現(xiàn)象，極易引起新的裂縫或者導(dǎo)致已有裂縫的擴展，對結(jié)構(gòu)的安全產(chǎn)生不利影響，故水工混凝土結(jié)構(gòu)工程中，溫度場的計算分析具有非常重要的作用。

為了避免常規(guī)有限元法追蹤動態(tài)裂紋時要求裂縫面與單元邊界一致、網(wǎng)格需不斷重劃分等缺點，考慮到擴展有限元法對于求解不連續(xù)問題的獨特優(yōu)勢，本文以擴展有限元法求解溫度場。對于水工混凝土帶縫結(jié)構(gòu)而言，裂縫內(nèi)介質(zhì)導(dǎo)熱系數(shù)相對周圍介質(zhì)較大或縫內(nèi)介質(zhì)溫度受邊界影響變化不大時，縫兩側(cè)溫度相等(稱等溫裂縫)，可按等溫裂縫處理，采用文獻[2]提出的結(jié)點加強方法處理；當(dāng)裂縫內(nèi)介質(zhì)的導(dǎo)熱系數(shù)相對較低時(稱絕熱裂縫)，可通過將裂縫的導(dǎo)熱系數(shù)取為較小值來模擬。文章首先分析裂縫結(jié)構(gòu)溫度場的逼近形式，然后在簡述熱傳導(dǎo)問題基本原理的基礎(chǔ)上，研究建立溫度場的擴展有限元法，并以十字板算例驗證了該方法的有效性，最后將其用于某面板堆石壩面板裂縫在降溫作用下的擴展分析。

1 溫度場的逼近形式

帶縫結(jié)構(gòu)溫度場的XFEM逼近形式與位移場類似[2-7]，可用下式表示

式中：x為點的坐標(biāo)向量；Ni(x),Nj(x)和Nk(x)為常規(guī)有限元形函數(shù)；Ti為結(jié)點溫度的常規(guī)自由度，ej和fk分別為溫度在各結(jié)點處的加強自由度；Φ(x)為貫穿單元結(jié)點加強函數(shù)、γ(x)是裂尖加強函數(shù)，而Φj和γk分別為二者的結(jié)點值；I為所有結(jié)點集合；J為所有與裂縫相交的單元結(jié)點集合；K為裂尖單元的結(jié)點集合。

裂縫的貫穿加強結(jié)點與位移模式中采用相同的加強函數(shù)Φ(x)，該函數(shù)連續(xù)但其一階導(dǎo)數(shù)不連續(xù)。等溫裂縫尖端的漸近溫度場可表示為為常數(shù)，故將選為裂尖加強函數(shù)[2]。絕熱裂縫尖端的漸近溫度場可表示為為常數(shù)，故將選為裂尖加強函數(shù)，與擴展有限元法的位移模式原理一致[3]。

2 熱傳導(dǎo)的基本原理

2.1 熱傳導(dǎo)的方程

考慮到對于均勻、各向同性固體，由于溫度升高而吸收的熱量必須等于從外界流入凈熱量與內(nèi)部的水化熱之和，也即熱量平衡，可知

式中：Q是單位時間單位體積內(nèi)發(fā)出的熱量，c和r分別為比熱和密度，t代表時間，a為導(dǎo)溫系數(shù)，a = λ / cρ。若無內(nèi)熱源，則溫度場的熱傳導(dǎo)方程為

2.2 邊值條件

為了確定所需溫度場，還須知道初始條件與邊界條件，二者合稱為邊值條件。通常認(rèn)為初始瞬時溫度分布為常數(shù)，即τ = 0時

邊界條件一般可分為四種：

（2）第二邊界條件：混凝土的表面熱流量q與時間函數(shù)已知

（3）第三邊界條件：假設(shè)通過混凝土表面的熱流量q與其表面的溫度T和氣溫aT的差成正比

式中：β是表面放熱系數(shù)。若β→∞時，轉(zhuǎn)化成第一邊界條件T=Ta；若β=0時，則轉(zhuǎn)化為絕熱狀況

（4）第四邊界條件：當(dāng)兩種不同的固體接觸時，如接觸良好，則在接觸面上的溫度和熱流均應(yīng)連續(xù)，邊界條件即為

2.3 XFEM離散方程

對于滿足(4～5)的不考慮內(nèi)熱源的二維溫度場問題，可根據(jù)變分原理等價于下述泛函的極值問題

將求解域劃分為有限個單元，泛函Ie為單元e的求解子區(qū)域ΔR上的泛函，表達如下

溫度場的計算網(wǎng)格及裂縫加強結(jié)點的選取與位移場和滲流場一致，加強結(jié)點仍分為裂縫貫穿加強與裂尖加強兩種方式[3]，兩種類型每個結(jié)點的自由度均為2個，部分雙重加強的結(jié)點自由度為3個，其他常規(guī)結(jié)點的自由度為1個。

將式（9）對單元各自由度求微分，得出離散方程如下

在式(12)與(14)中，Gia= Ni，Gib= NiΨi，Gic= Niφi。

將各單元結(jié)果加以集合，則對于求解區(qū)域的全部結(jié)點可得方程組

式中：{T }為所有結(jié)點的溫度陣列(包含加強溫度)；[H]為整體傳導(dǎo)矩陣，其元素Hij由與結(jié)點i有關(guān)的各單元hij與gij集合而成，[L]與[M]為時間項離散后的整體系數(shù)矩陣，其元素Lij與Mij分別由與結(jié)點i有關(guān)的各單元lij與mi集合而成，[F]為邊界條件形成的溫度荷載陣列；令第一類邊界上各結(jié)點的溫度等于給定值，求解方程(16)即可得到所有結(jié)點的溫度值。

表1 工況詳表Tab.1 Different conditions

圖1 十字形板模型及仿真圖Fig.1 Model and simulation diagrams of cross plate

3 算例分析

用第2部分介紹的方法及筆者在文獻[3]中關(guān)于擴展有限元法模擬裂縫擴展的方法，分析一帶有邊緣傾斜裂縫的十字形板在4種不同工況下的裂縫擴展過程，各工況的溫度和力學(xué)邊界條件詳見表1。十字形板模型的幾何尺寸及邊界條件如圖1-a所示，其中，L = 2 m，a = 0.4 m，裂縫采用絕熱加強形式，與水平軸成45°角，增量步長為0.3a；材料參數(shù)為E=2.184×105Pa，v = 0.3，α=1.67×10-5/℃。

四種工況下裂縫擴展11步后的路徑及兩種類型應(yīng)力強度因子分別如圖1-d和圖2所示。為方便比較，將Ι型與Π型裂紋尖端的應(yīng)力強度因子KΙ和KΠ均用K0=αE(πa0)0.5進行標(biāo)準(zhǔn)化，標(biāo)準(zhǔn)化之后的結(jié)果分別如表2以及圖2所示，可發(fā)現(xiàn)KΙ受邊界條件的影響較大； KΠ幾乎不受邊界條件影響，擴展一步之后趨于零；工況1下，裂紋擴展幾步之后KΙ和KΠ均近似為零，說明真實條件下的裂縫應(yīng)停止擴展；工況4下，裂紋擴展4步之后KΙ和KΠ趨于平穩(wěn)，說明真實條件下裂紋停止擴展或處于穩(wěn)定擴展?fàn)顟B(tài)。

圖2 各工況下標(biāo)準(zhǔn)化的SIFs隨裂縫擴展的變化圖Fig.2 Variations of normalized SIFs by crack propagation under different conditions

圖3 文獻[8]邊界元法計算結(jié)果Fig.3 Results of boundary element method in reference[8]

由圖1-d可知，溫度荷載對裂縫的擴展路徑影響非常大，只加溫度荷載和只加外力荷載的擴展路徑夾角大于90°，溫度和外力同時作用與單一外力或溫度作用下的裂縫擴展路徑差別也相當(dāng)大，可見裂縫問題中的溫度荷載是絕不能忽略的。圖1-b和1-c分別給出了在工況4下裂紋擴展3步的溫度場及放大1 000倍位移場。文獻[8]與本文算例模型一致，采用邊界元法計算，其中工況1,2,3,5分別與本文工況1,2,3,4一一對應(yīng)，文獻[8]各工況裂縫擴展路徑圖以及標(biāo)準(zhǔn)化之后的KΙ和KΠ分別如圖3所示，由圖可見，這與本文計算結(jié)果一致，說明了本方法的有效性。

表2 各工況下標(biāo)準(zhǔn)化的SIFs隨裂縫擴展的變化表Tab.2 Variations of normalized SIFs by crack propagation under different conditions

4 工程實例

某混凝土面板堆石壩壩頂長429 m，壩頂寬度為10 m，最大設(shè)計壩高為139 m，壩頂高程2 010 m，壩體填筑分墊層、過渡料、主堆石、次堆石等區(qū)，面板共分38塊，分縫長度約5 000 m，最大單塊長度為218 m，其中水上面板36塊；面板不等厚，厚度為0.3～0.7 m，寬度12 m，面板間設(shè)垂直縫，面板與趾板和高趾墻間設(shè)周邊縫，面板混凝土設(shè)計采用525#普通硅酸鹽水泥，為二級配混凝土，標(biāo)號為C25；該壩的正常蓄水位為2 005 m，設(shè)計及校核洪水位分別為2 005 m及2 008.28 m。經(jīng)幾年運行，堆石壩混凝土面板出現(xiàn)較嚴(yán)重的裂縫問題，且絕大部分為垂直于壩軸線且縱向分布于面板上的裂縫。多次的裂縫檢查發(fā)現(xiàn)面板混凝土受溫度影響顯著。面板在庫水位以上部分裂縫較多，且冬季產(chǎn)生的裂縫較多。2011年6月最新的現(xiàn)場檢查共發(fā)現(xiàn)水位以上大小裂縫135條，裂縫最長延伸至2 002.0 m高程水位以下，其中129 條裂縫為豎向裂縫，5條裂縫為橫向裂縫，裂縫數(shù)量比2010年10月統(tǒng)計的結(jié)果增加了12條，裂縫寬度在0.02 ～0.45 mm之間。因此有必要分析溫度變化對面板裂縫擴展的影響。下面以5#面板為例進行分析。

5#面板中間存在一條從2 005.5 m高程延伸至水位附近的縱向裂縫，檢測發(fā)現(xiàn)該裂縫在2 003.5 m高程處最深，達0.11 m，寬0.15 mm。2008年初出現(xiàn)持續(xù)的降溫過程，下面用擴展有限元法分析降溫過程2003.5 m高程處裂縫的擴展深度，為此，在2 003.5 m高程處沿壩軸線的面板法向截取一個剖面，通過建立該剖面的擴展有限元法簡化模型分析裂縫的擴展情況，剖面截取的示意圖如圖4所示，其中X向為面板寬度方向，Y向為面板厚度方向；注意到運行期面板過渡料區(qū)和墊層內(nèi)的溫度測值變幅較大，受外界氣溫影響大，而主堆石區(qū)和高程較低的過渡料區(qū)內(nèi)的溫度計測值的變幅則很小，故簡化模型不考慮堆石區(qū)溫度變化的影響。剖面的XFEM模型見圖5，模型共有960個單元，1 067個結(jié)點。

由于墊層料密實，擠壓墻與墊層料結(jié)合良好，運行期的擠壓墻壓碎后性能已經(jīng)基本與墊層相同，因此，在計算溫度應(yīng)力時僅考慮擠壓墻與面板間的接觸問題，擠壓墻與墊層之間的接觸問題不再考慮。通過在面板與擠壓墻之間設(shè)置無厚度Goodman接觸單元模擬二者的接觸問題。

圖4 面板法向剖面的截取圖Fig.4 Profile normal of the concrete panel

圖5 面板剖面的網(wǎng)格模型(單位: m)Fig.5 Mesh model of the panel profile

據(jù)面板壩所處地區(qū)的氣候特點，由氣溫觀測資料，選取旬平均氣溫最低的2008年1月下旬為典型工況之一進行分析，氣溫變化過程如圖6所示；由于相鄰面板之間溫度場影響較小，計算溫度場時，面板左右側(cè)取為絕熱，即僅考慮熱流沿面板法向流動，不考慮其縱向流動；面板混凝土的材料參數(shù)見表3，混凝土與空氣的熱交換系數(shù)β=2 009.28 kJ/m.d.℃，斷裂能Gf= 120 N/m。

接觸面單元勁度系數(shù)的選取[9]：（1）法向勁度系數(shù)，當(dāng)接觸面受壓時，為了防止兩邊二維單元在接觸處重疊嵌入，應(yīng)取較大數(shù)值，本節(jié)取107 kN/m3；若算出法向應(yīng)力為拉應(yīng)力，假定接觸面不能受拉，則取較小值，本節(jié)取10 kN/m3，使算出的拉應(yīng)力可忽略不計（2）切向勁度系數(shù)用適合增量法的經(jīng)驗公式，每級荷載更具接觸面單元的應(yīng)力狀態(tài)均用公式計算。

圖6 典型工況下氣溫過程線Fig.6 Temperature changes in typical working conditions

表3 面板混凝土材料參數(shù)表Tab.3 Material parameters of face concrete

通過對整個工況的溫度場計算及溫度計監(jiān)測資料可知，面板沿厚度方向從表面至底部受氣溫影響逐漸變?nèi)?，面板表面溫度與氣溫相同，最大溫差為-5.5℃，出現(xiàn)在2008年1月30日。

采用增量法求解冬季降溫期工況下(如圖6所示)面板裂縫的擴展過程，初始溫度場取年平均氣溫8.5℃，認(rèn)為此時面板內(nèi)部無應(yīng)力作用，所謂的增量是指某時刻溫度對比初始溫度場的增量，正負(fù)均可；若溫度增量為正，表示溫升，此時溫度應(yīng)力相比前一日有所減?。ɡ瓚?yīng)力為正），反之則溫度應(yīng)力增加。按照文獻[3]進行迭代計算，以時刻tn的溫度場為初始條件，得出時刻tn+1=tn+Δt的溫度場，以此類推，不斷求解新的溫度場和應(yīng)力場，并以最大主應(yīng)力準(zhǔn)則判斷裂縫在各時間步是否擴展，最終得出整個降溫期的擴展路徑。時間步取為1 d；通過反復(fù)試算，裂縫增量長取為0.01 m；計算得出裂縫的擴展路徑如圖7所示；2008年1月29日，裂縫達到最大深度0.27 mm；圖8為2008年1月29日裂縫穩(wěn)定后裂尖附近X向的局部應(yīng)力放大圖。

由圖7和圖8可知，溫降拉應(yīng)力使得面板縫向墊層延伸擴展，拉應(yīng)力在裂尖附近有一定的應(yīng)力集中現(xiàn)象，隨著溫度降低，裂縫不斷擴展，裂縫擴展深度為0.27 m，比檢測值大0.03 m(鉆芯取樣的檢測值為0.24 m)，可能是由于沒考慮三維情況及面板內(nèi)的鋼筋的緣故，擴展路徑與監(jiān)測結(jié)果較一致。

圖7 裂縫擴展路徑Fig.7 Crack path

圖8 裂縫穩(wěn)定后(2008年1月29日)裂尖附近的X向應(yīng)力圖(單位: MPa)Fig.8 Distribution of stress near the crack-tip in X direction when crack being stable on January 29, 2008

5 結(jié)論

鑒于溫度荷載對水工結(jié)構(gòu)此類大體積混凝土結(jié)構(gòu)的影響較大，易導(dǎo)致裂縫進一步擴展，產(chǎn)生不利影響，考慮到擴展有限元法處理不連續(xù)問題的獨特優(yōu)勢，故本文采用擴展有限元法求解水工帶縫結(jié)構(gòu)的溫度場，推導(dǎo)了用擴展有限元法求解不連續(xù)溫度場的離散方程，其中裂縫采用等溫形式進行模擬，即溫度連續(xù)但其一階導(dǎo)數(shù)不連續(xù)，并以帶縫十字形板的算例驗證了擴展有限元法求解溫度-應(yīng)力問題的正確性和有效性，最后將該方法用于某面板堆石壩面板裂縫在冬季連續(xù)降溫作用下的裂縫擴展分析，其裂縫擴展深度、路徑與實際監(jiān)測結(jié)果較一致。

[1] 龔召熊,張錫祥,肖漢江，等. 水工混凝土的溫控與防裂[M]. 北京:中國水利水電出版社,1999.

[2] Sukumar N, Chopp D L, Moes 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-47): 6 183-6 200

[3] 霍中艷,鄭東健.基于擴展有限元法的粘聚性裂縫模型的混凝土梁斷裂過程模擬[J]. 計算機輔助工程,2010,19（4）:29-33. HUO Z Y, ZHENG D J.Simulation of concrete cracking using cohesive crack model based on XFEM[J]. Computer Aided Engineering,2010,19(4):29-33.

[4] Dolbow J, Moes 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-52): 6 825-6 846.

[5] Khoei A R, Nikbakht M. An enriched finite element algorithm for numerical computation of contact friction problems[J]. International Journal of Mechanical Sciences,2007, 49(2): 183-199.

[6] Moes N, Cloirec M, Cartraud P, et al. A computational approach to handle complex microstructure geometries[J]. Computer Methods in Applied Mechanics and Engineering, 2003, 192(28-30): 3 163-3 177.

[7] Sukumar N, Chopp D L，Moes 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-47): 6 183-6 200.

[8] Prasad N N V, Aliabadi M H. Incremental crack growth in thermoelastic problems[J].International Journal of Fracture, 1994, 66:45-50.

[9] 殷宗澤. 土工原理[M]. 北京:中國水利水電出版社,2007.

Temperature field analysis of harbor engineering structure with cracks based on XFEM

LIN Xue-jun1, CHEN An1, HUO Zhong-yan2, CHEN Xu-dong3
(1. Taizhou Port and Shipping Management Bureau,Taizhou 317700,China; 2. Zhejiang Ocean University,Zhoushan 316000,China; 3. Zhengzhou University,Zhengzhou 450001,China)

For the long-term exposure to temperature load, certain parts of a concrete structure have more or less cracks during its construction and operation period. It is necessary to calculate and analyze the temperature fi eld. For the diffi culty of applying classical fracture criteria to the actual hydraulic engineering and simulating the process of cracking by conditional FEM, the XFEM was introduced into the analysis of the temperature fi eld in hydraulic structures in this paper. Firstly, the approximation form of the temperature fi eld of cracked structure was analyzed and provided, and then the discrete equation for solving discontinuous temperature fi eld was deduced, where the crack was simulated by isothermal form, namely, the temperature field is continuous but its first derivative is not continuous. The temperature fi eld by XFEM was established. And a numerical example of a cross shaped plate with slit was used to verify the validity of proposed method. Finally, this method was used to analyze the expansion of the cracks in the face slab of a concrete face rock-fi ll dam in winter. The results are consistent with the actual monitoring results.

harbor engineering concrete; structure with crack; temperature fi eld; XFEM

TU755；O 241.82

A

1005-8443(2017)03-0274-07

2016-11-25；

2017-02-16

浙江省交通廳科技計劃項目：港口碼頭健康檢測評估、修復(fù)加固技術(shù)水文水資源與水利工程科學(xué)國家重點實驗室開放基金（2015491411）

林學(xué)軍（1966-），男，浙江臺州人，高級工程師，主要從事港口航道工程建設(shè)與管理方面工作。

Biography:LIN Xue-jun(1966-),male,senior engineer.