剛性球體三維高速垂直自由入水載荷*

孫玉松，周穗華，張曉兵，孫玉明

（1. 海軍工程大學兵器工程學院，湖北 武漢 430033；

2. 91267 部隊，福建 福州 350000）

水上飛機、航天器、魚雷、水雷等入水均屬非周期性入水問題，有著廣泛的應用前景，是流體力學領域內研究的重要內容。自1929 年von Karman 提出忽略浮力，采用動量定理對楔形體以及平頭體入水載荷進行計算后，Wagner[1]將von Karman 的方法進行了更為詳細的推導，考慮到水面抬升，提出了基于小斜升角模型的近似平板理論。May 等[2]分析了入水初期流體流動特點，得出了鋼球入水初期的入水載荷系數。Korobkin 等[3]考慮結構體的彈性，在每一個時間步內首先求解結構體受力，而后根據結構體壓力分布求解其變形，是一種耦合度較高的計算方法。關于入水初期結構體載荷問題的理論研究基本遵循Wagner 漸進匹配近似理論，該理論主要針對可等效為二維結構體的入水模型，當流體三維流動特性較為明顯時，需引入修正系數對該方法進行修正[4-5]。宋保維等[6]和王永虎等[7]建立了六自由度的水雷水下彈道模型，依據不可壓縮流體非定常勢流理論和鏡像法求解入水載荷，并應用于大型彈體以約60 m/s 速度入水的計算中。孫士麗等[8]采用三維全非線性不可壓縮勢流理論方法研究了有限水深中非軸對稱體的斜向入水抨擊問題，對軸對稱體與非軸對稱體的垂直入水以及斜入水進行了模擬。

隨著數值方法的發展并在流體流動問題上的應用，工程上大量復雜流動問題得以解決，王健等[9]、馬慶鵬等[10]和朱珠等[11]使用基于網格的數值方法對結構體入水過程進行了仿真，并對入水載荷進行了分析。Oger 等[12]采用光滑粒子流體動力學方法（SPH），計及入水過程中流體的可壓縮性，計算了作用在剛體上的壓力。雖然數值方法已能有效解決大部分的工程問題，但對于從本質上解釋入水沖擊這一現象的意義有限。本文中所研究的剛性球體以160～240 m/s 的速度垂直入水，流體的運動是一個非定常的三維流動，且持續時間極短，撞水瞬間流體呈現出較強的彈性效應，并產生振動[13]。由于入水載荷極大且持續時間極短，重力作用和流體黏性作用均可忽略。液體不可壓縮，但需考慮其彈性，認為小幅度的彈性對于流體的不可壓縮性影響很小，在無黏不可壓流體流動模型的基礎上，采用微元邊界運動等效方法對運動邊界進行分段分析，計及入水過程中系統的動能損失，根據能量守恒，對剛性球體垂直入水初期流體的三維運動進行分析，求出剛性球體入水過程中的載荷。

1 剛性球體垂直入水動力學方程

剛性球體在入水初始時刻的位置示意圖如圖1 所示。

對于一個以速度Vp在Z 軸方向上運動的剛性球體來說，為了計算它的入水載荷，首先需要寫出其動力學方程。根據牛頓第二定律，剛性球體的動力學方程為：

圖 1 剛性球體位置示意圖Fig. 1 General view of the rigid sphere’s location

式中：m 為剛性球體的質量，F 為其入水阻力。

目前計算入水阻力的一個普遍采用的計算公式為：

式中：ρ 為水的密度，Av為剛性球體的最大截面積，Cd為與速度相關的阻力系數。

根據Charters 及Thomas 測量的阻力系數與入水速度的關系可以分為3 個區間，入水速度在亞聲速到Ma=0.5 區間時，Cd=0.384。式(2)在剛性球體對周圍流體形成完全擾動之后是有效的，但在擾動形成階段，該公式計算得到的結果嚴重偏小，Wagner 理論正是用于計算該過程的，自提出以來，圍繞該理論取得了較為豐碩的成果[14]，但基于Wagner 理論的方法在討論三維流動問題時往往是按照二維入水問題計算，然后添加修正系數[15]，本文針對Wagner 理論所存在的不足，采用剛性球體微元邊界運動等效的方法考慮三維效應，對剛性球體的入水載荷進行計算。

