孫紅碩， 楊 格,2， 吳 斌

(1. 武漢理工大學(xué) a. 土木工程與建筑學(xué)院；b. 硅酸鹽建筑材料國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430070；2. 天津大學(xué) 中國地震局工程綜合模擬與城鄉(xiāng)抗震韌性重點(diǎn)實(shí)驗(yàn)室，天津 300350)

鋼筋混凝土結(jié)構(gòu)中氯離子的侵蝕是造成鋼筋銹蝕的主要原因[1]?；炷林新入x子的擴(kuò)散過程可通過Fick第二定律來描述，通過氯離子擴(kuò)散系數(shù)來反映混凝土中氯離子擴(kuò)散的快慢。當(dāng)混凝土所處的外部環(huán)境發(fā)生變化時，如溫度、相對濕度等因素的改變[2,3]，以及混凝土本身的水化程度變化等[2]，都會導(dǎo)致混凝土中氯離子擴(kuò)散系數(shù)隨時間變化。此外，同一混凝土結(jié)構(gòu)截面所受應(yīng)力不同[4,5]、混凝土結(jié)構(gòu)中的非貫通裂縫[6～9]等因素都會導(dǎo)致氯離子擴(kuò)散系數(shù)隨空間位置發(fā)生變化[10]。

Collepardi等[11]基于Fick第二定律計(jì)算了氯離子擴(kuò)散方程的解析解，但其模型的氯離子擴(kuò)散系數(shù)是定值，沒有考慮隨時間和空間的變化，計(jì)算結(jié)果與實(shí)際偏差較大。余紅發(fā)等[12]采用隨時間變化的氯離子擴(kuò)散系數(shù)，其氯離子擴(kuò)散系數(shù)模型考慮了混凝土齡期、溫度等因素的綜合作用，推導(dǎo)了氯離子擴(kuò)散方程的解析表達(dá)式。Mejlbro等[13]考慮了混凝土齡期、養(yǎng)護(hù)條件、環(huán)境、溫度等因素對氯離子擴(kuò)散系數(shù)的影響，同樣推導(dǎo)出了氯離子擴(kuò)散方程的解析表達(dá)式。上述研究雖然給出了時變氯離子擴(kuò)散系數(shù)下的擴(kuò)散方程解析解，但其采用的溫度、相對濕度等氣候參數(shù)為恒定值，不符合實(shí)際情況。Hassan等[3]為研究氣候變化對氯離子擴(kuò)散的影響，通過蒙特卡羅采樣獲得每個月的溫度和相對濕度值，并計(jì)算對應(yīng)月份的氯離子擴(kuò)散系數(shù)，然后采用有限差分法對混凝土中氯離子擴(kuò)散進(jìn)行模擬。但有限差分法在氯離子擴(kuò)散系數(shù)變化速率較快時計(jì)算量非常大，不便于工程應(yīng)用，而且現(xiàn)有氯離子擴(kuò)散解析模型中擴(kuò)散系數(shù)不能考慮溫度、相對濕度等因素隨時間變化的情況[2]。因此，需要建立任意氯離子擴(kuò)散系數(shù)時變形式的解析模型，以便于考慮因溫度、相對濕度等多種因素變化導(dǎo)致氯離子擴(kuò)散系數(shù)隨時間變化的情況。

