999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

燃氣輪機發電機組余熱鍋爐低頻噪聲機理研究

2014-07-19 11:11:10唐虹云蔣偉康鐘振茂趙迎久
動力工程學報 2014年12期
關鍵詞:煙氣

唐虹云, 蔣偉康, 鐘振茂, 趙迎久

(1. 上海交通大學 振動、沖擊、噪聲研究所,上海 200240; 2. 華電重工股份有限公司,北京 100077)

燃氣輪機發電機組余熱鍋爐低頻噪聲機理研究

唐虹云1, 蔣偉康1, 鐘振茂2, 趙迎久2

(1. 上海交通大學 振動、沖擊、噪聲研究所,上海 200240; 2. 華電重工股份有限公司,北京 100077)

為研究燃氣輪機余熱鍋爐(HRSG)低頻噪聲的產生機理,采用大渦模擬(LES)結合FW-H積分方程法對爐內換熱器管陣的非定常流場和聲場進行數值模擬,研究了單個鰭片管、單個圓管以及圓管管陣的渦流噪聲特性,并與試驗結果進行了對比.結果表明:余熱鍋爐低頻噪聲主要源自換熱器管陣渦流噪聲;單根圓管的渦流噪聲具有明顯的低頻特征頻率,螺旋鰭片的存在降低了噪聲的特征頻率和總聲壓級;管陣中管間流場的相互作用增強,其渦流噪聲是具有一定帶寬的窄帶噪聲.

余熱鍋爐; 氣動噪聲; 大渦模擬; FW-H方程

隨著我國城市化進程的加快和生活用電需求的增長,以燃氣-蒸汽聯合循環為主的燃氣輪機電廠的建設地點越來越靠近居民區,在給人們的生活帶來便利的同時,燃氣輪機電廠的噪聲污染尤其是廠界低頻噪聲超標成為亟待解決的問題.燃氣輪機余熱鍋爐(HRSG)作為燃氣輪機電廠的主要設備之一,因其鍋爐本體聲輻射面積大且煙囪位置高而成為燃氣輪機電廠噪聲控制的重要對象.

燃氣輪機余熱鍋爐的爐內聲場十分復雜,國內外對該聲場產生機理的研究較少且結論并不一致.許志貴等[1]提出以燃氣輪機排氣噪聲為主的余熱鍋爐噪聲計算方法;Putnam[2]則通過案例分析,得出余熱鍋爐低頻噪聲由換熱器管陣渦流噪聲與聲腔模態耦合共振引起的結論.

筆者測量分析了某燃氣輪機電廠余熱鍋爐的廠界和機組各部件的輻射噪聲,發現影響廠界噪聲的主要成分是以63 Hz為中心的倍頻程帶.根據燃氣輪機排氣噪聲聲功率和鍋爐部件間傳聲損失,可以估算燃氣輪機排氣噪聲在鍋爐各部件處的聲壓級.根據實測的爐墻輻射噪聲聲壓級和爐墻的隔聲量可以獲得對應位置爐內噪聲的間接測量聲壓級,在以63 Hz為中心的倍頻程帶上,鍋爐內部噪聲間接測量聲壓級明顯高于燃氣輪機排氣噪聲的聲壓級.由此可知,低頻噪聲源不是燃氣輪機排氣噪聲,初步判斷為鍋爐內的換熱器管陣渦流噪聲.

管陣渦流噪聲問題可近似為圓柱繞流的氣動噪聲問題.國內外學者對不同雷諾數下的單圓柱、串列雙圓柱和錯排雙圓柱的漩渦脫落和輻射聲進行了大量的研究[3-8],結果表明漩渦脫落產生的脈動力是發聲的主要原因,漩渦脫落頻率與流噪聲的峰值頻率十分吻合.由于換熱器管通常帶有螺旋鰭片,并非理想的圓管,目前的研究主要針對螺旋鰭片對換熱性能的影響,關于螺旋鰭片對渦流噪聲影響的研究卻非常少.筆者采用大渦模擬(LES)結合FW-H積分方程的方法對爐內換熱器管陣的非定常流場和聲場進行了數值模擬,通過研究單個鰭片管、單個圓管和圓管管陣的漩渦脫落發聲機理及其噪聲特性,獲取渦流噪聲的特征頻率,分析換熱器管陣渦流噪聲的影響因素,并結合試驗結果,揭示余熱鍋爐低頻噪聲的主要聲源和產生機理.

