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

輪船-橋墩碰撞簡化荷載模型

2019-03-25 02:31:48宋彥臣王君杰尹海蛟劉慧杰
振動與沖擊 2019年5期
關鍵詞:船舶橋梁模型

宋彥臣, 王君杰, 尹海蛟, 劉慧杰

(1.同濟大學 土木工程學院,上海 200092;2. 青島市公路規劃設計院有限公司,山東 青島 266100)

船撞安全是跨航道橋梁設計中的一個重要問題,各國橋梁設計規范[1-4]對此都進行了規定。在這些規范中,均將船舶撞擊作用等效為靜力荷載。然而試驗[5-8]和碰撞數值模擬[9-12]表明,船與橋梁構件之間的接觸力接近于一個短時沖擊荷載,引起的動力效應對橋梁結構是重要的,不能被忽略。

沖擊碰撞試驗研究結果表明[13-15],采用數值模擬方法進行船橋碰撞分析是可信賴的。然而,在工程設計中采用復雜的碰撞數值模擬計算,對工程師來說是一個過高的要求和過重的負擔,因此建立考慮沖擊效應的簡化結構分析方法是船撞設計必須解決的問題。

沖擊譜方法是考慮結構沖擊效應的方法之一。Fan等[16-17]先后建立了駁船沖擊譜模型;Wang等[18]基于81條輪船撞擊力時程樣本,建立了輪船沖擊譜模型,并提出了SUM與IFM組合方法。然而,沖擊譜方法僅適用于線彈性結構,并不適用于結構非線性響應分析。

另一個可行的辦法是將船舶等效為一個質點與一個非線性彈簧,橋梁結構仍按照一般的有限元方法建模,即質點碰撞法。該方法最早在丹麥大帶橋[19-20]中提出,但由于非線性彈簧參數的確定缺少數據基礎,影響到其實用性;Consolazio等基于兩艘駁船的有限元模型,獲得了駁船撞擊力-撞深關系,通過與LS-DYNA精細化數值仿真對比,表明質點碰撞法具有很高的求解精度;Fan等[21]基于1艘12 000 DWT輪船的有限元模型,得到了輪船的擬靜力-撞深關系,并進行應變速率效應修正,通過與精細化有限元分析對比,研究結果亦表明質點碰撞模型具有很好的求解精度;Fan等[22]進一步根據兩艘5 000 DWT輪船的有限元模型,分析了被撞擊物厚度對撞擊力-撞深關系的影響。王君杰等[23]根據數值仿真獲得的81條輪船撞擊力-撞深關系樣本,考慮船艏結構的隨機性因素,建立了簡化的輪船撞擊力-撞深關系概率模型。

考慮結構沖擊效應的另一個方法是在橋梁結構上施加一個撞擊力時程,進行一般的動力響應分析。Consolazio等提出了簡化的駁船撞擊力時程模型;樊偉等[24-25]基于輪船船艏撞擊力-撞深關系,通過能量守恒與動量定理,假定撞深時程為二次多項式,得出了輪船撞擊力時程。然而,由于需要獲得輪船的船艏撞擊力-撞深關系,因此該方法不能避免精細化有限元數值仿真,影響到其實用性。歐洲規范建議了一種了簡化的船舶撞擊力時間過程確定方法。該規范將船舶碰撞分為彈性碰撞和非彈性碰撞兩類,并根據撞擊能量確定函數中的參數,但歐洲統一規范沖擊力時間過程的定義是概念性的,缺少必要的數據支撐。

根據上述討論,可以通過建立一定數量的代表性船舶有限元模型,進行碰撞數值模擬得到船舶撞擊力時程樣本,通過數理統計建立船舶撞擊橋梁的簡化沖擊力模型,并確定其參數是一個實用的途徑。

本文建立了5艘代表性輪船的數值模型,通過碰撞數值模擬計算分別獲得了45條撞擊力時間過程樣本。采用修正的半波正弦函數近似描述船舶撞擊力時程。通過統計分析,給出了基于樣本過程確定模型參數的方法,并得到了模型參數。以一座連續橋梁為例進行船撞動力反應分析,將簡化荷載模型的分析結果與碰撞數值模擬的結果進行對比,以分析簡化荷載模型計算結果的精度。

1 基于船-剛性墻碰撞數值模擬的撞擊力樣本過程

1.1 代表船舶

