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

顆粒增強(qiáng)金屬基復(fù)合材料細(xì)觀有限元建模方法的對(duì)比

2022-08-04 07:07:20劉駿華黃興民
機(jī)械工程材料 2022年4期
關(guān)鍵詞:復(fù)合材料有限元模型

劉駿華,張 娟,張 晨,黃興民

(西南交通大學(xué)1.力學(xué)與航空航天學(xué)院,應(yīng)用力學(xué)與結(jié)構(gòu)安全四川省重點(diǎn)實(shí)驗(yàn)室;2.材料科學(xué)與工程學(xué)院,成都,610031)

0 引 言

顆粒增強(qiáng)金屬基復(fù)合材料具有較高的強(qiáng)度、彈性模量、硬度以及良好的塑性,并且生產(chǎn)成本低廉,應(yīng)用前景廣泛。顆粒增強(qiáng)金屬基復(fù)合材料在承受載荷時(shí)的變形與失效等行為是非常復(fù)雜的,與基體和顆粒兩相材料的性能以及二者之間的多種相互作用密切相關(guān)。在研究微觀及細(xì)觀結(jié)構(gòu)對(duì)復(fù)合材料力學(xué)性能的影響中,采用有限元方法進(jìn)行數(shù)值模擬是一種非常有效的手段。想要進(jìn)行有效的模擬,必須先要建立合理的細(xì)觀有限元模型。目前,顆粒增強(qiáng)復(fù)合材料的細(xì)觀有限元模型主要有體胞模型和真實(shí)結(jié)構(gòu)模型2種。體胞模型的構(gòu)建基于復(fù)合材料的結(jié)構(gòu)特征,可分為單顆粒模型及多顆粒模型[1-2]。單顆粒模型由于其結(jié)構(gòu)簡(jiǎn)單,常用其來(lái)揭示基體、顆粒及界面之間的作用機(jī)理,但無(wú)法探究顆粒之間的相互關(guān)系[3]。對(duì)于多顆粒模型,可以假設(shè)顆粒是均勻分布的,也可以使用隨機(jī)序列吸附法生成隨機(jī)分布多顆粒模型[4]。隨機(jī)分布多顆粒模型相比于單顆粒模型雖然考慮了顆粒間的相互作用,并且一定程度上考慮了材料的不均勻性,但與復(fù)合材料的真實(shí)微觀結(jié)構(gòu)還有一定差異[5]。研究[6]表明,顆粒形貌越接近實(shí)際情況,模擬結(jié)果與試驗(yàn)結(jié)果的誤差越小,因此有學(xué)者發(fā)展了基于材料真實(shí)顯微圖像的建模方式,即真實(shí)結(jié)構(gòu)模型。真實(shí)結(jié)構(gòu)模型一般使用連續(xù)切片法建立,其原理是將試樣分成很多層,取得每一層的平面結(jié)構(gòu)圖像并輸入計(jì)算機(jī)中,然后利用三維重構(gòu)技術(shù),還原出試樣的三維真實(shí)結(jié)構(gòu)[7]。由于三維真實(shí)結(jié)構(gòu)模型的建立需要復(fù)雜的前期工作并且伴隨著龐大的計(jì)算量,因此有不少學(xué)者建立了二維真實(shí)結(jié)構(gòu)模型進(jìn)行模擬,得到的模擬結(jié)果與試驗(yàn)結(jié)果吻合得較好[8-11]。

體胞模型和真實(shí)結(jié)構(gòu)模型各有特點(diǎn),在模擬顆粒增強(qiáng)復(fù)合材料的力學(xué)行為時(shí)得到的結(jié)果也有差異,為了探究這2種模型的異同及特點(diǎn),作者以原位生成法制備的TiB2顆粒增強(qiáng)鐵基復(fù)合材料為研究對(duì)象,通過(guò)納米壓痕試驗(yàn)及有限元反演分析確定了鐵素體基體的冪硬化模型參數(shù),建立了二維的顆粒隨機(jī)分布的體胞模型和真實(shí)結(jié)構(gòu)模型,模擬材料在單軸拉伸時(shí)的力學(xué)行為,研究了2種模型有限元模擬結(jié)果的異同,并進(jìn)一步探究了邊界條件對(duì)模擬結(jié)果的影響,以期在一定程度上為顆粒增強(qiáng)鐵基復(fù)合材料微觀結(jié)構(gòu)的優(yōu)化設(shè)計(jì)提供參考。

