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

Pitzer活度系數模型研究與開發(fā)

2020-03-05 02:38:41韓莎莎鄭俊強孫曉巖項曙光
當代化工 2020年1期
關鍵詞:體系模型

韓莎莎 鄭俊強 孫曉巖 項曙光

摘 ?????要:基于Debye-Hückel(D-H)理論,Pitzer同時考慮了長程靜電相互作用能項、短程“硬核效應”的位能項以及三個離子間的相互作用能項據此修正了D-H理論,建立了適用于水電解質溶液體系(6摩爾離子強度以下)的活度系數模型。其中采用IAPWS-97物性方法計算純水的性質,進而對Pitzer活度系數模型進行研究、求解和開發(fā)。并采用實際工業(yè)生產中含電解質水溶液流程的應用實例進行模擬驗證,以及與Aspen plus化工模擬軟件在交互參數來源于相同文獻數據情況下的模擬結果相比,開發(fā)的模型在描述電解質水溶液體系的過量的吉布斯自由能、活度系數和焓值等熱力學性質方面的相對偏差低至0.1%左右,從而表明該模型的研究與開發(fā)可以應用于水電解質溶液體系的模擬、計算和優(yōu)化。

關 ?鍵 ?詞:電解質;活度系數;Pitzer模型;熱力學性質

中圖分類號:TQ013.1??????文獻標識碼:?A ??????文章編號: 1671-0460(2020)01-0221-07

Research and Development of Pitzer Activity Coefficient Model

HAN?Sha-shaZHENG?Jun-qiangSUN?Xiao-yanXIANG?Shu-guang

(Process Systems Engineering Institute, Qingdao University of Science and Technology, Shandong Qingdao 266042, China)

Abstract: Based on Debye-Hückel(D-H) theory, Pitzer took fully the long-range electrostatic interaction energy term, the potential energy term of the short-range "Hard-Core Effects" and the interaction energy between three-ions into account, and modified the D-H theory accordingly. An?activity coefficient model suitable for aqueous electrolyte solution system (below 6 mol ionic strength) was established. The properties of pure water were calculated by IAPWS-97 property method,and then the Pitzer activity coefficient model was researched, solved and developed.Using some application example of electrolyte-containing aqueous solution process in actual industrial production verified the results. Compared with the simulation results of chemical simulation software (Aspen plus)under the interactive parameters with the?same source of literature data, the relative deviations of the excess Gibbs free energy, the activity coefficient and enthalpies for aqueous electrolytes solution system that described by the developed model are all low to 0.1% or so. Therefore, the research and development of the model can be applied in?the simulation, calculation and optimization of aqueous electrolytes solution system.

Key words:Electrolytes;Activity coefficient; Pitzer model; Thermodynamic properties

自然界、生命體和工業(yè)過程中普遍存在著電解質溶液,是化工行業(yè)中的重要組成部分,也是眾多過程處理的對象,目前逐漸成為許多有機物和無機物反應的良好媒介,因此對電解質溶液的理論研究、電解質溶液的熱力學性質的研究及電解質過程模擬研究具有重要的工業(yè)實用價值和理論意義。

其中在電解質溶液理論及含電解質溶液體系的熱力學性質方面,Debye[1]、Meissner[2]、Bromley[3]、Chen[4]、陸小華[5]、左有祥[6,7]、Loehe[8]、李以圭[9]和杜艷萍[10]等都做出了很大的貢獻。目前Pitzer是用于計算水電解質溶液體系(尤其是離子強度為6摩爾以下的強電解質體系)的活度系數等熱力學性質較為準確的電解質活度系數模型,也是應用最為廣泛的電解質溶液理論。最初1973年,Pitzer修正了D-H理論[1],得到了經典的半經驗Pitzer模型[11],但適用的濃度較低。隨后為了擴大濃度適用范圍,用Margules方程修正了短程項,得到了Pitzer[12](1980年)模型。之后,Bromley[3](1973年)簡化的Pitzer模型、Pitzer[13](1975年)添加的靜電非對稱混合項、Fürst和Renon[14](1982年)研究的多種參數對模型用于1-1型電解質固液平衡的影響、李以圭[15,16](1986年)的Pitzer-Li方程、Simonson等[17](1986年)的Pitzer-Simonson方程、Kim等[18,19](1988年)回歸的高濃度體系參數、Clegg等[20,21](1992年)的Clegg-Pitzer模型、李以圭等[22,23](1994和1997年)的Li-Mather模型、Pitzer[24](1999年)以及Chen等[25](2008年)都對Pitzer模型做了相應的修正和完善。因此,參照Fortran語言編程如Zemaitis[26]中實現含電解質體系的模擬計算過程,也可通過Visual C++編程語言開發(fā)Pitzer模型,實現被已有的支持CAPE-OPEN標準的大型通用化工模擬軟件所調用,從而對工業(yè)中含電解質溶液過程進行設計、模擬、計算和優(yōu)化,更好地解決較復雜的工程問題。

