劉鴻瑨，胡欲立，楊 威，郝澤花

(西北工業(yè)大學(xué) 航海學(xué)院，陜西 西安，710072)

0 引言

我國(guó)海洋資源十分豐富，人們對(duì)海洋資源的探索開(kāi)發(fā)離不開(kāi)各種水下航行器[1]。水下熱滑翔機(jī)是一種新型無(wú)人自主水下航行器，其利用相變材料在固液相變過(guò)程中產(chǎn)生的體積變化，將海洋溫差能轉(zhuǎn)化為機(jī)械能，并完成水下的之字形運(yùn)動(dòng)。

國(guó)內(nèi)學(xué)者很早就對(duì)水下熱滑翔機(jī)展開(kāi)了研究，孔巧玲[2]、倪園芳[3]、任龍飛[4]等建立一維的水下熱滑翔機(jī)熱管換熱模型，采用Matlab 軟件，利用有限容積法對(duì)滑翔機(jī)上相變材料的固液相變過(guò)程換熱特性展開(kāi)了研究，計(jì)算過(guò)程中利用等效導(dǎo)熱系數(shù)的方法來(lái)表現(xiàn)凝固融化過(guò)程中自然對(duì)流作用的影響，計(jì)算過(guò)程并未擴(kuò)展到二維或三維情況。王延輝[5]和Yang 等[6]對(duì)溫差能驅(qū)動(dòng)的水下滑翔機(jī)展開(kāi)了動(dòng)力學(xué)分析與設(shè)計(jì)，并設(shè)計(jì)出了一種帶內(nèi)柔囊的換熱熱管，但是并未對(duì)熱管內(nèi)相變換熱過(guò)程溫度場(chǎng)及液相率的分布進(jìn)行深入研究。梁澤德等[7]使用正十六烷和石墨為9：1 組成的復(fù)合相變材料，研究對(duì)比了不同尺寸圓筒形和圓柱形豎直換熱容器在不同溫度下的相變換熱特性，但是并未考慮到壓力對(duì)相變過(guò)程的影響。李國(guó)道[8]以石蠟為研究對(duì)象，采用添加膨脹石墨的方法來(lái)提高石蠟換熱性能，并仿真計(jì)算了水平圓管內(nèi)的相變過(guò)程。楊銳[9]以海洋浮標(biāo)上低溫固液相變材料為研究對(duì)象，采用實(shí)驗(yàn)的方法研究了添加改進(jìn)劑對(duì)正十六烷相變溫度的影響。從前人研究?jī)?nèi)容來(lái)看，對(duì)水下熱滑翔機(jī)相變過(guò)程研究主要集中在數(shù)值仿真、熱管結(jié)構(gòu)以及相變材料幾方面。

在此基礎(chǔ)上，文中從融化溫度的角度，引入了壓力對(duì)相變過(guò)程的影響。建立了二維水下熱滑翔機(jī)換熱熱管中的相變換熱過(guò)程仿真模型。為方便模型建立以及網(wǎng)格繪制，同時(shí)也為得到自然對(duì)流作用對(duì)相界面的影響，使用FLUENT 軟件對(duì)水下熱滑翔機(jī)熱管內(nèi)相變過(guò)程溫度及液相率的變化進(jìn)行了仿真分析，對(duì)水下熱滑翔機(jī)溫差能動(dòng)力系統(tǒng)設(shè)計(jì)改進(jìn)具有積極意義。

1 相變換熱模型

1.1 物理模型

圖1 為熱管截面相變過(guò)程示意圖，管內(nèi)部存放相變材料，整個(gè)圓筒處于海水環(huán)境中。材料初始狀態(tài)為液態(tài)，當(dāng)管外海水低于材料相變溫度時(shí)，相變材料開(kāi)始凝固；當(dāng)管外海水溫度高于材料相變溫度時(shí)，管內(nèi)相變材料開(kāi)始融化。熱滑翔機(jī)所用相變材料為正十六烷，其物性參數(shù)受溫度影響較小。

為簡(jiǎn)化計(jì)算模型進(jìn)行如下處理：

圖1 熱管結(jié)構(gòu)示意圖Fig.1 Schematic diagram of heat pipe structure

1) 假設(shè)固液相變材料正十六烷各向同性，固液兩相的物性參數(shù)均為常數(shù)，不隨溫度發(fā)生變化，比如熱導(dǎo)熱系數(shù)、粘度等；

2) 相較于相變材料，管壁材料的導(dǎo)熱系數(shù)很大，管壁厚度可忽略；

