DIGITAL FILTER BASED METHOD FOR MEASURING THRUST RESPONSE TIME OF SATELLITE-BORNE MICRO-THRUSTER

Abstract

The present disclosure belongs to the technical field of space satellite propulsion, and particularly relates to a digital filter based method for measuring thrust response time of a satellite-borne micro-thruster. The method for measuring thrust response time in the present disclosure includes the following steps: S1: zeroing non-zero initial conditions of a torsional pendulum thrust measurement system to obtain an oscillating differential equation for a thrust measurement system after variable substitution; S2: connecting the digital filter in series behind the thrust measurement system after variable substitution to obtain an equivalent-sensitivity high-frequency thrust measurement system; S3: determining a system response of the equivalent-sensitivity high-frequency thrust measurement system; S4: determining and reading thrust response time of the satellite-borne micro-thruster from the system response; and S5: computing thrust to be measured by means of the system response inversely, and further confirming the thrust response time.

Claims

1. A digital filter based method for measuring thrust response time of a satellite-borne micro-thruster, comprising the following steps: S1: zeroing non-zero initial conditions of a torsional pendulum thrust measurement system to obtain an oscillating differential equation for a thrust measurement system after variable substitution; S2: connecting the digital filter in series behind the thrust measurement system after variable substitution to obtain an equivalent-sensitivity high-frequency thrust measurement system; S3: determining a system response of the equivalent-sensitivity high-frequency thrust measurement system; S4: determining and reading thrust response time of the satellite-borne micro-thruster from the system response; and S5: computing thrust to be measured by means of an equivalent-sensitivity high-frequency thrust measurement system inversely, and further confirming the thrust response time.

2. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 1, wherein in the step S1, a system response θ′(t) is measured by a displacement sensor under the action of the thrust f′(t) to be measured, variable substitution is carried out according to the non-zero initial conditions (θ.sub.0, {dot over (θ)}.sub.0), and by making
θ(t)=θ′(t)−{dot over (θ)}.sub.0t−θ.sub.0, {dot over (θ)}(t)={dot over (θ)}′(t)−{dot over (θ)}.sub.0 and {umlaut over (θ)}(t)={umlaut over (θ)}′(t), the oscillating differential equation for the thrust measurement system after variable substitution is obtained: $\stackrel{.Math.}{\theta }\left(t\right)+2{\mathrm{\zeta \omega }}_n\dot{\theta }\left(t\right)+{\omega }_n^2\theta \left(t\right)=\frac{L_f}{J}f\left(t\right)$ $f\left(t\right)=f^{\prime }\left(t\right)+f_0\left(t\right)$ $f_0\left(t\right)=-\frac{J}{L_f}\left(2{\mathrm{\zeta \omega }}_n{\dot{\theta }}_0+{\omega }_n^2{\dot{\theta }}_{0}t+{\omega }_n^2{\theta }_0\right)$ wherein f.sub.0(t) represents equivalent thrust related to the non-zero initial conditions, f′(t) represents the thrust to be measured, θ′(t) represents the system response measured by the displacement sensor under the action of the thrust, ζ represents a damping ratio, ω.sub.n represents a natural vibration frequency, J represents moment of inertia, L.sub.f represents a moment arm, θ.sub.0 represents an initial torsion angle, θ.sub.0 and {dot over (θ)}.sub.0 are constants, θ(t) represents a system response after variable substitution, and initial conditions θ(0)=0 and {dot over (θ)}(0)=0 are satisfied.

3. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 1, wherein in the step S2, a transfer function of the digital filter is, and $$\begin{gathered} D\left(s\right)=\frac{{\omega }_{n1}^2}{{\omega }_n^2}[1-2\left({\zeta }_1{\omega }_{n1}-{\mathrm{\zeta \omega }}_n\right)\frac{s+{\zeta }_1{\omega }_{n1}}{{\left(s+{\zeta }_1{\omega }_{n1}\right)}^2+{\omega }_{d1}^2})- \\ \frac{\left({\omega }_{n1}^2-{\omega }_n^2\right)-2\left({\zeta }_1{\omega }_{n1}-\zeta {\omega }_n\right){\zeta }_1{\omega }_{n1}}{{\omega }_{d1}}.Math.\frac{{\omega }_{d1}}{{\left(s+{\zeta }_1{\omega }_{n1}\right)}^2+{\omega }_{d1}^2}] \end{gathered}$$ a unit impulse response function of the digital filter is, $$\begin{gathered} d\left(t\right)=C_{\omega }^2\left[\delta \left(t\right)-2\left({\zeta }_1{C}_{\omega }-\zeta \right){\omega }_n{e}^{-{\zeta }_1{C}_{\omega }{\omega }_{n}t}\cos \left(\sqrt{1-{\zeta }_1^2}{C}_{\omega }{\omega }_{n}t\right)- \\ \frac{\left(C_{\omega }^2-1\right)-2\left({\zeta }_1{C}_{\omega }-\zeta \right){\zeta }_1{C}_{\omega }}{\sqrt{1-{\zeta }_1^2}{C}_{\omega }}{\omega }_n{e}^{-{\zeta }_1{C}_{\omega }{\omega }_{n}t}\sin \left(\sqrt{1-{\zeta }_1^2}{C}_{\omega }{\omega }_{n}t\right)\right] \end{gathered}$$ wherein ω.sub.n represents the natural vibration frequency, ζ represents the damping ratio, ω.sub.n1 represents a natural frequency of the equivalent-sensitivity high-frequency thrust measurement system, ω.sub.n1=C.sub.ωω.sub.n, C.sub.ω represents a frequency ratio, ζ.sub.1 represents a damping ratio of the equivalent-sensitivity high-frequency thrust measurement system, and ω.sub.d1=√{square root over (1−ζ.sub.1.sup.2)}ω.sub.n1 represents an improved vibration frequency.

4. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 3, wherein C.sub.ω>1.

5. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 3, wherein ζ.sub.1=0.7˜0.9.

6. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 1, wherein in the step S3, input of the digital filter is θ(τ), θ(τ) is computed according to a measured value θ′(τ) of the displacement sensor, and a system response θ.sub.1(t) output by the digital filter is computed: θ.sub.1(t)=∫.sub.0.sup.tθ(τ)d(t−τ)dτ, wherein τ is an integral variable.

7. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 1, wherein in the step S4, under the situation that the system response is 1% greater than or less than a steady-state system response, corresponding time is the thrust response time.

8. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 1, wherein in the step S5, a unit impulse response function h.sub.1(t) of the equivalent-sensitivity high-frequency thrust measurement system is, $h_1\left(t\right)=\frac{{\omega }_{n1}^2{L}_f}{k{\omega }_{d1}}{e}^{-{\zeta }_1{\omega }_{n1}t}\sin \left({\omega }_{d1}t\right)$ what is known is that output of the digital filter is θ.sub.1(t), and nominal thrust f(τ) is computed by means of an integral thrust equation for an equivalent-sensitivity high-frequency thrust measurement system: θ.sub.1(t)=∫.sub.0.sup.tf(τ)h.sub.1(t−τ)dτ.

9. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 8, wherein an estimated value F′(t) of the thrust f′(t) to be measured is: F′(t)=F(t)−f(t)−f.sub.0(t), wherein $f_0\left(t\right)=-\frac{J}{L_f}\left(2{\mathrm{\zeta \omega }}_n{\dot{\theta }}_0+{\omega }_n^2{\dot{\theta }}_{0}t+{\omega }_n^2{\theta }_0\right),$ f.sub.0(t) represents an equivalent thrust related to the non-zero initial conditions, F(t) represents an estimated value of the nominal thrust f(t), ζ represents a damping ratio, ω.sub.n represents a natural vibration frequency, J represents moment of inertia, L.sub.f represents a moment arm, θ.sub.0 represents an initial torsion angle, and θ.sub.0 and {dot over (θ)}.sub.0 are constants.

10. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 1, wherein the thrust response time comprises thrust increasing time and thrust decreasing time.

11. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 2, wherein the thrust response time comprises thrust increasing time and thrust decreasing time.

12. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 3, wherein the thrust response time comprises thrust increasing time and thrust decreasing time.

13. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 4, wherein the thrust response time comprises thrust increasing time and thrust decreasing time.

14. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 5, wherein the thrust response time comprises thrust increasing time and thrust decreasing time.

15. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 6, wherein the thrust response time comprises thrust increasing time and thrust decreasing time.

16. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 7, wherein the thrust response time comprises thrust increasing time and thrust decreasing time.

17. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 8, wherein the thrust response time comprises thrust increasing time and thrust decreasing time.

18. The method for measuring thrust response time of a satellite-borne micro-thruster according to claim 9, wherein the thrust response time comprises thrust increasing time and thrust decreasing time.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0027] FIG. 1 is a flow diagram of a method for measuring thrust response time in the present disclosure.

[0028] FIG. 2 is a schematic diagram of an equivalent-sensitivity high-frequency thrust measurement system in the present disclosure.

[0029] FIG. 3 is a variation graph of thrust increasing time.

[0030] FIG. 4 is a graph of measured values of system responses of an original thrust measurement system.

[0031] FIG. 5 is a variation graph of system responses after variable substitution.

[0032] FIG. 6 is a variation graph of system responses of a digital filter.

[0033] FIG. 7 is a graph of thrust and errors inversely computed under the situation that thrust increasing time is 25 ms.

[0034] FIG. 8 is a graph of thrust and errors inversely computed under the situation that thrust increasing time is 50 ms.

[0035] FIG. 9 is a graph of thrust and errors inversely computed under the situation that thrust increasing time is 100 ms.

[0036] FIG. 10 is a variation graph of system responses of a digital filter under the situation that thrust is gradually decreased.

[0037] FIG. 11 is a graph of thrust inversely computed under the situation that thrust is gradually decreased.

[0038] FIG. 12 is a graph of thrust errors inversely computed under the situation that thrust is gradually decreased.

DETAILED DESCRIPTION OF THE EMBODIMENTS

[0039] The specific implementations of the present disclosure will be described in detail below with reference to the accompanying drawings.

[0040] As shown in FIG. 1, a flow diagram of a method for measuring thrust response time in the present disclosure is provided, and an implementation process will be described in detail below.

[0041] 1. Establishing Relation Between Thrust Response Time and Time of Entering Steady State of Thrust Measurement System

[0042] Establishing a relation between thrust response time and time of entering a steady state of a thrust measurement system is a necessary condition for the thrust measurement system to measure given thrust response time.

[0043] 1) Thrust Response Time and Initial Conditions of Oscillating Differential Equation

[0044] Thrust increasing time and thrust decreasing time are required to be considered for thrust response time of a micro-thruster. Under the situation that small steady thrust is increased into large steady thrust, thrust response time may be represented by thrust increasing time, and under the situation that large steady thrust is decreased into small steady thrust, thrust response time may be represented by thrust decreasing time. The effects of initial conditions are certainly required to be considered for researching thrust response time.

[0045] An oscillating differential equation for a torsional pendulum thrust measurement system is

[00006]$\begin{array}{cc}{\stackrel{.Math.}{\theta }}^{\prime }\left(t\right)+2{\mathrm{\zeta \omega }}_n{\dot{\theta }}^{\prime }\left(t\right)+{\omega }_n^2{\theta }^{\prime }\left(t\right)=\frac{L_f}{J}{f}^{\prime }\left(t\right),& \left(1\right)\end{array}$

[0046] where f′(t) represents thrust to be measured, θ′(t) represents a system response (measured by a displacement sensor) under the action of thrust, ζ represents a damping ratio, ω.sub.n represents a natural vibration frequency, J represents moment of inertia, and L.sub.f represents a moment arm.

[0047] Under the situation that thrust is varied from one steady thrust to another steady thrust, the initial conditions are not zero but are

θ′(0)=θ.sub.0 and {dot over (θ)}′(0)={dot over (θ)}.sub.0  (2),

[0048] where under the situation that thrust is varied from non-zero steady thrust, it may be approximately considered that {dot over (θ)}′(0)≈0, and under the situation that thrust is varied from zero thrust, it may be approximately considered that θ′(0)≈0 and {dot over (θ)}′(0)≈0.

[0049] Obviously, the initial conditions of the oscillating differential equation for the thrust measurement system reflect from what initial thrust thrust is varied, for example, from a non-zero steady thrust or from zero thrust. Therefore, the initial conditions of the oscillating differential equation for the thrust measurement system are certainly required to be clearly defined for measuring the thrust response time.

[0050] 2) Establishing Relation Between Thrust Response Time and Time of Entering Steady State

[0051] An oscillating differential equation for a torsional pendulum thrust measurement system is

[00007]$\begin{array}{cc}{\stackrel{.Math.}{\theta }}^{\prime }\left(t\right)+2{\mathrm{\zeta \omega }}_n{\dot{\theta }}^{\prime }\left(t\right)+{\omega }_n^2{\theta }^{\prime }\left(t\right)=\frac{L_f}{J}{f}^{\prime }\left(t\right).& \left(3\right)\end{array}$

[0052] Under the situation that the initial conditions are not zero but are θ′(0)=θ.sub.0 and {dot over (θ)}′(0)={dot over (θ)}.sub.0, assuming that laplace transform of thrust is L[f′(t)]=F′(s), laplace transform of a system response and a derivative is

L[θ′(t)]=Θ′(s)  (4)

L[{dot over (θ)}′(t)]=sΘ′(s)−θ.sub.0  (5) and

L[{umlaut over (θ)}′(t)]=s.sup.2Θ′(s)−sθ.sub.0−{dot over (θ)}.sub.0)  (6),

[0053] laplace transform is carried out on two sides of the oscillating equation, and finishing is carried out to obtain

[00008]$\begin{array}{cc}{\Theta }^{\prime }-\left(s\right)=\frac{L_f}{J}{F}^{\prime }\left(s\right)\frac{1}{{\left(s+{\mathrm{\zeta \omega }}_n\right)}^2+{\omega }_d^2}+\frac{\left(s+{\mathrm{\zeta \omega }}_n\right){\theta }_0}{{\left(s+{\mathrm{\zeta \omega }}_n\right)}^2+{\omega }_d^2}+\frac{{\mathrm{\zeta \omega }}_n{\theta }_0+{\dot{\theta }}_0}{{\left(s+{\mathrm{\zeta \omega }}_n\right)}^2+{\omega }_d^2},& \left(7\right)\end{array}$

[0054] inverse laplace transform is carried out on two sides to obtain

