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

分子動力學模擬在復合固體推進劑研究中的應用

2018-01-04 02:57:08苗瑞珍高曉敏喬小平陳海洋郗明頗
兵器裝備工程學報 2017年12期
關鍵詞:體系研究

苗瑞珍,高曉敏,喬小平,陳海洋,郗明頗,孟 勇

(西安北方惠安化學工業有限公司, 西安 710302)

【彈藥工程】

分子動力學模擬在復合固體推進劑研究中的應用

苗瑞珍,高曉敏,喬小平,陳海洋,郗明頗,孟 勇

(西安北方惠安化學工業有限公司, 西安 710302)

概述了分子動力學模擬方法在復合固體推進劑的組分相容性、感度及安全性研究中的應用和分子動力學在研究固體推進劑中炸藥爆轟及燃燒反應機理方面的應用,為從分子層面解釋炸藥爆炸時的微觀反應機理提供了可能;分析了目前存在的問題及發展前景,可為相關科研人員提供理論參考。

分子動力學模擬;固體推進劑;性能;反應機理

1 前言

復合固體推進劑是一種以高分子粘合劑為基體,添加氧化劑、增塑劑和燃燒劑等填料制成的一種推進劑[1],隨著化學、材料等學科的不斷發展以及火箭技術的不斷進步,對固體推進劑的性能要求越來越高,高能、鈍感和低特征信號成為復合固體推進劑發展的主要目標[2]。性能優異的復合固體推進劑配方的研制,需要耗費大量的人力、物力和財力,不僅研制周期較長,同時實驗還具有一定的危險性。1993年,Cumming等[3]在25屆ICT國際會議上報道了HMX和PNMO相互作用的分子力學(MM)和分子動力學(MD)模擬計算工作,開創了近代計算化學對含能材料分子間相互作用理論研究的先河。美國伊利諾伊大學先進火箭發動機仿真中心(CSAR)自1997年開始就在美國能源部ASCI計劃的資助下進行了固體火箭發動機全系統仿真[4],其中一部分就涉及到復合固體推進劑的模擬研究[5]。近年來,隨著理論化學和計算化學以及計算機模擬技術的迅速發展,越來越多的研究者開始應用計算機模擬技術對固體推進劑的結構、性能和微觀反應機理進行研究以彌補實驗手段和理論研究在這方面的一些不足之處,以期為固體推進劑的配方設計,新材料設計、開發提供前期的理論預測與科學依據。

本文主要介紹了MD模擬方法在復合固體推進劑的性能及其炸藥微觀反應機理研究方面的應用進展,分析應用中存在的問題并展望未來的發展前景,以期為科研人員利用該方法開展固體推進劑研究提供參考。

2 在復合固體推進劑性能研究中的應用

2.1 在復合固體推進劑相容性研究中的應用

為了改善固體推進劑的低溫力學性能,降低其玻璃化轉變溫度,同時增加其柔韌性并優化加工性能,需要在配方中加入功能組分增塑劑,理想的增塑劑必須與粘結劑具有良好的相容性,才能保證推進劑具有優異的綜合性能[6-9]。此外也可以在炸藥中加入一定量相容性較好的高聚物粘結劑來改善其力學性能及加工性能。

根據高分子溶液理論,由粘合劑和增塑劑組成的混合體系可以視為高分子溶液,一般情況下溶解是溶質分子和溶劑分子互相混合作用的過程,兩種物質混合時,在恒溫恒壓下,這種過程能自發進行,形成相容體系的熱力學條件是:

ΔGM= ΔHM-TΔSM<0

(1)

其中:T是溶解時溫度;ΔSM、ΔGM和ΔHM分別代表混合熵,混合自由能和混合焓。

