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

類乘波構型機身的設計及鈍化對其性能的影響 ①

2020-08-01 00:54:06張艷華張登成鄭無計
固體火箭技術 2020年3期
關鍵詞:設計

羅 浩,張艷華,張登成,鄭無計

(空軍工程大學 航空工程學院,西安 710038)

0 引言

乘波體在高超聲速時具有高升力、高升阻比的優點,可以突破常規布局所面臨的“升阻比屏障”問題,因此成為高超聲速飛行器的首選布局。雖然各航天大國針對乘波構型進行的大量研究證明了其良好的高超聲速性能,但乘波構型也有自身的缺陷,如容積率低、偏離設計點性能下降、外形加工難等。因此,在理想乘波構型的基礎上,根據實用化要求對其進行改進是高超聲速飛行器氣動布局設計所必須考慮的,而且乘波構型的尖銳前緣在高超聲速時在面臨著嚴重的氣動熱環境,為了滿足防熱需求,對其邊緣實施鈍化是有效的途徑[1]。

在類乘波構型設計方面,范曉檣[2]初步建立了高超聲速二維進氣道的參數化設計方法,確定了進氣道參數化設計的5個基本幾何控制參數,為二維進氣道型面優化設計提供了可行的途徑。李曉宇[3]在二維設計的基礎上研究了三維效應對于前體流場及預壓縮性能的影響,提出了抑制三維側緣溢流影響、改善三維前體氣動性能的前體側緣設計思想。對于乘波體鈍化的研究,劉建霞[4-5]提出了一種新的邊緣鈍化方法,數值模擬和風洞試驗的結果表明這種鈍化方法在滿足防熱需求的同時氣動性能下降很小。文獻[6]分析了一致鈍化和非一致鈍化兩種典型鈍化方法的優缺點,指出了非一致鈍化方法的優越性。陳小慶[7]和李世斌[8]等研究了不同鈍化半徑和位置對乘波構型氣動力和氣動熱的影響,指出鈍化半徑大小的選擇要綜合考慮氣動力和氣動熱性能。文獻[9]提出了一種半徑可控參數化三維前緣鈍化方法,并且在脊形乘波體上進行了驗證,這種方法保證了鈍化前緣與上下表面光滑連接,鈍化可以參數化有利于程序實現。

現有的大部分文章針對乘波構型提出了一系列鈍化方法并研究了鈍化對氣動力和氣動熱性能的影響,推進了乘波構型的實用化發展,但因其固有的橫截面曲線復雜、容積率小和與發動機結合困難等缺陷,在航空航天領域的實際應用中都受到了一定的限制。本文設計了一種有較高實用性的類乘波構型高超聲速飛行器機體,在一定程度上解決了乘波構型在實際應用中存在固有缺陷的問題,并根據防熱需求和構型特點提出了一種鈍化方法,利用CFD方法對其鈍化后的氣動力和氣動熱性能進行了評估,研究可為類乘波構型高超聲速飛行器的設計和外形鈍化提供參考。

1 類乘波體構型的設計

1.1 二維流道前體設計

二維流道設計是三維構型設計的基礎,主要包括前體設計和后體設計。二維氣動外形的示意圖如圖1所示,飛行器的尺寸參數設定為,機身總長為2900 mm,前體長度為1300 mm,發動機段長度為700 mm,后體長度為900 mm,機身寬度為670 mm,機身高度為260 mm。下面將分別對前體和后體進行設計。

圖1 類乘波構型機體縱向截面氣動布局

前體的設計目標是在壓縮過程中使總壓損失和阻力盡量小的同時提高壓縮效果和流量捕獲率。由于楔面前體設計具有激波系簡單、出口流場均勻、熱防護容易實現和性能較好等優點,已經被應用到大多數高超飛行試驗方案和發動機地面試驗方案之中[10-11]。因此,本文采用基于斜激波理論通過二維楔面橫向擴展得到前體的方法進行前體構型的設計。

本文中飛行器在設計狀態下的飛行條件為H=26 km,Ma=6。借鑒經過飛行試驗[12]和現有較為成熟的前體設計方法[13],本文設計的前體采用三級轉折角。為保證較大的總壓恢復系數,采用各級轉折角等總壓恢復系數的原則來配置各級轉折角,給定前體的總壓恢復系數σ=0.9,則每級轉折角的斜激波總壓恢復系數σi為

(1)

(2)

式中σi為各級轉折角的總壓恢復系數;Mai1為第i道斜激波的波前馬赫數;β為激波角;k為氣體的絕熱指數。

