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

基于層析PIV的湍流邊界層展向渦研究

2016-06-22 14:47:03王洪平魏潤杰王晉軍
實驗流體力學 2016年2期

王洪平, 高 琪,*, 魏潤杰, 王晉軍

(1. 北京航空航天大學流體力學研究所 流體力學教育部重點實驗室, 北京 100191; 2. 北京立方天地科技有限責任公司, 北京 100083)

基于層析PIV的湍流邊界層展向渦研究

王洪平1, 高 琪1,*, 魏潤杰2, 王晉軍1

(1. 北京航空航天大學流體力學研究所 流體力學教育部重點實驗室, 北京 100191; 2. 北京立方天地科技有限責任公司, 北京 100083)

層析粒子圖像測速;湍流邊界層;旋渦識別;展向渦

0 引 言

湍流中的擬序結構被認為是湍流產生和發展的關鍵,因此針對擬序結構的研究一直是湍流研究的一個重要方向。近期的湍流研究主要關注擬序結構和湍流統計量之間的關系。Smits[1]在關于高雷諾數壁

湍流的綜述性文章中將擬序結構主要分為4類,包括近壁區條帶,馬蹄渦,大尺度結構(Large-Scale Motions,LSMs)和超大尺度結構(Very-Large-Scale Motions,VLSMs)。這里馬蹄渦是對稱以及不對稱馬蹄渦、Ω形渦和拱形渦的統稱,這些渦一般由1個或2個流向渦腿和1個展向渦頭組成,展向渦頭的旋轉方向和湍流邊界層的平均剪切方向一致[2]。Adrian[3]用平面粒子圖像測速(Particle Image Velocimetry, PIV)實驗證明了在湍流邊界層外區存在馬蹄渦,同時這些馬蹄渦的渦腿附著在壁面上并沿流向排列形成了流向渦包結構。目前,這種附著馬蹄渦模型已被廣泛接受[4]。但是,隨著數值計算以及實驗手段的不斷改進和完善,有學者提出了不同的意見。Yang[5]通過馬蹄渦發生器產生馬蹄渦并觀察馬蹄渦向下游發展的情況,發現由渦產生的噴射(Ejection)現象能夠導致二次馬蹄渦的產生,并且能夠誘導產生反向展向渦。Pirozzoli[6]通過條件平均進一步闡明了湍流中渦管和剪切層的相互作用,發現渦管基本分布在剪切層上并推斷前者可能是Kelvin-Helmholtz不穩定性的產物。Lozano-Durán[7]在追蹤擬序結構發展變化的過程中發現并不是所有的壁面附著結構都產生于壁面。Eitel-Amor[8]重做了Zhou[9]的數值實驗,指出以前的研究[9-10]可能過高地估計了馬蹄渦存在的生命周期以及再生的重要性和魯棒性,在湍流耗散環境中馬蹄渦會很快衰減,不足以維持湍流大尺度運動。

通過上面的討論可知,馬蹄渦雖然能夠成功地描述湍流邊界層中的一些現象,但仍有學者對其在充分發展湍流中的重要性提出了疑問。在本文中,定義正展向渦和湍流邊界層中馬蹄渦頭的方向一致,用“p”表示(Prograde Spanwise Vortices),負展向渦用“r”表示(Retrograde Spanwise Vortices),與湍流邊界層平均剪切方向相反。眾所周知,在湍流邊界層中存在和馬蹄渦渦頭相反的負展向渦,而且總是和正展向渦成對出現[2,6,11]。隨著雷諾數的升高,負展向渦的作用將變得更加明顯[2]。目前,對負展向渦的產生機理主要有3種解釋:正展向渦誘導局部剪切失穩[5],馬蹄渦的相互合并[12]以及Ω形渦的一部分[11]。然而,由于旋渦結構復雜的三維特性,目前這些解釋仍缺少三維流場的進一步實驗驗證。

