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

考慮軸向張力時變效應的圓柱體渦激振動響應特性研究

2019-05-27 02:24:34袁昱超薛鴻祥唐文勇
振動與沖擊 2019年9期
關鍵詞:模態(tài)振動結構

袁昱超,薛鴻祥,唐文勇

(1.上海交通大學 海洋工程國家重點實驗室,上海 200240;2.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心,上海 200240)

立管作為油氣生產系統(tǒng)的重要組成部分,起到將油氣資源由海底輸送至頂端平臺的作用,典型的平臺-頂張式立管-海床系統(tǒng)可簡化為圖1所示布置形式。影響立管動力響應的結構剛度可分為由立管固有屬性決定的彎曲剛度和由軸向張力提供的附加剛度兩種成分。在復雜波浪環(huán)境中,浮式平臺易發(fā)生垂蕩運動,誘發(fā)張緊器壓縮或拉伸,從而導致立管頂端張力隨時間波動。因此,考慮張力時變效應的細長圓柱體渦激振動相較張力恒定條件更接近于實際海洋環(huán)境。

近年來,Franzini等[1-2]開展了張力簡諧變化的小尺度立管模型試驗研究。Karniadakis等[3-4]基于二維切片理論借助CFD方法研究了結構彎曲剛度可變條件下圓柱體渦激振動問題。Park等[5]認為由時變張力引發(fā)的參數激勵會改變結構原有響應特征。Da Silveira等[6]發(fā)現張力時變條件下立管動力響應存在模態(tài)階躍。王東耀等[7]、Wu等[8]和Chen等[9]均指出平臺垂蕩引起的張力波動可能導致立管激發(fā)更高階模態(tài)及更大幅度動力響應。唐友剛等[10]研究了時變張力對于剪切流工況下TTR立管渦激振動響應的影響效應并得出結構響應存在對應0.5倍參激頻率的亞諧振成分。

圖1 平臺-頂張式立管-海床系統(tǒng)示意圖Fig.1 Sketch of platform-TTR-seabed system

流體力分解模型是廣泛用于預報圓柱體渦激振動的一類方法,詳見文獻[11-14]。以往流體力分解模型并未考慮時變張力對結構動力響應的影響。本文基于Wang等提出的時域流體力分解模型,在每一分析步之初更新結構剛度矩陣以計及張力時變效應。采用Franzini等公布的小尺度模型試驗結果驗證本文方法的有效性。結合另一38 m細長圓柱體模型,研究了28組工況下結構渦激振動響應特性,討論了時變張力的幅值和頻率對響應特性的影響規(guī)律。

1 渦激振動時域模型

本文選取Cartesian坐標系,其中x軸順應來流速度方向,y軸垂直于流速方向,z軸為圓柱體軸向。

1.1 渦激振動水動力載荷模型

基于流體力分解模型,渦激振動橫流向流體力載荷Fhydro,y可由式(1)表示,主要包含三個載荷成分,即激勵力Fv,阻尼力Fd和慣性力Fm。

(1)

式中:Cv為激勵力系數,A為響應位移幅值,D為直徑,f為響應頻率,V為遭遇流速,t為時間,ρf為流體密度,cf為水動力阻尼系數,Ca為附加質量系數,本文設定Ca為1.0。

圖2為激勵力系數Cv云圖,源于圓柱體受迫振蕩試驗數據[15]。該受迫振蕩試驗測得的St數約為0.193。

圖2 渦激振動激勵力系數云圖Fig.2 VIV excitation force coefficient contour

1.2 鎖定判定準則

本文選取無因次頻率帶寬[0.125,0.25]作為鎖定區(qū)間,該鎖定區(qū)間需根據式(2)對真實環(huán)境及試驗條件下St數的差異進行修正。

(fD/VSt)test=(fD/VSt)actual

(2)

