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

基于溝槽結構的泵噴推進器導管脈動壓力控制研究

2024-02-26 12:16:08葉金銘吳原潤孫大鵬鄒笑宇
船舶力學 2024年2期
關鍵詞:溝槽結構

葉金銘,吳原潤,孫大鵬,鄒笑宇

(海軍工程大學艦船與海洋學院,武漢 430033)

0 引 言

泵噴推進器是一種由轉子、定子和環狀導管構成的組合式推進裝置,在提高空泡的臨界航速、提升推進效率及降低輻射噪聲等方面具有明顯優勢。環狀導管的剖面為機翼型,根據不同航行體的工作環境、用途和性能需求,泵噴推進器的導管可分為加速型和減速型兩種[1]:加速型導管能夠將來流加速,減小轉子載荷,增加推進器推力,進而提高泵噴推進器效率;減速型導管能夠降低流體流入導管內流場的速度,提高轉子盤面處的壓強,進而有效推遲轉子葉片的空化,降低轉子的噪聲。

但泵噴推進器在實際應用中仍然存在一些問題,轉子葉梢間隙內的流場極為復雜[2-4],并時刻伴隨著渦結構的生成、發展、輸運和擴散,這些非定常的強漩渦流動不僅會引發梢渦空化,產生流動噪聲,還會使導管承受極大的脈動壓力載荷,產生結構振動噪聲。與此同時,過大的脈動壓力載荷對導管的結構強度和使用壽命產生影響,甚至會撕裂導管蒙皮。Hah[5]通過試驗研究發現在轉子尾緣后也會發生渦空化。Wu[6]使用透明材料加工制作了噴水推進器,并且在試驗中發現空化首先出現在梢渦中,隨著轉子負載增加,葉梢端面也會發生空化。軸流壓氣機在梢部流動控制方面已經得到了長久的研究,其中“機匣處理”技術是較為成熟的方式之一?!皺C匣處理”就是在壓氣機轉子梢部附近的內端壁處加工一定形式和數量的槽結構和縫結構,旨在減小葉尖處壓力面和吸力面的壓差,或者溝通葉片前緣和尾緣,減小葉片通道阻塞。Fujita[7]、Alone[8]和Ross[9]通過對不同類型的機匣處理結構進行試驗研究,發現不同類型的機匣處理結構均可以大幅度提高壓氣機的失速裕度,但會導致一定程度上的效率損失,且間隙越大,效率損失越多;Rolfes[10]針對低速壓氣機,通過數值計算和實驗研究了在兩種葉梢間隙尺寸下處理機匣結構對梢部流動的控制效果,研究結果顯示,在較大的葉梢間隙尺寸下,機匣處理結構抑制失速的效果更好,但效率相較于較小葉梢間隙尺寸時有所降低。葉金銘[11-12]等類比處理機匣技術在壓氣機中的應用,在泵噴推進器導管內壁上布置了一系列溝槽來實現對梢渦的控制,并初步驗證了溝槽結構對梢渦有控制效果;Cheng 等[13]研究了安裝在水翼頂端的懸垂凹槽對水翼性能的影響,發現在一定的間隙尺寸內能有效地抑制梢隙泄渦空化。為了更好地設計推進器導管,焦予秦等[14]針對導管螺旋槳進行了風洞試驗,測量了靜推力狀態和倒車狀態下導管內壁上的脈動壓力,并分析了其沿槳軸方向的分布規律和幅值特性;舒禮偉等[15]在循環水槽中進行了泵噴推進器導管脈動壓力測量,研究了轉子直徑和側斜角對導管脈動壓力的影響。

本文以安裝在艇后的某泵噴推進器為研究對象,在泵噴推進器轉子梢部對應的導管內壁處開設一定數量的軸向矩形溝槽,基于分離渦模擬方法并結合高質量結構化網格[16-17]分別對有無溝槽泵噴推進器進行數值模擬,通過數值計算結果分析溝槽結構對梢渦的控制效果以及對導管脈動壓力的控制效果。

