基于FE-SPH耦合算法的球體沖擊問題研究

李旦， 趙廷渝

(中國民用航空飛行學院飛行技術學院, 四川廣漢618307)

基于FE-SPH耦合算法的球體沖擊問題研究

李旦， 趙廷渝

(中國民用航空飛行學院飛行技術學院, 四川廣漢618307)

為充分利用SPH方法處理大變形問題與FEM方法處理沖擊問題時各自的優勢，提出了新型的三維SPH-FEM耦合算法，該算法在大變形區域采用SPH粒子離散，從而克服FEM方法單元畸變存在的問題，其他模擬區域使用FEM單元離散，既可以避免SPH邊界效應，也能保證耦合界面物理量具有連續性。使用該算法分別對球體沖擊土壤與水域的過程進行了三維數值模擬，選取含損傷的Johnson-Cook模型和Mie-Gruneisen狀態方程對土壤與水域進行參數計算，對比分析球體在土壤與水域中的速度、加速度以及有效應力數據結果，分析可知，此次模擬結果與其他已有仿真結果基本吻合，驗證了采用FE-SPH 耦合算法處理結構與流體相互作用問題是可行和有效的，為求解流固耦合問題提供了一條嶄新的途徑。

FE- SPH耦合算法；邊界條件；流固耦合；數值模擬

引言

諸如飛機水上迫降、返回艙回收以及波浪對海洋結構的作用等情況是工程實際中普遍存在流固耦合問題。由于該耦合效應具有非線性特性，利用傳統解析法進行求解具有一定的局限性，因此，數值模擬方法是解決此問題的重要手段[1]。

FEM方法采用拉格朗日格式求解結構域，ALE格式或者歐拉格式求解流體域，該方法雖然計算效率高，但需要劃分網格，尤其在拓撲結構變化很大時進行網格劃分困難，且后期計算處理較復雜，極有可能導致網格畸變[2]。因此，上述方法在模擬復雜自由表面流動與結構大變形相互作用問題時存在困難。SPH方法起初用于處理三維開放空間天體物理學問題，由于該方法不受網格束縛，在模擬大變形以及自由表面流動等問題時能有效避免網格畸變，且能在拉格朗日框架下對流體和結構同時進行求解，目前已廣泛應用于流體力學與連續固體力學。但是，在固體力學領域，用該方法模擬固體介質存在拉伸不穩定問題，計算精度和效率相比FEM方法較低。由于SPH方法中近似核函數沒有Kronecher delta張量性質，處理邊界條件比較困難[3]，因此，單獨采用SPH方法解決此類問題仍存在不足。

近年來，為解決SPH和FEM耦合界面的計算問題，FE-SPH耦合算法被提出并得到快速發展，其基本思想為：采用SPH方法處理局部大變形區域，FEM計算處理小變形區域，便于綜合利用有限元法計算精度好、效率高、易于處理復雜邊界與SPH方法模擬自由表面流動的優勢[4]。目前，FE-SPH耦合算法在工程領域得到廣泛應用，尤其是在數值模擬沖擊動力學問題方面，例如在處理高速沖擊及爆炸問題上普遍采用該耦合算法，同時，在材料加工處理方面，FE-SPH耦合算法也得到了應用。

1FE-SPH耦合算法基本概念

1.1三維FE-SPH耦合算法基本原理

三維FE-SPH耦合算法基本原理如圖1所示。Johnson G R提出將SPH粒子固結在有限元節點，通過固結的方式保證兩者之間傳遞相互作用力，且SPH粒子與對應單元節點加速度相同，從而實現SPH與FE耦合[5]。但是，在計算界面粒子應變率時該耦合算法沒有考慮鄰近有限元節點的影響，因此，耦合界面物理量的連續性無法得到保證。同時，該方法要求耦合界面節點數與粒子數相同，兩者位置需要重合，限制了耦合界面周圍粒子和單元劃分，導致粒子與單元密度不匹配，對計算精度產生影響。

圖1粒子和單元節點的連接示意圖

