Representation of light pressure resultant force and moment as a tensor series

TL;DR

This paper introduces a tensor series approach using Chebyshev polynomials to accurately model light pressure forces and moments on complex space structures, accounting for shadowing and reflections.

Contribution

It presents a novel tensor series method for calculating light pressure on complex geometries, including self-shadowing and reflection effects, with applications to space structures like solar sails.

Findings

Derived tensor-based expressions for light pressure forces and moments.

Validated the approach on flat, spherical, cylindrical bodies, and a space observatory.

Enabled improved analysis of space structure dynamics under radiation pressure.

Abstract

In this article we address a problem of determination of light pressure upon space structures of the complex geometric shape. For the surface element, we wrote a condition that this element can interact with light only from the front side, not from the back side. This condition To in the form of series of Chebyshev polynomials of the first kind. Chebyshev expansion lets us move to the series of tensors of increasing rank for the problem of determination of force and moment. We obtained expressions for the fiber method for determination of light pressure on space structures of complex geometry taking into account self--shadowing and reflections within the structure. We also give the expressions for tensor parametrization by the specularity coefficient in case of specular-diffuse reflection. For these expressions we calculated the principal moment and force upon two-sided flat solar sail, spherical and cylindrical bodies, and approximated light pressure upon a perspective space observatory Millimetron. The proposed expressions can be used in the ballistic analysis of solar sails and other space objects, which are significantly affected by the radiation pressure. Also, these results can be used to analyze the dynamics of movement of large-scale space structures around the center of gravity under the light pressure.

Equations

$$d\mathbf{F}=P(R)\left[-a_0\hat{\mathbf{n}}-a_1(\hat{\mathbf{n}}\cdot \hat{\mathbf{s}})\hat{\mathbf{s}}+a_2(\hat{\mathbf{n}}\cdot \hat{\mathbf{s}})\hat{\mathbf{n}}-2a_3(\hat{\mathbf{n}}\cdot \hat{\mathbf{s}}{)}^2\hat{\mathbf{n}}\right]dA,$$

$$a_0=\frac{\epsilon B\sigma T^4}{P(R)};$$

$$a_1=1-\rho s;$$

$$a_2=B\rho (1-s);$$

$$a_3=\rho s.$$

$$V(\hat{\mathbf{n}},\hat{\mathbf{s}})=\frac{\hat{\mathbf{n}}\cdot \hat{\mathbf{s}}-|\hat{\mathbf{n}}\cdot \hat{\mathbf{s}}|}{2}.$$

$$\displaystyle d\mathbf{F}=\frac{P(R)}{2}\bigg{(}-2a_{0}\hat{\mathbf{n}}-a_{1}(\hat{\mathbf{n}}\cdot\hat{\mathbf{s}}-|\hat{\mathbf{n}}\cdot\hat{\mathbf{s}}|)\hat{\mathbf{s}}+a_{2}(\hat{\mathbf{n}}\cdot\hat{\mathbf{s}}-|\hat{\mathbf{n}}\cdot\hat{\mathbf{s}}|)\hat{\mathbf{n}}$$

$$\displaystyle-a_{3}(\hat{\mathbf{n}}\cdot\hat{\mathbf{s}}-|\hat{\mathbf{n}}\cdot\hat{\mathbf{s}}|)^{2}\hat{\mathbf{n}}\bigg{)}dA,$$

$$\displaystyle d\mathbf{F}=\frac{P(R)}{2}\bigg{(}-2a_{0}\hat{\mathbf{n}}-a_{1}(\hat{\mathbf{n}}\cdot\hat{\mathbf{s}})\hat{\mathbf{s}}+a_{2}(\hat{\mathbf{n}}\cdot\hat{\mathbf{s}})\hat{\mathbf{n}}-2a_{3}(\hat{\mathbf{n}}\cdot\hat{\mathbf{s}})^{2}\hat{\mathbf{n}}+a_{1}|\hat{\mathbf{n}}\cdot\hat{\mathbf{s}}|\hat{\mathbf{s}}$$

$$\displaystyle-a_{2}|\hat{\mathbf{n}}\cdot\hat{\mathbf{s}}|\hat{\mathbf{n}}+2a_{3}(\hat{\mathbf{n}}\cdot\hat{\mathbf{s}})|\hat{\mathbf{n}}\cdot\hat{\mathbf{s}}|\hat{\mathbf{n}}\bigg{)}dA.$$

$$|\hat{\mathbf{n}}\cdot \hat{\mathbf{s}}|=\frac{2}{\pi }-\frac{4}{\pi }\sum _{n=1}^{\infty }\frac{(-1{)}^n{T}_{2n}(\hat{\mathbf{n}}\cdot \hat{\mathbf{s}})}{-1+4n^2}=$$

$$=-\frac{4}{\pi }\sum _{n=1}^{\infty }\sum _{k=0}^{n-1}\frac{(-1{)}^n(-1{)}^{k}n(2n-k-1)!}{(-1+4n^2)k!(2n-2k)!}{4}^{n-k}(\hat{\mathbf{n}}\cdot \hat{\mathbf{s}}{)}^{2(n-k)}.$$

$$|\hat{\mathbf{n}}\cdot \hat{\mathbf{s}}|=\sum _{m=1}^{\infty }{B}_m(\hat{\mathbf{n}}\cdot \hat{\mathbf{s}}{)}^{2m},$$

$$B_m=-\frac{(-1{)}^m{4}^{m+1}}{\pi (2m)!}\sum _{n=m}^{\infty }\frac{n(n+m-1)!}{(-1+4n^2)(n-m)!}.$$

$$N_{\max B}=\lfloor \frac{N_{\max }-1}{2}\rfloor ,$$

$$B_m\approx -\frac{(-1{)}^m{4}^{m+1}}{\pi (2m)!}\sum _{n=m}^{\lfloor \frac{N_{\max }-1}{2}\rfloor }\frac{n(n+m-1)!}{(-1+4n^2)(n-m)!}.$$

