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

芳綸短纖維增強橡膠復合材料界面力學性能表征及有限元模擬

2022-08-30 09:06:12劉霞楊曉翔高劍虹
福州大學學報(自然科學版) 2022年3期
關鍵詞:復合材料界面有限元

劉霞,楊曉翔,高劍虹

(1.福州大學機械工程及自動化學院,福建 福州 350108;2.泉州師范學院化工與材料學院,福建 泉州 362000)

0 引言

短纖維增強橡膠材料的應用領域廣泛,因其力學性能好、來源易獲取、輕量化等優點逐漸成為傳統材料的替代品[1],材料性能的精確預測為復合材料的安全使用提供了保障.由于細觀層面下的界面性能無法直接測量,建立細觀結構和復合材料宏觀力學性能之間的聯系顯得格外重要[2].

從短纖維增強復合材料的斷裂實驗來看[3],界面并非完全粘結:基體發生橫向斷裂、纖維從基體中拔出.在纖維與基體之間引入損傷模型能模擬纖維與基體之間的脫粘.界面損傷模型在細觀層面預測材料的損傷過程和力學性能得到了廣泛的應用[4-7].沈珉等[8]基于雙線性內聚力模型建立隨機植物短纖維復合材料二維代表體積單元,并研究了界面性能對有效彈性模量和拉伸行為的影響.段宇星等[9]基于剪滯理論建立短纖維增強橡膠材料的理論和數值模型,并對短纖維增強阻尼材料的彈性模型進行計算與擬合.

現階段,對纖維增強橡膠復合材料的性能探索以建立二維平面模型較主流,有限元技術的發展使建立復雜的細觀結構力學模型成為可能.Vaughan等[10-11]提出一種二維單胞復合材料模型建立的新算法,代表體積單元(representative volume element,RVE)單胞模型能保證纖維之間的分布距離跟實際細觀吻合.Qi等[12]基于Python-ABAQUS提出一種改進的隨機序列吸附算法(random sequence adsorption algorithm,RSA)來生成纖維,研究了具有一定長徑比的空間隨機分布短圓柱纖維增強復合材料的有效彈性性能.Ge等[13]提出一種基于迭代改進的生成高纖維體積分數(高達70%)隨機分布連續纖維增強復合材料細觀結構的新方法,克服了傳統隨機序列吸附法的干擾限制.該算法計算效率高,纖維生成隨機性好,但其RVE局限于二維平面且沒有考慮界面損傷.

本研究采用數值仿真/實驗相結合的方法,研究芳綸短纖維增強橡膠復合材料拉伸力學性能.基于Python-ABAQUS二次開發,建立含損傷界面的隨機纖維分布二維代表體積單元,采用cohesive單元表征界面性能,通過編寫的Python腳本程序對模型自動施加周期性邊界條件.最后通過ABAQUS模擬仿真RVE對拉伸載荷的力學響應,分析不同界面參數以及纖維含量對芳綸短纖維增強橡膠復合材料(aramid fiber-reinforced rubber composite,AFRC)軸向拉伸力學性能的影響.

1 實驗設計及過程

由于芳綸纖維高強度和高模量等優異的力學性能,在橡膠基體中少量加入即可提高制品的整體性能.為考察復合材料的宏觀力學性能,分別采用Kevlar對位芳綸短切纖維制備了4種配方的AFRC,見表1.表中配方代號F表示無添加纖維的橡膠基體,K表示添加Kevlar對位芳綸纖維.文中研究的Kevlar纖維長度為l=1.5 mm,試件中芳綸纖維的質量份數從1份增加到7份,而橡膠基體、炭黑和其余配合劑的摻加量均相同.實驗煉膠分二段混煉,分別在密煉機和開煉機上進行,各段混煉膠料停放時間不少于8 h,在平板硫化機上進行硫化,硫化溫度為143 ℃,時間為45 min.

表1 芳綸纖維增強橡膠復合材料配方Tab.1 Formula of aramid fiber-reinforced rubber composite

由于芳綸纖維的投放是無序且不定向的,在膠料中的分布并不完全呈現沿軸向分布[14].拉伸試樣所用的裁刀和裁片機符合國標GB/T2941[15]的要求,所制備的橡膠尺寸參考國家標準《硫化橡膠或熱塑性橡膠拉伸應力應變性能的測定》(GB/T528-2009[16])中“I型”啞鈴狀試樣的要求(如圖1所示).采用日本島津公司研發的AG-X plus島津電子萬能試驗機,在常溫(23 ℃)對啞鈴形試件進行單向拉伸.

