(武漢大學 水資源與水電工程科學國家重點實驗室，武漢 430072)

1 研究背景

鈣溶出性侵蝕是水泥基材料耐久性病害的一類，又稱溶蝕[1]，根據水泥基材料所處的環(huán)境的不同，可分為接觸溶蝕與滲透溶蝕[2]。接觸溶蝕是指水泥基材料在與水接觸過程中，水頭較低，水泥基材料內部未形成滲流作用或滲流作用效應可以忽略時，孔隙液中的鈣離子在濃度梯度作用下，不斷向外擴散，打破了固液兩相鈣濃度間的平衡關系，使得固相鈣不斷溶出的過程[3]。

在自然條件下，接觸溶蝕過程十分緩慢，溶蝕鋒面由表及里移動緩慢[4-5]。據報道，建筑物與水接觸100 a，溶蝕鋒面僅向前移動5～10 mm[6]，這一溶蝕特性對水泥基材料鈣溶蝕性能的研究增加了一定的難度。為了能在短時間內研究溶蝕對水泥基材料各方面性能的影響，加速溶蝕試驗和溶蝕數值模擬成了眾多學者研究溶蝕的2種方法。在數值模擬方面，國內外學者做了大量的研究。Reardon等[7]研究了水泥水化物與孔隙液中離子間的熱動力平衡關系；Berner[8]建立了孔隙液中鈣離子濃度與水基材料中鈣含量之間的固液平衡關系；Gérard等[9]修正了固液平衡關系，建立起簡化的一維鈣溶蝕數學模型；Nguyen等[10]結合尺度均質化分析建立了混凝土鈣離子浸出模型，引入孔隙曲折率來考慮骨料對溶蝕的影響；Wan等[11]通過試驗得到了波特蘭水泥基材料在6 mol/L的硝酸銨溶液中的固液平衡關系曲線，并據此建立了加速水泥凈漿溶蝕的一維數學模型[6]； Yu等[12]耦合化學動力學與化學熱力學提出了針對采用6 mol/L的硝酸銨溶液加速接觸溶蝕的一類新的模擬方法。

縱觀以上研究發(fā)現，水泥基材料接觸溶蝕數值模擬方法雖然較多，但大多以一維展開，求解方法以有限差分法為主，對于采用三維有限元方法模擬接觸溶蝕尚未見報道。對此，本文針對水泥基材料接觸溶蝕的三維有限元數值模擬展開研究，擬通過建立三維接觸溶蝕數學模型，采用FORTRAN語言設計有限元計算程序，通過數值計算研究各溶蝕特征量的時空分布規(guī)律。

2 三維鈣溶蝕數學模型

相關研究表明，水泥基材料在溶蝕過程中，鈣、硫、鉀、鋁等元素為主要溶出元素，而氫氧化鈣、水化硅酸鈣凝膠又是水泥的主要水化產物，因此，以鈣元素作為研究對象，可以較好地反映水泥基材料性能退化情況[13]。

參考現有的研究成果，本文根據鈣離子傳輸數學模型、固液平衡關系、孔隙率演化模型、有效擴散系數模型建立三維鈣溶蝕數學模型，并作如下假設：①孔隙液中的鈣離子與固相鈣含量之間瞬間達到平衡狀態(tài)，即忽略固體鈣溶解的時間效應;②不考慮孔隙液中鈣離子與固體骨架產生吸附與解吸作用;③不考慮溫度對溶蝕的影響;④材料各向同性。

圖1 表征體元內的質量守恒Fig.1 Conservation of mass in representative elementary volume

2.1 鈣離子傳輸數學模型

2.1.1 鈣離子傳輸控制方程

在求解域內任取一

個長、寬、高分別為Δx，Δy，Δz的六面體表征體元(REV)，如圖1所示。在域內僅考慮鈣離子擴散作用，根據表征體元內鈣離子的質量守恒關系，結合非克(Fick)第二定律可推導出水泥基材料溶蝕過程中鈣離子傳輸控制方程。

如圖1所示的表征體元，單位時間內x方向進入的鈣離子通量為

(1)

式中：Jx為x方向進入的鈣離子通量；Ax為表征體元側面面積；φ為孔隙率；Dx為x方向擴散系數；C為孔隙液中鈣離子濃度。

單位時間內x方向從表征體元流出的鈣離子通量為

(Jx+ΔJx)Ax=

(2)

