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

振蕩激波作用下受熱壁板主共振特性分析

2022-05-16 12:04:16葉柳青葉正寅
振動與沖擊 2022年9期
關鍵詞:振動

葉柳青, 葉正寅, 洪 正, 葉 坤

(西北工業大學 航空學院,西安 710072)

吸氣式超聲速/高超聲速飛行器因具有較高的潛在軍事與民用價值,已經成為了各國在航空航天領域投入巨大人力物力來研發的熱點和關鍵性技術之一[1-3]。而沖壓/超燃沖壓發動機[4-7]作為此類飛行器推進系統的核心,在一定程度上決定了飛行器最終能達到的技術高度。沖壓/超燃沖壓發動機是一類結構簡單的吸氣式發動機,它們直接利用空氣作為燃料氧化劑,使得飛行器具有連續可控的推力,具備更高的飛行靈活性和經濟性。在沖壓/超燃沖壓發動機內部存在激波、膨脹波等復雜的波系[8-10],而滿足結構質量和主動冷卻的設計要求需要采用薄壁結構, 因此激波主導流動中彈性壁板的動力學問題是吸氣式高速飛行器結構設計與優化的重要基礎科學問題。

近二十年來,許多學者針對激波主導流動中壁板的氣動彈性問題展開研究,取得了一系列研究成果。通過Von-Karman大變形板彎理論考慮壁板幾何非線性,并且根據CFD(computational fluid dynamics)代理模型和熱傳導理論,Miller等[11]建立了基于流固熱的氣動彈性分析模型,該模型可以對激波作用下的壁板氣動彈性響應進行長時間的計算。假定激波位置固定在二維壁板的中點不變,Visbal[12-13]通過Euler和Navier-Stokes方程來分別求解變形壁板上的無黏和黏性流場,利用雙向流固耦合的數值方法研究了二維彈性壁板在激波作用下的動力學特性,結果表明當激波強度較小時,氣動彈性穩定性顯著提高;而當激波強度較大時,臨界顫振動壓遠小于沒有激波作用時壁板的臨界顫振動壓,穩定性明顯降低。Boyer等[14]將Visbal的二維氣彈模型拓展到三維,研究了三維無黏流場中在激波的作用下彈性壁板的動力學響應,結果表明強激波的作用會增大壁板振動響應的幅值和頻率,而弱激波的作用將顯著增大臨界顫振動壓,即提高壁板的氣動彈性穩定性。

Brouwer[15]提出了可用當地活塞流理論來預測激波作用下彈性壁板表面的非定常氣動力。基于Von-Karman大變形板彎理論和當地活塞流理論,Ye等[16-17]建立了激波作用在二維壁板的任意位置時受熱壁板的氣動彈性穩定的理論分析模型,推導出了壁板發生氣動彈性失穩的邊界條件,并分析了激波強度、激波沖擊位置對臨界顫振動壓、振動響應幅值和頻率的影響。結果表明,隨著激波沖擊位置從壁板一端向另一端移動變化時,臨界顫振動壓、振動幅值和頻率都不是單調變化的。采用雙向流固耦合算法,李映坤等研究了無黏流場中斜激波沖擊作用下二維曲壁板的氣動彈性響應特性,著重分析了彎曲高度和動壓對系統振動特性的影響。結果顯示較小的彎曲高度可降低壁板顫振臨界動壓,而當彎曲高度進一步增大后,臨界顫振動壓迅速提高,并且準周期無規則運動狀態被激發出來。但是,這些研究都假定激波位置固定不變,目前僅有少量的文獻考慮激波位置隨時間發生改變。Dennis[18]采用風洞試驗對在快速移動的激波作用下彈性壁板的氣動彈性特性進行了初步的研究,結果表明,與激波位置固定不變的情況相比,彈性壁板在快速移動的激波作用下振動幅值顯著提高。

