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

微通道內連續合成十二烷基苯磺酸的響應面分析及混合過程模擬

2021-11-30 07:40:16孟維軍徐一鳴李平趙曦嚴佩蓉徐建鴻
化工進展 2021年11期

孟維軍,徐一鳴,李平,趙曦,嚴佩蓉,徐建鴻

(1 寧夏大學省部共建煤炭高效利用與綠色化工國家重點實驗室,銀川 寧夏 750000;2 清華大學化學工程系,化學工程聯合國家重點實驗室,北京 100084)

20 世紀40 年代后,十二烷基苯磺酸鹽以其良好的生物降解性、價格低廉等優勢在全世界范圍內被廣泛應用為洗滌劑的原料。烷基苯磺酸鹽和石油磺酸鹽在我國被廣泛使用[1]。目前應用較為廣泛的十二烷基苯(DDB)磺化反應器有釜式反應器、降膜反應器和噴射環流反應器,但這些反應器都存在著一定的問題,導致十二烷基苯磺酸產品成本過高。因此,高效、綠色、安全的磺化反應器是降低表面活性劑成本、促進陰離子表面活性劑供給側結構性改革的重要環節。微化工設備如微混合器、微反應器、微換熱器等通過精密加工技術制造的設備,以其1mm 以下的通道特征尺寸在近十幾年里受到普遍關注[2-4]。與傳統尺寸設備相比,微反應器和微混合器特征尺寸小、比表面積大、傳質傳熱能力比傳統化工設備大1~2個數量級,能夠在化工過程中提供更精確的溫度控制、更強的混合效果、更窄的停留時間分布以及更安全的工藝狀態[5]。因此,許多瞬間強放熱反應、有毒中間體反應、易爆炸反應都被應用于微化工設備中[6-10]。

十二烷基苯磺化反應是瞬間強放熱反應,并附帶平行副反應,而副產物對產品色澤影響嚴重,對于這一情況,人們希望磺化反應器能夠有足夠的傳熱性能、傳質性能和較短的停留時間,來控制局部的溫度、濃度對反應的影響。主凱等[11]以1,2-二氯乙烷作為溶劑,研究了DDB 在微反應器中用液相SO3進行磺化的過程特性,發現DDB磺化反應在微通道反應器內主要由傳質效果控制。Geng等[12]研究了微反應器連續化工藝和循環化工藝對液相三氧化硫DDB 磺化工藝的區別,通過優化流量、DDB 摩爾比、反應溫度等因素最終使得在保證產品活性物質含量高達97%的同時,停留時間為2.5s,老化時間為1h,遠遠快于現有工業過程。可見,微反應器可以有效地應用在DDB 磺化工藝中。但目前并沒有有關工藝參數對微通道反應器內磺化反應過程的交互影響的研究,各因素的交互影響對于反應過程的控制和優化有著重要意義。

由于特征尺寸的減小,流體在微通道反應器內受力情況由傳統尺寸化工設備中的慣性力主導向黏性力和表面張力主導轉變,因此流體的傳遞和混合特性也表現出一些新現象。關于如何高效地強化微尺度下液-液混合與傳質過程,許多研究者進行了大量的研究[13-17]。通過結構或外場作用對微通道中局部進行二次流混合強化是一種有效的強化手段,其中以T形撞擊流微混合器在工業化應用中最有前景[18]。Fonte等[19]采用平面激光誘導熒光(planar laser induced fluorescence,PLIF)技術對常規尺寸通道內對撞流的流體流動狀態和混合質量進行了研究,發現在當Re>104 時,兩股等動量圓射流流體撞擊面發生偏轉振蕩,不再維持相互分離的狀態。Madana和Ashraf[20]使用計算流體力學方法(computational fluid dynamics,CFD)對矩形截面T形微通道內的單相和多相流進行了分析,發現流態的轉變存在臨界Re,最終提出了一個入口幾何形狀的優化方案,以保證在低Re下實現流態的轉變。而對于微通道反應器來說,溶劑的加入能夠改變流體的密度、黏度和不同組分的溶解度,從而影響流體的Re和混合質量,但目前對于有溶劑條件下T形微通道反應器的對撞流流態和混合質量的研究還鮮有報道。