根據式(1)和式(2)可得，單位時間內在x方向從表征體元進入和流出的鈣離子通量的差值即為凈有流入量，如式(3)所示。

同理可得，單位時間內在y方向、z方向的凈有流入量，分別如式(4)和式(5)所示。

式中：Jy為y方向進入的鈣離子通量；Ay為表征體元側面面積；Dy為y方向擴散系數；Jz為z方向進入的鈣離子通量；Az為表征體元側面面積；Dz為z方向擴散系數。

在溶蝕過程中，固相鈣不斷向孔隙液中溶解，假設固相鈣的摩爾濃度為CS，單位時間內表征體元孔隙液內增加的鈣離子物質的量ΔQ1為

(6)

式中t為溶蝕時間。

水泥基材料為多孔介質，假設孔隙內充滿孔隙液，則表征體元中孔隙液的體積為φΔxΔyΔz，由此可得表征體元內孔隙液中鈣離子含量為φΔxΔyΔzC，單位時間內表征體元中鈣離子的物質變化量ΔQ2為

(7)

根據質量守恒定律，單位時間內x，y，z方向的鈣離子凈有流入量與孔隙液中因固相鈣溶解的鈣離子增加量之和等于單位時間內表征體元中鈣離子的物質變化量，即

(8)

需要指出的是式中Dx，Dy，Dz為鈣離子在孔隙液中不同方向的擴散系數，經體積均勻化后得到有效擴散系數，即

(9)

將式(9)代入式(8)整理可得鈣離子三維傳輸控制方程，即

(10)

2.1.2 定解條件

求解上述控制方程組，給定初始條件與邊界條件如下所述。

(1)初始條件：

C(x,y,z,t0)=C0。

(11)

式中C0為初始時刻孔隙液中鈣離子濃度。

(2)邊界條件：

①定通量邊界條件為

(12)

式中：q為邊界上已知鈣離子的擴散通量；nx，ny，nz為邊界的單位外法線方向余弦分量。

②離子交換邊界條件[4, 14]。在以往的研究中，鈣溶蝕數學模型邊界條件多采用第一類邊界條件，即0濃度邊界條件，這一類邊界條件在模型求解過程中，邊界各點孔隙液濃度始終保持為0。由固液平衡關系可知，邊界各點的固相鈣濃度也始終保持為0，即邊界各點瞬間完成鈣溶出過程，這顯然與事實不符；實質上，隨著靠近邊界側的求解域內的濃度梯度變化，邊界各點的孔隙液濃度也在不斷發(fā)生變化。因此，采用物質交換邊界條件更為符合實際情況，這類邊界條件認為經過邊界的鈣離子通量與邊界孔隙中鈣離子濃度C與環(huán)境水中鈣離子濃度C1之差成正比，其表達式如式(13)所示。采用這類邊界條件不僅可以反映邊界濃度變化，還可以通過邊界濃度求解鈣離子累積溶出量。本文采用該類邊界條件作為接觸邊界條件,即

式中：kT為邊界上孔隙液與環(huán)境水之間的傳質系數；C1為環(huán)境水中鈣離子濃度。

2.2 固液平衡關系

控制方程式的固相鈣濃度CS與液相鈣濃度C并非獨立變量，相關研究表明，CS與C兩者具有一定的熱動力平衡關系[7-8]。Buil’s等提出了一種便于計算的等溫平衡關系模型，Nakarai等進一步修正了Buil’s模型得到固液平衡曲線[4]，Berner通過試驗數據驗證了該模型的合理性[8]。本文引用該固液平衡曲線得到C與CS間的關系，即

式中：x1為C-S-H(水化硅酸鈣凝膠)開始快速溶解時孔隙液中的鈣離子濃度；x2為Ca(OH)2完全溶解時孔隙液中的鈣離子濃度；Csatu為孔隙液中鈣離子飽和濃度；CCSH為水泥基材料中C-S-H的摩爾濃度；CCH為水泥基材料中Ca(OH)2的摩爾濃度。

根據固液平衡關系，鈣離子傳輸控制方程中固相鈣濃度可轉化為孔隙液中鈣離子濃度，進一步得到含單一變量的控制方程，即

(15)

2.3 孔隙率演化模型

2.3.1 初始孔隙率計算

水泥基材料內部孔隙分為毛細孔與凝膠孔，即

φp0=φc0+φg0。

(16)

式中：φp0為水泥漿體的初始孔隙率；φc0為水泥漿體初始毛細孔隙率；φg0為水泥漿體初始凝膠孔隙率。φc0與φg0根據POWERS模型計算可得，即

