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

泵噴分布式脈動壓力激勵下泵噴艇體耦合系統振動聲輻射

2021-10-11 09:49:22黃修長蘇智偉師帥康饒志強華宏星
振動與沖擊 2021年18期

黃修長, 蘇智偉, 師帥康, 饒志強, 華宏星

(1.上海交通大學 振動、沖擊、噪聲實驗室,上海 200240;2.上海交通大學 機械系統與振動國家重點實驗室,上海 200240)

螺旋槳、泵噴等推進器在表面脈動壓力激勵下,一方面會產生結構振動直發聲以及流噪聲直發聲;另一方面會通過軸系、定子等傳遞激勵艇體,導致艇體產生強烈的聲輻射。因此,推進器-軸系-艇體的耦合振動聲輻射研究受到廣泛關注[1-3]。

目前已針對螺旋槳-軸系-艇體耦合系統的振動噪聲特性開展豐富的研究,明確了螺旋槳激勵下螺旋槳-軸系-艇體耦合系統的聲振響應特性[4-5]。但上述分析均采用在螺旋槳槳轂或0.7R處施加集中力(R為螺旋槳半徑),無法考慮分布式載荷的相位不同以及螺旋槳空間模態特性對響應的影響。為考慮上述問題,部分學者在分析螺旋槳-軸系耦合系統的響應或螺旋槳的流固耦合響應時,采用條帶法將隨機寬帶力映射到螺旋槳結構網格上進行施加[6],或利用面元法和有限元法(finite element method,FEM)耦合將計算得到的脈動壓力施加到結構網格[7-8]。

泵噴推進器作為一種組合推進器,由于具有較螺旋槳更優越的推進及靜音性能而被廣泛應用。泵噴推進器的激勵力更為復雜,其脈動壓力是在導管、定子、轉子吸力面和壓力面上呈復雜空間分布的力,在該空間分布力作用下會激勵起不同模態產生復雜的響應,其響應特性和響應機理更加復雜。目前針對泵噴推進器的研究主要集中在3個方面:一是水動力性能、表面脈動壓力激勵特性、流場特性等的研究[9-10];二是泵噴參數(如間隙、導管幾何參數和定轉子幾何參數等)對水動力學性能、脈動壓力等的影響[11-12];三是推進器脈動壓力激勵下艇體的聲輻射特性以及新型結構(如鋸齒后緣導管)的降噪效果[13]。但對于泵噴激勵下泵噴-軸系-艇體耦合系統振動聲輻射的研究鮮見公開報道。對于泵噴推進器激勵下的艇體振動聲輻射響應,泵噴推進器同時通過推進軸系激勵艇體以及通過導管-定子激勵艇體,這兩種途徑的傳遞特性、傳遞幅值相對大小等仍不明確。

本論文考慮suboff艇后泵噴推進器的激勵力為轉子、導管和定子上的分布式脈動壓力分布,通過插值的方法把流體脈動壓力映射到推進器結構濕表面,建立了泵噴-軸系-艇體的耦合動力學模型,獲得了其結構振動聲輻射響應,并分析了通過軸系激勵艇體、通過導管-定子激勵艇體這兩條傳遞途徑的振動聲輻射響應特性。

1 理論基礎

1.1 suboff艇后泵噴表面脈動壓力計算

不考慮推進器內部流場空化影響,可假設流體為不可壓縮黏性單相流,基于雷諾平均納維-斯托克斯(reynolds-averaged navier-stokes,RANS)的控制方程可寫為[14]

式中:i,j=1,2,3;ρ為流體密度;xi和xj分別為笛卡爾坐標系下3個方向的坐標分量;P為雷諾平均壓力;Uj和u′i分別為平均速度和脈動速度;μ為流體動力黏度系數;為雷諾應力項,為了使式(2)封閉,引入渦流黏度vt,雷諾應力可寫為

