自由落體下落速度振蕩現(xiàn)象的研究

李軼明

(中國石油大學(xué)(北京) 石油工程學(xué)院,北京 102249)

0 引 言

固體在流體中的下落過程是一個經(jīng)典的流體力學(xué)問題，同時在眾多工業(yè)領(lǐng)域中具有實際應(yīng)用價值。這一過程是流體與固體之間的受力耦合過程，會產(chǎn)生十分復(fù)雜的現(xiàn)象[1-2]，而物體的下落速度是描述問題的關(guān)鍵參數(shù)。物體在下落過程中，受到重力和流體阻力的作用，當(dāng)阻力和重力達到平衡時，降落速度將不會再增加。但這個速度是否穩(wěn)定則取決于物體形狀、運動狀態(tài)和物理特性、繞流場特性及流體物性等因素，其下落速度并不一定是一個固定值。Joseph和Liu[3]測量了圓柱以不同角度在不同流體中的自由下落過程，記錄了圓柱在水平和垂直方向上的運動軌跡，發(fā)現(xiàn)圓柱下落的過程與圓柱的尺寸、重量以及流體特性有關(guān)。Liu和Nelson等人[4]測量了不同流體中圓球沿傾斜平面的運動規(guī)律，實驗中發(fā)現(xiàn)與粘性流體不同的是在粘彈性流體中圓球出現(xiàn)了異常翻滾現(xiàn)象。Mollinger和Cornelissen等人[5]通過實驗研究了非牛頓流體中球形體下落過程中速度的振蕩現(xiàn)象，在球體下落過程中，其運動速度會不斷減低，但是當(dāng)速度降低到一定值之后突然增加，這與非牛頓流體的物性參數(shù)有關(guān)。

本文通過實驗和數(shù)值模擬的手段研究多種形狀的小型剛體在空氣中的下落過程，實驗?zāi)Ｐ驮趯ΨQ性上有較大差異，以驗證物體對稱性對下落速度穩(wěn)定性的關(guān)鍵作用。實驗在立式風(fēng)洞中進行，基于相對運動原理，物體相對于風(fēng)洞試驗段(即地面)的速度疊加與試驗段氣流相對(即相反)的速度，即可模擬物體下落時相對地面的速度。例如物體在立式風(fēng)洞向上的勻速氣流中對于試驗段處于懸浮(或稱懸停)狀態(tài)，就等同于物體在靜止空氣中相對于地面以該速度勻速下落的狀態(tài)；而如果物體對風(fēng)洞試驗段(即地面)存在相對運動，即相當(dāng)于該物體在靜止大氣中做變速降落運動。在下述討論中，本文將依這種相對運動的關(guān)系用立式風(fēng)洞的試驗現(xiàn)象來研究物體在近地面的靜止大氣中自由降落時的速度變化規(guī)律。這里強調(diào)近地條件，是說明物體的重力加速度不變，空氣的密度不變。

1 自由下落物體速度振蕩現(xiàn)象實驗研究

1.1風(fēng)洞立式試驗段

試驗在北京大學(xué)湍流與復(fù)雜系統(tǒng)國家重點實驗室一號邊界層風(fēng)洞立式試驗段中進行(見圖1)。試驗段為直徑100mm的有機玻璃圓管，可方便CCD拍攝，出口處裝有皮托管，通過差壓計測量來流速度[6-7]。為了降低氣流湍流度，試驗段底部裝有阻尼網(wǎng)，同時可防止實驗物落入風(fēng)洞中，在出口處裝有隔網(wǎng)防止實驗物飛出試驗段。風(fēng)洞風(fēng)速通過可控硅-交流電機調(diào)節(jié)，調(diào)節(jié)范圍0～25m/s。

1.2實驗?zāi)Ｐ?/p>

為了研究下落物體的下落特征，有針對性地選取了幾類典型形狀的物體進行實驗測量，包括：圓球、立方體、長方體、短圓柱體和圓盤。圓球為乒乓球，直徑40mm，重2.7g，最大截面積為1256.6mm2。其余模型由輕質(zhì)木材做成，模型參數(shù)如表1所示。圖2是試驗?zāi)Ｐ偷恼掌?/p>

