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

丙炔基雙酚A醚硼聚合物熱解過程的ReaxFF分子動力學模擬

2015-08-20 07:32:44高寧王一超劉育紅
化工學報 2015年4期
關鍵詞:優化結構模型

高寧,王一超,劉育紅

(西安交通大學化學工程與技術學院,陜西 西安 710049)

引 言

隨著科學技術的發展,耐高溫聚合物在交通、能源、建筑和航空航天等領域中的應用越來越廣泛[1-3]。傳統的耐高溫聚合物普遍存在熔融溫度高,成型工藝差等缺點。因此,設計新型具有良好工藝性能和優異耐熱性能的聚合物成為了研究人員關注的焦點[4-6]。從聚合物結構特性分析,分子鏈中含有炔基的樹脂因其固化前黏度小、固化過程無小分子產生、高溫殘炭率高等優點備受關注[7-8]。此外,在聚合物分子鏈中引入無機元素硼,因炭化過程中硼元素的促石墨化作用及高溫下B2O3保護層的形成,可明顯改善材料的耐高溫氧化性能[9-11]。基于上述兩方面的考慮,設計并制備了新型硼雜化丙炔基樹脂——丙炔基雙酚A 醚硼聚合物(PB)[12],其結構式如圖1所示。值得注意的是,PB 固化物的玻璃化轉變溫度(Tg)為416℃,在氮氣氛中800℃時的質量保持率可高達74%。PB 樹脂上述優異的耐熱性能為其作為高性能熱固性樹脂的應用提供了可能。

圖1 PB 的結構式Fig.1 Schematic structure of PB

PB 樹脂在高溫下發生熱分解,釋放出小分子氣體物質,結構從有機物逐步轉變為無定形態無機物。熱分解是由有機態向無機態轉化的重要步驟,也是影響熱解生成物結構的關鍵因素。因此,對PB 的熱解過程及其生成物結構進行研究十分必要。目前,采用熱重-紅外光譜聯用(TG-FTIR)法研究熱解過程較為普遍[13-14]。這一方法可以表征熱解生成的氣體種類,從而推斷出可能發生的反應。但由于熱解過程涉及的物理和化學變化較為復雜,TG-FTIR 法只局限于對熱解過程初步描述,難以準確探究PB樹脂在熱解過程中發生的結構變化。

分子模擬技術可以通過計算機在分子水平上呈現分子的構象變化并計算材料的性能,彌補實驗手段的不足,被稱為第3 種認知自然的方法。其中,反應分子動力學模擬(ReaxFF-MD)因其可模擬聚合物的成鍵、斷鍵過程而受到廣泛的關注。尤其在描述極端條件下的反應現象以及復雜的反應過程時,ReaxFF-MD 法可清晰呈現出物質的結構變化,揭示其反應機理。Jiang 等[15]采用ReaxFF-MD 法計算得出酚醛樹脂的主要熱解產物H2O 的反應活化能為(144±28)kJ·mol-1。這一數值與實驗值相符,證明ReaxFF-MD 在模擬熱解過程上是適用的。并且,通過ReaxFF-MD 法,Jiang 等提出H2O 的生成機理為β-H 的脫除,為酚醛樹脂的熱解理論提供了新的見解。由此可見,ReaxFF-MD 在描述聚合物的熱解過程,揭示其熱解機理等方面具有實驗分析方法不可比擬的優勢[16-17]。

本文在實驗分析的基礎上,采用ReaxFF-MD模擬方法,從微觀角度出發,系統研究PB 耐高溫樹脂的熱解過程。為探究聚合物的熱解過程及聚合物熱穩定性的成因提供新的借鑒。

1 實驗材料和方法

1.1 分析儀器

本實驗采用德國 NETZSCH 公司的型號為STA449-F3 的熱重分析儀及德國 Bruker 公司的Vertex70 型紅外光譜儀,升溫速率為10℃·min-1,氮氣氛流量40 ml·min-1。熱重分析(TGA)采用美國TA 公司的TA2950 熱重分析儀測定。測試條件:升溫范圍,室溫~800℃;升溫速率,10℃·min-1;氮氣氛流量為60 ml·min-1。

1.2 模型構建與優化

1.2.1 分子模型的建立[18-19]

