基于SBFE虛粒子邊界的結構入水數值模擬

李上明，趙梓斌

(中國工程物理研究院 總體工程研究所， 四川 綿陽 621000)

1 引言

在海洋和軍事工程中，常常會涉及到結構的出水和入水問題，許多學者做了大量的研究。譬如，陸宇等[1]采用LS-DYNA建立了水下垂直發射裝置頭罩外層內側內壓對頭罩開裂過程及沖擊力影響的有限元模型。黃鴻鑫等[2]通過結構入水的有限元模型分別研究了射彈的頭部形狀和質心位置對高速入水穩定性的影響。由于光滑粒子法(SPH)能夠模擬結構入水過程中自由液面的破碎和翻轉等問題，因此近些年來SPH算法逐步被廣泛應用于結構入水的數值模擬。SPH方法是由 Gingold 和 Monaghan[3]和Lucy[4]提出的一種粒子類方法。對于有限水域中結構入水的數值模擬，流域通常采用固壁邊界條件，但是當模擬的流域為無限域時，就會涉及到無反射邊界的應用。就目前看來，SPH方法的無反射邊界條件主要包括三大類：第一類為開放邊界條件[5-8]，第二類為海綿層邊界條件[9-13]，第三類為時間線插值邊界條件[14]。

開放邊界條件是指有進口和出口的固定區域，當粒子通過進口進入計算域時，則添加新的粒子，而當粒子通過出口離開計算域時，則將該粒子刪除。Morris等[5]采用周期性邊界條件模擬低雷諾數條件下的層流，但是周期性邊界允許粒子通過出口以后再經過入口進入計算域，這樣隨著時間的推進，擾動粒子的數量將會持續下降，會降低模擬的精度。為了提高模擬精度，Lastiwka等[6]在出口的上游設置流入區，在出口的下游設置流出區，并且根據實際的流速確定流入流出的粒子數量，這樣能有效地避免擾動粒子數量的下降，但是這種方法僅僅在流出域刪除了粒子，依然不能完全避免流出域產生的反射波對上游計算域的影響。Tafuni等[8]在計算域的上游或者下游定義了緩沖區，緩沖區內粒子的物理量等于其對應詭粒子的物理量。但是詭粒子的位置信息是緩沖區內粒子經過邊界對稱映射得到的，這樣就大大降低計算效率，尤其是邊界形狀比較復雜的時候。第二類無反射邊界條件是海綿層邊界條件。Lind[9]和Liu[10]在模擬波的傳播的時候采用了阻尼層的概念，即粒子的速度隨著與邊界的靠近而呈指數衰減，這樣就能消除一部分邊界上的反射波。Altomave等[11]在模擬長鳳波的產生和吸收的時候，同樣采用了阻尼層的概念，只是根據耗散海灘的需要，將粒子速度的衰減由指數型改為二次型。密度變化層的概念是由Gong等[12]提出的，通過在計算域外圍增加幾十層密度可變的粒子，粒子密度隨時間的變化率與其距邊界的距離相關，以此來消除固壁邊界的反射波。但是海綿層法在模擬無反射邊界時，由于需要在計算域的外面設置幾十層粒子，因此計算量較大，導致計算效率偏低，特別是對三維問題，計算量將迅速增加。最近， Pingping Wang[14]提出了時間域插值邊界條件。時間域插值邊界條件是一種基于特征線法的無反射邊界，通過在邊界設置3～4層的僅僅向下游傳遞干擾的邊界粒子來避免反射波的產生。邊界粒子的任一時刻的性質(速度、壓力等)，則是根據波傳播的速度，通過上游不同時刻的粒子采用拉格朗日插值來獲取的。但是時間插值類邊界條件需要存儲上一個或幾個時間步內的粒子信息，對計算內存要求較高。本文聯合SBFEM算法和SPH算法，提出了一種新的無反射邊界條件，即SBFE虛粒子邊界，該邊界具有計算精度高、效率快等優點。最后，我們采用2個結構入水算例驗證了該SBFE虛粒子邊界的準確性。需要特別說明的是，本文中的所有結果均采用FORTRAN 95自編程序計算所得。

2 SBFE虛粒子邊界

