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

流體黏性對受迫振蕩水平圓柱受力的影響

2020-07-01 05:06:58袁劍平毛鴻飛赫巖莉潘新祥賈寶柱
廣東海洋大學學報 2020年4期

袁劍平,毛鴻飛,赫巖莉,潘新祥,賈寶柱

流體黏性對受迫振蕩水平圓柱受力的影響

袁劍平1,毛鴻飛1,赫巖莉1,潘新祥2,賈寶柱2

(1. 廣東海洋大學海洋工程學院,廣東 湛江 524088;2. 廣東海洋大學海運學院,廣東 湛江 524088)

【目的】準確預測顯著黏性效應下結構物所受流體作用力。【方法】采用有限體積(FVM)方法構建兩相流數值波浪水槽模型。首先,通過半潛水平圓柱受波浪力數值結果與前人實驗數據的對比,驗證數值模型的準確性。其次,對不同流體黏性下,受迫振蕩淹沒水平圓柱上的流體作用力進行數值計算,通過黏性流結果和勢流結果的對比,分析流體黏性對圓柱受力的影響特征。最后,通過將黏性流和勢流理論壓力結果對比,并結合對圓柱周圍黏性渦量場分布特征分析,解釋不同流體黏性下,黏性流與勢流結果之間出現差異的原因。【結果】在不同流體黏性條件下,圓柱所受各倍頻流體作用力均出現一定程度上的差別,表明流體黏性的影響明顯。【結論】相比勢流理論,采用黏性流理論對受迫振蕩水平圓柱受力的影響的研究更為準確。

圓柱;受迫振蕩;流體作用力;數值模擬;波浪水槽

在流體與結構相互作用的流體動力學領域研究中,近自由面淹沒水平圓柱在水中受迫振蕩問題受到了關注和研究。基于勢流理論展開的相關研究為了簡化計算,將結構物上的流體作用力分解為固定物體上的波浪繞射力和振蕩運動結構物上的輻射力。Wu[1]基于勢流理論,考慮線性自由面條件,對近水面淹沒水平圓柱受迫振蕩受力問題進行了理論解析,并將圓柱受力分解到了不同倍頻上。Tyvand等[2]基于勢流理論,考慮線性自由面條件,對近水面水平圓柱初始入水沖擊問題進行了理論解析,分析了圓柱受力特征和自由面波動特征。上述基于勢流理論的研究僅考慮線性自由面條件,考慮非線性自由面條件的相關研究也相繼開展。Wu等[3]、Liu等[4]分別基于勢流理論,考慮非線性自由面條件,對近自由面淹沒水平圓柱受迫振蕩問題進行了計算,考察了輻射波浪分布情況,分析了輻射波的非線性。Kent等[5]基于勢流理論,考慮了非線性自由面條件,對淹沒水平圓柱垂向受迫振蕩的受力問題進行了數值研究,發現對于三倍頻垂向力,其結果區別于Wu[1]的解析結果。Guerber等[6-7]基于勢流理論建立了二維非線性數值水槽模型,并對淹沒水平圓柱受迫振蕩受力問題進行了數值計算研究,研究發現非線性自由面條件對二倍頻力有影響。

由于上述研究基于勢流理論忽略了流體黏性,無法對圓柱受力給出真實準確的計算結果;而僅考慮線性自由面條件也顯然把問題理想化。因此,本研究有必要使用基于黏性流理論對淹沒水平圓柱受迫振蕩受力問題開展研究。本研究在Teng等[8-9]和毛鴻飛等[10-11]分別對波浪作用于淹沒圓柱結構物和波浪對近水面結構物沖擊作用問題研究的經驗基礎上,采用基于黏性流理論建立的數值波浪水槽模型開展計算研究,模型的基本控制方程為Navier-Stokes方程,方程的數值離散采用有限體積方法,自由面捕捉采用流體體積函數方法,數值造波和消波采用松弛區和速度邊界結合方法,同時可以滿足長時間計算需求。首先,對半淹沒水平圓柱上的垂向波浪力進行了數值計算,與前人實驗數據對比,驗證數值模型對結構物上流體作用力計算的準確性。其次,對受迫振蕩水平圓柱受力進行計算,將流體作用力并分解至多倍頻項,考察其隨振幅變化的特征規律,并將多個黏性流體下的受力結果與勢流結果對比,分析黏性對各倍頻流體作用力的影響規律。最后,將基于黏性流理論和勢流理論得到圓柱上的壓力以及黏性力幅值進行對比,結合圓柱周圍黏性渦量場分布特征的分析,解釋不同流體黏性下流體作用力差別原因。

