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

引入DMD方法研究有/無控氣流分離的動態(tài)結構

2017-11-20 03:12:19洪樹立黃國平
航空學報 2017年8期
關鍵詞:模態(tài)結構分析

洪樹立, 黃國平

1.南京航空航天大學 能源與動力學院, 南京 210016 2.南京航空航天大學 江蘇省航空動力系統重點實驗室, 南京 210016

引入DMD方法研究有/無控氣流分離的動態(tài)結構

洪樹立1,2,*, 黃國平1,2

1.南京航空航天大學 能源與動力學院, 南京 210016 2.南京航空航天大學 江蘇省航空動力系統重點實驗室, 南京 210016

為分析非定常流動控制技術抑制分離流的機理,對彎曲擴壓通道的試驗模型進行了數值模擬,針對擴壓通道在無控和采用最佳射流頻率狀態(tài)下的計算結果引入了動力模態(tài)分解(DMD)技術進行分析。通過DMD技術能夠將包含時空信息的擴壓通道復雜流場進行分解,捕獲流場包含的動力信息和對應的擬序流動結構。將無控和有控流場分解的結果進行對比分析后表明:采用有效激勵措施時,和脫落渦頻率一致的渦系對流場的影響更加突顯,流場整體上表現得更加有序;非定常控制抑制了一部分渦的增長,使得各模態(tài)整體上更加穩(wěn)定;而有控流場占主導地位的渦系結構相比無控流場較為有序,且對主流區(qū)未形成明顯的直接影響。

流動分離; 非定常; 流動控制; 擬序結構; 動力模態(tài)分解(DMD)

航空動力系統壓縮部件通道內存在著較大的逆壓梯度,容易產生流動分離并形成大尺度的湍動渦,尤其在航空動力系統高負荷的發(fā)展下,使得流動更加惡化,影響了其效率和穩(wěn)定工作范圍。為了改善該現象,越來越多的研究者開始關注于采用流動控制的方法[1]。相比于定常流動控制,非定常流動控制往往用較小的能量注入起到相同的效果。盡管非定常流動控制取得了許多有效的成果并獲得了一些控制規(guī)律[2-4],然而對于其內在機理的認識還有欠缺。尤其是像航空動力系統的壓縮部件,分離流、湍動渦等常跨越了大尺度和頻率范圍[5-6],這種內部的復雜流動增加了對其抑制的難度,因此有必要進一步加強對非定常流動控制機理的研究。而對于各種渦系互相耦合的復雜流動,為研究分析帶來了困難。研究者發(fā)現,這種看似復雜無序的流動在時空的演化過程中具有某種有序的規(guī)律,存在著一些能夠反映流場內在規(guī)律的本質特征[7-8]。

為了得到這種流動的擬序結構,可以將時空流場包含的信息進行解耦。近年來研究較多的是本征正交分解(Proper Orthogonal Decomposition, POD),其基本思想是尋找一組基元的抓拍“Snapshot”流場,可被重構為這組正交基與各自的時間權重系數(隨時間分布的權重系數)的乘積之和。通過POD分解,能夠將原時空流場分解為具有一定能量等級的擬序結構。此方法常被用于各類的流場分析,近年來對于非定常流動控制機理的研究也有較好的應用[9-10]。而POD在提取流場時往往通過某種統計平均得到,丟失了系統的相位信息,具體反映在POD分解的某個模態(tài)依舊參雜了不同頻率的渦系結構,因而難以在動力學層面解析原流場。這使之在應用中也遇到了一些問題,例如不能很好地捕獲能量不高的小擾動[11];高階模態(tài)下摻雜的低頻信息被忽略[12]。

與POD不同的是,動力模態(tài)分解(Dynamic Mode Decomposition, DMD)是從動力層面對流場進行分解和提取,分解后得到的各階模態(tài)在時間上是互不相關的。DMD方法是由美國普林斯頓大學的Rowley等[13]和巴黎綜合理工大學的Schmid[14]首次提出的,其思路來源于非線性系統的Koopman分析。Koopman分析的思路是將對原系統的研究轉移到對Koopman算子的研究上,而該算子則包含了原動力系統的所有信息[13]。DMD方法就是Koopman分析的一種近似處理,這種技術提供了一種新的方法和視角去描述流動的擬序結構。

