• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    球形氣泡界面變化對(duì)尾渦性質(zhì)和尺寸的影響

    2017-10-13 06:07:25費(fèi)洋龐明軍
    化工學(xué)報(bào) 2017年9期
    關(guān)鍵詞:尾渦尾流氣泡

    費(fèi)洋,龐明軍

    ?

    球形氣泡界面變化對(duì)尾渦性質(zhì)和尺寸的影響

    費(fèi)洋,龐明軍

    (常州大學(xué)機(jī)械工程學(xué)院,江蘇省綠色過(guò)程裝備重點(diǎn)實(shí)驗(yàn)室,江蘇常州213164)

    利用計(jì)算流體力學(xué)法研究了中等Reynolds數(shù)下(25≤≤500)氣泡界面污染程度對(duì)其尾流的影響。借鑒圓球繞流和停滯帽模型,提出了一種模擬中等Reynolds數(shù)下受污染球形氣泡尾流的三維模型,氣泡界面污染程度取決于帽角()的大小,帽角越大表示氣泡表面污染程度越小。研究發(fā)現(xiàn):=25~200時(shí),污染程度的減小會(huì)減小尾渦長(zhǎng)度()、分離角()以及渦中心位置()和()的數(shù)值,但不會(huì)改變其與Reynolds數(shù)表征的關(guān)系;污染程度的減小會(huì)使=250~500時(shí)尾渦的三維特性減弱,使=350時(shí)有序脫落的尾渦的強(qiáng)度減小并最終使其不發(fā)生脫落,使=500時(shí)無(wú)規(guī)律脫落的尾渦的無(wú)序性減弱并最終使其不發(fā)生脫落。

    氣泡;界面;尾流;層流;流體動(dòng)力學(xué)

    引 言

    許多工業(yè)領(lǐng)域比如核能發(fā)電、化工、食品加工中都存在泡狀流動(dòng)。在這些流動(dòng)中,氣泡或多或少會(huì)受到污染,氣泡受污染后它的界面會(huì)發(fā)生變化,界面變化會(huì)使氣泡的尾流發(fā)生變化,進(jìn)而使氣泡的阻力系數(shù)和其他水動(dòng)力學(xué)特性發(fā)生變化。因此,深入研究氣泡界面變化對(duì)尾流的影響具有重要意義。為了研究氣泡界面污染程度對(duì)其水動(dòng)力學(xué)特性的影響,目前國(guó)內(nèi)外研究者均以在液相中加入不同濃度的表面活性劑作為污染源來(lái)開(kāi)展研究。表面活性劑分子作為污染源會(huì)吸附在氣泡表面上。對(duì)于一個(gè)上浮氣泡,表面活性劑分子在對(duì)流的影響下會(huì)從氣泡上部表面不斷地?cái)U(kuò)散到氣泡的尾部形成聚集現(xiàn)象[1-6]。

    目前有關(guān)氣泡尾流的研究主要集中于純凈氣泡,Bhaga等[7]實(shí)驗(yàn)研究了不同E?tv?s數(shù)()、Morton數(shù)()和Reynolds數(shù)()下氣泡的尾渦長(zhǎng)度、寬度和渦中心位置隨Reynolds數(shù)的變化關(guān)系;倪明玖[8]使用VOF方法研究了氣泡尾部尾渦的形成機(jī)理,發(fā)現(xiàn)尾渦的形成主要源于氣泡頂部與底部存在的壓差,同時(shí)也受氣泡尾部環(huán)流的影響;Tripathi等[9]模擬發(fā)現(xiàn)橢球形氣泡沿直線上升時(shí)的尾渦是對(duì)稱的,沿曲線上升時(shí)的尾渦是不對(duì)稱的;Gumulya等[10]模擬研究了橢球帽形氣泡和裙?fàn)顨馀莸奈矞u長(zhǎng)度隨雷諾數(shù)的變化規(guī)律。有關(guān)受污染氣泡的尾流研究,Veldhuis等[11]實(shí)驗(yàn)發(fā)現(xiàn)氣泡在純凈水和自來(lái)水中上浮時(shí)的尾流結(jié)構(gòu)明顯不同;Saito等[12]和Huang等[13-14]實(shí)驗(yàn)發(fā)現(xiàn)因表面活性劑在氣泡表面分布不均而產(chǎn)生的Marangoni效應(yīng)會(huì)對(duì)氣泡的尾流產(chǎn)生影響。由上述可知,目前對(duì)純凈氣泡尾流的研究相對(duì)較深入,但對(duì)受污染氣泡尾流的研究?jī)H局限于氣泡污染前后其尾流的變化情況,沒(méi)有涉及氣泡表面污染程度不同時(shí),其尾流的變化情況。因此,有必要詳細(xì)研究氣泡界面污染程度對(duì)其尾渦性質(zhì)、尺寸和脫落情況的影響。

    針對(duì)受污染氣泡的數(shù)值研究,主要有模型假設(shè)和真實(shí)模擬兩種方法。模型假設(shè)是基于對(duì)氣泡表面給予一定的假設(shè)進(jìn)行研究的,主要有均勻模型[15]、停滯帽模型[16-21]和摩擦系數(shù)模型[22]。均勻模型假設(shè)氣泡表面均勻污染,但這一假設(shè)不符合運(yùn)動(dòng)氣泡上表面活性劑的實(shí)際分布,故目前的模擬大多基于圖1所示的停滯帽模型,但停滯帽模型又難以對(duì)界面污染程度進(jìn)行定量控制。摩擦系數(shù)模型用摩擦系數(shù)來(lái)表示氣泡的污染,總體上與停滯帽模型類似。雖然近些年相關(guān)學(xué)者不使用任何模型假設(shè)實(shí)現(xiàn)了真實(shí)情況下受污染界面的模擬,但由于交界面處的輸運(yùn)方程需在伸縮與變形的交界面上求解,且還伴有網(wǎng)格的重構(gòu),故大多均停留在二維模擬階段[23-25]。然而Fleckenstein等[3]通過(guò)對(duì)比二維和三維的模擬結(jié)果發(fā)現(xiàn),兩種維度下氣泡的尾流情況大不相同。因此,想要更真實(shí)獲得受污染氣泡的尾流情況,就必須進(jìn)行三維模擬。

    本文使用計(jì)算流體力學(xué)法研究了中等Reynolds數(shù)下球形氣泡界面污染程度對(duì)其尾渦性質(zhì)、尺寸和脫落情況的影響。由于氣泡在中等Reynolds數(shù)下運(yùn)動(dòng)時(shí)其Marangoni數(shù)通常小于1(=¥/b,是氣體常數(shù),是吸附溫度,¥是氣泡界面上的極限濃度,是液相的動(dòng)力黏度,b是氣泡的運(yùn)動(dòng)速度),即氣泡表面兩種界面條件接近于階躍過(guò)渡[19-20],為方便研究,本文借鑒圓球繞流和停滯帽模型,提出了一種模擬中等Reynolds數(shù)下受污染球形氣泡尾流的三維簡(jiǎn)化模型,通過(guò)改變帽角大小來(lái)表示氣泡表面的污染程度,即帽角越大表示受污染程度越小。這一模型有效地避免了現(xiàn)有模型難以實(shí)現(xiàn)三維模擬以及單純停滯帽模型難以定量控制界面污染程度的缺點(diǎn)。

    1 數(shù)值計(jì)算

    1.1 幾何建模及邊界條件

    鑒于目前所研究的問(wèn)題,參考文獻(xiàn)[18-19,26]的建模方法,采用如圖2所示直徑為的圓球域作為計(jì)算區(qū)域,將直徑為的氣泡固定在圓球域中心,而流體相對(duì)氣泡發(fā)生流動(dòng),即,將-左半球設(shè)置為速度進(jìn)口邊界條件,進(jìn)口液相速度為,將-右半球設(shè)置為出流邊界條件(除了靜壓,其他所有參數(shù)均為零)。對(duì)于圓球域與氣泡的直徑比/,本文綜合考慮表1所示各學(xué)者[18-19,26]的取值,為同時(shí)保證模擬的準(zhǔn)確性和節(jié)省計(jì)算時(shí)間,取/=64。坐標(biāo)原點(diǎn)設(shè)在氣泡中心上,軸正方向?yàn)榱黧w流動(dòng)方向。為便于生成六面體網(wǎng)格,將氣泡表面和圓球域表面六等分。為考慮不同污染程度時(shí)氣泡的尾流情況,參考停滯帽模型,將氣泡分為圖3和圖4所示的迎流面和背流面兩部分,并將背流面分成圖4所示的21小塊,然后通過(guò)改變邊界條件來(lái)控制氣泡的污染程度。若帽角=0°(表示氣泡完全污染),將氣泡表面上所有的面都設(shè)為無(wú)滑移壁面邊界條件,即氣泡表面存在切應(yīng)力、不存在滑移速度,此時(shí)氣泡就類似于一個(gè)固體顆粒;若帽角=180°(表示氣泡未受污染),將氣泡表面上所有的面都設(shè)為無(wú)穿透壁面邊界條件,即氣泡表面存在滑移速度,不存在切應(yīng)力;若帽角在0°~180°之間,表示氣泡部分污染,氣泡表面上兩種邊界條件共存。為了便于建模,本文基于二分法,取部分污染氣泡的帽角90°,135°,157.5°和168.75°。若帽角90°,將氣泡的背流面,即圖4中的面6-26設(shè)為無(wú)滑移壁面邊界條件;若帽角135°,將圖4中氣泡背流面中的10-26設(shè)為無(wú)滑移壁面邊界條件;若帽角157.5°,將圖4中氣泡背流面中的18-26設(shè)為無(wú)滑移壁面邊界條件;若帽角168.75°,將圖4中氣泡背流面中的26設(shè)為無(wú)滑移壁面邊界條件。

    1.2 網(wǎng)格劃分

    對(duì)于網(wǎng)格劃分,本文采用六面體網(wǎng)格。由于球形氣泡和圓球域?yàn)檩S對(duì)稱結(jié)構(gòu),因此流域內(nèi)的總網(wǎng)格單元數(shù)由圓球域投影圓周上的網(wǎng)格份數(shù)和沿流域半徑上的網(wǎng)格份數(shù)決定。本文設(shè)氣泡投影圓周上的網(wǎng)格份數(shù)、流域投影圓周上的網(wǎng)格份數(shù)與流域半徑上的網(wǎng)格份數(shù)相等。由于氣泡表面的網(wǎng)格數(shù)也會(huì)對(duì)模擬結(jié)果產(chǎn)生一定影響,本文通過(guò)借鑒表1所示各學(xué)者[18-19,26]在流域投影圓周上的網(wǎng)格份數(shù),為同時(shí)保證模擬的準(zhǔn)確性和節(jié)省計(jì)算時(shí)間,取流域投影圓周上的網(wǎng)格份數(shù)為128個(gè),即氣泡表面的網(wǎng)格數(shù)為6144個(gè)、流域內(nèi)網(wǎng)格單元數(shù)約為78.64萬(wàn)個(gè)。有關(guān)氣泡投影圓周上網(wǎng)格份數(shù)的無(wú)關(guān)性檢驗(yàn)見(jiàn)2.1節(jié)。

    表1 域的尺寸與氣泡投影圓周上的網(wǎng)格數(shù)

    此外,為了更好地捕捉氣泡尾流的變化,本文進(jìn)行非均勻網(wǎng)格劃分,使靠近氣泡表面的網(wǎng)格較為密集,遠(yuǎn)離氣泡表面的網(wǎng)格較為稀疏。加密后氣泡表面第一層網(wǎng)格厚度與為滿足梯度要求而所需的邊界層厚度有關(guān),表2中羅列了Tomboulides等[27]給出的不同Reynolds數(shù)下所需邊界層厚度的需求值re以及本文加密后氣泡表面第1層網(wǎng)格的厚度use,可以看出實(shí)際值小于需求值,即這一比例下的第1層網(wǎng)格厚度滿足梯度要求。球形氣泡表面上的網(wǎng)格和流域中某一截面的網(wǎng)格分別如圖5和圖6所示。

    表2 不同Reynolds數(shù)下所需邊界層厚度的理論值和第1層網(wǎng)格厚度的實(shí)際值

    1.3 控制方程與數(shù)值方法

    從Kim等[28]給出的圓球繞流實(shí)驗(yàn)結(jié)果可以看出,當(dāng)Reynolds數(shù)大于800時(shí),尾流才會(huì)表現(xiàn)出湍流性質(zhì),并出現(xiàn)小尺度的渦結(jié)構(gòu)。因此,為節(jié)省計(jì)算時(shí)間,采用較為簡(jiǎn)單的層流模型進(jìn)行計(jì)算。此時(shí),需求解的連續(xù)性方程和動(dòng)量方程如式(1)和式(2)所示[29]

    ??=0 (1)

    式中,是液相速度,m·s-1;是靜壓,Pa;是液相動(dòng)力黏度,Pa·s。

    下文中的阻力系數(shù)由阻力求得,阻力由黏性力和壓力梯度力組成,如式(3)~式(6)所示[30]

    (4)

    (5)

    式中,D為阻力,N;DV為黏性力,N;DP為壓力梯度力,N;D為阻力系數(shù);為氣泡直徑,m;b為氣泡運(yùn)動(dòng)速度(這里等于進(jìn)口液相速度),m·s-1;為液相密度,kg·m-3;為氣泡表面的靜壓,Pa;x為氣泡表面切應(yīng)力的流向分量,Pa;為氣泡表面積,m2;x為氣泡表面法向單位矢量的流向分量。

    本文使用非穩(wěn)態(tài)時(shí)間格式迭代求解,求解時(shí)壓力-速度耦合采用SIMPLEC算法以加快收斂速度,單元中心的變量梯度采用最小二乘法離散,對(duì)流項(xiàng)采用具有三階精度的QUICK離散格式以減小偽擴(kuò)散,壓力項(xiàng)采用二階精度離散格式,時(shí)間項(xiàng)采用二階隱式離散格式。參考文獻(xiàn)[27],時(shí)間步長(zhǎng)的設(shè)置如表3所示。

    表3 不同Reynolds數(shù)下的時(shí)間步長(zhǎng)

    對(duì)于后文中的模擬結(jié)果,長(zhǎng)度尺度用氣泡直徑量綱1化,速度尺度用進(jìn)口液相速度量綱1化,時(shí)間尺度用/量綱1化。

    1.4 工況設(shè)置

    從圓球繞流的研究結(jié)果可知,在≤500的中等Reynolds數(shù)下,存在3種臨界Reynolds數(shù):1=212[27],2=277.5[31]和3=420[32]。當(dāng)≤210時(shí),為穩(wěn)態(tài)軸對(duì)稱尾流[33];當(dāng)=220~270時(shí),尾流由兩個(gè)平面對(duì)稱的發(fā)夾渦組成[33];當(dāng)=280~420時(shí),渦開(kāi)始脫落,但渦脫落的方向相同[31-32];當(dāng)>420時(shí),尾渦開(kāi)始無(wú)規(guī)律脫落[32]。因此本文將依據(jù)文獻(xiàn)劃分的3種區(qū)間,分別對(duì)受污染氣泡的尾流進(jìn)行分析和討論。此外,為對(duì)≤200時(shí)尾流隨Reynolds數(shù)的變化規(guī)律進(jìn)行深入分析,故在此區(qū)間內(nèi)取多個(gè)Reynolds數(shù)進(jìn)行模擬,具體工況設(shè)置如表4所示。

    表4 工況設(shè)置

    2 結(jié)果分析

    2.1 網(wǎng)格無(wú)關(guān)性檢驗(yàn)

    氣泡圓周上網(wǎng)格的份數(shù)會(huì)對(duì)計(jì)算的準(zhǔn)確性造成影響,本文分別取氣泡投影圓周上網(wǎng)格的份數(shù)為64、96、128、160和192份,并在最大Reynolds數(shù)=500下,對(duì)干凈氣泡和完全污染氣泡的阻力系數(shù)進(jìn)行監(jiān)測(cè),計(jì)算收斂時(shí)的監(jiān)測(cè)結(jié)果如表5所示??梢钥闯?,氣泡圓周上網(wǎng)格的份數(shù)為64份時(shí)的模擬結(jié)果與其他模擬結(jié)果相比有明顯偏差,本文為同時(shí)保證模擬結(jié)果的準(zhǔn)確性和節(jié)省時(shí)間花費(fèi),取氣泡圓周上網(wǎng)格的份數(shù)為128份。

    表5 氣泡圓周上不同的網(wǎng)格份數(shù)對(duì)干凈氣泡和完全污染氣泡的阻力系數(shù)的影響 (Re=500)

    2.2 結(jié)果準(zhǔn)確性檢驗(yàn)

    為確認(rèn)本文提出的簡(jiǎn)化模型和參數(shù)設(shè)置的準(zhǔn)確性,下文將模擬所得完全污染氣泡和干凈氣泡在不同下的阻力系數(shù)與Mei等[34-35]給出的公式[式(7)和式(8)]所對(duì)應(yīng)的阻力系數(shù)值進(jìn)行比較,如圖7所示。從中可以看出,完全污染氣泡的阻力系數(shù)與Mei等[34]給出的公式吻合得很好,干凈氣泡的阻力系數(shù)與Mei等[35]給出的公式有少量的偏差,這是由于真實(shí)情況下氣泡內(nèi)存在微弱的流動(dòng),流動(dòng)強(qiáng)度隨污染程度的減小而增加[36],而本簡(jiǎn)化模型忽略了氣泡內(nèi)部的流場(chǎng),但與真實(shí)情況相比偏差相對(duì)較小。

    界面完全污染球形氣泡阻力系數(shù)的計(jì)算表達(dá)式[34]

    界面完全干凈球形氣泡阻力系數(shù)的計(jì)算表達(dá)式[35]

    (8)

    有關(guān)結(jié)果準(zhǔn)確性的進(jìn)一步檢驗(yàn)見(jiàn)下文圖11。

    2.3 氣泡界面變化對(duì)其尾流的影響分析

    2.3.1=25~200時(shí)尾渦性質(zhì)隨帽角的變化 由于=25~200時(shí)氣泡尾流的表現(xiàn)形式相同,故下文僅對(duì)圖8所示=150時(shí),不同帽角下的尾流進(jìn)行分析。從圖中可以看出,隨著帽角的增加,尾渦的面積不斷減小,當(dāng)帽角≥157.5°時(shí),尾渦將完全消失。由此可見(jiàn),圖7中同一下干凈氣泡的阻力系數(shù)小于受污染氣泡的阻力系數(shù)正是由尾渦面積的減小和尾渦的消失造成的。此外,Takagi等[37]認(rèn)為氣泡受污染后減速的原因源于阻力系數(shù)的增加,由此可以推斷,尾渦面積的增大能降低氣泡在靜置液中的上升速度。這與通常所認(rèn)為的氣泡上浮減速源于交界面流動(dòng)性的降低并不矛盾,因?yàn)榻唤缑媪鲃?dòng)性的降低在本文中就是指帽角的減小,帽角的減小會(huì)使尾渦面積增大。這一點(diǎn)在能量守恒上是成立的,因?yàn)闅馀萆仙俣冉档秃蠊?jié)約下來(lái)的動(dòng)能被尾渦所消耗。對(duì)氣泡尾渦尺寸隨Reynolds數(shù)變化的具體分析見(jiàn)2.3.2節(jié)。

    此外,從圖8中還可以看出,對(duì)于有尾渦形成的帽角區(qū)間,帽角0°的完全污染氣泡,其尾流的表現(xiàn)形式與Tomboulides等[27]模擬的剛性顆粒的尾流表現(xiàn)形式相同,呈二維軸對(duì)稱形式;帽角=90°的半污染氣泡,其尾流也呈軸對(duì)稱形式;而帽角=135°的部分污染氣泡,尾流的表現(xiàn)形式與完全污染和半污染時(shí)的不同,有微弱的三維特性產(chǎn)生,為平面對(duì)稱形式。出現(xiàn)這一現(xiàn)象的原因與氣泡表面流向靜壓在帽角處的跳躍有關(guān),圖9給出了帽角=0°,90°和135°下氣泡表面壓力系數(shù)p沿流向分布,可以看出,帽角=135°時(shí)氣泡表面流向靜壓在帽角處有一個(gè)跳躍現(xiàn)象,它的存在使易受干擾的分離流線發(fā)生了變化,而帽角=0°和90°時(shí)泡表面流向靜壓在帽角處不存在這一跳躍現(xiàn)象。

    2.3.2=25~200時(shí)尾渦尺寸隨帽角的變化 為了定量分析氣泡受污染程度對(duì)其尾流的影響,Taneda[38]定義了如圖10所示用來(lái)描述尾流的相關(guān)參數(shù)(即尾渦長(zhǎng)度、渦中心位置和以及分離角),并預(yù)測(cè)了這些參數(shù)隨Reynolds數(shù)的變化情況。由于=250~500時(shí)流線呈明顯的三維特性,難以對(duì)尾渦準(zhǔn)確描述,且帽角≥157.5°時(shí),尾渦完全消失,故本文將用這些變量對(duì)=25~200范圍內(nèi)、帽角≤135°時(shí),氣泡的尾渦進(jìn)行深入分析,并給出其隨的變化曲線。

    圖11是帽角=0°、90°、135°時(shí)氣泡尾渦的相關(guān)尺寸:尾渦長(zhǎng)度/、渦中心位置/和/以及分離角隨Reynolds數(shù)的變化??梢钥闯?,帽角的變化只是改變了曲線上在數(shù)值上的大小,沒(méi)有對(duì)其隨Reynolds數(shù)的變化趨勢(shì)造成影響。

    2.3.3=250~500時(shí)尾渦性質(zhì)和尾渦尺寸隨帽角的變化 圖12為=250~500時(shí)不同帽角下的流線圖,可以看出:當(dāng)帽角=0°時(shí),流線雜亂無(wú)章;當(dāng)帽角帽角=90°時(shí),流線的混亂程度有所減小,對(duì)稱性有所增加;當(dāng)帽角帽角=135°時(shí),流線的三維特性已經(jīng)變得很弱,接近于二維特性;當(dāng)帽角≥157.5°時(shí),流線的三維特性消失,呈軸對(duì)稱形式。此時(shí),依舊能推斷出圖7中同一Reynolds數(shù)下干凈氣泡的阻力系數(shù)小于受污染氣泡的阻力系數(shù)正是由尾渦面積的減小和尾渦的消失造成的。氣泡受污染后尾渦面積會(huì)增大,使氣泡的阻力系數(shù)增加,進(jìn)而降低氣泡在靜置液中的上升速度。

    2.3.4=350~500時(shí)尾渦脫落情況隨帽角的變化 圓球繞流的研究結(jié)果已表明,當(dāng)280≤≤420時(shí),渦開(kāi)始脫落,但渦脫落的方向相同,仍然關(guān)于某一平面對(duì)稱[31-32];當(dāng)>420時(shí),尾渦開(kāi)始無(wú)規(guī)律脫落[32]。下文將在=350和=500這兩種情況下,通過(guò)氣泡尾部/=1處相差90°的兩監(jiān)測(cè)點(diǎn):點(diǎn)(/=1、/=0.2、/=0)以及點(diǎn)(/=1、/=0.2、/=0)處3個(gè)方向上的速度隨時(shí)間的變化來(lái)說(shuō)明氣泡界面變化對(duì)尾渦脫落情況的影響。兩個(gè)監(jiān)控點(diǎn)關(guān)于氣泡的位置如圖13所示,圖中表示氣泡直徑,坐標(biāo)原點(diǎn)位于氣泡中心。

    由于帽角0°~135°時(shí)氣泡尾部存在尾渦,為對(duì)=350時(shí)不同帽角下尾渦發(fā)生脫落時(shí)的對(duì)稱性進(jìn)行討論,圖14給出了=350,帽角=0°、90°和135°時(shí),氣泡尾部?jī)杀O(jiān)測(cè)點(diǎn)處流向速度隨時(shí)間的變化。可以看出:帽角=0°時(shí),兩監(jiān)測(cè)點(diǎn)處流向速度波動(dòng)的振幅相同,這說(shuō)明渦環(huán)脫落的方向相同,尾渦關(guān)于某一平面對(duì)稱,兩監(jiān)測(cè)點(diǎn)處流向速度波動(dòng)的頻率相同,這說(shuō)明渦環(huán)呈周期性脫落;帽角=90°時(shí),速度波動(dòng)的振幅明顯小于帽角=0°時(shí)速度波動(dòng)的振幅,這說(shuō)明此時(shí)脫落渦的強(qiáng)度已明顯減弱;帽角=135°時(shí),兩監(jiān)測(cè)點(diǎn)處流向速度沒(méi)有出現(xiàn)波動(dòng),為恒定值,這說(shuō)明污染程度較小的氣泡,其尾部不會(huì)產(chǎn)生渦脫落現(xiàn)象。

    由于帽角0°~135°時(shí)氣泡尾部存在尾渦,為對(duì)=500時(shí)不同帽角下尾渦發(fā)生脫落時(shí)的對(duì)稱性進(jìn)行討論,圖15給出了=500,帽角=0°,90°和135°時(shí),氣泡尾部?jī)杀O(jiān)測(cè)點(diǎn)處流向速度隨時(shí)間的變化。可以看出:當(dāng)帽角=0°時(shí),兩監(jiān)測(cè)點(diǎn)處流向速度波動(dòng)振幅不同,這說(shuō)明渦環(huán)脫落的方向不同,兩監(jiān)測(cè)點(diǎn)處流向速度波動(dòng)頻率不同,這說(shuō)明渦環(huán)隨機(jī)脫落,不再呈周期性脫落;當(dāng)帽角=90°時(shí),除了速度波動(dòng)的振幅減弱外,還可以看出速度波動(dòng)的隨機(jī)性有所減小,這表示污染程度減小后尾渦脫落的頻率與方向有趨向恒定的趨勢(shì)。與=350時(shí)相同,當(dāng)帽角=135°時(shí),兩監(jiān)測(cè)點(diǎn)處流向速度波動(dòng)消失,這表示污染程度較小的氣泡,其尾部不會(huì)產(chǎn)生渦脫落現(xiàn)象。

    3 結(jié) 論

    本文使用計(jì)算流體力學(xué)法研究了氣泡界面變化對(duì)其尾流的影響。借鑒圓球繞流和停滯帽模型,提出了一種模擬中等Reynolds數(shù)下受污染球形氣泡尾流的三維簡(jiǎn)化模型,氣泡界面污染程度取決于帽角的大小。研究發(fā)現(xiàn),對(duì)于不同污染程度的氣泡界面,其尾渦性質(zhì)與尺寸也不同,具體結(jié)論如下。

    (1)當(dāng)=25~200時(shí),隨著帽角的增加,尾渦的面積會(huì)不斷減小,當(dāng)帽角≥157.5°時(shí),尾渦將完全消失;對(duì)于尾渦尺寸,帽角的變化只會(huì)改變了尾渦相關(guān)尺寸(如尾渦長(zhǎng)度、尾渦中心位置和尾渦的分離角)上數(shù)值上的大小,不會(huì)對(duì)其隨Reynolds數(shù)的變化趨勢(shì)造成影響。

    (2)氣泡受污染后尾渦面積會(huì)增大,使氣泡的阻力系數(shù)增加,進(jìn)而降低氣泡在靜置液中的上升速度,并將節(jié)約下來(lái)的動(dòng)能被尾渦所消耗。

    (3)當(dāng)=250~500時(shí),氣泡表面污染程度的減小會(huì)使尾渦的三維特性減弱。

    (4)當(dāng)=350時(shí),氣泡表面污染程度的減小會(huì)使有序脫落的尾渦的強(qiáng)度減小并最終使其不發(fā)生脫落;當(dāng)=500時(shí),氣泡表面污染程度的減小會(huì)使無(wú)規(guī)律脫落的尾渦的無(wú)序性減弱并最終使其不發(fā)生脫落。

    符 號(hào) 說(shuō) 明

    CD——?dú)馀葑枇ο禂?shù) Cp——壓力系數(shù) D——圓球域直徑,m d——球形氣泡直徑,m Eo——E?tv?s數(shù) FD——阻力,N FDV——黏性力,N FDP——壓力梯度力,N h——?dú)馀葜行牡轿矞u中心的距離,m Lre——所需邊界層厚度,m Luse——?dú)馀荼砻娴谝粚泳W(wǎng)格的距離,m l——兩尾渦中心的距離,m Ma——Marangoni數(shù) Mo——Morton數(shù) nx——法向單位矢量的流向分量 p——靜壓,Pa R——?dú)怏w常數(shù),J·mol-1·K-1 Re——雷諾數(shù) S——?dú)馀荼砻娣e,m2 s——尾渦長(zhǎng)度,m T——溫度,K t——時(shí)間,s Δt——時(shí)間步長(zhǎng),s U——進(jìn)口液相速度,m·s-1 u——每個(gè)網(wǎng)格單元中的液相速度,m·s-1 ub——?dú)馀葸\(yùn)動(dòng)速度,m·s-1 α——方位角,(°) Γ¥——?dú)馀萁缑嫔系臉O限濃度,mol·kg-1 θ——帽角,(°) λ——?dú)馀葜行牡竭M(jìn)口或者壁面的距離,m μ——液相動(dòng)力黏度,Pa·s ρ——液相密度,kg·m-3 τx——表面切應(yīng)力流向分量,Pa φ——分離角,(°) 下角標(biāo) b——?dú)馀?re——需求值 use——使用值 x——流向

    References

    [1] TASOGLU S, DEMIRCI U, MURADOGLU M. The effect of soluble surfactant on the transient motion of a buoyancy-driven bubble[J]. Phys. Fluids, 2008, 20(4): 65-79.

    [2] ALKE A, BOTHE D. 3D numerical modeling of soluble surfactant at fluidic interfaces based on the VOF method[J]. FDMP, 2009, 5(4): 345-372.

    [3] FLECHENSTEIN S, BOTHE D. Simplified modeling of the influence of surfactants on the rise of bubbles in VOF-simulations[J]. Chem. Eng. Sci., 2013, 102(5): 514-523.

    [4] MURADOGLU M, TRYGGVASON G. Simulation of soluble surfactants in 3D multiphase flow[J]. J. Comput. Phys., 2014, 274: 737-757.

    [5] 裴明敬, 朱婷婷, 楊晶晶, 等. 環(huán)境友好型生物表面活性劑對(duì)浮選氣泡動(dòng)力學(xué)的影響研究[J]. 高校化學(xué)工程學(xué)報(bào), 2013, 27(6): 985-990. PEI M J, ZHU T T, YANG J J,. Study on effects of environment-friendly biosurfactants on bubble hydrodynamics in flotation column[J]. Journal of Chemical Engineering of Chinese Universities, 2013, 27(6): 985-990.

    [6] 劉艷艷, 李彥鵬, 朱婷婷. 表面活性劑對(duì)中尺度氣泡形狀及速度的調(diào)控研究[J]. 西安交通大學(xué)學(xué)報(bào), 2011, 45(10): 93-97. LIU Y Y, LI Y P, ZHU T T. Study on modulating shape and velocity of meso-scale bubble using surfactants[J]. Journal of Xi’an Jiaotong University , 2011, 45(10): 93-97.

    [7] BHAGA D, WEDER M E. Bubbles in viscous liquids: shapes, wakes and velocities[J]. J. Fluid Mech., 1981, 105: 61-85.

    [8] 倪明玖. 浮力作用下上升氣泡的變形和駐渦形成機(jī)理研究[J]. 工程熱物理學(xué)報(bào), 2009, 30(1): 76-80. NI M J. Bubble rising driven by buoyancy with deformation and standing vortex[J]. Journal of Engineering Thermophysics, 2009, 30(1): 76-80.

    [9] TRIPATHI M K, SAHU K C, GOVINDARAJAN R. Dynamics of an initially spherical bubble rising in quiescent liquid[J]. Nat. Commun., 2015, 6: 6268.

    [10] GUMULYA M, JOSHI J B, UTIKAR R P,. Bubbles in viscous liquids: time dependent behaviour and wake characteristics[J]. Chem. Eng. Sci., 2016 , 144: 298-309.

    [11] VELDHUIS C, BIESHEUWEL A, WIJNGAARDEN L V. Shape oscillations on bubbles rising in clean and in tap water[J]. Phys. Fluids, 2008, 20(4): 040705.

    [12] SAITO T, TORIU M. Effects of a bubble and the surrounding liquid motions on the instantaneous mass transfer across the gas-liquid interface[J]. Chem. Eng. J., 2015, 265: 164-175.

    [13] HUANG J, SAITO T. Influence of bubble-surface contamination on instantaneous mass transfer[J]. Chem. Eng. Technol., 2015, 38(11): 1947-1954.

    [14] HUANG J, SAITO T. Influences of gas-liquid interface contamination on bubble motions, bubble wakes, and instantaneous mass transfer[J]. Chem. Eng. Sci., 2017, 157: 182-199.

    [15] WANG Y P. A theoretical study of bubble motion in surfactant solutions[D]. Newark: New Jersey Institute of Technology, 1999.

    [16] SADHAL S S, JOHNSON R E. Stokes flow past bubbles and drops partially coated with thin films. Part 1. Stagnant cap of surfactant film-exact solution[J]. J. Fluid Mech., 1983, 126: 237-250.

    [17] FDHILA R B, DUINEVELD P C. The effect of surfactant on the rise of a spherical bubble at high Reynolds and Peclet numbers[J]. Phys. Fluids, 1996, 8(2): 310-321.

    [18] CUENOT B, MAGNAUDET J, SPENNATO B. The effects of slightly soluble surfactants on the flow around a spherical bubble[J]. J. Fluid Mech., 1997, 339: 25-53.

    [19] PALAPARTHI R, PAPAGEORGIOU D. Theory and experiments on the stagnant cap regime in the motion of spherical surfactant-laden bubbles[J]. J. Fluid Mech., 2006, 559: 1-44.

    [20] DUKHINA S S, KOVALCHUKB V I, GOCHEVCD G G,. Dynamics of rear stagnant cap formation at the surface of spherical bubbles in surfactant solutions at large Reynolds numbers under conditions of small Marangoni number and slow sorption kinetics[J]. Adv. Colloid Interface Sci., 2015, 222: 260-274.

    [21] DUKHINA S S, LOTFIBC M, KOVALCHUKB V I,. Dynamics of rear stagnant cap formation at the surface of rising bubbles in surfactant solutions at large Reynolds and Marangoni numbers and for slow sorption kinetics[J]. Colloid Surface A, 2016, 492: 127-137.

    [22] FLECHENSTEIN S, BOTHE D. Simplified modeling of the influence of surfactants on the rise of bubbles in VOF-simulations[J]. Chem. Eng. Sci., 2013, 102(5): 514-523.

    [23] GANESAN S, TOBISKA L. Arbitrary Lagrangian-Eulerian finite-element method for computation of two-phase flows with soluble surfactants[J]. J. Comput. Phys., 2012, 231: 3685-3702.

    [24] FUJIOKA H. A continuum model of interfacial surfactant transport for particle methods[J]. J. Comput. Phys., 2013, 234: 280-294.

    [25] KHATRI S, TORNBERG A K. An embedded boundary method for soluble surfactants with interface tracking for two-phase flows[J]. J. Comput. Phys., 2014, 256: 768-790.

    [26] MITTAL R, NAJJAR F M. Vortex dynamics in the sphere wake[C]// 30th Fluid Dynamics Conference, Fluid Dynamics and Co-located Conferences, Norfolk: AIAA J., 1999: 1-8.

    [27] TOMBOULIDES A G, ORSZAG S A. Numerical investigation of transitional and weak turbulent flow past a sphere[J]. J. Fluid Mech., 2000, 416: 45-73.

    [28] KIM H J, DURBIN P A. Observations of the frequencies in a sphere wake and of drag increase by acoustic excitation[J]. Phys. Fluids, 1988, 31(11): 3260-3265.

    [29] 張鳴遠(yuǎn), 景思睿, 李國(guó)君. 高等工程流體力學(xué)[M]. 北京: 高等教育出版社, 2012. ZHANG M Y, JING S R, LI G J. Advanced Engineering Fluid dynamics[M]. Beijing: Higher Education Press, 2012.

    [30] CROWE C T, SCHWARZKOPF J D, SOMMERFELD M,. Multiphase Flows with Droplets and Particles[M]. Boca Raton: Chemical Rubber Company Press, 2012.

    [31] NATARAJAN R, ACRIVOS A. The instability of the steady flow past spheres and disks[J]. J. Fluid Mech., 1993, 254(1): 323-344.

    [32] SAKAMOTO H, HANJU H. The formation mechanism and shedding frequency of vortices from a sphere in uniform shear flow[J]. J. Fluid Mech., 1995, 287: 151-171.

    [33] MAGARVEY R H, MACLATCHY C S. Vortices in sphere wakes[J]. Can. J. Phys., 1965, 43: 1649-1656.

    [34] MEI R, KLAUSNER J F, LAWRENCE C J. A note on the history force on a spherical bubble at finite Reynolds number[J]. Phys. Fluids, 1994, 6(1): 418-420.

    [35] MEI R. History force on a sphere due to a step change in the free-stream velocity[J]. Int. J. Multiphase Flow, 1993, 19: 509-525.

    [36] TASOGLU S, DEMIRCI U, MURADOGLU M. The effect of soluble surfactant on the transient motion of a buoyancy-driven bubble[J]. Phys. Fluids, 2008, 20(4): 65-79.

    [37] TAKAGI S, OGASAWARA T, MATSUMOTO Y. The effect of surfactant on the multiscale structure of bubbly flows[J]. Phil. Trans. R. Soc. A, 2008, 366: 2117-2129.

    [38] TANEDA S. Experimental investigation of the wake behind a sphere at low Reynolds number[J]. J. Phys. Soc. Jpn., 1956, 11(10): 1104-110.

    Influence of interface change for spherical bubble on vortex characteristic and size

    FEI Yang, PANG Mingjun

    (Jiangsu Key Laboratory of Green Process Equipment, School of Mechanical Engineering, Changzhou University, Changzhou 213164, Jiangsu, China)

    The numerical method is employed to investigate the influence of contaminated degree of bubble surface on its wakes for the spherical bubble under moderate Reynolds number (25≤≤500). By referencing the flow past a sphere and the stagnant cap model, one kind of three-dimensional model for contaminated spherical bubble under moderate Reynolds number is proposed.The interface contaminated degree is dependent on the magnitude of the cap angle. The larger the cap angel is, the slighter the bubble interface pollution is. The present results show that, for 25≤≤200, the magnitudes of the vortex length, the vortex center position distance to the bubble rear and the separation angle decrease with the decrease of the bubble surface contaminated degree but the distribution trends of those parameters against the Reynolds number are similar; for 250≤≤500, the decrease of the bubble surface contaminated degree weakens the three-dimensional property of vortexes, reduces the strength of orderly shedding vortexes until the shedding phenomenon disappears at=350, and reduces the disorder of shedding vortexes until the shedding phenomenon disappears too at=500.

    bubble; interface; wake; laminar flow; hydrodynamics

    10.11949/j.issn.0438-1157.20170290

    O 357.1

    A

    0438—1157(2017)09—3409—11

    2017-03-23收到初稿,2017-06-27收到修改稿。

    龐明軍。

    費(fèi)洋(1994—),男,碩士研究生。

    國(guó)家自然科學(xué)基金項(xiàng)目(51376026);江蘇省青藍(lán)工程項(xiàng)目。

    2017-03-23.

    Prof. PANG Mingjun, pangmj@cczu.edu.cn

    supported by the National Natural Science Foundation of China (51376026) and the Qinglan Project of Jiangsu Province.

    猜你喜歡
    尾渦尾流氣泡
    檸檬氣泡水
    欣漾(2024年2期)2024-04-27 15:19:49
    不同B-V頻率下的飛機(jī)尾渦數(shù)值模擬研究
    高空巡航階段的飛機(jī)尾渦流場(chǎng)演化特性研究
    SIAU詩(shī)杭便攜式氣泡水杯
    新潮電子(2021年7期)2021-08-14 15:53:12
    浮法玻璃氣泡的預(yù)防和控制對(duì)策
    冰凍氣泡
    基于激光雷達(dá)回波的動(dòng)態(tài)尾渦特征參數(shù)計(jì)算
    干擾板作用下飛機(jī)尾渦流場(chǎng)近地演變機(jī)理研究
    飛機(jī)尾流的散射特性與探測(cè)技術(shù)綜述
    錐形流量計(jì)尾流流場(chǎng)分析
    日本一二三区视频观看| 亚洲天堂av无毛| 国产成人午夜福利电影在线观看| 久久久久视频综合| 韩国高清视频一区二区三区| 青春草亚洲视频在线观看| 直男gayav资源| 国产综合精华液| 视频区图区小说| 午夜日本视频在线| h日本视频在线播放| 久久久欧美国产精品| 欧美人与善性xxx| 五月开心婷婷网| 国产成人精品久久久久久| 亚洲精品自拍成人| 精华霜和精华液先用哪个| 国产亚洲最大av| 2021少妇久久久久久久久久久| 欧美日韩亚洲高清精品| 国产精品.久久久| 一区二区三区精品91| 国产亚洲午夜精品一区二区久久| 成年人午夜在线观看视频| 尤物成人国产欧美一区二区三区| 亚洲色图av天堂| av国产免费在线观看| 亚洲av电影在线观看一区二区三区| 国产片特级美女逼逼视频| 精品亚洲成a人片在线观看 | 老熟女久久久| 自拍偷自拍亚洲精品老妇| 国产高清不卡午夜福利| 国产熟女欧美一区二区| 亚洲第一区二区三区不卡| 国产高潮美女av| 亚洲精品,欧美精品| 一二三四中文在线观看免费高清| 国产黄片视频在线免费观看| 欧美日韩国产mv在线观看视频 | 日韩在线高清观看一区二区三区| 日韩,欧美,国产一区二区三区| 欧美日韩综合久久久久久| 交换朋友夫妻互换小说| xxx大片免费视频| 精品久久久精品久久久| 一个人免费看片子| 午夜福利在线观看免费完整高清在| 又大又黄又爽视频免费| 六月丁香七月| 亚洲图色成人| 全区人妻精品视频| 亚洲av福利一区| 卡戴珊不雅视频在线播放| 国产亚洲一区二区精品| 亚洲在久久综合| 国产精品福利在线免费观看| 成人特级av手机在线观看| 国产精品嫩草影院av在线观看| 国产乱人偷精品视频| 国产在线免费精品| 26uuu在线亚洲综合色| 综合色丁香网| 成年av动漫网址| 一个人看视频在线观看www免费| 国产极品天堂在线| 噜噜噜噜噜久久久久久91| 少妇 在线观看| 久久这里有精品视频免费| 国产女主播在线喷水免费视频网站| 亚洲精品色激情综合| 午夜福利影视在线免费观看| av线在线观看网站| 肉色欧美久久久久久久蜜桃| 高清日韩中文字幕在线| videossex国产| 日日啪夜夜爽| 男女免费视频国产| 有码 亚洲区| 欧美xxxx性猛交bbbb| 青春草视频在线免费观看| 亚洲欧美精品自产自拍| 久久ye,这里只有精品| 国产乱来视频区| 亚洲无线观看免费| 视频区图区小说| 欧美日韩国产mv在线观看视频 | 97超视频在线观看视频| 国产黄色视频一区二区在线观看| 久久这里有精品视频免费| 精品人妻偷拍中文字幕| 色视频在线一区二区三区| 99热这里只有是精品在线观看| videossex国产| 久久久久久久久久久免费av| 国产免费一级a男人的天堂| 22中文网久久字幕| 少妇被粗大猛烈的视频| 日韩中文字幕视频在线看片 | 午夜老司机福利剧场| 一级二级三级毛片免费看| 久久97久久精品| 欧美极品一区二区三区四区| 亚洲色图av天堂| 久久久久久久久大av| 毛片女人毛片| 97超碰精品成人国产| 日本av免费视频播放| 超碰av人人做人人爽久久| 国产v大片淫在线免费观看| 国产精品99久久久久久久久| 最近中文字幕2019免费版| 菩萨蛮人人尽说江南好唐韦庄| 国产亚洲午夜精品一区二区久久| 高清欧美精品videossex| 欧美一级a爱片免费观看看| 在线观看av片永久免费下载| 国产精品精品国产色婷婷| 啦啦啦啦在线视频资源| 国产精品一区二区三区四区免费观看| 国产亚洲91精品色在线| 国产69精品久久久久777片| 久久鲁丝午夜福利片| 少妇人妻 视频| 插阴视频在线观看视频| 国产精品国产av在线观看| 黄片无遮挡物在线观看| 国产高清有码在线观看视频| 内地一区二区视频在线| 国产男女超爽视频在线观看| 一区二区三区精品91| 国产精品免费大片| 国产高清国产精品国产三级 | 最近中文字幕2019免费版| 草草在线视频免费看| 男女国产视频网站| 亚洲欧美中文字幕日韩二区| 欧美高清成人免费视频www| 亚洲av日韩在线播放| 欧美少妇被猛烈插入视频| 国内精品宾馆在线| 国产精品人妻久久久久久| 成年美女黄网站色视频大全免费 | 久久婷婷青草| 这个男人来自地球电影免费观看 | 夫妻性生交免费视频一级片| 在线观看av片永久免费下载| 午夜视频国产福利| 亚洲av国产av综合av卡| 伊人久久精品亚洲午夜| 亚洲精品久久久久久婷婷小说| 又大又黄又爽视频免费| 亚洲人与动物交配视频| 国产又色又爽无遮挡免| 国产伦精品一区二区三区视频9| 亚洲aⅴ乱码一区二区在线播放| 超碰97精品在线观看| 精品久久久噜噜| 汤姆久久久久久久影院中文字幕| 中文精品一卡2卡3卡4更新| 成人国产av品久久久| 嫩草影院新地址| 久久久色成人| 国产一区二区三区综合在线观看 | 亚洲四区av| 中文字幕亚洲精品专区| 亚洲欧美日韩东京热| 日本wwww免费看| av在线播放精品| 亚洲国产最新在线播放| 青春草视频在线免费观看| 国产av码专区亚洲av| 国产国拍精品亚洲av在线观看| 精品久久久精品久久久| 高清在线视频一区二区三区| 成人毛片60女人毛片免费| 日韩大片免费观看网站| 国产精品一二三区在线看| av专区在线播放| 少妇人妻一区二区三区视频| 欧美一区二区亚洲| 性色av一级| 夫妻性生交免费视频一级片| 搡女人真爽免费视频火全软件| av不卡在线播放| av黄色大香蕉| 国产女主播在线喷水免费视频网站| 国产精品精品国产色婷婷| 国产乱来视频区| 啦啦啦啦在线视频资源| 大香蕉久久网| 精品久久久久久电影网| 国产视频内射| 久热这里只有精品99| 日本与韩国留学比较| 成人影院久久| 亚洲精品一区蜜桃| 身体一侧抽搐| 久久久久久久精品精品| 伊人久久国产一区二区| 国产真实伦视频高清在线观看| 中国国产av一级| 久久精品国产a三级三级三级| 大香蕉久久网| 日韩av免费高清视频| 狂野欧美激情性bbbbbb| 久久这里有精品视频免费| av在线蜜桃| 一边亲一边摸免费视频| 人妻少妇偷人精品九色| 黄色日韩在线| 七月丁香在线播放| 欧美极品一区二区三区四区| 尾随美女入室| 国产成人午夜福利电影在线观看| 欧美日韩精品成人综合77777| 婷婷色麻豆天堂久久| 少妇人妻久久综合中文| 精品一区在线观看国产| 亚洲婷婷狠狠爱综合网| 啦啦啦视频在线资源免费观看| 日韩欧美精品免费久久| 80岁老熟妇乱子伦牲交| 岛国毛片在线播放| 国产乱人偷精品视频| 十分钟在线观看高清视频www | 国产精品爽爽va在线观看网站| 街头女战士在线观看网站| 免费播放大片免费观看视频在线观看| 美女国产视频在线观看| 有码 亚洲区| 我要看黄色一级片免费的| 尾随美女入室| 久久国产乱子免费精品| 国产老妇伦熟女老妇高清| 99久久精品一区二区三区| 国产永久视频网站| 天天躁日日操中文字幕| 欧美一级a爱片免费观看看| 日本色播在线视频| 亚洲av免费高清在线观看| 最后的刺客免费高清国语| 精品一区二区三区视频在线| 久久国产精品男人的天堂亚洲 | 亚洲国产精品一区三区| 高清毛片免费看| 内射极品少妇av片p| 最黄视频免费看| 麻豆成人午夜福利视频| 99国产精品免费福利视频| 免费高清在线观看视频在线观看| 深爱激情五月婷婷| 亚洲精品中文字幕在线视频 | 少妇的逼水好多| 男女无遮挡免费网站观看| 99视频精品全部免费 在线| 黑人高潮一二区| 成人影院久久| 精品午夜福利在线看| 日日摸夜夜添夜夜爱| av专区在线播放| 久久久久久久久久成人| 99热国产这里只有精品6| 97在线视频观看| 欧美高清性xxxxhd video| 成人免费观看视频高清| 亚洲精品国产成人久久av| 免费看av在线观看网站| 国产白丝娇喘喷水9色精品| 国产成人免费观看mmmm| 麻豆成人午夜福利视频| 99精国产麻豆久久婷婷| 国产精品人妻久久久影院| 婷婷色av中文字幕| 免费av不卡在线播放| 菩萨蛮人人尽说江南好唐韦庄| 国产69精品久久久久777片| 国产精品一区二区三区四区免费观看| 国产精品一区二区在线不卡| 性色av一级| av在线蜜桃| 一级二级三级毛片免费看| 国产精品久久久久久久久免| 亚洲精品视频女| 视频中文字幕在线观看| 99久久人妻综合| 日本免费在线观看一区| 亚洲,一卡二卡三卡| 亚洲av在线观看美女高潮| 精品人妻视频免费看| 嫩草影院入口| 久久久久国产网址| 国产成人91sexporn| 亚洲av综合色区一区| 伦理电影大哥的女人| 舔av片在线| 国产爽快片一区二区三区| 2022亚洲国产成人精品| 国模一区二区三区四区视频| 麻豆成人av视频| 国产精品秋霞免费鲁丝片| 美女内射精品一级片tv| 色综合色国产| 一边亲一边摸免费视频| av在线老鸭窝| 亚洲精品自拍成人| 在线免费观看不下载黄p国产| 伦精品一区二区三区| 一级a做视频免费观看| 麻豆精品久久久久久蜜桃| 在线看a的网站| 99九九线精品视频在线观看视频| 久久久精品免费免费高清| 欧美精品一区二区免费开放| 99热6这里只有精品| 人妻制服诱惑在线中文字幕| 草草在线视频免费看| 欧美日韩亚洲高清精品| 中文天堂在线官网| 国产 精品1| 欧美+日韩+精品| 黑人猛操日本美女一级片| 中文字幕亚洲精品专区| 久久精品国产自在天天线| 最近手机中文字幕大全| 免费观看av网站的网址| 校园人妻丝袜中文字幕| 亚洲最大成人中文| 自拍欧美九色日韩亚洲蝌蚪91 | 亚洲精品第二区| 国产一区二区三区综合在线观看 | 久久久久久伊人网av| 久久久久精品性色| 熟女电影av网| 丰满迷人的少妇在线观看| 亚洲三级黄色毛片| 亚洲熟女精品中文字幕| 韩国高清视频一区二区三区| 欧美日韩综合久久久久久| 欧美成人a在线观看| 中文字幕制服av| 国产精品偷伦视频观看了| 久久午夜福利片| 久久精品国产亚洲av涩爱| 免费不卡的大黄色大毛片视频在线观看| 高清在线视频一区二区三区| 22中文网久久字幕| 天天躁日日操中文字幕| 不卡视频在线观看欧美| 美女国产视频在线观看| 免费在线观看成人毛片| 自拍偷自拍亚洲精品老妇| 国产男人的电影天堂91| 大码成人一级视频| 九草在线视频观看| 国产男女超爽视频在线观看| 黄色日韩在线| 天天躁夜夜躁狠狠久久av| 国产精品av视频在线免费观看| 国产精品三级大全| 熟女人妻精品中文字幕| 国产av精品麻豆| 亚洲色图综合在线观看| 99热这里只有是精品50| 亚洲国产欧美在线一区| 蜜桃久久精品国产亚洲av| 婷婷色综合大香蕉| 熟女人妻精品中文字幕| 亚洲精品自拍成人| 男人和女人高潮做爰伦理| 身体一侧抽搐| 精品久久久久久久末码| 欧美日韩综合久久久久久| 精品一区在线观看国产| 秋霞伦理黄片| 亚洲av成人精品一二三区| 欧美极品一区二区三区四区| 日韩av不卡免费在线播放| 在线 av 中文字幕| 久久久久久久久久久免费av| 免费观看性生交大片5| tube8黄色片| 最近最新中文字幕免费大全7| 啦啦啦啦在线视频资源| 91久久精品电影网| 人妻少妇偷人精品九色| 91久久精品电影网| 国产在视频线精品| 精品国产一区二区三区久久久樱花 | 人妻 亚洲 视频| 深夜a级毛片| 妹子高潮喷水视频| av.在线天堂| 亚洲av免费高清在线观看| 久久鲁丝午夜福利片| 91久久精品国产一区二区成人| 赤兔流量卡办理| 国产高清国产精品国产三级 | 日本黄色日本黄色录像| 这个男人来自地球电影免费观看 | 伦精品一区二区三区| 波野结衣二区三区在线| 少妇人妻 视频| 97热精品久久久久久| 嫩草影院新地址| 久久久久久久大尺度免费视频| 纵有疾风起免费观看全集完整版| 国产亚洲午夜精品一区二区久久| 久久久久性生活片| 99视频精品全部免费 在线| 亚洲精品色激情综合| 男人和女人高潮做爰伦理| 在线精品无人区一区二区三 | 午夜福利影视在线免费观看| 黄色欧美视频在线观看| 国产 一区精品| 最后的刺客免费高清国语| 久久人妻熟女aⅴ| 91久久精品国产一区二区成人| 免费av不卡在线播放| 久久久精品94久久精品| 91狼人影院| 精品久久国产蜜桃| 国产精品不卡视频一区二区| 最近最新中文字幕免费大全7| av在线老鸭窝| 久久久成人免费电影| 老师上课跳d突然被开到最大视频| 成人黄色视频免费在线看| 国产精品蜜桃在线观看| 精品人妻视频免费看| 大香蕉97超碰在线| 亚洲成人av在线免费| 夫妻午夜视频| 少妇精品久久久久久久| 一级a做视频免费观看| 老女人水多毛片| 久久人人爽人人爽人人片va| 亚洲av在线观看美女高潮| 成人国产av品久久久| 国产精品一区二区性色av| 亚洲欧美清纯卡通| 国产高潮美女av| 国产中年淑女户外野战色| 国产亚洲av片在线观看秒播厂| 十分钟在线观看高清视频www | 国产精品国产三级国产av玫瑰| 国产免费又黄又爽又色| 搡女人真爽免费视频火全软件| 精品人妻一区二区三区麻豆| 五月伊人婷婷丁香| av在线老鸭窝| 一区在线观看完整版| 欧美高清性xxxxhd video| 亚洲精品自拍成人| 欧美日韩亚洲高清精品| 国产真实伦视频高清在线观看| 久久99热6这里只有精品| 亚洲久久久国产精品| 一级a做视频免费观看| 最近2019中文字幕mv第一页| 中国美白少妇内射xxxbb| 美女主播在线视频| 亚洲美女搞黄在线观看| 人人妻人人添人人爽欧美一区卜 | 色综合色国产| 精品一区在线观看国产| 看免费成人av毛片| 18禁在线播放成人免费| 国产男女超爽视频在线观看| 国精品久久久久久国模美| 国产av码专区亚洲av| 毛片一级片免费看久久久久| 天堂中文最新版在线下载| 亚洲国产最新在线播放| 久久精品人妻少妇| 在线免费十八禁| 欧美亚洲 丝袜 人妻 在线| 国产片特级美女逼逼视频| 国产黄频视频在线观看| 日韩一区二区三区影片| 22中文网久久字幕| 最近中文字幕高清免费大全6| 夫妻午夜视频| 亚洲av电影在线观看一区二区三区| 日日啪夜夜爽| 高清日韩中文字幕在线| 一级黄片播放器| av.在线天堂| 爱豆传媒免费全集在线观看| 男人和女人高潮做爰伦理| 中文资源天堂在线| 丰满乱子伦码专区| 99热6这里只有精品| 国产av码专区亚洲av| 欧美xxxx黑人xx丫x性爽| 国产成人精品一,二区| 身体一侧抽搐| 国内精品宾馆在线| 人妻制服诱惑在线中文字幕| 在线观看国产h片| 91精品一卡2卡3卡4卡| 久久久亚洲精品成人影院| 精品久久久久久久久av| 亚洲国产毛片av蜜桃av| 欧美精品一区二区免费开放| 777米奇影视久久| 在线观看一区二区三区激情| 国产真实伦视频高清在线观看| 欧美一区二区亚洲| 日韩一区二区视频免费看| av播播在线观看一区| 狂野欧美白嫩少妇大欣赏| 久久久久久久久久成人| 十分钟在线观看高清视频www | 性高湖久久久久久久久免费观看| 国产成人91sexporn| 亚洲va在线va天堂va国产| 国产伦理片在线播放av一区| 一本—道久久a久久精品蜜桃钙片| 高清在线视频一区二区三区| 午夜福利在线观看免费完整高清在| 亚洲精品国产成人久久av| 大又大粗又爽又黄少妇毛片口| 亚洲欧美日韩无卡精品| 美女内射精品一级片tv| 亚洲色图av天堂| 九色成人免费人妻av| 国产亚洲一区二区精品| 国产亚洲5aaaaa淫片| 卡戴珊不雅视频在线播放| 成人高潮视频无遮挡免费网站| 成人国产av品久久久| 最近最新中文字幕大全电影3| 精品人妻视频免费看| 亚洲av成人精品一区久久| 日韩大片免费观看网站| 高清视频免费观看一区二区| 久久久久久久精品精品| 我的老师免费观看完整版| 亚洲精品成人av观看孕妇| 日本与韩国留学比较| 欧美精品亚洲一区二区| 欧美精品人与动牲交sv欧美| 亚洲av不卡在线观看| 国内揄拍国产精品人妻在线| 另类亚洲欧美激情| 欧美性感艳星| 欧美xxⅹ黑人| videossex国产| 免费人妻精品一区二区三区视频| 亚洲国产色片| 亚洲av成人精品一二三区| 精品99又大又爽又粗少妇毛片| 亚洲国产欧美人成| 国产黄片美女视频| 午夜老司机福利剧场| 婷婷色麻豆天堂久久| 在线免费观看不下载黄p国产| 色婷婷av一区二区三区视频| av黄色大香蕉| 在线观看美女被高潮喷水网站| 国产成人freesex在线| 国产精品熟女久久久久浪| 久久人妻熟女aⅴ| 国产v大片淫在线免费观看| 国产高潮美女av| 内地一区二区视频在线| 亚洲欧美日韩无卡精品| 亚洲国产成人一精品久久久| 日本猛色少妇xxxxx猛交久久| av视频免费观看在线观看| 午夜福利视频精品| 波野结衣二区三区在线| 成年免费大片在线观看| 男女边摸边吃奶| 久久精品国产a三级三级三级| 欧美精品一区二区免费开放| 国产高清三级在线| av一本久久久久| 一个人免费看片子| 国产精品偷伦视频观看了| 韩国高清视频一区二区三区| 亚洲av欧美aⅴ国产| 国产乱来视频区| 亚洲欧美日韩东京热| 99热网站在线观看| 久久精品国产a三级三级三级| 午夜福利在线在线| 2021少妇久久久久久久久久久| 1000部很黄的大片| 亚洲欧美日韩卡通动漫| 成人综合一区亚洲| 91久久精品国产一区二区三区| 又黄又爽又刺激的免费视频.| 亚洲一级一片aⅴ在线观看| 欧美日韩亚洲高清精品| 22中文网久久字幕| 热re99久久精品国产66热6| 午夜福利高清视频| 国产免费一区二区三区四区乱码| 日本-黄色视频高清免费观看| 少妇猛男粗大的猛烈进出视频| 高清午夜精品一区二区三区| 六月丁香七月| 亚洲电影在线观看av| 丰满少妇做爰视频| 国产免费一级a男人的天堂| 视频中文字幕在线观看| www.av在线官网国产| 久久久久性生活片| 亚洲av欧美aⅴ国产| 中文字幕人妻熟人妻熟丝袜美| av国产精品久久久久影院| 人人妻人人澡人人爽人人夜夜|