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

帶襟翼的機翼尾渦合并數值計算

2017-04-19 09:31:38王志博
哈爾濱工業大學學報 2017年4期

王志博,孫 剛

(復旦大學 航空航天系,上海 200433)

帶襟翼的機翼尾渦合并數值計算

王志博,孫 剛

(復旦大學 航空航天系,上海 200433)

為獲得帶有襟翼的機翼尾渦的合并動力學過程,在驗證數值方法的基礎上,數值模擬帶有襟翼的機翼繞流尾流場,根據渦合并特征將合并過程劃分為誘導共轉階段、合并階段和軸對稱化階段.采用渦間距量綱—化尾流區域描述二渦誘導合并.變換弦長雷諾數、襟翼翼梢與機翼翼梢的間距、襟翼角度,改變襟翼翼梢渦與機翼翼梢渦的強度比,得到尾渦合并差異的特征.計算結果表明:隨著雷諾數的增加渦的強度增加,渦量的擴散程度減低,渦合并過程被推遲,空間誘導運動過程得到延長,渦系空間誘導運動增強,渦合并的雷諾數效應隨著雷諾數增加而減弱, 縮短渦間距,加速渦合并過程,但合并后的遠場渦尺寸和形態沒有顯著改變; 襟翼渦隨著襟翼角的減小而減弱,強度逐漸減弱的襟翼渦逐步被翼梢渦拉伸卷吸,微弱的襟翼渦系在近場中的合并過程完全改變,翼梢渦的運動軌跡并未受到誘導運動.翼梢渦合并的雷諾數效應表現為誘導運動過程的增強,尾流中強度小的渦系起不到明顯的誘導作用.

襟翼渦; 翼梢渦; 雷諾數; 渦合并; 環量

帶有翼面控制的航行體模型風洞試驗顯示,航行過程中機翼的翼梢渦和因襟翼打開導致翼面不連續而造成的渦系在近場的尾流中相互作用,這兩個渦系是尾流的主要成分.此外還包括從主機翼短艙結合部脫落的渦、尾翼翼梢渦、翼身組合體渦等各渦系,這些強度和尺寸各異的渦系在近場中最終合并為兩個對稱的集中渦.由于風洞中均勻來流的整流作用,將這些低動量高渦量的流體通過對流作用快速合并在一起.在近場尾流中,由于相互的誘導運動,脫體渦系發生了短波失穩包括已知的橢圓失穩、合作失穩和其他形式的失穩,渦在近場的互誘導作用使得渦核心區域渦量略有減低,渦發生了翻轉和擴散,渦系逐步趨向于合并.

相互誘導形成的渦系運動對尾流中最終形成的渦的尺寸和強度以及持續時間等均起著決定性作用,研究這些渦系的空間演化,可揭示尾跡的演化規律.通過文獻歸納,按照渦雷諾數的大小將渦的合并可分為準二維層流合并和三維短波橢圓合并兩種合并的形式.一般認為200

針對高雷諾數情況下的飛行器三維空間繞流,為了抽象出合理的尾渦結構的理論和實驗模型,Fabre等[3]推導兩個平行的有渦核結構的無限長渦絲的誘導失穩,并解釋了飛行器遠場尾跡中經常觀察的失穩形式.Crouch[4]建立了二渦對的共轉失穩模型,而Spalart[5]應用對轉二渦對描述飛機尾流的優化擾動,用最優擾動分析來描述實際觀察到的飛機尾渦構型是內襟翼的脫落渦和平尾翼梢渦的失穩可能性.Leweke等[6]研究了對轉渦對的三維短波失穩,失穩波長約為渦核的尺寸,通過細致的流場可視化和測速揭示了失穩渦對的空間結構.文獻[6]發現了各渦的擾動演化遵守一個運動匹配條件,每一個渦的耦合失穩被稱為合作失穩.Crow[7]構造了兩對對轉渦對空間演化,利用文獻[4]的對轉渦模型描述飛機遠場尾渦的失穩.