對于高分子體系, 如果異種分子間沒有相互作用(如氫鍵), 那么ΔHM值總是大于零的, 因此,混合焓項始終不利于兩者的混合, 能否均勻混合取決于熵項ΔSM的貢獻是否能克服混合焓項。實際中, 對于推進劑內高分子材料在混合過程中, 熵的增加非常小[10],因此ΔGM值的正負取決于ΔHM的大小。 Hildebrand等[11]通過研究表明物質間的相互作用能力決定于其內聚能密度, 因此引入了溶度參數的概念, 其定義為內聚能密度的平方根:

δ=(ΔE/V)1/2=[(ΔHV-RT)/V]1/2

(2)

ΔE、V、ΔHV分別為體系的內能、體積和蒸發熱。高分子材料混合過程中的焓變與高分子的溶度參數的關系為[20]:

ΔHM/V=(δ1-δ2)2φ1φ2

(3)

φ1、φ2分別為組分1、2 的體積分數. 由(3)式可見, 組分1、2 的溶度參數δ1、δ2越接近, ΔHM值越小, 體系相容性越好。因此, 溶度參數差值(Δδ)可以作為組分間相容性理論預測指標。

表1 采用不同方法所得溶度參數δ (J1/2·cm-3/2)[12]

表2 采用不同方法計算的Δδ(J1/2·cm-3/2) [12]

付一政等[12]用MD模擬方法預測了聚丁二烯(HTPB)和增塑劑葵二酸二辛酯(DOS)以及硝化甘油(NG)之間的相容性。表1為用MD模擬計算所得純物質的溶度參數和實驗值[13-15], 以及理論值(目前現有的從理論上計算高分子溶度參數的主要方法是Dunkel[16]發展的原子與基團貢獻法。表1同時給出了根據原子與基團貢獻法原理采用不同估算方法計算的值,其中δFedors、δvanKrevelen分別為采用Fedors[17]和van Krevelen[18]參數求得的溶度參數)。由表可得, 對于HTPB、DOS、NG 來說δFedors、δVK和δMD均與實驗值吻合得較好, 說明通過MD 模擬可以得到與實驗值比較吻合的溶度參數。

表2給了采用不同方法所得的Δδ值,有研究表明[19]對于分子間不存在強作用(如極性集團或者氫鍵作用)的高分子體系,兩種材料的Δδ只要滿足∣Δδ∣<(1.3-2.1) J1/2·cm-3/2, 表明兩者就可以相容。綜合表1 和2,不同的計算方法得到不同的溶度參數值,但是通過計算Δδ判斷共混物相容性的結果均一致, 即∣Δδ∣HTPB/DOS<(1.3-2.1) J1/2·cm-3/2, HTPB/DOS 共混物屬于相容體系, 而∣Δδ∣HTPB/NG的值較大,與相容性較好的HTPB/DOS體系相比, HTPB與NG的相容性較差,屬于不相容體系, 這一結論與實驗結果一致,由此說明MD模擬可以作為研究推進劑相容性的可靠方法。

李倩[20-21]則對高能推進劑中疊氮粘合劑與硝酸酯進行了 MD 模擬,得出了與實際相符的各組分的溶度參數;Li 等[22]從介觀尺度計算了 NEPE 推進劑粘合劑與增塑劑間的相分離情況;Abou-Rachid采用MD方法研究了固體推進劑的相容性[23];許曉娟,肖繼軍等[24]運用MD方法分別預測了以下5種體系相容性和穩定性大小的順序依次為:ε-CL-20/PEG>ε-CL-20/Estane5703>ε-CL-20/GAP>ε-CL-20/HTPB>ε-CL-20/F2314,與實驗結果一致;焦東明等[25]運用MD方法計算了HTPB粘合劑及增塑劑DOS、DOA、TOA、DBP和DOP的溶度參數,并以此來選擇合適的增塑劑,計算數值基本吻合實驗值;王歡等[26]運用MD方法研究了固體推進劑BTTN、NG、HMX、HTPB、DOS,EC等的相容性,結果表明,BTTN、NG均與EC相容,HTPB與DOS是相容體系,而HTPB/NG、HMX/BTTN和HMX/NG體系均為不相容體系。這些模擬計算為固體推進劑的配方設計與發展提供了強有力的理論指導和技術支持。

