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

應(yīng)用橋聯(lián)理論與有限元模型的航空復(fù)合材料結(jié)構(gòu)失效分析

2018-07-23 09:39:10蘆冠達(dá)黃爭(zhēng)鳴
航空學(xué)報(bào) 2018年6期
關(guān)鍵詞:復(fù)合材料有限元

蘆冠達(dá),黃爭(zhēng)鳴

同濟(jì)大學(xué) 航空與力學(xué)學(xué)院,上海 200092

復(fù)合材料因高比強(qiáng)度、高比模量、可設(shè)計(jì)性好等優(yōu)點(diǎn),在航空航天領(lǐng)域得到廣泛應(yīng)用,其極限承載能力是工程設(shè)計(jì)人員最為關(guān)注的一項(xiàng)性能指標(biāo)。Huang的橋聯(lián)模型[1-2],僅需組份材料的性能參數(shù),即可完成對(duì)復(fù)合材料彈塑性應(yīng)力場(chǎng)的計(jì)算,而且所有公式皆為解析式,十分方便于工程應(yīng)用。橋聯(lián)模型誕生不久,受邀參加了第一屆世界范圍破壞評(píng)估(WWFE-I)[3],預(yù)報(bào)精度雖超過了其他參評(píng)的細(xì)觀力學(xué)理論,但總體評(píng)分上尚不如一些宏觀理論[4-5]。

近年,橋聯(lián)模型在復(fù)合材料破壞分析和強(qiáng)度預(yù)報(bào)方面取得重要進(jìn)展,導(dǎo)出了添加纖維后基體中的各應(yīng)力集中系數(shù)[6],修正了層合板最終破壞判據(jù)[7],完善了纖維和基體界面開裂后對(duì)材料承載能力的影響[8],添加了層合板非致命破壞后的部分剛度衰減策略[9]。這些工作使橋聯(lián)理論預(yù)報(bào)WEEI-I的精度顯著提升,整體精度比賽會(huì)第一名Zinoviev理論的精度更高[10],更重要的是,組分材料的輸入?yún)?shù)與其原始性能完全相同,無需由任何復(fù)合材料的實(shí)驗(yàn)數(shù)據(jù)來反演確定輸入?yún)?shù)。除了組分材料的原始性能,僅需提供單向復(fù)合材料的橫向拉伸強(qiáng)度,以確定界面粘接強(qiáng)弱,這給實(shí)際應(yīng)用帶來極大方便。

有限元(FE)數(shù)值模擬已成為工程結(jié)構(gòu)尤其是復(fù)雜結(jié)構(gòu)分析不可或缺的工具,在復(fù)合材料領(lǐng)域也如此。但目前的商用有限元軟件,主要只提供基于宏觀破壞判據(jù)的復(fù)合材料強(qiáng)度計(jì)算模塊,需要使用者提供眾多復(fù)合材料測(cè)試數(shù)據(jù),這給用戶帶來很大不便。歐洲航天局在1992年開始開發(fā)ESAComp軟件,使用經(jīng)典層合板理論計(jì)算各層應(yīng)力值,并把初始層破壞(First Ply Failure)作為最終破壞判據(jù)[11]。2002年,F(xiàn)riebel等開發(fā)了Digimat材料用戶子程序[12],主要基于Mori-Tanaka細(xì)觀力學(xué)方法,但對(duì)復(fù)合材料的破壞分析,依然需要提供復(fù)合材料的實(shí)驗(yàn)數(shù)據(jù)。2015年,韓國(guó)研究人員Lee等運(yùn)用宏觀理論P(yáng)uck失效準(zhǔn)則與損傷有限單元法,通過子程序開發(fā),成功將該理論添加入ABAQUS的計(jì)算中,得到相關(guān)層合板算例的強(qiáng)度譜,并與實(shí)驗(yàn)結(jié)果作了對(duì)比[13]。這些工作豐富了復(fù)合材料結(jié)構(gòu)分析功能,但如何減少對(duì)復(fù)合材料實(shí)驗(yàn)數(shù)據(jù)的依賴尚有不足。

隨著Huang的橋連理論不斷發(fā)展與完善,如何將該理論與CAE(Computer Added Engineering)軟件結(jié)合分析工程結(jié)構(gòu)的破壞與強(qiáng)度問題,成為該理論發(fā)展的一大方向。2007年,Huang首次將橋聯(lián)模型結(jié)合商業(yè)有限元軟件ABAQUS,開發(fā)出用戶子程序UGENS[7,14]。本文在該工作基礎(chǔ)上,將UGENS程序改編,添加了上述橋聯(lián)模型的最新進(jìn)展,并應(yīng)用于WWEF-I考題分析和復(fù)雜層合板結(jié)構(gòu)T型接頭的有限元模擬與強(qiáng)度預(yù)報(bào),計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合良好。

