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

SLAB模型參數優化及其應用研究

2015-08-30 03:22:12曹井國門曉曄楊宗政天津科技大學海洋科學與工程學院天津300457
安全與環境工程 2015年5期
關鍵詞:優化模型

熊 發,曹井國,門曉曄,楊宗政(天津科技大學海洋科學與工程學院,天津300457)

SLAB模型參數優化及其應用研究

熊 發,曹井國,門曉曄,楊宗政
(天津科技大學海洋科學與工程學院,天津300457)

SLAB模型是由典型的淺層模型演變而來的,適用于重氣泄漏擴散的模擬。為提高該模型的適用性和精度,在天津臨港經濟區內進行了二氟一氯甲烷現場泄放試驗,試驗結果與理論分析顯示了較好的一致性;在此基礎上將目標函數與響應面最優化求解方法相結合,研究出一種新的模型參數優化方法,采用此方法對SLAB模型的大氣穩定度和地表粗糙度兩個環境參數進行優化;并以臨港經濟區為主體,應用優化后的SLAB模型對苯泄漏事故案例進行模擬,對事故危害區域進行劃分,分析不同季節典型氣候條件下事故的風險區域范圍變化。結果表明: SLAB模型的地表粗糙度優化結果為0.32 m,大氣穩定度等級為C;發生在各個季節的重氣泄漏事故都具有各自的特點,事故的影響范圍差異性較大;該模型的應用可以為相關事故的應急指揮工作提供有針對性的指導。

SLAB模型;重氣泄漏擴散模擬;參數優化;地表粗糙度;大氣穩定度;天津臨港經濟區

重氣在化工和石油化工行業應用相當常見,幾乎涉及生產、運輸、儲存等各個環節,例如苯、氯氣(Cl2)、硫化氫(H2S)、液化天然氣(LNG)等。重氣一旦發生泄漏,會由于分子質量大(如CO2)、突發化學變化、自身溫度低(如LNG)等原因形成重氣云團[1]。多數重質氣體屬于危險氣體,具有急性毒性、易燃易爆等特性,對人體和生態環境危害巨大。重氣泄漏擴散過程的相關研究開始于20世紀70年代[2],隨著毒性重氣泄漏事故的不斷增加,此方面的研究工作得到了更多關注。由于重氣泄漏試驗風險較大,用計算機模擬具有顯著優越性,國內外已經開發出大量重氣擴散的計算模型[3]。

SLAB模型是由美國能源部(United States Department of Energy)的勞倫斯利弗莫爾國家實驗室開發的,為適用于重氣泄漏擴散模擬的大氣擴散模型[4]。該模型是由典型的淺層模型演變而來,可對重氣瞬間泄放、短時間泄放、持續泄放以及泄放后的擴散過程進行數值模擬,并能夠處理地面蒸發池泄漏、水平噴射泄漏、垂直噴射泄漏和瞬間或短期蒸發池泄漏共計4種泄漏排放形式。該模型在使用時不考慮下墊面的地勢起伏,地表粗糙度水平一般小于氣云的垂直高度[5],并以守恒方程作為基本控制方程(包括質量守恒、動量守恒、能量守恒和組分守恒),通過穩態煙羽模式、瞬變煙團擴散模式以及兩者的組合來實現對各種不同泄漏形式的重氣擴散過程模擬[6]。SLAB模型的數學表達式和參數是影響模擬結果的兩個關鍵因素,尤其是模型參數的變化將會對模擬結果產生很大影響。因此,本文重點研究SLAB模型的參數優化方法,并以天津臨港經濟區為主體對該模型相關參數進行優化,進而探討該模型在化工園區的應用方式。

1 現場重氣泄放試驗

現場重氣泄放試驗在天津臨港經濟區內進行,泄放物質為二氟一氯甲烷(R22),試驗時間為2014 年11月3日上午,風向為西南偏南,相對濕度均值為17.6%,即時氣溫均值為18.4℃,風速均值為2.2 m/s,屬于二級風的范疇,大氣穩定度為C級,模擬儲罐泄漏事故發生后重質氣云的擴散過程,以為模型的參數優化提供現場觀測數據。

1.1試驗方法

試驗主要通過采集泄放開始后不同時間點的下風向各點位處的氣體樣品,并分析氣體樣品中的R22含量,其試驗方法如下:

1.1.1現場布置