2 剛性球體垂直入水載荷計算模型

剛性球體高速垂直入水時，背部將會形成空泡，假設入水深度過半球時流體與球體分離，未觸水部分無表面阻力。從球體觸水至觸水面發展完全這一撞水階段，由于撞水時間持續極短，流體黏性作用有限，因此將流體看作無黏流體，Lee 等[16]在研究高速入水條件下的彈道波的時候使用可壓縮波動方程和不可壓縮伯努利方程（非定常流動的伯努利方程），雖然看起來是前后矛盾的，但依然取得了非常好的結果。可以認為，在剛性球體高速入水的條件下，流體中有限的彈性對于流體的不可壓縮特性影響不大，本文模型是建立在無黏不可壓彈性流體基礎上的。

如圖2 所示是四分之一球體剖面由OMM′運動到O′NN′。

將弧邊MM′進行分割，任意兩相鄰點Mi 與Mi+1 組成微元邊界MiMi+1，MiMi+1 所代表的微元面在空間上的位置如圖3 所示。

圖 2 球體入水運動示意圖Fig. 2 General view of sphere’s micromovement

圖 3 弧段微元示意圖Fig. 3 General view of a microsegmental arc

圖中陰影部分即是弧段MiMi+1所代表的空間上的微元面si，弧段MiMi+1即是剛性球體截面圓的微元邊界，為了表達的方便，下面以微元邊界MiMi+1指代微元面si進行表述。

微元面s 的面積為dSi，表達式為：式中：δh 為微元邊界MiMi+1在豎直方向上的投影；θ 為OMi與OM 的夾角；R 為球體半徑。

過圓心O 分別與Mi、Mi+1連線并延長，與弧線NN'分別相交于Ni、Ni+1，如圖4 所示。由于流體無黏且不可壓，可認為弧段MiMi+1通過擴張到達NiNi+1。采用此等效方法，會使得水面附近的弧段在弧線NN'上無對應區域，當L 極小時，此部分弧長趨于0，可將其忽略。

設MiMi+1表面流體的等效擴張速度大小為vi，表達式為：

微元邊界MiMi+1在δt 時刻后擾動的流體區域示意圖如圖5 所示。

CiCi+1為在微元邊界MiMi+1作用下流體受迫運動產生的擾動波波面，波面的擴張速度為聲速c，波面CiCi+1后部的流體處于未擾動的狀態，速度為0。根據彈性理論，波的傳播速度為聲速c，則此刻波面距離運動的球心O 的距離Rci的表達式為：

式中：t 為系統時間；ti為微元邊界MiMi+1觸水時間。經過δt 后，新增擾動區域內流體質量δmθ的表達式為：

圖 4 微元邊界運動等效示意圖Fig. 4 General view of microboundary’smotion equivalent

圖 5 微元邊界擴張示意圖Fig. 5 General view of microsegmental arc’s expansion

當δt 足夠小時，可認為新增擾動區域內流體速度大小相等，設該區域內流體平均速度為vci，根據不可壓縮流體流動的連續性，可知vci的表達式為：

δmθ的動能增量δEθ可以表示為：

根據動量守恒，可知δmθ在垂直方向上的動量增量δMθ可以表示為：

設δt 時間內作用在δmθ上的力為Fci，平均值為 Fci，根據動量定理，則 Fci可以表示為：

在微觀層面，Fci是變化的，但作用面的位移速度為vci，則力Fci在δt 時間內做的功δWi可以采用平均力 來計算，即：

在δt 時間內，剛性球體的每一個微元邊界對流體做的功是流體動能增加量的兩倍，損失的能量將主要以波的形式不斷擴散出去。根據能量定律，在δt 時間內，剛性球體的動能損失量dE 應等于各個微元邊界對流體做功的總和，dE 的表達式為：

δt 時間內，剛性球體受到的合外力F 的表達式為：

剛性球體入水載荷，即入水加速度a 的表達式為：

隨著剛性球體入水深度的增大，以對流體擾動作用為主的Wagner 階段逐漸結束，本文所建立的模型便不再適用，后期剛性球體入水阻力F 的表達式為：

式中：Cd是隨速度變化的，在本文中，Cd=0.384。