圖1 風(fēng)洞立式試驗段

圓球、立方體、長方體、短圓柱體及圓盤

表1模型參數(shù)及平均下落速度測量結(jié)果

Table1Characterofmodelstesedinthewindtunnel

實驗?zāi)Ｐ湍Ｐ蛥?shù)及實驗結(jié)果準懸浮狀態(tài)平均風(fēng)速/(m·s?1)幾何尺寸體積/mm3質(zhì)量/g乒乓球8．5直徑40mm335102．7立方體7．8邊長28．94mm270002．7長方體約7．3邊長37．8mm、25mm、37．8mm357213．6短圓柱約6．8直徑40mm、高20mm251322．5圓盤約4．0直徑43．7mm、高10mm149981．5

1.3物體自由下落實驗觀察

首先將實驗?zāi)Ｐ椭糜陲L(fēng)洞垂直試驗段入口處，然后緩慢提高風(fēng)速，當(dāng)達到一定風(fēng)速后模型將升起進入試驗段并出現(xiàn)某種懸浮狀態(tài)。圖3為各實驗?zāi)Ｐ蛻腋≡谠囼灦沃械恼掌?/p>

從左至右：圓球、立方體、長方體、圓柱體及圓盤

實驗觀察到的這5種模型的懸浮狀態(tài)可以分為3類：

第1類，在一定風(fēng)速下，出現(xiàn)較穩(wěn)定的懸浮狀態(tài)。圓球與正方體模型屬于這一類，前者中心對稱，后者有三個正交的對稱軸。在一定風(fēng)速的氣流中，它們能較穩(wěn)定地懸浮于立式風(fēng)洞的試驗段中。圖4顯示的是風(fēng)洞試驗段中懸浮的正方體。觀察發(fā)現(xiàn)，懸浮的圓球和正方體模型會有一些方向和大小都是隨機的轉(zhuǎn)動和擺動。但相對來說，球模型轉(zhuǎn)動和擺動的幅度和速度要遠小于正方體模型。此外還觀察到處于懸浮狀態(tài)的正方體模型沒有出現(xiàn)一個頻度占優(yōu)的方位姿態(tài)，也就是說方位姿態(tài)是隨機的。由此不難想象，一個高空下落的正方物體可能以任一個方位姿態(tài)著地。

圖4 風(fēng)洞試驗段中懸浮的正方體

第2類，懸浮狀態(tài)較不穩(wěn)定。反映在風(fēng)洞實驗過程中是調(diào)整風(fēng)速使模型懸浮比較困難。在某一個風(fēng)速范圍內(nèi)，模型都會出現(xiàn)懸浮狀態(tài)，但卻維持時間不長就從試驗段中消失，或掉到試驗段入口，或沖到試驗段出口。本次試驗中，長方柱模型和短圓柱模型屬于這一類。實驗觀察發(fā)現(xiàn)在懸浮狀態(tài)下，這兩種模型自身都在做隨機的轉(zhuǎn)動和擺動，其幅度和速度都大于第1類的對稱性強的球形和正方形模型。而它們還有一個特點，就是處于懸浮狀態(tài)時，主要都是處于氣動阻力較大而不是較小的方位姿態(tài)，如長方體模型的迎風(fēng)面是面積大的矩形平面，而短圓柱則是圓平面。由此也可以推想，從高空落下的這類形狀的物體，以面積大的一面著地的概率較大。

第3類，懸浮狀態(tài)很不穩(wěn)定。圓盤模型就是屬于這一類。在風(fēng)洞立式試驗段中，在某一范圍的風(fēng)速下，都可能出現(xiàn)圓盤的懸浮狀態(tài)，但這種狀態(tài)保持的時間大都很短，如幾秒鐘，極少情況下也會較長，如10s以上。處于這種懸浮狀態(tài)時，模型會隨機地快速轉(zhuǎn)動和擺動，同時大部分時間是圓盤面處于大迎角而不是小迎角的狀態(tài)。可以理解為：因為圓盤在小迎角狀態(tài)是氣動不穩(wěn)定的，即氣動力矩會使其迎角增大進入大迎角的隨機擺動和轉(zhuǎn)動狀態(tài)。在自然界中，寬大的落葉是搖搖晃晃飄落下來的而不是以很小的迎角快速砸向地面也是同樣的機理。

