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

基于URANS與DDES方法的空腔近場噪聲數值研究

2016-11-24 06:36:47楊黨國王顯圣羅新福
振動與沖擊 2016年20期
關鍵詞:結構方法

劉 俊, 楊黨國, 王顯圣, 羅新福

(1.中國空氣動力研究與發展中心 空氣動力學國家重點實驗室,四川 綿陽 621000;2.中國空氣動力研究與發展中心 高速空氣動力研究所,四川 綿陽 621000)

?

基于URANS與DDES方法的空腔近場噪聲數值研究

劉 俊1, 楊黨國2, 王顯圣2, 羅新福2

(1.中國空氣動力研究與發展中心 空氣動力學國家重點實驗室,四川 綿陽 621000;2.中國空氣動力研究與發展中心 高速空氣動力研究所,四川 綿陽 621000)

采用基于SST(Shear-Stress Transport)湍流模式的URANS(Unsteady Reynolds-Averaged Navier-Stokes)和DDES (Delayed Detached Eddy Simulation)方法開展了馬赫數0.85的三維空腔非定常流動數值計算。計算結果表明:兩種方法得到的空腔底部靜壓、脈動壓力聲壓級和功率譜均與實驗及參考文獻結果具有良好的一致性;DDES在模擬流動失穩、小尺度結構等流動細節方面更具優勢,對高頻壓力脈動的捕捉也要優于URANS。通過對時均流場的分析,確定了模擬的空腔流動類型為過渡式流動,同時發現空腔內存在的復雜三維渦結構,并分析了這些渦結構對空腔流場特性的影響。

URANS;DDES;空腔噪聲;脈動壓力;渦結構

空腔流動是一種常見的流動現象。大自然中的山谷、城市中樓宇的間隙、行駛中開窗的汽車、大型客機的起落架艙、戰斗機的內埋彈艙、航天飛機隔熱板的縫隙等都可以簡化為空腔外形,當流體經過這些大大小小的空腔時,將演化出形式多樣的空腔流動現象。空腔流動可以看成后臺階流動和前臺階流動的組合,但是當前后臺階距離較近時,這兩種流動就會產生強烈的相互干擾,誘發更加復雜的流動現象。

研究發現影響空腔流動的因素包括空腔的形狀(方腔、圓柱腔等)、尺寸、長深寬比例等幾何因素以及馬赫數、雷諾數、迎角、邊界層厚度、邊界層形態等流動因素。在高馬赫數、高雷諾數的條件下,空腔中存在激波、剪切層、旋渦、壓縮/膨脹波、聲波等復雜流動結構及這些結構間的相互干擾。與這種非定常、非線性流動相伴而生的是劇烈的壓力脈動,試驗表明壓力脈動對應的聲壓級最高可達170 dB,極易誘發結構振動,導致疲勞失效。

為了更好地解決空腔流致噪聲等問題,近半個世紀以來,人們采用理論分析、風洞試驗和數值計算等手段,開展了大量的研究。在風洞試驗方面,研究人員主要采用表面脈動壓力傳感器[1-5]和PIV[6-7](Particle Image Velocimetry)等進行空腔噪聲及流動特性測量。其中,脈動壓力傳感器主要布置在物體表面,難以獲得空間的壓力變化信號,雖然其時間采樣頻率可達幾萬赫茲,但是空間分辨率不足;而PIV等非接觸測量手段能夠較清楚解析空間流場,具有較高的空間分辨率,但是由于采樣頻率低,時間解析精度差。為了精細模擬空腔內的多尺度湍流結構,需要同時具備較高的空間和時間分辨率,而數值計算具備這方面的優勢,因此在空腔流動研究方面數值計算手段受到越來越多的重視。

