Nonlinear dynamics near equilibrium points for a Solar Sail ` - - PowerPoint PPT Presentation

Nonlinear dynamics near equilibrium points for a Solar Sail

Ariadna Farr´ es ` Angel Jorba

ari@maia.ub.es angel@maia.ub.es

Universitat de Barcelona Departament de Matem` atica Aplicada i An` alisi

WSIMS – p.1/36

Contents

• Introduction to Solar Sails.
• Families of Equilibria.
• Periodic Motion around Equilibria.
• Reduction to the Centre Manifold.

WSIMS – p.2/36

What is a Solar Sail ?

• It is a proposed form of spacecraft propulsion that uses large membrane

mirrors.

• The impact of the photons emitted by the Sun onto the surface of the sail

and its further reflection produce momentum.

• Solar Sails open a new wide range of possible mission that are not

accessible for a traditional spacecraft.

WSIMS – p.3/36

Some Definitions

• The effectiveness of the sail is given by the dimensionless parameter β,

the lightness number.

• The sail orientation is given by the normal vector to the surface of the sail

(

n), parametrised by two angles, α and δ, where α ∈ [−π/2, π/2] and δ ∈ [−π/2, π/2].

α δ Sun-line

• n

x y z Ecliptic plane

WSIMS – p.4/36

Equations of Motion (RTBPS)

• We consider that the sail is perfectly reflecting. So the force due to the sail

is in the normal direction to the surface of the sail

n.

• Fsail = β ms

r2

ps

• rs,

n2 n.

• We consider the gravitational attraction of Sun and Earth: we use the

RTBP adding the radiation pressure to model the motion of the sail.

1 − µ µ

• FS
• FSail
• FE

Sun Earth

Y X t

WSIMS – p.5/36

Equations of Motion (RTBPS)

The equations of motion are:

¨ x = 2 ˙ y + x − (1 − µ)x − µ r3

ps

− µx + 1 − µ r3

pe

+ β 1 − µ r2

ps

• rs,

n2nx, ¨ y = −2 ˙ x + y − 1 − µ r3

ps

+ µ r3

pe

• y + β 1 − µ

r2

ps

• rs,

n2ny, ¨ z = − 1 − µ r3

ps

+ µ r3

pe

• z + β 1 − µ

r2

ps

• rs,

n2nz,

where,

nx = cos(φ(x, y) + α) cos(ψ(x, y, z) + δ), ny = sin(φ(x, y, z) + α) cos(ψ(x, y, z) + δ), nz = sin(ψ(x, y, z) + δ),

with φ(x, y) and ψ(x, y, z) defining the Sun - Sail direction in spherical coordinates (

rs = rps/rps ).

WSIMS – p.6/36

Equilibrium Points

• The RTBP has 5 equilibrium points (Li). For small β, these 5 points are

replaced by 5 continuous families of equilibria, parametrised by α and δ.

• For a fixed and small β, these families have two disconnected surfaces, S1

and S2. It can be seen that S1 is diffeomorphic to a sphere and S2 is diffeomorphic to a torus around the Sun.

• All these families can be computed numerically by means of a continuation

method.

WSIMS – p.7/36

Equilibrium Points

• 1
• 0.8
• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8 1

• 1.5
• 1
• 0.5

0.5 1 1.5 y x

• 0.04
• 0.03
• 0.02
• 0.01

0.01 0.02 0.03 0.04

• 1.02
• 1.01
• 1
• 0.99
• 0.98
• 0.97
• 0.96
• 0.95

y x

Equilibrium points in the {x, y}- plane

• 0.2
• 0.15
• 0.1
• 0.05

0.05 0.1 0.15 0.2

• 1.5
• 1
• 0.5

0.5 1 1.5 z x

• 0.015
• 0.01
• 0.005

0.005 0.01 0.015

• 1.015
• 1.01
• 1.005
• 1
• 0.995
• 0.99
• 0.985

z x

Equilibrium points in the {x, z}- plane

WSIMS – p.8/36

Some Interesting Missions

• Observations of the Sun provide information of the geomagnetic storms, as