1 有限元模型的建立

1.1 材料參數(shù)的確定

在有限元建模的過(guò)程中,需要輸入復(fù)合材料各相的材料參數(shù)。在采用原位生成的方法制備的TiB2顆粒增強(qiáng)鐵基復(fù)合材料中,鐵素體晶粒尺寸與純鐵中的鐵素體晶粒尺寸不同,而且晶粒中會(huì)有位錯(cuò)產(chǎn)生[12],這導(dǎo)致復(fù)合材料中的鐵素體與純鐵中的鐵素體在未受載時(shí)的初始狀態(tài)可能會(huì)有較大不同。若采用純鐵進(jìn)行單軸拉伸試驗(yàn)來(lái)獲得復(fù)合材料鐵素體基體的參數(shù)不合理,因此采用納米壓痕和有限元反演分析的方式獲得鐵素體基體的材料屬性。將TiB2顆粒視為彈性體,彈性模量為510 GPa,泊松比為0.25。對(duì)TiB2顆粒增強(qiáng)鐵基復(fù)合材料進(jìn)行納米壓痕試驗(yàn),測(cè)得基體的彈性模量為214 GPa,硬度為4.01 GPa,然后進(jìn)一步采用DAO等[13]推導(dǎo)的量綱函數(shù)及有限元反演分析,確定基體的單軸冪率硬化應(yīng)力-應(yīng)變響應(yīng)曲線(xiàn)。假設(shè)基體的應(yīng)力σ、應(yīng)變?chǔ)艥M(mǎn)足以下關(guān)系:

(2)

式中:E為基體的彈性模量;σy,εy分別為名義屈服強(qiáng)度與對(duì)應(yīng)的名義屈服應(yīng)變;εp為總應(yīng)變減去名義屈服應(yīng)變;n為冪硬化指數(shù);σr、εr分別為特征應(yīng)力與特征應(yīng)變。

將納米壓痕試驗(yàn)中采用的Berkvich壓頭等效成140.6°的圓錐壓頭,建立軸對(duì)稱(chēng)納米壓痕有限元模型[14]。參考ANTUNES等[15]和LEE等[16]的研究方法進(jìn)行有限元反演,得到鐵素體基體的單軸真應(yīng)力-真應(yīng)變曲線(xiàn)如圖1所示,鐵素體基體的冪硬化指數(shù)n為0.34,屈服強(qiáng)度σy為242.91 MPa,且由圖2可以看出,采用有限元反演分析方法得到壓痕試驗(yàn)過(guò)程中的力-位移曲線(xiàn)與試驗(yàn)結(jié)果相吻合,驗(yàn)證了有限元反演分析方法的準(zhǔn)確性。

圖1 采用有限元反演分析得到鐵素體基體的真應(yīng)力-真應(yīng)變曲線(xiàn)Fig.1 True stress-true strain curve of ferrite matrix by finite element inverse analysis

圖2 有限元反演分析方法得到壓痕試驗(yàn)過(guò)程中的力-位移曲線(xiàn)與試驗(yàn)結(jié)果的對(duì)比Fig.2 Comparison of force-displacement curves during indentation test obtained by finite element inverse analysis with experimental results

1.2 有限元模型的建立

對(duì)體積分?jǐn)?shù)20% TiB2顆粒增強(qiáng)鐵基復(fù)合材料建立二維代表性體積元模型。研究[11]表明,當(dāng)所選取的代表性體積元尺寸大于某一臨界值時(shí),其所代表的微結(jié)構(gòu)能夠反映材料微觀組織的基本信息。所選取的真實(shí)結(jié)構(gòu)模型尺寸為80 μm×80 μm。在復(fù)合材料的掃描電鏡(SEM)形貌中選擇尺寸為80 μm×80 μm合適區(qū)域,經(jīng)二值化處理轉(zhuǎn)化為矢量文件,然后導(dǎo)入ABAQUS軟件中進(jìn)行建模,實(shí)現(xiàn)基于微觀組織的有限元模型的建立,如圖3所示。測(cè)得所選區(qū)域真實(shí)顆粒的平均尺寸為3 μm。

