蔡琛芳,梁 彬

(中國航天空氣動力技術研究院,北京 100074)

0 引 言

自然界中樹葉和某些植物種子的飄落，撲克牌、紙片等輕質平板的掉落，以及水下氣泡上升等常見現(xiàn)象，都會呈現(xiàn)規(guī)律、美麗的運動軌跡。這些有趣的現(xiàn)象看似簡單，卻包含了復雜的流體力學和動力學機理，近年來引起了人們更多的關注和研究[1-11]。Belmont等[1]在三種不同流體中進行了紙片自由下落實驗，確定了紙片自由下落中擺動和翻滾兩種運動軌跡互相轉變的臨界弗勞德數(shù)Fr。Wang等[2-3]對二維平板的下落問題進行了實驗和計算研究，建立了相應的力學模型，解釋了其非定常氣動力機制。Wan等[4]、Tian和Shu[5-6]、Kumota和Mochizuki[7]分別采用實驗和數(shù)值模擬方法研究了平板自由下落中的流動渦結構和發(fā)展過程，并分析了其中的致力機理。Hu和Wang[8]在動力學系統(tǒng)方程的基礎上，引入微分方程分支理論，分析了平板自由下落的穩(wěn)定性問題。Zhong等[9-10]還使用立體視覺的方法，首次研究了三維、六自由度圓盤的自由下落，并引入了螺旋運動和平動兩種運動形式。這些研究結果表明，空氣中紙片或樹葉自由下落的雷諾數(shù)為1000左右，與較大昆蟲飛行的雷諾數(shù)相似[12]。研究平板自由下落問題，對樹葉和種子的飄落、昆蟲飛行等自然現(xiàn)象，以及仿生微型飛行器等工程應用的相關物理問題中的流體力學和動力學研究具有一定的指導意義。

綜上所述，以往的研究主要以實驗手段為主，并且均采用勻質平板作為研究對象。本文在上述研究的基礎上，針對輕質平板自由下落的簡化模型——二維平板自由下落這一經(jīng)典問題，進一步探討轉動慣量、初始角度和非勻質平板質心位置對下落運動的影響，采用耦合求解N-S方程和運動方程的方法數(shù)值模擬進行系統(tǒng)的研究，并分析其中的流體力學和動力學機制。

1 方 法

1.1 控制方程

平板自由下落過程中，受到氣動力和重力的共同作用，其運動受N-S方程和運動方程共同支配，兩者共同構成控制方程。

在二維慣性坐標系Oxy下，非定常不可壓N-S方程的無量綱形式為：

(1)

(2)

式中，u是無量綱速度，p是無量綱壓力，τ是無量綱時間，是梯度，2是拉普拉斯算子，Re是雷諾數(shù)。

二維平板在慣性坐標系下的運動方程如下：

(3)

式中，u、v分別是平板在x軸(水平方向)和y軸(垂直方向)上的運動速度；ω和θ是平板轉動速度和轉動角度；XA、YA和MA分別是x軸、y軸上的氣動力分量和氣動力矩；m是平板質量；I是平板轉動慣量；g是重力加速度。

1.2 計算方法

N-S方程的數(shù)值求解方法與文獻[13]中的方法相同，采用擬壓縮性算法求解[14-15]。算法和離散方法本文不再贅述。其中在求解動量方程時，引入了虛擬時間步來保證速度散度為0的不可壓條件。計算時采用剛性動網(wǎng)格技術，網(wǎng)格為貼體正交結構，計算過程中網(wǎng)格隨平板剛性運動。所用網(wǎng)格密度為60×49(徑向×周向)，時間步長Δτ=0.0005。

運動方程中需要通過求解N-S方程獲取氣動力和力矩項，而N-S方程中也需要通過求解運動方程獲取邊界條件，因此求解二維平板的自由下落時必須對兩者進行耦合求解。這里采用了 “弱耦合”求解[16]的方法，具體步驟如下：

1)假設已知平板tn時刻的狀態(tài)變量Λn=(un,vn,ωn,θn)，因此確定N-S方程的邊界條件，通過前文的求解方法數(shù)值求解N-S方程得到tn時刻的流場和氣動力、氣動力矩等信息；

2)采用預估校正歐拉法推進求解運動方程。首先用歐拉法估計tn+1時刻狀態(tài)變量的值(上標*表示估計值，f表示方程組(3)的右端函數(shù))：

(4)

3)修正tn+1時刻狀態(tài)變量的值，即tn+1時刻最終值：

(5)

4)根據(jù)解出的tn+1時刻的平板狀態(tài)變量，重復第1步的步驟。如此反復，推進求解。

綜上，可實現(xiàn)運動方程和N-S方程的耦合求解。

1.3 無量綱方法

2 數(shù)值模擬結果

結合之前的研究成果，本文選擇寬厚比β=14、無量綱質量m+=0.3857的平板在水中的下落作為研究對象，取ρ=1.0 g/cm3，υ=0.0089 cm2/s。根據(jù)估算的平均下落速度，Re=2094。

2.1 實驗結果的模擬

