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

均勻流中近壁面垂直流向振蕩圓柱水動力特性研究*

2013-09-25 03:06:34陳鎣付世曉許玉旺周青范迪夏
物理學報 2013年6期

陳鎣 付世曉 許玉旺 周青 范迪夏

1 引言

海底管道被廣泛應用于海底油氣的運輸.但由于鋪設過程中海床的不平整、海流和海浪對海床的沖蝕作用或由于管線的變形等原因都可能在海底管道與海底之間產生一定的間隙,從而發展成懸跨管道.懸跨的出現使其在海流作用下形成周期性的旋渦脫落,從而引起渦激振動.渦激振動是決定海底管道的使用壽命、引發懸跨疲勞失效的主要原因之一.過去幾十年中,國內外對于渦激振動的研究已取得諸多進展[1-5].

關于懸跨海底管道的問題,可以簡化為均勻來流下靠近壁面的靜止或振蕩圓柱的水動力特性問題.

已有許多學者在臨界雷諾數下進行過近壁面對靜止圓柱水動力特性影響的研究,如表1所示.總結其研究內容,主要集中在不同間隙比G/D(G為管道底部離海底的距離,D為管道直徑)和邊界層的厚度δ時升力、阻力系數的變化關系及泄渦抑制的臨界間隙比,即當圓柱與壁面之間的間隙比在一定范圍內時,隨間隙比減小,平均阻力快速減小而平均升力增大,且當間隙比小于某個值時,圓柱的泄渦將受到抑制.

關于平均阻力驟降問題,Zdravkovich[6]發現當圓柱間隙比小于邊界層厚度時阻力開始減小,所以認為阻力變化的主導因素是G/δ(δ為底邊界層的厚度)而不是間隙比G/D.但Hiwada等[7]在試驗中發現當δ/D=0.23時,阻力從G/D=0.5時就開始減小,而此時圓柱仍在壁面的邊界層之外.Buresti和Lanciotti[8]在不同δ/D和雷諾數下試驗則得出阻力大小的變化與G/D沒有單調的變化關系.Lei等[9]對不同間隙比與底邊界層厚度δ下的圓柱進行了研究,認為間隙比G/D與底邊界層的厚度δ均對圓柱水動力系數有很大影響,但沒有明確的相關性.Nishino等[10]在試驗中將底邊界設為與來流速度相同的可運動平板,從而可忽略底邊界層的影響,在G/D從0.5減小到0.35時仍有阻力減小的現象.可見,在平均阻力減小的主導原因問題上仍然存在爭議.

而對于圓柱泄渦受到完全抑制的臨界間隙比的研究卻有著更好的一致性,大部分研究結果顯示臨界間隙比大約為G/D<0.3,如Bearman和Zdravkovich[11],Grass等[12],Buresti和Lanciotti[8],Lei等[9].

以上研究均是針對靜止圓柱進行的,但在現實中,當海底管道在海流作用下發生渦激振動時,其水動力特性及泄渦抑制的臨界間隙比都將發生變化,因此有必要進一步了解振蕩圓柱的水動力特性.

在過去的幾十年中,已有許多關于自由邊界孤立圓柱強迫振蕩試驗的研究成果.早期Bishop和Hassan[13]在強迫振蕩試驗中發現了鎖定現象及力與相位角在某一臨界頻率突變的現象,Wiliiamson[14]通過分析振蕩圓柱體在不同振蕩頻率和振幅下的旋渦結構解釋了這一現象.Sarpkaya[15]在數據分析中,將升力分解為與速度同相位的分量和與加速度同相位的附加慣性力分量,利用這種分解方法可以解釋在振蕩過程中流體與結構物之間的能量轉換關系.Gopalkrishnan[16]在MIT的拖曳水池做了一系列強迫振蕩試驗,并將結果繪制成云圖,為SHEAR7,VIVANA等半經驗模型提供了較為完善水動力系數.

