亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        旋轉振蕩板尾流的控制研究1)

        2021-11-09 06:26:14陳國孝邵傳平
        力學學報 2021年7期
        關鍵詞:旋渦尾流脈動

        陳國孝 劉 喆 邵傳平,2)

        *(中國計量大學流體檢測與仿真研究所,杭州 310018)

        ? (國家管網集團西氣東輸公司南京計量研究中心,南京 210046)

        引言

        橋梁顫振是結構從氣流中吸取的能量大于結構阻尼所消耗的能量時發(fā)生的扭轉(或扭、彎耦合)振動,其振動的振幅隨時間不斷加大,形成發(fā)散的自激振動,使結構遭受直接破壞.提出適應現代越來越輕型化的新型橋梁的顫振控制方法,為今后橋梁建設提供技術儲備很有必要.

        Larsen 和Larose[1]研究了Tacoma 大橋的淺H 型截面在扭轉振蕩中旋渦的形成及運動過程,以及對截面上下兩側壓力分布的影響,提出了旋渦脫落模式造成顫振的機制.當橋面扭轉振動時,會對其尾流產生影響,形成各種新的旋渦脫落模式,這些旋渦模式與振幅、振頻等參數密切相關.反之,新的旋渦脫落模式往往帶來氣動升力、阻力和力矩的幅值和頻率的變化,從而影響物體的振動.對于有限厚度板,顫振過程與旋渦脫落過程相互耦合,因此顫振的流動控制與旋渦脫落控制具有密切聯(lián)系.

        橋梁顫振受多種因素影響,包括長寬比、結構與空氣密度比、橋梁剛度,結構阻尼以及端部條件等,進行系統(tǒng)性控制研究極為困難.因此進行合理的簡化十分必要.風引起的橋梁、高層大樓的振動,海洋石油平臺減振板在海流中的振動,都可以用板或薄矩形柱作為簡化模型,進行實驗和數值研究.

        顫振可分為單模態(tài)振動和多模態(tài)振動兩種.單模態(tài)顫振較為簡單,即工程結構以同一個振頻振動,且振幅往往沿展向變化較緩慢.Bourguet 等[2]研究了多模態(tài)激發(fā)的情況,發(fā)現在鎖頻區(qū)域附近展向不同位置的振頻不同,但是在兩個節(jié)點之間,渦激振動也是局部單頻振動事件,流體以局部單振頻輸送能量給結構體.Law 等[3-4]研究了懸索橋顫振尾流,發(fā)現在扭轉第一模態(tài)和一些橫彎模態(tài)振動之下,出現旋渦脫落;并發(fā)現扭轉顫振至少在振幅最大的中間三分之一展長上是展向相關的.因此在展向相關尺度較大的情況下,應用二維模型,采用最大振幅及該處振頻,進行模擬研究是合理的.

        Carberry 等[5]觀察到強迫振動與自由振動的渦脫模式、應力和相位的變化都具有很強的相似性.Morse 和Williamson[6]利用強迫振動實驗得到的升阻力結果進行自由振動預測,并與Govardhan 和Williamson[7]自由振動實驗的響應和渦脫落模式比較,他們觀察到強迫振動和自由振動的結果之間具有一致性.由于強迫振動影響參數少且易于實驗中控制,因此在顫振導數等研究中常用其取代自由振動實驗.

        工程中橋跨結構繞流的雷諾數Re=106~ 107,大多數風洞實驗達不到這么高的要求.Bruno 和Fransos[8]研究了Re對橋梁顫振導數的影響,得到公式

        式中,Fi為某顫振導數,Fi0為雷諾數無關項,ai,bi為常數,不同的顫振導數其數值不同,一般ai為1 的量級,bi為更小的量級.折減風速Ur=V∞/(feH),其中V∞,fe,H,分別為來流速度、強迫振蕩振頻和板厚.可以看出,隨Re增大,各導數趨近于和雷諾數無關的常數.

        Matsuda 等[9]研究了Re對橋梁顫振導數與顫振臨界風速的影響.用Re=2×104~ 4×104,3×105~1.0×106和1.5×106等3 種不同雷諾數模型實驗得到的顫振導數計算的顫振臨界風速分別為70 m/s,70 m/s 和75 m/s,三者差別不大.對于強迫振蕩問題,振幅與振頻(折算風速)的影響需要首先考慮,Re的影響與其相比是屬于第二位的.

        長期以來,人們對橋梁消除顫振做了很多研究.顫振的控制可分為結構控制和流動控制兩種類型.結構控制通過增大結構的機械阻尼或增大結構剛度(以抵消氣動負剛度效應)達到減振的目的.流動控制是在不改變結構基本形狀的前提下,施加小的局部裝置對流動進行干預,從而改變結構的氣動力分布,達到抑制顫振的目的.板的流動控制方法有:(1)在板上開縫,調節(jié)兩側的空氣壓差,抑制顫振和馳振[10-13];在橋梁中間沿展向開縫隙,實驗證明在一定范圍內顫振臨界風速隨縫隙寬度的增大而增大.(2)減少展向相關性,抑制顫振[14-16];在尾緣沿展向每隔一定距離設置凸起薄片,控制參數包括薄片寬度、高度和間距.優(yōu)化的薄片布置增加背壓,抑制顫振.(3)豎直隔板、護板、導流片[17];在橋面中間沿展向設置中央隔板,在前緣和尾緣設置導流板.(4)風嘴(wind fairing)[18],將矩形截面橋梁的前后緣削尖,使其趨近于流線體.(5)噴射氣體干擾流場[19-22].(6) 表面振動[23-24];機翼負壓面璧面周期性變形.(7)前緣和尾緣襟翼、升力面[25-28];在橋跨前緣和尾緣分別設置繞軸旋轉振蕩的被動控制襟翼,用襟翼的拍動抑制橋跨結構不穩(wěn)定性.(8)可變形機翼[29].(9)反饋控制[30].綜上所述,已有的顫振控制(包括反饋和非反饋主動控制及被動控制)方法主要采用的控制裝置大都與主體結構連接在一起,從流動控制觀點看,在很多情況下,最佳控制效果往往發(fā)生在控制裝置與主體結構分離之時.且很多用于機翼顫振的主動控制方法(特別是反饋控制方法),難以用于尺度大得多的橋梁顫振控制.從實際應用考慮,應盡量采用被動離體的方法進行顫振控制研究.

        根據文獻[5-6]的研究可以推論,如果某種方法可以在較大參數范圍內抑制強迫振蕩板的尾流旋渦脫落,并減小脈動升力和力矩,則這種方法一定能夠抑制板的顫振.

        本文采用窄條控制件方法,對強迫旋轉振蕩板尾流進行抑制.在板的兩側對稱地放置兩個窄條控制件,根據工程實際,控制件分別放在板的前、后緣以及中間位置進行研究.板寬厚比采用研究最多的塔科馬大橋的比例B/H=5.根據Shimada 和Ishihara[31]的綜述研究,寬厚比B/H=2.8~ 6 情況下,流體繞過柔性板,所激發(fā)板的振動,主要是扭轉顫振.

        控制件的方法由Staykowski 和Sreenivasan[32]提出.他們在靜止的主圓柱下游放置一個尺度很小的圓柱,在主圓柱雷諾數Re<100 條件下,成功抑制了主尾流的旋渦脫落.此后文獻[33-36]在鈍體上游放置小圓柱,通過小圓柱尾流對鈍體邊界層或者分離剪切層的影響,達到抑制鈍體旋渦脫落的目的.但這些僅限于對靜止鈍體,且雷諾數較低情況下的控制研究.強迫振蕩柱體尾流控制方面的研究較少,文獻[37-40]采用單窄條控制件對流向振蕩圓柱尾流控制做了研究.曹夢圓等[41]還用單窄條對橫向振蕩圓柱尾流旋渦脫落的影響做了研究,取得一定控制效果.

        1 實驗與數值模擬方法簡述及結果對比

        1.1 實驗方法

        實驗在中國計量大學回流式風洞中進行,風洞實驗段長2.0 m,寬和高均為0.6 m,可提供速度范圍為0.6~ 30 m/s,湍流度小于0.5% 的均勻來流.矩形板由亞克力板制成,長0.55 m,寬B=0.15 m,厚H=30 mm,寬厚比B/H=5.板內沿長度方向嵌有一個直徑25 mm 的有機玻璃圓軸,軸線平行于板長方向,軸的一端延伸出板外.風洞的前后側壁由透明有機玻璃構成,后壁開有一直徑25 mm 的圓形槽口.位于板外的圓軸垂直地穿過槽口,并與風洞外的旋轉振蕩裝置連接.振蕩裝置由伺服電機,轉盤和連桿機構構成,通過調節(jié)連桿機構與轉盤的連接位置改變矩形板振蕩幅度,通過調節(jié)伺服電機轉速改變板的振蕩頻率.

        窄條控制件由鋁合金制成,其長度0.6 m,厚度3 mm,寬度b=10 mm,窄條寬度與板的厚度比b/H=0.33.兩個相同的窄條對稱地放置于矩形板上下兩側,窄條長度方向與板旋轉軸平行.窄條表面(長×寬) 與來流方向的夾角對控制效果具有影響.根據Shao 和Wei[42]的研究,控制效果與窄條迎風面積成正比.本文選擇具有最佳控制效果的情況,即窄條面與來流垂直放置.窄條流向位置如圖1(a)所示,分為3 種情況,分別為x/B=?0.5,0 及0.5.窄條縱向位置以y/H=±1.3 為主,測量y/H變化對控制效果的影響.

        圖1 模型與實驗裝置示意圖Fig.1 Sketch of test model and experimental equipment

        應用煙線技術進行流動顯示,觀察不同工況下振蕩板尾流旋渦脫落狀態(tài).如圖1(a)所示,在矩形板上游一定位置豎直地放置一根直徑40 μm 的鉬絲,鉬絲兩端通過導線與充電電容器兩極連接.在鉬絲上涂抹分析級甘油,電容器放電加熱鉬絲,使甘油揮發(fā)成煙霧.來流吹過鉬絲,煙霧在展向中間斷面形成二維煙流,煙流繞過板,在下游清晰地顯示出尾渦結構.尾渦結構的演化狀況由高速像機 Photron FASTCAM Mini UX50 以500 FPS 的速度拍攝記錄.

        在矩形板下游x/H=12.5 處沿Y軸在y/H=?5.67~ 5.67 的區(qū)間內均勻布置35 個測量點,用DANTECT Streamline 熱線風速儀熱線探頭,逐點測量尾流脈動速度.每個點的采樣時間為10 s,采樣頻率為512 Hz,脈動信號頻率分辨度為0.1 Hz.表1 為振幅5°、振頻1 Hz 且控制件在尾緣時,其中某一測點的平均速度U、脈動速度均方根u′及功率譜主峰峰值P隨測量時間T(數據點數N)的變化.可知10 s采樣時間足以得到精確的統(tǒng)計量.

        表1 熱線測量時間對測量結果的影響Table 1 Influence of hot wire measurement time on measurement results

        1.2 數值模擬方法

        振蕩板尾流旋渦的演化過程及升力、阻力和扭轉力矩的變化,采用標準k-ε模型用Fluent 軟件進行數值模擬.用Gambit 軟件進行用網格劃分,整體網格如圖2(a)所示,計算區(qū)域的上、下邊界到板振蕩中心的距離為2B,入口和出口邊界與板旋轉軸線的距離分別為3B和10B.振蕩板周圍設置非結構網格加密區(qū),用于動網格計算;該區(qū)上下邊界到軸心距離為1.2B,上、下游邊界與軸心距離分別為B和 6B.通過UDF 函數定義板的強迫旋轉振蕩運動.

        圖2 計算區(qū)域網格劃分Fig.2 Grid generation of computational domain

        網格質量在很大程度上影響著Fluent 數值模擬結果的可靠性,本文采用的是非結構化網格,選擇了6 種密度的網格,對雷諾數Re=2800 的均勻流繞靜止矩形板的流動進行網格無關性驗證.

        表2 展示了板在平衡位置靜止時,不同網格所算出的時均阻力系數CD和Strouhal 數St值.隨著網格數的增加,St和CD變化越來越緩慢,逐漸穩(wěn)定.本文選用mesh 5 進行振蕩板繞流模擬.

        表2 網格密度對Re=2800 靜止矩形板繞流數值計算結果的影響Table 2 Influence of grid density on numerical results of Re=2800 static rectangular plat

        2 控制結果與分析

        本文研究的雷諾數Re=V∞H/v=2800,板的振幅β=0~ 10°,振頻feH/V∞=0~ 0.0857.其中V∞=1.4 m/s 為來流速度,H=0.03 m 為板的厚度,ν為空氣運動黏度.

        2.1 數值模擬渦量場與流動顯示結果的比較

        通過煙線顯示和值數模擬的渦量場及壓力場3 種方法顯示振蕩板的尾流狀態(tài).圖3 和圖4 為不同振頻振幅下,未加控制和施加控制件于板前緣、中央和尾緣等不同位置時,實驗與數值模擬結果的對比情況.受實驗場地限制,相機鏡頭距離風洞較近,為了拍攝盡可能大的尾流區(qū)域,鏡頭焦點在板下游一定位置.由于板與鏡頭有夾角,圖片靠近左、右邊緣處有一定變形.但總體上看,渦量場模擬結果與煙線實驗結果相似度很高.

        圖3 是振幅5°,無量綱振頻0.043 時控制件在不同流向位置的流動對比圖.在煙線圖中,振蕩板未加控制時(圖3(a1)),尾流中存在明顯的大尺度的旋渦脫落,脫落模式為2S,即一個振蕩周期內有兩個單渦從板兩側脫落下來.施加窄條控制件后,如圖3(b1)、圖3(c1)和圖3(d1),尾流結構發(fā)生較大變化,小尺度旋渦從每個窄條兩側脫落并流向下游.板尾流中不再出現大尺度旋渦脫落.從模擬的渦量場看,旋渦結構變化更加明顯.如圖3(b2)、圖3(c2)和圖3(d2)所示,從窄條上剛脫落的小尺度旋渦清晰可見,但在向下游移動過程中,這些小尺度渦衰弱、變模糊.窄條對壓力場也有較大影響.無控制時,如圖3(a3),板逆時針旋轉到最大角度,板下側靠近前緣處有較大的高負壓區(qū),下游大尺度旋渦結構中心的低壓區(qū)和旋渦外的區(qū)域差別明顯.施加窄條后,如圖3(b3)、圖3(c3)和圖3(d3)所示,板下側的負壓區(qū)有所減小,下游也不再有明顯的低壓結構.

        圖3 振幅β=5°,振頻feH/V∞=0.043 時,控制件橫向位置固定于y/H=± 1.3,流向位置x/B 變化時的流動對比圖Fig.3 Comparisons between the plate wakes without elements and with the elements at y/H=± 1.3 and different x/B positions.Plate oscillation amplitude β=5°,frequency feH/V∞=0.0428

        圖3 振幅β=5°,振頻feH/V∞=0.043 時,控制件橫向位置固定于y/H=± 1.3,流向位置x/B 變化時的流動對比圖(續(xù))Fig.3 Comparisons between the plate wakes without elements and with the elements at y/H=± 1.3 and different x/B positions.Plate oscillation amplitude β=5°,frequency feH/V∞=0.0428 (continued)

        圖4 是振幅增大到10°,無量綱振頻為0.054 的情況.無控制時,數值模擬的渦量場與流動顯示高度吻合,如圖4(a1)和圖4(a2)所示,板后有2P 模式的旋渦脫落,即一個振蕩周期內,板上下兩側各有一對轉向相反的渦脫落下來.與2S 模式單渦相比,2P 模式的對渦尺度更大.施加控制件后,如圖4(b1)、圖4(b2),圖4(c1)、圖4(c2),圖4(d1)和圖4(d2)所示,小渦從每個控制件兩側脫落并流向下游,板的近尾流沒有大渦.再向下游,小尺度渦已衰弱、變模糊,但是,其有合并重組,形成新的大尺度渦的趨勢.

        從壓力場看,無控制情況如圖4(a3)所示,板順時針轉動到最大角度,板上側前緣附近和下側尾緣附近都出現較大的負壓區(qū),形成順時針扭矩,對抗板向平衡位置恢復.在板近尾流中,存在明顯的大尺度高壓和低壓結構.低壓結構由大尺度旋渦形成.高壓結構是由于大尺度對渦結構移動較慢,外流受其阻擋產生.施加控制件后,如圖4(b3)、圖4(c3) 和圖4(d3)所示,板下側尾緣的負壓區(qū)近于消失;尾流中的低壓和高壓結構或者強度變弱,或者尺度變小.

        圖4 振幅β=10°和振頻feH/V∞=0.054 時,控制件橫向位置固定于y/H=± 1.3,流向位置x/B 變化時的流動對比圖Fig.4 Comparisons between the plate wakes without elements and with the elements at y/H=± 1.3 and different x/B positions.Plate oscillation amplitude β=5°,frequency feH/V∞=0.0428

        2.2 模擬與測量的脈動速度功率譜比較

        為了定量比較控制前后旋渦脫落強度變化,將分布于x/H=12.5 橫線上的每個監(jiān)測點得到的脈動速度進行頻譜分析,再對35 個監(jiān)測點的功率譜進行平均,得到平均功率譜.平均功率譜上的主峰峰值,代表旋渦脫落強度.需要說明的是,風洞來流中存在湍流度,其湍流脈動的能量廣泛分布于各個頻率上,形成背景寬譜,疊加在熱線測量的功率譜上.而數值模擬中來流沒有湍流度,不存在背景寬譜,因此兩者功率譜有一定差別.

        圖5 為振幅β=2.5°,振頻feH/V∞=0.043 時,未加控制和施加控制后,數值模擬和熱線測量的脈動速度功率譜比較.不加控制時的功率譜如圖5(a)所示,模擬和測量的功率譜上都存在一個主峰和若干個次峰,其中主峰頻率與板振頻相同,次峰頻率依次為板振頻的2,3,4 等整數倍.主峰是由于板的振動形成鎖頻旋渦脫落,在尾流脈動功率譜上的反映.由于流動的非線性性質,板的振蕩在尾流中產生高次諧波,分別形成2 次、3 次、4 次等諧波脈動,在功率譜上出現各個次峰.數值模擬的主峰和3 階諧波峰值略高于測量值,而2 階和4 階諧波峰值略低于測量值.由于數值模擬是理想系統(tǒng),能量高度集中在主頻和各階諧波頻率上,而風洞存在來流湍流度及壁面摩擦等因素,形成額外的系統(tǒng)耗散,可能引起能量在主頻、諧波頻率和其他頻率上的再分配.

        圖5 振幅β=2.5° 和振頻fe H/V∞=0.043 時,無控制件與控制件位于橫線y/H=±1.3 上不同x/B 位置時的尾流脈動速度功率譜比較Fig.5 Comparisons between power spectra of fluctuating velocities in the wakes without and with the control elements at y/H=±1.3 and different x/B positions,the plate oscillation amplitude β=2.5° and frequency fe H/V∞=0.043

        在橫線y/D= ±1.3 上,分別在x/B=?0.5,0,0.5 等位置施加控制件后,如圖5(b)~ 圖5(d),主峰峰值降低一半以上,各諧波峰值也大幅降低,這說明控制件對鎖頻旋渦脫落及尾流其他規(guī)則脈動具有抑制作用.模擬得到的控制與無控制功率譜主峰比分別為0.5,0.5 和 0.46,而測量得到的譜峰比分別為0.47,0.46 和0.43,兩者差別不大,說明模擬方法得到的控制效果基本反映實際情況.由于監(jiān)測位置距離板較遠,窄條上脫落的小渦衰減嚴重,加之相鄰監(jiān)測點間距與控制件尾流寬度相近,因此控制件旋渦脫落未在功率譜上形成尖峰.

        如圖6(a)所示,隨著振幅增大到5°時,未加控制尾流的功率譜主峰值大幅增大.各次諧波峰值也都有所增大,但沒有主峰增幅大.施加控制件后如圖6(b),主峰值大幅減小,各次諧波峰值也都減小.實驗測量的各個峰值均低于數值模擬值.控制效果方面,即控制與未控制時功率譜主峰比,實驗值為0.45,數值模擬為0.48,也相差不多.振幅增大到10°時,未控制尾流脈動速度功率譜主峰值相比5° 時并未增大,如圖7(a)所示.而測量結果證明,隨著振頻增大,主峰值是迅速減小趨勢(未給出圖),這可能是因為測量點離板有一定距離,高頻振蕩信號衰減比低頻快的緣故.施加控制后,功率譜主峰和各諧波峰值都有減小.控制與無控制主峰比,數值模擬和實驗測量均為0.78.

        圖6 振幅β=5° 和振頻fe H/V∞=0.043 時,無控制件與控制件位于橫線y/H=±1.3 上、x/B = 0 位置時的尾流脈動速度功率譜比較Fig.6 Comparisons between power spectra of fluctuating velocities in the wakes without and with the control elements at y/H=±1.3 and x/B = 0,the plate oscillation amplitude β=5° and frequency fe H/V∞=0.043

        圖7 振幅β=10° 和振頻fe H/V∞=0.054 時,無控制件與控制件位于橫線y/H=±1.3 上、x/B = 0 位置時的尾流脈動速度功率譜比較Fig.7 Comparisons between power spectra of fluctuating velocities in the wakes without and with the control elements at y/H=±1.3 and x/B = 0,the plate oscillation amplitude β=10° and frequency fe H/V∞=0.054

        2.3 旋渦脫落抑制效果

        設加控制件后實驗測量的尾流脈動速度平均功率譜的主峰值為P,未加控制時測量的功率譜主峰值為P0,則比值P/P0代表尾流旋渦脫落的實際抑制效果:當P/P0<1.0 時,施加控制件能夠減小鎖頻旋渦脫落的規(guī)則脈動能量,效果是正面的;而當P/P0>1.0 時,施加控制件會增大鎖頻旋渦脫落的規(guī)則脈動能量,效果是負面的.

        在圖8 中,兩控制件的橫向位置固定為y/H=±1.3,流向位置在x/B=?0.5,0,0.5 時的控制效果分別如圖8(a)~ 圖8(c)所示.

        圖8 控制件橫坐標y/H=±1.3 時不同振幅下尾流脈動速度功率譜主峰峰值比P/P0 隨振頻feH/V∞的變化Fig.8 Main peak ratio of power spectra of fluctuating velocities P/P0 versus non-dimensional oscillation frequency feH/V∞ for different oscillation amplitudes,with the elements at y/H=±1.3

        當控制件位于板前緣x/B=?0.5 時,振幅β=2.5°和5°時,在振頻范圍feH/V∞=0.005~ 0.08 內,控制件方法具有較好的控制效果.當振幅增大到β=7.5°時,控制效果減弱,但在上述振頻范圍內比值P/P0仍然低于1.0.當振幅增大到β=10°時,在feH/V∞=0.035~ 0.065 范圍內,控制件具有負面效果.

        當控制件位于板中央x/B=0 時,也有相似結果,振幅β=2.5°和5°時控制效果很好,β=7.5°也有正面控制效果,但β=10°時,在大多數振頻下都出現負面控制效果.當控制件位于板尾緣x/B=0.5 時,β=2.5°,5°,7.5° 和10°情況下,具有正面控制效果的振頻范圍分別為feH/V∞=0.007 5~ 0.08,0.01~ 0.075,0.02~ 0.085 和0.03~ 0.075.

        根據Shao 等[42]的研究結果,窄條控制效果與窄條迎風寬度b與板迎風面寬度H*之比b/H*有正相關性.當板與來流具有角度β時,其迎風寬度變?yōu)?/p>

        由該式計算的b/H*隨板攻角的變化如表3 所示.可以看出,隨著板攻角增大,b/H*快速減小,因此控制效果也快速減弱.板旋轉振蕩,其攻角可達振幅值,因此控制效果也隨振幅增大而變差.當β=10°時,b/H*減小到接近失效的臨界值.好在一般顫振振幅都是從小逐漸變大;如果小振幅時具有好的控制效果,則顫振振幅不會增大到不可控的程度.

        表3 窄條迎風寬度比 b/H* 隨攻角β 的變化Table 3 Variation of windward width ratio b/H*with angle of attack β

        2.4 數值模擬與實驗測量的尾流速度剖面對比

        振幅β=2.5°下,在下游x/H=12.5 測量的尾流速度剖面如圖9(a)所示.無控制(控制件位置x/B=∞)時,尾流速度虧損相對較小,尾流相對較窄,流速在y/H=±3 恢復到外流大小.施加控制件于板前緣x/B=?0.5 時,速度虧損加大,尾流變寬為y/H=?4~4.控制件位于板中央x/B=0 和板尾緣x/B=0.5 時,尾流寬度為y/H=?3.5~ 3.5,速度虧損低于控制件在前緣時,但仍高于無控制情況.控制件在板尾緣x/B=0.5 時,測量的速度剖面為W 型,其他控制件位置時速度剖面為V 型.如圖9(b)所示,數值模擬的速度剖面,在各控制件位置下都與實測剖面具有相同的形狀.無控制時,模擬的最大虧損與實測幾乎相同,但模擬的尾流寬度比實測窄.控制件位于前緣時,模擬的最大速度虧損值略大于實測,但尾流寬度兩者幾乎相同.控制件位于板中央時,模擬與實測速度剖面最為接近.控制件位于板尾緣時,模擬的W 型速度剖面,其兩個谷略低于實測值,而尾流中心的峰值略高于實測值,尾流寬度則與實測結果幾乎相同.

        圖9 振幅β=2.5°,振頻fe H/V∞=0.043 下,控制件在橫線y/H=±1.3 上不同x/B 位置時,位于板下游x/H=12.5 處的尾流速度剖面比較Fig.9 Comparisons between measured and simulated velocity profiles sampled at x/H=12.5 in the wakes with elements at y/H=±1.3 and different x/B positions,the oscillation amplitude β=2.5° and frequency fe H/V∞=0.043

        圖10 為振幅增大到β=7.5°的速度剖面情況.如圖10(a)所示,無控制時實測速度虧損較輕,施加控制后速度虧損加大,尤其當控制件位于板前緣時虧損最嚴重.控制件在板尾緣時,速度剖面為W 型,其他控制件位置下(包括無控制),速度剖面為V 型.如圖10(b)所示,數值模擬的各控制件位置下的速度剖面,在形狀和具體數值大小方面,都與實測結果很接近.

        圖10 振幅β=7.5°,振頻feH/V∞=0.054 下,控制件在橫線y/H=±1.3 上不同x/B 位置時,位于板下游x/H=12.5 處的尾流速度剖面比較Fig.10 Comparisons between measured and simulated velocity profiles sampled at x/H=12.5 in the wakes with elements at y/H=±1.3 and different x/B positions,the oscillation amplitude β=7.5° and frequency fe H/V∞=0.054

        圖11(a)為振幅β=2.5°時,在下游x/H=12.5 位置測量的尾流脈動速度均方根(RMS)沿y軸的分布情況.無控制時,尾流核心區(qū)脈動速度較大,分布呈M 狀,在y/H=±1.0 附近達到最高值,隨著離尾流中心線的距離變大,脈動速度迅速減小.加控制件于不同位置時,尾流核心區(qū)脈動速度減小.模擬的脈動速度均方根值分布如圖11(b)所示.無控制時脈動速度分布也呈M 狀,在y/H=±1.0 附近達到最高值,但比實測最高值低.施加控制件于前緣、中央、尾緣后,數值模擬的脈動速度均方根分布形狀均與實測情況相似,但具體數值均比實測值低.

        圖11 振幅β=2.5°,振頻fe H/V∞=0.043 下,控制件在橫線y/H=±1.3 上不同x/B 位置時,位于板下游x/H=12.5 處的尾流脈動速度均方根值比較Fig.11 Comparisons between measured and simulated r.m.s values of fluctuating velocities sampled at x/H=12.5 in the wakes with elements at y/H=±1.3 and different x/B positions,the oscillation amplitude β=2.5° and frequency fe H/V∞=0.043

        圖12(a)為振幅β=7.5°下,在下游x/H=12.5 位置測量的尾流脈動速度均方根沿y軸的分布情況.無控制件與控制件位于不同位置時,均方根值分布都呈M 狀,均方根值有控制時比無控制時低.如圖12(b)所示,數值模擬的脈動速度均方根分布形狀與實測情況一致,但在尾流核心區(qū)模擬的均方根數值低于實測值.實測尾流脈動速度高于模擬值,可能是風洞來流湍流度帶來的影響,具體情況有待進一步研究.

        圖12 振幅β=7.5°,振頻fe H/V∞=0.054 下,控制件在橫線y/H=±1.3 上不同x/B 位置時,位于板下游x/H=12.5 處的尾流脈動速度均方根值比較Fig.12 Comparisons between measured and simulated r.m.s values of fluctuating velocities sampled at x/H=12.5 in the wakes with elements at y/H=±1.3 and different x/B positions,the oscillation amplitude β=7.5° and frequency fe H/V∞=0.054

        2.5 升阻力與力矩控制效果

        由于缺乏實驗條件,振蕩板的升阻力和扭轉力矩無法實測得到.鑒于數值模擬與實驗測量在尾流旋渦流場、時均速度剖面、脈動速度功率譜和脈動速度均方根分布方面的吻合度良好,因此采用數值模擬方法將控制件縱坐標固定在板中央x/B=0 處,改變控制件與板中心軸的距離y/H,研究其對振蕩板升力、阻力及扭矩的控制效果.

        圖13(a)為振幅β=2.5°時,控制前后升力、阻力系數隨時間的變化.其中,C0,C10可以看出,控制前和控制后板的阻力系數脈動都很小,控制后平均阻力略有上升.控制前升力系數隨板的振蕩作周期性變動,控制前變動幅值較大,控制后變動幅值減小約一半.

        圖13 振幅β=2.5°,振頻feH/V∞=0.054 時,無控制件與控制件位于x/B=0,y/H=1.0 時板的升、阻力系數和力矩系數隨時間的變化Fig.13 Time series of force and torque coefficients of the plate in cases without elements and with the elements at x/B=0,y/H=0.5,and β=2.5°,feH/V∞=0.0536

        圖13(b)為控制前、后板的扭轉力矩隨時間的變化.控制前力矩系數作周期性變化,其中,Cm和Cm0變化幅值約為1.8,控制后變化幅值減為約1.0.

        圖14 為振幅β=5°時,窄條對板升力、阻力和力矩系數的控制效果.控制前、后板的阻力系數變化幅度都較小,控制后平均阻力略有上升.控制前板的升力系數和力矩系數變化幅度分別為2.0 和3.8,控制后分別減小為0.9 和2.3.

        圖14 振幅β=5°,振頻feH/V∞=0.054 下,無控制件與控制件位于x/B=0,y/H=±1.0 時板的升力、阻力系數和扭轉力矩系數隨時間的變化Fig.14 Time series of force and torque coefficients of the plate in cases without elements and with the elements at x/B=0,y/H=±1.0,and β=5°,feH/V∞=0.054

        板的扭矩和升力系數脈動幅度的減小,將直接減小顫振導數值,因此提高顫振臨界風速,對顫振抑制具有重要意義.

        將升力系數和力矩系數減去各自的時均值后,得到脈動升力系數和脈動力矩系數.對脈動升力和脈動力矩的時間序列分別求均方根,得到施加與未加控制時脈動升力系數均方根比值Cl/Cl0,和脈動扭矩系數均方根比值Cm/Cm0.當兩個比值均小于1 時,施加窄條對顫振具有正面抑制效果.

        圖15(a)是不同振幅下,有控制與無控制時的脈動扭矩均方根比值Cm/Cm0隨板中央控制件橫向位置y/H的變化情況,其中振頻固定為feH/V∞=0.054.工程中常用的中央控制板方法,是將控制板與橋面接觸或很靠近橋面安裝.但本文結果顯示,控制件很靠近板面時,脈動力矩控制效果不佳.隨著控制件離板面距離增加,脈動升力迅速下降,到y(tǒng)/H=1.0 時,降至最低值.再增大y/H,效果又迅速變差.Cm/Cm0最低值隨振幅的增大而增大,由β=2.5°時的0.6 上升到β=7.5°時的0.7.

        如圖15(b)所示,脈動升力均方根比Cl/Cl0隨控制件位置y/H的變化情況與脈動扭矩系數類似.窄條很靠近板面時,控制效果較差.隨著距離y/H增大,控制效果改善.窄條位置在y/H=0.8 附近時,控制效果最佳;再增大間距,控制效果又變差.β=2.5°,5°,7.5°時,Cl/Cl0最低值分別達到0.1,0.3 和0.5.

        圖15 固定振頻feH/V∞=0.054,不同振幅β 下,有控制和無控制時脈動力距系數均方根比Cm/Cm0 和脈動升力系數均方根比Cl/Cl0 隨中央控制件(x/B=0)橫向位置y/H 的變化圖Fig.15Cm/Cm0 and Cl/Cl0 versus y/H at fixed frequency fe H/V∞=0.054,where y/H is distance between the plate axis and the element on central line x/B=0;Cm/Cm0 (Cl/Cl) is root mean square ratio of fluctuating moment(fluctuating lift) with and without control

        圖16(a) 為振幅β=5°時,脈動扭矩比Cm/Cm0隨中央控制件橫向位置的變化.各振頻下,控制件最佳位置在y/H=0.85~ 1.0 之間.Cm/Cm0最低值隨振頻增大而上升,由feH/V∞=0.043 時的0.57 上升到feH/V∞=0.054,0.064 時的0.63 和0.84.

        圖16(b) 為脈動升力比Cl/Cl0隨控制件距離y/H的變化.當振頻較低,feH/V∞=0.043 時,控制件靠近板面時控制效果最好,隨著距離增大而變差.當振頻稍增大到feH/V∞=0.054,0.064 時,最佳控制效果的位置在y/H=0.85~ 1.0.3 個頻率下,Cl/Cl0最低值在0.2~ 0.25 之間,控制效果都很好.

        由圖15 和圖16 可知,具有較好控制效果的控制件y/H位置區(qū)域較狹窄.對于顫振來說,扭矩控制占首要地位,其次是升力.實際應用中,需根據橋跨結構固有頻率等因素,選擇合適的控制件位置.

        圖17(a)為固定振頻下,扭矩控制效果隨振幅的變化.當控制件很靠近板面y/H=0.75 時,在該頻率下,小振幅時控制效果不佳,但隨振幅增大,控制效果加強.當位置在y/H=0.875~ 1.0 時,扭矩抑制效果最好,且隨振幅改變,抑制效果變化不大.當y/H>1.13 時,抑制效果不佳.

        圖17(b)為固定振頻下,升力控制效果隨振幅的變化.y/H≤ 0.875 時,控制效果隨振幅增大而變的越來越好,到β=7.5°,Cl/Cl0低于0.25.y/H=1.0~1.125 時,Cl/Cl0在0.5~ 0.6 之間,幾乎不隨β改變.y/H>1.25,抑制效果不佳.

        圖17 振幅β 對板脈動扭矩比Cm/Cm0 和脈動升力比Cl/Cl0 控制效果的影響(振頻feH/V∞=0.054.中央控制件(x/B=0)在不同 y/H 位置)Fig.17 Influences of oscillation amplitude β on suppression of fluctuating torque and lift of the plate.(Oscillation frequency feH/V∞=0.0536,and the central elements (x/B=0) at different y/H positions)

        圖18(a)為脈動扭矩比隨feH/V∞的變化.當控制件很靠近板面y/H=0.75 時,脈動扭矩控制效果不是很好,且隨著振頻增大,效果越差.隨著y/H增大到0.875~ 1.0,Cm/Cm0曲線整體上下降很多,抑制效果變好,但Cm/Cm0隨振頻增大而增大的趨勢仍然存在.當y/H=1.125~ 1.5 時,仍有一定控制效果,且不隨振頻改變.工程橋梁的固有頻率都很低,加之顫振臨界風速很高,因此無量綱頻率很低.由于顫振氣動導數的存在,降低結構阻尼和結構剛度,致使橋跨結構固頻率變的更低.因此在低頻下具有好的控制效果,對顫振抑制仍具有重要意義.

        圖18 振頻feH/V∞對板脈動扭矩比Cm/Cm0 和脈動升力比Cl/Cl0 控制效果的影響(振幅β=5°,中央控制件(x/B=0)位于不同 y/H 位置)Fig.18 Influences of oscillation frequency feH/V∞ on suppression of fluctuating torque and lift of the plate (Oscillation amplitude β=5°,and the central elements (x/B=0) at different y/H positions)

        圖18(b)為脈動升力控制效果隨振頻的變化情況.控制件很靠近板面y/H=0.75 時,在低振頻下控制效果很好,但隨著振頻增大,效果越來越差.控制件在y/H=0.875~ 1.125 之間時,總體控制效果較好,特別是y/H=1.0 時,隨振頻增大,效果變得越來越好.

        3 尾流控制機理討論

        尾流控制機理與旋渦生成機理密不可分.靜止柱體尾流渦脫落的形成機制已經有多位作者進行了研究:文獻[43-44]在尾流速度剖面的線性穩(wěn)定性分析中,引入絕對不穩(wěn)定性和對流不穩(wěn)定性概念,認為近尾流中的絕對不穩(wěn)定性是旋渦脫落形成的主要因素;文獻[45-46]的研究證明,只有當近尾流中的絕對不穩(wěn)定性區(qū)足夠大時才能形成尾流的整體不穩(wěn)定和旋渦脫落;文獻[47-49]在柱體尾部引入噴射和抽吸,減小了絕對不穩(wěn)定性區(qū)并抑制了旋渦脫落.

        在振蕩板尾流中,除了絕對不穩(wěn)定性機制以外,還存在另一種旋渦生成機制,稱為信號放大機制[50]:振蕩運動在尾流中發(fā)出周期性擾動信號,該信號在向下游傳播過程中發(fā)展放大,形成旋渦脫落.在振蕩板尾流中,這兩種機制存在競爭:絕對不穩(wěn)定性機制占優(yōu)時形成非鎖頻旋渦脫落,其頻率與柱體振蕩頻率無關.信號放大機制占優(yōu)時形成鎖頻旋渦脫落,其頻率與板振蕩頻率相同.

        根據有關研究,靜止圓柱體尾流初次失穩(wěn)(從二維向三維轉變)的臨界雷諾數約為190,三次失穩(wěn)臨界雷諾數約為260,轉捩為湍流的臨界雷諾數低于400[51-52]).周期振蕩柱體尾流第三次失穩(wěn)的臨界雷諾數為377[53].而三次失穩(wěn)后,雷諾數再增大一個不是很大的定值,尾流將轉捩為湍流[54].在本文研究的雷諾數Re=2800 下,旋轉振蕩板的尾流為湍流.湍尾流的不規(guī)則脈動將使大尺度旋渦的能量向小渦遷移,并可能影響旋渦脫落的形成.施加控制件以后,將在兩個方面影響板尾流:一是改變板尾流的速度剖面,二是增強板尾流不規(guī)則脈動,因此在尾流穩(wěn)定性分析中,這兩個因素都需要考慮,才能正確反映施加控制件后尾流穩(wěn)定性的改變.

        為了研究信號放大問題,將無量綱化的尾流流速的兩個分量分解為平均量和脈動量

        將脈動速度u′,v′ 進一步分解為具有最大增長率的規(guī)則脈動項、其他規(guī)則脈動項和不規(guī)則脈動項3 部分.設規(guī)則擾動幅值變化率相對于其波動周期變化來說很慢,可按照多重尺度方法,引入時間慢變量:τ=εt,其中 ε 為一階小量,則

        式中u1,v1為有最大增長率的規(guī)則脈動幅值,uk和vk(k≥ 2)為其他規(guī)則脈動幅值,他們都是τ和空間坐標的函數.ω1=ω為最大增長的脈動圓頻率(主頻),ωk(k≥ 2)為其他規(guī)則脈動圓頻率.最后一項為不規(guī)則脈動.利用連續(xù)性方程,將動量方程中的壓力項消去,再將得到的單一方程各項分別乘以 eiωt再對時間t積分,認為τ獨立于t,有

        其中T相對于規(guī)則脈動周期和不規(guī)則脈動來說是足夠長的時間,從而使積分收斂;而相對于慢變量來說,εT又足夠短,各規(guī)則脈動幅值變化可忽略。如果將ω看作復數,則脈動幅值變化可化為ω的虛部,式(3)中的積分結果變?yōu)?/p>

        而不規(guī)則項積分后的結果,相當于湍流能譜(寬譜)在主頻率下的值等于0 (如果不為0,可吸收進最大增長的主頻脈動中).

        對不規(guī)則脈動的影響,引入渦黏系數(沒有其他假設),則擾動非線性項部分得到線性化

        式中vT1和vT2分別稱為第一和第二渦黏系數.在平行流假設下,設

        式中,k為波數,u和v為y坐標的函數.忽略V0,忽略U0隨流向坐標的變化,用連續(xù)性方程消去u,最后得到關于v的線性穩(wěn)定性方程

        平行流假設下渦黏系數僅為空間位置y的函數,不隨時間改變.

        根據有關研究[51,55],當采用飽和狀態(tài)(satuated state)流速作為基本流時,得到的主頻率與實驗旋渦脫落頻率很接近,而采用初始速度剖面時,得到的主頻率與實驗差別較大.因此穩(wěn)定性分析中采用飽和速度剖面;而兩個渦黏系數(假設為正實數),可以用飽和狀態(tài)下測量的脈動速度,以及經過濾波后得到的主頻脈動u1和v1,通過式(5)的關系推出.用重復測量結果計算的渦黏系數都非常接近,說明其具有唯一性.

        圖19(a)和圖19(d)分別為無控制和施加控制后各流向站位的速度剖面情況:無控制的速度剖面呈V 型,只有一個低谷;施加雙控制件后,速度剖面存在3 個低谷.

        圖19 振幅 β=2.5°,振頻 feH/V∞=0.043 時,無控制與有控制的尾流速度剖面與渦黏系數分布對比Fig.19 Comparisons of velocity profiles and eddy viscosities between the wakes without elements and with the elements at x/B=0.5,y/H=± 1.3 when β=2.5° and feH/V∞=0.0428

        圖19(b)和圖19(e)分別為無控制和施加控制后的尾流各流向站位的第一渦黏系數情況:無控制時,渦黏系數vT1主要分布在尾流核心區(qū)中y/H=?2~2 內,此外的區(qū)域接近于0,而且在靠近板尾緣時,vT1較低,隨著向下游其值增大.加控制后,vT1最高值增大到無控制時的9 倍,主要分布在y/H=?1~ ?3 和1~ 3 兩個區(qū)域內,其他區(qū)域接近0.

        圖19(c)和圖19(f)分別為無控制和施加控制后的第二渦黏系數在各流向站位的分布情況:無控制和有控制時的分布區(qū)域分別與第一渦黏系數無控制和有控制時的區(qū)域相同;施加控制后vT2最大值是無控制時最大值的4 倍.

        將y/H坐標歸一化為 η(?1≤η≤1),采用滿足邊界條件的車比雪夫組合函數系列φn

        作為基本函數系,將擾動量v展開成2N項級數之和并帶入穩(wěn)定性方程,將函數系的每個 φn乘以方程各項,然后對歸一化坐標進行積分,得到線性代數方程組.線性方程組系數矩陣各特征根的虛部,決定速度剖面的穩(wěn)定性.

        對每個x/H站位,對k=0.025~ 50,每隔0.025取一個波數,求出虛部最大的特征根ω=ωr+iωi,然后再取所有波數中的最大虛部ωimax,作為該x/H站位的最大擾動放大因子.對每個x/H站位,最大放大因子所對應的圓頻率ωrmax,再除以 2π,得無量綱頻率fmax.當擾動的無量綱頻率高于fmax時,該擾動在該x/H站位衰減或增長緩慢.

        圖20(a)為最大擾動放大因子ωimax隨流向站位x/H的變化.無控制時,考慮渦黏性計算的最大放大因子曲線,略低于不考慮渦黏性(在穩(wěn)定性方程中令vT1=vT2=0)時的計算曲線;兩條曲線隨x/H的增大而逐漸下降.施加控制后,不考慮渦黏性時計算的放大因子,在位靠近板尾緣的站位x/H=2.5~ 4 區(qū)間具有很高值,甚至高于未加控制時,與尾流抑制情況不太符合,且當x/H站位離開板較大距離后,放大因子仍保持在0.03 左右,不再減小.施加控制件時,考慮渦黏性計算的放大因子曲線,遠比無控制時考慮渦黏性計算的曲線低.相比不考慮渦黏性情況,考慮渦黏性計算的有控制時的放大因子曲線明顯更低,且隨著x/H站位離開板距離增大,放大因子趨近于0.渦黏系數的引入,體現了控制件方法將大渦能量向小渦和湍流脈動轉移,從而削弱旋渦脫落強度的作用.

        圖20 振幅β=2.5°,振頻feH/V∞=0.043 下,未加控制和加控制件于x/B=0.5,y/H=±1.3 時擾動最高放大因子ωimax 及擾動增長的最大頻率fmax 隨流向站位x/H 的變化曲線對比Fig.20 Maximum perturbation amplification factors and maximum frequencies of amplifying perturbations in the wakes without control and with the control elements at x/B=0.5,y/H=±1.3 when β=2.5° and feH/V∞=0.043

        圖20(b)為擾動放大的頻率曲線情況.無控制時,考慮和不考慮渦黏性影響的計算結果相近,擾動最高放大的頻率在0.12~ 0.14 之間,隨x/H站位的變化很小.施加控制后,不考慮和考慮渦黏性影響計算的擾動最高放大頻率范圍除很靠近板尾緣的小區(qū)域外,都遠低于無控制時,且隨著x/H站位距離的增大向更低頻收縮.考慮渦黏性計算的最高放大頻率曲線比不考慮渦黏性時更低,說明渦黏性對較高頻率的擾動具有更大的抑制作用.上述穩(wěn)定性分析結果,可以解釋施加控制件以后,通過改變尾流速度剖面,增大湍流尾流的渦黏性,從而削弱擾動放大因子,并縮小擾動放大的頻率范圍,進而減小旋渦脫落形成機會,削弱旋渦脫落強度的尾流控制機制.

        4 結論

        在旋轉振蕩板的兩側對稱地放置兩個相同的窄條控制件,研究其對板尾流的影響.控制件分別位于板的前緣、中央和尾緣,控制參數取窄條的橫向坐標y/H.在風洞中用煙線方法及熱線風速儀,對施加控制前后的旋轉振蕩板尾流進行流動顯示及尾流平均和脈動速度測量,并用標準k-ε模型對相應尾流進行數值模擬.不同振頻和振幅下,數值模擬的無控制和有控制渦量場,與流動顯示的尾流旋渦流場具有很高相似度.數值模擬的不同流向站位的平均速度剖面、脈動速度均方根分布和脈動速度功率譜,與測量結果在曲線形狀上很相似,具體數值整體上相差不大,吻合度良好.流動顯示圖說明,在一定橫向位置y/H下,對不同振幅和振頻的振蕩,前緣、中央和尾緣控制件都對板尾流大尺度旋渦脫落具有抑制效果.從控制件上脫落的小尺度旋渦在下游移動和進入板尾流過程中由清晰變得模糊,增強了板尾流的不規(guī)則脈動.控制前后尾流脈動速度功率譜主峰峰值比,代表旋渦脫落抑制效果.根據測量結果,當板振幅β≤ 7.5°時,在一定振頻范圍內,3 種流向位置的控制件都具有較好的尾流控制效果,尤其是中央控制件,有效果的頻率范圍較寬,且功率譜主峰比最低處僅為0.3 左右.數值模擬結果說明,當中央控制件位于一定y/H位置時,流體作用于板的脈動扭轉力矩均方根和脈動升力均方根比無控制時都有大幅下降.但控制效果很好的控制件位置區(qū)域范圍較窄.在振頻方面,也有一定限制:某些情況下,無量綱頻率太高,控制效果不佳.引入第一和第二渦黏系數,對施加控制前后板尾流速度剖面進行線性穩(wěn)定性分析,以探討控制機理.分析結果證明,施加控制件于一定位置后,尾流速度剖面形狀有較大改變,渦黏系數也成倍增大,最終結果,是使擾動最大放大因子隨流向位置的變化曲線比無控制時有很大降低,擾動放大的頻率范圍也大幅收窄.

        猜你喜歡
        旋渦尾流脈動
        新學期,如何“脈動回來”?
        家教世界(2023年25期)2023-10-09 02:11:56
        RBI在超期服役脈動真空滅菌器定檢中的應用
        小心,旋渦來啦
        大班科學活動:神秘的旋渦
        旋渦笑臉
        山間湖
        揚子江(2019年1期)2019-03-08 02:52:34
        地球脈動(第一季)
        飛機尾流的散射特性與探測技術綜述
        雷達學報(2017年6期)2017-03-26 07:53:06
        錐形流量計尾流流場分析
        水面艦船風尾流效應減弱的模擬研究
        91精品国产无码在线观看| 一性一交一口添一摸视频| 狠狠爱无码一区二区三区| 国产精品亚洲A∨无码遮挡| 亚洲午夜经典一区二区日韩| 欧美成人国产精品高潮| 亚洲国产精品sss在线观看av| 精品久久久久久午夜| 日本办公室三级在线观看 | 国产精品乱码一区二区三区| 亚洲中文字幕每日更新| 亚洲中文字幕熟女五十| 丰满精品人妻一区二区| 真人新婚之夜破苞第一次视频| 精品国产av无码一道| 国产99久久久国产精品免费| 狠狠躁夜夜躁人人爽超碰97香蕉| 性生交大片免费看淑女出招| 国产女人91精品嗷嗷嗷嗷| 凹凸世界视频a一二三| 成人中文乱幕日产无线码 | 人妻少妇乱子伦精品| 久久精品re| 亚洲av乱码国产精品观看麻豆| 欧美激情视频一区二区三区免费| 亚洲精品97久久中文字幕无码| 白丝美女被狂躁免费视频网站| 亚洲av手机在线播放| 精品国产一区二区三区av性色| 91精品国产综合成人| 国产一区二区毛片视频| 麻豆国产一区二区三区四区| 国产性一交一乱一伦一色一情| 特级毛片全部免费播放a一级| 精品人妻久久一日二个| 国产好大好硬好爽免费不卡| 日本中文字幕不卡在线一区二区| 国产三级韩三级日产三级| 亚洲欧美国产精品久久| 亚洲成a人片在线看| 日韩男女av中文字幕|