本文通過Design-Expert 實驗設計軟件設計了響應面實驗,以研究微反應器內溫度、摩爾比、體積流量對產物活性物質含量的交互影響。并通過ANSY-Fluent流體力學計算軟件對有溶劑條件下的T形微通道混合器對撞流流體進行了模擬,以表征微混合器中流體流動的狀態,通過對流體流動和混合機理的研究來更好地優化和控制微反應器內十二烷基苯的磺化工藝。

1 實驗部分

1.1 實驗原料及主要試劑

本實驗用到的十二烷基苯(DDB)[異構體混合物,阿法埃莎(中國)化學有限公司]與1,2-二氯乙烷(分析純,General-Reagent)配制成十二烷基苯溶液。SO3由30%發煙硫酸和五氧化二磷[98%,阿法埃莎(中國)化學有限公司]反應制得,溶于1,2-二氯乙烷制成質量分數為7.5%的液相SO3溶液。

實驗過程中配置SO3溶液后,通過使用0.1mol/L的NaOH 標準溶液滴定得到SO3溶液的物質的量濃度,再根據所需的SO3/DDB摩爾比配置DDB溶液。

1.2 反應機理

主反應式如圖1所示[21]。

圖1 反應機理

十二烷基苯(DDB)與SO3反應生成焦磺酸中間體,焦磺酸與DDB反應生成目標產物十二烷基苯磺酸(DBSA)以及主要副產物二芳基砜,副反應為可逆反應,主反應與副反應為平行反應。老化階段焦磺酸中間體與酸酐和DDB生成目標產物DBSA。

1.3 實驗流程

SO3溶液與DDB 溶液以1∶1 的流量分別通過兩個平流泵通入T形混合器內,然后進入反應器管道進行反應。混合器和反應器管道放置在水浴鍋中恒溫,產品接出后使用旋蒸儀除去產品中的二氯乙烷溶劑,旋蒸儀水浴溫度40℃。T 形混合器內徑1mm,微通道反應器內徑0.6mm,長度3000mm,材質均為304不銹鋼。實驗流程如圖2所示。

圖2 實驗流程

1.4 產品活性物質含量分析方法

依據GB/T 8447—2008、GB/T 5173—1995使用溴化十六烷基三甲基銨(CTAB)進行直接兩相滴定法測定產物中活性物質含量。依照式(1)計算。

式中,X為十二烷基苯磺酸產品活性物質質量分數,%;V為消耗的滴定液體積,mL;c為CTAB滴定液濃度,mol/L;M為十二烷基苯磺酸分子量,g/mol;m為取樣質量,g。

1.5 響應面實驗設計

使用Design-Expert(V10.0.4)中的Box-Behnken模塊設計實驗方案,采用三因素三水平中間點重復三次的設計方法,選取T 形混合器單側入口流量(A)、溫度(B)、SO3/DDB 摩爾比(C)3個因素為自變量,產品活性物質含量(X)為因變量,研究三個因素對產品活性物質含量的交互影響。T形混合器兩側入口流量比為1∶1,則單側入口流量分別為6mL/min、8mL/min、10mL/min 時對應的反應器內總流量為12mL/min、16mL/min、20mL/min,總流量對應的停留時間分別為4.23s、3.18s、2.54s。因素變量水平參考主凱等[11]和Geng等[12]的研究結果進行選取,見表1。

表1 實驗的因素水平表

中間點重復3次,共進行15次實驗,研究三個因素對產品活性物質含量的交互影響。實驗設計值及結果見表2。

表2 實驗設計與結果

2 結果與討論

2.1 響應面法優化產品活性物質含量

使用Design-Expert 軟件對結果進行三次回歸擬合,得到的產物活性物質質量分數Y(%)與流量(A,mL/min)、溫度(B,℃)、SO3/DDB 摩爾比(C)的多元回歸方程如式(2)所示。方差分析及可信度分析見表3、表4。

表4 響應面模型可信度分析

表3 響應面模型方差分析

由表3可知,回歸模型的F值為1696.67,且P值小于0.05,說明模型具有顯著意義。由表4 可知,模型的相關系數為0.9999,表明相關性良好;校正系數為0.9993,表明三次模型有效,且精密度高,可以用來分析結果。變異系數值為0.16%,說明實驗精確度高,擬合程度高。

2.2 因素交互作用分析