想要完美地解釋負展向渦的產生機理以及和正展向渦的關系,就必須得到時間解析的瞬時三維三分量(3D3C)速度場。隨著科學技術的發展和實驗手段的不斷完善,層析PIV技術[13-14]的出現使這種想法成為可能。從本質上講,層析PIV和平面PIV的原理相同,它通過多相機多視角成像,然后利用倍增代數重構技術(Multiplicative Algebraic Reconstruction Technique,MART)得到流場示蹤粒子的三維空間分布,進而互相關分析得到3D3C速度場。經過近10年的發展,層析PIV已經被用于測量高低雷諾數湍流邊界層[15-16]、零質量射流[17]、激波邊界層相互作用[18]、三角翼前緣渦[19]以及渦環撞壁實驗[20]等,得到了復雜流場的三維結構。

本文運用層析PIV技術測量了Reτ≈1700的湍流邊界層,并且從統計角度對正負展向渦進行了分析,希望對湍流邊界層中的擬序結構有更深入的理解。本文的組織結構如下:第1節和第2節分別介紹了層析PIV的具體實驗方法以及旋渦辨識的方法,這2部分是數據分析的基礎;第3節統計了正負展向渦沿法向的變化規律;第4節進一步討論了正負展向渦的空間位置對流場速度的影響;第5節討論了幾個典型的三維流場結構;最后給出總結和結論。

1 實驗方法

本實驗在北京航空航天大學大水洞完成,整個測量工況如圖 1所示,其中x,y,z分別代表流向、法向和展向。為了生成滿足要求的湍流邊界層,在大水洞中豎直安裝了一塊大小為7m×1m的有機玻璃板并且前緣為橢圓倒角,直徑為3mm的拌線布置在距平板前緣0.1m的位置。水洞底壁的寬度為1m,玻璃板被安裝在距側壁0.6m的位置,實驗時水位高度也為0.6m。測量體處在展向中心,距玻璃板前緣6.7m,下邊界距平板壁面的高度為5.5mm。流場中投放平均粒徑為20μm的空心玻璃微珠作為示蹤粒子,其密度為1.05×10-3g·mm-3。4臺分辨率為2456pixel×2048pixel的12位CCD相機呈“×”形布置在水洞一側,相機之間的夾角約為60°,各臺相機配置尼康f=45mm定焦體視鏡頭。為了減少由于視角引起的光學畸變,鏡頭的光圈f#被調到8和11之間,并且相機、鏡頭和物面滿足Scheimpflug條件。光源采用500mJ雙脈沖激光(鐳寶光電Vlite-500 PIV專用雙脈沖激光器),波長532nm,經過擴束光路被擴展成厚度為16mm的體光源。相機和激光之間通過同步器(北京立方天地SM-MicroPulse725同步器)進行控制,實驗時激光頻率為4Hz,2次曝光的時間間隔Δt=1.2ms。每臺相機采樣了300幅圖像,采樣時長約為38s,最終得到了150個速度場。

圖1 層析PIV實驗工況示意圖

層析PIV實驗前需要標定物理空間和像素空間的關系。本文采用多項式標定(x,z方向采用三次多項式,y向采用一次多項式),標定板在測量體厚度范圍內以2mm等間距移動并成像。初始得到的映射函數需要進行自修正[21],將映射函數的誤差控制在0.1pixel以內,具體處理過程可以參考文獻[22]。在層析PIV中,圖像前處理也是一個非常重要的環節,可以很大程度上提升重構的精度和效率。本文對原始圖像進行了時序最小背景剔除和空間滑動背景剔除,最終的圖像進行了3×3的高斯平滑,高斯核σ=0.5。經過圖像前處理后,粒子濃度約為0.06ppp(particle per pixel)。結合映射函數,空間粒子灰度場通過MART算法重構得到。本文進行10次MART迭代,重構體的大小為80mm×16mm×45mm,圖像放大率為0.055mm/pixel,因此重構體素(voxel)的個數為1454voxel×290voxel×818voxel。

重構后的粒子灰度場直接通過互相關分析確定粒子的位移。本文采用多層迭代體變形算法[23],最終窗口大小為48voxel×48voxel×48voxel,重疊率為75%。為了提高計算效率,第1層互相關采用快速傅里葉變換算法(Fast Fourier Transform, FFT),其余互相關在得到速度參考矢量后進行窗口預偏置,然后采用小位移直接互相關算法。在體變形時只采用線性插值,這已經能滿足計算精度的要求。另外,每一層互相關得到的速度場需要進行歸一化壞點剔除[24]和高斯濾波,才能成為下一層互相關迭代的參考速度矢量。最終得到的速度場在剔除壞點的基礎上進行了高斯平滑,平滑窗的大小為3×3×3,平滑強度σ=1。

