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

高超聲速風洞來流擾動測量及數據后處理技術研究

2019-11-07 10:52:34王俊鵬劉向宏趙家權
實驗流體力學 2019年5期
關鍵詞:模態測量

余 濤, 王俊鵬, 劉向宏, 趙家權, 吳 杰,2,*

(1. 華中科技大學 航空航天學院, 武漢 430074; 2. 不倫瑞克工業大學 流體力學所, 德國 不倫瑞克 38108)

0 引 言

風洞實驗是研究高超聲速空氣動力學的重要手段之一,因此高超聲速風洞是研究高超聲速空氣動力學不可或缺的地面實驗平臺。常規高超聲速風洞建成后,必須對流場的靜態和動態品質進行鑒定。其中,動態流場品質,即自由來流擾動的模態特征,是衡量風洞設計成功與否的重要指標。對于常規低速風洞,衡量其動態品質的參數是湍流度,即速度偏差的均方根[1]。相對而言,高超聲速風洞的流場結構復雜,其流場的動態品質定義更加困難。此外,諸多高超聲速風洞實驗發現:流場的動態品質對風洞實驗測量結果有嚴重影響,如層/湍流邊界層轉捩[2-5]、激波/邊界層干擾模式[6]、阻力系數測量[7]等,但是其中的作用機理仍不確定,有待于探索。尤其是高超聲速流動下的邊界層轉捩問題,不同類型來流擾動子模態直接決定了轉捩點位置及轉捩波的類型,如渦波模態通常引起第一模態(Tollmien-Schlichting wave)轉捩,而聲波模態則導致第二模態(Mack wave)轉捩[2, 8-10]。隨著我國高超聲速飛行器設計技術的迅速發展,研究人員對高超聲速風洞實驗的精度要求愈趨嚴格。因此,準確測量高超聲速風洞中自由來流擾動子模態對基于高超聲速風洞開展的相關基礎科學問題的研究至關重要。

國外對高超聲速風洞流場擾動的研究起步較早。1953年,Kovasznay使用熱線儀對超聲速來流擾動進行測量并建立了超聲速流動的擾動模態理論[11-12]。其主要思想是將一階小擾動引入到粘性、可壓縮、帶熱傳遞效應的Navier-Stokes方程,將方程線性化后對應的解定義為超聲速流動的擾動子模態,即渦波、熵波和聲波模態。其中,渦波是速度場中的旋轉分量的擾動;熵波則是在恒壓下熵、密度和溫度的變化;聲波模態主要描述在熵守恒情況下壓力、密度和溫度的擾動量。在這3種模態中,渦波和熵波滿足流線方程,因此可以通過流線追蹤到擾動的產生源;相比之下,聲波的產生機理更加復雜,且其可以跨越流線發展。此后,Morkovin對超聲速自由來流中的擾動來源展開了更為深入的研究[13-14]。在Kovasznay對來流擾動分類的基礎上,Morkovin進一步將聲波模態細化為渦流馬赫波和顫振馬赫波(Eddy Mach wave and Shivering Mach wave)。前者主要由超聲速湍流邊界層所產生的以偶極子和四極子為主的聲波輻射組成;后者則是由噴管上超聲速湍流邊界層的壓力間斷所引起的反射和衍射現象。此外,Morkovin和Laufer分別通過實驗研究發現,在超聲速風洞的自由來流中聲波模態是主要擾動源[14-15]。普渡大學Schneider教授提供的一張圖片形象地描述了高超聲速風洞中自由來流的擾動模態類型及演化特征,見圖1[16]。

圖1 高超聲速風洞自由來流的擾動[16]

