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

兩自由度平板大幅運動的氣動特性與穩定性的CFD研究

2021-06-10 00:29:09祝志文顏爽王欽華李加武
振動工程學報 2021年2期

祝志文 顏爽 王欽華 李加武

摘要: 為研究兩自由度薄平板大幅運動的氣動特性,評價其氣動穩定性,基于任意拉格朗日?歐拉描述法的動網格技術,通過有限差分法求解描述任意流變區域不可壓流動的控制方程,開展了不同折減風速下薄平板豎彎和扭轉運動繞流場的CFD (Computational Fluid Dynamics)模擬。研究認為,單自由度薄平板小幅運動的氣動力系統是線性和穩定的,即使單自由度大幅豎彎運動也是線性和氣動穩定的。但單自由度大幅扭轉運動的平板氣動力系統出現非線性,并隨折減風速的提高非線性變得顯著,且平板將進入氣動不穩定狀態。另外,大幅扭轉耦合不同豎彎振幅運動的平板,氣動力系統均為非線性并隨折減風速的提高越加顯著,而該非線性主要來自扭轉自由度的大幅運動;對該兩自由度耦合系統,當豎彎振幅較小和折減風速較高時,氣動力系統是不穩定的;但當豎彎振幅較大時,氣動力系統將是穩定的。

關鍵詞: 氣動穩定性; 薄平板; 大幅運動; CFD; 氣動力非線性

中圖分類號: U448.21+3; V211.3 ? ?文獻標志碼: A ? ?文章編號: 1004-4523(2021)02-0271-12

DOI:10.16385/j.cnki.issn.1004-4523.2021.02.007

引 ?言

橋梁跨徑的不斷增大使得結構的頻率和阻尼不斷降低,導致橋梁對自然風作用的敏感程度明顯增加,因而可能引起橋梁主梁的大幅運動。這些大幅運動,包括主梁的大幅渦激運動[1]和顫振[2],以及拉索的大幅振動[3?4]。對鈍體外形的橋梁主梁,或主梁因大幅運動產生較大的相對攻角效應,橋梁主梁的氣動力系統可能會因這種大幅運動而呈現顯著的非線性[5?6]。開展橋梁主梁大幅豎彎和扭轉運動下的氣動特性研究,有助于評價現有線性氣彈理論的應用條件和不足,有利于揭示橋梁大幅運動氣動力系統的非線性特征和影響因素,以及橋梁氣動穩定性的影響因素[7]。

在風洞中開展橋梁主梁大幅豎彎和扭轉運動試驗,存在試驗驅動機構設計和布置、測力和測壓、模型保護等諸多困難,而基于CFD開展這一研究,不僅能克服上述困難,還能方便地設置運動參數和邊界條件,具有風洞試驗無法相比的優勢。當基于CFD開展橋梁主梁的大幅運動模擬時,因主梁周圍空氣所占據的實際空間隨時間在發生顯著的變化,具有大幅運動的邊界,對此類顯著流變區域流體的描述,需采用任意拉格朗日?歐拉描述法[8]。下面以橋梁氣彈研究常采用的薄平板為研究對象,研究其大幅運動的氣動特征和氣動穩定性。

1 薄平板強迫運動繞流場控制方程

基于二維流體運動的任意拉格朗日?歐拉描述方程[8],通過建立計算域和物理域的隨時間變化的映射關系,可在不隨時間變化的二維計算域上,采用歐拉描述法實現對二維物理域上任意時變區域的流體流動的描述。根據橋梁主梁剛性斷面的事實,以及主梁豎彎和扭轉運動的特點,可簡化得到薄平板強迫運動繞流場控制方程。

