Ultra-long-distance power transmission and space-based energy capture

Research report

by Jacky Song

On behalf of the Project Elara Team

Emilliano Vilchis, Milo Schnak,

Charles Wang, Nicholas Schnorbus, Arnar Arnarsson, and

Luke Monteiro

Abstract

A prototype solar capture, microwave conversion, and transmission (SCMCT) system is proposed. System design considerations are described, and initial modeling of the SCMCT system with simulation results is shown. Further, a strategy for prototype testing to validate and extend upon current research is outlined. Approval of this strategy is requested to enable the continuation of this research.

1Overview

The Sun is a star, a nuclear furnace with a power output1 of \(3.828 \times 10^{26} \text{ W}\) [7]. Meanwhile, Earth's global power consumption2 in 2019 was \(2.61 \times 10^{12} \text{ W}\) [5], 14 orders of magnitude below the power output of the sun. Even the most powerful artificial fusion reactor could hardly match the 4.6 billion-year-old, 2 trillion-quadrillion-ton one at the center of the Solar System. However, Earth-based solar power, while highly successful, has not truly delivered on the full possibilities of harnessing the energy of the Sun.

The primary difficulties of space-based solar power are twofold. First, high upfront costs and the difficulties of launching space infrastructure present immense hurdles. Second, the technological capabilities that would allow the full realization of space-based power to the point it is worth the cost are not yet possible. Therefore, the primary goal of our present research is the development of a space-based power transmission system that provides:

It is without a doubt that the technologies involved in this system are far from reaching readiness. However, Abiri et. al. [1] have already demonstrated small-scale versions of this technology in low-Earth orbit. This research intends to supplement present achievements in the field with further advances in realizing space-based power.

2Theoretical analysis and results

Due to the difficulties of practical testing, a theoretical model of wireless power capture and transmission was first developed for initial analysis. The conceptual system consists of mirrors and satellites situated in highly eccentric geosynchronous orbits. The individually-spaced mirrors form a composite parabolic solar collector, concentrating sunlight to a fixed point. A secondary mirror at the foci of the composite collector, segmented and also of parabolic design, directs the focused light into specialized transmitter satellites. These satellites then transmit the power to Earth in the form of microwaves. The system is planned to be only one of many, with the intent to build a space-based stellar power collection swarm around geosynchronous and eventually higher orbits which may one day be able to provide all the energy for the Earth. A simplified diagram of the system is shown below:

Figure 1. An illustration of a portion of the conceptual SCMCT system

For the theoretical modelling, the special case of Maxwell's equations in free space, namely the 3D electromagnetic wave equation, was used. The electromagnetic wave equation is given by:

\(\displaystyle \frac{\partial^2}{\partial t^2} \left( \begin{array}{c} \mathbf{E}\\ \mathbf{B} \end{array} \right) = c^2 \nabla^2 \left( \begin{array}{c} \mathbf{E}\\ \mathbf{B} \end{array} \right)\)

The proposed power transmission system operates at a frequency of 1 GHz, which has been identified as a candidate frequency for microwave power transmission due to its low atmospheric attenuation [6]. On initial testing, it was found that numerical simulations of the electromagnetic wave equation exhibited rapid temporal oscillations that resulted in numerical instabilities. Therefore, the time-independent variant of the electromagnetic wave equation was used instead. As is well-known, this is obtained by the separation of variables of the electromagnetic wave equation, upon the simplifying assumption that the time variation of the electric field follows \(\mathbf{E} \sim \cos \omega t\). This results in the Helmholtz equation:

\begin{eqnarray*} \nabla^2 \mathbf{E} + k^2 \mathbf{E} = 0, & \nabla^2 \mathbf{} \mathbf{B} + k^2 \mathbf{B} = 0 & \end{eqnarray*}

