柔性仿羽毛結(jié)構(gòu)抑制邊界層轉(zhuǎn)捩的初步探索

葉正寅, 洪 正, 武 潔

(西北工業(yè)大學(xué), 西安 710072)

0 引 言

對于常規(guī)民用飛機(jī)而言，阻力增加1%就會在一次飛行中額外多消耗3 T的燃油，在相同起飛重量的情況下，會縮短120 km的航程[1]；對于相同起飛距離，升阻比提高1%，意味著可以增加14個乘客[2]。在構(gòu)成民用飛機(jī)阻力的成分中，表面摩擦阻力占到了50%[3]。層流邊界層的摩擦阻力遠(yuǎn)小于湍流邊界層。研究表明[3]: 飛機(jī)機(jī)翼、平尾、垂尾和短艙分別貢獻(xiàn)了總阻力的18%、4%、3%和3%。如果層流邊界層在這些部件表面上能夠達(dá)到20%、30%甚至40%的比例，則飛機(jī)阻力分別能夠減少8%、12%甚至16%。因此，推遲邊界層的轉(zhuǎn)捩在飛行器減阻方面具有很大的實(shí)際價值。

在目前已知的邊界層控制手段中，有吹/吸氣控制[4]、表面溝槽控制[5-7]、振蕩射流控制[8-9]、振蕩洛倫茲力控制[10]、展向變形的行波法控制[11-13]、流向變形控制[14-15]、展向粗糙帶控制[16]、邊界層燃燒[17]等。當(dāng)然，在航空領(lǐng)域的工程應(yīng)用方面，目前最主要的途徑仍然是通過設(shè)計(jì)層流機(jī)翼的方法推遲轉(zhuǎn)捩的發(fā)生[18]。為了實(shí)現(xiàn)層流的有效控制，邊界層的控制方法是當(dāng)前的重點(diǎn)研究方向。但是，一些手段在向?qū)嶋H應(yīng)用推廣的過程中存在很大的難度和代價，如展向行波法、洛倫茲力控制技術(shù)[13]。

層流的被動控制方法具有成本低的特點(diǎn)，上述溝槽控制轉(zhuǎn)捩的方法、粗糙表面的方法都屬于被動控制措施。事實(shí)上，人們甚至在高速流動中，也仍然在探索有效的被動控制方法[19-20]。

是否存在更直接的邊界層控制思路？自然界是人類尋求靈感的源泉，正如文獻(xiàn)[21]中所講，人類在發(fā)展飛行器的歷史進(jìn)程中，從鳥類飛行中得到過多種啟示,鳥類的羽毛經(jīng)過千萬年的進(jìn)化，仍然蘊(yùn)藏著許多奧秘。Tucher等人[22]測量了哈里斯“鷹”翼羽修剪前后在風(fēng)洞中自由滑翔時的最小阻力，發(fā)現(xiàn)在7.3 m/s到15 m/s的速度之間，未修剪的“鷹”阻力只有修剪后“鷹”阻力的70%～90%。Parfitt等[23]發(fā)現(xiàn)，許多企鵝物種喙周圍突出的羽毛區(qū)域平均減少了31%的阻力。

從我們對鳥類羽毛的直接觀測中不難發(fā)現(xiàn)，當(dāng)氣流順羽毛方向時流動是順利的，可當(dāng)氣流從側(cè)面流過羽毛時，羽毛的側(cè)緣就會卷起來阻止當(dāng)?shù)卣瓜蛄鲃?由于羽毛側(cè)緣卷起時基本上沿流向排列，在來流方向的迎風(fēng)面積很小，而沿展向方向的迎風(fēng)面積大)。因此，對于流經(jīng)表面的氣流而言，當(dāng)局部氣流沿展向流動時，在展向方向阻力很大，流向阻力小得多。而且離壁面高度越低，阻力會越大，當(dāng)高度超過一定值(對應(yīng)羽毛側(cè)緣卷起來的高度)時，阻力就沒有了。也就是說，羽毛對近物面氣流具有各向異性的阻力特征。而羽毛的流固耦合特性就是羽毛產(chǎn)生這種各向異性阻力的“激發(fā)開關(guān)”。