無反射邊界對結構入水、水下爆炸等領域的數值模擬具有很大的意義。如圖1所示，SBFE虛粒子邊界只需在SPH流體粒子的外層額外布置2～4層的SBFE虛粒子，虛粒子的層數應該滿足大于或者等于SPH粒子的光滑半徑。

圖1 SBFE虛粒子邊界示意圖Fig.1 SBFE virtual particle boundary

SBFE虛粒子的速度和加速度是通過核函數從附近的SPH粒子中插值獲得的，其插值公式為：

(1)

式(1)中：W表示光滑粒子法的核函數；N表示支持域內的粒子總數；下標i，j表示不同的粒子；m表示粒子質量；ρ表示粒子的密度。

進一步地，虛粒子的加速度值可以采用差分的方法來獲得，即：

(2)

式(2)中，Δt表示時間步長。

比例邊界單元示意圖如圖2所示。虛粒子之間可以構成一個個比例邊界單元，因此每一層邊界的無窮域動力剛度矩陣的剛度陣A和阻尼陣B能夠按照SBFEM的高頻近似方法獲得[15]。

需要特別說明的是，位于同一位置的SBFE虛粒子和比例邊界單元節點具有相同加速度值，所以比例邊界單元節點的加速度值可以通過同一位置的SBFE虛粒子的加速度值獲得。

進一步地就能求解每一個比例邊界單元節點的等效加速度，其表達式為：

(3)

式(3)中：a表示由各個虛粒子組成的加速度列向量；ns表示單元外法向的單位向量；N表示比例邊界單元的形函數。

圖2 比例邊界單元示意圖Fig.2 The diagram of SBFEM

采用高頻近似方法就能求解每層比例邊界單元各節點的壓力值，即：

(4)

其中，

(5)

式(5)中：p表示邊界節點壓力的列向量；p(i)(i=1，2，3，…，M)為輔助變量。

位于同一位置的SBFE虛粒子和比例邊界單元節點具有相同的壓力值，而位于比例邊界單元節點間的SBFE虛粒子的壓力，則可以采用線性插值的方式來獲得。如圖3所示，虛粒子m的壓力值可以通過比例邊界單元節點k和節點h的壓力值獲得，即：

(6)

圖3 SBFE虛粒子和比例邊界單元節點示意圖Fig.3 SBFE virtual particles and SBFE nodes

3 SBFE虛粒子邊界的計算效率和精度

我們驗證了SBFE虛粒子邊界的準確性和計算效率。算例中，水池的長0.2 m，寬0.1 m。楔形體的一條邊長為0.04 m，底升角為30°(如圖4所示)，以初速度3.15 m/s豎直落入水中，重力加速度為9.8 m/s2。另外設置一個同樣大小的水池，該水池的邊界采用固壁邊界條件。為了驗證SBFE虛粒子邊界的準確性，還需要一個長為0.4 m，寬為0.2 m的水池(見圖4(a))，以保證其在計算過程中，沒有反射波進入感興趣的區域。

圖4 二維楔形體入水模型示意圖Fig.4 The model of wedge water entry

表1表示不同邊界條件SPH方法的計算耗時。從表1中可以明顯的看出，SBFE虛粒子邊界能夠極大地減少計算域內的粒子數量，從而提高SPH方法模擬無窮域流體的計算效率。更重要的是， SBFE虛粒子邊界的計算耗時和固壁邊界的計算耗時幾乎相同。這表明，SBFE虛粒子邊界不會顯著增加計算耗時，具有計算效率高的優點。

表1 不同邊界條件的計算耗時

圖5表示水池底邊中點的壓力曲線。在固壁邊界條件下，我們能很明顯地發現壓力值存在2個波峰點，第1個波峰點是入射波和反射波疊加的峰值，其大小約為SBFE虛粒子邊界結果的2倍，而第2個波峰點則是反射波與底邊反射波，側邊反射波的多重疊加。從圖5中可以發現，無窮遠邊界中點的壓力值與SBFE虛粒子邊界的中點壓力值幾乎相同。這說明了SBFE虛粒子邊界能夠準確的模擬邊界的無反射特征。

圖5 水池底邊中點的壓力曲線Fig.5 The pressure of point 1

