CFD動態網格技術在水下航行體非定常操縱運動預報中的應用研究

李迎華，吳寶山，張 華

（中國船舶科學研究中心，江蘇 無錫214082）

1 引 言

對于水下滑翔機[1-2]這一類以剩余浮力為主要驅動方式的水下低速航行體而言，其在進行縱傾與浮力調節時，所作的運動是非定常的，因此對其非定常操縱運動的預報研究，是指導其設計開發的基礎。

由于水下低速航行體非定常操縱運動中其雷諾數低于臨界雷諾數，水動力系數會隨時間發生變化，因而采用傳統的拘束模型試驗[3-5]，以及基于準定常假設的數值方法[6]對其進行運動預報時，則需要對各種雷諾數下的水動力分別進行測試和計算，這顯然需要大量的工作，而且對于拘束模型試驗來說，其代價也是極其昂貴的。近年來，隨著商業軟件Fluent的不斷發展，基于RANS求解器的CFD數值預報方法已成為一種經濟、有效的手段。因此，研究非定常操縱運動的CFD數值模擬技術對水下低速航行體設計研究具有重要的意義。

對于數值模擬時水動力及運動參數需要實時求解的這類操縱運動，Fluent專門提供了一種研究手段—動網格技術。在模擬時將運動方程通過用戶自定義函數（UDF）導入至Fluent中，即可實現水動力和運動參數實時自主求解。因此，動網格技術適合于水下低速航行體非定常操縱運動的直接模擬。

目前動網格技術大多應用于周期性操縱運動的數值計算，如垂蕩、升沉等，所采用的動網格劃分形式一般是全非結構網格[7]，由于這種網格劃分形式在近壁面粘性邊界層內難以生成大展弦比的非結構網格，且網格數量較結構網格多出許多，因而將其推廣到非周期性操縱運動預報上有一定的局限性，因此，近年來，結合結構網格和非結構網格優勢的混合網格技術得到了重視和發展，張來平等[8]提出了三棱柱/四面體/矩形三維動態混合網格劃分形式，并用無粘流的圓球超聲速繞流算例進行了驗證，得到了與定常計算相一致的結果。

本文在此基礎上對三維動態混合網格模型進行了發展，提出了六面體/四面體三維動態混合網格模型，并以此建立了基于動網格技術的CFD數值預報方法。本文首先以帶攻角均勻直航的滑行狀態為研究對象，應用動網格技術對一典型水下低速航行體縮比模型主體及帶附體的水動力進行了計算以驗證其水動力的計算精度，為便于比較，同時給出了靜態網格計算結果，進而將其應用于一水下滑翔機自航模非定常滑翔運動的預報，以驗證動網格技術的應用能力。

2 基于動網格技術的數值計算方法

2.1 控制方程

時均的不可壓縮連續性方程和N-S方程（RANS方程）如（1）、（2）式所示：

2.2 湍流模型

對于本文研究的水下低速航行體而言，其運動過程中縱向雷諾數較低。RNG k-ε湍流模型提供了針對低雷諾數的有效粘性的微分解析式，且具備數值穩定性好和處理流線彎曲程度較大的流動更優等優點[9]，本文數值計算即采用RNG k-ε湍流模型。

2.3 幾何守恒律

動網格需滿足幾何守恒律，控制體積的時間導數由（3）式計算：

式中：nf是控制面的數目，是j面的面積矢量。

式中：δVj是控制面j在一個時間步長內的體積更新量。

2.4 邊界條件

（1）遠場邊界條件：采用零剪應力壁面邊界條件處理。

（2）物面條件：滿足壁面黏附條件，壁面處流體速度與運動邊界速度相同。

2.5 動態混合網格生成方法

本文的六面體/四面體動態混合網格模型如下建立：在近壁面處依次生成貼體六面體網格和四面體網格共同構成運動區域隨運動邊界一起運動 （四面體網格與六面體網格之間采用交界面過渡），外場生成四面體網格并時刻保持靜止，在運動區域和靜止區域之間采用四面體網格過渡構成變形區域。當物體位移較小時，采用“類彈簧原理”對變形區域網格節點進行松弛，當物體位移較大導致變形區域網格質量較差時，則在局部重新生成網格。

本文采用上述改進的動態混合網格模型對一典型水下低速航行體縮比模型表面及其周圍流動區域進行了網格劃分，此模型主體為細長回轉體，長2 150mm，直徑300mm，在主體中后部布置有一對后掠的水平翼，翼兩端間距1 200mm。網格劃分的具體形式為：主體C-O型（縱向-周向），附體C-O型（展向—弦向）；計算域大小為：沿航行體首前端、航行體徑向分別延伸3L距離，沿航行體尾后端延伸4L距離（L為航行體縱向長度）。單獨主體網格數65萬，帶附體網格數120萬。圖1、2為帶附體表面及近壁面流場區域網格劃分示意圖。為便于對比分析，也給出了靜態網格劃分示意圖，如圖3所示。單獨主體的靜態網格數量32萬，帶附體的靜態網格數量57萬。

