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

40%DCP溶液的熱分解模型

2017-05-23 00:27:31董澤陳利平陳網樺馬瑩瑩
化工學報 2017年5期
關鍵詞:方法模型

董澤,陳利平,陳網樺,馬瑩瑩

?

40%DCP溶液的熱分解模型

董澤1,陳利平1,陳網樺1,馬瑩瑩2

(1南京理工大學化工學院,江蘇南京210094;2索爾維投資有限公司(中國),上海201108)

準確的熱分解動力學模型有助于人們采取各種安全措施預防和控制物料熱失控導致的燃燒爆炸事故。以40%過氧化二異丙苯(DCP)的2,2,4-三甲基戊二醇二異丁酯(DIB)溶液為研究對象,運用差示掃描量熱儀(DSC)和絕熱量熱設備(VSP2)進行了量熱實驗,并采用TSS軟件(Thermal Safety Software)對數據進行動力學分析,建立了兩種分解模型:“級+級”模型(模型1)和“級+自催化”模型(模型2),采用Friedman法和非線性擬合方法求算其動力學參數。在運用所建立的兩種模型擬合曲線時,發現兩種模型對同種量熱模式數據擬合的相關系數非常接近,說明單一量熱模式在求算動力學上存在局限性。聯合采用基于動態掃描模式的DSC數據及基于絕熱模式的VSP2數據共同求算動力學,發現相對于模型2,模型1可以更好地反映分解過程,其兩步反應的活化能分別為115.5 kJ·mol-1和135.7 kJ·mol-1,指前因子的對數分別為28.3和31.6,反應級數分別為0.40和0.84。研究結果表明采用基于不同量熱模式的數據求算動力學有助于確定正確的動力學模型,從而獲得準確的動力學參數,并克服單一量熱模式下動力學求算的局限性。

熱分解;量熱模式;40%DCP溶液;動力學;分解模型

引 言

建立準確的動力學模型有助于人們了解物質分解的各種行為,有助于采取各種安全措施預防和控制物質熱失控事故,因而對物質熱安全性的研究具有重要的現實意義。求解分解反應動力學參數通常分為無模型方法和模型擬合方法。無模型方法是指在不知道確切動力學模型的情況下獲得活化能,當然,這并不意味著無模型方法不能用來確定反應模型,只要在整個反應進程中,活化能的值不發生較大變化,即可認為分解反應滿足用單步反應動力學描述的條件,從而可以方便地采用Malek方法[1]確定反應模型和指前因子。國際動力學委員會(ICTAC)建議采用包括Friedman方法[2]、Ozawa-Flynn-Wall方法[3-4]在內的等轉化率方法對上述條件進行有效性驗證[5-7]。而當反應不能僅用一個單步反應模型和一組Arrhenius參數描述時,即出現多步反應動力學情形時,動力學模型的求解只能采用模型擬合方法。模型擬合方法又可以分為線性模型擬合方法和非線性模型擬合方法。線性模型擬合方法[8-9]通過取對數等數學手段,把反應速率方程轉化為線性方程;該方法更適合用于單步動力學求算,多步反應往往很難被線性化。非線性模型擬合方法通過使擬合數據與測試數據間的差異最小化來實現,適用于單步和多步反應動力學求算。模型擬合方法的關鍵在于動力學模型的確定,特別是多步反應模型的情況,每一步都需要確定一個合適的反應模型。不恰當的反應動力學模型往往導致動力學求解的完全錯誤。

然而模型的建立通常十分困難,有時不同的反應模型都能較好地描述反應,對于多步反應動力學模型的建立就更是如此。針對這種情況,本文以40%過氧化二異丙苯(DCP)的2,2,4-三甲基戊二醇二異丁酯(DIB)溶液(以下簡稱40%DCP溶液)為研究對象開展動力學篩選研究,其中,DCP是一種典型的有機過氧化物,常用于各種不同的樹脂中來改善材料的物理性質[10],國內外學者基于不同的量熱設備研究了DCP、DCP溶液的熱分解性能及熱力學和動力學參數[11-15]。DIB則是一種常用的工業溶劑和環保增塑劑,能較好地溶解過氧化二異丙苯,保證DCP是以溶液形式進行分解的,且其沸點較高,產生的蒸氣壓較小,而不會對產氣的DCP分解產生影響。利用基于動態掃描模式的差示掃描量熱儀(DSC)和基于絕熱模式的絕熱量熱設備(VSP2),進行了量熱實驗,并且采用基于兩種模式的測試結果對反應動力學模型進行求解,以期篩選出正確的動力學模型,并擬合獲得相應的動力學參數。