(17)

式中：W/C為水灰比；α為水化度，可由式(18)計算。

(18)

根據水泥漿體在水泥基材料的體積占比，可計算水泥基材料的初始孔隙率，即

φ0=fp/mfm/cφp0。

(19)

式中：φ0為水泥基材料初始孔隙率;fp/m為水泥漿在砂漿中的體積分數；fm/c為砂漿在混凝土中的體積分數。

2.3.2 孔隙率演化

隨著溶蝕過程的進行，Ca(OH)2與C-S-H等水化產物不斷向孔隙液中溶解，使得孔隙體積不斷增大，因而材料孔隙率也不斷增大?？紤]到C-S-H溶解過程較為復雜，假設C-S-H以Ca(OH)2形式溶解[4, 6]，據此可建立孔隙率演化模型，即

(20)

式中：Δφ為孔隙率變化量；ΔCS為固相鈣濃度變化量；MCH為Ca(OH)2的摩爾質量；ρCH為Ca(OH)2的密度；CS0為初始固相鈣濃度。

假設Ca(OH)2的溶解主要造成毛細孔的增加，而C-S-H的溶解在引起凝膠孔增大的同時，也會造成毛細孔的增大[5]。引入分配系數γ來考慮C-S-H的溶解對毛細孔隙率的影響，因此，可用式(21)、式(22)分別描述毛細孔隙率與凝膠孔隙率的變化。

(21)

(22)

2.4 擴散系數演化模型

目前，有效擴散系數模型均建立在不含粗骨料的水泥凈漿或砂漿的基礎上，對于含有粗骨料的水泥基材料無法適用。相關研究表明[15]，粗骨料對于擴散系數的影響主要有2方面：一方面，粗骨料所在區(qū)域無擴散現象發(fā)生；另一方面，粗骨料與水泥漿間存在界面過渡區(qū)，會造成孔隙增大，進而加劇擴散。為此，相關研究引入孔隙曲折度來考慮粗骨料對擴散系數的綜合影響，建立水泥基材料的通用擴散系數模型，如式(23)所示。

(23)

3 有限元計算方法

3.1 計算格式

采用三維8節(jié)點等參單元離散空間域，試探函數為式(24)所示。

(24)

式中：C代表空間域上任意一點的濃度；Nm為等參單元的形函數；Cm為等參單元的第m個節(jié)點濃度；Hm為等參單元的第m個節(jié)點總水頭；M為8。

根據Galerkin有限元計算方法與隱式差分方法分別對空間域、時間域進行離散，得到有限元計算格式，如式(25)所示。

(25)

其中：

[DC]=

3.2 計算方法

3.2.1 協(xié)調質量矩陣處理

對于時間域離散，筆者曾采用過顯式差分法、中心差分法、隱式差分法3種方法，計算過程中發(fā)現：

(1)采用中心差分格式與顯式差分格式時，只能選取較小的時間步長，才能保證數值計算結果的穩(wěn)定性，若時間步長過大，單元部分節(jié)點將出現無意義的負濃度。另外，這2種格式均不能有效地減少數值振蕩問題。

(2)采用隱式差分格式，對減小數值振蕩具有較為顯著的效果，而且不會因時間步長的變大節(jié)點出現無意義的負濃度，但時間步長較小時，也會出現較大的數值振蕩現象。

綜合考慮，為避免在小時間步長下出現數值振蕩，筆者在隱式差分格式離散時間域的基礎上，采用集中質量矩陣代替協(xié)調質量矩陣，即對式(25)中的[DC]進行處理，具體處理方法詳見文獻[16]。

3.2.2 程序實現

采用FORTRAN語言進行程序設計，矩陣元素統(tǒng)一采用一維數組存儲，對矩陣應用稀疏矩陣行壓縮存儲方案，線性方程組求解采用超松弛預處理共軛梯度法(SSOR-PCG)。

4 數值驗證

本文分別對文獻[6]中一維加速溶蝕試驗、文獻[17]中二維加速溶蝕試驗開展三維數值模擬，通過數值模擬結果與試驗結果、有限差分法計算結果進行對比，驗證三維模型及有限元計算程序的正確性。

4.1 一維溶蝕模擬

