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

內(nèi)激勵型振蕩衰減流作用下輸流管道動力不穩(wěn)定分析

2021-02-06 11:25:06林震寰
振動與沖擊 2021年3期
關(guān)鍵詞:區(qū)域

張 挺,林震寰,林 通,張 恒

(1.福州大學(xué) 土木工程學(xué)院,福州 350116; 2.中國電建集團(tuán)貴陽勘測設(shè)計研究院有限公司,貴陽 550081)

輸流管道在工程中有著廣泛運用,當(dāng)管道內(nèi)的流速隨時間變化時,將產(chǎn)生脈動流,其中周期性脈動流最為常見。作為一種內(nèi)激勵的形式,會影響輸流管道的穩(wěn)定性,從而引起了很多學(xué)者和工程界的關(guān)注。

Chen[1]首次將脈動流表達(dá)式引入輸流管道橫向振動運動方程,采用HSU的方法和Bolotin的方法分別研究其作用下的動力不穩(wěn)定問題,之后Paidoussis等[2-3]對其作了進(jìn)一步地研究,證明了一定頻率和振幅的內(nèi)流脈動使原來不穩(wěn)定的輸流直管可能變得穩(wěn)定。Bolotin法常被作為結(jié)構(gòu)物在周期性或瞬態(tài)激勵下動力不穩(wěn)定區(qū)域的分析方法,例如:王杰方等[4]采用伽遼金變分法和Bolotin法,研究超空泡運動體的動力不穩(wěn)定區(qū)域和影響參數(shù)共振曲線的因素;Arani等[5]基于正弦剪切變形梁理論,采用微分求積法結(jié)合Bolotin法,研究得到脈動流作用下的雙壁納米碳管動力不穩(wěn)定區(qū)域。隨著對輸流管道認(rèn)知和數(shù)值方法的不斷深入及發(fā)展,Gorman等[6]、Azrar等[7]和Seo等[8]分別采用有限差分法、微分求積法、有限元法對內(nèi)流作用下輸流管道參數(shù)共振特性及動力不穩(wěn)定區(qū)域作了深入研究。Wang等[9-10]采用彈性伯努利-歐拉梁模型,研究了輸流碳納米管的振動和結(jié)構(gòu)失穩(wěn)分析。稅郎泉等[11]基于經(jīng)典流固耦合方程,采用二階Galerkin方法求解,研究了軸向周期外激勵對含有脈動流體兩端簡支輸流管道橫向振動穩(wěn)定性的影響。Panda等[12]采用多尺度法研究了在內(nèi)共振存在情況時兩端簡支輸流管道的主參數(shù)共振和組合共振。張計光等[13]同樣采用直接多尺度法研究了黏彈性輸流管在Winkler地基上的參數(shù)共振。以上研究均是針對流速無衰減的周期性脈動流作用下的輸流管道橫向振動而進(jìn)行的。

輸流管道在工程應(yīng)用中還可能發(fā)生一種極端現(xiàn)象,即當(dāng)輸水過程中由于斷電停泵或閥門關(guān)閉等原因,將出現(xiàn)水錘現(xiàn)象,管內(nèi)流體流速受流體黏性和管道摩擦阻力等因素影響,將呈現(xiàn)振蕩衰減的特性。對于水錘現(xiàn)象的研究,早期集中于對水錘波的模擬上[14],隨著研究的深入,學(xué)者們主要集中在考慮管道軸向振動與流體間的耦合作用(Fluid-Structure Interaction, FSI)及管壁黏滯作用的影響,在水錘激勵下對管道軸向振動做了大量研究[15-18]。然而,當(dāng)輸流管道中出現(xiàn)水錘現(xiàn)象時,不僅影響管道軸向振動特性,同時也會對管道的橫向振動特性及穩(wěn)定性產(chǎn)生顯著影響。就輸流管道橫向振動而言,水錘產(chǎn)生的振蕩衰減流激勵與流速無衰減的脈動流持續(xù)周期性激勵不同,前者的瞬態(tài)激勵將導(dǎo)致輸流管道的動力不穩(wěn)定性隨著時間的推移而發(fā)生變化,到目前為止還未見相關(guān)的研究報道。

