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

考慮老化時間影響的隔震支座橡膠本構Moony-Rivlin模型常數研究

2016-09-18 02:45:44李艷敏馬玉宏羅佳潤趙桂峰
振動與沖擊 2016年16期
關鍵詞:橡膠材料有限元模型

李艷敏, 馬玉宏, 羅佳潤, 趙桂峰

(1.廣州大學 工程抗震研究中心,廣州 510405;2. 上海筑房工程科技有限公司,上?!?00062;3.廣州大學 土木工程學院,廣州 510006)

?

考慮老化時間影響的隔震支座橡膠本構Moony-Rivlin模型常數研究

李艷敏1, 馬玉宏1, 羅佳潤2, 趙桂峰3

(1.廣州大學 工程抗震研究中心,廣州510405;2. 上海筑房工程科技有限公司,上海200062;3.廣州大學 土木工程學院,廣州510006)

目前Moony-Rivlin模型被廣泛應用于橡膠隔震支座的有限元分析中,其中,材料常數的標定是關鍵環節。為充分了解老化作用對橡膠材料常數的影響,采用最小二乘擬合法計算老化試驗后的Moony-Rivlin模型材料常數,擬合出材料常數隨老化時間的變化關系式,提出了修正Moony-Rivlin模型,并將擬合值和計算值進行對比。最后,根據上述方法確定的材料常數用Abaqus軟件對隔震橡膠支座老化前后的性能進行了模擬分析,并將有限元分析結果與試驗結果進行對比。結果表明:隨著老化時間的增加,Moony-Rivlin模型常數逐漸增大;老化前、后橡膠隔震支座的豎向剛度、水平剛度試驗測試與有限元模擬值的誤差均在可接受的范圍內,驗證了所建立修正本構模型的準確性,可為橡膠隔震支座老化性能的研究和橋梁等結構全壽命性能設計提供理論支撐。

Moony-Rivlin模型常數;橡膠;老化;有限元分析

早在19世紀,橡膠材料因其良好的黏彈性,被廣泛用作飛機、火車、汽車、船舶和建筑物的減震部件。但是,因為橡膠材料的復雜分子特性、材料和幾何的雙重非線性以及溫度、荷載時間、應變量等多種因素的復雜影響,所以橡膠材料力學性能的計算就顯得十分困難[1-3]。有限元分析方法的廣泛推廣為橡膠制品的工程模擬提供了廣闊的發展前景。而有限元模擬分析結果的準確性與所采用橡膠本構關系模型以及模型中材料常數的準確性有著密切聯系。隨著跨海橋梁的不斷發展,橡膠隔震支座在復雜海洋環境中的老化現象較為嚴重,但橡膠隔震支座有限元分析中普遍采用的橡膠本構模型Moony-Rivlin常數是否會發生變化,以及變化規律如何是目前急需解決的科學問題。

老化對橡膠材料性能的影響一直受到國內外學者的關注,并建立了許多不同的理論本構模型:RIVLIN[2]認為橡膠各向同性,且拉壓性質相同,將應變能密度函數改成了級數的形式,為Moony-Rivlin模型的發展提供了理論條件;NAMJOO等[4]在用數值模擬輪胎、土壤的相互作用時采用了三維有限元模型,其中采用Moony-Rivlin模型分析了不可壓縮橡膠輪胎的性能;余超等[5]在研究熱空氣和海水對橡膠的老化行為時,擬合出了兩種環境下力學性能指標與老化時間的關系;左亮等[6]采用理論估算確定Moony-Rivlin模型材料常數;黃建龍等[7]采用Ansys有限元分析軟件對橡膠材料的Moony-Rivlin模型和Yeoh模型進行分析對比研究;左亮等[8]在橡膠小變形范圍內,采用理論推導的方法,得到在相同硬度下橡膠Moony-Rivlin 模型材料系數的關系式;仲健林等[9]通過開展橡膠材料單軸拉伸試驗,并結合Exp-ln超彈性體本構模型和廣義黏彈性方法,提出了一種描述橡膠材料不同應變率下力學響應的超彈性本構模型。以上研究大多側重理論推導,缺少較好的試驗依據,難以充分考慮老化作用給橡膠材料性能產生的影響,更不能很好地解決老化作用引起的橡膠材料性能隨時間變化的問題。

