推板造波的光滑粒子流體動力學(xué)數(shù)值模擬探討

1、推板造波的光滑粒子流體動力學(xué)數(shù)值模擬探討    1 前言造波問題是船舶工程和海洋工程領(lǐng)域十分普遍的現(xiàn)象，是模擬船舶與海洋結(jié)構(gòu)物相互作用波浪環(huán)境的基礎(chǔ)。過去的幾十年中，造波問題的數(shù)值方法不斷得到更新和發(fā)展。最早對造波問題進行研究是基于勢流理論的邊界元方法(BEM)。勢流理論假設(shè)流體沒有旋度，沒有粘性，不可壓縮，流體的速度勢 滿足拉普拉斯方程。Longuet-Higgins 和Cokelet最早采用邊界元方法對二維水波問題進行研究。Baudic et al對于非線性水波的產(chǎn)生，傳播和吸收都做了詳細(xì)研究。賀五州和段文洋采用完全非線性邊界元方法對二維搖板造波問題

2、進行了時域計算。該類方法的優(yōu)點可以自然滿足質(zhì)量守恒，不存在能量耗散，計算量較小，但是對粘性效應(yīng)，自由表面發(fā)生翻轉(zhuǎn)、破碎現(xiàn)象還無法模擬。隨著計算機運算能力的不斷提高，很多學(xué)者發(fā)展和應(yīng)用基于直接求解N-S 或者雷諾平均N-S 方程的方法。該類方法最具挑戰(zhàn)性的問題是如何處理非穩(wěn)態(tài)非線性自由表面。目前，常用的方法是VOF 方法和Level Set方法，該類方法都是通過求解每個網(wǎng)格上的一個標(biāo)量函數(shù)來描述界面的幾何特征，標(biāo)量函數(shù)的輸運方程形式上完全相同，同屬于界面撲獲類算法，都有較強的拓?fù)浔磉_能力。但是目前該類數(shù)值方法對于自由表面，特別是伴隨有翻卷、破碎等現(xiàn)象的復(fù)雜自由表面問題，所取得的成就還很有限。目前

3、，一種被稱為無網(wǎng)格的方法引起人們關(guān)注。因為在處理大變形的情況下，不需要重新劃分網(wǎng)格。該類方法對于復(fù)雜自由表面問題具有較大的潛力。本文所采用的光滑粒子流體動力學(xué)(SPH)方法是無網(wǎng)格方法的一種，該方法的計算建立在一群離散的粒子上，函數(shù)以及導(dǎo)數(shù)的求解可以利用周圍粒子的相關(guān)值積分得到。某個粒子影響區(qū)域的大小稱為支持域。支持域內(nèi)越靠近中心粒子影響越明顯，支持域外粒子對計算結(jié)果沒有影響。SPH 模擬方法最早用于天體物理學(xué)，并由Monaghan將該方法引入自由表面流動問題。由于SPH 方法中是對每個粒子的運動軌跡進行模擬，這些粒子能夠自動滿足自由表面的運動學(xué)條件。目前，SPH 方法在自由表面流動問題上已經(jīng)

4、取得一些成就，例如Y.M.Lo& Shao用SPH 方法結(jié)合大渦模擬方法對近岸孤立波的爬升進行模擬。Colagrossi &Landrini采用SPH 方法對二維兩相流問題進行計算分析。Souto et al采用SPH 方法對液艙晃蕩的相位滯后問題進行研究，之后又實現(xiàn)了SPH 方法對液艙晃蕩橫搖力矩的模擬。國內(nèi)對SPH 方法的研究已經(jīng)開始，但是還處于學(xué)習(xí)和消化階段。雖然SPH 方法具有巨大潛力，但是該方法同時還存在一些問題有待解決，例如計算效率較低，數(shù)值結(jié)果不很穩(wěn)定，但是該方法編寫程序較容易，對于大變形問題和動邊界問題的處理方便。因此，本文嘗試采用SPH方法對二維推板造波問題進

5、行模擬，通過與全非線性邊界元和VOF 方法比較，對SPH 方法模擬造波問題的適用性進行研究。本文的主要內(nèi)容包括四個方面。首先對SPH 方法的基本原理進行簡單介紹，給出SPH方法模擬自由表面流動問題的控制方程，然后對求解控制方法的若干技巧進行了介紹，包括壓力與密度的顯式關(guān)系，人工粘性，邊界條件處理方法。然后給出二維推板造波問題計算結(jié)果，通過與全非線性邊界元和VOF 方法比較，得到了SPH 方法模擬該問題的精度特點，特別是當(dāng)波面出現(xiàn)翻轉(zhuǎn)現(xiàn)象時，SPH 方法都能夠準(zhǔn)確的反映出來。最后對SPH 方法計算造波問題的精度進行了總結(jié)，對存在的不足提出改進方法。2 基本原理2.1SPH 方法SPH 方法的基本

