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

有關近空間高超聲速飛行器邊界層轉捩和湍流的兩個問題

2017-04-28 04:52:37張涵信
空氣動力學學報 2017年2期

周 恒, 張涵信

(1.天津大學 力學系, 天津 300072; 2. 中國空氣動力研究與發展中心 國家計算流體力學實驗室, 北京 100191)

?

有關近空間高超聲速飛行器邊界層轉捩和湍流的兩個問題

周 恒1,*, 張涵信2

(1.天津大學 力學系, 天津 300072; 2. 中國空氣動力研究與發展中心 國家計算流體力學實驗室, 北京 100191)

和一般的飛行器一樣,在近空間飛行器的研制中,其邊界層的轉捩和湍流也是需要考慮的兩個重要問題。但即使是對一般的飛行器,“轉捩”和“湍流”也還是兩個歷經百年而仍未很好解決的問題,而對近空間飛行器來說,空氣動力學本身就還存在若干新的需要研究解決的基礎問題,邊界層的轉捩和湍流就更是沒有很好解決的問題。本文討論了兩個問題:1) 為增強對高超聲速飛行器邊界層轉捩預測的能力,需要開展哪些方面的研究工作及其困難;2) 是否有可能當飛行器飛行高度足夠大時,其邊界層就不會再有湍流問題?

近空間飛行器;轉捩;湍流;預測能力

1 轉捩預測問題

轉捩問題歷經百年的研究,在低速流方面已經有了不小的進展。而對高超聲速邊界層的轉捩問題,近年來的研究雖有一些進展,但仍存在嚴重的不足[1-2]。

可以肯定,整個轉捩過程可分為三個階段或作為三個問題來研究。1) 感受性問題,即邊界層外的擾動如何引發邊界層內以不穩定波形式出現的擾動?2) 被引發的不穩定波在邊界層內的演化。如果開始時幅值很小,例如一般認為,如果不穩定波的速度的幅值小于邊界層外緣速度的1%,則可以用線性理論;而大于1%,則要用非線性理論。3) 擾動演化到什么情況將觸發轉捩?即預測轉捩發生的判據。

以上的第二個問題研究的時間最長也最充分。特別是現在有了計算機和各種計算方法,只要知道擾動的初始條件,已經沒有原則上的困難了。

第三個問題原來是一個難題。長期以來人們提出多種邊界層內擾動演化的非線性理論,但始終無法和第三個問題,即演化到什么程度將觸發轉捩,聯系起來。在文獻[2]中,作者通過多種轉捩的直接數值模擬,發現了轉捩過程的機理是:擾動通過非線性作用對平均流進行修正,使得平均流剖面的線性穩定性發生變化。如果修正的結果是線性不穩定的參數區(例如線性穩定性理論的中性曲線所包含的區域)擴大,則會形成不穩定性的正反饋,從而觸發轉捩。他們并總結出,在技術人員喜歡用的半經驗轉捩預測方法eN法中,當不穩定波幅值達到1.5%時轉捩將發生。這可以認為第三個問題在一定程度上已經解決。

現今已有不少作者從事這方面的計算研究工作,例如Zhong等采用激波裝配法[4-9],Balakurma等則采用了激波捕捉法[10-11],然而他們對激波捕捉法的正確使用和激波后的擾動的非定常分析均存在缺陷。有文給出了線性擾動解,還用平行流假定作聲波、熵波和渦波分析,對高超聲速繞流這都是欠適當的。對有復雜激波系的問題,用激波裝配法來研究擾動通過激波經歷的變化是不現實的,而用激波捕捉法是否能可靠地反映這一變化則有待于新的能證實其可靠性的研究結果。

因此,為了解決轉捩預測問題,當前最需要做的基礎研究就是超聲速/高超聲速條件下的感受性問題。

前文已述,解決感受性問題的一個攔路虎就是飛行環境下的擾動是什么不清楚,為此就要做實地測量的飛行試驗。

從Zhong等已有的直接數值模擬結果來看,激波后的聲波在感受性問題中起主要作用,而激波前的聲波又在產生激波后的聲波中起主要作用。但在流體中,聲波的產生只有三種來源:一是由流體中的外來物體的體積周期性變化所引起;二是由外來物體對流體的周期性作用力所引起;三是由流體本身的運動產生了周期性的雷諾應力所引起。而在高空,實際上不存在外來物體(不算飛行器本身),因此可以認為不存在前兩種聲源。但第三種聲源激發聲波的效率很低,可以忽略不計。所以,飛行試驗的目標就是測出大氣中的以速度和溫度擾動出現的渦波和熵波。

