，，，，，

1. 航天工程大學(xué) 激光推進及其應(yīng)用國家重點實驗室，北京 101416 2. 西安衛(wèi)星測控中心，西安 710043

基于多顆微納衛(wèi)星的編隊、集群、星座，可以承擔(dān)更為復(fù)雜的空間任務(wù)[1]。星載微推進器是微納衛(wèi)星實現(xiàn)精確姿控、軌控、快速機動的基礎(chǔ)。目前，國際上已研制出多種新概念微推進器，如膠體微推進器[2]、激光燒蝕微推進器[3]、脈沖等離子體微推進器[4]、離子推進器等[5]，產(chǎn)生的推力通常在微牛至毫牛量級，甚至小到亞微牛量級[6]，給微推力的測量帶來一定的難度。

目前，推力測量基本都是基于力的動力效果，將推力轉(zhuǎn)換成測量臺架的振動幅值或轉(zhuǎn)動位移[7]。按照推進器是否與測量裝置固連，可將測量方法分為直接法與間接法[8]。直接法中，推進器固連在測量臺架的執(zhí)行元件上，推力直接轉(zhuǎn)換為測量臺架的振動幅值或轉(zhuǎn)動位移，典型裝置有扭擺結(jié)構(gòu)[9]、單擺結(jié)構(gòu)[10]、天平結(jié)構(gòu)[11]。直接法精度較高，但由于臺架承重推進器，導(dǎo)致固有頻率較小，因而響應(yīng)速度較慢，難以完成推進器設(shè)計者十分關(guān)注的動態(tài)推力的測量。相比之下，采用間接法測量時，推進器與測量臺架物理隔離，測量臺架固有頻率較高，能夠較快響應(yīng)推力變化，實現(xiàn)動態(tài)推力測量，同時還可避免供應(yīng)管路及測控電纜干擾，進行裝星后測量，典型裝置有吊擺結(jié)構(gòu)[12-13]、彈性盤結(jié)構(gòu)[14]、懸臂梁結(jié)構(gòu)[15-17]。

間接法測量結(jié)構(gòu)中，懸臂梁結(jié)構(gòu)簡單，測量時既作為推進器噴射羽流的承接元件，又作為形變元件，通過改變尺寸即可調(diào)整工作頻帶，是應(yīng)用最廣泛的力學(xué)結(jié)構(gòu)。目前應(yīng)用懸臂梁結(jié)構(gòu)時，主要基于靜力學(xué)模型，依據(jù)穩(wěn)態(tài)彎曲量與外力大小的線性關(guān)系對推進器穩(wěn)態(tài)推力進行測量[15-17]。而當(dāng)待測推力動態(tài)變化時，懸臂梁振幅較大，難以獲得穩(wěn)態(tài)彎曲量，無法解算推力，需要根據(jù)懸臂梁動力學(xué)模型，得到穩(wěn)態(tài)位移的時變值，進而確定推力測量方法。同時，由于噴射羽流的發(fā)散及反射氣流對來流干擾等因素影響，懸臂梁測量得到的沖擊力與實際推力存在一定誤差。文獻[17]進行了相應(yīng)的仿真分析，但受供氣軟管的影響，試驗中只能測量入口壓力0.35 MPa以下推力，且精度受限。需要采用精度較高的試驗方法，對沖擊力與推力關(guān)系進行研究，以量化測量誤差。

本文在懸臂梁動力學(xué)模型的基礎(chǔ)上，確定了系統(tǒng)傳遞函數(shù)，根據(jù)系統(tǒng)響應(yīng)特征，確定了動態(tài)推力測量方法。同時，通過試驗驗證了測量方法的正確性，并通過與扭擺測量結(jié)果進行對比，分析了測量誤差。

1 系統(tǒng)傳遞函數(shù)的建立

單位脈沖函數(shù)拉普拉斯變換為1，求得單位脈沖力下的系統(tǒng)響應(yīng)，做拉普拉斯變換即可求得系統(tǒng)傳遞函數(shù)。

如圖1所示，選擇懸臂梁未變形時的軸向，即各截面形心的連線作為x軸，選擇對稱面內(nèi)與x軸垂直的方向作為y軸。梁的截面尺寸為w×h，長度為l，密度為ρ。采用歐拉-伯努利梁理論，并將阻尼視為線性粘性阻尼。若阻尼系數(shù)為c，在單位長度推力f(x,t)作用下，懸臂梁彎曲振動方程[18]為：

