張 嶠，劉冬雨，羅 超，魏傳鋒

（中國空間技術(shù)研究院 載人航天總體部，北京 100094）

0 引言

“天宮二號”空間實驗室的實驗艙內(nèi)搭載了機(jī)械臂操作終端系統(tǒng)[1]。計劃在空間實驗室與載人飛船組合體飛行期間，航天員乘組與機(jī)械臂操作終端協(xié)同完成首次人機(jī)協(xié)同在軌維修技術(shù)試驗，為航天員與空間機(jī)器人的人機(jī)協(xié)同作業(yè)積累在軌經(jīng)驗。

在正式的在軌維修試驗之前，工程技術(shù)人員先期設(shè)計了機(jī)械臂操作終端在軌抓取小球的驗證試驗。在軌微重力環(huán)境下，密封艙內(nèi)不存在自然對流，載人環(huán)境的建立需通過艙內(nèi)強(qiáng)制通風(fēng)來保證，這將導(dǎo)致漂浮在艙內(nèi)的小球隨流場運(yùn)動，對機(jī)械臂操作終端抓取小球任務(wù)構(gòu)成干擾項。

本文建立了密封艙內(nèi)漂浮小球運(yùn)動的物理模型與數(shù)學(xué)模型，采用動網(wǎng)格方法，對小球的運(yùn)動速度和位移進(jìn)行數(shù)值模擬，得到典型工況下小球的運(yùn)動規(guī)律，并由此確立小球的設(shè)計參數(shù)范圍，作為機(jī)械臂操作終端系統(tǒng)試驗方案設(shè)計的依據(jù)。

1 計算模型

1.1 物性參數(shù)

實驗艙內(nèi)壓力體制基本均衡（標(biāo)稱1個標(biāo)準(zhǔn)大氣壓），艙內(nèi)通風(fēng)氣流速度不大于0.5 m/s，可將氣流簡化為不可壓縮流體。計算所用物性參數(shù)如表1所示。

表1 計算所用物性參數(shù)Table 1 Physical properties of the sealed cabin for calculation

1.2 邊界條件

通過分析，機(jī)械臂操作終端的機(jī)械手可對標(biāo)準(zhǔn)網(wǎng)球尺寸（φ=66 mm）的小球進(jìn)行抓取，考慮到小球幾何構(gòu)型的軸對稱性，將計算域簡化為圖1所示區(qū)域。

圖1 計算區(qū)域示意Fig.1 Computational model and domain

由于氣流已簡化為不可壓縮流體，所以可直接將左側(cè)流場入口定義為速度入口，右側(cè)定義為流場出口。進(jìn)一步將自由空間簡化為半徑0.3 m（約為小球半徑的10倍）的圓柱體計算域，邊界為無滑移壁面。

1.3 控制方程

不考慮化學(xué)反應(yīng)，對于單組分氣相工質(zhì)，流場控制方程如下[2-3]：

式中：ρ為流體密度；ui、uj分別為流體在x、y方向的速度矢量；p為流體壓強(qiáng)；為雷諾應(yīng)力；σij為黏性剪切應(yīng)力，

其中，μ為分子黏性系數(shù)，δij為 Kronecker Delta，當(dāng)i=j時δij=1，當(dāng)i≠j時δij=0。

使用有限體積法將控制方程組離散為網(wǎng)格單元上的線性代數(shù)方程[4]。當(dāng)計算對空間與時間分辨率要求不高時，直接采用雷諾時均模型（RANS），假定湍流中的流場變量由一個時均量和一個脈動量組成，引入Boussinesq假設(shè)，認(rèn)為湍流雷諾應(yīng)力與應(yīng)變成正比，則湍流模擬可歸結(jié)為對雷諾應(yīng)力與應(yīng)變之間的比例系數(shù)（即湍流黏性系數(shù)）的計算[3-5]。本文采用帶旋流修正的兩方程k-ε湍流模型對控制方程進(jìn)行封閉，并利用增強(qiáng)壁面函數(shù)法處理流動中的近壁區(qū)域[6-7]。

