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

球鼻船艏斜撞船舶舷側(cè)結(jié)構(gòu)的變形機(jī)理研究

2022-06-17 03:13:26王澤平胡志強(qiáng)
振動與沖擊 2022年11期
關(guān)鍵詞:船舶變形結(jié)構(gòu)

王澤平, 胡志強(qiáng), 劉 昆, 陳 剛

(1.上海交通大學(xué) 海洋工程國家重點(diǎn)實(shí)驗(yàn)室,上海 200240;2.哈爾濱工程大學(xué) 船舶工程學(xué)院,哈爾濱 150001;3.英國紐卡斯?fàn)柎髮W(xué) 工程學(xué)院,英國紐卡斯?fàn)?NE1 7RU;4.江蘇科技大學(xué) 船舶與海洋工程學(xué)院,鎮(zhèn)江 212003;5.中國船舶及海洋工程設(shè)計研究院,上海 200011)

隨著海上航行船舶數(shù)量的不斷增加,船舶發(fā)生碰撞事故的概率也逐漸增加。盡管很多先進(jìn)的避撞設(shè)備已經(jīng)應(yīng)用到船舶航行中,但船舶發(fā)生碰撞事故的風(fēng)險仍然存在;船舶碰撞事故會導(dǎo)致環(huán)境污染、經(jīng)濟(jì)損失甚至人員傷亡。有統(tǒng)計數(shù)據(jù)表明,船舶發(fā)生斜撞事故的概率要高于船舶正撞事故[1]。然而目前對于船舶碰撞事故的研究主要集中在船舶正撞場景,對于船舶斜撞場景開展的研究相對較少,因此需要針對船舶斜撞場景開展相關(guān)結(jié)構(gòu)耐撞性研究。

船舶碰撞的研究方法可以分為四種,包括經(jīng)驗(yàn)公式法,試驗(yàn)法,有限元數(shù)值仿真法和解析計算法。Minorsky[2]作為經(jīng)驗(yàn)公式法的開創(chuàng)者,通過對已發(fā)生的船舶碰撞事故結(jié)構(gòu)損傷進(jìn)行統(tǒng)計分析,進(jìn)而提出內(nèi)能吸收與結(jié)構(gòu)損傷體積的經(jīng)驗(yàn)公式。試驗(yàn)法分析船舶碰撞問題的準(zhǔn)確性較好,不過研究成本較高。Pedersen等[3]通過試驗(yàn)法研究船舶結(jié)構(gòu)的耐撞性;Amdahl等[5]通過開展模型試驗(yàn)研究船舶碰撞與擱淺時的結(jié)構(gòu)變形機(jī)理;Liu等[6-7]通過落錘沖擊試驗(yàn)研究了應(yīng)變率的影響。隨著計算機(jī)技術(shù)的發(fā)展,有限元數(shù)值仿真法被廣泛應(yīng)用于船舶碰撞問題的研究中。王自力等[8]通過數(shù)值仿真法研究了船舶碰撞動力學(xué)問題。劉敬喜等[9]通過數(shù)值仿真法研究了雙殼油船舷側(cè)結(jié)構(gòu)的耐撞性。劉昆等[10]采用全耦合技術(shù)分析了船體結(jié)構(gòu)的碰撞性能。此外,解析計算法是一種成本低且可以提供理想精度的方法,其通過提出結(jié)構(gòu)的理論變形模型,運(yùn)用塑性力學(xué)等理論,推導(dǎo)得到解析計算公式。Alexander[11]創(chuàng)造性地將解析計算法應(yīng)用到薄板結(jié)構(gòu)在軸向加載情況下的分析計算中。高振國等[12]通過解析計算法分析了強(qiáng)肋框結(jié)構(gòu)在面內(nèi)加載時的變形機(jī)理。Hong等[13]研究了船舶碰撞與擱淺過程中桁材結(jié)構(gòu)的變形機(jī)理,并提出了解析公式。

