第二類邊界條件下固液相變問題的改進型準穩態近似解*

楊小虎， 李少丹， 陳 凱

（武漢第二船舶設計研究所 熱能動力技術重點實驗室，武漢 430205）

引 言

固液相變現象廣泛存在于自然界和工業生產當中，如水（冰）的凝固（熔化）、冶金工業、熱能（冷量）存儲與利用[1]、電子設備或人體熱防護[2]以及生物醫學等領域[3].對固液相變傳熱問題的理論研究始于19 世紀奧地利物理學家Stefan 對極地冰熔化過程和土壤凍結問題的研究，故又稱為Stefan 問題.Stefan 問題是兩相耦合的具有移動邊界的非線性問題，且一般為瞬態問題，這使得其理論求解十分困難.目前，僅針對少數簡單情況下的Stefan 問題存在解析解，最典型的是Neumann 于1864 年首先給出的直角坐標系半無限大空間第一類邊界條件（恒壁溫）下一維Stefan 問題的精確解[4].

作為經典Stefan 問題的延伸，當邊界條件由第一類邊界條件變化為第二類邊界條件（恒熱流）時，其理論求解將變得十分復雜[5].Evans 等[6]采用Laplace 積分變換，并假定固液界面函數是時間的冪級數的形式，首次給出了恒熱流邊界條件下單相Stefan 問題的界面函數解析解表達式.Schiavone 等[7]進一步發展了這一方法并完整地給出了界面函數和溫度分布的表達式.然而，這種冪級數形式的解是無界函數，且必須無限求和才能得到真正收斂的解，這在實際計算中不太適用.El-Genk 等[8]于1979 年以微分方程的形式給出了這一問題的解析解，但其不便于直接使用，而且其正確性目前也存在爭議.

此外，當求解問題由直角坐標系變換到柱坐標系或球坐標系時，理論求解會變得更加困難.例如，針對工程中常見的圓柱體內（外）的相變問題，目前還沒有獲得精確解，實際處理中往往使用數值模擬計算.盡管數值模擬方法已經很成熟，但是使用起來比較麻煩，而且難以獲得通用的無量綱解，不利于實際問題的快速分析.因此，給出形式簡單的通用近似解析解對于實際問題的初步分析十分必要[9-11].目前，針對第二類邊界條件下Stefan 問題的近似求解方法很多，主要包括熱平衡積分近似法、攝動法、準穩態法等[4].熱平衡積分法數學求解過程復雜，特別是在柱坐標下，且對很多問題給出的解是微分方程組的形式，不便于實際使用.攝動法適用于Stefan 數Ste很小的情況，在Ste較大時精度較差，因此這兩種方法在實際使用時均有所限制[9].

準穩態法最早由Stefan 提出，并用于大潛熱固液相變問題的近似求解[4].Crank[12]利用這一方法求解了直角坐標系下邊界條件為?T/?x=C?T/?t（C為常數）的一維相變擴散問題，并推廣到圓柱和球坐標下.Lin 等[13]利用直角坐標系下的Neumann 精確解對準穩態法進行修正，使其具有更高的精度，并將這一修正方法推廣到柱坐標和球坐標恒壁溫邊界條件問題的求解.Megerlin 利用準穩態近似法求解了恒熱流邊界條件下直角坐標系下的熔化問題[4].Solomon[14]對Megerlin 的方法做出了改進，改進方法形式上與Goodman[15]提出的積分近似法一樣，只不過不需要滿足熱平衡條件.上述的所有準穩態近似方法均是從數學角度出發來進行分析求解，通過假定一個簡單形式的溫度分布函數并讓此函數滿足邊界條件和Stefan 界面條件，從而獲得界面函數和溫度分布的表達式.

總結而言，傳統的準穩態近似方法實際上是完全忽略了相變材料顯熱的影響，在近似假設中只考慮相變潛熱對傳熱過程的影響，這必然產生一定的誤差.此外，針對圓柱坐標系中恒熱流邊界條件下的固液相變問題的準穩態近似求解目前還未見文獻報道.本文將在準穩態近似方法的基礎上，從熱平衡角度出發推導給出一維平板和圓柱體外恒熱流熔化問題的近似解，通過這種方法獲得的近似解，比單純利用準穩態法從數學角度推導出來的近似解具有更高的精度.