2.6 動網格更新算法

動網格更新算法有彈簧近似光滑算法（Spring Based Smoothing）、動態分層算法（Dynamic Layering）和局部重構算法（Local Remeshing）。對于本文研究的四面體/六面體混合網格，選用彈簧近似光滑算法和局部重構算法作為動網格更新算法。

2.7 數值離散和求解

（1）時間項的離散：采用直接一階隱式離散。

（2）空間項的離散：其中擴散項以中心差分格式差分，對流項采用二階迎風格式。

應用SIMPLE法處理壓力速度耦合問題，離散方程以Gauss-Seidel迭代法求解。

2.8 時間步長的選取

由于時間項采用一階隱式離散，是無條件穩定的，因此時間步長的選取可以不用考慮離散格式的穩定性，理論上可以取較大的時間步長，但過大的時間步長將導致網格更新過快，進而出現負體積。本文從網格更新角度推導出時間步長計算公式如下：

不考慮控制面上的面積更新量則由（4）式可得對于單元j上的當地時間步長：

式中：αj=hj/Δhj，Δhj為單元j上的高度更新量，hj為單元j上更新前的高度。

全場時間步長應該取全動區域最小的即

式中：h1表示動區域最小的第一層網格高度；vel為邊界運動的合速度；α為一系數。

2.9 雙振蕩圓板驗證算例

為對動網格技術的可行性進行驗證，本文選取文獻[10]中的雙強迫振蕩圓板算例進行了計算，并與文獻給出的結果進行了對比。采用全非結構網格劃分形式對其進行網格劃分，網格數量12萬，如圖4所示。

采用彈簧近似光滑算法和局部重構算法作為動網格更新算法，將雙板振蕩速度以用戶自定義函數（UDF）形式導入到Fluent中，時間步長 Δt=0.1s。

本文計算得到的阻力系數變化周期T=25s，阻力系數幅值A=13.7；而文獻給出的阻力系數變化周期T=25s，阻力系數幅值A=14.2，可見本文采用的動網格技術可以得到相對比較精確的計算結果，因此可以應用于水下滑翔機非定常運動的預報。

3 水下滑翔機非定常運動的CFD計算研究

本文以中國船舶科學研究中心自主開發的“前哨”號水下滑翔機自航模[2]為對象，應用動網格技術對其非定常操縱運動進行了預報，以進一步驗證動網格技術的應用。

操縱運動預報的基礎在于水動力的預報，因此本文首先以帶攻角均勻直航的滑行狀態為研究對象，應用動網格技術對前述典型水下低速航行體縮比模型的主體及帶附體水動力的預報精度進行了驗證，為便于比較分析，同時給出了靜態網格計算結果。為分析水動力計算精度對運動預報的影響，本文還基于帶附體動網格計算的水動力對定常滑翔運動參數進行了分析。

3.1 基于動網格技術的水動力預報驗證

通過用戶自定義函數 （UDF）將水下低速航行體運動速度實時導入到Fluent中，時間步長Δt=0.002s，應用基于動態混合網格模型的動網格技術依次計算了雷諾數Re=2.4×106，攻角分別為6°、8°、10°、12°四種工況下主體及帶附體模型的水動力，得到了無因次水動力系數和水動力中心作用點位置lα′（lα′=-M′/Z′），并與風洞模型試驗結果[12]進行了對比。

（1）主體模型水動力計算結果

主體模型靜態網格和動網格計算結果如表1所示，水動力隨攻角變化曲線如圖5-8所示。

表1 水下航行體主體模型水動力計算結果（Re=2.4×106）Tab.1 The calculated results of the hydrodynamic forces for the main body of the underwater vehicle model（Re=2.4×106）

由計算結果可知，除個別工況外，主體模型靜態網格計算的軸向力系數相差16%以內，垂向力系數相差8%以內，俯仰力矩系數相差5%以內，水動力作用中心位置相差5%以內；動網格計算的軸向力系數相差5%以內，垂向力系數、俯仰力矩系數均相差10%以內，水動力作用中心位置相差13%以內。

（2）帶附體模型水動力計算結果及對運動預報的影響

帶附體模型操縱性水動力計算結果如表2所示，水動力計算結果隨攻角變化曲線如圖9-12所示（由于靜態網格在12°時失速，因此只給出前三個攻角對應的水動力值）。

由計算結果可知，除個別工況外，帶附體模型靜態網格計算的軸向力系數相差16%以內，垂向力系數和俯仰力矩系數相差14%以內，水動力作用中心位置相差5%以內；動網格計算的軸向力系數相差10%以內，垂向力系數相差9%以內，俯仰力矩系數相差13%以內，水動力作用中心位置lα′相差10%以內。

表2 水下航行體帶附體模型水動力計算結果（Re=2.4×106）Tab.2 The calculated results of the hydrodynamic forces for the main body with appendages of the underwater vehicle model（Re=2.4×106）

由主體及帶附體模型水動力計算結果可知，基于動態混合網格模型的動網格技術水動力預報結果與試驗吻合較好，且可以達到與靜態網格相類似的精度。