早期超聲速及高超聲速風洞的來流擾動測量以熱線風速儀為主。由于熱絲的直徑小(本文所使用Dantec 55P11熱絲直徑為5 μm),在超聲速及高超聲速流場測量中通常可以近似認為熱絲測量的流場脈動可以真實反映自由來流的脈動[17]。基于該假設,熱線風速儀在超聲速及高超聲速流場測量中得到了廣泛應用,如Smits使用熱線風速儀對馬赫數2.9的超聲速流進行測量[18]、Wu等人對馬赫數3的超聲速邊界層的均值和脈動剖面進行測量并與線化理論吻合[19]、Weiss等人對馬赫數2.5的自由來流進行流場擾動校核[20]等。通過熱線儀測量可得到流場的流量和總溫脈動,以及兩者的相關系數[11,18]。但是,流量和總溫脈動并不是典型的非定常流動特征量,也不足以推導Kovasznay所定義的3種擾動子模態[21]。Laufer通過實驗觀察發現,在超聲速風洞中噴管的湍流邊界層所產生的聲波模態通常為占優模態,并且通過熱線測量發現流量和總溫的擾動之間高度互相關[15, 20]。因此,Laufer假設來流擾動為單一噪聲場,各流場參數之間滿足等熵條件,進而可以推導其他流場參數的擾動特征。但是,Kovasznay認為在高超聲速風洞自由來流中等熵假設不具備普適性,建議采用3種擾動子模態互不相關假設[11]。隨著靜風洞的出現[16,22],風洞來流擾動主要來源于上游的渦波與熵波,Laufer提倡的聲波模態占優假設重新得到應用[23-25]。但是對于常規高超聲速噪聲風洞,分解來流擾動模態則需另辟蹊徑。Logan等提出了一種新型的擾動模態分解法[26, 27],假設流場參數的脈動量(速度、溫度和密度脈動)僅有兩個不同的來源:聲波模態

和其他模態,并認為聲波模態與其他擾動模態不相關。該方法整體低估了聲波模態對擾動模態的貢獻。Masutti對馮卡門研究所的H3馬赫數6風洞進行來流擾動模態標定,采用了雙熱線探頭與Pitot管聯合測量方法。在基于流場參數互不相關的假設下,發現聲波模態并非主要擾動,尤其在高雷諾數下熵波甚至高于聲波擾動模態的幅值[28]。Wu等人則通過單熱線探頭與Pitot管聯合測量法,對德國不倫瑞克工業大學的馬赫數3超聲速風洞進行擾動模態測量,發現聲波模態為來流的主要擾動模態,但是渦波與熵波模態也不可忽略[29-30]。之后,Schilden等人采用錐體探頭-DNS聯合法對德國不倫瑞克工業大學的馬赫數6高超聲速風洞進行了來流擾動模態的測量[31],在忽略渦波模態的前提下,結合實驗測量與直接數值模擬方法,獲得了來流中熵波與聲波模態的幅值。與Wu[29]之前所使用的測量法得到的結果進行對比,發現聲波模態在低頻范圍內吻合較好。Chaudhry 和 Candler采用DNS研究了自由來流擾動中聲波擾動經過皮托探頭產生的脫體激波后的傳遞函數[32],結果顯示自由來流中的壓力脈動經過脫體激波后脈動幅值反而增加,與Harvey以及Stainback[33]分別基于定常與非定常過程分析的結果有較大差距。

國內傳統的高超聲速風洞流場校測是將氣流的速度場(馬赫數分布的均勻度、梯度)和方向場(平均氣流偏角)作為標準,而忽略了流場的動態特征,如壓力、溫度和密度脈動等,相關方面的理論研究較為缺乏。直到近年來,國內學者開展高超聲速邊界層轉捩以及不穩定性研究,意識到來流擾動問題對邊界層轉捩波的重要性,方才重視風洞的來流擾動問題,如北京大學靜風洞使用Kulite壓力傳感器測量皮托脈動壓力[34]。雖然以上實驗結果具有重大意義,但是仍未提供自由來流的擾動模態信息,如渦波和熵波擾動。

整體而言,目前高超聲速風洞來流擾動的模態測量與后處理問題仍存在較大難度,國內外對于高超聲速風洞來流擾動的標定尚無定論。在實際情況中,不同的風洞由于設計、加工方案的不同,流場的擾動子模態所扮演的角色也不同,加上部分測量儀器本身的激波干擾問題,更進一步將擾動模態分析問題復雜化。本文基于Kovasznay的擾動模態分解理論,討論了不同類型高超聲速風洞的來流擾動子模態分析方法,并以德國不倫瑞克工業大學馬赫數6 Ludwieg式高超聲速風洞的實驗數據為例,展示了一種適用于常規高超聲速風洞的來流擾動子模態分解技術。

1 實驗設備、儀器及測量方法

1.1 馬赫數6 Ludwieg管風洞

