2.1 Geometry and coordinates

As shown in Fig. 2, let $x,y,z$ be directed along a set of dextral, orthogonal axes fixed in the reference frame centred at the Sun-Earth barycentre and rotating with the Earth. For all practical purposes, the origin may be taken as the centre of the Sun. The $x$ -axis is directed from the Sun to the Earth, the $z$ -axis is along the rotational angular velocity direction and the $y$ -axis is directed along the mutual perpendicular direction.

Figure 2.

Schematic diagram of the Sun, the Earth and the shade for showing the axis convention and the frame of reference. The shade is shown as a point for clarity. The figure is not to scale.

The shade is assumed to be a paraboloid of revolution and $q,d,o$ are associated with the shade’s principal axes with the origin at the apex and $o$ being associated with the axis of symmetry. The shade is treated as a surface; that is, a body of zero thickness. A schematic version is shown in Fig. 3. The equation for the shade surface takes the form

(1) \begin{align}o=\alpha \left({q}^{2}+{d}^{2}\right)\end{align}

where $\alpha $ is the curvature parameter. Let the shade have a radius $R$ so that it exists in the region

(2) \begin{align}{q}^{2}+{d}^{2}\le {R}^{2}\end{align}

Figure 3.

The shade. It is shown as opening outwards from the Sun since that is essential for stability, as will be demonstrated later. This figure is not to scale.

An integration can be used to show that the CoM lies at the point $(\mathrm{0,0},{o}_{B})$ where

(3) \begin{align}{o}_{B}=\frac{{\big(4{\alpha }^{2}{R}^{2}+1\big)}^{3/2}\left(6{\alpha }^{2}{R}^{2}-1\right)+1}{10\alpha \left[{\left(4{\alpha }^{2}{R}^{2}+1\right)}^{3/2}-1\right]}\end{align}

In the limit of small $\alpha $ , when the shade is reduced to a disk, one finds

(4) \begin{align}{o}_{B}=\frac{\alpha {R}^{2}}{2}\end{align}

Any point on the shade is uniquely identified by the two coordinates $q,d$ . The infinitesimal area element at point $(q,d)$ is given by

(5) \begin{align}\mathrm{d}A=\sqrt{4{\alpha }^{2}{q}^{2}+4{\alpha }^{2}{d}^{2}+1}\mathrm{d}q\mathrm{d}d\end{align}

Let the surface density – that is, mass per unit area of the material – be $\sigma $ so that the infinitesimal mass element is

(6) \begin{align}\mathrm{d}m=\sigma \sqrt{4{\alpha }^{2}{q}^{2}+4{\alpha }^{2}{d}^{2}+1}\mathrm{d}q\mathrm{d}d\end{align}

The reference configuration of the shade features the apex (origin of the qdo basis) at the point $({x}_{0},\mathrm{0,0})$ where ${x}_{0}$ is determined from the equilibrium force balance. Furthermore, the $o$ -axis of the shade is coincident with the $x$ -axis so that $o,q,d$ correspond to $x,y,z$ in the reference configuration (the shade is chosen to open out away from the Sun). In the displaced configuration, the apex is translated by distances ${\eta }_{x}$ , ${\eta }_{y}$ , and ${\eta }_{z}$ along the three axes. Further, to model the angular perturbations, the shade is rotated first through angle $\varphi $ about the $y$ -axis, then through angle $\theta $ about the new $z$ -axis and finally through angle $\psi $ about the new $x$ -axis (shade axis of symmetry).

2.2 Structure of the equations of motion

The translational equations of motion of the shade are given by

(7) \begin{align}\begin{array}{ll}m{\ddot{\eta }}_{x} {} ={F}_{x}\\[3pt] m{\ddot{\eta }}_{y} {} ={F}_{y}\\[3pt] m{\ddot{\eta }}_{z} {} ={F}_{z}\end{array}\end{align}

where $m$ is the mass of the shade and ${F}_{x,y,z}$ are the resultant forces along the three fundamental directions. At this step, the displacements ${\eta }_{x}$ , ${\eta }_{y}$ and ${\eta }_{z}$ can be arbitrary, although they will be assumed to be small in magnitude later on.

The rotational dynamics of the shade are given by Euler’s equations, with the angular velocity being related to the Euler angle derivatives through the usual Tait–Bryan conventions [Reference Greenwood24]. The angular momentum balance equations are given by

(8) \begin{align}\begin{array}{ll} {} {I}_{q}{\dot{\omega }}_{q}+({I}_{o}-{I}_{d}){\omega }_{o}{\omega }_{d}={\tau }_{q}\\[3pt] {} {I}_{d}{\dot{\omega }}_{d}+({I}_{q}-{I}_{o}){\omega }_{q}{\omega }_{o}={\tau }_{d}\\[3pt] {} {I}_{o}{\dot{\omega }}_{o}+({I}_{d}-{I}_{q}){\omega }_{d}{\omega }_{q}={\tau }_{o}\end{array}\end{align}

with $\tau $ used to denote the torque and the angular velocity being given by

(9) \begin{align}\omega =\big(\dot{\theta }\,\mathrm{sin}\,\varphi +\dot{\psi }\,\mathrm{cos}\,\theta\, \mathrm{cos}\,\varphi \big){\hat{\mathbf{n}}}_{x}+\big(\dot{\varphi }+\dot{\psi }\,\mathrm{sin}\,\theta \big){\hat{\mathbf{n}}}_{y}+\big(\dot{\theta }\,\mathrm{cos}\,\varphi -\dot{\psi }\,\mathrm{cos}\,\theta \,\mathrm{sin}\,\varphi \big){\hat{\mathbf{n}}}_{z}\end{align}

where the orthogonal set of unit vectors ${\hat{\mathbf{n}}}_{x}$ , ${\hat{\mathbf{n}}}_{y}$ and ${\hat{\mathbf{n}}}_{z}$ are fixed at the shade centre. When the angles and angular velocities are small in magnitude, these equations can be simplified considerably. Since only small perturbations of the angles from the reference configuration are realistic in this problem, this assumption is immediately made. Then, the gyroscopic coupling terms in Euler’s equations disappear, and one is left with

(10) \begin{align}{I}_{o}{\dot{\omega }}_{o} {} &={\tau }_{o}\nonumber\\[3pt] I{\dot{\omega }}_{q} {} &={\tau }_{q}\\[3pt] I{\dot{\omega }}_{d} {} &={\tau }_{d}\nonumber\end{align}

where the moments of inertia ${I}_{q}$ and ${I}_{d}$ are related as ${I}_{q}={I}_{d}=I$ from symmetry. Due to the perfect symmetry of the shade, the angle $\psi $ does not influence the dynamics of $\varphi $ and $\theta $ . Hence, it can be assumed that $\psi =0$ throughout. For any vector $\mathbf{V}$ , the associated components can be transformed from the $xyz$ basis to the $oqd$ basis as per the rotation matrix

