電磁力控制潛艇近水面興波繞流流場的數值研究

紀延亮，周本謀，黃亞冬

（南京理工大學 瞬態物理重點實驗室，南京 210094）

0 引 言

近現代的戰爭中以潛艇為代表的水下作戰平臺由于高度隱蔽性和優異的作戰性能，得到了各國的重視，相關方面的研究也一直經久不衰。譬如DAPRA使用Suboff[1]系列計算模型進行研究，得到了大量實驗數據，為水下航行器的外形設計提供了寶貴經驗。Serhat hosder等[2]使用風洞對該模型靜態和機動兩種情況的壁面湍流進行了實驗分析，得到不同偏航角度時尾部流動分離結構及壁面受力系數的分布情況。隨著數值算法的不斷完善，常規實驗難以實現的工況也能得到模擬，同時節約大量的研究經費。Gross等[3]以Suboff裸艇體為研究對象，從數值計算和實驗兩方面著手研究了不同雷諾數和攻角情況下艇身近壁區域的流動分離現象；Alin等[4]使用LES方法對潛艇機動過程中繞流流場結構進行了仿真研究，通過實驗對比說明其采用的數值方法可以很好地仿真各種工況的潛艇繞流流場。考慮到潛艇在水下幾百米的深水區域中活動，之前的研究主要基于無界繞流狀態，然而靠近水面位置，鈍體繞流會產生復雜的水面興波，興波的存在導致潛艇的阻力大大增加，同時穩定性變差。隨著潛艇作戰范圍的不斷擴展，近水面航行性能也需要關注，例如在復雜海況下，潛艇的穩定性受到波浪力的干擾會變差；潛艇潛射導彈時，需要上浮到近水面保持低速穩定航行。對于水面艦船的興波問題研究已有百余年歷史，早在19世紀八十年代，Kelvin就對靜水中的壓力點源興波進行了研究，后來基于Havelock源和Rankine源的格林函數求解得到了廣泛應用并且不斷發展，同時運用數值仿真方法使得興波問題分析時更加直觀全面[5-8]。潛艇近水面興波問題已有相關研究工作，Zhang Nan等[9-11]使用RANS模型對潛艇不同航行工況進行了系統的研究，將海底以及自由面的影響進行了完整分析；Tim Gourlay等[12]使用Havelock面元法研究了不同弗勞德數對應的近水面潛艇繞流興波問題，結果均顯示潛艇阻力較水下航行時有不同程度的增加。

由于近水面興波的不利影響較為突出，因此可以借助一些措施來改善這種現象，若改變外形結構，則會影響其整體水下機動性能，因此采用主動控制方法更加合適，常見的主動控制包括壁面吹吸法、側面注射法、電磁流動控制[13-15]等方法。其中電磁力具有場力的結構傳輸特性，能夠在不改變原來流場的邊界條件、不需要向弱導電流場（如鹽水、海水等）傳輸質量的情況下，比較方便地向流場傳輸動量、能量以及渦量，從而可以有效地改變和構造流體邊界層與流場的結構。Weier等[16]研究發現流向電磁力可以減小繞流物體的阻力，增加層流-湍流轉捩臨界雷諾數；張輝等[17]通過圓柱繞流實驗驗證了電磁力在減阻消渦方面的可行性。Liu等[18]模擬了無界繞流中潛艇的流場結構和受力情況，提出潛艇指揮臺結構是阻力和渦街脫落的主要來源，通過對指揮臺施加流向電磁力達到減阻消渦的目的。

本文使用VOF方法對潛艇近水面區域的繞流興波問題開展研究工作，基于高度函數法構建空間興波自由面，分析了潛艇在不同潛深下的興波特性和繞流流場結構，研究近水面興波擾動對潛艇整體受力的影響。通過局部電磁力控制，達到減小阻力，抑制尾流渦街，提高穩定性的目的，為近水面潛艇的機動性能和隱身能力的優化提供參考。

1 計算方法

1.1 兩相流控制方程

