黃曉婷, 孫鵬楠, 呂鴻冠, 殷曉瑞, 董嘉徐

(1.中山大學(xué) 海洋工程與技術(shù)學(xué)院, 廣東 珠海 519082;2.南方海洋科學(xué)與工程廣東省實驗室(珠海), 廣東 珠海 519082)

近年來,微型飛行器或潛水器的發(fā)展受到廣泛關(guān)注,飛行器或潛水器性能的改進離不開對翼型繞流問題的精確模擬。翼型表面流動如分離、附著等現(xiàn)象對翼型的動力性能具有重要影響。為保證在不同工況下飛行器或潛水器運行的可靠性,需要了解翼型繞流的流動特性和機理。這就要求對翼型繞流現(xiàn)象進行深入研究和分析。

翼型的幾何形狀對改善力學(xué)性能具有重要影響[1]?？梢园l(fā)現(xiàn),目前對于非對稱翼型的繞流現(xiàn)象的研究較少。其中,Anyoji等[2]經(jīng)過研究發(fā)現(xiàn),非對稱Ishii弧形翼型由于后緣附近的弧度,其壓力阻力要遠遠小于其他對稱翼型,在同等工況下,可獲得比對稱翼型更高的升力。因此,展開對Ishii非對稱翼型更深入研究,準確獲取流場信息及其運動特性,可為未來設(shè)計出良好動力性能的新型翼型提供新思路。

現(xiàn)有翼型繞流研究中,主流研究方法中主要有試驗法及數(shù)值模擬法。針對水翼繞流,季斌等[3]從試驗和數(shù)值模擬角度綜述水翼水動力特性。在飛行氣翼繞流探究中,Anyoji等[2]采用風(fēng)洞試驗探究及驗證Ishii翼型氣動性能。然而在試驗法中,試驗條件復(fù)現(xiàn)存在一定挑戰(zhàn),且邊界層湍流細節(jié)捕捉對儀器的精密性要求苛刻。此外,試驗法的成本高昂,且對多工況無法同時進行研究。因此現(xiàn)有研究大多采取數(shù)值模擬法展開研究。數(shù)值模擬法中，發(fā)展較為成熟的歐拉網(wǎng)格法較為常見。但對于網(wǎng)格法而言,翼型后緣狹窄尖角邊界層較薄,網(wǎng)格劃分邊界層存在一定難度[4],在網(wǎng)格離散中需要較高的網(wǎng)格構(gòu)建技術(shù)來滿足復(fù)雜的邊界條件。值得注意的是,近年來發(fā)展的無網(wǎng)格法對邊界條件、自由液面及大變形等問題[5]處理具有固有優(yōu)勢。

SPH方法作為無網(wǎng)格法中應(yīng)用廣泛的方法,在流體或固體力學(xué)問題處理中逐漸發(fā)揮重要作用[6]且對于自由液面變形動邊界問題的處理上SPH方法具有顯著優(yōu)勢[7]。目前,在水動力學(xué)領(lǐng)域上得到較多應(yīng)用[8-9]。Sun等[10]采用δ+-SPH模型對NACA0012翼型固定及自由運動狀態(tài)進行模擬,分析其水動力特性。Huang等[11]采用SPH方法對翼型繞流的運動性能進行探究;在空氣動力學(xué)領(lǐng)域,Shadloo等[12]驗證了WCSPH可以實現(xiàn)低雷諾數(shù)下不同攻角翼型繞流計算;Huang等[13]采用KGF-SPH模型,在低雷諾數(shù)下對氣翼的升阻力系數(shù)進行監(jiān)測,驗證該模型的準確性?？傻贸?SPH方法在翼型繞流問題研究上具有一定的優(yōu)勢及特色。但由于其本身的精度及穩(wěn)定性問題,SPH方法在翼型氣動性上的研究受到一定限制,翼型繞流中邊界層內(nèi)部流場特征獲取仍需進一步的方法改進和技術(shù)討論。