1 數值模型

采用Navier-Stokes方程對不可壓縮黏性流體的流動求解,其表示為:

為了精確計算湍流流動,模型結合了RNG模型,其對流輸運方程表示為:

式中,為湍動能,為耗散率,1ε、2ε、σσ為模型常數。

對自由面的捕捉采用流體體積函數法,單元密度和黏性系數表示為:

式中,為體積分數,氣相中= 0,液相中= 1,自由面內= 0 ~ 1,ρρ分別為液相和空氣相的密度,μμ分別為液相和空氣相的流體黏性系數。

為提高計算穩定性和界面分辨率,將人工壓縮項引入界面方程,即

數值水槽模型采用速度邊界法造波,結合松弛方法保證波浪穩定傳輸和消除反射,其設置如圖1所示。協助波浪生產、穩定傳播和吸收從結構物的一次反射浪的功能集成于松弛區I中,消除消波出口邊界的二次反射波的功能集成于松弛區II中。

松弛區內,以解析形式修正數值求解,松弛修正需要借助如下松弛函數:

式中,為函數的關系函數,松弛函數的分布如圖1所示。進而,松弛區內的數值修正可表達為:

式中,代表需要修正的物理量,即流體速度和體積分數,角標、分別表示解析值、計算值、目標值[13]。

構建數值波浪水槽模型首先需要對水槽和結構邊界條件進行初始設置:將“速度邊界條件”定義為“入口邊界”,即邊界上的流體速度始終按所需波浪理論給定,初始動態壓強定義為*/= 0;將“可自由進出邊界條件”定義為“頂邊界”,即初始速度定義為?u/= 0,初始動態壓強定義為*= 0;將“不可滑移邊界條件”定義為“出口邊界”、“底邊界”和“水槽中固定結構物邊界”,即初始速度定義為u= 0,初始動態壓強定義為*/= 0。

圖1 邊界條件

為實現強非線性波浪生產和傳播功能,數值波浪水槽模型速度邊界結合了五階Stokes波浪理論控制模塊。波面形式以及水平和垂向速度表達形式分別如式(10)、(11)和(12)所示,

式中,I為松弛區I長度,為階數,=為無量綱波幅,ab為轉化函數,為波數,為波幅,=(?I)?+為相位角,為波浪角頻率,為初相位,B為與水深和波數相關的無量綱參數,0和A為與水深和相關的無量綱參數。上述無量綱參數及相關函數具體表達式可參考Fenton五階Stokes波浪理論[14]。

通過上述數值方法獲得壓強和速度場分布后,結構物上的流體作用力可通過下式進行計算,

2 驗證計算

2.1 計算設置

為了驗證數值波浪水槽模型對結構物受流體作用力計算的準確性,根據Dixon等[15]的實驗室物理模型試驗,并參考前人[16-17]相關數值計算的設置,對半淺固定水平圓柱上垂向波浪力進行計算和對比驗證,本研究的驗證計算設置示意如圖2所示。為水平圓柱半徑,設置為0.125 m;為水深,設置為1.0 m,初始時刻軸心至于靜水面處。速度邊界的造波理論采用五階Stokes波浪理論,波浪周期為1.646 s;波幅設置為0.125 m。圓柱水平位置,松弛區I和松弛區II的長度均根據入射波波長的一定比例進行設置。