在進(jìn)一步考慮混凝土結(jié)構(gòu)應(yīng)力的不均勻分布和非貫通裂縫等因素導(dǎo)致氯離子擴(kuò)散系數(shù)隨空間變化的情況時，李鏡培等[6]根據(jù)混凝土結(jié)構(gòu)內(nèi)裂縫的深度，將混凝土結(jié)構(gòu)分為有裂縫區(qū)和無裂縫區(qū)，建立了考慮混凝土結(jié)構(gòu)在非貫通裂縫條件下的氯離子擴(kuò)散模型，但該模型無法考慮混凝土內(nèi)氯離子擴(kuò)散系數(shù)時變的影響。劉建文等[4]考慮了混凝土結(jié)構(gòu)截面線性分布應(yīng)力導(dǎo)致的氯離子擴(kuò)散系數(shù)隨空間的變化，同時也考慮了受混凝土齡期影響的氯離子擴(kuò)散系數(shù)隨時間的變化，建立了基于Fick第二定律的混凝土氯離子侵蝕模型，但該模型沒有辦法考慮溫度、相對濕度等因素時變的影響。上述解析模型只適用于指定的氯離子擴(kuò)散系數(shù)表達(dá)式，不適用于氯離子擴(kuò)散系數(shù)隨空間變化的其他形式，且不能考慮溫度、相對濕度等因素時變的影響。Song等[14]為研究修復(fù)混凝土保護(hù)層對混凝土結(jié)構(gòu)耐久性的影響，通過采用有限差分法中的Crank等[15]差分格式，計(jì)算了擴(kuò)散系數(shù)隨時間和空間同時變化的氯離子擴(kuò)散過程。采用有限差分法能夠同時考慮擴(kuò)散系數(shù)隨時間和空間變化的情況，但其在氯離子擴(kuò)散系數(shù)變化速率較快時很難確定合理的時間步長，且計(jì)算量非常大，不便于工程應(yīng)用。然而，在氯離子擴(kuò)散系數(shù)隨時間和空間均發(fā)生變化時，又很難實(shí)現(xiàn)建立適用于任意變化形式的解析模型。因此，需要高效簡便的求解方法，以解決有限差分法在氯離子擴(kuò)散系數(shù)變化速率較快時計(jì)算存在的問題。

為解決現(xiàn)有解析模型中氯離子擴(kuò)散系數(shù)不能考慮溫度、相對濕度等因素隨時間變化的不足，本文通過理論推導(dǎo)建立了任意氯離子擴(kuò)散系數(shù)時變形式的解析模型，該解析模型不僅能考慮混凝土齡期、氯離子結(jié)合能力和溫度等因素的影響，而且還可以考慮溫度、相對濕度等因素隨時間變化的情況。為解決氯離子擴(kuò)散系數(shù)關(guān)于空間、時間均發(fā)生變化時，擴(kuò)散系數(shù)隨時間變化較快會導(dǎo)致采用有限差分法計(jì)算量非常大，且很難確定合理時間步長的問題，本文提出采用等效氯離子擴(kuò)散系數(shù)代替時變氯離子擴(kuò)散系數(shù)方法。該方法可以避免采用有限差分法因擴(kuò)散系數(shù)隨時間的變化關(guān)系而減小時間步長的缺點(diǎn)，以解決有限差分法在氯離子擴(kuò)散系數(shù)隨時間變化較快時計(jì)算存在的問題。

1 氯離子擴(kuò)散解析模型

當(dāng)氯離子擴(kuò)散系數(shù)為關(guān)于時間的函數(shù)時，基于Fick第二定律的氯離子擴(kuò)散方程可表示為：

(1)

式中：C(x,t)為氯離子濃度，簡寫為C；t為結(jié)構(gòu)暴露于環(huán)境中的時間；x為計(jì)算點(diǎn)距混凝土表面的深度；D(t)為任意時變形式的氯離子擴(kuò)散系數(shù)；Cs為混凝土表面氯離子濃度；C0為混凝土中的初始氯離子濃度。

為了將式(1)化為常微分方程，引入ξ(x,t)(記為ξ)，并令

(2)

則有

(3)

(4)

(5)

將式(3)(5)代入式(1)，得

(6)

或

(7)

式(7)中的自變量只有唯一的ξ，于是便可將其寫成如下常微分方程：

(8)

式(8)相應(yīng)的初始條件和邊界條件為：(a)ξ=∞時，C=C0；(b)ξ=0時，C=Cs。

現(xiàn)在求解式(8)，并使其滿足初始條件(a)和邊界條件(b)。令

(9)

將式(9)代入式(8)，得

(10)

將式(10)分離變量，積分得

P=U1e-ξ2

(11)

即

(12)

將式(12)進(jìn)行積分，得

(13)

式中：U1，U2為積分常數(shù)，可根據(jù)定解條件(a)(b)確定。將邊界條件(b)代入式(13)，得

(14)

進(jìn)而得到

Cs=U2

(15)

將初始條件(a)代入式(13)，得

(16)

進(jìn)而得到

(17)

將U1和U2代入式(13)，得

(18)

進(jìn)而得到

(19)

即

(20)

式中：erfc(x)為關(guān)于變量x的互補(bǔ)誤差函數(shù)。

式(20)即為本文提出的氯離子擴(kuò)散解析模型。該氯離子擴(kuò)散解析模型中氯離子擴(kuò)散系數(shù)可為任意時變形式，以便于考慮溫度和相對濕度等因素導(dǎo)致的氯離子擴(kuò)散系數(shù)隨時間變化的情況，實(shí)現(xiàn)對若干年后氯離子濃度分布的更精確預(yù)測。

為進(jìn)一步驗(yàn)證本文提出的氯離子擴(kuò)散解析模型的正確性，當(dāng)模型僅考慮混凝土抗氯離子擴(kuò)散的能力隨時間增長逐漸增強(qiáng)的因素時，其t時刻的氯離子擴(kuò)散系數(shù)可表示為[16]：

(21)

式中：Dref為參照氯離子擴(kuò)散系數(shù)，即混凝土齡期t0(一般取28 d)時的氯離子擴(kuò)散系數(shù)；m為擴(kuò)散系數(shù)衰減率[12]。

將式(21)代入式(20)得

(22)

式(22)與文獻(xiàn)[12]中解析解完全相同，驗(yàn)證了本文提出的氯離子擴(kuò)散解析模型的正確性。

在本文提出的氯離子擴(kuò)散解析模型中，氯離子擴(kuò)散系數(shù)雖然可為任意時變形式，但并不適用于氯離子擴(kuò)散系數(shù)隨空間變化的情況。因此，當(dāng)擴(kuò)散系數(shù)隨時間、空間均發(fā)生變化時，需要采用有限差分方法進(jìn)行求解。

2 基于Crank-Nicolson有限差分格式計(jì)算變擴(kuò)散系數(shù)的氯離子擴(kuò)散

2.1 變擴(kuò)散系數(shù)時的傳統(tǒng)有限差分法求解

當(dāng)擴(kuò)散系數(shù)隨時間、空間均發(fā)生變化時，基于Fick第二定律的氯離子擴(kuò)散方程可以寫成

(23)

該方程可采用Crank-Nicolson[15]有限差分格式進(jìn)行求解。當(dāng)距混凝土表面的最遠(yuǎn)計(jì)算距離為l，擴(kuò)散總時長為tN時，計(jì)算t時刻距混凝土表面x處的氯離子濃度分布的方法如下：

(1)首先對時間和空間進(jìn)行網(wǎng)格剖分，構(gòu)造差分格式。其中剖分是將區(qū)域Ω={0≤x≤l,0≤t≤tN}用兩簇平行直線

(24)

(2)采用有限差分法的Crank-Nicolson格式可表示為[14]：

(25)