式中:k為雷諾時間平均下湍流動能;Cμ=0.09。采用剪切應力輸運(shear stress transfer,SST)k-ω湍流模型進行模擬。

1.2 泵噴-軸系-艇體耦合系統振動聲輻射模型

采用耦合有限元法開展泵噴-軸系-艇體耦合系統的動力學建模

式中:M和K分別為泵噴-軸系-艇體耦合系統的質量和剛度矩陣;C為阻尼矩陣,C=αM+βK;X為結構的位移向量;p為流場中的聲壓向量;F為泵噴的分布式脈動激勵力;Ma和Ka分別為流體的質量矩陣和剛度矩陣;H為流固耦合矩陣;∏建立了結構濕表面上的節點法向位移與結構振動位移向量X之間的關系。Ma,Ka和H的表達式分別為

式中:N和Nψ分別為結構和流體的插值形函數;ρa為流體密度;c為流體中波速。

由結構濕表面的法向振速viα,利用直接邊界元法(boundary element method,BEM)求得濕表面上的聲壓piα以及場點聲壓pj[15]

式中:J(ξ)為Jacobian矩陣;Ai為第i個單元的面積;(ξ)為第j個場點與泵噴-軸系-艇體耦合系統結構濕表面上第i個單元的第α個節點之間的距離;N(ξ)為邊界元中聲壓的插值形函數。

1.3 表面分布式脈動壓力和結構濕表面節點力映射

將計算流體動力學(computational fluid dynamics,CFD)模型中轉子、定子、導管表面單元的壓力以及泵噴結構濕表面的節點導出,通過徑向基函數插值將流體單元節點的壓力映射到泵噴結構濕表面節點上。假設空間存在一組點群xj及點群上的函數值fj(j=1,2,…,N),可構建各點之間的距離函數?(‖x-xj‖)與函數值之間的關系

式中:?(‖x-xj‖)為徑向基函數,它是插值點與目標點之間距離的函數;λj為徑向基函數的系數。通過xj可建立點群之間的λj系數求解線性方程組。空間任意點上的函數值f(x)則可通過式(8)求解。具體計算時首先找到泵噴結構濕表面單元中心附近的一組流體節點及壓力值,并求解徑向基函數插值系數,然后將泵噴結構濕表面單元中心代入式(8)中求解泵噴結構濕表面單元上的壓力分布。在求得泵噴結構濕表面單元的壓力分布后,將壓力乘以面積得到面單元上的集中力,將集中力按節點數平均分配到各個節點上,即可得到每個節點上的集中力分布。在每個時刻進行徑向基函數插值即可得到時域中的結構濕表面節點的三向集中力。利用傅里葉變換可得到頻率域結構濕表面節點的三向集中力。

2 數值模型

研究對象為suboff艇后泵噴推進器模型,如圖1所示。泵噴-軸系-suboff艇參數,如表1所示。殼體厚度為6 mm,并采用T型加強筋進行加強,縱向方向有4個艙壁;泵噴推進器通過定子支撐在殼體艉部。泵噴-軸系模型見圖1,推進軸系的長度為0.89 m、直徑為0.022 m,通過艉后軸承支撐在殼體艉部以及通過推力軸承支撐在最后一個隔艙壁的中心。艉后軸承的橫向和垂向剛度分別為4.5×109N/m,1.7×108N/m;推力軸承的縱向、橫向和垂向剛度分別為6.6×108N/m,2.3×108N/m和6.7×108N/m。殼體、加強筋、軸系、泵噴定子和導管的材料為鋼,其密度、彈性模量和泊松比分別為7 850 kg/m3,2.1×1011Pa,0.3。 轉子的材料為青銅,其密度、彈性模量和泊松比分別為8 790 kg/m3,1.1×1011Pa,0.3。 結構阻尼均取為0.001。

圖1 泵噴-軸系-艇體耦合系統Fig.1 Pump-jet-shafting-suboff system

