多個柔性梁V 型集群行為節能機理研究1)

張磊 敖雷 裴志勇

(武漢理工大學綠色智能江海直達船舶與郵輪游艇研究中心,武漢 430063)

(綠色智能江海直達船舶湖北省工程研究中心,武漢 430063)

引言

集群是動物較為普遍的行為,如大多數鳥類和魚類都會階段性地集群,這種聚集行為產生的有利作用稱為集群效應.生物學[1-4]從動物對資源需求和環境適應等方面闡述集群行為產生的原因,指出集群可以提高族群的防御力、捕食效率和繁殖率.然而,這些定性解釋并不能滿足人類對動物集群的認知需求,例如:遷徙過程中為什么大雁會階段性地使用一字形排列的集群方式;為什么大雁等體型較大的鳥類常采用V 型集群方式,而小型鳥類不采用.亟需定量研究動物集群行為,提升對這一自然現象的認知水平.

近幾十年國內外學者嘗試從受力和能量的角度對動物集群行為和集群效應的內在機理進行定量研究,采用的方法包括理論分析、試驗(活體試驗和仿生學試驗)研究和數值分析.首先,空氣動力學和流體動力學的相關理論被用來研究V 型集群方式下鳥群個體之間和菱形排列魚群個體之間的相對運動,指出集群中跟隨運動的個體可以利用前面個體產生的渦街從而節省能耗[5-6].然后,這種節能觀點也得到了試驗驗證.對V 型排列8 只鵜鶘的心率進行實時監測發現,群體飛行中個體心率較獨自飛行的鵜鶘減少了11.4%～14.5%[7];采用全球定位系統(GPS)對14 只朱鹮的飛行軌跡進行記錄發現,雖然個體位置不斷調整,但群體呈現出比較完整的V 型排列[8];測量鮭鱒魚在均勻流和圓柱尾流后的運動發現,魚可以利用圓柱尾流渦降低自身肌肉的參與度[9].由于活體試驗變量較難控制,且采集的數據(如心率和肌肉參與度等)只是集群節能觀點的側面證據,亟需能夠表征受力和能耗的直接數據.學者們開展了大量的仿生學試驗[10-17],根據仿生簡化方式的不同將這些試驗分為兩類,一類是多個柔性體流致振動的研究,如交錯排列的多個柔性體中領頭柔性梁的受力較單獨流致振動柔性梁的受力減少可以達到50%[13],本類研究沒有考慮柔性梁自主擺動的推進作用;第二類是多剛性體集群行為的研究,如縱向排列2 個剛性翼采用轉動或者升沉方式推進,根據不同的初始縱向間距會自然形成不同的集群模式,且合適的縱向間距下剛性翼的推進效率較單獨推進的剛性翼高[15,17],本類研究忽略了動物柔性的影響.為了更加準確地描述動物集群行為,需要開展多個自推進柔性體自形成集群運動的仿生學試驗,受制于目前的試驗條件,實施起來難度較大.

相對于試驗而言,數值方法求解動物集群這種多個自推進柔性體?流體介質相互作用問題更為經濟,得到的流場細節有利于開展機理分析.文獻[18-19]采用浸入邊界法對均勻流中并排2 個柔性梁的流致振動進行了模擬,本文作者基于Fluent 二次開發提出了流固耦合數值算法,對均勻流中并排2 個柔性梁的流致振動進行了模擬[20],研究均發現不同的橫向間距下柔性梁呈現同向和異向振動兩種模式.同樣,浸入邊界法被用于均勻流中縱向排列2 個柔性梁的流致振動計算[21],后梁受力和幅值較前梁大.文獻[14]采用面元法和振型疊加的聯合仿真算法對均勻流中不同排列(并排、縱向或交錯) 2 個柔性梁的流致振動問題進行討論,得出了不同排列方式下柔性體受力的相關數據.上述研究中柔性體固定在均勻流中,受到流體力產生被動變形,而自然界中鳥類或魚類不僅具有被動變形,它們還通過主動控制自身肌肉的形變(振翅或擺尾等)獲得動力.文獻[22]在流致振動前梁形成的尾流中引入自推進的柔性梁,采用浸入邊界法進行計算發現,后梁在前梁尾流中的運動呈現繞渦模式和穿渦模式,分別對應低阻力和高阻力的狀態.文獻[23-24]采用浸入邊界法對縱向多個柔性梁自推進的流固耦合問題進行求解,研究中柔性梁導邊做同振幅、同頻率的橫向指定運動,最終形成了緊密排列集群模式,多個柔性梁群體較單獨自推進梁具有更大的速度.同樣,浸入邊界法還用于并排2 個柔性體自推進的研究[25-26],根據橫向間距的不同呈現交替跟隨、交替領先和并駕齊驅的集群模式,其中交替跟隨模式下柔性梁的推進速度和效率較高.本文作者采用基于Fluent 二次開發的流固耦合數值算法對并排排列多個(小于4)自推進柔性梁進行了仿真[27],研究中的并排集群方式并不能降低能耗.