如圖1所示為大幅運動的薄平板,設OXY表示絕對參考系下的直角坐標系,并用oxy表示隨平板一起做剛性運動的動參考系下的直角坐標系,平板具有繞平衡位置的豎向自由度和繞轉動軸的扭轉自由度。設平板豎向運動速度為V0,繞通過其中點垂直于OXY平面軸的扭轉角速度為ω(t)=dψ(t)/dt。對薄平板上述給定的兩自由度運動,其對應的任意拉格朗日?歐拉描述方程可簡化為動參考系下關于原參變量u,v和p的N?S方程[9],也即:

3.2 平板豎彎和扭轉運動CFD模擬方法驗證

為驗證薄平板豎彎或扭轉運動CFD模擬數值方法的合理性,本文以非定常氣動力模擬為依據,通過有解析解的薄平板顫振導數Theodorsen解析解與本文識別結果的對比,驗證本文方法的合理性。對豎彎或扭轉運動的二自由度運動薄平板,其受到的非定常氣動力作用可用八個顫振導數來表達[11?12],即:

式中 ?k=2πfb/U_0為折減頻率;H_i^*和A_i^* (i=1?4)是顫振導數。

在不同折減風速下開展CFD計算,將分別得到作用在豎彎和扭轉運動薄平板上的氣動力,利用式(18)和(19)并通過最小二乘法可識別平板對應折減風速下的八個氣動導數。本文CFD模擬首先開展了網格無關性檢查,檢查對象是顫振導數中相對重要的A_2^*導數在折減風速為4,8和12上的值,物面第一個網格高度分別取0.0005b,0.001b和0.002b三種情況,將計算結果與Theodorsen解析解對比表明,在0.0005b和0.001b兩種情況下得到的結果與解析解一致性最好,也好于0.002b的情況。從非定常計算較大的計算量上考慮,本文CFD計算平板物面第一個網格高度取0.001b。

如圖4所示,將本文識別的顫振導數和平板Theodorsen解析解作了比較。可見在低折減風速下,二者吻合很好;隨著折減風速的增大,二者出現較小偏差。偏差的原因可能是解析解是基于無黏、無流動分離的勢流理論,但本文是真實的有黏空氣流動,且平板運動會導致流動分離,而黏性效應和流動分離可能隨折減風速的增大而影響顯著。本文識別結果與Theodorsen解析解[13]趨勢性非常一致,這說明本文數值模擬能有效地獲得運動薄平板上作用的氣動力,因而能用于研究運動薄平板的氣動力和能量特征。

3.3 單自由度小幅運動

為便于對比平板大幅運動的氣動力和能量特征,下面給出Vr=8下,平板分別作小幅豎彎(振幅為B/80)或小幅扭轉運動(振幅0.2°)的氣動力時程。圖5(a)和(b)分別是從多個氣動力周期中取出一個完整周期(每個周期氣動力完全相同),對應平板單自由度豎彎或單自由度扭轉運動氣動力系數時程曲線,可見無論是升力還是扭矩系數,其曲線具有明顯的諧波特性,其時程頻率等于扭轉強迫運動頻率。因此當平板單自由度小幅運動時,其氣動力系統具有良好的線性特性。

圖6為薄平板單自由度小幅扭轉運動在一個周期內4個關鍵時刻的風壓拓撲圖,對應從平衡位置出發后達到最大正攻角位置(圖6(a)),然后回到平衡位置(圖6(b)),再達到最大負攻角位置(圖6(c)),最后回到零攻角的平衡位置(圖6(d))。可見從一個狀態到另外一個狀態,壓力呈現規律性變化,反映了平板扭轉運動方向的改變。平板前緣為風壓值大和梯度高的區域,后緣為低風壓區且壓力變化遠不及前緣劇烈。

圖7(a)和(b)是在不同折減風速下豎彎運動的平板,一個周期內氣流分別在平板上表面各點、以及整個平板上所做的無量綱總功。可見氣流對平板表面所有點均做負功,也即平板表面均消耗氣流能量,前緣(L.E.)局部耗能最強,后緣(T.E.)耗能最弱。圖7(b)給出了Vr=8?60,豎彎運動平板一個周期內氣流作用總功隨折減風速的變化。從能量的角度可見,單自由度小幅豎彎運動的平板是氣動穩定的。

