• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于直曲組集算法的復(fù)雜液壓管路固有頻率分析

    2018-04-24 08:07:13張子駿趙通來劉永壽
    振動(dòng)與沖擊 2018年7期
    關(guān)鍵詞:直管固有頻率管路

    韓 濤, 劉 偉, 張子駿, 趙通來, 劉永壽

    (西北工業(yè)大學(xué) 力學(xué)與土木建筑學(xué)院, 西安 710129)

    航空液壓管道的流致振動(dòng)是液壓系統(tǒng)“跑冒滴漏”的重要原因之一。但航空管道系統(tǒng)走向、布局復(fù)雜。目前航空液壓管路的動(dòng)力學(xué)分析方法存在以下兩個(gè)問題:① 主要針對單個(gè)管路(包括直管、曲管等)的流固耦合振動(dòng)分析,難以反映管路局部系統(tǒng)或者整體系統(tǒng)的動(dòng)態(tài)品質(zhì),難以滿足航空工程復(fù)雜管路分析的需求;② 傳統(tǒng)液壓系統(tǒng)動(dòng)態(tài)特性分析將管路簡化為一維模型,偏重液壓系統(tǒng)功能性質(zhì)分析,但是難以反映空間管道結(jié)構(gòu)的響應(yīng)和振動(dòng)問題。隨著未來戰(zhàn)機(jī)減重和提高機(jī)動(dòng)性與生存力的要求,液壓系統(tǒng)向高壓化發(fā)展,這樣對航空管道的強(qiáng)度和穩(wěn)定性要求越來越高。因此,有必要對復(fù)雜管路系統(tǒng)流致振動(dòng)的分析方法展開進(jìn)一步研究。

    輸流單管(包括直管、曲管等)的流固耦合振動(dòng)問題已有大量的研究,Paisoussis等[1]推導(dǎo)出了輸流直管振動(dòng)方程并對其進(jìn)行穩(wěn)定性分析。Housner[2]研究了兩端支撐直管的動(dòng)力學(xué)特性,他采用Euler-Bernoulli梁模型,忽略重力、結(jié)構(gòu)阻尼、外部拉壓力和流體壓力效應(yīng),不考慮流體黏性、可壓縮性,獲得了輸流直管的橫向振動(dòng)控制方程。輸流曲管的研究過程中一直存在軸線可伸長和不可伸長兩種假定,其主要區(qū)別在于是否考慮軸向力作用。Chen[3-4]應(yīng)用Hamilton原理建立了軸線不可伸長曲管的運(yùn)動(dòng)方程。Misra等[5-6]分別基于軸線可伸長和不可伸長假定對曲管平面內(nèi)和平面外振動(dòng)進(jìn)行了分析。在此基礎(chǔ)上,很多學(xué)者對單個(gè)管路的流固耦合振動(dòng)特性進(jìn)行了深入研究,包括線性穩(wěn)定性研究[7-8]、非線性振動(dòng)研究[9-10]、模型的數(shù)值解法研究[11]以及管道的振動(dòng)控制研究[12-13]等,為輸流單管的流致振動(dòng)分析建立了較為完善的理論體系。

    然而,在實(shí)際工程應(yīng)用中,輸流管路系統(tǒng)的幾何形狀往往比較復(fù)雜。對于這些復(fù)雜管路,很難建立其振動(dòng)分析的解析方法,因此一些數(shù)值方法和半解析方法逐漸被研究者采用,如有限元法[14-15]、動(dòng)力剛度矩陣法[16-21]、傳遞矩陣法[22-23]、子結(jié)構(gòu)法[24]等。其中有限元法通用性最強(qiáng),但其計(jì)算量很大。動(dòng)力剛度矩陣法的計(jì)算精度不依賴于單元數(shù)量,對于尺寸很大的單元它依然能夠保證精度,因而可適用于復(fù)雜管路的振動(dòng)分析。Koo等采用Euler-Bernoulli梁理論建立了直管單元的動(dòng)力剛度矩陣,并采用傳遞矩陣法推導(dǎo)了曲管單元的動(dòng)力剛度矩陣,完成了復(fù)雜管路基于直管Euler梁模型的組集,為復(fù)雜管路的動(dòng)力學(xué)分析提供了重要思路。隨后,Koo 等采用同樣方法對液態(tài)金屬爐熱管系統(tǒng)的動(dòng)態(tài)特性進(jìn)行了研究。Shen等[25]研究了由周期復(fù)合材料結(jié)構(gòu)組成的復(fù)雜管路的減振問題,將管路中的曲管部分劃分為一系列直管單元來近似,基于直管Timoshenko梁模型完成組集。Dai等[26]在輸流直管控制方程中加入軸向組合力項(xiàng),結(jié)合傳遞矩陣法建立管路系統(tǒng)動(dòng)力剛度矩陣,并分析了流速、軸向組合力對單管和復(fù)雜管路頻率響應(yīng)函數(shù)的影響。但是,已有的研究中在構(gòu)建管路系統(tǒng)動(dòng)力剛度矩陣時(shí),都是基于單直管模型,即“直管-直管”組集,未曾實(shí)現(xiàn)“直管-曲管”組集[27]。另外,算法中對曲管部分常采用有限數(shù)目的直管單元來近似,再結(jié)合傳遞矩陣法來推導(dǎo),過程比較繁瑣,還可能導(dǎo)致動(dòng)力剛度矩陣維度增加。

    因此,本文對曲管采用曲管單元的精確模型,提出復(fù)雜管路系統(tǒng)模態(tài)分析的直曲組集算法。首先基于Euler-Bernoulli梁理論,在局部坐標(biāo)系下建立直管單元和曲管單元自由振動(dòng)的離散模式及動(dòng)力剛度矩陣;然后在全局坐標(biāo)系中實(shí)現(xiàn)直曲組集,建立復(fù)雜管路的系統(tǒng)動(dòng)力剛度矩陣和特征方程;最后采用所提算法分析管路布局對“Z”形管路固有頻率的影響規(guī)律,建立經(jīng)驗(yàn)公式,并加以試驗(yàn)驗(yàn)證。

    1 輸流直管的動(dòng)力剛度矩陣

    輸流直管模型,如圖1所示。管道長度為L,管道各處橫截面相同,彈性模量為E,剪切模量G,單位長度質(zhì)量為mp。流體相對管壁的流速恒定且為U,流體單位長度質(zhì)量為mf。定義一個(gè)直角坐標(biāo)系xyz,x軸沿著未變形的管道軸線,y軸,z軸均垂直于未變形的管道軸線。

    圖1 輸流直管模型

    基于以下假設(shè):① 忽略管道材料阻尼影響;② 忽略重力影響;③ 流體無黏、不可壓縮,那么輸流直管的流固耦合振動(dòng)方程可表示如下[28-29]

    EAp?2wx?x2-mfU?2wx?x?t-(mp+mf)?2wx?t2=0

    (1)

    EI?4wk?x4+(mfU2+PiAi)?2wk?x2+2mfU?2wk?x?t+

    (mp+mf)?2wk?t2=0, (k=y,z)

    (2)

    GJ?2φx?x2-Ix?2φx?t2=0

    (3)

    式中:wx為直管的軸向位移;wk為直管的橫向位移;φx為直管沿軸線方向的轉(zhuǎn)角;t為時(shí)間;Pi為內(nèi)壓;Ai為流體的橫截面積;EAp,EI,GJ分別為管道的拉伸剛度、彎曲剛度和扭轉(zhuǎn)剛度。

    式(1)~式(3)分別具有以下形式的解[30-31]

    wx(x,t)=wx(x)exp(iωt)

    wk(x,t)=wk(x)exp(iωt)

    φx(x,t)=φx(x)exp(iωt)

    (4)

    式中:ω為圓頻率;將其代入振動(dòng)方程式(1)~式(3),消去時(shí)間項(xiàng)exp(iωt),得到輸流直管在頻域內(nèi)的振動(dòng)方程,其解的形式可以表達(dá)為

    wx(x)=Aexp(kax)

    wk(x)=Bkexp(kbx)

    φx(x)=Cexp(kcx)

    (5)

    式中:ka,kb,kc分別為軸向波數(shù),橫向波數(shù)和扭轉(zhuǎn)波數(shù),由下列頻散方程決定

    (6)

    (mp+mf)ω2=0

    (7)

    (8)

    式(6)和式(8)有兩個(gè)根,式(7)有四個(gè)根,根據(jù)Euler梁理論,建立直管單元的離散模式,見表1。其中φy(x),φz(x)分別為沿著y軸,z軸方向的轉(zhuǎn)角;Nx(x)為軸向力;Ny(x),Nz(x)分別為剪力沿著y軸,z軸方向的分量;Mx(x)為沿著x軸方向的扭矩;My(x),Mz(x)為彎矩沿著y軸,z軸方向的分量。

    定義直管單元兩端的位移狀態(tài)矢量Ws和力狀態(tài)矢量Fs如下

    (9)

    (10)

    Cs=(A1A2By1By2By3By4

    Bz1Bz2Bz3Bz4C1C2)T

    (11)

    式中:上標(biāo)L為管道單元的左端;R為管道單元的右端;下標(biāo)s為直管單元,結(jié)合表1和式(9)~式(11),則有

    Ws=W1sCs

    (12)

    Fs=D2sCs

    (13)

    由式(12)和式(13)得到輸流直管的動(dòng)力學(xué)關(guān)系

    Fs=DsWs

    (14)

    表1 輸流直管自由振動(dòng)的離散模式

    2 輸流曲管的動(dòng)力剛度矩陣

    輸流曲管模型,如圖2所示。管道的曲率半徑為R,管道各處橫截面相同流體單位長度質(zhì)量為mf。定義一個(gè)曲線坐標(biāo)系xyz,x軸與未變形的曲管軸線相切,y軸垂直于未變形的曲管中線且處于曲管平面內(nèi),z軸垂直于未變形曲管所在的平面。

    圖2 輸流曲管模型

    根據(jù)Paidoussis的研究,對于輸流曲管,如果不考慮管道材料阻尼的影響,管內(nèi)流體無黏、不可壓縮,忽略重力的影響,那么其三維振動(dòng)方程可表示如下

    EI?4wy?s4+1R?3wx?s3+??s(AiPi-Nx)?wy?s+wxR+

    1R(AiPi-Nx)+mfU2?2wy?s2+1R?wx?s+1R+

    2mfU?2wy?t?s+1R?wx?t+(mp+mf)?2wy?t2=0

    (15)

    EI?4wz?s4-1R?2φx?s2-GJR?2φx?s2+1R?2wz?s2+

    ??s(AiPi-Nx)?wz?s+mfU2?2wz?s2+

    2mfU?2wz?t?s+(mp+mf)?2wz?t2=0

    (16)

    EIR?3wy?s3+1R?2wx?s2-??s(AiPi-Nx)+1R×

    (AiPi-Nx)?wy?s+wxR+mfU2R?wy?s+wxR-

    mfU?2wx?t?s-1R?wy?t-(mp+mf)?2wx?t2=0

    (17)

    -GJ?2φx?s2+1R?2wz?s2+EIRφxR-?2wz?s2+Ix?2φx?t2=0

    (18)

    式中:wx,wy,wz分別為曲管的切向位移、徑向位移以及垂直于曲管平面的位移;φx為曲管軸線方向的轉(zhuǎn)角;s為自然坐標(biāo);t為時(shí)間;Ai為流體的橫截面積;Pi為流體內(nèi)壓;管道的拉伸剛度、彎曲剛度、扭轉(zhuǎn)剛度分別用EAp,EI,GJ表示,則管道的軸向力Nx表達(dá)為

    Nx=EApε

    (19)

    在軸線不可伸長的情況下,ε=0;在軸線可伸長的情況下,ε=?wx?s-wyR。不難看出:式(15)和式(17)為曲管的面內(nèi)振動(dòng),式(16)和式(18)為曲管的面外振動(dòng)。

    式(15)~式(18)解的形式設(shè)為

    wx(s,t)=wx(s)exp(iωt)

    wy(s,t)=wy(s)exp(iωt)

    wz(s,t)=wz(s)exp(iωt)

    φx(s,t)=φx(s)exp(iωt)

    (20)

    將式(20)代回式(15)~式(18),消去時(shí)間項(xiàng)exp(iωt),可得到輸流曲管在頻域內(nèi)的振動(dòng)方程,略去非線性項(xiàng),在軸線可伸長情況下得到輸流曲管的自由振動(dòng)方程如下

    EI?4wy(s)?s4+1R?3wx(s)?s3+AiPi?2wy(s)?s2+1R?wx(s)?s-

    本文根據(jù)試驗(yàn)結(jié)果,通過對不同條件下干燥時(shí)間的分析并結(jié)合文獻(xiàn)[14]得出辣椒干燥時(shí)間的評價(jià)方法,如表2所示。

    EAPR?wx(s)?s-wy(s)R+mfU2?2wy(s)?s2+1R?wx(s)?s+

    2iωmfU?wy(s)?s+wx(s)R-(mp+mf)ω2wy(s)=0

    (21)

    EI?4wz(s)?s4-1R?2φx(s)?s2-GJR?2φx(s)?s2+1R?2wz(s)?s2+

    AiPi?2wz(s)?s2+mfU2?2wz(s)?s2+2iωmfU?wz(s)?s-

    (mp+mf)ω2wz(s)=0

    (22)

    EIR?3wy(s)?s3+1R?2wx(s)?s2+AiPiR+mfUR×

    ?wy(s)?s+wx(s)R+EAp?wx(s)?s-wy(s)R-

    iωmfU?wx(s)?s-wy(s)R+(mp+mf)×

    (23)

    -GJ?2φx(s)?s2+1R?2wz(s)?s2+EIR×

    φx(s)R-?2wz(s)?s2-Ixω2φx(s)=0

    (24)

    式(21)~式(24)具有如下形式的解:

    wx(s)=Aexp(kas)

    wy(s)=A′exp(kas)

    wz(s)=Bexp(kbs)

    φx(s)=B′exp(kbs)

    (25)

    式中:ka為面內(nèi)振動(dòng)的波數(shù),既包含軸向波數(shù),也包含面內(nèi)彎曲波數(shù),由式(26)決定

    deta11a12

    a21a22=0

    (26)

    式中:

    (AiPi+mfU2)/R2

    kb為面外振動(dòng)的波數(shù),既包含扭轉(zhuǎn)波數(shù),也包含面外彎曲波數(shù),由式(27)決定

    (27)

    式中:

    (mp+mf)ω2

    從式(26)和式(27)解得面內(nèi)振動(dòng)包含六個(gè)波數(shù),面外振動(dòng)也包含六個(gè)波數(shù),那么結(jié)合Euler梁理論,可以建立曲管單元的離散模式,如表2所示。表中系數(shù)αj,βj可結(jié)合式(26)和式(27)得到

    (28)

    (29)

    結(jié)合表2,同樣可以得到

    Wc=D1cCc

    (30)

    Fc=D2cCc

    (31)

    式中:Wc為位移狀態(tài)矢量,與式(9)相同,F(xiàn)c為力狀態(tài)矢量,與式(10)相同,下標(biāo)c為曲管單元,系數(shù)列陣為

    Cc=(A1A2A3A4A5A6

    B1B2B3B4B5B6)T

    (32)

    聯(lián)立式(30)~式(32)得到

    Fc=DcWc

    (33)

    3 基于動(dòng)力剛度矩陣法的直曲組集

    對于復(fù)雜管路系統(tǒng),其直曲組集計(jì)算過程,如圖3所示。按照“先分解再組集”的思路,將復(fù)雜管路系統(tǒng)分解為單管,依據(jù)Euler梁理論,建立單管的局部動(dòng)力剛度矩陣,然后運(yùn)用轉(zhuǎn)換矩陣建立全局坐標(biāo)系下單管的動(dòng)力學(xué)關(guān)系,再按照節(jié)點(diǎn)進(jìn)行組集,建立起系統(tǒng)的動(dòng)力剛度矩陣,結(jié)合邊界條件構(gòu)建系統(tǒng)特征方程,最后進(jìn)行求解。

    以管路模型為例進(jìn)行詳細(xì)說明,如圖4所示。首先將管路分解為單直管和單曲管,單管動(dòng)力剛度矩陣的建立過程如前文所述,第i個(gè)管道單元的動(dòng)力剛度矩陣記為Di,根據(jù)式(14)和式(33),單管在局部坐標(biāo)系下的動(dòng)力學(xué)關(guān)系為

    表2 輸流曲管自由振動(dòng)的離散模式

    Fi=DiWi

    (34)

    圖3 算法流程圖

    建立全局坐標(biāo)系x0y0z0和局部坐標(biāo)系x1y1z1,x2y2z2,…,xiyizi,…,xmymzm,m為管路系統(tǒng)中直管單元的數(shù)目。那么依據(jù)局部坐標(biāo)系xiyizi相對于全局坐標(biāo)系x0y0z0的方向余弦矩陣ti,可以得到管道單元i的轉(zhuǎn)換矩陣Ti,將管道單元i的力狀態(tài)矢量和位移狀態(tài)矢量轉(zhuǎn)換到全局坐標(biāo)系下

    Fig=TiFi

    (35)

    Wig=TiWi

    (36)

    圖4 管路組集示意圖

    將式(35)和式(36)代入式(34),得到

    (37)

    式中:下標(biāo)g為在全局坐標(biāo)系下,由于轉(zhuǎn)換矩陣Ti為正交矩陣,故管道單元i在全局坐標(biāo)下的動(dòng)力剛度矩陣為

    (38)

    那么,對于單元i,其在全局坐標(biāo)系下的動(dòng)力學(xué)關(guān)系為

    Fig=DigWig

    (39)

    將其分解為

    (40)

    同理,單元i+1的動(dòng)力學(xué)關(guān)系可表示為

    (41)

    管道單元i與單元i+1的交點(diǎn)記為節(jié)點(diǎn)i,其受力如圖5所示。在工程應(yīng)用中,直管單元與曲管單元在節(jié)點(diǎn)處相切或近似相切的情形比較常見,但也有不相切的情形,這兩種情形在全局坐標(biāo)系下的受力分析可以得到統(tǒng)一,也就是說,無論節(jié)點(diǎn)i在局部坐標(biāo)系下受力情況如何,都可以通過轉(zhuǎn)換矩陣將其所受力轉(zhuǎn)換為沿著全局坐標(biāo)軸x0,y0,z0方向的分量,節(jié)點(diǎn)i處的力矩分析也是如此。

    圖5 節(jié)點(diǎn)i受力分析圖

    在全局坐標(biāo)系下,根據(jù)該節(jié)點(diǎn)處的力平衡條件和位移連續(xù)性條件,可以得到

    (42)

    (43)

    (44)

    至此,便完成了“直管-曲管”的組集過程,依次類推,依據(jù)節(jié)點(diǎn)建立整個(gè)管路系統(tǒng)的動(dòng)力學(xué)關(guān)系如式(45)所示

    (45)

    FK=-KW(ω)

    (46)

    K為彈性系數(shù),節(jié)點(diǎn)i處受力平衡,得到

    (47)

    依據(jù)上式將式(46)疊加到式(45)中,即可得到含彈性支撐的系統(tǒng)動(dòng)力剛度矩陣。

    4 算 例

    4.1 計(jì)算驗(yàn)證

    直曲組集算法可用于含有“直管-曲管”典型管道元件的復(fù)雜管路計(jì)算,但是目前關(guān)于這類管路的理論成果較少,缺乏對比性,而關(guān)于單管的計(jì)算結(jié)果較多。如果所提算法是正確的,那么必然也可適用于單管的計(jì)算。所以先采用該法對單曲管進(jìn)行計(jì)算,并與文獻(xiàn)[18]進(jìn)行對比,對所提算法進(jìn)行初步驗(yàn)證。

    為了方便對比,管道材料參數(shù)及邊界條件設(shè)置與文獻(xiàn)[18]相同:楊氏模量E=210 GPa,截面慣性矩Ip=9.41×10-8m4,管截面面積Ap=4.05×10-4m2,流體截面面積Af=1.26×10-3m2,管道單位長度質(zhì)量mp=3.18 kg/m,流體單位長度質(zhì)量mf=1.26 kg/m,管道內(nèi)壓P=10 MPa,流體速度V=15 m/s,曲管的約束條件設(shè)為兩端固支。計(jì)算結(jié)果,如表3所示。

    從表3可知:直曲組集方法計(jì)算的一階、三階、五階頻率與文獻(xiàn)[18]給出的前三階頻率相近,相對誤差2%,但卻多出了兩階頻率,是參考文獻(xiàn)中所沒有的。這是因?yàn)槲墨I(xiàn)[18]中計(jì)算的是曲管的平面振動(dòng),而本文采用的是曲管的三維振動(dòng)模型,不僅包含面內(nèi)振動(dòng),也包含面外振動(dòng),為了進(jìn)一步說明,求得曲管的前五階振型,如表4所示??梢愿鼮榍宄乜闯?,一階、三階、五階頻率是曲管面內(nèi)振動(dòng)頻率,二階、四階頻率是曲管面外振動(dòng)頻率。采用所提算法求得的面內(nèi)振動(dòng)前三階頻率與參考文獻(xiàn)中平面振動(dòng)前三階固有頻率相吻合,初步證明所提算法是正確的。

    表3 曲管固有頻率對比結(jié)果

    表4 曲管的前五階振型

    為了進(jìn)一步證明直曲組集算法可用于復(fù)雜管路,分別采用該法與有限元法對“Z”形管路進(jìn)行計(jì)算。管路參數(shù)設(shè)置如下:楊氏模量E=200 GPa,泊松比υ=0.3,密度ρp=7.93×103kg/m3,管道外徑do=9.9 mm,壁厚t=1.2 mm,中間直管單元長520 mm,兩端直管單元的長度均為250 mm,曲管曲率半徑均為40 mm。管內(nèi)流體密度ρf=898 kg/m3,內(nèi)壓P=10 MPa,流速V=0.5 m/s,約束條件設(shè)為兩端固支。

    直曲組集計(jì)算結(jié)果,如圖6所示,圖6中下尖點(diǎn)對應(yīng)的橫坐標(biāo)即為管路的固有頻率,分別用W1,W2,W3,W4,W5來表示。然后采用有限元法計(jì)算了不同單元數(shù)目下“Z”形管路的前三階固有頻率。對比結(jié)果,如圖7所示。

    圖6 “Z”形管路的前五階頻率圖

    圖7 兩種方法計(jì)算結(jié)果對比

    從圖7可知,隨著劃分單元數(shù)目的增多,有限元法計(jì)算結(jié)果逐漸趨于穩(wěn)定,并與直曲組集計(jì)算結(jié)果逐漸吻合,說明直曲組集算法是正確的。而且,直曲組集算法計(jì)算時(shí)僅用了5個(gè)單元,而有限元法在計(jì)算單元數(shù)目達(dá)到30個(gè)時(shí)才趨于穩(wěn)定,說明直曲組集算法在保證計(jì)算精度的同時(shí),減少了計(jì)算單元數(shù)目,實(shí)現(xiàn)了大尺寸單元組集計(jì)算。

    4.2 管路布局對管路固有頻率的影響

    管路布局對管路的振動(dòng)特性有著重要影響,因?yàn)楣苈凡季职l(fā)生變化時(shí),管路的固有頻率也會(huì)隨之發(fā)生變化,而固有頻率是判斷管路穩(wěn)定性的一個(gè)重要參數(shù)。研究管路布局對管路固有頻率的影響規(guī)律,有助于管路設(shè)計(jì)時(shí),提高管路的穩(wěn)定性,避開激振頻率,預(yù)防發(fā)生共振等。所以,接下來以航空管路系統(tǒng)中常見的“Z”形管路為例,探討管路布局對其固有頻率的影響規(guī)律。

    “Z”形管A端、B端位置固定,如圖8所示。定義兩個(gè)幾何參數(shù):曲管至A端的距離l以及曲管的曲率半徑r。那么管路布局的變化可通過改變l和r的值來實(shí)現(xiàn)。為了方便尋找規(guī)律,采用控制變量法,即先給定曲管曲率半徑r,通過改變l來觀察其對“Z”形管路固有頻率的影響;然后再給定l,改變r(jià)值,觀察“Z”形管路固有頻率的變化規(guī)律。

    圖8 “Z”形管路布局圖

    管道材料參數(shù)設(shè)置同上,AB兩端點(diǎn)之間的橫向距離為600 mm,豎向距離為580 mm,管內(nèi)流體密度為ρf=1 000 kg/m3,內(nèi)壓Pi=21 MPa,流速U=12.48 m/s。根據(jù)航空管路設(shè)計(jì)標(biāo)準(zhǔn),曲管的曲率半徑不小于管道外徑的4倍,即r≥4do,故取r=40 mm,分別計(jì)算了l=80 mm,150 mm,220 mm,290 mm四種工況下管路的前三階固有頻率,如圖9所示。

    圖9 不同l值下管路的前三階頻率

    從圖9可知,在曲管曲率半徑一定的情形下,隨著曲管位置距離A端越遠(yuǎn),即l值越大,“Z”形管路的前兩階頻率均隨之減小,第三階頻率先減小后增大。

    一階固有頻率是判斷管路靜態(tài)失穩(wěn)的重要參數(shù),此處重點(diǎn)研究“Z”形管路的一階頻率隨l值的變化規(guī)律。定義一個(gè)比例系數(shù)α,α=l/yAB,yAB代表圖8中A端,B端的豎向距離,逐步改變l值,得到“Z”形管路的一階固有頻率隨α的變化,如圖10所示。

    從圖10可知,“Z”形管一階頻率W1隨著比例系數(shù)α的增大先減小后增大,在α=0.5時(shí)達(dá)到最小值,這與管路的幾何形狀有關(guān),管路A端、B端約束條件相同,當(dāng)曲管距離B端的距離達(dá)到與原來距離A端的距離相等(例如l=430 mm與l=150 mm,圖10中α=0.741與α=0.259)時(shí),管路幾何形狀是對稱的,故一階頻率是相同的,振型也是對稱的。觀察圖形可知:W1與α近似呈余弦函數(shù)關(guān)系,不妨假設(shè)

    W1=Acosωα+θ+B

    (48)

    為了進(jìn)一步確定式中待定系數(shù),定義:“Z”形管路幾何形狀呈中心對稱(即α=0.5)時(shí),管路的一階固有頻率為W0.5;α0=r0/yAB,r0表示給定的曲率半徑;當(dāng)α=α0時(shí),管路的一階固有頻率記為Wα0。那么結(jié)合圖10,可以將式(48)改寫為

    (49)

    圖10 一階固有頻率隨比例系數(shù)α的變化規(guī)律

    至此,只需將一組計(jì)算數(shù)據(jù)代入式(49),即可確定ω值,這樣就得到了“Z”形管路一階固有頻率關(guān)于比例系數(shù)α的表達(dá)式。為了驗(yàn)證式(49)是否正確,接下來對其進(jìn)行試驗(yàn)驗(yàn)證,如圖11所示。

    對比結(jié)果,如圖12所示。式(49)與計(jì)算結(jié)果基本吻合,與試驗(yàn)結(jié)果雖有誤差,但變化趨勢基本相符。而且通過數(shù)值分析發(fā)現(xiàn),試驗(yàn)結(jié)果與式(49)計(jì)算結(jié)果的相對誤差<5%,說明在誤差許可的范圍內(nèi),式(49)可用于“Z”形管路的設(shè)計(jì)計(jì)算。

    圖12 不同比例系數(shù)下的結(jié)果對比

    接下來研究曲管曲率半徑r發(fā)生變化時(shí)對“Z”形管路固有頻率的影響。給定l=290 mm,管道材料參數(shù)同上,分別計(jì)算了r=40 mm,80 mm,120 mm,160 mm四種工況下管路的前三階固有頻率,如圖13所示。

    圖13 不同r值下管路的前三階固有頻率

    由圖13可知,隨著曲管曲率半徑r的增大,“Z”形管路的前三階固有頻率均會(huì)隨之增大,一階頻率增大得較為緩慢,二階、三階增大地較為明顯。同樣,這里重點(diǎn)研究曲管曲率半徑對“Z”形管路一階固有頻率的影響規(guī)律。

    選取不同r值進(jìn)一步計(jì)算,“Z”形管路一階固有頻率的變化趨勢,如圖14所示。隨著曲管曲率半徑的增大,“Z”形管路的一階頻率增大速率越來越快。采用最小二乘法對計(jì)算結(jié)果分別進(jìn)行二次多項(xiàng)式擬合和三次多項(xiàng)式擬合,發(fā)現(xiàn)擬合曲線近乎重合,故可近似認(rèn)為“Z”形管路的一階固有頻率是曲率半徑r的二次函數(shù),即

    W1=a0+a1r+a2r2

    (50)

    對式(50)進(jìn)行試驗(yàn)驗(yàn)證,分別選取了r1=40 mm,r2=80 mm,r3=120 mm,r4=160 mm,r5=200 mm五種尺寸進(jìn)行試驗(yàn),結(jié)果如圖15所示。二次擬合結(jié)果與計(jì)算結(jié)果基本吻合,試驗(yàn)結(jié)果與擬合曲線的變化趨勢基本相同,但相對擬合結(jié)果偏小,這是因?yàn)楣艿郎腺N有傳感器等附件,給管道增添了附加質(zhì)量,所以導(dǎo)致試驗(yàn)結(jié)果偏小。數(shù)值分析顯示,擬合結(jié)果與試驗(yàn)結(jié)果的相對誤差<5%,如果考慮附加質(zhì)量的影響,相對誤差會(huì)更小,說明擬合結(jié)果是正確的,在誤差許可的范圍內(nèi),式(50)可用于“Z”形管路的設(shè)計(jì)計(jì)算。

    圖14 一階固有頻率隨曲率半徑的變化規(guī)律

    圖15 不同曲率半徑下的結(jié)果對比

    5 結(jié) 論

    利用管單元的動(dòng)力剛度矩陣,建立“直管-曲管”組合管道系統(tǒng)固有模態(tài)組集解法,通過驗(yàn)證與計(jì)算,得到以下結(jié)論:

    (1) 所提算法對曲管的計(jì)算結(jié)果與參考文獻(xiàn)吻合,相對誤差<2%。

    (2) 所提算法對“直管-曲管”組合管道系統(tǒng)是適用的,實(shí)現(xiàn)了大尺寸單元組集計(jì)算,較有限元法減少了計(jì)算單元數(shù)目。

    (3) 試驗(yàn)顯示,“Z”形管基頻關(guān)于布局參數(shù)的經(jīng)驗(yàn)公式誤差<5%,證明所擬公式正確,在誤差許可范圍內(nèi),可用于“Z”形管的設(shè)計(jì)計(jì)算。

    [1] PAISOUSSIS M P, LI G X. Pipes conveying fluid: a model dynamical problem[J]. Journal of Fluids & Structures, 1993, 7(2):137-204.

    [2] HOUSNER G W. Bending vibration of a pipe line containing flowing fluid[J]. Journal of Applied Mechanics, 1952, 19(3):205-208.

    [3] CHEN S S. Vibration and stability of a uniformly curved tube conveying fluid[J]. Journal of Acoustic Society of America, 1972, 51(Sup 1):223-232.

    [4] CHEN S S. Flow induced in-plane instabilities of curved pipes[J]. Nuclear Engineering and Design, 1972, 23(1): 29-38.

    [5] MISRA A K, PAIDOUSSIS M P, VAN K S. On the dynamics of curved pipes transporting fluid. Part I: Inextensible theory[J]. Journal of Fluids and Structures, 1988, 2(3): 221-244.

    [6] MISRA A K, PAIDOUSSIS M P, VAN K S. On the dynamics of curved pipes transporting fluid. PartII: Extensible theory[J]. Journal of Fluids and Structures, 1988, 2(3):245-261.

    [7] YANG Xiaodong, YANG Tianzhi, JIN Jiduo. Dynamic stability of a beam-model viscoelastic pipe for conveying pulsative fluid[J]. Acta Mechanica Solida Sinica, 2007, 20(4):350-356.

    [8] QIAN Qin, WANG Lin, NI Qiao. Vibration and stability of vertical upward-fluid-conveying pipe immersed in rigid cylindrical channel[J]. Acta Mechanica Solida Sinica, 2008, 21(5):331-340.

    [9] MENG Dan, GUO Haiyan, XU Sipeng. Non-linear dynamic model of a fluid-conveying pipe undergoing overall motions[J]. Applied Mathematical Modelling, 2011, 35(2):781-796.

    [10] WANG L, NI Q. A reappraisal of the computational modelling of carbon nanotubes conveying viscous fluid[J]. Mechanics Research Communications, 2009, 36(36):833-837.

    [11] LIANG Feng, WEN Bangchun. Forced vibrations with internal resonance of a pipe conveying fluid under external periodic excitation[J]. Acta Mechanica Solida Sinica, 2011, 24(6):477-483.

    [12] YU Dianlong, WEN Jihong, ZHAO Honggang, et al. Vibration reduction by using the idea of phononic crystals in a pipe-conveying fluid[J]. Journal of Sound and Vibration, 2008, 318(1/2):193-205.

    [13] 金基鐸, 楊曉東, 張宇飛. 固定約束松動(dòng)對輸流管道穩(wěn)定性和臨界流速的影響[J]. 振動(dòng)與沖擊, 2009, 28(6): 95-99.

    JIN Jiduo, YANG Xiaodong, ZHANG Yufei. Analysis of critical flow velocities of pipe conveying fluid under relaxation of boundary conditions[J]. Journal of Sound and Vibration, 2009, 28(6): 95-99.

    [14] JUNG D, CHUNG J. In-plane and out-of-plane motions of an extensible semi-circular pipe conveying fluid[J]. Journal of Sound & Vibration, 2008, 311(1/2):408-420.

    [15] OLSON L G, JAMISON D. Application of a general purpose finite element method to elastic pipes conveying fluid[J]. Journal of Fluids and Structures, 1997, 11(2): 207-222.

    [16] TANG Bin. Combined dynamic stiffness matrix and precise time integration method for transient forced vibration response analysis of beams[J]. Journal of Sound and Vibration, 2008, 309(3):868-876.

    [17] VIOLA E, RICCI P, ALIABADI M H. Free vibration analysis of axially loaded cracked Timoshenko beam structures using the dynamic stiffness method[J]. Journal of Sound & Vibration, 2007, 304(1/2):124-153.

    [18] 李寶輝, 高行山, 劉永壽,等. 輸液曲管平面內(nèi)振動(dòng)的波動(dòng)方法研究[J]. 固體力學(xué)學(xué)報(bào), 2012, 33(3): 1-7.

    LI Baohui, GAO Hangshan, LIU Yongshou, et al. In-plane vibration analysis of curved pipe conveying fluid with wave propagation method[J]. Chinese Journal of Solid Mechanics, 2012, 33(3): 1-7.

    [19] KOO G H, PARK Y S. Vibration analysis of a 3-dimensional piping system conveying fluid by wave approach[J]. International Journal of Pressure Vessels & Piping, 1996, 67(3):249-256.

    [20] KOO G H, YOO B. Dynamic characteristics of KALIMER IHTS hot leg piping system conveying hot liquid sodium[J]. International Journal of Pressure Vessels and Piping, 2000, 77(11):679-689.

    [21] KOO G H, PARK Y S. Vibration reduction by using periodic supports in a piping system[J]. Journal of Sound & Vibration,1998,210(1):53-68.

    [22] DONG M L, CHOI M J, OH T Y. Transfer matrix modelling for the 3-dimensional vibration analysis of piping system containing fluid flow[J]. Journal of Mechanical Science & Technology, 1996, 10(2):180-189.

    [23] JUNG D, CHUNG J, MAZZOLENI A, et al. Dynamic stability of a semi-circular pipe conveying harmonically oscillating fluid[J]. Journal of Sound & Vibration, 2008, 315(1):100-117.

    [24] ALDRAIHEM O J. Analysis of the dynamic stability of collar-stiffened pipes conveying fluid[J]. Journal of Sound & Vibration, 2007, 300(3/4/5):453-465.

    [25] SHEN H, WEN J, YU D, et al. The vibrational properties of a periodic composite pipe in 3D space[J]. Journal of Sound and Vibration, 2009, 328(1):57-70.

    [26] DAI H L, WANG L, QIAN Q, et al. Vibration analysis of three-dimensional pipes conveying fluid with consideration of steady combined force by transfer matrix method[J]. Applied Mathematics & Computation, 2012, 219(5):2453-2464.

    [27] IBRAHIM R A. Overview of mechanics of pipes conveying fluids—Part I: Fundamental studies[J]. ASME Journal of Pressure Vessel Technology, 2010, 132(3):1-32.

    [28] PAIDOUSSIS M P. Fluid-structure interactions: slender structures and axial flow[M]. Salt Lake City, UT:Academic Press, 2004.

    [29] 包日東, 梁峰. 兩端一般支承裂紋管道的動(dòng)力學(xué)特性[J]. 振動(dòng)與沖擊, 2016, 35(7):220-224.

    BAO Ridong, LIANG Feng. Dynamic characteristics of a cracked pipe conveying fluid with both elastically supported ends[J]. Journal of Vibration and Shock, 2016, 35(7):220-224.

    [30] 梁峰,包日東.微尺度輸流管道考慮熱效應(yīng)的流固耦合振動(dòng)分析[J].振動(dòng)與沖擊,2015,34(5):141-144.

    LIANG Feng, BAO Ridong. Fluid-structure interaction of microtubes conveying fluid considering thermal effect[J]. Journal of Vibration and Shock, 2015,34(5):141-144.

    [31] 王忠民, 鄒德志, 姜全友. 彈性地基上輸流管道主參數(shù)共振的主動(dòng)振動(dòng)控制[J]. 振動(dòng)與沖擊, 2016, 35(4):182-187.

    WANG Zhongmin, ZOU Dezhi, JIANG Quanyou. Active vibration control for principal parametric resonance of pipes conveying fluid resting on an elastic foundation[J]. Journal of Vibration and Shock, 2016, 35(4):182-187.

    猜你喜歡
    直管固有頻率管路
    基于水質(zhì)變化的供熱采暖管路設(shè)計(jì)
    現(xiàn)場測定大型水輪發(fā)電機(jī)組軸系的固有頻率
    液壓管路系統(tǒng)隨機(jī)振動(dòng)下疲勞分析
    硅鋼軋制過程中乳化液流量控制解耦研究及應(yīng)用
    山西冶金(2019年2期)2019-05-31 11:30:04
    2018年河南省各省轄市及直管縣(市)專利申請量統(tǒng)計(jì)表(1月)
    河南科技(2018年9期)2018-09-10 07:22:44
    2017年河南省各省轄市及直管縣(市)專利申請量統(tǒng)計(jì)表(12月)
    河南科技(2018年3期)2018-09-10 05:18:39
    2018年河南省各省轄市及直管縣(市)專利申請量統(tǒng)計(jì)表(3月)
    河南科技(2018年12期)2018-09-10 05:12:39
    對直管河道采砂管理的認(rèn)識與思考
    中國水利(2015年16期)2015-02-28 15:14:46
    總溫總壓測頭模態(tài)振型變化規(guī)律研究
    美航天服漏水或因管路堵塞
    太空探索(2014年4期)2014-07-19 10:08:58
    成人美女网站在线观看视频| 国产精品久久久久久av不卡| 久久精品影院6| av视频在线观看入口| 最近的中文字幕免费完整| 99热6这里只有精品| www日本黄色视频网| 日韩,欧美,国产一区二区三区 | 久久久成人免费电影| av在线天堂中文字幕| 亚洲最大成人手机在线| 成人综合一区亚洲| 国产毛片a区久久久久| 国产午夜精品久久久久久一区二区三区| 久久久精品94久久精品| 男女视频在线观看网站免费| 日本三级黄在线观看| 国产探花极品一区二区| 婷婷色麻豆天堂久久 | 麻豆成人av视频| 又粗又爽又猛毛片免费看| 久久人人爽人人爽人人片va| 精品欧美国产一区二区三| 午夜福利成人在线免费观看| 国产精品一区二区在线观看99 | 日韩欧美三级三区| 国产精华一区二区三区| 国产淫语在线视频| 国产一级毛片七仙女欲春2| 久久久久久久久大av| 国产精品精品国产色婷婷| 在线免费十八禁| 成人美女网站在线观看视频| 又粗又爽又猛毛片免费看| 麻豆国产97在线/欧美| 日韩亚洲欧美综合| 99热精品在线国产| 午夜福利网站1000一区二区三区| 成人特级av手机在线观看| 国产精品久久久久久精品电影小说 | 日韩成人伦理影院| 免费看a级黄色片| 日本猛色少妇xxxxx猛交久久| 丰满少妇做爰视频| 免费看av在线观看网站| 欧美97在线视频| .国产精品久久| 亚洲第一区二区三区不卡| 亚洲av男天堂| 久久久久网色| 啦啦啦韩国在线观看视频| 亚洲电影在线观看av| 亚洲一级一片aⅴ在线观看| 午夜福利在线观看吧| 直男gayav资源| 久久久久久久久久久免费av| 国产精品国产三级国产av玫瑰| 亚洲综合精品二区| 亚洲人成网站高清观看| 男女啪啪激烈高潮av片| 日韩成人av中文字幕在线观看| 日本黄色片子视频| 大香蕉久久网| 美女黄网站色视频| 最近的中文字幕免费完整| 小蜜桃在线观看免费完整版高清| 精品人妻视频免费看| 欧美一区二区亚洲| 中文字幕免费在线视频6| 不卡视频在线观看欧美| 26uuu在线亚洲综合色| 国产黄色小视频在线观看| 久久精品国产鲁丝片午夜精品| 一本久久精品| 九色成人免费人妻av| 深爱激情五月婷婷| 婷婷六月久久综合丁香| 爱豆传媒免费全集在线观看| 91久久精品国产一区二区成人| 亚洲国产欧美在线一区| 国产亚洲精品久久久com| 精品99又大又爽又粗少妇毛片| 国产精品电影一区二区三区| 日韩,欧美,国产一区二区三区 | 亚洲高清免费不卡视频| 舔av片在线| 国产精品久久视频播放| 人妻系列 视频| 又爽又黄无遮挡网站| 久久久精品大字幕| 黄色一级大片看看| 亚洲精品久久久久久婷婷小说 | av女优亚洲男人天堂| 麻豆av噜噜一区二区三区| 亚洲av成人精品一区久久| 欧美日本视频| 精品熟女少妇av免费看| 日本午夜av视频| 色综合站精品国产| 一级毛片我不卡| 亚洲精品456在线播放app| 精品国内亚洲2022精品成人| 亚洲成色77777| 亚洲真实伦在线观看| 精品午夜福利在线看| 国模一区二区三区四区视频| 麻豆乱淫一区二区| 秋霞在线观看毛片| 欧美另类亚洲清纯唯美| 久久久久久久亚洲中文字幕| 国产探花极品一区二区| 日韩一本色道免费dvd| 国产精品无大码| 变态另类丝袜制服| 一个人免费在线观看电影| 久久久久网色| 亚洲精品国产av成人精品| 中文字幕精品亚洲无线码一区| 日韩欧美在线乱码| 国产熟女欧美一区二区| 久久久久精品久久久久真实原创| 日韩强制内射视频| 男女视频在线观看网站免费| 午夜精品在线福利| 18禁动态无遮挡网站| 中文字幕av在线有码专区| 国产精品一区二区在线观看99 | 亚洲精华国产精华液的使用体验| 日本黄色视频三级网站网址| 久久精品久久久久久噜噜老黄 | 欧美高清成人免费视频www| 亚洲欧美精品自产自拍| 日本wwww免费看| 超碰av人人做人人爽久久| 亚洲精品国产av成人精品| 精品久久久久久成人av| 91狼人影院| 91午夜精品亚洲一区二区三区| 免费不卡的大黄色大毛片视频在线观看 | 久久鲁丝午夜福利片| 亚洲国产精品合色在线| 两性午夜刺激爽爽歪歪视频在线观看| 亚洲av二区三区四区| 成年女人看的毛片在线观看| 深夜a级毛片| 国产一区二区在线av高清观看| 一本一本综合久久| 干丝袜人妻中文字幕| 亚洲丝袜综合中文字幕| 最近中文字幕2019免费版| 中文字幕av成人在线电影| 日本免费a在线| 亚洲av不卡在线观看| 亚洲成人av在线免费| 国产成人精品婷婷| 亚洲欧美精品综合久久99| 小说图片视频综合网站| 久久综合国产亚洲精品| av线在线观看网站| 一区二区三区高清视频在线| 狂野欧美激情性xxxx在线观看| 极品教师在线视频| 一区二区三区高清视频在线| 男女视频在线观看网站免费| 18+在线观看网站| 26uuu在线亚洲综合色| 久久国内精品自在自线图片| 91精品一卡2卡3卡4卡| 国产一区二区在线观看日韩| 国产真实乱freesex| 国产伦理片在线播放av一区| 亚洲国产欧美在线一区| 精品久久久久久成人av| 2022亚洲国产成人精品| 99久久精品国产国产毛片| 国产午夜精品一二区理论片| 欧美日韩精品成人综合77777| 最近手机中文字幕大全| 两性午夜刺激爽爽歪歪视频在线观看| 99久久中文字幕三级久久日本| 在线免费观看不下载黄p国产| 国产探花极品一区二区| 男女边吃奶边做爰视频| 久久精品国产鲁丝片午夜精品| 成人毛片a级毛片在线播放| 能在线免费观看的黄片| 欧美日韩综合久久久久久| 亚洲精品国产成人久久av| 99久国产av精品| 日本黄大片高清| 视频中文字幕在线观看| 国产欧美日韩精品一区二区| 免费播放大片免费观看视频在线观看 | 蜜桃久久精品国产亚洲av| 国产精品一区二区性色av| 日本免费a在线| av在线老鸭窝| 日本五十路高清| 亚洲精品影视一区二区三区av| 可以在线观看毛片的网站| 亚洲成人中文字幕在线播放| 国产麻豆成人av免费视频| 色综合色国产| 日韩强制内射视频| 夜夜爽夜夜爽视频| 亚洲国产精品专区欧美| 亚洲一区高清亚洲精品| 人妻夜夜爽99麻豆av| 国产黄色小视频在线观看| 亚洲精品一区蜜桃| 日本色播在线视频| 简卡轻食公司| 欧美区成人在线视频| 狂野欧美激情性xxxx在线观看| 亚洲国产精品国产精品| 97人妻精品一区二区三区麻豆| 日本欧美国产在线视频| 亚洲成色77777| 99热网站在线观看| 成人漫画全彩无遮挡| 51国产日韩欧美| 欧美成人精品欧美一级黄| 日韩视频在线欧美| 天堂中文最新版在线下载 | 国产美女午夜福利| 午夜视频国产福利| 亚洲精品aⅴ在线观看| 国产成人精品婷婷| 午夜老司机福利剧场| 好男人在线观看高清免费视频| 久久久久性生活片| 国产亚洲精品久久久com| videos熟女内射| 在线观看66精品国产| 色综合亚洲欧美另类图片| 国产国拍精品亚洲av在线观看| 最近视频中文字幕2019在线8| 99九九线精品视频在线观看视频| 久久人人爽人人片av| 成人一区二区视频在线观看| 午夜福利网站1000一区二区三区| 久久精品综合一区二区三区| 综合色丁香网| 成人二区视频| www.av在线官网国产| 国产精品福利在线免费观看| 国产精品女同一区二区软件| 极品教师在线视频| av福利片在线观看| 日韩三级伦理在线观看| 亚洲av男天堂| 91在线精品国自产拍蜜月| 国产精品久久电影中文字幕| 伦理电影大哥的女人| 亚洲精品,欧美精品| av在线播放精品| 成人亚洲精品av一区二区| 国产精品乱码一区二三区的特点| 久久国产乱子免费精品| 美女xxoo啪啪120秒动态图| 久热久热在线精品观看| 成人午夜精彩视频在线观看| 村上凉子中文字幕在线| 国产精品久久久久久av不卡| www.av在线官网国产| 欧美性猛交黑人性爽| 丝袜喷水一区| 亚洲综合色惰| av线在线观看网站| 午夜a级毛片| 在线免费观看不下载黄p国产| 在线a可以看的网站| 日韩成人av中文字幕在线观看| 日韩视频在线欧美| 日韩三级伦理在线观看| 中文在线观看免费www的网站| 久久99蜜桃精品久久| 久久久久网色| 国产免费福利视频在线观看| 美女cb高潮喷水在线观看| 黄色配什么色好看| 午夜亚洲福利在线播放| 亚洲成av人片在线播放无| 亚洲精品乱久久久久久| 级片在线观看| 国产精品久久久久久精品电影| 午夜精品在线福利| 欧美三级亚洲精品| 日韩欧美精品免费久久| 国产精品1区2区在线观看.| 99九九线精品视频在线观看视频| 日韩欧美在线乱码| 九九热线精品视视频播放| 丝袜喷水一区| 国产精品一及| 国产在视频线在精品| 看非洲黑人一级黄片| 麻豆乱淫一区二区| 亚洲av中文av极速乱| 国产片特级美女逼逼视频| 亚洲成人av在线免费| 嫩草影院新地址| 日韩一区二区视频免费看| 免费不卡的大黄色大毛片视频在线观看 | 欧美极品一区二区三区四区| 中文字幕av成人在线电影| 日本一二三区视频观看| 熟女电影av网| 男人的好看免费观看在线视频| 亚洲色图av天堂| 最近视频中文字幕2019在线8| 一级毛片aaaaaa免费看小| 久久久久精品久久久久真实原创| 久久久国产成人免费| 免费无遮挡裸体视频| 亚洲欧美精品自产自拍| 99热这里只有精品一区| 欧美激情国产日韩精品一区| 久久久亚洲精品成人影院| 一夜夜www| 成人午夜精彩视频在线观看| 亚洲怡红院男人天堂| 成人性生交大片免费视频hd| 日韩国内少妇激情av| 黄片无遮挡物在线观看| 精品免费久久久久久久清纯| 久久久a久久爽久久v久久| 日本免费在线观看一区| 国产精品女同一区二区软件| 亚洲天堂国产精品一区在线| 老司机影院毛片| 欧美3d第一页| www.av在线官网国产| 中文精品一卡2卡3卡4更新| 韩国高清视频一区二区三区| 天天一区二区日本电影三级| 99热网站在线观看| 三级毛片av免费| av在线蜜桃| 亚洲精品国产av成人精品| 国产高清不卡午夜福利| 亚洲综合精品二区| 亚洲第一区二区三区不卡| 亚洲va在线va天堂va国产| 免费无遮挡裸体视频| av视频在线观看入口| 日本熟妇午夜| 精品免费久久久久久久清纯| www.av在线官网国产| 久久99热这里只频精品6学生 | 国产精品福利在线免费观看| 男人狂女人下面高潮的视频| 免费观看a级毛片全部| 美女高潮的动态| 自拍偷自拍亚洲精品老妇| 久久人人爽人人片av| 国产三级中文精品| 免费在线观看成人毛片| 日韩av不卡免费在线播放| 熟女人妻精品中文字幕| 又粗又爽又猛毛片免费看| 国产精品美女特级片免费视频播放器| 欧美一区二区精品小视频在线| 99久久精品国产国产毛片| 日本免费在线观看一区| 日本一二三区视频观看| 成人毛片60女人毛片免费| 久久久久网色| 日韩人妻高清精品专区| 亚洲人与动物交配视频| 亚洲精品久久久久久婷婷小说 | 1000部很黄的大片| 日韩 亚洲 欧美在线| 成人毛片a级毛片在线播放| 亚洲国产色片| 中文字幕av成人在线电影| 亚洲人与动物交配视频| 人人妻人人澡欧美一区二区| 99久国产av精品| 亚洲精品成人久久久久久| 老司机福利观看| 97人妻精品一区二区三区麻豆| 蜜桃亚洲精品一区二区三区| 久久久久精品久久久久真实原创| 级片在线观看| 中文精品一卡2卡3卡4更新| 大又大粗又爽又黄少妇毛片口| 少妇的逼好多水| 免费一级毛片在线播放高清视频| 国产亚洲午夜精品一区二区久久 | av黄色大香蕉| 国产一区亚洲一区在线观看| 一级毛片aaaaaa免费看小| 男女边吃奶边做爰视频| 91精品一卡2卡3卡4卡| 18禁在线播放成人免费| 级片在线观看| 天堂√8在线中文| 99久久精品一区二区三区| 一级毛片电影观看 | av又黄又爽大尺度在线免费看 | 国产成人91sexporn| 乱码一卡2卡4卡精品| 亚州av有码| 中文字幕久久专区| 日韩强制内射视频| 天天一区二区日本电影三级| 国产精品久久视频播放| 大话2 男鬼变身卡| 亚洲欧美一区二区三区国产| 黄片无遮挡物在线观看| 亚洲久久久久久中文字幕| www日本黄色视频网| 三级毛片av免费| 亚洲性久久影院| 欧美极品一区二区三区四区| 国产v大片淫在线免费观看| 黄色配什么色好看| 精品国产三级普通话版| 三级国产精品片| 亚洲国产高清在线一区二区三| 丝袜美腿在线中文| 中文资源天堂在线| 国产亚洲91精品色在线| 别揉我奶头 嗯啊视频| 日本午夜av视频| 国产一级毛片在线| 日本五十路高清| 亚洲精品国产av成人精品| 亚洲第一区二区三区不卡| 午夜爱爱视频在线播放| 黄片无遮挡物在线观看| 别揉我奶头 嗯啊视频| 亚洲av免费在线观看| 国产伦精品一区二区三区四那| 最近视频中文字幕2019在线8| 搞女人的毛片| 最后的刺客免费高清国语| 寂寞人妻少妇视频99o| 国产美女午夜福利| 91av网一区二区| 色尼玛亚洲综合影院| 国产成年人精品一区二区| 成人一区二区视频在线观看| 亚洲av一区综合| h日本视频在线播放| 春色校园在线视频观看| 日韩精品青青久久久久久| 国产精品伦人一区二区| 久久久久久久久大av| 国产探花在线观看一区二区| 深爱激情五月婷婷| 免费av不卡在线播放| 国产精品嫩草影院av在线观看| 日本黄色视频三级网站网址| 欧美zozozo另类| 国产av码专区亚洲av| 日韩亚洲欧美综合| av播播在线观看一区| 99热这里只有精品一区| 中文精品一卡2卡3卡4更新| 欧美高清成人免费视频www| 两个人视频免费观看高清| 99视频精品全部免费 在线| 国产真实乱freesex| 中文欧美无线码| 99久久人妻综合| 亚洲怡红院男人天堂| 成人综合一区亚洲| 在线观看美女被高潮喷水网站| 免费无遮挡裸体视频| 亚洲av免费高清在线观看| 神马国产精品三级电影在线观看| 精品久久久久久成人av| 久久精品影院6| 国产午夜精品论理片| 91久久精品国产一区二区成人| 免费观看在线日韩| 秋霞在线观看毛片| 亚洲美女视频黄频| 美女被艹到高潮喷水动态| 成人漫画全彩无遮挡| 国产黄片美女视频| 久99久视频精品免费| 少妇人妻一区二区三区视频| 国产 一区 欧美 日韩| 欧美性猛交黑人性爽| www.av在线官网国产| 大香蕉97超碰在线| 日产精品乱码卡一卡2卡三| 国产精品久久电影中文字幕| 国产av不卡久久| 麻豆成人av视频| 成人无遮挡网站| 日韩一本色道免费dvd| 中文字幕亚洲精品专区| 日韩大片免费观看网站 | 亚洲精品日韩在线中文字幕| 中文字幕亚洲精品专区| 天堂√8在线中文| 日韩欧美在线乱码| 一个人观看的视频www高清免费观看| 人人妻人人澡人人爽人人夜夜 | 韩国高清视频一区二区三区| 亚洲国产精品专区欧美| 日韩欧美在线乱码| 国模一区二区三区四区视频| 亚洲精品日韩在线中文字幕| 直男gayav资源| 国产精品一区二区在线观看99 | 国产亚洲一区二区精品| 亚洲综合色惰| 在线播放无遮挡| 毛片一级片免费看久久久久| 三级经典国产精品| 欧美xxxx性猛交bbbb| 免费播放大片免费观看视频在线观看 | 色综合色国产| 简卡轻食公司| 国内揄拍国产精品人妻在线| 秋霞伦理黄片| 久久韩国三级中文字幕| 日韩一本色道免费dvd| 蜜桃久久精品国产亚洲av| 精品久久国产蜜桃| 亚洲色图av天堂| 成人av在线播放网站| 久久精品国产自在天天线| 中文字幕久久专区| 狠狠狠狠99中文字幕| 麻豆久久精品国产亚洲av| 国产v大片淫在线免费观看| 国产精品一区二区三区四区久久| 特级一级黄色大片| 在现免费观看毛片| 3wmmmm亚洲av在线观看| 一级毛片aaaaaa免费看小| av视频在线观看入口| 国产高潮美女av| 精品99又大又爽又粗少妇毛片| 有码 亚洲区| 免费无遮挡裸体视频| 欧美xxxx性猛交bbbb| 草草在线视频免费看| 人人妻人人看人人澡| 亚洲国产最新在线播放| 国产三级中文精品| 99久国产av精品| 亚洲av.av天堂| 一级毛片aaaaaa免费看小| 一二三四中文在线观看免费高清| 国内精品美女久久久久久| 国产单亲对白刺激| 美女黄网站色视频| 亚洲在久久综合| 久久久久网色| 亚洲欧洲国产日韩| 日韩在线高清观看一区二区三区| 一本久久精品| 91av网一区二区| 可以在线观看毛片的网站| 国产成人a∨麻豆精品| 亚洲激情五月婷婷啪啪| 99久久成人亚洲精品观看| 狂野欧美激情性xxxx在线观看| 人体艺术视频欧美日本| 99热这里只有是精品在线观看| 97超视频在线观看视频| 一个人看视频在线观看www免费| 国产高清有码在线观看视频| 亚洲国产高清在线一区二区三| 久久精品国产99精品国产亚洲性色| 国产亚洲91精品色在线| 成年版毛片免费区| 久久久久性生活片| 国产欧美另类精品又又久久亚洲欧美| 人妻制服诱惑在线中文字幕| 国产麻豆成人av免费视频| 国产成人91sexporn| 国产伦一二天堂av在线观看| 亚洲最大成人手机在线| 国产不卡一卡二| 国产又色又爽无遮挡免| 91精品伊人久久大香线蕉| 91午夜精品亚洲一区二区三区| 中文字幕人妻熟人妻熟丝袜美| 免费看美女性在线毛片视频| 亚洲精品aⅴ在线观看| 精品久久久久久成人av| 激情 狠狠 欧美| 永久网站在线| 婷婷色av中文字幕| 亚洲欧洲日产国产| 插逼视频在线观看| 国产单亲对白刺激| 国产国拍精品亚洲av在线观看| 人人妻人人澡人人爽人人夜夜 | 女的被弄到高潮叫床怎么办| 搡女人真爽免费视频火全软件| 尾随美女入室| 国产精品熟女久久久久浪| 黄片wwwwww| 97超碰精品成人国产| 久久国内精品自在自线图片| 最新中文字幕久久久久| 日韩成人伦理影院| 久久久久网色| 亚洲国产精品专区欧美| 精品少妇黑人巨大在线播放 | 国产精品无大码| 秋霞伦理黄片| 亚洲精品国产av成人精品| 久久精品国产自在天天线| 国产精品不卡视频一区二区| 亚洲最大成人手机在线|