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

超聲速火箭橇流動特征和氣動力激勵振動分析

2022-11-08 09:10:44余元元王方元王彬孫建紅許常悅孫智
西北工業大學學報 2022年5期
關鍵詞:振動

余元元, 王方元, 王彬, 孫建紅,, 許常悅, 孫智

(1.南京航空航天大學 飛行器環境控制與生命保障工業與信息化部重點實驗室, 江蘇 南京 210016;2.南京航空航天大學 民航應急科學與技術重點實驗室, 江蘇 南京 211106)

在各種先進技術的推動下,超聲速和高超聲速飛行器設計受到了國內外研究機構的重點關注。例如,美國國防部21世紀初開始實施的NHFRP計劃,以及2019年發布的“高超聲速應用研究國家聯盟”項目[1]。開展高速飛行器設計研究需要開展大量的試驗以驗證新的理論與技術。對于飛行器的空氣動力學試驗而言,常見的試驗方式有風洞試驗、真實飛行試驗和地面滑行試驗等。由于受到空間限制,風洞僅能開展縮比尺寸高速飛行器試驗。真實尺寸高速飛行器的空中飛行試驗成本較高,不適合在設計階段開展。

火箭橇試驗系統屬于一種地面滑行試驗裝備。它以火箭發動機為動力裝置,推動裝有被試件的試驗平臺沿著高精度軌道高速滑行[2]。火箭橇試驗系統兼具風洞試驗和自由飛行試驗的優點。與風洞試驗相比,火箭橇試驗系統可以開展全尺寸試驗而非縮比模型,獲取的試驗數據更加真實可信。目前學者對于風洞流場品質的改善和試驗精度的提高已有過大量研究[3-5],而對于提高火箭橇試驗精度的研究較少。相比于自由飛行試驗,火箭橇試驗系統的整體研制周期更短,試驗效費比更高。因此,火箭橇試驗系統具備較好的實用價值和發展前景,可應用于機載設備的性能試驗、空氣動力學試驗等研究。

美國是最早開展超聲速火箭橇研究并著手建設的國家,設計并建造了二十多條用于火箭橇試驗的軌道。其中,最為著名的軌道為霍洛曼高速測試軌道[6]。美國以此為依托進行了多項試驗,積累了大量的數據,研究人員也針對這些數據提出了許多火箭橇系統理論。俄羅斯、英國及法國都相繼建成了自己的火箭橇試驗系統,主要用于彈射與導彈試驗[7]。

我國火箭橇試驗系統發展較晚,直至1997年才于湖北省建成第一條火箭橇滑軌。隨著我國襄北高精度火箭橇滑軌擴建工程的開展,滑軌長度已達六千余米,足以開展超聲速火箭橇試驗[8]。然而,大量的試驗數據表明,隨著運行速度的增加,火箭橇在滑軌上產生的不利振動會逐步加劇。該試驗場上的多次失敗試驗結果表明,這種不利振動會威脅到火箭橇試驗的安全[9]。火箭橇的振動來源主要涉及橇/軌接觸、火箭發動機推力波動以及高速運動帶來的氣動力致振。其中,氣動力致振是指火箭橇頭激波與地面、軌道與扣件發生碰撞反射,對橇體產生非線性的氣動力激勵振動。反射激波與邊界層相互作用會形成新的渦流[10],在一定程度上破壞理想試驗環境。此外,相較于旋成體部分,滑靴、支撐件等類方柱體結構會發生更為明顯的振動響應[11]。

