楊小虎， 李少丹， 陳 凱

（武漢第二船舶設(shè)計(jì)研究所 熱能動(dòng)力技術(shù)重點(diǎn)實(shí)驗(yàn)室，武漢 430205）

引 言

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

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

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

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

總結(jié)而言，傳統(tǒng)的準(zhǔn)穩(wěn)態(tài)近似方法實(shí)際上是完全忽略了相變材料顯熱的影響，在近似假設(shè)中只考慮相變潛熱對(duì)傳熱過(guò)程的影響，這必然產(chǎn)生一定的誤差.此外，針對(duì)圓柱坐標(biāo)系中恒熱流邊界條件下的固液相變問(wèn)題的準(zhǔn)穩(wěn)態(tài)近似求解目前還未見(jiàn)文獻(xiàn)報(bào)道.本文將在準(zhǔn)穩(wěn)態(tài)近似方法的基礎(chǔ)上，從熱平衡角度出發(fā)推導(dǎo)給出一維平板和圓柱體外恒熱流熔化問(wèn)題的近似解，通過(guò)這種方法獲得的近似解，比單純利用準(zhǔn)穩(wěn)態(tài)法從數(shù)學(xué)角度推導(dǎo)出來(lái)的近似解具有更高的精度.

本文主要結(jié)構(gòu)如下：首先，簡(jiǎn)單介紹了經(jīng)典的一維平板恒壁溫熔化問(wèn)題，并給出適合工程計(jì)算用的熔化分?jǐn)?shù)和無(wú)量綱換熱擬合關(guān)系式.然后，定量分析了準(zhǔn)穩(wěn)態(tài)近似法的相對(duì)誤差，并基于準(zhǔn)穩(wěn)態(tài)近似方法，從熱平衡角度出發(fā)推導(dǎo)給出恒熱流邊界下熔化問(wèn)題的無(wú)量綱關(guān)系式，并與其他近似解以及數(shù)值解進(jìn)行對(duì)比.針對(duì)圓柱體外的熔化問(wèn)題，本文給出了恒定熱流下的改進(jìn)型準(zhǔn)穩(wěn)態(tài)近似解，并與數(shù)值解進(jìn)行對(duì)比，最后給出了適合工程計(jì)算使用的無(wú)量綱關(guān)系曲線(xiàn)圖.

1 數(shù)學(xué)物理模型

1.1 直角坐標(biāo)系

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

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

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

液相區(qū)域傳熱方程：

1.2 柱坐標(biāo)系

類(lèi)似地，這里結(jié)合圖1(b）給出圓柱體外一維單相熔化問(wèn)題的數(shù)學(xué)描述.

液相區(qū)域傳熱方程：

固液界面條件：

2 無(wú)量綱解析解

2.1 一維平板熔化問(wèn)題

2.1.1 恒壁溫邊界條件

無(wú)量綱化是獲得具有普適性關(guān)系式的重要方法，這里定義問(wèn)題的主要無(wú)量綱參數(shù)如下：

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

可以看出，對(duì)Stefan 問(wèn)題的數(shù)學(xué)求解實(shí)際上只關(guān)心兩個(gè)變量：界面函數(shù)S和溫度分布θ.而在實(shí)際應(yīng)用中，最直接關(guān)心的是熔化體積分?jǐn)?shù)? 和換熱量.在直角坐標(biāo)系下，熔化體積分?jǐn)?shù)為

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

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

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

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

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

圖4 給出了不同Ste下，熔化分?jǐn)?shù)分別為0.2，0.5 和1 時(shí)的液相溫度分布，并與線(xiàn)性化溫度分布對(duì)比.可以看出，Ste越小，溫度越趨近于線(xiàn)性分布.定量地說(shuō)，在完全融化時(shí)刻，線(xiàn)性溫度分布假設(shè)在Ste為0.5，0.3 和0.1 時(shí)的最大誤差分別為5.7%，3.5%和1.2%.這是因?yàn)镾te是相變過(guò)程中，顯熱值與潛熱值的相對(duì)大小的度量.Ste越小，說(shuō)明顯熱相對(duì)于潛熱而言越小，當(dāng)小到顯熱可以忽略時(shí)，從壁面?zhèn)鬟f進(jìn)來(lái)的熱量幾乎全部被相變材料潛熱吸收，而被液相區(qū)域吸收的熱量可以忽略，因此整個(gè)液相區(qū)域的溫度梯度幾乎是相等的，也就是線(xiàn)性溫度分布.從數(shù)學(xué)的角度來(lái)講，就是熱擴(kuò)散方程式（1）中左邊的瞬態(tài)項(xiàng)可以忽略為0，也就是準(zhǔn)穩(wěn)態(tài)近似.

圖4 線(xiàn)性化溫度分布與精確值的對(duì)比Fig. 4 Comparison of approximation solutions and exact solutions

2.1.2 恒熱流邊界條件

上面提到的線(xiàn)性溫度分布假設(shè)實(shí)際上就是傳統(tǒng)的準(zhǔn)穩(wěn)態(tài)近似法的思想，針對(duì)恒定熱流邊界條件，對(duì)應(yīng)的線(xiàn)性溫度分布為

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

下面將定量對(duì)比這一改進(jìn)的準(zhǔn)穩(wěn)態(tài)近似解與其他近似解的精度，包括Evans 等[6]給出的無(wú)窮級(jí)數(shù)解，這里取其前三項(xiàng)：

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

圖5 恒壁溫平板單相熔化問(wèn)題數(shù)值解與精確解對(duì)比Fig. 5 Validation of the numerical method

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

圖6 幾種近似解與數(shù)值解的對(duì)比Fig. 6 Comparison of numerical results and approximation results

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

2.2 圓柱體外熔化問(wèn)題

2.2.1 恒壁溫邊界條件

對(duì)于圓柱體外的單相熔化問(wèn)題，分析方法與前面基本一樣.唯一不同的是，對(duì)于圓柱坐標(biāo)，準(zhǔn)穩(wěn)態(tài)近似溫度分布不再是線(xiàn)性分布，而是對(duì)數(shù)分布.這里，首先簡(jiǎn)單介紹圓柱體外恒壁溫邊界條件單相熔化問(wèn)題的準(zhǔn)穩(wěn)態(tài)近似解.

針對(duì)圓柱體外恒壁溫單相熔化問(wèn)題（圖1（b）），重新定義無(wú)量綱變量：

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

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

圖7 改進(jìn)型準(zhǔn)穩(wěn)態(tài)近似解與數(shù)值解的對(duì)比Fig. 7 Comparison of the improved approximation solutions and numerical results

圖8 改進(jìn)型準(zhǔn)穩(wěn)態(tài)近似解Fig. 8 Improved quasi-steady-state approximation solutions

3 結(jié) 論

參考文獻(xiàn)（ 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］ 李長(zhǎng)玉, 方彥奎, 劉福旭, 等. 熱防護(hù)服-空氣-皮膚熱傳導(dǎo)模型及其解析解[J]. 應(yīng)用數(shù)學(xué)和力學(xué), 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.