近場生成的渦系不滿足文獻[7]關于柱渦間距離遠大于渦的半徑的假設,也不滿足長波合作失穩的基本假設,也不是初始小擾動造成渦系失穩,近場渦系存在于機翼脫體的自由剪切流中,連通集中渦的剪切層起到了加速渦系合并的作用.航行器的尾流中存在大量的共轉渦合并的現象,尾渦合并不僅是遠場尾跡的形成的主要內容,而且與近場渦致阻力緊密相關,航行器布局設計應與渦致阻力之間通過渦合并動力學過程建立橋梁;遠場中存在的合作失穩則對飛機的尾流尺度效應產生影響,決定了機場的起降容積率.深入研究渦在近場的合并,將有助于揭示尾跡的形成機制和規模.Breitsamter[8]對飛機模型尾流風洞試驗進行了系統總結,指出飛機機翼的翼面不連續處的脫落尾渦在近場中以共轉和對轉合并是尾流結構形成的主要機制.Delbende等[9]擴展到三維螺旋形渦的合并,指出螺旋形渦的自誘導作用對渦的合并起到了延遲作用.Laurent等[10]指出共轉合并后的渦系不穩定性特征,Churchfield等[11]和Devenport等[12]對翼梢渦的合并進行了系統實驗和數值分析給出了翼尖渦初生的細節.Bristol等[13]對平行渦系誘導合并過程進行了模擬.

不同于經典的二維情況下的共轉渦合并,渦合并的三維效應顯示了獨特的三維渦拉伸,誘導共轉和剪切變形失穩等獨有的物理現象.調整舵面布局,例如調整內外襟翼的構型和布局,合理配比設計翼面積分布與后掠角、扭轉角等構型參數,制造合理的尾流渦系近場分布,延緩各類渦系共轉和對轉合并,對制造有益于航行性能的誘導下洗流場起到至關重要的作用,本文針對翼梢和襟翼脫落的渦系在尾流中的誘導共轉與合并運動特征,系統的描述尾流結構演化.

1 數值模型與驗證

1.1 計算設置

為研究繞流場中相鄰的兩個渦的相互誘導運動與合并,Bruin等[14-15]設計了如圖1所示的模型,并進行了系統的風洞試驗,測試了機翼的翼梢脫體渦與襟翼翼梢發放的渦在尾流中的相互作用的尾流速度場.試驗的弦長雷諾數為2.98×106,測量距機翼尾緣10倍的機翼展長至30倍的機翼展長范圍內若干橫截面的速度分布.可供開展數值模擬的驗證對比,試驗模型的詳細尺寸見表1.

圖1 實驗模型尺寸

表1 SWIM模型參數

計算采用有限體積法求解雷諾平均方程結合SSTk-ω湍流模型與近壁面處理模型,近壁面區域網格加密滿足給定來流速度條件下y+≈1,速度與壓力的耦合采用PISO算法,采用二階迎風對流格式離散動量方程,求解格式細節與離散數值處理方法在此不做詳述.比照風洞試驗,選取來流速度為60 m/s的常溫空氣作為計算介質,來流攻角為0°,遠場邊界采用速度入口,對應于風洞試驗中的來流湍流度為2.5%,遠場壓力出口和對稱面邊界,壁面采用了無滑移邊界條件.為獲得更為詳細的尾流結構,對尾流場的遠場采用200倍的機翼展長作為擴展計算域,并對尾流區域做了網格加密.

除了用于驗證風洞試驗的算例之外,計算中省略了風洞試驗中存在的支柱和導流罩,這些支撐結構對翼梢和襟翼后方的尾流流場造成的影響小,本文省略這些次要因素,網格無關性是通過渦核位置捕捉精度的比較來確定,最終選取網格為軸向×周向×展向 =450×100×50,網格劃分和邊界設置如圖2所示.

圖2網格劃分與邊界條件設置

1.2 數值計算結果驗證

本研究選取機翼展向距離導流罩0.15 m處的截面壓力測試數據與數值計算進行比較,壓力曲線分布比較如圖3所示.帶有襟翼的機翼增升效果顯著,機翼尾緣沒有顯著的壓力下降.

圖3 y/b1=0.5截面處壓力系數試驗值與計算值比較

Fig.3 Pressure coefficients comparison from numerical simulation and experiment aty/b1=0.5

比較圖4、5可知,采用PIV試驗獲得的尾渦位置與形態與數值計算結果基本一致,集中渦量的核心區位置略有差別,數值模擬襟翼渦與翼梢渦的形態和互誘導運動與實驗基本一致.本文建立的數值方法能夠準確模擬襟翼與翼梢形成的尾渦流動結構的演化.數值方法可獲得二渦在尾流中出現的合并現象,并描述近場渦系共轉渦對的合并造成渦的尺寸在尾流中迅速的擴大等三維誘導運動現象.