根據中國內河航運的現狀,本文選擇了5艘具有代表性的內河輪船,恒載噸位在500 DWT~5 000 DWT之間,其技術參數見表1。限于篇幅,本文只列出船舶的船艏結構圖,如圖1所示。

(a) 5000 DWT

(b) 3000-2 DWT

(c) 3000-1 DWT

(d) 1000 DWT

(e) 500 DWT

1.2 船-剛性墻碰撞模型

船橋碰撞過程極為復雜,涉及材料非線性、幾何非線性以及接觸非線性等問題,另外船艏結構與形狀、被撞物的形狀與剛度、撞擊速度船舶壓倉量不同均會對碰撞結果產生明顯影響。由于船舶與橋梁之間的相互作用,確定船舶撞擊荷載時應考慮橋梁結構的影響。然而沖擊荷載的上界可以通過假定橋梁為“剛性結構”獲得[26]。根據JCSS這一觀點,當船舶撞擊剛性墻時,應得到最大的撞擊荷載,因此本文的研究中采用了圖2所示的計算模型,即船舶-剛性墻正撞模型。從橋梁設計的角度來看,使用偏保守的荷載也符合橋梁設計的一般原則。

船舶撞擊力與船、橋墩的幾何形狀密切相關,然而綜合考慮這些因素的影響會使問題過于復雜,采用船-剛性墻碰撞計算模型可以使問題得到簡化,這也是采用這種計算模型的原因。

圖2 輪船正撞剛性墻面模型

1.3 船舶有限元模型

如前所述,船舶撞擊力受很多因素影響,這些因素使得船橋碰撞仿真分析較傳統數值仿真分析更為困難。為驗證數值仿真分析的可靠性,作者開展了船舶與鋼箱的縮尺模型沖擊試驗,深入討論了船舶碰撞數值仿真分析的建模方法與參數取值,結果表明采用數值仿真進行鋼箱碰撞分析具有較高的精度;在此基礎上,基于船舶設計標準與圖紙,建立了5艘輪船的精細化有限元模型,如圖3所示。

(a) 5000 DWT

(b) 3000-2 DWT

(c) 3000-1 DWT

(d) 1000 DWT

(e) 500 DWT

在船舶正撞剛性墻過程中,由于船體變形主要集中在船艏,船身耗能較小,在有限元網格方面,網格尺寸從船艏到船身逐漸增大;材料方面,船艏部分采用彈塑性本構模型,為提高計算效率,船身采用剛體材料。船舶鋼板鋼種為Q235,不同構件板厚在10~16 mm,單元類型采用殼單元,船艏網格尺寸約150 mm,如圖3所示。

采用線性隨動強化彈塑性本構關系描述船艏鋼材的非線性變形性質,采用Cowper-Symonds模型[27]考慮應變率效應的影響。線性隨動強化彈塑性本構方程[28]的屈服函數φ為

(1a)

(1b)

對于Q235鋼板,其材料參數取值見表2,其中E表示彈性模量(MPa),ν表示泊松比,其它參數意義同前。剛性墻材料參數按C40混凝土取值,如表2所示。

數值仿真采用顯式動力計算程序LS-DYNA進行。

1.4 計算工況和樣本

針對每艘船舶分別設置了9種撞擊工況,撞擊速度V0從1 m/s到5 m/s,為內河航道船舶常規航速范圍,工況設置如表3所示。

表3 撞擊速度工況設置

通過數值計算,得到了45條撞擊力時程樣本。典型的輪船撞擊力時程如圖4所示。從圖4中可以看出,輪船撞擊力時程整體上可以劃分為3個階段:初始階段,船舶撞擊力迅速上升;中期階段,船舶撞擊力出現波動;卸載階段,船舶與被撞剛性墻分離,撞擊力迅速衰減至零。

圖4 典型輪船撞擊力時程

2 修正半波正弦撞擊力時間過程模型

2.1 標準化撞擊力時間過程

記撞擊力時間過程為F(t),t∈[0,T],T為撞擊力樣本的持續時間。為觀察撞擊時間過程的規律,將時間軸無量綱化,即:

(2)

則撞擊力時間過程可以表示為F(τ),τ∈[0,1]。

(3)

定義無量綱撞擊力系數β(τ),

(4)

在式(3)和式(4)中,I0為撞擊力時程樣本沖量,即,

(5)

根據式(2)~式(4)與45條撞擊力時間過程樣本,45條無量綱系數β(τ)樣本過程如圖5所示。