(11) \begin{align}\left[\begin{array}{l}{V}_{o}\\ {V}_{q}\\ {V}_{d}\\ \end{array}\right]={\mathbf{R}}_{xy}^{oqd}\left[\begin{array}{l}{V}_{x}\\ {V}_{y}\\ {V}_{z}\\ \end{array}\right]\end{align}

where

(12) \begin{align}{\mathbf{R}}_{xyz}^{oqd}=\left[\begin{array}{c@{\quad}c@{\quad}c}1& \theta& -\varphi \\ -\theta & 1& 0\\ \varphi& 0& 1\end{array}\right]\end{align}

This matrix is orthogonal up to $O(\varphi ,\theta )$ and is antisymmetric. By using this matrix, one can resolve the angular derivatives to form the angular velocity vector. $\dot{\varphi }$ , which is about the $y$ -axis, can be resolved into the $oqd$ bases as $\dot{\varphi }\left(\theta {\hat{\mathbf{n}}}_{o}+{\hat{\mathbf{n}}}_{q}\right)$ . After neglecting the second-order term, this vector can be treated as being directed about the $q$ -axis. $\dot{\theta },$ on the other hand, is about the new $z$ -axis, which is the $d$ -axis. Resolving it as $\dot{\theta }\left(\varphi {\hat{\mathbf{n}}}_{x}+{\hat{\mathbf{n}}}_{z}\right)$ , it can just as well be treated as the old $z$ -axis. Thus, $\dot{\varphi }$ and $\dot{\theta }$ can interchangeably be assumed to be along the $y$ -axis or $q$ -axis and the $z$ -axis or $d$ -axis, respectively. This enables one to write as

(13) \begin{align}\begin{array}{ll}I\ddot{\varphi } {} ={\tau }_{y}\\ I\ddot{\theta } {} ={\tau }_{z}\end{array}\end{align}

On the left-hand side (LHS) of Equation (13), the angle derivatives are assumed to be about the body principal axes, while on the right-hand side (RHS) they have been treated as being about the space-fixed axes. Finally, for the angle $\psi $ , one can obtain in a similar manner that

(14) \begin{align}{I}_{o}\ddot{\psi }={\tau }_{x}\end{align}

Equations (13) and (14) follow from linearisation in both the small angles and torques around the equilibrium point and neglecting second- and higher-order terms following the use of the rotation matrix in Equation (12) for conversion between the axes.

The end result is intuitive; that is, in the limit of infinitesimal rotations, the composite three-dimensional rotation decouples into three separate two-dimensional rotations.

The axial as well as crosswise moments of inertia can be evaluated through elementary integration. The results are given by

(15) \begin{align}{I}_{o}=\frac{\pi \sigma }{60{\alpha }^{4}}\big[{(4{\alpha }^{2}{R}^{2}+1)}^{3/2}(6{\alpha }^{2}{R}^{2}-1)+1\big]\end{align}

and

(16) \begin{align}I=\frac{\pi \sigma }{\alpha }\left[\begin{array}{ll} {} \dfrac{{\left(4{\alpha }^{2}{R}^{2}+1\right)}^{3/2}\left(30{\alpha }^{4}{R}^{4}-6{\alpha }^{2}{R}^{2}+1\right)-1}{420{\alpha }^{3}}-2{o}_{B}\dfrac{{\left(4{\alpha }^{2}{R}^{2}+1\right)}^{3/2}\left(6{\alpha }^{2}{R}^{2}-1\right)+1}{60{\alpha }^{2}}\\[9pt] {} +{o}_{B}^{2}\dfrac{{\left(4{\alpha }^{2}{R}^{2}+1\right)}^{3/2}-1}{6\alpha }\end{array}\right]\end{align}

In the limit of $\alpha =0$ , the shade is reduced from a paraboloid to a disk; for small $\alpha $ , the following are obtained:

(17) \begin{align}{I}_{o}=\frac{\pi \sigma {R}^{4}}{2}+\frac{2\pi \sigma {\alpha }^{2}{R}^{6}}{3}\end{align}

and

(18) \begin{align}I=\frac{\pi \sigma {R}^{4}}{4}+\frac{5\pi \sigma {\alpha }^{2}{R}^{6}}{12}\end{align}

Recognising that $\pi {R}^{2}\sigma $ is the mass of the disk when $\alpha =0$ , a consistency check can be made.

2.3 Individual forces and torques

As already mentioned, the shade is subjected to five force fields, namely the Sun’s gravity field, Earth’s gravity field, radiation pressure, centrifugal force and Coriolis force. To make the calculation analytically tractable (up to quadratures), a linearisation process is used. Since the intended size of the shade (radius about 1000 km) is orders of magnitude smaller than the Earth-Sun distance ( $1.5\times 1{0}^{8}$ km), the force fields are linearised in the region containing the shade. Since shade locations close to the $x$ -axis are of interest, this linearisation is performed about the point $(x,y,z)=({x}_{0},\mathrm{0,0})$ . For any force field $\mathbf{F}\left(x,y,z\right)$ , one has at the point $({x}_{0},\mathrm{0,0})+(\xi ,\eta ,\zeta )$ , the expansion

(19) \begin{align}{F}_{x}\left({x}_{0}+\xi ,\eta ,\zeta \right)={F}_{x}\left({x}_{0},\mathrm{0,0}\right)+\frac{\partial {F}_{x}}{\partial x}\xi +\frac{\partial {F}_{x}}{\partial y}\eta +\frac{\partial {F}_{x}}{\partial z}\zeta \end{align}

and similar equations for ${F}_{y}$ and ${F}_{z}$ .

Next, each of the five forces acting on the shade is analysed. In each case, a computation is made of the force contribution on the infinitesimal element of the shade located at the point $(q,d)$ as per Equations (5) and (6). For brevity, the symbol ${w}_{0}$ is used to denote $\sqrt{4{\alpha }^{2}{q}^{2}+4{\alpha }^{2}{d}^{2}+1}$ . All computations have been performed using the software WXMaxima, a computer algebra system.

2.3.1 Sun’s gravity field

The force field has the form

(20) \begin{align}{\mathbf{F}}_{SG}=\frac{-G{M}_{S}\left({\hat{\mathbf{n}}}_{x}+y{\hat{\mathbf{n}}}_{y}+z{\hat{\mathbf{n}}}_{z}\right)}{{\left({x}^{2}+{y}^{2}+{z}^{2}\right)}^{3/2}}\end{align}

Performing the linearisation, the following is found for the force on the infinitesimal element:

(21) \begin{align} {} \mathrm{d}{\mathbf{F}}_{SG}&=\mathrm{d}q\mathrm{d}d\times \nonumber\\[4pt] {}& \quad {\hat{\mathbf{n}}}_{x}\left[-\dfrac{G{M}_{S}\sigma {w}_{0}}{{x}_{0}^{3}}\left({x}_{0}-2\alpha ({q}^{2}+{d}^{2})\right)+\dfrac{2G{M}_{S}\sigma {w}_{0}}{{x}_{0}^{3}}{\eta }_{x}\right]\nonumber\\[5pt] {}& \quad + {\hat{\mathbf{n}}}_{y}\left[-\dfrac{G{M}_{S}\sigma {w}_{0}}{{x}_{0}^{3}}{\eta }_{y}-\dfrac{G{M}_{S}\sigma {w}_{0}}{{x}_{0}^{3}}\alpha \left({q}^{2}+{d}^{2}\right)\theta \right]\\[5pt] {}& \quad + {\hat{\mathbf{n}}}_{z}\left[-\dfrac{G{M}_{S}\sigma {w}_{0}}{{x}_{0}^{3}}{\eta }_{z}+\dfrac{G{M}_{S}\sigma {w}_{0}}{{x}_{0}^{3}}\alpha \left({q}^{2}+{d}^{2}\right)\varphi \right]\nonumber\end{align}

2.3.2 Earth’s gravity field

The force field has the form

(22) \begin{align}{\mathbf{F}}_{EG}=-\frac{G{M}_{E}\left((x-{r}_{E}){\hat{\mathbf{n}}}_{x}+y{\hat{\mathbf{n}}}_{y}+z{\hat{\mathbf{n}}}_{z}\right)}{{\left({\left(x-{r}_{E}\right)}^{2}+{y}^{2}+{z}^{2}\right)}^{3/2}}\end{align}

and the linearisation yields the element-wise force

(23) \begin{align} \mathrm{d}{\mathbf{F}}_{EG}&=\mathrm{d}q\mathrm{d}d\times \nonumber\\[5pt] & \quad {\hat{\mathbf{n}}}_{x}\left[-\dfrac{G{M}_{E}\sigma {w}_{0}}{{\left({r}_{E}-{x}_{0}\right)}^{3}}\left(-\left({r}_{E}-{x}_{0}\right)-2\alpha ({q}^{2}+{d}^{2})\right)+\dfrac{2G{M}_{E}\sigma {w}_{0}}{{\left({r}_{E}-{x}_{0}\right)}^{3}}{\eta }_{x}\right]\nonumber\\[5pt] & {}\quad + {\hat{\mathbf{n}}}_{y}\left[-\dfrac{G{M}_{E}\sigma {w}_{0}}{{\left({r}_{E}-{x}_{0}\right)}^{3}}{\eta }_{y}-\dfrac{G{M}_{E}\sigma {w}_{0}}{{\left({r}_{E}-{x}_{0}\right)}^{3}}\alpha \left({q}^{2}+{d}^{2}\right)\theta \right]\\[5pt] & {}\quad + {\hat{\mathbf{n}}}_{z}\left[-\dfrac{G{M}_{E}\sigma {w}_{0}}{{\left({r}_{E}-{x}_{0}\right)}^{3}}{\eta }_{z}+\dfrac{G{M}_{E}\sigma {w}_{0}}{{\left({r}_{E}-{x}_{0}\right)}^{3}}\alpha \left({q}^{2}+{d}^{2}\right)\varphi \right]\nonumber\end{align}

2.3.3 Solar radiation pressure

This force arises from the impingement of the Sun’s photons on the shade surface. The photons themselves stream radially outward from the Sun. On every infinitesimal area of the shade, they transfer to the area element their component of momentum along the normal to the element. This transferred momentum is always along the outward normal of the area element.

The raw photon momentum is

(24) \begin{align}=\frac{k}{{r}^{2}}{\hat{\mathbf{e}}}_{r}\end{align}

where $k$ is a constant. At the point $(q,d)$ , the normal ${\hat{\mathbf{e}}}_{n}$ to the shade surface is

(25) \begin{align}{\hat{\mathbf{e}}}_{n}=\frac{-{\hat{\mathbf{n}}}_{o}+2\alpha {\hat{\mathbf{n}}}_{q}+2\alpha {\hat{\mathbf{n}}}_{d}}{{w}_{0}}\end{align}

which can be expressed in the $xyz$ bases by using the transpose of the rotation matrix Equation (12). The momentum component normal to the shade is given by the dot product of Equation (24) and Equation (25), and this component is directed along the outward normal, which is the negative of Equation (25) since the positive $o$ -axis is directed away from the Sun. After performing the computations, one finds the radiation pressure on the infinitesimal shade element to be

(26) \begin{align} {} \mathrm{d}{\mathbf{F}}_{RP}&=\frac{k}{{w}_{0}{x}_{0}^{3}}\mathrm{d}q\mathrm{d}d\times \nonumber\\[3pt] {}&\quad {\hat{\mathbf{n}}}_{x}\left[\begin{array}{ll} {} -{x}_{0}+4\alpha \left({q}^{2}+{d}^{2}\right)+2{\eta }_{x}+2\alpha q{\eta }_{}+2\alpha d{\eta }_{z}\nonumber\\[3pt] {} + \theta \left(-4\alpha q{x}_{0}+14{\alpha }^{2}{q}^{3}+14{\alpha }^{2}q{d}^{2}-3q\right)+\varphi \left(4\alpha d{x}_{0}-14{\alpha }^{2}{q}^{2}d-14{\alpha }^{2}{d}^{3}+3d\right) {} \end{array}\right]\nonumber\\[3pt] &\quad {}+ {\hat{\mathbf{n}}}_{y}\left[\begin{array}{ll} {} -2\alpha q{x}_{0}+8{\alpha }^{2}{q}^{3}+8{\alpha }^{2}q{d}^{2}+4\alpha q{\eta }_{x}+4{\alpha }^{2}{q}^{2}{\eta }_{y}+4{\alpha }^{2}qd{\eta }_{z}-\nonumber\\[3pt] {} \theta \left(-{x}_{0}+4{\alpha }^{2}{q}^{2}{x}_{0}-12{\alpha }^{3}{q}^{4}-12{\alpha }^{3}{q}^{2}{d}^{2}+10\alpha {q}^{2}+4\alpha {d}^{2}\right)\\[3pt] {} -\varphi \left(-4{\alpha }^{2}qd{x}_{0}+12{\alpha }^{3}{q}^{3}d+12{\alpha }^{3}q{d}^{3}-6\alpha qd\right) {} \end{array}\right]\nonumber\\[3pt] & \quad + {\hat{\mathbf{n}}}_{z}\left[\begin{array}{ll} {} -2\alpha d{x}_{0}+8{\alpha }^{2}{q}^{2}d+8{\alpha }^{2}{d}^{3}+4\alpha d{\eta }_{x}+4{\alpha }^{2}qd{\eta }_{y}+4{\alpha }^{2}{d}^{2}{\eta }_{z}-\nonumber\\[3pt] {} \theta \left(4{\alpha }^{2}qd{x}_{0}-12{\alpha }^{3}{q}^{3}d-12{\alpha }^{3}q{d}^{3}+6\alpha qd\right)-\nonumber\\[3pt] {} \varphi \left({x}_{0}-4{\alpha }^{2}{d}^{2}{x}_{0}+12{\alpha }^{3}{q}^{2}{d}^{2}+12{\alpha }^{3}{d}^{4}-4\alpha {q}^{2}-10\alpha {d}^{2}\right) {} \end{array}\right]\end{align}