引用文獻[6]中的加速水泥凈漿接觸溶蝕試驗資料，進行數值模擬。試件尺寸為40 mm×40 mm×20 mm，水灰比為0.35，在標準條件養(yǎng)護10個月后，采用6 mol/L的NH4NO3溶液進行加速水泥凈漿試塊溶蝕試驗。

假設初始時刻的試塊孔隙中充滿孔隙液，且孔隙液中鈣離子濃度為飽和濃度，如式(26)所示。基于本文所建立的數學模型，根據文獻中的試驗條件給出邊界條件，如式(27)所示，相應的邊界條件示意圖如圖2所示。

C(x,y,z,t0)=Csatu。

(26)

(27)

圖2 邊界條件示意圖Fig.2 Schematic diagram of boundary conditions

通過水泥水化模型計算可得到試塊中C-S-H與Ca(OH)2的摩爾含量，加速溶蝕過程中的固液平衡關系曲線其他參數取值與文獻[6]相同，如表1所示，模型計算其他參數如表2所示。

表1 固液平衡曲線參數Table 1 Parameters of solid-liquid equilibrium curve mol/m3

表2 模型計算基本輸入參數Table 2 The basic input parameters of FEM program

注:kT取值參考文獻[18]。

圖3(a)為本文模擬的溶蝕深度與文獻[6]試驗數據及有限差分法模擬的結果對比，從圖3(a)可知，本文模擬(FEM)結果與有限差分法模擬(FDM)結果具有較好的一致性，兩者與試驗結果在試驗前期吻合較好，后期具有一定的誤差；經分析，這可能是由于經過一段浸蝕時間后，溶蝕介質NH4NO3溶液濃度有所降低，減緩了試件溶蝕進程[6]。

圖3(b)—圖3(d)分別為本文計算結果與文獻采用有限差分法計算結果的對比情況。從孔隙液鈣離子濃度分布、固相鈣濃度分布及孔隙率分布3個方面對比可知，本文模型計算結果與文獻[6]模擬結果具有較好的一致性。

圖3 模擬溶蝕特征量與文獻[6]數據對比Fig.3 Comparison between modeling the spatial distributions of characteristic variables and the data from Ref.[6]

4.2 二維溶蝕模擬

引用文獻[17]中的接觸溶蝕試驗資料，對試驗中水泥凈漿試塊進行加速溶蝕數值模擬。試驗中水泥凈漿試塊尺寸為10 mm×10 mm×2 mm，水灰比為0.5，采用6 mol/L的NH4NO3溶液進行加速水泥凈漿溶蝕試驗。由于試件上下側與骨料接觸，本文假設試件這兩側不與NH4NO3溶液接觸，即認為該試驗為二維溶蝕試驗。邊界條件示意圖如圖4所示。

選取中間截面進行溶蝕深度分析，溶蝕深度隨時間變化關系如圖5所示，與文獻[17]試驗數據對比發(fā)現，本文模擬計算結果與試驗數據吻合良好。

圖4 邊界條件示意圖Fig.4 Schematic diagram of boundary conditions

圖5 模擬溶蝕深度與文獻[17]數據對比Fig.5 Comparison between modeling leaching depth and the experimental results from Ref.[17]

圖6為溶蝕24 h后各溶蝕特征量分布云圖，可以看出，固體鈣濃度由內及外逐漸降低，溶蝕程度逐漸增強，與文獻[17]試驗結果對比發(fā)現，固體鈣濃度分布與試驗所到溶蝕分布具有較好的一致性。孔隙液中鈣離子分布與固體鈣濃度分布規(guī)律基本一致。需要指出的是，孔隙率分布規(guī)律與固相鈣濃度分布規(guī)律正好相反，這是由于接觸溶蝕是一個由表及里的發(fā)展過程，與內部相比，外部固相鈣溶出量更大，孔隙體積增大程度更大。

圖6 模擬溶蝕特征量分布與文獻[17]試驗結果對比Fig.6 Comparison between modeling the spatial distributions of characteristic variables after 24 hours and the experimental results from Ref.[17]

選取接觸邊界節(jié)點進行分析，其溶蝕特征量隨時間變化關系如圖7所示。從圖7(a)可以看出，邊界節(jié)點隨著時間增加孔隙液鈣離子濃度逐漸減小并最終趨于穩(wěn)定；從圖7(b)可以看出，固相鈣含量變化呈現出相同的規(guī)律。這充分表明離子交換邊界條件能夠有效地避免采用“0濃度”邊界條件造成的問題，即接觸邊界節(jié)點孔隙液中鈣離子濃度始終為0，接觸邊界節(jié)點固相鈣溶出過程瞬間完成。

