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

水下細(xì)長體結(jié)構(gòu)“流-剛-彈耦合”特性及方法

2014-07-25 11:28:27劉志忠
艦船科學(xué)技術(shù) 2014年9期
關(guān)鍵詞:振動結(jié)構(gòu)

張 成,劉志忠

(中國艦船研究設(shè)計中心,湖北 武漢 400064)

水下細(xì)長體結(jié)構(gòu)“流-剛-彈耦合”特性及方法

張 成,劉志忠

(中國艦船研究設(shè)計中心,湖北 武漢 400064)

為對水下細(xì)長體的流固耦合特性進(jìn)行分析,本文建立細(xì)長體結(jié)構(gòu)的平面運(yùn)動學(xué)方程。為捕捉結(jié)構(gòu)的振動特性,采用雙漸進(jìn)法給出與之相應(yīng)的流體動力學(xué)方程,形成以“流-剛-彈耦合”為特點(diǎn)的方程組。在時域中推進(jìn)編程求解耦合方程組,得到結(jié)構(gòu)的剛體運(yùn)動軌跡、彈性振動情況。并與傳統(tǒng)有限元方法進(jìn)行對比,驗證本文方法的正確性及合理性。以此為理論依據(jù),對潛艇結(jié)構(gòu)的“流-剛-彈耦合”特性進(jìn)行分析。研究表明:水下細(xì)長結(jié)構(gòu)在水動壓力的作用下,流體、剛體運(yùn)動、彈性振動存在著耦合現(xiàn)象。本文將剛體運(yùn)動與彈性振動分別求出,直觀體現(xiàn)了“流-剛-彈耦合”這一現(xiàn)象。

流-剛-彈耦合;雙漸進(jìn)法;水彈性振動;流固耦合

0 引 言

水下結(jié)構(gòu)物在水動壓力載荷或爆炸沖擊載荷作用下,產(chǎn)生剛體運(yùn)動的同時伴隨著彈性變形,其中彈性變形包含整體的變形和局部的振動。對于水下結(jié)構(gòu)物來說,其剛體運(yùn)動及彈性振動的影響不同。剛體運(yùn)動直接影響結(jié)構(gòu)的運(yùn)動姿態(tài),但剛體運(yùn)動本身所消耗的能量并不會使結(jié)構(gòu)產(chǎn)生應(yīng)力或變形,也不會使結(jié)構(gòu)破壞,但彈性振動會影響結(jié)構(gòu)本身的應(yīng)力以及變形,振動過大還會產(chǎn)生塑性變形,甚至產(chǎn)生破口等破壞現(xiàn)象[1-4]。因此需要將水動力載荷及外力載荷共同作用下結(jié)構(gòu)產(chǎn)生的剛體運(yùn)動與彈性振動分開分析,盡量使外力做功轉(zhuǎn)化為剛體運(yùn)動的動能,減少彈性振動動能,必然有益于降低結(jié)構(gòu)的破壞及毀傷。

目前,在計算結(jié)構(gòu)上某一點(diǎn)的運(yùn)動時一般采用有限元法或邊界元法[5-10],這樣得到的軌跡包括整個結(jié)構(gòu)的剛體運(yùn)動和該點(diǎn)在結(jié)構(gòu)局部坐標(biāo)系下的彈性振動之和。需要研究的是,結(jié)構(gòu)上各點(diǎn)的彈性振動隨時間的變化情況,這就需要把這2種運(yùn)動分離開來,并能直觀地表示。為此,最容易想到的是,用合位移減掉質(zhì)心點(diǎn)的運(yùn)動軌跡即是各點(diǎn)的彈性振動。但是質(zhì)心點(diǎn)的運(yùn)動也包含了它本身的彈性振動。因此必須建立一個描述結(jié)構(gòu)剛體運(yùn)動的坐標(biāo)系和一個描述結(jié)構(gòu)各質(zhì)點(diǎn)局部變形的隨動坐標(biāo)系,剛體位移R加上局部位移r即構(gòu)成質(zhì)點(diǎn)的絕對位移,以此建立各質(zhì)點(diǎn)的運(yùn)動方程[11-13],并分別求解R和r即可將剛體運(yùn)動和彈性振動分開。

