Satellite design optimization for differential lift and drag applications

Utilizing differential atmospheric forces in the Very Low Earth Orbits (VLEO) regime for the control of the relative motion within a satellite formation is a promising option as any thrusting device has tremendous effects on the mission capacity due to the limited weight and size restrictions of small satellites. One possible approach to increase the available control forces is to reduce the mass of the respective satellites as well as to increase the available surface area. However, satellites of these characteristics suffer from rapid orbital decay and consequently have a reduced service lifetime. Therefore, achieving higher control forces is in contradiction to achieving a minimum orbital decay of the satellites, which currently represents one of the biggest challenges in the VLEO regime. In this work, the geometry of a given reference satellite, a 3UCubeSat, is optimized under the consideration of different surface material properties for differential lift and drag control applications while simultaneously ensuring a sustained VLEO operation. Notably, not only the consideration of sustainability but also the optimization with regard to differential lift is new in literature. It was shown that the advantageous geometries strongly depend on the type of gas-surface interaction and thus, two different final designs, one for each extreme type, are presented. In both cases, improvements in all relevant parameters could be achieved solely via geometry adaptions.


Introduction
In the recent past, interest in operating satellites at much lower altitudes than before, in the so-called VLEO regime (i.e. the entirety of orbits with a mean altitude < 450 km [1]), has increased due to a variety of advantages, which have been summarized by Crisp et al. [1].Similarly, the application of distributed satellite systems, which commonly are made up of several small satellites working together to achieve a shared goal, is nowadays prevalent [2].Due to their stringent volume and mass limitations, alternative solutions to the conventional propulsion methods to exert control forces are of great interest.In the VLEO regime, the use of the atmospheric forces from the prevailing residual atmosphere represents a promising solution.
At the Institute of Space Systems of the University of Stuttgart, this methodology is actively researched since 2018 [3].In the most recent publication, a planning tool for optimal three-dimensional formation flight maneuvers of satellites in VLEO using aerodynamic lift and drag via yaw angle deviations has been presented [4].As in the VLEO regime the large levels of orbital decay represent the major challenge to be overcome for a sustained operation to become reality, the planned trajectory is optimal in a sense that the overall decay during the maneuver is minimized.Thereby, the remaining lifetime of the satellites is maximized and the practicability and sustainability of the methodology significantly increased.Additionally, applying yaw angle deviations allows to simultaneously control the in-and out-of-plane relative motion via differential lift and drag.Throughout the work of Traub et al. [4], however, conventional 3UCubeSats augmented with two solar panels have been assessed.In parallel efforts, optimized shapes for VLEO satellites targeting a minimization of the atmospheric drag force and thus extension of operational lifetime have been developed.Thereby, the satellite geometry has been optimized via a novel 2D profile optimization tool [5].It was shown that the optimized satellite geometries offer pure passive lifetime extensions of up to 46 % compared to a GOCE like reference body.
In this article, which is based on the corresponding author's master's thesis [6] and builds upon said contributions, optimal satellite designs for differential lift and drag applications are presented.To the best of the authors' knowledge, such efforts have never been discussed in the literature.The article elaborates how the research objective leads to conflicting requirements, which need to be prioritised accordingly to arrive at an overall optimal design.All optimization steps are presented in detail and the corresponding considerations are discussed.In future efforts, it is foreseen to combine the optimal maneuver planning and the optimal satellite designs, which would represent a more holistic optimization approach.
The necessary fundamentals are briefly presented in Section 2. In Section 3, the for the calculations and optimization set parameters and constraints as well as the applied optimization steps are introduced.The design optimization for diffuse re-emission is conducted in Section 4, and the similar process for specular reflection can be found in Section 5.

Fundamentals
The forces acting on a satellite are the result of the interchange of momentum between the particles of the prevailing atmosphere and the satellite.Whereas the atmospheric forces attacking at the center of gravity are of interest within this work, the torque induced by the atmosphere was neglected and consigned to the attitude control system of the satellite.The specific drag force f D and specific lift force f L for a satellite of mass m can be calculated as shown in Eq. 1 and Eq.2: where ρ is the prevailing density, C D the drag coefficient, C L the lift coefficient, A re f the reference area, which corresponds to the projected area of the satellite perpendicular to the flow, and v rel the relative velocity to the local atmosphere.In Eq. 2, n denotes the surface normal vector of the surface under consideration.
In order to investigate the sustainability of the different satellite designs, their lifetime was estimated.To ensure comparability with the previous study [5], the same simplified equation for the lifetime of a satellite in a circular orbit was applied [7]: Here, t l is the lifetime of the satellite for a given orientation and no additional use of propulsion.H 0 is the atmospheric scale height, h 0 the base altitude and ρ 0 the density as shown in Tab. 1. µ E is the Earth's gravitational parameter and a = h 0 +R E with R E being the Earth's radius.Although Eq. 3 is very simplified due to the neglected changes in density with time and location, it is well suited for the purpose of this work, since all atmospheric changes have to be eliminated to only account for the influence of the satellite's design.

