Document [original]

Vol.:(0123456789)

1 3

CEAS Space Journal (2020) 12:397–410

https://doi.org/10.1007/s12567-020-00308-6

ORIGINAL PAPER

Orbit deployment anddrag control strategy forformation flight

whileminimizing collision probability anddrift

ZizungYoon1· YeerangLim1· SebastianGrau1· WalterFrese1· ManuelA.Garcia1

Received: 20 September 2018 / Revised: 26 February 2020 / Accepted: 10 March 2020 / Published online: 19 March 2020

Abstract

The compact form factor of nanosatellites or even smaller satellites makes them predestined for distributed systems such as

formations, constellations or large swarms. However, when it comes to orbit insertion of multiple satellites, these ride share

payloads have constrains in the deployment parameters such as sequence, direction, velocity and time interval. Especially

for formation flight missions without propulsion, where the satellites should minimize their relative distance drift either

passively or by atmospheric drag control, the initial ejection parameters must find a proper trade-off between collision

probability and relative drift. Hence, this article covers short-term (first orbit) collision analysis, long-term (30days) drift

analysis and atmospheric drag control strategy for long-term distance control of multiple satellites. The collision analysis

considers various orbit deployment parameters such as insertion direction and tolerance, orbital elements of insertion and

time span. To cover the parameter space, a Monte Carlo simulation was conducted to identify the impact of these param-

eters. It showed that for collision probability the major factor is the time span between two ejections and the precision of

the deployment vector. For long-term drift analysis, orbit perturbation such as atmosphere and J2 terms are considered. The

result showed that for drift minimizing, minimizing the along-track variation is more substantial than reducing the time span

between ejections. Additionally, a drag control strategy to reduce the relative drift of the satellites is described. The results

have been applied on the S-NET mission, which consists of four nanosatellites with the task to keep their relative distance

within 400 km to perform intersatellite communication experiments. The flight results for orbital drift show equal or better

performance (0.1–0.7km/day) compared to the worst-case simulation scenario, implying that orbit perturbation was chosen

correctly and all orbit injection tolerances were within specified range. The drag control maneuver showed good matching

to the flight results as well with a deviation for the maneuver time of approximately 10%.

Keywords Orbit deployment· Formation flight· Swarm· Collision risk· Orbit drift· Drag control

List of symbols

drel

Relative distance between two satellites

Drag coefficient

Critical distance

e Eccentricity

i Inclination

J2 Zonal harmonic term

Number of iterations

Number of collisions

Quantity of random deployment vectors

nED

Quantity of equal distributed deployment

vectors

nsat

Number of satellites

Collision probability

qref

Dynamic pressure

𝜌

Atmospheric density

a Semi-major axis

𝜃

Angular distance in along-track direction

𝛥tSS

Time span between two consecutively deployed

satellites

𝐯𝐃

Nominal deployment velocity vector

𝐯𝐑

Random deployment velocity vector

* Zizung Yoon

zizung.yoon@tu-berlin.de

Yeerang Lim

y[email protected]

Sebastian Grau

[email protected]

Walter Frese

walter.frese@tu-berlin.de

Manuel A. Garcia

[email protected]

1 Chair ofSpace Technology, Technische Universität Berlin,

Berlin, Germany

398 Z.Yoon et al.

1 3



𝐯𝐫𝐞𝐥

Unit vector of relative velocity

𝕍RD

Set of random deployment velocity vectors

𝕍ED

Set of equally distributed deployment velocity

vectors

Abbreviations

BC Ballistic coefficient

CW Clohessy–Wiltshire

DL Downlink

DLR Deutsches Zentrum für Luft- und Raumfahrt

EAGF Earth’s aspherical gravity field

ECI Earth centered inertial frame

GMAT General Mission Analysis Tool

ILRS International Laser Ranging Service

IMU Inertial measurement unit

ISL Intersatellite link

ISS International Space Station

LEO Low earth orbit

LVLH Local vertical local horizontal frame

MEMS Microelectromechanical systems

PDF Probability density function

RAAN Right ascension of ascending node

SNET S-band network of distributed nanosatellites

SSO Sun-synchronous orbit

TLE Two-line element set format

UL Uplink

USSA76 U. S. Standard Atmosphere 1976

1 Introduction

The number of launches of small satellites, especially Cube-

Sats and nanosatellites, has been increased exponentially

over the last years. The increasing number of formation and

swarm mission such as the Dove, Flock [1] or Lemur [2]

constellation were heavily contributing to this trend. Still

these small satellites are launched as secondary payload, and

thus have constrains in orbital parameters and deployment

conditions dictated either by the primary payload or launch

vehicle upper stage.

Especially for satellite swarms/formation with strict

requirements on relative distance (due to, e.g., communi-

cation, distributed measurement or proximity operation), a

proper deployment strategy is crucial to minimize collision

risk and set optimal (e.g., propellant saving) drift conditions.

Since propulsion is not accommodated in most CubeSat

missions, passively minimizing the relative drift by proper

initial deployment could reduce system complexity of the

satellites.

Hence, this article covers short-term (first orbit) colli-

sion analysis based on Monte Carlo simulation, long-term

(30days) drift propagation incorporating J2 and atmospheric

density, and atmospheric drag control strategy for long-term

distance control of multiple satellites. The analysis has been

applied on the S-NET mission, which consists of four nano-

satellites with the task to keep their relative distance within

400 km to perform (ISL) experiments [3]. The mission was

launched in February 2018 to an 580km (SSO).

2 Motivation

2.1 Considerations forpiggy‑back swarm

deployment

In general, the orbit insertion condition of share-ride mis-

sions is not as flexible as the primary payloads and is

constrained by the upper stage configuration and ejection

mechanism. Hence, for optimal deployment of multiple

spacecrafts, the following constraints must be considered:

• All piggy-back satellites of a launch are normally

deployed in a single orbit (e.g., Dove satellites). Re-

ignition capability of upper stage may allow for several

orbits.

• The attitude determination and control accuracy of upper

stage have to be considered when defining the ejection

direction. Normally attitude determination of upper stage

is performed by the integration of an IMU; thus, the atti-

tude information depends on angular random walk and

