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

顆粒增強橡膠復合材料有效力學性能預測分析

2017-08-01 00:03:28宿曉如馮春冬羅冬梅
兵器裝備工程學報 2017年7期
關鍵詞:有限元法力學性能復合材料

謝 悅,宿曉如,馮春冬,羅冬梅

(佛山科學技術學院 土木工程系, 廣東 佛山 528000)

?

顆粒增強橡膠復合材料有效力學性能預測分析

謝 悅,宿曉如,馮春冬,羅冬梅

(佛山科學技術學院 土木工程系, 廣東 佛山 528000)

采用幾種經典預測計算方法分析顆粒半徑大小、顆粒百分比、顆粒彈性模量等參數對隨機顆粒增強橡膠復合材料有限元模型有效力學性能的影響。結果表明:高百分比、小半徑、低彈性模量的顆粒有利于改善隨機顆粒增強橡膠復合材料的有效力學性能;所用方法中,經典Mori-Tanaka法、Halpin-Tsai法和以應變能法和廣義虎克定律為基礎的有限元法都有一定的局限性,而修正的經典Halpin-Tsai公式能有效改善經典公式預測的結果,另一種多尺度均質化方法的結果介于數值法與經典公式法之間,說明修正Halpin-Tsai法和均質化法的預測更為可靠。

顆粒;橡膠;復合材料;有效力學性能;均質化法

橡膠基復合材料因其彈性大、可設計性強、可以人工合成、價格低廉等優良性能廣泛應用于航空航天、土木建筑、交通運輸及體育器材等領域,通過添加纖維、顆粒等形狀和性能不一樣的增強相能夠進一步改善橡膠的性能。目前對于這類材料的研究一般是通過大量的物理實驗進行配方設計。李福強等[1]研究了未處理棉短纖維、尼龍短纖維和木質纖維素短纖維的用量對短纖維增強三元乙丙橡膠復合材料物理力學性能的影響;謝尊虎等[2]分析了提高硅橡膠各項性能的主要途徑和方法,指出了提高硅橡膠相關性能的發展方向;王作齡等[3]通過編譯橡膠試驗方法,對橡膠的實驗材料尺寸、試驗方法、試驗數據處理做了統一,為后續工作者提供試驗參考,齊海波等[4]通過試驗,利用電子顯微鏡對顆粒增強復合材料中顆粒長徑比的影響進行統計分析,推導出顆粒長徑比與其體積分數的關系,并將其運用到復合材料有效彈性模量計算的Eshelby等效夾雜理論及改進的自洽法中;研究結果表明:顆粒體積分數較低時,從Eshelby等效夾雜理論得到的結果具有較高的精度;體積分數較高時,改進的自洽法能作出較準確的預測;柏振海等[5]通過大量實驗找出復合材料有效彈性模量和體積分數的關系式,利用Voigt等應變假設和Reuss等應力假設成功引入J因子,預測各種不同復合材料體系的彈性模量。

隨著電子計算機的發展,數值模擬的優越性越來越顯著,通過建立合理的細觀力學模型分析復合材料的微觀結構對宏觀力學性能的影響已經成為復合材料設計的重要輔助手段,材料和力學工作者一直致力于研究直接用數值方法進行復合材料設計。黃乾鈺[6]利用Laplace變換與逆變換從基體的黏彈性出發推導了復合材料的黏彈性模型,沈珉等[7]采用細觀力學的Mori-Tanaka法研究非理想界面剛度對復合材料有效彈性模量的影響,姜劍等[8]在細觀力學基礎上利用數值模擬法分析了大變形條件下短纖維增強復合材料的力學性能,李慶[9-10]通過細觀力學有限元方法對炭黑顆粒填充橡膠復合材料的宏觀力學行為進行模擬仿真,重點分析單個圓形和方形炭黑填料粒子模型的變形場和應力場。上述結果表明有限元法對研究單個增強相復合材料的力學性能非常有效。