圖8(a)和(b)分別為單自由度小幅扭轉運動的平板,一個周期內氣流在平板上表面各點、以及整個平板所做的無量綱總功。從圖8(a)可見平板前緣為穩定區,其后到轉動軸(A.R)為不穩定區,也即氣流對平板做正功,平板從氣流中吸收能量,因此為激勵區;轉動軸到平板后緣氣流做負功,因而是耗散氣流能量的穩定區,此時平板上激勵區與穩定區共存。圖8(b)給出了Vr=8?60,一個扭轉周期內氣流對平板所做的總功,可見在所計算的折減風速范圍內,一個周期內氣流對扭轉運動的平板均做負功,也即平板消耗氣流的能量,因此單自由度小幅扭轉運動的平板也是氣動穩定的。

3.4 單自由度大幅運動

為研究薄平板單自由度大幅運動的氣動特性,需對豎彎或扭轉自由度方向設定遠大于上述單自由度小幅運動的幅值。根據文獻[7]對橋梁斷面氣動力非線性特性的研究結果,本文研究中,取大幅豎彎振幅為B/4,大幅扭轉振幅為20°,以研究其氣動力的線性特征和氣動穩定性。

3.4.1 單自由度大幅豎彎運動(LAH)

圖9(a)和(b)分別給出Vr=12和Vr=40時,分別從多個氣動力周期時程中,取出一個完整周期(每個周期氣動力完全相同),對應單自由度大幅豎彎運動的氣動力時程。可見雖然折減風速差別大,但升力和扭矩系數時程的相位差變化不明顯,特別是升力和扭矩系數時程曲線均具有較好的諧波特性,也即沒有觀察到明顯的氣動力非線性,因此可以認為,大幅豎彎運動的薄平板氣動力系統的非線性不明顯。

圖10(a)和(b)分別是在一個周期內,氣流在大幅豎彎運動平板上表面各點、以及對整個平板所做的無量綱總功。可見氣流對豎彎運動的平板表面所有點均做負功,且隨著折減風速的提高,功的絕對值雖越來越小,但仍為負。同時,在全部折減風速上,氣流對整個平板所做的無量綱總功也均為負。因此,在本文所考慮的較大折減風速范圍內,一個周期內氣流對大幅度豎彎運動的平板均做負功,也即消耗氣流的能量,因此單自由度大幅豎彎運動的平板是氣動穩定的。

3.4.2 單自由度大幅扭轉運動(LAP)

考察大幅扭轉運動的情況,圖11(a)和(b)分別為折減風速Vr=8和Vr=40時,分別從多個氣動力周期時程中,取出一個完整周期(每個周期氣動力完全相同),該周期內的氣動力系數時程曲線。與圖9相比,可見隨著折減風速的提高,升力和扭矩系數間相位逐漸變化,且時程曲線波形發生畸變,也即失去諧波特性,特別是高折減風速下氣動力時程曲線的畸變更顯著。對低和高折減風速下氣動扭矩系數時程的歸一化功率譜分析表明,在低折減風速下扭矩系數時程的頻率等于扭轉運動頻率,如圖12(a)所示;但高折減風速下扭矩系數時程雖然峰值頻率為扭轉運動頻率,但頻譜圖中明顯出現了倍頻成分,如圖12(b)所示,且倍頻成分的能量明顯小于扭轉運動頻率對應峰值能量,可見出現的非線性不是特別顯著。因此,當薄平板大幅扭轉運動時,其氣動力系統不再為線性,且隨著折減風速的提高,其非線性將越來越顯著。這可能是折減風速增大,繞薄平板流動的分離加劇,因此平板表面的流態隨著折減風速的提高在不斷改變,使得氣動力的非線性變得明顯。