圖5 無量綱系數樣本β(τ)

2.2 修正半波正弦模型

由圖5可知,樣本撞擊力時間過程復雜,但整體上可以分為三部分:第一部分為上升段,此時船艏主要發生彈性變形;第二部分為波動段,在此階段船艏部分發生明顯的屈曲和壓潰;第三部分為卸載段,在此階段船舶與剛性墻面逐漸分離,船艏少部分的彈性變形得到恢復。簡化沖擊荷載模型應在總體上應描述上述特點。沖擊荷載的常用簡化描述方法[29]包括三角脈沖、半波正弦脈沖以及矩形脈沖,其中半波正弦函數具有上升段與下降段,同時撞擊力的梯度不斷變化,最符合上述要求。由圖5可以觀察到,樣本撞擊力時間過程在中間段多數較為平坦,直接使用半波正弦模型并不適合,因此本文提出一個修正的半波正弦模型來描述撞擊力時間過程,即:

F(τ)=c(1+α1τ+α2τ2)sin(πτ) (0≤τ≤1)

(6)

根據樣本計算得到的沖量應等于由式(6)計算得到的沖量,因此可以得到

(7)

從而式(6)可以重寫為

(8a)

(1+α1τ+α2τ2)sin(πτ) (0≤τ≤1)

(8b)

式中:α1,α2為待定參數,c與α1,α2相關。由式(7)可知,根據式(8)計算得到的沖量恒等于根據樣本得到的沖量,而與α1,α2取值無關。

2.3 修正半波正弦模型的參數

已有研究中一般根據物理條件確定撞擊力時程的模型參數,即采用動能定理與動量定理。然而研究結果表明,當獲得準確的船舶的撞擊力-撞深關系時,這種方法確定的撞擊力時程可達到滿意的求解精度,而采用簡化的撞擊力-撞深關系確定的撞擊力時程,其求解精度較差。此時,可直接采用質點碰撞法計算橋梁結構響應。因此,本文在確定模型參數時未采用這種方法。

數理統計方法是確定模型參數的另一個重要手段,為確定修正半波正弦模型中的參數α1和α2,采用如下兩個原則:①樣本和模型關于力軸的形心位置相等;②樣本和模型關于形心的回轉半徑相等。這兩個原則的數學表達為

(9a)

(9b)

式中,tc表示樣本撞擊力時程的形心位置。式(9)在數學上保證了修正半波正弦撞擊力時程的總體形狀與幅值特征與撞擊力原始樣本保持相同。

根據式(9),可以得到α1和α2表達式如下

α1=

(10a)

α2=

(10b)

式中,A=1-4/π2,B=1-6/π2,τc和η定義如下

(11a)

(11b)

式中,I1表示撞擊力時程樣本對力軸的一次矩,I2表示撞擊力時程樣本對自身形心力軸的二次矩。

由式(3)、式(8)與式(11)可知,修正半波正弦模型需確定T、I0、τc及η的樣本取值。

2.3.1 持續時間T

根據撞擊力原始樣本確定的T樣本如表4所示。

通過觀察T樣本的規律,采用冪函數回歸T與撞擊速度V0的關系

(12)

根據式(12),對T的樣本關于V0進行回歸,γ1與α4的統計值如表5所示。

(13)

式中,Tsta表示T統計值,Tsam表示T樣本值。根據表4與式(12),T的統計誤差如表6所示。根據表6可知,式(12)統計誤差基本在±10%以內,各噸位船舶誤差均值在2%以內,具有較好的統計精度。

2.3.2 撞擊力沖量I0

根據樣本撞擊力時間過程計算得到撞擊力沖量見表7。

將撞擊力沖量I0表示為

I0=α3MV0

(14)

式中,M表示船舶滿載排水噸位(噸)。根據表7,α3的樣本如表8所示,可知α3的離散性很小,在1.03~1.12之間,為便于工程應用,α3按總體樣本均值取為1.07,如表8所示。

2.3.3τc和η

根據式(11),得到81條撞擊力時間過程樣本所對應的τc和η值,列于表9(a)與表9(b)中。

由表9(a)與表9(b),τc和η樣本值的離散性很小,樣本值范圍為[0.427,0.590]與[0.561,0.696],為便于工程應用,τc和η可取各船舶樣本均值,如表9所示。

根據表4、表7與表9,式(8)與式(10),典型的簡化撞擊力時程樣本如圖6所示。

