流固耦合問題的PD-SPH 建模與分析

姚學昊，黃 丹

(河海大學工程力學系，江蘇，南京 211100)

流固耦合(Fluid-structure interaction, FSI)問題廣泛存在于科學和工程領域，常涉及結構變形破壞以及流體破碎等復雜流動現象[1]，往往難以獲得解析解。隨著計算力學研究發展，任意拉格朗日-歐拉(Arbitrary lagrangian-eulerian, ALE)方法[2-3]、浸入邊界法(Immersed boundary method, IBM)[4]等各種數值方法被用于此類問題分析，然而在處理結構大變形、準確快速地追蹤自由表面和運動邊界等方面依然存在困難。

與網格類方法不同，光滑粒子動力學(Smoothed particle hydrodynamics, SPH)方法[1,5]是一種拉格朗日無網格法，它適合于追蹤自由表面和運動邊界，且在處理大變形問題時不會產生網格畸變。ANTOCI等[6]采用SPH 方法模擬流體與彈性結構的相互作用問題，獲得了較理想的結果。KHAYYER 等[7]和ZHANG 等[8]分別提出了一種不可壓縮SPH 方法以及多分辨率SPH 方法來提高流固耦合問題的求解精度與效率。但SPH 方法在模擬固體時仍需要一些修正技術以克服拉伸不穩定問題[1,5]等固有缺陷。近年來，研究人員開發了諸多SPH 與其他方法耦合的算法，如SPH 與有限元法(Finite element method, FEM)[9-10]、光滑點插值法(Smoothed point interpolation method, SPIM)[11]的耦合等。然而，上述算法主要針對含結構變形的流固耦合問題，在涉及固體損傷與破壞的相關問題中具有局限性。近場動力學(Peridynamics, PD)[12-14]是一種新興的非局部方法，采用空間積分代替傳統微分方程，在分析不連續力學問題時展現出獨特優勢，可應用于多種大變形和結構破壞[15-17]、流固耦合結構破壞[18-19]等模擬。采用PD-SPH 耦合方法來分析流固耦合問題，能同時發揮PD 方法模擬固體破壞和SPH 方法追蹤流體運動界面的優勢。

近來，部分學者開展了PD 與SPH 方法結合的初步研究。ZHOU 等[20]結合兩種方法提出SPD(Smoothed peridynamics)方法并應用于大變形與斷裂問題。FAN 等[21-22]提出一種PD-SPH 耦合模型模擬爆炸引起的土壤破碎問題；SUN 等[23]最近提出一種基于移動虛粒子的PD-SPH 耦合方法，并應用于流固耦合問題。

本文建立一種能求解流固耦合及結構破壞的新型PD-SPH 耦合模型。針對流-固界面處理，提出一種基于虛粒子和排斥力的耦合算法，既能防止流體粒子穿透交界面，又能對交界面處流體粒子進行邊界缺陷修正，提高計算精度。通過模擬流體作用下的彈性結構變形問題，驗證了耦合方法的可行性和有效性，并進一步開展了涉及結構斷裂破壞的復雜流固耦合問題全過程仿真分析。

1 數值模型

1.1 固體模型

1.1.1 常規態型PD 基本方程

如圖1 所示，將空間域R離散為一系列包含物性信息且具有相互作用f的PD 粒子，在任意時刻t，粒子Xi的運動方程可寫為：

圖1 PD 模型示意圖Fig. 1 Sketch for PD model

PD 理論分為鍵型(bond-based)PD 理論、常規態型(ordinary state-based)PD 理論以及非常規態型(non-ordinary state-based)PD 理論。由于鍵型與非常規態型PD 理論通常分別存在泊松比限制及數值不穩定問題，本文基于常規態型PD 理論進行固體域求解。此時，點對相互作用f可表示為：

1.2 流體模型

1.2.1 控制方程

1.3 流-固交界面處理方案

PD 和SPH 方法均通過粒子離散計算域，故可通過基于粒子-粒子接觸的耦合方案處理流-固交界面。圖3 為本文提出的基于虛粒子和排斥力的PDSPH 耦合方案示意圖。在該方案中，PD 粒子將作為兩種類型的耦合虛粒子(邊界虛粒子和內部虛粒子)參與SPH 計算。同時，反作用力將作用于PD粒子從而實現SPH 粒子對PD 粒子的影響。

該耦合方案具體計算流程如圖4 所示。

圖4 基于虛粒子和排斥力的PD-SPH 耦合方案流程Fig. 4 Flow chart for PD-SPH coupling scheme based on virtual particle and repulsive force

2 算例分析

2.1 靜水壓力作用下的鋁板變形

