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

陳德鵬

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

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

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

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

1 多物理場基本模型

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

1.1 濕擴散模型

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

(1)

式中,M為濕含量;Dm為濕擴散系數(shù),實際計算時,應對濕擴散系數(shù)進行修正,即用考慮Knudsen擴散影響的混凝土濕擴散系數(shù)Dmk替換Dm,Dmk=Dmkf,kf為Knudsen擴散影響系數(shù)[2];Qm為濕源或濕匯．

1.2 熱傳輸模型

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

(2)

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

1.3 力學平衡方程

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

εall=εe+εt+εsh

(3)

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

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

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

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

確定物理場模型后,可利用GUI界面進行模型選擇和尺寸、參數(shù)定義等,然后進行網(wǎng)格劃分、定義邊界條件等,從而可運行COMSOL進行計算．COMSOL在求解過程中,首先把結(jié)構(gòu)力學模型和通式微分方程組結(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é)果,濕度分布還需要首先進行應力轉(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 工程應用實例分析

某淺埋湖底隧道采用放坡明挖法施工,大開挖后隧道主體結(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ù),因此初始溫濕度分布均勻;② 鋼筋混凝土的比熱、導熱系數(shù)等熱物參數(shù)可利用改進的加權(quán)平均方法由混凝土、鋼筋自身參數(shù)值進行推定[15];③ 墻板內(nèi)不同長度上的溫濕度在距表面相同深度處的大小相同;④ 徑向混凝土板兩端自由伸縮;⑤ 因關(guān)注隧道服役期間的混凝土變形,故忽略水泥水化熱的熱變形效應;⑥ 模擬計算自成型后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)生影響,應屬于第三類邊界條件．實際模擬分析中,環(huán)境溫度和濕度此時應屬于溫濕度函數(shù)在結(jié)構(gòu)邊界上已知的情況,見圖5．實際計算中,每天的溫濕度分布由實驗值通過插值得到．

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

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

混凝土中初始溫度和濕度認為是均勻分布的,根據(jù)測定,隧道主體結(jié)構(gòu)混凝土初始溫度為12 ℃,初始濕度為3.5%(對應于RH=95%)[8]．

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

在實際應用算例的模擬計算中,所用到的混凝土基本濕熱物參數(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時較大的濕度變化應與此時期內(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)境溫濕度的波動相比,混凝土變形的變化存在滯后性,這符合導熱性差的混凝土的熱脹冷縮特點．

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

4.2.3 橫向伸縮縫位移

為了與光纖光柵監(jiān)測結(jié)果(伸縮縫位移)進行比較,首先將COMSOL和A-FEM模擬的應變結(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) 本文中的算例及工程實際應用實例的結(jié)果表明,基于多物理場耦合利用COMSOL計算混凝土內(nèi)部溫濕度分布、應力及變形比利用解析-有限元結(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ù)值模擬及應用[D].南京: 東南大學材料科學與工程學院,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ù)學模型及數(shù)值模擬[J].巖石力學與工程學報,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].建筑材料學報,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)