SPH 統一算法對自由流體沖擊彈性結構體問題模擬

第一作者沈雁鳴男，博士生,1984年生

通信作者陳堅強男，研究員,1966年生

郵箱:jq-chen@263.net

SPH統一算法對自由流體沖擊彈性結構體問題模擬

沈雁鳴，何琨，陳堅強，袁先旭

(中國空氣動力研究與發展中心，四川綿陽621000)

摘要：對自由流體沖擊彈性結構體涉及的自由界面運動、彈塑性變形、流固耦合作用等問題進行模擬時要求能有效計算介質運動變形與流固間耦合作用。利用光滑粒子流體動力學方法將流體、彈性體視為連續介質，用統一的N-S方程與虛擬壓縮狀態方程求解。基于初始粒子類型劃分提出新的介質接觸算法，通過在結構體表面增加結構體表面粒子，使流體與彈性體介質間相互作用實現穩定傳遞。計算結果表明，所用SPH統一算法與接觸算法能準確有效模擬自由流體沖擊彈性結構體變形回彈過程。

關鍵詞：流固耦合；變形回彈；SPH方法；接觸算法

收稿日期：2014-06-24修改稿收到日期：2014-08-07

中圖分類號:O353.4文獻標志碼：A

基金項目：國家自然科學基金重點項目(51239007)

Numerical simulation of free surface flow impacting elastic structure with SPH uniform method

SHENYan-ming,HEKun,CHENJian-qiang,YUANXian-xu(China Aerodynamics Research and Development Center, Mianyang 621000, China)

Abstract:Free surface flow impacting elastic structure is a common problem, usually referring to interface track, elastic or plastic structural deformation and fluid-structure interaction. To simulate these problems, it has to capture the deformation of interfacial flow and compute the fluid structure coupling interaction accurately. A smoothing particle hydrodynamics (SPH) uniform method where both fluid and solid phases are described by SPH was presented and the uniform continuous medium governing equations and week pressure state equation were solved. A new contact algorithm with the structure surface particles treated as virtual particles was developed. The virtual particles were defined in initial particle configuration, containing both fluid and structure particle characters and producing repulsive force to prevent penetration. Comparing the numerical results with experimental ones, it is demonstrated that the SPH uniform method with the new contact algorithm can be used to simulate the free surface flow impacting elastic structure problem easily and accurately.

Key words:fluid-structure interaction;deformation resilience; smoothed particle hydrodynamics; contact algorithm

自由流體沖擊彈性結構體即帶自由界面運動的流體與彈性結構體間發生的流體結構體耦合作用(fluid-structure interaction，FSI)廣泛存在于人類生產生活中，如船舶砰擊、波浪沖擊海洋平臺等。此類問題涉及自由界面運動、彈塑性變形、流固耦合作用等，對其進行數值模擬時要求能跟蹤流體介質界面運動變形及求解介質間耦合作用。

光滑粒子流體動力學方法[1](Smoothed Particle Hydrodynamics，SPH)作為成熟的無網格Lagrange方法廣泛用于沖擊動力學領域[2-4]，通過改造可用于運動變形及流固耦合問題求解。Oger等[5]通過流體粒子壓力投影方法獲得剛體表面壓力，實現剛體入水沖擊時流固耦合作用的數值計算；Gong等[6]則發展出新型固壁面壓力處理方法用于楔形體入水沖擊計算。以上計算中均將結構體視為剛體，未考慮結構體彈性變形。因沖擊強烈時結構體可能產生明顯變形。

為處理結構體彈性問題可利用SPH方法與有限元耦合計算。LS-DYNA[7]軟件中已嵌入SPH算法，可與有限元算法耦合處理各種流彈性問題。肖毅華等[8-9]已開展該項工作。該方法可提高彈性計算精度及計算效率，但SPH方法與FEM耦合計算時要求處理好單元與粒子間耦合算法及接觸算法。

具有材料強度的彈性體、水、空氣等均屬于廣義流體力學范疇，可用連續介質力學中守恒方程作為控制方程，而描述應力應變關系的本構方程、壓力狀態方程也具有通用性，因此稍加改進即可將所有氣液固介質統一用光滑粒子流體動力學方法求解。通過流固邊界處理、提高計算精度，便可獲得較準確的計算結果。Antoci等[10-11]先后已嘗試成功，Amini等[12]通過改進固液介質間接觸模型使交界面更平穩，Liu等[13]提出耦合動態流固邊界處理方法提高計算精度。該方法處理方式簡單有效，易編程實現，且方便程序并行改造，因此在流固耦合問題應用中逐漸增多。