2.3.4 Centrifugal force

The field is given by

(27) \begin{align}{\mathbf{F}}_{CF}={\mathrm{\Omega }}^{2}\left(x{\hat{\mathbf{n}}}_{x}+y{\hat{\mathbf{n}}}_{y}\right)\end{align}

and the elemental force follows as

(28) \begin{align}\begin{array}{ll} {} \mathrm{d}{\mathbf{F}}_{CF}=\mathrm{d}q\mathrm{d}d\times \\[5pt] {} {\hat{\mathbf{n}}}_{x}\left[\sigma {w}_{0}{\mathrm{\Omega }}^{2}\left({x}_{0}+\alpha ({q}^{2}+{d}^{2})\right)+\sigma {w}_{0}{\mathrm{\Omega }}^{2}{\eta }_{x}\right]+\\[8pt] {} {\hat{\mathbf{n}}}_{y}\left[\sigma {w}_{0}{\mathrm{\Omega }}^{2}{\eta }_{y}+\sigma {w}_{0}{\mathrm{\Omega }}^{2}\alpha \left({q}^{2}+{d}^{2}\right)\theta \right]+{\hat{\mathbf{n}}}_{z}\left[\sigma {w}_{0}{\mathrm{\Omega }}^{2}{\eta }_{z}-\sigma {w}_{0}{\mathrm{\Omega }}^{2}\alpha \left({q}^{2}+{d}^{2}\right)\varphi \right]\end{array}\end{align}

2.3.5 Coriolis force

The associated force follows from

(29) \begin{align}{\mathbf{F}}_{COR}=-2\mathbf{\Omega }\times \mathbf{v}\end{align}

and the elemental force is determined as

(30) \begin{align}\begin{array}{ll} {} \mathrm{d}{\mathbf{F}}_{COR}=2\sigma {w}_{0}\mathrm{\Omega }\mathrm{d}q\mathrm{d}d\times \\[5pt] {} {\hat{\mathbf{n}}}_{x}\left[{\dot{\eta }}_{y}+{x}_{p}\dot{\theta }\right]+{\hat{\mathbf{n}}}_{y}\left[-{\dot{\eta }}_{x}+\left(\dot{\varphi }{z}_{p}-\dot{\theta }{y}_{p}\right)\right]\end{array}\end{align}

where ${x}_{p}$ , ${y}_{p}$ and ${z}_{p}$ denote the $x$ , $y$ and $z$ coordinates from the CoM to a point on the shade. Since the rotation of the reference frame is about the $z$ -axis, ${\mathbf{F}}_{COR}$ does not have a $z$ -component.

The model has been set up in this section. In the following section, parametric studies are conducted.

3.1 Combined forces and torques

In this subsection, first, the combined forces and torques are given, based on the individual forces determined in the preceding subsection. The torques are always computed about the CoM of the shade. It is noted that the total forces and torques on the shade are obtained by integrating the element-wise expressions over the whole area of the shade. This integration takes the form

(31) \begin{align}force/torque=\underset{-R}{\overset{R}{\int }} \mathrm{d}q\underset{-\sqrt{{R}^{2}-{q}^{2}}}{\overset{\sqrt{{R}^{2}-{q}^{2}}}{\int }} \mathrm{d}df(q,d)\end{align}

and the two integration variables can be interchanged as well. If $f$ is antisymmetric in either variable, then the integration with respect to that variable will vanish. Hence, antisymmetric terms in $q$ and $d$ can be ignored while writing the force and torque expressions. With this stated, the following forms are found for the forces.

(32) \begin{align} {} {F}_{x} &=-\frac{G{M}_{S}\sigma {w}_{0}}{{x}_{0}^{3}}\left({x}_{0}-2\alpha ({q}^{2}+{d}^{2})\right)-\frac{G{M}_{E}\sigma {w}_{0}}{{\left({r}_{E}-{x}_{0}\right)}^{3}}\left({x}_{0}-2\alpha ({q}^{2}+{d}^{2})\right)+\frac{k}{{w}_{0}{x}_{0}^{3}}\left({x}_{0}-4\alpha ({q}^{2}+{d}^{2})\right)\nonumber\\ {}&\quad +\sigma {w}_{0}{\mathrm{\Omega }}^{2}\left({x}_{0}+\alpha ({q}^{2}+{d}^{2})\right)+{\eta }_{x}\left[\frac{2G{M}_{s}\sigma {w}_{0}}{{x}_{0}^{3}}+\frac{2G{M}_{E}\sigma {w}_{0}}{{\left({r}_{E}-{x}_{0}\right)}^{3}}-\frac{2k}{{w}_{0}{x}_{0}^{3}}+\sigma {w}_{0}{\mathrm{\Omega }}^{2}\right]+{\dot{\eta }}_{y}\left[2\sigma {w}_{0}\mathrm{\Omega }\right]\nonumber\\ {}&\quad +\dot{\theta }\left[2\sigma {w}_{0}\mathrm{\Omega }\left({q}^{2}+{d}^{2}\right)\right]\end{align}

(33) \begin{align} {} {F}_{y}&={\eta }_{y}\left[-\frac{G{M}_{S}\sigma {w}_{0}}{{x}_{0}^{3}}-\frac{G{M}_{E}\sigma {w}_{0}}{{\left({r}_{E}-{x}_{0}\right)}^{3}}+\frac{4k{\alpha }^{2}{q}^{2}}{{w}_{0}{x}_{0}^{3}}+\sigma {w}_{0}{\mathrm{\Omega }}^{2}\right]+{\dot{\eta }}_{x}\left[2\sigma {w}_{0}\mathrm{\Omega }\right]\nonumber\\[6pt] {} & \quad+ \theta \left[\begin{array}{ll} {} -\dfrac{G{M}_{S}\sigma {w}_{0}}{{x}_{0}^{3}}\alpha \left({q}^{2}+{d}^{2}\right)-\dfrac{G{M}_{E}\sigma {w}_{0}}{{\left({r}_{E}-{x}_{0}\right)}^{3}}\alpha \left({q}^{2}+{d}^{2}\right)\\ {} + \dfrac{k}{{w}_{0}{x}_{0}^{3}}\left(12{\alpha }^{3}{q}^{2}({q}^{2}+{d}^{2})-4{\alpha }^{2}{q}^{2}{x}_{0}-10\alpha {q}^{2}-4\alpha {d}^{2}+{x}_{0}\right)+\sigma {w}_{0}{\mathrm{\Omega }}^{2}\alpha \left({q}^{2}+{d}^{2}\right) {} \end{array}\right] {} \end{align}

