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

A柱渦動力學特性隨前窗傾角變化研究

2019-01-08 07:53:20袁海東楊志剛李啟良
同濟大學學報(自然科學版) 2018年12期
關鍵詞:模型

袁海東,楊志剛,3,李啟良

(1.同濟大學 上海地面交通工具風洞中心,上海 201804; 2.同濟大學 上海市地面交通工具空氣動力與熱環(huán)境模擬重點實驗室,上海 201804; 3.北京民用飛機技術研究中心,北京 102211)

運動中的汽車周圍存在兩種典型的流動分離形式[1].一種是在車體尾部,分離線垂直于流動方向,導致準二維分離,它的動能大部分通過湍流參混耗散,從而導致壓力損失和空氣阻力.另一種類型的分離,其特征是存在相對于來流傾斜的分離線,本質上是三維的.這種分離流動含有豐富的動能,渦核軸線大致沿著流線方向與其形成小夾角,并向下游輸送動能.類似的縱向渦可以在汽車的A柱和C柱區(qū)域產生,也存在于三角翼[2]和高速列車的外流場[3].

首先,從風噪聲的角度開始對A柱渦的拓撲結構進行了試驗研究[4]. A柱后面的分離包括一個起主導作用的主渦和一個位于A柱后面的較小尺度的二次渦.Haruna等[5]的研究證實了A柱渦的拓撲結構,表明分離的區(qū)域在風噪聲中起著重要作用,而并不是再附著流動的區(qū)域.Alam等[6-7]對A柱幾何形狀對側窗周圍流動的影響進行了參數化研究,表明側窗上的壓力波動大小主要取決于局部A柱弧度.為了深入研究A柱渦與風噪之間的關系,A柱渦的拓撲結構和動力學特性通過簡化的A柱模型進行了試驗和數值研究.Levy等[8]利用粒子圖像測速儀(PIV)測量了Dihedron模型的A柱渦,其研究證實,A柱渦是一個相對穩(wěn)定的縱向渦,由主渦和二次渦組成,局部最大不穩(wěn)定性存在于主渦、剪切層和二次分離線位置.后來的研究[9]表明,主要的壁面壓力脈動位于二次分離線位置,并且A柱渦不會自然破裂.Hoarau等[10]對A柱模型側壁脈動壓力場的時空特性進行了研究,指出了3個擾動形式是產生側壁壓力脈動的原因.這些擾動對應于錐形渦結構的長度尺度、寬度尺度和富含能量的高頻部分.

30° Dihedron模型[8]的A柱渦是相對穩(wěn)定的縱向渦,當前窗傾角為45°時[12]發(fā)生渦破裂.這顯示了前窗傾角與A柱渦穩(wěn)定性之間的關系,類似于三角翼[2]的情況. Gilliéron等[16]從理論上分析了渦破裂的物理機制,重點研究前窗傾角的影響.本文首先通過試驗和數值的方法研究了A柱渦的動力學演化過程隨前窗傾角的變化,然后分析了渦破裂對氣動阻力和壁面動力學特性的影響.

1 方法

1.1 Dihedron模型

在上述討論的研究中,Dihedron鈍體模型已被用于A柱渦的動力學特性研究.Dihedron模型代表了簡化車體的前部.它有3個固定的參數(長度L=0.38 m、寬度W=0.13 m、高度H=0.096 m),以及可變參數前窗傾角θ,見圖1.寬度和高度的比例與Ahmed 類車體模型相同[17].L0是與Dihedron模型相關的特征長度,其定義為模型橫截面的平方根,即L0=0.112 m.前窗傾角范圍從30°~90°,分別對應轎車、SUV、MPV、客車和卡車,間隔為15°.

1.2 風洞試驗

試驗在開口式回流風洞中進行,該風洞具有0.9 m長的測試段,其橫截面寬度為0.43 m,高度為0.28 m,如圖1所示.Dihedron模型安裝在測試段上,距離噴口為0.2 m.測試部分的湍流強度為0.5%,試驗阻塞比為10%.該模型由NACA 0020支撐柱支撐,距地面60 mm,以最大限度地減少Dihedron模型與地面之間的相互作用.Yuan等[18]描述了風洞的更多細節(jié).噴口的自由來流速度為U0=24 m·s-1,由Dihedron模型的特征長度定義的雷諾數為Re=0.12×106,這與文獻[8]研究相同.模型安裝保證為零偏航.

圖1 Dihedron模型試驗布置示意圖Fig.1 Schematic of experimental setup of a Dihedron model