近水面繞流問題涉及到水和空氣交界面，屬于兩相流問題，因此需要引入VOF模型進行求解。對于不可壓，密度不均勻的流場，含有表面張力的Navier-Stokes方程表示為：

其中：ρ為流體密度，μ 為動力粘度系數，Dij=（?iUj+?jUi）/2為形變張量，狄拉克分布函數 δS表示該表面張力項作用于自由面，σ和κ分別表示表面張力系數和曲率，n代表自由面法向量。

對于水—氣兩相流動，設單元中水的體積分數為c，則流場的密度和粘性可表示為：

其中：ρ和μ下標1，2分別代表水和空氣兩相，以c表示的對流方程改寫為：

設來流速度為U∞，對計算時間及流場參數進行無量綱化處理：

其中：帶有上標*的參數為實際流場計算值。計算時的各受力系數表示為：

其中：F為計算受力，S為潛艇濕表面積。

1.2 計算模型及邊界條件

潛艇的外形多種多樣，本文選取帶附體的回轉體結構為代表對該類潛行設備的繞流問題進行研究（圖1）。主體部分為軸向中心對稱的回轉體，頂部指揮臺為橢圓柱形。以艇身弦長L為參考長度，主體直徑為0.12L。指揮臺高度0.06L，橢圓面長半軸和短半軸尺寸分別為0.045L和0.02L，尾翼采用翼型設計。計算空間采用笛卡爾坐標系，定義坐標軸方向X，Y，Z分別為流向，橫向和垂直方向。指揮臺和頭部交錯排布電磁極，從而能夠產生沿壁面的流向電磁力。

圖1 潛艇模型Fig.1 Submarine model

初始時刻靜止水面坐標為Z=0，潛艇下潛過程中回轉體前駐點始終在（0，0，Z）位置移動。流場入口處為混合邊界條件，入口流向速度為U∞，壓力為0，流入流體組分保持恒定。出口以及其余四個側面為諾依曼邊界條件，即壓力、速度分量的法向梯度為零，以消除邊界的堵塞效應。潛艇表面為固壁邊界，計算時設為無滑移邊界條件。

1.3 空間和時間離散

計算區域如圖2所示，流場計算區域為4L×2L×2L（L為艇身長度），潛艇壁面初始網格劃分層數為N=9，即最小尺寸為L/29，計算開始后通過笛卡爾網格進行自適應動態劃分，流場內形成分層組織的八叉樹結構，算法采用Khokhlov的全線程數據結構，可以有效遍歷不同類型的單元網格。網格自適應依據局部渦量為準則，當滿足（7）式時網格被加密：

其中：a為網格長度，U為速度矢量，ε為介于0到1之間的常數。自適應判定周期為每一時間步，當局部參數不滿足該式時，網格將被粗糙化，本文設定ε為0.01。

圖2 空間網格劃分Fig.2 The scope of calculation field and the spatial division

采用時間分裂投影法對控制方程進行時間離散[19]，動量方程的對流項采用二階迎風的Bell-Colella-Glaz格式進行離散，這種格式對于CFL數小于1是穩定的。擴散項采用隱式Crank-Nicholson方法離散，它具有二階精度且無條件穩定。因此時間和空間離散均具有二階精度。

使用分段線性VOF方法構建兩相交界自由面，同時采用二階精度的高度函數估算曲率，結合平衡力表面張力離散方法，非平衡的自由面形狀有足夠的時間松弛為平衡形態，平衡形態以二階速率向精確值收斂。

1.4 電磁力源項

數值計算時電磁體積力以源項的形式加入到動量方程中，電磁力表現為電流密度矢量J和磁場強度B相互作用：

F=J×B

由歐姆定律得到：

J=σ（ E+U× ）B

其中：σ、E分別為電導率和電場強度。對于低速流動，磁場中帶電粒子流動所激勵的電流要遠遠小于外加電場，因此U×B可以忽略，電磁力的產生僅與外加電場和磁場有關，且假設電磁力沿壁面法向呈指數衰減。施加于表面的電磁力大小通過無量綱化的作用系數N表示，N為電磁力與流體慣性力的比值：

