流動隱身衣的均勻化設(shè)計與減阻特性*

王浩 姚能智 王斌 王學(xué)生

(華東理工大學(xué)機(jī)械與動力工程學(xué)院，上海 200237)

1 引言

2006 年，Pendry 等[1]證明了麥克斯韋方程組滿足空間坐標(biāo)變換的形式不變性，利用空間變換理論使光可以像流體一樣繞過物體傳播，實現(xiàn)了物體的隱身效果并提出光學(xué)隱身衣概念[2-4].隨后，眾多學(xué)者將變換理論拓展到了不同的領(lǐng)域，其中包括聲學(xué)[5-10]、熱學(xué)[11-24]、力學(xué)[25-28]等.流動隱身衣與光學(xué)隱身衣類似，流體在流經(jīng)隱身衣時也會主動繞道而行且不干擾流場.也正是其主動繞道而行的特性使得被隱身衣包裹的物體所承受的表面阻力顯著降低，這賦予了流動隱身衣在降低能耗、抵抗惡劣環(huán)境以及偽裝軍事目標(biāo)等方面的巨大潛力，十分具有研究意義.但是如何有效地調(diào)控流體繞過物體而不擾動流場目前還少有研究.

為此，Urzhumov 和Smith[29，30]利用達(dá)西定律和Brinkman-Stokes 方程所具有的空間坐標(biāo)變換形式不變性提出了流體動力學(xué)隱身衣(簡寫為流動隱身衣)的概念，并基于流動超材料設(shè)計理論構(gòu)造各向異性的滲透率實現(xiàn)了流體在多孔介質(zhì)中的隱身流動.然而，這種流動隱身衣因為需要負(fù)滲透率的多孔介質(zhì)被認(rèn)為難以實現(xiàn).隨后，Li 等[31]通過坐標(biāo)變換設(shè)計梯度深剖面結(jié)構(gòu)并利用實驗驗證了兩個可以用于聚集水波的環(huán)形裝置.Zou 等[32]受電磁波導(dǎo)隱身衣的啟發(fā)設(shè)計了一種梯度深剖面寬頻帶隱身衣，在實驗中實現(xiàn)了水波隱身.可是水波的傳播方式與流體的流動方式大相徑庭，這種水波隱身衣仍然無法有效調(diào)控流體流動.Tay 等[33]獨(dú)辟蹊徑，通過改變流域高度，不需要超材料即可實現(xiàn)流動隱身效果.Park 等[34]近期證明了流體動力學(xué)控制方程在一定約束條件下也同樣滿足空間坐標(biāo)變換的形式不變性，并基于此建立起變換流體動力學(xué).通過在流場中放置具有不同大小和方向的微米柱以獲得各向異性的動力粘度μ，成功調(diào)控了流體流動，完成了流動隱身衣實驗.基于變換流體動力學(xué)，Park 等[35]隨后又通過理論與實驗設(shè)計了具有不同功能的流動旋轉(zhuǎn)隱身衣、流動聚集隱身衣[36]等流動超材料.此外，Wang 等[37]與Boyko 等[38]利用流變學(xué)特性，通過施加溫度場、電場等方式均實現(xiàn)了流動隱身效果，但是這兩種方式都需要對流場持續(xù)輸入額外能量才能實現(xiàn).