1 數值計算方法及其驗證

1.1 湍流模型與數值離散

本文使用分離渦模擬(DES)方法[18-20],DES作為一種混合數值計算模型,它將LES方法和RANS方法的優點結合到了一起,將邊界層內的“附著”渦進行建模,利用RANS方法處理,并采用LES方法對邊界層外較大的“分離”渦進行數值模擬,因此DES 方法可以使用相對較少的計算資源取得相對精確的流場信息。在DES模擬中選用改進延時分離渦模擬(IDDES)進行數值計算,采用SIMPLE 算法完成壓力-速度耦合方程的求解,運動方程中的時間項通過二階隱式格式離散,對流項通過二階迎風格式進行離散,粘性力項通過中心差分格式進行離散,并采用SSTk-ω湍流模型對N-S 方程進行封閉。該模型對流場中渦結構的捕捉效果較好。

1.2 數值計算方法驗證

為了對數值計算方法的準確性進行驗證,本文選取NACA 16020 翼型的橢圓形水翼為研究對象,使用改進延時分離渦模擬(IDDES)方法計算梢渦的流場分布,并與試驗值[21]進行對比,來流速度為10 m/s,橢圓水翼的攻角為10°,最大弦長C=0.475 m,展長L=0.7125 m。原點在翼根弦長中心處,入口和出口與翼導邊和隨邊距離分別為2.5C、10.0C,計算域及坐標系如圖1 所示。進口設置為速度進口,出口設置為壓力出口,橢圓水翼表面設置為無滑移壁面。

應用O型網格對翼壁面附近的網格進行處理,并對梢渦軌跡所在位置進行加密,整個橢圓水翼計算域的總網格數為1000萬。計算穩定后,提取梢渦流場處的軸向速度和切向速度,圖2和圖3分別為試驗結果與計算結果在x方向3 個位置處的梢渦流場無因次軸向速度(Vx/U∞)和周向速度(Vt/U∞)沿z方向的分布情況。由圖2 和圖3 可以看出,IDDES 方法對梢渦流場的預報效果較好,尤其是在x/C=0.1 處梢渦空泡的初生位置附近,計算得到的渦核軸向速度和切向速度分布與試驗值也能很好地吻合。因此,上述對橢圓水翼梢渦流場計算結果與試驗結果的比較分析,驗證了數值計算方法的可靠性和適用性。

圖2 不同位置處的梢渦軸向速度分布Fig.2 Axial velocity in tip vortex at different axial positions

圖3 不同位置處的梢渦切向速度分布Fig.3 Tangential velocity in tip vortex at different axial positions

2 艇后泵噴推進器梢渦流場數值模擬

2.1 計算模型

以某水下航行體及其泵噴推進器為計算模型,如圖4 所示,該泵噴推進器采用的是前置定子和加速型導管,轉子與定子的葉剖面都是NACA 翼型,轉子為7 葉,定子為13 葉,轉子直徑為534.4 mm,導管與葉梢之間的最小間隙為4 mm,最大間隙為5.6 mm。

圖4 水下航行器及泵噴推進器模型Fig.4 Model of underwater vehicle and pumpjet propulsor

在導管內壁開設一定數量的矩型溝槽,由于以導管內壁母線為曲線,溝槽的徑向深度沿軸向會發生變化,將轉子葉梢前緣所對應的溝槽深度定義為溝槽深度H,轉子葉梢軸向長度L0為36.96 mm。溝槽結構的示意圖如圖5~6所示,溝槽深度H為12.8 mm,溝槽寬度B為9.6 mm,溝槽數量N為100。

圖5 溝槽結構示意圖Fig.5 Schematic diagram of groove structures

圖6 溝槽結構布置示意圖Fig.6 Schematic diagram of arrangement of groove structures