壓縮效果表現為壓力上升、速度下降。參考文獻[14]中對壓縮效果的要求,即進氣道出口的壓力大于1個大氣壓,高度H=26 km時大氣壓力為2152.06 Pa,也就是壓升pr≥47.1。

(3)

(4)

(5)

式中pr為壓升;δ為氣流轉折角;Ma1為波前馬赫數;Ma2為波后馬赫數。

流量捕獲率是衡量進氣道捕獲壓縮氣流能力的主要指標。為提高流量捕獲率,常用的方法是合理配置前體各級轉折角度使斜激波系在設計狀態下相交于進氣道唇口前緣處,這種設計方法被稱為“shock on lip”方法[14]。如圖2所示,在“shock on lip”狀態下,前體其他各級的長度可通過第一級長度和各級轉折角計算得到,其關系如下所示:

(6)

Ylip=Xliptanβ1

(7)

(8)

L1+L2+L3=1300 mm

(9)

式中Xlip為前體頂點到唇口的水平距離;Ylip為唇口點到前體上表面的垂直距離;Li為前體各級長度;δj為前體各級轉折角;βi為各級轉折的斜激波角度。

根據以上設計原則可得到前體設計的步驟,如圖3所示。經過這個流程最終可獲得既滿足進氣道進口壓縮效果又具有較大總壓恢復系數的前體各級轉折角和各級長度,最終得到總壓恢復系數σ=0.84,壓升pr=47.16,前體總轉折角δtotal=17.61°,前體各級轉折角分別為δ1=5.36°、δ2=5.57°、δ3=6.68°,前體各級長度分別為L1=546 mm、L2=286 mm、L3=468 mm。

從上面的計算過程中可以發現壓升和總壓恢復系數要求是相互矛盾的,這是因為壓縮效果和總壓恢復系數的主要影響因素是不同的,因此在進行前體設計時需要綜合考慮這兩者的要求,在兩者之間進行權衡處理。

1.2 二維流道后體設計

后體的設計目標是使高焓氣流在下表面盡可能膨脹的同時保證滿足升力、穩定性等氣動特性要求。本文的后體下表面采用斜切形式的單壁面三次曲線,考慮到實際要求后體尾部有一定的高度,為機身高度的0.1倍。后體的簡圖如圖4所示。

后體主要參數有后體長度L,高度H,尾部高度M,上表面初始膨脹角θ1,上表面出口切角θ2。后體下表面為三次曲線:

y=a1x3+a2x2+a3x+a4

(10)

其中,系數a1、a2、a3、a4由L、H、θ1、θ2確定。由于后體長度L和高度H已經確定,因此只需確定θ1和θ2就可以得到三次曲線。本文參考文獻[15]優化后的上表面初始膨脹角θ2和上表面出口切角θ2,取θ1=26.3°、θ2=5.5°。二維后體的型面參數分別為L=900 mm、H=260 mm、M=26 mm。

圖5給出了在設計條件下二維機身構型的壓力云圖,可見前體三級轉折面分別產生了3道激波,激波系結構保持較好、很好地匯聚到了一點。經過數值計算得出前體的總壓恢復系數σ=0.81,壓升pr=47.16,說明前體的設計達到了預期的效果,起到了很好的氣流壓縮作用。后體下表面的氣流從后體起始點開始膨脹,一直到遠場都處于膨脹區,后體產生了飛行器的大部分推力,但后體對飛行器的升力貢獻較小。

圖5 二維機體壓力云圖

1.3 側緣構型設計

本文采用的側緣構型為“上反”形式,這種構型截面簡單易加工、容積率高,但對橫向流動的抑制作用較差。側緣構型的設計參數為側緣后掠角α和側緣轉折角β,如圖6所示。為保證經過前體壓縮的氣流能完全進入進氣道內部,發動機段下表面的機身寬度和前體下表面的寬度保持一致,文獻[16]中X-43的構型尺寸,下表面的寬度為機身最大寬度的一半。結合上文給出的前體長度,可得到側緣后掠角;因為前體為三級轉折,故有3個對應的側緣轉折角。側緣構型具體參數為α=80.9°、β1=143.8°、β2=134.9°、β3=123.5°。

(a)Sweep angle of side edge (b)Turning angle of side edge

綜合上述各方面的設計最終確類乘波機身構型,如圖7所示。這種類乘波構型的橫截面曲線簡單、易于加工制造、機身容積大,易實現發動機與機身的一體化設計,是高超聲速飛行器機身可以使用的布局。