[00009]$\begin{array}{cc}{\theta }^{\prime }\left(t\right)=L^{-1}\left[{\Theta }^{\prime }\left(s\right)\right]=\frac{L_f}{J{\omega }_d}{\int }_0^t{f}^{\prime }\left(\tau \right)e^{-\zeta {\omega }_n\left(t-\tau \right)}\sin {\omega }_d\left(t-\tau \right)d\tau +{\theta }_0{e}^{-\zeta {\omega }_{n}t}\cos {\omega }_{d}t+\frac{{\mathrm{\zeta \omega }}_n{\theta }_0+{\dot{\theta }}_0}{{\omega }_d}{e}^{-\zeta {\omega }_{n}t}\sin {\omega }_{d}t,\mathrm{and}& \left(8\right)\end{array}$$\begin{array}{cc}{\theta }^{\prime }\left(t\right)={\theta }_1^{\prime }\left(t\right)+{\theta }_2^{\prime }\left(t\right)& \left(9\right)\end{array}$$\begin{array}{cc}{\theta }_1^{\prime }\left(t\right)=\frac{L_f}{J{\omega }_d}{\int }_0^{t}f\left(\tau \right)e^{-\zeta {\omega }_n\left(t-\tau \right)}\sin {\omega }_d\left(t-\tau \right)d\tau & \left(10\right)\end{array}$$\begin{array}{cc}{\theta }_2^{\prime }\left(t\right)=\sqrt{{\theta }_0^2+{\left[\frac{{\mathrm{\zeta \omega }}_n{\theta }_0+{\stackrel{.}{\theta }}_0}{{\omega }_d}\right]}^2}{e}^{-\zeta {\omega }_{n}t}\sin \left({\omega }_{d}t+\beta \right)\mathrm{and}& \left(11\right)\end{array}$$\begin{array}{cc}\tan \beta =\frac{{\omega }_d{\theta }_0}{{\mathrm{\zeta \omega }}_n{\theta }_0+{\dot{\theta }}_0}& \left(12\right)\end{array}$

are made,

[0055] where θ′.sub.1(t) includes a transient system response and a steady-state system response. For example, by substituting periodic thrust f′(t)=A.sub.f cos ωt (A.sub.f is a constant representing an amplitude of the periodic thrust, and ω represents a frequency of the periodic thrust) in θ′.sub.1(t), a system response of a thrust measurement system is

[00010]$\begin{array}{cc}{\theta }_1^{\prime }\left(t\right)=\frac{C}{\sqrt{1-{\zeta }^2}}{e}^{-\zeta {\omega }_{n}t}\cos \left({\omega }_{d}t-{\phi }_1\right)+C\cos \left(\omega t-{\phi }_2\right),& \left(13\right)\end{array}$$\begin{array}{cc}\mathrm{where}C=A_f\frac{L_f}{k}.Math.\frac{1}{\sqrt{{\left(1-R_{\omega }^2\right)}^2+{\left(2\zeta R_{\omega }\right)}^2}}\mathrm{and}{R}_{\omega }=\frac{\omega }{{\omega }_n},& \left(14\right)\end{array}$

[0056] where S=L.sub.f/k represents sensitivity of the thrust measurement system, k=Jω.sub.n.sup.2 represents a torsional stiffness coefficient, and R.sub.ω=ω/ω.sub.n represents a frequency ratio. Obviously, a first item of a system response in θ′.sub.1(t) is a transient system response, which represents a variation of a transient process with time, and a vibration frequency of the transient process is a vibration frequency of a thrust measurement system. A second term is a steady-state system response, which represents a variation of a steady-state process with time, and a vibration frequency of a steady-state process is a vibration frequency of periodic thrust. θ′.sub.2(t) is also a transient system response, a vibration frequency of which is a vibration frequency of a thrust measurement system.

[0057] In a word, an amplitude of a transient system response is gradually attenuated, an attenuation term is e.sup.−ζω.sup.n.sup.t, and assuming that corresponding time (corresponding time when a transient amplitude is attenuated to 1%) when the attenuation term is e.sup.−ζω.sup.n.sup.t=0.01 is defined as transient time or time t.sub.s of entering a steady state,

[00011]$\begin{array}{cc}{t}_s\approx \frac{4.605170}{{\mathrm{\zeta \omega }}_n}\approx \frac{4.605170}{2\pi \zeta }.Math.\frac{2\pi }{{\omega }_n}\approx 0.73\frac{T_n}{\zeta },& \left(15\right)\end{array}$

[0058] where T.sub.n=2π/ω.sub.n is a natural period of a thrust measurement system.

[0059] Time of entering a steady state or transient time of a thrust measurement system are inherent properties of the thrust measurement system, and the following conclusions are drawn:

[0060] (1) the greater a natural frequency of a thrust measurement system, the less a period, and the shorter time of entering a steady state;

[0061] (2) the less a damping ratio of a thrust measurement system, the longer time of entering a steady state; and

[0062] (3) only under the situation that thrust response time is longer than time of entering a steady state of a thrust measurement system, variation characteristics of thrust response time may be represented from a system response curve, or thrust response time may be determined and read from a system response curve.

[0063] For example, within a common damping ratio range 0.1≤ζ≤0.3, time of entering a steady state of a thrust measurement system is 3 times to 7 times a natural period, and assuming that thrust increasing time is t.sub.s, and a natural period of a thrust measurement system is T.sub.n,

t.sub.s≥7T.sub.n  (16).

[0064] Under the situation that thrust increasing time is required to be t.sub.s≤50 ms, a natural period of a thrust measurement system is required to be T.sub.n≤7.1 ms, that is, a linear frequency is f.sub.n≥140 Hz, which is difficult to realize by any thrust measurement system of any structure. It is indicated that in order to measure response time of rapid-varying thrust, it is not enough to increase a natural frequency of a thrust measurement system only from the design of a physical structure, and another way should be explored.

[0065] 2. Method for Zeroing Non-Zero Initial Conditions

[0066] During thrust response time measurement, initial conditions are not zero. However, during subsequent design of a digital filter and construction of an equivalent-sensitivity high-frequency thrust measurement system, since a model is established by means of a transfer function, non-zero initial conditions are certainly required to be zeroed.

[0067] An oscillating differential equation for a torsional pendulum thrust measurement system is

[00012]$\begin{array}{cc}{\stackrel{.Math.}{\theta }}^{\prime }\left(t\right)+2{\mathrm{\zeta \omega }}_n{\dot{\theta }}^{\prime }\left(t\right)+{\omega }_n^2{\theta }^{\prime }\left(t\right)=\frac{L_f}{J}{f}^{\prime }\left(t\right),& \left(17\right)\end{array}$

and

[0068] initial conditions are not zero but are

θ′(0)=θ.sub.0 and {dot over (θ)}′(0)={dot over (θ)}.sub.0  (18), where

[0069] θ.sub.0 and {dot over (θ)}.sub.0 are constants.

[0070] In order to guarantee that the initial conditions are zero, variable substitution is carried out, that is,

θ(t)=θ′(t)−{dot over (θ)}.sub.0t−θ.sub.0  (19)

{dot over (θ)}(t)={dot over (θ)}′(t)−{dot over (θ)}.sub.0  (20) and

{umlaut over (θ)}(t)={umlaut over (θ)}′(t)  (21) are made,

[0071] an oscillating equation is substituted, and finishing is carried out to obtain

[00013]$\begin{array}{cc}\stackrel{.Math.}{\theta }\left(t\right)+2{\mathrm{\zeta \omega }}_n\dot{\theta }\left(t\right)+{\omega }_n^2\theta \left(t\right)=\frac{L_f}{J}{f}^{\prime }\left(t\right)-2{\mathrm{\zeta \omega }}_n{\dot{\theta }}_0-{\omega }_n^2{\dot{\theta }}_{0}t-{\omega }_n^2{\theta }_0,& \left(22\right)\end{array}$

[0072] the initial conditions are substituted to obtain

θ(0)=0 and {dot over (θ)}(0)=0  (23), and

[0073] an oscillating differential equation for a thrust measurement system after variable substitution is rewritten as