Aerodynamic coefficients and gas-surface-interaction
In Eq. 1 and Eq. 2, the parameters m, A re f , C D , and C L are characteristics of the satellite's design.Whereas the first two can be directly derived and measured using the satellite or its model, the two aerodynamic coefficients are depending on the surface properties of the satellite and have to be calculated analytically or numerically.The two main parameters influencing the aerodynamic coefficients are the velocity of the reflected Incident Flux Incident Flux Figure 1: Specular reflection (left) and diffuse re-emission (right) [12].
particle, which is depending on the amount of exchanged energy between particle and wall, commonly described using energy accommodation coefficients, and the angle of reflection/reemission, which is additionally depending on the type of gassurface-interaction (GSI).There are two extreme types of GSI; specular reflections, where the angle between surface and reflected particle is equal to the angle between surface and incoming particle, and diffuse re-emissions, with a particle reemission according to a probabilistic velocity and direction distribution (see Fig. 1).
In general, specularly reflecting material are advantageous as they allow a deliberate deflection of the incoming particles.However, specularly reflecting materials are not yet available and remain subject of current research [8,9].Hence, the prevailing type of GSI in VLEO on the currently used materials is a diffuse or quasi-diffuse gas-surface interaction [10,11].

Accommodation coefficients
The thermal energy accommodation coefficient is a measure for the adaption of energy from the impinging particles E i to the energy E w , which the particles would have after reaching equilibrium with the satellite's wall temperature.E r corresponds to the actual energy of the reflected particles.Similarly, the momentum exchange between particle and surface is commonly described by the tangential momentum accommodation coefficient σ t and the normal momentum accommodation coefficient σ n as follows [13]: Here, τ i is the tangential momentum carried to the surface by the incident particle and τ r is the tangential momentum carried away from the surface by the reflected particle.The tangential momentum carried away from the surface by a diffusely reflected particle after reaching thermal equilibrium with the wall is ascribed to τ w and is per definition equal to zero.The normal momentum carried to the surface by the incident particle and normal momentum carried away from the surface by the reflected particle are p i and p r respectively, and p w is the normal momentum carried away from the surface by a diffusely reflected particle after reaching thermal equilibrium with the wall.

Differential lift and drag control method
In order to change the formation geometry, the method of differential lift and drag is generally understood to intentionally create differences in the acting aerodynamic forces [3], e.g. by rotating only one of the two satellites of the formation (see Fig. 2).Whereas differential drag effects the relative movement within the orbital plane, differential lift effects the out-of-plane relative movement.The maximum achievable differential specific forces are as shown in Fig. 3b.Within this work, it is assumed that all satellites within the formation share the same design.The following flight configurations 1 are defined, of which Fig. 3 gives an overview of: 1 In this article, the term "(flight) configuration" is used to describe the ori-• A nominal flight configuration experiencing as little drag as possible and no/negligible lift ( f D,min , see Fig. 3a left), • A maximum drag configuration for experiencing the highest possible drag ( f D,max , see Fig. 3a right), • A maximum lift configuration for experiencing the highest possible lift ( f L,max , see Fig. 3b).
With the aim of optimizing the satellite's design for a sustainable use of differential lift and drag control methods in VLEO, the specific optimization goals can be expressed as follows: In order to derive design recommendations, the different design characteristics are prioritized as they contribute to contradictory optimization goals.Considering that aerodynamic force control methods are not the means of choice for a time-critical maneuver due to the comparatively low absolute values of specific forces, the overall priority is to decrease the orbital decay.It can be assumed that the satellite spends most of its service life in the nominal flight configuration, and hence as little unwanted aerodynamic forces as possible are desired.Due to the prevailing condition of the contamination with atomic oxygen on currently used materials and a therefore diffuse gas-surfaceinteraction with generally little C L values, increasing the lift force is set to be the second priority within this work.
To summarize, the priorities underlying all decisions within this work are set as follows: with ↑ indicating an increase.However, these priorities are based on the considerations mentioned above, and naturally may differ for individual missions.

Optimization approach
The framework conditions set to pursue the defined goals alongside the utilized tools and optimization steps taken are presented in the following sub-sections.

Boundary conditions of the optimization
In this section, the chosen atmospheric conditions for the determination of the forces, as well as the chosen reference satellite and the derived geometry constraints for the design variations are briefly presented.entation of a satellite, which can be further defined by the angle between the satellite's longitudinal axis and the orbital plane (angle of sideslip -AOS) as well as the angle between the longitudinal axis and the plane formed by v rel and the orbit's normal vector (angle of attack -AOA).

Atmospheric conditions
The atmospheric conditions are set constant for all considerations within this work in order to ascribe changes in the specific forces to the design variations.Therefore, moderate space weather with a 10.7 cm flux F 10.7 of 140 sfu and a geomagnetic activity index of A p = 15 were chosen 2 .The data shown in Tab. 2 was obtained by averaging the output of the NRLMSISE-00 model [15] of one day per month of 2004 over the year for the chosen altitude of 350 km.Additionally, the satellite's wall temperature is set to T w = 300 K as commonly assumed in the regarding literature [14,16].

Reference satellite
All simulation results are related to a reference satellite under the same conditions in order to evaluate different satellite designs.Due to its universal applicability for various payloads and the size suitable for satellite formations, a 3UCubeSat with one panel each on the top and bottom is used as a reference 2 Detailed information about solar proxies and indices for space weather description can be found in the work of Doornbos [14].
Diffuse Re-emission α T,1 = 1.00 satellite (see Fig. 4).For compliance with the commonly utilized CubeSats [17], the mass of the satellite is set to be 5 kg.The reference satellite has a maximum length l max and height h max of 0.3 m each as well as a main body volume V MB of 0.003 m 3 with main body height and width of both 0.1 m.The width of the panel is equal to 0.003 m.
The dependence of the atmospheric forces on the surface properties of the satellite is taken into account by considering different energy accommodation coefficients α T and corresponding gas-surface-interaction models (GSIM).An overview of the applied surface properties is given in Tab. 3. The chosen gradation steps ∆ 1 and ∆ 2 are derived from the traditional surface materials used for satellites.α T,1 = 1.00 represents the worst case scenario for a deliberate deflection of the particle.The second value, α T,2 = 0.91 represents a value for the traditional surface materials [10,11] and α T,3 = 0.70 represents improved material properties for diffusely re-emitting materials [18].Therefore, ∆ 1 = 0.09 and ∆ 2 = 0.21, which is equally applied on the energy accommodation coefficients for the consideration of specular reflection as shown in Tab. 3. Additionally, this choice of α T allows comparison to other literature with similar gradation steps [19].