為解決上述問題,本文以Moony-Rivlin模型為例,根據老化試驗數據,采用最小二乘擬合法對橡膠在老化試驗過程中獲得的材料常數進行了分析,得出了橡膠材料常數隨老化時間的時變規律。最后,通過有限元模擬與試驗結果的對比分析,驗證本文研究成果的準確性,為今后深入研究復雜海洋環境中的橡膠隔震支座性能變化規律提供理論基礎,為近海隔震結構的設計提供參考。

1 Moony-Rivlin模型的本構關系

目前,國內外學者主要采用統計理論和唯象理論這兩種方法來研究橡膠材料本構關系中的應變能密度函數[1,10-12]。本文主要采用唯象理論。

唯象理論在處理橡膠材料時認為橡膠材料的變形都是各向同性的超彈性材料的均勻變形。隨著橡膠材料本構關系研究的不斷發展,建立了很多不同的本構模型,這些模型都屬于應變能函數的某種特殊形式。橡膠材料應變能函數有兩種表達方式,即用變形張量的三個不變量I1、I2和I3來表示的應變能函數W(I1,I2,I3);用主伸長比λ1、λ2和λ3的關系式來表示的W(λ1,λ2,λ3)。三個變形張量不變量與三個主伸長比的關系式如下:

I1=λ12+λ22+λ32

(1)

I2=λ12λ22+λ22λ32+λ12λ32

(2)

I3=λ12λ22λ32

(3)

λi=1+εi

(4)

式中:Ii為變形張量不變量;λi為定義的伸長比;εi為主軸方向的應變,i=1,2,3

橡膠材料這種超彈性材料的應變能函數轉化為多項式形式后,可由應變偏量能和體積應變兩部分之和構成,表示如下:

(5)

式中:N為多項式階數;J為彈性體積比;Di為定義材料的壓縮性,Di的值決定材料是否可壓縮,如果所有N=1都為0,則材料是完全不可壓縮。

對于多項式,無論N取何值,橡膠的初始剪切模量μ0、初始體積模量k0都取決于多形式一階(N=1)系數,即μ0=2(C10+C01),k0=2/D1。

對于完全多項式,如果N=1則只有線性部分的應變能保留下來,即Moony-Rivlin模型:

(6)

式中:W為應變勢能;I1、I2為變形張量;C10、C01、Di為材料常數;J為彈性體積比(代表模型的體積應變)。

2 橡膠材料老化試驗

在對天然橡膠隔震支座開展老化試驗的同時,設計了與支座相同材料的橡膠片,將支座與橡膠片放置在同樣的試驗環境下開展老化試驗,研究橡膠材料物理力學性能在老化作用下隨時間的變化規律。得到了大量的老化試驗數據,為本文的研究工作奠定了良好的基礎。

2.1試驗概況

在對橡膠支座進行老化試驗時,在同一個老化試驗箱里同時晾掛了橡膠片,且在試驗期間對橡膠片每隔2天取樣一次,以查看橡膠材料性能隨老化時間的變化情況。而進行老化試驗之前,先對橡膠片取樣并對12組橡膠隔震支座做了相應的基本性能測試試驗,然后,將其放置于80℃恒溫老化箱中進行20天(480 h)熱老化試驗,測試橡膠材料硬度、定伸應力、拉伸強度及扯斷伸長率在老化作用下隨時間的變化規律。橡膠隔震支座及橡膠片晾掛情況見圖1[13]。

圖1 橡膠隔震支座及橡膠塊體放置圖Fig.1 The place of rubber isolated bearing and rubber piece

2.2試驗結果

針對單獨老化480 h的定伸應力和硬度的測試數據,采用最小二乘擬合法計算相應的橡膠本構模型常數值。依照中位數的原則來選取數據,具體數據如圖2~圖5所示。

圖2 老化時間對橡膠50%、100%定伸應力的影響Fig.2 The influence of aging time on stretch stress at 50%、100% strain

3 用最小二乘法確定橡膠材料系數

3.1應力-應變關系推導

材料的主要特性通常由應力-應變關系來表征,而橡膠材料的應力應變關系可以通過應變能密度函數對其主伸長比求偏導來表示,此應力稱為Piola-Kirchhoff應力,應變稱為Cauchy-Green應變,其應力-應變形式如下[14]:

(7)

由式(1)~(6)得主軸力ti和主伸長比λi之間的關系:

(8)

式中:t1、t2、t3為三個主應力。

對于單軸拉伸試驗,有t2=t3=0,得到