圖3 真實(shí)結(jié)構(gòu)模型建立的3個(gè)階段Fig.3 Three stages of establishing real structure model: (a) SEM micromorphology; (b) results by binarization and (c) geometric model

體胞模型代表性體積元的建立則利用ABAQUS軟件進(jìn)行二次開(kāi)發(fā),采用Python語(yǔ)言編寫(xiě)程序,運(yùn)行程序建立同等顆粒尺寸、體積分?jǐn)?shù)的隨機(jī)分布圓形顆粒及不規(guī)則顆粒的體胞模型,如圖4所示。

圖4 隨機(jī)分布圓形顆粒及不規(guī)則顆粒的體胞模型Fig.4 Body cell models of circular particles with random distribution (a) and irregular particles with random distribution (b)

在ABAQUS軟件中進(jìn)行有限元模擬,將TiB2顆粒設(shè)為彈性體,采用CPS4R單元?jiǎng)澐志W(wǎng)格,平均單元尺寸為0.4 μm。代表性體積元模型一般可施加3種邊界條件[17],即同質(zhì)化邊界條件、對(duì)稱(chēng)性邊界條件和周期性邊界條件,如圖5所示。

圖5 建模可選擇的3種邊界條件Fig.5 Three alternative boundary conditions for modeling: (a) homogenization boundary condition; (b) symmetric boundary condition and (c) periodic boundary condition

1.3 應(yīng)力-應(yīng)變曲線(xiàn)的模擬及試驗(yàn)驗(yàn)證

(4)

式中:A,σij,εij分別為單元面積、應(yīng)力、應(yīng)變。

根據(jù)上述公式編寫(xiě)Python腳本得到采用一階均勻化方法計(jì)算的真應(yīng)力-真應(yīng)變曲線(xiàn)。等效宏觀方法將代表性體積元等效成宏觀模型,在對(duì)模型進(jìn)行單向拉伸加載時(shí),直接輸出拉伸方向的作用力和位移,然后除以模型的橫截面積和長(zhǎng)度,得到工程應(yīng)力-應(yīng)變曲線(xiàn)再轉(zhuǎn)化為真應(yīng)力-真應(yīng)變曲線(xiàn)。

按照GB/T 228—2010,在體積分?jǐn)?shù)20% TiB2顆粒增強(qiáng)鐵基復(fù)合材料上截取如圖6所示的拉伸試樣,用 MTS CMT5105 型電子萬(wàn)能試驗(yàn)機(jī)進(jìn)行室溫拉伸試驗(yàn),拉伸速度為 1 mm·min-1。

圖6 拉伸試樣的尺寸Fig.6 Size of tensile specimen

將用施加對(duì)稱(chēng)性邊界條件的真實(shí)結(jié)構(gòu)模型和體胞模型模擬出的結(jié)果與單軸拉伸試驗(yàn)結(jié)果進(jìn)行對(duì)比,結(jié)果如圖7所示。由圖7可以看出,采用2種方法得到不同模型的真應(yīng)力-真應(yīng)變曲線(xiàn)與試驗(yàn)結(jié)果較吻合,但采用等效宏觀方法計(jì)算得到的真應(yīng)力-真應(yīng)變曲線(xiàn)與試驗(yàn)結(jié)果更加吻合,相對(duì)誤差小于5%,因此后面采用等效宏觀方法計(jì)算真應(yīng)力-真應(yīng)變曲線(xiàn),同時(shí)也證明了建立的真實(shí)結(jié)構(gòu)模型與體胞模型的合理性。

圖7 對(duì)稱(chēng)性邊界條件下由真實(shí)結(jié)構(gòu)模型與體胞模型模擬得到數(shù)據(jù)經(jīng)等效宏觀方法和一階均勻化方法計(jì)算得到真應(yīng)力-真應(yīng)變曲線(xiàn)與試驗(yàn)結(jié)果的對(duì)比Fig.7 Comparison of true stress-true strain curves calculated by equivalent macro method and first-order homogenization method from simulation by real structural model (a) and body cell model (b-c) under symmetric boundary condition with test results: (b) circular particle with random distribution and (c) irregular particle with random distribution

2 結(jié)果與討論

2.1 邊界條件對(duì)模擬結(jié)果的影響

