仝 睿，付 浩，宋二祥

(清華大學(xué)土木工程系土木工程安全與耐久教育部重點試驗室，北京 100084)

淺層土體溫度的分布規(guī)律，直接影響土體中水分的蒸發(fā)、凍結(jié)等行為，從而影響到農(nóng)業(yè)生產(chǎn)、巖土工程等領(lǐng)域的研究建設(shè)[1-3]。在地表存在公路、鐵路等的覆蓋層時，冬季土體底部熱、頂部冷的溫度梯度，還會導(dǎo)致覆蓋效應(yīng)出現(xiàn)[4-7]。如何評價頂部存在覆蓋層時，土體溫度場的分布，對于探究覆蓋效應(yīng)的發(fā)生機理、實現(xiàn)凍土區(qū)路基的合理設(shè)計具有重要意義。

淺層土體溫度常常受太陽輻射、水分蒸發(fā)、大氣溫度等多方面因素影響。但如果將地表溫度變化視為正弦變化的溫度邊界條件，深處恒溫土層視為固定的溫度邊界條件，忽略傳熱外的其他作用因素，則土層溫度分布可以近似為一維瞬態(tài)傳熱問題。徐學(xué)祖等[8]在《凍土物理學(xué)》中的研究表明：“土體內(nèi)地溫隨時間的變化主要是熱量輸運的影響，數(shù)學(xué)上可歸納為熱傳導(dǎo)問題”。傳熱方程反映的熱量輸運情況，對地溫如何隨時間變化有著非常大的影響。因此，研究該一維傳熱方程的定解，對于估計實際土體內(nèi)溫度場的分布，具有重要的參考價值。

實際土體多為分層材料。每層土體內(nèi)，物理性質(zhì)會有波動，但近似為均質(zhì)材料可以簡化計算，且并不會導(dǎo)致較大的誤差。在《凍土物理學(xué)》中，即有將單一土層視為均質(zhì)材料，進行熱流計算的算例[8]。考慮到傳熱方程涉及的熱容量、傳熱系數(shù)、密度等參數(shù)主要受到含水率影響，當土體內(nèi)含水量隨時間變化不大時，將土體簡化為均質(zhì)材料進行分析是可行的。此外，本文隨后的推導(dǎo)表明，單一土層的分析，可以進一步推廣到分層土模型中。因此，對單一均質(zhì)土層的一維傳熱方程定解進行分析，對于工程問題中的多層土層溫度變化分析同樣具有參考價值。

關(guān)于一維傳熱方程的解析求解，此前已有一定研究。涂新斌等[9]，使用復(fù)變函數(shù)方法，分析了模型無限大，底部不存在恒溫點時，一維傳熱方程對應(yīng)的解析解。李翊神[10]使用杜哈梅積分法，研究了傳熱方程定解和波動方程定解之間的關(guān)系。王兆瑞等[11]使用傅里葉變換和積分函數(shù)表示方法，分析了一維傳熱方程的特解。林府標等[12]通過李代數(shù)，研究了一維廣義方程的解析解特點。左沖等[13]使用時域徑向積分法給出了積分形式的傳熱方程解。但是，目前尚未有人給出底部溫度恒定、頂部溫度變化邊界條件下，土體一維瞬態(tài)傳熱問題的解析解。

本文給出一維瞬態(tài)傳熱問題的解析解。通過與數(shù)值模擬及現(xiàn)場測試結(jié)果的對比驗證了解析解的正確性，進而給出了解析解在凍深估計、有關(guān)模型試驗中幾何相似關(guān)系的確定、“覆蓋效應(yīng)”數(shù)值模擬等方面的應(yīng)用。

1 一維瞬態(tài)傳熱問題解析解

1.1 解析解推導(dǎo)

土體或上覆材料的頂部，與大氣直接接觸，可近似用正弦曲線反映頂部空氣溫度每日或每年的周期性變化。而在底部較深處，則假定溫度基本恒定將，土近似為均質(zhì)材料，得到的控制方程及邊界條件為：