圖4 風洞試驗測得的x/b1=1.25處尾渦系位置與數值結果比較

圖5 風洞試驗測得x/b1=5.0處尾渦系位置與數值結果比較

2 尾流場結構特征

2.1翼梢渦與襟翼渦的誘導共轉階段x/Δ<20

在該階段中翼梢渦與襟翼渦各自具有獨立的渦系結構,可利用經典的蘭金渦理論模型描述渦核結構,隨著渦脫體后,受橫流效應和機翼脫體的自由剪切層影響,渦呈現軸對稱結構,渦的尺寸增長但渦量衰減.隨后兩個柱狀渦受到互相的誘導而翻轉,渦核軌跡形成了螺旋線,柱狀集中渦由于互誘導而逐漸失去軸對稱的特征,兩個渦逐漸靠近.

2.2 共轉渦的合并20

由于二渦間存在機翼尾流的自由剪切層,該剪切層起到了連接二渦的作用,剪切層卷起后,二渦的渦核靠近,渦首先出現橢圓失穩變形,渦相互靠近,直至不能區分兩個渦的邊界,無法辨識出渦核的各自位置.

2.3 結果渦的軸對稱化階段50

合并的渦在尾流中進一步軸對稱化,在軸對稱化的過程中,結果渦的尺寸保持一定.渦的結構逐步由橢圓形過渡到具有近似蘭金渦速度剖面的形狀.隨著尾渦在流場中的擴散,尾渦的尺寸達到極限尺寸,渦量的峰值緩慢衰減下降.

可見上述合并的動力學過程具有自相似的特征,即具有相似結構的襟翼渦和翼梢渦合并后的結果渦也具有與子渦相似的結構特征.尾流結構劃分如圖6所示.

圖6 帶有襟翼的機翼尾渦流場劃分

在第1個階段由于二渦的獨立性強,渦間相互作用以互誘導的翻轉運動為主,那么渦系可利用下述渦辨識方法得到渦核心運動軌跡.對于柱狀渦而言,沿著下游發展的渦系軸向渦量在任意橫截面的分布為

對應的環量為

單個渦的渦心位置為:

式中,積分上、下限以渦量的軸對稱分布區域而定.那么渦系的核心軌跡通過對截取有限個尾流中近場截面得到如圖7所示的互誘導翻轉運動.

圖7 尾流場中翼梢渦核與襟翼渦核的翻轉運動軌跡

采用求解雷諾平均方程的方法求解流場.雖然實驗中觀察到的渦在x/b1>10之后仍能存在高渦量的集中渦,但是由于過分數值耗散,本文中給出的量綱—化的渦量耗散較試驗顯著,圖8顯示了數值模擬給出的渦量在尾流中的耗散和擴散.在合并之前的階段,設定襟翼渦渦強度Γ1和機翼的渦強度Γ2,根據數值計算得到渦在近場運動的軌跡,忽略黏性耗散作用,根據渦管強度保持定理,在近場認為渦系的總的渦強度保持不變,那么有

圖8 襟翼渦與翼梢渦渦量衰減(U∞=60 m/s)

不考慮渦系的三維誘導作用造成的環量衰減,以平面渦系運動來近似模擬二渦的互誘導翻轉運動.渦系中的第m個渦受其他渦誘導產生的復速度為

兩邊同乘以Γm,兩邊對m求和得到:

積分得到:

那么復數的實部和虛部分別為常數,可寫出點渦系的重心位置為:

實際在尾流中不同的位置,處于機翼脫體的剪切層中的集中渦量迅速的耗散.故每一個流場剖面內的重心位置均不同,重心位置分布如圖9所示.重心位置的變化反映出二渦共轉誘導運動使渦處于不穩定狀態,相互靠近逐漸誘導擴散的特征.

3 渦間距對二渦合并影響