Geometry constraints
For comparability of the results within this work, all design variations have to comply with the following constraints regarding the geometry: • Constant volume: Main body volume remains constant.
• Height and length limitations: Maximum height and length of the reference satellite must not be exceeded.
Due to these geometric constraints, the mass of the satellite can be considered constant and is equal to 5 kg according to the reference satellite.Thus, the only variables in Eq. 1 and Eq. 2 are the aerodynamic coefficients and the reference area, and changes in the forces can be directly attributed to the design changes.

Utilized tools
The utilized open-source tools as well as an extension for a Matlab 2D profile optimization tool for the derivation of satellite geometries are introduced in the following sections.

Determining the atmospheric forces
Since numerous investigations of the influences of different design parameters are necessary, the calculation time for the aerodynamic forces per satellite model plays an important role.Therefore, the Matlab toolkit "ADBSat", which was developed at the University of Manchester [21], was used to study  the effects of different designs on the aerodynamic forces.It represents a means of calculating the acting forces with little computational effort.However, this method is not suitable for concave surfaces, since multiple reflections of the incoming particles cannot be taken into account.For selected and promising geometries as well as geometries causing multiple reflections, their actual performance was determined using the more computationally intensive but proven Direct Simulation Monte Carlo (DSMC) method, which is available within the gas-kinetic simulation framework "PICLas", developed at the University of Stuttgart [22].For the PICLas simulations, the necessary variable hard sphere parameters for the species listed in Tab. 2 were taken from Bird [23].A comparison of ADBSat and PICLas results can be found in Appendix A.

2D profile optimization
For optimizing and obtaining the satellite models, an already existing Matlab 2D optimization tool for reduced drag [5] was utilized and adapted in the scope of this work.This tool was developed for an analytical determination of a 2D profile with minimum drag for a given length and volume.A preceding determination of the type of volume derivation of the profile is necessary in order to comply with the given volume restrictions during the profile optimization.In the given optimization tool, three different volume derivation options are possible as shown on the left in Fig. 5. Whereas geometry option A is obtained by rotating the profile around the x-axis (longitudinal axis), geom-etry option B is obtained by mirroring the profile at the x-axis and then extruding symmetrically with the value of the largest cross-section width.Thereby, geometry option B has a squared rear side by default.Geometry option C is obtained by intersecting the extruded profile of option B with a cylinder with the given length and a radius equal to the largest cross-section width of the optimized profile.Additionally, the application of a tail geometry is possible, in which the height of the 2D profile is reduced again from a certain longitudinal position.Each individual profile optimization complies with the given volume restriction of V MB = 0.003 m 3 .Within this work, the 2D profile optimization tool was extended as presented in the following: Additional geometry options: The extrude option B was modified to consider the following two additional user inputs (see Fig. 5 right): 1. Extrusion height h mb : For the application of the optimization tool on a given structure such as the panels of this work's reference satellite, the height of the extrusion can be set as an input parameter.2. Nose width w N : Whereas it might be aerodynamically advantageous to have a sharp or pointy nose geometry, it could pose a problem for accommodating the payload and manufacturing of the satellite.Hence a definition of a nose width during the optimization process was implemented.
Lift optimization: Within the 2D optimization tool, a derivation of Sentman's equation was used [5].Since the Matlab tool is aiming on optimizing a 2D profile, it followed for Sentman's equation that the two direction cosine η, t = 0. Considering the definition of the area element on the profile as shown in Fig. 6, the direction cosine between the local x and y axis and the direction of the drag force is set as follows in the original tool: with By setting the two direction cosine according to the desired force direction of the lift force (see Fig. 6) to an optimization of a given profile with regard to increasing the lift force experienced by the profile is now possible as well.It should be noted, that therefore the sign of the original optimization function has to be inverted, since originally a decrease in the acting force was desired.

Optimization steps
The optimization process presented in this subsection is equally executed for the assumption of diffuse re-emission and specular reflection.In a first step, the design is varied with the aim of optimizing the design for solely differential drag or solely differential lift, in order to investigate the respective important design characteristics.In a second step, the design characteristics are combined for achieving a best possible trade-off according to the set priorities.
The optimization process for deriving a design optimal for differential drag is divided into the following steps: Drag 1a -Decreasing minimum drag (main body): The main body is optimized for minimum drag in the nominal flight configuration utilizing the extended Matlab 2D profile optimization tool (EPOT) while maintaining the maximum length of 30 cm and the main body volume of 0.003 m 3 .