20世紀80年代以來,人們采用數值手段研究了多種空腔構型在不同來流條件下的流動,并嘗試了大量的湍流模型,從二維和三維的非定常雷諾平均模擬[8-9,23](URANS)到大渦模擬[10](Large Eddy Simulation,LES),再到最近的直接數值模擬[11](Direct Numerical Simulation.DNS)。大量的研究結果表明,提升空腔流場和脈動壓力預測精度的關鍵仍在于對湍流的模擬。湍流直接影響剪切層的不穩定性及其在向下游傳播過程中的混合過程,以及剪切層對空腔后壁的撞擊,從而對空腔主要流動現象—流聲耦合作用產生影響,因此湍流模擬在空腔流動數值計算中起著決定性作用。傳統的RANS方法由于在時間平均過程中損失了湍流能譜,因此在模擬高頻壓力脈動方面存在嚴重缺陷。另一方面,LES和DNS所需的計算資源又超出了目前的承受能力,尤其是在高馬赫數和高雷諾數的情況下。近年來,在湍流模擬領域發展起來的脫體渦模擬[10,12-19](Detached Eddy Simulation,DES)方法,是一種RANS和LES的混合方法,這種方法吸收了RANS方法模擬近壁湍流對網格分辨率要求低以及LES方法在模擬分離區流動耗散小的優點,較好地模擬了大分離流動現象。對于空腔流動而言,剪切層和非定常渦等流動特征適合用DES中的LES模式解析,而近壁湍流和來流邊界層則可采用RANS方法模擬,因此采用DES方法能較好的解析空腔流動中的大部分流動特征。

為進一步探索適用于空腔流動的湍流模擬方法,本文采用URANS方法和延遲的脫體渦模擬(DDES)方法研究了馬赫數為0.85的三維空腔高速復雜流動特性,并對比分析了兩種方法在模擬空腔非定常流動、時均流動、壓力脈動等方面的表現。

1 數值方法

對于URANS方程,采用基于格點的有限體積法進行離散,對流項通過二階迎風格式重構得到,黏性項采用二階中心離散格式。湍流模型為Menter的k-w-SST兩方程模型。時間推進采用雙時間步的二階歐拉隱式格式,并采用多重網格方法加速收斂。

基于SST湍流模式的DES方法通過對湍動能耗散項的改造,實現RANS模式和LES模式的切換,改造后耗散項Dk為:

Dk=β*ρkω·FDES

FDES=max[(1-FSST)lRANS/lLES,1]

(1)

2 空腔模型及計算狀態

空腔長度L為0.508 m,寬度W和深度D均為0.102 m,長寬深比例為L∶W∶D=5∶1∶1。入口距離空腔前緣0.788 m,出口距離空腔后緣0.533 m,空腔模型見圖1。物面采用絕熱無滑移條件。計算網格如圖2所示,空腔分離區域需要滿足DES方法中 LES 模式的網格要求,因此在這一區域幾乎保持了網格各向同性,流向、法向和展向網格尺度相差在1.05倍以內。在空腔四周采用了減網格技術,使有限的網格更加集中于空腔內部及其周圍。網格總量為360萬,其中空腔內部及其附近網格為109萬,占比30%。時間步長為0.005倍無量綱時間,約為9.6 μs,共推進23 000 步。來流馬赫數為0.85,來流靜壓82 100 Pa,來流靜溫266.53 K。為了與試驗結果對比,在空腔底部沿流向設置了K20K29共10個壓力監測點。

圖1 空腔模型Fig.1 Cavity model

圖2 空腔流動計算網格Fig.2 Computational grids for the cavity flow

3 計算結果分析

3.1 非定常流動特征

采用Q準則(速度梯度的第二不變量)提取了空腔的瞬時三維渦結構(見圖3)。由于Kelvin-Helmholtz不穩定性,空腔上方剪切層易發生失穩形成旋渦。旋渦在向下游對流過程中不斷發展,尺寸和強度都得到加強,與后壁面相碰后,一部分旋渦回流進入空腔內部,一部分旋渦流入尾跡區。對比兩種方法的計算結果,DDES中,空腔內部及尾跡不斷有旋渦涌入,多尺度渦結構十分豐富,與實際流動相符。而URANS中僅存在一些大尺度渦。這主要是由于URANS方法在分離區耗散過大造成的。圖3中三維渦結構采用渦黏系數與層流黏性系數的比值染色,空腔內URANS的渦黏系數要比DDES大一個數量級左右。