(1)建立PB 分子鏈的結構模型。對于聚合度的選擇,要考慮軟件能夠承受的計算量以及模型是否合理。綜合考慮計算量及精確性,還要兼顧合理性與時效性,最終設定單分子鏈模型的聚合度為3,其模型結構如圖2(a)所示。可以看出,此時模型結構處于較高的能量狀態,而且遠離平衡狀態,因此需要對此結構進行松弛和優化。常用的優化方法主要有分子力學和分子動力學。首先,采用智能最小化的方法對此模型進行能量最小化以獲得合理的三維結構,如圖2(a)所示,分子力場為Dreiding,當能量收斂至0.0010 kcal·mol-1(1 kcal·mol-1=4.18 kJ·mol-1)時即停止優化。然后,將得到的分子模型在300 K 下進行20 ps 的NVT 動力學平衡,完全優化后的模型如圖2(b)所示。通過能量優化前后的對比圖片可以看出分子模型經過優化后結構會變得更合理,同時模型的能量也降低。

(2)定義一個三維盒子。在周期性邊界條件下,基于Dreiding 力場,對聚合物建立了三維周期邊界盒子,盒子尺寸為2.514 nm×2.514 nm×2.514 nm,其中放入8 條硼聚合物分子鏈。將聚合物的初始密度設定為0.88 g·cm-3(實驗值),其結構模型如圖3(a)所示。

圖2 能量優化前和能量優化之后的PB 分子鏈模型Fig.2 Snapshot of simulation of PB chain before and after equilibration

(3)再次采用智能最小化的方法對體系進行能量優化。然而,智能最小化方法采用的是分子力學方法,嚴格地講是在絕對零度下進行的,當體系能量較小時,通常能夠較為準確地優化出分子的最低能量從而得到分子的最優幾何構型。當分子體系較為龐大時,分子力學優化的結構可能只是體系的某個位置。為了避免這種情況,需要在較高的溫度下,對體系進行分子動力學優化,即先將分子“加熱”到較高溫度,讓分子結構在高溫下充分舒展和松弛,然后進行“冷卻”。因此,對模型進行退火模擬,令體系從300 K 為起點,步長為50 K,在101325 Pa下對每次的溫度進行5 ps 的NVT 動力學模擬。直至溫度達到700 K,接著再以50 K 為步長,將體系降至300 K,共重復進行10 次。至此,體系的空間結構得到了逐步優化。最后,在300 K 下進行40 ps的NVT 的分子動力學模擬,當體系的能量基本維持恒定時停止優化。此時,得到了如圖3(b)所示的完全優化之后的周期性邊界盒子。

圖3 在300 K 下能量優化前和能量優化之后的PB 周期邊界盒子(灰色:碳;紅色:氧;白色:氫;黃色:硼)Fig.3 Snapshot of simulation cell equilibrated at 300 K before equilibration and after equilibration(grey:C; red:O; white:H; yellow:B)

1.2.2 ReaxFF-MD 模擬 將優化完全的聚合物模型體系在300 K 下進行40 ps 的NVT 分子動力學平 衡計算后,分別以50、100 和200 K·ps-1的速率升溫至4300 K,通過觀察聚合物在熱解過程中發生的反應、小分子逸出的數量來研究熱解過程并探究其機理。積分算法采用Velocity Verlet 法,控溫方法為與熱浴耦合的方法(Berendsen 法),模擬步長為0.2 fs。

2 實驗結果與討論

2.1 ReaxFF 升溫過程模擬

首先探究不同升溫速率對PB 固化物熱解開始時間的影響。在此,分別模擬了PB 固化物在50、100 和200 K·ps-1升溫速率下聚合物的熱解過程。不同升溫速率下聚合物的個數隨模擬時間的變化如圖4所示。不同升溫速率下熱解開始的時間和溫度列于表1中。由表1可以看出,當升溫速率從50 K·ps-1升高至200 K·ps-1時,熱解反應的開始時間從14.0 ps 縮短到4.8 ps;熱解開始的溫度由700 K 升高到960 K。較高的升溫速率縮短了熱解反應的起始時間,促進了熱解的進行。升溫速率越高,熱解開始的溫度也越高。

表1 不同升溫速率硼聚合物熱解開始的時間和溫度Table 1 Pyrolysis time and temperature as function of heating rate

通過測試不同升溫速率下的DTG(圖5)曲線可知,隨著升溫速率的增加,熱解開始的溫度也升高。這是由于升溫較慢時,體系有充足的反應時間,即在低升溫速率下經歷的熱解時間增加,熱解開始的溫度較低。由此可見,ReaxFF 反應動力學模擬結果與實驗結論一致,進一步證明了ReaxFF 反應動力學在模擬硼聚合物熱解過程是適用的。

