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

    高水頭混流式水輪機減負荷瞬態(tài)流動特性研究

    2022-12-05 06:29:00孫龍剛徐卓飛郭鵬程鄭小波
    水利學(xué)報 2022年11期
    關(guān)鍵詞:導(dǎo)葉轉(zhuǎn)輪水輪機

    孫龍剛,徐卓飛,郭鵬程,2,鄭小波,2

    (1.西安理工大學(xué) 水利水電學(xué)院,陜西 西安 710048;2.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國家重點實驗室,陜西 西安 710048)

    1 研究背景

    水電由于具有快速高效的啟停響應(yīng),因而被用來調(diào)節(jié)電網(wǎng)容量及頻率的動態(tài)平衡。近年來清潔能源在電網(wǎng)中占比持續(xù)提高,以風(fēng)能和太陽能為代表的清潔能源的不定時并網(wǎng)對電網(wǎng)造成劇烈的不穩(wěn)定沖擊,此時水電機組必然經(jīng)歷快速且頻繁的工況轉(zhuǎn)換如減負荷過程以平衡電網(wǎng)參數(shù)[1-2]。水輪機運行模式的轉(zhuǎn)變,通常伴隨著回流、流動分離、空化、尾水管渦帶、葉道渦、卡門渦等不穩(wěn)定流動現(xiàn)象[3-6],從而在機組內(nèi)部誘發(fā)具有繁雜幅頻特征的壓力脈動。高水頭水輪機由于所受動態(tài)液體力更大,且水輪機減負荷過程中轉(zhuǎn)輪對來流的不適應(yīng)性增強,機組將承受更加復(fù)雜的不平衡作用力,由此產(chǎn)生的交變應(yīng)力載荷會進一步加劇機組關(guān)鍵部件的疲勞破壞和磨損,縮短機組壽命,降低機組水力性能,從而影響電站的安全穩(wěn)定運行[7-9]。因此,開展減負荷過程中高水頭水輪機動態(tài)響應(yīng)研究,對進一步促進水電與其他清潔能源的多能互補具有重要的現(xiàn)實意義。

    模型試驗與數(shù)值模擬是明確工況轉(zhuǎn)換過程中水輪機內(nèi)部不穩(wěn)定流動特性及其誘發(fā)水力振動的重要手段。Trivedi等[10]試驗研究了變轉(zhuǎn)速工況下功率急速上升和下降過程中高水頭混流式水輪機內(nèi)部的壓力脈動響應(yīng),發(fā)現(xiàn)動靜干涉頻率對應(yīng)的壓力脈動幅值在負荷變化期間相對較小,而在瞬態(tài)過程結(jié)束時提高了30倍。此外,快速啟停策略下水輪機轉(zhuǎn)輪葉片最大壓力脈動幅值是緩慢啟停策略幅值的1.5倍[11]。周勤等[12]通過數(shù)值方法研究水泵水輪機,發(fā)現(xiàn)甩負荷過程中轉(zhuǎn)輪進口回流是不穩(wěn)定壓力脈動及不均勻湍動能生成的重要原因,葉片扭矩和徑向力的最大幅值分別達初始值的10倍和60倍。Fu等[13]基于三維大渦模擬技術(shù)捕捉到水泵水輪機內(nèi)部的異常低頻壓力脈動,研究發(fā)現(xiàn)該種壓力脈動的誘發(fā)與水輪機內(nèi)部水錘現(xiàn)象、轉(zhuǎn)輪進口局部回流以及無葉區(qū)擋水環(huán)有關(guān)。基于一維管道與三維數(shù)值模擬耦合方法,陳秋華等[14]發(fā)現(xiàn)以部分負荷工況為初始工況,水泵水輪機飛逸過程中產(chǎn)生的轉(zhuǎn)輪徑向力和壓力脈動相對額定工況更為劇烈。Goyal等[15-16]通過水輪模型試驗發(fā)現(xiàn),甩負荷過程中水輪機尾水管內(nèi)形成的螺旋渦帶可分解為表征軸向運動的突進模式和表征軸向運動的旋轉(zhuǎn)模式,且突進模式比旋轉(zhuǎn)模式提前出現(xiàn),而增負荷過程中柱狀渦帶的形成主要受負荷變化過程中尾水管內(nèi)駐點、反向流動及回流區(qū)的直接影響。Hosseinimanesh等[17]的研究顯示,混流式水輪機運行在空載條件下,尾水管內(nèi)部的軸向反向流動和周向切向流動之間交界面上出現(xiàn)的不穩(wěn)定強剪切層,是造成渦結(jié)構(gòu)破裂和壓力波動向上游傳播至轉(zhuǎn)輪葉片的主要原因。Liu等[18]數(shù)值研究結(jié)果表明,水泵水輪機甩負荷過程中,尾水管壁面附近出現(xiàn)大量大尺度渦結(jié)構(gòu),引發(fā)低頻壓力脈動。Liu等[19]對原型水泵水輪機甩負荷條件下渦帶的形成過程進行了分析,發(fā)現(xiàn)尾水管中的渦帶出現(xiàn)在零力矩運行之前,隨著流量的減小,尾水管中的渦帶強度增加。作者前期研究顯示[20-21],水輪機甩負荷過程中,尾水管內(nèi)的渦旋強度顯著提升,而飛逸過程中的能量耗散主要發(fā)生在轉(zhuǎn)輪和尾水管內(nèi)部。此外,F(xiàn)u等[22-24]基于數(shù)值模擬方法研究了間隙、空化以及水聲學(xué)對水泵水輪機工況轉(zhuǎn)換過程的影響。

    綜上所述,水輪機工況轉(zhuǎn)換會顯著提升其內(nèi)部的壓力脈動幅值,且直接影響水輪機尾水管內(nèi)的渦流結(jié)構(gòu)形態(tài)及其演化過程。為進一步明確工況轉(zhuǎn)換條件下水輪機特征參數(shù)的動態(tài)響應(yīng)、水輪機內(nèi)部壓力脈動及渦流演化規(guī)律,本文基于動網(wǎng)格技術(shù)對一高水頭混流式模型水輪機負荷減小過程進行了瞬態(tài)數(shù)值研究,分析了水輪機壓力和速度演變規(guī)律、尾水管內(nèi)渦流演化特性,并建立轉(zhuǎn)輪軸向力與不穩(wěn)定壓力脈動之間的關(guān)聯(lián)。

    2 數(shù)值模擬策略

    2.1 高水頭模型水輪機本文以挪威科技大學(xué)(Norwegian University of Science and Technology-NTNU)公開的Francis99模型水輪機為研究對象,開展水輪機負荷減小過程的瞬態(tài)數(shù)值研究[25-27]。Francis99模型水輪機三維幾何結(jié)構(gòu)如圖1所示,由進口至出口分別為帶有14個嵌入式固定導(dǎo)葉的蝸殼、葉片數(shù)為28的活動導(dǎo)葉、由15個分流葉片和15個長葉片組成的轉(zhuǎn)輪以及彎肘型尾水管。該水輪機模型與原型轉(zhuǎn)輪出口直徑分別為0.349 m和1.778 m,對應(yīng)的原型機和模型機水頭分別為12.0 m和377 m。圖2所示為Francis99模型水輪機負荷減小過程中活動導(dǎo)葉開度隨時間變化曲線,由圖可知,負荷減小過程涉及到3種運行工況,分別為時段t=0至t=1.0 s,活動導(dǎo)葉開度為的9.84°的最優(yōu)工況(Best Efficiency Point-BEP),時段t=3.5至t=10.0 s,導(dǎo)葉開度為6.72°的部分負荷工況(Part Load-PL)以及時段t=1.0至t=3.5 s,導(dǎo)葉開度由9.84°線性減小至6.72°的負荷減小過程(BEP to PL)。此外,整個負荷變化過程中,水輪機轉(zhuǎn)輪的轉(zhuǎn)速均為333.0 r/min。

    圖1 Francis99模型水輪機

    圖2 活動導(dǎo)葉開度變化規(guī)律

    為了對數(shù)值方法及模擬控制進行驗證,NTNU提供了部分模型試驗結(jié)果,其中壓力及速度測量位置如圖3所示[26-27]。圖3中共布置了6個壓力測點,其中VL02位于活動導(dǎo)葉與轉(zhuǎn)輪之間的無葉區(qū),用以捕捉動靜干涉效應(yīng);轉(zhuǎn)輪葉片壓力面中間位置布置測點P42,出口側(cè)靠近上冠處布置測點P71,用來分析轉(zhuǎn)輪內(nèi)部的壓力脈動特征;同理,轉(zhuǎn)輪葉片吸力面靠近上冠處布置測點S51;尾水管錐管段間隔180°布置了兩個壓力測點DT5和DT6,用來分析尾水管壓力場隨負荷變化規(guī)律。此外,試驗中通過PIV(Particle Image Velocimetry)技術(shù)對尾水管錐管段3條線段L1、L2和L3進行了速度測量,其中L1與L2為垂直距離相差120 mm且與水輪機軸線垂直的水平線,其與x軸之間的角度均為79°;L3為與水輪機軸線重合的豎直線。速度測量中,對L1和L2上的28個位置進行測量,而L3上的測點數(shù)為19。

    圖3 壓力與速度測點布置示意圖

    2.2 導(dǎo)葉網(wǎng)格變形控制水輪機工況轉(zhuǎn)換是通過改變活動導(dǎo)葉開度進而控制水輪機流量來實現(xiàn)的,因此水輪機負荷改變數(shù)值研究中的一項重要內(nèi)容為按照一定的導(dǎo)葉運動規(guī)律動態(tài)調(diào)整活動導(dǎo)葉的開度且保證一定的網(wǎng)格精度。為了實現(xiàn)上述目標,本文采用一種幾何結(jié)構(gòu)布爾運算與網(wǎng)格重構(gòu)相結(jié)合的方法為活動導(dǎo)葉模擬域提供網(wǎng)格的自動更新,如圖4和圖5所示。當(dāng)導(dǎo)葉開度變化時,導(dǎo)葉區(qū)域的網(wǎng)格發(fā)生變形,求解器按照預(yù)設(shè)的網(wǎng)格質(zhì)量判斷準則,確定是否執(zhí)行網(wǎng)格重新劃分。導(dǎo)葉運動過程中,網(wǎng)格質(zhì)量判斷準則為同時滿足網(wǎng)格單元質(zhì)量大于0.3且網(wǎng)格偏斜角小于85°。若不滿足該條件,則首先獲取當(dāng)前時間步的導(dǎo)葉位置,然后經(jīng)過導(dǎo)葉與背景區(qū)域兩者之間的布爾運算獲得新的導(dǎo)葉位置,最后以相同的網(wǎng)格設(shè)置參數(shù)進行網(wǎng)格重新劃分。

    圖4 導(dǎo)葉運動控制流程圖

    圖5 活動導(dǎo)葉運動示意圖

    2.3 求解設(shè)置本文基于ANSYS Fluent進行水輪機負荷變化的數(shù)值計算,湍流模型采用SSTk-ω模型[28-29],該模型改進了近壁區(qū)低雷諾數(shù)流動的處理方式,能較好處理帶有逆壓梯度、流動分離的復(fù)雜流動問題。壁面邊界層則采用全y+壁面函數(shù)處理方式,這種方式的優(yōu)點在于使用混合壁面函數(shù),對精細網(wǎng)格仿效低y+壁面處理,而對粗糙網(wǎng)格仿效高y+壁面處理。通過剛體運動來處理轉(zhuǎn)輪域網(wǎng)格的旋轉(zhuǎn)運動,所有壁面均設(shè)置為無滑移,對流項采用二階迎風(fēng)格式,時間離散則采用高精度的二階格式,收斂標準設(shè)置為殘差小于0.0001。按照M?ssinger等[30]和Sotoudeh等[31]的研究,針對Francis99模型瞬態(tài)數(shù)值研究,瞬態(tài)時間步長設(shè)置為小于轉(zhuǎn)輪旋轉(zhuǎn)4°所用的時間能夠獲得較為可靠的數(shù)值結(jié)果。為了進一步提高數(shù)值精度,捕捉負荷變化過程中的壓力脈動及渦流特性,本文采用對應(yīng)于轉(zhuǎn)輪旋轉(zhuǎn)2°所用的時間,經(jīng)計算瞬態(tài)時間步長為0.001 s,且每個時間步內(nèi)迭代15次。首先進行BEP工況下的穩(wěn)態(tài)計算,并以此為初場進行BEP工況下的非穩(wěn)態(tài)計算,以轉(zhuǎn)輪旋轉(zhuǎn)4圈的非穩(wěn)態(tài)計算結(jié)果作為負荷減小瞬態(tài)計算的初始值。數(shù)值計算進口給定總壓,出口指定靜壓,這種邊界條件的設(shè)置有效保證了數(shù)值計算與模型試驗水輪機進出口壓差的一致性。按照NTNU提供的進口靜壓P1、流量Q和凈水頭H的試驗數(shù)據(jù),可通過式(1)至(3)獲得蝸殼進口總壓Pin和尾水管出口靜壓條件Pout,

    (1)

    (2)

    (3)

    式中:P1為蝸殼進口靜壓,Pa;Pout為尾水管出口靜壓,Pa;Q為流量,m3/s;ρ為密度,取999.8 kg/m3;g為重力加速度,取9.82 m/s2;Ain為蝸殼進口橫截面面積;Aout為尾水管出口橫截面面積;vin為蝸殼進口速度,m/s;vout為尾水管出口速度,m/s;z為水輪機進出口高程差,m。

    此外,本文對NTNU提供的進出口壓力數(shù)據(jù)進行非線性曲線擬合處理,取其擬合結(jié)果作為瞬態(tài)模擬的進出口邊界條件,如圖6中點劃線所示。

    圖6 進出口壓力邊界條件

    2.4 計算域離散及網(wǎng)格無關(guān)性驗證本文采用多面體網(wǎng)格對計算域進行離散,多面體網(wǎng)格的最大優(yōu)點在于其一個控制體周圍有大約10個相鄰單元,因而能以相對較少的控制體獲得更加精確的梯度和局部流動分離解[32-33]。圖7所示為計算域網(wǎng)格示意圖,其中固壁面如導(dǎo)葉、轉(zhuǎn)輪葉片等通過棱柱層網(wǎng)格進行邊界加密處理。

    圖7 計算域多面體網(wǎng)格示意圖

    為消除網(wǎng)格數(shù)量對計算結(jié)果的影響,本文采用美國機械工程師協(xié)會(ASME)推薦的網(wǎng)格收斂指數(shù)(GCI)進行網(wǎng)格離散誤差的估計[34-36]。GCI是一個具有95%置信區(qū)間、表示兩個對比網(wǎng)格中更密網(wǎng)格與漸進值之間距離的指標,并且預(yù)測進一步的網(wǎng)格細化對求解的影響。GCI網(wǎng)格無關(guān)性驗證需要3套不同數(shù)目的網(wǎng)格,分別為細密網(wǎng)格(Fine)、中等網(wǎng)格(Medium)和粗糙網(wǎng)格(Coarse),計算的近似相對誤差ea,外推相對誤差eext以及網(wǎng)格收斂指數(shù)GCI:

    (4)

    (5)

    (6)

    式中:φ為計算所選擇的關(guān)鍵變量;r為網(wǎng)格加密因子;p為采用定點迭代法計算的表觀級數(shù);下標1、2、3分別對應(yīng)于網(wǎng)格方案Fine、Medium、Coarse;上標21和32分別表示網(wǎng)格Fine相對于Medium以及網(wǎng)格Medium相對于Coarse的計算值。

    表1列出了BEP工況下數(shù)值計算離散誤差的計算過程及結(jié)果統(tǒng)計。其中,N1、N2、N3為三種不同密度下的網(wǎng)格數(shù),分別對應(yīng)Fine、Medium和Coarse網(wǎng)格方案。本文選擇水輪機水頭H(單位為m)和轉(zhuǎn)輪扭矩T(單位為N·m)作為網(wǎng)格無關(guān)性測試的關(guān)鍵變量。由表1統(tǒng)計結(jié)果知,三種密度網(wǎng)格的數(shù)值解以漸進形式收斂,表明網(wǎng)格加密有利于平均流場的求解。對Fine和Medium網(wǎng)格而言,計算的水頭不確定度分別為0.072%和0.18%,轉(zhuǎn)輪扭矩不確定度為0.066%和0.25%,表明Fine和Medium網(wǎng)格方案獲得的數(shù)值解誤差較小。為了平衡計算精度與計算資源之間的關(guān)系,本文最終選擇了Medium網(wǎng)格進行數(shù)值計算,網(wǎng)格總數(shù)目為1451萬,其中蝸殼和固定導(dǎo)葉網(wǎng)格數(shù)為292萬,活動導(dǎo)葉為284萬,轉(zhuǎn)輪為514萬,尾水管為361萬。圖8顯示了Medium網(wǎng)格方案BEP工況轉(zhuǎn)輪、尾水管的y+值分布云圖,由圖中可知,轉(zhuǎn)輪壁面最大y+值為11.893,尾水管壁面最大y+值為7.061,且轉(zhuǎn)輪壁面y+值主要范圍在2.5以下,而尾水管壁面y+值主要范圍在1.1左右,可以認為基本適用于SSTk-ω湍流模型。

    表1 數(shù)值計算離散誤差及不確定性統(tǒng)計表

    圖8 轉(zhuǎn)輪及尾水管壁面y+分布

    3 結(jié)果與討論

    3.1 外特性分析為驗證本文所采用數(shù)值策略的可靠性,圖9顯示了數(shù)值模擬與試驗測試獲得的水輪機外特性水頭H和流量Q的對比結(jié)果。

    圖9 負荷減小過程水輪機水頭與流量變化

    圖9結(jié)果顯示,數(shù)值預(yù)測的水輪機水頭與試驗測試趨勢平均值一致性較高。負荷減小過程中,水頭在BEP工況下保持恒定,導(dǎo)葉閉合后,水頭快速上升。t=3.5 s導(dǎo)葉停止閉合,水輪機開始運行在PL工況,此時水頭繼續(xù)上升,在t=4.2 s達到最高值,隨后緩慢下降。水頭變化結(jié)果表明,導(dǎo)葉的閉合對水輪機水頭有一定的提升,且導(dǎo)葉停止運動后這種影響仍會在一定時間內(nèi)持續(xù)。數(shù)值預(yù)測的流量與試驗測試結(jié)果之間存在一定差異,數(shù)值結(jié)果顯示為與導(dǎo)葉開度變化規(guī)律一致的線性變化,而試驗值則為非線性變化規(guī)律,且導(dǎo)葉閉合過程t=1.0 s至t=3.5 s時段內(nèi),試驗記錄的流量值變化率小于數(shù)值結(jié)果。造成試驗與數(shù)值之間差異的主要原因在于試驗測試的流量計存在一定的時間延遲[27]。盡管數(shù)值解與試驗測試存在一定差異,但可通過對比負荷減小過程中定導(dǎo)葉開度BEP階段、PL階段流量值與對應(yīng)穩(wěn)態(tài)試驗測試結(jié)果之間的差異來量化本文數(shù)值流量結(jié)果的可信度。負荷減小過程中,BEP與PL工況數(shù)值計算獲得的流量分別為0.193 24 m3/s和0.136 27 m3/s,對應(yīng)穩(wěn)態(tài)試驗條件流量分別為0.199 59 m3/s和0.139 62 m3/s,僅僅相差3.2%和2.4%,表明數(shù)值計算的流量準確性較高。此外,水輪機流量調(diào)節(jié)是通過改變活動導(dǎo)葉開度實現(xiàn)的,流量變化趨勢與導(dǎo)葉開度密切相關(guān),本文中數(shù)值模擬與試驗測試均顯示了隨導(dǎo)葉閉合流量逐漸減小的變化趨勢,因此可認為本文數(shù)值計算的流量變化規(guī)律具有較高的可信度。綜上所述,采用本文數(shù)值計算的模擬策略,能較好地反映負荷減小過程中水輪機外特性參數(shù)的動態(tài)響應(yīng)。

    3.2 壓力脈動分析為進一步分析負荷減小過程中水輪機內(nèi)部壓力脈動特性,圖10顯示了無葉區(qū)和尾水管內(nèi)VL02、DT5與DT6的壓力脈動系數(shù)Cp,其中Cp為按照式(7)計算的無量綱數(shù),

    (7)

    分析圖10結(jié)果可知,數(shù)值模擬獲得的無葉區(qū)壓力脈動變化規(guī)律與試驗測試結(jié)果吻合地較好,壓力系數(shù)在定導(dǎo)葉開度BEP及PL工況下顯示為均勻波動,而隨著導(dǎo)葉的閉合壓力系數(shù)線性減小。幅值方面,壓力系數(shù)數(shù)值結(jié)果相對于試驗值偏低,這與試驗測試過程中的機械振動、電磁振蕩以及數(shù)值模擬采用的網(wǎng)格、湍流模型和時間步長等有關(guān)。此外,由初始至終止時刻,數(shù)值模擬預(yù)測的壓力系數(shù)幅值保持一致,而試驗結(jié)果在負荷減小過程及PL工況下的壓力脈動幅值要更高。

    圖10 數(shù)值與試驗壓力系數(shù)對比

    尾水管內(nèi)的壓力系數(shù)結(jié)果顯示,數(shù)值模擬低估了尾水管錐管段的壓力脈動幅值,然而,其壓力系數(shù)變化趨勢與試驗值一致性較高,因而在定性分析尾水管內(nèi)的流動特性上具有較高的可靠性。BEP工況下,轉(zhuǎn)輪出口接近法向出流,數(shù)值壓力系數(shù)與試驗壓力系數(shù)平均值均保持為一個固定值,導(dǎo)葉閉合后,壓力系數(shù)迅速下降,隨后隨導(dǎo)葉開度的減小逐漸增加。與BEP工況相似,導(dǎo)葉停止運動進入PL工況,壓力系數(shù)首先迅速上升,隨后緩慢下降,數(shù)值結(jié)果在時刻t=7.2 s以后出現(xiàn)一個顯著的低頻、高幅值周期性脈動,而相同的周期性脈動在試驗過程中約在t=6.5 s后出現(xiàn)。低頻、高幅值周期性壓力系數(shù)的出現(xiàn),與部分負荷工況下尾水管內(nèi)螺旋形尾水管渦帶的出現(xiàn)有關(guān)。

    3.3 轉(zhuǎn)輪軸向水推力分析軸向水推力不僅是水輪機組推力軸承載荷的重要組成部分,而且其脈動特性對會直接引起廠房結(jié)構(gòu)的振動,因而是水輪機設(shè)計和運行過程中的一個重要參數(shù)[37]。為進一步分析軸向水推力相對水輪機負荷變化的響應(yīng)關(guān)系,圖11顯示了負荷減小過程中水輪機軸向水推力隨時間變化規(guī)律,其中正值表示數(shù)值模擬獲得的軸向水推力與參考坐標系z軸方向一致。圖中結(jié)果顯示,BEP工況下軸向水推力絕對值最高且保持在一個常數(shù)附近。軸向水推力變化與導(dǎo)葉運動規(guī)律比較一致,即隨著導(dǎo)葉的閉合而減小,且在導(dǎo)葉開始運動后,軸向水推力快速降低,而在導(dǎo)葉停止閉合后進入PL工況,軸向水推力首先快速上升,然后在一定時間范圍內(nèi)緩慢下降,最后做周期性脈動。

    圖11 負荷減小過程中軸向水推力隨時間變化曲線

    為確定影響軸向水推力的重要因素,本文通過相干分析方法建立軸向水推力與不穩(wěn)定壓力脈動之間的聯(lián)系[38],如圖12所示,圖中分別顯示了軸向水推力與無葉區(qū)測點VL02、轉(zhuǎn)輪域測點P42以及尾水管測點DT6的相干系數(shù)。圖12中橫坐標f/fn定義為一個無量綱頻率比,為計算頻率f與水輪機轉(zhuǎn)頻fn之間的比值。縱坐標為相干系數(shù),用來衡量兩個變量之間的相關(guān)程度,相干系數(shù)越大,表明兩個變量之間關(guān)聯(lián)程度越高。分析圖12可知,除較低頻率范圍及導(dǎo)葉、葉片通過頻率28fn和30fn附近,軸向力相對于DT6的相干性系數(shù)均高于VL02和P42。相干系數(shù)大意味著兩個信號之間的關(guān)聯(lián)程度高,因此認為軸向力Fz與尾水管中DT6的壓力脈動的相關(guān)度要顯著高于無葉區(qū)和轉(zhuǎn)輪內(nèi),表明尾水管內(nèi)部非穩(wěn)態(tài)壓力波動向上游的傳播是造成轉(zhuǎn)輪軸向力波動的重要因素。

    圖12 軸向力與測點VL02、P42、DT6壓力信號相關(guān)性

    3.4 尾水管內(nèi)速度變化情況分析螺旋形尾水管渦帶的形成,與錐管段內(nèi)的軸向速度及切向速度密切相關(guān)[39-40]。圖13和圖14分別為瞬態(tài)過程尾水管L1至L3上的軸向和切向速度分布,圖中軸向速度沿z軸正向為負,切向速度沿轉(zhuǎn)輪轉(zhuǎn)向為正,圖中橫坐標r/r0定義為一個無量綱距離比,為速度測量線L1、L2、L3上坐標位置相對于該線段中間位置距離的比值。由于NTNU未提供切向速度數(shù)據(jù),因此本文僅給出數(shù)值結(jié)果。為便于對不同工況的速度變量進行分析,本文將速度變化分為3個階段,其中BEP工況為階段Ⅰ,導(dǎo)葉開度變化過程BEP至PL為階段Ⅱ,PL工況為階段Ⅲ。

    圖13結(jié)果顯示,水平線L1及L2上數(shù)值預(yù)測的軸向速度與試驗結(jié)果比較吻合,豎直線L3上數(shù)值與試驗結(jié)果存在一定差異,但二者之間的變化趨勢比較一致。第Ⅰ階段,轉(zhuǎn)輪出口水流接近法向出流,軸向速度與主流方向相同,泄水錐下方L1與L2上的軸向速度高于尾水管錐管段壁面附近,而L3沿流向保持一致。導(dǎo)葉閉合的第Ⅱ階段,L1與L2上測得的軸向速度在-0.5

    圖13 軸向速度分布

    分析圖14結(jié)果可知,L1與L2上的切向速度變化比較一致,在第Ⅰ階段均與轉(zhuǎn)輪轉(zhuǎn)向一致,且轉(zhuǎn)軸中心處絕對值較小。第Ⅱ階段,轉(zhuǎn)軸中心附近發(fā)生切向速度的反轉(zhuǎn),且隨時間推移反向切向速度范圍逐漸擴大。第Ⅲ階段,轉(zhuǎn)軸中心附近的正向切向速度區(qū)隨導(dǎo)葉的閉合逐漸增加,而在這個正向速度區(qū)外側(cè),均顯示為反向切向速度。之后,切向速度也表現(xiàn)以轉(zhuǎn)軸為中心線,正向和負向交替出現(xiàn)的切向速度區(qū)。豎直線L3上的切向速度在第I和第Ⅱ階段變化較小,而在第Ⅲ階段后期切向速度出現(xiàn)方向的交替變化,且靠近轉(zhuǎn)輪中心測點經(jīng)歷更多的切向速度方向變化。

    圖14 切向速度分布

    3.5 尾水管渦帶形成分析由BEP工況過渡至PL工況,水輪機尾水管內(nèi)的壓力、速度均顯示了一定的低頻振蕩,這對尾水管渦帶的形成有一定影響。為進一步深入理解水輪機負荷變化過程中尾水管內(nèi)渦結(jié)構(gòu)演化及尾水管渦帶的形成過程,圖15至圖17分別顯示了階段Ⅰ至階段Ⅲ尾水管內(nèi)軸向速度云圖、速度矢量分布及以Q準則顯示的渦結(jié)構(gòu),圖中Q準則閾值q=700 s-2,軸向速度為負時,表示軸向速度方向與參考坐標系z軸相反,與主流方向一致。由于BEP工況不同時刻尾水管內(nèi)流場結(jié)構(gòu)比較接近,因此本文只給出了t=1.0 s時刻的結(jié)果。

    圖15 階段Ⅰ尾水管軸向速度及渦結(jié)構(gòu)

    圖16 階段Ⅱ尾水管軸向速度及渦結(jié)構(gòu)

    圖17 階段Ⅲ尾水管軸向速度及渦結(jié)構(gòu)

    圖15速度分布結(jié)果顯示,BEP工況下尾水管內(nèi)為幾乎與主流方向相同的純軸向流動的穩(wěn)定流場,軸向速度在錐管段中心附近較高,而在錐管段壁面附近較低。由于泄水錐的壁面效應(yīng),在其正下方出現(xiàn)范圍較小的局部回流區(qū)??梢暬瘻u結(jié)構(gòu)等值面結(jié)果未出現(xiàn)自由渦結(jié)構(gòu),僅僅顯示為附著于尾水管進口、泄水錐表面以及彎肘段的片狀渦結(jié)構(gòu),這主要是由于固壁面附近速度梯度較高所致,而且這種片狀渦結(jié)構(gòu)在整個模擬過程的不同時刻均存在。因此,BEP工況下尾水管內(nèi)流動光順,未誘發(fā)不穩(wěn)定渦結(jié)構(gòu)。

    第Ⅱ階段,水輪機活動導(dǎo)葉由9.84°線性關(guān)閉至6.72°,水輪機流量逐漸減小,轉(zhuǎn)輪進口沖角增大。分析圖16結(jié)果可知,t=1.8 s,尾水管內(nèi)的軸向速度由于流量的降低而減小,特別是轉(zhuǎn)軸附近區(qū)域,此時水流方向仍然為軸向流動,未出現(xiàn)回流及分離區(qū)。然而,尾水管錐管段內(nèi)出現(xiàn)了軸對稱的自由渦結(jié)構(gòu),表明水輪機負荷的減小,激發(fā)了尾水管內(nèi)部的不穩(wěn)定流動。t=2.7 s,隨著流量的減小,尾水管內(nèi)的軸向速度進一步降低,錐管段出口與肘管段進口中心區(qū)域出現(xiàn)由肘管段外側(cè)指向內(nèi)側(cè)的流動偏移區(qū)。此時,尾水管內(nèi)的自由渦結(jié)構(gòu)強度顯著增加,由泄水錐延伸至肘管段進口處。導(dǎo)葉繼續(xù)閉合至t=3.0 s,流動偏移強度迅速增強,其中軸向速度流線方向與尾水管軸線呈較大角度,其范圍擴大至彎肘段中部。

    此外,泄水錐下游出現(xiàn)反向軸向速度而形成回流區(qū),此時的渦強度相對t=2.7 s有所增強,但仍然保持為對稱結(jié)構(gòu)。導(dǎo)葉終止閉合時刻t=3.5 s,尾水管內(nèi)的軸向速度持續(xù)降低,流動偏移區(qū)域與回流區(qū)域范圍均顯著擴大,表明水輪機流量的減小使尾水管內(nèi)的水力不穩(wěn)定性進一步增強。然而,圖16中顯示的渦結(jié)構(gòu)沿軸向迅速收縮,而在徑向其渦核半徑增加。綜上所述,水輪機負荷變化過程中,尾水管內(nèi)的軸向速度逐漸降低,誘發(fā)不穩(wěn)定的流動偏移及形成回流區(qū)。可視化渦結(jié)構(gòu)首先顯示為強度較低的軸對稱結(jié)構(gòu),隨后渦強度沿軸向及徑向均顯著增強,最后沿軸向迅速收縮而沿徑向強度增加,而在負荷變化的整個過程中渦形態(tài)均顯示為對稱結(jié)構(gòu)。

    第Ⅲ階段,水輪機運行在定導(dǎo)葉開度PL工況下。t=3.7 s,尾水管內(nèi)的軸向速度分布與t=3.5 s時刻比較接近,但轉(zhuǎn)軸中心區(qū)域的回流強度提高且軸對稱渦結(jié)構(gòu)沿軸向收縮。t=5.5 s時,尾水管錐管段及彎肘段的低速區(qū)范圍顯著增加,對應(yīng)地,尾水管內(nèi)的流動偏移區(qū)域與回流區(qū)域范圍被擴大,且在轉(zhuǎn)軸中心左側(cè)形成兩個明顯的漩渦。然而,相對于t=3.7 s,該時刻的自由渦強度被顯著抑制,僅僅殘留附著在泄水錐表面的部分渦結(jié)構(gòu)。t=7.0 s,尾水管錐管段低速區(qū)內(nèi)出現(xiàn)了多個漩渦區(qū),渦結(jié)構(gòu)強度較上一時刻顯示了較大的差別,該時刻的渦強度有所增強,但渦形態(tài)演化為雙螺旋狀偏心結(jié)構(gòu)。t=8.5 s,錐管段流態(tài)紊亂程度加劇,轉(zhuǎn)軸中心的低速區(qū)不再關(guān)于尾水管中心線對稱,其在尾水管內(nèi)顯示為“S”狀,且靠近壁面的低速區(qū)處出現(xiàn)明顯的漩渦結(jié)構(gòu),對水流形成一定的擠壓作用而使壁面附近的速度提高。Q準則顯示的渦形態(tài)為強度較高、單一的螺旋狀偏心渦結(jié)構(gòu),為典型的部分負荷尾水管渦帶,t=8.5 s之后,該螺旋狀偏心渦帶做與轉(zhuǎn)輪轉(zhuǎn)向相同的旋轉(zhuǎn)運動且其形態(tài)保持固定。綜合以上分析可知,PL工況下,尾水管內(nèi)的回流范圍顯著增加,渦形態(tài)首先為軸對稱形,然后沿軸向劇烈收縮,隨后單一的柱狀渦結(jié)構(gòu)演化為強度升高的雙螺旋狀渦結(jié)構(gòu),最終該雙螺旋狀渦結(jié)構(gòu)合并為單一的螺旋狀偏心渦帶并保持其形態(tài)。

    4 結(jié)論

    本文采用基于多面體網(wǎng)格的動網(wǎng)格技術(shù)對Francis99高水頭混流式模型水輪機負荷減小過渡過程進行了數(shù)值計算,分析了負荷變化對水輪機外特性參數(shù)、不同位置壓力脈動及尾水管內(nèi)渦流結(jié)構(gòu)的影響,主要結(jié)論如下:

    (1)基于多面體網(wǎng)格技術(shù)和進出口壓力邊界條件的水輪機瞬態(tài)數(shù)值研究,能夠較可靠地模擬負荷變化過程中水輪機水頭、流量的變化規(guī)律。負荷變化過程中,無葉區(qū)壓力信號變化趨勢與導(dǎo)葉閉合規(guī)律一致。導(dǎo)葉閉合及停止運動的瞬間對尾水管壓力信號有較大的擾動作用。導(dǎo)葉停止運動進入部分負荷工況,尾水管內(nèi)逐漸產(chǎn)生低頻周期性壓力信號,且轉(zhuǎn)輪軸向力與尾水管內(nèi)壓力脈動關(guān)聯(lián)度最高,表明尾水管內(nèi)的壓力脈動向上游的傳播是影響軸向水推力形成的重要因素。

    (2)導(dǎo)葉閉合導(dǎo)致水輪機流量的降低,對尾水管內(nèi)的橫向流動及回流有一定的促進作用,加劇了尾水管內(nèi)的不穩(wěn)定流動。最優(yōu)工況下,水輪機內(nèi)部為光順流動,尾水管內(nèi)無自由渦結(jié)構(gòu)。負荷減小過程中,尾水管內(nèi)形成的強度較小的軸對稱渦帶結(jié)構(gòu),首先沿軸向拉伸,隨后沿軸向劇烈收縮但徑向尺寸增加。部分負荷工況下,軸對稱渦帶首先潰滅收縮,隨后演變?yōu)楦街谛顾F上的雙螺旋渦帶結(jié)構(gòu),最后雙螺旋渦帶合并為單一的強度較高的螺旋形尾水管渦帶。

    猜你喜歡
    導(dǎo)葉轉(zhuǎn)輪水輪機
    水輪機過流面非金屬材料的修復(fù)及防護
    大中型水斗式水輪機的關(guān)鍵技術(shù)
    詞語大轉(zhuǎn)輪
    ——“AABC”和“無X無X”式詞語
    寺廟里有座大書架——神奇的轉(zhuǎn)輪藏
    水輪機虛擬仿真動畫制作的研究
    壓氣機導(dǎo)葉用GTD222合金鑄造技術(shù)研究
    模具制造(2019年4期)2019-06-24 03:36:50
    我國第一臺分半鑄造的30萬千瓦水輪發(fā)電機轉(zhuǎn)輪實制成功
    蘭臺世界(2017年12期)2017-06-22 14:17:25
    水輪機過流部件改造與節(jié)能增效
    混流式模型水輪機導(dǎo)葉水力矩試驗
    一種新型的電子式導(dǎo)葉主令控制器的設(shè)計及應(yīng)用
    男人爽女人下面视频在线观看| 少妇被粗大的猛进出69影院| 亚洲在久久综合| 欧美日韩亚洲综合一区二区三区_| 欧美精品一区二区免费开放| 成人国语在线视频| 国产爽快片一区二区三区| 男人操女人黄网站| av.在线天堂| 免费高清在线观看日韩| 女人久久www免费人成看片| 亚洲成色77777| 老汉色av国产亚洲站长工具| 日本黄色日本黄色录像| 啦啦啦 在线观看视频| 成人免费观看视频高清| 亚洲第一区二区三区不卡| 满18在线观看网站| 亚洲人成电影观看| 国产亚洲最大av| 日韩免费高清中文字幕av| 丝袜在线中文字幕| 国产一区二区三区综合在线观看| 哪个播放器可以免费观看大片| 久久久久精品国产欧美久久久 | av卡一久久| 国产不卡av网站在线观看| 久久久精品区二区三区| 在线观看三级黄色| 免费观看a级毛片全部| 亚洲色图 男人天堂 中文字幕| 日韩精品有码人妻一区| 制服诱惑二区| 精品久久久久久电影网| 18禁裸乳无遮挡动漫免费视频| 极品人妻少妇av视频| √禁漫天堂资源中文www| 成人影院久久| 毛片一级片免费看久久久久| 下体分泌物呈黄色| 18禁动态无遮挡网站| 亚洲欧美一区二区三区国产| 精品一区二区三卡| 视频在线观看一区二区三区| 久久人人爽av亚洲精品天堂| 飞空精品影院首页| 大陆偷拍与自拍| 国产一区二区在线观看av| 日日摸夜夜添夜夜爱| 女的被弄到高潮叫床怎么办| 色94色欧美一区二区| bbb黄色大片| xxxhd国产人妻xxx| 亚洲国产精品国产精品| 大话2 男鬼变身卡| 熟女少妇亚洲综合色aaa.| 中文乱码字字幕精品一区二区三区| 亚洲欧美色中文字幕在线| 中文字幕另类日韩欧美亚洲嫩草| 高清av免费在线| 亚洲国产精品成人久久小说| 男女无遮挡免费网站观看| 男女下面插进去视频免费观看| 亚洲综合色网址| 不卡av一区二区三区| 国产一区有黄有色的免费视频| 欧美激情高清一区二区三区 | 麻豆乱淫一区二区| 日本黄色日本黄色录像| 在线观看免费午夜福利视频| 一级毛片我不卡| 91成人精品电影| 中文字幕人妻丝袜制服| 多毛熟女@视频| 国产精品一区二区精品视频观看| 欧美黑人欧美精品刺激| 精品一区二区三区四区五区乱码 | 亚洲欧美中文字幕日韩二区| 国产一卡二卡三卡精品 | 国产片特级美女逼逼视频| 巨乳人妻的诱惑在线观看| 亚洲欧美精品自产自拍| 国产精品三级大全| 侵犯人妻中文字幕一二三四区| 久久久国产欧美日韩av| 女人爽到高潮嗷嗷叫在线视频| 一本—道久久a久久精品蜜桃钙片| 天堂8中文在线网| 亚洲精品第二区| 丝袜喷水一区| 女人久久www免费人成看片| av在线观看视频网站免费| 高清av免费在线| 色视频在线一区二区三区| 老司机在亚洲福利影院| 欧美精品高潮呻吟av久久| 国产在线免费精品| 1024香蕉在线观看| 久久久国产欧美日韩av| 欧美97在线视频| 99热全是精品| 国产亚洲av高清不卡| 欧美 日韩 精品 国产| 飞空精品影院首页| 国产成人系列免费观看| 亚洲精品一区蜜桃| 久热这里只有精品99| 国产av一区二区精品久久| 精品人妻熟女毛片av久久网站| 日韩熟女老妇一区二区性免费视频| 国产男女超爽视频在线观看| 久久免费观看电影| 亚洲av男天堂| 国产男女内射视频| 一本大道久久a久久精品| 国产淫语在线视频| 波多野结衣一区麻豆| 国产成人免费无遮挡视频| 十八禁高潮呻吟视频| 色网站视频免费| 777米奇影视久久| 国产精品免费视频内射| 国产精品.久久久| 少妇人妻久久综合中文| 精品少妇一区二区三区视频日本电影 | 亚洲精品在线美女| 美女高潮到喷水免费观看| 女人精品久久久久毛片| 久久精品久久久久久久性| 亚洲精品美女久久av网站| 久久这里只有精品19| av不卡在线播放| 国产精品久久久人人做人人爽| 国产精品无大码| 一二三四中文在线观看免费高清| 国产淫语在线视频| 操出白浆在线播放| 最近手机中文字幕大全| 国产毛片在线视频| 精品一区二区免费观看| 韩国av在线不卡| 不卡视频在线观看欧美| 欧美日韩亚洲国产一区二区在线观看 | 女人爽到高潮嗷嗷叫在线视频| 精品国产一区二区三区四区第35| 1024视频免费在线观看| 人人妻人人添人人爽欧美一区卜| 91成人精品电影| 伊人亚洲综合成人网| 青青草视频在线视频观看| 成年av动漫网址| 久久人妻熟女aⅴ| 狂野欧美激情性xxxx| 免费黄色在线免费观看| 国产成人精品福利久久| 黑人猛操日本美女一级片| 黄色 视频免费看| 美女脱内裤让男人舔精品视频| 成人三级做爰电影| 久久久久久久国产电影| a级毛片在线看网站| av国产精品久久久久影院| 欧美日韩亚洲高清精品| 国产成人精品无人区| 亚洲情色 制服丝袜| 丝瓜视频免费看黄片| 天堂中文最新版在线下载| 黑人猛操日本美女一级片| 久久久久久久大尺度免费视频| 狠狠婷婷综合久久久久久88av| 国产亚洲午夜精品一区二区久久| 婷婷色av中文字幕| 欧美精品一区二区大全| 亚洲成人手机| 欧美黄色片欧美黄色片| 色综合欧美亚洲国产小说| 精品国产超薄肉色丝袜足j| 成人亚洲欧美一区二区av| 亚洲国产中文字幕在线视频| 97在线人人人人妻| 香蕉国产在线看| 丰满少妇做爰视频| 亚洲成av片中文字幕在线观看| av片东京热男人的天堂| 婷婷成人精品国产| 交换朋友夫妻互换小说| 咕卡用的链子| 欧美人与善性xxx| 免费不卡黄色视频| 97人妻天天添夜夜摸| 三上悠亚av全集在线观看| 亚洲成人免费av在线播放| 日本黄色日本黄色录像| 色综合欧美亚洲国产小说| 超碰97精品在线观看| 国产av码专区亚洲av| 久久影院123| 十八禁高潮呻吟视频| 午夜日本视频在线| 亚洲四区av| 91精品三级在线观看| 亚洲精品乱久久久久久| 久久人人97超碰香蕉20202| 高清黄色对白视频在线免费看| 久久综合国产亚洲精品| h视频一区二区三区| 中文字幕人妻熟女乱码| 侵犯人妻中文字幕一二三四区| 午夜福利免费观看在线| 亚洲情色 制服丝袜| 一个人免费看片子| 午夜老司机福利片| 秋霞在线观看毛片| 女性生殖器流出的白浆| 国产一区二区 视频在线| 午夜福利视频在线观看免费| 中文字幕人妻熟女乱码| 黄网站色视频无遮挡免费观看| 国产av码专区亚洲av| 久久人人爽av亚洲精品天堂| 秋霞伦理黄片| 色吧在线观看| 三上悠亚av全集在线观看| 欧美黑人精品巨大| 亚洲av成人精品一二三区| 日韩 亚洲 欧美在线| 久久久国产精品麻豆| 久久久精品94久久精品| 国产色婷婷99| 精品一区二区免费观看| 国产欧美日韩一区二区三区在线| 日韩av在线免费看完整版不卡| 丝袜美足系列| 国产一卡二卡三卡精品 | 飞空精品影院首页| 免费高清在线观看视频在线观看| 色婷婷久久久亚洲欧美| 大陆偷拍与自拍| 久久久久久久久免费视频了| 久久毛片免费看一区二区三区| 菩萨蛮人人尽说江南好唐韦庄| 国产精品99久久99久久久不卡 | 男人爽女人下面视频在线观看| 伦理电影大哥的女人| 国产精品三级大全| 欧美在线黄色| 99久久精品国产亚洲精品| 丰满饥渴人妻一区二区三| 亚洲欧美一区二区三区黑人| 亚洲国产最新在线播放| 欧美激情极品国产一区二区三区| 久久综合国产亚洲精品| 蜜桃在线观看..| 十分钟在线观看高清视频www| 欧美亚洲日本最大视频资源| 18在线观看网站| 日本欧美国产在线视频| 丝袜美腿诱惑在线| av卡一久久| 亚洲男人天堂网一区| 在线观看三级黄色| 午夜福利视频在线观看免费| 精品一区二区免费观看| 熟妇人妻不卡中文字幕| 哪个播放器可以免费观看大片| 中文字幕av电影在线播放| 精品亚洲乱码少妇综合久久| 午夜福利,免费看| 最黄视频免费看| av一本久久久久| 久久久久人妻精品一区果冻| 亚洲色图 男人天堂 中文字幕| 99久久综合免费| 女人高潮潮喷娇喘18禁视频| 国产成人精品久久久久久| 国产男女内射视频| 九色亚洲精品在线播放| www.av在线官网国产| 777米奇影视久久| 欧美 日韩 精品 国产| 国产伦人伦偷精品视频| 久久久久久人妻| 欧美最新免费一区二区三区| 久久久精品94久久精品| 久久人人97超碰香蕉20202| 香蕉丝袜av| 国产精品欧美亚洲77777| 国产精品久久久久久精品古装| 国产av精品麻豆| av网站在线播放免费| 国产精品久久久久久精品古装| 亚洲成人av在线免费| 亚洲免费av在线视频| 老司机在亚洲福利影院| 日韩视频在线欧美| 国产精品久久久av美女十八| 看免费成人av毛片| 国产一区二区在线观看av| 精品少妇久久久久久888优播| 纯流量卡能插随身wifi吗| 成年人免费黄色播放视频| 免费不卡黄色视频| 曰老女人黄片| 久久亚洲国产成人精品v| 免费在线观看黄色视频的| 亚洲自偷自拍图片 自拍| 国产成人午夜福利电影在线观看| 一区在线观看完整版| 女人爽到高潮嗷嗷叫在线视频| 久久久精品免费免费高清| 一区在线观看完整版| 午夜福利影视在线免费观看| 日韩一区二区视频免费看| 午夜av观看不卡| 午夜福利影视在线免费观看| 国产午夜精品一二区理论片| bbb黄色大片| 狂野欧美激情性bbbbbb| 日韩熟女老妇一区二区性免费视频| 国产亚洲欧美精品永久| 国产欧美日韩一区二区三区在线| 少妇被粗大猛烈的视频| 女的被弄到高潮叫床怎么办| 日韩伦理黄色片| 亚洲av综合色区一区| 亚洲精品久久久久久婷婷小说| 欧美日韩精品网址| 天天躁夜夜躁狠狠久久av| 欧美国产精品va在线观看不卡| 黑丝袜美女国产一区| 午夜福利视频精品| 看十八女毛片水多多多| 亚洲色图综合在线观看| 日本黄色日本黄色录像| 久久精品久久精品一区二区三区| 九草在线视频观看| 高清不卡的av网站| 国产精品.久久久| 啦啦啦在线免费观看视频4| 亚洲,欧美,日韩| 在线观看国产h片| 欧美变态另类bdsm刘玥| av一本久久久久| 国产av一区二区精品久久| 国产激情久久老熟女| 国产色婷婷99| av网站在线播放免费| 天天影视国产精品| 国产男人的电影天堂91| 精品一品国产午夜福利视频| 欧美国产精品一级二级三级| 1024香蕉在线观看| 综合色丁香网| 国产在线免费精品| 精品卡一卡二卡四卡免费| 色94色欧美一区二区| 日本爱情动作片www.在线观看| 免费在线观看完整版高清| 国产日韩一区二区三区精品不卡| 亚洲国产精品999| 国产日韩一区二区三区精品不卡| 尾随美女入室| 免费观看性生交大片5| 日本黄色日本黄色录像| 在线观看免费高清a一片| 伦理电影大哥的女人| 91精品国产国语对白视频| 国产一区二区激情短视频 | 国产 一区精品| 亚洲av综合色区一区| 亚洲美女黄色视频免费看| 天堂俺去俺来也www色官网| 国产极品天堂在线| 999精品在线视频| av免费观看日本| 国产成人午夜福利电影在线观看| 久久人人爽av亚洲精品天堂| 看免费成人av毛片| 9191精品国产免费久久| 男女之事视频高清在线观看 | 精品少妇一区二区三区视频日本电影 | 性少妇av在线| 亚洲欧美色中文字幕在线| 国产精品人妻久久久影院| 国产日韩欧美在线精品| 欧美另类一区| 亚洲,欧美精品.| 国产一区二区 视频在线| 中文字幕人妻丝袜一区二区 | 天堂中文最新版在线下载| 亚洲四区av| 午夜福利网站1000一区二区三区| 亚洲欧美日韩另类电影网站| 欧美日韩视频精品一区| av线在线观看网站| 少妇人妻久久综合中文| 国产色婷婷99| 久久久久人妻精品一区果冻| 亚洲国产日韩一区二区| 啦啦啦在线免费观看视频4| 成人免费观看视频高清| 又黄又粗又硬又大视频| 老鸭窝网址在线观看| 一二三四在线观看免费中文在| 又大又爽又粗| 18禁动态无遮挡网站| 国产精品秋霞免费鲁丝片| 丝瓜视频免费看黄片| 亚洲免费av在线视频| 国产色婷婷99| 在线观看国产h片| 久久精品国产亚洲av涩爱| 久久 成人 亚洲| 国产成人精品在线电影| 色婷婷久久久亚洲欧美| 18在线观看网站| 国产无遮挡羞羞视频在线观看| 观看av在线不卡| av线在线观看网站| 亚洲成国产人片在线观看| 久久狼人影院| 国产精品女同一区二区软件| 色94色欧美一区二区| 最近手机中文字幕大全| 午夜福利网站1000一区二区三区| 少妇 在线观看| 中文欧美无线码| 久久久久久久大尺度免费视频| 精品一区二区三区四区五区乱码 | 午夜免费男女啪啪视频观看| 两个人看的免费小视频| 成人免费观看视频高清| 欧美国产精品va在线观看不卡| 色精品久久人妻99蜜桃| 亚洲精品国产区一区二| 精品久久久久久电影网| 国产色婷婷99| 久久ye,这里只有精品| 国产一卡二卡三卡精品 | 久久亚洲国产成人精品v| 黄色怎么调成土黄色| 成人亚洲精品一区在线观看| 五月天丁香电影| 亚洲精品久久久久久婷婷小说| 亚洲国产欧美日韩在线播放| av片东京热男人的天堂| 精品人妻在线不人妻| 亚洲情色 制服丝袜| 国产精品.久久久| 老司机深夜福利视频在线观看 | 久久精品人人爽人人爽视色| 色婷婷久久久亚洲欧美| 久久99精品国语久久久| 在线亚洲精品国产二区图片欧美| 女人精品久久久久毛片| 亚洲人成网站在线观看播放| 久久97久久精品| 一级毛片黄色毛片免费观看视频| 又黄又粗又硬又大视频| 男的添女的下面高潮视频| 国产日韩欧美视频二区| 精品国产一区二区久久| 欧美在线黄色| 欧美日韩视频高清一区二区三区二| 人人妻人人澡人人爽人人夜夜| 免费日韩欧美在线观看| 婷婷色麻豆天堂久久| kizo精华| 国产男女内射视频| 日韩精品免费视频一区二区三区| 一级片'在线观看视频| 青春草视频在线免费观看| 国产一区二区 视频在线| 久久天堂一区二区三区四区| 久久久久精品国产欧美久久久 | 狂野欧美激情性bbbbbb| 欧美日韩福利视频一区二区| 亚洲欧美色中文字幕在线| 国产高清国产精品国产三级| 女人久久www免费人成看片| 不卡av一区二区三区| 欧美日韩成人在线一区二区| 伊人久久国产一区二区| 可以免费在线观看a视频的电影网站 | 亚洲美女视频黄频| 久久精品熟女亚洲av麻豆精品| 国产日韩一区二区三区精品不卡| 免费高清在线观看日韩| 久热爱精品视频在线9| 精品国产乱码久久久久久小说| 欧美精品高潮呻吟av久久| 久久影院123| 国产精品一区二区在线观看99| 中文乱码字字幕精品一区二区三区| 久久精品亚洲av国产电影网| 亚洲少妇的诱惑av| 国产人伦9x9x在线观看| 人人妻人人澡人人爽人人夜夜| 国产精品久久久久久精品电影小说| 欧美日韩视频精品一区| 一级爰片在线观看| 青青草视频在线视频观看| 人体艺术视频欧美日本| 一边摸一边做爽爽视频免费| 观看av在线不卡| 搡老岳熟女国产| 精品一区二区三区av网在线观看 | 久久 成人 亚洲| 男女下面插进去视频免费观看| 国产无遮挡羞羞视频在线观看| 最新在线观看一区二区三区 | 天天添夜夜摸| 秋霞在线观看毛片| 免费不卡黄色视频| 1024香蕉在线观看| 欧美乱码精品一区二区三区| 免费高清在线观看视频在线观看| av在线app专区| 久久久久精品国产欧美久久久 | 高清在线视频一区二区三区| 精品一区在线观看国产| 国产男女超爽视频在线观看| 久久人人爽人人片av| 亚洲av成人精品一二三区| 人妻 亚洲 视频| 91精品国产国语对白视频| 男的添女的下面高潮视频| 国产欧美日韩一区二区三区在线| 国产精品免费大片| 五月开心婷婷网| 婷婷色综合www| 亚洲av日韩在线播放| av女优亚洲男人天堂| 一级毛片黄色毛片免费观看视频| 国产精品秋霞免费鲁丝片| 欧美成人精品欧美一级黄| 国产麻豆69| 青春草视频在线免费观看| 日韩中文字幕视频在线看片| 国产淫语在线视频| 久久久久精品国产欧美久久久 | 操出白浆在线播放| 国产乱人偷精品视频| 爱豆传媒免费全集在线观看| 久久精品熟女亚洲av麻豆精品| 黑人欧美特级aaaaaa片| 久久精品人人爽人人爽视色| 超碰成人久久| 亚洲精品日韩在线中文字幕| 国产免费视频播放在线视频| 啦啦啦 在线观看视频| av免费观看日本| 电影成人av| 国产精品免费视频内射| 视频在线观看一区二区三区| 激情视频va一区二区三区| 十八禁人妻一区二区| 热re99久久精品国产66热6| 在线观看免费日韩欧美大片| 午夜精品国产一区二区电影| 18禁裸乳无遮挡动漫免费视频| 亚洲成人手机| 午夜福利网站1000一区二区三区| 啦啦啦啦在线视频资源| 国产成人91sexporn| 午夜免费观看性视频| 人人妻人人添人人爽欧美一区卜| 国产成人精品在线电影| a 毛片基地| 哪个播放器可以免费观看大片| 男女午夜视频在线观看| 日韩制服骚丝袜av| 日韩一区二区视频免费看| 亚洲色图综合在线观看| 精品免费久久久久久久清纯 | 一级毛片 在线播放| 亚洲三区欧美一区| 两性夫妻黄色片| 亚洲国产成人一精品久久久| 自拍欧美九色日韩亚洲蝌蚪91| 久热这里只有精品99| 久久人人爽av亚洲精品天堂| 亚洲成人av在线免费| 欧美日韩精品网址| 最黄视频免费看| 少妇精品久久久久久久| 国产精品一区二区在线不卡| 午夜免费男女啪啪视频观看| av卡一久久| 亚洲中文av在线| 精品久久久精品久久久| 80岁老熟妇乱子伦牲交| 婷婷色综合大香蕉| 欧美av亚洲av综合av国产av | 黄片小视频在线播放| 啦啦啦视频在线资源免费观看| 最近最新中文字幕大全免费视频 | 亚洲精品在线美女| 亚洲成国产人片在线观看| 精品国产一区二区三区四区第35| 青春草亚洲视频在线观看| 午夜免费男女啪啪视频观看| 在线观看免费午夜福利视频| 丰满饥渴人妻一区二区三| 七月丁香在线播放| 一区二区日韩欧美中文字幕| av卡一久久| 亚洲精品视频女| 悠悠久久av| 十分钟在线观看高清视频www| 亚洲第一青青草原| 日韩欧美精品免费久久| 午夜日本视频在线| 女人高潮潮喷娇喘18禁视频| 久久ye,这里只有精品|