time since lift-off. Typical values of attitude determina-

tion accuracy are few degrees in orbit.

• CubeSats and nanosatellites are mostly separated from

upper stage by a dispenser/container using a spring

mechanism. Compared to the classical separation ring

with pyro actuator, spring mechanism is beneficial in

controlling the relative ejection velocity and direction.

• For precise deployment, even the mechanical tolerance

of the dispenser mounted on the upper stage must be

considered.

Primarily, the initial deployment condition is set to minimize

collision risk with the upper stage or other spacecrafts of the

deployment sequence. Depending on mission requirements,

the secondary requirement would be to minimize the relative

drift of satellites (e.g., distributed nanosatellites of mission

(S-NET)) or maximize to achieve broad distribution in short

time (e.g., 3U Dove constellation by Planet Inc.).

2.2 Existing deployment strategies forformations/

swarms

Some deployment strategies have already been evaluated to

reach the objectives of each mission. For example, Kılıç [4]

analyzes the deployment strategy of a cluster of nanosatel-

lites. Also from the space debris mitigation aspect, proper

399

Orbit deployment anddrag control strategy forformation flight whileminimizing collision…

1 3

deployment is required to avoid fragmentation events [5].

A parametric study was utilized to test various deployment

parameters and thus analyze the dispersion of the 50 Cube-

Sats of the QB50 mission. Then, it is provided an insight

about the collision risk between the CubeSats, as well as, the

effect of the deployment strategy on the lifetime distribution

of CubeSats.

For swarm deployment, the first couple of orbits right

after the orbit injection are the most critical in terms of col-

lision risk. Even if satellites are equipped with propulsion,

especially for secondary payloads, a fixed waiting time for

initial activation (e.g., 30min for 6U CubeSats [6]) in order

to protect the main payload prohibits active orbit control

within that critical phase. Thus, the development of a right

deployment strategy is crucial to minimize collision between

the swarm members.

Dove satellites, a reference in this sector, have been

deployed in different ways. For example, 8 Flock 2E and

Flock 2E’ satellites were deployed from the (ISS) in May

2016. They were deployed in groups of two over four deploy-

ment windows by a spring deployer supplied by NanoRacks.

The system was attached to the Multi-Purpose Experiment

Platform of the Japanese Experiment Module. This is grap-

pled and set in deployment position and orientation by the

Japanese Remote Manipulator System [7].

3 Methodology

This section presents the propagator used and explains the

perturbation sources taken into account for the development

of this project. Afterward, the different methodologies con-

sidered for the collision and drift analyses are explained.

3.1 Orbital propagation

There are several classical special perturbation methods for

orbit propagation such as Cowell, Encke and DROMO. In

Cowell’s method, the equations are formulated in rectangu-

lar coordinates and integrated numerically. Encke’s method

approximates the orbit as an conic section. Since differences

grow with time, it becomes necessary to derive a new oscu-

lating orbit. This process is called rectification of the orbit.

The DROMO method is especially appropriated to propagate

complex orbits, e.g., comets and asteroids [8]. For the typi-

cal cases of this study (LEO satellites perturbed by atmos-

pheric drag and J2), Cowell’s method offers optimal results

in terms of time consumption and precision.

3.1.1 Cowell propagator

Cowell’s formulation for orbit propagation was initially pro-

posed by P. H. Cowell (Cambridge, UK). It is the most direct

approximation of all the perturbation’s methods available

nowadays, since the (Newtonian) perturbation forces are sim-

ply summed. The equation of motion in ECI system is given

in Eq.1:

where

and

are the accelerations caused by pertur-

bations. This work considers atmospheric drag which

requires velocity information. Therefore, Cowell’s equation

is reduced to a first-order system (Eq.2), which allows for a

broader class of integration methods. The used Cowell prop-

agator integrates the motion equations with a

Runge–Kutta–Fehlberg 7(8) numeric solver with variable

step size, where perturbations are included:

3.1.2 Orbital perturbation

For small satellites in LEO altitudes 400–800km, J2 per-

turbation of the Earth’s oblateness is the dominant force,

followed by lunar gravity and solar gravity. The atmospheric

density strongly depends on the solar activity and can vary

up to one order of magnitude [9]. The orbital perturbations

implemented and applied for obtaining the results of this

publication are J2 taken from (EAGF) and atmospheric drag

using USSA76 as density model.

Non-spherical Earth Harmonics

The Cowell simulation considered the zonal harmonic

term J2. The satellite’s acceleration due to the Earth’s non-

spheric gravity field EAGF can be expressed as a func-

tion derived from the gradient of the geopotential function

(Eq.3):

Atmospheric Drag

The acceleration experienced by the satellite due to

atmospheric drag is obtained by Eq.4:

(1)



=−𝜇

r3+apx(t,x,y,z,x,y,z)



=−𝜇y

r3+apy(t,x,y,z,x,y,z)



=−𝜇z

3+apz(t,x,y,z,x,y,z)

(2)

𝐗

[

𝐫

𝐯

]



𝐗

=[𝐯

−𝜇𝐫

3+𝐚p]

(3)

𝐚grav (𝐫,t)=∇

𝛷

(𝐫,t).

(4)

𝐚

𝐝=−

2𝜌v2

rel



𝐯𝐫𝐞𝐥

=−

𝜌v2

relBC−1

𝐯𝐫𝐞𝐥

400 Z.Yoon et al.

1 3

where BC is the ballistic coefficient (Eq.5), and



𝐯𝐫𝐞𝐥

the unit vector of velocity of the satellite relative to the

atmosphere:

In Eqs.4,5 and 6 ,

is the drag coefficient of the satellite,

A is its cross-sectional area normal to

𝐯𝐫𝐞𝐥

𝜌

is the local

atmospheric density,

is the spacecraft mass,

𝜔E

is the

angular velocity vector of the central body, and

𝐯w

stands

for the local wind velocity.

𝜌

is obtained from USSA76.

𝐯w

depends on the altitude, as empirical studies from satellites

have shown [10]. For altitudes between 200 and 300km,

𝐯w

might have up to 30% of

𝜔E

, but for altitudes above 500km,

it becomes reduced to a range where it can be neglected [10].

Additionally