From this point on, only the expressions for the electric field are written and it is implied that the same expressions are true for the magnetic field up to the replacement \(\mathbf{E} \rightarrow \mathbf{B}\). The Helmholtz equation was prepared in variational form (weak form) for obtaining a numerical solution by the finite element method. By multiplication with a test function \(\Phi (\mathbf{x})\), integrating over the simulation domain \(\Omega\), and simplifying via the Divergence Theorem, one obtains the general variational form:

\(\displaystyle k^2 \int_{\Omega} \Phi \cdot \mathbf{E} \hspace{0.17em} dV - \int_{\Omega} \nabla_J \mathbf{E} : \nabla_J \Phi \hspace{0.17em} dV + \int_{\partial \Omega} \Phi \cdot \nabla_J \mathbf{E} \hspace{0.17em} d \mathbf{A} = 0\)

Where : is the double dot product (tensor contraction for matrices). The boundary conditions must then be incorporated for a meaningful solution. As a starting case, the boundary conditions are given as follows:

\begin{eqnarray*} \frac{\partial \mathbf{E}}{\partial n} |_{\partial \Omega_P} = \frac{\partial \mathbf{E}}{\partial n} |_{\partial \Omega_W} = 0, & \Omega \in \mathbb{R}^3 & \\ \lim_{r \rightarrow \infty} \mathbf{E} = \epsilon & & \end{eqnarray*}

Where, respectively:

The open boundary condition \(\partial \Omega_B\) is approximated by a technique described on the Comsol finite-element software blog [3], in which a coordinate transform \(x \to \alpha \tanh x, y \rightarrow \beta \tanh y\) is applied to extend the domain out to infinity, where \(\alpha\) and \(\beta\) are respectively half the length and width of the domain bounding box. This approximately replicates an open boundary that allows outgoing electromagnetic waves to radiate outwards in all directions. However, the direct application of the coordinate transformation results in significant mathematical complexities due to the dependence of the representation of the variational form on the choice of coordinates. Thus tensors will be used instead to give expressions of physical laws that are invariant under coordinate transformations. From this point on, all expressions of the weak form will be given in tensors unless otherwise specified. All tensors are assumed to be within Euclidean space where upper and lower indices are equivalent, that is, \(A^i = A_i\). The Einstein summation convention is assumed, in which repeated indices are implicitly summed over, and all indices take the numerical values of \(i = 1, 2, 3\) unless otherwise specified.

To begin, the weak form may be expressed with tensors as:

\(\displaystyle k^2 \int_{\Omega} \Phi_i E^i \hspace{0.17em} dV - \int_{\Omega} \partial_j E^i \partial^j \Phi_i \hspace{0.17em} dV + \int_{\partial \Omega} \Phi_i \partial_j E^i \hspace{0.17em} dA^j = 0\)

One may reduce by one dimension to just consider the 2D case, as was done to simplify the problem for this initial stage of research:

\(\displaystyle k^2 \int_{\Omega} \Phi_i E^i \hspace{0.17em} dA - \int_{\Omega} \partial_j E^i \partial^j \Phi_i \hspace{0.17em} dA + \int_{\partial \Omega} \Phi_i \partial_j E^i \hspace{0.17em} dx^j = 0\)

