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

外源激勵作用下火箭發動機輸流管路振動特性研究

2024-02-01 01:57:22宮武旗
振動與沖擊 2024年2期
關鍵詞:模態振動

蘇 勇, 何 江, 張 淼, 宮武旗

(1. 西安交通大學 能源與動力工程學院,西安 710048; 2. 中國航天科技集團有限公司 西安航天動力研究所,西安 710100)

液體火箭發動機運行過程中,受煤油一級泵形成的強烈壓力波動和發動機主體振動的影響,與一級泵出口和燃燒室相連的煤油輸送管路(如圖1所示)會出現強烈的振動。可能造成結構失效破壞、發動機燃燒不穩定等問題,嚴重影響發動機可靠性。隨著大推力發動機的日益發展,管路振動問題愈發顯著,勢必嚴重制約大推力火箭發動機的研究[1]。因此開展煤油輸送管路振動研究對提升液體火箭發動機動力學認知水平、提高系統穩定性與可靠性具有積極意義。

圖1 發動機主體及管路結構示意圖Fig.1 Engine body and pipeline structure diagram

管路系統一般會存在泵、閥門、分支管路等產生的流體波動和系統基礎振動等外部激勵,因而會誘發管路結構振動,此類管路結構及流體還會相互耦合而產生流固耦合振動。看似穩定的管路系統在輸送流體的同時,會將振動能量沿管壁和管內流體傳播到系統的各個位置,造成管路及與之相連的其他元件和精密儀器的破壞,影響系統安全和管路動力系統的正常運行[2]。因此管路系統的研究是國內外眾多學者研究的熱點[3-6]。Gao等[7]結合飛機液壓管路系統的發展趨勢和振動控制技術要求,對液壓管路系統振動及控制技術進行了詳細綜述,并對振動控制技術在工程領域的應用提出了建議。Jia等[8]試驗研究了斜管中的氣固兩相流。將氣體作用下,斜管內固體顆粒沿斜管向下的流動分為了三種流型,為類似斜管的設計提供了理論依據。Adamkowski等[9]試驗研究了固定支撐的剛度對管道振動的影響,結果表明支撐剛度增大,管路振動加速度的頻率升高,同時使振動加速度頻率下的幅值減小。Liu等[10-11]采用流固耦合方法對天然氣出站管路振動問題進行研究,結果表明管路振動是由于流體壓力頻率與管道模態頻率接近導致。通過增大管徑、增加約束等方法,有效地對管路實現了減振。另外,El-Borgi等[12-14]研究了周期性分布的局部減振器對管道特定頻寬下的振動衰減作用。基于以上研究,可初步了解管道動力學特性以及管道流體流動、振動和脈動壓力相互之間關系等問題。然而,對流體引起管道振動機理和外源激勵對管道振動影響等方面的理解仍然不足。例如,Liu等通過單個測點的傅里葉變換頻譜結果分析流體激勵管道共振問題。僅從數個點的角度分析流激振動,很難確定流體引起結構振動的機理。Wu等[15-17]僅分析了管道內部流體與結構的相互作用,卻沒有考慮離心泵引入的流體波動和基礎振動,因此未確定外源激勵對管道振動的影響。

近十幾年來,流固耦合計算已成為解決復雜振動問題的有力工具[18-22],并已被用于輸流管道的振動研究[23-25]。Jiang等[26]采用單向流固耦合方法模擬離心泵中流動引起的機械振動和噪聲,對振動結果進行了可視化,闡明了共振噪聲產生和傳播的機制。Pittard等[27]利用數值模擬研究表明,管道振動加速度與流體流量之間呈近似二次關系,管道振動加速度與湍流引起的管壁壓力波動成正比。Du等[28]應用流固耦合計算分析了大推力液體火箭發動機燃燒室在強振動試驗條件下的結構振動,確定了燃燒室聲壓與振動信號的耦合特性,通過提高燃燒室結構的剛度對燃燒室聲振耦合模式進行解耦,使熱試驗時的振動加速度值降低了2/3。Zhang等[29]對存在空化流動的誘導輪進行了流致振動特性的數值研究,結果表明誘導輪中的空化引起的振動主要表現為一階彎曲,葉片的最大位移發生在葉尖前緣,從葉尖前緣到葉根逐漸減小。雖然流固耦合計算方法已被廣泛應用,然而對承高壓、受煤油一級泵形成的強烈壓力波動和發動機主體振動影響的煤油輸送管路的數值研究較少。