對以上案例進行分析可知,首先,浸入邊界法作為一種流固耦合邊界處理方法,主要通過與特定的計算流體力學方法(如格子玻爾茲曼方法)結合用于一些基礎研究,計算流體力學軟件中尚未接入浸入邊界法的相關模塊,使用該方法解決復雜科學問題或者工程問題的門檻較高;然后,動物集群流固耦合問題的相關研究方興未艾,集群模式也多集中在縱向或并排排列的2 個柔性體,缺少對多個柔性梁(大于2) V 型集群行為的研究.鑒于此,本文研究具有兩個目的,一是采用基于Fluent 二次開發的流固耦合數值方法對多個自推進柔性梁自形成的V 型集群過程進行模擬;二是對V 型集群中柔性梁個體的受力和運動參數進行討論,結合流場細節,揭示集群產生的原因和節能的內在機理,特別是對領頭柔性梁的節能機理進行闡述.

1 數值模型和基本原理

1.1 V 型集群數值模型

目前對動物集群行為多進行二維多個柔性體?流體介質相互作用流固耦合問題的研究[23-27],動物個體近似表達為二維柔性梁.多個柔性梁V 型排列集群計算模型如圖1 所示,x軸和y軸分別表示水平方向和垂直方向,多個二維柔性梁(即結構域Ωs)置于流體域Ωf中.本文柔性梁選取的最大數量為5,長為L,厚度為h.由于h?L,柔性梁可看作歐拉伯努利梁.柔性梁左端(即導邊Γl)固定并在y方向做指定的升沉運動y(t)=Acos(2πft),然后該升沉運動與流體的相互作用提供柔性梁水平方向自由運動的動力,最終形成穩定的V 型排列.多個柔性梁同相位升沉運動的幅值A和頻率f保持一致.導邊橫向(即y方向上)間距設置為D,運動過程中保持不變.柔性梁i和柔性梁j導邊縱向(即x方向上) 間距用Hij(t)表示,i,j=1,2,3,4,5.縱向間距初始值設置為Hij(t=0)=H0,在運動過程中縱向間距Hij(t)是隨時間變化的,由多個柔性體與流體之間的相互作用決定.

圖1 多個柔性體自推進V 型排列計算模型示意圖Fig.1 Schematic of multiple self-propelled flexible beams in V configuration

1.2 流固耦合問題控制方程

不可壓縮黏性流動力學控制方程包括質量連續性方程和動量方程,分別為

其中,uf是流體速度矢量,ρf是流體密度,pf表示壓力,μ是動力黏性系數,gf是外部體積力.為了求解流體動力學方程,還需要補充相應的初始條件和邊界條件,分別為

其中,式(2a)表示流體在初始時刻保持靜止,式(2b)表示流固耦合面Γi=Ωf∩Ωs上的速度邊界由柔性體的運動ui決定,式(2c)表示流體域邊界Γw為無滑移邊界條件.流體動力學方程采用有限體積法進行求解.

歐拉伯努利梁在x和y方向上的動力學控制方程分別為

其中,二維柔性梁線密度ρl=ρsh,ρs為柔性體密度,柔性梁位置矢量X(s,t)=(X(s,t),Y(s,t)),流體外力矢量F(s,t)=(Fx(s,t),Fy(s,t)),s表示沿著柔性梁軸線的拉格朗日坐標,柔性梁的楊氏模量為E,轉動慣量為I.運動過程中認為柔性梁自身長度保持不變,補充幾何約束條件如下

