一種自適應軸對稱FEM-SPH 耦合算法及其在高速沖擊模擬中的應用*

肖毅華，胡德安，韓 旭，楊 剛

(湖南大學汽車車身先進設計制造國家重點實驗室，湖南 長沙410082)

高速沖擊涉及材料的大變形、破碎和飛濺等現象，應用基于網格的方法對其進行數值模擬存在困難。拉格朗日網格方法會遭遇單元畸變而計算終止;單元侵蝕技術可以克服畸變問題，但會導致能量損失并且改變物質的幾何邊界。歐拉網格方法難于精確地捕捉物質界面和跟蹤材料的歷史數據。光滑粒子流體動力學(SPH)［1-2］作為一種無網格、拉格朗日粒子法，能克服基于網格的方法的缺陷。SPH 在處理大變形方面較有限元法(FEM)等拉格朗日網格方法有優勢，但在模擬小變形時的計算精度和效率都不及FEM，并且SPH 的邊界處理不如FEM 方便。對于高速沖擊問題，由于大變形通常只發生在有限的局部區域，耦合FEM 和SPH 是一種十分有效的模擬途徑。在絕大部分的小變形區域內用FEM 計算，在局部的大變形區域內用SPH 計算，這樣既能保證計算精度又能提高計算效率。目前，一些研究者［3-10］對FEM-SPH 耦合算法及其在高速沖擊問題中的應用進行了研究，并取得了一些有意義的結果。

本文中，提出一種自適應軸對稱FEM-SPH 耦合算法，并用于高速沖擊模擬。該算法自動將畸變單元轉化為粒子，既能避免單元畸變，又盡可能地減小SPH 計算區域。文獻［4-6，9］也給出了類似的算法，本文算法與它們的主要區別在于:提出一種新的耦合算法，實現單元和粒子計算的高精度耦合;提出采用最小內角轉化準則和單元分組轉化方式來實現單元向粒子的自動轉化。通過計算應力波傳播、泰勒桿撞擊問題，驗證算法的正確性;并通過模擬彈丸侵徹鋁板和混凝土問題，說明算法在高速沖擊模擬中的可靠性和實用性。

1 基本方程

為了求解軸對稱高速沖擊問題，采用如下形式的軸對稱SPH 方程

式中:ρ 為三維密度，η=2πrρ 為平面密度，m 為質量，v 為速度，σ 為應力，e 為質量內能，r 和z 分別為徑向和軸向坐標，W 為光滑函數。光滑函數選擇為B 樣條函數［11］。初始光滑長度

式中:αh為光滑長度放大因子，Δ 為粒子間距。光滑長度在計算中保持不變。為了穩定計算、消除數值振蕩，SPH 法采用文獻［12］中的結合型人工粘性。關于求解沖擊動力問題的軸對稱FEM 的基本方程，參考文獻［13］。

2 FEM 與SPH 的自適應耦合

為了實現FEM 和SPH 的自適應耦合，需要處理3 個關鍵問題:(1)同一物體中的單元與粒子的耦合;(2)單元與單元、單元與粒子以及粒子與粒子的接觸;(3)單元向粒子的自動轉化。下面介紹對上述3 個問題的處理方法。

2.1 單元-粒子耦合界面處理

本文中提出一種新的高精度耦合算法，圖1 為該算法的示意圖。耦合算法的關鍵是計算耦合界面附近單元與粒子間的相互作用力。本文算法計算相互作用力的基本思想是:一方面，將耦合界面附近的單元作為虛粒子包含到SPH 方程中，得到鄰近單元對粒子的作用力;另一方面，根據耦合界面附近的粒子應力確定作用在耦合界面上的單元邊上的等效面力，再將這些面力等效到單元節點上，最終得到鄰近粒子對單元節點的作用力。

耦合界面附近被粒子的影響域覆蓋的單元(以下稱耦合界面單元)作為虛粒子處理。每一個耦合界面單元作為一個虛粒子。圖2為虛粒子生成示意圖。虛粒子的變量在每一個時間步都根據所在單元的相關變量重新確定，質量和密度與對應單元的質量和密度相等，坐標和速度取為單元3 個節點的相應值的算術平均值，應力與單元的應力相等。虛粒子的半徑和光滑長度按下述方式確定