本文采用雙向流固耦合方法,研究了液氧煤油火箭發動機一級泵出口煤油輸送管路的振動特性。首先利用試驗得到仿真所需邊界條件,并對仿真結果的有效性進行驗證。之后分析引起管道振動的內部原因和外部原因。內部原因的研究內容是管道純流體流動的流固耦合研究;外部原因研究內容包括兩部分,分別是僅有外來流體壓力脈動激勵的流固耦合研究和包括所有外源激勵的熱試車工況流固耦合研究。綜合內外因素分析,得到管道振動原因及特性,為火箭發動機結構優化提供理論支撐。

1 數值方法和模型簡化

1.1 數值理論

1.1.1 流場計算

流場控制方程由Navier-Stokes方程、湍流輸運方程及狀態方程組成。其控制方程通用形式為

(1)

式中:Sφ為源項;U為三維坐標上的速度u,v,w;ρ為流體密度;φ為通用變量;Γφ為廣義擴散系數。

1.1.2 管道結構計算

瞬態動力學分析求解的基本運動方程為[30]

[M]u″+[C]u′+[K]u={F}

(2)

式中: [M]為質量矩陣; [C]為阻尼矩陣; [K]為剛度矩陣;u″為節點加速度向量;u′為節點速度向量;u為節點位移向量; {F}為流體激勵力。流固耦合計算時,{F}即為流固接觸面上流體的壓力脈動激勵,而模態計算時,若不考慮流體的壓力影響,則{F}為零。

1.2 計算方法

本文對流場的模擬方法選取雷諾平均方法,湍流模型為k-ε模型[31],流場壓力和速度耦合計算采用Coupled算法,動量和能量方程的離散均采用迎風格式,使用中心差分技術對擴散項進行離散。流體介質為液態煤油,計算過程中考慮為不可壓縮,且不考慮煤油在整個過程中的溫度變化。

固體結構基于有限元法求解,將連續的求解區域離散成有限個單元,在每個單元內假設位移函數形式,基于變分原理求得單元剛度矩陣,然后將單元內剛度矩陣組裝成總剛度矩陣,用于計算整體結構的頻率、模態、節點位移等,具體結構上某點的位移、內力可通過對周圍節點的插值來計算[32]。

流固耦合采用分區雙向流固耦合方式,即計算時,流體和結構分開計算,在流固接觸面上進行數據交換。流固接觸面上管道接收流體傳遞的激勵力,流體接收管道在流體激勵力作用下反饋的位移。流固耦合面上流體域和固體域滿足位移協調方程和力平衡方程

df=ds

(3)

τf=τs

(4)

式中:d為位移;τ為應力; 下標s和f分別為固體和流體。

本文基于以上理論公式和計算方法,應用ANSYS 2020版商業軟件進行雙向流固耦合計算。流場計算采用Fluent模塊,固體計算采用Transient Structural 模塊,流場和固體數據交換采用System Coupling 模塊。計算時間步長6.25×10-5s,計算時長0.11 s。

1.3 模型簡化

在管路安裝圖、發動機安裝詳圖等圖紙基礎上,結合現場調研收集的實際參數,建立管路三維模型。管路模型包含直管、彎管和波紋管三種管段。波紋管段包含多層波紋管,導流板,搖擺支撐和鎧裝裝置。管路簡圖如圖2所示,包含一個介質進口、三個出口和兩個位置的支撐。

圖2 管路結構簡圖Fig.2 Pipeline structure diagram

鑒于真實模型計算成本較高,在仿真前以剛度、質量等效原則,對模型進行適當簡化,簡化位置如圖3所示。將搖擺裝置中波紋管外層的鎧裝結構略去,波紋管8層簡化為1層,厚度減薄至原總厚度的3/4。搖擺裝置中結構之間的摩擦因數設置為0.03。在剛度計算時,設定波紋管段一端固定,在另一端分別施加軸向和徑向1 mm的位移,以計算對應的載荷力,利用載荷可得到相應的軸向剛度和彎曲剛度[33]。

圖3 模型簡化位置示意圖Fig.3 Model simplified position diagram

仿真值與試驗值對比如表1所示。軸向剛度與彎曲剛度的相對誤差分別為2.8%和2.9%,相對誤差較小表明簡化合理。鎧裝裝置和波紋管質量僅占整個波紋管段結構的質量約3.1%,占整個管路比例更小,因此質量變化的影響可以忽略。仿真驗證了簡化模型的有效性。