2.2 TG-FTIR 法表征PB 固化物的熱解過程

為了探究PB 固化物的高溫熱解過程,首先采用TG-FTIR 對熱分解過程中逸出的尾氣進行在線連續監測,結果見圖6。在整個熱分解過程,3964~3500 cm-1、1854~1500 cm-1處對應于水蒸氣的特征峰[14]及2360~2320 cm-1處屬于CO2的特征峰一直存在[13]。其中,水可能以自由水、結合水的形式存在,還可能在熱解反應的過程中生成。300℃之前裂解氣中的水峰可能由產物吸潮而吸收的自由水或 結合水揮發導致的,300℃之后裂解氣中的水應該是熱裂解反應生成的水分。從400℃開始,熱解尾氣的組成與之前有較大差別。從400℃的裂解氣的紅外譜圖中可以看到,裂解產物中除了水蒸氣和CO2,在3000~2600 cm-1和1500~1325 cm-1處出現了 甲烷的特征峰[20],在1850~1600 cm-1處出現了的雙鍵的振動峰,這對應于醛類、酮類及羧酸類等大分子物質[14]。此外,在1845~1500 cm-1處出現了對應于烯烴的特征峰,說明此時產物中的苯環已經開始發生降解[21]。600℃之后,烷烴和烯烴的紅外特征峰開始逐漸變小最終消失。在2250~2000 cm-1處代表CO 的特征峰[22],CO 的生成對應體系中的碳原子逐漸芳香化的過程[14,23]。

圖4 在50、100 和200 K·ps-1 升溫速率下的硼聚合物(C63H60O12B3)熱解過程的碎片分析Fig.4 Average number of fragments as function of heating rate at 50,100,and 200 K·ps-1

圖5 不同升溫速率時固化后PB 的DTG 曲線Fig.5 Derivative thermogravimetry curves of PB thermosets under different heating rates

圖6 PB 固化物的TG-FTIR 譜圖Fig.6 TG-FTIR spectra of cured PB

TG-FTIR 測試可以得出熱解過程生成的小分子種類,但無法探究熱解過程中PB 的結構變化。通過ReaxFF-MD 法,不僅可以驗證TG-FTIR 測試結果的可靠性,而且可獲得PB 在熱解過程中的斷鍵順序及結構變化。在此,以其中一個結構單元為研究對象,通過觀察升溫反應過程模擬中聚合物結構的變化有助于研究聚合物的熱解機理。首先,以100 K·ps-1的升溫速率,將聚合物體系升至4300 K。由圖4(b)可知,大約在7.3 ps(730 K)時,聚合 物開始降解。在降解前,可以觀察到聚合物鏈的數量并非一成不變,這是由高溫下聚合物鏈較強的活動性,不斷地與相鄰的分子鏈成鍵、斷鍵導致的。圖7描述了PB 固化物在100 K·ps-1熱解過程中斷鍵過程。隨著溫度的升高,大約在12 ps(1300 K)時,開始斷裂的是溴丙炔基雙酚A 醚的脂肪醚鍵,這與文獻[24]相一致;隨后,約在20 ps(2300 K) 時,結構單元中的另外一個脂肪醚鍵也隨之斷裂,雙酚A結構及雙酚A氧自由基結構在模擬體系中大量存在。此外,可以觀察到發生芳基醚鍵斷裂的反應較少,這是由于五電子的芳環不穩定;隨著溫度的進一步升高,約在28 ps(2800 K)時,雙酚A結構中的甲基脫除,此時有大量的甲基自由基[28 ps(b)]生成,甲基自由基可與氫結合生成裂解的主要產物CH4;當溫度進一步升高時,苯環開始斷裂,生成烯烴[31 ps(b)]和酮[36 ps(b)]等有機小分子。

圖7 PB 熱解過程不同時刻分子片段快照(灰色:碳;紅色:氧;白色:氫;黃色:硼)Fig.7 Snapshot of simulated pyolysis of PB at different time(grey :C; red:O; white:H; yellow:B)

可見,ReaxFF-MD 模擬得到的硼聚合物的熱解過程與TG-FTIR 的分析結果一致,均是由脂肪醚鍵的斷裂開始熱解反應,進而脫除甲基、生成甲烷小分子,隨后苯環開始裂解生成烯烴、羧酸類、酮類等有機化合物。因此,ReaxFF-MD 模擬熱解過程,不僅驗證了TG-FTIR 的結論,而且能夠更為直觀、從分子水平上“觀察”聚合物的結構變化。