2.2 計算域劃分與邊界條件設置

計算區域包含了外域和導管內域,導管內域是指導管前緣進口截面至尾緣出口之間的流場域,如圖7(a)所示。導管內域以外的流場域稱為外域,外域的大小以及邊界條件設置如圖7(b)所示,以艇長L為基礎尺寸,外域輪廓是半徑為1L的圓柱體,將艇首前側的外域邊界設置為速度進口,距艇首為1L;艇尾后側的外域邊界設置為壓力出口,距艇首為3L;外域的圓柱面邊界設置為對稱面。

圖7 計算域劃分形式Fig.7 Division of computational domain

導管內域又分為定子域、轉子域、葉頂域、后域等子域,如圖7(a)所示。當沒有溝槽結構時,葉頂域只包含導管內壁與轉子域之間的回轉體流場域,當溝槽結構存在時,葉頂域既包含了導管內壁與轉子域之間的回轉體流場域,又包含了軸向縫內部的流場域。在導管內域中,導管、定子和轉子壁面設置為無滑移壁面,對于IDDES 的壁面處理選擇Ally+壁面處理,最小允許壁面距離為1.0×10-6m;計算時間步長設置為2.0×10-4s,單位時間步長內轉子旋轉0.9°。

2.3 計算域網格劃分

在上述計算域劃分的基礎上,采用結構化網格對各區域進行網格劃分,外域網格和后域網格如圖8(a)和圖8(b)所示。為了保證網格的周期性或者對稱性,定子域、轉子域和葉頂域均可以被看作由單通道模型周期性旋轉得到,所以先進行單通道模型網格的劃分,再根據其周期性進行旋轉,得到其全通道模型網格,如圖8(c)~8(f)所示。

圖8 計算域網格網格Fig.8 Mesh of computation domain

外域、后域和定子域的網格數量分別為718萬、450萬和460萬。由于轉子梢部附近區域和葉頂間隙內的流場比較復雜,為了較為準確地對轉子梢渦流動進行計算,對轉子梢部附近區域以及葉頂間隙內區域進行了相應加密,設置四組不同網格密度的網格,對轉子域和葉頂域的網格進行無關性驗證,最終確定轉子域和有無溝槽葉頂域網格數量分別為1500萬、920萬和560萬。其中有無溝槽的葉頂域重疊部分網格相同,單個流道的軸向網格節點為90,周向網格節點為52,徑向網格節點為15;有溝槽葉頂域的溝槽的軸向網格節點為60,周向網格節點為28,徑向網格節點為25。得到的兩套計算網格總數分別為無溝槽的3698萬和有溝槽的4058萬。

2.4 梢渦控制效果研究

為了研究溝槽結構對梢隙流場中渦結構的影響,監測并比較有無溝槽狀態下泵噴推進器轉子葉梢端面、梢泄渦渦核和梢渦渦核的最小壓力變化。

首先對比有無溝槽狀態下葉梢端面壓力分布情況,圖9為同一時刻、同一轉子的葉梢端面壓力分布,從圖中可以看出,轉子葉梢端面前緣位置處的壓力極低,這主要是因為在轉子高速旋轉時,流體會在轉子葉梢前緣位置處發生劇烈的流動分離并形成梢隙分離渦,從而導致該位置壓力極低。在導管內壁開設溝槽結構后,可以發現溝槽結構不但提高了葉梢前緣位置處的最低壓力,還使整個葉梢端面上的平均壓力有所提高。同時繪制一個周期內轉子葉梢端面最低壓力的時域曲線,如圖10 所示,可以發現有溝槽狀態下轉子葉梢端面上的最低壓力在任意時刻均較無溝槽結構時有所提高。

圖9 轉子葉梢壓力分布Fig.9 Pressure distribution on rotor tips

圖10 轉子葉梢端面最低壓力Fig.10 Minimum pressure at tip face of rotor