此外,雷友峰等[11]通過分析復合材料細觀結構代表性體積元的力學響應,基于能量等效原理計算復合材料有效彈性模量,數值計算結果與部分實驗結果有較好的吻合度;高劍虹[12]利用該方法計算了短纖維增強橡膠復合材料的有效彈性模量。Charles L.等[13]比較了用不同方法計算短纖維增強復合材料有效剛度的差異,證明Mori-Tanaka法是預測短纖維復合材料有效彈性模量的理想方法。黃永霞等[14]基于均勻化理論對隨機夾雜分布的應力場計算模擬,得出Voronoi單元有限元法在計算多相復合材料時比普通位移有限元法效率高,并且能夠反映夾雜分布的隨機性。

汪文學等[15]明確指出了復雜應力狀態下材料主應力分量之間的非線性耦合對材料非線性行為的影響,并基于余應變能密度函數提出了新的、包含了非線性耦合效應的應力應變關系式。劉熠等[16]以遠場球對稱應力作用下非線性基體中含單個球形剛性粒子或單個球形孔洞為例,給出了非線性基體中 Eshelby等效夾雜方法的精確解,表明材料線性情況的Eshelby等效夾雜方法不能簡單地推廣應用于非線性基體中含剛性粒子或孔洞的力學分析中。

本文在ansys平臺上利用Monte-Carlo隨機投放法建立了隨機顆粒增強橡膠復合材料的有限元模型,分別用廣義虎克定律(有限元法)、應變能法[11-12]、Mori-Tanaka法[17]、Halpin-Tsai法[18]、修正的Halpin-Tsai[19]法和均質化法[20]六種方法研究顆粒增強橡膠復合材料中填充顆粒的半徑大小、百分含量、彈性模量等參數對復合材料有效力學性能的影響,通過比較各方法的計算結果,分析各個方法的差異,討論其優劣,為準確預測隨機顆粒增強橡膠復合材料的有效力學性能提供可靠方法。

1 有效彈性模量計算的理論基礎

1) 廣義虎克定律(有限元法)

以平面應力狀態為例,考慮線彈性情況,利用廣義虎克定律

{σij}=[D]{εij}

(1)

(2)

可得出:

(3)

(4)

2) 應變能法

雷友峰等[11]依據能量等效原理,利用細觀力學有限元法推導了以平均應變能密度和平均應變為變量的計算復合材料有效彈性模量的計算公式:

(5)

式中:∑Uij為復合材料單元平均應變能密度,ε為復合材料平均應變,V為復合材料單元體的總體積。

3)Mori-Tanaka法與Halpin-Tsai理論

Mori和Tanaka利用Eshelby夾雜理論和平均應力-應變的概念推導出一個計算短纖維增強復合材料有效力學性能的經典方法[17]。

設均質材料在其邊界上受到遠場均勻的應力σ0的作用,其本構關系為:

σ0=L0ε0

(6)

(7)

基體中應力的擾動部分為

(8)

利用Eshelby等效夾雜原理處理,得到夾雜相的應力為

(9)

式中:L1為夾雜相的彈性常數張量;ε*為夾雜的等效本征應變;σ′與ε′為由于單個夾雜的存在而相對于原本基體所引起的擾動應力和應變,采用Eshelby的推導結果有

ε′=Sε*

(10)

式中,S為Eshelby四階張量,根據平均場理論,最終得到復合材料的等效彈性模量為

L=L0(I+vfA)-1

(11)

式中,vf為夾雜相體積比,A為:

A={L0+(L1-L0)[vfI+(1-vf)S}-1(L0-L1)

(12)

由式(11)可以看出,Mori-Tanaka法得到的有效彈性模量除了與夾雜的體積比和彈性模量有關,還通過Eshelby張量體現夾雜之間的相互作用。但該公式無法體現夾雜的排列方式和夾雜尺寸大小的變化,不能很好體現夾雜之間的相互作用。

Halpin-Tsai經典理論源于利用連續介質理論推出的廣義自洽模型,對于體積含量較小的連續纖維增強復合材料有效彈性模量預測具有較高的準確率。該理論導出的通式:

(13)

可以用于包括彈性模量、剪切模量和泊松比在內的所有有效力學性能參數的計算。式中:P表示E11、E22、G12、γ12等不同的力學性能參數,f表示顆粒材料,m表示基質材料,η表示取不同力學性能計算得到的一個參數,ξ是一個與長徑比(l/d) 有關的材料參數,其定義見表1所示,表1同時定義了P代表的幾種典型的力學性能。

從Halpin-Tsai經典公式(13)可以看出,公式只包含了基體和夾雜材料的體積比和各自的材料參數,無法體現夾雜相的幾何形狀及夾雜相之間的相互作用。

表1 Halpin-Tsai經典公式P代表的力學性能

4) 修正的Halpain-Tsai法

劉平等[15]進一步探討 Halpin-Tsai修正混合率計算模型及其原理,并對其經驗擬合參數q進行理論識別,導出了著名的 Halpin-Tsai修正混合率公式,擴展了該公式的應用范圍:

0

(14)

4) 均質化法

多尺度均質化法將宏觀結構分解為無數周期性特征體積單元,引入攝動參變量ξ,利用攝動理論將與宏微觀變量有關的邊界值問題解耦為宏觀邊值問題和細觀邊值問題。在宏觀邊值問題中,組成Ω的非均勻材料被等效為一種均勻材料,而其等效材料性能通過在單胞上求解細觀邊值問題得到,細觀邊值問題的控制方程為均勻化方程,利用多尺度均勻化法能夠得到具有足夠精度的近似結果,同時還能理解不同的微結構材料性質對非均勻材料整體和局部響應的影響。

根據均質化原理,細觀尺度的位移展開量與宏觀位移之間滿足如下關系

(15)

(16)

羅冬梅等[20]通過引進新的解耦特征函數,得到一個齊次微分方程,只要利用周期性邊界條件的基本變形求得每種變形情況下的應變,就可以確定精確的特征函數,得到復合材料的有效彈性模量。

通過Kronecker函數

(17)

定義過渡函數

(18)

得到解耦特征函數

(19)

將式(19)代入方程(16),即可得到求解解耦特征函數的齊次微分方程:

(20)

從而得到均質化有效彈性模量

(21)

二維情況下,矩形截面的周期性位移邊界條件為:

(22)

式中,u1代表細觀坐標中任意的水平和豎直位移。

2 有限元建模

由于單顆粒夾雜代表性單元只適用于研究顆粒分布均勻的復合材料,對于顆粒分布不均勻的復合材料并不適用,所以,為了建立比較接近真實情況的代表性單元,考慮到顆粒之間的相互作用,本文利用Monte-Carlo隨機投放技術,先在指定的投放區域用ANSYS中的命令生成方框模型,調用數組SPH中的數據,使用WPAVE命令、SPH4命令和APDL參數化設計語言中的循環語句生成隨機顆粒模型,然后使用VSBV命令(體相減命令)將立方體區域減去顆粒區域即得到基體區域。最后用VGLUE(粘結體命令)將顆粒區域和基體區域相粘結,完成模型的構建,用來分析復合材料的有效力學性能。圖1為隨機顆粒增強橡膠復合材料的代表性單元,圖2為對應的有限元網格模型。

圖1 二維顆粒生成和投放模型

圖2 有限元網格模型

3 計算結果與討論

橡膠屬于大變形超彈性材料,不是常規的線彈性材料,彈性模量沒有固定值,本文的研究對像是丁腈橡膠,一般彈性模量為2~8 MPa,耐熱性120℃,耐寒性-40℃,耐候性佳,耐磨損,抗變形性好,不抗燃,儲存穩定年份5到10年。

橡膠的本構關系采用Mooney模型,其應變能函數為

W=C1(I1-3)+C2(I2-3)

(23)