$$\displaystyle d\mathbf{F}=\frac{P(R)}{2}\bigg{(}-2a_{0}\hat{\mathbf{n}}-a_{1}(\hat{\mathbf{n}}\cdot\hat{\mathbf{s}})\hat{\mathbf{s}}+a_{2}(\hat{\mathbf{n}}\cdot\hat{\mathbf{s}})\hat{\mathbf{n}}-2a_{3}(\hat{\mathbf{n}}\cdot\hat{\mathbf{s}})^{2}\hat{\mathbf{n}}+$$

$$\displaystyle a_{1}\sum\limits_{m=1}^{{N_{\max}}}B_{m}(\hat{\mathbf{n}}\cdot\hat{\mathbf{s}})^{2m}\hat{\mathbf{s}}-a_{2}\sum\limits_{m=1}^{{N_{\max}}}B_{m}(\hat{\mathbf{n}}\cdot\hat{\mathbf{s}})^{2m}\hat{\mathbf{n}}+2a_{3}\sum\limits_{m=1}^{{N_{\max}}}B_{m}(\hat{\mathbf{n}}\cdot\hat{\mathbf{s}})^{2m+1}\hat{\mathbf{n}}\bigg{)}dA.$$

$$\underset{p}{\underbrace{(\hat{\mathbf{n}}\cdot \hat{\mathbf{s}})\dots (\hat{\mathbf{n}}\cdot \hat{\mathbf{s}})}}\hat{\mathbf{n}}=(\underset{p+1}{\underbrace{\hat{\mathbf{n}}\otimes \hat{\mathbf{n}}\otimes \dots \otimes \hat{\mathbf{n}}}})\cdot \underset{p}{\underbrace{\hat{\mathbf{s}}\cdot \dots \cdot \hat{\mathbf{s}}}}={\mathcal{J}}_A^{p+1}\cdot \underset{p}{\underbrace{\hat{\mathbf{s}}\cdot \dots \cdot \hat{\mathbf{s}}}};$$

$$\underset{p}{\underbrace{(\hat{\mathbf{n}}\cdot \hat{\mathbf{s}})\dots (\hat{\mathbf{n}}\cdot \hat{\mathbf{s}})}}\hat{\mathbf{s}}=(\underset{p}{\underbrace{\hat{\mathbf{n}}\otimes \hat{\mathbf{n}}\otimes \dots \otimes \hat{\mathbf{n}}}}\otimes {\mathcal{E}}^2)\cdot \underset{p+1}{\underbrace{\hat{\mathbf{s}}\cdot \dots \cdot \hat{\mathbf{s}}}}={\mathcal{J}}_B^{p+2}\cdot \underset{p+1}{\underbrace{\hat{\mathbf{s}}\cdot \dots \cdot \hat{\mathbf{s}}}},$$

$$\displaystyle d\mathbf{F}=\frac{P(R)}{2}\bigg{(}-2a_{0}\hat{\mathbf{n}}-a_{1}\mathcal{J}_{B}^{3}\cdot\hat{\mathbf{s}}\cdot\hat{\mathbf{s}}+a_{2}\mathcal{J}_{A}^{2}\cdot\hat{\mathbf{s}}-2a_{3}\mathcal{J}_{A}^{3}\cdot\hat{\mathbf{s}}\cdot\hat{\mathbf{s}}+$$

$$a_1\sum _{m=1}^{N_{\max }}{B}_m{\mathcal{J}}_B^{2m+2}\cdot \underset{2m+1}{\underbrace{\hat{\mathbf{s}}\cdot \dots \cdot \hat{\mathbf{s}}}}-a_2\sum _{m=1}^{N_{\max }}{B}_m{\mathcal{J}}_A^{2m+1}\cdot \underset{2m}{\underbrace{\hat{\mathbf{s}}\cdot \dots \cdot \hat{\mathbf{s}}}}+$$

$$\displaystyle 2a_{3}\sum\limits_{m=1}^{{N_{\max}}}B_{m}\mathcal{J}_{A}^{2m+2}\cdot\underbrace{\hat{\mathbf{s}}\cdot\ldots\cdot\hat{\mathbf{s}}}_{2m+1}\bigg{)}dA.$$

$$d\mathbf{F}=P(R)\bigg{(}\mathcal{J}^{1}+\sum\limits_{n=2}^{{N_{\max}}}\mathcal{J}^{n}\cdot\underbrace{\hat{\mathbf{s}}\cdot\ldots\cdot\hat{\mathbf{s}}}_{n-1}\bigg{)}dA,$$

$${\mathcal{J}}^1=-a_0\hat{\mathbf{n}};$$

$${\mathcal{J}}^2=\frac{1}{2}{a}_2{\mathcal{J}}_A^2;$$

$$\displaystyle\mathcal{J}^{3}=\frac{1}{2}\bigg{(}-a_{1}\mathcal{J}_{B}^{3}-2a_{3}\mathcal{J}_{A}^{3}-B_{1}a_{2}\mathcal{J}_{A}^{3}\bigg{)};$$

$$\displaystyle\mathcal{J}^{n}=\frac{1}{2}\bigg{(}-B_{\frac{n-1}{2}}\frac{1-(-1)^{n}}{2}a_{2}\mathcal{J}_{A}^{n}+B_{\frac{n-2}{2}}\frac{1+(-1)^{n}}{2}\left(a_{1}\mathcal{J}_{B}^{n}+2a_{3}\mathcal{J}_{A}^{n}\right)\bigg{)},\ n>3;$$

$${\mathcal{J}}_A^n=\underset{n}{\underbrace{\hat{\mathbf{n}}\otimes \dots \otimes \hat{\mathbf{n}}}};$$

$${\mathcal{J}}_B^n=\underset{n-2}{\underbrace{\hat{\mathbf{n}}\otimes \dots \otimes \hat{\mathbf{n}}}}\otimes {\mathcal{E}}^2,$$