圖1 單元-粒子耦合算法示意圖Fig.1 Description of the element-particle coupling algorithm

圖2 虛粒子生成示意圖Fig.2 Description of the generation of imaginary particles

式中:rIP，i和hIP，i分別為虛粒子半徑和光滑長度，AE，I為虛粒子所在單元的面積，αh的含義同式(4)中的αh。

計算鄰近單元對粒子的作用力只需將虛粒子包含到SPH 計算中。如圖1 所示，粒子i 的影響域內包含實粒子(n1，…，n12)和虛粒子(n13，…，n19)。計算粒子i 時，實粒子和虛粒子都參與SPH 方程(1)～(3)的求和，SPH 方程的形式不發生改變。這樣，在動量方程(2a)和(2b)中，虛粒子的貢獻即為鄰近單元對粒子的作用力。另外，虛粒子也應用于粒子的應變率的計算。這可以克服SPH 近似的邊界缺陷，提高近似精度，使應變率和應力計算更準確。

圖3 為計算鄰近粒子對單元節點的作用力的示意圖。對于每一條位于耦合界面的單元邊，根據與它關聯的粒子的應力確定一個等效面力，并認為其均勻作用在單元邊上。等效面力的計算方式為

圖3 鄰近粒子對單元節點的作用力的計算Fig.3 Calculation of the forces on element nodes from adjacent particles

2.2 接觸界面處理

不同物體或者同一物體的不同部分間存在接觸作用，需要采用有效的算法進行處理，防止出現非物理的穿透。目前，已有一些較好的算法來處理各種類型的接觸，本文中基于已有接觸算法來進行接觸界面計算。對于單元與單元間的接觸，采用G.R.Johnson 等［14］提出的對稱接觸和滑移界面算法處理。對于單元與粒子間的接觸，參考G.R.Johnson 等［15］處理GPA 節點與單元接觸的算法處理。對于粒子與粒子間的接觸，參照J.C.Campbell 等［16］提出的粒子-粒子型罰接觸算法處理。

2.3 單元向粒子的自動轉化

為實現單元自動轉化為粒子，首先必須確定需要轉化為粒子的畸變單元(以下簡稱轉化單元)。本文中中采用三角形單元計算，其最小內角能很好地反映單元的畸變程度。因此，我們提出采用最小內角轉化準則來確定轉化單元，即:如果某一單元的最小內角小于指定的臨界值θc時，則該單元為轉化單元;否則該單元無需轉化。本文中，取θc=30°。

在確定轉化單元時，除了采用最小內角轉化準則，還結合了一種分組轉化方式。圖4 為單元分組轉化的示意圖。每個物體內的單元在計算的初始化階段分成單元組。在計算過程中，如果一個單元組中存在一個單元滿足最小內角轉化準則，則該組中的所有單元均作為轉化單元。使用分組轉化的好處是:快速地搜索到耦合界面單元，提高計算效率;同時，使耦合界面的幾何更連續和規則，提高計算精度和穩定性。單元分組的過程很簡單。首先找到包含物體的最小矩形域;然后，將矩形域劃分為參考尺寸為Δs的矩形子域。Δs的大小為

圖4 單元分組轉化Fig.4 Group-based conversion of elements

式中:κ 為大于1 的因數，它使單元組轉化生成的粒子的影響域不會超出其鄰接單元組覆蓋的范圍，本文計算取κ=1.5;Rmax為最大可能的粒子影響域半徑，按下式估計

式中:A0，max為初始時刻單元的最大面積。根據確定的Δs，計算沿徑向和軸向的矩形子域數目分別為

然后，可以計算矩形子域在徑向和軸向的精確尺寸為

將矩形子域編號為(Ir，Iz)，Ir和Iz分別表示矩形子域沿徑向和軸向的序號，取值分別為1 ～Nr和1～Nz。根據此編號方式，容易確定一個單元所在的矩形子域的編號為

式中:(rg，zg)為單元的重心坐標。由式(13)確定每個單元所在的矩形子域，同一矩形子域中的單元為一個單元組。