針對火箭橇的氣動力問題,國內外學者已經開展了一些相關的試驗和數值研究工作。Krupovage等[12]對3種雙軌超聲速火箭橇進行了多次試驗,發現氣動力是造成拉力改變的重要原因。Rigali等[13]開展了單軌超聲速火箭橇試驗,結果表明火箭橇地面試驗可以較好地還原低空飛行器的飛行環境。Zhang等[14-16]分析了火箭橇的數值邊界條件,并采用結構有限元方法對某型雙軌火箭橇進行數值模擬,計算結果與試驗結果吻合較好。張立乾等[17]對超聲速單軌火箭橇運行過程中的氣動變化進行了數值模擬,發現氣動力的非恒定性、火箭發動機所產生推力的偏心度都會引起橇身的振動。黨峰等[18]對某火箭橇一體化設計結構進行數值模擬,并討論了跨聲速工況時氣動外形對流場特性的影響。Lofthouse等[19]利用CFD軟件對霍洛曼高速測試軌道的窄軌火箭橇在空氣和氦氣環境下進行三維數值模擬。結果表明,垂直楔形滑橇產生的激波會沖擊到火箭橇體,激波產生的高溫高壓停滯區域會降低火箭橇性能,甚至會損壞火箭橇滑行試驗。Lumb等[20]設計了一種新的火箭橇試驗,并利用CFD進行效果測試。通過對運行軌道進行預熱,可以把火箭橇運行速度由Ma3提高到Ma5。Strike等[21]對位于固定在軌道上的火箭橇進行了來流Ma為2至5的風洞試驗,并改變了離地高度、滾轉角等參數,獲得了滑塊距離、來流馬赫數等參數對火箭橇氣動載荷的影響。

由于超聲速火箭橇常采用尖拱旋成體構型,橇體頭錐處形成的斜激波與地面、軌道和軌道扣件發生復雜的碰撞、反射以及相互作用。為了深入認識這些激波結構的特征及產生的氣動力影響,有必要逐步剝離反射物體(地面、軌道及其扣件)的影響,進而探索氣動力致振的機理。基于這些研究,可以為超聲速火箭橇的氣動減振設計提供理論依據,這就是本文的研究目的。

1 數值計算方法

1.1 控制方程組

超聲速火箭橇的滑行過程屬于可壓縮湍流問題。因此,超聲速火箭橇流場的控制方程為三維可壓縮Navier-Stokes方程。在笛卡爾坐標系中,包含連續性方程、動量方程和能量守恒方程的控制方程組可以寫作如下守恒形式:

式中:Sij=0.5(?ui/?xj+?uj/?xi)為應變率張量;T為靜溫;γ為空氣比熱比;μ為空氣的分子動力學黏性系數,可根據Sutherland公式求出;μT為湍流黏性系數;Pr為層流Prandtl數;PrT為湍流Prandtl數,常近似取定值0.92。

1.2 湍流模型

湍流黏性系數μT需要通過建模求解。本文采用基于SST兩方程湍流模型的尺度自適應模擬(SAS)方法[22]開展湍流計算。在SST兩方程湍流模型中,湍流黏性系數μT可通過求解湍動能k和渦量密度ω的輸運方程得到

湍流黏性系數μT可以通過(10)式求出

(10)

輸運方程中的各種參數詳見文獻[23]。為了實現SAS計算,需要在渦量密度方程(9)中添加SAS源項QSAS

(11)

式中:ζ2=3.51,σφ=2/3,C=2,κ=0.41;L為與湍流長度尺度

(12)

卡門長度尺度按(13)式計算

(13)

式中:Δ=(ΔxΔyΔz)1/3,CS=0.11。

1.3 計算方法

火箭橇沿著軌道的滑行過程為單一方向的滑行,故可通過指定運動函數的動網格方法對火箭橇的運行過程進行模擬。為了規定網格變形的范圍,需要在繪制計算域網格時指定運動部分。本文采用基于鋪層算法的動網格技術[23]對火箭橇的滑行過程進行模擬。作者前期的研究表明基于SAS方法的動網格技術可以用于超聲速火箭橇滑行流場的模擬[24]。

控制方程組(1)~(3)采用基于結構網格的有限體積方法進行求解。為了捕捉超聲速火箭橇流場中的激波等強間斷結構,采用二階精度的Roe通量差分裂格式對黏性項和對流項進行離散。為了保證時間的二階精度,時間推進采用子迭代技術,單個時間步長內進行20次迭代計算。在當前模擬中,火箭橇流場中的固壁設置為無滑移無穿透的絕熱壁,地面、軌道及軌道扣件處的固壁設置為滑移壁。其他邊界處設置壓力遠場邊界條件。

1.4 計算物理模型

