基于多物理場耦合的混凝土濕熱變形數(shù)值模擬

陳德鵬

(安徽工業(yè)大學(xué)建筑工程學(xué)院，馬鞍山243002)

(江蘇省建筑科學(xué)研究院有限公司高性能土木工程材料國家重點實驗室，南京210008)

混凝土材料自身的變形性能有其濕、熱變化本質(zhì)并存在耦合作用［1－2］，對于研究混凝土結(jié)構(gòu)的變形開裂及耐久性至關(guān)重要.混凝土的體積變化難免會造成混凝土及其結(jié)構(gòu)的變形開裂，使混凝土的滲透性大幅增加，混凝土的碳化、化學(xué)侵蝕、鋼筋銹蝕等耐久性能劣化過程更容易進(jìn)行，對混凝土耐久性極為不利［3－4］.混凝土是一種典型的多孔介質(zhì)，將多孔介質(zhì)的傳熱傳質(zhì)原理應(yīng)用于混凝土濕熱耦合傳輸及變形數(shù)值模擬計算，已成為混凝土研究的主要熱點之一［5－8］.可用于混凝土濕熱耦合數(shù)值模擬分析的方法有以下4 種:①基于多孔介質(zhì)熱、質(zhì)遷移的唯象熱力學(xué)方法;②以Luikov 傳熱傳質(zhì)耦合方程為基礎(chǔ)的數(shù)值分析方法;③有限元分析方法;④多物理場耦合分析軟件.熱流唯象方程同時考慮了熱流驅(qū)動力與質(zhì)流驅(qū)動力的影響，而質(zhì)流唯象方程只考慮了質(zhì)流驅(qū)動力的作用，沒有涉及熱流驅(qū)動力對質(zhì)流作用的交叉項，因而，未能全面揭示熱、質(zhì)交叉效應(yīng).因受試驗條件和許多材料參數(shù)的影響，唯象學(xué)分析方法難以在實際中應(yīng)用［9］.以Luikov傳熱傳質(zhì)耦合方程為基礎(chǔ)的數(shù)值分析方法，主要用于多孔介質(zhì)的傳熱傳質(zhì)分析，但在求解溫濕度耦合過程中可能會出現(xiàn)如下問題:①由于存在很難計算的復(fù)雜特征值而得不到結(jié)果;②在控制方程和邊界條件同時耦合時無能無力;③準(zhǔn)確性欠佳，難以反映出多孔介質(zhì)內(nèi)的溫濕度分布.用傳遞函數(shù)對耦合傳輸進(jìn)行求解的方法則相對可行［10］，但該解析法要求的條件較嚴(yán)格，傳遞函數(shù)的選擇較難，尤其是對于研究環(huán)境溫濕度條件變化情況下的混凝土內(nèi)部溫濕度分布及由此導(dǎo)致的混凝土變形，已很難利用傳遞函數(shù)進(jìn)行求解.此外，主流的商業(yè)有限元軟件還難以直接計算出混凝土中的濕度應(yīng)力及變形.近年來出現(xiàn)的多物理場耦合模擬計算軟件，易于將不同物理場耦合起來進(jìn)行模擬，已逐漸得到廣泛應(yīng)用.其中，基于偏微分方程對科學(xué)和工程問題進(jìn)行建模和仿真計算的交互開發(fā)環(huán)境系統(tǒng)COMSOL Multiphysics 已在多個領(lǐng)域發(fā)揮重要作用.

本文利用COMSOL 軟件對混凝土濕熱耦合變形進(jìn)行數(shù)值模擬，并結(jié)合算例，與實際實驗結(jié)果和解析－有限元結(jié)合解法(A－FEM)結(jié)果進(jìn)行比較，以探討利用COMSOL 軟件進(jìn)行濕熱耦合甚至濕－熱－力耦合作用下水泥基復(fù)合材料內(nèi)部溫濕度場、應(yīng)力場及耦合變形行為的數(shù)值模擬的可行性和有效性.

1 多物理場基本模型

利用COMSOL Multiphysics 軟件對水泥基復(fù)合的濕熱傳輸和耦合變形行為進(jìn)行模擬計算時，物理過程數(shù)學(xué)模型與解析－有限元結(jié)合解法中的濕熱耦合傳輸模型不同，有關(guān)解析－有限元結(jié)合解法的相關(guān)內(nèi)容詳見文獻(xiàn)［11］.

1.1 濕擴散模型