其中:C1和C2為橡膠的材料常數,I1和I2分別為橡膠的第一和第二應變不變量,為方便比較,取C1=7E/48,C2=E/48,將橡膠材料的初始彈性模量E=2 MPa代入其中得到:C1=0.292,C2=0.042,泊松比為0.499[7];顆粒的彈性模量為200 MPa,泊松比為0.3,本文主要討論有效彈性模量和有效泊松比的變化。

3.1 顆粒半徑的影響

設隨機投放的顆粒半徑分別為2、2.5、3、4、6 nm,所占面積百分比為15%,以常應變0.4進行加載運算,得到半徑與有效彈性模量之間的關系如圖3所示。

圖3 顆粒半徑對復合材料有效模量的影響

圖4 顆粒半徑對復合材料有效泊松比的影響

從圖3可以看出,(1)顆粒只在小半徑情況對彈性模量影響較大,隨半徑的增大,彈性模量逐漸趨于穩定。(2)計算方法對彈性模量的影響較為顯著,有限元法得到的結果偏大,修正Halpin-Tsai結果是幾種計算方法中最為保守的。與有限元法結果相比,兩者最大值相差36.5%,均質化法結果與應變能法結果較為接近。Mori-Tanaka法和Halpin-Tsai兩個經典理論都難以體現顆粒半徑的影響,圖中沒有顯示這兩種方法的結果,而修正的Halpin-Tsai法通過應力和應變的變化能夠體現顆粒半徑的影響,但結果比其他方法均小,偏于保守。(3)盡管顆粒彈性模量是橡膠基體材料的100倍,顆粒百分比也有15%,但由于橡膠材料的超彈性特性,復合材料的有效彈性模量增長并不明顯。

圖4為半徑變化對泊松比的影響,應變能法預測到的泊松比最大,為0.486,與均質化法得到的結果最接近。修正的Halpin-Tsai法預測的泊松比只有0.428,依然是所有預測結果中最保守的值,與最大值相差10.6%,有限元法得到的結果居中。總的來說,橡膠材料對大半徑顆粒不敏感,Mori-Tanaka法和Halpin-Tsai兩個經典理論不能體現顆粒半徑的影響,其他四種方法預測的結果差別不大。

3.2 顆粒百分比的影響

選取顆粒半徑為4 nm,顆粒彈性模量為200 MPa,顆粒百分比分別為5%、10%、15%、20%、25%、30%、40%,加載應變取0.4。顆粒百分比與有效彈性模量和有效泊松比的關系分別如圖5和圖6所示。

圖5 顆粒百分比對有效彈性模量的影響

圖6 顆粒百分比對有效泊松比的影響

圖5顯示顆粒百分比的增加使復合材料的有效彈性模量緩慢增大,其中修正的Halpin-Tsai法增長幅度最大,顆粒體積比為40%時,最大的有效彈性模量大約比基體材料增大80%,Halpin-Tsai得到的結果最小,由該方法算得的最大有效彈性模量只比基體材料增大51.5%。不同方法得到的最大有效彈性模量相差18.4%。有限元法和應變能法源自最基本的計算公式,其所用的應力-應變直接來自數值模擬,受網格尺寸、邊界條件處理、加載方式的選擇等因素的影響,其結果偏向修正的Halpin-Tsai法;經典Mori-Tanaka法在小體積比時的結果偏小,與Halpin-Tsai法接近,隨著顆粒體積比的增加,越來越與均質化法結果接近,處于幾種方法的中間,均質化法既考慮了顆粒之間的相互影響,也通過統計平均的方法考慮了不同加載方式的影響,結果比較可靠。圖6顯示Mori-Tanaka法預測的泊松比偏小,與最大泊松比相差28.3%,其他幾種方法的結果變化趨勢非常一致,且數值上也相差不大,說明Mori-Tanaka法在泊松比的預測上與其他方法有較大差異。

3.3 顆粒彈性模量的影響

