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

波浪與帶窄縫多箱體作用共振現象的模擬研究

2015-06-01 12:30:12寧德志蘇曉杰滕斌
海洋學報 2015年3期

寧德志,蘇曉杰,滕斌

(1.大連理工大學海岸和近海工程國家重點實驗室,遼寧大連 116024)

波浪與帶窄縫多箱體作用共振現象的模擬研究

寧德志1,蘇曉杰1,滕斌1

(1.大連理工大學海岸和近海工程國家重點實驗室,遼寧大連 116024)

針對波浪與帶有窄縫多箱體結構作用產生的流體共振問題,建立了基于域內源造波技術的二維非線性時域數值波浪水槽模型,其中自由水面滿足完全非線性運動學和動力學邊界條件,窄縫內流體引入人工阻尼來等效由于渦旋運動和流動分離引起的黏性耗散,計算域邊界采用高階邊界元進行離散。通過模擬三箱體間兩窄縫內相對波高變化,并與已發表的數值與實驗結果對比,驗證了本模型的準確性。同時通過大量的數值計算,分析了箱體數量對窄縫內水體共振頻率、共振波高以及對結構反射波高和透射波高的影響。

窄縫;流體共振;域內源造波技術;高階邊界元;非線性數值波浪水槽

1 引言

為了充分利用海洋空間資源,由多模塊組成的超大型浮體已在海洋工程領域得到利用[1],譬如用作海上機場、軍事基地、儲存器和海難救助點等,國內外學者已經開展了很多有關超大型浮體水動力特性的研究,如波浪沖擊[2],水彈性運動響應[3],非穩定外荷載作用下的動力響應[4]等。由于組成超大型浮體的各個模塊之間并不是無縫連接,大都存在相對模塊特征長度很小的窄縫,縫隙內的水體在某些頻率波浪作用下會發生共振現象,誘發很大的波浪爬高和荷載,進而對海洋結構的作業安全帶來很大的影響。

有關海洋結構物間窄縫內波浪水動力特性問題,國內外已經開展了許多相關理論、模型試驗和數值模擬研究,但大多數研究目前還只是局限于兩個結構或3個結構間縫隙水體共振問題。譬如,Miao等[5]采用漸近匹配法研究了帶狹縫二維雙箱的共振現象,給出了狹縫很小時雙箱的理論共振頻率。Saitoh等[6]對不同入射波浪作用下兩個方箱間窄縫的波高變化進行了試驗研究,發現窄縫內最大共振波高可以達到入射波高的5倍。Kristiansen和Faltinsen[7—8]對波浪與具有窄縫的二維固定及浮式箱體和固定式岸壁結構之間的相互作用分別開展了數值和試驗模擬,發現線性和非線性勢流數值結果都比共振條件下的窄縫內試驗波高要大,尤其是線性結果更大;并研究了共振條件下浮體3個自由度的運動響應與水動力系數。在三體結構雙縫隙方面,Iwata等[9]分別對不同入射波浪作用下3個方箱之間窄縫的波高變化進行了試驗研究,發現共振頻率與箱體的吃水深度、窄縫寬度和箱體個數成一定的函數關系。何廣華等[10]采用比例邊界有限元方法分別研究了波浪與三箱結構作用下窄縫內水體共振現象。Lu等[11—12]采用黏性流模型和改進的線性勢流模型對三體結構間窄縫內水體共振引起的波浪爬高和波浪荷載進行了研究,并發現在窄縫內自由水面上布置一定的人工阻尼時,勢流模型也可以得到與試驗和黏性流模型相一致的結果。通過上述研究發現,三箱體間窄縫內水體共振頻率明顯不同于雙箱浮體情況,甚至出現多個自振頻率等復雜情況,也進一步說明箱體數量對窄縫內水體運動特性具有重要影響,而目前對于3個及以上浮體間窄縫內水體共振問題的研究還很少。

