球體結構近水面運動過程附加質量研究*

張曉強，王冠，馮金富

(1.空軍工程大學 航空航天工程學院，陜西 西安　710038；2.中國人民解放軍95899部隊，北京　100085)

?

仿真技術

球體結構近水面運動過程附加質量研究*

張曉強1，王冠2，馮金富1

(1.空軍工程大學 航空航天工程學院，陜西 西安710038；2.中國人民解放軍95899部隊，北京100085)

摘要：基于CFD平臺，運用VOF法和動網(wǎng)格技術，針對球體結構近水面運動過程提出了一種新的附加質量求解方法，推導了相應的數(shù)學公式，構建了數(shù)值模型，設定不同的運動規(guī)律對方法進行了驗證分析，在此基礎上，計算了不同浸深狀態(tài)時球體結構的附加質量。結論表明，提出的求解方法正確有效；接近自由面過程中，球體結構的垂向和橫向附加質量均呈下降趨勢，頂端剛好接觸自由面時分別下降29.6%和16.2%，對球體結構穩(wěn)定性、安全性設計的影響不可忽略。

關鍵詞：球體結構；附加質量；自由液面；水動力；流固耦合；動網(wǎng)格；數(shù)值計算

0引言

球體結構廣泛應用于海洋工程、船舶工程等領域，如潛艇、深潛器、海洋營救裝置等[1-2]。目前，對此類結構物的經(jīng)濟性、穩(wěn)定性、安全性的設計要求日益增加，需要設計者研究球體結構的整體受力，明晰其水動力特性。因此，對于球體結構水動力特性尤其是近水面水動力特性的研究具有十分重要的工程意義。文獻[3]研究了有限深水中潛浮式球體的輻射與繞射問題，文獻[4]研究了有限深水中球形潛體產(chǎn)生的速度場及自由面上的位移場，但依然無法解決球體結構近水面運動過程中與水流體的運動耦合問題。

工程上通常采用附加質量方法進行解耦，附加質量是對流體中運動的物體改變運動狀態(tài)時受力特性的一種描述，是對潛體運動進行分析時的重要參數(shù)，對其動力外形布局設計、運動特性預報均具有指導價值?？偟膩碚f，附加質量求解有3類方法：①理論推導，原則上存在解析解，適于球體、柱體等簡單外形的物體[5]。②水池試驗，實驗方法復雜、耗時，特別當物體尺寸較大或線型復雜時，代價過高甚至無法承受，不利于初期研究[6-7]。③數(shù)值方法，隨著CFD(computational fluid dynamics)技術的發(fā)展，數(shù)值方法能夠進行無窮域任意外形、尺寸物體的附加質量計算，速度與精度受限于硬件配置和湍流模型的研究[8-9]。

目前，對于運動物體附加質量的研究，多集中于無窮域，即假設物體是在無窮的流體域中運動的，而且無窮域附加質量的經(jīng)典推導和結論已在實踐中得到廣泛應用[10-13]，而對于物體近水面運動的附加質量研究較少。無窮域附加質量的經(jīng)典結論及數(shù)值方法無法直接推廣到帶自由表面的問題中去[14]，尤其是對于各類水面附近運動物體的特性研究。

本文針對近水面運動球體結構，基于CFD平臺進行二次開發(fā)，提出了一種新的附加質量求解方法，并與理論計算結果進行了對比驗證，結論可應用于球體結構近水面運動的性能分析和預報。

1附加質量求解思路

球體結構近水面運動時主要受推力、重力、浮力、水動力等作用力，如圖1所示。

圖1　球體結構受力示意圖Fig.1　Schematic diagram of spherical structure's force

球體結構的動力學方程可寫為

(1)

把附加質量有關項移到方程右邊，可得

(2)

可以認為，方程右端后3項即為物體運動時所受的水作用力，以W表示即為

(3)

從W的組成部分可以看出，它包括了流體對運動物體的所有作用，考慮到CFD平臺下，物體以某一規(guī)律運動時，對物面邊界的壓強進行矢量積分就是W，提出以下求解策略：

垂向運動、浸深H情況下，以(v,a)運動時的流體作用力為Wa，

Wa=B-Fv-λ′v-λa,

(4)

以(v，-a)運動時的流體作用力為W-a，

W-a=B-Fv-λ′V+λa，

(5)

方程兩邊同時相減，可得附加質量的計算公式，

(6)

