卞曉琳，何 平，吳青柏，施燁輝

（1.清華大學(xué) 環(huán)境科學(xué)與工程系，北京 100084；2.北京交通大學(xué) 土建學(xué)院隧道及地下工程教育部工程研究中心，北京 100044；3.中國科學(xué)院寒區(qū)旱區(qū)環(huán)境與工程研究所 凍土工程國家重點實驗室，蘭州 730000）

1 引 言

青藏鐵路建設(shè)中廣泛采用了拋石氣冷路基、拋石護(hù)坡路基等工程措施，該類結(jié)構(gòu)能夠自然冷卻路基下的凍土層，從而達(dá)到減少甚至消除寒區(qū)路基融沉的效果[1－2]。拋石可以看作多孔介質(zhì)，其內(nèi)部流體的對流換熱為多孔介質(zhì)的傳熱傳質(zhì)問題。多孔介質(zhì)中的熱質(zhì)遷移問題是一個較為復(fù)雜的強(qiáng)非線性問題，無法求得解析解，通常采用數(shù)值方法進(jìn)行理論分析，如Neto等[3]利用Darcy流，通過數(shù)值計算，對拋石路堤的對流傳熱進(jìn)行了研究；Ahmed等[4]應(yīng)用顯式差分法，分析了不可壓縮流體通過多孔介質(zhì)的強(qiáng)迫對流效果和相應(yīng)的熱質(zhì)運(yùn)輸過程；米隆等[5]基于多孔介質(zhì)流體動力性理論，利用有限元方法，對青藏鐵路拋石路堤、拋石護(hù)坡降溫效果進(jìn)行了數(shù)值計算分析；賴遠(yuǎn)明等[6]根據(jù)多孔介質(zhì)中流體熱對流的連續(xù)性方程、動量方程和能量方程，引進(jìn)流函數(shù)和應(yīng)用伽遼金法，導(dǎo)出了多孔介質(zhì)對流換熱的有限元公式，并對普通道碴路基、拋石護(hù)坡路基和拋石路基的溫度場進(jìn)行了分析比較，結(jié)果表明普通道碴路基的修建將使路基下凍土的溫度升高，使凍土路基出現(xiàn)熱不穩(wěn)定，而拋石路基和拋石護(hù)坡路基均能降低凍土路基的溫度；姜凡等[7]應(yīng)用有限體積法，對某塊石路基內(nèi)的溫度和速度特性進(jìn)行數(shù)值分析，計算時對塊石層分別采用多孔介質(zhì)模型和“塊石”模型，結(jié)果表明冬季塊石層內(nèi)空氣速度場的形態(tài)為一個逆時針的渦包，夏季則為順時針渦包，采用“塊石”模型計算得到路基內(nèi)的整體溫度略低于多孔介質(zhì)模型的結(jié)果，速度則大于多孔介質(zhì)模型計算結(jié)果。

目前關(guān)于拋石護(hù)坡路基自然對流方面的數(shù)值模擬研究較多，而強(qiáng)通風(fēng)條件下拋石護(hù)坡路基內(nèi)強(qiáng)迫對流特性方面的研究則較少。為此，本文基于多孔介質(zhì)中流體熱對流的連續(xù)性方程、非達(dá)西流動量方程和能量方程，對強(qiáng)通風(fēng)條件下青藏鐵路典型拋石護(hù)坡路基內(nèi)溫度場和流速場的分布形態(tài)進(jìn)行理論分析，結(jié)合現(xiàn)場試驗結(jié)果，對拋石護(hù)坡路基的強(qiáng)迫對流特性進(jìn)行研究。

2 理論背景

將拋石層視為多孔介質(zhì)，忽略路堤長度方向的不均勻效應(yīng)，簡化為橫斷面方向的二維各向同性問題?；谝陨霞僭O(shè)，其內(nèi)部熱對流方式為非穩(wěn)態(tài)的非等溫滲流，控制方程為連續(xù)性方程、動量方程和能量方程。

連續(xù)性方程：

式中：vx、vy分別為空氣在 x、y方向上的速度分量。

動量方程：

考慮空氣不可壓縮且ρa(bǔ)是溫度的函數(shù)，為了簡化分析，使用Boussinesq近似，即只有重力項中的空氣密度是可變的，可表示為

式中：β為空氣的熱膨脹系數(shù)；ρ0、T0分別為空氣密度和穩(wěn)定的參考值。