根據(jù)式(25)逐步計(jì)算即可得到整個擴(kuò)散過程中的氯離子濃度分布，這一求解過程就是變擴(kuò)散系數(shù)時的傳統(tǒng)有限差分求解方法，后面簡稱為基于有限差分的傳統(tǒng)方法。其中，在氯離子擴(kuò)散系數(shù)為恒定值時，Crank-Nicolson格式的截?cái)嗾`差為o(τ2+h2)，因此需要限制步長以滿足計(jì)算精度。另外，在擴(kuò)散系數(shù)隨時間、空間發(fā)生變化時較大步長會引入新的誤差，往往需要更小的步長以減少因擴(kuò)散系數(shù)隨時間、空間發(fā)生變化引入的誤差，當(dāng)所取步長使得氯離子擴(kuò)散系數(shù)在該步長內(nèi)無變化可避免該誤差。

2.2 基于等效氯離子擴(kuò)散系數(shù)的有限差分法求解

通過式(20)可以看出，當(dāng)不考慮氯離子擴(kuò)散系數(shù)隨空間變化時，在基于Fick第二定律的氯離子擴(kuò)散模型中，隨時間變化的氯離子擴(kuò)散系數(shù)可由恒定的等效氯離子擴(kuò)散系數(shù)代替，即：

(26)

將恒定的等效氯離子擴(kuò)散系數(shù)Deq代替氯離子擴(kuò)散系數(shù)D(x)代入式(1)，其計(jì)算結(jié)果與本文所提氯離子擴(kuò)散解析模型計(jì)算結(jié)果一致，驗(yàn)證了上述方法的可行性。

前文介紹的有限差分法可以求解擴(kuò)散系數(shù)隨時間、空間均發(fā)生變化的情況，但通常需要較密的空間和時間網(wǎng)格劃分，導(dǎo)致計(jì)算量較大。為了減少時間上的網(wǎng)格數(shù)，假定氯離子擴(kuò)散系數(shù)隨空間變化時也可采用等效氯離子擴(kuò)散系數(shù)求解，本文擬采用等效氯離子擴(kuò)散系數(shù)代替前文的變擴(kuò)散系數(shù)。

假定氯離子擴(kuò)散系數(shù)可表示為：

D(t,x)=F(t)G(x)

(27)

式中：F(t)為氯離子擴(kuò)散系數(shù)中與時間相關(guān)的項(xiàng)；G(x)為與空間位置相關(guān)的項(xiàng)。

(28)

(29)

由于式(29)中的等效氯離子擴(kuò)散系數(shù)與時間t無關(guān)，可以避免采用有限差分法因擴(kuò)散系數(shù)隨時間的變化關(guān)系而減小時間步長的缺點(diǎn)，以解決有限差分法在氯離子擴(kuò)散系數(shù)變化速率較快時計(jì)算存在的問題。這一求解過程就是本文提出的基于等效氯離子擴(kuò)散系數(shù)的有限差分求解方法，后面則簡稱為基于有限差分的簡化方法。

3 算例與分析

3.1 算例一

廣州某混凝土簡支梁橋?qū)嵭陌?，板?80 mm?？紤]到簡支梁橋板服役狀態(tài)，根據(jù)文獻(xiàn)[17]假定板跨中截面受拉區(qū)最外側(cè)拉應(yīng)力水平為0.6，截面應(yīng)力分布為線性分布，受拉區(qū)截面高度150 mm，僅考慮混凝土板跨中受拉側(cè)氯離子一維擴(kuò)散，如圖1所示。參數(shù)取值：Cs=0.4%，C0=0，Dref=10-5.1mm2/s，m=0。以廣州氣候?yàn)槔?，廣州多年1～12月份平均氣溫為14.0，15.7，18.5，22.8，26.0，27.8，29.1，28.8，27.5，25.1，20.3，15.7 ℃；平均相對濕度為64.7%，73.3%，75.7%，77.4%，75.7%，79.0%，72.7%，74.3%，70.7%，62.9%，63.4%，61.4%[2]。

圖1 氯離子侵蝕示意/mm

任心波[18]通過試驗(yàn)?zāi)M海洋環(huán)境，給出了當(dāng)荷載類型為受拉時，拉應(yīng)力對混凝土的氯離子擴(kuò)散系數(shù)影響關(guān)系為：

D(δ)=D(0)(1-0.11δ+2.93δ2)

(30)

式中：δ=σ/f，δ為應(yīng)力水平，σ為混凝土截面所受應(yīng)力，f為混凝土的強(qiáng)度。

根據(jù)拉應(yīng)力對混凝土的氯離子擴(kuò)散系數(shù)影響關(guān)系可得混凝土彎曲受拉側(cè)擴(kuò)散系數(shù)隨深度變化函數(shù)為：

Dx=D0[1-0.11(0.6-0.004x)+

2.93(0.6-0.004x)2]

(31)

式中：x為距混凝土表面的距離；Dx為距混凝土表面x處的氯離子擴(kuò)散系數(shù)；D0為混凝土表面處的氯離子擴(kuò)散系數(shù)。則G(x)=[1-0.11(0.6-0.004x)+2.93(0.6-0.004x)2]。

根據(jù)本算例參數(shù)，并結(jié)合文獻(xiàn)[19,20]中溫度和相對濕度對氯離子擴(kuò)散系數(shù)的影響關(guān)系，計(jì)算出不同時間的氯離子擴(kuò)散系數(shù)值，假定該值為混凝土表面處對應(yīng)時間氯離子擴(kuò)散系數(shù)值D0(t)，則式(27)中F(t)=D0(t)。氯離子擴(kuò)散系數(shù)D(t,x)=F(t)G(x)。其中每個月按30 d計(jì)算，分別采用基于有限差分的簡化方法和基于有限差分的傳統(tǒng)方法計(jì)算第3年末的氯離子濃度分布。其中，時間步長τ=86400 s，空間長度l=150 mm，空間步長h=1 mm。

采用基于有限差分的傳統(tǒng)方法和基于有限差分的簡化方法分別計(jì)算得到第3年末的氯離子濃度分布見圖2。由式(25)可知，在每個計(jì)算步長內(nèi)氯離子擴(kuò)散系數(shù)無變化時，不會引入因氯離子擴(kuò)散系數(shù)變化產(chǎn)生的誤差，因此本算例中基于有限差分的傳統(tǒng)方法計(jì)算不會引入因氯離子擴(kuò)散系數(shù)隨時間變化產(chǎn)生的誤差。本算例以基于有限差分的傳統(tǒng)方法計(jì)算結(jié)果作為參考解，基于有限差分的簡化方法均方根誤差為0.0000066052%，驗(yàn)證了基于有限差分的簡化方法應(yīng)用于擴(kuò)散系數(shù)隨時間、空間均發(fā)生變化情況的可行性?；诖丝梢酝茢喑霎?dāng)擴(kuò)散系數(shù)隨時間、空間均發(fā)生變化時，該方法同樣可以避免采用有限差分法因擴(kuò)散系數(shù)隨時間的變化關(guān)系而減小時間步長的缺點(diǎn)。

圖2 氯離子濃度分布

本算例中均方根誤差計(jì)算公式為：

(32)

式中：M為總空間點(diǎn)數(shù)；C(i)為第i個空間點(diǎn)的氯離子濃度值，Cref(i)為參考解的第i個空間點(diǎn)的氯離子濃度。

3.2 算例二

為分析在氯離子擴(kuò)散系數(shù)隨時間變化較快的情況下，基于有限差分的傳統(tǒng)方法和基于有限差分的簡化方法的計(jì)算量和計(jì)算精度。參數(shù)取值：Cs=0.4%，C0=0，Dref=10-4.1mm2/s。根據(jù)參數(shù)取值并考慮氯離子擴(kuò)散系數(shù)衰減率m，可得第t時刻混凝土內(nèi)氯離子擴(kuò)散系數(shù)為：

(33)

通過第t時刻混凝土內(nèi)氯離子擴(kuò)散系數(shù)可以看出，氯離子擴(kuò)散系數(shù)衰減率越大氯離子擴(kuò)散系數(shù)隨時間變化越快。其中，基于有限差分的傳統(tǒng)方法，計(jì)算時間步長取2592000 s，m=0.64命名為組1；基于有限差分的傳統(tǒng)方法，計(jì)算時間步長取1 s，m=0.64命名為組2；基于有限差分的簡化方法，計(jì)算時間步長取2592000 s，m=0.64命名為組3；基于有限差分的傳統(tǒng)方法，計(jì)算時間步長取2592000 s，m=0.2命名為組4。所有組中空間步長都取1 mm，分別計(jì)算第3年末混凝土內(nèi)的氯離子濃度分布。以式(22)解析解為參考解，分別獲得所有組計(jì)算結(jié)果的均方根誤差(表1)。

表1 氯離子擴(kuò)散計(jì)算參數(shù)及氯離子濃度誤差

通過組1和組4計(jì)算所取參數(shù)和結(jié)果均方根誤差對比，相同時間步長下，氯離子擴(kuò)散系數(shù)隨時間變化快慢對基于有限差分法的傳統(tǒng)方法計(jì)算精度影響非常大。通過組1和組2計(jì)算所取參數(shù)和結(jié)果均方根誤差對比，當(dāng)氯離子擴(kuò)散系數(shù)隨時間變化相對較快時，基于有限差分的傳統(tǒng)方法通過減小時間步長雖然可提高計(jì)算精度，但是本算例中減少2592000倍的時間步長對精度提高并不顯著。其中，組2通過R 2017b版本的MATLAB計(jì)算時間已超過24 h，計(jì)算機(jī)處理器為Intel(R)Core(TM)i7-8700CPU@3.20GHz，如果在m=0.64情況下，要保證計(jì)算的均方根誤差達(dá)到5%以內(nèi)，至少需要將1 s的時間步長減小上萬倍，計(jì)算量非常大。因此在氯離子擴(kuò)散系數(shù)隨時間變化較快時，基于有限差分的傳統(tǒng)方法為保證計(jì)算精度需要采用非常小的時間步長，會導(dǎo)致計(jì)算量非常大。并且由于沒有氯離子擴(kuò)散系數(shù)隨時間變化速率與計(jì)算所需時間步長的精確關(guān)系，很難確定保證精度的合理時間步長。通過組1和組3計(jì)算所取參數(shù)和結(jié)果均方根誤差對比，在氯離子擴(kuò)散系數(shù)隨時間變化較快時，同樣步長下，基于有限差分的簡化方法計(jì)算精度比基于有限差分法的傳統(tǒng)方法計(jì)算精度高。這是由于氯離子擴(kuò)散系數(shù)隨時間變化越快，基于有限差分法的傳統(tǒng)方法為保證精度計(jì)算所需要的時間步長越小，本文所提方法將隨時間變化的氯離子擴(kuò)散系數(shù)轉(zhuǎn)化為恒定的氯離子擴(kuò)散系數(shù)值，使得本文所提方法在不減少計(jì)算步長的情況下仍能獲得較好的精度。驗(yàn)證了在氯離子擴(kuò)散系數(shù)隨時間變化速率較快時，本文所提方法一方面可避免基于有限差分法的傳統(tǒng)方法為保證計(jì)算精度需要采用極小時間步長的問題，另一方面可避免基于有限差分法的傳統(tǒng)方法為保證計(jì)算精度難以確定合理時間步長的問題。

4 結(jié) 論

(1)本文為解決現(xiàn)有解析模型中氯離子擴(kuò)散系數(shù)不能考慮溫度、相對濕度等因素隨時間變化的不足，通過理論推導(dǎo)建立了基于Fick第二定律的氯離子擴(kuò)散解析模型。該氯離子擴(kuò)散解析模型中氯離子擴(kuò)散系數(shù)可為任意時變形式，不僅能考慮混凝土齡期、氯離子結(jié)合能力和溫度等因素影響，而且還可以考慮溫度、相對濕度等因素隨時間變化的情況。通過與現(xiàn)有氯離子擴(kuò)散系數(shù)時變形式解析模型對比，驗(yàn)證了本文提出解析模型的正確性。

(2)當(dāng)考慮因混凝土結(jié)構(gòu)應(yīng)力的不均勻分布和非貫通裂縫等因素導(dǎo)致的氯離子擴(kuò)散系數(shù)同時隨空間變化的情況，即氯離子擴(kuò)散系數(shù)關(guān)于空間、時間均發(fā)生變化時，傳統(tǒng)有限差分法存在難以確定保證精度的合理時間步長以及計(jì)算量大的問題。本文提出采用等效氯離子擴(kuò)散系數(shù)代替時變的氯離子擴(kuò)散系數(shù)，在不減小時間步長的情況下仍能獲得較為滿意的精度。通過算例一和算例二分別驗(yàn)證了該方法在氯離子擴(kuò)散系數(shù)關(guān)于空間、時間均發(fā)生變化時的可行性，以及在計(jì)算氯離子擴(kuò)散方面的優(yōu)越性。