王永虎，石秀華

（1.中國(guó)民航飛行學(xué)院，四川 廣漢618307；2.西北工業(yè)大學(xué) 航海學(xué)院，西安710072）

自從上世紀(jì)20年代末水上飛機(jī)的水面降落問(wèn)題提出以來(lái)，結(jié)構(gòu)物入水沖擊的理論與試驗(yàn)研究得到了飛速發(fā)展，例如，水上飛機(jī)和宇宙飛船水面著陸、衛(wèi)星海面回收、空投雷彈入水沖擊、船舶在風(fēng)浪中砰擊和救生艇的海上拋落等[1，2].由于理論分析入水初始彈道比較困難，忽撲行為的研究也只局限于定性分析，所以入水初始彈道和忽撲在空投雷彈入水沖擊問(wèn)題中的研究就顯得尤為重要.20世紀(jì)60年代，WAUGH J G在假雷和實(shí)雷實(shí)驗(yàn)基礎(chǔ)上，對(duì)空投雷彈入水彈道和忽撲行為進(jìn)行了研究[2，3].為了研究空投水雷的入水以及水下運(yùn)動(dòng)行為機(jī)理，達(dá)到更準(zhǔn)確地布雷，MANN J L采用蒙特卡羅法獲得三維水雷模型 MINE6D的運(yùn)動(dòng)軌跡[4].CHU P C等建立時(shí)域3DOF運(yùn)動(dòng)方程以及回歸模型，對(duì)圓柱體跌落入水進(jìn)行了理論分析和實(shí)驗(yàn)研究[5，6].PARK Man－sung等[7]利用面元法數(shù)值計(jì)算了正切尖拱體以任意角度入水的沖擊載荷，然后對(duì)尖拱體的忽撲行為進(jìn)行了數(shù)值計(jì)算分析，但是在數(shù)值計(jì)算中忽略了入水空泡對(duì)入水彈道的影響，計(jì)算結(jié)果不可避免地存在誤差.由于空投雷彈斜入水的復(fù)雜性和特殊性，必須對(duì)其入水初始彈道進(jìn)行準(zhǔn)確分析以確保斜入水的航行安全性.

1 入水沖擊顯式動(dòng)力模型

用數(shù)理解析方法和數(shù)值方法研究入水初始彈道問(wèn)題時(shí)，必須把結(jié)構(gòu)動(dòng)力學(xué)方程和流體動(dòng)力學(xué)方程耦合起來(lái)求解，解析方法只適應(yīng)于簡(jiǎn)化模型，而數(shù)值模擬應(yīng)用范圍更廣.LS－DYNA是以顯式為主、隱式為輔的非線性動(dòng)力有限元仿真軟件，程序中的單元采用拉格朗日列式增量解法，特別適合求解各種二維、三維非線性動(dòng)力沖擊問(wèn)題，并與實(shí)驗(yàn)進(jìn)行無(wú)數(shù)次的對(duì)比，證實(shí)了模擬計(jì)算的可靠性[8].下面首先根據(jù)多物質(zhì)流固耦合方法，建立斜入水沖擊顯式動(dòng)力模型；然后數(shù)值模擬入水沖擊過(guò)程，最后應(yīng)用到入水初始彈道問(wèn)題的研究中.

1.1 固體域網(wǎng)格劃分

在入水沖擊的顯式動(dòng)力模型中，采用有限元方法離散入水結(jié)構(gòu)體（固體域），其有限元方程為

式中，M為總質(zhì)量矩陣；C為結(jié)構(gòu)阻尼系數(shù)；P為總體載荷矢量；F為單元應(yīng)力場(chǎng)等效節(jié)點(diǎn)力矢量組；H為總體結(jié)構(gòu)沙漏粘性阻尼力（t）為總體節(jié)點(diǎn)加速度矢量（t）為總體節(jié)點(diǎn)速度矢量.

為了真實(shí)模擬入水沖擊的液面隆起現(xiàn)象，并提高數(shù)值計(jì)算速度，求解算法上采用LS－DYNA程序中的多物質(zhì)ALE算法.由于無(wú)反射邊界只能施加到實(shí)體單元上，入水結(jié)構(gòu)體采用LS－DYNA單元庫(kù)中的六面體SOLID164體元模擬.為了與水域、空氣域耦合，固體域的網(wǎng)格劃分采用Lagrange網(wǎng)格算法.如果結(jié)構(gòu)體外形特殊（例如尖拱體），需采用布爾運(yùn)算劃分網(wǎng)格，對(duì)尖拱體進(jìn)行了八節(jié)點(diǎn)體元網(wǎng)格劃分.同時(shí)，為了節(jié)省計(jì)時(shí)，采用沙漏粘性阻尼方式控制單點(diǎn)積分引起沙漏模式，這對(duì)模擬液面隆起和飛濺等大變形現(xiàn)象十分有效.由于只討論入水初始彈道，所以將空投雷彈結(jié)構(gòu)外形簡(jiǎn)化為圖1所示的數(shù)值模型.

