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

    水下超空泡流動(dòng)問題數(shù)值計(jì)算研究綜述

    2023-07-10 02:26:54孫鐵志王世晟謝勃漢
    關(guān)鍵詞:模型研究

    孫鐵志 ,王世晟 ,謝勃漢 ,許 昊 *

    (1.大連理工大學(xué) 船舶工程學(xué)院,遼寧 大連,116024;2.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析優(yōu)化與CAE 軟件全國(guó)重點(diǎn)實(shí)驗(yàn)室,遼寧 大連,116024)

    0 引言

    超空泡技術(shù)是指通過頭部空化器、特殊頭型或氣體發(fā)生器等方法,在水下形成完全包裹航行體的空泡,以減小航行體整體阻力。采用水下超空泡技術(shù)的航行體,能夠?qū)崿F(xiàn)更高的航行速度,相對(duì)一般水下裝備有明顯的速度優(yōu)勢(shì)[1-2]。

    由于超空泡減阻的巨大潛在優(yōu)勢(shì),國(guó)內(nèi)外研究者已對(duì)超空泡發(fā)生、演化流動(dòng)機(jī)理和航行體水動(dòng)力特性等方面開展了廣泛研究。對(duì)于超空泡流動(dòng)的研究最初基于勢(shì)流理論開展,忽略黏性影響后得到了超空泡的形態(tài)及升阻力計(jì)算公式。但是實(shí)際超空泡流動(dòng)中存在明顯氣液交界面、多相流動(dòng)中的相變過程以及高雷諾數(shù)下的湍流特征,流場(chǎng)復(fù)雜,非線性和非定常特征明顯。通過勢(shì)流理論無法準(zhǔn)確預(yù)報(bào)這些特征帶來的影響,因此,試驗(yàn)研究成為了重要的研究方法[3-7]。然而,試驗(yàn)研究難以獲得泡內(nèi)全局壓力場(chǎng)、速度場(chǎng)、渦量場(chǎng)等流場(chǎng)細(xì)節(jié)和多相耦合作用特性。通過數(shù)值計(jì)算手段,基于界面捕捉方法獲得超空泡形態(tài)、相變模型模擬空化過程、湍流模型計(jì)算空泡表面脈動(dòng)及尾空泡破碎脫落,能夠計(jì)算水下超空泡流動(dòng)的復(fù)雜流場(chǎng)特征,獲得更全面細(xì)致的流場(chǎng)信息,為研究超空泡機(jī)理提供幫助。

    水下超空泡流動(dòng)作為一種復(fù)雜的流動(dòng)過程,計(jì)算方法對(duì)這一物理問題的有效預(yù)報(bào)十分重要。文中從理論模型方法和數(shù)值仿真方法出發(fā),介紹了水下超空泡數(shù)值計(jì)算的研究進(jìn)展。

    1 理論模型求解方法

    超空泡理論求解方法是基于理想流體勢(shì)流理論,通過數(shù)學(xué)方法處理復(fù)變函數(shù)的計(jì)算方法,能夠預(yù)報(bào)空泡形態(tài)及升阻力,但是理想流體假設(shè)忽略了黏性影響,對(duì)摩擦阻力預(yù)估不足。這部分方法包括自由流線理論、線性理論和獨(dú)立膨脹原理等。

    1.1 自由流線理論

    自由流線理論最早由Helmholtz 和Kirchhoff針對(duì)二維平板阻力問題提出。之后,經(jīng)過Levi-Civita[8]及Roshko[9]等人的發(fā)展,使自由流線理論能夠適用于一般二維曲面升阻力的計(jì)算。Wu[10]基于上述理論,提出了一種非零空化數(shù)條件下的自由流線理論,并給出了圖1 模型對(duì)應(yīng)的阻力D和升力L計(jì)算公式

    圖1 非零空化數(shù)自由流線模型[10]Fig.1 Free streamline model under non-zero cavitation number condition

    式中:b1、b2與A和B這2 點(diǎn)的復(fù)勢(shì)有關(guān);ρ為液體密度;b、β為與b1和b2有關(guān)的參數(shù);σ為空化數(shù);ε為與空化數(shù)有關(guān)的參數(shù);A1、A2和A3為和二維曲面有關(guān)的系數(shù)。

    圖1 為非零空化數(shù)自由流線模型。圖中,流場(chǎng)復(fù)勢(shì)W=qeiθ,q為速度幅值,θ為速度方向。AD,BD′處速度幅值歸一化后,可得相應(yīng)的無窮遠(yuǎn)處速度DE和D′E′為無限長(zhǎng)空泡的平行段,速度幅值q由1 降低到無窮遠(yuǎn)速度U。

    由于自由流線理論不對(duì)邊界條件進(jìn)行線性化處理,難以通過數(shù)學(xué)分析應(yīng)對(duì)復(fù)雜曲面邊界,因而在實(shí)際應(yīng)用中受限。此外,通過自由流線模型計(jì)算的空泡尾部閉合處存在駐點(diǎn),與其空泡邊界處絕對(duì)速度恒定的假設(shè)相悖。為此后續(xù)衍生出了多種尾流模型,包括Riabouchinsky[11]提出的對(duì)稱模型、Tulin[12]提出的單渦管及雙渦管模型(如圖2(a)所示)、Gilbarg[13]提出的回射流模型 (如圖2(b)所示)和Zhukovsky[14]提出的開式模型等。圖2 中:AB為水下平板;α為其傾角;O為平板上的駐點(diǎn)。雙渦管模型中,C、E為閉合渦管上的2 點(diǎn),渦管上速度由v0變?yōu)闊o窮遠(yuǎn)處速度v∞;在回射流模型中,閉合渦管上的2 點(diǎn)C、E分別表示回射流和空泡尾部駐點(diǎn)。

    圖2 空泡尾部不同閉合方式Fig.2 Different closure schemes on cavity tail

    1.2 線性理論

    為了解決自由流線理論在數(shù)學(xué)分析方面的難題,Tulin 等[15]針對(duì)二維小攻角水翼超空泡流動(dòng)問題,將物面邊界條件線性化,最先提出了線性理論(模型如圖3 所示),U∞為無窮遠(yuǎn)處速度,水翼前后緣坐標(biāo)分別為x=0 和x=s,水翼縱坐標(biāo)由函數(shù)y0表示,空泡上壁面和下壁面縱坐標(biāo)分別由ycu和ycl表示。模型升力L和阻力D計(jì)算公式如下

    圖3 線性理論模型Fig.3 Linear theory model

    式中,p為當(dāng)?shù)仂o壓。

    Geurst[16-17]隨后采用線性理論分別研究了局部及完全空化水翼。Auslaender[18]應(yīng)用線性空泡流理論,給出了在特定表面壓力分布條件下,近自由面超空泡水翼的升阻力及空泡形狀,圖4 為其理論模型示意圖。圖中:p∞為無窮遠(yuǎn)處壓力;h為水翼距自由液面深度;yc為空泡壁面坐標(biāo);自由液面及空泡壁面處的速度u(x)=0。

    圖4 Auslaender 的近自由面水翼模型[18]Fig.4 Auslaender’s near free surface hydrofoil model

    Parkin 等[19]進(jìn)一步研究了不同壓力分布下超空泡水翼的阻力特性。線性理論雖然解決了復(fù)雜二維曲面求解的問題,應(yīng)用范圍較自由流線理論更加廣泛,但是需要先對(duì)物面條件線性化,且當(dāng)空泡長(zhǎng)度接近物體長(zhǎng)度時(shí),計(jì)算值與試驗(yàn)結(jié)果差異明顯。董世湯[20]基于Geurst 的研究成果,將尾流模型修改為空泡末端不封閉的開式模型,解決了原模型在空泡長(zhǎng)度接近弦長(zhǎng)時(shí)的失效問題。

    1.3 獨(dú)立膨脹原理

    獨(dú)立膨脹原理最先由Logvinovich[21]提出,其認(rèn)為空泡各個(gè)橫截面相對(duì)于航行體軌跡的膨脹過程獨(dú)立于航行體之前或之后的運(yùn)動(dòng),只與無窮遠(yuǎn)處壓力及空泡內(nèi)部壓力、航行體速度、航行體穿過截面時(shí)刻的航行體尺寸及阻力有關(guān)。

    Logvinovich 推導(dǎo)的空泡截面擴(kuò)展方程為

    式中:k=k(σ)為系數(shù),微弱依賴于空泡數(shù),通常取常數(shù)k=4π/A,A≈2為經(jīng)驗(yàn)值;l為空泡長(zhǎng)度;S為空泡橫截面積;Δp為空泡內(nèi)外壓差;xh為空泡截面流向坐標(biāo)。

    之后,Serebrayakov[22]結(jié)合細(xì)長(zhǎng)體理論與獨(dú)立膨脹原理,進(jìn)行了非定??张萁孛娣匠痰挠?jì)算。Paryshev[23-24]引入了新的參數(shù)用以表征時(shí)間延遲,針對(duì)不同類型的空泡建立了線性穩(wěn)定性和震蕩理論,并提出了一種考慮回射流的空泡尾部閉合模型。張學(xué)偉等[25]基于獨(dú)立膨脹原理發(fā)展了定常和非定常自然及通氣超空泡計(jì)算方法。宋書龍等[26]通過計(jì)算各個(gè)截面的垂直方向,結(jié)合Riabouchinsky尾流模型,建立了任意軌跡三維超空泡形態(tài)計(jì)算方法。獨(dú)立膨脹原理區(qū)別于前2 種基于復(fù)變函數(shù)的二維方法,能夠進(jìn)行三維、定常及非定常超空泡形態(tài)計(jì)算,適用范圍更廣、具有一定的優(yōu)勢(shì)。

    Paryshev 得到的空泡截面方程為

    式中: τ=t-t1,t1是空化器穿過截面h的時(shí)刻;S0是空化器截面積;是空化器最大截面積處空泡的初始擴(kuò)張速度。

    理論模型方法能夠?qū)Τ张莸耐庑芜M(jìn)行快速計(jì)算,但由于其簡(jiǎn)化程度較大,難以對(duì)超空泡流中非線性現(xiàn)象有效求解。隨著計(jì)算機(jī)技術(shù)的發(fā)展,近年來多采用對(duì)Navier-Stokes 方程進(jìn)行空間和時(shí)間離散求解的計(jì)算流體力學(xué)(computational fluid dynamics,CFD)方法對(duì)超空泡問題進(jìn)行研究。

    2 數(shù)值仿真方法

    目前主要的流體動(dòng)力學(xué)計(jì)算方法包括有限體積法(finite volume method,FVM)、邊界元法(boundary element method,BEM)、光滑粒子流體動(dòng)力學(xué)方法(smoothed particle hydrodynamics,SPH)和格子玻爾茲曼方法(lattice Boltzmann method,LBM)。對(duì)于水下超空泡流動(dòng)問題,基于勢(shì)流理論的BEM 難以計(jì)算湍流影響[27-28];SPH 在求解帶自由液面多相流問題具有優(yōu)勢(shì),但求解自然空化問題時(shí)存在不足[29];LBM 是一種對(duì)于多相流問題具有天然優(yōu)勢(shì)的介觀尺度仿真方法,但目前與空化有關(guān)的研究主要集中于潰滅過程和游離泡群而非超空泡[30-33];FVM 因其通用性和豐富的物理模型在超空泡研究中被廣泛應(yīng)用,是目前的主流仿真方法。

    水下超空泡這一多相湍流問題的仿真涉及多相流模型、界面捕捉方法、湍流模型、熱力學(xué)模型、相變模型和聲學(xué)模型等。采取不同的模型和界面捕捉方法對(duì)預(yù)報(bào)的空泡形態(tài)、尾部泄氣方式、升阻力特性及流噪聲預(yù)報(bào)具有重要影響。其中,以多相流和湍流模型影響最為顯著[34]。

    2.1 多相流模型

    目前,采用FVM 的超空泡流動(dòng)仿真一般采用分散多相流模型,其中最常用的模型包括混合(Mixture)多相流模型、流體體積(volume of fluid,VOF)多相流模型以及非均質(zhì)歐拉-歐拉(Euler-Euler)多相流模型。Mixture 及VOF 模型采用均相流假設(shè),控制體單元內(nèi)混合相屬性由各相通過體積分?jǐn)?shù)加權(quán)平均得到,之后代入控制方程進(jìn)行求解。Mixture 模型假設(shè)流體中各相易混溶,適用于大范圍分散多相流的仿真[35],其混合相速度、密度和動(dòng)力黏度計(jì)算公式分別為

    式中:vm為混合相速度矢量;ρm為混合相密度;vi為i相速度矢量;αi為i相體積分?jǐn)?shù);ρi為i相密度;μm為混合相動(dòng)力黏度;μi為i相動(dòng)力黏度。

    相對(duì)于假設(shè)單元內(nèi)流體速度一致的VOF 模型,Mixture 模型允許單元內(nèi)各相存在速度差,但其捕捉空泡界面的能力相對(duì)較弱。相較VOF 和Mixture 模型,Euler-Euler 模型不采用均相流假設(shè),而是直接對(duì)屬性不同的各相及其對(duì)應(yīng)控制方程組分別求解,并包含相間相互作用仿真。

    在以上多相流模型的具體應(yīng)用中,陳鑫等[36]基于Mixture 模型開展了定常超空泡仿真,Javadpour 等[37]進(jìn)一步給出了錐形空化器航行體的升阻力系數(shù)經(jīng)驗(yàn)公式,并通過試驗(yàn)對(duì)比證明了仿真結(jié)果的準(zhǔn)確性。仲宵等[38]則基于Mixture 模型預(yù)報(bào)了錐形通氣航行體形成的超空泡內(nèi)部流場(chǎng)速度及渦量特征,研究了通氣空泡內(nèi)部射流區(qū)和回流區(qū)的流動(dòng)特性。Radzyuk 等[39]也采用相同方法研究了約束來流條件下航行體尾部超空泡形態(tài),并說明Mixture 模型能夠很好地預(yù)報(bào)尾部空泡特征。Gopala 等[40]則采用Mixture 模型研究了不同空化器形狀的影響。劉如石等[41]采用Mixture 模型研究了不同尾部結(jié)構(gòu)超空泡射彈的尾拍運(yùn)動(dòng)特性,發(fā)現(xiàn)射彈尾部形狀的差異會(huì)導(dǎo)致沾濕面產(chǎn)生的空泡形態(tài)不同,進(jìn)而影響彈體載荷特性。Thang[42]和Kim[43]等也都基于Mixture 模型預(yù)報(bào)了超空泡形態(tài)。上述研究表明Mixture 模型能夠?qū)Τ张輧?nèi)部的流場(chǎng)結(jié)構(gòu)、空泡尾部閉合和航行體升阻力進(jìn)行有效仿真,但無法克服相間混溶引起界面捕捉精度不足的問題。

    VOF 模型適用于由明顯交界面分割的不混溶多相流體[44]。它假設(shè)每個(gè)多相混合流體單元內(nèi)部各相具有相同的速度和壓力,其單元混合密度和混合動(dòng)力黏度與Mixture 模型一樣可通過式(11)和式(12)計(jì)算,單元混合相滿足N-S 方程?;赩OF模型和相界面壓縮界面捕捉技術(shù),Roohi 等[45]計(jì)算了自然空化水翼在云狀空化和超空化狀態(tài)下的流動(dòng),并通過空化動(dòng)態(tài)特性和水翼升力系數(shù)評(píng)估了計(jì)算的有效性。Pendar 等[46]基于類似模型仿真了半球頭航行體自然空化,并指出了空泡長(zhǎng)度和直徑之間的關(guān)系,討論了邊界層分離和回射流對(duì)空泡閉合和脫落的影響。Jiang 等[47]通過修改VOF模型中液相屬性,成功預(yù)報(bào)了表面活性減阻劑對(duì)通氣超空泡流的影響,發(fā)現(xiàn)減阻劑在不同流速下對(duì)通氣超空泡形態(tài)和阻力會(huì)產(chǎn)生不同的影響。對(duì)于自由液面的影響,向敏等[48]等通過VOF 仿真了通氣超空泡與自由面空氣卷吸耦合作用,討論了自由面對(duì)空泡形態(tài)、航行體水動(dòng)力參數(shù)的影響。在上述研究中,VOF 模型均表現(xiàn)出比Mixture 模型更強(qiáng)的超空泡相界面捕捉能力[49],Mixture 模型[40]和VOF 模型[50]計(jì)算的二維錐形體的空泡形態(tài)對(duì)比如圖5 所示。

    圖5 Mixture 與VOF 多相流模型對(duì)比Fig.5 Comparison between Mixture model and VOF model

    Euler-Euler 模型區(qū)別于上述2 種均相流模型,在控制方程組中各相分別擁有獨(dú)立的連續(xù)性方程、能量方程以及動(dòng)量方程。通過這種方法,能夠針對(duì)每一相設(shè)置獨(dú)立的計(jì)算模型,各相之間的相互作用通過相界面?zhèn)鬟f。

    Euler-Euler 模型的連續(xù)性方程為

    式中:γα為 α相體積分?jǐn)?shù);ρα為 α相密度;Uα為 α相速度矢量;SMSα為質(zhì)量源項(xiàng);為單位體積內(nèi)從β相向 α相的質(zhì)量流量。

    通過基于Euler-Euler 模型的數(shù)值仿真,周景軍等[51]獲得了不同航速及通氣率條件下典型主動(dòng)通氣超空泡航行體的泄氣規(guī)律。楊明等[52]基于Euler-Euler 模型建立了數(shù)值水洞,指出水洞洞壁及阻塞效應(yīng)會(huì)減小通氣超空泡尺寸。Zhou 等[53]采用此模型研究了楔形艏舵對(duì)超空泡形態(tài)的影響,發(fā)現(xiàn)舵后方形成的空泡和主空泡融合后,改變了空泡整體形態(tài),Euler-Euler 模型能夠捕捉復(fù)雜結(jié)構(gòu)下的相間交界面。Xu 等[54]采用同樣方法研究了通氣超空泡的參數(shù)化影響特征及非定常特性。上述研究中,相較于均相流模型,Euler-Euler 模型也能較好地捕捉超空泡的氣水交界面,如圖6 所示。

    圖6 Euler-Euler 模型計(jì)算的超空泡形態(tài)等值面圖[53]Fig.6 Iso-surface of supercavity morphology simulated by Euler-Euler model

    2.2 湍流模型

    湍流模型作為影響超空泡數(shù)值仿真的另一項(xiàng)重要影響因素,目前主要包括雷諾時(shí)均仿真(Reynolds averaged navier-stokes simulation,RANS)、大渦仿真(large eddy simulation,LES)和分離渦仿真(detached-eddy simulation,DES)。RANS 方法 是 對(duì)Reynolds 方程的時(shí)均化求解,具有計(jì)算量低、魯棒性強(qiáng)的優(yōu)點(diǎn),同時(shí)在超空泡主體形態(tài)仿真、航行體力學(xué)特性等方面具有良好的表現(xiàn)。例如,周景軍等[55]利用k-ε模型討論了不同空化器錐角對(duì)超空泡形態(tài)和水動(dòng)力特性的影響,提出了一種圓錐空化器和圓盤空化器組合的方式來增大升力系數(shù)的方法。Park等[56]基于RANS 模型計(jì)算了不同楔形角的二維空化器超空化流動(dòng)現(xiàn)象,并研究了空化數(shù)對(duì)空泡長(zhǎng)度的影響。Choi 等[57]則研究了超空泡航行體在不同雷諾數(shù)和空化數(shù)下的空泡形態(tài)及阻力,發(fā)現(xiàn)阻力系數(shù)與雷諾數(shù)關(guān)系不大。Lu 等[58]采用realizablek-ε方法,研究了不同攻角條件下的空泡形態(tài)和阻尼力特性,如圖7 所示。王威等[59]通過RNGk-ε方法指出自然超空泡向通氣超空泡演化過程中,總阻力會(huì)進(jìn)一步小幅度減小,說明通氣超空泡能改善航行體的阻力性能。

    圖7 RANS 方法預(yù)報(bào)的不同速度下的超空泡形態(tài)等值面圖[58]Fig.7 Iso-surface of supercavity morphology with different flow speed using RANS method

    不同于RANS 的時(shí)域平均,LES 方法通過空間濾波函數(shù)區(qū)別大尺度和小尺度渦結(jié)構(gòu)并分別求解,在空泡尾部泄氣、空泡噪聲等重要瞬態(tài)特性仿真時(shí)具有更高的求解精度。通過LES 方法,趙怡等[60]分析了水下超高速航行體在不同速度下超空泡形態(tài)的變化規(guī)律,并對(duì)比了空泡脫落過程的差異。Pendar 等[61-62]先后研究了球體、半球體、圓盤和圓錐空化器的超空泡流動(dòng),與試驗(yàn)結(jié)果對(duì)比后認(rèn)為相較于RANS 方法,LES 方法可以獲得更準(zhǔn)確的數(shù)值結(jié)果。張木等[63]通過LES 方法對(duì)不同空化器的空泡初生問題進(jìn)行了研究,并指出圓盤空化器的空泡初生于空化器頂端,具有良好的超空泡生成能力。Wang 等[64]基于均質(zhì)多相流模型和LES 方法預(yù)報(bào)了如圖8 所示的超空泡尾部閉合現(xiàn)象,討論了主動(dòng)通氣對(duì)動(dòng)態(tài)超空泡渦流形成、演化和脫落的影響。

    圖8 LES 模型預(yù)報(bào)超空泡尾部渦管及脫落現(xiàn)象[64]Fig.8 Tail vortex and shedding of supercavity using LES model

    DES 作為一種混合方法,兼顧了空泡湍流模擬的RANS 特征與LES 特征,能夠?qū)张菸膊炕厣淞骶_仿真,但也繼承了上述湍流模型各自的局限。Kunz 等[65]通過DES 與RANS 的對(duì)比,發(fā)現(xiàn)DES 展現(xiàn)了對(duì)空泡脫落過程的更精細(xì)仿真。Kinzel 等[66]也認(rèn)為DES 方法能捕獲更廣尺度的湍流,并基于DES 方法研究了通氣超空泡與射流的動(dòng)態(tài)相互作用[67]。Gao 等[68]采用DES 模型研究了通氣率對(duì)空泡形態(tài)的影響,發(fā)現(xiàn)空泡形態(tài)變化滯后于通氣率的變化。Zhou 等[69]對(duì)通氣超空泡內(nèi)部流場(chǎng)結(jié)構(gòu),采用DES 模型并結(jié)合水洞試驗(yàn)結(jié)果進(jìn)行對(duì)比,證明了這種模型能夠有效捕捉空泡邊界的波動(dòng)。Choe[70]對(duì)帶控制鰭的高速水下航行體的超空泡流動(dòng)進(jìn)行了計(jì)算,著重研究了在不同自由來流速度、通氣量和攻角時(shí)的水動(dòng)力特性,并捕捉到了如圖9 所示空泡壁面的非定常波動(dòng)現(xiàn)象。

    圖9 DES 模型相體積分?jǐn)?shù)等值面預(yù)報(bào)的尾鰭附近空泡壁面波動(dòng)[70]Fig.9 Cavity fluctuation near the tail fin predicted by phase volume fraction iso-surface using DES model

    總體而言,上述3 種湍流模型在超空泡數(shù)值仿真中均被廣泛采用。RANS 將控制方程進(jìn)行了時(shí)間平均,更適用于有限計(jì)算資源下超空泡整體形狀及航行體力學(xué)特性計(jì)算,但側(cè)重提供湍流的平均信息,對(duì)于預(yù)報(bào)超空泡尾部閉合及渦生長(zhǎng)脫落動(dòng)態(tài)過程精度不足。LES 通過對(duì)大尺度渦的直接求解,能夠捕捉到RANS 難以預(yù)報(bào)的高頻信息,在分離流動(dòng)及瞬態(tài)特性計(jì)算方面具有優(yōu)勢(shì),適用于超空泡尾部閉合及脫落問題的研究,但相對(duì)需要更多的計(jì)算資源。DES 結(jié)合上述2 種方法,在近壁面處采用RANS 而在自由流動(dòng)區(qū)域采用LES,降低了整體計(jì)算資源需求,同時(shí)在適當(dāng)?shù)木W(wǎng)格分辨率下也能夠較好地捕捉遠(yuǎn)離壁面位置超空泡流動(dòng)的非定常特征。

    2.3 空化模型

    在高速條件下,航行體肩部繞流區(qū)或尾部等低壓區(qū)會(huì)引起局部空化或自然超空泡現(xiàn)象,為了描述這一相變過程,需要在數(shù)值計(jì)算中引入空化模型??张輨?dòng)力學(xué)中的Rayleigh-Plesset 方程(R-P 方程)描述了液體中單個(gè)氣泡的生長(zhǎng)和潰滅過程,并為大多數(shù)空化模型奠定了理論基礎(chǔ),其方程如下

    式中:R為氣泡直徑: υL為氣泡周圍流體的運(yùn)動(dòng)黏度;γ為表面張力;pB為氣泡內(nèi)部壓力;ρL為液相密度。

    Kubota 等[71-72]提出了一種氣泡兩相流模型,將空泡流場(chǎng)視為密度變化較大的可壓縮黏性流體,描述了含大尺度渦的黏性流場(chǎng)與空泡之間的相互作用。區(qū)別于Kubota 模型,Merkle 模型[73]的推導(dǎo)建立在氣泡群而非單個(gè)氣泡,由混合密度導(dǎo)出相間質(zhì)量傳輸率。Kunz 等[74]提出了一種求解黏性兩相流的半解析模型,獲得了空泡形成過程中空化數(shù)、阻力系數(shù)等參數(shù)的變化規(guī)律。Singhal 等[75]針對(duì)繞二維水翼空化流,推導(dǎo)出了蒸汽與液體之間的相變率,提出了“全空化模型”。Schnerr 等[76-77]對(duì)R-P 方程進(jìn)行了簡(jiǎn)化,忽略了黏性效應(yīng)和表面張力的影響,并結(jié)合試驗(yàn)結(jié)果準(zhǔn)確預(yù)測(cè)了空泡生長(zhǎng)和潰滅,方程如下

    式中: αv為蒸汽相體積分?jǐn)?shù);ρv和 ρl分別為蒸汽相和液相密度;uj為流體速度的j方向分量;和分別為凝結(jié)源項(xiàng)質(zhì)量傳輸率和蒸發(fā)源項(xiàng)質(zhì)量傳輸率;。

    Owis 等[78]利用質(zhì)量分?jǐn)?shù)形式,計(jì)算了空泡魚雷上的多相流,計(jì)算結(jié)果與實(shí)測(cè)結(jié)果吻合良好。Gerber[79]建立了廣泛兩相空化流模型,將空化模型推廣到了更高的密度比。Zwart 等[80]在Kubota和Gerber 的基礎(chǔ)上,修正了質(zhì)量空化率方程中的蒸汽體積分?jǐn)?shù)項(xiàng)。

    2.4 聲學(xué)模型

    作為水下航行體的一項(xiàng)重要特性,超空泡航行的噪聲主要來源包括空泡周圍湍流壓力脈動(dòng)產(chǎn)生的噪聲、空泡氣-水界面處壓力波動(dòng)產(chǎn)生的噪聲和通氣氣體射流直接撞擊在空泡壁面產(chǎn)生的噪聲。對(duì)于噪聲問題,一般有解析求解、聲學(xué)類比方程及邊界積分等求解方法。Howe 等[81]聚焦通氣氣體直接撞擊空泡氣-水界面產(chǎn)生的噪聲,提出了一種將航行體頭部自噪聲壓力和空泡內(nèi)壁壓力分布相聯(lián)系的傳遞函數(shù),并準(zhǔn)確預(yù)測(cè)了氣體射流撞擊產(chǎn)生的自噪聲水平。張風(fēng)珍等[82]指出空泡尾部細(xì)小空泡的潰滅是超空泡航行體空化噪聲的主要來源,基于單空泡潰滅輻射噪聲方程并結(jié)合瑞麗分布,建立了空泡群潰滅輻射噪聲模型,分析了不同空泡數(shù)目和平均半徑條件下的空化噪聲波形和頻譜特性。Sun 等[83]采用FW-H 聲學(xué)類比法研究了調(diào)制通氣對(duì)超空泡聲學(xué)性能的影響,發(fā)現(xiàn)可以通過改變調(diào)制頻率的方法改變?cè)肼曋黝l。Ramesh等[84-85]通過邊界積分法對(duì)超空泡噪聲的傳播進(jìn)行了研究,發(fā)現(xiàn)相較于作用在超空泡表面的偶極子聲源,超空泡附近的湍流壓力脈動(dòng)產(chǎn)生的噪聲更大。Skidmore 等[86]基于單極子源輻射聲壓表達(dá)式仿真了脈動(dòng)超空泡腔內(nèi)壓力和近場(chǎng)輻射噪聲,與水洞試驗(yàn)結(jié)果進(jìn)行對(duì)比,證明了數(shù)值方法的準(zhǔn)確性。

    在上述超空泡航行體噪聲問題的研究中,基于公式推導(dǎo)的解析方法需要對(duì)模型進(jìn)行一定簡(jiǎn)化,一般只能分析一種噪聲源產(chǎn)生的噪聲。FW-H 聲學(xué)類比方程對(duì)3 個(gè)噪聲源項(xiàng)都明確定義了其物理意義,有助于分析噪聲的產(chǎn)生機(jī)理,同時(shí)能直觀地顯示空泡噪聲與流場(chǎng)動(dòng)態(tài)的關(guān)聯(lián),對(duì)于超空泡問題較為適用。相較于有限元方法,采用邊界積分方法仿真聲傳播問題能顯著提高計(jì)算效率,但對(duì)于超空泡內(nèi)部流動(dòng)的捕捉不夠精細(xì),導(dǎo)致在某些頻率無法求解正確的聲壓值。

    3 超空泡計(jì)算的典型影響因素

    3.1 自然空化數(shù)

    水下高速航行體會(huì)發(fā)生自然空化現(xiàn)象。當(dāng)空化數(shù)足夠小時(shí),自然空化形成的空泡能夠完全包裹航行體,進(jìn)而在不需要主動(dòng)通氣的情況下形成自然超空泡。

    在不同的自然空化數(shù)下,Passandideh-Frad 等[87]通過數(shù)值仿真指出楔形及圓盤空化器產(chǎn)生的超空泡長(zhǎng)度隨空化數(shù)減小而增長(zhǎng),Gugulothu[88]也給出了類似研究結(jié)果。金大橋等[89]研究了空化數(shù)對(duì)水下超空泡減阻特性的影響,發(fā)現(xiàn)隨著空化數(shù)的減小,壓差阻力系數(shù)、摩擦阻力系數(shù)和總阻力系數(shù)均減小,直至生成超空泡后基本不再變化。Jiang等[90]在不同自然空化數(shù)條件下分析了圓柱體的阻力系數(shù),同樣發(fā)現(xiàn)小空化數(shù)時(shí)圓柱阻力系數(shù)較小。張木等[91]研究了帶尾翼水下自然超空泡射彈的空泡形態(tài),發(fā)現(xiàn)空化器直徑一定時(shí),超空泡長(zhǎng)度和直徑隨空化數(shù)減小而增大。Kadivar 等[92]進(jìn)一步研究了不同楔形角錐形空化器的自然超空化流動(dòng),發(fā)現(xiàn)小角度空化器對(duì)空化數(shù)變化不敏感,而大楔形角空化器尾部超空泡長(zhǎng)度隨空化數(shù)減小迅速增大。從上述研究可以看出,與Logvinovich 提出的理論模型類似,無論是空化器還是鈍體尾部自然空化形成的超空泡,其形態(tài)均受到自然空化數(shù)直接且顯著的影響。

    3.2 通氣參數(shù)

    上述自然超空化需要高航速以達(dá)到所需的極小空化數(shù)。在實(shí)際應(yīng)用中,大部分航行體通過主動(dòng)通氣形成完全包裹航行體的穩(wěn)定空泡,即通氣超空泡。對(duì)于通氣超空泡而言,通氣參數(shù)是其重要影響參數(shù)之一。通氣率大小、通氣角度及周期性通氣等都會(huì)明顯影響空泡形態(tài)及空泡尾部閉合模式,進(jìn)一步影響超空泡航行體的受力情況[93-96]。

    關(guān)于通氣率與空泡形態(tài),楊武剛等[97]通過采用Mixture 模型的數(shù)值仿真,發(fā)現(xiàn)通氣率是影響空泡形態(tài)及穩(wěn)定性的最主要因素,而同一通氣率下的空化器部位通氣孔開孔數(shù)量對(duì)空泡形態(tài)影響不大。黃海龍等[98]研究了通氣角度的影響規(guī)律,發(fā)現(xiàn)隨著通氣角度的增加,空泡長(zhǎng)度減小,直徑增大。Rashidi 等[99]研究了不同通氣率下通氣超空泡的形態(tài)及內(nèi)部壓力,給出了空泡長(zhǎng)度、通氣率及空化數(shù)的關(guān)系。Nouri[100]和Kim[101]等的研究表明通氣率的增加會(huì)導(dǎo)致超空泡直徑和長(zhǎng)度增加。Cao 等[102]進(jìn)一步指出超空泡內(nèi)部壓力主要受靜水壓力影響,但是增加通氣率會(huì)增加空泡長(zhǎng)度,導(dǎo)致空泡閉合位置的逆壓梯度減小。Fronzeo 等[103]則發(fā)現(xiàn)通氣率變化后空泡內(nèi)部產(chǎn)生了多個(gè)高壓及低壓區(qū)域,表明在流動(dòng)細(xì)節(jié)分析中不適用空泡內(nèi)部恒壓假設(shè)。Wang 等[104]研究發(fā)現(xiàn),通氣率增加會(huì)導(dǎo)致如圖10 所示空泡內(nèi)部低速區(qū)擴(kuò)張、空泡壁面處渦量升高、渦帶拉伸更加明顯且規(guī)律。

    圖10 不同通氣率下空泡內(nèi)部渦量云圖[104]Fig.10 Vorticity field inside cavity under different ventilation rates

    Lindau 等[105]通過對(duì)通氣率增加正弦分量進(jìn)行調(diào)制,發(fā)現(xiàn)調(diào)制通氣能夠減小超空泡的脈動(dòng)情況。Sun 等[83]發(fā)現(xiàn)調(diào)制通氣可以改變空泡的脫落周期,并促進(jìn)發(fā)卡渦的形成及生長(zhǎng),同時(shí)能夠小幅提高航行體的升阻力性能,如圖11 所示。An 等[106]發(fā)現(xiàn)通氣口位置對(duì)于航行體的表面摩擦阻力影響顯著。Lv 等[107]則發(fā)現(xiàn)通氣率增加后,閉合模式由回射流轉(zhuǎn)變?yōu)殡p渦管,并認(rèn)為通氣增加后的流動(dòng)分離減小和超空泡不對(duì)稱性增加是導(dǎo)致這種轉(zhuǎn)變的重要原因。Pham 等[108]研究了高溫通氣對(duì)尾部閉合模式的影響,發(fā)現(xiàn)高溫通氣條件下閉合模式為雙渦管,但是提高來流速度后轉(zhuǎn)變?yōu)閱螠u管模式。

    圖11 不同調(diào)制頻率空泡尾部渦結(jié)構(gòu)[83]Fig.11 Vortex at cavity tail under different modulated frequencies

    總的來說,通氣參數(shù)的變化可影響空泡長(zhǎng)度及尾部閉合模式。其中,通氣率是較為重要的影響參數(shù),在一定范圍內(nèi)通氣率增加會(huì)使空泡長(zhǎng)度增加、尾部逆壓梯度減小,閉合模式也由回射流轉(zhuǎn)變?yōu)殡p渦管模式。

    3.3 來流條件

    不同的來流條件會(huì)對(duì)超空泡形態(tài)產(chǎn)生明顯的影響,其中弗勞德數(shù)Fr、來流攻角及周期性來流等條件對(duì)超空泡流動(dòng)特征和航行體升阻力特性具有重要研究意義[109-111]。

    Kim 等[112]仿真了不同來流速度和深度下超空泡航行體的超空泡形態(tài)和阻力特性,給出了阻力隨速度的變化曲線。Wang 等[113]在恒定通氣率條件下,研究了不同來流速度超空泡尾部泄氣機(jī)理,發(fā)現(xiàn)低Fr條件下,空泡尾部呈雙渦管模式,而高Fr條件下,空泡尾部呈周期性環(huán)狀脫落。Karn等[114]進(jìn)一步發(fā)現(xiàn),隨著來流速度的變化,還存在四渦管、脈動(dòng)及泡沫流等空泡尾部閉合模式。Jafarian 等[115]發(fā)現(xiàn)來流速度超過一定限度后,基本不再對(duì)空泡產(chǎn)生顯著影響。王威等[116]利用VOF模型及k-ε湍流模型,在航行體前部設(shè)置陣風(fēng)發(fā)生器產(chǎn)生周期性來流,發(fā)現(xiàn)周期性來流會(huì)使航行體的沾濕區(qū)域發(fā)生變化,進(jìn)一步導(dǎo)致空泡形態(tài)及升阻力產(chǎn)生如圖12所示的周期性變化。他們隨后研究了超空泡航行體在一定側(cè)滑角條件下的空泡形態(tài),給出了能夠減小沾濕面積的側(cè)滑角范圍[117]。Huang 等[118]進(jìn)一步研究了不同陣風(fēng)攻角及陣風(fēng)頻率下超空泡形態(tài)及壓力的變化情況,發(fā)現(xiàn)攻角及頻率的增加均會(huì)增強(qiáng)空泡波動(dòng)及泡內(nèi)壓力脈動(dòng)。Zou 等[119]則發(fā)現(xiàn)攻角會(huì)影響超空泡脫落頻率,并給出了脫落頻率與攻角的關(guān)系。Sun 等[120]研究了內(nèi)波對(duì)超空泡形態(tài)的影響,發(fā)現(xiàn)在周期性內(nèi)波的影響下,空泡脫落和壓力脈動(dòng)強(qiáng)度增加,脫落頻率與內(nèi)波頻率相關(guān)。Moltani 等[121]研究了近自由面的超空泡流動(dòng),發(fā)現(xiàn)空泡長(zhǎng)度隨來流速度的變化與遠(yuǎn)自由面區(qū)域不同,并且阻力系數(shù)隨深度增加而增大。

    圖12 由相體積分?jǐn)?shù)等值面表示的周期性來流對(duì)超空泡影響[116]Fig.12 Influence of periodic flow on supercavity shown by iso-surface of phase volume fraction

    來流速度同樣會(huì)影響超空泡尾部的閉合模式。許海雨等[122]研究發(fā)現(xiàn),Fr增加后,超空泡尾部由雙渦管模式轉(zhuǎn)變?yōu)榛厣淞髂J?空泡內(nèi)部壓力震蕩特性增強(qiáng);雙渦管閉合方式下,存在空泡表面及雙渦管2 條泄氣途徑,且兩者之間的比例只與Fr相關(guān),Fr增加則從表面泄漏氣體增加;回射流閉合方式下,氣體主要由空泡表面流出。

    3.4 幾何構(gòu)型

    超空泡航行體整體布局一般包括空化器、擴(kuò)張段、平行中體、尾部收縮段和控制舵翼。不同的形狀設(shè)計(jì)和幾何尺寸會(huì)明顯改變超空泡的形態(tài)及航行體阻力,這一問題也是近期國(guó)內(nèi)外研究關(guān)注的熱點(diǎn)之一[123-125]。

    陳鑫等[126]基于VOF 模型,對(duì)二維航行體通氣超空泡進(jìn)行了數(shù)值仿真,結(jié)果表明空化器設(shè)計(jì)及通氣流量需要與環(huán)境壓力相匹配,才能形成穩(wěn)定的空泡外形。王海斌等[127]采用均質(zhì)平衡流模型,同樣得到了空化器設(shè)計(jì)對(duì)超空泡幾何尺度存在明顯影響的結(jié)論。Jia 等[128]研究了不同空化器角度對(duì)空泡不對(duì)稱性及阻力特征,總結(jié)了空化器角度的影響規(guī)律。蔣增輝等[129]采用Mixture 模型對(duì)不同后體超空泡航行體的阻力及空泡形態(tài)進(jìn)行了預(yù)報(bào)。Ahn 等[130]研究了不同空化器構(gòu)型以及楔形空化器楔形角的影響。Mansour 等[131]則研究了不同頭型航行體超空泡形態(tài)和阻力系數(shù),如圖13 所示,并設(shè)計(jì)了新的頭部形狀,實(shí)現(xiàn)了更小的阻力系數(shù)。Prichard 等[132]通過數(shù)值計(jì)算結(jié)果,給出了空化數(shù)、阻力系數(shù)與空化器半徑之間的關(guān)系,并給出了特定航行體構(gòu)型下的最佳空化器直徑。

    圖13 不同頭部形狀對(duì)超空泡云圖形態(tài)影響[131]Fig.13 Influence of different head shapes on supercavity morphology

    Kim 等[133]在上述研究基礎(chǔ)上進(jìn)一步分析了空化器攻角a與空泡不對(duì)稱性之間的關(guān)系,發(fā)現(xiàn)無攻角條件下空泡主要受浮力作用影響而產(chǎn)生不對(duì)稱現(xiàn)象,而有攻角條件下浮力對(duì)空泡不對(duì)稱性產(chǎn)生的作用小于空化器的作用,如圖14 所示。綜合上述研究,幾何構(gòu)型主要影響了超空泡整體形態(tài)及升阻力系數(shù),但對(duì)于空泡閉合模式及脫落過程影響較小。

    圖14 不同空化器攻角下超空泡形態(tài)[133]Fig.14 Supercavity morphology under different cavitator attack angles

    4 結(jié)論

    文中綜述了國(guó)內(nèi)外進(jìn)行超空泡問題研究的數(shù)值計(jì)算方法及相關(guān)研究進(jìn)展。數(shù)值計(jì)算方法中理論模型求解方法包括自由流線理論、線性理論和獨(dú)立膨脹原理;數(shù)值仿真方法則以有限體積法為主,涉及多相流模型、空化模型、湍流模型和聲學(xué)模型。之后,圍繞著自然空化數(shù)、通氣參數(shù)、來流條件和幾何構(gòu)型4 個(gè)超空泡流動(dòng)計(jì)算重要影響因素,介紹了超空泡數(shù)值計(jì)算的研究進(jìn)展。得到的主要結(jié)論如下:

    1) 理論模型方法基于勢(shì)流理論求解超空泡形態(tài)及升阻力,但是勢(shì)流理論忽略了黏性影響,因而對(duì)摩擦阻力的預(yù)報(bào)存在偏差。解析方法中,自由流線理論和線性理論適用于二維問題,而獨(dú)立膨脹原理能夠完成三維非定常問題計(jì)算。

    2) 有限體積法是最為常用的水下超空泡數(shù)值仿真方法,涉及到多相流模型、空化模型、湍流模型以及聲學(xué)模型等。其中不同的多相流模型及空化模型對(duì)于特定問題可獲得較好的計(jì)算結(jié)果;湍流模型的選擇主要考慮對(duì)仿真結(jié)果時(shí)間和空間分辨率的需求,對(duì)于關(guān)注非定常流動(dòng)細(xì)節(jié)超空流的仿真,有必要采用求解精度更高的分離渦或大渦仿真等方法;聲學(xué)模型中聲學(xué)類比方程對(duì)超空泡流較為適用。

    3) 自然空化數(shù)、通氣參數(shù)、來流參數(shù)和幾何構(gòu)型是影響超空泡流動(dòng)特征的幾個(gè)主要參數(shù)。自然空化數(shù)降低往往增加相變強(qiáng)度;通氣參數(shù)則主要影響空泡的尺度及尾部閉合模式;來流速度主要影響尾部閉合模式,周期性來流則強(qiáng)迫空泡發(fā)生周期性波動(dòng);幾何構(gòu)型尤其是空化器構(gòu)型會(huì)對(duì)空泡整體形態(tài)產(chǎn)生重要影響。

    猜你喜歡
    模型研究
    一半模型
    FMS與YBT相關(guān)性的實(shí)證研究
    2020年國(guó)內(nèi)翻譯研究述評(píng)
    遼代千人邑研究述論
    重要模型『一線三等角』
    重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
    視錯(cuò)覺在平面設(shè)計(jì)中的應(yīng)用與研究
    科技傳播(2019年22期)2020-01-14 03:06:54
    EMA伺服控制系統(tǒng)研究
    新版C-NCAP側(cè)面碰撞假人損傷研究
    3D打印中的模型分割與打包
    亚洲精品一二三| 大香蕉97超碰在线| 亚洲婷婷狠狠爱综合网| 国精品久久久久久国模美| av卡一久久| 大片免费播放器 马上看| 久久6这里有精品| 最近中文字幕高清免费大全6| 免费观看a级毛片全部| 午夜福利高清视频| 亚洲精品aⅴ在线观看| 国产日韩欧美在线精品| 新久久久久国产一级毛片| 亚洲精品国产av蜜桃| 中文在线观看免费www的网站| 成人综合一区亚洲| 午夜免费观看性视频| 日韩一区二区视频免费看| 亚洲av不卡在线观看| 18禁动态无遮挡网站| 中文字幕人妻熟人妻熟丝袜美| 亚洲av中文字字幕乱码综合| 特级一级黄色大片| 人体艺术视频欧美日本| 亚洲欧美一区二区三区黑人 | 永久网站在线| 亚洲激情五月婷婷啪啪| 又粗又硬又长又爽又黄的视频| 一本久久精品| 日日啪夜夜爽| 亚洲av日韩在线播放| 免费电影在线观看免费观看| 丝袜美腿在线中文| 又爽又黄a免费视频| 看十八女毛片水多多多| www.色视频.com| 男男h啪啪无遮挡| 女的被弄到高潮叫床怎么办| 网址你懂的国产日韩在线| 婷婷色综合大香蕉| 在线观看美女被高潮喷水网站| 国产精品国产三级国产专区5o| 人妻 亚洲 视频| 纵有疾风起免费观看全集完整版| 亚洲国产欧美人成| 韩国高清视频一区二区三区| 一个人看视频在线观看www免费| 91aial.com中文字幕在线观看| 久久影院123| 午夜激情福利司机影院| 亚洲性久久影院| av播播在线观看一区| 久久ye,这里只有精品| 秋霞在线观看毛片| 亚洲最大成人中文| 99久久精品热视频| 亚洲四区av| 国产成人a区在线观看| 亚洲电影在线观看av| 最近最新中文字幕免费大全7| 一个人观看的视频www高清免费观看| 中文字幕亚洲精品专区| 免费大片黄手机在线观看| 五月伊人婷婷丁香| 丰满少妇做爰视频| 日韩视频在线欧美| 1000部很黄的大片| 嫩草影院精品99| 秋霞在线观看毛片| 人妻制服诱惑在线中文字幕| 亚洲怡红院男人天堂| 97热精品久久久久久| 极品少妇高潮喷水抽搐| 国内精品美女久久久久久| 18禁裸乳无遮挡动漫免费视频 | 日韩av在线免费看完整版不卡| 一级毛片电影观看| 欧美成人精品欧美一级黄| 亚洲最大成人中文| 亚洲国产精品专区欧美| 男人和女人高潮做爰伦理| 成人综合一区亚洲| 日韩一本色道免费dvd| 亚洲精品国产色婷婷电影| 国产中年淑女户外野战色| 人人妻人人看人人澡| 免费高清在线观看视频在线观看| 少妇裸体淫交视频免费看高清| 成人鲁丝片一二三区免费| 国产探花极品一区二区| 久久99热这里只有精品18| 欧美成人一区二区免费高清观看| 国产一区有黄有色的免费视频| 天堂中文最新版在线下载 | 尤物成人国产欧美一区二区三区| 精品熟女少妇av免费看| 夫妻午夜视频| 国产一级毛片在线| 成人亚洲精品一区在线观看 | 久久久久久国产a免费观看| 最近的中文字幕免费完整| 亚洲精品一二三| 欧美日韩在线观看h| 国产黄片视频在线免费观看| 欧美成人午夜免费资源| 永久网站在线| 精品一区在线观看国产| 欧美3d第一页| av天堂中文字幕网| 高清av免费在线| 99热6这里只有精品| 久久久久久久久久成人| 亚洲伊人久久精品综合| 最近手机中文字幕大全| 国产探花在线观看一区二区| 午夜福利网站1000一区二区三区| 国产亚洲最大av| 国产69精品久久久久777片| 免费看光身美女| 日韩不卡一区二区三区视频在线| 日韩亚洲欧美综合| 欧美高清成人免费视频www| 亚洲欧美日韩无卡精品| 精品久久久精品久久久| 99久久精品一区二区三区| 国产极品天堂在线| 在线观看免费高清a一片| 亚洲天堂国产精品一区在线| 精品人妻偷拍中文字幕| 一区二区三区四区激情视频| 又黄又爽又刺激的免费视频.| 亚洲精品aⅴ在线观看| 国产高清国产精品国产三级 | 看免费成人av毛片| 一本色道久久久久久精品综合| 男人爽女人下面视频在线观看| 午夜免费男女啪啪视频观看| 在线观看一区二区三区| 日韩免费高清中文字幕av| 联通29元200g的流量卡| 亚洲精品成人久久久久久| 18禁在线无遮挡免费观看视频| 久久久久久久亚洲中文字幕| 只有这里有精品99| 性色avwww在线观看| 一级毛片aaaaaa免费看小| 一区二区三区精品91| 日日摸夜夜添夜夜爱| 国产大屁股一区二区在线视频| 三级国产精品欧美在线观看| 国产精品爽爽va在线观看网站| 国产v大片淫在线免费观看| 成人亚洲欧美一区二区av| 亚洲精品色激情综合| 亚洲精品成人av观看孕妇| 亚洲最大成人中文| 2018国产大陆天天弄谢| 一级av片app| av天堂中文字幕网| 99九九线精品视频在线观看视频| 日日啪夜夜爽| 国产极品天堂在线| 性色avwww在线观看| 日韩欧美一区视频在线观看 | 纵有疾风起免费观看全集完整版| 91在线精品国自产拍蜜月| 老女人水多毛片| 天天躁夜夜躁狠狠久久av| 天美传媒精品一区二区| 99热这里只有是精品50| 深夜a级毛片| 丝袜美腿在线中文| 午夜免费观看性视频| 精华霜和精华液先用哪个| 一级毛片我不卡| 亚洲人成网站在线观看播放| 免费av毛片视频| 午夜爱爱视频在线播放| 人人妻人人澡人人爽人人夜夜| 九色成人免费人妻av| 男人狂女人下面高潮的视频| 日本黄大片高清| 久久人人爽av亚洲精品天堂 | 熟女电影av网| 成人漫画全彩无遮挡| 伊人久久国产一区二区| 日韩av在线免费看完整版不卡| 欧美zozozo另类| 精品一区二区三卡| eeuss影院久久| 国内精品宾馆在线| 岛国毛片在线播放| 日韩强制内射视频| 国产大屁股一区二区在线视频| 亚洲精品日韩在线中文字幕| 国产成人免费无遮挡视频| 免费少妇av软件| 中文欧美无线码| 日本一二三区视频观看| 国产一区二区三区综合在线观看 | 亚洲av免费高清在线观看| 高清视频免费观看一区二区| 我的女老师完整版在线观看| 又爽又黄无遮挡网站| 国产有黄有色有爽视频| 亚洲精品久久午夜乱码| 久久久久久久亚洲中文字幕| 精品国产乱码久久久久久小说| 五月天丁香电影| 十八禁网站网址无遮挡 | 视频区图区小说| 麻豆成人午夜福利视频| 大陆偷拍与自拍| 综合色av麻豆| 中文在线观看免费www的网站| 免费高清在线观看视频在线观看| 欧美xxxx黑人xx丫x性爽| 波多野结衣巨乳人妻| 亚洲怡红院男人天堂| 免费黄色在线免费观看| 国产有黄有色有爽视频| 免费观看性生交大片5| 青春草视频在线免费观看| 欧美日韩精品成人综合77777| 亚洲av免费高清在线观看| 搞女人的毛片| 国产一区二区亚洲精品在线观看| 新久久久久国产一级毛片| 亚洲精品,欧美精品| 国语对白做爰xxxⅹ性视频网站| 国产精品国产av在线观看| 一区二区av电影网| 人妻系列 视频| 夜夜看夜夜爽夜夜摸| 青青草视频在线视频观看| 亚洲人成网站高清观看| 国产精品av视频在线免费观看| 在线免费观看不下载黄p国产| 欧美最新免费一区二区三区| 久久99热这里只频精品6学生| 在线观看一区二区三区激情| 少妇人妻一区二区三区视频| 亚洲欧美一区二区三区国产| 欧美xxⅹ黑人| 成人无遮挡网站| 91精品一卡2卡3卡4卡| 高清日韩中文字幕在线| 欧美极品一区二区三区四区| 亚洲成人精品中文字幕电影| 老司机影院成人| 国产精品.久久久| 激情五月婷婷亚洲| 精品一区二区三卡| 亚洲精品乱久久久久久| 国产久久久一区二区三区| 日产精品乱码卡一卡2卡三| 亚洲综合精品二区| 免费av不卡在线播放| 久久综合国产亚洲精品| 18+在线观看网站| 亚洲精品色激情综合| 中文字幕亚洲精品专区| 精品人妻视频免费看| 中文乱码字字幕精品一区二区三区| 亚洲欧美一区二区三区国产| 亚洲经典国产精华液单| av国产精品久久久久影院| 乱码一卡2卡4卡精品| 三级经典国产精品| 我的女老师完整版在线观看| 国产亚洲精品久久久com| 亚洲人成网站高清观看| 汤姆久久久久久久影院中文字幕| 成人无遮挡网站| 免费看光身美女| 成人特级av手机在线观看| 真实男女啪啪啪动态图| 七月丁香在线播放| 直男gayav资源| 亚洲自拍偷在线| 国产人妻一区二区三区在| 欧美区成人在线视频| av国产精品久久久久影院| 久久久欧美国产精品| 爱豆传媒免费全集在线观看| 亚洲av在线观看美女高潮| 日本与韩国留学比较| 简卡轻食公司| 日韩欧美 国产精品| 好男人视频免费观看在线| 高清视频免费观看一区二区| 青春草亚洲视频在线观看| 免费少妇av软件| 少妇人妻久久综合中文| 高清毛片免费看| 国产精品福利在线免费观看| 国产美女午夜福利| 久久热精品热| 亚洲国产色片| 人妻夜夜爽99麻豆av| 亚洲婷婷狠狠爱综合网| 亚洲成色77777| 69人妻影院| 中文资源天堂在线| 亚洲最大成人av| 国产淫语在线视频| 精品午夜福利在线看| 久久精品国产鲁丝片午夜精品| 亚洲天堂av无毛| 国产中年淑女户外野战色| 两个人的视频大全免费| 男人爽女人下面视频在线观看| 美女脱内裤让男人舔精品视频| 一二三四中文在线观看免费高清| 欧美激情久久久久久爽电影| 久久6这里有精品| 亚洲性久久影院| 狂野欧美激情性xxxx在线观看| 中文字幕久久专区| 一边亲一边摸免费视频| 一级毛片aaaaaa免费看小| 爱豆传媒免费全集在线观看| 在线观看一区二区三区激情| 亚洲av福利一区| 观看免费一级毛片| 亚州av有码| av在线app专区| 日韩伦理黄色片| 亚洲精品456在线播放app| 久久99精品国语久久久| 青春草国产在线视频| 亚洲精品乱码久久久久久按摩| 在线 av 中文字幕| 午夜福利网站1000一区二区三区| av女优亚洲男人天堂| 亚洲第一区二区三区不卡| 国模一区二区三区四区视频| 99视频精品全部免费 在线| 色播亚洲综合网| 久久鲁丝午夜福利片| 成人漫画全彩无遮挡| 18+在线观看网站| 久久精品综合一区二区三区| 另类亚洲欧美激情| 日韩强制内射视频| 亚洲无线观看免费| 少妇裸体淫交视频免费看高清| 久久亚洲国产成人精品v| 久久久精品94久久精品| 国产成人午夜福利电影在线观看| 国产亚洲最大av| 午夜视频国产福利| 舔av片在线| 午夜免费男女啪啪视频观看| 听说在线观看完整版免费高清| 国产成人午夜福利电影在线观看| 欧美性猛交╳xxx乱大交人| 看黄色毛片网站| 成人特级av手机在线观看| 久久这里有精品视频免费| 九色成人免费人妻av| 日韩电影二区| 欧美少妇被猛烈插入视频| 久久久精品免费免费高清| 久久久国产一区二区| 国产免费又黄又爽又色| 色哟哟·www| 成人一区二区视频在线观看| 夫妻午夜视频| 欧美日本视频| 国产黄片视频在线免费观看| 亚洲色图av天堂| 久久韩国三级中文字幕| 最近中文字幕2019免费版| 丰满少妇做爰视频| 亚洲欧美一区二区三区黑人 | 午夜福利视频精品| 国产男女超爽视频在线观看| 国产在视频线精品| 麻豆久久精品国产亚洲av| 2021少妇久久久久久久久久久| 日本-黄色视频高清免费观看| 亚洲欧美清纯卡通| 日本黄大片高清| 精品一区二区三区视频在线| 久久精品人妻少妇| 国产亚洲一区二区精品| 欧美日韩精品成人综合77777| 欧美日韩亚洲高清精品| 亚洲美女视频黄频| 国产 一区 欧美 日韩| 爱豆传媒免费全集在线观看| 美女被艹到高潮喷水动态| 亚洲欧洲日产国产| 尾随美女入室| 亚洲欧美日韩卡通动漫| 大陆偷拍与自拍| 男插女下体视频免费在线播放| 肉色欧美久久久久久久蜜桃 | 国产精品偷伦视频观看了| 热re99久久精品国产66热6| 欧美日韩视频精品一区| 成人亚洲精品av一区二区| 搡女人真爽免费视频火全软件| 亚洲图色成人| 亚洲av福利一区| 99久久精品一区二区三区| 激情五月婷婷亚洲| 内射极品少妇av片p| 久久97久久精品| 亚洲经典国产精华液单| 亚洲四区av| 中文天堂在线官网| 99热全是精品| 国产一区二区在线观看日韩| 中文字幕人妻熟人妻熟丝袜美| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 亚洲精品,欧美精品| 婷婷色av中文字幕| 国产精品秋霞免费鲁丝片| videos熟女内射| 国产精品人妻久久久影院| 成人欧美大片| 少妇裸体淫交视频免费看高清| 九草在线视频观看| 日韩免费高清中文字幕av| av在线蜜桃| 亚洲欧美一区二区三区国产| 亚洲欧美日韩无卡精品| 国产在线一区二区三区精| 免费不卡的大黄色大毛片视频在线观看| 亚洲在线观看片| 亚洲欧洲日产国产| 久久久久久久久久成人| 一区二区三区免费毛片| 国产免费福利视频在线观看| 秋霞在线观看毛片| 在线观看一区二区三区| 一级爰片在线观看| 舔av片在线| 人妻制服诱惑在线中文字幕| 亚洲va在线va天堂va国产| 国产女主播在线喷水免费视频网站| 51国产日韩欧美| 街头女战士在线观看网站| 欧美日韩亚洲高清精品| 99久久中文字幕三级久久日本| 大又大粗又爽又黄少妇毛片口| 91久久精品国产一区二区三区| 一级二级三级毛片免费看| 乱系列少妇在线播放| 日本与韩国留学比较| 国产91av在线免费观看| 免费电影在线观看免费观看| 久久久精品免费免费高清| 久久综合国产亚洲精品| 亚洲av二区三区四区| 99久国产av精品国产电影| 亚洲精品亚洲一区二区| 中文字幕亚洲精品专区| 18禁裸乳无遮挡动漫免费视频 | 一区二区三区四区激情视频| 九九在线视频观看精品| 丝袜脚勾引网站| 自拍偷自拍亚洲精品老妇| 亚洲久久久久久中文字幕| 爱豆传媒免费全集在线观看| 亚洲精华国产精华液的使用体验| 国产免费视频播放在线视频| 寂寞人妻少妇视频99o| 免费黄色在线免费观看| 99热国产这里只有精品6| 丝袜喷水一区| 国产视频首页在线观看| 国产成人精品福利久久| 国产精品熟女久久久久浪| 亚洲精品久久久久久婷婷小说| 如何舔出高潮| 国产成年人精品一区二区| 免费av观看视频| 国产免费视频播放在线视频| 麻豆久久精品国产亚洲av| 建设人人有责人人尽责人人享有的 | 国产精品不卡视频一区二区| 在线观看免费高清a一片| 欧美xxⅹ黑人| 欧美丝袜亚洲另类| 国产精品成人在线| 老师上课跳d突然被开到最大视频| 亚洲国产欧美人成| www.av在线官网国产| 男插女下体视频免费在线播放| 亚洲av中文av极速乱| 久久久久久久国产电影| 自拍欧美九色日韩亚洲蝌蚪91 | 一级毛片aaaaaa免费看小| 一级毛片久久久久久久久女| 在线天堂最新版资源| 久久久久国产精品人妻一区二区| 99热国产这里只有精品6| 夜夜爽夜夜爽视频| 精品少妇久久久久久888优播| 精品国产一区二区三区久久久樱花 | 亚洲自拍偷在线| 亚洲国产成人一精品久久久| 亚洲伊人久久精品综合| 免费少妇av软件| 日本猛色少妇xxxxx猛交久久| 可以在线观看毛片的网站| av一本久久久久| 99re6热这里在线精品视频| 中国国产av一级| 久久精品国产亚洲av涩爱| 亚洲,一卡二卡三卡| 国产毛片a区久久久久| 七月丁香在线播放| 精品久久久噜噜| 黑人高潮一二区| 欧美97在线视频| 一级a做视频免费观看| 人妻 亚洲 视频| 中国国产av一级| 99re6热这里在线精品视频| 亚洲精品久久午夜乱码| 十八禁网站网址无遮挡 | 肉色欧美久久久久久久蜜桃 | 国产有黄有色有爽视频| 欧美激情在线99| 亚洲国产日韩一区二区| 免费观看无遮挡的男女| 亚洲色图av天堂| 国产免费一区二区三区四区乱码| 在线a可以看的网站| 插阴视频在线观看视频| 99视频精品全部免费 在线| 国产亚洲精品久久久com| 国产精品一区二区性色av| 下体分泌物呈黄色| 亚洲精品自拍成人| 人妻系列 视频| 欧美日韩视频精品一区| 在线亚洲精品国产二区图片欧美 | 一二三四中文在线观看免费高清| 亚洲精品aⅴ在线观看| 中文字幕制服av| 一级a做视频免费观看| 亚洲色图av天堂| 久久精品夜色国产| 美女国产视频在线观看| 欧美老熟妇乱子伦牲交| 美女cb高潮喷水在线观看| 国产精品无大码| 国产成人精品一,二区| 精品少妇久久久久久888优播| 三级国产精品欧美在线观看| 成人亚洲精品av一区二区| 精品久久国产蜜桃| 网址你懂的国产日韩在线| 亚洲精品第二区| 高清日韩中文字幕在线| 麻豆成人av视频| 日韩亚洲欧美综合| 日韩精品有码人妻一区| 另类亚洲欧美激情| 亚洲va在线va天堂va国产| 亚洲自拍偷在线| 成人综合一区亚洲| 少妇猛男粗大的猛烈进出视频 | 国产91av在线免费观看| 久久精品综合一区二区三区| 色网站视频免费| 国产精品一二三区在线看| 97超视频在线观看视频| 成人高潮视频无遮挡免费网站| 乱系列少妇在线播放| 久久久成人免费电影| 亚洲欧美成人精品一区二区| 亚洲欧美清纯卡通| 亚洲精品成人久久久久久| 人妻 亚洲 视频| 一区二区三区乱码不卡18| 日韩欧美精品v在线| 国产色婷婷99| 一个人观看的视频www高清免费观看| 亚洲av国产av综合av卡| 国产精品麻豆人妻色哟哟久久| 亚洲色图综合在线观看| 91午夜精品亚洲一区二区三区| 亚洲天堂国产精品一区在线| 亚洲精品久久午夜乱码| 久久久亚洲精品成人影院| 97在线人人人人妻| 精品少妇久久久久久888优播| 国产爱豆传媒在线观看| 日日啪夜夜爽| 久久久久久久精品精品| 天天躁夜夜躁狠狠久久av| 亚洲天堂国产精品一区在线| 久久精品久久精品一区二区三区| 男女啪啪激烈高潮av片| 国产综合懂色| 国产黄片视频在线免费观看| 男女边摸边吃奶| 久久久a久久爽久久v久久| 亚洲,一卡二卡三卡| 日日啪夜夜撸| 欧美一级a爱片免费观看看| 91精品伊人久久大香线蕉| 国产精品偷伦视频观看了| 联通29元200g的流量卡| 国模一区二区三区四区视频| 少妇人妻一区二区三区视频| 久久久精品免费免费高清| 高清欧美精品videossex| 少妇被粗大猛烈的视频|