在3種邊界條件下用真實(shí)結(jié)構(gòu)模型進(jìn)行沿水平方向的單軸拉伸有限元模擬,得到的拉伸真應(yīng)力-真應(yīng)變曲線(xiàn)如圖8所示。由圖8可以看出,不同邊界條件下模擬得到的真應(yīng)力-真應(yīng)變曲線(xiàn)基本重合,說(shuō)明邊界條件對(duì)模擬結(jié)果影響較小。

圖8 不同邊界條件下用真實(shí)結(jié)構(gòu)模型單軸拉伸模擬得到的真應(yīng)力-真應(yīng)變曲線(xiàn)Fig.8 True stress-true strain curves obtained by simulating uniaxial tension with real structure model under different boundary conditions

由圖9可以看出,邊界條件對(duì)約束位置的應(yīng)力、應(yīng)變場(chǎng)具有一定影響,由于同質(zhì)化邊界條件和對(duì)稱(chēng)性邊界條件對(duì)垂直加載方向的約束能力弱,使模型變形過(guò)程中沿垂直加載方向有一定的收縮,從而導(dǎo)致模型的最大應(yīng)力與最大等效塑性應(yīng)變位置具有較大的差異。周期性邊界條件下模型變形均勻,在邊界位置應(yīng)力場(chǎng)連續(xù)性較好,與復(fù)合材料內(nèi)部實(shí)際情況相吻合。對(duì)稱(chēng)性邊界條件與周期性邊界條件下的應(yīng)力與應(yīng)變?cè)茍D有較高的相似性,因此單軸拉伸試驗(yàn)?zāi)M時(shí)可用對(duì)稱(chēng)性邊界條件代替周期性邊界條件施加載荷。

圖9 不同邊界條件下用真實(shí)結(jié)構(gòu)模型模擬單軸拉伸真應(yīng)變1%時(shí)的Mises應(yīng)力云圖和等效塑性應(yīng)變?cè)茍DFig.9 Mises stress countour (a-c) and equivalent plastic strain countour (d-f) by simulating uniaxial tension with real structure model for true strian of 1% under different boundary conditions: (a,d) periodic boundary condition; (b,d) symmetric boundary condition and (c,f) homogenization boundary condition

2.2 2種建模方法下模擬結(jié)果的對(duì)比

由圖10可以看出,對(duì)2種模型模擬單軸拉伸得到的真應(yīng)力-真應(yīng)變曲線(xiàn)間差異較小,且均與試驗(yàn)結(jié)果吻合較好,相對(duì)誤差小于5%。根據(jù)真應(yīng)力-真應(yīng)變曲線(xiàn)計(jì)算得到的屈服強(qiáng)度及彈性模量如表1所示。

圖10 對(duì)稱(chēng)性邊界條件下用不同模型模擬單軸拉伸得到的真應(yīng)力-真應(yīng)變曲線(xiàn)與試驗(yàn)結(jié)果的對(duì)比Fig.10 Comparison of true stress-true strain curves by simulating uniaxial tension with different models under symmetric boundary condition with test results

表1 對(duì)稱(chēng)性邊界條件下用不同模型模擬單軸拉伸得到的彈性模量和屈服強(qiáng)度以及試驗(yàn)結(jié)果

由表1可以看出,對(duì)2種模型模擬單軸拉伸后得到的彈性模量及屈服強(qiáng)度與試驗(yàn)結(jié)果均相差較小,相對(duì)誤差小于10%。從工程應(yīng)用的角度來(lái)看,這兩種模型都可以較好地模擬復(fù)合材料在單軸拉伸過(guò)程中的應(yīng)力-應(yīng)變曲線(xiàn)。對(duì)真實(shí)結(jié)構(gòu)模型模擬單軸拉伸得到的彈性模量及屈服強(qiáng)度與試驗(yàn)結(jié)果間的相對(duì)誤差最小,分別為3.19%與5.63%,說(shuō)明真實(shí)結(jié)構(gòu)模型在彈性模量及屈服強(qiáng)度模擬方面的可靠性更高。