鑒于超聲速繞流存在頭激波的特點,為了減少頭激波引起的波阻,現有的超聲速火箭橇常采用尖拱旋成體構型。本文研究的超聲速火箭橇安裝在單軌上,橇體采用尖拱旋成體構型,如圖1a)所示。為簡化計算,數值模擬采用的火箭橇計算物理模型僅保留火箭橇試驗平臺的橇體、軌道及軌道扣件,去除下方的支撐結構及車身,如圖1b)所示。火箭橇橇體長度為2 170 mm,直徑為280 mm。軌道扣件長寬高分別為160,150和150 mm,沿軌道周期性安裝在地面,相鄰扣件間距為1 250 mm。火箭橇中軸線距離地面高度為563 mm。火箭橇流場數值模擬采用長方體計算域,計算域長、寬、高分別為28 328,4 500和4 000 mm。

圖1 超聲速火箭橇的物理模型圖

在超聲速火箭橇滑行過程中,旋成體構型的橇體頭錐處會形成附體激波。為了對比分析激波結構特征,本文還對2種特殊工況進行研究:①超聲速火箭橇橇體的無限空間內繞流;②僅包含地面的超聲速火箭橇橇體繞流。在第一種工況中,長方體計算域長、寬、高分別為41 239,26 320和26 880 mm;在第二種工況中,火箭橇中軸線距離地面高度也為563 mm,長方體計算域的長、寬、高分別為41 239,26 320和13 583 mm。

2 計算結果分析與討論

2.1 計算細節

為了便于討論分析,對研究的對象分別做如下標記:火箭橇橇體在不受限環境下的超聲速繞流記為“Case1”;火箭橇橇體在僅包含地面環境下的超聲速繞流記為“Case2”;帶軌道及軌道扣件的超聲速火箭橇滑行過程記為“Case3”。

3種研究對象模擬均采用結構網格,并且沿壁面法向對邊界層內的網格進行局部加密。Case1網格總數約為880萬;Case2網格總數約為770萬;在Case3中,由于采用基于鋪層算法的動網格技術模擬火箭橇的滑行過程,故針對運動區域和靜止區域分別構建結構網格。其中,運動區域網格總數約為170萬,靜止區域網格總數約為1 100萬。火箭橇壁面法向第一層網格間距的Δr+值約為O(10),滿足SAS方法中的近壁RANS計算要求。其中,Δr+表示利用黏性壁面單位規則化的壁面法向網格間距。所有模擬的時間步長均取定值5×10-6s,環境溫度取300 K,來流(滑行)馬赫數均取為2,來流(滑行)攻角均為零。此外,為了考察滑行速度對超聲速火箭橇氣動特征的影響,還對研究對象Case3模擬Ma2.5和Ma3運行速度下的火箭橇滑行過程。Case1的總迭代計算時間為0.05 s,Case2的總迭代計算時間為0.05 s,Case3的總迭代計算時間為0.04 s。

2.2 地面效應對超聲速火箭橇氣動特征的影響

當超聲速氣流掠過尖拱旋成體時,會在其頭部形成附體弓形斜激波[25]。激波屬于一種強間斷面,氣流在激波面前后的氣體參數會發生突然變化,滿足Rankine-Hugoniot關系。當火箭橇橇體處于不受限的超聲速氣流環境時,橇體的超聲速流場幾乎呈現對稱特征,橇體頭錐處的斜激波結構向遠方傳播直至耗散,橇體尾部會形成膨脹波,如圖2a)所示。當僅存在地面影響時,激波會在橇體下方發生規則反射,如圖2b所示。經過地面反射后的激波傳播至橇體下表面,發生激波/湍流邊界層相互作用,進而導致橇體表面局部壓力上升。由于耗散效應的存在,反射至橇體的激波經過二次碰撞后強度明顯減弱。對于超聲速火箭橇而言,向下傳播的頭激波會與軌道、軌道扣件以及地面發生碰撞。由于軌道扣件按照等間距方式安裝在火箭橇地面軌道兩側,故高速運動的頭激波會周期性地掠過軌道扣件,發生周期性碰撞。對比頭激波未作用于軌道扣件時刻和作用于軌道扣件時刻(如圖2c)~2d)所示),火箭橇橇體下壁面的壓力分布類似,這說明頭激波掠過軌道扣件時未產生反射至橇體的激波。然而,當頭激波與地面碰撞反射時,也會產生作用于橇體的反射激波。此外,當頭激波掠過軌道扣件時,會在軌道扣件后方產生局部高壓區。

圖2 火箭橇對稱截面上的瞬時靜壓分布

為了定量分析地面效應對火箭橇表面壓力分布的影響,圖3給出了沿著火箭橇上下表面中心子午線上的壓力系數Cp曲線。

