大體積混凝土濕熱耦合研究

高 鵬 李鵬輝 鮑克蒙

（1.新疆額爾齊斯河流域開發工程建設管理局，新疆 830000；2.清華大學 水沙科學與水利水電工程國家重點實驗室，北京 100084；2.清華大學 水利水電工程系，北京 100084）

隨著我國水利水電事業迅猛發展，我國將建設一大批大中型骨干水利水電工程，其中涉及大量的大壩、船閘、電站廠房等大體積混凝土結構.大體積水工混凝土結構，具有結構形式復雜、體積大、施工過程受環境影響大等的特點，其上常產生裂縫，影響結構安全.大體積混凝土裂縫的產生因素有結構承受外界如水壓、自重等結構荷載及溫度與濕度等環境變化［1].對于同一大體積水工混凝土結構，常常由于溫度與濕度的共同作用而引起混凝土結構產生由內而外或由外而內等不同型式的開裂現象.大體積混凝土結構所承受的結構外荷載已得到廣泛重視，若在結構設計中考慮充分，一般可以避免外荷載引起的結構性裂縫的產生.但對于混凝土結構產生的非結構性裂縫，人們更為關注的是溫度變形引起的溫度裂縫.對于這類裂縫常采用低熱水泥或摻加摻和料（如粉煤灰）以降低混凝土的絕熱溫升、采用預冷骨料或加冰拌和混凝土以控制混凝土的出機口溫度（如葛洲壩和三峽的7℃工程）、采用倉面噴霧與遮陽等增加壩面散熱與隔離太陽輻射、采用水管冷卻降低壩體內部溫度上升等溫控措施加以防止.但眾所周知，大多數混凝土表面裂縫的產生很大程度與混凝土表面的濕份喪失引起濕差應力（干縮應力）有關.實際上混凝土的溫度場和濕度場相互耦合會產生濕熱耦合效應：一方面混凝土濕份的遷移變化不但引起混凝土變形，同時也影響混凝土熱學性質的變化，從而影響溫度傳導過程和溫度應力；另一方面混凝土溫度場的變化也影響混凝土內濕份的遷移擴散，從而影響濕度分布和濕度應力.

對傳熱傳濕的研究，不同學者提出了不同的理論與模型［2－8]，這些理論總體上來看，可歸納為以下4種理論：①蒸發前沿模型理論：Sherwood等將材料內部分為濕區和干區兩部分，對于干燥過程，干區最初出現在材料的表面，并逐漸向內部的濕區推進.這一推進面稱為蒸發面，并且認為蒸發只在這一面上進行.在干區濕度主要通過水蒸氣的方式擴散，而在濕區濕度主要通過水滲透的形式擴散.②體積平均模型理論：Whitaker等人通過對控制點質量、動量和能量守恒方程采用容積平均化技術來建立宏觀控制方程，由于其考慮到了多孔介質內部濕份和能量的各種傳輸機制，因此，能較全面的反映多孔介質內部的傳熱傳質過程.但是，該模型的缺點是傳輸系數的確定極其困難.③非平衡熱力學模型理論：Luikov等［5－7]利用非平衡熱力學原理分析了多孔介質內部的傳熱傳質過程，建立了不穩定熱質輸運方程，并首次在方程中體現了傳熱對傳質過程的影響.該理論認為多孔介質的傳熱傳質，實際上是濃度場X、溫度場T和壓力場P相互作用、相互制約的過程.因此，利用能量守恒、質量守恒以及動量守恒原理便可以分別建立傳熱傳質的控制方程.該模型的缺點也表現在傳熱傳質參數的難以確定上.④蒸發凝結模型理論：Dayan等［8]認為蒸發并不僅僅發生在材料表面，在材料的內部也存在著水蒸氣的遷移，并伴隨有水蒸氣的蒸發和凝結.

1 混凝土濕熱耦合傳導基本方程

基于壩體混凝土材料的特點，本文提出以下基本假定：①將混凝土看作多孔介質，且在宏觀上是連續的、各向同性和均勻的.②在混凝土內部，存在固液氣三相狀態的水，且對傳導起主要作用的是液相和氣相形式的水.由于熱濕傳導過程相當緩慢，假定混凝土壩體內在局部始終處于固液氣三相熱動力學平衡狀態.③將水蒸氣看作理想氣體，可以利用理想氣體方程處理水蒸氣的各種參變量.④假定混凝土空隙內各處的水蒸氣壓力差是導致氣體水傳導的原因；毛細水壓力差是導致空隙內的毛細水傳導的原因.⑤假定壩體表面不存在液態形式的水散失，且只有氣體形式的水的散失（蒸發）.

以Luikov的傳熱傳質模型和Pedersen的傳熱傳濕模型［9－10]為基礎，結合混凝土傳熱傳濕的特點，根據能量守恒和質量守恒定律，建立壩體混凝土的濕熱耦合傳導控制方程如下：

式中，ρ為材料（如混凝土材料，基巖等）的密度；c為比熱容，單位J／（kg·K）；T為溫度；k為熱傳導系數（W／（m·K））；hlv為水的蒸發潛熱（J／kg）；δp為水蒸氣滲透性系數（kg／（m·s·Pa））；pv為水蒸氣壓強；λ為液態水（毛細水）的滲透性系數（kg／（m·s·Pa））；Pl為液態水（毛細水）壓強；QT為內熱源（如水化熱）；Qm為內濕源（如自干燥）.

如果將孔隙內水蒸氣的壓力看作溫度和含水量的函數，將孔隙內的毛細水吸收壓力近似看作含水量的單值函數，那么：

由pv＝pv（T，u），因此

由pl≈pl（u），因此

并且pv＝H（u）·pv，sat（T），代入式（1）和（2）整理可得：

濕熱耦合傳導方程存在三類邊界條件，它們的形式分別是：

第一類邊界條件：

第二類邊界條件：

第三類邊界條件：

結合混凝土壩體傳熱傳濕的特點，可以對這三類邊界條件加以細化：

1）基巖的側面，如果遠離壩體足夠遠，可近似為半無限問題，側面為對稱面，因此可用絕熱絕濕邊界處理.

2）基巖的底層，當取足夠深的巖體加以計算，認為其傳熱傳濕的過程已經完成，且受外界環境的影響可以忽略，因此可以看作絕熱絕濕邊界.

3）壩體和基巖的表面，其傳熱和傳濕邊界可以按情況不同看作第二或第三類邊界處理.

4）大壩蓄水后在上游壩面底部可看作第一類邊界條件.

2 有限元計算格式

利用Galerkin法推求傳熱傳濕方程的有限元格式［11－16]，熱傳導控制方程（5）對應的有限元求解格式為

濕傳導控制方程（6）對應的有限元格式為

如何求解上述兩個相互耦合的方程，以使求解過程既不會過于復雜又保證求解的精度，是壩體混凝土濕熱耦合有限元方程求解中的關鍵.本文采用了一種可被稱為交替變量的Galerkin法求解算法，通過交替變量將｛T｝和｛u｝由t時刻推求至t＋Δt時刻.其具體操作是，已知t0時刻的溫度值｛T｝0和濕度值｛u｝0，在時間由t推進到t＋Δt的一個步長內：

1）計算耦合方程（5）.將前一時間步計算得到的濕度值（或t0時刻的初始值）代入方程（5）中含｛u｝的項中；前一時間步計算得到的溫度值（或t0時刻的初始值）作為溫度的初值，這樣方程（5）化為僅含溫度變量的溫度擴散方程，在時間離散時僅計算Δt／2個步長，即取τ＝Δt／2，解方程可得到溫度向量｛T｝t＋Δt／2.

2）計算方程（6）.把溫度看作已知量，將第（1）步中計算得到的溫度結果｛T｝t＋Δt／2代入方程（6）中所有含｛T｝的項中.再用上一計算步計算得到的｛u｝t（或t0時刻的初始值）作為濕度初值，同樣計算半個步長Δt／2，求解可得到濕度解｛u｝t＋Δt／2.

3）計算方程（5）中溫度方程的后Δt／2個步長.將第（2）步中得到的｛u｝t＋Δt／2代入方程中含｛u｝的項中；以第（1）步中計算出的｛T｝＋Δt／2作為溫度的初值，同樣計算Δt／2步，可以得到t＋τ時刻的溫度值｛T｝t＋τ.

4）計算方程（6）中濕度方程的后Δt／2步.將第（3）步中得到的t＋τ時刻的溫度｛T｝t＋τ代入式中含｛T｝的項；并以第（2）步中計算得到的｛u｝t＋Δt／2作為初值，解方程可得到t＋τ時刻的濕度值｛u｝t＋τ.

重復第1步到第4步可以計算下一時間步的溫度值和濕度值.

此種解法的優點主要表現在：①既保證了方程原有的耦合特性，又使方程在局部求解中解耦.②有限元求解過程中，保持了總剛度矩陣的對稱性，從而可以獲得更快的求解速度和更小的內存消耗.③在計算中材料的屬性可以得到最快的更新，有利于解決材料參數的非線性問題.

3 算例分析

3.1 網格模型

算例采用如圖1所示的有限元計算網格.壩底高程120m，寬度98.4m；壩頂高程245m，寬度12m；下游壩坡1∶0.8；壩基計算區域取向上下游分別延伸100m，向地下延伸120m.

圖1 壩體網格圖

3.2 計算條件

圖2 壩頂網格及細部單元的位置

