999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

翼型空化起始對應(yīng)空化數(shù)及尺度效應(yīng)分析

2022-07-05 03:41:16常晟銘丁恩寶孫聰王超劉貴申
中國艦船研究 2022年3期
關(guān)鍵詞:效應(yīng)模型

常晟銘,丁恩寶,孫聰,王超,劉貴申

1 中國船舶科學(xué)研究中心,江蘇 無錫 214082

2 哈爾濱工程大學(xué) 船舶工程學(xué)院,黑龍江 哈爾濱 150001

0 引 言

空化是船海領(lǐng)域中經(jīng)常發(fā)生的一種現(xiàn)象,因包含了相變、可壓縮性、非定常等特點,長久以來一直是該領(lǐng)域的研究熱點。一般而言,空化會帶來一系列不利影響,空化起始和潰滅過程會使諸如魚雷、艇體、螺旋槳等海洋結(jié)構(gòu)物發(fā)生劇烈振動,進而引發(fā)噪聲、結(jié)構(gòu)破壞、效率降低等問題,所以也是海洋結(jié)構(gòu)物需要極力避免的現(xiàn)象[1]。雖然空化對海洋結(jié)構(gòu)物的危害一直是研究者們關(guān)注的問題,但對于其成因仍無定論。目前,廣為接受的表述是由Franc 和Michel 提出的,即空化是液體介質(zhì)在過高應(yīng)力作用下的碎裂過程,高溫和低壓的條件均會使液體介質(zhì)的水分子間的拉應(yīng)力劇烈增加,使液體發(fā)生碎裂并發(fā)生空化。在一定溫度下的臨界壓力便被稱為飽和蒸汽壓[2]。此外,空化的發(fā)生還離不開一種削弱其抗拉強度的物質(zhì),這種物質(zhì)被稱為“汽核”[3]。

近年來,許多研究者針對空化,特別是與翼型空化相關(guān)的機理問題進行了研究。例如,劉志輝[4]采用OpenFOAM 開源軟件,數(shù)值模擬二維水翼的非定常云空化,分析云空化形成過程及其分級斷裂現(xiàn)象,驗證了一種非線性湍流模型及其適用性。結(jié)果表明,該非線性湍流模型在研究梢渦空泡的過程中具有一定優(yōu)勢。劇冬梅等[5]采用實驗方法對Clark-Y 翼型的空化特性進行了研究,通過改變不同攻角,觀察空化數(shù)對翼型空化形式的影響機理,結(jié)果表明,隨著空化數(shù)的不斷減小,空化可以劃分為初始空化、片空化、云空化和超空泡這幾個類型。其中,云空化對翼型的水動力特性影響程度最大,亦即顯著增加了翼型的阻力,減小了翼型的升力,從而帶來一系列不利影響。丁恩寶等[6]采用STAR-CCM+軟件,分析不同尺度下翼型在流場域中的壓力場,以及通過改變流速確定翼型空化起始點對應(yīng)的臨界雷諾數(shù),分析尺度效應(yīng),結(jié)果表明,尺度較大的翼型其臨界雷諾數(shù)也較大。

從目前來看,雖然研究者采用了許多方法研究了翼型空化問題,但對于翼型空化尺度效應(yīng)的研究卻十分片面,僅是從臨界雷諾數(shù)下產(chǎn)生尺度效應(yīng)的角度進行了研究,而未引入空化模型。因此,關(guān)于空化起始點的尺度效應(yīng)問題還需要更深入的探究。

本文將以標準KCS 船使用的NACA 0012 翼型船舵作為參照,使用SSTk-ω 湍流模型和基于輸運方程的S-S 空化模型,通過改變環(huán)境壓強來改變空化數(shù),尋找空化起始對應(yīng)的臨界空化數(shù),然后對比不同縮比尺翼型空化起始點的流場、壓力場、表面壓力分布,研究不同攻角下翼型空化起始對應(yīng)的臨界空化數(shù)隨尺度變化的尺度效應(yīng),并分析相應(yīng)的機理。