圖13給出了一個扭轉運動周期內氣流對平板所做無量綱總功隨折減風速的變化,可見當折減風速小于16時,一個周期內氣流對平板所做總功為負,也即平板消耗氣流的能量;但當折減風速大于20后,氣流對上游側半個平板估的正功將大于對下游側半個平板所做負功,氣流對平板所做的總功由負變為正,此時一個周期內扭轉運動平板將從氣流中吸收能量,表明大幅扭轉運動平板已進入氣動不穩定狀態。

3.5 兩自由度耦合運動

從上面研究可知,當扭轉運動幅值角為20°且折減風速較高后,單自由度大幅扭轉運動的薄平板將進入氣動不穩定狀態。下面針對扭轉運動幅值角20°的情況,考察不同豎彎運動幅值平板兩自由度耦合系統的氣動特征,豎彎幅值分別采用小幅、中幅和大幅,對應的幅值分別是B/50,B/16和B/4,并設定耦合運動的豎彎和扭轉自由度相位角為0o(CFD可設定任意相位角)。分別開展不同折減風速的CFD模擬,比較作用在平板上的氣動力和一周期內氣流對平板做功的能量特征。

3.5.1 小幅豎彎大幅扭轉(SAH?LAP)

圖14(a)和(b)分別是折減風速Vr=10和Vr=40時,從各自多個氣動力周期時程中,取出一個完整周期(每個周期氣動力完全相同),對應小幅豎彎大幅扭轉運動平板一個周期內的氣動力系數時程曲線。可見耦合運動平板的氣動力系統表現出非線性,且隨著折減風速的提高,其非線性特性越來越顯著,其特征與圖11相似。

圖15給出了兩個折減風速下耦合運動薄平板,在一個周期內氣流在運動平板上表面各點所做的無量綱總功。可見每個自由度方向氣流作用的能量特征與單自由度類似,也即豎彎自由度方向平板表面點均消耗氣流能量,但扭轉自由度方向是迎風側上半個平板吸收氣流能量,下游側半個平板消耗氣流能量。當折減風速較高時,上游側半個平板吸收的能量明顯大于下游側半個平板消耗的能量,導致扭轉自由度方向吸收氣流能量。如綜合兩個自由度方向的總能量,在低折減風速下豎彎自由度消耗氣流的能量大于扭轉自由度吸收的能量,因此平板是氣動穩定的;但當折減風速較高后,豎彎自由度消耗氣流的能量小于扭轉自由度方向吸收的能量,此時兩自由度耦合運動系統是氣動不穩定的,如圖16所示。從圖15也可見,由于豎彎自由度振幅小,其消耗的氣流能量小,因此氣流在扭轉自由度方向作用的能量主導了兩自由度耦合運動系統氣流作用的總能量。

3.5.2 中幅豎彎大幅扭轉(MAH?LAP)

維持薄平板兩自由度運動的扭轉幅值角為20°,并將豎彎運動位移幅值增大到B/16。圖17給出了折減風速分別是Vr=10和Vr=40時,分別從多個氣動力周期時程中,取出一個完整周期(每個周期氣動力完全相同),對應一個周期內薄平板的升力和扭矩系數時程。可見即使低折減風速,薄平板氣動力的諧波特性變差,而隨著折減風速的提高,氣動力時程的非線性特性越加顯著。

圖18給出了高低兩個折減風速下耦合運動的薄平板,一個周期內氣流在運動平板上表面各點所做的無量綱總功。其曲線特征與圖15相似,也即豎彎自由度方向平板表面所有點均消耗氣流能量,但扭轉自由度方向是迎風側上半個平板吸收氣流能量,下游側半個平板消耗氣流能量。但由于豎彎自由度方向的振幅增大,豎彎運動消耗氣流的能量明顯增大,從兩個自由度方向的能量總和來看,在所有折減風速下,平板均消耗氣流的能量,如圖19所示。因此,以中幅豎彎大幅扭轉運動的薄平板是氣動穩定的。與前述小幅豎彎大幅扭轉運動情況相比,可見豎彎幅度的增大提高了平板的氣動穩定性。

