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

聚焦波作用下立柱高階波浪力特性研究

2020-10-27 01:03:20呂合媛盧文月
海洋工程 2020年5期
關鍵詞:優化

呂合媛,盧文月,李 欣

(上海交通大學 海洋工程國家重點實驗室,上海 200240)

海洋環境千變萬化,使得環境載荷十分復雜。隨著近些年來海洋環境不斷惡化,極端海況頻繁出現,波浪砰擊、高頻共振響應(ringing)等現象成為了海洋工程領域內非常受關注的話題。在極端波浪作用下,由于波浪及其相互作用具有非線性特性,結構物上的波浪荷載由入射譜峰值頻率附近的線性諧波分量和峰值頻率倍數的高階諧波分量組成[1],這些現象具有很強的高頻特性。海洋結構物多種多樣,無疑都會遭遇惡劣海洋環境,即面對強非線性條件,結構物所遭受破壞程度也會因波浪情況的不同而不同。關于海洋中波浪載荷的理論分析工具能較好的預測波浪力的前兩階成分,但對于嚴峻海況下的高頻成分預測結果與實際相差較大[2-3]。所以對總波浪力高階成分的研究是非常需要且有意義的。

最初關于波浪高頻特性的研究來源于對ringing現象的研究。ringing現象是多樣化的,會出現在錨鏈、平臺立柱、船體等部位,不斷有學者對實際工程中出現的相關現象提出分析模型。挪威石油管理局(NPD)也曾匯編了一份文件總結挪威運營商和一家挪威工程公司對于高階錨鏈響應分析的觀點,Natvig[4]在其基礎上進行深入的研究,提出新的分析模型,指出ringing產生過程中出現的非高斯波可能是誘發其產生的一重要因素。Davies等[5]對陡波作用下的模型和實尺度固定結構物和浮動張力腿平臺的ringing現象及其載荷特性進行了總結。Grue等[6]也通過不同尺度下的試驗研究,總結了各波浪參數對波浪力的影響,詳細解釋了高頻諧波成分出現的原因。Chaplin等[7]通過重復試驗,分析了作用于同一陡波的垂直圓柱在固定和底部由彈簧支撐兩種狀態下的試驗結果,得出了結構物的瞬態劇烈響應來源于波浪的高階成分的結論。另有相關試驗證實ringing現象就是高階諧波對結構物瞬態劇烈響應的激發,且通常是三階和頻分量[2, 6]。

極端波浪的激勵和響應等特性一直以來都是海上工程領域備受關注的話題,為了充分捕捉結構物在極端波浪中的總力和響應,必須分離并精確計算高階波載荷諧波,由此眾多學者開始展開一系列對于波浪高頻特性的研究。關于提取規則波高階諧波的研究多是利用傅里葉變換,但對于不規則波,傅里葉變換顯然不再適用[1]。Ma等[8]應用平均小波能量來分析諧波下波浪力的幅值,展開波浪水槽中垂直圓柱在聚焦波作用下的試驗,用平均波能分析了聚焦位置處確定頻率范圍內力和波高的前六階諧波分量,結合波能分析對六階分量進行對比研究,結果表明,該方法適用于瞬態或短時時空序列的諧波分析。Zang等[9]提出了一種適用于陡波作用下圓柱所受載荷分析的新方法,通過假設自由表面高度和聚焦波水平載荷在頻率和波陡上均存在廣義斯托克斯諧波序列,提取相互作用中較高的諧波,將一波群與其逆波群分別作用得到的力的時歷相結合,來提取總水動力載荷中的完整諧波結構,并通過垂直圓柱波浪試驗證明了此方法的優越性。此方法被稱作“逆相位分解方法”,后面也將利用此方法展開研究。

1 聚焦波理論模型

Chaplin[10]總結了實驗室內通過聚焦手段生成瞬態波的三種優化方式,分別為相速度法、逆傳播法和波群速法,對于聚焦波的模擬,普遍采用相位聚焦法,在物理水池或數值水池中使波浪能量在特定的地點與時刻匯聚,生成能量高度集中的強非線性波浪。由于不同波頻或波數的前進波在水中傳播時存在傳播速度不同的色散現象,在聚焦波的模擬中,需要使波浪能量定點匯聚,使各成分波在某一位置同時達到最大值,疊加形成聚焦波列。隨機波浪可以表示為無數個隨機余弦波浪的線性疊加:

(1)