1.4下落平均速度的測量

在風(fēng)洞立式試驗段中，實驗?zāi)Ｐ瓦M入接近懸浮的狀態(tài)時，物體受到的重力與阻力基本平衡，此時的風(fēng)速被認為是模型下落達到比較穩(wěn)定速度時的平均速度。這個速度可以通過風(fēng)洞立式試驗段出口的皮托管風(fēng)速計測得，結(jié)果如表1所示。實際上，對于對稱性好的圓球和正方物體模型，這個速度和它們自由下落的最大速度對應(yīng)，但對于對稱性不好的長方體、短圓柱特別是圓盤，表1的風(fēng)速只是對應(yīng)一種可能的懸停風(fēng)速，因而也是一種可能出現(xiàn)的(不是唯一的)下降的最大風(fēng)速。

1.5下落瞬時速度的概率分布

為了得到模型下落瞬時速度，詳細分析模型懸浮狀態(tài)出現(xiàn)時間段內(nèi)的所有圖像(每幀圖像拍攝時間間隔為1/15s)，通過獲取模型中心位置的坐標得到1/15s的平均速度，并近似認為這一速度值是瞬時速度。由于此速度為物體懸浮狀態(tài)下相對風(fēng)洞的運動速度，需疊加與風(fēng)洞氣流相對的速度得到物體的實際下落瞬時速度，最終將所有瞬時速度繪成箱式圖(見圖5)。

圖5 瞬時下落速度分布箱式圖

圖5反映出5種模型下落速度的分布規(guī)律，針對每個形狀的分布圖，上下兩點為最大值和最小值，中間的兩個長方形上沿代表速度概率分布為95％的速度值，下沿為5％，中間的水平直線代表分布函數(shù)的中值，而長方形中間的小正方形代表算數(shù)平均值。從瞬時速度值的平均值來看，球形、正方形、長方體、短圓柱體和圓盤的平均下落速度依次降低，這是由各模型的氣動外形不同引起的。除了平均下落速度，瞬時速度的概率分布可以更好地反映模型下落過程中的細節(jié)信息。從瞬時速度概率分布形式來看，球形和正方形由于其對稱性較好，同時在實驗中觀察到模型出現(xiàn)隨機性旋轉(zhuǎn)，其迎風(fēng)面積和風(fēng)阻系數(shù)等氣動參數(shù)隨時間的變化不大，球形和正方形模型下降速度概率分布函數(shù)的范圍較窄，最大值和最小值之差也較小。而由于正方形較球形來講不同姿態(tài)下的氣動特征稍有不同，這也使其速度概率分布函數(shù)上最大與最小值之間的差異較球形稍大。長方體和短圓柱的對稱性稍差，不同姿態(tài)角下的迎風(fēng)面積和阻力系數(shù)差距很大，從其速度概率分布函數(shù)來看，其寬度較對稱性較好的球形和正方形有較大增加，說明模型在下落過程中速度的震蕩較大，即模型以不同的姿態(tài)角和不同的速度下落，速度最大最小值之差可以達到1m/s。對于圓盤，最大最小值之間的差別很大，達到2.5m/s，這主要是由于其對稱性較差，圓盤不同姿態(tài)時的迎風(fēng)面積差距較大，可以帶來下落速度的較大差別，圓盤下落速度最小值出現(xiàn)在迎風(fēng)面積最大的姿態(tài)，而最大值則出現(xiàn)在圓盤迎角為0的姿態(tài)。由其下落速度概率分布函數(shù)的特征可知，兩種極限狀態(tài)并不是主要的狀態(tài)，也就是說圓盤不可能保持以這兩種姿態(tài)下落，實際下落是介于最大和最小值之間的狀態(tài)，由之前的實驗觀察可知，圓盤主要是以快速旋轉(zhuǎn)和擺動的狀態(tài)下落，并且是大迎角占優(yōu)，小迎角屬于不穩(wěn)定狀態(tài)。