1 數(shù)值計算方法

1.1 湍流模型

本文使用的模型為SSTk-ω 湍流模型[7-8],其中k,ω 的輸運方程分別為:

式中:Gk為由層流速度梯度產(chǎn)生的湍流動能;Gω為由ω 方程產(chǎn)生的湍流動能;k和ω 分別為湍流動能及耗散率; ρ為流體密度;ui為速度分量的時均值;xi和xj為位移分量;Tk和Tω分別為k和ω 的擴散率;Yk和Yω為由擴散產(chǎn)生的湍流;Dω為正交發(fā)散項。

1.2 空化模型

在研究空化問題時,普遍采用的是均質(zhì)混合流模型,主要思想是將空化發(fā)生時的液體和蒸汽兩相視為由球形汽泡及液體兩種成分組成的混合流,能夠?qū)⒍嘞噙B續(xù)介質(zhì)簡化成一個包含汽核的單相連續(xù)介質(zhì),進而通過控制方程求解。而Schnerr-Sauer 空化模型(以下稱S-S 空化模型)是一種基于均質(zhì)混合流的輸運方程模型[9],其中汽相體積分數(shù)的定義為

構(gòu)建S-S 空化模型源項的表達式為

2 計算模型

2.1 流場域劃分和工況設(shè)置

圖1 所示為劃分的計算域。其中,左側(cè)邊界為速度進口,右側(cè)邊界為壓力出口,上、下邊界為對稱平面。為滿足黏性流體力學(xué)條件,翼型表面的工況被設(shè)置為不可滑移壁面。計算域左側(cè)速度進口的邊界距離翼型左緣2 倍弦長(C),右側(cè)壓力出口的邊界距離翼型的右緣6 倍弦長,上、下兩個邊界距離翼型上緣和下緣均為2 倍弦長。流體密度和動力黏度系數(shù)為25 ℃的液態(tài)水參數(shù),設(shè)置流體密度為997.561 kg/m3,動力黏度系數(shù)為8.887 1×10?4Pa·s。

圖1 計算域的劃分Fig. 1 Division of computational domain

本文參照標準KCS 船所使用的NACA 0012型船舵,翼型的初始尺度與船舵的弦長一致,為5.5 m,而流場的流速對應(yīng)KCS 船的航速24 kn,即12.345 m/s。分別采用1:3,1:10 和1:20 縮尺比翼型進行尺度效應(yīng)的研究。在研究尺度效應(yīng)時,因參照了實船舵所選取的翼型,一般需要保證弗勞德數(shù)相等,而若要保證弗勞德數(shù)相等,不同尺度之間的翼型所對應(yīng)的流場流速會有一定差異,導(dǎo)致雷諾數(shù)的不同,并帶來尺度效應(yīng)。表1 給出了不同尺度翼型對應(yīng)的弦長及流場流速。

根據(jù)前人的研究經(jīng)驗,平板湍流的邊界層厚度需與其對應(yīng)的雷諾數(shù)及參照長度相關(guān),具體的數(shù)量關(guān)系為[10]

式中:h為平板湍流邊界層厚度;L為參照長度,本文即翼型的弦長;Re為對應(yīng)的雷諾數(shù)。

因縮尺比的變化會導(dǎo)致表1 所示的翼型雷諾數(shù)發(fā)生變化,根據(jù)式(9),所以對應(yīng)的邊界層厚度也會產(chǎn)生較大差異,進而對網(wǎng)格劃分帶來一定影響,即影響到近壁第1 層棱柱層厚度和棱柱層厚度。經(jīng)公式推導(dǎo),不同縮尺比翼型對應(yīng)的棱柱層厚度及近壁面第1 層棱柱層厚度如表2 所示。

表1 不同尺度翼型的工況設(shè)置Table 1 Working condition setting of different scale models of hydrofoil

表2 不同尺度翼型的邊界層參數(shù)設(shè)置Table 2 Boundary layer parameter setting of different scale models of hydrofoils