2.2 在復合固體推進劑感度及安全性研究中的應用

通常情況下分子中引發鍵的鍵級越小,鍵長越大,即越容易斷裂,引發分解和起爆,相應的其感度就越高。經典模擬過程不涉及鍵級,但可以給出鍵長的統計平均值。肖鶴鳴課題組以及國外的其他研究中曾使用引發鍵最大鍵長(Lmax)研究了含能材料的感度,所得結論相符實驗數據,說明釆用隨來關聯含能材料的感度是可靠的。

此外力學性能密切關聯感度,若推進劑體系的力學強度較小,就使得其剛性較小,柔性增強,即體系變“軟”,在體系受到外力作用時,可以有效緩沖和分散外力作用,減小炸藥顆粒之間的摩擦,使其內部應力分布均勻,從而減少“熱點”的形成,減小其感度,增加安全性。

結合能也是含能材料穩定性及安全性的一個重要參考指標,結合能(Ebind)為相互作用能(Einter)的負值,對于某一個兩組分體系A與B而言,A與B的之間的相互作用能等于平衡結構的總能量(Etotal)減去除掉A組分后該結構的能量(EB),再減去除掉B組分后該結構的能量(EA),即:

Ebind=-Einter=-(Etotal-EA-EB)

結合能是組分相互作用強弱的標志。結合能越大,表示組分間的相互作用越強,形成的體系越穩定,相容性越好,對結合能的貢獻主要是非鍵相互作用引起的,其中靜電作用力所占比例最大,范德華能量項所占比例最小。各組分間若存在強烈的分子間作用力,會增加了炸藥分子體系的力學穩定性以及熱力學穩定性,提高了共晶分子對機械外力的抗振性與耐熱性,也即降低了其撞擊感度以及熱感度,從而提高其安全性。

劉強等[27]運用MD方法計算了CL-20/TNT共晶不同晶面以及RDX/TNT共混體系的最大引發鍵鍵長,力學性能以及結合能與感度和其體系安全性的關系,由計算結果可知CL-20/TNT共晶體系中(100)晶面Lmax最大,拉伸模量、剪切模量和體積模量的排序依次為(3×2×1)>(001)>(010)>(100)>(120),而對于RDX/TNT共混體系來說,混合體系中RDX組分的引發鍵(N-O2)的Lmax均比純的RDX晶體短,說明TNT的加入會使體系感度下降,隨著體系中TNT組分含量的增加,其拉伸模量和剪切模量逐漸增大,體系剛性增加,感度增大,此結論與實驗結果能夠很好的吻合;肖鶴鳴課題組近年來采用MD方法模擬了多種含能材料及其復合物的結構性能與感度的關系并取得一定進展,豐富了高能復合物配方研制理論體系;付一政等人[28]采用MD模擬方法對CL-20、DNB、兩者的共混物及共晶的感度、結合能和力學性能進行了模擬計算,結果表明共晶和共混均會降低體系的感度,但共晶效果更加明顯;與共混物結構相比,共晶結構更加穩定,共晶和共混均可以改變復合體系的力學性能,降低體系的剛度,增加體系的柔性和安全性,但共混會使體系的力學性能劣化,與實驗結果一致。感度是衡量推進劑安全性的重要參數,MD模擬為感度研究提供了一種更簡單有效的方法。

2.3 在復合固體推進劑爆轟及燃燒反應機理中的應用