響應曲面圖能直觀反映各因素對因變量的影響以及兩兩因素之間的交互作用。圖3(a)為溫度和流量的響應曲面,圖3(b)為SO3/DDB 摩爾比和流量的響應曲面,圖3(c)為SO3/DDB 摩爾比和溫度的響應曲面。由圖3(a)可以看出,隨著摩爾比的增大,曲面斜率增大,產品活性物質含量受兩個因素的影響且都增加。其中,溫度和流量對產品活性物質含量的影響表現為:低流量時溫度對其影響不大,且隨著摩爾比的增大這一現象沒有改變;中流量時,溫度與因變量成正相關,隨著摩爾比的增大,溫度對其影響增大,但表現為低點降低,高點不變;高流量時,溫度對因變量的影響隨著摩爾比的增大而變強,但主要表現在高溫區。針對這些現象,分析是因為在T 形混合器中,流量的增大使得流體Re增大,對撞效果增強,同時在反應器長度不變的情況下,停留時間變短。因此在低流量時,由于混合效果不強,反應主要受傳質控制,受到溫度和摩爾比影響較小。中流量時,流量的增加使得停留時間變短,因此呈現為低溫區低點降低,高溫區高點不變。高流量時,混合效果增強,反應由傳質控制向反應控制轉變,且主反應相對SO3為二級反應,因此隨著摩爾比的增加,產物活性物質含量上升。由圖3(b)可以看出,隨著溫度的增加,曲面高度整體上升(曲線顏色變紅),即因變量數值整體增大,且整體趨勢由先減小后增大轉變為由低到高逐漸增大。結合圖3(a)現象,發現隨著溫度的增加,流量對因變量的影響發生了變化。分析認為當溫度升高到60~70℃時,反應物中有SO3氣體溢出,這一現象增強了反應器內的傳質和混合效果,因此在溫度為60~70℃時,之前由于停留時間減少而出現的中流量區的低點消失。從圖3(c)也可分析出以上結論。

圖3 響應曲面圖

3 混合過程模擬

通過響應面實驗發現,影響產品活性物質含量的主要因素為入口流量大小,認為該因素會影響混合器內混合效果,因此為表征微混合器內的流體流動和混合現象,使用CFD 的方法對微混合器內流體流動進行模擬。模擬了實驗中用到的十二烷基苯和SO3在二氯乙烷為溶劑條件下的混合狀態,與實驗所得結論進行對比。為使得混合器內混合現象的變化更加敏感,選用水和甲苯這兩種互不相溶液體作為主要混合流體,選用對這兩種組分都能溶解的乙醇液體作為溶劑,以表征不同條件下混合器內流體流動及混合狀態。

3.1 基本物理模型

文中針對實驗中用到的T形圓柱微通道混合器進行計算機模擬,采用SCDM軟件構造原尺寸三維幾何模型,如圖4 所示。T 形微混合器截面內直徑為1mm,入口與出口間管道長度24.5mm。兩股流體分別由兩個入口進入微混合器,入口邊界條件設為速度,出口設為壓力出口邊界條件。

圖4 T形混合器模型(單位:mm)

3.2 邊界條件控制方程

對于不同相態之間的模擬有許多模型,如流體體積函數、水平集、離散格子玻爾茲曼、標記粒子、界面追蹤等[22]。其中流體體積函數VOF(volume of fluid)能夠追蹤兩相界面移動變化,且復雜度小,精度高[23],被認為適用于微流體多相流的模擬。

VOF 模型關于多相流模擬計算的控制方程如下所示。

連續性方程如式(3)。

動量方程如式(4)。

動量方程中動量源項為界面張力FSV,Fluent中界面張力模型是連續界面力(continuum surface force,CSF)模型如式(5)。

式中,ρ為流體密度,kg/m3;t為時間,s;u為速度矢量,m/s;μ為流體黏度,Pa·s;p為壓力,Pa;T為溫度,K;FSV為界面張力,N/m;σ為界面張力系數;κ為界面曲率;δ為狄拉克函數;r為徑向坐標;n為單位法向量。

由于Fluent Database 中沒有十二烷基苯和二氯乙烷的相關物性,因此通過實驗測量了十二烷基苯和二氯乙烷的黏度,設置物性見表5。其余組分物理參數使用Fluent Database中給定的值。

表5 物性參數

3.3 邊界條件