(b) 5000 DWT, 4 m/s

3 簡化荷載模型的精度

3.1 橋梁信息

本文以一座連續梁橋為例,分別進行船橋接觸碰撞分析和基于撞擊力時程的動力響應分析,分析修正半波正弦簡化荷載模型的響應求解精度。

連續梁橋跨徑布置為80 m+140 m+140 m+80 m,主墩基礎采用14根D2.5 m鉆孔灌注樁,樁身混凝土采用C30混凝土,橋梁構造如圖7所示。

圖7 橋梁立面布置與橋墩群樁基礎

3.2 橋梁有限元模型

橋梁結構采用梁單元模擬;為進行接觸碰撞分析,承臺采用實體單元,采用剛體材料,材料參數按C30混凝土取值[30],E=3×104Mpa,ρ=2 600 kg/m3;橋梁上部結構混凝土為C40,參數值[30]:E=3.25×104Mpa,ρ=2 600 kg/m3,主梁在墩頂橫橋向約束;樁基礎采用C30混凝土;場地地質條件為粘土層,將土層分為7層,每層5 m,采用土彈簧模擬土體,考慮動力效應,土體水平抗力系數取m法[31]建議值的2倍[32]。橋梁的有限元模型如圖8所示。

3.3 分析工況

為分析修正半波正弦模型的精度,本文選取了50000 DWT,3000 DWT以及1000 DWT船舶,撞擊速度均為3 m/s,分別進行碰撞反應分析與時程動力響應分析,工況介紹如表10所示,表中含義解釋如下:IMP分析模型表示進行船舶-橋梁接觸響應分析;TH-OS分析模型中撞擊力時程由船舶正撞剛性墻模型得到,即為撞擊力原始樣本;TH-SS表示采用簡化的撞擊力時程,但其參數I0、T、τc和η取其樣本值;而TH-OS也表示采用簡化的撞擊力時程,其參數I0、T、τc和η取回歸統計值。

圖8 橋梁有限元模型

表10 分析工況

3.4 動力響應計算結果

根據表10所列工況,分別進行結構動力反應分析,結果如圖9所示,圖中涉及的符號含義解釋如下:Dp,撞擊向承臺位移;Mpile/Qpile,樁頂彎矩/剪力;Mpier/Qpier,墩底彎矩/剪力。根據圖9,可得出以下主要信息。

(a) 撞擊力

(b) Dp

(c) Mpile

(d) Qpile

(e) Mpier

(f) Qpier

5000 DWT輪船撞擊下,IMP與TH-OS工況撞擊力時程差異很小;相比TH-SE,TH-SS與IMP撞擊力時程的幅值更為接近;在結構響應方面,IMP與TH-OS工況各結構響應之間的差異也很小,各結構響應時程幾乎重合,表明采用剛性橋梁假定用于船橋碰撞分析的合理性;TH-SS與TH-SE工況各響應與IMP工況整體上符合良好,且5000 DWT輪船撞擊時TH-SS的求解精度優于TH-SE。

TH-OS的求解誤差為采用剛性橋梁假定所致;TH-SS的求解誤差在剛性橋梁假定的誤差基礎上,包含了采用修正半波正弦模型帶來的誤差;TH-SE的求解誤差則既包含了上述兩種誤差,同時也包含著模型參數的統計誤差導致的結構響應求解誤差。下文將分三個方面,分別討論剛性橋梁假定、修正半波正弦模型以及參數統計誤差帶來的求解誤差。

3.5 誤差分析

定義指標評價結構響應的求解精度

(15)

式中,rTH表示撞擊力時程模型(TH-OS,TH-SS與TH-SE)的結構響應峰值;rIMP表示實船撞擊下的結構響應峰值。根據表10所列工況,計算結果如表11所示。

根據圖9與表11,分別討論剛性橋梁假定、撞擊力模型化以及模型參數統計誤差對結構響應的影響。

3.5.1 剛性橋梁假設誤差分析

