考慮微觀結構影響的混凝土界面過渡區裂隙滲流-溶蝕耦合模型

王志良，張 躍，申林方，李 澤

(昆明理工大學建筑工程學院，云南，昆明 650500)

水泥漿與集料間的界面過渡區(ITZ)是混凝土材料的最薄弱環節[1]，具有孔隙率高、強度低、Ca(OH)2富集等特點，易成為裂隙的產生和擴展區。對于水工建筑物、隧道及橋梁等實際工程而言，混凝土結構長期與水接觸，其中的可溶性鈣與水中的Ca2+在濃度梯度的作用下不斷地溶解析出，導致混凝土結構孔隙率增大，物理力學特性降低。由于ITZ內的裂隙較為發育且Ca(OH)2含量相對較多，而Ca(OH)2的溶解性又遠大于其它水化產物，故ITZ內固相鈣的溶解速率較快[2?3]，溶蝕脫鈣現象也較為嚴重。因此，有必要建立混凝土ITZ裂隙的滲流-溶蝕耦合模型，并研究其固相鈣的溶解機制，為預測混凝土結構的服役性能提供理論依據。

針對水泥基材料的溶蝕問題，國內外學者開展了大量的研究工作，并取得了豐碩的成果。湯玉娟等[4]采用6 mol/L的氯化銨溶液對水泥砂漿內襯進行加速溶蝕試驗，并分析了其微觀孔隙結構的變化；方永浩等[5]采用室內試驗研究了含裂縫水泥基材料的滲透溶蝕過程，以及溶蝕前后的微觀結構變化；王曉梅等[6]研究了在水流作用下水泥基材料開裂表面的溶蝕特性，并建立了相應的Ca2+遷移和表面溶蝕模型；Segura等[7]基于Fick第二定律和未反應縮核模型對水泥砂漿的溶蝕脫鈣過程進行了數值模擬；Wan等[8]考慮固液平衡曲線的變化，建立了水泥漿體在6mol/L硝酸銨溶液中的加速溶蝕模型，并采用有限差分法對該模型進行了數值求解。然而，目前關于水泥基材料溶蝕特性的研究大多處于宏觀層面試驗以及溶蝕前后的微觀形貌觀測，極少涉及溶蝕過程中微觀結構的動態演變。另外，相比于普通的水泥基材料，ITZ的尺度較小，屬于微米級別，所以大部分學者主要對其形成原理、微觀結構特征及其對宏觀力學性質的影響展開研究，如：李冬等[9]建立了考慮混凝土微觀組分影響的宏觀力學性能理論預測模型，并指出混凝土的單軸抗拉強度受界面過渡區界面特性的影響顯著；Leemann等[10]基于微觀圖像分析，將靠近骨料表面并超過水泥漿體平均孔隙率15%的區域界定為ITZ；Kenny等[11?13]基于微觀試驗研究指出，ITZ的厚度分布在5μm～100μm；Barnes等[14]認為，隨著齡期的增長，ITZ內Ca(OH)2的排布狀態由平行于骨料發展為垂直于骨料表面。然而，ITZ的微觀試驗研究往往難以獲取ITZ裂隙溶蝕過程的動態演化，故數值模擬方法成為解決ITZ裂隙滲流-溶蝕問題的重要途徑。由于ITZ內部微觀結構的復雜性，傳統的宏觀數值計算方法(有限單元、有限差分等)難以考慮ITZ微觀形貌及復雜邊界條件的影響，研究Ca(OH)2溶解過程與流體滲流間的耦合作用機制。而基于分子動理論發展起來的格子Boltzmann方法，是一種宏觀離散、微觀連續的介觀模擬方法[15]，在許多傳統數值方法難以勝任的領域得到了廣泛應用，如多孔介質滲流、多物理化學場耦合作用等[15?18]。

鑒于此，本文在對混凝土ITZ進行微觀結構重構的基礎上，基于格子Boltzmann方法，采用雙分布函數分別模擬流體速度場和溶質濃度場的演化過程，建立了ITZ裂隙滲流-溶蝕耦合模型。根據兩個經典算例，分別驗證了計算模型在處理溶質對流-擴散及反應-擴散問題方面的有效性。最后，考慮混凝土ITZ微觀結構的影響，討論了不同初始流速、Ca(OH)2含量及Ca(OH)2排布狀態等因素對其裂隙滲流-溶蝕耦合作用機制的影響。

1 理論模型

1.1 裂隙滲流控制方程

假設ITZ裂隙中的流動處于層流狀態，且流體為不可壓縮的牛頓流體，則其流動過程滿足：

