COUPLED SIMULATION OF THE PROPULSION SYSTEM AND VEHICLE USING THE ESPSS SATELLITE LIBRARY

The paper documents the implementation and validation of the coupled simulation of the propulsion system and vehicle performed during the 4th development phase of the ESPSS (European Space Propulsion System Simulation) library running on the existing platform EcosimPro . This covers a signi¦cant update of the spacecraft propulsion system modeling: the Fluid §ow, Tanks and Combustion chamber components are updated to allow coupling to the vehicle£s motion, the Archimedes pressure coming from acceleration and rotations given by the vehicle or by any perturbation forces are taken into account, several new features are added to the Satellite library along with new components enabling full attitude control of a platform. A new powerful compact equation is presented for solving elegantly the Archimedes pressure coming from combined acceleration and rotation in the most general case (noncollinear). Eventually, a propulsion system is modeled to check the correct implementation of the new components especially those dealing with the e ̈ects of the mission on the propulsion subsystem.

The paper documents the implementation and validation of the coupled simulation of the propulsion system and vehicle performed during the 4th development phase of the ESPSS (European Space Propulsion System Simulation) library running on the existing platform EcosimPro .This covers a signi¦cant update of the spacecraft propulsion system modeling: the Fluid §ow, Tanks and Combustion chamber components are updated to allow coupling to the vehicle£s motion, the Archimedes pressure coming from acceleration and rotations given by the vehicle or by any perturbation forces are taken into account, several new features are added to the Satellite library along with new components enabling full attitude control of a platform.A new powerful compact equation is presented for solving elegantly the Archimedes pressure coming from combined acceleration and rotation in the most general case (noncollinear).Eventually, a propulsion system is modeled to check the correct implementation of the new components especially those dealing with the e¨ects of the mission on the propulsion subsystem.

INTRODUCTION
In the frame of the 4th development phase of the ESPSS library [13] running on the EcosimPro platform [4], the implementation and validation of the ¤coupled simulation of the propulsion system and vehicle¥ is discussed.
EcosimPro is a object-oriented Physical Simulation Modeling tool dedicated for system analyses.Its visual simulation interface is capable of solving various kinds of dynamic systems represented by writing equations and discrete events.It can be used to study both transients and steady states.For example, with the propulsion libraries ESPSS from European Space Agency (ESA), the user draws (and designs at the same time) the propulsion system with components of that speci¦c library consisting of tanks, lines, ori¦ces, thrusters, tees, etc.The user can enhance the design with components from the thermal library (heaters, thermal conductance, and radiators), the control library (analogue/digital devices), the electrical library, etc.
The present work follows the ¦rst step carried out in the previous phase regarding the ¤evolutionary behavior of components,¥ where a speci¦c Satellite library was set up for the simulation of the Attitude and Orbit Control System (AOCS) e¨ects on the propulsion subsystem.It included the §ight dynamic (orbit and attitude) capabilities for a full spacecraft with orbital and attitude perturbations.The present work focuses on a signi¦cant update of this Satellite platform propulsion system modeling and the coupling is performed by updating the Fluid §ow, Tanks and Combustion chamber components while taking into account the Archimedes pressure coming from acceleration and angular rates given by the vehicle or by any perturbation forces (and torques).
The paper presents ¦rst several new features added to the Satellite library as well as some new components for performing full attitude control of a platform.It is followed by a new compact equation for solving elegantly the Archimedes pressure and ¦nishes with some realistic validation test cases of bipropellant satellites in Low Earth Orbit (LEO).

ESPSS BACKGROUND
EcosimPro is a Physical Simulation Modeling tool developed for ESA by Empresarios Agrupados Internacional (Spain) since 1989.EcosimPro was a precursor and now, with its 28 years of careful growing, it belongs to the last generation of the common engineering tools after computer aided design (CAD) and integrated engineering analysis tools available on classical laptops.The kernel of EcosimPro is an expert solver of all the equations set in the di¨erent components of a system.Thanks to such expert solver, the tool allows to manipulate components like objects that can be independently further developed with more sophisticated equations.EcosimPro is based on a visual simulation tool for solving simple and complex physical processes that can be expressed in terms of equations (including ordinary di¨erential equations and di¨erential-algebraic equations) and discrete events.
Practically, the modeling of physical components is based on a basic ¤Ecosim-Pro language¥ (EL), an object-oriented programming language which is very similar to other conventional programming languages (Basic, C++) but is very powerful to write any equations and di¨erential equations for modeling continuous and discrete processes.EcosimPro employs a set of libraries containing various types of components (mechanical, electrical, pneumatic, hydraulic, etc.) which can be interconnected to model complex multidomain dynamic systems.The ESA ESPSS is a set of EcosimPro libraries written to model all aspects of a functional propulsion system.As a tool, ESPSS is relying on onedimensional (1D) §ow equations, thermodynamic relationships, and real §uid properties, there is no need for fudge factors; therefore, the results of the simulation could be considered as general as long as the §ow is 1D and homogeneous (either as monophase, or two-phase state or as a mixture).