根據現場風向,確定空地的上風向某點為泄放點,在泄放點處設置支架,將成品罐裝R22放置于支架上,調整泄放口高度為1.5 m,并以泄放點為起點,分別在下風向5 m、10 m、20 m、50 m、100 m處設置氣體樣品采集點。具體采樣點布設情況見圖1。

1.1.2樣品采集

分別在泄放過程開始后的5 min、10 min、20 min、30 min4個時間節點,對5個采樣點位處的氣體樣品進行采集,樣品采集點距離地面高度一律為0.5 m,泄放持續時間為15 min。

圖1 采樣點布設示意圖Eig.1 Schematic diagram of the sampling points

1.1.3樣品檢測

采用氣相色譜頂空進樣與ECD捕集連用的方法對氣體樣品進行檢測。首先用純R22氣體進樣,找到其出峰位置,并標定其峰面積與R22含量的對應關系,作出標準曲線;然后通過測定每個樣品的峰面積,即可得出其R22濃度。

1.2試驗結果及分析

各個采樣點位處的R22濃度隨時間的變化情況見圖2。

圖2 各個采樣點位處R22濃度隨時間的變化曲線Eig.2 Variation curves of R22 concentration at all sampling sites with time

由圖2可以看出:各個采樣點位處的R22濃度變化均為先升高后降低,曲線始末兩端均無限趨近于零,且在泄放過程持續階段變化趨勢較為平緩。這是由于氣體初始釋放時濃度較低,隨后慢慢富集形成重氣云團,重氣云團在重力作用下沉降,并在氣流作用下緩慢移動,與后續補充的重氣云團構成較為穩定的狀態,使固定點處氣體濃度趨于穩定;泄放過程停止后,缺少了后續氣體的補充,重氣云團在風力作用下逐漸向遠處擴散稀釋,最終濃度趨近于零。可見,各個采樣點位的R22濃度變化過程與氣體的泄放擴散過程基本相符,能夠反映氣體泄放—重力沉降—擴散稀釋的整個過程。

2 模型參數優化

模型參數的優化即是尋求一組參數,使模型的輸出與實際觀測數據之間按給定目標函數的度量方式達到最佳擬合[7],即通過模型參數識別、參數取值范圍確定、實測值獲取、目標函數確定等一系列過程完成對各個參數的率定和優選,從而提高模型的精確性,使其在實際應用中能夠更好地發揮作用。重氣泄漏擴散模型的參數一般包括物質理化參數、泄漏源參數、環境參數和系統參數四大類[4]。本文以天津臨港經濟區為主體,從參數是否可精確獲取和對模型預測結果的影響程度兩個方面考慮,對SLAB模型的大氣穩定度(Atmospheric Stability Class)和地表粗糙度(Surface Roughness Height)兩個環境參數進行優化。

2.1模型參數優化方法

本文對SLAB模型參數優化的具體方法是通過構建以模型實測值和預測值差值為基礎的目標函數,利用Design Export 8的試驗設計功能,計算參數不同取值下的目標函數值,分析參數取值與目標函數值之間的對應關系,并通過其最優化求解功能給出令目標函數值接近于零的最優化方案。具體步驟如下:

(1)目標函數構建。構建的目標函數為

式中:Fi0為實測值;Fi為模型預測值;x,y,…代表物質理化性質、泄漏源參數、氣象參數等客觀參數;θ1,θ2分別代表兩個需要優化的模型參數。

(2)參數取值范圍確定。地表粗糙度是從空氣動力學的角度出發,描述風速零點距離地表的平均高度值,反映的是地表對于風速減弱作用的影響,其大小取決于地表粗糙元的性質及流經地表的流體的性質[8]。經過長期研究,國內外學者已經得出多種確定不同地表類型的地表粗糙度的方法,并在分析總結大量試驗數據的基礎上,得出了各種地表類型不同條件下的地表粗糙度參考值[9]。化工園區內建筑物分布較為密集,但建筑物高度比城市建筑物低,同時有些園區內存在大量空地,其地表粗糙度取值應介于城市類型和草原類型之間,但考慮到多數化工園區為了運輸方便選址位于沿海地帶,地表粗糙度的取值會受其影響,因此本文將化工園區地表粗糙度參數的取值范圍確定為0.2~1 m。大氣穩定度指近地層大氣作垂直運動的強弱程度,是影響污染物在大氣中擴散的極重要因素。根據經典的Pasqiull穩定度劃分方法[10],可將大氣穩定度分為A(極不穩定)、B(中等不穩定)、C(弱不穩定)、D(中性)、E(弱穩定)、E(穩定)6種穩定度狀況。擴散問題應用領域的大氣穩定度通常以理查森數(Ri)或與之相聯系的參數為基礎所建立的穩定度分類法來劃分[11]。范紹佳等[12]給出了沿海地區理查森數Ri大氣穩定度推薦分類標準,據此分類標準對天津地區近一年內的大氣穩定度等級進行頻數統計,結果顯示天津地區的大氣穩定度等級按頻率由高到低排列依次為:C(39.72%)>D(35.16%)>B(10.24%)。由此確定出大氣穩定度等級優化前的初始取值范圍是B~D,經數字化處理后為2.0~4.0。

