C++ Question Advice for Calculating Vessel's Temperature?

[Mr Martian]

Hi All,

I have been stuck on something for a few weeks now, and I am just going round in circles. My goal is to calculate a temperature reading based on some defined parameters (like surface area, emmisivity, etc.) Right now I am not getting into any aerodynamics or atmospheric thermodynamics, I firstly need to be able to calculate temperature in vaccuum conditions from sunlight alone.

This is what I have (runs every timestep):

double sA = 12.0;    // surface area
double E = 0.5;        // emmisivity
double sigma = 5.670373e-8;    //Stefan Boltzmann constant

double radiation = E * sigma * pow(T, 4) * sA;
T += ((sun * sA) - radiation) * simdt;

The parameter "sun" is energy from the sun in W/m^2 (about 1300 at Earth's distance).
The parameter "T" is a header-defined value for temperature in degrees Kelvin.

My issues are, the temperature almost instantly climbs to around 460K when in sunlight (this is when orbiting Earth, so the value seems right-ish, based on my parameters) but I do not understand why this temperature increase is so rapid. When out of sunlight, the temperature will again radiply drop to below 100K in under a second. I would have imagined a realistic climb and drop in temperature respectively.

Additionally, when in sunlight, the output is extremely unstable and will wildly fluctuate and return -1.#INF if any time acceleration is applied at all.

I know I am probably just missing something obvious, or made some simple formula errors, but I could really use a fresh perspective on this, I'm starting to tear my hair out!

Thank you all!
MrMartian

 

[francisdrake]

I think, calculating the heat flow is only the first step. You are assuming a quasi-stationary condition, which is ok for the minimum and maximum temperature, which is an equilibrium after long time.

The heating up and cooling down of the material takes some time. The heat capacity is defined by the specific heat x the mass of the object. How fast the heat enters the object is defined by the heat conductivity.

This is further complicated, if the objects wall is multi-layered. Like an aluminum outer shell, insulation panel in between, then an inner steel shell. The science for this is called transient (instationary) heat transfer. To calculate it is a major task.

To simplify this for your application you might define a lump of mass with an average heat capacity and calculate per second (or either a bigger timestep) how much the incoming heat flow has - on average - changed the temperature of the mass.

 

[n72.75]

What you're calculating with your "T" variable is energy, not temperature. Let's rename that Q.

To calculate temperature you need: mass, m; heat capacity, c; and the energy you are calculating.

Temperature = m*c*Q;
Temperature = Q/(m*c);
Or better
Temperature = Q/(m*c(T));
Where c(T) is specific heat capacity as a function of temperature.

This may be helpful to look at: https://github.com/n7275/NASSP/blob...o/src_sys/PanelSDK/Internals/Thermal.cpp#L164

Also, be aware that because you are effectively solving a differential equation, you will have some maximum timestep lengths at which your thermal simulation is no longer stable, and will instead oscillate. https://en.wikipedia.org/wiki/Stiff_equation This is why NASSP has time aceleration limits (also emulating computers). We actually break each Orbiter timestep into multiple "systems timesteps" to maintain stability. If you see oscillating behavior that seems dependent on timestep length, that's what this is. Fixed timestep mode can help you find this.

 

[jedidia]

First off, what you're trying to do is pretty complex. I kinda did it once, but one major problem you have is that not all of your vessel is exposed to the sun. One side is shadowing itself. Figuring that even just approximately is non-trivial. I had very simple approximation going in which the vessel was just treated as a cuboid. From there I had to calculate its orientation to the sun, and which surfaces of the cuboid are receiving how much energy (usually, 3 surfaces will be receiving energy, while the remaining 3 are not, but the surfaces receiving energy will receive very different amounts depending on their angle to the sun and their size). This will involve some 3-dimensional trigonometry. Orbiters API has a couple of functions to help with that.

And then there's material properties, which I can't see in your equation at all. @fracisdrake has mentioned the most important ones (appart from emissivity): Heat capacity and heat flux, though in all honesty, heat flux isn't that important for a ballpark estimate, especially if you just abstract the entire hull to a cuboid (which I would recommend, at least in a first iteration... the alternatives start getting very complex very fast).

Heat capacity on the other hand is the thing you need to fix your instantaneous temperature change problem. It will determine how many kelvin your surfaces will heat up per energy received. If it has a high heat capacity, it will heat up slower than if it has a low heat capaacity.
THe same goes for cooling down. Heat capacity is effectively a measure of how good something is at storing heat energy. So a material with a high heat capacity with the same temperature as a low one will have much more heat stored up in it once it gets into shadow, ergo it will take longer for it to radiate all that energy off again. THat's essentially the crux of your problem, as far as I can tell.

Once you got that you can think about refinements, like maybe putting in some heat flux between the individual surfaces of your cuboid, and breaking up your cuboid into several cuboids or maybe turning it into an 8-sided (or more) cross-section if your hull has a cylindrical geometry, so you get more surfaces to play with.

 

[francisdrake]

T += ((sun * sA) - radiation) * simdt;

I think, a hack would be to put a limit to the temperature increase per timestep:
---
dT = ((sun * sA) - radiation) * simdt;

if (dT > 0)
T += (dT>0.1? 0.1 : dT);
else
T += (dT < 0.1? -0.1 : dT);
----
If we look on the heating-up curve (dT > 0), it checks if dT is larger than a threshold value (here: 0.1 K per timestep).
If yes it limits dT to +0.1 K per timestep. If not, because the equilibrium temperature is close, the calculated dT value is taken.

• The limit is a bit arbitrary, unless you dervice it from measurements or theoretical calculations.
• You get a linear ramp instead of an e-function aproaching the end temperature.
• It will still be oscilating around the end temperature, but only within a limited bandwidth.

 

[Thunder Chicken]

I teach engineering heat transfer. What you are trying to do is challenging even for engineers to do with fidelity. To what the others have said, if you are looking for realism, you need to account for the thermal capacitance of the spacecraft, transient analysis, and several other radiation heat fluxes that you haven't included, and you need to understand that the entire vessel isn't at a single temperature, so you you need to determine what temperature you want to determine.

What you call 'radiation' is the emission of the spacecraft surface at temperature T, which is fine as far as that goes.

So at a given instant, the net heat transfer rate into the vessel surface would be Q_net = ((sun * sA) - radiation), not T.

Now the problem is to determine how much does that rate of heat transfer change the temperature of the vessel. This requires a lot of understanding about the spacecraft materials, transient analysis, etc., and ultimately what temperature you are interested in, as different parts of the spacecraft will be a wildly different temperatures. For example, shuttle style silica tiles are extremely poor conductors of heat, so the only part of the tile that heats up / cools down quickly is the exposed surface, so the involved thermal mass is small, and the temperature on the interior of the vessel would hardly change. In this event Q_net would be close to zero, and the surface temperature would be the equilibrium temperature at steady state.

If you are just interested in determining the surface temperature for an adiabatic surface, you would achieve quasi-equilibrium very quickly, and you could get the temperature by setting your 'radiation' = sun*sA and solving for T:

T = pow(sun/(E*sigma),1/4)

This WILL change very fast, as soon as the surface is exposed to the sun. It seems that your value of about 460K is correct for this situation.

Also note that parts of the spacecraft exposed only to deep space, or to the Earth, will experience other equilibrium temperatures.

If you want to mimic a thermal transient assuming a uniform temperature mass, you could do something like this

Q_net = ((sun * sA) - radiation). //Note that Q_net is a rate of heat transfer

dT = (Q_net*simdt)/(mass*spec_heat)

T += dT

The product of mass and the specific heat would be the thermal capacitance of a uniform temperature object, and in your case your best bet would be to tweak those numbers until you get something that seems "realistic". This would cause the temperature change to lag and not be instantaneous, but it is a crude hack, not really rigorous for realistic heat transfer.

 

[Thunder Chicken]

What you're calculating with your "T" variable is energy, not temperature. Let's rename that Q.

To calculate temperature you need: mass, m; heat capacity, c; and the energy you are calculating.

Temperature = m*c*Q;

This isn't correct. Q (heat or energy) = m*c*dT, so dT = Q/(m*c), and this isn't the temperature, but the change in temperature per time step.

I think, a hack would be to put a limit to the temperature increase per timestep:
---
dT = ((sun * sA) - radiation) * simdt;

This is also incorrect. This wouldn't give you dT, it would give you the net energy per time step.

 

[n72.75]

This isn't correct. Q (heat or energy) = m*c*dT, so dT = Q/(m*c), and this isn't the temperature, but the change in temperature per time step.

While this is true, strictly speaking (and I should have been more clear), you can also use this method to calculate object temperature as long as you pick a reasonable value for heat capacity. I am making the assumption that energy is linearly proportional to temperature, with the equation.

More realistically, you would impliment a function, c(T), but I don't think assuming that the heat capacity is constant for your spacecraft hull is wildly unrealistic, at least as far as the kinds of physical simulation you can do in an Orbiter addon go.

This is effectively saying

[math]Q = \int_{0}^{T}cmdT[/math]

 

[Thunder Chicken]

While this is true, strictly speaking (and I should have been more clear), you can also use this method to calculate object temperature as long as you pick a reasonable value for heat capacity. I am making the assumption that energy is linearly proportional to temperature, with the equation.

More realistically, you would impliment a function, c(T), but I don't think assuming that the heat capacity is constant for your spacecraft hull is wildly unrealistic, at least as far as the kinds of physical simulation you can do in an Orbiter addon go.

This is effectively saying

[math]Q = \int_{0}^{T}cmdT[/math]

That's only true if you are starting at T = 0. Otherwise the [imath]{\Delta}T[/imath] increases linearly with Q.

And thanks for pointing out we have LaTex capability for equations! I was not aware of that.

 

[n72.75]

That's only true if you are starting at T = 0. Otherwise the [imath]{\Delta}T[/imath] increases linearly with Q.

And thanks for pointing out we have LaTex capability for equations! I was not aware of that.

To clarify, I have no expectation that Q be value that is representative of the real world, necessarily. It's probably sort-of close at some temperatures, but it's really just there for "book keeping" of energy in and out of the object. If c is constant, Q and T are linearly proportional, so if you pick a good c value for the material, and temperature range your simulated object will spend most of it's time in, this will work just fine. Otherwise you can take some specific heat curves from NIST or AIChE publications and use those.

 

[Thunder Chicken]

To clarify, I have no expectation that Q be value that is representative of the real world, necessarily. It's probably sort-of close at some temperatures, but it's really just there for "book keeping" of energy in and out of the object. If c is constant, Q and T are linearly proportional, so if you pick a good c value for the material, and temperature range your simulated object will spend most of it's time in, this will work just fine. Otherwise you can take some specific heat curves from NIST or AIChE publications and use those.

The temperature dependence of the specific heat here is of least concern. We're trying to get something to the OP that gets them something in the order-of-magnitude ballpark, not trying to add more decimals of accuracy to a faithful model. Just writing a correct and comprehensive transient energy equation for this system would require more information than OP has or is likely ever to get.

 

[Mr Martian]

Thank you all so much for all the info! @n72.75 thank you for the clarification... I had completely overlooked mass And thank you for the info on stiff equations. I have encountered this before but I was not aware there was a name for this. I have something which works better now, but still need to do more to safeguard agains this.

@Thunder Chicken thank you for the comprehensive breakdown. As you mentioned, I am really trying to get something more in the 'realistic-enough' ballpark. What I have now seems to achieve thatpretty well.

I am actually implementing this as part of an SDK, so this is for a struct which I have called a 'thermal area'. Users can define numerous thermal areas for their vessel and specify its properties, namely effective cross-sections (used for exposed area as @jedidia mentioned), surface area, emissivity, and specific heat. I know that there is a lot more to calculating something like a hull temperature, but I figure these parameters are enough for any users of the SDK to tweak around to define some zones on their vessels. The current safeguards I have in place to prevent wild fluctuations at higher simdt rate work well, prividing no rediculous values are provided (tiny mass and massive surface area etc.)

Now I suppose the next step would be the even more complicated task of adjusting for atmospheric heating and cooling!

 

[Mr Martian]

I think I'm going to have to ask for some advice again...

As I previously mentioned, I want to now calculate the temperature affect from the ambient atmosphere. After this, I will then try to calculate for dynamic pressure/friction heating.

I figure a good place to start would be to get the atmospheric temperature, so I created a simple function to calculate this, posted below:

double GetAtmTemp()
{
    if (!m_hRefbody) return NULL;
    if (!oapiPlanetHasAtmosphere(m_hRefbody)) return NULL;

    double p0 = oapiGetPlanetAtmConstants(m_hRefbody)->rho0;    // get atm density at mean radius
    return (vParent->GetAtmTemperature() * (vParent->GetAtmDensity() / p0));
}

In this code, m_hRefbody is the current surface reference body, vParent is the VESSEL class instantiating my API's class.

This function works fine for Venus, Earth and Mars. The temperature is uniform and does not change with day/night obviously, but it's good enough for an average atmospheric temperature, however, on Io, Titan, Triton etc. it returns a simillar temperature to Earth. Is this the right way to go about this? or is the GetAtmTemperature() function's accuracy dependant on a planet's atmosphere module? Has anyone tackled what I'm trying to accomplish before?

As always I truly appreciate any advice,
Thanks!

 

[Thunder Chicken]

I think I'm going to have to ask for some advice again...

As I previously mentioned, I want to now calculate the temperature affect from the ambient atmosphere. After this, I will then try to calculate for dynamic pressure/friction heating.

I figure a good place to start would be to get the atmospheric temperature, so I created a simple function to calculate this, posted below:

double GetAtmTemp()
{
    if (!m_hRefbody) return NULL;
    if (!oapiPlanetHasAtmosphere(m_hRefbody)) return NULL;

    double p0 = oapiGetPlanetAtmConstants(m_hRefbody)->rho0;    // get atm density at mean radius
    return (vParent->GetAtmTemperature() * (vParent->GetAtmDensity() / p0));
}

In this code, m_hRefbody is the current surface reference body, vParent is the VESSEL class instantiating my API's class.

This function works fine for Venus, Earth and Mars. The temperature is uniform and does not change with day/night obviously, but it's good enough for an average atmospheric temperature, however, on Io, Titan, Triton etc. it returns a simillar temperature to Earth. Is this the right way to go about this? or is the GetAtmTemperature() function's accuracy dependant on a planet's atmosphere module? Has anyone tackled what I'm trying to accomplish before?

As always I truly appreciate any advice,
Thanks!

Earth, Venus and Mars have Atmosphere.cfg files whereas the others do not. The planet configuration files for these also don't contain any temperature information, just pressure and density for drag effects. Orbiter may be defaulting to Earth's Atmosphere.cfg when one for the particular planet isn't available. You'll have to read up on what is done in Orbiter in the Doc folder.