表1 仿真與試驗剛度值比較Tab.1 Comparison of simulation and test stiffness values

1.4 網格劃分

網格無關性結果如圖4所示在額定工況下,使用入口-出口1、入口-出口3的壓降結果對流場進行網格無關性驗證。當網格數從410萬增加到630萬時,壓降幾乎沒有變化,后三種網格壓降的最大相對誤差僅1.4%,表明410網格數已滿足網格獨立性要求。以固體網格數對前10階模態頻率的影響進行網格無關性分析。網格數大于21萬后,管道模態頻率幾乎沒有變化。最終選用的流體網格數410萬,固體網格數54萬。

圖4 網格無關性分析Fig.4 Grid independence analysis

2 結果分析

2.1 仿真有效性驗證

熱試車過程中,主控室向試車平臺發送控制信號,試車平臺產生的試驗數據則通過傳感器傳輸至主控室,進行觀測記錄,如圖5所示。試驗測試儀器采用LMS公司開發的LMS Test.Lab系統,集成了高速多通道數據采集與試驗、分析、電子報告工具等功能,可進行包括數據采集、數字信號處理、旋轉機械分析、聲學等多種試驗。試驗中脈動壓力傳感器測量范圍0.05~0.50 MPa,振動加速度傳感器量程0~1 000g。數據測量過程自發動機點火指令發出開始,測量系統同步接收時統信號開始數據采集,仿真計算所需的邊界條件均源于主控室記錄的數據。

圖5 試車平臺簡介圖Fig.5 Test platform introduction diagram

在發動機熱試車試驗過程中,記錄振動加速度監測點1和壓力監測點4處的數據,位置如圖6所示。對比了仿真與試驗結果(如圖7、圖8所示)。由于試驗與仿真選取的測點位置難以完全一致,而且仿真很難完全模擬熱試車真實情況,導致壓力脈動和振動加速度與仿真結果存在輕微差異,但是仿真得到的流場壓力和振動加速度結果主要特征與試驗是一致的,驗證了數值分析模型的有效性。

圖6 監測點及支撐位置示意圖Fig.6 Monitoring point and fixed support position diagram

圖7 壓力頻譜比較結果Fig.7 Pressure spectrum comparison results

圖8 仿真和試驗振動加速度頻譜結果比較Fig.8 Comparison of simulation and experimental vibration acceleration spectrum results

2.2 管道模態

管道邊界條件包括4個固定支撐,分別為進、出口法蘭連接處、支撐1和支撐2,固定支撐位置見圖6。在支撐和搖擺裝置處包含M8和M32兩類螺栓預緊,預緊力分別為25 000 N和25 641 N,共有16個螺栓。兩處搖擺裝置及波紋管使用簡化后的模型結構。管道材料為結構鋼,密度7 830 kg/m3,彈性模量2×1011Pa,泊松比0.3。

比較了干模態、濕模態和預應力模態三種管道模態結果,其中干模態僅包含邊界條件,不包含流體質量和流體壓力;濕模態在邊界條件的基礎上,考慮流體質量的影響,但不包含流體壓力;預應力模態在邊界條件的基礎上,將流體的穩態壓力結果施加至流固接觸面,但不考慮流體質量影響。模態計算比較結果如圖9所示。濕模態計算中,流體在管路中靜止,得到的濕模態頻率相比干模態均出現不同程度的降低,與文獻[34]結果一致。本文研究的管道結構復雜,而且包含四處不同位置的固定支撐和多處螺栓預緊,多種因素影響下管道濕模態的低階頻率相比干模態變化較小。預應力模態是將流體穩態計算得到的壓力結果施加至流固接觸面,得到的低階模態頻率結果相比干模態和濕模態變化較大。綜合分析可知,流體質量對管道低階模態結果影響較小,流體穩態壓力對管道低階模態結果影響較大,而低階模態是著重關注點,因此采用預應力模態進行分析。預應力模態結果如圖10所示。

圖9 模態結果比較圖Fig.9 Modal result comparison diagram

圖10 預應力模態結果Fig.10 Prestressed modal results

2.3 引起管道振動內部因素

流體經過彎管、波紋管等結構會出現壓力波動[35],并最終以流體激勵力的形式作用在流固接觸面。當激勵力過大或其頻率與管道固有頻率一致時,會使管路產生振動。本節分析無外源激勵時,純流體流動引起的管道振動特征,邊界條件是火箭發動機熱試車的額定工況。