近年來,隨著船舶噸位的增加,船舶球鼻艏的安裝也更加普遍,例如補(bǔ)給船,穿梭油輪等。船舶球鼻艏結(jié)構(gòu)通常比較堅硬,發(fā)生船舶碰撞事故時更容易引起被撞船舷側(cè)結(jié)構(gòu)的損傷,進(jìn)而導(dǎo)致艙室進(jìn)水、漏油等事故的發(fā)生。為了減少因?yàn)榍虮囚甲矒魧?dǎo)致的損失,研究舷側(cè)結(jié)構(gòu)在球鼻艏撞擊場景下的結(jié)構(gòu)響應(yīng)具有重要意義。然而大多數(shù)已有的研究都集中在球鼻艏正撞舷側(cè)結(jié)構(gòu)的場景,比如,Lee等[14-15]研究了舷側(cè)外板在球鼻艏正撞場景下的結(jié)構(gòu)響應(yīng)。對于球鼻艏斜撞場景下舷側(cè)結(jié)構(gòu)耐撞性的研究較少。在斜撞場景下,舷側(cè)結(jié)構(gòu)的變形情況與撞擊角度相關(guān)性很大,舷側(cè)結(jié)構(gòu)的主要構(gòu)件,例如舷側(cè)外板、橫向肋板等構(gòu)件在斜撞場景下的變形模式與正撞場景是不同的,因此需要針對舷側(cè)結(jié)構(gòu)在球鼻艏斜撞場景下的變形機(jī)理開展研究。

1 解析計算基本理論

解析計算方法可以快速評估被撞船的耐撞性,通過應(yīng)用塑性力學(xué)理論,求解舷側(cè)結(jié)構(gòu)主要構(gòu)件在球鼻艏斜撞場景下的抵抗力。根據(jù)上限定理,外力對結(jié)構(gòu)的做功功率等于結(jié)構(gòu)內(nèi)部構(gòu)件的能量耗散率[16]

(1)

(2)

(3)

(4)

N0=σ0t

(5)

式中:σ0為流動應(yīng)力,t為舷側(cè)外板的厚度。

2 舷側(cè)外板的變形機(jī)理研究

2.1 舷側(cè)外板幾何變形模型

橢球形撞頭斜撞一塊邊長分別為2a和2b的矩形板,如圖1所示。該橢球形撞頭可以用如下方程的形式表示

(6)

圖1 橢球形撞頭斜撞矩形板示意圖Fig.1 A rectangular plate struck by an ellipsoid indenter

假設(shè)截面為橢圓的撞頭撞擊矩形板,如圖2所示。該橢圓截面是橢球形撞頭在xoy平面內(nèi)的截面,那么橢圓形截面可以用如下方程的形式表示

(7)

圖2 截面為橢球形的撞頭斜撞矩形板示意圖Fig.2 A rectangular plate struck by an indenter (b1+b2=2b)

矩形板的變形模式如圖3所示。在矩形板受壓變形的過程中,左側(cè)部分矩形板的應(yīng)變可以用如下方程表示

圖3 矩形板的變形模式示意圖Fig.3 Deformation model of the rectangular plate

因此,變形矩形板左側(cè)部分的應(yīng)變率是

(8)

式中:l1是矩形板左側(cè)部分未發(fā)生彎曲部位的瞬時長度;α1是矩形板左側(cè)部分的瞬時旋轉(zhuǎn)角度。

(9)

同樣,右側(cè)部分矩形板的應(yīng)變可以表示為

(10)

變形矩形板右側(cè)部分的應(yīng)變率是

(11)

將式(9)和(11)代入式(3)中,可以得到膜拉伸能量耗散率為

(12)

N0_sideplate=σ0tp

(13)

式中:N0_sideplate是塑性膜拉伸力;tp是矩形板的厚度。

橢球形撞頭的撞擊深度可以用下式表示

(14)

其中,b1,b2的長度可以分別表示為

b1=x1tanα1+l1cosα1

(15)

b2=x2tanα2+l2cosα2

(16)

式中,(x1,y1)和(x2,y2)是矩形板未彎曲部分與彎曲部分的交點(diǎn)。x1,x2分別為變形矩形板未彎曲部分與彎曲部分的交點(diǎn)在x軸方向上的坐標(biāo)值,如圖3所示。x1和x2可以表示為

(17)

(18)

由此可以得到橢球形撞頭的撞擊深度

(19)

橢球形撞頭的撞擊深度對瞬時旋轉(zhuǎn)角度求導(dǎo),可以得到撞擊速率為

(20)

由于矩形板的彎曲能量耗散占總能量耗散的比例很小,所以通常忽略彎曲能量耗散。因此,矩形板在橢球形撞頭撞擊時的瞬時阻力可以表示為

(21)

同樣,假設(shè)截面為橢圓的撞頭撞擊矩形板,該橢圓截面是橢球形撞頭在xoz平面內(nèi)的截面,那么橢圓形截面可以用如下方程的形式表示

(22)