動量守恒方程：

式中：ρ為流體密度；u為滲流流速；t為時間；P為壓力；υ為流體的運動粘滯系數。

1.2 裂隙溶質傳輸方程

溶質在裂隙中傳輸主要由對流和擴散控制，則其傳輸過程滿足[19]：

式中：C為溶質濃度；D為溶質擴散系數。

1.3 裂隙表面溶解機制

ITZ的Ca(OH)2含量非常高，且溶解度遠大于其它水化產物，故其溶蝕過程主要由Ca(OH)2的溶解造成，其它水化產物的溶解對孔隙率的影響可忽略[20]。為此，本文將Ca(OH)2溶解作為ITZ裂隙溶蝕的主要原因，其溶解過程如下：

假定裂隙溶液中的Ca2+和OH?全部來自于Ca(OH)2的溶解，則OH?濃度C(OH?)為Ca2+濃度C(Ca2+)的2倍。根據一階動力學反應模型[6,21]，Ca(OH)2的表面溶解過程可表示為：

式中：Cw為固-液界面處的Ca2+濃度；n為由壁面指向流體的法線方向，如圖1所示；kd為溶解速率系數；Keq為化學平衡常數，可表示為：

圖1 固-液界面溶解反應示意圖Fig.1 Diagram of the dissolution reaction at the solid-liquid interface

式中，下標eq表示為離子飽和狀態。

在裂隙表面處，Ca(OH)2的溶解反應使其固相體積減小，由式(5)可得到節點體積演化：

式中：Vw為固-液界面處濃度所對應的無量綱體積；A為反應表面面積；M為Ca(OH)2的摩爾體積。

假定初始時刻，Ca(OH)2固體節點的體積為1，則溶蝕過程中其體積不斷減小。當固體節點體積減小為0時，將其更新為流體節點。

2 格子Boltzmann模型

2.1 速度場演化方程

采用二維九速(D2Q9)的單松弛格子Boltzmann方程描述流體流動，其演化方程為：

式中：fi(x,t)和fieq(x,t)分別為t時刻、x位置處沿i方向的速度分布函數和平衡態分布函數；?t為時間步長；τf為無量綱松弛時間；e i為沿i方向的離散速度，其9個方向上的離散速度為：

采用Chapman-Enskog展開，可以由式(8)推導出宏觀方程式(1)、式(2)，并得到宏觀物理量與格子參數間的聯系：

2.2 濃度場演化方程

Ca2+在裂隙中的對流-擴散過程與流速呈線性關系，只需較少的離散速度即可表征其運移規律[19]，故采用二維五速(D2Q5)模型模擬溶質的運移過程，其演化方程為：

同理，由式(12)可推導出到溶質的對流擴散方程式(3)，并得到宏觀參數與格子參數間的關系如下：

2.3 邊界處理

采用具有二階精度的半反彈格式，來模擬裂隙滲流過程中流體與固體壁面間的相互作用，其速度分布函數可表示為[15]：

式中：xf為與壁面相鄰的流體節點；i為流體節點指向壁面節點的方向；?i為i的反方向。

對于Ca(OH)2表面的溶解反應邊界，采用Zhang等[22]提出的半反彈處理格式，其濃度分布函數可表示為：

本文采用壓力驅動裂隙內的流體流動，其出、入口處的壓力邊界采用Guo等[23]提出的非平衡態外推格式，即將邊界節點xb的分布函數分為平衡態和非平衡態兩部分。平衡態部分fieq(xb,t)可由式(10)求得；而非平衡態部分則由與之相鄰的流體節點xbf代替，故其分布函數可表示為[23]：

3 ITZ微觀結構及裂隙生成

假定ITZ區域由可溶物(Ca(OH)2)、不溶物(其它水化產物)、孔隙構成。為了建立ITZ裂隙的微觀結構模型，在構造區域(L×H)中間設置一條開度為a的裂隙，將整個區域分為ITZ固相區和裂隙區，其示意圖如圖2所示。本文采用Wang等[24]提出的四參數隨機生長法構建其微觀結構，實施步驟如下：

圖2 ITZ裂隙示意圖Fig.2 Schematic diagram of ITZ fracture

1)設計ITZ的構造范圍，并劃分裂隙區及固相區。在固相區內將Pc作為Ca(OH)2初始固體節點的分布概率，隨機生成初始固相生長核，其分布概率Pc不得超過Ca(OH)2的體積分數。

2)以Ca(OH)2初始固相生長核作為中心，根據不同方向的生長概率Pi向相鄰的8個節點隨機生長，如圖3所示。