1 理論概述

橋聯(lián)模型UGENS程序基礎(chǔ)詳見文獻(xiàn)[2],為節(jié)省篇幅,本文只簡(jiǎn)要列出各有關(guān)公式。

1.1 均值內(nèi)應(yīng)力計(jì)算

(1)

(2)

(3)

(4)

(5)

(6)

上述各式中Aij為橋聯(lián)矩陣A中的元素,A為[2]

(7)

其中:

(8)

應(yīng)用于層合板分析,有[2]:

(9)

(10)

(11)

(12)

(13)

進(jìn)一步將式(12)代入式(1)和式(2)中,可求出纖維和基體的內(nèi)應(yīng)力增量。

1.2 真實(shí)內(nèi)應(yīng)力

(14)

(15a)

(15b)

(16)

(17)

(18)

式中:

(19)

纖維中的點(diǎn)應(yīng)力均勻,其真實(shí)應(yīng)力與均值應(yīng)力相同。

1.3 破壞判據(jù)

當(dāng)前步纖維和基體中的真實(shí)應(yīng)力以及單層板分擔(dān)的外力按如下更新:

(20)

(21)

(σ)l=(σ)l-1+dσl

(22)

基于當(dāng)前真實(shí)應(yīng)力,若纖維或基體已經(jīng)破壞,則認(rèn)為復(fù)合材料達(dá)到了破壞。程序中,纖維采用修正的最大應(yīng)力破壞判據(jù)[15]:

(23)

(24)

一旦纖維破壞即式(23)成立,認(rèn)為層合板達(dá)到了致命破壞,對(duì)應(yīng)的外載定義為層合板的極限強(qiáng)度,程序終止計(jì)算。基體的破壞式(24)成立只是一種非致命破壞,需折減單層板的剛度。本文采用的剛度衰減格式為

(25)

存在纖維始終不破壞的情況,持續(xù)加載及每一層的非致命破壞,必然會(huì)導(dǎo)致層合板變形過大。為此,在程序中還施加了應(yīng)變控制條件:當(dāng)層合板最大應(yīng)變超過極限應(yīng)變時(shí),也視為發(fā)生了致命破壞,即需滿足

(26)

式中:h為粘結(jié)層厚度。

聚合物復(fù)合材料最大應(yīng)變一般不超過15%[10],一般層板結(jié)構(gòu)破壞時(shí)的極限應(yīng)變不超過10%。以下算例中均取10%作為極限應(yīng)變條件。

1.4 用戶子程序

上述理論發(fā)展均添加到UGENS子程序,與2007版本相比[14],只需額外多輸入一個(gè)單向復(fù)合材料的橫向拉伸強(qiáng)度,以便確定纖維和基體界面何時(shí)開裂,若缺省,認(rèn)為界面始終理想。基體的各應(yīng)力集中系數(shù),根據(jù)纖維和基體的原始性能參數(shù)計(jì)算得到。

2 單層板與層合板強(qiáng)度預(yù)報(bào)——WWFE分析

2.1 算例介紹

有限元建模均采用薄殼單元。只考慮雙軸向正應(yīng)力加載時(shí),在殼體結(jié)構(gòu)左下兩邊施以對(duì)稱邊界條件,如圖1所示;存在剪應(yīng)力的算例,則用左下節(jié)點(diǎn)鉸支的邊界條件,如圖2所示。以此達(dá)到整個(gè)殼體受均勻載荷的目的。

表1 4種復(fù)合材料的組份材料性能參數(shù)[17]Table 1 Constituent properties of four composites[17]

注:F和M分別表示纖維與基體。

表2 4種復(fù)合材料中基體材料的塑性參數(shù)[17]Table 2 Plastic parameters of matrices in four composites[17]

表3 4種基體材料的應(yīng)力集中系數(shù)及界面開裂參數(shù)Table 3 Stress concentration factors and cracking parameters of matrices in four composites

圖1 雙軸加載邊界條件[12]Fig.1 Boundary condition for biaxial loading[12]

圖2 受剪應(yīng)力邊界條件[12]Fig.2 Boundary condition for shear stress[12]