in the Geostorm Warning Mission.

Sun Earth x y z

0.01 AU 0.02 AU L1 ACE

Sail CME

• Observations of the Earth’s poles, as in the Polar Observer.

Sun Earth x z

L1

N S Summer Solstice Sail Sun Earth x z

L1

N S Winter Solstice Sail WSIMS – p.9/36

• C. McInnes, “ Solar Sail: Technology, Dynamics and Mission Applications.”, Springer-Praxis,

1999.

• D. Lawrence and S. Piggott, “ Solar Sailing trajectory control for Sub-L1 stationkeeping”,

AIAA 2005-6173.

• J. Bookless and C. McInnes, “Control of Lagrange point orbits using Solar Sail propulsion.”,

56th Astronautical Conference 2005.

• A. Farr´

es and `

• A. Jorba, “Solar Sail surfing along familes of equilibrium points.”, Acta

Astronautica Volume 63, Issues 1-4, July-August 2008, Pages 249-257.

• A. Farr´

es and `

• A. Jorba, “A dynamical System Approach for the Station Keeping of a Solar

Sail.”, Journal of Astronautical Science. ( to apear in 2008 )

WSIMS – p.10/36

From now on ...

We fix α = 0 and β = 0.051689.

Earth Sun

L1 L2 L3 SL1 SL2 SL3

• Here, we have 3 families of equilibrium points on the {x, z} - plane

parametrised by the angle δ.

• The linear behaviour for all these equilibrium points is of the type

centre×centre×saddle.

• We want to study the families of periodic orbits that appear around these

equilibrium points for a fixed δ.

• For practical reasons we focus on the region around SL1.

WSIMS – p.11/36

Family of equilibrium points around SL1 for α = 0 and β = 0.051689

• 0.015
• 0.01
• 0.005

0.005 0.01 0.015

• 1.005
• 1
• 0.995
• 0.99
• 0.985
• 0.98

z x Earth L1 SL1

WSIMS – p.12/36

Motion around the equilibrium points

• As we have said, the linear behaviour around the fixed point is

centre×centre×saddle.

• So up to first order the solutions around the fixed point are:

φ(t) = A0[cos(ω1t + ψ1) v1 + sin(ω1t + ψ1) u1] + B0[cos(ω2t + ψ2) v2 + sin(ω2t + ψ2) u2] + C0eλt vλ + D0e−λt v−λ

Where,

• ±iω1 eigenvalues with

v1 ± i u1 as eigenvectors.

• ±iω2 eigenvalues with

v2 ± i u2 as eigenvectors.

• ±λ eigenvalues with

vλ, v−λ as eigenvectors.

WSIMS – p.13/36

Motion around the equilibrium points

• We take the linear approximation to compute an initial periodic orbit for

each family. We then use a continuation method to compute the rest of the family.

• Planar family: A0 = γ and B0 = D0 = E0 = 0.
• Vertical family : B0 = γ and A0 = D0 = E0 = 0.
• We use a parallel shooting method to compute the periodic orbits.
• We have done this for different values of δ.

WSIMS – p.14/36

Planar Family of Periodic Orbits

• We have computed the planar family for δ = 0. (i.e. sail perpendicular to

Sun).

• 0.015
• 0.01
• 0.005

0.005 0.01 0.015

• 0.992
• 0.99
• 0.988
• 0.986
• 0.984
• 0.982
• 0.98

z x delta = 0 C x S S x S

WSIMS – p.15/36

Continuation of the Planar Family

• We have computed the planar family for δ = 0.001.
• 0.015
• 0.01
• 0.005

0.005 0.01 0.015

• 0.992
• 0.99
• 0.988
• 0.986
• 0.984
• 0.982
• 0.98

z x delta =10^-3 C x S S x S

WSIMS – p.16/36

Continuation of the Planar Family

• 0.015
• 0.01
• 0.005

0.005 0.01 0.015

• 0.992
• 0.99
• 0.988
• 0.986
• 0.984
• 0.982
• 0.98

z x C x S S x S

WSIMS – p.17/36

Planar Family of Periodic Orbits