能量方程：

式中：Ca為空氣的定壓比熱；分別為為介質(zhì)等效體積熱容和等效導(dǎo)熱系數(shù)。

空氣區(qū)域為流體區(qū)，采用如下控制方程（N-S）方程：

連續(xù)性方程為

計算中應(yīng)用顯熱容法，假設(shè)模型中含水介質(zhì)相變發(fā)生在溫度區(qū)間（Tm±ΔT*）。當(dāng)建立等效體積熱容時，應(yīng)考慮溫度間隔 ΔT*的影響，同時，假設(shè)介質(zhì)在正凍、未凍時的體積熱容Cf和Cu以及導(dǎo)熱系數(shù)λf和λu不取決于溫度，因此，簡化構(gòu)造出和的表達(dá)式：

式中：L為含水介質(zhì)單位體積相變潛能。

3 計算模型

3.1 工點概況

本研究以青藏鐵路北麓河試驗段拋石護(hù)坡路基結(jié)構(gòu)為原型。該試驗段位于青藏高原可可西里與風(fēng)火山之間，北麓河盆地南部，地貌形態(tài)屬北麓河沖洪積高原地貌。區(qū)域地勢開闊，地形略有起伏，沖溝發(fā)育，植被覆蓋也較好。該地區(qū)屬青藏高原干旱氣候區(qū)，寒冷干旱，空氣稀薄。一年內(nèi)凍結(jié)期長達(dá) 7～8個月，蒸發(fā)量大，年平均氣溫為-5.2 ℃，最高溫為23.2 ℃，最低溫為37.7 ℃，年平均降雨量為290.9 mm，相對濕度平均為57%；最大風(fēng)速為40 m/s，年平均風(fēng)速為4.10m/s，主導(dǎo)風(fēng)向為北西。

青藏鐵路北麓河試驗段為厚層地下冰地段，地下冰分布極為廣泛，鉆探資料表明，該地段全范圍內(nèi)均有含土冰層分布，體積含冰量大于50%，多為懸浮狀構(gòu)造。地層上部為粉、細(xì)砂夾碎塊石，下部活動層范圍主要為棕紅色粉質(zhì)黏土，上限附近有1.0～3.0 m厚的含土冰層，下部以棕紅色全風(fēng)化泥巖為主，局部有砂巖夾層。

3.2 計算模型

選取青藏鐵路北麓河試驗段拋石護(hù)坡路基結(jié)構(gòu)為計算斷面。拋石護(hù)坡厚度為80 cm，拋石粒徑為8～15 cm，路基填筑高度為3.7 m，地層材料包括拋石護(hù)坡、路堤填土、活動層以及多年凍土層。路基上方區(qū)域為流體區(qū)域，考慮路基左右護(hù)坡通風(fēng)條件不同，路堤計算邊界條件并不對稱，故本研究中仍以整體路基結(jié)構(gòu)為計算對象。

圖1為拋石護(hù)坡路基的計算模型。計算域以路基坡腳向外延伸各 30 m，天然地表上下各延伸30 m。圖中，區(qū)域①為路基砂礫料填土；區(qū)域②為細(xì)砂活動層；區(qū)域③、④分別為黏土活動層；區(qū)域⑤為黏土；區(qū)域⑥為強(qiáng)風(fēng)化泥巖；區(qū)域⑦為拋石護(hù)坡；區(qū)域⑧為流體區(qū)；A～N為區(qū)域邊界。

圖1 計算模型（單位：m）Fig.1 Calculation model (unit: m)

3.3 材料參數(shù)

研究工點各介質(zhì)的物理參數(shù)見表 1。表中，L為介質(zhì)的相變潛能；拋石層滲透率k=1.58×10-6；慣性阻力系數(shù)（非達(dá)西流的Beta因子）B=1 650.24 m-1，其他符號意義同前。

表1 路基結(jié)構(gòu)中各介質(zhì)的物理參數(shù)Table 1 Thermal parameters of media in embankments

3.4 邊界條件

考慮全球氣候變暖的影響，根據(jù)文獻(xiàn)[8]，取青藏鐵路未來50年年平均氣溫上升2.6 ℃，設(shè)初始年平均氣溫為-4.0 ℃，由于路基模型上表面溫度變化不僅與環(huán)境空氣有關(guān)，且受太陽輻射等復(fù)雜因素的綜合作用，故根據(jù)吸附面理論及北麓河氣象站的長期現(xiàn)場觀測資料，對計算域的熱邊界條件進(jìn)行如下假定。