圖1 I型啞鈴形狀試件形狀和尺寸(單位:mm)Fig.1 Shape and size of I-shaped dumbbell specimen (unit:mm)

采用1 mm·min-1的拉伸速率,對試件進行準靜態單軸拉伸實驗,結果見圖2.圖2中列出了纖維體積分數為0.486%、1.458%、2.430%和3.402%試件的實驗拉伸應力-應變曲線和將纖維/基體界面看作是理想粘結界面的有限元仿真應力-應變曲線.從圖2中可以看出,纖維體積含量對AFRC的拉伸應力-應變曲線有一定的影響,體現在纖維體積含量越高,同應變條件下AFRC的真實應力越高.K1試件的應力-應變曲線位置與K2試件基本重合,而K4的應力-應變曲線的位置較K2有明顯的變化,驗證了在一定的纖維含量范圍內,纖維含量的增加對AFRC的增強效果明顯,而對于小含量的纖維對AFRC的增強作用變化較微弱.

圖2 K系列試件受單軸拉伸載荷時的應力-應變曲線Fig.2 Stress-strain curve of K series specimen under uniaxial tensile load

2 有限元模型的建立

2.1 RVE單胞模型建立

復合材料在細觀層面的力學響應特性可反映其在宏觀層面的有效性能[17].復合材料在細觀層面并不是完全均勻的,但是纖維在復合材料的分布具有一定規律性,在統計學上稱其為周期性隨機分布,將包含足夠多的信息又能獨立代表復合材料的RVE擴展排布即可構成復合材料的整體結構.RVE是整體結構中具有代表的一部分,有些纖維部分在建模時難以避免會超出RVE的邊界,則需要將超出部分切割后,將其平移至對應邊的RVE內部[18].

采用Python語言編程建立二維短纖維隨機分布代表體積單元(RVE)模型.考慮芳綸纖維的體積含量、纖維的隨機位置分布以及纖維的大長徑比特性等細觀結構信息.程序編寫分為如下步驟:

1) 建立力學簡化模型,創建二維基體,其大小與RVE大小一致,并在該基體中創建第一根纖維.

2) 兩次調用Random函數分別完成對第一根纖維的平面旋轉和平面平移得到第二根纖維.

3) 利用判斷語句對投放纖維體積分數的判斷,未滿足設定纖維體積分數的投放時,則繼續投放纖維,直到達到設定參數.

4) 輸出模型.為了保證周期性邊界條件的順利自動實施,需要在RVE的對邊保證網格的一致性,這步可通過網格的復制得到保證.通過Python編寫的程序可獲得互不干涉的隨機分布纖維,建立隨機分布纖維增強橡膠復合材料的代表體積單元,模擬研究AFRC的單軸拉伸力學性能.建立了纖維體積分數為0.486%、1.458%、2.430%和3.402%的二維RVE模型,其中,纖維含量為3.402% AFRC的RVE有限元模型如圖3所示.針對本研究所建立的有限元模型,橡膠基體和芳綸纖維采用線性減縮積分單元(CPS4R);為保有限元模型能容納足夠多的纖維,模型尺寸取6 mm × 6 mm.

圖3 纖維含量為3.402% AFRC的有限元二維數值模型Fig.3 Two-dimensional finite element model of AFRC with 3.402% fiber

2.2 材料參數及界面表征

將隨機芳綸纖維增強橡膠復合材料的纖維/基體界面看作是非理想界面.雙線性內聚力模型能很好地描述應力和位移的非線性關系(如圖4所示),可用來表征AFRC的非理想界面.從圖2可以看出,采用理想界面(界面不發生脫粘)模擬材料的力學性能與真實的應力應變曲線有一定的偏差.利用ABAQUS軟件建立的RVE模型中,采用零厚度二維四節點cohesive單元(COH2D4)表征纖維與基體之間的非理想界面力學行為,界面單元的本構模型遵循雙線性內聚力模型,認為纖維與基體之間的界面性能皆相同,即界面單元的本構模型參數皆相同,有σmax=τmax,δ1=δ1n=δ1s,δd=δdn=δds.在模擬過程中基體與纖維不發生斷裂,模型僅在纖維與基體之間的界面發生破壞.