為了檢驗耦合方案的有效性、穩定性及計算精度，首先研究了靜水壓力作用下的鋁板變形問題[9]。如圖5 所示，在寬1.0 m 的水箱中，箱內水深H= 2.0 m，密度為 1000 kg/m3；水箱底部為厚度d= 0.05 m 的鋁板，其密度是 2700 kg/m3，彈性模量67.5 GPa，泊松比0.34。

圖5 靜水壓力下的鋁板示意圖Fig. 5 Sketch for aluminium plate under hydrostatic pressure

初始時刻，鋁板突然受到水柱的靜水壓力載荷而發生變形。經過短時間的初始振蕩，系統最終達到具有特定靜態變形的平衡狀態。根據理論解[9]，2.0 m 高水柱產生的靜水壓力載荷作用下鋁板跨中撓度大小為 -6.85×10-5m。在模擬計算中：人工粘性參數απ= 0.03 ； 參考聲速c0=190 m/s；光滑長度h=1.33Δx；近場范圍半徑 δ=3Δx；時間步長 Δt=1×10-6s；總模擬時間為1 s。流體粒子與固體粒子的初始間距一致，粒子空間分辨率分別取L/Δx=100、L/Δx=80和L/Δx=60，對應的粒子間距為0.01 m、0.0125 m 以及0.0167 m，此時SPH 粒子總數分別為21060、13648、7836，PD 粒子總數分別為525、340、195。本算例在配有主頻為2.7 GHz 的Intel CPU 的計算機上完成，計算總耗時分別為8 h 31 min、6 h 4 min、4 h 32 min。

圖6 給出了模擬得到的t=1 s 時不同空間分辨率(粒子間距)下流體壓力云圖以及放大1000 倍后的鋁板變形圖。可見，基于PD-SPH 的流固耦合求解方法在不同空間分辨率條件下均可得到光滑的流體壓力場和鋁板撓度場。

圖6 t = 1 s 時不同空間分辨率下流體壓力云圖和鋁板變形Fig. 6 The fluid pressure contours and the deflection of the aluminium plate at t = 1 s with different spatial resolutions

圖7 為三種不同空間分辨率下鋁板跨中撓度的時間歷程，表1 則給出了不同分辨率下跨中撓度數值解及其與解析解的相對誤差。結合圖7 與表1 可以看出，空間分辨率L/Δx=60時，由于粒子數較少，PD 和SPH 方法本身的計算精度較低；隨著空間分辨率的增加，PD-SPH 耦合方法的跨中撓度計算結果逐漸向靜態理論解析解收斂，并最終在L/Δx=100時，與解析解基本一致，相對誤差僅為3.07%。由此表明：本文提出的PD-SPH 耦合算法能有效、準確模擬準靜態流固耦合問題。

表1 不同空間分辨率下鋁板跨中撓度計算誤差Table 1 Computational error of mid-span deflection of aluminium plate with different spatial resolutions

圖7 鋁板跨中撓度Fig. 7 Mid-span deflection of aluminium plate

2.2 潰壩水流沖擊彈性板

如圖8 所示，寬為0.584 m 的水箱左側存在高H= 0.292 m，寬L= 0.146 m 的水柱，水柱在重力作用下倒塌產生潰壩流，并沖擊水箱中高he=0.08 m，厚度為w= 0.012 m 的彈性板，隨后發生劇烈的流固耦合作用。該算例中，水的密度為1000 kg/m3；彈性板的密度為 2500 kg/m3，彈性模量為1.0 MPa，泊松比為0.0。模擬中，人工粘性參數設為 απ= 0.02 ； 參考聲速c0=35 m/s；粒子間距Δx= 0.002 m，對應的SPH 粒子與PD 粒子總數分別為12744 和246；光滑長度h=1.5Δx；近場范圍半徑 δ=3Δx； 時間步長為 Δt=5×10-6s；總模擬時間為0.75 s。使用與算例2.1 相同的計算機配置，計算總耗時為2 h 20 min。

圖8 潰壩水流沖擊彈性板的初始條件示意圖Fig. 8 Sketch for initial conditions in dam-break flows impacting the elastic plate

為了對PD-SPH 結果進行定量驗證，在PDSPH 模擬中計算了彈性板頂端A點(如圖8 所示)的水平位移時間歷程，與其它已有文獻結果[29-32]的對比如圖9 所示。可以看出：在水流沖擊彈性板前期(0.4 s 前)，PD-SPH 方法的位移計算結果與其它方法所得結果吻合良好，且彈性板的位移最大值以及最大位移的出現時間基本一致，其中，PD-SPH 方法在0.244 s 時獲得彈性板水平位移最大值，約0.046 m；在0.4 s 以后，受自由表面融合現象的影響，不同方法得到的位移結果存在明顯差異，PD-SPH 模擬結果與粒子有限元法(PFEM)[29]、SPH-再生核粒子法(RKPM)[32]模擬結果更為接近。