本文需要對不同縮尺比和攻角的翼型空化起始所對應(yīng)的空化數(shù)進行研究。對應(yīng)的空化數(shù)定義如下:

式中: σ為空化數(shù);Pout為出口壓力,即環(huán)境壓強;ρw為液體密度;v為進口速度。

2.2 網(wǎng)格劃分

本文采用的網(wǎng)格為切割體網(wǎng)格和棱柱層網(wǎng)格相結(jié)合的形式,其中,翼型表面的邊界層設(shè)置為棱柱層網(wǎng)格形式,流場域網(wǎng)格設(shè)置為切割體網(wǎng)格形式。為提高計算精度,將翼型上表面前緣位置的網(wǎng)格和尾流區(qū)進行了一定程度的加密。最終,確定了如圖2 所示的網(wǎng)格劃分形式。當(dāng)改變翼型縮尺比時,需通過改變對應(yīng)網(wǎng)格基礎(chǔ)尺寸的方式,將整個計算域作放大或縮小處理。

圖2 網(wǎng)格劃分Fig. 2 Mesh division of computational domain

為驗證網(wǎng)格的收斂性,本文按照圖2 所示的網(wǎng)格劃分形式,通過改變網(wǎng)格的基礎(chǔ)尺寸,選擇了3 套不同加密形式的網(wǎng)格進行收斂性驗證。模型計算工況如下:弦長為5.5 m、原始尺度翼型攻角為5°、空化數(shù)為1.15。具體網(wǎng)格數(shù)及3 套網(wǎng)格劃分模型所對應(yīng)的升力系數(shù)CL和阻力系數(shù)CD如表3 所示。

表3 網(wǎng)格收斂性分析Table 3 Mesh convergence analysis

由表3 可見,相較于較為粗糙的Mesh-1 網(wǎng)格加密形式,Mesh-2 和Mesh-3 這兩種加密形式的網(wǎng)格模型二者的升力系數(shù)和阻力系數(shù)的值較接近,因此網(wǎng)格加密到Mesh-2 的程度已可以滿足對計算所需精度的要求,無需進一步細化網(wǎng)格。此外,本文還需對空化起始的機理作一定的研究,網(wǎng)格收斂性驗證中需對比翼型表面壓力分布。具體結(jié)果如圖3 所示。

由圖3 可見,3 種不同網(wǎng)格加密形式的翼型表面壓力分布結(jié)果之間差異并不顯著。因此,基于上述網(wǎng)格收斂性分析結(jié)果,本文最終確定了Mesh-2 網(wǎng)格加密形式。

圖3 不同網(wǎng)格加密形式的翼型表面壓力分布對比Fig. 3 Comparison of surface pressure distribution of hydrofoil with different forms of mesh refinement

2.3 模擬值與試驗值對比

為證明本文所用方法的合理性,將Mesh-2 網(wǎng)格加密形式應(yīng)用于攻角6°、弦長0.15 m 的NACA 66翼型,以及在雷諾數(shù)Re=8×105工況下通過數(shù)值模擬方式得到的升阻比(CL/CD)與試驗值[11]進行對比,如圖4(a)所示。圖4(b)給出的是空化數(shù)為1.41 時的空泡分布結(jié)果。

圖4 模擬值與試驗值的對比Fig. 4 Comparison between simulations and experimental results

由圖4 可見,在現(xiàn)有Mesh-2 網(wǎng)格劃分的基礎(chǔ)上,數(shù)值模擬結(jié)果與試驗值較接近,也就是說,數(shù)值模擬能夠精確模擬出翼型周圍的流場,并可解決與翼型空化起始相關(guān)的水動力機理問題。

3 計算結(jié)果

3.1 空化起始對應(yīng)空化數(shù)的確定方法