若響應頻率落在渦激振動激發(fā)帶寬內,鎖定發(fā)生,本文假定結構將被鎖定到渦激振動激勵中心對應的無因次頻率0.17,并根據式(2)進行修正,作為渦激振動主導激勵頻率。

當Cv為負時,水動力阻尼力作用于圓柱體,水動力阻尼系數可由式(3)計算。

(3)

當響應幅值和頻率超出試驗數據范圍時,采用Venugopal提出的阻尼模型[16],其中的經驗系數Chf,Csw及Clf推薦選取0.18,0.25及0.2。

1.3 結構動力響應求解

基于Euler-Bernoulli梁理論,橫流渦激振動動力學方程可由式(4)表示。

(4)

式中:m為結構單位長度質量,cs為結構阻尼系數cs=4πmfξ,ξ為結構阻尼比,E為彈性模量,I為慣性矩,T0為張力恒定成分,ΔT,PT和φ分別為時變張力的幅值、周期和初始相位。

圖3給出本文數值分析方法流程圖。

圖3 渦激振動時域分析流程圖Fig.3 Flow-chart of VIV time domain analysis

2 數值模型的試驗結果驗證

Franzini等針對長度L0=2.552 m,直徑D=0.022 m的立管模型開展了試驗研究,本節(jié)通過對比該試驗結果,驗證本文數值方法的有效性。試驗立管采用豎直放置形式,頂端施加40 N預張力,浸沒水中長度為2.257 m。模型頂部施加垂蕩幅值At/L0=1%的簡諧運動激勵。彎曲剛度為0.056 N·m2,質量比和浸沒重量分別為3.48和7.88 N/m。自由衰減試驗表明,該立管模型n階固有頻率fN,n為0.84nHz。本文采用有限元方法所得前三階固有頻率分別為0.85 Hz,1.70 Hz和2.57 Hz,與實測值吻合較好。

Franzini等給出了1個張力恒定工況和2個張力時變工況(ωT/ωCT=2和3)的試驗結果,其中ωT為時變張力圓頻率,ωCT為恒定張力時渦激振動響應圓頻率。本節(jié)對以上3組工況的試驗結果與數值結果進行對比,如圖4和圖5所示。

2.1 張力恒定工況

由圖4可知,當張力恒定時,渦激振動呈現出明顯的單一模態(tài)激發(fā)特征。預報所得位移幅值沿管長包絡線與試驗觀測總體較為吻合,僅在立管上半部分略大于實測值。兩者振幅峰值均出現在z/L=0.42附近,偏離立管中點。

監(jiān)測點(z/L0=0.43)處位移時歷服從正弦變化規(guī)律,預報和實測位移幅值均約為0.6D,幅值頻響譜與試驗結果吻合,該工況僅單一頻率被激發(fā),主導頻率接近于結構一階固有頻率。由此可知,約化速度VR,1=5.63時,渦激振動穩(wěn)定地激發(fā)該立管模型一階固有模態(tài)。

圖4 張力恒定工況結構響應對比Fig.4 Comparisons of structural response for constant tension case

2.2 張力時變工況

如圖5所示,張力時變工況下,立管渦激振動呈現明顯有別于張力恒定工況的響應特性,且受ωT/ωCT影響顯著。在ωT/ωCT=2工況中,幅值包絡線被顯著放大(約1倍)。Franzini等推斷這一放大現象是由于結構發(fā)生了由Mathieu不穩(wěn)定引起的參數共振。實測結果不再是標準的一階模態(tài)振型,尤其是在z/L位于[0.25,0.6]范圍內,一些高階模態(tài)明顯參與立管渦激振動。預報結果同樣顯示出高階模態(tài)被激發(fā)的響應特征,最大響應幅值略大于實測值。而對于ωT/ωCT=3工況,預報和實測響應幅值包絡線吻合較好,均展現出明顯的三階模態(tài)振型。預報值僅在z/L介于[0,0.3]內略大于試驗結果。預報與實測振幅峰值出現位置均偏離標準三階模態(tài)振型對應位置,呈現明顯的不對稱性。

