一種高容積阻力比水下熱滑翔機殼體外形設計

周徐斌 馬 捷

海洋工程國家重點實驗室上海交通大學，上海200240

0 引言

水下熱滑翔機是一種高效環保、浮力驅動、低噪聲的自主式水下運載器，可廣泛應用于海洋科學考察以及軍事領域。但是，由于海洋溫躍層的溫差較小（往往不超過20℃），溫差能的能量品質較低，導致水下熱滑翔機的熱機效率不僅較低，而且功率往往也不大，因此，需要尋找具有最低阻力特性的殼體外形，以使滑翔機擺脫溫差功率的限制，實現低功率運行，高續航周期的服務。同時，熱滑翔機作為一種水下運載器，不僅需要裝載自身航行所需的動力裝置、控制系統和大量傳感器，而且還要裝載為實現特殊功能和滿足工作需要的大量科研儀器設備，因此，所設計的滑翔機殼體外形必須具有較強的裝載能力，否則，將嚴重限制滑翔機的實際工作性能。

文獻［1-2］采用數值計算的方法對水下航行器，如潛艇、水下滑翔機等的低阻特性進行了研究，并且計算了在不同速度下殼體外形的阻力隨攻角的變化曲線以及其壓力云圖等水動力數據。本文旨在此基礎之上，通過分析不同殼體外形周圍繞流場的情況并總結其優缺點，綜合考慮水動力和裝載能力兩方面的性能，提出一種能夠克服滑翔機殼體外形的低阻特性與裝載能力之間矛盾的新型高性能殼體外形。

本文在水下熱滑翔機殼體外形的設計中，首先將提出一個無量綱的容積阻力比系數，該系數是衡量滑翔機殼體外形的低阻性能和裝載性能的綜合指標，該系數越大，說明殼體外形這兩方面的綜合性能越好，反之則越差。其次，本文將以4種常見殼體外形為基本出發點，通過數值計算的方法，分析出它們各自的優缺點。之后，再采用揚長避短、各取所長的設計思想，結合紡錘體外形的低阻力特性和橢球體頭尾外形的高裝載特性的優點，設計出一種具有高容積阻力比的水下熱滑翔機殼體外形，該外形在具有較好的低阻特性的同時又具有較強的裝載能力。

1 水動力數值計算的理論基礎

本文不考慮流體的熱效應，只關注流場的速度分布和阻力大小。本章將基于質量守恒、動量守恒定律，推導計算水動力學的控制方程。

1）質量連續方程

質量連續方程保證了流體流動的連續性，即單位時間內流體任何體積單元中質量的增加，等于同一時間間隔內流入該體積單元的凈質量，由于滑翔機周圍流體的速度較慢，不考慮流體密度的變化，所以

式中，i=1，2，3，分別代表水動力數值計算的三維正交坐標系中的3個維度；u為速度；x為流場空間點坐標。

2）動量連續方程

動量連續方程保證了經典力學中的動量定律，即數值計算的體積單元中流體的動量對時間的變化率等于外界作用在該體積單元上的各種外力之和。由于滑翔機周圍的流體為海水，因此，忽略流體膨脹粘性系數，并將流體粘性系數視為常數，則動量連續方程為：

式中，i，j=1，2，3，分別代表三維正交坐標系中的3個維度；υ為流體的運動粘度系數；p為壓力，Pa；ρ為密度，kg/m3；F為體積單元所受的外力，N；Sm為流體質量源（匯）。

3）湍流k-ε方程

本文采用雷諾平均（RANS）算法對滑翔機殼體周圍的湍流流動進行計算，即將湍流的N-S方程組對時間進行平均，濾掉高階小量，得到雷諾時均方程組。本文所討論的殼體周圍流體的運動雷諾數在105數量級，引入適用于非較高雷諾數情況下的湍流動能和湍流耗散率的標準k-ε方程使各方程封閉，如式（3）和式（4）所示：

式中，k為湍流動能，m2/s2，；ε為湍流耗散率，，；Gs和GT均為附加湍流生成項，其中gi為重力加速度在i方向上的分量，T為湍流質點的溫度，℃；μt為湍流粘性系數，，其中Ct為一個常數。相關的其他常數值按1972年Launder-Spalding的建議，取值如表 1所示［3］。

表1 k-ε方程相關參數Tab.1 Param etersof k-εfunction

4）壁面函數

如果將壁面作為光滑壁面處理，則流體不會因為摩擦受到來自壁面的阻力作用，這與實際情況不符，同時也會給流場的計算帶來誤差，尤其是在雷諾數不高而又使用k-ε模型的情況下，甚至還會造成較大誤差。為此，本文根據文獻［4］和文獻［5］的研究，引入壁面函數，即在壁面附近流體的速度分布公式，如式（5）所示：

