基于ALE，CEL和SPH方法的球形破片高速?zèng)_擊充液結(jié)構(gòu)對(duì)比研究

張 宇， 王彬文， 劉小川

(1.中國(guó)飛機(jī)強(qiáng)度研究所，西安 710065；2.結(jié)構(gòu)沖擊動(dòng)力學(xué)航空科技重點(diǎn)實(shí)驗(yàn)室，西安 710065)

1 引 言

易損性研究是提高飛機(jī)作戰(zhàn)生存力的基礎(chǔ)研究之一。其中油箱是飛機(jī)最大的易損結(jié)構(gòu)，約占飛機(jī)總體易損面積的75%[1]。飛機(jī)油箱毀傷包括水錘效應(yīng)和燃爆效應(yīng)，前者通過(guò)破壞油箱結(jié)構(gòu)完整性，后者通過(guò)引燃/引爆航空煤油，導(dǎo)致飛機(jī)失去作戰(zhàn)能力。

針對(duì)破片高速?zèng)_擊充液結(jié)構(gòu)研究，目前常用的研究方法有試驗(yàn)研究和數(shù)值模擬研究。其中試驗(yàn)研究需要大量的人力和物力，國(guó)外開(kāi)展了部分工作[2-8]。隨著計(jì)算機(jī)技術(shù)的不斷進(jìn)步，流固耦合問(wèn)題的數(shù)值仿真能力也逐漸完善。球形破片高速?zèng)_擊充液結(jié)構(gòu)是典型的大變形流固耦合過(guò)程，常用的分析方法包括隨機(jī)拉格朗日-歐拉法ALE(Arbitrary Lagrangian-Eulerian)、耦合歐拉-拉格朗日法CEL(Coupled Euler-Lagrange)和光滑粒子流體動(dòng)力學(xué)法SPH(Smoothed Particle Hydrodynamics)等[9]。

其中，Varas等[10]基于SPH和ALE方法，建立破片高速?zèng)_擊充液結(jié)構(gòu)有限元模型，研究了沖擊速度和充液比例對(duì)充液結(jié)構(gòu)毀傷效應(yīng)的影響，發(fā)現(xiàn)ALE方法中破片速度變化誤差更大，但在相同的網(wǎng)格密度下ALE方法得到的壓力精度更高；Sparks等[11]對(duì)比研究SPH和CEL方法模擬破片高速?zèng)_擊充液結(jié)構(gòu)，發(fā)現(xiàn)SPH方法可以用更少的單元實(shí)現(xiàn)對(duì)沖擊過(guò)程的模擬，但計(jì)算時(shí)間較長(zhǎng)；Artero-Guerrero等[12]采用CEL方法模擬破片高速?zèng)_擊碳纖維復(fù)合材料充液結(jié)構(gòu)；Pineda等[13]采用SPH-ALE結(jié)合的方法，研究了近壁面水中氣泡壓縮過(guò)程，并討論了SPH-ALE方法在這一特殊問(wèn)題上的局限性和前景；Sridhar等[14]針對(duì)金屬切削問(wèn)題，對(duì)比研究了ALE，CEL和SPH三種方法的可靠性；韓璐等[15]采用ALE方法研究了多枚破片高速?zèng)_擊對(duì)充液結(jié)構(gòu)的毀傷效應(yīng)；文獻(xiàn)[16-19]還基于ALE方法，研究充液比例對(duì)破片高速?zèng)_擊充液結(jié)構(gòu)的影響。此外，文獻(xiàn)[20-23]分別基于SPH和ALE方法對(duì)比研究復(fù)材平板入水沖擊、圓盤水漂、高速射流破土以及接觸爆炸等流固耦合問(wèn)題。

目前，對(duì)于破片高速?zèng)_擊充液結(jié)構(gòu)數(shù)值分析，學(xué)者們采用SPH，ALE和CEL三種方法開(kāi)展了大量的工作，并得到了相應(yīng)的結(jié)論。但是，這些研究基本基于某一種建模方法開(kāi)展，極少針對(duì)不同建模方法模擬精度進(jìn)行對(duì)比分析。數(shù)值模擬作為彌補(bǔ)試驗(yàn)的重要手段，必須合理選擇建模方法，因此需要了解常用建模方法的模擬精度及其差別。本文參考Varas等[5]試驗(yàn)結(jié)果，基于ALE，CEL和SPH三種方法，建立相應(yīng)的動(dòng)力學(xué)數(shù)值分析模型，對(duì)比沖擊過(guò)程中流體壓力、空腔尺寸、破片速度衰減、充液結(jié)構(gòu)變形以及計(jì)算成本等，研究三種方法的模擬精度及其差別，為破片高速?zèng)_擊充液結(jié)構(gòu)建模方法選取提供基礎(chǔ)。

