Modeling & Simulation of Aerospace Systems - Final Course Project

Bassissi, Colombo, De Luca, MSAS – Final project MSAS – Final project Bassissi Enrico, Colombo Alessandro, De Luca Maria Alessandra 1 Introduction and real system Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) satellite was designed and operated by the European Space Agency with the aim of mapping the Earth’s gravity field with high accuracy and spatial resolution [2]. For the purpose for which it was built, the spacecraft had to fly in an orbit very close to the planet’s surface and therefore a DFACs (Drag-Free and Attitude Control System) was mandatory to counteract the atmospheric drag which would led to faster orbital decay and perturbation in the gravity gradient measurement. A graphical representation is given in Fig. 1a, whose main components can be summarized in: the electrostatic gravity gradiometer, composed by 3 pairs of accelerometers in 3 orthogonal directions to measure the perturbations; the propulsion unit, made of a Xenon’s feed system and the dual Kaufman-type ion thrusters; the flow controller for the propellant mass flow rate. (a) GOCE spacecraft, image credit: ESA [1]. (b) GOCE gravitational field model. Figure 1: GOCE gravitational field model. The aim of the following analysis is to model and simulate of the Drag Free and Attitude Control System under nominal and off-nominal conditions; for this purpose some simplifications are carried out as described in the further paragraph. 2 The physical model In the following analysis only one accelerometer is considered, aligned with the direction of relative velocity of the s/c with respect to the atmosphere (assumed rigidly rotating with the Earth). The orbital motion is modeled with drag and J2 as main perturbation, considering that SRP is of several order of magnitude less effective at GOCE operative altitude and that the gravitational pull of the moon is affecting both the s/c and the internal subsystems together so does not have any relative effect on the DFACs. For the characterization of the atmosphere, the empirical model NRLMSISE-00 is used [3]. The perturbations acting on the spacecraft give origin to a relative motion between the proof mass in the accelerometer and the capacitor plates attached to the spacecraft itself. In particular, the mass moves due to a disparity between the overall forces acting on the spacecraft and the electrostatic forces generated in the accelerometer; because of this displacement xa , the AY 2020-21 – Prof. F. Topputo; TA: C. Giordano, V. Franzese 1 Bassissi, Colombo, De Luca, MSAS – Final project capacitance of the armatures change. This difference in capacitance produce a net outflow of current going toward an ideal operational amplifier that converts the current into an output voltage Vout , which is then used to restore the proof mass initial equilibrium condition by means of a PD controller. No resistance is considered in parallel to the capacitor in the OP-Amp and no border effects are considered for the capacitors. A graphical representation of the accelerometer is given in Fig. 2a. At this point, the output voltage Vout is also used as the input for the following control valve, modeled as an electromagnetic RL circuit coupled with a mechanical mass-spring-damper system. An open loop PI controller is introduced to appropriately adjust the amount of voltage to feed the current flowing in the spool, responsible of the electromagnetic force that counteracts the effect of a spring designed to be in equilibrium condition when the valve is closed [4]; in this way the orifice’s opening is controlled. Xenon is stored in a tank and discharged reaching choked flow condition in the valve cross section which acts as a throat. After going through a ionisation process in the first chamber of the electrical thruster, the Xe+ ions get accelerated by a series of grids with high potential difference, to be then neutralized by an ulterior electron flow and discharged to produce the desired thrust needed to compensate the drag, assumed aligned in the motion direction. Vc Vbias PD C1 Fext Fe1 Kfcv Cf xa PI Vout I Ki R L Vout Fe2 Vbias C2 mfcv Vc PD Xe0 (a) Accelerometer physical model. P, T xv Cfcv . m Acceleration grid V e- Xe+ Thrust (b) Physical models of flow control valve and ion thruster. Figure 2: Physical models of the DFACs principal subsystems. 3 The mathematical model Once the physical model is sketched, the mathematical model can be derived. In particular, the states chosen to describe the system are the keplerian elements [a, e, i, Ω, ω , θ] for the orbital part,the readout voltage and the proof mass displacement and velocity [Vout , xa , va ] for the accelerometer part and at the end the integral R of Vout , the current flowing in the spool and the displacement and the velocity of the valve [ Vout , I, xv , vv ]. In order to find the accelerometer characteristic equations, the following analysis is performed starting from a rest position of the seismic mass. Computing the KVL of the outer loop and obtaining the equivalent capacitance of the subsystem, the voltage drop Vx given by the suspended mass displacement can be determined, function also of the feeding voltage Vbias and the initial electrodes-mass gap g as in Eq. (1). Vx = Vbias xa g (1) The operational amplifier circuit can be solved to find the Vout evolution in time (Eq. (4a)), whose value is fed into the PD controller (characterized by Kpa and Kda gains) to produce the AY 2020-21 – Prof. F. Topputo; TA: C. Giordano, V. Franzese 2 Bassissi, Colombo, De Luca, MSAS – Final project control variable Vc used to restore the mass nominal position, as shown in Eq. (2). Vc = Kpa Vout + Kda V̇out (2) Indeed, after the PD, the voltage on the two capacitor plates is updated and, consequently, the electrostatic forces responsible of relocating the seismic mass, as shown in Eq. (3). 1 ∆V1 = Vbias − Vc − Vx 2 1 ∆V2 = Vbias + Vc + Vx 2 ∆V12 2(g − xa )2 ∆V22 = εA 2(g + xa )2 ⇒ Fe1 = εA (3a) ⇒ Fe2 (3b) In the end, the proof mass equation of motion is found by applying together the electrostatical forces and the external ones due to the perturbing accelerations of thrust and drag (Fext = ma (T − D)/Msc , where ma is the accelerometer mass and Msc the mass of the spacecraft). The overall state equations for the accelerometer are summarized in Eq. (4).  2 2  V̇out = − 1 2 Vbias va ε A g + xa (4a)   2 2  Cf (g − xa )2  Accelerometer: ẋa = va (4b)     1   ẍa = (4c) (Fext − Fe1 + Fe2 ) ma where ε and A are the accelerometer’s permittivity and section area. Concerning the valve subsystem, the current flowing in the spool is derived from the application of the output of the open loop PI controller to an equivalent RL circuit to take into account the inductance of the coil and the overall losses (Lv = 1 mH , Rv = 0.5 Ω). The movement of the valve, instead, can be written from a force balance on the valve pin mass due to the spring kf cv , the friction damping term c and the magnetic force induced by the current in the spool which is the responsible effect to keep the valve open. The state equations of the valve are given in Eq. (5). Control valve:  ˙ dt = Vout  ∫ Vout      I Rv 1   (Kpv Vout + Kiv ∫ Vout dt) −  I˙ = Lv Lv  ẋv = vv      1    ẍv = m [kf cv (10A0 + d0 − xv ) − Ki I − cvv ] f cv (5a) (5b) (5c) (5d) The model used for the valve is a circular orifice closed by a separator that moves perpendicular to the flow direction, thus the area actually open for the passage of fluid is calculated as the circular segment specified by the displacement xv .       r0 2 −1 xv − 10A0 −1 xv − 10A0 2 cos − 1 + sin 2 cos −1 (6) Av = 2 r0 r0 Assuming that the Xenon reaches always choked flow condition in the orifice, the mass flow rate can be described as only dependent on the valve area Av and the gas properties P2 , T2 , γ. The model to compute the thrust can be simplified assuming negligible the inlet Xenon speed, 100% flow ionisation and lumping the aforementioned process and the electromagnetic acceleration in a single expression knowing the difference potential ∆V between the grids, the Xenon ions mass mi and the electrons charge e− , as shown in Eq. (7). s   γ+1 r 2(γ−1) 2 γ 2∆V e− ṁ = P 2 Av ⇒ T = ṁ (7) γ+1 RXe T2 mi AY 2020-21 – Prof. F. Topputo; TA: C. Giordano, V. Franzese 3 Bassissi, Colombo, De Luca, MSAS – Final project The dynamics of the keplerian elements are described using Gauss Planetary Equations in the {t̂, n̂, ĥ} frame as shown in Eq. (8),where the tangential direction is aligned with the velocity of the spacecraft itself. The perturbing accelerations described in this reference frame, {at , an , ah }, are obtained as sum of the single contributions of J2, drag and thrust.   2a2 v   at (8a) ȧ =   µE        r sin θ 1   (8b)  ė = 2(e + cos θ)at − a an v      1    i̇ = r cos(ω + θ)ah (8c) h GPEs: 1   (8d) Ω̇ = r sin(ω + θ)ah    h sin i         1 1 r cos θ   an − r sin(ω + θ) cos i ah (8e) ω̇ = 2 sin θ at + 2e +   a e v h sin i          h r cos θ 1   θ̇ = − 2 sin θ at + 2e + an (8f)  2 r a ev 4 4.1 Numerical integration Linear stability analysis Once the mathematical model has been developed, the equations can be numerically integrated to determine the system response. The first thing to notice is that different domains with different dynamics are involved, and therefore the problem is stiff. To validate that, the overall system of differential equations is linearized around a reference condition, chosen as the initial one. In particular, the initial values of the keplerian elements are given, while the inner system is considered in rest position: all the states are set to null value, except for xv set equal to 10A0 + d0 since the valve is still closed and no thrust is produced. Therefore, the only coupling term between the GPEs and the DFACs ones is the drag force; however, by considering that one as a given input (i.e. the force already evaluated at initial time), the two sub-systems can be decoupled and the eigenvalues are computed separately, as shown in Fig. 3a. As evident from the figure, all the eigenvalues are placed in the origin or exhibit negative real part (including the two complex conjugates ones), and therefore the system can be considered stable. However the two sub-systems, as expected, show a completely different order of magnitude: the orbital dynamics is much slower, and therefore eigenvalues are much smaller compared to ones of the control part. Notice at the end that the previous results have been obtained by setting to zero the differential equation related to the true anomaly; indeed, since this quantity, for its nature, keeps increasing, it would lead to a different result with two complex conjugates poles in the right hand side and therefore to overall instability. 4.2 Integration scheme and step-size As a result of the afore-mentioned analysis, the integration schemes ode15s and ode23tb are considered, both implementing A-stable methods and therefore suitable for stiff problems. In order to compare them for the problem under analysis, the overall system of equations is integrated for 5 revolutions of the satellite using both the integration schemes. The tolerances are set equal in both cases and are differentiated according to the kind of state (smaller tolerances are used for the orbital mechanics (10−12 ) and for the mass dynamics (10−16 ), while slightly larger tolerances are adopted for the other dynamics); in this way each state can be accurately represented without increasing too much the computational time. The integration time needed AY 2020-21 – Prof. F. Topputo; TA: C. Giordano, V. Franzese 4 Bassissi, Colombo, De Luca, MSAS – Final project (a) Linearized system eigenvalues. (b) BDFs stability region [9] and ode15s step-size histogram. Figure 3: Linear stability analysis and ode15s features. is 1.4533 s with 5371 number of steps for ode15s and 59.697 s with 114064 steps for ode23tb. The results obtained are the same, but due to its higher computational cost and less efficiency ode23tb is discarded and ode15s is the method chosen for all the following analyses. The underlying idea implemented in ode15s is the one of the backward differences formulae (BDFs) from order 1 to 5. However this implicit method is A-stable only up to order 2: indeed, by increasing the order, the stability of the BDF’s worsens and the imaginary axis is slightly overcome; therefore the step-size should be accurately chosen [6]. Nevertheless, since most of the eigenvalues of the problem belong to the real axis, they remain in the stability region no matter which is the dimension of the step. Fig. 3b shows the steps-size used by ode15s to integrate the equations, together with the region of stability of the method [9]. As evident from the histogram the dimension of the step can vary quite a bit. For example, very small stepsize are mainly involved at the beginning of the integration to represent the transients; then, an automatic adjustment procedure is adopted for choosing subsequent steps and produce an accurate solution efficiently [7]. 5 The simulation framework Fig. 4 shows the flow chart adopted for the development of the project, with a particular insight on the simulation block, in which the main steps to determine the state variables are reported. 6 6.1 System response Nominal and Off-Nominal conditions The equations are firstly integrated in nominal conditions over a t-span of 5 orbits. Fig. 5 and R Fig. 6 show the evolution of the states (except for Vout , whose trend is the same of I) and the forces. As evident, using the nominal values, the DFACs is unable to keep up with the varying conditions of the operating environment,resulting in mismatch between thrust and drag. For the off-nominal conditions, the attention is held on some of the fallible components contained in the system. In particular three possible cases are analysed: erosion phenomena for the accelerating grids of the ion thruster, a leakage in the Xenon tank and a failure of the PI controller with subsequent recovery. In all cases the failure time is computed randomly. AY 2020-21 – Prof. F. Topputo; TA: C. Giordano, V. Franzese 5 Bassissi, Colombo, De Luca, MSAS – Final project Initialisation • Spacecraft system’s data • ICS (orbital & spacecraft systems) main.m Simulation • cond.type = 'Nominal'; simulation.m Linearisation • Reference conditions • Eigenvalue problem • Stability region Drag • NRLMSISE-00 atmospheric model • Relative velocity wrt air [a, e, i, Ω, ω, θ] • Orbital position to ECI frame States • Orbital • Accelerometer • Valve J2 • Position in ECEF frame • Earth coefficients [𝒙𝒂 , 𝒗𝒂 , 𝑽𝒐𝒖𝒕 ] • Readout voltage • PD controller • Electrostatic forces acting on the proof mass + [‫ 𝒕𝒖𝒐𝑽 ׬‬, I, 𝒙𝒗 , 𝒗𝒗 ] • PI controller • RL circuit • Electromechanical control valve • Variable throath area Sensitivity Analysis • Select parameters • ±60% of nominal value • std deviation analysis Simulation • cond.type = 'OffNominal’; Drag Acceleration Effect + Gauss Planetary Equations in ෠ frame [𝑡,Ƹ 𝑛, ො ℎ] Overall acceleration on the spacecraft J2 Acceleration Effect cond.type = 'Off-Nominal’ • cond.failure = {‘Grid Erosion’, ‘Xe Leakage’, ‘PI Failure & Recovery’} Xenon Tank • Mach = 1 discharge ሶ • Mass flow rate [𝒎] Optimisation • Select parameters • Lower/upper boundaries • Objective function • patternsearch.m + Ion Thruster • Acceleration grids • Thrust generation External forces acting on accelerometer proof mass Thrust Acceleration Effect + Figure 4: Simulation framework flow chart. Figure 5: Keplerian elements evolution in Nominal condition. A failure of the pressurization system of the Xenon tank produces a variation of the operating pressure. Assuming the system to be able to have a clear readure of the latter, it will automatically adjust the orifice opening trying to match the wanted mass flow rate by increasing Av and consequently xv , as shown in Fig. 7. Then a degradation of thruster performances is considered, here modeled as an sped up erosion phenomenon of the accelerating grids [8]. The achieved results show a similar behaviour to the previous off-nominal case. Lastly, assuming AY 2020-21 – Prof. F. Topputo; TA: C. Giordano, V. Franzese 6 + Bassissi, Colombo, De Luca, MSAS – Final project Figure 6: Control system states and Forces evolution in Nominal condition. a null output VI to be the result of malfunction in the controller, causes the fast run out of accumulated current in the spool and consequently closure of the valve and a drop of thrust to zero; assuming however the possibility of a recovery strategy, the system is able to recover its nominal behavior after a small transient. Figure 7: Forces evolution and control valve displacement in the three Off-Nominal conditions. AY 2020-21 – Prof. F. Topputo; TA: C. Giordano, V. Franzese 7 Bassissi, Colombo, De Luca, MSAS – Final project 6.2 Sensitivity analysis The sensitivity analysis is carried out to examine how uncertainties in some independent parameters contribute in the overall behaviour of the DFACs. The assessment of that is obtained through the evaluation of an objective function related to the compensation of the drag decay by means of the thrust; its formulation is inspired by the output error method [5] and defined as in Eq. (9), where the t-span considered is the one of one orbital period. T J= X 1 (T − D)2 2 var(T − D) (9) t=0 The variables of interest considered are 11: the gains Kpa , Kda , Ki , Kpv , Kiv , the valve spool mass and spring coefficient mf cv , kf cv , the variables of the Xenon P2 , T2 and the resistance and inductance introduced for the electro-mechanical modeling of the valve Rv , Lv ; on these quantities a perturbation of ±60% with respect to their nominal values is applied. The results are shown in Fig. 8, where the blue dots represent the evaluations of J for every variation of each parameter, the red dots represent the mean values and the yellow boxes the amplitudes of standard deviation with respect to the mean value for each parameter. The conclusions that can be drawn is that, out of the six optimisable parameters, actually only three are effectively modifying the system response. Concerning instead Rv and Lv , the parameters show opposite behaviour: the inductance has a weak impact on the objective function, whereas the resistance is actually much more influential. Furthermore among the other parameters under analysis, the system turn out to be very sensitive to changes in the valve spring stiffness and the operative pressure of the Xe. Figure 8: Standard deviation of main parameter on the objective function. 7 System optimization In order to achieve the best possible performance, with an optimised thrust with same amplitude and phase of the drag, patternsearch.m global optimiser is used. Starting from the previous sensitivity analysis,the optimisation is carried out on the higher std parameters, to which the system response is more sensitive, by setting proper bounds. However, considering kf cv and the pressure PXe as unmodifiable data, the optimization is performed on Kpa , KI , Kiv and Rv . The objective function and the t-span used are the same of the sensitivity analysis; the convergence AY 2020-21 – Prof. F. Topputo; TA: C. Giordano, V. Franzese 8 Bassissi, Colombo, De Luca, MSAS – Final project Table 1: Optimised Parameters. Nominal Optimised Kpa KI Kiv Rv [Ω] J 1e6 0.2 3 0.5 1223.47 2.1561e5 0.32507 4.508 0.25025 787.187 of the algorithm is reached for the values shown in Table 1, greatly reducing the value of J. Simulating now the response with those values, the result is shown in Fig. 9; excluding the initial transient, the maximum error between thrust and drag in absolute value is 0.0752 mN, compared to the 0.9201 mN of the nominal case, confirming the significant improvement obtained. Figure 9: System response with optimised parameters. AY 2020-21 – Prof. F. Topputo; TA: C. Giordano, V. Franzese 9 Bassissi, Colombo, De Luca, MSAS – Final project References [1] ESA GOCE (Gravity field and steady-state Ocean Circulation Explorer). 2002 EO Portal Directory, https://earth.esa.int/web/eoportal/satellite-missions/g/goce# LNAqF13c8Herb [2] Enrico Canuto, Luca Massotti. All-propulsion design of the drag-free and attitude control of the European satellite GOCE. 2009 Acta Astronautica, volume 64, 325–344 [3] J. M. Picone, A. E. Hedin,D. P. Drob, A. C. Aikin. NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues. 2002 Journal of Geophysical Research: Space Physics, volume 107, SIA 15-1-SIA 15-16 doi: 10.1029/2002JA009430 [4] N. D. Vaughan, J. B. Gamble. The Modeling and Simulation of a Proportional Solenoid Valve. 1996 ASME American Society of Mechanical Engineers, volume 118, 120–125. [5] Marco Lovera. Output-error model identification: linear time-invariant systems. 2020 Dipartimento di Scienze e Tecnologie Aerospaziali, Politecnico di Milano. [6] Lawrence Shampine, Mark Reichelt. The matlab ODE suite. 1997 SIAM Journal on Scientific Computing, volume 18, doi: 10.1137/S1064827594276424 [7] H.A. Watts. Starting step size for an ODE solver. 1982 Journal of Computational and Applied Mathematics, volume 9, 177-191 [8] Miguel Sangregorio, Kan Xie, Ningfei Wang, Ning Guo, Zun Zhang. Ion engine girds: Function, main parameters, issues, configurations, geometries, materials and fabrication methods. 2018 Chinese Journal of Aeronautics, volume 31, 1637-1639 [9] Nick Trefethen. Stability Regions of ODE Formulas. 2011 Mathworks, https: //www.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/23972/ versions/22/previews/chebfun/examples/ode/html/Regions.html AY 2020-21 – Prof. F. Topputo; TA: C. Giordano, V. Franzese 10