2.3.1 邊界條件

流場邊界條件如表2所示,流場瞬態計算以穩態收斂的結果為初場。管道邊界條件與模態計算時一致。

表2 額定工況流場邊界Tab.2 Flow field boundary of rated condition

2.3.2 純流體流動流固耦合流場結果

沿流向設置7個壓力監測點,位置見圖6。監測點的時域和頻譜結果如圖11所示,壓力時域和頻域的幅值均沿流向增大,表明沿流向壓力波動越來越劇烈。流場頻譜結果中,幅值較大的頻率為1 153 Hz和1 837 Hz。

圖11 流固耦合流場結果Fig.11 Fluid-solid coupling flow field results

2.3.3 純流體流動流固耦合振動結果

監測3個位置的振動加速度(見圖6)。監測點傅里葉變換的頻譜結果如圖12所示。結果包括120 Hz,228 Hz,241 Hz和335 Hz的振動加速度頻率,而無明顯的高頻振動。其中228 Hz和241 Hz兩者差值百分比僅為5.4%,由結構的非線性引起[36-38],應是同一頻率。振動加速度頻率120 Hz,241 Hz和335 Hz,分別對應3階(115 Hz)、6階(242 Hz)和9階(343 Hz)模態頻率,而流場的脈動壓力特征頻率(116 Hz,232 Hz和335 Hz)與管道3階,6階和9階固有頻率相近,因此,管道振動可能是由流體和管道共振產生。流場結果中,幅值較大的壓力頻率并未激起管道振動,原因是這些頻率與管道固有頻率并不相近,且產生的激振力未達到引起管道強迫振動的程度。

圖12 振動加速度頻譜Fig.12 Spectrum of vibration acceleration

雖然流場中存在與管道固有頻率和振動加速度頻率均相近的壓力頻率,但也存在并未激起管道振動的其他與管道固有頻率相近的壓力頻率,因此揭示流激振動的原因,需作進一步分析。首先提取流固接觸面上流場和管道結構在每個時間步中的數據,然后對數據進行可視化,分析流場壓力和振動加速度的相關性。流固接觸面上流場壓力和振動加速度的均值分布,如圖13所示。振動加速度較大的位置主要集中在中間管道和波紋管處。流場壓力沿流向逐漸減小,但在彎管、波紋管等位置出現明顯壓力變化。均值分布結果中,流場壓力和振動加速度無明顯相關性。

圖13 振動加速度和流場壓力均值分布Fig.13 The mean value distribution of vibration acceleration and flow field pressure

對流固接觸面上每個節點的振動加速度和流場壓力作傅里葉變換,得到3個主要振動加速度頻率(120 Hz,241 Hz和335 Hz)和相應壓力頻率的幅值分布,如圖14所示。3個振動加速度頻率中,121 Hz和241 Hz的較大幅值集中在中間管道和波紋管處,335 Hz的較大幅值集中在波紋管和靠近出口的直角彎管處。相應流場壓力頻率的幅值分布中,121 Hz較大的幅值集中在靠近出口管道、中間管道和波紋管處,241 Hz較大的幅值集中在中間管道和波紋管處,335 Hz的幅值在下游較大。對比可以看出,121 Hz和241 Hz頻率下的振動加速度和壓力均在中間管道處有較大的幅值,335 Hz頻率下的振動加速度和壓力均在下游波紋管和彎管處有較大的幅值。因為支撐2的作用,雖然121 Hz和335 Hz的流場壓力頻率在靠近出口附近有較大幅值,但是振動加速度頻率幅值在靠近出口附近并不大。振動加速度在流向第一個波紋管處有較大幅值,與管道振型和波紋管吸收變形的特性有關[39-40]。

圖14 振動加速度和流場頻率幅值分布Fig.14 Frequency amplitude distribution of vibration acceleration and flow field pressure

由圖14可知,流場壓力與振動加速度在相應頻率下的幅值分布有較明顯的相關性。綜合均值分布結果可知,管道振動與流場壓力頻率的幅值分布相關,與壓力的均值分布無關,表明管道自激振動是由流場壓力波動引起。流場壓力與振動加速度在相同頻率下幅值分布的相關性表明了流場壓力引起管道結構共振的原因。

