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

細編穿刺C/C復合材料熱導率數值模擬*

2017-06-19 19:09:41童明波王玉青朱書華
固體火箭技術 2017年3期
關鍵詞:復合材料有限元模型

梁 珩,童明波,王玉青,朱書華

(南京航空航天大學 飛行器先進設計技術國防重點學科實驗室,南京 210016)

細編穿刺C/C復合材料熱導率數值模擬*

梁 珩,童明波,王玉青,朱書華

(南京航空航天大學 飛行器先進設計技術國防重點學科實驗室,南京 210016)

根據細編穿刺復合材料的細觀和微觀結構,分別建立了纖維束和細編穿刺單胞有限元模型。采用周期性非絕熱溫度邊界條件,計算了纖維束和材料整體的等效熱導率。計算結果與經驗公式比較,具有高度的一致性。在此基礎上,進一步研究了纖維體積分數、基體和纖維熱導率對材料熱導率的影響。結果表明,隨著纖維含量的增加,材料兩個方向熱導率均有不同程度的下降,且差異逐漸減小,且基體對熱導率影響作用較大。文中采用的模型和周期性邊界條件與理論預期符合較好,為材料熱學和熱力耦合問題的分析提供了有用參考。

細編穿刺C/C復合材料;有限元;單胞;周期性邊界條件;熱導率

0 引言

隨著現代工業水平、國防技術的快速發展,C/C復合材料因其優異的比剛度、比強度和耐高溫性能,而廣泛用于航空航天等領域。C/C復合材料也從最初的兩向發展為三向、四向等多向編織C/C復合材料,從單一功能材料發展成為多功能材料[1]。復合材料性能的研究也逐漸得到了人們的重視。目前,針對C/C復合材料剛度和強度方面的相關研究較多,對熱物理性能的研究則鮮有報道。由于C/C復合材料結構復雜,工作環境溫度較高,通過試驗獲得材料熱物理性能費用較高,難度較大。因此,通過有限元分析,對材料的熱物理性能進行充分研究,能更好地完成復合材料結構設計工作。

對C/C復合材料的熱物理性能,國外已進行了大量的研究[2-4]。Rolfes R和Hammershmidt U[5]對纖維加強復合材料的橫向熱導率進行了深入研究,提出了多種經驗公式,并與試驗結果進行了詳細比對。Klett J W等[6]建立了一套有限元模型,對影響C/C復合材料的多種因素進行了分析,并準確預測了材料的橫向和縱向熱導率。Grujicic M等[7]針對不同基體和纖維的三維編織C/C復合材料建立了有限元模型,對其熱導率進行了對比分析。結果表明,纖維和基體對材料的熱導率都有重要的影響,纖維基體中的孔隙和裂紋會顯著下降材料熱導率。Hongzhou Li等[8]將周期性邊界條件與單胞模型相結合,建立了微觀和細觀兩級單胞模型,對平紋機織復合材料的熱導率進行了研究。Jian-Jun Gou等[9]針對三維四向編織復合材料,建立了一套有限元分析方法,對其熱導率進行了研究,并與試驗結果進行對比分析。

國內關于復合材料熱導率的研究資料則相對較少。程偉等[10]采用實驗和有限元法,研究了三維四向編織復合材料的等效熱物理性能,為后續的理論分析和數值模擬提供了實驗參考。Liu[11]和李典森等[12]分別對三維四向和五向編織復合材料的熱傳導性能進行了研究。然而,上述文獻中,對平行于熱傳導方向的側面均采用絕熱邊界條件,限制了胞體之間的熱量交換,與實際的情況不太相符。夏彪和盧子興[13]采用周期性的非絕熱溫度邊界條件和位移邊界條件,計算了三維四向和五向編織復合材料的整體等效熱傳導系數和熱膨脹系數,計算結果同已有文獻相比,與實驗值符合得更好。目前,關于C/C復合材料熱導率的研究,已經取得了一定的成果,但仍需深入的對合理的模擬方法進行研究,對指導材料結構與工藝設計,為工程選擇復合材料細觀結構提供依據,具有非常重要的意義。

本文針對細編穿刺C/C復合材料,分別建立了微觀和細觀有限元單胞模型,結合周期性邊界條件在單胞周圍施加非絕熱的溫度邊界條件,對材料的熱導率進行了預測。并進一步分析了纖維體積分數、基體和纖維熱導率對材料熱導率的影響規律。

1 分析方法

1.1 周期性邊界條件