使用Fluent 模擬T形微通道混合器中的穩態對撞流。流體區域使用ANSYS Mesh 中的狹縫與曲率(proximity and curvature)功能進行網格劃分,通過控制最大尺寸控制網格質量。模擬過程采用穩態求解,動量方程采用二階迎風格式,流場采用PISO壓力速度耦合的方法求解,壓力項關聯采用交錯壓力格式(PRESTO!),壁面無滑移,通過調整松弛因子使計算達到收斂,當各參數誤差值小于10-3認為收斂,未提及的設置保持默認。在軟件中,設置兩相間的界面張力以規定界面張力動量源項,設置水與甲苯的界面張力為0.0357。兩個入口皆為速度入口,出口為壓力出口,表壓為0。

3.4 網格無關性

如圖5 所示,通過改變最大尺寸控制網格質量,發現當網格數量由300000 增加至1000000 時,相同邊界條件下出口最大速度變化趨勢如圖5。發現范圍內網格數量對模擬結果幾乎無影響,為減小計算量并保證結果準確,最終選擇網格數量為433960的網格進行模擬。

圖5 網格無關性

3.5 入口雷諾數計算

計算兩相液體為水和甲苯、溶劑為乙醇時的入口雷諾數,混合液體的密度由式(6)計算。

式中,ρm為混合流體密度,kg/m3;ρi為組分i的密度,kg/m3;φi為組分i的體積分數。

混合液體黏度的計算使用Lobe 方法,因為兩股混合液體分別由兩種組分組成,且不考慮溫度影響,所以簡化為式(7)。

其中兩個參數如式(8)、式(9)。

運動黏度與動力黏度關系如式(10)。

雷諾數計算如式(11)。

式中,純液體黏度值較小的組分為A,另一黏度值較大的組分為B;α為混合物中組分的特性黏度參數,J/mol;ν為組分運動黏度,mm2/s;μ為液體動力黏度,mPa·s。

入口雷諾數計算見表6。

表6 入口雷諾數

3.6 密度分布云圖

流體中3 種組分的密度分別為:水998kg/m3,甲苯866kg/m3,乙醇790kg/m3。因此通過密度分布云圖可以看出流體混合效果。無溶劑條件下密度分布云圖如圖6所示,以乙醇為溶劑條件下密度分布云圖如圖7所示。

圖6 無溶劑條件下密度分布云圖

圖7 以乙醇為溶劑條件下密度分布云圖

從圖6 可以看出,當入口速度為0.12m/s 時(對應流量為6mL/min,入口雷諾數分別為119.4和177.34),由于兩種組分溶解性低,流體的雷諾數較小,對撞效果不明顯,T形對撞處流體呈現明顯的兩種密度分布,意味著流體混合效果不強。但隨著進入管道,體積減小,壓強增加,流體的分層在一段距離后消失,兩組分發生了混合,流體密度變化。隨后因為即將流出管道,壓力減小,兩種組分又分離,形成清晰的彈狀流(slug flow),也稱為泰勒流(Taylor flow),是微通道內低流速條件下的一種兩相流動形式,流型如圖6中的兩相以一連串子彈形狀周期性分布。而當流速到達0.15m/s 時(對應流量為7mL/min,入口雷諾數分別為149.25 和221.67),其他現象與之前相似,但中間的流體分層不再消失,當流速為0.17m/s(對應流量為8mL/min,入口雷諾數分別為169.15和251.23)時也表現出相似的現象。而當流速增加到0.19m/s(對應流量為9mL/min,入口雷諾數分別為199 和295.56) 和0.21m/s 時(對應流量為10mL/min,入口雷諾數分別為208.95 和310.34),中間再次出現這一現象。分析認為是由于雷諾數的變化,管道內慣性力與黏性力主導地位轉變,低雷諾數時(入口雷諾數范圍119.4~149.25 和177.34~221.67),由于對撞效果不強,所以對撞處產生了清晰的分層,但之后界面張力使得混合效果增強,分層消失;中雷諾數時(入口雷諾數范圍149.25~169.15和221.67~251.23),對撞效果依然不強,但相對較短的停留時間使得界面張力不能長時間作用在分層界面處,因此不再出現分層消失現象而是直接產生泰勒流;高雷諾數時(入口雷諾數范圍169.15~208.95和221.67~310.34),對撞效果強烈,但界面張力作用時間不足,因此出現短暫混合后,流體立即變為彈狀流。