[00014]$\begin{array}{cc}\stackrel{.Math.}{\theta }\left(t\right)+2{\mathrm{\zeta \omega }}_n\dot{\theta }\left(t\right)+{\omega }_n^2\theta \left(t\right)=\frac{L_f}{J}f\left(t\right),& \left(24\right)\end{array}$$\begin{array}{cc}f\left(t\right)=f^{\prime }\left(t\right)+f_0\left(t\right)\mathrm{and}& \left(25\right)\end{array}$$\begin{array}{cc}{f}_0\left(t\right)=-\frac{J}{L_f}\left(2{\mathrm{\zeta \omega }}_n{\dot{\theta }}_0+{\omega }_n^2{\dot{\theta }}_{0}t+{\omega }_n^2{\theta }_0\right),& \left(26\right)\end{array}$

[0074] where f.sub.0(t) represents equivalent thrust related to non-zero initial conditions.

[0075] Under the situation that thrust is varied from one steady thrust to another steady thrust, θ′(0)=θ.sub.0 and {dot over (θ)}′(0)≈0, and then

θ(t)=θ′(t)−θ.sub.0  (27).

[0076] Notes: a system response θ(t) after substitution is the same as a system response θ′(t) before substitution, a transient system response is the same as a steady-state system response, and only a system response curve moves up and down by a certain distance.

[0077] A method for zeroing non-zero initial conditions is as follows:

[0078] (1) A system response θ′(t) is measured by a displacement sensor under the action of thrust f′(t) to be measured, and variable substitution is carried out according to non-zero initial conditions (θ.sub.0,{dot over (θ)}.sub.0) to obtain

θ(t)=θ′(t)−{dot over (θ)}.sub.0t−θ.sub.0  (28),

[0079] where θ(t) represents a system response after variable substitution, and zero initial conditions, that is, θ(0)=0 and {dot over (θ)}(0)=0 are satisfied.

[0080] (2) Under the situation that thrust is varied from one steady thrust to another steady thrust (θ′(0)=θ.sub.0 and {dot over (θ)}′(0)≈0), a system response θ(t) after variable substitution and a system response θ′(t) before substitution have the same time of entering a steady state such that thrust response time may be determined and read from a system response θ(t) curve after substitution.

[0081] (3) According to a system response θ(t) after variable substitution, nominal thrust f(t) may be computed by means of an oscillating differential equation after variable substitution, such that thrust f′(t) to be measured is obtained:

[00015]$\begin{array}{cc}{f}^{\prime }\left(t\right)=f\left(t\right)-f_0\left(t\right),\mathrm{and}& \left(29\right)\end{array}$$\begin{array}{cc}{f}_0\left(t\right)=-\frac{J}{L_f}\left(2{\mathrm{\zeta \omega }}_n{\dot{\theta }}_0+{\omega }_n^2{\dot{\theta }}_{0}t+{\omega }_n^2{\theta }_0\right),& \left(30\right)\end{array}$

[0082] where f.sub.0(t) represents equivalent thrust related to non-zero initial conditions.

[0083] 3. Constructing Equivalent-Sensitivity High-Frequency Thrust Measurement System by Means of Digital Filter

[0084] In order to solve the problem that a natural frequency of a thrust measurement system is too low to satisfy the requirement of measuring response time of rapid-varying thrust, an equivalent-sensitivity high-frequency thrust measurement system is constructed by means of a digital filter, such that the natural frequency is greatly increased, and the requirement of measuring a response of rapid-varying thrust is satisfied.

[0085] As shown in FIG. 2, a schematic diagram of an equivalent-sensitivity high-frequency thrust measurement system in the present disclosure is shown. A method for constructing an equivalent-sensitivity high-frequency thrust measurement system by means of a digital filter is as follows:

[0086] (1) a digital filter is connected in series on a transfer function of a thrust measurement system after variable substitution to construct an equivalent-sensitivity high-frequency thrust measurement system;

[0087] (2) the constructed equivalent-sensitivity high-frequency thrust measurement system and a thrust measurement system before variable substitution have the same sensitivity, but a vibration frequency and a damping ratio may be adjusted through the design of the digital filter; and

[0088] (3) according to the measurement requirements of thrust increasing time and thrust decreasing time, the thrust increasing time and the thrust decreasing time are measured by adjusting the vibration frequency and the damping ratio.

[0089] 4. Method for Designing Digital Filter and Method for Computing System Response

[0090] 1) Method for Designing Digital Filter

[0091] Assuming that input of a thrust measurement system after variable substitution is nominal thrust f(t) (laplace transform is F(s)), and a system response thereof is θ(t) (computed according to a measured value θ′(t) of a displacement sensor, and laplace transform is Θ(s)). A transfer function of a digital filter is D(s), input of the digital filter is θ(t), and output of the digital filter is θ.sub.1(t) computed by means of θ(t), and laplace transform is Θ.sub.1(s)). A transfer function of a constructed equivalent-sensitivity high-frequency thrust measurement system is H.sub.1(s)=H(s)D(s), input of the equivalent-sensitivity high-frequency thrust measurement system is nominal thrust f(t), and output is θ.sub.1(t).

[0092] A transfer function of a torsional pendulum thrust measurement system is

[00016]$\begin{array}{cc}H\left(s\right)=\frac{L_f}{J}.Math.\frac{1}{s^2+2{\mathrm{\zeta \omega }}_{n}s+{\omega }_n^2}=\frac{L_f}{k}.Math.\frac{{\omega }_n^2}{s^2+2{\mathrm{\zeta \omega }}_{n}s+{\omega }_n^2},& \left(31\right)\end{array}$

[0093] where represents L.sub.f a moment arm, J represents moment of inertia, ζ represents a damping ratio, ω.sub.n represents a natural vibration frequency, ω.sub.d=√{square root over (1−ζ.sup.2)}ω.sub.n represents a vibration frequency, and k=Jω.sub.n.sup.2 represents a torsional stiffness coefficient.

[0094] A method for designing of a digital filter is to connect a digital filter in series to an output end of a thrust measurement system after variable substitution to form an improved torsional pendulum thrust measurement system. A transfer function of the improved torsional pendulum thrust measurement system is

[00017]$\begin{array}{cc}{H}_1\left(s\right)=\frac{L_f}{k_1}.Math.\frac{{\omega }_{n1}^2}{s^2+2{\zeta }_1{\omega }_{n1}s+{\omega }_{n1}^2},& \left(32\right)\end{array}$

[0095] where a natural vibration frequency of the improved torsional pendulum thrust measurement system is ω.sub.n1, a damping ratio is ζ.sub.1, and a torsional stiffness coefficient is k.sub.1.

[0096] Assuming that a transfer function of a digital filter is D(s),

[00018]$\begin{array}{cc}D\left(s\right)=\frac{H_1\left(s\right)}{H\left(s\right)}=\frac{\frac{L_f}{k_1}.Math.\frac{{\omega }_{n1}^2}{s^2+2{\zeta }_1{\omega }_{n1}s+{\omega }_{n1}^2}}{\frac{L_f}{k}.Math.\frac{{\omega }_n^2}{s^2+2{\mathrm{\zeta \omega }}_{n}s+{\omega }_n^2}}=\left(\frac{k}{k_1}\right)\left(\frac{{\omega }_{n1}^2}{{\omega }_n^2}.Math.\frac{s^2+2{\mathrm{\zeta \omega }}_{n}s+{\omega }_n^2}{s^2+2{\zeta }_1{\omega }_{n1}s+{\omega }_{n1}^2}\right).& \left(33\right)\end{array}$

