Dynamic-force extraction for micro-propulsion testing : Theory and experimental validation

A dynamic-force extraction, based on the least-squares method, is proposed for micro-propulsion testing. Having modeled the displacement oscillation of a micro-newton torsional pendulum, the time evolution of the dynamic force may be calculated if the stand constants are well calibrated. According to the linear characteristic of the motion equation, a reconstruction of the dynamic thrust reduces to solving linear equations. The simulation analysis shows that the error is affected by the sensor noise and the low-pass filter as well as the sampling rate. Validation experiments were performed showing that this method reconstructs the dynamic force well up to 8 Hz with an error less than 15 μN. The noise-induced error moreover varies little with frequency.


I. INTRODUCTION
Accurately controlling satellites at highly precise attitudes and/or orbits is vital in space missions requiring gravitationalwave or gravity-field measurements.To satisfy these requirements, micro-newton thrusters, such as cold gas engines and electromagnetic thrusters, were studied.Directly measuring such small thrusts is nearly impossible because the thrust is so miniscule that the device is unable to response to the changes it causes.Different kinds of thrust stands have been developed to accurately measure their performance.These stands are based on the torsion balance, [1][2][3][4] vertical pendulum, 5 indirect counterbalanced pendulum, 6 double pendulum balance, 7 and even the magnetically levitated balance. 8Thrust is deduced using balance displacement.The tools often employed for these measurements include laser interferometers, 9 fiber optic linear displacement system (LDS) sensors, 3 heterodyne interferometers, 7 capacitive sensors, 5 as well as linear variable differential transformer (LVDT) sensors, 1,2 all of which are susceptible to electromagnetic interference and radio frequency (RF) interactions when testing both laboratory DC and RF-powered plasma thrusters. 11ost of these balance displacement systems have great accuracy and resolution.For example, Jamison 1 measured thrusts as small as 86 nano-newtons and Soni 3 developed a stand with a resolution of 10 nN for steady force measurements.Jarrige 10 measured a cold gas thruster to a resolution of 20 nN and a measurement bandwidth (MBW) of 0.1 Hz.Increasing MBW decreases the resolution.
Nevertheless, the challenging performances demanded of high-precision thrusters are thrust level, accuracy, thrust noise, and dynamic response.1,12 The required a) Author to whom correspondence should be addressed: lifei@imech.ac.cn thrust range is between 0 and 3 mN with a resolution below 0.1 and 1 µN and a noise level of 0.1-1 µN/Hz 1/2 , in a frequency bandwidth range of 2-10 Hz. 12 To enlarge the MBW of the thrust stand, many researchers sought new methods to extract dynamic force.These methods may be divided into two categories: increasing eigenfrequency and dynamic modeling.
To increase the eigenfrequency of the stand, auxiliary equipment or a new design in controlling the stand system is required.Using two parallel metallic plates and a FP cavity, Canuto 12 developed a nano-balance with a 13.5-Hz resonance frequency and sub-micro-newton resolution.A large natural frequency made it suitable for dynamic measurements up to 2 Hz.Orieux 13 used a feedback control loop to increase the frequency of a pendulum balance up to 42 Hz despite its low sensitivity of 25 µN.This method is widely used for measuring thrust noise.Active control suppresses its resonance and reduces the noise by more than one order of magnitude, 7 which is significant for long-term stability.Hagiwara 14 developed a proportional-integral-derivative controller to suppress the displacement of the arm of the torsion balance during thrusts.Masuda 15 studied the sensitivity of the torsion pendulum.With strong magnetic damping, a feedback control system was used to achieve long torsional periods and high sensitivity for the balance at low frequency (e.g., 1 mHz).
The second category is dynamic modeling.By analyzing the equation of motion of the thrust stand, the time-dependent force can be calculated using mathematical methods and the displacement history of the balance.A typical method is based on derivation.During this process, the first differential (speed) and second differential (acceleration) of the displacement history, as well as the well-calibrated stand constants (e.g., inertia and spring coefficient) are used. 16,17Displacement-sensor samples of high frequency and high resolution are required.Therefore, this method has an obvious drawback in that the derivative may incur large errors arising from coupled noise in the displacement data.D'Souza 19 used wavelet and Fourier analyses to reduce noise and directly transformed the control differential equation into a difference equation.Once the displacement was measured, the corresponding force is obtained using this control equation.
In the current study, a method is proposed to reconstruct the dynamic force from the displacement data.It uses the property of superposition of the control equation solutions, a basic functional idea, and a least-squares (LS) method.Verification experiments were performed to study the applicability of this method as well as the relationship between error and noise.

