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

    串列布置三圓柱渦激振動頻譜特性研究1)

    2021-11-09 08:46:10涂佳黃譚瀟玲梁經(jīng)群
    力學(xué)學(xué)報 2021年6期
    關(guān)鍵詞:雷諾數(shù)升力振幅

    涂佳黃 胡 剛 譚瀟玲 梁經(jīng)群 張 平

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

    引言

    渦激振動現(xiàn)象廣泛存在于實際工程中,當(dāng)柱體結(jié)構(gòu)處于流場時,交替的渦脫落在結(jié)構(gòu)表面產(chǎn)生周期性的流體力,導(dǎo)致結(jié)構(gòu)產(chǎn)生振動.當(dāng)海洋立管[1-2]、橋梁結(jié)構(gòu)[3-4]、高層建筑[5-6]結(jié)構(gòu)振動頻率接近于自身固有頻率時,結(jié)構(gòu)會發(fā)生“鎖定”現(xiàn)象[7],誘發(fā)結(jié)構(gòu)產(chǎn)生較大的振動幅度,帶來災(zāi)難性后果.另一方面,強烈的結(jié)構(gòu)振動可以進(jìn)行能源轉(zhuǎn)換和利用[8-9].因此,有關(guān)渦激振動問題一直是研究者們所關(guān)注的熱點之一.

    目前,學(xué)者們對單圓柱渦激振動問題進(jìn)行了相關(guān)研究,取得了一系列成果[10-13].與單圓柱渦激振動問題相比,串列多圓柱渦激振動情況更加復(fù)雜,對實際工程的參考意義更大.Papaioannou 等[14]研究了不同間距比工況下串列雙圓柱渦激振動問題,發(fā)現(xiàn)間距比較小時,上游圓柱的鎖定區(qū)間范圍較單圓柱工況明顯擴(kuò)大且最大振幅增大,超過臨界間距比后,下游圓柱對上游圓柱振幅的影響較小.Prasanth 等[15]對Re=100 和L=5.5D工況串列雙圓柱進(jìn)行了數(shù)值模擬研究,發(fā)現(xiàn)下游圓柱的鎖定區(qū)間受到了上游圓柱的影響,明顯大于單圓柱渦激振動的現(xiàn)象.及春寧等[16]對低雷諾數(shù)下串列雙圓柱渦激振動中下游圓柱大振幅響應(yīng)的耦合機制進(jìn)行了全面分析,發(fā)現(xiàn)上游圓柱脫落漩渦激勵下游圓柱大振幅響應(yīng)的機理來自于脫落漩渦所誘發(fā)的位于下游圓柱的運動方向上低壓區(qū).Mysa 等[17]研究了L=(1.0~4.0)D,Re=100下串列雙圓柱渦激振動中對響應(yīng)起關(guān)鍵作用的因素,發(fā)現(xiàn)上游圓柱尾流與下游圓柱之間的耦合作用的強弱取決于橫流向流體力與圓柱位移之間的相位差,與下游圓柱運動同相的部分對其響應(yīng)起到了促進(jìn)作用.Mittal 等[18]通過數(shù)值模擬發(fā)現(xiàn)Re=100,L=5.5D時上游圓柱的振動響應(yīng)趨近于單圓柱渦激振動,而下游圓柱經(jīng)歷了大振幅振動.

    關(guān)于串列三圓柱渦激振動問題,圓柱體之間相互作用對振動影響很大,潛在的機理也值得更多的努力去研究.Igarashi 和Suzuki[19]通過試驗研究了L=(1.0~4.0)D,Re=10 900~39 200 的串列三圓柱繞流特性,根據(jù)下游圓柱上從上游圓柱分離的具有動態(tài)效應(yīng)的剪切層,將6 種尾流模態(tài)進(jìn)行了分類,并得到了發(fā)生雙穩(wěn)態(tài)現(xiàn)象的區(qū)域.Yu 等[20]通過研究發(fā)現(xiàn)三圓柱系統(tǒng)的振動響應(yīng)劇烈,橫流向最大振幅提高了25%.同時,結(jié)構(gòu)運動軌跡大多呈現(xiàn)為不規(guī)則形狀.Chen 等[21]對Re=100,L=(1.2~5.0)D工況下串列三圓柱進(jìn)行數(shù)值模擬研究,發(fā)現(xiàn)L=1.2D時,串列三圓柱均表現(xiàn)為馳振現(xiàn)象.指出馳振現(xiàn)象的3 個關(guān)鍵因素是平衡位置偏移、低頻振動以及漩渦脫落與圓柱之間的時機.Mahmoud 和Atef[22]對串列布置三圓柱系統(tǒng)的流致振動問題進(jìn)行了研究,著重分析了上中游兩靜止圓柱體間距比的變化對下游圓柱體動力響應(yīng)及其流場特性的影響.張志猛等[23]對Re=100 時上游圓柱固定串列三圓柱渦激振動進(jìn)行數(shù)值模擬研究,發(fā)現(xiàn)當(dāng)振幅較小時,上游圓柱的剪切層將三圓柱包裹在一起,尾流表現(xiàn)為單鈍體脫渦模式.當(dāng)振幅較大時,上游圓柱的旋渦重附著于下游圓柱上,尾流呈現(xiàn)出兩列并排的漩渦.Behara 等[24]對串列布置三圓柱進(jìn)行數(shù)值模擬分析,對于下游兩圓柱根據(jù)其振幅和頻率的響應(yīng),將約化速度劃分為3 個區(qū)間:初始激勵區(qū)、上鎖定區(qū)、下鎖定區(qū).涂佳黃等[25]對Re=100下剪切來流作用下串列三圓柱體雙自由度流致振動問題進(jìn)行了數(shù)值計算研究,隨固有頻率比的增大,上游圓柱順流向鎖定區(qū)間范圍會減小,而中下游圓柱雙向鎖定區(qū)間會擴(kuò)大.

    綜上所述,頻譜特性分析在多柱體結(jié)構(gòu)渦激振動問題中的研究討論較少,本文基于四步半隱式特征線分裂算子有限元方法,對均勻來流作用下串列布置三圓柱系統(tǒng)雙自由度渦激振動問題進(jìn)行了數(shù)值模擬,重點分析了圓柱體功率譜密度隨不同雷諾數(shù)、固有頻率比、約化速度的變化規(guī)律,對不同工況下功率譜密度的形態(tài)進(jìn)行對比,深入研究了功率譜密度對結(jié)構(gòu)動力響應(yīng)的影響并揭示了其內(nèi)在機理.為實際工程應(yīng)用提供參考依據(jù).

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

    1.1 流體控制方程

    基于ALE 方法下的不可壓縮黏性流體的N-S 方程的無量綱形式表達(dá)如下

    式中,xi為空間坐標(biāo);ui為流體速度;cj=uj-wj,cj為流體相對于網(wǎng)格移動速度的對流速度,wj為網(wǎng)格移動速度;p為流體壓力;t為時間;Re=U∞D(zhuǎn)/ν,Re為雷諾數(shù),D與U∞分別為特征長度尺寸與特征速度,ν為動力黏性系數(shù).

    1.2 固體控制方程

    彈性支撐結(jié)果運動體系可假設(shè)為之質(zhì)量-阻尼-彈簧系統(tǒng),考慮雙自由度運動的控制方程量綱歸一化的形式如下

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

    1.3 網(wǎng)格更新

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

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

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

    1.4 計算流程

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

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

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

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

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

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

    2 問題描述

    2.1 計算模型

    為了研究上游結(jié)構(gòu)尾流充分發(fā)展后對下游結(jié)構(gòu)振動響應(yīng)的影響,選取圓柱體結(jié)構(gòu)間距比Lx=5.5D.其他參數(shù)為:雷諾數(shù)Re=80,100,160,固有頻率比r=fn,x/fn,y=1.0,1.5,2.0,約化速度Ur=Ur,x=3~21,質(zhì)量比mr=2.0.計算域尺寸為[-30D,60D]×[-20D,20D],上游圓柱(upstream cylinder,UC)、中游圓柱(midstream cylinder,MC) 和下游圓柱(downstream cylinder,DC)的圓心位置分別為(-5.5D,0)、(0,0) 與(5.5D,0),本文計算域堵塞率B=5%,如圖1 所示.邊界條件設(shè)置如下,入口處采用速度入口:ux=U∞=1.0,uy=0;出口處采用壓力出口:p=0;上下邊界均為自由滑移邊界:?ux/?y=0,uy=0;圓柱表面均采用無滑移邊界:ux=uy=0.為了獲得較大的動力響應(yīng),不考慮結(jié)構(gòu)阻尼的影響.對于彈性支撐圓柱體結(jié)構(gòu),可簡化為質(zhì)量-彈簧系統(tǒng)模型,圓柱表面速度為

    圖1 計算模型和邊界條件Fig.1 Computational model and boundary conditions

    2.2 網(wǎng)格劃分

    流場計算域采用非結(jié)構(gòu)化三角形網(wǎng)格進(jìn)行劃分,其中圓柱附近及尾流區(qū)域進(jìn)行網(wǎng)格加密處理.表1為不同網(wǎng)格參數(shù)下串列三圓柱結(jié)構(gòu)渦激振動動力響應(yīng)特性的統(tǒng)計結(jié)果對比.由表1 可知:與粗網(wǎng)格GI相比,網(wǎng)格GII 的Xmax/D與Ymax/D的計算結(jié)果最大變化率分別為5.26%和4.00%;與密網(wǎng)格GIII 相比,其最大變化率分別下降為1.25%和0.85%.網(wǎng)格GII已滿足數(shù)值模擬結(jié)果網(wǎng)格收斂性要求和計算精度要求.綜上所述,本文所有算例采用的網(wǎng)格密度和分布情況與GII 類似,無量綱計算時間步長Δt=0.002.

    表1 網(wǎng)格獨立性驗證:串列布置三圓柱體流致振動計算結(jié)果(Re=100,r=1.0,Ur=6.0)Table 1 Grid independence test:The results for the three tandem circular cylinders at Re=100,r=1.0 and Ur=6.0

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

    本文分析了不同雷諾數(shù)、固有頻率比與約化速度工況下串列三圓柱體結(jié)構(gòu)體系雙自由度流致振動振幅特性、頻譜特性、流體力系數(shù)的變化情況.本文根據(jù)UC 的振動特性對約化速度進(jìn)行如下劃分:3 ≤Ur≤5(區(qū)間A);5 <Ur≤9(區(qū)間B);9 <Ur≤15(區(qū)間C1),Ur>15(區(qū)間C2).

    3.1 振幅特性

    由圖2 可知,在區(qū)間A 內(nèi),UC 在x和y兩個方向的振幅隨約化速度的增大而增大.在大振幅區(qū)間內(nèi),頻率比對UC 振幅的影響較大,具體表現(xiàn)為:當(dāng)r≤1.5 時,頻率比的變化對UC 的振動響應(yīng)的影響較小.各工況下橫流向振幅十分接近,在Ur=5 工況下達(dá)到最大振幅值后逐漸下降.頻率比較小時結(jié)構(gòu)剛度較大,圓柱順流向振動響應(yīng)微弱,x方向的振幅趨于0.當(dāng)r=2.0,Ur=5 時,UC 橫流向振動頻率fs,y=0.174(Re=80,100)和0.178(Re=160)接近其固有頻率,同時順流向振動頻率與固有頻率的比值在1 附近.此時UC 發(fā)生雙頻共振,其振動響應(yīng)顯著增強,順流向振幅遠(yuǎn)大于其他頻率比工況下的振幅,橫流向振幅也達(dá)到較大值.值得注意地是,當(dāng)Re=80時,在雙頻共振作用下的尾渦模式為2S 模式,漩渦強度逐漸減弱并在下游不遠(yuǎn)處會消失,如圖3(a)所示;Re=100 時的尾渦模式也呈2S 模式,但下游的漩渦重新卷起并脫落形成新的渦街;不同地是,Re=160時尾渦模式呈現(xiàn)P+S 模式,在較遠(yuǎn)的下游能保持較高的強度.雷諾數(shù)對UC 的振動響應(yīng)的影響較小,各工況下UC 的振幅變化規(guī)律類似.進(jìn)入C1區(qū)間后,UC 在x和y兩個方向的振幅逐漸趨于穩(wěn)定,x方向的振動響應(yīng)極其微弱.

    圖2 不同雷諾數(shù)與頻率比工況下,UC 在x 和y 兩個方向最大振幅隨約化速度的變化Fig.2 Variation of the maximum vibrating amplitudes of upstream cylinder with reduced velocity under different Reynolds number and frequency ratios

    圖3 Ur=5,r=2.0 時,流場瞬時渦量圖Fig.3 Instantaneous vorticity diagram of flow field at Ur=5 and r=2.0

    與UC 相比,雷諾數(shù)和頻率比對MC 的振幅的影響更加顯著.在區(qū)間A 內(nèi),MC 受到“弱鎖定”作用,如圖4 所示,Ur=4 時,x和y兩個方向的振幅取得第一個極大值,當(dāng)r=1.5 時順流向振幅接近MC 的順流向振幅最大值.特別地,當(dāng)r=2.0 時,MC 在x和y兩個方向的運動在區(qū)間A 內(nèi)的振幅很小,這種“屏蔽”現(xiàn)象可以用流場進(jìn)行解釋.如圖5 所示,當(dāng)Ur=4時,UC 在靠近MC 的方向產(chǎn)生漩渦,而漩渦脫落的位置與UC 的距離非常近,連續(xù)渦結(jié)構(gòu)的垂直距離較大,說明尾流發(fā)展過程中沒有渦撞擊的發(fā)生,漩渦直接從MC 和DC 的兩側(cè)通過,在UC 的下游不遠(yuǎn)處形成穩(wěn)定的兩列渦街.因此,MC 的漩渦脫落被顯著抑制,進(jìn)一步導(dǎo)致橫向動力響應(yīng)降低.在Ur=5 處,UC在P+S 尾渦模式中表現(xiàn)出較強的振動響應(yīng),渦結(jié)構(gòu)的垂直距離被拉大,從而對振動產(chǎn)生實質(zhì)性的抑制.Bao 等[29]也發(fā)現(xiàn)了類似的現(xiàn)象.在區(qū)間B 內(nèi),各工況下MC 的橫流向運動在Ur=7 附近發(fā)生頻率鎖定現(xiàn)象,此時y方向的振幅最大,如圖4 所示,且鎖定區(qū)間的范圍會隨著雷諾數(shù)的增大而擴(kuò)大.頻率比對y方向的振幅影響較小,各工況下的振幅比較接近.特別的,當(dāng)Re=160 時,在Ur=6 處的橫流向振幅,r=2.0工況Ymax/D=0.322 遠(yuǎn)小于1.073 (r=1.0) 和1.025(r=1.5).頻率比和雷諾數(shù)對x方向的振幅影響較大,Re<160 時,x方向的振幅隨頻率比的增大而增大,特別的是,當(dāng)Re=100,Ur=9 時,x方向的振幅與頻率比成反比關(guān)系.當(dāng)Re=160 時,r=1.0 工況下MC 在x方向上的振幅是其他頻率比工況下振幅的2~3 倍,但在Ur=7 時Xmax/D=0.073 與其他工況下的振幅接近.在區(qū)間C 內(nèi)(包括C1和C2),MC 逐漸退出鎖定區(qū)間,圓柱振動響應(yīng)逐漸減弱,因此MC在x和y方向的振幅逐漸減小并趨于穩(wěn)定.值得注意的是,當(dāng)r=1.0 時,在C1區(qū)間內(nèi)各雷諾數(shù)工況下MC 在x方向上的振動仍會保持相當(dāng)大的振幅,達(dá)到該工況下順流向振幅的最大值后逐漸下降趨于平穩(wěn).總的來說,MC 的動力響應(yīng)受UC 尾流振蕩的影響較大,相反,頻率比的影響較小.

    圖4 不同雷諾數(shù)與頻率比工況下,MC 在x 和y 兩個方向最大振幅隨約化速度的變化Fig.4 Variation of the maximum vibrating amplitudes of midstream cylinder with reduced velocity under different Reynolds numbers and frequency ratios

    圖5 Re=160,r=2.0 時,流場瞬時渦量圖Fig.5 Instantaneous vorticity diagram of flow field at Re=160 and r=2.0

    DC 振幅的變化趨勢與MC 類似,但DC 的“鎖定”區(qū)間的范圍比MC 更寬,且會隨著雷諾數(shù)的增大而明顯擴(kuò)大.在區(qū)間A 內(nèi)DC 在x和y兩個方向的振幅較小,Ur=4 時出現(xiàn)與MC 類似的“弱鎖定”現(xiàn)象,此時振幅出現(xiàn)小幅度的“上升-下降”過程,如圖6 所示.特別的是,當(dāng)r=2.0 工況DC 在區(qū)間A 內(nèi)的振幅波動很小.由于UC 和MC 尾流的共同“屏蔽”作用,當(dāng)Ur≤6 時,各工況下DC 的兩個方向的振幅均小于UC 和MC 的振幅;但當(dāng)Ur>6 時,DC 在兩個方向的振幅均大于UC 和MC 的振幅.由此可見,本文的臨界約化速度在6~7 之間,且受雷諾數(shù)與頻率比的影響較小[16,30].在區(qū)間B 內(nèi),DC 在y方向上突然劇烈振動,隨Ur增大,橫流向振幅取得最大值后下降,但仍保持較大的振幅值.但Re=160 時,在整個B 區(qū)間內(nèi)橫流向振幅均不斷增大.值得注意的是,與其他頻率比工況相比,r=2.0 時DC 進(jìn)入大振幅區(qū)間的約化速度Ur=7.在區(qū)間C1內(nèi),隨著Ur的增大,不同雷諾數(shù)和頻率比工況下DC 的橫流向振幅在取得最大值后逐漸減小.當(dāng)Re=160 時,r=1.5 與其他工況完全不同,此時橫流向振幅在整個區(qū)間C1內(nèi)均保持較大的數(shù)值.值得注意地是,當(dāng)r=2.0 時,橫流向振幅在Ur=14 處突然減小,Ymax/D=0.849(r=2.0)約為r=1.5 工況的53%(Ymax/D=1.583).在區(qū)間C 內(nèi)DC順流向振幅受頻率比的影響比橫流向振幅更敏感,具體地說:當(dāng)r=1.0 時,隨著Ur的增大,DC 的順流向運動在Ur=11 附近達(dá)到該工況下的最大振幅值后逐漸下降,此時DC 在C1區(qū)間內(nèi)的順流向振幅遠(yuǎn)大于其他頻率比工況下的振幅.當(dāng)r>1 時,不同雷諾數(shù)下DC 順流向運動在區(qū)間C1內(nèi)較弱,振幅較小且呈現(xiàn)逐漸下降的趨勢.特殊的是,當(dāng)Re=160,r=1.5時,DC 順流向振幅在Ur>12 后出現(xiàn)明顯的躍升,并達(dá)到該工況下順流向振幅最大值Xmax/D=1.12.C2區(qū)間內(nèi)DC 在兩個方向的振幅進(jìn)一步減小,逐漸趨于穩(wěn)定.值得注意地是,當(dāng)Re<160,r=1.0 時,DC 的橫流向振幅在C2區(qū)間內(nèi)出現(xiàn)增大的趨勢.同時,MC也具有相同的變化規(guī)律.

    圖6 不同雷諾數(shù)與頻率比工況下,DC 在x 和y 兩個方向最大振幅隨約化速度的變化Fig.6 Variation of the maximum vibrating amplitudes of downstream cylinder with reduced velocity under different Reynolds numbers and frequency ratios

    3.2 流體力系數(shù)分析

    Lx=5.5D工況下,MC 對UC 的影響極其微小,UC 的尾流能夠得到充分發(fā)展,導(dǎo)致其振動響應(yīng)類似于單圓柱工況[18].另一方面,受到UC 尾流和自身渦脫落的共同作用,使得MC 和DC 振動響應(yīng)發(fā)生顯著變化[15].

    由圖7(a) 可知,隨約化速度的增大,UC 在各工況下阻力系數(shù)的平均值(CD-mean)逐漸增大,頻率比對CD-mean 的影響很小,雷諾數(shù)越大CD-mean 越大.特別的,在r=2.0 工況下,當(dāng)Re=160,Ur=4 時CD-mean=2.405 遠(yuǎn)大于其他雷諾數(shù)工況,此時UC 順流向振幅Xmax/D=0.192 3 分別是0.015 8 (Re=80)和0.022 9(Re=100)的12 倍和8 倍.UC 的CD-mean在Ur=5 處達(dá)到最大值隨后逐漸降低并最終保持穩(wěn)定,CD-mean 隨頻率比的增大而增大,但與雷諾數(shù)成反比.區(qū)間A 內(nèi)MC 受UC 尾流的屏蔽作用,MC的CD-mean 幅值較小且小幅下降,此時頻率比小的工況下CD-mean 反而較大.隨著MC 進(jìn)入大振幅區(qū)間,CD-mean 顯著增大并在Ur=7 時達(dá)到最大值并逐漸減小趨于穩(wěn)定,此時CD-mean 受頻率比的影響很小.值得注意地是,MC 的CD-mean 穩(wěn)定時的幅值大約是UC 的一半.DC 的CD-mean 隨約化速度的變化十分復(fù)雜,同時雷諾數(shù)與頻率比對CD-mean 的影響也較為顯著.CD-mean 在區(qū)間A 和區(qū)間B 內(nèi)的波動較大,分別在Ur=4 和Ur=8 時達(dá)到兩個極值.在區(qū)間C 內(nèi),雷諾數(shù)越大DC 的CD-mean 越大,隨Ur增大,CD-mean 逐漸減小并趨于穩(wěn)定.

    圖7 不同雷諾數(shù)工況下,圓柱阻力系數(shù)平均值的變化Fig.7 Variation of mean value of drag coefficient under different Reynolds numbers

    UC 在各工況下升力系數(shù)的均方根值(CL-rms)隨約化速度的變化趨勢與CD-mean 變化類似,從圖8(a)可以看出,各工況下UC 的CL-rms 大致呈現(xiàn)先增大后降低直到平穩(wěn)的變化趨勢,并在Ur=4 處達(dá)到最大值.隨雷諾數(shù)的增大,UC 的CL-rms 變大,此時在y方向的振動響應(yīng)變化劇烈程度變強.但當(dāng)r≤1.5,Ur=5 時,雷諾數(shù)與UC 的CL-rms 卻成反比.MC 的CL-rms 在Ur=3 處就達(dá)到最大值,隨著Ur增大,CL-rms 迅速減小,頻率比和雷諾數(shù)對CL-rms的影響很小.在區(qū)間B 內(nèi),MC 的CL-rms 曲線變化的波動性明顯增強.當(dāng)Re=160 時,不同頻率比工況下CL-rms 曲線出現(xiàn)多個突變點.值得注意地是,當(dāng)Ur=6 時,MC 在r=2.0 工況下的CL-rms=0.081 遠(yuǎn)小于0.565(r=1.0)和0.577(r=1.5),導(dǎo)致此MC 的在y方向的振幅遠(yuǎn)小于其他頻率比工況下的振幅.在區(qū)間C 內(nèi),CL-rms 變化很小呈平穩(wěn)趨勢.“弱鎖定”作用使DC 的CL-rms 在Ur=4 處突然增大,隨后逐漸減小.當(dāng)5 <Ur<15 時,雷諾數(shù)和頻率比對DC 的CL-rms 的影響很強,DC 的CL-rms 曲線出現(xiàn)很強的波動.特別的,在Re=160 工況下,當(dāng)r=1.5,Ur=15時,MC 脫落的漩渦會與DC 上、下側(cè)的剪切層融合形成更強的漩渦,加劇DC 的振動響應(yīng),導(dǎo)致其升力系數(shù)的波動性會增強,如圖9 所示.i 時刻,DC 處于上方最大位移處,其尾流上方的漩渦基本已脫離,同時MC 尾流已傳播至其下方,使得DC 下表面受到較大吸力.在ii 時刻,DC 的尾流區(qū)會出現(xiàn)一個大漩渦,同時迎流面會受到MC 尾流的沖擊作用.隨后,DC 運動至下方最大位移處,其尾流上方的漩渦與MC 尾流已融合成一個整體,使得DC 上表面受到較大吸力.

    圖8 不同雷諾數(shù)工況下,圓柱升力系數(shù)均方根值的變化Fig.8 Variation of root mean square value of lift coefficient under different Reynolds numbers

    圖8 不同雷諾數(shù)工況下,圓柱升力系數(shù)均方根值的變化(續(xù))Fig.8 Variation of root mean square value of lift coefficient under different Reynolds numbers(continued)

    圖9 Re=160,r=1.5,Ur=15 時,DC 運動位移及升力時程曲線圖及對應(yīng)不同時刻渦量場云圖Fig.9 Time-history curve of displacement and lift coefficient of downstream cylinders and instantaneous pressure contours at different times when Re=160,r=1.5 and Ur=15

    3.3 流體力功率譜密度分析

    通過對圓柱體流體力系數(shù)功率譜密度(PSD)對比分析,發(fā)現(xiàn)圓柱流體力系數(shù)PSD 曲線在不同雷諾數(shù)工況下呈現(xiàn)不同的主峰幅值、頻譜成分及波動性,導(dǎo)致圓柱動力響應(yīng)特性發(fā)生較大變化.另外,研究發(fā)現(xiàn)r與Ur的變化對PSD 曲線有顯著影響,導(dǎo)致圓柱動力學(xué)響應(yīng)呈現(xiàn)明顯的差異性.本節(jié)針對r=1.0 和2.0 工況下三圓柱體流體力系數(shù)PSD 曲線的變化及其對圓柱的振動響應(yīng)的影響進(jìn)行詳細(xì)分析.

    在區(qū)間A 內(nèi),當(dāng)r=1.0 且Ur=3 時,各工況下UC 的流體力系數(shù)PSD 曲線出現(xiàn)兩階頻率峰值且曲線均比較光滑,隨著約化速度和雷諾數(shù)的增大,PSD曲線的波動性增強,如圖10(a) 所示,在一階頻率聚集了更多的能量,增強了UC 的振動響應(yīng).Ur增大到5 時,各雷諾數(shù)工況下升力系數(shù)PSD 曲線波動性減弱,曲線包含的頻率峰值增大至三階,有效增大了升力系數(shù)PSD 曲線包圍的面積,使得UC 獲得更多動能,在y方向上產(chǎn)生較大的振幅,具體分析見第3.4節(jié).在區(qū)間B 內(nèi),隨Ur增大,升力系數(shù)PSD 曲線波動性增強,但阻力系數(shù)PSD 曲線受約化速度的影響較小.各工況下升力系數(shù)PSD 曲線包含了主峰幅值與次峰幅值同等量級的多階頻率峰值,此時能量在次峰的分布比Ur=5 工況更加豐富,而主頻能量值更低,流體流向結(jié)構(gòu)的能量值減少,導(dǎo)致UC 的振動響應(yīng)比Ur=5 時有所減弱.隨Ur進(jìn)一步增大,升、阻力系數(shù)PSD 曲線趨于光滑,頻率峰值逐漸減少至二階,主峰幅值逐漸穩(wěn)定,次峰幅值逐漸減小直至可以忽略,UC 在區(qū)間C 內(nèi)的振動響應(yīng)逐漸減弱,在兩個方向的振幅趨于穩(wěn)定,最終PSD 曲線逐漸形成光滑的雙峰譜形態(tài).

    圖10 不同工況下,上游圓柱升力系數(shù)功率譜密度(紅線)與阻力系數(shù)功率譜密度(黑線)Fig.10 PSD of lift coefficient(red)and drag coefficient(black)of upstream cylinder under different case

    當(dāng)r增大至2.0 時,在區(qū)間A 內(nèi),如圖10(e) 和圖10(f)所示,各工況下流體力系數(shù)PSD 曲線的波動性很弱,隨著約化速度和雷諾數(shù)的增大,PSD 曲線的頻率峰值成分增大,UC 的運動變得復(fù)雜,其他特征與r=1.0 時類似.在區(qū)間B 內(nèi),各工況下圓柱流體力PSD 曲線均只出現(xiàn)三階頻率峰值,在Ur>7 后,PSD曲線的波動性逐漸減弱,主導(dǎo)頻率的能量值越來越低,UC 的動力學(xué)響應(yīng)也越來越弱,即圓柱振動響應(yīng)的強弱與主頻能量值成正比.需要說明的是,Ur=8時,雷諾數(shù)對升力系數(shù)PSD 曲線有顯著的影響.當(dāng)Re=160 時,升力系數(shù)PSD 曲線波動性與其他工況相比明顯增強,次峰及雜頻占據(jù)的能量值增多,此時各工況下升力系數(shù)PSD 曲線的主頻能量值幾乎相同.升力系數(shù)PSD 曲線的波動性對UC 運動的產(chǎn)生較大影響,Re=160 時UC 橫流向振幅Ymax/D=0.131 相較于Re=100 時的Ymax/D=0.366 下降了64%,同時比Ur=7 時的Ymax/D=0.705 下降了81%.研究發(fā)現(xiàn),當(dāng)r=1.0 時同樣具有類似規(guī)律.進(jìn)入?yún)^(qū)間C后,PSD 曲線隨約化速度的增大而變得光滑,雷諾數(shù)越小曲線的波動也越來越弱,最終形成光滑的兩階頻率峰值形態(tài).

    與UC 相比,MC 由于受到上游圓柱的尾流影響,流體力系數(shù)PSD 曲線隨雷諾數(shù)和約化速度的變化更加復(fù)雜.當(dāng)r=1.0 時,在區(qū)間A 內(nèi),不同工況下的升、阻力系數(shù)PSD 曲線均比較光滑,分別包含三階和兩階頻率峰值,但主頻能量值較小,圓柱振動響微弱.特別地,當(dāng)Ur=4 時,PSD 曲線的波動性隨雷諾數(shù)的增大而增強,如圖11(a)~圖11(c) 所示,PSD 曲線出現(xiàn)多階同等量級的頻率峰值,表明此時的泄渦頻率不位于鎖定區(qū)間.在區(qū)間B 內(nèi),當(dāng)Re<160,升力系數(shù)PSD 曲線的頻率峰值數(shù)量隨Ur的增大而增大,圓柱運動趨于復(fù)雜.Re=160 時,升力系數(shù)PSD 曲線波動性增強,并包含多階頻率峰值,主頻維持較高能量值,故此時MC 在y方向的運動能以較大的振幅持續(xù)到區(qū)間B 的末端.隨著Ur增大,阻力系數(shù)PSD 曲線變化較小,僅表現(xiàn)為波動性的差異,所以在區(qū)間B內(nèi)順流向振動相對平穩(wěn),振幅較區(qū)間A 有小幅度增大.特別地,當(dāng)Re=160 時,如圖11(c)所示,Ur=6 時阻力系數(shù)PSD 曲線主、次頻峰值的能量值接近,且能量在各雜頻的分布幅值也較多,圓柱順流向振動響應(yīng)較強.隨著Ur增大,阻力系數(shù)PSD 曲線的頻率幅值減少到三階(Ur=7),主頻能量值與Ur=6 工況相比差了一個數(shù)量級,順流向振幅Xmax/D=0.073 與Ur=6 工況比降低了三倍.流體力系數(shù)的PSD 曲線在區(qū)間C 后波動逐漸減弱,并逐漸形成光滑的三階頻率峰值曲線.

    圖11 不同工況下,中游圓柱升力系數(shù)功率譜密度(紅線)與阻力系數(shù)功率譜密度(黑線)Fig.11 PSD of lift coefficient(red)and drag coefficient(black)of midstream cylinder under different cases

    圖11 不同工況下,中游圓柱升力系數(shù)功率譜密度(紅線)與阻力系數(shù)功率譜密度(黑線)(續(xù))Fig.11 PSD of lift coefficient(red)and drag coefficient(black)of midstream cylinder under different cases(continued)

    當(dāng)r=2.0 時,在區(qū)間A 內(nèi)PSD 曲線的波動性隨雷諾數(shù)的增大而減弱,如圖11(f)所示,Ur=4 時MC的“弱鎖定”作用逐漸消失,流體經(jīng)過MC 時傳遞給圓柱的能量減少,導(dǎo)致圓柱的動力學(xué)響應(yīng)十分微弱.當(dāng)Ur=5 時,雷諾數(shù)對升、阻力系數(shù)PSD 的影響十分顯著.當(dāng)Re=100 時,阻力系數(shù)PSD 的主峰值頻率fscd=0.352 是升力系數(shù)PSD 的兩倍,從UC 脫落的漩渦直接從MC 上、下側(cè)通過,如圖12 所示,導(dǎo)致MC 的運動軌跡呈現(xiàn)為對稱的“8”字型,與Prasanth等[31]研究的結(jié)論一致.但當(dāng)Re=160 時,升、阻力系數(shù)PSD 曲線的各階峰值頻率完全重疊,UC 在一個周期內(nèi)脫落的孤立漩渦A 和一對反向的漩渦B+C,分別從MC 的上、下側(cè)向下游發(fā)展,由于上下側(cè)渦量不對稱,導(dǎo)致MC 的運動軌跡由“8”字形變成不對稱形狀.值得注意地是,上述工況下MC 受到UC 的“屏蔽”作用較大,且此時升、阻力曲線與位移的時程曲線的幾乎反向,導(dǎo)致MC 在Ur=5 時的振動響應(yīng)很弱.在區(qū)間B 內(nèi),MC 進(jìn)入鎖定區(qū)間,當(dāng)Ur接近8 時,PSD 曲線中會出現(xiàn)高階的頻率峰值,曲線波動性也大大減弱.此時MC 被完全鎖定,更多的能量由流體傳遞給圓柱,MC 獲得更多的動能而劇烈振動,x和y方向的振幅達(dá)到最大值.進(jìn)入?yún)^(qū)間C 后,Re<160 時升力系數(shù)的PSD 曲線的變化規(guī)律與r=1.0 時大體相同.但當(dāng)Re=160 時,在C1區(qū)間內(nèi),升力系數(shù)PSD曲線中出現(xiàn)了多階頻率峰值,包含的總能量值高,使得MC 產(chǎn)生較大的橫流向振幅.在C2區(qū)間中,升、阻力系數(shù)PSD 曲線波動性逐漸減弱,逐漸形成光滑的三階頻率峰值曲線.

    圖12 Re=100 與160 時中游圓柱流體力系數(shù)與位移時程曲線、運動軌跡及流場圖(r=2.0,Ur=5)Fig.12 Time history of the fluid force coefficient and displacement,the trajectory and the flow field of midstream cylinder at Re=100 and 160(r=2.0,Ur=5)

    DC 受到上中游兩圓柱的尾流共同作用,雷諾數(shù)、頻率比和約化速度對升、阻力系數(shù)PSD 的影響更加顯著.當(dāng)r=1.0 時,與MC 類似,PSD 曲線在Ur=4 時出現(xiàn)數(shù)量眾多的頻率峰值,主頻能量值較高,PSD 曲線的波動性隨雷諾數(shù)的增大而增強,如圖13(a)~圖13(c)所示,導(dǎo)致DC 在A 區(qū)間內(nèi)的運動在Ur=4 時最強烈.在區(qū)間B 內(nèi),隨Ur增大PSD 曲線中出現(xiàn)了更多數(shù)量的的頻率峰值,同時曲線的波動性也明顯增強,導(dǎo)致DC 的運動越趨復(fù)雜.特殊地,當(dāng)Re=100,Ur=7 時,升力系數(shù)PSD 曲線的頻率峰值比Ur=6 工況時更多,主頻能量值卻少了兩個數(shù)量級,致使DC 的橫流向振幅減小.當(dāng)Ur=8 時,雷諾數(shù)對PSD 曲線影響十分顯著,當(dāng)Re=100 時,如圖13(a)~圖13(c)所示,升力系數(shù)PSD 曲線光滑且中出現(xiàn)四階頻率峰值.當(dāng)Re=160 時,升力系數(shù)PSD 曲線由繁多的峰值構(gòu)成整體呈下降趨勢的曲線.研究發(fā)現(xiàn),上述的3 種不同的曲線形態(tài)是串列布置三圓柱體的PSD 曲線的最基本形態(tài).C 區(qū)間內(nèi),升力系數(shù)PSD曲線中出現(xiàn)多階頻率峰值,包含較高的能量值,導(dǎo)致下游圓柱產(chǎn)生強烈的動力學(xué)響應(yīng).當(dāng)Ur>13,PSD 曲線的波動性逐漸增強,雜頻占據(jù)的能量值增大,圓柱體獲取的能量值減少,導(dǎo)致下游圓柱的橫流向振幅逐漸減小.與UC 和MC 不同,當(dāng)Re=80 時,在Ur>15后,升力系數(shù)PSD 曲線中出現(xiàn)三階頻率峰值,主頻能量值隨Ur增大而變大,使得DC 橫流向振幅不斷增大.隨Ur增大,阻力系數(shù)PSD 曲線的主頻能量值逐漸降低,導(dǎo)致DC 順流向振幅減小.當(dāng)r=2.0 時,升力系數(shù)PSD 曲線隨約化速度的變化規(guī)律與r=1.0工況基本一致,如圖13(d)和圖13(f)所示,僅在頻率峰值數(shù)量及主頻能量值有細(xì)微差距.由此可見,流體力系數(shù)PSD 曲線隨約化速度的變化規(guī)律:Ur較小時PSD 曲線光滑,隨Ur增大PSD 曲線波動性增強,最終PSD 曲線逐漸光滑,與圓柱的運動軌跡的變化規(guī)律類似.

    圖13 不同工況下,下游圓柱升力系數(shù)功率譜密度(紅線)與阻力系數(shù)功率譜密度(黑線)Fig.13 PSD of lift coefficient(red)and drag coefficient(black)of downstream cylinder under different cases

    3.4 結(jié)構(gòu)頻譜機理分析

    結(jié)構(gòu)功率譜密度能反映單位頻帶內(nèi)激勵荷載功率隨頻率的分布,可以得到各個峰值頻率成分的幅值,從而衡量不同峰值頻率對結(jié)構(gòu)響應(yīng)的貢獻(xiàn),是結(jié)構(gòu)頻譜分析的重要手段,渦激振動中激勵荷載包括流體力、位移等[32].通過結(jié)構(gòu)功率譜密度分析,可以判斷激勵荷載信號響應(yīng)的強弱,從而對結(jié)構(gòu)動力學(xué)響應(yīng)進(jìn)行分析.圓柱的橫流向位移與升力系數(shù)功率譜密度曲線包圍的面積Ay和Acl分別代表位移與升力系數(shù)的平均功率值.本節(jié)以Re=100,r=1.0 工況為例,對結(jié)構(gòu)功率譜密度曲線和激勵荷載平均功率值與結(jié)構(gòu)振動的關(guān)系進(jìn)行分析.

    由圖14(a)可知,圓柱順流向位移平均功率值A(chǔ)y隨約化速度的變化趨勢與其在y方向的振幅變化趨勢類似.當(dāng)Ur=5 時,上游圓柱位移PSD 曲線出現(xiàn)fsy與3fsy二階頻率峰值,如圖15(a)所示,其主頻能量值較Ur=3 時更大,增大了位移平均功率值,激勵UC 發(fā)生大振幅振動.另一方面,此時渦脫落頻率fsy=0.18 與圓柱的固有頻率(fn,y=0.20)接近,UC 發(fā)生橫流向鎖定現(xiàn)象,產(chǎn)生較大振幅.當(dāng)Ur=5 和8 時,UC 的PSD 曲線表現(xiàn)為光滑的多峰譜特征,其運動軌跡表現(xiàn)為規(guī)則的“8”字形;但Ur=9 時,PSD 曲線的波動性增強并在主頻附近出現(xiàn)許多鋸齒形的窄帶峰值,削弱了主頻能量值,導(dǎo)致UC 運動的不規(guī)則性及隨機性增強,誘發(fā)產(chǎn)生小振幅順流向運動,如圖15(a)所示.進(jìn)一步研究發(fā)現(xiàn),在某約化速度區(qū)間內(nèi),圓柱橫流向振幅Ymax/D與位移信號平均功率值A(chǔ)y的大小成正比.

    圖14 Re=100,r=1.0 工況下圓柱位移與升力系數(shù)平均功率值隨約化速度的變化Fig.14 Variation of average power value of cylinder displacement and lift coefficient with reduced velocity at Re=100 and r=1.0

    通過對比發(fā)現(xiàn),在區(qū)間A 和B 內(nèi),升力系數(shù)平均功率值隨Ur的變化規(guī)律與升力系數(shù)均方根值的變化規(guī)律類似,且圓柱升力系數(shù)均方根值隨升力平均功率值A(chǔ)cl的增大而增大.與Ur=3 工況相比,如圖15(b)所示,當(dāng)Ur=5 時,UC 的升力系數(shù)PSD 曲線的頻譜成分更豐富,出現(xiàn)fscl,3fscl與5fscl三階頻率峰值,UC 的升力系數(shù)平均功率值增大,進(jìn)而增強了UC 的振動響應(yīng).Ur=9 工況下PSD 曲線的窄帶峰值數(shù)量明顯增多,且分布的頻帶更寬.在區(qū)間C 內(nèi),升力系數(shù)均方根值變化的波動加強,但維持了升力系數(shù)均方根值的變化趨勢.

    圖15 Re=100 與r=1.0 時,上游圓柱功率譜密度Fig.15 PSD of upstream cylinder at Re=100 and r=1.0

    需要注意的是,對不同約化速度區(qū)間的工況進(jìn)行頻譜分析時,必須優(yōu)先考慮振動頻率比(fs/fn,y)的影響.當(dāng)Ur=8 時,UC 升力系數(shù)平均功率值比Ur=3 工況更小,而此時UC 的Ymax/D達(dá)到0.271是Ur=3 工況下的近9 倍.這是由于Ur=8 時的主頻fscl=0.128 接近固有頻率fn,y=0.125,UC 發(fā)生橫流向共振并產(chǎn)生較大振幅.由此可見,對升力系數(shù)功率譜密度分析時,振動頻率比與激勵荷載平均功率值是影響結(jié)構(gòu)振動響應(yīng)強弱的兩個重要參數(shù).

    4 結(jié)論

    基于四步半隱式特征線分裂算子有限元方法,對串列布置三圓柱雙自由度結(jié)構(gòu)體系渦激振動問題進(jìn)行了數(shù)值模擬,分析了雷諾數(shù)、頻率比和約化速度的變化對結(jié)構(gòu)振動響應(yīng)及頻譜特性的影響,主要結(jié)論如下:

    (1)UC 兩個方向的振幅均比較小,雷諾數(shù)、頻率比對振幅的影響較小.MC 的鎖定區(qū)間的范圍隨著雷諾數(shù)的增大而擴(kuò)大,其動力響應(yīng)受頻率比的影響很小.本文的臨界約化速度在6.0 附近,且受雷諾數(shù)的影響較小.DC 的順流向振幅受頻率比和雷諾數(shù)的影響更為顯著.

    (2)UC 的CD-mean 與CL-rms 隨約化速度的變化出現(xiàn)先增大后降低直至平穩(wěn)的趨勢,受雷諾數(shù)和頻率比影響較小.由于MC 受到UC 尾流的影響,導(dǎo)致其CL-rms 與CD-mean 與UC 完全不同,約化速度較小時受雷諾數(shù)和頻率比的影響較大.DC 的CL-rms 受雷諾數(shù)和頻率比的影響十分顯著.中下游圓柱CD-mean的變化趨勢與UC 類似.

    (3) 圓柱流體力系數(shù)PSD 曲線的波動性隨雷諾數(shù)和頻率比的增大而增強,主峰能量值越大圓柱振動響應(yīng)越積極,PSD 曲線波動性越大,圓柱運動軌跡越不規(guī)則.當(dāng)升、阻力系數(shù)PSD 曲線的主峰重疊時,圓柱沿著非對稱軌跡運動.約化速度較小時PSD 曲線光滑,在大振幅區(qū)間內(nèi)PSD 曲線波動性增強且頻譜成分增多,最終PSD 曲線趨于光滑.

    (4)激勵荷載平均功率值隨Ur的變化趨勢與對應(yīng)的結(jié)構(gòu)動力響應(yīng)的變化類似.在同一約化速度區(qū)間內(nèi),結(jié)構(gòu)振動響應(yīng)的強弱與位移的平均功率值成正比.對升力系數(shù)功率譜密度分析時,區(qū)間A 和B 內(nèi)升力系數(shù)均方根值與升力平均功率值成正比.在其他區(qū)間振動頻率比對結(jié)構(gòu)振動響應(yīng)的影響更大.

    猜你喜歡
    雷諾數(shù)升力振幅
    高速列車車頂–升力翼組合體氣動特性
    無人機升力測試裝置設(shè)計及誤差因素分析
    基于自適應(yīng)偽譜法的升力式飛行器火星進(jìn)入段快速軌跡優(yōu)化
    基于Transition SST模型的高雷諾數(shù)圓柱繞流數(shù)值研究
    十大漲跌幅、換手、振幅、資金流向
    十大漲跌幅、換手、振幅、資金流向
    十大漲跌幅、換手、振幅、資金流向
    滬市十大振幅
    升力式再入飛行器體襟翼姿態(tài)控制方法
    失穩(wěn)初期的低雷諾數(shù)圓柱繞流POD-Galerkin 建模方法研究
    怎么达到女性高潮| 国产高清视频在线观看网站| 国产精品久久久av美女十八| 国产精品久久久久久亚洲av鲁大| 国产黄片美女视频| 精品久久久久久久久久久久久| 国产精品美女特级片免费视频播放器 | 岛国在线观看网站| 亚洲成人久久性| 日韩 欧美 亚洲 中文字幕| 亚洲精品美女久久av网站| 波多野结衣高清作品| av在线天堂中文字幕| 久久久久九九精品影院| 欧美丝袜亚洲另类 | x7x7x7水蜜桃| 国内精品久久久久精免费| 757午夜福利合集在线观看| 亚洲色图av天堂| 亚洲电影在线观看av| 国产精品乱码一区二三区的特点| 亚洲 欧美 日韩 在线 免费| 国产视频一区二区在线看| 欧美又色又爽又黄视频| 一进一出好大好爽视频| 精品久久久久久久毛片微露脸| 男人舔女人下体高潮全视频| 久久中文看片网| 淫秽高清视频在线观看| 中国美女看黄片| av片东京热男人的天堂| 黄色丝袜av网址大全| 午夜精品一区二区三区免费看| 久久久久免费精品人妻一区二区| 国产亚洲精品一区二区www| 不卡一级毛片| 久久婷婷人人爽人人干人人爱| 亚洲人成77777在线视频| 国产伦人伦偷精品视频| 一级片免费观看大全| 免费看十八禁软件| 欧美3d第一页| 成人18禁高潮啪啪吃奶动态图| 成人欧美大片| 成人特级黄色片久久久久久久| 成人av一区二区三区在线看| 国产在线观看jvid| 久久这里只有精品中国| 欧美黑人欧美精品刺激| 在线观看美女被高潮喷水网站 | 国产精品亚洲av一区麻豆| 欧美大码av| 免费无遮挡裸体视频| 黄色 视频免费看| 国产精品久久久人人做人人爽| 欧美性猛交╳xxx乱大交人| 一本久久中文字幕| 亚洲一卡2卡3卡4卡5卡精品中文| 国产精品久久久久久人妻精品电影| 国产区一区二久久| 看免费av毛片| 国产熟女xx| 精品一区二区三区视频在线观看免费| 中文字幕人成人乱码亚洲影| 国产欧美日韩精品亚洲av| 亚洲一码二码三码区别大吗| 一进一出抽搐动态| 黄色丝袜av网址大全| 欧美性长视频在线观看| 久久欧美精品欧美久久欧美| 国产高清激情床上av| 亚洲无线在线观看| 亚洲色图 男人天堂 中文字幕| 老司机午夜十八禁免费视频| 欧美最黄视频在线播放免费| 国产午夜精品论理片| 一二三四在线观看免费中文在| 精品少妇一区二区三区视频日本电影| 亚洲片人在线观看| 久久精品影院6| 中国美女看黄片| 国产97色在线日韩免费| 欧美国产日韩亚洲一区| 熟女少妇亚洲综合色aaa.| 久久午夜综合久久蜜桃| 欧美性长视频在线观看| 一本大道久久a久久精品| 天堂影院成人在线观看| 亚洲成人久久爱视频| cao死你这个sao货| 亚洲一码二码三码区别大吗| 久久香蕉激情| 亚洲男人的天堂狠狠| 日日夜夜操网爽| 女同久久另类99精品国产91| 97碰自拍视频| 国产91精品成人一区二区三区| 国内毛片毛片毛片毛片毛片| 十八禁人妻一区二区| 久久久国产成人精品二区| 国产探花在线观看一区二区| 免费无遮挡裸体视频| 最好的美女福利视频网| 男女床上黄色一级片免费看| 69av精品久久久久久| 最新在线观看一区二区三区| 又黄又粗又硬又大视频| 香蕉av资源在线| 日本一二三区视频观看| 日日摸夜夜添夜夜添小说| 国产精品久久久久久亚洲av鲁大| 男女视频在线观看网站免费 | 国产真人三级小视频在线观看| 久久亚洲真实| 久久人妻av系列| 欧美精品啪啪一区二区三区| 两个人的视频大全免费| 色在线成人网| 免费电影在线观看免费观看| 国产伦人伦偷精品视频| 国产成人aa在线观看| 黄色丝袜av网址大全| 两性午夜刺激爽爽歪歪视频在线观看 | 国产精品香港三级国产av潘金莲| 欧美日韩中文字幕国产精品一区二区三区| 每晚都被弄得嗷嗷叫到高潮| 婷婷精品国产亚洲av在线| 最近视频中文字幕2019在线8| 国产精品久久久久久人妻精品电影| www国产在线视频色| 99精品在免费线老司机午夜| 成人特级黄色片久久久久久久| 精品久久久久久,| 亚洲熟妇熟女久久| 俺也久久电影网| 亚洲国产精品久久男人天堂| 国产亚洲精品av在线| 夜夜看夜夜爽夜夜摸| 男人舔女人下体高潮全视频| av在线天堂中文字幕| 操出白浆在线播放| 夜夜躁狠狠躁天天躁| 欧美精品亚洲一区二区| 久久国产乱子伦精品免费另类| 国产一区二区激情短视频| 午夜福利视频1000在线观看| 日本三级黄在线观看| 不卡av一区二区三区| 国产av在哪里看| 女生性感内裤真人,穿戴方法视频| 国产成人影院久久av| 国产av一区在线观看免费| 少妇熟女aⅴ在线视频| 婷婷六月久久综合丁香| 日本一区二区免费在线视频| 国产免费男女视频| 国产高清视频在线观看网站| 国产v大片淫在线免费观看| 在线观看免费视频日本深夜| 欧美日韩亚洲综合一区二区三区_| 欧美另类亚洲清纯唯美| 久久精品人妻少妇| av超薄肉色丝袜交足视频| 99热这里只有是精品50| 天天躁狠狠躁夜夜躁狠狠躁| 成年人黄色毛片网站| 9191精品国产免费久久| 久久天堂一区二区三区四区| 欧美丝袜亚洲另类 | 999久久久国产精品视频| 精品免费久久久久久久清纯| 高潮久久久久久久久久久不卡| 国产精品久久久久久亚洲av鲁大| 日本 欧美在线| 怎么达到女性高潮| 午夜a级毛片| 好男人在线观看高清免费视频| 亚洲国产欧美网| 国产主播在线观看一区二区| 亚洲人成伊人成综合网2020| 亚洲精品在线美女| 我要搜黄色片| 精品国产乱子伦一区二区三区| 欧美日本亚洲视频在线播放| 99riav亚洲国产免费| 99国产精品99久久久久| 精品国产美女av久久久久小说| 欧美日韩一级在线毛片| 欧美精品亚洲一区二区| 九九热线精品视视频播放| 亚洲精品国产一区二区精华液| 亚洲熟女毛片儿| 欧美性猛交黑人性爽| 成熟少妇高潮喷水视频| 久久婷婷成人综合色麻豆| 国产一区二区三区在线臀色熟女| 男女那种视频在线观看| 97人妻精品一区二区三区麻豆| 日韩大码丰满熟妇| 欧美性长视频在线观看| 亚洲av日韩精品久久久久久密| 成人亚洲精品av一区二区| 三级国产精品欧美在线观看 | 免费在线观看成人毛片| 亚洲 欧美 日韩 在线 免费| 日韩av在线大香蕉| 一级毛片精品| 久久精品夜夜夜夜夜久久蜜豆 | 国产亚洲精品久久久久久毛片| av有码第一页| 99热这里只有是精品50| 久久精品aⅴ一区二区三区四区| 在线观看免费午夜福利视频| 国产精品久久久久久久电影 | 亚洲精品国产精品久久久不卡| 免费看日本二区| 人成视频在线观看免费观看| 亚洲午夜理论影院| 国产高清视频在线观看网站| 欧美另类亚洲清纯唯美| 天天躁夜夜躁狠狠躁躁| 色精品久久人妻99蜜桃| 不卡av一区二区三区| 哪里可以看免费的av片| 国产高清videossex| 国产三级中文精品| 欧美色视频一区免费| 757午夜福利合集在线观看| 18禁国产床啪视频网站| 亚洲男人的天堂狠狠| 久久精品人妻少妇| 国产精品 国内视频| 国产精品久久久久久久电影 | 欧美日韩乱码在线| 国产97色在线日韩免费| 日韩精品中文字幕看吧| 国产av一区二区精品久久| 嫁个100分男人电影在线观看| 啦啦啦韩国在线观看视频| 日韩欧美 国产精品| 午夜福利免费观看在线| av片东京热男人的天堂| 亚洲av五月六月丁香网| 久久精品国产亚洲av香蕉五月| 99riav亚洲国产免费| 看免费av毛片| 女人高潮潮喷娇喘18禁视频| 亚洲欧美一区二区三区黑人| 亚洲第一电影网av| 亚洲精品在线美女| 床上黄色一级片| 国产精品美女特级片免费视频播放器 | 成年免费大片在线观看| 成人国产综合亚洲| 国产欧美日韩一区二区三| 丝袜人妻中文字幕| 亚洲欧美日韩高清在线视频| 国产精品影院久久| 亚洲av第一区精品v没综合| 好男人在线观看高清免费视频| 国产高清有码在线观看视频 | 中文字幕高清在线视频| 国产三级中文精品| 午夜亚洲福利在线播放| 国产久久久一区二区三区| 欧美三级亚洲精品| 午夜福利免费观看在线| 熟女少妇亚洲综合色aaa.| 国产午夜福利久久久久久| 听说在线观看完整版免费高清| 亚洲国产精品久久男人天堂| 日本黄大片高清| 麻豆国产97在线/欧美 | 欧美 亚洲 国产 日韩一| 亚洲片人在线观看| 一本一本综合久久| 午夜日韩欧美国产| 欧美一区二区精品小视频在线| 国产精品av久久久久免费| 中出人妻视频一区二区| 动漫黄色视频在线观看| 亚洲狠狠婷婷综合久久图片| 色尼玛亚洲综合影院| 日韩免费av在线播放| 欧美黑人欧美精品刺激| 亚洲人成电影免费在线| 欧美绝顶高潮抽搐喷水| 97人妻精品一区二区三区麻豆| 男人的好看免费观看在线视频 | 亚洲人成77777在线视频| 国产精品亚洲av一区麻豆| 一二三四在线观看免费中文在| 国产一区二区在线av高清观看| 亚洲人成电影免费在线| 精品高清国产在线一区| 亚洲国产精品sss在线观看| 中国美女看黄片| 国产在线精品亚洲第一网站| 亚洲专区中文字幕在线| 精品免费久久久久久久清纯| 高清在线国产一区| 午夜视频精品福利| 熟妇人妻久久中文字幕3abv| 美女免费视频网站| www.www免费av| 男男h啪啪无遮挡| 国产精品av视频在线免费观看| 这个男人来自地球电影免费观看| 中文字幕久久专区| 日韩成人在线观看一区二区三区| 欧美黑人巨大hd| 在线播放国产精品三级| 12—13女人毛片做爰片一| 久久中文字幕人妻熟女| 国产av一区在线观看免费| 999久久久国产精品视频| 国产成人精品久久二区二区91| 人妻丰满熟妇av一区二区三区| av中文乱码字幕在线| 这个男人来自地球电影免费观看| 99精品在免费线老司机午夜| 美女大奶头视频| 91国产中文字幕| 欧美高清成人免费视频www| 久久精品夜夜夜夜夜久久蜜豆 | 国产成人啪精品午夜网站| 身体一侧抽搐| 午夜亚洲福利在线播放| 亚洲av成人精品一区久久| 欧美日本亚洲视频在线播放| 丝袜人妻中文字幕| 一区二区三区国产精品乱码| 国内揄拍国产精品人妻在线| 亚洲一卡2卡3卡4卡5卡精品中文| 欧美 亚洲 国产 日韩一| 美女高潮喷水抽搐中文字幕| 非洲黑人性xxxx精品又粗又长| 亚洲国产精品合色在线| 五月玫瑰六月丁香| 欧美中文日本在线观看视频| 亚洲成人免费电影在线观看| 最近视频中文字幕2019在线8| 伦理电影免费视频| 精品电影一区二区在线| 亚洲国产精品成人综合色| 国产激情偷乱视频一区二区| av欧美777| 一级作爱视频免费观看| 特级一级黄色大片| 老汉色av国产亚洲站长工具| 99久久综合精品五月天人人| 成人av一区二区三区在线看| 日本免费a在线| 夜夜躁狠狠躁天天躁| 午夜免费成人在线视频| 日韩大码丰满熟妇| 免费搜索国产男女视频| 精品国产乱子伦一区二区三区| 色在线成人网| 午夜精品在线福利| 两性午夜刺激爽爽歪歪视频在线观看 | 精品久久久久久久末码| 久久 成人 亚洲| 久久这里只有精品中国| 禁无遮挡网站| 一本大道久久a久久精品| 亚洲国产欧洲综合997久久,| 精品日产1卡2卡| 每晚都被弄得嗷嗷叫到高潮| 免费搜索国产男女视频| 一级毛片高清免费大全| 色综合亚洲欧美另类图片| 非洲黑人性xxxx精品又粗又长| 亚洲午夜精品一区,二区,三区| 国产午夜精品久久久久久| 午夜福利成人在线免费观看| 不卡一级毛片| 亚洲国产精品成人综合色| 亚洲国产欧美网| 在线免费观看的www视频| 欧美黄色淫秽网站| 三级国产精品欧美在线观看 | 十八禁人妻一区二区| 精品国产乱子伦一区二区三区| 久久精品人妻少妇| 91国产中文字幕| 老熟妇乱子伦视频在线观看| 亚洲中文av在线| 黄色丝袜av网址大全| 久久久久亚洲av毛片大全| 国模一区二区三区四区视频 | 国内精品一区二区在线观看| 亚洲男人天堂网一区| 最新在线观看一区二区三区| 妹子高潮喷水视频| 91av网站免费观看| 欧美一级a爱片免费观看看 | 国产97色在线日韩免费| 999久久久精品免费观看国产| 国产黄色小视频在线观看| 亚洲欧洲精品一区二区精品久久久| 午夜福利在线在线| 男人的好看免费观看在线视频 | 最近最新中文字幕大全免费视频| 国产在线精品亚洲第一网站| 亚洲自拍偷在线| 青草久久国产| 又黄又爽又免费观看的视频| 舔av片在线| 免费在线观看日本一区| 一进一出好大好爽视频| 91麻豆av在线| 日本五十路高清| 青草久久国产| 成人高潮视频无遮挡免费网站| 亚洲激情在线av| 在线看三级毛片| 亚洲精品中文字幕在线视频| 叶爱在线成人免费视频播放| 久久久久久久久免费视频了| 国产野战对白在线观看| 青草久久国产| 岛国在线观看网站| 国产精品爽爽va在线观看网站| 级片在线观看| 亚洲精品在线观看二区| 又大又爽又粗| 免费在线观看视频国产中文字幕亚洲| 亚洲av中文字字幕乱码综合| 午夜福利高清视频| 国产成人影院久久av| 人妻丰满熟妇av一区二区三区| 欧美在线黄色| 操出白浆在线播放| 一级黄色大片毛片| 一夜夜www| 看黄色毛片网站| 亚洲成人国产一区在线观看| 成年免费大片在线观看| 亚洲电影在线观看av| 搡老妇女老女人老熟妇| 免费一级毛片在线播放高清视频| 国产蜜桃级精品一区二区三区| 99久久综合精品五月天人人| 精品电影一区二区在线| 十八禁人妻一区二区| 久久久久国产精品人妻aⅴ院| 久久久久久久午夜电影| 久久人妻av系列| 色在线成人网| 午夜老司机福利片| 不卡一级毛片| 99久久无色码亚洲精品果冻| 欧美高清成人免费视频www| www.999成人在线观看| 一卡2卡三卡四卡精品乱码亚洲| 国产精品久久久久久亚洲av鲁大| 午夜精品在线福利| 免费无遮挡裸体视频| 91av网站免费观看| 天堂av国产一区二区熟女人妻 | 99国产精品一区二区蜜桃av| 欧美乱码精品一区二区三区| 视频区欧美日本亚洲| 脱女人内裤的视频| 黄色a级毛片大全视频| 色噜噜av男人的天堂激情| 首页视频小说图片口味搜索| 身体一侧抽搐| 男插女下体视频免费在线播放| 欧美不卡视频在线免费观看 | 欧美极品一区二区三区四区| 巨乳人妻的诱惑在线观看| 久久精品国产99精品国产亚洲性色| 我的老师免费观看完整版| 母亲3免费完整高清在线观看| 免费在线观看视频国产中文字幕亚洲| 少妇粗大呻吟视频| 亚洲18禁久久av| 国产高清videossex| www日本黄色视频网| 国产97色在线日韩免费| 男人舔女人的私密视频| 午夜福利视频1000在线观看| 在线a可以看的网站| 男女视频在线观看网站免费 | 一级毛片精品| 婷婷精品国产亚洲av在线| 亚洲人成电影免费在线| 国内精品一区二区在线观看| 91av网站免费观看| 久久人妻av系列| 精品久久久久久,| e午夜精品久久久久久久| 国内少妇人妻偷人精品xxx网站 | 男男h啪啪无遮挡| 亚洲av第一区精品v没综合| 日日摸夜夜添夜夜添小说| 一级毛片女人18水好多| 99久久99久久久精品蜜桃| 一本精品99久久精品77| 欧美日本亚洲视频在线播放| 琪琪午夜伦伦电影理论片6080| 美女 人体艺术 gogo| 久久伊人香网站| 国产黄片美女视频| 亚洲成人精品中文字幕电影| 婷婷精品国产亚洲av在线| 亚洲av成人不卡在线观看播放网| www.999成人在线观看| 狂野欧美白嫩少妇大欣赏| 久久国产精品影院| 深夜精品福利| 日本三级黄在线观看| 亚洲精品在线观看二区| 99riav亚洲国产免费| www日本在线高清视频| 五月玫瑰六月丁香| 九色国产91popny在线| 亚洲真实伦在线观看| 一级作爱视频免费观看| 色哟哟哟哟哟哟| 亚洲成av人片免费观看| 女同久久另类99精品国产91| 69av精品久久久久久| 无限看片的www在线观看| 国产黄色小视频在线观看| 黑人操中国人逼视频| 舔av片在线| 9191精品国产免费久久| 精品久久久久久久毛片微露脸| 精品高清国产在线一区| 一本综合久久免费| 久久香蕉国产精品| 亚洲精品av麻豆狂野| 一进一出好大好爽视频| 日本一本二区三区精品| 一区二区三区高清视频在线| 在线观看免费午夜福利视频| 国产三级中文精品| 国产精品美女特级片免费视频播放器 | 精品无人区乱码1区二区| 麻豆久久精品国产亚洲av| 国产真实乱freesex| 波多野结衣高清无吗| 亚洲成av人片免费观看| 一进一出抽搐gif免费好疼| 午夜精品在线福利| 亚洲成人久久性| 国产成人一区二区三区免费视频网站| 国产野战对白在线观看| 婷婷丁香在线五月| av欧美777| www.熟女人妻精品国产| 欧美黄色片欧美黄色片| 日日爽夜夜爽网站| 日韩欧美免费精品| 超碰成人久久| 麻豆av在线久日| 12—13女人毛片做爰片一| 国产成人一区二区三区免费视频网站| 在线视频色国产色| 免费在线观看成人毛片| 亚洲国产精品sss在线观看| 国产97色在线日韩免费| 搡老岳熟女国产| 国产区一区二久久| 亚洲乱码一区二区免费版| 国产单亲对白刺激| 在线观看一区二区三区| 午夜成年电影在线免费观看| 12—13女人毛片做爰片一| 久9热在线精品视频| 一区二区三区激情视频| 久99久视频精品免费| 精品人妻1区二区| 欧美丝袜亚洲另类 | 俺也久久电影网| 999久久久精品免费观看国产| 免费在线观看黄色视频的| 人人妻,人人澡人人爽秒播| 精品午夜福利视频在线观看一区| www日本黄色视频网| 国产成人啪精品午夜网站| 欧美黄色片欧美黄色片| 天堂av国产一区二区熟女人妻 | av片东京热男人的天堂| 欧美在线黄色| 法律面前人人平等表现在哪些方面| 伊人久久大香线蕉亚洲五| 18禁裸乳无遮挡免费网站照片| 亚洲成人免费电影在线观看| 国产精品爽爽va在线观看网站| 日日干狠狠操夜夜爽| 亚洲精品av麻豆狂野| 欧美一区二区国产精品久久精品 | 天天一区二区日本电影三级| 久久久久久久精品吃奶| 欧美不卡视频在线免费观看 | 亚洲国产日韩欧美精品在线观看 | 看黄色毛片网站| 色综合站精品国产| 成人精品一区二区免费| 嫩草影院精品99| 中文资源天堂在线| АⅤ资源中文在线天堂| 黄色毛片三级朝国网站| 叶爱在线成人免费视频播放| 日韩国内少妇激情av| 亚洲黑人精品在线| 日本黄色视频三级网站网址| 在线观看日韩欧美| 亚洲真实伦在线观看| 精品久久久久久成人av| 美女高潮喷水抽搐中文字幕|