基于反應力場的MD模擬方法是連接量子化學和分子動力學方法的橋梁,為大規模研究凝聚態性質及處理其中可能發生的化學反應過程提供了可能。近幾年反應力場在材料科學尤其對于固體推進劑中炸藥的爆轟及燃燒反應的研究中的應用越來越廣泛,不僅可以對新型炸藥的設計和使用過程提供理論依據,而且可以對分子原子級別的微觀物理化學現象進行解釋。其中用于研究反應機理的力場現在運用最多的是ReaxFF力場,它已經成功用于一些固體推進劑中炸藥的爆轟及燃燒反應的機理研究中。如三硝基甲苯(TNT),環三亞甲基三硝銨(RDX),季戊四醇四硝酸酯 (PETN),環四亞甲基四硝銨(HMX),1,3,5三胺基-2,4,5-三硝基苯(TATB),三過氧化三丙酮(TATP),硝基甲烷(NM)的沖擊起爆、沖擊點火和爆轟反應等,這不僅對含能材料的理論設計進行指導,而且對材料的使用、存儲、運輸安全問題都有重要的現實意義。

硝基甲烷其實并不是實際應用的炸藥,但是由于含有NO2基團,可以代表一類含有硝基的炸藥,鑒于其簡單的分子結構,且在沖擊作用下可以發生爆轟現象所以被廣泛用于炸藥反應的模擬中。Naomi Rom等[29]研究了壓力對硝基甲烷分解路徑的影響,發現在高壓和低壓下熱解引發反應存在一個競爭。Han等[30]研究了高壓下固相硝基甲烷的分解機理,發現在較低的溫度2000K和2500K時單個硝基甲烷分子會發生分子內的重排,生成亞硝酸酯(CH3ONO),接著發生C-N鍵的斷裂,而在較高的溫度3000K時兩個硝基甲烷分子間會發生質子轉移形成C3HNOOH和 CH2NO2,說明溫度會對反應路徑有較大的影響。Strachan等[31]研究了密度對環三亞甲基三硝胺(RDX)晶體熱解過程的影響,得到了隨著密度的變化各產物數量的變化規律以及低密度與高密度晶體分解機理的主要區別。Strachan等[32]研究了不同沖擊速度下RDX的分解反應,結果表明沖擊速度對RDX熱解反應機理有較大的影響。Michael F. Russo等[33]研究了不同當量比的1,5-二硝基縮二脲(DNB)與硝酸的燃燒反應,研究表明在一定的當量比時,會發生自燃反應,釋放出大量的熱能。M.R. Weismiller等[34]運用ReaxFF力場對硼烷氨氧化反應的動力學機理進行了研究。Adri C. T. van Duin等[35]研究了硼烷氨的脫氫反應和燃燒反應的機理,文中首先是運用量子力學計算得到的參數對硼烷氨的力場參數進行擬合,然后對硼烷氨單分子和多分子的熱解和燃燒反應進行了分子動力學模擬,得到產物第一分子H2的生成是分子內部的反應,反應的活化能是26.36kcal/mol,這與實驗值117 kJ/mol很接近[36]。這些模擬計算為固體推進劑的配方設計與發展提供了強有力的理論指導和技術支持。

3 結論

通過對已有的火炸藥品種進行分子動力學模擬研究,不僅可以彌補宏觀實驗手段研究尺度的不足,而且有助于認識火炸藥結構與性能之間的內在關系和作用機理,豐富火炸藥基礎理論體系,為其配方設計和工藝參數的選擇提供參考。

采用MD模擬方法對新型火炸藥品種的進行配方設計和性能預測,能夠大大縮短研制周期,降低成本,提高研發過程中的安全性。

由于一般推進劑中火炸藥的反應條件比較極端,反應速度非常快,用常規的實驗手段無法捕捉一些微觀信息,目前對火炸藥的分子動力學研究大部分僅僅是模擬仿真,缺少實驗驗證。由于計算機性能及模擬效率等因素的限制,模擬體系選取的原子數量有限,模擬規模較小,以及目前力場的不完善,仿真模型與實際情況存在較大的差距。

MD模擬方法在目前的應用中存在一定的局限性,但通過理論與實驗研究緊密結合,MD模擬方法必將成為研制高性能復合固體推進不可或缺的研究手段。

[1] 任務正, 王澤山.火炸藥理論與實踐 [M].北京: 中國北方化學工業總公司, 2001: 120-135.