通常，假定多孔水泥基材料中的濕傳輸為擴散形式，且濕度梯度是濕度(液態(tài)水和水蒸氣)的驅(qū)動力，根據(jù)Fick 定律和質(zhì)量守恒方程，對于一維傳輸有如下方程［2］:

式中，M 為濕含量;Dm為濕擴散系數(shù)，實際計算時，應(yīng)對濕擴散系數(shù)進(jìn)行修正，即用考慮Knudsen擴散影響的混凝土濕擴散系數(shù)Dmk替換Dm，Dmk=Dmkf，kf為Knudsen 擴散影響系數(shù)［2］;Qm為濕源或濕匯.

1.2 熱傳輸模型

通常采用Fourier 定律來描述水泥基材料中的熱傳輸，導(dǎo)熱系數(shù)則是已經(jīng)考慮了內(nèi)部對流等影響的表觀導(dǎo)熱系數(shù)(或稱名義導(dǎo)熱系數(shù))，則一維情況下熱傳輸模型可表示為

式中，Q 為熱源或熱匯，是由水泥水化放熱或其他非傳熱過程引起的熱量變化;T 為溫度;ρ 為水泥基材料表觀密度;cp為水泥基材料的比熱;λ 為水泥基材料的名義導(dǎo)熱系數(shù).

1.3 力學(xué)平衡方程

假定不同因素引起的混凝土的應(yīng)力均作用于材料整體，考慮到濕度孔隙壓力和熱應(yīng)力，本構(gòu)關(guān)系可以簡化為混凝土總應(yīng)變(εall)與不同因素引起的應(yīng)變之間的關(guān)系，即

式中，εe為外荷載引起的彈性應(yīng)變;εt為熱應(yīng)變;εsh為干縮應(yīng)變.

溫度變形可以根據(jù)溫度場利用混凝土材料的線性熱膨脹系數(shù)直接得到，而濕度分布引起的濕度變形很難如溫度的變形效應(yīng)一樣直接進(jìn)行處理.本文使用濕膨脹系數(shù)來解決濕度分布難以進(jìn)行變形計算的難題，并使得濕度變形計算像溫度變形計算一樣簡單.模仿熱膨脹系數(shù)，將濕膨脹系數(shù)定義為材料內(nèi)部單位濕度變化所導(dǎo)致的體積變化或長度變化.

從側(cè)吹浸沒燃燒熔池熔煉的長期性、穩(wěn)定性以及獲得經(jīng)濟指標(biāo)的先進(jìn)性等方面考察，其大規(guī)模工業(yè)化條件已經(jīng)成熟，再加上投資的節(jié)約性及配置的緊湊性，采用富氧側(cè)吹煤粉熔融還原工藝處理紅土鎳礦前景是樂觀的，該技術(shù)將是傳統(tǒng)紅土鎳礦工藝升級改造的首選技術(shù)。

2 基于COMSOL 的多物理場耦合模擬計算

瑞典COMSOL 公司開發(fā)的COMSOL Multiphysics 是一款大型的高級數(shù)值仿真軟件，因其高效的計算性能和杰出的多場直接耦合分析能力，被廣泛應(yīng)用于各個領(lǐng)域的科學(xué)研究以及工程計算［12］.在交互環(huán)境下，對于基于偏微分方程組的多物理耦合過程，COMSOL 不再需要編制復(fù)雜的偏微分方程組的求解器，而是利用其內(nèi)嵌的多種物理模型如化學(xué)反應(yīng)工程模型、熱傳導(dǎo)模型、結(jié)構(gòu)力學(xué)模型等.偏微分方程組模式是其功能最強大、最靈活的求解方法，有3 個數(shù)學(xué)應(yīng)用模式描述偏微分方程組:系數(shù)形式(coefficient form)、通式(general form)和弱形式(weak form).通常，通式模式適宜求解非線性問題，已足夠解決大部分多物理場耦合問題.本文圍繞混凝土中溫濕度場和濕熱耦合變形的濕－熱－力三場耦合問題，綜合采用通式微分方程組(描述熱傳導(dǎo)方程和濕擴散方程)和結(jié)構(gòu)力學(xué)模型進(jìn)行模擬計算.