為確定空化起始對應(yīng)的臨界空化數(shù),應(yīng)明確空化形式隨著空泡數(shù)的不斷減小,會逐漸發(fā)展為起始空泡、片空泡、云空泡、超空泡這幾種形式。根據(jù)式(10)的定義,當(dāng)環(huán)境壓強足夠大時,空泡數(shù)會很大,空化并不會發(fā)生。因此,在保證翼型周圍流速一定的情況下,設(shè)置足夠大的環(huán)境壓強,然后再不斷減小環(huán)境壓強,直至空化數(shù)縮小到起始空化發(fā)生之時,即可確定對應(yīng)工況下空化起始對應(yīng)的臨界空化數(shù)。圖5 所示為原始尺度翼型(弦長5.5 m 的NACA 0012 翼型)0°攻角時的空化起始對應(yīng)空化數(shù)的確定方式。

圖5 翼型空化起始對應(yīng)空化數(shù)的確定方式Fig. 5 Method for ascertaining the cavitation number corresponding to cavitation initiation of original scale hyrofoil

根據(jù)前人的經(jīng)驗[4],將汽相體積分數(shù)分布的最大值設(shè)置為0.1 時能夠清晰地捕捉到空化邊界,故本文也遵循此設(shè)置方法。由圖5(a)所示空化數(shù)為0.45 時的汽相分布結(jié)果不難看出,此時因環(huán)境壓強較大,液體難以發(fā)生汽化,空化并不會發(fā)生。而當(dāng)空化數(shù)減小至0.40 時(圖5(b)),因環(huán)境壓強過低,空化已發(fā)展為片空化形式,并非空化起始形式。最終,通過在0.40~0.45 的空化數(shù)范圍內(nèi)改變環(huán)境壓強,可以得到如圖5(c)所示的空化起始形式,最終確定空化起始對應(yīng)的臨界空化數(shù)為0.43。

3.2 空化起始對應(yīng)空化數(shù)隨攻角的變化

在流場流速保持一定時,若改變翼型的攻角,周圍的流場特性也會隨之發(fā)生一定的改變,并影響翼型空化起始對應(yīng)的空化數(shù)。本文以原始尺度翼型為例,給出了攻角與空化起始對應(yīng)空化數(shù)之間的關(guān)系,如圖6 所示。

圖6 原始尺度翼型空化起始對應(yīng)空化數(shù)與攻角的關(guān)系Fig. 6 Relationship between cavitation number and attack angle corresponding to cavitation initiation of original scale hydrofoil

由圖6 可見,原始尺度翼型的空化起始對應(yīng)的臨界空化數(shù)隨攻角的增加而不斷增大,二者的關(guān)系呈現(xiàn)為斗狀曲線,一般被稱為“空泡斗”[12]。在空泡斗的內(nèi)部,因空化數(shù)大于空化起始對應(yīng)的臨界空化數(shù),空化難以發(fā)生,所以此區(qū)域也被稱作“無空泡區(qū)”;而空泡斗的外部因空化數(shù)會小于空化起始對應(yīng)的臨界空化數(shù),在翼型上表面強吸力峰的位置形成空化,所以空泡斗的外部區(qū)域也被稱作“空泡區(qū)”。

圖7 所示為原始尺度翼型在不同攻角時對應(yīng)的空化起始點空泡形態(tài)分布圖。

由圖7 可見,隨著攻角的不斷增加,翼型空化起始位置逐漸向翼型前緣方向移動。在攻角為0°時,翼型上、下兩個表面均有空化起始,但在有攻角之后,空化起始的位置接近于翼型上表面前緣位置。從空化起始的發(fā)展程度來看,翼型攻角較小,即0°和5°時,空化起始的發(fā)展程度較小,但隨著翼型攻角的不斷增大,例如10°時,翼型附近流場發(fā)生了較大變化,從而影響到翼型上表面低壓區(qū)域的面積,因此空化起始的發(fā)展程度也會有所增強。為深入分析流場對空化起始點的作用機理,提取了5°和15°攻角下翼型周圍的流線圖(圖8)、速度矢量圖(圖9)和壓力云圖(圖10)。