選取顆粒半徑為4 nm,顆粒體積比為15%,加載應變取0.4,改變顆粒的彈性模量分別為200 MPa、400 MPa、600 MPa、800 MPa、1 000 MPa。顆粒彈性模量的變化對復合材料有效彈性模量和有效泊松比的影響如圖7和圖8所示。

圖7 顆粒彈性模量對有效彈性模量的影響

圖8 顆粒彈性模量對有效泊松比的影響

從圖7可以看出每一種方法都顯示復合材料的有效彈性模量在顆粒彈性模量較小時變化較大。顆粒彈性模量超過200 MPa之后,顆粒彈性模量的變化幾乎不再影響復合材料的有效彈性模量,有效彈性模量趨于常數。但不同的方法彈性模量趨近的值不一樣,修正的Halpin-Tsai法得到的值最大,其次是均質化法,有限元法和應變能法結果居中,Mori-Tanaka法得到的結果最小,與最大值相差30%。圖8也顯示顆粒彈性模量小于200 MPa時泊松比才有變化,超過200 MPa 后就趨近常數,Mori-Tanaka法的結果最大,幾乎接近橡膠基體的泊松比,較難體現體積比變化的影響,修正的Halpin-Tsai法則先下降然后逐漸趨于穩定,但其數值是所有結果中最小的。其他幾種方法的結果介于兩者之間,最大值與最小值相差12.2%,有較好的一致性。

4 結論

本文用六種不同的方法計算了隨機顆粒增強橡膠復合材料的有效彈性模量和有效泊松比,比較所得到的結果及其相互之間的差異,得出以下結論:

1) 顆粒增強橡膠復合材料有效彈性模量隨顆粒半徑的增大而逐漸減小并趨于穩定,變化有效泊松比的大小也產生類似的變化趨勢,證明小尺寸顆粒在提高材料的剛度方面作用更大。

2) 顆粒增強橡膠復合材料有效彈性模量隨顆粒百分比增大而增大,說明用小顆粒大體積比的增強相能有效改善材料的剛度;與其他方法相比,Mori- Tanaka法預測的有效泊松比過于偏小。

3) 顆粒增強橡膠復合材料有效彈性模量中顆粒彈性模量的影響只在比較小的情況下較明顯,隨著顆粒彈性模量的增大逐漸趨于穩定,但每種方法所得到的值有一定差異;泊松比與彈性模量的變化趨勢相互吻合。

4) 用于預測有效力學性能的幾種方法中,經典Halpin-Tsai和Mori-Tanaka法在全面體現顆粒特性影響方面有一定的局限性,有限元法和應變能法的結果容易受單元劃分和邊界條件的影響,修正的Halpin-Tsai法綜合考慮了經典理論與有限元結合的優勢,均質化法通過統計平均克服了有限元計算中的累積誤差,因此修正的Halpin-Tsai法和均質化法所得的結果具有更高的可靠性。

[1] 李福強,陳福林.幾種短纖維對三元乙丙橡膠/短纖維復合材料性能的影響[J].廣東橡膠,2010(8):7-10.

[2] 謝尊虎,曾凡偉,肖建斌.硅橡膠性能及其研究進展[J].特種橡膠制品,2011,32(2): 69-72.

[3] 王作齡(編譯).橡膠試驗方法(二十六)—摘自日本《ゴム試驗法》[J].橡塑資源利用,2011(1):36-48.

[4] 齊海波,呂寶華,王文華.顆粒增強復合材料彈性模量的統計分析方法[J].工程力學,2001(S):295-299.

[5] 柏振海,黎文獻,羅兵輝.一種復合材料彈性模量的計算方法[J].中南大學學報,2006,37(3):438-443.

[6] 黃乾鈺.聚合物基復合材料黏彈性性能預測及其應用研究[J].武漢理工大學大學學報,2015(4).

[7] 沈珉,郝培.顆粒增強復合材料非理想界面剛度和有效模量的理論估計[J].復合材料學報,2016(1):189-197.