表1 泵噴-軸系-suboff艇參數Tab.1 Parameters of pump-jet-shafting-suboff system

2.1 CFD模型和計算驗證

創建泵噴-軸系-suboff的CFD計算區域,如圖2(a)所示,將計算域劃分成靜止域和旋轉域,靜止域包含定子、艇體計算域和導管計算域,將轉子所在旋轉域設置為滑移網格,兩個計算域之間的交界面設置為interface邊界條件,進行數據傳遞。計算中遠前方入口設置為速度入口條件,遠后方出口設置為壓力出口條件,外圍遠場邊界與入口一樣為速度入口條件,轉子、定子及導管等固體表面采用無滑移固體壁面條件。部分進行流體計算的面網格,如圖2(b)所示。結合網格無關性分析結果,選擇合適的網格密度進行流域網格劃分。計算域網格共計2 760萬,其中:艇體域網格為1 000萬;泵噴定子域網格數量為660萬;泵噴轉子域網格為1 100萬。艇體表面和泵噴推進器均劃分精細邊界層網格,目標Y+值為1。

圖2 suboff艇后泵噴的計算域和CFD網格Fig.2 Computational domain and CFD model of the behind-hull suboff

計算工況為進速9.42 m/s,轉子轉速1 824 r/min,進行非定常計算,設置時間步長為4.568 7×10-5s,保證轉子每個時間步旋轉0.5°,計算的最高頻率為10 944 Hz,采樣頻率為21 888 Hz。數值計算的收斂性殘差滿足:連續性殘差收斂在10-5以下,3個方向速度殘差均在10-6以下,湍動能k殘差低于10-7。

CFD數值求解采用商業軟件ANSYS Fluent進行,采用SST k-ω湍流模型,壓力、速度耦合求解采用SIMPLEC(semi-implicit method for pressure-linked equations consistent)方法,其他項均采用二階迎風格式。

采用對suboff艇和泵噴單獨進行CFD計算驗證的方法對建立的模型進行驗證。對suboff艇體進行驗證時,保留圖2中的suboff網格,流域和網格密度與圖2相同。采用上述網格模型對進速3.051 m/s的工況進行計算,得到全附體suboff艇阻力(103.2 N)與實測阻力(102.3 N)誤差為0.9%[16],同時得到的槳盤面處的軸向伴流結果與試驗結果相似,如圖3所示。對泵噴進行驗證時,采用與圖2相同的流域、網格密度(其中:Grid2與圖2中網格數量相同;Grid3將轉子域中網格增加40%)和計算設置,對一五葉前置定子七葉轉子泵噴進行計算,如圖4所示。將計算結果與空泡水筒中的試驗結果進行對比,采用Grid2推力最大誤差為3.59%,最小誤差為0.23%,扭矩最大誤差為5.5%,效率的最大誤差為5.17%,精度較高,增加網格數量后(Grid3)增加了計算成本,精度提高不明顯,故采用Grid2滿足分析要求。

圖3 槳盤面處軸向伴流對比Fig.3 Comparison of axial velocity in the rotor plane

圖4 泵噴水動力性能計算模型和試驗驗證Fig.4 CFD model and experimental validation for propulsion performance of the pump-jet

2.2 結構動力學模型

在ABAQUS軟件中建立泵噴-軸系-艇體耦合系統的結構有限元模型,如圖5所示。殼體采用33 861個S4R單元進行模擬,加強筋采用4 354個B31單元進行模擬,艙壁采用4 362個S4R單元進行模擬。定子、轉子和導管采用22 641個C3D8實體單元進行模擬。推進軸系采用91個B31梁單元進行模擬。為考慮流固耦合,在泵噴-軸系-艇體耦合系統的結構有限元模型外建立圓柱形的水體,在結構有限元模型的濕表面和水體模型的表面之間建立壓力和速度耦合條件。其中圓柱形水體的長度和直徑分別為6 600 mm,2 500 mm(直徑為6倍導管最大直徑)。水體采用2 601 967個AC3D8聲學單元進行模擬。水的密度和體積模量分別為1 000 kg/m3和2×109Pa。在水體外表面施加無反射聲學邊界條件。將泵噴導管、定子和轉子濕表面以及艇體濕表面的速度導出到LMS VirtualLab中,利用直接邊界元法進行聲輻射的計算。結構單元和聲輻射單元的網格大小為5 mm,根據拇指法則,結果在600 Hz以內具有可信度[17]。