為了有效地克服原算法中存在的缺陷，在有限元節點處設置了與有限元節點對應的背景粒子，背景粒子具有SPH粒子的屬性，其質量、位置、速度、加速度與相應有限元節點保持一致，對SPH粒子進行數值積分時，有限元節點以背景粒子的形式加入SPH鄰近搜索算法中，通過背景粒子轉化方式將有限元節點納入SPH粒子的鄰近搜索列表，SPH使得耦合界面周圍的粒子和單元劃分變得相當自由，從而避免SPH方法中光滑核函數被邊界截斷，也能去除邊界效應，保證了耦合界面處物理量的連續性[6]。改進的耦合算法不僅實施簡單，也能使計算精度提高。

1.2三維FE-SPH耦合算法基本理論

FE-SPH自適應算法分為四類，本文主要對耦合算法進行分析研究，FE-SPH耦合算法示意圖如圖2所示。由于耦合算法要確定12個速度，計算流程相比其他算法更為繁瑣復雜，因此，確保方程組正確求解相當關鍵[7]。建立上述方程需要根據線動量守恒，角動量守恒以及速度、位置在耦合點匹配的假設條件，計算步驟分為以下三步：

(1)檢測四節點法向速度的增量。

(2)檢測主面內速度增量。

(3)檢測主面三節點存儲角動量部分的速度增量。

圖2有限元與光滑粒子耦合算法

步驟(1)中，總穿透量表達式：

δtotal=δadjust+δattach

步驟(2)中，光滑粒子與主面相對移動距離為：

在圖2(a)中，s′軸、粒子和主面相對運動方向一致，三個節點的速度變化量ΔU1、ΔU2、ΔU3也沿s′軸方向。四節點在主面內的速度增量形式與法向速度增量形式相同，即：

上述各步驟中，步驟(1)中三個軸向保持角動量守恒；步驟(2)在s′軸方向保持角動量守恒，保證節點在主面上的位置和速度得到匹配；步驟(3)保證x′、y′軸方向保持角動量守恒，其中x′、y′方向角動量變化量與主面內粒子點速度變化的大小相等方向相反，且法向動量守恒。通過上述計算最終得到以下三個方程：

ΔV3=-M1ΔV1+M2ΔV2/M3

其中，ΔV1、ΔV2、ΔV3表示主面內粒子點速度變化量，如圖2(b)所示。

綜合以上計算，可以得到有限單元與光滑粒子耦合時需要調整的參量值。但是，由于耦合界面處周圍節點的影響，計算結果并不能一次調整到位，一般需要迭代1～5次之后才能使方程組最終達到平衡。

2數值建模

本文利用SPH粒子模擬具有大變形的流體，FE單元模擬具有不規則形狀的固體，固體模型在軟件LS-DYNA中建立，同時完成網格劃分，流體部分模型在LS-PREPOST完成。

2.1FE-SPH建模

為了減少計算成本，本文只建立1/2模型，如圖3所示。圖3(a)所示球體采用默認的Belytshco-Lin-Tsay殼單元進行映射式有限元建模，并在SECTION_SHELL中定義殼單元算法、積分規律和相關屬性，球體視為剛體，故其旋轉自由度受到約束，通過關鍵字MAT_RIGID來定義[8]；圖3(b)所示為FE單元，通過設置單元密度來模擬土壤與水域；圖3(c)所示模型主要由SPH粒子組成，用來模擬具有大變形的流體或固體，控制其粒子疏密程度相當重要，粒子過多可能使后期計算時間過長，出錯概率增大，最終導致仿真失敗；通過初期建模與網格劃分，最終得到需要的模型，如圖3(d)所示。模型材料具體密度見表1。

圖3FE-SPH耦合模型

表1模型密度

2.2關鍵參數選擇

FE-SPH耦合算法中流體與結構交界面處理是一個難點，因此，在關鍵字BOUNDARY設置上需要仔細分析，選取合適的參數。例如關鍵字CONTACT、CONTROL、DEFINE、INITIAL、LOAD、SET等必須根據數值模擬要求賦值，本文中球體初始速度為4.88m·s-1，球體加載曲線比例因子為1，終止時間為0.1s。

水的狀態方程采用Mie-Gruneisen狀態方程[9]，表示為：

P=Pc+PT=A(μ)+B(μ)·eVO

3仿真結果分析

3.1沖擊過程分析