對于豎直放置的立管模型,由于結構自重存在,張力沿管長方向由頂端向底端線性遞減,結構固有模態(tài)為滿足類Bessel函數的非正交解集。因此,即使在均勻流條件下,結構整體振型在理論上也會呈現不對稱性。針對張力沿軸向分布不均勻的工況,Da Silveira等借助自編程序及Orcaflex軟件均得到與圖5所示趨勢一致的數值結果,即結構振型具有行波特性。

就監(jiān)測點(z/L0=0.43)處立管局部結構響應而言,ωT/ωCT=2的時變張力激勵使得預報與實測位移幅值相較張力恒定工況放大約一倍。試驗觀測存在位移振幅劇烈調制的現象且調制過程無明顯的規(guī)律性。對應的預報結果中雖然也存在一定的振幅調制,但就整個時間歷程而言位移響應表現得更為穩(wěn)定。預報所得位移幅值頻響譜與實測結果保持較高一致性,能量響應集中在ω/ωCT=1處并伴隨著少量成分出現在ω/ωCT=2,3,4,5等處。實測結果中其他頻率成分弱于預報結果所示。而對于ωT/ωCT=3工況,預報和實測響應位移時歷同樣呈現明顯的振幅調制特性,幅值在0.6D附近輕微波動。預報和實測頻響譜在ω/ωCT=1和3處均出現突出的能量響應峰值,且ω/ωCT=3處的能量峰值大于ω/ωCT=1處。ω/ωCT=2,4,5,6的頻率成分同時出現在預報和實測頻響譜中,只是實測頻譜中響應能量相對微弱。針對該試驗立管,Franzini等[17]在純時變張力(ωT/ωCT=3)激勵的試驗工況中發(fā)現穩(wěn)定的三階橫向振動響應,而無第一、二階響應成分,即這一約化速度下第三階模態(tài)響應伴隨著第一階渦激振動同時發(fā)生。這也在一定程度上解釋了為何該工況體現的頻響特征與下文大尺度圓柱體模型設計工況有所區(qū)別。

(a)At/L0=1%,ωT/ωCT=2,VR,1=5.8

(b)At/L0=1%,ωT/ωCT=3,VR,1=5.63

通過與模型試驗對比顯示,本文采用的數值預報方法可較有效地用于時變張力工況下圓柱體渦激振動響應分析。

3 時變軸向張力條件下渦激振動響應特性研究

與大多數真實海洋環(huán)境下細長結構物相比,小尺度立管具有一定的局限性,最顯著的一點在于渦激振動僅激發(fā)結構一階模態(tài)。為深入研究時變張力下渦激振動響應特性,有必要挑選另一尺度較大的圓柱體模型,使得在低流速條件下也可激發(fā)結構高階模態(tài)。本文旨在研究時變張力這單一要素對渦激振動的影響,故選取一水平放置圓柱體模型作為研究對象。

Trim等[18]公布了Norwegian Deepwater Programme (NDP)采用的38 m長圓柱體模型,模型橫向放置在試驗水池中,表1給出模型的主要參數。

3.1 計算工況

38 m細長圓柱體模型在不同端部預緊力下前幾階固有頻率分布如圖6所示,可以發(fā)現預緊力對固有頻率具有顯著的影響,其本質原因在于預緊力直接影響結構剛度。

渦激振動發(fā)生時,結構將被鎖定到距離主導頻率最近的結構固有頻率。而主導頻率僅與St數、圓柱體直徑和來流速度有關。因此,即便在不同預緊力下渦激振動主導頻率均保持恒定,但激發(fā)模態(tài)階數卻未必不變。考慮張力時變效應時,可根據變化幅值的相對大小將工況分為兩類。變化幅值較小時激發(fā)模態(tài)階數保持不變,而變化幅值較大時激發(fā)模態(tài)則有可能隨時間發(fā)生階躍。