其中：J0，B0分別為壁面上的電流密度和磁場強度值。

1.5 算法驗證

以潛艇模型Suboff為例，選取裸艇體對其無界繞流和水面興波問題進行算法驗證。無界繞流時對壁面網格進行不同層數的劃分，以艇身長度為參考的計算雷諾數為1.2×107。計算足夠時間周期后截取潛艇中剖面上的壓力系數與實驗結果[1]進行對比如圖3所示。不同的網格層數壓力變化趨勢一致，當網格劃分層數n=6和n=7時，較為粗糙的網格使得壓力系數在局部有較大的波動。提高劃分層數后，所得結果與實驗數值有較好的吻合，對潛艇近壁流場使用n=9的劃分層數時，可以滿足潛艇壁面附近的流場計算精度要求。

圖3 不同網格劃分層次的Suboff無界繞流壁面流向壓力分布Fig.3 The streamwise pressure distribution of the infinite deep flow around Suboff under various grid division

同時，為檢驗算法用于自由面興波的計算精度，對Suboff模型在水面航行的情況進行仿真計算，求得Fr在0.17和0.27時的阻力系數與實驗和文獻結果[9]進行對比。如表1所示，兩種Fr數下平均阻力系數與文獻中實驗值和模擬結果接近，最大誤差分別為3%和2.4%。這里認為該算法適用于自由面興波問題的研究。

表1 水面航行潛艇興波數值結果驗證Tab.1 Validation test of the wave-making resistance for the submarine on free surface

2 計算與結果分析

2.1 不同下潛深度流場分析

文獻[18]對潛艇無界繞流的流場結構及受力情況進行了詳細分析，本文以潛艇主體軸中心線為基準標定下潛深度h，以Fr=0.2為例進行潛艇近水面不同下潛深度的興波流場分析。從水面上方觀察，興波波型分布如圖4所示，當h=0.17時興波效應較為明顯，興波由兩部分組成：潛艇頭部興波和指揮臺興波，頭部由于離水面較遠，興波效果不明顯。指揮臺雖然尺寸小，但是更靠近水面，因此下游的大尺度興波以指揮臺興波為主。興波以一定的夾角向后方傳播，傳播的外側波動呈現出周期性起伏，而在指揮臺后方區域，由于尾部脫落渦的影響，興波波幅被抑制，后方三角區域內無明顯興波。圖5為潛艇不同潛深時中軸線上方的瞬時波高分布，越靠近水面繞流興波越明顯，三種潛深情況下指揮臺前后的興波波幅均為最大，進一步說明指揮臺是興波的主要來源。h=0.17L時，指揮臺后方興波的最大峰谷值分別為0.94%L和1.23%L，興波繼續向后傳播，波型出現了較為明顯的擾動，興波能量被耗散，波動幅度被明顯削弱。潛深較大時波型則無明顯擾動。

圖4 潛深h=0.17的興波波型Fig.4 Wave pattern for h=0.17

圖5 三種潛深潛艇中軸線上方的水面興波高度Fig.5 Wave patterns for different depths of submergence above the central axis