然而,大部分強迫振蕩試驗均是針對不考慮近壁面影響的自由邊界圓柱的,關于近壁面圓柱強迫振蕩的研究較少.Sumer等[17]對在均勻流和波浪作用下垂直流向振蕩圓柱在不同的間隙比下水動力特性進行了研究,發現平均阻力系數和最大升力系數隨振蕩幅值的增大而顯著增加;平均阻力系數隨間隙比減小而減小;由于壁面的作用,在G/D=1時升力系數與速度同相位分量與自由邊界時不同.Huang和Larsen[18]使用LES(large edd simulation)方法對在不同間隙比下圓柱在同一振蕩頻率和振幅下強迫振蕩時的情況進了數值模擬,發現主要是間隙比的減小導致了阻力的突然減小,且圓柱在G/D<0.3時泄渦并沒有得到完全抑制.

綜上所述,現有的關于海底管道水動力特性的研究主要是針對靜止圓柱,缺少關于近壁面振蕩圓柱現象與機理的研究.此外,現有的水動力系數數據一般都是由自由邊界圓柱強迫振蕩試驗得到的,缺少用于預報海底管道的渦激振動的水動力系數.因此進行近壁面圓柱強迫振蕩試驗對于了解海底管道的水動力特性以及對其渦激振動的預報均有十分重要的意義.

本文在距壁面不同間隙比G/D下,在拖曳水池中對均勻來流下的光滑圓柱進行垂直流向強迫振蕩試驗研究,并詳細分析了各水動力系數隨間隙比、振蕩頻率和振幅變化的規律.

2 試驗描述

2.1 試驗裝置

試驗在拖曳水池(如圖1所示)中進行,水池尺寸為192 m×10 m×42 m(長×寬×深),拖車最大速度為9 m/s.試驗裝置如圖2所示,其包括強迫振蕩裝置、圓柱和平板三部分,強迫振蕩裝置懸掛在拖車上,圓柱安裝在其中,平板通過鋼柱懸掛在拖車上且與圓柱的軸線平行.

強迫振蕩裝置由水平滑軌系統和豎直滑軌系統兩部分組成,可用程序控制伺服電機驅動模型前后或上下滑動.

平板為焊接在骨架上的尺寸為2.5 m×2.4 m×0.003 m(長×寬×厚)的不銹鋼板,通過四根支柱連接在拖車底部.此外,在平板的四個角用四根固定在拖車上的鋼絲牽引以防止其在水中運動時產生振動.

表1 近壁面對靜止圓柱水動力特性影響的研究總結

圖1 拖曳水池整體圖

圖2 試驗裝置示意圖

2.2 儀器設備

試驗中使用三分力儀來測量模型所受的力.利用伺服電機上的編碼器實時記錄模型的位移和速度.圖3為除平板外整個試驗裝置倒置時的照片.

圖3 模型與振蕩裝置整體圖

2.3 模型

模型為一個直徑0.25 m,長2 m的光滑聚丙烯圓柱,長粗比L/D=8,在圓柱的兩端各安裝一個三分力儀以測量其受到的力,具體布置如圖4所示.試驗中使用了兩個方法來消除模型端部連接處的邊界效應:首先,兩端連接處安裝了直徑為D的假體,以保證圓柱所受的為均勻流,在圓柱與假體之間有一個很小的(<1 mm)空隙以保證所測量到的力都是圓柱所受的;其次,兩塊直徑為3D的有機玻璃圓板分別安裝在假體末端,以消除豎直軌道對流場的影響.

圖4 模型布置

2.4 試驗內容

為了研究近壁面振蕩圓柱的水動力特性,本試驗在間隙比G/D分別為0.1,0.3,0.5,0.7,0.9時,使用0.1到0.24間的8組無因次振蕩頻率(?f0=f0D/U,其中f0為強迫振蕩頻率)及不同的無因次振幅Y0(Y0=A/D)對模型進行強迫振蕩.試驗中,拖車速度恒為0.8 m/s,與之相對應的雷諾數約為Re=2×105,各參數及坐標系的示意簡圖如圖5所示,具體工況如表2所示.

圖5 試驗中坐標及各參數示意圖

