姚鴻驍, 左 沖, 胡小飛, 姚偉岸*

(1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，大連 116024； 2.三一重工集團(tuán)有限公司，武漢 430100)

1 引 言

相變傳熱問題常見于多種工程實(shí)際問題中，如金屬鑄造、焊接、晶體生長、食品冷凍和相變儲(chǔ)能系統(tǒng)[1]等。因?yàn)橄嘧冞^程中存在隨時(shí)間變化的相變界面以及潛熱的吸收或者釋放，因此相變問題屬于非線性問題。只有少數(shù)一維或半無限問題可以獲得解析解，相變問題主要依靠數(shù)值方法求解。

處理移動(dòng)界面和潛熱變化的方法主要有固定域法和界面追蹤法兩種。固定域法將模擬區(qū)域視作一個(gè)整體，基于熱焓、等效熱容或者熱源項(xiàng)建立統(tǒng)一的傳熱控制微分方程，并通過溫度場的插值確定相變界面的位置[2,3]。但由于相變溫度附近參數(shù)存在不連續(xù)性甚至階躍，其非線性方程的模擬存在收斂性問題，尤其是相變潛熱數(shù)值較大時(shí)。界面追蹤法則不需要處理這種不連續(xù)性，而是分別在固相和液相兩個(gè)區(qū)域中求解瞬態(tài)傳熱問題，然后根據(jù)相變界面上的能量連續(xù)條件追蹤界面的移動(dòng)[4]，但界面追蹤法需要連續(xù)地對離散網(wǎng)格進(jìn)行重劃分。

空間離散的數(shù)值方法主要有有限元法、有限差分法、有限體積法、無網(wǎng)格法和邊界元法。其中邊界元法已經(jīng)成為求解傳熱問題的有效方法之一，如周楓林等[5]使用邊界元法分析了混凝土水壩的瞬態(tài)熱傳導(dǎo)問題；高效偉等[6]提出了一種三維等參管單元邊界元法，并將其應(yīng)用于熱傳導(dǎo)分析中；趙金軍等[7]提出了一種計(jì)算幾乎奇異邊界積分的邊界元方法并求解了三維非均質(zhì)熱傳導(dǎo)問題。相對其他方法,邊界元法的主要優(yōu)點(diǎn)是只需離散求解域的邊界。因此，邊界元法非常適合于求解移動(dòng)邊界問題和邊界未知的反問題，如左沖等[8]利用時(shí)域徑向積分邊界元法和界面追蹤法模擬了平面單相凝固問題；Wang等[9]將邊界元法與水平集法相結(jié)合模擬了具有復(fù)雜拓?fù)渥兓南嘧儐栴}；Yu等[10]提出了一種基于邊界元法的直接反演方法來辨識內(nèi)部空腔的形狀。

本文采用精細(xì)積分邊界元法和界面追蹤法求解兩相等溫相變問題。首先，分別在固液兩相中利用位勢問題的基本解建立積分方程，并采用徑向積分法[11]將其中的域積分轉(zhuǎn)化為邊界積分。然后，利用邊界元法和精細(xì)積分法[12]分別求得固液兩相的溫度場和邊界熱流密度。最后，基于相變界面上的能量平衡方程,采用界面追蹤法得到移動(dòng)界面的位置。

2 控制方程

如圖1所示的二維相變問題，相變界面Γif將整個(gè)區(qū)域Ωtotal分為固相Ωs和液相Ωl兩個(gè)區(qū)域。如忽略密度變化，則傳熱控制微分方程為

(x∈Ωs)(1)

(x∈Ωl)(2)

式中x=(x1,x2)，=?/?x1i+?/?x2j,T(x,t)為t時(shí)刻點(diǎn)x的溫度，c,k和分別為比熱容、熱導(dǎo)率和密度，f(x,t)為域內(nèi)熱源。

圖1 相變問題Fig.1 Phase change problem

邊界條件為

(x∈Γ1)

(3)

(x∈Γ2)

(4)

式中n為單位外法向量，q為熱流密度。

初始條件為

T(x,0)=T0(x)

(t=0)

(5)

相變界面Γif的溫度條件和能量平衡條件

T(x,t)=Tm

(x∈Γif)

(6)

(x∈Γif)

(7)

式中Tm為相變溫度，L為相變潛熱，Vn為相變界面的法向移動(dòng)速度，其由固液兩相的熱流得到。

分別在固相域和液相域中應(yīng)用徑向積分邊界元法求解，省略下標(biāo)s和l，固液兩相的控制方程可統(tǒng)一寫為

