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

王浩 姚能智 王斌 王學生

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

1 引言

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

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

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

2 均勻流動隱身衣的參數設計

在蠕動流情況下，對于不可壓縮的穩態流動，流體運動的連續性方程與動量守恒方程可表示為

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

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

為變換動力粘度，用于調控隱身衣區域流體流動;J-1，J-T分別為雅閣比矩陣的逆矩陣與逆矩陣的轉置.

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

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

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

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

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

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

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

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

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

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

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

3 數值模擬結果與討論

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

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

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

圖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.

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

圖4 (a) 壓強分布圖，黑色線條、紅色三角形與藍色三角形分別表示圖3(a)、圖3(c)與圖3(d)在 r ′=R2 處的壓強;(b) 速度分布圖，黑色線條、紅色線條與藍色線條分別表示圖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)，Fig.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)，Fig.3(c) (with inhomogeneous cloaks) and Fig.3(d) (with homogeneous cloaks) respectively.

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

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

圖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.

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

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

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

圖6 (a) 隱身效果，其中Ω 表示含隱身衣流場與純背景流場的差異程度;(b) 減阻效果，其中 ηreduce 表示隱身衣的減阻程度;(c)當流域雷諾數 R e=1 時，非均勻流動隱身衣(c1)與均勻流動隱身衣(c2)的速度分布;(d) 當流域雷諾數 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 .

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

4 結論

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

附錄A 流動偽裝裝置

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

圖A1 流動偽裝示意圖 (a) 僅含障礙物;(b) 隱身狀態;(c) 一個障礙物被偽裝成兩個;(d) 含兩個障礙物.從偽裝圖中可以看出，圖(c)紅色虛線外的流場與圖(d)兩個目標物體的流場一致，說明一個物體已經被隱藏并偽裝成兩個物體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.