為了探明該內(nèi)激勵型振蕩衰減流對輸流管道穩(wěn)定性的影響,以期在流速無衰減的周期性脈動流對輸流管道振動特性影響的研究基礎(chǔ)上更進(jìn)一步,本文引入指數(shù)衰減函數(shù)模擬水錘發(fā)生后管道中的流速變化過程,基于兩端支撐的輸流管道橫向運動微分方程,推導(dǎo)得到動力不穩(wěn)定區(qū)域的表達(dá)式,研究在流速振蕩衰減型的瞬態(tài)激勵下,兩端支撐輸流管道的不穩(wěn)定區(qū)間,并分析相關(guān)特征參數(shù)對不穩(wěn)定區(qū)域的影響。

1 運動微分方程及離散化

1.1 輸流管道橫向運動微分方程

如圖1所示,水平放置的單跨輸流管道,管道長度為L,管道軸線為x軸,管道橫向為y軸。忽略重力、軸向拉伸和流體壓力的影響,考慮管道黏彈性,其應(yīng)力-應(yīng)變關(guān)系滿足Kelvin-Voiget模型,則兩端支撐的輸流管道橫向運動的微分方程[3]可表述為:

圖1 兩端支撐輸流直管模型

(1)

式中:E*為包含管壁材料Kelvin-Voigt阻尼的彈性模量,Pa·s;I為管道橫截面慣性矩,m4;y為管道橫向位移,m;E為管道彈性模量,Pa;M和m分別為單位長度管內(nèi)流體質(zhì)量和單位長度管道質(zhì)量,kg;U為流體流速,m/s;c為流體黏性引起的摩擦因數(shù);g為重力加速度,m/s2。

引入無量綱變量:

ξ=x/L,η=y/L,χ=[EI(M+m)]-1/2cL2,

代入式(1),整理可得兩端支撐輸流管道無量綱橫向運動微分方程:

(2)

根據(jù)前期對水錘激勵下黏彈性輸流管道振動特性的研究成果[18],本文引入指數(shù)衰減項構(gòu)造隨時間變化的無量綱流速表達(dá)式,用于模擬管道發(fā)生水錘時管內(nèi)流速振蕩衰減特性,即:

u(τ)=u0(ε-μsinωτ)e-bωτ

(3)

式中:u0為管內(nèi)流體的初始流速,μ為流體流速的幅值,ω為流體流速的波動頻率,b為衰減系數(shù),ε為流動系數(shù)。

將式(3)代入式(2),可得:

β1/2u0μωcosωτe-bωτ-β1/2bωu0(ε-

cosωτe-bωτ+β1/2bωu0(ε-μsinωτ)e-bωτ+

(4)

可見,衰減項的引入會使得輸流管道的運動微分方程在對空間和時間的微分項系數(shù)上增加衰減系數(shù),這勢必將對輸流管道的動力不穩(wěn)定區(qū)域帶來一定影響,從而可進(jìn)一步分析流速衰減對輸流管道不穩(wěn)定區(qū)域的影響。

1.2 Galerkin離散

式(4)中包含了空間四階偏導(dǎo)和時間二階偏導(dǎo)項,本文采用二階Galerkin展開式,將輸流管道的橫向位移η(ξ,τ)表述為系統(tǒng)的廣義坐標(biāo)與管道振型函數(shù)的乘積之和,即:

(5)

式中:qr(τ)是離散系統(tǒng)的廣義坐標(biāo);φr(ξ)是同時滿足位移邊界條件和力邊界條件的管道無量綱振型函數(shù)。

將式(5)代入式(4)中,利用歐拉-伯努利梁的振型正交性,并在區(qū)間[0,1]上對ξ積分,可得:

(6)

式中:

M=F+2β1/2u0(ε-μsinωτ)e-bωτB;

K=Λ+[β1/2u0μωcosωτe-bωτ+β1/2bωu0(ε-

e-2bωτ-β1/2u0μωcosωτe-bωτ-β1/2bωu0(ε-

μsinωτ)e-bωτ-γ]C;

針對不同支撐的輸流管道模型,其特征方程特征值λr和振型函數(shù)φr(r=1,2)有所不同,從而使矩陣F、B、C、D和Λ有所不同。采用Bolotin方法對方程(6)應(yīng)用Fourier級數(shù)展開[3,19-20],即可得到內(nèi)激勵型振蕩衰減流作用下輸流管道動力不穩(wěn)定區(qū)域的數(shù)學(xué)表達(dá)式。

2 不穩(wěn)定區(qū)域的確定

2.1 第一動力不穩(wěn)定區(qū)域

對于第一動力不穩(wěn)定區(qū)域,將系統(tǒng)廣義坐標(biāo)q展開成如下形式:

(7)

式中:ak和bk為未知矩陣。

將式(7)代入式(6),經(jīng)化簡合并同類項后,可得:

(8)

式中:

G1=Θ1ak+Θ2bk;G2=-Θ2ak+Θ1bk;

G3=Θ3ak+Θ4bk;G4=Θ5ak+Θ6bk;

G5=Θ7ak+Θ8bk;G6=Θ9ak+Θ10bk;

G7=κ10Cak;G8=κ10Cbk;

Θ2=-κ1F-κ1κ5B;

Θ3=-κ1κ6B+κ7D-κ7C;

Θ4=-κ9D-κ8C+κ9C;

Θ5=κ1κ6B+κ7D-κ6C;

Θ6=κ9D+κ8C-κ9C;

Θ7=κ9D+κ8C-κ9C;

Θ8=-κ1κ6B+κ7D-κ7C;

Θ9=-κ9D-κ8C+κ9C;

Θ10=κ1κ6B+κ7D-κ7C;

κ5=2εβ1/2u0e-bωτ;κ6=β1/2u0μe-bωτ;

由式(8)可看出第一動力不穩(wěn)定區(qū)域由無數(shù)個參數(shù)波疊加而成,根據(jù)Paidoussis等[3]的研究表明當(dāng)k=1時,可得到足夠精度的結(jié)果,因此本文取k=1,并使同類項系數(shù)相等,可得關(guān)于a1和b1的齊次代數(shù)方程組,即

(9)

式中,

(κ3+κ4-κ5+κ7)C

Q12=-κ1F-κ6B-κ9D+(-κ8-κ9)C

Q21=κ1F+κ6B-κ9D+(-κ8+κ9)C

(κ3+κ4-κ5+κ7)C

式(9)周期解存在條件是齊次方程組的行列式等于零,即

|Q|=0

(10)

式(10)可得到振蕩衰減流激勵下不同支撐條件輸流管道的第一動力不穩(wěn)定區(qū)域。

2.2 第二動力不穩(wěn)定區(qū)域

同理,對于第二動力不穩(wěn)定區(qū)域的求解可將系統(tǒng)廣義坐標(biāo)q展開成如下形式:

(11)

將式(11)代入式(6),經(jīng)化簡合并同類項后,取k=2,可得代數(shù)方程組,

(12)

式中:

P11=Λ+κ2D+κ3C+κ4C-κ2C;

P12=-κ9D-κ8C+κ9C;

P13=2κ7B+κ7D-κ7C;

P21=-2κ9D-2κ8C+κ2C;

P23=-κ5B;

P31=κ6D-κ6C;

P32=κ5B+κ10Ca;

同樣,

|P|=0

(13)

時式(12)存在周期解。

由式(13)可得出振蕩衰減流激勵下不同支撐條件輸流管道的第二動力不穩(wěn)定區(qū)域。

3 不穩(wěn)定區(qū)域驗證

前述引入與時間和脈動周期有關(guān)的指數(shù)函數(shù)來描述振蕩衰減流,推導(dǎo)得到了輸流管道在振蕩衰減流激勵下的動力不穩(wěn)定區(qū)域表達(dá)式,若不考慮流速衰減,即流動系數(shù)ε=1和衰減系數(shù)b=0時,則方程退化為流速無衰減脈動流激勵下動力不穩(wěn)定區(qū)域方程。由于目前輸流直管在振蕩衰減流作用下的不穩(wěn)定區(qū)域研究未查閱到相關(guān)文獻(xiàn)成果,本節(jié)針對兩種不同支撐條件的輸流管道模型,選取兩個經(jīng)典案例即Chen[1]和Paidoussis[3]所做研究進(jìn)行對比,用于驗證本文所提出的模型。

3.1 兩端固支管道

對于兩端固支管道,其振型函數(shù)為:

φr=chλrξ-cosλrξ-

(14)

計算得到在ω/ω0-μ平面內(nèi),兩端固支輸流管道在不同初始流速u0與不同黏滯阻尼系數(shù)χ條件下第一動力不穩(wěn)定區(qū)域,如圖2所示(實線為χ=0,虛線為χ=0.2,點劃線為χ=0.5)。可見本文計算結(jié)果與Paidoussis等用脈動流作為激勵條件得到的計算結(jié)果(實心點)吻合良好。從圖中可以看出,隨著初始流速u0的增加,第一動力不穩(wěn)定區(qū)域向下移動,且不穩(wěn)定區(qū)域顯著增大;同時,隨著黏滯阻尼系數(shù)χ的增加,不穩(wěn)定區(qū)域起點向右偏移,且初始流速越小,χ對不穩(wěn)定區(qū)域的影響越大,不穩(wěn)定區(qū)域縮小地更快。

圖2 兩端固支輸流管道第一動力不穩(wěn)定區(qū)域隨平均流速u0和黏滯阻尼χ變化

3.2 兩端簡支管道

對于兩端簡支管道,其振型函數(shù)為:

φr=sinλrξ

(15)

式中:r=1,2,兩端簡支梁特征方程的特征值λ1=π,λ2=2π。本算例中各參數(shù)取值分別為:α=0,γ=10,β1/2=0.2,u0=6,χ=0,即表示在重力作用下,不考慮內(nèi)部管壁材料阻尼作用,忽略流體黏性引起摩擦的相對粗管。

計算得到在ω/ω0-μ平面內(nèi),兩端簡支輸流管道第一和第二動力不穩(wěn)定區(qū)域,如圖3所示。可見本文計算結(jié)果比Chen的結(jié)果(點劃線)不穩(wěn)定區(qū)域開口大,而與Paidoussis等[3]得到的計算結(jié)果(虛線)吻合良好。其原因是Chen[1]在推導(dǎo)輸流管道運動微分方程時未考慮縱向加速度項和管道橫向運動的影響,這充分表明了水流與管道的縱向加速度和管道的橫向運動對動力不穩(wěn)定區(qū)域的影響是不可忽略的。

圖3 兩端簡支輸流管道參數(shù)不穩(wěn)定區(qū)域

4 振蕩衰減流對不穩(wěn)定區(qū)域的影響

當(dāng)輸流管道中發(fā)生水錘現(xiàn)象,其流速會出現(xiàn)雙向衰減流的現(xiàn)象,這種現(xiàn)象反映在本文提出的衰減流表達(dá)形式上,與特征參數(shù)脈動幅值μ、衰減系數(shù)b、流速系數(shù)ε和時間τ有關(guān),本節(jié)以兩端簡支的輸流管道模型為例,分析輸流管道在振蕩衰減流激勵下動力不穩(wěn)定區(qū)域在空間和時間上的變化特性。為分析方便,其它系統(tǒng)參數(shù)分別取α=0,γ=0,β1/2=0.8,表示忽略重力影響下不考慮內(nèi)部管壁材料阻尼作用的相對細(xì)管。

4.1 振蕩衰減流特性

本文所引入振蕩衰減流的流速表達(dá)式,在衰減系數(shù)b、脈動幅值μ和流速系數(shù)ε的取值不同時,其流速的衰減形式有所不同。在給定脈動頻率ω(ω/ω0=1.57)和初始流速u0(u0=0.65)的條件下,可得到管道內(nèi)流體在不同衰減系數(shù)b、脈動幅值μ和流速系數(shù)ε的流速時程,如圖4(a)~(c)所示。從圖中可以看出隨著衰減系數(shù)b的增加,管內(nèi)流速衰減至零的速率增加,管內(nèi)流體恢復(fù)至靜止?fàn)顟B(tài)用時縮短(圖4(a));隨著脈動幅值μ的增加,其流速波峰波谷值隨之增加(圖4(b));同時,流速系數(shù)ε=1.0時,管道內(nèi)的流速雖然有衰減特性,但并未表現(xiàn)出雙向流的特性,仍處于周期性脈動流范疇,只有當(dāng)脈動幅值ε<1.0后,管道內(nèi)的流速逐漸呈現(xiàn)出雙向流的特性(圖4(c))。為了進(jìn)一步驗證所提出的流速表達(dá)式能較好地反映水錘激勵下,輸流管道內(nèi)水流速度呈現(xiàn)出的衰減雙向流特性,應(yīng)用文獻(xiàn)[18]所提出的黏彈性輸流管道模型計算得到本文系統(tǒng)參數(shù)條件下輸流管道內(nèi)的流速時程(ε=0.1,μ=1.1,b=0.02),繪于圖4(d)中,可見本文所提出的流速表達(dá)式與水錘發(fā)生后輸流管道內(nèi)的流速時程吻合良好。

(a)μ=1.1,ε=0.1

4.2 衰減系數(shù)b對不穩(wěn)定區(qū)域的影響

為了進(jìn)一步分析管內(nèi)流速衰減系數(shù)b對輸流直管的不穩(wěn)定區(qū)域的影響,當(dāng)時間τ=0.01時,計算得到兩組初始流速u0在ω/ω0-μ平面內(nèi),兩端簡支輸流直管第一和第二動力不穩(wěn)定區(qū)域,如圖5所示。從圖中可看出,隨著脈動幅值μ的增加,不穩(wěn)定區(qū)域范圍增大,管道不穩(wěn)定性增強;同時,隨著衰減系數(shù)b增加,第一動力和第二動力不穩(wěn)定區(qū)域起點會向下偏移,對比圖5(a)和(b)可見,初始流速u0越大,隨著衰減系數(shù)b的增加,第一動力不穩(wěn)定區(qū)域起點向下偏移越明顯,且不穩(wěn)定區(qū)域越大。以μ=1.1為例,在初始流速u0=0.65時,當(dāng)流速不衰減(b=0)情況下,第一動力不穩(wěn)定區(qū)域在1.74~2.23之間,當(dāng)流速衰減系數(shù)b=0.04時,第一動力不穩(wěn)定區(qū)域縮減至1.70~2.15之間,變化并不明顯;而在初始流速u0=1.30時,當(dāng)流速不衰減(b=0)情況下,第一動力不穩(wěn)定區(qū)域在1.51~2.52之間,當(dāng)流速衰減系數(shù)b增加至0.04時,第一動力不穩(wěn)定區(qū)域縮減至1.26~2.24之間。可見,初始流速對動力不穩(wěn)定區(qū)域的影響比衰減系數(shù)對動力不穩(wěn)定區(qū)域的影響更加顯著。

(a)u0=0.65

4.3 時間τ對不穩(wěn)定區(qū)域的影響

在輸流管道發(fā)生水錘時,流速是不斷衰減的,當(dāng)時間到達(dá)一定時刻,管內(nèi)流體將停止流動,因此,輸流管道的不穩(wěn)定區(qū)域也將隨著時間的推移而發(fā)生變化。為了分析特征參數(shù)時間τ對管道不穩(wěn)定區(qū)域的影響,圖6給出了脈動幅值μ=1.1時,在兩組不同初始流速條件下,不同衰減系數(shù)在ω/ω0-τ平面的動力不穩(wěn)定區(qū)域。可見,隨著時間τ的推移,不穩(wěn)定區(qū)域逐漸閉合,且隨著衰減系數(shù)b的增加,其不穩(wěn)定區(qū)域的閉合時間逐漸縮短,這表明隨著衰減系數(shù)b的增加,流速衰減越快,不穩(wěn)定區(qū)間閉合越快,當(dāng)管道內(nèi)流速衰減至0時,水錘過程結(jié)束,流體靜止,此時管道不穩(wěn)定區(qū)域消失。對比圖6(a)和(b)可得,隨著初始流速u0的增加,其不穩(wěn)定區(qū)域閉合時間會推遲;從圖中同樣可看出,隨著初始流速u0的增加,不穩(wěn)定區(qū)域增大。

(a)u0=0.65

5 結(jié) 論

本文引入指數(shù)衰減函數(shù)模擬水錘發(fā)生后管道中流體速度呈現(xiàn)的振蕩衰減特性,基于兩端支撐的輸流管道橫向運動微分方程,推導(dǎo)得到內(nèi)激勵型振蕩衰減流作用下動力不穩(wěn)定區(qū)域的表達(dá)式。在流速無衰減的脈動流激勵條件下,模擬得到的兩種不同支撐輸流管道不穩(wěn)定區(qū)域與前人數(shù)值結(jié)果吻合良好;為了進(jìn)一步驗證,將本文所提出的流速表達(dá)式結(jié)果與黏彈性輸流管道模型計算得到的流速時程進(jìn)行對比,結(jié)果表明,公式(3)能較好地反映水錘激勵下輸流管道內(nèi)水流速度呈現(xiàn)出的振蕩衰減特性。在此基礎(chǔ)上,分析了不同流速特征參數(shù),如脈動幅值μ、衰減系數(shù)b和初始流速u0,對兩端簡支輸流管道模型的動力不穩(wěn)定區(qū)域的影響。結(jié)果表明,隨著脈動幅值μ的增加,不穩(wěn)定區(qū)域范圍擴(kuò)大,隨著衰減系數(shù)b增加,不穩(wěn)定區(qū)域起點向下偏移;隨著時間τ的推移,不穩(wěn)定區(qū)域逐漸閉合,且隨著衰減系數(shù)b的增加,其不穩(wěn)定區(qū)域的閉合時間逐漸縮短;初始流速對動力不穩(wěn)定區(qū)域的影響比衰減系數(shù)對動力不穩(wěn)定區(qū)域的影響更加顯著。由于輸流直管在內(nèi)激勵型振蕩衰減流的作用下,其動力不穩(wěn)定性隨時間的推移而改變,因此研究其橫向振動響應(yīng)將是下一步的重點。

猜你喜歡
區(qū)域
分割區(qū)域
探尋區(qū)域創(chuàng)新的密碼
科學(xué)(2020年5期)2020-11-26 08:19:22
基于BM3D的復(fù)雜紋理區(qū)域圖像去噪
軟件(2020年3期)2020-04-20 01:45:18
小區(qū)域、大發(fā)展
商周刊(2018年15期)2018-07-27 01:41:20
論“戎”的活動區(qū)域
區(qū)域發(fā)展篇
區(qū)域經(jīng)濟(jì)
關(guān)于四色猜想
分區(qū)域
公司治理與技術(shù)創(chuàng)新:分區(qū)域比較
主站蜘蛛池模板: 欧美中出一区二区| 亚洲精品va| 国产偷倩视频| 毛片免费观看视频| h视频在线播放| 亚洲水蜜桃久久综合网站| 日韩在线视频网| 爱爱影院18禁免费| 国产91麻豆免费观看| 国产高清国内精品福利| a在线亚洲男人的天堂试看| 欧美专区在线观看| 国产亚洲精久久久久久无码AV| 伊在人亚洲香蕉精品播放| 91小视频在线播放| 久久香蕉国产线看精品| 久久精品视频一| 精品伊人久久久大香线蕉欧美| 国产精品久久久久久久伊一| 精品一区二区三区视频免费观看| 五月婷婷导航| 国产白浆在线| 国内精品小视频福利网址| 国产人成网线在线播放va| 亚洲一区二区日韩欧美gif| 国产香蕉97碰碰视频VA碰碰看| 999精品在线视频| 久久 午夜福利 张柏芝| 成人噜噜噜视频在线观看| 激情综合婷婷丁香五月尤物| 少妇精品在线| 成人综合在线观看| 97一区二区在线播放| 日本亚洲欧美在线| 久久精品aⅴ无码中文字幕| 成人av手机在线观看| 国产男女XX00免费观看| 奇米精品一区二区三区在线观看| 色综合色国产热无码一| 久久成人国产精品免费软件| 国产微拍一区二区三区四区| 99无码中文字幕视频| 在线va视频| 久久精品亚洲热综合一区二区| 亚洲av无码片一区二区三区| 日本91在线| 伊在人亚洲香蕉精品播放| 国产国模一区二区三区四区| 国产亚洲欧美在线中文bt天堂 | 亚洲天堂网在线播放| 中文字幕久久精品波多野结| 亚洲综合经典在线一区二区| 影音先锋丝袜制服| 亚洲三级视频在线观看| 永久成人无码激情视频免费| 色天天综合| 青草精品视频| 成年人福利视频| 91探花国产综合在线精品| 成人中文字幕在线| 日韩a级毛片| 麻豆国产在线不卡一区二区| 日本精品αv中文字幕| 国产无码在线调教| 国产波多野结衣中文在线播放| 91精品网站| 欧美国产成人在线| 国产美女91视频| 日本一区二区三区精品视频| 国产精品一区二区不卡的视频| 久久精品无码国产一区二区三区| 日本免费一区视频| 91精品视频网站| 综合人妻久久一区二区精品| 精品人妻AV区| 色婷婷亚洲综合五月| 2020国产精品视频| 福利国产微拍广场一区视频在线| 青青青国产在线播放| 无码综合天天久久综合网| 日本影院一区| 欧美国产日韩另类|