圖3 隨機生長方向示意圖Fig.3 Diagram of random growth direction

3)重復步驟2)的操作，直至Ca(OH)2的固體體積分數達到設計值。

4)在固相區的非Ca(OH)2范圍內，以一定的概率隨機生成不溶物的固相結構，直至體積孔隙率達到設定值。

5)在構造范圍內生成完可溶物(Ca(OH)2)、不溶物(其它水化產物)后，剩余的節點即為孔隙。

假定ITZ的構造范圍為L×H=200×50格子，在中心處設置一條裂隙，其開度為a=20個格子，固相區的孔隙率為?ITZ=0.15，Ca(OH)2體積含量為35%。采用上述方法，分別構建了Ca(OH)2水平生長、均勻生長和豎向生長的ITZ的微觀結構及相應裂隙。為了避免后續數值計算中，邊界處固相對流體流動的影響[19]，在裂隙兩端各增設了5個格子的孔隙，如圖4所示。圖中灰色物質為Ca(OH)2固體，黑色為不溶物，白色為孔隙。

圖4 Ca(OH)2不同生長方式的ITZ裂隙Fig.4 ITZ fractureswith different growth patterns of Ca(OH)2

4 模型驗證

4.1 半無限區域內溶質的對流-擴散問題

為了驗證本計算模型處理溶質對流-擴散問題的準確性，在半無限區域內建立了相應的計算模型。假定初始時刻計算域內的溶質濃度為0，流速恒定為u。隨后從左側注入濃度為C0的溶質，則一維的溶質對流-擴散方程可表示為：

針對L×H=200×10的計算網格，采用D2Q5模型進行了溶質對流-擴散的數值計算。為了便于分析，計算參數均采用格子單位。假定入口處的溶質濃度為C0=1.0，滲流流速為u=0.002，溶質的擴散系數為D=0.288。圖5為本文數值解與解析解的對比，從圖中可以看出，本文的數值解與解析解具有較好的一致性，這充分說明本計算模型處理溶質對流-擴散問題的有效性。

圖5 對流-擴散問題溶質濃度的數值解與解析解對比Fig.5 Comparison between numerical and analytic solutions of soluteconcentration for convection-diffusion problem

4.2 矩形區域內的反應-擴散問題

為了驗證本計算模型處理溶質反應-擴散問題的有效性，對于穩態下矩形區域內的反應-擴散問題進行了數值計算。假定計算域內的濃度分布由Laplace方程來描述，即：

計算模型如圖6所示，假定在上邊界處(y=H)發生一階化學反應；左邊界(x=0)注入恒定濃度為C0的溶質；右邊界(x=L)和下邊界(y=0)均為零通量邊界，其邊界條件如下：

圖6 穩態溶質反應-擴散問題示意圖Fig.6 Diagram of solutereaction-diffusion problem in steady state

在數值計算中，本算例的計算參數均采用格子單位。計算網格為L×H=200×50，假定初始時刻溶質濃度處于平衡狀態，C=Ceq=1.0，左側入口處的濃度C0=10.0，溶質的擴散系數為D=0.288。針對化學反應速率分別為kr=0.01、kr=0.05及kr=0.1 這3種情況，基于本文數值計算方法，求解了穩態溶質反應-擴散問題，其濃度場的本文數值解與解析解的對比，如圖7所示，圖中實線為解析解，虛線為本文數值解。由圖可知，本文數值解與解析解的吻合度非常高，兩者間的最大相對誤差僅為0.35%，這說明本文計算模型處理溶質反應-擴散問題的準確性。

圖7 反應-擴散問題溶質濃度的數值解與解析解對比Fig.7 Comparison between numerical and analytic solution of soluteconcentration for reaction-diffusion problem

5 分析討論

為了研究ITZ裂隙的滲流-溶蝕耦合作用機制，建立了如圖2(示意圖)及圖4(微觀結構圖)所示的計算模型。將ITZ裂隙的計算域L×H=1000μm×250μm 劃分200×50的網格。裂隙開度為a=100μm，模型中的計算參數取值[26]如表1所示。設在初始時刻裂隙內的Ca2+濃度處于飽和狀態，即C=Ceq，然后在左側入口處注入Ca2+濃度為0的軟水，同時在壓力驅動下裂隙發生滲流，并在固-液界面處發生Ca(OH)2的溶蝕現象。整個ITZ裂隙的初始孔隙率為?=0.515，隨著溶蝕的發生，其孔隙率逐漸增大。為了分析ITZ裂隙滲流-溶蝕間的耦合作用機制，采用滲透率表征其演化過程，相應的計算公式為：