(3)實測值和預測值獲取。本文的實測值是通過現場重氣泄放試驗獲取;預測值是通過SLAB模型計算得到。SLAB模型所需泄漏源參數和氣象參數與現場試驗環境保持一致,泄放物質的理化參數通過查閱相關資料獲得。

(4)制定模擬試驗方案。根據步驟(2)中確定的參數取值范圍,以現場重氣泄放試驗中B點在泄放開始后10 min時的實測數據為因變量,用Design Export 8的CCD模塊設計出9組兩因素多水平正交試驗,具體參數優化試驗方案見表1。

表1 參數優化試驗方案及結果Table 1 Scheme and result of the parameter optimization experiment

(5)回歸方程構建。應用Design Export 8的數據擬合功能,構建目標函數值關于待優化參數的回歸方程,并根據回歸方程的方差檢驗結果,得到精簡后的回歸方程。

(6)計算回歸方程最優解。利用Design Export 8的最優化求解功能,計算目標函數絕對值最小的點對應的參數值,即可得出模型的參數優化結果。

2.2模型參數優化結果

模型參數優化結果見表1。以大氣穩定度(S)和地表粗糙度(R)為自變量,目標函數值為因變量(P),建立目標函數值關于待優化參數的回歸方程,并對該回歸方程進行方差分析,檢驗回歸方程系數的顯著性,去除回歸方程中的不顯著項,得到如下精簡后的回歸方程:

令式(2)取值為零,利用Design Export 8的方程求解功能求得模型參數最優解為:地表粗糙度=0.32 m,大氣穩定度=2.51,目標函數Min= -0.000 124 211。根據實際情況,確定出的SLAB模型優化后參數為地表粗糙度為0.32 m,大氣穩定度等級為C,作為SLAB模型在化工園區內使用時的基礎參數值。

3 SLAB模型在天津臨港工業區的應用

3.1區域概況

天津臨港經濟區位于海河入海口南側灘涂淺海區,是通過圍海造地而形成的港口工業一體化的海上工業新城,地處中緯度歐亞大陸東岸,面對太平洋,季風環流影響顯著:冬季受蒙古冷高氣壓控制,盛行偏北風;夏季受西太平洋副熱帶高氣壓左右,多偏南風。該地區屬暖溫帶半濕潤大陸季風型氣候,有明顯由陸到海的過渡特點,年平均氣溫為14.4℃,7月最熱,月平均氣溫可達26℃,1月最冷,月平均氣溫為-4℃[13]。臨港經濟區內共包含化工生產類企業20余家,如天津大沽化工有限公司、樂金渤海化學有限公司、天津堿廠等,通過調研,該經濟區內及其周邊環境敏感點共有4個:臨港經濟區管委會、臨港生態餐廳、散貨交易中心和大沽炮臺。各環境敏感點的具體位置詳見圖3。

圖3 臨港經濟區環境敏感點的位置Eig.3 Location of each environmental sensitive point

3.2事故模擬

苯(Benzene,C6H6)是化工行業中最為常用的物料之一,是合成苯酚、苯乙烯、鹵代苯等芳香族化合物的重要原料,也是一種可燃、有毒、致癌的危險化學品。苯在常溫下為一種無色、有甜味的透明液體,其蒸氣密度比空氣重,屬于重氣的范疇。化工園區內關于苯突發環境污染事故較為多見,事故類型以罐體泄漏為主[14]。本文假設事故案例發生地為臨港經濟區內某化工廠,事故類型為液態苯儲罐泄漏,環境參數按照臨港經濟區全年氣象數據統計結果確定。苯泄漏事故的參數值詳見表2。