[2] 張曉宏, 莫紅軍.下一代戰術導彈固體推進劑研究進展 [J].火炸藥學報, 2007, 30(1): 24-27.

[3] CUMMING A, LEIPER G, ROBSON E.Molecular modelling as a tool to aid the design of polymer bonded explosives [J].Journal of Hazardous Materials 1993, 3: 23-26.

[4] DICKFL W A, FIEDLER R A, HEATH M T.Integrated simulation of solid propellant rockets [J].2000, 11: 21-24.

[5] DICK W A, HEATH M T, FIEDLER R A, et al.Advanced Simulation of Solid Propellant Rockets from First Principles [J].Journal of Hazardous Materials, 2005, 4: 10-13.

[6] EDGAR K J, BUCHANAN C M, DEBENHAM J S, et al.Advances in cellulose ester performance and application[J].Progress in Polymer Science, 2001, 26(9):1605-1688.

[7] MUTHIAH R, SOMASUNDARAN U I, VERGHESE T L, et al.Energetics and Compatibility of Plasticizers in Composite Solid Propellants[J].Defence Science Journal, 2013, 39(2).

[8] MATHIEU J, STUCKI H.Military high explosives [J].CHIMIA International Journal for Chemistry, 2004, 58(6): 383-389.

[9] JAWALKAR S N, MEHILAL, RAMESH K, et al.Studies on the effect of plasticiser and addition of toluene diisocyanate at different temperatures in composite propellant formulations.[J].Journal of Hazardous Materials, 2009, 164(2):549-554.

[10] SUN X Q, FAN X W, XUE-HAI J U, et al.Research Methods on Component Compatibility of Propellants[J].Chemical Propellants & Polymeric Materials, 2007, 152(4), 362-369.

[11] HILDEBRAND J H,SCOTT R L.The solubility of non-electrodytes[M].NewYork: Reinhold Publishing Corp., 1950: 424-434.

[12] 付一政, 劉亞青, 蘭艷花.端羥基聚丁二烯/增塑劑共混物相容性的分子動力學模擬[J].物理化學學報, 2009, 25(7):1267-1272.

[13] 楊月誠, 焦東明, 強洪夫,等.HTPB推進劑組分溶度參數的分子模擬研究[J].含能材料, 2008, 16(2):191-195.

[14] ABOU-RACHID H, LUSSIER L S, RINGUETTE S, et al.On the Correlation between Miscibility and Solubility Properties of Energetic Plasticizers/Polymer Blends: Modeling and Simulation Studies[J].Propellants Explosives Pyrotechnics, 2008, 33(4):301-310.

[15] 任玉立, 陳少鎮.關于硝化纖維素濃溶液的研究——體系溶度參數與相溶性[J].火炸藥, 1981(1):11-18.

[16] DUNKEL M.Calculation of intermolecular forces in organic compounds [J].PhysChem, 1928, 138: 42-54.

[17] FEDORS R F.A method for estimating both the solubility parameters and molar volumes of liquids [J].Polymer Engineering & Science, 1974, 14(2): 147-154.

[18] KREVELEN D W.Properties of polymers,3rd Edition [M].Elsevier Science Pub Co, 1990: 188-204.

[19] MANSON J A, SPERLING L H.Polymer Blends and Composites[M].Heyden, 1976.

[20] 李倩, 姚維尚, 譚惠民.疊氮粘合劑與硝酸酯溶度參數的分子動力學模擬[J].含能材料, 2007, 15(4):370-373.

[21] 李倩, 姚維尚, 譚惠民.分子動力學模擬疊氮熱塑性彈性體的楊氏模量及其與硝酸酯的溶度參數[J].火炸藥學報, 2007, 30(4):13-16.

[22] LI S, LIU Y, TUO X, et al.Mesoscale dynamic simulation on phase separation between plasticizer and binder in NEPE propellants[J].Polymer, 2008, 49(11):2775-2780.