The old coordinates are denoted \(x^i = \mathbf{x} = (x, y)\), and new coordinates denoted \(x^k = \mathbf{x}' = (x', y')\) where \(x' = x' (x)\) and \(y' = y' (y)\). A change of variables was applied on the weak form \(x^i \to x^k\). On the first integral term, the Kronecker delta was used to relabel indices from \(i \to k\), by the relationships \(\Phi_i = \delta_i^k \Phi_k, \quad E^i = \delta^i_k E^k\), from which one can substitute into the first integral term to obtain:

\(\displaystyle k^2 \int_{\Omega} \Phi_i E^i \hspace{0.17em} dA = k^2 \int_{\Omega} \delta_i^k \Phi_k \hspace{0.17em} \delta^i_k E^k \hspace{0.17em} dA\)

The two Kronecker deltas are contracted over both of their indices, so \(\delta_i^k \hspace{0.17em} \delta^i_k = \delta^i_i = \delta^1_1 + \delta^2_2 + \delta^3_3 = 3\). Therefore the first term further simplifies to \(\displaystyle 3 k^2 \int_{\Omega} \Phi_k \hspace{0.17em} E^k \hspace{0.17em} dA\). For the second term, the chain rule \(\partial_j = \dfrac{\partial x^k}{\partial x^j} \dfrac{\partial}{\partial x^k} = \partial_j x^k \hspace{0.17em} \partial_k\) may be used to rewrite derivatives in terms of the new coordinates \(x^k\). After substitution one finds:

\(\displaystyle \int_{\Omega} \partial_j E^i \partial^j \Phi_i \hspace{0.17em} dA = \int_{\Omega} (\partial_j x^k) \partial_k E^i (\partial^j x_k) \partial^k \Phi_i \hspace{0.17em} dA\)

Where, after rearrangement and using the Kronecker delta for relabeling the \(E\) and \(\Phi\) indices (as shown previously with the first term), becomes:

\(\displaystyle \int_{\Omega} (\partial_j x^k) (\partial^j x_k) \partial_k \delta^i_k E^k \partial^k \delta_i^k \Phi_k \hspace{0.17em} dA\)

Taking out the Kronecker deltas allows for further simplification, where \(\delta^i_k \delta_i^k = 3\) due to the double contraction previously explained:

\(\displaystyle \int_{\Omega} \delta^i_k \delta_i^k (\partial_j x^k) (\partial^j x_k) \partial_k E^k \partial^k \Phi_k \hspace{0.17em} dA = 3 \int_{\Omega} (\partial_j x^k) (\partial^j x_k) \partial_k E^k \partial^k \Phi_k \hspace{0.17em} dA\)

The third term may be simplified via the tensor transformation law \(dx^j = \dfrac{\partial x^j}{\partial x^k} dx^k = \partial_k x^j dx^k\), from which one may obtain:

\(\displaystyle \int_{\partial \Omega} \Phi_i \partial_j E^i \hspace{0.17em} dx^j = \int_{\partial \Omega} \Phi_i \partial_j E^i \partial_k x^j dx^k\)

Now substituting in \(\partial_j = \partial_j x^k \hspace{0.17em} \partial_k\) and the two Kronecker delta identities \(\Phi_i = \delta_i^k \Phi_k\) and \(E^i = \delta^i_k E^k\) the result is:

\(\displaystyle \int_{\partial \Omega} \delta_i^k \Phi_k \partial_j x^k \hspace{0.17em} \partial_k \delta^i_k E^k \hspace{0.17em} \partial_k x^j dx^k\)

Where, after another double contraction over the Kronecker deltas and rearranging, noting that \(\partial_j x^k \hspace{0.17em} \partial_k x^j = 1\), gives:

\(\displaystyle 3 \int_{\partial \Omega} \Phi_k \partial_j x^k \hspace{0.17em} \partial_k x^j \hspace{0.17em} \partial_k E^k \hspace{0.17em} dx^k = 3 \int_{\partial \Omega} \Phi_k \partial_k E^k \hspace{0.17em} dx^k\)

Altogether, after substitution of all simplified terms and ordering the terms such that the quasi-linear term is second due to software requirements, one obtains the weak form in the transformed coordinates:

\(\displaystyle - 3 \int_{\Omega} \partial_j x^k \hspace{0.17em} \partial^j x_k \hspace{0.17em} \partial_k E^k \partial^k \Phi_k \hspace{0.17em} dA + 3 k^2 \int_{\Omega} \Phi_k \hspace{0.17em} E^k \hspace{0.17em} dA + 3 \int_{\partial \Omega} \Phi_k \partial_k E^k \hspace{0.17em} dx^k = 0\)

After which dividing by three from all sides gives the most simplified general form:

\(\displaystyle - \int_{\Omega} \partial_j x^k \hspace{0.17em} \partial^j x_k \hspace{0.17em} \partial_k E^k \partial^k \Phi_k \hspace{0.17em} dA + k^2 \int_{\Omega} \Phi_k \hspace{0.17em} E^k \hspace{0.17em} dA + \int_{\partial \Omega} \Phi_k \partial_k E^k \hspace{0.17em} dx^k = 0\)

In typical vector calculus notation rather than tensor notation, this may alternatively be written as:

\(\displaystyle - \int_{\Omega} (\nabla_J \mathbf{x}' : \nabla_J \mathbf{x}') (\nabla \cdot \tilde{\mathbf{E}}) (\nabla \cdot \Phi) \hspace{0.17em} dA + k^2 \int_{\Omega} \Phi \cdot \tilde{\mathbf{E}}\, d \mathbf{A} + \int_{\partial \Omega} \Phi \cdot (\nabla \cdot \tilde{\mathbf{E}}) \hspace{0.17em} d \mathbf{x} = 0\)

When applying boundary conditions to the boundary integral, the only contribution is that of the radiative boundary condition, which reduces to the constant Dirichlet boundary condition \(\widetilde{\mathbf{E}} |_{\partial \Omega_B} = \epsilon \) and thus:

\(\displaystyle \int_{\partial \Omega} \Phi \cdot (\nabla \cdot \mathbf{E}) \hspace{0.17em} d \mathbf{x} = 0\)

And therefore the boundary integral vanishes. One may also show, after tedious derivation, that the Laplacian in a generalized coordinates transform \(x^i \to x^j\) takes the form:

\(\displaystyle \nabla^2 = \partial^i \partial_i x^j \partial_j + \partial^i x_j \partial_i x^j \partial^j \partial_j\)

So with the coordinate transforms \(x' = \alpha \tanh x\), \(y' = \beta \tanh y\), the transformed Laplacian becomes:

\(\displaystyle \nabla^2 = \alpha^2 {(1 - x'}^2) \left[ \frac{\partial}{\partial x'} \left( {(1 - x'}^2) \frac{\partial}{\partial x'} \right) \right] + \beta^2 {(1 - y'}^2) \left[ \frac{\partial}{\partial y'} \left( {(1 - y'}^2) \frac{\partial}{\partial y'} \right) \right]\)

Thus the Helmholtz equation in the aforementioned transformed coordinates has the form:

\(\displaystyle \alpha^2 {(1 - x'}^2) \left[ \frac{\partial}{\partial x'} \left( {(1 - x'}^2) \frac{\partial}{\partial x'} \right) \right] + \beta^2 {(1 - y'}^2) \left[ \frac{\partial}{\partial y'} \left( {(1 - y'}^2) \frac{\partial}{\partial y'} \right) \right] \tilde{\mathbf{E}} (x', y') + k^2 \tilde{\mathbf{E}} (x', y') = 0\)

The geometry of the simulation was parametrized by four constants, the simulation volume width \(D\) and length \(L\), the primary reflector radius \(R_1\), secondary reflector radius \(R_2\), opening gap radius \(b\) (where the opening is located at the rear of the primary reflector dish), and the primary reflector anchor point \(d_1\). The primary collector and secondary mirror were respectively parametrized as follows:

\begin{eqnarray*} x_1 (t) = d_1 - \frac{1}{2 R_1} t^2, \quad y_1 (t) &= t, &t \in [- R_1, - b] \cup [b, R_1] \\ x_2 (t) = d_1 - \frac{R_1}{2} + \frac{1}{2 R_1} t^2, \quad y_2 (t) &= t, &t \in [-R_2, R_2] \end{eqnarray*}

The bounding box of the simulation were parametrized as follows, for which \(t \in [0, 1]\):

\(\displaystyle \begin{align*} x_T (t) &= Lt - \frac{L}{2}, \quad y_T (t) = \frac{D}{2}\\ x_B (t) &= Lt - \frac{L}{2}, \quad y_B (t) = - \frac{D}{2}\\ x_L (t) &= - \frac{L}{2}, \quad y_L (t) = Dt - \frac{D}{2}\\ x_R (t) &= \frac{L}{2}, \quad y_R (t) = Dt - \frac{D}{2} \end{align*}\)

Simulations were done in FreeFEM++ software [4], with a custom FreeFEM-based solver. This solver was written in the FreeFEM language and utilizes FreeFEM++'s built-in finite-element solver and mesh discretizer, but otherwise is completely custom and tailored to the specific simulation parameters. A solution in transformed \((x', y')\) coordinates was first obtained, before the built-in polynomial interpolation was used to transform the solution to \((x, y)\) coordinates, as shown in the subsequent figures.

Figure 2. Solution in transformed coordinates \((x', y')\)

Figure 3. Solution in physical-space \((x, y)\) coordinates

The simulation was intended as a proof-of-concept and thus it should be emphasized that the simulation results were not of specific significance and should not be intended as research data. However, several distinctive and possibly promising features did emerge out of the initial simulation. The electric field strength being highest approximately between the parabolics, and decaying quickly outside of the immediate vicinity, corresponds with the expected result. However, it should be noted that the rather weak field variation, in which the differences of the highest and lowest field strengths is merely \(5 \times 10^{- 3} \text{N/C}\), does not provide sufficient support for these claims. Thus these are merely deduction from observations and cannot be proven at present without further testing.

As a separate test of the fidelity of the results, a custom FreeFEM-based solver was tested against a finite-difference simulation in a similar but much simpler Helmholtz boundary-value problem. A solver written with the Python findiff library [2], a library that implements finite difference methods for PDEs, was tested, alongside a slightly different finite-difference solver that used findiff for simply derivative discretization. The two solvers differed at most by \(6.64 \times 10^{- 13} \text{N/C}\) and as such the rest of the finite-difference solution results are given by that of the primary findiff solver. The validation test took the form of solving the scalar-valued Helmholtz equation on the unit square \(\Omega = [0, 1] \times [0, 1]\), with the following boundary conditions, which can be shown to yield a non-trivial solution:

\begin{eqnarray*} E (x, 0) = 3 & & E (x, 1) = 10\\ E (0, y) = 6 & & E (1, y) = 1 \end{eqnarray*}

A FreeFEM-based solver, similar although not identical to the aforementioned fully-featured vector solver, was compared against the findiff solver. The results are shown below:

Figure 4. Numerical solution obtained by FreeFEM-based (finite element) solver

Figure 5. Numerical solution obtained by findiff (finite-difference) solver

Comparison of the numerical solutions from the two solvers showed excellent agreement and the completion of an important intermediary diagnostic for the primary FreeFEM-based solver. Diferences between the two solutions showed almost zero difference on a number of test points evaluated from both solutions, as shown below:

Figure 6. Comparison between the FEM and FDM solvers on various test points

The results suggest that the solver shows promise, at least to the extent of present testing. However, the author would like to note that a further validation test of the vector-valued Helmholtz equation has not been successfully completed as of the time of writing. Difficulties with performing an inverse coordinate transformation from the finite-difference numerical solution have yet to be resolved, among other issues. Two different methods have been attempted, that of polynomial interpolation as well as a neural network fit, but both do not meet research standards and much further work is necessary. In addition, it was found that the order of integrals was a significant factor in deciding the result of the solution, as well as the sign convention used, where there were significantly different values when changed. It is believed to be a result of the conventions of the FreeFem++ software, but should continue to be investigated.

3Planned upcoming research

Further theoretical work, including a more thorough treatment of the power collection system and theoretical modelling of the power transmission system, is necessary for advancing the research. More detailed studies of orbits and spacecraft design are also hoped for in upcoming work. In addition, construction and testing of physical prototypes to augument computational simulations is also essential. Prototypes are intended to be assembled via 3D printed components treated with electroless plating to form a smooth, metallic coating, after which specific tests may be performed. The verification of such prototypes in an experimental setting is highly anticipated to be carried out in the near future.

4Disclosures

The author would like to acknowledge the assistance of using LLMs as a source of valuable knowledge and an educational tool. No part of this paper was taken in whole or in part from content provided by the use of such tools. In addition, the author acknowledges the great assistance provided by the FreeFEM++ software documentation, including its many examples which included example solvers of the Helmholtz equation. Another excellent reference was the website of Nicolás Guarín-Zapata, whose work on finite-difference methods on infinite domain was invaluable.3

This paper is released into the public domain and the author withdraws all rights to this work, in accordance Project Elara guidelines. For more details about the project, please see https://elaraproject.github.io/

5Acknowledgements

On behalf of the Project and as the author of this paper, I would like to acknowledge the numerous contributions that went into this research. The support and collaborative effort of the entire Project team was essential to all of the research's successes. Milo played a key role in providing engineering insight, Luke in mathematical and physical analysis, Charles in offering fresh suggestions and ideas, Max (Emilliano) in continual work on research, and Arnar and Nicholas for their immense enthusiasm. I would also like to thank Professor Persans and Professor Newberg of Rensselaer Polytechnic for their support, knowledge, and their generously-offered resources, and Professor Tysoe for valuable information on the suggestion of electronless plating. I thank the support of my family and friends, including Gabrielle for her valuable suggestions of emphasizing the distinctive features of the project. Finally, I thank Dr. Harold White of Limitless Space Institute and Dr. Andrew Bramante of Greenwich High School for their invaluable mentorship and guidance.

Footnotes

1. According to the IAU defined value of the solar luminosity \(L_{\odot}\), after rounding. Back ↩

2. Found through averaging the provided 2019 global power consumption figure of 22,848 TWh, as stated by the IEA, over one year. Back ↩

3. See Finite differences in infinite domains https://nicoguaro.github.io/posts/infinite_fdm/. Back ↩

Bibliography

[1]

Abiri, B., Arya, M., Bohn, F., Fikes, A., Gal-Katziri, M., Gdoutos, E., Goel, A., Gonzalez, P. E., Kelzenberg, M., Lee, N., Marshall, M. A., Roy, T., Royer, F., Warmann, E. C., Vaidya, N., Vinogradova, T., Madonna, R., Atwater, H., Hajimiri, A., and Pellegrino, S., A lightweight space-based solar power generation and transmission satellite, https://arxiv.org/abs/2206.08373 (2022).

[2]

Baer, M., Findiff software package, https://github.com/maroba/findiff (2018).

[3]

Frei, W., Using perfectly matched layers and scattering boundary conditions for wave electromagnetics problems, https://www.comsol.com/blogs/using-perfectly-matched-layers-and-scattering-boundary-conditions-for-wave-electromagnetics-problems (2015).

[4]

Hecht, F. New development in freefem++. J. Numer. Math. 20, 3-4 (2012), 251–265. https://freefem.org/.

[5]

IEA. Electricity information: overview. Technical Report, International Energy Agency, Paris, 2021. https://www.iea.org/reports/electricity-information-overview.

[6]

ITU-R. Attenuation by atmospheric gases and related effects. Recommendation ITU-R P.676-13, International Telecommunications Union (Radiocommunication Sector), 2022. https://www.itu.int/rec/R-REC-P.676.

[7]

Mamajek, E. E., Prsa, A., Torres, G., Harmanec, P., Asplund, M., Bennett, P. D., Capitaine, N., Christensen-Dalsgaard, J., Depagne, E., Folkner, W. M., Haberreiter, M., Hekker, S., Hilton, J. L., Kostov, V., Kurtz, D. W., Laskar, J., Mason, B. D., Milone, E. F., Montgomery, M. M., Richards, M. T., Schou, J., and Stewart, S. G., Iau 2015 resolution b3 on recommended nominal conversion constants for selected solar and planetary properties, https://arxiv.org/abs/1510.07674 (2015).