黃苗苗，張 楠，卜淑霞，朱愛軍，張 華

（中國船舶科學(xué)研究中心水動力學(xué)重點實驗室，江蘇無錫214082）

0 引 言

無論對于水面艦船還是水下航行體，快速性一直都是綜合航行性能的關(guān)鍵組成部分之一。經(jīng)過長期的研究和發(fā)展，通過線型優(yōu)化降低阻力的空間和幅度越來越小，進一步的減阻必須借助于新型減阻技術(shù)。海洋魚類體表結(jié)構(gòu)經(jīng)過億萬年的自然優(yōu)化選擇，進化形成了許多特異結(jié)構(gòu)和優(yōu)異減阻功能，具有非凡的低阻高效游動能力。仿生減阻研究的重要啟示首先來自于對海豚的研究[1-2]。1936年，英國生物學(xué)家James Gray[1]基于海豚快速游動現(xiàn)象提出了著名的格雷悖論，從而引發(fā)了之后幾十年對海洋高速魚類仿生減阻研究的熱潮[3-7]。

國外在上世紀六十年代就開展了仿生柔性面的減阻研究，包括實驗及數(shù)值計算方法研究[8-11]。我國從上世紀末開始柔性面的理論和試驗研究，如中國船舶科學(xué)研究中心的張效慈[12]開展了柔性面減阻降噪的理論分析；黃微波等[13]采用CFD 模擬研究了定常隨行波作用下的渦街形式。整體來看，與國外相比我國在仿生柔性面減阻方面的研究較少，特別是數(shù)值計算方面有待深入研究。

Carpenter 等[14]在研究海豚皮膚時發(fā)現(xiàn)：海豚表皮光滑，在皮膚表面下隱藏著一些纖維結(jié)節(jié)，當(dāng)海豚游動時，隨著滑過海豚體表的水流剪切力的增大，海豚皮膚逐漸由光滑轉(zhuǎn)變成具有一定幾何形狀的非光滑形態(tài)，從而實現(xiàn)減阻目的。也就是說，活體海豚皮膚的柔性是自主依靠肌肉上下垂直抖動形成表面有規(guī)律的動態(tài)起伏。而這方面的研究比較欠缺，因此，本文開展了鯨豚表皮微尺度柔性變形的非定常數(shù)值計算研究。

1 生物觀測

在柔性面計算參數(shù)的選取時，本文參考了國外文獻中的海豚表皮觀測結(jié)果。雖然用肉眼觀察比較模糊，但是在很多齒鯨（如海豚）的皮膚表面確實分布著波紋狀的規(guī)則溝嵴。這些溝嵴較淺，在海豚死后會因為皮膚松弛而模糊或消失不見，同時很難在皮膚的組織切片中觀察到，文獻中采用快速固化方法對活體齒鯨亞目動物的皮膚嵴進行了壓印和測量[15]。

圖1 海豚背部皮膚嵴[15]Fig.1 Surface cutaneous ridges on the dorsal thorax of a dolphin[15]

由于這種皮膚嵴的尺寸很小，都是微米（μm）量級，觀測活體鯨豚尚需使其處于安靜狀態(tài)，并借助于特別的處理方法。對于研究仿生減阻的人，我們關(guān)心的是這種皮膚嵴在鯨豚高速游動過程中是如何變化發(fā)揮作用的，抑或是不動即可減阻？本來通過實時觀測就能直接回答這個問題，但是目前的測試水平尚達不到，因此本文開展了這方面的數(shù)值研究工作來進行探索。

2 計算方法

文中采用計算效率高，預(yù)報精度也能滿足工程實用的RANS 方法對柔性表皮開展阻力流場的數(shù)值模擬，數(shù)學(xué)控制方程包括連續(xù)性方程和Reynolds方程[16]，湍流模型采用k-ε模型。

文中數(shù)值計算涉及的柔性變形依靠動網(wǎng)格技術(shù)實現(xiàn)。對于動網(wǎng)格，在任意一控制體?V（邊界是運動的）上，標(biāo)量?積分形式的守恒方程為

守恒方程中的時間導(dǎo)數(shù)項采用一階后向差分格式寫作