每組工況至少采集120 s的數據,其中包括大約8 s的初始段,3 s的過渡段,100 s的振蕩段以及5 s的結束段.圖6是一個典型的簡諧振蕩工況阻力時歷原始數據(工況G/D=0.1,A/D=0.05,?f0=0.220).

表2 試驗工況表

圖6 工況G/D=0.1,A/D=0.05,無因次頻率為0.22時阻力時歷原始數據

3 基本方程與數據處理

當圓柱受外部力在垂直來流方向上做簡諧運動時,其受到的力中同時包含Strouhal項與強迫振蕩項.圓柱在垂直來流方向如以下形式運動:

其中,Y0為振蕩振幅,f0為振蕩頻率.

其所受到的垂直來流方向升力可表示為

其中,Lm為平均升力,Lo為振蕩升力幅值,φ0為振蕩升力相位角,Ls為Strouhal升力幅值,φs為Strouhal升力相位角.

在來流阻力方向,由于受近壁面的影響,阻力在1倍振蕩頻率 f0上也有分量,可表示為

其中,Dm為平均阻力,Don為振蕩阻力在nf0上的幅值,ψn振蕩阻力在nf0上的相位角,Ds為Strouhal阻力的幅值,ψs為Strouhal阻力的相位角.當圓柱進行強迫振蕩時,振蕩頻率項占主導,所以忽略Strouhal頻率項,再將(2)和(3)式中的每一項進行無因次化后,可得升力系數表達式為

阻力系數表達式為

其中CLm和CLo分別為平均升力系數和振蕩升力系數幅值,CDm和Don分別為平均阻力系數和振蕩阻力系數在nf0上的幅值,如圖7所示.

一個波形x(t)可以表示為如下傅里葉級數形式:

其中a0,an,bn分別為

n即表示以上系數為nf0上計算得到的分量,計算升力系數時取1,計算阻力系數時取2.

對于升力系數,可得到:

根據(10)和(11)式所求得的振蕩升力系數幅值和相位角,可將振蕩升力系數如圖8所示進行分解,其中與圓柱運動速度同相位的項表示為

將升力系數中與圓柱加速度同相位的項表示為

關于附加質量系數,按其定義可表示為

其中,MA0為附加質量,V為圓柱的體積.

對(1)式兩次微分后可得圓柱在升力方向的加速度為

根據牛頓第二定律,由(13)和(15)式可將附加質量MA0表示為

將(16)式代入(14)式即可得附加質量系數為

對于阻力系數,可得到:

由于圓柱沒有在來流方向振蕩,不對于振蕩阻力系數幅值進行繼續分解.

圖7 升力、阻力系數平均值、振蕩值示意圖

4 試驗結果與討論

本文對原始數據進行低通濾波濾除高頻噪聲后,使用以上介紹的方法處理得到所需要的各系數,對比了不同間隙比下各水動力系數隨無因次頻率和振幅的變化關系.

4.1 間隙比G/D對平均阻力的影響

圓柱在不同間隙比G/D下平均阻力系數CDm隨無因次頻率的變化關系如圖9所示.可以看出,對于振幅均為0.3的G/D=0.5,0.7,0.9三個工況,平均阻力系數隨G/D的減小而減小.這一規律同樣適用于使用不同振蕩振幅的G/D=0.1和0.3兩個工況.對于靜止圓柱,平均阻力系數CDm是一個常數,其值隨G/D減小而減小;當圓柱發生強迫振蕩時,CDm將隨無因次振蕩頻率而變化,但其在各自振蕩頻率下隨G/D減小而驟減的規律仍然適用.此外,隨G/D的減小平均阻力系數在鎖定區域附近產生的峰值逐漸消失,這與Sumer等[17]觀察到的現象相似.

圖8 升力系數向量的分解示意圖

圖9 不同間隙比G/D下平均阻力系數

