一維傳熱方程瞬態問題解析解及其應用

仝 睿，付 浩，宋二祥

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

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

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

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

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

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

1 一維瞬態傳熱問題解析解

1.1 解析解推導

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

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

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

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

對于頂部與底部的固定值部分，存在對應特解：

對于頂部的正弦邊界條件，設對應的溫度解形式為：

根據特解需要滿足的傳熱方程式(1)，解出對應的系數之間的關系：

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

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

假如初始值條件選取為：

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

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

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

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

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

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

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

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

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

其中：

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

1.2 近似解及其討論

式(10)給出的精確解，表達形式較為復雜?？梢詫⑻厥饨馐?8)，作為精確解式(11)的近似。近似解的相對誤差大小為：

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

2 解析解的數值驗證

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

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

3 結合現場試驗數據對解析解進行驗證

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

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

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

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

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

3.1 振幅衰減驗證

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

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

圖5 路基振幅-深度關系圖Fig. 5 Amplitude vs depth diagram of subgrade

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

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

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

圖6 邊坡振幅-深度關系圖Fig. 6 Amplitude vs depth diagram of slope

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

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

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

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

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

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

這一估計數值與式(23)的實際計算數值存在較大的差距。存在差距的原因將在第四節最后進行解釋。

3.2 相位滯后性驗證

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

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

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

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

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

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

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

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

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

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

4 不同性質材料交界面的反射函數推導

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

初始值條件設為：

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

假設上方入射溫度方程為：

假設入射后溫度方程為：

觀察到熱流量：

如果只存在入射現象，則可以得到方程：

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

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

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

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

化簡為：

從而可以解出：

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

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

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

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

這一比率略大于實際值59.48%?？紤]到實際石子與石子之間不是完全貼合，石砟層的傳熱系數應略小于按照空氣、石子并聯估算的式(24)，式(40)得到的比率略高于實際值，屬于較為合理的結果。

5 反射函數的數值模擬驗證

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

初始值條件設為：

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

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

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

6 解析解的應用

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

假設某地區地大氣溫度變化遵循公式：

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

求解T＞273.15 K，得到：

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

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

6.2 使用解析解估計室內試驗的“尺寸效應”

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

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

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

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

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

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

考慮交界面處的反射現象計算得到的日溫度變化，經過瀝青道面衰減為：

年溫度變化衰減為：

對應的，修正后大氣溫度變化條件近似為：

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

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

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

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

其余模擬參數按照文獻[4]進行設置。

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

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

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

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

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

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

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

7 結論

本文針對一維土柱的瞬態問題，求解了對應的解析解。通過數值模擬、現場試驗兩種方式驗證了解析解的正確性。本文還推導了解析解在固定邊界處、不同材料交界面上存在的反射現象。通過對推導結果的分析，發現：

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

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

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