使用NovaSensor固態(tài)低壓傳感器測量Dihedron模型表面的靜態(tài)壓力.在每個測量位置,以500 Hz采樣率采集4 000個壓力數據.將表面壓力p與參考壓力p0之間的壓力差除以動態(tài)壓力,以給出壓力系數Cp,如等式(1)所示.在地板下安裝六分量天平以測量Dihedron模型的氣動阻力.阻力系數CD表示為等式(2).

(1)

(2)

式中:ρ為空氣密度;Fx為作用在Dihedron模型上的氣動力的x方向分量;U0為噴口的自由來流速度;A為Dihedron模型的正投影面積.

1.3 數值計算

為了捕捉非定常的流動特征、解決振動和噪聲問題,混合RANS / LES(雷諾時均納維-斯托克斯方程/大渦模擬)方法得到迅速發(fā)展,并且分離渦模擬(DES)是最受歡迎的方法之一.DES 97[19]的原始策略為三維非定常數值解,亞格子網格尺度模型在網格密度對于LES足夠精細的區(qū)域中起作用,并且在其余區(qū)域使用雷諾平均模型.最近的改進[20],如延遲DES(DDES)[21]和改進的DDES(IDDES)[22]增加了額外的功能來解決網格誘導分離(GIS)和對數層不匹配(LLM)的弱點,使DES在復雜流動解決方案中更加強大. DES的可行性已經在空氣動力學[23]和氣動聲學[24]方面得到證明.

所有模擬均在西門子公司開發(fā)的商業(yè)代碼STAR-CCM+11.06.010-R8中完成.該解決方案使用穩(wěn)定的RANS結果進行初始化.時間步長為2×10-5s.在30個對流單元之后采集30個對流單元(L/U0)的時間平均統(tǒng)計結果.計算域的長、寬、高為15L、11W、12H.在模型前4L位置施加自由來流速度U0的速度入口,入口邊界處的湍流由湍流強度和湍流黏度比確定,分別為0.005和10.在Dihedron模型后面的下游10L處施加0 Pa標準壓力出口條件.其余邊界設置為對稱邊界條件,并且Dihedron模型被置于計算域橫截面的中心.Dihedron模型表面采用無滑移壁面邊界條件.

表1 粗細網格的空間分辨率Tab.1 Spatial resolution for coarse and fine grids

圖2 圍繞Dihedron模型的計算網格Fig.2 Computational grid around Dihedron model and zoom of boundary layer

1.4 計算有效性

通過對比30° Dihedron模型的計算和試驗結果來驗證計算的有效性.表2對比了阻力系數的試驗值和計算值,其中風洞試驗的阻力系數通過公式(3)[25]進行了修正.

(3)

式中:Cdcor為阻力系數的修正值;Cdm為阻力系數的測量值;ΔCdHB為受水平浮力影響的阻力修正值;εS為對射流擴張和偏移的修正;εN為對于噴口阻塞效應的修正;εC為對于收集口阻塞效應的修正.由表2中結果可以看出,3套網格與修正后的試驗值相差在5%以內,網格1和網格3結果接近,并且與網格2差距比較大.

表2 阻力系數的計算值與試驗值的對比Tab.2 Drag coefficient from simulations compared with experimental results

注:ΔCD1是計算值和試驗值之間的差異; ΔCD2是網格1和網格2分別與網格3之間的差異.

圖3顯示了模型表面的平均壓力系數,圖中x和z為水平和豎直方向坐標,原點位于模型對稱面前端頂點.由圖3可知,網格1和網格3與試驗結果吻合較好,網格2與試驗結果差異較大,由此可以得出,在近壁區(qū)域的網格加密對于準確捕獲表面壓力是非常有必要的.

a 正面b 頂面c 背面d 側面z/H=0.4

圖3 縱向對稱線和側面的平均壓力系數與試驗值的比較Fig.3 Comparison of mean pressure coefficient on longitudinal symmetry lines and side surface with experimental values

圖4對比了A柱渦橫剖面上的平均速度投影和流線圖的試驗與計算結果.其中,試驗結果參考了文獻[8]中的PIV測量結果,測量平面為文獻中的平面5,V為速度在測量平面內的分量.盡管3套網格都成功預測了主渦和二次渦,但網格2獲得的主渦尺寸和位置與試驗結果不符,網格3較網格2更接近試驗結果.