$$\mathbf{r}\times \mathbf{a}={\mathcal{R}}^2\cdot \mathbf{a}.$$

$${\mathcal{R}}^2=\left(\begin{array}{ccc}0& -r_3& r_2\\ r_3& 0& -r_1\\ -r_2& r_1& 0\end{array}\right).$$

$$\underset{p}{\underbrace{(\hat{\mathbf{n}}\cdot \hat{\mathbf{s}})\dots (\hat{\mathbf{n}}\cdot \hat{\mathbf{s}})}}({\mathcal{R}}^2\cdot \hat{\mathbf{n}})=(\underset{p}{\underbrace{\hat{\mathbf{n}}\otimes \hat{\mathbf{n}}\otimes \dots \otimes \hat{\mathbf{n}}}}\otimes {\mathcal{R}}^2\cdot \hat{\mathbf{n}})\cdot \underset{p}{\underbrace{\hat{\mathbf{s}}\cdot \dots \cdot \hat{\mathbf{s}}}}={\mathcal{L}}_A^{p+1}\cdot \underset{p}{\underbrace{\hat{\mathbf{s}}\cdot \dots \cdot \hat{\mathbf{s}}}};$$

$$\underset{p}{\underbrace{(\hat{\mathbf{n}}\cdot \hat{\mathbf{s}})\dots (\hat{\mathbf{n}}\cdot \hat{\mathbf{s}})}}({\mathcal{R}}^2\cdot \hat{\mathbf{s}})=(\underset{p}{\underbrace{\hat{\mathbf{n}}\otimes \hat{\mathbf{n}}\otimes \dots \otimes \hat{\mathbf{n}}}}\otimes {\mathcal{R}}^2)\cdot \underset{p+1}{\underbrace{\hat{\mathbf{s}}\cdot \dots \cdot \hat{\mathbf{s}}}}={\mathcal{L}}_B^{p+2}\cdot \underset{p+1}{\underbrace{\hat{\mathbf{s}}\cdot \dots \cdot \hat{\mathbf{s}}}}.$$

$$\displaystyle d\mathbf{M}=\mathbf{r}\times d\mathbf{F}=\mathcal{R}^{2}\cdot d\mathbf{F}=\frac{P(R)}{2}\bigg{(}-2a_{0}\mathcal{R}^{2}\cdot\hat{\mathbf{n}}-a_{1}\mathcal{L}_{B}^{3}\cdot\hat{\mathbf{s}}\cdot\hat{\mathbf{s}}+a_{2}\mathcal{L}_{A}^{2}\cdot\hat{\mathbf{s}}-$$

$$2a_3{\mathcal{L}}_A^3\cdot \hat{\mathbf{s}}\cdot \hat{\mathbf{s}}+a_1\sum _{m=1}^{N_{\max }}{B}_m{\mathcal{L}}_B^{2m+2}\cdot \underset{2m+1}{\underbrace{\hat{\mathbf{s}}\cdot \dots \cdot \hat{\mathbf{s}}}}-a_2\sum _{m=1}^{N_{\max }}{B}_m{\mathcal{L}}_A^{2m+1}\cdot \underset{2m}{\underbrace{\hat{\mathbf{s}}\cdot \dots \cdot \hat{\mathbf{s}}}}+$$

$$\displaystyle 2a_{3}\sum\limits_{m=1}^{{N_{\max}}}B_{m}\mathcal{L}_{A}^{2m+2}\cdot\underbrace{\hat{\mathbf{s}}\cdot\ldots\cdot\hat{\mathbf{s}}}_{2m+1}\bigg{)}dA,$$

$${\mathcal{L}}_A^n=\underset{n-1}{\underbrace{\hat{\mathbf{n}}\otimes \dots \otimes \hat{\mathbf{n}}}}\otimes {\mathcal{R}}^2\cdot \hat{\mathbf{n}};$$

$${\mathcal{L}}_B^n=\underset{n-2}{\underbrace{\hat{\mathbf{n}}\otimes \dots \otimes \hat{\mathbf{n}}}}\otimes {\mathcal{R}}^2.$$

$$d\mathbf{M}=P(R)\bigg{(}\mathcal{L}^{1}+\sum\limits_{n=2}^{{N_{\max}}}\mathcal{L}^{n}\cdot\underbrace{\hat{\mathbf{s}}\cdot\ldots\cdot\hat{\mathbf{s}}}_{n-1}\bigg{)}dA,$$

Full text

1. Introduction

This work111This work was prepared during development of BMSTU-Sail Space Experiment Rachkin et al (2011); bms (2014) is related to the calculation of the light pressure force on space structures. The phenomenon of light pressure with radiations by J. Maxwell based on his theory of electromagnetism Maxwell (1873). The light radiation pressure was found by P. Lebedev Lebedew (1901) in the series of experiments. The first one who suggested using the force of light pressure to fly in space was F. Tsander in the mid 20-s of the 20 century Tsander (1969). Regarding solar sails, there is a sufficiently complete overview of the current state of solar sails development in the book McInnes (2004). Up to the present moment several experiments regarding solar sail technology were conducted in space, including Znamya-1 and Znamya-2 Raykunov et al (2009), Nanosail-D2 Alhorn et al (2011); Johnson et al (2011), Lightsail-1 Ridenoure et al (2015), IKAROS Tsuda et al (2011); Kawaguchi (2014).

Solar radiation pressure (SRP) is affecting celestial bodies providing various effects Burns et al (1979, 2014). One of these effects is the phenomenon called Yarkovsky effect Neiman et al (1965); Öpik (1951); Beekman (2006). It is originated from the anisotropy of heat radiation emission from heated surface of celestial body Radzievskii (1952) and provides the additional moment which causes acceleration or deceleration of spinning Paddack (1969); Rubincam (2000) and also affecting the orbital parameters Rubincam (1995). For small bodies (e.g. asteroids) it can increase the spinning velocity up to the disintegration by centrifugal forces Radzievskii (1954). There are several good analytical models of SRP on celestial bodies, including Vokrouhlický’s model Vokrouhlický and Farinella (1998); Vokrouhlický (1999); Vokrouhlický and Čapek (2002), other models by Katasev and Kulikova (1981); Hartmann et al (1999); Bottke et al (2006) etc.