圖2 驗證計算示意圖

2.2 波浪力驗證對比

圖3為采用本文數值模型計算獲得的圓柱所受垂向波浪力數值結果與Dixon等[15]實驗數據的對比情況。圖中,無量綱垂向波浪力用z/(π2)表示,時間的無量綱形式表示為/(為入射波周期),考察一個周期內垂向波浪力的時間歷程對比情況。通過觀察可知,本研究數值結果與前人實驗數據總體上吻合良好,這表明,該數值模型對結構物受力計算具有良好的準確性。

圖3 一個周期內圓柱上的垂向波浪力

3 振蕩圓柱上各倍頻流體作用力

3.1 算例設置

本研究的黏性流數值計算根據Wu[1]考慮線性自由面條件的勢流理論的解析研究條件,以及Guerber等[6-7]完全非線性的勢流理論數值模擬研究條件進行設置,計算域如圖4所示。圓柱分別在垂向方向和水平方向受迫振蕩運動,振蕩幅值用表示,圓柱半徑為10 cm,水深設置為4.0 m。圓柱初始淹沒深度,即圓心與靜水面距離為0.3 m,為了消除反射波的影響,將松弛區II設置在水槽兩側,長度取2.5。

圖4 計算域示意圖

3.2 不同黏性下各倍頻流體作用力對比

將本研究計算所得不同黏性下的黏性流數值結果與勢流結果進行對比。用于對比的勢流結果包括前人考慮線性自由面條件的理論解析結果,以及本文考慮非線性自由面條件采用高階邊界元方法計算所得到的數值結果,該數值計算方法參見Zhou等[18]。

為了獲得圓柱上各頻率的流體作用力分量,將流體作用力水平和垂向方向分量按照傅里葉級數形式進行分解,表達為:

式中,(0)為力的均值,(m)為倍頻力,(m)為倍頻分量的相位,下標和分別表示水平和垂向方向分量。考察各倍頻流體作用力隨振幅變化特征,并進行對比。為了對比和分析黏性影響,需排除其他參數的比例影響,將各倍頻流體作用力按照無量綱形式(m)/(π22) 表示,將振幅按照無量綱形式/表示。

圖5為圓柱在振蕩方向上所受一倍頻流體作用力的對比情況。從圖中結果可見,對于小黏性工況,當振幅較小時,黏性流結果與勢流結果比較接近;而當振幅較大時,出現黏性流結果逐漸小于勢流結果的現象,且二者差距逐漸增大。對于大黏性工況,當振幅較小時,黏性流結果明顯大于勢流結果;當振幅較大時,無量綱的黏性流結果隨著振幅增大出現先逐漸減小,而后逐漸增大的現象。兩種不同黏性下的黏性流結果和勢流結果之間均有明顯差別,這說明一倍頻流體作用力受到了流體黏性變化的顯著影響。此外,自由面條件不同的勢流結果非常接近,這表明線性自由面條件和非線性自由面條件對一倍頻流體作用力的影響不顯著。

圖5 振蕩方向上的一倍頻流體作用力

Fig. 5 The first harmonic components of the hydrodynamic forces in the oscillation direction

圖6為圓柱在振蕩方向上所受流體作用力的均值以及高倍頻量的對比情況。從圖中結果可見,相比之下,當振幅較大時,小黏性工況下的黏性流結果和考慮非線性自由面條件的勢流結果總體上更為接近。大黏性工況下的黏性流結果與其他各組結果存在不同程度的差別,這表明,該情況下圓柱振蕩方向上的流體作用力均值和高倍頻量受到流體黏性變化不同程度的影響。此外,當振幅較大時,兩組勢流結果之間也存在一定差別,其原因為自由面條件的不同;也正是因為都考慮了非線性自由面條件,本研究小黏性工況黏性流結果與本文的勢流結果更加接近,相比前人考慮線性自由面條件的勢流結果在大振幅下有偏差。