因此,本文克服以往SPH方法難以模擬黏性邊界層流動的不足,采用δ+-SPH模型,實現(xiàn)黏性流動下翼型繞流的模擬,展開翼型特性分析。在SPH計算中,高雷諾數(shù)下翼型繞流模擬存在三大難點,主要應(yīng)用以下方法進行解決:①針對流體繞過翼型時,由于負壓產(chǎn)生的張力不穩(wěn)定現(xiàn)象,采用張力不穩(wěn)定性控制(tensile instability control,TIC)技術(shù),實現(xiàn)Ishii翼型黏性繞流的壓力場和速度場等穩(wěn)定數(shù)值模擬。②為了避免粒子的拉格朗日結(jié)構(gòu)引起的計算精度下降現(xiàn)象,在每個時間步初始時刻,施加粒子修正技術(shù)(particle shifting technique,PST),確保粒子的均勻分布。③為實現(xiàn)邊界層黏性力的精確計算,采用自適應(yīng)粒子細化(adaptive particle refinement,APR)技術(shù),減少粒子數(shù)量,提高計算效率,實現(xiàn)粒子高分辨率同時確保計算精度。

1 δ+-SPH模型

在δ+-SPH模型中,考慮流體的真實物理黏性,在動量方程中施加黏性項,以精確模擬邊界層的流動。再有,在連續(xù)方程中添加密度耗散項[14],以避免高頻的壓力振蕩。此外,為了避免傳統(tǒng)SPH模型中由于粒子的拉格朗日結(jié)構(gòu)引起的計算精度降低現(xiàn)象,δ+-SPH模型針對粒子進行位移修正,引入粒子位移修正技術(shù)(particle shifting technique,PST),確保在計算過程粒子分布的均勻性。因此,基于流體弱可壓假設(shè)下,δ+-SPH模型[15]最終寫為

(1)

式中：i和j表示支持域內(nèi)目標粒子及其對應(yīng)的鄰近粒子；fi是質(zhì)量力,本文計算均取fi=0；ρi,mi,Vi,ui,ri分別代表粒子i的密度、質(zhì)量、體積、速度和位移；光滑核函數(shù)Wij采用Wendland的C2核函數(shù)[16],光滑長度h=2Δx,其中Δx是初始粒子間距；耗散項Di和πi的具體參數(shù)可參考文獻[15]。Fij為壓力項,在本文中采用張力不穩(wěn)定控制技術(shù)進行處理。

為確保粒子的均勻分布,在δ+-SPH模型中,動量方程計算得出的粒子位移ri需要施加人工位移修正量δri,δri具體定義為

(2)

2 張力不穩(wěn)定性控制技術(shù)

在傳統(tǒng)SPH模型中,由于負壓會導(dǎo)致張力不穩(wěn)定現(xiàn)象[18],流場信息將無法很好呈現(xiàn)。基于Sun等[15]的研究結(jié)果,為更好實現(xiàn)漩渦附著、分離、脫落細節(jié)的捕捉,展開對Ishii翼型流場規(guī)律分析,本文采用TIC技術(shù)改進壓力梯度項。對壓力項Fij進行改進[15]，如(3)式所示

(3)

DF表示液面區(qū)域自由液面粒子及相鄰粒子的集合。

3 拉格朗日擬序結(jié)構(gòu)可視化

在歐拉算法的模擬結(jié)果中,流場特征量(如渦量)通?；跉W拉框架下的瞬態(tài)速度場計算得出,難以給出流體介質(zhì)在某個時間段內(nèi)的輸運特征,而本文采用拉格朗日擬序結(jié)構(gòu)對流場中漩渦結(jié)構(gòu)進行可視化,可更加清晰地獲取流場細節(jié)。拉格朗日擬序結(jié)構(gòu)通過計算有限時間李亞普諾夫指數(shù)(finite-time Lyapunov exponent,FTLE)實現(xiàn)。在FTLE云圖中呈現(xiàn)的脊線為拉格朗日擬序結(jié)構(gòu)[19]。