[0097] In order to shorten time of entering a steady state and improve a capability of measuring thrust response time, a natural vibration frequency is required to be increased, and by making ω.sub.n1=C.sub.ωω.sub.n (C.sub.ω>1) and k.sub.1=k (having the same sensitivity),

[00019]$\begin{array}{cc}D\left(s\right)=\frac{{\omega }_{n1}^2}{{\omega }_n^2}.Math.\frac{s^2+2{\mathrm{\zeta \omega }}_{n}s+{\omega }_n^2}{s^2+2{\zeta }_1{\omega }_{n1}s+{\omega }_{n1}^2},\mathrm{and}& \left(34\right)\end{array}$$\begin{array}{cc}{H}_1\left(s\right)=\frac{L_f}{k}.Math.\frac{{\omega }_{n1}^2}{s^2+2{\zeta }_1{\omega }_{n1}s+{\omega }_{n1}^2}.& \left(35\right)\end{array}$

[0098] A transfer function of a digital filter is

[00020]$\begin{array}{cc}\begin{array}{c}D\left(s\right)=\frac{{\omega }_{n1}^2}{{\omega }_n^2}\left[1-\frac{2\left({\zeta }_1{\omega }_{n1}-{\mathrm{\zeta \omega }}_n\right)s}{s^2+2{\zeta }_1{\omega }_{n1}s+{\omega }_{n1}^2}-\frac{\left({\omega }_{n1}^2-{\omega }_n^2\right)}{s^2+2{\zeta }_1{\omega }_{n1}s+{\omega }_{n1}^2}\right]\\ =\frac{{\omega }_{n1}^2}{{\omega }_n^2}\left[1-\frac{2\left({\zeta }_1{\omega }_{n1}-{\mathrm{\zeta \omega }}_n\right)s}{{\left(s+{\zeta }_1{\omega }_{n1}\right)}^2+{\omega }_{d1}^2}-\frac{\left({\omega }_{n1}^2-{\omega }_n^2\right)}{{\left(s+{\zeta }_1{\omega }_{n1}\right)}^2+{\omega }_{d1}^2}\right]\end{array},& \left(36\right)\end{array}$

[0099] where ω.sub.d1=√{square root over (1−ζ.sub.1.sup.2)}ω.sub.n1 represents an improved vibration frequency. Through further simplification, a transfer function of a digital filter is

[00021]$\begin{array}{cc}D\left(s\right)=\frac{{\omega }_{n1}^2}{{\omega }_n^2}[1-2\left({\zeta }_1{\omega }_{n1}-{\mathrm{\zeta \omega }}_n\right)\frac{s+{\zeta }_1{\omega }_{n1}}{{\left(s+{\zeta }_1{\omega }_{n1}\right)}^2+{\omega }_{d1}^2})-\frac{\left({\omega }_{n1}^2-{\omega }_n^2\right)-2\left({\zeta }_1{\omega }_{n1}-{\mathrm{\zeta \omega }}_n\right){\zeta }_1{\omega }_{n1}}{{\omega }_{d1}}.Math.\frac{{\omega }_{d1}}{{\left(s+{\zeta }_1{\omega }_{n1}\right)}^2+{\omega }_{d1}^2}].& \left(37\right)\end{array}$

[0100] Assuming that a unit impulse response function of a digital filter is d(t),

[00022]$\begin{array}{cc}\begin{array}{c}d\left(t\right)=L^{-1}\left[D\left(s\right)\right]\\ =\frac{{\omega }_{n1}^2}{{\omega }_n^2}[\delta \left(t\right)-2\left({\zeta }_1{\omega }_{n1}-{\mathrm{\zeta \omega }}_n\right)e^{-{\zeta }_1{\omega }_{n1}t}\cos \left({\omega }_{d1}t\right)-\\ \frac{\left({\omega }_{n1}^2-{\omega }_n^2\right)-2\left({\zeta }_1{\omega }_{n1}-{\mathrm{\zeta \omega }}_n\right){\zeta }_1{\omega }_{n1}}{{\omega }_{d1}}{e}^{-{\zeta }_1{\omega }_{n1}t}\sin \left({\omega }_{d1}t\right)]\end{array},& \left(38\right)\end{array}$

and

[0101] a unit impulse response function of a digital filter is

[00023]$\begin{array}{cc}d\left(t\right)=C_{\omega }^2\left[\delta \left(t\right)-2\left({\zeta }_1{C}_{\omega }-\zeta \right){\omega }_n{e}^{-{\zeta }_1{C}_{\omega }{\omega }_{n}t}\cos \left(\sqrt{1-{\zeta }_1^2}{C}_{\omega }{\omega }_{n}t\right)-\frac{\left(C_{\omega }^2-1\right)-2\left({\zeta }_1{C}_{\omega }-\zeta \right){\zeta }_1{C}_{\omega }}{\sqrt{1-{\zeta }_1^2}{C}_{\omega }}{\omega }_n{e}^{-{\zeta }_1{C}_{\omega }{\omega }_{n}t}\sin \left(\sqrt{1-{\zeta }_1^2}{C}_{\omega }{\omega }_{n}t\right)\right].& \left(39\right)\end{array}$

[0102] Input of the digital filter is θ(τ) (computed θ′(τ) according to a measured value of a displacement sensor), output of the digital filter is θ.sub.1(t), then

θ.sub.1(t)=∫.sub.0.sup.tθ(τ)d(t−τ)dτ  (40), and

[0103] therefore, a system response of an equivalent-sensitivity high-frequency thrust measurement system is obtained.

[0104] A method for designing a digital filter is as follows:

[0105] (1) by taking equal sensitivity and increase in natural vibration frequency as objectives, a transfer function of a digital filter is designed;

[0106] (2) according to the transfer function of the digital filter, a unit impulse response function of the digital filter is obtained through inverse laplace transform; and

[0107] (3) a system response of the digital filter is computed by means of a system response of a thrust measurement system after variable substitution.

[0108] 2) Method for Computing System Response of Digital Filter

[0109] Output of a digital filter is required to be computed by means of a unit impulse response function thereof and through a numerical integration method.

[0110] Assuming that the unit impulse response function of the digital filter is d(t),

[00024]$\begin{array}{cc}d\left(t\right)=C_{\omega }^2\left[\delta \left(t\right)-2\left({\zeta }_1{C}_{\omega }-\zeta \right){\omega }_n{e}^{-{\zeta }_1{C}_{\omega }{\omega }_{n}t}\cos \left(\sqrt{1-{\zeta }_1^2}{C}_{\omega }{\omega }_{n}t\right)-\frac{\left(C_{\omega }^2-1\right)-2\left({\zeta }_1{C}_{\omega }-\zeta \right){\zeta }_1{C}_{\omega }}{\sqrt{1-{\zeta }_1^2}{C}_{\omega }}{\omega }_n{e}^{-{\zeta }_1{C}_{\omega }{\omega }_{n}t}\sin \left(\sqrt{1-{\zeta }_1^2}{C}_{\omega }{\omega }_{n}t\right)\right].& \left(41\right)\end{array}$$\begin{array}{cc}d\left(t\right)=d_1\left(t\right)+d_2\left(t\right),& \left(42\right)\end{array}$$\begin{array}{cc}{d}_1\left(t\right)=C_{\omega }^2\delta \left(t\right),\mathrm{and}& \left(43\right)\end{array}$$\begin{array}{cc}{d}_2\left(t\right)=C_{\omega }^2\left[-2\left({\zeta }_1{C}_{\omega }-\zeta \right){\omega }_n{e}^{-{\zeta }_1{C}_{\omega }{\omega }_{n}t}\cos \left(\sqrt{1-{\zeta }_1^2}{C}_{\omega }{\omega }_{n}t\right)-\frac{\left(C_{\omega }^2-1\right)-2\left({\zeta }_1{C}_{\omega }-\zeta \right){\zeta }_1{C}_{\omega }}{\sqrt{1-{\zeta }_1^2}{C}_{\omega }}{\omega }_n{e}^{-{\zeta }_1{C}_{\omega }{\omega }_{n}t}\sin \left(\sqrt{1-{\zeta }_1^2}{C}_{\omega }{\omega }_{n}t\right)\right]& \left(44\right)\end{array}$