II. MICRO-THRUSTER BALANCE SYSTEM
As in former studies, [1][2][3][4] the thrust stand used here is based on a torsion pendulum (Fig. 1).The measurement range of this stand is 1-3000 µN for steady-force measurements with an accuracy of ∼1 µN under ambient noise. 18ts main components included a pendulum arm, calibration coil, damper, displacement sensor, and two flexural pivots (Riverhawk, 6004-600).Each pivot has a torsional spring rate of 1.92 × 10 −4 Nm/deg.The displacement sensor is the most important component of this balance system.A capacitive sensor (Fogale, Mc900) having a range of 100 µm and an accuracy of 2 nm was used.This pendulum holds an axial load of up to 2 kg.The 50-cm length pendulum arm was used to balance gravity and converts the thrust into a displacement.Similar to reports in the literature, 20,21 the electromagnetic calibration technique was used here.A coil was fixed on the base plate opposite a magnet fixed to the arm.This coil is connected with a high-precision current source making a magnetcoil system supplying a known force from 2 to 10 000 µN.A damper, comprising a circular permanent magnet and a copper plate, operates based on Lenz's law but also acts as a counterweight.Measurement and control systems and homemade software were developed to record and analyze the sensor signal, the on/off switching of the thruster, the change in distance between arm and sensor, and the magnet-coil gap.The gap ensures a reliable measurement of micro-thrusts in a vacuum chamber.

A. Equation of motion of the pendulum
When the torsion pendulum undergoes small-angle swings, the motion of the pendulum arm follows the simple harmonic equation of forced vibration.The control equation is 22,23 where θ represents the angular position relative to the initial position, L is the distance between pivot center and thrust action point, I (kg m 2 ) is the moment of inertia of the pendulum, c (N s m) is the damping constant, and k (N m/rad) is the spring constant.Setting and substituting into Eq.( 1) yields where ζ is the coefficient of damping, ω n is the natural frequency, and ω d is the frequency of damped motion.For small deflections in our cases, the small-angle approximation may be used to deduce the angular deflection θ ≈ x/L s , where x represents the displacement of the pendulum, and L s represents the distance between the pivot center and the displacement sensor.
From theory, the amplitude of the pendulum obeys the amplitude-frequency function given by the magnification factor, 2,24 with relative frequency Setting the coefficient of damping to 0.1, Fig. 2 plots the log-log relationship between the rate of amplitude amplification and relative frequency.In Ref. 25, Saulson found a frequency dependence to damping.Nevertheless, the coefficient of damping is treated as constant, independent of frequency as assumed in Refs.16 and 19.For our situation, it is an appropriate assumption that is validated by experiments.
As can be seen, the amplitude ratio changes slightly at low frequencies and increases rapidly to a maximum value at the resonance frequency.At higher frequencies, the amplitude FIG. 2. Amplitude of the pendulum response versus relative frequency at ζ = 0.1.decreases rapidly; for example, it is 1% when the frequency is 10 times the natural frequency (Fig. 2).