圖3 沿火箭橇上下表面中心子午線上的瞬時壓力系數分布

由于頭激波的存在,火箭橇頭錐處的壓力明顯高于其他區域。當火箭橇橇體處于不受限環境時,頭錐后方的橇體表面壓力系數分布較為均勻,這與橇體表面流動未受到周圍環境影響有關。當僅存在地面的影響時,頭激波經過地面反射作用于x≈4.1D附近的橇體表面,反射激波使得該處的Cp值明顯抬升。當地面安裝有軌道及軌道扣件時,軌道碰撞反射的激波使得x≈2.2D附近的Cp值得到抬升,地面反射的激波使得x≈5.7D附近的Cp值也得到抬升。需要特別注意的是,在2.2D

反射激波在橇體表面作用的位置并非單獨的點,而是沿著橇體表面的一個曲線。因此,橇體下表面的局部壓力波動勢必導致橇體升力的變化。表1給出了不同工況下的超聲速火箭橇受力特征。

表1 超聲速火箭橇的平均受力特性

可以看出,反射激波的存在使得火箭橇產生升力,這將加劇火箭橇試驗平臺車身與軌道接觸引起的振動。此外,反射至橇體的激波還會導致火箭橇阻力略微上升。圖4給出了橇體下表面的瞬時摩擦阻力系數Cf分布。Cf的定義為

(14)

(15)

式中:τw為壁面摩擦應力;n代表橇體表面法向;ρ∞為環境空氣密度;U∞為火箭橇來流(滑行)速度。從圖4中可以看出,3個橇體的頭錐處Cf值明顯高于其他位置。根據Rankine-Hugoniot關系不難得知,這與頭錐位于頭激波波后有關。在Case2和Case3中,如圖4b)~4d)所示,由于反射激波與橇體湍流邊界層發生相互作用,在橇體下表面形成弧形作用區域。與Case1相比,這些弧形作用區域處的Cf值略有下降,而作用區域下游的Cf值則明顯上升。Zuo等[27]也發現了類似的現象。從整體而言,反射激波對橇體產生的氣動力激勵會使得橇體的表面摩擦阻力增大,這與表1得出的阻力變化趨勢一致。

圖4 火箭橇下表面上的瞬時摩擦因數分布

激波與湍流邊界層相互作用在影響壁面摩擦因數的同時,也對近壁處的渦量生成產生影響。壁面渦通量(ωBVF)可以反映壁面的渦動力學特征[28]。ωBVF的定義為

(16)

當雷諾數較高時,ωBVF可近似改寫成

(17)

公式(16)和(17)中:ω為偽渦矢量,n表示壁面法向矢量。從公式(17)中可以看出,ωBVF可以反映壁面的雙重動力學影響[29]:一是外部流場對壁面的脹壓效應;二是壁面的剪切效應,也即近壁附近渦量產生的機制。

圖5給出了火箭橇下表面的瞬時壁面渦通量幅值分布。在3種工況中,橇體頭錐處均有明顯散布的|ωBVF|峰值,這與橇體頭激波的脹壓影響有關。在Case1工況中,橇體頭錐下游的|ωBVF|值幾乎為零;在Case2工況中,地面反射激波對橇體中部產生脹壓影響,進而導致該位置處出現|ωBVF|峰值;在Case3工況中,橇體上出現多處|ωBVF|峰值,這與頭激波與軌道、軌道扣件以及地面發生碰撞產生復雜的反射激波有關。此外,頭激波與軌道扣件碰撞對橇體表面的ωBVF值分布影響較弱,這進一步印證了頭激波與軌道扣件碰撞后產生的擾動波屬于強度較弱的壓縮波。由此可見,反射激波與橇體的湍流邊界層發生相互作用是橇體近壁處渦量產生的重要原因。需要注意的是,激波的脹壓效應對跨越激波的渦量變化影響可以忽略,激波引起的渦量增加則與流體拉伸效應有關[30]。這進一步驗證了ωBVF反映的上述動力學影響機制,即剪切效應是渦量產生的重要原因。

圖5 火箭橇下表面上的瞬時壁面渦通量幅值分布