本文以細(xì)長對稱柱體為研究對象,建立潛體的運(yùn)動學(xué)方程,同時結(jié)合雙漸進(jìn)法(又叫DAA法)建立流體的動力學(xué)方程[14],研究無限水域“流固耦合”效應(yīng)。研究成果可用于潛射導(dǎo)彈、潛艇、水下管道、深海立柱等海洋工程結(jié)構(gòu)的設(shè)計中。

1 基本理論

1.1 結(jié)構(gòu)運(yùn)動方程

圖1 坐標(biāo)系示意圖Fig.1 The coordinate system

(1)

系統(tǒng)動能表達(dá)式如下:

(2)

系統(tǒng)勢能表達(dá)式如下:

(3)

其中M為結(jié)構(gòu)總質(zhì)量。

物面上任意點(diǎn)的坐標(biāo)用(x,γ)表示,其中γ為截面自z軸順時針方向轉(zhuǎn)角(見圖2),a(x)為x截面的半徑,φ(x)為結(jié)構(gòu)物面的斜率見圖3,作用于面元dA=(a(x)/cosφ(x))dxdγ上的壓力p(x,γ),剪應(yīng)力τm(x,γ)和τγ(x,γ)。令δu′和δw′表示這個面沿x方向和z方向的虛位移(見圖3)。

圖2 微元面Fig.2 The infinite surface

圖3 軸對稱體微段dxFig.3 The infinite length ‘dx’ of the symmetric body

將虛位移分解為dA垂直方向和剪切方向并乘以相應(yīng)的應(yīng)力便得到外力在這個單元面上所做的虛功:

(4)

式中:結(jié)構(gòu)質(zhì)心產(chǎn)生水平虛位移δX投影到u和w方向分別為δXsin?和-δXcos?;質(zhì)心產(chǎn)生垂向虛位移δZ投影到u和w方向分別為δZcos?和δZsin?;產(chǎn)生轉(zhuǎn)動虛位移δ?在u和w方向形成的虛位移為-δ?a(x)cosγ和δ?x。

將總動能K、彈性變形能U和外力所做的虛功δW應(yīng)用到非有勢力做功的哈密爾頓原理[15]。

(5)

式中各變分量δX,δZ,δ?,δu,δw都是獨(dú)立的。可以得出:

(6)

(7)

(8)

(9)

(10)

式(6)和式(7)中,結(jié)構(gòu)的剛體運(yùn)動X,Z與流體壓力p有關(guān)。壓力p是結(jié)構(gòu)與流體交界面處的物面壓力,該壓力p不僅隨流場運(yùn)動發(fā)生變化,而且與結(jié)構(gòu)的物面運(yùn)動(含結(jié)構(gòu)剛體運(yùn)動、彈性振動)息息相關(guān)。在下文DAA法的推導(dǎo)中可進(jìn)一步看到物面壓力p與結(jié)構(gòu)的剛體運(yùn)動、彈性振動有關(guān),所以式(6)和式(7)體現(xiàn)了結(jié)構(gòu)的運(yùn)動與流體壓力之間的流固耦合特性。

式(8)中剛體運(yùn)動參數(shù)θ與彈性振動參數(shù)u和w、流體壓力p有關(guān),而壓力p與結(jié)構(gòu)的物面運(yùn)動(含結(jié)構(gòu)剛體運(yùn)動參數(shù)X,Z,θ,彈性振動參數(shù)u和w)息息相關(guān)。水動壓力把剛體運(yùn)動X,Z,θ和彈性振動u,w耦合在一起。這些耦合現(xiàn)象的本質(zhì)仍然是流固耦合,為了對結(jié)構(gòu)的運(yùn)動(剛體運(yùn)動和彈性振動)進(jìn)行更深入解析,本文提出“流-剛-彈耦合”這一說法。