3) 液態(tài)工質(zhì)溫度分布不均會(huì)使其密度分布不均勻，進(jìn)而產(chǎn)生自然對(duì)流作用，在仿真計(jì)算過(guò)程中密度選項(xiàng)選擇布辛奈斯克(Boussinesq)近似，其作用是在動(dòng)量方程中浮力項(xiàng)中，將原本為常數(shù)的密度改為與溫度有關(guān)的函數(shù)。

由于水下熱滑翔機(jī)熱管為細(xì)長(zhǎng)圓管，溫度隨管長(zhǎng)變化不明顯，故截取橫截面分析，使用二維結(jié)構(gòu)網(wǎng)格。網(wǎng)格用ICEM CFD 繪制，圓面直徑為24 mm。經(jīng)網(wǎng)格無(wú)關(guān)性檢驗(yàn)，取網(wǎng)格單元數(shù)量9 396，節(jié)點(diǎn)數(shù)量9 289 時(shí)，計(jì)算能保證較高的精度和較快的計(jì)算速度。將網(wǎng)格導(dǎo)入FLUENT 18.0，選擇二維分離隱式雙精度求解器，并將圓截面中心點(diǎn)作為觀(guān)測(cè)節(jié)點(diǎn)，監(jiān)測(cè)該點(diǎn)溫度變化。熱滑翔機(jī)在水下以緩慢速度滑翔，其熱管與周?chē)Ｋg為對(duì)流換熱，假設(shè)滑翔機(jī)運(yùn)行速度為0.28 m/s[10]，經(jīng)計(jì)算取對(duì)流換熱系數(shù)為324 W/(m2·K)，仿真計(jì)算時(shí)間步長(zhǎng)取0.1 s。

正十六烷在固液相變過(guò)程體積變化大，且固態(tài)密度大于液態(tài)密度，其融點(diǎn)剛好在海洋溫躍層溫差范圍內(nèi)(4～25℃)，比較適合作為溫差能熱機(jī)的相變工質(zhì)。文中采用正十六烷為相變材料，其物性參數(shù)如表1 所示。

1.2 數(shù)學(xué)模型

FLUENT 軟件中的Solidification/Melting 模型[11]可以用來(lái)仿真材料相變過(guò)程，模型中能量方程采用焓法[12]求解，在固態(tài)、液態(tài)和固液糊狀區(qū)建立統(tǒng)一的能量方程，不需要跟蹤相界面的移動(dòng)，以簡(jiǎn)化計(jì)算過(guò)程。

表1 正十六烷物性參數(shù)Table 1 Physical parameters of n-hexadecane

式中：H為單位體積相變材料的總焓；ρ為材料密度；υ為流體速度；S為動(dòng)量源項(xiàng)；k為導(dǎo)熱系數(shù)；ν為運(yùn)動(dòng)粘度；p為壓強(qiáng)。

焓孔隙率技術(shù)將糊狀區(qū)域(部分凝固區(qū)域)作為多孔介質(zhì)，每個(gè)微元內(nèi)的孔隙度等于該孔內(nèi)的液體分?jǐn)?shù)。在動(dòng)量方程添加動(dòng)量源項(xiàng)

式中：f為液相率，在液體區(qū)f=0，在固體區(qū)f=1，在糊狀區(qū)域f在0～1 之間變化；ε為常數(shù)，作用是防止仿真計(jì)算過(guò)程中式(4)分母變成為0，通常取一個(gè)極小數(shù)值0.001；Amush為糊狀區(qū)常數(shù)，一般取1×105。

FLUENT Solidification/Melting 模型中液相率f的定義如下

對(duì)于純物質(zhì)，其融化溫度lT等于其凝固溫度sT，但在計(jì)算過(guò)程中為得到較好的收斂性，2 個(gè)溫度之間通常取一個(gè)極小的差值，文中在后續(xù)計(jì)算過(guò)程中差值取0.5 K。

焓法求解能量方程的過(guò)程中，液相率是用溫度定義的，故須建立液相率和工質(zhì)焓的關(guān)系。在計(jì)算過(guò)程中，工質(zhì)焓定義為

式中：L為材料相變潛熱；href為參考焓值；Tref為參考溫度；cp為相變材料的比定壓熱容。

根據(jù)以上公式建立焓和溫度的關(guān)系。對(duì)于相變過(guò)程中材料的導(dǎo)熱系數(shù)采用如下表達(dá)形式

式中，ks、kl分別為材料固相和液相的導(dǎo)熱系數(shù)。

材料比熱容示意見(jiàn)圖2，具體定義如下