圖2 是本文計算模型，頂部為剛性無變形的平板固壁，底部為柔性面，兩者為平行狀態(tài)。左側(cè)設(shè)為水流入口，右側(cè)設(shè)為出口，計算中出入口均設(shè)為周期性邊界條件。計算域底部邊長記為L，計算域大小為L×L×5L，頂、底部距離Δz=5L，以防止頂部固壁邊界層影響下層流場。本文計算模型采用結(jié)構(gòu)化網(wǎng)格，網(wǎng)格總數(shù)為400萬。

圖2 計算域及柔性面上的網(wǎng)格分布Fig.2 Computational domain and mesh pattern for compliant surface

本文選用剛性波來代替海豚可能不變形的波狀規(guī)則溝嵴，用駐波和流向隨行波來代替海豚皮膚可能的柔性變形，尚不考慮其表面的彈性影響。

剛性波面表達式為

駐波柔性面表達式為

流向隨行波柔性面表達式為

3 模擬結(jié)果

為了探尋可能的仿生減阻方案，本文對剛性的及柔性的波狀壁面進行了阻力數(shù)值模擬，作為對比，同步計算了剛性平板。計算工況Re范圍為106～107。

3.1 計算精度檢驗

3.1.1 與實驗結(jié)果比較

本文計算方法與試驗結(jié)果給出的平板壁面無量綱速度分布曲線如圖3 所示，圖中圓點為Abra?ham的經(jīng)典實驗結(jié)果[17]，實線為本文數(shù)值計算結(jié)果，兩者總體吻合良好。

3.1.2 與經(jīng)典擬合公式比較

采用Schlichting 平板湍流邊界層壁面局部摩擦系數(shù)近似公式（7）[18-19]，計算本文剛性平板局部摩擦系數(shù)，并與數(shù)值模擬結(jié)果對比。如表1所示，兩者差異基本在1%之內(nèi)，說明本文的數(shù)值計算方法是可行的。

圖3 平板壁面無量綱速度試驗對比Fig.3 Comparison of dimensionless velocity on plate between the simulation and experiment results

表1 平板摩擦阻力系數(shù)對比結(jié)果Tab.1 Comparison of plate friction coefficients

3.2 流場分析

針對不同壁面表達式（4）～（6），選取入口速度5 m/s時的流場開展對比分析。下文中涉及的流向剖面為本文計算域中y=0.5L處的xz平面，如圖4所示。

3.2.1 剛性波面ζ = acoskx

圖5給出了對應(yīng)的速度云圖。右向流動的邊界層內(nèi)，由波峰→波谷→波峰形成的壁面凹陷，產(chǎn)生了流道周期性收縮擴散區(qū)，對應(yīng)了近壁面流場的高速與低速區(qū)。故速度等值線呈波狀分布，且與波面基本上呈反相。其中波谷→波峰面為迎流物面，它的壓力顯然應(yīng)該高于相鄰波峰→波谷的背流面。于是產(chǎn)生流向的壓差阻力，這是平板繞流不存在的。

圖4 流向剖面示意圖Fig.4 Streamwise section

圖5 剛性波面流向剖面速度云圖Fig.5 Velocity magnitude of rigid wave surface

圖6 則是這一凹陷區(qū)的速度矢量圖。對于右向流動，它清楚顯示了順時針旋轉(zhuǎn)的駐渦存在。這一駐渦將使凹陷區(qū)相當(dāng)部分壁面速度梯度比平板小，甚至梯度變負。這樣，即使在波峰處，速度梯度高于平板，但其總粘性阻力將比平板低。由于順時針駐渦的存在，流向最大流速不在剛性波面的波谷處，而是在波谷稍后的下游處，所以速度波狀等值線的“谷”仍然對應(yīng)壁面波峰處，而“峰”則對應(yīng)壁面波谷稍后的下游處。

圖6 剛性波面近壁區(qū)速度矢量圖Fig.6 Velocity vector of rigid wave surface

圖7 駐波柔性面流向剖面速度場云圖Fig.7 Velocity magnitude of active compliant surface with standing waveform