基巖為花崗巖，壩體混凝土和基巖的材料基本參數見表1.壩體和基巖的初始溫度都為30℃，環境溫度在計算過程中保持15℃；壩體的初始含水量為0.094 5kg／kg，基巖的初始含水量為0.004 395kg／kg，環境空氣的相對濕度在計算過程中恒為20%；考慮壩體和巖體表面的熱輻射作用，日光的日平均照射強度為67W／m2；平均風速取為5.0m／s，在當前風速下的對流換熱系數取為32W／（m2·s），不考慮毛細作用，不考慮混凝土的水化作用，考慮熱傳導與濕傳導的影響；時間步長為1d，計算60d.

表1 材料熱學參數

3.3 計算結果

計算得到60d后壩頂附近的溫濕度場如圖3～4所示.圖4中壩體內部濕度約為0.094 5kg／kg，等濕線間間隔0.015kg／kg，向外遞減.圖5～6分別為圖2中的細部處，距離壩面不同深度處的各點的溫度和濕度過程線圖.從以上的溫濕度場圖中不難看出，壩體內的溫度受到環境的影響較大，從初始時的30℃降低到20℃附近，而壩體內部的含水量仍然保持接近初始含水量0.094 5kg／kg，并且壩體的濕度梯度主要集中在壩面附近.

圖3 壩頂附近溫度場圖

從圖5的溫度過程線可以看出壩體表面的溫度迅速下降至環境溫度15℃附近，然而由于日光照射的影響，壩面達到吸放熱平衡時的溫度仍略高于15℃.從圖6的壩面附近的含水量變化過程線可以得到：距離壩面45cm處含水量的變化過程線已基本趨平，距離壩面20cm處的含水量的下降也已非常小，60d內僅下降了0.01kg／kg左右，然而壩面附近的含水量下降的較為劇烈，平均每10d就下降0.01kg／kg，這與溫濕度的實際規律是非常吻合的.

4 結 論

將混凝土看作多孔介質，基于Luikov傳熱傳質模型和Pedersen傳熱傳濕模型，結合水工大體積混凝土的特點，建立了以溫度和含水量為基本變量、含內熱源（如水化熱）和內濕源（如自干燥）的混凝土濕熱傳導耦合模型，采用交替變量的Galerkin法實現了混凝土濕熱耦合方程的有限元求解.該方法可用于的混凝土壩濕熱耦合仿真分析，為混凝土壩保溫保濕以防止混凝土壩開裂提供理論依據.

［1]杜 彬.混凝土大壩保溫保濕新技術［M].北京：中國水利水電出版社，2012.

［2]Whitaker S.Diffusion and Dispersion in Porous Media［J].AIChE J，1967，13：1066－1085.

［3]Whitaker S.Advances in the Theory of Fluid Motion in Porous Media［J].Ind.Eng.Chem，1969，61：14－28.

［4]Whitaker S.Simultaneous Heat，Mass and Momentum Transfer in Porous Media：A Theory of Drying［J].Advances in Heat Transfer，1977（13）：119－203.

［5]Luikov A V.Theory of Energy and Mass Transfer.Prentice－Hall［M].［s.l.]：Englewood Cliffs N.J.，1961.

［6]Luikov A V.Heat and Mass Transfer in Capillary－Porous Bodies［M].New York：Pergamon Press，1966.

［7]Luikov A V.System of Differential Equation of Heat and Mass Transfer in Capillary－porous Bodies［J].Heat Mass Transfer，1975，18：1－14.

［8]Dayan A.Self－similar Temperature，Pressure and Mois－ture Distributions Within an Intensely Heated Concrete Slab［J].Heat Mass Transfer，1985，25：1469－1476.

［9]Pedersen C R.A Transient Model For Analyzing The Hygrothermal Behavior of Building Constructions［M].Thermal Insulation Laboratory，Technical University of Denmark，1990：14－19.

［10]Rode C，Burch D M.Empirical Validation of a Transient Computer Model for Combined Heat and Moisture Transfer，Thermal Performance of the Exterior Envelopes of Building VI，December 4－8，1995，Clearwater Beach，FL，1995，283－295.

［11]陸金甫，關 治.偏微分方程數值解［M].2版.北京：清華大學出版社，2003.

［12]關 治，陸金甫.數值分析基礎［M].北京：高等教育出版社，1998.

［13]焦修剛.混凝土熱濕耦合傳導及結構抗裂仿真［D].北京：清華大學，2005.

［14]蘇向輝.多層多孔結構內熱濕耦合遷移特性研究［D].南京：南京航空航天大學，2003.

［15]虞維平.含濕多孔介質傳熱傳質的綜合理論及其實驗研究［D].北京：清華大學，1984.

［16]劉光廷，黃達海.混凝土濕熱傳導與濕熱擴散特性的實驗研究［J].三峽大學學報：自然科學版，2002，24（1）：15－23.