從Wagner 階段結束到可使用經典入水阻力表達式之間的過渡階段流體流動較為復雜，這里進行近似處理，當采用經典入水阻力表達式計算得到的球體入水阻力Fc大于本文模型計算得到的入水阻力FR時，便使用經典的入水阻力表達式計算入水載荷。

綜上所述，剛性球體垂直入水載荷表達式為：

3 模型驗證

剛性球體高速入水是強非線性過程，涉及固、液、氣三相的運動，采用多介質ALE 方法對剛性球體入水沖擊過程進行數值計算。

為了節省計算量，根據對稱性，采用四分之一模型對入水過程進行計算，計算域網格采用八節點六面體單元，空氣和水均采用ALE 網格，剛性球體采用Lagrangian 網格，四分之一計算域有限元模型如圖6 所示。

定義球體為剛體，水和空氣均選用空材料模型，水的狀態方程采用Grüneisen 方程，空氣則采用線性多項式方程。Grüneisen 方程的壓力表達式為：

式中：C=1 480 m/s 為介質中聲速；S1、S2、S3為沖擊波輸入參數，S1=2.56，S2=-1.986，S3=0；μ 為介質壓縮比，μ=ρ/ρ0-1，ρ0為常溫狀態下水初始密度，ρ 為水當前密度；γ0=0.493 4 為Grüneisen 初系數；a 為Grüneisen 系數修正項，為1.393 7；E 為體積內能。

圖 6 有限元模型Fig. 6 Finite element model

線性多項式（Linear_polynomial）狀態方程的壓力表達式為：

式中：μ=ρ/ρ0-1，E 為體積內能，C0=C1=C2=C5=C6=0，C3=C4=0.4。

設置過球體剖面的兩個面為對稱邊界面，其余面為非反射邊界面。采用罰函數方法進行流固耦合。球體半徑為0.2 m，質量為261 kg，入水速度為200 m/s。經過計算得到球體入水載荷曲線如圖7所示。

對于有限元仿真結果的驗證一般采用試驗的方法，但大型球體高速入水所需動能極大，實驗室內難以達到相應的試驗條件，且在高速入水條件下，結構體入水載荷難以準確測量，而入水空泡的外形容易觀測，入水空泡與能量轉化量直接相關，仿真模型的準確性可以通過對比空泡的外形來驗證[17,18]。在入水空泡的研究方面，May 等[19]通過大量的試驗觀察得到入水初期入水空泡輪廓除頭部以外可近似為拋物線，顧建農等[20]對此進行了進一步驗證，拋物線的方程為：

式中：h 為入水深度；CD為結構體入水阻力系數；db為入水結構體特征尺寸，本文中為球體直徑。

在當前入水條件下，基于有限元模型和May的空泡模型得到的入水空泡輪廓對比結果如圖8所示。

可知除觸水部分外，兩種方法計算得到的空泡輪廓一致，因此使用該有限元方法和參數計算剛性球體入水沖擊過程是準確可行的，可使用該有限元方法對式（16）所表示的理論方法進行驗證。

有限元方法與理論方法分別得到的剛性球體入水載荷曲線對比結果如圖9 所示。

理論方法與有限元方法分別求得的入水載荷峰值相差15.2%，入水載荷峰值出現的時間相差32%。為進一步驗證理論方法的合理性，下面使用兩種方法分別對質量為391.5 kg 以及入水速度為240 m/s時鋼球的入水載荷進行計算，結果見圖10 和圖11。

此入水條件下，入水載荷峰值相差14.4%，入水載荷峰值出現的時間相差29.4%。

此入水條件下，入水載荷峰值相差18.7%，入水載荷峰值出現的時間相差18.2%。

圖 8 入水空泡輪廓對比圖Fig. 8 Comparison of the cavities

圖 9 入水載荷對比圖Fig. 9 Comparison of water-entry impact

圖 10 球體質量為391.5 kg 時入水載荷對比圖Fig. 10 Comparison of water-entry impact when the sphere’s mass is 391.5 kg

圖 11 入水速度為240 m/s 時入水載荷對比圖Fig. 11 Comparison of water-entry impact when initial velocity is 240 m/s

由圖9～11 可知，本文所建立的理論模型得到入水載荷峰值結果偏小，峰值出現時間偏晚，由于該理論模型未考慮入水沖擊過程中產生的水堆，水堆的出現使得球體周圍等效水面抬升，導致入水沖擊載荷峰值時間的提前。同時由于理論模型均勻化了流體內部復雜的壓力波動，因此計算結果也較為平緩，峰值寬度也較大。