3.2.2 駐波波面ζ = acoskxcosωt

圖7 為穩(wěn)定計算結(jié)果中選取的一個周期內(nèi)速度云圖。自t=0 時刻（圖7（a））起：節(jié)點A、B 之間波面從波峰降至平面（t=T/4），再降至波谷（T/2），一直在讓出流道空間，對應(yīng)z 向便存在較大低速影響區(qū)；而此時相鄰的節(jié)點BC 波面自波谷至平面（T/4），再至波峰（T/2），不斷壓縮流道空間，從而z向為較小的高速影響區(qū)。故其速度等值波狀線與運動的駐波面反相。

雖然波峰、波谷不斷在交替變化，但造成壁面高低不平必然會產(chǎn)生壓差阻力。若駐波最大波高與剛性波一致，在t＞0 時駐波壁面的高低不平程度低于剛性波壁面，可以預(yù)期對應(yīng)的壓差阻力將比剛性波小。

圖8 給出了駐波面的峰→谷→峰凹陷區(qū)，同樣存在順時針旋轉(zhuǎn)的駐渦，這意味著其摩擦阻力應(yīng)低于平板。

3.2.3 流向隨行波ζ = acos( kx - ωt)

對于隨行波壁面上任意兩個質(zhì)點，它們對應(yīng)的z向無窮遠處均為來流u，而它們自身僅在z向上下運動，規(guī)律完全一樣，只差相位。這就注定了邊界層流場的速度等值線，應(yīng)該有隨行波面類似形狀，并且是同相的。圖9的速度云圖證實了這種分析。

在大地坐標(biāo)系中觀察波峰→波谷→波峰形成的凹陷區(qū)。對于剛性波面，位置、大小都是不變的；對于駐波面，位置不變，波高呈周期性正負變化。雖然隨行波以波形的波速c=流向移動，但波面質(zhì)點x位置固定不變。只是位于“波峰→波谷”波面上的質(zhì)點總是以有規(guī)律的正z向速度驅(qū)趕流體，使“峰→谷”波面以流向速度c擠壓凹陷區(qū)中的流體空間，而凹陷區(qū)中的“谷→峰”波面上質(zhì)點以一定規(guī)律沿負z向速度運動，使谷→峰波面以流向速度c讓出空間，接納上游被驅(qū)趕來的流體，如圖10速度矢量所示。

圖8 駐波柔性面近壁區(qū)速度矢量圖Fig.8 Velocity vector of active compliant surface with standing waveform

圖9 隨行波柔性面流向剖面速度場云圖Fig.9 Velocity magnitude of active compliant surface with traveling waveform

這一物理過程表明，隨著流向波速由零（剛性波面）開始增加，“谷→峰”波面將由第3.2.1 節(jié)中迎流面角色轉(zhuǎn)變?yōu)椴粩嘧尦隽鲌隹臻g的角色，壓力分布由高變低；而“峰→谷”波面由第3.2.1 節(jié)中背流面角色轉(zhuǎn)變?yōu)椴粩囹?qū)趕、擠壓流體向下游的角色，壓力分布由低變高。從而使壓差阻力由剛性波面的值不斷減小，乃至變?yōu)橥屏Α?/p>

這一物理過程同時表明，流向波速由零增加，意味著凹陷區(qū)中近壁面流向速度增加，并逐漸減弱原剛性波面凹陷區(qū)中的駐渦強度。只要波速足夠大，駐渦強度會低于駐波、剛性波，乃至消失，并不斷使壁面正向速度呈梯度增加。所以可預(yù)計到摩擦阻力將不斷增加。

圖10 隨行波柔性面近壁區(qū)速度矢量圖Fig.10 Velocity vector of active compliant surface with traveling waveform

3.3 阻力的數(shù)值預(yù)報

對文中四種面板開展阻力模擬對比，涉及的無量綱參數(shù)如下：

式中，δv為壁面粘性單位，λ為柔性面波動的波長，a為柔性面波動的波幅，T為波動周期，u為流速，c+為隨行波無量綱行波波速。