λ即為球體在浸深h時的垂向附加質量。同理，橫向運動情況下，不考慮浮力B，亦可得式(6)，由于計算域的對稱性，沿2個橫向坐標軸的運動實質上是相同的，文中均作橫向討論。式中，只有監(jiān)控受力和加速度，排除了其他參數(shù)的影響。由以上過程可知，如果在無窮域，λ為常量，那么只要得出定常運動時與加速運動時的受力之差，也可得出物體的附加質量，這也是文獻[7-8]中作者通過數(shù)值計算物體附加質量的思路，可以看做是本文計算方法的一種特殊情況。

2數(shù)學模型

由前面的推導過程可知，本文的計算策略同樣適用于無窮域，而無窮域球體的附加質量存在理論解析解，也有比較可靠的試驗資料，為驗證方法的正確性，以無窮域球體的附加質量作為計算參考。

本文的數(shù)值模擬包含水和空氣2種流體，所以把流動當作兩相流，自由面為水與空氣的交界面，整個過程是一個非定常過程，為了能夠精確地跟蹤描述水空自由面，采用成熟的VOF(volume of fluid)方法處理，主要基于以下優(yōu)點[15]：①VOF法基于相函數(shù)的概念，追蹤的是網(wǎng)格單元中流體體積的變化，而非追蹤自由液面流體質點的運動，能夠很好地描述自由面的波動、翻卷、噴濺，界面的銳利程度較高。②相對于MAC(marker and cell)方法，VOF法可以處理自由面重入等強非線性現(xiàn)象，所需計算時間更短、存儲量更少，對于硬件的要求尤其是內(nèi)存大小相對較低。

2.1流場控制方程[15-16]

流動速度遠遠小于音速，所以認為水和空氣均為不可壓縮流體，故采用以下三維控制方程：

▽·U=0,

(7)

(8)

式中：U表示流體質點的速度矢量；ρ，μ分別為流體的密度和動力黏性系數(shù)；g為重力加速度；FSV為表面張力的等價體積力形式，借助于相函數(shù)F獲得。

VOF法基本原理是利用計算網(wǎng)格單元中流體體積量的變化和網(wǎng)格單元本身體積的比值函數(shù)F來確定自由面的位置和形狀。跟蹤相與相之間的界面是通過求解一相或多相的容積比率的連續(xù)方程來完成的。下標q=1時代表水，下標q=2時代表空氣，對第q相，這個方程如下：

(9)

方程右端的源項為0，于是求解方程可寫為

(10)

主相容積比率的計算基于如下的約束：

(11)

考慮水的黏性，利用穩(wěn)定性較好的湍流模型“標準k-ε模型”對控制方程進行封閉。其模型方程為

湍流動能方程k:

Gk+Gb-ρε+Sk；

(12)

擴散方程ε:

(13)

式中：Gk表示由層流速度梯度而產(chǎn)生的湍流動能；Gb為由浮力產(chǎn)生的湍流動能；C1ε，C2ε為常量；σk和σε為k方程和ε方程的湍流Prandtl數(shù)；Sk和Sε為用戶定義的源項。

2.2網(wǎng)格劃分

整體計算域為20 m×10 m×10 m的長方體，球體直徑D=1 m，如圖2所示?？紤]到運動物體的曲面外形，而非結構網(wǎng)格能較好地適應曲面的變化，故模型的面網(wǎng)格采用三角網(wǎng)格，計算域采用四面體和六面體；為了同時保證計算速度和精度，減少動網(wǎng)格的數(shù)量，將整個計算域劃分為靜止區(qū)域和變化區(qū)域，劃分網(wǎng)格時在質量優(yōu)先的準則下對靜止區(qū)域使用四面體網(wǎng)格，對變化區(qū)域使用六面體網(wǎng)格[17]，在計算過程中，靜止區(qū)域的網(wǎng)格保持不變，變化區(qū)域的網(wǎng)格根據(jù)邊界的運動進行實時更新。

圖2　計算域示意圖Fig.2　Schematic diagram of computational region

為了兼顧模型曲面網(wǎng)格的精細度與總體計算網(wǎng)格數(shù)量，模型壁面網(wǎng)格尺寸為0.03 m，外圍邊界網(wǎng)格尺寸為0.3 m，最終網(wǎng)格數(shù)量在149萬，網(wǎng)格質量絕大部分超過0.7，如圖3所示。

圖3　網(wǎng)格劃分結果Fig.3　Mesh generation

2.3動態(tài)網(wǎng)格法

