一種基于懸臂梁結構的動態推力測量方法

，，，，，

1. 航天工程大學 激光推進及其應用國家重點實驗室，北京 101416 2. 西安衛星測控中心，西安 710043

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

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

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

本文在懸臂梁動力學模型的基礎上，確定了系統傳遞函數，根據系統響應特征，確定了動態推力測量方法。同時，通過試驗驗證了測量方法的正確性，并通過與扭擺測量結果進行對比，分析了測量誤差。

1 系統傳遞函數的建立

單位脈沖函數拉普拉斯變換為1，求得單位脈沖力下的系統響應，做拉普拉斯變換即可求得系統傳遞函數。

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

(1)

圖1 懸臂梁結構模型Fig.1 Cantilever structure model

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

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

(2)

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

(3)

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

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

(4)

模態函數φi(x)為：

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

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

(5)

其中，

(6)

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

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

(7)

圖2 懸臂梁測量系統簡化結構Fig.2 Simplified diagram of cantilever measuring system

2 推力測量方法

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

2.1 階躍響應特性分析

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

(8)

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

(9)

(10)

式(10)表明，yi(ls,∞)與f0為線性關系，線性系數為i階增益值Ki(ls,lf)。根據式(10)，系統穩態位移為:

(11)

由式(11)，推力會引起穩態彎曲量變化，穩態彎曲量大小與推力為線性關系，線性系數為系統增益值，系統增益值為各階子系統增益值之和。

2.2 參數標定及恒定力解算方法

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

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

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

圖4 求取穩態位移的末端均值法Fig.4 The end mean method for obtaining the steady state displacement

2.3 動態推力解算方法

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

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

(12)

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

(13)

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

(14)

基于確定的動態推力解算及參數標定方法，可測量步驟如下：

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

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

本質上，消除高階噪聲影響后，懸臂梁結構可視為固有頻率較高、阻尼比較小的典型的二階欠阻尼系統進行推力測量。由二階欠阻尼系統頻率響應特性可得，當推力器脈沖頻率高于懸臂梁固有頻率10倍以上時，穩態輸出近似恒定[19]，因而仍可采用本文方法測量平均推力。

3 試驗驗證

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

3.1 試驗平臺

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

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

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

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

圖6 基于懸臂梁結構的測量系統Fig.6 Measuring system based on cantilever structure

3.2 試驗結果

首先通過施加電磁標定力標定出系統增益值，然后控制推進器輸出推力，得到系統位移響應，采用末端均值法求取時變穩態位移值，結合標定出的狀態增益值，求取時變推力。

3.2.1 系統增益值標定

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

圖7 標定力逐漸增大、減小時懸臂梁響應Fig.7 Response of cantilever when gradually increasing and reducing the calibration force

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

3.2.2 推力測量結果

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

圖8 推力作用下系統響應曲線Fig.8 System response curve under thrust

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

圖9 推力響應與求取的穩態位移對比Fig.9 Comparison of thrust response and the calculated steady state displacement

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

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

4 誤差分析

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

4.1 沖擊力測量誤差

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

(15)

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

(16)

4.2 推力與沖擊力關系

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

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

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

按照沖擊力變化規律，將距離分為3個區段：

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

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

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

通過減壓閥調整入口壓力，得到平穩段位置基本保持不變，說明噴口形狀一定時，沖擊力平穩段位置也是一定的。得到沖擊力與距離關系后，以下通過扭擺來測量實際推力大小。

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

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

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

圖13 標定力響應與推力響應對比Fig.13 Response comparison of calibration force and thrust

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

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

確定US1及US2后，可得懸臂梁結構測量推進器階躍力誤差US：

(17)

5 結束語

本文確定了基于末端均值法求時變穩態位移，進而計算動態推力的方法，得到的主要結論如下：

1) 懸臂梁測量系統在結構上為無窮多個欠阻尼二階子系統的并聯，系統增益值與剛度系數大小為倒數關系。

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

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

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

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