1 試劑及測試條件

1.1 試劑

DIB的純度為98.5%,密度為0.941 g·ml-1,沸點為280℃,CAS號為6846-50-0。DCP的純度為98%,密度為1.56 g·ml-1,熔點在40℃左右,CAS號為80-43-3。生產商都為阿拉丁試劑公司。

1.2 測試條件

所用DSC型號為200F3,由德國Netzsch公司生產,樣品池為高壓316不銹鋼坩堝,樣品質量(8.000±0.020) mg,溫度范圍30~200℃,動態溫升速率分別用1、2、4、8 K·min-1。

所用VSP2由美國Fauske公司開發生產,采用薄壁316不銹鋼球為測試球,實驗所用測試球的質量為39.75 g,容積為115 ml,所用樣品質量為44.07 g,溫度范圍30~300℃。升溫程序為:30~100℃,2 K·min-1線性升溫;100~300℃用加熱-等待-搜索模式,搜索的靈敏度為0.1 K·min-1 [16]。

2 結果與討論

2.1 DSC結果及分析

DSC動態模式下得到的40%DCP溶液測試曲線見圖1、圖2,其中,圖1是放熱速率與溫度的曲線(原始測試曲線),圖2是放熱速率與時間的曲線(已取基線,并截取了放熱段,由于圖1曲線有重合且不夠清晰,故附加了圖2)。從圖2中明顯觀察到40%DCP溶液的放熱曲線包含兩個相互疊加的放熱峰,第1個放熱峰尚未結束,第2段放熱已經開始。在圖1中,主峰(第2個分解峰)的起始分解溫度隨溫升速率的增大而升高,其范圍為131~140℃,放熱量的平均值為340 J·g-1,與文獻中DCP的分解溫度放熱量非常接近[17-21];次峰(第1個分解峰)的起始溫度也隨溫升速率增大而升高,其范圍為105~115℃,放熱量的平均值為7 J·g-1(具體見表1)。為了了解次峰的形成原因,在2 K·min-1升溫速率下,對DIB溶劑進行了動態DSC測試,結果見圖3,發現DIB在120~170℃范圍內有一段微弱的放熱過程,其起始溫度為122.5℃,比40%DCP溶液中次峰的起始分解溫度略高,比放熱量為5.4 J·g-1,與40%DCP溶液中次峰的放熱量比較接近,從而判斷次峰是由DIB放熱引起,相關結果見表2。

表1 40% DCP溶液動態DSC實驗數據

表2 98.5%DIB溶劑動態DSC實驗數據

按照ICTAC的建議,采用無模型的Friedman法[22]求取動態模式下40%DCP溶液的活化能值,該方法中活化能會隨著溫度與轉化率的變化而變化[23-26],其求解公式見式(1)

式中,為轉化率;a為指前因子,s-1;為理想氣體常數,J·(mol·K)-1;()為機理函數;a為活化能,kJ·mol-1;a為溫度,K;為時間,s。

由于反應最初和結束階段實驗的誤差較大,因此僅分析范圍為0.10~0.90的情況,經計算,活化能的值為115~145 kJ·mol-1,且隨著反應進度呈現很大波動,說明40%DCP溶液的分解反應機理不符合單步動力學模型,具體結果見圖4。

2.2 VSP2結果

40%DCP溶液的VSP2實驗曲線見圖5,經過一段時間的加熱-等待-搜索循環后,樣品在121.5℃測試到放熱信號,之后溫度和壓力持續上升;在155.2℃時,溫度和壓力陡然上升,溫升速率突然變大,經過1.6 min在241.5℃左右放熱結束;樣品的最大溫升速率達到10.3℃·min-1,最大壓升速率達到2.65 MPa·min-1。

2.3 TSS軟件求解動力學

