鐘新然， 戴　瑛， 陳　熹

(同濟大學(xué) 航空航天與力學(xué)學(xué)院，上海 200092)

?

混合微分變換法在混凝土溶質(zhì)運移問題中的應(yīng)用

鐘新然， 戴瑛， 陳熹

(同濟大學(xué) 航空航天與力學(xué)學(xué)院，上海 200092)

摘要：微分變換法是一種基于泰勒公式的求解方法，在求解包含間斷條件的溶質(zhì)運移偏微分方程(組)時會失效.采用將有限差分法與微分變換法相結(jié)合的求解思路，可使問題得到解決.利用該方法得到混凝土SO42-的反應(yīng)性溶質(zhì)運移方程以及CO2-SO42-和CO2-Cl-的全耦合溶質(zhì)運移方程組的近似解表達(dá)式并將結(jié)果與有限元結(jié)果進(jìn)行比較，證明了該方法的有效性.

關(guān)鍵詞：混凝土； 溶質(zhì)運移； 微分變換法； 有限差分法

混凝土是一種多孔介質(zhì)，當(dāng)它作為地下結(jié)構(gòu)的材料時，土壤或地下水中的CO2,SO42-,Cl-等化學(xué)物質(zhì)的侵蝕會導(dǎo)致混凝土力學(xué)性能劣化[1-4]，從而影響混凝土的承載力.因此，為了保證地下結(jié)構(gòu)的使用安全，有必要對混凝土中的溶質(zhì)運移過程進(jìn)行分析.

1997年，Gospodinov等[4]通過相關(guān)實驗標(biāo)定SO42-的擴散系數(shù)，采用不穩(wěn)定差分格式計算時間步長，獲得二維硫酸鹽溶質(zhì)運移方程的數(shù)值解.2001年，竇紅[5]將有限體積和有限元進(jìn)行結(jié)合，求解非線性對流擴散問題并證明其數(shù)值解的穩(wěn)定性.2009年，左曉寶等[6]對一維硫酸鹽侵蝕混凝土的化學(xué)-力學(xué)(M-C)耦合方程進(jìn)行Matlab自主編程，得到數(shù)值解.2012年，Bonakdar等[7]結(jié)合閉合形式的誤差函數(shù)及有限元對氯離子擴散方程進(jìn)行計算，分別獲得考慮反應(yīng)項為1階以及2階的近似解析解.

混凝土的溶質(zhì)運移過程通常是由偏微分方程(組)描述的[4-7]，而微分變換法(differential transform method，DTM)是近20年發(fā)展起來的一種求解偏微分方程的方法，它由趙家奎于1986年首次提出[8]，此后， Arikoglu等[9]、 ?zdemir等[10]、 Kangalgil等[11]和Shahed[12]等將其用于求解多種偏微分方程.在此，本文亦將其運用于求解混凝土中的溶質(zhì)運移問題.

微分變換法是一種基于泰勒級數(shù)的求解方法，但在處理某些含間斷條件的偏微分方程時，會出現(xiàn)級數(shù)的待定系數(shù)值不確定的問題，如混凝土溶質(zhì)運移問題中初始條件間斷時，級數(shù)的待定系數(shù)取值不唯一，從而使后續(xù)的求解無法進(jìn)行.為解決該問題，本文利用Chu等[13]的混合微分變換法(hybrid differential transform and finite difference method, hybrid DTM and FDM)的求解思路，將問題在空間上按差分格式離散，在時間域上仍采用DTM進(jìn)行求解，解決了間斷條件導(dǎo)致的系數(shù)取值不確定問題.采用這種求解方法，得到混凝土SO42-的反應(yīng)性溶質(zhì)運移以及CO2-SO42-和CO2-Cl-的全耦合溶質(zhì)運移方程的近似解表達(dá)式.

1微分變換法