3.5.3 兩自由度耦合大幅運動(LAH?LAP)

進一步增大豎彎運動幅值到B/4,使得兩個自由度方向均為大幅運動,也即大幅耦合運動。圖20給出了折減風速分別是Vr=10和Vr=40時,分別從多個氣動力周期時程中,取出一個完整周期(每個周期氣動力完全相同),對應薄平板一個周期內大幅耦合運動的升力和扭矩系數時程。可見氣動力的非線性明顯,也隨折減風速的增大而顯著。與圖17相比,相同的折減風速下氣動力的非線性并沒有明顯變化。

圖21為低折減風速(Vr=10)和高折減風速(Vr=40)下大幅耦合運動的薄平板,一個運動周期內氣流在運動平板上表面各點所做的無量綱總功。豎彎和扭轉運動氣流作用能量沿平板表面的分布與圖18稍有不同,可能是兩自由度大幅耦合運動導致。但因豎彎自由度方向的大幅運動,導致豎彎自由度方向消耗氣流的能量大幅度增加,使得一個周期內氣流對平板表面所有點均做負功,因此氣流作用在運動平板上的總能量為負,并明顯大于圖16小幅豎彎大幅扭轉和圖19中幅豎彎大幅扭轉運動的情況,如圖22所示,因此大幅耦合運動的平板是氣動穩定的,原因是豎彎自由度方向的振幅增大,顯著提高了其氣動穩定性。

4 結 ?論

本文基于任意拉格朗日?歐拉描述法,開展了薄平板豎彎和扭轉強迫運動繞流場CFD模擬,研究了薄平板大幅運動的氣動特征和穩定性,得到下述結論:

1) 單自由度小幅運動的薄平板,其氣動力系統是線性的,也是氣動穩定的。

2) 單自由度大幅豎彎運動薄平板是氣動穩定的,且氣動力系統非線性不明顯;但單自由度大幅扭轉運動薄平板氣動力系統是非線性的,其非線性隨折減風速提高而顯著,且折減風速較高后,其氣動力系統將變得不穩定。

3) 大幅扭轉的耦合運動平板,氣動力系統均為非線性,其非線性主要來自扭轉自由度,非線性隨折減風速的提高更顯著。

4) 大幅扭轉的耦合運動平板氣動穩定性,取決于豎彎自由度消耗的氣流能量與扭轉自由度吸收氣流能量的相對大小。當豎彎振幅較小時和折減風速較高時,氣動力系統是不穩定的;當豎彎自由度方向振幅較大時,其氣動力系統將是穩定的。

參考文獻:

[1] Zhu Q, Xu Y L, Zhu L D, et al. Vortex-induced vibration analysis of long-span bridges with twin-box decks under non-uniformly distributed turbulent winds[J]. Journal of Wind Engineering & Industrial Aerodynamics, 2018,172: 31-41.

[2] Malík Josef. Sudden lateral asymmetry and torsional oscillations in the riginal Tacoma suspension bridge[J]. Journal of Sound and Vibration, 2013,332(15):3772-3789.

[3] 祝志文.超臨界雷諾數下拉索順風向自激力特性研究[J].振動工程學報,2014, 27(3):377-384.

Zhu Zhiwen. Characteristics of along wind self-excited forces of a cable under super-critical Reynolds number[J]. Journal of Vibration Engineering, 2014, 27(3):377-384.

[4] 祝志文,陳 魏,李健朋,等.多塔斜拉橋加勁索渦激振動實測與時域解析模態分解研究[J].中國公路學報,2019,32(10):247-256.