矩形板受到撞擊時的瞬時阻力計算方法與上文提出的方法相同,可以表示為如下形式

(23)

Haris等的研究發(fā)現(xiàn)橢球形撞頭撞擊矩形板的瞬時阻力與上文提出的矩形板受到兩個方向上的橢圓形截面撞頭撞擊時的瞬時阻力關(guān)系可以表示為

(24)

在斜撞過程中,橢球形撞頭和矩形板之間的摩擦力可以表示為

FS=μ·Fp·sinβ

(25)

由于摩擦產(chǎn)生的能量耗散可以表示為

(26)

2.2 舷側(cè)外板破裂的預(yù)報

在碰撞過程中,如果外板發(fā)生破裂,那么外板的抵抗力會出現(xiàn)明顯下降,因此能夠較為準(zhǔn)確地預(yù)測出外板的破裂十分有必要。通常臨界斷裂應(yīng)變被當(dāng)做板受到面外載荷作用時發(fā)生膜拉伸的最大應(yīng)變,即認(rèn)為當(dāng)板發(fā)生膜拉伸的應(yīng)變值達(dá)到臨界應(yīng)變值時,板出現(xiàn)破裂。這一方法被廣泛使用,因此本文也采用這種方法預(yù)測外板的破裂。低碳鋼的塑性拉伸應(yīng)變值通常在0.2~0.35,在船舶的初步設(shè)計階段,設(shè)計師可以根據(jù)具體情況確定這一數(shù)值的大小。將式(8)代入式(19)中,可以得到舷側(cè)外板最大膜拉伸應(yīng)變和球鼻艏撞深之間的關(guān)系

(27)

式中,εm為外板的臨界斷裂應(yīng)變值。根據(jù)式(27)可以得到舷側(cè)外板發(fā)生破裂時球鼻艏的撞擊深度。

2.3 舷側(cè)外板破裂后的解析計算公式

隨著碰撞過程的進(jìn)行,球鼻艏撞深增加,舷側(cè)外板將會發(fā)生破裂。Wang等分析了板的變形過程,提出了預(yù)測板破裂后的阻力計算公式

(28)

式中:n為板破裂后的裂紋條數(shù),n>2;l為裂縫的長度;φ為球鼻艏頂角的一半;μ為摩擦因數(shù)。

3 桁材變形模型

為了研究桁材結(jié)構(gòu)在不同斜撞角度下的變形機(jī)理,首先對桁材的斜撞過程進(jìn)行了數(shù)值模擬。選取三種典型的斜撞角度,數(shù)值模擬了球鼻艏結(jié)構(gòu)在30°、45°和60°三種場景下斜撞桁材的過程。桁材在30°、45°和60°場景下的有限元模型、變形模型以及橫截面的變形過程如圖4~6所示。

(a) 有限元模型

(b) 變形模型

(c) 橫截面變形過程圖4 桁材在數(shù)值模擬中的變形(β=30°)Fig.4 Deformation of the web girder in numerical simulation(β=30°)

(a) 有限元模型

(b) 變形模型

(c) 橫截面變形過程圖5 桁材在數(shù)值模擬中的變形(β=45°)Fig.5 Deformation of the web girder in numerical simulation(β=45°)

(a) 有限元模型

(b) 變形模型

(c) 橫截面變形過程圖6 桁材在數(shù)值模擬中的變形 (β=60°)Fig.6 Deformation of the web girder in numerical simulation

從圖4~6中可以看出,在不同的斜撞角度下,桁材腹板的折疊高度比是相似的,因此根據(jù)數(shù)值模擬的結(jié)果提出了桁材的理論變形模型,如圖7和圖8所示。其中,b1為受壓桁材結(jié)構(gòu)的一側(cè)長度,另一側(cè)長度為b2=b-b1。

圖7 桁材撞擊后的理論變形模型Fig.7 The deformation model of the web girder after collision

圖8 橫截面受壓折疊變形過程(黑點(diǎn)代表塑性鉸所在位置)Fig.8 Folding deformation process of the cross section (the black dots represent the location of the plastic hinges)

根據(jù)所提出的桁材理論變形模型,通過推導(dǎo)可以得到適用于不同斜撞角度的桁材變形阻力計算公式。限于篇幅,詳細(xì)推導(dǎo)過程請參閱Wang等[17]9-10。推導(dǎo)后得到的桁材結(jié)構(gòu)在形成第1個褶皺時的瞬時阻力公式為

(29)