微分變換法假設(shè)待求函數(shù)可展為泰勒級數(shù)形式，從而將原控制微分方程及初、邊界條件轉(zhuǎn)化為級數(shù)形式表示的方程，通過求解由各階級數(shù)的待求系數(shù)構(gòu)成的代數(shù)方程組，最終獲得級數(shù)形式的近似解[8].

以一維溶質(zhì)運移問題為例，設(shè)待求的化學(xué)物質(zhì)質(zhì)量濃度ρ(x,t)可以表示為[14]

(1)

式中：U(k,h)稱為ρ(x,t)的微分變換函數(shù)，與ρ(x,t)有如下關(guān)系：

(2)

待求的原函數(shù)及其微分變換函數(shù)具有加法、乘積、微分等性質(zhì)[14]，利用上述性質(zhì)，可將控制微分方程以及初始、邊界條件轉(zhuǎn)化為關(guān)于U(k,h)的遞推公式，求得U(k,h)的各項值并代入式(1)，即為ρ(x,t)的級數(shù)形式的近似解.

2含間斷條件的溶質(zhì)運移問題

一維溶質(zhì)運移問題的控制方程具有如下形式[15]：

(3)

式中：D為擴散系數(shù).

對于待求區(qū)域，在t=0時刻，x=0處，化學(xué)物質(zhì)即將由此開始滲入，在x>0處,化學(xué)物質(zhì)尚未到達(dá)，初始條件可寫為

(4)

初始端及末端邊界條件為

(5)

注意到式(4)的初始條件(含式(5)的邊界條件)在x=0處是間斷的.若用常規(guī)的DTM進(jìn)行求解，設(shè)待求化學(xué)物質(zhì)質(zhì)量濃度ρ(x,t)仍可表示為式(1)的級數(shù)形式，且與U(k,h)滿足式(2)的關(guān)系.將式(1)代入初始條件(4)，得到

(6)

而將式(1)代入邊界條件式(5)后，得到

(7)

顯然，當(dāng)ρ0≠0時，由式(6)與式(7)可知，U(k,h)在k=h=0處的取值是矛盾的，因此，對于這種間斷問題，常規(guī)的微分變換法是失效的.

3基于有限差分的微分變換法

為了解決上述初值點間斷現(xiàn)象，本文采用混合微分變換法對混凝土溶質(zhì)運移問題進(jìn)行求解，其求解步驟如下[13,16-17]：

考慮一般情況下的控制微分方程

(8)

初始條件為式(4)，邊界條件為式(5).

在空間沿x方向等分成n段，可得n+1個節(jié)點(含首末端).對域內(nèi)每一個節(jié)點，都有

(9)

設(shè)域內(nèi)兩相鄰節(jié)點的間距為Δx，則i點處的函數(shù)ρ(x,t)關(guān)于x的各微分項可由差分[18]獲得

i=1,2,…,n-1

(10)

差分后的控制微分方程化為

(11)

參照式(1),(2)的微分變換關(guān)系，在任意第i個節(jié)點上建立微分變換關(guān)系，有

(12)

(13)

將式(12)代入式(11)并結(jié)合微分變換的基本性質(zhì)，可得關(guān)于微分變換函數(shù)Ui(k)的遞推公式

(14)

初始條件和邊界條件相應(yīng)地轉(zhuǎn)化為

(15)

(16)

由式(15)和(16)可知，初始點間斷問題不再影響到系數(shù)的取值.利用式(15)和(16)，以及Ui(k)的遞推關(guān)系式(14)可以求得Ui(k)，將Ui(k)代入式(12)即可得到各節(jié)點處未知函數(shù)的近似表達(dá)式.

4數(shù)值算例

4.1SO42-反應(yīng)性運移問題

SO42-在混凝土中反應(yīng)性溶質(zhì)運移方程[7]為

(17)

式中：ρS為SO42-質(zhì)量濃度；D為擴散系數(shù)；R(ρS)為反應(yīng)項，當(dāng)R(ρS)僅取1階反應(yīng)項時[7]

(18)

式中：K為反應(yīng)速率.

初始條件:

ρS|t=0=ρS,ini,x=0;ρS|t=0=0,x∈(0,L]

(19)

邊界條件:

(20)

基本計算參數(shù)[7]為D=2×10-12m2·s-1，K=5×10-9s-1，L=0.05 m，ρS,ini=ρS,0=25 g·m-3.對域內(nèi)劃分節(jié)點見圖1.

圖1　節(jié)點劃分(單位：m)

由式(11)可知，式(17)在空間點上有

i=1,2,…,9

(21)

其中ρS,i=ρS,i(t)，為在域內(nèi)x=xi點的SO42-質(zhì)量濃度值，在邊界上，即i=0和10，ρS,i(t)的值可由邊界條件確定.

由式(12),(13)可知，ρS,i(t)及其微分變換函數(shù)Wi(k)可以表示為

(22)

根據(jù)式(21)，利用原函數(shù)與微分變換函數(shù)的基本性質(zhì)，可得微分變換函數(shù)Wi(k)滿足

(k+1)Wi(k+1)=

(23)

由初始條件式(19)和邊界條件式(20)，并結(jié)合式(23)，可計算得到各個節(jié)點的Wi(k)值，代入式(22)后，可得到各個節(jié)點位置處SO42-質(zhì)量濃度值ρS,i(t)的近似表達(dá)式.

分別取t=100，200，300 d，計算SO42-的侵蝕狀況，與有限元結(jié)果的對比見圖2.由圖2可知，近似解的結(jié)果與有限元的結(jié)果是一致的.此外，隨著時間的增長，SO42-逐步向遠(yuǎn)端侵蝕，對于一個給定位置，其化學(xué)物質(zhì)濃度也隨時間增加.

圖2　t=100,200,300 d時SO42-的侵蝕深度

由于近初始端的侵蝕情況嚴(yán)重，選取x=0.005 m，考察該處SO42-的質(zhì)量濃度在時間域上的變化情況，并與有限元結(jié)果進(jìn)行比較，見圖3.由圖3可見，在運移初期，近似解與有限元解存在偏差，在60 d以后兩者就幾乎完全重合了.運移初期的偏差主要是由運移時間所致，SO42-要經(jīng)過一段時間才能到達(dá)x=0.005 m處，在此之前，該點處的SO42-的質(zhì)量濃度為0.而用混合微分變換法進(jìn)行求解時，是通過初始條件和邊界條件，在所有離散點的整個時域t≥0上進(jìn)行的，得到的是級數(shù)解，對溶質(zhì)的“滯后”現(xiàn)象無法做到完全符合，不過這種偏差會隨離散點的增加而減小.

圖3　x=0.005 m處SO42-侵蝕深度的結(jié)果對比

4.2兩種溶質(zhì)的全耦合運移問題

兩種化學(xué)物質(zhì)作用下的耦合溶質(zhì)運移控制微分方程組[19]為

(24)

式中：ρⅠ和ρⅡ，DⅠ和DⅡ分別為兩種化學(xué)物質(zhì)的質(zhì)量濃度以及擴散系數(shù)；RⅠ為阻滯因子；v為對流速度；K為兩種化學(xué)物質(zhì)的反應(yīng)系數(shù).

初始條件:

(25)

邊界條件:

(26)

控制方程式(24)在空間各點上的形式為

(27)

(28)

(29)

由于算例中SO42-的擴散系數(shù)DⅠ較小，取略長的時間t=500，750，1 000 d，考察SO42-的擴散情況，并與有限元計算結(jié)果進(jìn)行比較，見圖4.圖5為x= 0.005 m處，SO42-的質(zhì)量濃度在時間域上的變化情況，及其與有限元計算結(jié)果的比較.

圖5　x=0.005 m處SO42-侵蝕深度的結(jié)果對比耦合)

Fig.5Concentration of SO42-at 0.005 m

圖6為CO2的擴散情況，由于其擴散系數(shù)DⅡ相對較大，時間分別取t=20，35，50 d.圖7為x= 0.005 m處，CO2的質(zhì)量濃度在時間域上的變化情況，及其與有限元結(jié)果的比較.

圖4和6顯示近似解與有限元的結(jié)果基本一致，且隨著時間的增長，SO42-和CO2的逐步向遠(yuǎn)端侵蝕；隨深度的增加，侵蝕速度很快減小.

圖6　t = 20, 35, 50 d時CO2的侵蝕深度耦合)

Fig.6Concentration of CO2ast=20, 35, 50 d

由圖5和7可知，x= 0.005 m處SO42-和CO2的質(zhì)量濃度隨時間而增加，與4.2節(jié)中的情況類似，近似解對溶質(zhì)的“滯后”現(xiàn)象無法做到完全符合，相對于CO2的擴散情況，由于SO42-的擴散系數(shù)較CO2低兩個數(shù)量級，精度相對較差，導(dǎo)致質(zhì)量濃度的計算結(jié)果偏差略大.

圖7　x=0.005 m處CO2侵蝕深度的結(jié)果對比耦合)

Fig.7Concentration of CO2at 0.005 m

取時間t=20，35，50 d，考察Cl-和CO2的擴散情況，同時考察x= 0.005 m處，Cl-和CO2的質(zhì)量濃度隨時間的變化情況，見圖8～11，有限元計算結(jié)果也繪于相應(yīng)的圖中.

圖8　t=20, 35, 50 d時Cl-的侵蝕深度

圖9　x=0.005 m處Cl-侵蝕深度的結(jié)果對比

圖10　t = 20, 35, 50 d時CO2的侵蝕深度

圖11　x=0.005 m處CO2侵蝕深度的結(jié)果對比

5結(jié)論

本文利用基于有限差分的微分變換法，解決了混凝土溶質(zhì)運移問題中，初邊值條件間斷導(dǎo)致的普通微分變換法求解失效的問題，并用該方法求解了混凝土中SO42-的反應(yīng)性溶質(zhì)運移以及CO2-SO42-和CO2-Cl-的全耦合溶質(zhì)運移問題，得到了問題的近似解，通過與有限元結(jié)果的比較，證實了該方法的有效性.

參考文獻(xiàn):

[1]Hope B B, Ip A K C. Chloride corrosion threshold in concrete[J]. ACI Materials Journal, 1987, 84(4): 306.

[2]Loo Y H, Chin M S, Tam C T,etal. A carbonation prediction model for accelerated carbonation testing of concrete[J]. Magazine of Concrete Research, 1994, 46(168): 191.

[3]Somuah S K, Boah J K, Leblanc P,etal. Effect of sulfate and carbonate ions on reinforcing steel corrosion, as evaluated using AC impedance spectroscopy[J]. ACI Materials Journal, 1991, 88(1): 49.

[4]Gospodinov P N, Kazandjiev R F, Partalin T A,etal. Diffusion of sulfate ions into cement stone regarding simultaneous chemical reactions and resulting effects[J]. Cement and Concrete Research, 1999, 29(10): 1591.

[5]竇紅. 對流擴散方程的一種顯式有限體積-有限元方法[J]. 應(yīng)用數(shù)學(xué)與計算數(shù)學(xué)學(xué)報, 2001, 15(2): 45.

DOU Hong. A explicit finite volume-finite element method for nonlinear convection-diffusion problems[J]. Communication on Applied Mathematics and Computation, 2001, 15(2): 45.

[6]左曉寶, 孫偉. 硫酸鹽侵蝕下的混凝土損傷破壞全過程[J]. 硅酸鹽學(xué)報, 2009, 37(7): 1063.

ZUO Xiaobao, SUN Wei. Full process analysis of damage and failure of concrete subjected to external sulfate attack[J].

Journal of the Chinese Ceramic Society, 2009, 37(7): 1063.