[23] ABOU-RACHID H, LUSSIER L, RINGUETTE S, et al.On the correlation between miscibility and solubility properties of energetic plasticizers/polymer blends: modeling and simulation studies [J].Propellants, Explosives, Pyrotechnics, 2008, 33(4): 301-310.

[24] 許曉娟, 肖繼軍, 黃輝,等.ε-CL-20基PBX結構和性能的分子動力學模擬HEDM理論配方設計初探[J].中國科學:, 2007, 37(6):556-563.

[25] 焦東明, 楊月誠, 強洪夫,等.HTPB固體推進劑增塑劑選取分子模擬研究[J].化學研究與應用, 2009, 21(6):805-809.

[26] 王歡, 孫治丹, 張常山, 等.固體推進劑組分相容性的分子動力學模擬[J].火炸藥學報, 2016, 39(5): 69-73.

[27] 劉強.RDX 與 CL-20 及其共晶和復合體系的 MD 模擬研究[D].南京: 南京理工大學, 2014.

[28] 付一政, 康志鵬, 郭志婧, 等.共晶和共混對 CL-20/DNB 感度和熱解機理影響的 MD 模擬[J].含能材料, 2017, 25(2): 94-99.

[29] ROM N, ZYBIN S V, VAN DUIN A C T, et al.Density-Dependent Liquid Nitromethane Decomposition: Molecular Dynamics Simulations Based on ReaxFF[J].The Journal of Physical Chemistry A, 2011, 115(36): 10181-10202.

[30] LIU J, CHAMBREAU S D, VAGHJIANI G L.Thermal Decomposition of 1, 5-Dinitrobiuret (DNB): Direct Dynamics Trajectory Simulations and Statistical Modeling[J].The Journal of Physical Chemistry A, 2011, 115(28): 8064-8072.

[31] STRACHAN A, KOBER E M, VAN DUIN A C, et al.Thermal decomposition of RDX from reactive molecular dynamics[R].California Inst of Tech Pasadena Materials and Processes Simulation Center, 2005.

[32] STRACHAN A, VAN DUIN A C T, CHAKRABORTY D, et al.Shock waves in high-energy materials: The initial chemical events in nitramineRDX[J].Physical Review Letters, 2003, 91(9): 098301.

[33] RUSSO JR M F, BEDROV D, SINGHAI S, et al.Combustion of 1, 5-Dinitrobiuret (DNB) in the Presence of Nitric Acid Using ReaxFF Molecular Dynamics Simulations[J].The Journal of Physical Chemistry A, 2013, 117(38): 9216-9223.

[34] WEISMILLER M R, RUSSO JR M F, VAN DUIN A C T, et al.Using molecular dynamics simulations with a ReaxFF reactive force field to develop a kinetic mechanism for ammonia borane oxidation[J].Proceedings of the Combustion Institute, 2013, 34(2): 3489-3497.

[35] WEISMILLER M R, DUIN A C T, LEE J, et al.ReaxFF reactive force field development and applications for molecular dynamics simulations of ammonia borane dehydrogenation and combustion[J].The Journal of Physical Chemistry A, 2010, 114(17): 5485-5492.

[36] BROWN C M, JACQUES T L, HESS N J, et al.Dynamics of ammonia borane using neutron scattering[J].Physica B: Condensed Matter, 2006, 385: 266-268.

ApplicationofMolecularDynamicsMethodinResearchofSolidPropellants

MIAO Ruizhen, GAO Xiaomin, QIAO Xiaoping, CHEN Haiyang, XI Mingpo, MENG Yong

(Xi’an North Huian Chemical Industry Co., Ltd., Xi’an 710302, China)

The application progress of molecular dynamics simulation method inresearch of component compatibility, sensitivity and safety of composite solid propellants at home and abroad is summarized. In addition, the application of molecular dynamics in the study of detonation and combustion reactions in solid propellants is also summarized.It is possible to explain the microscopic reaction mechanism of explosive explosion at the molecular level.And the existing problems in applications and development prospects are analyzed. The references are provided for scientific researchers in this field to develop the research work of propellants and explosives using this method.