同樣地,為了求解柔性梁的動力學方程還需要補充合適的初始條件和邊界條件為

其中,式(4a)和(4b)分別表示初始時刻柔性梁的位置和運動狀態,本文柔性梁初始時刻無形變且保持靜止.式(4c)表示柔性梁導邊Γl做指定的升沉運動,yi是升沉運動平橫位置在oxy坐標系中的坐標值(見圖1).式(4d)表示流體作用在流固耦合面上的力pf(s,t),nf是流體力單位矢量,指向流體外,該流體力通過求解流體動力學方程獲得.柔性梁x軸方向上的動力學方程為二階微分方程.柔性梁導邊在y軸方向上做無相位差的升沉運動y(t)=Acos(2πft),即可將研究的柔性梁看作是導邊具有相同升沉運動激勵的懸臂梁,其y軸方向上的動力學方程可采用廣義支座激勵的振型疊加法進行處理.首先,將懸臂處支座激勵看作是等效載荷的作用,等效載荷處理方法參見文獻[28] 第17 章;然后,由于密度、尺寸、楊氏模量和轉動慣量等參數都相同,則各懸臂梁振型模態都是一致的.柔性梁的微分方程采用4 階龍格庫塔法進行求解.

上述兩場(即流場和結構場)通過耦合面Γi相互作用,從式(3)、式(2b)和式(4d)可以看出,流場產生的作用于結構物上的流體力將改變結構場的動力學條件,反過來結構場運動狀態(即速度)的改變提供影響流場的邊界條件,流固耦合計算中需要保持耦合面上兩場速度和受力的連續性,表達式分別為

其中,ns是結構力單位矢量,指向結構外.根據如式(5a)耦合面上的運動采用動網格技術更新計算域.

1.3 松耦合方法

采用分區算法中的松耦合方法對多個柔性體?流體介質相互作用流固耦合問題進行求解.在單個時間步內,分別對流場和結構場進行一次計算,然后兩場之間數據也只交換一次,具體計算流程如圖2所示,?為表示任意流體參量的通用變量.從時間步tn到tn+1的具體計算流程敘述如下:

圖2 松耦合方法流程圖Fig.2 Flow chart of the proposed loosely-coupled scheme

(1)時間步tn已知參量uf=uf,n,X(s,t)=X(s,tn),?X(s,t)/?t=?X(s,tn)/?t;

(2)時間步tn到tn+1:

(a)求解流體動力學方程積分獲得作用在柔性梁上的流體力pf,n;

(b)求解結構動力學方程獲得tn+1 時刻柔性體的位置矢量X(s,tn+1);

(c)根據結構物位置矢量采用動網格技術更新流體計算域.

2 驗證實例

為了驗證本文數值方法的有效性,對文獻[23]中縱向2 個柔性體自推進的集群行為進行模擬.柔性體導邊做同相位的指定升沉運動y(t)=Acos(2πft),選擇特征長度L、速度Uref=2πAf和流體密度ρf無因次化參量,則柔性體的線密度無因次量表示為=ρl/(ρfL),雷諾數為Re=ρfUrefL/μ,彎曲剛度為K=選擇=0.2,Re=200,K=0.8,A/L=0.5,H0/L=8 的工況進行計算.H0表示2 個柔性體之間的初始縱向間距.為了消除流體域無滑移壁面Γw阻塞效應的影響,選擇足夠大的流體計算域,x和y方向上尺寸分別為[?60L,30L] × [?20L,20L].為了減少網格量采用混合網格對流體域進行空間離散,如圖3 所示.為了準確地獲取柔性梁附近區域的詳細流場信息,網格需要進行局部加密,圖3 中紅色方框內局部放大圖顯示了柔性梁周圍非結構網格加密情況,其中最小網格邊長為?Lmin=0.02L.流場和結構場的計算時間步長選擇統一值?t=0.01 s(約為0.001T),T=1/f為柔性梁指定升沉運動的周期.

圖3 計算域空間離散示意圖Fig.3 Spatial discretization of the fluid domain

