Designing Thermal Systems and Material Selections for Columbia's LIONESS CubeSat
- Gilbert Tohme
- Dec 26, 2024
- 10 min read
Updated: Mar 18
This project focuses on the Columbia Space Initiative's (CSI) CubeSat (cube satellite) Mission team of designing its first satellite, LIONESS (Line Imaging Orbiter for Nanosatellite-Enabled Spectrographic Surveys). It will host an integral field unit (IFU) spectrograph, inspired by the ground-based models of Columbia University's Schiminovich Astronomy and Instrumentation Lab where we will study diffuse gases in nearby galaxies. This mission is supported by the University Nanosat Program and the CubeSat Launch Initiative, with a funded launch scheduled by 2027.
I collaborated with a sub-team of ~10 students within the CubeSat’s Structures/Thermal sub-team. My current work centers on designing thermal simulation and processing material selection for both 1U and 6U CubeSat configurations.
Material Selection Research
Upon first joining the CSI's CubeSat as part of my M.S. journey at Columbia, I first worked with my fellow team members to conduct thermal component research. The selection of materials for the CubeSat is critical for ensuring reliable performance under extreme space conditions. We have considered both passive and active thermal mitigation methods.
Passive methods, which are reliable and low-power, involve using materials that naturally dissipate heat without manual intervention. Active cooling methods, which are reserved for situations requiring precise temperature control, will require careful regulation to balance power consumption and performance.
Research suggested that materials must exhibit high thermal conductivity, low mass, and robust durability to withstand the harsh conditions of space. Specifically, they must endure UV radiation, extreme temperature fluctuations, and mechanical stresses. Lightweight, conductive materials like aluminum are excellent candidates due to their favorable properties.
Additionally, we are accounting for differences between 1U and larger configurations, such as 6U. The lighter mass of a 1U CubeSat may lead to greater temperature fluctuations, while its lower thermal inertia could potentially allow for better thermal conductivity and less structural stress over time. This analysis ensures that our material selection is tailored to both the specific size and operational needs of the CubeSat.
Python Simulations
To better understand the thermal dynamics of our CubeSat, we utilized Python simulations to model the system under simplified assumptions. The primary assumption was that the CubeSat follows a standard 1U configuration with dimensions of 10 cm x 10 cm x 10 cm. Intensive thermal properties, such as thermal conductivity, node mass fraction, emissivity, and absorptivity, were assumed to remain constant. For a 10-node simulation, we uniformly distributed the mass across all nodes, presuming a constant density throughout the CubeSat. Shown below is the math we conducted based on this:
Thermal Node Modeling for 1U CubeSat: Heat Transfer Dynamics and Assumptions
The chalkboard notes describe the approach for modeling heat transfer in a 1U CubeSat, focusing on thermal dynamics and heat distribution within a simplified node-based system.
Thermal Assumptions and Setup
The CubeSat is assumed to be positioned between the Sun and Earth, with heat transfer dominated by solar radiation Q_sun, Earth's albedo Q_albedo, and emitted radiation Q_emi,out. However, the temperature 𝑇 that defines the satellite's initial heat transfer conditions is unclear and may need clarification. Node T_33 is assigned no mass and no heat transfer, possibly to simplify the model at the central point.
Mass and Node Considerations Nodes T_i2, T_i3, and T_i4 represent mass dimensions, but in the future, we may need to add more nodes to better balance the distribution of mass and dimensions. Our model assumes the CubeSat's axis of symmetry corresponds to its hottest cross-section, simplifying the thermal analysis. Heat Transfer Model
The diagram we drew visually outlines our CubeSat surface as follows:
Top Surface: Q_sun enters the CubeSat, while Q_emi,out exits it.
Bottom Surface: Both Earth's radiation Q_earth and Q_albedo enter, while Q_emi,out exits.
Right and Left Surfaces: Q_emi,out is emitted outward on both sides.
The node orientations are as follows:
Outer Box Corner Nodes: T_11, T_15, T_55, T_51
Outer Box Edge-Aligned Nodes: T_12, T_14, T_25, T_45, T_54, T_52, T_41, T_21
Outer Box Center Nodes: T_13, T_35, T_53, T_31
Inner Box Corner Nodes: T_22, T_24, T_44, T_42
Inner Box Center Nodes: T_23, T_34, T_42, T_32
Central Node: T_33, at the very center of both boxes, experiences symmetric heat transfer in all four directions but no net heat transfer because of symmetry assumptions.
Assumptions
The Heat Transfer between horizontally symmetric nodes, such as T_32 and T_34, is negligible due to symmetry.
The Heat Transfer between vertically aligned nodes, such as T_23 and T_43, is not negligible and must be calculated.
Equation Simplifications
We are using symmetry to reduce the complexity of the heat transfer equations. For instance:
T_14 = T_12, which allows a factor of 2 to simplify the -Q_cond term.
The equations for nodes T_12 and T_14 can be used due to symmetry, reducing our total calculations.
We can further reduce the equation count by collapsing the representation of T_25 to T_15 and T_24 to T_14.
However, this current model has notable limitations. It does not account for complex 2D heat dynamics, as it essentially treats the CubeSat as a thermally homogenous "brick." Moreover, it lacks the ability to incorporate heat generation within the system and presents challenges when attempting to generalize results for other configurations or conditions. Despite these limitations, the Python simulations, alongside our manual calculations, provided valuable insights into how heat is distributed across the CubeSat. These efforts will serve as the foundation for future work, including more advanced simulations. Our Python code is shown below:
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
from scipy import interpolate
import pandas as pd
from scipy.linalg import cholesky, solve_triangular
# Constants
sigma = 5.67e-8 # Stefan-Boltzmann constant
# Orbital Parameters
H = 400 # altitude, must be greater than ~300 for calcs [km]
R = 6370 # Earth radius [km]
h = (H+R)/R # function of relative orbit radius
incl = 51.6 # inclination [deg]
beta = incl*np.pi/180 # angle of inclination of orbit wrt solar vector [rad]
beta_max = np.pi/2 - np.arccos(1/h)
P = 5420 # period of orbit [s]
# Heat Transfer Coefficients and Constants
S = 1370 # solar irradiance [W]
Fb0 = 1 # view factor of unity
Fbp = 0.8853 # view factor
albedo = 0.26 # albedo
T_earth = 288 # temperature of Earth [K]
emis_earth = 0.6 # emissivity of Earth
# Satellite Specifications
alpha = 0.7263 # weighted absorptivity of spacecraft
A = 0.06 # frontal area [m^2]
Aslice = (1/3)*((2*0.1*0.3)+(2*0.1*0.2)) # node slice exterior surface area [m^2]
m = 18 # mass of spacecraft[kg]
emis_sat = 0.5798 # weighted emissivity of satellite
cp = 896 # heat capacity of Al 6061-t6 [J/kg.K]
T_init = 10 # initial temperature [C]
L = 0.025 # length between nodes [m]
k = 150 # thermal conductivity of AL 6061 [W/m.K]
# Define Orbital Position, phi, and Eclipse Angle Start and Ends, phi_es and phi_ee respectively
phi = np.arange(0, 6*2*np.pi, 1/10000)
phi_es = np.pi - np.arccos(np.sqrt(h**2 - 1)/(h*np.cos(beta))) # angular position eclipse starts [rad]
phi_ee = np.pi + np.arccos(np.sqrt(h**2 - 1)/(h*np.cos(beta))) # angular position eclipse ends [rad]
# Define Orbital Position, phi, and Eclipse Angle Start and Ends, phi_es and phi_ee respectively
phi = np.arange(0, 6*2*np.pi, 1/10000)
phi_es = np.pi - np.arccos(np.sqrt(h**2 - 1)/(h*np.cos(beta))) # angular position eclipse starts [rad]
phi_ee = np.pi + np.arccos(np.sqrt(h**2 - 1)/(h*np.cos(beta))) # angular position eclipse ends [rad]
# Define Orbital Position, phi, and Eclipse Angle Start and Ends, phi_es and phi_ee respectively
phi = np.arange(0, 6*2*np.pi, 1/10000)
phi_es = np.pi - np.arccos(np.sqrt(h**2 - 1)/(h*np.cos(beta))) # angular position eclipse starts [rad]
phi_ee = np.pi + np.arccos(np.sqrt(h**2 - 1)/(h*np.cos(beta))) # angular position eclipse ends [rad]
# Define view factors as a function of time
t = ViewFactors['t']
Fa = ViewFactors['Albedo View Factor']
VF = interpolate.interp1d(t, Fa)
# Create Heat Flux Vectors
def Q_S(t):
phi = t*(2*np.pi * 60)/(P)
if ((phi_es < (phi%(2*np.pi))) & ((phi%(2*np.pi)) < phi_ee)):
return 0
else:
return Q_S0
def Q_A(t):
phi = t*(2*np.pi * 60)/(P)
if ((phi%(2*np.pi))< phi_es) & ((phi%(2*np.pi)) < np.pi):
return Q_A0*VF(t)
elif ((phi_ee < (phi%(2*np.pi))) & ((phi%(2*np.pi)) > np.pi)):
return Q_A0*VF(t)
else:
return 0
def Q_E(t):
return Q_E0
# Sanity Check Functions When Compared to Steps Taken Above
phi = np.arange(0, 3*2*np.pi, 1/10000)
t = phi*P/(2*np.pi * 60)
qs = []
time = 0
for time in t:
qs.append(Q_S(time))
time += 1
qa = []
time = 0
for time in t:
qa.append(Q_A(time))
time += 1
qe = []
time = 0
for time in t:
qe.append(Q_E(time))
time += 1
qtot = np.array(qs)+np.array(qa)+np.array(qe)
fig, ax = plt.subplots(figsize=(10,5))
lw = 2
plt.rc('font', size=12)
ax.plot(t, qs, linewidth=lw, label = 'Solar Radiation Input [W]')
ax.plot(t, qa, linewidth=lw, label = 'Albedo Radiation Input [W]')
ax.plot(t, qe, linewidth=lw, label = 'Earth Radiation Input [W]')
ax.plot(t, qtot, linewidth=lw, label = 'Total Radiation Input [W]')
ax.set(xlabel='Time in Orbit [Minutes]',
ylabel='Radiation [W]',
title='Total Radiation Received By Satellite Over Time')
ax.grid()
plt.legend(loc=0)
plt.show()This first bit of code outputs the following graph:
Next:
first_orbit_albedo = np.array(qa)[t<=P]
len(first_orbit_albedo[first_orbit_albedo==0])
Output: 59651
len(t)
Output: 188496
(5420*19883/62832)/60
Output: 28.585715346744756
min(qtot), max(qtot)
Output: (7.208156145995331, 80.65207087707532)
# Create vector field for system of transient temperature ODEs without active cooling
def temp_ODE_vectorfield(w, t, p):
"""
Defines the differential equations for a 1D rod approximation of LIONESS using 5 coupled nodes.
Arguments:
w : vector of the state variables:
w = [T1, T2, T3, T4, T5, y1, y2, y3, y4, y5]
t : time
p : vector of the parameters:
p = Q_S, Q_A, Q_E, A, emis_sat, sigma, T_earth, m, cp, k, L
"""
T1, T2, T3, T4, T5, y1, y2, y3, y4, y5 = w
Q_S, Q_A, Q_E, A, emis_sat, sigma, T_earth, m, cp, k, L = p
# Create f = (T1', T2', T3', T4', T5', y1, y2, y3, y4, y5):
f = [y1,
y2,
y3,
y4,
y5,
(1/(m*cp))*(Q_S(t) + (k*A/L)*(T2-T1)-(emis_sat*A*sigma*(T1**4))),
(1/(m*cp))*((k*A/L)*(T3-T2)+(k*A/L)*(T1-T2)-(emis_sat*Aslice*sigma*(T2**4))),
(1/(m*cp))*((k*A/L)*(T4-T3)+(k*A/L)*(T2-T3)-(emis_sat*Aslice*sigma*(T3**4))),
(1/(m*cp))*((k*A/L)*(T5-T4)+(k*A/L)*(T3-T4)-(emis_sat*Aslice*sigma*(T4**4))),
(1/(m*cp))*(Q_A(t) + Q_E(t)+((k*A/L)*(T4-T5))-(emis_sat*A*sigma*(T5**4))),]
return f
# Create vector field for system of transient temperature ODEs without active cooling
def temp_ODE_vectorfield(w, t, p):
"""
Defines the differential equations for a 1D rod approximation of LIONESS using 5 coupled nodes.
Arguments:
w : vector of the state variables:
w = [T1, T2, T3, T4, T5, y1, y2, y3, y4, y5]
t : time
p : vector of the parameters:
p = Q_S, Q_A, Q_E, A, emis_sat, sigma, T_earth, m, cp, k, L
"""
T1, T2, T3, T4, T5, y1, y2, y3, y4, y5 = w
Q_S, Q_A, Q_E, A, emis_sat, sigma, T_earth, m, cp, k, L = p
# Create f = (T1', T2', T3', T4', T5', y1, y2, y3, y4, y5):
f = [y1,
y2,
y3,
y4,
y5,
(1/(m*cp))*(Q_S(t) + (k*A/L)*(T2-T1)-(emis_sat*A*sigma*(T1**4))),
(1/(m*cp))*((k*A/L)*(T3-T2)+(k*A/L)*(T1-T2)-(emis_sat*Aslice*sigma*(T2**4))),
(1/(m*cp))*((k*A/L)*(T4-T3)+(k*A/L)*(T2-T3)-(emis_sat*Aslice*sigma*(T3**4))),
(1/(m*cp))*((k*A/L)*(T5-T4)+(k*A/L)*(T3-T4)-(emis_sat*Aslice*sigma*(T4**4))),
(1/(m*cp))*(Q_A(t) + Q_E(t)+((k*A/L)*(T4-T5))-(emis_sat*A*sigma*(T5**4))),]
return f
# Create vector field for system of transient temperature ODEs without active cooling
def temp_ODE_vectorfield(w, t, p):
"""
Defines the differential equations for a 1D rod approximation of LIONESS using 5 coupled nodes.
Arguments:
w : vector of the state variables:
w = [T1, T2, T3, T4, T5, y1, y2, y3, y4, y5]
t : time
p : vector of the parameters:
p = Q_S, Q_A, Q_E, A, emis_sat, sigma, T_earth, m, cp, k, L
"""
T1, T2, T3, T4, T5, y1, y2, y3, y4, y5 = w
Q_S, Q_A, Q_E, A, emis_sat, sigma, T_earth, m, cp, k, L = p
# Create f = (T1', T2', T3', T4', T5', y1, y2, y3, y4, y5):
f = [y1,
y2,
y3,
y4,
y5,
(1/(m*cp))*(Q_S(t) + (k*A/L)*(T2-T1)-(emis_sat*A*sigma*(T1**4))),
(1/(m*cp))*((k*A/L)*(T3-T2)+(k*A/L)*(T1-T2)-(emis_sat*Aslice*sigma*(T2**4))),
(1/(m*cp))*((k*A/L)*(T4-T3)+(k*A/L)*(T2-T3)-(emis_sat*Aslice*sigma*(T3**4))),
(1/(m*cp))*((k*A/L)*(T5-T4)+(k*A/L)*(T3-T4)-(emis_sat*Aslice*sigma*(T4**4))),
(1/(m*cp))*(Q_A(t) + Q_E(t)+((k*A/L)*(T4-T5))-(emis_sat*A*sigma*(T5**4))),]
return fThis part of the code setting up a vector field for a system of transient temperature ODEs outputs the following graph:
fig, ax = plt.subplots(figsize=(10,5))
T_diff = T1_t-T5_t
ax.plot(t, T_diff)
ax.set(xlabel='Time in Orbit [Minutes]',
ylabel='Temperature [K]',
#xlim=[0, 100],
ylim=[min(T_diff)-0.5, max(T_diff)+0.5],
title='Temp Difference Between Sun-Facing and Earth-Facing Surfaces [K]')
ax.grid()
plt.show()The temperature difference after computing the ODEs is shown as this:
def temp_cooling_ODE_vectorfield(w, t, p):
"""
Defines the differential equations for a 1D rod approximation of LIONESS using 5 coupled nodes.
Arguments:
w : vector of the state variables:
w = [T1, T2, T3, T4, T5, y1, y2, y3, y4, y5]
t : time
p : vector of the parameters:
p = Q_S, Q_A, Q_E, A, emis_sat, sigma, T_earth, m, cp, k, L
"""
T1, T2, T3, T4, T5, y1, y2, y3, y4, y5 = w
Q_S, Q_A, Q_E, A, emis_sat, sigma, T_earth, m, cp, k, L = p
# Create f = (T1', T2', T3', T4', T5', y1, y2, y3, y4, y5):
f = [y1,
y2,
y3,
y4,
y5,
(1/(m*cp))*(Q_S(t) + (k*A/L)*(T2-T1)-(emis_sat*A*sigma*(T1**4))),
(1/(m*cp))*((k*A/L)*(T3-T2)+(k*A/L)*(T1-T2)-(emis_sat*Aslice*sigma*(T2**4))),
(1/(m*cp))*((k*A/L)*(T4-T3)+(k*A/L)*(T2-T3)-(emis_sat*Aslice*sigma*(T3**4))),
(1/(m*cp))*((k*A/L)*(T5-T4)+(k*A/L)*(T3-T4)-(emis_sat*Aslice*sigma*(T4**4))),
(1/(m*cp))*(Q_A(t) + Q_E(t)+((k*A/L)*(T4-T5))-(emis_sat*A*sigma*(T5**4))-50),]
return f
# ODE solver parameters
abserr = 1.0e-8
relerr = 1.0e-6
phi = np.arange(0, 5*2*np.pi, 1/10000)
t = phi*P/(2*np.pi * 60)
# Initial Conditions
T1 = 273.15+T_init
T2 = 273.15+T_init
T3 = 273.15+T_init
T4 = 273.15+T_init
T5 = 273.15+T_init
y1 = 0
y2 = 0
y3 = 0
y4 = 0
y5 = 0
# Pack up the parameters and initial conditions:
p = [Q_S, Q_A, Q_E, A, emis_sat, sigma, T_earth, m, cp, k, L]
w0 = [T1, T2, T3, T4, T5, y1, y2, y3, y4, y5]
# Call the ODE solver.
wsol = odeint(temp_cooling_ODE_vectorfield, w0, t, args=(p,),
atol=abserr, rtol=relerr)
fig, ax = plt.subplots(figsize=(10,5))
plt.rc('font', size=12)
T1c_t = wsol[:, 0]
T2c_t = wsol[:, 1]
T3c_t = wsol[:, 2]
T4c_t = wsol[:, 3]
T5c_t = wsol[:, 4]
ax.plot(t, T1c_t, label='Sun-Facing Surface Temp [K]')
ax.plot(t, T2c_t,label='Inside Sun-Facing Temp [K]')
ax.plot(t, T3c_t,label='Satellite Center Temp Temp [K]')
ax.plot(t, T4c_t,label='Inside Earth-Facing Temp [K]')
ax.plot(t, T5c_t,label='Earth-Facing Surface Temp [K]')
ax.set(xlabel='Time in Orbit [Minutes]',
ylabel='Temperature [K]',
#xlim=[0, 100],
#ylim=[260, 287],
title='Temperatures Within Satellite As a Function of Orbital Period With 25 W Cooling Capability')
ax.grid()
plt.legend(loc=0)
plt.show()
This piece of code outputs this:
where:
-20.9*(70-51-(13-17.75)*0.001*83)**(-1)
Output: -1.0776389909380357
-(283e3/8.314)*(1/1500-1/298)-(257.3e3/(8.314*298))
Output:-12.319588155931157
#Cholesky Decompisition
# Define a symmetric positive-definite matrix `A` representing thermal conduction
# Simplified example for a 3-node system
A = np.array([[4, -1, 0],
[-1, 4, -1],
[0, -1, 3]])
# Define a vector `b` representing heat sources or boundary conditions
b = np.array([100, 200, 150])
# Perform Cholesky decomposition: A = L * L^T
L = cholesky(A, lower=True)
# Solve L * y = b using forward substitution
y = solve_triangular(L, b, lower=True)
# Solve L^T * x = y using back substitution
x = solve_triangular(L.T, y, lower=False)
print("Solution x:", x)
Solution x: [45.12195122 80.48780488 76.82926829]What's Next
Our next steps involve advancing our simulation efforts by transitioning from Python-based models to a more sophisticated tool, Thermal Desktop. This software will allow us to conduct more accurate and comprehensive thermal simulations, accounting for complex heat generation, multi-dimensional dynamics, and environmental interactions in space. By leveraging Thermal Desktop’s capabilities, we aim to refine our understanding of thermal behavior and develop effective mitigation strategies.
Additionally, we will evaluate and select appropriate thermal management techniques, both passive and active, to optimize heat dissipation and ensure the CubeSat's functionality. This work will also guide us in automating the generalization of the heat distribution model, enabling quicker identification of critical thermal spots. These efforts will bring us closer to achieving a robust thermal design capable of enduring the extreme conditions of space.
Comments