2.3.1 動力學模型及求解原理 采用TSS軟件(Thermal Safety Software)對40%DCP溶液分解的動力學進行求解。該高等熱分析軟件應用廣泛、功能強大,主要通過對物質分解過程的模擬來驗證所建立動力學模型的正確性,并從擬合結果中得到模型相關的動力學參數。軟件中包含有多種類型的分解模型,如-order模型、自催化模型、Erofeev模型(適用于顆粒固體的表面分解動力學)、Jander模型(適用于分散控制的動力學)等[27]。

對于-order模型,反應速率的表達式見式(2)

其中,為級反應的反應級數。

對于自催化模型,反應速率的表達式見式(3),右側第1項表示自催化模型中引發反應的速率,第2項表示自催化反應的反應速率。令=1/2,即自催化因子[28-30],其具體表述見式(4)。式中,E為引發反應與自催化反應活化能的差值。綜合式(3)~式(5),自催化的反應速率可由式(6)表示。

(4)

E=1-2(5)

令=2,=2

其中,1、2分別表示引發反應和自催化反應的指前因子,1、2為各自的活化能,1、2為反應的反應級數。

熱平衡方程的公式為

在有了動力學模型后,將其與熱平衡方程式(7)結合。其中,轉化率和樣品溫度s是由實驗測得;在動態測試中(q為任意時刻的放熱速率),在絕熱測試中=(T-0)/Dad(T是放熱開始后任意時刻的溫度,0是起始分解溫度,Dad是放熱段的總溫升)。對于絕熱的VSP2來說,體系內不存在溫差從而不需要考慮熱傳導,只需擬合∞即可。而對于動態DSC,∞可以通過測試直接得到。因此,動力學參數0、、等只需在給出初始猜測值的基礎上進行非線性擬合即可獲得。這部分工作主要借助TSS完成。

2.3.2 40%DCP溶液的動力學研究 在對40%DCP溶液的DSC和VSP2數據單獨進行動力學建模時,發現有兩種模型可以較好地描述測試結果,分別是“級+級”模型(模型1)和“級+自催化”模型(模型2)。其中模型1中的兩個級反應分別是A(DIB)→B和E(DCP)→F;模型2中級反應指的是A(DIB)→B,自催化反應指的是E(DCP)→F。

用兩種模型分別對動態DSC和絕熱VSP2數據進行非線性擬合,得到4組動力學參數和4組擬合的動力學曲線,見表3和圖6~圖9。

表3 動力學參數模擬結果

結果表明:兩種模型對DSC數據的相關系數分別為0.9992和0.9994;對VSP2數據的擬合相關系數分別為0.9781和0.9726。顯然,單從擬合的相關系數上來看,這兩種模型都很適合描述40%DCP溶液在動態模式和絕熱模式下的分解。然而,理論上來說,物質分解的模型不會因為量熱模式不同而改變,因而兩種模型中必定有一個是不完全符合實際分解過程的。另外,兩種模型對同種量熱模式下測試數據的相關系數非常接近,這從一個側面說明了利用單一量熱模式的測試數據研究物質分解動力學是存在局限性的。

為此,采用同種模型同時對絕熱數據和動態數據共同求算動力學,以進一步篩選正確的動力學模型并求算準確的動力學參數。這種方法需要用一種模型同時對絕熱和動態數據進行擬合,通過不斷改變模型中的動力學參數,使擬合結果達到最優,即擬合曲線同時與兩種量熱模式測試曲線的相關系數最大,并獲得相應的動力學參數,擬合曲線見圖10~圖13。經計算,模型1和模型2擬合的相關系數分別為0.9798和0.9096,顯然,模型1的擬合結果明顯優于模型2,因此,可以認為40%DCP溶液分解過程應遵循模型1,具體動力學參數見表3。

從6組動力學參數對比中,發現模型1在對不同數據進行擬合后的結果是比較穩定的,動力學參數相差不大,而模型2的參數結果卻出現了較大的波動,這也從側面反映出模型1是40%DCP分解的準確動力學模型,與文獻[11-14]中純DCP或DCP溶液的分解模型相符。

3 結 論

(1)DSC的結果表明,40%DCP溶液的分解由兩個峰組成,第1個峰由溶劑的變化引起,第2個峰為DCP的分解。溶液中DCP分解的起始分解溫度為(135.0±5.0)℃,峰溫范圍為(162.0±8.0)℃,放熱量為(342.0±6.0) J·g-1。VSP2絕熱量熱實驗結果表明,40%DCP溶液的起始分解溫度為121.5℃,最大溫升速率對應的溫度為155.2℃。