但這是非常困難的任務。首先,如果用超聲速/高超聲速飛行器作為儀器載體,則儀器只能感受飛行器頭激波后的擾動。要還原成大氣中原來的擾動,是要解一個反問題,而反問題的解是不唯一的。其次,儀器的傳感器要有非常高的靈敏度。例如,一般估計對馬赫數為20的飛行器,速度背景擾動的幅值為萬分之一的量級。而如果馬赫數為6時,同樣的擾動對應的無量綱幅值就是萬分之三的量級。因此,測量時儀器就要有能分辨速度的萬分之一或萬分之三的變化的能力。而頻率范圍則要高達幾十萬赫茲,對溫度傳感器的要求也類似。如果要用亞聲速的飛行器作為儀器載體,則對儀器傳感器的靈敏度和頻響范圍的要求低得多,但如何將該飛行器送上高空后能有一段亞聲速的平飛段是一難題,而且亞聲速的飛行器本身對大氣流動的干擾會反映到傳感器上,如何排除?有人認為在30 km以下似乎可以用探空氣球作為儀器載體。但探空氣球的存在改變了周圍相當大的一片區域的流場,儀器的傳感器能否伸出到離氣球足夠遠的地方?

通常,在飛行器的設計中,地面的風洞實驗起了很大的作用,但要將地面風洞作為工程問題的轉捩預測的手段,是不夠的。如前所述,人們早就認識到對超聲速流而言,小擾動可分解為三種擾動。風洞的實驗段收縮比加大后,三種擾動更小,于是就認為風洞氣流幾乎沒有聲、熵、渦擾動。這在最初只是從減小背景擾動幅值的結果來考慮的。但現在知道,從感受性的角度看,聲波擾動、渦波擾動及熵波擾動的作用是不同的。同樣幅值的聲波和渦波導致的轉捩結果是不同的[3]。而上面的討論中已提到,高空大氣中本質上不存在聲波。而靜風洞的出現則在不自覺的條件下正好滿足了這一要求(靜風洞的出現比文獻[3]中所紀錄的聲波擾動和渦波擾動導致的小鈍頭錐上的轉捩線定性上就不同的研究結果要早得多)。但如果要模擬真實的轉捩,則靜風洞中的剩余擾動必須和真實情況下的背景擾動一致(無量綱化后一致),而我們實際上并不知道實際擾動是什么(如上所述,需要由飛行試驗測定)。而且即使知道了,要在風洞中實現(至少是對轉捩起作用的擾動頻率段的譜要一致)也是很大的難題。

此外,還有其他問題使得現有的風洞不能完全滿足預測工程問題中的轉捩的要求。這是因為:1)N-S方程無量綱化后,和不可壓縮流的情況不同,方程的系數并不是常數,而還有無量綱化后的粘性系數和普朗特數。粘性系數和溫度有關,而且對高超聲速流來說,實際飛行中邊界層的溫度還可能很高,從而粘性系數還要受氣體分子內能被激發、離解、甚至化學反應的影響。因此,要做到模型和實際飛行器的無量綱方程中的粘性系數一致,只有二者的溫度場完全一致才行,而這在現有的風洞中是不可能做到的。2)模型和實際飛行器的表面溫度還要一致,因為邊界層內的不穩定波的演化受溫度邊界條件的影響是很大的,而這一項的要求顯然也不可能在實驗中實現。因此,做到風洞試驗和工程實際的完全相似目前是做不到的。

既然這樣,那風洞實驗還能起什么作用呢?我們認為,目前它能起的最大作用是能提供有較完備參數的實驗結果,可以作為檢驗轉捩預測方法的依據。即將轉捩預測方法先直接面對風洞和模型的參數,而不是針對實際飛行器的參數,以檢驗其預測轉捩的結果是否和風洞實驗結果相符。但這樣做,就要求給出風洞實驗段的背景擾動的頻率譜及模型壁面溫度(不是實驗開始前的溫度,而是實驗過程中的溫度)分布。

2 湍流問題