式中，k為卡爾曼數，k=0.4；A為常數，A=5.5；y為流體質點距壁面的距離，m；u*為流體質點相對壁面的摩擦速度，m/s，，其中 τ0為壁面剪應力，N。

式中，當＞1 000時，該范圍內壁面流體的作用力 τ=ρu*2，N。

2 對4種常見殼體外形的綜合性能分析

2.1 4種常見殼體外形的幾何模型

本文研究的水下熱滑翔機殼體外形的總體主尺寸總長1.5m，最大直徑0.2m。一方面，由于這2個參數與水動力性能計算的殼體外形的周圍繞流流場的特征有關［6］，另一方面，為便于分析比較，本文在水下熱滑翔機殼體外形的研究中將這2個參數的值保持不變，4種常見殼體外形的幾何模型如圖1所示，從球體首尾外形殼體到紡錘體殼體，外形逐漸變得細長尖削。

2.2 4種常見殼體外形的低阻特性分析

2.2.1 計算邊界與網格

數值計算的流體域是一個長4.5m、寬0.5m、高0.5m的長方體域，流場的來流邊界距滑翔機殼體的頭部1.4m，去流邊界距尾部1.6m。相對滑翔機機身的尺寸（長1.5m，橫截面的最大直徑0.2 m）而言，流體域的尺寸足夠大，確保了在計算滑翔機水動力性能時流場的進流和尾流都能充分發展，從而保證了數值計算的可靠性［7］。來流邊界設置為固定流速0.5m/s的邊界條件，去流邊界為可自由流出的全壓為零的邊界條件，如圖2所示（以橢球體首尾外形為例）。

圖1 4種常見的殼體外形設計方案Fig.1 Gallery of four common shellshape designs

圖2 數值計算流場域示意圖Fig.2 The flowing area for numerical calculation

結構化網格具有生成速度快、網格質量好、收斂性較高等優點，因此，本文對流體域采用正交六面體結構化網格來進行劃分，如圖3所示。以橢球體首尾外形殼體的計算模型為例，共計3 531 000個網格。圖3（a）為網格劃分完畢的橢球體首尾外形殼體的流場域正視圖，圖3（b）為側視圖。

圖3 數值計算流場域網格劃分示意圖Fig.3 Mesh of the flowing area applied in the numerical calculation

流體域的進流遠端和去流遠端采用粗網格，如圖3所示的淺色線邊長（10～16mm）部分。由于這部分流場與殼體的相互作用以及流動變化已經較小，因而粗網格既可以保證足夠的精度，又可以減少網格總數，達到減輕計算負荷的效果。在壁面附近及流動發生剝離和渦旋形成的地方，則采用精細網格，如圖3所示的深色網格線邊長（3.5～7mm）部分。一方面，要保證滑翔機殼體外形的復雜型線能夠被準確刻畫，另一方面，也能保證計算精度，以準確模擬這些區域的流體流動狀態［8］。

同時，為保證網格質量，將不同邊長的網格與網格之間的大小幅度變化率控制在1～1.1之間，粗細網格之間采用等比漸變的模式調整網格幅度，整個流體域內的網格長寬比（AspectRatio）也控制在1～6以內，最大長寬比為5.3∶1，沒有極端畸變的網格。

采用有限體積法離散控制方程。由于采用高階的六面體網格，對動量方程、湍流方程等采用二階迎風格式離散，壓力和速度的耦合迭代則采用SIMPLEC 算法［9］。

2.2.2 4種滑翔機殼體外形的阻力特性

在相同的收斂條件下（設質量平衡方程和動量守恒方程的最小相對收斂誤差為10-7），采用STREAM軟件和FLUENT軟件求解器對第2.2.1小節中劃分好網格和設置好邊界條件的流場模型進行求解，得到水動力學數據，如表2所示。通過比對用2種軟件分別計算第2.2.1小節所建立的流場域網格模型的數據結果，來驗證兩者的結果是否具有較高的一致性。

對采用2種求解器分別進行數值計算所得到的各外形的阻力數據進行比較可以發現：在壓差阻力計算中，對橢球體首尾外形的求解偏差最大，為3.16%，對圓錐體首尾外形的求解偏差最小，為0.54%；在摩擦阻力計算中，對紡錘體外形的求解偏差最大，為4.23%，對圓錐體首尾外形的求解偏差最小，為0.26%。總而言之，分別采用FLUENT和STREAM進行數值計算得到的阻力數據的偏差均小于5%，結果表現出了較高的一致性，說明數值計算的結果是可靠的。