式中：u為單元體溫度；ρ 為材料密度；c為材料比熱容；k為材料導(dǎo)熱系數(shù)；T0為頂部溫度均值；T1為頂部溫度正弦變化部分振幅；Tbot為深層土體恒定溫度；Tini為初始狀態(tài)的溫度分布。

式(1)可根據(jù)數(shù)學(xué)物理方程中的參數(shù)分離法進行求解[14]：

式中：B1、B2為待定系數(shù)；λ 為變量分離法選取的常數(shù)。

對于頂部與底部的固定值部分，存在對應(yīng)特解：

對于頂部的正弦邊界條件，設(shè)對應(yīng)的溫度解形式為：

根據(jù)特解需要滿足的傳熱方程式(1)，解出對應(yīng)的系數(shù)之間的關(guān)系：

又因為x=0 時，需符合頂部溫度邊界條件，可以得到對應(yīng)Tsin(ωx)的特解為：

從而，針對非齊次邊界條件(2)，對應(yīng)的特解為：

假如初始值條件選取為：

則這組特解反映的是頂部正弦溫度邊界條件，會形成振幅不斷衰減的正弦溫度波動，向土體下方傳播。

當波動傳播到土體底部時，特解式(8)中，取x=h，得到的邊界條件無法和底部u(h,t)=Tbot的邊界條件相協(xié)調(diào)，說明特解式(8)并不能同時符合上、下端邊界條件。此時需要如圖1 所示，疊加函數(shù)式(10)，才能滿足底部的邊界條件。

圖1 疊加示意圖Fig. 1 Composition of functions

而函數(shù)式(10)又會導(dǎo)致頂部邊界條件無法滿足。將反射函數(shù)不斷疊加下去，最終得到了同時滿足上、下邊界條件，級數(shù)形式的精確解：

這兩個無窮級數(shù)，都可以通過比較判別法證明絕對收斂。以第一個級數(shù)和為例：

從而可以知道，無窮級數(shù)V(x,t)收斂。

注意到精確解(式(11))中，不斷疊加的級數(shù)項，反映了溫度波動變化傳播到土體底部時，存在反射現(xiàn)象。反射波動回到土體頂部時，又會再次反射。而且式(10)中，與應(yīng)力波、光波等波動類似，在固定端反射時，都出現(xiàn)了“半波損失”的現(xiàn)象。

如果給定的初始值條件與式(9)存在差異，不妨設(shè)：

將w(x)以齊次通解式(3)為基底進行傅里葉展開，可以得到對應(yīng)的齊次解：

其中：

土體或覆蓋層頂部的空氣溫度變化，對應(yīng)的頂部溫度邊界條件應(yīng)當包含每日溫度變化、每年溫度變化兩組正弦變化函數(shù)。分別對應(yīng)：

1.2 近似解及其討論

式(10)給出的精確解，表達形式較為復(fù)雜。可以將特殊解式(8)，作為精確解式(11)的近似。近似解的相對誤差大小為：

取一組典型數(shù)值進行估計。土的密度為1600 kg/m3，比熱容1200 J/(kg·K)，傳熱系數(shù)k=1.2 W/(m·K)，h=20 m，溫度年波動T1=20 K，土體底部溫度Tbot=278.15 K，此時按式(17)計算相對誤差上限約為0.0024%。這意味著，針對自然界中的土體傳熱問題，式(8)可以很好地近似實際溫度變化。

2 解析解的數(shù)值驗證

由圖2 可見，位于0.1 倍、0.2 倍、0.5 倍波長位置處的土體，在經(jīng)過1～2 個周期后數(shù)值模擬的溫度已與解析解接近一致，由此可以看出解析解作為穩(wěn)態(tài)解的正確性。此外，初始時刻溫度的差異正說明初始溫度場的影響需要一段時間來平衡。

圖2 解析解和模擬結(jié)果對比Fig. 2 Comparison between analytical solution and simulation results