氣溫變化規(guī)律：

式中：t為時間變量（h），當(dāng)t=0對應(yīng)的初始時間為 7月20日；假定通過JK邊的熱流密度為q =0.06 W/m2，邊界AMJ和DNK視為絕熱。

由于拋石護(hù)坡中空氣與外界連通，其內(nèi)部溫度場和速度場直接受到外界通風(fēng)條件的影響，所以拋石層與空氣的接觸面EF和GH設(shè)置為流體可通過界面，在這些界面上施加溫度荷載的同時，外界風(fēng)可通過這些界面穿過拋石層。此外，為了使計算與實際更為接近，計算中將風(fēng)速邊界施加于距離坡腳30 m處的AB邊界上，即取模型邊界AB為風(fēng)速的入口邊界，氣溫按式（11）變化。

對北麓河氣象站觀測資料分析，計算工點地面以上10.85 m高度標(biāo)準(zhǔn)風(fēng)速變化規(guī)律可近似為

根據(jù)“綜合冪次律”理論[9]，考慮到試驗段的主導(dǎo)風(fēng)向為NW，故入口邊界AB的風(fēng)速可近似表達(dá)為

式中：V為任意高度風(fēng)速；V10.85為10.85 m高度的標(biāo)準(zhǔn)風(fēng)速；y、y10.85分別為V、V10.85對應(yīng)的高度；α值由試驗獲得，本次計算取0.16。

根據(jù)北麓河氣象站的實際觀測資料，熱初始條件按式（12），不考慮升溫的天然地表溫度方程作為區(qū)域②～⑥的上邊界條件進(jìn)行反復(fù)多年計算，逐時段求解，直至得到穩(wěn)定的溫度場為止；此時年變化層以下穩(wěn)定場基本保持穩(wěn)定，年變化層以上相同位置的溫度值在同一時間逐年相同，本模型取5月31日溫度場作為該區(qū)計算的溫度初始條件。速度初始條件取不考慮能量方程式（4）和式（8），定溫條件下以式（16）為邊界條件，計算得到的5月31日穩(wěn)定速度場作為計算流體區(qū)速度初始條件。

4 計算結(jié)果與分析

4.1 模型的合理性驗證

以青藏鐵路北麓河試驗段拋石護(hù)坡路基結(jié)構(gòu)斷面為研究對象，對本文所建立的計算模型的合理性進(jìn)行驗證。根據(jù)現(xiàn)場觀測資料，對計算域的熱邊界條件作如下假定：

該試驗段于2004年9月施工完畢，計算中以路基在2004年10月5日實測地溫數(shù)據(jù)作為計算的熱初始邊界，以式（18）為速度邊界計算得到的相應(yīng)時刻溫度速度場作為初始速度條件。

圖2、3為青藏鐵路北麓河試驗段拋石護(hù)坡路基結(jié)構(gòu)斷面在夏季8月30日和冬季12月27日溫度場模擬與實測結(jié)果對比圖。本文中溫度場實測資料由中國科學(xué)院寒區(qū)旱區(qū)環(huán)境與工程研究所提供，從圖中可以看出，模擬所得溫度場與實測值分布規(guī)律基本一致，路堤的修筑使下部凍土上限有所抬升，存在地溫升高的現(xiàn)象。圖4為青藏鐵路北麓河試驗段拋石護(hù)坡路基結(jié)構(gòu)斷面在環(huán)境風(fēng)速最大時刻的速度分布對比示意圖。從圖中可以看出，左側(cè)拋石護(hù)坡中空氣在外界大氣風(fēng)場的強(qiáng)迫作用下，基本沿拋石護(hù)坡斜向上運(yùn)動，靠近護(hù)坡外表面的空氣流速大，護(hù)坡下部速度較小，模擬空氣流速在 1.4×10-3～13.1 m/s之間，實測空氣流速在1.6×10-3～14.6 m/s之間，模擬結(jié)果與實測結(jié)果基本一致。

圖2 拋石護(hù)坡路基夏季溫度場（單位：℃）Fig.2 Temperature distributions of the crushed-rock slope embankment in summer (unit: ℃)

圖3 拋石護(hù)坡路基冬季溫度場（單位：℃）Fig.3 Temperature distributions of the crushed-rock slope embankment in winter (unit: ℃)