湍流黏性比νt/ν,表示亞格子尺度(模擬的湍流)對整體耗散水平的貢獻量.其中νt為湍流黏性,ν為空氣的運動黏性.如圖5所示,該值隨網格細化而降低.根據文獻[26],網格3獲得的湍流黏性比小于10,滿足本文研究的要求.綜合以上試驗與計算結果的對比,選擇網格3為本文的網格生成策略.

a 試驗b 網格1c 網格2d 網格3

圖4時均流線與速度云圖
Fig.4Streamlinescalculatedbymeanflowfieldandnormalizedmeanvelocity

b 網格1c 網格2d 網格3

2 結果與討論

2.1 流動分離和局部坐標系的定義

由圖6中速度梯度的二階不變量Q[27]的等值面圖可知,A柱鋒利的邊緣強制附近的流動發(fā)生分離.流動分離對前窗傾角非常敏感,明顯地區(qū)分為3種不同的流動分離形式.對于30° Dihedron模型,一個柱狀縱向渦沿著A柱邊緣形成并延伸直到模型頂部發(fā)生破裂.對于45°的前窗傾角,A柱渦的破裂提前到縱向渦形成的初始段.繼續(xù)增大前窗傾角,在整個A柱區(qū)域形成離散的展向小渦.

渦核軸線是流動旋轉所圍繞的軸線.渦核軸線精確描述了渦結構的位置并與其他渦結構進行區(qū)分.基于時均速度場,根據Haimes和Sujudi[28-29]提出的提取方法計算渦核,如圖7所示.α和β分別表示渦核軸線與來流方向的夾角在豎直和水平平面的投影.基于渦核軸線建立的局部坐標系如圖8所示.選取一條直線讓渦核盡量分布在這條直線的周圍.這條直線被模型的上下表面截得的線段定義為A柱渦的特征長度Lv.局部坐標系的原點設在這條線段的起點,x軸與這條線段重合.局部坐標系是在原始坐標系的基礎上分別繞y0軸和新的z1軸旋轉-α和β得到的.下文的分析是基于局部坐標系展開的.

a θ=30°b θ=45°c θ=60°d θ=75°e θ=90°

圖6瞬時的Q=4×106等值面圖
Fig.6Iso-surfaceofQ=4×106

a θ=30°b θ=45°c θ=60°d θ=75°e θ=90°

2.2 渦破裂

沿A柱渦主渦渦核軸線靜壓系數變化曲線如圖9所示.這些壓力系數曲線同樣可以分成3類.對于45° Dihedron模型,在A柱渦的整個長度尺度內存在明顯的逆壓梯度,暗示了渦破裂的發(fā)生.而對于30° Dihedron模型,壓力系數曲線表現(xiàn)為多峰狀,表明渦破裂得到抑制.對于前窗傾角更大的另外3個模型,壓力系數變化不大并維持在較高的水平.

a 局部坐標系b 后處理面與線

圖9 靜壓系數沿主渦渦核軸線的演化Fig.9 Evolution of static pressure coefficient along primary vortex core axis

A柱渦沿渦核軸線的對流速度Vx/U0可以通過提取A柱渦的特征線沿x方向的速度分量獲得,如圖10所示.對于不同的前窗傾角,A柱渦的對流速度有明顯的區(qū)別.對于30°模型,A柱渦在x/Lv=0.1~0.8區(qū)域對流速度都高于自由來流速度,在x/Lv=0.1處為最大值1.4.較大的對流速度表明A柱渦是一個富含動能的縱向渦,對于動量和渦量的輸運效率較高.而對于較大的前窗傾角60°、75°和90°,A柱渦對流速度幾乎消失.一個明顯的過渡階段存在于前窗傾角為45°的模型中.

圖10 軸向速度沿主渦渦核軸線的演化Fig.10 Evolution of normalized axial velocity along primary vortex core axis

圖11和12顯示了A柱渦不同位置的旋轉速度和軸向速度剖面,揭露了A柱渦空間不對稱的本質特征.這一現(xiàn)象類似于三角翼的前緣渦和Ahmed類車體的C柱渦,但明顯不同于自由旋轉射流.模型的側表面和A柱附近的剪切層是產生非對稱性的主要原因.對于不同前窗傾角的Dihedron模型,軸向速度剖面表現(xiàn)為不同的形狀,比如30°模型為射流狀,而60°模型為尾跡狀,而對于處于過渡階段的45°模型兩種速度剖面形狀同時存在.30°傾角的模型旋轉速度表現(xiàn)為蘭金渦的特征,與較大傾角的模型區(qū)別比較明顯.從圖10~12可以看出,對于30°傾角模型,A柱附近的速度高于自由來流速度20%,在關于A柱渦[8]和C柱渦[30]的研究中同樣得到過類似的結論.圖11所示,對于60°傾角模型軸向速度存在明顯的駐點,這暗示著縱向渦已經破裂.

從圖14所示的時均速度云圖和流線圖可以看出,對于30° Dihedron模型,A柱渦包含剪切層、主渦和二次渦,這和實車試驗結果[5]以及過去關于30° Dihedron模型的研究結論[8]一致.而對于較大前窗傾角的Dihedron模型以及45°模型中截面高度高于x/Lv=0.4的位置,A柱渦中的二次渦消失.由于二次渦是由主渦誘導產生的,因此二次渦的消失表明主渦的強度減弱或者不能穩(wěn)定存在.

a θ=30°

b θ=45°

c θ=60°

圖11 沿直線z=0的軸向和旋轉速度剖面Fig.11 Axial and azimuthal velocity profiles along lines z=0

a θ=30°

b θ=45°

c θ=60°

圖12 沿直線y=0的軸向和旋轉速度剖面Fig.12 Axial and azimuthal velocity profiles along lines y=0

a θ=30°b θ=45°c θ=60°d θ=75°e θ=90°

圖13 不同傾角模型在平面x/Lv為0.2、0.4、0.6、0.8上的旋流數Fig.13 Ratio of azimuthal and axial velocity on planes at x/Lv=0.2, 0.4, 0.6, and 0.8

a θ=30°b θ=45°c θ=60°

d θ=75°e θ=90°

圖14 不同傾角模型在平面x/Lv為0.2、0.4、0.6、0.8上的時均速度云圖和流線圖Fig.14 Streamlines calculated by mean field on section planes at x/Lv=0.2, 0.4,0.6, and 0.8 colored by normalized mean w

圖15顯示了不同前窗傾角的Dihedron模型在不同高度截面上x方向的瞬時渦量云圖.對于30°Dihedron模型,沿著A柱高度方向,剪切層中不斷產生的渦量注入到主渦中,并在主渦中積累.沿著A柱高度方向,主渦的尺寸逐漸增大,類似的描述在levy等[8]的研究中也出現(xiàn)過.對于30° Dihedron模型,在x/Lv=0.4的位置主渦開始破裂為離散的小渦.對于前窗傾角更大的Dihedron模型,主渦的破裂更明顯,A柱渦變成由剪切層和離散的展向小渦組成.

a θ=30°b θ=45°c θ=60°

d θ=75°e θ=90°

圖15 不同傾角模型在平面x/Lv為0.2、0.4、0.6、0.8上的瞬時渦量云圖Fig.15 Transient vorticity on section planes at x/Lv=0.2, 0.4, 0.6, and 0.8

為了解釋A柱渦如何破裂為離散的小渦,圖16同時顯示了x方向的瞬時渦量和流線圖,u為沿渦核軸線的對流速度.對于30° Dihedron模型,沿著A柱高度方向,剪切層中產生的渦量不斷注入并聚集在主渦中.同時A柱渦沿主渦軸線較高的對流速度具有較強的渦量輸運效率,這使得A柱渦主渦中的渦量保持平衡.A柱渦以相對穩(wěn)定的縱向渦結構存在,并在近壁區(qū)域誘導產生二次渦.對于45°Dihedron模型,前窗傾角增大使得A柱渦沿渦核軸線的對流速度降低,對于渦量的輸運效率降低.但是剪切層不斷產生的渦量使得A柱渦中渦量不斷聚集,為了保持渦量平衡,A柱渦只能破裂為離散的小渦,同時二次渦消失.圖17描述了A柱渦存在的兩種不同形態(tài),分別為穩(wěn)定的縱向渦和不穩(wěn)定的展向渦.類似地,Wu等[34]用渦量平衡理論解釋了圓柱繞流的流動形態(tài)隨雷諾數的變化.

a θ=30°b θ=45°

圖16不同傾角模型在平面x/Lv為0.2、0.4、0.6、0.8上的瞬時渦量和瞬時流線
Fig.16Transientvorticityonsectionplanesatx/Lv=0.2,0.4,0.6,0.8andtransientstreamlinescoloredbyu

2.3 氣動阻力

圖18中的試驗和計算結果表明,Dihedron模型的氣動阻力隨著前窗傾角增大而明顯增大.從計算結果可以看出,氣動阻力主要來自壓差阻力.圖19顯示了模型表面的壓力云圖試驗和計算結果的對比.模型表面的壓力系數試驗和計算結果吻合較好,其中45°模型差異稍大,根據上面的分析可知,對于45°模型來說,A柱渦的流動形態(tài)處在過渡階段,模型尾部的壓力隨著前窗傾角的變化沒有明顯的差異(這里并沒有給出).由此可知,Dihedron模型氣動阻力的差異主要來自前窗表面壓力的差異.

隨著前窗傾角增大,較大的前窗表面壓力是阻力增加的主要原因.在模型的頂面、底面和側面存在較大的分離區(qū)域,這些分離區(qū)也隨著前窗傾角的增加而增大.較大的分離區(qū)阻礙了前方來流順利地流向模型尾部,動量在模型前部積聚從而導致了阻力的增加.另外,隨著前窗傾角的增大,A柱渦的對流速度降低,按照上面的分析,A柱渦的流動形態(tài)由穩(wěn)定的縱向渦轉換成不穩(wěn)定的展向離散小渦,A柱渦的動量輸運效率降低,從而導致動量在模型前部的聚集,這也是導致阻力增加的原因.

a A柱渦

b 穩(wěn)定的縱向渦c 不穩(wěn)定的展向渦

圖18 時均阻力系數的試驗值和計算值對比Fig.18 Mean drag coefficient from experimental and numerical results

氣動阻力的明顯差異表明流動結構存在明顯區(qū)別.關于Ahmed類車體的C柱渦進行過類似的討論[17,35].增大后背傾角,阻力增加并在30°時發(fā)生阻力驟降.這與C柱渦的動力學特性有關.富含動量的C柱渦帶走尾部的動量降低背壓產生阻力,當后背傾角大于30°時,C柱渦破裂成離散小渦,動量被留在尾部提升背壓,從而阻力降低.關于隨前窗傾角變化氣動阻力和A柱渦動力學特性的討論還比較少.

Gilliéron等[36]從理論上分析了三維鈍體傾斜角度對縱向渦形成的影響.如圖20所示,在渦核軸線和再附線之間靠近壁面的位置的合成速度V是由主流速度Vx和縱向渦誘導速度Vθ組成的.對于A柱渦和C柱渦,再附線分別位于背風側和迎風側,因此合成速度V分別趨向于遠離渦核軸線和靠近渦核軸線.這導致了相比于A柱渦,C柱渦有更充足的動量供應從而能維持為穩(wěn)定的縱向渦,而A柱渦更容易破裂為離散的展向小渦.

a 試驗,θ=30°b 計算,θ=30°c 試驗,θ=15°d 計算,θ=45°e 試驗,θ=60°f 計算,θ=60°g 試驗,θ=75°h 計算,θ=75°i 試驗,θ=90°j 計算,θ=90°

圖19 時均表面壓力系數的試驗值與計算值對比Fig.19 Mean surface pressure coefficient from experimental and numerical results

a A柱渦

b C柱渦圖20 縱向渦破裂趨勢示意圖[35]Fig.20 Schematic of longitudinal vortexes breakdown trend[35]

從減阻的角度來說,A柱渦盡量維持為穩(wěn)定的縱向渦,更多地將模型前部的動量輸運到模型尾部;而C柱渦盡量破裂為離散的展向渦,從而把動量留在模型尾部,通過降低模型前部的壓力,增加模型尾部的壓力,從而減小阻力.但是,根據前面的分析,從減阻角度對縱向渦動力學特性的期望和自然狀態(tài)下縱向渦的動力學特性存在明顯的差異.現(xiàn)有的一些定常的或非定常吹吸氣控制在A柱渦和C柱渦的控制上取得明顯的效果[37-38].因此,在A柱和C柱位置聯(lián)合使用吹吸氣控制是一種潛在高效的減阻方案,同時保證流量平衡.這還需要進一步的研究.

2.4 壁面動力學特性

現(xiàn)有的研究[6]已經表明,A柱附近的流動是影響車內舒適性的重要氣動噪聲源.車身表面的湍流壓力脈動和聲學壓力脈動是產生車內氣動噪聲的兩種壓力激勵,目前有很多關于兩種激勵對車內噪聲貢獻相對大小的研究.圖21顯示了湍流壓力的脈動和表面時均流線.壓力脈動主要出現(xiàn)在分離區(qū)域靠近再附線的位置,這和現(xiàn)有的試驗研究吻合[5].通常,對于30° Dihedron模型,壓力脈動最大值出現(xiàn)在二次分離線和再附線之間,Levy等[9]得出了類似的結論.不同前窗傾角模型中,45°模型側面的壓力脈動最大,這證實了該區(qū)域渦破裂的發(fā)生.