流場的品質和參數通過激光多普勒測速(LaserDopplerVelocimetry,LDV)進行標定,表 1給出了流場的一些基本信息。實驗時水溫為18℃,對應的運動粘性系數ν為1.055mm2/s。當自由來流速度U∞約為0.413m/s時,流場的湍流度為0.97%。通過Musker方法[25]計算得到壁面摩阻速度uτ為15.57mm/s,邊界層厚度δ為120mm,對應的基于壁面摩阻速度和動量厚度的雷諾數分別為1769和5176。根據公式用壁面粘性尺度(WallUnits,WU)對物理量進行無量綱化。

(1)

因此,粘性尺度的測量體大小為1176WU×235WU×661WU。速度場的法向測量范圍為100

表1 湍流邊界層的基本參數

2 旋渦及渦心辨識方法

2.1 三維流場的旋渦辨識

在湍流邊界層中,存在渦量的地方可能存在旋渦,也可能是剪切層,因此,必須從湍流背景流場中提取真正的旋渦。根據Robinson[26]關于旋渦的定義,旋渦通常會以一定的速度對流,需要考慮伽利略不變性(Galilean Invariant),但在湍流邊界層中旋渦的對流速度往往與所處位置的流動結構有關。目前的旋渦識別方法大多通過分析局部速度梯度張量來實現,例如Q準則、Δ準則、λ2準則以及λci準則,這些判據之間存在一定的聯系,某些情況下是相互等價的[27]。本文采用基于當地脈動速度u的λci準則,即旋渦強度判據,該判據用局部速度梯度的復特征值的虛部來表示旋轉的強度。由于分析的是脈動速度,湍流邊界層當地時均速度U(僅考慮流向速度)已經從流場中剔除。考慮到PIV數據存在噪音以及空間分辨率的問題,速度梯度的計算往往存在誤差。同時,旋渦強度在湍流邊界層中并不是均勻分布,而是與法向高度y+有關。因此,閾值的設置不僅要考慮誤差,還要考慮y+的影響。根據其他學者的研究成果[2,11],本文閾值的設定如公式(2)所示。

(2)

(3)

2.2 三維流場中展向渦位置

本文對展向渦的定位采用了類似文獻[28]的方法,即找出展向渦所對應的連通區域,以該連通區域的質心作為三維流場中展向渦的位置。具體計算過程為:

3 展向渦沿法向分布規律

3.1 展向渦強度

(a) 代表平均值

(b) 代表脈動值

3.2 展向渦對流速度

(a)

(b)

圖3 正負展向渦的流向對流速度(a)和法向對流速度(b)沿高度的變化。(a)中的黑色虛線代表速度剖面,(b)中的黑色虛線代表V+=0。

Fig.3 (a) The mean streamwise velocity of prograde and retrograde spanwise vortices as a function ofy+,theblackdottedlinerepresentsthemeanvelocityprofile. (b)Themeanwall-normalvelocityofprogradeandretrogradespanwisevorticesasafunctionofy+,theblackdottedlinerepresentsV+=0.

(a) 正展向渦 (b) 負展向渦 (c) 全流場

圖4 流向速度和法向速度的聯合概率密度分布

Fig.4 The joint probability density function of streamwise velocity and wall-normal velocity. (a) prograde spanwise vortices, (b) retrograde spanwise vortices, (c) the whole flow field.

4 正負展向渦的空間關系

通過第3節的分析可知:正負展向渦在法向的區別主要集中在強度、空間流場結構以及對流速度上,本節主要討論正負展向渦在x-y平面內的空間關系。以正展向渦的渦心為中心,在半徑為200WU的區域內統計周圍正負展向渦出現的概率,在展向z允許有±2Δz+(±19.4WU)的偏差。在統計過程中,同一個展向渦如果同時被不同的中心渦選中,則選擇其中距離最近的一個。圖5為統計得到的結果,圖5(a)表示中心為正展向渦時周圍正展向渦的分布情況,而圖5(b)表示中心是正展向渦,負展向渦在其周圍的分布情況,橫縱坐標分別表示離中心渦的流向和法向距離,云圖代表周圍渦出現的概率,顏色越紅說明在該位置出現的概率越高。在圖5中,正展向渦基本上沿流向水平排列,間距約為60WU,由于樣本數量有限導致統計的分辨率較低。對于圖5(b)中的負展向渦,基本上分布在正展向渦的流向和法向。Natrajan[11]在文章中展示了法向分布的這種情況,但并沒有展示流向分布這種情況,在實際流場中,后者往往也是存在的,其出現的概率幾乎和法向出現的概率相等。另外,仔細觀察負展向渦法向出現的位置可以發現:在正展向渦底部時其偏向上流,而在正展向渦上部時偏向下游,這與Natrajan[11]的結果一致。