從圖7可以發現,溶劑的引入可以很好地提升混合效果,流體在T形對撞處不再形成明顯的分層。進入管道后液體立刻混合,一段時間后發生變化。當流速為0.12m/s、0.15m/s、0.17m/s 時(對應入口雷諾數分別為89.46 和113.32、111.82 和141.64、126.73和160.53),進入管道后流體先均勻混合,后逐漸分離,但形成尺寸小于管道尺寸的不規則的液滴,當流速增加到0.19m/s和0.21m/s時(對應入口雷諾數分別為149.09 和188.86、156.55 和198.30),流體進入管道后逐漸分離,形成規則的以密度從大到小排布在管道由內到外的三層流體。

水和甲苯這兩種互不相溶液體在引入了溶劑后,流態發生明顯改變,由原本的彈狀流轉變為分布更加均勻的流型。從密度分布云圖來看,處于代表三種液體密度值(水998kg/m3,甲苯866kg/m3,乙醇790kg/m3)中間區域的密度分布明顯增多,這代表了三種組分處于混合狀態的占比增高。從流型來看,流型由彈狀流轉變為特征尺寸小于管道的不規則液滴,相接觸面積增加,增強了混合強度。圖7中流型呈現出不同的表現形式,認為是隨著雷諾數的增大,流體受力向慣性力為主導轉變,同時隨著溶劑的引入,減小了兩相的界面張力作用,因此更加加強了慣性力對流型的影響,呈現為形成特征尺寸小于管道尺寸的液滴,在經過了碰撞區域一段距離之后,由于液滴間界面張力的作用,小液滴再次分開。由于流速的增大,慣性力使得密度更高的液滴有更高的速度而分布在管道更內測,液滴也由密度從大到小由內而外地分布在管道內,而低流速時就沒有這種情況。

3.7 壓力分布云圖

壓力分布云圖如圖8 所示。從圖中可以看出,流速為0.21m/s時,有溶劑條件下和無溶劑條件下管道內壓力都是從入口到出口遞減,區別在于有溶劑條件下壓力的變化更規律,分析是因為無溶劑條件下形成了明顯的泰勒流,導致管內壓力分布不均。

圖8 流速為0.21m/s時壓力分布云圖

3.8 DDB與SO3混合狀態表征

為表征實驗過程中微混合器內DDB與SO3的混合狀態,模擬了不同流量下體積分數為0.2DDB+0.8 二氯乙烷和0.2SO3+0.8 二氯乙烷的混合狀態,所得DDB體積分數分布云圖和SO3體積分數分布云圖如圖9所示。

從圖9 中可以看出,DDB 和SO3進入T 形混合器后只分布在管道的各自入口側,這一現象與實驗所得結論吻合,即過程受傳質效果影響。且單側入口流量為6mL/min、7mL/min、8mL/min時,DDB和SO3在進入T 形混合器一段距離后體積分數分布出現局部升高,而按照入口的體積分數計算,完全混合后DDB 和SO3的體積分數應同為0.1,這種現象說明進入T形混合器一段距離后DDB和SO3有一部分處于混合不完全的狀態。而隨著流速的增大,單側入口流量為9mL/min、10mL/min時,這一現象有所緩解(即圖中紅色區域變少),也與實驗所得當流速增大時混合效果變強的結論符合。認為是混合器內流體受力情況受入口流速大小影響,隨著流速增大,流體受到的慣性力更強,不同相之間碰撞更強烈,相對減少了表面張力的作用,也就增強了混合效果。

圖9 二氯乙烷為溶劑條件下DDB和SO3體積分數分布云圖

4 結論

通過響應面實驗,研究了流量、溫度、摩爾比對以二氯乙烷為溶劑條件下微通道反應器內十二烷基苯磺化過程的交互影響。通過CFD 計算流體力學方法對T形微通道混合器內流體流動和傳質現象進行表征,模擬了與實驗條件一致的以二氯乙烷為溶劑時DDB和SO3的混合狀態,研究了不同流速下T形混合器內的混合狀態。以水和甲苯作為互不相溶的兩相流體,以乙醇作為溶劑,研究了有無溶劑條件下不同流速對流體混合效果的影響。