本文將在前人研究基礎上,建立自由水面滿足完全非線性邊界條件的二維時域數值波浪水槽模型,采用域內源造波方法產生入射波浪并布置前置阻尼層消除二次反射影響,實現在較小計算域內進行長時間模擬;參照Lu等[11]的方法在窄縫水體自由水面上引入常人工阻尼系數來等效共振條件下黏性耗散。通過模擬波浪與3箱體間雙窄縫內水體運動問題,并與Iwata等[9]實驗結果及Lu等[12]的CFD數值結果進行對比,證明了本模型的準確性。進而通過大量數值計算研究了箱體數量(最多至6個)對窄縫內水體共振頻率、共振波高以及對結構迎浪側反射波高和背浪側透射波高的影響,總結一般性規律,為如超大型浮體這樣由多模塊組成的具有窄縫的海洋結構水動力分析提供參考。

2 數學模型

考慮單向規則波浪與具有窄縫的多固定箱體相互作用問題,其布置如圖1所示。建立二維笛卡爾坐標系Oxz,z=0位于靜水面上,且z軸向上為正,x軸右方向為正。計算域包含自由水面Гf和固體邊界ГN(包括水底Гd和箱體Гb)。波浪由控制垂直源造波面(Гs)的流量密度產生。圖中h為水槽靜水深,W為箱體寬度,D為箱體吃水深度,Wg為兩箱體間窄縫的寬度,N為箱體個數,本研究中取值為2、3、4、5、6。考慮問題一般性,本文中假定各個箱體的寬度、吃水深度和縫隙寬度均一致,圖1中迎浪側e點、背浪側f點和各個箱體件窄縫(gap)內波高變化規律是本文研究的重點。在流體無黏、不可壓縮和流動無旋的假定下,整個流域的速度可用速度勢的梯度來描述,上述問題的控制方程為由速度勢滿足的泊松(Poisson)方程[13],即

式中,q*(xs,z,t)=2Vδ(x-xs)為造波源強度;造波位置x=xs(本文均取xs=0);V為流體質點水平速度,本文給定二階Stokes速度解析解。

在自由水面上,滿足完全非線性動力學和運動學邊界條件,本文中采用混合歐拉-拉格朗日方法更新自由水面,利用物質導數,并在計算域的上游和下游區域的自由水面分別布置人工阻尼層來吸收從結構物反射回來的波浪與出流波浪[14],在窄縫內布置一常參數人工阻尼來近似由于渦和分流引起的黏性耗散[15]。自由水面邊界條件可以寫成以下形式:

式中,η代表自由水面的鉛垂位移,g是重力加速度,X0=(x0,0)是指水質點初始靜止時的位置,阻尼系數

用于計算域邊界的兩個阻尼層,x1和x2分別是左右兩側阻尼層的起點位置;L為阻尼層長度,本文取1.5倍波長(即1.5λ);ω是波浪角頻率。阻尼系數μ2用于窄縫內自由水面,其數值根據試驗黏性耗散來確定;k是波數,滿足如下線性色散方程關系

在固定的結構表面和底面邊界上,流體法向速度為0,滿足固壁不可滲透邊界條件。由于本研究在時域內進行,自由水面滿足靜初始條件,即起始時刻速度勢和波面均為0。

對于滿足上述控制方程和邊界條件的定解問題,在整個流域內對速度勢應用格林第二定理,可轉換如下邊界積分方程[13]:

式中,p=(x0,z0)為源點;q=(x,z)為場點;C為固角系數;Ω代表整個流域;G是簡單格林函數,考慮到水底鏡像,可以表示為如下形式:

式中,r1為p和q兩點距離,

r2為p和q關于水底鏡像之間距離,

本文用三節點高階邊界元離散計算域成一些曲線單元,單元內任一點的幾何坐標和速度勢等物理量可以用如下二次形狀函數hi(ξ)插值得到,

式中,ξ代表固有坐標,取值范圍(-1.0,1.0)。