另外,觀察空腔中間截面的密度梯度云圖(見圖4),DDES結果中可清晰地看到空腔上方剪切層失穩形成渦結構以及渦結構向下游流動并與空腔后壁相碰等典型空腔流動現象。而URANS中由于耗散過大,流動不穩定性明顯被抑制,剪切層從空腔前緣發展至后緣仍未失穩。通過分析兩種方法在空腔非定常流動特征方面的表現,可以發現,DDES與URANS相比,在精細模擬流動失穩、小尺度渦結構等流動細節方面更具優勢。

圖3 瞬時空腔三維多尺度渦結構Fig.3 Instantaneous cavity three-dimensional multi-scale vortex structures

圖4 瞬時空腔中間截面密度梯度Fig.4 Instantaneous contour of density gradient on cavity middle plane

3.2 時均流動特征

對于空腔這種強非定常流動,隨著時間的推進,流場的變化十分劇烈,從瞬時結果中難以分析出流動的一般規律。為了定量地分析流場,對流場進行時間平均處理,從變化的流場中提取出一些不變的流動特征。

對于中截面二維時均流場(如圖5所示),URANS和DDES方法得到類似的流場結構,空腔內部均存在兩個未完全分離的渦結構,前緣附近的渦占據的面積較大,后緣附近的渦較小,其中URANS中的后緣渦更小。對于開式空腔流動,空腔中主要存在一個大的回旋渦;而閉式空腔流動在前后緣各有一個渦,同時流線在空腔中部附壁。因此,從渦結構形態上講,計算得到的空腔流動結構介于開式和閉式流動類型之間,屬于過渡式流動。另外,背景采用時均渦黏性系數染色,同樣可以看到URANS的渦黏系數要明顯大于DDES。

圖5 空腔中截面流線Fig.5 Streamlines on cavity middle plane

圖6展示了空腔底部中心線上壓力分布。由于沒有靜壓實驗數據,僅與參考文獻[20]中部分結果進行了對比。從對比的結果來看,本文計算的結果與其他CFD (Computational Fluid Dynamics)結果量值基本一致,變化趨勢也十分相似,尤其是本文的DDES計算結果跟FOI的HYBO以及Alenia 的EARSM-DES結果符合良好。除了靠近后壁部分(x/L從0.95~1區間),空腔底部壓力分布跟亞跨聲速過渡式流動類型的壓力曲線(見圖7)保持一致,在x/L<0.4區間,壓力系數基本維持在0附近,在0.4

圖6 空腔底部壓力分布對比[20]Fig.6 Comparison of pressure coefficient on cavity floor[20]

圖7 基于流向壓力分布的空腔流動分類[5]Fig.7 Classification for cavity flows based on streamwise pressure distribution[5]

采用Q準則提取時均流場三維旋渦結構(見圖8)。從URANS和DDES的結果來看流動基本保持對稱,但是呈現明顯的三維效應。空腔內外形成了多種旋渦結構,包括邊緣處的流向邊緣渦、尾流區的尾跡渦以及腔內的馬蹄渦和流向渦。計算得到的馬蹄渦(見圖8和圖10(a))靠近空腔的后壁,后壁區是空腔流動比較重要的區域,剪切層與固壁碰撞、腔內流體溢出等對空腔流動起重要作用的現象均發生在此區域。多篇文獻[6,21-22]均發現了在后壁區的馬蹄渦結構(如圖9所示)。在自由來流中加入障礙物,一般在障礙物周圍會產生類似馬蹄渦結構。對于空腔流動,由于流動向下游運動過程中,流線向下偏轉,后壁阻擋了這部分流動從而導致馬蹄渦的形成。分析發現,除邊緣渦卷走部分腔內流體外,大部分流入空腔的流體均通過馬蹄渦溢出空腔,因此該馬蹄渦對于及時排除進入腔內的流動,降低腔內壓力脈動具有十分重要的作用。

此外,DDES計算結果還捕捉到了一對轉向相反的沿流向發展的旋渦(見圖10(b))。PENG[12]在其計算的三維空腔流動中也發現了類似的流向旋渦(見圖11)。 觀察空腔底面壓力系數海拔面(見圖12),可以看到海拔面沿橫向基本保持一致,但是在上升斜坡上出現了兩個凹槽,凹槽剛好位于流向渦正下方。我們知道,在旋渦的卷吸作用下旋渦下方一般會形成較明顯的低壓區,如三角翼大攻角條件下,前緣渦的卷吸作用會使背風區產生大面積的低壓區,從而形成非線性渦升力。因此,空腔底部靜壓出現的這種局部低壓區可能與流向旋渦有著十分密切的關系。