其中,N為組成波數目,an、kn、ωn、εn分別為第n個組成波的波幅、波數、圓頻率和隨機相位,為了使所有組成波在特定位置xc及特定時刻tc同時達到最大值,隨機相位εn可以設為εn=knxc-ωntc。這樣,波浪的分布可以描述為:

(2)

采用新波理論推演聚焦波群來模擬海洋環境中的極端波浪。Tromans等[11]提出了能夠模擬給定波峰值的瞬態波“新波”理論。該理論用來模擬給定目標波峰高度的聚焦波,優點為不受組成波數目的影響,通過實現設定的目標波峰高度,并將波浪譜能量按照相關比例分配到不同的頻率區間內,進而求得各組成波波幅的高度,無論組成波數目多少,在目標位置處聚焦波波浪的波峰值均為事先指定值,且該方法不需要很長的隨機波列[12-13]。新波理論模型是將波面看做一個高斯隨機過程,把歷史記錄中最不利的波面情況通過數學方法表達出來。單位波峰的標準型波面可以表達為波面升高的自相關函數形式,即:

(3)

式中:δω為能量譜的頻率間隔,S(ωn)是組成波對應的波浪譜,σ2是波浪譜S(ωn)的方差。通過ηN(t)=ηNΔr(t)可以獲得與給定波峰值一致的瞬態波波浪升高的時歷表達。波峰波谷的不對稱性依賴于所選取波浪譜的形態。每個組成波的波幅值按照下述公式求解:

(4)

則單向入射波聚焦模型可以表示為:

(5)

圖1 局部波浪參數描述Fig. 1 Description of local wave parameters

聚焦波是一種瞬時大波,屬于非平穩過程,采用局部波浪參數對其進行描述。如圖1所示,聚焦波的波峰值Hc為靜水面到最大波峰高度的垂向距離;波谷Ht為峰前波谷到靜水面的距離;波峰前后波谷之前的時間間距為聚焦波的谷谷周期Ttt,上升周期和下降周期分別為波面從峰前波谷到波峰、從波峰到峰后波谷的時間,聚焦波的理論局部波數k和局部波長L均通過谷谷周期計算得到,波陡定義為波峰與谷谷周期的比值Hc/Ttt;KR為散射參數,定義為立柱等效半徑R與波長L的比值;KC數定義為πA/R。

2 數值模型及準確性驗證

2.1 模型參數

數值計算波浪水池模型如圖2所示,水池長10 m,高1.25 m,初始自由液面水深為0.9 m,自由液面以上高度為0.35 m,在距離水池入口5 m處放置虛擬浪高儀,水池尾端1.5 m為“消波區域”,結構物放置于水池5 m處。數值計算假設流體不可壓縮、黏性系數為常數,控制方程為連續性方程、N-S方程和自由液面控制方程。數值波浪水池入口端邊界條件設為速度入口,另一端出口設為壓力出口,底部與立柱表面設為壁面(無滑移),水池兩側面設為對稱面。

圖2 數值波浪水池Fig. 2 Numerical wave tank

參照真實海洋平臺所處海域情況設計工況,采用波浪參數為有效波高Hs=13.4 m,譜峰周期Tp=14.7 s,峰高因子γ=2.4的JONSWAP譜,并以200的縮尺比進行換算。結構物參照真實海洋平臺立柱以200的縮尺比換算得來,為一個橫斷截面為正方形的立柱,表1給出了結構物的具體尺寸。

表1 結構物尺寸Tab. 1 Size of structure

選擇STARCCM+軟件中的VOF波模塊造聚焦波,是一種純數值造波技術,需要在波浪入口處導入各個成分波浪的參數,包括每個成分波的波幅、周期和相位信息。成分波的參數參照之前所述方法計算得來。數值水池“消波區域”通過在邊界施加控制方程,整個流域按照此方程進行控制。

為了得到效率與質量較高的數值模擬計算,進行數值計算之前需要進行網格和時間步長的敏感性分析。經驗證,在純波浪作用下,水池寬度對波浪的計算沒有影響,為提高計算效率,波浪校核階段建立“偽二維”數值波浪水池,即將水池寬度設置為很小量(y方向),這里設置為0.001 m。在數值波浪水池中,從自由液面附近到整個計算域進行了五次網格密度的過渡,表2給出了網格劃分的區域位置,以及每塊區域內網格的尺寸為基準值的百分比。分別設置基礎尺寸為0.003 m、0.006 m、0.009 m、0.012 m、0.015 m來測試網格的敏感性,由圖3可以看出基礎尺寸越小計算精度越高,0.003 m、0.006 m的網格計算精度均較高且相差不大,考慮到0.006 m的基礎網格計算時間短于0.003 m,所以最終選擇0.006 m的網格進行計算。