以轉子半徑R為參考基礎尺寸,設轉子葉梢尾緣的軸向位置為0,在其前后不同軸向距離Δx的位置處建立一系列垂直于轉子軸線的橫截面,各橫截面的直徑與該位置處導管內壁的直徑一致,其中部分橫截面位置如圖11 所示。待計算穩定后,取同一時刻兩種泵噴推進器在同一角度位置葉梢附近的各橫截面壓力分布,如圖12所示。從圖12中可以看出,在導管內壁開設軸向溝槽結構后,轉子尾緣后的渦核壓力有明顯提升,轉子尾緣前的渦核壓力也有一定提升。

圖11 橫截面示意圖Fig.11 Diagram of cross sections

圖12 各橫截面上渦核壓力分布比較Fig.12 Comparison of vortex core pressure distributions on each cross section

由于轉子尾緣后的梢泄渦螺距角較小,為了便于觀察和分析軸向溝槽結構對梢泄渦渦核壓力的影響,建立一系列通過轉子軸線的軸向剖面,剖面之間的夾角為2°,如圖13 所示??梢钥闯?,在轉子梢部尾緣附近及順著梢泄渦流動方向的軸向剖面上,有軸向溝槽結構時的渦核壓力比無軸向溝槽時的渦核壓力有明顯升高。

圖13 各軸向剖面上渦核壓力分布比較Fig.13 Comparison of vortex core pressure distributions on each axial section

為了詳細分析軸向溝槽結構對梢泄渦渦核壓力的提升效果,待兩種結構的泵噴推進器轉子梢部流動計算結果穩定后,取二者一個轉子旋轉周期內各橫截面上的最低渦核壓力的時均值和時域最小值進行比較,如圖14所示。從圖14可以看出,在導管內壁開設軸向溝槽結構后,各橫截面上的最低渦核壓力均有提升,其中在轉子梢部尾緣附近及以后的橫截面(Δx≥-0.01R)上,渦核壓力提升最明顯,在葉梢導邊附近的橫截面(Δx≤-0.12R)上,渦核壓力提升效果次之,在葉梢端面中部附近位置的橫截面(-0.11R≤Δx= -0.02R)上,渦核壓力提升效果較小。

圖14 在各橫截面位置處最低渦核壓力的時均值和時域最小值比較Fig.14 Comparison of the time-averaged and time-domain minimum values of the minimum vortex core pressure at each cross-section position

同時采用基于Pressure 標量的等值面法將梢隙流場中的渦核低壓體積進行可視化,取Pressure為-120 000 Pa、-150 000 Pa 和-180 000 Pa,即可得到導管內流場中的渦結構形態以及低壓區域體積的變化情況,如圖15所示。從圖15中可以看出,導管內壁開設軸向溝槽結構后,遠離葉梢的梢泄渦所在位置的低壓區域體積明顯變小,這說明軸向溝槽結構確實可以抑制梢泄渦發展,顯著減小梢泄渦的低壓區范圍,有利于控制梢泄渦空化。

圖15 渦結構位置及低壓區域體積對比Fig.15 Comparison of vortex structure locations and low pressure area volumes

2.5 脈動壓力控制效果研究

通過前面的分析可以看出,軸向溝槽結構能夠提高泵噴推進器轉子梢部端面最低壓力和轉子梢泄渦渦核的壓力,從而抑制轉子梢泄渦空化。

轉子梢泄渦與導管內壁距離較近,梢泄渦渦核處的低壓會引起梢泄渦附近的導管內壁區域壓力較低。由于軸向溝槽結構能夠提高轉子梢泄渦渦核壓力,梢泄渦渦核附近的導管內壁最低壓力也會明顯增大,兩種泵噴推進器在同一時刻的導管內壁壓力分布如圖16 所示。從圖16 可以看出,無軸向溝槽結構時,轉子梢泄渦附近的導管內壁壓力較低,在同一軸向位置處,導管內壁壓力在周向上的不均勻程度較大;有軸溝槽結構時,轉子梢泄渦附近的導管內壁壓力明顯升高,在同一軸向位置處,導管內壁壓力在周向上的不均勻程度明顯降低。