2 建模方法

2.1 ALE建模方法

ALE方法與拉格朗日網(wǎng)格重新分區(qū)法緊密相關(guān)，保證網(wǎng)格的運(yùn)動(dòng)獨(dú)立于物質(zhì)，避免網(wǎng)格畸變。該方法在每個(gè)時(shí)間步先進(jìn)行拉格朗日運(yùn)動(dòng)計(jì)算，隨后重新劃分網(wǎng)格進(jìn)行下一步的計(jì)算。采用ALE方法離散的質(zhì)量、動(dòng)量和能量守恒控制方程分別見(jiàn)式(1～3)。

(1)

(2)

(3)

2.2 CEL建模方法

CEL方法一般使用拉格朗日法離散固體，使用歐拉法離散流體，兩者的計(jì)算一開(kāi)始是完全分開(kāi)的，但在界面處拉格朗日單元和歐拉單元互為邊界條件，并采用圖1的方法進(jìn)行耦合計(jì)算，保證滿足界面協(xié)調(diào)條件，解決復(fù)雜的流-固耦合問(wèn)題和大變形問(wèn)題。

圖1 CEL方法耦合界面信息傳遞[24]

2.3 SPH建模方法

SPH方法作為一種典型的拉格朗日形式的無(wú)網(wǎng)格方法，其用一系列包含材料性質(zhì)和邊界條件等信息的粒子描述系統(tǒng)運(yùn)動(dòng)狀態(tài)。SPH方法最早用于三維開(kāi)放空間的天體物理學(xué)研究，后廣泛應(yīng)用于大變形的動(dòng)力學(xué)問(wèn)題。采用SPH方法離散的質(zhì)量、動(dòng)量和能量守恒控制方程見(jiàn)式(4～6)。

(4)

(5)

(6)

3 有限元數(shù)值建模

為驗(yàn)證模型的有效性，建立與參考文獻(xiàn)[5]中試驗(yàn)?zāi)Ｐ屯耆恢碌挠邢拊Ｐ?。試?yàn)?zāi)Ｐ腿鐖D2所示，其中充液比例100%，箱體尺寸為750 mm×150 mm×150 mm，采用2.5 mm厚6063-T5鋁合金平板，球形破片采用4340鋼彈，直徑12.5 mm，沖擊速度為900 m/s?？紤]模型對(duì)稱性，為提高計(jì)算效率，建立1/2模型，計(jì)算時(shí)長(zhǎng)0.6 ms。

圖2 Varas試驗(yàn)件

為保證三種建模方法得到的有限元模型網(wǎng)格的一致性，本文基于CATIA軟件建立三維模型，并將其導(dǎo)入Hypermesh軟件得到網(wǎng)格單元，分別導(dǎo)入ABAQUS有限元軟件(CEL方法)和 LS -DYNA 有限元軟件(ALE和SPH方法)進(jìn)行接觸、邊界和載荷等前處理。

針對(duì)本文模型，Artero-Guerrero等[12]發(fā)現(xiàn)當(dāng)網(wǎng)格尺寸為2 mm時(shí)，能在合理的計(jì)算時(shí)間成本內(nèi)獲得較高的計(jì)算精度。同時(shí)學(xué)者們發(fā)現(xiàn)在射彈沖擊充液箱體過(guò)程中，當(dāng)射彈、沖擊點(diǎn)附近箱體結(jié)構(gòu)及液體網(wǎng)格尺寸約為1∶1∶1時(shí)，能較好地避免發(fā)生計(jì)算泄漏問(wèn)題(Leakage problems)[10,11]。為提高計(jì)算效率，本文將沖擊主要影響區(qū)域網(wǎng)格加密為 2 mm，且射彈、沖擊點(diǎn)附近箱體結(jié)構(gòu)及液體網(wǎng)格尺寸比例為1∶1∶1。

最終充液結(jié)構(gòu)劃分為23034個(gè)殼單元，球形破片劃分為752個(gè)實(shí)體單元，水劃分為214668個(gè)單元。

3.1 ALE模型

ALE模型中，為有效模擬水的流動(dòng)及水與結(jié)構(gòu)之間的相互作用，必須在四周建立空氣域，其中空氣域的大小必須保證箱體結(jié)構(gòu)變形后仍然在Euler網(wǎng)格區(qū)域內(nèi)[12]。