SATELLITE LIBRARY
The satellite library contains the main components needed to build a satellite AOCS with the major interactions and perturbations that occur in §ight [5,6]: Moon, Sun, and other planets£ perturbations, forcing the implementation and use of thrusters for the North-South station keeping of a GEO (Geostationary Earth Orbit) satellite (the addition of planets allows analyses relative to interplanetary travels, §yby, etc.); Earth §atness or so-called ¤J2¥ perturbations to be used for the heliosynchronous satellites; Sun pressure interaction with the solar arrays (SAs) of the satellite; drag forces from upper Earth atmosphere; gravity gradient torques; and magnetic torques from the interaction of the Satellite magnetic moment with the Earth magnetic ¦eld.
In addition, a control of the satellite attitude can be added (using the existing components of the existing CONTROL library) for performing an Earth Pointing with actions on a set of reaction wheels (RWs) or on a set of thrusters.Finally, the components of the library are shown in Fig. 1.
Three ports have been designed for the transfer of information between the components of the library: (1) Force port: this is a multifunctional port allowing inputs from a set of thrusters, RWs, SAs, Drag Areas (DAs), Magnetotorquers (MTs), and gravity gradient (GG) with the port directions set at IN for the satellite frame and OUT for all other components.The variables of the port are of type SUM in order to automatically account for all mass §ow rates, forces, moment, angular momentum, and power coming from all connected components; (2) State port: this is a multipurpose port allowing the attitude and orbit control, three-dimensional (3D) visualization, as well as the needed inputs for the SAs, DAs, GG, and Tanks.For the port directions, there are set at OUT for the satellite frame and IN for all other components; and (3) InteractionsFluids port: this port allows to communicate interactively the mass, the location of the center of mass (CoM), and the inertia matrix of each liquid §uid in the tanks component.

Frame Component
The major component of the Satellite library has been improved in the 4th phase by adding the perturbation from the planets (from Mercury to Pluto) in order to get the possibility to enable interplanetary §ights.Its description is recalled here for completeness: the frame component is in charge of solving the §ight dynamic with input data of the initial orbit and with nongravity (i.e., nonbody) forces vectors coming from the thrusters (including reaction control thrusters (RCTs), combustion chamber nozzle, or nozzle extension from other ESPSS library), SA, DA.The perturbation forces coming from Moon, Sun, and all other planets gravity are directly solved inside the component (actually, the complete set of equations is used for computing the perturbation in order to solve accurately any §y-by or orbit around any planet).
The ¦rst main equation of the component is based on the momentum equation as follows: where � v is the velocity vector of the satellite frame with respect to the Earth centered inertial (ECI) frame; G is the gravitational constant; M is the mass of the focus body (Earth); � r is the radius vector from the focus body center; � T is the resultant thrust vector applied on the satellite CoM, in ECI; � P is the perturbation and interaction force vector in ECI frame; and m the instantaneous mass of the satellite that is corrected dynamically by the thrusters consumption if any.
The second main equation of the component is based on the dynamic equation for the attitude of satellite around its CoM as follows: The source terms come from the torques and perturbations due to the thrusters, MT, GG, and RWs, and from the perturbation torques due to the other mobile parts of the satellite (SA) and with � Ÿ being the instantaneous rotation with respect to the inertial orientation frame of the satellite and the angular momentum with respect to the satellite of the RWs � H RWi and of the other mobile parts � H i , but with all their coordinates and derivative written in the satellite axis.The angular momentum � H takes into account the possible ¤evolutionary behavior¥ of the §uid in the tanks through the inertia matrix update.
The attitude angles of the satellite must be known in order to properly set the local data like the thrust with respect to the satellite frame into the ECI frame.Those can be given by the Cardan angles (yaw, pitch, roll) of the satellite axis with respect to the orbital frame.The attitude angles of the orbital frame can be given by the Euler angles (precession, nutation, proper rotation) with respect to the ECI frame.But in order to ensure robustness in the solutions, the equations dealing with the attitude angles are based upon the quaternion theory.Thus, the general equation of the quaternion theory is to be solved simultaneously with the other previous equations.
The integration of Eq. ( 1) produces at each time t the quaternion Q(t) that enables to retrieve (with suitable conversion matrixes) the attitude angles without any risk of singularities with unit quaternion.As the quaternion of the satellite with respect to the orbital frame is available (given by the state port), it can be used in an external control loop of the attitude of the satellite for an Earth pointing command by setting to zero its imaginary components.