圖2 不同區(qū)域比熱容示意圖Fig.2 Schematic diagram of specific heat capacity of different regions

1.3 壓力對(duì)相變過(guò)程的影響

熱滑翔機(jī)相變材料在融化時(shí)體積膨脹會(huì)壓縮蓄能器中的氣體，其動(dòng)力系統(tǒng)相當(dāng)于一個(gè)壓力系統(tǒng)，而壓力變化又會(huì)改變相變材料的融化溫度，因此研究壓力對(duì)相變過(guò)程的影響很有必要。

Ma 等[13]通過(guò)分析壓力對(duì)相變材料和液壓油中氣體溶解度的影響來(lái)研究不同壓力條件下動(dòng)力系統(tǒng)的存儲(chǔ)效率；Xia 等[14]通過(guò)在能量方程加入壓力做功的源項(xiàng)，建立了一維相變換熱模型，發(fā)現(xiàn)由于壓力做功相對(duì)于相變潛熱非常小，壓力對(duì)相變過(guò)程的影響也非常小。

材料發(fā)生固液相變時(shí)在融化溫度處是一個(gè)兩相平衡的狀態(tài)，依據(jù)克拉配龍方程[15]，物質(zhì)在固液兩相共存時(shí)，其熔化溫度和壓力存在如下關(guān)系

式中，ρl和ρs分別為材料液相和固相的密度，進(jìn)一步積分可得

式中，下角標(biāo)1、2 分別代表相變材料的不同溫度壓力狀態(tài)。

依據(jù)上式即可計(jì)算得到材料在不同壓力條件下的融化溫度，如圖3 所示。

圖3 融化溫度隨壓力變化曲線(xiàn)Fig.3 Curve of melting temperature versus pressure

2 仿真結(jié)果及分析

2.1 融化過(guò)程

計(jì)算中設(shè)置管內(nèi)的相變材料初始溫度為16℃，此時(shí)低于正十六烷的融化溫度，正十六烷處于凝固狀態(tài)。管外設(shè)置為恒壁溫邊界條件，海水溫度為25℃，熱管從周?chē)Ｋ形諢崃?，溫度升高，?dāng)溫度高于正十六烷融化溫度時(shí)，正十六烷開(kāi)始融化過(guò)程。

從圖4 和圖5 可以看出，隨著時(shí)間的推進(jìn)，熱管外層相變材料最先開(kāi)始融化，并緩慢地向熱管中心推進(jìn)。隨著壓力的上升，熱管內(nèi)相變材料融化速度明顯變慢。在4 000 s 時(shí)5 Mpa 壓力條件下，已有大約85%的相變材料完成融化；而相同時(shí)間30 Mpa 壓力條件下，只有不到30%相變材料變成液態(tài)。這是因?yàn)椋S著壓力升高，材料融點(diǎn)升高，相變材料需要吸收更多的熱量來(lái)達(dá)到相變溫度，在換熱條件相同的條件下就需要花費(fèi)更長(zhǎng)的時(shí)間。融化過(guò)程中，相界面開(kāi)始能保持為圓形，后面由于液體部分的增加，自然對(duì)流作用使得相界面形狀發(fā)生了改變。

圖5 4 000 s 時(shí)融化過(guò)程液相率分布Fig.5 Liquid fraction distribution of melting process at 4 000 s

圖6 為融化過(guò)程中不同壓力條件下液相率的變化，開(kāi)始時(shí)刻材料全部為固態(tài)，此時(shí)液相率的值為0。當(dāng)此值變?yōu)? 時(shí)，表明熱管內(nèi)材料已全部變?yōu)橐簯B(tài)。從圖中可以看出，壓力越高，曲線(xiàn)越靠近橫軸，液相率變化梯度也越小。在0.1，5，10，15，25，30 MPa 壓力下，材料完全融化所需時(shí)間分別約為64，72，86，101，179，271 min，其中5，15 MPa 條件下的融化時(shí)間比0.1 MPa 條件下融化時(shí)間分別長(zhǎng)了12%，34%，而35 MPa 條件下，在300 min 時(shí)融化過(guò)程仍未全部完成。

圖6 不同壓力下液相率變化曲線(xiàn)Fig.6 Curves of liquid fraction at different pressure

圖7 為不同壓力下中心點(diǎn)溫度的變化。初始時(shí)刻，中心點(diǎn)溫度上升很快，此時(shí)熱量主要以顯熱的形式存儲(chǔ)在材料中，當(dāng)達(dá)到融化溫度時(shí)，材料溫度不再上升，熱量以相變潛熱的形式存儲(chǔ)，材料開(kāi)始進(jìn)行融化過(guò)程。當(dāng)材料融化為液態(tài)后，由于外界溫度仍然高于此時(shí)液態(tài)材料，溫度繼續(xù)上升。