圖4 將計算得到的2 個柔性梁速度和縱向間距的時歷曲線與浸入邊界法的結果進行對比可知,兩種方法計算得到的曲線吻合良好.在t/T=0 至2 區間內,縱向排列2 個柔性梁保持相同的速度水平移動,接著后梁突然加速,前后梁的縱向間距逐漸變小直到最后達到一個穩定的間距((H?L)/L=5).

圖4 縱向排列柔性梁導邊運動參數時歷曲線對比:(a) 水平運動速度;(b) 縱向間距Fig.4 Comparison of time history of (a) streamwise velocity and(b) longitudinal separation distance of two tandem beams between the present solutions and Ref.[23]

此外,對數值計算的網格和時間步長無關性進行說明.由于柔性梁動力學方程求解中使用了振型疊加的處理方式,還需要開展模態截斷分析(即模態無關性分析).另外選擇兩套空間離散網格(最小網格邊長分別為?Lmin=0.01L和?Lmin=0.04L)、兩個時間計算步長(?t=0.002 s 和?t=0.05 s)和兩個振型模態數(jmax=3 和jmax=10)進行計算.圖5 對比了集群運動穩定后的柔性梁水平運動速度,仿真計算結果隨著網格尺寸和時間步長的減小單調收斂,隨著柔性梁模態數的增加單調收斂.同時,當網格尺寸選擇?Lmin=0.02L,時間步長?t=0.01 s,模態數jmax=5 時,仿真結果滿足計算的精度需求,在后面V 型集群行為的模擬中采用相同的設置.

圖5 柔性梁導邊水平運動速度時歷曲線對比:(a) 網格無關性分析;(b) 時間步長無關性分析;(c) 模態截斷分析Fig.5 Comparison of streamwise velocity curves of the leading edge in(a) grid independence study,(b) time-step independence study,and(c) mode truncation study

圖5 柔性梁導邊水平運動速度時歷曲線對比:(a) 網格無關性分析;(b) 時間步長無關性分析;(c) 模態截斷分析(續)Fig.5 Comparison of streamwise velocity curves of the leading edge in(a) grid independence study,(b) time-step independence study,and(c) mode truncation study (continued)

3 數值模擬結果及分析

對多個柔性體V 型排列進行仿真,用來模擬自然界中動物V 型集群行為.文獻[29-30]中指出自然界中動物的彎曲剛度系數近似為O(1),這里仍然選用第三節驗證案例中的柔性體,參數為=0.2,Re=200,K=0.8,A/L=0.5,柔性體的數量分別為3 個和5 個.

3.1 3 個柔性梁的穩定V 型排列

3 個柔性梁橫向間距設置為D=0.475L,初始縱向間距設置為H12(t=0)=H13(t=0)=H0=2L.展示了V 型排列3 個自推進柔性梁之間水平速度差值Uij(t)=ui(t)?uj(t)和縱向間距Hij(t)=Xi(s=0,t)?Xj(s=0,t)的時歷曲線,其中ui(t)和Xi(s=0,t)分別代表柔性梁導邊的水平速度和位置.領頭柔性梁1 與柔性梁2 和3 的水平速度差值在t/T=4 后逐漸趨于穩定,繞著零點振蕩,幅值為0.5Uref.柔性梁2 和3 與領頭柔性梁1 之間的縱向間距先增大然后快速減小,t/T=4 后基本保持一致且位于[?0.91L,?1.16L]的區間內作穩定的余弦變化.