1 換熱器管陣渦流噪聲模型的分析

1.1 工程模型簡化

由于鍋爐內存在流、熱、聲三者的耦合,且管陣本身是由帶螺旋翅片的鰭片管組成的,其幾何特征復雜,對工程對象進行完整建模計算需要大量的資源,故對分析模型進行如下簡化:

(1) 根據煙氣組分,采用氣體狀態方程計算煙氣的物理特性,其各項數據與空氣的對應數據相差不超過5%,因此計算中采用空氣代替煙氣;

(2) 假定鍋爐內的流場和溫度場是穩定的,每個熱交換器管陣內部的煙氣熱力學參數均勻不變;

(3) 將三維鰭片管管陣簡化為二維圓管管陣,建模仿真流程如圖1所示;

(4) 由于爐內不同部位管陣的設計參數和工質參數不同,因此選取最有代表性的位于最底層、溫度最高且流速最快的高壓過熱器和位于最頂層、溫度最低且流速最慢的除氧蒸發器為研究對象,計算其特征頻率,進而估算爐內整體的噪聲特征頻率范圍.

圖1 三維鰭片管管陣簡化為二維圓管管陣建模仿真流程圖

1.2 數值模擬

數值模擬采用計算流體力學(CFD)結合聲類比的計算方法[9-11].該算法核心是流體動力近場和聲學遠場的計算,前者通過求解湍流脈動為后者提供聲源信息,后者根據聲源信息計算遠場聲學量.筆者采用大渦模擬和FW-H積分方程外推法求解流場和聲場.

1.2.1 大渦模擬方法

大渦模擬方法運用濾波函數對N-S方程進行處理.

(1)

(2)

式中:上標“-”表示大尺度分量;u為流速;ui、uj分別為流體在xi、xj方向上的分量;p為壓強;ρ為流體密度;v為流體動力黏度.

(3)

有多種模型可以模擬亞格子尺度雷諾應力,此處采用Smargorinsky模型.

(4)

(5)

(6)

將小于網格尺度的漩渦從流場中過濾,只計算大渦,小渦的解則由亞格子模型模擬.此方法比直接數值模擬(DNS)效率更高,消耗系統資源更少,比雷諾平均(RANS) 方法更精確.

1.2.2 FW-H積分方程

FW-H積分方程由Lighthill方程推導得出,在氣動聲學領域具有重要意義.

(7)

式中:c0為聲速;ρ0和p0分別為聲學計算參考密度和壓強;v為垂直于積分面的流體速度;V為積分面移動速度;σ(f)為Dirac函數;τij為剪切應力;H(f)為Heaviside函數;Tij為Lighthill應力張量.

方程右端的3項分別對應單極子、偶極子和四極子聲源.該計算方法在本例中用于固定邊界、低馬赫數流動的計算條件,可忽略單極子和四極子聲源的作用,重點考慮偶極子聲源.

1.3 分析模型與邊界條件

對管徑d=40 mm的單根換熱器圓管建立10d×20d×2.5d的三維計算域和10d×20d的二維計算域,圓管距入口端、出口端和邊壁的長度分別為5d、15d和5d.三維模型中,入口采用速度邊界條件,出口采用自由流邊界條件,徑向兩側和軸向兩端均采用周期性邊界條件,圓管表面采用無滑移壁面邊界條件.鰭片管的有限元模型如圖2(a)所示,翅片寬度為1 mm,翅片高度為15 mm,管徑、管長、計算域尺寸和邊界條件與單根圓管的三維模型一致.根據廠方提供的換熱器管陣信息,取12×64原始管陣的局部陣列12×10作為分析對象(沿流動方向共12橫排,每排10根換熱器管),建立管徑d均為40 mm的管陣二維模型.管陣為等邊三角形錯排式,管間距為2.1d,計算域為25d×200d,流動方向和邊界條件與單根圓管二維模型一致,有限元模型如圖2(b)所示.為了準確捕捉邊界層流動信息,所有模型的近壁面第一層網格高度均為5×10-3d,保證壁面y+在1的量級.

1.4 分析工況

對鰭片管三維模型、圓管三維模型和圓管二維模型,采用高壓過熱器處和除氧蒸發器處2種煙氣的熱力學參數,并在5 m/s、10 m/s、15 m/s、20 m/s和25 m/s 5種流速下分別進行仿真計算.所選擇的仿真流速涵蓋了絕大多數種類余熱鍋爐的主要工況,其爐內煙氣流速遠低于常規電站鍋爐.