進一步分析正負展向渦的空間位置對流場速度的影響。如圖6(a)所示,以正展向渦為中心,負展向渦的位置可以用他們之間的連線與流向的夾角θ表示(半徑變化較小,一般在100個粘性尺度以內,如圖5(b)所示)。θ的取值范圍為0°到360°,180°和0°表示負展向渦在正展向渦的上下游,而90°和270°表示在上層和下層。連線中心的脈動速度可以由插值得出,用該速度可以比較容易和準確地判斷出流場速度與展向渦空間位置的關系。圖6(b)給出了連線中心的脈動速度(u+,v+)隨角度θ的變化規律,圓點代表流向速度u+,三角代表法向速度v+。隨著角度的變化,流向速度呈現出正弦曲線的變化規律,而法向速度表現出與之相對的余弦變化趨勢。當負展向渦分布在正展向渦的流向,即θ為0°或180°時,流向速度的平均值幾乎為0,而法向速度的平均值卻最大,約為1.2;當負展向渦分布在正展向渦的法向時,即θ為90°或270°時,流向速度表現出高低速區域,而法向速度卻維持在0左右,如圖6(a)中的紅藍箭頭所示。該規律完全背離了圖4所描述的速度之間的依賴關系,高低速區的產生無法簡單地用掃掠和噴射來解釋。需要注意的是,這種和渦結構有很強關系的高低速區域的尺度都比較小,與圖5中展向渦之間的距離相當。一般來說,噴射產生于近壁區并逐漸向上運動,相反,掃掠產生于較高位置朝壁面運動[7]。這種大范圍相對運動容易產生大尺度結構,而中小尺度的流動結構更多是由于不穩定性和旋渦運動造成的。

(a)

(b)

Fig.5 The distribution of prograde (a) and retrograde (b) spanwise vortices around a prograde vortex

(a)

圖6 (a)高低速區域與正負展向渦位置關系示意圖,(b)速度與正負展向渦空間位置的關系

Fig.6 (a) The schematic diagram about the relationship between the low (high) momentum region and the spatial arrangement of spanwise vortices, (b) The corresponding variations inu+(v+)andθ

5 三維流場

6 結 論

本文通過層析PIV測量了Reτ≈1700的湍流邊界層,并統計分析了100

(1) 對于旋渦強度,正負展向渦表現出不同的性質。隨著法向高度的增加正展向渦逐漸變弱而負展向渦基本保持不變。

(2) 流向速度和法向速度呈現出很強的相關性,導致流場中的大尺度高低速區域成為雷諾應力-uv的主要來源。

(3) 正展向渦在流向呈現出等間距分布,而負展向渦大多分布在正展向渦的流向和法向,但并沒有統計這種分布在流向持續的長度。

(4) 在小尺度范圍內(約為100WU),高低速區域并不受噴射和掃掠機制的控制,而是更多和渦結構的空間位置以及局部不穩定性有關。

[1] Smits A J, Mckeon B J, Marusic I. High-Reynolds number wall turbulence[J]. Annual Review of Fluid Mechanics, 2011, 43(1): 353-375.

[2] Wu Y, Christensen K T. Population trends of spanwise vortices in wall turbulence[J]. J Fluid Mech, 2006, 568: 55-76.

[3] Adrian R J, Meinhart C D, Tomkins C D. Vortex organization in the outer region of the turbulent boundary layer[J]. J Fluid Mech, 2000, 422: 1-54.

[4] Adrian R J. Hairpin vortex organization in wall turbulence[J]. Physics of Fluids, 2007, 19(4).

[5] Yang W, Meng H, Sheng J. Dynamics of hairpin vortices generated by a mixing tab in a channel flow[J]. Exp Fluids, 2001, 30(6): 705-722.

