CN114218750B - 基于数字滤波器的星载微推力器推力响应时间测量方法 - Google Patents

基于数字滤波器的星载微推力器推力响应时间测量方法 Download PDF

Info

Publication number
CN114218750B
CN114218750B CN202111352455.6A CN202111352455A CN114218750B CN 114218750 B CN114218750 B CN 114218750B CN 202111352455 A CN202111352455 A CN 202111352455A CN 114218750 B CN114218750 B CN 114218750B
Authority
CN
China
Prior art keywords
thrust
frequency
response
measurement system
response time
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202111352455.6A
Other languages
English (en)
Other versions
CN114218750A (zh
Inventor
洪延姬
冯孝辉
金星
叶继飞
崔海超
冯高平
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Original Assignee
Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Peoples Liberation Army Strategic Support Force Aerospace Engineering University filed Critical Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Priority to CN202111352455.6A priority Critical patent/CN114218750B/zh
Publication of CN114218750A publication Critical patent/CN114218750A/zh
Priority to US17/964,466 priority patent/US20230150696A1/en
Application granted granted Critical
Publication of CN114218750B publication Critical patent/CN114218750B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
    • B64G1/36—Guiding or controlling apparatus, e.g. for attitude control using sensors, e.g. sun-sensors, horizon sensors
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H17/00—Networks using digital techniques
    • H03H17/02—Frequency selective networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/12—Timing analysis or timing optimisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Remote Sensing (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Chemical & Material Sciences (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Combustion & Propulsion (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明属于航天卫星推进技术领域，具体涉及一种基于数字滤波器的星载微推力器推力响应时间测量方法。本发明的推力响应时间测量方法包括步骤，S1：对扭摆推力测量系统的非零初始条件进行归零处理,得到变量代换后的推力测量系统的振动微分方程；S2：在变量代换后的推力测量系统后串联一个数字滤波器,得到等效灵敏度高频推力测量系统；S3：确定等效灵敏度高频推力测量系统的系统响应；S4：从系统响应判读星载微推力器的推力响应时间；S5：利用系统响应反演计算待测推力，进一步确认推力响应时间。本发明的等效灵敏度高频推力测量系统的固有频率大幅提高并极大缩短了进入稳态时间，实现了快变推力响应时间的测量。

Description

基于数字滤波器的星载微推力器推力响应时间测量方法
技术领域
本发明涉及航天卫星推进技术领域，特别是一种基于数字滤波器的星载微推力器推力响应时间测量方法。
背景技术
在卫星精确控制场合，对微推力器提出了快速调节推力要求，需要将推力响应时间控制在指定范围内。推力响应时间可用推力上升时间和推力下降时间表示，当推力从较小稳定推力变化到较大稳定推力时，推力响应时间采用推力上升时间表示；当推力从较大稳定推力变化到较小稳定推力时，推力响应时间采用推力下降时间表示，此时，推力响应时间是微推力器的考核指标之一。
推力响应时间的测量难点在于：推力变化频率远大于推力测量系统的固有频率，或推力测量系统的固有频率太低，不满足快变推力响应时间测量要求。首先，由于用于微推力测量的推力测量系统的高灵敏度要求，推力测量系统的刚度系数很小，使得固有频率降低；其次，由于推力器搭载在推力测量系统上，推力测量系统的转动部件(推力器、配重、结构梁等)的转动惯量很大，也使得固有频率降低；最后，如果测量10ms级推力响应时间，则要求推力测量系统的固有线频率为100Hz以上，这是目前任何(扭摆、吊摆、倒摆、弯振)推力测量系统，在物理结构上都难以实现的。
现有的星载微推力器推力响应时间测量方法主要缺陷如下：
(1)不了解推力响应时间与推力测量系统的系统响应之间内在联系。实际上，推力测量系统的系统响应，包括暂态系统响应和稳态系统响应，只有推力响应时间大于推力测量系统的暂态持续时间(或进入稳态时间)，才能从推力测量系统的系统响应反映出推力响应时间。
(2)针对快变推力响应时间的动态推力计算方法不成熟，难以反映准确的推力上升过程或推力下降过程。在动态推力计算中，以差商作为导数的近似值反向求解振动微分方程(振动方程反问题求解)，一是误差大：差商作为导数的近似值，一阶导数误差很大，二阶导数误差更大；二是收敛性差：步长较大时截断误差增大，步长较小时舍入误差增大。
发明内容
有鉴于此，本发明的目的在于解决快变推力响应时间的测量问题。
为解决上述技术问题，本发明提供了一种基于数字滤波器的星载微推力器推力响应时间测量方法，包括以下步骤：
S1：对扭摆推力测量系统的非零初始条件进行归零处理,得到变量代换后的推力测量系统的振动微分方程；
S2：在变量代换后的推力测量系统后串联一个数字滤波器,得到等效灵敏度高频推力测量系统；
S3：确定等效灵敏度高频推力测量系统的系统响应；
S4：从系统响应判读星载微推力器的推力响应时间；
S5：利用等效灵敏度高频推力测量系统反演计算待测推力，进一步确认推力响应时间。
步骤S1中，在待测推力f′(t)作用下系统响应θ′(t)由位移传感器测量得到，根据非零初始条件进行变量代换，令 得到变量代换后的推力测量系统的振动微分方程为：
f(t)＝f′(t)+f0(t)
其中，f0(t)是与非零初始条件相关的等效推力，f′(t)为待测推力，θ′(t)为该推力作用下由位移传感器测量得到的系统响应，ζ为阻尼比，ωn为固有振动频率，J为转动惯量，Lf为力臂；θ0表示初始扭转角，θ0和为常数，θ(t)为变量代换后系统响应，满足初始条件θ(0)＝0和
步骤S2中，数字滤波器的传递函数为：
数字滤波器的单位脉冲响应函数为
其中，ωn为固有振动频率，ζ为阻尼比，ωn1为等效灵敏度高频推力测量系统的固有频率，ωn1＝Cωωn，Cω为频率比，ζ1为等效灵敏度高频推力测量系统的阻尼比，为改进后振动频率。其中，Cω＞1，ζ1＝0.7～0.9。
步骤S3中，数字滤波器的输入为θ(τ)，通过位移传感器测量值θ′(τ)计算得到θ(τ)，数字滤波器输出的系统响应θ1(t)为：其中,τ为积分变量。
步骤S4中，当系统响应与稳态系统响应比较相差1％时所对应的时间为推力响应时间。
推力响应时间包括推力上升时间和推力下降时间。
步骤S5中，等效灵敏度高频推力测量系统的单位脉冲响应函数h1(t)为
已知数字滤波器的输出为θ1(t)，利用等效灵敏度高频推力测量系统的推力积分方程：计算名义推力f(τ)。
待测推力f′(t)的估计值F′(t)为：F′(t)＝F(t)-f0(t),其中，F(t)是名义推力f(t)的估计值，ζ为阻尼比，ωn为固有振动频率，J为转动惯量，Lf为力臂，θ0表示初始扭转角，θ0和为常数；当系统响应与稳态系统响应比较相差1％时所对应的时间为推力响应时间。
有益效果
本发明通过在固有频率低的推力测量系统上串联数字滤波器，构建了等效灵敏度高频推力测量系统，采用数字滤波器的设计方法，使得等效灵敏度高频推力测量系统的固有频率大幅提高并极大缩短了进入稳态时间，实现了快变推力响应时间的测量，并且等效灵敏度高频推力测量系统与原推力测量系统的灵敏度相等，实现了等灵敏度测量。通过数字滤波器的设计，可调节振动频率和阻尼比，根据推力上升时间和推力下降时间的测量需求，通过调节振动频率和阻尼比，实现快变推力上升时间和推力下降时间的测量。
附图说明
图1是本发明推力响应时间测量方法的流程图。
图2是本发明等效灵敏度高频推力测量系统示意图。
图3是推力上升时间的变化图。
图4是原推力测量系统的系统响应测量值图。
图5是变量代换后系统响应的变化图。
图6是数字滤波器的系统响应变化图。
图7是推力上升时间25ms下反演计算的推力和误差图。
图8是推力上升时间50ms下反演计算的推力和误差图。
图9是推力上升时间100ms下反演计算的推力和误差图。
图10是推力逐渐下降时数字滤波器的系统响应变化图。
图11是推力逐渐下降时反演计算的推力图。
图12是推力逐渐下降时反演计算的推力误差图。
具体实施方式
下面将参考附图对本发明的具体实施方式进行详细说明。
如图1所示，是本发明推力响应时间测量方法的流程图，下面对实施过程进行详细说明。
1.建立推力响应时间与推力测量系统的进入稳态时间关系
建立推力响应时间与推力测量系统的进入稳态时间关系，是推力测量系统能够测量给定推力响应时间的必要条件。
1)推力响应时间与振动微分方程初始条件
微推力器的推力响应时间不但需要考虑推力上升时间还要考虑推力下降时间。推力由较小稳定推力上升到较大稳定推力时，推力响应时间可用推力上升时间表示；推力由较大稳定推力下降到较小稳定推力时，推力响应时间可用推力下降时间表示。研究推力响应时间必须考虑初始条件的影响。
扭摆推力测量系统的振动微分方程为
式中，f′(t)为待测推力，θ′(t)为该推力作用下系统响应(由位移传感器测量得到)，ζ为阻尼比，ωn为固有振动频率，J为转动惯量，Lf为力臂。
当推力从一个稳定推力变化到另一个稳定推力时，初始条件不为零，为
式中，若是从一个非零稳定推力开始变化，可近似认为若是从零推力开始变化，可近似认为θ′(0)≈0和
显然，推力测量系统的振动微分方程的初始条件，反映了推力是从怎样初始推力开始变化的，例如，从非零稳定推力开始变化，还是从零推力开始变化。因此，测量推力响应时间必须明确推力测量系统的振动微分方程的初始条件。
2)建立推力响应时间与进入稳态时间关系
扭摆推力测量系统的振动微分方程为
如果初始条件不为零，为θ′(0)＝θ0和设推力的拉普拉斯变换为L[f′(t)]＝F′(s)，系统响应及导数的拉普拉斯变换为
L[θ′(t)]＝Θ′(s) (4)
振动方程两边取拉普拉斯变换，整理得
两边取拉普拉斯反变换，可得
令
θ′(t)＝θ′1(t)+θ′2(t) (9)
在θ1′(t)中包括暂态系统响应和稳态系统响应。例如，在θ1′(t)中代入周期推力f′(t)＝Af cosωt(Af为常数表示周期推力的振幅，ω为周期推力的频率)，推力测量系统的系统响应，为
其中
式中，S＝Lf/k为推力测量系统的灵敏度，为扭转刚度系数，Rω＝ω/ωn为频率比。显然，在θ1′(t)中系统响应的第一项为暂态系统响应，代表暂态过程随着时间变化，暂态过程的振动频率为推力测量系统的振动频率；第二项为稳态系统响应，代表稳态过程随着时间变化，稳态过程的振动频率为周期推力的振动频率。θ2′(t)也为暂态系统响应，其振动频率为推力测量系统的振动频率。
总之，暂态系统响应的振幅是逐渐衰减的，衰减项为设衰减项时所对应时间(暂态振幅衰减到1％所对应时间)，定义为暂态持续时间或进入稳态时间ts，有
式中，Tn＝2π/ωn为推力测量系统的固有周期。
推力测量系统的进入稳态时间或暂态持续时间，是推力测量系统的固有属性，具有以下结论：
(1)推力测量系统的固有频率越大则周期越小，进入稳态时间越短；
(2)推力测量系统的阻尼比越小，进入稳态时间越大；
(3)推力响应时间只有大于推力测量系统的进入稳态时间条件下，才能从系统响应曲线表现出推力响应时间的变化特点，或才能从系统响应曲线判读推力响应时间。
例如，在常用阻尼比0.1≤ζ≤0.3范围内，推力测量系统的进入稳态时间为固有周期的3～7倍，设推力上升时间为ts，推力测量系统的固有周期为Tn，则有
ts≥7Tn (16)
如果要求推力上升时间ts≤50ms，则要求推力测量系统的固有周期为Tn≤7.1ms，也就是线频率为fn≥140Hz，这是任何结构的推力测量系统都难以做到的。说明测量快变推力响应时间，仅考虑从物理结构设计上提高推力测量系统的固有频率是远远不够的，需要另辟蹊径。
2.非零初始条件的归零处理方法
在推力响应时间测量中初始条件不为零，但是，在后续数字滤波器设计和构建等效灵敏度高频推力测量系统中，由于利用传递函数建立模型，必须将非零初始条件进行归零处理。
扭摆推力测量系统的振动微分方程为
初始条件不为零，为
式中，θ0和为常数。
为了保证初始条件为零，进行变量代换，令
代入振动方程，整理得
代入初始条件，可得
将变量代换后的推力测量系统的振动微分方程，改写为
f(t)＝f′(t)+f0(t) (25)
式中，f0(t)是与非零初始条件相关的等效推力。
当推力从一个稳定推力变化到另一个稳定推力时，有θ′(0)＝θ0和可得
θ(t)＝θ(t)-θ0 (27)
说明：代换后系统响应θ(t)与代换前系统响应θ′(t)，暂态系统响应和稳态系统响应都相同只是系统响应曲线上下移动一定距离。
非零初始条件的归零处理方法为：
(1)在待测推力f′(t)作用下系统响应θ′(t)是由位移传感器测量得到，根据非零初始条件进行变量代换，为
式中，θ(t)为变量代换后系统响应，满足初始条件为零θ(0)＝0和条件。
(2)当推力从一个稳定推力变化到另一个稳定推力时(θ′(0)＝θ0和)，变量代换后系统响应θ(t)与代换前系统响应θ′(t)的进入稳态时间相同，因此，从代换后系统响应θ(t)曲线可判读推力响应时间。
(3)根据变量代换后系统响应θ(t)，利用变量代换后振动微分方程，可计算名义推力f(t)，从而求得待测推力f′(t)为
f′(t)＝f(t)-f0(t) (29)
式中，f0(t)是与非零初始条件相关的等效推力。
3.利用数字滤波器构建等效灵敏度高频推力测量系统
为了解决推力测量系统固有频率过低，不能满足快变推力响应时间测量的难题，采用数字滤波器构建等效灵敏度的高频推力测量系统，大幅提高固有频率，满足快变推力响应测量的需求。
如图2所示，是本发明等效灵敏度高频推力测量系统示意图。采用数字滤波器构建等效灵敏度高频推力测量系统的方法为：
(1)在变量代换后的推力测量系统的传递函数上，串联一个数字滤波器，构造等效灵敏度高频推力测量系统；
(2)所构造的等效灵敏度高频推力测量系统，灵敏度与变量代换前的推力测量系统相同，但是，通过数字滤波器的设计，可调节振动频率和阻尼比；
(3)根据推力上升时间和推力下降时间的测量需求，通过调节振动频率和阻尼比，实现推力上升时间和推力下降时间的测量。
4.数字滤波器设计方法与系统响应计算方法
1)数字滤波器设计方法
设变量代换后的推力测量系统的输入为名义推力f(t)(拉普拉斯变化为F(s))，其系统响应为θ(t)(由位移传感器测量值θ′(t)计算得到，拉普拉斯变化为Θ(s))；数字滤波器的传递函数为D(s)，数字滤波器的输入为θ(t)，数字滤波器的输出为θ1(t)(由θ(t)计算得到，拉普拉斯变化为Θ1(s))；所构造的等效灵敏度高频推力测量系统的传递函数为H1(s)＝H(s)D(s)，等效灵敏度高频推力测量系统的输入为名义推力f(t)，输出为θ1(t)。
扭摆推力测量系统的传递函数为
式中，Lf为力臂，J为转动惯量，ζ为阻尼比，ωn为固有振动频率，为振动频率，为扭转刚度系数。
数字滤波器设计方法是采用数字滤波器与变量代换后的推力测量系统的输出端串联，构成改进的扭摆推力测量系统。改进后扭摆推力测量系统的传递函数为
式中，改进后扭摆推力测量系统的固有振动频率为ωn1，阻尼比为ζ1，扭转刚度系数为k1。
设数字滤波器的传递函数为D(s)，则有
为了减小进入稳态时间提高推力响应时间的测量能力，需要增大固有振动频率，令ωn1＝Cωωn(Cω＞1)和k1＝k(灵敏度相等)，有
数字滤波器的传递函数为
式中，为改进后振动频率。进一步简化，可得数字滤波器的传递函数为
设数字滤波器的单位脉冲响应函数为d(t)，则有
可得数字滤波器的单位脉冲响应函数为
数字滤波器的输入为θ(τ)(通过位移传感器测量值θ′(τ)计算得到)，设数字滤波器的输出为θ1(t)，则有
从而，获得等效灵敏度高频推力测量系统的系统响应。
数字滤波器的设计方法为：
(1)以灵敏度相等和增大固有振动频率为目标，设计数字滤波器的传递函数；
(2)根据数字滤波器的传递函数，通过拉普拉斯反变换获得数字滤波器的单位脉冲响应函数；
(3)利用变量代换后的推力测量系统的系统响应，计算数字滤波器的系统响应。
2)数字滤波器的系统响应计算方法
数字滤波器的输出量需要采用其单位脉冲响应函数，利用数值积分方法进行计算。
设数字滤波器的单位脉冲响应函数为d(t)，有
令
d(t)＝d1(t)+d2(t) (42)
代入以下方程
可得
式中，可根据系统响应θ(t)直接计算，的计算需要采用数值积分方法。
为了计算的积分方程，设时间采样步长为h，ti＝ih(i＝0,1,2,…)，有
由于[t0＝0,ti＝ih]划分为i等分，对于τj＝jh(j＝0,1,2,…,i)，令函数
g(ti,τj)＝θ(τj)d2(ti-τj) (50)
的计算方法为采用梯形积分公式起步计算，再交替采用3点和4点辛普森积分公式计算，具体为：
(1)当子区间数目n＝1时，采用梯形积分公式起步计算。
(2)当子区间数目n＝2时，采用3点(2个区间)辛普森积分公式计算。
(3)当子区间数目n＝3时，采用4点(3个区间)辛普森积分公式计算。
(4)当子区间数目n＝2k(k＝2,3,4,…)时，采用3点辛普森积分公式计算。
用2k+1个节点ti＝ih(i＝0,1,2,…,2k)，将积分区间[t0,t2k]划分为2k等分，子区间长度为h＝t2k/(2k)，每个子区间上反复使用3点辛普森积分公式，可得复化辛普森公式为
(5)当子区间数目n＝2k+1(k＝2,3,4,…)时，交替采用3点和4点辛普森积分公式计算。
用2k+2个节点ti＝ih(i＝0,1,2,…,2k+1)，将积分区间[t0,t2k+1]划分为2k+1等分，子区间长度为h＝t2k+1/(2k+1)，在前面的2k-2(k＝2,3,4,…)个子区间反复使用3点辛普森积分公式，在后面2k-2、2k-1、2k、2k+1等4个点使用4点辛普森积分公式，可得复化辛普森公式为
在数字滤波器的系统响应计算中，变量代换前的推力测量系统的固有频率ωn和阻尼比ζ是已知的，根据推力响应时间测量需要，选取等效灵敏度高频推力测量系统的固有频率ωn1＝Cωωn(通过选取频率比Cω＞1实现)，选取阻尼比ζ1＝0.7～0.9进一步达到减小进入稳态时间的目的。
5.利用等效灵敏度高频推力测量系统反演计算推力方法
由于等效灵敏度高频推力测量系统的固有频率大幅提高，可实现快变推力响应时间的测量。
一方面，从数字滤波器的系统响应，也就是等效灵敏度高频推力测量系统的系统响应，可间接判读推力上升时间和推力下降时间；另一方面，根据数字滤波器的系统响应(也就是等效灵敏度高频推力测量系统的系统响应)，利用等效灵敏度高频推力测量系统的传递函数和单位脉冲响应函数，反演计算推力，可以进一步确认推力上升时间和推力下降时间。此时，推力响应时间的判读更为合理和可靠。
1)建立等效灵敏度高频推力测量系统的推力积分方程
等效灵敏度高频推力测量系统的传递函数为
等效灵敏度高频推力测量系统的单位脉冲响应函数h1(t)为
已知数字滤波器的输出为θ1(t)，设名义推力f(τ)，等效灵敏度高频推力测量系统的推力积分方程为
上述方程称为推力积分方程，根据推力积分方程，利用θ1(t)可反演计算名义推力f(τ)。
2)推力积分方程离散化计算推力
将等效灵敏度高频推力测量系统的单位脉冲响应函数h1(t)代入推力积分方程，可得
已知等效灵敏度高频推力测量系统的系统响应为[ti,θ1(ti)](i＝0,1,2,3,…)，t0＝0时θ1(t0)＝0，采样时间步长为h，ti＝ih，令θ1(ti)＝θ1i和F(ti)＝Fi。对于任意给定ti＝ih(i＝1，2,3…)，推力积分方程组为
式中，用F(τ)代替f(τ)是表示F(τ)为推力f(τ)的估计值。
将区间[0,ti]划分为i等分(i＝1,2,3,…)，节点为τj＝jh(j＝0,1,…,i)，被积分函数为
满足
采用复化辛普森离散化方法，将推力积分方程组离散化为推力线性方程组，需要交替使用3点和4点辛普森积分公式，分以下5种情况讨论。
(1)当i＝1时，根据梯形积分公式起步计算，可得
可得
(2)当i＝2时，根据3点辛普森积分公式，可知
可得
a20F0+a21F1＝Cfθ12 (66)
(3)当i＝3时，根据4点辛普森积分公式，可知
可得
a30F0+a31F1+a32F2＝Cfθ13 (68)
(4)当i＝2k(k＝2,3,…)时，根据3点辛普森积分公式，可得
可得
(5)当i＝2k+1(k＝2,3,…)时，在2k+2个点中，前面的i＝0,1,…,2k-3,2k-2等点采用3点辛普森积分公式；最后的4个点i＝2k-2,2k-1,2k,2k+1等采用4点辛普森积分公式，可得
即有
(6)计算得到名义推力f(τ)的估计值F(τ)，设待测推力f′(τ)的估计值为F′(τ)，可得待测推力估计值为
F′(t)＝F(t)-f0(t) (80)
如果在步骤(1)至步骤(5)中，线性方程组表示为矩阵形式，则有
推力的复化辛普森计算方法，需要交替使用3点和4点辛普森公式，计算方法构造复杂，但是由于复化辛普森公式的截断误差与时间步长h4成正比，计算精度大幅提高。
如图1所示，本发明基于数字滤波器的星载微推力器推力响应时间测量方法，包括以下步骤：
S1：对扭摆推力测量系统的非零初始条件进行归零处理,得到变量代换后的推力测量系统的振动微分方程；
S2：在变量代换后的推力测量系统后串联一个数字滤波器,得到等效灵敏度高频推力测量系统；
S3：确定等效灵敏度高频推力测量系统的系统响应；
S4：从系统响应判读星载微推力器的推力响应时间。
S5：利用等效灵敏度高频推力测量系统反演计算待测推力，可以进一步确认推力响应时间。
下面结合一个具体的实施例进行详细说明。
实施例
已知扭摆推力测量系统的扭转刚度系数为k＝0.075N·m/rad，振动频率为ωd＝0.25rad/s，阻尼比为ζ＝0.2，固有频率为ωn＝0.255155rad/s(固有线频率fn＝0.040609Hz，固有周期为24.6s)，转动惯量为J＝1.152002kg·m2，力臂为Lf＝0.4m。
在扭摆转动部件中，搭载3kg推力头和3kg配重其转动惯量为2×3×0.42＝0.96kg·m2，再加上阻尼器、标定力装置、横梁等，转动部件的转动惯量为J＝1.152002kg·m2。
采用该推力测量系统，利用所提出的推力响应时间的数字滤波器测量方法，测量推力响应时间分别为25ms、50ms、100ms的情况。
初始条件为：θ0＝50μrad，
由于灵敏度为S＝Lf/k，初始条件为θ0＝50μrad和对应的推力为：这表示推力从9.375μN开始变化。
1.建立推力响应时间与推力测量系统的进入稳态时间关系
该推力测量系统进入稳态时间为
因此，该推力测量系统只能测量推力响应时间大于89.8s的情况，远不能满足推力响应时间为25ms、50ms、100ms的测量要求。
2.非零初始条件的归零处理方法
在待测推力f′(t)作用下系统响应θ′(t)是由位移传感器测量得到。此处，为了检验推力响应时间的数字滤波器测量方法，取推力函数形式为：
f′(t)＝20[1-exp(-t/τ)](单位：μN)
式中，τ为时间常数表征推力上升时间大小，推力上升时间为tf≈5τ。由f′(t)根据变量代换前振动微分方程，计算出θ′(t)作为位移传感器的测量值。
为了研究推力响应时间分别为25ms、50ms、100ms的情况，推力的时间常数分别取τ＝5ms、10ms、20ms。如图3所示，①、②、③分别为推力响应时间为25ms、50ms、100ms时，推力上升时间的变化(推力从9.375μN开始变化到29.375μN，推力增量为f′(t))。
如图4所示，①、②、③分别为推力响应时间为25ms、50ms、100ms时，原推力测量系统的系统响应测量值。
根据初始条件θ0＝50μrad和进行变量代换，变量代换后系统响应为：
如图5所示，①、②、③分别为推力响应时间为25ms、50ms、100ms时，变量代换后系统响应的变化。
变量代换后，推力测量系统的振动微分方程为
其中，f(t)＝f′(t)+f0(t)
式中，f0(t)是与非零初始条件相关的等效推力。
3.利用数字滤波器构建等效灵敏度高频推力测量系统
在变量代换后，推力测量系统的传递函数上，串联一个数字滤波器，构造等效灵敏度高频推力测量系统。
如图2所示，在变量代换后，推力测量系统的传递函数为H(s)，其输入为待求名义推力f(t)，其输出为计算得到的θ(t)；所串联的数字滤波器的传递函数为D(s)，数字滤波器的输入为θ(t)，数字滤波器的输出为θ1(t)(由θ(t)计算得到)；所构造的等效灵敏度高频推力测量系统的传递函数为H1(s)＝H(s)D(s)，等效灵敏度高频推力测量系统的输入为名义推力f(t)，输出为θ1(t)。
4.数字滤波器设计方法与系统响应计算方法
为了减小进入稳态时间提高推力响应时间的测量能力，需要增大固有振动频率，令ωn1＝Cωωn(Cω＞1)和k1＝k(灵敏度相等)，数字滤波器的输出为θ1(t)，则有：
式中，d(t)为数字滤波器的单位脉冲响应函数，由数字滤波器设计得到。
采用数字滤波器的系统响应计算方法，计算得到等效灵敏度高频推力测量系统的系统响应θ1(t)。
选取等效灵敏度高频推力测量系统的频率比为Cω＝1000，则固有频率ωn1＝255.155rad/s(固有频率由变量代换前的0.255155rad/s，提升为变量代换后的255.155rad/s，即固有周期变为24.6ms)，并且选取阻尼比ζ1＝0.7，变量代换后推力测量系统进入稳态时间为：
如图6所示，①、②、③分别为数字滤波器的系统响应(由图5数据计算得到)。由图6可判读所施加推力的推力响应时间分别为25ms、50ms、100ms，并且推力响应时间25ms是刚好能够判读的情况。
5.利用等效灵敏度高频推力测量系统反演计算推力方法
已知数字滤波器的输出为θ1(t)，设名义推力f(τ)，等效灵敏度高频推力测量系统的推力积分方程为
根据推力积分方程，利用θ1(t)可反演计算名义推力f(τ)。
根据所构造的推力积分方程，采用数值积分离散化方法建立推力线性方程组，求解名义推力f(τ)的估计值F(τ)，设待测推力f′(τ)的估计值为F′(τ)，可得待测推力估计值为
F′(t)＝F(t)-f0(t)
待测推力f′(τ)估计值F′(τ)的误差为
式中，推力函数形式为
f′(t)＝20[1-exp(-t/τ)](单位：μN)
式中，推力的时间常数分布取τ＝5ms、10ms、20ms。
如图7所示，为推力上升时间25ms下反演计算推力和误差。从图7可清晰确认推力时间响应为25ms，并且20ms以后推力误差小于0.5％。
如图8所示，为推力上升时间50ms下反演计算推力和误差。从图8可清晰确认推力时间响应为50ms，并且20ms以后推力误差小于0.5％。
如图9所示，为推力上升时间100ms下反演计算推力和误差。从图9可清晰确认推力时间响应为100ms，并且20ms以后推力误差小于0.5％。
6.推力下降时间测量的情况
下面举例说明在推力逐渐下降时，推力下降时间测量情况。
取推力函数形式为：f′(t)＝-20[1-exp(-t/τ)](单位：μN)
式中，当推力的时间常数取τ＝5ms、10ms、20ms时，推力下降时间分别为τf＝25ms、50ms、100ms。
将上述推力代入以下扭摆推力测量系统的振动微分方程
初始条件为：θ0＝150μrad，
式中，阻尼比为ζ＝0.2，固有频率为ωn＝0.255155rad/s，转动惯量为J＝1.152002kg·m2，力臂为Lf＝0.4m。采用常微分方程数值计算方法，求得θ′(t)作为系统响应的测量值。
初始条件θ0＝150μrad和表示推力从28.125μN开始逐渐下降，减小的增量为20μN(推力从28.125μN开始减小到8.125μN)。
如图10所示，为数字滤波器的系统响应变化(由θ′(t)计算得到)。从图10可判读推力下降时间分别为25ms、50ms、100ms，并且推力下降时间25ms是刚好能够判读的情况。
如图11所示，为反演计算的推力(由θ1(t)计算得到)。从图11可清晰确认推力下降时间分别为25ms、50ms、100ms。
如图12所示，为反演计算的推力误差。在上述三种情况下，在20ms以后，推力误差都小于1％。
综上，本发明的主要技术特征和优点为：
(1)针对目前不了解推力响应时间与推力测量系统的系统响应之间内在联系的现状，本发明建立了推力响应时间与推力测量系统的进入稳态时间关系。首先，指明了推力响应时间测量中，推力从一个稳定推力开始变化，意味着推力测量系统的初始条件不为零；其次，推力测量系统的进入稳态时间是推力测量系统的固有属性，推力测量系统的固有频率越大则进入稳态时间越短，推力响应时间只有大于推力测量系统的进入稳态时间条件下，才能从系统响应曲线判读推力响应时间；最后，指明了10ms级推力响应时间测量，要求推力测量系统的固有频率为100Hz以上，这是目前任何(扭摆、吊摆、倒摆、弯振)推力测量系统在物理结构上都难以实现的，因此，需要另辟蹊径。
(2)针对目前快变推力响应时间的动态推力计算方法不成熟和不准确的现状，本发明通过固有频率低的推力测量系统上串联数字滤波器，构建等效灵敏度高频推力测量系统，采用数字滤波器设计方法，使得等效灵敏度高频推力测量系统的固有频率大幅提高并极大缩短了进入稳态时间，实现了快变推力响应时间的测量。首先，根据等效灵敏度高频推力测量系统的系统响应曲线判读推力响应时间；其次，通过推力反演计算由推力曲线确认推力响应时间，使得推力响应时间评估结论更加准确和可靠。
(3)星载微推力器推力响应时间的数字滤波器测量方法，所构建的等效灵敏度高频推力测量系统的固有频率提升500-1500倍，阻尼比调节范围由0.1-0.3变为0.7-0.9，大幅减小了进入稳态时间，实现了10ms级推力响应时间的测量。
以上所述仅为本发明的具体实施方式而已，并不用来限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种基于数字滤波器的星载微推力器推力响应时间测量方法，其特征在于，包括以下步骤：
S1：对扭摆推力测量系统的非零初始条件进行归零处理,得到变量代换后的推力测量系统的振动微分方程；
S2：在变量代换后的推力测量系统后串联一个数字滤波器,得到等效灵敏度高频推力测量系统；
S3：确定等效灵敏度高频推力测量系统的系统响应；
S4：从系统响应判读星载微推力器的推力响应时间，当系统响应与稳态系统响应比较相差1％时所对应的时间为推力响应时间；
S5：利用等效灵敏度高频推力测量系统反演计算待测推力，进一步确认推力响应时间，
等效灵敏度高频推力测量系统的单位脉冲响应函数h1(t)为
已知数字滤波器的输出为θ1(t)，利用等效灵敏度高频推力测量系统的推力积分方程：计算名义推力f(τ)，
其中，为扭转刚度系数，J为转动惯量，Lf为力臂，ωn为固有振动频率，ωn1为等效灵敏度高频推力测量系统的固有频率，ζ1为等效灵敏度高频推力测量系统的阻尼比，为改进后振动频率。
2.根据权利要求1所述的星载微推力器推力响应时间测量方法，其特征在于，步骤S1中，在待测推力f′(t)作用下系统响应θ′(t)由位移传感器测量得到，根据非零初始条件进行变量代换，令
得到变量代换后的推力测量系统的振动微分方程为：
f(t)＝f′(t)+f0(t)
其中，f0(t)是与非零初始条件相关的等效推力，f′(t)为待测推力，θ′(t)为该推力作用下由位移传感器测量得到的系统响应，ζ为阻尼比，ωn为固有振动频率，J为转动惯量，Lf为力臂；θ0表示初始扭转角，θ0和为常数，θ(t)为变量代换后系统响应，满足初始条件θ(0)＝0和
3.根据权利要求1所述的星载微推力器推力响应时间测量方法，其特征在于，步骤S2中，数字滤波器的传递函数为
其中，ωn为固有振动频率，ζ为阻尼比，ωn1为等效灵敏度高频推力测量系统的固有频率，ωn1＝Cωωn，Cω为频率比，ζ1为等效灵敏度高频推力测量系统的阻尼比，为改进后振动频率。
4.根据权利要求3所述的星载微推力器推力响应时间测量方法，其特征在于，Cω>1。
5.根据权利要求3所述的星载微推力器推力响应时间测量方法，其特征在于，ζ1＝0.7～0.9。
6.根据权利要求1所述的星载微推力器推力响应时间测量方法，其特征在于，步骤S3中，数字滤波器的输入为θ(τ)，通过位移传感器测量值θ′(τ)计算得到θ(τ)，并且计算数字滤波器输出的系统响应θ1(t)为：其中,τ为积分变量。
7.根据权利要求1所述的星载微推力器推力响应时间测量方法，其特征在于，待测推力f′(t)的估计值F′(t)为：F′(t)＝F(t)-f0(t),其中，f0(t)是与非零初始条件相关的等效推力，F(t)是名义推力f(t)的估计值，ζ为阻尼比，ωn为固有振动频率，J为转动惯量，Lf为力臂，θ0表示初始扭转角，θ0和为常数。
8.根据权利要求1-7任一项所述的星载微推力器推力响应时间测量方法，其特征在于，推力响应时间包括推力上升时间和推力下降时间。
CN202111352455.6A 2021-11-16 2021-11-16 基于数字滤波器的星载微推力器推力响应时间测量方法 Active CN114218750B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202111352455.6A CN114218750B (zh) 2021-11-16 2021-11-16 基于数字滤波器的星载微推力器推力响应时间测量方法
US17/964,466 US20230150696A1 (en) 2021-11-16 2022-10-12 Digital filter based method for measuring thrust response time of satellite-borne micro-thruster

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111352455.6A CN114218750B (zh) 2021-11-16 2021-11-16 基于数字滤波器的星载微推力器推力响应时间测量方法

