全時域初值推進邊界元法分析氯鹽侵蝕下混凝土結構服役壽命

蔣瓊明，楊綠峰，2，陳 正

（1.廣西大學 土木建筑工程學院；工程防災與結構安全教育部重點實驗室，南寧530004；2.廣西壯族自治區 住房和城鄉建設廳，南寧530028）

暴露在氯鹽環境下的混凝土結構，由于氯離子的侵入使鋼筋表面的氯離子濃度達到臨界濃度而導致鋼筋銹蝕，從而影響到混凝土結構的耐久性和安全性，因此，預測氯離子在混凝土中的濃度分布是評價混凝土結構服役壽命的一個重要環節。雖然氯離子侵入混凝土包括了擴散作用、滲透作用、電化學遷移等綜合因素的結果，但在實際的應用中通常認為擴散起主導作用，并采用基于Fick第二定律的氯離子擴散模型進行分析。余紅發等［1－3］提出了能同時考慮氯離子的結合能力、氯離子擴散系數的時間依賴性和結構微缺陷影響的混凝土結構服役壽命解析模型；Tang等［4］討論了使用過于簡單的數學解分析混凝土服役壽命可能帶來的誤差，并給出了一個改進的解析分析方法。雖然解析方法計算量少，使用簡單，但由于解析解常受到規則邊界條件和材料均勻性的局限，所以實際上常常需要借助于有限元等數值方法。Han［5］考慮結合氯離子和蒸發水對氯離子擴散系數的影響，利用有限元方法分析了混凝土中氯離子濃度的分布；Var Dimitri等［6］應用有限差分法對氯離子引起混凝土中鋼筋腐蝕的初始時間進行了評估。楊綠峰等［7］提出了氯離子擴散場的計算長度，并據此建立了混凝土中氯離子擴散分析的補償長度理論［8］及服役壽命分析的邊界元法［9］。有限元法適應性強，但在計算過程中通常在時間域和空間域中需要同時采用細密離散網格，計算量較大、效率不高。邊界元法的未知量僅存在于邊界結點上，可以使待求解問題的維數降低一維，未知量較少，計算效率得到很大提高。

使用邊界元法解決瞬態擴散問題時，可采用與時間相關［7，9］及與時間不相關［11］2種基本解。楊綠峰［7］和郭力［12］分別采用與時間相關和與時間不相關的基本解建立了混凝土中氯離子擴散分析的邊界元法。前者可以用較少的時間步得到精度較高的計算結果［8］，所以更適合于在氯離子擴散分析中應用，而且可以選擇采用2種不同的時間進程解決時域積分問題［13］：分時段逐步推進法和全時域初值推進法，前者將時間域離散后逐段建立邊界元法計算格式，并以前一時間段的計算結果作為下一時間段的初始值；后者總是從初始時刻開始積分，不需要對每個時間步的結果求域內積分，保證了邊界元法計算模型僅在擴散域邊界上離散，提高計算效率［14］。

本文采用與時間相關的基本解，并利用氯離子擴散場補償長度理論［8］建立了混凝土中氯離子擴散分析的邊界元法，分別利用分時段逐步推進法和全時域初值推進法作為邊界元法的時間進程，比較了2種時間進程在氯離子擴散問題中的特點，結果表明：基于全時域初值推進法建立混凝土中氯離子擴散分析的邊界元法，可以避免分時段逐步推進法在氯離子擴散域內產生的積分項和域內積分問題，從而簡化計算格式，顯著提高計算效率，而且能夠保證邊界元模型中時間域離散網格和計算格式的收斂性，提高計算結果的穩定性。在此基礎上，進一步研究了混凝土結構的服役壽命，分析了氯離子擴散維數對混凝土結構腐蝕過程與服役壽命的影響。

1 混凝土中氯離子擴散分析的邊界元法

基于Fick第二定律的二維擴散方程可以表示為：

其邊界條件和初始條件為：

式中：C為混凝土中氯離子濃度（%，氯離子質量占混凝土質量的百分比，本文中氯離子濃度均為此單位）；q為氯離子的濃度梯度，是C的方向導數；t為混凝土暴露于氯離子環境中的時間（a）；D為混凝土中氯離子擴散系數（mm2／a）；Cs、qs分別為混凝土表面氯離子濃度、濃度梯度；C0為混凝土內部初始氯離子濃度；Γ1、Γ2和Ω分別為混凝土的第一、二類邊界和域內。

將混凝土中氯離子擴散的時間域［t0，tn］離散成n個子域［tk－1，tk］，（k＝1，2，…，n），且：

依據文獻［7］，可以針對時間子域 ［tk－1，tk］，將式（1）轉化為如下積分方程：

式中：AP為與源點P點位置有關的系數，當點P位于氯離子擴散域內部時，AP＝1，反之，當源點P位于擴散域邊界上時，AP＝，此時β表示擴散域邊界上P點兩側切線的內夾角；為tk時刻混凝土中源點P處的氯離子濃度；C＊、q＊是氯離子擴散問題的基本解［7］。

可以看出，式（4）中每一項積分的起始點都是時間子域的初始時刻tk－1，并將前一個時間子域的計算結果作為該子域的初始值，屬于分時段逐步推進法。由于等號右側第二項屬于域內積分，因而該邊界元模型通常需要開展域內離散，不僅增加計算量，而且有可能導致計算誤差積累，降低計算精度，消減了邊界元法的優勢。

2 全時域初值推進的氯離子擴散分析邊界元法

為了克服分時段逐步推進法存在的問題，這里采用全時域初值推進法建立混凝土中氯離子擴散分析的邊界元法模型，即總是從氯離子擴散的時間域［t0，tn］的初始時刻t0開始，對于時域中的任意時刻tk，可建立邊界積分方程：

上式等號右側的域內積分項只存在于時間域的初始時刻t0，為了避免域內離散，需要消除該項，此時可以利用

代替式（1）、（2）中的氯離子濃度C，建立混凝土中氯離子擴散分析的控制方程及其邊界條件和初始條件：

其邊界條件和初始條件為：

根據上述兩式，可以建立關于CI的氯離子擴散分析的邊界積分方程：

如此即可消除積分方程中的域內積分項，從而充分發揮邊界元法的優勢。

在式（9）基礎上，建立混凝土中氯離子擴散分析的邊界元法計算模型。首先根據補償長度理論［8］確定計算模型的幾何尺寸，并沿該模型的邊界離散出N個線性單元，同時將氯離子擴散的時間域［t0，tn］離散成n個等步長子域，可得邊界元法控制方程：

式中，｛CI和 ｛qI分別表示ts時刻單元結點上的氯離子濃度CI及其梯度qI。

其中，［N］＝ ［N1N2］表示單元形函數矩陣；C和q分別由C＊和q＊在時間子域上積分得到：

這里：

將各邊界單元的控制方程集成總體方程，據此可求得邊界上各結點處的CI和qI，進而可得到域內點的氯離子濃度值。

3 氯離子侵蝕下混凝土結構服役壽命分析的邊界元法

處于氯鹽環境下混凝土結構的服役壽命是指結構建成開始使用到結構失效的時間過程，包括誘導期、發展期和失效期3個時期。其中，誘導期是指混凝土內鋼筋表面氯離子濃度達到臨界值所需的時間；發展期是指從鋼筋銹蝕到混凝土保護層脹裂所需的時間；失效期是指從保護層開裂到結構失效所需的時間。文獻中通常將誘導期作為氯鹽環境下鋼筋混凝土結構的服役壽命［9］，而將發展期和失效期作為結構的安全儲備。因此這里將誘導期壽命近似視為混凝土結構的服役壽命。

根據混凝土結構中鋼筋腐蝕的閾值濃度Ccr和鋼筋表面某點 （xd，yd）處t時刻的氯離子濃度C（xd，yd，t），可以建立混凝土結構腐蝕的極限方程：

當g （Ccr，C）＝0，即C（xd，yd，T）＝Ccr時，鋼筋開始銹蝕，稱T為氯鹽腐蝕下混凝土結構的服役壽命。

4 算例分析

算例1：使用文獻［7］中的算例，邊長為1 000mm的鋼筋混凝土方板，厚度為100mm，將其中的4個面用環氧樹脂封閉，沿厚度方向留下2個相鄰面AB和AD暴露在含氯離子的溶液中，如圖1所示?；炷恋?個暴露面上的表面氯離子濃度Cs＝1%，初始氯離子濃度為0，混凝土的氯離子擴散系數D為1.0×10－12m2／s。分別采用文獻［7］和本文中2種時間進程的邊界元法格式計算混凝土板內部對角線處的氯離子濃度，并比較2種時間進程對計算結果的影響。

計算擴散時間分別為10a、30a時混凝土中的氯離子濃度，根據混凝土中氯離子擴散分析的補償理論［7－8］，邊界元計算模型的長度Lt：

擴散時間為10、30a時的補償長度分別取為L10＝70mm、L30＝110mm，兩者都小于構件的實際長度，這里取計算長度建立邊界元法計算模型，利用全時域初值推進邊界元法計算混凝土板對角線AC上的氯離子濃度值。計算模型中，時間域被均勻劃分為5個子域。為了分析模型的收斂性，將模型的全部邊界分別離散為8、16和32個等長度單元，并將計算結果與有限元相比較，有限元模型的離散網格為10mm×10mm，時域離散步長為Δt＝0.5a。結果如圖2所示。從中可以看出，全時域初值推進邊界元法計算結果與有限元法吻合良好，將邊界離散為16個邊界單元時即可獲得較好的計算精度。

圖1 氯離子二維擴散計算模型

圖2 離散邊界元數量對計算結果的影響

為了對比全時域初值推進法和分時段逐步推進法對邊界元法計算結果和計算效率的影響，這里分別利用兩種時間進程方法計算擴散時間為30a時，混凝土板對角線AC上坐標為（40mm，40mm）處的氯離子濃度，并根據有限元法結果分析計算精度和計算效率。2種方法取同樣的BEM離散模型，即將邊界離散為16個邊界單元。另外，由于分時段逐步推進法含有域內積分項，需將擴散域離散為10mm×10mm的域內單元，并采取高斯積分的方法加以處理。計算結果如圖3、表1所示，從中可以看出時間域劃分的時段數量對2種方法計算結果的精度和效率有明顯影響。對于全時域初值推進法，邊界元法計算用時少，且計算結果隨著時段數量的增加而快速收斂。而分時段逐步推進法耗時隨著時段數量的增加而顯著增加，而且在時段不超過5個時，計算結果快速向有限元結果靠近，但當時域離散的時段繼續增加時，誤差反而快速增大，計算結果不收斂。

圖3 2種時間進程的計算精度比較

表1 2種時間進程的計算耗時比較

根據對BEM模型的分析，分時段逐步推進法產生誤差的原因在于式（4）右側第2項的域內積分。采用常規的高斯積分法計算域內積分時，由于該積分核函數C＊在源點處的值會隨著時段增加、時間步長減小而迅速增大。擴散系數D取單位值時，C＊值在源點附近的值如圖4所示，從中可以看出，積分核函數在源點處的值隨著時間步長的減小而快速增大，且增大部分主要集中在源點附近很小的區域內，造成函數曲線越來越不光滑，從而導致常規的高斯積分法產生誤差。此時，即便增加積分點，也難以有效提高計算精度，如圖5所示。

盡管可以采取細分域內離散網格或強行令積分核函數在源點處的值為零［10］等方法解決該問題，但前者會弱化邊界元法不需要域內離散的優點，增大計算工作量；后者并沒有全面解決核函數在源點附件區域過分增大的問題。

圖4 取不同時間步長時的基本解特征

圖5 取不同高斯積分點個數對相對誤差的影響

算例2 考慮一個截面尺寸500mm×500mm的方形截面混凝土柱，其鋼筋配置情況如圖6所示，鋼筋的保護層厚度為40mm。混凝土柱的4個側面都暴露于氯鹽環境中，設混凝土柱表面氯離子濃度Cs＝0.8%，初始氯離子濃度0，臨界氯離子濃度為0.22%［15］，混凝土擴散系數為1.0×10－12m2／s。這里采用本文建立的全時域初值推進邊界元法預測該混凝土柱的服役壽命。

圖6 混凝土柱截面及配筋

在柱的角部存在氯離子二維擴散現象，所以可以從角部取出一個邊長為Lt＝140mm的混凝土試塊，建立邊界元法計算模型，其中Lt為計算長度，按照50a結構服役時間和式（17）確定，氯離子沿AE邊、AF邊進入混凝土試件內，如圖7（a）所示；而在遠離拐角的中部將發生氯離子一維擴散，因此從邊中部取出長度為Lt＝140mm的混凝土試塊，建立邊界元法計算模型，氯離子沿GH邊進入混凝土試件內，如圖7（b）所示。這里利用全時域初值推進邊界元法計算兩個模型中鋼筋周圍的氯離子濃度，并和有限元法計算結果相比較，結果如圖8所示。有限元法計算模型中，氯離子擴散域采用10mm×10mm的離散單元網格，采用Δt＝0.5a的步長離散時域；邊界元法計算模型中，全部邊界離散為40個等長單元，并采用和有限元法相同的時域離散方式。從圖8可以看出，本文建立的全時域初值推進邊界元法與有限元計算結果基本吻合，具有較高的計算精度。而且，盡管二個模型中鋼筋與構件表面的距離都是40mm，但氯離子濃度及其隨時間變化的規律有明顯區別，二維擴散區鋼筋表面的氯離子濃度明顯高于一維擴散區，而且隨著擴散時間的增長，二者之間的差別逐漸增大。

圖7 混凝土柱中氯離子擴散計算模型

圖8 一、二維擴散條件下鋼筋表面的氯離子濃度

根據式（16），利用全時域初值推進邊界元法分別計算了一維、二維擴散區鋼筋銹蝕時混凝土柱的服役壽命，并將計算結果與有限元法相比較，詳見圖9。圖中結果再次表明，本文方法與有限元法結果基本吻合，具有很高的計算精度。而且可以看出，二維擴散條件下，混凝土結構的服役壽命較一維擴散條件大大降低，因此實際設計中，采用二維氯離子擴散模型進行計算更符合真實情況，才能保證鋼筋混凝土結構的耐久性。

圖9 混凝土柱的服役壽命預測值

5 結 語

建立了混凝土中氯離子擴散分析的全時域初值推進邊界元法，研究了分時段逐步推進邊界元法計算結果不隨時域細分而收斂的原因。與分時段逐步推進邊界元法相比，全時域初值推進邊界元法可以消除氯離子擴散域內的積分項，避免域內離散，不僅簡化了計算格式，而且提高了迭代過程的穩定性和計算結果的精度，減少了計算耗時。研究結果表明，氯離子擴散維數對氯離子濃度分布和混凝土結構服役壽命有顯著影響，應在結構設計中加以考慮。

［1］余紅發，孫偉，鄢良慧，等.混凝土使用壽命預測方法的研究I—理論模型［J］.硅酸鹽學報，2002，30（6）：686－690.Yu H F，Sun W，Yan L H，et al.Study on prediction of concrete servive life I－theoretical model［J］.Journal of the Chinese Ceramic Society，2002，30（6）：686－690.

［2］余紅發，孫偉，鄢良慧，等.混凝土使用壽命預測方法的研究II—模型驗證與應用［J］.硅酸鹽學報，2002，30（6）：691－695.Yu H F，Sun W，Yan L H，et al.Study on prediction of concrete service life II －model's examination and application ［J］.Journal of the Chinese Ceramic Society，2002，30（6）：691－695.

［3］余紅發，孫偉，鄢良慧，等.混凝土使用壽命預測方法的研究III—混凝土使用壽命的影響因素及混凝土壽命評價［J］.硅酸鹽學報，2002，30（6）：696－701.Yu H F，Sun W，Yan L H，et al.Study on prediction of concrete service life III －evaluation of influencing factors and service life ［J］.Journal of The Chinese Ceramic Society，2002，30（6）：696－701.

［4］Tang L P，Joost G.On the mathematics of timedependent apparent chloride diffusion coefficient in concrete［J］.Cement and Concrete Research，2007，37（4）：589－595.

［5］Han S H.Influence of diffusion coefficient on chloride ion penetration of concrete structure［J］.Construction and Building Materials，2007，21：370－378.

［6］Var Dimitri V，Trapper P A.Probabilistic evaluation of initiation time of chloride induced corrosion ［J］.Reliability Engineering ＆ System Safety，2008，93：364－372.

［7］楊綠峰，陳正，王燚，等.混凝土中氯離子二維擴散分析的邊界元法［J］.硅酸鹽學報，2009，37（7）：1110－1117.Yang L F，Chen Z，Wang Y，et al.Boundary element method for analysis of two－dimensional chloride diffusion in concrete［J］.Journal of the Chinese Ceramic Society，2009，37（7）：1110－1117.

［8］楊綠峰，陳正，劉鴻亮，等.混凝土中氯離子時變擴散的邊界元分析［J］.土木建筑與環境工程，2011，33（4）：73－79.Yang L F，Chen Z，Liu H L，et al.Boundary element method for time－dependent chloride diffusion and the chloride distribution in concrete ［J］.Journal of Civil，Architectural ＆ Environmental Engineering，2011，33（4）：73－79.

［9］陳正，楊綠峰，馮慶革，等.高性能混凝土的氯離子擴散及服役壽命研究［J］.建筑材料學報.2010，13（2）：222－227.Chen Z，Yang L F，Feng Q G，et al.Study of chloride diffusion of high－performance concrete and its service life by the boundary element method ［J］.Journal of Building Material，2010，13（2）：222－227.

［10］Mohammadi M，Hematiyan M R，Marin L.Boundary element analysis of nonlinear transient heat conduction problems involving non－homogenous and nonlinear heat sources using time－dependent fundamental solutions［J］.Engineering Analysis with Boundary Elements，2010，34：655－665.

［11］Singh K M，Kalra S.Time integration in the dual reciprocity boundary element analysis of transient diffusion ［J］.Engineering Analysis with Boundary Elements，1996，18：73－102.

［12］Guo L，Chen T，Gao X W.Transient meshless boundary element method for prediction of chloride diffusion in concrete with time dependent nonlinear coefficients［J］.Engineering Analysis with Boundary Elements，2012，36：104－111.

［13］Feng Z X，Tang S W，Li Z X.A new time marching BE scheme for transient heat transfer with a thermal moving boundary ［J］.Engineering Analysis with Boundary Elements，1998，22：71－75.

［14］Wrobe L C，Aliabad M H.The boundary element method：applications in thermo－fluids and acoustics［M］.Chichester：Johu Wiley＆Sons，2002.

［15］王曉舟，金偉良，金立兵.基于廣義擴散方程海工混凝土結構耐久壽命概率評估方法［J］.海洋工程.2009，27（3）：85－90.Wang X Z，Jin W L，Jin L B.A probabilistic approach of life prediction of marine concrete structures based on GDF［J］.The Ocean Engineering，2009，27（3）：85－90.