1.4 動網(wǎng)格模型

使用動網(wǎng)格模型模擬流場形狀由于邊界運(yùn)動而隨時間改變的問題。網(wǎng)格的更新過程根據(jù)每個迭代步中邊界的變化情況自動完成[5,8]。系統(tǒng)考慮式(1)～式(3)，在任意一個控制體中，廣義標(biāo)量Φ的積分守恒方程為

式中：u為流體速度矢量；ug為移動網(wǎng)格的網(wǎng)格速度；Γ為擴(kuò)散系數(shù)；SΦ為源項；?V代表控制體V的邊界。

方程中的時間導(dǎo)數(shù)項，可以用一階后向差分格式表示為

其中n和n+1代表不同的時間層。

其中dV/dt是控制體的時間導(dǎo)數(shù)。為了滿足網(wǎng)格守恒定律，控制體的時間導(dǎo)數(shù)可以表示為

其中：nf是控制體積的面網(wǎng)格數(shù)；Aj為面j的面積矢量；

其中δVj為控制體積中面j在時間間隔Δt內(nèi)掃過的空間體積。

微重力環(huán)境下，忽略小球所受重力，來流對小球的作用力驅(qū)使小球隨動，作用在小球壁面上的力Fi包含壓強(qiáng)效應(yīng)和黏性剪切力2部分，即

式中Fv即為黏性力，由近壁函數(shù)求得，i代表空間矢量分量。由于來流速度較低，所以壓強(qiáng)效應(yīng)與黏性剪切力的比例相當(dāng)，黏性力不可忽略。

由小球表面的作用力，可推知當(dāng)前時刻小球在流場中的加速度，則下一時刻小球的運(yùn)動速度及位移可表示為：

將小球視為剛體，通過編寫用戶自定義函數(shù)（UDF），調(diào)用Define_CG_Motion宏，對小球的運(yùn)動方式進(jìn)行描述[4-5]。

與方形結(jié)構(gòu)相比，小球構(gòu)型不甚規(guī)則，為了避免小球移動過程中出現(xiàn)負(fù)體積網(wǎng)格，采用彈簧順滑動網(wǎng)格模型實時更新流場，使用非結(jié)構(gòu)網(wǎng)格對整個計算域進(jìn)行離散，初始網(wǎng)格數(shù)量為27 200。

1.5 離散格式

對動量方程使用二階迎風(fēng)格式進(jìn)行離散，對湍流方程使用一階迎風(fēng)格式進(jìn)行離散。動網(wǎng)格計算中采用一階迎風(fēng)格式對時間項（式(6)）進(jìn)行離散，步長為0.005 s，共計算20 s。非穩(wěn)態(tài)流場壓力-速度的耦合求解采用PISO算法[9-10]。

2 計算結(jié)果與分析

2.1 流場特性分析

將小球質(zhì)量m=58 g（標(biāo)準(zhǔn)網(wǎng)球重量）、來流速度v=0.1 m/s的工況定義為基準(zhǔn)工況，對基準(zhǔn)工況進(jìn)行計算，不同時刻下的計算域網(wǎng)格如圖2所示?？梢钥闯?，計算網(wǎng)格隨著時間的推移而重構(gòu)：網(wǎng)格邊界移動前由“彈簧”組成的系統(tǒng)處于平衡狀態(tài)；在網(wǎng)格邊界節(jié)點發(fā)生位移后，節(jié)點位移形成的力破壞了“彈簧系統(tǒng)”原有的平衡，系統(tǒng)經(jīng)過調(diào)整將達(dá)到新的平衡，即在新的位置上重新獲得受力平衡，據(jù)此經(jīng)過迭代計算，得到使各節(jié)點上的合力等于0的、新的網(wǎng)格節(jié)點位置。當(dāng)運(yùn)動邊界的位移遠(yuǎn)遠(yuǎn)大于網(wǎng)格尺寸時，為了避免因網(wǎng)格畸變過大導(dǎo)致計算不收斂，須將局部網(wǎng)格重新劃分或合并，用新的網(wǎng)格代替原來的網(wǎng)格。