圖7 機身布局圖

2 類乘波體構型的鈍化

高超聲速飛行時,由于空氣劇烈壓縮飛行器會受到嚴重的氣動加熱,從上節類乘波構型機體的氣動布局圖可以看出,第一級外壓楔形體具有尖銳的前緣,側壁面和上、下壁面的連接處也是尖銳的,這些部位的熱腐蝕會非常嚴重,因此有必要對其進行鈍化處理。現有的鈍化方法都是針對分開乘波體上下表面的一個尖銳邊緣,本文中飛行器機體的側壁面和上、下壁面的連接處分別產生了一個尖銳的邊緣,在進行鈍化時需要綜合考慮兩個邊緣。

2.1 前緣鈍化

圖8 鈍化方法示意圖

2.2 邊緣鈍化

上邊緣線和下邊緣線不在同一個平面上,且下邊緣線是一條三維空間曲線,為降低鈍化的難度,可參考文獻[9]中的鈍化方法將復雜的三維曲面設計問題簡化為二維平面設計問題,即在各個橫向站立位置設計一系列的鈍化曲線,然后對其放樣得到三維鈍化曲面。對于上邊緣線和下邊緣線如果分開進行鈍化,在同一橫截面得到的兩條鈍化曲線融合設計的難度比較大,放樣后可能造成側面形狀凹凸不平,因此有必要結合側壁面對上、下邊緣線進行綜合鈍化處理。為保證鈍化曲線與上下表面光滑連接,本文中采用增加材料的鈍化方法,曲線類型采用三次Bézier曲線。三次Bézier曲線共需要4個控制點控制曲線形狀,如圖9所示。

4個控制點間的關系由式(11)表示:

(11)

2個端點處的曲率為

(12)

圖9 三次Bézier曲線

在P1、P4點處鈍化曲線分別于下表面和上表面相切,由此可得到關于δ1和δ2的方程組:

(13)

此處向量叉乘定義為

a×b=(a1,a2)×(b1,b2)=a1b2-a2b1

(14)

式中P1、P2、P3、P4為對應點的坐標向量;V1、V2分別為P1、P2點處的切矢量;k1、k2分別表示P1、P2點處的曲率;δ1、δ2表示變量。

通過求解式(12)和式(13)可以確定控制點P2和P3的坐標向量,控制點P1和P4的坐標向量已經確定,其分別位于橫向截面與上、下邊緣線的交點,從而可以得到一條光滑的鈍化曲線。

針對本文中的類乘波機身構型,以上邊緣線為例,具體的鈍化方法分為以下幾個步驟:

綜上所述,體外沖擊波聯合扶他林乳膠劑可更好地發揮體外沖擊波的物理生物學效應和扶他林消炎鎮痛的效果,并產生疊加效應,從而顯著提高OA患者的膝關節功能,減輕膝關節疼痛癥狀。但在臨床應用過程中沖擊波的沖擊強度、頻率以及治療周期等參數還有待規范,并且本研究觀察周期短,其長期療效及治療機制仍有待進一步研究。

(1) 首先確定各橫向站位處的鈍化半徑。文獻[19]中提出了一種用經驗公式確定鈍化半徑的非一致鈍化方法,可由駐點處的鈍化半徑和當地側緣后掠角確定橫向站位處的鈍化半徑,如式(15)所示:

r=rs(cosα)2.2

(15)

式中r為橫向站位處的鈍化半徑;rs為駐點處的鈍化半徑;α為當地側緣后掠角。

由類乘波體的構型特點可知,本文中的駐點位于第一級楔形體的尖銳前緣,前緣處的鈍化半徑已經確定為20 mm,故各個橫向站位處的鈍化半徑可以計算得到。

(2)在橫向站位處,上表面線和側表面線構成角γ,做出角γ的角平分線,在角平分線上做垂線,垂線分別于上表面線和側表面線相交,取垂線的長度d為2倍鈍化半徑r,即d=2r。

(3)垂線與上表面線的交點到邊緣點的距離為h,將側表面線沿著上表面線方向向外移動長度為h的距離,延長上表面線和下表面線構成新的側表面線。在新的上表面線和側表面線構成的夾角上重復步驟(2)。

(4)采用上述的三次Bézier曲線實施鈍化。步驟(3)中得到的垂線與上表面線和側表面線的交點分別對應三次Bézier曲線中的點P1和點P4,根據三次Bézier曲線的確定方法可以得到上邊緣線的鈍化曲線,即圖10中的紅色曲線1。