此外，根據表2中的數據，可知殼體模型的壓差阻力起主導作用，相對而言，摩擦阻力的影響則較小。因摩擦阻力與雷諾數有關，在雷諾數保持不變的基礎上，本文將重點分析壓差阻力。

表2 STREAM和FLUENT求解各種殼體外形阻力系數計算結果比對表Tab.2 Com parison of drag coefficients calcu lation resu lts for d ifferen t shell shapesby STREAM and FLUENT

1）阻力系數

利用數值計算得到的結果，通過下式，可計算各滑翔機殼體外形的阻力系數C［10］：

drag

式中，Fp為壓差阻力，N；Ff為摩擦阻力，N；為流體的動壓，Pa；Sw為濕表面積，m2。

阻力系數是衡量滑翔機殼體外形阻力性能最重要、最直接的一個無量綱指標，因為它與水下航行物體在航行過程中所受的阻力成正比［11］。

根據數值模擬計算結果，利用式（7）計算各滑翔機殼體外形的阻力系數，整理后如表3所示。阻力系數由低到高排列分別為紡錘體外形、橢球體首尾外形、圓錐體首尾外形、球形首尾外形，可見外形越是尖削細長，阻力系數便越低。

2）速度矢量圖分析

表3 4種常見殼體外形阻力系數表Tab.3 Drag coefficientof four comm on shell shapes

流場的速度矢量圖表示穩態情況下流體域內各個節點的平均速度矢量，表征的是數值計算得到的流場各個體積單元流體質點的流速方向和大小。充分了解這些速度信息后就可以了解流場的整體情況，從而分析各種殼體外形的水動力性能發生差異的原因。根據阻力系數計算結果可知，低阻特性最好的是紡錘體外形，最差的是球形頭尾外形。因此，利用速度矢量圖對這兩個極端情況進行分析，尋找滑翔機的阻力成因，可以為本文的設計提供思路。

本文主要觀察流場域中壁面附近的流動，這部分流場與殼體之間的相互作用最為強烈，而在流場域中遠離壁面的地方，由于遠離邊界層，且雷諾數不大，區域中顯示為層流。

球形頭尾外形的速度矢量圖如圖4所示。在球形頭部，因來流受到阻滯，出現了錐形降速區（淺色區域）。由伯努利定理可知，這部分流體若失去速度，則該區域的壓力會上升，從而產生壓差阻力。

圖4 球形首尾外形殼體周圍流場速度矢量圖Fig.4 Vector graph of the overall flow field around the shellwith ballhead and tail

由于尾部為球形，曲率變化急促，尾流還來不及過渡，使得原本緊貼壁面流動的邊界層由于突然脫離壁面，一部分尾流會因動能耗盡且受到其它較快流體的擠壓作用而成為倒流來填充相應的緊靠尾部的空間，因此，在尾部的三角形區域內可能會發生嚴重的邊界層分離現象，這可以通過觀察尾部速度矢量圖來驗證，如圖5所示。由圖可發現，在球形外殼的尾跡中出現了2個渦旋中心，這將導致較大誘導阻力的產生［12］。

紡錘體外形的速度矢量圖如圖6所示。由圖可看出，流動無論是在頭部還是在尾部始終都是平緩過渡。其中首部低速區的范圍最為狹窄，降速幅度也不大，在尾跡中，既沒有渦旋中心，也沒有發生邊界層分離現象，幾乎沒有較大的誘導阻力的產生，因此，其總體阻力最小。

圖6 紡錘體外形殼體周圍流場速度矢量圖Fig.6 Vectorgraph of the flow field around the spindle shell

根據以上分析比較發現，尾部細長，且外形呈流線型，像紡錘體這樣的二次曲率連續的型線的殼體在阻力特性方面具有較大優勢，繞流平緩過渡；而型線較鈍、尾部曲率變化較為急促的外殼則容易產生較大的壓差阻力，從而使得低阻特性變差。

但如果計算這4種殼體外形的有效容積（可用于滑翔機裝載的這部分容積結果整理如表4所示）又可發現：紡錘體外形的有效容積最低，這是因其形狀過于細長；而球形首尾外形的低阻特性雖然最差，但卻具有最大的有效容積，較尖削的圓錐體首尾外形和橢球體首尾外形的有效容積則居中。

表4 4種常見殼體外形有效容積Tab.4 Effective capacity of the four comm on shell shapes

綜上所述，可總結出以下矛盾：呈流線型且尾部細長的殼體外形雖然具有較佳的低阻特性，但有效容積卻較小。

2.3 容積阻力比系數

