平板自由下落的數(shù)值模擬

蔡琛芳,梁 彬

(中國航天空氣動(dòng)力技術(shù)研究院,北京 100074)

0 引 言

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

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

1 方 法

1.1 控制方程

平板自由下落過程中，受到氣動(dòng)力和重力的共同作用，其運(yùn)動(dòng)受N-S方程和運(yùn)動(dòng)方程共同支配，兩者共同構(gòu)成控制方程。

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

(1)

(2)

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

二維平板在慣性坐標(biāo)系下的運(yùn)動(dòng)方程如下：

(3)

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

1.2 計(jì)算方法

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

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

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

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

(4)

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

(5)

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

綜上，可實(shí)現(xiàn)運(yùn)動(dòng)方程和N-S方程的耦合求解。

1.3 無量綱方法

2 數(shù)值模擬結(jié)果

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

2.1 實(shí)驗(yàn)結(jié)果的模擬

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

圖1 擺動(dòng)運(yùn)動(dòng)下落軌跡Fig.1 Trajectories of fluttering

表1 擺動(dòng)運(yùn)動(dòng)平均速度Table 1 Average velocities of fluttering

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

2.2 轉(zhuǎn)動(dòng)慣量的影響

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

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

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

(a)高頻擺動(dòng)

(b)低頻擺動(dòng)

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

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

2.3 質(zhì)心的影響

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

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

圖5 質(zhì)心偏移時(shí)下落運(yùn)動(dòng)軌跡Fig.5 Trajectories of freely falling plates (l>0)

圖6 質(zhì)心偏移動(dòng)時(shí)無量綱速度u+-v+曲線Fig.6 u+ versus v+ (l>0)

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

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

(a)運(yùn)動(dòng)軌跡

(b)u+-v+曲線圖7 擺翻運(yùn)動(dòng)Fig.7 Flutter-tumbling

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

3 討 論

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

圖8 擺動(dòng)運(yùn)動(dòng)時(shí)無量綱速度和氣動(dòng)力變化曲線Fig.8 Dimensionless velocities and aerodynamic force components of fluttering

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

圖9 擺動(dòng)運(yùn)動(dòng)時(shí)平板受力情況示意圖Fig.9 Plates attitude,forces,and velocities of fluttering

圖10 翻滾運(yùn)動(dòng)時(shí)平板受力情況示意圖Fig.10 Plates attitude,forces,and velocities of tumbling

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

同理，擺翻運(yùn)動(dòng)看似復(fù)雜，但受力情況與擺動(dòng)和翻滾運(yùn)動(dòng)類似，這里不再分析，簡述如下：在擺動(dòng)運(yùn)動(dòng)部分氣動(dòng)力作用點(diǎn)距質(zhì)心較遠(yuǎn)，氣動(dòng)力矩作用大于轉(zhuǎn)動(dòng)慣性作用，造成了擺動(dòng)運(yùn)動(dòng)；在翻滾運(yùn)動(dòng)部分氣動(dòng)力作用點(diǎn)距質(zhì)心較近，氣動(dòng)力矩作用小于轉(zhuǎn)動(dòng)慣性作用，造成了翻滾運(yùn)動(dòng)。

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

圖11 非對稱擺動(dòng)運(yùn)動(dòng)時(shí)平板受力情況示意圖Fig.11 Plates attitude,forces,and velocities of asymmetrical-fluttering

4 結(jié) 論

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

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

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