觀察兩種潛深的水下流場結構，圖6為t=6時刻艇身周圍不同截面的三維渦量分布，截取位置為0.3-1.2L，間隔距離0.2L。指揮臺的脫落渦街主要影響艇身上方區域，尾渦沿流向脫落后擾動先增強后逐漸衰減，艇身后方則以艇艉脫落渦街為主。當h=0.17時，指揮臺后方0.5L截面附近的渦量變化最為劇烈，較h=0.27時有所提前，且受水面起伏的波浪影響，指揮臺后方脫落的渦街有豎直向上發展的趨勢；而當h=0.27時，水下脫落的渦街受水面影響較小，渦量沿流向的發展相對滯后且豎直向的抬升不明顯。通過上述對比可知，興波使得靠近水面的繞流流場渦量發展更迅速，這就導致指揮臺尾部區域的流場變化更加劇烈，影響了潛艇整體的受力分布。對指揮臺附近的局部流場結構進行分析，選取t=6時刻的流場分布如圖7、圖8所示，與h=0.27時的繞流流場結構相比，h=0.17時的流場結構更加復雜，指揮臺頂部前緣的流動分離被抑制，后方會形成傾斜向上的低速低壓帶狀區域。這主要是由于指揮臺興波豎直方向的波動導致近水面區域流體的速度、壓力產生強烈擾動，當潛艇靠近水面時，這種擾動對壁面產生影響。與此同時指揮臺的尾流與其相互疊加，使得后方尾流區域的三維流場結構更加復雜。而隨著繼續下潛，興波效應不明顯，興波面豎直方向的擾動對潛艇繞流流場結構的影響也越來越小。

圖6 不同潛深水下三維渦量分布Fig.6 Underwater 3-D vorticity for different depths of submergence

圖7 不同潛深指揮臺附近的流向速度分布Fig.7 The distribution of streamwise velocity around appendage for various depth of submergence

圖8 不同潛深指揮臺附近的壓力系數分布Fig.8 The distribution of the pressure coefficient around appendage for various depth of submergence

對不同潛深的潛艇受力系數進行分析，計算濕表面積為整個潛艇表面，取值約為0.361 3。總阻力系數為Cfd，包括壓差阻力和摩擦阻力兩部分，不同潛深的平均阻力系數如表2所示。選取三種潛深潛艇繞流計算穩定后的阻力系數取均值與無界繞流情況下的阻力系數進行對比，距離水面越近，所受阻力越大，并且增長速率也有所加快。尤其是h=0.17時，阻力系數的增幅達到了18.4%，可見水面興波對近水面潛行的潛艇總阻力影響是非常明顯的。在h=0.27的位置時，雖然興波效果不再明顯，但是阻力系數仍有4.1%的增加。

表2 不同潛深阻力對比Tab.2 Comparison of the resistance for various depth of submergence

圖9表示了h=0.17的潛艇繞流計算過程中受力系數的變化情況。阻力，橫向力及垂向力系數分別表示為Cfd、Cfh和Cfv。阻力系數在t=2后就趨于穩定，圍繞其均值上下波動。同時橫向力和垂向力系數有較為強烈的波動，且二者均呈現出無序性。由此可見，靠近水面航行時，潛艇受到的各方向力的擾動有所增加，說明興波對潛艇的穩定性有明顯影響，表現為較大幅度的不規則波動，這些問題為近水面的活動帶來了不利影響。

圖9 h=0.17時潛艇各方向受力計算收斂曲線Fig.9 Force coefficient convergence history for h=0.17

3 電磁力流動控制效果

上文分析了潛艇水面及近水面航行時的水面興波及整體受力情況，本節針對Fr=0.2，潛深h=0.17的潛艇繞流采取近壁面電磁力控制的方式，研究艇身受力及周圍流場變化情況。已知潛行時指揮臺和潛艇頭部均有水面興波產生，這里對指揮臺和潛艇頭部的繞流流動進行控制。與整體加力相比，局部作用力的施加在實際過程中更加易于實現，同時，潛艇壁面的摩擦阻力也不會增加過多。設電磁力同時作用于潛艇頭部和指揮臺兩部分，作用時刻始于t=4。電磁力作用后，水面興波波型有所變化，圖10為N=30時中軸線上方水面興波波型分布，指揮臺后方波面擾動被抑制，重新呈現出有規則的起伏波動。

圖10 N=30時潛深h=0.17的水面波高對比Fig.10 Comparison of the wave pattern for h=0.17,N=30