B. Basic functional idea and the LS method
Based on the functional idea and the LS method, a new method for dynamic-thrust measurements is proposed.The unknown thrust can be reconstructed by solving Eq. ( 5).The basic functional idea used is to assume a force function expressed as {f τ }; here the curly braces signify a sequence of numbers with f τ denoting any one of these numbers.That is, f τ defines the function that ensures Eq. ( 5) is expressible in theory and is here denoted by x τ .The expression for x τ contains f τ , which means that x τ is represented by f τ .The theoretical solution of the displacement expressed by f τ and the measured experimental values are used in the variance formula to find the function f τ that makes the variance a minimum.That f τ is the force that generates the displacement.Thus, this method seeks a functional solution for the force.
For a self-contained description of the method, the wellknown LS method is described briefly here.The variance R of a function φ j (e.g., the response of the balance to excitation), from the known data y i (e.g., measured data of the displacement sensor), is defined as where a j are unknown coefficients (e.g., the dynamic force we want to measure), m represents the number of samples (e.g., number of displacement data), and n represents the number of a j .Then seeking the extremum of the multivariate function, we have . . .
This set of linear equations for a j is in a typical form to which the LS method applies.Measuring the dynamic force is achieved by adding the pendulum control equation, Eq. ( 5), and then applying the LS method, which is described next.unit step input is known [Eq.( 12)], it can be put into the LS method, Eq. ( 10), then used to replace the expression ϕ j , and obtain Eq. ( 15) by linear rearrangement.By solving this linear system of equations, the dynamic force may be finally reconstructed.

C. Dynamic-force extraction method
In more detail, the method by which to extract the dynamic-force is based on discretization of the equations and the superposition of solutions of the governing equation.In Fig. 4, the continuous force can be assumed to be a sequence of superimposed forms of quasi-impulse actions f (t) = f τ δ(t − τ), where f τ is the size of the impulse at t = τ and δ(t − τ) is the unit impulse function.
Let the function x(τ) represent the theoretical solution at time τ, the time the impulse acts.An arbitrary displacement function takes the form In practice, instead of the impulse function, a narrow rectangular impulse is used, having the benefit of overcoming interval limitations.In a damped balance system, the effect of the quasiimpulse is equivalent to a pair of time-inverted step functions.At time zero, the unit step function effect x(t) is 24 where ζ, ω n , and ω d are the coefficient of damping, the natural frequency, and the frequency of damped vibration, respectively.The quasi impulse becomes Substituting Eq. ( 13) into Eq.( 11) gives where x(t) is the displacement, f τ is the force at time τ, m is the number of points between the initial and final times.Substituting Eq. ( 14) into Eq.( 10) gives where x i is the displacement function caused by the unit impulse at time t i with τ its time of duration.Here f i denotes the size of the impulse, i.e., the dynamic force; our target, which can be more positive or negative.

IV. EXPERIMENTAL A. Experimental setup
To validate the dynamic-force extraction method, experiments were performed using the torsion balance (Fig. 1).During these experiments, the key issue is how to exert a known force on the balance for reference in comparison with the reconstructed forces.An electromagnetic calibration device [Fig.5; labeled (8) in Fig. 1] was used to generate the external excitation.It has two components, a permanent magnet and a coil.The flat cylindrical magnet, made of rubidium nickel alloy, has a high magnetic energy density.Its superficial magnetic field intensity is about 4000 G.The circular hole in the core of the magnet makes it easier to fix on the side of the pendulum arm.A special copper coil of radius 12.5 mm is placed facing parallel to the magnet (Fig. 5); its position and angle of pitch are precisely controlled using a 5-axis kinematic mount.Its center is coincident with the axis of symmetry of the magnet.A distance of 230 µm is maintained between magnet and coil.
With a resistance of about 0.6 Ω, this coil is connected in series with a 50-Ω resistance.A current-source-type signal generator is used to generate the different kinds of waveforms, corresponding to different dynamic forces.To avoid electric heating, the current is kept below 200 mA at an accuracy of better than 0.1 mA.After careful calibration, an oscilloscope is used to monitor the voltage/current of the coil for the known force.
The current-force coefficient can be measured using a high precision weighing balance.By placing the permanent magnet horizontally, the electromagnetic force between coil and permanent magnet can be measured by weighing the balance when the current changes.The rate of change between sensor output and force was deduced to be Aided by this electromagnetic force generator, validation experiments were conducted with an arbitrary dynamic force.Three excitation waveforms were used to generate known dynamic forces, including square, sinusoidal, and sawtooth waveforms.The frequency range of these waveforms was from 0.05 to 10 Hz in intervals of 0.5 Hz and the amplitudes ranged from 120 to 300 µN.The sampling rate of the displacement sensor was set at 100 Hz.

