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

    底部加熱肥皂泡上準(zhǔn)二維湍流的數(shù)值模擬*

    2022-11-09 05:52:56賀嘯秋熊永亮彭澤瑞

    賀嘯秋,熊永亮,徐 順,彭澤瑞,陳 波

    (華中科技大學(xué) 航空航天學(xué)院 力學(xué)系,武漢 430074)

    引言

    湍流運(yùn)動(dòng)廣泛存在于自然界中,從星系中的星際介質(zhì)運(yùn)動(dòng),到杯子中被攪拌的咖啡流動(dòng),雖然時(shí)空尺度差距懸殊但都屬于湍流流動(dòng)[1].流體的湍流運(yùn)動(dòng)對傳熱傳質(zhì)具有重要的影響[2].當(dāng)湍流流場中流體的運(yùn)動(dòng)僅具有兩個(gè)獨(dú)立的空間自由度時(shí),稱之為二維湍流[3-6].相比三維湍流,二維湍流的結(jié)構(gòu)更易于觀察與分析,且在直接數(shù)值模擬(DNS)所需要的計(jì)算資源顯著減小的同時(shí)又保留著湍流的許多關(guān)鍵特征[3].另一方面,二維湍流并不是三維湍流進(jìn)行簡單維度壓縮的結(jié)果,它也具有與三維湍流完全不同的顯著差別[3].其中最核心的差別是最早由Kraichnan[7]、Leith[8]和Batchelor[9]發(fā)現(xiàn)的反向能量串級(jí),即在充分發(fā)展的二維湍流中,在特定尺度范圍內(nèi)存在著能量從尺寸較小的渦向尺寸較大的渦反向傳遞的現(xiàn)象.在研究海洋或者大氣中的大尺度運(yùn)動(dòng)時(shí),海洋或大氣的側(cè)向尺寸遠(yuǎn)大于高度方向的尺寸,在研究較大尺度的流動(dòng)問題時(shí)可以忽略掉流體在高度方向的運(yùn)動(dòng)而作二維流動(dòng)處理[10-11].因而,在地球物理與行星物理學(xué)中,二維湍流具有重要的研究價(jià)值[12].

    嚴(yán)格意義下的二維湍流在實(shí)驗(yàn)室或自然界中并不存在[3],研究二維湍流主要依賴于DNS[3]或者在實(shí)驗(yàn)條件下借助肥皂膜及電磁力抑制形成準(zhǔn)二維流動(dòng)[13].眾所周知,肥皂薄膜的厚度只有數(shù)微米,但其側(cè)向尺寸可以達(dá)到數(shù)十厘米以上,肥皂薄膜被廣泛用于準(zhǔn)二維湍流實(shí)驗(yàn)[12].Kellay[13]在實(shí)驗(yàn)中研究了加熱肥皂泡中的準(zhǔn)二維湍流.實(shí)驗(yàn)中一個(gè)半球形的肥皂泡被固定在基座上,基座以水浴的方法對肥皂泡底部進(jìn)行加熱,肥皂泡上的流體在浮力的驅(qū)動(dòng)下形成復(fù)雜的湍流運(yùn)動(dòng)[14-15].肥皂泡上的流體運(yùn)動(dòng)與大氣流動(dòng)都具有準(zhǔn)二維流動(dòng)特征,且具有相似的幾何特征,這使得學(xué)者以肥皂泡為簡化模型開展了大氣流動(dòng)中尺度結(jié)構(gòu)特征的研究[13-19].熱帶氣旋是大氣中的一種特殊中尺度結(jié)構(gòu),對人類社會(huì)具有重要的影響,是多種極端天氣與自然災(zāi)害的直接原因[10].Seychelles 等[14]在實(shí)驗(yàn)中觀察到肥皂泡上出現(xiàn)了能長時(shí)間獨(dú)立存在的島渦,島渦的運(yùn)動(dòng)規(guī)律與大氣中的熱帶氣旋類似,都具有超擴(kuò)散的特征[14,16].此外,島渦的強(qiáng)度隨時(shí)間變化的規(guī)律與大氣中的熱帶氣旋類似[15].基于這些特點(diǎn),Meuel 等[17]基于肥皂泡上島渦與颶風(fēng)的運(yùn)動(dòng)特征提出了預(yù)測大氣氣旋運(yùn)動(dòng)軌跡的模型.在湍流方面,研究者們發(fā)現(xiàn)肥皂泡上的流動(dòng)具有其特有的曲面準(zhǔn)二維熱對流特征,通過實(shí)驗(yàn)與DNS,可利用這一標(biāo)準(zhǔn)模型開展湍流流動(dòng)機(jī)理的研究與湍流理論的考察[16,19-21].例如通過對肥皂泡的旋轉(zhuǎn),發(fā)現(xiàn)了旋轉(zhuǎn)效應(yīng)產(chǎn)生了新的溫度小尺度脈動(dòng)標(biāo)度關(guān)系[18].因而研究肥皂泡上的湍流熱對流不僅具有重要的應(yīng)用價(jià)值,同時(shí)也可以拓展與豐富學(xué)界對湍流物理機(jī)制的認(rèn)識(shí).

    借助赤平極射投影法,Xiong 等[20]首先實(shí)現(xiàn)了具有半球面幾何的二維流場DNS.他們計(jì)算了Rayleigh 數(shù)為1×108的肥皂泡流場,由于缺乏能量耗散,所得到的計(jì)算結(jié)果并不穩(wěn)定[20].Meuel 等[15]通過實(shí)驗(yàn)研究了肥皂泡上的漩渦結(jié)構(gòu)特征,并與DNS 計(jì)算結(jié)果進(jìn)行了對比,發(fā)現(xiàn)兩者符合得較好.此后,Bruneau 等[19]在DNS 代碼中引入了能量耗散項(xiàng),使得計(jì)算得到的流場能夠達(dá)到統(tǒng)計(jì)穩(wěn)態(tài),計(jì)算獲得的流場與實(shí)驗(yàn)觀察到的肥皂泡上的流動(dòng)具有相同的特征,同時(shí)還發(fā)現(xiàn)DNS 計(jì)算結(jié)果對一定的外部耗散并不敏感,他們還通過DNS 結(jié)果發(fā)現(xiàn)了肥皂泡對流具有Bo59標(biāo)度律,這與實(shí)驗(yàn)中肥皂泡上觀察到的結(jié)果一致.最近,Meuel 等[18]同時(shí)使用實(shí)驗(yàn)與DNS 方法研究了旋轉(zhuǎn)肥皂泡上的溫度脈動(dòng),他們的研究中實(shí)驗(yàn)與DNS 獲得了非常相近的結(jié)果.經(jīng)過近十年的發(fā)展,肥皂泡上熱對流的DNS 已日趨成熟,相比實(shí)驗(yàn)方法,使用DNS 研究肥皂泡上的湍流熱對流也展露出較大的優(yōu)勢.首先,目前的實(shí)驗(yàn)手段在獲得肥皂泡上的空間特征結(jié)構(gòu)關(guān)聯(lián)上還具有一定難度,難以同時(shí)測量整個(gè)流場[13-14,17].對肥皂泡上的準(zhǔn)二維湍流的深入研究有必要獲得包含全部尺度特征的完整流場信息,而采用DNS 方法可以滿足這一需求[22].其次,實(shí)驗(yàn)方法也難以精確地在曲面幾何中測量如能量級(jí)串、結(jié)構(gòu)函數(shù)等湍流的小尺度特征[23-24].而DNS 獲得的流場包含全部尺度的特征,可以使用編程后處理的方法高效快捷地計(jì)算湍流的小尺度特征.另外,使用DNS 方法可以方便地在更大、更明確的參數(shù)空間內(nèi)研究肥皂泡曲面空間上的熱對流,并規(guī)避了如非Boussinesq 效應(yīng)、不均勻加熱邊界、材質(zhì)與外界擾動(dòng)等實(shí)驗(yàn)中常見的誤差來源.本文將詳細(xì)介紹這一套簡便求解肥皂泡上流場的并行計(jì)算DNS 方法,以及針對球面二維湍流的分析方法.

    文中通過該方法對多個(gè)高Rayleigh 數(shù)下的熱對流開展了DNS,并探索了肥皂泡上的湍流雙級(jí)串現(xiàn)象與小尺度脈動(dòng)的標(biāo)度關(guān)系.

    1 模型介紹與控制方程

    本文研究的是在底部加熱的半球形肥皂泡上的準(zhǔn)二維湍流.肥皂薄膜的厚度遠(yuǎn)小于肥皂泡的半徑R,可以忽略流體在肥皂薄膜厚度方向的運(yùn)動(dòng),使用二維半球面來近似肥皂泡的幾何外形.赤道是肥皂泡與基座接觸的邊界,幾何外形為一個(gè)半徑為R的圓周.首先,需要建立坐標(biāo)系描述肥皂泡上流體微元的空間位置.三維直角坐標(biāo)系Oxyz以二維半球面的球心O為原點(diǎn),ex,ey和ez分別為x,y和z方向上的單位基矢,其中ez方向與重力方向相反.在直角坐標(biāo)系下,肥皂泡用曲面方程x2+y2+z2=R2(z≥0)描述,赤道的曲線方程為x2+y2=R2.采用球坐標(biāo)系Or?θ,其中r≡R可方便地將該問題簡化為二維,球坐標(biāo)系的單位基矢用er,e?和eθ表示.與球坐標(biāo)類似,地理坐標(biāo)系?*Oθ*能更加直觀地描述肥皂泡上任意一點(diǎn)的空間位置,經(jīng)度定義為?*=?,緯度定義為θ*=90?-θ.肥皂泡頂部緯度為 90?(θ*=90?),而在赤道緯度為 0?(θ*=0?),地理坐標(biāo)系的基矢為e?*=e?,eθ*=-eθ.對于肥皂泡上的速度場u有如下定義:

    肥皂泡上流場的唯一邊界是赤道,邊界條件為恒溫?zé)o滑移(T=T0,u=0)邊界條件.若肥皂泡的初始溫度遠(yuǎn)小于肥皂泡的赤道溫度,隨著時(shí)間的推進(jìn),在赤道附近的流體,通過赤道從外界吸收熱量,使其溫度升高密度降低.此時(shí),流體在浮力的驅(qū)動(dòng)下向肥皂泡頂部運(yùn)動(dòng),肥皂泡中產(chǎn)生熱對流.數(shù)值模擬中,Boussinesq 近似常被用來描述流體溫度與密度的關(guān)系.當(dāng)流體運(yùn)動(dòng)至遠(yuǎn)離赤道的位置時(shí),熱量向周圍的低溫空氣擴(kuò)散,導(dǎo)致溫度下降密度逐漸上升,受到的浮力逐漸減弱.觀察實(shí)驗(yàn)中的肥皂泡流場可以發(fā)現(xiàn),由于肥皂泡既向冷空氣散熱又與冷空氣產(chǎn)生摩擦,肥皂泡上的湍流熱對流不斷持續(xù)并達(dá)到了統(tǒng)計(jì)平衡態(tài)[13-18].為了使DNS 中肥皂泡進(jìn)入統(tǒng)計(jì)平衡態(tài),外冷卻項(xiàng)ST和外摩擦項(xiàng)Fu被分別加入到能量方程與動(dòng)量方程中,用于描述流體與環(huán)境冷空氣的散熱與摩擦,S為外冷卻系數(shù),而F是外摩擦因數(shù).以往的研究表明,外冷卻項(xiàng)和外摩擦項(xiàng)可以較好地描述肥皂泡與實(shí)驗(yàn)環(huán)境之間的能量交換[19-21].因此,流體在肥皂泡上運(yùn)動(dòng)的控制方程可寫為Ra建立了流體受到的浮力與黏性力之間的一個(gè)關(guān)系,Pr則反映了動(dòng)量傳輸及熱量傳輸速率之比,對于特定的物質(zhì),Pr是流體本身的屬性,此處為一定值.

    考慮到需要在球坐標(biāo)系下求解控制方程(4),則在肥皂泡的頂點(diǎn)處引入了奇異性.而在直角坐標(biāo)系下無法對這一曲面進(jìn)行二維描述,采用三維描述并進(jìn)行離散則降低了曲面的光滑性.通過赤平極射投影法引入的投影空間坐標(biāo)系x?Oy?可以規(guī)避以上全部缺點(diǎn).首先,定義肥皂泡的頂點(diǎn)關(guān)于赤道平面的鏡像點(diǎn)為B點(diǎn),B點(diǎn)在無量綱直角坐標(biāo)系Ox?y?z?中的坐標(biāo)為(0,0,-1).令A(yù)點(diǎn)為肥皂泡上任意一點(diǎn),存在唯一直線經(jīng)過A點(diǎn)與B點(diǎn),直線AB與赤道平面的交點(diǎn)為C,C點(diǎn)就是A點(diǎn)的赤平極射投影點(diǎn).圖 1 展示了赤平極射投影幾何關(guān)系,其中橙色平面為赤道平面,A是肥皂泡上任意一點(diǎn),B為肥皂泡頂點(diǎn)關(guān)于赤道平面的鏡像點(diǎn),C為投影點(diǎn).投影坐標(biāo)系Ox?y?的原點(diǎn)為肥皂泡的球心,x?軸與x?軸重合,y?軸與y?軸重合.由簡單的幾何關(guān)系可以推導(dǎo)出,無量綱的直角坐標(biāo)系向投影空間的正變換is為

    圖1 赤平極射投影法示意圖Fig.1 An illustration of the stereographic projection

    而對應(yīng)的逆變換為

    在計(jì)算空間內(nèi),基矢不再是常向量,而隨著空間坐標(biāo)變化而變化.通過基矢的變換(12)與逆變換(13)可以方便地求得在投影空間與物理空間的矢量之間的轉(zhuǎn)換公式.例如,速度在投影空間與物理空間的分量具有以下關(guān)系:

    使用基矢的逆變換(13)代入到速度定義(14)中,即可獲得速度的逆變換公式:

    在投影空間中求解控制方程(4)還需要計(jì)算矢量算符在投影空間中的形式.變換is的Jacobi 矩陣為

    2 數(shù)值方法

    圖2 交錯(cuò)網(wǎng)格示意圖Fig.2 The scheme of the staggered mesh

    圖3 并行算法示意圖Fig.3 The scheme of the parallel algorithm

    3 結(jié)果討論

    3.1 算例信息與瞬時(shí)流場

    本文一共設(shè)置了三個(gè)算例用以分析Rayleigh 數(shù)對流場的影響,三個(gè)算例的Ra分別為3×107,3×109和3×1011,詳細(xì)信息見表 1.所有算例的Pr固定為7,因?yàn)榉试硭奈锢硇再|(zhì)與水近似.以往的研究表明S? =F?=0.06是一個(gè)合理的選擇[19,21],本文計(jì)算中也采用該值.為了消除網(wǎng)格分辨率對結(jié)果的影響,對每個(gè)算例都嘗試了多個(gè)網(wǎng)格分辨率:256×256,512×512,1 024×1 024,1 536×1 536,2 048×2 048和2 560×2 560.圖 4 顯示了算例2 在不同網(wǎng)格分辨率下,系統(tǒng)總擬渦能隨時(shí)間的變化曲線.從圖中可以觀察到,不同的網(wǎng)格分辨率下肥皂泡上的總擬渦能都在 50個(gè)無量綱時(shí)間單位內(nèi)達(dá)到了統(tǒng)計(jì)穩(wěn)態(tài).在本文所有的數(shù)值模擬中,每個(gè)算例都保證足夠的計(jì)算時(shí)長,使流場達(dá)到統(tǒng)計(jì)穩(wěn)態(tài),并對結(jié)果進(jìn)行長時(shí)間的統(tǒng)計(jì)計(jì)算.隨著網(wǎng)格分辨率的提升,總擬渦能的統(tǒng)計(jì)平均值在緩慢增加,但當(dāng)網(wǎng)格分辨率超過2 048×2 048時(shí),DNS 已可以獲得網(wǎng)格獨(dú)立結(jié)果.得益于并行計(jì)算效率的提升,本文首次使用2 048×2 048的網(wǎng)格分辨率開展了肥皂泡上熱對流的DNS,同時(shí)也將DNS 的Rayleigh 數(shù)提升到Ra=3×1011.

    圖4 不同網(wǎng)格分辨率下,算例2 的總擬渦能隨時(shí)間的變化曲線(左)和擬渦能總量的統(tǒng)計(jì)值隨網(wǎng)格分辨率變化曲線(右)Fig.4 The temporal evolution of the total enstrophy on the bubble with different mesh resolutions(left)&the variation of the mean total enstropy with the mesh resolutions(right)

    表1 算例參數(shù)信息Table 1 Information for the simulated cases

    通過觀察瞬態(tài)流場云圖可以直觀地分析肥皂泡上流場的宏觀特點(diǎn).圖 5 展示了三個(gè)算例的瞬時(shí)溫度場T,瞬時(shí)動(dòng)能場u2/2和瞬時(shí)擬渦能場ω2/2.從圖 5 中可以清晰地觀察肥皂泡上非常卷曲的羽流以及尺度較大的島渦.在經(jīng)典的Rayleigh-Bénard 對流中,流場由冷熱羽流所自組織的大尺度環(huán)流(large scale circulation)所主導(dǎo).當(dāng)對流槽具有特定的縱橫比時(shí),在合適的Rayleigh 數(shù)控制下,流場中還會(huì)出現(xiàn)角渦.肥皂泡沒有冷邊界,流場中不存在冷羽流;而肥皂泡同時(shí)也沒有側(cè)向邊界及其形成的角渦,肥皂泡上的相干結(jié)構(gòu)主要由熱羽流及其演化形成.從溫度云圖上可以觀察到,邊界層上產(chǎn)生了大量的羽流,羽流在浮力的作用下向更高緯度運(yùn)動(dòng).而同一時(shí)刻的擬渦能云圖則顯示羽流的形成與運(yùn)動(dòng)過程中產(chǎn)生了不同尺寸的渦,一些渦相互融合后形成了尺度較大的島渦,島渦主要居于較高緯度,與實(shí)驗(yàn)中觀察到的情況相同[15,17].從動(dòng)能云圖中可以發(fā)現(xiàn)流場的動(dòng)能集中在一些尺寸較大的島渦中.

    圖5 瞬時(shí)流場圖:(a)無量綱瞬時(shí)溫度場,Ra=3×107;(b)無量綱瞬時(shí)溫度場,Ra=3×109;(c)無量綱瞬時(shí)溫度場,Ra=3×1011;(d)無量綱瞬時(shí)動(dòng)能場,Ra=3×107;(e)無量綱瞬時(shí)動(dòng)能場,Ra=3×109;(f)無量綱瞬時(shí)動(dòng)能場,Ra=3×1011;(g)無量綱瞬時(shí)擬渦能場,Ra=3×107;(h)無量綱瞬時(shí)擬渦能場,Ra=3×109;(i)無量綱瞬時(shí)擬渦能場,Ra=3×1011Fig.5 The instantaneous flow field:(a)dimensionless instantaneous temperature field,Ra=3×107;(b)dimensionless instantaneous temperature field,Ra=3×109;(c)dimensionless instantaneous temperature field,Ra=3×1011;(d)dimensionless instantaneous kinetic energy field,Ra=3×107;(e)dimensionless instantaneous kinetic energy field,Ra=3×109;(f)dimensionless instantaneous kinetic energy field,Ra=3×1011;(g)dimensionless instantaneous enstrophy field,Ra=3×107;(h)dimensionless instantaneous enstrophy field,Ra=3×109;(i)dimensionless instantaneous enstrophy field,Ra=3×1011

    而隨著Ra的增加,瞬時(shí)流場出現(xiàn)了明顯變化.當(dāng)Ra=3×107時(shí),羽流和渦的尺寸較大,動(dòng)能分布在較大的流動(dòng)結(jié)構(gòu)中,而且可以從溫度云圖中清晰的觀察到邊界層.當(dāng)Ra增加到3×109時(shí),羽流更加細(xì)小與卷曲,而渦的尺寸也明顯變小,同時(shí)羽流和渦的數(shù)量有明顯增加,動(dòng)能集中在更加細(xì)長的流動(dòng)結(jié)構(gòu)中.另外,溫度邊界層的厚度大幅減小,從溫度云圖中已經(jīng)無法清晰觀察到邊界層.而當(dāng)Ra=3×1011時(shí),羽流和渦變得更加小,而其數(shù)量更進(jìn)一步增加時(shí),動(dòng)能云圖顯示流動(dòng)結(jié)構(gòu)的尺寸將繼續(xù)減小.

    3.2 動(dòng)能與擬熱能的波數(shù)譜與通量

    流場量的波數(shù)譜和通量表征了湍流場不同尺度的渦之間的相互作用,對于三維流場或二維平面流場,計(jì)算波數(shù)譜和通量需要進(jìn)行空間Fourier 變換.肥皂泡上的流場具有曲面幾何,需要專門的數(shù)學(xué)工具來計(jì)算波數(shù)譜和通量.任意球?qū)ΨQ空間函數(shù)Ψ 均可以展開為球諧函數(shù)的級(jí)數(shù):

    在肥皂泡上尋找Bo59標(biāo)度律需要考慮曲面幾何帶來的影響.由于肥皂泡上不同緯度上的重力投影隨著緯度而變化,這是Bo59假設(shè)中所不包含的特點(diǎn).這暗示只有在重力隨緯度變化較小的區(qū)域內(nèi),有可能發(fā)現(xiàn)Bo59標(biāo)度律.對于緯度為 θ*,重力加速度在球面切向的投影為gcosθ*,對之求 θ*的導(dǎo)數(shù)為gsinθ*,這說明越接近赤道重力投影隨緯度的變化越小.然而,在導(dǎo)出Bo59標(biāo)度律的理論中假定邊界條件不產(chǎn)生影響,因而只有在足夠遠(yuǎn)離赤道才能發(fā)現(xiàn)Bo59標(biāo)度律.結(jié)合以上幾點(diǎn)考慮,Bo59標(biāo)度律只在一定緯度范圍內(nèi)有效.

    本文三個(gè)算例的動(dòng)能與擬熱能的波數(shù)譜如圖 6 所示,圖中虛線為Bo59標(biāo)度參考線.觀察圖 6 中的動(dòng)能波數(shù)譜可以發(fā)現(xiàn),在5≤kR≤100的范圍內(nèi),Bo59標(biāo)度律準(zhǔn)確地描述了不同Ra下的動(dòng)能譜隨波數(shù)的標(biāo)度規(guī)律.相反的是,圖 6 中的擬熱能的波數(shù)譜則表明,在本文所考慮的三個(gè)Ra條件下,肥皂泡的擬熱能波數(shù)譜標(biāo)度特征沒有觀察到Bo59標(biāo)度律.在肥皂泡上流場中,邊界層中流體的溫度最高,而當(dāng)流體隨著羽流上升離開邊界層后,溫度快速下降.He 等[21]發(fā)現(xiàn),邊界層的平均溫度與脈動(dòng)溫度均遠(yuǎn)大于肥皂泡其余部分,這暗示邊界層對肥皂泡上溫度場具有很強(qiáng)的影響.圖 6 中的擬熱能波數(shù)譜包含了整個(gè)肥皂泡上溫度場的脈動(dòng)信息,溫度邊界條件對波數(shù)譜的影響無法忽略,所以擬熱能波數(shù)譜相比動(dòng)能包含了更多邊界的影響而沒有Bo59標(biāo)度律.隨著Ra的增加,在小波數(shù)的范圍內(nèi),動(dòng)能譜和擬熱能有較明顯的減小,這與瞬時(shí)流場中羽流和渦的尺寸隨著Ra增加而明顯變小的規(guī)律相符合.對于圖 6 中所有Ra的波數(shù)譜,隨著波數(shù)k的增加,動(dòng)能譜與擬熱能譜逐漸減小.這表明肥皂泡上尺寸越大的流動(dòng)結(jié)構(gòu)具有更大的動(dòng)能與擬熱能,與從瞬時(shí)流場觀察到的現(xiàn)象相同.

    圖6 動(dòng)能與擬熱能的波數(shù)譜:(a)動(dòng)能;(b)擬熱能Fig.6 The wave number spectra for the kinetic energy and entropy:(a)the kinetic energy;(b)the entropy

    通量(flux)代表了湍流場中不同尺度的渦之間的物理量傳遞,大于0的通量表明物理量從大尺度的渦向小尺度的渦級(jí)串.反之當(dāng)通量小于 0時(shí),物理量從小尺度的渦向大尺度渦逆級(jí)串.Kraichnan[7]、Leith[8]和Batchelor[9]的二維湍流理論預(yù)言了雙級(jí)串的現(xiàn)象:一個(gè)尺度范圍中,動(dòng)能能級(jí)逆級(jí)串而擬渦能保持正向級(jí)串.圖 7 顯示了整個(gè)肥皂泡上動(dòng)能、擬渦能、擬熱能和浮力項(xiàng)的通量,圖中包含全部三個(gè)算例的結(jié)果.在三個(gè)不同的Ra條件下,圖 7 中的 ΠE曲線都具有相同的規(guī)律,在較小波數(shù)的區(qū)間ΠE>0,隨著波數(shù)的增加,ΠE快速減小到負(fù)值.在 ΠE達(dá)到最小值之后,ΠE隨著波數(shù)的增加逐漸增加,重新達(dá)到ΠE>0,此時(shí) ΠE隨波數(shù)的變化較小.三個(gè)算例中的流場均出現(xiàn)了動(dòng)能能級(jí)逆級(jí)串的現(xiàn)象,隨著Ra的增加,能級(jí)擬級(jí)串的尺度區(qū)間逐漸變大.

    通過觀察圖 7 中的擬渦能通量曲線可以發(fā)現(xiàn),三個(gè)算例的 ΠW始終大于0,肥皂泡上擬渦能始終保持從大渦向小渦的傳遞方向.在不同的Ra下,ΠW在較小波數(shù)的范圍內(nèi)保持穩(wěn)定,但隨著波數(shù)的增加,ΠW先快速增加,達(dá)到最大值之后又快速減少.三個(gè)算例的擬熱能的通量同樣始終大于0,這說明溫度也始終是由較大的渦傳遞給較小的渦.從瞬時(shí)溫度場可以觀察到,羽流和渦在赤道處產(chǎn)生時(shí)尺寸較大,在向高緯度運(yùn)動(dòng)的過程中破碎消散,這一過程符合擬熱能通量所描述的逆熱能傳遞特征.當(dāng)Ra=3×1011與Ra=3×109時(shí),ΠT隨著波數(shù)的增加單調(diào)緩慢減小.而在Ra=3×107的條件下,擬渦能通量曲線隨著波數(shù)的增加先增加后減小,出現(xiàn)一個(gè)極大值.

    浮力通量是對整個(gè)肥皂泡上的浮力場Teθ*進(jìn)行球諧分解后計(jì)算獲得,計(jì)算方法與動(dòng)能的通量計(jì)算方法相同.由于肥皂泡上的流體運(yùn)動(dòng)是由浮力所驅(qū)動(dòng),浮力通量具有特別的重要性.圖 7 中的浮力通量曲線顯示三個(gè)算例的ΠBuoy始終大于0.不同Ra下的浮力通量曲線具有相同的變化規(guī)律.在區(qū)間kR≤5中,ΠBuoy隨著波數(shù)增加保持穩(wěn)定.在kR≥10的條件下,浮力通量隨著波數(shù)的增加而增加,達(dá)到最大值之后又快速減小.ΠBuoy的上升暗示一個(gè)在較小尺度上存在局部的湍流動(dòng)能注入,而這個(gè)湍流動(dòng)能注入來源是動(dòng)能能量逆級(jí)串:小尺度渦向大尺度渦注入了湍流動(dòng)能.隨著Ra的增加,ΠBuoy顯著減小.由于Ra的增加會(huì)使得羽流的尺寸快速減小,而且小尺寸的熱羽流受到的浮力更小,于是浮力通量在整個(gè)波數(shù)范圍上都減小了.

    圖7 肥皂泡上的動(dòng)能通量、擬渦能通量、擬熱能通量和浮力的通量:(a)動(dòng)能;(b)擬渦能;(c)擬熱能;(d)浮力Fig.7 The fluxes of kinetic energy,enstrophy,entropy and buoyancy on the soap bubble:(a)the kinetic energy;(b)the enstrophy;(c)the entropy;(d)the buoyancy

    3.3 結(jié)構(gòu)函數(shù)

    Kolmogrov 在1941 年提出各向同性均勻湍流假設(shè),并使用結(jié)構(gòu)函數(shù)來描述不同尺度渦之間的能量交換.結(jié)構(gòu)函數(shù)被定義在物理空間中,更加直觀地揭示出流場脈動(dòng)的特點(diǎn),并且能夠給出更高階的流場信息.n階速度結(jié)構(gòu)函數(shù)的定義如下:

    x是空間中一點(diǎn)的位置矢量,而d為距離矢量,d=|d|,n為正整數(shù),〈·〉算符代表系綜平均.按照同樣的方法,可以定義n階溫度結(jié)構(gòu)函數(shù):

    圖8 結(jié)構(gòu)函數(shù)計(jì)算方法Fig.8 The scheme of the distance for calculating the structure function

    由于溫度場與速度場的具有連續(xù)性,當(dāng)d足夠小時(shí),溫度與速度結(jié)構(gòu)函數(shù)的標(biāo)度系數(shù)退化為n以滿足Taylor 展開特性.在計(jì)算結(jié)構(gòu)函數(shù)時(shí),需要對物理量的差值進(jìn)行時(shí)間,經(jīng)度和緯度上的平均.所有達(dá)到統(tǒng)計(jì)平衡態(tài)之后的時(shí)間步都在時(shí)間平均的范圍內(nèi),而經(jīng)度平均的范圍在0?≤?*≤360?.為了能在消除邊界層的影響同時(shí)保證計(jì)算的準(zhǔn)確性,緯度平均的范圍被限制在10?≤θ*≤90?.圖 9 顯示了緯度與經(jīng)度方向上的從2 到9 階結(jié)構(gòu)函數(shù),其中緯度方向上的結(jié)構(gòu)函數(shù)在左邊,經(jīng)度方向上的結(jié)構(gòu)函數(shù)在右邊.將緯度方向上的溫度或速度結(jié)構(gòu)函數(shù)與經(jīng)度方向上的結(jié)構(gòu)函數(shù)相對比,可以發(fā)現(xiàn)兩者的標(biāo)度規(guī)律非常接近,但是數(shù)值有一些微小的差別.這與Bo59標(biāo)度律成立的前提條件相符合:速度或溫度的差值僅與兩點(diǎn)之間的距離相關(guān),與方向無關(guān).

    圖9 緯度與經(jīng)度方向上2 到9 階的溫度和速度結(jié)構(gòu)函數(shù):(a)緯度方向的溫度結(jié)構(gòu)函數(shù);(b)緯度方向的速度結(jié)構(gòu)函數(shù);(c)經(jīng)度方向的溫度結(jié)構(gòu)函數(shù);(d)經(jīng)度方向的速度結(jié)構(gòu)函數(shù)Fig.9 The temperature and velocity structure functions in the latitude and longitude directions,n=2~9:(a)the temperature structure functions in the latitude direction;(b)the velocity structure functions in the latitude direction;(c)the temperature structure functions in the longitude direction;(d)the velocity structure functions in the longitude direction

    從圖 9 可以觀察到,速度結(jié)構(gòu)函數(shù)在d/R≤O(0.01)的范圍內(nèi),標(biāo)度系數(shù)為n.而在0.01≤d/R≤0.1的區(qū)間中,Bo59標(biāo)度律較好地描述了速度結(jié)構(gòu)函數(shù)的標(biāo)度特征.在階數(shù)n逐漸增加的情況下,速度結(jié)構(gòu)函數(shù)始終能在同一區(qū)間內(nèi)較好地滿足Bo59標(biāo)度律,但受到湍流間歇的影響,其指數(shù)在高階時(shí)并非線性增加.

    另一方面,觀察溫度結(jié)構(gòu)函數(shù)曲線上可以發(fā)現(xiàn),溫度結(jié)構(gòu)函數(shù)在d/R≤O(0.01)的區(qū)間內(nèi)滿足而在 0.01≤d/R≤0.1的區(qū)間內(nèi)較好地符合Bo59標(biāo)度律.由于排除了邊界層的影響,結(jié)構(gòu)函數(shù)得到的結(jié)果與波數(shù)譜不同.

    4 結(jié)論與展望

    本文介紹了底部加熱的半個(gè)肥皂泡這一新的準(zhǔn)二維湍流,同時(shí)詳細(xì)介紹了其DNS 方法.計(jì)算中針對肥皂泡二維球面的曲面幾何特點(diǎn),借助球極投影,將物理空間的控制方程轉(zhuǎn)換到投影坐標(biāo)系中,降低了離散化控制方程的難度.此后還介紹了球面上利用球諧分解計(jì)算波數(shù)譜與通量,以及計(jì)算結(jié)構(gòu)函數(shù)的分析方法.最后,針對Ra=3×107,Ra=3×109,Ra=3×1011三個(gè)算例,本文計(jì)算了動(dòng)能和擬熱能的波數(shù)譜,以及動(dòng)能、擬熱能、擬渦能和浮力的通量.另外也計(jì)算了Ra=3×1011算例的溫度與速度結(jié)構(gòu)函數(shù).不同Ra的結(jié)果都表明了肥皂泡上存在二維湍流中雙級(jí)串的現(xiàn)象.湍流動(dòng)能波數(shù)譜較好地滿足了Bo59標(biāo)度律的理論預(yù)示,但擬熱能由于主要受邊界的影響較大,其波數(shù)譜中沒有觀察到足Bo59標(biāo)度律.但通過湍流結(jié)構(gòu)函數(shù),在去除了赤道邊界的高溫邊界層影響后,無論速度還是溫度結(jié)構(gòu)函數(shù),都在一定范圍較好地符合Bo59的理論預(yù)示.未來可進(jìn)一步通過增加Rayleigh 數(shù)范圍,研究Rayleigh 數(shù)與能量注入尺度以及不同維度上浮力作用間的聯(lián)系.

    99久久无色码亚洲精品果冻| 男女视频在线观看网站免费 | 国产69精品久久久久777片 | 精品人妻1区二区| 哪里可以看免费的av片| 欧美一区二区精品小视频在线| 成人国产综合亚洲| 91大片在线观看| 99精品久久久久人妻精品| 99精品欧美一区二区三区四区| 999久久久精品免费观看国产| 人妻丰满熟妇av一区二区三区| 免费在线观看黄色视频的| 国产成人精品久久二区二区91| 欧美黄色片欧美黄色片| 天天躁狠狠躁夜夜躁狠狠躁| 国产一区二区三区视频了| 亚洲av片天天在线观看| 又爽又黄无遮挡网站| a级毛片在线看网站| 亚洲无线在线观看| 国产av一区二区精品久久| 久久久久久九九精品二区国产 | av在线播放免费不卡| 欧美一区二区精品小视频在线| 国产成+人综合+亚洲专区| 欧美高清成人免费视频www| 中文在线观看免费www的网站 | 一二三四在线观看免费中文在| 欧美高清成人免费视频www| 国产又黄又爽又无遮挡在线| 国产成人欧美在线观看| 国内精品久久久久久久电影| 青草久久国产| 少妇裸体淫交视频免费看高清 | 麻豆国产97在线/欧美 | 老汉色av国产亚洲站长工具| 精品人妻1区二区| 日韩精品中文字幕看吧| 一级毛片高清免费大全| 精品久久蜜臀av无| 欧美人与性动交α欧美精品济南到| 色播亚洲综合网| 国产黄a三级三级三级人| 他把我摸到了高潮在线观看| 亚洲精品一区av在线观看| 12—13女人毛片做爰片一| 亚洲va日本ⅴa欧美va伊人久久| 两个人看的免费小视频| 亚洲免费av在线视频| tocl精华| 最近最新免费中文字幕在线| 波多野结衣高清作品| 男女下面进入的视频免费午夜| 成人国产一区最新在线观看| 国产亚洲精品久久久久久毛片| 精品高清国产在线一区| 50天的宝宝边吃奶边哭怎么回事| 国产片内射在线| 99久久综合精品五月天人人| 一级作爱视频免费观看| 69av精品久久久久久| www.999成人在线观看| 国产在线观看jvid| 悠悠久久av| 色老头精品视频在线观看| 日韩大尺度精品在线看网址| 免费在线观看亚洲国产| 久久久国产成人免费| 少妇人妻一区二区三区视频| 久久国产乱子伦精品免费另类| 国产探花在线观看一区二区| 精品国产乱子伦一区二区三区| 毛片女人毛片| 日日干狠狠操夜夜爽| 精品国产乱码久久久久久男人| 免费高清视频大片| 亚洲精品国产精品久久久不卡| 在线观看一区二区三区| 免费在线观看影片大全网站| 亚洲自偷自拍图片 自拍| 国产激情偷乱视频一区二区| 精品一区二区三区av网在线观看| 在线播放国产精品三级| 校园春色视频在线观看| 法律面前人人平等表现在哪些方面| 狂野欧美白嫩少妇大欣赏| 岛国在线免费视频观看| 午夜免费观看网址| 午夜免费观看网址| 国产精品久久久人人做人人爽| 久久久国产欧美日韩av| 国产一级毛片七仙女欲春2| 国语自产精品视频在线第100页| av福利片在线观看| 欧美激情久久久久久爽电影| 天天一区二区日本电影三级| 精品一区二区三区av网在线观看| 国内久久婷婷六月综合欲色啪| √禁漫天堂资源中文www| 国产精品永久免费网站| 老司机福利观看| 老司机福利观看| 一卡2卡三卡四卡精品乱码亚洲| 两人在一起打扑克的视频| 一本一本综合久久| 国产精华一区二区三区| 日本 欧美在线| 69av精品久久久久久| 一边摸一边抽搐一进一小说| 久久天躁狠狠躁夜夜2o2o| 一夜夜www| 国产精品久久久久久人妻精品电影| a级毛片a级免费在线| 亚洲国产看品久久| av在线播放免费不卡| 婷婷六月久久综合丁香| 丰满人妻一区二区三区视频av | 床上黄色一级片| 丰满人妻一区二区三区视频av | av视频在线观看入口| 亚洲人成77777在线视频| 麻豆国产97在线/欧美 | 日韩大码丰满熟妇| 免费搜索国产男女视频| 叶爱在线成人免费视频播放| 一夜夜www| 两个人视频免费观看高清| 欧美日韩乱码在线| 免费高清视频大片| 丰满人妻一区二区三区视频av | 日韩高清综合在线| 中亚洲国语对白在线视频| 亚洲黑人精品在线| 国产精品久久视频播放| 久9热在线精品视频| 亚洲成a人片在线一区二区| 欧美另类亚洲清纯唯美| 亚洲第一电影网av| 看片在线看免费视频| 亚洲av片天天在线观看| 欧美性长视频在线观看| 欧美成人午夜精品| 97碰自拍视频| 国产成人系列免费观看| 国产亚洲精品第一综合不卡| 天天躁狠狠躁夜夜躁狠狠躁| 久久天堂一区二区三区四区| 国产黄a三级三级三级人| 日韩欧美免费精品| 亚洲国产精品成人综合色| 久久久久九九精品影院| 久久午夜亚洲精品久久| 在线观看日韩欧美| www.精华液| 高清在线国产一区| 一级毛片女人18水好多| 亚洲自偷自拍图片 自拍| 国产精品国产高清国产av| 国产精品1区2区在线观看.| 欧美高清成人免费视频www| 给我免费播放毛片高清在线观看| 老鸭窝网址在线观看| 国产精品久久久人人做人人爽| www.自偷自拍.com| 午夜福利在线在线| 久久精品人妻少妇| 国语自产精品视频在线第100页| 又黄又粗又硬又大视频| 国产伦在线观看视频一区| 国产亚洲精品一区二区www| 国产亚洲精品av在线| 国产激情欧美一区二区| 最近最新中文字幕大全电影3| 久久久久久大精品| 国产一区二区在线av高清观看| 国产午夜精品久久久久久| 亚洲av美国av| 亚洲精品美女久久久久99蜜臀| 日韩三级视频一区二区三区| 美女黄网站色视频| 啦啦啦免费观看视频1| 国产男靠女视频免费网站| 免费在线观看影片大全网站| 国产亚洲精品久久久久久毛片| 免费在线观看视频国产中文字幕亚洲| 午夜精品一区二区三区免费看| 天堂影院成人在线观看| 成人亚洲精品av一区二区| 99热6这里只有精品| 久久久国产成人免费| 免费看十八禁软件| 久久天躁狠狠躁夜夜2o2o| 99久久无色码亚洲精品果冻| cao死你这个sao货| 黄频高清免费视频| АⅤ资源中文在线天堂| 免费在线观看完整版高清| 在线观看一区二区三区| 波多野结衣高清无吗| 久久久久久九九精品二区国产 | 国内少妇人妻偷人精品xxx网站 | 日本免费a在线| 免费在线观看完整版高清| 国产精品久久久人人做人人爽| 国产蜜桃级精品一区二区三区| 免费在线观看完整版高清| 免费看a级黄色片| 2021天堂中文幕一二区在线观| 国产精品久久久久久亚洲av鲁大| www.自偷自拍.com| 亚洲电影在线观看av| 国产精品1区2区在线观看.| 激情在线观看视频在线高清| 亚洲国产精品合色在线| 看黄色毛片网站| 色综合婷婷激情| 欧美最黄视频在线播放免费| 亚洲天堂国产精品一区在线| 日韩欧美国产在线观看| 亚洲在线自拍视频| 精品国产乱子伦一区二区三区| 午夜精品一区二区三区免费看| 久久久精品大字幕| 成人一区二区视频在线观看| 可以在线观看毛片的网站| 国产视频内射| 欧美日韩乱码在线| 国产精品一及| 波多野结衣高清作品| 88av欧美| videosex国产| 一个人免费在线观看电影 | 亚洲最大成人中文| 曰老女人黄片| 亚洲全国av大片| 国产高清激情床上av| 97人妻精品一区二区三区麻豆| 国产精品久久久av美女十八| www日本黄色视频网| 久久人人精品亚洲av| 国产野战对白在线观看| 日韩成人在线观看一区二区三区| 黄色 视频免费看| 在线观看一区二区三区| 国产私拍福利视频在线观看| 精品乱码久久久久久99久播| 国产成人精品久久二区二区91| 久久久久国内视频| www国产在线视频色| 国产av麻豆久久久久久久| 婷婷精品国产亚洲av| 亚洲成人精品中文字幕电影| 久久精品综合一区二区三区| 午夜精品在线福利| 久久精品成人免费网站| 久久这里只有精品19| 久久久久久久久免费视频了| 欧美日韩精品网址| 国产在线观看jvid| 久久中文字幕一级| 99精品久久久久人妻精品| 99热这里只有是精品50| 国产亚洲精品久久久久久毛片| 亚洲成人国产一区在线观看| 亚洲欧美一区二区三区黑人| 两人在一起打扑克的视频| 看片在线看免费视频| 亚洲精品在线美女| 亚洲精品久久成人aⅴ小说| 免费人成视频x8x8入口观看| 亚洲精品国产一区二区精华液| 亚洲一区二区三区不卡视频| 天天躁夜夜躁狠狠躁躁| 久久精品影院6| 国产精华一区二区三区| 亚洲第一电影网av| 小说图片视频综合网站| 久久这里只有精品中国| 在线观看美女被高潮喷水网站 | 宅男免费午夜| 亚洲国产欧美一区二区综合| 欧美日韩福利视频一区二区| 日韩欧美国产在线观看| 成人永久免费在线观看视频| 久久久久久九九精品二区国产 | 国产精品久久久av美女十八| 久久久久久久午夜电影| 成人午夜高清在线视频| 毛片女人毛片| 法律面前人人平等表现在哪些方面| 亚洲 欧美一区二区三区| 男男h啪啪无遮挡| 国产伦在线观看视频一区| 日韩中文字幕欧美一区二区| 女人高潮潮喷娇喘18禁视频| 日韩成人在线观看一区二区三区| 国产欧美日韩一区二区三| 99re在线观看精品视频| 亚洲精品美女久久久久99蜜臀| 日本 欧美在线| 99热6这里只有精品| 国产aⅴ精品一区二区三区波| 看片在线看免费视频| 老熟妇乱子伦视频在线观看| 人妻夜夜爽99麻豆av| 18禁黄网站禁片免费观看直播| 香蕉av资源在线| 91成年电影在线观看| 欧美日韩乱码在线| 亚洲一卡2卡3卡4卡5卡精品中文| 99热这里只有是精品50| 国产精品亚洲美女久久久| 一级黄色大片毛片| 伊人久久大香线蕉亚洲五| 国产又色又爽无遮挡免费看| 久久精品影院6| 亚洲国产精品成人综合色| 香蕉国产在线看| 成熟少妇高潮喷水视频| 午夜激情福利司机影院| www.精华液| 欧美极品一区二区三区四区| 中文字幕最新亚洲高清| 亚洲av成人一区二区三| 国产黄片美女视频| 少妇粗大呻吟视频| 99re在线观看精品视频| 中文字幕人妻丝袜一区二区| 国产精品久久视频播放| 正在播放国产对白刺激| 国产亚洲精品一区二区www| 三级男女做爰猛烈吃奶摸视频| 少妇粗大呻吟视频| 国产av一区二区精品久久| 免费高清视频大片| 大型av网站在线播放| 久久久久国产一级毛片高清牌| 免费在线观看日本一区| 国产精品免费视频内射| 伊人久久大香线蕉亚洲五| 欧美午夜高清在线| 久久天躁狠狠躁夜夜2o2o| 欧美久久黑人一区二区| av免费在线观看网站| 久久欧美精品欧美久久欧美| 岛国视频午夜一区免费看| 日韩精品免费视频一区二区三区| 淫秽高清视频在线观看| 嫁个100分男人电影在线观看| 精品日产1卡2卡| 精品高清国产在线一区| 成人三级黄色视频| av有码第一页| 黄片大片在线免费观看| 成年免费大片在线观看| 在线十欧美十亚洲十日本专区| 波多野结衣高清无吗| 18禁国产床啪视频网站| 一个人免费在线观看电影 | 日韩欧美精品v在线| 午夜亚洲福利在线播放| 国产日本99.免费观看| 人妻丰满熟妇av一区二区三区| 国语自产精品视频在线第100页| 少妇人妻一区二区三区视频| 床上黄色一级片| 精品午夜福利视频在线观看一区| 长腿黑丝高跟| 国产伦人伦偷精品视频| 热99re8久久精品国产| 国产伦一二天堂av在线观看| 男人舔女人的私密视频| av有码第一页| 精品一区二区三区四区五区乱码| 淫妇啪啪啪对白视频| 欧美性长视频在线观看| 国产野战对白在线观看| 好看av亚洲va欧美ⅴa在| 国产精品永久免费网站| 两性午夜刺激爽爽歪歪视频在线观看 | 人妻丰满熟妇av一区二区三区| 欧美乱码精品一区二区三区| 我的老师免费观看完整版| 男人舔女人下体高潮全视频| 久久久国产精品麻豆| 俺也久久电影网| 国产v大片淫在线免费观看| 国产探花在线观看一区二区| 在线a可以看的网站| 一级a爱片免费观看的视频| 97人妻精品一区二区三区麻豆| 久久中文看片网| 香蕉丝袜av| 午夜老司机福利片| 50天的宝宝边吃奶边哭怎么回事| 欧美不卡视频在线免费观看 | 日韩精品免费视频一区二区三区| 他把我摸到了高潮在线观看| 国产三级黄色录像| 91国产中文字幕| 精品久久久久久久久久免费视频| 成人永久免费在线观看视频| 久久婷婷成人综合色麻豆| 99热只有精品国产| 国产精品久久久久久人妻精品电影| 欧美黑人巨大hd| 母亲3免费完整高清在线观看| 可以免费在线观看a视频的电影网站| 99国产精品99久久久久| 超碰成人久久| 无遮挡黄片免费观看| 中文字幕熟女人妻在线| 亚洲成人久久爱视频| 亚洲精品av麻豆狂野| 国产午夜福利久久久久久| 日韩欧美三级三区| 18美女黄网站色大片免费观看| 免费在线观看成人毛片| 啦啦啦观看免费观看视频高清| 国产亚洲精品久久久久久毛片| 色综合亚洲欧美另类图片| 一级片免费观看大全| 色av中文字幕| 又黄又粗又硬又大视频| 此物有八面人人有两片| 国产一区二区在线av高清观看| 免费看日本二区| 亚洲色图 男人天堂 中文字幕| 成人永久免费在线观看视频| www日本在线高清视频| av天堂在线播放| 舔av片在线| 狂野欧美激情性xxxx| 两个人的视频大全免费| 五月玫瑰六月丁香| 欧美色视频一区免费| 久久国产精品人妻蜜桃| 在线a可以看的网站| 国产一区在线观看成人免费| 国产精品自产拍在线观看55亚洲| av在线天堂中文字幕| xxx96com| 国产爱豆传媒在线观看 | 黄色 视频免费看| 久久精品国产清高在天天线| 欧美在线一区亚洲| 久久人妻av系列| 好看av亚洲va欧美ⅴa在| 久久久国产欧美日韩av| 亚洲av第一区精品v没综合| 久99久视频精品免费| 色老头精品视频在线观看| 国产成人系列免费观看| 久久精品国产亚洲av高清一级| 久久精品国产99精品国产亚洲性色| 特级一级黄色大片| 亚洲av成人不卡在线观看播放网| 免费高清视频大片| 亚洲五月天丁香| 法律面前人人平等表现在哪些方面| 欧美日韩一级在线毛片| 亚洲欧美精品综合一区二区三区| 在线永久观看黄色视频| 久久这里只有精品19| 91在线观看av| 欧美av亚洲av综合av国产av| 精品国产乱码久久久久久男人| 夜夜躁狠狠躁天天躁| 18禁观看日本| 午夜亚洲福利在线播放| 伦理电影免费视频| 99re在线观看精品视频| 18美女黄网站色大片免费观看| 亚洲人成网站高清观看| 久久久久亚洲av毛片大全| 男插女下体视频免费在线播放| 久久 成人 亚洲| 伊人久久大香线蕉亚洲五| 嫩草影院精品99| bbb黄色大片| 中文字幕av在线有码专区| 18禁观看日本| 91国产中文字幕| 窝窝影院91人妻| 国产精品久久久久久亚洲av鲁大| 国产精品精品国产色婷婷| 久久久久久亚洲精品国产蜜桃av| 国产99白浆流出| 国产三级黄色录像| 欧美色欧美亚洲另类二区| 亚洲va日本ⅴa欧美va伊人久久| 国产欧美日韩精品亚洲av| e午夜精品久久久久久久| 悠悠久久av| 免费观看人在逋| 人妻久久中文字幕网| 老司机靠b影院| 亚洲18禁久久av| 精品熟女少妇八av免费久了| 婷婷精品国产亚洲av| 日本精品一区二区三区蜜桃| 99在线视频只有这里精品首页| av免费在线观看网站| 中文字幕人妻丝袜一区二区| 午夜激情福利司机影院| 一本大道久久a久久精品| 99re在线观看精品视频| av片东京热男人的天堂| av免费在线观看网站| 久久精品影院6| 精品久久蜜臀av无| 亚洲专区中文字幕在线| 51午夜福利影视在线观看| 精品久久久久久久久久免费视频| 91字幕亚洲| 亚洲一卡2卡3卡4卡5卡精品中文| 香蕉久久夜色| 男女下面进入的视频免费午夜| 操出白浆在线播放| 特级一级黄色大片| 国产精品,欧美在线| 日本撒尿小便嘘嘘汇集6| 国产精品98久久久久久宅男小说| 国产成年人精品一区二区| 国产av一区在线观看免费| 少妇粗大呻吟视频| 男女之事视频高清在线观看| 午夜久久久久精精品| 亚洲国产高清在线一区二区三| 国产精品自产拍在线观看55亚洲| 欧美中文综合在线视频| 久久久久国产精品人妻aⅴ院| 最近视频中文字幕2019在线8| 亚洲国产精品久久男人天堂| 婷婷精品国产亚洲av在线| 成人av在线播放网站| 又黄又粗又硬又大视频| 亚洲精品美女久久av网站| 99热这里只有是精品50| 天天躁狠狠躁夜夜躁狠狠躁| 午夜福利视频1000在线观看| 精品少妇一区二区三区视频日本电影| 欧美一区二区国产精品久久精品 | 久久久久精品国产欧美久久久| 日韩中文字幕欧美一区二区| 久久久国产精品麻豆| 熟女少妇亚洲综合色aaa.| 男女午夜视频在线观看| 欧美成人一区二区免费高清观看 | 母亲3免费完整高清在线观看| 亚洲人成网站高清观看| 女人高潮潮喷娇喘18禁视频| 最近最新免费中文字幕在线| 一区二区三区国产精品乱码| 亚洲精品国产精品久久久不卡| 美女扒开内裤让男人捅视频| 老司机靠b影院| 久久亚洲真实| 国产精品99久久99久久久不卡| 国产亚洲欧美在线一区二区| 校园春色视频在线观看| 日韩中文字幕欧美一区二区| 久久精品成人免费网站| 欧美高清成人免费视频www| 久久天堂一区二区三区四区| av福利片在线| 一进一出抽搐gif免费好疼| 国产成人精品无人区| 成年女人毛片免费观看观看9| 欧美日韩黄片免| 国产主播在线观看一区二区| 久9热在线精品视频| 久久婷婷成人综合色麻豆| 制服丝袜大香蕉在线| 日日爽夜夜爽网站| 国产三级黄色录像| 国产精品 欧美亚洲| 免费看a级黄色片| 成人手机av| 999久久久精品免费观看国产| 舔av片在线| 亚洲av成人精品一区久久| 美女高潮喷水抽搐中文字幕| 一卡2卡三卡四卡精品乱码亚洲| 老汉色∧v一级毛片| a级毛片在线看网站| 国产亚洲精品久久久久久毛片| 香蕉丝袜av| 脱女人内裤的视频| av有码第一页| 久久久久精品国产欧美久久久| av中文乱码字幕在线| 日韩高清综合在线| 熟女电影av网| 国产高清有码在线观看视频 | 三级男女做爰猛烈吃奶摸视频| 日韩三级视频一区二区三区| 熟女少妇亚洲综合色aaa.| 久久精品夜夜夜夜夜久久蜜豆 | 免费看美女性在线毛片视频| 国产伦一二天堂av在线观看| 琪琪午夜伦伦电影理论片6080| 国产精品一区二区免费欧美| 久久精品影院6| 色尼玛亚洲综合影院| 午夜福利18| 欧美成人免费av一区二区三区| 日本免费一区二区三区高清不卡| 男插女下体视频免费在线播放| 亚洲国产日韩欧美精品在线观看 | 日韩av在线大香蕉| 男女下面进入的视频免费午夜| www.熟女人妻精品国产| 日本五十路高清|