2.2 計(jì)算結(jié)果

根據(jù)表1和表2的輸入數(shù)據(jù)并應(yīng)用表3中的參數(shù),對(duì)WWFE的14個(gè)考題進(jìn)行了分析計(jì)算。為省篇幅,這里只給出最終結(jié)果的得分情況,如表4所示。評(píng)估組織者共設(shè)置了5大類125個(gè)評(píng)分點(diǎn)。當(dāng)理論預(yù)報(bào)值與實(shí)驗(yàn)值的偏差在±10%之內(nèi)時(shí)得分為A,預(yù)報(bào)值與實(shí)驗(yàn)值偏差在±10%~±50%之間得分為B,超出±50%時(shí)得分為C[5]。根據(jù)評(píng)分標(biāo)準(zhǔn),本文有限元求解A+B之和的得分比例為86.4%,遠(yuǎn)高于WWEF評(píng)估中得分最高的Zinoviev理論76.8%的A+B得分比例。證實(shí)了二維橋聯(lián)理論UGENS子程序預(yù)報(bào)復(fù)合材料強(qiáng)度特性的有效性。

表4 WWFE評(píng)估匯總表Table 4 Summary of results for WWFE

注:RC1~RC5依次表示單層板破壞強(qiáng)度分析、層合板初始破壞分析、層合板最終破壞分析、層合板應(yīng)力應(yīng)變曲線分析、定性特征分析。

3 T型接頭受彎失效分析

3.1 材料參數(shù)

T型接頭是復(fù)合材料結(jié)構(gòu)連接的一種常見形式,在航空領(lǐng)域應(yīng)用廣泛,也受到相關(guān)學(xué)者的廣泛研究。為進(jìn)一步驗(yàn)證二維橋聯(lián)理論UGENS子程序的功效,這里對(duì)文獻(xiàn)[18]報(bào)告的T型接頭在彎曲載荷下的極限破壞進(jìn)行了有限元模擬,并與實(shí)驗(yàn)數(shù)據(jù)對(duì)比。

文獻(xiàn)[18]給出了單層板的力學(xué)性能參數(shù)和層合板幾何數(shù)據(jù),但未能提供組分材料的性能參數(shù)。本文分析中,纖維材料參數(shù)取自文獻(xiàn)[19],基體材料參數(shù)取自文獻(xiàn)[20],不計(jì)基體塑性變形,詳見表5。

根據(jù)文獻(xiàn)[18],T型接頭由不同結(jié)構(gòu)層合板共固化而成,本文對(duì)兩個(gè)層合板的結(jié)合處以及曲率劇變(即90°彎角)處導(dǎo)致的富樹脂區(qū)均添加了一層ABAQUS自帶的Cohesive單元。

幾何上,Cohesive單元類似無厚度的平面單元,該單元上可承受3種作用力,即正應(yīng)力tn、剪應(yīng)力ts和tt,分別定義[21]為

(27)

表5 T700/QY9611組分材料的力學(xué)性能[19-20]Table 5 Mechanical properties of T700 fiber and QY9611 matrix[19-20]

圖3 Cohesive單元的線彈性-線軟化本構(gòu)模型Fig.3 Linearly elastic-softening constitutive model for a cohesive element

在圖3所示的線彈性-線軟化本構(gòu)模型中,當(dāng)δ=δ0時(shí)材料屈服,當(dāng)δ=δmax時(shí)材料開裂。假設(shè)粘結(jié)層厚度為h,對(duì)應(yīng)的3個(gè)應(yīng)變?yōu)?/p>

(28)

當(dāng)δ<δ0時(shí),材料為線彈性,此時(shí)有

(29)

當(dāng)δ0≤δ<δmax時(shí),材料處于軟化區(qū),此時(shí)有

(30)

式中:D為損傷因子,其大小由圖3確定。

由于粘結(jié)層破壞往往并非單一開裂破壞模式所致,應(yīng)考慮混合模式的開裂準(zhǔn)則。本文采用ABAQUS提供的二次應(yīng)力判據(jù)和B-K準(zhǔn)則,二次應(yīng)力判據(jù)將正應(yīng)力引起的Ⅰ型斷裂與剪應(yīng)力引起的Ⅱ型斷裂結(jié)合起來,簡(jiǎn)化為一個(gè)應(yīng)力-位移曲線,其中損傷起始的臨界點(diǎn)由式(31)確定,當(dāng)滿足式(31)時(shí)裂紋開始產(chǎn)生[21]。