本文選取ΔT=500 N和4 000 N分別代表張力變幅較小和較大條件,14個時變張力頻率涵蓋高頻區(qū)及低頻區(qū),時變張力初始相位取φ=0。設計的28組計算工況參數見表2,ΔT/T0和ωT/ωCT分別為張力無因次變化幅值和頻率。流速均為0.4 m/s。

3.2 時變與恒定軸向張力結構響應對比

本節(jié)選取ΔT/T0=0.1,ωT/ωCT=1/20工況,對張力時變及恒定條件下渦激振動響應特性進行分析。采用連續(xù)小波變換(CWT)方法獲取圓柱體中點處(z=19 m)響應位移時頻分布云圖,母小波選取Morlet小波。位移幅值頻響譜根據快速傅里葉變換(FFT)得到。張力時變及恒定條件下圓柱體中點響應對比如圖7所示,其中前三行分別為端部張力、響應位移和響應頻率時歷,第四行為位移幅值頻響譜。

圖7 時變與恒定張力下結構響應對比(z=19 m)Fig.7 Response comparison with time-varying and constant tension (z=19 m)

由圖7可知,張力恒定時,渦激振動具有明顯的單模態(tài)激發(fā)特征,響應幅值和頻率在時間維度上保持不變。位移幅值頻響譜顯示,單一被激發(fā)的頻率成分為2.14 Hz。根據文獻[18],與之對應的試驗值為2.134 Hz,兩者吻合度較高。

當張力時變時,渦激振動幅值隨時間呈現規(guī)律性變化,其變化周期基本等于時變張力周期(即10 s)。伴隨著張力的變化,能量將由外界傳入渦激振動系統(tǒng)。根據能量守恒原理,響應位移趨向于變大以消散輸入的能量,而輸入能量大致與時變張力變化速度的平方成正比,故位移幅值變化應遵循張力變化速度絕對值的變化趨勢。因此,在一個變化周期內,響應位移按照類似于“山峰”的形狀調制兩次,對應張力變化的兩次加速和減速過程。如圖中淺藍色標記所示,最小幅值出現的時刻點并非嚴格對應于張力變化速度VT=0,而是略滯后于后者,即存在遲滯現象。

對于張力時變工況,響應頻率同樣展現出時變特性,隨張力的增大而變高,反之亦然。由圖6可知,端部預緊力位于4 500~5 500 N時,激發(fā)模態(tài)階數始終保持第三階,而結構固有頻率與張力存在正相關關系,因此響應頻率的變化與張力變化契合。根據位移幅值頻響譜,時變張力下渦激振動為多頻疊加響應,其中能量最強的峰頻仍為2.14 Hz,對應恒定張力工況下渦激振動主導頻率,同時激發(fā)了一些其他頻率成分,激發(fā)區(qū)間為[1.8,2.4]Hz,且主導頻率強度較張力恒定工況明顯變弱。

3.3 時變張力幅值效應

張力時變效應對渦激振動響應的影響由變化幅值和頻率共同決定,本節(jié)挑選ωT/ωCT=1和1/50代表時變張力頻率較高和較低兩種情況,研究時變張力幅值效應,圓柱體中點處結構響應對比如圖8和圖9所示,左右兩列分別對應ΔT/T0=0.1和0.8。

圖8 ωT/ωCT=1時結構響應對比(z=19 m)Fig.8 Response comparison with ωT/ωCT=1 (z=19 m)

圖9 ωT/ωCT=1/50時結構響應對比(z=19 m)Fig.9 Response comparison with ωT/ωCT=1/50 (z=19 m)