圖7 邊界節(jié)點孔隙液鈣離子濃度、固相鈣濃度 隨時間變化關系曲線Fig.7 Time-histories of calcium ion concentration in pore fluid and solid phase

通過數值模擬結果與試驗數據、有限差分法模擬結果對比可知，各項溶蝕特征量數據吻合良好，本文模型及程序的正確性得到驗證；采用離子交換邊界條件能夠有效避免“0濃度”邊界條件造成的問題。

圖8 試件5有限元模型Fig.8 FEM model of specimen 5

5 算 例

5.1 模型簡介

目前溶蝕試驗研究大多以一維、二維展開，三維較少涉及，因此，本文給出一個三維溶蝕的數值算例，研究三維溶蝕下的各溶蝕特征量的空間分布規(guī)律。考慮一個圓柱體水泥凈漿試件，尺寸為Φ20 mm×20 mm，將其置于6 mol/L的NH4NO3溶液中，使試件與溶液充分接觸；水泥化學組分及其他模型計算基本輸入參數取值見表1、表2；試件的有限元模型如圖8所示。

圖9 各截面溶蝕特征量模擬計算結果Fig.9 Modeling results of characteristic variables at each section

5.2 結果分析

為了研究溶蝕的空間分布規(guī)律，選取y=0 mm和z=10 mm 2個截面展開分析。經過40 d溶蝕后，各截面的各溶蝕特征量模擬計算結果云圖如圖9所示。從圖9可以看出，無論是徑向還是軸向，各溶蝕特征量的空間分布規(guī)律與上文算例呈現出一致的規(guī)律。

為了研究各溶蝕特征量隨時間的變化規(guī)律，在試件徑向取4個特征點，分別距離表層為1，2，3，5 mm，編號依次為①,②,③,④；在試件軸向，以試件頂面圓心為基點，從上至下取4個特征點，分別距離表層為1，2，3，5 mm，編號依次為⑤,⑥,⑦,⑧。根據計算可得，各溶蝕特征量與時間的變化關系，如圖10所示?？梢钥闯?，隨著時間的變化，特征點的各溶蝕特征量變化速率均呈現出由快變慢的總體規(guī)律，說明各處的溶蝕速率由快變慢；從圖10(b)可以看出，固相鈣濃度在7 000 mol/m3左右時，溶蝕速率開始變慢，這是由于此時Ca(OH)2溶解完成，C-S-H開始溶解，而C-S-H前期的溶解速率遠遠小于Ca(OH)2的溶解速率。無論是徑向還是軸向，各個特征點開始溶蝕的時刻均不相同，呈現出外部早于內部的總體規(guī)律，體現出接觸溶蝕是一個由表及里的過程。需要指出的是，在徑向與軸向上，距表層距離相同的特征點的溶蝕開始時間并不相同，徑向明顯早于軸向，據分析，這是由于所選取的徑向特征點所在處的濃度梯度大于軸向特征點，從而造成孔隙液中鈣離子徑向擴散速率大于軸向擴散速率。

圖10 各截面特征點孔隙液中鈣離子濃度、固相鈣濃度、 孔隙率隨時間變化的變化關系曲線Fig.10 Time-histories of calcium ion concentration in pore fluid and solid phase and porosity at feature points

6 結 論

(1)在隱式差分格式離散時間域的基礎上，用集中質量矩陣代替協(xié)調質量矩陣可以消除對時間步長選擇的限制，能夠有效避免數值振蕩現象。

(2)稀疏矩陣行壓縮存儲方案雖然具有內存消耗較小的特點，但在形成總體剛度矩陣時卻十分耗時，對大型有限元計算較為不利。

(3)本文建立了三維鈣溶蝕模型，通過對加速接觸溶蝕過程數值模擬，經過結果對比驗證，表明該模型能夠用于模擬三維狀態(tài)下發(fā)生的各維度的接觸溶蝕，具有一定的通用性。

(4)采用離子交換邊界條件，得到邊界濃度變化規(guī)律，解決了采用“0濃度邊界條件”造成接觸邊界各點瞬間完成溶出過程的問題。

(5)通過開展三維溶蝕數值試驗，結果表明：在三維狀態(tài)下，接觸溶蝕發(fā)展呈現出由表及里，各處溶蝕速率由快到慢，徑向溶蝕速率大于軸向溶蝕速率的總體特征。