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

    低雷諾數(shù)下串列布置雙圓柱渦激振動(dòng)特性研究1)

    2022-03-19 01:54:36涂佳黃黃林茜何永康呂海宇梁經(jīng)群
    力學(xué)學(xué)報(bào) 2022年1期
    關(guān)鍵詞:振動(dòng)

    涂佳黃 黃林茜 何永康 呂海宇 梁經(jīng)群

    (湘潭大學(xué)土木工程與力學(xué)學(xué)院,湖南湘潭 411105)

    引言

    渦激振動(dòng)現(xiàn)象廣泛存在于實(shí)際工程中,如海洋立管[1]、橋梁拉索[2]、高層建筑[3]等.彈性支撐結(jié)構(gòu)體置于一定速度的流場(chǎng)中,會(huì)在其兩側(cè)產(chǎn)生交替脫落的旋渦,并誘發(fā)結(jié)構(gòu)表面產(chǎn)生流體脈動(dòng)力,引起結(jié)構(gòu)振動(dòng).當(dāng)彈性支撐結(jié)構(gòu)的振動(dòng)頻率與渦脫落頻率基本重合時(shí),引發(fā)共振,從而對(duì)結(jié)構(gòu)造成不可逆的破壞[4].因此,對(duì)于渦激振動(dòng)問(wèn)題一直是研究的熱點(diǎn)之一,并取得了一系列成果[5-8].

    在實(shí)際工程中,結(jié)構(gòu)往往以結(jié)構(gòu)群的形式出現(xiàn),因此對(duì)圓柱群渦激振動(dòng)問(wèn)題進(jìn)行研究具有重要意義.影響渦激振動(dòng)的主要參數(shù)包括雷諾數(shù)[9]、剪切率[10]、質(zhì)量比[11]、間距比[12]、頻率比[13]等.近些年來(lái)學(xué)者們關(guān)于各個(gè)參數(shù)對(duì)圓柱渦激振動(dòng)的影響開(kāi)展了大量的研究.Papaioannou 等[14]對(duì)不同間距比工況下串列雙圓柱與單圓柱工況渦激振動(dòng)問(wèn)題進(jìn)行了研究,當(dāng)間距比較小時(shí),上游圓柱所對(duì)應(yīng)的鎖定區(qū)間較單圓柱工況明顯擴(kuò)大,但當(dāng)間距比超過(guò)臨界間距比時(shí),上游圓柱基本不受下游圓柱影響,其振幅與單圓柱工況相似.Mittal 和Kumar[15]對(duì)Re=100 工況下,串列和錯(cuò)列布置雙圓柱體流致振動(dòng)進(jìn)行了數(shù)值模擬,當(dāng)圓柱排布方式從串列布置轉(zhuǎn)變?yōu)椴⒘胁贾脮r(shí),雙圓柱結(jié)構(gòu)運(yùn)動(dòng)軌跡將從“8”字形轉(zhuǎn)變?yōu)椤皺E圓”形.及春寧等[16]對(duì)Re=100 條件下,串列雙圓柱在不同情況下的流致振動(dòng)進(jìn)行了數(shù)值模擬研究,發(fā)現(xiàn)無(wú)論上游圓柱固定或者振動(dòng),下游圓柱的振幅會(huì)強(qiáng)于上游圓柱.Qin 等[17]對(duì)串列布置雙圓柱體進(jìn)行了實(shí)驗(yàn)研究,在不同的頻率比工況下分別發(fā)現(xiàn)了鎖定、渦激發(fā)以及馳振等不同的振動(dòng)機(jī)制.陳文曲等[18]運(yùn)用分塊耦合法對(duì)下游圓柱渦致振動(dòng)的升阻力特性進(jìn)行了模擬研究,發(fā)現(xiàn)結(jié)構(gòu)表面升力的變化和渦脫落模態(tài)以及初始渦脫落時(shí)刻之間的內(nèi)在聯(lián)系.鄒琳等[19]運(yùn)用動(dòng)網(wǎng)格技術(shù)對(duì)Re=100 條件下串列雙圓柱渦激振動(dòng)響應(yīng)進(jìn)行數(shù)值模擬研究,發(fā)現(xiàn)在不同間距比工況下,下游圓柱升阻力系數(shù)在折減速度Ur=7.2 附近會(huì)出現(xiàn)“跳躍現(xiàn)象”,并出現(xiàn)阻力峰值.陳威霖等[20]對(duì)小間距比下串列雙圓柱尾流耦合機(jī)制進(jìn)行了模擬研究,發(fā)現(xiàn)由于串列雙圓柱之間平衡位置差的調(diào)和當(dāng)間距比L/D=1.1 時(shí),在較廣的折減速度區(qū)間內(nèi)圓柱會(huì)發(fā)生大振幅運(yùn)動(dòng).段松長(zhǎng)等[21]研究了錯(cuò)列角度對(duì)雙圓柱渦激振動(dòng)的影響,當(dāng)折減速度較小時(shí),柱體升阻力系數(shù)隨錯(cuò)列角度的增大而增大,當(dāng)錯(cuò)列角度固定時(shí),折減速度的改變對(duì)升力影響較大而對(duì)阻力影響較小.Ryan 等[22]對(duì)層流范圍內(nèi)臨界質(zhì)量比進(jìn)行了探究,當(dāng)Re=40~ 95 范圍內(nèi)存在臨界質(zhì)量比,且其值隨著雷諾數(shù)的增大而減小.Chen 等[23]對(duì)間距比L/D=1.2~ 5.0 范圍內(nèi)串列三圓柱進(jìn)行數(shù)值模擬研究,并分析了平衡位置偏移、低頻振動(dòng)以及漩渦脫落與圓柱之間的時(shí)機(jī)3 個(gè)關(guān)鍵因素對(duì)馳振現(xiàn)象的影響.

    綜上所述,目前研究者們對(duì)于影響多圓柱渦激振動(dòng)參數(shù)的探討大多集中在取兩到3 個(gè)變量,且多集中在對(duì)雷諾數(shù)、間距比、剪切率等關(guān)鍵參數(shù)影響的研究,而對(duì)不同頻率比的影響研究較少.然而,在實(shí)際工程中,結(jié)構(gòu)會(huì)受到多種環(huán)境參數(shù)變量的影響.因此,本文基于四步半隱式特征線(xiàn)分裂算子有限元方法對(duì)層流條件下串列雙圓柱雙自由度渦激振動(dòng)問(wèn)題進(jìn)行了模擬研究,主要分析了頻率比(r=1.0,1.5,2.0),間距比(L/D=2.5,3.5,5.5),剪切率(k=0.0,0.05,0.1)以及折減速度(Ur=3~ 20) 4 個(gè)關(guān)鍵參數(shù)對(duì)結(jié)構(gòu)動(dòng)力響應(yīng)、尾渦特性、頻譜特性、相位特征與能量傳遞方面的影響,并揭示其內(nèi)在機(jī)理,為實(shí)際工程應(yīng)用提供參考依據(jù).

    1 流固耦合數(shù)值方法

    1.1 流體控制方程

    基于任意-拉格朗日歐拉方法下不可壓縮黏性流體N-S 方程的無(wú)量綱形式表達(dá)如下[24]

    式中,xi與xj分別為i,j方向空間坐標(biāo);ui與p為流體速度與壓力;cj=uj-wj,其中wj為網(wǎng)格移動(dòng)速度,cj為流體相對(duì)于wj的對(duì)流速度;t為時(shí)間;雷諾數(shù)Re=UcD/ν,其中D與Uc分別為特征長(zhǎng)度尺寸與中心來(lái)流速度,ν為流體動(dòng)力黏性系數(shù).

    1.2 固體控制方程

    彈性支撐圓柱結(jié)構(gòu)運(yùn)動(dòng)體系可假設(shè)為質(zhì)量-阻尼-彈簧系統(tǒng),雙自由度運(yùn)動(dòng)結(jié)構(gòu)的控制方程無(wú)量綱形式如下

    式中,X與Y分別為結(jié)構(gòu)順流向和橫流向無(wú)量綱位移;ξ為結(jié)構(gòu)阻尼系數(shù),為了得到最大結(jié)構(gòu)位移響應(yīng),取ξ=0;Ur,x=Uc/(fn,xD)和Ur,y=Uc/(fn, yD)分別為x方向(順流向)與y方向(橫流向)折減速度,其中fn,x和fn,y分別為彈性支撐圓柱體結(jié)構(gòu)的順流向和橫流向自然頻率;分別為阻力系數(shù)和升力系數(shù);mr=m/(ρD2)為結(jié)構(gòu)折合質(zhì)量,其中m為單位長(zhǎng)度結(jié)構(gòu)質(zhì)量,ρ為流體密度.本文采用Newmark-β時(shí)間積分法求解結(jié)構(gòu)動(dòng)力方程.

    1.3 網(wǎng)格更新

    為了避免網(wǎng)格失效,本文對(duì)Laplace 方程的邊值問(wèn)題方法進(jìn)行了改進(jìn),采用改進(jìn)的Laplace 方程對(duì)網(wǎng)格坐標(biāo)進(jìn)行更新[25]

    式中,Si為網(wǎng)格節(jié)點(diǎn)i方向位移;Гm和Гf分別為網(wǎng)格動(dòng)邊界和固定邊界;gi是運(yùn)動(dòng)邊界上的節(jié)點(diǎn)位移;τ是網(wǎng)格形變控制參數(shù),表達(dá)式如下

    式中,Δe為計(jì)算網(wǎng)格單元的面積(或體積);Δmin和Δmax分別為網(wǎng)格單元中最小與最大的面積(或體積).本文采用Galerkin 有限元方法求解該Laplace 方程的邊值問(wèn)題.

    1.4 計(jì)算流程

    本文采用分區(qū)迭代方法求解鈍體結(jié)構(gòu)渦激振動(dòng)問(wèn)題,該算法的計(jì)算流程如下:

    (1) 求解流體控制方程:運(yùn)用基于四步半隱式CBS 穩(wěn)定化有限元方法求解流體控制方程式(1)和(2),獲得t(n +1)時(shí)刻流場(chǎng)速度與壓力,從而得到流體作用于結(jié)構(gòu)上的流體力.

    (2)求解固體控制方程:將流體力施加到結(jié)構(gòu),以Newmark-β時(shí)間積分方法求解結(jié)構(gòu)運(yùn)動(dòng)控制方程式(3),得到t(n +1)時(shí)刻結(jié)構(gòu)的動(dòng)力響應(yīng).

    (3) 網(wǎng)格更新:采用Galerkin 有限元方法求解Laplace 方程式(4),獲得網(wǎng)格節(jié)點(diǎn)位移及網(wǎng)格速度,并更新網(wǎng)格坐標(biāo).

    (4)返回第一步開(kāi)始計(jì)算t(n +2)時(shí)刻,并循環(huán)至系統(tǒng)達(dá)到穩(wěn)定為止.

    本文的數(shù)值方法已運(yùn)用于渦激振動(dòng)相關(guān)問(wèn)題求解過(guò)程中,并能得到較好的數(shù)值結(jié)果[26-28],驗(yàn)證其正確性與可靠性.

    2 問(wèn)題描述

    2.1 計(jì)算模型

    基于上述方法本文對(duì)串列雙圓柱流致振動(dòng)問(wèn)題進(jìn)行了數(shù)值模擬研究,其計(jì)算模型如圖1 所示,兩圓柱在橫流向與順流向均可自由振動(dòng),選取雷諾數(shù)Re=100,其他參數(shù)選取為:間距比L/D=2.5,3.5,5.5;剪切率k=0.0,0.05,0.1;固有頻率比r=fn,x/fn,y=1.0,1.5,2.0;折減速度Ur=Ur,x=3~ 20;質(zhì)量比mr=7.85.計(jì)算域尺寸為[-30D,60D] × [-10D,10D],上游圓柱(upstream cylinder,UC)的圓心位置為(0,0),下游圓柱(downstream cylinder,DC)的圓心位置為(L/D,0).本文計(jì)算域堵塞率B=5%,邊界條件設(shè)置如下:模型入口設(shè)定為速度入口ux=Uc+ky(Uc=1.0),uy=0;模型出口設(shè)定為壓力出口p=0;上下邊界均為自由滑移邊界?ux/?y=0,uy=0;結(jié)構(gòu)表面為無(wú)滑移壁面u x=uy=0.對(duì)于彈性支撐圓柱體結(jié)構(gòu),為了獲得較大的結(jié)構(gòu)動(dòng)力響應(yīng),故不考慮阻尼影響,將其簡(jiǎn)化成質(zhì)量-彈簧系統(tǒng)模型,其量綱歸一化時(shí)間步長(zhǎng)設(shè)為Δt=0.002.

    圖1 計(jì)算模型與邊界條件Fig.1 The computational model and boundary conditions

    2.2 網(wǎng)格劃分

    本文流場(chǎng)計(jì)算域采用非結(jié)構(gòu)化三角形網(wǎng)格劃分法,為了更準(zhǔn)確的刻畫(huà)出圓柱尾渦的形成與發(fā)展,故對(duì)圓柱附近及尾流區(qū)域進(jìn)行網(wǎng)格加密處理.由表1 可知,與密網(wǎng)格GIII 相比,網(wǎng)格GI 的Xmax/D與Ymax/D的計(jì)算結(jié)果最大變化率分別為6.25%和4.96%;網(wǎng)格GII 的最大變化率分別下降為0.33%和0.06%.由于網(wǎng)格GII 已滿(mǎn)足數(shù)值模擬結(jié)果網(wǎng)格收斂性要求和計(jì)算精度要求,且若采用密網(wǎng)格GIII,雖能獲得更高的計(jì)算精度,但計(jì)算耗時(shí)會(huì)大幅度增加,效率降低.綜上所述,本文所有算例采用的網(wǎng)格密度和分布情況與GII 類(lèi)似.

    表1 網(wǎng)格獨(dú)立性驗(yàn)證:串列布置雙圓柱流致振動(dòng)計(jì)算結(jié)果(L/D=2.5,r=1.0,k=0.0, Ur=6.0)Table 1 Grid independence test:the results for the two tandem circular cylinders at L/D=2.5,r=1.0,k=0.0, Ur=6.0

    3 計(jì)算結(jié)果與分析

    3.1 振幅特性

    本文對(duì)不同剪切率,固有頻率比,間距比以及折減速度工況下串列布置雙圓柱雙自由度流致振動(dòng)振幅變化情況進(jìn)行了分析.由圖2 可知,圓柱橫流向振幅要遠(yuǎn)大于其順流向振幅,且UC 的共振區(qū)間明顯要寬于DC 的共振區(qū)間.整體上看,頻率比較小工況(r=1.0 與r=1.5)對(duì)UC 振幅值變化影響較小,其幅值隨折減速度的變化較為接近,但當(dāng)r增大到2.0 時(shí),UC 在兩個(gè)自由度方向的振幅明顯增大.一方面,當(dāng)頻率比較大時(shí),x,y方向的結(jié)構(gòu)剛度不同,使其容易受到流體脈動(dòng)力的影響,進(jìn)而導(dǎo)致脈動(dòng)力引發(fā)的振動(dòng)越大.另一方面,可能由于此時(shí)圓柱流向振動(dòng)和橫向振動(dòng)同時(shí)達(dá)到共振有關(guān),雙頻共振的出現(xiàn)使得頻率比為2.0 時(shí),圓柱振幅明顯增大.特別的是,UC 在兩個(gè)自由度方向達(dá)到最大振幅值的折減速度不同,而DC 基本同步.另外,由于來(lái)流以及間隙流的干擾,使得UC 與DC 兩者“鎖定區(qū)間”的折減速度范圍存在差值,UC 較DC 更早進(jìn)入共振狀態(tài),同時(shí)也會(huì)更早退出.

    圖2 不同頻率比與剪切率工況下,串列雙圓柱雙自由度最大振幅隨折減速度的變化(L/D=2.5)Fig.2 The variation of the maximum vibrating amplitude of tandem double cylinders with reduced velocity under different frequency ratios and shear rates in two degrees of freedom (L/D=2.5)

    圖2 不同頻率比與剪切率工況下,串列雙圓柱雙自由度最大振幅隨折減速度的變化(L/D=2.5) (續(xù))Fig.2 The variation of the maximum vibrating amplitude of tandem double cylinders with reduced velocity under different frequency ratios and shear rates in two degrees of freedom (L/D=2.5) (continued)

    3.1.1L/D=2.5

    當(dāng)折減速度(Ur=3~ 4) 較小時(shí),兩圓柱在x,y方向僅微弱振動(dòng),振幅接近于零.當(dāng)Ur=5 時(shí),振幅開(kāi)始“啟動(dòng)”,UC 在橫流向的振幅值迅速增大,逐漸進(jìn)入鎖定狀態(tài).然而,順流方向僅在r=2.0 工況下圓柱體結(jié)構(gòu)會(huì)發(fā)生明顯的振幅,在其它頻率比(r=1.0,1.5)工況下,UC 基本不振動(dòng),如圖2 所示.由圖3 可知,r=1.0 與r=1.5 工況所受到的流場(chǎng)時(shí)均壓力基本相同,UC 迎流面受到來(lái)流的沖擊作用受到正壓,其尾部由于渦流的脫落產(chǎn)生吸力受到負(fù)壓作用.當(dāng)r=2.0 時(shí),UC 受到的正壓面積更廣,與此同時(shí)DC 的前駐點(diǎn)附近開(kāi)始出現(xiàn)正壓,兩圓柱尾部負(fù)壓面積均增大,這可能是該工況比其它兩個(gè)工況振幅更大的原因,如圖3(c)所示.隨Ur進(jìn)一步的增大,圓柱開(kāi)始進(jìn)入“鎖定區(qū)間”(Ur=5~ 10),在此區(qū)間內(nèi)圓柱的振動(dòng)頻率與渦脫落頻率無(wú)限接近,從而誘發(fā)共振現(xiàn)象.由圖2 觀(guān)察到兩圓柱在橫流向的振幅值迅速增大,并在不同的折減速度處取得極值.另外,UC 在不同自由度方向上被“完全鎖定”的折減速度不同,其橫流向振幅最大值在Ur=7 處取得,其值為0.82(k=0.05,r=2.0).然而,在Ur=5 時(shí)順流向振幅達(dá)到最大值,其值為0.19(k=0.0,r=2.0).對(duì)于下游圓柱則沒(méi)有出現(xiàn)錯(cuò)位現(xiàn)象,DC 在兩個(gè)自由度方向上被“完全鎖定”的折減速度相同,其橫流向上最大振幅值達(dá)到了1.36(k=0.0,r=2.0)是UC 的1.7 倍,順流向上最大振幅值為0.24(k=0.0,r=2.0)是UC 的1.3 倍.當(dāng)Ur≥9 時(shí),UC 退出鎖定區(qū)間,在兩個(gè)流向上的振幅均趨于穩(wěn)定,橫流向最大振幅值維持在0.2 附近,而順流向最大振幅值已趨于零.然而,DC 橫流向的振幅在Ur=9 時(shí),還維持著較大的振幅值1.0.隨著Ur的繼續(xù)增大DC 在x方向的振幅值為0,在y方向的振幅由1.0 遞減到0.2 左右,當(dāng)Ur≥16 時(shí),UC 與DC 在兩個(gè)流向的振幅值均趨于0,兩圓柱體結(jié)構(gòu)不再發(fā)生振動(dòng).

    圖3 k=0.0 與Ur=7 時(shí),不同頻率比下流場(chǎng)無(wú)量綱時(shí)均壓力圖Fig.3 Dimensionless time-averaged pressure diagram of flow field at different frequency ratios at k=0.0 and Ur=7

    3.1.2L/D=3.5

    當(dāng)間距比增大到3.5 時(shí),UC 與DC 在兩個(gè)自由度上最大振幅值變化規(guī)律與L/D=2.5 工況大致相同,UC 在Ur=6 處Ymax/D值最大,其最大值為0.79(k=0.05,r=2.0),而DC 在Ur=8 處Ymax/D值最大,其最大值為1.33(k=0.1,r=2.0),如圖4 所示.其中,UC 的振幅變化趨勢(shì)較為規(guī)律,而DC 較為特殊.特殊在于:均勻來(lái)流工況下,DC 在x與y方向上最大振幅值出現(xiàn)了“平臺(tái)期”,其“鎖定區(qū)域”圖形類(lèi)似于梯形,在Ur=7 與8 時(shí)結(jié)構(gòu)會(huì)發(fā)生最大振幅,如圖4(a)所示.剪切來(lái)流工況下(k=0.05),頻率比較小時(shí)(r=1.0,1.5)DC 在橫流向上的振幅值隨Ur的增大規(guī)則的增大,然而在r=2.0 工況下,Ymax/D在Ur≤7 時(shí)緩慢的增加,但當(dāng)Ur=8 時(shí),Ymax/D會(huì)從0.5(Ur=7)跳躍到1.47(Ur=8),DC 在此頻率比工況振幅啟動(dòng)開(kāi)關(guān)要晚于另外兩個(gè)頻率比工況,如圖4(b) 所示.由圖5 可觀(guān)察到,在r=1.5 工況下,DC 在Ur=7 與Ur=8 時(shí)尾部的負(fù)壓面積相當(dāng),如圖5(b)所示,使得DC 在橫流向上的振幅相同,出現(xiàn)了“平臺(tái)期”.然而,在r=2.0 工況下,DC 在Ur=7 時(shí)尾部負(fù)壓區(qū)面積接近于0,當(dāng)折減速度增長(zhǎng)到8 時(shí),其尾部負(fù)壓區(qū)面積增大并超過(guò)了其他頻率比工況,導(dǎo)致結(jié)構(gòu)大振幅現(xiàn)象會(huì)延時(shí)出現(xiàn).當(dāng)剪切率進(jìn)一步增大到k=0.1 時(shí),UC 在小頻率比工況下(r=1.0,1.5)橫流向振幅被完全鎖定的折減速度區(qū)間有所擴(kuò)大(Ur=6~ 8),如圖4(c)所示.此時(shí)DC 在r=1.5 工況下無(wú)論是橫流向振幅還是順流向振幅都較平緩.當(dāng)r=1.0 與r=2.0 時(shí),其在橫流向上的振幅值在Ur≤7 時(shí)都較平緩,然而在Ur=8 處會(huì)突然跳躍到一個(gè)極大值,這與k=0.05 工況相似,但不同的是,在k=0.1 工況下,DC 在r=1.0 時(shí)Xmax/D最大會(huì)增長(zhǎng)到0.5(Ur=8),結(jié)構(gòu)順流方向振動(dòng)會(huì)達(dá)到最大值.

    圖4 不同頻率比與剪切率工況下,串列雙圓柱雙自由度最大振幅隨折減速度的變化(L/D=3.5)Fig.4 The variation of the maximum vibrating amplitude of tandem double cylinders with reduced velocity under different frequency ratios and shear rates in two degrees of freedom (L/D=3.5)

    圖5 不同頻率比工況下,流場(chǎng)無(wú)量綱時(shí)均流場(chǎng)壓力圖(k=0.05)Fig.5 Dimensionless time-averaged pressure diagram of flow field at different frequency ratios (k=0.05)

    3.1.3L/D=5.5

    當(dāng)間距比進(jìn)一步增大到5.5 時(shí),兩圓柱之間的間距超過(guò)臨界間距比,此時(shí)UC 的振動(dòng)情況與單圓柱工況相似,而DC 受到上游圓柱尾流的影響明顯增強(qiáng)[12].圖6 給出了L/D=5.5 時(shí)不同剪切來(lái)流工況下UC 與DC 雙自由度振幅變化情況.由圖可觀(guān)察到UC 較DC 振幅變化更加規(guī)律,其振幅變化主要集中在4<Ur<9 范圍內(nèi).隨著剪切率的增大,不同固有頻率比對(duì)DC 影響較大,對(duì)UC 影響較小.這是由于對(duì)于串列雙圓柱繞流,上游圓柱的渦激振動(dòng)是一種共振現(xiàn)象,然而下游圓柱受到共振作用的同時(shí),還受到UC 尾流與圓柱的相互作用,隨著剪切率與間距比的增大,尾流與圓柱的相互作用逐漸起到主導(dǎo)作用,使得不同頻率比對(duì)圓柱振幅的影響越發(fā)明顯.Ur≤7時(shí),DC 在r=2.0 工況下的Ymax/D值一直小于r=1.5 與r=1.0 工況.但是當(dāng)Ur=8 時(shí),r=2.0 工況下的橫流向最大振幅值會(huì)超過(guò)另兩個(gè)工況取得最大值,又會(huì)在Ur=9 時(shí)迅速減小.與均勻來(lái)流工況相比,DC 在剪切來(lái)流中振幅變化范圍有所擴(kuò)大,在UC 尾流的影響下,DC 出現(xiàn)了次峰,且剪切率越大,UC 尾流與DC 相互作用越強(qiáng),次峰越明顯,如圖6(b)與圖6(c)所示.相較于k=0.05 工況,剪切率增大到k=0.1 時(shí),DC 在順流向的幅值變化呈波動(dòng)狀,且頻率比越小波動(dòng)越明顯.DC 在橫流向上,r=1.0 與r=2.0 工況均會(huì)在Ur=8 處取得極大值,然而r=1.5 工況則處于平穩(wěn)的波動(dòng)狀,沒(méi)有明顯的峰值,與L/D=3.5 工況類(lèi)似,如圖4(c)所示.由于雙圓柱體結(jié)構(gòu)間距較大,UC 受到DC 的影響較小,UC 尾渦脫落完整,DC 上行渦拉長(zhǎng)在遠(yuǎn)尾流區(qū)脫落,如圖7 所示.而下行渦由于DC 的劇烈振動(dòng)受到抑制,在整個(gè)鎖定區(qū)間內(nèi),圓柱脫渦情況較為接近,使得此工況下圓柱振幅大小接近.

    圖6 不同頻率比與剪切率工況下,串列雙圓柱雙自由度最大振幅隨折減速度的變化(L/D=5.5)Fig.6 The variation of the maximum vibrating amplitude of tandem double cylinders with reduced velocity under different frequency ratios and shear rates in two degrees of freedom (L/D=5.5)

    圖7 r=1.5,k=0.1 時(shí)不同折減速度下的同一時(shí)刻流場(chǎng)瞬時(shí)渦量圖Fig.7 Instantaneous vorticity diagram of flow field at the same time at r=1.5 and k=0.1

    圖7 r=1.5,k=0.1 時(shí)不同折減速度下的同一時(shí)刻流場(chǎng)瞬時(shí)渦量圖 (續(xù))Fig.7 Instantaneous vorticity diagram of flow field at the same time at r=1.5 and k=0.1 (continued)

    3.2 相位特性變化

    圓柱渦激振動(dòng)中升力與位移的相位差變化體現(xiàn)了流體力做功的正負(fù).當(dāng)結(jié)構(gòu)的振動(dòng)頻率接近渦脫落頻率時(shí)會(huì)發(fā)生同步現(xiàn)象,升力幅值突然變大并且其與結(jié)構(gòu)位移間的相位會(huì)有180°的跳躍[29-31].

    3.2.1L/D=2.5

    當(dāng)L/D=2.5 時(shí),雙圓柱體結(jié)構(gòu)升力與位移之間的相位差在折減速度Ur=3~ 6 時(shí)接近于零,如圖8 所示,說(shuō)明此時(shí)圓柱結(jié)構(gòu)在一個(gè)振動(dòng)周期內(nèi)能量?jī)魝鬏斄繛榱?隨著折減速度進(jìn)一步增大,在Ur=6~ 9 區(qū)間內(nèi)升力與位移之間的相位差完成了從同相到反相的轉(zhuǎn)變,相位從0°跳躍到180°,能量傳遞在此區(qū)間內(nèi)變得復(fù)雜起來(lái),并且這個(gè)轉(zhuǎn)變區(qū)間正對(duì)應(yīng)著圓柱振幅鎖定區(qū)間,由此說(shuō)明相位轉(zhuǎn)變與結(jié)構(gòu)振動(dòng)幅度有著密切的關(guān)系,如圖2 所示.整體來(lái)說(shuō),剪切率的改變對(duì)升力與位移相位差的影響不明顯,UC 與DC 僅在相位轉(zhuǎn)變區(qū)間(Ur=6~ 9)較為雜亂而在其它區(qū)間內(nèi)都較為規(guī)則.當(dāng)間距比較小時(shí),UC 與DC 結(jié)合類(lèi)似一個(gè)整體,在來(lái)流作用下,兩者之間的相位差呈現(xiàn)同步變化.另一方面,頻率比的變化對(duì)UC 與DC 的影響主要在于:在相位轉(zhuǎn)變區(qū)間,頻率比越小,結(jié)構(gòu)升力與位移之間的相位差完成從同相到反相的轉(zhuǎn)變速度越快.然而,頻率比越大在鎖定區(qū)間圓柱振幅越大,如圖2 所示.由此可以總結(jié)出:隨頻率比的增大會(huì)逐漸加強(qiáng)對(duì)結(jié)構(gòu)渦激振動(dòng)效應(yīng)的影響,能量從流體傳遞到柱體的速度越緩慢,且過(guò)程會(huì)更復(fù)雜.

    圖8 在不同剪切率與頻率比下,串列雙圓柱升力與橫流向位移相位差變化圖(L/D=2.5)Fig.8 Diagram of phase difference between lift and transverse displacement of double cylinders in series at different shear rate and frequency ratio (L/D=2.5)

    以L(fǎng)/D=2.5,k=0.0,r=2.0 工況為例,當(dāng)Ur=6,8,10 時(shí),UC 上升力與位移間相位差分別處于同相、相位轉(zhuǎn)變以及反相區(qū)間.圖9 為3 個(gè)折減速度工況下升力時(shí)程曲線(xiàn)與位移時(shí)程曲線(xiàn)圖,由圖9(a)可知,當(dāng)Ur=6 時(shí),升力與位移同相變化,升力時(shí)程曲線(xiàn)變化幅值明顯大于位移時(shí)程曲線(xiàn),且橫流向位移時(shí)程曲線(xiàn)變化呈現(xiàn)標(biāo)準(zhǔn)的正弦曲線(xiàn)而升力位移時(shí)程曲線(xiàn)雖然同為正弦曲線(xiàn),但其每一個(gè)峰值由一大一小兩個(gè)次峰組成.當(dāng)Ur=8 時(shí),升力與位移處于相位轉(zhuǎn)變的復(fù)雜期,升力時(shí)程曲線(xiàn)呈現(xiàn)明顯的調(diào)制振幅,橫流向時(shí)程曲線(xiàn)雖為調(diào)制振幅但在各個(gè)瞬時(shí)階段總的振幅值相差不大,且位移幅值變化區(qū)間要大于升力時(shí)程曲線(xiàn),如圖9(b)所示.當(dāng)Ur=10 時(shí),升力與位移完成了相位轉(zhuǎn)變,升力時(shí)程曲線(xiàn)與位移時(shí)程曲線(xiàn)反相且均會(huì)呈現(xiàn)規(guī)則的正弦曲線(xiàn),升力變化幅值區(qū)間再次大于橫流向位移變化幅值區(qū)間.由此可以總結(jié)出:升力與位移時(shí)程曲線(xiàn)在同相與反相區(qū)間會(huì)呈規(guī)則的正弦曲線(xiàn),然而在相位轉(zhuǎn)變期間,能量傳遞紊亂,升力與位移時(shí)程曲線(xiàn)均會(huì)發(fā)生調(diào)制振動(dòng).

    圖9 在Ur=6,8,10 時(shí),上游圓柱升力與橫流向位移時(shí)程曲線(xiàn)圖(k=0.0,r=2.0)Fig.9 The time history curve of lift force and transverse displacement of upstream cylinder at Ur=6,8 and 10 (k=0.0,r=2.0)

    3.2.2L/D=3.5

    當(dāng)間距比L/D=3.5 時(shí),隨著兩圓柱之間間隙的增大,上游圓柱的尾流將對(duì)下游圓柱造成一定的影響.如圖10 所示,在均勻來(lái)流工況下,UC 與DC 上升力與位移相位差的轉(zhuǎn)變與L/D=2.5 工況相似,都主要在“鎖定區(qū)間”完成相位的轉(zhuǎn)變.剪切來(lái)流下UC 在Ur=4 處出現(xiàn)了“弱鎖定”現(xiàn)象,特殊的是隨著剪切率增大到0.1 時(shí),UC 與DC 在相位轉(zhuǎn)變區(qū)間(Ur=7~ 8)出現(xiàn)了“平臺(tái)期”,在這個(gè)“平臺(tái)期”內(nèi),隨折減速度的變大,升力與位移相位差不發(fā)生改變,流體與柱體之間能量傳遞相同.如圖10(c)所示,“平臺(tái)期”僅出現(xiàn)在r=1.5 與r=2.0 工況中,UC 在此兩個(gè)頻率比的相位差接近而DC 隨頻率比的增大升力與位移的相位差增大.通過(guò)對(duì)“平臺(tái)期”附近折減速度工況進(jìn)行加算發(fā)現(xiàn)從Ur=6 增長(zhǎng)到Ur=7 的過(guò)程中,各工況呈規(guī)律的線(xiàn)性增長(zhǎng),而從Ur=8 增長(zhǎng)到Ur=9的過(guò)程中,會(huì)在Ur=8 處出現(xiàn)跳躍增長(zhǎng).較為特殊的是,在小頻率比r=1.0 時(shí),UC 相位差會(huì)從40°跳躍到120°,而DC 相位差的變化較小,Ur=7 與Ur=8兩個(gè)工況的相位差幾乎一致.

    圖10 在不同剪切率與頻率比下,串列雙圓柱升力與橫流向位移相位差變化圖(L/D=3.5)Fig.10 Diagram of phase difference between lift and transverse displacement of double cylinders in series at different shear rate and frequency ratio (L/D=3.5)

    以L(fǎng)/D=3.5,k=0.1 工況為例,DC 橫流向振幅在Ur=8 時(shí)取得最大值,結(jié)構(gòu)進(jìn)入“鎖定狀態(tài)”發(fā)生共振現(xiàn)象.圖11 給出了不同頻率比工況下DC 在共振狀態(tài)的升力時(shí)程曲線(xiàn)與位移時(shí)程曲線(xiàn).當(dāng)r=1.0 時(shí),升力與橫流向位移時(shí)程曲線(xiàn)均為調(diào)制振幅,且橫流向位移時(shí)程曲線(xiàn)呈現(xiàn)為更為特殊的“拍”頻.當(dāng)升力曲線(xiàn)幅值達(dá)到最大值時(shí),位移曲線(xiàn)幅值則達(dá)到最小值,兩者之間會(huì)發(fā)生錯(cuò)位振動(dòng).當(dāng)r=1.5 時(shí),升力與位移時(shí)程曲線(xiàn)調(diào)制幅度變小,時(shí)程曲線(xiàn)變化情況越發(fā)規(guī)律,位移時(shí)程曲線(xiàn)幅度明顯要大于升力時(shí)程曲線(xiàn).當(dāng)頻率比進(jìn)一步增大到2.0 時(shí),橫流向位移時(shí)程曲線(xiàn)已恢復(fù)成規(guī)律的正弦曲線(xiàn),而升力時(shí)程曲線(xiàn)則呈現(xiàn)一個(gè)大周期包含3 個(gè)小周期的特殊調(diào)制振幅.

    圖11 在不同頻率比工況下,下游圓柱升力與橫流向位移時(shí)程曲線(xiàn)圖(k=0.1,Ur=8)Fig.11 The time history curve of lift force and transverse displacement of downstream cylinder at Ur=6,8 and 10 (k=0.1, Ur=8)

    3.2.3L/D=5.5

    當(dāng)間距比進(jìn)一步增大到5.5 時(shí),此時(shí)兩圓柱之間的距離超過(guò)了臨界間距比,上游圓柱與來(lái)流之間的相互作用類(lèi)似于單圓柱體,而下游圓柱既受到來(lái)流作用,又受到上游圓柱尾流的干擾,使得其橫流向位移與升力相位差隨折減速度的變化越發(fā)復(fù)雜.由圖12 可知剪切率的不同對(duì)UC 影響較小而對(duì)DC 影響較大.均勻來(lái)流工況下,如圖12(a) 所示,UC 上升力與位移間的相位差變化較為規(guī)律,主要在“鎖定區(qū)間”完成了相位的轉(zhuǎn)變,且頻率比越大,轉(zhuǎn)變的速度越緩慢.DC 的不同之處是:當(dāng)Ur=8 時(shí),DC 升力與位移相位差增長(zhǎng)到150°附近,而在Ur=9處相位差跌落到120°出現(xiàn)了能量的“反哺”,隨著折減速度的增大又再次回到180°.在剪切來(lái)流當(dāng)中,當(dāng)折減速度較小時(shí)(Ur=3),上下兩圓柱存在30°的初始相位差,隨著Ur的增大,升力與位移間的相位差先減小到0°隨后在“鎖定區(qū)間”(Ur=6~ 8)完成從同相到反相的轉(zhuǎn)變.特別的是,當(dāng)DC 退出“鎖定區(qū)間”后,在Ur=9 處,升力與位移間的相位差會(huì)從180°減小到120°左右,最后再隨折減速度的增大又重新回到180°.從圖12(b)與圖12(c)可以觀(guān)察到,剪切率越大(k=0.1),下游圓柱相位差重新回到180°的折減速度越大.當(dāng)k=0.05 時(shí),r=1.5 與r=2.0 工況下的DC 相位差變化同步,不同的是r=1.0 工況下的相位差在小幅度上升后又會(huì)在Ur=11 處降低到120°左右,最后在Ur=16 處徹底完成相位的轉(zhuǎn)變.當(dāng)剪切率增大到1.0 時(shí),UC 與DC 僅在r=1.5 頻率比工況下出現(xiàn)了“平臺(tái)期”.對(duì)于上游圓柱而言,當(dāng)Ur≥9 時(shí),3 個(gè)頻率比工況均已完成相位的轉(zhuǎn)變,但特殊的是r=1.0 工況較另兩個(gè)頻率比工況要更加不穩(wěn)定.對(duì)于下游圓柱,當(dāng)Ur=8 時(shí),3 個(gè)頻率比工況相位差值較大,r=1.0 與r=2.0 工況均是在到達(dá)一個(gè)極大值后在Ur=9 時(shí)相位差值回落且隨著折減速度的增大最后在較大折減速度下重新上升到180°.然而,對(duì)于r=1.5 工況沒(méi)有出現(xiàn)相位差回落現(xiàn)象,其升力與位移間的相位差隨著折減速度的增大緩慢增長(zhǎng),并最終在Ur=18 時(shí)徹底完成相位的轉(zhuǎn)變.

    圖12 在不同剪切率與頻率比下,串列雙圓柱升力與橫流向位移相位差變化圖(L/D=5.5)Fig.12 Diagram of phase difference between lift and transverse displacement of double cylinders in series at different shear rate and frequency ratio (L/D=5.5)

    由圖12(c) 可觀(guān)察到當(dāng)圓柱退出鎖定區(qū)間后,UC 在Ur=8 處完成相位轉(zhuǎn)變后,隨著折減速度的增大,將會(huì)一直維持反相變化.圖13 給出了L/D=5.5,k=0.1,r=1.0,Ur=14 工況UC 與DC 的升力時(shí)程曲線(xiàn)與橫流向位移時(shí)程曲線(xiàn)圖.由圖13(a) 可知,UC 的橫流向位移時(shí)程曲線(xiàn)為一個(gè)大周期包含兩個(gè)小周期的調(diào)制振幅,而升力時(shí)程曲線(xiàn)變化類(lèi)似于正弦曲線(xiàn),但每臨近兩個(gè)峰值幅值有細(xì)小的差別.另一方面,DC 升力時(shí)程曲線(xiàn)呈現(xiàn)特殊的調(diào)制振幅,波峰變化相同,而波谷則為一低一高兩種波況,且兩者之間最小振幅值為接近1.7 倍大小的關(guān)系,如圖13(b)所示.

    圖13 Ur=14 時(shí),串列雙圓柱升力與橫流向位移時(shí)程曲線(xiàn)圖(k=0.1,r=1.0)Fig.13 The time history curve of lift force and transverse displacement of upstream cylinder and downstream cylinder at Ur=14 (k=0.1, r=1.0)

    3.3 流體力功率譜密度分析

    在3.2 節(jié)中通過(guò)對(duì)串列雙圓柱的升力與橫流向位移相位差變化情況進(jìn)行對(duì)比分析,探究了不同工況下圓柱結(jié)構(gòu)與流體之間能量轉(zhuǎn)換的關(guān)系.本節(jié)將以r=1.5 工況為例進(jìn)一步對(duì)下游圓柱流體力系數(shù)功率譜密度(PSD)在不同間距比以及剪切率工況下的主峰幅值、頻譜成分及波動(dòng)性進(jìn)行分析,深入探究結(jié)構(gòu)能量變化情況及內(nèi)在力學(xué)機(jī)理.

    如圖14 所示,在均勻來(lái)流工況下,當(dāng)Ur=4 時(shí),不同工況下的升阻力PSD 曲線(xiàn)均較為光滑,主頻能量值較小,圓柱振動(dòng)微弱.當(dāng)折減速度增長(zhǎng)到6 時(shí),在小間距比工況下(L/D≤3.5)流體力PSD 曲線(xiàn)增大至3 階,升力PSD 曲線(xiàn)包圍的面積有效增大使得DC 獲得更多的動(dòng)能,在橫流向開(kāi)始產(chǎn)生較大的振幅.然而,在L/D=5.5 工況下,流體力系數(shù)PSD 曲線(xiàn)次峰及雜頻占據(jù)的能量值增多,下游圓柱運(yùn)動(dòng)更加復(fù)雜.當(dāng)Ur為7~ 8 時(shí),此時(shí)下游圓柱位于鎖定區(qū)間,PSD 曲線(xiàn)出現(xiàn)多階同等量級(jí)的頻率峰值,波動(dòng)性增強(qiáng),能量由流體大量流向圓柱體,DC 劇烈振動(dòng).當(dāng)Ur增大到9 時(shí),在L/D≤3.5 工況下,下游圓柱流體力PSD 曲線(xiàn)波動(dòng)性減弱,能量已基本完成從流體傳入結(jié)構(gòu)物的轉(zhuǎn)變.然而在L/D=5.5 工況下,流體力PSD 曲線(xiàn)出現(xiàn)了多個(gè)雜頻,且能量均勻分布在各雜頻上,使得能量出現(xiàn)“反哺”現(xiàn)象,能量會(huì)從結(jié)構(gòu)傳回到流體中,如圖12 所示.

    圖14 不同間距比工況下,下游圓柱升力系數(shù)功率譜密度(紅線(xiàn))與阻力系數(shù)功率譜密度(黑線(xiàn)) (k=0.0,r=1.5)Fig.14 The power spectral density of lift coefficient (red) and drag coefficient (black) of downstream cylinder under different spacing ratios (k=0.0,r=1.5)

    與均勻來(lái)流相比,剪切來(lái)流工況下,各工況流體力PSD 曲線(xiàn)波動(dòng)性增強(qiáng),所包含的多階頻率峰值增多,能量變化劇烈相對(duì)應(yīng)的下游圓柱振動(dòng)情況越發(fā)劇烈,如圖15 所示.在Ur=4 工況下,升力PSD 曲線(xiàn)與阻力PSD 曲線(xiàn)主峰頻率重合,使得能量大量由流體傳入圓柱體結(jié)構(gòu),在升力與橫流向相位差圖中出現(xiàn)了“弱鎖定”現(xiàn)象.隨折減速度的增大(Ur≤8),升力系數(shù)PSD 曲線(xiàn)出現(xiàn)的次峰數(shù)目增多且構(gòu)成整體呈下降趨勢(shì)的曲線(xiàn),包含的總能量較多,DC 逐漸進(jìn)入鎖定狀態(tài)并振動(dòng)劇烈.值得注意的是,當(dāng)L/D≥3.5 時(shí),如圖15(b)與15(c)所示,Ur=7 與Ur=8 工況下流體力系數(shù)PSD 曲線(xiàn)出現(xiàn)的頻率峰值數(shù)目、整體曲線(xiàn)變化趨勢(shì)以及各階頻率峰值所對(duì)應(yīng)的能量值都相同.從圖16 觀(guān)察到,此時(shí)兩個(gè)折減速度工況在一個(gè)周期內(nèi)的能量變化曲線(xiàn)其高、低峰值以及曲線(xiàn)變化形態(tài)都相同,這就意味著折減速度從7 增至8 的過(guò)程中圓柱能量變化相同,這可能是升力與橫流向位移相位差出現(xiàn)“平臺(tái)期”的原因,如圖10(c)與圖12(c)所示.

    圖15 不同間距比工況下,下游圓柱升力系數(shù)功率譜密度(紅線(xiàn))與阻力系數(shù)功率譜密度(黑線(xiàn)) (k=0.1,r=1.5)Fig.15 The power spectral density of lift coefficient (red) and drag coefficient (black) of downstream cylinder under different spacing ratios(k=0.1,r=1.5)

    圖16 在Ur=7 與Ur=8 工況下串列雙圓柱的能量曲線(xiàn)(L/D=3.5,k=0.1,r=1.5)Fig.16 Energy curve of tandem double cylinders under the of Ur=7 and Ur=8 cases (L/D=3.5,k=0.1,r=1.5)

    圖16 在Ur=7 與Ur=8 工況下串列雙圓柱的能量曲線(xiàn)(L/D=3.5,k=0.1,r=1.5) (續(xù))Fig.16 Energy curve of tandem double cylinders under the of Ur=7 and Ur=8 cases (L/D=3.5,k=0.1,r=1.5) (continued)

    4 結(jié)論

    基于四步半隱式特征線(xiàn)分裂算子有限元方法,對(duì)串列布置雙圓柱雙自由度結(jié)構(gòu)體系渦激振動(dòng)問(wèn)題進(jìn)行了數(shù)值模擬,分析了間距比、剪切率、頻率比以及折減速度四個(gè)關(guān)鍵影響因素的變化對(duì)結(jié)構(gòu)振動(dòng)響應(yīng)及流場(chǎng)的影響,主要結(jié)論如下.

    (1)UC 的共振區(qū)間明顯寬于DC 的共振區(qū)間,UC 在兩個(gè)自由度方向達(dá)到最大值的折減速度不同,而DC 基本同步.由于來(lái)流與間隙流的干擾,使得UC 與DC“鎖定區(qū)間”的折減速度范圍存在差值,UC 較DC 更早進(jìn)入也更早退出共振狀態(tài).固有頻率比與剪切來(lái)流的變化對(duì)DC 影響較大,對(duì)UC 影響較小.

    (2)兩圓柱主要在鎖定區(qū)間完成相位的轉(zhuǎn)變.頻率比越大對(duì)圓柱渦激振動(dòng)影響越大,能量從流體傳遞到柱體的速度越緩慢與復(fù)雜,圓柱完成從同相到反相的速度越緩慢.特別的是,在L/D≥3.5,r≥1.5的剪切來(lái)流工況下會(huì)出現(xiàn)相位差的“平臺(tái)期”.此“平臺(tái)期”出現(xiàn)在鎖定區(qū)間Ur=7 與Ur=8 工況下,此時(shí)隨著折減速度的增大,升力與位移之間的相位差不變.當(dāng)間距比超過(guò)臨界值時(shí),對(duì)DC 相位轉(zhuǎn)變影響較大,下游圓柱在完成相位轉(zhuǎn)變后會(huì)隨著折減速度的增大出現(xiàn)相位差回落現(xiàn)象,并且剪切率越大,回落的折減速度區(qū)間越廣.

    (3)在均勻來(lái)流工況下,升阻力PSD 曲線(xiàn)頻率呈兩倍關(guān)系,然而在剪切來(lái)流工況下,則基本重合.當(dāng)間距比較大時(shí),折減速度越大,流體力PSD 曲線(xiàn)出現(xiàn)的雜頻越多,由于能量均勻分布在各雜頻上,使得能量出現(xiàn)“反哺”現(xiàn)象.在剪切來(lái)流工況下,當(dāng)L/D≥3.5 時(shí),Ur=7 與Ur=8 工況中流體力系數(shù)PSD 曲線(xiàn)高度相似,同時(shí)兩工況的能量變化也相同,這可能是在升力與位移相位差圖中出現(xiàn)“平臺(tái)期”的原因.

    猜你喜歡
    振動(dòng)
    振動(dòng)的思考
    某調(diào)相機(jī)振動(dòng)異常診斷分析與處理
    振動(dòng)與頻率
    This “Singing Highway”plays music
    具非線(xiàn)性中立項(xiàng)的廣義Emden-Fowler微分方程的振動(dòng)性
    中立型Emden-Fowler微分方程的振動(dòng)性
    基于ANSYS的高速艇艉軸架軸系振動(dòng)響應(yīng)分析
    船海工程(2015年4期)2016-01-05 15:53:26
    主回路泵致聲振動(dòng)分析
    UF6振動(dòng)激發(fā)態(tài)分子的振動(dòng)-振動(dòng)馳豫
    帶有強(qiáng)迫項(xiàng)的高階差分方程解的振動(dòng)性
    中国国产av一级| 国产精品99久久99久久久不卡| 国产日韩欧美视频二区| 欧美变态另类bdsm刘玥| 满18在线观看网站| 多毛熟女@视频| 亚洲欧美一区二区三区久久| 97在线人人人人妻| 欧美日韩综合久久久久久| 久久九九热精品免费| 美女扒开内裤让男人捅视频| 欧美国产精品va在线观看不卡| 日韩制服丝袜自拍偷拍| 中文字幕人妻丝袜一区二区| 亚洲av日韩精品久久久久久密 | 91精品伊人久久大香线蕉| 欧美日韩黄片免| 久久毛片免费看一区二区三区| 亚洲国产精品一区二区三区在线| av又黄又爽大尺度在线免费看| 99国产精品99久久久久| 国产成人欧美| 亚洲中文av在线| 男人操女人黄网站| 久久 成人 亚洲| 丝袜在线中文字幕| 视频区欧美日本亚洲| 久9热在线精品视频| 国产高清videossex| 亚洲国产精品成人久久小说| av视频免费观看在线观看| 欧美精品一区二区大全| xxx大片免费视频| 日韩免费高清中文字幕av| 99久久人妻综合| 一二三四社区在线视频社区8| 国产精品国产三级专区第一集| 丰满人妻熟妇乱又伦精品不卡| 一边摸一边抽搐一进一出视频| 波多野结衣av一区二区av| 成在线人永久免费视频| e午夜精品久久久久久久| 午夜日韩欧美国产| 美女视频免费永久观看网站| 新久久久久国产一级毛片| avwww免费| 纵有疾风起免费观看全集完整版| 黑人猛操日本美女一级片| 欧美精品亚洲一区二区| 午夜福利,免费看| 一区二区三区精品91| 国产视频首页在线观看| 国产成人欧美| 国产av精品麻豆| 国产成人91sexporn| 精品卡一卡二卡四卡免费| 自线自在国产av| 看十八女毛片水多多多| 十分钟在线观看高清视频www| 欧美乱码精品一区二区三区| 国产精品久久久av美女十八| 国产成人精品无人区| 色综合欧美亚洲国产小说| 五月天丁香电影| 男人添女人高潮全过程视频| 超色免费av| 在线观看国产h片| 久久人人爽人人片av| 国产又色又爽无遮挡免| 日本午夜av视频| 亚洲精品av麻豆狂野| 少妇粗大呻吟视频| 日韩熟女老妇一区二区性免费视频| 欧美日韩福利视频一区二区| 国产在线一区二区三区精| 免费高清在线观看日韩| 777米奇影视久久| 欧美精品啪啪一区二区三区 | 黄片小视频在线播放| 欧美亚洲 丝袜 人妻 在线| 精品第一国产精品| 好男人电影高清在线观看| 亚洲成色77777| 99久久精品国产亚洲精品| 天天躁日日躁夜夜躁夜夜| 精品亚洲成国产av| 丝袜美腿诱惑在线| 亚洲天堂av无毛| 亚洲成人国产一区在线观看 | 国产av精品麻豆| 国产福利在线免费观看视频| 在线观看免费日韩欧美大片| 欧美黄色淫秽网站| 成人18禁高潮啪啪吃奶动态图| 亚洲精品av麻豆狂野| 日本欧美国产在线视频| 亚洲精品中文字幕在线视频| 免费黄频网站在线观看国产| 国产又色又爽无遮挡免| 国产男女超爽视频在线观看| 天堂中文最新版在线下载| 国产黄频视频在线观看| 精品欧美一区二区三区在线| 国产99久久九九免费精品| 秋霞在线观看毛片| 国产精品.久久久| 极品人妻少妇av视频| 精品久久久精品久久久| 国产午夜精品一二区理论片| 日本一区二区免费在线视频| 欧美精品啪啪一区二区三区 | 亚洲久久久国产精品| 日本午夜av视频| 丝袜美足系列| 欧美在线黄色| 欧美精品一区二区免费开放| 女人爽到高潮嗷嗷叫在线视频| 搡老岳熟女国产| 久久免费观看电影| 91精品伊人久久大香线蕉| 在线观看人妻少妇| a级毛片在线看网站| 青春草亚洲视频在线观看| 亚洲精品在线美女| 国产有黄有色有爽视频| 欧美精品av麻豆av| 黑人巨大精品欧美一区二区蜜桃| cao死你这个sao货| 啦啦啦在线免费观看视频4| 乱人伦中国视频| 男女之事视频高清在线观看 | 国产精品麻豆人妻色哟哟久久| 大型av网站在线播放| 黑人猛操日本美女一级片| 亚洲精品中文字幕在线视频| 亚洲欧美一区二区三区黑人| 黄色a级毛片大全视频| 自线自在国产av| 成年美女黄网站色视频大全免费| 久久热在线av| 国产亚洲av片在线观看秒播厂| 国产福利在线免费观看视频| 国精品久久久久久国模美| videos熟女内射| 久久免费观看电影| 午夜福利视频精品| 成人午夜精彩视频在线观看| 18禁黄网站禁片午夜丰满| 咕卡用的链子| a级片在线免费高清观看视频| 国产熟女午夜一区二区三区| 久久久久网色| 国产主播在线观看一区二区 | 99国产精品99久久久久| 亚洲精品国产av蜜桃| 熟女av电影| 亚洲av欧美aⅴ国产| 爱豆传媒免费全集在线观看| 韩国高清视频一区二区三区| 99精国产麻豆久久婷婷| 精品一区在线观看国产| av天堂在线播放| 亚洲国产欧美一区二区综合| 99国产精品99久久久久| 久久久欧美国产精品| 99精品久久久久人妻精品| 婷婷色综合www| 免费看不卡的av| bbb黄色大片| 日本欧美国产在线视频| 国产精品二区激情视频| 国产成人欧美在线观看 | 久久99一区二区三区| 在线av久久热| 日日夜夜操网爽| 欧美成人精品欧美一级黄| 18禁国产床啪视频网站| 久久久精品国产亚洲av高清涩受| 久热爱精品视频在线9| 在线天堂中文资源库| 成人国产av品久久久| 国产高清国产精品国产三级| 欧美+亚洲+日韩+国产| 亚洲七黄色美女视频| 久久热在线av| 男人添女人高潮全过程视频| 日本午夜av视频| 日本欧美国产在线视频| 黄色一级大片看看| 久久ye,这里只有精品| 丝瓜视频免费看黄片| 国产女主播在线喷水免费视频网站| 50天的宝宝边吃奶边哭怎么回事| 亚洲第一av免费看| 精品国产一区二区久久| 黄网站色视频无遮挡免费观看| 波野结衣二区三区在线| 91老司机精品| 亚洲精品久久成人aⅴ小说| 后天国语完整版免费观看| 欧美日韩国产mv在线观看视频| 另类亚洲欧美激情| 18禁裸乳无遮挡动漫免费视频| 国产精品一二三区在线看| 天天添夜夜摸| 人人妻人人添人人爽欧美一区卜| 国产又爽黄色视频| 日本一区二区免费在线视频| 午夜两性在线视频| 日日夜夜操网爽| 日韩视频在线欧美| 国产成人免费观看mmmm| 久久久久久久国产电影| 亚洲人成电影观看| 男人爽女人下面视频在线观看| 捣出白浆h1v1| 国产亚洲欧美精品永久| 久久精品国产a三级三级三级| 久久久久久久久久久久大奶| 国产成人精品久久二区二区免费| 成人午夜精彩视频在线观看| 久久久久视频综合| 成人国产一区最新在线观看 | 午夜视频精品福利| 亚洲一区二区三区欧美精品| 大片免费播放器 马上看| 国产精品国产av在线观看| 日本wwww免费看| 少妇 在线观看| 国产精品一区二区在线观看99| 国产熟女午夜一区二区三区| 亚洲国产精品一区三区| 美女国产高潮福利片在线看| 午夜精品国产一区二区电影| 亚洲,欧美,日韩| 丝袜人妻中文字幕| 青春草亚洲视频在线观看| 激情视频va一区二区三区| 尾随美女入室| 欧美人与性动交α欧美精品济南到| 看免费av毛片| 成年av动漫网址| 午夜福利免费观看在线| 老司机午夜十八禁免费视频| 久久精品国产综合久久久| 国产精品麻豆人妻色哟哟久久| 男女免费视频国产| 国产精品 欧美亚洲| 日韩一区二区三区影片| 视频区欧美日本亚洲| 国产高清国产精品国产三级| av欧美777| 久久这里只有精品19| 国产视频一区二区在线看| 国产熟女欧美一区二区| 亚洲国产欧美网| 免费在线观看完整版高清| 婷婷丁香在线五月| 男女午夜视频在线观看| 日本一区二区免费在线视频| 老司机靠b影院| 色精品久久人妻99蜜桃| a级毛片黄视频| 夫妻性生交免费视频一级片| 国产成人免费无遮挡视频| 欧美变态另类bdsm刘玥| 亚洲国产中文字幕在线视频| 久久国产精品男人的天堂亚洲| 在线观看免费午夜福利视频| 婷婷色av中文字幕| 超色免费av| 男人爽女人下面视频在线观看| 国产精品一区二区在线观看99| 久久人人爽av亚洲精品天堂| 国产不卡av网站在线观看| 久久久久国产一级毛片高清牌| 成年av动漫网址| 在线观看免费日韩欧美大片| 麻豆乱淫一区二区| 国产精品 国内视频| 咕卡用的链子| 久久国产精品大桥未久av| xxxhd国产人妻xxx| 亚洲国产精品999| 成人三级做爰电影| 亚洲av成人不卡在线观看播放网 | 人人澡人人妻人| av欧美777| 一本—道久久a久久精品蜜桃钙片| 女人被躁到高潮嗷嗷叫费观| 成年人黄色毛片网站| 手机成人av网站| 天天添夜夜摸| 9191精品国产免费久久| 成人免费观看视频高清| 一区二区三区四区激情视频| 一二三四社区在线视频社区8| av一本久久久久| 国产真人三级小视频在线观看| 少妇人妻 视频| 国产精品久久久久久精品古装| 色精品久久人妻99蜜桃| 国产精品免费视频内射| 国产精品 国内视频| 欧美性长视频在线观看| 狠狠婷婷综合久久久久久88av| 国产1区2区3区精品| 亚洲,欧美精品.| 夫妻午夜视频| 两个人免费观看高清视频| 久久久久国产一级毛片高清牌| 亚洲国产精品国产精品| 老司机午夜十八禁免费视频| 日本五十路高清| 亚洲一卡2卡3卡4卡5卡精品中文| 麻豆乱淫一区二区| 只有这里有精品99| 精品国产超薄肉色丝袜足j| 伊人亚洲综合成人网| 热99国产精品久久久久久7| 老司机影院毛片| 真人做人爱边吃奶动态| 日本黄色日本黄色录像| 校园人妻丝袜中文字幕| 中文字幕人妻丝袜一区二区| 校园人妻丝袜中文字幕| 国产一区二区三区综合在线观看| 日本av手机在线免费观看| 久久亚洲国产成人精品v| www.精华液| 巨乳人妻的诱惑在线观看| 十分钟在线观看高清视频www| 乱人伦中国视频| 黄色片一级片一级黄色片| 晚上一个人看的免费电影| 免费人妻精品一区二区三区视频| 欧美成人午夜精品| 18在线观看网站| 久久久久国产精品人妻一区二区| 亚洲一区二区三区欧美精品| 国产视频一区二区在线看| 一级,二级,三级黄色视频| 老司机在亚洲福利影院| 电影成人av| 亚洲欧美中文字幕日韩二区| 亚洲av欧美aⅴ国产| 热99国产精品久久久久久7| 肉色欧美久久久久久久蜜桃| 亚洲自偷自拍图片 自拍| 国产男女超爽视频在线观看| av电影中文网址| 99精国产麻豆久久婷婷| 又黄又粗又硬又大视频| 精品国产乱码久久久久久男人| 午夜91福利影院| 午夜福利在线免费观看网站| av电影中文网址| 国产熟女欧美一区二区| 国产精品熟女久久久久浪| 成年美女黄网站色视频大全免费| 欧美精品啪啪一区二区三区 | 成年人黄色毛片网站| 丰满饥渴人妻一区二区三| 欧美日韩成人在线一区二区| 亚洲精品中文字幕在线视频| 欧美日韩一级在线毛片| 国产深夜福利视频在线观看| 久久久久久亚洲精品国产蜜桃av| 精品亚洲成a人片在线观看| 日韩制服骚丝袜av| tube8黄色片| 国产成人欧美在线观看 | 人人澡人人妻人| 人人妻人人爽人人添夜夜欢视频| 亚洲精品久久午夜乱码| 九草在线视频观看| 亚洲成人手机| 国产亚洲精品久久久久5区| 亚洲自偷自拍图片 自拍| 成年美女黄网站色视频大全免费| 曰老女人黄片| 无遮挡黄片免费观看| 可以免费在线观看a视频的电影网站| 天天躁狠狠躁夜夜躁狠狠躁| 国产男人的电影天堂91| av视频免费观看在线观看| 伦理电影免费视频| 99re6热这里在线精品视频| 91老司机精品| 国产亚洲欧美在线一区二区| 亚洲国产最新在线播放| 国产伦人伦偷精品视频| 亚洲精品日本国产第一区| 99国产精品一区二区蜜桃av | 99香蕉大伊视频| 亚洲专区中文字幕在线| 亚洲av男天堂| 日韩 欧美 亚洲 中文字幕| 国产精品三级大全| 人妻一区二区av| 高潮久久久久久久久久久不卡| 男女床上黄色一级片免费看| 91精品伊人久久大香线蕉| 精品久久久久久久毛片微露脸 | 一二三四在线观看免费中文在| 一级,二级,三级黄色视频| 美女国产高潮福利片在线看| 老司机影院毛片| 欧美日韩视频高清一区二区三区二| 少妇精品久久久久久久| 青草久久国产| 亚洲欧美一区二区三区久久| 黑人欧美特级aaaaaa片| 午夜激情av网站| 成人亚洲精品一区在线观看| 国产高清videossex| 涩涩av久久男人的天堂| 欧美亚洲 丝袜 人妻 在线| 在线 av 中文字幕| 久久女婷五月综合色啪小说| 国产精品一区二区免费欧美 | 午夜福利视频在线观看免费| 免费观看a级毛片全部| 国产黄色视频一区二区在线观看| 性色av乱码一区二区三区2| 男男h啪啪无遮挡| 十八禁网站网址无遮挡| 国产成人啪精品午夜网站| 国产亚洲av片在线观看秒播厂| 99精品久久久久人妻精品| 国产视频一区二区在线看| cao死你这个sao货| 女性被躁到高潮视频| 高清黄色对白视频在线免费看| 亚洲精品一卡2卡三卡4卡5卡 | 久久久久久亚洲精品国产蜜桃av| 午夜福利乱码中文字幕| 国产一区二区 视频在线| 国产精品九九99| 国语对白做爰xxxⅹ性视频网站| 国产免费福利视频在线观看| 久久国产亚洲av麻豆专区| 1024视频免费在线观看| 久久精品久久久久久久性| 久9热在线精品视频| 99久久人妻综合| 久久久久视频综合| 亚洲精品一二三| 久久99热这里只频精品6学生| 亚洲成色77777| 热99国产精品久久久久久7| 久久久久久久久久久久大奶| 999精品在线视频| 精品久久久久久电影网| 亚洲国产看品久久| 黄片播放在线免费| 色网站视频免费| 美女国产高潮福利片在线看| 欧美精品一区二区免费开放| 午夜免费成人在线视频| 午夜久久久在线观看| 欧美黄色片欧美黄色片| 好男人视频免费观看在线| 满18在线观看网站| 另类精品久久| 亚洲伊人久久精品综合| 久久狼人影院| 校园人妻丝袜中文字幕| 中文字幕另类日韩欧美亚洲嫩草| 一级a爱视频在线免费观看| 国产精品成人在线| av线在线观看网站| 国产亚洲欧美精品永久| 精品少妇黑人巨大在线播放| 一本久久精品| 视频在线观看一区二区三区| 亚洲 国产 在线| 中文字幕色久视频| 熟女av电影| 亚洲中文日韩欧美视频| 亚洲av欧美aⅴ国产| 婷婷成人精品国产| 十分钟在线观看高清视频www| 人妻 亚洲 视频| 亚洲色图综合在线观看| 久久精品国产亚洲av高清一级| 欧美成人午夜精品| 一区二区av电影网| 亚洲国产看品久久| 一级毛片女人18水好多 | 在现免费观看毛片| 成年人黄色毛片网站| 亚洲少妇的诱惑av| 亚洲欧洲精品一区二区精品久久久| 精品人妻熟女毛片av久久网站| 一本久久精品| 嫁个100分男人电影在线观看 | 美女扒开内裤让男人捅视频| 捣出白浆h1v1| 欧美精品啪啪一区二区三区 | 国产成人91sexporn| 久久精品aⅴ一区二区三区四区| 丰满迷人的少妇在线观看| 人人妻人人澡人人爽人人夜夜| www日本在线高清视频| 国产黄频视频在线观看| 亚洲成人免费av在线播放| 欧美精品av麻豆av| 天天添夜夜摸| 麻豆乱淫一区二区| av片东京热男人的天堂| 国产欧美日韩一区二区三 | 中文乱码字字幕精品一区二区三区| 免费高清在线观看视频在线观看| 丝瓜视频免费看黄片| 男女边吃奶边做爰视频| 九草在线视频观看| 中文乱码字字幕精品一区二区三区| 国产一区二区三区av在线| av在线老鸭窝| 人人妻人人添人人爽欧美一区卜| 国产男人的电影天堂91| 一本一本久久a久久精品综合妖精| 91精品伊人久久大香线蕉| 成在线人永久免费视频| 亚洲欧美一区二区三区久久| 午夜福利,免费看| 天天操日日干夜夜撸| 国产片特级美女逼逼视频| 视频在线观看一区二区三区| 蜜桃国产av成人99| 日韩一卡2卡3卡4卡2021年| 精品亚洲成a人片在线观看| 色综合欧美亚洲国产小说| h视频一区二区三区| 青春草亚洲视频在线观看| 中文字幕人妻丝袜一区二区| 久久久久久久久免费视频了| 亚洲中文av在线| 男人爽女人下面视频在线观看| 美女主播在线视频| 国产主播在线观看一区二区 | 91精品三级在线观看| 欧美国产精品va在线观看不卡| 午夜av观看不卡| 亚洲欧美色中文字幕在线| 精品少妇久久久久久888优播| 黄频高清免费视频| 晚上一个人看的免费电影| 日本猛色少妇xxxxx猛交久久| 成人免费观看视频高清| 亚洲欧美日韩另类电影网站| 后天国语完整版免费观看| 亚洲专区中文字幕在线| 日韩一卡2卡3卡4卡2021年| 在线观看一区二区三区激情| 国产亚洲av片在线观看秒播厂| 别揉我奶头~嗯~啊~动态视频 | 欧美国产精品va在线观看不卡| 在线看a的网站| 国产精品99久久99久久久不卡| 国产成人影院久久av| 国产一级毛片在线| 久久久精品国产亚洲av高清涩受| 国产野战对白在线观看| 久久久久国产一级毛片高清牌| 在线观看免费视频网站a站| 在线观看国产h片| 国产欧美亚洲国产| 亚洲国产最新在线播放| 国产成人精品无人区| 一边摸一边做爽爽视频免费| 观看av在线不卡| 精品人妻1区二区| 中文精品一卡2卡3卡4更新| 80岁老熟妇乱子伦牲交| 国产一区二区 视频在线| 五月开心婷婷网| 又黄又粗又硬又大视频| a级毛片在线看网站| 蜜桃国产av成人99| 99久久精品国产亚洲精品| 国产日韩欧美在线精品| 9色porny在线观看| 99国产精品一区二区蜜桃av | 男女高潮啪啪啪动态图| 欧美日韩综合久久久久久| 亚洲精品国产区一区二| 午夜久久久在线观看| 国产av精品麻豆| 黄色 视频免费看| 一级毛片电影观看| 中文字幕另类日韩欧美亚洲嫩草| 少妇粗大呻吟视频| 伦理电影免费视频| 51午夜福利影视在线观看| 成年av动漫网址| 久久热在线av| 天天操日日干夜夜撸| 十八禁高潮呻吟视频| 黄色怎么调成土黄色| 男的添女的下面高潮视频| 男女午夜视频在线观看| 啦啦啦中文免费视频观看日本| 1024视频免费在线观看| 天天躁狠狠躁夜夜躁狠狠躁| 黄色毛片三级朝国网站| 嫁个100分男人电影在线观看 | 久久青草综合色| 久久久久久久久久久久大奶| 91成人精品电影| 国产成人免费无遮挡视频| 国产精品偷伦视频观看了| 日韩视频在线欧美| 久久精品国产综合久久久|