水及周圍空氣采用歐拉單元，并進(jìn)行共節(jié)點(diǎn)(Ele Tol_DupNod)處理，其中空氣劃分為298506個(gè)單元，整個(gè)模型共536960個(gè)單元。

通過(guò)關(guān)鍵字*ALE_MULTI-MATERIAL_GROUP將水和空氣兩種ALE單元進(jìn)行耦合。流固耦合采用關(guān)鍵詞*CONSTRAINED_LAGRANGE _IN_SOLID設(shè)置接觸，并通過(guò)*MAT_009-NULL和*EOS狀態(tài)方程描述水的運(yùn)動(dòng)，通過(guò)*MAT_009-NUL和*EOS狀態(tài)方程描述空氣的運(yùn)動(dòng)。模型如圖3(a)所示。

3.2 CEL模型

CEL模型中，同樣需要建立空氣域?？諝庥虻拇笮⊥珹LE模型一致。水及周圍空氣采用歐拉單元，并進(jìn)行共節(jié)點(diǎn)處理。各類單元數(shù)量與ALE模型完全一致。

通過(guò)對(duì)水及空氣相應(yīng)模型進(jìn)行Part_Edit_Type_Eulerian處理，得到歐拉單元，并基于Load_Create Predefined Field_Other_Material assignment設(shè)置歐拉單元中水區(qū)域和空氣區(qū)域。流固耦合接觸采用軟件自帶的通用接觸(General contact)即可。通過(guò)EOS狀態(tài)方程和Dynamic Viscosity黏性系數(shù)描述水和空氣的運(yùn)動(dòng)。模型如圖3(b)所示。

3.3 SPH模型

在SPH模型中，水的運(yùn)動(dòng)不需要空氣模型。因此基于Hypermesh軟件得到的網(wǎng)格單元，通過(guò)Mesh_SPH Generation操作，在拉格朗日單元中心位置生成相應(yīng)的SPH單元。

通過(guò)關(guān)鍵詞*BOUNDARY_SPH_SYMMETRY_PLANE設(shè)置SPH單元的對(duì)稱面，并基于關(guān)鍵詞*CONTACT_AUTOMATIC_NODES_TO_SURFACE定義流固耦合接觸。通過(guò)*MAT_009-NULL和*EOS狀態(tài)方程描述水的運(yùn)動(dòng)過(guò)程。模型如圖3(c)所示。

圖3 三種有限元模型

3.4 材料參數(shù)

三種模型均采用Johnson-Cook本構(gòu)模型和失效模型模擬充液結(jié)構(gòu)的破壞變形，球形破片采用彈塑性模型，具體材料參數(shù)列入表1。水和空氣具體材料參數(shù)列入表2。

表2 水及空氣材料參數(shù)[5]

4 模型驗(yàn)證

文獻(xiàn)[5]研究了充液結(jié)構(gòu)特定位置(圖4中點(diǎn)PTn和點(diǎn)PTf)的壓力-時(shí)間歷程和球形破片沖擊過(guò)程中形成的空腔尺寸。

圖4 Varas試驗(yàn)件簡(jiǎn)化模型

當(dāng)球形破片沖擊速度為900 m/s時(shí)，基于ALE，CEL和SPH三種方法，建立的有限元模型得到的點(diǎn)PTn和點(diǎn)PTf的壓力-時(shí)間歷程如圖5所示。圖6分別給出了0.028 ms,0.084 ms和 0.140 ms 時(shí)刻試驗(yàn)獲得的空腔尺寸及數(shù)值仿真得到的空腔尺寸。結(jié)合圖5和圖6可以看出，基于ALE，CEL和SPH方法數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果吻合度較好，誤差基本在15%以內(nèi)，證明建立的有限元模型的準(zhǔn)確性。

圖5 壓力-時(shí)間曲線

圖6 0.028 ms,0.084 ms和0.140 ms時(shí)空腔尺寸(單位：mm)

5 三種方法對(duì)比分析

5.1 充液結(jié)構(gòu)特定位置壓力

球形破片入水后產(chǎn)生的沖擊波壓力場(chǎng)在破片運(yùn)動(dòng)前方以半球狀向四周擴(kuò)張，所經(jīng)之處產(chǎn)生沖擊壓力載荷峰值，之后流體快速膨脹形成空腔，沖擊波后壓力急劇衰減。以點(diǎn)PTn和點(diǎn)PTf作為特定位置點(diǎn)，三種方法得到的壓力-時(shí)間歷程曲線如圖5所示。其中點(diǎn)PTn位于網(wǎng)格加密區(qū)，其附近網(wǎng)格尺寸小于點(diǎn)PTf附近網(wǎng)格尺寸。表3給出了相應(yīng)的壓力載荷峰值及其對(duì)應(yīng)時(shí)間節(jié)點(diǎn)。