(a) 鰭片管有限元模型

(b) 管陣有限元模型

根據電廠實測煙氣流速進行仿真計算.煙氣參數見表1.瞬態計算時間步長為0.000 5 s,時間步數為4 000,噪聲分析頻率為1 000 Hz,頻率分辨率為0.5 Hz.

表1 煙氣熱力學參數

2 渦流噪聲機理與影響因素

2.1 鰭片管與圓管仿真結果的對比

圖3給出了高壓過熱器處煙氣參數下鰭片管和圓管模型的截面渦量云圖.由圖3可見,管壁均存在明顯的漩渦脫落現象.圖4給出了高壓過熱器處煙氣參數下三維圓管在不同流速下的聲壓級.由圖4可知,漩渦脫落產生的氣動噪聲存在明顯的低頻特征頻率.

(a) 鰭片管模型截面渦量分布

(b) 圓管三維模型截面渦量分布

(c) 圓管二維模型截面渦量分布

Fig.3 Vorticity distribution in middle section of different models at high pressure superheater area for flow velocity of 10 m/s

進一步比較鰭片管和圓管在不同煙氣參數和不同流速下的噪聲特征頻率和量級.圖5給出了2種煙氣參數下鰭片管和圓管的特征頻率隨流速的變化.由圖5可以看出,鰭片管和圓管的特征頻率均與流速呈良好的線性關系,煙氣參數的差別對其影響甚微;單一圓管的三維模型和二維模型在計算特征頻率時具有高度的一致性,為三維管陣簡化為二維管陣提供了依據;螺旋鰭片的存在使氣動噪聲的特征頻率移向低頻.

圖4 高壓過熱器處煙氣參數下圓管三維模型在不同流速下的聲壓級

Fig.4 Sound pressure level of the 3D tube model at high pressure superheater area for different flow velocities

圖5 不同煙氣參數下鰭片管和圓管的特征頻率隨流速的變化

Fig.5 Variation of characteristic frequency of finned-tube and single tube with flow velocity for different flue gas parameters

進一步分析螺旋鰭片的移頻效應.圓管特征頻率f的定義[12]為

(8)

式中:Sr為斯特勞哈爾數.已有研究結果表明,當雷諾數Re為1×103~1×105時,Sr可近似作為常數,故理論上可認為當管徑d一定時,特征頻率與流速呈線性關系,仿真煙氣的Re為2 500~12 500,滿足理論假設,且圓管和鰭片管的數值模擬結果也與這一結論吻合,故在余熱鍋爐正常工作的煙氣流速范圍內可將移頻效應看做螺旋鰭片帶來的等效直徑的影響.各流速下的等效直徑如圖6所示.由圖6可知,隨著流速的變化,等效直徑在1.27d左右變化,偏差在2%以內,具有較好的一致性,從工程應用的角度,可將各流速下圓管特征頻率按等效直徑折算得到對應鰭片管的特征頻率,為鰭片管管陣簡化為圓管管陣提供依據.不同規格的螺旋鰭片受幾何形狀和尺寸的影響而具有不同的等效直徑,筆者側重于提供等效直徑移頻這一工程應用方法,并不對其他影響因素進行詳細探討.

圖6 不同流速下鰭片管的等效直徑

圖7給出了2種煙氣參數下鰭片管和圓管的噪聲總聲壓級隨流速的變化.由圖7可知,噪聲總聲壓級與流速呈正相關,同時煙氣參數的差別對其有明顯影響,低溫氣體總聲壓級高出高溫氣體3~7 dB.由于二維模型在計算時存在聲相關長度(Sound Correlation Length)的選擇問題,因此單一圓管的三維模型和二維模型在計算噪聲聲壓級時存在較大差別,該參數決定了聲源的模擬尺度,進而決定了噪聲量級,但工程計算中無法預知該參數的精確值,故二維模型在計算噪聲量級時僅具備相對參考價值.螺旋鰭片的存在降低了氣動噪聲的總聲壓級,鰭片管與圓管總聲壓級相差10~12 dB.

圖7 不同煙氣參數下鰭片管和圓管的噪聲總聲壓級隨流速的變化

Fig.7 Variation of total sound pressure level of finned-tube and single tube with flow velocity for different flue gas parameters

2.2 管陣渦流噪聲仿真計算結果及分析

圖8給出了高壓過熱器處煙氣參數下管陣的截面渦量云圖.由圖8可知,管陣間存在明顯的漩渦脫落現象,但管間漩渦的產生和脫落規律與單管有區別:最接近入口的底層管子形成規律的尾渦,左右交替脫落;順著流向,越靠后的管子,其漩渦脫落越復雜,普遍存在前端管子漩渦脫落附著在后端管子上,經過湍流發展后再次脫落的現象.

圖8 高壓過熱器處煙氣參數下管陣模型的截面渦量分布

Fig.8 Vorticity distribution in middle section of tube array for high pressure superheater area

進一步研究管陣的噪聲特性.沿流向取出各橫排的換熱器管進行聲學計算.圖9以管陣第1排和第12排為例給出其噪聲頻譜.由圖9可知,各橫排換熱器管的噪聲均存在低頻特征頻率,且順著流向,橫向管排的特征頻率逐漸由95.9 Hz過渡到97.9 Hz,后排換熱器管在特征頻率附近(5 Hz的窄帶頻段)還存在多個次峰(與主峰峰值相差不到3 dB),即單一特征頻率鈍化為94.9~97.9 Hz的特征頻段,這一現象是由管間流場的相互作用引起的,周邊換熱器管漩渦脫落附著導致的脈動力波動,渦量云圖特征明顯.對比各橫排的噪聲量級可知,第1排雖然特征頻率明顯,但總聲壓級低于后排10 dB以上,這是因為后排管陣的湍流強度明顯增強,可見管陣間看似不規律的漩渦脫落在窄帶特征頻段處加強了氣動噪聲.

(a) 第1排換熱管的聲壓級

(b) 第12排換熱管的聲壓級

Fig.9 Sound pressure level of tube array at different horizontals for high pressure superheater area

除氧蒸發器管陣的計算結果也有類似規律,但由于該管陣處的氣體流速只有高壓過熱器管陣處的1/2,其噪聲量級相差15 dB以上,對爐內低頻噪聲貢獻量少,在此不再贅述.

2.3 數值模擬結果與試驗結果的對比

根據前文論述的螺旋鰭片的移頻效應,高壓過熱器圓管管陣的特征頻段94.9~97.9 Hz經移頻可得鰭片管管陣的特征頻段73.6~76.1 Hz,與試驗測得的噪聲頻譜特征頻段74~78 Hz高度吻合.

圖10給出了數值模擬得到的管陣渦流噪聲與試驗測試得到的爐內再生噪聲(即鍋爐內部噪聲間接測量值減去燃氣輪機排氣噪聲)的對比結果.由圖10可知,二者在倍頻程譜趨勢上高度吻合,均存在63 Hz的峰值特征頻段,可進一步證明爐內低頻噪聲主要由管陣渦流噪聲引起.因數值模擬的二維特性,仿真得到的管陣渦流噪聲的絕對量級與試驗結果不具備可比性.

圖10 渦流噪聲仿真結果與試驗結果的對比

Fig.10 Comparison of vortex noise between simulation results and experimental measurements

3 結 論

(1) 余熱鍋爐低頻噪聲主要源自換熱器管管陣渦流噪聲.煙氣流過密排管陣時,在鰭片管表面形成規律的漩渦脫落,漩渦脫落導致的脈動力以偶極子源的方式產生明顯的低頻噪聲.

(2) 單根圓管的噪聲具有明顯的特征頻率,螺旋鰭片的存在降低了噪聲的特征頻率和總聲壓級.

(3) 由于管間流場的干涉,管陣中單管渦流噪聲是具有單一頻率的簡諧聲波,但管陣渦流噪聲是具有一定帶寬的窄帶噪聲.

(4) 優化換熱管陣的設計,降低渦流噪聲可以減弱噪聲源;也可以在余熱鍋爐集中排煙段采取消聲降噪措施,以減少渦流噪聲從煙囪的排放.

[1] 許志貴,侯曉東. 燃機余熱鍋爐噪聲計算與控制[J].余熱鍋爐, 2007(1): 6-11.

XU Zhigui, HOU Xiaodong. Noise computation and control of the heat recovery steam generator[J]. Heat Recovery Steam Generator, 2007 (1): 6-11.

[2] PUTNAM R A. Low frequency noise from heat recovery steam generators[C]//Proceedings of Inter-noise 2006. Honolulu, Hawaii, USA:I-INCE, 2006.

[3] COX J S, RUMSEY C L, BRENTNER K S,etal. Computation of vortex shedding and radiated sound for a circular cylinder: subcritical to transcritical Reynolds numbers[J]. Theoretical and Computational Fluid Dynamics, 1998, 12: 233-253.

[4] 麻偉巍,蔣運幸,余宇軒. 雙圓柱繞流的數值模擬[J]. 復旦學報:自然科學版,2013,52(5):592-600.

MA Weiwei, JIANG Yunxing, YU Yuxuan. The numerical simulation of flow around two circular cylinders[J]. Journal of Fudan University: Natural Science,2013,52(5):592-600.

[5] 劉景偉,郭海燕,趙婧. 等直徑串列雙圓柱體繞流的數值模擬[J]. 中國海洋大學學報,2013,43(12):92-97.

LIU Jingwei, GUO Haiyan, ZHAO Jing. Numerical simulation of flow around two tandem circular cylinders of equal diameter[J]. Periodical of Ocean University of China, 2013,43(12):92-97.

[6] 王云飛,陳浮,李國占,等. 不同間距下并列圓柱繞流的大渦模擬[J]. 工程熱物理學報,2013,34(11):2039-2042.

WANG Yunfei, CHEN Fu, LI Guozhan,etal. Large eddy simulation of flow around side-by-side circular cylinders at various spacing[J]. Journal of Engineering Thermophysics, 2013,34(11):2039-2042.

[7] 唐科范,FRANKE J. 用解NS方程和FW-H積分的混合方法計算圓柱繞流噪聲[J]. 水動力學研究與進展, 2009, 24(2): 190-199.

TANG Kefan, FRANKE J. Numerical simulation of noise induced by flow around cylinder using the hybrid method with the solutions of NS equation and FW-H integration[J]. Chinese Journal of Hydrodynamics, 2009, 24(2): 190-199.

[8] 龍雙麗,聶宏,許鑫. 不同雷諾數下圓柱繞流氣動噪聲數值模擬[J]. 聲學技術,2011,30(2): 111-116.

LONG Shuangli, NIE Hong, XU Xin. Numerical simulation of noise induced by flow around a cylinder at different Reynolds number[J]. Technical Acoustics, 2011, 30(2): 111-116.

[9] 徐俊偉,吳亞鋒,陳耿. 氣動噪聲數值計算方法的比較與應用[J]. 噪聲與振動控制,2012(4):6-10.

XU Junwei, WU Yafeng, CHEN Geng. Comparison and application on the aero-acoustics numerical computing methods[J]. Noise and Vibration Control, 2012(4): 6-10.

[10] 余雷,宋文萍,韓忠華,等. 基于混合RANS/LES方法與FW-H方程的氣動聲學計算研究[J]. 航空學報,2013,34(8):1795-1805.

YU Lei, SONG Wenping, HAN Zhonghua,etal. Aeroacoustic noise prediction using hybrid RANS/LES method and FW-H equation[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(8): 1795-1805.

[11] 魏志,王玉璋. 閥體后90°圓形彎管內部流場的數值模擬[J]. 動力工程學報,2013,33(8):595-599.

WEI Zhi, WANG Yuzhang. Numerical simulation on inner flow in a 90° circular-sectioned bend with fore-end valve[J]. Journal of Chinese Society of Power Engineering, 2013,33(8):595-599.

[12] BLEVINS R D. Flow-iduced vibration[M]. 2nd ed. New York: Van Nostrand Reinhold Company Inc., 1990.

Research on Formation Mechanism of Low Frequency Noise in a Heat Recovery Steam Generator

TANGHongyun1,JIANGWeikang1,ZHONGZhenmao2,ZHAOYingjiu2

(1. Institute of Vibration, Shock and Noise, Shanghai Jiaotong University, Shanghai 200240, China; 2. Huadian Heavy Industries Co., Ltd., Beijing 100077, China)

To investigate the formation mechanism of low frequency noise in a heat recovery steam generator (HRSG), numerical simulation was carried out to analyze the non-stationary flow field and aero-acoustic noise in the heat exchanger tube array using large eddy simulation (LES) method combined with Ffowcs Williams-Hawkings (FW-H) equation. Features of vortex noise were studied respectively for sing tube, finned-tube and tube array, and the results were compared with that of actual measurements. Results show that the low frequency noise of HRSG mainly comes from the vortex shedding in the heat exchanger tube array. The vortex noise of single tube has obvious characteristic frequency, and the spiral fin reduces both the characteristic frequency and total sound level of finned tube. For tube arrays, the vortex noise is of the narrowband type with a certain bandwidth, due to enhanced interaction of flow fields between the tubes.

heat recovery steam generator; aerodynamic noise; large eddy simulation; FW-H equation

1674-7607(2014)12-0978-07

TK229,TB533

A

470.30

2014-03-03

2014-04-02

唐虹云(1989-),男,重慶人,碩士研究生,研究方向為:流體噪聲. 蔣偉康(通信作者),男,教授,博士,電話(Tel.):021-34206332-820;E-mail:wkjiang@sjtu.edu.cn.

猜你喜歡
煙氣
固體吸附劑脫除煙氣中SOx/NOx的研究進展
化工管理(2022年13期)2022-12-02 09:21:52
水泥窯高溫高塵SCR煙氣脫硝系統的應用
SCR煙氣脫硝噴氨優化改造研究
煙氣爐在干熄焦溫風干燥中的應用
山東冶金(2019年2期)2019-05-11 09:12:16
基于參數自整定模糊PID的SCR煙氣脫硝控制
測控技術(2018年2期)2018-12-09 09:00:52
基于煙氣脫硫儀控系統的設計與應用
電子測試(2018年1期)2018-04-18 11:52:15
燒結煙氣脫硫工藝技術分析
燒結煙氣中二氧化硫的測試方法
基于非分散紫外吸收法的便攜式煙氣分析儀在煙氣二氧化硫監測中的應用研究
煙氣導向式雙室熔煉爐
主站蜘蛛池模板: 国产自无码视频在线观看| 国产成人久久综合777777麻豆 | 国产黑丝视频在线观看| 国产打屁股免费区网站| av色爱 天堂网| 亚洲福利片无码最新在线播放| 亚洲视频四区| 国产精品亚洲一区二区三区z| 国产真实乱子伦视频播放| 久久国产精品娇妻素人| 国产成人免费高清AⅤ| 国产第一页亚洲| 九九九精品成人免费视频7| 国产一区免费在线观看| 嫩草国产在线| 国产精品太粉嫩高中在线观看| 亚洲欧州色色免费AV| 71pao成人国产永久免费视频| 国产在线八区| 亚洲无线一二三四区男男| 91欧洲国产日韩在线人成| 国产欧美日本在线观看| 亚洲成人高清在线观看| 夜夜高潮夜夜爽国产伦精品| 99久久这里只精品麻豆| 美女毛片在线| 中文字幕在线永久在线视频2020| 欧美一区二区三区香蕉视| 成人日韩欧美| 国产精品美女网站| 伊人久久青草青青综合| 欧美一区二区三区国产精品| 伊人久久影视| 五月天综合婷婷| 国产精品成人一区二区不卡| 国产国产人在线成免费视频狼人色| 亚洲无码A视频在线| 国产91透明丝袜美腿在线| 日韩第一页在线| 亚洲制服中文字幕一区二区| 欧美在线综合视频| 欧美激情视频在线观看一区| 色噜噜狠狠色综合网图区| 全部免费特黄特色大片视频| 欧洲熟妇精品视频| 91日本在线观看亚洲精品| 日韩无码视频网站| 九九九久久国产精品| 日韩美毛片| 欧美视频在线第一页| 国产在线观看成人91| 四虎永久免费在线| 国产一区二区视频在线| 国产本道久久一区二区三区| 露脸真实国语乱在线观看| 88国产经典欧美一区二区三区| 波多野结衣第一页| 日韩黄色大片免费看| 麻豆精品视频在线原创| 亚洲av无码人妻| 欧美中文字幕一区| 九九热精品在线视频| 亚洲日本www| 日韩国产无码一区| 日本黄色不卡视频| 亚洲成人高清无码| 99久久免费精品特色大片| 国产一区二区三区精品久久呦| 91免费精品国偷自产在线在线| 亚洲中文字幕23页在线| av一区二区无码在线| 国产乱论视频| 欧美有码在线| 夜夜爽免费视频| 午夜丁香婷婷| 老熟妇喷水一区二区三区| 曰AV在线无码| 国产午夜精品鲁丝片| 欧美精品一二三区| 成人综合网址| 波多野结衣一区二区三区AV| AV不卡在线永久免费观看|