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

超空泡運動體圓柱薄殼動力屈曲及可靠性分析

2014-08-11 14:50:07王杰方安偉光宋向華
振動與沖擊 2014年8期
關鍵詞:區域

王杰方, 安偉光, 宋向華

(哈爾濱工程大學 航天工程系, 哈爾濱 150001)

超空泡運動體圓柱薄殼動力屈曲及可靠性分析

王杰方, 安偉光, 宋向華

(哈爾濱工程大學 航天工程系, 哈爾濱 150001)

首先將超空泡運動體模擬成受動態軸向載荷作用的圓柱薄殼,推導了結構的動力穩定性微分方程和動力不穩定區域,然后考慮動態軸向載荷的隨機性,采用有限步長迭代法將給出的動力屈曲失穩的多個安全余量方程線性化,并利用逐步搜索法找出有效的安全余量方程,最后結合逐步等效平面法計算了艙段動力屈曲的可靠性指標。通過算例分別分析了載荷頻率、速度和載荷比例系數這三個隨機參數的變化對動力屈曲可靠性的影響。計算結果為如何選擇載荷頻率、速度和載荷比例系數的安全范圍提供了理論依據。

動力屈曲;可靠性;有限步長迭代法;逐步搜索法;逐步等效平面法

超空泡運動體是在低壓汽化或者人工通氣形成的超空泡中高速航行的運動體,其頭部空化器和水相互接觸所產生的阻力隨速度的二次方增長,結構所承受的軸向載荷非常大,容易發生屈曲,另外,空泡形狀和尺寸的不穩定性會導致載荷隨時間變化,對于承受動態軸向載荷作用的圓柱薄殼結構,一般來說,結構只有軸向振動,但是,當軸向載荷的擾動頻率與結構橫向的固有頻率之間的比值存在某種關系時,圓柱薄殼結構的橫向振動的振幅迅速增大,導致結構將喪失動力穩定性。因此,研究超空泡運動體的動力屈曲問題是十分必要的。

目前,國內外對超空泡運動體動力屈曲問題的研究主要有:Ruzzene[1]將超空泡運動體模擬成彈性的軸對稱殼體,采用Bolotin方法分析了各種軸向屈曲載荷和不同的殼體結構形式下的動態屈曲穩定性;Edward 等[2]用有限元法和多學科優化方法改進了結構的動態進了結構的動態穩定性;施連會等[3]用Bolotin方法對水下航行體的動力穩定性問題進行了數值計算,并分析了載荷參數和航行體變化參數對主動力不穩定區域的影響規律;宋向華等[4]對超空泡射彈的截錐形結構的動力穩定性進行了計算。運動體的靜力屈曲可靠性分析的文獻主要有:顧永維等[5]將超空泡運動體簡化為變截面梁,研究了其抗彎穩定可靠性的問題;周凌等[6]將運動體的環向加肋艙段簡化為變厚度的圓柱薄殼,對其進行了屈曲可靠性分析,同時,周凌等[7]對超空泡運動體強度和穩定性兩失效模式組成的串聯系統進行了非概率可靠性分析。

現有資料中,分析超空泡運動體的動力屈曲和靜力屈曲可靠性的文獻較多,但是分析動力屈曲可靠性的文獻極少,而實際情況下,與載荷相關的幾何參數、物理參數、流場參數等都存在隨機性,結構的動力不穩定區也存在隨機性,因此,在確定性動力屈曲分析的基礎上,研究超空泡運動體動力屈曲的可靠性是十分必要的。

文中將有限步長迭代法和逐步等效平面法相結合,將功能函數線性化,并計算結構動力屈曲的可靠度指標。有限步長迭代法對于非線性程度較高的極限狀態方程收斂速度快、迭代精度高,而等效平面法能得到體系對應的線性功能函數,有利于進一步結合其他子系統計算總系統的可靠度指標。另外,由于圓柱薄殼的動力不穩定區的數目眾多,即動力屈曲的安全余量方程眾多,而實際影響結構安全的只是那些離載荷頻率較近的不穩定區域,稱為有效不穩定區,因此,文中還將給出搜索有效安全余量方程的逐步搜索法。

1 超空泡運動體受力分析及動力穩定性分析

1.1 超空泡運動體圓柱薄殼艙段受力分析

對于超空泡運動體水平向前運動時的動力屈曲分析,其受力可以簡化為軸向的均布載荷,即頭部阻力和尾部推力,二者大小相等,作用在頭部空化器的阻力[8]為