(1)

圖1 懸臂梁結(jié)構(gòu)模型Fig.1 Cantilever structure model

式中：y(x,t)為x處撓度;E為楊氏模量;I為截面二次矩;ρl為單位長度質(zhì)量。

懸臂梁振動是由無窮多個主振動的疊加：

(2)

式中：φi(x)、Qi(t)分別為i階主振動的模態(tài)函數(shù)、時間函數(shù)。i階固有振動頻率ωi為：

(3)

式中：ω1稱為振動基頻，ωi(i≥2)稱為振動高頻。βl由頻率方程確定：

cos(βl)cosh(βl)+1=0

(4)

模態(tài)函數(shù)φi(x)為：

φi(x)= cos(βix)-cosh(βix)+

ξi[sin(βix)-sinh(βix)]

(5)

其中，

(6)

為研究方便，將噴射羽流沖擊力等效為作用于lf位置處的集中載荷f0δ(x-lf)。懸臂梁初始為靜止?fàn)顟B(tài)，在位置x=lf處，存在瞬間作用的脈沖力f(lf,t)，其瞬間作用沖量為S，單位長度的沖量為Sδ(x-lf)。

將阻尼系數(shù)c表示為c=2ρlζiωi，ζi為各階振動的阻尼比大小，可得系統(tǒng)傳遞函數(shù)G(x,s)為：

(7)

圖2 懸臂梁測量系統(tǒng)簡化結(jié)構(gòu)Fig.2 Simplified diagram of cantilever measuring system

2 推力測量方法

求得傳遞函數(shù)后，根據(jù)階躍推力f0δ(x-lf)下的系統(tǒng)響應(yīng)特性，確定推力測量方法。

2.1 階躍響應(yīng)特性分析

輸入f0δ(x-lf)時，R(s)=f0/s，系統(tǒng)輸出為：

(8)

根據(jù)式(8)，易得f(lf,t)=f0δ(x-lf)作用下，傳感器測量位置x=ls處系統(tǒng)響應(yīng)為:

(9)

(10)

式(10)表明，yi(ls,∞)與f0為線性關(guān)系，線性系數(shù)為i階增益值Ki(ls,lf)。根據(jù)式(10)，系統(tǒng)穩(wěn)態(tài)位移為:

(11)

由式(11)，推力會引起穩(wěn)態(tài)彎曲量變化，穩(wěn)態(tài)彎曲量大小與推力為線性關(guān)系，線性系數(shù)為系統(tǒng)增益值，系統(tǒng)增益值為各階子系統(tǒng)增益值之和。

2.2 參數(shù)標(biāo)定及恒定力解算方法

根據(jù)式(11)，輸出恒定標(biāo)定力f0，懸臂梁振動達到穩(wěn)態(tài)時，測量出y(ls,∞)，即可標(biāo)定出K(ls,lf)。標(biāo)定出K(ls,lf)后，測量出待測推力的y(ls,∞)后，即可計算出f0大小。但實際測量中，推力作用較短或變化較快時，系統(tǒng)穩(wěn)態(tài)誤差較大，無法直接得到y(tǒng)(ls,∞)，需要確定y(ls,∞)的求解方法。

圖3 系統(tǒng)階躍響應(yīng)低通濾波前后信號變化Fig.3 Signal change before and after low pass filtering of system step response

濾波后，消除1階動態(tài)分量影響，即可得到y(tǒng)(x,∞)。動態(tài)分量中，e-ζiωit為衰減項，由于懸臂梁阻尼比較小，通常在10-3量級，在小的時間區(qū)間內(nèi)，動態(tài)分量接近等幅振蕩。選取1階響應(yīng)振動周期Td1作為區(qū)間長度，對濾波后系統(tǒng)響應(yīng)yFilter(x,t)采用末端求均值來估算y(x,∞)，如圖4所示，推力作用在P2點停止，通過Td1確定區(qū)間起始點P1，對P1到P2區(qū)間內(nèi)位移求平均值可估算出P2點對應(yīng)的y(x,∞)。由于阻尼的作用，動態(tài)分量并不是等幅振蕩，而是按幾何級數(shù)衰減，P1到P2區(qū)間內(nèi)均值與y(x,∞)會存在一定誤差，顯然，推力作用時間tf越大，振幅衰減程度越大，P1到P2區(qū)間內(nèi)幅值差異越小，采用末端均值法估算y(x,∞)精度越高。