圖4 拋石護(hù)坡路基斷面左側(cè)護(hù)坡內(nèi)速度場Fig.4 Velocity distribution of the pore air in crushed-rock side slope

從上述結(jié)果中也可以看出，模擬溫度場和實測溫度場的分布形態(tài)略有差異，整體而言，本文建立的計算模型及參數(shù)的獲取方法是合理的，能夠反映出拋石護(hù)坡路基在強(qiáng)通風(fēng)條件下的主要特性。

4.2 溫度場分析

圖5為拋石護(hù)坡路基修筑完成2年后的溫度分布。由圖中可以看出，拋石層主要以導(dǎo)熱為主，越遠(yuǎn)離環(huán)境，則凍土層內(nèi)的溫度變化受環(huán)境的影響越小。此外，由于強(qiáng)風(fēng)的方向為從左向右，所以路堤左側(cè)護(hù)坡附近空氣與填土的換熱能力高于右側(cè)，故左側(cè)填土接近大氣來流方向的溫度變化較右側(cè)快，左側(cè)等溫線的斜率略高于右側(cè)。由圖 5(a)可見，路基在7月15日最高溫度升至10 ℃以上，最低溫度為-2.5 ℃，路基中部以下產(chǎn)生了一個-2.5 ℃的低溫凍土核，其溫度低于周圍土體約-1.5 ℃。天然地表下，凍土的最大融深坐標(biāo)y=-1.5 m，而路基中線處y=0.7 m，表明拋石護(hù)坡路基的修筑，使多年凍土上限提升了約 2.2 m。由圖 5(b)可見，路基在次年冬季1月15日已全部完成回凍，但天然地表以下局部仍有0 ℃等溫線存在，表明拋石護(hù)坡路基下部土體回凍速度較天然地表下部土體更快。

圖6為拋石護(hù)坡路基修筑完成50年后的溫度分布。由圖中曲線看出，隨著全球氣候變暖，凍土的整體溫度有所升高，多年凍土上限（0 ℃等溫線）有所下降。由圖6(a)可見，夏季7月15日最高氣溫升至13 ℃，最低氣溫為-2 ℃，路基中線以下土體的內(nèi)部產(chǎn)生了“似眼球狀”融化夾層，由圖6(b)可見，盡管冬季外界溫度較低，而“似眼球狀”融化夾層仍然存在，且范圍有所增大，表明該融化夾層貫穿于拋石護(hù)坡路基修筑完成50年的整個凍融區(qū)，不利于路基的穩(wěn)定。

圖5 修筑完成后第二年溫度分布（單位：℃）Fig.5 Temperature distributions of the crushed-rock slope embankment after 2 years of the construction (unit：℃)

圖6 修筑完成后第50年溫度分布（單位：℃）Fig.6 Temperature distributions of the crushed-rock slope embankment after 50 years of the construction (unit: ℃)

4.3 速度場分析

強(qiáng)通風(fēng)條件下模型中的空氣運(yùn)動主要以強(qiáng)迫對流為主，故不同時刻空氣運(yùn)動規(guī)律差別不大，本文僅以環(huán)境風(fēng)速最大時刻的速度場為例對拋石護(hù)坡路基的速度場進(jìn)行分析。

圖7為本次計算環(huán)境風(fēng)速最大時刻拋石護(hù)坡路基的速度分布示意圖，由于拋石層內(nèi)外空氣運(yùn)動差異較大，所有圖中矢量只表示流體的運(yùn)動方向，不代表大小。由圖可見，計算域的風(fēng)向自左向右，故流通區(qū)矢量的方向以向右為主，迎風(fēng)拋石護(hù)坡層中空氣運(yùn)動方向大致為沿護(hù)坡斜向上，背風(fēng)拋石護(hù)坡層中空氣運(yùn)動方向以從下到上運(yùn)動為主。圖8為左側(cè)拋石護(hù)坡層細(xì)部的速度分布圖。由圖可以看出，拋石層中空氣運(yùn)動較為零亂，空氣在拋石縫隙處繞流拋石，整體而言，拋石孔隙中空氣運(yùn)動形式為“繞流”。拋石層表面空氣速度最大，內(nèi)部較小，空氣速度分布區(qū)間為 1.24×10-3～12.8 m/s，與現(xiàn)場試驗測得的風(fēng)速區(qū)間基本一致。

圖7 拋石護(hù)坡路基速度場分布Fig.7 Velocity vectors at maximum wind