國內外學者對板的非線性共振問題已開展了較多的研究工作。在20世紀90年代,Okajima等[19]對懸臂矩形板的共振頻率與曲率和厚度變化的關系進行了研究。在均勻分布的載荷作用下,Du等[20]對阻尼夾層圓板的非線性超諧共振問題進行了研究。基于Von-Karman 大變形理論,Xue等[21]建立了磁電彈性薄板的非線性無阻尼強迫振動的數學模型,并采用改進的L-P(lindstedt-poincare)法對問題進行求解,詳細地分析了板的厚度、外激勵力、壓電材料、壓磁材料以及磁電彈材料等對系統主共振區間、彈簧剛度特性和振幅跳躍現象的影響。基于多尺度法,張小廣等[22-24]對矩形板和圓板的主共振問題進行了大量的研究。在矩形板方面,張小廣等對四邊固支約束的陶瓷-金屬材料功能梯度矩形板的主共振問題進行了理論研究。在圓板方面,Hu等研究了靜磁場作用下導電旋轉圓板的主共振問題,并詳細分析了磁感性強度、轉速以及靜載等參數對頻率和振幅的影響。近年來,馬冰冰等[24]對在常磁場引起的靜載荷作用下鐵磁圓板的主共振問題進行了研究,并分析了不同電磁參量對共振振幅的影響。采用一種簡化的氣動模型[25]對作用于板上的氣動載荷進行描述,李文強等[26]對氣動載荷作用下旋轉運動導電圓板的主共振問題進行了研究。以上關于壁板主共振問題的研究工作中并沒有考慮激波的作用。

實際上,沖壓/超燃發動機內部的激波串往往存在振蕩特征[27-29],這種激波串振蕩會給發動機壁面帶來嚴酷的壓力脈動載荷。就作者目前所知,比較缺乏對振蕩激波作用下壁板的主共振特性的研究。本文基于Von-Karman 大變形理論和當地活塞流理論,采用 Galerkin 方法建立了振蕩激波作用下壁板振動的動力學模型,通過龍格-庫塔法對非線性動力學方程進行數值計算獲得系統非線性振動響應,詳細分析了振蕩激波強度、激波振蕩幅值、溫度、激波振蕩的中心位置、來流馬赫數對系統幅頻響應特性的影響,尤其是振動突跳與雙穩態特性。

1 振蕩激波作用下受熱壁板的非線性動力學方程

1.1 控制方程

當三維平壁板在流向上的尺寸遠小于另一個方向的尺寸時,它可以簡化為二維平壁板模型,這里僅考慮二維各向同性材料的壁板。振蕩激波流場中的二維平壁板,邊界條件為兩端簡支,如圖1所示。壁板的長度為l,厚度為h,其中壁板厚度遠小于壁板長度。壁板的上表面受振蕩激波流場中非定常氣動載荷的作用,入射斜激波前的氣流密度、速度和馬赫數分別為ρu,l,Uu,l,Mau,l,反射斜激波后的氣流密度、速度和馬赫數分別是ρu,r,Uu,r,Mau,r。壁板下表面受空腔壓力的作用,空腔壓力為pd。基于Hamilton原理和Von-Karman大變形板彎理論,壁板在振蕩激波作用下的氣動彈性方程[30]為

(1)

1.2 熱應力

在飛行器高馬赫數(一般大于2.2)飛行時,氣動加熱的影響一般不容忽視,還應該考慮由于氣動加熱產生的熱應力[31-33]。假設壁板受熱達到穩態并且溫度均勻分布,設T0為初始溫度,由溫升ΔT=T-T0引起的面內熱應力為

(2)

式中,α為材料的熱膨脹系數。

1.3 氣動力

在分析振蕩激波作用下受熱壁板的主共振特性中,能夠有效地預測激波振蕩流場中的非定常氣動載荷是至關重要的。活塞理論在20世紀50年代被Lighthill[34]和Ashely[35]等提出后,由于具有簡潔的表達式,且具有較好的實用性,被廣泛應用于超聲速、高超聲速非定常氣動力的計算。該理論假設當馬赫數(Ma)?1時,機翼所產生的擾動沿翼面法向傳播,而機翼表面各點間的相互影響很小,就像氣缸中活塞所產生的的擾動傳播一樣。在等熵假設條件下,通過動量定理和等熵關系式可以得到經典活塞理論計算公式為

(3)

式中:P為翼面表面任意一點的壓力;P∞為來流壓強;a∞為來流音速;γ為空氣比熱比;νn為翼面法向速度,具體表達式如下

(4)

式中,U∞為來流速度。

當來流馬赫數較小時,一般采用修正的線性活塞理論[36]

(5)

