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

    定常黏性外流氣動(dòng)力分量的同源性及流場(chǎng)斷層掃描診斷原理

    2019-12-30 05:25:04鄒舒帆吳介之
    關(guān)鍵詞:尾流氣動(dòng)力升力

    鄒舒帆, 吳介之

    (北京大學(xué) 工學(xué)院 湍流與復(fù)雜系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室, 北京 100871)

    0 引 言

    減阻是外流空氣動(dòng)力學(xué)的基本任務(wù)之一,大飛機(jī)減阻是當(dāng)代航空技術(shù)的一個(gè)迫切需求。減阻的前提是對(duì)阻力主要組分做出正確的物理解釋和準(zhǔn)確的理論預(yù)測(cè)。在傳統(tǒng)流場(chǎng)診斷和減阻方案中采用的原則是線性分解、各個(gè)擊破和線性疊加,并據(jù)此選定合適的阻力分類方式,追溯阻力各組分的物理來源,尋找對(duì)應(yīng)的減阻措施[1]。例如,采用改變壁面條件來控制摩擦阻力[2],采用吹吸氣抑制流場(chǎng)分離來控制壓差阻力[3],采用翼梢小翼來控制誘導(dǎo)阻力,等等。在實(shí)際工程應(yīng)用中則把對(duì)各個(gè)阻力的相應(yīng)減阻方案結(jié)合起來以達(dá)到總體減阻的目標(biāo)。

    但是,如我們最近的研究[4]評(píng)述的,當(dāng)前國內(nèi)外的減阻思路存在若干根本性弊病,普遍的問題包括:阻力分類不唯一,阻力成分占比預(yù)測(cè)不精確,總體減阻結(jié)果也并非各減阻措施的簡(jiǎn)單線性疊加。

    其實(shí),對(duì)于任意構(gòu)型的定常繞流(包括雷諾時(shí)均定常湍流),在亞聲速情況下只有兩種阻力分解方式具有嚴(yán)格的理論基礎(chǔ)[5]:第一種把物面應(yīng)力積分得到的總阻力D分解成壓差阻力Dp和摩擦阻力Df,其物理意義是清晰而無歧義的;第二種則用流場(chǎng)動(dòng)量平衡的體積分得到誘導(dǎo)阻力Di和型阻DP。這里的流場(chǎng)體積分可以轉(zhuǎn)換為外邊界的面積分,再通過選取合適的控制面簡(jiǎn)化為適當(dāng)尾流截面的積分。第一種分解方式是普適的,對(duì)任何物體在任何流動(dòng)條件下都成立,但無法將阻力與實(shí)際流場(chǎng)中的結(jié)構(gòu)聯(lián)系起來,故需轉(zhuǎn)向基于動(dòng)量平衡積分的第二種分解方式??墒牵捎贑FD和EFD都無法分別單獨(dú)測(cè)量Di和DP,它們中至少一個(gè)必須先給出精確的理論定義 (另一個(gè)可由總阻力推出)才知道要測(cè)量或計(jì)算什么,其可靠性完全取決于理論的可靠性。理論好,則定義合理可信;若理論過于簡(jiǎn)化,則定義不準(zhǔn),減阻措施打不中要害,而毛病何在卻是根據(jù)定義算不出測(cè)不到的?,F(xiàn)在工程中常用的誘導(dǎo)阻力公式均基于線化尾流模型,引入Trefftz平面假設(shè),通過展向升力[6]或者尾流面上流向渦量分布[7]來確定誘導(dǎo)阻力。這些公式對(duì)小攻角大展弦比薄翼的附著流基本適用,但對(duì)真實(shí)的復(fù)雜流動(dòng)卻偏差很大。為此,Zou等[4]對(duì)該問題進(jìn)行了深入的研究,提出了三維不可壓縮定常黏性流動(dòng)中基于尾流積分的氣動(dòng)力分量精確理論。

    這個(gè)理論仍基于傳統(tǒng)的線性思維,卻揭示了兩個(gè)無法回避的問題。一是在給定構(gòu)型和流動(dòng)條件下,型阻和誘導(dǎo)阻力不再是純數(shù)而是尾流截面位置的函數(shù)。Zou等[4]從理論和數(shù)值上確認(rèn)了該變化的趨勢(shì),并查明了其隨尾流截面位置變化的物理根源。型阻和誘導(dǎo)阻力客觀性的喪失,使得尋找合適阻力分量測(cè)量位置、確認(rèn)減阻措施的有效性成了無法解決的難題。二是升力、型阻和誘導(dǎo)阻力是物理上同源的,可分別采用Lamb矢量的體積分或其矩的面積分來定義。Lamb矢量是所有氣動(dòng)力的核心。這個(gè)觀察在Wu等[8]的書中已初現(xiàn)端倪,但未系統(tǒng)追究。它的結(jié)論是:若要通過線性思維減少某阻力分量,必會(huì)對(duì)其余氣動(dòng)力分量產(chǎn)生影響。氣動(dòng)力分量的非客觀性和同源性使得原有的基于線性分解、各個(gè)擊破、線性疊加的傳統(tǒng)流動(dòng)診斷和減阻方案無法應(yīng)對(duì)真實(shí)非線性流場(chǎng)。

    針對(duì)上述問題,我們提出了兩種不同的解決方案。方案一保留傳統(tǒng)氣動(dòng)力分解框架,給出附著流的誘導(dǎo)阻力中場(chǎng)近似。該方案兼顧了阻力分解的精確性和公式的易用性,適用于常規(guī)的大飛機(jī)氣動(dòng)設(shè)計(jì)。方案二突破傳統(tǒng)的阻力分解框架,構(gòu)造面向任意復(fù)雜構(gòu)型氣動(dòng)力診斷的新型斷層掃描技術(shù),可稱為空氣動(dòng)力學(xué)CT,其理念來源于現(xiàn)代空氣動(dòng)力學(xué)理論。Wu等[9]述評(píng)了關(guān)于黏性可壓縮復(fù)雜流動(dòng)的現(xiàn)代空氣動(dòng)力學(xué)基礎(chǔ)理論及其CFD/EFD檢驗(yàn),其中提出了現(xiàn)代空氣動(dòng)力學(xué)理論在近場(chǎng)可以分為三個(gè)不同的層次:速度壓力層次、結(jié)構(gòu)層次以及因果層次。所有這些層次的理論相互豐富,而非相互排斥。對(duì)于復(fù)雜流動(dòng),不同層次的理論應(yīng)該聯(lián)合使用,并與 CFD/EFD 緊密結(jié)合。結(jié)合CFD/EFD的現(xiàn)實(shí)情況可知,非線性近場(chǎng)得到的數(shù)據(jù)是最為準(zhǔn)確的,也是利用有限的CFD/EFD數(shù)據(jù)進(jìn)行局部動(dòng)力學(xué)診斷和優(yōu)化的最佳場(chǎng)所。若需要,甚至可以將截面切入物體內(nèi)部.這種觀念上的解放,意味著不再糾結(jié)于尋找“客觀”的阻力分類方式,而是轉(zhuǎn)向研究流場(chǎng)結(jié)構(gòu)與氣動(dòng)力的直接聯(lián)系并追溯其在壁面的產(chǎn)生根源。

    本文介紹筆者的課題組基于上述認(rèn)識(shí)的研究成果。為敘述方便完整起見,本文第 1節(jié)回顧Zou等[4]的研究成果,第2節(jié)系統(tǒng)闡述氣動(dòng)力分量的同源性及其本征機(jī)理,第3節(jié)提出空氣動(dòng)力學(xué)CT技術(shù),展示新的進(jìn)展。本文僅限于不可壓流。

    1 定常黏性不可壓外流的氣動(dòng)力分解理論及緊致中場(chǎng)近似

    考慮一個(gè)繞任意飛行器的三維定常黏性流動(dòng),均勻來流速度為U=Uex,飛行器B的外邊界為?B。圖 1定義了固定于飛行器的坐標(biāo)系(x,y,z),流向坐標(biāo)x的原點(diǎn)位于機(jī)翼前緣,并列出了下面將用到的標(biāo)記:控制面∑包圍了流體區(qū)域Vf以及飛行器B,因此V=Vf∪B。流體域外邊界?Vf的單位法向量指向流體外部。在大雷諾數(shù)下,∑上的黏性力可以忽略。控制面∑僅在尾流處切割渦量場(chǎng),其前面和側(cè)面上的流動(dòng)是無黏無旋的。尾流中的渦量ω≠0僅僅在尾流面上很小的一個(gè)區(qū)域Wv?W內(nèi)存在。本文所有物理量均采用來流速度U、機(jī)翼弦長(zhǎng)c以及機(jī)翼面積S無量綱化。

    圖1 流體分析區(qū)域以及標(biāo)記

    首先回顧當(dāng)前常用的定常不可壓層流誘導(dǎo)阻力和型阻公式(時(shí)均湍流中的公式只需略加修改),再簡(jiǎn)述精確的一般理論,探討氣動(dòng)力客觀性喪失的物理根源,最后提出保留傳統(tǒng)氣動(dòng)力分解框架的誘導(dǎo)阻力中場(chǎng)近似。

    1.1 常用型阻和誘導(dǎo)阻力公式

    采用尾流積分(W-積分)定義的型阻是清晰且明確的[10-11]:

    (1)

    其中P=p+ρ|u|2/2是總壓,p、ρ和u分別是壓強(qiáng)、密度和速度。式(1)適用于任何飛行器的定常繞流。與此相反,誘導(dǎo)阻力Di的尾流積分表達(dá)式卻并不唯一。在當(dāng)前的工程應(yīng)用中,有兩個(gè)簡(jiǎn)單的Di定義公式。一個(gè)是經(jīng)典的Prandtl[6,12]公式,無黏且DP=0,

    (2)

    其中Λ和An分別是展弦比和升力展向分布的Fourier展開系數(shù)。對(duì)于橢圓升力分布,δ=0,從而得到單個(gè)機(jī)翼的最小誘導(dǎo)阻力。另一個(gè)是Maskell公式[7],可用于黏流:

    (3)

    其中ωx是流向渦量,2ψx=-ωx,r=yey+zez是遠(yuǎn)場(chǎng)Trefftz平面T上的位置矢量。公式(2)和公式(3)均基于線化尾流模型,在尾流截面上只保留流向渦量ω=(ωx,0,0)。該假設(shè)使得公式(2)和公式(3)在小攻角大展弦比薄翼的附著流中基本適用[13],但Zou等[4]指出它們?cè)谡鎸?shí)復(fù)雜流動(dòng)中有較大偏差。

    過去30年中,如何定義一般流動(dòng)情況下的誘導(dǎo)阻力一直是個(gè)非常活躍的研究領(lǐng)域[5,14-22]。在這些討論中,Di的另一個(gè)通用公式常作為研究的起點(diǎn)。令u=U+v=(U+u′,v,w),有[5,16-22]

    (4)

    其中k=|v|2/2是擾動(dòng)動(dòng)能。然而,公式(4)僅僅意味著誘導(dǎo)阻力可以表示為被尾流對(duì)流帶走的擾動(dòng)動(dòng)能減去流向擾動(dòng)動(dòng)量,物理含義含糊,而且被積函數(shù)在尾流中彌散。這個(gè)公式不適用于實(shí)際的尾流測(cè)量?,F(xiàn)有的理論推導(dǎo)(可壓縮流動(dòng))主要采用攝動(dòng)方法,Maskell公式(3)則被證明是個(gè)一階近似[16]??偟膩砜?,在非線性黏性流動(dòng)中如何精確定義誘導(dǎo)阻力并解釋其物理含義,一直沒有顯著的進(jìn)展,這也導(dǎo)致了對(duì)整個(gè)誘導(dǎo)阻力預(yù)測(cè)理論“長(zhǎng)期存在的不滿”(Spalart[22])。

    1.2 型阻與誘導(dǎo)阻力的精確定義

    為精確定義基于尾流積分的阻力分量,通常從動(dòng)量方程出發(fā),方程(4)就是其中的一個(gè)結(jié)果。其實(shí),誘導(dǎo)阻力的概念和理論源于Prandtl[6],但很少人注意到在同一篇論文中,Prandtl 已經(jīng)給出了定常無黏流升阻力的一般渦力理論,其誘導(dǎo)阻力公式僅是這個(gè)理論在升力線近似下簡(jiǎn)化到最低階的解析表示。Wu等[8]在Prandtl渦力理論的基礎(chǔ)上,將其進(jìn)一步拓展到了黏性流動(dòng)中,并采用體積分以及尾流截面積分,分別給出了誘導(dǎo)阻力和型阻的精確的一般定義式:

    (5)

    ≡DP0+DP1

    (6)

    公式(5)是一個(gè)具有明確物理意義的緊致積分,Marongiu等[13]已用數(shù)值計(jì)算驗(yàn)證了它的正確性。它可以嚴(yán)格地轉(zhuǎn)換為公式(4),說明誘導(dǎo)阻力的體積分定義和W上的非緊致積分等價(jià)。公式(6)則是從公式(1)直接通過導(dǎo)數(shù)矩變換得來,兩者也等價(jià)。此式表明,型阻的物理載體是物體定常運(yùn)動(dòng)產(chǎn)生的Lamb矢量場(chǎng)的矩,積分區(qū)域是緊致的。這對(duì)于物理理解、精確計(jì)算測(cè)量和流動(dòng)診斷具有極其重要的作用。

    為了實(shí)際計(jì)算和測(cè)量的方便,積分區(qū)域應(yīng)越小越好。型阻的一般定義式(6)滿足該需求,而誘導(dǎo)阻力的一般定義式(5)可進(jìn)一步轉(zhuǎn)換為尾流截面積分,但留下明顯依賴于x的非緊致項(xiàng):

    (7)

    對(duì)于上述嚴(yán)格的三維黏性尾流截面積分理論,Zou等[4]的論文中采用兩個(gè)典型的定常流計(jì)算結(jié)果——雷諾Re=1×105、攻角α=4°的繞Xt=1.0橢圓機(jī)翼的附著流以及Re=5×105、α=20°的繞后掠角χ=76°細(xì)長(zhǎng)三角翼的分離流——分別驗(yàn)證了其在附著流和分離流中的準(zhǔn)確性與有效性。具體細(xì)節(jié)可以參閱該論文,在此我們只關(guān)注型阻DP和誘導(dǎo)阻力Di與總阻力D隨x的演化關(guān)系,如圖2所示。由公式(6)+公式(7)計(jì)算出的總阻力與壁面應(yīng)力積分計(jì)算出的結(jié)果符合得很好,兩者最大的相對(duì)誤差小于2%。我們發(fā)現(xiàn)在近尾流區(qū)域,Di和DP都是x相關(guān)的。橢圓機(jī)翼的Di和DP在機(jī)翼尾緣隨x變化較大,但隨后便基本保持不變。但是,三角翼的DP一直單調(diào)上升,而Di最初占據(jù)了D的絕大部分并在后緣附近達(dá)最大值,然后單調(diào)減少,這導(dǎo)致在中場(chǎng)以后DP占了D的絕大部分。已有理論證明[14,23-25],隨x→∞,有Di→0,DP→D。

    (a)橢圓機(jī)翼,Λ=7,Re=1×105, α=4°

    (b)三角翼,χ=76°,Re=5×105, α=20°

    1.3 氣動(dòng)力分量x依賴性的物理根源

    1.2節(jié)從理論和數(shù)值兩個(gè)方面確認(rèn)了:氣動(dòng)力分量總是具有x依賴性;這種依賴性在附著流中比較弱,即在近尾流比較明顯,而在遠(yuǎn)尾流基本可以認(rèn)為不存在;但在定常分離流中卻很強(qiáng):在近尾流和遠(yuǎn)尾流都有著明顯的體現(xiàn)。

    進(jìn)而,Zou等[4]找到Di和DP的x依賴性根源。在運(yùn)動(dòng)學(xué)上,對(duì)公式(5)沿x方向求導(dǎo),可得:

    (8)

    這表明誘導(dǎo)阻力隨x的變化取決于尾流截面上的擾動(dòng)Lamb矢量的x分量積分。若將擾動(dòng)速度和渦量分解為v=vxex+vπ,ω=ωxex+ωπ,下標(biāo)π代表著尾流W截面上的切向向量,有

    ex·(v×ω)=ex·(vπ×ωπ)

    (9)

    因此,只要在W(x)上vπ×ωπ≠0,從尾流平面上算得的誘導(dǎo)阻力肯定是x依賴的。在動(dòng)力學(xué)上,x依賴性可從物體做功和流體區(qū)域Vf的動(dòng)能變化求得:

    (10)

    現(xiàn)在的問題是:在可接受的誤差范圍內(nèi)是否存在一個(gè)尾流區(qū)間,使dDi/dx=-dDp/dx?0?這個(gè)關(guān)鍵判據(jù)可以用來檢驗(yàn)誘導(dǎo)阻力和型阻是否可以被視為與尾流截面位置無關(guān)的常數(shù)。1.2節(jié)的計(jì)算結(jié)果表明,這個(gè)問題的答案取決于具體物體構(gòu)型和流動(dòng)條件。粗略地說,完全附著流中的阻力分量可以視為常數(shù),而分離流中的阻力分量必然是尾流截面位置的函數(shù)。后者涵蓋了比前者更廣泛的物體構(gòu)型和流動(dòng)條件,包括繞機(jī)翼機(jī)身組合體到整個(gè)飛機(jī)的定常流動(dòng)。

    1.4 誘導(dǎo)阻力的緊致中場(chǎng)近似

    前已說明,x依賴性的存在使得選取合適的位置測(cè)量誘導(dǎo)阻力和型阻成了一個(gè)關(guān)鍵問題。由式(10)可知,在此尾流位置上必須忽略耗散。根據(jù)上述分析,這在近尾流和遠(yuǎn)尾流都是不可能的,唯一的選擇是兩者之間的弱非線性區(qū)域,稱之為中場(chǎng)尾流區(qū)。在其中的尾流平面W上,誘導(dǎo)阻力所固有的核心非線性效應(yīng)即擾動(dòng)Lamb矢量v×ω應(yīng)該保留下來。最終可把公式(7)簡(jiǎn)化為比公式(2)和公式(3)準(zhǔn)確的緊致中場(chǎng)近似:

    (11)

    (a)橢圓機(jī)翼

    (b)三角翼

    2 合力組分的物理同源性與縱橫分解

    回顧誘導(dǎo)阻力公式(5)和型阻公式(6),并考慮渦力理論中升力的公式

    L=L0+L1

    (12)

    可以看出,型阻的線性項(xiàng)DP0和線性升力L0一樣,來自線化Lamb矢量U×ω;型阻的非線性項(xiàng)DP1與非線性升力L1以及誘導(dǎo)阻力Di一樣,來自非線性Lamb矢量v×ω。它們之間的關(guān)系如圖4所示。

    圖4 合力各組分的物理同源性示意圖

    這個(gè)觀察表明,定常物體運(yùn)動(dòng)產(chǎn)生的Lamb矢量場(chǎng)u×ω及其演化是理解所有氣動(dòng)力分量的關(guān)鍵。誘導(dǎo)阻力與升力和型阻本征相關(guān),說明對(duì)它們必須做一體化考慮,很難分割開來實(shí)施優(yōu)化。以三角翼尾流為例,如圖5所示,在x=1.2截面,尾渦對(duì)不僅提供了升力,也同時(shí)提供了Di和Dp。

    為了深入理解升阻力各分量的同源性和本征關(guān)聯(lián),我們從 Lamb矢量的縱橫分解以及其發(fā)展演化入手。眾所周知,任意一個(gè)分片可微的向量場(chǎng)f可以分解為一個(gè)標(biāo)量場(chǎng)的梯度和一個(gè)向量場(chǎng)的旋度之和,即著名的Helmholtz分解定理:

    f=φ+×ψ,·ψ=0

    (13)

    φ和ψ稱為f的標(biāo)量勢(shì)和向量勢(shì)。它們代表的流動(dòng)過程分別稱為縱過程和橫過程,因此 Helmholtz 分解也叫縱橫分解。如果勢(shì)函數(shù)本身也是可測(cè)的物理量,則稱此縱橫分解為自然的,十分可貴。對(duì)于 Lamb矢量l=ω×u,Wu等[27]指出它有自然的 Helmholtz 分解

    (a)ωx

    (b)ρex·(v×ω)

    (c)ρez·(u×ω)

    (d)0.5ρr·(u×ω)

    l=l‖+l⊥=-P-×(νω)

    l‖=-P,l⊥=-×(νω)

    (14)

    可見,縱Lamb矢量l‖與總壓梯度平衡,總壓就是Lamb矢量的標(biāo)量勢(shì)。正因如此,可以將型阻的面積分轉(zhuǎn)換成Lamb的矩的積分.對(duì)上述方程兩邊在Vf中進(jìn)行積分,將梯度項(xiàng)變換成面積分,并分別在x方向投影得到:

    (15)

    (16)

    這里τ=μω×n是切應(yīng)力。因此,縱Lamb矢量l‖的x分量體積分代表著壁面壓差阻力與流場(chǎng)型阻的差。這個(gè)體積分也可以寫為-Di‖。橫 Lamb矢量l⊥的體積分完全提供了壁面摩擦力,它又是誘導(dǎo)阻力的橫分量。因此,Lamb矢量的縱橫部分將壁面的應(yīng)力與流域/尾流積分的定義簡(jiǎn)潔而明確的聯(lián)系起來,如圖6所示。

    圖6 基于Lamb矢量縱橫分解的總阻力及其分量示意圖

    Wu等[27]還觀察到,Lamb矢量完全由橫分量l⊥以及其全場(chǎng)效應(yīng)所驅(qū)動(dòng) (非定常時(shí)還包括歷史積累效應(yīng))。正是這種非線性的全場(chǎng)效應(yīng)導(dǎo)致了包括湍流在內(nèi)的各種旋渦流動(dòng)的極為復(fù)雜的演化。橫分量l⊥變化引起了縱分量l‖的被動(dòng)變化,它類似于壓力梯度的行為。然而,與壓力本身不同的是,由于l‖=-P中包含動(dòng)能,所以它受擴(kuò)散方程的控制。該現(xiàn)象與型阻和黏性擴(kuò)散之間的內(nèi)稟物理關(guān)聯(lián)一致。

    3 空氣動(dòng)力學(xué)斷層掃描技術(shù)

    第1節(jié)指出了傳統(tǒng)阻力分解的種種困難,說明只對(duì)附著流可以建立合理的中場(chǎng)近似。第2節(jié)發(fā)現(xiàn)在所有的氣動(dòng)力分量在物理上同源于 Lamb矢量的體積分或矩積分,升力以及阻力分量無法分割開來實(shí)施優(yōu)化。若要考慮更加復(fù)雜的構(gòu)型(如全機(jī)或嶄新氣動(dòng)概念飛行器)的氣動(dòng)設(shè)計(jì),就需要突破傳統(tǒng)的阻力分解框架,構(gòu)建新型的技術(shù)路線。本節(jié)旨在為這個(gè)根本性的創(chuàng)新提供一個(gè)堅(jiān)實(shí)的基礎(chǔ)。

    3.1 空氣動(dòng)力學(xué)斷層掃描診斷的理論基礎(chǔ)

    和對(duì)升力來源的認(rèn)識(shí)發(fā)展一樣,研究阻力也需要突破線性分解、各個(gè)擊破、線性疊加的傳統(tǒng)思維和用速度、動(dòng)能解釋誘導(dǎo)阻力的籠統(tǒng)認(rèn)識(shí),樹立整體的、相互作用的思維,從速度、動(dòng)能追問到具體的流動(dòng)結(jié)構(gòu)及其產(chǎn)生機(jī)制?;谶@一考量,結(jié)合CFD/EFD的現(xiàn)實(shí)情況可知,非線性近場(chǎng)得到的數(shù)據(jù)是最為準(zhǔn)確的,也是利用有限的CFD/EFD數(shù)據(jù)進(jìn)行局部動(dòng)力學(xué)診斷和優(yōu)化的最佳場(chǎng)所。因此,研究的對(duì)象應(yīng)由速度、動(dòng)能轉(zhuǎn)變?yōu)閷?dǎo)致升阻力的緊致 Lamb矢量場(chǎng)的縱橫部分在近場(chǎng)的非線性演化特性,并追溯其產(chǎn)生機(jī)制。如圖7所示,尾流截面甚至可以切入物體內(nèi)部,將Lamb矢量場(chǎng)及其相關(guān)邊界量隨流向x或展向y的截面積分變化(可拓展到任選方位的截面族)作為診斷和溯源的重要線索。

    圖7 近場(chǎng)空氣動(dòng)力學(xué)CT截面示意圖

    這種通過掃描、切片得到一系列流動(dòng)結(jié)構(gòu)的截面數(shù)據(jù)并進(jìn)行局部流場(chǎng)結(jié)構(gòu)的動(dòng)力學(xué)分析,并結(jié)合邊界渦量流追溯其產(chǎn)生機(jī)制的技術(shù),我們稱為空氣動(dòng)力學(xué)斷層掃描技術(shù) (亦即空氣動(dòng)力學(xué)CT 技術(shù))。

    為此,我們首先寫出一般不可壓定常黏性流動(dòng)的空氣動(dòng)力學(xué)CT渦力公式為:

    (17)

    其中:

    分別為外邊界項(xiàng)以及物面邊界項(xiàng),σ=μ?ω/?n,m=2或3,m為空間維數(shù)。與通常的不可壓縮渦力理論[28]相比,出現(xiàn)了由于切割物體產(chǎn)生的線積分項(xiàng)FSB。

    若把基于尾流截面積分的氣動(dòng)力公式推廣至空氣動(dòng)力學(xué)CT中,也會(huì)出現(xiàn)相應(yīng)的邊界積分,在此不再贅述。

    上述兩種理論屬于結(jié)構(gòu)層次[9],它只關(guān)注流動(dòng)結(jié)構(gòu)如何決定氣動(dòng)力,并不關(guān)注流動(dòng)結(jié)構(gòu)如何形成。后一個(gè)問題涉及復(fù)雜的因果關(guān)系鏈條,其最主要和最基本的機(jī)理是渦量如何在壁面產(chǎn)生[9],因此還需要考慮基于邊界渦量流的合力積分公式,它在三維流中是:

    (18)

    其中:

    σp=n×p,σvis=μ(n×)×ω

    邊界渦量流σ=σp+σvis(簡(jiǎn)寫為 BVF)衡量了單位面積的壁面在單位時(shí)間內(nèi)由于無滑移條件產(chǎn)生并擴(kuò)散進(jìn)入流體內(nèi)部的渦量大小。它代表著由于黏附性條件而通過切向壓力梯度產(chǎn)生渦量的機(jī)制,因此也展現(xiàn)了兩個(gè)基本過程——縱橫過程——在邊界上的黏性與線性耦合。壁面是產(chǎn)生渦量的地方,因此也是最初形成 Lamb矢量的地方。雖然由于黏附條件而在靜止邊界上有u=0從而l=0,但其縱場(chǎng)和橫場(chǎng)部分并不分別為0,而是在壁面有:

    ρl⊥+μ×ω=0,ρl‖+p=0

    (19)

    Wu等[27]已經(jīng)推導(dǎo)出不可壓縮流動(dòng)的Lamb矢量輸運(yùn)方程。該方程有一個(gè)擴(kuò)散項(xiàng)μ2l的體積分,可以轉(zhuǎn)換成邊界積分μn·l,該形式與BVF相似且密切相關(guān)??梢缘玫絃amb矢量的縱橫分量在壁面上與BVF 之間的關(guān)系:

    ρn×l⊥=-ρn×l‖=-μ(n×)×ω=n×p=σ

    (20)

    ρn·l⊥=-ρn·l‖=-μ(n×)·ω=n·p

    (21)

    可見,橫向Lamb矢量在壁面的切向分量即為BVF,法向分量為法向壓力梯度。在壁面產(chǎn)生的橫向Lamb矢量首先通過擴(kuò)散進(jìn)入流體內(nèi)部靠近邊界的黏性子層中,因?yàn)樵搮^(qū)域有著很強(qiáng)的黏性作用。

    具體應(yīng)用空氣動(dòng)力學(xué)CT時(shí),可遵循Wu等[9]的建議:面對(duì)特定的流動(dòng)問題,首先考察全局流場(chǎng),然后更多地關(guān)注局部結(jié)構(gòu)和過程,最后查明其物理原因。為了獲得流場(chǎng)多層次的視野,需要把所有可能層次的理論有機(jī)地結(jié)合起來使用,且與CFD/EFD得到的精確數(shù)據(jù)結(jié)合。為了展示空氣動(dòng)力學(xué)CT的具體應(yīng)用,下一節(jié)將以第2節(jié)中的細(xì)長(zhǎng)三角翼分離流為例,在相同的計(jì)算參數(shù)下,對(duì)其升阻力進(jìn)行流場(chǎng)診斷,并據(jù)此提出新的概念性的增升減阻方案。

    3.2 基于空氣動(dòng)力學(xué)CT的流場(chǎng)診斷

    先從全局流場(chǎng)的角度來檢查三角翼在不同攻角情況下升阻力的變化曲線,如圖8所示。圖8(a)表明,升力基本是線性增長(zhǎng),而阻力按二次曲線增長(zhǎng)。因此升阻比L/D一開始急劇增長(zhǎng),在α=4°時(shí)達(dá)到了最大值。此后L/D一直下降并逐漸平緩。

    為了定量的探討氣動(dòng)力隨攻角的變化關(guān)系,取α=4°以及α=20°兩個(gè)代表性的算例進(jìn)行分析。在結(jié)構(gòu)層次的理論中,氣動(dòng)力的變化與渦息息相關(guān)。流動(dòng)結(jié)構(gòu)在空間中并不是均勻分布的,參照Yang等[29]的方法,將流場(chǎng)分為三個(gè)區(qū):區(qū)域I(Region I)是下表面以下區(qū)域,區(qū)域III(Region III)是脫體渦區(qū),區(qū)域II(Region II)是兩者中間的上表面近壁區(qū)。

    (a)升阻力系數(shù) (b)升阻力

    圖9 流場(chǎng)分區(qū)示意圖(背景為α=20°時(shí)細(xì)長(zhǎng)三角翼x/c=0.8截面ωx云圖)

    當(dāng)然更好的方式是將脫體渦單獨(dú)隔離出來計(jì)算。但數(shù)值結(jié)果表明,只要保證上表面近壁區(qū)處于區(qū)域II 中,其值就很穩(wěn)定。連接集中渦與近壁層的剪切層對(duì)升阻力貢獻(xiàn)都較小。值得注意的是,這種分區(qū)并不意味著各區(qū)域的效果是獨(dú)立的:近壁區(qū)產(chǎn)生的渦量從翼尖分離進(jìn)入流場(chǎng)形成自由剪切層,它進(jìn)而卷繞成集中的脫體渦;而集中渦又會(huì)反過來影響上翼面近壁區(qū)的流場(chǎng),誘導(dǎo)出二次渦。這種渦結(jié)構(gòu)的產(chǎn)生和誘導(dǎo)都是相互的,因此空間上的分解是對(duì)結(jié)果的描述,而不是原因。真實(shí)的物理原因需要進(jìn)一步采用動(dòng)理學(xué)的公式追溯流場(chǎng)結(jié)構(gòu)產(chǎn)生的因果關(guān)系。

    現(xiàn)在采用渦力理論進(jìn)行分析。圖10 顯示了近壁區(qū)域 (Region I+II)與脫體渦區(qū)域 (Region III)在不同攻角下升力隨著流向的變化曲線。圖中每個(gè)區(qū)域的升力采用渦力公式(17)計(jì)算,三個(gè)區(qū)域相加得到的升力與采用壁面積分得到的升力曲線幾乎完全一致,充分驗(yàn)證了公式的正確性和計(jì)算結(jié)果的精確性。

    (a)α=4° (b)α=20°

    α=4°時(shí),升力系數(shù)一直沿流向持續(xù)增加,尾緣部分的增速只稍微降低。且由于不存在脫體渦,近壁區(qū)域提供了全部升力。α=20°時(shí),升力系數(shù)在到達(dá)尾緣之前一直持續(xù)增加,但在尾緣部分其增速幾乎為0。與小攻角不同,脫體渦也提供了一部分升力。x較小時(shí),前緣渦與近壁區(qū)比較接近,兩個(gè)區(qū)域?qū)狭Φ呢暙I(xiàn)不能完全區(qū)分開來。但總體來說,近壁區(qū)提供了80%左右的升力,而脫體渦提供了 20%左右的升力。隨x增大,前緣渦不斷加強(qiáng),提供的升力逐漸線性化并占總升力的更高比例。到了尾緣附近,前緣渦提供的升力仍然穩(wěn)定增長(zhǎng),而近壁區(qū)提供的升力則略有下降。最終前緣渦給整個(gè)機(jī)翼提供了大約 40%的升力,近壁區(qū)提供了大約 60%的升力。這與 Yang等[29]的結(jié)論定性上一致。

    除了沿流向引入不同的橫流切片外,還可以沿展向引入不同的流向切片觀察升力的展向分布,如將機(jī)翼分為寬度為y=0.05的條帶,對(duì)每條條帶進(jìn)行積分。由于各條帶的面積不同,內(nèi)翼面區(qū)域顯然積分面積更大,為避免誤解,可以考慮單位面積的升力系數(shù)Cl(y)/S(y)隨y的變化曲線,如圖11所示。

    (a)α=4° (b)α=20°

    α=4°時(shí),單位面積產(chǎn)生的渦升力隨著|y|的增大而增大,在外翼面區(qū)域達(dá)到最大值。但當(dāng)α=20°時(shí),Cl(y)/S(y)隨著|y|的增大先增加后減小,峰值在|y|≈0.1附近,靠近前緣渦的渦心位置??倖挝幻娣e渦升力曲線的變化來源于近壁區(qū)與脫體渦產(chǎn)生的渦升力之間的競(jìng)爭(zhēng)。脫體渦產(chǎn)生的Cl(y)/S(y)隨著|y|的增大先保持不變,然后急劇增大,并且在外翼面達(dá)到最大值。曲線急劇增長(zhǎng)的位置與脫體渦的位置是相關(guān)的。而近壁區(qū)產(chǎn)生的Cl(y)/S(y)隨著|y|的增大先輕微變大,然后迅速變小。其曲線變化的位置與脫體渦曲線的變化位置是一致的。

    通過解剖兩個(gè)特征攻角下的流動(dòng)結(jié)構(gòu)對(duì)升力系數(shù)貢獻(xiàn)沿流向和展向的變化,可以發(fā)現(xiàn)造成三角翼升阻比下降的主要流動(dòng)結(jié)構(gòu)位于尾外翼面處。由于在兩個(gè)算例中機(jī)翼下表面均未發(fā)生分離,可以斷定上翼面的結(jié)構(gòu)是產(chǎn)生升阻比下降的主要原因。

    為了進(jìn)一步鎖定該流動(dòng)結(jié)構(gòu),還需考慮尾流截面型的氣動(dòng)力公式。公式(12)中第一項(xiàng)線性升力L0可轉(zhuǎn)換為尾流截面積分形式:

    (22)

    圖12是L0的被積函數(shù)在靠近機(jī)翼尾緣的橫截面上分布的云圖。由圖可見,上表面的|y|≈0.2處,正是主渦誘導(dǎo)的二次渦結(jié)構(gòu)在近壁區(qū)產(chǎn)生負(fù)升力的地方,與渦力理論的結(jié)果一致。正是由于該結(jié)構(gòu)的存在,使得近壁區(qū)的升力不再繼續(xù)增長(zhǎng)。

    圖12 α=20°,細(xì)長(zhǎng)三角翼x/c=0.9處線性升力L0云圖

    圖13展示了α=20°細(xì)長(zhǎng)三角翼上表面BVF型公式(18)給出的阻力和升力被積函數(shù)云圖。顯然,二次渦是在壁面上產(chǎn)生負(fù)升力的機(jī)制,而該二次渦也同樣產(chǎn)生了大量阻力。圖14進(jìn)一步給出代表型阻的總壓損失沿展向變化。顯然,二次渦不僅產(chǎn)生負(fù)升力還造成總壓損失的大峰值。

    (a)-ex·(x×σp)/2

    (b)ex·(x×σp)/2

    圖14 型阻總壓損失C(P∞-P)在x/c=0.6截面沿機(jī)翼邊界的分布

    3.3 基于空氣動(dòng)力學(xué)CT的概念性設(shè)計(jì)

    若僅從線性的設(shè)計(jì)思路出發(fā),要增升減阻只需簡(jiǎn)單的切掉尾緣附近的翼尖,使得二次渦沒有產(chǎn)生的機(jī)會(huì)。但對(duì)于大雷諾數(shù)下大攻角分離流來說,壁面上壓強(qiáng)產(chǎn)生的力占主導(dǎo),此時(shí)對(duì)于任意形狀的平板可以很簡(jiǎn)單的證明如下關(guān)系式:

    (23)

    因此線性設(shè)計(jì)思路無法幫助我們提高升阻比。

    (a)Case 4

    (b)Case 5

    (c)Case 6

    圖16 不同曲面三角翼截面曲線

    Case 4為圓弧曲線,Case 5采用曲率更大的兩個(gè)半圓弧拼接,Case 6采用樣條曲線。仍按α=20°、Re=5×105的流動(dòng)參數(shù)計(jì)算。計(jì)算的升阻力系數(shù)以及升阻比結(jié)果如表1所示。

    Case 4的升力增加,阻力有少許增加,升阻比增加了一些,但變化不大,說明這個(gè)方法可能有效。Case 5的升力減小了14%,但阻力減少了34%,升阻比提高了31%,這是由于其沿展向的大曲率基本上把將二次渦完全去掉了。這個(gè)改進(jìn)太過激進(jìn)。我們吸收了 Case 4和Case 5的優(yōu)點(diǎn),重新設(shè)計(jì)了 Case 6的曲線??梢钥闯?,升力基本沒變,但阻力減小了10%,從而使得升阻比也增加了10%。

    表1 不同曲面三角翼氣動(dòng)力系數(shù)表

    三個(gè)算例的升力基本沒變,且σp的貢獻(xiàn)遠(yuǎn)大于σvis,因此可以采用邊界型渦力公式來追蹤σp以闡述阻力減小導(dǎo)致升阻比增加的原因。圖15同時(shí)展示了邊界型公式中阻力被積函數(shù)-ex·(x×σp)/2的云圖。與圖13(a)對(duì)比,Case 4的阻力分布變化不大,在尾緣處仍有很大的阻力貢獻(xiàn)。雖然機(jī)翼前部和中部的阻力減少,但總阻力基本不變。Case 5中的阻力明顯減小,只在尾緣部分區(qū)域中出現(xiàn)阻力貢獻(xiàn),在尾緣外翼部基本沒有阻力。在翼中部還出現(xiàn)了明顯的推力,這說明該設(shè)計(jì)明顯減少了σp,從而減少了對(duì)阻力的貢獻(xiàn)。美中不足的是σp減少的太多,也導(dǎo)致了升力的減小。Case 6中阻力與圖13(a)對(duì)比,在尾緣處的阻力貢獻(xiàn)區(qū)域有所減小,在翼型中段的阻力減少得更多,而且分布更加均勻。這也證明了這種設(shè)計(jì)能通過減少壁面切向壓力梯度沿阻力方向的投影,進(jìn)而減少阻力。三種翼型下表面均是附著流動(dòng),其壓力以及壓力梯度分布均可用勢(shì)流解來估計(jì),與原方案的阻力系數(shù)差別不大。

    當(dāng)然,以上探索還只是基于對(duì)機(jī)翼和流場(chǎng)的初步認(rèn)識(shí),尚未提供優(yōu)化設(shè)計(jì)的系統(tǒng)理論方案。這也是今后需要進(jìn)一步探究和深化的方向。

    4 結(jié) 論

    傳統(tǒng)空氣動(dòng)力學(xué)流場(chǎng)診斷和減阻方案往往采用線性分解、各個(gè)擊破的思路。Zou等[4]發(fā)現(xiàn)上述思路并不能有效地處理真實(shí)復(fù)雜流動(dòng)。我們發(fā)現(xiàn)物體定常運(yùn)動(dòng)產(chǎn)生的Lamb矢量場(chǎng)ω×u及其演化,是制約所有氣動(dòng)力分量的關(guān)鍵物理載體。Lamb矢量可自然的分解為縱橫兩個(gè)部分。橫Lamb矢量體積分的x分量(即誘導(dǎo)阻力的橫分量)完全提供了壁面摩擦力;縱Lamb矢量體積分的x分量(即誘導(dǎo)阻力的縱分量)提供了壓差阻力和型阻之差。這意味著采用線性思維,單獨(dú)減少特定阻力分量的思路,會(huì)不可避免的影響其余氣動(dòng)力分量,很難達(dá)到優(yōu)化升阻比的目的。

    本文突破傳統(tǒng)阻力分解框架,提出面向任意復(fù)雜構(gòu)型的空氣動(dòng)力學(xué)CT原理,允許截面切入物體,給出了可切入物體內(nèi)部的渦力公式和尾流截面公式。通過對(duì)大攻角三角翼的流場(chǎng)進(jìn)行斷層掃描,找到了產(chǎn)生負(fù)升力和型阻相對(duì)應(yīng)的主要渦結(jié)構(gòu),并可以追溯到特定壁面區(qū)域的邊界渦量流?;谝陨险J(rèn)識(shí), 提出了概念性的增加升阻比的設(shè)計(jì)思路與參考構(gòu)型,并證實(shí)了該診斷原理的適用性和合理性。

    這個(gè)新的設(shè)計(jì)理念,對(duì)于復(fù)雜構(gòu)型流場(chǎng)診斷理論如何在實(shí)際工程應(yīng)用中發(fā)揮作用,也有著重要的啟發(fā)作用。常用的基于壓力速度層次的流場(chǎng)優(yōu)化策略,雖然方便易行,還能與EFD有效結(jié)合,但往往無法揭示物理根源;基于結(jié)構(gòu)層次的流場(chǎng)優(yōu)化策略,對(duì)CFD來說是合適的,卻由于EFD很難獲得三維流場(chǎng)而無法進(jìn)行驗(yàn)證??磥恚枰獙讉€(gè)層次的理論有機(jī)結(jié)合起來,對(duì)CFD數(shù)據(jù)進(jìn)行有針對(duì)性的挖掘,鎖定具體的影響氣動(dòng)力的關(guān)鍵流場(chǎng)結(jié)構(gòu)與壁面物理來源。這個(gè)認(rèn)識(shí)也能用于指導(dǎo)EFD對(duì)壁面和流場(chǎng)關(guān)鍵區(qū)域的測(cè)量,反過來佐證CFD計(jì)算和理論的正確性,并在此基礎(chǔ)上進(jìn)行下一步的優(yōu)化設(shè)計(jì)。這就不僅是建立CFD或EFD與理論的單向連接,而是建立起CFD/EFD/理論的相互雙向互動(dòng),充分發(fā)揮流體力學(xué)研究中每種手段的優(yōu)勢(shì),共同進(jìn)行流場(chǎng)診斷與優(yōu)化。

    應(yīng)當(dāng)看到,以近場(chǎng) Lamb矢量為主攻對(duì)象的空氣動(dòng)力學(xué) CT 技術(shù)才剛剛起步,還面臨新的挑戰(zhàn)。以定常不可壓縮流動(dòng)為例,氣動(dòng)力各分量不只取決于尾流中的分布渦量,而是由渦量場(chǎng)及其誘導(dǎo)的速度場(chǎng)共同決定,后者是前者的Biot-Savart積分, 有非局域效應(yīng)。Lamb矢量場(chǎng)既有散度,又有旋度,比渦量場(chǎng)內(nèi)容豐富復(fù)雜。其線性部分導(dǎo)致的L0和DP0公式僅提供了指導(dǎo)設(shè)計(jì)的一些基本線索,還需仔細(xì)研究其非線性部分的作用。對(duì)不同幾何構(gòu)型與流動(dòng)條件, 有旋有散的非線性Lamb矢量場(chǎng)及其矩的各分量有何特殊分布、相互如何制約、如何優(yōu)化, 是個(gè)嶄新課題。

    致謝:非常感謝張錫金教授和張淼博士啟發(fā)我們進(jìn)行此項(xiàng)研究,感謝史一蓬教授、蘇衛(wèi)東教授、白鵬研究員、李周復(fù)研究員、錢戰(zhàn)森研究員、張日葵博士、劉羅勤博士、高安康以及康林林對(duì)本研究極為有益的討論與建議。

    猜你喜歡
    尾流氣動(dòng)力升力
    高速列車車頂–升力翼組合體氣動(dòng)特性
    飛行載荷外部氣動(dòng)力的二次規(guī)劃等效映射方法
    無人機(jī)升力測(cè)試裝置設(shè)計(jì)及誤差因素分析
    基于自適應(yīng)偽譜法的升力式飛行器火星進(jìn)入段快速軌跡優(yōu)化
    側(cè)風(fēng)對(duì)拍動(dòng)翅氣動(dòng)力的影響
    飛機(jī)尾流的散射特性與探測(cè)技術(shù)綜述
    升力式再入飛行器體襟翼姿態(tài)控制方法
    錐形流量計(jì)尾流流場(chǎng)分析
    水面艦船風(fēng)尾流效應(yīng)減弱的模擬研究
    高速鐵路接觸線覆冰后氣動(dòng)力特性的風(fēng)洞試驗(yàn)研究
    男人舔奶头视频| 国内精品美女久久久久久| 午夜福利在线观看吧| 在线播放无遮挡| 美女xxoo啪啪120秒动态图| 欧美高清性xxxxhd video| 国内揄拍国产精品人妻在线| 又爽又黄无遮挡网站| 国产真实伦视频高清在线观看| 国产亚洲一区二区精品| 2021天堂中文幕一二区在线观| 国产麻豆成人av免费视频| 亚洲av国产av综合av卡| 国产人妻一区二区三区在| 国产乱人视频| 熟女人妻精品中文字幕| 亚洲精品视频女| 久久精品久久久久久久性| 深夜a级毛片| 成人亚洲精品av一区二区| 久久精品久久久久久噜噜老黄| 亚洲精品乱码久久久v下载方式| 日韩制服骚丝袜av| 午夜久久久久精精品| 黑人高潮一二区| 秋霞伦理黄片| av在线天堂中文字幕| 亚洲在久久综合| 在线观看美女被高潮喷水网站| 亚洲,欧美,日韩| 久久韩国三级中文字幕| 欧美日韩综合久久久久久| 欧美成人一区二区免费高清观看| 精品熟女少妇av免费看| 久99久视频精品免费| 亚洲精品久久久久久婷婷小说| 国产单亲对白刺激| 久久99热这里只有精品18| 亚洲人成网站高清观看| 色哟哟·www| 国产午夜精品论理片| 一个人观看的视频www高清免费观看| 亚洲成人av在线免费| 亚洲四区av| 国产极品天堂在线| 亚洲av不卡在线观看| 久久久久久伊人网av| 欧美潮喷喷水| 免费看美女性在线毛片视频| 久久久久久久久久久免费av| 免费观看的影片在线观看| 日韩av不卡免费在线播放| 超碰97精品在线观看| 国产成人91sexporn| 国产欧美日韩精品一区二区| 成人漫画全彩无遮挡| 一个人看的www免费观看视频| 日本欧美国产在线视频| 22中文网久久字幕| 伦精品一区二区三区| 69人妻影院| 夜夜看夜夜爽夜夜摸| 99热6这里只有精品| 搡女人真爽免费视频火全软件| 免费少妇av软件| av在线天堂中文字幕| 老司机影院毛片| 亚洲欧美中文字幕日韩二区| 亚洲国产精品sss在线观看| 91精品国产九色| 寂寞人妻少妇视频99o| 久久精品久久精品一区二区三区| 麻豆成人午夜福利视频| 午夜福利视频精品| a级毛片免费高清观看在线播放| 久久久成人免费电影| 国内精品美女久久久久久| 免费黄频网站在线观看国产| 欧美日韩视频高清一区二区三区二| 在线天堂最新版资源| 国产精品一二三区在线看| 亚洲国产精品成人久久小说| 国产单亲对白刺激| 最近中文字幕2019免费版| 一二三四中文在线观看免费高清| 美女xxoo啪啪120秒动态图| 插逼视频在线观看| 日本黄色片子视频| 亚洲av成人av| 黄色日韩在线| 麻豆av噜噜一区二区三区| 亚洲自偷自拍三级| 欧美不卡视频在线免费观看| 国产黄色视频一区二区在线观看| 菩萨蛮人人尽说江南好唐韦庄| 亚洲欧美成人精品一区二区| 亚洲欧美清纯卡通| 国产不卡一卡二| 国产成人午夜福利电影在线观看| 97热精品久久久久久| 国产午夜精品久久久久久一区二区三区| 亚洲精品一区蜜桃| 免费高清在线观看视频在线观看| 亚洲av在线观看美女高潮| freevideosex欧美| 国产美女午夜福利| 黄片wwwwww| 看黄色毛片网站| 午夜福利在线观看吧| 久久99蜜桃精品久久| 国产探花在线观看一区二区| 亚洲高清免费不卡视频| 免费高清在线观看视频在线观看| 成人美女网站在线观看视频| 国产av国产精品国产| 亚洲精品乱码久久久久久按摩| 久久精品国产鲁丝片午夜精品| 在线观看美女被高潮喷水网站| 欧美精品一区二区大全| 亚洲精品国产av蜜桃| 亚洲av成人av| 国产有黄有色有爽视频| 中文字幕免费在线视频6| 精品一区在线观看国产| 国产av码专区亚洲av| 在线 av 中文字幕| 亚洲美女视频黄频| 国产伦在线观看视频一区| 熟女人妻精品中文字幕| 亚洲精品乱久久久久久| 三级经典国产精品| 99久国产av精品| 嫩草影院精品99| 国产高潮美女av| 精品人妻一区二区三区麻豆| 最近中文字幕高清免费大全6| 日韩一本色道免费dvd| 免费黄频网站在线观看国产| 国产精品久久久久久av不卡| 亚洲国产精品sss在线观看| 久久久久久久久中文| 老司机影院成人| 一级毛片 在线播放| 网址你懂的国产日韩在线| 国产精品蜜桃在线观看| 嫩草影院精品99| 亚洲av中文av极速乱| av专区在线播放| 久久久久网色| 国产精品国产三级国产av玫瑰| 免费av不卡在线播放| 亚洲av电影不卡..在线观看| 校园人妻丝袜中文字幕| 国产黄频视频在线观看| 国产成人午夜福利电影在线观看| 精品人妻熟女av久视频| 一边亲一边摸免费视频| 我的老师免费观看完整版| 亚洲在线自拍视频| 黄色一级大片看看| 白带黄色成豆腐渣| 国产91av在线免费观看| 观看美女的网站| 亚洲真实伦在线观看| 18+在线观看网站| 99热全是精品| 噜噜噜噜噜久久久久久91| 日韩av不卡免费在线播放| 99久久精品一区二区三区| 国产精品美女特级片免费视频播放器| 少妇的逼好多水| 综合色av麻豆| videossex国产| 亚洲婷婷狠狠爱综合网| 国模一区二区三区四区视频| 欧美成人一区二区免费高清观看| 亚洲国产最新在线播放| or卡值多少钱| eeuss影院久久| 熟女电影av网| 天堂网av新在线| 最近2019中文字幕mv第一页| 国产高清国产精品国产三级 | 日韩在线高清观看一区二区三区| 少妇熟女欧美另类| 综合色丁香网| 日本wwww免费看| 毛片一级片免费看久久久久| 午夜精品国产一区二区电影 | 国产欧美日韩精品一区二区| 久久久久久久久久人人人人人人| 又大又黄又爽视频免费| 六月丁香七月| 免费黄网站久久成人精品| av福利片在线观看| 色综合亚洲欧美另类图片| 国产精品久久视频播放| 中文字幕免费在线视频6| www.av在线官网国产| 最新中文字幕久久久久| 亚洲精品影视一区二区三区av| 久久精品国产亚洲av天美| 国内精品一区二区在线观看| 免费高清在线观看视频在线观看| 2018国产大陆天天弄谢| 精品久久久久久久末码| 午夜精品国产一区二区电影 | 亚洲不卡免费看| 国产老妇女一区| 国产精品一及| 一二三四中文在线观看免费高清| 别揉我奶头 嗯啊视频| 性色avwww在线观看| 丝袜喷水一区| 一区二区三区高清视频在线| 国产视频首页在线观看| 麻豆av噜噜一区二区三区| 亚洲最大成人中文| 婷婷六月久久综合丁香| 欧美成人一区二区免费高清观看| 午夜精品一区二区三区免费看| 亚洲精品自拍成人| 精品一区在线观看国产| 亚洲成人久久爱视频| 国产午夜精品论理片| 亚洲欧洲国产日韩| 超碰av人人做人人爽久久| 麻豆乱淫一区二区| 免费看光身美女| 精品久久久久久电影网| 亚洲最大成人手机在线| 国产淫语在线视频| 国产黄频视频在线观看| 久久午夜福利片| 精品99又大又爽又粗少妇毛片| 久久久久久久久久人人人人人人| 亚洲av男天堂| 午夜爱爱视频在线播放| 不卡视频在线观看欧美| 3wmmmm亚洲av在线观看| 欧美xxxx黑人xx丫x性爽| av黄色大香蕉| 午夜福利高清视频| 超碰av人人做人人爽久久| 免费黄网站久久成人精品| 久久亚洲国产成人精品v| 亚洲精品亚洲一区二区| 99九九线精品视频在线观看视频| 特大巨黑吊av在线直播| 免费无遮挡裸体视频| 日韩不卡一区二区三区视频在线| 亚洲国产精品成人综合色| 精品99又大又爽又粗少妇毛片| 亚洲国产高清在线一区二区三| 久久精品综合一区二区三区| 女人被狂操c到高潮| 国产精品.久久久| 丰满人妻一区二区三区视频av| 能在线免费看毛片的网站| 日本三级黄在线观看| 人人妻人人看人人澡| 久久久久久久久久人人人人人人| 亚洲怡红院男人天堂| 国语对白做爰xxxⅹ性视频网站| 欧美 日韩 精品 国产| 女人久久www免费人成看片| 热99在线观看视频| 日韩亚洲欧美综合| 日本一本二区三区精品| 99久久九九国产精品国产免费| 亚洲欧美成人综合另类久久久| av卡一久久| 国产91av在线免费观看| 欧美三级亚洲精品| 精品久久久久久久久av| 国产精品蜜桃在线观看| 国产高清不卡午夜福利| 又粗又硬又长又爽又黄的视频| 国产中年淑女户外野战色| 国产探花极品一区二区| 日韩制服骚丝袜av| 男人和女人高潮做爰伦理| 亚洲在线观看片| freevideosex欧美| 国产黄片视频在线免费观看| 国产在视频线在精品| 1000部很黄的大片| 一夜夜www| 日韩欧美国产在线观看| 天堂中文最新版在线下载 | 乱系列少妇在线播放| 天堂av国产一区二区熟女人妻| 五月伊人婷婷丁香| 午夜福利视频1000在线观看| 日韩一区二区视频免费看| 国产午夜精品久久久久久一区二区三区| 男女视频在线观看网站免费| 国产日韩欧美在线精品| 欧美精品一区二区大全| 偷拍熟女少妇极品色| 日韩一区二区三区影片| 成年版毛片免费区| 日韩欧美精品v在线| 有码 亚洲区| 国模一区二区三区四区视频| 亚洲一区高清亚洲精品| 日韩欧美国产在线观看| 搞女人的毛片| 天天躁夜夜躁狠狠久久av| 色吧在线观看| 久久6这里有精品| 一级黄片播放器| 又黄又爽又刺激的免费视频.| 国产伦一二天堂av在线观看| 欧美人与善性xxx| 91久久精品电影网| 黄片wwwwww| 综合色av麻豆| 大香蕉97超碰在线| 青春草国产在线视频| 亚洲天堂国产精品一区在线| 欧美日本视频| 国产免费福利视频在线观看| 九九久久精品国产亚洲av麻豆| videossex国产| 欧美97在线视频| 2022亚洲国产成人精品| 国模一区二区三区四区视频| 久久久欧美国产精品| 亚洲欧美精品专区久久| 一区二区三区乱码不卡18| 亚洲丝袜综合中文字幕| 日韩伦理黄色片| 一本久久精品| 国产有黄有色有爽视频| 婷婷色综合www| 免费av不卡在线播放| 两个人视频免费观看高清| 久久久久免费精品人妻一区二区| 激情 狠狠 欧美| 国产伦精品一区二区三区四那| 亚洲18禁久久av| 老女人水多毛片| 少妇被粗大猛烈的视频| 成年人午夜在线观看视频 | 91精品伊人久久大香线蕉| 国产v大片淫在线免费观看| 中文在线观看免费www的网站| h日本视频在线播放| 男女下面进入的视频免费午夜| 99久国产av精品国产电影| 国产成人一区二区在线| 在线免费观看的www视频| 亚洲成人久久爱视频| 国产成人freesex在线| 久久这里只有精品中国| 国产人妻一区二区三区在| 久久久久国产网址| 国产欧美另类精品又又久久亚洲欧美| 日韩一区二区三区影片| 美女脱内裤让男人舔精品视频| av在线天堂中文字幕| videos熟女内射| 亚洲av中文av极速乱| 午夜福利在线观看吧| 80岁老熟妇乱子伦牲交| 菩萨蛮人人尽说江南好唐韦庄| 国产黄色视频一区二区在线观看| 网址你懂的国产日韩在线| 丰满乱子伦码专区| 内地一区二区视频在线| 大话2 男鬼变身卡| 搡老妇女老女人老熟妇| 免费人成在线观看视频色| 男人舔女人下体高潮全视频| 在线免费观看不下载黄p国产| 国产精品一区二区三区四区免费观看| 91久久精品国产一区二区三区| 中文字幕久久专区| 国产精品av视频在线免费观看| 午夜精品在线福利| 亚洲成色77777| 嫩草影院入口| 全区人妻精品视频| 亚洲最大成人手机在线| 日本黄大片高清| 在线 av 中文字幕| 日韩一本色道免费dvd| 最近2019中文字幕mv第一页| 男人舔奶头视频| 少妇熟女欧美另类| 精品人妻视频免费看| 嘟嘟电影网在线观看| 美女大奶头视频| 亚洲精品中文字幕在线视频 | 亚洲精品影视一区二区三区av| 高清视频免费观看一区二区 | 国产精品一区二区在线观看99 | 国产成人一区二区在线| 午夜亚洲福利在线播放| 大香蕉久久网| 六月丁香七月| 天堂中文最新版在线下载 | 午夜爱爱视频在线播放| freevideosex欧美| 亚洲精品456在线播放app| 亚洲熟女精品中文字幕| 精品99又大又爽又粗少妇毛片| 亚洲av成人精品一二三区| 别揉我奶头 嗯啊视频| 亚洲综合色惰| 精品一区二区三卡| 免费av观看视频| 91久久精品国产一区二区三区| 亚洲精华国产精华液的使用体验| 免费看美女性在线毛片视频| 国产精品一区二区三区四区免费观看| av国产久精品久网站免费入址| 欧美区成人在线视频| 高清欧美精品videossex| 日本黄大片高清| 久99久视频精品免费| 亚洲av国产av综合av卡| 91精品国产九色| 又大又黄又爽视频免费| 国产色爽女视频免费观看| 国产成人精品婷婷| 赤兔流量卡办理| 偷拍熟女少妇极品色| 成人高潮视频无遮挡免费网站| 嫩草影院精品99| 不卡视频在线观看欧美| 成人无遮挡网站| 永久网站在线| 97精品久久久久久久久久精品| 国产亚洲5aaaaa淫片| 婷婷色av中文字幕| 精品久久久久久电影网| 欧美97在线视频| 精品一区在线观看国产| 成人高潮视频无遮挡免费网站| 大片免费播放器 马上看| 欧美区成人在线视频| 好男人在线观看高清免费视频| 久久久亚洲精品成人影院| 99热全是精品| 最新中文字幕久久久久| 国产成人a区在线观看| 亚洲av福利一区| 成人午夜精彩视频在线观看| 久久久国产一区二区| 国产色婷婷99| 我的女老师完整版在线观看| 久久精品国产亚洲av天美| 国内揄拍国产精品人妻在线| 搡女人真爽免费视频火全软件| 丰满少妇做爰视频| 国产成人精品福利久久| av在线天堂中文字幕| 午夜精品在线福利| 亚洲av在线观看美女高潮| 欧美日韩在线观看h| 日韩精品有码人妻一区| 国产黄色视频一区二区在线观看| 国产精品久久久久久久久免| 亚洲精品自拍成人| 午夜激情福利司机影院| 中文字幕av成人在线电影| 99久国产av精品国产电影| 99视频精品全部免费 在线| 午夜老司机福利剧场| 99久久九九国产精品国产免费| 亚洲精品国产成人久久av| 午夜久久久久精精品| 五月伊人婷婷丁香| 日本黄色片子视频| 神马国产精品三级电影在线观看| 老司机影院成人| 亚洲成人久久爱视频| 国产免费福利视频在线观看| 观看美女的网站| 看黄色毛片网站| 天堂√8在线中文| 国产成人a区在线观看| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 国产免费福利视频在线观看| 卡戴珊不雅视频在线播放| 亚洲性久久影院| 中国国产av一级| 男插女下体视频免费在线播放| 日本熟妇午夜| 日韩成人av中文字幕在线观看| 亚洲精品一二三| 国产视频首页在线观看| 少妇的逼水好多| 国产成人精品婷婷| 精品国内亚洲2022精品成人| 日韩三级伦理在线观看| 国产色爽女视频免费观看| 高清日韩中文字幕在线| 精品久久久久久久久久久久久| 日韩精品有码人妻一区| 特大巨黑吊av在线直播| 亚洲欧美精品专区久久| 国产精品久久久久久精品电影小说 | 菩萨蛮人人尽说江南好唐韦庄| 欧美xxxx性猛交bbbb| 国产成人freesex在线| av福利片在线观看| 内地一区二区视频在线| 我的老师免费观看完整版| 嫩草影院精品99| 日日啪夜夜爽| 日本黄大片高清| 嫩草影院入口| 欧美高清成人免费视频www| 久久这里只有精品中国| 国内精品宾馆在线| 色综合色国产| 国产精品爽爽va在线观看网站| 十八禁国产超污无遮挡网站| 搡女人真爽免费视频火全软件| 久99久视频精品免费| 99视频精品全部免费 在线| 超碰av人人做人人爽久久| 最近最新中文字幕大全电影3| 五月伊人婷婷丁香| 亚洲精品成人av观看孕妇| 久久亚洲国产成人精品v| 国产不卡一卡二| 色视频www国产| 欧美三级亚洲精品| 久久韩国三级中文字幕| 高清视频免费观看一区二区 | 亚洲aⅴ乱码一区二区在线播放| 天堂网av新在线| 亚洲自拍偷在线| 晚上一个人看的免费电影| 亚洲av男天堂| 天堂影院成人在线观看| 黄片无遮挡物在线观看| 久久精品久久久久久久性| av国产免费在线观看| 亚洲精品一区蜜桃| 日韩三级伦理在线观看| 久久久久久久大尺度免费视频| freevideosex欧美| 亚洲精品中文字幕在线视频 | 成人亚洲精品一区在线观看 | 精品午夜福利在线看| 国产色爽女视频免费观看| 国产av不卡久久| 精品久久久久久久人妻蜜臀av| 综合色丁香网| 国产男女超爽视频在线观看| 亚洲欧美中文字幕日韩二区| 2022亚洲国产成人精品| 亚洲av不卡在线观看| 色综合亚洲欧美另类图片| 联通29元200g的流量卡| 青春草亚洲视频在线观看| 大香蕉久久网| 禁无遮挡网站| 中文字幕av成人在线电影| 国产淫片久久久久久久久| 亚洲自拍偷在线| 尤物成人国产欧美一区二区三区| 久久人人爽人人爽人人片va| 在线播放无遮挡| 日本av手机在线免费观看| 少妇高潮的动态图| 亚洲va在线va天堂va国产| 午夜福利在线观看免费完整高清在| 日韩亚洲欧美综合| 婷婷色麻豆天堂久久| 国产麻豆成人av免费视频| 亚洲成人精品中文字幕电影| 九九爱精品视频在线观看| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 国产淫片久久久久久久久| 亚洲精品视频女| 国产精品久久久久久久久免| 菩萨蛮人人尽说江南好唐韦庄| 成人综合一区亚洲| 嫩草影院精品99| 国产成人精品福利久久| 美女cb高潮喷水在线观看| 男女啪啪激烈高潮av片| 啦啦啦啦在线视频资源| 国产伦精品一区二区三区视频9| 天堂网av新在线| 青春草亚洲视频在线观看| 色综合色国产| 亚洲四区av| 国产av在哪里看| 亚洲精品乱码久久久v下载方式| 亚洲av福利一区| 国产午夜精品论理片| 亚洲精品影视一区二区三区av| 国产免费一级a男人的天堂| 欧美丝袜亚洲另类| 欧美日韩精品成人综合77777| 免费黄频网站在线观看国产| 亚洲最大成人av| 99久久九九国产精品国产免费| 久久国内精品自在自线图片| 三级国产精品欧美在线观看| 99热这里只有是精品在线观看| 又大又黄又爽视频免费| 精品人妻一区二区三区麻豆| 国产成人免费观看mmmm| 色尼玛亚洲综合影院| 高清日韩中文字幕在线| 最近最新中文字幕大全电影3| 亚州av有码|