這樣在任一邊界單元內,物理量和幾何量都可以通過形狀函數插值得到,也即,邊界單元都是等參的。積分方程(4)經高階邊界元離散后,可以表示成如下形式:

式中,Ne1、Ne2、Ne3分別為自由水面、物面和出入流邊界及造波源面上劃分的單元個數,J(ξ)是聯系大地坐標和固有坐標的雅可比行列式。

由于本文采用了高階邊界元方法,自由水面邊界條件式(2)中用到的速度勢的空間導數項進而可以表示成如下形式:

式中,n=(nx,nz)為單位法向量,指出水體為正。

最后把未知量都移到方程(9)的左側,該方程組就可以寫出如下矩陣乘積形式:

式中,X是未知的速度勢和速度勢法向導數,A為空間系數矩陣,F為由已知的速度勢和速度在邊界上積分得到的列向量。

計算中認為當前時刻物面上的速度勢法向導數和自由水面上的速度勢是已知的,根據積分方程計算當前時刻物面上的速度勢和自由水面上的速度勢法向導數,然后應用四階Runga-Kutta法,根據自由水面條件式(2)計算下一時刻的水質點位置和自由水面上的速度勢,再用二次形狀函數在舊單元上插值求得新節點上的物理量來對自由水面網格重新劃分,重新應用積分方程計算下一時刻物面上的速度勢和自由水面上的速度勢法向導數。這樣計算周而復始,直到計算結束[16—17]。

3 數值計算及討論

3.1 模型準確性和穩定性

作為算例,本文以Iwata等[9]的實驗來進行數值模擬波浪與具有窄縫的三箱體的相互作用的研究,驗證本文數學模型的準確性。這里選用的實驗參數為水槽靜水深h=0.5 m,箱體的寬度W=0.5 m,吃水深度D=0.252 m,入射波高H0=0.024 m,箱體間窄縫寬度0.05 m。在數值模型中,計算域長度取7.5倍波長,在水槽的左右兩端各布置1.5倍波長的阻尼層,造波源位于x=0,箱體1左側面邊界位于距離造波源2.5倍波長的位置,然后依次按Wg調整箱體2和3的位置。通過開展數值收斂性實驗,自由水面上每個波長布置15個單元,窄縫內布置2個單元,計算域垂向邊界和造波源面分別布置10個單元,箱體側面邊界上均布置6個單元,底面邊界上均布置12個單元;時間步長Δt=T/60 s,每個算例模擬30個周期。

圖2給出了窄縫寬度Wg為0.05 m情況下兩個窄縫中心位置無量綱波高Hg/H0隨入射波波數kh的變化關系,及本文數值結果與實驗數據[9]、黏性流模型數值結果[12]的比較。由于窄縫內水體黏性與波浪條件無關,故在某一波浪頻率下通過取一系列不同人工黏性系數進行數值模擬,并與實驗數據進行對比,確定窄縫中間水面的人工黏性系數取為μ2=0.03;數值模擬中,波高Hg以測試點的穩定時間段內波面時間序列中波峰值與波谷絕對值和的平均來計算。從圖中可以看出,箱體1與箱體2之間窄縫內數值模擬波面高在kh=1.35時達到最大,也即窄縫內流體發生共振,波高為入射波高的4.6倍。箱體2與箱體3之間窄縫內數值模擬波面高在kh=1.4時達到最大,波高為入射波高的4.4倍。同時發現窄縫1內出現2個共振頻率,這與單縫隙內的共振規律是不同的,也進一步說明浮體數量對窄縫內水體運動規律的重要影響。整體上兩種數值結果與實驗數據均符合的很好,在個別位置處甚至本文結果比黏性流模型結果與實驗數據吻合的更好,說明所建立模型在取得合適的人工黏性系數情況下可以準確模擬多箱體窄縫內流體共振問題。

