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

    可壓縮壁湍流熱力學量統(tǒng)計特性分析

    2023-06-28 00:48:18傅亞陸袁先旭劉朋欣余明
    航空學報 2023年9期
    關鍵詞:壁溫馬赫數(shù)法向

    傅亞陸,袁先旭,劉朋欣,余明

    空氣動力學國家重點實驗室,綿陽 621000

    可壓縮壁湍流是航空航天領域中廣泛存在的流動現(xiàn)象,深入理解并掌握其統(tǒng)計規(guī)律具有重要的科學意義和工程價值[1-2]??蓧嚎s湍流比不可壓縮湍流更加復雜,不僅包含各種形態(tài)的渦結構,還包含膨脹和壓縮等復雜波系結構。在可壓縮湍流中,動能-內能通過壓力做功互相轉化,速度、壓力、密度、溫度和熵等運動學/動力學量和熱力學量脈動存在強相互耦合[3],因此研究熱力學平均量和脈動量對于理解湍流中動能-內能之間的轉換、建立壁模型等方面具有重要意義[4]。

    隨著計算性能的飛速提升和數(shù)值方法的發(fā)展和進步[5-6],高精度直接數(shù)值模擬(Direct Numerical Simulation,DNS)方法已成為研究湍流機理的有效手段[7-9]。目前國內外學者已經(jīng)對平板邊界層[10-13]、槽道[14-17]和圓管[18]等規(guī)范壁湍流開展了廣泛的DNS 研究,并取得較大進展。然而目前大多數(shù)研究僅針對速度、溫度和壓力脈動的統(tǒng)計特性,而對密度、熵以及各熱力學量脈動間的相關性的研究還十分有限。

    Coleman 等[14]對馬赫數(shù)為1.5 和3.0 的可壓縮等溫壁槽道湍流開展了直接數(shù)值模擬,分析了溫度、密度和壓力脈動的均方根、流向和展向兩點相關量、聯(lián)合概率密度分布等統(tǒng)計特性,結論表明壓縮性對熱力學量脈動的影響不能忽略,等溫邊界條件顯著影響了壁面法線方向的平均密度和溫度梯度,最高溫度/最低密度出現(xiàn)在槽道中心線處;隨馬赫數(shù)增大,壁面?zhèn)鳠崧孰S之增大,近壁條帶結構相干性增強。2011 年,Wei 和Pollard[19]研究了馬赫數(shù)為0.2、0.7 和1.5 條件下,速度脈動和溫度、密度、壓力等熱力學量的脈動均方根、平坦度、偏斜度隨中心線馬赫數(shù)變化的規(guī)律,并進一步分析了馬赫數(shù)對密度和壓力脈動的相關系數(shù)以及脈動壓力梯度和脈動渦量梯度相關系數(shù)的影響規(guī)律,發(fā)現(xiàn)二階統(tǒng)計量在內區(qū)和馬赫數(shù)相關,但在外區(qū)幾乎不受馬赫數(shù)影響,且近壁大尺度結構受馬赫數(shù)影響明顯[20]。

    2013 年,Donzis 和Jagannathan[4]對DNS 數(shù)據(jù)進行分析,研究了可壓縮均勻各向同性湍流的壓力、密度和溫度脈動的統(tǒng)計量,采用小擾動分析研究了脈動量、互相關和譜的標度規(guī)律,發(fā)現(xiàn)壓力和密度脈動的概率密度函數(shù)在低湍流馬赫數(shù)時呈負偏態(tài),與不可壓縮結果一致,而在高湍流馬赫數(shù)時呈正偏態(tài);雷諾數(shù)效應對單點統(tǒng)計量來說并不顯著。Gerolymos 和Vallet[21]研究了低雷諾數(shù)可壓縮槽道湍流中的熱力學量脈動,并推導了熱力學通量的輸運方程,研究了統(tǒng)計量、相關系數(shù)、平衡分析,發(fā)現(xiàn)溫度、密度、壓力的歸一化脈動均方根在槽道任意位置量級均相同,在低馬赫數(shù)條件下亦是如此。對熱力學脈動量、互相關進行研究發(fā)現(xiàn),近壁區(qū)壓力和溫度相關性較弱,但在槽道中心線處增大;整個槽道內,壓力和密度脈動、壓力和溫度脈動之間呈正相關,近壁區(qū)密度和溫度呈負相關。除此之外,熵脈動和溫度脈動存在很強的正相關,密度脈動存在很強的負相關,在法向范圍y+∈(1, 100)范圍內(其中y+為黏性尺度下的法向坐標),熵脈動與溫度脈動的相關系數(shù)接近1,表明了二者的強相關性。速度-熱力學量脈動相關性和中心線馬赫數(shù)幾乎無關。Gerolymos 和Vallet[22]進一步研究了熱力學量脈動相關系數(shù)之間的關系,結果表明,基本熱力學量的二階矩以及無量綱脈動熵與特征馬赫數(shù)呈正相關,且這些無量綱量具有相同的量級。

    以上文獻調研表明,目前對可壓縮壁湍流中熱力學量脈動的系統(tǒng)研究依然較少,壓縮性對熱力學量脈動的影響尚不明晰。因此,本文通過分析不同槽道中心線馬赫數(shù)MaC=0.88~6.15 以及不同壁溫比(壁溫和恢復溫度之比)條件下的DNS 數(shù)據(jù),研究可壓縮槽道湍流中熱力學量脈動統(tǒng)計特性,進而給出壓力、密度、溫度、熵脈動的統(tǒng)計矩,以及各物理量之間的相關系數(shù)。

    1 物理模型和數(shù)值方法

    本文研究的物理模型為可壓縮槽道湍流,如圖1 所示,其中Ub、ρ0、T0分別為來流速度、來流密度、來流溫度,Tw為壁面溫度,u、v、w為速度分量。數(shù)據(jù)來源本文作者在前期建立的DNS 數(shù)據(jù)庫,具體計算設置和結果的對比驗證請參考文獻[1],在此簡要介紹如下。

    圖1 物理模型示意圖Fig. 1 Schematic diagram of physical model

    流動的控制方程為完全氣體的可壓縮Navier-Stokes 方程,用張量形式可表示為

    式中:t為時間;ρ、p、E分別為密度、壓力和單位質量總能;ui(i=1, 2, 3)為xi方向的速度分量;τij為分子黏性應力;f1為體積力項;?為熱匯項;qj為分子熱傳導引起的熱流。

    分子黏性應力滿足如下線性本構關系:

    式中:μ為分子黏性系數(shù),由Sutherland 公式給出。

    分子熱傳導引起的熱流qj為

    式中:T表示溫度;導熱系數(shù)k=μcp/Pr,cp為定壓比熱,Pr=0.7 為普朗特數(shù)。

    熱力學量滿足如下完全氣體狀態(tài)方程:

    式中:R為氣體常數(shù);cv為定容比熱;γ=1.4 為比熱比;e為內能;s為熵。此外,動量方程中增加了均勻分布的體積力項f1以抵消黏性力做功,在能量方程中增加了熱匯項?以抵消黏性耗散產(chǎn)生的熱流,保持體平均溫度不變。

    體平均溫度Tb定義為

    式中:α為熱流控制參數(shù);Tr代表來流馬赫數(shù)Ma0時邊界層恢復溫度。

    其中:r為恢復系數(shù)。

    對于本文中的可壓縮槽道湍流,流向和展向均為周期邊界條件,上下壁面處速度滿足無滑移和不可穿透條件,溫度滿足等溫條件。

    本文共考察7 組不同的算例,物理參數(shù)如表1所示。其中平均馬赫數(shù)Mab、槽道中心線處馬赫數(shù)MaC和摩擦雷諾數(shù)Reτ的定義為

    表1 物理參數(shù)Table 1 Physical parameters

    式中:Ub為平均速度;UC為槽道中心線處速度;TC為槽道中心線處溫度;h為半槽道高度;ρw為壁面密度;μw為黏性系數(shù);uτ為摩擦速度。

    基于壁面總切應力τw(黏性和雷諾切應力之和)、壁面密度ρw和黏性系數(shù)μw,uτ可定義為

    在本文的算例中,各組算例的摩擦雷諾數(shù)Reτ均接近500,平均馬赫數(shù)范圍為0.77~4.79,槽道中心線馬赫數(shù)范圍為0.88~6.15,包含了從亞聲速到高超聲速的參數(shù)范圍;A1~A5 為近似絕熱壁條件的流動,C1~C2 為高超聲速冷壁流動。

    直接數(shù)值模擬采用Li 等開發(fā)的基于有限差分的OpenCFD-SC 開源程序[23]。其中對流項在Steger-Warming 矢通量分裂后采用七階迎風格式離散,黏性項采用六階中心差分格式離散,時間推進采用3 階TVD 龍格-庫塔格式。

    計算網(wǎng)格劃分如圖2 所示,展向z和流向x采用均勻網(wǎng)格,法向y采用雙曲正切函數(shù)進行拉伸(在壁面附近密集)。計算參數(shù)如表2 所示。(Lx,Ly,Lz)為計算域尺寸,(Nx,Ny,Nz)為網(wǎng)格數(shù)。需要說明的是,算例C2 中Lx=6πh,這是因為在較低壁溫下流向尺度更大,過小的計算域長度會引入非物理周期性。網(wǎng)格分辨率的選取已在文獻[1]中做了說明,結果表明黏性尺度下的網(wǎng)格間距足以分辨湍流小尺度脈動。

    表2 計算參數(shù)Table 2 Simulation parameters

    圖2 網(wǎng)格劃分示意圖Fig. 2 Schematic diagram of grid

    為了更好地理解物理量的脈動,本文采用系綜平均統(tǒng)計方法,在流向、展向和時間方向作平均。對于流動物理量φ,其系綜平均為,對應的湍流脈動為φ′,即

    脈動均方根為φ′rms。采用摩擦速度uτ、壁面總切應力τw、壁面密度ρw無量綱化的變量用上標“+”表示。

    2 熱力學量脈動的瞬時場結構

    圖3~圖6 分別為算例A1、A5 和C1 的溫度、熵、密度和壓力脈動的瞬時分布,分別對應了低馬赫數(shù)和高馬赫數(shù)絕熱壁流動,以及高馬赫數(shù)冷壁流動。為了更清晰地展示流場的三維結構,本文選取了與壁面垂直的沿流向和展向分布的x-y和y-z平面和近壁處y+=30 的x-z平面上的流場分布。

    圖3 溫度脈動T′/T0 的瞬時分布Fig. 3 Instantaneous structures of temperature fluctuations T′/T0

    低馬赫數(shù)流動中,近壁x-z平面上溫度脈動T′(見圖3)主要表現(xiàn)為沿流向分布的條帶結構,與流向速度脈動類似。隨馬赫數(shù)的增加和壁溫的降低,溫度脈動中沿流向正負交替分布的行波包絡結構逐漸顯現(xiàn),與速度脈動類似,這是壓縮性效應在流場相干結構中的直接體現(xiàn)。此外,在近壁處的剪切層所引起的較強的湍流耗散會使局部溫度升高,對應薄層的強溫度脈動。隨著法向高度的增加,相干結構的尺度逐漸增加(y-z平面)。在外區(qū)存在大尺度結構(x-y平面),其中附著在壁面上的結構與流向夾角約為45°,而非附著結構則更傾向于各向同性的分布。

    熵脈動s′(見圖4)在x-z平面上以大尺度的“斑塊”結構為主,沒有明顯地表現(xiàn)出條帶結構,表明在近壁區(qū)溫度脈動與熵脈動的相關性較弱。在高馬赫數(shù)流動中,近壁處的熵脈動同樣存在沿流向正負交替的行波包絡結構,表明壓縮性效應同樣體現(xiàn)在熵脈動中,并可能使溫度與熵脈動相關性增強。在外區(qū),熵脈動的相干結構與溫度脈動相干結構分布類似,表明二者相關性較強。

    圖4 熵脈動s′/s0 的瞬時分布Fig. 4 Instantaneous structures of entropy fluctuations s′/s0

    密度脈動ρ′(見圖5)在外區(qū)與溫度和熵脈動的相干結構類似,但符號相反。在近壁處,密度脈動強度較高,但是沿流向拉長的條帶結構并不明顯,而是表現(xiàn)為流向和展向特征尺度相差不大的在平面上“各向同性”的相干結構,與壓力脈動(見圖6)類似。根據(jù)狀態(tài)方程,等溫壁面附近壓力和密度脈動是強相關的,因此相干結構分布類似這一現(xiàn)象也是合理的。在外區(qū),壓力脈動與其他熱力學量的脈動均存在明顯差別,因為壓力不僅是熱力學量,也是動力學量,其相干結構同樣受到速度脈動和渦量脈動的影響。

    圖5 密度脈動ρ′/ρ0 的瞬時分布Fig. 5 Instantaneous structures of density fluctuations ρ′/ρ0

    圖6 壓力脈動p′/p0 的瞬時分布Fig. 6 Instantaneous structures of pressure fluctuations p′/p0

    整體而言,溫度脈動T′、熵脈動s′和密度脈動ρ′在外區(qū)的相干結構分布類似;在近壁區(qū)的低馬赫數(shù)流動中,上述熱力學量脈動結構的分布形式差別較大,而在高馬赫數(shù)流動中,行波包絡形式的結構均較為明顯,表明壓縮性效應的增強對熱力學量脈動及其互相關性具有較大影響,在流動建模中需考慮這一因素。

    3 結果與討論

    本節(jié)重點討論熱力學的平均量和脈動量的單點統(tǒng)計矩,包括脈動均方根、偏斜度、平坦度,以及各熱力學量之間及其與速度脈動的互相關系數(shù)。

    3. 1 平均分布

    圖7 為熱力學量的平均值沿法向的分布。圖7(a)中的溫度分布Tˉ(y)由中心線處平均溫度無量綱化。由于壁溫效應主要體現(xiàn)在近壁處,圖7 中的橫坐標采用了對數(shù)坐標。對于算例A1~A5,壁面處溫度梯度很低,在y/h<0.02 范圍內平均溫度幾乎不變,而在該法向高度以上,平均溫度單調降低。算例C1 和C2 則在近壁區(qū)存在一個溫度的極大值點,在靠近壁面處溫度隨法向高度的增加而升高,遠離壁面處的變化趨勢相反。需要說明的是,即使高馬赫數(shù)流動中流場中平均溫度的最大值與槽道中心線處相差很大,這一溫度梯度效應對統(tǒng)計量影響也并不顯著。下文的研究將進一步揭示,平均溫度梯度對湍流統(tǒng)計特性的影響主要體現(xiàn)在近壁區(qū)。圖7(b)為平均密度沿法向高度的變化,由體平均密度無量綱化。由于流場中沒有明顯的壓力梯度,平均壓力沿法向高度的變化很小,因此平均密度和平均溫度幾乎呈反比。

    圖7 熱力學平均量沿法向分布Fig. 7 Profiles of mean thermodynamic variables

    圖7(c)為平均熵sˉ(y)沿法向高度的變化。由于熵的絕對值沒有意義,本文將壁面處的平均熵取為0,即考察平均熵與壁面平均熵之差。對于絕熱壁算例A1~A5,壁面處的平均熵最高,隨法向高度增加而單調遞減;冷壁算例C1 和C2 則先增加后減小,在近壁處存在一個極大值點。平均熵的分布與平均溫度分布在趨勢上具有一致性。

    由于流場中沒有明顯的壓力梯度,平均壓力沿法向的變化很小。由平均動量方程可知:

    其中,黏性項的影響很小,因此壓力沿法向高度的變化主要源于雷諾正應力的法向分量。圖7(d)所示為平均壓力與壁面平均壓力差的法向分布,由壁面平均切應力無量綱化。該結果與文獻[1]中雷諾應力法向正應力分量的分布和變化趨勢相一致。該結果也表明,對于無明顯壓力梯度的流動,平均壓力的分布主要取決于速度脈動,而非熱力學量。

    3. 2 脈動均方根

    對于流動物理量的脈動量φ′,其均方根(RMS)的定義為

    溫度、密度、熵和壓力脈動的均方根如圖8 所示,分別由壁面平均溫度Tw、密度ρw、熵s0、壓力pw歸一化。

    圖8 歸一化的均方根熱力學量脈動沿壁面法向分布Fig. 8 Profiles of RMS thermodynamic variables fluctuations

    溫度脈動均方根T′RMS/Tw如圖8(a)所示。由于壁面處等溫條件的限制,其壁面脈動強度為0。隨法向高度的增加,溫度脈動在近壁處達到峰值,并在外區(qū)y>0.2h范圍內單調遞減。峰值所對應的法向位置隨馬赫數(shù)的增加而逐漸遠離壁面。此外,在冷壁條件下,近壁處的溫度脈動表現(xiàn)出“雙峰值”的特點,這是因為平均溫度分布是非單調變化的。對于算例C2,其峰值位置約位于y=0.05h處,該峰值上下的平均溫度梯度較大,加之緩沖區(qū)的強湍流脈動,使得在峰值上下的近壁處出現(xiàn)了這一雙峰值的特點。這一現(xiàn)象在冷壁條件湍流邊界層中是廣泛存在的[17,24]。隨中心線馬赫數(shù)MaC的增加,外區(qū)y>0.2h范圍內的溫度脈動強度顯著增強。由于馬赫數(shù)是反映流體動能和內能量級之比的參數(shù),不難推斷,隨馬赫數(shù)的增加,黏性對湍動能的耗散所引起的內能脈動增強,且壓力對體膨脹運動做功增強,使得溫度脈動隨馬赫數(shù)增強。

    密度脈動均方根ρ′RMS/ρw如圖8(b)所示。與溫度不同,密度在壁面處不需直接滿足第一類條件,因此壁面處的密度脈動不為0。隨中心線馬赫數(shù)MaC增加和壁溫的降低,壁面處密度脈動強度顯著增加。隨法向高度的增加,密度脈動強度先減小后增大,其最小值在近壁區(qū),且隨馬赫數(shù)的增加和壁溫的降低而逐漸外移。在外區(qū)y>0.2h范圍內,密度脈動強度隨法向距離的增加而增加。與壁面脈動不同,冷壁條件的算例C2 在外區(qū)比算例C1 更低。從物理意義上,流體微團內的密度脈動對應流體微團的體膨脹運動,因此密度脈動強弱直接反映了湍流脈動中壓縮性效應的強弱。因此不難推斷,在近壁處,壓縮性效應隨馬赫數(shù)的增加和壁溫的降低而增強,這與Yu和Xu[17]對真實壓縮性效應的研究結論是一致的;在外區(qū),壓縮性的強弱同樣與馬赫數(shù)呈正相關,而隨壁溫的降低,壓縮性效應更弱。

    熵的脈動強度如圖8(c)所示。與溫度脈動類似,熵是不直接與運動學和動力學相關的物理量。對于完全氣體,熵值變化主要與2 個因素有關:一是由黏性耗散引起動能向內能的轉換這一非可逆過程引起;二是由壓力對流體微團的強壓縮作用引起,類似于激波對流體的非等熵壓縮作用。在壁面附近,雖然黏性耗散作用很強,但是黏性耗散的脈動(反映在溫度脈動中)強度不高;而壓縮性效應隨馬赫數(shù)的增加和壁溫的降低顯著增強,因此熵脈動強度增加。在外區(qū),熵脈動隨法向高度、馬赫數(shù)和壁溫的變化與密度和溫度類似,在此不再贅述。

    壓力脈動強度如圖8(d)所示。由壁面平均壓力pw無量綱化的壓力脈動強度隨馬赫數(shù)和壁溫的變化趨勢較為復雜。對于最低馬赫數(shù)的算例A1,壓力脈動強度沿法向高度的分布先增大后減小,與不可壓縮湍流的結果一致,具體驗證請參考文獻[25]。對于算例A2~A5,隨馬赫數(shù)的增加,壁面處壓力降低,而接近槽道中心線處壓力升高。特別是對于算例A5,在外區(qū)y>0.4h至槽道中心線處壓力脈動強度幾乎不變。隨壁溫的降低,壓力脈動強度進一步減小,但在外區(qū)壓力脈動強度幾乎為常數(shù)這一現(xiàn)象沒有明顯改變。與溫度、密度和熵等物理量不同,壓力既是熱力學量,也是動力學量,是建立動量和能量之間關系的橋梁。Yu 等[26]的研究對壓力脈動泊松方程進行了分解,發(fā)現(xiàn)由剪切運動所引起的壓力脈動與馬赫數(shù)無關,而由體膨脹運動所引起的壓力脈動則與馬赫數(shù)有較大的關聯(lián)。在黏性尺度無量綱化后,壁面和槽道中心線處壓力脈動與馬赫數(shù)呈正相關。本文以外尺度無量綱化,因此這一趨勢沒有得到直接體現(xiàn),但是不難總結,壓縮性效應對壓力脈動的影響十分顯著。

    3. 3 概率密度分布

    概率密度函數(shù)(Probability Density Function,PDF)能夠更為完整地反映各物理量在各強度脈動的分布[4]。本文對PDF 預乘了脈動的平方,該曲線與橫軸所圍成的面積即代表了脈動的方差,直接反映了不同強度脈動對均方根的貢獻。

    圖9 為y+=30 的預乘PDF 分布。其中溫度、熵、密度和壓力脈動均表現(xiàn)為較強的正負不對稱性。對于溫度和熵脈動,正值脈動分布范圍較小,在約為1.2T′RMS和1.2s′RMS處達到峰值,且高于2.4T′RMS和2.4s′RMS的脈動對均方根的貢獻可以忽略不計;與之相比,負值脈動分布范圍更廣,在約為-1.8T′RMS和-1.8s′RMS處達到峰值,在脈動值低于-4.2T′RMS和-4.2s′RMS范圍內的極端事件對均方根的貢獻較小。整體而言,溫度和熵脈動的預乘PDF 分布隨馬赫數(shù)變化不大,而對壁溫變化更為敏感。特別是冷壁算例C2,其強負值脈動的極端事件更為頻繁,對脈動均方根的貢獻更大。

    圖9 y+=30處熱力學量脈動的預乘概率密度函數(shù)分布Fig. 9 Pre-multiplied probability density functions of thermodynamic fluctuations at y+=30

    密度脈動的預乘PDF(圖9(c))則與溫度和熵脈動幾乎呈現(xiàn)關于縱坐標軸線對稱的分布:負值脈動分布范圍較小,但預乘PDF 值更高;正值脈動分布范圍較大,極端事件對脈動均方根的貢獻更高。與溫度和熵不同的是,密度脈動的PDF分布對馬赫數(shù)和壁溫均很敏感。隨馬赫數(shù)的增加,負值脈動的分布范圍變窄,預乘PDF 的峰值升高,表明極端事件發(fā)生的概率更低,對脈動均方根的貢獻更小,而降低壁溫則起到了相反的作用。正值脈動的預乘PDF 分布則隨馬赫數(shù)的升高和壁溫的降低而表現(xiàn)為分布范圍更廣、峰值更低的特點,即極端事件頻率更高、貢獻更大。3.2 節(jié)對密度脈動強度的討論表明,隨馬赫數(shù)的升高和壁溫的降低,近壁處的壓縮性效應更強。圖9(c)中的PDF 分布同樣證實了這一結論的正確性——其預乘PDF 值升高表明流場中密度脈動不僅是速度脈動對平均密度梯度的輸運引起,強壓縮事件所引起的局部流體密度升高這一效應同樣不可忽略。

    壓力脈動p′的預乘概率密度分布同樣表現(xiàn)為正負不對稱的特點。Pumir[27]的研究表明,在不可壓湍流中,壓力脈動的PDF 分布fp′呈負偏斜,且正方向尾部接近高斯分布。對于本文的算例,隨著馬赫數(shù)的升高和壁溫的降低,負值脈動的預乘PDF 分布變窄,且正、負壓力脈動趨向于對稱。Donzis 和Jagannathan[4]的研究表明,可壓縮湍流中正壓力脈動的概率增大,這與本文中所給出的結果定性上是一致的,表明可壓縮與不可壓縮湍流中的壓力脈動統(tǒng)計特性存在較大區(qū)別。

    圖10 為外區(qū)y=0.2h處的預乘PDF 分布。與近壁區(qū)y+=30 處的結果相比,溫度脈動和熵脈動在正值和負值區(qū)域上的分布對稱性更強,且隨馬赫數(shù)和壁溫的變化不大。密度脈動和壓力脈動則表現(xiàn)處與y+=30 處一致的非對稱性,且隨馬赫數(shù)和壁溫的變化趨勢與y+=30 處相同,差別在于極端事件對脈動強度的貢獻更小。

    圖10 y=0.2h處熱力學量脈動的預乘概率密度函數(shù)分布Fig. 10 Pre-multiplied probability density functions of thermodynamic fluctuations at y=0.2h

    3. 4 偏斜度和平坦度

    本節(jié)討論熱力學量的三階和四階統(tǒng)計矩,即偏斜度S(φ′)和平坦度F(φ′)(其中φ′為任一流動物理量),其定義為

    式中:偏斜度S(φ′)反映了隨機變量分布的不對稱性,完全對稱分布的三階矩為0;F(φ′)則反映了隨機變量的間歇性,其中高斯分布的間歇因子為3.0,高于或低于該值表示該隨機變量偏離高斯分布。

    圖11 為各算例熱力學量脈動的偏斜度沿壁面法向分布。整體而言,溫度、密度、熵、壓力脈動幾乎在全法向高度均不滿足正態(tài)分布。其中溫度脈動的偏斜度隨馬赫數(shù)和壁溫變化不大。在y<0.2h范圍內變化較大,其偏斜度為負,幅值較大,表明近壁處的溫度脈動為負的概率更高。在y>0.2h的外區(qū),偏斜度從負值逐漸變?yōu)檎担砻髟谕鈪^(qū)溫度脈動為正的概率更高。熵脈動的偏斜度的變化趨勢與溫度脈動類似,但是S(s′)隨馬赫數(shù)和壁溫變化更為敏感,其中算例C1 和C2 在全法向高度內均為負。

    圖11 偏斜度沿壁面法向y/h 分布Fig.11 Distributions of skewness factors, plotted against y/h

    密度脈動的偏斜度隨馬赫數(shù)的增加而升高,且隨壁溫的降低而降低。對于較低馬赫數(shù)的算例A1 和A2,密度脈動在靠近壁面處為正,而在遠離壁面處為負,這一趨勢與溫度脈動相反。而在較高馬赫數(shù)的算例中,密度脈動在全法向高度范圍內均為正。在3.3 節(jié)的討論中已經(jīng)指出,由于在較高馬赫數(shù)流動中存在較強的壓縮事件,密度脈動的正極端事件比負極端事件發(fā)生更為頻繁且強度更高,因此隨馬赫數(shù)的增加,密度脈動的偏斜度越來越高。

    對于壓力脈動,低馬赫數(shù)算例A1 和A2 在壁面處的偏斜度幾乎為0,這與不可壓縮湍流中的結論相一致[28]。而在壁面處的其他法向位置,壓力脈動的偏斜度均為負值。隨馬赫數(shù)的增加和壁溫的降低,壁面處的偏斜因子逐漸變?yōu)檎担ㄏ蛭恢玫钠倍纫搽S馬赫數(shù)的增加和壁溫的降低而逐漸升高,由低馬赫數(shù)的負值逐漸增加為正值。這一變化規(guī)律與密度脈動偏斜度隨馬赫數(shù)和壁溫的變化規(guī)律類似,這是由于高馬赫數(shù)流動中的強壓縮事件會同時引起正密度脈動和正壓力脈動。

    圖12 為熱力學量脈動的平坦度沿壁面法向分布。高斯分布的平坦因子為3.0。平坦度越大,間歇性越強,出現(xiàn)極端事件的概率也就越大。對于溫度脈動,在靠近壁面y<0.006h范圍內,其平坦度均高于3.0,表明即使是在等溫條件的約束下,靠近壁面處的溫度脈動同樣具有高于高斯分布的間歇性。在0.01h<y<0.03h范圍內,F(xiàn)(T′)明顯升高,且在平均溫度最大值位置處達到峰值。這一現(xiàn)象是顯而易見的,因為在平均脈動最大值位置處溫度梯度為0,較弱的溫度脈動對總溫度脈動均方根的貢獻很弱,而較強的溫度脈動則貢獻很強,進而表現(xiàn)為較強的間歇性。在0.1h<y<0.5h范圍內的外區(qū),F(xiàn)(T′)幾乎為常數(shù),且略低于3.0,表明此時溫度脈動的間歇性較弱。在槽道中心線附近,F(xiàn)(T′)值略有抬升。熵脈動的平坦度變化趨勢與溫度脈動類似,在近壁y<0.006h和外區(qū)0.1h<y<0.5h范圍內F(s′) 與F(T′) 具體數(shù)值相差也不大,而在0.01h<y<0.03h范圍內的峰值明顯升高,特別是對于高馬赫數(shù)的流動。整體而言,F(xiàn)(s′) 和F(T′)對馬赫數(shù)和壁溫的變化不敏感。

    圖12 平坦度沿壁面法向y/h 分布Fig.12 Distributions of flatness factors, plotted against y/h

    密度脈動的平坦度F(ρ′)隨法向高度的變化趨勢與F(s′)和F(T′)一致,而各組算例的結果之間存在明顯差異。在近壁y<0.006h范圍內F(ρ′)隨馬赫數(shù)變化不大,而隨壁溫的降低明顯升高,表明近壁區(qū)的壓縮性效應使得密度脈動的間歇性更強。在y>0.1h范圍內,F(xiàn)(ρ′)隨馬赫數(shù)的升高和壁溫的降低均有明顯增強,表明壓縮性效應增強了密度脈動的間歇性。

    壓力脈動的平坦度F(p′)與其他熱力學量有顯著差異。對于絕熱壁流動,壓力脈動的間歇性在全法向位置的變化不大(算例A1 槽道中心線附近的結果除外),而在冷壁條件下,y<0.02h范圍內的近壁處壓力脈動間歇性明顯升高,與密度脈動的平坦度在近壁處的變化趨勢一致。

    3. 5 熱力學量的互相關系數(shù)

    對于量熱完全氣體,溫度、熵、密度和壓力之間由氣體的狀態(tài)方程相關聯(lián)。在等熵流動中若存在小擾動,上述物理量之間脈動由等熵關系確定,其相關系數(shù)也為1.0。為了研究可壓縮壁湍流中上述物理量之間的脈動關系,本文進一步對熱力學量脈動的互相關進行考察,2 個流動物理量脈動φ′1和φ′2之間的相關系數(shù)定義為[29]

    各熱力學量脈動互相關系數(shù)R(ρ′,T′),R(ρ′,s′),R(ρ′,p′) ,R(s′,T′),R(s′,p′),R(p′,T′)沿法向(外尺度y/h)分布,如圖13 所示。以下首先討論密度、壓力和溫度三者的互相關。

    圖13 熱力學量脈動間相關系數(shù)沿外尺度y/h的分布Fig. 13 Distribution of correlation coefficients between thermodynamic fluctuations, plotted against outer scaling y/h

    密度和溫度脈動的互相關系數(shù)R(ρ′,T′)(見圖13(a))在壁面附近變化較為劇烈,在靠近壁面時為正,在近壁區(qū)迅速下降為負值,在y>0.1h的區(qū)域內達到-0.8~-0.9 的強負相關,且這一相關性隨法向位置、馬赫數(shù)和壁溫的變化不大。

    密度與壓力脈動的互相關系數(shù)R(ρ′,p′) 如圖13(c)所示,二者在壁面是完全相關的,即相關系數(shù)為1.0。根據(jù)完全氣體狀態(tài)方程不難推斷,由于等溫邊界條件的影響,即T′w=0,密度脈動和壓力脈動直接滿足p′=ρ′RTˉ這一代數(shù)關系。隨法向高度的增加,二者相關系數(shù)逐漸降低,在y>0.1h的區(qū)域隨法向高度變化不大。對比不同算例,不難發(fā)現(xiàn),密度和壓力脈動的相關性隨馬赫數(shù)的增加而減弱,而隨壁溫的降低而增強。

    壓力和溫度脈動的互相關系數(shù)R(p′,T′)如圖13(f)所示。二者在近壁處相關性較強,而在y>0.1h的外區(qū)相關系數(shù)幾乎為0,表明壓力脈動對溫度脈動的作用很小。

    下面討論熵脈動與其他物理量之間的互相關。熵的控制方程相對簡單,在RANS 和LES 方程中需模型化的未知項最少[21],因此討論熵脈動具有較強的應用價值。但其物理意義不如其他熱力學量明確,因此需要借助其他熱力學量對其進行理解。但是除近壁區(qū)外,熵脈動與溫度脈動的相關系數(shù)R(s′,T′)幾乎為1.0,表明二者具有很強的互相關性,與Gerolymos[21]和Kovasznay[30]等的研究結果一致。由熵的定義可知,熵的全微分可以展開為溫度和密度的全微分之和,即

    由于溫度和密度脈動在外區(qū)相關性很強(如圖13(a)所示),因此熵脈動與溫度脈動之間同樣存在相關性。進一步說,熵脈動與密度脈動和壓力脈動的相關系數(shù)R(s′,ρ′)和R(s′,p′)與R(ρ′,T′)和R(p′,T′)分布類似——R(s′,ρ′)在近壁處變化較大,在外區(qū)約為0.9;R(s′,p′)在外區(qū)幾乎為0,表明在平均意義上,壓力脈動對熵脈動無直接作用。

    3. 6 速度脈動對熱力學量的輸運

    本節(jié)進一步討論速度場對熱力學量脈動輸運的影響,主要分析密度、溫度和熵脈動與流向和法向速度脈動的相關系數(shù),統(tǒng)計結果如圖14所示。

    圖14 熱力學量脈動與流向速度和法向速度脈動的相關系數(shù)沿y/h 的分布Fig. 14 Distribution of correlation coefficients between thermodynamic fluctuations and streamwise and wall normal velocity fluctuations, plotted against y/h

    其中速度與溫度脈動的相關系數(shù)R(u′,T′)(見圖14(a))在近壁區(qū)值為正,在外區(qū)值為負;R(v′,T′)(見圖14(b))在近壁區(qū)值為負,在外區(qū)值為正。這一現(xiàn)象與Gerolymos 和Vallet[21]的結論一致。如3.1 節(jié)所討論,這是由于近壁區(qū)平均溫度隨法向高度增加而升高,而外區(qū)平均溫度隨法向高度增加而減小所導致的。壁湍流中的速度脈動主要源于上拋和下掃事件。在近壁區(qū),壁面附近的低溫低速流體向上運動,遠離壁面的高溫高速流體向下運動,因此R(u′,T′)為正、R(v′,T′)為負。而在外區(qū),上拋和下掃事件引起面附近的高低速流體向上運動,遠離壁面的低溫高速流體向下運動,因此R(u′,T′)為負、R(v′,T′)為正。

    強雷諾比擬中指出,流向速度與溫度脈動具有強相關性,即R(u′,T′)≈1。然而大量數(shù)值模擬結果表明,二者的相關系數(shù)的絕對值一般不高于0.9,且與壁溫有關。本文的速度和溫度脈動相關系數(shù)為R(u′,T′)≈-0.7。法向速度與溫度脈動之間的相關系數(shù)的幅值R(v′,T′)與R(v′,s′)相接近,約為0.4。此外,馬赫數(shù)和壁溫對R(u′,T′)和R(v′,T′)的影響并不顯著。

    根據(jù)3.5 節(jié)的討論,溫度脈動與密度脈動呈現(xiàn)較強的負相關,而與熵脈動呈現(xiàn)較強的正相關,因此R(u′,ρ′) 和R(v′,ρ′) 與R(u′,T′) 和R(v′,T′)的變化趨勢相反,而R(u′,s′)和R(v′,s′)與之變化趨勢相同。這一結果表明,RANS 和LES 中的湍流熱通量的湍流模型可以等效地用于建立密度和熵脈動通量的湍流模型。

    4 結 論

    本文對可壓縮壁湍流直接數(shù)值模擬數(shù)據(jù)庫進行了統(tǒng)計分析,研究了不同中心線馬赫數(shù)MaC=0.88~6.15 以及不同壁溫比Tw/Tr=1.0~0.2 條件下溫度、密度壓力和熵等熱力學量脈動的相干結構和統(tǒng)計特性,得出以下主要結論:

    1) 在高馬赫數(shù)冷壁條件下,近壁溫度、密度、熵和壓力均表現(xiàn)為沿流向正負交替的行波包絡結構,這是由于壓縮性效應引起的。

    2) 溫度和熵脈動統(tǒng)計特性隨馬赫數(shù)和壁溫變化趨勢類似,主要由黏性耗散產(chǎn)熱和壓力對體膨脹運動做功所引起。其中脈動均方根溫度脈動隨馬赫數(shù)增大而增強,并且冷壁條件下,近壁處脈動表現(xiàn)出“雙峰值”特點。脈動的預乘PDF分布受馬赫數(shù)影響不大,但對壁溫變化更為敏感。偏斜度和平坦度隨馬赫數(shù)和壁溫變化不大,近壁處的負脈動概率更高,但在y>0.2h的外區(qū),正脈動概率更高。

    3) 密度脈動與壓縮性效應直接相關,隨馬赫數(shù)的增加和壁溫的降低,近壁處壓縮性效應增強,脈動均方根顯著增大;外區(qū)壓縮性強弱與馬赫數(shù)呈正相關,但與壁溫呈負相關。脈動的PDF分布對馬赫數(shù)和壁溫均很敏感。偏斜度隨馬赫數(shù)的增加而升高,隨壁溫的降低而降低。平坦度在近壁y<0.006h范圍內隨馬赫數(shù)變化不大,而隨壁溫的降低明顯升高;在y>0.1h,壓縮性效應增強了密度脈動的間歇性,其隨馬赫數(shù)的升高和壁溫的降低均有明顯增強。

    4) 溫度脈動、熵脈動和密度脈動具有較強的相關性,而這些熱力學量與壓力脈動的相關性較弱。湍流運動對溫度、熵和密度脈動的輸運作用類似。

    猜你喜歡
    壁溫馬赫數(shù)法向
    高馬赫數(shù)激波作用下單模界面的Richtmyer-Meshkov不穩(wěn)定性數(shù)值模擬
    爆炸與沖擊(2024年7期)2024-11-01 00:00:00
    落石法向恢復系數(shù)的多因素聯(lián)合影響研究
    一維非等熵可壓縮微極流體的低馬赫數(shù)極限
    載荷分布對可控擴散葉型性能的影響
    機組啟動過程中溫度壓力控制分析
    壁溫對氣化爐操作的指導
    降低鄒縣發(fā)電廠#6爐屏式過熱器管壁溫度
    低溫狀態(tài)下的材料法向發(fā)射率測量
    落石碰撞法向恢復系數(shù)的模型試驗研究
    不透明材料波段法向發(fā)射率在線測量方法
    亚洲成a人片在线一区二区| 亚洲精品色激情综合| 赤兔流量卡办理| 精品日产1卡2卡| 乱系列少妇在线播放| 最近手机中文字幕大全| 国产高清有码在线观看视频| 免费不卡的大黄色大毛片视频在线观看 | 久久亚洲精品不卡| 国产爱豆传媒在线观看| 不卡视频在线观看欧美| 色av中文字幕| 国产私拍福利视频在线观看| 国产精品一区二区性色av| 亚洲av五月六月丁香网| 在线国产一区二区在线| 欧美极品一区二区三区四区| 国产黄色视频一区二区在线观看 | 99久久精品一区二区三区| 老师上课跳d突然被开到最大视频| 最新在线观看一区二区三区| 国产精品久久视频播放| 亚洲18禁久久av| 最新中文字幕久久久久| 人妻夜夜爽99麻豆av| 婷婷精品国产亚洲av| 欧美日韩国产亚洲二区| 日本黄色片子视频| 男人和女人高潮做爰伦理| 搡女人真爽免费视频火全软件 | 伦精品一区二区三区| 亚洲精品国产av成人精品 | 成熟少妇高潮喷水视频| 在线免费十八禁| 日韩高清综合在线| 久久6这里有精品| or卡值多少钱| 国产欧美日韩精品亚洲av| 国产亚洲91精品色在线| 五月伊人婷婷丁香| 2021天堂中文幕一二区在线观| 亚洲第一区二区三区不卡| 插逼视频在线观看| 国产成人freesex在线 | 婷婷精品国产亚洲av在线| 亚洲成人精品中文字幕电影| 久久久久久九九精品二区国产| 日本黄色视频三级网站网址| 日本-黄色视频高清免费观看| 熟妇人妻久久中文字幕3abv| 久久亚洲国产成人精品v| 久久久精品大字幕| 99热这里只有是精品50| 老司机午夜福利在线观看视频| 国产精品久久久久久久久免| 国产中年淑女户外野战色| 三级毛片av免费| 国产美女午夜福利| av女优亚洲男人天堂| 亚洲av熟女| 天堂动漫精品| 自拍偷自拍亚洲精品老妇| 国产精品福利在线免费观看| 亚洲欧美日韩无卡精品| 国产精品久久久久久亚洲av鲁大| 色哟哟哟哟哟哟| 国产免费一级a男人的天堂| 亚洲第一区二区三区不卡| 搞女人的毛片| 伦精品一区二区三区| 亚洲成av人片在线播放无| 国产在线男女| 国产精品电影一区二区三区| 中文字幕久久专区| 日韩欧美国产在线观看| 国产av一区在线观看免费| 免费黄网站久久成人精品| 亚洲av不卡在线观看| 成年女人毛片免费观看观看9| 少妇人妻精品综合一区二区 | 看非洲黑人一级黄片| 不卡一级毛片| 国产爱豆传媒在线观看| ponron亚洲| 欧美又色又爽又黄视频| 两个人的视频大全免费| 亚洲欧美精品自产自拍| 国产爱豆传媒在线观看| 日本免费一区二区三区高清不卡| 午夜福利在线观看免费完整高清在 | 69av精品久久久久久| 国产精品久久久久久av不卡| 天美传媒精品一区二区| 超碰av人人做人人爽久久| 偷拍熟女少妇极品色| 一夜夜www| 我的女老师完整版在线观看| 神马国产精品三级电影在线观看| 欧美一区二区精品小视频在线| 日韩精品青青久久久久久| 国产精品久久久久久亚洲av鲁大| 久久久久久九九精品二区国产| 免费在线观看影片大全网站| 91久久精品电影网| 中文字幕av成人在线电影| 一区二区三区高清视频在线| 日韩一本色道免费dvd| 国产精品一区二区性色av| 国产高清不卡午夜福利| 国产极品精品免费视频能看的| 国产极品精品免费视频能看的| 99国产极品粉嫩在线观看| 搡老岳熟女国产| 久久精品国产自在天天线| avwww免费| 五月伊人婷婷丁香| www.色视频.com| 天堂动漫精品| 男女边吃奶边做爰视频| 成人漫画全彩无遮挡| 狂野欧美激情性xxxx在线观看| 欧美性猛交黑人性爽| a级毛片a级免费在线| 国产高清有码在线观看视频| 97热精品久久久久久| 国模一区二区三区四区视频| 日本欧美国产在线视频| 国产人妻一区二区三区在| 免费在线观看成人毛片| 秋霞在线观看毛片| 久久精品国产亚洲网站| 国产中年淑女户外野战色| 高清毛片免费看| 精品久久久久久久久久免费视频| 特级一级黄色大片| 又粗又爽又猛毛片免费看| 欧美+日韩+精品| 亚洲av五月六月丁香网| 欧美区成人在线视频| 成人亚洲欧美一区二区av| 欧美+亚洲+日韩+国产| eeuss影院久久| 免费电影在线观看免费观看| 九九热线精品视视频播放| 亚洲国产精品成人综合色| 亚洲精品国产av成人精品 | 亚洲精品一卡2卡三卡4卡5卡| 国产老妇女一区| 寂寞人妻少妇视频99o| 自拍偷自拍亚洲精品老妇| 亚洲欧美成人综合另类久久久 | 成人鲁丝片一二三区免费| 久久人人爽人人片av| 久久草成人影院| 成熟少妇高潮喷水视频| 日日摸夜夜添夜夜添小说| 国产伦在线观看视频一区| 欧美不卡视频在线免费观看| 亚洲成人久久性| 欧美bdsm另类| 99久久精品热视频| 久久综合国产亚洲精品| 精品免费久久久久久久清纯| 一级毛片我不卡| 桃色一区二区三区在线观看| 男人狂女人下面高潮的视频| 国产成人a区在线观看| 国产亚洲精品综合一区在线观看| 亚洲中文字幕日韩| 国产精品福利在线免费观看| 香蕉av资源在线| 中文字幕熟女人妻在线| 日本欧美国产在线视频| 人妻久久中文字幕网| 波多野结衣巨乳人妻| 黄色配什么色好看| 精品熟女少妇av免费看| 深夜a级毛片| 国产一区亚洲一区在线观看| 国产高清有码在线观看视频| 成人av在线播放网站| 日韩欧美国产在线观看| 国产伦一二天堂av在线观看| 久久久久久九九精品二区国产| 一进一出好大好爽视频| 观看美女的网站| 国产片特级美女逼逼视频| 国产精品久久久久久久久免| 日日摸夜夜添夜夜添小说| 人人妻人人澡人人爽人人夜夜 | 亚洲经典国产精华液单| 亚洲精品在线观看二区| 香蕉av资源在线| 成年女人永久免费观看视频| 91久久精品国产一区二区成人| 久久天躁狠狠躁夜夜2o2o| 一级黄片播放器| 人妻丰满熟妇av一区二区三区| 免费无遮挡裸体视频| 欧美+亚洲+日韩+国产| 丰满人妻一区二区三区视频av| 一夜夜www| 能在线免费观看的黄片| 级片在线观看| 白带黄色成豆腐渣| 自拍偷自拍亚洲精品老妇| 免费高清视频大片| 夜夜爽天天搞| 成年免费大片在线观看| 午夜视频国产福利| 插逼视频在线观看| www日本黄色视频网| 久久久久久大精品| 免费在线观看成人毛片| 99久久精品一区二区三区| 日本黄大片高清| 婷婷六月久久综合丁香| 69人妻影院| 在线看三级毛片| 国产伦精品一区二区三区四那| 亚洲国产高清在线一区二区三| 18禁黄网站禁片免费观看直播| 综合色av麻豆| 亚洲欧美中文字幕日韩二区| 丝袜喷水一区| 三级男女做爰猛烈吃奶摸视频| www日本黄色视频网| 激情 狠狠 欧美| 人妻少妇偷人精品九色| 免费黄网站久久成人精品| 国产伦一二天堂av在线观看| 国产高清视频在线播放一区| 天堂动漫精品| 最后的刺客免费高清国语| 女的被弄到高潮叫床怎么办| 91在线观看av| 天天躁夜夜躁狠狠久久av| 亚洲av免费在线观看| 性插视频无遮挡在线免费观看| 老司机福利观看| 欧美+日韩+精品| 国产av在哪里看| 哪里可以看免费的av片| 亚洲av免费高清在线观看| 尤物成人国产欧美一区二区三区| 看免费成人av毛片| 此物有八面人人有两片| 99热只有精品国产| 欧美成人a在线观看| 精品99又大又爽又粗少妇毛片| 日韩精品中文字幕看吧| 国产一区二区在线观看日韩| av在线播放精品| 少妇人妻精品综合一区二区 | 日本一二三区视频观看| 国产91av在线免费观看| 国产成人一区二区在线| 亚洲激情五月婷婷啪啪| 黑人高潮一二区| 久久久久国内视频| 99热这里只有精品一区| 大香蕉久久网| av国产免费在线观看| a级毛片a级免费在线| 亚洲五月天丁香| 特大巨黑吊av在线直播| 中国美白少妇内射xxxbb| 成人永久免费在线观看视频| 亚洲激情五月婷婷啪啪| 国产成人福利小说| 国产麻豆成人av免费视频| 日日摸夜夜添夜夜添小说| 99久久久亚洲精品蜜臀av| 国产一区二区在线av高清观看| 亚洲欧美清纯卡通| 欧洲精品卡2卡3卡4卡5卡区| 老熟妇乱子伦视频在线观看| 午夜a级毛片| 色播亚洲综合网| 欧美成人精品欧美一级黄| 亚洲国产精品合色在线| av国产免费在线观看| 狠狠狠狠99中文字幕| 一卡2卡三卡四卡精品乱码亚洲| 成人高潮视频无遮挡免费网站| eeuss影院久久| 毛片一级片免费看久久久久| 欧美高清性xxxxhd video| 免费无遮挡裸体视频| 草草在线视频免费看| 在线观看午夜福利视频| 老师上课跳d突然被开到最大视频| 伦理电影大哥的女人| 少妇被粗大猛烈的视频| 搡老妇女老女人老熟妇| 一级黄片播放器| 欧美高清成人免费视频www| 神马国产精品三级电影在线观看| 久久精品综合一区二区三区| 人妻少妇偷人精品九色| 99久久精品一区二区三区| 久久热精品热| 国产国拍精品亚洲av在线观看| 男插女下体视频免费在线播放| 久久精品国产鲁丝片午夜精品| 国产精品1区2区在线观看.| 欧美bdsm另类| 日本a在线网址| 精品久久久久久久久久久久久| 久久精品91蜜桃| 国产国拍精品亚洲av在线观看| 99热全是精品| 亚洲最大成人中文| 一进一出好大好爽视频| 高清毛片免费观看视频网站| 亚洲在线自拍视频| 午夜福利在线观看免费完整高清在 | 天堂av国产一区二区熟女人妻| 欧美人与善性xxx| 精品久久久久久久久亚洲| 69人妻影院| 国产精品伦人一区二区| 国产精品乱码一区二三区的特点| 国内揄拍国产精品人妻在线| 亚洲精品国产av成人精品 | 国产精品久久久久久精品电影| 亚洲国产精品成人综合色| 老女人水多毛片| 在线观看av片永久免费下载| 日韩欧美免费精品| 亚洲av电影不卡..在线观看| 欧美区成人在线视频| 日韩三级伦理在线观看| 简卡轻食公司| 在线播放无遮挡| 久久精品国产自在天天线| 午夜福利在线观看吧| 亚洲人成网站高清观看| 亚洲av中文av极速乱| 一级毛片久久久久久久久女| 色av中文字幕| 国产高清不卡午夜福利| 长腿黑丝高跟| 国产日本99.免费观看| 国产精品乱码一区二三区的特点| 欧美bdsm另类| 高清毛片免费观看视频网站| 日日摸夜夜添夜夜添小说| 欧美日本亚洲视频在线播放| 最近中文字幕高清免费大全6| 色尼玛亚洲综合影院| 两个人的视频大全免费| 丰满乱子伦码专区| 老司机午夜福利在线观看视频| av天堂中文字幕网| 天堂√8在线中文| 中出人妻视频一区二区| 亚洲专区国产一区二区| 国产精品一区二区三区四区久久| 99热这里只有是精品在线观看| 少妇高潮的动态图| 九九久久精品国产亚洲av麻豆| 一级毛片我不卡| 美女被艹到高潮喷水动态| 婷婷精品国产亚洲av在线| 91久久精品国产一区二区成人| 国产 一区 欧美 日韩| 亚洲精品亚洲一区二区| 亚洲av熟女| 成人综合一区亚洲| 国产三级在线视频| 国产成人a∨麻豆精品| 国产精品一区二区免费欧美| 99国产极品粉嫩在线观看| 久久韩国三级中文字幕| 中文字幕人妻熟人妻熟丝袜美| 蜜臀久久99精品久久宅男| 99国产精品一区二区蜜桃av| 精品久久久久久久久久久久久| 日本免费一区二区三区高清不卡| 天美传媒精品一区二区| 十八禁国产超污无遮挡网站| 男人舔奶头视频| 97热精品久久久久久| 91在线观看av| 一a级毛片在线观看| 国产人妻一区二区三区在| 亚洲国产精品合色在线| АⅤ资源中文在线天堂| 极品教师在线视频| 日本黄色视频三级网站网址| 欧美性猛交╳xxx乱大交人| 亚洲第一电影网av| 亚洲精品亚洲一区二区| 亚洲国产精品成人综合色| 成人综合一区亚洲| 91av网一区二区| 国产淫片久久久久久久久| 久久久久国产网址| 天堂动漫精品| 国产v大片淫在线免费观看| 国产女主播在线喷水免费视频网站 | 国产女主播在线喷水免费视频网站 | 色5月婷婷丁香| 国产精品久久久久久精品电影| 日韩 亚洲 欧美在线| 欧美最新免费一区二区三区| 亚洲最大成人av| 国产在线精品亚洲第一网站| 亚洲欧美日韩无卡精品| 男人的好看免费观看在线视频| 亚洲欧美成人精品一区二区| 高清日韩中文字幕在线| 日本黄色片子视频| 中文亚洲av片在线观看爽| 日韩欧美免费精品| 久久久久久久亚洲中文字幕| 亚洲自拍偷在线| 国内精品宾馆在线| 日韩,欧美,国产一区二区三区 | 久久久久久久亚洲中文字幕| 五月伊人婷婷丁香| 亚洲欧美中文字幕日韩二区| 久久鲁丝午夜福利片| 久久精品综合一区二区三区| 久久久久久久久久久丰满| 亚洲真实伦在线观看| 午夜免费激情av| 伊人久久精品亚洲午夜| 长腿黑丝高跟| 国产久久久一区二区三区| av天堂在线播放| av女优亚洲男人天堂| h日本视频在线播放| 激情 狠狠 欧美| 国产成人a区在线观看| 久久人人精品亚洲av| videossex国产| 99精品在免费线老司机午夜| 一进一出好大好爽视频| 国产69精品久久久久777片| 国产熟女欧美一区二区| 天美传媒精品一区二区| 亚洲三级黄色毛片| 热99在线观看视频| 成人av在线播放网站| 99视频精品全部免费 在线| 乱人视频在线观看| 午夜精品一区二区三区免费看| 成人特级黄色片久久久久久久| 免费在线观看影片大全网站| av卡一久久| 国产亚洲精品av在线| 黑人高潮一二区| 亚洲国产精品久久男人天堂| 午夜影院日韩av| 亚洲国产精品成人久久小说 | 九九在线视频观看精品| 一个人免费在线观看电影| 亚洲在线自拍视频| 婷婷六月久久综合丁香| 国产单亲对白刺激| 一进一出抽搐gif免费好疼| 一区二区三区免费毛片| 日本熟妇午夜| 激情 狠狠 欧美| 国产色爽女视频免费观看| 精品一区二区免费观看| 国产精品一二三区在线看| av在线蜜桃| 97热精品久久久久久| 国产精品乱码一区二三区的特点| 91久久精品国产一区二区三区| 国产大屁股一区二区在线视频| 国产 一区精品| 51国产日韩欧美| 人妻夜夜爽99麻豆av| 日本一二三区视频观看| 色综合色国产| 亚洲七黄色美女视频| 久久久久久久久中文| 黑人高潮一二区| 直男gayav资源| 久久人人精品亚洲av| 最近的中文字幕免费完整| 精品一区二区三区视频在线观看免费| 欧美zozozo另类| 赤兔流量卡办理| 少妇被粗大猛烈的视频| 六月丁香七月| 国产又黄又爽又无遮挡在线| av天堂中文字幕网| 国产一区二区在线观看日韩| 久久人人爽人人爽人人片va| 欧美性猛交╳xxx乱大交人| 能在线免费观看的黄片| 夜夜看夜夜爽夜夜摸| 一区二区三区免费毛片| 伊人久久精品亚洲午夜| 最近视频中文字幕2019在线8| 午夜久久久久精精品| 18禁黄网站禁片免费观看直播| 成年免费大片在线观看| 此物有八面人人有两片| 99热6这里只有精品| 亚洲激情五月婷婷啪啪| 国产成人aa在线观看| 亚洲av中文av极速乱| 国产人妻一区二区三区在| 国产女主播在线喷水免费视频网站 | 久久人人爽人人爽人人片va| 又黄又爽又刺激的免费视频.| 色噜噜av男人的天堂激情| 如何舔出高潮| 搡老岳熟女国产| 国产亚洲欧美98| 村上凉子中文字幕在线| 美女 人体艺术 gogo| 12—13女人毛片做爰片一| 18禁裸乳无遮挡免费网站照片| av天堂中文字幕网| 毛片一级片免费看久久久久| 特级一级黄色大片| 国产黄色小视频在线观看| 欧美人与善性xxx| 黄色视频,在线免费观看| 在线国产一区二区在线| 亚洲精品国产成人久久av| 亚洲丝袜综合中文字幕| 97超视频在线观看视频| 国产成人a区在线观看| www.色视频.com| 高清毛片免费看| 日本 av在线| 欧美高清性xxxxhd video| 最新中文字幕久久久久| 亚洲av不卡在线观看| 美女大奶头视频| 极品教师在线视频| 精华霜和精华液先用哪个| 老司机午夜福利在线观看视频| 国产 一区精品| 亚洲av不卡在线观看| eeuss影院久久| 蜜臀久久99精品久久宅男| 少妇熟女欧美另类| 嫩草影视91久久| 欧美色欧美亚洲另类二区| 蜜臀久久99精品久久宅男| 国产日本99.免费观看| 日本五十路高清| 久久99热6这里只有精品| 国产精品福利在线免费观看| 成年免费大片在线观看| 高清午夜精品一区二区三区 | 毛片女人毛片| 久久精品91蜜桃| 在线免费十八禁| 日韩av不卡免费在线播放| 在线免费观看不下载黄p国产| 男人舔奶头视频| 成人鲁丝片一二三区免费| 97人妻精品一区二区三区麻豆| 欧美极品一区二区三区四区| av国产免费在线观看| 成年av动漫网址| 日本 av在线| 天堂影院成人在线观看| 免费不卡的大黄色大毛片视频在线观看 | 日本黄色片子视频| 美女大奶头视频| 校园春色视频在线观看| 男女下面进入的视频免费午夜| 亚洲图色成人| 91麻豆精品激情在线观看国产| 国产又黄又爽又无遮挡在线| 国产精品一区二区三区四区久久| 欧美一级a爱片免费观看看| 国产亚洲91精品色在线| 波多野结衣巨乳人妻| 欧美一区二区亚洲| 伊人久久精品亚洲午夜| 卡戴珊不雅视频在线播放| 在线a可以看的网站| 老司机影院成人| 久久精品人妻少妇| 精品无人区乱码1区二区| 亚洲美女黄片视频| 一卡2卡三卡四卡精品乱码亚洲| 草草在线视频免费看| 亚洲av成人精品一区久久| 国产一区二区亚洲精品在线观看| 亚洲精华国产精华液的使用体验 | 观看美女的网站| 在线观看66精品国产| 久久午夜亚洲精品久久| 亚洲av中文av极速乱| 亚洲欧美精品综合久久99| 色综合亚洲欧美另类图片| 久久久精品大字幕| 国产成人福利小说| 亚洲人成网站高清观看| 菩萨蛮人人尽说江南好唐韦庄 | 在线观看av片永久免费下载| 99久久精品国产国产毛片| 成人亚洲欧美一区二区av| 麻豆国产av国片精品| 身体一侧抽搐| 毛片女人毛片| 在线观看免费视频日本深夜| 午夜精品在线福利| 日产精品乱码卡一卡2卡三| 国产精品久久电影中文字幕| 五月伊人婷婷丁香| 啦啦啦啦在线视频资源|