𝐯w

, which is dominantly eastward wind, mainly

changes the i for high inclination orbits. Inclination change

is not primary relevant for the drag analysis.

3.2 Collision analysis

As Florijn points out in [11], the reduction of the collision

probability can be achieved through orbit design and deter-

mined by the satellite deployment method. The method he

utilized to calculate the collision probability is based on the

propagation of the nominal position of the satellites with an

error ellipsoid that describes the position and velocities by

a three-dimensional (PDF). It is concluded that the PDF-

based method is more suitable when a large number of sat-

ellites are going to be evaluated and the Monte Carlo-based

method to calculate collision probability is more suitable for

a few number of satellites due to its intrinsic computational

intensity.

Patera [12] calculates the probability of collision using a

nonlinear approach. The proposed method is applicable to

low-velocity space vehicle encounters that involve nonlin-

ear relative motion. According to him, “the total probability

of collision over a specified time is obtained by integrating

the collision rate over the appropriate time interval.” The

accuracy of this nonlinear model was verified with the com-

parison of a Monte Carlo simulation of 6000 runs, where

the collision probability was only 2% greater than the result

obtained by the Monte Carlo approach.

A Monte Carlo method is used to compute the collision

of various deployment scenarios. In a Monte Carlo statistical

method, the variable of interest is computed from a random

input variable. In this study,

𝕍RD ∋𝐯𝐑

represents the random

set of deployment velocity vectors of each satellite w.r.t the

upper stage. The random variation from the ideal (

𝐯𝐃

) is

derived from the sum of the following uncertainties:

(5)

⋅

(6)



𝐯𝐫𝐞𝐥 =𝐯−

𝝎

E×𝐫−𝐯w.

•

edis

: variation of ejection velocity vector due to the

ejection mechanism (generally spring or linear actua-

tor) and mechanical tolerance of dispenser itself

•

emech

: tolerance of the mechanical alignment between

dispenser and upper stage

•

eus

: error in attitude control of upper stage

Assuming normal distribution, the random deployment

vector

𝐯𝐑

is given by:

with

1𝜎

variation

Therefore, for every satellite,

|𝕍RD|=nR

. The permutation

of each satellite pair

i,j=1…n

sat)

results in

number of iterations that is required to obtain an accept-

able amount of values for the minimum (

drel

), which is the

parameter of interest. The bigger the

, the more accurate

is the result of the Monte Carlo approach. To find a trade-off

between computing time and accuracy, the following equa-

tion is used to obtain the order for

(ni)

[13]:

𝜀

is the percentage error of the mean. s is the calculated

standard deviation,



is the calculated mean, and

z𝛼∕2

is the

z-score associated with

𝛼∕2

where

𝛼

is the significance level.

Finally, collision is counted if

drel

≤

within orbit

propagation. Therefore, the collision probability is sim-

ply obtained by

3.3 Drift analysis

The study of drift can therefore estimate the lifetime of a

mission, its different phases or the time between maneuvers

for missions with active orbit control capabilities. If the drift

can be estimated and minimized by controlling the initial

deployment condition, even formation flight and proxim-

ity operations without propulsion could be feasible. Thus,

herein the drift analysis investigated the effect of initial orbit

insertion condition and orbit perturbation. For the direction

of deployment, the three main directions to analyze corre-

spond to the three axis of (LVLH) reference system:

(7)

𝐯𝐑

∼N(𝐯

𝐃

,𝜎

(8)

𝜎2=edis +emech +eus.

(9)

⋅

(

−1

)

∕

(10)

(

z𝛼∕2

100s

𝜀



Y)2

(11)

401

Orbit deployment anddrag control strategy forformation flight whileminimizing collision…

1 3

• Zenith deployment impulses in radial direction change

the shape (e) and in-plane orientation (

𝜔

and

𝜈

) of the

satellite, but not the size or plane of the orbit.

• Along-track deployment along-track impulses produce

a change in the shape and size of the orbit (a, i and e),

and also in the orientation in-plane (

𝜔

and

𝜈

), but do not

modify the plane.

• Cross-track deployment finally, impulses in cross-track

direction will affect mainly the plane of the orbit, chang-

ing i and

𝛺

To analyze the long-term drift behavior depending on initial

condition of orbit insertion and orbit perturbation, a Monte

Carlo method was applied. Similar to the collision analysis

method, the (

nED

) of the

𝐕

𝐃

,𝐄𝐃

is selected. The set of

𝐕

𝐃

,𝐄𝐃

is distributed around the nominal deployment velocity

𝐯𝐃

Each vector is propagated with Cowell propagator (Eq.1),

which returns

nED

propagated position sets for each satel-

lite. These sets are combined to obtain the relative distances

among the satellites. There are

nED ⋅(

−1

)

∕

combina-

tions for each pair of satellites

i,j=1…n

sat)

and each

drift between two satellites after a given amount of time.

3.4 Differential drag control

Differential drag is caused when the ballistic coefficients

(Eq.5) of the spacecraft in a formation are not equal. The

magnitude of differential drag depends on the difference in

ballistic coefficients and also the altitude of the spacecraft

formation. Differential drag is a favorable method to com-

pensate orbital energy difference, especially for small satel-

lites who do not carry active propulsion. Mean motion of

satellites could be matched by decreasing larger semi-major

axes values using higher drag.

After S-NET satellites were launched in February, con-

sistent secular drifts are observed in relative distances

between satellites. As anticipated from Clohessy–Wiltshire

equations [14], these drifts result from velocity component

in along-track direction. Expressed in Keplerian orbit ele-

ments, those linear increases are caused by discrepancies

in satellite’s orbital energies (i.e., semi-major axis dispari-

ties). Since mean motion of the satellites is a function of

semi-major axis only, those differences affect the satellite’s

angular velocity.

3.4.1 Differential drag controller

Since differential drag affects the in-plane dynamics only, six

relative states in the LVLH frame could be reduced to a single

angular distance

𝜃

and its derivative



𝜃

. The angle

𝜃

is defined

w.r.t. the reference satellite as shown in Fig.1. As shown,



𝜃

is equivalent to the satellite’s mean motion for circular orbits.

The reference satellite is the one which has the smallest