Brouwer等已經證實了當地活塞流理論預測激波主導流場中氣動力的可行性。在當地活塞流理論中,相應的自由來流流動參數用當地參數來代替。本文研究振蕩激波作用下彈性壁板的主共振特性,主要考慮結構非線性,因此這里采用當地1階活塞流理論,得到壁板上表面的氣動力如下。

入射斜激波前的氣動載荷為

(6)

反射斜激波后的氣動載荷為

(7)

(8)

反射斜激波后壁板的上下壓差為

(9)

1.4 無量綱化

引入無量綱參數

將上述無量綱參數代入式(1),則無量綱化后的控制方程為

(10)

由于壁板兩端簡支,相應的無量綱邊界條件為

(11)

實際上,本文將壁板兩端的約束簡化為固定鉸支座,如圖1所示。因此除了式(11)中所表現的邊界條件(即壁板的兩端可以轉動,并且沿著垂直方向不可以移動)以外,壁板兩端沿著水平方向也是不可以移動的。

圖1 受振蕩激波作用的二維壁板

1.5 Galerkin離散

采用Galerkin方法對式(10)進行離散,將位移函數展開成各階諧波模態的疊加。對于簡支邊界條件,位移函數為

(12)

將式(12)代入到式(10),對方程各項同時乘以sin(jπx)(j=1,2,…,∞),并沿著板長積分。

這里假設壁板上表面的激波將以壁板上某一點為中心作簡諧振動,入射斜激波沖擊在壁板上的位置記為ξk(ξk為無量綱量,ξk=xk/l),則ξk=ξ0+Acos(ωτ),ξ0為中心位置,A為激波振蕩的幅值,ω為激波振蕩的頻率。由于振蕩激波的存在,壁板左右兩端壓差不同,采用Galerkin方法對氣動力項處理如下

(13)

(14)

其中

Q=Q1+Q2+Q3+Q4+Q5,

2 模型驗證

關于受熱壁板在振蕩激波作用下的主共振問題,在第1章中已經建立了相應的動力學分析模型。本文采用4階Runge-Kutta法對此動力學模型(即式(14))進行數值積分,取無量綱時間步長為Δτ=0.001,初值為W0(ξ)=1,觀察壁板上順氣流75%長度(ξ=0.75)處的時間歷程響應,分析系統的幅頻響應特性,從而探討振蕩激波作用下受熱壁板的主共振特性。

表1 不同激波角下的靜壓比和動壓比

2.1 數值方法驗證

圖2 本文結果與Dowell的研究結果比較

圖3 本文結果與Ye等的研究結果比較

2.2 模態收斂性分析

為了選擇合適的模態數目求解振蕩激波作用下受熱壁板的非線性氣動彈性響應,有必要研究不同模態數目對響應結果的影響,即:收斂性分析。選取激波角為σ=18°,假定激波圍繞壁板中點作簡諧振動,即ξ0=0.500,激波振蕩幅值為A=0.100,振蕩頻率為ω=10,圖4顯示模態數目分別取N=2、4、6時,振動幅值隨動壓的變化。其中:WA為彈性壁板的振動幅值。從圖4可知,用2階諧波模型(N=2)算得的幅值與用4、6階諧波模型(N=4、6)算得的幅值基本一致。由此可見,當來流動壓較小時,受熱壁板在振蕩激波作用下作強迫振動,描述這種情況下的壁板變形至少取2階諧波模態。為保證計算的精度,在下面所有的分析中,均采用6階諧波模型來進行計算。

圖4 模態數目對振動幅值的影響

3 振蕩激波作用下受熱壁板主共振特性分析

采用4階Runge-Kutta法對系統動力學方程進行數值計算,得到系統的幅頻響應特性曲線,探討振蕩激波作用下受熱壁板的主共振特性,尤其是振動突跳與雙穩態特性。其中,幅頻響應特性曲線分別采用升頻掃描和降頻掃描計算[37]獲得,即采用前一個外激勵頻率激勵下響應的穩態解作為下一個外激勵頻率計算的初始條件進行數值計算。本文重點關注主共振區,因此主要截取壁板1階固有頻率附近的幅頻響應曲線進行分析。壁板的第j階無量綱固有頻率[38]為ωj=(jπ)2。影響振蕩激波作用下受熱壁板主共振特性的因素有很多,這里分別詳細地分析了振蕩激波強度、激波振蕩幅值、溫度、激波振蕩的中心位置、來流馬赫數以及下壁面背壓對壁板主共振特性的影響規律。