迄今為止,在近空間以高超聲速飛行的飛行器是否必須考慮湍流問題還沒有一個可靠的說法。實際飛行實驗的結果非常少。經常被引用的只有幾篇文章,都還是20世紀七、八十年代的。文獻[12]是關于一個長近4 m、半錐角為5°的小鈍頭錐,從約100 km高空返回地面的過程中做的試驗。從100 km高空到約13 km高空,其飛行馬赫數都在20左右。結果發現一直降到30 km左右,才開始在尾部出現湍流。隨著高度不斷降低,轉捩線往前移,即湍流區不斷擴大。文獻[13]是1974年發表的,研究對象為長67.23 cm、鼻錐半徑16.09 cm、底部直徑48.69 cm的9°鼻錐。由拖曳火箭把它送上天,以5 km/s、迎角12.5°再入飛行,錐面有泰氟隆涂層,實驗結果用粗糙度稍大的轉捩準則處理。當飛行高度是36.58 km以上,氣流全為層流;30.48 km時,湍流區已擴展至鼻端,即高度為36.58 km以下會有湍流的影響。文獻[14]也有類似的結果。后來美國根據自由飛和風洞實驗結果,提出了不同粗糙度轉捩準則的轉捩公式[15]。我們要求近空間飛行器和航天飛行器一樣,基本上是光滑的,即實際粗糙度比上述實驗的粗糙度小,則可以合理地認為,開始發生轉捩的高度至少應在40 km到45 km間。2015年又有幾篇有關美國HIFiRE-1計劃的分析試驗報告,但其飛行高度值僅達到稍高于25 km。

型號設計中總提出一個問題:什么高度就可以不考慮湍流和轉捩影響,特別是近空間滑翔飛行器?從上面的分析來看,似乎可以認為飛行高度在45km以上就可以不考慮轉捩和湍流問題,但實際上問題也許不是這么簡單。

首先,為了驗證我們的結論,我們努力收集了國外飛行器的飛行結果,包括雙子星座、Apollo、美國各種航天飛機、蘇聯的飛船、英國的Hermes航天器、日本的小型飛船、美國的HTV-2以及中國的飛船,特別是哥倫比亞失事調查、復飛以及HTV-2前部的破壞等,都沒有發現是因為在40~50 km以上沒有考慮轉捩和湍流所導致的。最近美國正在改造世界最大的高超聲速風洞(九號),其一個重要目的是能夠做邊界層自然轉捩的實驗,他們把對應的飛行高度上限定為50 km。 以上是從已有的實際情況總結出來的結論,或許我們也可以從另一個角度對此問題做一些分析,即從氣體稀薄程度來判斷是否會有湍流。

本文作者曾在另一篇文章中[16]指出,對不同的問題,判斷是否要考慮氣體的稀薄效應要用不同的標準,否則按早年的觀點,高度在80 km以下的范圍,對飛行器的設計來說,都可以用連續介質模型。2005年黃章峰等在他們的一篇文章[17]中,報告了他們所做的一個在30 km高空、以馬赫數4.5飛行的平板邊界層的湍流直接數值模擬。結果發現,所得到的湍流和低速邊界層的湍流有幾乎完全相同的特征。例如,在壁面區存在低速條紋,其間距約為100個粘性長度,此間距會隨離壁面的距離增加而逐漸加大。在湍流近壁區,存在相干結構,其特征也和低速邊界層中的湍流相似。因此,我們可以假定在近壁區出現低速條紋,且條紋間距約為100粘性長度,是湍流邊界層的普適現象。在判斷是否有湍流問題時,判斷氣體稀薄與否的參考長度也許就應取這一長度。

我們也可以先假定連續介質模型可用,對在更高空處飛行的平板做類似的計算,但要計算多種工況則工作量較大,而采用以下的從層流解來估計則要方便得多。

設在某一高度有一平板以0°迎角及某一速度飛行,則可按相似性解得到層流邊界層的解,從而得到相應的粘性長度。

例如,對位于10 000 m和30 000 m高空的平板層流邊界層,采用相似性解,得到馬赫數為0.85和4時,位于離前緣2 m處的粘性長度l+見表1。

表1 不同高空、不同馬赫數下的粘性長度Table 1 Viscous length at different altitude and different Mach number

由以上結果看,馬赫數和壁面的溫度條件對粘性長度的影響不大,而高度則影響較大。高度從10 000 m增至30 000 m,l+就增加為原來的約4.8倍。