近年來DMD方法已經成功運用于一些特征較為明顯的流動結構。Schmid[15]運用DMD方法提取了射流中與頻率相關的占主導地位的動力學模態(tài);Zhang等[16]將POD和DMD方法運用于圓柱繞流,指出對于流動特征單一的單圓柱繞流,盡管POD方法也捕獲了卡門渦街,但是其高階模態(tài)的物理意義難以解釋,而DMD分解得到的高階模態(tài)就可以認為是卡門渦街的“諧波”;而對于雙圓柱繞流,其具有2個重要渦街,POD捕獲的渦結構摻雜了其他渦系,而DMD清晰地捕獲了2個渦系;Seena和Sung[17]研究了空腔流,通過DMD分析發(fā)現當來流附面層達到一定條件時,會引起共振現象(Hydrodynamic Resonances),使得空腔內附面層的結構和來流附面層的結構在頻率和波數上都一致;Nastase等[18]研究了三維的葉狀射流,通過DMD分析發(fā)現Kelvin-Helmholtz結構是其主要形態(tài)。目前來說,DMD方法主要應用于分析典型的周期性非定常流動,而對較為復雜的非定常擬序結構流動的應用還很少見到報道。

本文采用DMD方法,對擴壓通道內流場在無控和采用非定常流動控制下的狀態(tài)進行分析研究,提取擴壓通道內受控前后流場的擬序結構,并有效使用DMD分析結果探討了非定常激勵對擬序分離流的控制機理。

1 DMD基本原理與實現方法

采用DMD方法可以對海量的試驗數據或數值模擬結果進行處理分析,其算法推導過程參見文獻[13-14],本文主要介紹其基本原理和實現方法。

1.1 時空速度場的正切近似

vj+1=Avj

(1)

式中:A為系統矩陣,反映了前后2個時刻流場之間的關系。對于非線性系統,上述假設是正切近似的,而對于線性系統則嚴格成立,因此時空流場可以表述為

(2)

1.2 系統矩陣的物理意義

實際上式(2)的形式和狀態(tài)空間向量矩陣的表達式形式一樣,只是將對時間的導數換成了對時間的推進,這相當于將多自由度系統的自由振動分析應用在流場的動力分析上[19]。

圖1 二自由度無阻尼系統[19]Fig.1 No-damping system of two-degree-of-freedom[19]

對于流場,式(2)中A就包含了原時空流場的動力信息,其特征值反映了各流場空間點的速度增益(或阻尼)和比率變化的頻率信息;A的特征向量反映了相應向量下的流場信息,相當于多自由度系統的固有振型。

1.3 系統矩陣降維處理與分解

一般來說A的維度較高(依賴于流場數據點),因而很難直接進行求解,而應用Arnoldi方法[20],可以將A進行降維處理,把求解A的特征值和特征向量轉換到求解A的伴隨矩陣S上,S的表達式為

(3)

伴隨矩陣S只有最后一列是未知的,該列與前N-1個流場乘積的線性組合就反映了后一時刻的流場,可以得到

(4)

S=T-1ΛT

(5)

定義矩陣Φ為

(6)

則式(4)可以變?yōu)?/p>

(7)

(8)

圖2 DMD具體過程[21]Fig.2 Detailed process of DMD[21]

DMD技術將時空流場結構做線化近似,并進行降維處理,通過提取流場的主要動力信息,為研究復雜流場和不穩(wěn)定性觸發(fā)的擬序流動結構提供了認識相關機理的工具。

2 典型通道內氣流分離及其控制

2.1 數值模擬方法

本文引入DMD方法,嘗試分析典型通道內氣流分離的動態(tài)擬序結構,及非定常激勵對其流場結構的作用。為此,這里以一個前期有較多研究基礎的拐彎擴壓通道作為流場研究對象。在擴壓通道彎角處設置了葉片,并針對葉背處出現的二次流采用無源脈沖射流控制技術進行抑制。利用葉盆(對應于大氣高壓一側)和葉背處(對應于通道內低壓一側)的壓差驅動流體經過葉片縫隙,并通過縫隙中轉子的運動形成非定常有質量脈沖射流,最大射流速度為35 m/s,具體結構如圖3所示。