2 數(shù)值模擬

為了更加詳細研究物體在下落過程中的運動行為，采用數(shù)值模擬手段對實驗結(jié)果進行補充研究。主要應(yīng)用FLUENT軟件中的k-ε模型[8-9]求解不可壓Navier-Stokes方程組得到速度和壓力場分布，對物體下落過程中物體周圍的流體流動狀態(tài)進行研究，同時得到物體在空氣中下落姿態(tài)特征、運動軌跡及受力情況。為了簡化問題，采用二維數(shù)值計算，選擇正方形作為研究對象，采用了兩種數(shù)值模擬方法。

第1種模擬方法模擬了固定的正方形模型在不同來流方向時的繞流流場，計算域尺寸2m×2m，正方形尺寸為0.2m×0.2m，位于計算域中心。正方形周圍網(wǎng)格劃分采用正方形網(wǎng)格，較遠處為三角形網(wǎng)格，網(wǎng)格均采用固定網(wǎng)格，計算正方形在不同姿態(tài)下受到的力和力矩。取正方形一個邊(見圖7左上圖，最下面的邊)的垂直線作為基準，與來流夾角為0°～45°，間隔5°，共10個角度。圖7給出了其中4個角度(0、15°、30°及45°)。網(wǎng)格節(jié)點數(shù)為3000左右。

第2種數(shù)值模擬采用了更大的計算域，尺寸為20m×20m，計算中采用了動態(tài)網(wǎng)格技術(shù)，用于模擬正方形的運動過程。相對較大的計算域使得物體有足夠的運動空間，對于下落速度不穩(wěn)定的情況，可以得到更長時間段內(nèi)的運動軌跡(因為運動網(wǎng)格技術(shù)不能處理物體運動到計算域之外的情況，一旦正方形移動出計算域，計算將停止)。正方形周圍的網(wǎng)格劃分被加密，可以提高計算精度同時減少遠處梯度較小處的計算量，總網(wǎng)格節(jié)點數(shù)為10000左右。

2.1動態(tài)網(wǎng)格技術(shù)

如上所述，為了計算正方形的運動過程，采用了動態(tài)網(wǎng)格技術(shù)[10]，圖6為初始狀態(tài)和計算過程中的網(wǎng)格分布的局部圖(t=0、5、10及15s四個時刻)。在初始狀態(tài)時，正方形位于計算域中心，通過用戶自定義文件給出正方形的質(zhì)量和扭轉(zhuǎn)慣量，正方形在重力、流體壓力及粘性應(yīng)力作用下，發(fā)生運動。每隔一個計算時間步長，通過六自由度受力分析可以得到正方形在X方向和Y方向上的速度以及Z方向上的轉(zhuǎn)動角速度。根據(jù)速度參數(shù)可以得到正方形下一時刻的位置，網(wǎng)格會被自動重新劃分，并進行下一步計算。

2.2邊界條件及求解過程

計算域入口(底部)為速度邊界條件，出口為壓力出口邊界條件，兩側(cè)為無滑移壁面邊界條件，正方形模型四邊為固壁邊界條件。由于作了二維簡化，來流風(fēng)速不能以實驗得到的風(fēng)速作為參考，通過調(diào)整入口邊界條件，同時監(jiān)測物體阻力系數(shù)和y方向速度的變化，當(dāng)y方向的加速度接近0時認為正方形在此速度下會進入懸浮狀態(tài)，由此確定入口速度邊界條件的速度大小為20.5m/s。

對于第1種模擬方法，采用定常求解方法，模擬各種姿態(tài)下正方形所受的力和扭矩。對于第2種方法，為了模擬正方形在計算域中的運動過程，采用了時間非定常求解方法，計算時間間隔為0.5ms。

圖6 動態(tài)網(wǎng)格

圖8為正方形受到的阻力系數(shù)和扭矩系數(shù)隨角度的變化規(guī)律。阻力系數(shù)在0°的時候最小，隨著角度的增加，阻力系數(shù)會隨之增加。當(dāng)角度為45°時，阻力系數(shù)達到極值。由于正方形的對稱性，超過45°之后的受力情況會以45°為對稱軸呈對稱分布。當(dāng)角度為90°時，受力情況與0°相同。

圖8同時給出扭矩系數(shù)的分布規(guī)律，當(dāng)正方形逆時針旋轉(zhuǎn)一定角度之后(小于45°)，物體扭矩系數(shù)為負值，即受到一個使之順時針旋轉(zhuǎn)的扭矩，扭矩絕對值大小隨著角度的增加而增加。當(dāng)角度達到15°左右時，達到最大值。當(dāng)角度進一步增加之后，扭矩的絕對值下降，當(dāng)角度為45°時，扭矩恢復(fù)為接近0。超過45°之后，由于正方形的對稱性，受到的力矩方向發(fā)生變化，絕對值以45°為對稱軸呈對稱分布，即一旦角度超過45°，正方形將受到一個使之偏轉(zhuǎn)角度增加的力矩。

圖7 速度場分布(以最大速度歸一處理)

圖8 阻力系數(shù)及扭矩系數(shù)圖

2.3運動軌跡及姿態(tài)數(shù)值模擬

除了不同姿態(tài)角下正方形周圍流場的計算模擬以外，采用動態(tài)網(wǎng)格技術(shù)對處于20m×20m計算域中的正方形的運動規(guī)律進行了數(shù)值模擬研究。正方形的初始位置在計算域中心，初始角度為0°。圖9顯示了正方形在50s的計算時間段內(nèi)角度隨時間的變化規(guī)律。當(dāng)T=0時，雖然根據(jù)穩(wěn)態(tài)計算結(jié)果(見圖7)正方形的受力為對稱受力，受到的扭矩也接近零，即正方形能夠保持這個姿態(tài)，但是在瞬態(tài)流場計算過程中，氣流的不穩(wěn)定性可以使正方形發(fā)生微小的角度偏轉(zhuǎn)。

根據(jù)前文對扭矩的計算可知，當(dāng)正方形姿態(tài)發(fā)生變化后，會受到一個反向的扭矩使之回到初始狀態(tài)(角度為0的狀態(tài))。但是由于慣性的作用，當(dāng)正方形姿態(tài)角為0時，仍然存在一個角速度使之越過角度為0這個姿態(tài)，隨之受到反向力矩使角速度絕對值降低，直至降為0，這時正方形的姿態(tài)角達到極值；而后角速度改變方向，正方形的姿態(tài)角再次向角度為0的狀態(tài)變化。從模擬結(jié)果(見圖9)可以看出正方形姿態(tài)度隨時間呈現(xiàn)出上述震蕩變化過程，其振幅在10s之內(nèi)會出現(xiàn)逐漸增加的現(xiàn)象，但是在接下來的過程中這一幅值會穩(wěn)定在27°～32°之間。

圖10為正方形在計算域中的軌跡，圖11為計算過程中t=10、15、20和25s 4個時刻流場中的速度分布圖，由最大速度25m/s進行歸一化處理。正方形的初始姿態(tài)為0°，阻力系數(shù)最小，這時正方形受到的重力要稍大于阻力，因此存在向下的速度，從圖10得到的正方形運動軌跡可以看出這是一個下落過程。隨著正方形姿態(tài)角的變化，尤其是當(dāng)達到30°左右時，由于阻力系數(shù)的升高，阻力與重力大體平衡，從軌跡圖上可以看到Y(jié)值的振蕩規(guī)律，即正方形進入懸浮狀態(tài)。

圖9 正方形姿態(tài)角度隨時間的變化

圖10 正方形的運動軌跡

圖11 速度分布圖

3 討論與結(jié)論