圖5 泵噴-軸系-艇體耦合系統結構振動聲輻射模型Fig.5 Vibro-acoustic model of coupled pump-jet-shafting-suboff system

3 仿真結果與分析

3.1 表面脈動壓力

取出某個轉子葉片表面0.7R,0.8R,0.9R和相應轉子區域對應的導管內壁面某點的脈動壓力,其頻譜結果如圖6所示。由圖6可知,葉片的脈動壓力以軸頻30.4 Hz及其倍頻為主,導管上的脈動壓力以葉頻273.7 Hz(9倍軸頻)和2倍葉頻547.4 Hz為主。總體而言,導管上的脈動壓力幅值比轉子葉片上監測點的脈動壓力幅值大1~2個數量級。

圖6 轉子不同半徑處和導管的表面脈動壓力頻譜Fig.6 Spectrum of pressure on the different radius of rotor surface and the duct

轉子縱向、橫向和垂向合力的時程曲線和頻譜,如圖7所示。由圖7可知,由于不同轉子葉片上的脈動壓力相位不同,疊加后各個方向的脈動力在葉頻及其倍頻處出現較大的峰值,在軸頻及其倍頻處也表現為峰值,縱向脈動力的幅值大于垂向和橫向脈動力。隨著頻率的增加,脈動力的幅值呈衰減趨勢。

圖7 轉子三向脈動力的時程曲線和頻譜Fig.7 Time-history,spectrum of pulsation force in three directions for the rotor

圖8給出了CFD計算得到的表面脈動壓力結果和插值后導入到ABAQUS中結構濕表面脈動壓力結果的對比,給出了轉子壓力面和導管內壁面的結果。由圖8可知,采用徑向基函數插值后,脈動壓力的分布和幅值均能夠保證足夠的精度。

圖8 分布式脈動壓力的CFD結果和FEM結果比較Fig.8 Comparison of fluctuation pressure distribution on surface of propeller,duct between CFD and mapped result in FEM

3.2 模態結果與分析

對泵噴-軸系-艇體耦合系統動力學模型開展模態分析,得到典型模態結果如圖9、表2所示。

表2 耦合系統模態頻率列表 Tab.2 List of modal frequencies of the coupled system Hz

圖9 泵噴-軸系-艇體耦合系統典型模態Fig.9 Typical modal shapes of the coupled pump-jet-shafting-suboff system

下面計算泵噴-軸系-艇體耦合系統的振動聲輻射響應。考慮4種載荷施加方法:一是將導管、轉子和定子上的脈動壓力施加到結構網格上(‘Pumpjet’);二是僅將轉子上的脈動壓力施加到結構網格上(‘Ro-tor’);三是對轉子葉片上的脈動壓力進行求和,將每個葉片上的三向集中力施加到轉子葉片0.7R處(‘0.7R’);四是對整個轉子的脈動壓力求和,將整個轉子的三向集中力施加到槳轂中心(‘Hub’)。

3.3 通過軸承的傳遞力結果與分析

