基于Lattice-Boltzmann方法的水泥基材料溶液傳輸過程模擬

李東遙, 單鈺涵, 韓宇棟, 丁小平, 侯東偉,3

(1.上海交通大學 土木工程系, 上海 200240; 2.中冶建筑研究總院有限公司, 北京 100088;3.上海市公共建筑和基礎(chǔ)設(shè)施數(shù)字化運維重點實驗室, 上海 200240)

混凝土結(jié)構(gòu)通常處于氣、水、土及其耦合的服役環(huán)境中.外界的侵蝕性介質(zhì),通過擴散、對流、滲透等輸運機制到達材料內(nèi)部,引發(fā)溶蝕、化學腐蝕及電化學作用等物理、化學過程,最終導致材料損傷、破壞,這是引發(fā)混凝土材料與結(jié)構(gòu)耐久性問題的主要原因[1-2].在海洋環(huán)境中,氯離子侵入鋼筋表面,與鋼筋鈍化膜結(jié)合并擴散,使鋼筋脫鈍[3].脫鈍后的鋼筋在氧氣和水的作用下發(fā)生電化學腐蝕,進而引起鋼筋銹蝕和結(jié)構(gòu)脹裂.鋼筋混凝土結(jié)構(gòu)的壽命一般以氯離子到達鋼筋表面并達到臨界濃度為判據(jù),因此氯離子隨海水在混凝土中的傳輸過程決定了鋼筋混凝土結(jié)構(gòu)的耐久性能.海水在混凝土中的傳輸,通常采用水分傳輸系數(shù)或氯離子擴散系數(shù)來統(tǒng)一描述[4-6].

混凝土中,水蒸氣和液態(tài)水的傳輸方式不同.水蒸氣在濃度梯度的作用下向混凝土內(nèi)部擴散,而液態(tài)水則在壓力梯度下發(fā)生滲透[7].對于理想狀態(tài)下單次干燥或濕潤過程,可通過水分擴散方程予以描述[8-11];對于氯離子傳輸過程,傳統(tǒng)方法是使用考慮其與孔隙表面的結(jié)合效應(yīng)建立的改進擴散方程[12]描述.但水分遷移系數(shù)和氯離子擴散系數(shù)都只是一個唯象的綜合性參數(shù),并不能反映其微觀作用機理;且現(xiàn)有研究還沒有反映物理本質(zhì)的、對溶液傳輸機制進行描述的方法.實質(zhì)上,海水是一種“溶液”：從微觀機理上,溶液中作為溶劑的水和作為溶質(zhì)的各類鹽離子雖然同時向混凝土內(nèi)部遷移,但是其遷移規(guī)律不同,應(yīng)被看作相互區(qū)別但又相互耦合的2個遷移過程[5].

為了考察溶液在水泥基材料中的傳輸過程,本文首先采用自生長模型建立了水泥基材料細觀結(jié)構(gòu)模型,然后基于格子-玻爾茲曼(lattice-Boltzmann)方法,計算模擬水泥基多孔介質(zhì)中水分的多相傳輸和離子的遷移、擴散過程,初步考察以氯鹽溶液為代表的“水分-離子”雙重遷移機制,并分析材料參數(shù)的影響規(guī)律.

1 Lattice-Boltzmann方法簡介

Lattice-Boltzmann方法(LBM)是近幾十年發(fā)展的一種在微觀尺度上求解流體力學方程的離散化方法.微觀流體可看作空間微元中一系列的流體微粒,比分子尺度大,比宏觀尺度小.在LBM中,計算區(qū)域被離散成標準化網(wǎng)格,時間被離散成時間步,而微粒的運動則采用Boltzmann方程或其離散形式描述.LBM模型通常包括格子(即離散速度集合)、平衡態(tài)分布函數(shù)及其演化方程3個部分.演化方程通常表示為：

fi(x+ciδt,t+δt)-fi(x,t)=Ωi(x,t)

(1)