Publications (2)

Publication Number Publication Date
CN114218750A CN114218750A (zh) 2022-03-22
CN114218750B true CN114218750B (zh) 2024-06-25

Family

ID=80697232

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111352455.6A Active CN114218750B (zh) 2021-11-16 2021-11-16 基于数字滤波器的星载微推力器推力响应时间测量方法

Country Status (2)

Country Link
US (1) US20230150696A1 (zh)
CN (1) CN114218750B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116773138B (zh) * 2023-08-23 2023-12-19 国科大杭州高等研究院 冷气微推力响应时间测量系统及方法
CN117288377B (zh) * 2023-09-05 2024-04-30 国科大杭州高等研究院 扭摆式微推力测量在线标定装置
CN116930668B (zh) * 2023-09-15 2023-12-05 国科大杭州高等研究院 电推力器响应时间测量的检测系统与运行方法
CN117346926B (zh) * 2023-12-06 2024-04-09 国科大杭州高等研究院 微推力测量方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE19829289C2 (de) * 1998-06-30 2001-12-06 Siemens Ag Verfahren zur Berechnung der Koeffizienten eines nichtrekursiven digitalen Filters
US9115662B1 (en) * 2009-07-10 2015-08-25 The Boeing Company Health-adaptive reaction control system
CN102689685B (zh) * 2012-06-01 2014-08-06 哈尔滨工程大学 基于在线可用功率的动力定位船推进器负载限制方法
CN104990666B (zh) * 2015-05-27 2020-03-20 中国人民解放军战略支援部队航天工程大学 一种基于比例回归法的二阶振动测量系统的系统参数标定方法
CN108254145B (zh) * 2017-12-29 2020-06-26 苏州东菱智能减振降噪技术有限公司 一种实现多振动台同步振动的控制方法
CN108829946A (zh) * 2018-05-29 2018-11-16 中国人民解放军战略支援部队航天工程大学 一种基于动态补偿技术的推力计算方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
一种基于动态补偿技术的微小稳态推力还原方法;周伟静 等;航空学报;摘要，第0-5节 *