(1)

其中,ρw為流體密度,Ac是空化器的橫截面積,Cx是空化器的阻力系數,V是航行體的運動速度。零攻角時,圓盤空化器的阻力系數為

Cx=Cx0(1+σ)

(2)

Cx0=0.5+1.81(φ/360-0.25)-2(φ/360-0.25)2

(3)

(4)

其中,φ為空化器錐角,圓盤空化器的錐角180°,σ為空化數,pc為空泡內壓力,p∞為環境壓力,其表達式為

p∞=ρwgH+p

(5)

其中,H為航行深度,p為標準大氣壓。

運動體在高速航行過程中,空泡形狀和尺寸的不穩定性會導致軸壓幅值p0隨時間變化,本文將這一動態的軸向載荷簡化為[1]

p(t)=p0+pδ(t)=p0+(1+δcosθt)

(6)

(7)

其中,p0為不隨時間變化的部分,pδ(t)為隨時間變化的部分,δ為擾動載荷的比例系數,是擾動載荷的幅值與載荷p0的比值,θ為擾動載荷的角頻率,dn為圓盤空化器的直徑。

1.2 動力穩定性微分方程

根據符拉索夫的彈性殼體理論,適用于細長圓柱薄殼的殼體穩定平衡方程應滿足:

(1) 平衡方程中,在圓周切線方向考慮橫向剪切力的影響;

(2) 幾何方程中,在彎曲變形的幾何關系中計入切向位移的影響;

(3) 物理方程中,中面力中計入了彎曲變形的影響,而中面外的力中則計入了中面應變的影響。

聯立平衡方程、幾何方程和物理方程,并引入函數Φ=Φ(α,β,t),得到用Φ表示的圓柱薄殼艙段的穩定平衡方程為[9]

(2+1)222Φ-

(8)

令Φ(α,β,t)=f(t)sinnαcoskβ(n=iπR/L)(薄殼兩端沒有徑向位移),i,k取整數值,分別為軸向和周向的半波數,代入式(8)得到圓柱薄殼動力穩定性微分方程為

(9)

其中,

g(n,k)=

從式(9)的量剛分析得,無彎矩狀態下的軸向內力N1的量綱應為N/m,即p(t)的量綱為N,文獻[3]中p(t)的量綱是N/m2,是錯誤的。

整理上式得

(10)

1.3 動力不穩定區

上述方程即mathieu方程,它有一個重要的性質:當它的系數(Ωn,k,γ,θ)存在某些關系時,方程式具有無限增長的解,這些無限增長的解在參數平面上布滿了許多區域,即動力不穩定區,根據Bolotin方法,圓柱薄殼的前三階動力不穩定區域邊界為[10]

(11)

(12)

(13)

其中,θ11和θ12,θ21和θ22,θ31和θ32分別為第一、第二和第三階動力不穩定區域的上邊界和下邊界。

文獻[3]給出了以速度V為變量的動力不穩定區域,是二維的動力不穩定區域,由式(7)可知,其他參數取定值時,p0取決于運動體的航行速度V,而p0決定了mathieu方程中的參數Ωn,k,γ,載荷比例系數δ決定了mathieu方程中的參數γ,因此,下文中給出的以V和δ兩者同時作為變量的三維動力不穩定區域能更全面的反映不穩定區的信息。

2 圓柱薄殼動力屈曲可靠性分析

2.1 動力屈曲安全余量方程

當超空泡體動態軸向載荷的頻率落入動力不穩定區間時,運動體將發生參數共振。圓柱薄殼動力屈曲可靠性分析即計算動態軸向載荷的頻率不落入結構動力不穩定區域的概率。從式(9)和式(11)~(13)可知,每組(i,k)分別對應著三個動力不穩定區域,圓柱薄殼艙段存在多個動力不穩定區域,因此,圓柱薄殼動力屈曲可靠性分析相當于多個安全余量方程的可靠性分析。每組(i,k)對應的第一、二、三階動力不穩定區域的功能函數形式一致,即

(14)

(15)

(16)

上述三式中,gi(X)<0,則動態軸向載荷的角頻率落入了動力不穩定區域;結構發生屈曲失效;gi(X)>0,則結構穩定;gi(X)=0,則結構有周期解,處于失效和可靠的臨界狀態。

2.2 動力屈曲可靠性分析方法