圖3給出了波數kh=1.40,窄縫寬度Wg=0.05 m,入射波高H0=0.024 m情況下t=26T和30T時的整個計算域波面分布,圖中虛線分界處分別對應3個箱體所在的位置。從圖中可以看出,兩個時刻的波面曲線已經完全重合,包括箱體1前的反射波,箱體3后的透射波和3個箱體之間的窄縫內波面;并且兩端的波面基本趨于0,說明兩端阻尼層吸收波浪的效果很理想;箱體1前的波浪形成了穩定的立波,說明其反射回去的波浪透過造波源被前端阻尼層完全吸收,而對入射波浪沒有產生影響。以上現象說明本模型的模擬結果已達到穩定。

圖2 兩窄縫中無因次波高Hg/H0與波數kh間的關系Fig.2 Distribution of non-dimensional wave height versus incident wave number at two narrow gaps

圖3t=26T和30T兩個時刻的水槽波面分布Fig.3 Snapshot of wave elevation along the wave flume att=26Tand 30T

3.2 數值結果

下面仍以上述工況為例,保持各參數不變,只是改變箱體數量,進一步分析箱體數量對窄縫內水體共振頻率、波浪爬高以及結構前反射波高和結構后透射波高的影響。

圖4和圖5是箱體個數N分別為4和5時各個窄縫內無因次波高Hg/H0隨入射波波數kh的變化關系。各窄縫的主共振頻率相仿,其所對應的共振波高在處于中間位置的窄縫內達到最大。在大于主共振頻率的某一頻率處,還會有次共振現象發生,且處于兩端位置的窄縫內較中間位置的窄縫內次共振現象更為明顯,尤為突出的是迎浪側的窄縫內,在圖5的窄縫1中高共振頻率所對應的共振波高已大于其低共振頻率所對應的波高。其原因可能是各窄縫內自振頻率下的波浪會向兩側傳播并與迎浪側透射的波浪相互作用,處于中間位置的窄縫由于對稱關系,從窄縫兩側反射來的波浪疊加達到最大,而處于邊界的窄縫由于其相對位置關系則會發生復雜的波浪干涉現象。

圖4 4箱體時各窄縫內波高隨波數的變化關系Fig.4 Dimensionless wave height againstkhat various gaps for four boxes

圖6和圖7分別為不同數量箱體時窄縫1(即迎浪側第一個窄縫)和窄縫N-1(背浪側第一個窄縫)內無因次波高隨波數的分布情況。可以看出隨著箱體數量的增加,主頻共振發生時,窄縫內波高減小,對應的共振頻率向低頻偏移。并且箱體數量較多時,窄縫會在多個頻率發生共振,且高共振頻率處的波高隨箱體個數的增加而增大,窄縫1內在箱體數N為5和6時已大于低共振頻率處的波高,窄縫5內高頻共振波高也大于低頻共振波高。各個工況下窄縫N-1內主共振頻率都不小于窄縫1內主共振頻率,而窄縫1和N-1內相鄰箱體數對應的次頻共振頻率的差值近似常數(Δkh≈0.1)。

圖5 5箱體時各窄縫內波高隨波數的變化關系Fig.5 Dimensionless wave height againstkhat various gaps for five boxes

圖6 窄縫1內無因次波高隨波數的變化關系Fig.6 Dimensionless wave height againstkhat gap 1

圖7 窄縫N內無因次波高隨波數的變化關系Fig.7 Dimensionless wave height againstkhat gapN

表1和表2列出了不同箱體數情況下各個窄縫內對應共振頻率和共振波高分布情況。可以看出,除了得到與圖4、圖5、圖6和圖7相同的規律外,還會發現對于同一箱體個數N,各窄縫內的低共振頻率會從迎浪側向背浪側逐漸增大,對應的共振波高則是中間窄縫內最大,向迎浪側和背浪側兩側減小;而高共振頻率則是處于中間位置窄縫內最大,向迎浪側和背浪側兩側減小,對應的共振波高則是中間位置窄縫內最小,向兩側增大。

