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

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

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

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

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

1 剛性球體垂直入水動力學(xué)方程

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

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

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

式中：m 為剛性球體的質(zhì)量，F(xiàn) 為其入水阻力。

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

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

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

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

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

如圖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 所示。由于流體無黏且不可壓，可認(rèn)為弧段MiMi+1通過擴張到達NiNi+1。采用此等效方法，會使得水面附近的弧段在弧線NN'上無對應(yīng)區(qū)域，當(dāng)L 極小時，此部分弧長趨于0，可將其忽略。

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

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

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

式中：t 為系統(tǒng)時間；ti為微元邊界MiMi+1觸水時間。經(jīng)過δt 后，新增擾動區(qū)域內(nèi)流體質(zhì)量δmθ的表達式為：

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

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

當(dāng)δt 足夠小時，可認(rèn)為新增擾動區(qū)域內(nèi)流體速度大小相等，設(shè)該區(qū)域內(nèi)流體平均速度為vci，根據(jù)不可壓縮流體流動的連續(xù)性，可知vci的表達式為：

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

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

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

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

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

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

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

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

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

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

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

3 模型驗證

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

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

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

式中：C=1 480 m/s 為介質(zhì)中聲速；S1、S2、S3為沖擊波輸入?yún)?shù)，S1=2.56，S2=-1.986，S3=0；μ 為介質(zhì)壓縮比，μ=ρ/ρ0-1，ρ0為常溫狀態(tài)下水初始密度，ρ 為水當(dāng)前密度；γ0=0.493 4 為Grüneisen 初系數(shù)；a 為Grüneisen 系數(shù)修正項，為1.393 7；E 為體積內(nèi)能。

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

線性多項式（Linear_polynomial）狀態(tài)方程的壓力表達式為：

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

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

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

式中：h 為入水深度；CD為結(jié)構(gòu)體入水阻力系數(shù)；db為入水結(jié)構(gòu)體特征尺寸，本文中為球體直徑。

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

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

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

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

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

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

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

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

圖 10 球體質(zhì)量為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 可知，本文所建立的理論模型得到入水載荷峰值結(jié)果偏小，峰值出現(xiàn)時間偏晚，由于該理論模型未考慮入水沖擊過程中產(chǎn)生的水堆，水堆的出現(xiàn)使得球體周圍等效水面抬升，導(dǎo)致入水沖擊載荷峰值時間的提前。同時由于理論模型均勻化了流體內(nèi)部復(fù)雜的壓力波動，因此計算結(jié)果也較為平緩，峰值寬度也較大。

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

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

4 變參數(shù)分析

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

4.1 變質(zhì)量分析

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

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

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

4.2 變體積分析

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

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

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

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

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

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

4.3 變速度分析

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

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

圖 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 結(jié) 論

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

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

（2）入水載荷峰值與剛性球體質(zhì)量呈反比關(guān)系。隨著剛性球體質(zhì)量的增大，峰值持續(xù)時間變化較??；

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

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