Drag 1b -Decreasing minimum drag (panel):
The panel is optimized for minimal drag in the nominal flight configuration using the EPOT, considering a given extrusion height of 30 cm.The trivial solution, omitting the panel, is not considered within this work as it would greatly reduce the achievable control forces.Thus, the size of the panel remains l max • h max .
Drag 2 -Increasing maximum drag: The loss in maximum drag due to the first and second optimization steps has to be kept as small as possible.Considering the optimum angle for the maximum of C D (see Fig. 7) and the set geometry restrictions, the following approaches can be derived: • Increasing number of area elements with C D,max : For an AOS of 90°, A re f cannot be further increased due to the given l max and h max and thus, the overall C D of the whole body needs to be maximized.A possible solution for increasing C D is implementing the largest perpendicular area to the flow in the maximum drag configuration.Within this work, this was achieved by changing the panel position and approximating the satellite shape with vertical and horizontal area elements (see Figs. 8a-8c).
• Multiple reflections: For a consideration of the whole geometry and not only single area elements (as within ADBSat), a geometry promoting multiple reflections poses another possible approach to increase the drag.This approach was implemented in the design by overlaying the optimized 2D profile with a zigzag curve for particle-trapping concave areas (see Figs. 8d-8e).
Opposite to the design optimization process for differential drag, where a high difference of f D,min and f D,max but generally lower absolute values of f D,min and f D,max are desired, the goal for optimizing the design for differential lift only consists of increasing f L,max .The process for optimizing the design for differential lift was divided into the two steps: Lift 1a -Increasing maximum lift: For achieving higher lift forces, the profile obtained using the EPOT lift optimization function (see Fig. 8f and 8g) was implemented in the design.
Lift 1b -Eliminating lift reducing areas: When considering the reference satellite with its frontal area it is apparent, that for an AOS ≈ 45°(i.e.maximum C L value for the panels, see Fig. 7), the frontal area and the velocity vector form an angle ≈ −45°.Therefore, the frontal area experiences lift as well, but in the opposite direction as desired.
Hence, when dealing with different profiles than the optimized one, areas experiencing lift in the opposite direction than desired must be decreased.
Each design was generally tested in a variety of AOA and AOS in order to determine the orientation achieving the highest drag and lift force, referred to as the maximum drag and maximum lift configuration, respectively.The same was applied to the reference satellite.Therefore, in order to compare the optimization in drag and lift forces, the respective configurations with the maximum achievable values are compared to each other.

Design optimization -diffuse re-emission
The results of a design variation are generally related to the reference satellite with the same surface properties, if not stated otherwise.The gas-surface-interaction model applied for the ADBSat calculations within this section is Sentman's model [24].

Differential drag optimization
Within the design considerations for the use of differential drag, the influence of general shape and tail geometry, surface structures and panel position has been analyzed, and the results are presented in the following sections.

General shape and tail geometry
As investigated by Hild [20] and Walsh [25], it can be advantageous for the overall lifetime of the satellite to apply a tail geometry on slender bodies.In the following, the effect of a tail geometry for the three different volume derivation options A, B and C is investigated in order to find the optimum combination for an increased lifetime.At the same time, their influence on the achievable maximum drag and thereby the overall differential drag is assessed.The influence of the tail length on the lifetime in nominal flight configuration and experienced drag in maximum drag configuration as well as the resulting differential drag for different α T can be seen in Fig. 9, Fig. 10 and Fig. 11.The data shown was obtained using ADBSat.
It is apparent, that the optimum tail length generally depends on the volume derivation option and the surface properties.Whereas the geometry option A is promising for high α T , the increase in lifetime for geometry options B and C generally increases with decreasing α T .The three following geometry options and tail lengths are optimal concerning the lifetime for the given cases of surface properties: for α T = 1.00 and α T = 0.91 geometry A with 23 % tail length and for α T = 0.70 geometry A with 25 % tail length.However, as shown in Fig. 10, geometry option B experiences the higher drag forces in maximum drag configuration as expected due to its vertical side surfaces, which are perpendicular to the flow in maximum drag configuration and thus correspond to the highest C D values.Even though geometry option A can reach a lifetime increase of up to 13%, the negative effect on the maximum achievable drag by the rounded side surfaces leads to an overall loss in differential drag compared to the reference satellite as visible in Fig. 11.The increase in lifetime with geometry option B is smaller than for A and C, but due to the generally higher maximum drag, an increase in differential drag compared to the reference satellite is possible.However, as the priority of this work is to increase the lifetime, options for increasing the maximum drag while maintaining the geometry characteristics of geometry option A were investigated.