圖16 泵噴推進器導管內壁壓力分布Fig.16 Pressure distributions on the inner wall of the duct of pumpjet propulsor

在轉子旋轉的過程中,梢泄渦系的空間位置會隨轉子的轉動而發生變化,導管內壁面上低壓區域的空間位置也會隨之發生周期性變化,這導致了導管內壁面上的壓力出現大幅度的周期性脈動。轉子梢泄渦核壓力越低,導管內壁最低壓力也會隨之降低,轉子旋轉過程中,導管內壁壓力脈動幅度越大。為了研究溝槽結構對導管內壁面脈動壓力的控制效果,在轉子葉梢上方的導管內壁面上布置了兩組監測點,當導管內壁有溝槽結構時,A 組監測點在溝槽周向正中間位置,B 組監測點在相鄰兩溝槽周向正中間位置,其周向位置如圖17所示。

圖17 兩組監測點在導管內壁周向上分布Fig.17 Two groups of monitoring points distrib?uted circumferentially on the inner wall of the duct

以轉子葉梢導邊正上方導管內壁處的軸向位置為基點(軸向坐標設置為0),隨邊正上方導管內壁處的軸向位置為L0,在導管內壁上建立A 組的9 個監測點,各點在導管上的軸向分布如圖18 所示,無溝槽與有溝槽對應監測點的軸向位置相同。由于A組監測點在溝槽正中間位置,所以部分監測點位于溝槽內壁中。

圖18 導管內壁脈動壓力A組監測點Fig.18 Group A monitoring points of pulsating pressure on inner wall of duct

因B 組監測點在周向上位于相鄰兩溝槽正中間位置,所以各監測點均分布在導管內壁上,B 組中的監測點數量與A 組相同,各點的軸向位置與A 組中對應監測點的軸向位置也相同,如圖19所示。

圖19 導管內壁脈動壓力B組監測點Fig.19 Group B monitoring points of pulsating pressure on inner wall of duct

待計算結果穩定后,提取一個轉子旋轉周期內各監測點上的脈動壓力時域曲線,如圖20所示,可以看出各監測點的壓力脈動呈現出明顯的葉頻特性。在導管內壁開設溝槽結構后,A 組中各監測點的壓力脈動幅度均有明顯減小,在轉子葉梢導邊之前的區域,脈動幅度降低量較小,如圖20(a)中A1時域曲線所示;在轉子葉梢中前位置至轉子葉梢尾緣附近的脈動幅度降低量最大,如圖20(a)中A3~A6時域曲線所示。B 組檢測點中,在轉子葉梢導邊之前區域,有溝槽結構的導管內壁脈動壓力脈動幅度降低量較小,如圖20(b)中B1 時域曲線所示;在轉子葉梢導邊與隨邊之間,有無溝槽結構的導管內壁脈動壓力脈動幅度沒有明顯改變,如圖20(b)中B2~B5 時域曲線所示;導管內壁開設溝槽結構后,在轉子葉梢隨邊之后區域,導管內壁上的壓力脈動幅度有明顯減小,如圖20(b)中的B6~B9時域曲線所示。

圖20 導管內壁脈動壓力時域曲線Fig.20 Time domain curve of fluctuating pressure on inner wall of duct

對導管內壁脈動壓力時域曲線進行Fourier分析,即可得到各監測點的脈動壓力的頻域特征,如圖21 所示。從圖21(a)中可以看出,在導管內壁開設溝槽結構后,A 組監測點上脈動壓力的各階葉頻幅值均有明顯減小,在梢泄渦被抑制區域(A3~A9 區域),各階葉頻幅值減小量最大。從圖21(b)可以看出,在導管內壁開設溝槽結構后,監測點中B3、B4 的一階葉頻幅值有所增大,但二階和三階葉頻的幅值均有減小,且這兩點的脈動壓力時域曲線的波動幅度基本一致;其余各點的各葉頻幅值均有所減小,其中轉子梢部隨邊之后的B6、B7、B8、B9四點最為明顯。

