基于FLUENT 的水下熱滑翔機相變過程仿真分析

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

(西北工業大學 航海學院，陜西 西安，710072)

0 引言

我國海洋資源十分豐富，人們對海洋資源的探索開發離不開各種水下航行器[1]。水下熱滑翔機是一種新型無人自主水下航行器，其利用相變材料在固液相變過程中產生的體積變化，將海洋溫差能轉化為機械能，并完成水下的之字形運動。

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

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

1 相變換熱模型

1.1 物理模型

圖1 為熱管截面相變過程示意圖，管內部存放相變材料，整個圓筒處于海水環境中。材料初始狀態為液態，當管外海水低于材料相變溫度時，相變材料開始凝固；當管外海水溫度高于材料相變溫度時，管內相變材料開始融化。熱滑翔機所用相變材料為正十六烷，其物性參數受溫度影響較小。

為簡化計算模型進行如下處理：

圖1 熱管結構示意圖Fig.1 Schematic diagram of heat pipe structure

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

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

3) 液態工質溫度分布不均會使其密度分布不均勻，進而產生自然對流作用，在仿真計算過程中密度選項選擇布辛奈斯克(Boussinesq)近似，其作用是在動量方程中浮力項中，將原本為常數的密度改為與溫度有關的函數。

由于水下熱滑翔機熱管為細長圓管，溫度隨管長變化不明顯，故截取橫截面分析，使用二維結構網格。網格用ICEM CFD 繪制，圓面直徑為24 mm。經網格無關性檢驗，取網格單元數量9 396，節點數量9 289 時，計算能保證較高的精度和較快的計算速度。將網格導入FLUENT 18.0，選擇二維分離隱式雙精度求解器，并將圓截面中心點作為觀測節點，監測該點溫度變化。熱滑翔機在水下以緩慢速度滑翔，其熱管與周圍海水之間為對流換熱，假設滑翔機運行速度為0.28 m/s[10]，經計算取對流換熱系數為324 W/(m2·K)，仿真計算時間步長取0.1 s。

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

1.2 數學模型

FLUENT 軟件中的Solidification/Melting 模型[11]可以用來仿真材料相變過程，模型中能量方程采用焓法[12]求解，在固態、液態和固液糊狀區建立統一的能量方程，不需要跟蹤相界面的移動，以簡化計算過程。

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

式中：H為單位體積相變材料的總焓；ρ為材料密度；υ為流體速度；S為動量源項；k為導熱系數；ν為運動粘度；p為壓強。

焓孔隙率技術將糊狀區域(部分凝固區域)作為多孔介質，每個微元內的孔隙度等于該孔內的液體分數。在動量方程添加動量源項

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

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

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

焓法求解能量方程的過程中，液相率是用溫度定義的，故須建立液相率和工質焓的關系。在計算過程中，工質焓定義為

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

根據以上公式建立焓和溫度的關系。對于相變過程中材料的導熱系數采用如下表達形式

式中，ks、kl分別為材料固相和液相的導熱系數。

材料比熱容示意見圖2，具體定義如下

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

1.3 壓力對相變過程的影響

熱滑翔機相變材料在融化時體積膨脹會壓縮蓄能器中的氣體，其動力系統相當于一個壓力系統，而壓力變化又會改變相變材料的融化溫度，因此研究壓力對相變過程的影響很有必要。

Ma 等[13]通過分析壓力對相變材料和液壓油中氣體溶解度的影響來研究不同壓力條件下動力系統的存儲效率；Xia 等[14]通過在能量方程加入壓力做功的源項，建立了一維相變換熱模型，發現由于壓力做功相對于相變潛熱非常小，壓力對相變過程的影響也非常小。

材料發生固液相變時在融化溫度處是一個兩相平衡的狀態，依據克拉配龍方程[15]，物質在固液兩相共存時，其熔化溫度和壓力存在如下關系

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

式中，下角標1、2 分別代表相變材料的不同溫度壓力狀態。

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

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

2 仿真結果及分析

2.1 融化過程

計算中設置管內的相變材料初始溫度為16℃，此時低于正十六烷的融化溫度，正十六烷處于凝固狀態。管外設置為恒壁溫邊界條件，海水溫度為25℃，熱管從周圍海水中吸收熱量，溫度升高，當溫度高于正十六烷融化溫度時，正十六烷開始融化過程。

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

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

圖6 為融化過程中不同壓力條件下液相率的變化，開始時刻材料全部為固態，此時液相率的值為0。當此值變為1 時，表明熱管內材料已全部變為液態。從圖中可以看出，壓力越高，曲線越靠近橫軸，液相率變化梯度也越小。在0.1，5，10，15，25，30 MPa 壓力下，材料完全融化所需時間分別約為64，72，86，101，179，271 min，其中5，15 MPa 條件下的融化時間比0.1 MPa 條件下融化時間分別長了12%，34%，而35 MPa 條件下，在300 min 時融化過程仍未全部完成。

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

圖7 為不同壓力下中心點溫度的變化。初始時刻，中心點溫度上升很快，此時熱量主要以顯熱的形式存儲在材料中，當達到融化溫度時，材料溫度不再上升，熱量以相變潛熱的形式存儲，材料開始進行融化過程。當材料融化為液態后，由于外界溫度仍然高于此時液態材料，溫度繼續上升。

圖7 中心節點溫度變化曲線Fig.7 Temperature change curves at the center node

對比圖6、圖7 及液相率分布圖可知，當中心點完成融化時，熱管內的整個融化過程并未全部完成，這正好與云圖結果相吻合。出現這種情況，是因為固態正十六烷比液態正十六烷密度大，在融化過程中，熱管中心的固態正十六烷會下沉，同時由于液體部分自然對流的作用，熱管內相變材料最后融化的位置是在中心節點的下方。

2.2 凝固過程

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

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

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

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

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

圖10 凝固過程液相率變化曲線Fig.10 Change curves of liquid fraction during solidification process

圖11 為凝固過程中心節點的溫度變化。對比圖10、圖11 可以看出，在中心節點凝固完成時，整個熱管內相變材料也都完成了凝固。

圖11 凝固過程中心節點溫度變化曲線Fig.11 Temperature change curves at the center node during solidification process

3 結束語

計算分析表明，壓力會改變材料的融化溫度進而影響材料相變過程。壓力升高使得材料融化溫度上升，在相同溫度條件下完全融化所需時間會變長，但是凝固所需時間反而變短。如果融點高于了海水暖水層溫度，相變材料將無法完成相變過程，熱滑翔機將不能穩定工作。文中每個算例仿真計算過程中壓力都是取一定值進行計算，實際相變過程中壓力是受體積變化的影響，并處于一個連續變化的狀態，如何實現壓力變化和相變過程耦合仿真是一個值得研究的方向。下一步的工作重點將是建立壓力和相變過程的耦合仿真，使模型更接近于實際的滑翔機上的相變過程；搭建滑翔機相變過程仿真平臺，實測相變過程的壓力變化和相變時間，并與仿真結果對比分析。