本文將流體、彈性體視為連續介質，用統一的N-S方程與虛擬壓縮狀態方程求解。基于初始粒子類型劃分提出新的介質接觸算法，即在結構體表面增加表面粒子，使流體與彈性體介質間相互作用，實現穩定傳遞。介紹SPH統一算法下連續介質控制方程、本構方程、狀態方程、人工應力穩定項及新提出的流體與結構體交界面接觸算法，利用典型算例驗證所提SPH統一算法及流體結構體交界面接觸算法。

1數學物理模型

1.1連續介質控制方程及狀態方程

具有材料強度的彈性體、水、空氣等均屬于廣義流體力學范疇，可用連續介質力學中守恒方程作為控制方程，即

(1)

(2)

由此，可利用光滑粒子法公式對以上方程進行離散及插值轉換，獲得方程組的光滑粒子離散格式，即

(3)

(4)

在兩相流動或流固耦合計算中，雖mi≠mj，ρi≠ρj，但式(3)、(4)仍可有效避免不同介質質量及密度差所致插值誤差。

總應力σαβ主要有壓力項P及偏應力Sαβ構成，即

σαβ=-Pδαβ+Sαβ

(5)

式中：P為壓力值或正應力，由狀態方程給出；Sαβ為偏應力，流體力學中也稱粘性應力項ταβ。

應力與應變成比例，比例系數為動力粘滯率即粘性系數μ，單位Pa·s，與時間無關

(6)

對彈性體，偏應力項稱剪應力項，材料力學中剪應力與剪應變關系非單純比例關系，而與時間相關，比例系數為材料剪切模量，單位Pa。剪應變率滿足

(7)

由式(6)、及連續介質控制方程組可將彈性體與流體耦合一起計算。

正應力P通過狀態方程中粒子密度、內能等狀態值計算獲得。然而，在不可壓縮連續介質問題中，介質實際狀態方程限制時間步長取值大小，即時間步長不能太小。如何準確高效計算壓力項為不可壓縮計算的主要任務。大部分液體、固體可作為“不可壓縮”連續介質，由物理意義講，均為可壓縮，僅壓縮率非常低。介質密度變化與聲速平方成正比，即

(8)

式中：B=ρ0c2/γ，c為實際聲速值，ρ為水密度。

引入人工壓縮率概念[14]，即將不可壓縮介質或微可壓縮介質視為可壓縮介質，適當放大壓縮性，調整聲速取值，放寬時間步長限制。對水流，由于采用人工壓縮率后壓縮變化仍較小，基本不改變水流狀態。

1.2人工應力穩定項

Swegle等[15]發現SPH方法存在不穩定性問題并進行證明。利用SPH方法模擬彈性體彈塑性變形時，存在明顯拉伸應力不穩定性。為此，Gray等[16]在動量方程中添加人工應力穩定項，即

(9)

式中：f=W(rij)/W(Δx)，Δx為粒子間距。

(10)

(11)

據經驗，n=4，ε=0.3。

1.3流體與結構體交界面接觸算法

該算法為流彈性計算中的難點。對流體與彈性體耦合計算，須實現兩種不同介質相互作用力的有效傳遞。通過插值計算判斷彈性體與流體交界面粒子進而獲取交界面空間信息，再通過積分求解流固作用。該方法在復雜形狀下適用性有限，如當彈性體具有拐角或曲線形狀時則無法通過簡單插值判斷交界面。本文發展的新接觸算法，基于初始粒子類型劃分，在結構體表面增加粒子，該種粒子具有流體、彈性體粒子兩種屬性，分別參與流體、彈性體粒子間插值計算。通過利用該種粒子使流體與彈性體介質間相互作用實現穩定傳遞。