molecular dynamics simulations; solid propellant;properties; reaction mechanism

2017-08-15;

2017-08-30

苗瑞珍(1990—),女,碩士,主要從事含能材料結構與性能研究。

10.11809/scbgxb2017.12.010

本文引用格式:苗瑞珍,高曉敏,喬小平,等.分子動力學模擬在復合固體推進劑研究中的應用[J].兵器裝備工程學報,2017(12):40-44.

formatMIAO Ruizhen, GAO Xiaomin, QIAO Xiaoping, et al.Application of Molecular Dynamics Method in Research of Solid Propellants[J].Journal of Ordnance Equipment Engineering,2017(12):40-44.

TJ55

A

2096-2304(2017)12-0040-05

(責任編輯周江川)

猜你喜歡
體系研究
FMS與YBT相關性的實證研究
2020年國內翻譯研究述評
遼代千人邑研究述論
構建體系,舉一反三
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
探索自由貿易賬戶體系創新應用
中國外匯(2019年17期)2019-11-16 09:31:14
EMA伺服控制系統研究
新版C-NCAP側面碰撞假人損傷研究
如何建立長期有效的培訓體系
現代企業(2015年1期)2015-02-28 18:43:18
“曲線運動”知識體系和方法指導
主站蜘蛛池模板: 囯产av无码片毛片一级| 国产毛片网站| 国产精品自在自线免费观看| 亚洲福利一区二区三区| 午夜啪啪网| 国产成人av一区二区三区| 欧美A级V片在线观看| 中文字幕欧美日韩高清| 精品国产99久久| 四虎影视国产精品| 欧美在线中文字幕| 欧美三级视频网站| 亚洲精品视频免费看| 亚洲天堂日韩在线| 国产成人麻豆精品| 色呦呦手机在线精品| 国产在线无码一区二区三区| 91无码视频在线观看| 爽爽影院十八禁在线观看| 亚洲无码37.| 亚洲91在线精品| 国产精品爽爽va在线无码观看| 国产亚洲精品自在久久不卡| 免费看美女自慰的网站| 国产成人1024精品下载| 国产黄色爱视频| 91精品国产91久久久久久三级| 日韩性网站| 国产成人一级| 国产在线日本| 久久久久九九精品影院| 成人久久精品一区二区三区| 谁有在线观看日韩亚洲最新视频| 五月天久久综合国产一区二区| 欧美亚洲欧美| 综合社区亚洲熟妇p| 91极品美女高潮叫床在线观看| 精品国产成人a在线观看| 亚洲成人免费看| 亚洲第一黄片大全| 又大又硬又爽免费视频| 精品视频一区二区观看| 日本一区高清| 在线观看91香蕉国产免费| 又黄又爽视频好爽视频| 国产成人综合亚洲网址| 一本大道视频精品人妻| 久久午夜夜伦鲁鲁片不卡| 天天综合色网| 日韩视频福利| 免费视频在线2021入口| 狠狠色狠狠综合久久| 91青青草视频| 成人国产一区二区三区| 欧美区在线播放| 免费人成在线观看成人片| 国产熟睡乱子伦视频网站| 成人午夜天| 手机精品福利在线观看| 亚洲日韩欧美在线观看| 国产99免费视频| 在线亚洲精品自拍| 成人午夜视频网站| 99这里精品| 狠狠色丁香婷婷| 天天综合色天天综合网| 人人艹人人爽| 18禁高潮出水呻吟娇喘蜜芽| 国产微拍一区| 久久伊伊香蕉综合精品| 日韩黄色在线| 91麻豆国产在线| 欧美一级片在线| 伊人久久综在合线亚洲2019| аv天堂最新中文在线| 亚洲第一黄片大全| 亚洲天堂首页| 日本人妻一区二区三区不卡影院 | 99热这里只有精品免费| 亚洲一区二区三区香蕉| 99视频精品全国免费品| 性欧美在线|