圖21 導管內壁脈動壓力頻域幅值對比Fig.21 Comparison of frequency domain amplitude of pulsating pressure on inner wall of duct

由此可見,導管內壁的軸向溝槽結構對導管內壁壓力脈動基本沒有負面影響,且在導管內壁中的大多數區域能有效降低脈動壓力的幅值。導管脈動壓力既可以激發流動噪聲,也是引起導管結構振動噪聲的激勵源。因此,通過導管內壁的軸向溝槽結構降低導管內壁的壓力脈動輻值,有望成為降低推進器導管流動噪聲和結構振動噪聲的有效途徑。

2.6 泵噴推進器水動力性能比較

航速為Vs時不帶推進器的裸艇體阻力記為R0,由于泵噴推進器與艇體間存在相互作用,會使艇體阻力變大,此時的艇體阻力記為Rs,當推進器推力T與艇體阻力Rs平衡時可以使用下式來計算泵噴推進器的效率ηa:

式中,n為泵噴推進器的轉速,Q為轉子的扭矩。

由于在CFD計算時,艇體阻力Rs和推進器推力T可能略有差別,本文中裸艇體的阻力R0=6215 N,無軸向縫狀態下,艇體阻力Rs=10 051.92 N,推進器推力T=10 043.57 N;有軸向縫狀態下,艇體阻力Rs=9974.74 N,推進器推力T=9907.2 N。無軸向縫狀態和有軸向縫狀態時,艇體阻力Rs和推進器推力T分別相差0.08%和0.68%,雖然艇體阻力Rs和推進器推力T存在差別,但是相差很小,可以通過將公式(1)進行修正來計算泵噴推進器效率ηa,修正后的泵噴推進器效率計算公式如下:

泵噴推進器的轉子推力系數KTb、扭矩系數KQ、推進器推力系數KTa的計算公式分別如下:

式中,Tb為泵噴推進器轉子的推力,Q為轉子的扭矩,ρ為水的密度,n為泵噴推進器的轉速,Vs為艇的航速,D為泵噴推進器轉子直徑。

泵噴推進器水動力性能計算結果如表1所示。在導管內壁開設軸向縫之后,泵噴推進器的KTb、KQ、KTa均有所減小,ηa略有增加,其中KTb減小了2.33%,KQ減小了1.05%,KTa減小了1.38%,ηa增加了0.46%,這說明軸向縫結構對泵噴推進器的水動力性能影響較小,在推進器效率上還有一定的改善作用。

表1 泵噴推進器水動力性能對比Tab.1 Comparison of hydrodynamic performance of pumpjet propulsors with and without grooves

3 結 論

本文借鑒航空發動機中的壓氣機“處理機匣”技術,在泵噴推進器導管內壁開設了一定數量的軸向溝槽結構,用于減弱泵噴推進器的梢渦,從而減小泵噴推進器導管內壁上的脈動壓力。結合分離渦模擬(DES)方法和滑移網格方法,分別對有無溝槽結構的泵噴推進器進行了數值模擬,通過對計算結果進行對比分析,研究了溝槽結構對梢渦渦核壓力、導管內壁脈動壓力和水動力性能的影響,得到了如下結論:

(1)在導管內壁開設溝槽結構后,溝槽結構不但提高了葉梢前緣位置處的最低壓力,還使整個葉梢端面上的平均壓力有所提高,能有效抑制流體在轉子葉梢前緣位置處發生流動分離并形成梢隙分離渦。

