亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        考慮非平衡效應的過冷水滴凝固特性

        2017-11-22 01:12:34肖光明杜雁霞王橋郭龍王茂
        航空學報 2017年2期
        關鍵詞:結(jié)冰水滴界面

        肖光明, 杜雁霞, 王橋, 郭龍, 王茂

        中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室, 綿陽 621000

        考慮非平衡效應的過冷水滴凝固特性

        肖光明, 杜雁霞*, 王橋, 郭龍, 王茂

        中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室, 綿陽 621000

        非平衡凝固是過冷條件下水滴凝固過程的重要現(xiàn)象。本文針對飛機結(jié)冰過程過冷水滴的非平衡凝固效應,發(fā)展了改進的凝固特性預測模型及數(shù)值計算方法,并自行搭建了實驗系統(tǒng),開展了所建過冷水滴凝固模型與數(shù)值預測方法的實驗驗證。研究表明,所發(fā)展的改進模型可有效表征水滴過冷階段的非平衡凝固效應,因而對冷水滴凝固速率的預測有較好的改進;當過冷度為0 ℃時,過冷模型退化為傳統(tǒng)模型?;谒ǚ椒ǎ_展了過冷度及冷卻條件對水滴凝固特性的影響分析,獲得了不同條件下水滴凝固過程的溫度分布及相界面變化特征。研究表明,過冷度越大或水滴尺度越小,凝固速率相對越高;在考慮非平衡凝固效應的條件下,過冷水滴凝固速率要高于不考慮非平衡凝固效應的工況。相關研究可為結(jié)冰熱力學模型的改進,以及結(jié)冰特性的精細化預測提供參考。

        過冷水滴; 結(jié)冰; 非平衡效應; 凝固; 改進模型

        飛機結(jié)冰是影響飛行安全的重要隱患之一,也是過冷水滴撞擊于飛機表面并發(fā)生凝固的一種特殊相變現(xiàn)象[1]。由于過冷條件的存在,因而飛機結(jié)冰具有顯著的非平衡凝固特征[2]。當水滴低于凝固點溫度以液態(tài)形式存在時,往往會形成亞穩(wěn)平衡態(tài)[3-4]。此時,只要施加一個較小的擾動即可觸發(fā)凝固,并使其回到穩(wěn)定的平衡態(tài)[5],而能量的波動、界面、雜質(zhì)、振動等均是觸發(fā)亞穩(wěn)態(tài)液體發(fā)生凝固的擾動源[5-6]。

        過冷水滴的凝固通常分為2個典型階段[2,6]。第1階段由形核開始,水滴從熱力學非平衡態(tài)過渡到熱力學平衡態(tài)的階段,即枝晶形成階段,也有研究者稱其為部分凝固階段[7]。對于撞擊于飛機表面的過冷水滴,由于機體表面提供了異相形核的條件,因此在此階段形核過程從界面逐步發(fā)展到整個水滴,使水滴由液相變成冰水共存的模糊相,水滴溫度也由過冷態(tài)上升到凝固溫度所處的平衡態(tài)[6-7]。第2階段為完全凝固階段,即液固相界面由固相向液相推進直至凝固完成的過程[7-8]。由于撞擊引起的異相形核作用,飛機結(jié)冰過程水滴凝固的第1階段顯著快于第2階段[2,6]。鑒于結(jié)冰物理過程的復雜性,目前大多關于過冷水滴結(jié)冰的研究均把重點放在凝固第2階段即相界面的推進過程上。由于第2階段以平衡凝固為主,因此,多數(shù)研究者基于平衡凝固的相關理論與方法開展了結(jié)冰特性的預測研究, 如基于Enthalpy-Porosity法研究凝固過程的液/固相變行為[9-10],但該方法主要針對平衡凝固過程,無法表征過冷水滴凝固的第1階段即非平衡凝固過程的影響特征。

        近年來,隨著飛機結(jié)冰預測精度要求的提高,結(jié)冰過程的非平衡凝固效應引起了研究者越來越多的關注。Worster[11]、Ellen[12]等提出了平面生長理論,針對飛機結(jié)冰凝固過程的2個階段,采用了不同的預測方法。Feuillebois等[6]研究了非平衡凝固對結(jié)冰第2階段的初始物性參數(shù)的影響特性。Blake等[2]在考慮形核過程的基礎上,基于FLUENT二次開發(fā)發(fā)展了過冷水滴凝固特性的預測方法。可以看出,由于過冷條件下水滴的非平衡凝固效應對后期結(jié)冰速率及結(jié)冰特征有著重要影響而受到了越來越多的關注[13-14]。但由于凝固過程的復雜性,非平衡凝固特性預測方法的相關研究目前還較為缺乏,預測精度也有待提高,因而非平衡凝固規(guī)律特征的相關研究也較為薄弱。本文針對過冷水滴結(jié)冰的非平衡凝固效應,發(fā)展了能表征非平衡凝固效應的過冷水滴凝固特性預測方法,并自行搭建了實驗系統(tǒng)開展所建方法的實驗驗證。相關研究可為結(jié)冰熱力學模型的改進,以及結(jié)冰特性的精細化預測提供參考。

        1 考慮非平衡效應凝固模型的建立

        對單個過冷水滴凝固特性的研究有助于深入揭示飛機結(jié)冰過程過冷水滴凝固的物理特性[3]。為便于觀測和實驗,本文以單個水滴為對象,研究過冷水滴結(jié)冰過程非平衡凝固現(xiàn)象的共性特征。凝固實驗研究表明,過冷水滴的凝固過程由枝晶形成和相界面推進2個典型階段構(gòu)成。在枝晶形成階段,水滴由透明態(tài)轉(zhuǎn)變?yōu)槟:龖B(tài),如圖1(a)所示;在相界面推動階段,固/液相界面由冷卻面即底面向頂部平行推進直至凝固完成,如圖1(b)~圖1(d)所示。

        為建立相應的數(shù)理模型,圖2顯示了過冷水滴凝固過程的簡化示意圖。

        在凝固初始階段即枝晶形成階段,過冷水滴溫度由過冷態(tài)上升至凝固點,并伴隨潛熱的部分釋放,是典型的非平衡凝固階段;凝固第2階段為由熱擴散驅(qū)動的界面推進階段,也是在等溫下進行的平衡凝固階段。盡管研究表明凝固第1階段相對于第2階段的時間較短,但由于結(jié)冰條件的不同,使第1階段結(jié)束時形成了第2階段的不同初始條件,從而影響了凝固的后續(xù)特征。因此,過冷水滴整個凝固過程的有效預測應綜合考慮第1階段的非平衡凝固效應及第2階段的平衡凝固效應。針對過冷水滴的凝固兩個典型階段的特點,如何在相界面推進過程的預測中考慮第1階段的非平衡效應,是有效預測結(jié)冰全過程凝固特性需解決的重要問題。

        圖1 過冷水滴凝固的典型階段
        Fig.1 Typical freezing stages of supercooled water droplet

        圖2 過冷水滴凝固過程的簡化示意圖
        Fig.2 Schematic of supercooled droplet freezing process

        考慮到在凝固第1階段完成并形成混合態(tài)的過程中過冷水滴已經(jīng)有了潛熱的部分釋放,因此,在凝固第2階段即相界面推進過程的預測中應考慮凝固第1階段的影響。將無量綱過冷度ε表示為

        ε=cpsΔT/L

        (1)

        式中:ΔT=Tf-Tsupercooled,Tf為凝固溫度,Tsupercooled為過冷態(tài)溫度;L為相變潛熱;cps為固相定壓比熱容。

        凝固第1階段結(jié)束時混合態(tài)中未發(fā)生相變的液相分數(shù)可表示為

        (2)

        式中:cpl為液相定壓比熱容。

        因此,混合態(tài)的液/固相變潛熱可表示為[6]

        Lmix=Lf1

        (3)

        混合態(tài)的物性參數(shù)可表示為

        ρmix=ρlf1+ρs1-f1

        (4)

        cpmix=cplf1+cps1-f1

        (5)

        λmix=λlf1+λs1-f1

        (6)

        式中:ρ、cp和λ分別為密度、定壓比熱容和熱導率;下標mix、l和s則分別對應混合態(tài)液相和固相物性參數(shù)。表1分別給出了液相水和固相冰對應的物性參數(shù)[6]。

        表1 冰和水的材料物性參數(shù)[6]Table 1 Physical properties of ice and water[6]

        當獲得凝固第1階段凝固過程的物性參數(shù)后,包含非平衡凝固的相變問題即可轉(zhuǎn)化為平衡凝固問題。為借鑒平衡凝固條件下液/固相變傳熱的預測方法,作以下幾點假設:

        1) 由于撞擊形成的異相形核作用,凝固第1階段時間顯著小于第2階段。

        2) 忽略水滴凝固過程因體積膨脹引起的水滴變形。

        3) 界面推進過程的液相區(qū)實質(zhì)為已在凝固第1階段發(fā)生了部分相變的混合態(tài),因此凝固第2階段液相區(qū)物性參數(shù)由混合態(tài)參數(shù)代替。

        基于上述假設并借鑒Enthalpy-Porosity方法,則凝固過程的傳熱控制方程可描述為

        (7)

        ρgβT-Tref+S

        (8)

        (9)

        式中:u為液體流動速度矢量;p為液體壓力;h為熱焓;μ為液體黏性系數(shù);g為重力加速度矢量;β為液體熱膨脹系數(shù),且滿足Boussinesq近似;Tref為參考溫度;矢量S為凝固第2階段兩相區(qū)的源項,可描述為[10]

        S=-C1-f2u

        (10)

        其中:C為界面推進過程兩相共存區(qū)的特征參數(shù);f2為結(jié)冰過程第2階段未凝固的液相分數(shù),可表示為

        (11)

        其中:Tl和Ts分別為液/固相變過程完全融化和完全凝固時的溫度。

        2 過冷水滴凝固特性的數(shù)值分析

        圖3顯示了直徑為4 mm、高度為2.5 mm,過冷度為-10 ℃條件下水滴凝固過程的溫度分布及相區(qū)分布特性預測結(jié)果,t為水滴凝固時間??梢钥闯觯谙嘟缑嫱七M的初始階段,在熱擴散的作用下,相界面由固相區(qū)向液相區(qū)平行推進。隨著凝固過程的發(fā)展,相界面與等溫線逐步向水滴頂部彎曲。

        圖3 水滴凝固過程的溫度場及相界面變化特性
        Fig.3 Temperature and interface evolution characteristics of droplet freezing process

        圖4 過冷度對相界面變化速率的影響特性
        Fig.4 Effect characteristics of supercooled temperature on freezing rate of phase interface

        圖4顯示了過冷度分別為-15、-10、-5、0 ℃條件下相界面變化速率的比較(縱坐標為相界面推進高度ΔH與水滴半徑R的比值)。可以看出,在相同冷卻面溫度條件下,隨著水滴過冷度的增大,液/固相變的驅(qū)動力增加,相界面移動速度相應增大;當過冷度為0 ℃時,過冷模型退化為傳統(tǒng)Enthalpy-Porosity模型。

        3 實驗驗證

        為驗證所發(fā)展方法的有效性,本項目自行搭建了實驗系統(tǒng)并開展了相應算法的實驗研究。實驗系統(tǒng)如圖5所示,由半導體制冷系統(tǒng)、Agilent 34970A多點數(shù)據(jù)采集儀、ANV TF100 PID溫度監(jiān)視器、溫度控制器、MotionXtra HG-100K高速攝像機、FLIR E60紅外熱像儀、LED無影光源系統(tǒng)、工控機及調(diào)壓電源組成。實驗時,采用滴管在冷表面產(chǎn)生不同尺度的水滴;開啟半導體制冷裝置,由PID溫度控制系統(tǒng)將制冷裝置冷卻面溫度維持在0 ℃以下實驗所需的溫度條件,使水滴緩慢冷卻至過冷態(tài)直至水滴完全凝固?;诮缑孀粉櫟姆椒?,采用高速攝像機記錄獲得水滴凝固過程相界面隨時間的變化特性;同時,采用紅外熱像儀記錄水滴凝固過程溫度的實時變化特性。每組工況重復3次,取每組平均值作最終實驗值。

        在水滴凝固實驗中,在過冷期間的形核與枝晶生長階段,水滴逐步由透明態(tài)向模糊態(tài)過渡;當水滴達到平衡溫度,固/液界面逐步出現(xiàn)并由液相區(qū)向固相區(qū)移動。如果以固/液界面出現(xiàn)的起始時刻作為凝固第2階段的起點,用高速相機即可記錄相界面隨時間的變化特征。圖6顯示了圖3相同條件下過冷水滴凝固第2階段相界面隨時間的變化特性??梢钥闯觯嘟缑孀钕仍诶鋮s面產(chǎn)生,并逐步由固相區(qū)向液相區(qū)平行推進,與計算相界面移動特性相似。

        圖5 水滴凝固過程的實驗系統(tǒng)
        Fig.5 Experimental setup for droplet freezing

        圖6 實驗凝固過程相界面隨時間的變化特征
        Fig.6 Experimental solidification evolution characteristics of phase interface vs time

        圖7為本文發(fā)展的數(shù)值模型與傳統(tǒng)模型獲得的界面變化計算結(jié)果與實驗結(jié)果的比較??梢钥闯?,在相同條件下,由于考慮了過冷水滴在過冷階段非平衡凝固效應的影響,改進的預測模型的計算與實驗結(jié)果的吻合程度優(yōu)于傳統(tǒng)Enthalpy-Porosity預測模型,且界面推進速率要快于不考慮非平衡效應的工況。改進后的模型與實驗結(jié)果吻合得更好,表明了所建模型的有效性。

        圖8顯示了半徑R分別為2 mm和3 mm兩種不同尺寸水滴,在過冷度為-10 ℃條件下的計算與實驗相界面移動特性的比較。從計算與實驗結(jié)果的比較可以看出,總體而言,計算與實驗相界面隨時間的移動特性吻合較好,但由于計算中忽略的凝固過程水滴體積膨脹效應引起的變形及冒尖現(xiàn)象在凝固后期更為顯著,因而在凝固前期計算相界面與實驗相界面推進速率的吻合程度要優(yōu)于后期。同時還可以看出,水滴尺度越小或冷卻面溫度越低,凝固速率相對越高。

        為考查所建模型對凝固過程溫度變化特性預測的有效性,圖9顯示了半徑為2 mm的水滴在過冷度為-10 ℃條件下界面推進過程的溫度變化特性紅外測試圖像及溫度曲線??梢钥闯?,界面推進階段,水滴由過冷態(tài)上升至平衡態(tài),水滴溫度維持在0 ℃附近并保持相對恒定,當凝固完成,水滴由平衡溫度逐步降低并趨于冷卻面溫度。

        圖7 水滴凝固過程相界面變化特性的實驗與計算對比
        Fig.7 Comparison of experimental and simulated phase interface evolution characteristics of droplet freezing

        圖8 不同水滴直徑凝固過程相界面變化特性比較
        Fig.8 Phase interface evolution characteristics of droplet freezing with different diameters

        圖9 水滴凝固過程溫度變化特性
        Fig.9 Temperature evolution characteristics of droplet freezing process

        圖10 水滴凝固過程溫度變化特性的實驗與計算對比
        Fig.10 Comparison of experimental and simulated temperature evolution characteristics of droplet freezing

        圖10則顯示了過冷度分別為-10 ℃和-14 ℃ 條件下,直徑為4 mm水滴凝固界面推進過程中水滴表面計算溫度與實驗紅外溫度變化特性的比較??梢钥闯觯捎跍y量誤差及計算模型忽略變形效應假設帶來的誤差,計算與實驗溫度變化歷程在凝固后期存在一定程度的偏差,但總體變化趨勢基本相似。且在凝固期間,與一般的液/固相變過程相同,計算與實驗水滴溫度均維持在平衡溫度即0 ℃附近,較好反映了水滴相變過程的恒溫特性。

        4 結(jié) 論

        1) 基于Enthalpy-Porosity模型,發(fā)展了過冷水滴結(jié)冰特性預測模型及數(shù)值計算方法,并自行搭建了水滴凝固實驗臺,開展了數(shù)值方法的實驗驗證。所發(fā)展方法計算結(jié)果與實驗結(jié)果吻合較好,當過冷度為0 ℃時,過冷模型退化為傳統(tǒng)Enthalpy-Porosity模型,表明所建方法的有效性。

        2) 相對于傳統(tǒng)方法,所發(fā)展的過冷水滴結(jié)冰特性預測模型及數(shù)值方法能有效表征非平衡凝固效應,可用于過冷水滴凝固特性的預測,從而將傳統(tǒng)基于平衡凝固的Enthalpy-Porosity模型拓展至非平衡凝固研究領域。

        3) 基于所建方法,開展了水滴凝固特性的影響因素分析,獲得了不同過冷條件下水滴凝固過程的溫度分布及相界面變化特征。研究表明,過冷態(tài)溫度越大、水滴尺度越小或冷卻面溫度越低,水滴的相變速率越高;在考慮非平衡凝固效應的條件下,過冷水滴凝固速率要高于不考慮非平衡凝固效應的工況。

        本研究將進一步考慮水滴結(jié)冰過程的體積膨脹效應,為結(jié)冰熱力學模型的改進以及結(jié)冰特性的精細化預測提供參考。

        [1] 王洪偉, 李先哲, 宋展. 通用飛機結(jié)冰適航驗證關鍵技術及工程應用[J]. 航空學報, 2016, 37(1): 335-350.

        WANG H W, LI X Z, SONG Z. Key airworthiness validation technologies for icing of general aviation aircraft and their engineering application[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(1): 335-350 (in Chinese).

        [2] BLAKE J, THOMPSON D, RAPS D, et al. Simulating the freezing of supercooled water droplets impacting a cooled substrate[J]. AIAA Journal, 2015, 53(7): 1725-1739.

        [3] FUMOTO K, KAWANAMI T. Study on freezing characteristics of supercooled water droplets impacting on solid surfaces[J]. Journal of Adhesion Science and Technology, 2012, 26(4-5): 463-472.

        [4] KING W D. Freezing rates of water droplets[J]. Journal of the Atmospheric Sciences, 1975, 32(2): 403-408.

        [5] TABAKOVA S, FEUILLEBOIS F, RADEV S. Freezing of a suspended supercooled droplet with a heat transfer mixed condition on its outer surface[C]//1st International Conference on Applications of Mathematics in Technical and Natural Sciences, 2009, 1186(1): 240-247.

        [6] FEUILLEBOIS F, LASEK A, CREISMEAS P, et al. Freezing of a subcooled liquid droplet[J]. Journal of Colloid and Interface Science, 1995, 169(1): 90-102.

        [7] BURTNETT E. Volume of fluid simulations for droplet impact on dry and wetted hydrophobic and superhydrophobic surfaces[D]. Mississippi: Mississippi State University, 2012.

        [8] JUNG S, DORRESTIJN M, RAPS D, et al. Are superhydrophobic surfaces best for icephobicity?[J]. Langmuir, 2011, 27(6): 3059-3066.

        [9] VOLLER V, PRAKASH C. A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems[J]. International Journal of Heat and Mass Transfer, 1987, 30(8): 1709-1719.

        [10] CRANK J. Free and moving boundary problems[M]. Oxford: Oxford University Press, 1987.

        [11] WORSTER M G. Solidification of fluids[M]. Cambridge: Cambridge University Press, 2000.

        [12] ELLEN N, JACCO M H, EDWIN W, et al. Aircraft icing in flight: Effects of impact of supercooled large droplets[C]//29th Congress of the Aeronautical Sciences, 2014.

        [13] TABAKOVA S, FEUILLEBOIS F. On the solidification of a supercooled liquid Droplet lying on a surface[J]. Journal of Colloid and Interface Science, 2004, 272(1): 225-234.

        [14] SZIMMAT J. Numerical simulation of solidification process in enclosures[J]. Heat Mass Transfer, 2002, 38(4): 279.

        (責任編輯: 李明敏)

        URL:www.cnki.net/kcms/detail/11.1929.V.20161205.1640.002.html

        Freezingcharacteristicsofsupercooledwaterdropletinconsiderationofnon-equilibriumeffect

        XIAOGuangming,DUYanxia*,WANGQiao,GUOLong,WANGMao

        StateKeyLaboratoryofAerodynamics,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China

        Non-equilibriumeffectisanimportantphenomenoninfreezingofsupercooledwaterdropletinaircrafticingprocess.Basedontheenthalpy-porositymodel,anumericalpredictionmethodforfreezingofsupercooledwaterdropletisdeveloped.Theexperimentalsystemfordropletfreezingisbuiltupandseveralexperimentsareperformedtovalidatethenumericalmethodproposed.Theresultsindicatethatthedevelopedmodelisvalidandcanbeusedtopredictthefreezingcharacteristicsofsupercooledwaterdroplet.Basedontheimprovedfreezingmodel,theinfluenceofthedegreeofsupercooledandcoolingconditionsonthecharacteristicsofsupercooleddropletareanalyzed.Whenthedegreeofsupercooledisdecreasedtozero,thedevelopedmodeldegeneratestothetraditionalmodel.Thegreaterorsmallerthedegreeofsupercooledorsmallerthedropletsis,therelativelyhigherfreezingrateis.Inconsiderationoftheeffectofnon-equilibriumconditions,thefreezingandmovingrateofinterfaceishigherthanthetraditionalmodel.Relatedresearchcanprovideimportantreferenceforimprovingicingthermodynamicmodelandrefiningthepredictionmethodforicingaccretion.

        supercooledwaterdroplet;icing;non-equilibriumeffect;solidification;improvedmodel

        2016-08-24;Revised2016-10-25;Accepted2016-11-23;Publishedonline2016-12-051640

        s:NationalNaturalScienceFoundationofChina(51308531,11672322);NationalBasicResearchProgramofChina(2015CB755800)

        .E-mailyanxiadu@163.com

        2016-08-24;退修日期2016-10-25;錄用日期2016-11-23; < class="emphasis_bold">網(wǎng)絡出版時間

        時間:2016-12-051640

        www.cnki.net/kcms/detail/11.1929.V.20161205.1640.002.html

        國家自然科學基金 (51308531,11672322); 國家“973”計劃 (2015CB755800)

        .E-mailyanxiadu@163.com

        肖光明, 杜雁霞, 王橋, 等. 考慮非平衡效應的過冷水滴凝固特性J. 航空學報,2017,38(2):520703.XIAOGM,DUYX,WANGQ,etal.Freezingcharacteristicsofsupercooledwaterdropletinconsiderationofnon-equilibriumeffectJ.ActaAeronauticaetAstronauticaSinica,2017,38(2):520703.

        http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

        10.7527/S1000-6893.2016.0309

        V211.3

        A

        1000-6893(2017)02-520703-07

        猜你喜歡
        結(jié)冰水滴界面
        水滴大變樣
        “水滴”船
        通體結(jié)冰的球
        國企黨委前置研究的“四個界面”
        當代陜西(2020年13期)2020-08-24 08:22:02
        冬天,玻璃窗上為什么會結(jié)冰花?
        基于FANUC PICTURE的虛擬軸坐標顯示界面開發(fā)方法研究
        魚缸結(jié)冰
        水滴瓶
        人機交互界面發(fā)展趨勢研究
        手機界面中圖形符號的發(fā)展趨向
        新聞傳播(2015年11期)2015-07-18 11:15:04
        国内自拍视频在线观看h| 色哟哟网站在线观看| 国产日产高清欧美一区| 色婷婷精品综合久久狠狠| 亚洲av手机在线观看| 国色天香社区视频在线| 性一交一乱一伦一色一情孩交| 亚洲国产精品久久久天堂不卡海量 | 欧美综合天天夜夜久久| 玩两个丰满老熟女| 人妻系列无码专区久久五月天| 亚洲女同一区二区三区| 一本久久综合亚洲鲁鲁五月天| 国内a∨免费播放| 91精品综合久久久久m3u8| 亚洲精品天堂日本亚洲精品 | 中文字幕无线码| 国内精品一区视频在线播放| 精品一区二区三区女同免费| 日本高清一级二级三级| 久久人妻内射无码一区三区| 亚洲A∨无码国产精品久久网| 亚洲精品一区二区三区日韩 | 亚洲av无码专区亚洲av伊甸园| 日韩欧美一区二区三区中文精品| 色婷婷亚洲十月十月色天| 国产流白浆视频在线观看| 99热这里有精品| 国产精品99久久久精品免费观看| 亚洲福利网站在线一区不卡| 国产精品免费一区二区三区四区| 色伦专区97中文字幕| 亚洲欧美日韩精品香蕉| 曰日本一级二级三级人人| 奇米影视第四色首页| 久久夜色撩人精品国产小说| av天堂手机一区在线| 麻豆国产精品久久人妻| 日韩精品一区二区亚洲av| 亚洲av激情久久精品人| 日韩不卡的av二三四区|