表2 網格劃分方案Tab. 2 Mash scheme

圖3 不同基礎尺寸下模擬對比Fig. 3 Simulation comparison under different foundation sizes

2.2 波浪校核

利用新波理論推演出一系列單向入射的聚焦波群,聚焦位置為距離入口5 m處,聚焦時間為5 s,波幅為1的標準聚焦波型,選擇0.02 m,0.04 m,0.06 m,0.07 m作為波峰值來制定工況,表3為工況參數。

表3 工況參數Tab. 3 Parameters of cases

尤其是以聚焦波為代表的極端波浪的傳播過程是非線性的,在聚焦目標位置處,各組成波的相位和幅值與目標值不可避免會存在差別。為了提高聚焦波的模擬質量,采用相位—幅值優化方法對波浪信號進行優化。通過比較聚焦目標位置測量譜與目標譜的差距,調整波浪輸入信號的相位及幅值信息[14]。迭代優化過程如下:

φi+1=φi+(φT-φMi)

(6)

(7)

其中,φi為初始輸入的波浪信號或上一次輸入的波浪信號的相位角,φi+1為新的信號相位角,φT為目標波浪時歷組成波成分的相位角,φMi為試驗得到的波浪時歷組成波的相位角,a為波浪組成成分的幅值,其上下標具有與φ相同的含義,迭代優化過程中均是多退少補的過程,盡管在迭代優化過程中高階波浪成分均作為自由組成波處理,但是由于每次優化都是基于上一次模擬結果進行,因而相位—幅值優化方法內在的考慮了波浪非線性效應。

圖4 工況2優化結果Fig. 4 Optimization results of NO.2 work conditions

圖5 工況2優化后模擬波浪譜與波浪靶譜對比 Fig. 5 Comparison diagram of simulated wave spectrum and wave target spectrum after optimization of NO.2 work conditions

作為示例,圖4為工況2的優化過程,圖5為工況2優化后模擬波浪譜與波浪靶譜對比。

圖6為各個工況波浪優化結果,圖6(a)為設定4個工況的優化結果,圖6(b)為設定工況所對應增加180°相位角的波列優化結果。從中可以看出相位—幅值優化方法使得聚焦波最終模擬結果與目標時歷更為接近,優化效果較好。0.02 m,0.04 m,0.06 m校核結果相比0.07 m效果較好,原因為波幅為0.07 m工況波陡較大,會出現波浪破碎現象,校波過程較為困難。

圖6 工況1~4優化結果Fig. 6 Optimization results NO.1~NO.4 work conditions

2.3 模型準確性驗證

將數值水池寬度擴大至1 m與模型試驗波浪水槽一致。計算立柱的三種固定方式在聚焦波作用下所受波浪力,三種固定方式分別為:完全固定、單自由度縱蕩、單自由度垂蕩。立柱的運動設置通過在模擬樹中添加DFBI模塊,在此模塊中選擇所涉及的運動方式來完成。數值計算涉及網格運動,應用重疊網格法進行建模。重疊網格部分的網格參數按照自由液面附近網格尺寸設置,初始化后,重疊網格與數值水池網格交界處會自動過渡處理,圖7為算例初始化后三維網格模型。

圖7 三維網格模型Fig. 7 Three-dimensional mesh model

選擇工況2進行模型準確性驗證,圖8為固定立柱在物理試驗與數值模擬下所受縱向與垂向波浪力對比,圖9為單自由度垂蕩運動狀態下物理試驗與數值模擬中立柱所受縱向與垂向波浪力及垂蕩運動對比。可以看出各結果的整體趨勢一致,但數值有略微差別,分析此差別是由校波結果與試驗環境引起的,可以得出此數值計算模型具有準確性。

圖8 固定立柱在物理試驗與數值模擬下所受x與z方向波浪力對比Fig. 8 Comparison of wave forces in x and z directions on fixed columns in physical test and numerical simulation