(2)建立了兩種分解模型,模型1為“級+級”,模型2為“級+自催化”。在動態DSC數據動力學擬合的結果中,模型1、2的擬合相關系數分別為0.9992、0.9994。盡管擬合相關系數高,卻并不能說明所篩選的動力學模型是可靠的,表明基于單一量熱模式的測試數據在進行動力學計算中存在局限性。VSP2的動力學分析結果也同樣驗證了這一結論。

(3)同時對動態和絕熱測試結果進行動力學計算的結果表明,模型1(相關系數為0.9798)明顯優于模型2(相關系數為0.9096),說明同時采用兩種量熱模式的測試數據分析動力學是可行的,可以避免單一量熱模式在模型篩選上的局限性問題。

(4)40%DCP溶液熱分解的動力學可以表述為:

反應1

A(DIB)→B (級反應)

d/d=exp(28.3)exp(-115500/)(1-)0.40

反應2

E(DCP)→F (級分解反應)

d/d=exp(31.6)exp(-135700/)(1-)0.84

References

[1] MALEK J. The kinetic analysis of non-isothermal data[J]. Thermochimica Acta, 1992, 200(92): 257-269.

[2] SBIRRAZZUOLI N. Is the Friedman method applicable to transformations with temperature dependent reaction heat?[J]. Macromolecular Chemistry and Physics, 2007, 208(14): 1592-1597.

[3] OZAWA T. A new method of analyzing thermogravimetric data[J]. Bulletin of the Chemical Society of Japan, 1965, 38(11): 1881-1886.

[4] FLYNN J H, WALL L A. General treatment of the thermogravimetry of polymers[J]. Journal of Research of the National Bureau of Standards- Physics & Chemistry, 1966, 70A(6): 487-522.

[5] VYAZOVKIN S, BURNHAM A K, CRIADO J M,. ICTAC kinetics committee recommendations for performing kinetic computations on thermal analysis data[J]. Thermochimica Acta, 2011, 520(1/2): 1-19.

[6] VYAZOVKIN S. Modification of the integral isoconversional method to account for variation in the activation energy[J]. Journal of Computational Chemistry, 2001, 22(2): 178-183.

[7] BUDRUGEAC P. Differential non-linear isoconversional procedure for evaluating the activation energy of non-isothermal reactions[J]. Journal of Thermal Analysis and Calorimetry, 2002, 68(1): 131-139.

[8] PEREZMAQUEDA L A, CRIADO J M, GOTOR F J,. Advantages of combined kinetic analysis of experimental data obtained under any heating profile[J]. Journal of Physical Chemistry A, 2002, 106(12): 2862-2868.

[9] PEREZMAQUEDA L A, SANCHEZJIMENEZ P E, CRIADO J M. Combined kinetic analysis of solid-state reactions: ?a powerful tool for the simultaneous determination of kinetic parameters and the kinetic model without previous assumptions on the reaction mechanism[J]. Journal of Physical Chemistry A, 2006, 110(45): 12456-62.

[10] WU K W, HOU H Y, SHU C M. Thermal phenomena studies for dicumyl peroxide at various concentrations by DSC[J]. Journal of Thermal Analysis & Calorimetry, 2006, 83(1): 41-44.

[11] SOMMA I D, MAROTTA R, ANDREOZZI R,. Dicumylperoxidethermal decomposition in cumene: development of a kinetic model[J]. Industrial & Engineering Chemistry Research, 2012, 51(22): 601-609.

[12] LV J Y, CHEN L P, CHEN W H,. Kinetic analysis and self-accelerating decomposition temperature (SADT) of dicumylperoxide[J]. Thermochimica Acta, 2013, 571(20): 60-63.

[13] WU S H, SHYU M L, YETPOLE I,. Evaluation of runaway reaction for dicumyl peroxide in a batch reactor by DSC and VSP2[J]. Journal of Loss Prevention in the Process Industries, 2009, 22(6): 721-727.

[14] TALOUBA I B, BALLAND L, MOUHAB N,. Kinetic parameter estimation for decomposition of organic peroxides by means of DSC measurements[J]. Journal of Loss Prevention in the Process Industries, 2011, 24(4): 391-396.