采用大渦模擬進行數值仿真,邊界條件根據試驗給定,進口為大氣條件,出口設定背壓,保證進口馬赫數Ma=1,采用雙時間步長進行加速求解,時間步長為1×10-5s,計算域及坐標設置如圖4所示。通過對比不同狀態(tài)下脫落渦頻率、相對總壓損失系數、平均壓力等參數表明試驗和數值模擬結果具有較好的吻合度,因此本文采用的數值模擬具有一定的可信度,更多相關的試驗和數值模擬介紹可以參考文獻[22-23]。

圖3 擴壓通道結構圖Fig.3 Structure diagram of divergent channel

圖4 擴壓通道計算域Fig.4 Computational domain of divergent channel

2.2 整體控制效果

圖5 控制效果對比Fig.5 Contrast of control effect

3 基于DMD的彎曲通道流場分離及其控制分析

基于DMD原理,對z方向(垂直于xOy平面)上的渦量Ωz進行DMD分解,Ωz的表達式為

(9)

式中:uy和ux分別為y和x方向上的速度分量。有控流場和無控流場均選取相同的時間步長Δt和時間步數N,時間步長Δt取為1×10-4s,時間步數N取為300,共記0.03 s,約為7個分離渦準周期。

3.1 全局能量和特征值標準偏差分析方法

特征值標準偏差的定義為

(10)

3.2 非定常控制對全局能量分布的影響

圖6給出了有控和無控流場全局能量和頻率的對應關系(頻率為0的模態(tài)可以看作對時間平均的流場,其全局能量相比其他各階高出一個量級,圖中并未給出)。從圖中可以看出,隨著頻率的增加,全局能量大致呈遞減趨勢。

無控流場的全局能量峰值所在頻率為373 Hz,有控流場全局能量峰值位置為331 Hz,相比無控流場其峰值點的頻率更加接近脫落渦頻率,說明盡管無控流場中脫落渦占了主導地位,但是摻雜了其他一些對流場有重要影響的渦系結構。此外有控流場全局能量最大值相比其他頻率以及無控流場各頻率對應的全局能量要高出一個數量級,表明施加非定常控制后,通過不同模態(tài)間的能量轉移最終使得流場中占主導地位的渦系結構更加突出,表明流場更加有序。

圖6 有/無控流場的全局能量分布 Fig.6 Distributions of global energy norm of flowfield with and without control

3.3 非定常控制對整體穩(wěn)定性的影響

圖7給出了流場施加非定常控制前后特征值的分布(圖中λRe和λIm分別表示特征值的實部和虛部),無控流場DMD分解得到的特征值有46個模態(tài)落在單位圓外,103個落在單位圓內,受控后有39個模態(tài)落在單位圓外,126個落在單位圓內,表明非定常控制抑制了一部分渦的增長,呈衰減趨勢的渦系增多。

無控流場特征值標準偏差σ=0.059,有控流場特征值標準偏差σ=0.023,說明在施加非定常控制之后,各模態(tài)在整體上更加接近穩(wěn)定的固定周期運動。

圖7 有/無控流場的特征值分布 Fig.7 Distributions of flow field eigenvalues with and without control

總的來說,通過DMD技術,可以讓更多動力模態(tài)處于單位圓內,并使得特征值標準偏差減少,這都說明經恰當周期激勵控制(F=1.0),流場能變得更為有序。

3.4 非定常控制對擬序結構的影響

采用DMD方法,可以把流場各渦系的動力學行為和相應的空間擬序結構相結合,以便更好地分析和描述原流場的特征。圖8為彎曲擴壓通道的分離流場和施加恰當周期激勵控制(F=1.0)后流場的不同模態(tài)的渦量云圖。對于無控流場,將頻率為373、651、2 522 Hz對應的模態(tài)分別記為MODE A、MODE B和MODE C。對于有控流場,將采用DMD方法后各相近頻率為331、643、2 517 Hz 對應的模態(tài)分別記為MODE A′、MODE B′和MODE C′。

從DMD的結果可以看出,渦的空間尺度隨著頻率的增加有所減小。并且采用DMD方法,能夠捕獲有控和無控流場的主導渦系結構,即MODE A′和MODE A。對比圖8(a)和8(b)中有控和無控的分解結果可以看出,MODE A′反映的空間渦系結構較MODE A更加均勻有序,并且MODE A中的渦結構較為混亂,在向通道下游演變的過程中影響到了主流區(qū),這一點也反映在MODE B′和MODE B中。