Zhu Zhiwen. CHEN Wei, LI Jianpeng,et al. Field observation of vortex-induced vibration of stiffening cables in a multi-tower cable-stayed bridge with application of analytical mode decomposition[J]. China Journal of Highway and Transport, 2019,32(10):247-256.

[5] Wu T, Kareem A. A nonlinear analysis framework for bluff-body aerodynamics: A Volterra representation of the solution of Navier-Stokes equations[J]. Journal of Fluids and Structures, 2015,54:479-502.

[6] 祝志文,石亞光,李健朋.扁平箱梁基于Volterra理論的氣動力非線性特性研究[J].中國公路學報,2018, 31 (1): 74-81.

Zhu Zhiwen, Shi Yaguang, Li Jianpeng. Research on nonlinear aerodynamics of flat steel box girder based on Volterra theory[J]. China Journal of Highway and Transport , 2018, 31 (1): 74-81.

[7] 祝志文,顧 明,陳政清.兩自由度系統氣彈響應特性研究[J].工程力學,2008,25(8):23-30.

Zhu Zhiwen, Gu Ming, Chen Zhengqing. Study on response characteristics of two-DOF aeroelastic system[J]. Engineering Mechanics, 2008, 25(8):23-30.

[8] Filipovic Nenad, Mijailovic Srboljub, Tsuda Akira, et al. An implicit algorithm within the arbitrary Lagrangian-Eulerian formulation for solving incompressible fluid flow with large boundary oscillations[J]. Computer Methods in Applied Mechanics and Engineering, 2006,195(44-47): 6347-6361.

[9] Zhu Zhiwen, Chen Zhengqing, Gu Ming. CFD based simulations of flutter characteristics of ideal thin plates with and without central slot[J]. Wind and Structures, 2009, 12(1): 1-19.

[10] Tau E Y. A second-order Projection method for incompressible Navier-Stokes equations in arbitrary domains[J]. Journal of Computational Physics, 1994, 115:147-152.

[11] 祝志文,陳政清.數值模擬橋梁斷面的顫振導數和顫振臨界風速[J].中國公路學報,2004,15(4):41-50.

Zhu Zhiwen, Chen Zhengqing. Numerical simulations for aerodynamic derivatives and critical flutter velocity of bridge deck[J]. China Journal of Highway and Transport, 2004, 15(4):41-50.

[12] 祝志文,顧 明,陳政清.基于3211多階躍激勵CFD模型的顫振導數識別研究[J].振動工程學報,2008, 21(1):18-23.

Zhu Zhiwen, Gu Ming, Chen Zhengqing. Identification of flutter derivatives based on 3211 multi-step excitation of CFD model[J]. Journal of Vibration Engineering, 2008,21(1):18-23.

[13] Scanlan R H, Tomko J J. Airfoil and bridge deck flutter derivatives[J]. Journal of Engineering Mechanics, ASCE, 1971,97(6): 1171-1737.

CFD investigation on aerodynamics and stability of a two degree-of-freedom thin plate undergoing large-amplitude oscillation

ZHU Zhi-wen1, YAN Shuang1, WANG Qin-hua1, LI Jia-wu2

(1. Key Laboratory of Structure and Wind Tunnel of Guangdong Higher Education Institutes (Shantou University),

Shantou 515063, China; 2. Highway School, Chang'an University, Xi'an 710064, China)