Surface structures
Since the main body optimization with the EPOT led to side areas with different orientations to the flow and therefore smaller C D values, the maximum drag experienced is smaller compared to the reference satellite with completely vertical side surfaces.In order to increase the C D values of the side surfaces in the maximum drag configuration while maintaining the overall optimized profile shape obtained by the EPOT, a rasterized cross section was added to the optimized profile in order to approximate the curved surface by horizontal and vertical surface elements.Therefore, the frontal circular area, which was swept along the optimized profile, is rasterized using the midpoint circle algorithm 3 to two different raster step sizes as shown in Fig. 8a and Fig. 8b.Here, Raster 1 has a step size of r f ront /14.3 , with r f ront being the radius of the frontal area, which is a result of the EPOT for given volume restriction, l max , α T and geometry option.With this step size, the geometry constraint of constant main body volume remains met.Raster 2 has a step size of r f ront /10.Another possible means to increase the maximum drag is to use surface structures promoting multiple reflections.Therefore, the geometries shown in Fig. 8d and Fig. 8e, which are derived by overlaying the optimized profile for geometry option A with a zigzag-curve with an extrema distance d ex of 5 mm (design MultRef 1) and 10 mm (design MultRef 2) as well as extrema of ±5 mm, were investigated further.
The results for the above presented satellite designs are listed in Tab. 4. Here, the lifetime was estimated in the nominal flight configuration and the maximum drag in its respective configuration.The data shown was obtained using PICLas and is referring to the data from the reference satellite obtained using PICLas as well.It can be seen, that the loss in maximum    drag is smaller for the design with additional vertical side surfaces than the designs promoting multiple reflections.The force distribution of the different designs compared to the reference satellite is shown in Fig. 12 for α T = 1.00 and obtained using PICLas.As the area element's C D value is the highest for the vertical surfaces, these surfaces experience the highest drag forces.However, the effect of multiple reflections for the tested design is lower than expected and the additional decrease in vertical side surface area leads to an overall loss in experienced drag.The expected effect of multiple reflections increasing the drag in maximum drag configuration for the tested design did not appear.However, the results for the minimum drag in nominal flight configuration differed up to 25 % compared to the results obtained using ADBSat.Therefore, the effect of multiple reflections did occur, but not in the desired way.It can be assumed, that the utilization of structures promoting multiple reflections for an increased maximum drag force is generally possible, although the tested design is not suitable in this respect.The rasterized cross section on the other side did not lead to a loss in the generally increased lifetime by the optimized profile, and additionally increased the effective side area and thereby the differential drag.In the case of a chosen geometry option C, the rasterized cross section can be applied as well for increasing the vertical side surfaces (see Fig. 8c).

Panel position
Taking up the theoretical approach for further increasing the maximum drag, changing the panel position as shown in Fig. 13 poses another possibility to increase the surface perpendicular to the flow, which is marked in orange.This parameter study of the panel position was performed on a main body of geometry option A, since the geometry option C has limited width to move the panel due to the sharp nose, and since changing the panel position on geometry option B does not change the amount of surface elements perpendicular to the flow.The influence of the panel position on lifetime, maximum drag and differential drag can be seen in Fig. 14 and was estimated using ADBSat, by gradually moving the panel from the centered position to the edge of the frontal circular area.It can be seen, that as expected, the maximum drag can be increased, but simultaneously the lifetime decreases due to the greater A re f in nominal flight configuration.Hence, changing the panel position is not a recommended design variation for the use of differential drag control methods in the given case.

Differential lift optimization
The lift values compared and shown in this section are referring to the optimum angle of attack for the respective satellite geometry.Hence, the compared lift values represent the maximum possible lift, but may be experienced at different AOSs.
Optimizing the profile of the satellite using the EPOT's lift optimization function led to a simple geometry as presented in Fig. 8f for a given main body height.However, here the side surface of the main body and of the panel form a different angle with the macroscopic velocity vector.Hence, the geometry can be further optimized by increasing the main body height to the original height of the panel as shown in Fig. 8g.In order to validate the EPOT's output, other geometries (e.g. the designs built for the investigations within the previous section) have been evaluated as well, but none of them were able to increase the lift such as design Lift 2, when compared to the reference satellite.However, since design Lift 2 contains a sharp nose geometry, possibly posing a problem towards the accommodation of the payload, the consequences of a mandatory frontal area were investigated as well as shown in Fig. 15.In order to obtain the ADBSat results for the change in lift force, the utilized satellite models were obtained extruding the with the EPOT optimized profile considering the respective main body heights h MB and nose widths w N .It is apparent, that for a given nose width, the overall frontal area increases with increasing main body height, and thus the lift opposite to the desired lift direction increases as well.Hence, the greater the given nose width, the smaller the optimum main body height and achievable increase in lift force.In order to achieve a high increase in lift force, a frontal nose area should be avoided.If not possible, the main body height has to be chosen according to the applied nose width for a maximum increase in lift force.

Discussion on material as design factor
Within this section, the influence of improved surface properties under the assumption of diffuse re-emission is discussed.In Fig. 16, the data obtained for the considerations of a frontal nose area, geometry option B, and a varying main body height as presented in Section 4.2 for different α T is shown with respect to the data from the reference satellite with α T = 1.00.Here, the comparison to the reference satellite of α T = 1.00 was chosen in order to visualize the effect of varying α T and varying the geometry.It can be seen, that changing the reference satellite's α T from 1.00 to 0.91 increases the lift more than optimizing the geometry with h N = 0 for α T = 1.00.So not only does the change from diffuse re-emission to specular reflection enable greater atmospheric forces than a design variation, but also a reduction of α T for diffusely re-emitting materials.Especially for high α T , adapting the material is therefore a more powerful design variation than optimizing the geometry with regard to increased control forces.As Fig. 16 shows, the effect of a change in geometry (e.g. from h N = 0 cm to h N = 2 cm) increases with decreasing α T .

Performance of the optimized designs
The final design recommendation for the assumption of a diffuse re-emission is presented in the following and its actual performance was evaluated using PICLas.After investigating several design derivations of the reference satellite, a satellite geometry with a revolved lifetimeoptimized 2D profile including an additionally added rasterized cross-section swept along the 2D profile as a main body, and an extruded lifetime-optimized 2D profile for the panel was shown to be advantageous for decreased minimum drag and increased differential drag.The design (Raster 1) is shown in Fig. 8a with 2D profiles optimized for α T = 1.00.
To increase the experienced lift force in the desired direction, it was shown to be advantageous to increase the side surface area experiencing lift, as well as decreasing the frontal area, which otherwise experiences a lift force in the opposite direction than desired.The design recommendation (Lift 2), which applies for all the tested α T within the diffuse re-emission section, is shown in Fig. 8g.