圖9 垂蕩運動狀態下物理試驗與數值模擬中所受x與z方向波浪力及垂蕩運動對比Fig. 9 Comparison of wave forces in x and z directions and heave motion under heave motion condition in physical test and numerical simulation

3 結果分析

3.1 逆相位分解法

逆相位反演法是基于斯托克斯理論呈現的一種通過改變波群相位并對時歷結果進行線性組合來分離高階波浪力的方法[1]。根據經典的斯托克斯波理論,結構物水平水動力載荷為:

F=A1f11cosφ+A2(f20+f22cos2φ)+A3(f31cosφ+f33cosφ)+A4(f42cosφ+f44cosφ)

(8)

其中,Ai為波浪力中第i階幅值,系數fmn代表波浪與波浪力之間的傳遞函數,φ=ωt+φ0為入射波成分的相位。逆相位分解法分離波浪力中奇次諧波和偶次諧波,對于給定的波浪條件,需要兩個來自相差180°相位角的反向入射波時歷。對于初始選擇的信號,會生成一列在聚焦位置出現波峰的波浪時歷,因此反向信號會生成一列在聚焦位置出現波谷的波浪時歷。此時的水動力載荷為:

Fπ=A1f11cos (φ+π)+A2[f20+f22cos2(φ+π)]+A3[f31cos (φ+π)+f33cos (φ+π)]+

A4[f42cos (φ+π)+f44cos (φ+π)]

(9)

對兩列時歷進行基本運算可以得到:

(10)

(11)

即可將波浪力中奇偶次成分分離。

3.2 波浪力分離處理

對數值模擬得到的波浪力時歷按照逆相位分解法進行處理,首先對對應兩列時歷進行峰值位置核對,來消除時間上的偏差引起的誤差,然后按照上述方法將對應兩列數據加和處理得到偶次諧頻波浪力,差減處理得到奇次諧頻波浪力,最后分別還原各階波浪力對應的時歷。圖10~21為工況1、2縱向與垂向波浪力、高頻分離結果以及各階波浪力時歷(圖中e為實際波浪力,1為一階波浪力,2為二階波浪力,3為三階波浪力,4為四階波浪力)。工況3、4的結果不再進行展示。

圖10 工況1,固定狀態下Fx、高頻分離結果及各階波浪力時歷Fig. 10 NO.1, Fx, high-frequency separation results and wave force of each order in a fixed state

圖11 工況1,固定狀態下Fz、高頻分離結果及各階波浪力時歷Fig. 11 NO.1, Fz, high-frequency separation results and wave force of each order in a fixed state

圖12 工況1,單自由度縱蕩運動下Fx、高頻分離結果及各階波浪力時歷Fig. 12 NO.1, Fx, high-frequency separation result and wave force of each order under single freedom of surge motion

圖13 工況1,單自由度縱蕩運動狀態下Fz、高頻分離結果及各階波浪力時歷Fig. 13 NO.1, Fz, high-frequency separation result and wave force of each order under single freedom of surge motion

圖14 工況1,單自由度垂蕩運動狀態下Fx、高頻分離結果及各階波浪力時歷Fig. 14 NO.1, Fx, high-frequency separation result and wave forces of each order under single freedom of heave motion

圖15 工況1,單自由度垂蕩運動狀態下Fz、高頻分離結果及各階波浪力時歷Fig. 15 NO.1, Fz, high-frequency separation result and wave forces of each order under single freedom of heave motion

圖16 工況2,固定狀態下Fx、高頻分離結果及各階波浪力時歷Fig. 16 NO.2, Fx, high-frequency separation results and wave force of each order in a fixed state

圖17 工況2,固定狀態下Fz、高頻分離結果及各階波浪力時歷Fig. 17 NO.2, Fz, high-frequency separation results and wave force of each order in a fixed state

圖18 工況2,單自由度縱蕩運動狀態下Fx、高頻分離結果及各階波浪力時歷Fig. 18 NO.2, Fx, high-frequency separation result and wave force of each order under single freedom of surge motion

圖19 工況2,單自由度縱蕩運動狀態下Fz、高頻分離結果及各階波浪力時歷Fig. 19 NO.2, Fz, high-frequency separation result and wave force of each order under single freedom of surge motion

圖20 工況2,單自由度垂蕩運動狀態下Fx、高頻分離結果及各階波浪力時歷Fig. 20 NO.2, Fx, high-frequency separation result and wave forces of each order under single freedom of heave motion