(8)

3 精細(xì)積分邊界元法

基于式(8)的積分方程為

(9)

(10)

式中r(x,y)為源點(diǎn)y到場點(diǎn)x的距離。

需要利用徑向積分法[11]將式(9)的域積分轉(zhuǎn)化為邊界積分，其中?T(x,t)/?t為未知量，可用徑向基函數(shù)和坐標(biāo)分量的多項(xiàng)式將其近似表示為

(11)

(12)

(13)

式中Dy為N維橫向量，其第j個(gè)分量為

(14)

(15)

(16)

(17)

將邊界Γ離散為Ne個(gè)邊界單元和Nb個(gè)邊界點(diǎn)，域內(nèi)配置Ni個(gè)內(nèi)部點(diǎn)，N=Nb+Ni。源點(diǎn)y取遍所有邊界點(diǎn)和內(nèi)部點(diǎn)，由式(12)得出式(18)，即

(18)

(19)

然后，利用精細(xì)積分法[12]求解式(19)，如果r在時(shí)間步[tk,tk + 1]內(nèi)是關(guān)于時(shí)間的線性函數(shù)，則式(19)的求解如下。

Xk + 1=E[Xk+B-1r0+r1]-

B-1(r0+Δtr1)-r1

(20)

式中 Δt=tk + 1-tk，E=exp(BΔt)

r0=r(tk)，r1=[r(tk + 1)-r0]/Δt

4 界面追蹤法

為了追蹤相變界面的移動(dòng)位置，界面追蹤法首先需要確定界面的移動(dòng)速度和方向。界面移動(dòng)速度由式(7)計(jì)算得到，而界面移動(dòng)方向由式(21)計(jì)算[14]

(21)

式中 下標(biāo)j為移動(dòng)界面上的j號節(jié)點(diǎn)，下標(biāo)i-1和i為j號節(jié)點(diǎn)相鄰的兩個(gè)單元，l為單元長度。

時(shí)間步[tk,tk + 1]內(nèi)界面追蹤法的步驟如下。

(1) 指定一個(gè)界面最大移動(dòng)距離Δs。

(2) 根據(jù)初始條件或者上一時(shí)間步數(shù)值結(jié)果得到tk時(shí)刻固相域和液相域的邊界熱流qs(tk)和ql(tk)，跟據(jù)式(7)計(jì)算速度V(tk)={Vj(tk)}。

(3) 利用式(22)計(jì)算時(shí)間步長[15]

(22)

式中 ‖·‖∞為向量的無窮范數(shù)。

(4) 預(yù)測tk + 1時(shí)刻的界面移動(dòng)速度為Vp(tk + 1)= 0.7V(tk)。

(5) 利用式(23)計(jì)算移動(dòng)界面上節(jié)點(diǎn)預(yù)測位置

(23)

(6) 更新tk + 1時(shí)刻的網(wǎng)格，在固相域和液相域分別利用精細(xì)積分邊界元法求解未知溫度場和熱流密度qs(tk + 1)和ql(tk + 1)，再依據(jù)式(7)計(jì)算速度V(tk + 1)={Vj(tk + 1)}。

(7) 收斂性檢測，如收斂條件滿足

(24)

則進(jìn)入步驟(8)，否則令Vp(tk + 1)=V(tk + 1)，回到步驟(5)。

(8) 步驟(5)得到的移動(dòng)界面位置和步驟(6)得到的溫度場即為tk + 1時(shí)刻的計(jì)算結(jié)果，對移動(dòng)界面進(jìn)行網(wǎng)格重劃分，并根據(jù)需要增加或者刪減內(nèi)部節(jié)點(diǎn)。令tk=tk + 1，回到步驟(1)。

5 數(shù)值算例

利用四個(gè)數(shù)值算例來驗(yàn)證本文提出的數(shù)值方法。數(shù)值計(jì)算通過Matlab編程實(shí)現(xiàn)，界面追蹤法終止迭代的容許差限統(tǒng)一為ε=0.001。

算例1考慮一個(gè)一維半無限凝固相變問題。初始溫度為T0=2 ℃，x1=0處邊界的溫度保持為Tw=-10 ℃，其他邊界保持絕熱。相變材料的熱物理參數(shù)列入表1。該問題可用圖2所示的二維1×0.1(m2)矩形區(qū)域的相變問題近似。

表1 相變材料熱物理參數(shù)Tab.1 Thermo -physical properties of the phase change materials

圖2 算例1的邊界元計(jì)算模型Fig.2 BEM model for Example 1