表1 模型參數Table 1 Model parameters

式中：Q為裂隙面滲流量；dP為壓力差。

在初始時刻ITZ裂隙的滲流為Poiseuille流，其平均流速為：

5.1 滲流流速影響

為了研究滲流流速對裂隙滲流-溶蝕耦合作用機制的影響，針對Ca(OH)2體積含量為35%且均勻生長的ITZ裂隙，當初始平均流速uˉ0分別為4.167×10?6m/s、8.333×10?6m/s、1.250×10?5m/s時，進行了其滲流-溶蝕過程的數值計算。圖8為裂隙的整體孔隙率與相對平均流速的時程演化曲線。由圖可知，在初始階段由于軟水的突然注入導致Ca2+濃度顯著下降，使得壁面處的溶蝕速度加快。隨著時間的推移，裂隙內的Ca2+濃度逐漸趨于穩定，溶蝕速率逐漸減慢，其整體孔隙率和相對平均流速的變化趨勢也逐漸變緩。同時，對于初始平均流速較大的裂隙，由于其Ca2+濃度的傳輸速率較快，溶蝕過程迅速趨于穩定，故受隙寬影響的相對平均流速也表現出相似的趨勢?？傮w而言，滲流流速越快，裂隙的溶蝕速率越快，整體孔隙率和相對平均流速的變化趨勢也相應的越快。此外，圖9為溶蝕2.5×104s后裂隙中線處Ca2+濃度與滲流流速的分布圖。由圖可知，在裂隙中線位置處，Ca2+濃度隨著滲流流速的增加而減小，且在入口附近濃度變化較快，變化趨勢隨著向裂隙內部的推進而逐漸變緩。同時由于受到Ca(OH)2溶蝕裂隙壁面微觀結構變化的影響，滲流流速沿著裂隙長度方向并非均勻分布，而呈現出不規則的波動，同時流速的波動也對Ca2+濃度產生了一定的影響，在流速變化較大的位置處，Ca2+濃度的波動也較為劇烈。

圖8 整體孔隙率與相對平均流速的時程曲線Fig.8 Timehistory curvesof overall porosity and relative average flow velocities

圖9 t=2.50×104 s時裂隙中線處Ca2+濃度與滲流流速分布Fig.9 Distribution of Ca2+concentration and velocity at the central part of fracture at t=2.50×104 s

5.2 Ca(OH)2含量影響

為了分析Ca(OH)2含量對ITZ裂隙滲流-溶蝕特性的影響，在初始平均流速uˉ0為4.167×10?6m/s時，針對Ca(OH)2均勻生長且體積含量ω分別為35%、45%及50%這3種工況下的滲流-溶蝕過程進行了數值計算。圖10為不同Ca(OH)2含量情況下裂隙相對滲透率隨時間的變化趨勢，從圖中可以看出，在溶蝕初期，由于注入軟水使得裂隙內的Ca2+濃度顯著下降，溶蝕速率加快，相對滲透率的增長速率也較快。隨著時間的推移，裂隙相對滲透率的變化趨勢逐漸趨于穩定，且K/K0(ω=50%)>K/K0(ω=45%)>K/K0(ω=35%)。圖11為裂隙中截面(x=L/2)不同位置處的Ca2+濃度的時程曲線。由圖可知，距離壁面不同位置處Ca2+濃度的時程變化具有相似的變化規律，僅僅在壁面附近的濃度稍大。在初始階段，受注入軟水的影響ω=35%的ITZ裂隙中Ca2+濃度下降幅度最大、ω=45%次之、ω=50%最小，主要是由于Ca(OH)2含量較高的ITZ裂隙，其Ca(OH)2溶蝕速率較快，致使裂隙內部的Ca2+濃度相對較高。隨著裂隙滲流-溶蝕耦合作用的進行，3者的濃度變化趨勢漸趨穩定，且穩定后的Ca2+濃度變化規律為C(Ca2+)(ω=50%)>C(Ca2+)(ω=45%)>C(Ca2+)(ω=35%)?？傮w而言，3種Ca(OH)2含量情況下，裂隙溶蝕演化特征的差異主要是由于Ca(OH)2含量越高，壁面處其與流體的接觸面積越大，發生溶蝕的概率也相應增加，裂隙內Ca2+濃度分布也就越高，如果Ca2+得不到有效的運移，它又反過來降低Ca(OH)2的溶蝕速率。故裂隙的相對滲透率受控于Ca(OH)2含量與Ca2+濃度的綜合作用，即，裂隙內的Ca2+能夠得到及時運移，Ca(OH)2溶解加速，則ITZ內Ca(OH)2含量即使較高，也能夠快速溶解；反之，則會抑制Ca(OH)2溶解。