對通過各個軸承的傳遞力進行分析。推力軸承和艉后軸承在不同載荷施加方式下的多向傳遞力,如圖10所示。對于推力軸承,通過其傳遞的三向力對于‘Pumpjet’和‘Rotor’基本重合,‘0.7R’將使結果偏大;‘Hub’將使212.9 Hz以前的結果偏小。推力軸承縱向傳遞力的最大響應峰值為軸頻30.4 Hz。5個顯著的峰值分別為209.9 Hz(軸系一階縱振),212.9 Hz(7倍軸頻),243.3 Hz(8倍軸頻)和273.7 Hz(9倍軸頻,即葉頻)以及547.5 Hz(18倍軸頻,即2倍葉頻)。軸系一階縱振的峰值頻率顯著和相鄰的峰值具有可比性。一些小的響應峰值出現在334.6 Hz(艇體一階縱振),360.4 Hz,391.3 Hz(軸系二階縱振),486.6 Hz和577.9 Hz。對于推力軸承垂向傳遞力,‘0.7R’和‘Hub’施加方式在葉頻273.7 Hz處出現峰值,但是在‘Pumpjet’和‘Rotor’情況下并未出現。主要的峰值頻率為30.4 Hz,118.6 Hz(軸系一階彎曲),121.7 Hz(4倍軸頻),243.3 Hz(8倍軸頻),365.0 Hz(12倍軸頻),415.2 Hz(在泵噴處出現最大變形的軸系彎曲模態),517.0 Hz(17倍軸頻)和547.5 Hz(18倍軸頻)。推力軸承橫向傳遞力和垂向傳遞力具有類似的規律,主要的峰值頻率為30.4 Hz,118.6 Hz,121.7 Hz,152.1 Hz(5倍軸頻),243.3 Hz,304.1 Hz(10倍軸頻),365.0 Hz,416.7 Hz(在泵噴處出現最大變形的軸系彎曲模態)和517.0 Hz。除了軸頻,4個最大的響應峰值出現在121.7 Hz,243.3 Hz,365.0 Hz和415.2 Hz;小的峰值出現在273.7 Hz,334.6 Hz(11倍軸頻)。推力軸承縱向傳遞力的幅值比橫向和垂向傳遞力大一個數量級。

圖10 推力軸承和艉后軸承的力傳遞Fig.10 Force transmission by thrust bearing and stern tube bearing

對于艉后軸承傳遞力,可得到類似結論。垂向傳遞力的顯著峰值出現在30.4 Hz,118.6 Hz,121.7 Hz,152.1 Hz,182.5 Hz(5倍軸頻),212.9 Hz(7倍軸頻),243.3 Hz,304.1 Hz,365.0 Hz,416.7 Hz和547.5 Hz。橫向傳遞力的顯著峰值出現在30.4 Hz,118.6 Hz,121.7 Hz,152.1 Hz,182.9 Hz,212.9 Hz,243.3 Hz,304.1 Hz,365.0 Hz,416.7 Hz和517.0 Hz。在艉后軸承處的垂向和橫向傳遞力比推力軸承處的垂向和橫向傳遞力大一個數量級。

在上述峰值頻率中,大部分為泵噴推進器脈動壓力的軸頻及其倍頻,以及結構的模態,如艇體一階縱振模態、軸系的一階和二階縱振模態,但是‘Hub’無法激勵起豐富的脈動壓力峰值頻率和模態頻率。

3.4 水下聲輻射結果與分析