CFD數(shù)值計算中，整個計算域由網(wǎng)格組成，物體的運動意味著物面邊界的移動，即運動條件需要定義在邊界條件上，使用自定義函數(shù)(user-defined function，UDF)對其進行描述，UDF是使用C語言編寫的函數(shù)，使用Fluent.Inc提供的DEFINE宏定義，被求解器動態(tài)地加載。

隨著邊界的運動，部分網(wǎng)格被不斷壓縮，部分網(wǎng)格被不斷拉伸，壓縮和拉伸產(chǎn)生的畸變會影響計算的精度，畸變達到一定程度還會產(chǎn)生負體積網(wǎng)格而導致計算失敗，這種情況下就需要用到彈簧光順、動態(tài)層鋪、局部網(wǎng)格重構等網(wǎng)格更新方法。

彈簧光順法不改變網(wǎng)格的拓撲結構，能保證網(wǎng)格質量，適用于小位移問題；其他2種方法改變了網(wǎng)格的拓撲結構， 新網(wǎng)格的變量需要從舊網(wǎng)格中插值

得到。對于大位移和強切變問題僅靠節(jié)點松弛不能保證網(wǎng)格質量，而且有可能出現(xiàn)網(wǎng)格相交的情況，根據(jù)網(wǎng)格特點，采用彈簧光順法和局部重構法相結合的方法來解決。物面邊界運動后，首先使用彈簧光順法調(diào)整網(wǎng)格，當網(wǎng)格的品質參數(shù)如光順性、正交性、節(jié)點分布特性，高寬比和單元尺寸改變到一定程度，其畸變率超過設定值(本文為0.75)時，使用局部網(wǎng)格重構法，通過預置參數(shù)如網(wǎng)格的最小、最大長度，網(wǎng)格單元的偏斜度、尺度函數(shù)參數(shù)等控制網(wǎng)格移動后的品質。

3仿真校驗

本文驗證計算的浸深位置為自由水面距離球心2D(y=2 m)，根據(jù)工程經(jīng)驗，可以認為超過此浸深時，自由水面對球體結構附加質量的影響可以忽略不計。在UDF中設定不同的運動規(guī)律，速度v=[1,2]與加速度a=[3，-3，10，-10]進行自由組合，前者單位為m/s，后者單位為m/s2，即物體分別以(v,a)=(1,3)，(1,-3)，(1,10)，(1,-10)，(2,3)，(2,-3)，(2,10)，(2,-10)等組合進行運動。數(shù)值計算結束后，利用前述方法進行后處理，附加質量計算結果如表1所示，其中時間為球體的運動時間，在本文設定的運動規(guī)律下，球體短時間運動所產(chǎn)生的位移相對其浸深為小量，可以認為球體的浸深位置近似不變。多步計算平均結果與理論參考值相比的誤差如表2～4所示，平均步數(shù)n就代表對前n個時間步的計算結果進行平均。表1～4中，y代表垂向運動，z代表橫向運動，va[1,3,-3]代表結果是利用(v，a)=(1,3)，(1,-3)2組運動規(guī)律下的監(jiān)控數(shù)據(jù)進行計算的，其他以此類推。

表1　不同組合運動規(guī)律下的垂向和橫向附加質量計算結果

表2　不同平均步數(shù)下垂向附加質量的相對誤差(v=1)

表3　不同平均步數(shù)下垂向附加質量的相對誤差(v=2)

表4　不同平均步數(shù)下橫向附加質量的相對誤差

從表1可以看出，不同組合運動規(guī)律下，隨著計算時間的增加，速度增量的影響變大，附加質量的波動也隨之變大，故不宜過多。表2～4中，一步計算的相對誤差最大為0.700 02%，多步平均的計算結果波動比較平緩，相對誤差最大為3.46%，雖然計算結果普遍小于理論值，但誤差基本上都在可接受范圍內(nèi)，同時，也說明工程上可以選擇時間代價較小的一步計算結果。

為便于直觀對比，以垂向運動為例，繪制圖4和圖5，圖4為同一組加速度下不同初始速度的計算結果對比，圖5為同一初始速度下不同加速度的計算結果對比，其中，F(xiàn)low time代表流動時間。從圖4中容易看出，初始速度較大者的計算結果波動大，較小的相對平穩(wěn)。從圖5可以看出，加速度較大者波動較小，更接近理論值。

圖4　相同加速度下的垂向附加質量Fig.4　Vertical added mass in the same accelerated velocity

圖5　相同初速度下的垂向附加質量Fig.5　Vertical added mass in the same initial velocity

綜上，可以認為，本文的計算策略與所選的速度、加速度無關聯(lián)，是正確有效的，而選擇初始速度較小加速度較大的組合運動規(guī)律能夠降低計算結果的波動程度。