雖然335 Hz頻率下的壓力幅值比121 Hz和241 Hz大,但引起的振動幅值卻是最小的。相反121 Hz頻率下的壓力幅值最小,但引起的振動加速度幅值最大,表明管道低階模態更易被激發,因此應重點關注低階模態頻率的振動。

2.4 引起管道振動外部因素

在2.3節的流固耦合仿真中,進出口邊界及支撐等四處均設置為固定,且不包括上游流體壓力脈動激勵。發動機實際運行環境復雜,管道邊界受到系統振動影響,并非為完全固定支撐,且進口包含上游流體壓力脈動激勵。因此進一步分析了進口壓力脈動激勵和基礎振動對管道振動的影響。

2.4.1 僅有外來流體壓力脈動流固耦合

(1) 邊界條件

管道邊界條件與模態計算一致,流場進口邊界條件為試驗穩定段脈動壓力數據加上平均壓力,出口為質量流量,如表3所示,流場計算以穩態結果為初場。進口壓力脈動激勵頻譜如圖15所示,包含334 Hz,910 Hz和3 639 Hz三個主要特征頻率。

表3 流場邊界條件Tab.3 Flow field boundary conditions

圖15 進口脈動壓力激勵頻譜圖Fig.15 Inlet pressure fluctuation spectrum

(2) 流固耦合振動結果

振動加速度頻譜結果如圖16所示,結果中主要包括120 Hz,228 Hz,241 Hz和335 Hz振動加速度頻率,無明顯高頻振動。振動加速度頻率和幅值結果與2.3節無明顯區別,表明進口壓力脈動激勵對振動加速度結果影響較小。

圖16 振動加速度頻譜Fig.16 Spectrum of vibration acceleration

2.4.2 熱試車工況流固耦合

(1) 邊界條件

管道邊界條件中,兩處支撐和進出口連接法蘭為振動邊界,數據采用額定工況下穩定試驗段的振動加速度,頻譜如圖17所示,包含335 Hz,909 Hz和3 635 Hz等頻率。由于出口法蘭與支撐2位置相近,因此出口法蘭采用支撐2處的振動數據。管道其余邊界條件與模態計算一致,流場邊界條件與2.4.1節一致。流場計算以穩態結果為初場,進行雙向流固耦合計算。

圖17 振動加速度頻譜結果Fig.17 The results of vibration acceleration spectrum at the support

(2) 流固耦合振動結果

振動加速度頻譜結果如圖18所示, 結果中主要包括116 Hz,211 Hz,243 Hz,327 Hz和3 637 Hz振動加速度頻率,其中211 Hz和243 Hz認為是同一頻率。振動加速度頻率結果與2.4.1節相比,增加了高頻3 626 Hz,低階頻率(小于1 000 Hz)和相應的幅值無明顯變化。結果表明高頻振動與邊界處振動相關,但邊界處振動對低階振動影響較小。

通過提取流固接觸面的節點信息,得到整個流固接觸面上振動加速度的均值分布,如圖19所示。由圖19可知,振動加速度較大的位置主要集中在進出口法蘭和兩個波紋管之間,與2.3節的結果不同。另外,中間管道的振動加速度均值較小,主要由于兩個波紋管的吸振作用。不同振動加速度頻率下的幅值分布表明,低頻下較大的幅值分布范圍更廣,而高頻下較大的幅值集中在基礎振動源附近。4個位置的基礎振動中不包含116 Hz和241 Hz兩個頻率,且兩個頻率下的幅值分布與2.3節相比發生明顯變化,因此應是與基礎振動發生耦合,使原有的幅值分布結果發生改變。同樣,335 Hz下頻率的幅值分布結果也是自激振動與基礎振動的共同作用結果。高頻下的振動加速度幅值分布可以看出,高頻基礎振動在管道結構中會迅速衰減,傳播能力較弱,因此高頻基礎振動的影響也較小。

圖19 振動加速度結果Fig.19 Vibration acceleration results

綜合所有情況的振動結果可以得出,低頻振動主要由流體自激產生,與外源激勵無關,但會與外源激勵出現一定程度的耦合。高頻振動源于基礎振動,與流體自激無關。可針對不同振動加速度頻率的產生原因制定相應的減振措施。

從圖20的應力均值分布結果可以看出,應力較大的位置集中在兩處固定支撐、彎管和波紋管處。應力較大位置是容易出現結構失效的危險位置,應重點關注。

圖20 應力均值分布結果Fig.20 Stress mean distribution results

3 結 論