無論載荷頻率落入哪一個動力不穩定區域,結構都會失去穩定性,因此,動力屈曲的可靠性分析等同于多個失效函數串聯的可靠性分析。對于串聯系統而言,對各失效函數失效概率貢獻最大的驗算點,對結構體系失效概率的貢獻也最大。因此,近似計算串聯系統的失效概率時,可利用有限步長迭代法將各失效模式的功能函數在各自的驗算點處線性化,并采用等效平面法計算結構的動力屈曲指標。相比一般的迭代法,有限步長迭代法[11]對于非線性程度較高的極限狀態方程收斂速度快、迭代精度高,而等效平面法[11]是將那些經過有限步長迭代法線性化后的功能函數逐步等效,最終得到一個線性功能函數,這個線性功能函數中的常數項即為結構的動力屈曲可靠性指標,且這個線性功能函數有利于計算系統的可靠度指標。

對于第i個失效函數,有限步長迭代法在相互獨立的正態隨機變量空間中的迭代公式為

(17)

(18)

(j=1,2,…,n)

(19)

那么,經有限步長迭代法線性化的功能函數為

其中,

對于串聯系統,其失效概率采用逐步等效平面方法計算。逐步等效平面法是通過兩兩極限狀態面的依次等效,逐步減少極限狀態面的數目,最后剩下一個極限狀態面時,即為系統的功能函數。由式(20)中Ai恒大于0,那么,考慮兩個失效函數的情形,其線性化功能函數的等效函數為

(21)

其中,(B11,B12,…,B1n)和(B21,B22,…,B2n)均為單位向量。

兩兩等效平面對應的可靠指標公式為

βe=βp=-Φ-1[Φ(-β1)+

Φ(-β2)-Φ2(-β1,-β2;ρ12)]

(22)

其中,Φ2(-β1,-β2;ρ12)計算方法[9]如下:

(23)

(24)

B1=-A1(-β1+A1)

(25)

(26)

等效后各隨機變量對可靠度指標的敏感系數為

(27)

將敏感系數歸一化

(28)

兩個失效模式的等效線性功能函數為

(29)

2.3 逐步搜索法

有效不穩定區域是指對動力屈曲可靠性指標有主要影響的那些不穩定區域,從2.1節可知,圓柱薄殼艙段存在多個動力不穩定區域,精確計算超空泡運動體圓柱薄殼艙段的動力屈曲可靠度的過程應為:對每組(i,k)取值所對應的前三階動力不穩定區域都給出其安全余量方程,再線性化這些安全余量方程,最后利用逐步等效平面法來計算動力屈曲可靠度。顯然,對于那些遠離載荷頻率θ的動力不穩定區域,其可靠性指標很大,當某一不穩定區的可靠性指標β≥βc(βc為有效可靠性指標)時,可以認為這一不穩定區對系統的動力屈曲可靠性指標的影響極小,分析它對系統可靠度的影響是完全不必要的。

本文給出如圖1所示的逐步搜索法的計算框圖,來篩選有效不穩定區域。圖1中,取軸向最大半波數為a,周向最大半波數為b,搜索出a×b組振型對應的3×a×b個不穩定區中,影響動力屈曲可靠性指標的有效不穩定區域。圖1中,計算B1j,β1,B2j,β2,B3j,β3,采用式(17)~(20)的有效步長迭代法。搜索出有效的不穩定區域后,采用式(21)~(29)計算圓柱薄殼的動力屈曲可靠性指標即可。

圖1 逐步搜索法Fig.1 Step by step searching method

3 超空泡運動體圓柱薄殼艙段數值算例

超空泡運動體航行深度H=10 m,流場密度ρw=1 000 kg/m3,對于自然超空泡,空泡內的飽和蒸汽壓pc=2 350 Pa(20℃),標準大氣壓p=101 325 Pa。圓柱薄殼艙段的幾何參數:R=0.2 m,L=4 m,h=3 mm,dn=0.2 m。材料物理參數:E=209 GPa,ρ=7 850 kg/m3,μ=0.3。

3.1 動力屈曲計算結果及分析

由1.2節可知,i=0時,p*→∞,γ→0,即載荷比例系數δ→0,相當于受靜力載荷作用,不必進行動力屈曲及可靠性分析;計算表明,k=0時,f=θ/2π的最低值都在以上2 500 Hz,因此也不必進行動力屈曲及可靠性分析。