3.1 振蕩激波強度對受熱壁板主共振特性的影響

首先分析振蕩激波強度的影響,固定來流馬赫數為Ma=3.5,來流動壓為λu,l=30,選取不同的激波角σ為17°、18°、20°、22°、24°,由斜激波關系式計算得到的入射斜激波前與反射斜激波后的無量綱動壓比和無量綱靜壓比見表1。假定激波圍繞著壁板中點作簡諧振蕩,即ξ0=0.500,激波振蕩的幅值為A=0.010。溫升為ΔT/ΔTcr=0.1。采用4階Runge-Kutta法求解系統動力學方程,得到不同振蕩激波強度下受熱壁板的幅頻響應特性曲線,如圖5所示。

其中,圖5實線為激波振蕩頻率不斷增大的正向掃頻曲線,虛線為激波振蕩頻率不斷減小的反向掃頻曲線,后文幅頻響應曲線中實線和虛線的定義與此處一致,之后不再贅述。對于同一種算例來說,其升頻掃描得到的幅頻響應特性曲線與降頻掃描得到的幅頻響應特性曲線在兩邊的部分(即遠離共振區部分)應該是重合的。根據這個特點可從圖中區分各條曲線。

由圖5可知,當激波振蕩頻率接近主共振頻率時,系統振動幅值急劇增大,同時,在共振峰處,正向掃頻和反向掃頻過程均存在明顯的振動突跳和雙穩態現象,在升頻過程突跳點A與降頻過程突跳點B之間,形成雙穩態區間。從圖5(a),隨著振蕩激波的強度不斷增大,系統共振峰幅值不斷增大,主共振幅頻響應曲線不斷整體右移,主共振頻率不斷增大。對于弱激波而言,σ=17°,系統不存在振動突跳和雙穩態現象;當激波強度增大即σ=18°時,系統存在明顯的振動突跳和雙穩態現象,共振峰處的雙穩態區間對應的頻率為ω∈[10.3,13.0];當σ=20°時,升頻跳躍點較之前右移,而降頻跳躍點大幅右移,使得雙穩態區間變窄,此時共振峰處的雙穩態區間對應的頻率為ω∈[14.2,15.1];繼續增大激波強度,當σ=22°,主共振幅頻響應曲線較之前整體右移,但振動跳躍和雙穩態現象再次消失;當進一步增大激波強度,考慮強激波σ=24°時,系統出現顯著的振動突跳和雙穩態現象,不僅共振峰值和主共振頻率顯著提高,雙穩態區間寬度大幅變寬,此時共振峰處的雙穩態區間對應的頻率為ω∈[19.9,25.1]。

(a)

3.2 激波振蕩幅值對受熱壁板主共振特性的影響

分別取激波振蕩幅值為A=0.005、0.010、0.030、0.050、0.100、0.200,來流馬赫數為Ma=3.5,激波角為σ=18°,溫升為ΔT/ΔTcr=0.10,中心位置為ξ0=0.50,繪制系統的幅頻響應特性曲線如圖6所示。

由圖6可知,隨著激波振蕩幅值的增大,系統共振峰值不斷增大,并且發生上升跳躍和下降跳躍的頻率均隨著激波振蕩幅值的增大而增大。當激波振蕩幅值較小時,即A=0.005, 系統不存在振動突跳和雙穩態現象;當激波振蕩幅值增大時,即A=0.01時,系統出現明顯的振動突跳和雙穩態現象,共振峰處的雙穩態區間對應的頻率為ω∈[11.1,13.0];當A=0.03時,升頻跳躍點和降頻跳躍點都像右移,雙穩態區間變寬,此時共振峰處的雙穩態區間對應的頻率為ω∈[13.4,16.1];當激波振蕩幅值進一步增大時,即A=0.050、0.10,從圖6還可知,隨著激波振蕩幅值的增大,升跳躍點和降頻跳躍點不斷右移,但雙穩態區間不斷減小直至消失。當A=0.050時,共振峰處的雙穩態區間對應的頻率為ω∈[16.0,16.9];當A=0.100時,振動突跳和雙穩態現象消失。而當繼續增大激波振動幅值,即A=0.200時,從圖中可以看出,系統出現顯著的振動突跳和雙穩態現象,并且雙穩態區間寬度將急劇變寬,此時共振峰處的雙穩態區間對應的頻率為ω∈[20.1,29.4]。