圖7 原始尺度翼型空泡起始點空泡形態(tài)分布隨攻角的變化Fig. 7 Variation of the cavitation shape distribution of the cavitation initial point of the original scale hydrofoil with attack angle

圖8 不同攻角翼型周圍的流線分布Fig. 8 Streamline distribution of hydrofoil at different attack angles

由圖8 可見,翼型周圍的繞流場整體上比較平穩(wěn),但其前緣及尾緣附近的流線有一定程度的扭曲。相較于5°攻角時的情況,15°攻角時的翼型前緣流線扭曲程度更高,分布密度也更大,使得前緣附近流場的流速更大,流速較大的區(qū)域也明顯增加,如圖9 所示。

圖9 不同攻角翼型的速度矢量分布Fig. 9 Velocity vector distribution of hydrofoil at different attack angles

根據(jù)伯努利原理(Bernoulli's theorem),在同一流線上流速越大的位置壓力越小。因此,一方面,翼型前緣附近流場的流速較大會導(dǎo)致壓強有所降低,更容易低于水的飽和蒸汽壓,即空化更容易起始于翼型前緣位置的附近;另一方面,在15°攻角時,翼型前緣附近較大流速所在的區(qū)域較大(圖9),使得翼型前緣附近的低壓區(qū)域面積較大(圖10),從而極大地增加了此時翼型空化起始的發(fā)展程度。

圖10 不同攻角翼型的壓力云圖Fig. 10 Pressure contours of hydrofoil at different attack angles

根據(jù)科恩達效應(yīng)(Coanda effect),翼型前緣附近的流體會貼近翼型輪廓流動,受到一個向心力的作用,其大小朝向翼型方向,來源于流場中流體壓力的合力。這證明了翼型前緣附近的流場壓力較小,容易發(fā)生空化的這一結(jié)論。

3.3 空化起始對應(yīng)空化數(shù)的尺度效應(yīng)

隨著翼型尺度的變化,翼型雷諾數(shù)也會有一定的變化,進而發(fā)生尺度效應(yīng),使得空化起始對應(yīng)空化數(shù)有所變化。為探究此時翼型的尺度效應(yīng),本文選擇了1:3,1:10,1:20 縮尺比翼型進行分析,得到了不同攻角翼型時空化起始對應(yīng)的空化數(shù),如圖11 所示。

由圖11 可見,1:3 縮尺比翼型并未受到尺度效應(yīng)的影響,說明無論攻角如何變化,翼型空化起始對應(yīng)的空化數(shù)均與原始尺度翼型的相同。但是,隨著翼型縮尺比的變小,翼型空化起始對應(yīng)的空化數(shù)也均有所降低,且翼型縮尺比越小,空化數(shù)降幅越大。相較于其他攻角的情況,15°攻角翼型空化起始對應(yīng)的空化數(shù)與原始尺度翼型的差別非常大,說明此時的尺度效應(yīng)最明顯。圖12 給出了15°攻角翼型空化起始對應(yīng)的空化數(shù)隨縮尺比的變化趨勢。

圖11 不同尺度翼型空化起始對應(yīng)空泡數(shù)與攻角的關(guān)系Fig. 11 Relationship between cavitation number and attack angle corresponding to cavitation initiation of different scale hydrofoils

圖12 15°攻角翼型空化起始所對應(yīng)空泡數(shù)的尺度效應(yīng)Fig. 12 Scale effect of cavitation initiation corresponding to cavitation number at fifteen degrees attack angle

對于翼型空化起始所對應(yīng)的空化數(shù),改變翼型縮尺比,其受到的尺度效應(yīng)影響隨之變化,空化起始位置也隨翼型縮尺比的變化而變化。圖13給出了5°,10°,15°不同攻角翼型的空化起始位置受尺度效應(yīng)的影響情況。

由圖13 可見,各攻角下,隨著翼型縮尺比發(fā)生變化,翼型空化起始位置也都因尺度效應(yīng)而改變。具體而言,在翼型縮尺比縮小的情況下,翼型空化起始位置逐漸向其后緣方向偏移,而空化起始的發(fā)展程度也隨之逐漸增強。因此,在研究水翼時,特別是模型試驗研究中,翼型尺度效應(yīng)的影響不容忽視,采用的水翼模型尺度也不可與實體尺度相差過大。

圖13 不同攻角翼型空化起始位置隨尺度的變化Fig. 13 Variation of cavitation initiation position of hydrofoils with scale at different attack angles

為探究尺度效應(yīng)對翼型空化起始對應(yīng)的空化數(shù)及其位置的影響機理,本文提取了15°攻角下原始尺度翼型和1:20 縮尺比翼型的流線圖、速度矢量圖和壓力云圖進行對比,分別如圖14、圖15和圖16 所示。

由圖14 可見,對于原始尺度和1:20 縮尺比翼型,在15°攻角下二者的空化起始位置、流線密集程度和形狀扭曲程度幾乎相同,未受到尺度效應(yīng)的影響。但是,在尺度效應(yīng)的影響下,1:20 縮尺比翼型上表面的第1 條流線與翼型之間的間距明顯大于原始尺度翼型。相同的結(jié)果同樣在圖15 所示的速度矢量分布中得到反映,即原始尺度翼型周圍流體流動與翼型上表面的貼合程度大于1:20縮尺比翼型的情況。

造成上述現(xiàn)象的原因是,采用式(9)計算邊界層厚度后,1:20 縮尺比翼型邊界層的相對厚度要大于原始尺度翼型的邊界層厚度。二者邊界層厚度之間的差異,一方面使得1:20 縮尺比翼型附近流場受到黏性的影響較大,特別是貼近翼型表面的流體分子受到較大的黏性剪應(yīng)力作用,阻礙了流體流動;另一方面,會影響到控制方程的求解及邊界層的轉(zhuǎn)捩過程,以及不同縮尺比翼型周圍流場的分布,通過伯努利原理使翼型周圍壓力分布產(chǎn)生如圖16 所示的差異,帶來了一系列的尺度效應(yīng),最終影響了空化起始對應(yīng)的空化數(shù)及其位置和形態(tài)。

圖16 不同尺度翼型15°攻角下的壓力云圖Fig. 16 Pressure contours of different scale hydrofoils at fifteen degrees attack angle

4 結(jié) 論

本文采用了STAR-CCM+軟件,采用SSTk-ω湍流模型和S-S 空化模型數(shù)值模擬了翼型空泡特性,分析了不同縮尺比翼型的空化起始對應(yīng)空化數(shù)受尺度效應(yīng)的影響,得到了如下的結(jié)論:

1) 對于相同縮尺比翼型,隨著攻角的增大,空化起始對應(yīng)空化數(shù)有所增大,二者的關(guān)系呈現(xiàn)為一個斗狀曲線,也稱“空泡斗”。