(31)

鑒于該T型板結(jié)構(gòu)為整體共固化而成,本文將粘結(jié)層取為與純基體的力學(xué)性能參數(shù)相同更為科學(xué),但文獻(xiàn)[18]采用了另一組粘結(jié)層參數(shù),一并列于表6中。分別應(yīng)用兩組粘結(jié)層參數(shù)進(jìn)行了有限元模擬。

圖4 COHESIVE單元B-K準(zhǔn)則圖解[21]Fig.4 Schematic of B-K criterion used for COHESIVE elements[21]

表6 粘結(jié)層的力學(xué)性能[18]Table 6 Mechanical properties of cohesive layer[18]

3.2 結(jié)構(gòu)參數(shù)

T型構(gòu)件由兩個(gè)L型腹板、底板、中心富樹脂區(qū)3部分構(gòu)成,各部尺寸詳見圖5。底板3的厚度為4 mm,含32層鋪層,L型腹板1和2的厚度都為1.5 mm,由13層單層板制成。各部分鋪層角詳見表7。

圖5 T型構(gòu)件結(jié)構(gòu)圖[18]Fig.5 Structural diagram for T-joint[18]

表7 T型構(gòu)件鋪層參數(shù)[18]Table 7 Lamination parameters of T-joint[18]

3.3 有限元建模

結(jié)構(gòu)示意如圖6,在離底板下底面80 mm處施加500 N的水平力,在底部?jī)蓚?cè)距中心90 mm處設(shè)置固支邊界。

層板用S4R殼單元,中央三角富樹脂區(qū)用C3D8R實(shí)體單元。在曲率巨變或應(yīng)力集中區(qū)域進(jìn)行網(wǎng)格加密。表7中的各子結(jié)構(gòu)彼此接觸面上插入COH3D8 Cohesive單元。

圖6 T型構(gòu)件有限元模型Fig.6 FE modelling for T-joint

3.4 結(jié)果與討論

圖7為T型結(jié)構(gòu)受載343 N時(shí)的位移云圖。達(dá)到182 N時(shí),右側(cè)L型層合板在與富樹脂層連接的部分單元發(fā)生基體的壓縮破壞,程序?qū)Πl(fā)生非致命破壞的單元進(jìn)行剛度折減;當(dāng)荷載為343 N時(shí),由于局部應(yīng)變過大,左側(cè)L型腹板與底板連接區(qū)域發(fā)生致命破壞,層合板單元(圖8)的應(yīng)變超過臨界值,此時(shí)計(jì)算終止。

圖7 T型接頭破壞時(shí)的位移云圖Fig.7 Displacement contour of T-joint under failure load

模擬的T型接頭在加載點(diǎn)處的荷載-位移曲線如圖9所示,極限荷載為343 N,而原文獻(xiàn)中的實(shí)驗(yàn)結(jié)果為288.5 N,兩者相差54.5 N,預(yù)報(bào)誤差為18.9%,低于原模擬結(jié)果19.8%的誤差率。由于極限破壞由臨界應(yīng)變條件控制,進(jìn)一步降低臨界應(yīng)變至5%,預(yù)報(bào)的極限載荷為307.2 N,與實(shí)驗(yàn)對(duì)比誤差更小。因?yàn)楸纠欣w維和基體皆為脆性材料,取5%的臨界應(yīng)變似更合理。這說明,將橋聯(lián)理論通過UGENS用于復(fù)雜復(fù)合材料結(jié)構(gòu)的破壞與強(qiáng)度模擬十分有效,剛度與實(shí)驗(yàn)結(jié)果相比略小,相對(duì)誤差為8.21%。原因可能是取自其他文獻(xiàn)的組分材料彈性性能并非與實(shí)際的完全一致。不出所料,若將粘結(jié)層的強(qiáng)度參數(shù)取自文獻(xiàn)[18],預(yù)報(bào)的復(fù)合材料強(qiáng)度會(huì)明顯降低,見圖9,這也表明,在接觸面引入Cohesive單元是恰當(dāng)?shù)摹H舾鹘Y(jié)構(gòu)通過粘接而成,粘接面接觸不良以及非一致性,將導(dǎo)致粘接強(qiáng)度明顯降低,對(duì)復(fù)合材料結(jié)構(gòu)的承載能力也會(huì)產(chǎn)生顯著影響。

圖8 T型接頭破壞位置及形式Fig.8 Failure modes and locations of T-joint