圖9 A 點水平方向位移時程Fig. 9 Time history of horizontal displacement of point A

圖10 為不同時刻下，潰壩水流沖擊彈性板產生的流體飛濺與結構變形。可以看出：0.16 s 時，水流沖擊彈性板并沿板的左側向上爬升，同時，彈性板在沖擊力作用下出現彎曲變形；0.26 s 時，水流完全覆蓋彈性板左側，使其產生較大變形；0.42 s 時，水流撞擊右側壁面形成射流，彈性板則出現較大幅度的回彈現象。同時，圖10 還給出了文獻中PFEM 方法[29]和SPH 方法[30]的模擬結果，比較看出：每個時刻3 種方法得到的流體自由表面形狀和彈性板變形都是相近的。此外，PD-SPH方法也計算得到了光滑的流體壓力場，并且與PFEM結果相似。由此表明：本文提出的PD-SPH 耦合方法對劇烈流固耦合問題的模擬是可行的。

圖10 不同時刻彈性板變形和流體飛濺的PFEM[29]、SPH[30]和PD-SPH 結果的比較Fig. 10 Comparison of PFEM[29], SPH[30] and PD-SPH results of plate deformation and fluid splashing at different times

2.3 流體作用下的結構破壞

為了驗證本文PD-SPH 方法在流固耦合破壞問題模擬方面的可行性，本節模擬結構在潰壩流體作用下的運動、破壞過程。如圖11 所示，水箱長0.8 m，高0.4 m，其內部水柱和擋板的幾何尺寸、材料參數均與算例2.2 相同。為了實現擋板在水流作用下的斷裂破壞，假定擋板材料強度極限較低。本算例中，SPH 和PD 粒子總數分別為14267 和258，總模擬時間為0.8 s，計算總耗時為2 h 27 min。

圖11 流體作用下結構破壞的初始條件示意圖Fig. 11 Sketch for initial conditions for structural failure under fluid action

圖12 為不同時刻下結構變形破壞和流體飛濺的模擬結果，圖13 則給出了結構質心A(如圖11 所示)的坐標隨時間的變化。在0.18 s 前，流體運動過程以及結構變形與算例2.2 完全一致；0.18 s 時，結構底部開始斷裂(如圖12(a)所示)，質心A的y坐標時程曲線出現瞬時震蕩(如圖13(b)所示)；在水流的持續作用下，裂紋不斷擴展，產生的應力波也在水體內部傳播，導致流體壓力場振蕩，0.19 s 時裂紋呈現明顯的V 字型，上部擋板僅剩1 個粒子與底部連接(如圖12(b)所示)；最終在0.2 s 時，上部擋板與底部完全斷開(如圖12(c)所示)；隨后，擋板在水流作用下向水箱右側運動，于0.318 s 時撞擊水箱底部(如圖12(d)和圖13(b)所示)，并出現反彈現象(接觸力由式(18)計算)；經過幾次接觸回彈，擋板開始沿水箱底部滑動，0.524 s時擋板撞擊水箱右側(如圖12(f)和圖13(a)所示)，并開始向左滑動。最終，擋板將在水箱底部保持靜止狀態。本文模擬所得擋板運動過程與文獻[33]基本一致，由此表明：本文提出的PD-SPH 方法可應用于流固耦合破壞問題求解，且能夠較好地預測流體運動過程和結構變形破壞以及破壞后部分結構的運動。

圖12 不同時刻結構變形破壞和流體飛濺的模擬結果Fig. 12 Simulation results of structural deformation and failure and fluid splashing at different times

圖13 結構質心A 的坐標時間歷程Fig. 13 Time history of the coordinate of centroid A

3 結論

本文建立了一種新的PD-SPH 耦合模型求解流固耦合作用下的結構破壞問題，主要工作和結論如下：

(1)分別采用PD 方法與SPH 方法離散固體域和流體域，可充分發揮兩種方法的優勢。通過基于虛粒子和排斥力的耦合算法處理流-固交界面，能夠對流固界面處的流體粒子施加速度邊界，修正核函數缺陷，提高計算精度。

(2)采用PD-SPH 耦合方法模擬了靜水壓力作用下的鋁板變形問題以及潰壩水流沖擊彈性板問題，所得結構變形和流體運動過程與解析解或文獻結果吻合較好，驗證了本文PD-SPH 耦合方法在模擬復雜流固耦合問題方面的適用性。

(3)對流體沖擊作用下的結構破壞過程分析表明：所提出的PD-SPH 耦合方法易于實現流固耦合作用下結構從變形到破壞乃至破壞后部分結構運動的全過程仿真模擬，可為流-固耦合結構破壞問題的研究提供新參考。