通過對邊界層的大量研究[24-29]，人們發(fā)現(xiàn)存在著一種共性特點(diǎn)，即：邊界層從層流轉(zhuǎn)化為湍流的轉(zhuǎn)捩過程，必然出現(xiàn)Λ渦和發(fā)卡渦，而這兩種渦的“腿”就是一種傾斜的局部流向渦，它會誘導(dǎo)氣流的展向流動，如果羽毛規(guī)則的柔性邊緣因?yàn)檎瓜蛄鲃佣砥?，就會對局部展向流動有阻滯作用并抑制Λ渦和發(fā)卡渦的發(fā)展，從而抑制轉(zhuǎn)捩的發(fā)生。

如前面所提到，研究者們利用被動或主動的措施來減小邊界層的阻力，但大多是針對湍流區(qū)的減阻。然而，由于湍流邊界層的阻力遠(yuǎn)大于層流邊界層，減阻的效果沒有層流控制來得更直接有效。鳥類高效的飛行效率必然對應(yīng)著其羽毛表面有著極低的飛行阻力，從而可以推斷層流流動在羽毛表面是占主導(dǎo)的。那么，觀察到的羽毛表面的各向異性阻力特征是否對邊界層從層流到湍流的轉(zhuǎn)捩有抑制效果呢？在國家數(shù)值風(fēng)洞項(xiàng)目的資助下，作者對柔性羽毛進(jìn)行了流固耦合特性分析，旨在探索羽毛的這種柔性結(jié)構(gòu)特征在抑制邊界層轉(zhuǎn)捩中的效果。

1 羽毛模型

羽毛的結(jié)構(gòu)和排列方式非常復(fù)雜，建立等價于真實(shí)羽毛的模型是非常困難甚至不可能的。因此，抓住所關(guān)注的主要特征對模型進(jìn)行簡化是十分必要的。如引言中提到的，在本文中，羽毛的柔性體現(xiàn)在存在展向流動時羽毛側(cè)緣會卷起，從而阻礙該展向流動，而沒有展向流動時羽毛則不會卷起即不起作用。此外，考慮到真實(shí)羽毛的特點(diǎn)，展向流動應(yīng)當(dāng)可以穿過羽毛，即羽毛不能完全抑制展向流動。對于羽毛表面的這種各向異性的阻力特征，本文采用體積力的形式對其進(jìn)行建模。為簡單起見，羽毛側(cè)緣卷起時對物面法向的影響忽略不計(jì)。

1.1 展向阻力模型

對于卷起的羽毛側(cè)緣給當(dāng)?shù)卣瓜蛄鲃訋淼淖枇?，有以下幾點(diǎn)考慮：

(1)展向脈動速度越大，那么這種阻力也越大；

(2)越靠近表面，羽毛越密，則阻力也隨之增大；

(3)這種阻力只能使得當(dāng)?shù)卣瓜蛎}動速度減小，但不能改變其方向。

羽毛阻礙展向脈動流動產(chǎn)生的阻力，概念上與黏性阻力類似，則基于(1)、(2)兩點(diǎn)且類比黏性力的公式，可將沿展向的阻力?；癁椋?/p>

(1)

為了確定σy的值，本文假設(shè)在離表面最近、高度為h0的控制點(diǎn)處，即展向脈動抑制作用最強(qiáng)的位置，展向脈動在一個計(jì)算時間步內(nèi)恰好被衰減到0。于是，在該位置處利用動量定理近似有

(ρfy)h=h0VΔt=0-ρv′V

(2)

式中：V是離散控制體體積，Δt是計(jì)算的時間步長。將式(1)代入式(2)，可得:

(3)

更一般的，在任意高度h處，根據(jù)動量定理整理可得:

(4)

式中：下標(biāo)tn、tn+1分別指代當(dāng)前時刻和下一時刻。式(4)表明，使用式(3)來確定σy，則任意高度處的展向脈動速度相對衰減率只與高度成反比，與速度大小無關(guān)。

假設(shè)有一均勻流動的展向脈動速度固定為v′=1，網(wǎng)格和參數(shù)采用下文中的設(shè)置，但壁面處理為滑移壁面。這樣，流動的變化僅受到模型公式(1)的影響。圖1給出了1、10、100時間步后任一位置處v′沿壁面高度的分布情況。從圖1中可以看到，模型對展向流動的衰減作用,越靠近壁面越強(qiáng)。即使在離壁面稍遠(yuǎn)的位置處，多次的衰減也將使展向流動趨于滯止。顯然，模型只能削減速度的大小直到v′=0，然后模型不再起作用，速度方向也不會因此改變。

圖1 展向脈動在模型作用下的衰減特征

1.2 流向阻力模型

羽毛卷起的側(cè)緣除了主要對展向脈動流動起抑制作用外，也會對沿流向的流動產(chǎn)生或多或少的阻礙。因此，從完整性的角度考慮，也需要對羽毛卷起帶來的流向阻力進(jìn)行建模。

對于流向阻力的體積力模型，有以下幾點(diǎn)考慮：

(1)流向阻力的公式與展向形式保持一致；

(2)展向脈動流動越強(qiáng)，則卷起的特征越顯著，由此帶來的流向阻力也越大。

與展向阻力模型保持形式一致，則流向阻力的模型可寫為:

(5)

式中：fx為流向的體積力，σx是待確定的流向阻力系數(shù)，u是瞬時速度在流向的分量。與展向阻力系數(shù)σy相關(guān)聯(lián)，σx的確定按如下公式:

(6)

2 數(shù)值方法和計(jì)算設(shè)置

2.1 控制方程

使用守恒形式的三維N-S方程來描述流動，形式如下:

(7)

式中：Q=(ρ,ρu,ρv,ρw,ρe)T是守恒變量，v、w是瞬時速度在展向和壁面法向方向上的分量，e是氣體的總能，S=(0,ρfx,ρfy,0,uρfx+vρfy)T是方程的源項(xiàng)，羽毛的各向異性阻力模型正是在此處起作用。Fc、Gc、Hc是分別沿x、y、z方向的對流通量，F(xiàn)v、Gv、Hv則是相應(yīng)的黏性通量。

使用半離散的格心型有限差分方法來數(shù)值求解N-S方程。為求得控制點(diǎn)處的通量導(dǎo)數(shù)，先基于格心處已知的原始變量，通過五階DCS(Dissipative Compact Scheme)插值得到界面上的原始變量的左迎風(fēng)值和右迎風(fēng)值。然后，對于對流通量，根據(jù)得到的左右迎風(fēng)值，使用Roe格式計(jì)算界面處的對流通量。對于黏性通量，左右迎風(fēng)值取算數(shù)平均作為界面值，然后由通量公式得到。得到界面上的對流通量和黏性通量后，采用六階緊致中心格式來計(jì)算通量導(dǎo)數(shù)。最后，利用顯式的四步二階龍格-庫塔法來計(jì)算下一時刻格心處的流場變量值。五階DCS和格心型有限差分法的具體公式，可參閱文獻(xiàn)[30]，本文中DCS中控制耗散的參數(shù)取為0.5。

2.2 計(jì)算設(shè)置

計(jì)算域如圖2 所示。基本流為Blasius相似性解，在計(jì)算域入口處需引入擾動。初始擾動由3個T-S(Tollmien-Schlichting)波組成，具體形式如下

圖2 計(jì)算域示意圖

(8)

入口處擾動的參數(shù)選取參考文獻(xiàn)[31]中的設(shè)置，如表1所示。

表1 入口擾動參數(shù)