2.3 ReaxFF 模擬PB 固化物熱解產物的生成機理

通過ReaxFF 反應動力學模擬PB 熱解過程,主要產物有CH4、H2O、CO 和H2(圖8),前3 種產物的在TG-FTIR 分析中得到了證明。其中,熱解生成的 H2是非極性分子,吸收峰很弱,所以在TG-FTIR 譜圖中無法體現。此外,TG-FTIR 譜圖中顯示熱解過程中有CO2分子生成,但在模擬過程中無CO2生成。導致這種差異性的原因是在熱解過程中模擬體系里沒有O2分子,且聚合物結構中沒有羧基、酯基的結構[25]。

圖8 PB 固化物熱解過程生成的小分子CH3?、CH4、CO、H2 和H2O 的分布Fig.8 Average number of fragments of as function of time during pyrolysis of PB thermosets:CH3?,CH4,CO,H2 and H2O

通過追蹤產物CH4的生成過程,可得到4 種CH4的主要生成途徑(RPW m1—4)(圖9)。由此圖所列反應式可知,CH4主要是由CH3?和裂解過程生成的碎片結構中的氫原子結合而成,此結論驗證了上述推測。其中,氫原子主要來源于羥基[包括雙酚A 結構中的羥基(RPW m1)和活潑氫(雙酚A結構中苯環的鄰位氫(RPW m2)、裂解產物重排生成的苯乙烯的α-H(RPW m3)、裂解產物重排生成的炔烴的α-H(RPW m4)],而CH3?主要來源于雙酚A 結構中甲基的脫離。

H2O 分子的3 種主要的生成機理(RPW w1—3)可參考圖10。由圖10可見,雙酚A 結構中的—OH與H?之間的結合以及—OH 之間的結合是生成H2O分子的主要途徑(RPW w2 和RPW w3)。由雙酚A中的—OH 生成H2O 的反應機理在環氧樹脂中也有提及[26]。此外,H2O 的生成還存在分子間的反應,即與硼原子相連的兩個—OH 反應脫去水分子。在熱解過程中,頻繁出現C—B—O 的三元環,RPW w1同樣也是此三元環的一種主要的生成途徑。

圖11列出了模擬熱解過程中H2的5 種主要的生成途徑(RPW h1—5)。從圖11中可以觀察到,H?基團與碎片中的活潑氫發生撞擊是生成H2的主要機理(RPW h1、RPW h2、RPW h5)。RPW h3 為雙酚A 結構—OH 中的氫與甲基中的氫結合生成分子H2。RPW h4 為分子內反應,即甲基中的兩個氫原子形成H2,但在整個熱解階段,此類型的反應出現的頻率較低。

同時,通過對CO 分子產生過程的追蹤,可以獲得CO 分子的5 種主要的生成機理(RPW c1—5)(圖12)。CO 是由側鏈末端羰基或熱解產物中含羰基化合物脫羰基產生的。CO 產物的生成帶走了分子結構中的氧,有助于聚合物的成炭過程,很好地驗證了實驗推測的結果。

圖9 ReaxFF-MD 模擬過程中CH4 的主要生成機理Fig.9 Typical mechanism of CH4 formation observed in ReaxFF MD

圖10 ReaxFF-MD 模擬過程中H2O 的主要生成機理Fig.10 Typical mechanism of H2O formation observed in ReaxFF MD

圖11 ReaxFF-MD 模擬過程中H2 的主要生成機理Fig.11 Typical mechanism of H2 formation observed in ReaxFF MD

圖12 ReaxFF 模擬過程中CO 的主要生成機理Fig.12 Typical mechanism of CO formation observed in ReaxFF MD

3 結 論

本文采用將實驗分析與ReaxFF-MD 相結合方法,研究了PB 固化物的熱解過程并分析了其炭化物的結構。通過觀察升溫過程中PB 固化物結構的變化,可得到其熱解過程中的斷鍵順序以及小分子 生成機理。研究發現硼聚合物熱解過程中首先發生脂肪醚鍵的斷裂,然后雙酚A 中—CH3脫除生成CH4,當溫度進一步升高時,苯環裂解生成烯烴和酮等有機小分子。在整個熱分解過程中,生成的小分子主要有H2、H2O、CO 和CH4。裂解過程中的CH4主要由甲基和裂解過程生成的碎片結構中的活潑氫結合生成,H2主要由H?與碎片中的活潑氫撞擊生成的,雙酚A 結構中的—OH 和氫的結合是生成H2O 分子的主要途徑,而CO 的生成對應裂解的碎片結構中羰基的脫除。模擬結果與實驗結論吻合良好,同時從分子角度描述了PB 固化物熱解過程。

[1]Liu Jingang(劉金剛),Shen Dengxiong(沈登雄),Yang Shiyong(楊士勇).Recent research and development of matrix resins for high temperature polymer composites [J].Aerospace Mater.Technol.(宇航材料工藝),2013,44 (4):8-13

[2]Hirano K,Asami M.Phenolic resins—100 years of progress and their future [J].React.Funct.Polym.,2013,73 (2):256-269

[3]Cheng R,Zhou Q,Ni L,Chen Y,Wang J.Synthesis and thermal property of boron-silicon-acetylene hybrid polymer [J].J.Appl.Polym.Sci.,2011,119(1):47-52

[4]Chung T C M,Jeong Y,Chen Q,Kleinhammes A,Wu Y.Synthesis of microporous boron-substituted carbon (B/C) materials using polymeric precursors for hydrogen physisorption [J].J.Am.Chem.Soc.,2008,130(21):6668-6669

[5]Devapal D,Packirisamy S,Sreejith K J,Ravindran P V,George B K.Synthesis,characterization and ceramic conversion studies of borosiloxane oligomers from phenyltrialkoxysilanes [J].J.Inorg.Org.Polym.Mater.,2010,20(4):666-674

[6]Zhou Q,Mao Z,Ni L,Chen G,Zhan M S.Chemistry and application of dicarbadodecaborane and its derivatives [J].J.Appl.Polym.Sci.,2007,104(4):2498-2503

[7]Zhang Shijie(張世杰),Zhang Wei(張煒),Guo Yalin(郭亞林),Bai Xia(白俠).Research progress in a new ablate resistant material-polyarylacetylene resin [J].Thermosetting Resin(熱固性樹脂),2007,22(6):42-46

[8]Zhang Dexiong(張德雄),Zhang Yan(張衍).An overview of polyarylacetylene as the resin matrix for high temperature composites [J].Journal ofSolid Rocket Technol.(固體火箭技術),2001,24(1):53-59

[9]Liu Y H,Qiang J P,Jing X L.Synthesis and properties of a novel hyperbranched borate [J].J.Polym.Sci.Part A:Polym.Chem.,2007,45(15):3473-3476

[10]Zhong D H,Sano H,Uchiyama Y,Kobayashi K.Effect of low-level boron doping on oxidation behavior of polyimide-derived carbon films [J].Carbon,2000,38(8):1199-1206

[11]Wang J G,Guo Q G,Liu L,Song J R.Study on the microstructural evolution of high temperature adhesives for graphite [J].Carbon,2002,40(13):2447-2452

[12]Gao Ning (高寧).Synthesis,characterization and molecular dynamics simulations of a novel dipropargyl ether of bisphenol A based boron-containing polymer[D].Xi’an:Xi’an Jiaotong University,2014

[13]Atkinson J L,Vyazovkin S.Thermal properties and degradation behavior of linear and branched poly(L-lactide)s and poly(L-lactide-co-glycolide)s [J].Macromol.Chem.Phys.,2012,213(9):924-936

[14]Wang Shurong(王樹榮),Liu Qian(劉倩),Zheng Yun(鄭赟),Wen Lihua(文麗華),Luo Zhongyang(駱仲泱),Cen Kefa(岑可法).Mechanism study of biomass pyrolysis by thermogravimetric analysis coupled with infrared spectroscopy [J].J.Eng.Thermophys(工程熱物理學報),2006,27(2):351-353

[15]Jiang D E,van Duin A,Goddard W A,Dai S.Simulating the initial stage of phenolic resin carbonizationviathe ReaxFF reactive force field [J].J.Phys.Chem.A,2009,113(25):6891-6894

[16]Tack J L,Ford D M.Thermodynamic and mechanical properties of epoxy resin DGEBF crosslinked with DETDA by molecular dynamics [J].J.Mol.Graphics Modell.,2008,26(8):1269-1275

[17]Ho S R,Aoyagi M.Molecular-cluster-assembly method for analysis of high-dimensional structures ofp-cresol chains in phenolic polymers [J].Macromol.Theory Simul.,2013,22(9):443-61

[18]Andrew R.Molecular Modelling Principles and Applications [M].Harlow:Pearson Education Limited,2001:66-70

[19]Tang Bing (唐兵).Conformational statistics and computer simulation of restricted polymer chains [D].Chengdu:Sichuan University,2002

[20]Salema A A,Afzal M T,Motasemi F.Is there synergy between carbonaceous material and biomass during conventional pyrolysis? A TG-FTIR approach [J].J.Anal.Appl.Pyrol.,2014,105(1):217-226

[21]Su S P,Wilkie C A.The thermal degradation of nanocomposites that contain an oligomeric ammonium cation on the clay [J].Polym.Degrad.Stabil.,2004,83(2):347-362

[22]Chen Y F,Chen Z Q,Xiao S Y,Liu H B.A novel thermal degradation mechanism of phenol-formaldehyde type resins [J].Thermochimica Acta,2008,476(1/2):39-43

[23]Maschio G,Koufopanos C,Lucchesi A.Pyrolysis,a promising route for biomass utilization [J].Bioresource Technol.,1992,42(3):219-231

[24]Itoh M,Iwata K,Ishikawa J I,Sukawa H,Kimura H,Okita K.Various silicon-containing polymers with Si(H)―C≡C units [J].J.Polym.Sci.Part A:Polym.Chem.,2001,39(15):2658-2669

[25]Chenoweth K,Cheung S,van Duin A,Goddard W A,Kober E M.Simulations on the thermal decomposition of a poly(dimethylsiloxane) polymer using the ReaxFF reactive force field [J].J.Am.Chem.Soc.,2005,127(19):7192-7202

[26]Diao Z J,Zhao Y M,Chen B,Duan C L,Song S.ReaxFF reactive force field for molecular dynamics simulations of epoxy resin thermal decomposition with model compound [J].J.Anal.Appl.Pyrol.,2013,104(1):618-624

猜你喜歡
優化結構模型
一半模型
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
論《日出》的結構
主站蜘蛛池模板: 永久免费av网站可以直接看的 | 国产一区三区二区中文在线| 玖玖精品在线| 成人在线天堂| 丁香亚洲综合五月天婷婷| 亚洲人成网站色7799在线播放| 91精品国产自产在线观看| 亚洲色图欧美激情| 黄色网站不卡无码| 91口爆吞精国产对白第三集| 91在线一9|永久视频在线| 亚洲无码高清免费视频亚洲| 亚洲中文字幕无码mv| 亚洲国产成熟视频在线多多| 在线精品自拍| 国产成人免费视频精品一区二区| 日韩第一页在线| 精品福利视频导航| 在线看国产精品| 好久久免费视频高清| 亚洲精品欧美日韩在线| 欧美精品一二三区| 国产爽爽视频| a级毛片免费网站| 中日韩欧亚无码视频| 国产不卡网| 欧美精品影院| 四虎精品国产AV二区| 伊人激情久久综合中文字幕| 伊人久久大香线蕉影院| 亚洲伊人久久精品影院| 国产精选自拍| 99性视频| 很黄的网站在线观看| 四虎成人精品在永久免费| 亚洲天堂精品在线| 四虎永久在线| 国产国产人免费视频成18| 性69交片免费看| 美女被躁出白浆视频播放| 99视频精品全国免费品| 国产网友愉拍精品| 亚洲日本中文字幕天堂网| 国产综合精品日本亚洲777| 人妻少妇久久久久久97人妻| 国产拍揄自揄精品视频网站| 国产天天色| 在线观看精品国产入口| 国产午夜在线观看视频| 国产在线观看一区二区三区| 国产95在线 | 99视频有精品视频免费观看| 亚洲精品你懂的| 国产欧美日韩在线一区| www.亚洲天堂| 亚洲a级在线观看| 精品国产免费观看| 亚洲无码A视频在线| 欧美亚洲国产日韩电影在线| 99国产在线视频| 亚洲Av综合日韩精品久久久| 女同久久精品国产99国| 一级毛片高清| 中国精品久久| 成人福利在线免费观看| 青青青草国产| 中文字幕色在线| 亚洲永久色| 午夜激情婷婷| 制服丝袜国产精品| 日本国产一区在线观看| 亚洲欧洲自拍拍偷午夜色| a在线亚洲男人的天堂试看| 91伊人国产| 极品国产在线| 一本大道东京热无码av| 免费人欧美成又黄又爽的视频| 91精品国产一区| 看国产一级毛片| 亚洲精品色AV无码看| 免费中文字幕在在线不卡| 在线五月婷婷|