圖1為本文對流體結構體耦合作用問題處理方式，其中方形粒子為離散的流體粒子，參與流體粒子質量及動量方程計算；菱形粒子為固壁邊界粒子，受完全約束條件，只通過流體粒子進行壓力差值，壓力值參與流體粒子動量方程計算；三角粒子為彈性體內部粒子，只參與彈性體粒子及彈性體邊界粒子質量、動量方程計算；外層星形粒子為彈性體邊界粒子，與彈性體粒子插值計算時具備彈性體粒子屬性，與流體粒子插值計算時視為固壁粒子，靜水壓力由周圍流體粒子壓力值插值獲得。由此能較便捷處理各類復雜流體-剛體/結構體耦合問題，對結構體外形復雜程度無限制條件，且較易擴展至三維。

圖1　對流體結構體耦合作用問題處理方式 Fig.1 Contact algorithm for fluid-structure interaction

為防止流體粒子穿透結構體，本文參考固壁邊界處理方法，對結構體表面粒子施加Lennard-Jones排斥力[17]，排斥力表達式為

(12)

2結果與分析

為提高計算精度及穩定性，利用密度最小二乘插值、人工粘性與人工熱流及速度平均項等方法。利用潰壩沖擊及彈性懸臂梁振動等單介質運動變形問題進行算例考核，并進一步模擬彈性攔壩板泄流及潰壩沖擊彈性板后變形回彈過程。

2.1潰壩流動

計算模型參數為：初始時刻水柱高H=1 m，水面寬L=2 m；水箱高D=3 m，長d=5.366 m；水粒子數5 000個，粒子間距Δx=0.02 m，光滑長度h=1.33Δx，計算時間步長Δt=10-5s，狀態方程中水聲速值為14(gH)0.5。不同時刻水柱變形數值模擬結果見圖2。由圖2看出，水柱在重力作用下迅速垮塌，遇垂直固壁后產生猛烈撞擊，水頭沿垂直壁面向上爬升并向外翻滾，產生的浪頭與后部水面再次發生拍擊產生二次浪花飛濺。

圖2　潰壩問題不同時刻水體壓力圖 Fig.2 Time evolution of water pressure

距垂直板1.653 m處水位高度變化見圖3。由圖3看出，本文計算所得水位高度與試驗結果在大部分時間段相同。

圖3　垂直面前方1.653 m處水位高度與時間關系 Fig.3 Water level at Δx= 1.653 m front vertical wall

2.2彈性懸臂梁

懸臂彈性梁問題為典型的彈性運動驗證算例，初始時刻給定彈性梁一定運動速度，可觀測獲得彈性梁振動周期、振幅及衰減系數。取彈性梁長L=0.2 m，直徑H=0.02 m，一端固定于壁中。給自由端初速度使整個彈性梁產生振動，其在不同位置的速度滿足

(13)

f(x)=[cos(kl)+cosh(kl)][cosh(kx)-cos(kx)]+

[sin(kl)-sinh(kl)][sinh(kx)-sin(kx)]

(14)

圖4　懸臂彈性梁彈性振動過程壓力云圖 Fig.4 Time evolution of the pressure of elastic cantilever

將所得彈性梁振動周期、最大振幅與理論值進行比較。本文的A=0.025 7 m，T=0.289 s，與文獻結果Ar= 0.025 m，Tr=0.287 7 s[10-11]相符，較理論計算值Aa=0.023 m，Ta=0.254 s略大。

2.3彈性攔壩板

文獻[10]通過實驗研究彈性板在隨時間變化的水壓作用下變形，并與計算進行比較。圖5為該試驗二維示意圖，圖中幾何參數為：A=0.1 m；H=0.14 m；L= 0.079 m；s=0.005 m。

圖5　彈性攔壩板算例示意圖 Fig.5 Initial configuration of elastic gate

本文利用光滑粒子統一算法，結合結構體表面粒子處理流體與交界面接觸算法，對彈性攔壩板進行模擬。粒子間距Δx=0.001 m，時間步長Δt=2×10-6s，不考慮表面張力作用；彈性板密度1.1×103kg/m3，體積模量k=1.67×107Pa，剪切模量=3.57×106Pa。水密度1×103kg/m3，聲速值c=48.3 m/s。共15 088個粒子，其中水粒子14 000，固壁邊界粒子708，彈性板粒子300，彈性邊界粒子80。不同時刻流場及彈性板變形計算結果見圖6。由圖6看出，在水壓作用下，彈性板自由端逐漸偏移且發生彈性變形。其中彈性板固定端附近彎曲變形較嚴重，下半部變形較小，近似呈剛體運動；彈性板在水壓作用下形成的水流出口先逐漸增大到最大值后保持一定時間不變，而后隨水位降低水壓不斷減小，使受力減小、彈性彎曲程度降低，水流出口又逐漸減小。對應時刻中計算結果與試驗在水位高度及彈性變形上基本符合。