圖6 不同激波振蕩幅值下受熱壁板的幅頻響應曲線

3.3 溫度對振蕩激波作用下受熱壁板主共振特性的影響

分別取無量綱溫升為ΔT/ΔTcr=0、0.1、0.5、1.0,來流馬赫數為Ma=3.5,激波角為σ=18°,中心位置為ξ0=0.500,激波振蕩幅值為A=0.010,繪制系統的幅頻響應特性曲線如圖7所示。

從圖7可知,隨著溫升的不斷增大,系統共振峰值不斷增大。當不考慮溫升時,即ΔT/ΔTcr=0,系統不存在振動突跳和雙穩態現象;當考慮溫升時,即ΔT/ΔTcr=0.1、0.5、1.0時,系統存在明顯的振動突跳和雙穩態現象;從圖7還可知,雙穩態區間的寬度隨著溫升的增大幾乎不改變,發生上升跳躍和下降跳躍的頻率隨著溫度的不斷增大而減小,系統雙穩態區間位置不斷向左移,即系統主共振區對應的頻率逐漸減小,這說明了提高溫度對系統產生“剛度弱化效應”。

(a)

3.4 激波振蕩的中心位置對受熱壁板主共振特性的影響

分別取激波振蕩的無量綱中心位置為ξ0=0.100、0.300、0.45、0.500、0.700、0.900,來流馬赫數為Ma=3.5,激波角為σ=18°,溫升為ΔT/ΔTcr=0.1,激波振蕩幅值為A=0.010,繪制系統的幅頻響應特性曲線如圖8所示。

圖8 不同中心位置下受熱壁板的幅頻響應曲線

從圖8可知,激波振蕩的中心位置對系統共振峰值有很大的影響,當激波圍繞壁板中點振蕩時,即ξ0=0.50時,系統共振峰值最大。隨著中心位置向兩端移動時,系統共振峰值不斷減小。相反地,隨著中心位置向兩端移動時,主共振頻率不斷增大。當激波圍繞壁板中點振蕩時,即ξ0=0.50時,系統出現明顯的振動突跳和雙穩態現象;而當中心位置移動到ξ0=0.45時,雙穩態區間的寬度急劇變窄;當中心位置移動到ξ0=0.1、0.9、0.3、0.7時,振動突跳和雙穩態現象消失。由此可見,將激波振蕩中心位置往兩端移動時可有效地抑制振動突跳和雙穩態現象。

3.5 來流馬赫數對受熱壁板主共振特性的影響

分別取來流馬赫數為Ma=3.0、3.5、3.8、4.0、4.5,來流動壓為λu,l=30,激波角為σ=20°,溫升為ΔT/ΔTcr=0.1,中心位置為ξ0=0.500,激波振蕩幅值為A=0.01,繪制系統的幅頻響應特性曲線如圖9所示。

由圖9可知,隨著來流馬赫數不斷增大,系統共振峰幅值不斷增大,主共振幅頻響應曲線不斷整體右移,主共振頻率不斷增大。從圖9可知,來流馬赫數對受熱壁板主共振特性的影響與振蕩激波強度、激波振蕩幅值對受熱壁板主共振特性的影響類似。當來流馬赫數較小時,Ma=3.0,系統不存在振動突跳和雙穩態現象;當馬赫數增大,即Ma=3.5,系統存在明顯的振動突跳和雙穩態現象;當來流馬赫數增大,升頻跳躍點較之前右移,而降頻跳躍點大幅右移,雙穩態區間寬度不斷變窄直至消失;而來流馬赫數較大時,即Ma=4.5,系統出現顯著的振動突跳和雙穩態現象,不僅共振峰值和主共振頻率顯著提高,雙穩態區間寬度大幅變寬。

圖9 不同來流馬赫數下受熱壁板的幅頻響應曲線

4 結 論