本文通過雙向流固耦合數值方法,研究了液氧煤油發動機一級泵出口管路振動特性。結合管道模態、瞬態流場和振動結果,得到引起管道振動的原因及不同外源激勵對管道振動的影響。主要結論如下:

(1) 管道低頻振動由流體自激的壓力波動引起。無外源激勵時的管內流體壓力頻率接近管道固有頻率,而且振動加速度頻率與流場壓力和管道模態相近的頻率一致。另外,無外源激勵的流場壓力與振動加速度在低頻下的幅值分布結果有明顯相關性,因此管道低頻振動由流場自激的壓力波動引起。

(2) 管道外來流體壓力脈動對管道振動無明顯影響。僅有外來流場壓力脈動激勵的振動加速度結果與無外源激勵結果無明顯差異,表明外來流體壓力脈動激勵對振動結果無明顯影響。

(3) 管道高頻振動源于邊界傳遞的基礎振動。熱試車工況中的高頻振動結果,表明高頻振動與基礎振動相關。另外,高頻振動源沿管道會迅速衰減,且基礎振動對低頻振動加速度結果無明顯影響。

(4) 熱試車工況下,管道主要包括120 Hz,241 Hz,335 Hz和3 640 Hz 4個振動加速度頻率。另外應力較大的位置集中在兩處固定支撐、彎管和波紋管處,是容易出現疲勞失效的危險位置,應重點關注。

猜你喜歡
模態振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
噴水推進高速艇尾部振動響應分析
This “Singing Highway”plays music
振動攪拌 震動創新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
車輛CAE分析中自由模態和約束模態的應用與對比
國內多模態教學研究回顧與展望
高速顫振模型設計中顫振主要模態的判斷
航空學報(2015年4期)2015-05-07 06:43:35
基于HHT和Prony算法的電力系統低頻振蕩模態識別
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 日本尹人综合香蕉在线观看| 91精品啪在线观看国产60岁 | 亚洲第一区精品日韩在线播放| 99精品免费在线| 国产国产人成免费视频77777| 亚洲伊人天堂| 亚洲无码精彩视频在线观看| 91精品啪在线观看国产| 大陆精大陆国产国语精品1024| yy6080理论大片一级久久| 精品欧美视频| 欧美日韩午夜| 91成人在线观看| 亚洲天堂网在线观看视频| 日本在线欧美在线| 三级国产在线观看| 多人乱p欧美在线观看| 九色在线观看视频| 狠狠五月天中文字幕| 亚洲第一黄片大全| 韩日免费小视频| 精品国产www| 五月天香蕉视频国产亚| 亚洲天堂日本| 午夜不卡视频| 99在线视频精品| 一级全黄毛片| 亚洲无线一二三四区男男| 久久久受www免费人成| 精久久久久无码区中文字幕| 99re在线视频观看| 成年免费在线观看| 欧美激情视频一区| 国产色婷婷| 久久熟女AV| 国产精品刺激对白在线| 欧美成人亚洲综合精品欧美激情| 精品成人一区二区三区电影| 国产毛片高清一级国语 | 中文字幕欧美日韩| 国产成人艳妇AA视频在线| 亚洲区第一页| 国产门事件在线| 国产女主播一区| 99这里精品| 999在线免费视频| 国产微拍一区| 亚洲欧美成人| 精品国产自在现线看久久| 日韩欧美国产成人| 在线色国产| 欧美伊人色综合久久天天| 2019年国产精品自拍不卡| h网站在线播放| 色久综合在线| 女人一级毛片| 日韩无码真实干出血视频| 久久一日本道色综合久久| 欧美日韩中文国产va另类| 1024你懂的国产精品| 亚洲福利视频一区二区| 免费观看国产小粉嫩喷水| 国外欧美一区另类中文字幕| 欧美日韩福利| 亚州AV秘 一区二区三区| 欧美日韩高清在线| 五月婷婷精品| 国产日韩丝袜一二三区| 国产精品片在线观看手机版| 久久综合结合久久狠狠狠97色| 国产精品区视频中文字幕| 福利视频一区| 亚洲制服中文字幕一区二区| 国产精品yjizz视频网一二区| 永久免费无码成人网站| 欧美成人综合视频| 综合色在线| 国产成人av一区二区三区| 91香蕉视频下载网站| 亚洲熟妇AV日韩熟妇在线| 午夜毛片免费观看视频 | 日本不卡在线视频|