圖6　流場和彈性板變形情況 Fig.6 Comparison of experimental results with SPH simulation

圖7　彈性攔壩板自由端位移曲線 Fig.7 Horizontal and vertical displacements of the free end

彈性板自由度在x，y方向位移變化曲線見圖7。由圖7看出，計算結果與試驗結果符合度較好，x，y方向位移從初始時刻開始增加，約0.15 s達最大值后緩慢下降。因計算所得彈性板變形規律與試驗結果基本一致，彈性板彎曲程度直接決定水流排泄速度，故彈性板內水壩水位高度變化計算結果與試驗結果也應較符合，此由圖8水位隨時間變化曲線比較結果也可獲得印證。

圖8　彈性板后水位高度隨時間變化曲線 Fig.8 Water level behind the gate

2.4潰壩沖擊彈性板

為確定流體彈性體剪切模量在流彈性計算中的影響，本文利用潰壩沖擊彈性板問題進行計算比較。具體計算參數為：容器內左側初始時刻存放水體寬L= 0.146 m，高H=2L，容器總寬W=4L。在水體前方L處有一長Lp=0.08 m，寬W=0.012 m彈性板密度1.1×103kg/m3。水密度為1×103kg/m3，取聲速c= 34.1m/s。粒子間距Δx=0.002 m，粒子總數11 810個，其中水粒子10 658，彈性板粒子240個(含彈性板邊界粒子84個)，時間步長Δt=1×10-5s。本文用三種彈性模量6×106N/m2，1.2×107N/m2及2.38×107N/m2，并用剛體模型進行比較。

具有不同彈性模量彈性板受水體沖擊后計算結果見圖9。由圖9看出，剪切模量大小會直接決定彈性體受水體作用過程中彎曲變形程度，剪切模量越小彈性板最大彎曲幅度越大。此由圖10不同剪切模量下彈性板振動幅度時間歷程曲線比較結果可獲得證明。所有彈性板基本均在2.4 s左右達最大彎曲幅度后逐漸回彈，在0.92 s出現左側彎曲極大值。

圖9　不同彈性模量對潰壩沖擊影響比較，自上而下：6.0×10 6 N/m 2，1.2×10 7 N/m 2，2.38×10 7 N/m 2及剛體 Fig.9 Breaking dam on a elastic baffles with different elastic modulus：6.0×10 6 N/m (first row), 1.2×10 7 N/m 2 (second row), 2.38×10 7 N/m 2(third row) and rigid baffle(last row)

圖10　彈性板振動幅度時間歷程曲線比較 Fig.10 Time history of the displacement of elastic plate

總之，不同剪切模量對彈性體變形均有明顯影響，但對水體沖擊變形過程影響較小，而彈性體變形與水體流動規律亦較相似。尤其將彈性體替換為剛體計算時所得水流狀態及壓力云圖與剪切模量為2.38×107N/m2時基本一致。常規合金材料剪切模量為1×109N/m2量級，因此在流體動壓不特別大、且材料受沖擊面與材料橫截面比值不很大時，基本可視其為剛體。

3結論

(1)通過用改進的光滑粒子流體動力學統一算法，將流體、彈性體視為連續介質，用N-S方程及虛擬壓縮狀態方程對自由流體沖擊彈性結構體過程統一求解，并用人工應力穩定項抑制彈性變形中出現的拉伸應力不穩定性現象。

(2)為實現流體與彈性體耦合計算中介質間相互作用有效傳遞，本文基于初始粒子類型劃分，發展出新的接觸算法模型，通過在彈性結構體表面增加結構體表面粒子，使粒子具有流體、彈性體粒子兩種屬性，分別參與流體、彈性體粒子間的插值計算。用結構體表面粒子使計算中流體與彈性體介質間作用得到準確穩定傳遞。