Differential drag and lift
The recommendation for differential drag and the recommendation for differential lift do not share many geometric characteristics.Since a sustainable operation in VLEO is the goal, design Lift 2 is not suitable for the general use for the relative motion control by means of the atmospheric forces.Only for α T = 0.70, no loss in lifetime compared to the reference satellite with α T = 0.70 is achieved with design Lift 2. However, another solution for materials with higher α T is of interest.
Raster 1 was shown to be the best trade-off, since the vertical side surfaces for increasing the maximum drag simultaneously have a positive effect on increasing the maximum lift force and therefore, an improvement in all critical parameters was achieved.The flow-fields in all three flight configurations for the reference satellite and design Raster 1 are given in Figs.17a-17f.It can be seen, that the 2D profile optimization of main body and panel reduces the shock occurring in Fig. 17b for x < 0 m.The reduction of opposed lift experienced by the frontal area is visible in Fig. 17f, where the shock preceding the frontal area of the recommended design is smaller compared to the reference satellite.The respective specific forces and lifetimes for the designs Raster 1 and Lift 2 are given in Tab. 5. Due to the symmetry, the minimum lift is equal to zero.Since the altitude loss per maneuver plays an important role, the value for drag in the maximum lift configuration f D,L is given as well.

Design optimization -specular reflection
The satellite models within this section were generally evaluated with the Cercignani-Lampis-Lord (CLL) model of ADB-Sat.As the required input parameters are α n and σ t and not α T , they were derived from the chosen α T given in Tab.Table 5: Data for the reference satellite, design Raster 1 and Lift 2, α T = 1.00.Data obtained using PICLas, with * marked data obtained using ADBSat.

Differential drag optimization
Within this section, the influence of general shape, tail geometry, main body height, and frontal area on the experienced drag forces is presented.Sections 5.1.2and 5.1.3differ from the respective sections under the design optimization process of diffuse re-emission, as they are linked to the results of the investigation of general shape and tail geometry.

General shape and tail geometry
In a first step, the best possible volume derivation of the optimized 2D profile and the appropriate tail length (if any) were investigated.An overview of the lifetime in the nominal flight configuration influenced by the three geometry derivation options A, B and C and the tail length is given in Fig. 18.The data shown was obtained using ADBSat.For very low energy accommodation coefficients, the geometry option B is the most advantageous for a high increase in lifetime.With increasing energy accommodation coefficient, geometry option C is the more promising volume derivation option.From Fig. 18a and Fig. 18b it can be seen, that for the mentioned advantageous volume derivation options the design with no tail geometry is the most promising design respectively.

Main body height
The volume derivation option B is promising regarding a long lifetime of the satellite.Since the original 2D optimization tool does not investigate different main body heights for option B, their influence was investigated further.Figure 19 gives an overview of the lifetime in nominal flight configuration, maximum drag in maximum drag configuration and differential drag depending on the main body height.The profiles used were obtained utilizing the EPOT with a given extrusion height, which is set equal to the main body height h MB .The data shown was obtained using ADBSat.It can be seen, that especially for very low energy accommodation coefficients, a high main body height has great potential for increasing the overall lifetime.Additionally, an increased main body height is advantageous for increasing the experienced drag in maximum drag configuration and thus also advantageous for a high differential drag force.The with regard to a long lifetime and high differential drag optimized satellite, obtained using geometry option B, is depicted as design SpecOpt in Fig. 8h.1.9 1.95

Consequences of practical considerations
As aerodynamically promising design SpecOpt is, the practicability however of accommodating the payload can be lowered due to the sharp nose.A more practicable solution would be to include a frontal area.Therefore, the following nose geometries were considered: a constant nose width of w N = 3 cm, a constant nose area equal to a third of the reference satellite's frontal area A N,0 , a constant nose area equal to the half of A N,0 as well as a constant nose area equal to two thirds of A N,0 .An exemplary depiction of the tested satellite geometries can be found as design SpecPrac in Fig. 8i.
The dependence of lifetime in nominal flight configuration and differential drag for the considered nose geometries and different α T is shown in Fig. 20.The utilized satellite models were derived using the EPOT considering h MB and w N , and the data shown was obtained using ADBSat.It can be seen that for a given w N regardless of h MB , a smaller main body height is more advantageous regarding the lifetime, as an increase in h MB simultaneously increases the frontal area and therefore the drag.If the frontal area remains constant for different h MB , i.e. a decreasing nose width w N for increasing main body height h MB , the influence of α n and σ t as well as the achievable improvement in lifetime and differential drag decreases with increasing frontal area.Therefore it is recommended to choose a main body height as small as possible for a mandatory given nose width.However, if the value of the nose width is not of direct importance, but a frontal area is required, it is recommended to keep the nose width as small as possible.As visible when comparing the scales of the y-axes of Fig. 19a and Fig. 20a, any type of frontal area leads to a significant loss in lifetime and thereby also reduces the achievable increase in differential drag compared to the reference satellite.

Differential lift optimization
The suggested profile by the 2D optimization tool is equal to the one as derived in Section 4.2 and thus, the influence of different frontal areas was investigated again for the respective α T .Therefore, the satellite models from Section 4.2 were evaluated for α T = [0.00;0.09; 0.30] using ADBSat.Figure 21 shows, that whereas the absolute value of specific lift force also depends on the energy accommodation coefficient, no dependence of the percentage increase in lift with respect to the reference satellite can be seen.Additionally it is visible that the optimal main body height is depending on the nose width similarly to the equivalent investigation results in the diffuse re-emission section.For a high increase in experienced lift force due to solely geometry design factors, a frontal area should be avoided or, if not possible, the nose width should be kept as small as possible.