semi-major axis, because drag could only generate negative

acceleration and reduces the orbital energy.

The main idea of differential drag control is to accelerate

each satellite’s



𝜃

until it becomes the same as the reference

value, by increasing cross-sectional area and adjusting



𝜃



𝜃

zero when the satellite faces the same amount of drag as the

reference and has a nonzero value when the satellite is in high-

drag mode [1, 15]. To derive



𝜃

caused by differential drag,

induced

𝛥V

should be examined first. Since the cross-sectional

area and corresponding ballistic coefficient in high-drag mode

BCH

are different from the low-drag value

BCL

𝛥V

during

the discretized time period

𝛥t

could be derived as in Eq.12.

To simplify the formulation, dynamic pressure

qref

is set as an

average value from the reference satellite:

Considering the secular along-track term in CW equa-

tions only, mean relative velocity

y

could be approximately

derived as a function of

𝛥V

For the special case of circular orbits, angular acceleration



𝜃

could be converted from relative velocity

y

by using the

reference satellite’s mean semi-major axis at epoch

aref

Required time interval to maintain high-drag mode can be

obtained by dividing initial angular velocity offset with

angular acceleration above:

(12)

𝛥

2𝜌v2

rel(BCL−BCH)𝛥t

ref

(BC

−BC

)𝛥t

(13)



≈−3y

=−3

△

V=3q

ref

(BC

−BC

)𝛥t

(14)



𝜃

=y

a

ref

𝛥t=

ref

a

ref

(BCH−BCL)

(15)

𝛥

t=−

(



𝜃−



𝜃ref)



𝜃

Fig. 1 Relative state

𝜃

in along-track direction

402 Z.Yoon et al.

1 3

4 Analysis setup

4.1 Mission S‑NET overview

The mission goal of S-NET is to demonstrate intersatellite

communication with distributed nanosatellites. The mis-

sion consists of four satellites, to test multi-hop communi-

cation with different protocols and routing algorithm. Four

satellites allow for redundancy on satellite level, since the

ISL can also be performed with some restrictions on three

satellites.

To reduce system complexity of the space segment, a

propulsion system was omitted, and therefore, the rela-

tive formation will be controlled by the initial separa-

tion parameters. The goal is to keep the satellites within

a range of

400 km

to enable a stable ISL for at least four

months. Therefore, the separation direction, angle, veloc-

ity and sequence will be the primary starting values for

the formation drift. A summary of mission parameters is

provided in Table1.

The four S-NET satellites were accommodated in dis-

pensers (Fig.2) and inserted into a

580 km

SSO via Soyuz/

Fregat. The re-ignition capability of Fregat upper stage

allowed for a circular orbit even for secondary payloads.

The nominal orbital elements for the upper stage at time of

ejection are specified in Table2. These values were used for

both collision and drift analysis. Only for true anomaly of

collision analysis, a different value of

30◦

was used. Besides,

for both analyses, the considered perturbations for the propa-

gation with Cowell are listed in Table3.

The considered distance from one dispenser to the next

dispenser’s CoG is

0.5 m

, and the time interval between

consecutive separations

𝛥tSS

was set to

10s

according to the

simulation results in this article.

Finally, the characteristics of the nominal deployment

vector

𝐯𝐃

that generate both the random set of deploy-

ment vectors

VD,R

for the collision analysis and the equal

Table 1 Overview of mission

S-NET [16]Parameter Value

No. of satellites 4

Orbit height 580 km SSO

Launch date 2018-02-01

Mass 8.8 kg per SC

Volume 240

240

mm3

Communication UHF (satellite bus)

S-band (payload)

Attitude Control Three axis with MEMS sensor arrays,

reaction wheels, 3

magnetorquers

Payload SLink TRX for ISL, UL & DL Laser reflector for high precision

position measurement

Ground station UHF: Berlin, Svalbard, Backnang S-band: Berlin, DLR Neustrelitz

Fig. 2 Dispensers for S-NET mounted on Fregat upper stage. See

Fig.3 for dispenser orientation and naming conventions

Table 2 Initial nominal classical orbital elements of the upper stage

(

COE

) for the drift analysis

Element Value

Semi-major axis (

aus

)

6971 ±8

Eccentricity (

eus

)

0.0001 ±0.0016◦

Inclination (

ius

)

97.75 ±0.12◦

Argument of perigee (

𝜔us

)

205.02 ±150◦

Right ascension of ascending node (

𝛺us

)

226.25◦

True anomaly (

𝜈us,i

)

120◦

403

Orbit deployment anddrag control strategy forformation flight whileminimizing collision…

1 3

distributed set of deployment vectors for the drift analysis

VD,ED

are summarized in Table4.

4.2 Collision analysis

From the nominal deployment parameters, the mean mini-

mum relative distance between two consecutive satellites

is obtained as 13.258 m with a standard deviation of 3.371

𝛼

∕2

= 2.58 (99% level of confidence) and

𝜀

= 1%. Thus,

according to Eq.10, at least 4305 iterations are required to

achieve 99% confidence level that the calculated value of the

minimum

drel

between two satellites is within 1% of the true

minimum distance’s value. Finally, a sample size of 250.000

provides an accurate level for the real mean and

is set to

500. The critical distance

is set to

0.30 m

. To compare the

effects of perturbations, two propagation methods were used

to run the simulation: One refers to the Keplerian propaga-

tion without perturbations and the other uses the Cowell

propagator to integrate the orbit, inclusive of J2 and atmos-

pheric drag perturbations. Also, different deployment direc-

tions and separation sequences were considered to figure out

the optimal deployment strategy.

4.2.1 Deployment direction

The direction of deployment is defined by the attitude con-

trol of the upper stage. The Euler rotation sequence from

XYZ (along-track/cross-track/nadir) to the deployment

direction assumed in this paper is pitch, followed by yaw

and then by roll. The reference direction is positive along-

track where the pitch is

0◦

and yaw

0◦

(P0Y0). Hence, to

point a deployment toward the zenith direction, the rotation

of the upper stage is

90◦

in pitch and

0◦

in yaw (P90Y0),

as per schematic in Fig.3. The cross-track direction is

(P0Y90), respectively. The nominal deployment is shown