式中，F(xiàn)t為總阻力，F(xiàn)p為其中壓力積分值，F(xiàn)τ為其中的剪切應(yīng)力積分值，s為面板表面積。結(jié)果處理中柔性面與平板均取實際面積，駐波形式的柔性面面積做周期變化，因此取最大與最小波幅的面積平均值。受力分析均采用計算穩(wěn)定后的結(jié)果，其中流向隨行波及駐波柔性面均選取計算穩(wěn)定后兩個周期內(nèi)的時歷監(jiān)測值作為最終計算取值。

為了便于比較，表2～4中所有算例的波陡全部一樣。

3.3.1 剛性波面ζ = acoskx

從表2可以看出，剛性波面摩擦阻力系數(shù)Cτ如同平板一樣，隨來流速度增加而下降。由于第3.2.1節(jié)所述順時針駐波的存在，相同來流速度下，剛性波面摩擦阻力系數(shù)Cτ均比平板小。

邊界層內(nèi)來流經(jīng)剛性波面，所產(chǎn)生的壓差阻力系數(shù)Cp隨速度增加而增加。但剛性波面總的阻力系數(shù)Ct=Cp+Cτ仍然在所有算例中都高于平板，所以剛性波面相對于平板是增阻的。

3.3.2 駐波波面ζ = acoskxcosωt

表3同樣顯示，隨著來流速度提高，駐波波面的摩擦阻力系數(shù)與平板一樣，是下降的，由于第3.2.2節(jié)所述駐渦的存在，Cτ均比平板小。

邊界層來流經(jīng)駐波波面后，同樣產(chǎn)生了壓差阻力系數(shù)Cp。而駐波波面的總阻力系數(shù)Ct仍然均高于對應(yīng)的平板。所以，駐波波面是增阻的。

表3 平板與駐波變形的柔性面阻力對比計算Tab.3 Resistances of rigid plate and active compliant surface with waveform acoskxcosωt

3.3.3 流向隨行波波面ζ = acos( kx - ωt)

表4 給出的是隨行波面與平板阻力系數(shù)的比較。根據(jù)第3.2.3 節(jié)流動物理過程及機理分析，將重要參數(shù)波速c+=與其他參數(shù)一起列于表4中。

表4 平板與流向隨行波變形的柔性面阻力對比計算Tab.4 Resistances of rigid plate and active compliant surface with waveform acos( kx - ωt)

（1）Case 12與Case 13、Case 19與Case 20是兩組除波速外，其他條件完全一樣的隨行波。計算值比較表明，無論來流是低速還是高速，波速提高，Cτ增加，Cp下降。

（2）Case 12到Case 20算例，結(jié)合表3剛性波面情況，可說明：當(dāng)c+從零開始增加，隨行波摩擦阻力系數(shù)將從低于平板值增加，約在c+≈0.93 時與平板持平，隨后超過平板；而隨行波的壓差阻力系數(shù)將從剛性波Cp＞0開始逐步下降，約在c+≈0.73時開始變?yōu)镃p＜0，并繼續(xù)減小。

（3）隨著c+由零增加，隨行波總阻力系數(shù)Ct=Cp+Cτ從大于平板開始，逐步下降，在C+≈0.63 左右時，達到平板水平；C+繼續(xù)增加，Ct開始小于平板，并在C+≈1.61,Ct≈0；C+繼續(xù)增加，Ct＜0。對應(yīng)的減阻百分數(shù)η（正值增阻、負值減阻）：C+＜0.63,η ＞0;C+≈0.63,η ≈0;C+＞0.63,η ＜0。

（4）表2～4數(shù)據(jù)說明，海豚表皮僅當(dāng)以流向隨行波形式做精細生物運動時才能實現(xiàn)減阻。3.3.4 減阻效率

（1）有用功率

前面的數(shù)值模擬是通過給定流向隨行波這種柔性壁面的運動形式w( x,t )=?來求解流場，研究阻力成因。如圖11所示，在實際中，總需要力F?( x,t )來驅(qū)動柔性表皮以實現(xiàn)這種運動w( x,t )。無論是海豚還是仿生體，在減阻時，一定會尋求效率的最大化。若來流速度（或前進速度）為u時，表皮做隨行波運動，可減阻Δf，則獲得的有用功率為N0= u ?Δf。