同一個雷諾數下,采用渦間距量綱—化尾流場截面位置,可完整的描述尾渦合并的過程,改變襟翼翼梢與機翼翼梢的間距Δ對二渦合并的影響,通過間距調整使二渦合并的過程得到改變,擴展計算了襟翼翼梢與機翼翼梢不同間距情況下的二渦合并的情況, 圖10顯示了采用x/Δ作為量綱—化參數描述不同渦間距的二渦合并的過程,上述劃分適用于描述尾渦合并的3個階段.如圖10所示,縮短襟翼翼梢與機翼翼梢的間距,加速翼梢渦與襟翼渦的合并過程,使得合并渦的尺寸得到快速的增長并迅速達到極限尺寸.

圖9 襟翼渦與翼梢渦共轉誘導運動的重心軌跡

如圖11在x/Δ=60處的尾流場截面處,渦已經處于合并后的軸對稱化階段,對于柱狀軸對稱渦而言,結果渦的尺寸并不隨著渦間距減小而顯著的增加,但渦量的峰值具有較大差異.表明量綱一的參數x/Δ無法描述軸對稱化后的渦量擴散的過程.

圖10 不同翼梢間距的近場尾渦合并

圖11 x/Δ=60處合并后渦的渦量剖面

4 渦合并的雷諾數效應

隨著來流速度的增加,尾流中集中渦的強度急劇的增加,合并過程必然隨著渦強度增加而產生差別.對于全尺度的實際流動而言,渦的耗散過程是極為緩慢的,但是對于模型尺度的渦而言,渦的耗散相對較快.本文改變翼的展長對比發現渦的合并依賴于渦間距離和渦的初始尺寸與形態.對于較低的雷諾數而言,集中渦總是尺寸快速增長,因此渦總會在近場首先合并.隨著雷諾數增加,渦的合并過程被推遲到更大的x/Δ值,因此近場的渦互誘導運動現象與低雷諾數略有不同,即使弦長雷諾數大于105后,尾流中渦的合并過程繼續隨著雷諾數的增長而被推遲.本文選取對應于風洞試驗的標準算例對應風速的0.5倍和2.0倍進行模擬,進而分析渦系運動軌跡與合并雷諾數效應.圖12、13 顯示在同一個截面位置不同量級的雷諾數情況下兩個渦的合并狀態.對比圖12、13可知隨著渦的強度增加渦量的擴散程度減低,合并被推遲,獨立渦量的空間互誘導運動得到延長.

圖12 Δ=75 mm,x/Δ=20 處軸向渦量等值線云圖

圖13 Δ=25 mm,x/Δ=20 處軸向渦量等值線云圖

在與風洞實驗對應的數值算例中發現,由于襟翼渦的渦量集中,初始強度約是翼梢渦的2倍,那么翼梢渦受到襟翼渦誘導而發生翻轉卷曲運動,如圖14所示隨著雷諾數的增加而逐漸增強,渦的卷曲效應增加,由圖14可知翼梢渦翻轉軌跡曲率隨雷諾數而增大.在圖14中顯示二渦強度比在流向趨于接近的同時,二渦的誘導作用形成的軌跡趨于一致.

圖14 襟翼渦與翼梢渦近場軌跡(x/Δ<20)

5 襟翼渦強度減弱時的渦合并

通過改變襟翼角度,使得二渦的強度比在初始生成時減小,分析二渦強度變化對合并過程的影響.渦間誘導效應決定于渦強度的大小和二者強度之比Γ1/Γ2.在上述驗證模型算例的基礎上調整襟翼角度分別為20°、10°、5°,隨著襟翼角的增加,襟翼渦的強度不斷的增大,襟翼渦對翼梢渦的誘導作用變得越來越顯著,結果渦的尺寸有所增長.襟翼角為5°時,處于自由剪切尾流中的襟翼渦受到機翼的翼梢渦的拉伸而變形,形成快速的橢圓失穩,逐漸拉伸成為剪切層被翼梢渦卷起,最終形成自由剪切層與機翼翼梢渦合并,翼梢渦并不發生翻轉運動.隨著襟翼角度的增加,至10°時,襟翼渦的強度增大,襟翼渦與機翼翼梢渦的互誘導運動增強,翼梢渦出現了翻轉運動,襟翼渦與翼梢渦合并后的渦尺寸顯著的增大.在如圖15所示的尾流云圖可觀察到襟翼20°時二渦呈現了互誘導翻轉運動.

圖15 尾流中渦量分布云圖對比(Δ=100 mm,截面位置x/Δ=2.5至15間隔2.5)

6 結 論