[8] 姜劍,顧伯勤,張斌.大變形條件下短纖維增強橡膠基復合材料的力學行為[J].南京工業大學學報(自然科學版),2016,38(4):101-104.

[9] 李慶.周期性邊界條件下炭黑增強橡膠基復合材料有效彈性性能數值模擬[J].福州大學學報(自然科學版),2013(1):97-103.

[10]李慶,楊曉翔.顆粒增強橡膠細觀力學性能二維數值模擬[J].應用力學學報,2012,29(5):607-612.

[11]雷友峰,魏德明,高德平.細觀力學有限元法預測復合材料宏觀有效彈性模量[J].燃氣渦輪試驗研究,2003,1(3):11-16.

[12]高劍紅,吳奕錦,陳佳彬.短纖維橡膠復合材料的非線性有限元研究[J].泉州師范學院學報,2014,32(2):31-35.

[13]CHARLES L.,TUCKER III,ERWIN LIANG.Stiffness predictions for unidirectional short-fiber composites:Review and evaluation[J].Composites Science and Technology,1999,59,655-671.

[14]黃永霞,郭然,李偉.顆粒增強復合材料有效模量的Voronoi單元有限元法分析[J].重慶大學學報,2016(5):63-72.

[15]汪文學,高雄善裕.正交各向異性復合材料板的非線性彈性應力應變關系[J].固體力學學報,1991,12(3):353-358.

[16]劉熠,黃筑平,王仁.關于非線性基體中 Eshelby 等效夾雜方法適[J].力學學報,1997,29(4):506-512.

[17]MORI T,TANAKA K.Average stressin matrix and average elastic energy of materials with misfitting inclusions[J].Acta Metallurgica, 1973(21):571-574.

[18]HALPIN J C,KARDOS J L.The Halpin-Tsai equations:A review[J].Polym Eng Sci,1976(16):344-352.

[19]劉平,萬澤青.復合材料等效彈性模量預測方法的改進[J].揚州大學學報,2007,10(1):21-23.

[20]羅冬梅,汪文學,高雄善裕,等.確定均質化法中精確周期性邊界條件的新解法及其在復合材料剛度預測中的應用[J].機械強度,2006,28(4):217-523.

(責任編輯 楊繼森)

Prediction of Effective Mechanical Properties of Rubber Composites Reinforced with Particles

XIE Yue, SU Xiao-ru, FENG Chun-dong, LUO Dong-mei

(Department of Civil Engineering, Foshan University, Foshan 528000, China)

The influence of radius, volume fraction, particle’s elastic modulus and the other parameters on the effective mechanical properties of composite materials is analyzed by several classical calculation methods. The results show that it is beneficial to improve the effective mechanical properties of random particles reinforced rubber composite by using the particles with high volume fraction, small radius and low elastic modulus. The results of the finite element method based on the strain-energy method and the generalized Hooke’s law are greater than those from other methods, which is effected by the mesh division and boundary condition in the finite element method. The results from the traditional Mori-Tanaka method and Halpin-tsai theory are limited for some factors, and they cannot consider the influence of the sizes of particles and their inter-action. The modified Halpin-Tsai formula can improve the results from traditional Halpin-Tsai theory effectively, and the results from multi-scale homogenization methods are always among of the proposed methods, and it is proved that the homogenization method and the modified Halpin-Tsai formula are more reliable.

particle; rubber composites; effective mechanical properties; homogenization method

10.11809/scbgxb2017.07.030

2017-01-17;

2017-02-26

國家自然科學基金項目(10772047/A020206,11172066/A020305);廣東省自然科學基金項目(S2011010004874);佛山科學技術學院大學生創新創業資助項目“佛山市高校和醫院科研基礎平臺”(2016AG100341)

謝悅(1992—),碩士,主要從事土木工程結構的分析與研究。

羅冬梅(1965—),女,博士,教授,主要從事納微米復合材料的力學性能研究。