3 結(jié)合現(xiàn)場試驗數(shù)據(jù)對解析解進行驗證

本課題組此前在蘭新鐵路武威段411+600 區(qū)域路堤土層中不同深度處布設(shè)溫、濕度傳感器和凍脹計以獲得相應(yīng)的測試數(shù)據(jù)[15]。

本課題組將監(jiān)測儀器布置在0.1 m、0.4 m、1.0 m、1.6 m、2.5 m、3.3 m 處，以獲得對應(yīng)深度處的溫濕度數(shù)據(jù)?？紤]到2019/06/01～2021/06/01這一時間段內(nèi)測得的溫度數(shù)據(jù)和2017/11/11～019/06/01 這一時間段內(nèi)測得的溫度數(shù)據(jù)相差不大，本文選用后一時間段內(nèi)的溫度數(shù)據(jù)，對解析解加以驗證。該時間段內(nèi)的路基土體溫度數(shù)據(jù)如圖3 所示，邊坡土體溫度數(shù)據(jù)如圖4 所示。

圖3 路基土體溫度Fig. 3 Temperature of subgrade soil

圖4 邊坡土體溫度Fig. 4 Temperature of slope soil

觀察圖3 和圖4 可以發(fā)現(xiàn)，土體不同深度處的溫度都近似以正弦函數(shù)的形式周期性變化。正弦函數(shù)的振幅隨著深度的增加不斷減小，相位隨深度增加也存在明顯的“滯后性”現(xiàn)象。

3.1 振幅衰減驗證

根據(jù)不同深度處土體的最高溫、最低溫計算溫度變化的振幅，進而根據(jù)振幅的與深度的關(guān)系，對解析解進行驗證。

路基中土體溫度的振幅如圖5 所示。使用指數(shù)擬合公式進行回歸，得到如下公式：

圖5 路基振幅-深度關(guān)系圖Fig. 5 Amplitude vs depth diagram of subgrade

該回歸公式的R2為0.9881，回歸的效果較好。根據(jù)此前測得的土體密度1600 kg/m3，比熱容按照1200 J/(kg·K)進行估計，回歸參數(shù)對應(yīng)的土體傳熱系數(shù)為：

該取值在粉土正常傳熱系數(shù)范圍0.6 W/(m·K)～1.4 W/(m·K)之內(nèi)。

邊坡中土體溫度的振幅如圖6 所示。

圖6 邊坡振幅-深度關(guān)系圖Fig. 6 Amplitude vs depth diagram of slope

使用指數(shù)擬合公式進行回歸，得到下式：

該回歸公式的R2為0.989。根據(jù)此前測得的土體密度1600 kg/m3、比熱容按照1200 J/(kg·K)進行估計，回歸參數(shù)對應(yīng)的土體傳熱系數(shù)為：

與根據(jù)路基中溫度變化情況回歸出的土體傳熱系數(shù)極為接近。

根據(jù)式(19)和式(21)可以得到，路基頂部的石砟層，將大氣溫度變化衰減為：

根據(jù)鐵道部數(shù)據(jù)，石砟表觀密度在2700 kg/m3左右，堆積密度1500 kg/m3左右。據(jù)此推算，路基頂部石砟層孔隙率在44.4%左右。石砟層厚度按照平均值0.65 m 計算。石砟比熱容按照1000 J/(kg·K)計算?？諝鈱?dǎo)熱系數(shù)為0.023 W/(m·K)，石砟導(dǎo)熱系數(shù)取為1.28 W/(m·K)。估算得到的等價導(dǎo)熱系數(shù)為：

根據(jù)這個估算的年溫度變化衰減為：

這一估計數(shù)值與式(23)的實際計算數(shù)值存在較大的差距。存在差距的原因?qū)⒃诘谒墓?jié)最后進行解釋。

3.2 相位滯后性驗證

將不同傳感器測得最高溫的時間點，與傳感器的深度進行回歸，從而對解析解進行驗證。