式中：x為格點位置;ci為流體粒子離散速度集合,i為離散速度維數(shù);δt為離散時間步長;t為當前時間步;fi為以速度ci運動的速度分布函數(shù);Ωi為碰撞算子,表示微粒間的碰撞對速度分布函數(shù)的影響.

Ωi的表達式為：

(2)

LBM由分子動力學方法發(fā)展而來,其基本思想是基于粒子碰撞和Boltzmann統(tǒng)計來描述微觀物理過程,因此不依賴于宏觀連續(xù)性假設(shè),可以模擬不規(guī)則形狀介質(zhì)中微觀-細觀尺度上的復雜流動問題.針對具體的物理場景,其無量綱參數(shù)依照約定選取,而物理場景中各個物理機制的重要性和主次分別,則根據(jù)物理過程和材料特性,通過構(gòu)建相應(yīng)的演化模型來實現(xiàn).近年來,隨著LBM理論的發(fā)展,學者們建立了多維度、多相流問題的計算方法[13-17].基于此,本文采用多種LBM模型,考察水泥基材料中的溶液傳輸問題.

2 水泥基材料微觀結(jié)構(gòu)建模

2.1 微觀結(jié)構(gòu)特征

水泥熟料為形狀不規(guī)則,并根據(jù)用途和性能要求,具有一定的粒徑分布的顆粒材料.水泥顆粒是C3S、C2S、C3A、C4AF等主要成分構(gòu)成的固溶體,其中C3S占45%～60%.

在微觀尺度上的水化硅酸鈣(C-S-H)凝膠中,由于顆粒間水分通過擴散作用由外而內(nèi)發(fā)生水化過程,在水化顆粒內(nèi)外層生成的C-S-H凝膠中存在低密度和高密度組分,且在較大顆粒的內(nèi)部存在未水化部分[18-21].而根據(jù)Jennings等[19]的理論,未水化的C3S平均密度為3.21g/cm3,水化后的C-S-H凝膠密度從1.75g/cm3至2.20g/cm3不等,此外Ca(OH)2密度為2.24g/cm3.綜合來看,水化后的水泥基材料可看作密度由內(nèi)而外逐漸減小的固相顆粒混合而成的多孔復合材料.因此,不能使用一般的多孔材料微觀結(jié)構(gòu)生成方法,而應(yīng)該參考水泥水化過程,將水灰比和水化度作為輸入量,建立具有一定密度梯度和孔隙率的水泥微觀結(jié)構(gòu)模型.

2.2 自生長法建模

本文提出了一種新的水泥基多孔介質(zhì)細觀模型生成方法,即基于Matlab軟件,參考Wang等[22]提出的多孔介質(zhì)模型生成步驟,改進其中單個顆粒的生成過程為水泥顆粒水化過程.全部步驟(見圖1)如下：

(1)將粒徑分布曲線離散為各個粒徑對應(yīng)的顆粒數(shù)量(圖1(a)).

(2)以粒徑為限制條件,在隨機形狀的多邊形未水化內(nèi)核外,生成水泥水化產(chǎn)物(圖1(b)).

(3)水化產(chǎn)物的數(shù)量達到預設(shè)后,得到單個水泥顆粒水化后的模型(圖1(c)).

(4)將單個顆粒的水化水泥顆粒模型隨機投放至二維區(qū)域中(圖1(d)、(e)).

圖1 水泥微觀結(jié)構(gòu)模型生成步驟Fig.1 Cement microstructure model generation steps

下面使用數(shù)學描述圖1(b)中的單個水泥顆粒的水化過程,即水泥顆粒模型自生長公式：

(3)

式中：di為生長出的第i個微元的密度;fρ沒有統(tǒng)一的表達式形式,而是隨著模型參數(shù)變化,用于控制各等級粒徑的顆粒體現(xiàn)出的不同水化情況;rg為最終粒徑;rg,min、rg,max為所有粒徑的最大值和最小值;imax為待水化的所有微元的總數(shù);KW/C與水灰比成正相關(guān),這里直接取水灰比的值,KW/C=1.0的情況下不進行特殊處理,與Jennings等[19]的文獻數(shù)據(jù)一致,在KW/C<1.0時,近似認為水泥顆粒邊緣的最小密度隨水灰比線性減小,KW/C=0時,水泥顆粒為未水化狀態(tài);Kα與水化度成正相關(guān),這里直接取水化度的值,水化度越大,水化深度越深,在水化度為0時,水化深度為0.本模型中水化度統(tǒng)一取1.0,意為標準養(yǎng)護 28d 后的水泥水化度.