4近水面球體結構的附加質量

為考察近水面球體結構附加質量隨浸深的變化規(guī)律，定義浸深H為自由水面距離球心的長度，利用本文的方法計算浸深0.5D，0.6D，0.8D，1.0D，1.5D時球體的附加質量。

依據(jù)前述方法，不影響速度前提下盡可能選擇較大的加速度，設定[1,10,-10]，[2,10,-10]，[4,10,-10]3組運動規(guī)律來計算球體在近水面不同浸深下的附加質量，如圖6，7所示，圖中Step1代表一步計算結果，垂向、橫向附加質量的變化曲線對比如圖8所示。

圖6　隨浸深變化的垂向附加質量Fig.6　Depth varying curve of vertical added mass

圖7　隨浸深變化的橫向附加質量Fig.7　Depth varying curve of lateral added mass

圖8　垂向與橫向附加質量對比圖Fig.8　Comparison of vertical added mass and lateral added mass

從圖6，7可以看出，不同運動規(guī)律下一步計算結果之間吻合較好，再次說明設定的運動規(guī)律不影響計算結果；球體接近水面的過程中，雖然一直保持全浸沒狀態(tài)，浸濕外形不變，但由于自由面的存在，附加質量不再是常量，而是隨浸深呈下降趨勢，而且越接近自由面，變化越劇烈；球體剛好浸沒在水面下時，相較于深水時垂向、橫向附加質量已經(jīng)分別下降了29.6%和16.2%。

從圖8中可以看出，橫向與垂向的附加質量變化規(guī)律基本一致，均隨浸深減小呈下降趨勢，橫向附加質量在幅值上大于垂向，在球體頂端接觸自由面時達到最大，兩者相差19%。

5結束語

本文基于CFD平臺，綜合運用VOF和動網(wǎng)格技術，提出了一種近水面運動物體附加質量的求解方法，并參考球體的理論結果進行了驗證。在此基礎上，以球體為例，研究了其不同浸深下附加質量的變化規(guī)律，結果表明，由于自由水面的存在，雖然保持全浸沒狀態(tài)，但近水面時球體的垂向和橫向附加質量均呈下降趨勢，尤其在接近自由水面時更明顯，不利于近水面運動過程的控制。結論對于球體結構的穩(wěn)定性、安全性設計具有指導價值，對其他形狀的潛器也具有一定參考意義。本文未考慮自由面的波動，波浪對近水面物體附加質量的影響尚待進一步研究。

參考文獻：

[1]杜明明.兩層流體中潛浮式結構物的水動力特性研究[D].上海:上海交通大學,2007.

DU Ming-ming.Hydrodynamic Characteristics of Submerged Structures in a Two-Layer Fluid[D]. Shanghai:Shanghai Jiaotong University,2007.

[2]劉濤.深海載人潛水器耐壓球殼設計特性分析[J].船舶力學，2007，11(2):214-220．

LIU Tao.Reaearch on the Design of Spherical Pressure Hull in Manned Deep-Sea Submersible[J].Journal of Ship Mechanics,2007,11(2):214-220.

[3]Lopes DBS. Hydrodynamic Coefficients of a Submerged Pulsating Sphere in Finite Depth[J]. Ocean Engineering, 2002,29(1):1391-1398.

[4]Matiur Rahman. Simulation of Diffraction of Ocean Waves by a Submerged Sphere in Finite Depth[J]. Applied Ocean Research,2001,23(1):305-317.

[5]黃景泉，張宇文.魚雷流體力學[M].西安：西北工業(yè)大學出版社，1989.

HUANG Jing-quan, ZHANG Yu-wen. Torpedo Fluid Mechanics[M].Xi′an:Northwestern Polytechnic University Press,1989.

[6]李遠林，伍曉榕.非線性橫搖阻尼的實驗確定[J].華南理工大學學報：自然科學版，2002，30 (2):79-82.

LI Yuan-lin,WU Xiao-rong.Experimental Determination of Nonlinear Roll Damping[J].Journal of South China University of Technology：Natural Science ed,2002,30(2):79-82.

[7]David Valentin,Alexandre Presas,Eduard Egusquiza,et al.Experimental Study on the Added Mass and Damping of a Disk Submerged in a Partially Fluid-Filled Tank with Small Radial Confinement[J].Journal of Fluids and Structures,2014,50(1):1-17.

[8]WANG K,JIN Sh,LIU G.Numerical Modeling of Free-Surface Flows with Bottom and Surface-Layer Pressure Treatment[J].Journal of Hydrodynamics, 2009, 21(3):352-359.