TH-OS具有很好的求解精度,總體誤差均值為2%,表明采用剛性橋梁假定確定船舶撞擊力是合理的。具體來說,承臺位移、樁身內力的求解精度高,最大誤差為2%,如CGB/1000Qpile,;墩底彎矩Mpier、剪力Qpier求解誤差相對稍差,最大誤差為8%,如CGB/1000Qpier。主要原因在于撞擊點以下的結構峰值響應均發生在撞擊持續時間內,此時IMP工況船舶撞擊力與TH-OS極為接近,因而求解精度很高;撞擊點以上的結構峰值響應均發生于撞擊結束時刻,如圖9(e)與圖9(f)所示,由于船舶與結構的相互作用,IMP與TH-OS撞擊力卸載速率不同,從而造成撞擊結束后墩底內力的不同,因而求解精度稍差。

3.5.2 模型化誤差分析

采用修正半波正弦模型描述船舶撞擊力引起的響應求解誤差較小。對比TH-OS與TH-SS可知,各工況由于撞擊力修正半波正弦模型化導致的最大誤差均值在6%以內,如CGB/1000;TH-OS與TH-SS所有工況總體誤差均值相同,表明采用修正半波正弦模型近似描述輪船撞擊力時程樣本是合理的。

3.5.3 參數統計誤差分析

修正半波正弦模型參數統計誤差導致結構響應的誤差在可接受范圍之內。對比TH-SS與TH-SE,模型參數統計誤差導致的結構響應求解誤差均值在11%以內,如CGB/5000;對于CGB/3000與CGB/1000,模型參數統計導致的結構響應誤差較小,誤差均值分別為4%與3%。總體來看,模型參數統計誤差導致的結構響應誤差在可接受范圍之內,因此模型參數的統計值是合理的。

3.6 誤差機理分析

橋梁結構在船舶撞擊作用下的運動方程可以表示為

表11 精度分析

(16)

式中:M、C、K分別表示質量矩陣、阻尼矩陣與剛度矩陣;u表示結構位移系向量;s表示荷載空間分布向量。進行線性變換u=Φq,并對運動方程進行解耦可得

(17)

式中:Φ表示結構振型矩陣,qj表示振型坐標;ξj、ωj、γj分別表示振型阻尼比、振型圓頻率以及振型參與系數,j=1,2,…,N,N表示振型數量。

由式(17)可知,振型參與系數可作為一種參考指標,用于評價結構各階振型參與所關心響應的貢獻大小。通過進行模態分析,所選連續梁橋的振型參與系數如圖10所示。

圖10 振型參與系數

由圖10可知,連續梁橋主要控制振型頻率范圍為0.98~2.5 Hz。根據表11可知,CGB/5000 TH-SS求解精度較高,誤差均值為1%,TH-SE求解精度稍差,誤差均值為10%。通過對CGB/5000中IMP、TH-OS、TH-SS與TH-SE各工況撞擊力進行FFT頻譜分析,結果如圖11所示。

圖11 撞擊力頻譜成分

根據圖11,可得出:CGB/5000 TH-SS、TH-SE撞擊力在1.0~2.2 Hz以內的荷載成分幅值高于IMP,而在2.2~2.5 Hz范圍的荷載成分幅值則低于IMP。由表11可知,CGB/5000 TH-SS的誤差均值為1%,由此可知在1.0~2.5 Hz范圍內的振型反應誤差經過相互抵消,結構總反應求解精度較好;TH-SE撞擊力在1.0~2.2 Hz以內的荷載成分幅值低于TH-SS,在2.2~2.5 Hz范圍的荷載成分幅值與TH-SS相當,最終TH-SE的求解精度低于TH-SS。由此可知,修正半波正弦撞擊力模型的求解精度與橋梁結構的動力特性有關。

為從一般性角度討論修正半波正弦撞擊力模型的適用性,令qj(t)=γjδj(t),由式(17)可得

(18)

(19)

采用45條TH-OS撞擊力時程與45條TH-SE撞擊力時程,通過求解式(18)與式(19),得到各工況反應比r的譜值,并將各船舶9個撞擊速度得到的r譜值進行平均,如圖12所示。

圖12 振型反應比譜值

根據圖12,可得出以下信息:

(1)當f<0.5 Hz時,r基本在1.0左右,表明TH-OS與TH-SE得到的振型反應基本相同,主要是對于低頻單自由度體系,最大響應主要受撞擊荷載沖量控制;

(2)當0.5 Hz≤f<4.0 Hz時,0.8

(3)當f>4.0 Hz時,采用TH-SE簡化撞擊力模型得到的單自由度體系響應誤差超過20%,隨著單自由度體系固有頻率增加,誤差逐漸增大。