采用上述方式確定轉化單元后，將這些單元從單元列表中刪除，并且增加新的粒子。新增粒子的變量的確定方式與前述虛粒子相同。與虛粒子不同的是，這些粒子產生后，在后續計算中其變量通過SPH 方程計算。

轉化單元變成粒子后，單元區域與粒子區域的耦合界面發生變化，需更新耦合界面單元以用于耦合計算。圖4 說明了使用分組轉化時如何快速地搜索到耦合界面單元。例如，單元組(I，J)中的單元A 滿足最小內角轉化準則，則單元組(I，J)中的所有單元轉化為粒子，而單元組(I，J)的鄰接單元組中的單元作為耦合界面單元。根據單元組的編號規則，單元組(I，J)的鄰接單元組容易確定為

轉化單元變成粒子后，單元區域表面也需更新。圖5 為圖4 中單元A 所在的單元組。該單元組表面包含位于單元區域表面的線段N1N2，N2N3，…，N5N6和內部的線段N6N7，N7N8，…，N24N1。當該單元組內的單元轉化為粒子后，線段N1N2，N2N3，…，N5N6從單元區域表面中刪除，而線段N6N7，N7N8，…N24N1增加到單元區域表面中去。單元區域表面更新后，單元邊與粒子間的關聯也相應更新。與位于線段N1N2，N2N3，…，N5N6上的單元邊關聯的粒子被釋放;由連接到線段N6N7，N7N8，…，N24N1的單元轉化生成的粒子需要關聯到該線段上的單元邊。確定新生成的粒子與單元邊的關聯時，存在圖5 中類似單元A 和B 的2 種典型情況。單元A 只有1 條邊(N6N7)位于新增的單元區域表面上，由單元A 生成的粒子與邊N6N7關聯，并且與該邊固連，不允許與其發生分離或相對滑移。單元B 有2 條邊(N8N9和N9N10)位于新增的單元區域表面上，由單元B 生成的粒子同時與邊N8N9和N9N10關聯，但只與兩者中距離最近的邊固連。為實現粒子與單元邊的固連，本文中應用G.R.Johnson 等［15］連接GPA 節點與單元邊的算法。

圖5 單元區域表面和單元邊與粒子間關聯的更新Fig.5 Update of surfaces of element regions and association between element sides and particles

3 數值算例

3.1 彈性桿中的應力波傳播

如圖6(a)所示，半徑1 mm、長20 mm 的彈性桿以50 m/s 速度垂直撞擊剛性面。彈性桿的材料參數為:密度ρ=7.85 g/cm3，楊氏模量E=206 GPa，泊松比ν=0.3。由于彈性桿具有大長徑比，本問題近似為一維應力問題，可以得到解析解，為驗證算法提供參考。采用耦合算法和有限元法計算此問題，并將計算結果與解析解進行對比。圖6(b)為耦合計算的模型，彈性桿的上、下兩半分別采用1 000 個粒子和1 000 個單元離散。有限元法采用2 000 個單元離散彈性桿，單元的分布形式與耦合計算模型相同。有限元求解采用LS-DYNA 軟件。

圖7 為計算得到的3 μs 時沿軸線r=0.5 mm 的應力分布。本文耦合算法的計算結果與DYNA 計算結果符合很好，耦合界面處的應力很光滑，應力波準確地傳播通過耦合界面。耦合算法計算的應力波幅值和前沿位置與一維解析解也基本相符，只是由于采用二維軸對稱計算，存在彌散效應，所以應力有一定振蕩。由上述計算結果可見，本文提出的耦合算法是比較精確的。

圖6 彈性桿撞擊剛性面的模型Fig.6 Model of elastic bar impacting a rigid wall

圖7 3 μs 時沿軸線r=0.5 mm 的應力分布Fig.7 Stress distributions along the axis of r=0.5mm at 3 μs

3.2 泰勒桿

