基于FE-SPH耦合算法的球體沖擊問(wèn)題研究

李旦， 趙廷渝

(中國(guó)民用航空飛行學(xué)院飛行技術(shù)學(xué)院, 四川廣漢618307)

基于FE-SPH耦合算法的球體沖擊問(wèn)題研究

李旦， 趙廷渝

(中國(guó)民用航空飛行學(xué)院飛行技術(shù)學(xué)院, 四川廣漢618307)

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

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

引言

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

FEM方法采用拉格朗日格式求解結(jié)構(gòu)域，ALE格式或者歐拉格式求解流體域，該方法雖然計(jì)算效率高，但需要?jiǎng)澐志W(wǎng)格，尤其在拓?fù)浣Y(jié)構(gòu)變化很大時(shí)進(jìn)行網(wǎng)格劃分困難，且后期計(jì)算處理較復(fù)雜，極有可能導(dǎo)致網(wǎng)格畸變[2]。因此，上述方法在模擬復(fù)雜自由表面流動(dòng)與結(jié)構(gòu)大變形相互作用問(wèn)題時(shí)存在困難。SPH方法起初用于處理三維開(kāi)放空間天體物理學(xué)問(wèn)題，由于該方法不受網(wǎng)格束縛，在模擬大變形以及自由表面流動(dòng)等問(wèn)題時(shí)能有效避免網(wǎng)格畸變，且能在拉格朗日框架下對(duì)流體和結(jié)構(gòu)同時(shí)進(jìn)行求解，目前已廣泛應(yīng)用于流體力學(xué)與連續(xù)固體力學(xué)。但是，在固體力學(xué)領(lǐng)域，用該方法模擬固體介質(zhì)存在拉伸不穩(wěn)定問(wèn)題，計(jì)算精度和效率相比FEM方法較低。由于SPH方法中近似核函數(shù)沒(méi)有Kronecher delta張量性質(zhì)，處理邊界條件比較困難[3]，因此，單獨(dú)采用SPH方法解決此類(lèi)問(wèn)題仍存在不足。

近年來(lái)，為解決SPH和FEM耦合界面的計(jì)算問(wèn)題，F(xiàn)E-SPH耦合算法被提出并得到快速發(fā)展，其基本思想為：采用SPH方法處理局部大變形區(qū)域，F(xiàn)EM計(jì)算處理小變形區(qū)域，便于綜合利用有限元法計(jì)算精度好、效率高、易于處理復(fù)雜邊界與SPH方法模擬自由表面流動(dòng)的優(yōu)勢(shì)[4]。目前，F(xiàn)E-SPH耦合算法在工程領(lǐng)域得到廣泛應(yīng)用，尤其是在數(shù)值模擬沖擊動(dòng)力學(xué)問(wèn)題方面，例如在處理高速?zèng)_擊及爆炸問(wèn)題上普遍采用該耦合算法，同時(shí)，在材料加工處理方面，F(xiàn)E-SPH耦合算法也得到了應(yīng)用。

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

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

三維FE-SPH耦合算法基本原理如圖1所示。Johnson G R提出將SPH粒子固結(jié)在有限元節(jié)點(diǎn)，通過(guò)固結(jié)的方式保證兩者之間傳遞相互作用力，且SPH粒子與對(duì)應(yīng)單元節(jié)點(diǎn)加速度相同，從而實(shí)現(xiàn)SPH與FE耦合[5]。但是，在計(jì)算界面粒子應(yīng)變率時(shí)該耦合算法沒(méi)有考慮鄰近有限元節(jié)點(diǎn)的影響，因此，耦合界面物理量的連續(xù)性無(wú)法得到保證。同時(shí)，該方法要求耦合界面節(jié)點(diǎn)數(shù)與粒子數(shù)相同，兩者位置需要重合，限制了耦合界面周?chē)Ｗ雍蛦卧獎(jiǎng)澐?，?dǎo)致粒子與單元密度不匹配，對(duì)計(jì)算精度產(chǎn)生影響。

圖1粒子和單元節(jié)點(diǎn)的連接示意圖

為了有效地克服原算法中存在的缺陷，在有限元節(jié)點(diǎn)處設(shè)置了與有限元節(jié)點(diǎn)對(duì)應(yīng)的背景粒子，背景粒子具有SPH粒子的屬性，其質(zhì)量、位置、速度、加速度與相應(yīng)有限元節(jié)點(diǎn)保持一致，對(duì)SPH粒子進(jìn)行數(shù)值積分時(shí)，有限元節(jié)點(diǎn)以背景粒子的形式加入SPH鄰近搜索算法中，通過(guò)背景粒子轉(zhuǎn)化方式將有限元節(jié)點(diǎn)納入SPH粒子的鄰近搜索列表，SPH使得耦合界面周?chē)牧Ｗ雍蛦卧獎(jiǎng)澐肿兊孟喈?dāng)自由，從而避免SPH方法中光滑核函數(shù)被邊界截?cái)?，也能去除邊界效?yīng)，保證了耦合界面處物理量的連續(xù)性[6]。改進(jìn)的耦合算法不僅實(shí)施簡(jiǎn)單，也能使計(jì)算精度提高。

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

FE-SPH自適應(yīng)算法分為四類(lèi)，本文主要對(duì)耦合算法進(jìn)行分析研究，F(xiàn)E-SPH耦合算法示意圖如圖2所示。由于耦合算法要確定12個(gè)速度，計(jì)算流程相比其他算法更為繁瑣復(fù)雜，因此，確保方程組正確求解相當(dāng)關(guān)鍵[7]。建立上述方程需要根據(jù)線動(dòng)量守恒，角動(dòng)量守恒以及速度、位置在耦合點(diǎn)匹配的假設(shè)條件，計(jì)算步驟分為以下三步：