球體沖擊土壤與水域過程基本相似，因此選取沖擊土壤過程分析研究。圖4所示為球體沖擊土壤歷程圖，從中可以看出FE單元與SPH粒子模擬土壤受沖擊后作用后變形效果明顯不同。由于球體與SPH粒子直接碰撞，SPH粒子相比FE單元變形量明顯要大，產生粒子濺射現象，效果更為直觀，可以看出SPH方法在模擬自由表面流動問題時極具優勢[10]。但是，同樣可以看出FEM單元在模擬土壤固體結構動力問題上有其獨特優勢，體現出有限元法求解結構動力問題時計算精度好、效率高、易于處理復雜邊界的優點。

圖4球體沖擊土壤歷程圖

圖5所示為SPH粒子受力云圖，從俯視圖中可以看出，SPH粒子在沖擊力作用下受力區域向四周延伸，其中心區域受力最大，變形量最大，粒子濺射現象最為明顯。

圖5粒子受力云圖

圖6所示為FE單元模擬土壤受力云圖，可以清晰體現土壤各個部分受力大小，充分體現了有限元法計算精度高，效率高的特點，數據結果對后期分析處理土壤結構動力問題相當重要[11]。

圖6FEM單元受力云圖

3.2數據結果分析

首先對球體撞擊土壤與水域時速度、加速度數據進行濾波處理，得到最終的球體速度、加速度時間歷程曲線，再分別選取FE單元與SPH粒子的一個單元得到有效應力隨時間變化曲線。通過對比球體撞擊土壤與水域計算結果，可以發現FE-SPH耦合算法在模擬沖擊動力學問題上的優越性，球體速度時間歷程曲線如圖7與圖8所示。

圖7球體沖擊土壤速度時間歷程曲線

圖8球體沖擊水域速度時間歷程曲線

圖7與圖8所示分別為球體沖擊土壤與水域速度時間歷程曲線，對比可以看出，在初始速度、能量、加載方式等因素相同的條件下，0.1s時間內球體撞擊水域條件下速度首先減到最小值(圖中取絕對值)，并逐漸趨于零，受其密度、粘性等因素影響，沖擊力比較均勻，速度變化比較平滑[12]。相比較而言，從圖7曲線變化趨勢可知，由于土壤密度較大，沖擊過程中阻力明顯高于水域，導致速度減小率更大，這與其他仿真模擬結果基本吻合。

球體加速度時間歷程曲線如圖9與圖10所示。

圖9球體沖擊土壤加速度時間歷程曲線

圖10球體沖擊水域加速度時間歷程曲線

比較圖9與圖10加速度時間歷程曲線，發現球體沖擊土壤過程中加速度峰值明顯高于水域條件下加速度峰值，最大值達到540.41m/s2，這與土壤密度較大有直接關系。同時，球體撞擊土壤與水域都出現二次加速度峰值，但水域條件下速度明顯要快，這是由于材料變形涉及沖擊波、復雜的材料特性，材料屬性不同，變形效應會有差別[13]。SPH粒子所產生的沖擊波間接作用于FE單元，使其發生形變，并儲存一部分能量，當FE單元吸收的能量重新釋放時，會使球體產生反向加速度，而這種結構動力學效應在撞擊土壤時更為明顯，本次模擬結果與其他已有數值仿真結果基本吻合。

SPH粒子與FE單元有效應力隨時間變化歷程曲線如圖11所示。

圖11有效應力時間歷程曲線

分別選取耦合界面處一個SPH粒子與對應的單元節點，經過數據處理得到如圖11所示的曲線，可以看出SPH粒子承受的有效應力高于FE單元，在0.05s之后，對FE單元施加的有效應力趨于恒定值，這與單元材料屬性、結構動力特性有關，仿真結果與已有數值結果基本吻合[14]。

4結束語

對于諸如民機水上迫降、返回艙回收、高速沖擊及爆炸問題，由于材料變形極大，因此是數值模擬中的一類難題。同時，材料變形涉及沖擊波、復雜的材料特性等因素，想要成功實現數值模擬比較困難[13]。綜合FEM法與SPH方法在模擬材料大變形問題時各自的優點，即FEM方法計算效率高和SPH方法模擬大變形能力強的優點，發展為比較成熟的FE-SPH耦合算法，該算法為高速沖擊及爆炸問題進行數值模擬提供了一種高效、準確的途徑。在對FE-SPH耦合算法不斷完善的基礎上，開發實用的可市場化的數值仿真軟件，進一步將其應用于工程實際問題，這是該研究方向后續發展的重點[14]。