λ22=λ32=1/λ1

(9)

對于不可壓縮材料則有:

I3=λ12λ22λ32=1

(10)

由式(9)、(10)得到不可壓縮橡膠材料的主應力和主伸長比與不變量的關系式如下:

(11)

式中:C10=?W/?I1、C01=?W/?I2。

3.2材料常數的確定

為了確定材料常數C10和C01,假定橡膠是完全不可壓縮材料。根據試驗所得到的定伸應力和應變數據應用最小二乘擬合法求出參數C10和C01。首先,令:

A1=2(λ12-1/λ1),B1=t1/A1,E1=1/λ1

則函數式(11)可改寫成:

B1=C10+E1C01

(12)

根據試驗數據λi和ti求出相應的A1、B1和E1,并得到老化前后的應力-應變關系,見圖6、圖7。

根據最優平方逼近式的正規方程組[14]:

(13)

式中:(,)是內積,基函數(x)=xj,cj是系數(下同)。所逼近的近似多項式為:

(14)

式中:基函數φj(x)=xj。

用矩陣的形式表示為:

式中:(φi,φj)為內積,(i=0,1,2…n;j=0,1,2…n)。

由上述應力-應變數據,分別計算出老化前與老化后的相關數據。然后,由矩陣方程組:

(16)

計算出老化0~480 h工況下的橡膠材料常數C10和C01,并根據計算結果擬合出橡膠材料常數隨老化時間的變化情況,詳見圖6、圖7。

圖6 老化時間對橡膠材料常數C10的影響Fig.6 The influence of aging time on rubber material constant C10

圖7 老化時間對橡膠常數的影響Fig.7 The influence of aging time on rubber material constant absolute value

式(6)中的D1值可以根據橡膠的初始體積模量k0=2/D1來確定,k0為橡膠材料的初始體積模量。而橡膠支座的的體積模量通常根據表2來確定。

表中,E0為橡膠彈性模量,近似等于3G;G為橡膠剪切模量;E∞為橡膠體積彈性模量;k為與橡膠硬度有關的彈性模量修正系數。

由圖4可知橡膠硬度的初始值為43,老化480 h(小時)后的硬度值為47,根據表1中的數據采用插值方法計算出相應的體積模量,并計算出相應的D1值和1/D1值,計算結果如表3所示。

表1 材料常數擬合值和計算值的誤差對比

表2 支座橡膠材料性能參數[15]

表3 老化工況下計算的壓縮性常數

由表2的數據來看,老化前、后橡膠材料的壓縮性常數變化不大,后續研究中忽略壓縮性常數D1的變化,僅考慮材料常數C10和C01的變化。

因此,隨著老化時間t的變化,應變能密度函數公式近似擬合為:

W=0.375 1e0.000 5t(I1-3)+

(17)

式中:C10=0.375 1e0.000 5t,C01=-0.202 2e0.000 5t,t為試驗時間(h)。

實際使用時間treal與試驗時間t可以根據下列公式進行轉換:

(18)[13]

式中,Ea反應活化能(kJ/mol·K),取95 kJ/mol;R為氣體常數(取8.3[J/mol·K]),其與氣體類型無關,僅與量綱有關;Treal為實際環境中的絕對溫度(K),取293 K;Ttest為熱氧老化試驗的絕對溫度(K),取353 K;Treal為實際老化時間;t為試驗時間。

由式(17)可見,考慮老化影響后橡膠Moony-Rivlin模型常數會隨著老化時間的變化而發生變化,本構模型能夠體現老化時間的變化規律,可為橡膠隔震支座老化性能的研究提供理論支撐。

4 有限元分析與試驗結果對比分析[13]

4.1老化前有限元分析與試驗結果對比

以橡膠支座試驗結果作為與有限元分析對比的參考標準,有限元分析中參數采用上述擬合公式算出的Moony-Rivlin參數。

由于對橡膠隔震支座豎向剛度試驗時未記錄力-位移曲線,本文只展示有限元模擬的力-位移曲線。豎向力取隔震支座下封板表面193個節點的豎向反力之和,位移取隔震支座上封板表面193個節點豎向位移的平均值,模型圖見圖8。老化前后的豎向剛度見圖9,水平剛度則取5#支座的實測曲線與有限元分析結果進行對比,見圖10。

圖8 橡膠隔震支座有限元模型圖Fig.8 The finite element model of isolation rubber bearing