由圖11可以看出,顆粒分布與形狀的不同導(dǎo)致基體和顆粒的應(yīng)力場(chǎng)和塑性應(yīng)變場(chǎng)具有明顯的差異。體胞模型對(duì)顆粒的形狀與分布進(jìn)行了一定近似處理,而真實(shí)結(jié)構(gòu)模型中的顆粒具有更多的棱角,顆粒與基體之間有較多的接觸面積,因此相比于近似化的體胞模型其內(nèi)部應(yīng)力集中區(qū)域較多,這也是導(dǎo)致真實(shí)結(jié)構(gòu)模型下在相同應(yīng)變時(shí)的應(yīng)力略高于體胞模型的原因。復(fù)合材料基體高應(yīng)力區(qū)域多在沿加載方向的顆粒之間,并且顆粒間距越小,其中間的基體等效應(yīng)力越大;不同模型中基體較大的等效塑性應(yīng)變主要集中在高應(yīng)力區(qū)域,并且等效塑性應(yīng)變?cè)茍D與各自的等效應(yīng)力云圖具有相似性;2種模型中顆粒的高應(yīng)力區(qū)域多位于顆粒的尖角處和顆粒間距較小處。由表2可知:真實(shí)結(jié)構(gòu)模型顆粒的最大等效應(yīng)力分別為隨機(jī)分布圓形顆粒體胞模型和隨機(jī)分布不規(guī)則顆粒體胞模型的2.22倍和2.23倍;真實(shí)結(jié)構(gòu)模型基體的最大等效塑性應(yīng)變分別為隨機(jī)分布圓形顆粒體胞模型和隨機(jī)分布不規(guī)則顆粒體胞模型的1.38倍和2.58倍。由此可見(jiàn),2種模型模擬顆粒增強(qiáng)金屬基復(fù)合材料局部微觀等效應(yīng)力場(chǎng)及應(yīng)變場(chǎng)的結(jié)果有明顯差異,這是由于在顆粒增強(qiáng)金屬基復(fù)合材料實(shí)際服役過(guò)程中,基體常因產(chǎn)生較大塑性變形而失效,而顆粒則容易在其高應(yīng)力區(qū)域的尖角處發(fā)生斷裂。因此,在對(duì)復(fù)合材料的斷裂、失效等行為進(jìn)行模擬時(shí),真實(shí)結(jié)構(gòu)模型和體胞模型的模擬結(jié)果會(huì)存在較大的差異。

圖11 對(duì)稱(chēng)性邊界條件下用不同模型模擬單軸拉伸1%時(shí)基體的Mises應(yīng)力云圖和等效塑性應(yīng)變?cè)茍D以及顆粒Mises應(yīng)力云圖Fig.11 Mises stress countour (a-c) and equivalent plastic strain countour (d-f) of matrix (a-f) and Mises stress countour of particles (g-i) by simulating uniaxial tension with different models for true strain of 1% under symmetric boundary condition:(a, d, g) body cell model of circular particle with random distribution; (b,e,h) body cell model of irregular particle with random distribution and (c,f,i)real structure model

表2 對(duì)稱(chēng)性邊界條件下用不同模型模擬單軸拉伸真應(yīng)變?yōu)?%時(shí)基體的最大等效應(yīng)力、最大等效塑性應(yīng)變及顆粒的最大等效應(yīng)力

3 結(jié) 論

(1) 對(duì)原位生成法制備的TiB2顆粒增強(qiáng)鐵基復(fù)合材料建立了顆粒隨機(jī)分布的體胞模型和真實(shí)結(jié)構(gòu)模型,采用ABAQUS軟件根據(jù)不同模型對(duì)單軸拉伸進(jìn)行有限元模擬后,采用等效宏觀方法獲取的真應(yīng)力-真應(yīng)變曲線(xiàn)與試驗(yàn)結(jié)果基本吻合,相對(duì)誤差小于5%,驗(yàn)證了這2種模型的合理性。

(2) 同質(zhì)化邊界條件、對(duì)稱(chēng)性邊界條件和周期性邊界條件下用真實(shí)結(jié)構(gòu)模型模擬拉伸得到的真應(yīng)力-真應(yīng)變曲線(xiàn)基本重合,說(shuō)明邊界條件對(duì)模擬結(jié)果影響較小;對(duì)稱(chēng)性邊界條件下的應(yīng)力與應(yīng)變?cè)茍D與周期性邊界條件相似,在單軸拉伸試驗(yàn)?zāi)M時(shí)可用對(duì)稱(chēng)性邊界條件代替周期性邊界條件施加載荷。