路基中土體取得最高溫的時間點如圖7 所示。其中縱軸的數(shù)值(時間/d)，代表傳感器自2017 年11 月11 日開始測量數(shù)據(jù)以來，經(jīng)過了多少天測得了最高溫數(shù)據(jù)。使用線性擬合公式進行回歸，回歸得到如下公式：

圖7 路基最高溫時間點-深度關(guān)系圖Fig. 7 Maximum temperature time vs depth diagram of subgrade

該回歸公式的R2為0.9883，回歸的效果較好。根據(jù)此前測得的土體密度1600 kg/m3，比熱容按照1200 J/(kg·K)進行估計，回歸參數(shù)對應(yīng)的土體傳熱系數(shù)為：

這個取值在粉土正常傳熱系數(shù)范圍0.6 W/(m·K)～1.4 W/(m·K)之內(nèi)。

邊坡中土體取得最高溫的時間點如圖8 所示。使用線性擬合公式進行回歸，回歸得到公式：該回歸公式的R2為0.9793。根據(jù)此前測得的土體密度1600 kg/m3，比熱容按照1200 J/(kg·K)進行估計，回歸參數(shù)對應(yīng)的土體傳熱系數(shù)為：

圖8 邊坡最高溫時間點-深度關(guān)系圖Fig. 8 Maximum temperature time vs depth diagram of slope

這個取值在粉土正常傳熱系數(shù)范圍0.6 W/(m·K)～1.4 W/(m·K)之內(nèi)。

值得注意的是，根據(jù)相位滯后性反推的土體傳熱系數(shù)，不論路基和邊坡，都明顯高于根據(jù)振幅衰減情況反推得到的土體傳熱系數(shù)。

觀察到溫度變化的傳播函數(shù)與波的傳播函數(shù)類似，且式(11)推導(dǎo)了固定溫度邊界產(chǎn)生的反射波存在，及其對應(yīng)的“半波損失”現(xiàn)象。筆者有一個猜想：“溫度變化傳播到不同性質(zhì)材料的交界面上，也會發(fā)生反射和入射現(xiàn)象，帶來振幅的改變”，筆者將在后文對這一猜想加以證明。

4 不同性質(zhì)材料交界面的反射函數(shù)推導(dǎo)

考慮兩種不同材料交界面處情況，其中，上方材料厚度為x1、密度為ρ1、熱容量為c1、傳熱系數(shù)為k1，下方材料厚度為x2、密度為ρ2、熱容量為c2、傳熱系數(shù)為k2，溫度邊界條件設(shè)為：

初始值條件設(shè)為：

參考式(7)，從t=0 時刻開始，頂部正弦溫度變化逐漸向下傳輸。而當溫度傳輸?shù)浇唤缑嫔蠒r，需要在交界面上同時滿足溫度的連續(xù)性條件，以及熱流量的平衡條件。

假設(shè)上方入射溫度方程為：

假設(shè)入射后溫度方程為：

觀察到熱流量：

如果只存在入射現(xiàn)象，則可以得到方程：

由式(35)中的第一個溫度連續(xù)性方程可知，未知參數(shù)A1=1。代入第二個熱流量平衡方程，可以看到，除非兩種材料的c、ρ、k值相等，不然熱流量平衡方程肯定無法滿足。

與振動波、光波類似，這里可以取cρk這一乘積值作為表征材料對熱量傳輸能力的參數(shù)。

當兩種材料的值cρk不同時，為同時滿足交界面處的溫度連續(xù)性方程、熱流量平衡方程，就需要引入反射波。假設(shè)反射波方程為：

從而在交界面處，得到方程：

化簡為：

從而可以解出：

從解式(39)可以看出，反射溫度函數(shù)的振幅一定小于入射溫度函數(shù)。當溫度函數(shù)從cρk值比較大的材料傳入cρk值比較小的材料時，入射后溫度函數(shù)的振幅會大于入射前溫度函數(shù)的振幅，同時反射溫度函數(shù)的振幅系數(shù)為正。而當溫度函數(shù)從cρk比較小的材料傳入cρk值比較大的材料時，入射后溫度函數(shù)的振幅會小于入射前溫度函數(shù)的振幅，同時反射溫度函數(shù)的振幅系數(shù)為負，存在“半波損失”現(xiàn)象，如圖9 所示。