B. Parameters calibration
One needs to know all the parameters including ζ, ω d , and ω n before Eq. ( 15) can be solved.To calibrate these parameters, modeling and measuring the vibration of the torsional pendulum is a feasible scheme.If f (t) = 0, Eq. ( 5) becomes a homogeneous differential equation.The measured displacement of the free vibration may be compared to its theoretical solution, where A and ϕ are determined by the initial conditions.The solution can be factored into an attenuation x = e −ζt and an oscillation x = cos(ω d t + ϕ).The frequency ω d is obtained easily from the vibration period.The decay rate ζ can be deduced by taking the logarithm of the attenuation, ln(x) = −ζt, Hence, ζ is found using two or more points, with maximum displacements normally in different periods.In our experimental conditions, the values of these parameters were ζ = 0.69, ω n = 10.59 rad/s, and ω d = 7.66 rad/s.

V. RESULTANT ANALYSIS A. Typical reconstruction results
To verify the dynamic-force extraction method, three different waveforms for the coil current (corresponding to force) were applied to the torsional pendulum (Fig. 6).The initial phase of the input force is non-zero for better universality.The frequency of these input force waveforms was 0.1 Hz; the amplitude was 120 µN.The original displacement data and the reconstruction results (Fig. 7) show the responses of the damped pendulum to the various forces.In particular, in Fig. 7(b), damped oscillations had continued for more than 5 s.Hence, only the steady force lasting for about 10 s can be extracted using the conventional method.The time-dependence of the reconstructed force was found to coincide basically with the sawtooth input force [Fig.7(c)], indicating that this method is feasible for dynamic-force extraction.

System noise analysis
Apparently, noise in the displacement data has impacted on the reconstruction accuracy of the dynamic force.The inherent system noise, including balance and sensor, must be studied.A static experiment with no force loading was recorded for system-noise analysis.A fast Fourier transform (FFT) (Fig. 8) shows that the amplitude of noise varies with frequency.The downward sloping line is the linear fitting of the noise at low frequency (below 1 Hz) extrapolating the expected 1/f noise.At low frequency, the noise spectral density is mainly composed of signals from seismic motion. 16At frequencies larger than the natural frequency (about 1.1 Hz), the noise is suppressed by the pendulum response (Fig. 2).In Fig. 8  there is a clear peak around 10.5 Hz corresponding to the filter of the digital-to-analog converter for the displacement sensor.

Effect of noise and frequency
To study the effect of noise on the reconstructed force, a simulation analysis was performed in which noise was added to the theoretical displacement.The specific steps are as follows: First, a theoretical motion is simulated using Eq. ( 14) given the balance was excited by a sinusoidal force at a specific frequency.Then random white noise is added at the desired energy level to this theoretical motion and used as input data for the reconstruction program.Finally, the reconstructed force FIG. 8. FFT analysis of the system noise.Note that a low-pass filter (f c = 9.71 Hz) was used in the reconstruction process to reduce the displacement noise at high frequencies.More details and analysis concerning the filter is given in the Appendix.
is compared with the initial sinusoidal force to analyze the accuracy.
A frequency range from 0.1 to 10 Hz for the sinusoidal driving force was used and the amplitude for the random noise was varied from 0 to 30 mV (0%-3%).Similar to the experimental data, artificial noisy data used to reconstruct dynamic force and measurement error may be calculated using the known sinusoidal forces.Figure 9 plots the relationship between reconstruct error (standard deviation) and noise level.The first feature to note is that error increases as noise levels rise for frequencies less than 8 Hz.At zero noise, the reconstructed force error is about 2 µN, whereas it is about 7.5 µN for 10 mV noise and 14.5 µN for 20 mV noise.The error deviation also increases as noise strengthens.It tends to zero in the absence of noise and less than 15 µN for a 20-mV noise level.The second feature of the error-noise relationship is that error increases notably when the frequency exceeds 8 Hz for all noise levels.This is caused by the low pass filter mentioned above.Here, the cutoff-frequency of the filter used is 9.71 Hz.It attenuates the input displacement signal with frequencies above 8 Hz (Fig. 12).Within the range of 8 Hz, the error remains nearly stable while the deviation in error increases slightly with frequency.The red stars in Fig. 9 mark the experimental results as for the results from the sinusoidal waveform [Fig.7(a)].Its trend with frequency coincides with the simulation results.For frequencies below 8 Hz, the error of the measured dynamic force is about 10-14 µN and insensitive to frequency.