6、理論包括兩部分：核近似和粒子近似。核近似的思想來自函數(shù)f (x)的積分表達式(1)式中為包含x的積分域。h 是核函數(shù)的光滑長度，由于核函數(shù)W 不是 函數(shù)，故(3)式只是近似式，(3)式這種表示函數(shù)的近似方法稱為核近似(Kernel approximation)。核函數(shù)有多種形式，其中三次樣條核函數(shù)，高斯核函數(shù)，五次樣條核函數(shù)應(yīng)用得最廣泛。所有核函數(shù)需要滿足一些特性，如正交性，緊支性，非負(fù)性，單調(diào)性等等。求解流體動力學(xué)方程的關(guān)鍵除了求函數(shù)本身，更重要的是如何用核近似來計算函數(shù)導(dǎo)數(shù)。根據(jù)流體系統(tǒng)是由具有獨立質(zhì)量，占有獨立空間的有限個粒子組成的假設(shè)，可以將連續(xù)形式的函數(shù)核近似用粒子近似表達。2.2

7、SPH 形式的流體動力學(xué)方程水波問需要求解的方程包括質(zhì)量守恒方程和動量守恒方程。其中i F 表示質(zhì)量力，在水波問題求解時i F 表示重力項， 表示某一矢量方向。接下來對求解上述方程的一些技巧進行介紹。3 數(shù)值解法若干技巧3.1 人工可壓縮性在SPH 方法中認(rèn)為流體系統(tǒng)具有弱可壓縮性。流場密度和壓力的關(guān)系可以通過顯式關(guān)系直接表示出來，這種壓力的確定方法就稱為人工可壓縮性。P = C ，對于水 = 7， 0 C 為聲速，可取流體區(qū)域中粒子最大速度的10 倍。這是目前SPH 方法模擬水波問題所普遍采用的形式。3.2 人工粘性在SPH 方法的計算中，由于粒子容易發(fā)生穿透，重疊現(xiàn)象，使得核近似計算結(jié)果不

8、穩(wěn)定，為了避免以上現(xiàn)象發(fā)生，Monaghan在動量方程中加入人工粘性項，主要作用可以避免粒子相互穿透，使得整個區(qū)域的粒子排列更加有序，可以較好的改善計算穩(wěn)定性。3.3 邊界處理當(dāng)粒子位于計算區(qū)域內(nèi)部時稱為自由粒子。當(dāng)粒子的支持域沒有被邊界截斷， (3)式和(5)式成立。當(dāng)粒子位于邊界區(qū)域附近時，粒子支持域被邊界截斷，此時 (3)式和(5)式不再成立，邊界附近處就要采用特殊的方法進行處理。Monaghan3提出過固壁邊界的處理辦法。為了保證固壁邊界不可穿透，在邊界處安放一些固壁邊界粒子，當(dāng)自由粒子靠近固壁粒子時，固壁粒子將施加外力使自由粒子遠(yuǎn)離固壁邊界來保證自由粒子不穿透固壁。但該方法一方面缺乏

9、物理意義，另一方面當(dāng)自由粒子速度較高時，固壁邊界往往會施加過大的外力而使整個粒子系統(tǒng)不穩(wěn)定。本文采用一種鏡像粒子方法。在邊界外側(cè)生成一些鏡像粒子，鏡像粒子與對應(yīng)自由粒子到邊界的距離相同，鏡像粒子與所對應(yīng)的自由粒子的速度切向相同，法向相反，如圖1 所示。這樣可以保證邊界處的自由粒子沒有密度缺損，而且可以保證自由粒子和鏡像粒子的速度累加效應(yīng)在固壁邊界法向上為零，滿足固壁邊界不可穿透條件。3.4 步進方法在解問題(8)(11)時，通常采用的步進方法主要有三種，改進Euler 法，蛙跳 (leap-frog)法和四階的龍格-庫塔法。但是通過大量的計算發(fā)現(xiàn)，這些顯式的步進格式對計算結(jié)果影響并不明顯，因此