表1 各窄縫對應的低共振波數和共振波高(kh,Hg/H0)Tab.1 Low resonant wave number and wave height at various gaps(kh,Hg/H0)

表2 各窄縫對應的高共振波數和共振波高(kh,Hg/H0)Tab.2 High resonant wave number and wave height at various gaps(kh,Hg/H0)

圖8 迎浪側e點波高隨波數的變化關系Fig.8 Wave height at the weather-side pointeagainstkh

圖8和圖9分別給出了迎浪側e點和背浪側f點無因次波高隨波數的變化關系。從圖8中可以看出,在窄縫1內主頻共振頻率處,迎浪側e點波高出現最小值,且箱體數N越大,對應的極值越小(N=6時,最小的Hg/H0=1.4),在低頻處無因次波高接近于2,相當于波浪線性化并被全反射,而在高頻處無因次波高會大于2,這是由于波浪的非線性增強,高階諧波貢獻所致。從圖9可以看出,在窄縫N內主頻共振頻率處,背浪側f點波高也出現明顯的峰值,但其變化規律與圖8相反,隨著箱體數N的增大,透射波高的峰值也增大(N=6時,最大的Hg/H0=0.78),在高頻處,透射波高接近于0。

圖9 背浪側f點波高隨波數的變化關系Fig.9 Wave height at the lee-side pointfagainstkh

4 結論

本文基于域內源造波的時域高階邊界元方法建立波浪與具有窄縫的多箱體結構相作用的完全非線性數值水槽模型,對不同個數箱體時窄縫內流體的共振頻率、共振波高及結構迎浪側波高和背浪側透射波高等進行了模擬研究。通過與已發表實驗數據和數值結果進行對比驗證,表明本文所建立數學模型可以準確模擬波浪與具有窄縫的多箱體相互作用過程,且在較小計算域內可長時間模擬得到穩定的結果,沒有在入射邊界發生二次反射現象。研究發現:箱體寬度、吃水、窄縫寬度一定時,不同數量箱體,窄縫內水體發生共振的頻率是不同的且箱體個數較多時會有多個共振頻率。隨著箱體數量的增加,主頻共振發生時,窄縫內波高減小,主頻共振頻率值減小;次頻共振頻率處的波高隨箱體個數的增加而增大。對于一給定箱體數N的情況,各窄縫內的主頻共振頻率會從迎浪側向背浪側逐漸增大,對應的共振波高則是中間窄縫內最大,向迎浪側和背浪側兩側減小;而次共振頻率則從處于中間位置窄縫內最大,向迎浪側和背浪側兩側減小,對應的共振波高則是中間位置窄縫內最小,向兩側增大。迎浪側波高在窄縫1內共振發生時最小,背浪側波高則在窄縫N-1內共振發生時出現峰值;隨著箱體個數的增加,共振發生時迎浪側波高減小,背浪側波高增大。

[1] 繆國平,劉應中.征服海洋之夢——超大型海洋浮式結構物[J].自然雜志,1996,18(1):26-30.

Miao Guoping,Liu Yingzhong.A dream to conquer the ocean——super large floating ocean structures[J].Ziran Zazhi,1996,18(1):26-30.

[2] Yashimoto H,Ohmatsu S,Hoshino K,et al.Slamming load on a very large floating body with a shallow drift[J].Journal of Marine Science and Technology,1997,2:163-172.

[3] Kagemoto H,Fujino M,Murai M.Theoretical and experimental predictions of the hydro-elastic response of a very large floating structure in waves[J].Applied Ocean Research,1998,20:135-144.

[4] Qiu L,Liu H.Three-dimensional time-domain analysis of very large floating structures subjected to unsteady external loading[J].Journal of Offshore Mechanics and Arctic Engineering,2006,129(1):21-28.

[5] Miao G P,Ishida H,Saitoh T.Influence of gaps between multiple floating bodies onwaveforces[J].China Ocean Engineering,2000,14(4):407-422.