泵噴-軸系-艇體耦合系統的水下輻射聲功率(sound radiation power,SRP)和均方振速(mean-square velocity,MSV)級,如圖11所示。不同的載荷施加方式下輻射聲功率和均方振速的變化規律和傳遞力類似。可見轉子上分布式脈動壓力引起的振動聲輻射將占主導作用,主要的峰值頻率出現在30.4 Hz,91.3 Hz(3倍軸頻),121.7 Hz,152.1 Hz,182.5 Hz(6倍軸頻),209.9 Hz,212.9 Hz,243.3 Hz,273.7 Hz,304.1 Hz,334.6 Hz(艇體一階縱振),365.0 Hz,395.4 Hz(13倍軸頻),416.7 Hz,456.2 Hz(15倍軸頻),486.6 Hz(16倍軸頻),517.0 Hz,547.5 Hz和577.9 Hz。 一些小的峰值出現在209.9 Hz,235.7 Hz和252.4 Hz。并且在艇體一階縱振模態334.6 Hz處的響應比相鄰峰值均大。傳遞力中在軸系一階彎曲模態118.6 Hz出現的峰值在輻射聲功率和均方振速中未得到明顯體現,這是由于軸系彎曲模態通過艉后軸承傳遞的力以及通過推力軸承傳遞到艙壁的力都難以有效激勵起艇體的振動聲輻射。推力軸承傳遞力在軸系一階縱振模態209.9 Hz出現的峰值在輻射聲功率和均方振速中仍具有較大的峰值。為進一步研究不同方向脈動壓力激勵起的艇體聲輻射,取出泵噴的縱向脈動壓力施加到轉子、定子和泵噴表面,將得到的輻射聲功率和同時在泵噴表面施加三向脈動壓力的結果進行比較,如圖11(c)所示,可見在2倍軸頻以上,聲輻射主要由縱向脈動壓力主導。

圖11 輻射聲功率和均方振速Fig.11 Sound radiation power level and mean-square velocity level

給出4個典型頻率30.4 Hz(軸頻),243.3 Hz(8倍軸頻),273.7 Hz(葉頻),334.6 Hz(11倍軸頻)的聲輻射指向性,如圖12所示。場點位于水平面內,是以原點為中心直徑50 m的圓。軸頻處指向性為橢圓,首尾方向的輻射噪聲較大;8倍軸頻處為蝴蝶形,首尾方向和正橫方向的輻射聲壓較小。除葉頻外,其他3個頻率處‘Pumpjet’和‘Rotor’基本重合。葉頻處各種載荷施加方式的指向性差異較大,這是由于聲輻射同時包含了艇體、泵噴轉子、泵噴定子和導管的輻射聲壓,導管本身脈動壓力在葉頻處較大,也產生直接聲輻射,并且和轉子、定子等輻射聲壓疊加,而其他施加方式下導管的振動較小。

圖12 4個頻率處的指向性圖Fig.12 Sound radiation directivity pattern at four frequencies

為衡量導管和定子上的脈動壓力對水下聲輻射的貢獻,僅將導管和定子上的脈動壓力施加到結構網格上,計算輻射聲功率級和均方振速級,結果如圖13所示。由圖13可知,導管和定子上的脈動壓力對艇體水下聲輻射的主要貢獻體現在葉頻和2倍葉頻處,這是由于導管上的脈動壓力在這兩個頻率處為顯著峰值。另一方面,結合圖9的模態結果可知,與這兩個激勵頻率接近的結構模態也被激發,出現峰值。分別在整個泵噴推進器表面、只有轉子表面、只有導管-定子表面施加脈動激勵力,計算得到各自的輻射聲功率和均方振速,并進行比較。可知,在單獨導管-定子脈動壓力激勵下,存在一些突出線譜峰值,量級和轉子脈動壓力激勵下的輻射聲功率和均方振速量級可比擬,但它們對整體聲輻射的貢獻可以忽略。在兩者同時激勵時,均方振速和輻射聲功率級并沒有顯著增大。這是由于在導管-定子和轉子脈動壓力同時激勵下,相位差會導致相互抵消。為驗證抵消效應,取出艇體和導管上某監測點的法向速度響應,如圖14所示。由圖14可知,在整個泵噴推進器脈動壓力激勵時,法向速度會比導管-定子、轉子分別施加脈動壓力時小,這表明,同時施加脈動壓力會由于相位差相互抵消,不同位置處的相互抵消情況會有所差別。

圖13 導管-定子、轉子分別施加脈動壓力時響應Fig.13 SRPand MSV under distributed pressure applied on duct-stator,rotor separately

圖14 導管-定子、轉子分別施加脈動壓力時法向振速響應Fig.14 Normal velocity under distributed pressure applied on duct-stator,rotor separately