張力高頻變化條件下,ΔT/T0=0.1時,時頻分析及快速傅里葉變換結果均顯示結構動力響應特性與恒定張力工況大致相似,不同之處在于峰頻2.14 Hz的強度變弱,且額外出現0.14 Hz被激發(fā)的現象;而ΔT/T0=0.8時,結構響應變得相對復雜,位移時歷曲線呈現不穩(wěn)定性,響應頻率也隨著時間發(fā)生明顯的變化,能量最強的頻率成分為1.99 Hz,該頻率值非渦激振動主導頻率2.14 Hz,位移幅值顯著放大,超過恒定張力工況下對應值的兩倍,達到0.055 m,這一響應放大現象將在下文中著重討論。

張力低頻變化條件下,隨著ΔT/T0的增加,主激發(fā)頻率(2.14 Hz)的強度變弱,頻率激發(fā)帶寬相應變寬。由于更多亞諧振頻率成分參與渦激振動,ΔT/T0增大時響應位移和頻率時歷呈現出更為復雜的規(guī)律,結構響應具有明顯的周期性,遵循的周期與時變張力周期一致。

本節(jié)涉及的四個計算工況對應的響應位移時空分布如圖10所示,第一、二行分別對應ωT/ωCT=1和1/50,左右兩列對應ΔT/T0=0.1及0.8。張力高頻變化條件下,ΔT/T0較小時,結構響應幾乎與恒定張力工況一致;而當ΔT/T0較大時,整體振動形態(tài)仍大致滿足三階振型,但響應位移時空分布不具有明顯的規(guī)律。大幅值與小幅值振動交替出現,整個動力響應過程變得不再穩(wěn)定。張力低頻變化條件下,變化幅值效應呈現出新的特性。ΔT/T0較小時,結構在張力達到極值時被激發(fā)出有別于恒定張力工況的非標準三階振型;而ΔT/T0較大時,結構響應出現直觀的模態(tài)階躍現象。張力高頻變化條件下,由于周期太短,結構無法在不同模態(tài)間形成完整的轉換過程,故模態(tài)階躍并未出現。

圖10 不同ΔT/T0下響應位移時空分布Fig.10 Temporal-spatial distributions of displacement with different ΔT/T0

3.4 時變張力頻率效應

時變張力頻率對圓柱體動力響應的影響同樣顯著。本節(jié)選取ΔT/T0=0.1和0.8兩組工況分別討論張力變幅較小和較大時變化頻率對渦激振動響應的影響。結構中點處動力響應對比情況如圖11和圖12所示,由左至右分別對應ωT/ωCT=2,1/4和1/60。

小ΔT/T0條件下,當ωT/ωCT=2(即PT=0.25 s)時,結構響應并不遵循0.25 s的變化周期,與恒定張力工況相比,振幅調制現象明顯,位移顯著變大,極值達到0.03 m;當ωT/ωCT=1/4和1/60時,動力響應的周期性逐漸顯現,峰頻2.14 Hz始終占據主導,其他頻率成分同時被激發(fā),這些被額外激發(fā)的亞諧振頻率成分的具體數值嚴格遵循ωCT±kωT(k=1,2,3,……)的規(guī)律。

大ΔT/T0條件下,當ωT/ωCT=2時,響應位移幅值急劇增大,響應能量集中在1.99 Hz,而非一直保持的2.14 Hz,可認為此工況已誘發(fā)時變張力與渦激振動的結構聯合共振,詳見下文;當ωT/ωCT=1/60時,周期性重新出現,結構動力響應回歸穩(wěn)定狀態(tài);而對于ωT/ωCT=1/4工況,響應位移和頻率均顯得較為無序,位移幅值較小,各激發(fā)頻率間的強弱關系也并不固定,可認為其處于高頻共振區(qū)與低頻穩(wěn)定區(qū)之間的過渡區(qū)域。ΔT/T0較大時,同樣存在著ωCT±kωT的頻率成分可能被激發(fā)的情況。由圖10可知,時變張力頻率對渦激振動響應的另一影響效應體現在高頻變化抑制模態(tài)階躍的發(fā)生。