采用DMD方法,能夠較為明顯地反映出施加恰當周期激勵控制(F=1.0)后流場的變化,并且可顯示出各主要動力模態(tài)的流動結構及控制前后差異。因此,DMD方法可以為分析復雜擬序流動提供有效的技術途徑。

圖8 動力模態(tài)渦量分布Fig.8 Vorticity distributions of dynamic modes

4 結 論

本文以彎曲擴壓通道內分離流為研究對象,采用大渦模擬方法對脈沖射流抑制擴壓通道內流動分離及施加周期性激勵控制后的流場進行了數值仿真,并引入DMD動力模態(tài)分解法對計算的結果進行分析,得出如下結論:

1) 采用DMD分析方法可以有效捕獲流場的動力信息并且顯示其對應模態(tài)的空間渦系結構,分解所得的各階模態(tài)隨著頻率的增加對應的空間尺度有所減小。

2) 通過DMD分解的頻譜表明,相比于無控流場,有控流場分解所得與脫落渦頻率一致的渦系結構所占的主導地位更加突顯,反映出流場更加有序。

3) 非定常控制抑制了一部分渦的增長,呈衰減趨勢的渦系增多,各模態(tài)在整體上更加穩(wěn)定。

4) 無控流場占主導地位的渦系結構較為混亂,該渦在向下游演變的過程中會更多地影響到主流區(qū),而有控流場占主導地位的渦系結構比較有序,且沿通道向下游發(fā)展影響到的主流相對減少。

[1] RIVIR R B, SONDERGAARD R, BONS J P, et al. Passive and active control of separation in gas turbines: AIAA-2000-2235[R]. Reston, VA: AIAA, 2000.

[2] HERGT A, MEYER R, ENGEL K. Experimental investigation of flow control in compressor cascades: GT-2006-90415[R]. New York: ASME, 2006.

[3] GREENBLATT D, WYGNANSKI I J. The control of flow separation by periodic excitation[J]. Progress in Aerospace Sciences, 2000, 36(7): 487-545.

[4] WYGNANSKI I J. The variables affecting the control of separation by periodic excitation: AIAA-2004-2505[R]. Reston, VA: AIAA, 2004.

[5] SEIICHI I, MASATO F, KENICHIRO I, et al. Vortical flow structure and loss generation process in a transonic centrifugal compressor impeller: GT-2007-27791[R]. New York: ASME, 2007.

[6] HASHMI S, QIAO W Y, CHEN P P. Prediction of hub corner stall characteristics of a highly loaded low speed single stage fan[J]. Journal of Thermal Science, 2011, 20(2): 106-114.

[7] ADRIAN R J, MEINHART C D, TOMKINS C D. Vortex organization in the outer region of turbulent boundary layer[J]. Journal of Fluid Mechanics, 2000, 422: 1-54.

[8] DENNIS D J C, NICKELS T B. Experimental measurement of large-scale three-dimensional structures in a turbulent boundary layer. Part 1: Vortex packets[J]. Journal of Fluid Mechanics, 2011, 673: 180-217.

[9] FENG L H, WANG J J, PAN C. Proper orthogonal decomposition analysis of vortex dynamics of a circular cylinder under synthetic jet control[J]. Physics of Fluid, 2011, 23: 014106.

[10] 朱劍鋒, 黃國平, 傅鑫, 等. 基于POD方法的彎曲擴壓通道分離流控制的時空特性分析[J]. 航空學報, 2014, 35(4): 921-932.