It is necessary to study the effects of SRP for the precise orbit propagation. There are several studies related to research of solar radiation effects on GPS satellites Bar-Sever and Russ (1997); Bar-Sever and Kuang (2004); Rodriguez-Solano et al (2012); Fliegel et al (1992); Fliegel and Gallini (1996); Springer et al (1999) It is important to analyze the angular movement around the center of inertia for space structures in distant space far from Earth’s atmosphere in low gravity gradient conditions because moment from light radiation can affect the attitude of spacecrafts Kinzel (2005).

It is possible to create a more exact model of SRP by numerical integration of SRP equations around the surface represented as a mesh of simple elements Ziebart (2004).

During the development of the theory of light pressure acting on space objects some authors got the result according to which in the equations for the resultant force and moment of light pressure on arbitrary body it is possible to separate the description of structure’s surface from its orientation Rios-Reyes and Scheeres (2004, 2005); Rios-Reyes (2006); Rios-Reyes and Scheeres (2007); Scheeres (2007); McMahon and Scheeres (2010, 2014, 2015). These authors developed a so-termed tensor approach for the description of light pressure. Tensor approach is also being developed in the works Jing et al (2012, 2014). This work was done using the tensor approach for the description of the light pressure force, and it is a continuation of the author’s previous research Nerovnyi and Zimin (2014); Zimin and Nerovnyy (2015); Zimin and Nerovnyi (2016).

2. Light pressure on an infinitesimal surface element

As it was shown in work Rios-Reyes and Scheeres (2004) it is possible to write analytical expressions for the resultant force and moment of light pressure forces on the spacecraft surface under certain assumptions. Now we will derive the expression for the resultant vector and resultant moment of an arbitrary convex structure.

Let us write the expression for the light pressure force on the infinitesimal surface element in general form:

[TABLE]

where generalized optical parameters are given as the following:

[TABLE]

In this case, the following symbols are used:

• •

$P(R)=q_{0}/c$ — light pressure on a flat diffuse area normal to the incident radiation at a distance $R$ from the Sun, which is assumed to be a point source, and the light flux at the distance $R$ is $q_{0}$, where $c$ is a speed of light in vacuum;

• •

$\hat{\mathbf{n}}$ — local unit vector of normal to the surface of space structure;

• •

$\hat{\mathbf{s}}$ — unit vector defining the direction from the light source to the surface in the spacecraft coordinate system;

• •

$\varepsilon$ — emissivity of the surface element;

• •

$\sigma$ — Boltzmann constant;

• •

$T$ — temperature of surface element;

• •

$B$ — coefficient equal to $2/3$ for diffuse scattering surface (makes sense only for the materials without anisotropy of optical parameters depending on direction along the surface), it is often called as Lambertian coefficient;

• •

$\rho$ — total reflectivity of the surface material (spectral reflectivity integrated by wavelength from 0 to $\infty$);

• •

$s$ — specularity coefficient showing which part of the radiation is reflected glossy, and what – diffusely.

In the other studies (e.g. Rios-Reyes and Scheeres (2004)) there was an implicit assumption that the lightning configuration does not change with time i.e. there is no change in the set of lighted and shadowed surfaces of spacecraft. In the subsequent derivations for the convex structure we assume that the structure can take an arbitrary orientation relative to the incident solar radiation. Taking into account of this condition is similar to replacement in (1) of expression $\hat{\mathbf{n}}\cdot\hat{\mathbf{s}}$ by the following relationship $V(\hat{\mathbf{n}},\hat{\mathbf{s}})$:

[TABLE]

The above statement can be interpreted as a visibility function in the works Scheeres (2007). (McMahon and Scheeres, 2010, equation (15)).

In point of fact, if $\hat{\mathbf{n}}\cdot\hat{\mathbf{s}}>0$, it turns out that the area is illuminated from the internal side which cannot be. Let us write an expression (1) with this new relationship:

[TABLE]

and after transformations,

[TABLE]

Since the range of the $\hat{\mathbf{n}}\cdot\hat{\mathbf{s}}$ is $[-1;1]$, then $|\hat{\mathbf{n}}\cdot\hat{\mathbf{s}}|$ may be expanded in the well-known series of Chebyshev polynomials of the first kind:

[TABLE]

The main purpose of Chebyshev expansion is the representation of non-analytical visibility function as a power series of $\hat{\mathbf{n}}\cdot\hat{\mathbf{s}}$. After this, we can utilize the method originated from works by D.J. Scheeres, L. Rios-Reyes, and others Rios-Reyes and Scheeres (2004) by transforming of multiple scalar and vector products into the assembly of dyadic for $\hat{\mathbf{n}}$ and several scalar products of $\hat{\mathbf{s}}$. The original method of D.J. Scheeres and others provides an analytical separation of attitude from the description of the surface, so one can calculate some tensor parameters once in the beginning and then calculate the resultant force and moment using simple scalar products by orientation vector which allows using this method in the restricted calculation environments (e.g. in spacecraft’s onboard computer) or as a part of dynamics model of spacecraft.

After transformations we can write the series (3) as follows:

[TABLE]

where

[TABLE]

The series in the expression (5) diverge, however for any finite number of terms of the series in the final summation (4) we have a convergent sequence. Hereafter for the infinite sums we will consider that they either constrained by the maximum number of terms $N_{\max}$, or renormalization has been made if appropriate. Analysis of possible renormalization of these series is beyond the scope of this article. If we constrain the number of terms in the series (4) by $N_{\max}$, the upper bound for $B_{m}$ series (5) will be as follows:

[TABLE]

where $\lfloor x\rfloor$ is a floor function of real $x$. We can write the approximate equation:

[TABLE]

By substituting expression (4) in (2) we can obtain an expression for the infinitesimal force of light pressure as a series form:

[TABLE]

3. The resultant vector and moment of light pressure forces

Using the approach from work Rios-Reyes and Scheeres (2004) we can write the following relations:

[TABLE]

where $\mathcal{J}_{A}^{p+1}$ и $\mathcal{J}_{B}^{p+2}$ – some tensors in Euclidean space with a rank of $p+1$ и $p+2$, respectively.

Let us write (7) in the new tensor notation:

[TABLE]

By grouping members of infinite sums in (8) with each other and by introducing additional notation we can obtain an expression for the radiation pressure force on unit area which receives only the radiation that falls on the front side:

[TABLE]

where

[TABLE]

where $\mathcal{E}^{2}$ — unit tensor of second rank.

Let introduce the tensor $\mathcal{R}^{2}$ as follows:

[TABLE]

Obviously, the tensor $\mathcal{R}^{2}$ from (10) is skew-symmetric and can be represented in the following matrix form:

[TABLE]

For the infinitisemal moment we will use the following equations:

[TABLE]

By writing (7) in the new tensor notation we can get:

[TABLE]

where

[TABLE]

After transformations similar to those we did in the derivation of the expression of the force, we can obtain an expression for the moment from infinitesimal force of light pressure:

[TABLE]

where

[TABLE]

By integrating expression (9) and (14) over the entire surface let us introduce the following values:

[TABLE]

where $n>0$.

Finally we can get the expressions for the resultant force and moment of light pressure on a body with convex shape:

[TABLE]

4. Light pressure on the structure of the arbitrary shape

Considering the problem of the calculation of the radiation pressure force the structure with complex geometry taking into account self-shadowing and reflections in structure, we should note the fact that this problem is related to the problem of radiative heat transfer in a complex structure, which has been shown multiple times by different authors, it cannot be solved in unified general form for the structure of arbitrary configuration. Howell et al (2015). Nevertheless, there are some approximate methods to calculate the radiative heat transfer in complex structures from the view factor algebra to the Monte–Carlo methods (ibid.).

This problem is also close to the problems faced in computer graphics in the analysis of global illumination Liu et al (2004); Tsai and Shih (2006); Kristensen et al (2005); Mäki-Patola (2003); Sloan et al (2002). In these works authors used a method of calculation of global lightning with precomputed bidirectional radiation distribution function (BRDF) for each element of the surface. After the BRDF calculation, one can dynamically change the ambient light luminance map and calculate the brightness of each surface element. Moreover, if the construction of BRDF is a complex computational problem, the process of restoring of luminance of surface of the object with existing BRDF is less resource consuming.

Given the above, we present the main assumptions used for the development of the calculation method for determination of the light pressure on the structure of complex geometry.

Assumption 1**.**

We assume that the expressions of infinitesimal force and moment on the body with the complex shape can also be written in a form similar to (9) and (14) respectively.

This statement will be briefly called as the assumption of similarity of notation.

Assumption 2**.**

Since some surface elements can be illuminated not by direct sunlight but by light reflected and re-emitted by structure, we will assume that the vector of direction from the sun $\mathbf{s}$ is not a constant unit vector and it is dependent on the position of infinitesimal area $dA$.

This statement will be briefly called as an assumption of variable light flux. Vector $\mathbf{s}$ will be called as a vector of a local light load.

Assumption 3**.**

BRDF will be approximated by a second-rank tensor $\mathbf{s}=\mathcal{D}^{2}\cdot\hat{\mathbf{s}}$, where the unit vector $\hat{\mathbf{s}}$ is constant for the entire surface.

This proposal will be called as a distribution assumption.

Assumption 4**.**

During calculation of BRDF we assume that the vector $\mathbf{s}$ is zero if there are no such vector $\hat{\mathbf{s}}$ for which the area $dA$ is been illuminated.

Assumption 5**.**

When taking into account the light pressure from thermal radiation from structure we assume that the magnitude of force is also associated with the BRDF, $\mathbf{f}\propto\mathcal{D}^{2}\cdot\hat{\mathbf{n}}$.

In the distribution assumption it is possible to introduce the high–order approximation, however, we decided to simplify the problem by limiting the approximation by the tensor of a second rank. Thus, the need to use higher-order approximation remains still open. It should be noted that the expressions (9) and (14) are already contain features, according to which the area element $dA$ does not receive radiation incident from the internal side.

By analogy with section 3 we can transform the vector expressions in (9) and (14) assuming approximation

[TABLE]

We can obtain the following expressions:

[TABLE]

Since the tensors $\tilde{\mathcal{J}}_{A}$, $\tilde{\mathcal{J}}_{B}$, $\tilde{\mathcal{L}}_{A}$, $\tilde{\mathcal{L}}_{B}$ are introduced as well as tensors $\mathcal{J}_{A}$, $\mathcal{J}_{B}$, $\mathcal{L}_{A}$, $\mathcal{L}_{B}$ respectively, the final expressions for the infinitesimal force and moment for the surface area element $dA$ will look like the follows:

[TABLE]

where

[TABLE]

By integrating (20) and (21) across the surface $A$ we can obtain the final expression for the resultant force and moment of light pressure on the spacecraft with complex geometric shape:

[TABLE]

where

[TABLE]

where $n>0$.

Expressions (22) — (38) should be supplemented by the procedure of calculation of tensor approximation $\mathcal{D}^{2}$:

[TABLE]

In case of convex shape $\mathcal{D}^{2}=\mathcal{E}^{2}$. In this work, we will not discuss the calculation techniques for determination of BRDF, some of them can be derived from works by Kristensen et al (2005) and Tsai and Shih (2006) by the transition from spherical harmonics representation of BRDF to the tensor representation. The above definitions are aimed to create the equations of resultant force and moment in which the parameters of the structure are analytically separated from its attitude towards the sun, eq. (35) and (36). For the given equations we will develop the approximation method of SRP in the section 5 and show the numerical example results in the section 8.