圖11 ΔT/T0=0.1時結構響應對比(z=19 m)Fig.11 Response comparison with ΔT/T0=0.1 (z=19 m)

圖12 ΔT/T0=0.8時結構響應對比(z=19 m)Fig.12 Response comparison with ΔT/T0=0.8 (z=19 m)

4 聯合激勵共振討論

時變張力和渦激振動聯合激勵下,結構可能發(fā)生共振,響應位移幅值將顯著增大。出于結構安全考慮,共振狀態(tài)是海洋細長結構物在工程設計階段希望極力避免的,因此本節(jié)將就這一可能出現的共振現象展開討論。

圖13給出各工況下最大位移分布,圖中藍色實線標注出恒定張力工況對應值。當ωT/ωCT大于1/2時(高頻區(qū)),最大響應位移普遍大于恒定張力工況,且極值出現在ωT/ωCT=2處,與文獻[2]中試驗結果一致。由此可推斷出ωT=2ωCT是時變張力與渦激振動聯合激勵下結構最易發(fā)生共振的核心頻率。而ωT=2ωCT正是觸發(fā)結構Mathieu失穩(wěn)的最危險條件,這種由于Mathieu失穩(wěn)造成的結構共振伴隨著渦激振動鎖定現象的發(fā)生被進一步放大。因此,本文認為此類聯合激勵下發(fā)生的共振可定義為Mathieu型VIV共振。當ωT/ωCT處于高頻區(qū)某一固定值時,響應位移的放大效應隨著ΔT/T0的增大而越發(fā)明顯,這一規(guī)律可以基于能量守恒原理給出相應解釋。當ωT相同時,時變張力向渦激振動系統(tǒng)輸入的能量與ΔT/T0正相關,故ΔT/T0較大時,響應位移趨向于增大更多以消除更多輸入的能量。圖中A/D>20的結果超出普遍認知,表明在極端工況條件下Mathieu型共振對結構響應的放大效應極為顯著。文獻[9]基于尾流振子模型研究此類極端工況時也得出與圖13所示數量級相當的預報結果。

圖13 不同工況下響應位移最大值Fig.13 Maximum displacements under different cases

隨著ωT由共振中心2ωCT向兩側偏離,最大位移也將逐漸減小并趨向于張力恒定工況對應值。ωT/ωCT小于1/2時的結果在圖13左側的紅色方框內被局部放大。由于ΔT/T0較大時,高頻共振區(qū)與低頻穩(wěn)定區(qū)之間存在著結構響應較為混亂的過渡區(qū)域,故ωT/ωCT位于[0.05,0.25]時,ΔT/T0=0.8時的最大響應位移小于ΔT/T0=0.1對應值。如圖8和圖12所示,Mathieu型共振發(fā)生時,主導結構響應的激發(fā)頻率不再為渦激振動鎖定頻率2.14 Hz,而是偏移至其附近的一個新值,約為2 Hz,因此可認為Mathieu型共振誘發(fā)了失諧現象。

5 結 論

本文基于渦激振動流體力分解模型,對其預報軸向張力恒定和時變工況的有效性進行了試驗結果驗證,并以某38 m圓柱體模型為研究對象,討論了時變張力對渦激振動響應的影響,得到如下結論:

(1)張力時變條件下結構渦激振動響應具有明顯的分時特征,呈現出振幅調制、遲滯、頻率轉換及多模態(tài)響應疊加等有別于張力恒定工況的新特性。

(2)隨著無因次時變張力幅值ΔT/T0的增大,激發(fā)頻率帶寬變寬,主導頻率強度變弱。由于更多的亞諧振頻率成分參與,響應位移及頻率時歷表現出更加復雜的規(guī)律性。張力低頻變化條件下,大ΔT/T0工況中結構響應存在模態(tài)階躍現象。

