• <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 06:26:20程友良焦慎俐
    力學(xué)學(xué)報 2021年7期
    關(guān)鍵詞:受迫振動方柱柱體

    程友良 焦慎俐 許 強 劉 鑫

    (華北電力大學(xué)(保定)動力工程系,河北保定 071000)

    引言

    柱體繞流這一現(xiàn)象普遍存在于自然界,最早注意到這一現(xiàn)象的是意大利著名畫家達·芬奇.法國科學(xué)家Benard 于1908 年進行了圓柱繞流實驗,最早在實驗中發(fā)現(xiàn)繞流圓柱后出現(xiàn)的渦旋脫落現(xiàn)象.1911 年馮·卡門教授對圓柱繞流后的渦旋脫落現(xiàn)象進行了理論推導(dǎo),確定了渦旋生成的條件[1].1925 年,Relf 等[2]通過在均勻水流中噴射油的方法,對Re=1100、自轉(zhuǎn)速比為0~ 10.4 時旋轉(zhuǎn)圓柱的流場進行攝影,得到了5 種轉(zhuǎn)速比流場的概貌.1959 年,Taneda[3]對圓柱繞流的尾流和渦街進行實驗研究,結(jié)果表明當Re>45 時尾流處出現(xiàn)卡門渦街,且渦在下游發(fā)展處形成了二次卡門渦街.基于實驗和數(shù)值模擬的方法,Tamura 等[4-7]通過對方柱繞流研究,發(fā)現(xiàn)相較于尖角方柱,圓角r/D=1/6 方柱的平均阻力和升力有一定的減小.Kawai[8]通過方棱柱和矩形棱柱的風(fēng)洞實驗,指出圓角、切角、凹角3 種處理方法中,圓角化處理對方棱柱的氣動彈性不穩(wěn)定性抑制效果最好,風(fēng)致振動的幅值隨著圓角參數(shù)的增大而減小.Hu 等[9]基于不同轉(zhuǎn)角半徑方柱的近尾跡實驗研究,指出當r/D從0 增大到1/2 時,脫落渦的最大強度逐漸減弱,與渦旋相關(guān)的環(huán)流逐漸減小50%,斯特勞哈爾數(shù)St增大約60%,渦旋對流速度隨尾跡寬度增加約25%,渦旋形成長度和尾跡閉合長度增加近一倍.Shadaram 等[10]通過風(fēng)洞試驗對矩形柱體的繞流進行研究,結(jié)果表明尾流處的湍流強度與方柱的長寬比有關(guān),但在高雷諾數(shù)下斯特勞哈爾數(shù)(St)基本不變.Kumar 等[11]基于實驗研究指出圓角參數(shù)對繞流體的流動特性有顯著影響,且r/D值越高,尾跡附近的脫落越均勻,并可以抑制圓柱體可能發(fā)生的失穩(wěn).工程上通過對方柱角部的處理(如圓角、切角、凹角等),達到降低流致振動,減小流動阻力和脈動升力的效果[12-16].Zhao 等[17]通過數(shù)值研究指出,在Re=200 的情況下,當圓角參數(shù)r/D=0.1 時,圓角完全抑制了馳振;當圓角參數(shù)r/D=0.01 時,對馳振的抑制作用不大;當圓角參數(shù)r/D為0.03 和0.05 時,圓角有效地抑制了馳振,馳振幅度極小;振動頻率是位移的主頻,渦脫落頻率是升力系數(shù)的主頻.Rocchio 等[18]基于大渦模擬指出圓角化柱體的流動特征與尖角柱體的流動特性有明顯的差異,即使r/D的值很小,平均回流區(qū)的長度也會顯著增加.文獻[19-24]基于柱體繞流,更好地分析高層建筑、架空線路在風(fēng)中的性能,改進結(jié)構(gòu),提高安全系數(shù).

    受迫振動作為一種簡化的運動形式相對容易實現(xiàn),并且流致振動和受迫振動的流場信息在一定程度上具有相似性,許多學(xué)者將研究重心放在受迫振動問題上.Bishop 等[25]較早進行了圓柱受迫振動的研究,對圓柱流向以及橫向的受迫振動進行了實驗.Bearman 等[26-27]對低雷諾數(shù)下橫流振動圓柱的流動進行研究指出最大振幅出現(xiàn)在鎖定區(qū)域的下限附近,而從靜止狀態(tài)發(fā)展起來的柱體振幅和已經(jīng)振動柱體的自由流速度之間的最終響應(yīng)沒有明顯的差別;當振幅大于0.6D時,觀察到不對稱的渦旋脫落模式.Sarpkaya[28]曾提出,渦激振動研究的目的是根據(jù)受迫振動的結(jié)果來預(yù)測自激振動的運動和受力情況,或是通過自激振動的研究結(jié)果來預(yù)測受迫振動的升力、阻力系數(shù)等情況.Yokoi 等[29]對三維圓柱受迫振動進行數(shù)值模擬和實驗,發(fā)現(xiàn)了沿展向方向渦旋脫落的不同時性,并指出三維效應(yīng)對渦脫模態(tài)的影響.平煥等[30]對低雷諾數(shù)下的直圓柱和波浪型圓柱橫向受迫振動進行了數(shù)值模擬研究,通過改變圓柱振動的頻率和振幅,來分析二者所受升力、阻力及鎖定區(qū)間,并且觀察到鎖定區(qū)間內(nèi)二者尾渦模式有所不同.Ma 等[31]基于受迫振動實驗,對方柱脈動橫向力特性與氣動彈性不穩(wěn)定現(xiàn)象之間的聯(lián)系進一步分析研究,對非定常力特性與氣動彈性振動幅值之間的關(guān)系提出了新的見解.劉俊等[32]通過實驗研究分析指出渦激振動遲滯現(xiàn)象通常伴隨振動的幅值階躍和頻率躍遷.Huang 等[33]指出無論運動方向如何,如果流體與振動圓柱之間的相對速度的最大值和最小值保持相同,則可以產(chǎn)生類似的尾跡結(jié)構(gòu),同時指出流動與結(jié)構(gòu)之間能量傳遞的變化征兆不一定是尾跡模式的改變,而是渦旋脫落時間的改變.Qu 等[34]提出非線性耦合尾跡振子模型,它是能夠同時再現(xiàn)自由振動和強迫振動實驗的模型,但它存在一個重要挑戰(zhàn)即對于頻率有一定的限制,這是Zhang 等[35]基于已建立的硬質(zhì)涂層圓柱殼非線性強迫共振特性的非線性解析公式,指出NiCoCrAlY+YSZ 硬質(zhì)涂層有利于抑制柱體的過度振動.Chen 等[36]通過受迫振動法,指出柔性振動管上升力主要取決于管道本身在提升方向上的振幅,而受管道振動的自由度、旋轉(zhuǎn)方向和主振動方向的影響較小.Wang 等[37]通過對均勻流和線性剪切流中柔性長圓柱體的單頻自由振動與規(guī)定正弦運動的剛性圓柱體強迫振動進行對比研究,指出兩種情況下的水動力系數(shù)有很強的相似性,強迫振動與自由振動的尾跡截面非常相似.Chen 等[38]基于小間隙直徑比近壁圓筒在Re=200 時的流致橫向振動進行的數(shù)值模擬結(jié)果,指出運動振幅一般隨折減速度的增大而先增大后減小,沒有明顯的滯回轉(zhuǎn)變;當小間隙直徑比減小時,最大振幅減小,出現(xiàn)最大振幅的折合速度增大,且從圓柱上部脫落的順時針渦流仍然很強,而從圓柱下部脫落的逆時針渦流逐漸減弱,渦街向上偏轉(zhuǎn)較小.本課題組曾在超臨界雷諾數(shù)條件下,通過使用大渦模擬方法,研究了波浪型圓柱和折線形圓柱等變截面立式圓柱的繞流問題以尋找降低渦激振動的最優(yōu)解.通過對繞流時的均方根升力系數(shù)對比,發(fā)現(xiàn)折線形圓柱的減振效率大于波浪型圓柱[39].

    最近,杜曉慶等[40]針對不同風(fēng)攻角分析了方柱圓角化處理后的流場作用機理,結(jié)果表明,方柱在繞流時穩(wěn)定性較差,圓角化有助于改善其穩(wěn)定性.但他們并未結(jié)合柱體的受迫振動進行研究分析,前人對于圓角化受迫振動柱體的深入研究也未曾有過.故本文將進一步研究圓角化對受迫振動方柱繞流特性的影響機理.為了分析圓角化對受迫振動方柱造成影響的流場作用機理,本文采用了Fluent 軟件用戶自定義函數(shù)(UDF)中的DEFINE_ CG_MOTION 宏進行編程,并采用了流場計算域區(qū)域劃分以及動網(wǎng)格技術(shù)中的動態(tài)層技術(shù),在Re=200 時,考慮方柱截面不同圓角化的影響,研究了均勻流作用下5 種圓角化r/D=1/2,1/4,1/5,1/8 和0 受迫振動方柱的繞流特性,以期為受迫振動方柱繞流穩(wěn)定性的提高提供參考依據(jù).

    1 數(shù)值模型

    1.1 基本控制方程

    對于不可壓縮的黏性牛頓流體,其控制方程為Navier?Stokes 方程.在二維直角坐標系下,連續(xù)性方程為

    動量方程為

    式中u和v分別為x和y方向的速度分量,p為壓力,ρ 為流體的密度,t為流體運動時間,v為流體的運動黏度系數(shù).

    1.2 動網(wǎng)格守恒方程

    對于邊界移動的任意控制體積上的一般標量的守恒方程可以寫為

    式中,V為大小和形狀都隨時間變化的控制體積,?V為控制體積的運動邊界,ug為運動網(wǎng)格的運動速度,u為流體速度矢量,Γ 為耗散系數(shù),S為源項,? 為單位質(zhì)量某種物理量.

    對式(4)中時間導(dǎo)數(shù)項用一階向后差分公式可得

    式中,上標n和n+1 分別代表當前和下一步的時間層,第n+1 個時間步的控制體積由下式得到

    為了滿足網(wǎng)格的守恒定律,控制體積隨時間的導(dǎo)數(shù)可由下式求得

    式中,nf為控制體積上面網(wǎng)格的數(shù)量,Aj為面j的面積矢量.

    式中,?Vj是控制體上面j在時間Δt內(nèi)掃出來的空間體積.

    1.3 無量綱參數(shù)

    無量綱參數(shù)的使用是為了把不同條件下得到的結(jié)果進行類比.受迫振動涉及到的無量綱參數(shù)有流動參數(shù)和流固耦合參數(shù),下面將分別對這些參數(shù)進行介紹.

    1.3.1 雷諾數(shù)

    雷諾數(shù)是描述流體流動最基本的參數(shù),雷諾數(shù)表示了流體微元慣性力與黏性力的比值

    式中,ρ 和U分別為來流流體的密度與速度,D為柱體的特征長度,μ 為流體的動力黏度系數(shù).

    雷諾數(shù)的大小不同,所表示的物理意義也不同.黏性力傾向于使流體中的擾動衰減,而慣性力傾向于使流體中的擾動增加,因此雷諾數(shù)的大小決定了流動是處于層流還是湍流,一般當雷諾數(shù)小于2300時,流動為穩(wěn)定層流;當雷諾數(shù)大于2300 時,流動為不穩(wěn)定湍流.在圓柱繞流問題中,雷諾數(shù)的大小還會影響到柱體尾流渦旋的脫落形態(tài).

    1.3.2 斯特勞哈爾數(shù)

    在流體流過柱體后會出現(xiàn)一系列交替的渦旋,斯特勞哈爾數(shù)是用于表征其渦旋脫落特性的重要相似準則數(shù),也即渦旋脫落的無量綱頻率,用于表征渦旋脫落的快慢,是描述非定常性的重要特征參數(shù),其表達式為

    式中,fs為柱體渦旋脫落頻率,D為柱體的特征長度,U為來流流體的速度.

    St數(shù)的大小與雷諾數(shù)、柱體形狀以及表面粗糙度等因素有關(guān),其中最主要的影響因素為雷諾數(shù).在較大的雷諾數(shù)范圍內(nèi)斯特勞哈爾數(shù)的值沒有太大的變化(在亞臨界雷諾數(shù)區(qū)域內(nèi),St數(shù)穩(wěn)定在0.2 左右)[41].本文Re數(shù)為200,St數(shù)也基本可視為0.2 不變.

    1.3.3 升力、阻力系數(shù)

    在流體流經(jīng)柱體過程中,可以將柱體受到的總作用力分解為與來流方向垂直的升力和與來流方向平行的阻力兩部分.升力的產(chǎn)生原因是柱體尾流渦旋的脫落,升力的值隨時間周期性變化,穩(wěn)定工況下時均值為0.繞流所產(chǎn)生的阻力包括壓差阻力和摩擦阻力,摩擦阻力的產(chǎn)生是由于流體的黏性.壓差阻力的產(chǎn)生是由于流體流經(jīng)柱體產(chǎn)生了邊界層分離,在柱體后方渦旋脫落形成負壓區(qū),而在柱體前方,流體由于受到阻礙速度減小壓力增加,柱體前后壓力一大一小產(chǎn)生壓差.雖然繞流升力與阻力產(chǎn)生的原因不同,但是二者均是由摩擦力與壓力構(gòu)成.柱體繞流升力、阻力積分表達式如下

    式中,FD和FL分別為阻力與升力,D為特征長度,p和 τ 分別為壓力和摩擦力.

    對阻力、升力進行無量綱化可得阻力系數(shù)和升力系數(shù)

    1.3.4 無量綱振幅

    無量綱振幅為柱體振幅與柱體特征長度之比,是計算受迫振動的重要無量綱參數(shù)

    式中,A為柱體振幅.

    1.3.5 無量綱頻率

    無量綱頻率為柱體受迫振動的驅(qū)動頻率與靜止柱體自然瀉渦頻率之比

    式中,fs為柱體渦旋脫落頻率,fe為柱體受迫振動的驅(qū)動頻率.

    無量綱頻率在受迫振動中是十分關(guān)鍵的無量綱數(shù),它的改變會影響到鎖定區(qū)間以及尾渦模態(tài).

    1.4 計算域及模型網(wǎng)格

    本文中雷諾數(shù)為200,計算域為40D×20D的矩形,方柱中心距離上游入口距離為10D,距離下游出口為30D,距離上下對稱邊界距離為10D,其中D為方柱棱邊邊長.流體從左向右流動,左側(cè)上游為速度入口,右側(cè)下游為壓力出口,上下兩側(cè)邊界為對稱邊界,方柱為無滑移壁面.流場劃分網(wǎng)格時采用結(jié)構(gòu)網(wǎng)格,由于方柱進行受迫振動需要采用到動網(wǎng)格,為避免計算域內(nèi)的所有網(wǎng)格都運動而影響計算效率,在距方柱中心左右3D處各設(shè)置一個滑移面來分割靜止網(wǎng)格與動網(wǎng)格,如圖1 所示.

    圖1 計算域示意圖Fig.1 Schematic diagram of the calculation domain

    網(wǎng)格的運動通過Fluent 中的UDF 實現(xiàn),利用剛體的運動宏(DEFINE _CG_MOTION)把方柱的速度傳遞給網(wǎng)格,網(wǎng)格進一步通過動態(tài)層法進行更新,在上下邊界處對網(wǎng)格進行分裂和合并.在方柱周圍的網(wǎng)格劃分較密,遠離方柱處網(wǎng)格劃分稀疏,取壁面法向第一層網(wǎng)格高度為0.004D,網(wǎng)格增長率為1.1.在對控制方程進行離散時,速度與壓力耦合采用SIMPLEC 算法,各物理量參數(shù)殘差控制在10?5.

    通過把方柱的尖角優(yōu)化為圓角從而對其進行被動控制,圓角特征參數(shù)為r/D,其中r為圓角半徑.本文選r/D=1/2,1/4,1/5,1/8 和0,5 種不同的圓角參數(shù),對其進行柱體繞流模擬,并分析不同圓角參數(shù)對柱體繞流的影響.隨后對圓角方柱進行受迫振動的數(shù)值模擬,并將模擬結(jié)果進行比較,觀察圓角化對其受迫振動的影響.

    圓角化時,計算區(qū)域與尖角方柱計算域一致.采用相同的棱邊節(jié)點數(shù)和柱體壁面法向第一層網(wǎng)格高度,以保證模擬結(jié)果不受網(wǎng)格劃分方式的影響.同樣的,為了方便動網(wǎng)格的生成與計算,將整個計算域通過“interface”邊界分為3 部分,中間部分為動網(wǎng)格和圓角化方柱所在的區(qū)域,此處網(wǎng)格尺寸較小,柱體壁面處網(wǎng)格較密.如圖2 所示為r/D=1/4 的網(wǎng)格劃分情況.

    圖2 圓角方柱網(wǎng)格劃分情況Fig.2 Grid division of rounded square column

    2 計算結(jié)果及分析

    本文對方柱橫向受迫振動進行研究,方柱橫向受迫振動的方程為

    式中y為方柱橫向運動的位移,對其求導(dǎo)可得到速度方程

    2.1 網(wǎng)格無關(guān)性驗證

    分別取單條棱邊上節(jié)點數(shù)為60,90,120,三種網(wǎng)格劃分方式觀察對結(jié)果的影響;另取節(jié)點數(shù)為90時,壁面法向第一層網(wǎng)格高為0.002D,0.01D進行無關(guān)性驗證,并與文獻[41-45]實驗和數(shù)值模擬結(jié)果做對比進行正確性驗證,如表1 所示.

    從表1 中可以看出,當單條棱邊節(jié)點數(shù)為90 時,其阻力系數(shù)和St與其他實驗和數(shù)值模擬結(jié)果吻合較好,且網(wǎng)格高度為0.002D和0.004D模擬所得結(jié)果相近.在確保正確性的前提下,考慮到計算時間,故選擇棱邊節(jié)點數(shù)為90,第一層網(wǎng)格高度為0.004D的網(wǎng)格劃分方式進行后續(xù)模擬計算.

    表1 靜止方柱擾流結(jié)果比較Table 1 Comparison of perturbation results of stationary square column

    2.2 模型驗證

    在驗證網(wǎng)格無關(guān)性與邊界條件設(shè)置正確性后,對圓柱r/D=1/2 受迫振動進行模擬以驗證模型及自編UDF 的準確性.圓柱受迫振動與自激振動均包含流向振動和橫向振動,而在實際渦激振動中流向振動振幅相對于橫向振動振幅普遍較小,鑒于大多數(shù)學(xué)者將研究重點放在橫向振動上,因此本文也將進行柱體橫向受迫振動的模擬研究.

    對圓柱受迫振動的研究主要是在不同的振幅和振動頻率下對其升力、阻力系數(shù),鎖定區(qū)間以及尾渦脫落模式進行分析.在流體力學(xué)問題的研究中,為使研究的結(jié)果更加直觀和準確,通常引入無量綱參數(shù)進行分析研究,本文引入無量綱振幅A*與無量綱頻率F,由前述提及的St公式以及St值可知fs=0.386.本文對無量綱振幅為0.2 下不同無量綱頻率的受迫振動進行模擬,圖3 為不同頻率下的升力、阻力系數(shù)歷時曲線.

    由圖3 可知,在不同的無量綱頻率下升力、阻力系數(shù)歷時曲線變化趨勢有所不同.在圖3(a)中可以看出升力、阻力系數(shù)存在明顯的“差拍振動”現(xiàn)象,數(shù)值上波動較大存在周期性.隨著無量綱頻率的增加振動進入鎖定狀態(tài),此時升力、阻力系數(shù)曲線變?yōu)橐?guī)律的正弦周期曲線,“拍”現(xiàn)象消失,升力、阻力系數(shù)的數(shù)值要小于靜止繞流時的數(shù)值.隨著無量綱頻率的繼續(xù)增大,受迫振動仍處于鎖定區(qū)間內(nèi),升力、阻力系數(shù)的大小隨著無量綱頻率的增大而增大.當無量綱頻率進一步增大超過鎖定區(qū)間上限的無量綱頻率時,升力、阻力系數(shù)曲線再一次出現(xiàn)“拍”現(xiàn)象,且由于無量綱頻率的增加升力、阻力系數(shù)曲線振動更加劇烈,升力系數(shù)幅值也明顯增大.將本文模擬得到的鎖定區(qū)間范圍與Bearman[26-27]研究給出同工況下鎖定區(qū)間范圍相比,可以得到基本一致的結(jié)果.F=0.75,1.15 處于鎖定區(qū)間外,F=0.85,1.05 處于鎖定區(qū)間內(nèi),同時證明了本文自編程序和模擬的正確性.文獻[42,45-46]給出關(guān)于柱體繞流三維尾跡的Re的范圍,方柱在Re≤ 300 左右,圓柱Re≈ 180,而本文選取Re=200,并根據(jù)以上驗證的準確性可知,文本選取二維是可行的.

    圖3 Re=200,A*=0.2 時受迫振動的升力、阻力系數(shù)歷時曲線Fig.3 Diachronic curve of lift and drag coefficient of forced vibration when Re=200 and A*=0.2

    2.3 圓角化對靜止方柱繞流的影響

    圖4 為靜止方柱在不同圓角參數(shù)下的阻力系數(shù)和升力系數(shù)的歷時曲線.從圖4 中可以看出,由于圓角的存在,方柱繞流時阻力系數(shù)會有明顯的減小,r/D=0,1/8,1/5,1/4 時隨著圓角半徑的增加,減阻效果越來越好,在選取圓角參數(shù)為r/D=1/4 時,減阻效果最佳;r/D=1/2 時較于圓角化方柱r/D=1/8,1/5,1/4,其阻力系數(shù)增加,減阻效果下降;同樣,當方柱圓角化r/D=0,1/8,1/5,1/4 時,其升力系數(shù)也有明顯的減小,但隨著圓角的增大其減小程度不如阻力系數(shù)明顯;r/D=1/2 時較于圓角方柱,其升力系數(shù)增加,減阻效果較差.Tamura 等[4]曾對圓角方柱和尖角方柱的氣動特性進行了繞流模擬對比,其結(jié)果顯示圓角方柱的氣動阻力較尖角方柱明顯減小,與本文模擬結(jié)果相符,從而再次驗證了本文模擬結(jié)果的正確性.

    圖4 不同圓角半徑下的繞流參數(shù)Fig.4 Flow parameters under different fillet radii

    圖5 為r/D=1/4 和r/D=0 時的繞流渦量圖,從圖5 中可以看出當方柱存在圓角時其脫落渦旋之間的間隔變小,也即渦旋脫落的頻率加快,但其尾部渦旋脫落寬度有所變窄,渦量在數(shù)值上也較小于尖角方柱,因此阻力減小.綜上所述,對方柱進行圓角化可以降低阻力而升力降低較小,故能提高其繞流穩(wěn)定性.

    圖5 方柱繞流渦量圖Fig.5 Vorticity of the flow around the square column

    2.4 圓角化對受迫振動方柱的影響機理分析

    通過前文的模擬,圓角參數(shù)為r/D=1/4 的方柱有著較小的升力、阻力系數(shù),故對方柱進行圓角化可有效減阻.本節(jié)將通過模擬不同圓角參數(shù)下的受迫振動方柱,分析圓角化對受迫振動方柱繞流的影響機理.

    2.4.1 圓角化能改變其鎖定區(qū)間范圍

    (1)圓角參數(shù)r/D=1/4

    在無量綱振幅A*為0.2 和0.4 的工況下對不同無量綱頻率的受迫振動進行模擬,觀察其升力、阻力系數(shù)歷時曲線及鎖定區(qū)間.

    圖6 A*=0.2 時圓角方柱不同無量綱頻率下的升力、阻力系數(shù)歷時曲線Fig.6 Duration curve of the lifting resistance coefficient under different dimensionless frequency of square columns with rounded corners when A*=0.2

    圖7 A*=0.4 時圓角方柱不同無量綱頻率下的升力、阻力系數(shù)歷時曲線Fig.7 Duration curve of the lifting resistance coefficient under different dimensionless frequency of square columns with rounded corners when A*=0.4

    由圖6 和圖7 可知,圓角參數(shù)為r/D=1/4,在A*=0.2 工況下,當F=0.85 時振動為差拍振動,F=0.9 時受迫振動就已經(jīng)進入鎖定區(qū)間;當F=1.15 時仍處于鎖定,F=1.2 時又進入了差拍振動狀態(tài).由此可以推斷出A*=0.2 下圓角方柱的受迫振動鎖定區(qū)間.同理,在A*=0.4 工況下,當F=0.75 時振動為差拍振動,F=0.8 時受迫振動就已進入鎖定區(qū)間;當F=1.2 時處于鎖定,F=1.25 時又進入了差拍振動狀態(tài).由以上模擬結(jié)果可知,隨著圓角方柱受迫振動振幅的增大,其鎖定區(qū)間范圍也隨之增加.圓角化方柱同時具有方柱和圓柱的某些性質(zhì),在繞流過程中流場的變化不及尖角圓柱劇烈,且在振動過程中也不存在由馳振向渦激振動轉(zhuǎn)變的情況,可提高其振動穩(wěn)定性[41-44].

    (2)圓角參數(shù)r/D=1/8

    同樣對無量綱振幅A*為0.2 和0.4 的工況下對不同無量綱頻率的受迫振動進行模擬.

    結(jié)合圖8 和圖9 可知,r/D=1/8 圓角方柱在A*=0.2 工況下當F=0.8 時振動為差拍振動,F=0.85 時受迫振動就已經(jīng)進入鎖定區(qū)間;當F=1.2 時處于鎖定狀態(tài),F=1.25 時又已進入差拍振動狀態(tài).在A*=0.4 工況下當F=0.7 時振動為差拍振動,F=0.75 時受迫振動已進入鎖定區(qū)間;當F=1.2 時處于鎖定狀態(tài),F=1.25 時又進入了差拍振動狀態(tài).

    由此可知,r/D=1/8 參數(shù)下的圓角方柱其鎖定區(qū)間范圍同r/D=1/4 類似,隨著振幅比的增加而增寬,一方面說明了圓角方柱流場性質(zhì)相似,另一方面相互驗證了模擬結(jié)果的準確性.圓角半徑越大,其鎖定區(qū)間范圍越小,振動穩(wěn)定性越好.

    (3)圓角參數(shù)r/D=0

    同樣對r/D=0 參數(shù)下的圓角方柱,無量綱振幅A*為0.2 和0.4 的工況下對不同無量綱頻率的受迫振動進行模擬.當A*=0.2 時,在F=0.8 時振動表現(xiàn)為“差拍振動”,當F=0.85 時就已經(jīng)進入了鎖定狀態(tài),相對于A*=0.4 時較為提前;在F=1.25 時,仍處于鎖定區(qū)間內(nèi),鎖定區(qū)間的下限較A*=0.4 工況后移.方柱在A*=0.2 下的鎖定區(qū)間長度較A*=0.4 時寬,且由圖10 和圖11 可知,各個無量綱頻率下的阻力系數(shù)均較A*=0.4 時小.在鎖定區(qū)間內(nèi),隨著無量綱頻率的增加阻力系數(shù)也會相應(yīng)增加.由前述運動方程可知,當A較小時,方柱運動速度也較小,流場在相同的時間內(nèi)變化不劇烈,方柱在流場內(nèi)受到的阻力也就隨之降低.隨著無量綱頻率的增加,方柱運動速度幅值增加,同時在流場中受到的阻力也增加.

    圖8 r/D=1/8,A*=0.2 時圓角方柱不同無量綱頻率下的升力、阻力系數(shù)歷時曲線Fig.8 Duration curve of lifting resistance coefficient under different dimensionless frequency of rounded square column with r/D=1/8 and A*=0.2

    圖9 r/D=1/8,A*=0.4 時圓角方柱不同無量綱頻率下的升力、阻力系數(shù)歷時曲線Fig.9 Duration curve of lifting resistance coefficient under different dimensionless frequency of rounded square column with r/D=1/8 and A*=0.4

    圖10 r/D=0,A*=0.2 時圓角方柱不同無量綱頻率下的升力、阻力系數(shù)歷時曲線Fig.10 Duration curve of lifting resistance coefficient under different dimensionless frequency of rounded square column with r/D=0 and A*=0.2

    圖11 r/D=0,A*=0.4 時圓角方柱不同無量綱頻率下的升力、阻力系數(shù)歷時曲線Fig.11 Duration curve of lifting resistance coefficient under different dimensionless frequency of rounded square column with A*=0.4,r/D=0

    當A*=0.4 時,在F=0.85 時升力、阻力系數(shù)均呈不規(guī)則的振動狀態(tài),阻力系數(shù)的“差拍振動”現(xiàn)象較升力系數(shù)明顯.這是由于方柱振動既受到了給定驅(qū)動頻率的作用,又受到了靜止方柱自然瀉渦頻率的影響.當無量綱頻率繼續(xù)增大到0.9 時,這種現(xiàn)象消失,升力、阻力系數(shù)變?yōu)檩^規(guī)則的等幅值振動狀態(tài),當無量綱頻率增加到1.25,升力、阻力系數(shù)曲線再次出現(xiàn)拍動現(xiàn)象,由此可以推斷出方柱在A*=0.4 下的鎖定區(qū)間.從圖11 中還可以看出,當方柱處在非鎖定區(qū)間時,其升力、阻力系數(shù)幅值較大,造成這一現(xiàn)象的原因是方柱前緣后緣都有尖銳的棱邊,其渦旋脫落點是固定的,會使得本就處于“差拍振動”的流場變得更加復(fù)雜.與Bearman[26-27]研究的同雷諾數(shù)同無量綱振幅下圓柱的鎖定區(qū)間相比,方柱的整體鎖定區(qū)間更偏向于無量綱頻率大于1 處.

    2.4.2 圓角化能改變升力、阻力大小

    由上所述,圓角參數(shù)為r/D=1/5 時,其基本性質(zhì)及無量綱參數(shù)變化趨勢會與前兩種圓角化r/D=1/4,1/8 方柱相似,本小節(jié)著重對比不同圓角參數(shù)下的升力、阻力系數(shù)均值及最值的大小.

    由圖11 可知,在各個對應(yīng)工況下,尖角方柱的升力系數(shù)和阻力系數(shù)均大于圓角方柱的升力、阻力系數(shù).且隨著方柱圓角半徑的增加,升力、阻力系數(shù)數(shù)值也越來越小,對方柱進行圓角化被動控制的效果也越加明顯.

    不同圓角參數(shù)下的阻力系數(shù)曲線趨勢較為一致,由圖11 可知尖角與圓角的鎖定區(qū)間差別較大,圓角方柱的鎖定區(qū)間隨著振幅比的增加而有所增加,與圓柱的變化趨勢較為一致.對比圖12(a)和圖12(c)不同圓角參數(shù)下的阻力系數(shù)曲線可知,在鎖定區(qū)間內(nèi)阻力系數(shù)的數(shù)值相差較小,振幅比的提高使得圓角減阻作用減弱.另可觀察到,尖角方柱的鎖定區(qū)間偏向于無量綱頻率大于1 處,而經(jīng)過圓角化后的方柱,其鎖定區(qū)間整體左移,且鎖定區(qū)間范圍有所減小.

    在升力系數(shù)曲線圖12(b)和圖12(d)中,可以明顯觀察到尖角方柱在鎖定區(qū)間內(nèi)其升力系數(shù)是逐漸減小的,而優(yōu)化過后的圓角方柱隨著圓角半徑的增大和振幅比的提高,其升力系數(shù)在鎖定區(qū)間內(nèi)逐漸增大.在低振幅比和低圓角半徑下,升力系數(shù)值雖呈下降趨勢,但下降斜率不及尖角方柱的大,升力系數(shù)的數(shù)值也較小.由于方柱自身形狀的原因,在其發(fā)生振動時會在低無量綱頻率下先發(fā)生馳振,馳振振動劇烈,在橫向振動時代表振動劇烈程度的升力系數(shù)自然也會較大,隨著無量綱頻率的增大,振動由馳振向渦激振動轉(zhuǎn)變,升力系數(shù)也隨之減小.等到振動完全轉(zhuǎn)變?yōu)闇u激振動時,升力系數(shù)達到最小值,在這之后隨著無量綱頻率的增加,升力系數(shù)又隨之變大.馳振通常發(fā)生于帶有尖銳棱邊的柱體,例如方柱和三角形柱等,同樣會造成結(jié)構(gòu)體的安全性問題.本文通過對方柱進行圓角化處理,可有效減弱馳振的發(fā)生.且通過圓角化的處理,使方柱的振動性質(zhì)更接近于圓柱和六邊形柱體,提高了振動的穩(wěn)定性和安全性[41].

    圖12 不同圓角參數(shù)下的升力、阻力系數(shù)圖Fig.12 Lifting resistance coefficients under different fillet parameters

    2.4.3 圓角化能改變尾渦脫落范圍

    由圖13 所示,當方柱靜止時尾流處出現(xiàn)卡門渦街,尾渦模態(tài)呈現(xiàn)為經(jīng)典的2S 模態(tài).當方柱處于受迫振動的鎖定區(qū)間內(nèi)時,尾渦模態(tài)與靜止方柱相似同為2S,渦與渦之間的間距受自身振動頻率的影響,頻率越大間距越小.當方柱受迫振動還未進入鎖定區(qū)間時,其尾流場較為復(fù)雜,渦旋之間的流向間距變大,尾渦模態(tài)表現(xiàn)為P+S 模態(tài).當方柱受迫振動遠離鎖定區(qū)間時,隨著自振頻率的增大渦旋之間的流向間距進一步減小,而橫向間距逐漸拉大,在下游較遠處形成兩列渦旋.而圓柱分別在相同工況下受迫振動時,尾渦脫落模式始終呈現(xiàn)為2S 模態(tài),其脫落渦旋的間隔與自振頻率有關(guān).由圖13 還可知,圓柱的渦旋脫落點并不固定,方柱的渦旋脫落點固定在后緣棱邊處,方柱渦旋脫落后,會在方柱后緣處形成較大的負壓區(qū),負壓區(qū)的出現(xiàn)會使得速度變化劇烈,這是導(dǎo)致方柱和圓柱受迫振動差別較大的主要原因之一.

    圖13 等值線渦量圖Fig.13 Contour line vorticity diagram

    圖14 為圓角化受迫振動方柱在A*=0.2 時的渦量圖.由圖14 可知,圓角方柱不論是在進入鎖定區(qū)間之前還是在離開鎖定區(qū)間后,其尾渦脫落模式均為2S,與同工況下受迫振動圓柱的尾渦脫落模式相似.尖角方柱未進入鎖定區(qū)間時,其脫落模式為P+S 且尾跡寬度較大;在高無量綱頻率下,二者渦旋脫落模式均為2S,但圓角方柱的尾跡寬度較尖角方柱窄且脫落渦旋之間的間隔較大,渦旋脫落頻率較慢,負壓區(qū)有所減小和后移,流場較為穩(wěn)定.由于圓角的存在,其渦旋脫落點下移至圓角后緣處,進而使得分離層更易相遇,減弱了流場的擾動,因此圓角越大減阻效果越好.相較之下,方柱圓角化后尾渦流場情況更接近圓柱.綜上所述,對方柱進行圓角化可改變方柱在振動時受到的外力,有效減弱振動時流場的復(fù)雜程度及振動劇烈程度,雖不及圓柱流場穩(wěn)定,但仍可大幅提高柱體振動安全性.

    3 結(jié)論

    對Ansys Fluent 進行二次開發(fā),并對流場計算域進行區(qū)域劃分以便利用動網(wǎng)格技術(shù)中動態(tài)層法實現(xiàn)柱體受迫振動,從而實現(xiàn)對受迫振動柱體繞流流場的流固耦合模擬.在Re=200 時,對均勻流作用下圓角化r/D=1/2,1/4,1/5,1/8 和0 受迫振動方柱的繞流進行數(shù)值模擬研究,分析了圓角化受迫振動方柱的升力、阻力系數(shù)、尾流渦量和鎖定區(qū)間的變化規(guī)律,澄清了圓角化對受迫振動方柱穩(wěn)定性的影響機理.得到如下主要結(jié)論:

    圖14 r/D=1/4 時受迫振動方柱尾流渦量圖Fig.14 Vorticity diagram of forced vibration square column wake when r/D=1/4

    (1)與尖角方柱相比,圓角化方柱的升力、阻力系數(shù)更小,且隨著圓角半徑的增加,升力、阻力系數(shù)會進一步減小,不同圓角參數(shù)下阻力系數(shù)曲線的趨勢較為一致;圓角化方柱振動時的尾渦流場與圓柱相似,低振幅比下圓角方柱的渦旋脫落模式均為2S 模態(tài),尾流渦旋形狀更加規(guī)則,渦旋脫落尾跡更窄,流場劇烈程度更小;圓角方柱鎖定區(qū)間范圍有所減小,且圓角半徑越大,鎖定區(qū)間越小,圓角化方柱鎖定區(qū)間更偏向關(guān)于F=1 對稱.

    (2)更進一步地,圓角化方柱隨著圓角半徑的增大和振幅比的提高,其升力系數(shù)在鎖定區(qū)間內(nèi)逐漸增大;圓角化方柱可有效減弱方柱振動時的馳振效應(yīng),使方柱的振動與圓柱的振動相似,大大提高柱體振動時的安全穩(wěn)定性.

    本文通過二維數(shù)值計算,在低雷諾數(shù)Re=200下,對不同圓角化方柱進行研究分析,為研究柱體二維向三維流動轉(zhuǎn)變對應(yīng)的臨界雷諾數(shù)提供一定參考.應(yīng)用本文方法可進行多體模型以及多自由度的模擬.對于工程實際如海底電纜等低流速下的安全穩(wěn)定性具有一定參考意義.由于實際工程中雷諾數(shù)往往較高,故可把雷諾數(shù)增大,采用合適的湍流模型,對不同截面柱體的受迫振動進一步研究分析.本文中數(shù)值計算的柱體受迫振動均是在垂直于來流的橫向進行的,而在實際工程領(lǐng)域,柱體不僅在橫向有運動,同時在來流方向也會有運動,因此進行多自由度的受迫振動研究更具有工程實際意義.

    猜你喜歡
    受迫振動方柱柱體
    上游切角倒角小間距比串列方柱大渦模擬研究
    振動與沖擊(2022年7期)2022-05-04 05:18:38
    串列多方柱氣動特性的試驗研究
    振動與沖擊(2021年9期)2021-05-17 05:32:12
    不同倒角半徑四柱體繞流數(shù)值模擬及水動力特性分析
    海洋工程(2021年1期)2021-02-02 02:48:12
    基于多介質(zhì)ALE算法的柱體高速垂直入水仿真
    利用自制教具探究受迫振動與共振
    談擬柱體的體積
    外注式單體液壓支柱頂蓋與活柱體連接結(jié)構(gòu)的改進
    2017年中考數(shù)學(xué)模擬試題(十)
    方柱繞流中聚乙烯熔體流變行為的數(shù)值模擬
    中國塑料(2015年10期)2015-10-14 01:13:20
    永磁行星齒輪傳動系統(tǒng)受迫振動分析
    亚洲精品aⅴ在线观看| 在线免费观看不下载黄p国产| a级毛色黄片| 午夜精品一区二区三区免费看| 人妻制服诱惑在线中文字幕| 亚洲精品,欧美精品| 大码成人一级视频| 99久国产av精品国产电影| 国产91av在线免费观看| 人人妻人人爽人人添夜夜欢视频 | 亚洲国产精品国产精品| 人人妻人人爽人人添夜夜欢视频 | 美女cb高潮喷水在线观看| 国产爽快片一区二区三区| 熟女av电影| 91久久精品电影网| 少妇的逼好多水| 六月丁香七月| 99久久人妻综合| 内射极品少妇av片p| 亚洲精品乱码久久久久久按摩| 少妇人妻久久综合中文| 亚洲国产精品成人综合色| 午夜激情福利司机影院| 国产色婷婷99| 纵有疾风起免费观看全集完整版| 国产成人精品婷婷| 五月天丁香电影| av在线app专区| 国产毛片a区久久久久| 国产黄频视频在线观看| 午夜免费鲁丝| 中文字幕亚洲精品专区| 女人久久www免费人成看片| 午夜福利在线在线| 少妇裸体淫交视频免费看高清| 欧美xxxx性猛交bbbb| 亚洲在线观看片| 夫妻性生交免费视频一级片| 99久久中文字幕三级久久日本| 一二三四中文在线观看免费高清| 欧美成人午夜免费资源| 大陆偷拍与自拍| 亚洲欧美成人综合另类久久久| 观看美女的网站| 青春草国产在线视频| 成人午夜精彩视频在线观看| 一区二区三区精品91| 97在线视频观看| 国内精品宾馆在线| 免费人成在线观看视频色| 夜夜看夜夜爽夜夜摸| 久久久成人免费电影| 69av精品久久久久久| 嘟嘟电影网在线观看| 色5月婷婷丁香| 男女边摸边吃奶| 日本猛色少妇xxxxx猛交久久| 春色校园在线视频观看| 精品一区二区三卡| 亚洲精品乱码久久久久久按摩| 在线观看三级黄色| 精品一区二区三区视频在线| 国产精品国产av在线观看| 欧美精品一区二区大全| 亚洲最大成人av| 欧美xxxx性猛交bbbb| 日韩一区二区视频免费看| 国产女主播在线喷水免费视频网站| 亚洲av免费在线观看| 久久久久九九精品影院| 男人添女人高潮全过程视频| 日韩av不卡免费在线播放| 又粗又硬又长又爽又黄的视频| 国产亚洲精品久久久com| 精品亚洲乱码少妇综合久久| 国产av码专区亚洲av| 草草在线视频免费看| 天天躁日日操中文字幕| 亚洲av成人精品一区久久| 你懂的网址亚洲精品在线观看| 少妇人妻久久综合中文| 插阴视频在线观看视频| 成人一区二区视频在线观看| 亚洲一区二区三区欧美精品 | 国产欧美亚洲国产| 26uuu在线亚洲综合色| 观看美女的网站| 免费看日本二区| 久久亚洲国产成人精品v| 91午夜精品亚洲一区二区三区| 精品久久国产蜜桃| 成人鲁丝片一二三区免费| 中国美白少妇内射xxxbb| 国产亚洲5aaaaa淫片| 在线观看国产h片| av女优亚洲男人天堂| 精品视频人人做人人爽| 91久久精品电影网| 在线观看美女被高潮喷水网站| 国产精品一二三区在线看| a级一级毛片免费在线观看| 身体一侧抽搐| 国产一区二区亚洲精品在线观看| 成人鲁丝片一二三区免费| 99re6热这里在线精品视频| 日日啪夜夜撸| 国产av国产精品国产| 中文字幕亚洲精品专区| 国产毛片a区久久久久| 国产精品av视频在线免费观看| 久久午夜福利片| 欧美极品一区二区三区四区| 亚洲成色77777| 国产精品久久久久久精品电影| 18禁裸乳无遮挡免费网站照片| 在线观看人妻少妇| 亚洲成色77777| 亚洲av.av天堂| 免费观看的影片在线观看| 中文在线观看免费www的网站| 欧美三级亚洲精品| 久久久精品94久久精品| 少妇 在线观看| 久久亚洲国产成人精品v| 亚洲熟女精品中文字幕| 国产在视频线精品| 日韩电影二区| 国产黄色免费在线视频| 国精品久久久久久国模美| 精品99又大又爽又粗少妇毛片| 亚洲成色77777| 国产成人免费观看mmmm| 婷婷色综合www| 欧美97在线视频| 交换朋友夫妻互换小说| av国产精品久久久久影院| 亚洲av在线观看美女高潮| 国产老妇女一区| 国产成年人精品一区二区| 成年人午夜在线观看视频| a级一级毛片免费在线观看| 久久99精品国语久久久| 免费看不卡的av| 亚洲精品日韩在线中文字幕| 亚洲,一卡二卡三卡| 国产久久久一区二区三区| 成人免费观看视频高清| 丰满人妻一区二区三区视频av| 国产精品嫩草影院av在线观看| 欧美一级a爱片免费观看看| 天天躁夜夜躁狠狠久久av| 欧美性感艳星| 亚洲不卡免费看| 国产精品久久久久久久电影| 深爱激情五月婷婷| 亚洲天堂国产精品一区在线| freevideosex欧美| 亚洲婷婷狠狠爱综合网| 国产午夜精品一二区理论片| av免费观看日本| 好男人视频免费观看在线| 日韩视频在线欧美| 3wmmmm亚洲av在线观看| 下体分泌物呈黄色| 久久女婷五月综合色啪小说 | 午夜精品国产一区二区电影 | 啦啦啦中文免费视频观看日本| 久久精品国产鲁丝片午夜精品| 国产熟女欧美一区二区| 18禁在线播放成人免费| 老女人水多毛片| 91午夜精品亚洲一区二区三区| 在线观看美女被高潮喷水网站| 九九爱精品视频在线观看| 欧美日韩一区二区视频在线观看视频在线 | 亚洲伊人久久精品综合| 亚洲成人久久爱视频| 麻豆成人av视频| 草草在线视频免费看| 久久久久网色| 黄色配什么色好看| 精品久久国产蜜桃| 国产在线一区二区三区精| 永久免费av网站大全| 青春草亚洲视频在线观看| 蜜桃久久精品国产亚洲av| 国产免费福利视频在线观看| 久久99热这里只频精品6学生| 亚洲精品一区蜜桃| 日日啪夜夜撸| 男男h啪啪无遮挡| 乱码一卡2卡4卡精品| 久久热精品热| 日本爱情动作片www.在线观看| 亚洲精品乱码久久久v下载方式| 天堂中文最新版在线下载 | 在线观看人妻少妇| 日本色播在线视频| 国产美女午夜福利| 在线a可以看的网站| 一边亲一边摸免费视频| 丝瓜视频免费看黄片| 成人毛片60女人毛片免费| 久久综合国产亚洲精品| 男女无遮挡免费网站观看| 免费高清在线观看视频在线观看| 午夜福利视频1000在线观看| 国产一区二区三区综合在线观看 | 最近中文字幕高清免费大全6| 精品人妻一区二区三区麻豆| 国产精品偷伦视频观看了| 97热精品久久久久久| 午夜视频国产福利| 亚洲av.av天堂| 少妇熟女欧美另类| 国产午夜福利久久久久久| a级毛片免费高清观看在线播放| 干丝袜人妻中文字幕| .国产精品久久| 国国产精品蜜臀av免费| av在线app专区| 亚洲精品日韩在线中文字幕| 大香蕉97超碰在线| 尾随美女入室| 亚洲精品乱久久久久久| 性色av一级| 最近中文字幕2019免费版| 91久久精品国产一区二区三区| 国产欧美日韩一区二区三区在线 | 色视频www国产| 成年版毛片免费区| 欧美性猛交╳xxx乱大交人| 97在线人人人人妻| 日本午夜av视频| 久久久色成人| 欧美日韩综合久久久久久| 国产高清国产精品国产三级 | 日本免费在线观看一区| 亚洲av一区综合| tube8黄色片| 国产 一区精品| 春色校园在线视频观看| 精品少妇黑人巨大在线播放| 日本猛色少妇xxxxx猛交久久| 国产爽快片一区二区三区| 中文字幕av成人在线电影| 卡戴珊不雅视频在线播放| 久热这里只有精品99| 日本色播在线视频| 免费在线观看成人毛片| 永久免费av网站大全| 免费观看无遮挡的男女| 麻豆成人av视频| 五月玫瑰六月丁香| 亚洲美女搞黄在线观看| 亚洲自拍偷在线| 国产精品熟女久久久久浪| 免费大片黄手机在线观看| 日本猛色少妇xxxxx猛交久久| 九九爱精品视频在线观看| 边亲边吃奶的免费视频| 日韩免费高清中文字幕av| 欧美老熟妇乱子伦牲交| 色哟哟·www| 国产在线一区二区三区精| av在线观看视频网站免费| 午夜福利视频1000在线观看| 日韩,欧美,国产一区二区三区| 国产一区亚洲一区在线观看| 99热这里只有精品一区| 国产一区二区三区av在线| 午夜福利视频1000在线观看| 国产男女超爽视频在线观看| 女人十人毛片免费观看3o分钟| 日韩欧美一区视频在线观看 | 精品人妻偷拍中文字幕| 毛片女人毛片| 高清午夜精品一区二区三区| 亚洲精品影视一区二区三区av| 国产成人一区二区在线| a级毛片免费高清观看在线播放| 在线观看人妻少妇| 亚洲综合色惰| 国产久久久一区二区三区| 国产精品成人在线| 久久精品夜色国产| 在线观看av片永久免费下载| 老司机影院成人| 亚洲av二区三区四区| 欧美一区二区亚洲| 亚洲精品一区蜜桃| 99热国产这里只有精品6| 少妇熟女欧美另类| 男女边摸边吃奶| 免费观看的影片在线观看| 国产精品.久久久| 免费大片18禁| 人人妻人人爽人人添夜夜欢视频 | 久久精品人妻少妇| 乱系列少妇在线播放| 精品国产一区二区三区久久久樱花 | 日韩制服骚丝袜av| 国产爽快片一区二区三区| 国产黄频视频在线观看| 国产精品国产三级国产专区5o| 欧美三级亚洲精品| 男插女下体视频免费在线播放| 三级经典国产精品| 一区二区三区四区激情视频| 久久精品国产亚洲av涩爱| 男女边吃奶边做爰视频| 日韩电影二区| 一边亲一边摸免费视频| 国产淫语在线视频| 国产黄a三级三级三级人| 男女那种视频在线观看| 99久久中文字幕三级久久日本| 国产爽快片一区二区三区| 国产真实伦视频高清在线观看| 国产精品99久久99久久久不卡 | 色婷婷久久久亚洲欧美| 大话2 男鬼变身卡| 在线 av 中文字幕| 亚洲精品aⅴ在线观看| 美女cb高潮喷水在线观看| 国产精品一区二区在线观看99| 日本爱情动作片www.在线观看| 欧美老熟妇乱子伦牲交| 欧美性感艳星| 特级一级黄色大片| 少妇裸体淫交视频免费看高清| 欧美高清成人免费视频www| 麻豆国产97在线/欧美| av在线老鸭窝| 精品午夜福利在线看| 国产大屁股一区二区在线视频| 一区二区av电影网| 伊人久久国产一区二区| 精品久久久久久电影网| 亚洲av在线观看美女高潮| 在线观看一区二区三区| 欧美高清性xxxxhd video| 精华霜和精华液先用哪个| 超碰97精品在线观看| 午夜日本视频在线| 人妻系列 视频| 国产永久视频网站| 国产成人aa在线观看| 两个人的视频大全免费| 乱系列少妇在线播放| 国产高潮美女av| 美女cb高潮喷水在线观看| 欧美区成人在线视频| 国产成人a∨麻豆精品| 不卡视频在线观看欧美| 久久久久久久精品精品| www.av在线官网国产| 校园人妻丝袜中文字幕| 九草在线视频观看| 欧美成人a在线观看| 欧美+日韩+精品| 在线观看国产h片| 97超碰精品成人国产| 国产又色又爽无遮挡免| 国产精品偷伦视频观看了| 中文乱码字字幕精品一区二区三区| 全区人妻精品视频| 男插女下体视频免费在线播放| 人妻制服诱惑在线中文字幕| 午夜福利视频精品| 天天一区二区日本电影三级| 免费观看无遮挡的男女| 欧美最新免费一区二区三区| 免费观看性生交大片5| 欧美+日韩+精品| 青青草视频在线视频观看| 久久久久久久久久久丰满| 又大又黄又爽视频免费| 一区二区三区四区激情视频| 成人高潮视频无遮挡免费网站| 免费av毛片视频| 蜜桃亚洲精品一区二区三区| 制服丝袜香蕉在线| 一区二区三区免费毛片| 欧美日本视频| 国产亚洲午夜精品一区二区久久 | 有码 亚洲区| 人妻 亚洲 视频| www.色视频.com| 国产大屁股一区二区在线视频| 色吧在线观看| 日本三级黄在线观看| 九九久久精品国产亚洲av麻豆| 国产久久久一区二区三区| 22中文网久久字幕| 国产 精品1| 91午夜精品亚洲一区二区三区| 午夜激情福利司机影院| 美女脱内裤让男人舔精品视频| 性色avwww在线观看| 精品国产一区二区三区久久久樱花 | 视频区图区小说| 国产精品国产三级国产专区5o| 国产黄频视频在线观看| 香蕉精品网在线| 亚洲激情五月婷婷啪啪| 久久久久性生活片| 久久久久久九九精品二区国产| 国产成人aa在线观看| 国产黄色免费在线视频| 国产片特级美女逼逼视频| 亚洲经典国产精华液单| 卡戴珊不雅视频在线播放| 久久女婷五月综合色啪小说 | 成人美女网站在线观看视频| 亚洲欧洲日产国产| 欧美潮喷喷水| 国产精品女同一区二区软件| 亚洲在久久综合| 国产亚洲午夜精品一区二区久久 | 一区二区三区四区激情视频| h日本视频在线播放| 国产精品精品国产色婷婷| 人人妻人人爽人人添夜夜欢视频 | 成人毛片60女人毛片免费| 欧美zozozo另类| 特级一级黄色大片| 亚洲高清免费不卡视频| 国产有黄有色有爽视频| 亚洲va在线va天堂va国产| 身体一侧抽搐| 免费黄色在线免费观看| 国产精品一区二区三区四区免费观看| 欧美日韩国产mv在线观看视频 | 国产午夜精品久久久久久一区二区三区| 色视频www国产| 黄片wwwwww| 久久人人爽人人片av| 日韩欧美精品v在线| 欧美日韩亚洲高清精品| 亚洲国产av新网站| 久久久久久久久久久丰满| 不卡视频在线观看欧美| 丰满乱子伦码专区| 人妻 亚洲 视频| 一级毛片久久久久久久久女| 亚洲av二区三区四区| 91在线精品国自产拍蜜月| 白带黄色成豆腐渣| 久久国内精品自在自线图片| 国产爱豆传媒在线观看| 99热全是精品| 成人国产麻豆网| 又黄又爽又刺激的免费视频.| 久久99热这里只有精品18| 日本一本二区三区精品| 国产高清国产精品国产三级 | 两个人的视频大全免费| a级毛色黄片| 3wmmmm亚洲av在线观看| 成人亚洲精品av一区二区| 国产精品一及| 少妇的逼水好多| 中文精品一卡2卡3卡4更新| 夜夜爽夜夜爽视频| 嫩草影院精品99| 看免费成人av毛片| 男男h啪啪无遮挡| 国产午夜福利久久久久久| 51国产日韩欧美| 天天躁夜夜躁狠狠久久av| 成人二区视频| 免费av毛片视频| 国产淫语在线视频| 成人特级av手机在线观看| 一个人看视频在线观看www免费| 精品久久久久久久末码| 在线a可以看的网站| 国产精品久久久久久精品电影小说 | 国内少妇人妻偷人精品xxx网站| 建设人人有责人人尽责人人享有的 | 国产伦在线观看视频一区| 狂野欧美激情性xxxx在线观看| 国产成人精品一,二区| 男女国产视频网站| 国产毛片a区久久久久| 成人漫画全彩无遮挡| 成年版毛片免费区| videossex国产| 伦理电影大哥的女人| 天美传媒精品一区二区| 全区人妻精品视频| 99热这里只有是精品在线观看| 蜜臀久久99精品久久宅男| 欧美丝袜亚洲另类| 在线观看一区二区三区| 联通29元200g的流量卡| 亚洲欧洲日产国产| 亚洲精品一二三| 男男h啪啪无遮挡| 两个人的视频大全免费| 亚洲精品影视一区二区三区av| 岛国毛片在线播放| 国产片特级美女逼逼视频| 国产一区二区三区av在线| 欧美成人精品欧美一级黄| 亚州av有码| av在线亚洲专区| 欧美zozozo另类| 成人免费观看视频高清| 黄色视频在线播放观看不卡| 欧美亚洲 丝袜 人妻 在线| 亚洲四区av| 久久99蜜桃精品久久| 国产日韩欧美亚洲二区| 久久ye,这里只有精品| 国产午夜精品一二区理论片| 大片免费播放器 马上看| 国产永久视频网站| 午夜精品国产一区二区电影 | 一级黄片播放器| 草草在线视频免费看| 中文字幕制服av| 精品午夜福利在线看| 久久久久久久久久成人| 大香蕉97超碰在线| 日韩av不卡免费在线播放| 内地一区二区视频在线| 久久精品综合一区二区三区| h日本视频在线播放| 久久久久九九精品影院| 久久99蜜桃精品久久| 最近最新中文字幕大全电影3| 欧美性猛交╳xxx乱大交人| 欧美日本视频| 亚洲欧美清纯卡通| 亚洲真实伦在线观看| 精品亚洲乱码少妇综合久久| 亚洲国产欧美在线一区| 久久鲁丝午夜福利片| 精品一区二区三区视频在线| 精品国产乱码久久久久久小说| 国产一级毛片在线| 2018国产大陆天天弄谢| 五月开心婷婷网| 亚洲欧美中文字幕日韩二区| 国产成人精品久久久久久| 天堂俺去俺来也www色官网| 日日撸夜夜添| 久久亚洲国产成人精品v| 一二三四中文在线观看免费高清| 精品人妻偷拍中文字幕| 久久久久久久久久久丰满| 噜噜噜噜噜久久久久久91| 国产精品99久久久久久久久| 久久久色成人| 久久久久精品久久久久真实原创| 午夜视频国产福利| 日韩成人av中文字幕在线观看| 天堂中文最新版在线下载 | 肉色欧美久久久久久久蜜桃 | 亚洲欧美一区二区三区黑人 | 亚洲欧美成人精品一区二区| 国产成人a∨麻豆精品| 亚洲av国产av综合av卡| 国产综合懂色| 免费看a级黄色片| 久久久久久久亚洲中文字幕| av国产精品久久久久影院| 欧美精品人与动牲交sv欧美| 亚洲精品第二区| 精品人妻视频免费看| 日韩电影二区| 欧美成人a在线观看| 亚洲av中文av极速乱| 91精品国产九色| 精品久久国产蜜桃| 国产高清不卡午夜福利| 99热这里只有是精品在线观看| 日韩制服骚丝袜av| 午夜精品国产一区二区电影 | 亚洲婷婷狠狠爱综合网| 亚洲最大成人中文| 亚洲最大成人av| 五月开心婷婷网| 亚洲最大成人手机在线| 亚洲av国产av综合av卡| 99热6这里只有精品| 黄片wwwwww| av在线老鸭窝| 日韩强制内射视频| a级毛色黄片| 日韩不卡一区二区三区视频在线| 国产高清三级在线| 亚洲最大成人av| 国产淫语在线视频| 色视频在线一区二区三区| 黄色日韩在线| 69av精品久久久久久| 天堂中文最新版在线下载 | 能在线免费看毛片的网站| 国内精品美女久久久久久| 性插视频无遮挡在线免费观看| 人妻夜夜爽99麻豆av| 欧美3d第一页| 亚洲精品aⅴ在线观看| 热re99久久精品国产66热6| 亚洲怡红院男人天堂| 国产探花在线观看一区二区| 国产国拍精品亚洲av在线观看| 国产精品不卡视频一区二区| 一二三四中文在线观看免费高清| 永久网站在线| 国产精品久久久久久精品电影小说 |