取12個橡膠隔震支座性能試驗值的平均值與有限元模擬值進行對比。支座豎向剛度的試驗結果為188.66 kN/mm,有限元模擬結果為184.56 kN/mm,誤差為2.17%;水平剛度的試驗結果為0.189 kN/mm(未考慮溫度修正),0.193 kN/mm (考慮溫度修正),有限元模擬結果為0.181 7 kN/mm,誤差分別為3.86%、5.85%,可見,模擬精度比較高。因此,下文研究中按照同樣方法確定Moony-Rivlin本構模型參數。

4.2老化后試驗結果與有限元模擬對比

與前述方法類似,老化后的豎向剛度有限元分析結果見圖9,老化后水平剛度測試值與有限元對比見圖11。

老化后隔震支座豎向剛度的試驗數據:233.79 kN/mm,有限元模擬老化后的豎向剛度為210.16 kN/mm,誤差為10.11%,這個誤差基本認為在接受的范圍內;老化后水平剛度試驗數據為0.215 4 kN/mm(未考慮溫度修正),0.211 4 kN/mm(考慮溫度修正),有限元模擬數據為0.234 9 kN/mm,誤差分別為9.09%、11.12%。同時,可以看出:有限元模擬老化后的豎直剛度比老化前的剛度增大了13.87%;而有限元模擬老化后的水平剛度比老化前剛度增大了29.28%。

4.3誤差分析

(1) 本文中采用的Moony-Rivlin本構模型參數是根據試驗數據計算并擬合之后的結果,試驗數據本身存在一定的誤差。

(2) 用有限元模擬時,采用的是理想模型,沒有考慮實際生產過程及試驗工裝方面的影響,也會造成一定的誤差。

(3) 有限元模擬老化后的橡膠支座性能分析時,沒有充分考慮橡膠厚度對材料參數的影響,認為其老化程度相同,這也是產生誤差的一個原因。

5 結 論

本文分析了老化作用對橡膠材料常數C10和C01的影響,根據老化試驗數據,采用最小二乘擬合法計算材料常數C10、C01和|C01|,并擬合出老化時間與材料常數C10和材料常數|C01|的關系,且計算了擬合值與計算值的誤差范圍;進而得到了考慮了老化作用影響的橡膠材料修正Moony-Rivlin本構模型。最后利用修正橡膠本構模型,采用有限元軟件Abaqus對橡膠隔震支座進行了模擬分析,并與試驗結果進行對比。主要結論如下:

(1) 老化等因素會對橡膠Moony-Rivlin模型常數產生影響;

(2) 隨著老化時間的增加,橡膠Moony-Rivlin模型常數隨著老化時間的增加而增大;

(3) 利用修正的Moony-Rivlin橡膠本構模型進行有限元模擬與試驗數據吻合良好,證明所建立的本構模型合理準確;

(4) 所建立的修正Moony-Rivlin橡膠本構模型可以充分考慮老化時間的影響,可為橡膠隔震支座老化性能的研究和橋梁等結構全壽命性能設計提供理論支撐。

本文已經對橡膠隔震支座橡膠本構模型Moony-Rivlin模型常數隨老化時間的變化規律開展了研究,并根據試驗數據擬合出了材料常數隨老化時間的變化關系式,得出了一些規律。但是,隨著隔震支座在跨海橋梁等近海工程中的推廣,由于復雜的海洋環境,仍有一些問題需要更深入的研究:

(1)考慮海蝕環境對隔震支座橡膠本構模型Moony-Rivlin模型常數的影響,結合試驗結果,計算并總結出模型常數隨海蝕的變化規律;

(2)綜合考慮老化時間和老化厚度兩方面的因素,得出更為接近實際老化工況的變化規律;

(3)考慮老化和海蝕的綜合作用,來擬合材料常數隨環境的變化關系式。

[1] GENT A N.Engineer with rubber how to design rubber: components[M]. Munich:Hanser Publishers,2000:259-275.

[2] RIVLIN R S.Large elastic deformations of isotropic materials[C]//Philosophical Transactions of the Royal Society of London, Physical Sciences and Engineer. London,UK:1948.

[3] 鄭明軍,王文靜,陳政南,等.橡膠Moony-Rivlin模型力學性能常數的確定[J].橡膠工業,2003,50(8):462-463.