(3)考慮張力時變效應時,除渦激振動鎖定頻率外,對應ωCT±kωT的亞諧振頻率成分可能被激發(fā)。張力高頻變化時,結構響應與張力恒定工況相似;而對于張力低頻變化的情況,結構響應具有與時變張力周期對應的周期性。較高的時變張力頻率會抑制結構響應在各模態(tài)間的相互轉換。

(4)時變張力與渦激振動聯合激勵下,結構會在時變張力頻率ωT=2ωCT附近發(fā)生Mathieu型VIV共振,響應位移幅值將被明顯放大,且放大程度與ΔT/T0正相關。Mathieu型共振可能誘發(fā)結構響應失諧,結構將被鎖定到一個新的激發(fā)頻率而非渦激振動主導頻率。

本文研究結論有助于加深對復雜海洋環(huán)境下海洋細長結構物渦激振動問題的理解和認識,在特定條件下誘發(fā)的Mathieu型VIV共振現象需要引起足夠重視。

猜你喜歡
模態(tài)振動結構
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
中立型Emden-Fowler微分方程的振動性
論《日出》的結構
國內多模態(tài)教學研究回顧與展望
創(chuàng)新治理結構促進中小企業(yè)持續(xù)成長
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
UF6振動激發(fā)態(tài)分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 久久久精品国产SM调教网站| 久久这里只有精品免费| 中文字幕无码av专区久久| 国产精品深爱在线| 国产美女91呻吟求| 亚洲国产精品人久久电影| 在线精品视频成人网| 狠狠色丁香婷婷综合| 国产精品9| 亚洲综合片| 久久男人视频| 波多野结衣无码中文字幕在线观看一区二区 | 国产在线精彩视频二区| 亚洲欧洲自拍拍偷午夜色| 久久精品亚洲热综合一区二区| 免费一级毛片在线观看| 四虎在线高清无码| 香蕉视频在线精品| 欧美人人干| 毛片一级在线| 免费一级无码在线网站| 久久五月视频| 色呦呦手机在线精品| swag国产精品| 亚洲人成成无码网WWW| 欧美精品成人| 久久久精品久久久久三级| 精品无码一区二区三区在线视频| 亚洲精品视频免费观看| 在线观看亚洲精品福利片| 狠狠干综合| 无码aⅴ精品一区二区三区| 第九色区aⅴ天堂久久香| 都市激情亚洲综合久久| 国产精品欧美在线观看| 波多野结衣无码AV在线| 三区在线视频| 亚洲欧美日韩天堂| 成人无码一区二区三区视频在线观看 | 亚洲国产清纯| 亚洲成a人在线播放www| 中文字幕久久亚洲一区| 综合久久五月天| 98超碰在线观看| 欧美精品伊人久久| 熟女成人国产精品视频| 一区二区三区国产| 成人免费一级片| 丁香六月综合网| 无码中文AⅤ在线观看| 日韩精品毛片人妻AV不卡| h视频在线播放| 亚洲三级电影在线播放| 国产Av无码精品色午夜| 91精品专区国产盗摄| 欧美国产精品不卡在线观看| 97在线观看视频免费| 国产jizz| 又黄又爽视频好爽视频| 国产91色在线| 在线观看视频99| 亚洲欧州色色免费AV| 2021国产乱人伦在线播放| 九九这里只有精品视频| 精品在线免费播放| 国产熟女一级毛片| 亚洲愉拍一区二区精品| 自偷自拍三级全三级视频| 午夜小视频在线| 欧美成人影院亚洲综合图| 国产精品第一区| 久久精品欧美一区二区| 成人毛片在线播放| 日韩小视频网站hq| 亚洲天堂视频在线观看免费| 亚洲国产成人精品一二区| 老司机久久精品视频| 久久综合伊人77777| 五月天久久婷婷| 国产区人妖精品人妖精品视频| 久久一色本道亚洲| 亚洲—日韩aV在线|