計(jì)算域的三維尺寸(以入口處邊界層的位移厚度無量綱化)和網(wǎng)格的參數(shù)如表2所示。展向和流向網(wǎng)格都是平均分布的，沿壁面法向，網(wǎng)格點(diǎn)按如下公式進(jìn)行拉伸。

(9)

表2 計(jì)算域尺寸和網(wǎng)格參數(shù)

3 結(jié)果分析與討論

首先，模擬了平板上的邊界層轉(zhuǎn)捩流動，一方面可以驗(yàn)證所采用的數(shù)值方法的可靠性，另一方面提供了一個沒有加入羽毛模型的基準(zhǔn)參考。然后，加入羽毛模型后再進(jìn)行平板邊界層轉(zhuǎn)捩的模擬，對比基準(zhǔn)結(jié)果就可以得到羽毛對邊界層轉(zhuǎn)捩的影響。

3.1 平板邊界層轉(zhuǎn)捩

盡管湍流是隨機(jī)且復(fù)雜的，但仍是有跡可循的。無量綱的壁面高度定義為：

(10)

(11)

圖3所示為無量綱的邊界層速度分布。在入口x=0處，速度分布是典型的層流剖面，而在x=300處，速度分布已經(jīng)演化為典型的湍流剖面。這說明使用的數(shù)值方法和計(jì)算設(shè)置可以反映層流到湍流轉(zhuǎn)捩的過程。

圖3 層流區(qū)和湍流區(qū)的速度剖面

壁面摩擦系數(shù)Cf的定義為：

(12)

對于層流，摩擦系數(shù)為

(13)

式中：δ是當(dāng)?shù)氐倪吔鐚游灰坪穸?。Ducros等人[33]提出了一種經(jīng)驗(yàn)公式來描述湍流邊界層摩擦系數(shù)的空間演化規(guī)律，具體形式為：

(14)

式中：l′=125δ0。圖4中給出了平板上轉(zhuǎn)捩流動，層流以及Ducros經(jīng)驗(yàn)公式的壁面摩擦系數(shù)。從圖中可以看到，大約從x=80開始，平板邊界層的摩擦系數(shù)開始快速偏離層流對應(yīng)的值，該點(diǎn)可視為層流向湍流轉(zhuǎn)捩的起始點(diǎn)。摩擦系數(shù)在x=190達(dá)到峰值，之后流動進(jìn)入完全湍流狀態(tài)。湍流區(qū)的摩擦系數(shù)與經(jīng)驗(yàn)公式的結(jié)果所在的區(qū)間是一致的。

圖4 沿壁面的摩擦系數(shù)分布

平板邊界層從層流到湍流的演化過程中，存在一些典型的大尺度相干結(jié)構(gòu)。利用Q準(zhǔn)則，圖5呈現(xiàn)了轉(zhuǎn)捩過程中不同階段的渦結(jié)構(gòu)，用距離壁面的高度進(jìn)行了著色。為了更好地顯示效果，展向通過復(fù)制延伸了一倍?？梢钥吹睫D(zhuǎn)捩開始后先是形成Λ渦，然后是Λ渦演化成發(fā)卡渦，進(jìn)而生成環(huán)狀渦，并最終發(fā)展成環(huán)狀渦鏈結(jié)構(gòu)。對典型渦系結(jié)構(gòu)的清晰捕捉也體現(xiàn)了數(shù)值模擬的可靠性。

(a)t=1.6T0

3.2 羽毛模型對平板轉(zhuǎn)捩的影響

為了研究羽毛的各向異性阻力特征對轉(zhuǎn)捩的影響，將第1節(jié)中建立的羽毛阻力模型加入N-S方程，其余設(shè)置均保持相同。流向阻力模型中的控制參數(shù)m考慮三種情況：(1)m=∞，對應(yīng)羽毛帶來的流向阻力忽略不計(jì)的情形；(2)m=1.5，對應(yīng)流向阻力比展向阻力小的情形；(3)m=1，對應(yīng)流向和展向阻力大小相等的情形。實(shí)際上，羽毛卷起的側(cè)緣能夠影響的高度是有限的。根據(jù)前面的平板邊界層模擬結(jié)果，模型起作用的范圍限制在z