圖6 振蕩方向上力的均值和高倍頻量

Fig. 6 The mean and the high harmonic components of the hydrodynamic forces in the oscillation direction

圖7為圓柱在水平振蕩情況下所受垂向力均值和二倍頻垂向力的對比。圖中結果可見,當振幅較小時,不同黏性下的黏性流結果與勢流結果均比較接近;而隨著振幅增大,兩次黏性流結果與勢流結果之間的逐漸出現不同程度的差別,相比之下,大黏性工況下的黏性流結果與勢流結果的差別相對更大一些。這表明,當振幅較大時,圓柱水平振蕩情況下所受垂向力均值和二倍頻垂向力均受到流體黏性變化的顯著影響。此外,兩組勢流結果之間的差別原因仍為自由面條件不同。需要注意的是,水平振蕩情況下圓柱所受垂向力均值并不為零,這是由于該情況下輻射波并不以振蕩中心對稱。

圖7 水平振蕩下垂向力均值和二倍頻垂向力

Fig. 7 The mean vertical forces and the second harmonic components of the vertical forces

3.2 不同黏性下各倍頻流體作用力對比

為了解釋不同黏性的黏性流工況與本文考慮非線性自由面勢流工況計算所得圓柱上流體作用力之間差別原因,并揭示黏性影響機理。以圓柱水平方向受迫振蕩,不同黏性條件,/= 0.40,1.25和1.75工況為例,將流體作用力分解為壓力和黏性力,從而進行對比,并對圓柱附近黏性渦量場分布特征進行分析。

圖8為將圓柱所受流體作用力分解后的壓力和黏性力時間歷程對比情況。從對比結果可見,/= 0.40工況下,對于壓力,兩組黏性流結果與勢流結果的幅值很接近,說明流體黏性對壓力影響不顯著。/= 1.25工況下,對于壓力,兩組黏性流結果的幅值均小于勢流結果。/= 1.75工況下,對于壓力,小黏性工況下的壓力幅值小于勢流下的壓力,而大黏性工況下的壓力幅值大于勢流下的壓力。在本文所考察的工況下,對于黏性力,小黏性工況下的黏性力遠小于壓力,相比之下,大黏性工況下的黏性力對總流體作用力的貢獻相對較大。

圖9為水平壓力為負極值時圓柱附近黏性流渦量場。圖中,“Vorticity”為渦量,可表達為=(?????)/(/),其中為流速大小,為圓柱直徑;“F”為慣性力,箭頭指向表示流體相對加速度的方向;“S”和“L”分別代表小黏性工況和大黏性工況,對應的數值為無量綱振幅/值。從圖9-a和9-d可見,兩種黏性小振幅下,圓柱附近無明顯的渦旋產生。從圖9-b、9-c和9-e可見,隨著振幅增大,貼體渦旋出現圓柱加速度方向相反一側。從圖9-f可見,渦旋運動相比圓柱運動滯后性更強,此刻的渦旋位置發生明顯改變。同時,可以判斷,隨著振幅和黏性增大,渦旋壓力在流體作用力中的作用和占比也會發生變化。另外,需要注意的是,圓柱上下渦旋均為明顯的非對稱分布。

圖9 水平壓力為負極值時的圓柱周圍渦量場