圖9 載荷-位移曲線對(duì)比Fig.9 Comparison of load-displacement curves

3.5 T-型接頭設(shè)計(jì)

3.5.1 不同組分材料的影響

由于發(fā)生致命破壞的是L型腹板,只需對(duì)該構(gòu)件進(jìn)一步考慮。首先,比較3種不同組分材料體系,其力學(xué)性能參數(shù)如表8所示。

將L型腹板用上述3種材料替換,底板性能依然由表5提供,由此得到的T-型接頭載荷-位移曲線如圖10所示,不同材料體系(用于L腹板)的極限破壞模式如表9所示。

表8 3種復(fù)合材料的組分材料性能參數(shù)Table 8 Constituent properties of three composites

圖10 4種組分材料載荷-位移曲線對(duì)比Fig.10 Comparison of load-sdisplacement curves of T-joints made of 4 constituent systems

表9 4種組分材料體系的極限破壞模式Table 9 Ultimate failure modes of four constituent systems

圖10表明,兩種碳纖維材料體系預(yù)報(bào)的結(jié)構(gòu)剛度相同,但強(qiáng)度明顯有異,這是由于AS4纖維的拉壓強(qiáng)度遠(yuǎn)低于T700纖維的水平。而兩種玻璃纖維/環(huán)氧樹脂材料預(yù)報(bào)的模量與強(qiáng)度均相近,源于兩者采用了相同性能的纖維。綜合而言,碳纖維樹脂基體組合比玻璃纖維/樹脂基體組合有更高的模量與強(qiáng)度,而原結(jié)構(gòu)所用的T700/QY9611材料腹板相比其他3種材料有更強(qiáng)的抗彎曲載荷能力。

3.5.2 不同鋪層角

下面考慮L型腹板不同鋪層角對(duì)結(jié)構(gòu)承載能力的影響。假定層數(shù)、厚度與材料性能不變,設(shè)計(jì)了3種不同鋪層角,如表10所示。

依然取材料T700/QY9611,結(jié)構(gòu)致命破壞模式如表11所示,預(yù)報(bào)的不同腹板鋪層T型結(jié)構(gòu)載荷-位移曲線如圖11所示。

從圖11中可以看出,4種不同鋪層的T型接頭在彎曲載荷下的力學(xué)響應(yīng)差異很大,說明鋪層方案對(duì)復(fù)合材料結(jié)構(gòu)優(yōu)化影響顯著。按剛度大小排序依次是:鋪層3>鋪層1>原鋪層>鋪層2;按極限承載能力排序依次是:鋪層3>鋪層1>原鋪層>鋪層2。由此可見,鋪層3的抗彎能力較原鋪層的抗彎載荷與剛度都有顯著提高,為4種鋪層設(shè)計(jì)中的最優(yōu)方案。

表10 T型接頭中3種L型腹板鋪層參數(shù)Table 10 Parameters of three lay-ups of L-web in T-joint

表11 4種鋪層結(jié)構(gòu)極限破壞模式Table 11 Ultimate failure modes of four lay-ups

圖11 4種鋪層的載荷-位移曲線對(duì)比Fig.11 Comparison of load-displacement curves with 4 kinds of lay-ups

4 結(jié) 論

本文將最新發(fā)展的橋聯(lián)理論與有限元軟件ABAQUS結(jié)合,通過編寫UGENS殼單元子程序,實(shí)現(xiàn)了對(duì)復(fù)合材料層合板結(jié)構(gòu)在復(fù)雜載荷下的漸進(jìn)破壞模擬和強(qiáng)度預(yù)報(bào)。通過對(duì)WWFE中14道考題的模擬預(yù)報(bào),其結(jié)果相較實(shí)驗(yàn)的精確度提高到86.4%,超過該競(jìng)賽評(píng)分最高的Zinoviev理論(76.8%)。在此基礎(chǔ)上,對(duì)復(fù)合材料復(fù)雜結(jié)構(gòu)——T型接頭進(jìn)行了有限元極限分析,結(jié)果與有關(guān)文獻(xiàn)中的測(cè)試報(bào)告吻合良好(誤差8.21%),并進(jìn)一步對(duì)T型接頭中腹板的組分材料與鋪層方案進(jìn)行了優(yōu)化設(shè)計(jì)。應(yīng)用本文介紹的程序和方法,可以對(duì)更復(fù)雜復(fù)合材料或復(fù)合材料與其他材料組合結(jié)構(gòu)進(jìn)行破壞分析、強(qiáng)度預(yù)報(bào)和設(shè)計(jì)選型。

