• <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 建模方法研究
    久久影院123| 亚洲怡红院男人天堂| av福利片在线观看| 美女视频免费永久观看网站| 免费av中文字幕在线| 亚洲av日韩在线播放| 国产色爽女视频免费观看| 欧美成人一区二区免费高清观看| 视频区图区小说| 精品国产一区二区三区久久久樱花 | av视频免费观看在线观看| 欧美精品亚洲一区二区| 国产亚洲91精品色在线| 午夜福利高清视频| 人妻一区二区av| 老熟女久久久| 人妻少妇偷人精品九色| a级毛片免费高清观看在线播放| 国产精品无大码| 成人毛片60女人毛片免费| 我的老师免费观看完整版| 精品亚洲乱码少妇综合久久| 精品久久久久久电影网| 久久影院123| 大陆偷拍与自拍| 亚洲国产毛片av蜜桃av| 亚洲国产精品成人久久小说| a级毛片免费高清观看在线播放| 汤姆久久久久久久影院中文字幕| 亚洲无线观看免费| 美女视频免费永久观看网站| 欧美zozozo另类| 99久久精品热视频| 制服丝袜香蕉在线| 蜜桃在线观看..| 亚洲av免费高清在线观看| 精品一区在线观看国产| 高清黄色对白视频在线免费看 | 男的添女的下面高潮视频| 三级经典国产精品| 熟女av电影| 少妇精品久久久久久久| 蜜桃久久精品国产亚洲av| 久久久久久久大尺度免费视频| 久久久久精品性色| 能在线免费看毛片的网站| 高清av免费在线| 久久99热这里只有精品18| 国产精品国产三级专区第一集| 免费av不卡在线播放| 少妇精品久久久久久久| 欧美精品亚洲一区二区| 一区二区三区免费毛片| 欧美少妇被猛烈插入视频| 我的女老师完整版在线观看| 午夜精品国产一区二区电影| 美女cb高潮喷水在线观看| 国精品久久久久久国模美| 国产成人a∨麻豆精品| 91在线精品国自产拍蜜月| 丰满乱子伦码专区| 国产精品国产三级国产av玫瑰| 成人黄色视频免费在线看| 国产女主播在线喷水免费视频网站| 色综合色国产| 国产毛片在线视频| 人妻 亚洲 视频| 在线观看一区二区三区激情| 亚洲丝袜综合中文字幕| 久久久久久久亚洲中文字幕| 久久久久久久久久久丰满| 在线亚洲精品国产二区图片欧美 | 高清午夜精品一区二区三区| 网址你懂的国产日韩在线| 欧美+日韩+精品| a级一级毛片免费在线观看| 亚洲第一av免费看| av线在线观看网站| 特大巨黑吊av在线直播| 少妇熟女欧美另类| 久久久色成人| av卡一久久| 大片免费播放器 马上看| 久久97久久精品| 亚洲精品国产成人久久av| 自拍欧美九色日韩亚洲蝌蚪91 | 国产乱人视频| 国产高潮美女av| 亚洲欧美日韩东京热| 在线播放无遮挡| 精品一品国产午夜福利视频| 国产一区有黄有色的免费视频| 亚洲精品456在线播放app| 妹子高潮喷水视频| 国产中年淑女户外野战色| 成年人午夜在线观看视频| 国产av国产精品国产| 午夜免费鲁丝| 偷拍熟女少妇极品色| av一本久久久久| 免费av中文字幕在线| 亚洲精品国产色婷婷电影| 亚洲色图综合在线观看| 七月丁香在线播放| 亚洲精品,欧美精品| 日本wwww免费看| 国产精品国产三级专区第一集| 午夜精品国产一区二区电影| 2021少妇久久久久久久久久久| 美女主播在线视频| 成人免费观看视频高清| 国产淫语在线视频| 国产成人aa在线观看| 国产精品免费大片| 国产 一区精品| 亚洲欧美日韩东京热| 亚洲国产毛片av蜜桃av| 特大巨黑吊av在线直播| 国产精品人妻久久久久久| 久久久久久久久久久免费av| 亚洲性久久影院| 肉色欧美久久久久久久蜜桃| 观看av在线不卡| 一级二级三级毛片免费看| 在现免费观看毛片| 一级黄片播放器| 又爽又黄a免费视频| 欧美极品一区二区三区四区| 下体分泌物呈黄色| 成人无遮挡网站| av又黄又爽大尺度在线免费看| 99久久精品热视频| 欧美日韩在线观看h| 国产av码专区亚洲av| 国产91av在线免费观看| 在线观看免费高清a一片| 五月玫瑰六月丁香| 国产亚洲5aaaaa淫片| 大香蕉久久网| 成人综合一区亚洲| 一级av片app| 欧美xxxx黑人xx丫x性爽| 国产精品人妻久久久影院| 婷婷色麻豆天堂久久| 插逼视频在线观看| 免费久久久久久久精品成人欧美视频 | 国产淫语在线视频| 中文字幕免费在线视频6| 久久午夜福利片| 精华霜和精华液先用哪个| 亚洲精品一二三| 97热精品久久久久久| av国产免费在线观看| 一级毛片aaaaaa免费看小| 亚洲成人手机| 亚洲精品久久久久久婷婷小说| 一个人免费看片子| 啦啦啦视频在线资源免费观看| 欧美精品人与动牲交sv欧美| 国产色婷婷99| 性色avwww在线观看| 久久久久精品久久久久真实原创| 在线观看免费视频网站a站| 亚洲成人手机| 女人十人毛片免费观看3o分钟| 偷拍熟女少妇极品色| 久久ye,这里只有精品| 国产高清有码在线观看视频| 久久6这里有精品| 国产成人免费无遮挡视频| 久久 成人 亚洲| 毛片一级片免费看久久久久| av在线播放精品| 亚洲,一卡二卡三卡| 亚洲欧美日韩卡通动漫| 亚洲国产最新在线播放| 久久国内精品自在自线图片| 夜夜骑夜夜射夜夜干| 我要看黄色一级片免费的| 最后的刺客免费高清国语| 毛片一级片免费看久久久久| 国产免费一区二区三区四区乱码| 久久综合国产亚洲精品| 久久国内精品自在自线图片| 下体分泌物呈黄色| 在线免费观看不下载黄p国产| 乱系列少妇在线播放| 黄色视频在线播放观看不卡| 少妇 在线观看| 亚洲av不卡在线观看| 美女视频免费永久观看网站| 秋霞伦理黄片| 欧美老熟妇乱子伦牲交| 国产一区二区三区综合在线观看 | 国产av精品麻豆| 波野结衣二区三区在线| 大香蕉久久网| 全区人妻精品视频| 久久久久久久久久人人人人人人| 国产精品久久久久久久电影| 日本-黄色视频高清免费观看| 日日啪夜夜撸| 91aial.com中文字幕在线观看| 日韩伦理黄色片| 国产亚洲欧美精品永久| 丰满乱子伦码专区| 人人妻人人添人人爽欧美一区卜 | 啦啦啦中文免费视频观看日本| 大片免费播放器 马上看| 日本一二三区视频观看| 亚洲va在线va天堂va国产| 中文字幕精品免费在线观看视频 | 成年美女黄网站色视频大全免费 | 免费高清在线观看视频在线观看| 免费大片黄手机在线观看| 网址你懂的国产日韩在线| 午夜福利高清视频| 日本av免费视频播放| 国产精品av视频在线免费观看| 又粗又硬又长又爽又黄的视频| 日本av手机在线免费观看| av卡一久久| 色视频www国产| 直男gayav资源| 久久国产精品男人的天堂亚洲 | 干丝袜人妻中文字幕| 中文字幕亚洲精品专区| 日本av手机在线免费观看| 激情五月婷婷亚洲| 18禁在线播放成人免费| 性色avwww在线观看| 乱系列少妇在线播放| 超碰97精品在线观看| 久久精品国产a三级三级三级| 久久国产精品男人的天堂亚洲 | 尤物成人国产欧美一区二区三区| 国产精品麻豆人妻色哟哟久久| 久久久精品免费免费高清| 午夜免费男女啪啪视频观看| 精品酒店卫生间| 边亲边吃奶的免费视频| 日本色播在线视频| 日韩中文字幕视频在线看片 | 一级爰片在线观看| 亚洲,欧美,日韩| 女性生殖器流出的白浆| 日韩av不卡免费在线播放| 午夜日本视频在线| 国产亚洲欧美精品永久| 国产精品熟女久久久久浪| 久久精品熟女亚洲av麻豆精品| 身体一侧抽搐| 高清日韩中文字幕在线| 天堂8中文在线网| 六月丁香七月| 久久久久精品性色| 秋霞在线观看毛片| 91精品国产九色| 中文字幕人妻熟人妻熟丝袜美| 久热这里只有精品99| 看免费成人av毛片| 亚洲精品自拍成人| 久久精品国产自在天天线| 少妇人妻久久综合中文| 国产成人精品婷婷| 久久 成人 亚洲| 日本欧美国产在线视频| 国产无遮挡羞羞视频在线观看| 亚洲美女搞黄在线观看| 男人添女人高潮全过程视频| 久久亚洲国产成人精品v| 少妇高潮的动态图| 亚洲av国产av综合av卡| 久久久久精品久久久久真实原创| 在线观看三级黄色| 深爱激情五月婷婷| 久久97久久精品| 国产av国产精品国产| 九色成人免费人妻av| 婷婷色综合大香蕉| 纵有疾风起免费观看全集完整版| 一二三四中文在线观看免费高清| 肉色欧美久久久久久久蜜桃| 看免费成人av毛片| 日韩电影二区| 亚洲不卡免费看| 国产亚洲一区二区精品| 在线亚洲精品国产二区图片欧美 | 亚洲第一区二区三区不卡| 免费人妻精品一区二区三区视频| 久久99热这里只有精品18| 夫妻午夜视频| 国产免费视频播放在线视频| 国产精品一区www在线观看| 极品教师在线视频| 国产色婷婷99| 日韩,欧美,国产一区二区三区| 黄片wwwwww| 色婷婷av一区二区三区视频| 国产精品秋霞免费鲁丝片| 亚洲av电影在线观看一区二区三区| av国产久精品久网站免费入址| 免费看av在线观看网站| 午夜免费观看性视频| 黄色一级大片看看| 国产欧美日韩精品一区二区| 国产又色又爽无遮挡免| 日本黄大片高清| 在线观看一区二区三区| 狂野欧美激情性xxxx在线观看| 男人和女人高潮做爰伦理| 国产精品伦人一区二区| 久久这里有精品视频免费| 街头女战士在线观看网站| 晚上一个人看的免费电影| 日日啪夜夜爽| 伦精品一区二区三区| 在线观看av片永久免费下载| 亚洲av男天堂| 一级毛片久久久久久久久女| 美女cb高潮喷水在线观看| 欧美高清性xxxxhd video| 黄色一级大片看看| 国产精品福利在线免费观看| 色哟哟·www| 久久午夜福利片| 丝瓜视频免费看黄片| 免费人成在线观看视频色| h日本视频在线播放| 国产一区二区三区av在线| 日韩强制内射视频| 在线 av 中文字幕| 国产日韩欧美在线精品| 欧美最新免费一区二区三区| 亚洲熟女精品中文字幕| 多毛熟女@视频| 国产亚洲精品久久久com| 最近最新中文字幕大全电影3| 女性生殖器流出的白浆| 26uuu在线亚洲综合色| 国产av国产精品国产| 亚洲精品成人av观看孕妇| 高清在线视频一区二区三区| 日日摸夜夜添夜夜添av毛片| 国产精品一区二区性色av| 久久久久精品久久久久真实原创| 舔av片在线| 国产淫片久久久久久久久| 日日啪夜夜爽| 国产在线视频一区二区| 久久久午夜欧美精品| 午夜激情久久久久久久| 一本久久精品| 老女人水多毛片| 波野结衣二区三区在线| 九九久久精品国产亚洲av麻豆| 亚洲精品国产成人久久av| 国产片特级美女逼逼视频| 中文字幕亚洲精品专区| 欧美成人午夜免费资源| 一个人免费看片子| 欧美激情极品国产一区二区三区 | 欧美日韩综合久久久久久| 少妇人妻一区二区三区视频| 精品久久久久久久久av| 国语对白做爰xxxⅹ性视频网站| 久热久热在线精品观看| 午夜免费男女啪啪视频观看| 亚洲第一区二区三区不卡| 国产伦精品一区二区三区视频9| 一级二级三级毛片免费看| 人妻少妇偷人精品九色| 免费黄频网站在线观看国产| 免费av不卡在线播放| 国产 一区精品| 久久青草综合色| 麻豆乱淫一区二区| 精品一品国产午夜福利视频| 全区人妻精品视频| 一级毛片 在线播放| 精品久久久精品久久久| 欧美成人精品欧美一级黄| 18+在线观看网站| 久久99精品国语久久久| 国产高清有码在线观看视频| 九九爱精品视频在线观看| 国产欧美亚洲国产| 婷婷色综合大香蕉| 久久国产亚洲av麻豆专区| 亚洲av不卡在线观看| 久久久精品免费免费高清| 边亲边吃奶的免费视频| 国国产精品蜜臀av免费| 国产视频内射| 看十八女毛片水多多多| 国产精品国产三级国产av玫瑰| 在现免费观看毛片| 尾随美女入室| 哪个播放器可以免费观看大片| 嫩草影院入口| 一级爰片在线观看| 这个男人来自地球电影免费观看 | 国产亚洲午夜精品一区二区久久| 能在线免费看毛片的网站| 免费看日本二区| 国产成人午夜福利电影在线观看| 亚洲欧美日韩无卡精品| av卡一久久| 亚洲av成人精品一区久久| 永久网站在线| 国产高清国产精品国产三级 | 最近中文字幕2019免费版| 男女下面进入的视频免费午夜| 在线观看免费高清a一片| 尤物成人国产欧美一区二区三区| 人人妻人人澡人人爽人人夜夜| 色婷婷av一区二区三区视频| 性色avwww在线观看| 美女中出高潮动态图| 91午夜精品亚洲一区二区三区| 国精品久久久久久国模美| 精品一区二区三卡| 亚洲精品日韩av片在线观看| 插阴视频在线观看视频| 国产一区有黄有色的免费视频| 久久精品人妻少妇| 久久人人爽av亚洲精品天堂 | 亚洲一区二区三区欧美精品| 午夜免费鲁丝| 97在线人人人人妻| 日韩人妻高清精品专区| 亚洲一级一片aⅴ在线观看| 人妻系列 视频| 国产精品人妻久久久久久| 亚洲怡红院男人天堂| 春色校园在线视频观看| 日韩电影二区| 蜜臀久久99精品久久宅男| 91久久精品国产一区二区成人| 亚洲图色成人| 亚洲激情五月婷婷啪啪| av卡一久久| 亚洲精品一二三| 久久久久久伊人网av| 国产中年淑女户外野战色| 亚洲av中文字字幕乱码综合| freevideosex欧美| 啦啦啦视频在线资源免费观看| 直男gayav资源| 成人毛片60女人毛片免费| 噜噜噜噜噜久久久久久91| 美女视频免费永久观看网站| 黄色日韩在线| 有码 亚洲区| 如何舔出高潮| 欧美区成人在线视频| 亚洲欧洲日产国产| 狂野欧美白嫩少妇大欣赏| 国产欧美亚洲国产| 多毛熟女@视频| 高清在线视频一区二区三区| 欧美精品一区二区大全| 日韩成人av中文字幕在线观看| 免费黄网站久久成人精品| 天堂中文最新版在线下载| 少妇的逼好多水| 午夜免费鲁丝| 日日啪夜夜爽| 身体一侧抽搐| 日韩强制内射视频| 国产精品99久久99久久久不卡 | 国产精品一二三区在线看| 夫妻性生交免费视频一级片| 在线观看国产h片| 国产无遮挡羞羞视频在线观看| 欧美激情极品国产一区二区三区 | 99热网站在线观看| 美女高潮的动态| 亚洲欧美一区二区三区黑人 | 中国美白少妇内射xxxbb| 国产成人91sexporn| 色视频www国产| 在线看a的网站| 国产av国产精品国产| 亚洲av中文字字幕乱码综合| 久久久久精品久久久久真实原创| 在线 av 中文字幕| 九草在线视频观看| 国产男人的电影天堂91| 妹子高潮喷水视频| 亚洲av欧美aⅴ国产| 在线观看美女被高潮喷水网站| h视频一区二区三区| 一区二区av电影网| 毛片女人毛片| 国产成人freesex在线| 国产亚洲欧美精品永久| 欧美高清性xxxxhd video| 男女边摸边吃奶| 亚洲综合精品二区| 岛国毛片在线播放| 欧美国产精品一级二级三级 | 国产片特级美女逼逼视频| 亚洲欧美日韩东京热| 高清av免费在线| 丝瓜视频免费看黄片| 男人舔奶头视频| 欧美日韩国产mv在线观看视频 | 大片电影免费在线观看免费| 亚洲av在线观看美女高潮| 欧美日韩视频精品一区| 国产精品99久久99久久久不卡 | 最后的刺客免费高清国语| 国产av一区二区精品久久 | 久久久久人妻精品一区果冻| 久久精品国产鲁丝片午夜精品| 国产一区二区三区av在线| 精品一区二区三区视频在线| 亚洲成人av在线免费| 国产av国产精品国产| 欧美老熟妇乱子伦牲交| 高清日韩中文字幕在线| 欧美成人精品欧美一级黄| 欧美精品人与动牲交sv欧美| 91精品伊人久久大香线蕉| 国产在线男女| 亚洲av男天堂| 久久久久久久国产电影| 在线观看免费日韩欧美大片 | 久久国产精品大桥未久av | 18禁裸乳无遮挡动漫免费视频| 极品教师在线视频| av网站免费在线观看视频| 王馨瑶露胸无遮挡在线观看| 伦精品一区二区三区| 欧美国产精品一级二级三级 | 亚洲美女搞黄在线观看| 亚洲欧美精品专区久久| 亚洲欧美中文字幕日韩二区| 黑丝袜美女国产一区| 精品少妇久久久久久888优播| 成年女人在线观看亚洲视频| 日本vs欧美在线观看视频 | 久久久久性生活片| 成人毛片a级毛片在线播放| av国产久精品久网站免费入址| 高清视频免费观看一区二区| 大码成人一级视频| 国产深夜福利视频在线观看| 在线观看免费日韩欧美大片 | 国产精品.久久久| 亚洲国产精品一区三区| 日韩电影二区| 高清日韩中文字幕在线| 亚洲av不卡在线观看| av在线观看视频网站免费| 99久国产av精品国产电影| 网址你懂的国产日韩在线| 青春草亚洲视频在线观看| 日韩电影二区| 国产精品国产三级国产av玫瑰| 纯流量卡能插随身wifi吗| 久久99精品国语久久久| 亚洲国产精品成人久久小说| 人体艺术视频欧美日本| 久久 成人 亚洲| kizo精华| a级毛片免费高清观看在线播放| 少妇精品久久久久久久| 亚洲欧美成人综合另类久久久| 亚洲av综合色区一区| 婷婷色麻豆天堂久久| 最近的中文字幕免费完整| 国产高潮美女av| 黄色视频在线播放观看不卡| 中文天堂在线官网| 国产淫片久久久久久久久| 久久久久久久久久人人人人人人| 久久ye,这里只有精品| 联通29元200g的流量卡| 久久99热这里只有精品18| 国产一级毛片在线| 成人黄色视频免费在线看| 妹子高潮喷水视频| 欧美日韩在线观看h| 婷婷色av中文字幕| 一区在线观看完整版| 久热久热在线精品观看| 国产免费福利视频在线观看| 你懂的网址亚洲精品在线观看| 免费久久久久久久精品成人欧美视频 | 精品久久久噜噜| 日韩大片免费观看网站| 纯流量卡能插随身wifi吗| 深夜a级毛片| 色视频www国产| 国产精品国产三级专区第一集| 国产精品久久久久久精品电影小说 | 日产精品乱码卡一卡2卡三| 麻豆精品久久久久久蜜桃| 国产午夜精品久久久久久一区二区三区| 国产精品一二三区在线看| 久久精品人妻少妇| 狂野欧美白嫩少妇大欣赏| 亚洲精品日韩在线中文字幕| 久久人人爽人人片av| 国产黄频视频在线观看| 人妻系列 视频| 中文欧美无线码| 热re99久久精品国产66热6| 亚洲天堂av无毛| 美女视频免费永久观看网站| 免费观看a级毛片全部| 国产精品蜜桃在线观看| 精品国产三级普通话版|