(1)檢測(cè)四節(jié)點(diǎn)法向速度的增量。

(2)檢測(cè)主面內(nèi)速度增量。

(3)檢測(cè)主面三節(jié)點(diǎn)存儲(chǔ)角動(dòng)量部分的速度增量。

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

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

δtotal=δadjust+δattach

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

在圖2(a)中，s′軸、粒子和主面相對(duì)運(yùn)動(dòng)方向一致，三個(gè)節(jié)點(diǎn)的速度變化量ΔU1、ΔU2、ΔU3也沿s′軸方向。四節(jié)點(diǎn)在主面內(nèi)的速度增量形式與法向速度增量形式相同，即：

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

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

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

綜合以上計(jì)算，可以得到有限單元與光滑粒子耦合時(shí)需要調(diào)整的參量值。但是，由于耦合界面處周?chē)?jié)點(diǎn)的影響，計(jì)算結(jié)果并不能一次調(diào)整到位，一般需要迭代1～5次之后才能使方程組最終達(dá)到平衡。

2數(shù)值建模

本文利用SPH粒子模擬具有大變形的流體，F(xiàn)E單元模擬具有不規(guī)則形狀的固體，固體模型在軟件LS-DYNA中建立，同時(shí)完成網(wǎng)格劃分，流體部分模型在LS-PREPOST完成。

2.1FE-SPH建模

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

圖3FE-SPH耦合模型

表1模型密度

2.2關(guān)鍵參數(shù)選擇

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

水的狀態(tài)方程采用Mie-Gruneisen狀態(tài)方程[9]，表示為：

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

3仿真結(jié)果分析

3.1沖擊過(guò)程分析

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

圖4球體沖擊土壤歷程圖

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

圖5粒子受力云圖

圖6所示為FE單元模擬土壤受力云圖，可以清晰體現(xiàn)土壤各個(gè)部分受力大小，充分體現(xiàn)了有限元法計(jì)算精度高，效率高的特點(diǎn)，數(shù)據(jù)結(jié)果對(duì)后期分析處理土壤結(jié)構(gòu)動(dòng)力問(wèn)題相當(dāng)重要[11]。

圖6FEM單元受力云圖

3.2數(shù)據(jù)結(jié)果分析

首先對(duì)球體撞擊土壤與水域時(shí)速度、加速度數(shù)據(jù)進(jìn)行濾波處理，得到最終的球體速度、加速度時(shí)間歷程曲線，再分別選取FE單元與SPH粒子的一個(gè)單元得到有效應(yīng)力隨時(shí)間變化曲線。通過(guò)對(duì)比球體撞擊土壤與水域計(jì)算結(jié)果，可以發(fā)現(xiàn)FE-SPH耦合算法在模擬沖擊動(dòng)力學(xué)問(wèn)題上的優(yōu)越性，球體速度時(shí)間歷程曲線如圖7與圖8所示。

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

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

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

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

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

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

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

SPH粒子與FE單元有效應(yīng)力隨時(shí)間變化歷程曲線如圖11所示。

圖11有效應(yīng)力時(shí)間歷程曲線

分別選取耦合界面處一個(gè)SPH粒子與對(duì)應(yīng)的單元節(jié)點(diǎn)，經(jīng)過(guò)數(shù)據(jù)處理得到如圖11所示的曲線，可以看出SPH粒子承受的有效應(yīng)力高于FE單元，在0.05s之后，對(duì)FE單元施加的有效應(yīng)力趨于恒定值，這與單元材料屬性、結(jié)構(gòu)動(dòng)力特性有關(guān)，仿真結(jié)果與已有數(shù)值結(jié)果基本吻合[14]。

4結(jié)束語(yǔ)

對(duì)于諸如民機(jī)水上迫降、返回艙回收、高速?zèng)_擊及爆炸問(wèn)題，由于材料變形極大，因此是數(shù)值模擬中的一類(lèi)難題。同時(shí)，材料變形涉及沖擊波、復(fù)雜的材料特性等因素，想要成功實(shí)現(xiàn)數(shù)值模擬比較困難[13]。綜合FEM法與SPH方法在模擬材料大變形問(wèn)題時(shí)各自的優(yōu)點(diǎn)，即FEM方法計(jì)算效率高和SPH方法模擬大變形能力強(qiáng)的優(yōu)點(diǎn)，發(fā)展為比較成熟的FE-SPH耦合算法，該算法為高速?zèng)_擊及爆炸問(wèn)題進(jìn)行數(shù)值模擬提供了一種高效、準(zhǔn)確的途徑。在對(duì)FE-SPH耦合算法不斷完善的基礎(chǔ)上，開(kāi)發(fā)實(shí)用的可市場(chǎng)化的數(shù)值仿真軟件，進(jìn)一步將其應(yīng)用于工程實(shí)際問(wèn)題，這是該研究方向后續(xù)發(fā)展的重點(diǎn)[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-),男,甘肅天水人,碩士生,主要從事飛機(jī)水上追降數(shù)值方法方面的研究,(E-mail)1171085449@qq.com; 趙廷渝(1965-),男,重慶人，教授,碩士,主要從事民用航空發(fā)動(dòng)機(jī)性能及控制技術(shù)方面的研究,(E-mail)2494238988@qq.com

1673-1549(2017)02-0053-06

10.11863/j.suse.2017.02.11

O351

A