圖6給出了加入模型后計(jì)算得到的全流場的渦系結(jié)構(gòu)圖，上一小節(jié)中未加模型的結(jié)果也一并給出作為對比?？梢钥吹皆诩尤胗鹈淖枇δＰ秃?，開始形成大尺度渦的位置較無模型情況明顯靠后。無模型情況下，x=100附近就已經(jīng)形成了發(fā)卡渦，而加入模型且m=∞時，x=300處發(fā)卡渦才出現(xiàn)，并且轉(zhuǎn)捩到湍流過程中形成的渦系要稀疏得多。不同的m值得到的結(jié)果有顯著的差異。m=1.5與m=∞相比，開始形成典型渦結(jié)構(gòu)的位置相差不大，但后續(xù)過程中形成的渦系結(jié)構(gòu)要豐富一些。m=1情況下，開始形成渦的位置較m=∞明顯提前，x=200附近已經(jīng)可以看到形成的發(fā)卡渦了。此外，之后形成的渦系結(jié)構(gòu)相較之下也要豐富得多。

圖6 沿平板的渦系結(jié)構(gòu)(從上到下：無模型，m=∞,m=1.5,m=1)

圖7 展向脈動速度均方根沿平板的分布(從上到下：無模型，m=∞,m=1.5,m=1)

圖8 平均湍動能沿平板的分布(從上到下：無模型，m=∞,m=1.5,m=1)

圖9給出了這幾種情況下的壁面摩擦系數(shù)。加入羽毛模型后，即使在轉(zhuǎn)捩發(fā)生以后，壁面摩擦系數(shù)也沒有增大很多，保持在與層流對應(yīng)的摩擦系數(shù)相當(dāng)?shù)乃?。特別地，在m=1的情況下，轉(zhuǎn)捩之后的摩擦系數(shù)不但沒有增大，反而下降得比層流時還要小。湍流區(qū)的摩擦系數(shù)比層流區(qū)的摩擦系數(shù)還要小，這與通常的認(rèn)知是相悖的。實(shí)際上，在加入羽毛模型且不忽略流向阻力情況下，平板的阻力來源除了壁面的摩擦力外，還有羽毛阻礙流向流動產(chǎn)生的反作用力。因此，平板的總阻力系數(shù)應(yīng)當(dāng)寫為:

圖9 不同m值得到的摩擦系數(shù)分布

(15)

阻力系數(shù)的第I部分對應(yīng)的正是壁面摩擦系數(shù)，第II部分對應(yīng)的是羽毛阻礙流向流動帶來的反作用力系數(shù)。

m=∞對應(yīng)的是流向阻力小到可以忽略的情形，因此阻力僅來自壁面摩擦。對于m=1.5和m=1，圖10給出了總阻力系數(shù)及其兩部分系數(shù)的分布。轉(zhuǎn)捩開始后生成的大尺度渦結(jié)構(gòu)激發(fā)羽毛模型，可以看到，由此帶來的反作用力阻力開始迅速增長并在湍流區(qū)維持較高水平。從總阻力系數(shù)來看，轉(zhuǎn)捩之后阻力系數(shù)相比層流時顯著增大，這樣就與湍流區(qū)阻力比層流區(qū)阻力大的認(rèn)知是一致的了。以阻力開始突增的點(diǎn)為轉(zhuǎn)捩開始的位置，則m=1.5對應(yīng)的轉(zhuǎn)捩起始位置在x=240附近，而m=1對應(yīng)x=160附近。此外，轉(zhuǎn)捩之后與不加模型的平板阻力相比，m=1.5的總阻力仍然小一些，而m=1時反作用力帶來的阻力占據(jù)了主導(dǎo)，總阻力不降反而增加了不止一倍。