[6] Saitoh T,Miao G P,Ishida H.Theoretical analysis on appearancecondition of fluid resonance in a narrow gap between two modulesof very large floating structure[C]//Proceedings of the Third Asia-Pacific Workshop on Marine Hydrodynamics,Shanghai,China,2006:170-175.

[7] Kristiansen T,Faltinsen O M.Studies on resonant water motion between a ship and a fixed terminal in shallow water[J].Journal of Offshore Mechanics and Arctic Engineering,2009,131:021102.

[8] Kristiansen T,Faltinsen O M.A two-dimensional numerical andexperimental study of resonant coupled ship and piston-modemotion[J].Applied O-cean Research,2010,32:158-176.

[9] Iwata H,Saitoh T,Miao G P.Fluid resonance in narrow gaps ofvery large floating structure composed of rectangular modules[C]//Proceedings of the Fourth International Conference on Asian and Pacific Coasts,Nanjing,China,2007:815-826.

[10] 何廣華,滕斌,李博寧,等.應用比例邊界有限元法研究波浪與帶狹縫三箱作用的共振現象[J].水動力學研究與進展,2006,21(3):418-424.

He Guanghua,Teng Bin,Li Boning,et al.Research on the hydrodynamic influence from the gaps between threeidentical boxes by a scaled boundary finite element method[J].Journal of Hydrodynamics,2006,21(3):418-424.

[11] Lu L,Cheng L,Teng B,et al.Numerical investigation of fluid resonance in two narrow gaps of three identical rectangular structures[J].Applied Ocean Research,2010,32(2):177-190.

[12] Lu L,Teng B,Cheng L,et al.Modelling of multi-bodies in close proximity under water waves——fluid resonance in narrow gaps[J].Science China Physics,Mechanics&Astronomy,2011,54(1):16-25.

[13] Ning D Z,Teng B,Eatoack Taylor R,et al.Numerical simulation of non-linear regular and focused waves in an infinite water-depth[J].Ocean Engineering,2008,35(8/9):887-899.

[14] Tanizawa K.Long time fully nonlinear simulation of floating body motions with artificial damping zone[J].The Society of Naval Architects of Japan,1996,180:311-319.

[15] Kim Y W.Artificial damping in water wave problems:Ⅰ.Constant damping[J].International Journal of Offshore and Polar Engineering,2003,13(2):88-93.

[16] 周斌珍,寧德志,滕斌.造波板運動造波實時模擬[J].水動力學研究與進展,2009,24(4):1-12.

Zhou Binzhen,Ning Dezhi,Teng Bin.Real-time simulation of waves generated by a wave maker[J].Journal of Hydrodynamics,2009,24(4):1-12.

[17] 陳麗芬,寧德志,滕斌,等.潛堤上波流傳播的完全非線性數值模擬[J].力學學報,2011,43(5):834-843.

Chen Lifen,Ning Dezhi,Teng Bin,et al.Full nonlinear numerical simulation for wave-current propagation over a submerged bar[J].Chinese Journal of Theoretical and Applied Mechanics,2011,43(5):834-843.

Numerical study of fluid resonance induced by wave action on multi-boxes with narrow gaps

Ning Dezhi1,Su Xiaojie1,Teng Bin1

(1.State Key Laboratory of Coastal and Offshore Engineering,Dalian University of Technology,Dalian 116024,China)

Based on wave generation technique by a inner-domain source,a two-dimensional nonlinear numerical wave flume over time is developed to investigate the fluid resonance induced by the interaction between wave and multi-objects with narrow gaps.In the numerical model,the fully nonlinear kinematic and dynamic boundary conditions are set for the instantaneous free surface;the artificial damping is introduced into the gap to simulate the viscous dissipation due to vortex motion and flow separation;the computational domain is discretized using higher-order boundary elements.The proposed model is validated by the published experimental and numerical data of the relative wave height at two narrow gaps of three boxes.Numerical experiments are performed to study the following:the effects of the number of the boxes on the resonant frequency,wave height at various gaps and the reflected and transmitted wave heights of the objects.