2) 隨著翼型攻角的增大,翼型空化起始位置不斷向翼型前緣的方向移動,空化起始的發(fā)展程度隨之增強。其原因在于,隨著翼型攻角的增大,翼型前緣流線的扭曲程度有所增強,分布密集程度也更大,前緣附近流場流速也有所增大。根據(jù)伯努利原理和科恩達效應(yīng),流場的差異作用于壓力場,使得翼型空化起始受到影響。

3) 隨著翼型尺度的縮小,不同縮尺比翼型空化起始對應(yīng)的空化數(shù)均有所降低,尺度越小,翼型空化起始對應(yīng)空化數(shù)降低的幅度越大,即受尺度效應(yīng)的影響越大。而翼型空化起始位置隨尺度的減小向翼型后緣方向偏移,發(fā)展程度隨尺度的減小而有所增強。造成尺度效應(yīng)的原因是,翼型周圍邊界層的相對厚度因不同雷諾數(shù)帶來的差異而有所變化,影響了翼型表面流場受到的黏性剪應(yīng)力作用,以及控制方程的求解和邊界層轉(zhuǎn)捩的過程。

針對翼型空化的尺度效應(yīng)問題,未來還需要進一步開展如下研究:

1) 本文所述翼型空化研究僅限于二維層面,而空化往往是一個三維及非定常的問題,特別是對于三維尺度的翼型空化問題,在三維翼型端部還會形成一個內(nèi)卷的流向渦,發(fā)生梢渦空化。上述三維特性通過二維方法是無法研究的,未來應(yīng)在三維基礎(chǔ)上研究翼型空化問題。