ZHENG Mingjun,WANG Wenjing,CHEN Zhengnan,et al.Determination for mechanical constants of rubber Moony-Rivlin model[J]. Rubber Industry, 2003,50(8):462-463.

[4] NAMJOO M,GOLBAKHSSELN H. Numerical simulation of tire/soil interaction using a verified 3D finite element model [J].J Cent South Unlv,2014,21:817-821.

[5] 余超,文慶珍,余宏偉,等.丁苯橡膠在熱空氣和海水中老化性能的比較[J].合成橡膠工業,2010(1):56-59.

YU Chao,WEN Qingzhen,YU Hongwei, et al.Comparision of properties of styrene-butadiene rubber aged in sea water and in hot air[J].China Synihetc Rubber Ndusiry,2010(1):56-59.

[6] 左亮,肖緋雄.橡膠Moony-Rivlin模型常材料系數的一種確定方法[J].機械制造,2008,46(7):38-40.

ZUO Liang,XIAO Feixiong.One method of determination for material constants of rubber Moony-Rivlin model[J]. Machine Building, 2008,46(7):38-40.

[7] 黃建龍,解廣娟,劉正偉.基于Moony-Rivlin和Yeoh模型的超彈性橡膠材料有限元分析[J].橡塑技術與設備,2008,34(12):22-26.

HUANG Jianlong,XIE Guangjuan,LIU Zhengwei. Finite element analysis of super-elastic rubber materials based on the Moony-Rivlin and Yeoh model[J]. Rubber/Plastics Technology and Equipment, 2008, 34(12):22-26.

[8] 左亮,肖緋雄.橡膠Moony-Rivlin模型材料系數對軸向剛度影響分析[J].彈性體,2008,18(3):54-56.

ZUO Liang,XIAO Feixiong.Analysis of materials coefficient of rubber Moony-Rivlin model impacting on the axial stiffness[J].China Elastomerics, 2008,18(3):54-56.

[9] 仲健林,任杰,馬大為.基于Exp-ln模型與廣義黏彈性理論的橡膠本構模型及其應用研究[J].振動與沖擊,2015,34(19):150-156.

ZHONG Jianlin,REN Jie,MA Dawei.Constitutive model and its application for rubber material based on Exp-In model and generalized viscoelastic theory[J].Journal of Vibration and Shock,2015,34(19):150-156.

[10] 朱武.某重型卡車疊層橡膠彈簧有限元分析及穩健設計[D].湖南:湖南大學,2011.

[11] 卜繼玲,黃友劍.軌道車輛橡膠彈性元件設計計算算法[M].北京:中國鐵道出版社,2010.

[12] 劉萌,王青青,王國權.橡膠Moony-Rivlin模型中材料常數的確定[J].橡膠工業,2011,58(4):241-245.

LIU Meng, WANG Qingqing, WANG Guoquan. Determination for material constants of rubber Moony-Rivlin model[J].Rubber Industry, 2011,58(4):241-245.

[13] 羅佳潤.海蝕環境下橡膠隔震支座性能劣化規律研究[D].廣州:廣州大學,2014.

[14] 黃慶專.Moony-Rivlin模型及其系數最小二乘解法[J].熱帶農業科學,2009,29(5): 20-24.

HUANG Qingzhuan.Rubber Moony-Rivlin model and its modulus’s least-square solution[J].Chinese Journal of Tropical Agriculture, 2009,29(5): 20-24.

[15] GB 20688.1—2007.橡膠支座:第3部分:建筑隔震橡膠支座[S] .北京:中國標準出版社,2007.

The effect of aging on the material constant of the rubber isolator’s constitutive model Moony-Rivlin

LI Yanmin1, MA Yuhong1, LUO Jiarun2, ZHAO Guifeng3

(1. Engineering Seismic Research Centre of Guangzhou University, Guangzhou 510405, China;2. Shanghai Construction Engineering Technology Co., Ltd., Shanghai 200062, China;3. School of Civil Engineering of Guangzhou University, Guangzhou 510006, China)