圖4 求取穩(wěn)態(tài)位移的末端均值法Fig.4 The end mean method for obtaining the steady state displacement

2.3 動態(tài)推力解算方法

當(dāng)推力變化時，穩(wěn)態(tài)分量y(x,∞)為時變值，動態(tài)分量也為時變值，但由于阻尼比不變，采用末端均值法仍可消除動態(tài)分量影響，得到時變y(x,∞)。因而，對于系統(tǒng)響應(yīng)中任意時刻的y(x,∞)求解(等效于推力作用在該時刻停止)，向前以Td1為區(qū)間進行積分求均值，即可估算出該時刻y(x,∞)。例如，推力作用在圖4中P4點停止時，在P3到P4區(qū)間內(nèi)求均值可估算出P4點y(x,∞)，根據(jù)式(11)，即可求解出tf時刻推力。

綜合以上分析，tf時刻y(x,∞)估計值可表示為:

(12)

實際中，位移傳感器測量得到的響應(yīng)數(shù)據(jù)為離散點，傳感器采樣率為ωs，Td1區(qū)間內(nèi)數(shù)據(jù)點個數(shù)m=ωsTd1，tf時刻對應(yīng)數(shù)據(jù)個數(shù)n=ωstf，若濾波后系統(tǒng)響應(yīng)位移數(shù)組為yFilter(i)，則y(x,∞)估算值為:

(13)

則t時刻推力f0(t)估計值為:

(14)

基于確定的動態(tài)推力解算及參數(shù)標(biāo)定方法，可測量步驟如下：

1)標(biāo)定系統(tǒng)增益值K(ls,lf)。通過施加長時間恒定標(biāo)定力f0，測量y(x,∞)，計算出K(ls,lf)。

2)對系統(tǒng)響應(yīng)進行濾波，通過末端均值法求出時變y(x,∞)，根據(jù)式(14)，計算f0(t)。

本質(zhì)上，消除高階噪聲影響后，懸臂梁結(jié)構(gòu)可視為固有頻率較高、阻尼比較小的典型的二階欠阻尼系統(tǒng)進行推力測量。由二階欠阻尼系統(tǒng)頻率響應(yīng)特性可得，當(dāng)推力器脈沖頻率高于懸臂梁固有頻率10倍以上時，穩(wěn)態(tài)輸出近似恒定[19]，因而仍可采用本文方法測量平均推力。

3 試驗驗證

通過搭建試驗平臺，對測量方法進行驗證。

3.1 試驗平臺

測量試驗平臺由推進器、控制模塊及測量模塊組成，如圖5所示。冷氣推進器樣機輸出推力在微牛至毫牛之間。測量模塊由懸臂梁、電動位移臺、標(biāo)定力產(chǎn)生裝置(磁鐵與線圈)、位移傳感器組成，主要作用是將沖擊力轉(zhuǎn)變?yōu)槲灰菩盘?；控制模塊由位移臺控制箱、傳感器控制箱及前端放大器、電源、計算機、遙控開關(guān)、無線壓力接收模塊組成，主要用于控制測量模塊及推進器工作，接收測量模塊輸出的位移信號以及推進器輸出的入口壓力值。

建立測量模型時將推力等效為集中載荷，因而推進器噴口面積要小。測量模塊及推進器如圖6所示，整個推進器集成于一個平板上，噴嘴直徑為1 mm。噴口前壓力值表征推進器工況的重要參數(shù)，無線壓力變送器測量并每隔1 s發(fā)送一次噴口前壓力值，計算機接收并顯示壓力值。