1.2 考慮可壓縮性的水動壓力計算方法

前面通過推導(dǎo)得出潛體結(jié)構(gòu)的運(yùn)動微分方程,下面討論潛體結(jié)構(gòu)表面水動壓力的計算方法。在計算流固耦合問題時,采用考慮可壓縮性的勢函數(shù)微分方程式[14]:

DAA法就是從上述控制方程的基礎(chǔ)上經(jīng)推導(dǎo)得到。DAA法的基本方程如下[2]:

1)一階雙漸近法方程DAA1

2)二階雙漸近法方程DAA2

(11)

式中:ps為流體中的散射壓力;pi為入射波壓力;Mf為流體質(zhì)量矩陣;Ωf為流體頻率矩陣;Af為流體單元的面積矩陣;uI為入射波速度;x為結(jié)構(gòu)位移;G為坐標(biāo)轉(zhuǎn)換矩陣。

(12)

式中:P=ps+pi,即總的流場動壓力值;φ為流場速度勢。

(13)

式(13)中的流場散射壓力ps通過求解DAA方程式(11)得到,入射波壓力pi作為已知條件,流場總速度勢φ即可得解。在得到流場速度勢的基礎(chǔ)上,可以通過伯努利方程得到流場動壓力計為Pd。

(14)

(15)

式中:V為結(jié)構(gòu)運(yùn)動速度。求解式(15)即得到考慮非線性效應(yīng)的流場動壓力P。

2 方程組的求解

在文獻(xiàn)[14]中對DAA法求流場動壓力P進(jìn)行了詳細(xì)介紹,下面對結(jié)構(gòu)運(yùn)動方程作簡要介紹。

因為結(jié)構(gòu)動力學(xué)響應(yīng)中通常低頻成分是主要的,從計算精度考慮,允許采用較大的時間步長。因此,本文采用應(yīng)用最為廣泛的無條件穩(wěn)定的隱式算法Newmark方法[16]。它采用下列假設(shè):

(16)

(17)

3 方法對比驗證

本文列出了水下細(xì)長結(jié)構(gòu)的水彈性振動耦合方程和水動壓力方程。為驗證本文方法的合理性,通過對細(xì)長圓柱殼在水下運(yùn)動的求解與傳統(tǒng)有限元方法求解的結(jié)果進(jìn)行對比。

圖4 柱體模型Fig.4 The cylinder model

圖5 柱體一階模態(tài)(頻率21 Hz)Fig.5 The first mode of the cylinder(the frequency 21 Hz)

圓柱殼半徑為r=0.5 m,長度L=6 m,殼厚h=5 m,彈性模量取E=1.1 E11,泊松比ε取0.3,總質(zhì)量M=11 386 kg,轉(zhuǎn)動慣量J=28 246 kg·m2。坐標(biāo)系取在其質(zhì)心點(diǎn)處,沿水平方向有3.5 m/s初速度,沿垂直方向有10 m/s的初速度,如圖4所示。圖5為圓柱殼網(wǎng)格剖分情況和一階自由振動模態(tài)。

圖6 質(zhì)心軌跡曲線Fig.6 The trace of the barycenter

圖7 偏角時歷曲線Fig.7 The history curve of the offset angle

圖8 垂向速度時歷曲線Fig.8 The history curve of the vertical velocity

圖9 橫速度時歷曲線Fig.9 The history curve of the horizontal velocity

通過圖6~圖9可以看出,2種解法得出的圓柱殼剛體運(yùn)動軌跡幾乎完全一致,證明本文計算方法與傳統(tǒng)結(jié)構(gòu)動力學(xué)計算方法得出的結(jié)果吻合的很好。但是傳統(tǒng)有限元解法很難直觀得出結(jié)構(gòu)上各點(diǎn)的彈性振動軌跡,在傳統(tǒng)解的剛體運(yùn)動軌跡里包含了彈性振動的成分,由于模型的彈性振動十分微小,故在圖6~圖9中的傳統(tǒng)解中未體現(xiàn)出彈性振動的成分。采用本文計算方法求解得出的彈性振動曲線如圖10所示。