確定物理場模型后，可利用GUI 界面進(jìn)行模型選擇和尺寸、參數(shù)定義等，然后進(jìn)行網(wǎng)格劃分、定義邊界條件等，從而可運行COMSOL 進(jìn)行計算.COMSOL 在求解過程中，首先把結(jié)構(gòu)力學(xué)模型和通式微分方程組結(jié)合起來進(jìn)行轉(zhuǎn)換，轉(zhuǎn)變?yōu)榻y(tǒng)一的通式形式的微分方程組，然后進(jìn)行求解得到溫度場、濕度場和位移場，實現(xiàn)三場的耦合求解.利用COMSOL 分析多物理場耦合問題，可以避免一些松散耦合解法求解多物理場問題的誤差，實現(xiàn)了多物理場的過程耦合，得出的物理過程的數(shù)值解答更加準(zhǔn)確［13］.COMSOL 的求解結(jié)果可以用等勢線、曲線、云圖及動畫等多種方式來表達(dá).

3 算例分析

試驗數(shù)據(jù)取自Hundt 開展的長達(dá)3年的混凝土變形試驗，6 個混凝土試件的尺寸均為2.4 m ×0.4 m ×0.4 m［14］;另外，利用解析－有限元結(jié)合解法對Hundt 試驗情況進(jìn)行數(shù)值模擬得到的結(jié)果，也作為對比數(shù)據(jù).在解析－有限元結(jié)合解法求解中，首先將濕熱耦合傳輸偏微分方程、邊界條件和初始條件進(jìn)行Laplace 變換，再利用傳遞函數(shù)將耦合方程組轉(zhuǎn)變?yōu)楸容^容易求解的四階偏微分方程，然后在Laplace 域中進(jìn)行求解，最后再對Laplace 域中求解結(jié)果進(jìn)行Laplace 逆變換，就可得到時間和空間上的溫濕度分布.根據(jù)溫濕度分布的解析結(jié)果，濕度分布還需要首先進(jìn)行應(yīng)力轉(zhuǎn)化過程，然后通過有限單元法計算混凝土在溫濕度耦合作用下的變形及規(guī)律［11］.

Hundt 試驗的試驗條件如圖1所示，試件側(cè)表面(不含端面)經(jīng)處理成為絕熱絕濕狀態(tài)保證其一維擴散，試件一端密封置于80 ℃的溫度環(huán)境中，另一端不作密封處理置于To=20 ℃、RH=45%的環(huán)境中.試件成型后對混凝土試件不同位置的溫度進(jìn)行了測定，并繼續(xù)對經(jīng)28 d 標(biāo)準(zhǔn)養(yǎng)護(hù)后的混凝土試件進(jìn)行濕度分布及變形測試.圖中，M0為試件內(nèi)部初始濕含量，T0為內(nèi)部初始溫度.

圖1 Hundt 試驗試件及邊界條件示意圖

基于多物理場耦合利用COMSOL 進(jìn)行模擬計算與利用解析－有限元結(jié)合解法進(jìn)行計算分析中需要用到的混凝土參數(shù)見表1.數(shù)值模擬計算結(jié)果與試驗結(jié)果見圖2和圖3.

表1 混凝土濕熱耦合變形計算所需參數(shù)

圖2 混凝土試件長度變化

圖3 混凝土內(nèi)部濕含量(550 d 時)

從圖2和圖3可看出，混凝土試件長度和含濕量的數(shù)值模擬結(jié)果與試驗值符合很好，說明通過該方法對混凝土材料及結(jié)構(gòu)進(jìn)行變形模擬分析是可行的.圖3中曲線的形狀可從以下2 個方面進(jìn)行解釋:①受端部環(huán)境濕度影響，混凝土內(nèi)部濕流向端面2 擴散;②由于附近高溫環(huán)境的影響，密封端濕流向端面2 方向遷移.因此，試件中部混凝土飽和后的濕含量比初始時還要高.

4 工程應(yīng)用實例分析

某淺埋湖底隧道采用放坡明挖法施工，大開挖后隧道主體結(jié)構(gòu)采用現(xiàn)澆混凝土跳倉法施工，每隔60 m 設(shè)置豎向變形縫，豎向變形縫間混凝土跳倉施工段內(nèi)設(shè)置3 個施工縫，形成長度為15 m 的4個施工段，施工段內(nèi)按底板－中隔墻－側(cè)墻－頂板的順序進(jìn)行混凝土澆注.利用本文提出的基于COMSOL 多物理場耦合數(shù)值模擬方法對隧道側(cè)墻混凝土的變形進(jìn)行模擬計算，數(shù)值模擬結(jié)果與主體結(jié)構(gòu)側(cè)墻某伸縮縫處位移的光纖光柵監(jiān)測結(jié)果、解析－有限元結(jié)合解法計算結(jié)果進(jìn)行比較.