圖21 工況2,單自由度垂蕩運動狀態下Fz、高頻分離結果及各階波浪力時歷Fig. 21 NO.2, Fz, high-frequency separation result and wave forces of each order under single freedom of heave motion

由以上波浪力結果顯示此次數值模擬計算較為穩定,且每組對應波列作用下的波浪力結果有良好的對稱性。工況4為所有工況中最陡的一個,波浪優化結果較差于其他工況,其波浪力結果相比于前三個工況較為不穩定,兩列波浪力對稱性略差,可以推斷波浪優化校核結果對波浪力結果有一定的影響,也可以認識到對于較陡波浪情況下數值模擬的計算不易控制,對于較陡波浪數值模擬的準確性有待繼續考究。

由高頻分離結果可以發現奇次成分中一階以上成分相較于一階成分非常小,在奇階次圖中幾乎難以發現,偶階次圖中可以清晰分辨出二階及更高階成分,則說明高階成分相比于一階成分為非常小量。同時可以發現二階成分前仍有一段頻率較低的成分存在,這是由于逆相位分解法推演過程中二階成分中含有和頻與差頻成分(式(11)),圖中頻率較低部分為差頻成分[12]。

縱觀數據處理結果,各工況下單自由度縱蕩運動狀態下各階波浪力均小于固定狀態,且高階力明顯較小;單自由度垂蕩運動狀態下縱向波浪力與固定狀態相差甚微,而垂向各階波浪力明顯高于固定狀態,且高階成分突出。

由各階波浪力時歷圖可以觀察到,高階波浪力時歷主要集中分布在波浪聚焦發生前后,越接近聚焦時間,波浪力越大,并且從工況1至工況4,高階成分逐漸變得明顯突出,說明高階波浪力大小與波浪自身的條件有關,能量越集中的地方就更易出現高階成分,高階波浪力的出現使得波浪力峰值更為尖銳,這會導致結構物受力變大、增速較大,結構物受力過于突然有可能會導致非常嚴重的后果。

3.3 各階波浪力統計值

圖22 各階波浪力平均值統計Fig. 22 Statistics of the average values of wave forces of each order

圖中可以看出一般情況下波浪力階次越高平均值越小,偶爾會出現高階波浪力大于低階波浪力的情況。

圖22(a)、(c)、(e)為各狀態下各階縱向波浪力平均值統計,對于縱向波浪力,固定狀態下各階值均為最大;其次是單自由度垂蕩運動狀態,此狀態一階波浪力與固定狀態下相差甚小,但高階力明顯小于固定狀態;單自由度縱蕩運動狀態下各階均值最小,且高階力有明顯的下降趨勢。可以得出固定結構物承受的縱向載荷更為嚴重;結構物的縱向運動會一定程度上削減縱向各階波浪力成分,且高階成分明顯減小;結構物的垂蕩運動對縱向波浪力一階成分影響較小,高階成分明顯減小。

圖22(b)、(d)、(f)為各狀態下各階垂向波浪力平均值統計,對于垂向波浪力,單自由度垂蕩運動狀態下各階均值最大;其次是固定狀態;單自由度縱蕩運動下各階均值最小,且高階力有明顯的下降趨勢。聚焦波在工況4波幅為0.07 m,第二節波浪校核時效果較差,出現過波浪破碎現象,此時呈現出該工況下各階波浪力不穩定。可以得出結構物的縱蕩運動可以一定程度上削減垂向各階波浪力成分,垂蕩運動反而增加了垂向各階波浪力成分。

整體看來單自由度縱蕩運動對縱向與垂向各階波浪力均有緩和作用,且縱向波浪力各階成分有明顯減小;而單自由度垂蕩運動對縱向各階波浪力影響較小,對垂向各階波浪力產生了加劇作用,且高階分量有明顯增加,進一步增加了垂向波浪力。

4 結 語

基于立柱在聚焦波作用下的數值計算,講述了一種利用改變波浪信號來分離垂直立柱所受高階力的方法。介紹了聚焦波的模擬與分離高階波浪力的具體理論方法,以STARCCM+軟件為工具,分別計算了立柱在固定、單自由度縱蕩和單自由度垂蕩三種運動狀態下所受波浪力,詳細講述了高階波浪力分離方法的應用。經系統的研究,得出以下結論:

1) 波浪較穩定的狀態下,數值模擬也有較好的穩定性,每組對應波列作用下的波浪力有良好的對稱性。