圖10 不同時刻柱振動曲線Fig.10 The vibration curves at the different time

圖11 柱體質(zhì)心彈性振動時歷曲線Fig.11 The history curve of the elastic vibration at the center of the cylinder

圖10表示柱體不同時刻的振動曲線,其中每條曲線表示某一時刻柱體軸線的彈性變形情況。可以看出,該工況下柱體彈性振動主要以1階和2階彈性振動為主。圖11所示為圓柱殼質(zhì)心點(diǎn)處彈性振動的軌跡,該點(diǎn)處的振動軌跡接近于正弦曲線。這主要是由于柱體為對稱結(jié)構(gòu),質(zhì)心處于對稱中心的位置,并且該工況下結(jié)構(gòu)的彈性振動表現(xiàn)為1階彈性振動,故質(zhì)心點(diǎn)處的振動時歷曲線接近正弦曲線是合理的。該曲線振動周期為0.064左右,即其振動頻率為15.6左右,比1階固有頻率小,說明水動壓力產(chǎn)生的附加質(zhì)量降低了柱體的振動頻率。

圖12 柱體頂部彈性振動時歷曲線Fig.12 The history curve of the point on the top of the cylinder

圖13 柱體尾部彈性振動時歷曲線Fig.13 The history curve of the point on the end part of the cylinder

剛彈耦合方法的優(yōu)越性體現(xiàn)在采用傳統(tǒng)的結(jié)構(gòu)動力學(xué)計算方法不能將剛體運(yùn)動與彈性振動分開,而采用本文的方法可將彈性振動單獨(dú)提取出來。

4 潛艇模型在均勻來流作用下的運(yùn)動

以潛艇模型在均勻來流作用下的運(yùn)動為例,研究非均勻截面結(jié)構(gòu)的剛體運(yùn)動和彈性振動。為了減小計算量選取潛艇縮比模型進(jìn)行計算,潛艇縮比模型長L=6 m,最大直徑D=0.7 m,厚h=5 mm,彈性模量E=2.1 E11,總質(zhì)量M=1 574.6 kg,轉(zhuǎn)動慣量I=4 000 kg·m2,垂直潛艇方向均勻來流速度10 m/s,潛艇初始時刻靜止。計算模型如圖14所示。

圖14 潛艇模型Fig.14 The submarine model

圖15給出了均勻流作用下的潛艇運(yùn)動。流體力作用點(diǎn)在潛艇結(jié)構(gòu)質(zhì)心之上,產(chǎn)生一順時針的力矩。使?jié)撏г谒邪l(fā)生偏轉(zhuǎn),直到潛艇的軸線方向與來流方向一致。

采用剛體彈性體耦合的方法計算潛艇運(yùn)動可提取潛艇結(jié)構(gòu)的彈性振動如圖16~19所示。圖16表示柱體不同時刻的振動曲線,其中每條曲線表示某一時刻柱體軸線的彈性變形情況。從圖16可以看出,在本文設(shè)置的工況作用下結(jié)構(gòu)的彈性振動以四階以上彈性振動為主,也有部分時刻以2、3階彈性振動為主。圖17~圖19分別表示潛艇的頭部、中部、尾部各點(diǎn)振動隨時間的變化情況。其中頭部和質(zhì)心點(diǎn)處振動的幅度相近,而尾部的振動幅度較大,這與潛艇的結(jié)構(gòu)特征情況是相吻合的。由圖17中頻域曲線可以看出在該工況下,結(jié)構(gòu)上各質(zhì)點(diǎn)的振動頻率成分分別為3.2 Hz,18.1 Hz,43.0 Hz,表明在水動壓力作用下結(jié)構(gòu)的振動以低頻振動為主。