圖7 展示了計算穩定后一個升沉運動周期內典型時刻3 個柔性體的集群行為以及相應的渦量和壓力云圖.如圖7(a)～圖7(d)所示,y方向上,柔性梁導邊從其升沉運動的平衡位置(如柔性梁1,Y1(s=0,t=0)=0)向最大負幅值運動(如柔性梁1,=?A,然后向上運動再次經過平衡位置直到t/T=1 時達到最大正幅值處.x方向上,柔性梁2 和3 齊頭并進緊隨領頭柔性梁1 后,這一現象與圖6(b)紅色方框內放大的縱向間距時歷曲線一致.觀察形成的渦量圖,其中渦根據相應梁的編號來命名,渦量值用特征值ωref=Uref/L無因次化,正負號(即 +,-)分別表示渦繞逆時針和順時針旋轉.t/T=1/4 時刻3 個柔性梁上表面導邊處邊界層分離形成的負渦1?,2?和3?向x正向運動,隨后在t/T=3/4 時脫離柔性梁在尾流后混合形成負的卡門渦群(1?+3?和2?).同理,t/T=3/4 時刻3 個柔性梁下表面導邊處邊界層分離形成的正渦1 +,2+和3+向x正向運動,隨后在t/T=1+1/4 時脫離柔性梁在尾流后混合形成正的卡門渦群(1++2+和3 +).由此可見,一個完整升沉運動周期內3 個柔性梁形成的V 型集群模式呈現出一對完整的卡門渦街群.附錄視頻1 展示了3 個柔性體V 型集群一個周期內完整的渦量變化.圖7(e)～圖7(h)給出相應典型時刻的流場壓力云圖,遠場參考壓力設置為零,特征壓力用來無因次化計算得到的壓力值,紅色和藍色分別表示最大正壓力值2 和最小負壓力值?2.為了弄清3 個自推進柔性梁自形成V 型集群行為產生的原因,嘗試結合渦量云圖和流場壓力云圖給出解釋.邊界層分離以及渦產生的地方形成低壓區,例如圖7(a)中,t/T=0/4 時刻3 個柔性體上表面處邊界層發生分離,同時柔性體1 下表面分離的邊界層脫離其尾部,這四個區域如圖7(e)所示均為負壓區,而柔性體下表面具有較大的正壓力,另外,該時刻柔性體從升沉運動平衡位置向下運動產生了向上變形,柔性體上下表面壓力差在水平方向上形成了沿著x負方向的合力,與柔性體水平運動方向相同,該合力推動柔性體加速運動.t/T=1/4 時刻柔性體達到升沉運動的最大負振幅位置處開始向上運動,此時柔性體具有向下的變形(見圖7(f)),結合上表面的負壓和下表面的較大的壓力形成了向著x正方向的合力,該合力與柔性體水平運動方向相反,該阻力使得柔性體減速運動.采用同樣的方法也可以對t/T=1/2 和t/T=3/4 時刻的渦量云圖和壓力云圖進行分析.由以上分析可知,自推進多個柔性梁之間通過流體介質相互作用(如渦的相互混合)產生不同的流場分布,不同的壓力分布和柔性梁形變共同決定了柔性梁的水平受力,從而影響多個柔性體之間的相互運動.圖8 給出了一個升沉運動周期內3 個柔性梁水平受力Fx的時歷曲線,Fx<0 時水平力為推力,柔性梁加速運動;反之,柔性梁作減速運動.后排柔性梁2 和3 受到的推力和阻力幅值較柔性梁1 大.

圖6 3 個柔性梁V 型排列運動參數時歷曲線:(a) 速度差值Uij(t)=ui(t)?uj(t);(b) 縱向間距Hij(t)=Xi(s=0,t)?Xj(s=0,t) (i,j=1,2,3)Fig.6 Time history of (a) velocity difference Uij(t)=ui(t)?uj(t),and(b) longitudinal separation distance Hij(t)=Xi(s=0,t)?Xj(s=0,t) (i,j=1,2,3) of three beams in V configuration

圖7 3 個柔性梁V 型排列流場細節圖:(a～d) 渦量云圖;(e～h) 壓力云圖.特征值ωref=Uref/L 和pref=ρf 分別被用來無因次化計算得到的渦量和壓力Fig.7 Fluid details:(a～d) vortices contour and (e～h) pressure contour of three beams in V configuration.The vortices and pressures are normalized by ωref=Uref/L,pref=ρf ,respectively

圖8 3 個柔性梁水平受力時歷曲線,特征值Fref=(1/2)ρf 被用來無因次化計算得到的合力Fig.8 Hydrodynamic force experienced by three flexible beams in the horizontal direction during one heaving motion period.The forces are normalized by Fref=(1/2)ρf

3.2 5 個柔性梁的穩定V 型排列

5 個柔性梁橫向間距為D=0.475L,初始縱向間距為H12(t=0)=H13(t=0)=H24(t=0)=H35(t=0)=H0=2L.圖9 顯示了V 型排列5 個自推進柔性梁之間水平速度差值和縱向間距的時歷曲線.t/T=4 后,第二排的柔性體2 和3 與領頭柔性體1 之間的速度差值以及第三排柔性體4 和5 與第二排柔性體2 和3 的速度差值逐漸趨于周期性變化,變化周期為3T.同時如圖9(c)和9(d)所示,V 型排列三排柔性體之間的縱向間距時歷曲線也呈周期為3T規律性波動,第二排與領頭柔性體之間的縱向間距位于[?0.98L,?1.76L]區間內,第三排與第二排柔性體之間的縱向間距較小,位于[?0.45L,?1.03L]區間內.這說明5 個柔性體在自推進集群運動過程中雖然整體一直保持V 型排列,但個體之間的相互位置不斷發生著調整,該結論與文獻[8]中對14 只朱鹮V 型集群的動物活體試驗結果得出的結論相符合.

圖10 選取展示了t/T=5.5～8.5 內典型時刻5 個柔性體的集群行為以及相應的渦量云圖.如圖10(a)～圖10(d)所示,y方向上,柔性梁導邊從其升沉運動的最大負幅值(t/T=5.5)向上運動經過平衡位置(t/T=5.75)達到最大正幅值(t/T=6.0)處,然后向下運動經過平衡位置(t/T=6.25),圖10(e)～圖10(h)以及圖10(i)～圖10(l) 柔性體導邊依次重復上述運動.x方向上,第二排柔性體與領頭柔性體縱向間距先變大而后變小,而第三排柔性體與第二排柔性體縱向間距經歷了先變小然后變大的過程,這一現象與圖9(c)和9(d)中縱向間距時歷曲線一致.觀察5 個柔性體V 型集群運動產生的渦量圖可以看出,一個升沉運動周期T內,柔性體上(下)表面邊界層分離依次產生一對負(正)渦,渦在向尾流移動過程中受到多個柔性體之間通過流體介質產生的相互作用而擠壓混合形成了復雜的渦量云圖.對比3 個升沉周期相同時刻(如圖10(a),圖10(e)和圖10(i))的渦量圖發現,渦相互作用的模式基本相似,但由于多個柔性體之間相互位置的差別,造成渦受擠壓或者混合位置以及大小都不盡相同,這些差別將會導致柔性體受力的不同.5 個柔性體V 型集群運動在三個指定升沉運動周期內完整的渦量變化詳見附錄視頻2.

圖9 5 個柔性梁V 型排列運動參數時歷曲線:(a,b) 速度差值Uij(t)=ui(t)?uj(t);(c,d) 縱向間距Hij(t)=Xi(s=0,t)?Xj(s=0,t) (i,j=1,2,3,4,5)Fig.9 Time history of (a,b) velocity difference Uij(t)=ui(t)?uj(t),and(c,d) longitudinal separation distance Hij(t)=Xi(s=0,t)?Xj(s=0,t)(i,j=1,2,3,4,5) of five beams in V configuration

圖10 5 個柔性梁V 型排列流場渦量云圖.特征值ωref=Uref/L 被用來無因次化計算得到的渦量Fig.10 Vortices contour of five beams in V configuration.The vortices are normalized by ωref=Uref/L

圖11(a)～圖11(l)給出了相應時刻的流場壓力云圖,根據3.1 節對壓力場的分析方法,邊界層分離以及渦產生的地方形成低壓區,結合柔性體的形變,可以對作用在柔性體上流體力進行說明,這里不再詳述.圖12 給出了5 個柔性體水平方向受力的時歷曲線,第三排柔性體4 和5 受到的推力(Fx<0)幅值最大,第二排柔性體2 和3 次之,領頭柔性體最小;領頭柔性體受到的阻力(Fx>0)幅值最小,后排柔性體較大.

圖11 5 個柔性梁V 型排列流場壓力云圖.特征值pref=ρf 被用來無因次化計算得到的壓力Fig.11 Pressure contour of five beams in V configuration.The pressures are normalized by pref=ρf

圖12 5 個柔性梁水平受力時歷曲線,特征值Fref=(1/2)ρf L 被用來無因次化計算得到的合力Fig.12 Hydrodynamic force experienced by five flexible beams in the horizontal direction during three heaving motion periods.The forces are normalized by Fref=(1/2)ρf L

3.3 V 型集群性能分析

上述內容對多(3 和5)個柔性體的V 型集群運動進行了描述,接下來引入三個主要參數對多柔性體自推進集群運動性能進行討論[23-25].集群運動水平方向平均推進速度Ui,平均輸入功率Pi和推進效率ηi分別為

其中,下標i表示柔性體的編號,Tc表示多柔性體在升沉運動自推進作用下形成的穩定集群運動周期.觀察圖6 和圖9 的運動參數時歷曲線可知,3 個和5 個柔性體V 型集群水平速度差值和縱向間距曲線分別呈現T和3T的周期性變化,即3 個和5 個柔性體形成的V 型集群運動周期分別為T和3T.需要注意的是,式(6b)中平均輸入功率包括兩部分,一部分產生供應升沉運動需要的能量,一部分提供水平運動的能量.引入文獻[23]中單柔性體自推進運動的數據進行對比,結果見表1.在多個柔性體形成穩定的集群運動后個體水平方向平均推進速度保持一致,即集群具有一個固定的平均推進速度(U=Ui).從表1 可以看出,在柔性體輸入功率基本相同的情況下,V 型排列集群中個體自推進獲得的水平方向速度較單柔性體提升比例超過15%,雖然5 個柔性體集群中第二排柔性體2 和3 需要更大的輸入功率,但是其推進效率較單個柔性體提升比例超過14.5%.該對比包含兩方面的內容:(1)在相同能量的輸入下,集群中個體可以獲得更高的推進速度;(2)即使集群中某些個體存在高能量輸入的情況,但其推進的效率仍然較高.本文的計算結果給出了V 型排列集群行為具有降低能耗的直接數據證據.此外值得注意的是,本文研究的V 型集群中領頭柔性梁的運動速度和效率均較單獨柔性梁高,且是集群中效率最高的個體.這一結論打破了集群中領頭“羊”位置最耗能的認知.

表1 多個自推進柔性體性能參數比較Table 1 Propulsive performance of multiple self-propelled flexible beams

接下來嘗試對V 型集群行為節能機理進行闡述.文獻[23-25,31]等對集群中后排柔性體節能的原因做出了解釋,認為后排柔性體可以利用前排柔性體尾流中的渦,從本文3.1 節和3.2 節對集群流場(渦量云圖和壓力云圖)的分析可以得出同樣的結論,即前排柔性體尾渦在后排柔性體前端形成了低壓區,造成柔性體前后具有更大的壓差(如圖10(f)和圖10(j),圖11(f)和圖11(j)中柔性體3 周圍的渦量圖和壓力場所示),從而為水平運動提供更大的推力.

這里需要重點對領頭柔性梁節能的機理進行闡述,圖13(a)對比了領頭柔性梁的水平受力時歷曲線,圖13(b)給出了領頭柔性體尾部和導邊橫向位移差值(即橫向形變,Δy=Y(s=L,t)?Y(s=0,t))時歷曲線的對比,可以看出,V 型集群領頭柔性梁水平方向受力的幅值較單柔性體大,且橫向形變幅值較小.圖14分別給出了單柔性體和3 柔性體V 型集群的流場細節圖.t/T=5.25 時刻,3 柔性體V 型集群運動領頭柔性梁尾部脫落的渦受到第二排柔性梁2 的限制,沿著柔性梁2 的下表面運動,與單柔性體壓力云圖對比可知,其尾渦形成的低壓區下移(如圖14(d) 和14(e)中紅色方框所示),且領頭柔性梁的尾部上表面處壓力較大,這導致t/T=5.25 時刻集群中領頭柔性體的最大橫向形變較單柔性體小.t/T=5.25 和t/T=5.3 時刻分別是單柔性體受到推力和3 個柔性體中領頭柔性體受到推力最大的時刻,對比圖14(a),圖14(d)和圖14(c),圖14(f)可以看出,由于領頭柔性體尾渦與第二排柔性梁2 和3 導邊邊界層分離渦的相互作用,集群領頭柔性體上下表面壓差增大,所以雖然t/T=5.3 時刻領頭柔性體的橫向形變較小,但是其水平方向的推力仍然比t/T=5.25 時刻單柔性體的推力略大.從上述分析得出,集群中領頭柔性體與后排柔性體通過流體介質發生相互作用(如領頭柔性梁尾渦與后排柔性梁邊界層分離渦相互作用造成流場變化),一方面使得柔性體橫向形變幅值變小,另一方面造成柔性體上下表面壓力差變大.在這兩方面因素(更小的迎水面和更大的推力)的共同作用下,領頭柔性體獲得較高的運動速度和推進效率.

圖13 領頭柔性梁對比:(a) 水平受力時歷曲線;(b) 橫向形變時歷曲線.特征值Fref=(1/2)ρf L 被用來無因次化計算得到的合力Fig.13 The comparison of the leading beams:(a) hydrodynamic force in the horizontal direction and (b) lateral deformation during one heaving motion period.The forces are normalized by Fref=(1/2)ρf L

圖14 單柔性體和3 個柔性梁流場細節圖:(a～c) 渦量云圖;(d～f) 壓力云圖.特征值ωref=Uref/L 和pref=ρf 分別被用來無因次化計算得到的渦量和壓力Fig.14 Fluid details of single beam and 3 beams:(a-c) vortices contour and (d-f) pressure contour.The vortices and pressures are normalized by ωref=Uref/L,pref=ρf ,respectively

4 結論

本文采用基于Fluent 二次開發的數值方法求解多個柔性體?流體介質相互作用的V 型集群問題.首先,使用該算法對縱向2 個柔性體的二維集群問題進行計算,將結果與浸入邊界法仿真結果進行對比說明松耦合方法的可靠性.然后,對自然界中常見的動物V 型集群行為進行模擬,分別計算了3 柔性體和5 柔性體升沉運動自推進的二維集群行為.在4 個升沉運動周期(即t/T=4)后,多個柔性體都形成了穩定的V 型排列,柔性體間橫向間距保持一定,縱向間距在一定區域范圍內呈現周期性變化,如5 柔性體V 型集群中第二排柔性體與領頭柔性體縱向間距穩定在[?0.98L,?1.76L]區間內.該計算結果與文獻[8]中朱鹮V 型集群的動物活體試驗結論一致,即動物在自推進集群運動過程中雖然整體一直保持V 型排列,但個體之間的相互位置不斷發生著調整.引入水平方向平均推進速度、輸入功率和推進效率三個參數對集群中個體推進性能進行對比可知,不僅后排柔性梁的速度和推進效率得到提升,領頭柔性梁也得到大幅提升,其中推進速度較單獨柔性體自推進速度提升比例超過15%,推進效率提升比例也超過14.5%,該計算結果給出了V 型排列集群行為降低能耗的定量直接證據.值得注意的是,該V 型集群運動中領頭柔性梁是集群中速度和效率最高的個體.這一結果打破了集群中領頭“羊”位置最耗能的普遍認知.

為了弄清集群行為產生的原因和不同位置個體的節能機理,對數值仿真得到的二維集群流場細節進行了詳細分析.多個柔性體通過流體介質發生相互作用(如邊界層分離形成的渦相互擠壓混合)造成壓力場的變化,在此壓力場和柔性體自身形變的共同作用下,多個柔性體逐漸形成了穩定的V 型集群排列.V 型集群中不同位置的個體采用不同策略提升推進性能.后排柔性梁節能機理與其他研究者給出的結論一樣,即其可利用前排個體的尾流渦進而節省能耗;本文還對集群中領頭柔性梁的節能機理進行了闡述,即在自身產生的尾渦與后排柔性體的相互作用下,領頭柔性梁自身形變變小和所受推力幅值變大.

為了使多個柔性體的集群運動更加接近自然界中豐富多樣的動物集群現象,研究具有不同長度和彎曲剛度的多個柔性梁采用不同指定運動(升沉、轉動或兩者相結合)方式自推進進而自發形成的集群運動行為將是下一步的工作方向.