對于本文所選橋梁案例來說,結構的主要振型基本在2.5 Hz以內,由圖12可知,在此范圍內結構的振型響應的精度可以保持在±20%以內,與表11算例分析結果相符。實際上要全面評價修正半波正弦撞擊力模型的求解精度,需選擇足夠數量具有代表性的橋梁進行船橋碰撞動力響應分析,然而這需要一個不斷積累的過程。

4 結 論

本文建立了5艘不同噸位輪船的有限元模型,采用正撞剛性墻模型,獲得了45條不同撞擊速度下的撞擊力時程樣本過程,通過本文的討論,得出的主要結論如下:

(1)建立了修正半波正弦簡化撞擊力模型,給出了模型的參數確定方法,并確定了參數取值。

(2)算例分析表明,本文建立的修正半波正弦簡化撞擊力模型(TH-SE)平均計算誤差在10%以內,具有良好的求解精度,因此具有較好的工程實用價值。

(3)本文荷載模型建立在剛性橋梁基礎之上,橋例分析表明,剛性橋梁假定造成的響應求解誤差在8%以內,因此采用該假定用于確定船舶撞擊荷載是可以接受的。

(4)隨著橋梁算例的不斷積累,本文所建立的簡化撞擊力模型的誤差能得到更具體的解釋,荷載模型和參數的合理性將得到進一步的檢驗,荷載模型的適用范圍也將更明確。這需要一個不斷積累的過程與更深入的研究探討。

猜你喜歡
船舶橋梁模型
一半模型
計算流體力學在船舶操縱運動仿真中的應用
《船舶》2022 年度征訂啟事
船舶(2021年4期)2021-09-07 17:32:22
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
船舶!請加速
手拉手 共搭愛的橋梁
句子也需要橋梁
高性能砼在橋梁中的應用
3D打印中的模型分割與打包
主站蜘蛛池模板: 成人福利免费在线观看| 亚洲不卡影院| 欧美精品二区| 亚洲va欧美va国产综合下载| 国产一区二区三区精品欧美日韩| 99久久国产精品无码| 久久综合结合久久狠狠狠97色 | 免费高清a毛片| 国产伦片中文免费观看| 成人伊人色一区二区三区| 亚洲一级色| 日本尹人综合香蕉在线观看| 中文字幕精品一区二区三区视频 | 精品国产免费观看一区| 91久久夜色精品国产网站| 午夜在线不卡| 国产尹人香蕉综合在线电影| 色婷婷综合激情视频免费看| 亚洲青涩在线| 国产免费久久精品44| 日韩精品亚洲精品第一页| 国产凹凸视频在线观看| 日韩欧美国产中文| 91美女视频在线观看| 欧美一级视频免费| 午夜国产精品视频| h视频在线观看网站| 免费看黄片一区二区三区| 九九热在线视频| 亚洲女同欧美在线| 国内毛片视频| 日韩a级片视频| 深夜福利视频一区二区| 亚洲精品第一在线观看视频| 热这里只有精品国产热门精品| 精品国产美女福到在线不卡f| 亚洲色中色| 无码中文字幕精品推荐| 久草视频精品| 精品国产网| 伊人激情综合| 亚洲全网成人资源在线观看| 伊人成人在线视频| 久久视精品| 黄网站欧美内射| 91精品人妻一区二区| 亚洲成av人无码综合在线观看| 天堂亚洲网| 91久久国产综合精品女同我| 热伊人99re久久精品最新地| 久久久久国产精品嫩草影院| 四虎影视8848永久精品| 国产小视频免费| 日韩毛片在线视频| 免费在线观看av| 91免费国产在线观看尤物| 91小视频在线观看| 秋霞国产在线| 欧美日韩高清在线| 呦女亚洲一区精品| 狠狠色狠狠色综合久久第一次| 久草视频中文| 在线观看免费人成视频色快速| 国产香蕉在线| 99在线视频免费| 亚洲视频无码| 精品久久香蕉国产线看观看gif| 日韩毛片基地| 91亚洲影院| 国产91透明丝袜美腿在线| 亚洲v日韩v欧美在线观看| 中文字幕无线码一区| 中国丰满人妻无码束缚啪啪| 亚洲v日韩v欧美在线观看| 久青草免费在线视频| 国产人成在线视频| 国产理论最新国产精品视频| 91无码人妻精品一区| 国产精品免费p区| 亚洲天堂免费观看| 亚洲第一视频免费在线| 国产免费久久精品99re丫丫一|