在非均質材料的細觀有限元分析中,單胞模型被廣泛用來計算材料的等效參數[14]。為獲得較準確的結果,邊界條件選取的合理性至關重要。對于機織復合材料,假設纖維在基體中是均勻且按一定規則分布,則其細觀結構具有明顯的周期性,宏觀結構即可視為相同的單胞結構按照一定的周期性排列而成。本文研究的細編穿刺材料為三向細編穿刺結構,織物的X-Y向為緞布炭纖維疊層,Z向為穿刺的炭纖維束。材料具有明顯的平移對稱性,選取最小平移對稱單位作為該材料的代表性體積單元進行研究,如圖1所示。

在平移對稱變換下,模型的溫度分布也應滿足平移對稱條件。假設滿足平移對稱條件的單胞為六面體,如圖2所示,0≤x≤a,0≤y≤b,0≤z≤h。周期性溫度邊界條件的推導過程如下所示:

根據模型的平移對稱性,單胞中任意一點P′(x′,y′,z)都可通過其他單胞上的對應點P(x,y,z)在x、y、z方向上經過若干次平移變換得到,即

(1)

其中,i,j,k=0,±1,±2,…表示點P沿著x、y、z軸平移到P′的距離。

對于點P′(x′,y′,z)和P(x,y,z)處的溫度場,存在如下關系:

(2)

將式(1)代入式(2)中,得

(3)

根據平移對稱性,六面體的三對平行面都可看作是前一個單胞的面沿垂直于該面的坐標軸,通過一次平移變換到后一個單胞中。因此,可得到六面體三隊平行面之間的溫度約束條件:

(4)

其中,下標x=a和x=0等代表對應面上的普通節點,并不包括單胞邊上的點和頂點。

由于單胞的邊同時屬于兩個相鄰的平面,則每個邊上的點應同時滿足兩個面的位移約束方程。因此,需重新推導邊界條件。給六面體單胞的所有邊和頂點進行編號,如圖2所示。

為得到單胞各邊之間的約束條件,可將單胞的12條邊分成3組:平行于x軸的邊Ⅸ、Ⅹ、Ⅺ和Ⅻ;平行于y軸的邊Ⅴ、Ⅵ、Ⅶ和Ⅷ,以及平行于z軸的邊Ⅰ、Ⅱ、Ⅲ和Ⅳ。可看出,不同組的邊是相互獨立的,每一組內的任意一條邊都可由組內其他邊通過平移變換轉化得到。以邊Ⅰ為例,邊Ⅱ、Ⅲ、Ⅳ相對于邊Ⅰ的變換方式分別為(i=1,j=0,k=0),(i=1,j=1,k=0)和(i=0,j=1,k=0)。其他兩組邊具有相似的變換方式。因此,可推導出各組之間的溫度約束條件:

(5)

(6)

(7)

以上推導的邊同樣不包含單胞的頂點,因為任意一個單胞的頂點都同屬于3條正相交的邊,且任意一個頂點都可由其他頂點通過平移變換得到。例如,頂點7就可通過頂點1通過變換方式(i=1,j=0,k=1)得到,由此可推導出各個頂點之間的溫度約束關系:

(8)

1.2 宏觀本構關系與等效參數的求解

本文所采用的單胞穩態熱傳導分析控制方程為

(9)

式中T為溫度;λx、λy、λz、λxy等分別為各個方向的熱導率。

由控制方程(9)和單胞表面的邊界條件,即可求解出單胞的溫度分布,x、y、z方向上的熱通量與溫度梯度的關系,由式(10)給出:

(10)

在邊界條件中,引入宏觀熱通量作為系統的一個獨立自由度,則可用參考點上的集中熱流來表示宏觀熱通量,這個宏觀熱通量與集中熱流之間的關系可通過能量公式等效得到。例如,將x方向上的集中熱流Qx加載在單胞的相應參考點上,其他的自由度都不約束,則在單位時間內,這個熱流產生的能量為

QxΔt=ρcVTx

(11)

式中ρ為材料密度;c為材料比熱容;V為單胞體積;Tx為x方向上溫度變量。

公式兩邊同時關于t求導,可得

(12)

又因為

(13)

所以

(14)

(15)

將式(15)代入式(14)中,可得:

(16)

所以,可推出:

(17)

基于上述公式,用Python語言編寫了溫度周期性邊界條件,以顯示約束方式施加到模型對應的節點上。在纖維束單胞模型中,施加不同方向的集中熱流,計算得到纖維束不同方向的熱導率λ1、λ2、λ3。再將這些結果作為細編穿刺單胞模型的材料參數進行輸入,計算得到細編穿刺材料各個方向的熱導率λx、λy、λz。