圖2(a)給出了i=1,k=1時,依賴于V、δ和f=θ/2π這三個參數的前三階動力不穩定區域;圖2(b)給出了i=1,k=1時,前三階不穩定區域寬度變化圖。

圖2(a) i=1,k=1時動力不穩定區域Fig.2(a) Dynamic unstable regions(i=1,k=1)

圖2(b) i=1,k=1時動力不穩定區域寬度Fig.2(b) The width of dynamic unstable regions(i=1,k=1)

結論1(由圖2得到):

速度和載荷比例系數增大,不穩定的載荷頻率值減小,且第一階不穩定載荷頻率減小的速度明顯大于后兩階;第一階動力不穩定區域寬度明顯大于后兩階的寬度。因此,對于承受動態軸向載荷且高速運動的超空泡運動體而言,必須充分考慮第一階不穩定區域對于結構的動力穩定性的影響,避免發生參數共振。

圖3給出i=2,k=1時,依賴V、δ和f=θ/2π這三個參數的前三階動力不穩定區域;圖4給出i=3,k=1時的前三階動力不穩定區域。

圖3 i=2,k=1時動力不穩定區域Fig.3 Dynamic unstable regions(i=2,k=1)

圖4 i=3,k=1時動力不穩定區域Fig.4 Dynamic unstable regions(i=3,k=1)

結論2(由圖2(a)、圖3、圖4得到):

若周向半波數k不變,軸向半波數i增大時,圓柱薄殼艙段的前三階不穩定區整體上移,且不穩定區域的寬度逐漸減小,進一步計算k取2或3或更大的數值,且i增大時也能得到相同的結論,只是前三節動力不穩定區整體上移的幅度比k=1時要??;但是,若取軸向半波數i不變,周向半波數k增大時,并不能得到以上結論。因此,在圖1的逐步搜索法中,以周向半波數k作為外循環,軸向半波數i作為內循環是合理的。

3.2 動力屈曲可靠性計算結果及分析

經過試算,高階振型對應的動力不穩定區域不僅頻率值較高而且寬度較小,實際結構在運動過程中,其載荷的頻率值達不到那么高,這些不穩定區域對結構的動力穩定性影響很小,可以忽略。計算表明,圖1的逐步搜索法中,取a=5,b=5,βc=5.0是合理的。以f、V和δ作為基本隨機變量進行動力屈曲可靠性分析,考慮這三個隨機變量相互獨立且服從正態分布,它們的隨機參數見表1。

表1 隨機變量取值

輸入表1中的數據,采用圖1所示的逐步搜索法,搜索得到的有效的不穩定區域數目為1,圓柱薄艙段的動力屈曲可靠度指標β=3.266 1。

為了研究三個隨機變量的變異程度對動力屈曲可靠性指標的影響,現以表1中隨機參數為基準,設為X0,各隨機參數的均值不變,將要研究的隨機變量的變異系數取0.75X0-2.0X0,其余隨機變量的變異系數仍為X0,得到如表2所示的動力屈曲可靠性指標隨各隨機參數變異系數的變化值,“( )”中的數字為采用逐步搜索法得到的有效不穩定區域的個數。

表2 動力屈曲可靠性指標隨變異系數的變化值

結論3(由表2得到):

① 載荷比例系數的變異系數對動力屈曲可靠性指標的影響極小,可以忽略;

② 軸向載荷圓頻率的變異系數越大,動力屈曲可靠性指標越小;

③ 速度的變異系數越大,動力屈曲可靠性指標越小。

為了研究三個隨機變量的均值對動力屈曲可靠性指標的影響,現以表1中隨機參數為基準,設為Y0,各隨機參數的變異系數不變,將要研究的隨機變量的均值取表3中所示的值,其余隨機變量的均值仍為Y0,得到如表3所示的動力屈曲可靠性指標隨各隨機參數均值的變化值,“()”中的數字為采用逐步搜索法得到的有效不穩定區域的個數。

表3 動力屈曲可靠性指標隨均值的變化值

結論4(由表3得到):

① 載荷比例系數的均值對動力屈曲可靠性指標的影響極小,可以忽略;

② 軸向載荷圓頻率的均值位于110 Hz,150 Hz附近時,動力屈曲可靠性指標相對較大,而位于100 Hz,130 Hz附近時,動力屈曲可靠性指標相對較小,這說明軸向載荷圓頻率取某些范圍內的值時,結構動力屈曲可靠性指標相對較大,結構較安全,且這些范圍并不連續。