圖8 左側(cè)拋石護(hù)坡路基速度場分布（細(xì)部）（單位：m/s）Fig.8 Velocity distribution of the pore air in crushed-rock side slope at maximum wind（unit: m/s）

5 結(jié) 論

（1）強(qiáng)通風(fēng)條件下，拋石護(hù)坡路基對保護(hù)多年凍土具有積極作用。路基修筑完成2年的溫度計算結(jié)果表明，拋石護(hù)坡路基的存在使夏季多年凍土上限提高了約 2.2 m，冬季拋石護(hù)坡路基下部土體回凍速度較天然地表下部土體更快。

（2）降溫作用主要集中在護(hù)坡附近有限范圍之內(nèi)，對路基中部的降溫作用相對較弱，拋石護(hù)坡對凍土路基本體的保護(hù)作用有限。從長期降溫效果來看，由于全球氣候變暖的影響，強(qiáng)通風(fēng)條件下拋石護(hù)坡路基中線以下土體的內(nèi)部可能產(chǎn)生“似眼球狀”融化夾層，不利于路基的穩(wěn)定。

（3）迎風(fēng)拋石護(hù)坡層中空氣運(yùn)動方向大致為沿護(hù)坡斜向上，背風(fēng)拋石護(hù)坡層中空氣運(yùn)動方向以從下到上運(yùn)動為主，拋石層內(nèi)空氣的運(yùn)動形式為“繞流”，拋石層表面空氣速度最大，內(nèi)部較小，空氣速度分布區(qū)間為 1.24×10-3～12.8 m/s，與現(xiàn)場試驗測得的風(fēng)速區(qū)間基本一致。

[1] SUN Zhi-zhong, MA Wei, LI Dong-qing.In situ test on cooling effectiveness of air convection embankment with crushed rock slope protection in permafrost regions[J].Journal of Cold Regions Engineering, 2005, 19(2): 38－51.

[2] 吳青柏, 趙世運(yùn), 馬巍, 等.青藏鐵路塊石路基結(jié)構(gòu)的冷卻效果監(jiān)測分析[J].巖土工程學(xué)報, 2005, 27(12):1386－1390.WU Qing-bai, ZHAO Shi-yun, MA Wei, et al.Monitoring and analysis of cooling effect of block-stone embankment for Qinghai-Tibet railway[J].Chinese Journal of Geotechnical Engineering, 2005, 27(12): 1386－1390.

[3] NETO H L, QUARESMA J N, COTTA R M.Natural convection in three-dimensional porous cavities: Integral transform method[J].International Journal of Heat andMass Transfer, 2002, 45(14): 3013－3022.

[4] AHMED N, SUNDADA D K.Nonlinear flow in porous media[J].Journal of the Hydraulic Division, ASCE,1969, 95: 1847－1857.

[5] 米隆, 賴遠(yuǎn)明, 張克華.凍土通風(fēng)路基溫度場的三維非線性分析[J].冰川凍土, 2002, 24(6): 765－769.MI Long, LAI Yuan-ming, ZHANG Ke-hua.3D nonlininear analyses for temperature field of ventilative embankment in cold regions[J].Journal of Glaciolgy and Geocryology, 2002, 24(6): 765－769.

[6] 賴遠(yuǎn)明, 張魯新, 徐偉澤.青藏鐵路拋石路基的溫度特性研究[J].冰川凍土, 2003, 25(6): 291－296.LAI Yuan-ming, ZHANG Lu-xin, XU Wei-ze.Temperature features of broken rock mass embankment in the Qinghai-Tibetan railway[J].Journal of Glaciolgy and Geocryology, 2003, 25(6): 291－296.

[7] 姜凡, 劉石, 王海剛, 等.凍土碎（塊）石路堤傳熱與對流特性數(shù)值計算分析[J].中國科學(xué)（D輯）, 2003,33(增刊): 133－144.JIANG Fan, LIU Shi, WANG Hai-gang, et al.Numerical analysis on the heat transfer and convection characteristics of broken rock mass embankment[J].Science in China (Series D), 2003, 33(Supp.): 133－144.

[8] 秦大河.中國西部環(huán)境演變評估綜合報告[M].北京:科學(xué)出版社, 2002.

[9] 貝爾.多孔介質(zhì)流體動力學(xué)[M].北京: 中國建筑工業(yè)出版社, 1983.