[6]Pirozzoli S. Flow organization near shear layers in turbulent wall-bounded flows[J]. Journal of Turbulence, 2011: N41.

[7] Lozano-Durán A, Jiménez J. Time-resolved evolution of coherent structures in turbulent channels: characterization of eddies and cascades[J]. J Fluid Mech, 2014, 759: 432-471.

[8] Eitel-Amor G, ?rlü R, Schlatter P, et al. Hairpin vortices in turbulent boundary layers[J]. Physics of Fluids, 2015, 27(2): 025108.

[9] Zhou J, Adrian R J, Balachandar S, et al. Mechanisms for generating coherent packets of hairpin vortices in channel flow[J]. J Fluid Mech, 1999, 387: 353-396.

[10] Kim K, Sung H J, Adrian R J. Effects of background noise on generating coherent packets of hairpin vortices[J]. Physics of Fluids, 2008, 20(10): 105107.

[11] Natrajan V K, Wu Y, Christensen K T. Spatial signatures of retrograde spanwise vortices in wall turbulence[J]. J Fluid Mech, 2007, 574: 155-167.

[12] Adrian R J, Balachandar S, Lin Z C. Spanwise growth of vortex structure in wall turbulence[J]. KSME International Journal, 2001, 15(12): 1741-1749.

[13] Gao Q, Wang H P, Shen G X. Review on development of volumetric particle image velocimetry[J]. Chinese Science Bulletin, 2013, 58(36): 4541-4556.

[14] Elsinga G E, Scarano F, Wieneke B, et al. Tomographic particle image velocimetry[J]. Exp Fluids, 2006, 41(6): 933-947.

[15] Elsinga G E, Adrian R J, Van Oudheusden B W, et al. Three-dimensional vortex organization in a high-Reynolds-number supersonic turbulent boundary layer[J]. J Fluid Mech, 2010, 644: 35-60.

[16] Gao Q. Evolution of eddies and packets in turbulent boundary layers[D]. Minnesota: The University of Minnesota, 2011.

[17] 高琪, 王洪平. 層析PIV技術及其合成射流測量[J]. 中國科學:技術科學, 2013, (7): 828-835.

Gao Q, Wang H P. The tomographic PIV and it's application on synthetic jet measurement[J]. Science China: Technological Sciences, 2013, (7): 828-835.

[18] Humble R A, Elsinga G E, Scarano F, et al. Three-dimensional instantaneous structure of a shock wave/turbulent boundary layer interaction[J]. J Fluid Mech, 2009, 622: 33-62.

[19] 王成躍, 高琪, 魏潤杰, 等. 三角翼前緣渦三維流動顯示和層析PIV測量[C]. 第八屆全國流體力學學術會議, 2014.

Wang C Y, Gao Q, Wei R J, et al. 3D flow visualization and tomographic particle image velocimetry for vortex breakdown of a delta wing[C]. The 8th National Congress of Fluid Mechanics, 2014.

[20] Violato D, Ianiro A, Cardone G, et al. Investigation on circular and chevron impinging jets by IR thermography and time-resolved tomographic PIV[C]. ASME-JSME-KSME 2011 Joint Fluids Engineering Conference, 2011.

[21] Wieneke B. Volume self-calibration for 3D particle image velocimetry[J]. Exp Fluids, 2008, 45(4): 549-556.

[22] Gao Q, Wang H P, Wang J J. A single camera volumetric particle image velocimetry and its application[J]. Science China: Technological Sciences, 2012, 55(9): 2501-2510.

[23] Scarano F. Iterative image deformation methods in PIV[J]. Meas Sci Technol, 2002, 13(1): R1-R19.

[24] Westerweel J, Scarano F. Universal outlier detection for PIV data[J]. Exp Fluids, 2005, 39(6): 1096-1100.

[25] Kendall A, Koochesfahani M. A method for estimating wall friction in turbulent wall-bounded flows[J]. Exp Fluids, 2008, 44(5): 773-780.

[26] Robinson S K. Coherent motions in the turbulent boundary layer[J]. Annual Review of Fluid Mechanics, 1991, 23(1): 601-639.

[27] Chakraborty P, Balachandar S, Adrian R J. On the relationships between local vortex identification schemes[J]. J Fluid Mech, 2005, 535: 189-214.