4 無窮水域楔形體入水數值模擬

楔形體幾何結構如圖6所示，楔形體的高度h為0.01 m，底升角α為20°，結構的密度為600.3 kg/m3，結構的質量為0.18 kg。結構分別以初速度100 m/s、80 m/s、60 m/s、40 m/s、20 m/s、5 m/s落入水中。水池的大小為0.2 m×0.2 m。其中一個水池的邊界采用SBFE虛粒子邊界來模擬無窮水域的情形，另一個水池的邊界則采用固壁邊界條件來模擬有限域的結構入水問題。SPH粒子的初始間距為1e-3 m，整個流體域被離散為40 401個粒子。固體結構的初始粒子間距同樣采用1e-3 m，整個固體域被離散為297個粒子。在數值模擬過程中，時間步長均為5e-7 s，總共模擬5e-3 s，模擬總步數為10 000步。

圖7表示楔形體結構以100 m/s速度入水時流體的壓力云圖。圖7中，左側壓力云圖表示SBFE虛粒子邊界的模擬結果，右側壓力云圖表示固壁邊界的模擬結果。

圖6 楔形體幾何結構圖Fig.6 The geometry of two-dimensional wedge

圖7 楔形體結構入水壓力云圖Fig.7 The pressure distribution of wedge water entry

從圖7中可以看出，當壓力波到達邊界時，能穿過SBFE虛粒子邊界傳向無窮遠，而固壁邊界則會將壓力波反射重新進入計算域，從而影響整個流域的壓力場。圖7表明SBFE虛粒子邊界能夠模擬壓力波的透射過程，從而實現準確高效模擬無窮域中結構入水問題。

進一步的，我們模擬了剛體結構以100 m/s的速度入水的速度衰減曲線，如圖8，結構入水深度曲線如圖9所示。從圖9中可以看出，在結構入水的初期，即入水砰擊階段，SBFE虛粒子邊界的模擬結果和固壁邊界的模擬結果幾乎相同。這是因為此時反射波仍然沒有進入計算域的緣故。但是，在結構入水模擬的后期，固壁邊界的速度衰減比SBFE虛粒子邊界的速度衰減更快，這是因為固壁邊界的反射波抵達了結構底部，對結構施加了一個額外阻力的緣故。

圖8 楔形體結構入水的速度衰減曲線Fig.8 The velocity of the wedge in water entry

圖9 楔形體結構的入水深度曲線Fig.9 The depth of the wedge in water entry

不同初速度的結構入水深度曲線如圖10。從圖10可以看出，隨著初速度的提升，結構入水深度迅速增加。表2展示了2種不同邊界條件下結構在t=5 ms時的入水深度。

圖10 不同初速度的結構入水深度曲線Fig.10 The depth of the wedge with different velocities in water entry

表2 結構不同入水初速度的入水深度Table 2 The depth error of the wedge with different velocities in water entry

在表2中，入水初速度越大，2種邊界條件下模擬結果的相對誤差越大。這表明對于無窮水域中結構高速入水問題，固壁邊界反射波對模擬結果的影響已經不能再忽略。這是因為入水初速度越大，流體域中壓力波的最大值越大，固壁邊界的反射波也會越強烈，這樣就會導致流體域內的壓力場變動越劇烈，從而對無窮域內結構入水的數值模擬結果產生不可忽略的影響。因此，在SPH方法中模擬無窮域中結構低速入水時，采用固壁邊界可以得到理想的結果，但是當結構入水速度提高時，則需要進一步考慮邊界反射波的影響。SBFE虛粒子邊界能很好地模擬邊界波的透射過程，從而能夠準確地模擬無窮域中結構入水的問題。

5 結論

1) 本文提出的SBFE虛粒子邊界能夠在SPH框架下準確地模擬邊界壓力波的透射過程。

2) 在結構入水砰擊階段，結構入水速度急劇減小。SBFE虛粒子邊界和固壁邊界的模擬結果幾乎相同，對于僅僅關注砰擊載荷的數值模擬，可以采用固壁邊界條件。

3) 采用SPH方法模擬無窮水域中的結構入水時，低速入水，可以采用固壁邊界條件，高速入水，采用無反射邊界。