取每組G/D中?f0=0.14時的CDm繪制成CDm隨G/D變化的曲線,如圖10所示.CDm隨G/D減小而減小,在0.3<G/D<0.7時變化趨勢較快,當G/D足夠大時將趨于一個定值.本文中的平板邊界層約為δ/D=0.144,所以在G/D大于0.3之后的工況圓柱均在平板邊界層δ外面而不受其影響,可見平均阻力驟減現象主要是由G/D的減小造成的,這與Hiwada等[7],Nishino等[10],Huang和Larsen[18]的研究結論相似.

圖10 平均阻力系數隨不同間隙比G/D變化曲線(?f0=0.14)

4.2 間隙比G/D對振蕩阻力系數的影響

圓柱在不同間隙比G/D下在1f0上振蕩阻力系數的幅值CDo1隨無因次頻率的變化關系如圖11所示.CDo1隨振蕩頻率的增加而增加,在f?0=0.2—0.22左右產生峰值.對于振幅均為0.3的G/D=0.5和0.7兩個工況,CDo1隨G/D的減小而增大,當G/D=0.3時,雖然振幅不相同,但CDo1繼續增大.但G/D=0.1時,CDo1的值小于其他所有工況,這是因為其無因次振幅最小.對于靜止圓柱,在G/D<0.3時,由于泄渦受到抑制,振蕩升阻力系數也隨之減小,但對振蕩圓柱卻沒有這一現象,這可能是由于圓柱振蕩使近壁面對泄渦的抑制作用減弱而導致的.對于自由邊界的強迫振蕩圓柱,Gopalkrishnan[16]的結果中CDo1的值很小,但在本試驗中由于近壁面的影響而使其變大.

圖11 不同間隙比G/D下振蕩阻力系數在1f0上幅值

圓柱在不同間隙比G/D下在2f0上振蕩阻力系數的幅值CDo2隨無因次頻率的變化關系如圖12所示,振蕩阻力系數幅值隨無因次頻率的增加而逐漸增加,除G/D=0.1和0.3的工況外,其余有近壁面的工況均和自由邊界振蕩時一樣出現了兩個峰值,但近壁面振蕩圓柱出現第一個峰值時的振蕩頻率要比自由邊界時更高,這是由于近壁面的存在導致St數增大,從而使鎖定區間往更高頻的地方移動.此外,對比圖11和圖12,CDo1普遍要大于CDo2,但隨著G/D的增加,兩者之間的差距逐漸減小.

圖12 不同間隙比G/D下振蕩阻力系數在2f0上幅值

4.3 間隙比G/D對振蕩升力系數的影響

圓柱在不同間隙比G/D下振蕩升力系數幅值CLo隨無因次頻率的變化關系如圖13所示.除G/D=0.1的工況外,其余工況均隨無因次振蕩頻率的增加而增加,在低頻區增幅較小,但在鎖定區增幅急劇變大,在高頻區有近壁面的工況數值要高于自由邊界時的情況.對于振幅均為0.3的G/D=0.5,0.7,0.9三個工況,振蕩升力系數幅值隨G/D的減小而增大,但由于振幅不同,這一規律并不適用于G/D=0.1和0.3兩個工況.當G/D=0.3時,振蕩升力系數幅值與其他更大間隙比的工況有類似的趨勢和數值,可見這時圓柱的泄渦并沒有得到抑制.當G/D繼續減小為0.1時,振蕩升力系數幅值在低頻區與其他工況相近,但從鎖定區開始,其結果就遠遠小于間隙比更大的工況.有兩個原因造成了這一現象,其一,G/D=0.1時離壁面很近,泄渦受到抑制比較嚴重;其二,振蕩幅值A/D=0.05相比于其他工況都要小很多.

圖13 不同間隙比G/D下振蕩升力系數幅值

4.4 間隙比G/D對升力相位角的影響