首先采用實驗結果對數(shù)值模擬方法進行驗證。圖1和表1給出了二維平板自由下落軌跡的數(shù)值模擬結果和實驗結果。其中實驗結果來源于Anderson等[2]使用擬二維平板在水槽中的實驗：H=0.081 cm、β=14、I=0.026 g·cm2、Re=1147。數(shù)值模擬中采用了相應的無量綱參數(shù)。

圖1 擺動運動下落軌跡Fig.1 Trajectories of fluttering

表1 擺動運動平均速度Table 1 Average velocities of fluttering

同時，為了研究平板初始角度對下落運動的影響，分別選取初始角度0°、45°、90°和135°進行了計算。根據(jù)數(shù)值模擬結果，其他參數(shù)不變的情況下，不同初始角度平板進入穩(wěn)定下落運動的時間長短不一，但并不影響其穩(wěn)定后的運動。在后面的數(shù)值模擬當中，都采用了更快進入穩(wěn)定、規(guī)律的下落運動狀態(tài)的初始角度(一般為45°或135°)。

2.2 轉動慣量的影響

為研究轉動慣量對平板下落運動的影響，在其他參數(shù)不變情況下，對無量綱轉動慣量I+=0.001～5的平板進行了數(shù)值模擬。圖2和圖3分別顯示了幾個典型狀態(tài)下落過程運動軌跡和無量綱速度變化曲線。

數(shù)值模擬結果顯示，無量綱轉動慣量I+≈0.95是下落運動轉變的臨界點。當平板轉動慣量較小(I+<0.95，圖2a、圖2b時，下落運動呈現(xiàn)規(guī)律的擺動現(xiàn)象。隨著轉動慣量的增大,擺動的頻率減小、周期增大。當平板轉動慣量較大(I+>0.95，圖2c)時，下落運動轉變?yōu)橐?guī)律的翻滾，轉動慣量變化對其運動影響并不明顯。

(a)高頻擺動 (b)低頻擺動 (c)翻滾圖2 平板自由下落運動軌跡Fig.2 Trajectories of freely falling plates

(a)高頻擺動

(b)低頻擺動

(c)翻滾圖3 無量綱速度變化曲線Fig.3 Dimensionless velocities

圖2和圖3顯示，兩種運動都呈現(xiàn)了周期性的速度和位移變化，但區(qū)別在于擺動運動的速度和位移變化是對稱的，而翻滾運動在水平方向和角度方向上體現(xiàn)了單側的偏向性。兩種運動的特點在水平速度和垂直速度u+-v+曲線(圖4)中得到明顯體現(xiàn)：擺動運動為“∞”形軌跡，且周期越長“∞”形軌跡越大；翻滾運動則呈單側的“O”形軌跡。

2.3 質心的影響

前文的研究當中，認為二維平板為勻質，即質心位于平板中點。這里首先選取I+=0.0323的擺動下落狀態(tài)，對非勻質,即不同的質心位置(設質心與平板中點距離為l，l=0～0.25L),平板進行了數(shù)值模擬研究，如圖5和圖6所示。

圖4 無量綱速度u+-v+曲線Fig.4 u+ versus v+

圖5 質心偏移時下落運動軌跡Fig.5 Trajectories of freely falling plates (l>0)

圖6 質心偏移動時無量綱速度u+-v+曲線Fig.6 u+ versus v+ (l>0)

結果顯示，當質心位置稍偏離平板中點時(如l=0.01L，圖5a)，平板的擺動運動出現(xiàn)了較小的不對稱，這一不對稱也同樣體現(xiàn)在速度的周期性變化中(圖6)。隨著質心位置逐漸偏離(如l=0.04L，圖5b)，擺動運動的不對稱越來越大，擺動在一側的特征越來越不明顯，而另一側反之；當l=0.15L時(圖5d)，一側的擺動已經(jīng)消失，此時的平板下落運動已經(jīng)近似翻滾運動。這個單側擺動逐漸消退的變化過程在u+-v+曲線有明顯體現(xiàn)(圖6)——隨著質心位置的偏離，“∞”形軌跡的左側逐漸萎縮直至消失，退化為“o”形軌跡。特別是在擺動運動和翻滾運動轉變的臨界點附近(如l=0.13L)，u+-v+曲線上雖然還殘存著“∞”形運動特征(圖6)，但從運動軌跡上已經(jīng)看不出搖擺運動(圖5c)。也就是說，在質心位置偏離平板中點時，下落運動是一種特殊的擺動運動，根據(jù)其運動特征可稱為“非對稱擺動”。

圖7(a)顯示了I+=1.0、質心偏移l=0.01L時的平板自由下落運動軌跡。可見在翻滾運動的狀態(tài)下，發(fā)生一定的質心偏移將會衍生出一種擺動和翻滾交替出現(xiàn)的復雜運動，根據(jù)其運動特征可稱為“擺動-翻滾混合運動”，可簡稱“擺翻運動”。擺翻運動的u+-v+曲線呈現(xiàn)雙“∞”形軌跡(圖7b)。