在形成第1個褶皺時的平均阻力公式為

(30)

根據(jù)上限定理,當(dāng)平均阻力最小時,能量耗散也最小,即:

(31)

將式(30)代入式(31)中,可以得到:

(32)

同樣,當(dāng)桁材結(jié)構(gòu)形成第n個褶皺時(n≥2),瞬時阻力公式為

(33)

式中:δ′為球鼻艏在斜撞方向上的位移距離;t為腹板的厚度。

桁材上帶板的阻力公式可以采用梁理論進(jìn)行計算

(34)

式中:a是帶板寬度的一半;b1+b2是帶板的長度;tf是帶板的厚度。

4 橫向肋板變形模型

球鼻艏斜撞橫向肋板時,假設(shè)橫向肋板受到長度

為d的線荷載作用,切向碰撞力使橫向肋板上邊緣產(chǎn)生水平位移,垂向碰撞力使橫向肋板產(chǎn)生二維折疊變形,所提出的橫向肋板理論變形模型,分別如圖9和圖10所示。

圖9 橫向肋板碰撞前后的橫截面Fig.9 The cross section of transverse frame before and after collision

圖10 橫向肋板受壓變形模式Fig.10 Deformation mode of the crushing transverse frame

根據(jù)所提出的橫向肋板理論變形模型,通過推導(dǎo)可以得到適用于不同斜撞角度的切向力和垂向力的計算公式。限于篇幅,詳細(xì)推導(dǎo)過程請參閱Wang等[17]10-12。推導(dǎo)后得到的切向力和垂向力的計算公式如下

(35)

(36)

式中:δ′為球鼻艏在斜撞方向上的位移距離;t為肋板的厚度。

5 數(shù)值仿真

本文通過開展數(shù)值仿真分析來驗(yàn)證所提出解析公式的合理性與準(zhǔn)確性。

5.1 有限元模型

在數(shù)值仿真中,選取一艘雙殼油輪作為被撞船,一艘60 000 t船舶的橢球形球鼻艏作為撞擊船艏。撞擊船艏和被撞船舷側(cè)結(jié)構(gòu)的有限元模型如圖11所示,表1列出了被撞船的主尺度。

圖11 油輪舷側(cè)結(jié)構(gòu)和撞擊船艏的有限元模型Fig.11 Finite element model of the side structure of the tanker and the bow

表1 被撞船主尺度表Tab.1 Dimentions of the struck ship

使用有限元軟件LS_DYNA進(jìn)行數(shù)值仿真,使用Johnson-Cook GTN材料模型[18]。撞擊船艏的材料屬性定義為剛體,被撞船舷側(cè)結(jié)構(gòu)的材料參數(shù)如表2和表3所示。撞擊速度定義為3 m/s,摩擦因數(shù)定義為0.3[19]。有限元模型中采用的是四節(jié)點(diǎn)四邊形的Belytschko-Tsay單元。為了在保證數(shù)值模擬精度的同時有合理的計算時間,有限元模型的網(wǎng)格尺寸設(shè)置為200 mm。被撞船舷側(cè)兩端采用固定邊界條件;本文研究的是船舶斜撞過程中相對危險的情況,即撞擊船在碰撞過程中不發(fā)生偏轉(zhuǎn),被撞船舷側(cè)結(jié)構(gòu)在較大撞深情況下的損傷變形機(jī)理。在數(shù)值模擬中,被撞船和撞擊船之間的接觸設(shè)置為面面接觸。

表2 GTN模型參數(shù)Tab.2 GTN parameters

表3 Johnson-Cook模型參數(shù)Tab.3 Johnson-Cook parameters

5.2 碰撞場景定義

船舶碰撞時斜撞角度的定義如圖12所示,圖13為三種典型的碰撞位置,兩種碰撞場景的定義分別如表4和圖14所示。工況1中的碰撞位置為位置1,工況2中的碰撞位置為位置2,隨著碰撞的進(jìn)行,撞擊船艏會經(jīng)過位置3。

圖12 斜撞角度β的定義Fig.12 Definition of the collision angle β

圖13 三種典型碰撞位置Fig.13 Three typical collision positions

表4 不同碰撞場景中的碰撞位置Tab.4 Definition of collision positions in different collision scenarios

(a) 工況1

(b) 工況2圖14 兩種斜撞場景的定義Fig.14 Definition of two oblique collision cases

5.3 碰撞過程中參與構(gòu)件的確定