1)采用渦間距值進行量綱—化,可劃分含襟翼渦的翼梢尾流結構,根據運動特征可劃分為誘導共轉階段、共轉渦的合并、結果渦的軸對稱化等3個主要階段.但在對襟翼渦與翼梢渦的合并過程中未觀察到二渦互誘導翻轉運動軌跡超過180°的情況,說明襟翼渦系與翼梢渦系的三維空間拉伸失穩存在極限.并且在同一個襟翼角和雷諾數下,雖然襟翼渦和翼梢渦的間距不同,但合并之后的結果渦具有相近的尺寸和形態.

2)襟翼渦的強度遠小于翼梢渦的強度時,在近場襟翼渦被拉伸變形并被翼梢渦卷起,而翼梢渦的運動軌跡并未受到襟翼渦的誘導運動,隨著襟翼角度的增加,襟翼渦強度增大,二渦形成了相互誘導運動與共轉合并.

3)當Re<106時,隨著Re的增加,可觀察到尾流中一定程度的Re效應.隨著Re的增大,渦的強度增加,渦量的擴散程度減低,合并過程被推遲,空間誘導運動得到延長,Re效應減低.

[1] PAULO J, FERREIRA S, PEREIRA J. Reynolds number dependence of two-dimensional laminar co-rotating vortex merging[J]. Theoretical and Computational Fluid Dynamics, 2005, 19(1):65-75. DOI: 10.1007/s00162-004-0154-0.

[2] CERRETELLI C, WILLIAMSON C. The physical mechanism for vortex merging [J]. Journal of Fluid Mechanics, 2003, 475(1):41-77. DOI: 10.1017/S0022112002002847.

[3] FABRE D, JACQUIN L, LOOF A. Optimal perturbations in a four-vortex aircraft wake in counter-rotating configuration [J]. Journal of Fluid Mechanics, 2002, 451(1):319-328.

[4] CROUCH J D. Instability and transient growth for two trailing-vortex pairs [J]. Journal of Fluid Mechanics, 1997, 350(1):311-330.

[5] SPALART R. Airplane trailing vortices[J]. Annual Review of Fluid Mechanics,1998,30:107-138.

[6] LEWEKE T, WILLIAMSON C H K. Cooperative elliptic instability of a vortex pair [J]. Journal of Fluid Mechanics, 1998, 360(1): 85-119.

[7] CROW S. Stability theory for a pair of trailing vortices[J]. AIAA Journal, 1970, 8(12):2172-2179. DOI: 10.2514/3.6083.

[8] BREITSAMTER C. Wake vortex characteristics of transport aircraft [J]. Progress in Aerospace Sciences, 2011, 47(2):89-134.DOI: 10.1016/j.paerosci.2010.09.002.

[9] DELBENDE I, PITON B, ROSSI M. Merging of two helical vortices [J]. European Journal of Mechanics-B/Fluids, 2015, 49(2):363-372.DOI: 10.1016/j.euromechflu.2014.04.005.

[10]LAURENT J, DAVID F, DENIS S, et al. Unsteadiness, instability and turbulence in trailing vortices[J]. Comptes Rendus Physique, 2005, 6(5):399-414. DOI:10.1016/j.crhy.2005.05.007.

[11]CHURCHFIELD J, BLAISDELL A. Numerical simulations of a wingtip vortex in the near field[J]. Journal of Aircraft, 2012, 46(1):230-243.DOI:10.2514/1.38086.

[12]DEVENPORT J, RIFE C, LIAPIS I. The structure and development of a wing-tip vortex[J]. Journal of Fluid Mechanics, 1996, 312(4):67-106. DOI: 10.1017/S0022112096001929.

[13]BRISTOL R L, ORTEGA J M, MARCUS P S. On cooperative instabilities of parallel vortex pairs[J]. Journal of Fluid Mechanics, 2004, 517(517):331-358. DOI:10.1017/S0022112004001016.

[14]BRUIN A. Test report for wake survey tests behind SWIM model geometry in DNWLST and DNW-LLF wind tunnels: NLR-TR-2001-183 [R]. Netherlands: DNW-LLF/ LST wind tunnels, 1999.

[15]BRUIN A,Ganzevles F. Data analysis of wake survey tests behind SWIM model geometry in DNW-LST and DNW-LLF wind tunnels: NLR-TR-2001-201 [R]. Netherlands: DNW-LLF/ LST wind tunnels, 2001.