表3 壓力峰值及其對(duì)比

結(jié)合圖5和表3可以看出，三種方法均能有效模擬充液結(jié)構(gòu)內(nèi)部特定位置處的壓力載荷及其對(duì)應(yīng)時(shí)間；其中在相同的單元密度下，CEL方法計(jì)算的壓力載荷峰值誤差整體較小，ALE方法次之，SPH方法最大，且SPH方法得到的壓力載荷峰值均偏?。坏獵EL方法在壓力峰值后的計(jì)算結(jié)果震蕩程度最大(曲線后半段)；對(duì)比時(shí)間節(jié)點(diǎn)，發(fā)現(xiàn)SPH方法獲得的壓力載荷峰值時(shí)間點(diǎn)較晚，說(shuō)明SPH方法通過(guò)粒子之間相互作用傳遞力導(dǎo)致一定延遲；對(duì)比點(diǎn)PTn和點(diǎn)PTf載荷峰值平均誤差，可以看到點(diǎn)PTn平均誤差遠(yuǎn)小于點(diǎn)PTf，說(shuō)明單元尺寸越小，壓力載荷峰值計(jì)算精度越高，但對(duì)載荷峰值時(shí)間節(jié)點(diǎn)基本沒(méi)有影響；其中SPH方法計(jì)算精度隨單元尺寸變化影響最大，因此采用SPH方法時(shí)應(yīng)謹(jǐn)慎考慮單元尺寸的選取。

5.2 空腔尺寸

球形破片高速?zèng)_擊充液結(jié)構(gòu)產(chǎn)生水錘效應(yīng)，包括空腔的擴(kuò)張和坍塌階段。選取0.028 ms，0.084 ms 和0.140 ms三個(gè)時(shí)刻，得到相應(yīng)位置的空腔尺寸如圖6所示。表4給出了三種方法數(shù)值分析結(jié)果與試驗(yàn)結(jié)果的誤差值。

表4 空腔尺寸及其對(duì)比

通過(guò)對(duì)比可以看出，在相同單元密度下，CEL方法建立的有限元模型計(jì)算誤差整體而言最大，且得到的空腔尺寸均大于試驗(yàn)結(jié)果；ALE和SPH方法得到的空腔尺寸整體而言與試驗(yàn)結(jié)果誤差較小，其中ALE方法得到的空腔尺寸均大于試驗(yàn)結(jié)果，SPH方法得到的空腔尺寸均小于試驗(yàn)結(jié)果。因此，研究破片高速入射形成的空腔特性或者水錘效應(yīng)的空腔的擴(kuò)張和坍塌階段，ALE和SPH方法建立的有限元模型分析精度更高。

5.3 破片速度變化

沖擊過(guò)程中，球形破片動(dòng)能通過(guò)摩擦阻力和壓差阻力分別傳遞給箱體壁板和水。忽略重力影響，基于牛頓第二定律，球形破片在無(wú)限大水中運(yùn)動(dòng)方程[25]為

(7)

將式(7)進(jìn)行積分，得到破片速度衰減與時(shí)間之間關(guān)系為

vp=1/(βt+1/vo)

(8)

基于三種建模方法得到的球形破片速度衰減-時(shí)間變化曲線如圖7所示，同時(shí)給出了球形破片理論速度衰減-時(shí)間變化曲線?？梢钥闯?，三種方法得到的球形破片速度-時(shí)間曲線變化形式基本一致，且和理論值擬合程度較好；ALE方法在整個(gè)沖擊階段得到的速度曲線比較平滑，但整體而言速度偏小，誤差在10%以內(nèi)；CEL方法在破片初始沖擊時(shí)刻速度變化波動(dòng)較大，在沖擊后半段與理論值擬合最好；SPH方法在破片初始沖擊時(shí)刻擬合效果最好。因此，對(duì)于破片速度的模型，三種方法的精度基本一致，均能滿足工程需求。

圖7 破片速度-時(shí)間曲線

5.4 結(jié)構(gòu)變形

球形破片沖擊過(guò)程中，破片和水錘效應(yīng)共同作用，引起充液結(jié)構(gòu)大面積變形/破片。因此，結(jié)構(gòu)變形/損傷預(yù)測(cè)精度是充液結(jié)構(gòu)毀傷程度的重要指標(biāo)。

在本文模型中，計(jì)算時(shí)間足夠長(zhǎng)，保證獲得的結(jié)構(gòu)變形為最終塑性變形。圖8給出了ALE方法得到的前壁板和后壁板最大位移節(jié)點(diǎn)的位移時(shí)間曲線。可以看出，在3 ms時(shí)，壁板變形已經(jīng)趨于 穩(wěn)定。