本文所研究的高超聲速來流擾動模態測量與分析基于德國不倫瑞克工業大學串列式馬赫數6高超聲速風洞完成(實驗部分由本文通信作者在德國不倫瑞克工業大學完成)。該高超聲速風洞是一座Ludwieg式風洞,其實驗段口徑為500 mm,風洞的有效運行時間為80 ms。圖2為馬赫數6 Ludwieg式風洞示意圖[35]。風洞的高壓儲氣段與Laval噴管通過快速控制閥分開。在風洞啟動前,儲氣段內儲存著高溫高壓空氣,而控制閥下游的部分則通過真空泵抽成真空。在開啟快速控制閥門的瞬間,會產生一系列的非定常膨脹波,該膨脹波以聲速向儲氣段的上游行進;與此同時,該膨脹波驅動管內的氣體,達到對應的儲氣段啟動馬赫數(該馬赫數取決于儲氣段與Laval喉道部分的面積比)。當膨脹波到達儲氣段底端后,再次被反射回來。反射膨脹波到達快速控制閥附近時,快速控制閥關閉,風洞的運行也同時結束。在快速控制閥的下游,受膨脹波驅動的高壓氣體會沿著Laval噴管膨脹加速,氣流在到達實驗段時加速到設計馬赫數。更多關于該風洞的信息可參考文獻[35-41]。

圖2 德國不倫瑞克工業大學馬赫數6 Ludwieg式高超聲速風洞[35]

Fig.2Mach6LudwiegwindtunnelatTechnicalUniversityofBraunschweiginGermany[35]

1.2 實驗測量儀器及數據處理方法

在高超聲速風洞中,來流擾動脈動頻率高、流動環境復雜,對來流擾動測量技術有較高要求。目前,廣泛應用于高超聲速風洞自由來流擾動測量的儀器主要是熱線風速儀、皮托探頭以及聚焦激光差分干涉儀。其中,熱線風速儀可以獲取來流的總溫和流量脈動,其動態響應頻率可達百千赫茲以上;壓力皮托管測量的是激波后總壓脈動,雖然動態響應頻率更高,但是其所測得的信號并非自由來流中的真實壓力脈動;聚焦激光差分干涉儀采用的是非介入式測量方式,可以測量自由來流中的密度脈動。本文使用了恒溫熱線風速儀(Constant Temperature Anemometer, CTA)和皮托管,見圖3和4,結合兩種不同儀器的測量數據,通過擾動模態分析,獲取不同擾動子模態的特征。

圖3 熱線風速儀系統

圖4 皮托壓力探頭

熱線儀的測量原理是熱絲產生的熱量與周圍氣流對流換熱平衡。對于超聲速及高超聲速流動,目前廣泛采用以下經驗公式描述該過程:

(1)

其中,Nu為Nusselt數,I為通過熱絲的電流,R為熱絲的電阻,k為熱絲周圍流體熱傳遞系數,L為熱絲長度,T為溫度,η是總溫恢復系數。下角標w表示熱絲的狀態,0則表示駐點流動狀態。對于熱絲而言,X和Y通常為常數,其主要取決于熱絲本身的特性。f(τ)和g(τ)分別是過熱比τ的函數;n為雷諾數指數,通常在0.40~0.55之間。根據經驗,當來流馬赫數大于1.3時,Nu與馬赫數無關;此外,當熱絲雷諾數大于20時,總溫恢復系數η為常數(0.94±0.001)[42]。由于熱線儀的電壓輸出是流量ρu和總溫T0的函數,可表達為:

E=E(ρu,T0)

(2)

對式(2)進行如下轉換:

E′=km(ρu)′+kT0(T0)′

(3)

其中,

(4)

(5)

km為流量敏感系數,kT0為總溫敏感系數,這兩個重要的參數都是熱絲過熱比的函數,需要通過校核進行標定。在實際數據的分析過程中,需要對式(3)進行變換,通常的做法是將式(3)進行平方后整理為均方根值的形式:

(6)

在使用熱線儀進行測量前,通常需要對熱線儀進行校核,確定熱線的敏感系數。在低速流動中,通常將流動視為不可壓縮流動,且對環境溫度的考慮不多。對于超聲速及高超聲速自由來流而言,可壓縮性與溫度的影響不可忽略。為了校核方便,通常將式(1)表達為以下形式:

E2=F+M(ρu)n

(7)

其中,F和M是過熱比和總溫的函數。對式(2)進行變化得到:

(8)

其中,Ra是熱絲在環境溫度下的電阻。對于式(7)和(8),可以將過熱比函數F和M與熱線儀參數以及周圍流場參數聯系起來:

(9)

(10)

在式(9)和(10)中,D為熱絲的直徑,μ為流體的動力粘性系數。進一步可以將熱線儀的流量敏感系數和總溫敏感系數進行如下表達:

(11)

(12)

其中,a、b為常數,n為式(1)中的雷諾指數。認為熱絲周圍的流體熱傳遞系數以及動力粘性系數滿足如下關系式:

(13)

(14)

其中,下角標ref表示參考值,對于理想氣體,ak=bμ=0.768[11]。為了獲取熱線儀的流量和總溫敏感系數,通常的做法是在某一特定總溫下,改變來流的雷諾數,對流量的敏感系數進行校核;對于溫度的敏感系數,則需要在特定流量下通過改變來流的總溫來獲取。對于常規超聲速及高超聲速風洞,改變來流的總溫后,往往涉及到流量的變化。因此,本文選擇了Smits發展的方法[18],即保持風洞的總溫不變,在不同的過熱比狀態下改變風洞自由來流的單位雷諾數,建立起熱線風速儀輸出電壓與流場自由來流流量的關系式,如圖5所示。基于圖5中的校核曲線,可以同時獲取f(τ)和g(τ)的函數分布。在確定了f(τ)和g(τ)的函數分布后,根據式(11)和(12)可最終獲取流量和總溫敏感系數的分布,如圖6所示。由圖6可知,在低過熱比條件下,熱線儀對總溫變化更為敏感;隨著過熱比的增加,熱線儀對流量變化更為敏感。

圖5 熱線儀校核曲線

圖6 熱線儀流量與總溫敏感系數比較

皮托探頭測量激波后的總壓脈動,為了與熱線測量的數據進行耦合,需要將總壓脈動轉化為自由來流的靜壓脈動。由于不同的擾動子模態經過激波后均產生新的聲波脈動,即使通過直接數值模擬,也只能對極為簡化的情況進行分析。在本實驗中,采用Chaudhry等人最近直接數值模擬所獲得的傳遞函數[32]。

2 高超聲速來流擾動模態及其分解方法

通過熱線風速儀與皮托探頭測量,可獲取流量、總溫以及壓力脈動的信息,但是以上信息并未與擾動子模態進行關聯。為了定量標定擾動子模態,需要進行擾動模態離解分析。

Kovasznay在1950年首先發展了使用熱線儀測量超聲速來流中擾動子模態的技術,并且基于小擾動理論發展了擾動模態理論,對超聲速流動中的小擾動進行定性與定量的描述[11-12]。根據不同流場參數的特征,Kovasznay認為超聲速流場中的擾動是3種擾動子模態的疊加,即渦波ω、熵波θ和聲波σ擾動,各子模態擾動的定義如下所示:

(15)

(16)

(17)

上角標′表示流場參數的脈動,γ為比熱比,Ma為來流馬赫數,u、T、p分別為速度、溫度和壓力,nx為聲波的方向余弦。通過以上定義可見,來流擾動模態通常由多個流場參數以及它們之間的相關系數決定。首先,根據本實驗測量儀器所能獲取的流場參數的定義,可得:

m=ρu

(18)

(19)

p=ρRT

(20)

T為靜溫,cp為定壓比熱容。由于已經假設自由來流方向與熱線放置方向垂直,可認為除了流向以外,其他方向的速度脈動為零。對式(18)~(20)進行變化可得:

(21)

(22)

(23)

其中,

(24)

β=α(γ-1)Ma2

(25)

將式(21)~(23)左右變量進行更換后,代入到式(15)~(17),即可得到擾動子模態與實驗測量參數之間的關系式:

(26)

(27)

(28)

根據式(26)~(28)可知,在小擾動假設條件下,聯合熱線風速儀和皮托探頭測量,可最終獲得自由來流的擾動子模態幅值。

3 實驗數據與分析

3.1 皮托管測量