Expressions (22) — (39) give us a full statement of the approximate method for determination of resultant vector and moment of light radiation pressure upon a structure with complex geometry taking into account self-shadowing, multipath reflection, a variability of the optical parameters and temperature over the surface of the structure. The reduced pressure from reflected light is incorporated into the BRDF as far as tensor $\mathcal{D}^{2}$ can generate the non-unit vector $\mathbf{s}$ of light load. The self-shadowing is considering in two ways: in the Chebyshev expansion of visibility function and BRDF $\mathcal{D}^{2}$. The surface properties can be specular, diffuse, specular diffuse, as well as other rotationally symmetrical scattering law other than the Lambert law can be used.

Note that due to the physical meaning of the BRDF it cannot create infinite local light load vector from the finite sun vector load as this would mean the infinite concentration of energy in the structure which is impossible to achieve. The components of the tensor $D^{2}_{ij}$ always have a finite value for the whole area $A$. On the other hand the approximation tensor $\mathcal{D}^{2}$ is included in the expression for the $\tilde{\mathcal{J}}_{A}$, $\tilde{\mathcal{J}}_{B}$, $\tilde{\mathcal{L}}_{A}$, $\tilde{\mathcal{L}}_{B}$ in a variable number of times, which depends on the rank of approximation. Therefore, the question of convergence of the series in terms (35) and (36) remains still open. the

5. The method of approximating of the light pressure by tensor series

The resulting expressions for the resultant force and moment in the series form allows us to offer a method of approximation of light radiation pressure.

Let us assume that in some way (e.g. by Monte–Carlo ray tracing) we obtained the values of normalized force $\mathbf{f}^{(i)}=\mathbf{F}^{(i)}/P(R)$ and normalized moment $\mathbf{m}^{(i)}=\mathbf{M}^{(i)}/P(R)$ for the given set of different orientations $\hat{\mathbf{s}}^{(i)}$ in the amount of $N$. Superscript $(i)$ indicates the number of orientation case for which $\mathbf{F}$ and $\mathbf{M}$ were calculated by the different method (e.g. by Monte–Carlo), from $1$ to $N$. We need to get approximated values of tensor components in the expansion for $\mathbf{F}$ and $\mathbf{M}$. To do this, we define the vectors for the unknown components of the tensors $\tilde{\mathcal{I}}^{n}$ и $\tilde{\mathcal{K}}^{n}$ (hereinafter we indicate vector and matrix dimensions in brackets):

[TABLE]

where $M$ – number of terms in the expansion of force and moment (for simplicity we assume that this number is the same for the calculation of force and moment). In the above two equations subscript is a vector component, superscript is the order of series expansion.

Let us define the vectors of free terms:

[TABLE]

The matrix of coefficients will be as the follows:

[TABLE]

For vectors $\mathbf{f}$ and $\mathbf{m}$ let us introduce the following overdefined system of linear equations:

[TABLE]

For $\mathbf{j}$ and $\mathbf{k}$ we can find approximations $\tilde{\mathbf{j}}$ and $\tilde{\mathbf{k}}$ by the least squares method:

[TABLE]

where + is a pseudo–inverse operator for square matrix.

After calculation of approximated values of $\mathbf{f}$ and $\mathbf{m}$ we can reconstruct the approximated tensors $\tilde{\mathcal{I}}^{n}$ and $\tilde{\mathcal{K}}^{n}$ using (42) and (43) respectively. After this it is possible to determine the approximated value of the resultant force and moment for the arbitrary orientation $\hat{\mathbf{s}}$ of light source around structure using the expressions (35) and (36) respectively.

6. Parametrization of the resultant vector and moment in case of specular-diffuse reflection

For the ballistic analysis, it is often necessary to study the impact of the change of parameters. That is necessary to make a parametrization for our problem. Parametrization problem can be viewed in two ways: as a problem of geometric parameterization and as a problem of parametrization of the optical parameters. In this work, we will make parametrization by specularity coefficient $s$.

Suppose we obtain an approximation for the tensors describing the force and moment of light pressure in pure diffuse and pure specular cases. For the tensors (37) of the first rank, we write the expression for diffuse and specular cases, taking the temperature distribute on the surface constant:

[TABLE]

For the tensors (37) of second rank we will get:

[TABLE]

therefore, assuming that the tensors $\tilde{\mathcal{J}}_{A}^{2}$ are constant at different $s$, we obtain

[TABLE]

For the tensors (37) of third rank let us consider diffuse and specular cases:

[TABLE]

then approximating expression for the third rank will be as follows:

[TABLE]

In the general case of arbitrary rank we will get:

[TABLE]

where we get

[TABLE]

Similarly, we obtain expressions for the tensors (38):

[TABLE]

the resultant force and moment will be determined by the expressions (35) and (36).

The expressions (49) — (55) represent the approximation of tensors $\tilde{\mathcal{I}}$ and $\tilde{\mathcal{K}}$ for specular-diffuse surface when we know the components of tensors $\tilde{\mathcal{I}}$ and $\tilde{\mathcal{K}}$ in two cases: where the reflection is completely diffuse ($s=0$) and completely specular ($s=1$).

7. Analytical examples

In the analytical examples below the light source orientation vector $\hat{\mathbf{s}}$ is defined by two angles $\alpha$ and $\beta$ as follows:

• •

$\alpha\in[0,2\pi]$ — angle between unit vector $\hat{\mathbf{e}}_{1}$ of axis $Ox_{1}$ and projection of vector $\hat{\mathbf{s}}$ on the plane $Ox_{1}x_{3}$;

• •

$\beta\in[-\frac{\pi}{2};\frac{\pi}{2}]$ — angle between plane $Ox_{1}x_{3}$ and vector $\hat{\mathbf{s}}$.

The components of vector $\hat{\mathbf{s}}$ can be written as follows:

[TABLE]

7.1. Two-sided specular-diffuse flat solar sail