泰勒桿實驗是采用金屬圓桿高速正撞擊剛性平面，經常被用來研究材料的本構模型和驗證數值算法的準確性。采用本文的自適應耦合算法對泰勒桿實驗進行模擬，將模擬結果與實驗結果對比，驗證本文算法的準確性。本文中模擬的泰勒桿實驗為:直徑D0=7.6 mm、長度L0=25.4 mm 的阿姆克鐵圓桿以221 m/s 的速度正撞擊剛性面。模擬中，將阿姆克鐵視為理想彈塑性材料，使用Mie-Grüneisen 狀態方程計算材料壓力，相關的材料特性參數取自文獻［17］。圓桿在初始時刻采用2 680 個單元離散。

圖8 圓桿在不同時刻的變形Fig.8 Deformation of the cylinder at different time

圖8 為自適應耦合算法計算的圓桿在10、20、40 和60 μs 等4 個時刻的變形圖。由該圖可見，圓桿下端的單元隨著變形的發展逐漸轉化為粒子;轉化形成的粒子區域具有很好的連續性，粒子區域和單元區域間的耦合界面較規則。計算在60 μs 時終止，此時圓桿幾乎不再變形。在60 μs 時刻的變形圖中，對比了計算的圓桿變形形狀與實驗結果［18］(黑點表示)，兩者符合較好。計算的圓桿體變形后的特征尺寸為:長度Lc=19.1 mm，撞擊端直徑Dc=13.9 mm，離撞擊端0.2L0處的寬度Wc=8.6 mm;相應的實驗結果［18］依次為Le=19.8 mm，De=13.7 mm，We=8.6 mm;兩者的平均相對誤差(ε=(|Lc-Le|/Le+|Dc-De|/De+|Wc-We|/We)/3)不到2%。計算的圓桿長度略小于實驗結果，其誤差相對于其他2 個特征尺寸的誤差要大，這與文獻［17］給出的計算結果相似，可能需要采用更合理的材料模型來提高模擬結果的準確性。綜合上述計算結果可見，本文的自適應耦合算法和相應的程序是正確可靠的。

3.3 鋁板侵徹

考慮直徑12.9 mm、長88.9 mm、φ(頭部曲率半徑與直徑之比)為3.0 的卵形4340 鋼彈以不同速度正侵徹厚26.3 mm、邊長304 mm 的6061-T651 鋁板。采用自適應耦合算法和SPH 法對此問題進行對比計算。計算中，彈體視為理想彈塑性體;靶體等效為304 mm 直徑的圓靶，采用Johnson-Cook 粘塑性模型和Mie-Grüneisen 狀態方程。彈體的材料參數為:ρ=7.83 g/cm3，E=206 GPa，ν=0.3，Y0=1.43 GPa。靶體的材料參數見表1。靶體的基本參數(ρ、E 和ν)取自文獻［19］，Johnson-Cook 模型參數A、B 和n 通過最小二乘法擬合實驗測得的應力應變曲線［19］得到，其余參數參考文獻［20］中相同材料的參數。自適應耦合計算的初始模型由11 444 個三角形單元組成;SPH 計算模型由10 963 個粒子組成，彈/靶接觸界面采用粒子-粒子接觸算法處理。

表1 靶體材料參數Table 1 Material properties of target

圖9 給出了2 種算法計算的5 種不同初始沖擊速度vi下彈體的剩余速度vr，并與實驗結果［19］進行了對比。自適應耦合算法和SPH 的結果與實驗結果都符合較好。在5 種工況下，自適應耦合法與實驗結果的平均相對誤差為3.4%，SPH 法與實驗結果的平均相對誤差為2.9%，兩者的計算精度接近。

圖10 給出了初始沖擊速度為508 m/s 工況下2 種算法計算的相同時刻的變形結果(只給出了彈體和靶體中心區域的變形)，2 種算法給出的變形很接近。應用自適應耦合算法時，只有彈/靶接觸界面附近的局部區域內的少數單元轉化為粒子。在相同的計算條件下，到200 μs 時，自適應耦合算法耗時20 min，SPH 方法耗時71 min，前者的計算效率大大高于后者。對于其他工況，也有相同的結論。

圖9 彈體剩余速度的比較Fig.9 Comparison of residual velocities of the projectile

圖10 508 m/s 沖擊速度下的變形比較Fig.10 Comparison of deformations at impact velocity of 508 m/s

3.4 混凝土侵徹

考慮長143.7 mm、直徑25.4 mm、φ=3.0 的卵形彈以606 m/s 的速度正侵徹厚178 mm、抗壓強度48 MPa 的混凝土板。采用自適應耦合算法和固定耦合算法進行計算。自適應耦合計算的初始模型采用1 124 個三角形單元離散彈體和33 600 個三角形單元離散靶體，計算過程中將畸變單元轉化為粒子。固定耦合計算模型采用1 124 個三角形單元離散彈體和33 600 個粒子離散靶體，計算中無需將單元轉化為粒子。彈體考慮為彈塑性等向硬化材料，材料參數為:ρ=8.3 g/cm3，E=200 GPa，ν=0.3，Y0=1.72 GPa，切線模量ET=15 GPa。混凝土材料模型采用HJC 模型，材料參數參考文獻［21］。

圖11 為彈體的速度時間歷程。自適應耦合算法與固定耦合算法的計算結果很接近，計算的剩余速度分別為464 和466 m/s，均與實驗結果［22］(449 m/s)符合較好。圖12 為2 種方法計算的4 個時刻的變形。在每一變形圖中，軸線右邊為自適應耦合計算結果，左邊為固定耦合計算結果。從圖中可以看到，彈體在侵徹過程中變形很小，混凝土靶體被侵徹后發生嚴重破壞，在彈體入射面和射出面出現材料的飛濺和脫落。這些現象與物理實驗中觀察到的現象相符。圖13 為計算得到的混凝土靶的損傷圖像。2 種算法給出的損傷的分布趨勢大致相同:在彈體通道的周圍形成連續的完全損傷區域;在遠離彈體通道處，出現損傷裂紋，并擴展到靶板的邊緣。上述結果表明:兩種算法均能較好地模擬彈體對混凝土靶的侵徹毀傷過程。從計算效率上看，自適應耦合算法和固定耦合算法分別耗時59 和167 min，前者的計算效率遠高于后者。

圖11 彈體的速度時間歷程Fig.11 Velocity histories of the projectile

圖12 自適應耦合和固定耦合計算的變形比較Fig.12 Comparison of deformations calculated by adaptive and fixed coupling

圖13 混凝土損傷圖像的比較Fig.13 Comparison of the damage of the concrete plate

4 結 論

提出了一種自適應軸對稱FEM-SPH 耦合算法并開發了相應的程序。通過對應力波傳播和泰勒桿問題的計算驗證了算法和程序的正確性，并將該算法應用于彈體侵徹鋁板和混凝土板的模擬中，獲得了與實驗相符的計算結果。提出的自適應算法既克服了單元畸變，又盡可能地減小了SPH 計算區域。它充分利用了FEM 和SPH 各自的優勢，計算精度和效率都較好，非常適合于高速沖擊問題模擬。另外，本文算法雖是針對二維軸對稱問題的，但其具有很好的擴展性，其中涉及的單元-粒子耦合算法、單元向粒子的轉化算法等都很容易拓展到三維問題。在后續研究中，將基于本文思想，進一步開發三維自適應FEM-SPH 耦合算法，從而滿足更多實際工程問題的模擬需求。

［1］ Lucy L B.A numerical approach to the testing of the fission hypothesis［J］.The Astronomical Journal，1977，82(12):1013-1024.

［2］ Gingold R A，Monaghan J J.Smoothed particle hydrodynamics:Theory and application to non-spherical stars［J］.Monthly Notices of the Royal Astronomical Society，1977，181(2):375-389.

［3］ 武玉玉，何遠航，李金柱.耦合方法在超高速碰撞數值模擬中的應用［J］.高壓物理學報，2005，19(4):385-389.WU Yu-yu，HE Yuan-hang，LI Jin-zhu.Application of the coupling method in simulating the hypervelocity impact［J］.Chinese Journal of High Pressure Physics，2005，19(4):385-389.

［4］ 王吉，王肖鈞，卞梁.光滑粒子法與有限元的耦合算法及其在沖擊動力學中的應用［J］.爆炸與沖擊，2007，27(6):522-528.WANG Ji，WANG Xiao-jun，BIAN Liang.Linking of smoothed particle hydrodynamics method to standard finite element method and its application in impact dynamics［J］.Explosion and Shock Waves，2007，27(6):522-528.