本文皮托壓力測量過程中采用的是Kulite XCQ 系列壓力傳感器,其有效頻率最高達40 kHz。實驗中最低來流單位雷諾數Reunit為8.35×106/m,最大為15.6×106/m。圖7為馬赫數6自由來流中不同單位雷諾數下的皮托壓力脈動分布,實驗結果顯示在低雷諾數下皮托壓力ppitot脈動幅值約為1.8%。整體而言,該風洞的流場質量較馮卡門研究所馬赫數6風洞和普渡大學馬赫數6風洞的噪聲狀態低[28,35]。此外,隨著單位雷諾數的增加,皮托壓力脈動幅值下降,說明該高超聲速風洞屬于典型噪聲風洞[43],其噪聲主要來源為Laval噴管表面湍流邊界層所產生的噪聲輻射。

圖7 自由來流中皮托管壓力脈動

但是,皮托管測量的是激波后的總壓脈動,不能直接視為自由來流的壓力脈動,我們這里使用Chaudhry等直接數值模擬所獲得的傳遞函數[32]。 Chaudhry等使用直接數值模擬對具有迎角的聲波擾動穿越激波過程進行分析,探明了自由來流中壓力脈動穿越激波的變化特征,并建立了來流擾動與激波后總壓脈動之間的關系(圖8),可見皮托探頭所產生的激波前后壓力脈動比值同時也是頻率的函數。 由于本文使用的壓力傳感器的動態響應僅達到40 kHz,圖8顯示在40 kHz以下,激波前后壓力脈動的比值約為0.38~0.4(之后所使用的熱線儀數據也僅僅使用45 kHz以下信號)。根據以上關系,可以估算皮托總壓脈動與自由來流靜壓脈動的均方根比值為0.632。

圖8 皮托總壓與馬赫數6自由來流壓力脈動傳遞函數(聲波迎角為120°)[44]

Fig.8TransferfunctionbetweenpitottotalpressureandMach6freestreampressuredisturbance(sonicattackangleis120°)[44]

3.2 熱線測量

為了確保實驗數據具有統計意義,熱線儀的測量過程中熱絲探頭的位置與皮托探頭的位置保持一致,風洞運行的雷諾數也保持相同。受限于熱絲的最大溫度,實驗過程中使用的最大過熱比僅達到0.55。通過熱絲校核,得到流量和總溫脈動的均方根值,如圖9所示。流量的脈動隨著來流雷諾數的增加而逐漸降低,與皮托總壓脈動的趨勢相似;但是,總溫的脈動在低雷諾數下跳躍較大,且隨著雷諾數增加,其幅值也增加。

圖9 馬赫數6自由來流流量和總溫脈動

根據第2節擾動模態分解的方法,進行各擾動子模態的計算,獲得圖10所示結果。聲波擾動模態遠高于渦波和熵波擾動,占據了總擾動的60%以上;相比之下,熵波和渦波擾動較低,約各占總體擾動量的15%。由于該Ludwieg式管風洞流場結構極為簡單,沒有復雜的閥門機構,渦波擾動基本可以忽略;至于熵波擾動,本次實驗過程中僅僅加熱至430 K,且嚴格控制風洞運行間歇的穩流時間,故其幅值預期也較低。

圖10 自由來流擾動子模態分析

此外,為了進一步證實該風洞屬于噪聲為主要擾動的風洞,對熱線輸出電壓信號與敏感系數之間的關系進行了線性擬合,如圖11所示。擬合結果顯示,獲得的熱線電壓輸出與熱絲敏感系數存在良好的線性吻合,且隨著來流雷諾數的增加,擬合直線的斜率減小,與Laufer判定聲波占優脈動風洞的特征一致[15]。由于聲波占優假設不是本文的重點,這里不再闡述,更多細節可參考文獻[15]。

圖11 熱線儀輸出電壓與敏感系數線性擬合

Fig.11Linearfittingofoutputvoltageandsensitivitycoefficientofhot-wire

4 總結與討論

本文基于德國不倫瑞克工業大學馬赫數6 Ludwieg式高超聲速風洞的自由來流試驗數據,介紹了高超聲速風洞中自由來流擾動模態的測量與分析問題,通過聯合使用熱線風速儀和皮托壓力探頭測量技術,并結合直接數值模擬結果,進行來流擾動模態的離解分析,最終獲取了高超聲速流動下不同擾動子模態的幅值。實驗結果顯示,該風洞為典型的噪聲風洞,其聲波模態高達擾動總模態的69%,熵波和渦波模態各約占15%。該實驗技術為國內諸多高超聲速風洞的流場擾動測量提供了思路,為基于高超聲速風洞開展的實驗提供了借鑒與參考。

