基于SBFE虛粒子邊界的結(jié)構(gòu)入水?dāng)?shù)值模擬

李上明，趙梓斌

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

1 引言

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

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

2 SBFE虛粒子邊界

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

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

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

(1)

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

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

(2)

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

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

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

進(jìn)一步地就能求解每一個(gè)比例邊界單元節(jié)點(diǎn)的等效加速度，其表達(dá)式為：

(3)

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

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

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

(4)

其中，

(5)

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

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

(6)

圖3 SBFE虛粒子和比例邊界單元節(jié)點(diǎn)示意圖Fig.3 SBFE virtual particles and SBFE nodes

3 SBFE虛粒子邊界的計(jì)算效率和精度

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

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

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

表1 不同邊界條件的計(jì)算耗時(shí)

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

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

4 無窮水域楔形體入水?dāng)?shù)值模擬

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

圖7表示楔形體結(jié)構(gòu)以100 m/s速度入水時(shí)流體的壓力云圖。圖7中，左側(cè)壓力云圖表示SBFE虛粒子邊界的模擬結(jié)果，右側(cè)壓力云圖表示固壁邊界的模擬結(jié)果。

圖6 楔形體幾何結(jié)構(gòu)圖Fig.6 The geometry of two-dimensional wedge

圖7 楔形體結(jié)構(gòu)入水壓力云圖Fig.7 The pressure distribution of wedge water entry

從圖7中可以看出，當(dāng)壓力波到達(dá)邊界時(shí)，能穿過SBFE虛粒子邊界傳向無窮遠(yuǎn)，而固壁邊界則會(huì)將壓力波反射重新進(jìn)入計(jì)算域，從而影響整個(gè)流域的壓力場。圖7表明SBFE虛粒子邊界能夠模擬壓力波的透射過程，從而實(shí)現(xiàn)準(zhǔn)確高效模擬無窮域中結(jié)構(gòu)入水問題。

進(jìn)一步的，我們模擬了剛體結(jié)構(gòu)以100 m/s的速度入水的速度衰減曲線，如圖8，結(jié)構(gòu)入水深度曲線如圖9所示。從圖9中可以看出，在結(jié)構(gòu)入水的初期，即入水砰擊階段，SBFE虛粒子邊界的模擬結(jié)果和固壁邊界的模擬結(jié)果幾乎相同。這是因?yàn)榇藭r(shí)反射波仍然沒有進(jìn)入計(jì)算域的緣故。但是，在結(jié)構(gòu)入水模擬的后期，固壁邊界的速度衰減比SBFE虛粒子邊界的速度衰減更快，這是因?yàn)楣瘫谶吔绲姆瓷洳ǖ诌_(dá)了結(jié)構(gòu)底部，對結(jié)構(gòu)施加了一個(gè)額外阻力的緣故。

圖8 楔形體結(jié)構(gòu)入水的速度衰減曲線Fig.8 The velocity of the wedge in water entry

圖9 楔形體結(jié)構(gòu)的入水深度曲線Fig.9 The depth of the wedge in water entry

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

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

表2 結(jié)構(gòu)不同入水初速度的入水深度Table 2 The depth error of the wedge with different velocities in water entry

在表2中，入水初速度越大，2種邊界條件下模擬結(jié)果的相對誤差越大。這表明對于無窮水域中結(jié)構(gòu)高速入水問題，固壁邊界反射波對模擬結(jié)果的影響已經(jīng)不能再忽略。這是因?yàn)槿胨跛俣仍酱?，流體域中壓力波的最大值越大，固壁邊界的反射波也會(huì)越強(qiáng)烈，這樣就會(huì)導(dǎo)致流體域內(nèi)的壓力場變動(dòng)越劇烈，從而對無窮域內(nèi)結(jié)構(gòu)入水的數(shù)值模擬結(jié)果產(chǎn)生不可忽略的影響。因此，在SPH方法中模擬無窮域中結(jié)構(gòu)低速入水時(shí)，采用固壁邊界可以得到理想的結(jié)果，但是當(dāng)結(jié)構(gòu)入水速度提高時(shí)，則需要進(jìn)一步考慮邊界反射波的影響。SBFE虛粒子邊界能很好地模擬邊界波的透射過程，從而能夠準(zhǔn)確地模擬無窮域中結(jié)構(gòu)入水的問題。

5 結(jié)論

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

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

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