圖10 不同Ca(OH)2含量時相對滲透率的時程曲線Fig.10 Time history curvesof therelative permeability under different Ca(OH)2 contents

圖11 裂隙中截面不同位置處的Ca2+濃度時程曲線Fig.11 Time history curvesof Ca2+ concentration at different positionsof thecentral section of fracture

5.3 Ca(OH)2排布狀態影響

為了研究Ca(OH)2排布狀態對ITZ裂隙滲流-溶蝕特性的影響，針對如圖4所示的Ca(OH)23種排布狀態(水平生長、均勻生長和豎向生長)，在初始平均流速為uˉ0=8.333×10?6m/s時，進行了裂隙滲流-溶蝕過程的數值計算。圖12為3種Ca(OH)2排布狀態下ITZ裂隙整體孔隙率與相對滲透率的時程演化曲線。從圖中可以看出：在溶蝕前期，Ca(OH)2水平生長整體孔隙率變化相對較小，也就是溶蝕過程相對較慢。這主要是由于滲流初期水平生長的Ca(OH)2與流體接觸面積較大，致使其溶解釋放出大量的Ca2+(如圖13所示)，從而抑制Ca(OH)2后續的溶解反應。隨著裂隙滲流過程的進行，Ca(OH)2水平生長的裂隙，其整體孔隙率發展較快，并逐漸超越另外2種排布狀態。這主要是由于Ca(OH)2水平生長與裂隙滲流方向一致，溶解后其連通性較好，致使相對滲透率發展較快(如圖12(b)所示)，從而更利于Ca2+的運輸及溶蝕的發展。同時，裂隙的相對滲透率在溶蝕早期主要受Ca(OH)2溶解釋放出的Ca2+影響，待穩定后的3種排布狀態的規律為：水平生長最大，均勻生長次之，豎向生長最小。從滲流場的微觀流線圖(圖14)來看，Ca(OH)2水平生長裂隙的內部流線迂回度最小，而豎向生長則由于固相不溶物的影響，其迂回度最大，從而影響了裂隙整體的滲流流速和溶蝕進程。故Ca(OH)2微觀結構上的排布狀態，是影響ITZ裂隙滲流-溶蝕特性的一個重要因素。

圖12 不同Ca(OH)2排布狀態時整體孔隙率與相對滲透率的時程曲線Fig.12 Time history curves of overall porosity and relative permeability under different arrangements of Ca(OH)2

圖13 t=3.0×104 s時不同Ca(OH)2排布狀態裂隙中線處的Ca2+濃度分布Fig.13 Ca2+ concentration distribution at the central part of fracture under different arrangementsof Ca(OH)2 at t=3.0×104 s

圖14 溶蝕后ITZ裂隙局部滲流流線圖Fig.14 Local streamline diagram of ITZ fracture after dissolution

6 結論

基于格子Boltzmann方法，建立了考慮微觀結構影響ITZ裂隙滲流-溶蝕耦合演化過程的數值模型，并結合算例驗證了本數值模型在處理對流-擴散及反應-擴散問題方面的準確性。在此基礎上，討論了ITZ裂隙不同初始滲流流速、Ca(OH)2含量及Ca(OH)2排布狀態等因素，對其裂隙滲流-溶蝕耦合作用機制的影響，得到以下結論：

(1)ITZ裂隙的初始滲流流速越快，壁面的溶蝕速率越快，其整體孔隙率增大越快。而流體中的Ca2+濃度隨著滲流流速的增加而減小，且在入口附近濃度變化較快，其變化趨勢隨著向裂隙內部的推進而逐漸變緩。

(2)Ca(OH)2含量越高，在壁面處其與流體的接觸面積越大，溶蝕現象越易發生，導致裂隙內Ca2+濃度也相應的增加，這又反過來抑制Ca(OH)2的進一步溶蝕。故ITZ裂隙的滲流-溶蝕過程受控于Ca(OH)2含量與Ca2+濃度的綜合作用。

(3)對于Ca(OH)2不同排布狀態的ITZ裂隙，待溶蝕穩定后其相對滲透率為：水平生長最大，均勻生長次之，豎向生長最小。這主要是由于Ca(OH)2水平生長裂隙的內部流線迂回度最小，而豎向生長則由于固相不溶物的影響，則迂回度最大，從而影響了裂隙整體的滲流流速和溶蝕進程。