圖4 界面的雙線性內聚力模型Fig.4 Elastic-softening cohesive zone model for the interface

界面單元處于線彈性階段(0≤δ≤δ1),應力隨著分離位移增大而增大.當界面應力達到損傷初始準則后,界面損傷開始.采用的內聚力單元的初始損傷判據如下:

(1)

式中:σn和τt分別是界面的法向應力和切向應力;[σn]和[τt]分別是界面的法向脫粘強度和切向脫粘強度.

界面進入損傷演化階段(δ1≤δ≤δd),應力隨著分離位移增大而減小,直至為零,此時界面完全失效.界面損傷參數(scalar stiffness degradation,SDEG)可描述為界面損傷程度,損傷參數SDEG可表示為:

(2)

式中:δ為損傷過程中界面的最大分離位移.

芳綸纖維可視為各向同性彈性材料,其彈性模量為70 GPa,泊松比為0.3.纖維在基體材料的分布位置是隨機且不定向的.

基體是均勻各向同性的超彈性材料.根據唯象理論,橡膠類材料的應變能函數W是右柯西-格林變形張量C=FTF的函數,即:

W=W(C)=W(I1,I2,I3)

(3)

式中:F=?x/?X為變形梯度;x和X分別為歐拉構型與拉格朗日構型中的位移矢量;Ii為應變不變量.

采用經典的超彈性本構模型Ogden N3[19]描述橡膠基體的應力-應變行為,Ogden N3模型可表示為:

(4)

利用ABAQUS有限元軟件中的材料參數回歸模塊,得到表1中F0無纖維試件的材料參數如下:μ1=8.116 MPa,α1=1.975,μ2=5.725 MPa,α2=2.388,μ3=5.467 MPa,α3=4.495.由于纖維的面內分布特點,假設二維數值模型為平面應力問題.

2.3 周期性邊界條件的建立

對細觀有限元模型合理施加邊界條件有利于獲得準確的細觀力學響應.周期性邊界條件的實施必須滿足對應邊界面的應力連續和位移連續[20].AFRC的二維RVE幾何表述見圖5,對RVE模型施加周期性邊界條件,在相對平行的兩個對面上以及在各個頂點滿足:

圖5 AFRC的二維RVE幾何表述Fig.5 Geometrical frame of two-dimentional RVE for AFRC

ul2=ul1+(uA1-uA0);ul4=ul3+(uA2-uA0)

(5)

式中:u表示邊界上節點的位移向量;下標l1,l2,l3,l4代表4條邊;A0,A1,A2,A3代表4個角點.

基于周期性邊界條件理論與ABAQUS二次開發程序的編寫規則,編寫對有限元模型自動添加周期性邊界的腳本,并通過“equation”實現式(5)在ABAQUS的用戶界面的自動運行.約束A0點在x和y的位移,約束A1點在y的位移,約束A2點在x的位移,并通過在A2點施加y向位移載荷以實現對代表體積單元的單向拉伸模擬.

3 有限元分析

通過RVE模型可計算不同纖維體積含量的AFRC應力應變曲線,其中AFRC的真實應力由加載方向的約束反力除以變形后的表面積計算得到,真實應變ε的計算公式為:

(6)

式中:u為加載施加的位移值;L為RVE的邊長.

通過式(7)可以計算得到AFRC在單向拉伸時的應力應變曲線.通過取曲線在0.002應變下的正切模量來近似AFRC的有效彈性模量E[21].

3.1 網格密度對有效彈性模量的影響

對于體積分數為4.302%的RVE幾何模型建立了不同的網格密度,圖6給出了AFRC的有效彈性模量與單元網格密度的關系圖.從圖6中的試算結果可以看出,隨著網格密度的增加,有效彈性模量的值逐漸趨于平緩,這說明在一定合適范圍內的網格密度條件下,網格密度對仿真的結果沒有影響.考慮計算機的計算時長和計算收斂性,網格密度不宜過細.在網格密度大于132 × 132時材料的有效彈性模量時預測偏差小于5%,所以本研究在建立RVE模型時,其網格密度取值為132 × 132,基體單元對應的尺寸為0.055 mm.對于芳綸纖維這種大長徑比且直徑非常微小的材料,在仿真時需保證結果的收斂性,需要在纖維與基體連接處進行網格的細化來保證計算結果的收斂性.