圖8 節(jié)點(diǎn)位移

圖9分別給出了充液結(jié)構(gòu)前壁板、后壁板和側(cè)壁板的變形情況。表5給出了試驗(yàn)與三種數(shù)值分析方法得到的壁板最大變形及其誤差。基于圖9和表5可以看出，ALE，CEL和SPH三種方法預(yù)測(cè)的位移曲線變化趨勢(shì)基本一致，與試驗(yàn)結(jié)果擬合較好，均能有效的預(yù)測(cè)充液結(jié)構(gòu)的變形程度；同時(shí)三種方法對(duì)前壁板的預(yù)測(cè)精度遠(yuǎn)高于后壁板，最大誤差為3.8%。這是因?yàn)楹蟊诎逶谄破謴厍?，沖擊波已經(jīng)導(dǎo)致其產(chǎn)生一定的預(yù)應(yīng)力場(chǎng)，故破片沖擊后壁板的物理過(guò)程遠(yuǎn)遠(yuǎn)比前壁板復(fù)雜，導(dǎo)致數(shù)值模擬精度下降。整體而言，CEL方法對(duì)充液結(jié)構(gòu)前壁板和后壁板最大變形的預(yù)測(cè)精度最高。因此，在表征充液結(jié)構(gòu)破片高速?zèng)_擊下的變形模式研究，CEL方法預(yù)測(cè)精度最高，SPH方法次之。

圖9 結(jié)構(gòu)變形曲線

表5 結(jié)構(gòu)最大變形及其對(duì)比

5.5 計(jì)算成本

在流固耦合動(dòng)力學(xué)模型中，時(shí)間成本是必須考慮的問(wèn)題。本文采用DELL Precision 7920工作站開(kāi)展相關(guān)計(jì)算。處理器為Inter(R) Xeon(R) Gold 6148 CPU @2.40GHz(2處理器)，安裝內(nèi)存64 GB。

在相同的單元尺寸下，均提交24核并行計(jì)算。計(jì)算0.5 ms時(shí)，ALE模型計(jì)算總時(shí)長(zhǎng)為2286 s，CEL模型計(jì)算總時(shí)長(zhǎng)為2488 s，SPH模型計(jì)算總時(shí)長(zhǎng)為4324 s?？梢钥闯?，ALE方法和CEL方法計(jì)算時(shí)長(zhǎng)基本一致；同時(shí)考慮到ALE模型和CEL模型中存在空氣模型，總體單元數(shù)量是SPH模型的兩倍多，進(jìn)一步說(shuō)明SPH方法計(jì)算效率低，但SPH方法建模方式相對(duì)簡(jiǎn)便，且后處理相對(duì)簡(jiǎn)單。因此在建模方法選取上，應(yīng)綜合考慮前后處理和計(jì)算時(shí)長(zhǎng)問(wèn)題。

6 結(jié) 論

基于Varas試驗(yàn)結(jié)果和模型，建立相應(yīng)的ALE，CEL和SPH有限元模型，分別對(duì)比研究三種方法對(duì)流體壓力變化、形成的空腔尺寸、破片速度衰減變化和充液結(jié)構(gòu)變形等預(yù)測(cè)精度，并分析相應(yīng)的計(jì)算時(shí)間成本，得到如下結(jié)論。

(1) ALE，CEL和SPH三種方法均能有效模擬破片高速?zèng)_擊充液結(jié)構(gòu)的流固耦合動(dòng)力學(xué)過(guò)程，且數(shù)值分析結(jié)果與試驗(yàn)結(jié)果吻合較好，誤差在可接受范圍內(nèi)。

(2) ALE方法預(yù)測(cè)的空腔尺寸精度較高；CEL方法預(yù)測(cè)的流體壓力、破片速度衰減和充液結(jié)構(gòu)變形精度較高；SPH方法預(yù)測(cè)的空腔尺寸和破片速度衰減精度較高。

(3) 當(dāng)網(wǎng)格尺寸一致時(shí)，SPH方法計(jì)算時(shí)長(zhǎng)約為ALE和CEL方法的兩倍，但SPH方法前后處理更加簡(jiǎn)便，且該算法不會(huì)出現(xiàn)負(fù)體積等計(jì)算問(wèn)題。

因此，針對(duì)具體的研究問(wèn)題，應(yīng)充分考慮研究對(duì)象、前后處理和計(jì)算時(shí)長(zhǎng)問(wèn)題，選擇合適的建模方法。