(34) \begin{align} {} {F}_{z}&={\eta }_{z}\left[-\frac{G{M}_{S}\sigma {w}_{0}}{{x}_{0}^{3}}-\frac{G{M}_{E}\sigma {w}_{0}}{{\left({r}_{E}-{x}_{0}\right)}^{3}}+\frac{4k{\alpha }^{2}{d}^{2}}{{w}_{0}{x}_{0}^{3}}+\sigma {w}_{0}{\mathrm{\Omega }}^{2}\right]\nonumber\\[6pt] {} &\quad +\varphi \left[\begin{array}{ll} {} \dfrac{G{M}_{S}\sigma {w}_{0}}{{x}_{0}^{3}}\alpha \left({q}^{2}+{d}^{2}\right)+\dfrac{G{M}_{E}\sigma {w}_{0}}{{\left({r}_{E}-{x}_{0}\right)}^{3}}\alpha \left({q}^{2}+{d}^{2}\right)\\[4pt] {} +\dfrac{k}{{w}_{0}{x}_{0}^{3}}\left(-12{\alpha }^{3}{d}^{2}({q}^{2}+{d}^{2})+4{\alpha }^{2}{d}^{2}{x}_{0}+10\alpha {d}^{2}+4\alpha {q}^{2}-{x}_{0}\right)+\sigma {w}_{0}{\mathrm{\Omega }}^{2}\alpha \left({q}^{2}+{d}^{2}\right) {} \end{array}\right] {} \end{align}

One can observe that the coefficients in ${F}_{y}$ and ${F}_{z}$ are identical and/or negatives of each other, as is expected on account of the system symmetry.

For the torques, the following forms are obtained:

(35) \begin{align} {\tau }_{y}&={\eta }_{z}\left\{\begin{array}{ll} {} \dfrac{G{M}_{S}\sigma {w}_{0}}{{x}_{0}^{3}}\left[\alpha \left({q}^{2}+{d}^{2}\right)-{o}_{B}\right]+\dfrac{G{M}_{E}\sigma {w}_{0}}{{\left({r}_{E}-{x}_{0}\right)}^{3}}\left[\alpha \left({q}^{2}+{d}^{2}\right)-{o}_{B}\right]\nonumber\\[8pt] {} +\dfrac{k}{{w}_{0}{x}_{0}^{3}}\left[-4{\alpha }^{3}{d}^{2}\left({q}^{2}+{d}^{2}\right)+4{\alpha }^{2}{d}^{2}{o}_{B}-2\alpha {d}^{2}\right]+{\mathrm{\Omega }}^{2}\sigma {w}_{0}{x}_{0}\left[-\alpha \left({q}^{2}+{d}^{2}\right)+{o}_{B}\right] {} \end{array}\right\}\\[10pt] &\quad +\varphi \left\{\begin{array}{ll} {} \dfrac{G{M}_{S}\sigma {w}_{0}}{{x}_{0}^{3}}\big[-3{\alpha }^{2}{\left({q}^{2}+{d}^{2}\right)}^{2}+3\alpha \left({q}^{2}+{d}^{2}\right){o}_{B}+3{d}^{2}+{x}_{0}\left(\alpha ({q}^{2}+{d}^{2})-{o}_{B}\right)\big]\\[12pt] {}+ \dfrac{G{M}_{E}\sigma {w}_{0}}{{\left({r}_{E}-{x}_{0}\right)}^{3}}\big[-3{\alpha }^{2}{\left({q}^{2}+{d}^{2}\right)}^{2}+3\alpha \left({q}^{2}+{d}^{2}\right){o}_{B}+3{d}^{2}+\left({r}_{E}-{x}_{0}\right)\left(-\alpha ({q}^{2}+{d}^{2})+{o}_{B}\right)\big]\\[12pt] {}+ \dfrac{k}{{w}_{0}{x}_{0}^{3}}\Bigg[\begin{array}{ll} {} 12{\alpha }^{4}{d}^{2}{\left({q}^{2}+{d}^{2}\right)}^{2}-12{\alpha }^{3}{o}_{B}\left({q}^{2}{d}^{2}+{d}^{4}\right)+6\alpha {d}^{2}{o}_{B}-3{d}^{2}\\[6pt] {}+ {x}_{0}\left(-4{\alpha }^{3}{d}^{2}({q}^{2}+{d}^{2})+4{\alpha }^{2}{d}^{2}{o}_{B}-2\alpha {d}^{2}\right) {} \end{array}\Bigg]\\[14pt] +{\mathrm{\Omega }}^{2}\sigma {w}_{0}{x}_{0}\left[-\alpha \left({q}^{2}+{d}^{2}\right)+{o}_{B}\right] {} \end{array}\right\} {} \end{align}

and

(36) \begin{align}{\tau _z} &= {\eta _y}\left\{ \begin{array}{l}\dfrac{{G{M_S}\sigma {w_0}}}{{x_0^3}}\left[ { - \alpha \left( {{q^2} + {d^2}} \right) + {o_B}} \right] + \dfrac{{G{M_E}\sigma {w_0}}}{{{{\left( {{r_E} - {x_0}} \right)}^3}}}\left[ { - \alpha \left( {{q^2} + {d^2}} \right) + {o_B}} \right] \\[10pt]+ \dfrac{k}{{{w_0}x_0^3}}\left[ {4{\alpha ^3}{q^2}\left( {{q^2} + {d^2}} \right) - 4{\alpha ^2}{q^2}{o_B} + 2\alpha {q^2}} \right] + {{\rm{\Omega }}^2}\sigma {w_0}{x_0}\left[ {\alpha \left( {{q^2} + {d^2}} \right) - {o_B}} \right]\end{array} \right\}\nonumber\\[10pt]& \quad + \theta \left\{ \begin{array}{l}\dfrac{{G{M_S}\sigma {w_0}}}{{x_0^3}}\big[ { - 3{\alpha ^2}{{\left( {{q^2} + {d^2}} \right)}^2} + 3\alpha \left( {{q^2} + {d^2}} \right){o_B} + 3{q^2} + {x_0}\left( {\alpha \left( {{q^2} + {d^2}} \right) - {o_B}} \right)} \big] \\[10pt]+ \dfrac{{G{M_E}\sigma {w_0}}}{{{{\left( {{r_E} - {x_0}} \right)}^3}}}\big[ { - 3{\alpha ^2}{{\big( {{q^2} + {d^2}} \big)}^2} + 3\alpha \big( {{q^2} + {d^2}} \big){o_B} + 3{q^2} + \big( {{r_E} - {x_0}} \big)\big( { - \alpha \big( {{q^2} + {d^2}} \big) + {o_B}} \big)} \big] \\[10pt]+ \dfrac{k}{{{w_0}x_0^3}}\,\,\left[ \begin{array}{l}12{\alpha ^4}{q^2}{\left( {{q^2} + {d^2}} \right)^2} - 12{\alpha ^3}{o_B}{q^2}\left( {{q^2} + {d^2}} \right) + 6\alpha {q^2}{o_B} - 3{q^2} \\[10pt]+ {x_0}\left( { - 4{\alpha ^3}{q^2}\left( {{q^2} + {d^2}} \right) + 4{\alpha ^2}{q^2}{o_B} - 2\alpha {q^2}} \right)\end{array} \right]\\[15pt]+ {{\rm{\Omega }}^2}\sigma {w_0}{x_0}\left[ { - \alpha \left( {{q^2} + {d^2}} \right) + {o_B}} \right]\end{array} \right\}\end{align}