2 有限元模型

2.1 纖維束單胞幾何模型

在纖維束理論分析中,假設纖維束為纖維在基體中均勻規則分布的單向纖維增強復合材料,纖維體積分數即為其堆積密度。假設纖維在基體中的分布為正六邊形分布,如圖3所示。據此假設,可建立正六棱柱單胞模型。然而,正六棱柱單胞模型復雜的外形,使得施加邊界條件相對困難。所以,將正六棱柱模型拓展為圖3中矩形區域作為纖維束單胞,建立厚度方向為單位長度的矩形纖維束單胞模型,并施加周期性溫度邊界條件。由于不同的纖維體積含量會導致纖維束熱傳導率的不同,本文分別選取10%~85%纖維體積含量的纖維束進行研究。

(18)

材料采用T300炭纖維,單絲半徑r=3.8 μm ,帶入式(18)中求得a=16.18 μm,b=9.34 μm。

2.2 細編穿刺單胞幾何模型

細編穿刺織物是機織碳布和正交非織造三向織物的組合形式,如圖5所示。材料的成型工藝是將炭纖維經整經、穿綜、卷緯等工序在織布機上織X-Y向互相交織的平面碳布,再根據織物Z向間距和尺寸,確定Z向鋼針矩陣,碳布逐層穿刺,通過加壓控制保證碳布層密度和Z向垂直度,達到設計高度后,用炭纖維置換Z向鋼針形成織物。

由于復合材料平面結構是由一層經紗和一層緯紗相互交織而成,因此在其面內具有嚴格的周期性。在厚度方向上,由于細編穿刺材料是由多層壓制穿刺而來,假設各層在厚度方向上均勻鋪設,忽略表層的影響,則材料內部各層在3個方向上均具有周期性。為了方便有限元模擬,假設Z向纖維束的橫截面近似圓形,X-Y向纖維束的橫截面近似為橢圓形,且都為密實實體,各向纖維束之間的空隙全部由基體填充。假設X-Y向纖維交織方式均按式(19)給定的余弦公式排列。建立的細編穿刺單胞有限元模型如圖6所示。

(19)

經電鏡掃描(SEM)分析,得到單胞的尺寸a=b=2.60 mm,h=0.44 mm,X-Y向纖維橫截面尺寸w=0.61mm,c=0.20 mm,纖維間距u=v=0.69 mm,Z向纖維橫截面半徑r=0.29mm,長度與單胞厚度尺寸相同。通過計算得出細編穿刺單胞模型中,纖維束體積含量Vy=48.48%。由于纖維是各向異性材料,而纖維束在碳布中呈交織的彎曲狀態,因此需在纖維束中定義局部坐標系,使Z軸始終與纖維束方向平行。

3 結果與分析

3.1 纖維束單胞模型計算結果比較

在纖維束中,軸向的熱導率λ3可通過混合法則計算得到:

(20)

式中λfa和λm分別為纖維軸向熱導率和基體的熱導率。

纖維束橫向熱導率的情況則相對復雜,為簡化計算,假設纖維束為橫觀各向同性材料,即

(21)

同時為方便有限元模擬,以及計算準則的需要,提出如下假設[5]:

(1)纖維的橫截面為圓形;

(2)所有的纖維都完全相同;

(3)纖維均勻的分布在基體中;

(4)基體的熱導率是各向同性的,纖維的熱導率在面內是各向同性的。

基于上述假設,1982年,Rayleigh[5,15]在三角級數理論基礎上,提出了一種經驗公式用于計算纖維束橫向熱導率:

(22)

其中

(23)

S4=0.032 350 24

(24)

Hermans[5,16]在保留了Rayleigh理論的第一要素后,提出一種計算準則用于計算纖維束橫向熱導率。

(25)

式中λft為纖維橫向熱導率。

式(25)僅針對纖維與基體為同心圓的復合材料,因此也被稱為自洽公式,如圖7所示。

計算采用的纖維與基體的參數如表1所示。

表1 碳基體和T300纖維的熱導率[17]