Let us consider a flat solar sail with specular–diffuse surface properties which are different on the opposite sides. The coordinate system origin is placed in the geometric center of the solar sail. The front side (side 1) has reflection coefficient $\rho_{1}$ and specularity parameter $s_{1}$. The other side (side 2) has reflectivity $\rho_{2}$ and specularity $s_{2}$. The area of the solar sail is $A$.

The components of tensors $\mathcal{I}$ and $\mathcal{K}$ can be calculated by the following equations:

[TABLE]

By setting the maximum number of terms in the series as $N_{\max}=6$ we can obtain the following analytical expressions for the components of the principal force and moment of light pressure:

[TABLE]

To check the correctness of these statements we made a comparison with the known solutions for specular and diffuse solar sails Forward (1989); McInnes (2004); Polyakhova (2011).

At the first, we selected a two-sided solar sail, one side of which is completely specular, the other is absorbing, i.e. $\rho_{1}=1,\ \rho_{2}=0,\ s_{1}=s_{2}=1$. Solar sail area was set to 1. We calculated the non-zero components of the resultant vector of light pressure on a solar sail, namely projections on axises $Ox_{1}=Ox$ and $Ox_{3}=Oz$, depending on the angle of inclination $\alpha$ considering $\beta=0$. The calculation results are shown in Fig. 1. As can be seen from the figure, there is a fairly good convergence of the approximate solution to the exact solution. On this and subsequent figures (fig. 1, 2 and 3) the range of angles $\alpha\in[0;\pi]$ corresponds to the lightning of side 2, and range $\alpha\in[\pi;2\pi]$ corresponds to the lightning of side 1, and for $\alpha=\pi/2$ and $\alpha=3\pi/2$ the light strikes perpendicular to the surface.

Then we selected the two-sided solar sail, one side of which is diffuse, $B=2/3$, the other – absorbing, i.e. $\rho_{1}=1$, $\rho_{2}=0$, $s_{1}=s_{2}=0$. Solar sail area was set to 1. We calculated the non-zero components of the resultant vector of light pressure on a solar sail, namely projections on axises $Ox_{1}=Ox$ and $Ox_{3}=Oz$, depending on the angle of inclination $\alpha$ considering $\beta=0$. The calculation results are shown in Fig. 2. As can be seen from the figure, there is a fairly accurate convergence of the approximate solutions to the exact solution.

For the solar sail with unit area with parameters $\rho_{1}=1,\ \rho_{2}=1,\ s_{1}=1,\ s_{2}=0$ we calculated the projections of principal force and moment on the axises $Ox_{1}=Ox$ and $Ox_{3}=Oz$ considering $\beta=0$ ( fig. 3). At the indicated chart values of ​​force are divided by $P(R)$.

7.2. Specular–diffuse sphere

Let us consider a sphere of radius $R_{0}$ with a homogeneous specular-diffusive surface, the reflection coefficient of which is equal to $\rho$ and the degree of specular reflection is $s$.

The expressions for the components of tensors $\mathcal{I}$ and $\mathcal{K}$ have the following form:

[TABLE]

We obtained the following analytical expressions considering $N_{\max}=6$:

[TABLE]

For $\rho=1$ and $s=1$ we get:

[TABLE]

7.3. Specular–diffuse cylinder

Let us consider the cylinder with following parameters:

• •

$\rho_{0}$ — reflectance of the envelope;

• •

$\rho_{1}$ — reflectance of the butt surface $+x_{3}$;

• •

$\rho_{2}$ — reflectance of the butt surface $-x_{3}$;

• •

$s_{0}$ — specularity coefficient of the envelope;

• •

$s_{1}$ — specularity coefficient of the butt surface $+x_{3}$;

• •

$s_{2}$ — specularity coefficient of the butt surface $-x_{3}$;

• •

$R_{1}$ — radius of the cylinder;

• •

$H$ — height of the cylinder.

The expressions for the components of tensors $\mathcal{I}$ and $\mathcal{K}$ have the following form:

[TABLE]

For the number of terms of the series $N_{\max}=6$ we obtained analytical dependences for the principal force and moment of light radiation pressure:

[TABLE]

On the figures 4 and 5 there are graphs of components of principal force and moment of light radiation pressure upon a cylinder with the following parameters: $\rho_{1}=\rho_{2}=0$, $\rho_{0}=1$, $s_{1}=s_{2}=0$, $s_{0}=1$, the height and radius are equal to 1m, depending on angle $\alpha$ considering $\beta=0$. In these figures for the angles $\alpha=\pi/2$ and $\alpha=3\pi/2$ light strikes perpendicularly to the one of the ends of cylinder.

8. Numerical example of $\mathbf{F}$ and $\mathbf{M}$ approximation

8.1. Model spacecraft

As a numerical example we selected the problem of approximation of resultant force and moment of light radiation pressure upon a perspective space observatory “Millimetron” (“Spectr-M”, Fig. 6)). The purpose of the Millimetron is to study various astronomical objects in the Universe at the wavelength range 60 $\mu$m — 20 mm with an unprecedented sensitivity in the single-dish observation mode and an extremely high spatial resolution as an element of a ground-space very long baseline interferometry (VLBI) system Kardashev et al (2000); Wild et al (2008).

To achieve that, the mission will be equipped with a 10-m diameter deployable primary mirror with the high surface accuracy of the reflective surface (RMS$\leq 10\mu m$), and even more large-sized sun shields for passive cooling. As well as high requirements to accuracy and stability of attitude control system are defined (1” and 0.2” respectively). Hence, а precise estimation of light pressure, as the main factor affecting the movement of the center of mass and around the center of mass, is so important.

The antenna and cryogenic instruments should be cooled down to 4.5 K by a combination of passive cooling with radiation shields and active cooling with mechanical coolers. These requirements led to a bunch of engineering challenges, connecting with the development of a reliable deployment system, light weight antenna structures with high thermo-elastic distortion stability and high accuracy of the reflective surface.