(3)通過二維潰壩及彈性懸臂梁算例考核，證明SPH統一算法對自由流動界面追蹤及連續介質彈性變形計算的有效性，并進行彈性攔壩板泄流、潰壩沖擊彈性板變形回彈模擬。結果表面，本文所用SPH統一算法及介質接觸算法能準確、有效、便捷模擬自由流動沖擊彈性結構體問題，計算結果與試驗結果吻合度較好。

參考文獻

[1]Gingold R A, Monaghan J J. Smoothed particle hydro-dynamics: theory and application to non-spherical stars[J]. Monthly Notices R. Astronomy Soc., 1977, 181: 375-389.

[2]許慶新,沈榮瀛,周海亭. SPH方法在剪切式碰撞能量吸收器中的應用[J]. 振動與沖擊,2007,26(9): 108-111.

XU Qing-xin, SHEN Rong-ying, ZHOU Hai-ting. Applying sph method to shearing energy absorbers[J]. Journal of Vibration and Shock, 2007, 26(9): 108-111.

[3]鄭曉梅,徐志宏,湯文輝,等. 基于二階Godunov格式的SPH方法模擬彈塑性波的傳播[J].振動與沖擊,2011,30(2): 60-64.

ZHENG Xiao-mei, XU Zhi-hong, TANG Wen-hui, et al. Simulation of elastic-plastic wave propagation based on smoothed particle hydrodynamics method with second order Godunov scheme [J]. Journal of Vibration and Shock, 2011, 30(2): 60-64.

[4]初文華,張阿漫,明付仁,等. SPH-FEM耦合算法在爆炸螺栓解鎖分離過程中的應用[J]. 振動與沖擊,2012,31(23):197-202.

CHU Wen-hua, ZHANG A-man, MING Fu-ren, et al. Application of three-dimensional SPH-FEM coupling method in unlocking process of an explosion bolt [J]. Journal of Vibration and Shock, 2012, 31(23): 197-202.

[5]Oger G, Doring M, Alessandrini B,et al. Two-dimensional SPH simulations of wedge water entries[J]. Journal of Computational Physics, 2006(213): 803-822.

[6]Gong K, Liu H, Wang B L. Water entry of a wedge based on SPH model with an improved boundary treatment[J]. Journal of Hydrodynamics, 2009, 21(6): 750-757.

[7]John O H. Theoretical manual[M]. Livermore Software Technology Corporation(LSTC), 1998.

[8]肖毅華,韓旭,胡德安. 流體與結構相互作用問題的FE-SPH 耦合模擬[J]. 應用力學學報,2011, 28(1):13-18.

XIAO Yi-hua, HAN Xu, HU De-an. A coupling algorithm of finite element and smoothed particle hydrodynamics[J]. Chinese Journal of Computational Physics,2011, 28(1): 13-18.

[9]陳小波. 近海風機結構體系環境荷載及動力響應研究[D]. 大連：大連理工大學, 2011.

[10]Antoci C, Gallati M, Sibilla S.Numerical simulation of fluid- structure interaction by SPH[J].Computers and Structures, 2007, 85: 879-890.

[11]Rafiee A, Thiagarajan K P. An SPH projection method for simulating fluid-hypoelastic structure interaction[J].Comput Methods Appl Mech. Engrg., 2009, 198: 2785-2795.

[12]Amini Y, Emdad H. Farid M. A new model to solve fluid-hypo-elastic solid interaction using the smoothed particle hydrodynamics (SPH) method [J]. European Journal of Mechanics B/Fluid,2011, 30: 184-194.

[13]Liu Mou-bin, Shao Jia-ru, Li Hui-qi. Numerical simulation of hydro-elastic problems with smoothed particle hydro-dynamics method[J]. Journal of Hydrodynamics,2013,25(5):673-682.

[14]Monaghan J J. Simulating free surface flows with SPH[J]. Journal of Computational Physics, 1994, 110: 399-406.

[15]Swegle J W, Attaway S W. On the feasibility of using smoothed particle hydrodynamics for underwater explosion calculations[R]. New Mexico,US: SAND95-0311, Sandia, Albuquerque, 1995.

[16]Gray J P, Monaghan J J, Swift R P. SPH elastic dynamics [J]. Comput. Methods Appl. Mech. Eng., 2001, 190: 6641-6662.

[17]Liu G R, Liu M B. Smoothed particle hydrodynamics: a meshfree particle method [M]. World Scientific Publishing Co Pte Ltd,2003.