(a)運動軌跡

(b)u+-v+曲線圖7 擺翻運動Fig.7 Flutter-tumbling

當質心位置偏離平板中點較遠時(如l=0.25L)，平板保持豎直穩(wěn)定的定常下落，未出現(xiàn)速度和位移的周期性變化。

3 討 論

首先以質心位于平板中心、I+=0.0323的擺動狀態(tài)為例，對平板下落過程進行氣動力分析，圖8顯示了該狀態(tài)下一個擺動運動周期中的無量綱速度、角度和氣動力隨時間的變化曲線。圖8中可見，當平板擺動到最左側時(tw=0)，存在一定的垂直速度而水平速度為0，角度約65°，此時平板狀態(tài)相當于攻角25°,來流速度等于垂直速度。此時氣動力方向偏右上，作用點偏平板右側，相對質心產(chǎn)生正(逆時針)的氣動力矩(圖8)。同時，平板受到重力作用，氣動力和重力的共同作用下使平板向右下方平動的同時逆時針轉動(tw=0～0.25)。tw=0.25～0.5時同理也能產(chǎn)生相似的作用。反之，當tw=0.5～1時，平板的狀態(tài)正好相反，產(chǎn)生了相反的作用。圖9將這一作用過程詳細直觀地進行了描述，給出了典型位置的總氣動力大小、方向和作用位置(括號中的值代表總氣動力距平板質心力臂，下同)。在整個擺動過程中由于平板轉動慣量較小，氣動力矩產(chǎn)生較大的轉動加速度，平板的轉動能夠較快地響應氣動力矩變化，循環(huán)往復形成周期擺動。

圖8 擺動運動時無量綱速度和氣動力變化曲線Fig.8 Dimensionless velocities and aerodynamic force components of fluttering

再來分析翻滾運動的情況，圖10以I+=1.0時為例顯示了翻滾下落時的平板受力情況。平板向左側運動的過程中，具有較大的水平速度，此時氣動力略大于重力且方向相反，平板在慣性作用下保持平動且加速了順時針轉動，這一過程同擺動運動中相似(圖9)。當平板運動到最左端時，受力情況也與擺動運動中類似(圖9)，但由于平板轉動慣量較大，逆時針方向的氣動力矩僅僅減緩了平板的順時針轉動，因此平板仍然保持順時針轉動并逐漸回到上一狀態(tài)，如此往復形成周期性翻滾運動。

圖9 擺動運動時平板受力情況示意圖Fig.9 Plates attitude,forces,and velocities of fluttering

圖10 翻滾運動時平板受力情況示意圖Fig.10 Plates attitude,forces,and velocities of tumbling

接下來對非對稱擺動的自由下落情況進行受力分析。圖11中以圖5(b)中非對稱擺動下落時的狀態(tài)為例顯示了平板受力變化情況，其中基本過程與擺動運動一致，不同點在于：由于質心位置偏移，不僅重力的作用點發(fā)生了變化，同時導致了氣動力相對質心力矩的左右不對稱。在平板向右側運動的過程中，氣動力作用點距質心較近，氣動力矩較小作用微弱，故擺動的速度、距離和時間較長；反之平板向左側運動時，氣動力作用點距質心較遠，氣動力矩較大作用明顯，故擺動過程較短(圖11)。這個過程周期變化形成了非對稱擺動現(xiàn)象。

同理，擺翻運動看似復雜，但受力情況與擺動和翻滾運動類似，這里不再分析，簡述如下：在擺動運動部分氣動力作用點距質心較遠，氣動力矩作用大于轉動慣性作用，造成了擺動運動；在翻滾運動部分氣動力作用點距質心較近，氣動力矩作用小于轉動慣性作用，造成了翻滾運動。

最后，質心偏離平板中心較遠時的定常下落可視為一種非對稱擺動或擺翻運動的極限情況：擺動較短的過程(或擺動部分)消失，擺動較長的過程(或翻滾部分)增至無限。此狀態(tài)中平板保持豎直下落，實際攻角和轉動速度為0，氣動力和重力作用于一條直線上，平板在合力作用下保持定常下落狀態(tài)。

圖11 非對稱擺動運動時平板受力情況示意圖Fig.11 Plates attitude,forces,and velocities of asymmetrical-fluttering

4 結 論

本文利用動網(wǎng)格技術，耦合求解N-S方程和運動方程，對二維平板的自由下落問題進行了研究，數(shù)值模擬出了以往實驗研究結果[1,2,5,11]中出現(xiàn)的兩種規(guī)律運動：擺動運動和翻滾運動。主要結果如下：

1)無量綱轉動慣量和平板質心位置是決定平板自由下落運動方式的重要參數(shù)。其中，無量綱轉動慣量I+≈0.95是擺動運動和翻滾運動轉變的臨界點；

2)非勻質平板質心位置偏離中點時，平板自由下落出現(xiàn)新的、衍生于擺動運動和翻滾運動的規(guī)律下落運動，根據(jù)其運動特點可稱為非對稱擺動和擺翻運動。