本文主要是根據Pitzer修正的水電解質溶液體系活度系數計算模型[13](1975年模型)進行開發(fā)并通過對一些應用實例的模擬計算并驗證結果對該開發(fā)的Pitzer活度系數模型進行分析、討論和評價。

1 ?Pitzer活度系數模型

下文將對Pitzer活度系數模型作詳細闡述。

Pitzer活度系數模型以水電解質活度系數模型為基礎,區(qū)別于其他物性方法的是它與其他活度系數模型無重疊。它能夠準確地計算高達6 mol離子強度以上的電解質水溶液的行為,但不能用于存在除水之外的任何其他溶劑或混合溶劑的體系,任何非水的分子組分都應被認為是溶質,分子溶質存在與否也都可滿足,并全被看作亨利組分。超臨界氣體的溶解度是采用亨利定律來模擬計算的;汽相逸度系數采用Redlich-Kwong-Soave狀態(tài)方程,并且所有其他汽相性質也都假定是理想的,但遇到氣相中存在締合行為時,例如羧酸類(RCOOH)或氫氟酸(HF)體系,Redlich-Kwong-Soave則不適用。對于羧酸類,選擇一個非電解質活度系數模型與Hayden-O‘Connell或Nothnagel方法相結合;對于HF體系則選擇專用的ENRTL-HF或WILS-HF物性方法。

本文選取的Pitzer活度系數模型即水電解質活度系數模型計算的系統總過剩吉布斯能的基本方程見公式(1)。

(1)

式中:GE—剩余吉布斯自由能;

B、C、θψ—模型中交互參數;

??????R —氣體常數;

??????m —質量摩爾濃度,mol/kg;

??????Z —電荷數。

其中長程靜電相互作用項為:

(2)

即可推出活度系數模型:

(3)

其中:n—摩爾數;

*?—非對稱。

2 ?模型求解步驟

在進料物流的溫度T,壓力P和質量摩爾濃度mi或組成zi已知條件下,其各個種類(包括所有分子和離子組分)的活度系數的詳細求解步驟如下:

2.1 ?計算各類參數

Pitzer活度系數模型用于計算電解質體系所需要用到的參數主要包括:陽離子-陰離子間的參數β(0)β(1)β(2)β(3)C?φ;陽離子-陽離子間的參數θcc';陰離子-陰離子間的參數θaa';陽離子1-陽離子2-共同的陰離子三個離子間的參數ψcc'a;陰離子1-陰離子2-共同的陽離子三個離子間的參數ψcaa'以及分子-分子間、分子-離子間的參數β(0)、β(1)C?φ;還有靜電非對稱混合效應項Фij。以上所有參數都可參考Pitzer相關系列文獻[26-30],從而建立一個小型數據庫,供開發(fā)的模型模擬計算使用。

(1)計算陽-陰離子間的交互參數BijCij

因對于不同價型的電解質類型,Bij的計算公式也不相同,其中類似于1-n型(1-1,1-2和2-1等)電解質,BijCij都是針對帶不同電荷的離子即陰陽離子間的,其計算公式如下:

???(4)

針對n-m型的電解質并且n大于1、m也大于1(如2-2,2-3和3-4型等),Bij的計算公式如下:

???(5)

其中α1取值2.0;α2取值12.0;α3取值1.4。

相對應的求導得到Bij'(對于1-n型電解質):

???(6)

求導得到Bij'(對于2-2型電解質):

??(7)

????????????(8)

(2)計算帶同種電荷的離子間交互參數θij

???(9)

λ為短程離子間效應描述為離子強度函數。

(3)計算三元交互作用參數ψijk

(10)

Λ是忽略任何離子強度依賴的三離子交互項。

其中(2)和(3)計算的參數是多電解質組分溶液中可能存在的雙離子和三離子間的相互作用。對于單一的電解質溶液,它們的值都為零。

(4)計算靜電非對稱混合效應項Фij

(11)

求導得:

(12)

對于參數X的計算,根據體系類型不同,其計算方法也不同,當不存在非對稱的混合校正項時:

(13)

當靜電非對稱混合修正項用多項式表達時為:

(14)

當靜電非對稱混合修正項用積分形式表達時為:

(15)

2.2 ?計算各組分活度系數

長程作用項:

????(16)

其中Debye-Hückel參數Aφ、離子強度Ix為:

?????(17)

水的密度為:

?(18)

以上式子中:k —玻爾茲曼常數;

NA—阿伏伽德羅常數;

Qe—電子電荷;

??—介電常數;

M?—相對分子質量,kg/kmol。

fI)進行求導得:

?(19)

以上各類參數將其代入(3)式就可計算出體系中各個種類的活度系數。離子種類的活度系數模型:

?(20)

對于水,類似地計算活度系數的對數如下:

(21)

(22)

由于計算出的是以質量摩爾濃度mi為單位標度的活度系數值,可由下列公式將其轉化為以摩爾分數xi為單位標度:

(23)

離子類:

(24)

水分子:

(25)

2.3 ?計算Gm*EHm*E

在上一節(jié)中求出lnγi之后,即可求:

???????(26)

由于計算Hm*E需要對活度系數表達式求導,為了避免復雜的求導過程,可以采用以下方式計算:

?????????(27)

???(28)

很明顯看出,這里采用的是對Gm*E進行差分法求導計算:

?(29)

其中T2-T1可以取一個很小的值0.001。

Hm*EGm*E計算完成之后,再按照常規(guī)體系的總焓和總吉布斯自由能計算模型得出總焓值和總吉布斯自由能值,既而得到總熵值。

3 ?模型開發(fā)

通過Visual studio 2005平臺采用C++語言進行編程實現,具體開發(fā)過程如下。

(1)定義變量

//定義組分列表

vector<int>compIdList;

vector<BOOL>henryList;

//定義組分組成

vector<double>xm;

//定義參數數組集

vector<vector<double>>B0ij,B1ij,B2ij, B3ij,Cij,θij,ψijk;

(2)函數類型

//創(chuàng)建計算活動系數所需的相互作用參數數組。

voidInitialParameters();

//計算長程靜電相互作用項

CalcLongRange(double T,vector<double>&x,vector<double>& LR);

//計算短程位能項

CalcShortRange(double T,vector&x,vector& SR);

//計算三離子交互項

CalcThreeIonInteraction (double T,vector&x,vector& TII);

//計算各個種類的活度系數

GetActivityCoefficient(double T,vector&x,vector& ac);

GetActivityCoefficient(CVariant&compIdList, CVariant& composition, double P, double T, VARIANT* value);

//計算Hm*E

GetExtraEnthalpyDepature(CVariant&compIdList, CVariant& composition, double P, double T, double * extraEnthalpy);

(3)對亨利組分的處理

//初始化一個數組,對應于所有的組分,如果組分為亨利組分,則對應的位置為,否則為0

henryList=vector(count,0);