[15] ZANG N. Thermal stability analysis of dicumylperoxide[J]. Lecture Notes in Electrical Engineering, 2012, 137: 279-286.

[16] ASKONAS C F, BURELBACH J P, LEUNG J C. The versatile VSP2: a tool for adiabatic thermal analysis and vent sizing applications[C]// 28th Annual North American Thermal Analysis Society (NATAS) Conference. Orlando, Florida, 2000: 4-6.

[17] WU S H, WANG Y W, WU T C,. Evaluation of thermal hazards for dicumyl peroxide by DSC and VSP2[J]. Journal of Thermal Analysis & Calorimetry, 2008, 93(1): 189-194.

[18] HOU H Y, LIAO T S, DUH Y S,. Thermal hazard studies for dicumyl peroxide by DSC and TAM[J]. Journal of Thermal Analysis & Calorimetry, 2006, 83(1): 167-171.

[19] WU K W, HOU H Y, SHU C M. Thermal phenomena studies for dicumyl peroxide at various concentrations by DSC[J]. Journal of Thermal Analysis & Calorimetry, 2006, 83(1): 41-44.

[20] LU K T, CHU Y C, CHEN T C,. Investigation of the decomposition reaction and dust explosion characteristics of crystalline dicumylperoxide[J]. Journal of Hazardous Materials, 2008, 161(1): 246-56.

[21] VALDES O J R, MORENO V C, WALDRAM S P,. Experimental sensitivity analysis of the runaway severity of dicumyl peroxide decomposition using adiabatic calorimetry[J]. Thermochimica Acta, 2015, 617(4): 510-513.

[22] FRIEDMAN H L. Kinetics of thermal degradation of char-forming plastics from thermogravimetry. Application to a phenolic plastic[J]. Journal of Polymer Science Part C Polymer Symposia, 2007, 6(1): 183-195.

[23] VYAZOVKIN S. Kinetic concepts of thermally stimulated reactions in solids[J]. International Reviews in Physical Chemistry, 2010, 19(1): 45-60.

[24] VYAZOVKIN S. On the phenomenon of variable activation energy for condensed phase reactions[J]. New Journal of Chemistry, 2000, 24(11): 913-917.

[25] VYAZOVKIN S. Chapter 13-Isoconversional Kinetics// Handbook of Thermal Analysis & Calorimetry[M]. Amsterdam: Elsevier, 2008, 5(8): 503-538.

[26] VYAZOVKIN S, SBIRRAZZUOLI N. Isoconversional kinetic analysis of thermally stimulated processes in polymers[J]. Macromolecular Rapid Communications, 2006, 27(18): 1515-1532.

[27] KOSSOY A, KOLUDAROVA E. Specific features of kinetics evaluation in calorimetric studies of runaway reactions[J]. Journal of Loss Prevention in the Process Industries, 1995, 8(4): 229-235.

[28] RODUIT B, HARTMANN M, FOLLY P,. Prediction of thermal stability of materials by modified kinetic and model selection approaches based on limited amount of experimental points[J]. Thermochimica Acta, 2014, 579(5): 31-39.

[29] RODUIT B, HARTMANN M, FOLLY P,. Thermal decomposition of AIBN, Part B: Simulation of SADT value based on DSC results and large scale tests according to conventional and new kinetic merging approach[J]. Thermochimica Acta, 2015, 355: 6-24.

[30] MOUKHINA E. Thermal decomposition of AIBN, Part C: SADT calculation of AIBN based on DSC experiments[J]. Thermochimica Acta, 2015, 621: 25-35.

Thermal decomposition model for solution of 40% dicumyl peroxide

DONG Ze1,CHEN Liping1,CHEN Wanghua1,MA Yingying2

(1School of Chemical Engineering, Nanjing University of Science & Technology, Nanjing 210094, Jiangsu, China;2Solvay (China) Co., Ltd., Shanghai 201108, China)