在SPH計算程序中采用逆向FTLE展開對FTLE計算,以便能清晰反映漩渦結(jié)構(gòu)邊界特征[20]。逆向FTLE計算方程為

(4)

4 多級粒子分辨率技術(shù)

翼型的首尾部往往較為狹窄,翼型繞流中表面邊界層內(nèi)湍流流動的捕捉對于網(wǎng)格法和無網(wǎng)格法均具有一定挑戰(zhàn)性。因此,本文將利用多級粒子分辨率技術(shù)[21],將翼型表面離散成分辨率高的固定虛粒子,實現(xiàn)精確捕捉湍流流動細節(jié),提高邊界黏性力的計算精度。多級粒子分辨率技術(shù)原理如圖1所示。

圖1 多級粒子分辨率技術(shù)原理示意圖[21]

主要原理闡述如下:不同分辨率的粒子定義為不同的粒子集合。首先定義翼型附近區(qū)域為細化區(qū),在細化區(qū)邊緣定義過渡區(qū),如圖1b)所示,充當(dāng)其邊界條件,獲取進入細化區(qū)域母粒子的信息。過渡區(qū)域的寬度設(shè)置要大于核函數(shù)的半徑,防止其邊緣處核函數(shù)截斷而造成計算誤差。在計算中,參與流場SPH計算的粒子稱為SPH粒子;活性關(guān)閉,不參與SPH計算的粒子為非活性粒子。在細化區(qū),母粒子分裂,如圖1a)所示,在正方形4個頂點處形成更小粒子即子粒子,子粒子活性激活,參與流場的SPH計算,即為SPH粒子。而母粒子在撕裂之后,作為非活性粒子不再參與SPH計算,但仍然隨著流場運動,通過從周圍細化的子粒子插值獲得流場信息如速度、壓力等。相反的,在過渡區(qū)域,母粒子為SPH粒子,子粒子屬性從周圍SPH粒子即母粒子插值得出。采用Shepard插值方法,非活性粒子的物理量φ從周圍SPH粒子完成插值,插值公式為

(5)

經(jīng)過此過程,2層分辨率粒子完成了相互之間的信息交換。多級粒子分辨率技術(shù)對局部流場粒子加密,在不同的分辨率區(qū)域之間,過渡區(qū)域及細化區(qū)域并不隨著粒子的移動而移動。在計算中,當(dāng)母粒子撕裂而來的子粒子流出細化區(qū)及過渡區(qū)域時,不發(fā)生粒子合并而將直接被刪除,而母粒子流出細化區(qū)域時,將正常參與SPH流場計算?；谝陨狭Ｗ蛹毣枷?在多層分辨率粒子間,保持粒子尺度一樣,每層粒子的過渡區(qū)域流場信息從下一層插值取得,從而實現(xiàn)不同層次間粒子信息交換,實現(xiàn)多級分辨率的SPH數(shù)值計算,提高流場計算精度。在下文將首先通過圓柱繞流基準算例驗證在δ+-SPH中添加多級粒子分辨率技術(shù)的有效性,從而應(yīng)用于翼型繞流問題的研究中。

5 基準算例驗證

圓柱繞流問題是流體力學(xué)中的經(jīng)典問題,是檢驗數(shù)值算法的經(jīng)典算例。選取雷諾數(shù)Re=3 000,人工聲速c0=15 m/s時的圓柱繞流模擬算例,與試驗結(jié)果對比以檢驗δ+-SPH模型的數(shù)值精度。

在該算例中,施加多級粒子分辨率技術(shù),采用5層分辨率不同的粒子。圓柱表面的粒子間距為D/Δx=400,最遠處粒子間距為D/Δx=25。在計算中,粒子總數(shù)為1.45×106。假設(shè)在計算域中粒子分辨率均取為D/Δx=400,粒子總數(shù)將達到7.68×107。可見,多級粒子分辨率技術(shù)對重點區(qū)域流場進行加密,在保證計算精度條件下,有效減少了粒子總數(shù),從而減少計算時間,提高計算效率,對實現(xiàn)SPH方法在工程上應(yīng)用具有重要意義。

