柔性仿羽毛結構抑制邊界層轉捩的初步探索

葉正寅, 洪 正, 武 潔

(西北工業大學, 西安 710072)

0 引 言

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

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

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

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

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

通過對邊界層的大量研究[24-29]，人們發現存在著一種共性特點，即：邊界層從層流轉化為湍流的轉捩過程，必然出現Λ渦和發卡渦，而這兩種渦的“腿”就是一種傾斜的局部流向渦，它會誘導氣流的展向流動，如果羽毛規則的柔性邊緣因為展向流動而卷起，就會對局部展向流動有阻滯作用并抑制Λ渦和發卡渦的發展，從而抑制轉捩的發生。

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

1 羽毛模型

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

1.1 展向阻力模型

對于卷起的羽毛側緣給當地展向流動帶來的阻力，有以下幾點考慮：

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

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

(3)這種阻力只能使得當地展向脈動速度減小，但不能改變其方向。

羽毛阻礙展向脈動流動產生的阻力，概念上與黏性阻力類似，則基于(1)、(2)兩點且類比黏性力的公式，可將沿展向的阻力模化為：

(1)

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

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

(2)

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

(3)

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

(4)

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

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

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

1.2 流向阻力模型

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

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

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

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

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

(5)

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

(6)

2 數值方法和計算設置

2.1 控制方程

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

(7)

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

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

2.2 計算設置

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

圖2 計算域示意圖

(8)

入口處擾動的參數選取參考文獻[31]中的設置，如表1所示。

表1 入口擾動參數

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

(9)

表2 計算域尺寸和網格參數

3 結果分析與討論

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

3.1 平板邊界層轉捩

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

(10)

(11)

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

圖3 層流區和湍流區的速度剖面

壁面摩擦系數Cf的定義為：

(12)

對于層流，摩擦系數為

(13)

式中：δ是當地的邊界層位移厚度。Ducros等人[33]提出了一種經驗公式來描述湍流邊界層摩擦系數的空間演化規律，具體形式為：

(14)

式中：l′=125δ0。圖4中給出了平板上轉捩流動，層流以及Ducros經驗公式的壁面摩擦系數。從圖中可以看到，大約從x=80開始，平板邊界層的摩擦系數開始快速偏離層流對應的值，該點可視為層流向湍流轉捩的起始點。摩擦系數在x=190達到峰值，之后流動進入完全湍流狀態。湍流區的摩擦系數與經驗公式的結果所在的區間是一致的。

圖4 沿壁面的摩擦系數分布

平板邊界層從層流到湍流的演化過程中，存在一些典型的大尺度相干結構。利用Q準則，圖5呈現了轉捩過程中不同階段的渦結構，用距離壁面的高度進行了著色。為了更好地顯示效果，展向通過復制延伸了一倍。可以看到轉捩開始后先是形成Λ渦，然后是Λ渦演化成發卡渦，進而生成環狀渦，并最終發展成環狀渦鏈結構。對典型渦系結構的清晰捕捉也體現了數值模擬的可靠性。

(a)t=1.6T0

3.2 羽毛模型對平板轉捩的影響

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

圖6給出了加入模型后計算得到的全流場的渦系結構圖，上一小節中未加模型的結果也一并給出作為對比。可以看到在加入羽毛的阻力模型后，開始形成大尺度渦的位置較無模型情況明顯靠后。無模型情況下，x=100附近就已經形成了發卡渦，而加入模型且m=∞時，x=300處發卡渦才出現，并且轉捩到湍流過程中形成的渦系要稀疏得多。不同的m值得到的結果有顯著的差異。m=1.5與m=∞相比，開始形成典型渦結構的位置相差不大，但后續過程中形成的渦系結構要豐富一些。m=1情況下，開始形成渦的位置較m=∞明顯提前，x=200附近已經可以看到形成的發卡渦了。此外，之后形成的渦系結構相較之下也要豐富得多。

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

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

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

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

圖9 不同m值得到的摩擦系數分布

(15)

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

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

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

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

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

4 結 論

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

1)準流向渦是邊界層從層流轉捩到湍流過程中典型的大尺度渦結構[34]，羽毛對展向流動的抑制作用可以削弱準流向渦，從而延緩了轉捩的發生。即使邊界層已經發展到了湍流狀態，因為流向渦的削弱，導致湍流脈動強度顯著降低，壁面摩擦系數降低到與層流相當的水平。

2)若考慮羽毛側緣卷起對流向流動帶來額外阻力，相較于無流向阻力時，模型帶來的流向阻力使得轉捩開始得更早且湍流脈動更強。但在湍流區的摩擦阻力卻更小，甚至低于層流阻力。這是因為額外的流向阻力修正了來流方向的平均流動，使得壁面附近的速度變化得非常緩慢，速度梯度的減小了帶來了壁面摩擦阻力的減小。此外，修正的速度型出現了拐點，導致流動的穩定性降低，從而促進了湍流的發展。

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

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