z=x2+ix3=reiφ,N2=aFe-iψ+akeiψ

N3=Faeiψ+e-iψak

H=a(cosψ+2λsinψ)(0.5-F)

C=(k-1)(cosψ-2λsinψ)a2exp[2λ(ψ-π)]

D=(1-k)a3exp[2λ(ψ-π)]

式中:a為纖維的半徑,可取單位值;開裂角ψ由下述方程求解:

R(eiφ)=(eiφ-eiψ)0.5+iλ(eiφ-e-iψ)0.5-iλe-iφ

H0=(1-k)(1+4λ2)sin2ψ

J1=kG0-1-2(1-k)ξexp(2λψ)cosψ

J2=2(1-k)ξexp(2λψ)sinψ

J3=2(1-k)ξ·

exp(2λψ)[J1cosψ-J2sinψ]/J2

需要指出的是,ξ=1對(duì)應(yīng)奇異開裂,此時(shí)開裂角ψ無解。然而,總可以微調(diào)纖維/基體性能(其測(cè)試存在離散性)使得ξ≠1。

基體的臨界von Mises等效應(yīng)力計(jì)算公式為

猜你喜歡
復(fù)合材料有限元
金屬復(fù)合材料在機(jī)械制造中的應(yīng)用研究
纖維素基多孔相變復(fù)合材料研究
新型有機(jī)玻璃在站臺(tái)門的應(yīng)用及有限元分析
基于有限元的深孔鏜削仿真及分析
基于有限元模型對(duì)踝模擬扭傷機(jī)制的探討
民機(jī)復(fù)合材料的適航鑒定
復(fù)合材料無損檢測(cè)探討
磨削淬硬殘余應(yīng)力的有限元分析
TiO2/ACF復(fù)合材料的制備及表征
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 久久香蕉国产线看观看式| 青青热久免费精品视频6| 欧美日本在线| Jizz国产色系免费| 久久五月视频| 亚洲乱亚洲乱妇24p| 午夜视频www| 日韩欧美91| 无码有码中文字幕| 国产精品午夜福利麻豆| 无码AV日韩一二三区| 日本91视频| 亚洲中文字幕无码mv| 毛片网站在线看| 国产三级韩国三级理| 真实国产乱子伦视频| 色婷婷久久| 成人福利在线观看| 国产亚洲精品自在线| 久久96热在精品国产高清| 成人亚洲国产| 久久国产精品夜色| 国产aⅴ无码专区亚洲av综合网| 免费高清自慰一区二区三区| 三级毛片在线播放| 国产真实乱子伦视频播放| 亚洲精品视频在线观看视频| 国产视频你懂得| 国产精品内射视频| 亚洲欧美精品日韩欧美| 在线观看亚洲成人| 亚洲首页国产精品丝袜| 欧美成人h精品网站| 在线免费看片a| 伊人激情综合网| 国产成人亚洲综合A∨在线播放| 成年片色大黄全免费网站久久| 天天躁狠狠躁| 91青青视频| 午夜日b视频| 黄色网站在线观看无码| 久久免费成人| 在线观看免费黄色网址| 天天色天天操综合网| 久久亚洲精少妇毛片午夜无码| 黄色网站不卡无码| 国产精品99r8在线观看| 亚洲成肉网| av一区二区无码在线| 波多野结衣一区二区三视频 | 99久久精品视香蕉蕉| 一级成人a毛片免费播放| 国产三区二区| 麻豆AV网站免费进入| 狠狠色噜噜狠狠狠狠色综合久| 毛片久久网站小视频| 2020国产在线视精品在| 夜精品a一区二区三区| 精品人妻AV区| 日韩欧美国产综合| 欧美性色综合网| 久久特级毛片| 国产精品久久久久久久久| 亚洲av无码片一区二区三区| 国产亚洲一区二区三区在线| 日韩精品无码免费专网站| 久久久受www免费人成| 美女无遮挡免费视频网站| 亚洲无线视频| 日本三级欧美三级| 亚洲天堂日韩在线| 性欧美在线| 老司机aⅴ在线精品导航| 久久久精品久久久久三级| 91 九色视频丝袜| 国内精品久久久久久久久久影视| 久久这里只有精品23| 国产99视频在线| 欧美国产日韩在线播放| 毛片视频网址| 国产尤物视频在线| 在线观看国产网址你懂的|