for(inti=0;i

if(henryComps.GetLongAt(i)==compIdList[j])

{henryList[j]=1; break;}//以上是亨利組分的識別,識別出的亨利組分和可凝組分的交互作用參數都設為0。

(4)計算純水的性質

CComPtrutilitiesManager;

//公用工程的接口,計算純水摩爾體積和焓值

utilitiesManager->GetMoleWaterVolumeByTP(T,P,&moleWaterVolume);

HRESULT hr=utilitiesManager->GetMoleEnthalpyByTP(T,P,value);

4 ?模型驗證與對比

以實際工業(yè)生產中的含電解質流程實例對開發(fā)的Pitzer模型進行模擬驗證,體系中涉及的汽相的性質采用RK-Soave狀態(tài)方程計算,液相的熱力學性質則采用本文開發(fā)的Pitzer模型進行模擬。最后采用自身開發(fā)的數據庫的基礎物性以及Pitzer系列文獻中[26-30]的交互作用參數值(包括1-1型、2-1型等電解質間參數),將此開發(fā)模型的模擬計算結果(OPEN)與Aspen plus 7.2化工模擬軟件的模擬結果(AP)進行對比分析,進一步驗證開發(fā)模型的可行性和準確性。

4.1 ?強電解質水溶液體系

強電解質水溶液體系包括單一強電解質和多組分電解質水溶液體系,大多是水鹽體系。下面選取表1中的三種常見模擬體系:氫氧化鉀水溶液體系、氯化鎂水溶液體系以及苛性鈉溶液體系實例進行模擬計算,由于這些體系中模擬組分涉及的均是強電解質,所以在此不予考慮水自身的電離影響(因為水本身是一種弱電解質,電離極其微弱)。其中體系中純水組分的性質如水的液相摩爾體積、液相焓值等都是采用自己開發(fā)的IAPWS-97物性方法(水的蒸汽表)進行計算的,其進料溫度均設為35.5 ℃,壓力為101 325 Pa,總進料量100 kmol·h-1,各個體系所包含的模擬組分及組成如表1所示,根據第二節(jié)模型的計算步驟和第三節(jié)經編程開發(fā)實現的模型進行模擬計算,計算出的各個組分的液相活度系數值結果對比見表2所示。

各個體系的焓、熵和吉布斯自由能值結果對比見表3。

4.2 ?含分子溶質的弱電解質水溶液體系

含分子溶質的弱電解質水溶液體系包括:帶有酸性氣體的鹽水溶液,帶有酸性氣體的乙醇胺、二乙醇胺、二甘醇胺和甲基二乙醇胺等。本小節(jié)選取易揮發(fā)性弱電解質體系NH3-H2O-CO2體系進行模擬計算,體系中純水的性質也是采用自己開發(fā)的IAPWS-97物性方法(水的蒸汽表)進行計算的并且本體系中要考慮水的電離(因有弱酸根離子存在),其體系進料溫度為60 ℃,壓力為1.0 atm,體系中的模擬組分及進料量如表4所示,通過物料平衡、化學反應平衡以及相平衡聯立計算(期間多次調活度系數模型)可以得出可能存在的所有種類(包括分子和離子)具體的摩爾分數,再取一組與Aspen Plus相同的組成根據第二節(jié)模型的計算步驟和第三節(jié)經編程開發(fā)實現的模型進行模擬計算,其各個組分的活度系數模擬計算結果對比見表5,體系的焓、熵和吉布斯自由能值模擬結果對比見表6。

根據表2-3和表5-6的模擬結果表明:在與化工模擬軟件Aspen Plus相同的進料條件下,各種二元或三元的交互作用參數都參考文獻中的數據,針對4.1中的三個強電解質水溶液體系的熱力學性質的偏差都很小,基本上是一致的;4.2中的易揮發(fā)性弱電解質體系的計算因有弱電解質電離及弱酸根離子的水解等復雜化學反應的存在,其各種性質的計算結果相對偏差較大一點,但也都控制在0.1%左右。

綜上,本文研究并開發(fā)的Pitzer活度系數模型的模擬結果在離子強度為6摩爾以下的電解質體系方面和Aspen Plus軟件的計算結果基本上是一致的,從而表明,采用我們自身的數據庫以及純水的計算方法來驗證開發(fā)的模型是準確的,因此模型的開發(fā)也基本上是正確的。

5 ?結果與討論

基于Debye-Hückel(D-H)理論,根據Pitzer同時考慮長程靜電相互作用能項、短程“硬核效應”的位能以及三個離子間的相互作用能項修正的D-H理論,建立并開發(fā)適用于強電解質水溶液體系和易揮發(fā)弱電解質溶液體系的模型,可以得出以下結論:

(1)選取的計算純水性質的物性方法(自身開發(fā)的IAPWS-97方法)正確;

(2)該模型在描述電解質溶液的非理想性時綜合考慮了長程項、短程項以及離子作用項等,是目前較為準確地計算低濃度下強電解質水溶液的熱力學性質的模型。

(3)通過與Aspen plus軟件的模擬結果的對比分析,采用本文開發(fā)的Pitzer活度系數模型在計算電解質水溶液體系的熱力學性質的相對偏差總體低至0.1%左右。

綜上所述,該模型的建立、求解及開發(fā)的計算結果正確,Pitzer模型是目前應用最多的,但由于計算過程中,參數較多,計算較為復雜,適用的溫度、壓力及濃度范圍也有限,因此為后人在電解質方面的半經驗模型的研究打下基礎。

參考文獻:

[1]Debye P,Hückel E. The Theory of electrolytes.1. Lowering of freezing point and related phenomena[J]. PhysikalischeZeitschrift, 1923, 24: 185-206.

[2]Meissner H P, Tester J W.Activity Coefficients of Strong Electrolytes in Aqueous Solutions[J]. Aiche Journal, 1972, 18(3):661-662.

[3]Bromley L A. Thermodynamic properties of strong electrolytes in aqueous solutions[J]. Aiche Journal, 1973, 19(2):313-320.

[4]C.-C. CHEN,H.I. BRITT, J.F. BOSTON, L.B. EVANS. Local Compositions Model for Excess Gibbs Energy of Electrolyte Systems: Part I: Single Solvent, Single Completely Dissociated Electrolyte Systems[J]. AIChE J., 1982, 28 (4):588-596.

[5]陸小華,王延儒,時均. 電解質溶液的熱力學進展[J]. 化工進展,1988, 7(3):13.

[6]Zuo Y.X., Guo T. M. Model of thermodynamic properties of electrolyte solution [J]. Journal of Petroleum University?(Natural Science Edition), 1989(04):107-115.

[7]左有祥, 郭天民. 電解質溶液的分子熱力模型[J]. 化工學報, 1990, 41(1):1-9.

[8]Loehe J R, Donohue M D. Recent advances in modelingther- modynamic properties of aqueous strong electrolytesystems [J]. AIChEJ.,1997, 43(1):180-195.

[9]Li Y. G. Electrolyte solution theory[M]. Beijing: Tsinghua University Press, 2005.

[10]杜燕萍, 劉斌, 侯海云,張東東,白海浪,杜麗穎. 電解質溶液熱力學模型的研究進展[J]. 當代化工, 2016, 45(11): 2632-2637.

[11]K. S. PITZER, Thermodynamics of electrolytes. I. Theoretical basis and general equations[J]. J. Phys. Chem., 1973, 77(2): 268-277.

[12]Pitzer K S. Electrolytes from diluted solutions to fuse salts[J]. J Am Chen Soc, 1980,102(9): 2902-2906.

[13]Pitzer K S. Thermodynamics of electrolytes. V. effects of higher-order electrostatic terms[J]. Journal of Solution Chemistry, 1975, 4(3): 249-265.

[14]Fürst, W., H. Renon. Effects of the Various Parameters in the Application of Pitzer's Model to Solid-Liquid Equilibrium Preliminary Study for Strong 1-1 Electrolytes[J]. Ind. Eng. Chem. Process Des.Dev., 1982, 21(3): 396-400.

[15]li y g, pitzer k s. thermodynamics of aqueous sodium chloride solutions at high temperatures and pressures(i)thermodynamic properties over 373-573 k and 0.1-100 MPA[J]. Chinese Journal of Chemical Engineering, 1986, 1(2): 138-149.

[16]李以圭, Pitzer K S. 高溫高壓下氯化鈉水溶液的熱力學(Ⅱ)——在823 K及100 MPa條件下NaCI-H2O體系的熱力學性質[J]. 化工學報, 1987, 37(2): 249-255.

[17]Pitzer K S , Simonson J M . Thermodynamics of multicomponent, miscible, ionic systems: theory and equations[J]. The Journal of Physical Chemistry B, 1986, 90(13): 3005-3009.

[18]Kim H T , Frederick W J . Evaluation of Pitzer ion interaction parameters of aqueous electrolytes at 25.degree.C. 1. Single salt parameters[J]. Journal of Chemical & Engineering Data, 1988, 33 (2): 177-184.

[19]Kim H T, Frederick W J. Evaluation of Pitzer ion interaction parameters of aqueous mixed electrolyte solutions at 25°C. Part 2: Ternary mixing parameters[J]. Journal of Chemical & Engineering Data, 1988, 33 (3): 278-283.

[20]Clegg S L, Pitzer K S. Thermodynamics of Multicomponent, Miscible, Ionic Solutions: Generalized Equations for Symmetrical Electrolytes[J]. The Journal of Physical Chemistry, 1992, 96(8): 3513-3520.

[21]Clegg S L, Pitzer K S, Brimblecombe P . Thermodynamics of multicomponent, miscible, ionic solutions. Mixtures including unsymmetrical electrolytes[J]. Journal of Physical Chemistry, 1992, 96 (23): 9470-9479.

[22]Li Y G , Mather A E . The Correlation and Prediction of the Solubility of Carbon Dioxide in a Mixed Alkanolamine Solution[J]. Industrial & Engineering Chemistry Research, 1994, 33(8): 2006-2015.

[23]Li Y, Mather A E. Correlation and Prediction of the Solubility of CO2?and H2 S in Aqueous Solutions of Triethanolamine[J]. Industrial & Engineering Chemistry Research, 1995, 34(7): 2545-2550.

[24] Pitzer K S , Wang P , Rard J A , et al. Thermodynamics of Electrolytes. 13. Ionic Strength Dependence of Higher-Order Terms; Equations for CaCl2?and MgCl2[J]. Journal of Solution Chemistry, 1999, 28(4): 265-282.

[25]Chen C C, Song Y, Bollas G M. Formulation and Behavior of a Symmetric Electrolyte Nrtl Model for Gibbs Energy of Electrolyte Systems[C]. Aiche Meeting. 2008.

[26]Zemanitis J F, Clark D M ,Rafal M, et al. Handbook of Aqueous Electrolyte Thermodynamics: Theory & Application[M]. The American Institute of Chemical Engineers, 1986.

[27]Pitzer K S, Mayorga G. Thermodynamics of electrolytes. II. Activity and osmotic coefficients for strong electrolytes with one or both ions univalent[J]. The Journal of Physical Chemistry, 1973, 77 (19): 2300-2308.

[28]Pitzer K S, Mayorga G. Thermodynamics of electrolytes. III. Activity and osmotic coefficients for 2-2 electrolytes[J]. Journal of Solution Chemistry, 1974, 3(7): 539-546.

[29]Pitzer K S, Kim J J. Thermodynamics of electrolytes. IV. Activity and osmotic coefficients for mixed electrolytes[J]. Journal of the American Chemical Society, 1974, 96(18): 5701-5707.

[30]Pitzer, K.S., J.R. Peterson, and L. F. Silvester, Thermodynamics of Electrolytes. IX. Rare Earth Chlorides, Nitrates, and Perchlorates[J]. Solution Chem., 1978, 7(1): 45-56.

猜你喜歡
體系模型
一半模型
構建體系,舉一反三
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
探索自由貿易賬戶體系創(chuàng)新應用
中國外匯(2019年17期)2019-11-16 09:31:14
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
如何建立長期有效的培訓體系
“曲線運動”知識體系和方法指導
“三位一體”德育教育體系評說
中國火炬(2010年7期)2010-07-25 10:26:09
主站蜘蛛池模板: 欧美成人手机在线视频| 亚洲综合中文字幕国产精品欧美| 无码精品一区二区久久久| 日韩欧美国产三级| 久久国产精品嫖妓| 激情网址在线观看| 欧美精品另类| 91亚洲视频下载| 国产在线日本| 中日韩欧亚无码视频| 久久77777| 天天婬欲婬香婬色婬视频播放| 成年女人a毛片免费视频| 亚洲午夜福利在线| lhav亚洲精品| 国产99视频精品免费观看9e| 五月天香蕉视频国产亚| 国产AV毛片| 蜜桃视频一区二区三区| 婷婷在线网站| 欧美国产综合视频| 久草网视频在线| 亚洲狼网站狼狼鲁亚洲下载| 国产三级视频网站| 一本久道热中字伊人| 欧美一级特黄aaaaaa在线看片| 亚洲不卡网| 999国产精品永久免费视频精品久久| 亚洲永久视频| 国产在线精品美女观看| 国产精品一区二区国产主播| 亚洲天堂日韩av电影| 高清欧美性猛交XXXX黑人猛交 | 免费Aⅴ片在线观看蜜芽Tⅴ| 欧美亚洲国产精品第一页| 国产精品久久久久久久久| 亚洲天堂区| 国产第一页亚洲| 中文一区二区视频| 米奇精品一区二区三区| A级毛片无码久久精品免费| 午夜视频在线观看免费网站| 色综合中文| 试看120秒男女啪啪免费| 精品视频一区二区三区在线播| 免费观看三级毛片| 亚洲乱码视频| 国产美女自慰在线观看| 91小视频版在线观看www| AV无码无在线观看免费| 国产人妖视频一区在线观看| 国产成人精品无码一区二| 丝袜高跟美脚国产1区| 2024av在线无码中文最新| 东京热高清无码精品| 伊人久久青草青青综合| 国产福利拍拍拍| 天堂av高清一区二区三区| 永久天堂网Av| 亚洲综合狠狠| 欧美性精品| 91激情视频| a色毛片免费视频| 四虎影视无码永久免费观看| 91香蕉视频下载网站| 91啦中文字幕| 夜夜操天天摸| 国产三级国产精品国产普男人| 国产精品浪潮Av| a毛片在线免费观看| 亚洲精品天堂自在久久77| 国产真实乱人视频| 片在线无码观看| 99视频在线精品免费观看6| 色天堂无毒不卡| 在线精品亚洲国产| 有专无码视频| 国产久操视频| 国产h视频在线观看视频| 在线观看欧美国产| 国产女人爽到高潮的免费视频| 国产不卡在线看|