（2）輸入功率

不失一般性，對于流向隨行波可以在任意一時刻t0取一個波長λ 內(nèi)的波面來研究仿生物體輸入功率Ni。

為實現(xiàn)流向隨行波，此時波面上任意一質(zhì)點必須具有z向速度w= aω sin( )kx - ωt0。而這樣變速運動，波面應(yīng)該具有數(shù)值解給出的P( )x,t0分布。所以海豚或仿生物體必須具有驅(qū)動力F( x,t0)dx = P( x,t0)ds ?cos θ( x,t0)= P( x,t0)dx 來維持波面存在。即F( x,t0)= P( x,t0)。

圖11 隨行波柔性面受力分析示意圖Fig.11 Diagram of forces on traveling waveform surface

（3）減阻效率η0

定義

式中，Δf ( λ )為一個波長內(nèi)的減阻值。

表4 給出了在Case 1～Case 9 算例中，c+=0.75 時，減阻效率最大。應(yīng)該注意的是，表4 結(jié)果僅說明：海豚及仿生物體不會或不應(yīng)該一味追求減阻率；確實存在最佳組合，使流向隨行波的減阻效率最高；對于任何一個預(yù)期前進速度，特定物體（已有a+，λ+）應(yīng)該有一個最佳波速c+，使得節(jié)能效率最大，對于仿生物體來說將有一個權(quán)衡的最佳組合（c+，a+，λ+），使得減阻率與減阻效率均比較理想。

4 結(jié) 論

基于前人對海豚表皮微尺度觀測結(jié)果，本文提出了仿生柔性表皮三種可能的數(shù)學(xué)模型：剛性波形、駐波波形及流向隨行波形。通過數(shù)值計算進行了減阻預(yù)報和機理分析，獲得了以下結(jié)論：

（1）RANS方程及動網(wǎng)格技術(shù)較好地描述了二維波形壁面的邊界層流場，預(yù)報了阻力成分；

（2）剛性波、駐波、低波速的流向隨行波摩擦阻力系數(shù)Cτ小于平板的機理是：對右向來流，在波峰→波谷→波峰形成的凹陷區(qū)存在順時針旋轉(zhuǎn)的駐渦，這降低了壁面右向流速度梯度，甚至部分改變了梯度符號。

（3）流向隨行波摩擦阻力系數(shù)Cτ隨波速增加而增加（由小于平板轉(zhuǎn)而大于平板）的機理是：隨行波波面質(zhì)點始終無流向速度，而波面卻以波速在擠壓近壁面流體進行流向運動，從而使原本低波速時存在的駐渦逐漸消失，并使壁面速度呈梯度增加。

（4）流向隨行波壓差阻力系數(shù)Cp隨波速增加逐漸減小，甚至變向成推力的機理是：波速增加，隨行波凹陷區(qū)峰→谷波面質(zhì)點向上運動的速度隨之增加，擠壓流體加劇而提升波面壓力，谷→峰波面質(zhì)點向下運動速度也隨之增加，更快速地讓出空間而使該波面壓力進一步下降。

（5）剛性波、駐波的總阻力系數(shù)Ct都大于平板，流向隨行波隨波速增加Ct一直不斷下降，由大于平板轉(zhuǎn)為小于平板。海豚柔性表皮，唯有流向隨行波形式才有可能減阻，仿生航行體MEMS技術(shù)應(yīng)該要有流向隨行波功能。

（6）針對不同的來流速度（或者前進速度）具有一定參數(shù)特征（a+，λ+）的流向隨行波，存在最佳波速，使得減阻效率最大，進化使海豚巡航時會遵循此原則。而對η0尋優(yōu)，對流向隨行波仿生體進行設(shè)計具有重要參考價值。

（7）從本質(zhì)上講，流向隨行波是推力發(fā)生器，本文對隨行波Cp機理的闡述應(yīng)該同樣適用于蛇形運動等具有隨行波鰭的水下生物。

致謝：在此對惠昌年研究員在本文編寫過程中提出的寶貴意見表示感謝。