本文主要結構如下：首先，簡單介紹了經典的一維平板恒壁溫熔化問題，并給出適合工程計算用的熔化分數和無量綱換熱擬合關系式.然后，定量分析了準穩態近似法的相對誤差，并基于準穩態近似方法，從熱平衡角度出發推導給出恒熱流邊界下熔化問題的無量綱關系式，并與其他近似解以及數值解進行對比.針對圓柱體外的熔化問題，本文給出了恒定熱流下的改進型準穩態近似解，并與數值解進行對比，最后給出了適合工程計算使用的無量綱關系曲線圖.

1 數學物理模型

1.1 直角坐標系

固液相變問題可分為凝固和熔化兩類，其數學描述和求解過程是類似的，本文以熔化過程為例進行說明.根據固相區域是否存在溫度梯度，熔化問題又可以分為單相熔化問題和兩相熔化問題.熔化過程中主要的熱量是以相變潛熱的形式被吸收的，因此這里不考慮固相區域初始過冷的情形，而僅分析初始溫度為相變溫度的情況，也就是單相熔化問題.

在直角坐標系下，單相熔化問題也就是一維平板的熔化問題，其物理模型見圖1(a）.初始時整個區域均處于初始溫度Ti（Ti=Tm），Tm為相變溫度.t＞0 時，邊界x=0 處的溫度突然升高到Tw（Tw＞Tm）或者給定恒定加熱熱流q′′，熔化過程隨即開始，固液界面s(t)由左向右移動.由于固相區域始終處于相變溫度，故不需要求解，這里僅對液相區域進行求解.如無特殊說明，下文中所有變量和參數均默認為是液相的.假定相變材料的物性參數保持為恒定值，該問題的數學描述如下.

圖1 傳導型固液相變物理模型：（a）平板單相熔化問題；（b）圓柱體外單相熔化問題Fig. 1 Schematics of single phase Stefan problems: (a) the semi-finite slab; (b) the cylinder

液相區域傳熱方程：

1.2 柱坐標系

類似地，這里結合圖1(b）給出圓柱體外一維單相熔化問題的數學描述.

液相區域傳熱方程：

固液界面條件：

2 無量綱解析解

2.1 一維平板熔化問題

2.1.1 恒壁溫邊界條件

無量綱化是獲得具有普適性關系式的重要方法，這里定義問題的主要無量綱參數如下：

其中，ε 為誤差函數；λ 是熔化常數，由超越方程式（16）決定.由式（16）可以看出，熔化常數實際上僅僅是Ste的函數，這里將其定義為λ=f(Ste).

可以看出，對Stefan 問題的數學求解實際上只關心兩個變量：界面函數S和溫度分布θ.而在實際應用中，最直接關心的是熔化體積分數? 和換熱量.在直角坐標系下，熔化體積分數為

式（18）和（21）給出了表征一維平板單相熔化過程的最主要的無量綱關系式，現在的問題是如何獲得f(Ste)和g(Ste)的表達式.由于f(Ste)=λ 是由超越方程式（16）決定的，理論上講是無法獲得其解析表達式的.而顯式表達式顯然是更便于實際使用的，特別是在參數化設計中.Carslaw 等[16]對誤差函數ε(λ)的級數展開取一階近似，得到了f(Ste)的近似表達式.然而，這一近似表達式僅對小的Ste（Ste＜0.1）成立，當Ste較大時，誤差較大.因此，這里通過擬合關系式來獲得f(Ste)和g(Ste)的表達式.

這里以設備熱管理中常用的典型的石蠟相變材料n-eicosane[17-18]為例，可計算得其在50 ℃過熱壁面溫度下的Ste為0.47.對于低熔點金屬相變材料，其Ste一般比同等溫差下的石蠟小.對于大多數電子設備熱管理問題，溫差一般不會超過50 ℃.因此，可以認為在這一類問題中Ste≤0.5，這里僅在此范圍內對f(Ste)和g(Ste)進行擬合.

圖2 f(Ste)和g(Ste)擬合曲線Fig. 2 Fitting curves of f(Ste) and g(Ste)

此外，由式（18）和（21）可以得出，Nu與? 的乘積僅僅是Ste的函數.圖3 給出了Nu·?-Ste變化曲線，可以看出Nu·?≈1.且Ste越小時，Nu·? 越小且越接近于1.這一結論實際上可以通過一個簡單的假設得到.假定在液相區域內溫度成線性分布，即θ=1-X/S(τ)，則

圖3 Nu·? 隨Ste 的變化Fig. 3 Variation of Nu·? with Ste

圖4 給出了不同Ste下，熔化分數分別為0.2，0.5 和1 時的液相溫度分布，并與線性化溫度分布對比.可以看出，Ste越小，溫度越趨近于線性分布.定量地說，在完全融化時刻，線性溫度分布假設在Ste為0.5，0.3 和0.1 時的最大誤差分別為5.7%，3.5%和1.2%.這是因為Ste是相變過程中，顯熱值與潛熱值的相對大小的度量.Ste越小，說明顯熱相對于潛熱而言越小，當小到顯熱可以忽略時，從壁面傳遞進來的熱量幾乎全部被相變材料潛熱吸收，而被液相區域吸收的熱量可以忽略，因此整個液相區域的溫度梯度幾乎是相等的，也就是線性溫度分布.從數學的角度來講，就是熱擴散方程式（1）中左邊的瞬態項可以忽略為0，也就是準穩態近似.

圖4 線性化溫度分布與精確值的對比Fig. 4 Comparison of approximation solutions and exact solutions

2.1.2 恒熱流邊界條件

上面提到的線性溫度分布假設實際上就是傳統的準穩態近似法的思想，針對恒定熱流邊界條件，對應的線性溫度分布為

將式（31）代入式（29）可以得到

下面將定量對比這一改進的準穩態近似解與其他近似解的精度，包括Evans 等[6]給出的無窮級數解，這里取其前三項：

由于針對這一問題目前還沒有公認的解析解，這里給出數值解作為參考.數值解由商業軟件FLUENT 求解獲得，數值算法采用經典的焓方法，并且忽略液相區域的自然對流.計算網格和時間步長的選擇也經過了無關性驗證.作為數值算法有效性的驗證，圖5 給出了恒壁溫條件下平板單相融化問題的數值解（Ste=0.47），并與理論解做對比.可以看出兩者吻合得非常好，說明數值解具有足夠的精確度.下面所有的數值計算將采用同樣的方法.

圖5 恒壁溫平板單相熔化問題數值解與精確解對比Fig. 5 Validation of the numerical method

圖6 對比了幾種不同近似解的精度，分別給出了Ste為0.1，0.3 和0.5 時固液界面位置S隨無量綱時間Fo的變化曲線.以數值計算結果作為精確值的基準參考，定量對比不同近似解在S等于1 時對應的Fo數的大小相對于精度值的偏差，見表1.不難看出，Goodman[15]的積分近似解與數值計算結果吻合得最好，相對誤差在1%以內.Evans 等[6]的級數解（取前三項）在Ste較大時，誤差較大，達到-13.8%.El-Genk 等[8]給出的解在所有Ste情況下都與數值解偏離較多，在9% ～ 10%之間.傳統的不考慮顯熱的準穩態近似解偏差也較大，在Ste從0.1 增加到0.5 時，其計算誤差從-3.6%增加到-16.9%.本文提出的改進型準穩態結果與數值解吻合較好，其精度僅次于Goodman[15]的積分近似解，在Ste從0.1 增加到0.5 時，其計算誤差在1.2% ～ 3.8%之間.與積分近似解（式（38））不同的是，本文給出的是界面函數關于Ste和Fo數的顯式表達式，比隱式表達式使用起來更加方便.

圖6 幾種近似解與數值解的對比Fig. 6 Comparison of numerical results and approximation results

表1 幾種近似解誤差對比（以S=1 時的Fo 作為對比指標）Table 1 Comparison of different approximate solutions (Fo as the index for S=1)

2.2 圓柱體外熔化問題

2.2.1 恒壁溫邊界條件

對于圓柱體外的單相熔化問題，分析方法與前面基本一樣.唯一不同的是，對于圓柱坐標，準穩態近似溫度分布不再是線性分布，而是對數分布.這里，首先簡單介紹圓柱體外恒壁溫邊界條件單相熔化問題的準穩態近似解.

針對圓柱體外恒壁溫單相熔化問題（圖1（b）），重新定義無量綱變量：

與式（41）不同的是，此時的壁面溫度Tw(t)不再是常數，而是時間的函數，且滿足邊界條件：

圖7 將式（48）與數值解進行對比，可以看出，改進型準穩態近似解與數值解十分接近，Ste越小，近似解與數值解的相對誤差越小.由于式（48）是關于S的超越方程，無法獲得解析解，這里以曲線圖的形式給出其在不同Ste下隨無量綱時間Fo數變化的曲線，以供實際計算參考，如圖8.具體使用時，由外徑內徑比確定最大的Smax=γ，并在在圖上畫水平線，與相應的Ste的曲線相交，該交點對應的橫坐標即為完全熔化時間，對應的曲線即為熔化曲線，如圖7 中紅色曲線（γ=2，Ste=0.3）對應的界面函數為圖8 中的紅色曲線部分.

圖7 改進型準穩態近似解與數值解的對比Fig. 7 Comparison of the improved approximation solutions and numerical results

圖8 改進型準穩態近似解Fig. 8 Improved quasi-steady-state approximation solutions

3 結 論

參考文獻（ References ） ：

［1］ZHAO B C, LI T X, GAO J C, et al. Latent heat thermal storage using salt hydrates for distributed building heating: a multi-level scale-up research[J].Renewable and Sustainable Energy Reviews, 2020, 121（7）: 109712.

［2］ 李長玉, 方彥奎, 劉福旭, 等. 熱防護服-空氣-皮膚熱傳導模型及其解析解[J]. 應用數學和力學, 2021, 42（2）: 162-169.（LI Changyu, FANG Yankui, LIU Fuxu, et al. A thermal protective clothing-air-skin heat conduction model and its analytical solution[J].Applied Mathematics and Mechanics, 2021, 42（2）: 162-169.（in Chinese））

［3］FLEISCHER A S.Thermal Energy Storage Using Phase Change Materials: Fundamentals and Applications[M].Berlin: Springer, 2015.

［4］CRANK J.Free and Moving Boundary Problems[M]. Oxford: Clarendon Press, 1984.

［5］TAO L. On free boundary problems with arbitrary initial and flux conditions[J].Zeitschrift für Angewandte Mathematik und Physik(ZAMP), 1979, 30（3）: 416-426.

［6］EVANS G, ISAACSON I, MACDONALD J. Stefan-like problems[J].Quarterly of Applied Mathematics, 1950,8（3）: 312-319.

［7］SCHIAVONE P, CONSTANDA C, MIODUCHOWSKI A.Integral Methods in Science and Engineering[M].Berlin: Springer Science & Business Media, 2012.

［8］EL-GENK M S, CRONENBERG A W. Solidification in a semi-infinite region with boundary conditions of the second kind: an exact solution[J].Letters in Heat and Mass Transfer, 1979, 6（4）: 321-327.

［9］QU P, ZHANG C, LIAO X, et al. A new computation method for solidification process in a finite, initially overheated slab[J].Journal of Thermal Science, 1992, 1（4）: 272-277.

［10］BELI G E. Solidification of a liquid about a cylindrical pipe[J].International Journal of Heat and Mass Transfer, 1979, 22（12）: 1681-1686.

［11］YANG X H, LIU J. A novel method for determining the melting point, fusion latent heat, specific heat capacity and thermal conductivity of phase change materials[J].International Journal of Heat and Mass Transfer, 2018,127（Part B）: 457-468.

［12］CRANK J. Diffusion with rapid irreversible immobilization[J].Transactions of the Faraday Society, 1957, 53:1083-1091.

［13］LIN S, JIANG Z. An improved quasi-steady analysis for solving freezing problems in a plate, a cylinder and a sphere[J].Journal of Heat Transfer, 2003, 125（6）: 1123-1128.

［14］SOLOMON A D. The applicability and extendability of Megerlin’s method for solving parabolic free boundary problems[J].Moving Boundary Problems, 1977, 1: 187-202.

［15］GOODMAN T R. The heat-balance integral and its application to problems involving a change of phase[J].Journal of Fluids Engineering, 1958, 80（2）: 335-342.

［16］CARSLAW H S, JAEGER J C.Conduction of Heat in Solids[M]. Oxford: Oxford Clarendon Press, 1959.

［17］RUKH S, PASHA R A, NASIR M A. Heat transfer enhancement of round pin heat sinks usingn-eicosane as PCM: an experimental study[J].Heat and Mass Transfer, 2019, 55（2）: 309-325.

［18］ZENG Y, DONG J, KHODADADI J M. Thermal coupling-decoupling mechanism of heat transfer across van der Waals interfaces inn-eicosane[J].International Journal of Heat and Mass Transfer, 2021, 164（10）: 120603.