黃章峰的DNS顯示,從層流轉捩為湍流時,壁面摩擦系數增為原來的3~4倍。其他人對不可壓邊界層做的轉捩直接數值模擬,也得到類似的結果。如果取3,則對應于湍流邊界層,粘性長度單位約相當于層流的(1/3)1/2≈0.58倍。因此,對湍流邊界層來說,粘性長度和高度變化的規律也不變,在高度從10 000 m增加到30 000 m時,也是增為原來的4.8倍,或高度每增加10 000 m,粘性長度約增加約2.2倍。在比30 000 m更高的高空,這一增加的比例變化不會有量級的變化。

但是從10 000 m高空開始,高度只要每增加10 000 m,大氣分子自由程就增為原來的約4.7倍。因此,高度每增加10 000 m,湍流邊界層近壁區的低速條紋間距和當地分子自由程之比就要降為原來的約2.2/4.7=0.47倍。當這個比值降到一定的程度,低速條紋就無法維持了,即不會有湍流了。

10 000 m高空分子自由程約為0.22μ,30 000 m高空時約為4.7μ。對300 K壁溫情況:馬赫數為4時,在靠近壁面處,溫度升高很少,當地分子自由程和來流的分子自由程基本相同。而相應地湍流邊界層低速條紋間距(0.58×100l+)則分別約為1.322 mm和6.322 mm,間距和分子自由程之比分別為6009和1345。

如果馬赫數達20,則對定溫壁(300 K),在層流邊界層厚度的1/8處,溫度約為來流溫度的11倍[18](對湍流邊界層應更高一些)。相應地,對30 000 m高空,該處的分子自由程可達52μ,或0.052 mm。一個低速條紋間距對應僅121個分子自由程。

按一般湍流研究結果看,湍流邊界層近壁區的低速條紋是一種“大”結構,充分發展湍流中還應有比它尺度至少小一個量級的小結構。而要形成一個結構,其尺度又不能太小,必須包含足夠多的分子,或者說,其尺度應該至少比分子自由程大一個量級。因此,低速條紋的間距至少應該比分子自由程大兩個量級,或約100倍。否則,至少連續介質模型下的那種湍流是否能存在要存疑。但這也許還不能完全否定可能存在湍流。事實上,也許在一般的湍流和完全沒有湍流之間存在一個過渡區,就像從層流到湍流存在一個轉捩過程那樣。在那個范圍內,可能存在既非層流又非充分發展的湍流。

同樣地分析,到40 000 m高空時,一個低速條紋僅對應約56個分子自由程。從物理上講,應已不能形成通常意義下的湍流。

如果近空間飛行器采取的是耐高溫材料防熱,則壁面接近于絕熱情況。馬赫數為8時,壁面溫度即可達來流溫度的約11倍(層流相似性解),即在40 km高空已不會有湍流。

由于長航程的高超聲速飛行器一般都要飛行于較高的高度,例如至少是40 km或更高,飛行馬赫數也較大,因此可能的確不需要考慮湍流和轉捩的問題。

但是,不能無條件地應用這一結論。

實際上,前面的分析都是針對0°迎角飛行的平板而做的,其隱含的假定就是邊界層內的壓力和環境壓力相同。對有迎角平板,情況就不一樣了。

例如,假定平板以速度v,迎角θ飛行,則按牛頓理論,平板迎風面單位面積所受壓力約為:

在海平面,空氣單位重量約為1.29kg/m3,大氣壓約為10 339kg/m2。在30 km高空,空氣單位重量和大氣壓力分別降為約0.0189kg/m3和119kg/m2。而在50 km高空,空氣單位重量約為0.000 89kg/m3(以上數據不是精確值,因為高空的大氣溫度、密度和壓力都是不斷變動的)。如果馬赫數為15,則其速度約為4500 m/s,代入式(1),得

如果要求和30 km高空處大氣壓相當,即約為119kg/m2,則Sin2θ=0.0647,或θ=14.7°。即迎角達到14.7°時,其邊界層內的壓力就和在30 km高空以同樣速度飛行但迎角為0°的平板一樣,也就有可能發生轉捩。

對于像錐體那樣的飛行體,即使有迎角,在迎風面的壓力升高顯然遠小于平板的情況,因此在飛行試驗中未發現高度超過30 km時發生轉捩也是可以理解的。但對有升力面的飛行器,則在升力面上,要根據實際邊界層中的壓力和溫度,看邊界層中的分子自由程和按連續介質模型所得湍流邊界層中低速條紋間距(約100個粘性長度)相比的值。可以考慮暫時將這一臨界值定為1%,大于它則多半不能維持湍流,也就不存在轉捩問題了。