(編輯 張 紅)

Numerical simulation of trailing vortex merge in a flapped wing wake

WANG Zhibo, SUN Gang

(Department of Aeronautics and Astronautics, Fudan University, Shanghai 200433, China)

A validated numerical simulation modeling of flapped wing wake is carried out by using RANS discretion. The flapped wing wake region is divided into three regions with significant different dynamic features. Three regions including co-rotation field, merging field and axial symmetry far field are induced. A non dimensional parameter of vortex interval is used to describe above regions. The Reynolds number, vortex interval, flap angle are used as control variables to change ratio of circulation of flap vortex and tip vortex. Reynolds effect is found as prolonged co-rotation field and enhanced rotation motion in wake. Decreased vortices interval leads an accelerated merging in near field but a limited dimension of resulted vortex. Small flap angle leads a weak flap vortex. The dominated tip vortex stretches and filament flap vortex in free share layer. The merging process is totally changed. The tip vortex trajectory does not expand. Reynolds number effect is the enhancement of induction in tip vortices merge. Induction of vortices with low intensity is not obvious in the wake.

flap vortex; wingtip vortex; Reynolds number; vortex merge; circulation

10.11918/j.issn.0367-6234.201510028

2015-10-12

國家重點基礎研究發展計劃(2014CB744802)

王志博(1983—),男,博士研究生; 孫 剛(1966—),男,教授,博士生導師

孫 剛,Gang_sun@fudan.edu.cn

O355

A

0367-6234(2017)04-0073-07

主站蜘蛛池模板: 中文字幕日韩视频欧美一区| 蜜桃视频一区二区| 亚洲性日韩精品一区二区| 漂亮人妻被中出中文字幕久久| 国产黑丝视频在线观看| 全部毛片免费看| 三上悠亚一区二区| 国产激情第一页| 精品一区国产精品| 国产成人精品综合| 国产综合色在线视频播放线视| 亚洲va视频| 亚洲国产欧美国产综合久久| 国产精品视频观看裸模| 国产chinese男男gay视频网| 色悠久久综合| 女人一级毛片| a级免费视频| 国产xx在线观看| 亚洲愉拍一区二区精品| 第一页亚洲| 国产美女无遮挡免费视频网站| 亚洲高清在线天堂精品| 久久人妻系列无码一区| 午夜性爽视频男人的天堂| 丝袜高跟美脚国产1区| 成年人国产网站| 国产第一页屁屁影院| 毛片免费试看| 亚洲综合18p| 久久国产精品夜色| 青青操国产视频| 亚洲日韩AV无码精品| 欧美日韩高清| 99r在线精品视频在线播放| 九九视频在线免费观看| 亚洲一区精品视频在线 | 亚洲天堂免费观看| 538国产在线| 四虎成人精品在永久免费| 香蕉精品在线| 伊在人亚洲香蕉精品播放| 欧美成人免费午夜全| 亚洲一级毛片在线观| 亚洲日本一本dvd高清| 青青操视频在线| 国产美女在线观看| 亚洲成人高清无码| 欧美日一级片| 成人综合在线观看| 亚洲香蕉伊综合在人在线| 99精品在线视频观看| 国产最爽的乱婬视频国语对白| 一级毛片无毒不卡直接观看| 成人毛片免费在线观看| 亚洲大学生视频在线播放| 亚洲av无码牛牛影视在线二区| 91在线精品麻豆欧美在线| 亚洲精品第1页| 美女无遮挡免费网站| 亚洲中文字幕在线观看| 国产精品自在自线免费观看| 小蝌蚪亚洲精品国产| 国产女主播一区| 国产成人av大片在线播放| 伊人久久精品无码麻豆精品 | 嫩草国产在线| 久久大香香蕉国产免费网站| 无码网站免费观看| 亚洲中文字幕97久久精品少妇| 国产在线观看91精品亚瑟| 久久香蕉国产线看观看式| 欧美在线中文字幕| 亚洲第一视频网| 久久福利网| 国产精品手机在线播放| 亚洲日产2021三区在线| 亚洲欧美在线看片AI| 2020国产免费久久精品99| 国产在线高清一级毛片| 丝袜美女被出水视频一区| 国产亚洲精|