圖8 空腔時均流場三維流動結構Fig.8 Time-averaged cavity three-dimensional flow structures

圖9 后緣馬蹄渦結構Fig.9 Horse-shoe vortex near the cavity trailing edge

圖10 空腔內部馬蹄渦和流向渦(DDES)Fig.10 Horse-shoe vortex and streamwise vortex in cavity interior (DDES)

圖11 x/L=0.7截面速度向量[12]Fig.11 Velocity field at x/L=0.7 plane[12]

圖12 空腔底面壓力系數海拔面(DDES)Fig.12 Altitude of pressure coefficient for the cavity floor

3.3 非定常脈動壓力特征

選取位于空腔前部和后部的K21、K28兩個典型測點(見圖1)與試驗結果對比,可以看到,本文計算的脈動壓力聲壓級的幅值以及峰值頻率與試驗均符合較好。其中,URANS和DDES在低頻部分(見圖13(a)和圖13(b))的脈動聲壓頻譜與試驗結果基本重合,然而在高頻區(見圖13(c))就會發現URANS衰減過快,而DDES符合湍流衰減的-5/3次律。

對于空腔底部脈動壓力聲壓級分布(見圖14),兩種方法均得到類似的規律:沿流動方向,聲壓級曲線先降后升,最低點位于x/L=0.15附近,最高點位于空腔后壁區域。通過與試驗及其他CFD結果比較,可以看到本研究計算的結果與試驗及參考文獻規律保持一致,但是采用CFD(包括本文及參考文獻[20])得到的聲壓級均略高于試驗結果3~4 dB。LAWSON[24]等人也發現了這一問題,他們認為,數值計算和試驗結果的差異可能是由于數值方法未能夠精確模擬空腔入口湍流邊界層及其相干結構造成的,而事實上,空腔入口湍流邊界層是影響空腔噪聲的重要因素。

圖13 脈動壓力功率譜對比Fig.13 Comparison of pressure fluctuations power spectral

圖14 空腔底部聲壓級分布對比[20]Fig.14 Comparison of distributions of OASPL at the bottom of cavity[20]

4 結 論

(1) 采用基于SST模式的URANS和DDES方法研究了馬赫數0.85的三維空腔非定常流動,數值計算結果與與試驗結果吻合較好,驗證了數值方法的可行性和正確性。

(2) URANS方法在分離區耗散大,DDES通過在分離區切換到LES模式,降低耗散,能更加精細模擬空腔內的流動失穩、小尺度結構和壓力脈動等流動細節。

(3) 兩種方法在計算空腔脈動壓力特性方面,均能分辨出不同離散頻率下的噪聲,并獲得與試驗結果較接近的聲壓級分布規律,尤其是DDES方法在對高頻部分的噪聲捕捉的更為準確。

(4) 通過分析空腔時均流場特性,發現空腔流動中存在邊緣渦、尾跡渦以及位于腔內的馬蹄渦及流向渦等三維流動結構,并推斷空腔底部出現的局部低壓區可能與流向渦存在一定的內在關聯。

[1] PLENTOVICH E B, STALLINGS R L, TRACY M B. Experimental cavity pressure measurements at subsonic and transonic speeds: static-pressure results[R]. NASA TP-3358, 1993.

[2] STALLINGS R L, PLENTOVICH E B, TRACY M B, et al. Measurements of store forces and moments and cavity pressures for a generic store in and near a box cavity at subsonic and transonic speeds[R]. NASA TM-4611, 1995.

[3] STALLINGS R L, FLOYD J, WILCOX B. Experimental cavity pressure distributions at supersonic speeds[R]. NASA TP-2683, 1987.

[4] TRACY M B, PLENTOVICH E B. Characterization of cavity flow fields using pressure data obtained in the Langley 0.3-meter transonic cryogenic tunnel[R]. NASA TM-4436, 1993.