不同纖維體積含量的纖維束熱導率有限元預測結果如表2所示。有限元預測結果與法則計算結果對比如圖8所示。可看出,有限元預測纖維束軸向和橫向的熱導率結果與法則計算結果高度重合。隨著纖維含量的增加,纖維束軸向熱導率λ3線性升高。這是因為纖維軸向熱導率比基體高,纖維束軸向熱量傳導主要由纖維承擔。而纖維束橫向熱導率λ1=λ2則逐漸降低,主要因為纖維橫向熱導率相比基體低得多,使得基體在纖維束橫向熱量傳導中占主導地位,隨著纖維體積含量的增加,基體含量減少,使得橫向熱導率相應降低。70%纖維體積含量的纖維束單胞模型熱通量分布圖如圖9所示。

表2 不同纖維含量體積的纖維束熱導率

為進一步驗證周期性邊界條件的可靠性,參照文獻[18]建立瀝青基C/C復合材料單胞模型,如圖10所示。采用周期性邊界條件,忽略界面相的影響,對T300纖維束進行了熱導率的模擬。模擬結果如表3所示。由于沒有考慮纖維與基體界面相的影響,模擬得到的結果比試驗值偏大。綜合經驗公式和試驗對比,認為本文所建模型和周期性溫度邊界條件可較準確地預測纖維束軸向和橫向熱導率。

參數FEM/[W/(m·K)]試驗結果/[W/(m·K)]誤差/%λ1=λ288.4374.3018.98λ3289.61253.8614.08

3.2 纖維含量對細編穿刺復合材料熱導率的影響

為研究纖維體積含量對細編穿刺復合材料熱導率的影響,采用70%纖維體積含量的纖維束,通過改變纖維束經密和緯密改變纖維體積含量,建立如圖11所示的一組模型。

不同纖維體積含量的細編穿刺復合材料的熱導率如圖12所示。隨著纖維束體積含量的增加,單胞面內熱導率和Z向熱導率均有不同程度的下降。這主要是由于C/C 復合材料作為一元復合材料,其熱導率與纖維含量和組元熱導率有關,滿足混合規則。所以,熱導率相對較大的基體體積含量越少,材料的熱導率就會越小。其次,由于基體熱導率小于纖維軸向熱導率卻遠大于纖維橫向熱導率。使得復合材料主要沿纖維束軸向方向進行熱量的傳導,纖維束橫向方向則對熱量的傳導起到一定的阻礙作用。因為細編穿刺復合材料Z向纖維束含量遠小于X-Y向纖維束含量,X-Y向纖維束對縱向熱量的傳導阻礙作用較大。所以,Z向熱導率要小于面內熱導率。從圖12還可看出,隨著纖維體束積含量的增加,面內熱導率與Z向熱導率差距逐漸縮小,這為材料的設計提供一種實用方法,當材料需要在2個方向有不同的熱導率比例時,可通過調節材料經密緯密來實現。

圖11 纖維束體積含量不同的細編穿刺模型

Fig.11 Fine weave pierced C/C composite of different yarn volume fraction

3.3 細編穿刺復合材料熱通量分布

以48.48%纖維束體積含量的單胞為例,在X方向施加集中熱流 ,熱通量分布圖如圖13所示。可看出,X向纖維束熱通量較大,熱通量主要集中在經紗和緯紗的交織處,以及周圍的基體。纖維中熱通量最大值為1 739 W/m2,周圍基體的熱通量最大值為1 459 W/m2。

圖14為在單胞Z向施加集中熱流的熱通量分布圖。熱量主要集中在Z向纖維束及周圍的基體上,Z向集中熱流下纖維束熱通量最大值為26.85 W/m2,基體為25.17 W/m2。

由圖13和圖14可知,細編穿刺復合材料的熱載荷主要由相應方向的纖維束承擔,同時纖維束周圍基體也承擔了較多熱載荷。

3.4 參數研究

復合材料的熱導率是由纖維、基體和機織方式共同決定的。采用不同的纖維和基體,復合材料的熱導率也會有一定的差別。因此,需要對復合材料的敏感性進行研究。通過分別改變纖維絲和基體的熱導率,對細編穿刺碳碳復合材料單胞模型進行有限元模擬,模擬結果如圖15和圖16所示。

圖15中,纖維絲熱導率不變,將基體的熱導率變為原來的1/10為7.04 W/(m·K)。可看出,隨著碳基體熱導率的減小,單胞整體熱導率也較大程度下降。當基體熱導率小于纖維熱導率時,隨著纖維含量的增加,單胞面內和縱向熱導率增加。由于細編穿刺復合材料增加了Z向的纖維束,所以縱向熱導率與面內熱導率相差不多。