Once again, symmetries consistent with the setup are observed. The expression for ${\tau _x}$ is not displayed since it does not affect the system stability, as previously discussed.

3.2 System eigenvalues

The system of governing equations given by Equation (7) and Equation (13) with forces and torques as described by Equations (32)–(36) is considered. This is by definition a linear system since linearisation was employed during the derivation of Equations (32)–(36). The parameter values used are as given in Table 1.

Table 1.

Parameters and values used in the model

Table 1. Long description

A table titled ‘Parameters and values used in the model’ contains seven rows and three columns. The columns are labeled ‘Quantity’, ‘Symbol’, and ‘SI unit value’. The table lists various astronomical and physical parameters along with their symbols and SI unit values. The rows are as follows: Row 1: Quantity, Universal gravitational constant; Symbol, G; SI unit value, 6.67 times 10 to the power of negative 11 newton meter squared per kilogram squared. Row 2: Quantity, Mass of Sun; Symbol, M subscript S; SI unit value, 1.98 times 10 to the power of 30 kilograms. Row 3: Quantity, Mass of Earth; Symbol, M subscript E; SI unit value, 5.97 times 10 to the power of 24 kilograms. Row 4: Quantity, Earth’s orbiting radius; Symbol, r subscript E; SI unit value, 1.49 times 10 to the power of 11 meters. Row 5: Quantity, Angular velocity; Symbol, omega; SI unit value, 1.991 times 10 to the power of negative 7 radians per second. Row 6: Quantity, Radius of the shade; Symbol, R; SI unit value, 10 to the power of 6 meters. Row 7: Quantity, Solar radiation pressure constant; Symbol, k; SI unit value, 1.01 times 10 to the power of 7 newtons.

The equilibrium position ${x_0}$ is found by equating to zero the terms not involving perturbations in Equation (32). In the absence of radiation pressure, ${x_0}$ reduces to the Lagrange point ${L_1}$ for a shade. When the radiation pressure is included, the equilibrium is shifted closer to the Sun. The eigenvalues associated with the equilibrium point are examined to determine which translational motions and which rotational motions are stable or unstable.

The two important parameters used to describe the shade are $\sigma $ and $\alpha $ . A $\sigma $ value of $765{\,\textrm{kg}}/{\textrm{k}}{{\textrm{m}}^2}$ results in a radiation pressure that can generate a force to balance the force due to Sun’s gravity for normal incidence. A sheet of paper typically weighs 50+ tons per ${\textrm{k}}{{\textrm{m}}^2}$ , so this density is two orders of magnitude lighter than paper. For $\alpha $ , a value of 0.001 means that the 1000 km radius shade has an axial extent of 1000 km.

With $\alpha = 0$ , two stable (centre, purely imaginary eigenvalue pair) and three unstable (saddle, positive and negative real eigenvalue pairs) modes are seen for all values of $\sigma $ . The two stable modes correspond to translational $y$ - and $z$ -displacements, while the unstable modes correspond to translational $x$ -displacement as well as rotational motions. For a given $\sigma $ , when $\alpha $ is increased, a change in stability of the rotational modes is observed as $\alpha $ is increased beyond a cutoff. This is shown in Fig. 4, where after fixing the value of $\sigma $ at $20,000{\,\textrm{kg}}/{\textrm{k}}{{\textrm{m}}^2}$ and $40,000\frac{{{\textrm{kg}}}}{{{\textrm{k}}{{\textrm{m}}^2}}},$ respectively, in Fig. 4a and b, $\alpha $ is varied from 0 to 0.0005, and the real parts of the eigenvalues associated with a system equilibrium point are plotted. Since the eigenvalues arise in pairs with the same real parts, in Fig. 4a and b, there are only 5 loci for the loci plotted for the real parts of the 10 eigenvalues in each of the plots.

Figure 4.

Real parts of the eigenvalues as the curvature parameter $\alpha $ is varied. There are 10 eigenvalues but there are only 5 loci, as more than 1 eigenvalue has the same real part. Simulation results for two different values of $\sigma $ have been shown.

The cases depicted in Fig. 4 are based on simulations involving paraboloid shades with uniform mass distribution. For all $\alpha $ values, a pair of real eigenvalues having the values $\pm 2.1 \times {10^4}$ is noted. This corresponds to the saddle point associated with the $x$ -displacement. The fixed points associated with the $y$ - and $z$ -displacements are centres and the corresponding eigenvalues are purely imaginary eigenvalues, with a zero real part, as shown in the diagram. In the case of $\sigma $ = $20,000{\,\textrm{kg}}/{\textrm{k}}{{\textrm{m}}^2}$ in Fig. 4a, for $\alpha $ values below $1.4 \times {10^{ - 4}}$ , a second pair of positive and negative real eigenvalues is observed to emerge. This pair, which actually occurs twice, corresponds to the two rotational modes. As $\alpha $ is increased from zero, the loci of the positive and negative eigenvalues are found to approach each other until they turn into a purely imaginary pair at $\alpha = 0.00014$ and remain so for higher values of $\alpha $ . A very similar trend is observed for $\sigma $ = $40,000{\,\textrm{kg}}/{\textrm{k}}{{\textrm{m}}^2}$ depicted in Fig. 4b, with the difference being that the cutoff value of $\alpha $ is higher ( $\alpha = 2.8 \times {10^{ - 4}}$ ) for stabilisation of the rotational modes.

In Fig. 5, the cutoff value of $\alpha $ at which the rotational modes become stable is shown for different values of $\sigma $ . It can be seen that heavier shades need greater curvature to stabilise the rotational modes. This makes sense since gravity has a destabilising effect on these modes while radiation pressure has a stabilising one.

Figure 5.

Transition curve in the $\sigma $ – $\alpha $ plane.