[5] MAUREEN B, TRACY M B, PLENTOVICH E B. Cavity unsteady-pressure measurements at subsonic and transonic speeds[R]. NASA TP-3669, 1997.

[6] MURRAY N E, UKEILEY L S. Flow field dynamics in open cavity flows[R]. AIAA—2006-2428.

[7] ATVARS K, KNOWLES K, RITCHIE S A, et al. Experimental and computational investigation of an‘open’ transonic cavity flow[J]. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 2009, 223(4): 357-368.

[8] SHIEH C M, MORRIS P J. Comparison of two-and three-dimensional turbulent cavity flows[R]. AIAA—2001-0511.

[9] KIM H J, ARADAG S, KNIGHT D D. Two and three dimensional simulations of supersonic cavity flow[R]. AIAA—2006-2431.

[10] NAYYAR P, BARAKOS G N, BADCOCK K J, et al. Analysis and control of transonic cavity flow using DES and LES[R]. AIAA—2005-5267.

[11] ROWLEY C W, COLONIUS T, BASU A J. On self-sustained oscillations in two-dimensional compressible flow over rectangular cavities[J]. Journal of Fluid Mechanics, 2002, 455: 315-346.

[12] PENG S H. Simulation of turbulent flow past a rectangular open cavity using DES and unsteady RANS[R]. AIAA—2006-2827.

[13] LI Z S, HAMED A, BASU D. Numerical simulation of sidewall effects on the acoustic field in transonic cavity[R]. AIAA—2007-1456.

[14] 劉健. 采用 DES 類方法研究基本起落架非定常大范圍分離流動[D]. 北京:清華大學, 2011.

[15] 陳龍,伍貽兆,夏健. 基于DES的高雷諾數空腔噪聲數值模擬[J]. 計算力學學報,2011,28(5): 749-753.

CHEN Long, WU Yizhao, XIA Jian. High reynolds number cavity acoustic numerical simulation using DES[J]. Chinese Journal of Computational Mechanics, 2011,28(5): 749-753.

[16] 陳都. 基于脫體渦模擬方法的空腔流場計算研究[D]. 南京:南京航空航天大學, 2012.

[17] 譚玉婷. 空腔非定常流動特性的數值模擬研究[D]. 南京:南京航空航天大學,2012.

[18] 張寶兵. 空腔流動的機理模擬和控制[D]. 南京:南京航空航天大學,2012.

[19] 司海青,王同光. 邊界條件對三維空腔流動振蕩的影響[J]. 南京航空航天大學學報,2006,38(5):595-599.