a θ=30°b θ=45°c θ=60°d θ=75°e θ=90°

圖21 壁面流線和壓力脈動Fig.21 Surface streamlines calculated by mean wall shear stress and root mean square of pressure coefficient on Dihedron model

在汽車行業(yè),相比于遠場噪聲車內噪聲受到更多的關注.湍流壓力脈動可以激勵側窗振動從而向車內輻射噪聲.除此之外,根據渦聲理論[39],A柱渦本身可以產生空氣動力學噪聲.前側窗區(qū)域無論湍流壓力脈動還是聲學壓力脈動都會隨前窗傾角變化而變化.現(xiàn)有的研究傾向性地認為相比于湍流壓力脈動,聲學壓力脈動對車內的噪聲貢獻量更大[40-41].因此,前窗傾角間接合理地影響車內噪聲水平.

3 結論

(1)本文通過分離渦模擬和風洞試驗研究了Dehidron模型前窗傾角對A柱渦動力學演化的影響,風洞試驗結果證實了數值方法的有效性.

(2)數值計算結果表明,隨著前窗傾角的增大,A柱渦發(fā)生渦破裂,A柱渦由穩(wěn)定的縱向渦轉化為不穩(wěn)定的展向渦.前窗傾角增大降低了A柱渦沿渦核軸線的渦量輸運效率,受主渦中渦量平衡的影響,A柱渦破裂為離散的展向小渦.

(3)試驗和計算結果都表明,隨著前窗傾角的增加,Dihedron模型的阻力增大,這與A柱渦由縱向渦向展向渦的轉化有關.進一步分析表明,減阻對類車體周圍縱向渦動力學特征的訴求和自然狀態(tài)下縱向渦的動力學特征矛盾.

(4)模型壁面動力學特征分析表明,前窗傾角間接合理地影響著車內的氣動噪聲水平.

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 精品少妇人妻无码久久| 欧美日韩北条麻妃一区二区| 国产激爽大片高清在线观看| 中国毛片网| 亚洲国产精品人久久电影| 伊人网址在线| 任我操在线视频| 久久伊人操| 精品亚洲欧美中文字幕在线看| 国产成人亚洲精品无码电影| 欧美色99| 女人18毛片久久| 国产精品19p| 2019国产在线| 亚洲福利片无码最新在线播放 | 日本午夜精品一本在线观看| 亚洲高清中文字幕| 国产va免费精品观看| 国产美女一级毛片| 国产精品女在线观看| 毛片在线看网站| 黄色福利在线| 国产va视频| 精品人妻无码区在线视频| 国内嫩模私拍精品视频| 国产精品第一区| 亚洲欧洲日产国码无码av喷潮| 五月婷婷伊人网| 毛片在线播放a| 自拍亚洲欧美精品| 日本高清有码人妻| 精品人妻一区无码视频| 久久综合色88| 91青青草视频在线观看的| 欧美第二区| 亚洲精品无码抽插日韩| 亚洲天堂.com| 伊人查蕉在线观看国产精品| 亚洲日韩精品无码专区97| 夜夜操天天摸| 国产内射一区亚洲| www.91中文字幕| 亚洲婷婷在线视频| 中文字幕精品一区二区三区视频| 3344在线观看无码| 性欧美久久| 国产视频久久久久| 久久综合九色综合97网| 欧美日韩激情在线| 久久毛片网| 最新无码专区超级碰碰碰| 国产在线一二三区| 欧美午夜在线观看| 久久黄色影院| 999精品视频在线| 国产视频 第一页| 精品视频一区二区三区在线播| 国产主播一区二区三区| 精品国产香蕉伊思人在线| 国产精品天干天干在线观看| 91亚洲精选| 无码精品国产VA在线观看DVD| 免费毛片在线| 亚洲国产天堂久久综合| 日韩AV手机在线观看蜜芽| 在线观看国产精品日本不卡网| 精品色综合| 国产成人精品在线1区| 久久超级碰| 在线播放91| 成人免费网站久久久| 中文精品久久久久国产网址 | 欧美精品影院| 国产毛片网站| a亚洲视频| 欧美第二区| 国产午夜福利在线小视频| 国产成人亚洲精品蜜芽影院| 最新国产在线| 欧美成a人片在线观看| 91免费观看视频| 成人伊人色一区二区三区|