• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于RANS方法的導(dǎo)管推進(jìn)器梢隙流動(dòng)數(shù)值模擬

    2015-08-24 05:47:22阮華張志榮辛公正
    中國(guó)艦船研究 2015年5期
    關(guān)鍵詞:弦長(zhǎng)槳葉湍流

    阮華,張志榮,辛公正

    基于RANS方法的導(dǎo)管推進(jìn)器梢隙流動(dòng)數(shù)值模擬

    阮華,張志榮,辛公正

    中國(guó)船舶科學(xué)研究中心,江蘇無(wú)錫214082

    以導(dǎo)管推進(jìn)器為研究對(duì)象,采用雷諾時(shí)均納維斯托克斯(RANS)方法對(duì)其梢隙流動(dòng)進(jìn)行數(shù)值模擬研究。通過對(duì)網(wǎng)格類型、湍流模型的適用性研究以及對(duì)梢部流場(chǎng)的研究探討,初步建立了基于RANS梢隙流動(dòng)的數(shù)值模擬方法;數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好。對(duì)比分析發(fā)現(xiàn),結(jié)構(gòu)化網(wǎng)格與非結(jié)構(gòu)化網(wǎng)格相比能捕捉到梢隙流動(dòng)中更加細(xì)節(jié)的流場(chǎng)信息,如壁面邊界層流動(dòng)等,更適合于梢隙流動(dòng)的數(shù)值模擬。3種湍流模型SST k-ω,RNG k-ε及RSM的計(jì)算結(jié)果基本一致,都能有效模擬梢隙流動(dòng)。通過間隙區(qū)域流場(chǎng)分析發(fā)現(xiàn),梢隙流動(dòng)的驅(qū)動(dòng)力主要是葉面與葉背之間的壓差,受壁面邊界層流動(dòng)的影響。流體進(jìn)入間隙時(shí)流動(dòng)分離形成間隙分離渦,間隙泄漏流穿過間隙與吸力面?zhèn)攘黧w相互作用,卷起形成梢隙渦,在約37.5%弦長(zhǎng)位置處形成并附著在槳葉壁面發(fā)展,大約在75%弦長(zhǎng)位置與槳葉分離進(jìn)入尾流場(chǎng)中。研究獲得了梢隙渦的起始、發(fā)展、脫落的變化過程。

    導(dǎo)管推進(jìn)器;梢隙流動(dòng);網(wǎng)格類型;湍流模型;機(jī)理分析

    期刊網(wǎng)址:www.ship-research.com

    引用格式:阮華,張志榮,辛公正.基于RANS方法的導(dǎo)管推進(jìn)器梢隙流動(dòng)數(shù)值模擬[J].中國(guó)艦船研究,2015,10(5):83-91.

    RUAN Hua,ZHANG Zhirong,XIN Gongzheng.Numerical simulation of the tip leakage flow in a ducted propulsor based on theRANSmethod[J].Chinese Journalof Ship Research,2015,10(5):83-91.

    0 引言

    導(dǎo)管槳、泵噴等組合推進(jìn)器與敞開式螺旋槳相比,在操縱性能和高速性能方面有明顯的優(yōu)勢(shì)。此類推進(jìn)器通常由轉(zhuǎn)子、定子和導(dǎo)管組成,導(dǎo)管與轉(zhuǎn)子葉梢之間存在間隙。推進(jìn)器工作時(shí),轉(zhuǎn)子槳葉葉面與葉背之間存在壓差,在壓差的作用下槳葉壓力面?zhèn)壬也繀^(qū)域的流體徑向流動(dòng),穿過間隙到吸力面?zhèn)?,并與槳葉間主流相互作用,卷起形成大強(qiáng)度的梢隙泄出渦,常稱之為梢隙渦,這種流動(dòng)稱為梢隙流動(dòng)。

    目前,已有多個(gè)國(guó)家海軍的潛艇采用了泵噴推進(jìn)方式,在泵噴推進(jìn)器中,梢隙間隙流形成的不穩(wěn)定梢隙渦不僅會(huì)使推進(jìn)器效率降低,產(chǎn)生的激振力還容易在渦核內(nèi)低壓區(qū)產(chǎn)生空化,空泡一旦發(fā)生,噪聲會(huì)顯著增加,而低噪聲對(duì)于潛艇的生存至關(guān)重要。船用導(dǎo)管槳也面臨梢隙渦空化帶來(lái)的導(dǎo)管剝蝕和空泡噪聲危害。實(shí)驗(yàn)研究發(fā)現(xiàn),梢隙渦空化的起始先于槳葉表面空化的起始,因此,研究梢隙流動(dòng)對(duì)于減小梢隙渦強(qiáng)度、推遲或抑制渦空化起始、降低噪聲和提高空泡起始航速具有重要意義。

    研究表明,梢隙區(qū)域流動(dòng)過程非常復(fù)雜,為了搞清內(nèi)部流動(dòng)機(jī)理,世界上很多學(xué)者致力于這方面的研究。目前,研究梢隙流動(dòng)的手段主要有實(shí)驗(yàn)測(cè)量與數(shù)值模擬2種。在實(shí)驗(yàn)測(cè)量方面,Judge和Oweis等[1-3]在空泡水筒中對(duì)導(dǎo)管螺旋槳梢隙區(qū)域流場(chǎng)進(jìn)行了激光多普勒測(cè)速(LDV)和粒子成像測(cè)速(PIV)測(cè)量,其對(duì)速度場(chǎng)的分析揭示了梢部區(qū)域漩渦的復(fù)雜相互作用和很強(qiáng)的非定常特性。Kim等[4]通過實(shí)驗(yàn)發(fā)現(xiàn):間隙大小對(duì)軸流式噴水推進(jìn)裝置性能影響很大,間隙小的效率高;此外,對(duì)于等間隙比實(shí)船性能預(yù)報(bào),需要考慮間隙作用的尺度效應(yīng)。Miorini等[5-6]利用高精度的平面PIV觀測(cè)噴水泵轉(zhuǎn)子梢隙流動(dòng)的發(fā)展過程,得到了詳細(xì)的瞬時(shí)和相平均葉梢流動(dòng)的結(jié)構(gòu)和梢隙渦的發(fā)展過程;并用三維PIV研究了噴水泵轉(zhuǎn)子葉梢間隙梢隙渦流動(dòng)結(jié)構(gòu)和湍流動(dòng)力特性。Tan等[7]對(duì)噴水泵轉(zhuǎn)子梢隙渦空泡與壓力面隨邊區(qū)域的空泡進(jìn)行高速攝影,揭示了兩者之間復(fù)雜的相互作用過程,并解釋了空泡是如何影響噴水泵轉(zhuǎn)子的性能。

    在數(shù)值模擬方面,劉登成等[8]采用RANS方法結(jié)合SST k-ω模型比較了不同網(wǎng)格劃分形式對(duì)梢隙流動(dòng)結(jié)果的影響,結(jié)果表明,結(jié)構(gòu)化網(wǎng)格模擬梢隙流動(dòng)有一定的優(yōu)勢(shì)。Yu等[9]使用RANS方法研究梢隙流動(dòng)對(duì)導(dǎo)管槳敞水性能的影響,發(fā)現(xiàn)梢隙流動(dòng)影響了導(dǎo)管和槳葉葉梢的環(huán)量,改變了導(dǎo)管推力和槳推力。Wu等[10]采用RANS和大渦模擬(LES)2種數(shù)值方法對(duì)噴水推進(jìn)器葉梢區(qū)域混雜多種流動(dòng)結(jié)構(gòu)的湍流流場(chǎng)進(jìn)行了數(shù)值模擬,結(jié)果表明,LES與RANS方法相比能獲得更詳盡的流場(chǎng)信息。Cahuzac等[11]對(duì)NACA 65葉片轉(zhuǎn)子梢隙區(qū)域流動(dòng)進(jìn)行了LES數(shù)值模擬,結(jié)果顯示,LES可較全面地重現(xiàn)模型實(shí)驗(yàn)結(jié)果。Hah等[12-13]采用URANS和LES方法研究了噴水泵轉(zhuǎn)子梢隙泄漏流動(dòng)結(jié)構(gòu)和湍流特征,發(fā)現(xiàn)LES數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合更好,對(duì)渦的運(yùn)動(dòng)軌跡和非定常壓力場(chǎng)分布的計(jì)算更為準(zhǔn)確。

    盡管對(duì)梢隙流動(dòng)已有較多研究,但是大部分是從宏觀角度研究梢隙流動(dòng)對(duì)推進(jìn)器性能的影響,對(duì)梢隙流動(dòng)細(xì)節(jié),梢隙泄渦的產(chǎn)生、發(fā)展過程及機(jī)理,梢隙泄渦空化初生等尚未深入研究,無(wú)法預(yù)報(bào)梢隙渦空化起始,數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果相比還有比較大的偏差,梢隙流動(dòng)數(shù)值模擬的方法還不夠成熟。

    由于實(shí)驗(yàn)研究對(duì)于設(shè)備要求高,因此實(shí)驗(yàn)難度大;而隨著CFD技術(shù)的逐漸成熟和計(jì)算機(jī)性能的提升,其為研究梢隙流動(dòng)提供了強(qiáng)大的手段。

    本文將以導(dǎo)管螺旋槳為研究對(duì)象,采用RANS方法對(duì)導(dǎo)管槳的梢隙流動(dòng)進(jìn)行定常數(shù)值模擬,通過對(duì)網(wǎng)格類型、湍流模式的適用性研究以及對(duì)梢部流場(chǎng)的分析探討,初步建立梢隙流動(dòng)的數(shù)值模擬方法。

    1 數(shù)值計(jì)算方法

    由于梢隙流動(dòng)主要以漩渦的形式存在,受流體粘性和邊界層的影響顯著,目前大多采用粘性數(shù)值計(jì)算方法。粘性方法大致可分為3類:直接數(shù)值模擬方法、LES方法和RANS方法。直接數(shù)值模擬最為準(zhǔn)確,但由于對(duì)網(wǎng)格數(shù)量和計(jì)算機(jī)性能的要求太高,目前無(wú)法用于工程計(jì)算。LES方法模擬漩渦流動(dòng)較為準(zhǔn)確,能模擬到很多小尺度的漩渦,但計(jì)算量也比較大。RANS方法已被廣泛應(yīng)用于工程實(shí)踐中,可以有效模擬梢隙流動(dòng)這種大尺度的漩渦流動(dòng),常用來(lái)在推進(jìn)器設(shè)計(jì)過程中分析多方案的優(yōu)劣,效率較高。因此,本文采用RANS方法進(jìn)行數(shù)值模擬研究。

    1.1控制方程

    對(duì)非穩(wěn)態(tài)的N-S方程作時(shí)間平均,將物理量表示為平均量和脈動(dòng)量2部分,得到不可壓縮流體雷諾時(shí)均方程:

    式中:t為時(shí)間;uˉi為流體時(shí)均速度;u'i為流體脈動(dòng)速度;ρ為流體密度;p為壓強(qiáng);fi為質(zhì)量力;v為流體運(yùn)動(dòng)粘性系數(shù);i=1,2,3;j=1,2,3。動(dòng)量方程引入了脈動(dòng)量乘積時(shí)間平均值,稱為雷諾應(yīng)力。這種方法只計(jì)算大尺度平均流動(dòng),而所有湍流脈動(dòng)對(duì)平均流動(dòng)的影響體現(xiàn)到雷諾應(yīng)力τij中。由于在控制方程中引入變量雷諾應(yīng)力造成了方程不封閉,要使方程組封閉,需對(duì)雷諾應(yīng)力采用Boussinesq假設(shè),引入新的湍流模型。本文采用SST k-ω,RNG k-ε和RSM這3種湍流模型計(jì)算。

    1.2研究對(duì)象、網(wǎng)格及數(shù)值方法

    1.2.1研究對(duì)象

    本文研究的導(dǎo)管槳來(lái)自于文獻(xiàn)[1-3],文獻(xiàn)中針對(duì)此導(dǎo)管槳內(nèi)流場(chǎng)進(jìn)行了實(shí)驗(yàn)研究。導(dǎo)管槳的主要參數(shù)和幾何模型分別如表1和圖1所示。槳葉各半徑處剖面弦長(zhǎng)C=0.446DP恒定,導(dǎo)管為等直徑圓管,導(dǎo)管長(zhǎng)度范圍是-0.178~0.091m,槳葉剖面形式為NACA66mod厚度形式和a=0.8拱弧線型。

    表1 導(dǎo)管槳主要參數(shù)Tab.1 Them ain geometry of ducted p ropeller

    圖1 導(dǎo)管槳幾何模型Fig.1 The ducted propellergeometrymodel

    該實(shí)驗(yàn)在空泡水筒中進(jìn)行,由于導(dǎo)管外壁與空泡水筒內(nèi)壁完全貼合,水筒中水流全部經(jīng)過槳盤面,為模擬空泡水筒內(nèi)壁的影響,將導(dǎo)管延長(zhǎng)至槳盤面前方0.385m處。

    1.2.2計(jì)算域大小

    計(jì)算域劃分為2個(gè)不同的區(qū)域(圖2):包含螺旋槳的區(qū)域P-zone和槳外區(qū)域O-zone。P-zone的軸向范圍是-0.05~0.05 m,半徑為0.309 6 m。入口距離槳盤面0.385m,出口距槳盤面10DP,徑向直徑6DP。

    圖2 計(jì)算區(qū)域劃分Fig.2 The computationaldomain

    文獻(xiàn)[1-3]中模型實(shí)驗(yàn)時(shí)導(dǎo)管與水筒內(nèi)壁完全貼合,導(dǎo)管處于流場(chǎng)外邊界,流動(dòng)類似管道流動(dòng),因此加長(zhǎng)導(dǎo)管模擬水筒內(nèi)壁作用。為簡(jiǎn)化對(duì)敞開式水筒的模擬,采用上述計(jì)算域,計(jì)算中根據(jù)實(shí)驗(yàn)進(jìn)速系數(shù)給定導(dǎo)管入口進(jìn)流速度(給定流量),以保證與實(shí)驗(yàn)工況一致,而導(dǎo)管外的流場(chǎng)對(duì)于導(dǎo)管內(nèi)的流場(chǎng)不起作用,這與實(shí)驗(yàn)條件吻合。

    1.2.3網(wǎng)格劃分方法

    O-zone區(qū)域采用全結(jié)構(gòu)化網(wǎng)格劃分方法。P-zone采用2種網(wǎng)格劃分方式:一種為非結(jié)構(gòu)化網(wǎng)格劃分方法;另一種為混合網(wǎng)格,即0.8 r以上區(qū)域網(wǎng)格為結(jié)構(gòu)化而內(nèi)半徑區(qū)域網(wǎng)格為非結(jié)構(gòu)化。計(jì)算區(qū)域網(wǎng)格如圖3所示。

    圖3 計(jì)算域網(wǎng)格劃分Fig.3 The computational grid

    1.2.4數(shù)值方法及邊界條件

    采用Fluent求解器、多參考坐標(biāo)系(MRF),數(shù)值計(jì)算采用RANS定常模擬結(jié)合SST k-ω湍流模型,壓力速度耦合計(jì)算采用SIMPLE算法。邊界條件設(shè)置:入口為速度進(jìn)口邊界,給定均勻來(lái)流速度條件,出口為壓力出口邊界,導(dǎo)管、槳轂、螺旋槳為壁面無(wú)滑移邊界條件,不同區(qū)域之間采用交界面設(shè)置。

    2 計(jì)算結(jié)果與分析

    2.1非結(jié)構(gòu)化網(wǎng)格依賴性分析

    首先針對(duì)外部區(qū)域O-zone網(wǎng)格對(duì)計(jì)算結(jié)果的影響進(jìn)行了研究,當(dāng)該區(qū)域網(wǎng)格數(shù)大于100萬(wàn)時(shí)網(wǎng)格變化對(duì)計(jì)算結(jié)果的影響已經(jīng)很小,因此,選取該區(qū)域網(wǎng)格數(shù)為150萬(wàn)并固定不變。而對(duì)重點(diǎn)關(guān)注區(qū)域P-zone則進(jìn)行了較為詳細(xì)的網(wǎng)格收斂性研究。采用3套網(wǎng)格進(jìn)行計(jì)算,網(wǎng)格數(shù)分別為155萬(wàn)、340萬(wàn)和560萬(wàn)。來(lái)流速度U=5.92m/s,槳葉轉(zhuǎn)速1 200 r/min,進(jìn)速系數(shù)J=0.971。并將計(jì)算的槳葉推力系數(shù)Kt和轉(zhuǎn)矩系數(shù)Kq與實(shí)驗(yàn)值進(jìn)行對(duì)比分析,如表2所示。文獻(xiàn)[1]的模型實(shí)驗(yàn)中,在此進(jìn)速系數(shù)下,槳葉推力系數(shù)Kt=0.31,轉(zhuǎn)矩系數(shù)Kq=0.056。

    表2的計(jì)算結(jié)果表明:P-zone網(wǎng)格數(shù)在340萬(wàn)以上時(shí)計(jì)算結(jié)果基本不再變化,此時(shí)得到的數(shù)值結(jié)果可視為收斂的解;與模型實(shí)驗(yàn)結(jié)果的比較分析,可知Kt,Kq偏差較小,說(shuō)明建立的模型合理有效,能正確反映實(shí)驗(yàn)工況條件。

    表2 網(wǎng)格收斂性分析Tab.2 The convergence analysisof grids

    圖4為采用Q準(zhǔn)則的渦量等值面圖,其中Q=1.0×106,在葉端面靠近壓力面?zhèn)群臀γ鎮(zhèn)纫约皹~導(dǎo)邊附近區(qū)域渦量較大,說(shuō)明這些區(qū)域存在較強(qiáng)的剪切或渦漩流動(dòng),圖中葉背側(cè)長(zhǎng)管狀渦量等值面即是從葉背發(fā)展脫落的梢隙渦。圖5所示為梢隙區(qū)域流線,顯示了梢隙渦形成的大致過程。流體在壓力面?zhèn)犬a(chǎn)生徑向運(yùn)動(dòng),繞過間隙流動(dòng)到吸力面?zhèn)?,與葉背側(cè)主流相互作用形成漩渦并附著在槳葉表面,隨著梢隙流動(dòng)沿弦長(zhǎng)下游不斷增強(qiáng),梢泄渦不斷發(fā)展,渦結(jié)構(gòu)愈發(fā)明顯,渦強(qiáng)也不斷增大。在弦長(zhǎng)中下游位置處,梢隙渦脫離槳葉表面進(jìn)入尾流中。

    圖4 渦量等值面(Q=1.0×106)Fig.4 The vorticitymagnitude isosurface(Q=1.0×106)

    圖5 梢隙區(qū)域流線Fig.5 The stream line in the tip region

    圖6所示為沿弦長(zhǎng)方向不同位置的壓力分布。從圖中可知:葉面與葉背之間的壓差是流體從壓力面?zhèn)冗\(yùn)動(dòng)到吸力面?zhèn)鹊尿?qū)動(dòng)力;由梢隙渦運(yùn)動(dòng)的軌跡可知,渦核的壓力在0.75弦長(zhǎng)位置附近處達(dá)到最低,往下游由于渦結(jié)構(gòu)的變化渦心壓力有所增加。圖7為葉梢剖面0.75弦長(zhǎng)位置間隙區(qū)域的速度矢量圖。該圖揭示了間隙內(nèi)流動(dòng)的特征,壓力面?zhèn)攘黧w在葉面壓差作用下首先往徑向流動(dòng)然后進(jìn)入間隙內(nèi),同時(shí)出現(xiàn)分離流動(dòng),形成漩渦。間隙內(nèi)大部分流體穿過間隙到吸力面?zhèn)龋c吸力面?zhèn)攘黧w相混合,向內(nèi)卷起形成漩渦,即梢泄渦。

    圖6 弦長(zhǎng)方向切片壓力分布Fig.6 The pressure of slicesalong the chord

    圖7 葉梢0.75弦長(zhǎng)位置間隙速度矢量Fig.7 The vector in the tip clearance at0.75 chord p lane

    通過上述流場(chǎng)分析可知,非結(jié)構(gòu)化網(wǎng)格可模擬出典型的梢隙流動(dòng)現(xiàn)象,但是對(duì)流場(chǎng)細(xì)節(jié)信息捕捉不足,比如未計(jì)算出間隙區(qū)葉梢剖面和導(dǎo)管內(nèi)壁面的邊界層,而邊界層的厚度與分離對(duì)梢隙流動(dòng)有很大的影響,同時(shí)模擬的間隙內(nèi)分離渦和梢隙渦的結(jié)構(gòu)也比較模糊。其原因是在如此小的間隙區(qū)域內(nèi)很難劃出精細(xì)化的網(wǎng)格,要想達(dá)到能模擬邊界層流動(dòng)的網(wǎng)格精度,網(wǎng)格數(shù)量將會(huì)很大。而在遠(yuǎn)離壁面區(qū)域,網(wǎng)格分布稀疏,渦核內(nèi)網(wǎng)格分布數(shù)量少,且渦核區(qū)域速度梯度大,流場(chǎng)變化劇烈,對(duì)于梢隙渦的模擬誤差將會(huì)增大。加上非結(jié)構(gòu)網(wǎng)格本身數(shù)值耗散性較大,因此,采用非結(jié)構(gòu)化網(wǎng)格模擬梢隙流動(dòng)是不合適的。

    2.2采用混合網(wǎng)格計(jì)算

    梢隙流動(dòng)重點(diǎn)關(guān)注的是間隙區(qū)域,為提高計(jì)算效率,對(duì)P-zone區(qū)域采用混合網(wǎng)格劃分方法,在槳葉0.8r以上為結(jié)構(gòu)化網(wǎng)格,0.8r以下為非結(jié)構(gòu)化網(wǎng)格,內(nèi)半徑網(wǎng)格按照上述非結(jié)構(gòu)網(wǎng)格生成方法生成,網(wǎng)格數(shù)為360萬(wàn)并固定不變。外半徑0.8r以上區(qū)域表示為P-zone-o,劃分3種不同的網(wǎng)格,網(wǎng)格的主要變化在葉梢剖面梢隙區(qū)域,如表3所示,其余區(qū)域平滑過渡。

    表3 3種梢隙區(qū)域網(wǎng)格Tab.3 Three grid types in tip region

    圖8所示為3套網(wǎng)格計(jì)算的在葉梢剖面0.75弦長(zhǎng)位置處切面吸力面?zhèn)鹊乃俣仁噶繄D。從圖中可以看出:C-Grid間隙區(qū)網(wǎng)格數(shù)量少、尺寸大,未能模擬到明顯的邊界層;M-Grid只模擬到了葉梢剖面的邊界層;F-Grid網(wǎng)格在壁面附近網(wǎng)格更小更密,因此很好地模擬到了間隙區(qū)葉梢剖面和導(dǎo)管內(nèi)壁面的邊界層,同時(shí)與非結(jié)構(gòu)化網(wǎng)格相比在梢隙渦渦核區(qū)分布了更多的計(jì)算網(wǎng)格。速度場(chǎng)分析說(shuō)明,只有在梢隙區(qū)布置相當(dāng)數(shù)量且在壁面附近較密的網(wǎng)格才能模擬到導(dǎo)管內(nèi)壁與葉梢剖面的邊界層流動(dòng),導(dǎo)管內(nèi)壁邊界層比葉梢剖面的邊界層更薄,因此需在導(dǎo)管內(nèi)壁面布置更小尺度的網(wǎng)格,才能模擬到導(dǎo)管邊界層。

    圖8 3套網(wǎng)格計(jì)算的葉梢0.75弦長(zhǎng)位置截面間隙速度矢量Fig.8 The vector in the tip clearance at0.75 chord plane of three grid types

    圖9和圖10分別為結(jié)構(gòu)化網(wǎng)格與非結(jié)構(gòu)化網(wǎng)格計(jì)算的渦量等值面圖。從圖中可以看出,在Q=5.0×106時(shí),非結(jié)構(gòu)化網(wǎng)格與結(jié)構(gòu)化網(wǎng)格相比,梢隙渦渦量范圍小,沒有計(jì)算出槳葉表面和導(dǎo)管內(nèi)壁面的邊界層剪切流動(dòng)渦量,而邊界層流動(dòng)通常與渦的強(qiáng)度和大小密切相關(guān)。因?yàn)榻Y(jié)構(gòu)化網(wǎng)格在間隙區(qū)域內(nèi)的網(wǎng)格分布比非結(jié)構(gòu)化網(wǎng)格精細(xì),其在葉端與導(dǎo)管內(nèi)壁面的邊界層內(nèi)布置了一定數(shù)量的網(wǎng)格,梢隙渦渦核區(qū)也分布了相當(dāng)數(shù)量的網(wǎng)格,對(duì)渦量以及渦結(jié)構(gòu)的計(jì)算更為準(zhǔn)確。而非結(jié)構(gòu)化網(wǎng)格稀疏且網(wǎng)格數(shù)值耗散較大,未能捕捉到明顯的間隙區(qū)域內(nèi)的渦量,渦量計(jì)算偏小。

    圖9 結(jié)構(gòu)化網(wǎng)格渦量等值面(Q=5.0×106)Fig.9 The vorticitymagnitude isosurface(Q=5.0×106)of structured grid

    圖10 非結(jié)構(gòu)化網(wǎng)格渦量等值面(Q=5.0×106)Fig.10 The vorticitymagnitude isosurface(Q=5.0×106)of unstructured grid

    圖11 壓力和流線沿弦向變化Fig.11 The pressure and stream line ofslices

    圖11所示為弦長(zhǎng)方向不同位置截面的壓力分布。在葉片的吸力面?zhèn)?,沿著梢隙渦脫落的方向在渦核區(qū)出現(xiàn)了明顯的低壓區(qū),同時(shí)在壓力面?zhèn)攘黧w進(jìn)入間隙內(nèi),葉梢剖面壁面也出現(xiàn)了低壓區(qū)且在弦長(zhǎng)中下游位置處達(dá)到最低。切面內(nèi)間隙內(nèi)的流線顯示,梢隙渦大約在37.5%弦長(zhǎng)位置處起始發(fā)展,沿著弦長(zhǎng)往下游不斷發(fā)展,大約在75%弦長(zhǎng)位置與槳葉分離進(jìn)入尾流場(chǎng)。通過對(duì)葉梢剖面0.75弦長(zhǎng)位置截面間隙速度矢量圖的分析發(fā)現(xiàn),流體在從壓力面產(chǎn)生徑向流動(dòng)進(jìn)入間隙區(qū)時(shí),由于幾何曲率突變出現(xiàn)了流動(dòng)分離,其與間隙流相互作用形成了漩渦,如圖12所示。

    圖12 葉梢0.75弦長(zhǎng)位置間隙速度矢量Fig.12 Thevector in the tip clearance at0.75 along the chord plane

    圖13和圖14所示為數(shù)值計(jì)算的梢部漩渦Q準(zhǔn)則(Q=5.0×106)識(shí)別與實(shí)驗(yàn)空泡對(duì)梢隙渦的顯示。對(duì)比實(shí)驗(yàn)結(jié)果發(fā)現(xiàn),間隙區(qū)的低壓區(qū)以及梢隙渦的渦結(jié)構(gòu)和運(yùn)動(dòng)軌跡捕捉吻合較好,大約在75%弦長(zhǎng)位置處與槳葉壁面分離進(jìn)入尾流場(chǎng)中。

    圖13 梢部區(qū)域漩渦(計(jì)算)Fig.13 The low pressure zone of vortex in tip region(Cal.)

    圖14 空泡對(duì)漩渦的顯示(實(shí)驗(yàn))Fig.14 The visualization ofvortex in the tip region through cavitation(Exp.)

    由以上分析可知,結(jié)構(gòu)化網(wǎng)格對(duì)梢隙流動(dòng)的模擬比非結(jié)構(gòu)化網(wǎng)格更準(zhǔn)確,可捕捉到更多的流場(chǎng)細(xì)節(jié)。尤其是可獲得導(dǎo)管內(nèi)壁與葉梢剖面邊界層的流動(dòng)與分離,計(jì)算得到壁面邊界層的渦量,這對(duì)于研究梢隙渦的發(fā)展變化具有重要作用;同時(shí),采用結(jié)構(gòu)化網(wǎng)格劃分在梢隙渦渦核內(nèi)分布更多的網(wǎng)格,與非結(jié)構(gòu)化網(wǎng)格相比對(duì)渦結(jié)構(gòu)的計(jì)算更為準(zhǔn)確,結(jié)構(gòu)化網(wǎng)格計(jì)算得到的間隙渦和梢隙渦渦結(jié)構(gòu)及運(yùn)動(dòng)軌跡與實(shí)驗(yàn)結(jié)果吻合較好。

    2.3湍流模式的影響

    上述計(jì)算僅采用了SST k-ω湍流模型進(jìn)行計(jì)算,為研究不同湍流模型對(duì)模擬梢隙流動(dòng)的影響,分別采用RNG k-ε湍流模型結(jié)合標(biāo)準(zhǔn)壁面函數(shù)以及RSM湍流模型進(jìn)行了計(jì)算。

    圖15所示為采用3種不同湍流模型計(jì)算的梢隙區(qū)域渦量流線以及沿弦長(zhǎng)方向截面的壓力分布圖。從圖中可以看出,3種湍流模型計(jì)算的結(jié)果基本相似。圖16所示為3種湍流模型計(jì)算的葉梢0.75弦長(zhǎng)位置截面間隙區(qū)的流線。從中可以看出,圖16(a)和圖16(c)中間隙區(qū)內(nèi)的流動(dòng)差別不大,圖16(b)中間隙內(nèi)部有更加明顯的間隙渦,并有較完整的渦結(jié)構(gòu),即SST k-ω與RSM湍流模型的流場(chǎng)結(jié)果基本一致,而RNG k-ε湍流模型在靠近隨邊的梢隙區(qū)計(jì)算得到了更強(qiáng)的渦量。

    圖15 3種湍流模型計(jì)算的梢部流場(chǎng)Fig.15 The flow field of three turbulencemodels

    圖16 3種湍流模型計(jì)算的葉梢0.75弦長(zhǎng)位置間隙區(qū)的速度矢量Fig.16 The stream line in the tip clearance region at0.75 chord p lane of three turbulencemodels

    由以上分析可知,這3種湍流模型對(duì)計(jì)算梢隙流動(dòng)的差別較小,說(shuō)明3種湍流模型都能對(duì)梢隙流動(dòng)進(jìn)行有效模擬。

    3 結(jié)論

    本文基于RANS方法開展了導(dǎo)管槳梢隙流動(dòng)的數(shù)值模擬研究,對(duì)不同網(wǎng)格數(shù)量、網(wǎng)格類型及湍流模型進(jìn)行系列計(jì)算分析,得出以下主要結(jié)論:

    1)計(jì)算結(jié)果表明,采用結(jié)構(gòu)化網(wǎng)格更適合于對(duì)流場(chǎng)變化劇烈且以漩渦運(yùn)動(dòng)為主要形式的梢隙流動(dòng),同時(shí),在間隙區(qū)域需要布置更精細(xì)的網(wǎng)格才能模擬得到梢隙流動(dòng)更加全面的精細(xì)流場(chǎng)信息。

    2)通過對(duì)梢隙流動(dòng)的流場(chǎng)分析發(fā)現(xiàn),梢隙流動(dòng)的驅(qū)動(dòng)力主要是葉面與葉背之間的壓差,其受壁面邊界層流動(dòng)的影響。流體進(jìn)入間隙時(shí)流動(dòng)分離形成間隙分離渦,間隙泄漏流穿過間隙與吸力面?zhèn)攘黧w相互作用,卷起形成梢隙渦,在約37.5%弦長(zhǎng)位置處形成并附著于槳葉壁面發(fā)展,大約在75%弦長(zhǎng)位置處與槳葉分離進(jìn)入尾流場(chǎng)中。

    3)通過計(jì)算獲得了梢隙渦的起始、發(fā)展、脫落的變化過程。

    4)3種湍流模型的計(jì)算結(jié)果顯示,梢隙流動(dòng)的數(shù)值模擬結(jié)果差別不大,SST k-ω與RSM湍流模型數(shù)值模擬的結(jié)果基本一致,RNG k-ε湍流模型結(jié)合壁面函數(shù)在間隙區(qū)能模擬到更加明顯的間隙渦,說(shuō)明3種湍流模型都能對(duì)梢隙流動(dòng)進(jìn)行有效模擬。

    通過本文的研究,建立了基于RANS的梢隙流動(dòng)的數(shù)值模擬方法,且經(jīng)過了實(shí)驗(yàn)的初步驗(yàn)證。

    [1] JUDGE CQ,OWEISG F,CECCIO SL,etal.Tip leakage vortex inceptiononaductedrotor[C/OL]// CAV2001:Session A 6.001[2001-06-06].http:// core.ac.uk/download/pdf/9412518.pdf.

    [2]OWEISG,CECCIO S L,CHESNAKAS C,et al.Tip leakage vortex(TLV)variability from a ducted propeller under steady operation and its imp lications on cavitation incep tion[C]//Proceedings of Fifth International Symposium on Cavitation(CAV2003).Osaka,Japan,2003.

    [3] OWEISG F,CECCIOSL.Instantaneousand time-averaged flow fields ofmultip le vortices in the tip region of a ducted propulsor[J].Experiments in Fluids,2005,38(5):615-636.

    [4]KIM M C,CHUN H H.Experimental investigation into the performance of the axial-flow-type waterjet according to the variation of impeller tip clearance[J]. Ocean Engineering,2007,34(2):275-283.

    [5] MIORINIR L,WU H X,KATZ J.The internal structure of the tip leakage vortexwithin the rotor ofan axial waterjet pump[J].Journal of Turbomachinery,2012,134(3):031018-1-031018-12.

    [6] MIORINIR L,WU H X,TAN D,etal.Three-dimensional structure within the tip leakage vortex ofan axial waterjetpump[C]//Proceedingofthe ASME-JSME-KSME 2011 Joint Fluid EngineeringConference(AJK-Fluids 2011).Hamamastu,Shizuoka,Japan,2011.

    [7]TAN D,LI Y C,MIORINI R,et al.Role of large scale cavitating vortical structures in the rotor passage of an axial waterjet pump in performance breakdown[C]//Proceedings of 30th Sym posium on Naval Hyd rodynam ics.Hobart,Tasmania,Australia,2014.

    [8] 劉登成,張志榮,黃國(guó)富,等.帶前置定子導(dǎo)管槳梢隙流動(dòng)數(shù)值分析[C]//第二十一屆全國(guó)水動(dòng)力學(xué)研討會(huì)暨第八屆全國(guó)水動(dòng)力學(xué)學(xué)術(shù)會(huì)議暨兩岸船舶與海洋工程水動(dòng)力學(xué)研討會(huì)文集.無(wú)錫:水動(dòng)力學(xué)研究與進(jìn)展編委會(huì),2008:695-703.

    [9] YU L,GREVE M,DRUCKENBROD M,et al.Numerical analysis of ducted propeller performance under open water test condition[J].Journal of Marine Science and Technology,2013,18:381-394.

    [10] WU H X,MIORINIR L,KATZ J.Analysis of turbulence in the tip region of a waterjet pump rotor[C]// Proceedings of the ASME 2010 3rd Joint US-European Fluids Engineering Summer Meeting and 8th International Conference on Nanochannels,Microchannels,and Minichannels.American Society ofMechanical Engineers,2010,1:699-711.

    [11] CAHUZAC A,BOUDET J,JACOB M C,et al. Large-eddy simulation of a rotor tip-clearance flow[C]//17th AIAA/CEAS Aeroacoustics Con ference(32nd AIAA Aeroacoustics Conference).Portland,Oregon,2011.AIAA 2011-2947.

    [12]HAH C.Investigation of turbulent tip leakage vortex in an axial water jet pump with large eddy simulation[C]//Proceedings of the ASME 2012 Fluids Engineering Summer Meeting(FEDSM 2012).American Society ofMechanicalEngineers,2012.

    [13]HAH C,KATZ J.Investigation of tip leakage flow in an axial waterjet pump with large eddy simulation[C]//30th Symposium on Naval Hydrodynam ics.Hobart,Tasmania,Australia,2014.

    [責(zé)任編輯:易基圣]

    Num ericalsim u lation of the tip leakage flow in a ducted p ropu lsor based on the RANSmethod

    RUAN Hua,ZHANGZhirong,XINGongzheng
    China Ship Scientific Research Center,Wuxi214082,China

    In this paper,a ducted propulsor is taken as the objectof the research,where the RANSsimulation method isused tostudy the tip leakage flow of thepropulsor.Through the applicability study of the numerical grid type and the turbulencemodel,and by analyzing the flow field at the tip region,preliminary numerical methods based on RANS for tip leakage flow are finally proposed,and the corresponding results agreewell with the experimental data.Through the comparison of unstructured and structuredmesh's results,it is observed that structured mesh can capture more flowing detail(for example,the boundary flow).Therefore,the structured mesh ismore appropriate for the numerical simulation of tip leakage flow.Meanwhile,three turbulencemodels,SST k-ω,RNG k-ε,and RSM are used,and their numerical results are essentially similar,which suggests thatall three can be app lied in the tip leakage flow simulation.Through the analysis of tip clearance region's flow field,it is found that tip leakage flow's driving power is primary the p ressure difference between the front and back side of the rotor,and the flow is affected by the wall boundary flow. When the flow enters the tip gap,flow separation is occurred,and a gap separation vortex is formed.The tip clearance flow goes across the gap and interactswith the passage flow,finally developing into Tip Leakage Vortex(TLV).TLV's incep tion position is about 37.5%cord length and attaches to the rotor wall,and then sheds from the rotorwall atapproximately 75%cord length into thewake flied.The numerical results clearly demonstrate the changing processof the tip leakage vortex inception,development,and shedding.

    ducted propulsor;tip leakage flow;grid type;turbulentmodel;mechanism analysis

    U661.1

    ADO I:10.3969/j.issn.1673-3185.2015.05.014

    2014-12-22網(wǎng)絡(luò)出版時(shí)間:2015-10-8 11∶15

    阮華(通信作者),男,1989年生,碩士。研究方向:船舶計(jì)算流體力學(xué)。E-mail:ruanhuahust@163.com

    張志榮,男,1966年生,博士,研究員。研究方向:計(jì)算流體力學(xué)。E-mail:zhangzr@cssrc.com.cn

    網(wǎng)絡(luò)出版地址:http∶//www.cnki.net/kcms/detail/42.1755.TJ.20151008.1115.034.htm l

    猜你喜歡
    弦長(zhǎng)槳葉湍流
    探究奇偶旋翼對(duì)雷達(dá)回波的影響
    淺談圓錐曲線三類弦長(zhǎng)問題
    弦長(zhǎng)積分的極限性質(zhì)與不等式
    立式捏合機(jī)槳葉結(jié)構(gòu)與槳葉變形量的CFD仿真*
    重氣瞬時(shí)泄漏擴(kuò)散的湍流模型驗(yàn)證
    弦長(zhǎng)積分的極限性質(zhì)與不等式
    直升機(jī)槳葉/吸振器系統(tǒng)的組合共振研究
    立式捏合機(jī)槳葉型面設(shè)計(jì)與優(yōu)化①
    “青春期”湍流中的智慧引渡(三)
    “青春期”湍流中的智慧引渡(二)
    亚洲五月色婷婷综合| 精品久久久久久久毛片微露脸| 男人舔女人的私密视频| 色老头精品视频在线观看| 亚洲av中文字字幕乱码综合 | 亚洲av片天天在线观看| 久久久久久久久免费视频了| 国产区一区二久久| 国产精品九九99| 日本五十路高清| 精品一区二区三区av网在线观看| 亚洲第一欧美日韩一区二区三区| 一本一本综合久久| 国产精品av久久久久免费| 日韩欧美一区二区三区在线观看| 国产精品久久久久久亚洲av鲁大| av免费在线观看网站| 伦理电影免费视频| 18禁黄网站禁片午夜丰满| 亚洲黑人精品在线| 欧美激情极品国产一区二区三区| 国产主播在线观看一区二区| 久久精品91蜜桃| 老司机福利观看| 别揉我奶头~嗯~啊~动态视频| 欧美成人性av电影在线观看| 亚洲精品美女久久久久99蜜臀| 久久精品人妻少妇| 身体一侧抽搐| 91老司机精品| 国产极品粉嫩免费观看在线| 一级片免费观看大全| 国产黄a三级三级三级人| 一进一出好大好爽视频| 人人妻人人澡欧美一区二区| 激情在线观看视频在线高清| 窝窝影院91人妻| 在线播放国产精品三级| 午夜精品久久久久久毛片777| 免费人成视频x8x8入口观看| 亚洲精品美女久久av网站| 中文字幕人成人乱码亚洲影| 老汉色∧v一级毛片| 亚洲最大成人中文| 看黄色毛片网站| 国产aⅴ精品一区二区三区波| 亚洲成国产人片在线观看| 亚洲欧美激情综合另类| а√天堂www在线а√下载| 在线免费观看的www视频| 黄网站色视频无遮挡免费观看| 欧美精品亚洲一区二区| 亚洲最大成人中文| 国产亚洲精品综合一区在线观看 | 男人操女人黄网站| 十八禁人妻一区二区| a级毛片在线看网站| 欧美日韩亚洲综合一区二区三区_| 久久香蕉国产精品| 欧美激情 高清一区二区三区| 99re在线观看精品视频| 97人妻精品一区二区三区麻豆 | 亚洲av五月六月丁香网| 国产精品乱码一区二三区的特点| 国产亚洲精品一区二区www| 久久热在线av| 九色国产91popny在线| 国产97色在线日韩免费| 精品久久蜜臀av无| 久久这里只有精品19| 日韩一卡2卡3卡4卡2021年| 99在线视频只有这里精品首页| 免费观看人在逋| 久久久久久大精品| 成人亚洲精品一区在线观看| 香蕉av资源在线| 日韩欧美一区二区三区在线观看| 久久伊人香网站| 十分钟在线观看高清视频www| 麻豆国产av国片精品| 高清毛片免费观看视频网站| 色播亚洲综合网| 国产精品九九99| 少妇粗大呻吟视频| 少妇 在线观看| 少妇熟女aⅴ在线视频| 国产精品乱码一区二三区的特点| 国产亚洲精品av在线| 午夜老司机福利片| 久久久久九九精品影院| 久久久久国内视频| 久久亚洲真实| 亚洲 欧美 日韩 在线 免费| 九色国产91popny在线| 国产精品一区二区免费欧美| 男女那种视频在线观看| 精品国产乱子伦一区二区三区| 亚洲av电影不卡..在线观看| 一卡2卡三卡四卡精品乱码亚洲| 精品午夜福利视频在线观看一区| 国产v大片淫在线免费观看| 成年人黄色毛片网站| 久久久久久久精品吃奶| 在线天堂中文资源库| av欧美777| 男女午夜视频在线观看| 日韩精品青青久久久久久| 精品国产亚洲在线| 午夜影院日韩av| 最好的美女福利视频网| 国产精品综合久久久久久久免费| 法律面前人人平等表现在哪些方面| 精品久久蜜臀av无| 色尼玛亚洲综合影院| 两个人视频免费观看高清| 亚洲熟妇中文字幕五十中出| 国产一卡二卡三卡精品| 日本 av在线| 国产不卡一卡二| 成人三级黄色视频| 少妇被粗大的猛进出69影院| 最好的美女福利视频网| 亚洲自拍偷在线| 99国产极品粉嫩在线观看| 久久久久久九九精品二区国产 | 精品久久久久久成人av| 国产精品久久久久久人妻精品电影| 亚洲,欧美精品.| 欧美又色又爽又黄视频| 国产黄片美女视频| 一级黄色大片毛片| 亚洲国产欧美网| 欧美激情久久久久久爽电影| 老司机深夜福利视频在线观看| 久久国产精品影院| 国产精品国产高清国产av| 欧美激情久久久久久爽电影| 在线观看舔阴道视频| 亚洲 欧美一区二区三区| www国产在线视频色| 人人澡人人妻人| 我的亚洲天堂| 国产精品野战在线观看| 久久久久久久久免费视频了| 高清在线国产一区| 精品久久蜜臀av无| 亚洲免费av在线视频| 在线观看日韩欧美| 国产区一区二久久| 欧美中文日本在线观看视频| 色综合欧美亚洲国产小说| 50天的宝宝边吃奶边哭怎么回事| 中文资源天堂在线| 又黄又粗又硬又大视频| 激情在线观看视频在线高清| 大型黄色视频在线免费观看| 一a级毛片在线观看| 婷婷精品国产亚洲av在线| 美女免费视频网站| 国产私拍福利视频在线观看| 欧美日韩中文字幕国产精品一区二区三区| 免费搜索国产男女视频| 老司机午夜十八禁免费视频| 久久 成人 亚洲| 搡老熟女国产l中国老女人| 国产黄片美女视频| 午夜福利在线在线| 99热只有精品国产| 老鸭窝网址在线观看| 中文字幕人妻熟女乱码| 法律面前人人平等表现在哪些方面| 国产主播在线观看一区二区| 草草在线视频免费看| 精品电影一区二区在线| 亚洲精品在线美女| 丝袜在线中文字幕| 久久精品人妻少妇| 一级a爱片免费观看的视频| 精品久久久久久,| 美女午夜性视频免费| 国产成人系列免费观看| 天堂动漫精品| 国产av一区二区精品久久| 午夜a级毛片| 琪琪午夜伦伦电影理论片6080| 国产精品久久久人人做人人爽| 一进一出好大好爽视频| 亚洲国产欧洲综合997久久, | 中文字幕人妻熟女乱码| 精品久久久久久久末码| 国产欧美日韩一区二区精品| 91国产中文字幕| 两性夫妻黄色片| 精品日产1卡2卡| 亚洲男人的天堂狠狠| 亚洲人成网站高清观看| 国产精品电影一区二区三区| 久久精品国产综合久久久| 好看av亚洲va欧美ⅴa在| 亚洲一区二区三区不卡视频| 美国免费a级毛片| 女性被躁到高潮视频| 曰老女人黄片| 天天一区二区日本电影三级| 中文亚洲av片在线观看爽| 国产一区二区三区视频了| 日韩欧美国产一区二区入口| 啦啦啦免费观看视频1| av欧美777| 国产三级黄色录像| 色精品久久人妻99蜜桃| 91国产中文字幕| 别揉我奶头~嗯~啊~动态视频| 国产成人一区二区三区免费视频网站| xxxwww97欧美| 久久伊人香网站| 午夜免费鲁丝| 免费在线观看日本一区| 国产成+人综合+亚洲专区| 精品福利观看| 黄色成人免费大全| 免费在线观看影片大全网站| 此物有八面人人有两片| 亚洲第一av免费看| 在线永久观看黄色视频| 久久久久久国产a免费观看| 真人一进一出gif抽搐免费| 国内久久婷婷六月综合欲色啪| 熟妇人妻久久中文字幕3abv| 成人av一区二区三区在线看| 男女下面进入的视频免费午夜 | 精品久久久久久久毛片微露脸| 母亲3免费完整高清在线观看| 欧美色欧美亚洲另类二区| 黄片播放在线免费| 窝窝影院91人妻| 亚洲av成人一区二区三| 侵犯人妻中文字幕一二三四区| 色在线成人网| 搡老岳熟女国产| 国产成年人精品一区二区| 欧美激情极品国产一区二区三区| 国产一区二区激情短视频| 人成视频在线观看免费观看| 两性夫妻黄色片| 久99久视频精品免费| 可以免费在线观看a视频的电影网站| 黄色毛片三级朝国网站| 夜夜看夜夜爽夜夜摸| 亚洲真实伦在线观看| 欧美黄色片欧美黄色片| 日本a在线网址| 麻豆成人午夜福利视频| 老司机午夜福利在线观看视频| 日本三级黄在线观看| 真人做人爱边吃奶动态| 亚洲中文字幕日韩| 午夜免费成人在线视频| 亚洲精品色激情综合| 久久久久国内视频| 一卡2卡三卡四卡精品乱码亚洲| 国产精品亚洲美女久久久| 国产一级毛片七仙女欲春2 | 日本精品一区二区三区蜜桃| 国产亚洲欧美在线一区二区| 性色av乱码一区二区三区2| 青草久久国产| 欧美色欧美亚洲另类二区| 一本精品99久久精品77| 淫妇啪啪啪对白视频| 午夜两性在线视频| 日韩欧美 国产精品| 婷婷精品国产亚洲av在线| 丁香欧美五月| 国产又色又爽无遮挡免费看| 久久人妻福利社区极品人妻图片| 精品欧美国产一区二区三| 99久久99久久久精品蜜桃| 麻豆成人av在线观看| 韩国精品一区二区三区| 亚洲片人在线观看| 国产午夜精品久久久久久| 久久亚洲真实| 国产亚洲欧美98| 亚洲在线自拍视频| 黑人巨大精品欧美一区二区mp4| 嫩草影视91久久| 性色av乱码一区二区三区2| 性欧美人与动物交配| 免费在线观看视频国产中文字幕亚洲| 久久狼人影院| 在线观看午夜福利视频| 欧美大码av| 精品少妇一区二区三区视频日本电影| 午夜福利一区二区在线看| 国产精品久久久久久人妻精品电影| 欧美日韩亚洲国产一区二区在线观看| 美女国产高潮福利片在线看| 久热这里只有精品99| 午夜福利在线观看吧| 久久午夜综合久久蜜桃| 精品国产亚洲在线| 久久久水蜜桃国产精品网| 久99久视频精品免费| 国产成人精品久久二区二区免费| 婷婷精品国产亚洲av在线| 精品国产乱码久久久久久男人| 欧美中文日本在线观看视频| 在线十欧美十亚洲十日本专区| 丁香欧美五月| 母亲3免费完整高清在线观看| www.熟女人妻精品国产| 国产视频一区二区在线看| 宅男免费午夜| 两性夫妻黄色片| 久久 成人 亚洲| 免费看a级黄色片| 亚洲av第一区精品v没综合| 国产亚洲欧美精品永久| 国产精品一区二区精品视频观看| 色综合亚洲欧美另类图片| 99riav亚洲国产免费| 久久青草综合色| 亚洲av成人不卡在线观看播放网| 淫秽高清视频在线观看| 国产97色在线日韩免费| 白带黄色成豆腐渣| 91麻豆精品激情在线观看国产| 一级片免费观看大全| 黄色丝袜av网址大全| 国产一级毛片七仙女欲春2 | 一卡2卡三卡四卡精品乱码亚洲| 婷婷精品国产亚洲av| 亚洲精品色激情综合| 国产真实乱freesex| 成人国产综合亚洲| АⅤ资源中文在线天堂| 这个男人来自地球电影免费观看| 久久午夜综合久久蜜桃| 老汉色∧v一级毛片| 听说在线观看完整版免费高清| 久久久国产成人精品二区| 精品福利观看| 99精品欧美一区二区三区四区| 99久久国产精品久久久| 久久精品亚洲精品国产色婷小说| 日本 av在线| 在线观看免费日韩欧美大片| 一级毛片女人18水好多| 亚洲av成人不卡在线观看播放网| 黄频高清免费视频| 午夜福利成人在线免费观看| 一二三四在线观看免费中文在| 十八禁网站免费在线| 国内精品久久久久久久电影| 在线十欧美十亚洲十日本专区| 国产欧美日韩一区二区三| 中文字幕人妻熟女乱码| 男男h啪啪无遮挡| 亚洲第一av免费看| www.www免费av| 啦啦啦韩国在线观看视频| 男女午夜视频在线观看| 亚洲五月天丁香| 国产精品电影一区二区三区| 亚洲自拍偷在线| 视频在线观看一区二区三区| 国内精品久久久久精免费| 婷婷精品国产亚洲av在线| 久久亚洲真实| 精品国内亚洲2022精品成人| 亚洲美女黄片视频| 啪啪无遮挡十八禁网站| 午夜免费鲁丝| 亚洲精品一区av在线观看| 色精品久久人妻99蜜桃| 日日摸夜夜添夜夜添小说| 18禁黄网站禁片免费观看直播| 俺也久久电影网| 人人妻人人看人人澡| avwww免费| 后天国语完整版免费观看| 亚洲国产欧美网| 午夜福利一区二区在线看| 在线观看免费午夜福利视频| 两性午夜刺激爽爽歪歪视频在线观看 | 亚洲在线自拍视频| 欧美+亚洲+日韩+国产| 满18在线观看网站| 看免费av毛片| 巨乳人妻的诱惑在线观看| 欧美在线一区亚洲| 啦啦啦观看免费观看视频高清| 精品久久久久久,| 禁无遮挡网站| 精品国产国语对白av| 男女做爰动态图高潮gif福利片| 国产真实乱freesex| 欧美激情久久久久久爽电影| 哪里可以看免费的av片| 亚洲第一欧美日韩一区二区三区| 999久久久精品免费观看国产| 国内毛片毛片毛片毛片毛片| 久久人妻av系列| 国产精品久久久久久亚洲av鲁大| 欧美黑人巨大hd| 一级毛片精品| 欧美黑人精品巨大| 国产极品粉嫩免费观看在线| 亚洲av五月六月丁香网| 美女高潮到喷水免费观看| 亚洲精品国产一区二区精华液| 国产高清视频在线播放一区| 日韩av在线大香蕉| 久久国产亚洲av麻豆专区| 一区福利在线观看| 日韩有码中文字幕| 国产精品一区二区精品视频观看| 18禁国产床啪视频网站| 亚洲成国产人片在线观看| 久久精品91蜜桃| 色播亚洲综合网| 欧美乱码精品一区二区三区| 欧美午夜高清在线| 久久久久久亚洲精品国产蜜桃av| 日日爽夜夜爽网站| 999久久久国产精品视频| 亚洲精品国产精品久久久不卡| 国产极品粉嫩免费观看在线| 99国产精品99久久久久| 此物有八面人人有两片| 露出奶头的视频| 非洲黑人性xxxx精品又粗又长| 国产成人精品久久二区二区免费| 国产精品影院久久| 中文字幕另类日韩欧美亚洲嫩草| 久久久久国产一级毛片高清牌| 美女国产高潮福利片在线看| 日韩欧美一区视频在线观看| 久久精品人妻少妇| 午夜福利一区二区在线看| tocl精华| 国产免费男女视频| 午夜免费观看网址| 精品无人区乱码1区二区| 欧美精品啪啪一区二区三区| 亚洲av美国av| 久久久久久亚洲精品国产蜜桃av| 少妇 在线观看| 色综合站精品国产| 国产精品久久电影中文字幕| 每晚都被弄得嗷嗷叫到高潮| 黄频高清免费视频| 久久久久久久久免费视频了| 国产三级在线视频| 亚洲精品久久国产高清桃花| 亚洲中文日韩欧美视频| 国产午夜精品久久久久久| www.精华液| 国产精品综合久久久久久久免费| а√天堂www在线а√下载| 国产精品一区二区免费欧美| 午夜影院日韩av| 桃红色精品国产亚洲av| 91九色精品人成在线观看| 久久精品国产清高在天天线| 亚洲欧美激情综合另类| 又紧又爽又黄一区二区| 一个人免费在线观看的高清视频| 久久伊人香网站| 神马国产精品三级电影在线观看 | 亚洲九九香蕉| 国产在线精品亚洲第一网站| 无遮挡黄片免费观看| 一本久久中文字幕| 91九色精品人成在线观看| 欧美激情 高清一区二区三区| 两人在一起打扑克的视频| 亚洲五月色婷婷综合| 欧美性猛交╳xxx乱大交人| 最新美女视频免费是黄的| 成人国产一区最新在线观看| 久久天躁狠狠躁夜夜2o2o| 精品国产国语对白av| 国产91精品成人一区二区三区| 欧美丝袜亚洲另类 | 国产精品久久久久久亚洲av鲁大| 视频在线观看一区二区三区| 韩国av一区二区三区四区| 亚洲精品色激情综合| 色综合欧美亚洲国产小说| 日本免费一区二区三区高清不卡| 男女视频在线观看网站免费 | 一级毛片高清免费大全| 黄色女人牲交| 首页视频小说图片口味搜索| 99久久久亚洲精品蜜臀av| av超薄肉色丝袜交足视频| 欧美乱码精品一区二区三区| 岛国视频午夜一区免费看| 精品日产1卡2卡| 嫩草影视91久久| 久久性视频一级片| 老熟妇乱子伦视频在线观看| 亚洲无线在线观看| 观看免费一级毛片| 一a级毛片在线观看| 非洲黑人性xxxx精品又粗又长| 中文亚洲av片在线观看爽| 两个人视频免费观看高清| 欧美性长视频在线观看| 国产高清有码在线观看视频 | 99热6这里只有精品| 亚洲精品中文字幕在线视频| 麻豆av在线久日| 日韩精品中文字幕看吧| 国产精品自产拍在线观看55亚洲| 啦啦啦韩国在线观看视频| 日韩欧美三级三区| 丰满人妻熟妇乱又伦精品不卡| 人人妻人人澡欧美一区二区| 免费在线观看黄色视频的| 国产一级毛片七仙女欲春2 | 国产精品免费视频内射| av视频在线观看入口| 在线观看舔阴道视频| 美女高潮到喷水免费观看| 欧美一级a爱片免费观看看 | 国产精品一区二区三区四区久久 | 黄网站色视频无遮挡免费观看| 男女视频在线观看网站免费 | 婷婷精品国产亚洲av| 满18在线观看网站| 熟妇人妻久久中文字幕3abv| 黄色视频,在线免费观看| 麻豆一二三区av精品| 久久久水蜜桃国产精品网| 亚洲人成网站高清观看| 亚洲激情在线av| 国产高清videossex| 日韩欧美一区视频在线观看| 色av中文字幕| 在线视频色国产色| 首页视频小说图片口味搜索| 免费看a级黄色片| x7x7x7水蜜桃| 一级a爱视频在线免费观看| 欧美性猛交黑人性爽| 波多野结衣巨乳人妻| 一级毛片高清免费大全| 亚洲男人天堂网一区| 久久国产精品人妻蜜桃| 久久中文看片网| 国产aⅴ精品一区二区三区波| 香蕉av资源在线| 国产aⅴ精品一区二区三区波| 香蕉av资源在线| 最新美女视频免费是黄的| 18禁国产床啪视频网站| 制服人妻中文乱码| 精品乱码久久久久久99久播| 男女那种视频在线观看| 国产成年人精品一区二区| 精品一区二区三区av网在线观看| 深夜精品福利| 日韩大尺度精品在线看网址| 久久精品国产清高在天天线| 亚洲人成网站高清观看| 琪琪午夜伦伦电影理论片6080| 国产亚洲精品一区二区www| av视频在线观看入口| 精品少妇一区二区三区视频日本电影| 一区二区三区精品91| 国产精品亚洲一级av第二区| 中亚洲国语对白在线视频| 婷婷亚洲欧美| 亚洲精品中文字幕一二三四区| 国产精品亚洲美女久久久| av在线天堂中文字幕| 91在线观看av| 天天添夜夜摸| 韩国av一区二区三区四区| 视频区欧美日本亚洲| 50天的宝宝边吃奶边哭怎么回事| 香蕉av资源在线| svipshipincom国产片| 在线观看一区二区三区| 成人亚洲精品av一区二区| 欧美在线黄色| 色av中文字幕| 一区福利在线观看| 日韩欧美在线二视频| 国产高清激情床上av| 精品久久久久久成人av| 亚洲精品久久成人aⅴ小说| www.精华液| 在线永久观看黄色视频| 正在播放国产对白刺激| 宅男免费午夜| 黄片播放在线免费| 亚洲av成人av| 中文字幕av电影在线播放| 女性生殖器流出的白浆| 午夜日韩欧美国产| 国产精品精品国产色婷婷| 亚洲国产日韩欧美精品在线观看 | 久久久久久国产a免费观看| 国产99久久九九免费精品| 国产欧美日韩精品亚洲av| 免费观看精品视频网站| 变态另类丝袜制服| 一级毛片女人18水好多| 久久久水蜜桃国产精品网| 中文字幕av电影在线播放| 99精品欧美一区二区三区四区|