本方法旨在精細控制水泥顆粒的密度分布.式(3)引入三角函數(shù)項近似水泥密度的排布,可認為與Fick第二定律的雙向擴散模擬結(jié)果近似.此外,本方法通過fρ對不同等級顆粒的密度分布進行單獨設(shè)置.根據(jù)陳長久等[23]的研究成果,在28d標準養(yǎng)護后,水泥顆粒水化深度約為3.5μm,因此粒徑在7μm以下的水泥顆粒將被完全水化,但由于反應(yīng)物的濃度梯度仍然存在,這些顆粒還將持續(xù)進行生長和擴散.因此參考更大粒徑的顆粒密度分布,減小其中央最大密度,并確保最終模型整體密度與水泥石密度仍然一致;而對于粒徑在7μm 以上的顆粒,為統(tǒng)一模型最終整體密度與水泥石密度,需要增大大粒徑顆粒的水化深度,直到模型整體密度滿足要求.

自生長法生成的水泥基材料模型的密度(ρ)分布如圖2所示.考慮到投放各等級粒徑的顆粒數(shù)為正整數(shù)，估算模型邊界為1420μm×1420μm。不同預設(shè)水灰比下，均可得到如圖2所示的水泥基材料模型，從而進行擴散、滲流等問題的數(shù)值模擬。

圖2 水泥基材料模型的密度分布Fig.2 Meso-structure model of cement

2.3 模型參數(shù)計算

2.3.1水灰比驗證

設(shè)置式(3)中的控制參數(shù)KW/C,即可調(diào)整水泥細觀模型的預設(shè)水灰比.根據(jù)Pommersheim等[24-25]的研究,對于單個水泥顆粒,在原始邊界內(nèi)外的產(chǎn)物體積之比約為1.0∶2.2.因此,對模型中的微元密度di進行統(tǒng)計分析,得到計算水灰比：

(4)

式中：dun為未水化水泥顆粒的密度;nun為模型中未水化部分的微元個數(shù);nhy為模型中水化后的微元個數(shù).

基于水泥顆粒生長參數(shù),隨機生成100個水泥基材料二維模型,得到的預設(shè)水灰比和計算水灰比及其誤差,如表1所示.

表1 預設(shè)水灰比和計算水灰比及其誤差

由表1可以看出,預設(shè)水灰比與計算水灰比的誤差較小.因此在后文中使用0.40、0.45、0.50這3個預設(shè)水灰比,建立具有代表性的3個水泥模型,進行模擬計算.

2.3.2孔隙率統(tǒng)計

模型生成后,遍歷模型中的孔隙,從而得到孔隙率φ,其表達式為：

(5)

式中：nd=0為孔隙微元的個數(shù),即模型中密度為0的微元個數(shù);ntotal為模型整體的微元個數(shù).

使用Image J軟件對模型中的毛細孔進行統(tǒng)計,并使用分水嶺(watershed)算法對連通的孔隙進行分離,以孔隙的二維面積表征孔隙大小,并根據(jù)圓形面積公式反推其等效孔徑.各模型的堆積孔隙率與最可幾孔徑如表2所示.

表2 各模型的堆積孔隙率與最可幾孔徑

由表2可以看出,當預設(shè)水灰比為0.45時,模型孔徑顯著增大;當水預設(shè)灰比為0.50時,模型整體孔隙率φ顯著上升.其原因在于,預設(shè)水灰比從0.40上升到0.45,對于水泥模型的主要影響為水化度增加,顆粒增大,從而導致堆積孔相應(yīng)增大;預設(shè)水灰比從0.45上升到0.50,水泥顆粒水化度已達到最大,此時水泥顆粒本身變化不大,但多余的水分形成孔隙,導致模型整體孔隙率上升.