in the schematic in Fig.3, which is the baseline to compare

other deployment parameters. The nominal ejection veloc-

ity vector

𝐯𝐃

has a magnitude of

1.5 ms−1

and point toward

zenith (P90Y0).

4.2.2 Deployment order

The deployment order changes the initial distances of sat-

ellites. Linear A configuration (P3–P1–P4–P2) means that

S1 is deployed from position P3, then S2 from P1, then S3

from P4 and finally S4 from P2. Thus, the initial distance

on the upper stage is maximized. Linear B configuration

(P1–P2–P3–P4) implies that the directly neighboring satel-

lites are ejected beginning from positions P1, P2, P3 and

finally P4. A summary is given in Table5.

4.3 Drift analysis

For the drift analysis, all propagations (upper stage and

deployed satellites) are carried out with the same set of

perturbation forces that includes J2 and are propagated

with Cowell. Atmospheric drag was not considered since

the effect is negligible for the altitude of

580 km

. The sim-

ulations have been carried out for the three main direction

zenith, along-track and cross-track of the nominal deploy-

ment, to achieve results for

30 d

from insertion. Since the

three main directions provide the worst and best cases for

Table 3 Considered perturbations for the S-NET collision and drift

analysis

Perturbation Model

Earth’s aspherical gravity field

(EAGF)

J2 spherical harmonic

Atmospheric drag (AD) U.S. Standard Atmosphere 76

Table 4 Deployment velocity vector and tolerance for S-NET

Element Value

Magnitude of

𝐯𝐃

1.5 ms−1±3%(1

𝜎

)

Direction of

𝐯𝐃

90◦±2◦(1𝜎)

0.30 m

Table 5 Setup for the collision simulation

Propagator Duration Perturb.

Direction Order

Kepler

6400 s

no 500 P90Y0 (zenith) A

P90Y90 (zenith) B

Cowell

(Runge

Kutta 7)

6400 s

J2 USSA76 500 P90Y0 (zenith) A

P90Y90 (zenith) B

Fig. 3 Schematic of the upper stage’s attitude for deployment toward

zenith direction. Nominal deployment conditions:

COEA

, Linear A,

P90Y0

404 Z.Yoon et al.

1 3

drift, parametrical exploration of those combinations was

not done in this work. Especially for a nearly circular orbit,

the result can be compared with the Clohessy–Wiltshire

equations for relative motion.

To estimate the behavior of the four satellites after

deployment with reduced computation load, a set of

equally distributed worst-case deployment velocity vec-

tors has been defined, based on the tolerance provided in

Table4.

The set of deployment velocity vectors

𝕍ED

has been

reduced by defining the maximum deviations from the

nominal values for each case. Taking the full set of ran-

dom vectors

𝕍RD

would exceed acceptable computation

time, especially for 35days of drift propagation. Thus,

𝕍ED

is defined as:

and illustrated in Fig.4.

Thus, it results in 17 vectors and covers the most devi-

ated cases from the nominal vector

𝐯𝐃

as:

(16)

𝕍ED =

𝐯

𝐃∪

𝐯

�

𝐄𝐃 ∪

𝐯

��

𝐄𝐃

• The nominal deployment vector

𝐯𝐃

, refer to Table4 for

values

•

𝐯′

𝐄𝐃

: magnitude:

1.5 ms−1+3%(1𝜎)

, angle offset

𝜃=2◦(1𝜎)

w.r.t

𝐯𝐃

, azimuth variation:

0, 45 …360◦

•

𝐯′′

𝐄𝐃

: magnitude:

1.5 ms−1−3%(1𝜎)

, angle offset

𝜃=2◦(1𝜎)

w.r.t

𝐯𝐃

, azimuth variation:

0, 45 …360◦

A summary of drift simulation setup is given in Table6.

4.4 Drag control maneuver planning

For S-NET, S4 is at the very front and S2 has the slowest

angular velocity. Obviously, S2 will be the one with high-

drag mode, which has the maximum cross-sectional area.

The attitude of the following satellite is assumed to be

tilted

45◦

, and the leading satellite is assumed to continu-

ously tumble with respect to the velocity direction, which

generates 17% of area difference.

Maximum cross-sectional area for S-NET is

√3

times

larger than the square cross section.

Although the area difference is small between two satel-

lites, it is still sufficient for S-NET mission because secular

drifts are gradual compared to general cases.

To calculate required

𝛥t

for drag control maneuver to stop

relative drift, atmospheric density is set as a rough average

based on the ephemeris, which is

𝜌=2.5 ×10−14 kg∕m3

, and

mean dynamic pressure of the reference satellite, i.e., S4, is

derived as a function of reference velocity:

The drag coefficient depends on numerous factors such as

shape of the satellite, momentum exchange with molecular

particles, satellite surface temperature and chemical reac-

tions at the satellites surface. Since satellites do not incor-

porate methods of measuring it, the numerical value must

be estimated. Empirical values from orbital measurements

given by [17] indicate a good estimation of

Cd=2.3

for the

cube-shaped S-NET satellites at 580km altitude. The only

unknown parameters in the equation are semi-major axis

and the reference velocity. Since two semi-major axes are

almost identical to each other, semi-major axes value should

be precisely determined. For instance, Brouwer mean con-

version oscillates in

102

m level, and the semi-major axis

difference is in

100

m order. To get accurate semi-major axes,

(17)

ref =

𝜌×v2

rel

Fig. 4 Sketch of deployment vectors

𝕍ED

used for the drift analysis

Table 6 Setup for the drift

simulation Propagator Duration Perturb.

Direction Order

Cowell (Runge Kutta 7) 35 d J2 17 P90Y0 (zenith) A

P0Y90 (cross) B

P0Y0 (along-track)

405

Orbit deployment anddrag control strategy forformation flight whileminimizing collision…

1 3

(GMAT) propagation was conducted with 50x50 EGM96

gravitational model and the resulting positional differences

were averaged over time.

The a difference (values given in Table7) causes mean

motion difference of approximately

1.182 ×10−9s−1

. This

again results in

3.062 ×10−3

rad angular drift in along-track

after 30days. This is equivalent to

21 km

range drift, which