圖7 中心節(jié)點(diǎn)溫度變化曲線(xiàn)Fig.7 Temperature change curves at the center node

對(duì)比圖6、圖7 及液相率分布圖可知，當(dāng)中心點(diǎn)完成融化時(shí)，熱管內(nèi)的整個(gè)融化過(guò)程并未全部完成，這正好與云圖結(jié)果相吻合。出現(xiàn)這種情況，是因?yàn)楣虘B(tài)正十六烷比液態(tài)正十六烷密度大，在融化過(guò)程中，熱管中心的固態(tài)正十六烷會(huì)下沉，同時(shí)由于液體部分自然對(duì)流的作用，熱管內(nèi)相變材料最后融化的位置是在中心節(jié)點(diǎn)的下方。

2.2 凝固過(guò)程

為研究壓力對(duì)正十六烷凝固過(guò)程的影響，仿真計(jì)算時(shí)設(shè)置熱管內(nèi)相變材料初始溫度為25℃，管外設(shè)置為恒壁溫邊界條件，海水溫度設(shè)置為16℃。圖8、圖9 分別為凝固過(guò)程1 000 s 和4 000 s時(shí)熱管內(nèi)液相率分布云圖。

圖8 1 000 s 凝固過(guò)程液相率分布Fig.8 Liquid fraction distribution of solidification process at 1 000 s

對(duì)比圖8 和圖9 可知，材料最先從外部開(kāi)始凝固并逐漸向內(nèi)部發(fā)展，壓力越高凝固速率越快。同一壓力條件下，周向上凝固速率相同，整個(gè)過(guò)程相界面都保持為圓形。

圖9 4 000 s 凝固過(guò)程液相率分布Fig.9 Liquid fraction distribution of solidification process at 4 000 s

圖10 為凝固過(guò)程中液相率的變化。初始時(shí)刻材料全部處于融化狀態(tài)，此時(shí)液相率的值為1。當(dāng)此值變?yōu)? 時(shí)，表明材料已全部凝固。從圖中可知，隨著壓力升高，曲線(xiàn)越靠近縱軸，下降梯度也越來(lái)越大，材料凝固速率加快。在0.1，5，10，15，25，30，35 MPa 壓力下，正十六烷完全凝固所需時(shí)間為333，252，203，168，124，113，101 min，其中5，15 MPa 壓力下凝固速率分別比0.1 MPa條件下凝固速率提高了24%和50%。

圖10 凝固過(guò)程液相率變化曲線(xiàn)Fig.10 Change curves of liquid fraction during solidification process

圖11 為凝固過(guò)程中心節(jié)點(diǎn)的溫度變化。對(duì)比圖10、圖11 可以看出，在中心節(jié)點(diǎn)凝固完成時(shí)，整個(gè)熱管內(nèi)相變材料也都完成了凝固。

圖11 凝固過(guò)程中心節(jié)點(diǎn)溫度變化曲線(xiàn)Fig.11 Temperature change curves at the center node during solidification process

3 結(jié)束語(yǔ)

計(jì)算分析表明，壓力會(huì)改變材料的融化溫度進(jìn)而影響材料相變過(guò)程。壓力升高使得材料融化溫度上升，在相同溫度條件下完全融化所需時(shí)間會(huì)變長(zhǎng)，但是凝固所需時(shí)間反而變短。如果融點(diǎn)高于了海水暖水層溫度，相變材料將無(wú)法完成相變過(guò)程，熱滑翔機(jī)將不能穩(wěn)定工作。文中每個(gè)算例仿真計(jì)算過(guò)程中壓力都是取一定值進(jìn)行計(jì)算，實(shí)際相變過(guò)程中壓力是受體積變化的影響，并處于一個(gè)連續(xù)變化的狀態(tài)，如何實(shí)現(xiàn)壓力變化和相變過(guò)程耦合仿真是一個(gè)值得研究的方向。下一步的工作重點(diǎn)將是建立壓力和相變過(guò)程的耦合仿真，使模型更接近于實(shí)際的滑翔機(jī)上的相變過(guò)程；搭建滑翔機(jī)相變過(guò)程仿真平臺(tái)，實(shí)測(cè)相變過(guò)程的壓力變化和相變時(shí)間，并與仿真結(jié)果對(duì)比分析。