基于變換理論設(shè)計隱身衣已經(jīng)有了令人欣喜的結(jié)果.然而在此之前，無論是在光學(xué)[1，2]、熱學(xué)[12，18]還是流體動力學(xué)[34，36]中，大多數(shù)基于坐標(biāo)變換設(shè)計的隱身衣都具有非均勻各向異性的設(shè)計參數(shù)，這使得我們在實際制備隱身衣時充滿了挑戰(zhàn).而本文基于變換流體動力學(xué)求解得到流動隱身衣的設(shè)計參數(shù)，利用等效介質(zhì)理論與積分中值定理將設(shè)計參數(shù)簡化，得到了流動隱身衣的均勻各向異性設(shè)計參數(shù).通過數(shù)值模擬證明了使用該均勻參數(shù)設(shè)計的均勻流動隱身衣，在不同流場中具有與非均勻流動隱身衣一樣的隱身效果.利用均勻參數(shù)設(shè)計制備流動隱身衣將大大降低實際制備難度，同時這一簡化方法也可以拓展到其他領(lǐng)域，如力學(xué)[20]、熱學(xué)[12，18]、光學(xué)[1，4]、電磁學(xué)[2]與聲學(xué)[5，6]等.此外，本文基于流動隱身衣在不同流場中的適用性設(shè)計了一種流動偽裝裝置，它可以將任意物體的流場偽裝成多個其他物體的流場，即可以將隱藏物體的流場分布偽裝成由其他物體所引起的期望流場分布.這種流動偽裝裝置實現(xiàn)了流場中任意物體的流場偽裝，為設(shè)計流動偽裝提供了新方法.最后，當(dāng)雷諾數(shù)不斷增加，計算流域脫離蠕動流進(jìn)入非蠕動流(層流范圍內(nèi))時，本文定量對比分析了均勻流動隱身衣的隱身效果.又考慮到流動隱身衣的實際設(shè)計意義，本文分析對比了均勻流動隱身衣所具有的減阻效率.結(jié)果表明，本文所設(shè)計的均勻流動隱身衣在非蠕動流中仍然具有較好的隱身性能與減阻性能.

2 均勻流動隱身衣的參數(shù)設(shè)計

在蠕動流情況下，對于不可壓縮的穩(wěn)態(tài)流動，流體運(yùn)動的連續(xù)性方程與動量守恒方程可表示為

其中u，μ，p分別表示速度矢量、流體動力粘度與壓強(qiáng).研究結(jié)果[18，34]表明，在蠕動流中方程(1)和方程(2)的確滿足空間坐標(biāo)變換的不變性.根據(jù)變換流體動力學(xué)，通過引入雅閣比空間變換矩陣(J)可以將方程(1)和方程(2)從虛擬空間x(x，y，z)變換到物理空間x′(x′，y′，z′)，于是有

其中u′=Ju/det(J)，

為變換動力粘度，用于調(diào)控隱身衣區(qū)域流體流動;J-1，J-T分別為雅閣比矩陣的逆矩陣與逆矩陣的轉(zhuǎn)置.

為實現(xiàn)流動隱身，需要對虛擬空間施加流動隱身所需的坐標(biāo)變換.因此流動隱身衣的空間坐標(biāo)變換映射關(guān)系式可以寫成如下形式:

其中R1，R2分別為流動隱身衣的內(nèi)外半徑，如圖1(a)所示.根據(jù)該映射關(guān)系，虛擬空間中心處一點(diǎn)被拉伸成物理空間中的一個半徑為R1的圓柱障礙物(r′ ＜R1)，而背景域(r′ ＞R2)的空間保持不變，作用效果如圖1(b) 所示.

圖1 模型示意圖 (a) 邊界條件;(b) 空間坐標(biāo)變換Fig.1.Schematic models:(a) Boundary conditions;(b) the coordinate transformation.

實現(xiàn)這種作用效果的關(guān)鍵之處，即虛擬空間中的流體流經(jīng)具有均勻各向異性粘度張量的流動隱身衣時，流體將被壓縮進(jìn)預(yù)定義的物理空間中.而這種變換效果則是由(5)式所示的變換粘度引起的，該變換粘度可以表示為[18，34]

值得注意的是，該變換粘度為不均勻各向異性張量，而根據(jù)流動超材料的等效介質(zhì)理論[34-36]，依靠微米柱的空間排布規(guī)則可使得流體在流經(jīng)每一層微米柱時都會有不同的流速，繼而等效為各向異性的粘 度張量，也即μeff=(〈u0〉s/〈upillar〉s)·μ，其 中upillar與u0分別代 表有無 微米柱 的速度 場，〈·〉s表示平均值，μ為背景流體粘度值.具體的制備方法可以在文獻(xiàn)[34]中找到.

不可忽略的是(7)式中的變換粘度在實驗室中往往很難設(shè)計，這是因為μ′是隨半徑r′不斷變化的.與熱學(xué)[18，24]、流體動力學(xué)[34-36]等領(lǐng)域相同，可以采用一種簡化方法得到簡化的變換粘度μ′′，將其表示為