利用小型立式風(fēng)洞對五種典型形體的剛性模型進行了其在上升氣流中懸浮特性的實驗，主要目的之一是研究物體在近地大氣中自由下落達到重力與空氣阻力平衡時下降速度的特性。實驗觀察、測量和分析表明，空間對稱性好的物體具有穩(wěn)定的懸停特性，表明這類物體有穩(wěn)定的自由下降最大速度；而空間對稱性不好的物體，則在上升氣流中呈現(xiàn)不穩(wěn)定的懸停狀態(tài)，表明這類物體具有不確定的最終下降速度。換言之，從高空(這里假定是幾百米的近地高空)落物來看，這類具有不穩(wěn)定懸停狀態(tài)的物體的最終落地速度是不確定的。它有一個理論極限值，即始終以最小阻力姿態(tài)下落時達到的最大速度。但實際的自由落體運動中這個速度是很難達到的，因為最小阻力的方位姿態(tài)多是氣動不穩(wěn)定的，在自由下降的過程中任何擾動都將發(fā)散性地改變其空間姿態(tài)，增大了阻力，降低了速度。

用小型的立式風(fēng)洞和小模型做懸浮實驗?zāi)軌蚨ㄐ院湍撤N程度定量地反映不同形體的物體在上升氣流中的懸浮特性和由此推斷出的自由落體在重力和氣動阻力達到平衡時的下降速度特性。但由于風(fēng)洞截面積和模型截面積之比為6，風(fēng)洞壁影響了懸浮狀態(tài)的物體產(chǎn)生的橫向運動，此外試驗段高度也不夠，不能更好地反映懸浮不穩(wěn)定的運動歷程。所以實驗的一些定量結(jié)果只有參考意義。

二維正方形物體周圍流場以及正方形的阻力系數(shù)、扭矩和運動軌跡等數(shù)值試驗結(jié)果表明：正方形不存在姿態(tài)角不變的平衡位置，在氣流小擾動的作用下，正方形發(fā)生微小的角度偏轉(zhuǎn)，在與角速度相反的扭矩和慣性共同作用下，正方形物體姿態(tài)角出現(xiàn)正負振蕩現(xiàn)象，振蕩幅值最終穩(wěn)定在27～32°。

致謝:感謝北京大學(xué)工學(xué)院湍流與復(fù)雜系統(tǒng)國家重點實驗室提供的實驗設(shè)備。

參考文獻:

[1]Graf W.Hydraulics of sediment transport[M].New York: McGraw-Hill,1971.

[2]周光炯,嚴宗毅,等.流體力學(xué)[M].北京: 高等教育出版社,1993.

[3]Joseph D D,Liu Y J.Orientation of long bodies falling in a viscoelastic liquid[J].J Rheol,1993,37: 1-22.

[4]Liu Y J,Nelson J,Feng J,et al.Anomalous rolling of spheres down an inclined plane[J].Non-Newtonian Fluid Mechanics,1993,50(2/3): 305-330.

[5]Mollinger A M,Cornelissen E C,van den Brule B H A A.An unexpected phenomenon observed in particle settling: oscillating falling spheres[J].J Non-Newtonian Fluid Mech,1999,86: 389-393.

[6]Merzkirch,Wolfgang.Flow visualization[M].London: Academic Press,1987.

[7]Pankhurst R C,Holder D W.Wind-tunnel technique[M].London: Sir Isaac Pitman & Sons,1952.

[8]Hinze J O.Turbulence[M].New York: McGraw-Hill,1975.

[9]Launder B E,Spalding D B.Lectures in mathematical models of turbulence[M].London: Academic Press,1972.

[10]Snyder D O,Koutsavdis E K,Anttonen J S R.Transonic store separation using unstructured cfd with dynamic meshing[R].AIAA-2003-3913,2003.

作者簡介:

李軼明(1975-)，男，吉林省吉林市人，講師。研究方向：實驗流體力學(xué)及石油工程井筒多相流。通信地址：中國石油大學(xué)(北京)石油工程學(xué)院油氣井工程系(102249)。E-mail: ymli@cup.edu.cn