Solar Arrays Component
For the current release of the library, either ¦xed SAs on the satellite body or an automatic orientation of the SA with respect to the Sun is considered.The component is linked to the Frame component (state port) in order to get the Satellite to Sun vector and other information.
The classic equations of Sun pressure interaction on the SA are used in this component with input coe©cients of re §ection, specular re §ection, and absorption.During eclipses of Sun by the Earth, those interactions are null.
The output of the SAs are the 3D forces due to the Sun pressure, their moments, and the power produced by the solar cells according to the input data of the e©ciency and area of the solar cells.
As for other components, the SA component is written in terms of vector to allow a number of SAs (up to 12, for instance) that are all described in terms of location, orientation of the axis of rotation, and orientation of the canonical normal to the solar cells, size.

Thrusters Component
When activated, this component outputs the thrust vector with respect to the satellite, the moment induced by the thrusters£ forces, and the mass §ow rates (for monopropellant as well as for bipropellant).In case of electric propulsion, the electric power consumption is provided as well.

Reaction Wheels Component
Once activated, this component outputs the angular momentum vector with respect to the satellite as well as the electric power consumption of all RWs considered.This component has been improved in the 4th phase by adding the operational limits in terms of maximum torque and angular rate or rpm limits.

Drag Array Component
This component is a new component of the 4th phase.It outputs the drag force coming from the dynamic pressure in §ight.The Earth atmosphere type is taken from ESA ECSS (European Cooperation for Space Standartization) standard [7] with several cases depending on the Sun activity.

Gravity Gradient Component
This component is also a new component of the 4th phase.Based on the input from the state port coming from the satellite frame component and the inertia matrix [Inertia] of the whole body (including §uids), this component provides the moment vector with respect to the satellite without any forces (pure torque): where R is the distance from CoM to the Earth center; � e Earth is the unit vector from the CoM to the Earth; and � r is the local vector from the CoM.

Tanks Component
The classic relation for Archimedes£ pressure dp = −ργ(z −z 0 ) was already taken into account in the existing Flow1D tank model of the original ESPSS [3] but this was too restrictive.In space, the satellite is rotating around its CoM (and some couples are produced by thrusters and RWs or windmill e¨ect by the solar pressure on the SA).Further particular forces, e. g., originating from thrusters or solar pressure on the SAs, produce acceleration.Therefore, the Archimedes£ pressure needs to take into account the e¨ect of the acceleration and the e¨ect of the rotations even when not collinear to the acceleration.

A new equation for the Archimedes£ pressure
If � a represents the acceleration of the CoM in the inertial frame, m� g all the body forces, and � F the surface forces like the reaction to the gravity, the thrust forces from the thrusters, the Sun pressure force, etc., one gets to The NavierStokes equation is the Newton£s law of momentum ¤ρ� a abs = ρ� g + � f¥ adapted for Newtonian §uids (with a control volume ¦xed in the noninertial space considered) where � a abs is the absolute acceleration of a point (located at � r from the CoM) in an inertial frame and � f are all the surface forces per unit of volume (those are derived from � F above).Once we have � Ÿ being the rotation around CoM, � r and � V rel being the location and velocity of the §uid particle with respect to the noninertial frame, and µ being the viscosity, one can write the instantaneous momentum equation as Considering a static behavior in the noninertial frame, i. e., � V rel = 0, it remains for the point � r: In the literature, it was not possible to ¦nd this kind of simple relation; so, no references can be linked to it.
The pressure di¨erence between two points is the integral of the scalar product � ∇p • � ds: Note that centripetal acceleration (see [8]) Hence, if the density ρ is constant, one gets the new equation: .
To the authors£ knowledge, this compact formulation for the Archimedes£ pressure di¨erence in an accelerating-rotating frame could not be found in the literature.