實際上，(8)式在使用時為達(dá)到隱身效果還需要乘上補(bǔ)償系數(shù)，該補(bǔ)償系數(shù)可以通過優(yōu)化獲得.經(jīng)過簡化，變換粘度的切向分量便不再隨 半徑r′變化.但是，變換粘度的徑向分量卻仍取決于半徑r′，這意味著仍然需要非均勻的設(shè)計參數(shù).雖然這種經(jīng)過簡化的設(shè)計參數(shù)可以通過多層流體與微米柱交替的結(jié)構(gòu)實現(xiàn)在流動隱身衣中以階梯式的方式離散化[34])，但是這仍對制備工藝要求較高.

為了進(jìn)一步降低制備難度，可以采用等效介質(zhì)理論與積分中值定理處理簡化過的變換粘度.經(jīng)過計算，可以得到均勻的流動隱身衣，其平均變換粘度為

其中Rm=(R1+R2)/2，ΔR=R2-R1，經(jīng)過積分可以得到

求解定積分還需要代入流動隱身衣的內(nèi)徑R1與外徑R2分別作為積分的下限與上限.然而，r′=R1為(10)式的一個奇點(diǎn)，這導(dǎo)致無法將R1作為下限代入(10)式求解徑向平均粘度.實際上，在自然界中也并不存在粘度為無窮大的流體.不過由于在r′=R1處為無 滑移固 體壁面(固體壁 面可等效為粘度無窮大的流體)，流體無法穿透該壁面，所以此處的等效粘度已經(jīng)天然滿足流動隱身衣對徑向粘度值無窮大的要求.于是在求解平均變換粘度時，便可以取壁面R1附近的半徑作為積分下限，即NR1(R2/R1＞N ＞1)，代入求得定積分結(jié)果為

圖2 (a)徑向粘度與(b)切向粘度隨 r ′/R1 的分布Fig.2.Distributions of (a) radial viscosity and (b) azimuthal viscosity with r ′/R1 .

3 數(shù)值模擬結(jié)果與討論

根據(jù)得到的變換粘度值等式(7)式，(12)式和(13)式，可以使用有限元軟件COMSOL Multiphysics 中提供的自定義偏微分方程(Partial Differential Equation，PDE)計算模塊驗證我們的設(shè)計結(jié)果.其中，COMSOL 的離散化方案與線性算法采用Galerkin 有限元方法.

如圖1(a)所示，流動隱身衣的內(nèi)外徑分別為R1=1.5 cm，R2=3 cm (基于此計算得到N≈1.1).計算流域尺寸為L×H ×D，即是流域長10 cm×高10 cm×深50 μm，其中深度為流域沿z軸的尺寸.由于流域的深度遠(yuǎn)遠(yuǎn)小于流域的長度與寬度，即D?L和H，該流場的流動為Hele-Shaw 流動[39]，可以被視為二維流動.流場沿y方向自下向上流動，且進(jìn)出口壓差 Δp=1 kPa (p1＞p2) ;流場的左右兩邊界均為無滑移壁面.此外，計算流域的流體選用室溫293.15 K 的水，該流體介質(zhì)的動力粘度μ=1×10-3Pa · s，密度ρ=9 97.1kg/m3.

圖3 所示為均勻背景流場下的速度分布圖.圖3(a)中給出了背景流場，流場中流線與等壓線都平行且均勻分布.當(dāng)圓柱形障礙物放入流場后，如圖3(b)所示，流場中原本平行均勻分布的流線與等壓線由于障礙物的干擾變得扭曲.在圖3(b)的基礎(chǔ)上將非均勻流動隱身衣添加在障礙物周圍，如圖3(c)所示，在障礙物“穿”上隱身衣之后，位于隱身衣之外(r′ ＞R2)的流線與等壓線又一次恢復(fù)到圖3(a) 的狀態(tài)，達(dá)到了中心障礙物隱身的效果，這與Park 等[34]設(shè)計的非均勻流動隱身衣的作用效果一致.同樣地，當(dāng)使用本文設(shè)計的均勻流動隱身衣替換障礙物周圍的非均勻流動隱身衣時，如圖3(d)所示，隱身衣外的流場也恢復(fù)到圖3(a)的狀態(tài)，這表示均勻流動隱身衣具有與非均勻流動隱身衣相同的隱身效果.

圖3 均勻流場下的速度分布圖，黑色線條為流線，白色線條為等壓線 (a) 不含障礙物;(b) 含障礙物;(c) 含障礙物與非均勻隱身衣;(d) 含障礙物與均勻隱身衣Fig.3.Velocity distributions superimposed with streamlines (black color) and isobars (white color) for uniform flow fields:(a) Obstacle absent;(b) obstacle existent;(c) obstacle and inhomogeneous cloak existent;(d) obstacle and homogeneous cloak existent.

為了進(jìn)一步定量化對比驗證所設(shè)計的均勻流動隱身衣，可以在r′=R2處取若干點(diǎn)，在背景流場、含非均勻流動隱身衣的流場與含均勻流動隱身衣的流場中分別對這些點(diǎn)的壓強(qiáng)進(jìn)行比較.這些點(diǎn)在隱身衣作用前后壓強(qiáng)偏差值的大小，即可用來表示流動隱身衣作用后對背景流場的干擾程度.如圖4(a)所示，在計算流域r′=R2處取24 個點(diǎn)進(jìn)行對比，可以看到這些點(diǎn)在非均勻流動隱身衣與均勻流動隱身衣作用下的壓強(qiáng)均與背景流場中該點(diǎn)處的壓強(qiáng)重合，這表示我們設(shè)計的均勻流動隱身衣與非均勻流動隱身衣一樣具有不干擾背景流場的效果.此外，如圖4(b)所示，可以取y=0 處流場的速度分布圖與背景流場的流速vin進(jìn)行比較，通過比較確定流動隱身衣外 (|x/R2＞1|) 的速度與背景流場的速度值相同，也可證明我們所設(shè)計的均勻流動隱身衣具有流動隱身效果.

圖4 (a) 壓強(qiáng)分布圖，黑色線條、紅色三角形與藍(lán)色三角形分別表示圖3(a)、圖3(c)與圖3(d)在 r ′=R2 處的壓強(qiáng);(b) 速度分布圖，黑色線條、紅色線條與藍(lán)色線條分別表示圖3(a)、圖3(c)與圖3(d)在 y=0 處的速度分布Fig.4.(a) Pressure distributions.Black line，red triangle and blue triangle indicate the pressure at r ′=R2 in Fig.3(a) (only background)，F(xiàn)ig.3(c) (with inhomogeneous cloaks) and Fig.3(d) (with homogeneous cloaks) respectively.(b) Velocity distributions.Black line，red line and blue line indicate the velocity at y=0 in Fig.3(a) (only background)，F(xiàn)ig.3(c) (with inhomogeneous cloaks) and Fig.3(d) (with homogeneous cloaks) respectively.

上述數(shù)值模擬驗證都是基于均勻背景流場進(jìn)行的，為了進(jìn)一步驗證所設(shè)計的均勻流動隱身衣在非均勻流場中的可行性，將計算流域的進(jìn)出口邊界條件進(jìn)行修改.分別將進(jìn)出口邊界的左半邊與右半邊作為新的進(jìn)出口邊界條件并將其余邊界均設(shè)置為無滑移壁面，進(jìn)出口壓力與其他參數(shù)保持不變，得到的數(shù)值模擬流場速度分布圖如圖5 所示.

如圖5(a)所示為非均勻背景流場，流場內(nèi)的流線與等壓線不再均勻分布.將圓柱形障礙物與非均勻流動隱身衣“放入”到該非均勻流場中，如圖5(b)所示，同均勻流場一樣，流動隱身衣外(r′ ＞R2)的流線、等壓線再一次恢復(fù)到圖5(a)背景流場中的分布狀態(tài).如果采用所設(shè)計的均勻流動隱身衣替換圖5(b)中的非均勻流動隱身衣，如圖5(c)所示，背景流場中的流線、等壓線分布狀態(tài)與圖5(a)所示相同，這表示本文設(shè)計的均勻流動隱身衣一樣可以實現(xiàn)非均勻流場中障礙物的流動隱身效果.