在數(shù)值計(jì)算過(guò)程中，采用LS－DYNA的中心差分時(shí)間積分的顯式方法，計(jì)算系統(tǒng)各節(jié)點(diǎn)在第n時(shí)間步時(shí)刻的加速度a、速度v和位移s分別為

式中，F(xiàn)ext為節(jié)點(diǎn)外力和體力矩陣，F(xiàn)int為節(jié)點(diǎn)內(nèi)力矩陣.

圖1 空投雷彈斜入水沖擊數(shù)值模型及計(jì)算域

1.2 流體域網(wǎng)格劃分

空投雷彈入水沖擊涉及到三相介質(zhì)：固體域、空氣域和水域.空氣域和水域均采用LS－DYNA中的空材料模式＊MAT＿NULL來(lái)描述，即用本構(gòu)模型和狀態(tài)方程來(lái)同時(shí)描述流體材料，通過(guò)Gruneisen狀態(tài)方程確定壓力－體積的關(guān)系式：

式中，ρ0為材料密度；vc為沖擊波速度；E0為單位體積內(nèi)能；S1、S2、S3、γ0為材料常數(shù)；μ為密度變化率；α為一階體積修正.在LS－DYNA中采用＊EOS＿GUNEISEN關(guān)鍵字定義，其中水和空氣2種流體采用Euler網(wǎng)格建模，建模時(shí)水和空氣介質(zhì)均使用映射方法劃分網(wǎng)格.流體材料本構(gòu)模型和狀態(tài)方程具體參數(shù)分別見(jiàn)表1，本文數(shù)值模型采用cm－g－μs單位制，流體材料本構(gòu)模型和狀態(tài)方程具體參數(shù)分別見(jiàn)表1.

表1 流體材料的本構(gòu)模型和狀態(tài)方程的相關(guān)參數(shù)

1.3 流固耦合ALE算法

本文采用ALE算法，結(jié)合Euler算法與Lagrange算法的優(yōu)點(diǎn)，結(jié)構(gòu)采用Lagrange單元算法，流體采用Euler/ALE多物質(zhì)單元算法.ALE算法先執(zhí)行一個(gè)或幾個(gè)Lagrange時(shí)步計(jì)算，單元網(wǎng)格隨材料流動(dòng)而產(chǎn)生變形.然后，為了保持變形后的物質(zhì)邊界條件，對(duì)內(nèi)部單元重新劃分網(wǎng)格，將變形網(wǎng)格中的單元變量（密度、能量、應(yīng)力張量等）和接點(diǎn)速度矢量輸運(yùn)到重分后的新網(wǎng)格中，最后執(zhí)行ALE時(shí)步計(jì)算.特點(diǎn)是耦合面隨著Lagrange單元和Euler單元形成的共同邊界產(chǎn)生運(yùn)動(dòng)變形，Euler單元也就可以隨著結(jié)構(gòu)大范圍移動(dòng)，保證了數(shù)值計(jì)算順利進(jìn)行，且大大節(jié)約計(jì)時(shí).

LS－DYNA程序中的多物質(zhì)ALE－Lagrange算法可以傳遞ALE網(wǎng)格中的流體材料和Lagrange結(jié)構(gòu)體間的接觸力，能方便地通過(guò)＊CONSTRAINED＿LAGRANGE＿IN＿SOLID關(guān)鍵字把流體和固體單元進(jìn)行耦合，且建模時(shí)流體網(wǎng)格和固體網(wǎng)格可以交叉重疊.通過(guò)＊SECTION＿SOILD＿ALE關(guān)鍵字來(lái)定義單元算法類(lèi)型并標(biāo)識(shí)相關(guān)單元算法.為了更接近模擬無(wú)限水域的分析情況，在流體單元的邊界上定義無(wú)反射邊界條件來(lái)簡(jiǎn)化入水沖擊模型.

1.4 網(wǎng)格劃分尺寸控制

流體網(wǎng)格密度對(duì)數(shù)值結(jié)果是有影響的，理論上講網(wǎng)格劃分越細(xì)，數(shù)值計(jì)算結(jié)果越接近真實(shí)值，但實(shí)際中細(xì)密的網(wǎng)格劃分會(huì)造成數(shù)值計(jì)算中誤差累積、結(jié)果偏離和計(jì)時(shí)延長(zhǎng).考慮到計(jì)算機(jī)允許的計(jì)算能力和工程應(yīng)用的精度要求，合理確定網(wǎng)格密度和計(jì)算規(guī)模是非常必要的.

為了確定最佳的網(wǎng)格密度，并捕捉應(yīng)力場(chǎng)的最大梯度變化，這里采用一個(gè)無(wú)量綱參量來(lái)描述網(wǎng)格密度，定義Euler－Lagrange單元網(wǎng)格尺寸比為

式中，LEu為Euler單元網(wǎng)格特征尺寸；LLag為L(zhǎng)agrange單元網(wǎng)格特征尺寸.一般而言，尺寸比的選取以在仿真過(guò)程中數(shù)值收斂性好，計(jì)算結(jié)果較精確為準(zhǔn).這里對(duì)Rm＝0.4～0.9的情況進(jìn)行了數(shù)值計(jì)算，計(jì)算結(jié)果見(jiàn)圖2，圖中Cd為入水沖擊阻力系數(shù).

圖2 不同Rm的計(jì)算精度曲線

由圖2可以看出，隨著尺寸比的減少，曲線震蕩明顯減弱，而且入水沖擊載荷峰值也在減少，當(dāng)Rm＝0.6時(shí)曲線趨于平穩(wěn)，當(dāng)Rm＝0.4時(shí)得到最穩(wěn)定的曲線.數(shù)值計(jì)算效益和經(jīng)濟(jì)性不僅要考慮結(jié)果的精確度，而且要看計(jì)算消耗的機(jī)時(shí)和存取容量.雖然Rm＝0.4的結(jié)果較精確，但計(jì)時(shí)較長(zhǎng)和存儲(chǔ)容量驟然增大.因此，這里在空投雷彈入水沖擊數(shù)值計(jì)算中，Rm取0.5～0.6較合理，故本文后續(xù)的彈道分析均采用Rm＝0.6.

2 計(jì)算結(jié)果與分析

為了便于模型的建立和高效的數(shù)值求解，這里主要分析空投雷彈二維入水運(yùn)動(dòng)學(xué)行為，即入水初始彈道，將結(jié)構(gòu)體定義為剛性體以減少自由度的數(shù)量，每次數(shù)值仿真的固體域的尺寸和網(wǎng)格劃分方法都不變.為了保證入水彈道軌跡和射流飛濺大部分出現(xiàn)在流體域中，根據(jù)入水初始彈道時(shí)間和水平距離確定流體域尺寸，這樣會(huì)導(dǎo)致每次數(shù)值計(jì)算中網(wǎng)格數(shù)量不一致，但不影響數(shù)值仿真結(jié)果的準(zhǔn)確性.

尖拱體斜入水的初始條件和物理參數(shù)為：根據(jù)空投雷彈的實(shí)際質(zhì)量，設(shè)定尖拱體質(zhì)量為18.4kg，局部坐標(biāo)系中質(zhì)心的轉(zhuǎn)動(dòng)慣量在K文件中設(shè)定為＊DIM，Ivert，ARRAY，6，1，1，，，、＊SET，IVERT（1，1，1），5189233、＊SET，IVERT （4，1，1），2312606、＊SET，IVERT（6，1，1），5189233，單位為cm－g－μs.數(shù)值計(jì)算得到各情況下的入水初始彈道，并通過(guò)二維圖來(lái)表示.

首先，為了研究入水角對(duì)忽撲行為的影響，在初始入水速度為100m/s，初始攻角為0，對(duì)入水角分別以10°、20°和30°的初始彈道進(jìn)行模擬.采用尖拱頂點(diǎn)運(yùn)動(dòng)軌跡來(lái)形象顯示水下初始彈道，見(jiàn)圖3所示.可以看出，入水角越小，運(yùn)動(dòng)軌跡越偏離初始軌道；也表明入水角越小，尖拱體越容易產(chǎn)生忽撲行為甚至彈跳，與文獻(xiàn)[7]采用面元法的數(shù)值結(jié)果吻合，究其原因是入水角小的尖拱體在入水過(guò)程中受到的縱傾力矩大.