圖9 不同材料交界面反射示意圖Fig. 9 Reflection in interface between different materials

從本節(jié)的推導(dǎo)中可以看出，溫度函數(shù)的振幅變化會受到不同材料交界面的影響。而相位變化，則不會受到不同材料交界面的影響。因此，式(28)、式(29)根據(jù)相位變化回歸得到的熱傳導(dǎo)系數(shù)，更接近于真實的土體平均熱傳導(dǎo)系數(shù)。

而式(26)估算得到的比率，明顯小于實際衰減后的振幅比率，也是因為沒有考慮交界面處反射波。如果將交界面處的反射波納入考慮之中，將土體傳熱系數(shù)按照式(28)估計為0.989 W/(m·K)。重新估算得到的衰減后比率為：

這一比率略大于實際值59.48%。考慮到實際石子與石子之間不是完全貼合，石砟層的傳熱系數(shù)應(yīng)略小于按照空氣、石子并聯(lián)估算的式(24)，式(40)得到的比率略高于實際值，屬于較為合理的結(jié)果。

5 反射函數(shù)的數(shù)值模擬驗證

使用COMSOL 對解得的反射函數(shù)進行驗證。上方材料厚度為x1=1 m、密度為ρ1=1200 kg/m3、熱容量為c1=1400 J/(kg·K)、傳熱系數(shù)為k1=1.0 W/(m·K)，下方材料厚度為x2=1 m、密度為ρ2=1400 kg/m3、熱容量為c2=1600 J/(kg·K)、傳熱系數(shù)為k2=1.0 W/(m·K)。溫度邊界條件設(shè)為：

初始值條件設(shè)為：

模擬到10.25 d 時，交界面附近的數(shù)值解、解析解如圖10 所示。補充了忽略熱流量平衡方程得到的無反射波參考解作為對比。

圖10 中，仿真結(jié)果與理論分析結(jié)果基本一致，而無反射對應(yīng)解則存在明顯偏差。

圖10 仿真結(jié)果與理論分析結(jié)果對照圖Fig. 10 Comparison of simulation result and analytical result

6 解析解的應(yīng)用

6.1 使用解析解對冬季凍深進行估計

假設(shè)某地區(qū)地大氣溫度變化遵循公式：

考慮完全裸露均質(zhì)土體，忽略衰減較多的日溫侵入。在地下x/m 處對應(yīng)的溫度為：

求解T＞273.15 K，得到：

以西安地區(qū)為例，文獻[15]測得西安地區(qū)粉質(zhì)粘土比熱容為1231.6 J/kg·K，傳熱系數(shù)為1.9885 W/(m·K)，密度為1959.184 kg/m3。西安地區(qū)平均溫度T0為288.65 K，年溫度變化振幅約為18.5 K。根據(jù)式(45)估算，凍深為：

查閱相關(guān)文獻[16]知，西安地區(qū)冬季觀測到的凍土深度為0.45 m。估計的凍深與實際的凍深相比，誤差為20%?？紤]到實際情況土層并非是均勻介質(zhì)，且土體內(nèi)水分遷移、化學(xué)反應(yīng)等因素對于土體的溫度分布會產(chǎn)生影響，這一誤差較為合理。20%的誤差，可以滿足估測的需要。

6.2 使用解析解估計室內(nèi)試驗的“尺寸效應(yīng)”

在使用小尺寸樣本對實際大尺寸情況進行試驗?zāi)M時，一個常見的問題是：頂部溫度條件如何選取，才可以模擬實際自然界中的溫度變化情況。

例如，室內(nèi)試驗裝置高度1 m，用于模擬實際10 m 深的土體溫度變化。則對應(yīng)的，模擬年溫度變化的室內(nèi)試驗頂部溫度變化周期，應(yīng)當選取為365/100=3.65 d。