圖2 雷諾數(shù)Re=9 500圓柱繞流流場特征

圖2中選取圓柱上表面的δ+-SPH計算圖展示圓柱繞流的漩渦發(fā)展過程。通過漩渦流線圖與試驗結(jié)果[22]對比,可以發(fā)現(xiàn)試驗結(jié)果與δ+-SPH結(jié)果吻合良好。δ+-SPH模型與多級粒子分辨率技術(shù)相結(jié)合,漩渦發(fā)展細節(jié)能較好被模擬,δ+-SPH結(jié)合多級粒子分辨率技術(shù)可以較為準確模擬繞流現(xiàn)象中分離、附著現(xiàn)象。SPH中拉格朗日擬序結(jié)構(gòu)能清晰展示漩渦演化過程,呈現(xiàn)清晰的內(nèi)部細節(jié),與試驗結(jié)果基本保持一致,說明拉格朗日擬序結(jié)構(gòu)獲取流場特征的高精度?；谝陨嫌懻?可見SPH方法在漩渦結(jié)構(gòu)捕捉等方面具有獨特優(yōu)勢。在下節(jié)中,本文將采用多級粒子分辨率技術(shù)對計算域粒子進行處理,結(jié)合δ+-SPH計算結(jié)果中的拉格朗日擬序結(jié)構(gòu)展開對翼型繞流現(xiàn)象研究。

6 Ishii翼型繞流

6.1 Ishii翼型繞流的初始工況和邊界條件設(shè)置

Ishii翼型具有翼厚比低、前緣半徑小等特征,翼型形狀[23]如圖3所示。本文模擬將Ishii翼型固定在矩形流場中,流場左側(cè)為流入邊界,來流速度設(shè)為U0=1 m/s,右側(cè)為流出邊界,上下邊界為自由滑移壁面邊界。翼型的上表面施加無滑移邊界條件。

圖3 攻角α=6°時的Ishii翼型形狀

對Ishii翼型表面的無量綱量升阻力系數(shù)測量計算,獲取翼型的氣動力特性。對阻力系數(shù)CD,升力系數(shù)CL定義為：

(6)

式中，fD和fL分別定義為繞流Ishii翼型的阻力和升力,在SPH程序計算中,這兩者的計算基于流體粒子與結(jié)構(gòu)粒子之間相互作用力的平衡算得。在本文翼型模擬中,弦長L=1 m。

6.2 多級粒子分辨率技術(shù)對計算效率的提升

為提高計算效率,采用6層不同的粒子分辨率,在翼型表面的最高粒子分辨率為L/Δx=800,最低粒子分辨率即距翼型最遠處取L/Δx=25。在本計算中粒子總數(shù)為6.40×105,如果采用L/Δx=800對全流場進行計算,粒子總數(shù)將達到3.07×107,可見采用多級粒子分辨率技術(shù),粒子總數(shù)減少將近43倍,除去不同尺度粒子之間的插值耗時,多級粒子分辨率技術(shù)可將計算效率提高約20倍。

6.3 雷諾數(shù)Re=3 000下繞流模擬結(jié)果和驗證

首先,在雷諾數(shù)Re=3 000、攻角α=0°條件下,模擬了Ishii翼型繞流特性。阻力系數(shù)的計算結(jié)果如圖4所示。FVM結(jié)果與SPH結(jié)果基本吻合,驗證了SPH方法的計算精度。

圖4 雷諾數(shù)Re=3 000, 攻角α=0°Ishii翼型阻力曲線

6.4 雷諾數(shù)Re=30 000下繞流模擬結(jié)果和分析

6.4.1 升阻力系數(shù)

隨后,在雷諾數(shù)Re=30 000、攻角α=0°,3°,6°條件下,模擬了Ishii翼型繞流的特性。α=6°時,翼型阻力及升力系數(shù)的SPH計算結(jié)果與FVM結(jié)果對比如圖5所示。在粒子初始階段到歷經(jīng)劇烈抖動之后,最終CL在0.60～0.85之間動態(tài)波動。同樣地,曲線可發(fā)現(xiàn)阻力CD系數(shù)在0.044上下波動。