圖16中,基體和纖維橫向熱導率保持不變,軸向熱導率提高到200 W/(m·K)。由圖16可知,隨著軸向熱導率的增加,單胞整體各個方向的熱導率均有不同程度的增加,但增加幅度較小。由于纖維橫向熱導率依然小于基體熱導率,單胞面內熱導率依然隨纖維含量的增加呈下降趨勢。Z向熱導率則小幅上揚。

4 結論

(1)本文所采用的周期性邊界條件保證了在單元胞體之間的熱量交換,與實際情況較為相符,預測結果與經驗公式吻合程度較高。為后續復合材料的結構性能研究及優化奠定了一定的研究基礎。

(2)針對不同機織密度的單胞進行有限元分析,研究表明,隨著纖維含量的增加,材料兩個方向熱導率均有不同程度的下降,且差異逐漸減小,為材料的設計提供一種實用方法。

(3)針對不同纖維含量的單胞的材料分析表明,材料對基體熱導率敏感性較高。當基體熱導率大于纖維橫向熱導率時,材料Z向熱導率和面內熱導率隨著纖維含量的增加而下降,其中面內熱導率下降趨勢較快。反之,熱導率隨纖維含量增加呈上升趨勢。

[1] 石振海,李克智,李賀軍,等.航天器熱防護材料研究現狀與發展趨勢[J].材料導報,2007,21(8):15-18.

[2] Zou M Q,Yu B M,Zhang D M,et al.Study on optimization of transverse thermal conductivities of unidirectional composites [J].Journal of Heat Transfer,2003,125(2):980-987.

[3] Goo N S,Woo K.Measurement and prediction of effective thermal conductivities of woven fabric composites[J].International Journal of Modern Physics B,2003,17(8-9):1808-1813.

[4] Woo K,Goo N S.Thermal conductivity of carbon-phenolic 8-harness satin weave composites[J].Composite Structures,2004,66(1-4):521-526.

[5] Rolfes R,Hammerschmidt U.Transverse thermal conductivity of CFRP laminates: a numerical and experimental validation of approximation formulae[J].Composites Science and Technology,1995,54 (1):45-54.

[7] Klett J W,Edie D D.Flexible towpreg for the fabrication of high thermal conductivity carbon/carbon composites[J].Carbon,1995,33(10):1485-1503.

[6] Grujicic M,Zhao C L,Dusel E C,et al.Computational analysis of the thermal conductivity of the carbon-carbon composite materials[J].Journal of Materials Science,2006,41(24):8244-8256.

[8] Li H Z,Li S G,Wang Y C.Prediction of effective thermal conductivities of woven fabric composites using unit cells at multiple length scales[J].Journal of Materials Research,2011,26(3):384-394.

[9] Gou J J,Zhang H,Dai Y J,et al.Numerical prediction of effective thermal conductivities of 3D four-directional braided composites[J].Composite Structures,2015,125:509-519.

[10] 程偉,趙壽根,劉振國,等.三維四向編織復合材料等效熱特性數值分析和試驗研究[J].航空學報,2002,23(2):102-105.

[11] Liu Z G,Zhang H G,Lu Z X.Investigation on the thermal conductivity of 3-dimensional and 4-directional braided composites[J].Chinese Journal of Aeronautics,2007,20(4):327-331.

[12] 李典森,盧子興,劉振國,等.三維五向編織復合材料導熱性能的有限元分析[J].航空動力學報,2008,23(8):1455-1460.

[13] 夏彪,盧子興.三維編織復合材料熱物理性能的有限元分析[J].航空學報,2011,32(6):1040-1049.

[14] 潘清.細觀力學可視參數化建模及一種新型擴展有限元方法研究[D].南京:南京航空航天大學,2012.

[15] Lord Rayleigh.On the influence of obstacles arranged in rectangular order upon the properties of a medium[J].Phil.Mag,1982,34:481-502.

[16] Hermans J J.The elastic properties of fiber reinforced materials when the fibers are aligned[C]//Proc.Konigl.Nederl.A kad.Weteschappen,1967,B70,l.

[17] 張沫.多向編織復合材料熱物理性能研究[D].哈爾濱:哈爾濱工業大學,2008.

[18] 支佳運.瀝青基碳/碳復合材料熱導率的理論研究[D].北京:華北電力大學,2015.

(編輯:薛永利)

Simulation on effective thermal conductivity of fine weave pierced C/C composite

LIANG Heng,TONG Ming-bo,WANG Yu-qing,ZHU Shu-hua

(Ministerial Key Discipline Laboratory of Advanced Design Technology of Aircraft,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,China)