Abstract: In order to investigate the aerodynamic characteristics and aerodynamic stability of the two degree-of-freedom (2-DOF) thin plate undergoing large-amplitude oscillation, the governing equations for arbitrary-deformed fluid domain in the two-dimensional incompressible form are numerical solved based on the Arbitrary-Lagrangian-Eulerian description and the finite difference method in moving-grid system. Then, computational fluid dynamics (CFD) is employed to simulate the flow field around the thin plate undergoing heave and/or pitch oscillation at various reduced wind speeds. The research finds that the aerodynamic system around the thin plate is linear and stability when the thin plate experiences small-amplitude vibration in heave or pitch, and it is also linear and stable even when the thin plate undergoes large-amplitude heave motion. However, when the thin plate undergoes large-amplitude pitch motion, the aerodynamic system of the thin plate presents nonlinearity which will become significant and even turn into instability with the increase of reduced wind speed. Meanwhile, when the thin plate experiences 2-DOF oscillation with large amplitude in pitch coupled with different amplitude in heave, nonlinearity will also be presented and be more significant with the increase of reduced wind speed. It is noted that such nonlinearity is generated from the large-amplitude oscillation in pitch. For the coupled vibration system with 2-DOF, when the vibration amplitude is small while the reduced wind speed is high, its aerodynamic system is instable. While when the vibration amplitude in heave is large, the aeroelastic system will be stable.

Key words: aerodynamics stability; thin plate; large-amplitude oscillation; CFD; aerodynamic nonlinearity

作者簡介: 祝志文(1968-),男,教授,博士生導師。電話: 13574876655; E-mail: zhuzw@shu.edu.cn

主站蜘蛛池模板: 91网在线| 国产在线八区| 色135综合网| 欧美日韩国产精品综合 | 欧洲亚洲一区| 人妖无码第一页| 国产杨幂丝袜av在线播放| 欧美成人免费一区在线播放| 精品少妇人妻一区二区| 蜜芽一区二区国产精品| 九九精品在线观看| 亚洲成在线观看 | 狠狠色狠狠色综合久久第一次| 欧美一区日韩一区中文字幕页| 中文字幕天无码久久精品视频免费| 一级不卡毛片| 午夜视频免费一区二区在线看| 99re66精品视频在线观看| 操操操综合网| 国产女人18毛片水真多1| 日韩成人午夜| 欧美日韩中文字幕在线| 欧美色视频在线| 亚洲天堂免费在线视频| 免费不卡在线观看av| 女人18毛片水真多国产| 91蝌蚪视频在线观看| 国产三级精品三级在线观看| 亚洲欧美精品在线| 久久九九热视频| 亚洲成人黄色在线| 欧美亚洲国产精品久久蜜芽| 国产超碰一区二区三区| 91在线一9|永久视频在线| 国产又爽又黄无遮挡免费观看 | 亚洲清纯自偷自拍另类专区| 婷婷色狠狠干| 国产欧美在线观看精品一区污| 色偷偷综合网| 伊人国产无码高清视频| 在线播放精品一区二区啪视频| 欧洲成人免费视频| 人人爽人人爽人人片| 成人夜夜嗨| 亚洲侵犯无码网址在线观看| 国产成人精品在线| 91久久精品日日躁夜夜躁欧美| 成人福利在线视频| 91国语视频| 亚洲日本中文字幕乱码中文| 99精品伊人久久久大香线蕉| 久久精品国产999大香线焦| 免费国产无遮挡又黄又爽| 久久人人妻人人爽人人卡片av| 在线一级毛片| 久久精品欧美一区二区| 九色视频在线免费观看| 久久亚洲国产一区二区| 国产精品美女网站| 一级毛片高清| 无码一区中文字幕| 一级做a爰片久久毛片毛片| 欧美精品伊人久久| 久青草网站| 亚洲va欧美va国产综合下载| 在线国产91| 亚洲一区二区三区在线视频| 久久综合亚洲鲁鲁九月天| 亚洲三级a| 亚洲人成网7777777国产| 青青青亚洲精品国产| 国产日本一区二区三区| 欧美日韩国产成人高清视频| 亚洲全网成人资源在线观看| 久久久久无码国产精品不卡| 国产极品美女在线播放| 亚洲综合色吧| 国产成人精品一区二区三区| 2021国产乱人伦在线播放| 91麻豆国产精品91久久久| 免费看美女自慰的网站| 久久午夜夜伦鲁鲁片不卡 |