懸臂梁(尺寸為210 mm×40 mm×0.3 mm，材料為304不銹鋼)與線圈、位移傳感器探頭通過同一個光學(xué)平板固定在電動位移臺上，電動位移臺用于調(diào)節(jié)噴口與懸臂梁之間的距離。距離不同會引起噴管內(nèi)氣體膨脹狀態(tài)的不同，導(dǎo)致作用在懸臂梁上的氣體沖擊力大小的不同[17]。測量時，懸臂梁的寬度必須能包容全部射流區(qū)域。由式(3)，懸臂梁振動頻率與寬度無關(guān)，為了包容全部射流區(qū)域，可將寬度做得很大。但寬度增大，會引起系統(tǒng)增益值增大，導(dǎo)致位移量減小，測量靈敏度降低。因而理論上，寬度選擇應(yīng)為滿足包容全部射流區(qū)域的最小值。冷氣噴管出口擴張角度為12°，距離變化最大值為100 mm，計算得寬度最小值為36.3 mm。為留有一定裕度，選擇寬度40 mm。粘貼在懸臂梁上的磁鐵(Φ6×1.5 mm，質(zhì)量0.29 g，材料為N35)與線圈(30匝)構(gòu)成標(biāo)定力產(chǎn)生裝置。經(jīng)過電子天平標(biāo)定，線圈電流0.1A時，電磁力大小為0.511 mN。位移傳感器為電容式位移傳感器(分辨力30 nm，量程1 mm)，可以直接以數(shù)字量形式輸出懸臂梁響應(yīng)位移大小。

圖5 基于懸臂梁結(jié)構(gòu)的推力測量平臺結(jié)構(gòu)Fig.5 Structure map of thrust measurement platform based on cantilever structure

圖6 基于懸臂梁結(jié)構(gòu)的測量系統(tǒng)Fig.6 Measuring system based on cantilever structure

3.2 試驗結(jié)果

首先通過施加電磁標(biāo)定力標(biāo)定出系統(tǒng)增益值，然后控制推進器輸出推力，得到系統(tǒng)位移響應(yīng)，采用末端均值法求取時變穩(wěn)態(tài)位移值，結(jié)合標(biāo)定出的狀態(tài)增益值，求取時變推力。

3.2.1 系統(tǒng)增益值標(biāo)定

標(biāo)定中，設(shè)置位移測量臂長ls=106 mm，采樣率ωs=1041.7 Hz，標(biāo)定力f0力臂lf=174 mm時，控制線圈電流依次輸出0.1 A、0.2 A及0.3 A，對應(yīng)f0為0.511 mN、1.022 mN及1.533 mN(斥力)，然后逐漸減小到0.1 A時，系統(tǒng)響應(yīng)如圖7所示。

圖7 標(biāo)定力逐漸增大、減小時懸臂梁響應(yīng)Fig.7 Response of cantilever when gradually increasing and reducing the calibration force

圖7中，標(biāo)定力作用下，系統(tǒng)經(jīng)過動態(tài)振動，最終穩(wěn)定在新的平衡位置，穩(wěn)定時間約為50 s。卸載標(biāo)定力后，系統(tǒng)回到初始平衡位置，表明系統(tǒng)零漂較小。不同大小f0穩(wěn)態(tài)位移y(ls,∞)分別為25.3101 μm、60.6212 μm、75.6327 μm，根據(jù)式(11)，可得系統(tǒng)增益值K(ls,lf)=4.953×10-2。

3.2.2 推力測量結(jié)果

推進器作用力臂lf=135 mm，噴口入口壓力0.45 MPa。懸臂梁與噴口距離24 mm時，系統(tǒng)響應(yīng)如圖8所示。推力作用下，懸臂梁開始彎曲振動，平均位置基本穩(wěn)定在90 μm左右。在20 s之前，振幅整體上呈現(xiàn)衰減趨勢；20 s之后，振幅基本穩(wěn)定，但波動較大。推進器停止工作后，懸臂梁開始自由振動，振幅逐漸衰減，最終穩(wěn)定在初始平衡位置。推力停止輸出后，懸臂梁振幅按指數(shù)規(guī)律衰減，而推進器工作時，振幅波動明顯大于環(huán)境噪聲，說明振幅的波動是由推力輸出的不穩(wěn)定引起的。雖然存在振幅的波動，但平衡位置并沒有發(fā)生變化，說明推力輸出不穩(wěn)定為隨機突變式，相當(dāng)于脈沖力的作用，因而不會引起平衡位置的改變。對圖8進行頻譜分析得到ωd1=6.401 6 Hz，ωd2=40.879 1 Hz，設(shè)置濾波時截止頻率為8 Hz。