為解決第2.2小節中出現的低阻特性與有效容積之間的矛盾，更客觀、合理地衡量滑翔機外殼的綜合性能，根據實際設計過程中的需要，本文引入了一個無量綱的“容積阻力比”系數作為衡量水下熱滑翔機殼體外形的可行性和綜合性能指標。由于該系數的分子代表裝載能力的高低，分母表征的是低阻性能，所以該系數的值越大，就說明該殼體既具有較大的裝載量，同時又具有較低的阻力。

式中，Ve為所設計殼體的有效容積，m3；Vs為標準外殼容積［13］，m3，其物理意義是在確定主尺度后外殼可能擁有的最大容積，可由下式計算：

式中，Ls為殼體的最大縱向長度，m；Bs為殼體的最大橫向寬度，m；Hs為殼體的最大垂向高度，m。

3 新型滑翔機殼體外形

根據式（8）計算得到的各殼體外形的容積阻力比系數如表5所示。由表可見，橢球體首尾外形的容積阻力比系數最大，而紡錘體外形的容積阻力比系數則最小，這說明就綜合性能而言，紡錘體外殼是最差的一種殼體設計方案。本文希望通過合理的設計，使外殼外形在阻力特性和裝載能力這兩方面都具有較好的表現。

表5 各殼體外形容積阻力比系數Tab.5 Capacity to d rag ratio of different shell shapes

在殼體外形設計中，本文充分發揮紡錘體外形和橢球體首尾外形各自的優勢，采取將兩者的優勢結合起來的辦法設計了一種新型的水下熱滑翔機殼體外形（圖7）。該殼體外形由3部分組成，其首部為一個長軸為0.2m，短軸為0.1m的半橢球體，尾部為紡錘體從拐點（即最大半徑0.1m處）至尾端點的部分，中間段為將新的首部和尾部連接起來的圓柱體。這種殼體外形克服了紡錘體總體過于細長、有效容積小的缺點，但又最大程度地保留了紡錘體外形中起關鍵作用的低阻部分。

通過SOLIDWORKS軟件計算，新型殼體外形的有效容積為0.038m3。

圖7 新型水下熱滑翔機殼體外形正視圖Fig.7 Frontview ofnew designed shell shape of the ocean thermalglider

在阻力特性方面，對新型殼體外形利用STREAM軟件進行了求解，并且使其邊界條件和收斂條件與第2節的水動力計算的數值保持一致，計算得到的新型水下熱滑翔機殼體外形的摩擦阻力、壓差阻力、總阻力和濕表面積的數據如表6所示。以表6中的數據為基礎，根據式（7）計算新型殼體外形的阻力系數，可得其結果為0.006。

表6 新型水下熱滑翔機殼體外形水動力數據Tab.6 H yd rodynam ic data of the new designed shape shellof ocean therm alglider

最后，根據表6的水動力數據和式（8），計算出了新型水下熱滑翔機殼體外形的容積阻力比系數，為1.200 1，其要比橢球體首尾外形高6.5%，比圓錐首尾外形高17.06%，比球體首尾外形高13.99%，比紡錘體外形高22.6%。

4 結 論

本文的出發點是利用CFD數值計算的方法來分析各種殼體外形的水動力性能，首先分別利用STREAM和FLUENT軟件對具有相同邊界條件和收斂條件的數值計算模型進行了求解，兩者數值解的最大相對誤差不超過5%，表現出了較高的一致性，說明數值解是可靠的。通過對數值計算得到的速度矢量圖和阻力系數等進行分析比較，發現外形呈流線型且尾部尖而細長的殼體的低阻特性雖然較好，但有效容積小、裝載能力較差。

在對水下熱滑翔機殼體外形進行綜合性能分析時，不僅要考慮殼體外形的低阻特性，還要考慮殼體外形的裝載能力，因此，本文提出了一種采用容積阻力比系數來進行綜合比較的方法。同時，還設計了一種半橢球首、紡錘體尾的新型殼體外形。與4種常見的滑翔機殼體外形相比，新型殼體外形具有最大的容積阻力比系數，是一種綜合性能較高、低阻特性和裝載能力較為平衡的外形設計方案，可作為水下熱滑翔機殼體外形設計的一種新選擇。

由于水下熱滑翔機也屬水下潛器的一種，在軍事應用中就相當于一艘微型潛艇，因而本文有關水下熱滑翔機殼體外形的水動力分析的研究對水下潛艇的設計而言也具有一定的借鑒意義。首先，本文采用正交六面體結構化網格對全流場域劃分的方法也可以應用到潛艇殼體的水動力分析之中，它具有劃分速度快、網格質量高、計算效率高等特點。其次，本文提出的容積阻力比系數既考慮了容積裝載效率，又考慮了低阻特性，因潛艇同樣也要考慮這兩方面的綜合性能，所以該系數也可以為某些潛艇的殼體外形設計提供參考。