Case of tanks and components of the Fluid Flow 1D library
Inside tanks, the free surface should be considered for one of the points; hence, the Archimedes pressure at � z b with respect to the point � z free (a single point belonging to the free surface is enough) is given by where � F is the sum of all nonbody forces acting on the CoM, i. e., surface or contact forces: thrust, perturbation forces due to Sun pressure, drag, etc. and also including the reaction force to the weight for the case of simulation on ground.
The free surface of the tank is given by a 3D discretization of the tank using the equation like (2) for the location of the liquid phase until the volume reaches the current given volume of liquid.
The shape of the free surface is a paraboloid having an axis parallel to � Ÿ.Its deviation from � Ÿ depends on the strength of the speci¦c force � F /m.For the tubes, a similar equation is used: .

It is to be highlighted that in this equation
where tau is the time constant.Finally, the equation used for the Archimedes pressure is the following: Those dynamic vectors are set as global for every component of ESPSS, in the existing Fluid Flow 1D library, so that there is no need to add any links between the tubes and the frame for being able to use Eq.(3).

LogicCmdThr Component
This new component allows to perform the active stabilization with thrusters.Its input data are the force, location, and orientation of the thrusters, and some data regarding the stabilization accuracy (angular rate limit and attitude angles limit).
The logic for the stabilization follows the principle of the Dead-zone control with hysteresis: for each axis, the angular rate drives the evolution of the attitude angle.Once a switching limit (that is, a combination of angle and angular rate) is reached, a torque should be executed in order to control the attitude by reversing the sign of the angular rate.This kind of control is sketched in Fig. 2 in the phase space of one axis and without any perturbations: the limit cycle is the loop ABCDA inde¦nitely.
The command of the thrusters is performed with [matCV], a rectangular con¦guration matrix (its columns are the torque components produced by each thruster activated at 100%) according to So, the pseudoinverse of the rectangular con¦guration matrix allows selecting each thruster£s level for a given set of torques.Because the thrusters can only produce thrust in one direction, a vector from the null space of the rectangular matrix [matCM] should be added in order to get only null or positive levels for each thruster: Figure 2 Attitude control Dead-zone with hysteresis (torque = 7.6 • 10 −6 N•m; inertia = 0.0325 kg•m 2 ; and angular rate limit = 0.0057 deg/s)

LimitNb Component
Linked to the previous component, this component limits automatically the number of thrusters ¦red simultaneously and, in case of need, sequences the ¦ring between the thrusters.Its input data are the maximum number of thrusters simultaneously ¦red and the period of sequencing.

Statistics Component
Linked to the previous component, this component is exclusively used for transparently controlling the thrusters (directly from the input port) while keeping trace of the use of di¨erent thrusters.

SATELLITE LIBRARY VALIDATION CHECKS
The validation of the library including §ight dynamic validation, attitude control validation, and tank Archimedes pressure validation is described in [9].

SIMULATION CASES FOR COUPLED PROPULSION SYSTEM AND VEHICLE
In order to highlight the capabilities of the SATELLITE library and, particularly, the coupled Simulation of the Propulsion System and Vehicle, various realistic examples are presented.

Propulsion System of a Generic Satellite with Vehicle Coupling by Induced Acceleration and Angular Rate
The model of the propulsion system considered is shown in Fig. 3.It comprises 4 main tanks, several valves, ¦lters, tees and tubes, and 4 RCTs.
The new features of the satellite library allow the connection of the force vector from each RCT to the Frame component, from each tank to the Frame component in order to set the vehicle mass, acceleration, inertia matrixes, and angular rate.Those variables are taken into account for the Archimedes£ pressure in the tanks with right location of their free surface and into the whole manifold and, ¦nally, at the RCS (Reaction Control System) injector.
The main results of the simulation are shown in Fig. 4 where a thrust pulse of 3 s is performed.During the thrust pulse, the mass decreases according to the tanks depletion (coming from the thrusters consumption through the tubing manifold) and after the pulse, the acceleration cancels while the angular rate around the transverse axis of the vehicle remains constant.

