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

陳德鵬

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

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

混凝土材料自身的變形性能有其濕、熱變化本質(zhì)并存在耦合作用[1-2],對于研究混凝土結(jié)構(gòu)的變形開裂及耐久性至關(guān)重要．混凝土的體積變化難免會造成混凝土及其結(jié)構(gòu)的變形開裂,使混凝土的滲透性大幅增加,混凝土的碳化、化學(xué)侵蝕、鋼筋銹蝕等耐久性能劣化過程更容易進行,對混凝土耐久性極為不利[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èi)的溫濕度分布．用傳遞函數(shù)對耦合傳輸進行求解的方法則相對可行[10],但該解析法要求的條件較嚴格,傳遞函數(shù)的選擇較難,尤其是對于研究環(huán)境溫濕度條件變化情況下的混凝土內(nèi)部溫濕度分布及由此導(dǎo)致的混凝土變形,已很難利用傳遞函數(shù)進行求解．此外,主流的商業(yè)有限元軟件還難以直接計算出混凝土中的濕度應(yīng)力及變形．近年來出現(xiàn)的多物理場耦合模擬計算軟件,易于將不同物理場耦合起來進行模擬,已逐漸得到廣泛應(yīng)用．其中,基于偏微分方程對科學(xué)和工程問題進行建模和仿真計算的交互開發(fā)環(huán)境系統(tǒng)COMSOL Multiphysics已在多個領(lǐng)域發(fā)揮重要作用．

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

1 多物理場基本模型

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

1.1 濕擴散模型

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

(1)

式中,M為濕含量;Dm為濕擴散系數(shù),實際計算時,應(yīng)對濕擴散系數(shù)進行修正,即用考慮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ù)),則一維情況下熱傳輸模型可表示為

(2)

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

1.3 力學(xué)平衡方程

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

εall=εe+εt+εsh

(3)

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

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

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é)模型進行模擬計算．

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

3 算例分析

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

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

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

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

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

圖2 混凝土試件長度變化

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

從圖2和圖3可看出,混凝土試件長度和含濕量的數(shù)值模擬結(jié)果與試驗值符合很好,說明通過該方法對混凝土材料及結(jié)構(gòu)進行變形模擬分析是可行的．圖3中曲線的形狀可從以下2個方面進行解釋:① 受端部環(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è)墻-頂板的順序進行混凝土澆注．利用本文提出的基于COMSOL多物理場耦合數(shù)值模擬方法對隧道側(cè)墻混凝土的變形進行模擬計算,數(shù)值模擬結(jié)果與主體結(jié)構(gòu)側(cè)墻某伸縮縫處位移的光纖光柵監(jiān)測結(jié)果、解析-有限元結(jié)合解法計算結(jié)果進行比較．

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

根據(jù)混凝土材料及跳倉施工的特點,計算中采用如下假定:① 混凝土為多孔介質(zhì)且內(nèi)部均勻連續(xù),因此初始溫濕度分布均勻;② 鋼筋混凝土的比熱、導(dǎo)熱系數(shù)等熱物參數(shù)可利用改進的加權(quán)平均方法由混凝土、鋼筋自身參數(shù)值進行推定[15];③ 墻板內(nèi)不同長度上的溫濕度在距表面相同深度處的大小相同;④ 徑向混凝土板兩端自由伸縮;⑤ 因關(guān)注隧道服役期間的混凝土變形,故忽略水泥水化熱的熱變形效應(yīng);⑥ 模擬計算自成型后60 d始,根據(jù)同配比混凝土室內(nèi)試驗及現(xiàn)場測試，受水泥水化影響的混凝土變形平均值取為230×10-6．

根據(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)混凝土和外包防水層的綜合效果較好,可認為無水分交換,且由于隧道側(cè)墻混凝土在一定深度湖底的土層中溫度相對恒定,因此可假定為絕濕絕熱邊界．

混凝土中初始溫度和濕度認為是均勻分布的,根據(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é)合解法對某隧道混凝土墻板進行模擬計算，得到的內(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可見,距離隧道內(nèi)表面較遠時,混凝土內(nèi)部溫濕度變化對環(huán)境溫濕度變化不敏感．與環(huán)境溫濕度的波動相比,混凝土變形的變化存在滯后性,這符合導(dǎo)熱性差的混凝土的熱脹冷縮特點．

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

4.2.3 橫向伸縮縫位移

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

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

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

5 結(jié)論

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

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

)

[1]Shoukry S N,William G W,Downie B,et al.Effect of moisture and temperature on the mechanical properties of concrete[J].ConstructionandBuildingMaterials,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].ConcreteInternational,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].CementandConcreteResearch,2007,37(5): 759-766.

[5]De Borst R.Challenges in computational materials science: multiple scales,multi-physics and evolving discontinuities[J].ComputationalMaterialsScience,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]//TheFourthBiotConferenceonPoromechanics.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].InternationalJournalofHeatandMassTransfer,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].AppliedMathematicsandComputation,2013,219(13): 7275-7285.

[9]Jensen O M,Hansen P F.Autogenous deformation and RH-change in perspective[J].CementandConcreteResearch,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].InternationalJournalofHeatandMassTransfer,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].InternationalJournalofCivilEngineering,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-hydro-mechanical model of saturated porous media and numerical modelling[J].ChineseJournalofRockMechanicsandEngineering,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].JournalofBuildingMaterials,2007,10(2): 127-131.(in Chinese)