圖6 網格密度對AFRC宏觀有效彈性模量的影響Fig.6 Effect of mesh density on the macroscopic effective modulus

3.2 內聚力模型參數對AFRC應力應變曲線的影響

圖7中給出不同界面剛度K、不同初始裂紋產生位移δ1、不同完全破壞位移δd下內聚力模型的5組界面模型數據.分別應用到數值模型中進行有限元仿真,得到了纖維體積含量為3.402%的AFRC的界面參數對應力應變曲線的影響,如圖8所示.

圖7 內聚力模型參數Fig.7 Parameters of cohesive zone model

圖8 界面參數對有效彈性模量的影響Fig.8 Influence of interface parameters on effective modulus

對比S1和S2,在K和δd取值相同、δ1取值不同的情況下,可以看出界面的斷裂能(即曲線與橫軸所圍的面積)對應力應變曲線有很大的影響,著重體現在大斷裂能導致曲線在初始彈性階段非線性特征不明顯,表現為需要大應力才能使曲線進入非線性階段.S1的δ1較大,對應拉伸曲線的最大應力也越大.

對比S3和S4,在斷裂能相同的情況下,S3的界面剛度比S4的界面剛度要大,在應力應變曲線中,S3曲線在初始階段的斜率未與S4重合,界面剛度小則對應的應力應變曲線的斜率小,反映材料在彈性階段的有效彈性模量小,可以說明界面剛度對復合材料的有效彈性模量的估計有一定的影響.

對比S2和S3,保證其他參數不變,只改變界面的δd,在完全失效位移增大一倍的情況下,曲線在非線性階段兩者沒有明顯的分離.特別地,S2和S3有相同的界面剛度值,兩者對應的應力應變曲線在線性階段是重合的,可以推斷界面剛度相同的復合材料其有效彈性模量相同.

綜上所述,以及對比S1、S3和S5,推斷出內聚力模型的界面剛度K只影響曲線的線性階段,對于材料在非線性階段的影響取決于δ1、δd的值.

3.3 界面剛度對AFRC彈性模量的影響

為了研究界面性質對AFRC軸向拉伸性能的影響,模擬相應體積含量的AFRC的拉伸力學行為,選取20~1×106MPa·mm-1之間共15個不同界面剛度值,計算每個界面剛度值所對應的有效彈性模量.圖9給出4種纖維體積分數對應的不同界面剛度下的AFRC軸向拉伸的界面剛度-彈性模量曲線,可以看出,界面剛度對AFRC的有效彈性模量有很大的影響.選取合適的界面剛度對預測AFRC的力學性能有很大的影響.同時,界面剛度的增大使得AFRC的有效彈性模量增大.

圖9 AFRC的剛度-彈性模量曲線Fig.9 E-K curves of AFRC

在同一體積含量下AFRC的有效彈性模量隨著界面剛度的提高,AFRC的軸向拉伸模量有所提高,說明界面剛度對AFRC軸向拉伸有效性能有一定的影響.界面剛度值較大時,可看作是界面的理想粘結情況下,界面剛度值趨近于無窮大,曲線呈現水平穩定值.對于曲線在小界面剛度下,曲線呈現非線性特性,且呈現隨著界面剛度的增大,材料的有效彈性模量也增加.由此可以推斷,界面剛度和宏觀有效性能之間有單調遞增的關系.

3.4 纖維體積含量對AFRC彈性模量的影響

從圖9中可以看出,不同體積分數的曲線皆在一點聚集后分散.將聚集的一點稱為臨界交點,把臨界交點對應AFRC的有效彈性模量稱為臨界有效彈性模量.在本研究所制備的AFRC的臨界界面剛度,按照圖9中曲線擬合計算為50 MPa·mm-1.不同的復合材料的加工工藝所對應的臨界界面剛度不同,可以推斷對纖維與基體之間的界面進行處理可提高材料的有效性能.