6.3 使用解析解對“覆蓋效應(yīng)”模擬條件進行修正

此前針對“覆蓋效應(yīng)”現(xiàn)象的數(shù)值模擬，頂部邊界往往直接設(shè)置為大氣溫度變化[4,17-18]；或者設(shè)置為固定溫度[19]。實際情況下，頂部道面也會對大氣溫度變化有折減和滯后作用?？紤]這一情況，根據(jù)頂部道面情況，對覆蓋效應(yīng)的頂部邊界條件進行修正。

例如，瀝青道面常見厚度為0.2 m，瀝青密度為1200 kg/m3、比熱容為1300 J/(kg·K)、導(dǎo)熱系數(shù)為0.7 W/(m·K)，下方土層密度為1600 kg/m3、比熱容為1200 J/(kg·K)，導(dǎo)熱系數(shù)為1 W/(m·K)。

此前本課題組的論文[4]中，選擇頂部模擬條件為：年平均溫度為283.15 K、日溫度變化振幅為10 K、年溫度變化振幅為15 K。

考慮交界面處的反射現(xiàn)象計算得到的日溫度變化，經(jīng)過瀝青道面衰減為：

年溫度變化衰減為：

對應(yīng)的，修正后大氣溫度變化條件近似為：

除了頂部溫度邊界條件外，從第1 節(jié)也可以看出，模擬需要設(shè)置合理的初始值條件。實際土體經(jīng)過多年頂部溫度作用，土體內(nèi)溫度分布已經(jīng)近似于式(9)，因此，模擬的初始值應(yīng)按照式(9)進行設(shè)置，才可以得到較為接近實際情況的水熱遷移結(jié)果。

從式(10)可知，模型高度需要滿足：

才可以忽略底部固定邊界條件產(chǎn)生的反射波影響。因此，將模型高度修改為10 m。

修正前，模型的初始值條件為全域288.15 K。修正后的初始值條件設(shè)計為：

其余模擬參數(shù)按照文獻[4]進行設(shè)置。

修正前、后冬季溫度對比如圖11 所示。

圖11 修正前、后溫度對比圖Fig. 11 Temperature comparison before and after modification

從圖11 中可以看出，修正后土體溫度波動更小，冬季頂部溫度更高。

修正前、后含水量對比如圖12 所示。

圖12 修正前、后含水量對比圖Fig. 12 Water volumetric content comparison before and after modification

從圖12 中可以看出，修正后，土體頂部的水分聚集更少。

修正前、后的結(jié)果對比表明：修正前的模型，夸大了路基土體頂部的溫度變化幅度，得到的水分遷移量偏大，對“覆蓋效應(yīng)”影響的預(yù)測偏大。

7 結(jié)論

本文針對一維土柱的瞬態(tài)問題，求解了對應(yīng)的解析解。通過數(shù)值模擬、現(xiàn)場試驗兩種方式驗證了解析解的正確性。本文還推導(dǎo)了解析解在固定邊界處、不同材料交界面上存在的反射現(xiàn)象。通過對推導(dǎo)結(jié)果的分析，發(fā)現(xiàn)：

(1) 土體頂部溫度正弦變化帶來的土體溫度改變與應(yīng)力波存在一定相似性，會以波動的形式向下傳播；但波動振幅在傳播過程中會以指數(shù)形式不斷衰減。

(2) 溫度正弦變化的傳播在固定邊界處存在帶有“半波損失”的反射現(xiàn)象；在不同材料交界面處也存在反射現(xiàn)象。

(3) 該解析解存在多種應(yīng)用場景?？梢愿鶕?jù)該解析解，估計土體某一深度處的溫度變化情況；估計冬季凍結(jié)深度；計算室內(nèi)試驗滿足“尺寸效應(yīng)”的合理邊界條件；計算土體頂部以上存在覆蓋層時，對應(yīng)的頂部溫度邊界條件。