鯨豚表皮微尺度柔性變形的減阻原理研究

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

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

0 引 言

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

國外在上世紀(jì)六十年代就開展了仿生柔性面的減阻研究，包括實驗及數(shù)值計算方法研究[8-11]。我國從上世紀(jì)末開始柔性面的理論和試驗研究，如中國船舶科學(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ù)的選取時，本文參考了國外文獻(xiàn)中的海豚表皮觀測結(jié)果。雖然用肉眼觀察比較模糊，但是在很多齒鯨（如海豚）的皮膚表面確實分布著波紋狀的規(guī)則溝嵴。這些溝嵴較淺，在海豚死后會因為皮膚松弛而模糊或消失不見，同時很難在皮膚的組織切片中觀察到，文獻(xiàn)中采用快速固化方法對活體齒鯨亞目動物的皮膚嵴進(jìn)行了壓印和測量[15]。

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

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

2 計算方法

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

文中數(shù)值計算涉及的柔性變形依靠動網(wǎng)格技術(shù)實現(xiàn)。對于動網(wǎng)格，在任意一控制體?V（邊界是運(yùn)動的）上，標(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ī)則溝嵴，用駐波和流向隨行波來代替海豚皮膚可能的柔性變形，尚不考慮其表面的彈性影響。

剛性波面表達(dá)式為

駐波柔性面表達(dá)式為

流向隨行波柔性面表達(dá)式為

3 模擬結(jié)果

為了探尋可能的仿生減阻方案，本文對剛性的及柔性的波狀壁面進(jìn)行了阻力數(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 流場分析

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

3.2.1 剛性波面ζ = acoskx

圖5給出了對應(yīng)的速度云圖。右向流動的邊界層內(nèi)，由波峰→波谷→波峰形成的壁面凹陷，產(chǎn)生了流道周期性收縮擴(kuò)散區(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)部分壁面速度梯度比平板小，甚至梯度變負(fù)。這樣，即使在波峰處，速度梯度高于平板，但其總粘性阻力將比平板低。由于順時針駐渦的存在，流向最大流速不在剛性波面的波谷處，而是在波谷稍后的下游處，所以速度波狀等值線的“谷”仍然對應(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ū)。故其速度等值波狀線與運(yùn)動的駐波面反相。

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

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

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

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

在大地坐標(biāo)系中觀察波峰→波谷→波峰形成的凹陷區(qū)。對于剛性波面，位置、大小都是不變的；對于駐波面，位置不變，波高呈周期性正負(fù)變化。雖然隨行波以波形的波速c=流向移動，但波面質(zhì)點x位置固定不變。只是位于“波峰→波谷”波面上的質(zhì)點總是以有規(guī)律的正z向速度驅(qū)趕流體，使“峰→谷”波面以流向速度c擠壓凹陷區(qū)中的流體空間，而凹陷區(qū)中的“谷→峰”波面上質(zhì)點以一定規(guī)律沿負(fù)z向速度運(yùn)動，使谷→峰波面以流向速度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ū)中的駐渦強(qiáng)度。只要波速足夠大，駐渦強(qiáng)度會低于駐波、剛性波，乃至消失，并不斷使壁面正向速度呈梯度增加。所以可預(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é)流動物理過程及機(jī)理分析，將重要參數(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 左右時，達(dá)到平板水平；C+繼續(xù)增加，Ct開始小于平板，并在C+≈1.61,Ct≈0；C+繼續(xù)增加，Ct＜0。對應(yīng)的減阻百分?jǐn)?shù)η（正值增阻、負(fù)值減阻）：C+＜0.63,η ＞0;C+≈0.63,η ≈0;C+＞0.63,η ＜0。

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

（1）有用功率

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

（2）輸入功率

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

為實現(xiàn)流向隨行波，此時波面上任意一質(zhì)點必須具有z向速度w= aω sin( )kx - ωt0。而這樣變速運(yùn)動，波面應(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ù)期前進(jìn)速度，特定物體（已有a+，λ+）應(yīng)該有一個最佳波速c+，使得節(jié)能效率最大，對于仿生物體來說將有一個權(quán)衡的最佳組合（c+，a+，λ+），使得減阻率與減阻效率均比較理想。

4 結(jié) 論

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

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

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

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

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

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

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

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

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