基于ALE，CEL和SPH方法的球形破片高速沖擊充液結構對比研究

張 宇， 王彬文， 劉小川

(1.中國飛機強度研究所，西安 710065；2.結構沖擊動力學航空科技重點實驗室，西安 710065)

1 引 言

易損性研究是提高飛機作戰生存力的基礎研究之一。其中油箱是飛機最大的易損結構，約占飛機總體易損面積的75%[1]。飛機油箱毀傷包括水錘效應和燃爆效應，前者通過破壞油箱結構完整性，后者通過引燃/引爆航空煤油，導致飛機失去作戰能力。

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

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

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

2 建模方法

2.1 ALE建模方法

ALE方法與拉格朗日網格重新分區法緊密相關，保證網格的運動獨立于物質，避免網格畸變。該方法在每個時間步先進行拉格朗日運動計算，隨后重新劃分網格進行下一步的計算。采用ALE方法離散的質量、動量和能量守恒控制方程分別見式(1～3)。

(1)

(2)

(3)

2.2 CEL建模方法

CEL方法一般使用拉格朗日法離散固體，使用歐拉法離散流體，兩者的計算一開始是完全分開的，但在界面處拉格朗日單元和歐拉單元互為邊界條件，并采用圖1的方法進行耦合計算，保證滿足界面協調條件，解決復雜的流-固耦合問題和大變形問題。

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

2.3 SPH建模方法

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

(4)

(5)

(6)

3 有限元數值建模

為驗證模型的有效性，建立與參考文獻[5]中試驗模型完全一致的有限元模型。試驗模型如圖2所示，其中充液比例100%，箱體尺寸為750 mm×150 mm×150 mm，采用2.5 mm厚6063-T5鋁合金平板，球形破片采用4340鋼彈，直徑12.5 mm，沖擊速度為900 m/s?？紤]模型對稱性，為提高計算效率，建立1/2模型，計算時長0.6 ms。

圖2 Varas試驗件

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

針對本文模型，Artero-Guerrero等[12]發現當網格尺寸為2 mm時，能在合理的計算時間成本內獲得較高的計算精度。同時學者們發現在射彈沖擊充液箱體過程中，當射彈、沖擊點附近箱體結構及液體網格尺寸約為1∶1∶1時，能較好地避免發生計算泄漏問題(Leakage problems)[10,11]。為提高計算效率，本文將沖擊主要影響區域網格加密為 2 mm，且射彈、沖擊點附近箱體結構及液體網格尺寸比例為1∶1∶1。

最終充液結構劃分為23034個殼單元，球形破片劃分為752個實體單元，水劃分為214668個單元。

3.1 ALE模型

ALE模型中，為有效模擬水的流動及水與結構之間的相互作用，必須在四周建立空氣域，其中空氣域的大小必須保證箱體結構變形后仍然在Euler網格區域內[12]。

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

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

3.2 CEL模型

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

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

3.3 SPH模型

在SPH模型中，水的運動不需要空氣模型。因此基于Hypermesh軟件得到的網格單元，通過Mesh_SPH Generation操作，在拉格朗日單元中心位置生成相應的SPH單元。

通過關鍵詞*BOUNDARY_SPH_SYMMETRY_PLANE設置SPH單元的對稱面，并基于關鍵詞*CONTACT_AUTOMATIC_NODES_TO_SURFACE定義流固耦合接觸。通過*MAT_009-NULL和*EOS狀態方程描述水的運動過程。模型如圖3(c)所示。

圖3 三種有限元模型

3.4 材料參數

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

表2 水及空氣材料參數[5]

4 模型驗證

文獻[5]研究了充液結構特定位置(圖4中點PTn和點PTf)的壓力-時間歷程和球形破片沖擊過程中形成的空腔尺寸。

圖4 Varas試驗件簡化模型

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

圖5 壓力-時間曲線

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

5 三種方法對比分析

5.1 充液結構特定位置壓力

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

表3 壓力峰值及其對比

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

5.2 空腔尺寸

球形破片高速沖擊充液結構產生水錘效應，包括空腔的擴張和坍塌階段。選取0.028 ms，0.084 ms 和0.140 ms三個時刻，得到相應位置的空腔尺寸如圖6所示。表4給出了三種方法數值分析結果與試驗結果的誤差值。

表4 空腔尺寸及其對比

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

5.3 破片速度變化

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

(7)

將式(7)進行積分，得到破片速度衰減與時間之間關系為

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

(8)

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

圖7 破片速度-時間曲線

5.4 結構變形

球形破片沖擊過程中，破片和水錘效應共同作用，引起充液結構大面積變形/破片。因此，結構變形/損傷預測精度是充液結構毀傷程度的重要指標。

在本文模型中，計算時間足夠長，保證獲得的結構變形為最終塑性變形。圖8給出了ALE方法得到的前壁板和后壁板最大位移節點的位移時間曲線?？梢钥闯?，在3 ms時，壁板變形已經趨于 穩定。

圖8 節點位移

圖9分別給出了充液結構前壁板、后壁板和側壁板的變形情況。表5給出了試驗與三種數值分析方法得到的壁板最大變形及其誤差?；趫D9和表5可以看出，ALE，CEL和SPH三種方法預測的位移曲線變化趨勢基本一致，與試驗結果擬合較好，均能有效的預測充液結構的變形程度；同時三種方法對前壁板的預測精度遠高于后壁板，最大誤差為3.8%。這是因為后壁板在破片侵徹前，沖擊波已經導致其產生一定的預應力場，故破片沖擊后壁板的物理過程遠遠比前壁板復雜，導致數值模擬精度下降。整體而言，CEL方法對充液結構前壁板和后壁板最大變形的預測精度最高。因此，在表征充液結構破片高速沖擊下的變形模式研究，CEL方法預測精度最高，SPH方法次之。

圖9 結構變形曲線

表5 結構最大變形及其對比

5.5 計算成本

在流固耦合動力學模型中，時間成本是必須考慮的問題。本文采用DELL Precision 7920工作站開展相關計算。處理器為Inter(R) Xeon(R) Gold 6148 CPU @2.40GHz(2處理器)，安裝內存64 GB。

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

6 結 論

基于Varas試驗結果和模型，建立相應的ALE，CEL和SPH有限元模型，分別對比研究三種方法對流體壓力變化、形成的空腔尺寸、破片速度衰減變化和充液結構變形等預測精度，并分析相應的計算時間成本，得到如下結論。

(1) ALE，CEL和SPH三種方法均能有效模擬破片高速沖擊充液結構的流固耦合動力學過程，且數值分析結果與試驗結果吻合較好，誤差在可接受范圍內。

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

(3) 當網格尺寸一致時，SPH方法計算時長約為ALE和CEL方法的兩倍，但SPH方法前后處理更加簡便，且該算法不會出現負體積等計算問題。

因此，針對具體的研究問題，應充分考慮研究對象、前后處理和計算時長問題，選擇合適的建模方法。