Based on the meso and micro structures of the fine weave pierced C/C composite,finite element models of a single braiding yarn and C/C composites were established, respectively.By using the periodic non-adiabatic temperature boundary condition, the effective thermal conductivity of yarns and composites were calculated.The calculation results were highly consistent with empirical formula.Furthermore, the law of effective thermal conductivity by material parameters(i.e.fiber volume fraction,fiber and matrix thermal conductivity)was studied. With the increase of the fiber fraction,the transverse and longitudinal effective thermal conductivities decrease with reduced discrepancy.The thermal conductivity of the matrix has much more effects on the global effective thermal conductivity.Based on the FE model and boundary conditions proposed in this paper,the calculation results are in good agreement with the theory, which provides a useful reference for the analysis of thermal and thermo-mechanical coupling.

fine weave pierced C/C composites;FEM;unit cell;periodic boundary condition;thermal conductivity

2016-01-08;

2016-02-14。

國家自然科學基金(11302105);中國航空科學基金(2014ZF52074)。

梁珩(1989—),男,博士,研究方向為飛行器結構設計。E-mail: leunghengh@sina.com

童明波,男,教授,研究方向為飛行器結構設計。E-mail:tongw@nuaa.edu.cn

V258

A

1006-2793(2017)03-0364-08

10.7673/j.issn.1006-2793.2017.03.017

猜你喜歡
復合材料有限元模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
民機復合材料的適航鑒定
復合材料無損檢測探討
電子測試(2017年11期)2017-12-15 08:57:13
3D打印中的模型分割與打包
磨削淬硬殘余應力的有限元分析
TiO2/ACF復合材料的制備及表征
應用化工(2014年10期)2014-08-16 13:11:29
基于SolidWorks的吸嘴支撐臂有限元分析
RGO/C3N4復合材料的制備及可見光催化性能
主站蜘蛛池模板: 麻豆精品视频在线原创| 久久久久国产精品熟女影院| 亚洲成在线观看| 亚洲欧洲天堂色AV| 中文字幕无码av专区久久| 麻豆AV网站免费进入| 国产精品一区二区不卡的视频| 亚洲国产天堂久久综合| 国产剧情无码视频在线观看| 高清乱码精品福利在线视频| 日本欧美成人免费| 日本一区二区三区精品国产| 在线免费a视频| 亚洲h视频在线| 在线观看欧美精品二区| 国产精品免费p区| 日本午夜三级| 在线亚洲小视频| 狠狠综合久久| 亚洲天堂视频网| 久久五月天综合| 久久久黄色片| 精品无码一区二区在线观看| 国产免费怡红院视频| 3344在线观看无码| 成年人福利视频| 久久99国产乱子伦精品免| 成人综合久久综合| 国产你懂得| 99热国产这里只有精品9九| 精品久久久久久成人AV| 91福利在线观看视频| 亚洲视频三级| 少妇精品在线| 五月婷婷伊人网| 色综合成人| 成人噜噜噜视频在线观看| 亚洲a级毛片| 国产第一页屁屁影院| 色窝窝免费一区二区三区| 超薄丝袜足j国产在线视频| 国产免费羞羞视频| 亚洲人成在线精品| 曰韩免费无码AV一区二区| 国产www网站| 亚洲无码日韩一区| 亚洲AV无码乱码在线观看裸奔| 999国产精品| 久久a毛片| 亚洲男人的天堂在线| 天天爽免费视频| 无码精油按摩潮喷在线播放| 香蕉久人久人青草青草| 国产av无码日韩av无码网站 | 伊人久久综在合线亚洲2019| 欧美一区二区人人喊爽| 亚洲午夜福利精品无码不卡 | 97免费在线观看视频| 在线观看91香蕉国产免费| 亚洲高清日韩heyzo| 不卡无码网| 不卡色老大久久综合网| 午夜免费视频网站| 欧美一区二区精品久久久| 亚洲中文精品人人永久免费| 亚洲无限乱码一二三四区| 国内自拍久第一页| 欧美一级在线播放| 第一页亚洲| 国内嫩模私拍精品视频| 欧美爱爱网| 中国毛片网| 亚洲午夜天堂| 国产精品网曝门免费视频| 最新日韩AV网址在线观看| 高清无码一本到东京热| 伊人91在线| 一级片免费网站| 久久综合九九亚洲一区| 国产成人精品一区二区不卡| 欧美A级V片在线观看| 在线观看无码a∨|