10、本文采用了較簡單的改進Euler 法，其取值與粒子數(shù)和聲速有關(guān)，一般粒子數(shù)目增大或聲速增大，時間步長t取值需減小。4 數(shù)值結(jié)果與分析推板造波問題是研究波浪與結(jié)構(gòu)物相互作用的典型例子。目前對于該方面的研究通常還只限于小波幅非線性波的模擬，但是對于大幅強非線性波模擬，特別是伴隨有翻轉(zhuǎn)、破碎現(xiàn)象的造波問題，傳統(tǒng)的數(shù)值方法還面臨極大的挑戰(zhàn)。本文所計算的推板造波模型如圖2 所示，在平坦池底的左側(cè)安放有造波機，右側(cè)為一垂直固壁。由于本文主要模擬SPH 方法短期造波過程，假設(shè)波浪沒有到達右側(cè)墻壁，不存在波浪反射問題。有關(guān)避免波浪反射，增加阻尼區(qū)的研究正在進行中。坐標(biāo)原點建立在初始時刻造波機與靜止水面的交點，

11、x 軸水平向右為正， y 軸垂直向上為正。推板運動速度u a t w = sin ，位移s a t a w = ? cos + ，a為推板運動幅值， 表示推板運動頻率。水深hd 1.0 米，水池長度= 20 f l 米。4.1 弱非線性波的模擬當(dāng)取推板運動頻率 = 3.14秒-1，運動周期T = 2 = 2秒，推板運動幅值a = 0.1米，由有限水深的線性波色散關(guān)系可以估算所造波數(shù)k = 2 =1.01米-1，波長 = 6.25米，并由傳遞函數(shù)21 可得生成的波高4 sinh2 /(2 sinh 2 ) 0.0990 = + = d d d a a kh kh kh 米， 波陡0.032 1

12、20 0 a = < ， = 0.16 < 1 2 d h 。可以判斷生成波屬于有限水深弱非線性波。對于該問題的計算，主要是驗證SPH 方法模擬造波的精度特點。整個流體區(qū)域分布400*20=8000個粒子，每個粒子空間大小V = 2.5e ? 4 i 米2，每個粒子光滑長度= 0.05 i h 米，核函數(shù)采用三次樣條核函數(shù)，人工粘性系數(shù)= 0.06 ，時間步長t = 2.0e ? 4秒，步進50000 步，在Core23.0 的PC 機上計算時間5.1 小時。圖3 給出SPH 方法在不同時刻推板附近區(qū)域的粒子分布結(jié)果。由圖中的結(jié)果可以看出，SPH 方法計算的自由表面非常光順，在推板

13、和底部邊界附近沒有出現(xiàn)粒子發(fā)散和穿透現(xiàn)象，表明鏡像粒子的邊界處理辦法是穩(wěn)定可靠的。其中i f 稱為體積分?jǐn)?shù)，表示i 粒子影響域內(nèi)粒子空間與影響域面積的比值， i f 越大說明該粒子周圍的粒子數(shù)越多， i f 越小說明該粒子周圍的粒子數(shù)越少。j V 表示i 粒子支持域內(nèi)第j 個粒子占據(jù)的空間， i h 表示第i 個粒子的光滑長度， 表示與核函數(shù)有關(guān)的光滑長度伸展系數(shù)，當(dāng)核函數(shù)取三次樣條函數(shù)時 = 2.0。通過大量的計算表明，當(dāng)體積分?jǐn)?shù) 0.75 i f 就可以認(rèn)為，該粒子位于自由表面上。圖4 給出了BEM 方法與SPH 方法結(jié)果的比較。由于自由表面某些區(qū)域的粒子存在相互擠壓，根據(jù)以上判斷標(biāo)準(zhǔn)對少

14、數(shù)區(qū)域不能準(zhǔn)確分辨存在空缺，但是從總體效果來看，該方法還是能夠反映出整個波面形狀的發(fā)展趨勢。由自由表面形狀的比較結(jié)果可以發(fā)現(xiàn)，在模擬線性規(guī)則波時SPH 與BEM 方法的計算結(jié)果吻合得較好，體現(xiàn)了SPH 方法在弱非線性自由表面波的模擬上具有較高的精度。4.2 強非線性波的模擬當(dāng)增加推板運動頻率 = 6.28秒-1，推板運動幅值保持a = 0.1米，運動周期T =1秒，由有限水深的線性波色散關(guān)系可以估算所造波長 = 1.56米，波數(shù)k = 4.02米-1，同樣由傳遞函數(shù)可得生成波幅0.199 0 a = 米，波陡2 0.255 1 7 0 a = > ， = 0.64 > 1 2 d