Performance of the optimized designs
The final design recommendation for the assumption of a specular reflection is presented in the following and its actual performance was evaluated using PICLas.

Differential drag and lift -separate consideration
Especially the extruded optimized 2D profile with an extrusion height equal to the main body height of 30 cm was shown to drastically increase the lifetime and thus the differential drag for low α T .The resulting design (SpecOpt) is shown in Fig. 8h with a 2D profile optimized for α T = 0.00.This design is presented to increase the understanding of advantageous design characteristics prior to deriving a realistic design recommendation.
Similar to the case of diffuse re-emission, it was shown to be advantageous to increase the total side surface area by inclination, as well as eliminating a frontal area.Regardless of the tested GSIM and α T , the design optimal for differential lift is given in Fig. 8g, as presented in Section 4.2.

Differential drag and lift
Opposite to the designs within the section for diffuse reemission, the two previous satellite design recommendations share many geometry characteristics.The only difference between these two is the slightly curved profile, which is the result of the 2D optimization tool with regard to an increased lifetime.As design SpecOpt has a higher increase in lifetime (first priority) compared to design Lift 2, but almost the same increase in lift (second priority), design SpecOpt is also well suited for the use of both, differential lift and drag control methods.However, while design SpecOpt is the theoretical optimal design for the use in VLEO in a satellite formation, it not only poses difficulties for the accommodation of the payload, but also for other factors such as manufacturing or attitude control.Hence, a more practicable but less optimal solution is presented and discussed.A minimum required frontal area is assumed with a nose width of w N = 3 cm4 .In order to keep the loss in lifetime as small as possible, the main body height should remain as   small as possible as well and therefore it is set equal to the reference satellite with h MB = 10 cm.The new design (SpecPrac) is depicted in Fig. 8i with a 2D profile optimized for α T = 0.00.It was obtained using the EPOT geometry B under consideration of given w N and h MB .The panel is optimized for an increased lifetime as well using geometry option B and a given total height of h max = 30 cm.The according flow-fields for the reference and the design SpecPrac are given in Fig. 22.In Tab. 6, the corresponding data is given.The lifetime decreasing effect of a frontal area in the case of the practicable design recommendation can be seen in Fig. 22b, where the deflection of the impinging particles by 180°leads to an accumulation of particles ahead of the satellite.However, the reduction of the total frontal area compared to the reference led to an increase in lifetime.In Fig. 22d the particles are reflected in a wider range with an (in the picture) downward shift due to the increased curvature of the main body's side surface.The overlapping region of particles reflected by the panels and particles reflected by the rearward part of the side surface can be seen in the range around −0.5 m < y < 0.3 m.In this case, the non-symmetrical deflection in the maximum drag configuration leads to a lift force, which is comparatively small, but nonetheless can build up to a considerable effect over time.As the main body's frontal surface area is smaller than the reference's, it was possible to increase the lift force (see Fig. 22f).

Conclusion
In this work, the design of a given reference satellite has been optimized with regard to the use of differential lift and drag control methods and a sustainable application in orbit.The influence of design variations was generally evaluated using ADB-Sat, and the performance of the promising designs was verified using the DSMC code PICLas.
For materials with diffuse re-emission, the design optimized for differential drag and the design optimized for differential lift did not share many geometric characteristics.Hence, when deriving a design optimized for both differential lift and drag, a priority of either lift or drag has to be set.However, for all given surface properties, an optimization with regard to the use of differential lift and drag control methods was possible with the primary objective of increasing the lifetime.Even for the case of complete accommodation α T = 1.00, differential drag and lift could be increased by 2 % and 3 % respectively, while at the same time increasing the lifetime under the set conditions by 12 %.Despite the possible increase of differential forces, the differential lift remains two orders of magnitude smaller compared to the achievable differential drag.Hence, a differential lift maneuver necessary for influencing the out-of-plane motion cannot be performed without a significant loss in altitude.Generally, the design variations had more effect for specular materials than for materials reflecting particles diffusely.However, the theoretical optimal design for a satellite with specular materials entails little practicability.A more practicable solution was presented, which still achieved an increased lifetime of 200 %, increased differential drag by 4 % and differential lift by 1 % compared to the reference satellite with the same surface properties.The same optimized design compared to the reference satellite with a diffuse re-emitting material achieved an increase of 172 % in lifetime, 116 % in differential drag and 2133 % in differential lift.Therefore, the more promising approach for realizing relative motion control in a satellite formation by differential lift and drag for high α T is an improvement of the currently available surface properties, as only for specularly reflecting materials the differential lift and drag are of same order of magnitude.
However, the optimized 3D designs presented and discussed throughout the article are derived from optimal 2D profiles and thus, the direct optimization of the 3D bodies for possibly further increasing the performance is currently still pending.Additionally, the combination of different surface properties in one satellite design was not considered and represents a possible continuation of this work.Applying a diffusely reflecting frontal area and specular side surfaces for a combination of sustainable application in orbit and increased forces for the use in relative motion control within the formation is one example of the possibilities for mixed surface properties.The improvement of the materials is the most influential design parameter in order to make the sustainable implementation of thrust-free control methods possible and represents an important step towards future space systems.Hence, advancing material research is of great interest for the application of VLEO satellites and satellite formations.