Bipropellant Propulsion System with Vehicle Coupling by Induced Acceleration and Angular Rate
The model of the propulsion system considered is shown in Fig. 5.It comprises two main tanks in parallel for monomethylhydrazine (MMH) and nitrogen tetroxide (NTO), one pressure vessel for the Helium gas, a pressure regulator, several valves, nonreturn valve, tees and tubes, and a main rocket engine.The new features of the satellite library allow the connection of the force vector from the main rocket engine nozzle extension to the Frame component, from each tank to the Frame component in order to set the vehicle mass, acceleration, inertia matrixes, and angular rate (the moving mass of helium is neglected in this model).Those variables are taken into account for the Archimedes£ pressure in the tanks with right location of their free surface and into the whole manifold and, ¦nally, at the main rocket engine injector.
The main results of the simulation are shown in Fig. 6 where a thrust pulse of the main rocket engine lasts for 8 s.
The derivative of three vectors � F , � Ÿ, and − −−− → OCoM used for the §uid §ow 1D components are starting with some delay with respect to the original one given by the frame component of the Satellite library, the time constant being set to 0.1 s.During the thrust pulse, the mass decreases according to the tanks depletion (coming from the rocket engine consumption through the tubing manifold of MMH and NTO) and after the pulse, the acceleration cancels while the angular rate around one of the transverse axes of the vehicle remains constant at a high rate because the main rocket engine thrust is intentionally deviated along X axis by +70 mm relative to the CoM.
Those results show that the new feature regarding the coupling between propulsion and vehicle are well managed and the propulsion system can be simulated with a great realism.

Active Attitude Control with Thruster
The model of the propulsion system considered is shown in Fig. 7.It comprises 4 thrusters for performing thrust in one main direction and with components along a transverse axis for performing a full three axis attitude control.The logic presented in the previous section is used for the assessment of the system, taking into account the fact that only one thruster can be ¦red every second and that a delay of 1 s is needed before using a thruster.
Figure 8 shows a zoom (during 150 s) that the constraints imposed to the propulsion system are well respected during the simulation: according to the value of the switch function of the Dead-zone control with hysteresis, none or one of the thrusters are activated (level between 30% and 100%) for a duration of 1 s and, when needed, thrust pulses are successively performed every second.The overall result of the stabilization within ±0.2 • for the three-axis roll, pitch, yaw is shown in their phase space (Fig. 9b).One can notice that only a few very small overshoots occur for the control of the roll axis, while for pitch and yaw, the attitude goal is ful¦lled entirely on this example.

CONCLUDING REMARKS
The paper has presented in detail the progress made during the 4th phase of the work performed on the ESPSS library within the existing tool EcosimPro .The signi¦cant advance is related to the implementation of the Archimedes pressure for the most general case of noncollinear acceleration and angular rate for the ESPSS §uid components.Thanks to the use of a new very compact equation, this could be performed in a much easier form than in the previous complex formulations.This has been successfully implemented for the components needed for a propulsion system.Two examples of propulsion systems have been presented with highlights the coupling between propulsion vehicle.
Among the progress, new components have been implemented successfully into the Satellite library allowing to simulate with ¤good¥ realism typical spacecraft missions, particularly of interest for propulsion applications as the attitude control for a small satellite with thrusters.

Figure 1
Figure 1 Palette of components for the Satellite library of ESPSS with coupled components in other ESPSS libraries , � r is the distance between the point considered − − → OM and the real current CoM.Because in a dynamic system the three vectors � F /m, � Ÿ, and the vector from the reference origin O to the CoM noted − −−− → OCoM may change every time, a very powerful technique for fast integration of the system dynamic vs. time is to use additionally the derived variables − For example, one uses −−−−→ Omega derived from � Ÿ by The derivative of three vectors � F , � Ÿ, and − −−− → OCoM used for the §uid §ow 1D components are almost identical to the original one given by the Frame component of the Satellite library.

Figure 4
Figure 4 Center of mass, angular rate, acceleration, and mass of the coupled propulsion system with a thrust pulse 3 s

Figure 6
Figure 6 Center of mass, rate, acceleration, mass of the coupled bipropellant propulsion system

Figure 7 Figure 8
Figure 7 Propulsion system for active AOCS

Figure 9
Figure 9 Active stabilization of attitude within 0.2 • for two orbits: (a) orbit and attitude angles; and (b) phase spaces