Periodic Orbits for δ = 0.

• 1
• 0.995
• 0.99
• 0.985
• 0.98
• 0.975
• 0.97
• 0.965
• 0.96
• 0.955-0.06
• 0.04
• 0.02 0 0.02

0.04 0.06

• 1
• 0.5

0.5 1 x y

• 1 -0.995
• 0.99
• 0.985
• 0.98
• 0.975
• 0.97
• 0.965
• 0.03
• 0.02
• 0.01 0 0.01

0.02 0.03

• 0.015
• 0.01
• 0.005

0.005 0.01 0.015 z x y z

Periodic Orbits for δ = 0.01.

• 1 -0.995-0.99-0.985-0.98-0.975-0.97-0.965
• 0.03
• 0.02
• 0.01

0.01 0.02 0.03

• 0.015
• 0.01
• 0.005

0.005 0.01 0.015 z x y z

• 1.005-1-0.995
• 0.99
• 0.985
• 0.98
• 0.975
• 0.97
• 0.965
• 0.96
• 0.04
• 0.03
• 0.02
• 0.01

0.01 0.02 0.03 0.04

• 0.015
• 0.01
• 0.005

0.005 0.01 0.015 z x y z

WSIMS – p.18/36

Planar Family of Periodic Orbits

Familly for δ = 0

• 1-0.995
• 0.99
• 0.985
• 0.98
• 0.975
• 0.97
• 0.965
• 0.96
• 0.955
• 0.05
• 0.04
• 0.03
• 0.02
• 0.01

0.01 0.02 0.03 0.04 0.05

• 0.015
• 0.01
• 0.005

0.005 0.01 0.015 z x y z

Familly for δ = 0.01

• 1.005-1-0.995
• 0.99
• 0.985
• 0.98
• 0.975
• 0.97
• 0.965
• 0.96
• 0.04
• 0.03
• 0.02
• 0.01

0.01 0.02 0.03 0.04

• 0.015
• 0.01
• 0.005

0.005 0.01 0.015 z x y z

WSIMS – p.19/36

Vertical Family of Periodic Orbits

• 0.9818
• 0.9816
• 0.9814
• 0.9812
• 0.981
• 0.9808
• 0.9806
• 0.9804
• 0.9802
• 0.98
• 0.9798
• 0.0008
• 0.0006
• 0.0004
• 0.0002

0.0002 0.0004 0.0006 0.0008

• 0.03
• 0.02
• 0.01

0.01 0.02 0.03 z delta = 0 x y z

• 0.9818
• 0.9816
• 0.9814
• 0.9812
• 0.981
• 0.9808
• 0.9806
• 0.9804
• 0.9802
• 0.98
• 0.9798
• 0.0008
• 0.0006
• 0.0004
• 0.0002

0.0002 0.0004 0.0006 0.0008

• 0.03
• 0.02
• 0.01

0.01 0.02 0.03 z delta = 0.001 x y z

• 0.9818
• 0.9816
• 0.9814
• 0.9812
• 0.981
• 0.9808
• 0.9806
• 0.9804
• 0.9802
• 0.98
• 0.9798
• 0.0008
• 0.0006
• 0.0004
• 0.0002

0.0002 0.0004 0.0006 0.0008

• 0.03
• 0.02
• 0.01

0.01 0.02 0.03 z delta = 0.005 x y z

• 0.9818
• 0.9816
• 0.9814
• 0.9812
• 0.981
• 0.9808
• 0.9806
• 0.9804
• 0.9802
• 0.98
• 0.9798
• 0.001
• 0.0008
• 0.0006
• 0.0004
• 0.0002

0.0002 0.0004 0.0006 0.0008 0.001

• 0.03
• 0.02
• 0.01

0.01 0.02 0.03 z delta = 0.01 x y z

WSIMS – p.20/36

Reduction to the Centre Manifold

Using an appropriate linear transformation, the equations around the fixed point can be written as,

˙ x = Ax + f(x, y), x ∈ R4, ˙ y = By + g(x, y), y ∈ R2,

where A is an elliptic matrix and B an hyperbolic one, and