2) 對于翼型空化的尺度效應(yīng)研究,空化起始的尺度效應(yīng)僅是其中一小部分??栈凑掌浒l(fā)展階段包括了起始空化、片空化、云空化和超空化,這些階段發(fā)生點對應(yīng)的空化數(shù)及空化形態(tài)均因雷諾數(shù)的差異會帶來一系列尺度效應(yīng),這也是未來需要著重研究的方面。

猜你喜歡
效應(yīng)模型
一半模型
鈾對大型溞的急性毒性效應(yīng)
懶馬效應(yīng)
場景效應(yīng)
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
應(yīng)變效應(yīng)及其應(yīng)用
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
偶像效應(yīng)
主站蜘蛛池模板: 久久综合五月婷婷| 人妻一本久道久久综合久久鬼色| 日韩毛片免费视频| 午夜不卡视频| a毛片在线免费观看| 香蕉eeww99国产在线观看| 国产福利微拍精品一区二区| 色妞www精品视频一级下载| 国产日韩欧美一区二区三区在线| 毛片在线播放网址| 国产福利一区视频| 伊人狠狠丁香婷婷综合色| 免费又爽又刺激高潮网址| 免费黄色国产视频| 久久亚洲国产最新网站| 国产精品久久久精品三级| 国产精品免费久久久久影院无码| 在线亚洲精品自拍| 国产无人区一区二区三区| 亚洲第一成年免费网站| 国产不卡在线看| 日韩色图在线观看| 国产精品久久自在自线观看| 一级爱做片免费观看久久| 欧美精品三级在线| 91国内视频在线观看| 国产欧美日韩91| 国产精品一区二区无码免费看片| av免费在线观看美女叉开腿| 国产精品一区二区在线播放| 亚洲女人在线| 蜜桃臀无码内射一区二区三区| 中文无码日韩精品| 国产高清在线精品一区二区三区| 久久国产精品娇妻素人| 青青草原偷拍视频| 亚洲精品第一页不卡| 国产视频一二三区| 中文字幕 欧美日韩| 亚洲国产综合自在线另类| 日韩a级毛片| 性欧美在线| 亚洲无限乱码一二三四区| 四虎AV麻豆| 国产精品美女自慰喷水| 日韩一区精品视频一区二区| 丁香婷婷激情网| 台湾AV国片精品女同性| 亚洲美女久久| 找国产毛片看| 色噜噜中文网| 欧美色图久久| 日韩经典精品无码一区二区| 国产人人干| 久久人搡人人玩人妻精品| 久久精品人妻中文系列| 国产噜噜噜| 青青青国产视频| 欧美日韩在线成人| 中文字幕天无码久久精品视频免费 | 欧美国产日产一区二区| 国产成人一区二区| 中文字幕佐山爱一区二区免费| 中国国产A一级毛片| 2020精品极品国产色在线观看| 久久久久亚洲精品成人网| 精品视频91| 亚洲第一黄色网址| 欧美不卡视频在线观看| 人妻精品久久无码区| 日本高清在线看免费观看| 国产麻豆va精品视频| 国产亚洲精| 亚洲精品国产综合99| 午夜福利视频一区| 潮喷在线无码白浆| 伊人网址在线| 999国内精品视频免费| 丁香五月婷婷激情基地| 亚洲综合欧美在线一区在线播放| 成人精品区| 国产菊爆视频在线观看|