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

    旋轉(zhuǎn)肥皂泡熱對流能量耗散與邊界層特性的數(shù)值模擬*

    2022-10-27 02:59:12賀嘯秋熊永亮彭澤瑞徐順
    物理學(xué)報 2022年20期
    關(guān)鍵詞:熱對流肥皂泡赤道

    賀嘯秋 熊永亮? 彭澤瑞 徐順

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

    2) (工程結(jié)構(gòu)分析與安全評定湖北省重點實驗室,武漢 430074)

    將底部加熱的半個肥皂泡作為一個新的熱對流模型,結(jié)合了肥皂泡固有的球面與準(zhǔn)二維特征,由此有助于理解行星大氣流動中的復(fù)雜物理機制與熱對流特性.本文使用直接數(shù)值模擬方法計算了旋轉(zhuǎn)肥皂泡上的湍流熱對流,研究了肥皂泡上的溫度與黏性邊界層以及擬熱能和動能耗散規(guī)律.結(jié)合肥皂泡上溫度場與速度場特征,分別根據(jù)溫度脈動均方根最大值以及速度脈動邊界處斜率延長線與最大值交點提出了肥皂泡上溫度與黏性邊界層的識別方法.研究發(fā)現(xiàn),當(dāng)肥皂泡從邊界吸收能量時,擬熱能耗散與動能耗散均集中在邊界層中,肥皂泡上的溫度邊界層與黏性邊界層厚度與瑞利數(shù) Ra 存在明確的標(biāo)度關(guān)系.相比經(jīng)典Rayleigh-Bénrad對流(RB 對流)模型,溫度標(biāo)度指數(shù)具有較為接近的結(jié)果,但速度標(biāo)度指數(shù)存在一定的差異.此外,在混合區(qū),均方根溫度(T*)隨緯度(θ)具有近似 T*~θ0.5 的標(biāo)度關(guān)系,這與RB 對流模型及其相應(yīng)的理論預(yù)測一致.最后通過能量平衡方程發(fā)現(xiàn),肥皂泡上擬熱能內(nèi)耗散率 和動能內(nèi)耗散率 比擬熱能外耗散率和動能外耗散率 大1 個量級,擬熱能與動能的內(nèi)部耗散率在邊界層中具有支配地位,隨著肥皂泡旋轉(zhuǎn)速率的增加,熱羽流難以輸運到高緯度地區(qū),進一步降低了擬熱能與動能外耗散率的影響.

    1 引言

    海洋的平均深度約為數(shù)千米,但縱向長度可以達到數(shù)千公里;同樣,地球大氣圈的厚度只有數(shù)十公里,但是在經(jīng)度與緯度方向的尺寸達到了數(shù)千公里,因而大尺度的海洋和大氣流動具有一定的二維流動特征[1,2].先前的研究表明,二維湍流可以較好地近似描述宏觀大尺度的大氣與海洋流動規(guī)律[2].類似地,肥皂薄膜的宏觀尺寸可以達到厘米量級,但是其厚度只在微米量級,因此流體在肥皂薄膜中的湍流運動非常好地近似了二維湍流,與行星大氣與海洋的大尺度流動具有一定相似性[3,4],實驗中廣泛地將肥皂薄膜應(yīng)用于二維湍流的研究[5].Wu等[6]使用肥皂薄膜研究了二維Couette 流動,并驗證了被動標(biāo)量在二維湍流中的擴散所遵守的標(biāo)度率.Kellay等[7]則使用豎直肥皂薄膜驗證了二維湍流的能量逆級串規(guī)律.近十年來,曲面肥皂薄膜[8]也被應(yīng)用于實驗研究中: 底部加熱的肥皂泡作為一個較獨特的熱對流模型受到了很多學(xué)者的關(guān)注[3,4,9,10-14].Kellay[8]成功將半個肥皂泡長時間穩(wěn)定地固定在一個水浴加熱的銅制基座上.通過旋轉(zhuǎn)銅質(zhì)基座,可以實現(xiàn)肥皂泡繞對稱軸進行旋轉(zhuǎn)[13].肥皂泡被底座加熱,肥皂水在浮力的驅(qū)動下進行對流輸運;當(dāng)加熱足夠強時,肥皂泡上就產(chǎn)生了湍流熱對流[9].

    在實驗中,肥皂泡上出現(xiàn)了能長時間存在的島渦,學(xué)者發(fā)現(xiàn)島渦的運動軌跡與颶風(fēng)軌跡均滿足超擴散關(guān)系[9].實驗結(jié)果表明肥皂泡上島渦的運動軌跡特征以及強度與隨機觀察的一百多個颶風(fēng)的軌跡以及強度之間具有非常類似的特征[9,12],學(xué)者依據(jù)這些特征建立了快速預(yù)報颶風(fēng)軌跡錐的模型[3,12].肥皂泡模型刻畫了流體在二維曲面中的湍流熱對流運動,同時也有助于深入了解二維湍流熱對流的物理機制[4,14].現(xiàn)階段,對肥皂泡大尺度特征[3,9,12]和小尺度特征[4,10,13,14]都開展了一定研究.Bruneau等[4]使用直接數(shù)值模擬(DNS)方法研究了靜止肥皂泡上的湍流熱對流,發(fā)現(xiàn)了肥皂泡上能量串級符合Bolgiano-Oberbeck 的標(biāo)度規(guī)律.而He等[14]使用DNS 方法研究了旋轉(zhuǎn)對肥皂泡上流場與溫度場的影響,獲得了旋轉(zhuǎn)對肥皂泡上小尺度脈動的影響規(guī)律.

    在三維湍流系統(tǒng)中,能量由大尺度注入湍流,并在慣性區(qū)間內(nèi)發(fā)生能量級串,最終在Kolmogrov尺度以下被消耗殆盡[15,16].而在二維情況下,湍流中存在雙向能量級串: 湍流動能與擬渦能沿不同的方向傳遞[2].能量耗散對湍流場中的當(dāng)?shù)孛}動具有重要的影響.慣性區(qū)間內(nèi)的動能級串的過程中,二階速度結(jié)構(gòu)函數(shù)的標(biāo)度直接受到動能耗散率的影響[17,18],其中為速度場,為運動黏度.在湍流熱對流中,溫度是主動標(biāo)量,溫度脈動在流體中的傳播也具有級串形式[19],其標(biāo)度規(guī)律與擬熱能耗散率直接相關(guān)[20],其中為熱擴散系數(shù),而為溫度場.同時,能量耗散的高低對湍流熱對流系統(tǒng)的全局熱通量具有調(diào)節(jié)作用[21].于是,研究能量耗散對揭示湍流熱對流的物理機制具有重要的作用.例如,經(jīng)典的Rayleigh-Bénrad 熱對流(RB 對流)包含一層厚度為的流體和上下兩個水平的恒溫邊界: 下邊界溫度高于上邊界溫度,溫差為,熱能被流體從下邊界輸運至上邊界[22,23].RB 對流的控制參數(shù)為瑞利數(shù)Ra和普朗特數(shù)Pr:

    其中是重力加速度的大小,是熱膨脹系數(shù).人們最關(guān)心Ra與Pr如何影響RB 對流的熱輸運效率[22,23].熱輸運效率用努塞爾特數(shù)Nu表征,即通過這一層流體的無量綱熱通量:

    式 中是選取的特征速 度.Nu與和有以下準(zhǔn)確關(guān)系[21]:

    式中,〈·〉V代表對流體體積求空間與時間平均.GL理論基于(5)式和(6)式推導(dǎo)出了Nu與Re隨Ra和Pr變化的公式,并得到了大量研究的驗證[24-28].GL 理論認為,在湍流熱對流中能量耗散通常集中在邊界層.事實上,Grossmann 與Lohse[24-26]在描述RB 對流中的能量耗散時,將邊界層中的作用與湍流體中的作用分開討論.人們在進一步發(fā)展GL 理論的過程中,將羽流的影響加入了這一物理圖像,并將羽流歸納為邊界層流體脫離邊界而產(chǎn)生的運動[27,28].許多數(shù)值仿真[29,30]與實驗研究[31-33]都驗證了這一物理圖像的正確性.最近,人們又利用區(qū)分研究邊界層與湍流體中的平均能量耗散這一思想將GL 理論推廣至水平熱對流[34]、雙標(biāo)量RB 對流[35]、內(nèi)加熱對流[36]以及傾斜RB 對流[37].對于以上提及的多個熱對流模型,GL 理論的推廣理論成功得到了Nu和Re隨Ra和Pr變化的規(guī)律,并被多個相互獨立的實驗和數(shù)值研究所驗證[34-37].這說明了GL 理論具有較好的通用性,GL 理論所描述的物理圖像在湍流熱對流中廣泛存在.然而,據(jù)作者所知,肥皂泡模型的研究中尚未有專門和系統(tǒng)的工作投入于研究能量耗散與邊界層特性.所以本文希望借助GL 理論提出的物理圖像,對肥皂泡上熱對流的能量耗散與邊界層特性進行深入的研究,并分析肥皂泡旋轉(zhuǎn)所帶來的影響,從而對肥皂泡模型的物理特點做進一步的揭示.

    由于RB 對流簡單的邊界與幾何特征中蘊含了熱對流的本質(zhì)規(guī)律,目前對熱對流的研究主要集中在經(jīng)典的RB 對流模型,人們對RB 對流中能量的輸運與耗散做了深入細致的研究,取得了豐富的成果[22,23].在經(jīng)典RB 對流模型的研究中,旋轉(zhuǎn)所帶來的影響也是熱點方向之一,即所謂的旋轉(zhuǎn)Rayleigh-Bénrad 對流(RRB 對流).旋轉(zhuǎn)作用由羅斯比數(shù)的倒數(shù) 1/Ro表征,1/Ro越小,旋轉(zhuǎn)作用越弱;反之 1/Ro越大,旋轉(zhuǎn)作用越強.對于旋轉(zhuǎn)速度為的RRB 對流,羅斯比數(shù)Ro的定義為

    RRB 對流中,隨著 1/Ro的逐漸增加,流場依次進入弱旋轉(zhuǎn)態(tài)、中旋轉(zhuǎn)態(tài)與強旋轉(zhuǎn)態(tài)[38].在弱旋轉(zhuǎn)態(tài)中,流場中旋轉(zhuǎn)的作用基本可以忽略[39].而當(dāng)流場處于中旋轉(zhuǎn)態(tài)時,由于Ekman 泵吸效應(yīng),Nu隨著1/Ro增大而有一定程度的上升[40-43].在強旋轉(zhuǎn)態(tài)中,流體的運動被抑制,Nu隨著 1/Ro的增大而劇烈下降[44].在靠近上下水平邊界處,存在一層非常薄的溫度邊界層,RB 對流的主要熱阻源于溫度邊界層[22].而RB 對流中能量的輸運與耗散均與邊界層有重要的關(guān)系[29,30,33].在RB 對流中,大尺度環(huán)流(largescale circulation,LSC)和角渦是兩個獨特的流動特征,對Nu與Re有著重要的影響.黃茂靜和包蕓[45]利用DNS 方法研究了二維與三維方腔在不同Ra下溫度邊界層的厚度變化.她們指出由于二維與三維情況下角渦數(shù)量的不同,導(dǎo)致了二者在溫度分布變化上的差異,但二維與三維流動中邊界層厚度都與Ra存在明確的標(biāo)度關(guān)系[45].在二維較高Ra(Ra=1013)下,包蕓等[46]發(fā)現(xiàn)Nu和LCS 受Ra的影響.他們發(fā)現(xiàn),兩者隨Ra的變化規(guī)律有較強的關(guān)聯(lián)性,而GL 理論可以部分準(zhǔn)確地預(yù)測Nu隨較高Ra的變化規(guī)律[46].另外,包蕓等[47]研究了Ra較高的條件下(Ra=1010),溫度邊界層厚度和Nu隨Pr變化的規(guī)律.其中溫度邊界層的厚度隨Pr的變化并不明顯[47].在較低Pr數(shù)下,Nu隨Pr上升而增加;而當(dāng)Pr數(shù)較高時,難以觀察到Nu隨Pr的變化[47].Wang等[48]通過實驗測量了準(zhǔn)二維系統(tǒng)中邊界層的溫度脈動的空間分布規(guī)律,并發(fā)現(xiàn)平均溫度分布與邊界層厚度存在一個普適的標(biāo)度關(guān)系,并指出邊界層方程對理解邊界層脈動提供了一個理論框架.另外,Zhou 和Xia[49]對下邊界的溫度邊界層進行了高精度測量,他們發(fā)現(xiàn)脈動溫度均方根(RMS)的最大值距離下邊界的距離,是一個可以非常好的描述溫度脈動統(tǒng)計特征的空間尺度.

    LSC 與黏性邊界層有著密切關(guān)系,Sun等[50]利用粒子圖像測速(particle image velocimetry,PIV)方法實驗測量了在大尺度環(huán)流作用下的邊界層特征,發(fā)現(xiàn)黏性邊界層厚度隨Re的標(biāo)度指數(shù)是-0.5.Zhou等[51]也分析了溫度與黏性邊界層剖面曲線,他們發(fā)現(xiàn)在使用適當(dāng)?shù)奈锢砹窟M行縮放后,RB 對流的溫度與黏性邊界層剖面與經(jīng)典Prandtl-Blasius(PB)層流邊界層符合得較好.最近,方明衛(wèi)等[52]對溫度邊界層剖面曲線與PBP (Prandtl-Blasius-Pohlhausen)解進行了擬合,找到了擬合參數(shù)隨Ra與Pr的變化規(guī)律.而在邊界層外的混合區(qū)中,脈動溫度的均方根(root mean square,RMS)剖面曲線與到邊界的距離滿足冪律函數(shù)的形式[50].Sun等[50]測量的脈動溫度RMS 曲線與Adrian 的理論研究[53]所預(yù)測的結(jié)果一致.而最近的實驗研究[54]也獲得了滿足冪律函數(shù)形式的脈動溫度RMS 分布曲線.

    對于RB 對流的能量耗散,學(xué)術(shù)界也取得了很多重要研究進展.在實驗研究方面,He等[31,32]首先成功測量出RB 對流中空間不同位置的時間平均擬熱能耗散率代表求時間平均.他們發(fā)現(xiàn),在溫度邊界層中受平均溫度的主導(dǎo),而在邊界層外的湍流體中受脈動溫度的主導(dǎo)[31,32].Ni等[55]使用拉格朗日粒子追蹤法直接測量了二階速度結(jié)構(gòu)函數(shù),并利用結(jié)構(gòu)函數(shù)間接獲得了當(dāng)?shù)貏幽芎纳⒙蚀韺臻g某處局部體積 δV求空間與時間平均.他們發(fā)現(xiàn),隨Ra的變化可以用冪律函數(shù)進行標(biāo)度[55].數(shù)值研究方面,Zhang等[29]研究了水(Pr=5.3)與空氣(Pr=0.7)中,動能耗散率和擬熱能耗散率的統(tǒng)計特性、空間分布規(guī)律以及與Ra的標(biāo)度關(guān)系.而Xu等[30]研究了極低Pr的條件下(Pr=0.025)擬熱能耗散率在湍流體與邊界層中的統(tǒng)計特征.他們發(fā)現(xiàn)盡管邊界層區(qū)域較小,但其對總的擬熱能與動能耗散占到支配作用,并建立了和與Ra的標(biāo)度關(guān)系[29,30].在理論研究方面,Petschel等[56]從能量耗散的角度定義了耗散邊界層,并推廣到了任意流動.

    本文借助RB 對流中所建立起的經(jīng)典研究方法與物理概念,利用DNS 可以直接獲得肥皂泡上完整的速度場和溫度場信息的優(yōu)勢[57],結(jié)合肥皂泡模型的特點,提出了肥皂泡溫度與黏性邊界層厚度的定義方法: 溫度邊界層的厚度為脈動溫度RMS 最大值所處緯度距離赤道的大地距離;而黏性邊界層厚度為緯度方向脈動速度RMS 曲線近赤道線性部分的延長線與其最大值的交叉緯度距離赤道的大地距離.同時本文也對肥皂泡模型中的能量耗散與邊界層特征進行分析與研究,比較和RB 對流模型的異同.相比RB 對流,肥皂泡熱對流模型中沒有明顯的LSC 與角渦的存在,且浮力強弱隨著緯度而發(fā)生變化,其球面特征還使得不同緯度上空間的熱通量強度不同.通過驗證RB 對流中物理規(guī)律是否可以用于描述肥皂泡模型的特征,既可以促進充分認識肥皂泡模型的物理機制,同時也檢驗RB 對流模型中理論的通用性.本文的結(jié)構(gòu)如下: 第2 節(jié)詳細介紹DNS 所使用的數(shù)學(xué)模型;第3 節(jié)選取了多個不同數(shù)量級的Ra和羅斯比數(shù)開展數(shù)值研究,探索了肥皂泡上的能量輸運和流場特征,Nu與Re隨Ra以及 1/Ro的變化規(guī)律,同時還討論了肥皂泡上邊界層和能量耗散特征;第4 節(jié)進行總結(jié)分析.

    2 DNS 數(shù)值模擬方法

    2.1 數(shù)值模型

    其中u和T分別為量綱一的速度場與溫度場.肥皂泡控制方程的量綱一參數(shù)為Ra,Pr和Ro,它們的定義分別為

    1/Ro的大小與轉(zhuǎn)速成反比: 1/Ro越大,旋轉(zhuǎn)的作用就越大.

    與經(jīng)典的Rayleigh-Bénard 方程相比,肥皂泡上流場的控制方程多了兩項: 外冷卻項-ST和外摩擦項-Fu.S和F分別為量綱一冷卻系數(shù)與摩擦系數(shù),外冷卻項與外摩擦項表征了肥皂泡模型與外部環(huán)境中冷空氣產(chǎn)生熱交換與摩擦.當(dāng)S=0且F ?=0時,肥皂泡模型的能量方程和Oberbeck-Boussinesq 方程組中的能量方程一樣.由于肥皂泡只有赤道一條恒溫?zé)徇吔?能量單向通過赤道流入肥皂泡,使肥皂泡上的溫度隨著時間的增加而不斷增加,這一過程為非平衡過程.最終肥皂泡上的溫度處處與赤道溫度相同時,肥皂泡停止從赤道吸收能量而達到平衡態(tài),此時流體失去浮力驅(qū)動形成熱寂[4,14].冷卻項使肥皂泡上的平均溫度始終小于赤道溫度,在適當(dāng)?shù)腟取值范圍下,可使極地平均溫度在數(shù)量級上遠低于赤道溫度,近似為0.肥皂泡的赤道與頂部的溫差近似保持恒定,使熱對流在恒定的Ra下可以一直進行下去,并使肥皂泡進入統(tǒng)計穩(wěn)平衡態(tài)[4,14].而如果S ?=0 且F=0 時,肥皂泡的DNS 隨著仿真時間步的推進而無法避免地產(chǎn)生爆炸[11,4,14].肥皂泡上的湍流熱對流屬于二維湍流,存在能量逆級串的現(xiàn)象[11,4,14].能量逆級串使流體從赤道吸收的熱能以羽流相互作用的方式形成大尺度的經(jīng)向流[2,11],最終不斷快速匯聚的經(jīng)向動能使得CFL 條件永遠無法得到滿足從而使DNS 仿真失敗[11].向控制方程中添加外摩擦項,其意義在于抑制肥皂泡上大尺度流動結(jié)構(gòu)的能量不受限制的積累,使肥皂泡進入統(tǒng)計平衡態(tài).在以往的二維湍流DNS 的研究中,添加外摩擦項也是一種普遍采用的方法[2].以往的研究已經(jīng)表明,S=F=0.06 及其所在的量級可以較好地模擬肥皂泡上的統(tǒng)計平衡熱對流[4,14],此時北極部分的平均量綱一溫度相比赤道量綱一溫度T0=1 有多個數(shù)量級的差異,可近似為0.同時本文也會對兩者帶來的外耗散影響進行討論,以驗證外冷卻項與外摩擦項對流動的作用不會影響系統(tǒng)的物理機制.

    以肥皂泡的球心為原點建立量綱一三維笛卡爾坐標(biāo)系 (x1,x2,x3),x1和x2坐標(biāo)軸均在赤道平面.三個方向的基矢分別為e1,e2和e3.e3的方向為從原點指向肥皂泡的球頂.肥皂泡模型的邊界為赤道,赤道上保持恒溫?zé)o滑移邊界條件(u=0與T=1).肥皂泡的曲面方程為

    而描述赤道的曲線方程為

    量綱一三維笛卡爾坐標(biāo)系 (x1,x2,x3) 所對應(yīng)的球坐標(biāo)系為(r,?,ψ) :r為徑向長度,?為方位角,ψ為極角.球坐標(biāo)系三個方向的單位基矢分別為er,e?和eψ.為了方便描述數(shù)值計算結(jié)果,引入地理坐標(biāo)系中的經(jīng)度與緯度來描述肥皂泡上任意一點的空間位置: 經(jīng)度坐標(biāo)與方位角定義相同,也用?表示;緯度坐標(biāo)是極角ψ的余角,用θ=π/2-ψ表示.在赤道上緯度θ為0°,而極角ψ=90°;在肥皂泡的頂點緯度θ為9 0°,而極角ψ=0°.經(jīng)度方向基矢與方位角方向基矢相同,為e?.緯度方向基矢與極角方向基矢方向相反,eθ=-eψ.

    2.2 數(shù)值方法

    控制方程(8)是向量形式,在求解時需要根據(jù)坐標(biāo)系確定算符的具體形式.由于肥皂的幾何外形為二維半球曲面,使用笛卡爾坐標(biāo)系求解控制方程(8)需要額外的方法處理曲面,增加了求解難度.如果使用球坐標(biāo)系求解,眾所周知,球坐標(biāo)系的NS 方程在極點處引入了奇異性,因而必須通過其他方式回避這一問題.球極投影坐標(biāo)系 (x,y) 可以解決肥皂泡曲面幾何所帶來的難題.球極投影的幾何關(guān)系如圖1 所示.肥皂泡球面上任意A點與南極點B之間的連線在赤道平面上唯一相交于點C,通過該一一映射關(guān)系可以將球面上所有點映射在一個平面圓上.而球極投影坐標(biāo)系的原點x=y=0也處于肥皂泡的球心.x方向與y方向分別與笛卡爾坐標(biāo)系中x1方向和x2方向重合.

    圖1 不同坐標(biāo)系下的計算域與球極投影平面投影示意圖Fig.1.Computational domain under different coordinate system and the illustration of the stereographical projection.

    通過簡單的幾何關(guān)系推導(dǎo),可以得到笛卡爾坐標(biāo)系 (x1,x2,x3) 與球極投影坐標(biāo)系 (x,y) 的變換為

    而逆變換為

    為獲得球極投影坐標(biāo)系中的矢量算符的具體形式,首先需要計算球極投影坐標(biāo)系的拉梅系數(shù).x和y兩個方向的拉梅系數(shù)可寫為

    由于兩個方向的拉梅系數(shù)相同,下文定義L=Lx=Ly.球極投影坐標(biāo)系在x與y方向上的單位基矢分別是hx與hy,分別為

    速度矢量在球極投影坐標(biāo)系中可以寫為

    其中u和v分別為速度矢量在x和y方向的投影.速度矢量在笛卡爾坐標(biāo)系、球坐標(biāo)系和地理坐標(biāo)系中的各個分量為

    由于肥皂泡為二維半球面,徑向速度恒為0:ur ≡0.而根據(jù)基矢關(guān)系可以知道,uθ=-uψ.

    根據(jù)曲線坐標(biāo)系中梯度算符的通式,可以寫出?算符在球極投影坐標(biāo)系中的形式:

    而拉普拉斯算符 Δ 在球極投影坐標(biāo)系中為

    拉格朗日導(dǎo)數(shù)為

    在獲得以上單位基矢與算符的轉(zhuǎn)換關(guān)系以后,可以將球極投影坐標(biāo)系內(nèi)的算符代入到矢量形式的控制方程中,從而獲得球極投影坐標(biāo)系下的控制方程.連續(xù)性方程可以寫為

    動量方程為

    能量方程為

    在球極投影坐標(biāo)系 (x,y) 中,肥皂泡的曲面方程為x2+y2≤1,是一個半徑為1 的圓.將計算域設(shè)定為邊長為1.01 的正方形,使用交錯結(jié)構(gòu)化網(wǎng)格來離散計算域: 溫度與壓強結(jié)點在網(wǎng)格的中心,而速度結(jié)點位于網(wǎng)格邊的中點.網(wǎng)格結(jié)點分布如圖2 所示,其中i和j分別為x和y方向上的結(jié)點編號.

    圖2 網(wǎng)格結(jié)點分布圖Fig.2.Node distribution inside a grid cell.

    邊界條件使用懲罰法來實現(xiàn),時間離散格式為二階Gear 格式,而非線性項使用三階Murmanlike 格式來離散,更詳細的數(shù)值方法與驗證可參考文獻[4,14].

    2.3 算例設(shè)置

    表1 列出了本文肥皂泡模型DNS 算例的全部參數(shù)信息.由于肥皂水中肥皂的濃度很低,其物性與水非常接近,算例的Pr固定為7 不變.算例的Ra的變化范圍為3×106—3×109,小Ra下采用1024×1024 的數(shù)值解析度,在最高Ra下網(wǎng)格分辨率達到了 2 048×2048,以往的研究表明這一數(shù)值解析度可以充分獲得網(wǎng)格獨立解[4,11,14].

    表1 算例詳細信息Table 1.The detailed information of the cases.

    2.4 統(tǒng)計穩(wěn)態(tài)

    在數(shù)值實驗中,肥皂泡的初始狀態(tài)為靜止態(tài),即u=0 且T=0.隨著計算的進行,時間步不斷向前迭代,肥皂泡上的流體開始受熱并運動,產(chǎn)生湍流熱對流逐漸達到統(tǒng)計穩(wěn)態(tài)[4,14].此時肥皂泡上的湍流得到了充分發(fā)展,肥皂泡吸收的能量與其耗散的能量達到了動態(tài)平衡.在這樣的條件下,流場量可以分解為平均值與脈動值之和,如溫度T:

    在下文,符號〈·〉表示時間和經(jīng)度空間平均.肥皂泡在經(jīng)度?方向具有旋轉(zhuǎn)對稱性,流場變量在統(tǒng)計上并不依賴經(jīng)度坐標(biāo)?.于是,在控制參數(shù)確定的情況下,所有流場變量的平均值是緯度θ的函數(shù).

    3 結(jié)果與討論

    3.1 瞬態(tài)流場

    對所研究的不 同Ra都考慮了4 個 1/Ro: 0,0.1,1 和 10,以便分析旋轉(zhuǎn)作用的影響.這里首先從最直觀的流場瞬時云圖來展示不同Ra與不同旋轉(zhuǎn)角速度 1/Ro給流場帶來的影響.圖3 給出了不同Ra與 1/Ro條件下溫度場T在肥皂泡上的瞬時分布云圖.當(dāng)Ra=3×106時,在沒有旋轉(zhuǎn)的情況下(1/Ro=0),可以從瞬時溫度場中清晰地觀察到羽流的形狀特點: 羽流尺寸大,而數(shù)量相對少.隨著Ra升高至Ra=3×109,肥皂泡表面的羽流的尺寸大幅減小,且更加卷曲,同時羽流數(shù)量大幅增加.在肥皂泡高速旋轉(zhuǎn)的情況下(1/Ro=10),羽流的形狀與尺度大小同未旋轉(zhuǎn)的情況下相比并無明顯變化,但是羽流造成的影響難以上升到肥皂泡頂部區(qū)域.

    圖3 不同 1/Ro 與 Ra 下,T 在肥皂泡上的分布 (a) 算例 Ra=3×106,1/Ro=0;(b) 算例 Ra=3×109,1/Ro=0;(c) 算例 Ra=3×109,1/Ro=10Fig.3.The instantanous temperature field with different 1/Ro and Ra : (a) The case of Ra=3×106,1/Ro=0;(b) the case of Ra=3×109,1/Ro=0;(c) the case of Ra=3×109,1/Ro=10.

    同時,圖4 給出了在不同Ra與 1/Ro設(shè)置下,肥皂泡上的瞬時動能的對數(shù) l gEk=lg(1/2u2).在肥皂泡沒有旋轉(zhuǎn)的情況下,在Ra=3×106時,可通過區(qū)分動能集中的區(qū)域識別出流動結(jié)構(gòu)的尺寸較大;而當(dāng)Ra=3×109時,流動結(jié)構(gòu)的尺寸明顯變小,同時變得更加破碎.如果保持Ra大小不變,在高速旋轉(zhuǎn)時 1/Ro=10,肥皂泡上靠近赤道的部分出現(xiàn)了大尺寸且動能集中的流動結(jié)構(gòu),這是旋轉(zhuǎn)帶來的穩(wěn)定經(jīng)向流.另一方面,在肥皂泡的頂部附近,流動結(jié)構(gòu)的動能較小,但尺寸相比沒有旋轉(zhuǎn)的情況更大.

    圖4 不 同 1/Ro 與 Ra 下,log Ek 在肥皂泡上的分布 (a)算例 Ra=3×106,1/Ro=0;(b)算例 Ra=3×109,1/Ro=0;(c)算例 Ra=3×109,1/Ro=10Fig.4.The instantanous filed of log Ek with different 1/Ro and Ra : (a) The case of Ra=3×106,1/Ro=0;(b) the case of Ra=3×109,1/Ro=0;(c) the case of Ra=3×109,1/Ro=10.

    3.2 Nu 與Re

    在經(jīng)典RB 熱對流模型中,熱能從下邊界注入,從上邊界流出,Nu被用來表征流體對流輸運熱能的效率.但在肥皂泡模型中,只有赤道邊界:熱量通過這條邊界流入肥皂泡,然后通過外部和內(nèi)部耗散在肥皂泡上消耗殆盡.本文使用Nu來表征肥皂泡赤道上的對流傳熱效率,Nu同時也代表了肥皂泡從外界吸收能量的效率,Nu越大意味著肥皂泡通過對流從外界吸收熱能的能力越強.于是,Nu的定義為赤道處的無量綱熱通量,其形式如下:

    在熱對流中,溫度較高的流體受到浮力的作用從赤道向肥皂泡的頂部運動,將熱能與動能輸運至肥皂泡內(nèi)部.但是,可以證明在經(jīng)度方向上的平均速度為0.首先寫出球坐標(biāo)系下的平均連續(xù)性方程:

    圖6 給出了Nu與Re隨Ra變化的規(guī)律.從圖6可以看出,Nu隨Ra變化近似滿足標(biāo)度率~Ra13,這與許多RB 對流中所觀察到的結(jié)果非常接近.圖6 中結(jié)果表明,1/Ro的變化對Nu的影響非常小,只有在Ra較低(Ra≤3×107),并且1/Ro=

    圖6 Nu(上)和 Re(下)隨 Ra(左)和 1/Ro(右)的變化規(guī)律 (a) 在 1/Ro 恒定的條件下,Nu 隨 Ra 的標(biāo)度規(guī)律;(b) 在 Ra 恒定的條件下,Nu 隨 1/Ro 的變 化規(guī) 律;(c) 在 1/Ro 恒定的條件下,Re 隨 Ra 的標(biāo)度規(guī) 律;(d) 在 Ra 恒定的 條件 下,Re 隨 Ra 的變化規(guī)律Fig.6.The variation of Nu(up) and Re(down) with Ra(left) and 1/Ro(right): (a) The scaling behavior of Nu with Ra in condition of fixed 1/Ro;(b) the variation of Nu with 1/Ro in condition of fixed Ra;(c) the scaling behavior of Re with Ra in condition of fixed 1/Ro;(d) the variation of Re with 1/Ro in condition of fixed Ra.

    10 的條件下,Nu有很微弱的上升.綜合全部結(jié)果可以發(fā)現(xiàn),1/Ro對Nu標(biāo)度率的影響可以忽略.從定義公式(33)可以知道,Nu的大小主要取決于邊界層中的平均溫度分布曲線.圖5 數(shù)據(jù)表明,越靠近赤道,平均溫度隨著旋轉(zhuǎn)速度的變化越微弱[14].這暗示科里奧利力的作用隨著緯度的減小而減小.如果將速度u=u?e?+uθeθ代入動量方程的科里奧利力項,則可以得到:

    圖5 不同 1/Ro 和 Ra 下,平均溫度 〈T〉 隨緯度 θ 的變化規(guī)律 (a) 算例 Ra=3×106;(b) 算例 Ra=3×107;(c) 算例Ra=3×108;(d) 算例Ra=3×109Fig.5.The variation of mean temperature 〈T〉 with the latitude θ for the different 1/Ro and Ra : (a) The cases of Ra=3×106;(b)the cases of Ra=3×107;(c) the cases of Ra=3×108;(d) the cases of Ra=3×109.

    當(dāng)緯度為0 時,eθ=e3所以e3×eθ=0;另一方面,在赤道上有e3×e?=er;所以赤道上科里奧利力項可以寫為.由于肥皂泡二維半球面的幾何特征,徑向上的力對流體的運動沒有任何作用.在赤道上,科里奧利項的作用為0,而且在靠近赤道的區(qū)域內(nèi),科里奧利力的作用都很弱.這就揭示了旋轉(zhuǎn)對Nu的影響非常微小的原因.肥皂泡上Nu隨 1/Ro的變化規(guī)律與旋轉(zhuǎn)RB 對流有較明顯的不同.旋轉(zhuǎn)RB 對流中Nu隨著 1/Ro的逐漸增加呈現(xiàn)先升后降的變化規(guī)律: 在合理的參數(shù)范圍內(nèi),1/Ro可以使Nu最大增加30%[41,40].在適當(dāng)?shù)男D(zhuǎn)速度下,旋轉(zhuǎn)RB 對流中的羽流變?yōu)樨Q直排列的渦,產(chǎn)生Ekman 抽吸效應(yīng),將邊界層中的流體向湍流體中輸運,使Nu有明顯的上升[43,42].Ekman抽吸效應(yīng)的效率隨著Pr數(shù)的上升而得到加強,而隨著Ra數(shù)的上升Nu被削弱[40].但是只有在三維流場中才會產(chǎn)生Ekman 抽吸效應(yīng),所以肥皂泡上Nu上升的機理與RB 對流具有較大不同.

    另一方面,Re隨著Ra的變化曲線也具有冪律函數(shù)的形式.當(dāng)肥皂泡靜止或者旋轉(zhuǎn)速度較小時(1/Ro≤1),Re滿足標(biāo)度率Re~Ra0.48,而 且Re隨1/Ro的變化可以忽略.但當(dāng) 1/Ro=10 時,在較低Ra的條件下,Re出現(xiàn)了明顯的下降,而在較高Ra的條件下,Re隨 1/Ro的變化依然較小,這導(dǎo)致Re的標(biāo)度系數(shù)從 0.48 上升到 0.51.Re同時受到Ra與1/Ro的影響,與隨緯度分布的規(guī)律直接相關(guān)[14].在之后討論黏性邊界層的小節(jié)會展示,當(dāng)1/Ro=10時,緯度方向速度的脈動受到極大的抑制,max()也 隨之下降.但 隨著Ra上升,m ax() 下降幅度逐漸縮小,使Re的降幅也變小.

    3.3 均方根溫度分布與溫度邊界層

    RB 對流的溫度邊界層厚度δT可以根據(jù)平均溫度〈T〉來定義,也可以根據(jù)脈動溫度的均方根來定義[49].RB對流的實驗發(fā)現(xiàn),當(dāng)非??拷舷缕桨鍟r,流體的平均溫度〈T〉與到邊界的距離d之間的變化曲線表現(xiàn)出較好的線性規(guī)律.對線性分布區(qū)域做擬合線并且延長,當(dāng)擬合線達到湍流體區(qū)域的平均溫度時,其與邊界的距離可定義為邊界層的厚度[50,49].但是,肥皂泡對流模型中只有一個熱邊界,所以在邊界的法線方向上流場并不對稱,且靜止肥皂泡上的經(jīng)向流非常微弱[14],這導(dǎo)致流場中不存在一個有明確物理定義的湍流體區(qū)域.RB 對流還可以用另外一種方式定義δT: 脈動溫度的均方根達到最大值處到邊界的距離[29,50,49].因此,本文根據(jù)肥皂泡的曲面外形特點,采取以下方式定義溫度邊界層厚度: 均方根脈動溫度T*是θ的函數(shù),當(dāng)T*達到最大值時,θ=θT.緯度θT距離赤道的大地距離(即半球面上由赤道到緯度θT的最短距離: 圓周角為θT的大圓弧長)定義為溫度邊界層厚度δT.在無量綱坐標(biāo)系下,肥皂泡的半徑R ≡1,所以δT和θT之間的關(guān)系為θT=δT/R=δT.在肥皂泡上,緯度θT的物理含義與δT相同,代表了流體到邊界(赤道)的距離,兩者的差別在于量綱.

    圖7 的左邊給出了Ra=3×106,Ra=3×107,Ra=3×108和Ra=3×109的靜止肥皂泡上T*隨緯度的分布曲線,并且用虛線標(biāo)注了θT,即邊界層的厚度.圖7 的右邊則給出不同羅斯比數(shù)下邊界層厚度δT與Ra的關(guān)系.圖7 中的數(shù)據(jù)顯示,δT隨Ra變化呈指數(shù)標(biāo)度,而 1/Ro對溫度邊界層厚度的影響并不明顯.δT與Ra之間的標(biāo)度關(guān)系為δT~Ra-0.32.和3.2 節(jié)討論的原因一樣,動量方程中科里奧利的作用項隨著緯度的降低而減小,在赤道處變?yōu)?.而邊界層厚度很薄,而且隨著Ra的上升而快速下降,所以邊界層厚度受旋轉(zhuǎn)作用的影響非常微弱.值得一提的是,在RB 對流模型的研究文獻中,溫度邊界層厚度的標(biāo)度規(guī)律δT~已經(jīng)多次被研究者發(fā)現(xiàn)[45,49].在肥皂泡模型中,δT與Ra的標(biāo)度關(guān)系與RB 對流非常相似.并且,肥皂泡的δT標(biāo)度系數(shù)與Nu的標(biāo)度系數(shù)完全相反,這也符合RB 對流中的特征.從幾何角度觀察,溫度邊界層靠近赤道且厚度很薄,流動尺度很小.之前的研究也表明,在邊界層中受無滑移邊界條件的影響,流體的運動速度較低且溫度的脈動也很微弱[4,14],于是人工散熱與人工摩擦項的作用在邊界層中不明顯.同時,邊界層中浮力的投影與重力的夾角可以近似為1 80°,與RB 對流中的情況相同,這就說明了肥皂泡在邊界層中熱能輸運的物理機制與RB 對流很相近,使得δT隨著Ra的變化規(guī)律也與RB 對流接近.同時,由于溫度邊界層之外,溫度梯度較低;因而對于這一系統(tǒng),溫度邊界層的厚度也恰好與系統(tǒng)邊界層的溫度梯度成反比關(guān)系,因而其與Nu成反比關(guān)系.

    圖7 (a)均方根溫度剖面分布與溫度邊界層定義示意圖;(b)不同 1/Ro 下,溫度邊界層厚度隨 Ra 變化規(guī)律Fig.7.(a) The variation of T* profile with Ra;(b) the thickness of the thermal boundary layers with Ra and 1/Ro.

    圖8為將不同算例的T*曲線用δT與max(T*)進行歸一化后的結(jié)果.在保持 1/Ro恒定的條件下,圖8 的每個子圖中,Ra從 3×106增加至 3×109,提升了3 個數(shù)量級.結(jié)果顯示,當(dāng)旋轉(zhuǎn)作用較小時(1/Ro≤1),即使Ra增加了3 個數(shù)量級,不同Ra的T*歸一化曲線相互重合得較好,表現(xiàn)出較強的自相似性.當(dāng)旋轉(zhuǎn)作用達到最大時(1/Ro=10),在邊界層中與靠近邊界層的部分((Rθ/δT)≤2),不同Ra下的T*歸一化曲線依然重合得較好.產(chǎn)生這種現(xiàn)象的原因在于科里奧利力在赤道附近作用微弱,如3.2 節(jié)所討論.在離邊界層較遠的區(qū)域(2 ≤Rθ/δT≤4),除Ra=3×106外,其他Ra下的T*歸一化曲線依然展現(xiàn)出較好的自相似性.在離邊界層非常遠的區(qū)域(Rθ/δT≥4),只有Ra=3×108與Ra=3×109的T*歸一化曲線重合得較好.由此可以發(fā)現(xiàn),在科里奧利力作用明顯的條件下,隨著Ra上升,T*的歸一化曲線的重合度逐漸增加;而當(dāng)緯度越靠近赤道,T*曲線的重合度越明顯.在較強的旋轉(zhuǎn)作用下,圖中Ra=3×106在高緯度的差異可能是由于隨著Ra的降低,其羽流尺度相對更大,當(dāng)羽流延展到遠離邊界層區(qū)域時,科里奧利力作用逐漸增加,從而造成了羽流結(jié)合科里奧利力導(dǎo)致的溫度分布變化.隨著Ra上升,熱羽流尺寸相對減小,這一耦合科里奧利力篡改溫度分布的能力減弱,因而曲線的重合度更高.這些都表明,科里奧利力項的作用會改變T*的空間分布.另外,圖8的結(jié)果也表明,m ax(T*) 是一個非常合適的歸一化因子.m ax(T*) 直接描述了溫度脈動的強弱,與Ra相關(guān),于是可以較好地表征T*的自相似特性[49].

    圖8 不 同 1/Ro 和 Ra 條件下,使 用 δT 與 m ax(T*) 進行歸一化之后的 T*曲線 (a) 1/Ro=0;(b) 1/Ro=0.1;(c) 1/Ro=1;(d)1/Ro=10Fig.8.The RMS temperature distribution normalized by δT and m ax(T*) for the different 1/Ro and Ra : (a) 1/Ro=0;(b) 1/Ro=0.1;(c) 1/Ro=1;(d) 1/Ro=10.

    圖9 采用對數(shù)坐標(biāo)軸顯示了不同Ra下T*在邊界層外側(cè)的隨θ變化的曲線,每個子圖的結(jié)果具有不同的 1/Ro并且使用冪律函數(shù)對T*進行了擬合,擬合結(jié)果用不同線型的黑線展示.對比擬合結(jié)果與T*曲線可以發(fā)現(xiàn),在本文研究的羅斯比倒數(shù)范圍內(nèi)(0 ≤1/Ro≤10),邊界層外的T*(θ) 都具有冪律函數(shù)的形式.在 1/Ro保持不變的條件下,盡管Ra變化范圍達到了3 個數(shù)量級,但T*都在一定緯度范圍內(nèi)近似滿足T*~θξ.當(dāng)肥皂泡旋轉(zhuǎn)作用較弱時(1/Ro≤1),近似滿足T*~θξ的區(qū)間超過了一個數(shù)量級,擬合指數(shù)ξ與非常接近.但當(dāng)肥皂泡高速旋轉(zhuǎn)時(1/Ro=10),近似滿足T*~θξ的區(qū)間長度縮短,小于一個數(shù)量級,并且ξ與的差距明顯增加.對于所有算例,在θ≥1 時,流體非常靠近肥皂泡頂,脈動溫度快速下降[14].于是T*(θ)不再具有冪律函數(shù)的形式.肥皂泡上的T*在邊界層外的空間分布與RB 對流也具有一定的相似性.RB 對流中,邊界層外可以分為對流區(qū)(或混合區(qū))和對流核心區(qū)[58].Adrian[53]由慣性力與黏性力相互平衡的假定,推導(dǎo)得出對流區(qū)T*(d) 具有冪律函數(shù)的形式:T*~.多項RB 對流的實驗研究也發(fā)現(xiàn)了T*的空間分布曲線具有冪律函數(shù)形式.然而Adrian 基于慣性力與浮力平衡的假設(shè)預(yù)測T*(d)同樣可以具有對數(shù)函數(shù)的形式T*~lnd,并在實驗中得到驗證.最近,何宇昊與夏克青[54]使用高精度的實驗方法測量了RB 對流槽中剪切主導(dǎo)區(qū)與羽流主導(dǎo)區(qū)的T*(d).剪切主導(dǎo)區(qū)中,大尺度環(huán)流帶來的橫向平均速度使流場處于強剪切狀態(tài)中.而羽流主導(dǎo)區(qū)中,大量的羽流相互混合聚集并從熱邊界上升至冷邊界.他們的研究表明,剪切主導(dǎo)區(qū)內(nèi)T*~,而羽流主導(dǎo)區(qū)內(nèi)T*~lnd.比較RB 對流和肥皂泡模型的結(jié)果可以發(fā)現(xiàn),肥皂泡模型與RB 對流在邊界層外具有較大的物理差異.由于大尺度環(huán)流非常微弱,于是肥皂泡上不存在與剪切主導(dǎo)區(qū)中類似的流動狀態(tài).另一方面,羽流在肥皂泡的赤道上隨機生成并不斷運動,羽流的混合隨機發(fā)生在肥皂泡上不同的區(qū)域,羽流主導(dǎo)區(qū)中類似的流動狀態(tài)也不存在于肥皂泡上.但是,肥皂泡上添加了外摩擦項,引入了較大的外摩擦黏性力;而隨著θ的上升,在肥皂泡緯度方向上重力的投影越來越小.這些都暗示T*~θξ的分布規(guī)律可能取決于慣性力與黏性力相互平衡機制,但是詳細物理圖像還需要更加深入的研究.

    圖9 不同 1/Ro 和 Ra 條件下,T* 的剖 面曲線 (a) 1/Ro=0;(b) 1/Ro=0.1;(c) 1/Ro=1;(d)1/Ro=10Fig.9.The profile of T* for the different 1/Ro and Ra : (a) 1/Ro=0;(b) 1/Ro=0.1;(c) 1/Ro=1;(d) 1/Ro=10.

    3.4 緯度脈動速度均方根的空間分布與黏性邊界層

    圖10 不 同 Ra 與 1/Ro 條件下,均方根緯度速度 緯度剖面曲線 (a) Ra=3×106;(b) Ra=3×107;(c) Ra=3×108;(d)Ra=3×109Fig.10.The profile of the RMS velocity in the latitude direction for the differnet Ra and 1/Ro : (a) Ra=3×106;(b)Ra=3×107;(c) Ra=3×108;(d) Ra=3×109.

    當(dāng)肥皂泡上流體的湍流運動得到充分發(fā)展之后,流場的平均速度相比脈動速度非常微弱[14].之前的章節(jié)已經(jīng)證明,緯度方向的平均速度〈uθ〉=0;靜止肥皂泡上的大尺度環(huán)流非常微弱,其經(jīng)度方向的平均速度同樣非常微弱[14].因而,不便于通過平均速度場來定義黏性邊界層.本文參考RB 對流的相關(guān)研究文獻[50,51],基于緯度方向脈動速度的均方根來定義黏性邊界層厚度δu.圖11的左側(cè)示例了通過曲線來定義黏性邊界層厚度的方法.對于肥皂泡對流模型,隨θ變化的曲線在靠近赤道的部分具有線性特征.將線性部分做擬合線并延長,當(dāng)擬合線達到脈動緯度速度的均方根最大值 m ax() 時,緯度為θu.θu距離赤道的大地距離R×θu就是黏性邊界層的厚度δu.

    圖11 (a)黏性邊界層厚度 δu 的定義方法(示意算例Ra=3×107,1/Ro=0);(b)不 同 1/Ro 下 δu 隨Ra變化規(guī)律Fig.11.(a) The definition of the viscous boundary layer thicknesses δu with the example case of Ra=3×107 and 1/Ro=0;(b) the variation of δu with Ra for the different 1/Ro.

    圖11 的右邊給出了δu隨Ra的變化規(guī)律.對于肥皂泡熱對流模型,δu與Ra的關(guān)系呈現(xiàn)冪律函數(shù)的形式.而隨著Ra增加,黏性邊界層的厚度快速減小.對于靜止肥皂泡,1/Ro=0,δu近似滿足標(biāo)度率δu~Ra-0.2.當(dāng) 1/Ro≤1 時,肥皂泡上旋轉(zhuǎn)作用較弱,δu隨 1/Ro的變化可以忽略,此時δu隨Ra變化的標(biāo)度率保持恒定.當(dāng) 1/Ro=10 時,肥皂泡高速旋轉(zhuǎn),Ra相對較低時(Ra≤3×108),δu明顯下降.但隨著Ra升高到 3×109,1/Ro對δu的影響可以忽略.結(jié)果,當(dāng) 1/Ro=10 時,δu隨Ra變化的標(biāo)度指數(shù)從-0.2 變?yōu)?0.14.RB 對流中,由于大尺度環(huán)流的存在,使用水平脈動速度的均方根曲線來定義δu.在早期的文獻中,Xin 和Xia[59]首先在圓柱形對流槽中測量得到δu~Ra-0.25,而Qiu 和Xia[60]測得δu~Ra-0.22.這與肥皂泡中的結(jié)果較為接近.最近,孫超等[50]在使用PIV 方法測量了矩形對流槽的黏性邊界層,并得到δu~Ra-0.37.

    3.5 擬熱能耗散率和動能耗散率

    令單位質(zhì)量流體所包含的擬熱能和動能分別為ET=T2/2與Ek=.將能量方程和動量方程兩邊分別乘以T或者u,就可以得到ET與動能Ek的控制方程:

    圖12 不同 Ra 與 1/Ro 條件下,擬熱能內(nèi)耗散率 在緯度方向的分布 (a) Ra=3×106;(b) Ra=3×107;(c)Ra=3×108;(d) Ra=3×109.子圖為邊界層附近的放大Fig.12.The distribution of the internal thermal energy dissipation rate in the latitude direction for the different 1/Ro and Ra :(a) Ra=3×106;(b) Ra=3×107;(c) Ra=3×108;(d) Ra=3×109.The insets are the zoom-in for the boundary layers.

    隨著Ra的上升,δT快速減小,溫度邊界層中和隨θ下降得更快.肥皂泡在邊界層中的平均溫度梯度最大[14],于是擬熱能內(nèi)耗散集中在溫度邊界層中,而在肥皂泡的其他區(qū)域,擬熱能內(nèi)耗散較小.在邊界層中,科里奧利力項的作用非常小,1/Ro的上升對擬熱能耗散率的影響可以忽略.

    圖14 給出了不同Ra和 1/Ro參數(shù)組合下,和隨緯度θ變化的曲線,并且使用黑色和紅色虛線表示不同 1/Ro下黏性邊界層的厚度δu.因為1/Ro≤1 的條件下,δu隨 1/Ro的變化可以忽略,于是圖14 使用黑色虛線標(biāo)識δu(1/Ro=0),而使用紅色虛線標(biāo)識δu(1/Ro=10).可以看出,隨著Ra數(shù)上升,δu(1/Ro=10) 與δu(1/Ro=0)之間的差別不斷減小,最終在Ra=3×109時相互重合.從圖14 還可以發(fā)現(xiàn)動能外耗散率比內(nèi)耗散率要小一個量級,這暗示肥皂泡與外部的相互作用不會給湍流運動產(chǎn)生最本質(zhì)的影響.此外,〈Tu3〉在邊界層中隨緯度快速增長,達到最大值后緩慢下降,在肥皂泡頂處為0.且隨著Ra的上升,邊界層變薄,〈Tu3〉的峰值越來越靠近赤道.在 1/Ro≤0.1 的條件 下,〈Tu3〉隨 1/Ro并沒有變化.而 1/Ro=1 時,〈Tu3〉在黏性邊界層外隨θ的上升有更明顯的下降.在1/Ro=10 時,〈Tu3〉在黏性邊界層外隨θ的增加下降非???當(dāng)θ≥45°時,〈Tu3〉的值約等于0.而在黏性邊界層 內(nèi),〈Tu3〉的變化更加復(fù) 雜.在Ra=3×106或Ra=3×107的情況下,〈Tu3〉在邊界層內(nèi)的峰值有明顯的上升.但隨著Ra增加到3×108或 3×109時,〈Tu3〉峰值不再隨 1/Ro變化.

    圖13 不同 Ra 與 1/Ro 條件下,擬熱能內(nèi)耗散率 在緯度方向的分布 (a) Ra=3×106;(b) Ra=3×107;(c)Ra=3×108;(d) Ra=3×109.內(nèi)插圖為邊界層附近的放大Fig.13.The distribution of the internal thermal energy dissipation rate in the latitude direction for the different 1/Ro and Ra :(a) Ra=3×106;(b) Ra=3×107;(c) Ra=3×108;(d) Ra=3×109.The insets are the zoom-in for the boundary layers.

    動能內(nèi)耗散率的平均在邊界層中有小幅的增長,之后隨著θ的增加快速下降,在黏性邊界層中形成一 個尖銳的峰值.當(dāng)θ≥30°時,變 得非 常接 近于 0.隨 著Ra的增 加,在黏性邊界層中隨θ上升下降更加劇烈,峰值變得越發(fā)尖銳,且不斷接近赤道.圖14 中可以觀察到低速旋轉(zhuǎn)(1/Ro≤ 1)對在黏性邊界層內(nèi)外的分布沒有明顯影響.但當(dāng)肥皂泡高速旋轉(zhuǎn)時(1/Ro=10),在黏性邊界層中的峰值內(nèi)有明顯的下降,但在黏性邊界層外又有輕微的增加.綜合圖14 中的結(jié)果可以發(fā)現(xiàn),在黏性邊界層中的值遠大于邊界層之外的值.

    圖14 不同 Ra 與 1/Ro 條件下,動能耗散率以及浮力項在 緯度 方向 的分 布 (a),(b),(c) Ra=3×106;(d),(e),(f)Ra=3×107;(g),(h),(i) Ra=3×108;(j),(k),(l) Ra=3×109.內(nèi)插圖為邊界層附近的放大Fig.14.The distribution of the kinetic energy dissipation rate and the buoyancy term in the latitude direction for the different 1/Ro and Ra : (a),(b),(c) Ra=3×106;(d),(e),(f) Ra=3×107;(g),(h),(i) Ra=3×108;(j),(k),(l) Ra=3×109.The insets are the zoom-in for the boundary layers.

    4 結(jié)論

    通過簡單加熱肥皂泡可以獲得一個新的熱對流模型,該模型具有典型的準(zhǔn)二維球面特征,相比RB 對流不存在由幾何邊界奇點形成的角渦,其浮力作用也隨緯度不同而變化.本文提出了加熱肥皂泡的數(shù)值模擬方法,并通過直接數(shù)值模擬研究了不同Ra與 1/Ro下旋轉(zhuǎn)肥皂泡上的邊界層特征和能量耗散規(guī)律.結(jié)合肥皂泡上的流動特征提出了相應(yīng)的Nu,Re與溫度和黏性邊界層厚度的定義.其中δT定義為m ax(T*) 所在緯度距離赤道的大地距離,而δu是曲線近赤道線性部分延長至 m ax() 時所處緯度距離赤道的大地距離.DNS 數(shù)值結(jié)果表明加熱肥皂泡上的湍流熱對流具有以下特點.

    1)當(dāng)肥皂泡無旋轉(zhuǎn)運動時,Nu與Re與Ra的標(biāo)度關(guān)系為Nu~Ra0.30和Re~Ra0.48;盡管熱對流模型與RB 對流不同,但其標(biāo)度關(guān)系與標(biāo)度指數(shù)相近,這進一步表明這一標(biāo)度關(guān)系對熱對流的普適性.旋轉(zhuǎn)對Nu的影響可以忽略,但當(dāng)旋轉(zhuǎn)作用很強時(1/Ro=10),旋轉(zhuǎn)在低Ra時對Re有一定的抑制,使得Re與Ra的標(biāo)度關(guān)系變化為Re~Ra0.50.

    2)溫度邊界層的厚度δT隨Ra的標(biāo)度關(guān)系為δT~Ra-0.32,且不隨羅斯比數(shù)的變化而變化.在溫度邊界層外側(cè),均方根溫度的剖面曲線與緯度θ的近似關(guān)系為而隨著肥皂泡旋轉(zhuǎn)作用的增加,滿足標(biāo)度關(guān)系的緯度區(qū)間逐漸縮短.

    3)黏性邊界層的厚度δu隨Ra的標(biāo)度關(guān)系為δu~Ra-0.20.只有當(dāng)旋轉(zhuǎn)作用足夠強(1/Ro=10),并且Ra相對較小(Ra≤3×108)時,δu出現(xiàn)明顯下降,使得與Ra標(biāo)度關(guān)系變?yōu)棣膗~Ra-0.14.

    猜你喜歡
    熱對流肥皂泡赤道
    “雪龍2”號過赤道,救下一船人!
    軍事文摘(2024年2期)2024-01-10 01:58:44
    十二月·肥皂泡泡
    距赤道越遠越喜愛黃色
    最熱的地方不在赤道
    不銹鋼板傳熱特性試驗研究
    座艙空氣非定常流動特征及數(shù)值模擬策略2012年研究進展報告
    肥皂泡為什么是圓形?
    吹肥皂泡
    基于空氣流體動力學(xué)的高速列車制動盤散熱性能模擬
    最熱的地方不在赤道
    奧秘(2015年3期)2015-09-10 07:22:44
    99re在线观看精品视频| 人人妻人人澡人人看| 国产无遮挡羞羞视频在线观看| 精品久久久精品久久久| 真人做人爱边吃奶动态| 日本vs欧美在线观看视频| 精品国产乱码久久久久久男人| 在线观看一区二区三区激情| 亚洲精品在线美女| 国产精品免费视频内射| 咕卡用的链子| 在线观看免费视频日本深夜| 在线观看免费视频网站a站| 女人久久www免费人成看片| 亚洲人成电影观看| 午夜久久久在线观看| 色播在线永久视频| 伊人久久大香线蕉亚洲五| 在线永久观看黄色视频| 亚洲黑人精品在线| 女人被狂操c到高潮| 国产精品电影一区二区三区 | 精品久久久久久久毛片微露脸| a在线观看视频网站| 久久久久国产一级毛片高清牌| 一级毛片高清免费大全| 丝瓜视频免费看黄片| 欧美激情 高清一区二区三区| 国产成人精品久久二区二区91| 免费在线观看亚洲国产| 免费在线观看日本一区| 老汉色∧v一级毛片| 另类亚洲欧美激情| 日韩大码丰满熟妇| 国产人伦9x9x在线观看| 国产视频一区二区在线看| 美女午夜性视频免费| 极品教师在线免费播放| 十分钟在线观看高清视频www| 丁香六月欧美| 成人特级黄色片久久久久久久| 在线观看免费视频网站a站| 欧美日韩乱码在线| 十分钟在线观看高清视频www| 啦啦啦视频在线资源免费观看| 日韩中文字幕欧美一区二区| 国产不卡av网站在线观看| а√天堂www在线а√下载 | 亚洲欧美色中文字幕在线| 久久香蕉国产精品| 亚洲人成电影免费在线| 久久青草综合色| 侵犯人妻中文字幕一二三四区| 国产又色又爽无遮挡免费看| 天天操日日干夜夜撸| 可以免费在线观看a视频的电影网站| 我的亚洲天堂| 少妇 在线观看| 亚洲国产精品sss在线观看 | a级毛片在线看网站| 国产单亲对白刺激| 欧美精品啪啪一区二区三区| 色在线成人网| 亚洲七黄色美女视频| 女人精品久久久久毛片| 不卡av一区二区三区| 婷婷精品国产亚洲av在线 | 国产精品久久久久成人av| 国产深夜福利视频在线观看| 国产精品自产拍在线观看55亚洲 | 免费黄频网站在线观看国产| av有码第一页| 操美女的视频在线观看| 日韩免费高清中文字幕av| 女人被躁到高潮嗷嗷叫费观| 巨乳人妻的诱惑在线观看| 黄网站色视频无遮挡免费观看| 午夜亚洲福利在线播放| 十八禁网站免费在线| 精品一区二区三卡| 丝瓜视频免费看黄片| 少妇粗大呻吟视频| 亚洲熟女毛片儿| 国产精品综合久久久久久久免费 | 美女高潮到喷水免费观看| 国产有黄有色有爽视频| 国产av又大| 国产亚洲欧美98| 高清在线国产一区| 免费一级毛片在线播放高清视频 | 99riav亚洲国产免费| 日本wwww免费看| 欧美激情高清一区二区三区| 欧美日韩精品网址| 成人精品一区二区免费| 丰满人妻熟妇乱又伦精品不卡| 久久久精品免费免费高清| 亚洲一区中文字幕在线| 国产成人免费无遮挡视频| 亚洲国产欧美网| 中文字幕精品免费在线观看视频| 一本大道久久a久久精品| 乱人伦中国视频| 日韩视频一区二区在线观看| 777米奇影视久久| 超碰成人久久| av线在线观看网站| 叶爱在线成人免费视频播放| 成人国语在线视频| 国产精品98久久久久久宅男小说| 午夜福利视频在线观看免费| 欧美av亚洲av综合av国产av| 99re在线观看精品视频| 国产精品一区二区在线观看99| 精品人妻在线不人妻| 色老头精品视频在线观看| 757午夜福利合集在线观看| 亚洲成人免费电影在线观看| 国产淫语在线视频| 男女床上黄色一级片免费看| 51午夜福利影视在线观看| 黄频高清免费视频| 这个男人来自地球电影免费观看| 最新在线观看一区二区三区| 一区在线观看完整版| 色94色欧美一区二区| 亚洲一区二区三区欧美精品| 夜夜爽天天搞| 精品一区二区三区四区五区乱码| 人人妻人人添人人爽欧美一区卜| 777久久人妻少妇嫩草av网站| 精品国产超薄肉色丝袜足j| 日韩欧美国产一区二区入口| avwww免费| 国产一区二区三区综合在线观看| 侵犯人妻中文字幕一二三四区| 国产成人av激情在线播放| 国产精品久久久久成人av| 极品少妇高潮喷水抽搐| 两个人免费观看高清视频| 午夜福利免费观看在线| 在线观看免费午夜福利视频| 熟女少妇亚洲综合色aaa.| 欧美最黄视频在线播放免费 | 国产欧美日韩一区二区精品| 国产成人免费观看mmmm| xxxhd国产人妻xxx| 真人做人爱边吃奶动态| 在线观看免费高清a一片| 精品福利观看| 免费在线观看影片大全网站| 亚洲精品中文字幕在线视频| 美女视频免费永久观看网站| 人人妻,人人澡人人爽秒播| 丁香欧美五月| 久久国产精品男人的天堂亚洲| 黄色视频不卡| 母亲3免费完整高清在线观看| 99热国产这里只有精品6| 大码成人一级视频| 久久国产精品人妻蜜桃| 欧美激情久久久久久爽电影 | 国产色视频综合| 国产熟女午夜一区二区三区| 国产成人欧美在线观看 | 在线观看免费视频网站a站| 老司机影院毛片| 欧美精品一区二区免费开放| 制服诱惑二区| 精品午夜福利视频在线观看一区| 亚洲第一欧美日韩一区二区三区| 亚洲熟妇熟女久久| 国产亚洲精品久久久久5区| av网站免费在线观看视频| 后天国语完整版免费观看| 99香蕉大伊视频| 少妇的丰满在线观看| 久久久久久久久免费视频了| 黑人操中国人逼视频| 久久精品亚洲精品国产色婷小说| 欧美 亚洲 国产 日韩一| 9191精品国产免费久久| 九色亚洲精品在线播放| 精品国产亚洲在线| 国产精品影院久久| 18禁观看日本| 中文字幕人妻丝袜一区二区| 欧美激情高清一区二区三区| 黄色片一级片一级黄色片| 久久影院123| 国产成人免费观看mmmm| 高清黄色对白视频在线免费看| 久99久视频精品免费| 国产欧美日韩一区二区精品| 无人区码免费观看不卡| 夜夜夜夜夜久久久久| 在线观看免费高清a一片| 又紧又爽又黄一区二区| 欧美成狂野欧美在线观看| 精品国产亚洲在线| 国产精品国产av在线观看| 母亲3免费完整高清在线观看| 免费在线观看黄色视频的| 看免费av毛片| 婷婷成人精品国产| 精品乱码久久久久久99久播| 亚洲aⅴ乱码一区二区在线播放 | 色老头精品视频在线观看| 在线观看66精品国产| 久久久国产成人免费| 极品人妻少妇av视频| 丝袜美腿诱惑在线| 18禁国产床啪视频网站| 国产精品九九99| 村上凉子中文字幕在线| 村上凉子中文字幕在线| 亚洲精品粉嫩美女一区| 人人澡人人妻人| 日日夜夜操网爽| 又黄又粗又硬又大视频| 国产成人系列免费观看| 国产成人精品久久二区二区91| 操美女的视频在线观看| 9191精品国产免费久久| 日韩欧美国产一区二区入口| 久久久国产精品麻豆| 中文字幕人妻熟女乱码| 亚洲欧美日韩高清在线视频| 黄色成人免费大全| av欧美777| 黄频高清免费视频| 成人三级做爰电影| 欧美日韩精品网址| 一本一本久久a久久精品综合妖精| 午夜精品久久久久久毛片777| 国产成人免费无遮挡视频| 欧美成人午夜精品| 国产一区有黄有色的免费视频| 欧美日本中文国产一区发布| 亚洲九九香蕉| 交换朋友夫妻互换小说| 欧美日韩视频精品一区| 欧美国产精品一级二级三级| 久久 成人 亚洲| 欧美精品啪啪一区二区三区| 精品国产超薄肉色丝袜足j| 久久精品国产亚洲av香蕉五月 | 熟女少妇亚洲综合色aaa.| 一级毛片女人18水好多| 超碰成人久久| 国产男靠女视频免费网站| 老司机影院毛片| 欧美日韩福利视频一区二区| 亚洲精品成人av观看孕妇| 露出奶头的视频| 国产高清国产精品国产三级| 免费少妇av软件| tocl精华| 又紧又爽又黄一区二区| 日本a在线网址| 黄色a级毛片大全视频| 搡老岳熟女国产| 制服诱惑二区| 老司机在亚洲福利影院| 亚洲九九香蕉| 在线国产一区二区在线| 一本一本久久a久久精品综合妖精| 欧美丝袜亚洲另类 | 欧美久久黑人一区二区| 欧美精品高潮呻吟av久久| 国产色视频综合| 欧美人与性动交α欧美软件| 国产免费现黄频在线看| 欧美精品av麻豆av| 99re在线观看精品视频| 久久午夜亚洲精品久久| 久热爱精品视频在线9| 国产一区二区三区在线臀色熟女 | 一级毛片精品| 一a级毛片在线观看| 免费人成视频x8x8入口观看| 黄色片一级片一级黄色片| 少妇裸体淫交视频免费看高清 | 91麻豆av在线| 国产高清国产精品国产三级| 免费一级毛片在线播放高清视频 | 成人国语在线视频| 欧美大码av| 咕卡用的链子| 国产精品久久久av美女十八| 国产精品一区二区在线观看99| 一级作爱视频免费观看| 久久影院123| av电影中文网址| 精品少妇一区二区三区视频日本电影| 久久亚洲精品不卡| 亚洲精品自拍成人| 性少妇av在线| 日日爽夜夜爽网站| 国产欧美日韩一区二区精品| 亚洲情色 制服丝袜| 精品少妇一区二区三区视频日本电影| 巨乳人妻的诱惑在线观看| 黄色视频不卡| 亚洲一码二码三码区别大吗| 中文字幕色久视频| 国产人伦9x9x在线观看| 手机成人av网站| 91字幕亚洲| 99国产精品免费福利视频| 亚洲精华国产精华精| 欧美丝袜亚洲另类 | 99久久综合精品五月天人人| 女同久久另类99精品国产91| 亚洲av日韩精品久久久久久密| 久久久久视频综合| 又黄又爽又免费观看的视频| 俄罗斯特黄特色一大片| 日日爽夜夜爽网站| 一二三四社区在线视频社区8| 午夜久久久在线观看| 高潮久久久久久久久久久不卡| av线在线观看网站| 国产成人一区二区三区免费视频网站| 国产高清视频在线播放一区| 午夜精品在线福利| 侵犯人妻中文字幕一二三四区| 免费女性裸体啪啪无遮挡网站| 中文字幕av电影在线播放| av在线播放免费不卡| 99久久精品国产亚洲精品| 国产99白浆流出| 国产精品久久视频播放| 欧美精品亚洲一区二区| 亚洲aⅴ乱码一区二区在线播放 | 波多野结衣av一区二区av| 麻豆国产av国片精品| 中文亚洲av片在线观看爽 | 婷婷成人精品国产| 久久午夜亚洲精品久久| 亚洲精品中文字幕在线视频| 后天国语完整版免费观看| 麻豆乱淫一区二区| 桃红色精品国产亚洲av| 老司机在亚洲福利影院| 日韩人妻精品一区2区三区| 亚洲国产精品合色在线| 日韩欧美免费精品| 91在线观看av| 亚洲av电影在线进入| 脱女人内裤的视频| 十八禁网站免费在线| 久久国产精品人妻蜜桃| 亚洲男人天堂网一区| 正在播放国产对白刺激| 国产精品免费一区二区三区在线 | 无遮挡黄片免费观看| 老司机午夜十八禁免费视频| 99久久人妻综合| 国产欧美日韩综合在线一区二区| 女人被躁到高潮嗷嗷叫费观| 日韩有码中文字幕| 女人被躁到高潮嗷嗷叫费观| 国产欧美日韩综合在线一区二区| 成人国产一区最新在线观看| 国产成人一区二区三区免费视频网站| 亚洲伊人色综图| 操出白浆在线播放| 男人舔女人的私密视频| 精品久久久久久电影网| 亚洲精品在线观看二区| 国产区一区二久久| 天天影视国产精品| 高清av免费在线| 搡老岳熟女国产| 可以免费在线观看a视频的电影网站| 麻豆乱淫一区二区| 欧美日韩中文字幕国产精品一区二区三区 | 热99re8久久精品国产| 手机成人av网站| 久久精品影院6| 好男人在线观看高清免费视频| 一区二区三区国产精品乱码| 国内精品美女久久久久久| 十八禁网站免费在线| 国产淫片久久久久久久久 | 色综合欧美亚洲国产小说| 日韩中文字幕欧美一区二区| 中文字幕熟女人妻在线| 精品99又大又爽又粗少妇毛片 | 淫妇啪啪啪对白视频| 村上凉子中文字幕在线| 俄罗斯特黄特色一大片| xxx96com| 色噜噜av男人的天堂激情| 嫁个100分男人电影在线观看| 国产久久久一区二区三区| 男女那种视频在线观看| 色在线成人网| 国产麻豆成人av免费视频| 网址你懂的国产日韩在线| 色老头精品视频在线观看| 最近视频中文字幕2019在线8| 高潮久久久久久久久久久不卡| 搡老岳熟女国产| 久9热在线精品视频| 欧美日韩综合久久久久久 | 毛片女人毛片| 成人高潮视频无遮挡免费网站| 天堂动漫精品| 人妻久久中文字幕网| 97碰自拍视频| 亚洲最大成人手机在线| 国内精品一区二区在线观看| 欧美日韩中文字幕国产精品一区二区三区| 亚洲第一电影网av| 午夜福利18| 亚洲国产色片| 婷婷精品国产亚洲av| 欧美日韩中文字幕国产精品一区二区三区| 99久久综合精品五月天人人| 久久精品人妻少妇| 日韩欧美精品免费久久 | 久久久久国产精品人妻aⅴ院| 九色成人免费人妻av| 婷婷丁香在线五月| 亚洲一区二区三区色噜噜| 一个人观看的视频www高清免费观看| 久久久久久久久中文| 日本精品一区二区三区蜜桃| 人妻久久中文字幕网| 在线观看av片永久免费下载| 亚洲精华国产精华精| 国产亚洲精品一区二区www| www.999成人在线观看| 很黄的视频免费| 亚洲色图av天堂| 亚洲av电影在线进入| 一级毛片高清免费大全| 国产成年人精品一区二区| 日韩欧美在线二视频| 最好的美女福利视频网| 俺也久久电影网| 变态另类成人亚洲欧美熟女| 久久精品91无色码中文字幕| 两个人看的免费小视频| 美女大奶头视频| 中文字幕精品亚洲无线码一区| 亚洲电影在线观看av| 日韩欧美国产一区二区入口| 久久中文看片网| 亚洲欧美日韩高清专用| 99国产精品一区二区三区| 国产一区二区亚洲精品在线观看| 国内久久婷婷六月综合欲色啪| 日韩人妻高清精品专区| 最近最新免费中文字幕在线| 亚洲va日本ⅴa欧美va伊人久久| 欧美性猛交黑人性爽| 日日摸夜夜添夜夜添小说| 一级a爱片免费观看的视频| 少妇人妻精品综合一区二区 | 偷拍熟女少妇极品色| 国产精品电影一区二区三区| 久久久久久久久大av| 人妻久久中文字幕网| 久久人人精品亚洲av| 在线免费观看不下载黄p国产 | 亚洲国产精品久久男人天堂| 成人三级黄色视频| 狂野欧美白嫩少妇大欣赏| 小说图片视频综合网站| 久久6这里有精品| 国产精品野战在线观看| 欧美bdsm另类| 亚洲五月婷婷丁香| 欧美绝顶高潮抽搐喷水| 欧美一区二区国产精品久久精品| 精品一区二区三区av网在线观看| 婷婷六月久久综合丁香| 亚洲欧美精品综合久久99| 男女午夜视频在线观看| 久9热在线精品视频| 动漫黄色视频在线观看| 18美女黄网站色大片免费观看| 国产毛片a区久久久久| 欧美日韩乱码在线| 免费av观看视频| 亚洲av成人av| 久久久久久久亚洲中文字幕 | 亚洲国产欧洲综合997久久,| 国产高清videossex| 草草在线视频免费看| 国产单亲对白刺激| 国产野战对白在线观看| 国产视频内射| 美女大奶头视频| 国产成人啪精品午夜网站| 久久久成人免费电影| 高清在线国产一区| 久久99热这里只有精品18| 日韩大尺度精品在线看网址| 日韩欧美免费精品| 国产亚洲精品久久久久久毛片| bbb黄色大片| 搡老岳熟女国产| 久久人妻av系列| 免费av毛片视频| 午夜视频国产福利| 亚洲av熟女| 亚洲国产高清在线一区二区三| 亚洲人成电影免费在线| eeuss影院久久| 精品一区二区三区av网在线观看| 偷拍熟女少妇极品色| 小说图片视频综合网站| 99热6这里只有精品| 久久亚洲精品不卡| 国产精品野战在线观看| 欧美在线一区亚洲| 亚洲电影在线观看av| 国产亚洲av嫩草精品影院| 午夜福利在线观看免费完整高清在 | 91在线精品国自产拍蜜月 | 国产精品精品国产色婷婷| 91九色精品人成在线观看| 欧美性感艳星| 无限看片的www在线观看| 少妇裸体淫交视频免费看高清| 久久精品国产亚洲av涩爱 | 中文字幕人成人乱码亚洲影| 亚洲熟妇中文字幕五十中出| 美女高潮的动态| 色视频www国产| av欧美777| 国产高清三级在线| 男人舔奶头视频| 可以在线观看的亚洲视频| 亚洲精品在线观看二区| 少妇的丰满在线观看| 搡老岳熟女国产| 色综合欧美亚洲国产小说| 美女大奶头视频| 人妻久久中文字幕网| 国产伦精品一区二区三区视频9 | 91在线观看av| a级毛片a级免费在线| 国产激情偷乱视频一区二区| 亚洲狠狠婷婷综合久久图片| 亚洲国产日韩欧美精品在线观看 | 日韩欧美 国产精品| 亚洲成人久久爱视频| 99国产极品粉嫩在线观看| 国产亚洲欧美98| 2021天堂中文幕一二区在线观| 成人性生交大片免费视频hd| 好男人电影高清在线观看| 精品久久久久久,| 99久久无色码亚洲精品果冻| 欧美精品啪啪一区二区三区| 一个人免费在线观看电影| 无限看片的www在线观看| 听说在线观看完整版免费高清| 99久久99久久久精品蜜桃| 亚洲av日韩精品久久久久久密| 真人做人爱边吃奶动态| 噜噜噜噜噜久久久久久91| 日韩欧美在线乱码| 在线观看美女被高潮喷水网站 | h日本视频在线播放| 国产精品久久久久久精品电影| 色av中文字幕| 九九久久精品国产亚洲av麻豆| 少妇的丰满在线观看| 精品人妻1区二区| 欧美日韩一级在线毛片| 999久久久精品免费观看国产| 青草久久国产| 国内精品美女久久久久久| 成年人黄色毛片网站| 不卡一级毛片| 好男人在线观看高清免费视频| 亚洲 欧美 日韩 在线 免费| 欧美激情在线99| 99视频精品全部免费 在线| 亚洲精华国产精华精| 亚洲中文字幕一区二区三区有码在线看| 非洲黑人性xxxx精品又粗又长| 一a级毛片在线观看| 免费在线观看日本一区| 一二三四社区在线视频社区8| 少妇的逼水好多| 免费电影在线观看免费观看| 免费一级毛片在线播放高清视频| 国产黄色小视频在线观看| 麻豆成人av在线观看| tocl精华| 最好的美女福利视频网| 久久伊人香网站| 97超视频在线观看视频| 99精品在免费线老司机午夜| 有码 亚洲区| 好看av亚洲va欧美ⅴa在| 午夜亚洲福利在线播放| 中文亚洲av片在线观看爽| 在线观看日韩欧美| 亚洲国产精品999在线| 男女那种视频在线观看| 综合色av麻豆| 欧美日韩瑟瑟在线播放| 亚洲av免费高清在线观看| 日本与韩国留学比较| 久久伊人香网站| 麻豆国产97在线/欧美| 亚洲国产日韩欧美精品在线观看 |