(3) 用真實(shí)結(jié)構(gòu)模型和體胞模型模擬單軸拉伸得到的真應(yīng)力-真應(yīng)變曲線(xiàn)間差異較小,且均與試驗(yàn)結(jié)果吻合較好,相對(duì)誤差小于5%,得到的彈性模量及屈服強(qiáng)度與試驗(yàn)結(jié)果均相差較小,相對(duì)誤差小于10%,且真實(shí)結(jié)構(gòu)模型模擬單軸拉伸得到的彈性模量與屈服強(qiáng)度的誤差小于體胞模型,說(shuō)明真實(shí)結(jié)構(gòu)模型在彈性模量及屈服強(qiáng)度模擬方面的可靠性更高。

(4) 用不同模型模擬得到復(fù)合材料基體與顆粒在單軸拉伸時(shí)的局部微觀等效應(yīng)力場(chǎng)及應(yīng)變場(chǎng)有明顯差異,真實(shí)結(jié)構(gòu)模型模擬顆粒的最大等效應(yīng)力分別為隨機(jī)分布圓形顆粒體胞模型和隨機(jī)分布不規(guī)則顆粒體胞模型的2.22倍和2.23倍,基體的最大等效塑性應(yīng)變分別為這2種體胞模型的1.38倍和2.58倍。

猜你喜歡
復(fù)合材料有限元模型
一半模型
重要模型『一線(xiàn)三等角』
重尾非線(xiàn)性自回歸模型自加權(quán)M-估計(jì)的漸近分布
民機(jī)復(fù)合材料的適航鑒定
復(fù)合材料無(wú)損檢測(cè)探討
3D打印中的模型分割與打包
磨削淬硬殘余應(yīng)力的有限元分析
TiO2/ACF復(fù)合材料的制備及表征
基于SolidWorks的吸嘴支撐臂有限元分析
RGO/C3N4復(fù)合材料的制備及可見(jiàn)光催化性能
主站蜘蛛池模板: 国产在线小视频| 大香伊人久久| 国产麻豆精品久久一二三| 美女毛片在线| 无码福利日韩神码福利片| 欧美国产视频| 尤物特级无码毛片免费| 久久中文无码精品| 国产精品无码影视久久久久久久 | 国产办公室秘书无码精品| 国产成人高清亚洲一区久久| 免费高清自慰一区二区三区| 久久国产精品影院| 免费国产福利| 高潮爽到爆的喷水女主播视频| av在线无码浏览| 国产成人欧美| 国产尤物视频网址导航| 精品国产自在现线看久久| 国产在线观看精品| 手机永久AV在线播放| 久久中文字幕2021精品| 在线毛片网站| 欧美成人手机在线视频| 亚洲 欧美 日韩综合一区| 在线综合亚洲欧美网站| 精品精品国产高清A毛片| 亚洲IV视频免费在线光看| 无码免费视频| 国产美女91视频| 亚洲欧美人成人让影院| 国产91特黄特色A级毛片| 午夜日本永久乱码免费播放片| 五月丁香在线视频| 91po国产在线精品免费观看| 一级毛片免费高清视频| 国产高清精品在线91| 亚洲综合18p| 国产日韩丝袜一二三区| 40岁成熟女人牲交片免费| 91国内外精品自在线播放| 亚洲精品国产日韩无码AV永久免费网| 国产黄在线观看| 狠狠色丁香婷婷| 日本精品视频一区二区| 五月天久久综合| 91麻豆国产视频| 97se亚洲综合在线| 国产精品成人久久| 欧美午夜网| 亚洲丝袜中文字幕| 少妇露出福利视频| 欧美亚洲日韩不卡在线在线观看| 无码日韩精品91超碰| 色婷婷天天综合在线| 中文无码日韩精品| 亚洲色图欧美| 成人一区在线| 亚洲一级无毛片无码在线免费视频| 久久午夜夜伦鲁鲁片不卡| 色婷婷在线影院| 国产麻豆另类AV| 久久免费观看视频| 精品91自产拍在线| 国产肉感大码AV无码| 99草精品视频| 亚洲欧美成人网| 一级做a爰片久久免费| 动漫精品啪啪一区二区三区| 99精品影院| 尤物国产在线| 国产精品主播| 国产成人毛片| 91小视频在线观看| 亚洲看片网| 熟妇丰满人妻| 久久久国产精品无码专区| a级毛片网| 国产成人精品男人的天堂| 亚洲视频免| 日韩欧美国产另类| 亚洲综合片|