本文基于Von-Karman 大變形理論和當地一階活塞流理論, 建立了振蕩激波作用下受熱壁板振動的動力學模型,通過4階龍格-庫塔法對非線性動力學模型進行求解,得到了振蕩激波作用下受熱壁板的幅頻特性響應曲線,并詳細討論了振蕩激波強度、激波振蕩幅值、溫度、激波振蕩的中心位置、來流馬赫數對系統振動突跳和雙穩態特性的影響。主要結論如下:

(1)受熱壁板在振蕩激波作用下存在振動突跳和雙穩態現象等典型的非線性動力學行為,并且系統的幅頻響應特性曲線在雙穩態區表現出明顯的“硬特性”。

(2)隨著振蕩激波強度、激波振蕩幅值和來流馬赫數增大,系統共振峰值不斷單調增大,而雙穩態區間由不存在到逐漸變大,再變小直至消失,然后再急劇變大。

(3)溫度增大會對系統產生“剛度弱化效應”,使主共振幅頻響應曲線整體左移,但對系統的雙穩態區間影響較小。隨著溫度的增大,系統共振峰值不斷增大。

(4)當激波圍繞壁板中點振蕩時,系統共振峰值最大。隨著中心位置向兩端移動時,系統共振峰值不斷減小,而主共振頻率不斷增大。將激波振蕩中心位置往兩端移動時可有效地抑制振動突跳和雙穩態現象。本文主要是從數值的角度研究了振蕩激波作用下受熱壁板的主共振特性,在我們下一步的研究工作中,將從試驗的角度研究受熱壁板在振蕩激波流場中的主共振特性。

猜你喜歡
振動
振動的思考
科學大眾(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
帶有強迫項的高階差分方程解的振動性
主站蜘蛛池模板: 国产在线专区| 国产女主播一区| 久久夜夜视频| 国产97视频在线观看| 手机永久AV在线播放| 51国产偷自视频区视频手机观看 | 黄色在线网| 亚洲第一成人在线| 波多野结衣中文字幕一区二区| 国产精品视频999| 亚洲成人在线网| 亚洲国产欧美国产综合久久 | 国产玖玖视频| 色婷婷在线播放| 伊人久久久久久久久久| 欧美成在线视频| 亚洲系列中文字幕一区二区| 亚洲一道AV无码午夜福利| 国产视频a| 麻豆a级片| 亚洲精品你懂的| 免费看美女自慰的网站| 中文字幕资源站| 中文天堂在线视频| 无码精品国产VA在线观看DVD| 亚洲欧洲日韩国产综合在线二区| 国产一级精品毛片基地| 国产激爽大片高清在线观看| 国产女人在线| a级高清毛片| 色综合天天视频在线观看| 亚洲嫩模喷白浆| 91香蕉国产亚洲一二三区| 国产永久在线观看| 色综合久久无码网| 国产亚洲视频免费播放| 午夜一级做a爰片久久毛片| 亚洲综合片| 国产高清无码麻豆精品| 色哟哟国产精品| 一级毛片在线免费视频| 欧美亚洲国产视频| 精品欧美一区二区三区久久久| 国产综合在线观看视频| 亚洲色大成网站www国产| 国产麻豆福利av在线播放| 久久综合激情网| 四虎国产成人免费观看| 欧美视频在线不卡| 精品午夜国产福利观看| 黄色网址手机国内免费在线观看| 国产中文在线亚洲精品官网| 日本黄色不卡视频| 欧美高清日韩| 爽爽影院十八禁在线观看| 东京热av无码电影一区二区| 久久久久国产精品熟女影院| 色亚洲激情综合精品无码视频| 亚洲成AV人手机在线观看网站| 在线观看国产精美视频| 黄色不卡视频| 亚洲第一区在线| 99久久性生片| 欧美三级日韩三级| 国产网站免费观看| 国产丝袜无码一区二区视频| 一级高清毛片免费a级高清毛片| 国产女人在线视频| 国产色伊人| 18禁影院亚洲专区| AV在线麻免费观看网站 | 免费xxxxx在线观看网站| 丰满的少妇人妻无码区| 亚洲欧美h| 国产aⅴ无码专区亚洲av综合网| 亚洲嫩模喷白浆| 国产丝袜第一页| 亚洲人成网7777777国产| 国产精品亚洲а∨天堂免下载| 欧美v在线| 国产va在线| 国产伦片中文免费观看|