圖8 推力作用下系統(tǒng)響應(yīng)曲線Fig.8 System response curve under thrust

采用末端均值法，根據(jù)濾波后響應(yīng)曲線求取y(ls,∞)隨時間的變化與響應(yīng)數(shù)據(jù)的對比如圖9所示。y(ls,∞)整體上保持平穩(wěn)，并存在較小的波動。y(ls,∞)平均值為-92.013 5 μm，標(biāo)準差為0.979 6 μm。

圖9 推力響應(yīng)與求取的穩(wěn)態(tài)位移對比Fig.9 Comparison of thrust response and the calculated steady state displacement

由于K(ls,lf)=4.953×10-2，根據(jù)式(11)，計算得到的推力隨時間變化曲線如圖10所示。推力開始作用的初始0.5 s，推力逐漸增大，而后逐漸平穩(wěn)。計算得推力平均值為2.584 mN，標(biāo)準差為0.028 mN，推力相對于平均值的最大變化不超過0.080 mN。

圖10 推力隨時間變化曲線Fig.10 Thrust curve in variation of time

4 誤差分析

采用非接觸測量方式時，由于發(fā)散及反射氣流對來流干擾等因素影響，實際上懸臂梁測量得到的是尾流沖擊力，與實際推力存在一定誤差。將作用在懸臂梁上的沖擊力的測量誤差記為US1，沖擊力與實際推力誤差記為US2，推力測量誤差記為US。

4.1 沖擊力測量誤差

首先來分析US1大小。式(14)中，K(ls,lf)標(biāo)定誤差來源包括y(ls,∞)求取誤差，以及標(biāo)定力f0本身的輸出誤差。y(ls,∞)誤差等于位移傳感器測量誤差0.050%。設(shè)置相同采樣率及懸臂梁參數(shù)進行理論計算，濾波所對應(yīng)的幅值計算誤差為-0.470%。綜合位移測量及濾波誤差，可認為y(ls,∞)誤差為:

(15)

標(biāo)定力誤差為1.210%，標(biāo)定及沖擊力解算中均需要求解y(ls,∞)，則認為:

(16)

4.2 推力與沖擊力關(guān)系

在相同的測量環(huán)境及推進器工況條件下，采用試驗手段，以精度較高的扭擺系統(tǒng)的測量結(jié)果為推力基準，與懸臂梁測量結(jié)果進行對比，分析得到US2。

顯然，當(dāng)噴口與懸臂梁距離變化時，沖擊力大小也是不同的。在2～100 mm之間，保持入口壓力0.45 MPa，以回進程方式，通過電位移臺逐漸增大、減小距離，每隔2 mm測量5次取平均值，得到的沖擊力及誤差值如圖11所示，隨著距離變化，沖擊力大小變化明顯，最大值約為最小值的3倍。

圖11 推力測量值隨噴口與懸臂梁距離的變化Fig.11 Measured thrust in variation of distance between nozzle and cantilever

按照沖擊力變化規(guī)律，將距離分為3個區(qū)段：

1) I區(qū)段為沖擊力上升段(2～24 mm)，沖擊力隨著距離的增大而增大，尤其是10～14 mm區(qū)間內(nèi)，沖擊力快速增大，較小的距離差值也會引起沖擊力的明顯變化，因而測量誤差也相對較大。

2) II區(qū)段為沖擊力平穩(wěn)段(24～42 mm)，沖擊力隨距離變化較小，平均值為2.597 mN，變異系數(shù)為0.49%。在此區(qū)間內(nèi)測量時，能夠測量得到較為穩(wěn)定的沖擊力值，然后得到該區(qū)間內(nèi)沖擊力與實際推力的關(guān)系。測量推力時，首先測量得到?jīng)_擊力后，根據(jù)確定的關(guān)系即可計算出實際推力。

3) III區(qū)段為減小段(42～100 mm)，沖擊力隨著距離的增大而減小，且區(qū)間內(nèi)誤差較大，可能的原因是距離較大時，尾流擴散時易受隨機誤差影響。

通過減壓閥調(diào)整入口壓力，得到平穩(wěn)段位置基本保持不變，說明噴口形狀一定時，沖擊力平穩(wěn)段位置也是一定的。得到?jīng)_擊力與距離關(guān)系后，以下通過扭擺來測量實際推力大小。

