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

陳德鵬

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

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

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

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

1 多物理場基本模型

利用COMSOL Multiphysics 軟件對水泥基復合的濕熱傳輸和耦合變形行為進行模擬計算時，物理過程數學模型與解析－有限元結合解法中的濕熱耦合傳輸模型不同，有關解析－有限元結合解法的相關內容詳見文獻［11］.

1.1 濕擴散模型

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

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

1.2 熱傳輸模型

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

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

1.3 力學平衡方程

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

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

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

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

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

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

確定物理場模型后，可利用GUI 界面進行模型選擇和尺寸、參數定義等，然后進行網格劃分、定義邊界條件等，從而可運行COMSOL 進行計算.COMSOL 在求解過程中，首先把結構力學模型和通式微分方程組結合起來進行轉換，轉變為統一的通式形式的微分方程組，然后進行求解得到溫度場、濕度場和位移場，實現三場的耦合求解.利用COMSOL 分析多物理場耦合問題，可以避免一些松散耦合解法求解多物理場問題的誤差，實現了多物理場的過程耦合，得出的物理過程的數值解答更加準確［13］.COMSOL 的求解結果可以用等勢線、曲線、云圖及動畫等多種方式來表達.

3 算例分析

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

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

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

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

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

圖2 混凝土試件長度變化

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

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

4 工程應用實例分析

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

4.1 數值模擬分析建模

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

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

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

圖5 隧道內部年度月平均溫濕度

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

混凝土中初始溫度和濕度認為是均勻分布的，根據測定，隧道主體結構混凝土初始溫度為12 ℃，初始濕度為3.5%(對應于RH=95%)［8］.

4.2 模擬計算結果及分析

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

4.2.1 溫濕度分布

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

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

表2 數值模擬中所用的混凝土參數

圖6 混凝土內部溫濕度模擬計算結果

4.2.2 濕熱耦合變形

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

圖7 混凝土濕熱耦合變形的COMSOL 數值模擬結果

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

4.2.3 橫向伸縮縫位移

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

由圖8可知，混凝土變形的COMSOL 數值模擬結果比解析－有限元結合解法結果更加接近光纖光柵實測結果，變形發展規律與光纖光柵監測結果基本一致，較準確地反映了混凝土在真實結構中的變形發展情況.

圖8 光纖光柵監測結果與數值模擬結果比較

5 結論

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

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

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］陳德鵬.基于多孔介質濕熱傳輸理論的混凝土濕熱耦合變形數值模擬及應用［D］.南京:東南大學材料科學與工程學院，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，Fazio 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］盛金昌.多孔介質流－固－熱三場全耦合數學模型及數值模擬［J］.巖石力學與工程學報，2006，25(S1):3028－3033.

Sheng Jinchang.Fully coupled thermo－hydro－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］.建筑材料學報，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)