The Moony-Rivlin model is widely applied in finite element analysis of rubber isolators, in which a key point is to determine material parameters of the model. In order to study the effect of aging on material constant of rubber intensively, the material constants of Moony-Rivlin model were obtained through investigating aging test result of rubber by the Least-square method. The relationship between the material constants and the aging time was derived and a modified model was achieved. The fitting values were compared with the calculating ones. Finally, the performance of the rubber isolators was simulated and studied by using the Abaqus software and the material constants were determined from the above method. The results of finite element analysis were compared with the test results. The results show that: Moony-Rivlin model constant increases gradually with aging time; the errors between the test values and the finite element simulation values for the vertical stiffness and the horizontal stiffness are acceptable, whether before aging or aged, which proves the accuracy of the modified constitutive model. This work provides theoretical support for the performance study of rubber isolated bearings under the aging environment and the life-cycle performance design about bridges and other structures.

Moony-Rivlin model constant; rubber; aging; finite element analysis

國家重點基礎研究發展計劃973項目(2011CB013606);國家自然科學基金項目(51578170);國家自然科學基金高鐵聯合基金重點項目(U1334209);長江學者和創新團隊發展計劃項目(IRT13057);廣州市屬高校科技計劃項目(1201421152)

2015-12-07修改稿收到日期:2016-03-03

李艷敏 女,碩士生,1992年生

馬玉宏 女,博士,研究員,1972年生E-mall:849502749@qq.com

P315.966

A

10.13465/j.cnki.jvs.2016.16.026

猜你喜歡
橡膠材料有限元模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
一種鞋底橡膠材料
橡膠工業(2016年12期)2016-02-23 16:41:08
橡膠材料單軸拉伸疲勞壽命預測的有限元分析
橡膠工業(2015年7期)2015-08-29 06:33:40
一種能消除擠出口模上的滯留物的橡膠材料
橡膠工業(2015年11期)2015-08-01 09:08:54
一種用于橡膠材料自修復的微膠囊的制備方法
橡膠工業(2015年4期)2015-07-29 09:17:20
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 久久久久亚洲AV成人人电影软件| 欧美亚洲一区二区三区在线| 色综合婷婷| 亚洲福利视频网址| 日韩午夜伦| 国产成人资源| 一本久道热中字伊人| 亚洲无线一二三四区男男| 东京热av无码电影一区二区| 日韩黄色大片免费看| 久久国产精品麻豆系列| 毛片网站在线看| 91小视频版在线观看www| 亚洲欧美一区二区三区图片| 成人福利在线观看| 国产福利免费在线观看| 国产欧美又粗又猛又爽老| 午夜天堂视频| 亚洲AV无码乱码在线观看裸奔| 亚洲综合狠狠| 任我操在线视频| 国产一级二级在线观看| 欧美成人精品一区二区| 日韩色图在线观看| 在线观看亚洲成人| 成年女人18毛片毛片免费| 免费在线视频a| 亚洲AV无码精品无码久久蜜桃| 亚洲天堂视频在线免费观看| 亚洲综合香蕉| 亚洲成av人无码综合在线观看 | 老司机午夜精品网站在线观看| 国产91特黄特色A级毛片| 中文字幕在线日本| 熟妇人妻无乱码中文字幕真矢织江| 青青草原国产精品啪啪视频| 二级毛片免费观看全程| 在线观看亚洲国产| 国产精品成人第一区| 国产精品免费露脸视频| 四虎影视国产精品| 国产亚洲精久久久久久久91| 免费人成网站在线高清| 色偷偷一区二区三区| 日本高清免费不卡视频| 欧美一级在线| 国产精品三区四区| 色视频国产| 97狠狠操| 国产av一码二码三码无码| 国产精品毛片一区| 伊人大杳蕉中文无码| 18黑白丝水手服自慰喷水网站| 亚洲第一福利视频导航| 99久久免费精品特色大片| 亚洲精品无码AⅤ片青青在线观看| 国产又爽又黄无遮挡免费观看 | 国产在线一区视频| 四虎综合网| 在线观看精品自拍视频| 999国内精品久久免费视频| 91蜜芽尤物福利在线观看| 日本手机在线视频| 国产大片黄在线观看| 韩国v欧美v亚洲v日本v| 亚洲欧洲日韩综合| 精品视频第一页| 国产精品不卡片视频免费观看| 丁香婷婷激情综合激情| 色综合热无码热国产| 99精品热视频这里只有精品7| 久久这里只有精品8| 四虎国产在线观看| 国产精品短篇二区| 久久福利片| 五月婷婷中文字幕| 亚洲国产欧洲精品路线久久| 波多野结衣一二三| 57pao国产成视频免费播放 | 亚洲精品视频网| 婷婷久久综合九色综合88| 东京热av无码电影一区二区|