表2 苯泄漏事故的參數值Table 2 Parameters of Benzene leakage

3.2.1云團運動軌跡分析

參照上文中確定的模型優化參數,用修正后的SLAB模型對苯泄漏后的擴散過程進行模擬,分別截取泄放開始后1 160 s、1 970 s、2 420 s和3 090 s 4個時刻,對苯蒸汽云團的運動軌跡以及隨時間的變化規律進行數值模擬,其模擬結果見圖4。

圖4 不同時刻苯濃度空間分布圖Eig.4 Spatial distribution of Benzene concentration at different moments

由圖4可以看出,苯泄漏后擴散的過程較為緩慢,這是由于液態苯經水平噴射過程泄漏后揮發成苯蒸汽,苯蒸汽在重力作用下沉降,并在地面附近大量富集并形成重氣云團;之后重氣云團在氣流的驅動力作用下向下風向移動,并在擴散過程中被不斷稀釋,使得重氣云團向側風向擴張,同時其濃度降低。此次模擬事故中,苯在泄漏后1 160 s時刻到達臨港經濟區管委會附近,之后越過臨港經濟區西北邊界線,繼續向下風向擴散;隨后在1 970 s時刻到達大沽炮臺,但此時重氣云團的濃度已經較低;在泄漏過程停止后(1 800 s),近源端的污染物濃度急速下降,污染區域的擴張也進一步放緩;最后空氣中的重氣云團在氣流的作用下進一步稀釋擴散,在3 090 s時刻重氣云團移動軌跡中間出現低濃度區間,近源端的污染物也已擴散殆盡。由于在泄放持續進行過程中,氣云在穩態煙羽擴散模式下擴散,氣云中污染物的量基于下風向距離(x)變量守恒,氣云前端仍處于側風向擴張的階段,所以此階段氣云的形狀呈扇形[如圖4(a)];當泄放過程停止后,氣云的擴散模式變為瞬變煙團擴散,氣云中污染物的總量保持不變,氣云在風力作用下向下風向移動,同時由于大氣的湍流卷吸作用,污染物濃度持續下降,導致氣云由扇形向帶狀轉變[見圖4(b)至(d)]。

3.2.2事故影響范圍分析

苯的急性毒性為中等水平,因此本文選用適用于急性毒性水平不高的化學污染事故緊急響應計劃指南(ERPG)作為風險區域確定的標準。通過查閱美國環境保護署(EPA)官方網站[15],得到苯對應的風險等級ERPG-1、ERPG-2和ERPG-3分別為50 ppm、150 ppm和1 000 ppm,并應用SLAB模型對假設的苯儲罐泄漏事故在不同季節條件下進行數值模擬,在其他條件不變的情況下,得到苯泄漏事故不同風險等級的危害縱深,見表3。

表3 不同季節條件下苯泄漏事故不同風險等級的危害縱深Table 3 Comparison of the hazard depth of Benzene leakage in different seasons

由表3可以看出:季節變化引起的苯泄漏事故危害縱深波動較為明顯,不同風險級別對應的危害縱深之間差別很大,但變化規律較為一致:ERPG-1對應的不同季節條件下苯泄漏事故危害縱深在3 660~5 164 m范圍內變化,ERPG-2對應的危害縱深變化范圍為1 965~2 650 m,ERPG-3對應的危害縱深變化范圍較小,為384~652 m。4個季節按苯泄漏事故危害縱深從大到小排序為:冬季>秋季>春季>夏季。