［1］盧錦國，梁中剛，吳方良，等.水下航行體回轉水動力數值計算研究［J］.中國艦船研究，2011，6（6）：8-12，27.

LU JG，LIANG ZG，WU F L，etal.Numerical calculation on hydrodynamic performance of the submerged vehicle in turning motion［J］.Chinese Journal of Ship Research，2011，6（6）：8-12，27.

［2］馬冬梅，馬錚，張華，等.水下滑翔機水動力性能分析及滑翔姿態優化研究［J］.水動力學研究與進展，2007，22（6）：703-708.

MA D M，MA Z，ZHANG H，et al.Hydrodynamic analysis and optimization on the gliding attitude of the underwater glider［J］.Journal of Hydrodynamics，2007，22（6）：703-708.

［3］柏鐵朝，梁中剛，周軼美，等.潛艇操縱性水動力數值計算中湍流模式的比較與運用［J］.中國艦船研究，2010，5（2）：22-28.

BAO T C，LIANG ZG，ZHOU Y M，et al.Comparison and application of turbulencemodes in submarinemaneuvering hydrodynamic forces computation［J］.Chinese Journalof Ship Research，2010，5（2）：22-28.

［4］須賀一彥，CRAFT T J，LACOVIDESH.汎用的な解析的壁関數モデル（第1報、滑面～粗面亂流に対応した流れ場のモデル）［C］//日本機械學會論文集B編.東京，2005.

KAZUHIKO S，CRAFT T J，LACOVIDESH，et al.A generalized analytical wall-function for turbulence（1st Report，a flow field model for smooth and rough wall turbulence）［C］//The Japanese Society of Mechanical Engineers.Tokyo，2005.

［5］CRAFT T J，GERASIMOV A V，LACOVIDESH，et al.Progress in the generalization of wall-function treatments［J］.International Journal of Heat Fluid Flow，2002，23：148-160.

［6］吳利紅，俞建成，封錫盛.水下滑翔機器人水動力研究與運動分析［J］.船舶工程，2006，28（1）：12-16.

WU L H，YU JC，FENG X S.Hydrodynamic research and motion analysis of AUG［J］.Ship Engineering，2006，28（1）：12-16.

［7］熊耀宇，段宏，姜治芳，等.數值模擬現代水面艦船帶自由面的湍流繞流場［J］.中國艦船研究，2010，5（6）：16-20，32.

XIONG Y Y，DUAN H，JIANG Z F，et al.Numerical simulation ofmodern surface ship in free surface turbulence flow［J］.Chinese Journal of Ship Research，2010，5（6）：16-20，32.

［8］REN H L，LIUW X.Calculation ofhydrodynam ic coefficients of floating body with complex wetted surface［J］.Journalof Ship Mechanics，2008，12（3）：359-367.

［9］操盛文，吳方良.尺度效應對全附體潛艇阻力數值計算結果的影響［J］.中國艦船研究，2009，4（1）：33-37，42.

CAO SW，WU F L.Investigation of scaling effects on numerical computation of submarine resistance［J］.Chinese Journal of Ship Research，2009，4（1）：33-37，42.

［10］王沖，劉巨斌，張志宏，等.水下滑翔機沿縱剖面滑行時水動力特性計算與分析［J］.艦船科學技術，2009，31（1）：134-136，163.

WANG C，LIU JB，ZHANG Z H，et al.Numerical research on dynamic characteristic of underwater glider when it runs in longitudinal section［J］.Ship Science and Technology，2009，31（1）：134-136，163.

［11］吳方良，吳曉光，馬運義，等.潛艇實艇阻力預報方法研究［J］.中國艦船研究，2009，4（3），28-32.

WU F L，WU XG，MA Y Y，etal.Method of predicting the total submerged resistance of submarines［J］.Chinese Journal of Ship Research，2009，4（3）：28-32.

［12］羅寅，張阿漫，朱永凱，等.結合亞格子模型的Lattice-Boltzmann算法［J］.中國艦船研究，2010，5（1）：10-13.

LUO Y，ZHANG A M，ZHU Y K，et al.Combination of the Lattice-Boltzmann method and the sub-grid model［J］.Chinese Journal of Ship Research，2010，5（1）：10-13.

［13］WEBB D C，SIMONETTIP J，JONESC P.SLOCUM：an underwater glider propelled by environmental energy［J］.IEEE Journal of Oceanic Engineering，2001，26（4）：447-452.