在5.2節(jié)中確定了碰撞場景,可以結(jié)合撞擊球鼻艏結(jié)構(gòu)的輪廓圖和被撞船舷側(cè)的結(jié)構(gòu)圖確定參與斜撞過程中的船體構(gòu)件。圖15(a)是工況1,撞擊過程通過將撞擊球鼻艏結(jié)構(gòu)的輪廓圖和被撞船舷側(cè)的結(jié)構(gòu)圖結(jié)合起來進(jìn)行模擬。在碰撞開始時,與球鼻艏接觸的構(gòu)件有舷側(cè)外板和桁材,所以舷側(cè)結(jié)構(gòu)的抵抗力由這兩個構(gòu)件產(chǎn)生,舷側(cè)外板的變形區(qū)域?yàn)镃DHI;隨著碰撞的進(jìn)行,橫向肋板也參與到舷側(cè)結(jié)構(gòu)抵抗力的計算中,舷側(cè)外板的變形區(qū)域?yàn)镃EGI,如圖15(b)所示。圖15(c)和(d)是工況2,在碰撞開始時,與球鼻艏接觸的構(gòu)件有舷側(cè)外板和橫向肋板,舷側(cè)結(jié)構(gòu)的抵抗力由這兩個構(gòu)件產(chǎn)生,舷側(cè)外板的變形區(qū)域?yàn)锳HJK,如圖15(c)所示;隨著碰撞的進(jìn)行,舷側(cè)外板的變形區(qū)域?yàn)镕GJK,如圖15(d)所示。

(a) 工況1

(b) 工況1

(c) 工況2

(d) 工況2圖15 碰撞過程中變形區(qū)域及參與構(gòu)件的確定Fig.15 Determination of deformation area and participating components during collision

6 數(shù)值驗(yàn)證與討論

使用有限元軟件LS_DYNA對兩種斜撞場景進(jìn)行數(shù)值仿真,將得到碰撞過程中球鼻艏撞擊深度與碰撞力和能量耗散之間的關(guān)系曲線;解析方法中的能量耗散是通過對膜拉伸能量耗散率積分得到,將數(shù)值仿真結(jié)果與解析計算結(jié)果進(jìn)行對比,如圖16~圖19所示。

圖16 碰撞場景1中碰撞力-撞深曲線對比結(jié)果Fig.16 The comparison results of resistance-depth curves in case 1

圖17 碰撞場景1中能量-撞深曲線對比結(jié)果Fig.17 The comparison results of energy-depth curves in case 1

圖18 碰撞場景2中碰撞力-撞深曲線對比結(jié)果Fig.18 The comparison results of resistance-depth curves in case 2

圖19 碰撞場景2中能量-撞深曲線對比結(jié)果Fig.19 The comparison results of energy-depth curves in case 2

從圖16~圖19的對比結(jié)果中可以看出,數(shù)值仿真方法和解析方法得出的結(jié)果吻合良好。隨著碰撞過程的進(jìn)行,碰撞力和能量耗散曲線呈上升趨勢。當(dāng)碰撞場景1和2中球鼻艏的撞擊深度達(dá)到3.2 m時,舷側(cè)外板發(fā)生破裂;舷側(cè)外板發(fā)生破裂后抵抗力會下降,能量耗散曲線的增長速度也出現(xiàn)下降,解析方法和數(shù)值仿真方法都能捕捉到這一特點(diǎn)。在舷側(cè)外板破裂初期,解析方法和數(shù)值仿真方法得到的碰撞力-撞深曲線存在一些偏差,這是因?yàn)樵谙蟼?cè)外板破裂初期,膜拉伸能量耗散仍然是舷側(cè)外板最主要的能量耗散形式,舷側(cè)外板的變形模式不會在很短的時間內(nèi)轉(zhuǎn)化為撕裂變形模式,這是解析方法和數(shù)值仿真方法得到的結(jié)果之間存在偏差的主要原因??傮w來說,解析方法和數(shù)值仿真方法得到的結(jié)果對比良好,這也驗(yàn)證了本文所提出的解析公式的合理性。

7 結(jié) 論