圖15 潛艇剛體運(yùn)動Fig.15 The rigid movement of the submarine

圖16 不同時刻潛艇振動曲線Fig.16 The vibration of the submarine at different times

圖17 潛艇頭部彈性振動頻域與時域曲線Fig.17 The vibration curve of the top submarine in the frequency domain and the time domain

圖18 潛艇質(zhì)心彈性振動時歷曲線Fig.18 The vibration curve of the barycenter of the submarine

圖19 潛艇尾部彈性振動時歷曲線Fig.19 The vibration curve of the afterbody of the submarine

5 結(jié) 語

本文在計算流體力時,采用考慮可壓縮性的雙漸近法,用于捕捉結(jié)構(gòu)的彈性振動特性,同時將結(jié)構(gòu)的剛體運(yùn)動和彈性振動分離開來,將流固耦合系統(tǒng)中結(jié)構(gòu)的彈性振動特性很直觀的展現(xiàn)出來。得到結(jié)論如下:

1)在水動壓力的作用下,不僅結(jié)構(gòu)與流體之間存在耦合,而且結(jié)構(gòu)的剛體運(yùn)動與彈性振動之間也存在耦合現(xiàn)象。

2)殼體結(jié)構(gòu)在水下運(yùn)動時,其振動特性表現(xiàn)為以低頻振動為主。

3)結(jié)構(gòu)在水中運(yùn)動時的總體彈性振動的動能遠(yuǎn)小于剛體運(yùn)動的動能,因此水結(jié)構(gòu)在水中穩(wěn)定運(yùn)行時,不會產(chǎn)生塑性變形或破壞。

本文為了簡化問題,選取細(xì)長對稱殼體為研究對象,對于艦船等結(jié)構(gòu)不太適合。在以后的研究工作中還將考慮結(jié)構(gòu)的非對稱性,以及考慮非細(xì)長結(jié)構(gòu)的“流剛彈耦合”問題,并將對爆炸沖擊載荷作用下的“流剛彈耦合”現(xiàn)象進(jìn)行分析。

[1] 寧建國,王成,馬天寶.爆炸與沖擊動力學(xué)[M].北京:國防工業(yè)出版社,2010.

[2] 黃曉明,朱錫,牟金磊,等.圓柱殼在水下爆炸作用下鞭狀響應(yīng)試驗研究[J].哈爾濱工程大學(xué)學(xué)報,2010,31(10):1278-1285.

HUANG Xiao-ming,ZHU Xi,MU jin-lei,et al.Study on the whipping response of a stiffened cylindrical shell in an underwater explosion[J].Journal of Harbin Engineering University,2010,31(10):1278-1285.

[3] 董海,劉建湖,吳有生,等.水下爆炸氣泡脈動作用下細(xì)長加筋圓柱殼的鞭狀響應(yīng)分析[J].船舶力學(xué),2007,11( 2):250-258.

DONG Hai,LIU Jian-hu,WU You-sheng,et al.Whipping response analysis of slender stiffened cylindrical shell subjected to underwater explosion with bubble pulse[J].Journal of Ship Mechanics,2007,11(2):250-258.

[4] 李玉節(jié),張效慈,吳有生,等.水下爆炸氣泡激起的船體鞭狀運(yùn)動[J].中國造船,2001,42(3):1-7.

LI Yu-jie,ZHANG Xiao-ci,WU You-sheng,et al.Whipping response of ship hull induced by underwater explosion bybble[J].Shipbuilding of CHINA,2001,42(3):1-7.

[5] GU M X,WU Y S,XIA J Z.Time domain analysis of non-linear hydroelastic response of ships[J].In:Proc.Of 4th PRADS,Varna,Bulgaria,989.

[6] 王朝暉,夏錦祝,吳有生.波浪引起的彈性船體響應(yīng)的時域數(shù)值模擬[J].中國造船,1995,131(4):91-96.