(1)微通道反應器內十二烷基苯磺化過程中,產品活性物質質量分數在流量下受溫度和摩爾比的影響不大,隨著流量增加到8~10mL/min 范圍內,產品活性物質質量分數提高的同時,受到溫度和摩爾比的影響更明顯,且當溫度提升到60~70℃范圍時,由于SO3氣體的溢出,增強了管道內流體的混合效果,此溫度下低流量條件下的產品活性物質質量分數也有一定升高。認為產品活性物質質量分數主要受到傳質效果的影響,當傳質效果達到一定強度時這一過程開始向反應控制轉變。

(2)發現低流速下T形混合器內確實存在DDB和SO3的混合強度較弱的問題,DDB和SO3在進入T形混合器一段距離后,分別分布在管道的各自入口側。而隨著流速的增大,DDB和SO3完全分離的部分減少,相互混合的部分增多,與實驗所得結論吻合。

(3)發現由于流速同時影響T形碰撞區的碰撞強度和流體在管道內的停留時間,因此在界面張力為主導的微通道內,碰撞強度和界面張力的作用時間都能夠強化混合效果,而溶劑的引入能夠很好地強化兩種互不相溶流體的混合效果和管道內分布情況,且混合效果的變化與響應面實驗結果的討論吻合,混合效果隨流量的增大呈現先減弱后加強的變化趨勢。研究了微通道內十二烷基苯磺化過程不同因素的交互作用,對磺化過程的影響因素有了進一步的認識,對磺化工藝的優化有重要意義。同時通過CFD計算流體力學方法對T形微通道混合器的模擬,發現了溶劑條件對流體流態和混合機理的影響以及不同流速下流體混合狀態變化的規律。

主站蜘蛛池模板: 中文字幕在线观| 国产精品性| 久久人体视频| 91久久天天躁狠狠躁夜夜| 欧美精品啪啪一区二区三区| 伊人国产无码高清视频| 国产99视频精品免费观看9e| 国产 在线视频无码| 亚洲色大成网站www国产| 国产二级毛片| 亚洲一区色| 一区二区欧美日韩高清免费 | 在线欧美日韩国产| 思思热精品在线8| 亚洲精品日产精品乱码不卡| 青草精品视频| 国产v精品成人免费视频71pao| 黄色在线不卡| 一区二区三区在线不卡免费| 国产综合精品日本亚洲777| a毛片免费在线观看| 欧美有码在线| 精品久久高清| 久久综合干| 在线看片中文字幕| 午夜精品影院| 综合色婷婷| 国产精品第一区在线观看| 免费毛片a| 中文字幕久久波多野结衣| 综合人妻久久一区二区精品 | 久久免费看片| 无码人中文字幕| 青青青国产精品国产精品美女| 亚洲无线国产观看| 国产区成人精品视频| 欧美亚洲欧美| 亚洲欧美国产五月天综合| 国产成人精品免费视频大全五级| 日韩精品久久久久久久电影蜜臀| 992tv国产人成在线观看| 一本大道无码高清| 欧美亚洲香蕉| 手机精品视频在线观看免费| 婷婷五月在线视频| 国产精品美女免费视频大全 | 高清码无在线看| 无码啪啪精品天堂浪潮av| 中日无码在线观看| 午夜福利视频一区| 久久亚洲黄色视频| 国产成人盗摄精品| 激情视频综合网| 色综合五月| 五月天婷婷网亚洲综合在线| 久久国产精品嫖妓| 久久免费视频播放| 国产色婷婷| 亚洲精品桃花岛av在线| 91九色视频网| 国产精品久久久久久久久久98| 欧美日韩v| 日韩天堂网| 国产区网址| 亚洲精品无码抽插日韩| 熟妇丰满人妻| 国产日韩欧美精品区性色| 国产激爽大片在线播放| 国产无码精品在线| 伊人丁香五月天久久综合 | 九九视频在线免费观看| 熟妇无码人妻| 狠狠五月天中文字幕| 91久久大香线蕉| 亚洲国产欧美国产综合久久| 白丝美女办公室高潮喷水视频| 国产精品久久久久久搜索 | 欧美日韩国产在线播放| 欧洲精品视频在线观看| 欧美精品亚洲精品日韩专区| 女人av社区男人的天堂| 欧美国产日产一区二区|