圖11 m=∞時的平均流向速度剖面

圖12所示為m=1對應(yīng)的不同位置的平均速度剖面。x=100、x=200處的速度分布仍與層流剖面相近，而在x=300、x=400處，速度剖面明顯被修正。與m=∞不同的是，修正的速度分布并不是像湍流剖面那樣逐漸變得飽滿。特別地，修正的速度剖面在壁面處的梯度比層流剖面明顯要小，這對應(yīng)了圖10中湍流區(qū)的壁面摩擦系數(shù)比層流還小的情形。壁面處緩慢的速度變化主要是由于m=1時模型帶來的流向阻力削弱了平均流向流動，而這種削弱作用在壁面處最為強(qiáng)烈。另一方面，這也使得平均流動的速度分布與典型的層流速度分布和湍流速度分布有本質(zhì)的差別。從圖12中可以看到，層流速度分布或者湍流速度分布是不存在拐點(diǎn)的，而當(dāng)m=1時，x=300、x=400處的速度分布顯然是存在拐點(diǎn)的。拐點(diǎn)的出現(xiàn)使得流動的穩(wěn)定性變?nèi)酰瑪_動更容易激發(fā)。這可能就是前面圖6～圖8m=1的結(jié)果中渦系結(jié)構(gòu)更加豐富、脈動強(qiáng)度顯著增強(qiáng)的原因所在。正因?yàn)橥牧鞯脑鰪?qiáng)，在大概z>1的高度范圍內(nèi)，x=300、x=400處平均速度分布的修正效果明顯。

圖12 m=1時的平均流向速度剖面

4 結(jié) 論

從觀察發(fā)現(xiàn)，羽毛表面流向和展向流動的阻力特征并不相同。本文以體積力的形式對這種各向異性阻力進(jìn)行了建模，并利用直接數(shù)值模擬的手段研究了其對平板邊界層轉(zhuǎn)捩的影響。研究得到如下主要結(jié)論：

1)準(zhǔn)流向渦是邊界層從層流轉(zhuǎn)捩到湍流過程中典型的大尺度渦結(jié)構(gòu)[34]，羽毛對展向流動的抑制作用可以削弱準(zhǔn)流向渦，從而延緩了轉(zhuǎn)捩的發(fā)生。即使邊界層已經(jīng)發(fā)展到了湍流狀態(tài)，因?yàn)榱飨驕u的削弱，導(dǎo)致湍流脈動強(qiáng)度顯著降低，壁面摩擦系數(shù)降低到與層流相當(dāng)?shù)乃健?/p>

2)若考慮羽毛側(cè)緣卷起對流向流動帶來額外阻力，相較于無流向阻力時，模型帶來的流向阻力使得轉(zhuǎn)捩開始得更早且湍流脈動更強(qiáng)。但在湍流區(qū)的摩擦阻力卻更小，甚至低于層流阻力。這是因?yàn)轭~外的流向阻力修正了來流方向的平均流動，使得壁面附近的速度變化得非常緩慢，速度梯度的減小了帶來了壁面摩擦阻力的減小。此外，修正的速度型出現(xiàn)了拐點(diǎn)，導(dǎo)致流動的穩(wěn)定性降低，從而促進(jìn)了湍流的發(fā)展。

3)盡管模型的流向阻力帶來了摩擦阻力的進(jìn)一步降低，但總阻力卻是增加的。流向阻力較大時，總的阻力甚至比典型湍流區(qū)阻力還大。因此，羽毛在展向和流向方向真實(shí)的力學(xué)特性關(guān)系對于減阻來說至關(guān)重要。

本文初步探究了羽毛對平板邊界層轉(zhuǎn)捩的影響，旨在為邊界層層流控制提供一種新的思路，希望可以拋磚引玉，與廣大研究者們共同深入探究。本文的模型及計(jì)算等方面仍存在很多不足之處，后續(xù)研究工作仍在進(jìn)行當(dāng)中。