為了更直觀說明渦旋對圓柱表面壓力的影響,將/= 0.40和1.75下,水平壓力達到負極值時圓柱表面壓強分布進行對比,如圖10所示。其中,圖10-a為對圓柱表面環向角度的定義,Angel表示環向角度。從圖10-b壓強對比結果可見,曲線的凹陷出現在貼體渦旋對應位置(圖9),且在大黏性工況下,曲線凹陷范圍較大,曲線在水平方向上的積分面積較大。這表明,在小黏性,/= 1.75工況下,由于渦旋壓強始終小于周圍壓強,則其減小了圓柱表面正壓側表面壓強。由于渦旋壓力的幅值相對較小,黏性流下的壓力幅值小于勢流下的壓力幅值,無量綱一倍頻力呈現隨振幅增大而逐漸減小的趨勢(/= 0.80 ~ 1.75);在大黏性,/= 1.75工況下,渦旋壓力幅值相對較大,并在圓柱所受壓力占比中起主導作用,使黏性流下圓柱所受壓力幅值增大至大于勢流結果,因而無量綱一倍頻力呈現隨振幅增大而逐漸增大的趨勢(/= 1.50 ~ 1.75)。

通過上述分析可以得知,對于一倍頻流體作用力,小黏性工況,大振幅下,由于渦旋壓力的存在,圓柱所受壓力幅值相比勢流結果較小,黏性流下一倍頻力隨振幅增大逐漸小于勢流結果,二者之間差距也逐漸增大。大黏性工況,小振幅下,由于黏性力的貢獻相對較大,黏性流一倍頻流體作用力大于勢流結果,而隨著振幅增大,貼體渦旋由起減小圓柱所受壓力作用轉為增大壓力的作用,因而黏性流下一倍頻力無量綱結果呈現先減小后增大的趨勢。對于振蕩方向上的均值力和高倍頻力,大黏性下黏性流和勢流結果之間有明顯差別,這是因為此時黏性力對流體作用力貢獻較大,即黏性力的均值和高倍頻分量對總的流體作用力所造成的影響。水平振蕩下的垂向力均值和二倍頻垂向力,在圓柱大幅水平振蕩情況下,受到自由面波動的影響,圓柱上下所產生的渦旋為非對稱分布,導致了兩組黏性流結果與勢流結果有不同程度的差別。

4 結論

本研究對流體黏性對自由面下受迫振蕩水平圓柱所受流體作用力特征的影響開展了黏性流數值計算研究,研究結論總結如下:

(1)對于振蕩方向上的一倍頻力,兩組黏性流結果與勢流結果均有明顯差別,其原因為大幅振幅下渦旋的產生對于圓柱表面壓力造成的影響,以及大黏性工況下黏性力對流體作用力的貢獻較大。

(2)對于振蕩方向上力的均值和高倍頻力,大黏性工況結果與其他各組結果有明顯差別,其原因為黏性力的貢獻較大,黏性力的均值和高倍頻分量對總流體作用力有影響。

(3)對于水平振蕩下的垂向力均值和二倍頻垂向力,大幅振蕩下,兩組黏性流結果均與勢流結果有差別,并且大黏性工況下差別更明顯,其原因為在自由面影響下,圓柱上下渦旋分布不對稱對圓柱所受垂向壓力有影響。

[1] WU G X. Hydrodynamic forces on a submerged circular cylinder undergoing large-amplitude motion[J]. Journal of Fluid Mechanics, 1993, 254: 41-58.

[2] TYVAND P A, MILOH T. Free-surface flow due to impulsive motion of a submerged circular cylinder[J]. Journal of Fluid Mechanics, 1995, 286: 67-101.

[3] WU G X, EATOCK TAYLOR R. Time stepping solutions of the two-dimensional nonlinear wave radiation problem[J]. Ocean Engineering, 1995, 22(8): 785-798.

[4] LIU Y M, ZHU Q, YUE D K P. Nonlinear radiated and diffracted waves due to the motions of a submerged circular cylinder[J]. Journal of Fluid Mechanics, 1999, 382: 263-282.

[5] KENT C P, CHOI W. An explicit formulation for the evolution of nonlinear surface waves interacting with a submerged body[J]. International Journal for Numerical Methods in Fluids, 2007, 55(11):1019-1038.

[6] GUERBER E, BENOIT M, GRILLI S T, et al. Modeling of fully nonlinear wave interactions with moving submerged structures[C/OL]. Proceedings of the International Offshore and Polar Engineering Conference, 2010, 529-536.