[9]傅慧萍，李杰.附加質量CFD計算方法研究[J].哈爾濱工程大學學報，2011，32(2):148-152.

FU Hui-ping,LI Jie.Numerical Studies of Added Mass Based on the CFD Method[J].Journal of Harbin Engineering University,2011, 32(2):148-152.

[10]李剛，段文洋，郭志彬.復雜構型潛器附加質量的研究[J].哈爾濱工業(yè)大學學報，2010，42(7):1145-1148.

LI Gang, DUAN Wen-yang, GUO Zhi-bin.Added Mass of Submerged Vehicles with Complex Shape[J]. Journal of Harbin Institute of Technology, 2010,42(7):1145-1148.

[11]林超友，朱軍．潛艇近海底航行附加質量數(shù)值計算[J].船舶工程，2003，25(1):26-29.

LIN Chao-you, ZHU Jun.Numerical Computation of Added Mass of Submarine Maneuvering with Small Clearance to Sea-Bottom[J]. Ship Engineering,2003,25(1):26-29.

[12]Rajeev K Jaiman, Manoj K Parmar, Pardha S Gurugubelli.Added Mass and Aeroelastic Stability of a Flexible Plate Interacting with Mean Flow in a Confined Channel[J].Journal of Applied Mechanics,2014,81(1):1-9.

[13]朱仁傳,郭海強,繆國平,等.一種基于CFD理論船舶附加質量與阻尼的計算方法[J].上海交通大學學報,2009,43(2):198-203.

ZHU Ren-chuan, GUO Hai-qiang, MIAO Guo-ping,et al.A Computational Method for Evaluation of Added Mass and Damping of Ship Based on CFD Theory[J].Journal of Shanghai Jiaotong University,2009, 43(2):198-203.

[14]王恭義.物體在帶自由表面的流體中作任意運動的附加質量[D].大連:大連理工大學,2007.

WANG Gong-yi. Added Mass for Bodies Moving Arbitrarily in Fluid with Free Surface[D].Dalian:Dalian University of Technology,2007.

[15]車得福，李會雄.多相流及其應用[M].西安：西安交通大學出版社，2007:546-549.

CHE De-fu,LI Hui-xiong. Multiphase Flow and Its Application[M].Xi′an:Xi′an Jiaotong University Press,2007:546-549.

[16]CHEN X F,LI Y C,WANG Y X,et al.Numerical Simulation of Wave Interaction with Perforated Caisson Breakwaters[J].China Ocean Engineering, 2003, 17(1):33-43.

[17]褚江.非結構動網(wǎng)格生成方法研究[D].南京:南京理工大學,2006.

CHU Jiang.Research of the Generation of Dynamic Unstructured Meshes[D].Nanjing:Nanjing University of Science and Technology,2006.

Added Mass of Spherical Structure during Its Moving Process Towards Water Surface

ZHANG Xiao-qiang1, WANG Guan2, FENG Jin-fu1

(1.AFEU，Aeronautics and Astronautics Engineering College,Shaanxi Xi’an 710038, China；2. PLA,No.95889 Troop，Beijing 100085,China)

Abstract:Aiming at the added mass of spherical structure during its moving process towards water surface, a new method is proposed by using volume of fluid(VOF) model and dynamic mesh technology based on computational fluid dynamics (CFD) platform. The mathematical formula is derived, and then, numerical model is built, and the validation and analysis of the proposed method is promoted through the results of different motion laws. At last, the added mass of spherical structure in different depth is computed; the conclusions shows that, the proposed method proves correct; the vertical and horizontal added mass of spherical structure are declining during its moving towards free surface, when the top of spherical structure touches the free surface, the vertical added mass drops by 29.6% and horizontal drops by 16.2%, and the effect of the changes on its designation of stabilization and safety shouldn't be neglected.

Key words:spherical structure；added mass；free surface；hydrodynamic force；fluid-structure interaction；dynamic mesh；numerical computation

*收稿日期：2015-10-23；修回日期：2015-11-18

作者簡介：張曉強(1984-)，男，河南南陽人。工程師，博士生，主要從事新概念武器設計與仿真。

通信地址：710038陜西省西安市灞橋區(qū)東里村霸陵路1號二部軍械教研室E-mail：18001283156@189.cn

doi:10.3969/j.issn.1009-086x.2016.02.034

中圖分類號：TB126；U661

文獻標志碼：A

文章編號：1009-086X(2016)-02-0215-07