ZHU J F, HUANG G P, FU X, et al. Spatiotemporal characteristics analysis for controlling flow separation in divergent curved channels with POD method[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(4): 921-932 (in Chinese).

[12] ZHANG W, WANG Y, QIAN Y H. Stability analysis for flow past a cylinder via lattice Boltzmann method and dynamic mode decomposition[J]. Chinese Physics B, 2015, 24(6): 378-384.

[14] SCHMID P J. Dynamic mode decomposition of numerical and experimental data[J]. Journal of Fluid Mechanics, 2010, 656: 5-28.

[15] SCHMID P J. Application of the dynamic mode decomposition to experimental[J]. Experiments in Fluids, 2011, 50(4): 1123-1130.

[16] ZHANG Q S, LIU Y Z, WANG S F. The identification of coherent structures using proper orthogonal decomposition and dynamic mode decomposition[J]. Journal of Fluids and Structures, 2014, 49: 53-72.

[17] SEENA A, SUNG H J. Dynamic mode decomposition of turbulent cavity flows for self-sustained oscillations[J]. International Journal of Heat & Fluid Flow, 2011, 32(6): 1098-1110.

[18] NASTASE I, MEALEM A, EI HASSAN M. Image processing analysis of vortex dynamics of lobed jets from three-dimensional diffusers[J]. Fluid Dynamics Research, 2011, 43: 065502.

[19] 胡海巖. 機械振動基礎[M]. 哈爾濱: 哈爾濱工業(yè)大學出版社, 2004.

HU H Y. Mechanical vibration foundation[M]. Harbin: Harbin Institute of Technology Press, 2004 (in Chinese).

[20] GREENBAUM A. Iterative methods for solving linear system[M]. Philadelphia: SIAM, 1997.

[21] SCHMID P J, LI L, JUNIPER M P, et al. Applications of the dynamic mode decomposition[J]. Theoretical and Computational Fluid Dynamics, 2011, 25(1): 249-259.

[22] 朱劍鋒, 黃國平, 傅鑫, 等. 一種控制氣流分離的無源微脈沖射流技術研究[J]. 航空學報, 2013, 34(8): 1757-1767.

ZHU J F, HUANG G P, FU X, et al. Investigation of technology for controlling flow separation by micro-pulsed-jet without external device[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(8): 1757-1767 (in Chinese).

[23] 朱劍鋒, 黃國平, 傅鑫, 等. 脈沖射流控制彎曲擴壓管道流動分離的特點[J]. 航空動力學報, 2015, 30(12): 2942-2948.

ZHU J F, HUANG G P, FU X, et al. Characteristic of controlling flow separation in divergent curved channels by pulsed jet[J]. Journal of Aerospace Power, 2015, 30(12): 2942-2948 (in Chinese).

[24] 潘翀, 陳皇, 王晉軍. 復雜流場的動力學模態(tài)分解[C]//第八屆全國實驗流體力學學術會議. 北京: 中國學術期刊電子出版社, 2010: 77-82.

PAN C, CHEN H, WANG J J. Dynamic modal decomposition of complex flow field[C]//The Eighth National Conference on Experimental Fluid Mechanics. Beijing: China Academic Journal Electronic Publishing House, 2010: 77-82 (in Chinese).

(責任編輯: 王嬌)

*Correspondingauthor.E-mail:hong_0815@163.com

IntroducingDMDmethodtostudydynamicstructuresofflowseparationwithandwithoutcontrol

HONGShuli1,2,*,HUANGGuoping1,2

1.CollegeofEnergyandPowerEngineering,NanjingUniversityofAeronauticsandAstronautics,Nanjing210016,China2.JiangsuProvinceKeyLaboratoryofAerospacePowerSystem,NanjingUniversityofAeronauticsandAstronautics,Naning210016,China

Toanalyzethemechanismofdepressingflowseparationwithunsteadyflowcontroltechnology,anumericalsimulationfortheexperimentalmodelofthedivergentchanneliscarriedout.Dynamicmodedecomposition(DMD)technologyisadoptedtostudytheflowfieldofthecurveddivergentchannelwithandwithoutpulsedjetcontrol.WithDMDtechnology,thecomplexflowfieldofthedivergentchannelcontainingspatialandtemporalinformationcanbedecomposedhierarchically,anddynamicalinformationaswellasspatialcoherentstructurecorrespondingtothevortexcanbecapturedandrendered.Acomparisonofthedecomposedflowfieldswithandwithoutcontrolshowsthatwitheffectiveexcitation,thecoherentstructurewiththefrequencyapproximatingthefrequencyofthesheddingvortexbecomesmoredominantintheinitialflowfield,andtheoverallflowfieldturnsoutbetomoreordered.Somecoherentstructures,decomposedfromflowfieldwithoutcontrol,aresuppressedtomakeallmodesmoresteady.Thedominantstructureofthecontrolledflowfieldhasnoobviousinfluenceonthemainflow.

flowseparation;unsteady;flowcontrol;coherentstructure;dynamicmodedecomposition(DMD)

2016-10-21;Revised2016-11-18;Accepted2016-12-15;Publishedonline2017-01-121117

URL:www.cnki.net/kcms/detail/11.1929.V.20170112.1117.003.html

NationalNaturalScienceFoundationofChina(51176072)

2016-10-21;退修日期2016-11-18;錄用日期2016-12-15; < class="emphasis_bold">網絡出版時間

時間:2017-01-121117

www.cnki.net/kcms/detail/11.1929.V.20170112.1117.003.html

國家自然科學基金 (51176072)

.E-mailhong_0815@163.com

洪樹立, 黃國平. 引入DMD方法研究有/無控氣流分離的動態(tài)結構J. 航空學報,2017,38(8):120876.HONGSL,HUANGGP.IntroducingDMDmethodtostudydynamicstructuresofflowseparationwithandwithoutcontrolJ.ActaAeronauticaetAstronauticaSinica,2017,38(8):120876.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.120876

V231.3

A

1000-6893(2017)08-120876-08

猜你喜歡
模態(tài)結構分析
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
隱蔽失效適航要求符合性驗證分析
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
電力系統不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統及其自動化發(fā)展趨勢分析
論《日出》的結構
國內多模態(tài)教學研究回顧與展望
創(chuàng)新治理結構促進中小企業(yè)持續(xù)成長
基于HHT和Prony算法的電力系統低頻振蕩模態(tài)識別
由單個模態(tài)構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 91欧美在线| 国产亚洲精品97在线观看| 亚洲天堂视频网| 久久窝窝国产精品午夜看片| 小说 亚洲 无码 精品| 97精品国产高清久久久久蜜芽| 日韩美毛片| 日韩国产黄色网站| 日韩专区欧美| 72种姿势欧美久久久久大黄蕉| 女人毛片a级大学毛片免费| 国产精品密蕾丝视频| 欧美高清国产| 日韩欧美91| 欧美一级专区免费大片| 亚洲国产在一区二区三区| 99ri国产在线| 精品福利视频导航| 99中文字幕亚洲一区二区| 伊人91在线| 午夜啪啪福利| 午夜综合网| 视频国产精品丝袜第一页| 亚洲无码精品在线播放| 免费在线播放毛片| 欧美午夜视频| 永久天堂网Av| 亚洲一区二区三区国产精华液| 国产微拍一区| 久久综合色视频| 久久久四虎成人永久免费网站| 日韩欧美国产综合| 男人天堂亚洲天堂| 91精品啪在线观看国产60岁 | 亚洲第一香蕉视频| 三上悠亚精品二区在线观看| 亚洲精品麻豆| 亚洲第一视频网| 九色综合伊人久久富二代| 欧美性天天| 1769国产精品视频免费观看| 免费激情网址| 欧美日本在线播放| 亚洲精品第五页| 欧美日本在线观看| 五月激激激综合网色播免费| 国产欧美精品专区一区二区| 亚洲成年人网| 精品视频一区在线观看| 亚洲国产天堂久久综合| 日本91视频| 少妇高潮惨叫久久久久久| 国产精品亚洲αv天堂无码| 四虎永久在线精品国产免费| 国产第一页屁屁影院| 免费99精品国产自在现线| 一边摸一边做爽的视频17国产| 中文字幕在线一区二区在线| 麻豆国产精品一二三在线观看| 亚洲黄色网站视频| 波多野结衣无码中文字幕在线观看一区二区| 精品久久久久久中文字幕女| 欧美性精品| 欧美啪啪一区| а∨天堂一区中文字幕| 精品一区国产精品| 国产精品手机在线观看你懂的| 亚洲精品免费网站| 国产激爽爽爽大片在线观看| a毛片免费看| 园内精品自拍视频在线播放| 成人夜夜嗨| 美女无遮挡免费视频网站| 毛片大全免费观看| 国产亚洲一区二区三区在线| 国产00高中生在线播放| 久久久91人妻无码精品蜜桃HD| 亚洲天堂高清| 国产午夜一级毛片| 992tv国产人成在线观看| 香港一级毛片免费看| 成人亚洲天堂|