③ 速度的均值小于400 m/s時,動力屈曲可靠性指標較大,結構相對較為安全,而速度的均值大于450 m/s時,結構會失效,因此,速度取某些范圍內的值時,結構較為安全。

4 結 論

從超空泡運動體的動力屈曲分析及可靠性分析的算例得出如下結論:

(1) 速度和載荷比例系數增大,不穩定的載荷頻率值減小,且第一階不穩定載荷頻率減小的速度明顯大于后兩階;

(2) 周向半波數不變,軸向半波數增大,圓柱薄殼艙段的前三階不穩定區整體上移,且不穩定區域的寬度逐漸減小;

(3) 速度和軸向載荷圓頻率的變異系數越大,動力屈曲可靠性指標越小,載荷比例系數的變異系數對動力屈曲可靠性指標的影響極小,可以忽略;

(4) 載荷比例系數的均值對動力屈曲可靠性指標的影響很小,可以忽略,軸向載荷圓頻率和速度的均值取某些范圍內的值時,結構動力屈曲可靠性指標相對較大,結構較安全。

綜上所述,本文的研究為如何選擇載荷頻率、速度和載荷比例系數的安全范圍提供了理論依據。還有利于計算總系統的可靠度指標。

[ 1 ] Ruzzene M. Dynamic buckling of periodically stiffened shells: application to supercavitating vehicles [J].International Journal of Solids and Structures.2004,41(3-4): 1039-1059.

[ 2 ] Edward A, Vipperla V,Ramana G,et al. Structural response and optimization of a supercavitating torpedo [J].Finite Elements in Analysis and Design.2005,41(6):563-582.

[ 3 ] 施連會,王安穩. 軸向載荷下超空泡航行體動力穩定性的數值研究[J]. 振動與沖擊,2011,30(2):55-59. SHI Lian-hui, WANG An-wen.Numerical study on dynamic stability of supercavitating vehicles subjected to axial loads [J]. Journal of Vibration and Shock, 2011, 30(2): 55-59.

[ 4 ] 宋向華,安偉光,劉明. 軸向周期載荷下超空泡射彈的動力穩定性分析[J]. 哈爾濱工程大學學報,2012,33(10):1238-1243. SONG Xiang-hua,AN Wei-guang,LIU Ming. Dynamic stability analysis of supercavitating projectile subjected to axial periodic load [J]. Journal of Harbin Engineering University, 2012,33(10):1238-1243.

[ 5 ] 顧永維,安偉光,安海. 水下高速運動體的抗彎穩定可靠性分析[J]. 哈爾濱工程大學學報,2008,29(7):683-686. GU Yong-wei, AN Wei-guang,AN Hai. Analyzing of buckling resistance of an underwater elasticbody experiencing high speed motion [J]. Journal of Harbin Engineering University, 2008, 29(7): 5683-686.

[ 6 ] 周凌,安偉光,安海. 水下通氣超空泡航行體結構屈曲可靠性分析[J]. 工程力學,2010,27(9):248-256. ZHOU Ling, AN Wei-guang,AN Hai.Structural buckling reliability analysis of underwater ventilated supercavitating vehicles [J]. Engineering Mechanics, 2010, 29(7): 248-256.

[ 7 ] 周凌,安偉光,安海.超空泡運動體強度與穩定性的非概率可靠性分析[J]. 哈爾濱工程大學學報,2009,30(4):362-367. ZHOU Ling,AN Wei-guang,AN Hai.Non-probabilistic reliability analysis of supercavitating vehicles based on structure strength and buckling[J].Jouranl of Harbin Engineering University,2009,30(4):362-367.

[ 8 ] Ahn S S, Ruzzene M. Optimal design of cylindrical shells for enhanced buckling stability:Application to supercavitating underwater vehicles[J]. Finite Elements in Analysis and Design,2006,42(11):967-976.

[ 9 ] 符拉索夫 瓦 札.殼體的一般理論[M].北京:人民教育出版社,1960:262-271.

[10] 符華,鮑洛金.彈性體系的動力穩定性[M].北京:高等教育出版社,1960:26-31.

[11] 貢金鑫. 工程結構可靠度計算方法[M].大連:大連理工大學出版社,2003:263-268.

Dynamic buckling and reliability analysis of a cylindrical thin shell for supercavitating vehicles