narrow gap;fluid resonance;inner-domain source generation technique;higher-order boundary element;nonlinear numerical wave flume

O353.2

A

0253-4193(2015)03-0126-08

寧德志,蘇曉杰,滕斌.波浪與帶窄縫多箱體作用共振現象的模擬研究[J].海洋學報,2015,37(3):126—133,

10.3969/j.issn.0253-4193.2015.03.013

Ning Dezhi,Su Xiaojie,Teng Bin.Numerical study of fluid resonance induced by wave action on multi-boxes with narrow gaps[J].Haiyang Xuebao,2015,37(3):126—133,doi:10.3969/j.issn.0253-4193.2015.03.013

2014-01-22;

2014-04-21。

國家自然科學基金項目(51179028,51222902,51221961);教育部新世紀優秀人才支持計劃(NCET-13-0076);中央高校基本科研業務費專項資金(DUT13YQ104)。

寧德志(1975—),男,黑龍江省五常市人,教授,博士生導師,主要從事非線性波浪及其與結構物相互作用的研究。E-mail:dzning@dlut.edu.cn

主站蜘蛛池模板: 国产成人超碰无码| 亚洲第一黄色网址| 欧美日本在线| 欧洲亚洲一区| 97国产精品视频自在拍| 91精品国产福利| 国产97视频在线| 精品国产成人av免费| 亚洲欧美日韩中文字幕一区二区三区| 国产精品亚洲五月天高清| 免费观看无遮挡www的小视频| 欧美精品1区2区| 日韩乱码免费一区二区三区| 免费无码一区二区| 久久久久国色AV免费观看性色| 四虎永久免费地址| 久久久国产精品免费视频| 欧美日韩一区二区三| 国产成人精品一区二区三在线观看| 超碰精品无码一区二区| 精品久久人人爽人人玩人人妻| 国产成人综合欧美精品久久| 在线免费亚洲无码视频| 99热最新网址| 婷婷亚洲视频| 伊人色在线视频| 免费a级毛片视频| 色成人亚洲| 丝袜国产一区| 精品天海翼一区二区| 欧美一级大片在线观看| 国产sm重味一区二区三区| 91成人在线免费观看| 国产尤物在线播放| 国产a在视频线精品视频下载| 欧美特黄一级大黄录像| 久久鸭综合久久国产| 69免费在线视频| 亚洲日韩高清在线亚洲专区| 国产在线自乱拍播放| 欧美视频在线播放观看免费福利资源| 亚洲欧美自拍中文| 免费又黄又爽又猛大片午夜| 精品午夜国产福利观看| 91亚洲国产视频| 久久人人妻人人爽人人卡片av| 九九这里只有精品视频| 亚洲国产中文在线二区三区免| 最新国产午夜精品视频成人| 香蕉国产精品视频| 国产欧美精品专区一区二区| 国产日韩精品一区在线不卡| 精品欧美视频| 中日韩一区二区三区中文免费视频| 欧美人与性动交a欧美精品| 国产美女在线免费观看| 亚洲美女操| 尤物精品视频一区二区三区 | 国产一级二级在线观看| 久久精品人人做人人爽97| a亚洲视频| 欧美无专区| yjizz国产在线视频网| 少妇精品网站| 亚洲欧洲日韩久久狠狠爱| 91香蕉视频下载网站| 中文字幕免费播放| 国产精品视频猛进猛出| 日韩成人午夜| 亚洲国产欧美国产综合久久| 国产精品久久久久久久久kt| 婷婷亚洲天堂| 欧美a在线看| 精品国产一二三区| 欧美国产在线一区| 国产成人在线小视频| 色天天综合| 理论片一区| 欧美精品xx| 视频二区国产精品职场同事| 日韩精品久久无码中文字幕色欲| 欧美特级AAAAAA视频免费观看|