(5)對下邊緣線實施鈍化。在下邊緣線上重復上述步驟(1)、(2)、(4),即可得到下邊緣的鈍化曲線,即圖10中的紅色曲線2。

圖10 側緣鈍化方法示意圖

通過這個方法可在不同的橫截面站立位置分別得到連接側壁面和上、下壁面的光滑鈍化曲線,接著對鈍化曲線族進行放樣得到鈍化曲面。這種鈍化方法同時對上下邊緣線進行了處理,不僅增大了機身容積,也保證了鈍化曲線和上下邊緣的光滑連接。如圖11所示,黑色曲線表示原機身構型橫截面,紅色曲線表示鈍化后的側緣曲線,鈍化前后縱向截面的形狀對比和鈍化后的氣動布局分別如圖12和圖13所示。可以看出鈍化后機身容積有所增加。乘波構型容積率的計算公式定義為式(16),其中,V和S分別為乘波構型的容積和面積,不同文獻對于S有不同的取法,一般取為乘波構型的表面浸潤面積或者投影面積[21-22]。

圖11 鈍化前后側緣對比

圖12 縱向截面形狀對比

圖13 鈍化后氣動布局

本文參考文獻[22]將取S為乘波構型的投影面積。類乘波構型鈍化前后的機身容積、投影面積和容積率如表1所示。可見,鈍化后類乘波構型的容積和容積率都有所增加,其中容積增加了8.7%,容積率增加了3.6%。

(16)

表1 鈍化前后容積和容積率對比

3 類乘波體構型的鈍化

3.1 計算網格及數值方法

由于物理模型是對稱的且本文只研究氣動力和氣動熱特性,為節省計算時間,僅對半模的外流場進行數值計算。高超聲速飛行時,飛行器壁面流動復雜,模擬所有流動現象難度大,本文在計算中忽略了化學反應、高溫氣體效應等因素。氣動熱的數值計算對網格的依賴性高,本文的計算網格由ICEM軟件生成,鈍化前的網格為六面體結構網格;由于鈍化后側面是曲面,表面復雜因此采用非結構網格,在壁面附近生成棱柱網格來模擬附面層。根據文獻[23]中提出的高超聲速氣動熱數值計算壁面網格準則計算第一層網格高度,為提高邊界層的計算精度,在壁面附近進行加密處理,壁面法向第一層網格高度為5.9×10-6m,鈍化前后的計算網格分別如圖14和圖15所示。

圖14 鈍化前的計算網格

圖15 鈍化后的計算網格

網格數量會影響數值計算的時間和精度,網格數量的增加雖然會使流場模擬更加精細,但是在滿足流場結構正確后,再增加網格數量也不會提高計算的精確性,還會增大計算量,因此進行了網格無關性分析。鈍化前、后分別生成了3組網格,在馬赫數Ma=6,攻角α=4°條件下,分別對3組網格計算了鈍化前和鈍化后縱向對稱面處的駐點溫度,鈍化前后的網格數量和計算結果分別如表2、表3所示。可見鈍化前、后3組網格計算結果的誤差都不超過5%,為了節省計算機資源,本文選取鈍化前和鈍化后的網格數量分別為260萬、300萬。

表2 鈍化前3組結構網格和計算結果

表3 鈍化后3組非結構網格和計算結果

流動控制方程為可壓流的三維雷諾平均N-S方程,采用RNGk-ε方程湍流模型和壁面函數來模擬湍流,采用AUSM+格式進行空間離散,時間項采用LU-SGS隱式方法,壁面條件為無滑移等溫的邊界條件,遠場邊界為自由來流,計算條件如表4所示。

表4 計算條件參數

3.2 算法驗證

氣動熱的計算精度與湍流模型、離散格式和計算網格等密切相關,為了驗證數值方法的可靠性,本文與文獻[24]研究的鈍化雙錐體模型的風洞實驗結果進行了對比驗證。為了確保驗證的準確性和節約計算時間,根據不同的空間尺寸和幾何特征對鈍化雙錐體模型的關鍵部位進行局部加密,并進行網格的無關性驗證。來流條件為:Ma=12.43,p∞=18.27 Pa,T∞=107 K,Tw=297 K,對網格數量分別為300萬、350萬、400萬的3組網格劃分方法的熱流系數進行計算,在鈍化雙錐體X/L=1.2處的計算結果分別為0.112、0.117、0.121,3組網格之間的最大相對誤差均不超過10%,考慮到計算機的計算能力,最終采用網格數量為200萬的網格劃分方法。圖16(a)、(b)給出了實驗和數值計算的壓力系數和熱流系數分布,其中壓力系數和熱流系數分別定義為