are made.

[0111] The following equation

θ.sub.1(t)=∫.sub.0.sup.tθθ(τ)d(t−τ)dτ  (45) is substituted to

obtain

θ.sub.1(t)=ϑ.sub.1(t)+ϑ.sub.2(t)  (46),

ϑ.sub.1(t)=∫.sub.0.sup.tθ(τ)d.sub.1(t−τ)dτ=∫.sub.0.sup.tθ(τ)C.sub.ω.sup.2δ(t−τdτ=C.sub.ω.sup.2θ(t)  (47), and

ϑ.sub.2(t)=∫.sub.0.sup.tθ(τ)d.sub.2(t−τ)dτ  (48),

[0112] where ϑ.sub.1(t) may be directly computed according to a system response θ(t), and ϑ.sub.2(t) is required to be computed through a numerical integration method.

[0113] In order to compute an integral equation for ϑ.sub.2(t), assuming that a time sampling step is h, t.sub.i=ih(i=0, 1, 2, . . . ), and

ϑ.sub.2(t.sub.i)=∫.sub.0.sup.t.sup.iθ(τ)d.sub.2(t.sub.i−τ)dτ  (49).

[0114] Since [t.sub.0=0, t.sub.i=ih] is divided into i equal parts, for τ.sub.j=jh(j=0, 1, 2, . . . , i), a function

g(t.sub.i,τ.sub.j=θ(τ.sub.j)d.sub.2(t.sub.i−τ.sub.j)  (50) is established.

[0115] A method for computing ϑ.sub.2(t) is to use a trapezoidal integral formula to start computation, and then to use three-point and four-point simpson integral formulas for computation, which is specifically as follows:

[0116] (1) Under the situation that there are n=1 sub-intervals, the trapezoidal integral formula is used to start computation.

[00025]$\begin{array}{cc}{\vartheta }_2\left(t_1\right)=\frac{h}{2}\left[g\left(t_1,{\tau }_0\right)+g\left(t_1,{\tau }_1\right)\right]& \left(51\right)\end{array}$

[0117] (2) Under the situation that there are n=2 sub-intervals, the three-point (two intervals) simpson integral formula is used for computation.

[00026]$\begin{array}{cc}{\vartheta }_2\left(t_2\right)=\frac{h}{3}\left[g\left(t_2,{\tau }_0\right)+4g\left(t_2,{\tau }_1\right)+g\left(t_2,{\tau }_2\right)\right]& \left(52\right)\end{array}$

[0118] (3) Under the situation that there are n=3 sub-intervals, the four-point (three intervals) simpson integral formula is used for computation.

[00027]$\begin{array}{cc}{\vartheta }_2\left(t_3\right)=\frac{3h}{8}\left[g\left(t_3,{\tau }_0\right)+3g\left(t_3,{\tau }_1\right)+3g\left(t_3,{\tau }_2\right)+g\left(t_3,{\tau }_3\right)\right]& \left(53\right)\end{array}$

[0119] (4) Under the situation that there are n=2k (k=2, 3, 4, . . . ) sub-intervals, the three-point simpson integral formula is used for computation.

[0120] An integral interval [.sub.0, t.sub.2k] is divided into 2k equal parts by means of 2k+1 nodes t.sub.i=ih (i=0, 1, 2, . . . , 2k), each sub-interval has a length of h=t.sub.2k/(2k), and the three-point simpson integral formula is repeatedly used in each sub-interval such that a compound simpson formula may be obtained

[00028]$\begin{array}{cc}{\vartheta }_2\left(t_{2k}\right)=\frac{h}{3}\left[g\left(t_{2k},{\tau }_0\right)+4\underset{i=1}{\overset{k}{.Math.}}g\left(t_{2k},{\tau }_{2i-1}\right)+2\underset{i=1}{\overset{k-1}{.Math.}}g\left(t_{2k},{\tau }_{2i}\right)+g\left(t_{2k},{\tau }_{2k}\right)\right].& \left(54\right)\end{array}$

[0121] (5) Under the situation that there are n=2k+1 (k=2, 3, 4, . . . ) sub-intervals, the three-point and four-point simpson integral formulas are alternately used for computation.

[0122] An integral interval [t.sub.0, t.sub.2k+1] is divided into 2k+1 equal parts by means of 2k+2 nodes t.sub.i=ih (i=0, 1, 2, . . . , 2k+1), each sub interval has a length of h=t.sub.2k+1/(2k+1), the three-point simpson integral formula is repeatedly used in first 2k−2 (k=2, 3, 4, . . . ) sub-intervals, and the four-point simpson integral formula is used in last four points, that is, 2k−2, 2k−1, 2k and 2k+1 such that a compound simpson formula may be obtained

[00029]$\begin{array}{cc}{\vartheta }_2\left(t_{2k+1}\right)=\frac{h}{3}\left[g\left(t_{2k+1},{\tau }_0\right)+4\underset{i=1}{\overset{k-1}{.Math.}}g\left(t_{2k+1},{\tau }_{2i-1}\right)+2\underset{i=1}{\overset{k-2}{.Math.}}g\left(t_{2k+1},{\tau }_{2i}\right)\right]+\frac{17h}{24}g\left(t_{2k+1},{\tau }_{2k-2}\right)+\frac{9h}{8}g\left(t_{2k+1},{\tau }_{2k-1}\right)+\frac{9h}{8}g\left(t_{2k+1},{\tau }_{2k}\right)+\frac{3h}{8}g\left(t_{2k+1},{\tau }_{2k+1}\right).& \left(55\right)\end{array}$

[0123] During computation of a system response of a digital filter, a natural frequency ω.sub.n and a damping ratio ζ of a thrust measurement system before variable substitution are known, a natural frequency ω.sub.n1=C.sub.ωω.sub.n of an equivalent-sensitivity high-frequency thrust measurement system (realized by selecting a frequency ratio C.sub.ω>1) is selected according to the requirement of measuring thrust response time, and by selecting the damping ratio ζ.sub.1=0.7˜0.9, time of entering a steady state is further shortened.

[0124] 5. Method for Inversely Computing Thrust by Means of Equivalent-Sensitivity High-Frequency Thrust Measurement System

[0125] Since a natural frequency of an equivalent-sensitivity high-frequency thrust measurement system is greatly increased, response time of rapid-varying thrust may be measured.

[0126] On the one hand, from system response of a digital filter, that is, system response of an equivalent-sensitivity high-frequency thrust measurement system, thrust increasing time and thrust decreasing time may be indirectly determined and read. On the other hand, according to the system response of the digital filter (that is, the system response of the equivalent-sensitivity high-frequency thrust measurement system), thrust may be inversely computed by means of a transfer function and a unit impulse response function of the equivalent-sensitivity high-frequency thrust measurement system such that the thrust increasing time and the thrust decreasing time may be further confirmed. In this way, the thrust response time may be determined and read more reasonable and reliable.

[0127] 1) Establishing Integral Thrust Equation for Equivalent-Sensitivity High-Frequency Thrust Measurement System

[0128] A transfer function of an equivalent-sensitivity high-frequency thrust measurement system is

[00030]$\begin{array}{cc}{H}_1\left(s\right)=\frac{L_f}{k}.Math.\frac{{\omega }_{n1}^2}{s^2+2{\zeta }_1{\omega }_{n1}s+{\omega }_{n1}^2},& \left(56\right)\end{array}$

and

[0129] a unit impulse response function h.sub.1(t) of the equivalent-sensitivity high-frequency thrust measurement system is

[00031]$\begin{array}{cc}{h}_1\left(t\right)=L^{-1}\left[H_1\left(s\right)\right]=\frac{{\omega }_{n1}^2{L}_f}{k{\omega }_{d1}}{e}^{-{\zeta }_1{\omega }_{n1}t}\sin \left({\omega }_{d1}t\right).& \left(57\right)\end{array}$

[0130] What is known is that output of the digital filter is θ.sub.1(t), and assuming that nominal thrust is f(τ), the integral thrust equation for the equivalent-sensitivity high-frequency thrust measurement system is

θ.sub.1(t)=∫.sub.0.sup.tf(τ)h.sub.1(t−τ)dτ  (58).

[0131] The above equation is called integral thrust equation, and according to the integral thrust equation, the nominal thrust f(τ) may be inversely computed by means of θ.sub.1(t).

[0132] 2) Computing Thrust by Means of Integral Thrust Equation in Discretization Manner