[7] GUERBER E, BENOIT M, GRILLI S T, et al. A fully nonlinear implicit model for wave interactions with submerged structures in forced or free motion[J]. Engineering Analysis with Boundary Elements, 2012, 36(7): 1151-1163.

[8] TENG B, MAO H F, LU L. Viscous effects on wave forces on a submerged horizontal circular cylinder[J]. China Ocean Engineering, 2018, 32(3): 245-255.

[9] TENG B, MAO H F, NING D Z, et al. Viscous numerical examination of hydrodynamic forces on a submerged horizontal circular cylinder undergoing forced oscillation[J]. Journal of Hydrodynamics, 2019, 31(5): 887-899.

[10] 毛鴻飛, 陳洪洲. 非完全淹沒水平圓柱上波浪力特征的數值模擬[J]. 水科學進展, 2019, 30(5): 749-759.

[11] 毛鴻飛, 李芳成, 吳光林,等. 基于黏性流理論對平板受波浪沖擊的兩相流數值研究[J]. 廣東海洋大學學報, 2019, 39(4): 73-80.

[12] WELLER H G, TABOR G, JASAK H, et al. A tensorial approach to computational continuum mechanics using object oriented techniques[J]. Computers in Physics, 1998, 12(6): 620-631.

[13] JACOBSEN N G, FUHRMAN D R, FREDS?E J. A wave generation toolbox for the open-source CFD library: OpenFoam? [J]. International Journal for Numerical Methods in Fluids, 2012, 70(9): 1073-1088.

[14] FENTON J D. A fifth-order Stokes theory for steady waves[J]. Journal of Waterway, Port, Coastal and Ocean Engineering-asce, 1985, 111(2): 216-234.

[15] GRAHAM DIXON A, GREATED C A, SALTER S H. Wave forces on partially submerged cylinders[J]. Journal of the Waterway Port, Coastal and Ocean Division, 1979, 105(4): 421-438.

[16] WESTPHALEN J, GREAVES D M, WILLIAMS C K, et al. Extreme Wave Loading on Offshore Wave Energy Devices using CFD: A Hierarchical Team Approach[C/OL]. Proceedings of European Wave and Tidal Energy. Uppsala, Sweden, 2009.

[17] HU Z, CAUSON D M, MINGHAM C G, et al. Numerical simulation of floating bodies in extreme free surface waves[J]. Natural Hazards and Earth System Sciences, 2011, 11(2): 519-527.

[18] ZHOU B Z, NING D Z, TENG B, et al. The numerical simulation of fully nonlinear deep-water waves[J]. Acta Oceanologica Sinica, 2011, 33(1): 27-35.

Viscous Effects on Hydrodynamic Forces on a Horizontal Circular Cylinder under Forced Oscillation

YUAN Jian-ping1, MAO Hong-fei1, HE Yan-li1, PAN Xin-xiang2, JIA Bao-zhu2

(1.,,524088,;2.,,524088,)

【Objective】To accurately predict the hydrodynamic force of a circular structure under viscous fluid effects.【Method】A viscous fluid numerical wave tank model was established by using the finite volume method. Firstly, the accuracy of the numerical model was validated against available experimental data for the wave forces on a semi-submerged circular cylinder. Secondly, the viscous numerical calculations for different fluid viscosities were carried out. The comparisons of the hydrodynamic forces on the circular cylinder predicted by the viscous fluid model and the potential flow model were conducted to show the viscous effects on the hydrodynamic forces.Finally, by studying the components of the pressure forces and the viscous shear forces, and by analyzing the vorticity fields around the cylinder, the reasons for the differences between the results based on the viscous fluid theory and the potential flow theory were explained.【Result】There are different degrees of differences among the harmonic components of the hydrodynamic forces on the horizontal circular cylinder under different viscosities.【Conclusion】The numerical results predicted by the viscous fluid model are more accurate than those from the potential flow models.