圖11為電磁力強度N=30時水下渦量等值線分布，與圖6相比，指揮臺后方流向X渦量有所增加，然而Y、Z兩方向的渦量則被明顯抑制，說明豎直向和橫向的擾動減小。艇身壁面流體由于電磁力的作用，動量增加，壁面剪切運動增強，因此Y、Z渦量值有所增加。從圖12中可見加力前后潛艇周圍流場速度和壓力的變化，從速度分布圖可以看出流向電磁力將能量導入近壁流體中，使得近壁流體動量增加，壁面周圍形成了一層速度較大的剪切流動。加力后指揮臺頂部及后方的流動分離被抑制，斜上方的低壓條帶也有明顯減小，同時指揮臺后方和艇艉區域的壓力有明顯升高。

圖11 N=30時水下渦量分布Fig.11 The underwater 3-D vorticity at N=30

電磁力作用后艇身各方向受力變化如圖13所示，當N為10和30時，減阻效果相似，這主要是由于電磁力的增加后，雖然壓差阻力出現下降，但是潛艇壁面附近的流動速度增加，相應的摩擦阻力也會上升，二者綜合之后，使得兩種工況下的減阻效果相似。隨著電磁力進一步增強，壓差阻力的下降值將明顯大于壁面摩阻的增量，因此總阻力還是呈下降的趨勢，當N=50時阻力最多下降了24.1%。圖13（b）、（c）中反映了控制后橫向力和垂向力的變化情況，當N超過30后，橫向力的波動被很好地抑制住。豎直方向的垂向力也隨著電磁力強度的增加而逐漸向0值靠近。由此可見，電磁力將能量注入近壁流體，流體產生定向運動，這種定向運動改變了潛艇周圍的流場結構，削弱了上方水面波浪產生的擾動影響，阻力和穩定性均有所改善。

圖13 不同電磁力作用系數對各向受力系數的影響Fig.13 Variation of spatial force coefficient for different N

上文已知指揮臺是主要的興波源，并且靠近水面的指揮臺周圍的流場結構更加復雜。若潛艇所受的擾動主要集中在指揮臺附近，則理論上對其單獨控制也會得到不錯的減阻減振效果。這里作為對比，將潛艇頭部和指揮臺進行單獨控制，研究各部位的控制效果。定義指揮臺部位的控制模式為A，單獨控制頭部為模式B，作用系數選取N=30，計算結果如圖14所示。可以看出，減阻方面模式B要明顯優于模式A，同時橫向和垂向受力控制也是模式B起主導作用。主要由于B模式誘導出的流體對整個艇身形成包覆，因此對壁面的實際控制范圍要遠遠大于模式A，同時也說明近水面指揮臺引起的興波結構對潛艇整體產生了干擾。

圖14 不同電磁力作用位置對各向受力系數的影響Fig.14 Variation of spatial force coefficient for different controlling locations

4 結 語

本文對潛艇近水面航行時的興波現象和艇身周圍的流場結構進行了研究，分析了Fr=0.2時不同下潛深度對應的興波波幅以及艇身周圍繞流流場結構的差異。當潛艇靠近水面時，興波效應增強，指揮臺后方受到水面波浪起伏的影響，渦量變化更加劇烈，近壁區域的流場結構更加復雜，尾部脫落發卡渦與水面興波疊加，使得水面波動能量被耗散，規則的興波傳播被破壞，同時艇身阻力顯著增加，整體穩定性降低。使用近壁流向電磁力對頭部和指揮臺控制后，潛艇后方的發卡渦被抑制，隨著電磁力強度的增加，潛艇總阻力明顯下降，當N達到50時，潛艇橫向受力波動也得到了很好的抑制。對潛艇頭部和指揮臺進行獨立控制，對比發現頭部控制效果要好于指揮臺控制，說明指揮臺引起的興波對潛艇整個艇身產生了影響。綜上所述，電磁力的控制作用可以消除興波的影響，優化各類潛行設備近水面航行時的受力和穩定性。

參 考 文獻：

[1]Huang T T,Liu H L,Groves N C.Experiments of the DARPA(Defense Advanced Research Projects Agency)Suboff Program[J].Experiments of the Darpa Suboff Program,1989.