臨港地區春季以東南風為主導風向,風速與年均值相比較大,發生在春季的泄漏擴散事故其影響范圍分布更為狹長,污染范圍擴張速度更快,且下風向敏感點較多,應做好事故的防范與應急準備工作;夏季主導風向同樣是以東南風為主,風速與春季相比較小,高溫高濕是夏季的一大特點,當有重氣泄漏擴散事故發生時,污染物在高溫環境中分子活性更大,加快了其蒸發速率以及側風向和鉛垂風向上的運動速率,從而縮短了其事故影響距離,但包括苯在內的多數重質氣體都具有易燃易爆性,夏季的高溫特性在增加了其擴散稀釋速率的同時也增加了污染物燃燒爆炸的可能性,因此在事故應急過程中,應注意防止引發燃爆事故等次生危害;秋季的氣候特點與春季類似,其主導風向較春季略微偏東,如發生危險品突發泄漏事故更易對臨港管委會和散貨交易中心等環境敏感點產生直接影響;冬季是氣溫最低的時期,臨港地區冬季的平均氣溫在0℃左右,主導風向為西北風,冬季的低溫環境造成重氣氣體泄漏后擴張較為緩慢,能給相關人員留出更多的應急響應時間,但其影響距離卻是幾個季節發生的同類型事故中最大的,雖然主導風向是朝向渤海的,但冬季經常性的大氣逆溫現象會為污染物的擴散稀釋造成極大困難。

綜上分析可見,臨港地區發生在各個季節的重氣泄漏事故都具有各自的特點,事故的影響范圍差異性較大。

4 結 論

本文在天津臨港經濟區以二氟一氯甲烷(R22)為泄放介質,進行了現場重氣泄放試驗,并應用氣相色譜法對不同時刻下風向不同距離點位處采集的氣體樣品進行目標氣體濃度分析,得到了R22泄放擴散過程濃度變化曲線,其與理論分析結果具有較好的一致性;在此基礎上研究出了一種新的模型參數優化方法,并采用此方法對大氣穩定度和地表粗糙度兩個關鍵參數進行優化,優化結果為地表粗糙度為0.32 m,大氣穩定度等級為C;通過將優化后的SLAB模型在臨港經濟區進行模擬應用,達到了預期效果。此外,該模型還可以與各種環境應急監測[16]手段相結合,為相關事故的應急指揮工作提供輔助決策信息。

[1]張啟平,麻德賢.危險物泄漏擴散過程的重氣效應[J].北京化工大學學報,1998,25(3):86-90.

[2]Britter R E.Atmospheric dispersion of dense gas[J].Ann.Rev. Fluid Mech.,1989,21:317-344.

[3]黃文宏,章保東,包其富,等.重質氣體泄漏擴散模型研究綜述[J].浙江化工,2009,40(7):18-22.

[4]Ermak D L.User's Manual for SLAB:An Atmospheric Dispersion Model for Denser-Than-Air Releases[M].USA:Lakes-Environmental,1997.

[5]孫召賓,鄭洪波,張樹深.重氣體擴散模型分類及工業應用模型分析比較[J].中國人口資源與環境,2011,21(1):435-438.

[6]孫召賓.危險化學品泄放事故后果計算模型的研究及應用[D].大連:大連理工大學:2012.

[7]劉毅,陳吉寧,杜鵬飛.環境模型參數優化方法的比較[J].環境科學,2002,23(2):1-6.

[8]周艷蓮,孫曉敏,朱治林.幾種典型地表粗糙度計算方法的比較研究[J].地理研究,2007,26(5):887-896.

[9]布和敖斯爾.全球環境變化研究的區域模式與遙感和GIS方法[M].呼和浩特:內蒙古教育出版社,2005:239.

[10]Pasquill E.The estimation of the dispersion of windborne material[J].Meteo.Mag.,1961,90:33-49.

[11]Mohan M,Siddiqui T A.Analysis of various schemes for the estimation of atmospheric stability classification[J].Atmospheric Environment,1998,21(32):3775-3781.

[12]范紹佳,林文實,蘇雄暉,等.理查森數Ri在沿海近地層大氣穩定度分類中的應用[J].熱帶氣象學報,1999,15(4):370-375.

[13]天津臨港經濟區管理委員會.臨港簡介[EB/OL].[2014-1-3]. http://lg.ykdz.gov.cn/linganggaikuang/linganggongyequjianjie/.

[14]杜紅巖,王延平,盧均臣.2012年國內外石油化工行業事故統計分析[J].中國安全生產科學技術,2013,9(6):184-188.

[15]NIOSH.Immediately Dangerous To Life or Health[EB/OL]. [2014-5-7].NIOSH Publications and Products,http://www.cdc.gov/niosh/idlh/idlhintr.html.

[16]趙文喜,陳素寧.突發性環境污染事故的環境風險評價與應急監測[J].安全與環境工程,2009,16(4):14-17.

Parameter Optimization of SLAB Model and Its Application