［5］ 張志春，強洪夫，高巍然.一種光滑粒子流體動力學-有限元法轉換算法及其在沖擊動力學中的應用［J］.西安交通大學學報，2011，45(1):105-110.ZHANG Zhi-chun，QIANG Hong-fu，GAO Wei-ran.Conversion of 3D distorted finite elements into SPH particles during impact dynamic deformation［J］.Journal of Xi’an Jiaotong University，2011，45(1):105-110.

［6］ Johnson G R.Linking of Lagrangian particle methods to standard finite element methods for high velocity impact simulations［J］.Nuclear Engineering and Design，1994，150(2/3):265-274.

［7］ Attaway S W，Heinstein M W，Swegle JW.Coupling of Smooth particle hydrodynamic with finite element method［J］.Nuclear Engineering and Design，1994，150(2/3):199-205.

［8］ De Vuyst T，Vignjevic R，Campbell J C.Coupling between meshless and finite element methods［J］.International Journal of Impact Engineering，2005，31(8):1054-1064.

［9］ Sauer M.Adaptive kopplung des netzfreien SPH-verfahrens mit finiten elementen zur berechnung von impaktvorg-?ngen［D］.Munich:Universit?t der Bundeswehr München，2000.

［10］ Xiao Y H，Han X，Hu D A.A coupling algorithm of finite element method and smoothed particle hydrodynamics for impact computations［J］.Computers，Materials＆Continua，2011，23(1):9-34.

［11］ Monaghan J J，Lattanzio J C.A refined particle method for astrophysical problems［J］.Astronomy and Astrophysics，1985，149(1):135-143.

［12］ Johnson G R.Artificial viscosity effects for SPH impact computations［J］.International Journal of Impact Engineering，1996，18(5):477-488.

［13］ Johnson G R，Stryk R A，Holmquist T J，et al.Numerical algorithm in a Lagrangian hydrocode［R］.WL-TR-1997-7039.Wright Laboratory，1997.

［14］ Johnson G R，Stryk R A.Symmetric contact and sliding interface algorithms for intense impulsive loading computations［J］.Computer Methods in Applied Mechanics and Engineering，2001(35/36):4531-4549.

［15］ Johnson G R，Stryk R A，Beissel S R，et al.An algorithm to automatically convert distorted finite elements into meshless particles during dynamic deformation［J］.International Journal of Impact Engineering，2002，27(10):997-1013.

［16］ Campbell J C，Vignjevic R，Libersky L D.A contact algorithm for smoothed particle hydrodynamics［J］.Computer Methods in Applied Mechanics and Engineering，2000，184(1):49-65.

［17］ Libersky L D，Petschek A G，Carney T C，et al.High strain Lagrangian hydrodynamics:A three-dimensional SPH code for dynamic material response［J］.Journal of Computational Physics，1993，109(1):67-75.

［18］ Johnson G R，Holmquist T J.Evaluation of cylinder-impact test data for constitutive model constants［J］.Journal of Applied Physics，1988，64(8):3901-3910.

［19］ Piekutowski A J，Forrestal M J，Poormon KL，et al.Perforation of aluminum plates with ogive-nose steel rods at normal and oblique impacts［J］.International Journal of Impact Engineering，1996，18(6):877-887.

［20］ Kmetyk L N，Yarrington P.CTH analyses of steel rod penetration into aluminum and concrete targets with comparison to experiment data［R］.SAND 94-1498.Sandia National Laboratories，1994.

［21］ Holmquist T J，Johnson G R.A computational constitutive model for concrete subjected to large strains，high strain rates and high pressure［C］∥Murphy M J，Backofen J E.14th International Symposium on Ballistics.USA:American Defense Preparedness Association，1993:591-600.

［22］ Hanchak S J，Forrestal M J，Young E R，et a1.Perforation of concrete slabs with 48 MPa and 140 MPa unconfined compressive strengths［J］.International Journal of Impact Engineering，1992，12(1):1-7.