在臨界交點的弱界面剛度區域,即非理想界面的剛度小于臨界交點值,隨著纖維體積分數的增加有效彈性模量減小,即小含量纖維對AFRC的增強作用較小.在強界面剛度區域,即非理想界面的剛度大于臨界交點值,纖維體積含量的增大會增加材料的有效彈性模量.但是界面的增加作用不是無限增加的,當界面剛度超過一定值,材料的有效彈性模量會趨于水平線,這種情況下AFRC的界面可看做是理想粘結.根據圖2模擬的理想粘結情況下的應力應變曲線,在仿真中將界面處假設成理想粘結對預測的有效彈性模量以及模擬曲線的應力應變曲線有一定的偏差.

對于纖維含量為0.486%的AFRC,界面剛度對有效彈性模量的影響小,強界面剛度區域對材料的有效彈性模量的增強作用較弱界面區域僅增加27.12%.芳綸的大長徑比特性使得纖維在軸向受力能傳遞應力,但是由于芳綸的直徑小,使得小含量的纖維比大含量的纖維所承受的受力面積少.界面剛度特性反映了一定的界面粘結程度,小含量纖維相比大含量纖維在界面處承受應力傳遞的作用范圍少.圖9中纖維含量為1.458%的AFRC有效彈性模量在界面剛度范圍內增加了86.63%.同樣地,纖維含量為2.430%和3.402%的AFRC有效彈性模量在界面剛度范圍分別增加了192.07%和385.35%.特別地,臨界交點的存在可以說明在進行復合材料的加工工藝中要規避弱界面剛度區域.通過對纖維表面的化學處理[22]或物理粘附[23]改善纖維與基體之間的浸潤性;或使用偶聯劑[24]可增加宏觀復合材料的有效彈性模量.

3.5 界面參數的確定

由實驗給出的纖維含量為0.486%、1.458%、2.430%和3.402%四種不同體積分數下的材料的應力應變曲線,可以得到其宏觀有效彈性模量,在剛度-彈性模量曲線上可找到對應界面剛度的數值,實驗測得的有效彈性模量在圖9用空心點表示.AFRC的有效彈性模量值唯一確定該材料的界面剛度值.

第3.2小節分析了界面參數對材料應力應變曲線的影響,說明界面剛度只影響應力應變曲線在彈性階段的斜率,通過求解實驗應力應變曲線中的有效彈性模量與圖9的剛度-彈性模量曲線比對,可確定在纖維含量、纖維尺寸、纖維隨機分布狀態都跟有限元模型相同情況下的AFRC的界面剛度的具體值.再調試ABAQUS中δ1的值,得到一組應力應變的仿真曲線,與實驗的應力應變曲線進行對比,找到形狀吻合實驗曲線的初始裂紋位移δ1;再調試δd的值,直到找到吻合實驗應力應變曲線在非線性階段的最終裂紋脫粘位移δd.

擬合曲線在彈性階段很好地貼近實驗曲線后,通過模擬應力應變曲線可以在非線彈性階段確定δ1和δd的值.S5組的結果表明在確定K后,δ1和δd對曲線的影響較小.可通過擬合曲線在線性階段的最大值與實驗曲線的契合程度確定δ1的值.S2和S3組指出擬合曲線與實驗曲線在非線性階段的形狀擬合有一定的關系.可通過擬合曲線與實驗曲線的貼合程度確定δd的值.同樣地,采用本文相同的方法可確定纖維體積含量為0.486%、1.458%、2.430%的AFRC界面CZM模型的參數,其具體參數如表2.通過表2預估出的具體的界面剛度K、初始裂紋產生位移δ1和最終失效位移δd,可模擬在同等纖維含量下、纖維分布隨機、纖維尺寸以及相同加工工藝下的AFRC的其他力學性能.

表2 材料界面的內聚力模型參數確定Tab.2 Cohesive zone model parameter determination of material interface

圖10給出4種不同纖維含量通過內聚力模型表征非理想界面有限元仿真和實驗得到的應力應變曲線.在同一應變條件下,圖10中仿真與實驗曲線的應力值偏差最大不超過5%,可以認為,利用ABAQUS軟件仿真得到的4種不同纖維含量的材料的應力應變曲線能很好地跟實驗曲線擬合.

圖10 應力應變曲線仿真與實驗結果對比Fig.10 Comparison of stress-strain curve values and experimental results

4 結語

1) 建立了隨機纖維分布的二維代表體積單元有限元模型,在滿足纖維隨機分布的同時還能保證各纖維之間互不干涉,添加周期性邊界條件的數值模型能更符合真實材料.