3 結束語

上面討論了研制近空間飛行器需要研究考慮的兩個空氣動力學問題。給出了它們的研究現狀,提出了為解決它們需要開展的研究及可能遇到的困難。由于問題的復雜,也由于已有的研究成果很少,因此提出的看法和建議不能說是完全正確的,也不能說是很全面的。隨著研究的開展,可能還會發現新的問題,或者發現文中的觀點不準確甚至不正確,這都是正常的。如果我們的觀點或建議對今后的研究中能起到一定的促進作用,我們就很滿意了。

[1]Berin J J, Cummings R M. Critical hypersonic aerothermodunamic phenomena[J]. Annual Review of Fluid Mechanics, 2006, 38: 129-157.

[2]周恒, 蘇彩虹, 張永明. 超聲速/高超聲速邊界層的轉捩機理及預測[M]. 北京: 科學出版社, 2015: 1-49, 88-96.

[3]Pope A, Goin V L. High speed wind tunnel testing[M]. United Kingdom: Wiley, 1965.

[4]Zhong X. Leading-edgereceptivity to free-stream disturbance waves for hypersonic flow over aparabola[J]. Journal of Fluid Mechanics, 2001, 441: 315-367.

[5]Ma Y, Zhong X. Receptivity of asupersonic boundary layer overaflatplate. Part1. Waves tructures and interactions[J]. Journal of Fluid Mechanics, 2003, 488: 31-78.

[6]Ma Y, Zhong X. Receptivity of asupersonic boundary layer overaflatplate. Part2. Recep-tivity to free stream sound[J]. Journal of Fluid Mechanics, 2003, 488: 79-121.

[7]Ma Y, Zhong X. Receptivity of asupersonic boundary layer overaflatplate. Part3. Effects of different types of free-stream disturbances[J]. Journal of Fluid Mechanics, 2005, 532: 63-109.

[8]Ma Y, Zhong X. Boundary-layer receptivity of Mach 7.99 flow overabluntcone to free-stream acoustic waves[J]. Journal of Fluid Mechanics, 2006, 556(1): 55-103.

[9]Zhong X, Wang X. Directnumerical simulation on there ceptivity, instability and transition of hypersonic boundary layers[J]. Annual Review of Fluid Mechanics, 2012, 44(1): 527-561.

[10]Balakumar P. Receptivity of asupersonic boundary layer to acoustic disturbances[J]. AIAA J., 2009, 47(5): 1069-1078.

[11]Balakumar P, Kegerise M A. Receptivity of hypersonic boundary layers over straight and flared cones[J]. AIAA J., 2010, 53.

[12]Wright R L, Zoby E V. Flight measuement of boundary layer transition of a 50 half angle cone at free stream Mach number of 20 (reentry)[R]. NASA TM X-2253, 1971.

[13]Hayer D T, Herskovitz S B, Lennon J F, et al. An ablation technique for enhancing reentry antenna performance, flight test results[R]. FACRL-TR-74-0572, ADA012250, 1974.

[14]Hayer D T, Herskovitz S B, Lennon J F, et al. Flight test data comparing electron attachment by ablation products and by liquid injection. AIAA-75-181[R]. Reston: AIAA, 1975.

[15]Batt R Q, Legner H H. A review of roughness-induced nosetiptransition[J]. AIAA J., 1983, 21(1).

[16]周恒, 張涵信. 空氣動力學的新問題[J]. 中國科學(中文版), 2015, 40(10): 104709.

[17]Huang Zhangfeng, Zhou Heng, Luo Jisheng. Direct numerical simulation of a supersonic turbulent boundary layer on a flat plate and its analysis[J]. Science in China G, 2005, 48(5): 626-640.

[18]Anderson J D Jr. Hypersonic and high-temperature gas dynamics[M]. Second Edition, AIAA Education Series, 2006, 293-296.

Two problems in the transition and turbulence for near space
hypersonic flying vehicles

Zhou Heng1,*, Zhang Hanxin2

(1.MechanicalSystemofTianjinUniversity,Tianjin300072,China;2.NationalLaboratoryforComputationalFluidDynamics,ChinaAerodynamicsResearchandDevelopmentCenter,Beijing100191,China)