[0133] A unit impulse response function h.sub.1(t) of an equivalent-sensitivity high-frequency thrust measurement system is substituted into an integral thrust equation to obtain

[00032]$\begin{array}{cc}{C}_f{\theta }_1\left(t\right)={\int }_0^{t}f\left(\tau \right)e^{-{\zeta }_1{\omega }_{n1}\left(t-\tau \right)}\sin {\omega }_{d1}\left(t-\tau \right)d\tau ,& \left(59\right)\end{array}$$\begin{array}{cc}{C}_f=\frac{k{\omega }_{d1}}{{\omega }_{n1}^2{L}_f},{\omega }_{d1}=\sqrt{1-{\zeta }_1^2}{\omega }_{n1}\mathrm{and}{\omega }_{n1}=C_{\omega }{\omega }_n.& \left(60\right)\end{array}$

[0134] What is known is that the system response of the equivalent-sensitivity high-frequency thrust measurement system is [t.sub.i, θ.sub.1(t.sub.i)](i=0, 1, 2, 3, . . . ), θ.sub.1(t.sub.0)=0 under the situation that t.sub.0=0, a sampling time step is h, t.sub.i=ih, and θ.sub.1(t)=θ.sub.1i and F(t.sub.i)=F.sub.i. For any given t.sub.i=ih (i=1, 2, 3 . . . ), an integral thrust equation is

C.sub.f=θ.sub.1i=∫.sub.0.sup.ihF(τ)e.sup.−ζ.sup.1.sup.ω.sup.n1.sup.(ih−τ)sin ω.sub.d1(ih−τ)dτ  (61),

[0135] where replacing f(τ) with F(τ) means that F(τ) is an estimated value of thrust f(τ).

[0136] An interval [0, t.sub.i] is divided into i equal parts (i=1, 2, 3, . . . ), a node is τ.sub.j=jh(j=0, 1, . . . , i), and an integrated function is

G(t.sub.i,τ.sub.j)=F(τ.sub.j)e.sup.−ζ.sup.1.sup.ω.sup.n1.sup.(t.sup.i.sup.−τ.sup.j.sup.)sin ω.sub.d1(t.sub.i−τ.sub.j)  (62),

and

G(t.sub.i,τ.sub.0)=G(t.sub.i,0)=F.sub.0e.sup.−ζ.sup.1.sup.ω.sup.n1.sup.(ih)sin ω.sub.d1(ih)  (63) and

G(t.sub.i,τ.sub.i)=F(τ.sub.i)e.sup.−ζ.sup.1.sup.ω.sup.n1.sup.(t.sup.i.sup.−τ.sup.i.sup.)sin ω.sub.d1(t.sub.i−τ.sub.i)=0  (64) are satisfied.

[0137] Three-point and four-point simpson integral formulas are required to be used for discretizing an integral thrust equation group into a linear thrust equation group through a compound simpson discretization method, which will be discussed by means of the following five situations.

[0138] (1) Under the situation that i=1, computation is started according to a trapezoidal integral formula to obtain

[00033]$\begin{array}{cc}{C}_f{\theta }_{11}={\int }_0^{h}F\left(\tau \right)e^{-{\zeta }_1{\omega }_{n1}\left(h-\tau \right)}\sin {\omega }_{d1}\left(h-\tau \right)d\tau =\frac{h}{2}{e}^{-{\zeta }_1{\omega }_{n1}h}\sin \left({\omega }_{d1}h\right)F_0& \left(65\right)\end{array}$$\mathrm{and}$$a_{10}{F}_0=C_f{\theta }_{11}\mathrm{and}{a}_{10}=\frac{h}{2}{e}^{-{\zeta }_1{\omega }_{n1}h}\sin {\omega }_{d1}h.$

[0139] (2) Under the situation that i=2, according to a three-point simpson integral formula, it can be known that,

[00034]$\begin{array}{cc}{C}_f{\theta }_{12}={\int }_0^{2h}F\left(\tau \right)e^{-{\zeta }_1{\omega }_{n1}\left(2h-\tau \right)}\sin {\omega }_{d1}\left(2h-\tau \right)d\tau =\frac{h}{3}{e}^{-{\zeta }_1{\omega }_{n1}\left(2h\right)}\sin {\omega }_{d1}\left(2h\right)F_0+\frac{4h}{3}{e}^{-{\zeta }_1{\omega }_{n1}\left(h\right)}\sin {\omega }_{d1}\left(h\right)F_1& \left(66\right)\end{array}$$\mathrm{and}$$a_{20}{F}_0+a_{21}{F}_1=C_f{\theta }_{12},$$\begin{array}{cc}{a}_{20}=\frac{h}{3}{e}^{-{\zeta }_1{\omega }_{n1}\left(2h\right)}\sin {\omega }_{d1}\left(2h\right)\mathrm{and}{a}_{21}=\frac{4h}{3}{e}^{-{\zeta }_1{\omega }_{n1}\left(h\right)}\sin {\omega }_{d1}\left(h\right)& \left(67\right)\end{array}$

[0140] may be obtained.

[0141] (3) Under the situation that i=3, according to a four-point simpson integral formula, it may be known that

[00035]$\begin{array}{cc}{C}_f{\theta }_{13}={\int }_0^{3h}F\left(\tau \right)e^{-{\zeta }_1{\omega }_{n1}\left(3h-\tau \right)}\sin {\omega }_{d1}\left(3h-\tau \right)d\tau =\frac{3h}{8}{e}^{-\zeta {\omega }_{n1}\left(3h\right)}\sin {\omega }_{d1}\left(3h\right)F_0+\frac{9h}{8}{e}^{-{\zeta }_1{\omega }_{n1}\left(2h\right)}\sin {\omega }_{d1}\left(2h\right)F_1+\frac{9h}{8}{e}^{-{\zeta }_1{\omega }_{n1}\left(h\right)}\sin {\omega }_{d1}\left(h\right)F_2& \left(68\right)\end{array}$$\mathrm{and}$$a_{30}{F}_0+a_{31}{F}_1+a_{32}{F}_2=C_f{\theta }_{13}$$\begin{array}{cc}{a}_{30}=\frac{3h}{8}{e}^{-{\zeta }_1{\omega }_{n1}\left(3h\right)}\sin {\omega }_{d1}\left(3h\right),& \left(69\right)\end{array}$$a_{31}=\frac{9h}{8}{e}^{-{\zeta }_1{\omega }_{n1}\left(2h\right)}\sin \omega \left(2h\right)\mathrm{and}$$a_{32}=\frac{9h}{8}{e}^{-{\zeta }_1{\omega }_{n1}\left(h\right)}\sin {\omega }_{d1}\left(h\right)$