2) 聚焦波模擬優化過程中利用相位—幅值優化法,在波浪穩定的狀態下,取得了良好的優化結果,但對于波陡較大、波浪狀態嚴峻的情況,數值模擬與物理試驗經對比更易于破碎,優化過程難以實現,導致與目標值相差較大。

3) 逆相位分解方法可以有效分離波浪力的高階成分。

4) 文中分別計算了結構物三種狀態下所受波浪力,不同運動狀態下立柱所受波浪力有顯著不同。整體看來,固定結構物承受的縱向波浪載荷最嚴重,單自由度垂蕩運動結構物承受的垂向波浪載荷最嚴重;單自由度縱蕩運動對縱向與垂向各階波浪力均有緩和作用,且縱向波浪力各階成分有明顯減小;而單自由度垂蕩運動對縱向各階波浪力影響較小,對垂向各階波浪力產生了加劇作用,且高階分量有明顯增加,進一步增加了垂向波浪力。

研究過程中不難發現,數值計算在模擬真實波浪時仍存在一些不足,對于較難問題仍需繼續研究,來提高數值計算的準確性。

猜你喜歡
優化
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
PEMFC流道的多目標優化
能源工程(2022年1期)2022-03-29 01:06:28
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
圍繞“地、業、人”優化產業扶貧
今日農業(2020年16期)2020-12-14 15:04:59
事業單位中固定資產會計處理的優化
消費導刊(2018年8期)2018-05-25 13:20:08
4K HDR性能大幅度優化 JVC DLA-X8 18 BC
幾種常見的負載均衡算法的優化
電子制作(2017年20期)2017-04-26 06:57:45
主站蜘蛛池模板: 97人人做人人爽香蕉精品| 日本精品影院| 97超级碰碰碰碰精品| 免费高清a毛片| 91色老久久精品偷偷蜜臀| 91成人精品视频| 国产精品久久久久无码网站| 国产色婷婷视频在线观看| 99热线精品大全在线观看| 毛片免费高清免费| 国产成人亚洲综合a∨婷婷| 亚洲男人在线天堂| 国产视频只有无码精品| 55夜色66夜色国产精品视频| 亚洲福利视频网址| 国产91色在线| 久久综合丝袜长腿丝袜| 玖玖精品视频在线观看| 真实国产精品vr专区| 亚洲乱码视频| 成年人视频一区二区| 久久久久国产精品嫩草影院| 国内丰满少妇猛烈精品播| 国产日韩欧美一区二区三区在线| 国产欧美日韩综合在线第一| 亚洲男人的天堂久久香蕉网| 久久久波多野结衣av一区二区| 国产精品国产三级国产专业不 | 国产一级精品毛片基地| 69av免费视频| 久久亚洲高清国产| 欧美一区二区啪啪| 国产日韩久久久久无码精品| 丁香婷婷在线视频| 亚洲国产天堂在线观看| 国产自视频| 国产精品免费p区| 高清无码一本到东京热| 久久久久国产一区二区| 欧美成人第一页| 久爱午夜精品免费视频| 亚洲人成人无码www| 亚洲人成网站在线观看播放不卡| 成人精品亚洲| 喷潮白浆直流在线播放| 成年女人a毛片免费视频| 区国产精品搜索视频| 高清无码不卡视频| 亚洲AⅤ综合在线欧美一区| 欧美三级视频网站| 国产一区二区网站| 亚洲色无码专线精品观看| 亚洲男人在线天堂| 人妻少妇久久久久久97人妻| A级毛片高清免费视频就| jijzzizz老师出水喷水喷出| 国产精品护士| 国产精品三级专区| 亚洲色成人www在线观看| 亚洲人成网站色7799在线播放 | 日韩高清在线观看不卡一区二区| 亚洲国产综合精品一区| 免费无遮挡AV| 老司机久久精品视频| 激情网址在线观看| 亚洲欧美综合精品久久成人网| 美女国产在线| 欧美有码在线| 宅男噜噜噜66国产在线观看| 麻豆AV网站免费进入| 日韩毛片基地| 国产一在线| 麻豆AV网站免费进入| 国产欧美日韩va另类在线播放| 国产成人综合亚洲欧美在| 91无码网站| 午夜日韩久久影院| 高潮爽到爆的喷水女主播视频| 国产精品内射视频| 国产亚洲精品97在线观看| 国产精品爆乳99久久| 国产97视频在线观看|