For the research and development of near space flying vehicles, also as the same for conventional flying vehiclessuch as airplanes, the problems of transition and turbulence of the boundary layers are two important issues must be taken into consideration. However, even for conventional flying vehicles, these two problems are still not fully resolved, even though the investigations have been lasted for more than 100 years already. For near space flying vehicles, not only the related aerodynamics has its own unsolved fundamental scientific problems, let along the problems of transition and turbulence. In this paper, two related problems are focused on: 1) In order to enhance our capability of predicting the transition of the boundary layer of a hypersonic flying vehicle, what kinds of research work should we do and what difficulties we might face? 2) Would it be possible that there would be no problem of turbulence for its boundary layer if the attitude of the flying vehicle is sufficiently high?

near space flying vehicles; transition; turbulence; capability of predicting

0258-1825(2017)02-0151-05

2017-02-17;

2017-03-23

周恒*(1929-),男,中國科學院院士,研究方向:邊界層轉捩,流動穩定性非線性理論. E-mail: hzhou1@tju.edu.cn

周恒, 張涵信. 有關近空間高超聲速飛行器邊界層轉捩和湍流的兩個問題[J]. 空氣動力學學報, 2017, 35(2): 151-155.

10.7638/kqdlxxb-2017.0016. Zhou H, Zhang H X. Two problems in the transition and turbulence for near space hypersonic flying vehicles[J]. Acta Aerodynamica Sinica, 2017, 35(2): 151-155.

V211.3

A doi: 10.7638/kqdlxxb-2017.0016

主站蜘蛛池模板: 久久久亚洲国产美女国产盗摄| 亚洲福利一区二区三区| 成人看片欧美一区二区| 日本国产精品一区久久久| 国产微拍精品| 黄色成年视频| 久996视频精品免费观看| 精品福利国产| 热久久综合这里只有精品电影| 夜夜操天天摸| 国产自在线拍| 国产成年女人特黄特色大片免费| 国产一区免费在线观看| 色呦呦手机在线精品| 欧美亚洲激情| 国产迷奸在线看| 午夜性刺激在线观看免费| 国产免费高清无需播放器| 久久人人爽人人爽人人片aV东京热| 91在线播放免费不卡无毒| 72种姿势欧美久久久久大黄蕉| 无码一区二区波多野结衣播放搜索| 欧美日韩中文国产va另类| 日韩在线2020专区| 亚洲无码高清免费视频亚洲| 久久这里只有精品66| 99精品影院| 国产精品亚洲片在线va| 五月婷婷伊人网| 一区二区三区高清视频国产女人| 91po国产在线精品免费观看| 少妇极品熟妇人妻专区视频| 天堂在线视频精品| 自慰高潮喷白浆在线观看| 精品三级在线| 日韩在线网址| 久久亚洲天堂| 亚洲人成影院午夜网站| 一边摸一边做爽的视频17国产 | 天堂网亚洲综合在线| 久久精品视频一| 日本精品中文字幕在线不卡| 国产福利小视频高清在线观看| 国产无码高清视频不卡| 真实国产乱子伦高清| 欧美伊人色综合久久天天| a级毛片在线免费| 日韩精品高清自在线| 一本色道久久88综合日韩精品| 亚洲系列无码专区偷窥无码| 丝袜无码一区二区三区| 操国产美女| 在线观看的黄网| 国产亚洲视频中文字幕视频| A级毛片高清免费视频就| 国产爽歪歪免费视频在线观看 | 国产日本欧美亚洲精品视| 欧美在线一级片| 日韩第九页| 动漫精品中文字幕无码| 国产亚洲精品自在线| 国产麻豆永久视频| 欧美97色| 国产高清又黄又嫩的免费视频网站| 麻豆国产精品视频| 国产成人a毛片在线| 中国一级特黄大片在线观看| 久久精品女人天堂aaa| www.亚洲国产| 91无码人妻精品一区| 一级毛片免费观看久| 国产精品视频猛进猛出| 久久香蕉国产线看观看精品蕉| 国产中文一区a级毛片视频| 婷婷六月激情综合一区| 久久国产亚洲偷自| 一级毛片免费播放视频| 色有码无码视频| 久久久久久久久18禁秘| 97se亚洲| 久久精品这里只有精99品| 欧美日韩午夜视频在线观看|