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

        ?

        剪切流場(chǎng)下液滴碰撞的流變特性

        2015-04-24 05:50:00王程遙張程賓陳永平
        關(guān)鍵詞:液膜液滴流場(chǎng)

        王程遙 張程賓 陳永平,3 張 林

        (1東南大學(xué)能源與環(huán)境學(xué)院,南京 210096)(2中國(guó)工程物理研究院激光聚變研究中心,綿陽(yáng) 621900)(3揚(yáng)州大學(xué)水利與能源動(dòng)力工程學(xué)院,揚(yáng)州 225127)

        剪切流場(chǎng)下液滴碰撞的流變特性

        王程遙1張程賓1陳永平1,3張 林2

        (1東南大學(xué)能源與環(huán)境學(xué)院,南京 210096)(2中國(guó)工程物理研究院激光聚變研究中心,綿陽(yáng) 621900)(3揚(yáng)州大學(xué)水利與能源動(dòng)力工程學(xué)院,揚(yáng)州 225127)

        基于VOF液/液相界面追蹤方法,建立了不可壓縮水/油單乳液液滴動(dòng)力學(xué)模型并進(jìn)行數(shù)值求解,模擬研究了剪切流場(chǎng)條件下2個(gè)相同體積的液滴在碰撞過(guò)程中的相互作用及變形行為.觀察了液滴碰撞過(guò)程中液滴的運(yùn)動(dòng)軌跡,并對(duì)相應(yīng)的內(nèi)在機(jī)理進(jìn)行了分析.在剪切流場(chǎng)作用下,兩液滴的碰撞過(guò)程分為接近、碰撞、分離3個(gè)階段.由于碰撞過(guò)程中液滴間的相互擠推作用,液滴分離后,液滴間的側(cè)向質(zhì)心間距Δy/a增大.此外,分析了液滴碰撞過(guò)程中毛細(xì)數(shù)對(duì)液滴間相互作用的影響.兩液滴在碰撞靠近過(guò)程中,在碰撞區(qū)中心處產(chǎn)生一個(gè)高壓區(qū),隨著毛細(xì)數(shù)Ca從0.2增加到0.4,界面擠壓變形越明顯,液滴變形系數(shù)D也從0.32增加到0.51.

        液滴;碰撞;剪切流場(chǎng);VOF方法

        液滴動(dòng)力學(xué)是多相流學(xué)科和乳液工程技術(shù)的重要研究領(lǐng)域,液滴碰撞過(guò)程的流變行為是乳化、共混等工業(yè)過(guò)程[1-3]中主動(dòng)調(diào)控的重要理論基礎(chǔ).因此,準(zhǔn)確掌握碰撞過(guò)程中液滴的流變行為、內(nèi)在機(jī)理及其影響規(guī)律不僅對(duì)液滴動(dòng)力學(xué)理論研究方向具有重要的學(xué)術(shù)價(jià)值,而且對(duì)于有效地控制液滴工藝流程具有重要的現(xiàn)實(shí)意義.

        目前,國(guó)際上已對(duì)單個(gè)液滴的流變行為開(kāi)展了較為深入的研究.20世紀(jì)30年代,Taylor等[4-5]就對(duì)單個(gè)液滴在流場(chǎng)中變形與破裂行為開(kāi)展了實(shí)驗(yàn)研究.隨后,Rallison等[6]、Youngren等[7]對(duì)單一液滴的形變行為進(jìn)行了數(shù)學(xué)建模和數(shù)值模擬.Stone[8]則對(duì)單液滴在流場(chǎng)中流變行為的實(shí)驗(yàn)與理論工作進(jìn)行了總結(jié),綜述了液滴界面動(dòng)力學(xué)方面的研究進(jìn)展,著重介紹了流體物性參數(shù)和不同的流場(chǎng)類(lèi)型等因素對(duì)液滴形變和液滴破碎的影響.最近,Salkin等[9]對(duì)液滴在流場(chǎng)中遇到障礙物的破碎行為進(jìn)行了實(shí)驗(yàn)研究.Ulloa等[10]實(shí)驗(yàn)研究了剪切流場(chǎng)的受限程度對(duì)液滴形變與破碎的影響.劉向東等[11]采用VOF方法對(duì)剪切流場(chǎng)中的單個(gè)雙乳液液滴進(jìn)行了液滴形變及破碎特性的研究,結(jié)果表明液滴的穩(wěn)定形變隨著毛細(xì)數(shù)的增大而加劇,并且在給定的雷諾數(shù)下存在一個(gè)臨界毛細(xì)數(shù),小于臨界毛細(xì)數(shù)時(shí),液滴只發(fā)生穩(wěn)定變形,反之液滴破碎.

        相比于單個(gè)液滴在外流場(chǎng)中的流變行為,2個(gè)或者更多液滴在外流場(chǎng)的碰撞過(guò)程中存在液滴間的相互作用,使得其流變機(jī)理變得更加復(fù)雜.目前,已有少部分學(xué)者在此領(lǐng)域開(kāi)展了相關(guān)的實(shí)驗(yàn)研究,針對(duì)液滴在剪切流場(chǎng)中碰撞過(guò)程的運(yùn)動(dòng)軌跡進(jìn)行了可視化觀測(cè)[12-13],但受到量測(cè)手段的限制,對(duì)液滴在剪切流場(chǎng)中相互碰撞過(guò)程的形變與碰撞的耦合機(jī)理以及壓力變化特性尚未充分揭示.為此,本文基于VOF(volume of fluid)液/液相界面追蹤方法,建立了不可壓縮水/油單乳液液滴動(dòng)力學(xué)模型并進(jìn)行數(shù)值求解,模擬研究了在剪切流場(chǎng)作用下2個(gè)相同體積的球形液滴在碰撞過(guò)程中的變形、相互作用及壓力演化.

        1 數(shù)學(xué)模型

        雙板平行剪切流場(chǎng)下2個(gè)半徑a=100 μm的油滴在基質(zhì)水中的初始位置如圖1所示,流動(dòng)方向?yàn)閤方向,速度梯度方向?yàn)閥方向,2個(gè)液滴關(guān)于中心點(diǎn)對(duì)稱(chēng),x方向與y方向的初始質(zhì)心距分別為Δx和Δy.2塊平板在2個(gè)平行平面上分別以相同速度U=0.02 m/s向相反的方向滑動(dòng),產(chǎn)生剪切率G=2U/H的恒定剪切流場(chǎng)(H為兩板間距),作用于液滴從而產(chǎn)生運(yùn)動(dòng)變形.剪切流場(chǎng)的左右兩側(cè)邊界采用周期性邊界條件.

        在剪切流場(chǎng)下,液滴的碰撞現(xiàn)象實(shí)質(zhì)上是一種液/液相界面的運(yùn)動(dòng)與變形所產(chǎn)生的宏觀結(jié)果.本文采用保持質(zhì)量守恒性較好的VOF方法[14]來(lái)描述相界面的位置及其運(yùn)動(dòng).VOF方法就是通過(guò)研究網(wǎng)格中每一相的體積比函數(shù)αi來(lái)確定介質(zhì)面、

        圖1 雙板平行剪切流場(chǎng)下兩個(gè)液滴初始位置示意圖

        自由面、間斷面以及其他內(nèi)部運(yùn)動(dòng)界面,能夠捕獲各種不規(guī)則形狀的自由表面,可以較好地模擬乳液液滴在流場(chǎng)中的變形、碰撞等復(fù)雜的變化.定義連續(xù)相m、離散相d兩相的體積比函數(shù)分別為αm和αd,則αm+αd=1.

        在整個(gè)流場(chǎng)計(jì)算區(qū)域內(nèi),滿足以下輸運(yùn)方程:

        (1)

        整個(gè)流場(chǎng)區(qū)域滿足以下連續(xù)性方程:

        (2)

        動(dòng)量方程

        (3)

        其中

        ρ=∑αiρi

        (4)

        μ=∑αiμi

        (5)

        式中,v為流體速度;ρ和μ分別為流體的密度和動(dòng)力黏度;張力源項(xiàng)Fvol采用CSF(continuum surface force)方法[15]計(jì)算,即

        (6)

        式中,m相和d相之間的界面曲率為

        (7)

        在簡(jiǎn)單剪切流中,液滴碰撞變形主要受到慣性力、黏性力及相界面張力的作用.為表征各相參數(shù)及各種作用力間的相對(duì)大小,本文引入以下參數(shù):Reynolds數(shù)(Re=ρmGa2/μm)表征慣性力與黏性力之間的相對(duì)大小;Capillary數(shù)(Ca=μmGa/σ)表征黏性力和表面張力之間的相對(duì)大小,其中,界面張力系數(shù)σ=10-5N/m;黏度比λ=μd/μm;密度比β=ρd/ρm.μd=μm=0.001 Pa·s,ρd=ρm=1 000 kg/m3.為表征液滴碰撞過(guò)程中的變形程度,引入Taylor變形系數(shù)[1-2]D=(L-B)/(L+B),其中L和B分別為發(fā)生變形后的橢圓形液滴的長(zhǎng)軸與短軸.本文將初始液滴半徑a作為特征長(zhǎng)度,剪切率的倒數(shù)G-1作為特征時(shí)間,則無(wú)量綱時(shí)間t*=t/G,液滴的初始位置參數(shù)為(Δx0/a, Δy0/a).

        對(duì)圖1兩平行平板之間的流場(chǎng)進(jìn)行矩形網(wǎng)格劃分.流體流動(dòng)控制方程(即方程(1)~(3))采用有限體積差分法進(jìn)行數(shù)值求解.在求解過(guò)程中,壓力場(chǎng)和速度場(chǎng)的耦合采用SIMPLE算法,動(dòng)量方程采用一階迎風(fēng)差分格式離散,基于分段式界面計(jì)算(PLIC)方法的Geo-Reconstruct格式對(duì)液/液界面進(jìn)行重構(gòu).控制微分方程經(jīng)離散成為線性代數(shù)方程組后,采用Gauss-Seidel方法進(jìn)行迭代求解.迭代計(jì)算中所采用的欠松弛因子分別為:壓力項(xiàng)0.2,密度項(xiàng)0.5,體積力源項(xiàng)0.5,動(dòng)量源項(xiàng)0.2.采用不同的網(wǎng)格尺寸進(jìn)行網(wǎng)格獨(dú)立性檢測(cè),以確保所得的數(shù)值解均為網(wǎng)格獨(dú)立的解.時(shí)間步長(zhǎng)取0.1 ms,當(dāng)每個(gè)迭代時(shí)間步的流體速度殘差小于10-3時(shí),認(rèn)為該時(shí)間步迭代計(jì)算收斂.

        本文對(duì)Guido等[12]開(kāi)展的液滴碰撞實(shí)驗(yàn)進(jìn)行了數(shù)值模擬對(duì)比,如圖2所示.由圖可見(jiàn),液滴碰撞演化過(guò)程的數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好,表明本文所建立的模型準(zhǔn)確可靠.

        (a) 液滴碰撞的實(shí)驗(yàn)圖片

        (b) 數(shù)值模擬的結(jié)果

        2 結(jié)果分析與討論

        2.1 碰撞過(guò)程中液滴間的相互作用

        為研究剪切流場(chǎng)下2個(gè)液滴碰撞過(guò)程的流變特性,本文對(duì)Re=0.2,Ca=0.2,H=2 mm條件下2個(gè)相同大小的液滴在剪切流場(chǎng)中的碰撞過(guò)程進(jìn)行數(shù)值模擬.液滴在x方向與y方向的初始質(zhì)心距分別為Δx0/a=4, Δy0/a=0.8.圖3給出了液滴在剪切流場(chǎng)中經(jīng)歷碰撞過(guò)程的運(yùn)動(dòng)軌跡.由于液滴間的相互作用,液滴間的y向距離先增大后減小,隨后穩(wěn)定于一個(gè)值,最終的y向距離比初始的距離大.由圖可知,隨著剪切流場(chǎng)開(kāi)始作用,液滴先變形為橢圓狀,隨后在流場(chǎng)拖曳力的作用下,液滴相互靠近并碰撞,直至液滴最終相互繞過(guò)對(duì)方,完成碰撞.在碰撞過(guò)程中,受液滴間強(qiáng)烈的相互作用影響,液滴產(chǎn)生了劇烈形變,以致在2個(gè)液滴碰撞作用處分別形成了一個(gè)平滑界面.隨后,兩液滴以流場(chǎng)中心為對(duì)稱(chēng)軸旋轉(zhuǎn)并相互繞過(guò)對(duì)方后分離,當(dāng)液滴完全分開(kāi)后,液滴形狀恢復(fù)為碰撞前的橢圓狀.

        (a) 液滴碰撞過(guò)程時(shí)間序列圖

        (b) 兩液滴相對(duì)位置軌跡

        為了更好地理解液滴碰撞過(guò)程中的流變行為,圖4給出了與圖3相對(duì)應(yīng)時(shí)間點(diǎn)的碰撞過(guò)程壓力分布時(shí)間序列圖,圖中的壓力是以101 325 Pa為基準(zhǔn)的相對(duì)壓力值.如圖所示,隨著施加外剪切流場(chǎng),液滴開(kāi)始變形,在表面張力的作用下,在變形液滴長(zhǎng)軸方向的端部區(qū)域出現(xiàn)高壓區(qū),從而在長(zhǎng)軸端部的界面內(nèi)外形成較高的壓力梯度.根據(jù)Young-Laplace方程Δp=2σ/R,界面兩端的壓力梯度與界面曲率半徑R呈反比關(guān)系,因此,在長(zhǎng)軸端部高壓力梯度區(qū)所對(duì)應(yīng)的液滴相界面的界面曲率半徑較小(見(jiàn)圖4(b)).隨后液滴相互靠近,液滴之間液膜的壓力開(kāi)始變大,促使液膜排液,液膜厚度逐漸減小,當(dāng)液膜最薄時(shí),液膜間的壓力達(dá)到最高值(見(jiàn)圖4(c)).此時(shí),膜壓的增大使得液膜產(chǎn)生擠壓作用,使液滴碰撞面附近產(chǎn)生大的變形,液滴的相對(duì)接觸面變平甚至凹陷.當(dāng)液滴開(kāi)始分開(kāi),液滴間的液膜厚度增加,同時(shí)伴隨著液滴的回縮,流場(chǎng)中心的壓力快速下降,變?yōu)樨?fù)壓(見(jiàn)圖4(d)~(e)).在壓力下降至極低值時(shí),低壓區(qū)產(chǎn)生抽吸效應(yīng),在y方向上產(chǎn)生了一個(gè)拖曳力,減小了該階段中液滴y向的質(zhì)心距(見(jiàn)圖4(e)).同時(shí),液滴的形狀出現(xiàn)短暫不對(duì)稱(chēng)的現(xiàn)象,液滴相對(duì)的端部出現(xiàn)高曲率的尖端.在Guido等[12]的實(shí)驗(yàn)中也觀察到了類(lèi)似的抽吸現(xiàn)象.最終,液滴完全分開(kāi),如圖4(f)所示.

        圖4 液滴碰撞過(guò)程壓力分布圖

        圖5給出了碰撞過(guò)程中液滴變形系數(shù)的變化趨勢(shì).隨著液滴相互靠近,變形量開(kāi)始變大;當(dāng)液滴相互對(duì)齊擠壓時(shí),液滴的形變量最大;隨之液滴相互越過(guò)對(duì)方,當(dāng)液滴開(kāi)始分開(kāi)時(shí),出現(xiàn)負(fù)壓抽吸效應(yīng),液滴被短暫拉長(zhǎng),形變量達(dá)到第2個(gè)最大值;當(dāng)液滴完全分開(kāi)后,液滴形狀恢復(fù),隨后基本保持不變.

        (a) 液滴碰撞過(guò)程時(shí)間序列圖

        (b) 液滴形變趨勢(shì)

        2.2 Ca數(shù)的影響

        為研究不同的Ca數(shù)對(duì)液滴碰撞的影響,本文在H=1.0mm,Ca=0.2,0.3,0.4的條件下對(duì)液滴碰撞過(guò)程進(jìn)行數(shù)值模擬研究,結(jié)果如圖6所示.由圖可知,Ca數(shù)的影響主要體現(xiàn)在液滴的形變上,Ca數(shù)表征了促使液滴變形的黏性剪切力與抵抗液滴變形的表面張力的比值.Ca數(shù)越大,黏性剪切力的相對(duì)值就增大,從而導(dǎo)致液滴的變形程度增大.同時(shí),隨著Ca數(shù)的增大,剪切力增大,碰撞過(guò)程中,剪切流場(chǎng)對(duì)2個(gè)液滴的外圍擠推作用增強(qiáng),強(qiáng)化了對(duì)碰撞面的擠壓作用,因此,兩液滴碰撞面處的界面內(nèi)凹現(xiàn)象變得愈加明顯.

        圖6 不同Ca數(shù)下液滴的形變

        3 結(jié)論

        1) 在剪切流場(chǎng)作用下,兩液滴碰撞過(guò)程分為接近、碰撞、分離3個(gè)階段.碰撞過(guò)程中,液滴先變形為橢圓狀,隨后相互靠近、碰撞,然后液滴相互繞過(guò)對(duì)方,最終,兩液滴分離并恢復(fù)至碰撞前的橢圓狀.由于碰撞過(guò)程中液滴間的相互擠推作用,液滴分離后,垂直于液滴碰撞方向的質(zhì)心間距增大.

        2) 碰撞過(guò)程中,隨著液滴相互靠近,液滴的變形量開(kāi)始變大;當(dāng)液滴碰撞并相互擠壓時(shí),液滴的形變量最大;隨之液滴相互繞過(guò)對(duì)方并分開(kāi),此時(shí),由于低壓區(qū)的抽吸拖曳效應(yīng),兩液滴相對(duì)的端部出現(xiàn)尖端,液滴形變量達(dá)到第2個(gè)最大值;當(dāng)液滴完全分開(kāi)后,液滴形狀恢復(fù).

        3) 碰撞過(guò)程中,液滴碰撞區(qū)域處會(huì)形成平滑界面.兩液滴碰撞靠近過(guò)程中,在碰撞區(qū)中心處產(chǎn)生一個(gè)高壓區(qū),使得兩液滴間的液膜對(duì)碰撞的界面產(chǎn)生擠壓作用,碰撞界面變平,并且隨著毛細(xì)數(shù)Ca從0.2增加到0.4,界面擠壓變形越明顯,液滴變形系數(shù)也從0.32增加到0.51.

        References)

        [1]Jafari S M, Assadpoor E, He Y H, et al. Re-coalescence of emulsion droplets during high-energy emulsification [J].FoodHydrocolloids, 2008, 22(7): 1191-1202.

        [2]Lobo L, Svereika A. Coalescence during emulsification: 2. role of small molecule surfactants [J].JournalofColloidandInterfaceScience, 2003, 261(2): 498-507.

        [3]Ramic A J, Hudson S D, Jamieson A M, et al. Temporary droplet-size hysteresis in immiscible polymer blends [J].Polymer, 2000, 41(16): 6263-6270.

        [4]Taylor G I. The viscosity of a fluid containing small drops of another fluid [J].ProceedingsoftheRoyalSocietyA:MathematicalPhysicalandEngineeringSciences, 1932, 138(834): 41-48.

        [5]Taylor G I. The formation of emulsions in definable fields of flow [J].ProceedingsoftheRoyalSocietyA:MathematicalPhysicalandEngineeringSciences, 1934, 146(858):501-523.

        [6]Rallison J M, Acrivos A. A numerical study of deformation and burst of a viscous drop in an extensional flow [J].JournalofFluidMechanics, 1978, 89(1):191-200.

        [7]Youngren G K, Acrivos A. On the shape of a gas bubble in a viscous extensional flow [J].JournalofFluidMechanics, 1976, 76(3):433-442.

        [8]Stone H A. Dynamics of drop deformation and breakup in viscous fluids [J].AnnualReviewofFluidMechanics, 1994, 26: 65-102.

        [9]Salkin L, Courbin L, Panizza P. Microfluidic breakups of confined droplets against a linear obstacle: the importance of the viscosity contrast [J].PhysicalReviewE, 2012, 86: 036317-1-036317-6.

        [10]Ulloa C, Ahumada A, Cordero M L. Effect of confinement on the deformation of microfluidic drops [J].PhysicalReviewE, 2014, 89: 033004-1-033004-6.

        [11]劉向東, 陳永平, 張程賓, 等. 雙乳液液滴形變及破碎特性[J]. 工程熱物理學(xué)報(bào), 2012, 33(3): 457-459. Liu Xiangdong, Chen Yongping, Zhang Chengbin, et al. Deformation and breakup characteristics of double emulsion drop [J].JournalofEngineeringThermophysics, 2012, 33(3): 457-459.(in Chinese)

        [12]Guido S, Simeone M. Binary collision of drops in simple shear flow by computer-assisted video optical microscopy [J].JournalofFluidMechanics, 1998, 357:1-20.

        [13]Tretheway D C, Muraoka M, Leal L G. Experimental trajectories of two drops in planar extensional flow [J].PhysicsofFluids, 1999, 11(5): 971-981.

        [14]Hirt C W, Nichols B D. Volume of fluid (VOF) method for the dynamics of free boundaries [J].JournalofComputationalPhysics, 1981, 39(1): 201-225.

        [15]Brackbill J U, Kothe D B, Zemach C. A continuum method for modeling surface tension [J].JournalofComputationalPhysics, 1992, 100(2): 335-354.

        Rheological behavior of interactive drops in shear flow

        Wang Chengyao1Zhang Chengbin1Chen Yongping1,3Zhang Lin2

        (1School of Energy and Environment, Southeast University, Nanjing 210096, China)(2Research Center of Laser Fusion, China Academy of Engineering Physics, Mianyang 621900, China)(3School of Hydraulic, Energy and Power Engineering, Yangzhou University, Yangzhou 225127, China)

        Based on VOF (volume of fluid) liquid-liquid interface tracking method, a model for incompressible water-oil emulsion drop hydrodynamics is developed and numerically solved to study the hydrodynamic interaction and deformation behavior between two equal-sized drops during the collision in a shear flow. The motion trajectory during the collision is observed and the corresponding underlying hydrodynamics are analyzed. Under the shear flow, drop behavior during the collision can be divided into three stages: approach, collision and separation. Due to the hydrodynamic interaction between the drops, the final lateral separation of drop mass centers Δy/aincreases after the drop separation. In addition, the effects of the capillary number on the hydrodynamic interaction between two interactive drops are also analyzed. During the collision, a high pressure region is observed at the center of the flow domain. As the capillary number increases from 0.2 to 0.4, the interface deformation is more obvious and the deformation coefficientDof the drop increases from 0.32 to 0.51.

        drop; collision; shear flow; volume of fluid method

        10.3969/j.issn.1001-0505.2015.02.020

        2014-12-07. 作者簡(jiǎn)介: 王程遙(1988—),女,博士生;陳永平(聯(lián)系人),男,博士,教授,博士生導(dǎo)師, ypchen@seu.edu.cn.

        江蘇省自然科學(xué)基金資助項(xiàng)目(BK20130009,BK20130621)、中國(guó)工程物理研究院科學(xué)技術(shù)發(fā)展基金資助項(xiàng)目(2012A0302015).

        王程遙,張程賓,陳永平,等.剪切流場(chǎng)下液滴碰撞的流變特性[J].東南大學(xué)學(xué)報(bào):自然科學(xué)版,2015,45(2):309-313.

        10.3969/j.issn.1001-0505.2015.02.020

        TK124

        A

        1001-0505(2015)02-0309-05

        猜你喜歡
        液膜液滴流場(chǎng)
        考慮軸彎曲的水潤(rùn)滑軸承液膜建模方法
        高空高速氣流下平板液膜流動(dòng)與破裂規(guī)律
        大型空冷汽輪發(fā)電機(jī)轉(zhuǎn)子三維流場(chǎng)計(jì)算
        液膜破裂對(duì)PCCS降膜的影響*
        液滴間相互碰撞融合與破碎的實(shí)驗(yàn)研究
        噴淋液滴在空氣環(huán)境下的運(yùn)動(dòng)特性
        轉(zhuǎn)杯紡排雜區(qū)流場(chǎng)與排雜性能
        基于HYCOM的斯里蘭卡南部海域溫、鹽、流場(chǎng)統(tǒng)計(jì)分析
        基于瞬態(tài)流場(chǎng)計(jì)算的滑動(dòng)軸承靜平衡位置求解
        豎直窄矩形通道內(nèi)彈狀流中液膜特性研究
        日韩精品人妻少妇一区二区| 国产熟人av一二三区| 国产成人精品一区二区三区av| 久久久噜噜噜噜久久熟女m| 国产青春草在线观看视频| 手机在线观看亚洲av| 久久亚洲av午夜福利精品西区| 亚洲一区二区精品在线| 人妻少妇偷人精品免费看| 老女老肥熟女一区二区| 成年女人a毛片免费视频| 国产无遮挡又爽又刺激的视频老师 | av人妻在线一区二区三区| 隔壁人妻欲求不满中文字幕| 中文字幕午夜精品一区二区三区 | 精品亚洲少妇一区二区三区| 69天堂国产在线精品观看| 久久久久成人精品免费播放| 亚洲国产成a人v在线观看| 加勒比特在线视频播放| av中国av一区二区三区av| 亚洲国产中文字幕在线视频综合| 夜夜骚久久激情亚洲精品| 欧美成人看片一区二区三区尤物 | 久久精品国产精品亚洲艾| 亚洲熟女少妇一区二区三区青久久| 亚洲国产精品情侣视频| 亚洲精品粉嫩美女一区| 一个人看的www片免费高清视频| 国产国语亲子伦亲子| 亚洲愉拍99热成人精品热久久| 越南女子杂交内射bbwbbw| 66lu国产在线观看| 精品综合久久久久久8888| japanese色国产在线看视频| 成年女人18毛片毛片免费| 成年人视频在线播放麻豆| 久久狼人国产综合精品| 日韩精品 在线 国产 丝袜| 特黄大片又粗又大又暴| 朝鲜女子内射杂交bbw|