The deployment concept is based on the previously implemented structure in the Radioastron project Kardashev (1997). To achieve the required high surface accuracy of the primary mirror after the deployment and compensate reflecting surface distortions under the action of space environment factors, each petal is composed of a supporting spatial frame and three independent reflecting panels. The panels are mounted on the supporting frame via precise cryogenic actuators to adjust their reflecting surface. An active surface control system is used to control all reflecting panels of the mirror by using wave-front sensing. Additionally, a high modulus carbon fiber reinforced plastic (CFRP), providing a lightweight structure with high thermal stability, has been chosen for the material of the primary mirror.

The Millimetron space mission is an approved project as part of Russian Space Program with the provisional launch date after 2025. It is being led by the Astro Space Center of the P.N. Lebedev Physical Institute of the Russian Academy of Science in cooperation with numerous Russian and international organizations. The main industrial partner is the Reshetnev Information Satellite Systems Company.

The geometric model of the observatory is shown in Fig. 7. It should be noted that configuration of solar panels is still undecided and thus we analyze the linear design of solar panels in the model.

8.2. Description of Monte–Carlo raytracing

As a raytracing software, we used the Tracer program which was developed in Bauman MSTU Leonov (2012).

This software is used for analysis of radiation heat transfer in complex space structures. This software can use both specular, diffuse and specular–diffuse model of surface parameters.

To calculate the light pressure force we modified the modules which responsible for the processing of absorption of the light ray and its re-radiation (diffuse or specular or combined). Geometrical model of the structure is represented as a set of triangular elements. Representation of surface by simple elements can also be useful for calculation of the matrix of inertia of structure Dobrovolskis (1996) in the case of analysis of its angular motion. See Simonelli et al (1993) for the reconstruction of asteroid’s shape using a mesh of simple elements.

In the case of ray absorption the program records the energy of incident beam $I_{ik}$, where $i$ and $k$ – number of triangle element and number of incident beam for this element, respectively. Then, the program calculates the components of force and moment from the absorbed light as follows:

[TABLE]

where $A_{i}$ — area of triangle element; $\hat{\mathbf{n}}_{i}$ — normal to triangle element; $\hat{\mathbf{s}}_{ik}$ — direction vector for the incident beam; $\mathbf{r}_{i}$ — position vector of the current triangle element, pointed to the center of triangle.

After processing of beam absorption event the software calculates the set of new beams which has to be emitted from the element $A_{i}$. The outgoing beams (fig. 8) were considered similar to the incoming beams:

[TABLE]

where $m$ — number of outgoing beam.

The main vector and main moment of light pressure can be calculated as follows:

[TABLE]

where $N_{i}$ is a number of triangle elements; $N_{ik}$ — number of beams which were absorbed by the element with number $i$; $N_{im}$ — number of beams which were emitted by the element with number $i$.

The data on the main vector and main moment together with the original data on the direction of the external radiation vector were processed in the program implemented in GNU Octave software package Eaton et al (2015). The program did the numerical evaluation of expressions (47) and (48) and reconstructed tensors $\mathcal{I}$ and $\mathcal{K}$ and provided the necessary charts.

8.3. Approximation results

The principal force and moment of light pressure were calculated at 60 different orientations of the incident light. The resultant moment was calculated about the origin, located in the center of the outer edge of the outer heat shield.

In all cases, the estimated number of cast rays was equal to 1,000,000. It was assumed that the spacecraft radiators remain stationary when changing the orientation about the Sun, and solar panels are rotated to ensure the minimum angle of incidence of the light. All calculations neglected components of the spacecraft structure within the outer heat shield, as in these cases of orientation, when radiation from the sun falls into the bowl of the reflector, it is not compatible with the continued operation of the spacecraft as a telescope.

Fig. 9 — 12 represent the calculation results for absolute value of resultant force and moment. Fig. 9 — 10 correspond to the case when the entire spacecraft is fully specular, and Fig. 11 — 12 — when the entire surface of spacecraft is diffuse. The angle ranges $[-90\degree,+90\degree]$ corresponds to the ordinary orientation of space observatory. As it follows from the figures, with the increasing of numbers of terms in series (35) and (36) the force and moment approximation accuracy are also increasing. Depending on the task the number of terms in the series (35) and (36) may be restricted to achieve the acceptable level of approximation accuracy. Fig. 13 and 14 represent the components of principal force and moment in the specular and diffuse cases, respectively. In this model, the approximation of the resultant force converges fast enough to the known values ​​obtained by stochastic simulation. In the case of the approximation of the moment, it may require a greater number of terms for a more precise approximation. It is worth noting that the largest error of approximation of the principal moment is observed in the areas corresponding to the abnormal orientation of the test spacecraft.

Using the results of the calculation for specular and diffuse cases, we can obtain approximations for the following three specular–diffuse cases, namely $s=0{.}75$ (Fig. 15), $s=0{.}5$ (Fig. 16) moreover, $s=0{.}25$ (Fig. 17), and compare them with the results from the raytracing The calculation results are presented in the figures below. As can be seen from Fig. 15 — 17, this method provides a good approximation comparing to the Monte - Carlo results, and it may be used in the calculation of the principal force and moment of the light pressure on the spacecraft with arbitrary shape, whose surface has specular–diffusive properties.

9. Conclusion

We got the analytical expression for the radiation pressure force on a body with the arbitrary geometrical shape. The proposed method can be used to analyze the dynamics of the center of mass movement in and around the center of mass of space objects, which are significantly affected by the radiation pressure. The proposed method allows a parameterization of the properties of the body, at least by the degree of specularity.

10. Acknowledgements

Authors would like to thank assistant professor Marchevsky I.K. from department “Applied Mathematics” of BMSTU, assistant Kotsur O.S. from department “Aerospace Systems” of BMSTU and assistant Goncharov D.A. from department “Theoretical Mechanics” of BMSTU for their valuable advice and discussions which lead to the creation of explained method.