4.1 數(shù)值模擬分析建模

根據(jù)隧道主體結(jié)構(gòu)施工時的豎向施工縫和分層施工高度，取15 m ×1 m ×3 m 的墻板作為分析對象，如圖4所示.

圖4 混凝土墻板尺寸示意圖(單位:m)

在隧道結(jié)構(gòu)內(nèi)側(cè)(x =0)，假定空氣與混凝土直接接觸(建設(shè)期內(nèi))，空氣中溫濕度將對混凝土內(nèi)部溫度和濕度產(chǎn)生影響，應(yīng)屬于第三類邊界條件.實際模擬分析中，環(huán)境溫度和濕度此時應(yīng)屬于溫濕度函數(shù)在結(jié)構(gòu)邊界上已知的情況，見圖5.實際計算中，每天的溫濕度分布由實驗值通過插值得到.

圖5 隧道內(nèi)部年度月平均溫濕度

在隧道外側(cè)(x=l)，自防水結(jié)構(gòu)混凝土和外包防水層的綜合效果較好，可認(rèn)為無水分交換，且由于隧道側(cè)墻混凝土在一定深度湖底的土層中溫度相對恒定，因此可假定為絕濕絕熱邊界.

混凝土中初始溫度和濕度認(rèn)為是均勻分布的，根據(jù)測定，隧道主體結(jié)構(gòu)混凝土初始溫度為12 ℃，初始濕度為3.5%(對應(yīng)于RH=95%)［8］.

4.2 模擬計算結(jié)果及分析

在實際應(yīng)用算例的模擬計算中，所用到的混凝土基本濕熱物參數(shù)與利用解析－有限單元結(jié)合解法相同［11］，由于采用了“線性濕膨脹系數(shù)”，原解析－有限元結(jié)合解法中的變形相關(guān)參數(shù)無需再使用，見表2.

4.2.1 溫濕度分布

利用基于多物理場耦合的混凝土濕熱耦合變形模擬方法和解析－有限元結(jié)合解法對某隧道混凝土墻板進(jìn)行模擬計算，得到的內(nèi)部溫濕度分布情況見圖6.圖中，“x/l =0.4”表示距離內(nèi)側(cè)表面相對深度為0.4 處的位置.

從圖6可看出，在隧道內(nèi)側(cè)為時變環(huán)境邊界條件、外側(cè)為絕熱絕濕條件的情況下，混凝土內(nèi)部溫濕度分布的變化滯后于實際環(huán)境溫濕度的變化，且隨距混凝土表面深度的不同而存在差異.距離內(nèi)表面越近，試件內(nèi)溫濕度的波動幅度也越大;而隧道外側(cè)混凝土內(nèi)的溫濕度變化波動較小.濕度和溫度變化曲線的發(fā)展趨勢相反.圖6(b)中，240 d 時較大的濕度變化應(yīng)與此時期內(nèi)混凝土內(nèi)部水分因環(huán)境溫度較高而蒸發(fā)比較活躍有關(guān).

表2 數(shù)值模擬中所用的混凝土參數(shù)

圖6 混凝土內(nèi)部溫濕度模擬計算結(jié)果

4.2.2 濕熱耦合變形

COMSOL 法和解析－有限元結(jié)合解法的濕熱耦合變形結(jié)果如圖7所示.

圖7 混凝土濕熱耦合變形的COMSOL 數(shù)值模擬結(jié)果

由圖7可見，距離隧道內(nèi)表面較遠(yuǎn)時，混凝土內(nèi)部溫濕度變化對環(huán)境溫濕度變化不敏感.與環(huán)境溫濕度的波動相比，混凝土變形的變化存在滯后性，這符合導(dǎo)熱性差的混凝土的熱脹冷縮特點.

4.2.3 橫向伸縮縫位移

為了與光纖光柵監(jiān)測結(jié)果(伸縮縫位移)進(jìn)行比較，首先將COMSOL 和A－FEM 模擬的應(yīng)變結(jié)果乘以伸縮縫間墻體長度得到伸縮縫間墻體變形值.

由圖8可知，混凝土變形的COMSOL 數(shù)值模擬結(jié)果比解析－有限元結(jié)合解法結(jié)果更加接近光纖光柵實測結(jié)果，變形發(fā)展規(guī)律與光纖光柵監(jiān)測結(jié)果基本一致，較準(zhǔn)確地反映了混凝土在真實結(jié)構(gòu)中的變形發(fā)展情況.