圖2 基準(zhǔn)工況中網(wǎng)格隨時間的變化Fig.2 Sequence of gird variations in the basic calculation case

與圖2對應(yīng)的速度幅值云圖如圖3所示。在t=0時刻，小球處于準(zhǔn)穩(wěn)態(tài)，其左右兩側(cè)的受力均勻，沒有運(yùn)動速度；隨著自由來流撞擊小球，小球被流場驅(qū)動。通過圖3(b)、3(c)可以看出，由于自由來流的黏性特征，小球背風(fēng)面將出現(xiàn)由逆壓梯度帶來的不斷脫落的旋渦區(qū)。圖4給出了圖3對應(yīng)的流線圖，展現(xiàn)了小球背風(fēng)區(qū)的旋渦形狀。脫落的旋渦將會導(dǎo)致背風(fēng)面壓力不穩(wěn)定，因此小球表面受力將會呈現(xiàn)周期特性。

圖3 基準(zhǔn)工況中不同時刻速度幅值云圖Fig.3 Velocity contours in the basic case

圖4 基準(zhǔn)工況t=10 s時刻背風(fēng)區(qū)流線圖Fig.4 Streamlines in the basic case at t=10 s

2.2 小球運(yùn)動特性分析

2.2.1 計算工況

為系統(tǒng)研究小球在自由來流流場中的運(yùn)動規(guī)律，下面分別以小球質(zhì)量及來流速度為變化參數(shù)，對表2所示的10組工況進(jìn)行仿真分析。

表2 計算工況設(shè)置Table 2 Parameter settings for the calculation cases

2.2.2 小球運(yùn)動特性

小球速度隨時間變化規(guī)律與受力變化規(guī)律一致，在受力轉(zhuǎn)入振蕩狀態(tài)之前，速度曲線斜率較大；當(dāng)受力轉(zhuǎn)入振蕩狀態(tài)之后，速度斜率下降，基本呈現(xiàn)勻加速特性。將速度曲線積分可得到位移曲線。由于小球運(yùn)動速度不斷加快，所以位移隨時間推進(jìn)而快速增大。

對工況1～工況5進(jìn)行系列計算，小球受力、運(yùn)動速度和運(yùn)動位移隨時間變化規(guī)律分別如圖5、圖6和圖7所示?？梢钥闯觯鹤饔迷谛∏虮砻娴牧Υ嬖谝粋€明顯的啟動峰，來流速度越大，啟動峰越明顯；當(dāng)小球剛剛啟動時，作用力較大，亦即加速度較大；隨著小球在流場中被“帶著”移動，小球所受作用力呈逐漸下降的趨勢，這是由于小球自身慣性所產(chǎn)生的滯后效應(yīng)而導(dǎo)致的。還可以看出，隨著來流速度的增加，小球表面作用力迅速轉(zhuǎn)捩到振蕩狀態(tài)。結(jié)合2.1節(jié)的流場特性分析可知，振蕩受力正是由于小球背風(fēng)區(qū)旋渦脫落導(dǎo)致的。

對工況6～工況10進(jìn)行系列計算，小球受力、運(yùn)動速度和運(yùn)動位移隨時間變化分別如圖8、圖9和圖10所示。經(jīng)過對比可以看出，500 g小球的受力特性與58 g小球的基本相同，均呈現(xiàn)出啟動峰-穩(wěn)態(tài)-振蕩狀態(tài)3個分部區(qū)域；不同之處在于500 g小球的慣性比58 g小球大，因此受力曲線轉(zhuǎn)入振蕩狀態(tài)的時間點后移，小球運(yùn)動滯后特性更加顯著，導(dǎo)致振蕩周期變長。

圖5 58 g小球受力隨時間變化的計算結(jié)果Fig.5 The force on the floating ball of 58 g

圖6 58 g小球速度隨時間變化的計算結(jié)果Fig.6 Velocity characteristics of the floating ball of 58 g