circular cylinder; forced oscillation; hydrodynamic forces; numerical simulation; wave tank

TV139.2+6

A

1673-9159(2020)04-0116-08

10.3969/j.issn.1673-9159.2020.04.016

2020-04-10

廣東省自然科學基金(51979045,51479017);國家裝備預研基金項目(6142204190711);國家國防科工局穩定支持課題(JCKYS2019604SXJQR-02);廣東省教育廳重點領域專項(2019KZDZX1024);廣東省教育廳高校青年創新人才項目(2019KQNCX045);廣東海洋大學科研啟動費項目(120602-R19024);廣東省“沖一流”省財政專項資金建設項目(231419010)

袁劍平(1979-),男,碩士,高級工程師,研究方向為船舶智能控制系統開發。E-mail: yjp_103@163.com

毛鴻飛(1985-),男,博士,講師,研究方向為波浪與結構物的相互作用。E-mail:maohongfei-gdou@qq.com

袁劍平,毛鴻飛,赫巖莉,等. 流體黏性對受迫振蕩水平圓柱受力的影響[J].廣東海洋大學學報,2020,40(4):116-123.

(責任編輯:劉朏)

主站蜘蛛池模板: 国产69囗曝护士吞精在线视频| 成人亚洲视频| 园内精品自拍视频在线播放| 一本大道东京热无码av| 亚洲AV无码精品无码久久蜜桃| 最新国产精品第1页| 欧美亚洲激情| 激情六月丁香婷婷四房播| 成人综合久久综合| 91无码国产视频| 91网址在线播放| 精品精品国产高清A毛片| 国产成人毛片| 国产在线观看一区二区三区| 国产成人亚洲无吗淙合青草| 欧美激情成人网| 99热这里只有精品国产99| 东京热高清无码精品| 欧美无专区| 99久久精品免费视频| 国产最新无码专区在线| 91国内视频在线观看| 婷婷色婷婷| 亚洲有码在线播放| 精品视频在线观看你懂的一区| 51国产偷自视频区视频手机观看| 日韩国产另类| 午夜一级做a爰片久久毛片| 亚洲国产一区在线观看| 久久精品国产亚洲AV忘忧草18| 国产丰满成熟女性性满足视频| www.狠狠| 亚洲一级毛片| 国产激情无码一区二区三区免费| 精久久久久无码区中文字幕| 911亚洲精品| 色综合久久88色综合天天提莫 | 国产熟女一级毛片| 亚欧成人无码AV在线播放| 精品欧美日韩国产日漫一区不卡| 伊人无码视屏| 免费观看无遮挡www的小视频| 亚洲永久视频| 成人一区专区在线观看| 2022国产91精品久久久久久| 幺女国产一级毛片| 国产精品刺激对白在线| 国产经典在线观看一区| 一区二区三区高清视频国产女人| 手机成人午夜在线视频| a欧美在线| 国产成人欧美| 四虎影视无码永久免费观看| 国产精品丝袜在线| 欧美人在线一区二区三区| 狂欢视频在线观看不卡| 亚洲三级a| 伊人国产无码高清视频| 中文字幕欧美日韩高清| 国产免费黄| 亚洲国产理论片在线播放| 在线a视频免费观看| 日韩毛片免费视频| 丁香五月激情图片| 国产成人综合亚洲网址| 国产白浆在线| 国产网友愉拍精品| 伊人91视频| 亚洲乱码视频| 中文字幕日韩欧美| 91久久夜色精品| 国产精品免费电影| 亚洲狼网站狼狼鲁亚洲下载| 亚洲天堂久久久| 99久久99视频| 欧美一区二区福利视频| 中国毛片网| 亚洲无限乱码| 99精品国产自在现线观看| 91色在线视频| 成人一级免费视频| 国产最新无码专区在线|