(2)溝槽結構可以顯著提高梢渦渦核壓力,減小梢隙泄渦和梢渦位置處低壓區域的體積,即溝槽結構可以抑制梢渦空化,當發生梢渦空化時,溝槽結構也可以減小梢渦空化體積。

(3)導管內壁的溝槽結構能有效降低導管內壁的脈動壓力幅值,降低引起導管結構振動的激勵源。

(4)導管內壁開設軸向縫使推進器的效率略有增加,說明軸向縫能夠在對泵噴推進器推進性能影響較小的情況下較為明顯地控制轉子的梢部流動,為泵噴推進器的空化抑制提供了新途徑。

猜你喜歡
溝槽結構
基于數值模擬的2種條紋溝槽減阻特性對比分析
一種具有多形式鋼片結構的四季胎
輪胎工業(2021年10期)2021-12-24 17:23:35
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
一種低噪音的全路況輪胎
輪胎工業(2020年9期)2020-03-01 18:58:44
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
二層溝槽織構對機床導軌表面潤滑特性的影響
論《日出》的結構
溝槽爆破參數優化及成本分析
Influence of machining parameters on groove surface morphology of condenser for heat column
機床與液壓(2015年3期)2015-11-03 07:02:03
主站蜘蛛池模板: 国产亚卅精品无码| 日本午夜视频在线观看| 亚洲国产综合精品中文第一| 美女内射视频WWW网站午夜| 国产va在线观看免费| 无码视频国产精品一区二区| 亚洲成a人片77777在线播放| 国产电话自拍伊人| 亚洲国产AV无码综合原创| 又黄又湿又爽的视频| 国产一级一级毛片永久| 免费看美女自慰的网站| 国产系列在线| 国产精品免费电影| 国产十八禁在线观看免费| 亚洲无码在线午夜电影| 国产欧美日韩在线一区| 噜噜噜久久| 久久精品国产一区二区小说| 老司机精品一区在线视频| 露脸国产精品自产在线播| 国产毛片一区| 国产精品成| 国产粉嫩粉嫩的18在线播放91| 99精品免费欧美成人小视频 | 日本国产精品一区久久久| 91视频99| 日韩欧美国产成人| 亚洲成人免费看| 亚洲精品中文字幕午夜| 亚洲午夜久久久精品电影院| 99久视频| 国产精品久久精品| 国产成人精品优优av| 亚洲狠狠婷婷综合久久久久| 97在线碰| 国产精品女在线观看| 国产精品久久久久久久久kt| 国产情精品嫩草影院88av| 亚洲青涩在线| 国产精品嫩草影院视频| 欧美69视频在线| 美女高潮全身流白浆福利区| 午夜国产小视频| 国产精品第一区| 巨熟乳波霸若妻中文观看免费| 在线国产资源| 国产成人乱无码视频| 无码精品一区二区久久久| 久久semm亚洲国产| 福利国产在线| 国产成人精品在线1区| 精品福利视频导航| 亚洲国产精品久久久久秋霞影院| 国产麻豆永久视频| 国产在线精品香蕉麻豆| 超清无码熟妇人妻AV在线绿巨人| 一区二区偷拍美女撒尿视频| 2020国产精品视频| 国产一级二级在线观看| 国产一在线观看| 午夜视频免费试看| 国产精品网拍在线| 精品久久国产综合精麻豆| 波多野结衣无码中文字幕在线观看一区二区 | 99久久精品国产综合婷婷| 伊人久久综在合线亚洲2019| 国产在线视频欧美亚综合| 国产视频久久久久| 日韩黄色大片免费看| 精品福利视频网| 1769国产精品视频免费观看| 久久精品丝袜| 玖玖精品在线| 亚洲人成成无码网WWW| 一本大道视频精品人妻| 香蕉eeww99国产在线观看| 亚洲日本www| 午夜性刺激在线观看免费| 99re热精品视频国产免费| 亚洲国内精品自在自线官| 四虎国产永久在线观看|