在火箭助推下,火箭橇通過加速實現超聲速滑行。由于火箭助推器推力并非完全均勻,火箭橇的超聲速滑行過程屬于變加速過程。為了對比認識滑行速度對超聲速火箭橇流動特性的影響,把超聲速火箭橇的變加速滑行過程簡化為不同滑行速度下的勻速運動。圖6給出了不同滑行馬赫數下的超聲速火箭橇速度場分布。

圖6 不同滑行馬赫數下的超聲速火箭橇速度場

可以看出,超聲速火箭橇的尾跡存在明顯的非定常特征,這與反射激波誘導橇體表面的渦量生成有關。隨著滑行馬赫數的增加,火箭橇頭錐處的弓形斜激波激波角變小。當Ma=2時,火箭橇的尾跡向地面偏斜。然而,隨著Ma的增加,火箭橇的尾跡向上抬升,且呈現明顯的上下周期擺動。經過分析可知火箭橇尾跡的周期擺動波長λ約等于軌道扣件間距。不難得知,這些擾動波的傳播速度與火箭橇的滑行速度一致。因此,這些擾動波的傳播頻率fw可以通過(18)式進行計算

(18)

式中,U∞為超聲速火箭橇滑行速度。當火箭橇的滑行馬赫數分別為2,2.5和3時,這些尾跡擾動波的波動頻率分別為554,680和816 Hz。

2.3 超聲速火箭橇的氣動力激勵振動

與列車軌道類似,火箭橇的軌道扣件也是等間距安裝在地面上。當超聲速火箭橇滑行時,頭激波會周期性地掠過軌道扣件,頭激波碰撞到軌道、軌道扣件以及地面時均會產生反射激波。這些反射激波作用在橇體的位置會發生周期性變化,對橇體產生周期性的氣動力激勵,進而導致火箭橇產生振動。

為了獲取火箭橇氣動力激勵振動的特征頻率,在火箭橇橇體下表面設置了3個探測點,如圖4d)所示。圖7給出了這幾個探測點處的壓力信號曲線。

圖7 火箭橇下表面上的壓力信號

圖中,p∞為環境壓力。當Ma=2時,3個探測點處的壓力信號存在明顯的周期性波動;當Ma=2.5時,P1探測點處的壓力信號幾乎不存在波動,這意味著此時P1點位于反射激波作用區域上游;當Ma=3時,P2點處的壓力信號也幾乎不存在波動,這說明反射激波作用的區域隨著滑行速度的增加而向下游繼續移動。由此可見,滑行速度增大,頭部激波強度增大,激波角度減小,相應的反射激波角也減小,其影響區域后移,影響區域的壓力波動頻率變高,波動值變化范圍減小。

圖8 壓力信號的功率譜密度分析

反射激波對火箭橇橇體產生的氣動力激勵振動可以通過對壓力信號作功率譜密度(PSD)分析得出。圖8為不同滑行馬赫數下3個壓力信號的PSD分析曲線,曲線峰值對應的Strouhal數(Sr)即為氣動力激勵振動的特征頻率。Sr的定義為

(19)

式中:f為頻率。從圖8中可以看出,Ma=2,2.5和3時的特征Strouhal數分別為0.224,0.220和0.225(對應的頻率分別為554,680和816 Hz)。當滑行馬赫數為2時,襄北試驗場的實測氣動力激勵振動Sr數為0.22,這表明當前計算結果具有較好的可信性。由上文分析可以知道,氣動力激勵振動的頻率與尾跡擾動波的波動頻率一致,這意味著火箭橇氣動力振動對尾跡波動具有鎖頻現象。這種鎖頻現象也常見于振翅射流射流和圓柱周期振蕩流[31-32]。此外,PSD曲線中還可以看出,反射激波引起的火箭橇氣動力激勵振動存在倍頻關系,這說明氣動力振動存在諧聲現象,這在Lamb的超聲速火箭橇試驗中得到了印證[33]。振翅射流和圓柱周期振蕩流中的鎖頻現象也往往在功率譜曲線上出現倍頻關系[31-32]。因此,有理由認為氣動力振動的諧聲現象與氣動力振動對尾跡波動的鎖頻現象密切相關。

3 結 論