圖8 光纖光柵監(jiān)測結(jié)果與數(shù)值模擬結(jié)果比較

5 結(jié)論

1)根據(jù)混凝土的多孔介質(zhì)特點，基于多物理場耦合和多孔介質(zhì)濕熱傳輸原理提出了可利用COMSOL 進(jìn)行計算的濕－熱－力多物理場模型.經(jīng)算例和工程實例模擬計算證實了多物理場模型的正確性和COMSOL 軟件進(jìn)行混凝土耦合變形計算的可行性.

2)本文中的算例及工程實際應(yīng)用實例的結(jié)果表明，基于多物理場耦合利用COMSOL 計算混凝土內(nèi)部溫濕度分布、應(yīng)力及變形比利用解析－有限元結(jié)合解法進(jìn)行求解更加方便，模擬計算結(jié)果也更加準(zhǔn)確.

References)

［1］Shoukry S N，William G W，Downie B，et al.Effect of moisture and temperature on the mechanical properties of concrete［J］.Construction and Building Materials，2011，25(2):688－696.

［2］陳德鵬.基于多孔介質(zhì)濕熱傳輸理論的混凝土濕熱耦合變形數(shù)值模擬及應(yīng)用［D］.南京:東南大學(xué)材料科學(xué)與工程學(xué)院，2007.

［3］Mehta P K.Durability—critical issues for the future［J］.Concrete International，1997，19(7):27－32.

［4］Hossain K M A，Lachemi M.Strength，durability and micro－structural aspects of high performance volcanic ash concrete［J］.Cement and Concrete Research，2007，37(5):759－766.

［5］De Borst R.Challenges in computational materials science:multiple scales，multi－physics and evolving discontinuities［J］.Computational Materials Science，2008，43(1):1－15.

［6］Bazant Z P，Yu Q，Li G H，et al.Modeling of creep and hygrothermal deformations of concrete:intriguing consequences of nano－porosity［C］//The Fourth Biot Conference on Poromechanics.Columbia，USA，2009:3－16.

［7］Tariku F，Kumaran K，F(xiàn)azio P.Transient model for coupled heat，air and moisture transfer through multilayered porous media［J］.International Journal of Heat and Mass Transfer，2011，53(15/16):3035－3044.

［8］Sykora J，Sejnoha M，Sejnoha J.Homogenization of coupled heat and moisture transport in masonry structures including interfaces［J］.Applied Mathematics and Computation，2013，219(13):7275－7285.

［9］Jensen O M，Hansen P F.Autogenous deformation and RH－change in perspective［J］.Cement and Concrete Research，2001，31(12):1859－1865.

［10］Chang W J，Weng C I.Analytical solution to coupled heat and moisture diffusion transfer in porous materials［J］.International Journal of Heat and Mass Transfer，2000，43(19):3621－3632.

［11］Chen D，Qian C，Liu C.A numerical simulation approach to calculating hygrothermal deformation of concrete based on heat and moisture transfer in porous medium［J］.International Journal of Civil Engineering，2010，8(4):287－296.

［12］CnTech.COMSOL Multiphysics 全球第一款真正的多物理場耦合分析軟件［EB/OL］.(2008－09－05)［2012－02－01］.http://www.cntech.com.cn/newsflash/2008/comsol.html.

［13］盛金昌.多孔介質(zhì)流－固－熱三場全耦合數(shù)學(xué)模型及數(shù)值模擬［J］.巖石力學(xué)與工程學(xué)報，2006，25(S1):3028－3033.

Sheng Jinchang.Fully coupled thermo－h(huán)ydro－mechanical model of saturated porous media and numerical modelling［J］.Chinese Journal of Rock Mechanics and Engineering，2006，25(S1):3028－3033.(in Chinese)

［14］Hundt J.W?rme－und feuchtigkeit leitung in beton unter einwirkung eines temperaturf?lles，technical report 256［R］.Berlin:Deutscher Ausschu? für Stahlbeton，1975.

［15］陳德鵬，錢春香，王輝，等.水泥基材料比熱容測定及計算方法的研究［J］.建筑材料學(xué)報，2007，10(2):127－131.

Chen Depeng，Qian Chunxiang，Wang Hui，et al.Research on determination and calculation method of specific heat capacity of cement－based materials ［J］.Journal of Building Materials，2007，10(2):127－131.(in Chinese)