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

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

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

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

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

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

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

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

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

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

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

11、x 軸水平向右為正， y 軸垂直向上為正。推板運(yùn)動(dòng)速度u a t w = sin ，位移s a t a w = ? cos + ，a為推板運(yùn)動(dòng)幅值， 表示推板運(yùn)動(dòng)頻率。水深hd 1.0 米，水池長度= 20 f l 米。4.1 弱非線性波的模擬當(dāng)取推板運(yùn)動(dòng)頻率 = 3.14秒-1，運(yùn)動(dòng)周期T = 2 = 2秒，推板運(yùn)動(dòng)幅值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 ?？梢耘袛嗌刹▽儆谟邢匏钊醴蔷€性波。對(duì)于該問題的計(jì)算，主要是驗(yàn)證SPH 方法模擬造波的精度特點(diǎn)。整個(gè)流體區(qū)域分布400*20=8000個(gè)粒子，每個(gè)粒子空間大小V = 2.5e ? 4 i 米2，每個(gè)粒子光滑長度= 0.05 i h 米，核函數(shù)采用三次樣條核函數(shù)，人工粘性系數(shù)= 0.06 ，時(shí)間步長t = 2.0e ? 4秒，步進(jìn)50000 步，在Core23.0 的PC 機(jī)上計(jì)算時(shí)間5.1 小時(shí)。圖3 給出SPH 方法在不同時(shí)刻推板附近區(qū)域的粒子分布結(jié)果。由圖中的結(jié)果可以看出，SPH 方法計(jì)算的自由表面非常光順，在推板

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 個(gè)粒子占據(jù)的空間， i h 表示第i 個(gè)粒子的光滑長度， 表示與核函數(shù)有關(guān)的光滑長度伸展系數(shù)，當(dāng)核函數(shù)取三次樣條函數(shù)時(shí) = 2.0。通過大量的計(jì)算表明，當(dāng)體積分?jǐn)?shù) 0.75 i f 就可以認(rèn)為，該粒子位于自由表面上。圖4 給出了BEM 方法與SPH 方法結(jié)果的比較。由于自由表面某些區(qū)域的粒子存在相互擠壓，根據(jù)以上判斷標(biāo)準(zhǔn)對(duì)少

14、數(shù)區(qū)域不能準(zhǔn)確分辨存在空缺，但是從總體效果來看，該方法還是能夠反映出整個(gè)波面形狀的發(fā)展趨勢。由自由表面形狀的比較結(jié)果可以發(fā)現(xiàn)，在模擬線性規(guī)則波時(shí)SPH 與BEM 方法的計(jì)算結(jié)果吻合得較好，體現(xiàn)了SPH 方法在弱非線性自由表面波的模擬上具有較高的精度。4.2 強(qiáng)非線性波的模擬當(dāng)增加推板運(yùn)動(dòng)頻率 = 6.28秒-1，推板運(yùn)動(dòng)幅值保持a = 0.1米，運(yùn)動(dòng)周期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 ?？梢猿醪焦烙?jì)生成波屬于深水強(qiáng)非線性波，生成的波浪會(huì)發(fā)生翻轉(zhuǎn)破碎現(xiàn)象。對(duì)于該問題的計(jì)算，主要是驗(yàn)證SPH 方法在模擬強(qiáng)非線性波破碎波的精度特點(diǎn)?？紤]到模擬波長減小，將計(jì)算區(qū)域長度取為=10 f l 米，既保證在模擬時(shí)間內(nèi)波浪不會(huì)發(fā)生反射，又可以減少粒子數(shù)目，減少計(jì)算量。當(dāng)整個(gè)流體區(qū)域分布200*20=4000 個(gè)粒子，其它因素保持不變的情況下，步進(jìn)25000 步計(jì)算耗時(shí)2.5 小時(shí)。圖5 給出SPH 方法從初始狀態(tài)到3.0 秒時(shí)，每隔0.5 秒推板附近區(qū)域粒子分布結(jié)果。由圖中結(jié)果可以看出，SPH 方法計(jì)算出的波高度大約為0.4 米，波長在1.5 米至1.8 米之間。波面形狀并不光順，但是自

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

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

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