XIONG Ea,CAO Jingguo,MEN Xiaoye,YANG Zongzheng
(College of Marine Science and Engineering,Tianjin University of Science and Technology,Tianjin 300457,China)

SLAB model is evolved from the typical shallow layer model and suitable for heavy gas leakage and diffusion simulation.In order to improve the applicability and accuracy of the model,this paper conducts the monochlorodifluoromethane releasing experiment in Tianjin Harbor Economic Area and the result is quite consistent with the theory.Based on the result,the paper combines the objective function and the response surface optimization algorithm and develops a new model parameter optimization method.Then the paper applies the method to optimizing 2 parameters of the SLAB model,namely,the surface roughness and the level of atmospheric stability.In Tianjin Harbor Economic Area,the paper applies the optimized SLAB model to simulating the benzene leakage case,divides the accident hazard zones and analyzes the size change of these hazard zones under typical seasonal conditions.The results show that the surface roughness of SLAB model is 0.32 m,and the atmospheric stability class is C.The heavy gas leak accidents that occur in each season have their own characteristics with obvious differences in hazard zones under influence.So the model can provide corresponding guidance for emergency command of accidents.

SLAB model;heavy gas leakage diffusion simulation;parameter optimization;surface roughness height;atmospheric stability class;Tianjin Harbor Economic Area

X51;X937

A

10.13578/j.cnki.issn.1671-1556.2015.05.004

1671-1556(2015)05-0020-05

2015-01-05

2015-05-31

國家重大科學儀器專項(2012YQ060165)

熊 發(1989—),男,碩士研究生,主要研究方向為安全事故模擬及預測。E-mail:faxiong@126.com

楊宗政(1974—),男,博士,教授,主要從事三廢處理及資源化利用方面的研究。E-mail:yzz3520@163.com

猜你喜歡
優化模型
一半模型
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 成人在线观看不卡| 成人在线亚洲| 日韩欧美国产成人| 色偷偷一区| 狠狠色婷婷丁香综合久久韩国| 欧美亚洲综合免费精品高清在线观看| 青青草原国产| 国内熟女少妇一线天| 亚洲一区波多野结衣二区三区| 久久亚洲国产视频| 国产精品一老牛影视频| 最新无码专区超级碰碰碰| av在线手机播放| 亚洲美女久久| 免费高清毛片| 国产精品99r8在线观看| 国产色婷婷| 欧美成在线视频| 亚洲成在线观看| 狠狠综合久久久久综| 国产成人精品视频一区视频二区| 国产精品综合色区在线观看| 国产精品30p| 国产又黄又硬又粗| 丁香婷婷在线视频| 99视频精品在线观看| 久久五月天国产自| 四虎精品国产AV二区| 欧美一区二区丝袜高跟鞋| 国产在线欧美| 在线视频精品一区| 国产人人射| 国产综合在线观看视频| 3344在线观看无码| 久久综合激情网| 婷婷综合色| 国产jizz| 国产精品护士| 亚洲精品欧美日本中文字幕 | 亚洲av色吊丝无码| 亚洲精品国产成人7777| 99视频在线精品免费观看6| 四虎在线观看视频高清无码| 久久永久免费人妻精品| 久久成人18免费| 国产成人免费视频精品一区二区| 国产91小视频在线观看 | 91黄视频在线观看| 在线视频亚洲欧美| 精品国产黑色丝袜高跟鞋| 亚洲一区精品视频在线| 亚洲第一在线播放| 嫩草在线视频| 中文字幕无线码一区| 亚洲成aⅴ人片在线影院八| 在线国产欧美| 伊人久久综在合线亚洲2019| 亚洲日韩国产精品综合在线观看| 欧美色图第一页| 欧美午夜网站| 久久毛片网| 九九热精品视频在线| 亚洲另类国产欧美一区二区| 伊人成人在线| 日韩久久精品无码aV| 国产又色又刺激高潮免费看| 国产小视频免费| 色欲色欲久久综合网| 久久一本日韩精品中文字幕屁孩| AV天堂资源福利在线观看| 国产95在线 | 91精品人妻互换| 视频在线观看一区二区| 四虎免费视频网站| 国产成人狂喷潮在线观看2345| 国产欧美日韩专区发布| 日韩欧美中文在线| 毛片久久久| 中日韩一区二区三区中文免费视频 | 国产91色| 四虎永久在线精品影院| 欧美区一区|