is similar to actual behavior from measurement data pre-

sented in Fig.15. Once all the variables are prepared, mean

dynamic pressure and the angular acceleration values could

be achieved as follows:

and thus,

And the angular velocity of each satellites is obtained by:

Required

𝛥t

estimation results to stop the drift are as given

in Table8.

5 Results formission S‑NET

5.1 Results forcollision analysis

5.1.1 Cowell propagation

The Cowell propagation results are plotted in Fig.5. The

histogram shows that

90%

of the collision probability occurs

within the first

10 min

after the deployment process starts.

The next potential collision window is after half orbital

(18)

ref =1

⋅𝜌⋅

(√

𝜇

aref

−𝜔E⋅aref ⋅cosi

=7.29 ⋅10−7kg

(19)



𝜃

ref



aref

(BCH−BCL)=1.06 ⋅10−15 1

(20)



𝜃

√

𝜇

period and then again after one full orbit, where the satellites

approach the initial state. After one orbit, the cumulative

probability reaches close to 1.

5.1.2 Components ofthedeployment vectors thatcause

collision

Here, the collision cases are examined in detail. Figure6

shows the distribution and mean values of azimuth, elevation

and magnitude of the deployment velocity vectors that cause

a collision for the satellite combination S3S4. The bottom

figure verifies the influence of deploying with two different

velocity magnitudes on the collision probability. From this

graph, it is observed that it is more probable for a collision

Table 7 Semi-major axis of S4 and S2, based on flight data

Sat4 a (m) Sat2 a (m)

𝛥

a (m)

6.296e+06 6.296e+06

−5.029

Table 8 Required

𝛥t

for drag control maneuver to stop relative drift

Sat1 Sat2 Sat3

Required (day) 7.29 12.84 6.00

Fig. 5 Collision histogram and cumulative distribution of nominal

deployment over a orbital period with J2 and atmospheric drag

Fig. 6 Mean value of the components of the deployment velocity vec-

tors that cause a collision, including J2 and atmospheric drag

406 Z.Yoon et al.

1 3

to occur when the deployment velocity’s magnitude of the

second deployed nanosatellite is greater than the first one.

5.1.3 Effect ofdeployment configuration

The linear A (P3–P1–P4–P2) configuration with an attitude

of P90Y0 was originally proposed under the assumption

that a greater distance between the slots of the satellite’s

dispenser would generate a lower collision probability. How-

ever, the linear B (P1–P2–P3–P4) configuration (P90Y0)

was tested in order to verify that assumption. From the

results provided by the total cumulative collision graph,

which is presented in graph Fig.7, it is observed that the

Linear A drives a

0.2048%,

while the

of Linear B

(P90Y0) is

0.1788%

. Therefore, the distance between slots

is not the only parameter affecting collision probability, it

is also influenced by the order in which they are deployed

and the attitude of the upper stage. For example, a simula-

tion with Linear A configuration and upper stage’s attitude

of P90Y90 provides a

0.1748%

, while Linear B with

P90Y90 drives a

0.1800%

, showing that a greater dis-

tance between the slots has more influence when the upper

stage’s dispenser is oriented perpendicular to the orbital

path.

5.1.4 Effect ofdeployment time span

Among the simulated deployment parameters, the span

between the satellites insertion

𝛥tSS

results in the great-

est variation of collision. Figure8 provides the collision

probabilities w.r.t

𝛥tSS

for two cases: without perturbation

(Keplerian orbit) and with perturbation (J2 and atmospheric

drag). The figure shows for

𝛥tAB =5

s a collision probability

difference of 0.23% (0.80% vs. 1.03%) only. According to

the results, the propagation method and the perturbation

effects are negligible. After

10 s

converges to less than

1%.

5.1.5 Effect ofdeployment direction

The effect of the deployment direction is analyzed by varying

the yaw and pitch angle depicted in Fig.9. The red colored

zone highlights the highest

and occurs in the cross-track

direction (P0Y90) with value of

0.364%

. In contrast, the low-

est

0.204%

occurs in zenith direction (P90Y0) which is

Fig. 7 Comparison between Linear A and Linear B deployment con-

figurations. Vertical red line marks one orbital period

Fig. 8 Effect on overall collision probability w.r.t time spans

𝛥tAB

after one orbit with perturbation and without perturbation

Fig. 9 Collision probability for six deployment directions includ-

ing J2 and atmospheric drag,

𝛥tSS

=10s. Cross-track (P0Y90)

 = 0.364%, along-track (P0Y0)

 = 0.29%, zenith (P90Y0)

=0.204%

407

Orbit deployment anddrag control strategy forformation flight whileminimizing collision…

1 3

colored in dark blue. The along-track (P0Y0) is intermediate

with

=0.29%. Despite variations, the effect of deploy-

ment direction is rather small, since the majority of colli-

sions occur within the first 600s (see also Fig.7), where the

relative motion effect is not fully unfolded yet.

5.1.6 Effect oforbital elements oftheupper stage

To analyze the effect of the orbit characteristic on

, the

orbital elements e and a of insertion orbit were varied in

discrete steps. Figure10 plots the results for a range of

e=0…0.015

and

a=6871 …7071

km. It shows a varia-

tion of only 0.01% (min: 0.20%, max: 0.21%) for a nearly

circular orbit. The dashed white rectangle highlights the tol-

erance window of the upper stage’s orbit insertion tolerance

specification (values also given in Table2). In this range, the

only changes

±0.002%

, indicating that a and e introduce

negligible variability for nominal orbit insertion.

5.2 Drift behavior

Simulations were run for the four S-NET satellites with time

span between satellites set to

𝛥t

=10s, which resulted in suffi-

ciently small

=0.204%. For each direction, a set of point-

ing errors was added to the deployment vector, as discussed

in Section4.1. The results are presented in three subsec-

tions, depending upon the deployment direction employed.

Since the purpose was to identify the extrema and along-

track results in maximum drift (worst case), a parametrical

exploration was not part of this paper.

5.2.1 Along‑track deployment

The drift produced after along-track deployment is the larg-

est among the three directions considered. Figure11 shows

the maximum drift obtained combining the sets of initial

vectors for each combination of satellite pairs. The drift in

this case is very large, reaching almost

780 km

for this maxi-

mum-drift case. The maximum drift corresponds in all cases

to combinations with maximum difference in deployment

speed (one satellite with an initial

1.45 ms−1

and the other

one with

1.55 ms−1

), which provoke the maximum variation

in a and therefore in the orbital period.

5.2.2 Cross‑track deployment

The results for cross-track direction are obviously very dif-

ferent compared with along-track. Cross-track deployment

causes mainly a change in an orbit’s plane (i) and not size

(a). In Fig.12, the maximum drift for each combination is

shown, reaching about

21 km

drift after 30days. Again, the

results are very similar for the six pair combinations.

Fig. 10 Collision probability w.r.t variation of e and a shows a varia-

tion of 0.01%. White dotted rectangle represents orbit insertion nomi-

nal tolerance window

Fig. 11 Maximum drift for each combination for along-track deploy-

ment

Fig. 12 Maximum drift for each combination for cross-track deploy-

ment

408 Z.Yoon et al.

1 3

5.2.3 Zenith deployment

Observing Fig.13, the maximum drift obtained for zenith

deployment is of the same order as for cross-track, namely

less than

22 km

for the worst case. Obviously, the closer the

deployment order, the smaller the drift. The drift is caused

by the along-track component of

𝕍ED

. It is thus important to

minimize this component in the launch, to comply with the

requirements of the mission.

Finally, Fig.14 shows the maximum drift after 30days

of propagation by varying

𝛥tSS

. It presents an almost linear

variation of drift versus time span; thus, to minimize the

drift, it is necessary to deploy the satellites with minimum

time span as possible. Thus, a trade-off between

(Fig.8)

and drift (Fig.14) must be found by defining proper

𝛥tSS

5.2.4 Relative distance inorbit

The relative drift among the satellites was measured in orbit

using the signal run time of the S-band radio during ISL

sessions. The measurement was correlated with TLE data

and ILRS data to increase the accuracy and is illustrated in

Fig.15. Interestingly, the maximum value of measured

drel

occurs between satellite S4 and S2 and coincides with the

simulation result for zenith deployment, given in Fig.13.

All other satellite distances are smaller than the simulation

case. Overall, it implies that the tolerance assumptions made

for

𝐯𝐃

(Table4) were realistic and all launch requirements

have been fulfilled.

Fig. 13 Maximum drift for each combination for zenith deployment

Fig. 14 Effect of different deployment time spans on total drift

obtained after 30 days of propagation (J2 and drag), for the worst

between satellites 1 and 4 in zenith deployment

Fig. 15 Simulation (worst case) vs. orbit data of relative drift shows

that drift value is less or equal (S4-S2) compared to simulation

results. This indicates that simulation assumptions regarding orbital

perturbation and ejection vector tolerances were correct and that the

real initial condition was within the specified tolerances

Fig. 16 Flight results of relative distance control. Nonlinear drift (red

box) is a result of differential drag control started at day 190. Since

SNET-A and D were controlled, the not involved pair B-C remains

linear

409

Orbit deployment anddrag control strategy forformation flight whileminimizing collision…

1 3

5.3 Drag control results

The graph in Fig.16 shows the real relative distance plot.

The experiment has been started at day 190. The most

ahead flying S-NET A was put into nadir pointing mode,

to allow other satellites to catch up. SNET-D is the slowest

satellite, and it is following SNET-A and SNET-C which

are randomly tumbling. The plots are from two-line ele-

ments analysis. The acceleration is

−3.8 m∕d2

between

SNET-D and SNET-C and

−

3.0 m

∕

between SNET-D

and SNET-A from the plot fitting. Average acceleration is

about

−

3.4 m

∕

. Due to some operational constraints, the

accumulated drag control was performed for 247.8h during

151days. The expected acceleration for this case can be

averaged as:

The calculated and scaled acceleration is consistent with the

in-orbit result with an error of about 10%. This error could

be cased by

estimation error, variable atmospheric den-

sity, irregular tumbling and attitude control error.

Non-constant atmospheric density during propagation is

surmised as a cause for vestigial secular behavior.

6 Conclusion

An analysis on various deployment conditions to develop an

orbit deployment and drag control strategy to minimize the

collision and drift rate for distributed space systems (e.g.,

swarms or formations) is presented. This analysis includes

the impact of deployment conditions such as orbital ele-

ments, time span, deployment direction and order for mul-

tiple satellites on collision probability and relative drift.

Due to the large number of relevant parameters to define

the insertion condition, a Monte Carlo simulation was con-

ducted. Only perturbations due to J2 and atmospheric drag

(

𝜌

from USSA76) were considered. The parameter space was

limited to a SSO LEO swarm mission assuming a near cir-

cular orbit, and as a baseline four satellites where simulated.

The upper stage was assumed to provide flexible pointing

and deployment configuration. At the end, the analysis was

applied for the distributed nanosatellite mission S-NET. The

following aspects can be concluded from the analysis:

• For

analysis, the variation of e and a is mostly negli-

gible for nearly circular orbits (Fig.10).

• For

analysis, the relevant time window to consider is

in the order of few orbits after insertion (Fig.7). In fact,

approximately 90% of all collisions events occur within

the first 600s (Figs.5, 9). After one orbit, the probability

(21)

y

=

𝜃⋅aref =−55.3

⋅

247.8h

24h⋅151

=−3.78

of collision approaches zero due. Thus, the deployment

direction does not have significant impact for

since the

effect of relative motion is not fully unfolded for the first

600s.

• For

analysis, the major factor is the time span

𝛥tSS

and

the precision of the deployment vector. Even the effect

of drag and J2 perturbation can be neglected due to the

short time window of relevance (Fig.8).

• Long-term (few months) drift analysis considering J2 and

atmospheric drag showed very good matching to flight

results. The effect of sun and moon gravitation and solar

pressure were not analyzed in this work.

• Cross-track and zenith deployment do not show significant

differences in drift rate. However, cross-track deployment

can generate a drift in the orbital plane due to differential

(RAAN) rate, which could be more difficult to correct via

differential drag. RAAN control could be achieved via dif-

ferential lift, as introduced in [18].

• Along-track component does not significantly reduce

Pc,

but is major cause of long-term drift. So these three main

directions result in the extrema boundaries for drift. So to

minimize drift, minimizing the along-track component is

the key, more than reducing the time span

𝛥tAB

(Figs.12

and 14 ).

• The value of

to obtain BC has been estimated from [17],

since direct measurement is not feasible. Even though drag

control maneuver in order to stop the relative drift rate

showed good matching to the flight results with a deviation

for the maneuver time of approximately 11% (Section5.3),

this was good enough to apply the drag control maneuver

successfully in orbit.

• Table9 contains recommended deployment configuration,

direction, time span and resulting collision probability and

maximum drift rate for the mission S-NET.

Acknowledgements Open Access funding provided by Projekt DEAL.

We are grateful to all project participants for their contributions, espe-

cially the colleagues from the German Aerospace Center (DLR) for

their technical advices and management support. Special thanks to the

company IQ wireless for their great contribution to the S-band com-

munication technology and Astro- und Feinwerktechnik Ltd., for their

high precision deployment system, without their work, it would not be

Table 9 Selected parameters for the S-NET mission

Parameter Value

Deployment configuration Linear B (P1–P2–P3–P4)

Deployment direction P90Y0 (Zenith)

Time span

𝛥tSS

10 s

Collision probability

0.1788%

Maximum drift (S1S4)

22 km∕35

days

410 Z.Yoon et al.

1 3

possible to minimize the relative position drift of the S-NET cluster.

Also we are very pleased to have been working with Exolaunch and

Roscosmos, thanks for a perfect launch and a precise orbit insertion.

The project S-NET is implemented on behalf of the Federal Ministry

of Economy and Technology (ger. Bundesministerium für Wirtschaft

und Technologie, BMWT) and the German Aerospace Center (DLR)

under grant no. 50YB1225.

Open Access This article is licensed under a Creative Commons Attri-

bution 4.0 International License, which permits use, sharing, adapta-

tion, distribution and reproduction in any medium or format, as long

as you give appropriate credit to the original author(s) and the source,

provide a link to the Creative Commons licence, and indicate if changes

were made. The images or other third party material in this article are

included in the article’s Creative Commons licence, unless indicated

otherwise in a credit line to the material. If material is not included in

the article’s Creative Commons licence and your intended use is not

permitted by statutory regulation or exceeds the permitted use, you will

need to obtain permission directly from the copyright holder. To view a

copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

References

1. Foster, C., Hallam, H., Mason, J.: Orbit determination and differ-

ential-drag control of planet labs cubesat constellations. In: AIAA

(ed.) AIAA Astrodynamics Specialist Conference (2015)

2. Hand, E.: CubeSats promise to fill weather data gap. Science

350(6266), 1302–1303 (2015). https ://doi.org/10.1126/scien

ce.350.6266.1302

3. Yoon, Z., Frese, W., Bukmaier, A., Briess, K.: System design of

an s-band network of distributed nanosatellites, CEAS Space J.

6(1). https ://doi.org/10.1007/s1256 7-013-0058-1

4. Kılıç, Ç, Scholz, T., Asma, C.: Deployment strategy study of

qb50 network of cubesats, IEEE—(2013). https ://doi.or g/10.1109/

RAST.2013.65813 49.http://ieeex plore .ieee.org/docum ent/65813

49/

5. Oltrogge, D.: Debris mitigation strategies, standards and best prac-

tices with sample application to CubeSats. In: Improving Space

Operations Workshop, San Antonio, TX (2013)

6. SLO, C.P.: 6u cubesat design specification, Tech. rep., California

Polytechnic State University (2016).http://www.cubes at.org/

7. Blau, P.: Flock 1/2 - planet labs earth observation satellites, Space-

flight 101.http://space fligh t101.com/flock /

8. Peláez, J., Hedo, J.M., de Andrés, P.R.: A special perturbation

method in orbital dynamics. Celest. Mech. Dyn. Astron. 97(2),

131–150 (2007)

9. Hallmann, W., Ley, W., Wittmann, K.: Handbook of Space Tech-

nology. Wiley, Hoboken (2008)

10. King-Hele, D.G., Walker, D.: Upper-atmosphere zonal winds from

satellite orbit analysis. Tech. rep, ROYAL AIRCRAFT ESTAB-

LISHMENT FARNBOROUGH (ENGLAND) (1982)

11. Florijn, D.: Collision analysis and mitigation for distributed space

systems, Master’s thesis, Delf Universita of Technology (2015)

12. Patera, R.P.: Satellite collision probability for nonlinear relative

motion. J. Guid. 68, 1–26 (2003). https ://doi.org/10.1016/j.actaa

stro.2016.04.031

13. Oberle, W.: Monte Carlo simulations: number of iterations and

accuracy. US Army Research Laboratory (2015)

14. Clohessy, W., Wiltshire, R.: Terminal guidance system for satellite

rendezvous. J. Aerosp. Sci. 27(9), 653–658 (1960). https ://doi.

org/10.2514/8.8704

15. Foster, C., Mason, J., Vittaldev, L., Vivek an, L., Beukelaers, V.,

Stepan, L., Zimmerman, R.: Differential drag control scheme for

large constellation of planet satellites and on-orbit results. In: 9th

International Workshop on Satellite Constellations and Formation

Flying (2017)

16. Yoon, Z., Frese, W., Briess, K.: Design and implementation of a

narrow-band intersatellite network with limited onboard resources

for IoT. Sensors 19(19), 4212 (2019)

17. Moe, K., Moe, M., Wallace, S.: Improved satellite drag coeffi-

cient calculations from orbital measurements of energy accom-

modation. J. Spacecr. Rockets 35(3), 266 (1998). https ://doi.

org/10.2514/2.3350

18. Traub, C., Herdrich, G.H., Fasoulas, S.: Influence of energy

accommodation on a robust spacecraft rendezvous maneuver

using differential aerodynamic forces. CEAS Space J. 12(1),

43–63 (2020). https ://doi.org/10.1007/s1256 7-019-00258 -8

Publisher’s Note Springer Nature remains neutral with regard to

jurisdictional claims in published maps and institutional affiliations.