亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        應用加權(quán)緊致非線性格式的VFE-2鈍前緣三角翼轉(zhuǎn)捩模擬*

        2016-10-10 02:31:19王光學王圣業(yè)王東方鄧小剛
        國防科技大學學報 2016年4期
        關鍵詞:三角翼迎角吸力

        王光學,王圣業(yè),王東方,鄧小剛

        (1.國防科技大學 航天科學與工程學院, 湖南 長沙 410073;2.中山大學 物理學院, 廣東 廣州 510006)

        ?

        應用加權(quán)緊致非線性格式的VFE-2鈍前緣三角翼轉(zhuǎn)捩模擬*

        王光學1,2,王圣業(yè)1,王東方1,鄧小剛1

        (1.國防科技大學 航天科學與工程學院, 湖南 長沙410073;2.中山大學 物理學院, 廣東 廣州510006)

        為研究前緣轉(zhuǎn)捩對鈍前緣三角翼渦結(jié)構(gòu)的影響,采用高階精度加權(quán)緊致非線性格式和γ-Reθ轉(zhuǎn)捩模型對VFE-2中等半徑鈍前緣三角翼進行數(shù)值模擬。將計算結(jié)果與試驗結(jié)果進行詳細對比,結(jié)果表明:鈍前緣三角翼的前緣分離渦發(fā)生在翼尖下游,在特定雷諾數(shù)下其具體發(fā)生位置受轉(zhuǎn)捩因素影響,采用全湍流模型計算會推遲分離,而耦合轉(zhuǎn)捩模型后的計算結(jié)果和試驗結(jié)果吻合良好。運用耦合轉(zhuǎn)捩模型方法,對鈍前緣三角翼渦結(jié)構(gòu)隨迎角變化進行模擬。計算結(jié)果與試驗結(jié)果吻合,表明在較小的迎角下,前緣不會產(chǎn)生分離誘導渦;隨迎角不斷增大,分離誘導渦在三角翼后緣附近產(chǎn)生并向上游移動。

        轉(zhuǎn)捩模型;三角翼;高精度格式;加權(quán)緊致非線性格式;分離誘導渦

        現(xiàn)代戰(zhàn)斗機和導彈多采用三角翼布局,以獲得良好的飛行品質(zhì)和機動性能。在不大的迎角下,三角翼在上翼面形成前緣分離渦。前緣分離渦能提供渦升力,提高飛行器的氣動特性,但同時也使流動變得復雜。

        20世紀50年代以來對三角翼氣動特性的研究一直持續(xù)至今,其中較有影響力的是國際渦流動試驗(Vortex Flow Experiment, VFE)。20世紀80年代開展的第一次國際渦流動試驗(VFE-1),主要關注65°尖前緣三角翼構(gòu)型[1]。該構(gòu)型產(chǎn)生“典型”的渦結(jié)構(gòu),包括起支配作用的主渦和弱的二次渦。對尖前緣三角翼構(gòu)型的數(shù)值模擬也同時發(fā)展。早期,代數(shù)湍流模型以及低階精度數(shù)值方法,無法較好模擬該類流動。而近年來隨著一方程和二方程湍流模型尤其是高精度數(shù)值方法的發(fā)展,尖前緣三角翼流動已能夠被很好地模擬[2-3]。

        隨著高超聲速飛行器的研制,鈍前緣三角翼受到廣泛關注。鈍前緣三角翼的渦結(jié)構(gòu)不同于尖前緣,分離點在翼尖下游某處,其具體位置受多種因素影響,包括雷諾數(shù)、迎角和前緣半徑等[4-5]。2001年Hummel[6]提出開展第二次國際渦流動試驗(VFE-2),其主要目的之一就是為鈍前緣三角翼的計算提供準確的驗證數(shù)據(jù)。近年來,以該試驗為依據(jù),國際上許多學者對65°鈍前緣三角翼構(gòu)型開展了數(shù)值模擬研究。Schutte[7]基于TAU代碼,采用單方程(Spalart-Allmaras, SA)模型和Wilcoxk-ω模型對該構(gòu)型進行了數(shù)值計算,并表明SA模型比k-ω模型更接近試驗值;Fritz[8]基于FLOWer代碼,采用了SA模型、Wilcoxk-ω模型和雷諾應力模型(Reynolds Stress Model, RSM)模型,而結(jié)論卻與前者相反,即k-ω模型較好;Crivellini等[9]對比兩者結(jié)果,發(fā)現(xiàn)很難找出造成兩者結(jié)論相反的因素是湍流模型還是數(shù)值方法,并指出應該采用高階精度數(shù)值方法以盡量避免數(shù)值方法帶來的不確定因素。另外,鈍前緣三角翼前緣分離渦的形成常常伴隨轉(zhuǎn)捩的發(fā)生。Fritz[8]采用固定轉(zhuǎn)捩技術成功模擬了鈍前緣三角翼渦結(jié)構(gòu),但轉(zhuǎn)捩位置固定使該方法具有很大的局限性。

        1 高精度數(shù)值模擬方法

        1.1計算模型

        圖1 VFE-2三角翼模型Fig.1 Model configuration of VFE-2 delta wing

        采用VFE-2 65°后掠角三角翼模型。如圖1所示,該模型分為4部分:前緣、平板部分、后緣及整流罩[10]。前緣部分提供4種可替換外形,重點關注中等半徑鈍前緣外形,同時采用尖前緣外形做對比研究。

        1.2計算網(wǎng)格

        計算網(wǎng)格自主生成,網(wǎng)格節(jié)點數(shù)約600萬,網(wǎng)格結(jié)構(gòu)為多塊對接網(wǎng)格,共27塊。圖2給出了表面網(wǎng)格。網(wǎng)格拓撲采用C-H型,以避免翼尖奇性軸的產(chǎn)生。計算區(qū)域的遠場邊界取為50倍根弦長。壁面的第一層網(wǎng)格達到了1.0×10-6弦長,網(wǎng)格在背風區(qū)和各個剪切層附近均進行了適當?shù)募用?,以保證分離區(qū)、附面層內(nèi)和剪切層的數(shù)值模擬精度。

        圖2 三角翼計算網(wǎng)格Fig.2 Computational mesh for delta wing

        1.3高精度數(shù)值方法簡介

        加權(quán)緊致非線性格式(Weighted Compact Nonlinear Scheme, WCNS)首先由Deng等[11]在2000年提出。其后,不同學者[12-14]發(fā)展了多種形式的WCNS,并在廣泛的流動問題中證明了格式的優(yōu)勢。采用WCNS系列中的一種典型五階顯式離散格式WCNS-E-5。WCNS-E-5由于其低耗散、高魯棒和優(yōu)秀的自由流和渦保持特性,被廣泛應用于各種實際流動問題的高精度數(shù)值模擬中。時間格式方面,均為定常流動,采用隱式上下對稱高斯塞德爾(Lower Upper Symmetric Gauss Seidel, LU-SGS)方法。

        對于湍流模型,采用Menter′sk-ωSST模型,其守恒形式為:

        (1)

        (2)

        其中源項及系數(shù)具體參見文獻[15]。

        對于前緣轉(zhuǎn)捩問題,采用γ-Reθ模型,其守恒形式為:

        (3)

        (4)

        其中源項、系數(shù)以及與SST模型的耦合參見文獻[16]。

        另外需要指出的是,在計算網(wǎng)格導數(shù)及雅克比時,平臺中采用了滿足幾何守恒律的對稱守恒網(wǎng)格導數(shù)方法(Symmetrical Conservative Metric Method, SCMM)[17],有利于提高高精度有限差分方法的魯棒性并減小數(shù)值誤差。

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

        2.1尖和鈍前緣三角翼渦結(jié)構(gòu)對比

        鈍前緣三角翼的渦結(jié)構(gòu)與傳統(tǒng)尖前緣三角翼的有很大差異。圖3為Ma=0.4,Re=6×106,α=13.3°條件下尖前緣和中等半徑鈍前緣三角翼渦結(jié)構(gòu)對比。計算均采用全湍流模型,時間推進采用LU-SGS方法。圖3左側(cè),尖前緣三角翼展現(xiàn)了典型的分離誘導前緣渦的特點:來流從翼尖開始分離形成主渦,主渦在翼面再附后卷起外側(cè)附著流形成二次渦。圖3右側(cè),鈍前緣三角翼渦結(jié)構(gòu)較為復雜。來流流過鈍翼尖并未發(fā)生分離,隨著下游半展長的增加而前緣半徑恒定,機翼在流向上變得越來越尖。前緣相對鈍度減小,使來流在翼尖下游某處發(fā)生分離。

        圖3 尖和鈍前緣三角翼流動對比Fig.3 Contrast of sharp-edge and blunt-edge flows

        對于鈍前緣三角翼,前緣可被分為3部分:附著區(qū)、分離區(qū)和過渡區(qū)。附著區(qū),來流未發(fā)生分離;分離區(qū),來流在前緣分離產(chǎn)生典型的分離誘導渦,并伴隨二次渦結(jié)構(gòu);過渡區(qū),來流在前緣尚未產(chǎn)生分離,但在下游,邊界層受逆壓梯度(前緣吸力)影響而分離形成內(nèi)渦,該渦強度較弱,在接近尾緣處受主渦影響而消失。

        圖4為尖前緣和中等半徑鈍前緣三角翼典型站位處的壓力分布計算值與試驗值對比。試驗主要在NASA Langley研究中心的LTPT風洞[10]中進行,包括測力和測壓試驗。在x/cr=0.2處,尖前緣三角翼的吸力峰位于80%當?shù)卣归L處,即渦核位置,并且計算值與試驗吻合很好,而鈍前緣三角翼的吸力峰位于駐點處,流動為附著流。在x/cr=0.4處,鈍前緣三角翼已產(chǎn)生分離誘導渦,而計算結(jié)果中尚未出現(xiàn)吸力峰,即計算的前緣分離點推遲了。在x/cr=0.6處,鈍前緣三角翼出現(xiàn)了兩個吸力峰,其中內(nèi)渦強度較弱,位于60%當?shù)卣归L處;計算結(jié)果中同樣出現(xiàn)了兩個吸力峰,但均比試驗結(jié)果靠外。在x/cr=0.8和x/cr=0.95處,鈍前緣三角翼的內(nèi)渦強度變得更弱并逐漸消失;而計算的吸力峰仍比試驗值靠外。以上典型站位處的壓力分布對比表明,采用全湍流模型可以成功模擬尖前緣三角翼渦結(jié)構(gòu),這與文獻[2-3]的結(jié)論一致,但在模擬鈍前緣三角翼時出現(xiàn)了分離推遲現(xiàn)象。

        前緣分離渦可提高機翼上表面吸力,即產(chǎn)生渦升力。而前緣鈍度增加會使前緣分離渦推后和減弱。試驗測得尖、鈍前緣的法向力系數(shù)CN分別為0.570和0.517,可以看出,前緣鈍度增加將使機翼升力減小。另外,計算所得尖、鈍前緣的法向力系數(shù)CN分別為0.545和0.461,尖前緣三角翼計算結(jié)果和試驗值吻合較好,而鈍前緣三角翼計算結(jié)果偏小,仍是由分離推遲造成。

        2.2轉(zhuǎn)捩對鈍前緣三角翼渦結(jié)構(gòu)影響

        Fritz[8]指出鈍前緣三角翼的前緣分離對渦黏性以及數(shù)值黏性十分敏感,在采用全湍流模型(包括SAE模型、Wilcox′sk-ω模型和RSM模型)模擬時,均出現(xiàn)了分離推遲的現(xiàn)象。

        本節(jié)研究Ma=0.4,Re=6×106和α=13.3°條件下,有/無轉(zhuǎn)捩模型對模擬鈍前緣三角翼渦結(jié)構(gòu)的影響。圖5為有/無轉(zhuǎn)捩模型時三角翼表面壓力分布對比??梢钥吹剑徊捎棉D(zhuǎn)捩模型所計算的分離發(fā)生位置比采用轉(zhuǎn)捩模型的明顯靠下游,表明了鈍前緣三角翼的分離發(fā)生位置對于渦黏性十分敏感。

        圖4 尖前緣和鈍前緣三角翼典型站位處表面壓力計算值與試驗值對比Fig.4 Comparisons of surface pressure with experiments for sharp-edge and blunt-edge delta wing at typical stations

        圖5 有/無轉(zhuǎn)捩模型三角翼表面壓力云圖對比Fig.5 Contrast of surface pressure contours of with and without transition model

        圖6為x/cr分別為0.2,0.4,0.6,0.8和0.95五個典型站位處表面壓力計算值和試驗值的對比。試驗表明,前緣分離渦的起始位置在x/cr為0.2~0.4。在x/cr=0.2處,前緣并未產(chǎn)生分離渦,是否采用轉(zhuǎn)捩模型結(jié)果并無差別。在x/cr=0.4處,不采用轉(zhuǎn)捩模型的壓力分布中未出現(xiàn)吸力峰,即前緣未發(fā)生分離;而采用轉(zhuǎn)捩模型,前緣開始出現(xiàn)分離誘導渦。在x/cr分別為0.6,0.8和0.95三個站位處,采用轉(zhuǎn)捩模型的計算結(jié)果和試驗吻合很好,而不采用轉(zhuǎn)捩模型,計算的吸力峰靠外。

        圖6 有/無轉(zhuǎn)捩模型,三角翼典型站位處壓力分布計算值與試驗值對比Fig.6 Comparisons of surface pressure obtained with experiments for blunt-edge delta wing at typical stations of with and without transition model

        綜合圖5和圖6表明,基于高精度WCNS,通過耦合轉(zhuǎn)捩模型可很好模擬鈍前緣三角翼渦結(jié)構(gòu),而采用全湍流模型會延遲主渦的發(fā)生。分離推遲現(xiàn)象是由于全湍流模型在三角翼前緣層流區(qū)高估渦黏性造成的。由于γ-Reθ模型并不模擬實際的物理過程[16],而相關試驗也未研究前緣轉(zhuǎn)捩問題,因此前緣轉(zhuǎn)捩對分離渦影響的物理機理將在今后采用大渦模擬方法進行研究。

        2.3鈍前緣三角翼渦結(jié)構(gòu)隨迎角變化

        在特定的來流雷諾數(shù)和馬赫數(shù)下,鈍前緣三角翼的分離渦發(fā)生位置以及整體渦結(jié)構(gòu)主要受迎角影響。本節(jié)采用上節(jié)耦合轉(zhuǎn)捩模型的方法,探究鈍前緣三角翼渦結(jié)構(gòu)隨迎角變化的規(guī)律。為與風洞試驗對比,計算條件選擇:Ma=0.4,Re=3×106,α分別為10.2°,13.3°和15.3°。風洞試驗在德國宇航中心的TWG風洞中完成[18]。

        圖7(a)為WCNS的下轉(zhuǎn)捩模型計算得到的三角翼表面壓力分布,圖7(b)為試驗采用壓敏漆技術得到表面壓力分布[18]。計算結(jié)果與試驗符合得很好。在迎角為10.2°時,整個三角翼前緣均為附著流動,并未產(chǎn)生分離誘導渦,僅在接近尾緣處存在較弱的內(nèi)渦。迎角增大到13.3°時,在x/cr=0.4附近前緣流動出現(xiàn)分離,并形成較強的分離誘導渦。迎角繼續(xù)增大到15.3°時,前緣分離點提前至x/cr=0.2附近。對于中等半徑鈍前緣三角翼,隨著迎角增大,前緣分離誘導渦從無到有,逐漸向上游移動并不斷增強,使渦升力逐漸提高。同時計算的法向力系數(shù)CN分別為0.336,0.468和0.548,也印證了該結(jié)論。

        (a) 計算流體力學方法(a) Computational fluid dynamics

        (b) 壓敏漆方法(b) Pressure sensitive paint圖7 不同迎角下鈍前緣三角翼表面壓力云圖Fig.7 Surface pressure contours of blunt-edge delta wing at variation of the angle-of-attack

        3 結(jié)論

        1)尖前緣三角翼的分離誘導渦發(fā)生在翼尖,基于WCNS的全湍流SST模型能很好模擬;鈍前緣三角翼的分離誘導渦發(fā)生在翼尖下游某處,對渦黏性十分敏感,采用基于WCNS的全湍流模型會出現(xiàn)分離推遲現(xiàn)象。

        2) 基于高精度WCNS,在SST湍流模型上耦合γ-Reθ轉(zhuǎn)捩模型,可以很好模擬中等半徑鈍前緣三角翼的流場。

        3) 在特定雷諾數(shù)和馬赫數(shù)下,中等半徑鈍前緣三角翼的前緣分離渦的產(chǎn)生存在臨界迎角。小于該迎角,前緣不會產(chǎn)生分離誘導渦;超過該迎角,分離誘導渦在三角翼后緣附近產(chǎn)生,并隨迎角增大而向上游移動,同時法向力也隨之增大。

        References)

        [1]Drougge G.The international vortex flow experiment for computer code validation [C]//Proceedings of ICAS-Proceedings, Jerusalem, 1988, 1: 9-1-9-23.

        [2]王光學, 鄧小剛, 劉化勇, 等. 高階精度格式WCNS在三角翼大攻角模擬中的應用研究[J]. 空氣動力學學報, 2012, 30(1): 28-33.

        WANG Guangxue, DEND Xiaogang, LIU Huayong, et al. Application of high-order scheme(WCNS) at high angles of incidence for delta wing [J]. Acta Aerodynamica Sinica, 2012, 30(1): 28-33.(in Chinese)

        [3]王光學, 鄧小剛, 王運濤, 等. 三角翼渦破裂的高精度數(shù)值模擬[J]. 計算物理, 2012, 29(4): 489-494.

        WANG Guangxue, DEND Xiaogang, WANG Yuntao, et al. High-order numerical simulation of vortex breakdown on delta wing [J]. Chinese Journal of Computational Physics, 2012, 29(4): 489-494.(in Chinese)

        [4]Luckring J M.Initial experiments and analysis of blunt-edge vortex flows for VFE-2 configurations at NASA Langley, USA [J]. Aerospace Science and Technology, 2013, 24(1): 10-21.

        [5]Luckring J M, Hummel D J.What was learned from the new VFE-2 experiments [J]. Aerospace Science and Technology, 2013, 24(1): 77-88.

        [6]Hummel D J.The international vortex flow experiment(VFE-2): background, objectives and organization [J]. Aerospace Science and Technology, 2013, 24(1): 1-9.

        [7]Schutte A, Ludeke H.Numerical investigations on the VFE-2 65-degree rounded leading edge delta wing using the unstructured DLR TAU-Code [J]. Aerospace Science and Technology, 2013, 24(1): 56-65.

        [8]Fritz W.Numerical simulation of the peculiar subsonic flow-field about the VFE-2 delta wing with rounded leading edge[J]. Aerospace Science and Technology, 2013, 24(1): 45-55.

        [9]Crivellini A, D′Alessandro V, Bassi F.High-order discontinuous Galerkin RANS solution of the incompressible flow over a delta wing [J]. Computer & Fluids, 2013, 88(1): 663-677.

        [10]Chu J, Luckring J M. Experimental surface pressure data obtained on 65°delta wing across Reynolds number and Mach number ranges [R].NASA TM 4645, 1996.

        [11]Deng X, Zhang H. Developing high-order weighted compact nonlinear schemes[J]. Journal of Computational Physics, 2000, 165(1): 22-44.

        [12]Nonomura T, Fujii K. Effects of difference scheme type in high-order weighted compact nonlinear schemes[J]. Journal of Computational Physics, 2009, 228(10): 3533-3539.

        [13]Liu H, Deng X G, Mao M L. High-order behaviors of weighted compact fifth-order nonlinear schemes[J]. AIAA Journal, 2007, 48(8): 2093-2097.

        [14]Deng X G, Mao M L, Tu G H, et al. Extending weighted compact nonlinear schemes to complex grids with characteristic-based interface conditions [J]. AIAA Journal, 2010, 48(12): 2840-2851.

        [15]Menter F R, Kuntz M, Langtry R B. Ten years of industrial experience with the SST turbulence model[J]. Turbulence, Heat and Mass Transfer 4, 2003: 625-632.

        [16]Langtry R B, Menter F R. Correlation-based transition modeling for unstructured parallelized computational fluid dynamic codes[J]. AIAA Journal, 2009, 47(12): 2894-2906.

        [17]Deng X G, Min Y B, Mao M L, et al. Further study on geometric conservation law and application to high-order finite difference schemes with stationary grids [J]. Journal of Computational Physics, 2013, 239: 90-111.

        [18]Konrath R, Klein C, Engler R,et al.Analysis of PSP results obtained for the VFE-2 65°delta wing configuration at sub- and transonic speeds [C]//Proceedings of 44th AIAA Aerospace Sciences Meeting and Exhibit, 2006: 60.

        Numerical simulation on VFE-2 rounded leading edge delta wing using weighted compact nonlinear scheme

        WANG Guangxue1,2, WANG Shengye1, WANG Dongfang1, DENG Xiaogang1

        (1.College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China; 2.School of Physics, Sun Yat-sen University, Guangzhou 510006, China)

        In order to evaluate the influence of transition on the vortex structure of delta wing, a numerical simulation of VFE-2 rounded leading edge delta wing was carried out by using a high-order scheme-weighted compact nonlinear scheme and theγ-Reθtransition model. A comparisons between calculated results and experiment data indicate that the leading edge vortex begins at a certain distance of the wing apex and the transition has great influence on the onset of leading edge vortex. Using turbulence model without transition, the leading edge separation is delayed much, while with transition model the calculated results show a good agreement with experiment data. With transition model, numerical simulation on VFE-2 rounded leading edge delta wing at variation of the angle-of-attack was carried out. The calculated results which agree well with experiment data show that at a low angle-of-attack, there is no separation-induced leading-edge vortex, but with the increase of angle-of-attack the leading edge separation displays being closer to the trailing edge and moves upstream.

        transition model; delta wing; high-order scheme; weighted compact nonlinear scheme; separation-induced vortex

        10.11887/j.cn.201604002http://journal.nudt.edu.cn

        2016-03-01

        國防科學技術大學科研計劃資助項目(ZDYYJCYJ20140101)

        王光學(1976—),男,重慶忠縣人,博士研究生,E-mail:wgx111@sina.com;鄧小剛(通信作者),男,教授,博士,博士生導師,E-mail:xgdeng2000@vip.sina.com

        TN95

        A

        1001-2486(2016)04-008-06

        猜你喜歡
        三角翼迎角吸力
        深水大型吸力錨測試技術
        ROV在海上吸力樁安裝場景的應用及安裝精度和風險控制
        化工管理(2022年11期)2022-06-03 07:08:24
        三角翼機翼搖滾主動控制多學科耦合數(shù)值模擬
        航空學報(2021年12期)2022-01-10 07:56:14
        連續(xù)變迎角試驗數(shù)據(jù)自適應分段擬合濾波方法
        前緣和轉(zhuǎn)軸影響翼搖滾特性的數(shù)值模擬*
        深水吸力樁施工技術研究
        CY—06三角翼無人機
        航空模型(2016年10期)2017-05-09 06:22:13
        不同后掠角三角翼的靜態(tài)地面效應數(shù)值模擬
        超強吸力
        少年科學(2015年7期)2015-08-13 04:14:32
        失速保護系統(tǒng)迎角零向跳變研究
        科技傳播(2014年4期)2014-12-02 01:59:42
        亚洲av午夜福利精品一区二区 | 波多野无码AV中文专区| 日本精品人妻一区二区三区| 自拍偷区亚洲综合第一页| 蜜桃18禁成人午夜免费网站| av综合网男人的天堂| 亚洲国产激情一区二区三区| 乱人伦视频69| 亚洲中文字幕熟女五十| 久久午夜av一区二区| 久久久久亚洲精品无码系列| 国模丽丽啪啪一区二区| 美丽人妻被按摩中出中文字幕 | 亚洲h在线播放在线观看h| av大片在线无码免费| 久草视频在线这里只有精品| 99视频一区二区日本| 国产黑丝美腿在线观看| 人妻 色综合网站| 特级毛片a级毛片在线播放www| 国产AV无码无遮挡毛片| 国产av一区二区三区天美| 亚洲av无码国产精品色| 女人被狂躁高潮啊的视频在线看| 色综合久久无码中文字幕app| 亚洲免费成年女性毛视频| 国产一区二区三免费视频| 国产国语亲子伦亲子| 国产乱人伦在线播放| 韩国主播av福利一区二区| 精品日韩在线观看视频| 激情内射人妻1区2区3区| 一个人看的视频www免费| 久久精品国产亚洲一区二区| 国内偷拍第一视频第一视频区| 久久久精品人妻一区二区三区妖精| 亚洲成a∨人片在线观看不卡 | 国产无套粉嫩白浆在线观看| 欧美老熟妇欲乱高清视频| 日韩无码尤物视频| 婷婷开心五月亚洲综合|