may be obtained.

[0142] (4) Under the situation that i=2k(k=2, 3, . . . ) according to the three-point simpson integral formula, it may be obtained that

[00036]$\begin{array}{cc}{C}_f{\theta }_{1i}={\int }_0^{\mathrm{ih}}F\left(\tau \right)e^{-{\zeta }_1{\omega }_{n1}\left(\mathrm{ih}-\tau \right)}\sin {\omega }_{d1}\left(\mathrm{ih}-\tau \right)d\tau =\frac{h}{3}{e}^{-{\zeta }_1{\omega }_{n1}\left(\mathrm{ih}\right)}\sin {\omega }_{d1}\left(\mathrm{ih}\right)F_0+\frac{4h}{3}\underset{j=1}{\overset{k}{.Math.}}{e}^{-{\zeta }_1{\omega }_{n1}\left[i-\left(2j-1\right)\right]h}\sin {\omega }_{d1}\left[i-\left(2j-1\right)\right]{\mathrm{hF}}_{2j-1}+\frac{2h}{3}\underset{j=1}{\overset{k-1}{.Math.}}{e}^{-{\zeta }_1{\omega }_{n1}\left(i-2j\right)h}\sin {\omega }_{d1}\left(i-2j\right){\mathrm{hF}}_{2j}& \left(70\right)\end{array}$$\mathrm{and}$$a_{i0}{F}_0+\underset{l=1}{\overset{2k-1}{.Math.}}{a}_{\mathrm{il}}{F}_l=C_f{\theta }_{1i}\left(i=2k;k=2,3,.Math.\right),$$\begin{array}{cc}{a}_{i0}=\frac{h}{3}{e}^{-{\zeta }_1{\omega }_{n1}\left(\mathrm{ih}\right)}\sin {\omega }_{d1}\left(\mathrm{ih}\right),& \left(71\right)\end{array}$$\begin{array}{cc}{a}_{\mathrm{il}}=\frac{4h}{3}{e}^{-{\zeta }_1{\omega }_{n1}\left(i-l\right)h}\sin {\omega }_{d1}\left(i-l\right)h\left(l=1,3,.Math.,2k-1\right)\mathrm{and}& \left(72\right)\end{array}$$\begin{array}{cc}{a}_{\mathrm{il}}=\frac{2h}{3}{e}^{-{\zeta }_1{\omega }_{n1}\left(i-l\right)h}\sin {\omega }_{d1}\left(i-l\right)h\left(l=2,4,,2k-2\right)& \left(73\right)\end{array}$

may be obtained.

[0143] (5) Under the situation that i=2k+1 (k=2, 3, . . . ), in 2k+2 points, the three-point simpson integral formula is used for first i=0, 1, . . . , 2k−3, 2k−2 equal points, and the four-point simpson integral formula is used for last four points i=2k−2, 2k−1, 2k, 2k+1, and may be obtained,

[00037]$$\begin{gathered} \begin{array}{cc}{C}_f{\theta }_{1i}={\int }_0^{\mathrm{ih}}F\left(\tau \right)e^{-{\zeta }_1{\omega }_{n1}\left(\mathrm{ih}-\tau \right)}\sin {\omega }_{d1}\left(ih-\tau \right)d\tau = \\ \frac{h}{3}{e}^{-{\zeta }_1{\omega }_{n1}\left(\mathrm{ih}\right)}\sin {\omega }_{d1}\left(ih\right)F_0+\frac{4h}{3}\underset{j=1}{\overset{k-1}{.Math.}}{e}^{-{\zeta }_1{{\omega }_{n1\left[i-(i-\left(2j-1\right)\right]}}_h}\sin {\omega }_{d1}\left[i-\left(2j-1\right)\right]{\mathrm{hF}}_{2j-1}+ \\ \frac{2h}{3}\underset{j=1}{\overset{k-2}{.Math.}}{e}^{-{\zeta }_1{\omega }_{n1}\left(i-2j\right)h}\sin {\omega }_{d1}\left(i-2j\right)h{F}_{2j}+\frac{17h}{24}{e}^{-{\zeta }_1{\omega }_{n1}\left[i-\left(2k-2\right)\right]h}\sin {\omega }_{d1}\left[i-\left(2k-2\right)\right]h{F}_{2k-2}+\frac{9h}{8}{e}^{-{\zeta }_1{\omega }_{n1}\left[i-\left(2k-1\right)\right]h}\sin {\omega }_{d1}\left[i-\left(2k-1\right)\right]h{F}_{2k-1}+\frac{9h}{8}{e}^{-{\zeta }_1{\omega }_{n1}\left(i-2k\right)h}\sin {\omega }_{d1}\left(i-2k\right)h{F}_{2k}\mathrm{that}\mathrm{is},a_{\mathrm{i0}}{F}_0+\underset{l=1}{\overset{2k}{.Math.}}{a}_{\mathrm{il}}{F}_l=C_f{\theta }_{1i}\left(i=2k+1;k=2,3,.Math.\right),& \left(74\right)\end{array} \end{gathered}$$$\begin{array}{cc}{a}_{i0}=\frac{h}{3}{e}^{-{\zeta }_1{\omega }_{n1}\left(\mathrm{ih}\right)}\sin {\omega }_{d1}\left(\mathrm{ih}\right),& \left(75\right)\end{array}$$\begin{array}{cc}{a}_{il}=\frac{4h}{3}{e}^{-{\zeta }_1{\omega }_{n1}\left(i-l\right)h}\sin {\omega }_{d1}\left(i-l\right)h\left(l=1,3,.Math.,2k-3\right),& \left(76\right)\end{array}$$\begin{array}{cc}{a}_{il}=\frac{2h}{3}{e}^{-{\zeta }_1{\omega }_{n1}\left(i-l\right)h}\sin {\omega }_{d1}\left(i-l\right)h\left(l=2,4,.Math.,2k-4\right),& \left(77\right)\end{array}$$\begin{array}{cc}{a}_{\mathrm{il}}\frac{17h}{24}{e}^{-{\zeta }_1{\omega }_{n1}\left(i-l\right)h}\sin {\omega }_{d1}\left(i-l\right)hl=2k-2\mathrm{and}& \left(78\right)\end{array}$$\begin{array}{cc}{a}_{\mathrm{il}}=\frac{9h}{8}{e}^{-{\zeta }_1{\omega }_{n1}\left(i-l\right)h}\sin {\omega }_{d1}\left(i-l\right)hl=2k-1,2k.& \left(79\right)\end{array}$

[0144] (6) An estimated value F(τ) of nominal thrust f(τ) is computed, and assuming that the estimated value of the thrust F′(τ) to be measured is f′(τ), the estimated value of the thrust to be measured is obtained:

[00038]$\begin{array}{cc}{F}^{\prime }\left(t\right)=F\left(t\right)-f_0\left(t\right),\mathrm{and}& \left(80\right)\end{array}$$\begin{array}{cc}{f}_0\left(t\right)=-\frac{J}{L_f}(2{\mathrm{\zeta \omega }}_n{\stackrel{.}{\theta }}_0+{\omega }_n^2{\stackrel{.}{\theta }}_{0}t+{\omega }_n^2\stackrel{.}{{\theta }_0).}& \left(81\right)\end{array}$