但是,本文在后處理方法上假設流場脈動量之間互不相關,在小擾動假設前提下成立,有一定的局限性。此外,該實驗數據后處理方法與Ali[45]及Wu等[29,39]所使用的方法不同,且在子模態擾動幅值上呈現較大差距,其原因還需要進一步探究與驗證。本實驗研究還明確了自由來流壓力脈動在擾動模態計算過程中的重要性。

傳統的通過皮托探頭數據推導來流靜壓的方法,如Harvey及Stainback[33]分別基于定常與非定常過程的轉換,均與目前直接數值模擬的結果呈現極大差距,甚至是相反的趨勢,直接影響了擾動子模態的幅值。因此建議未來采用非介入式的光學測量方法對該問題進行深究。

最后,由于本文所使用測量方法與常規單一測量方法(如PIV或者熱絲測量)不同,使用兩種不同的測量技術在不同車次的風洞中測量所引入的誤差也更大,如何定量確定本方法的誤差范圍也將是接下來的研究重點。

猜你喜歡
模態測量
把握四個“三” 測量變簡單
滑動摩擦力的測量和計算
滑動摩擦力的測量與計算
測量的樂趣
車輛CAE分析中自由模態和約束模態的應用與對比
測量
國內多模態教學研究回顧與展望
高速顫振模型設計中顫振主要模態的判斷
航空學報(2015年4期)2015-05-07 06:43:35
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 亚洲人在线| 欧美精品一区二区三区中文字幕| 91精品国产综合久久香蕉922| 精品国产一区二区三区在线观看| 2020精品极品国产色在线观看 | 99久久免费精品特色大片| 福利姬国产精品一区在线| 日韩a在线观看免费观看| 在线精品欧美日韩| 老色鬼欧美精品| 人妻精品久久无码区| 福利一区在线| 国产精品2| 国产自产视频一区二区三区| 一级片一区| 亚洲欧美国产高清va在线播放| 亚洲乱强伦| 国产成人夜色91| 久久久久久尹人网香蕉| 在线观看国产小视频| 国产夜色视频| 三级国产在线观看| 国产精品无码一二三视频| 亚洲黄网视频| 亚洲美女一级毛片| 伊人五月丁香综合AⅤ| 国产精品永久不卡免费视频| 国产青榴视频| 国产精品一区二区在线播放| 日本在线国产| 国产va免费精品观看| 91精品专区国产盗摄| 亚洲午夜国产片在线观看| 青青久久91| 波多野结衣AV无码久久一区| 国产一区二区三区在线观看视频 | 国产综合欧美| 天天综合色网| 国产又色又刺激高潮免费看| 18禁色诱爆乳网站| 无码区日韩专区免费系列| 久久semm亚洲国产| 色婷婷在线播放| 亚洲动漫h| 国产美女91视频| 国产日韩欧美一区二区三区在线 | 色综合久久88| 国产精品一区在线麻豆| 国产无遮挡猛进猛出免费软件| 亚洲AV无码久久天堂| 在线播放真实国产乱子伦| 国产精品yjizz视频网一二区| 精品视频第一页| 国产在线专区| 国产精品福利导航| 中国美女**毛片录像在线| h视频在线播放| 日韩天堂视频| 久久精品91麻豆| 午夜精品一区二区蜜桃| 亚洲人精品亚洲人成在线| 在线精品亚洲一区二区古装| 日韩不卡免费视频| 亚洲高清无在码在线无弹窗| 午夜激情福利视频| 免费人成网站在线观看欧美| 久久久久人妻一区精品色奶水| 在线中文字幕日韩| 日韩国产综合精选| av在线5g无码天天| 日本在线亚洲| 22sihu国产精品视频影视资讯| 亚洲天堂伊人| 91破解版在线亚洲| 狼友视频一区二区三区| 香蕉视频在线观看www| 福利在线免费视频| 在线欧美a| 国产99视频在线| 欧洲高清无码在线| www中文字幕在线观看| 又污又黄又无遮挡网站|