CITE THIS NOTEBOOK: Structural stability of a lightsail for laser-driven interstellar flight by Dan-Cornelius Savu and Andrew Higgins. Wolfram Community OCT 20 2022.
ORIGINAL ARTICLE: Dan-Cornelius Savu, Andrew J. Higgins (2022), Structural stability of a lightsail for laser-driven interstellar flight. Acta astronautica.:376-393. https://doi.org/10.1016/j.actaastro.2022.09.003.
Article abstract: The structural stability of a lightsail under the intense laser flux necessary for interstellar flight is studied analytically and numerically. A sinusoidal perturbation is introduced into a two-dimensional thin-film sail to determine if the sail remains stable or if the perturbations grow in amplitude. A perfectly reflective sail material that gives specular reflection of the laser illumination is assumed in determining the resulting loading on the sail, although other reflection models can be incorporated as well. The quasi-static solution of the critical point between shape stability and instability is found by equating the bending moments induced on the sail due to radiation pressure with the restoring moments caused by the strength of the sail material and the tension applied at the edges of the sail. From this quasi-static solution, analytical expressions for the critical value of elastic modulus and boundary tension magnitude are found as a function of sail properties (e.g., thickness) and the amplitude and wave number of the initial sinusoidal perturbation. These same expressions are also derived from a more formal variational energy (virtual work) approach. A numerical model of the complete lightsail dynamics is developed by discretizing the lightsail into rectangular finite elements. By introducing torsional and rectilinear springs between the elements into the numerical model, a hierarchy of models is produced that can incorporate the effects of bending and applied tension. The numerical models permit the transient dynamics of a perturbed lightsail to be compared to the analytic results of the quasi-static analysis, visualized as stability maps that show the rate of perturbation growth as a function of sail thickness, elastic modulus, and applied tension. The analytic theory is able to correctly predict the stability boundary found in the numerical simulations. The stiffness required to make a thin lightsail stable against uncontrolled perturbation growth appears to be unfeasible for known materials, however, a relatively modest tensioning of the sail (e.g., via an inflatable structure or spinning of the sail) is able to maintain the sail shape under all wavelengths and amplitudes of perturbations.

Introduction

This notebook simulates the dynamics of a thin rectangular lightsail while undergoing acceleration due to radiation pressure applied by a laser beam of high (1–100
GW/
2
m
) intensity. The results of these simulations are compared to analytic expressions for the critical material stiffness and tension required to maintain the material stability of the lightsail. The question of lightsail structural stability arises due to the large laser loads that the sail is required to sustain in order to achieve relativistic speeds in a timely manner: given the large radiation pressure loads, will an ideally thin lightsail retain its structural integrity and remain functional despite the inevitable presence of perturbations causing non-uniform loading of the lightsail? These perturbations were simulated by sinusoidally deforming the discretized lightsail from its ideal flat shape prior to the beginning of the numerical simulation. If, within the simulation runtime, the maximum amplitude of the lightsail did not double from its original perturbation amplitude, then the lightsail configuration was termed stable oscillatory (oscillatory due to the fact that the maximum amplitude of stable lightsail configurations displays quasi-periodic behavior), with the lightsail configuration being termed potentially unstable otherwise. Although this notebook only explicitly features the tension and torsion (TnT) numerical model, any of the three numerical models described in the article can be used in this notebook to perform a lightsail simulation. That any of the three numerical models can be used in one notebook is achieved through the reduction of the more general TnT model into one the simpler numerical models. To produce the multiple simulations required to construct the stability maps found in the paper, an addendum Mathematica notebook was used to call this notebook in a loop while passing different lightsail parameter values for simulation upon each iteration (see Mathematica notebook attachment at the end of this post). For a more complete explanation of the approach used to generate the numerical lightsail models, the reader is encouraged to consult Section 3 of this notebook’s corresponding article.

Preamble

Packages

No additional packages are used throughout this notebook.

Units and Constants : all calculations are done in SI MKS units