圖5 雷諾數(shù)Re=30 000下的升阻力系數(shù)時歷曲線

表1給出了本文計算結(jié)果及文獻中試驗結(jié)果[23]與其他數(shù)值結(jié)果[24](2-D Lam、2-D RANS)的對比。可以看出,多種算法的升阻力系數(shù)值基本吻合,SPH結(jié)果與2-D Lam計算結(jié)果吻合最佳。通過比較圖6中攻角α=0°,3°,6°的升力系數(shù)可以看出,攻角變大、升力系數(shù)增加。此外,攻角增大引起黏性漩渦的非周期性脫落,使得升力系數(shù)隨時間產(chǎn)生振蕩。

表1 升阻力結(jié)果驗證

圖6 不同攻角α下升力系數(shù)時歷曲線

6.4.2 SPH計算黏性繞流問題的速度場與壓力場

圖7展示了雷諾數(shù)Re=30 000、攻角α=6°條件下,Ishii翼型繞流的速度場和粒子分布情況。

圖7 Ishii翼型繞流的速度場和粒子分布

對比t=0.7 s時刻,施加與不施加PST時的粒子分布,可見：不施加PST時,翼型周圍粒子分布十分不規(guī)則,且速度場呈現(xiàn)出不穩(wěn)定;而施加PST,翼型表面邊界層內(nèi)流體粒子分布十分均勻,保證了速度場計算的精度。

在SPH模擬中,實現(xiàn)穩(wěn)定的壓力場計算至關(guān)重要。然而在傳統(tǒng)SPH算法中,常常存在壓力振蕩及張力不穩(wěn)定性現(xiàn)象。在δ+-SPH模型中,通過在連續(xù)方程中添加密度耗散項[14],可得到光滑的壓力場。針對張力不穩(wěn)定現(xiàn)象,本文采用張力不穩(wěn)定性控制技術(shù)進行解決。為驗證其有效性,圖8中給出Ishii翼型繞流模擬中施加TIC與不施加TIC的壓力云圖?？梢?不施加TIC時,翼型周圍流場中出現(xiàn)數(shù)值空洞,導(dǎo)致計算結(jié)果失真;施加TIC時,Ishii翼型周圍壓力分布均勻,且沒有數(shù)值噪聲。

圖8 壓力場云圖(t=4.076 s,Re=30 000)

以上對比說明,δ+-SPH模型結(jié)合TIC技術(shù)可以有效消除壓力噪聲,進而較好實現(xiàn)翼型的壓力數(shù)值變化預(yù)報。此外,觀察壓力云圖可以發(fā)現(xiàn),翼型首部附近區(qū)域存在巨大負壓區(qū),這是引起張力不穩(wěn)定現(xiàn)象的主要原因。因此,采取張力不穩(wěn)定性控制技術(shù)是準確獲取流場特征的重要前提。

6.4.3α=6°時Ishii翼型繞流中漩渦結(jié)構(gòu)分析

圖9 翼型周圍流場中的漩渦結(jié)構(gòu)(t=12.326 s,Re=30 000)

圖9給出了翼型周圍流場中的拉格朗日擬序結(jié)構(gòu)和渦量場?？梢?翼型表面附近高分辨率粒子的布置使得邊界層內(nèi)湍流細節(jié)能能夠清晰地捕捉,反映出多級粒子分辨率技術(shù)在高雷諾數(shù)下黏性繞流問題中的重要作用。相比于渦量場(網(wǎng)格法較常使用該方法捕捉漩渦結(jié)構(gòu)),在SPH框架中捕捉的拉格朗日擬序結(jié)構(gòu)更能清晰重現(xiàn)渦結(jié)構(gòu)的輪廓和內(nèi)部細節(jié)。此外,相較于渦量場云圖,FTLE云圖受閾值的影響不敏感,更加便于進行流場細節(jié)分析。這也是SPH粒子法處理翼型繞流問題的優(yōu)勢之一。