升力相位角為圓柱運動方向與升力之間的夾角,其決定圓柱與流體之間能量傳遞的方向,當0≤φ0≤π時,能量從流體傳向圓柱,反之當-π≤φ0≤0時能量從圓柱傳向流體,其結果如圖14所示.從圖14可看出升力相位角分上下兩個“分支”,但其在物理意義上是連續的,即升力對圓柱的作用從抑制到激勵的過渡.在低振蕩頻率時,相位角均在“下分支”且均為負值.隨著頻率的增加,相位角絕對值增大,直到達到π左右后突變為正值“跳躍”到“上分支”,之后隨頻率的增加相位角值突然減小,Williamson[14]用泄渦模式從2P到2S的突變來解釋這一現象,最后相位角的值降到0左右.

在“下分支”中,對于振幅均為0.3的G/D=0.5,0.7,0.9三個工況,相位角的數值比較相近,但G/D=0.5時發生能量傳遞方向突變的無因次頻率較另兩個工況大,且隨著間隙比繼續減小為0.3和0.1,這一頻率繼續增大.可見發生能量傳遞方向突變的無因次頻率隨G/D的減小而增大.此外,在振蕩頻率較小時,不同G/D時的相位角數值比較相近,隨頻率增大,相位角的絕對值隨G/D減小相應減小.

在“上分支”中,對于自由邊界的情況下,相位角可分為三個區間,在0.125≤?f0≤0.182和0.227≤?f0≤0.243兩個區間內為正值(激勵區),在兩者之間的區間內為負值(阻尼區),而對于有近壁面的情況則沒有這樣的區間劃分,其所有值均為正.在“上分支”中,近壁面的存在使相位角的值均比自由邊界時的大.

圖14 不同間隙比G/D下升力相位角

圖15 不同間隙比G/D下升力系數與速度同相位分量

升力系數與速度同相位分量CLV0隨無因次頻率的變化關系如圖15所示,正值表示能量從流體傳向圓柱.從低頻率區開始隨頻率增加,CLV0逐漸增加并從負值過渡到正值.對于自由邊界的情況,CLV0在f?0=0.17左右達到峰值后開始下降,隨著頻率繼續增加,在f?0=0.2時降到最小值,之后繼續增加回到正值;和相位角相對應,正值激勵區主要在0.125≤f?0≤0.182和0.227≤f?0≤0.243兩個區間.在低頻率區,對于振幅均為0.3的G/D=0.5,0.7,0.9三個工況,CLV0隨G/D的減小而減小,當G/D下降到0.3和0.1時,兩者的CLV0基本相同,但仍比G/D=0.5時有所減小.可見,在負值區CLV0隨G/D的增大而增大.此外,隨著G/D減小,CLV0值從負過渡到正的無因次頻率也隨之增加.在有近壁面時,出現峰值的頻率也比自由邊界時提高,并且CLV0值從負過渡到正值后隨沒有再出現負值.圖15中Sumer等[17]的結果也可得到相似規律.

從以上結論可看出,近壁面的存在對圓柱的能量傳遞有很大的影響,自由邊界圓柱強迫振蕩所得到的水動力系數并不能用來預報海底管道的渦激振動.

圖16 不同間隙比G/D下升力系數與加速度同相位分量

4.5 間隙比G/D對附加質量系數的影響

升力系數與加速度同相位分量CLA0及附加質量系數CM0如圖16和圖17所示.與加速度同相位的升力包含兩項:勢流附加質量力以及由泄渦引起的與加速度同相位的力.對于圓柱,前者是一個可由勢流理論解出的常量1;而后者是與泄渦強度息息相關,從而與圓柱強迫振蕩的無因次頻率有關系.從圖17中可看出,對于附加質量系數CM0,除G/D=0.1的工況外,其余工況均有相似的趨勢,即在低頻率區為負值且基本維持不變,隨振蕩頻率的增加CM0逐漸增大,最后在0.22≤f?0≤0.24內基本保持穩定.可見CM0的值只有在某個特定的振蕩頻率區間內才是定值,正是由于在不同的振蕩頻率下泄渦引起的與加速度同相位的力不同所造成的.在鎖定區域之前的低頻率區,振幅相同的G/D=0.5,0.7和0.9三個工況中的CM0值隨G/D減小而減小,而當G/D繼續減小為振幅不同的0.3和0.1兩個工況時,這一規律仍然適用.此外,圖16中,關于升力系數與加速度同相位分量的結果與Sumer等[17]的G/D=1的結果比較接近;圖17中附加質量系數隨G/D減小而減小的規律與Yamamoto等[20]用勢流理論得到的結論相似.