本文采用基于鋪層算法的動網格技術對超聲速火箭橇的高速滑行過程進行數值模擬。為了深入認識地面效應對超聲速火箭橇氣動特性的影響,對比分析了橇體在不受限環境和僅存在地面效應下的超聲速繞流。當前工作還對不同滑行速度下的火箭橇流場拓撲特征以及氣動力激勵振動進行了研究。通過對計算結果的分析與討論,有如下結論:

1) 超聲速火箭橇的頭激波傳播至軌道和地面時發生激波反射,頭激波與軌道扣件發生碰撞產生強度較弱的壓縮波,導致火箭橇表面壓力出現不同幅度的振蕩。

2) 反射激波的存在使得超聲速火箭橇產生升力,且阻力略有增加。升力的產生將進一步加劇火箭橇車體與軌道之間的碰撞振動。此外,反射激波使得橇體表面的渦量增加,進而使得其尾跡區流動呈現非定常波動。

3) 反射激波對超聲速火箭橇產生氣動力激勵振動。通過對壓力信號的功率譜分析可知,該振動的頻率和尾跡的波動頻率一致,即氣動力振動對尾跡波動具有鎖頻現象。此外,功率譜密度分析表明氣動力激勵振動存在倍頻關系,這說明氣動力激勵振動存在與鎖頻現象密切相關的諧聲現象。

猜你喜歡
振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
某調相機振動異常診斷分析與處理
大電機技術(2022年5期)2022-11-17 08:12:48
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
This “Singing Highway”plays music
具非線性中立項的廣義Emden-Fowler微分方程的振動性
中立型Emden-Fowler微分方程的振動性
基于ANSYS的高速艇艉軸架軸系振動響應分析
船海工程(2015年4期)2016-01-05 15:53:26
主回路泵致聲振動分析
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
帶有強迫項的高階差分方程解的振動性
主站蜘蛛池模板: 第九色区aⅴ天堂久久香| 成年人国产视频| 国产特级毛片| 青青草国产在线视频| 亚洲免费毛片| 国产精品亚洲专区一区| 91蜜芽尤物福利在线观看| 国产97色在线| 最近最新中文字幕免费的一页| 国产在线日本| 国产黑丝视频在线观看| 在线观看av永久| 污污网站在线观看| 国产高清无码麻豆精品| 少妇精品网站| 91精品国产综合久久不国产大片| 中文字幕色在线| 97国内精品久久久久不卡| 亚洲性影院| 婷婷六月激情综合一区| 国产成年无码AⅤ片在线 | 无码啪啪精品天堂浪潮av| 免费无码网站| 国产超碰一区二区三区| 91最新精品视频发布页| 久久无码av三级| 国产一级妓女av网站| 欧美自慰一级看片免费| 四虎在线观看视频高清无码| 真人高潮娇喘嗯啊在线观看 | 国产精品成人AⅤ在线一二三四| 91在线视频福利| 色吊丝av中文字幕| 中日韩一区二区三区中文免费视频| 国产免费精彩视频| JIZZ亚洲国产| 超碰免费91| 亚洲午夜福利精品无码| 91尤物国产尤物福利在线| 日韩人妻精品一区| 波多野衣结在线精品二区| 在线国产你懂的| 人妻无码中文字幕第一区| 三级欧美在线| 精品视频在线观看你懂的一区 | 乱人伦99久久| 色哟哟国产成人精品| 亚洲久悠悠色悠在线播放| 日韩精品无码免费专网站| 亚洲国产精品美女| 亚洲国产欧美国产综合久久| 高清精品美女在线播放| 日本人又色又爽的视频| 不卡无码网| 亚洲欧美日韩精品专区| 亚洲三级视频在线观看| 亚洲AV一二三区无码AV蜜桃| 中文字幕亚洲综久久2021| 国产在线小视频| 99热免费在线| 亚欧成人无码AV在线播放| 亚洲欧美成人影院| 欧美啪啪视频免码| 国产人妖视频一区在线观看| 五月婷婷亚洲综合| 综合色88| 久草国产在线观看| 99久久99视频| 日韩av资源在线| 免费观看欧美性一级| 不卡视频国产| 欧美一区二区三区国产精品| 亚洲黄色视频在线观看一区| 国产免费精彩视频| 国产精品jizz在线观看软件| 亚洲美女视频一区| 亚洲精品第1页| 91成人在线观看| 亚洲永久色| 好吊色妇女免费视频免费| 亚洲国产精品国自产拍A| 色综合久久久久8天国|