從圖9可見,在翼型形狀及壓力梯度影響下,Ishii翼型漩渦的生成、附著、分離等主要發(fā)生在翼型上緣?；诙嗉壛Ｗ臃直媛始夹g(shù),可以發(fā)現(xiàn),尺度很小的漩渦也能清晰捕捉,再次說明了多級粒子分辨率技術(shù)的有效性。

7 拓展應(yīng)用:拍動翼型繞流問題

為進一步展示SPH在處理運動翼型繞流問題中的優(yōu)勢,本節(jié)采用δ+-SPH模型對拍動水翼問題進行模擬和驗證。在網(wǎng)格法中,對運動翼型繞流問題的模擬容易引起網(wǎng)格變形而降低計算精度。即使采用重疊網(wǎng)格法,也常常因涉及不同網(wǎng)格間的插值而使算法變得復(fù)雜。相比而言,得益于SPH的拉格朗日粒子特性,對于大幅運動邊界問題的模擬較為方便。

本節(jié)選擇頭部半圓尾部尖角的拍動翼型(具體細節(jié)參考文獻[25])進行SPH計算與驗證;采用無量綱參數(shù)包括雷諾數(shù)Re、斯特勞哈爾數(shù)Sr、幅值無量綱系數(shù)AD表征流體黏性、翼型拍動頻率、翼型拍動振幅對翼型繞流模擬的影響。拍動幅值無量綱系數(shù)AD=2A/D,其中,A為翼型拍動最大幅值,D取為翼型頭部半圓直徑。

圖10給出了雷諾數(shù)Re=440,斯特勞哈爾數(shù)Sr=0.035,無量綱幅值系數(shù)AD=1.47條件下,SPH模擬的漩渦結(jié)構(gòu)與試驗結(jié)果[25]的對比。在t1時刻,SPH結(jié)果與試驗結(jié)果[25]吻合較好,體現(xiàn)了δ+-SPH模型在運動翼型繞流模擬上的精度和獨特優(yōu)勢。觀察SPH不同時刻的結(jié)果,可見拉格朗日擬序結(jié)構(gòu)能呈現(xiàn)拍動翼周期拍動時漩渦演化過程,這為未來翼型推進性能與漩渦關(guān)系的研究提供了新思路。

圖10 拍動水翼繞流尾跡的試驗結(jié)果[25]與SPH結(jié)果

8 結(jié) 論

本文采用δ+-SPH方法對不同雷諾數(shù)和攻角下的翼型繞流問題進行了數(shù)值模擬研究,分析了升阻力特性和流場漩渦特征,并將SPH方法從固定翼繞流拓展應(yīng)用到了拍動翼繞流的模擬中。研究結(jié)果表明:

1) SPH方法在模擬真實黏性邊界層和漩渦運動的流體動力學(xué)問題中具有競爭力。不同于網(wǎng)格法,SPH方法易于實現(xiàn)拉格朗日擬序結(jié)構(gòu)的捕捉,動態(tài)呈現(xiàn)漩渦結(jié)構(gòu)演化過程,因此在未來探究漩渦與流體動力特性關(guān)系上具有重要優(yōu)勢。

2) 采用TIC技術(shù)及PST技術(shù)能夠有效提高翼型繞流的速度場和壓力場計算精度,并實現(xiàn)對升阻力系數(shù)較高精度的預(yù)報;采用多級粒子分辨率技術(shù)能夠有效減少粒子總數(shù),縮短計算時間,從而便于在工程中應(yīng)用。

3) SPH方法不需要生成網(wǎng)格,能自動追蹤運動和變形的流固界面,而且該方法的拉格朗日和無網(wǎng)格粒子特性使其在模擬運動翼型繞流問題上具有先天優(yōu)勢和巨大潛力。在未來研究中,利用該優(yōu)勢,SPH方法在翼型繞流的渦振或顫振問題研究中有望發(fā)揮重要作用。