圖17 不同間隙比G/D下附加質量系數

4.6 不同振蕩幅值對升阻力系數的影響

圖18 為不同無因次振幅A/D對升力和阻力系數的影響的對比圖.可以看出,CDm,CDon,CLo均隨無因次振幅A/D的增大而增大.對于振蕩升力和阻力系數,是因為當圓柱朝壁面運動時,由于間隙比突然減小而導致圓柱靠近壁面一側的流速增大,使圓柱表面壓強減小,所以產生了升力極值,更大的振蕩幅值會產生更小的間隙,從而使振蕩值隨振幅增加而變大.這一結論與Bishop和Hassan[13]及Sarpkaya[15]的關于自由邊界振蕩圓柱的結果相似.對于平均阻力系數CDm,Sarpkaya[15]認為是由于A/D增大使圓柱在來流方向的投影面積增大,從而導致圓柱在來流方向受到的平均阻力增加.

圖18 G/D=0.5時,無因次振幅A/D對各系數的影響 (a)平均阻力系數;(b)1f0上振蕩阻力系數;(c)2f0上振蕩阻力系數;(d)振蕩升力系數幅值

5 結論

本文在不同間隙比G/D下,對均勻來流下垂直流向強迫振蕩光滑圓柱的水動力系數進行了試驗研究.首先對模型試驗進行了詳細描述,然后對各系數及試驗數據處理方法進行了說明,詳細分析了各水動力系數隨間隙比、振蕩頻率和振幅變化的規律,得到以下結論.

1)對于近壁面振蕩圓柱,平均阻力系數在低無因次頻率區與靜止圓柱的結果相近,且在間隙比G/D從0.7到0.3有一個快速減小的過程,此時的圓柱完全在平板邊界層之外,這意味著近壁圓柱阻力驟降現象主要是由間隙比減小引起的.

2)振蕩圓柱泄渦受到完全抑制的臨界間隙比要小于靜止圓柱.

3)在近壁面的作用下,升力相位角從負突變到正的振蕩頻率增加,且間隙比G/D越小,突變頻率越大.由于近壁面的影響,升力相位角在高頻率區不再出現負值且比自由邊界時大,在低頻率區則比自由邊界時小.可見近壁面的存在對振蕩圓柱的能量傳遞有著重要的影響,自由邊界圓柱強迫振蕩所得到的水動力系數不能用來預報海底管道的渦激振動.

4)對于振蕩圓柱,附加質量系數只有在一定的頻率范圍內才是定值.在低頻率區域,附加質量系數的絕對值隨間隙比G/D減小而增大.

5)圓柱在進行強迫振蕩時,其平均阻力系數、振蕩阻力系數和振蕩升力系數均隨無因次振幅A/D的增加而增大.

[1]Sarpkaya T 2004 J.Fluid Struct.19 389

[2]Naudascher E 1987 J.Fluid Struct.1 265

[3]Chen Z H,Fan B C,Zhou B M,Li H Z 2007 Chin.Phys.Soc.16 1077

[4]Xia Y,Lu D T,Liu Y,Xu Y S 2009 Chin.Phys.Lett.26 034702

[5]Xu W H,Du J,Yu J X,Li J C 2011 Chin.Phys.Lett.28 124704

[6]Zdravkovich M M 1985 Appl.Ocean Res.7 197

[7]Hiwada M,Mabuchi I,Kumada M,Iwakoshi H 1986 Trans.Jpn.Soc.Mech.Eng.B 52 25

[8]Buresti G,Lanciotti A 1992 J.Wind Eng.Ind.Aerodyn.41 639

[9]Lei C,Cheng L,Kavanagh K 1999 J.Wind Eng.Ind.Aerodyn.80 263