A set of illustrations is provided in Figs 6 and 7a and b to further elucidate the effects of radiation pressure and gravity with regard to stabilisation of rotational modes. The centre of pressure is computed by using the formula $\,{\vec {\!r}_{COP}} = {{\smallint \,\vec {\!r}d{F_{rp}}} \over {\smallint d{F_{rp}}}}$ , where $\,{\vec {\!r}_{COP}}$ is the centre of radiation pressure on the paraboloid shade, $d{F_{rp}}$ is the differential force due to radiation pressure at at a position $\vec r$ on the shade. The net resultant radiation force will pass through the centre of pressure located at $\,{\vec {\!r}_{COP}}$ .

Figure 6.

An aligned configuration of the paraboloid shade, for depiction of the locations of the CoP and the CoM and the torque due to the radiation pressure.

Figure 7.

A tilted configuration of the paraboloid shade, for depiction of the effects of gravitational forces from the Sun and the Earth on the shade and the resulting destabilising torque on the paraboloid shade.

As illustrated in Fig. 6, in an aligned configuration of a perfectly symmetric shade, the net resultant force is directed along the Sun-Earth axis and the $\,{\vec {\!r}_{COP}}$ and $\,{\vec {\!r}_{COM}}$ lie on the same axis. Here, $\,{\vec {\!r}_{COM}}$ is the CoM of the shade. The result is a net zero torque. In this configuration, the gravitational force due to the Sun and the Earth also does not result in a non-zero torque.

As illustrated in Fig. 7a, in a tilted configuration, the lower portion of the shade is subjected to higher radiation pressure compared to the upper portion of the shade owing to the closer orientation of the outward normal with the incident photon direction. The result is the movement of the centre of pressure (CoP) towards the lower portion of the shade. Given the net thrust $\,{\vec {\!F}_{rp,net}}$ passes through this location, the result is a counter-clockwise (CCW) torque about the CoM, which remains fixed. In this case, this is a restoring torque which provides stabilisation. A similar mechanism can be observed if the shade is tilted the other way, with the result being a clockwise (CW) restoring torque due to the radiation pressure force.

The effects of gravitational forces on the rotational mode stabilisation are illustrated through Fig. 7b. The gravitational forces comprise those due to the Sun’s gravity and the Earth’s gravity. When the shade is tilted in the configuration shown in Fig. 7b, the upper portions and the lower portions of the shade are at different distances from the Earth and the Sun, resulting in different gravitational forces, with the result being a torque in the CW direction for both the Sun’s and the Earth’s gravitational forces. This results in destabilisation of the shade.

4.1 Material reflectivity degradation

A set of investigations has been carried out to understand the effect of degradation of material reflectivity on the stability analysis of the shade. The results obtained are presented here.

The radiation pressure constant $k$ depends on the surface reflectivity of the shade, which is expected to change over the operational lifetime of the shade owing to different factors. A few simulations using optimal and slightly diminished values of $k$ for two different values of $\sigma $ were carried out to investigate the effect of material reflectivity loss on the equilibrium position ${x_0}$ . For these simulations, the radiation pressure constant $k$ was computed by using $\left( {1 - f} \right){k_0}$ , where ${k_0}$ is the nominal radiation pressure coefficient and $f$ is the fractional reflectivity loss. The results obtained for the two cases $\sigma = 20,000{\,\textrm{kg}}/{\textrm{k}}{{\textrm{m}}^2}$ , $\alpha = 0.0005{\textrm{k}}{{\textrm{m}}^{ - 1}}$ and $\sigma = 40,000{\,\textrm{kg}}/{\textrm{k}}{{\textrm{m}}^2}$ , $\alpha = 0.0005{\textrm{k}}{{\textrm{m}}^{ - 1}}$ (both stable regions) are summarised in Table 2.

From the results presented in Table 2, it can be inferred that due to the degradation of reflectivity, the net effect of the different forces results in a net sunward force at the nominal equilibrium position due to the reduced radiation pressure. To compensate for this force imbalance and realise a state of equilibrium, the position of the paraboloid shade needs to be shifted towards the Earth as outlined in Table 2. The equilibrium position shifts linearly for the reflectivity loss, with a $35{\,\textrm{km}}$ shift for $1{{\% }}$ reflectivity loss for $\sigma = 20,000{\,\textrm{kg}}/{\textrm{k}}{{\textrm{m}}^2}$ and around $17{\,\textrm{km}}$ shift for $\sigma = 40,000{\,\textrm{kg}}/{\textrm{k}}{{\textrm{m}}^2}$ .

The control architecture to compensate for this equilibrium shift involves the coordinated use of the peripheral panels to generate decoupled forces and moments as outlined in Section 4. The relevant control mode for this is the generation of a pure x-direction force that can be achieved by pitching up of two peripheral panels located ${180^ \circ }$ apart. Partially tilted peripheral panels in the reference configuration can help in adjusting the radiation pressure force in either direction in such cases. They can be rotated towards the Sun-facing orientation, increasing the effective radiation pressure contribution and restoring the outward force. Given the peripheral panels have an area of ${A_p}$ , the total shade area is $A$ , and considering the case where the nominal tilted configuration makes an angle $\beta $ with the Sun-Earth line, the maximum additional force which can be generated is $\frac{{{A_p}}}{A}\left( {1 - {\textrm{co}}{{\textrm{s}}^2}\beta } \right){F_{{\textrm{rad}}}}$ . Thus, to compensate for the reflectivity loss $f$ , we should have to satisfy $\frac{{{A_p}}}{A}\left( {1 - {\textrm{co}}{{\textrm{s}}^2}\beta } \right) \ge f$ .

4.2 Mass budget of the control system

Another stability constraint required to be considered is the mass budget of the control system. The proposed control architecture consists of actuators and hinges at the periphery, which together introduce additional mass towards the periphery. As a result, the CoM ${o_B}$ is shifted axially towards the rim of the paraboloid shade.

Figure 13.

Real parts of eigenvalues for investigation of peripheral mass distribution on paraboloid shade stability, as the curvature parameter $\alpha $ is varied. Results of simulations for three different values of the rim mass enhancement factor $f$ have been shown for areal density $\sigma = 20,000{\,\textrm{kg}}/{\textrm{k}}{{\textrm{m}}^2}$ .

To investigate the effect of ${o_B}$ shift on the shade stability, the mass distribution over the paraboloid shade is made non-uniform through the use of two parameters, namely ${r_{rim}}$ and $f$ . Here, $f$ is the rim mass enhancement factor and ${r_{rim}}$ is the rim location. The CoM ${o_B}$ is then calculated using the formula in Equation (37).

