張磊 敖雷 裴志勇

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

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

引言

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

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

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

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

1 數(shù)值模型和基本原理

1.1 V 型集群數(shù)值模型

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

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

1.2 流固耦合問題控制方程

不可壓縮黏性流動力學(xué)控制方程包括質(zhì)量連續(xù)性方程和動量方程,分別為

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

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

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

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

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

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

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

其中,ns是結(jié)構(gòu)力單位矢量,指向結(jié)構(gòu)外.根據(jù)如式(5a)耦合面上的運動采用動網(wǎng)格技術(shù)更新計算域.

1.3 松耦合方法

采用分區(qū)算法中的松耦合方法對多個柔性體?流體介質(zhì)相互作用流固耦合問題進(jìn)行求解.在單個時間步內(nèi),分別對流場和結(jié)構(gòu)場進(jìn)行一次計算,然后兩場之間數(shù)據(jù)也只交換一次,具體計算流程如圖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)求解流體動力學(xué)方程積分獲得作用在柔性梁上的流體力pf,n;

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

(c)根據(jù)結(jié)構(gòu)物位置矢量采用動網(wǎng)格技術(shù)更新流體計算域.

2 驗證實例

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

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

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

圖4 縱向排列柔性梁導(dǎo)邊運動參數(shù)時歷曲線對比:(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]

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

圖5 柔性梁導(dǎo)邊水平運動速度時歷曲線對比:(a) 網(wǎng)格無關(guān)性分析;(b) 時間步長無關(guān)性分析;(c) 模態(tài)截斷分析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 柔性梁導(dǎo)邊水平運動速度時歷曲線對比:(a) 網(wǎng)格無關(guān)性分析;(b) 時間步長無關(guān)性分析;(c) 模態(tài)截斷分析(續(xù))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 數(shù)值模擬結(jié)果及分析

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

3.1 3 個柔性梁的穩(wěn)定V 型排列

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

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

圖6 3 個柔性梁V 型排列運動參數(shù)時歷曲線:(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 型排列流場細(xì)節(jié)圖:(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 個柔性梁的穩(wěn)定V 型排列

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

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

圖9 5 個柔性梁V 型排列運動參數(shù)時歷曲線:(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)給出了相應(yīng)時刻的流場壓力云圖,根據(jù)3.1 節(jié)對壓力場的分析方法,邊界層分離以及渦產(chǎn)生的地方形成低壓區(qū),結(jié)合柔性體的形變,可以對作用在柔性體上流體力進(jìn)行說明,這里不再詳述.圖12 給出了5 個柔性體水平方向受力的時歷曲線,第三排柔性體4 和5 受到的推力(Fx<0)幅值最大,第二排柔性體2 和3 次之,領(lǐng)頭柔性體最小;領(lǐng)頭柔性體受到的阻力(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 型集群性能分析

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

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

表1 多個自推進(jìn)柔性體性能參數(shù)比較Table 1 Propulsive performance of multiple self-propelled flexible beams

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

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

圖13 領(lǐng)頭柔性梁對比:(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 個柔性梁流場細(xì)節(jié)圖:(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 結(jié)論

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

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

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