圖5 非均勻流場下的速度分布圖，黑色線條為流線，白色線條為等壓線 (a) 不含障礙物;(b) 含障礙物與非均勻隱身衣;(c) 含障礙物與均勻隱身衣Fig.5.Velocity distributions superimposed with streamlines (black color) and isobars (white color) for non-uniform flow fields:(a) Obstacle absent;(b) obstacle and inhomogeneous cloak existent;(c) obstacle and homogeneous cloak existent.

實際上，流動隱身衣除了可以在流場中隱藏目標(biāo)以外，還可以實現(xiàn)流場的偽裝分布.如圖A1 所示，當(dāng)對指定障礙物進(jìn)行流動隱身的同時在隱身衣周圍放置其他偽裝物體，那么隱身衣外流場的表現(xiàn)僅與偽裝物體有關(guān)，于是可以將指定障礙物隱蔽并將原始流場分布偽裝為其他物體的流場分布.

根據(jù)上述分析驗證，可以確定所設(shè)計的均勻流動隱身衣在蠕動流中是有效的.雖然在設(shè)計之初并沒有考慮對流項(ρu·?u)，但是為了進(jìn)一步確認(rèn)所設(shè)計的均勻流動隱身衣在雷諾數(shù)更大的非蠕動流中是否有效，將在控制方程(2)和方程(4)中加入對流項后做進(jìn)一步分析驗證.而對于定量評定隱身衣的隱身效果，仍然采用與圖4(a)類似的辦法.在背景流場、含非均勻流動隱身衣的流場與含均勻流動隱身衣的流場中，分別在隱身衣外徑R2處取n個點(diǎn)，計算其均方差Ω，可以表示為

其中pc，pb分別表示含隱身衣流場與不含隱身衣流場中R2處的壓強(qiáng)，計算得出的均方差Ω就可以表示放入隱身衣后的流場與純背景流場的差異程度.如圖6 所示，為不同雷諾數(shù)下計算得到的均方差，隨著雷諾數(shù)逐漸增加進(jìn)入非蠕動流階段后，無論是均勻隱身衣還是非均勻隱身衣，計算得到的Ω都在持續(xù)增加，這表示兩種隱身衣的隱身效果都在逐漸降低.這主要是因為此時控制方程中的對流項在起作用，通過方程(1)—(5)(忽略對流項)設(shè)計出的流動隱身衣已經(jīng)不能完美適用.此外，通過數(shù)值模擬得到當(dāng)雷諾數(shù)Re=1 與Re=10 時兩種流動隱身衣的速度分布圖，如圖6(c)和圖6(d)所示.雖然此時的Ω相比于蠕動流已經(jīng)增加很多，但流動隱身衣也仍然具有一定適用性.

圖6 (a) 隱身效果，其中Ω 表示含隱身衣流場與純背景流場的差異程度;(b) 減阻效果，其中 ηreduce 表示隱身衣的減阻程度;(c)當(dāng)流域雷諾數(shù) R e=1 時，非均勻流動隱身衣(c1)與均勻流動隱身衣(c2)的速度分布;(d) 當(dāng)流域雷諾數(shù) R e=10 時，非均勻流動隱身衣(d1)與均勻流動隱身衣(d2)的速度分布Fig.6.(a) Cloaking effects，in which Ω presents the difference between flow fields with cloaks and pure background flow fields;(b) Drag reduction effects，in which η reduce denotes the degree of drag reduction of cloaks;(c) Velocity distributions of inhomogeneous cloaks (c1) and homogeneous cloaks (c2) at R e=1 ;(d) Velocity distributions of inhomogeneous cloaks (d1) and homogeneous cloaks (d2) at R e=10 .