WANG Jie-fang, AN Wei-guang, SONG Xiang-hua

(Department of Aerospace Engineering, Harbin Engineering University, Harbin 150001, China)

Firstly, a supercavitating vehicle was modeled as a cylindrical thin shell loaded with a axial dynamic load. The dynamic stability differential equations and the unstable regions of the shell were deduced. Secondly, the safety margin equations for its dynamic buckling were given and linearized with the limited step length iteration method considering the randomness of axial load. Then, the step-by-step searching method was proposed to search those effective safety margin equations. Finally, the reliability index of the dynamic buckling was calculated with the step-by-step equivalent plane method. Through numerical examples, the influences of change of load frequency, velocity and proportional coefficient of load on the dynamic buckling reliability were analyzed, respectively. The calculation results provided a theoretical basis for choosing the safety range of load frequency, velocity and proportional coefficient of load.

dynamic buckling; reliability; limited step length iteration method; step-by-step searching method; step-by-step equivalent plane method

國家科技部國際合作專項(2012DFR00070)

2013-01-28 修改稿收到日期:2013-05-22

王杰方 女,博士,1987年生

O342

A

10.13465/j.cnki.jvs.2014.08.005

猜你喜歡
區域
分割區域
探尋區域創新的密碼
科學(2020年5期)2020-11-26 08:19:22
基于BM3D的復雜紋理區域圖像去噪
軟件(2020年3期)2020-04-20 01:45:18
小區域、大發展
商周刊(2018年15期)2018-07-27 01:41:20
論“戎”的活動區域
敦煌學輯刊(2018年1期)2018-07-09 05:46:42
區域發展篇
區域經濟
關于四色猜想
分區域
公司治理與技術創新:分區域比較
主站蜘蛛池模板: 国产美女精品在线| 男女性午夜福利网站| 在线精品自拍| 国产精品lululu在线观看| 国产小视频a在线观看| 国产麻豆永久视频| 韩日无码在线不卡| 亚洲有无码中文网| 直接黄91麻豆网站| 久久天天躁夜夜躁狠狠| 伊人五月丁香综合AⅤ| 国产在线欧美| 91精品国产一区自在线拍| 欧美一级大片在线观看| 欧美激情综合一区二区| 久久精品一卡日本电影| 国产黄色视频综合| 性色在线视频精品| 亚洲一区二区日韩欧美gif| 亚洲国产成人久久77| 亚洲A∨无码精品午夜在线观看| 亚洲日本www| 人妻一区二区三区无码精品一区| 精品偷拍一区二区| 国产97视频在线| 91福利一区二区三区| 久久精品欧美一区二区| 亚洲中文无码h在线观看| 久久久久88色偷偷| 亚洲天堂免费| 欧美日韩在线亚洲国产人| 青青草国产一区二区三区| 国产成人一级| 亚洲男人的天堂在线观看| 亚洲天堂2014| 噜噜噜久久| 日韩毛片免费| 伊人久久婷婷| 婷婷成人综合| 日韩二区三区无| 国产啪在线91| 伊人查蕉在线观看国产精品| 日韩黄色精品| 成人在线亚洲| 亚洲综合第一页| 色噜噜综合网| 伊人久久婷婷五月综合97色| 色综合网址| 亚洲天堂精品在线观看| 日本在线国产| 黄色免费在线网址| 九九久久99精品| 亚洲国产看片基地久久1024| 嫩草国产在线| 又爽又大又黄a级毛片在线视频 | 欧美一级特黄aaaaaa在线看片| 99精品这里只有精品高清视频| 国产网友愉拍精品| 亚洲中文在线看视频一区| 一本久道热中字伊人| 天堂av综合网| 国产精品久久久久鬼色| 91麻豆精品国产91久久久久| 香蕉蕉亚亚洲aav综合| 欧美黄网在线| 成人日韩欧美| 亚洲国内精品自在自线官| 黄色网站不卡无码| www.亚洲天堂| 欧美一区中文字幕| 亚洲国产成熟视频在线多多| 噜噜噜久久| 98精品全国免费观看视频| 亚洲欧美日韩成人高清在线一区| 全裸无码专区| 欧美成人看片一区二区三区| 亚洲欧美日韩成人高清在线一区| 超薄丝袜足j国产在线视频| 日韩欧美国产另类| 999福利激情视频| 在线国产91| 国产精品性|