圖7 58 g小球位移隨時間變化的計算結(jié)果Fig.7 Displacement characteristics of the floating ball of 58 g

圖8 500 g小球受力隨時間變化的計算結(jié)果Fig.8 Calculated force on the floating ball of 500 g against time

圖9 500 g小球速度隨時間變化的計算結(jié)果Fig.9 Calculated velocity of the floating ball of 500 g against time

圖10 500 g小球位移隨時間變化的計算結(jié)果Fig.10 Calculated displacement of the floating ball of 500 g against time

2.3 分析結(jié)果

對2.2.2節(jié)中小球運(yùn)動的關(guān)鍵數(shù)據(jù)進(jìn)行提取，形成表3。對比兩種小球的運(yùn)動速度及位移變化曲線得出，在運(yùn)動末期500 g小球的速度及位移均約為58 g小球的1/7?？梢钥闯?，當(dāng)最大來流速度為0.5 m/s時，標(biāo)準(zhǔn)網(wǎng)球重量的小球在10 s時的運(yùn)動位移將超過400 mm，此種情況下機(jī)械臂操作終端捕獲小球的成功率將大大降低。綜合考慮機(jī)械臂操作終端的運(yùn)動速度、機(jī)械手抓取速度及視覺相機(jī)視場特性， 用于在軌試驗的小球的質(zhì)量應(yīng)不小于500 g，抓取小球試驗的有效時長應(yīng)不長于20 s。

表3 小球運(yùn)動關(guān)鍵數(shù)據(jù)Table 3 Key parameters of ball motion

3 在軌試驗

2016年10月—11月，“神舟十一號”航天員乘組駐留“天宮二號”空間實驗室期間，機(jī)械臂操作終端完成了2次抓取500 g漂浮小球試驗任務(wù)，均一次性成功。圖11展示了漂浮小球從靜止到運(yùn)動10、20 s的圖像。將小球尺寸（φ=66 mm）作為標(biāo)準(zhǔn)刻度進(jìn)行對比可以看出，小球在10 s時刻的運(yùn)動距離為54.9 mm，20 s時刻的運(yùn)動距離為185.5 mm，與圖10、表3中500 g小球在0.5 m/s的艙內(nèi)來流速度下的仿真預(yù)示值吻合程度較高，初步驗證了本文計算結(jié)果的正確性。

圖11 小球在軌運(yùn)動圖像Fig.11 Movement of floating ball in-orbit

圖12給出機(jī)械臂操作終端在軌抓取漂浮小球的始、末畫面，試驗取得成功。首次人機(jī)協(xié)同在軌維修技術(shù)試驗，為航天員與空間機(jī)器人的協(xié)同作業(yè)積累了在軌經(jīng)驗，為提高空間站機(jī)械臂精細(xì)操作能力和控制魯棒性提供了在軌基礎(chǔ)方法和參數(shù)。

圖12 機(jī)械臂操作終端抓取漂浮小球試驗圖像Fig.12 Pictures of floating ball caught by space manipulator terminal

4 結(jié)束語

本文采用CFD動網(wǎng)格方法，對密封艙內(nèi)漂浮小球在軌隨來流的運(yùn)動規(guī)律進(jìn)行了數(shù)值模擬，得到了小球的受力、運(yùn)動速度和位移隨時間變化的規(guī)律。研究結(jié)果表明：小球在載人密封艙內(nèi)通風(fēng)流場中的運(yùn)動呈現(xiàn)非線性特性；質(zhì)量較大的小球在流場作用下轉(zhuǎn)入振蕩狀態(tài)的時間點后移，小球運(yùn)動滯后，振蕩周期變長。

通過研究，工程技術(shù)人員得到了小球的設(shè)計參數(shù)范圍，作為機(jī)械臂操作終端系統(tǒng)試驗方案設(shè)計的依據(jù)?！吧裰凼惶枴焙教靻T乘組駐留“天宮二號”空間實驗室期間，機(jī)械臂操作終端完成了抓取漂浮小球的試驗任務(wù)，初步驗證了本文計算結(jié)果的正確性。