扭擺原理如圖12所示，推進器與配重在扭轉(zhuǎn)橫梁的兩端平衡。旋轉(zhuǎn)橫梁與轉(zhuǎn)軸固定，轉(zhuǎn)軸通過上下兩個撓性樞軸(用于提供回復(fù)力)與固定部件相連。推進器輸出推力時，推力轉(zhuǎn)化為橫梁的線位移，隨著位移的增大，推力與樞軸回復(fù)力逐漸趨于平衡。位移傳感器測量橫梁穩(wěn)態(tài)線位移，結(jié)合電磁標(biāo)定力標(biāo)定出的扭轉(zhuǎn)剛度系數(shù)，即可解算出推力大小。

圖12 扭擺測量推力原理Fig.12 Thrust measurement principle of torsion balance

標(biāo)定得到，扭擺有阻尼振動周期為11.910 s，振動頻率為8.40×10-2Hz，阻尼比為0.283。相比之下，扭擺系統(tǒng)的阻尼比、振動頻率均比懸臂梁大2個數(shù)量級。施加1.063 mN標(biāo)定力后，控制推進器輸出推力，扭擺系統(tǒng)響應(yīng)如圖13所示，推力作用下，扭擺響應(yīng)特性為二階系統(tǒng)典型階躍響應(yīng)。

圖13 標(biāo)定力響應(yīng)與推力響應(yīng)對比Fig.13 Response comparison of calibration force and thrust

結(jié)合標(biāo)定出的剛度系數(shù)，同時考慮測量臂長及推力力臂比例關(guān)系在內(nèi)，并考慮扭擺測量誤差，解算得推力大小為2.707 mN。保持測量條件不變，經(jīng)過5次測量，計算得推力平均值為2.707 mN，標(biāo)準差為9.85×10-3mN。因而有，懸臂梁測量結(jié)果小于扭擺測量結(jié)果，平均值相差-4.06%。也就是說，在平穩(wěn)段測量得到的沖擊力與實際推力誤差為-4.064%，因而US2=-4.064%。

同時，對比圖13與圖8可得，扭擺系統(tǒng)振動頻率較小，對于高頻擾動具有一定的抑制作用，而懸臂梁振動頻率較大，響應(yīng)較快，能夠響應(yīng)推力輸出不穩(wěn)定波動。

確定US1及US2后，可得懸臂梁結(jié)構(gòu)測量推進器階躍力誤差US：

(17)

5 結(jié)束語

本文確定了基于末端均值法求時變穩(wěn)態(tài)位移，進而計算動態(tài)推力的方法，得到的主要結(jié)論如下：

1) 懸臂梁測量系統(tǒng)在結(jié)構(gòu)上為無窮多個欠阻尼二階子系統(tǒng)的并聯(lián)，系統(tǒng)增益值與剛度系數(shù)大小為倒數(shù)關(guān)系。

2) 懸臂梁系統(tǒng)阻尼比較小，位移響應(yīng)中的動態(tài)分量在較小的時間區(qū)間內(nèi)接近等幅振蕩。通過低通濾波消除高階動態(tài)分量的影響，以一階振蕩周期為區(qū)間長度求平均值可消除動態(tài)分量影響，得到精度較高的穩(wěn)態(tài)位移時變值。

3) 懸臂梁與噴口之間距離對沖擊力影響較大，隨著距離的增大，沖擊力先增大后減小。在一定距離區(qū)間內(nèi)，會出現(xiàn)一段平穩(wěn)區(qū)，在此區(qū)間內(nèi)進行測量，可以得到穩(wěn)定的沖擊力數(shù)值。入口壓力0.45 MPa時，在距離為24～40 mm區(qū)間內(nèi)，出現(xiàn)沖擊力穩(wěn)定區(qū)，推力平均值為2.597 mN，變異系數(shù)為0.49%。

4) 與扭擺測量結(jié)果對比得到，沖擊力測量值小于實際推力。推力測量誤差由沖擊力測量誤差及沖擊力與真實推力之間的誤差組成。

結(jié)合仿真分析，研究射流與懸臂梁之間其他作用形式對測量結(jié)果影響，以及進一步細化誤差分析，是下一步的研究重點。