2) 通過對AFRC的非理想界面參數的研究表明,有效彈性模量與界面剛度有單調遞增的關系,即隨著界面剛度的提高,AFRC的有效彈性模量增加.在一定纖維體積含量范圍內,纖維體積含量對材料有效彈性模量的增加是有條件的:在弱界面剛度區域,纖維體積含量的增加會導致AFRC有效彈性模量的降低;在強界面剛度區域,纖維體積含量的增加使得AFRC有效彈性模量增加.

3) 通過預估的界面剛度K,初始裂紋位移δ1,最終失效位移δd等界面參數模擬得到的應力應變曲線,與纖維含量0.486%、1.458%、2.430%和3.402% AFRC的實驗應力應變曲線吻合.

猜你喜歡
復合材料界面有限元
國企黨委前置研究的“四個界面”
當代陜西(2020年13期)2020-08-24 08:22:02
民機復合材料的適航鑒定
基于FANUC PICTURE的虛擬軸坐標顯示界面開發方法研究
復合材料無損檢測探討
電子測試(2017年11期)2017-12-15 08:57:13
人機交互界面發展趨勢研究
手機界面中圖形符號的發展趨向
新聞傳播(2015年11期)2015-07-18 11:15:04
磨削淬硬殘余應力的有限元分析
TiO2/ACF復合材料的制備及表征
應用化工(2014年10期)2014-08-16 13:11:29
基于SolidWorks的吸嘴支撐臂有限元分析
RGO/C3N4復合材料的制備及可見光催化性能
主站蜘蛛池模板: 91毛片网| 一本大道视频精品人妻| 国产大片喷水在线在线视频| 五月六月伊人狠狠丁香网| 久久免费视频播放| 国产高清毛片| 日韩国产 在线| 日本成人在线不卡视频| 国产手机在线ΑⅤ片无码观看| 国产91透明丝袜美腿在线| 久青草网站| 国产无人区一区二区三区| 99精品欧美一区| 国内精品久久九九国产精品 | 国产精品一区二区无码免费看片| 亚洲色婷婷一区二区| 亚洲天堂视频网| 91精品啪在线观看国产91九色| 国产内射一区亚洲| 亚洲成人在线网| 亚洲综合久久成人AV| 亚洲不卡av中文在线| 99久久精品久久久久久婷婷| 香蕉色综合| 97亚洲色综久久精品| 一区二区无码在线视频| 国产H片无码不卡在线视频| 国产精品视频a| 国产中文一区二区苍井空| 国产乱人乱偷精品视频a人人澡| 不卡午夜视频| 国产91高跟丝袜| 亚洲国产午夜精华无码福利| 亚洲国产欧美自拍| 国产在线拍偷自揄拍精品| 亚洲成人精品久久| 国产日本一线在线观看免费| 米奇精品一区二区三区| 国产va在线| 国产chinese男男gay视频网| 国产黑人在线| 亚洲免费黄色网| 亚洲精品亚洲人成在线| 在线播放国产一区| 91久久夜色精品国产网站| 久久国产精品波多野结衣| 国产香蕉97碰碰视频VA碰碰看| igao国产精品| 国产美女精品在线| 亚洲欧美日韩中文字幕一区二区三区| 亚洲男女在线| 国产情侣一区| JIZZ亚洲国产| 亚洲欧美色中文字幕| 欧美成人手机在线观看网址| aa级毛片毛片免费观看久| 午夜精品影院| 国产欧美日韩一区二区视频在线| 国产在线精品99一区不卡| 日韩毛片在线播放| 婷婷伊人久久| 国产乱子伦一区二区=| 国产粉嫩粉嫩的18在线播放91| 国产精女同一区二区三区久| 日韩一级毛一欧美一国产| 欧美成人h精品网站| 国产精品无码AⅤ在线观看播放| 亚洲一级毛片免费看| 亚洲精品综合一二三区在线| 国产精品视频免费网站| 国产精品永久在线| 亚洲开心婷婷中文字幕| 五月婷婷亚洲综合| 欧亚日韩Av| 欧美综合一区二区三区| 天天综合网色中文字幕| 欧美激情伊人| 天天综合网色中文字幕| 大学生久久香蕉国产线观看 | 成人国产小视频| 亚洲精品桃花岛av在线| 999国产精品|