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

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

2. 91267 部隊(duì)，福建 福州 350000）

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

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

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

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

對(duì)于一個(gè)以速度Vp在Z 軸方向上運(yùn)動(dòng)的剛性球體來(lái)說(shuō)，為了計(jì)算它的入水載荷，首先需要寫(xiě)出其動(dòng)力學(xué)方程。根據(jù)牛頓第二定律，剛性球體的動(dòng)力學(xué)方程為：

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

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

目前計(jì)算入水阻力的一個(gè)普遍采用的計(jì)算公式為：

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

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

2 剛性球體垂直入水載荷計(jì)算模型

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

如圖2 所示是四分之一球體剖面由OMM′運(yùn)動(dòng)到O′N(xiāo)N′。

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

圖 2 球體入水運(yùn)動(dòng)示意圖Fig. 2 General view of sphere’s micromovement

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

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

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

過(guò)圓心O 分別與Mi、Mi+1連線(xiàn)并延長(zhǎng)，與弧線(xiàn)NN'分別相交于Ni、Ni+1，如圖4 所示。由于流體無(wú)黏且不可壓，可認(rèn)為弧段MiMi+1通過(guò)擴(kuò)張到達(dá)NiNi+1。采用此等效方法，會(huì)使得水面附近的弧段在弧線(xiàn)NN'上無(wú)對(duì)應(yīng)區(qū)域，當(dāng)L 極小時(shí)，此部分弧長(zhǎng)趨于0，可將其忽略。

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

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

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

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

圖 4 微元邊界運(yùn)動(dòng)等效示意圖Fig. 4 General view of microboundary’smotion equivalent

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

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

δmθ的動(dòng)能增量δEθ可以表示為：

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

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

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

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

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

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

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

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

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

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

3 模型驗(yàn)證

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

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

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

式中：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ù)修正項(xiàng)，為1.393 7；E 為體積內(nèi)能。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

由圖9～11 可知，本文所建立的理論模型得到入水載荷峰值結(jié)果偏小，峰值出現(xiàn)時(shí)間偏晚，由于該理論模型未考慮入水沖擊過(guò)程中產(chǎn)生的水堆，水堆的出現(xiàn)使得球體周?chē)刃嫣?，?dǎo)致入水沖擊載荷峰值時(shí)間的提前。同時(shí)由于理論模型均勻化了流體內(nèi)部復(fù)雜的壓力波動(dòng)，因此計(jì)算結(jié)果也較為平緩，峰值寬度也較大。

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

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

4 變參數(shù)分析

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

4.1 變質(zhì)量分析

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

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

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

4.2 變體積分析

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

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

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

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

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

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

4.3 變速度分析

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

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

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

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

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

5 結(jié) 論

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

（1）基于無(wú)黏不可壓流體流動(dòng)模型，考慮流體彈性，采用邊界運(yùn)動(dòng)等效方法對(duì)運(yùn)動(dòng)邊界進(jìn)行分段分析對(duì)剛性球體入水載荷進(jìn)行計(jì)算是可行的；

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

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

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