15、h ?？梢猿醪焦烙嬌刹▽儆谏钏畯姺蔷€性波，生成的波浪會發(fā)生翻轉(zhuǎn)破碎現(xiàn)象。對于該問題的計算，主要是驗證SPH 方法在模擬強非線性波破碎波的精度特點?？紤]到模擬波長減小，將計算區(qū)域長度取為=10 f l 米，既保證在模擬時間內(nèi)波浪不會發(fā)生反射，又可以減少粒子數(shù)目，減少計算量。當(dāng)整個流體區(qū)域分布200*20=4000 個粒子，其它因素保持不變的情況下，步進25000 步計算耗時2.5 小時。圖5 給出SPH 方法從初始狀態(tài)到3.0 秒時，每隔0.5 秒推板附近區(qū)域粒子分布結(jié)果。由圖中結(jié)果可以看出，SPH 方法計算出的波高度大約為0.4 米，波長在1.5 米至1.8 米之間。波面形狀并不光順，但是自

16、由表面發(fā)生翻轉(zhuǎn)破碎的現(xiàn)象還不清晰?？梢钥闯鯯PH 方法模擬的最大波高處往往不連續(xù)，發(fā)生了比較模糊的破碎現(xiàn)象，而BEM 方法自由表面始終保持光順沒有發(fā)生破碎現(xiàn)象。在推板附近波浪形成階段，兩種方法計算的自由表面形狀符合得較好，當(dāng)波浪傳播大約一個波長后，SPH 模擬波面與BEM 計算結(jié)果存在相位差，而且隨著波浪傳播距離增大相位差也增大。這是由于SPH 方法模擬這種較陡峭的波峰會發(fā)生破碎現(xiàn)象，使得生成的波長增大波速增加的結(jié)果。為了能夠更加清晰的觀察到翻轉(zhuǎn)破碎現(xiàn)象，將粒子總數(shù)增加到400*40=16000，在1.5 秒至1.8 秒時刻對推板附近區(qū)域進行放大。當(dāng)粒子數(shù)增加后在1.5 秒時，SPH 方法計算

17、的波形開始發(fā)生翻轉(zhuǎn)，在1.6 秒和1.7 秒可以清晰的觀察到波面翻轉(zhuǎn)的發(fā)展過程，在1.8 秒時翻轉(zhuǎn)現(xiàn)象結(jié)束。BEM 方法在所對應(yīng)時刻還無法對該現(xiàn)象進行清晰的捕捉。因此，為了更好的對于波浪細(xì)節(jié)進行描述，必須使該區(qū)域內(nèi)的粒子數(shù)目增加到一定范圍。本文推薦能夠?qū)Σɡ朔D(zhuǎn)進行清晰描述的粒子半徑/8 0 r a i ，其中0 a 為浪高。為了對該頻率和幅值下波浪發(fā)生的翻轉(zhuǎn)現(xiàn)象進行更加充分驗證，本文還采用FLUENT軟件對該問題進行了計算。FLUENT 軟件對微分方程的求解采用有限體積法，自由表面捕捉采用VOF 方法。左側(cè)推板采用動網(wǎng)格方法，推板運動方式編寫UDF 程序加以控制，運動方式與SPH 方法完全相

18、同。圖8 給出了不同時刻兩種方法的自由表面的比較結(jié)果。其中SPH方法的粒子數(shù)為16,000，VOF 方法的網(wǎng)格數(shù)為200,000。圖中可以看出在各時刻，SPH 與VOF 方法的計算結(jié)果總體上都吻合的很好，波面形狀保持了較好的一致性。只是在自由表面發(fā)生翻轉(zhuǎn)破碎后，復(fù)雜的自由表面形狀還有待實驗的驗證，例如射流發(fā)生的時間和形狀，氣穴生成的大小和位置等等。5 結(jié)論與展望本文采用SPH 方法對推板造波問題進行模擬研究。當(dāng)所造波是線性波時，由于BEM 方法已經(jīng)能夠?qū)υ搯栴}進行高效、準(zhǔn)確的計算，因此SPH 與該方法完全吻合的計算結(jié)果，可以充分證明SPH 在模擬線性波時具有較高的精度。當(dāng)模擬非線性波，特別是破碎波的發(fā)展過程時，BEM 方法就顯得力不從心，但是SPH 卻能非常方便的實現(xiàn)波浪翻轉(zhuǎn)，破碎，重入等復(fù)雜的自由表面現(xiàn)象。通過與VOF 方法的比較，兩種方法所計算出的波形在不同時刻都具有較好的一致性，只是波浪破碎后所產(chǎn)生的射流，