WANG Zhao-hui,XIA Jin-zhu,WU You-sheng.Time domain numerical simulation of wave-induced responses of elastic ships[J].Shipbuilding of CHINA,1995,131(4):91-96.

[7] 顧學(xué)康,段文洋.船體在波浪中的非線性水動壓力[J].船舶力學(xué),2001,5(6):27-39.

GU Xue-kang,DUAN Wen-yang.Nonlinear hydrodynamic pressures on ship hulls in waves[J].Journal of Ship Mechanics, 2001,5(6):27-39.

[8] SIMA C,ZHANG X C,WU Y S. Applications of boundary element method in viscous fluid-structure coupling motion problems[J].J. Ship Mech.,2002(3):35-39.

[9] 陳徐均,吳有生,崔維成.海洋浮體二階非線性水彈性力學(xué)分析——二階力對浮體振動時間響應(yīng)的影響. 船舶力學(xué)[J].2003,7(2):11-20.

CHEN Xu-jun,WU You-sheng,CUI Wei-cheng.Second order nonlinear hydroelastic analyses of floating bodies—effect of second-order forces to time response of vibration[J].Journal of Ship Mechanics,2003,7(2):11-20.

[10] DU S X.A general theory for the hydroelastic response of a structure manoeuvring in viscous fluid[J].J.Ship Mech.,1999(3):21-34.

[11] 陳偉民,李敏,鄭仲欽,等.細(xì)長軸對稱體的水彈性振動特性分析[J].中國科學(xué), 2010, 40(9):1165-1173.

CHEN Wei-min,LI Min,ZHENG Zhong-qin,et al. Hydroelastic dynamic characteristics of a slender axis-symmetric body[J].Physics, Mechanics & Astronomy,2010,40(9):1165-1173.

[12] 夏錦祝.細(xì)長浮體的水彈性力學(xué)理論[D].無錫:中國船舶科學(xué)研究中心,1994.

[13] 陳徐均,董克,潘小強(qiáng),等.海洋浮體水彈性力學(xué)研究歷史與現(xiàn)狀[J].解放軍理工大學(xué)學(xué)報(自然科學(xué)版),2003,4(6):41-49.

CHEN Xu-jun,DONG Ke,PAN Xiao-qiang,et al.Existing hydroelasticity theories of ocean floating bodies[J].Journal of PLA University of Science and Technology,2003,4(6):41-49.

[14] 姚熊亮,孫士麗,陳玉,等.非線性雙漸進(jìn)法應(yīng)用于水中結(jié)構(gòu)瞬態(tài)運(yùn)動的研究[J].振動與沖擊, 2010,29(10):9-15.

YAO Xiong-liang,SUN Shi-li,CHEN Yu,et al.Transientm otions of sub merged structures with nonlinear fluid-structure interaction method[J].Journal of ibration and shock,2010,29(10):9-15.

[15] 鷲津久一郎. 彈性和塑性力學(xué)中的變分法[M].北京:科學(xué)出版社,1984.

[16] 王勖成,邵敏.有限單元法基本原理和數(shù)值方法(2版)[M].清華大學(xué)出版社,1997,448-454.

Research on the fluid-rigid-elastic coupling characteristics of the underwater slender structures and its method

ZHANG Cheng,LIU Zhi-zhong

(China Ship Development and Design Center,Wuhan 430064,China)

To analyze the fluid structure coupling property for the underwater slender structures, in the paper it gives the plane kinematic equations;to get the vibration characteristics, the DAA method is adopted giving the corresponding hydrodynamic equations. And then the equations set characterized by ‘rigid-elastic coupling’and ‘fluid-structure coupling’ are formed. According to the specialties of the movements of the slender structures in water, the rational simplified methods and decoupled methods are given. The coupled equations are solved in the time domain to get the rigid movement of the structure, the elastic vibration and the pressure distribution. The results are compared with that of the traditional methods to validate the correctness and rationality of the method in the paper. Based on the method in the paper, the rigid-elastic coupling characteristics and fluid structure coupling characteristics are analyzed. The conclusion indicates that the rigid movements and the elastic vibration of the underwater slender structures are coupled and these two kinds of movement are solved separately to reflect the coupling phenomenon.