同時可以看出本文所建立的理論方法得到的入水載荷有一個明顯的拐點，這是由于隨著入水深度的增加，剛性球體對周圍流體的擾動已不是入水阻力形成的主要原因，剛性球體入水過程中阻力所做的功主要轉化為排開流體的動能，入水阻力主要成因的轉換使得入水過程存在一個過渡的區域，本文對過渡區域進行了簡化，當采用經典入水阻力方程計算得到的入水阻力較大時，便采用經典入水阻力方程計算入水載荷，由于計算模型的不同，導致了拐點的出現。

但總體來看，理論方法和有限元方法計算結果一致性較好，該理論模型能夠反映出剛性球體入水過程中載荷的主要影響因素。

4 變參數分析

為探究影響入水載荷的因素，本文采用控制變量法，計算不同質量，不同體積，不同入水速度條件下的入水載荷并進行對比分析。

4.1 變質量分析

對體積和入水速度相同，質量分別為130.5、261 和391.5 kg 的剛性球體入水過程中的入水載荷進行計算，并對結果進行比較，結果如圖12所示。

由圖12 可知，隨著剛性球體質量的增加，入水載荷峰值將顯著降低，峰值持續時間變化較小。不同質量條件下剛性球體入水載荷峰值如表1 所示。

由表1 可知，其它條件不變時，入水載荷峰值與剛性球體質量呈反比關系。

4.2 變體積分析

將密度和入水速度相同，半徑分別為0.1、0.2 和0.3 m 的剛性球體代入模型中進行計算，由于體積（半徑）的改變，密度相同的條件下，剛性球體質量也會發生變化，計算得到的3 個剛性球體入水載荷如圖13 所示。

圖 12 入水載荷對比圖Fig. 12 Comparison of water-entry impact

表 1 入水載荷峰值對比Table 1 Comparison of the water-entry peak impact

由圖13 可知，隨著體積的增大，剛性球體入水載荷峰值降低，且峰值持續時間顯著增大。不同體積條件下入水載荷峰值如表2 所示。

表 2 入水載荷峰值對比Table 2 Comparison of the water-entry peak impact

由表2 可知，在剛性球體密度不變的條件下，其入水載荷峰值與半徑成反比例關系。

4.3 變速度分析

將質量和體積均相同，入水速度v 分別為160、200 和240 m/s 的剛性球體代入計算模型，得到3 種不同入水速度條件下剛性球體入水載荷如圖14 所示。

由圖14 可知隨著入水速度的增大，入水載荷峰值的增大比較顯著，但峰值持續時間會有所縮短。不同入水速度條件下剛性球體入水載荷峰值如表3 所示。由表3 可知，入水載荷峰值與入水速度的平方近似成線性關系。

圖 13 入水載荷對比圖Fig. 13 Comparison of water-entry impact

圖 14 入水載荷對比圖Fig. 14 Comparison of water-entry impact

表 3 入水載荷峰值對比Table 3 Comparison of the water-entry peak impact

5 結 論

本文中基于無黏不可壓流體流動模型，考慮流體彈性，采用微元邊界運動等效方法對運動邊界進行分段分析，計及入水過程中系統的動能損失，根據能量守恒定理，對剛性球體高速垂直自由入水過程中流體的三維流動進行了理論分析，建立了基于無黏不可壓彈性流體的剛性球體垂直高速入水載荷計算模型，并驗證了該方法的可行性。基于此模型，進一步分析了剛性球體質量、半徑以及入水速度對入水載荷的影響，結果表明：

（1）基于無黏不可壓流體流動模型，考慮流體彈性，采用邊界運動等效方法對運動邊界進行分段分析對剛性球體入水載荷進行計算是可行的；

（2）入水載荷峰值與剛性球體質量呈反比關系。隨著剛性球體質量的增大，峰值持續時間變化較小；

（3）入水載荷峰值與等密度剛性球體的半徑成反比例關系。隨著剛性球體體積的增大，入水峰值降低，且峰值持續時間顯著增長；

（4）入水載荷峰值與剛性球體入水速度的平方近似成線性關系。隨著入水速度的增大，峰值持續時間會有所縮短。