盡管隨著雷諾數(shù)的增加，均勻流動隱身衣的隱身性能要比非均勻隱身衣差，如圖6(a)所示.但是，流動隱身衣具有的現(xiàn)實意義更偏向于降低隱身對象所受到的流動阻力.因此，為定量分析均勻流動隱身衣的減阻性能，采用阻力計算公式[40]計算隱身障礙物的表面阻力.如圖6(b)所示，隨著雷諾數(shù)的不斷增加，雖然已經(jīng)進(jìn)入非蠕動流階段，但是兩種流動隱身衣的隱身效率仍然較高，并且均勻流動隱身衣的減阻效率要更好.根據(jù)圖6，當(dāng)雷諾數(shù)增加時，本文所設(shè)計的均勻流動隱身衣相比于非均勻流動隱身衣雖然在隱身效果上有所降低，卻擁有了更佳的減阻能力.

4 結(jié)論

從變換流體動力學(xué)出發(fā)，可以得到傳統(tǒng)流動隱身衣的設(shè)計參數(shù)，而該設(shè)計參數(shù)的非均勻性對實際制備流動隱身衣造成了較大的阻礙.基于此，本文提出了一種將非均勻流動超材料簡化為均勻流動超材料的方法，該簡化方法利用等效介質(zhì)理論與積分中值定理，可以得到一組均勻各向異性流動隱身衣的設(shè)計參數(shù).通過有限元仿真軟件的數(shù)值模擬，驗證了基于均勻設(shè)計參數(shù)設(shè)計的流動隱身衣與非均勻流動隱身衣具有相同的隱身效果，而且同樣適用于多種流場.這種簡化非均勻超材料的方法可以拓展到其他領(lǐng)域，采用該方法得到的均勻設(shè)計參數(shù)將大大降低實際制備超材料的難度，進(jìn)一步推動超材料的實際應(yīng)用.此外，基于所設(shè)計均勻流動隱身衣對不同流場的適用性，提出了一種可以改變流場分布的流動偽裝裝置.最后，通過對比分析發(fā)現(xiàn)，均勻流動隱身衣的隱身效果隨著雷諾數(shù)增加而不斷降低，盡管不能完美適用于層流流動，但對于緩慢層流流動狀態(tài)也仍具有一定適用性.通過分析被隱身對象的表面阻力特性發(fā)現(xiàn)，均勻流動隱身衣的減阻效率相比于非均勻流動隱身衣要更高.隨著雷諾數(shù)的不斷增加，均勻或非均勻流動隱身衣的減阻性能均有所降低，但仍然都擁有較高的減阻效率，且均勻流動隱身衣的減阻性能更穩(wěn)定.

附錄A 流動偽裝裝置

圖A1(a)所示為原始流場的分布圖，流域中心為圓柱形障礙物，此時根據(jù)流場分布狀態(tài)可以輕易地判斷出障礙物的位置信息.而當(dāng)中心障礙物“套上”流動隱身衣后，如圖A1(b)所示，此時再根據(jù)隱身衣外流場的分布狀態(tài)就不能確定障礙物的位置.在此基礎(chǔ)上將流動隱身衣左右分別添加一個圓柱形物體作為偽裝物體，如圖A1(c))和(圖A1(d)所示，此時流動偽裝裝置外(紅色虛線外的區(qū)域)流場的分布狀態(tài)僅與左右兩個偽裝物體的位置和大小有關(guān)，再根據(jù)流場分布狀態(tài)很難辨別出中心障礙物的位置信息.這樣一種具有偽裝功能的裝置，可以將原始由中心障礙物所產(chǎn)生的流場分布變成偽裝物體所產(chǎn)生的流場分布.

圖A1 流動偽裝示意圖 (a) 僅含障礙物;(b) 隱身狀態(tài);(c) 一個障礙物被偽裝成兩個;(d) 含兩個障礙物.從偽裝圖中可以看出，圖(c)紅色虛線外的流場與圖(d)兩個目標(biāo)物體的流場一致，說明一個物體已經(jīng)被隱藏并偽裝成兩個物體Fig.A1.Schematic of hydrodynamic camouflage:(a) Only one obstacle exists;(b) cloaking;(c) one obstacle camouflaged as two obstacles;(d) two obstacles exist.As can be seen from the camouflage devices，the flow fields outside the red dashed line in panel(c) coincide with the flow fields of the two-target objects in panel (d)，indicating that one object has been hidden and camouflaged as two objects.