[2]Hosder S,Simpson R L.Experimental investigation of unsteady flow separation on a maneuvering axisymmetric body[J].Journal of Aircraft,2015,44(4)：1286-1295.

[3]Gross A,Kremheller A,Fasel H F.Simulation of flow over suboff bare hull model[C].49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition(AIAA),2011-290.

[4]Alin N,Fureby C,Svennberg U,et al.Large eddy simulation of the transient flow around a submarine during maneuver[C].International Conference on Numerical Ship Hydrodynamics,2013.

[5]Tarafder M S,Suzuki K.Numerical calculation of free-surface potential flow around a ship using the modified Rankine source panel method[J].Ocean Engineering,2008,35(5-6)：536-544.

[6]Zhang Baoji.Shape optimization of bow bulbs with minimum wave-making resistance based on rankine source method[J].Journal of Shanghai Jiaotong University：science,2012,17(01)：65-69.

[7]Yong Z,Zhi Z,Li Z.Ship hull optimization based on wave resistance using wavelet method[J].Journal of Hydrodynamics,2015,27(2)：216-222.

[8]Peng H,Ni S,Qiu W.Wave pattern and resistance prediction for ships of full form[J].Ocean Engineering,2014,87(5)：162-173.

[9]Zhang N,Ying L M,Yao H Z,et al.Numerical simulation of free surface viscous flow around submarine[J].Journal of Ship Mechanics,2005,9(3)：29-39.

[10]張 楠,沈泓萃,姚惠之.潛艇近海底與近水面繞流數值模擬研究[J].船舶力學,2007,11(4)：498-507.Zhang N,Shen H C,Yao H Z.Numerical simulation of flow around submarine operating close to the bottom or near surface[J].Journal of Ship Mechanics,2007,11(4)：498-507.

[11]張 楠,沈泓萃,姚惠之.潛艇阻力與流場的數值模擬與驗證及艇型的數值優化研究[J].船舶力學,2005,9(1)：1-13.Zhang N,Shen H C,Yao H Z.Validation of numerical simulation on resistance and flow field of submarine and numerical optimization of submarine hull form[J].Journal of Ship Mechanics,2005,9(1)：1-13.

[12]Gourlay T,Dawson E.A havelock source panel method for near-surface submarines[J].Journal of Marine Science and Application,2015,14(3)：215-224.

[13]Mamori H,Iwamoto K,Murata A.Effect of the parameters of traveling waves created by blowing and suction on the relaminarization phenomena in fully developed turbulent channel flow[J].Physics of Fluids,2014,26(1)：L73-L76.

[14]Majumder S,Sanyal D.Relaminarization of axisymmetric turbulent flow with combined axial jet and side injection in a pipe[J].J Fluids Eng,2010,132(10)：101101.

[15]Berger T W,Kim J,Lee C,et al.Turbulent boundary layer control utilizing the Lorentz force[J].Phys Fluids,2000,12(3)：631-649.

[16]Weier T,Gerbeth G,Mutschke G,et al.Control of flow separation using electromagnetic forces[J].Flow Turbul Combust,2003,71(1-4)：5-17.

[17]張 輝,范寶春,陳志華,等.Lorentz力對弱導電流體中圓柱受力的影響[J].水動力學研究與進展,2007,22(6)：766-773.Zhang H,Fan B C,Chen Z H,et al.The influences of Lorentz forces on the drag and lift coefficient of a circular cylinder in a low-conducting electrolytic fluid[J].Journal of Hydrodynamics,2007,22(6)：766-773.

[18]Liu H,Zhou B,Liu Z,et al.Numerical simulation of flow around a body of revolution with an appendage controlled by electromagnetic force[J].Proceedings of the Institution of Mechanical Engineers Part G Journal of Aerospace Engineering,2013,227(2)：303-310.

[19]Popinet S.Gerris：a tree-based adaptive solver for the incompressible Euler equations in complex geometries[J].Journal of Computational Physics,2003,190(2)：572-600.