4 結 論

開展了suboff艇后泵噴推進器脈動壓力激勵下泵噴-軸系-艇體耦合系統的振動聲輻射響應計算和分析。利用CFD模型獲得了艇體泵噴推進器的脈動壓力。采用耦合有限元法和邊界元法獲得了分布式脈動壓力激勵下泵噴-軸系-艇體耦合系統的傳遞力、聲輻射響應。得到結論如下:

(1)泵噴-軸系-艇體耦合系統的振動聲輻射以轉子脈動壓力激勵下的響應占主導,而其中2倍軸頻以上以縱向脈動壓力的貢獻占主導。

(2)泵噴-軸系-艇體耦合系統振動聲輻射響應中,顯著的峰值包含泵噴脈動壓力的軸頻及其倍頻成分(以軸頻、葉頻最大),也包含軸系一階縱振和二階縱振、艇體縱向振動等結構模態。

(3)導管的脈動壓力以葉頻和2倍葉頻為主,導管-定子脈動壓力激勵下振動聲輻射響應峰值出現在上述頻率處;且響應量級和轉子表面脈動壓力激勵處于同一量級;但導管-定子、轉子表面同時施加脈動壓力時,由于相位相互抵消,均方振速和輻射聲功率級并沒有顯著增大。

主站蜘蛛池模板: 国产97视频在线| 美女毛片在线| 亚洲毛片网站| 国产95在线 | 欧美精品1区2区| 国产地址二永久伊甸园| 亚洲中文字幕23页在线| 国产一在线| 亚洲成人一区二区| 色婷婷亚洲综合五月| 国产99在线| 啦啦啦网站在线观看a毛片| 国产精品免费久久久久影院无码| 毛片免费视频| 国产专区综合另类日韩一区| 欧美啪啪网| 一区二区理伦视频| 国产精品爆乳99久久| 国产在线观看精品| 九九久久精品免费观看| 成人看片欧美一区二区| 美女裸体18禁网站| 青青青国产免费线在| 国产亚洲精久久久久久久91| 日韩无码视频播放| 国产午夜一级毛片| 午夜成人在线视频| 高清无码不卡视频| 成人韩免费网站| 九九热精品视频在线| 亚洲日本中文字幕乱码中文| 伊人久久福利中文字幕| 亚洲成人一区二区三区| www.youjizz.com久久| 亚洲国产成人久久精品软件| 欧美色视频日本| 一区二区三区成人| 日韩欧美国产精品| 精品1区2区3区| 成人va亚洲va欧美天堂| 国内a级毛片| 72种姿势欧美久久久大黄蕉| 免费激情网站| 亚洲国产清纯| 国产精品久久久精品三级| 99尹人香蕉国产免费天天拍| 亚洲中文字幕无码爆乳| 亚洲综合日韩精品| 欧美 亚洲 日韩 国产| AV在线天堂进入| 日韩第一页在线| 亚洲成综合人影院在院播放| 天天摸天天操免费播放小视频| a天堂视频| 51国产偷自视频区视频手机观看| 亚洲最新网址| 毛片大全免费观看| 免费Aⅴ片在线观看蜜芽Tⅴ| 欧美日韩久久综合| 综合成人国产| 国产精品久久精品| 狠狠亚洲五月天| 美女无遮挡免费视频网站| 国产成在线观看免费视频 | 亚洲AV无码精品无码久久蜜桃| 国产情精品嫩草影院88av| 日韩高清一区 | 国产肉感大码AV无码| 国产欧美日本在线观看| 亚洲无码高清一区| 久久黄色影院| 91最新精品视频发布页| 国产日韩欧美一区二区三区在线| 亚洲中文无码h在线观看| 亚洲三级色| 凹凸精品免费精品视频| 毛片视频网| 亚洲最大综合网| 国产高清色视频免费看的网址| 青草娱乐极品免费视频| 欧美日韩高清在线| 毛片久久久|