Appendix A. Comparison of ADBSat and PICLas results
In this section, the results obtained using ADBSat and PI-CLas for the recommended designs presented in the previous sections are compared.It can be seen, that the results of ADBSat generally comply with the ones obtained using PICLas.The biggest difference is visible for design MultRef 1 and 2 for the minimum drag.However, design MultRef 1 and 2 are the designs created in order to promote multiple reflections.Here, the minimum drag obtained using PICLas is 25 % higher for design MultRef 1 and 12 % higher for design MultRef 2 than estimated with ADBSat.If only geometries suitable for the panel method are regarded, the values for the specific forces obtained using PICLas are generally smaller than the ones calculated with ADBSat.However, the biggest difference is only a 4 % smaller maximum drag of PICLas compared to ADBSat.The generally smaller specific forces within PICLas can be explained due to small differences in the edges of the models.The surface mesh used within ADBSat is better suited to model complex surfaces as the triangular shape of the mesh elements allows more narrow and sharp edges.On the other side, sharp edges or small radii within the mesh for PICLas can pose a problem for the 3D mesh generation.Here, rounded and sharp edges are approximated with a greater radius, as a smaller radius increases the cells and thereby the necessary amount of simulated particles and thus simulation time.Altogether it can be said, that for geometries without multiple reflections the results of ADBSat and PICLas do not differ more than 4 % and the differences can be traced back to differences in the mesh.

Appendix B. Derivation of the ADBSat input parameters
For the application of the CLL model in ADBSat, the normal energy accommodation coefficient α n and the tangential momentum accommodation coefficient σ t are required and therefore are derived from the given α T in Tab. 3. The following derivation is based on transformations of Sentman's equation by Hild [20] and the consideration of a 2D profile.
Since the kinetic energy E = 1 2 mv 2 is proportional to the squared momentum I 2 = m 2 v 2 under the assumption of constant mass, Eq. 4 can be also expressed as with the momentum I consisting of its normal and tangential components p and τ The normal and tangenatial energy accommodation coefficient α n and α t hence can be defined as The parameter τ i and p i as properties of the impinging particle depend on the properties of the atmosphere and are proportional to the following terms [20]: where T w is the wall temperature of the satellite and T i is the temperature of the incident particles [20].Finally, the momentum components of the reflected particle have to be expressed.Therefore, Eq.B.1 can be transformed into Inserting the above derived parameters into the following equations from Hild for τ r and p r [20] τ the tangential and normal momentum carried away from the wall by the particle can be described using the parameters α T and g: and thus, α n and σ t can be obtained utilizing Eq. 5 and Eq.B.3 for given α T and g.However, it should be noted that τ i , p i and p w are determined using Sentman's equation in a 2D simplified version, which is valid for one area element with a given orientation regarding the velocity vector.Hence, the values for τ i , p i and p w as well as τ r and p r change with different orientation of the considered area element.Assuming a constant energy and momentum accommodation over the satellite, the respective parameters may still be used on one representative area element in order to derive the normal energy accommodation coefficient and tangential momentum accommodation coefficient.The derived input parameters for the use of the CLL model according to the given α T in Tab. 3 can be found in Tab.B.7 and were used for all the CLL ADBSat calculations within this work.

Figure 2 :Figure 3 :
Figure 2: Depiction of a maximum differential drag force.flight direction vector of lift force vector of drag force

Figure 6 :
Figure 6: Definition of the 2D optimization tool function values according to Hild [20].

Figure 7 :
Figure 7: Aerodynamic coefficients for a flat plate depending on the angle ψ between flow and surface normal according to Sentman's model.

Figure 11 :
Figure 11: Dependence of differential drag on general geometry and tail length.

Figure 12 :Figure 13 :
Figure 12: Force distribution shown in head on view in maximum drag configuration for α T = 1.00.

Figure 14 :
Figure 14: Influence of the panel position on lifetime, maximum drag and differential drag referred to the panel position at x = 0, α T = 0.91.

Figure 15 :
Figure 15: Influence of nose width w N and main body height h MB on the lift force, α T = 1.00

Figure 16 :
Figure 16: Dependence of lift increase on the main body height and surface properties compared to the reference satellite with α T = 1.00.

Figure 17 :
Figure 17: Simulated particle density n in the flow-field around the reference satellite and the recommended design Raster 1 for α T = 1.00.

Figure 18 :
Figure 18: Increase in lifetime ∆t L depending on general shape and tail length l T .

Figure 19 :
Figure 19: Increase in lifetime, drag, and differential drag depending on the main body height h MB .

Figure 20 :Figure 21 :
Figure 20: Lifetime and differential drag depending on the main body height for different frontal areas.

Figure 8
gives an overview over the designs whose results are listed in Fig. A.23.

Figure 22 :Figure A. 23 :
Figure 22: Simulated particle density n in the flow-field around the optimal design SpecOpt and the more practicable version SpecPrac for α T = 0.00.

Table 2 :
Atmospheric data used for all calculations in this work: temperature of the impinging particle T i , molecular mass M, molecular speed ratio s, and particle density n i of species i.

Table 3 :
Gradation of the surface properties.

Table 4 :
Data for the tested designs with surface structures, α T = 1.00.
+ erf(γs)).(B.6)τ w is per definition equal to zero, and the normal momentum p w carried away after reaching thermal equilibrium with the wall is proportional to

Table B .
7: Derived CLL input parameters according to the specified α T gradation in Tab. 3.