孫紅碩， 楊 格,2， 吳 斌

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

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

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

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

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

1 氯離子擴散解析模型

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

(1)

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

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

(2)

則有

(3)

(4)

(5)

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

(6)

或

(7)

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

(8)

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

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

(9)

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

(10)

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

P=U1e-ξ2

(11)

即

(12)

將式(12)進行積分，得

(13)

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

(14)

進而得到

Cs=U2

(15)

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

(16)

進而得到

(17)

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

(18)

進而得到

(19)

即

(20)

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

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

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

(21)

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

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

(22)

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

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

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

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

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

(23)

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

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

(24)

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

(25)

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

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

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

(26)

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

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

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

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

(27)

式中：F(t)為氯離子擴散系數(shù)中與時間相關的項；G(x)為與空間位置相關的項。

(28)

(29)

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

3 算例與分析

3.1 算例一

廣州某混凝土簡支梁橋?qū)嵭陌?，板?80 mm?？紤]到簡支梁橋板服役狀態(tài)，根據(jù)文獻[17]假定板跨中截面受拉區(qū)最外側(cè)拉應力水平為0.6，截面應力分布為線性分布，受拉區(qū)截面高度150 mm，僅考慮混凝土板跨中受拉側(cè)氯離子一維擴散，如圖1所示。參數(shù)取值：Cs=0.4%，C0=0，Dref=10-5.1mm2/s，m=0。以廣州氣候為例，廣州多年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]通過試驗模擬海洋環(huán)境，給出了當荷載類型為受拉時，拉應力對混凝土的氯離子擴散系數(shù)影響關系為：

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

(30)

式中：δ=σ/f，δ為應力水平，σ為混凝土截面所受應力，f為混凝土的強度。

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

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

2.93(0.6-0.004x)2]

(31)

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

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

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

圖2 氯離子濃度分布

本算例中均方根誤差計算公式為：

(32)

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

3.2 算例二

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

(33)

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

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

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

4 結 論

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

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