[1] JOHNSON G R.Linking of Lagrangian particle methods to standard finite element methods for high velocity impact computations[J].Nuclear Engineering and Design,1994,150(2-3):265-274.

[2] JOHNSON G R,STRYK R A,BEISSEL S R.SPH for high velocity impact computations[J].Computer Methods in Applied Mechanics and Engineering,1996,139(1-4):347-373.

[3] MONAGHAN J J.SPH without a tensile instability[J].Journal of Computational Physics,2000,159(2):290-311.

[4] VIGNJEVIC R,DE VUYST T,CAMPBELL J C.A frictionless contact algorithm for meshless methods[C]∥Proceedings of International Conference on Computational & Experimental Engineering and Sciences(ICCES 2007),Miami,USA,November 27-29,2007,3(2):107-112.

[5] JOHNSON G R,STRYK R A.Conversion of 3D distorted elements into meshless particles during dynamic deformation[J].International Journal of Impact Engineering,2003,28(9):947-966.

[6] BONET J,KULASEGARAM S.Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations[J].International Journal for Numerical Method in Engineering,2000,47(6):1189-1214.

[7] JOHNSON G R,COOK W H.Fracture characteristics of three metals subjected to various strains,strain rates,temperatures and pressures[J].Engineering Fracture Mechanics,1985,21(1):31-48.

[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] BEISSEL S R,GERLACH C A,JOHNSON G R.Hypervelocity impact computations with finite elements and meshfree particles[J].International Journal of Impact Engineering,2006,33(1-12):80-90.

[10] JOHNSON G R,COOK W H.A constitutive model and data for metals subjected to large strains,high strain rates and high temperatures[C]∥Proceedings of the 7th International Symposium on Ballistics,Hague,Netherlands,April 19-21,1983:541-547.

[11] HU D A,HAN X,XIAO Y H,et al.Research developments of smoothed particle hydrodynamics method and its coupling with finite element method[J].Chinese Journal of Theoretical and Applied Mechanics,2013,45(5):639-652.

[12] XIAO Y H,HU D A,HAN X.A coupling algorithm of finite element and smoothed particle hydrodynamics[J].Chinese Journal of Computational Physics,2011,28(2):219-225.

[13] YANG G,LIANG C,LIU P,et al.Numerical simulation of ricochet problem of projectile penetrating into concrete target based on 3d FE-SPH adaptive coupling algorithm[J].Engineering Mechanics,2013,30(9):276-282.

[14] SCHWER L.Optional strain-rate forms for the Johnson-cook constitutive model and the role of the parameter epsion_0[C]∥Proceedings of the 6th European LS-DYNA Users’ Conference,Gothenburg,June 15-17,2007:1-17.

Study on Sphere Shock Problem Based on FE-SPH Coupling Algorithm

LIDan,ZHAOTingyu

(Civil Aviation Flight University of China, Guanghan 618307, China)

In order to make full use of SPH method to deal with the problem of large deformation and FEM, the new SPH-FEM coupling algorithm is proposed in this paper. In this algorithm, SPH particle discretization is used in the large deformation area to overcome the FEM distortion problem, and FEM unit discretization are used in other simulation area, both to avoid SPH boundary effects, but also to ensure the physical continuity of the coupling interface. The three-dimensional numerical simulation of the impact of the sphere on the soil and water is carried out. The Johnson-Cook model and the Mie-Gruneisen equation of state are used to calculate the soil and water parameters, and the velocity of the sphere in the soil and water , acceleration and effective stress data are compared. It can be seen that the simulation results are in good agreement with other existing simulation results, which proves that it is feasible and effective to deal with the structure-fluid interaction problem by using FE-SPH coupling algorithm. Coupling problem provides a new way.

FE-SPH coupling algorithm; boundary condition; fluid-solid coupling; numerical simulation

2016-12-02

李 旦(1992-),男,甘肅天水人,碩士生,主要從事飛機水上追降數值方法方面的研究,(E-mail)1171085449@qq.com; 趙廷渝(1965-),男,重慶人，教授,碩士,主要從事民用航空發動機性能及控制技術方面的研究,(E-mail)2494238988@qq.com

1673-1549(2017)02-0053-06

10.11863/j.suse.2017.02.11

O351

A