An accurate thermal decomposition model helps people take measures of prevention and controlling the burning and explosion accidents, which are caused by thermal run away of materials. This paper mainly studied a sample: 40%DCP solution (dicumyl peroxide dissolved in 2,2,4-trimethyl-1,3-pentanediol diisobutyrate). Experiments were performed with two devices, differential scanning calorimeter (DSC) and vent sizing package 2 (VSP2), analyzed the kinetics by TSS (Thermal Safety Software) and established two decomposition models: model 1 “-order followed with-order” and model 2 “-order followed with autocatalysis”, then, the kinetic parameters were estimated by Friedman method and non-linear simulation method. From simulated curves, both models described the decomposition curve for dynamic calorimetry mode or adiabatic calorimetry mode well, which illustrated the limitation of kinetics study for single calorimetric mode. Therefore, paper proposed a method that estimated the parameters based on both dynamic calorimetry mode and adiabatic calorimetry mode decomposition curves, found that only model 1 explained the decomposition well. Therefore, 40%DCP decomposition can be expressed as model 1, in which the activation energy for two steps are 115.5 kJ·mol-1and 135.7 kJ·mol-1, the natural logarithm of pro-exponential factor are 28.3 and 31.6, and reaction order are 0.40 and 0.84. This study proved that estimation of kinetic parameters with two different calorimetry modes can help determine right kinetic model, obtain accurate kinetic parameters, and overcome the limitation of kinetics study with single calorimetric mode.

thermal decomposition; calorimetry mode; 40%DCP solution; kinetics; decomposition model

10.11949/j.issn.0438-1157.20161516

TQ 013.1;TQ 013.2

A

0438—1157(2017)05—1773—07

陳利平。

董澤(1991—),男,博士研究生。

2016-10-27收到初稿,2016-12-22收到修改稿。

2016-10-27.

Prof. CHEN Liping, clp2005@hotmal.com

猜你喜歡
方法模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
學習方法
3D打印中的模型分割與打包
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 日韩精品高清自在线| 亚洲欧州色色免费AV| 亚洲成人黄色在线| 岛国精品一区免费视频在线观看 | 久久久久国产精品熟女影院| 无码区日韩专区免费系列| 青草午夜精品视频在线观看| 国产视频 第一页| 欧美成人影院亚洲综合图| 欧美自慰一级看片免费| 成人av手机在线观看| 国产丝袜一区二区三区视频免下载| 国产精品综合色区在线观看| 无码人妻免费| 国产91精品最新在线播放| 免费女人18毛片a级毛片视频| 人妻丝袜无码视频| 国产1区2区在线观看| 欧美日韩一区二区在线播放| 九九九国产| 国产黄网永久免费| 欧美性精品| 亚洲一区二区三区在线视频| 99er这里只有精品| 欧美一级片在线| 一级毛片在线播放| 91精品在线视频观看| 亚洲国产欧美国产综合久久| 欧美成人精品一级在线观看| 中文字幕在线免费看| 99久久成人国产精品免费| 久久久久亚洲精品成人网| 无码免费的亚洲视频| 欧美国产日韩在线观看| 91精品国产一区| 国产在线小视频| 亚洲专区一区二区在线观看| 久久先锋资源| 色呦呦手机在线精品| 影音先锋丝袜制服| 欧美一区二区三区不卡免费| 亚洲欧美成人综合| 欧美精品高清| 久久久久青草线综合超碰| 青青国产成人免费精品视频| 亚洲美女久久| 在线精品亚洲国产| 亚洲日本精品一区二区| 亚洲一道AV无码午夜福利| 91丨九色丨首页在线播放| 国模视频一区二区| 亚洲欧洲日韩国产综合在线二区| 国产成人精品亚洲77美色| 国产成人无码久久久久毛片| 国产浮力第一页永久地址| 国产1区2区在线观看| 日韩av手机在线| 在线视频97| 亚洲成人精品| 无码精油按摩潮喷在线播放| a亚洲天堂| 国产精品污视频| 三上悠亚一区二区| 亚洲嫩模喷白浆| 波多野结衣在线se| 一区二区三区国产| 青草视频网站在线观看| 国产凹凸一区在线观看视频| 国产精品久久久久久搜索| 亚洲天堂视频网站| 色婷婷啪啪| 国产本道久久一区二区三区| 亚洲综合久久成人AV| 免费观看成人久久网免费观看| 国产拍在线| 国产在线一区二区视频| 国产在线视频欧美亚综合| 日韩视频福利| 精品五夜婷香蕉国产线看观看| 国产又粗又猛又爽视频| 久久久久久高潮白浆| 亚洲AV一二三区无码AV蜜桃|