format:XIE Yue, SU Xiao-ru, FENG Chun-dong, et al.Prediction of Effective Mechanical Properties of Rubber Composites Reinforced with Particles [J].Journal of Ordnance Equipment Engineering,2017(7):142-147.

O343.5;TJ04

A

2096-2304(2017)07-0142-06

本文引用格式:謝悅,宿曉如,馮春冬,等.顆粒增強橡膠復合材料有效力學性能預測分析[J].兵器裝備工程學報,2017(7):142-147.

猜你喜歡
有限元法力學性能復合材料
Pr對20MnSi力學性能的影響
云南化工(2021年11期)2022-01-12 06:06:14
正交各向異性材料裂紋疲勞擴展的擴展有限元法研究
Mn-Si對ZG1Cr11Ni2WMoV鋼力學性能的影響
山東冶金(2019年3期)2019-07-10 00:54:00
民機復合材料的適航鑒定
復合材料無損檢測探討
電子測試(2017年11期)2017-12-15 08:57:13
INCONEL625+X65復合管的焊接組織與力學性能
焊接(2015年9期)2015-07-18 11:03:53
三維有限元法在口腔正畸生物力學研究中發揮的作用
TiO2/ACF復合材料的制備及表征
應用化工(2014年10期)2014-08-16 13:11:29
集成對稱模糊數及有限元法的切削力預測
RGO/C3N4復合材料的制備及可見光催化性能
主站蜘蛛池模板: 99精品免费欧美成人小视频| 国产精品爆乳99久久| 亚洲欧美另类视频| 色偷偷男人的天堂亚洲av| 国产日本一线在线观看免费| 国产特级毛片aaaaaa| 亚洲国产精品成人久久综合影院| 国产一级在线播放| 精品综合久久久久久97| 亚洲欧美国产视频| 精品在线免费播放| 少妇精品在线| 91在线播放免费不卡无毒| 美女高潮全身流白浆福利区| 国产高清在线观看| 又猛又黄又爽无遮挡的视频网站| 一级毛片基地| 久久伊人操| 国产剧情国内精品原创| 一级一级一片免费| 久久永久视频| AV片亚洲国产男人的天堂| m男亚洲一区中文字幕| 四虎影视无码永久免费观看| 麻豆精品国产自产在线| 欧美www在线观看| 欧洲熟妇精品视频| 欧美午夜一区| 成人精品亚洲| 少妇露出福利视频| 欧美精品一区在线看| 伊人激情综合网| 日韩精品高清自在线| 99久久精品无码专区免费| 亚洲av无码人妻| 伊人成人在线| 香蕉视频在线观看www| 国产精品久线在线观看| 福利姬国产精品一区在线| 亚洲欧美人成人让影院| 欧美另类图片视频无弹跳第一页| 亚洲天堂网2014| 国产精品人人做人人爽人人添| 91小视频在线| а∨天堂一区中文字幕| 91偷拍一区| 久青草网站| 日本高清视频在线www色| 亚洲欧洲日本在线| 午夜精品久久久久久久无码软件 | 538国产视频| jizz国产视频| 国产Av无码精品色午夜| 日韩福利视频导航| 欧美不卡在线视频| 久久黄色一级片| 国产成人一区二区| 免费国产高清视频| 亚洲AV一二三区无码AV蜜桃| 国产精品国产主播在线观看| 亚洲自偷自拍另类小说| 无码内射在线| 日韩在线第三页| 在线观看免费黄色网址| 黄色网站不卡无码| 丁香婷婷激情网| 国产地址二永久伊甸园| 免费日韩在线视频| 欧美一区国产| 理论片一区| 国模私拍一区二区| 亚洲欧美综合精品久久成人网| 精品视频在线一区| 国产人人干| 91网红精品在线观看| 人妻精品久久久无码区色视| 日韩激情成人| 91色在线观看| 国产在线精品99一区不卡| 日韩高清欧美| 国产欧美性爱网| 成人亚洲国产|