f(0, 0) = g(0, 0) = 0 and Df(0, 0) = Dg(0, 0) = 0.

• We want to obtain y = v(x), with v(0) = 0, Dv(0) = 0, the local

expression of the centre manifold.

• The flow restricted to the invariant manifold is

˙ x = Ax + f(x, v(x)).

WSIMS – p.21/36

Approximating the Centre Manifold

To find y = v(x) we substitute this expression on the differential equations. So v(x) must satisfy,

Dv(x)Ax − Bv(x) = g(x, v(x)) − Dv(x)f(x, v(x)).

(1) We take,

v(x) =  

|k|≥2

v1,kxk,

• |k|≥2

v2,kxk   , k ∈ (N ∪ {0})4,

its expansion as power series. The left hand side is a linear operator w.r.t v(x) and the right hand side is non-linear.

WSIMS – p.22/36

Approximating the Centre Manifold

The left hand side of equation (1),

L(x) = Dv(x)Ax − Bv(x),

diagonalizes if A and B are diagonal. In particular, if A = diag(iω1, −iω1, iω2, −iω2) and B = diag(λ, −λ) then,

L(x) =     

• |k|≥2

(iω1k1 − iω1k2 + iω2k3 − iω2k4 − λ)v1,kxk

• |k|≥2

(iω1k1 − iω1k2 + iω2k3 − iω2k4 + λ)v2,kxk      .

WSIMS – p.23/36

Approximating the Centre Manifold

The right hand side of equation (1),

h(x) = g(x, v(x)) − Dv(x)f(x, v(x)),

can be expressed as,

h(x) =  

|k|≥2

h1,kxk ,

• |k|≥2

h2,kxk  

T

,

where hi,k depend on vi,j in a known way (i = 1, 2).

• It can be seen that for a fixed degree |k| = n, the hi,k depend only on the

vi,j such that |j| < n.

WSIMS – p.24/36

Approximating the Centre Manifold

Now we can solve equation (1) in an iterative way, equalising the left and the right hand side degree by degree. We have to solve a diagonal system at each degree. Notice:

• It is important to have a fast way to find the hi,k to get up to high degrees.
• We do not recommend to expand f(x, y) y g(x, y), and then compose

with y = v(x). One should find other alternative ways, faster in terms of computational time.

• The matrixes A and B don’t have to be diagonal, but then one must solve

a larger linear system at each degree.

WSIMS – p.25/36

On the efficient computation of hi,j

We recall that the equations of motion for α = 0 are,

¨ x = 2 ˙ y + x − κs x − µ r3

ps

− κe x + 1 − µ r3

pe

+ κsail z(x − µ) r3

psr2

, ¨ y = −2 ˙ x + y − κs r3

ps

+ κe r3

pe

• y + κsail zy

r3

psr2 ,

¨ z = − κs r3

ps

+ κe r3

pe

• z − κsail r2

r3

ps

,

where κs = (1 − µ)(1 − β cos3 α), κe = µ, κsail = β(1 − µ) cos2 α sin α.

• To expand the equations of motion we use the Legendre polynomials.

WSIMS – p.26/36

On the efficient computation of hi,j

For example:

• 1/rps can be expanded as,
• n≥0

cnTn(x, y, z),

where the Tn(x, y, z) are homogeneous polynomials of degree n that are computed in a recurrent way.

Tn = 2n − 1 n xTn−1 − n − 1 n (x2 + y2 + z2)Tn−2,

with T0 = 1,

T1 = x.

WSIMS – p.27/36

On the efficient computation of hi,j

• The functions f(¯

x, ¯ y) and g(¯ x, ¯ y) can be computed in a reccurrent way

as they are found after applying a linear transformation to the expansion of the system.

• Composing these recurrences with v(¯

x) we can compute the expansions

• f f(¯

x, v(¯ x)) and g(¯ x, v(¯ x)) in a recurrent way and so for the hi,j.

For example:

T0 = 1, T1 = x(¯ x, v(¯ x)), Tn = 2n − 1 n x(¯ x, v(¯ x))Tn−1 − n − 1 n