Effect of cut-off frequency f c
Figure 9 shows that for a given noise level, the extractedforce error almost remains constant from 0.1 to 8 Hz at fixed f c of 9.71 Hz.For different f c , this error changes.As shown in Fig. 10, for noise of 17.5-mV (approximate experimental noise), the force errors are about 8.5 µN and 19 µN for f c of 8.0 and 11.2 Hz, respectively.The reconstructed force error increases rapidly as f c increases.

Effect of sampling rate
The sampling rate of the displacement sensor may also influence the reconstruction accuracy.Subject to the same level of noise (10 mV, 1%), four sampling rates (100, 200, 500, and 1000 s −1 ) were used in a simulation analysis to obtain the response (Fig. 11).The reconstruction results became better as the sampling rate increased.This is because at a high sampling rate, more displacement elements are used for the summation in Eq. ( 9), thereby reducing the equivalent noise.Nevertheless, the computation time of the reconstruction increased FIG.11.Reconstruction result at different sampling rates.
significantly with the higher sampling rates.Therefore, it is best to reduce the sensor noise and choose an appropriate sampling rate.

VI. CONCLUSION
A mathematical method was proposed to extract the timedependent force for micro-propulsion.A functional idea and LS method were used to reduce the measurements needed to solve the system of linear equations.Its main characteristic is that no auxiliary equipment is required.The MBW of the conventional torsion balance is improved significantly.
The dynamic-force extraction method and calibration of the stand parameters (i.e., ζ, ω d , and ω n ) were studied.Validation experiments were performed using an electromagnetic device.The dynamic force was found to be well extracted from the oscillation of the balance with the reconstructed force coinciding with the known input force.Further simulation analysis shows that the reconstructed error increases as sensor noise or f c increases.Moreover, the reconstruction results become better as sampling rates increase but at the expense of computation times.For the low-pass filter ( f c = 9.71 Hz) used, the MBW is up to 8 Hz and the error for the measured dynamic force is about 10-14 µN at a sampling rate of 100 Hz. decrease considerably the error of the reconstructed force.A lower f c corresponds to a lower force error as also confirmed by Fig. 10.Nevertheless, a lower f c loses dynamic information of the force.Therefore, f c must be optimized to reach a compromise between dynamic-force accuracy and bandwidth.The filter parameter values chosen were PF = 40 rad, SF = 100 rad, MAP = 0.05 dB, MAS = 30 dB, and FS = 100 s −1 corresponding to f c = 9.71 Hz.Hence, noise for frequencies above 8 Hz was removed (Fig. 12).

Figure 3
Figure 3 summarizes the basic process involved in reconstructing the dynamic force.Being a linear equation, Eq. (5) has the superposition property and has a theoretical solution for each special case.If the pendulum response x(t) to the

FIG. 4 .
FIG. 4. Sketch illustrating the discretization of a continuous function.

FIG. 12 .
FIG. 12. Contribution of 17.5-mV white noise to the reconstructed force with/without filter.