圖16中,X為表面各點到頂點的距離;L為鈍化雙椎體的特征長度。可以看出,實驗值和數值計算結果吻合度較高、變化趨勢基本一致,故認為本文所采用的數值計算方法是可靠的。

(a)Pressure coefficient (b)Heat flux coefficient

4 類乘波體構型的鈍化

用上面的數值方法分別對類乘波構型鈍化前后的氣動力和氣動熱特性進行了計算分析,以研究鈍化對類乘波構型高超聲速飛行器機體氣動性能的影響。

4.1 鈍化對氣動力性能的影響

圖17給出了在馬赫數Ma=6條件下,鈍化前后類乘波構型的升力系數、阻力系數和升阻比隨攻角的變化趨勢。

圖17 鈍化對氣動力性能的影響

從升力系數曲線可以看出,鈍化前、后類乘波構型的升力系數都隨攻角增加幾乎呈線性增長,鈍化后升力系數略微減小。從阻力系數曲線可以看出,鈍化前、后類乘波構型的阻力系數都隨攻角增加呈非線性增長,這是因為攻角增加改變了激波的附體特性,導致激波阻力非線性增長;鈍化對阻力系數產生了較大影響,鈍化后阻力系數明顯增大。從升阻比曲線可以看出,鈍化前后升阻比隨攻角的增大都是先增大后略有減小,數值計算的結果表明,最大升阻比由鈍化前的5.26減小到鈍化后的4.58,減小了約13%。雖然鈍化使得升阻比減小,但鈍化后升阻比仍然保持在較高的數值,最大升阻比可到達4,說明這種鈍化方法使類乘波構型機身有良好的氣動特性。

圖18給出了馬赫數Ma=6,攻角α=4°時,鈍化前、后縱向對稱面的壓力云圖對比。如圖18(a)所示,鈍化前激波很好地貼附在下表面,下表面前體產生了高壓區,下表面壓力遠大于上表面,上下表面產生了明顯的壓力差,這是升力的主要來源。如圖18(b)所示,鈍化后類乘波構型前緣上表面的氣流受到了擾動,使激波不再完全依附在前緣線上,導致下表面的高壓氣體沿著鈍化前緣泄漏到上表面,降低了類乘波構型的氣動性能;但是激波并未完全脫體,下表面高壓區仍然十分明顯,乘波效應較為顯著。

(a)Sharp model (b)Blunt model

從圖19可以看出,鈍化后產生了脫體激波,激波的位置發生了改變,導致下表面的高壓氣體泄漏到上表面,使上表面邊緣壓力增大,減小了升力,但鈍化前后下表面的高壓氣流分布都很均勻。由圖20可見,鈍化前上下表面產生了均勻的壓力差,提供了穩定的升力,鈍化后上下表面的壓力差略有減小,鈍化對類乘波構型邊緣處的流場造成了影響,導致邊緣附近上下表面都產生了壓力波動,但鈍化后在展向0≤z/z0≤0.7范圍內,上下表面壓力分布仍較為均勻。

圖19 橫截面壓力云圖

圖20 橫截面展向壓力分布

4.2 鈍化對氣動熱性能的影響

圖21分別給出了Ma=6,攻角α=4°時,鈍化前后類乘波構型前部側壁面的熱流密度云圖。由圖21可見,鈍化前側壁面邊緣的熱流密度較大,鈍化后除駐點附近區域其余部分的熱流密度明顯減小;數值計算的結果表明,壁面最大熱流密度由鈍化前的6.16×106W/m2減小到鈍化后的1.44×106W/m2,減小了約77%,鈍化使熱流密度急劇減小,極大改善了氣動熱性能。邊緣的熱流密度不到駐點熱流密度的一半。鈍化對熱流密度的影響較大,提高了類乘波構型的氣動熱性能。

(a)Sharp model (b)Blunt model

從上面的數值計算和分析可知,鈍化降低了類乘波構型的氣動力性能,極大地改善了氣動熱性能,雖然氣動力性能有所下降,但是升阻比仍然保持在較高的數值,乘波特性較為顯著。

5 結論

本文設計了一種類乘波構型高超聲速飛行器機體,根據熱防護需求對其進行了修形鈍化,通過CFD的方法分別計算了鈍化前后類乘波構型的氣動力和氣動熱性能,分析了邊緣鈍化對類乘波構型氣動性能的影響,研究結果表明:

(1)按照等總壓恢復系數設計的前體較好地平衡了總壓恢復系數和壓縮效果要求,結合后體設計和側緣設計得到的類乘波構型高超聲速飛行器有良好的氣動性能。

(2)鈍化對類乘波構型阻力的影響大于對升力的影響,鈍化后升力系數略微減小,但阻力系數明顯增大;這是由于鈍化改變了激波結構,使下表面的高壓氣流向上表面泄漏,升力減小,激波阻力增加。鈍化可有效降低壁面最大熱流密度,很好地提高氣動熱性能。

(3)在進行飛行器鈍化時要權衡氣動力和氣動熱性能,本文提出的鈍化方法既有效地減小了最大熱流密度也保證了較高的升阻比,較好地兼顧了類乘波構型的氣動力和氣動熱性能。

猜你喜歡
設計
二十四節氣在平面廣告設計中的應用
河北畫報(2020年8期)2020-10-27 02:54:06
何為設計的守護之道?
現代裝飾(2020年7期)2020-07-27 01:27:42
《豐收的喜悅展示設計》
流行色(2020年1期)2020-04-28 11:16:38
基于PWM的伺服控制系統設計
電子制作(2019年19期)2019-11-23 08:41:36
基于89C52的32只三色LED搖搖棒設計
電子制作(2019年15期)2019-08-27 01:11:50
基于ICL8038的波形發生器仿真設計
電子制作(2019年7期)2019-04-25 13:18:16
瞞天過海——仿生設計萌到家
藝術啟蒙(2018年7期)2018-08-23 09:14:18
設計秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設計叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
從平面設計到“設計健康”
商周刊(2017年26期)2017-04-25 08:13:04
主站蜘蛛池模板: 激情午夜婷婷| 欧美一级99在线观看国产| 波多野结衣在线se| 四虎亚洲精品| 91人妻日韩人妻无码专区精品| 精品国产免费第一区二区三区日韩| 中文字幕在线视频免费| 一区二区自拍| 茄子视频毛片免费观看| 精品一区二区三区视频免费观看| 欧美日本在线| 99久久精品久久久久久婷婷| 亚洲欧美在线看片AI| P尤物久久99国产综合精品| AV色爱天堂网| 色天天综合| 日韩一级二级三级| 91福利免费视频| 欧美另类图片视频无弹跳第一页| 国产aaaaa一级毛片| 国产精品一区二区不卡的视频| 日本亚洲最大的色成网站www| 中文字幕有乳无码| 欧美一级在线看| 国产精品自在自线免费观看| 久久人午夜亚洲精品无码区| 小13箩利洗澡无码视频免费网站| 欧美成人二区| 日本日韩欧美| 亚洲天堂区| 国产va免费精品观看| 国产一级毛片网站| 亚洲成年网站在线观看| 日本三级欧美三级| 国产精品太粉嫩高中在线观看| 亚洲日韩精品无码专区| 国产第二十一页| 国产精品第一区| 亚洲国产成人久久精品软件 | 制服丝袜一区| 亚洲天堂视频在线播放| 在线观看国产精品第一区免费 | 欧美精品亚洲二区| 免费av一区二区三区在线| 亚洲AⅤ综合在线欧美一区| 香港一级毛片免费看| 亚洲成a人片77777在线播放| 日韩视频福利| 伊人久热这里只有精品视频99| 无码免费的亚洲视频| 国产精品吹潮在线观看中文| 国产原创第一页在线观看| 久久青青草原亚洲av无码| 国产产在线精品亚洲aavv| 日韩欧美视频第一区在线观看| 日本欧美在线观看| 亚洲天堂免费| 亚洲天堂网在线视频| 色噜噜狠狠狠综合曰曰曰| 一级毛片在线直接观看| 91精品视频网站| 日韩大片免费观看视频播放| 色综合天天视频在线观看| 999国产精品| 夜夜操国产| 一区二区三区毛片无码| 国产午夜福利在线小视频| 国产美女免费| 国产福利一区在线| 视频一本大道香蕉久在线播放| 五月婷婷中文字幕| 国产呦视频免费视频在线观看| 国内老司机精品视频在线播出| 最新精品久久精品| 亚洲无码电影| 亚洲小视频网站| 女人av社区男人的天堂| 亚洲日韩精品综合在线一区二区 | 久久精品欧美一区二区| 成人精品在线观看| 真实国产精品vr专区| 免费人成在线观看成人片|