• x(¯

x, v(¯ x))2 + y(¯ x, v(¯ x))2 + z(¯ x, v(¯ x))2 Tn−2.

WSIMS – p.28/36

Validation Test

• Given an initial condition v0, we denote v1 and ˜

v1 to the integration at time t = 0.1 of v0 on the centre manifold and the complete system respectively.

• The error behaves as: |˜

v1 − v1| = chn+1, where h is the distance to the

• rigin of v0.
• If we consider the centre manifold up to degree 8:

h |˜ v1 − v1| n + 1 0.04 2.3643547906724647e − 15 0.08 1.2618898774811476e − 12 9.059923 0.16 6.9534006796827247e − 10 9.105988 0.32 3.9879163406855978e − 07 9.163700

WSIMS – p.29/36

Results for δ = 0

We have computed the reduction of to the centre manifold around Sub-L1 up to degree 32. (it takes 17min of CPU time)

• After this reduction we are in a four dimensional phase space

(x1, x2, x3, x4).

• We fix a Poincar´

e section x3 = 0 to reduce the system to a three dimensional phase space.

• We have taken several initial conditions and computed their successive

images on the Poincar´ e section.

WSIMS – p.30/36

Results for δ = 0

• 0.8
• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8 1

• 1.5
• 1
• 0.5

0.5 1 1.5 x4 = 0.01

• 0.8
• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8

• 1.5
• 1
• 0.5

0.5 1 1.5 x4 = 0.05

• 0.8
• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8

• 1.5
• 1
• 0.5

0.5 1 1.5 x4 = 0.11

• 0.8
• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8

• 1.5
• 1
• 0.5

0.5 1 1.5 x4 = 0.17

WSIMS – p.31/36

Results for δ = 0

• 1.5
• 1
• 0.5

0.5 1 1.5-0.8

• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8 1 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 x4 = 0.01

• 1.5
• 1
• 0.5

0.5 1 1.5-0.8

• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 x4 = 0.05

• 1.5
• 1
• 0.5

0.5 1 1.5-0.8

• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 x4 = 0.11

• 1.5
• 1
• 0.5

0.5 1 1.5-0.8

• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 x4 = 0.17

WSIMS – p.32/36

Results for δ = 0 (for a fixed energy level)

• 1
• 0.5

0.5 1

• 1
• 0.5

0.5 1 h = 0.09

• 1
• 0.5

0.5 1

• 1
• 0.5

0.5 1 h = 0.11

• 1
• 0.5

0.5 1

• 1
• 0.5

0.5 1 h = 0.13

• 1
• 0.5

0.5 1

• 1
• 0.5

0.5 1 h = 0.15

WSIMS – p.33/36

Results for δ = 0.05

• 1.5
• 1
• 0.5

0.5 1 1.5

• 1.5
• 1
• 0.5

0.5 1 1.5 x4 = 0.28

• 1.5
• 1
• 0.5

0.5 1 1.5

• 1.5
• 1
• 0.5

0.5 1 1.5 x4 = 0.49

• 1.5
• 1
• 0.5

0.5 1 1.5

• 1
• 0.5

0.5 1 1.5 x4 = 0.7

• 1.5
• 1
• 0.5

0.5 1 1.5

• 1
• 0.5

0.5 1 1.5 x4 = 0.84

WSIMS – p.34/36

Results for δ = 0.05

• 1.5
• 1
• 0.5

0.5 1 1.5

• 1.5
• 1
• 0.5

0.5 1 1.5 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 x4 = 0.28

• 1.5
• 1
• 0.5

0.5 1 1.5

• 1.5
• 1
• 0.5

0.5 1 1.5 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 x4 = 0.49

• 1
• 0.5

0.5 1 1.5-1.5

• 1-0.5

0 0.5 1 1.5 0.5 0.6 0.7 0.8 0.9 1 1.1 x4 = 0.7

• 1
• 0.5

0.5 1 1.5-1.5

• 1-0.5

0 0.5 1 1.5 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 x4 = 0.84

WSIMS – p.35/36

The End

Thank You !!!

WSIMS – p.36/36