Unit definitions below are done for inputting different units into functions.
When a value is stated in a non-MKS unit (meter-kilogram-seconds), multiply by the unit below to put the value into MKS.
If you want to put a value into a non-MKS unit, divide by the unit below.
◼
  • Length
    • In[]:=
    km=1000;cm=
    -2
    10
    ;mm=
    -3
    10
    ;µm=
    -6
    10
    ;nm=
    -9
    10
    ;AU=149597870700;
    ◼
  • Time
    • In[]:=
    minute=60;hour=60minute;day=24hour;year=365.25day;​​ms=
    -3
    10
    ;µs=
    -6
    10
    ;ns=
    -9
    10
    ;
    ◼
  • Power
    • In[]:=
    kW=1000;MW=1000kW;GW=1000MW;
    ◼
  • Pressure
    • In[]:=
    kPa=
    3
    10
    ;MPa=
    6
    10
    ;GPa=
    9
    10
    ;bar=
    5
    10
    ;kbar=
    5
    10
    ;Mbar=
    11
    10
    ;
    ◼
  • Definition of Degrees in radians
    • In[]:=
    Degrees=
    2π
    360
    ;
    ◼
  • Constants
    • In[]:=
    c0=299792458.;(*ThespeedoflightinvacuuminSIunits.*)
    In[]:=
    g0=9.8066;(*ThestandardgravityofEarthinSIunits.*)
    In[]:=
    LY=c0*year;(*AlightyearinSIunits.*)

    Setting Global Membership and Range Parameter Assumptions

    No additional global assumptions made throughout this notebook.
    In[]:=
    (*$Assumptions={};*)

    Graphics Settings

    In[]:=
    SetOptions[{ListPlot,Plot,ListLinePlot,ListLogLogPlot},​​LabelStyle{FontFamily"Times New Roman",FontSize18,FontColorBlack},​​Frame{{True,True},{True,True}},FrameTicksStyle{{{FontFamily"Times New Roman",FontSize18,FontColorBlack},Directive[FontOpacity0,FontSize0]},{{FontFamily"Times New Roman",FontSize18,FontColorBlack},Directive[FontOpacity0,FontSize0]}},​​AxesTrue,​​AxesStyle{Dashed,Thin},​​AspectRatio1/GoldenRatio,​​ImageSizeScaled[0.6]];

    Auto-Hide Input Cell Function

    To automatically collapse an input cell upon execution, simply write AutoCollapse[] in a newline—within the cell—right after the end of the input cell you want to collapse (code taken from https://mathematica.stackexchange.com/questions/680/how-to-keep-input-cells-hidden-after-evaluating-notebook):
    In[]:=
    AutoCollapse[]:=(If[$FrontEnd=!=$Failed,SelectionMove[EvaluationNotebook[],All,GeneratedCell];​​FrontEndTokenExecute["SelectionCloseUnselectedCells"]])

    Lightsail Simulation (User) Input Parameters

    This section constitutes the only section where the user is free to change variable numerical inputs in order to study a different lightsail numerical simulation configuration. A few variables that have expressions assigned to them are placed in this section for the sake of conceptual convenience. The user-input simulation and lightsail bulk properties in this section will be used to generate the equations of motion of the discrete lightsail model in the next section.
    ​
    This notebook explicitly features only the equations of motion of the tension and torsion (TnT) numerical lightsail model. However, because the TnT model is the more general of the three numerical models described in the article, the user can readily reduce the TnT model to either the torsion lightsail model (no rectilinear springs and no boundary tension) or the rigid lightsail model (no rectilinear springs, no boundary tension, and no torsion springs). Changing between numerical models can be done by the user through the assignment of different integer values to the numModel variable. The numModel possible input values and their corresponding numerical model are the following: 0 for the rigid lightsail model, 1 for the torsion lightsail model, and 2 for the TnT lightsail model. Note that the article originally made use of three separate notebooks for each lightsail numerical model, and the programming accommodation that was made to include all three models in one notebook results in small quantitative differences between the numerical output of this notebook and the results described in the article.
    In[]:=
    n=50;(*Thisvariablesetsthenumberofplatesusedforthediscretization.Theauthorsdonotrecommendrunningasimulationwithmorethan50elements.*)
    In[]:=
    tfinal=1;(*Totalsimulationruntime,[s].Theauthorsdonotrecommendrunningasimulationforlongerthan1second.*)
    In[]:=
    numModel=2;(*Thisvariableallowstheusertochangebetweenthethreelightsailnumericalmodels.Thepossibleinputsare0=RigidModel,1=TorsionModel,and2=TnTModel.*)

    Lightsail Geometry

    In[]:=
    thick=10µm;(*Lightsailthickness,[m].*)
    In[]:=
    L=1;(*Lightsaillength,[m].Typicallysetto1meter.*)
    In[]:=
    W=1;(*Lightsailwidth,[m].Typicallysetto1meter.*)

    Lightsail Material Properties

    In[]:=
    density=1000;(*Lightsailmassdensity,kg
    3
    m
    .*)
    In[]:=
    M=density*thick*W*L;(*Lightsail(full)mass,[kg].*)
    The Young’s modulus value used throughout the TnT simulations in the article was 5 GPa. For the torsion model, two sample modulus values are 1480 GPa and 0.108 GPa for a stable and an unstable (torsion model) lightsail configuration, respectively.
    In[]:=
    EModulus=5GPa;(*LightsailmaterialYoung'smodulus,[Pa].*)

    Initial Sinusoidal Perturbations Parameters

    The initial sinusoidal perturbations are of the form
    amplitude*Sin
    2*π*x
    λ
    
    .
    In[]:=
    ν=3/2;(*Initialsinusoidalperturbationsmodenumber,[-].*)
    In[]:=
    λ=N[L/ν];(*Initialsinusoidalperturbationswavelength,[m].*)
    In[]:=
    amplitude=0.01mm;(*Initialsinusoidalperturbationsamplitude,[m].*)

    Incident Light Properties

    In[]:=
    I0=10GW;(*Incidentlightintensity,[GW/
    2
    m
    ].*)
    In[]:=
    Pressure=
    I0
    c0
    ;(*Incidentradiationpressure,[Pa].Thespeedoflightisc0=299792458m/s.*)
    In[]:=
    acc=
    2*Pressure
    M
    ;(*Accelerationofaflat,perfectlyreflectivelightsail,[m/
    2
    s
    ].Variabletobeusedlaterinthedatapost-processing.*)

    Applied Boundary Tension Magnitude

    Sample boundary tension magnitude values are 0.01288 N and
    3.085*
    -6
    10
    N for a stable and an unstable (TnT model) lightsail configuration, respectively. If either the rigid lightsail model or the torsion lightsail model is chosen for simulation (numModel = 0 or numModel = 1, respectively), then the Tension variable will be assigned the value 0 in the next section.
    In[]:=
    Tension=0.01288;(*Appliedboundarytensionmagnitude,[N].*)

    Numerical Method Implementation

    In this section, the equations of motion of a lightsail discretized using
    n
    rigid slice elements are solved numerically using the Wolfram Mathematica NDSolve[] function. Based on the bulk lightsail properties defined by the user in the above section, the lightsail (pseudo) finite-element model properties such as lightsail slice element geometry, rectilinear spring stiffness, etc. are first computed. These computed discrete lightsail properties are then used to generate the full set of ordinary differential equations (ODEs) governing the motion of the numerical lightsail as well as the sinusoidal perturbation initial conditions. The ODEs as well as the initial conditions are then passed as arguments to NDSolve[] for numerical solution.

    Lightsail (Pseudo) Finite-Element Model Properties

    Lightsail Slice Element Geometry and Inertia

    In[]:=
    l=L/n;(*Lightsailsliceelementlength,[m].*)
    In[]:=
    m=density*thick*W*l;(*Lightsailsliceelementmass,[kg].*)
    In[]:=
    IG=(1/12)*m*(l^2);(*Lightsailsliceelementmomentofinertiaaboutitscenterofmass,kg*
    2
    m
    .*)
    In[]:=
    areaC=thick*W;(*Lightsailslicecross-sectionalarea,[
    2
    m
    ].*)

    Lightsail Slice Element Radiation-Driven Force

    In[]:=
    f0=2*
    I0
    c0
    *W*l;(*Lightsailsliceelementradation-drivenforcenominalmagnitude,[N].Nominalasintheforcemagnitudebeforethe
    2
    cos
    θ
    i
    modifierduetoinclination(appliedintheEOMsfoundbelow).*)

    Torsional Springs Stiffness

    To study only the effects of boundary tension upon lightsail structural stability, the torsional spring constant was typically set to zero when using the TnT model.
    In[]:=
    IArea=(1/12)*W*(thick^3);(*Lightsailsliceelementsecondmomentofareainthedirectionoftheinternalbendingmoments,[
    4
    m
    ].*)
    In[]:=
    KT=IfnumModel==0||numModel==2,0,N
    EModulus*IArea
    l
    ;(*Torsionalspringsconstant,[N*m].*)
    In[]:=
    (*KT=
    EModulus*IArea
    l
    ;*)

    Rectilinear Springs Stiffness

    In[]:=
    ls=0;(*Lengthofthespringsatrest,[m].Alwayssettozero.*)
    In[]:=
    KS=(n-1)*
    EModulus*areaC
    L
    ;(*Rectilinearspringsconstant,[N/m].*)
    The variable Tension is assigned the value 0 when the TnT model is not used for simulation.

    Lightsail Numerical Model Equations of Motion (EOMs)

    Using the discrete lightsail properties above, the equations of motion of the system are generated below. For a more detailed account as to their Lagrangian-based derivation, please see the notebook’s respective journal article.

    Constraint Equations

    Matrix Form Coefficients

    Sinusoidal Perturbation Initial Conditions

    Below is a display of the initial sinusoidal perturbation. The red dots represent the center of mass of each rigid element while the blue curve represents the ideal (analytical) sinusoidal shape.

    NDSolve

    Graphical Representation(s) of the Numerical Solution

    Following the numerical computation of the full dynamics of the lightsail, this section contains three graphical representations of the numerical lightsail behavior: static temporal snapshots of the lightsail, a plot of maximum amplitude of the lightsail as simulation time progresses, and two animation panels (each set in different non-inertial frames). As promised in the previous section, posk[] /. solns is used throughout the codes generating the graphical representations.

    Lightsail Simulation Static Snapshots

    The function below is used to generate temporal snapshots of the numerical lightsail. To better appreciate the shape dynamics of the lightsail, the snapshots are set in a reference frame attached to the center of mass of the numerical lightsail.

    Maximum Amplitude Growth Plot

    Animation Panel

    Following the static figures generated above, below are two animations of the motion of the discretized lightsail. As with the numerical lightsail snapshots, in order to better appreciate the lightsail shape dynamics, both animations are visualized in non-inertial (accelerating) reference frames. The first animation is set in a reference frame attached to a flat lightsail and the second animation is set in a reference frame attached to the center of mass of the discretized lightsail.

    Flat Lightsail Frame Animation

    Center of Mass (CoM) Frame Animation

    Animation Video Export

    The code commented out below can be used to generate and export high-quality animation videos of the lightsail numerical simulation (in the CoM reference frame). The high-quality animation is first generated frame by frame using a Do[] function with each frame being saved to a user-specified directory in .MX format. The .MX frames are then imported and assigned to the animationVideoFrames variable as members of a single-level list. The animationVideoFrames variable is then exported as a .mp4 video file. The individually exported frames collectively have a storage size of up to 8 MB and the .mp4 video file export has a storage size of up to 3 MB (when running simulations for no more than 1 second). If the code in this subsubsection is run on a device with 8 GB of RAM or less, then the authors recommend all other user applications besides this Wolfram Mathematica notebook to be closed to ensure no memory overload.