在計(jì)算過程中，相變界面每移動(dòng)0.02 m,則在固相域增加2個(gè)邊界節(jié)點(diǎn)和3個(gè)內(nèi)部節(jié)點(diǎn)，相應(yīng)在液相域減少2個(gè)邊界節(jié)點(diǎn)和3個(gè)內(nèi)部節(jié)點(diǎn)。圖3給出了不同時(shí)刻相變界面位置與解析解[16]對比的相對誤差，吻合良好?？梢钥闯觯分別取0.004 m，0.006 m和0.008 m時(shí),最大誤差分別是0.63%，0.89%和1.15%，計(jì)算結(jié)果均具有很高精度。

圖3 算例1的界面位置相對誤差Fig.3 Relative errors of interface position for Example 1

算例2考慮一個(gè)無限角域凝固問題。初始溫度為T0=0.3 ℃，x1=0和x2=0處邊界的溫度保持為Tw=-1 ℃，其他邊界保持絕熱。相變材料的熱物理參數(shù)列入表1。該問題可用圖4所示的一個(gè)二維1×1(m2)正方形區(qū)域的凝固相變問題來近似。

圖4 算例2的邊界元計(jì)算模型Fig.4 BEM model for Example 2

表2列出了Δs=0.02 m時(shí)，采用本文方法模擬得到的不同時(shí)刻相變界面與x1=1和x1=x2兩條直線交點(diǎn)的x2坐標(biāo)值以及對應(yīng)的解析解[17]，其中括號內(nèi)為相對誤差，最大相對誤差分別為 0.96% 和0.44%。雖然界面追蹤法需要迭代求解相變界面的位置，但本文方法能夠快速收斂得到結(jié)果，如對算例2，每個(gè)時(shí)間步的平均迭代次數(shù)為2.9，最大次數(shù)為4。

表2 算例2相變界面位置與誤差Tab.2 Interface positions and relative errors for Example 2

雖然銅管域內(nèi)不發(fā)生相變，但需對銅管進(jìn)行瞬態(tài)熱傳導(dǎo)問題數(shù)值分析。分別在銅管域和相鄰的相變材料域應(yīng)用邊界元法，得到兩組微分方程?；阢~管-相變材料界面上的溫度連續(xù)條件和熱流平衡條件，這兩組微分方程可聯(lián)立形成一組新微分方程組。利用精細(xì)積分法求解得到溫度場和邊界熱流。

圖5 算例3的邊界元計(jì)算模型Fig.5 BEM model for Example 3

基于對稱性，選取裝置截面的1/8進(jìn)行數(shù)值計(jì)算。圖6給出了不同時(shí)刻相變界面的移動(dòng)位置，同時(shí)圖7給出了固相分?jǐn)?shù)的變化。以ANSYS軟件的數(shù)值模擬結(jié)果作為參考解，對比結(jié)果表明,本方法具有很高的模擬精度。

圖6 算例3的相變界面位置Fig.6 Phase change interface position for Example 3

圖7 算例3的固相分?jǐn)?shù)變化Fig.7 Evolution of solid fraction for Example 3

算例4考慮如圖8所示區(qū)域的相變問題。初始溫度為T0=-10 ℃，外邊界的溫度保持為Tw=25 ℃。相變材料的熱物理參數(shù)列入表1。

圖8 算例4的邊界元計(jì)算模型Fig.8 BEM model for Example 4

圖9給出了不同時(shí)刻相變界面的移動(dòng)位置，并與ANSYS軟件數(shù)值結(jié)果進(jìn)行了對比?？梢钥闯鱿嘧儏^(qū)域發(fā)生了拓?fù)渥兓疚姆椒ㄈ匀荒軌蛴行揖_地模擬此類相變問題。

圖9 算例4的相變界面位置Fig.9 Phase change interface position for Example 4

6 結(jié) 論

本文將精細(xì)積分邊界元法與界面追蹤法相結(jié)合求解二維雙相相變問題。數(shù)值算例的結(jié)果表明，本文提出的方法能夠有效地模擬相變問題，并且具有精度高和迭代次數(shù)少的優(yōu)點(diǎn)。對于相變區(qū)域發(fā)生拓?fù)渥兓膯栴}本文方法仍然是有效的，能夠得到精確的模擬結(jié)果。當(dāng)相變界面最大移動(dòng)距離Δs取不同數(shù)值時(shí)，本文方法依然能夠得到穩(wěn)定的數(shù)值解，表明本文方法具有很好的數(shù)值穩(wěn)定性。