(37) \begin{align}o_{B}=\frac{\displaystyle\int_{r\leq r_{\mathrm{rim}}}\alpha\left(q^{2}+d^{2}\right)w_{0}\,\mathrm{d}q\,\mathrm{d}d+\int_{r_{\mathrm{rim}}\leq r\leq R}(1+f)\alpha\left(q^{2}+d^{2}\right)w_{0}\,\mathrm{d}q\,\mathrm{d}d}{\displaystyle\int_{r\leq r_{\mathrm{rim}}}w_{0}\,\mathrm{d}q\,\mathrm{d}d+\int_{r_{\mathrm{rim}}\leq r\leq R}(1+f)w_{0}\,\mathrm{d}q\,\mathrm{d}d}\end{align}

This effectively means that the annular area marked by ${r_{rim}}$ has a larger areal density $\left( {1 + f} \right)\sigma $ , which results in the uneven mass distribution.

Simulations have been carried out for different values of ${r_{rim}}$ and $f$ for shades with different areal density $\sigma $ . The numerical studies have revealed a few important things as discussed below.

• • Though the absolute ${o_B}$ values increased rapidly with $\alpha $ (increasing from ${o_B} \approx 264{\,\textrm{km}}$ at $\alpha = 0.0005$ to ${o_B} \approx 5990{\,\textrm{km}}$ at $\alpha = 0.01$ ), the fractional impact of rim mass enhancement factor $f = 0.05$ from a uniform to non-uniform mass distribution on ${o_B}$ shift is around $\approx 0.184{{\% }}$ over the range of $\alpha $ . For $f = 0.03$ , the fractional ${o_B}$ shift comes down to $\approx 0.12{{\% }}$ .

• • The rotational stability structure of the shade has been found to be very sensitive to the shift in the ${o_B}$ position (mass distribution). The evolution of eigenvalues for $\sigma = 20,000{\,\textrm{kg}}/{\textrm{k}}{{\textrm{m}}^2}$ for a range of rim mass enhancement values $f$ is depicted in Fig. 13. In all these cases, the enhanced mass is concentrated in the outermost $2{{\% }}$ of the shade. In the uniform-mass case $\left(\, {f = 0} \right)$ , the rotational mode is stable for all $\alpha $ above $1.4 \times {10^{ - 4}}$ . For nonzero $f$ , an instability appears within the otherwise-stable region, as can be observed in Fig 13b and c. A pair of eigenvalues has non-zero real parts over a finite range of $\alpha $ , producing an ‘unstable bubble’ in the stability map.

• • The influence of the areal density on the stability regime of the paraboloid shade for non-uniform mass distribution was also investigated through another set of simulations. A representative figure depicting the stability behaviour for four areal density values for $f = 0.03$ is illustrated in Fig. 14. It can be observed from the plots in Fig. 14 that the region of stability for rotational modes expands with an increase in areal density $\sigma $ . The instability bubble observed for $\sigma = 28,000{\,\textrm{kg}}/{\textrm{k}}{{\textrm{m}}^2}$ slowly diminishes and disappears at $\sigma = 32,500{\,\textrm{kg}}/{\textrm{k}}{{\textrm{m}}^2}$ .

Figure 14.

Real parts of eigenvalues for investigation of peripheral mass distribution on paraboloid shade stability, as the curvature parameter $\alpha $ is varied. Results of simulations for four different values of the areal density $\sigma $ have been shown for the rim mass enhancement factor $f = 0.03$ .

It can be inferred from the above set of simulation results that the mass distribution of the control hardware remains a critical design parameter for the architecture, and the control system mass budget must be specified exactly as part of any operational implementation. Compensation strategies can include optimal choice of $\alpha $ and $\sigma $ based on stability maps as well as use of counter-masses on the anti-Sun side of the shade close to the apex. This can be used for restoring the wider stability range of the uniform-density case.

4.3 Non-linear recovery for larger angular perturbations

As mentioned in Section 2.2, the stability analysis presented in this study is based on a rigid body model of the shade. Second, only small perturbations of the different Euler angles were assumed in Section 2.2, which was then followed by the derivation of the different force and torque equations. This has been carried out through a linearisation of the equations of motion about the equilibrium configuration, which is valid in the limit of small angular perturbations. However, to verify that the linear stability conclusions remain quantitatively valid for moderate perturbations of practical interest, a representative single degree-of-freedom rotational simulation was conducted and is presented here.

The investigated case is for the parameters $\sigma = 40,000{\,\textrm{kg}}/{\textrm{k}}{{\textrm{m}}^2}$ and $\alpha = 5 \times {10^{ - 4}}{\textrm{k}}{{\textrm{m}}^{ - 1}}$ . As depicted in Fig. 4b, stabilised rotational modes occur for this choice. The linearised rotational dynamics for this state can be said to be governed by an oscillator with a natural frequency of $\omega = 1.88 \times {10^{ - 4}}$ rad/s, which is the imaginary part of the corresponding eigenvalue. Physics wise, this corresponds to the rotational mode near the equilibrium, which can be modeled as follows:

(38) \begin{align}{\ddot{\theta}} + \omega^{2}\theta = 0\end{align}

Here, $\theta $ denotes the libration angle of the shade’s symmetry axis relative to the equilibrium (Sun-pointing) orientation, considered for a single transverse rotation axis. From the rotational symmetry, the single degree-of-freedom treatment considered for the simulations here is representative for both the pitch and roll modes.

A representative nonlinear extension of Equation (38) can be obtained by replacing the linear restoring torque with the form ${\omega ^2}{\textrm{sin}}\theta {\textrm{cos}}\theta $ , which reduces to the linear case at small angles. While the considered single degree-of-freedom model does not capture the full coupling present in the multi degree-of-freedom model of the original system, this analysis serves as a representative test of whether the period of oscillation is significantly altered by the nonlinear ${\textrm{sin}}\theta {\textrm{cos}}\theta $ dependence of the restoring torque.

As depicted in Fig. 15, for initial perturbations of ${5^ \circ }$ and ${10^ \circ }$ , the nonlinear oscillation period departs from the linear prediction by $0.25{{\% }}$ and $0.85{{\% }},$ respectively. It is observed that the trajectory remains within ${0.17^ \circ }$ and ${1.32^ \circ }$ of the linear solution over multiple oscillation periods. However, the deviations between the linear solutions and the non-linear oscillations grow with amplitude. This becomes clear in the plots depicted in Fig. 15 beyond ${20^ \circ }$ . The results confirm that the linearised stability conclusions are quantitatively reliable throughout a moderate perturbation range till ${15^ \circ }$ . A complete nonlinear analysis with a multi degree-of-freedom model, including translation–rotation coupling and pitch and roll coupling, is identified as an important direction for future investigation.

Figure 15.

Evolution of the recovery from different initial angular perturbations ranging from ${5^ \circ }$ to ${30^ \circ }$ . Dashed lines are used to depict the non-linear recovery trajectories and solid lines are used to depict the trajectories following the solution of the linearised eigenvalue problem.