SI Haiqing, WANG Tongguang. Influence of boundary condition on 3-D cavity flow-induced oscillations[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2006,38(5):595-599.

[20] PENG S H. M219 Cavity flow. DESider—A european effort on hybrid RANS-LES modelling[M]. Springer,2009,103:270-285

[21] CROOK S D, LAU T C W, KELSO R M. Three-dimensional flow within shallow, narrow cavities[J]. Journal of Fluid Mechanics, 2013, 735: 587-612.

[22] KHANAL B, KNOWLES K, SADDINGTON A J. Computational study of flowfield characteristics in cavities with stores[J]. Aeronautical Journal, 2011, 115(1173): 669-681.

[23] 余培汛,白俊強,郭博智. 剪切層形態對開式空腔氣動噪聲的抑制[J]. 振動與沖擊,2015,34(1):156-164.

YU Peixun, BAI Junqiang, GUO Bozhi, et al. Suppression of aerodynamic noise by altering the form of shear layer in open cavity[J]. Journal of Vibration and Shock, 2015,34(1):156-164.

[24] LAWSON S J, BARAKOS G N. Review of numerical simulations for high-speed, turbulent cavity flows[J]. Progress in Aerospace Sciences, 2011, 47(3): 186-216.

Numerical simulation of near-field cavity noise by URANS and DDES

LIU Jun1, YANG Dangguo2, WANG Xiansheng2, LUO Xinfu2

(1. State Key Laboratory of Aerodynamics, China Aerodynamics Research and Development Center, Mianyang 621000, China;2. High Speed Aerodynamics Research Institute, China Aerodynamics Research and Development Center, Mianyang 621000, China)

Numerical simulations of three-dimensional unsteady cavity flow at Mach number 0.85 were conducted by using the delayed detached eddy simulation (DDES) and the unsteady RANS (URANS) methods based on the SST turbulence model. The obtained static pressure, sound pressure level (SPL), and power spectral of pressure fluctuations at the cavity bottom by the both methods agreed well with experimental and other CFD results. Compared with URANS, DDES is more suitable for modeling flow details such as shear layer instability and fine scale turbulence structures, and also for high frequency pressure fluctuation capture. The simulated cavity flow can be identified as transitional flow by checking the time averaged flow characteristics. Meanwhile, complicated three-dimensional flow structures were also captured by both DDES and URANS, the effects of which were discussed.

unsteady Reynolds-averaged Navier-Stocks (URANS); delayed detached eddy simulation (DDES); cavity noise; pressure fluctuation; vortex structure

國家973項目(613240);自然科學基金(11402286);空氣動力學國家重點實驗室研究基金(SKLA20140302)

2015-06-18 修改稿收到日期:2015-10-14

劉俊 男,碩士,助理研究員,1986年10月生

楊黨國 男,博士,副研究員,1980年1月生

E-mail:yangdg-cardc@163.com

V211; TB84

A

10.13465/j.cnki.jvs.2016.20.025

猜你喜歡
結構方法
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
學習方法
論《日出》的結構
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
主站蜘蛛池模板: 免费av一区二区三区在线| 亚欧美国产综合| 国产凹凸视频在线观看| 99精品免费在线| 日韩高清在线观看不卡一区二区| 亚洲天堂在线免费| 二级特黄绝大片免费视频大片| 无码'专区第一页| 欧美性猛交一区二区三区 | 小说 亚洲 无码 精品| 成人国产精品网站在线看| 亚洲熟妇AV日韩熟妇在线| 手机在线看片不卡中文字幕| 一本大道视频精品人妻| 國產尤物AV尤物在線觀看| 国产区免费| 国产精品久线在线观看| 香蕉国产精品视频| 永久毛片在线播| 日韩小视频在线观看| 国产午夜在线观看视频| 国产无遮挡猛进猛出免费软件| 久久96热在精品国产高清| 在线视频亚洲欧美| 熟女视频91| 色综合中文综合网| 日本AⅤ精品一区二区三区日| 国产三级成人| 精品久久香蕉国产线看观看gif| www.国产福利| 亚洲视频a| 日韩麻豆小视频| a免费毛片在线播放| 日韩美女福利视频| 亚洲人精品亚洲人成在线| 91人妻日韩人妻无码专区精品| 中文字幕伦视频| 欧美日韩导航| 免费人成又黄又爽的视频网站| 婷婷在线网站| 国产h视频免费观看| 国产成人精品高清在线| 一区二区在线视频免费观看| 欧美精品导航| 超清无码一区二区三区| 99精品在线视频观看| 国产日韩欧美在线视频免费观看| 国产精品亚洲一区二区三区在线观看| 香蕉国产精品视频| 国产又粗又爽视频| 日本欧美午夜| 精品国产免费观看| 欧美一级夜夜爽www| 亚洲永久精品ww47国产| 国产99视频精品免费视频7| 日韩高清欧美| 2021国产精品自产拍在线| 成人国产精品一级毛片天堂| 热99re99首页精品亚洲五月天| 亚洲综合片| 欧洲成人免费视频| 亚洲第一成年人网站| 欧美成a人片在线观看| 2021天堂在线亚洲精品专区| 国产精品免费福利久久播放| 欧美日韩国产成人高清视频| 欧美特黄一级大黄录像| 97视频在线观看免费视频| 玖玖精品在线| 日本午夜精品一本在线观看| 尤物成AV人片在线观看| 日日摸夜夜爽无码| 又黄又湿又爽的视频| 精品无码一区二区三区在线视频| 欧美日一级片| 国产欧美精品专区一区二区| 91色爱欧美精品www| 久久久久88色偷偷| 久久大香伊蕉在人线观看热2 | 一级毛片在线直接观看| 97人人模人人爽人人喊小说| 亚洲大学生视频在线播放 |