[28] Del Alamo J C, Jiménez J, Zandonade P, et al. Self-similar vortex clusters in the turbulent logarithmic region[J]. J Fluid Mech, 2006, 561: 329-358.

(編輯:張巧蕓)

Study of spanwise vortices in turbulent boundary layer flow based on tomographic PIV

Wang Hongping1, Gao Qi1,*, Wei Runjie2, Wang Jinjun1

(1. Fluid Mechanics Key Laboratory of Ministry of Education, Beijing University of Aeronautics and Astronautics, Beijing 100191, China; 2. MicroVec., Inc, Beijing 100083, China)

tomographic PIV;turbulent boundary layer;vortex identification;spanwise vortices

1672-9897(2016)02-0059-08

10.11729/syltlx20150086

2015-06-08;

2015-07-30

國家自然科學基金(11472030,11327202,11490552)

WangHP,GaoQ,WeiRJ,etal.StudyofspanwisevorticesinturbulentboundarylayerflowbasedontomographicPIV.JournalofExperimentsinFluidMechanics, 2016, 30(2): 59-66. 王洪平, 高 琪, 魏潤杰, 等. 基于層析PIV的湍流邊界層展向渦研究. 實驗流體力學, 2016, 30(2): 59-66.

V211.7;O

A

王洪平(1987-),男,四川廣元人,博士。研究方向:實驗流體力學。通信地址: 北京航空航天大學航空科學與工程學院流體力學研究所(100191)。E-mail: h.p.wang@foxmail.com。

*通信作者E-mail:qigao@buaa.edu.cn

主站蜘蛛池模板: 国产精品自拍露脸视频| 亚洲美女视频一区| 国产福利免费视频| 久久精品亚洲热综合一区二区| 成人福利免费在线观看| 国产精品女人呻吟在线观看| 国产91av在线| 在线高清亚洲精品二区| 国产精品亚欧美一区二区| 亚洲黄网视频| 日韩在线成年视频人网站观看| 伊人久久大线影院首页| 欧类av怡春院| 亚洲欧洲自拍拍偷午夜色无码| 最新国产麻豆aⅴ精品无| 好紧好深好大乳无码中文字幕| 丝袜久久剧情精品国产| 日本一本正道综合久久dvd| 国产精品一线天| 成年人福利视频| 久久久受www免费人成| 99视频只有精品| 91亚洲精品第一| 青青草原国产一区二区| 91在线无码精品秘九色APP| 特级毛片免费视频| 日韩中文精品亚洲第三区| 天天色综网| 国产一区二区三区在线精品专区 | 欧美激情视频一区二区三区免费| 国产精品刺激对白在线| 国产99精品久久| 欧美天天干| 亚洲va在线观看| 亚洲综合一区国产精品| 久久影院一区二区h| 久久国产av麻豆| 婷婷六月综合网| 狠狠做深爱婷婷久久一区| 欧美日韩高清| 国产第一页第二页| 国产黄网永久免费| 精品无码人妻一区二区| 午夜电影在线观看国产1区| 性做久久久久久久免费看| 国产在线精品人成导航| 国产激情第一页| 中文字幕精品一区二区三区视频| 波多野结衣中文字幕一区二区| 欲色天天综合网| 亚洲av综合网| 亚洲精品欧美重口| 欧美午夜小视频| 亚洲人在线| 国产av色站网站| 欧美天堂久久| 黄色网址免费在线| 国产AV毛片| 大香伊人久久| 亚洲婷婷六月| 久草网视频在线| 国产成人精品视频一区二区电影 | 国产凹凸视频在线观看| 欧美精品成人一区二区在线观看| 自拍欧美亚洲| 一区二区三区四区日韩| 香蕉综合在线视频91| 2021国产在线视频| 99ri精品视频在线观看播放| 国产尤物jk自慰制服喷水| 欧美综合区自拍亚洲综合天堂| 欧美丝袜高跟鞋一区二区| 久久公开视频| 亚洲第一色网站| 国产农村妇女精品一二区| 国产网站免费看| 色爽网免费视频| 亚洲Av综合日韩精品久久久| 亚洲男人的天堂网| 日韩人妻精品一区| 中文字幕无线码一区| 91探花在线观看国产最新|