3 水分、離子與溶液傳輸過程的LBM實現(xiàn)

3.1 水分傳輸過程的LBM模型

氣、液兩相的物質(zhì)傳輸及相變過程,屬于單組分多相流模型.對于復雜邊界,應(yīng)使用Shan等[26]的偽勢模型進行多相模擬,其演化方程的形式為：

(6)

式中：τ為弛豫時間,與流體黏度有關(guān);ωi為離散速度i的權(quán)系數(shù);E為外力項,在邊界處即作為邊界條件.

平衡態(tài)分布函數(shù)的形式為：

(7)

式中：ρ為流體的宏觀密度;u為流體的宏觀速度;cs為格子聲速.

本研究采用的離散速度格式為D2Q9模型,即二維模型中具有9個離散速度：

(8)

(9)

流體的宏觀密度ρ和宏觀速度u確定如下：

(10)

(11)

式中：F(x,t)為流體粒子所受的作用力,主要包括流-流作用力Fc和流-固作用力Fads,其表達式為：

(12)

(13)

式中：Gc、Gads分別代表流-流作用系數(shù)和流-固作用系數(shù),決定了表面張力和固相的潤濕性;s(x,t)為標示函數(shù),在固體微元處取1,在流體微元處取0;Ψ(x,t)為序參量,表示x位置上t時刻流體的物質(zhì)狀態(tài),出于計算穩(wěn)定性的考慮,在本研究中取：

(14)

式中：ρ0為臨界密度,即氣相和液相的分界密度;Ψ0為常量.

流體格點上的壓力值P可由非理想狀態(tài)方程求得：

(15)

3.2 離子擴散過程的LBM模型

當研究鋼筋混凝土結(jié)構(gòu)在海工環(huán)境下的腐蝕過程時,通常采用氯離子擴散系數(shù)予以表征,這其實是假定了混凝土中的孔隙水處于飽和狀態(tài),且僅發(fā)生環(huán)境中的氯離子向內(nèi)部擴散的過程.因此,僅考慮氯離子的擴散作用,忽略其他作用,LBM中的平衡態(tài)分布函數(shù)形式如下：

feq=ωiρ

(16)

在模型的半固體格點處需要采用部分反彈格式.經(jīng)過比較,本研究決定采用Walsh等[27]改良的部分回彈方法.該方法在碰撞步執(zhí)行前,重定向了需要回彈的速度分量,提高了并行計算的效率,且保證了在單個微元處的質(zhì)量守恒,具有一定的優(yōu)越性.其碰撞步演化方程為：

(17)

式中：0≤ns≤1是一個連續(xù)的變量,代表微元內(nèi)傳輸介質(zhì)的平均相對密度.

在模型中,水泥的細觀結(jié)構(gòu)僅包含毛細孔,而凝膠孔的存在體現(xiàn)為單個微元的密度減小.假設(shè)溶液可進入固相網(wǎng)格,其飽和濃度與該網(wǎng)格孔隙度成正比.因此將平衡態(tài)分布函數(shù)形式修正為：

feq(x,t)=(1-ns)ωiρ(x,t)

(18)

3.3 溶液傳輸過程的LBM模型

由于離子僅能通過孔隙液進行擴散,因此通過雙分布函數(shù)的LBM模型進行計算,將水分傳輸與離子擴散同時進行.根據(jù)對流擴散方程(ADE),流場給定的條件下,擴散的平衡態(tài)分布函數(shù)可取為與濃度、對流速度均線性相關(guān)的形式：

(19)

4 結(jié)果與討論

4.1 水分傳輸過程

使用0.40、0.45、0.50這3個水灰比的模型進行計算.模型左側(cè)區(qū)域初始條件設(shè)置為液相密度,右側(cè)初始條件為氣相密度.為保證計算穩(wěn)定,設(shè)定密度比為20∶1.模型左側(cè)邊界與右側(cè)邊界設(shè)置為恒定壓力差以模擬滲透壓的作用,均采用Zou-He壓力邊界條件[28].模型上下側(cè)邊界則采用循環(huán)邊界條件.模型內(nèi)部的非流體微元統(tǒng)一設(shè)置為密度大于1.80g/cm3的微元,避免產(chǎn)生過多的封閉孔隙.不同水灰比的水泥基材料中水分傳輸過程如圖3所示.

由圖3可見,在水分傳輸初期,即1001個時間步左右，液相(圖中黑色部分)由于壓力差以及毛細作用,迅速充滿模型表面附近的孔隙.與此同時，可以看到在模型內(nèi)部出現(xiàn)局部凝結(jié)的水滴.這是氣相在蒸汽壓作用下向小毛細孔中擴散，隨后累積、附著在固相微元(圖中灰色部分)表面形成的，隨著傳輸進行到中后期，液相在壓力差作用下,優(yōu)先通過較大的孔隙通道,而后進入較小孔,并與包含凝結(jié)水的孔隙連通.模擬結(jié)果體現(xiàn)了液相與氣相相互轉(zhuǎn)化并共同傳輸?shù)倪^程.縱向比較不同水灰比的模擬結(jié)果發(fā)現(xiàn),在較高水灰比的模型中,水分傳輸過程發(fā)展較快.

圖3 不同水灰比的水泥基材料中水分傳輸過程Fig.3 Moisture transmission process in cementitious materials with different water-cement ratios

4.2 氯離子擴散過程

仍然選取上述3個模型進行計算.模型的初始條件為：除左側(cè)和右側(cè)設(shè)置氯離子濃度C為1.0mol/L外,其他位置均為0mol/L,且無外力作用;模型孔隙中已完成飽水過程,以計算氯離子在多孔介質(zhì)孔隙水中的擴散過程.模型的左側(cè)及右側(cè)均為恒定濃度,設(shè)置方法為借助通量(濃度×速度)守恒得到邊界處未知的離散速度分布.模型上側(cè)和下側(cè)為循環(huán)邊界條件.計算結(jié)果如圖4所示.

由圖4可見,在擴散初期,模型左側(cè)邊界附近聚集有一定量的氯離子;隨著擴散的進行,邊界處的氯離子通過孔隙通道持續(xù)向內(nèi)部擴散,濃度前鋒逐漸推進;擴散進行到后期,模型中大部分孔隙被氯離子溶液充滿,此時的主導過程為氯離子從右側(cè)溶出,與左側(cè)的補充達到動態(tài)的平衡,模型整體的氯離子濃度分布不再發(fā)生變化,系統(tǒng)形成穩(wěn)態(tài)擴散.縱向比較3個模型中的擴散過程發(fā)現(xiàn),在擴散的初始階段,水灰比大的模型擴散更快;隨著時間的進行,介質(zhì)中形成穩(wěn)定的濃度梯度場,不同模型中的濃度分布趨于一致,均接近線性分布.

圖4 不同水灰比模型的氯離子擴散過程Fig.4 Chloride ion diffusion process of different water-cement ratio models

4.3 溶液傳輸過程

仍選取上述3個模型進行模擬,右側(cè)不設(shè)置濃度邊界條件,以模擬氯離子與水分同步擴散.上下兩側(cè)設(shè)置為循環(huán)邊界條件,其余邊界條件相同.計算結(jié)果如圖5所示.

由圖5可見：在水分尚未傳輸?shù)侥Ｐ蛢?nèi)部時,氯離子的傳輸主要在邊界處進行,且迅速達到飽和狀態(tài),與水分傳輸幾乎保持同步;隨著時間的延長,水分在模型孔隙中已經(jīng)以蒸汽的形式先一步到達較大的孔隙,而氯離子只能隨液相水擴散,因此落后于水分遷移;參考圖3可知,在約10001個時間步后,水分已經(jīng)基本充滿模型內(nèi)部,而此時溶液中的氯離子仍然在擴散;80001個時間步時,水灰比為0.50的模型率先完成傳輸過程,達到穩(wěn)態(tài).

圖5 不同水灰比模型的溶液傳輸過程Fig.5 Solution transfer process of different water-cement ratio models

4.4 擴散系數(shù)計算

由于擴散系數(shù)的計算需要穩(wěn)定的擴散濃度梯度場,因此考察飽和介質(zhì)中的離子濃度場變化過程.以濃度分布的變化率小于0.1%為穩(wěn)態(tài)擴散的標志,得到穩(wěn)態(tài)擴散所需的時間約為500000個時間步長.取第520001個時間步長時沿擴散方向的氯離子擴散速度和氯離子濃度.根據(jù)氯離子的有效擴散系數(shù)公式,可以計算該模型的氯離子擴散系數(shù)(Deff)：

Deff=JL/Δn

(20)

式中：J為穩(wěn)態(tài)擴散中的氯離子通量;L為擴散距離;Δn為模型兩側(cè)的濃度差.

不同水灰比模型的氯離子擴散系數(shù)的計算結(jié)果如表3所示.

表3 不同水灰比模型的氯離子擴散系數(shù)

將文獻收集到的水泥凈漿的氯離子擴散系數(shù)[29-31],與本文的數(shù)值模擬結(jié)果進行對比,如圖6所示.

由圖6可見：盡管同樣是對水泥凈漿使用電通量法測定氯離子擴散系數(shù),但不同文獻的測定結(jié)果仍有一定的差別,這是由材料制備的離散性和電通量法本身的誤差造成的.通過對比可以發(fā)現(xiàn),LBM的模擬結(jié)果與試驗測得的氯離子擴散系數(shù)在同一數(shù)量級,與特定試驗的測量數(shù)據(jù)接近,同時隨水灰比的變化規(guī)律近似;在水灰比為0.45和0.50的情況下,模擬結(jié)果呈現(xiàn)出與實測結(jié)果非常接近的趨勢,但在水灰比為0.40時,與實測結(jié)果吻合稍差.原因在于,在水化度不變的前提下,水灰比低于一定的數(shù)值會導致毛細孔大幅減少,其作為模型滲透性的影響作用減弱.此時,水泥顆粒本身的致密性起到了決定作用.理想模型控制了水泥的水化度不變,在單個顆粒水化完成后再投放進模型,因此,低水灰比的理想模型保留了一定的滲透性.而在實際情況中,室內(nèi)試驗統(tǒng)一采用養(yǎng)護28d的試件.相對于高水灰比的試件,低水灰比的試件中,水泥顆粒水化程度較低,且在水化前就已經(jīng)處于更為致密的狀態(tài),因此抗?jié)B性能較好.

圖6 不同文獻中水泥凈漿的氯離子擴散系數(shù)Fig.6 Chloride ion diffusion coefficient in different references

5 結(jié)論

(1)水泥基材料水分傳輸?shù)哪M中,在蒸汽壓作用下，液相和氣相在材料內(nèi)部孔隙中傳輸，并相互連通,達到飽和狀態(tài).水灰比大的模型中,水分傳輸過程較快.

(2)水泥基材料氯離子擴散的模擬中,由于模型兩側(cè)濃度差恒定,氯離子濃度前鋒不斷推進,直到形成穩(wěn)定的濃度梯度場.水灰比大的模型中,這一過程發(fā)生得更早.

(3)在水分與離子共同遷移的溶液傳輸模擬中,水分的傳輸過程在時間和空間上都先于離子擴散過程，這是因為離子擴散作用需要在液相環(huán)境中才能進行.

(4)使用擴散系數(shù)描述氯離子的穩(wěn)態(tài)擴散,在水灰比為0.45、0.50時,模型與試驗結(jié)果相符,水灰比為0.40時模擬得到的擴散系數(shù)偏大,這是由于模型生成時預先統(tǒng)一規(guī)定了水化度,而沒有考慮低水灰比時水化產(chǎn)物的水化度偏低所致.