[7]Bonakdar A, Mobasher B, Chawla N. Diffusivity and micro-hardness of blended cement materials exposed to external sulfate attack[J]. Cement and Concrete Composites, 2012, 34(1): 76.

[8]趙佳奎. 微分變換及其在電路中的應(yīng)用[M]. 武漢：華中理工大學(xué)出版社, 1988.

ZHAO Jiakui. The application of differential transform method in the circuit problem[M]. Wuhan: Huazhong University of Science and Technology Press, 1988.

[9]Arikoglu A, Ozkol I. Solution of boundary value problems for integro-differential equations by using differential transform method[J]. Applied Mathematics and Computation, 2005, 168(2): 1145.

[10]?zdemir ?, Kaya M O. Flapwise bending vibration analysis of a rotating tapered cantilever Bernoulli—Euler beam by differential transform method[J]. Journal of Sound and Vibration, 2006, 289(1): 413.

[11]Kangalgil F, Ayaz F. Solution of linear and nonlinear heat equations by differential transform method[J]. Seluk Journal of Applied Mathematics, 2007, 8(1):75.

[12]Shahed M E. Positive solutions of boundary value problems fornth order ordinary differential equations[J]. Electronic Journal of Qualitative Theory of Differential Equations, 2008, 15(1):1.

[13]Chu H P, Chen C L. Hybrid differential transform and finite difference method to solve the nonlinear heat conduction problem[J]. Communications in Nonlinear Science and Numerical Simulation, 2008, 13(8): 1605.

[14]Chen K, Huei Ho S. Solving partial differential equations by two-dimensional differential transform method[J]. Applied Mathematics and Computation, 1999, 106(2): 171.

[15]余紅發(fā), 孫偉. 混凝土氯離子擴散理論模型[J]. 東南大學(xué)學(xué)報: 自然科學(xué)版, 2006, 36(2): 68.

YU Hongfa, SUN Wei. Model research on chlorine ion diffusion in concretes[J]. Journal of Southeast University: Nature Science, 2006, 36(2): 68.

[16]Beilina L, Johnsson C. Hybrid FEM/FDM method for an inverse scattering problem[C]∥Numerical Mathematics & Advanced Applications. Milan: Springer, 2003:545-556.

[18]張文生. 科學(xué)計算中的偏微分方程有限差分法[M]. 北京: 高等教育出版社, 2006.

ZHANG Wensheng. Finite difference methods for partial differential equations in science computation[M]. Beijing: Higher Education Press, 2006.

[19]Liang M T, Lin S M. Modeling the transport of multiple corrosive chemicals in concrete structures: synergetic effect study[J]. Cement and Concrete Research, 2003, 33(12): 1917.

Application of Hybrid DTM and FDM in Solute Transport Problem of Concrete

ZHONG Xinran, DAI Ying, CHEN Xi

(College of Aerospace Engineering and Applied Mechanics, Tongji University, Shanghai 200092， China)

Abstract：The differential transform method (DTM) ceases to be effective when solving the solute transport equation with discontinuous condition. In order to solve this problem, hybrid DTM and finite difference method(FDM) were adopted to solve the solute transport equation of SO42-, coupled CO2- SO42-and coupled CO2-Cl- of concrete, and the approximate solutions were obtained. The comparison of the results with those obtained by FEM shows that the proposed method is effective to obtain the approximate solution of solute transport problem.

Key words：concrete; solute transport; finite difference method(FDM); differential transform method(DTM)

收稿日期：2015-06-23

基金項目：國家“九七三”重點基礎(chǔ)研究發(fā)展計劃(2011CB013800)

通訊作者：戴瑛(1966—)，女，教授，博士生導(dǎo)師，工學(xué)博士，主要研究方向為計算力學(xué).E-mail:ydai@#edu.cn

中圖分類號：TU37

文獻(xiàn)標(biāo)志碼：A

第一作者： 鐘新然(1992—)，男，博士生，主要研究方向為混凝土多場耦合. E-mail:101839@#edu.cn