同時從計算結果也可看出：數值計算結果與試驗仍存在一定差距，這在很大程度上是因為實際流動處于層流和湍流的混合區，而本文采用適用于充分發展湍流（Re>1×107）的湍流模型去解算，因而會帶來一定的誤差。

為分析水動力計算精度對運動預報的影響，本文基于上述動網格計算的水動力和試驗測得的水動力，分別計算了給定雷諾數Re=2.4×106下帶附體模型穩定滑翔時不同剩余浮力系數k（k=ΔB/B，ΔB為剩余浮力，B為水下靜均衡后的浮力） 對應的潛浮角 χ （χ=θ-α，θ為縱傾角，α 為攻角）、水平速度分量uξ和垂直速度分量uζ，計算結果曲線如圖13-15所示。

由計算結果來看：動網格預報的結果與試驗預報結果相比，潛浮角相差0.4°以內，水平速度和垂直速度相差0.007m/s以內，表明了上述動網格水動力計算誤差對運動預報的影響比較小。

3.2 基于動網格技術的非定常運動的CFD預報

對“前哨”號水下滑翔機自航模垂直面內滑翔運動的預報可采用如圖16所示的兩類坐標系：地面坐標系E-ξζ和運動坐標系o-xz，其中原點o位于水下航行體浮心處，ox軸指向首部為正，oz軸垂直ox軸指向下方為正。圖中X，Z，M分別為此水下航行體所受的軸向、垂向水動力和縱傾力矩。

此水下滑翔機垂直面內非定常滑翔運動方程可簡化為如（7）～（12）式所示的形式：

式中：W（t）為水下滑翔機重量；xg（t）為水下滑翔機重心縱向坐標；zg（t）為水下滑翔機重心垂向坐標；B為水下滑翔機所受浮力。

將滑翔運動方程（7）～（12）進行時間一階向前差分離散，并把離散后的方程以UDF形式導入到Fluent中，即可實現Fluent中每一時間步的水動力和運動參數實時自主求解。其求解流程如圖17所示。

本文采用基于動態混合網格模型的動網格技術對“前哨”號水下滑翔機非定常下潛運動中的一段進行了數值模擬，結果如圖18、19所示。

從數值計算結果與試驗的比較來看：基于動網格技術數值方法的預報結果與試驗符合較好，且變化趨勢一致。

4 結論和展望

本文基于傳統的動態混合網格模型，發展了一種改進的三維動態混合網格模型，以此建立了非定常操縱運動數值預報方法，并以帶攻角均勻直航的滑行狀態為研究對象，應用動網格技術對一典型水下低速航行體縮比模型主體及帶附體的水動力進行了計算，結果與風洞模型試驗結果吻合較好，從而驗證了基于動態混合網格模型的動網格技術對操縱運動水動力的預報精度；本文進而將其應用于一水下滑翔機自航模非定常運動的預報，預報結果與試驗符合較好，初步驗證了動網格技術的應用能力和本文建立的方法的有效性。同時動網格技術具有更廣闊的應用前景，主要表現在：

（1）動網格技術可以應用于水下多體相對運動，潛艇的坐底、靠岸等的數值模擬。

（2）動網格技術可以應用于非定常操縱運動過程中流體動力的精細研究。

[1]馬冬梅,馬 錚.水下滑翔機運動及水動力特性研究[R].無錫:中國船舶科學研究中心科技報告,2007.

[2]李 龍,張 華,陳季軍,馬 錚,吳寶山.水下滑翔機自航模演示試驗研究[C].中國造船工程學會學術論文集,2008,No.200804.

[3]Ma Ling,Cui Weicheng.Simulation of dive motion of a deep manned submersible[J].Journal of Ship Mechanics,2004,8(3):31-38.

[4]張 華,吳寶山.扁平型潛水器空間大機動數學模型研究[C].中國造船工程學會學術論文集,2004,No.200405.

[5]謝俊元,馬 嶺,胡 震.載人深潛器運動操縱模擬[J].中國造船,2006,47(2):62-69.

[6]王長濤,俞建成,吳利紅,封錫盛.水下滑翔機器人運動機理仿真與實驗[J].海洋工程,2007,25(1):64-69.

[7]黃昆侖,龐永杰,蘇玉民,朱 軍.潛器線性水動力系數計算方法研究[J].船舶力學,2008,12(5):697-703.

[8]張來平,張涵信.三維動態混合網格技術及非定常計算方法[C].全國流體力學青年研討會論文集,2005.

[9]Fluent 6.3.26 User’s Guide.Fluent Inc.[CP].

[10]顧 正.二維單圓柱、雙圓柱繞流問題及三維振蕩板運動的數值模擬[D].上海:上海交通大學,2007.

[11]施生達.潛艇操縱性[M].北京:國防工業出版社,1995.

[12]倪歆韻,潘子英.低速水下航行體風洞模型試驗[R].無錫:中國船舶科學研究中心科技報告,2008.