fluid-rigid-elastic coupling;DAA;hydroelastic vibration;fluid structure interaction

2013-09-13;

2013-11-13

張成(1987-),男,碩士,工程師,研究方向為結(jié)構(gòu)動力學(xué)。

U664.113

A

1672-7649(2014)09-0025-07

10.3404/j.issn.1672-7649.2014.09.005

猜你喜歡
振動結(jié)構(gòu)
振動的思考
噴水推進(jìn)高速艇尾部振動響應(yīng)分析
《形而上學(xué)》△卷的結(jié)構(gòu)和位置
This “Singing Highway”plays music
論結(jié)構(gòu)
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結(jié)構(gòu)的應(yīng)用
模具制造(2019年3期)2019-06-06 02:10:54
振動攪拌 震動創(chuàng)新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
論《日出》的結(jié)構(gòu)
創(chuàng)新治理結(jié)構(gòu)促進(jìn)中小企業(yè)持續(xù)成長
主站蜘蛛池模板: 国产性精品| 国产精品观看视频免费完整版| 一级高清毛片免费a级高清毛片| 一区二区日韩国产精久久| 91精品视频网站| 国产99在线| 国产免费a级片| 日韩人妻精品一区| 欧美yw精品日本国产精品| 在线无码九区| 亚洲综合第一页| 最新国产在线| 色噜噜狠狠色综合网图区| 国产毛片片精品天天看视频| 萌白酱国产一区二区| 在线国产综合一区二区三区| 8090成人午夜精品| 色综合婷婷| 91精品福利自产拍在线观看| 久久久久人妻一区精品色奶水| 亚洲精品第一页不卡| 91视频青青草| 色九九视频| 国产欧美中文字幕| 香蕉国产精品视频| 亚洲日韩第九十九页| 欧美黄网站免费观看| 97国产在线视频| 亚洲动漫h| 精品国产美女福到在线直播| 2020最新国产精品视频| 成人亚洲国产| 亚洲AV成人一区二区三区AV| 亚洲欧洲综合| 青青青国产视频手机| 1769国产精品免费视频| 国产精品欧美日本韩免费一区二区三区不卡 | 国产欧美日韩精品第二区| 丝袜国产一区| 91麻豆国产在线| 国产美女精品人人做人人爽| a级毛片在线免费| 一区二区欧美日韩高清免费| 欧美视频免费一区二区三区| 国产精品专区第一页在线观看| 日本少妇又色又爽又高潮| 亚洲精品视频网| 无码免费视频| 亚洲成人网在线播放| 久久 午夜福利 张柏芝| 亚洲人人视频| 国产h视频免费观看| 天堂av综合网| 国产欧美在线| 亚洲男女在线| 精品少妇人妻无码久久| 国产小视频在线高清播放| 成人无码一区二区三区视频在线观看| 国产精品天干天干在线观看 | 女高中生自慰污污网站| 久久77777| 午夜a级毛片| 亚洲嫩模喷白浆| 久久久久亚洲av成人网人人软件| 国产资源站| 国产成人精彩在线视频50| 国产内射在线观看| 久久美女精品| 在线高清亚洲精品二区| 国产色网站| 日本人妻一区二区三区不卡影院| 成人免费午夜视频| 久996视频精品免费观看| 五月婷婷导航| 色屁屁一区二区三区视频国产| 日韩国产一区二区三区无码| 欧美亚洲一区二区三区导航| 亚洲va欧美ⅴa国产va影院| 麻豆精品视频在线原创| 国产精品 欧美激情 在线播放| 精品無碼一區在線觀看 | 亚洲AV色香蕉一区二区|