Also Published As

Publication number Publication date
CN114218750A (zh) 2022-03-22
US20230150696A1 (en) 2023-05-18

Similar Documents

Publication Publication Date Title
CN114218750B (zh) 基于数字滤波器的星载微推力器推力响应时间测量方法
CN102620729B (zh) 一种机抖激光陀螺惯性测量单元数字滤波器设计方法
CN109960831B (zh) 用于扭摆系统的微推力平滑降噪优化还原方法
CN105974415A (zh) 一种机载sar方位空变运动误差的高精度补偿方法
CN105387859A (zh) Mems传感器组合温度漂移误差补偿方法
EP4375122A1 (en) Temperature compensation method and apparatus based on direct current charging base
CN106052668B (zh) 一种大量程硅微陀螺仪非线性数字补偿方法
Uchiyama et al. Measurement of instantaneous flow rate through estimation of velocity profiles
CN109426143B (zh) 负载转矩估算方法、系统、机电控制系统、方法及电机
JP2867477B2 (ja) オンライン設備の寿命予測方法
CN104655876A (zh) 一种恒加速度和振动复合输入情况下的线加速度计校准方法
CN116026328A (zh) 微惯导的零偏滞回效应补偿模型的构建方法和补偿方法
CN109298376B (zh) 一种基于标准电能表组的电能量值传递的方法及系统
Kumar et al. A method to determine wall shear stress from mean profiles in turbulent boundary layers
CN110967045B (zh) 光纤传感器的温度效应误差补偿系统与设计方法
CN210981427U (zh) 一种燃油液位检测装置
CN111929585B (zh) 电池电荷状态计算装置、方法、服务器及介质
EP0429472A1 (de) Verfahren zur verbesserung des nordungsergebnisses.
JP2987635B2 (ja) ディジタル信号処理方法及びその装置
CN112964336B (zh) 空气流量标定设备及其标定方法
US20240192250A1 (en) Method for improving calibration efficiency of accelerometer
CN111504551B (zh) 一种基于最小二乘复指数法的应变式力矩仪带宽扩展方法
JP2003515117A (ja) 慣性測定システム
CN112231890B (zh) 基于扭摆测量系统的高平稳电推力器的推力测评方法
CN118258532A (zh) 一种车辆运行阻力系数的估算方法、装置及车辆

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant