王永虎，石秀華

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

自從上世紀(jì)20年代末水上飛機(jī)的水面降落問題提出以來，結(jié)構(gòu)物入水沖擊的理論與試驗(yàn)研究得到了飛速發(fā)展，例如，水上飛機(jī)和宇宙飛船水面著陸、衛(wèi)星海面回收、空投雷彈入水沖擊、船舶在風(fēng)浪中砰擊和救生艇的海上拋落等[1，2].由于理論分析入水初始彈道比較困難，忽撲行為的研究也只局限于定性分析，所以入水初始彈道和忽撲在空投雷彈入水沖擊問題中的研究就顯得尤為重要.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ù)值方法研究入水初始彈道問題時(shí)，必須把結(jié)構(gòu)動(dòng)力學(xué)方程和流體動(dòng)力學(xué)方程耦合起來求解，解析方法只適應(yīng)于簡(jiǎn)化模型，而數(shù)值模擬應(yīng)用范圍更廣.LS－DYNA是以顯式為主、隱式為輔的非線性動(dòng)力有限元仿真軟件，程序中的單元采用拉格朗日列式增量解法，特別適合求解各種二維、三維非線性動(dòng)力沖擊問題，并與實(shí)驗(yàn)進(jìn)行無數(shù)次的對(duì)比，證實(shí)了模擬計(jì)算的可靠性[8].下面首先根據(jù)多物質(zhì)流固耦合方法，建立斜入水沖擊顯式動(dòng)力模型；然后數(shù)值模擬入水沖擊過程，最后應(yīng)用到入水初始彈道問題的研究中.

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算法.由于無反射邊界只能施加到實(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ì)算過程中，采用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來描述，即用本構(gòu)模型和狀態(tài)方程來同時(shí)描述流體材料，通過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ù)分別見表1，本文數(shù)值模型采用cm－g－μs單位制，流體材料本構(gòu)模型和狀態(tài)方程具體參數(shù)分別見表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)體間的接觸力，能方便地通過＊CONSTRAINED＿LAGRANGE＿IN＿SOLID關(guān)鍵字把流體和固體單元進(jìn)行耦合，且建模時(shí)流體網(wǎng)格和固體網(wǎng)格可以交叉重疊.通過＊SECTION＿SOILD＿ALE關(guān)鍵字來定義單元算法類型并標(biāo)識(shí)相關(guān)單元算法.為了更接近模擬無限水域的分析情況，在流體單元的邊界上定義無反射邊界條件來簡(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ǎng)格密度，定義Euler－Lagrange單元網(wǎng)格尺寸比為

式中，LEu為Euler單元網(wǎng)格特征尺寸；LLag為L(zhǎng)agrange單元網(wǎng)格特征尺寸.一般而言，尺寸比的選取以在仿真過程中數(shù)值收斂性好，計(jì)算結(jié)果較精確為準(zhǔn).這里對(duì)Rm＝0.4～0.9的情況進(jìn)行了數(shù)值計(jì)算，計(jì)算結(jié)果見圖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ì)算得到各情況下的入水初始彈道，并通過二維圖來表示.

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

然后，對(duì)質(zhì)量為18.4kg尖拱體以固定的30°入水角且零攻角情況下的斜入水進(jìn)行數(shù)值仿真，選擇入 水 初 速 度v0分 別 為 30.48 m/s，60 m/s 和100m/s，分析初始速度對(duì)入水初始彈道的影響，見圖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)初始軌跡的影響，見圖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ì)，起到抑制的作用.所以，在小角度斜入水過程中應(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ì)入水初始彈道幾乎沒有什么影響，空投雷彈都是沿著近似相同的彈道運(yùn)行；③空投雷彈自身的質(zhì)量對(duì)入水初始彈道也有影響，在其他條件相同的條件下質(zhì)量越小越容易產(chǎn)生忽撲行為；④對(duì)帶有攻角姿態(tài)而言，攻角大小對(duì)初始彈道有不同的影響，特別針對(duì)小角度斜入水情況，可以考慮采用適當(dāng)攻角來彌補(bǔ)小角度帶來的負(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.