然后，對(duì)質(zhì)量為18.4kg尖拱體以固定的30°入水角且零攻角情況下的斜入水進(jìn)行數(shù)值仿真，選擇入 水 初 速 度v0分 別 為 30.48 m/s，60 m/s 和100m/s，分析初始速度對(duì)入水初始彈道的影響，見(jiàn)圖4所示，圖中顯示的是二維Oxy坐標(biāo)系中頂點(diǎn)運(yùn)動(dòng)軌跡.可以看出，雖然入水速度不同，但只要攻角和入水角不變，尖拱頂點(diǎn)的運(yùn)動(dòng)軌跡是相同的，即入水角對(duì)入水初始彈道的影響微乎其微.圖5給出了不同質(zhì)量的尖拱體在v0＝100m/s和相同時(shí)間情況下的入水初始彈道二維圖，可以看出，隨著入水體質(zhì)量的增加，在相同時(shí)間內(nèi)，質(zhì)量小的尖拱體入水初始彈道曲率較大，從而很容易產(chǎn)生忽撲或彈跳行為.

圖3 以相應(yīng)入水角入水的初始彈道

圖4 入水速度不同的入水初始彈道

圖5 質(zhì)量不同的入水初始彈道

最后，空投雷彈與流體剛剛接觸的時(shí)候具有一定的入水攻角，這里數(shù)值分析了初始攻角對(duì)入水運(yùn)動(dòng)初始軌跡的影響，見(jiàn)圖6所示.分別選取了3種情況進(jìn)行模擬仿真，即入水攻角α0分別為－5°、0°和5°，入水初始速度為100m/s，入水角為固定的30°.

從仿真結(jié)果可以看出，在入水初速度和入水角速度一定的條件下，入水初始攻角對(duì)入水初始運(yùn)動(dòng)軌跡和姿態(tài)角變化影響較大.正攻角入水使雷彈往拉平方向運(yùn)動(dòng)，加劇忽撲甚至出現(xiàn)彈跳行為；負(fù)攻角入水使雷彈具有俯沖下潛的趨勢(shì)，起到抑制的作用.所以，在小角度斜入水過(guò)程中應(yīng)注意入水攻角的影響作用.

圖6 初始攻角不同的水下初始彈道

3 結(jié)論

采用ANSYS/LS－DYNA大型非線性顯式動(dòng)力程序進(jìn)行空投雷彈斜入水初始彈道模擬研究，實(shí)現(xiàn)入水空泡生成和擴(kuò)展現(xiàn)象以及不正常入水造成的忽撲行為模擬，考慮了初始條件和物理?xiàng)l件對(duì)斜入水初始彈道的影響，得到以下結(jié)論：①入水角對(duì)入水初始彈道影響很大，如果入水角很小可能產(chǎn)生忽撲行為甚至彈跳行為，這對(duì)空投雷彈而言是不利的影響因素；②入水初速度對(duì)入水初始彈道幾乎沒(méi)有什么影響，空投雷彈都是沿著近似相同的彈道運(yùn)行；③空投雷彈自身的質(zhì)量對(duì)入水初始彈道也有影響，在其他條件相同的條件下質(zhì)量越小越容易產(chǎn)生忽撲行為；④對(duì)帶有攻角姿態(tài)而言，攻角大小對(duì)初始彈道有不同的影響，特別針對(duì)小角度斜入水情況，可以考慮采用適當(dāng)攻角來(lái)彌補(bǔ)小角度帶來(lái)的負(fù)面影響，此入水初始彈道規(guī)律將指導(dǎo)后續(xù)的入水彈道定性或定量分析.

[1]WANG Y，SHI X，WANG P.Dynamical response analysis of incautious water entry of UUV based on exact body shape approach[C].7th World Congress on Intelligent Control and Automation.Chongqing，China：IEEE，2008：4 870－4 875.

[2]WAUGH J G.Water－entry pitch modeling[J].Journal of Hydronautics，1968，2（2）：87－92.

[3]WAUGH J G，STUBSTAD G W.Hydroballistics modeling，AD A007529[R].1975.

[4]MANN J L.Deterministic and stochastic modeling of the water entry and descent of three－dimensional cylinderical bodies[D].USA：Massachusetts Institute of Technology，2005.

[5]CHU P C，F(xiàn)AN C W.3Drigid body impact burial prediction model[J].Advances in Mechanics，2004，5：43－52.

[6]CHU P C，F(xiàn)AN C W.Pseudo－cylinder parameterization for mine impact burial prediction[J].Journal of Fluids Engineering，2005，127：1 515－1 520.

[7]PARK Man－sung，JUNG Young－rae.Numerical study of impact force and ricochet behavior of high speed water－entry bodies[J].Computers ＆ Fluids，2003，32：939－951.

[8]HALLQUIST J O.LS－DYNA theoretical manual[M].Livermore，California：Livermore Software Technology Corporation，2006.