[10]Nishino T,Roberts G T,Zhang X 2007 Phys.Fluids 19 025103

[11]Bearman P W,Zdravkovich M M 1978 J.Fluid Mech.89 33

[12]Grass A J,Raven P W J,Stuart R J,Bray J A 1984 J.Energ Resour-ASME 106 70

[13]Bishop R E D,Hassan A Y 1964 Proc.Roy.Soc.Lond.A 277 51

[14]Williamson C H K,Roshko A 1988 J.Fluid Struct.2 355

[15]Sarpkaya T 1978 ASCE J.Waterway Port Coast.Ocean Div.104 275

[16]Gopalkrishnan R 1993 Ph.D.Dissertation(Cambridge,MA,USA:MIT)

[17]Sumer B M,Fredsoe J,Jensen B L,Christiansen N 1994 J.Waterway.Port.C.-ASCE 120 233

[18]HuangZY,LarsenCM201029thInternationalConferenceonOcean,Offshore and Arctic Engineering Shanghai,China,2010 20006

[19]Roshko A,Steinolfson A,Chattoorgoon V 1975 Proceedings of the 2nd U.S.National Conference on Wind Engineering Research Fort Collins,USA,1975 IV 15

[20]Yamamoto T,Nath J H,Slotta L S 1974 ASCE J.Waterway Port.Coast.Ocean Engng.Div.100 345

主站蜘蛛池模板: 婷婷激情五月网| 亚洲免费福利视频| 成年人免费国产视频| 国内精品视频区在线2021| 日本精品中文字幕在线不卡| 国产丝袜无码精品| 国产真实乱了在线播放| 亚洲欧美国产视频| 亚洲天堂日韩在线| 亚洲国产精品无码久久一线| 免费看av在线网站网址| 国内自拍久第一页| 无码免费的亚洲视频| A级毛片高清免费视频就| 免费久久一级欧美特大黄| 欧美成人免费一区在线播放| 色综合成人| 国产产在线精品亚洲aavv| 午夜少妇精品视频小电影| 国产欧美日韩18| 激情影院内射美女| 国产十八禁在线观看免费| 国产波多野结衣中文在线播放| 熟妇丰满人妻| 欧美色香蕉| 亚洲手机在线| 午夜色综合| 国产精品自在拍首页视频8| 在线播放91| 婷婷午夜影院| 国产精品成人一区二区不卡| 久久久久无码精品| yjizz国产在线视频网| 中国丰满人妻无码束缚啪啪| 高清无码一本到东京热| 97人人做人人爽香蕉精品| 精品国产成人国产在线| 思思99热精品在线| 亚洲精品无码久久毛片波多野吉| 国产精品冒白浆免费视频| 人人澡人人爽欧美一区| 国产三级毛片| 亚洲AV电影不卡在线观看| 国产人在线成免费视频| 中字无码av在线电影| 欧美人人干| 国产亚洲欧美另类一区二区| 国产白浆在线| 精品无码国产一区二区三区AV| 在线观看无码av五月花| 欧美精品v欧洲精品| 久久五月视频| 91久久青青草原精品国产| 久久国语对白| 欧美中文一区| 国产乱人激情H在线观看| 亚洲二三区| 视频二区国产精品职场同事| 免费观看国产小粉嫩喷水| 国产成人精品男人的天堂下载 | 天堂在线视频精品| 99re这里只有国产中文精品国产精品 | 亚洲综合色婷婷中文字幕| 国产微拍精品| AV在线天堂进入| 日韩在线视频网| 亚洲综合天堂网| 久久国产成人精品国产成人亚洲 | 国产青榴视频| 99热国产在线精品99| 内射人妻无码色AV天堂| 久久这里只精品国产99热8| 国产丰满成熟女性性满足视频| 国产亚洲美日韩AV中文字幕无码成人| 国产丝袜啪啪| 亚洲va视频| 国产成人精品一区二区三区| 美女被躁出白浆视频播放| 亚洲第一成年网| 91无码国产视频| 特级精品毛片免费观看| 日韩a级毛片|