本文研究了船舶球鼻艏結(jié)構(gòu)斜撞船舶舷側(cè)時,舷側(cè)結(jié)構(gòu)主要構(gòu)件的變形機(jī)理,提出了舷側(cè)外板、桁材和橫向肋板三種構(gòu)件的理論變形模型,運(yùn)用塑性力學(xué)理論推導(dǎo)得到了三種構(gòu)件在不同斜撞角度下變形過程中的解析計算公式。定義了兩種典型的斜撞場景,通過有限元軟件LS_DYNA進(jìn)行了數(shù)值仿真,數(shù)值仿真結(jié)果驗(yàn)證了解析計算結(jié)果的準(zhǔn)確性,解析計算公式可以很好地捕捉到舷側(cè)外板發(fā)生破裂后碰撞力和能量耗散增長速度下降這一特點(diǎn)。

球鼻艏也是船舶常見的設(shè)計形式,船舶斜撞事故發(fā)生的概率要高于船舶正撞事故,而目前針對船舶碰撞的研究主要集中在船舶正撞場景,對于船舶斜撞場景的研究是不充分的,因此,研究舷側(cè)結(jié)構(gòu)在船舶球鼻艏結(jié)構(gòu)斜撞場景下的變形機(jī)理是很有意義的。本文對舷側(cè)結(jié)構(gòu)在球鼻艏斜撞場景下的變形機(jī)理進(jìn)行了研究,所得到的解析公式對于船舶舷側(cè)結(jié)構(gòu)耐撞性的評估具有較好的工程應(yīng)用價值。

猜你喜歡
船舶變形結(jié)構(gòu)
計算流體力學(xué)在船舶操縱運(yùn)動仿真中的應(yīng)用
《船舶》2022 年度征訂啟事
船舶(2021年4期)2021-09-07 17:32:22
《形而上學(xué)》△卷的結(jié)構(gòu)和位置
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
船舶!請加速
論結(jié)構(gòu)
中華詩詞(2019年7期)2019-11-25 01:43:04
“我”的變形計
例談拼圖與整式變形
會變形的餅
論《日出》的結(jié)構(gòu)
主站蜘蛛池模板: 亚洲 欧美 日韩综合一区| 香蕉色综合| 456亚洲人成高清在线| 国产精品欧美在线观看| 国产女人在线观看| 亚洲高清资源| 不卡网亚洲无码| 青青操国产视频| 9966国产精品视频| 欧美国产精品拍自| 日韩毛片基地| 亚洲成a∧人片在线观看无码| 熟女成人国产精品视频| 综合亚洲网| 99精品热视频这里只有精品7| 亚洲精品无码久久毛片波多野吉| 91精品国产丝袜| 亚洲欧洲一区二区三区| 日本人妻一区二区三区不卡影院 | 波多野结衣AV无码久久一区| 久久精品国产亚洲麻豆| 欧美一区中文字幕| 无码aⅴ精品一区二区三区| 午夜少妇精品视频小电影| 国产成人一级| 在线看国产精品| 精品亚洲欧美中文字幕在线看| 中文字幕 日韩 欧美| 久久精品这里只有国产中文精品 | 亚洲国语自产一区第二页| 亚洲国产欧美国产综合久久| 亚洲人在线| 嫩草国产在线| 亚洲一级色| 亚洲AⅤ综合在线欧美一区| 性欧美在线| 中文无码日韩精品| 72种姿势欧美久久久久大黄蕉| 国产97视频在线| 久夜色精品国产噜噜| 激情无码字幕综合| 国产成人亚洲无吗淙合青草| 精品1区2区3区| 午夜福利在线观看成人| 中国毛片网| 亚洲天堂网2014| 成人字幕网视频在线观看| 国产毛片高清一级国语 | 国产亚洲欧美在线中文bt天堂| 国产18在线| 狠狠综合久久| 婷婷99视频精品全部在线观看 | 午夜福利无码一区二区| 怡红院美国分院一区二区| 国产无码精品在线播放| 日韩欧美国产中文| 理论片一区| 被公侵犯人妻少妇一区二区三区| 综合色天天| 福利在线不卡| 超碰91免费人妻| 婷婷激情亚洲| 亚洲日本中文综合在线| 精品国产电影久久九九| 人妻21p大胆| 91视频日本| 国产日韩欧美黄色片免费观看| 亚洲国产午夜精华无码福利| 国产成人凹凸视频在线| 成人午夜网址| 99视频精品在线观看| 青草视频久久| 伊人五月丁香综合AⅤ| 国产精鲁鲁网在线视频| 国内精品九九久久久精品 | 在线日本国产成人免费的| av在线无码浏览| 国产欧美日韩精品第二区| 久久精品国产亚洲AV忘忧草18| 国产丝袜无码精品| 青草91视频免费观看| 凹凸国产分类在线观看|