CODEWITHSUNDEEP
pythonmathsimulationphysicsgame-physics
Calrizan
Calrizan

Reputation: 13

Uncontrollable Simulation Oscillation

This is likely a math problem as much as it is a programming problem, but I seem to be encountering severe oscillations in temperature in my class method "update()" when warp is set for a high value (1000+) in the code below. All temperatures are in Kelvin for simplicity.

(I am not a programmer by profession. This formatting is likely unpleasant.)

import math

#Critical to the Stefan-Boltzmann equation. Otherwise known as Sigma
BOLTZMANN_CONSTANT = 5.67e-8

class GeneratorObject(object):
    """Create a new object to run thermal simulation on."""
    def __init__(self, mass, emissivity, surfaceArea, material, temp=0, power=5000, warp=1):
        self.tK = temp                                                  #Temperature of the object.     
        self.mass = mass                                                #Mass of the object.
        self.emissivity = emissivity                                    #Emissivity of the object. Always between 0 and 1.
        self.surfaceArea = surfaceArea                                  #Emissive surface area of the object.
        self.material = material                                        #Store the material name for some reason.
        self.specificHeat = (0.45*1000)*self.mass                       #Get the specific heat of the object in J/kg (Iron: 0.45*1000=450J/kg)
        self.power = power                                              #Joules/Second (Watts) input. This is for heating the object.
        self.warp = warp                                                #Warp Multiplier. This pertains to how KSP's warp multiplier works.

    def update(self):
        """Update the object's temperature according to it's properties."""
        #This method updates the object's temperature according to heat losses and other factors.
        self.tK -= (((self.emissivity * BOLTZMANN_CONSTANT * self.surfaceArea * (math.pow(self.tK,4) - math.pow(30+273.15,4))) / self.specificHeat) - (self.power / self.specificHeat)) * self.warp

The law used is the Stefan-Boltzmann law for calculating black-body heat losses:

Temp -= (Emissivity*Sigma*SurfaceArea*(Temp^4-Amb^4))/SpecificHeat)

This was ported from a KSP plugin for quicker debugging. Object.update() is called 50 times per second.

Would there be a solution to preventing these extreme oscillations that doesn't involve executing the code multiple times per step?

Upvotes: 0

Views: 327

Answers (2)

Robert Dodier
Robert Dodier

Reputation: 17576

I guess KSP = Kerbal Space Platform, so I gather this is a problem in game physics. If so maybe an approximation with the same qualitative behavior is sufficient. Maybe an exponential curve which starts at the initial temperature and falls to the ambient temperature is enough. Pick the decay constant by matching the heat transfer at the initial time.

Sometimes an approximation is good enough. I don't know if this is one of those situations.

Upvotes: 0

Hristo Iliev
Hristo Iliev

Reputation: 74395

Your integration scheme is bad as already hinted by @Beta and @tom10. The integration timestep is self.warp units of time, i.e. self.warp seconds since your work with physical units. This is not the way things are done. You should first convert the equation to a dimensionless format by expressing each term in some sort of computational units. For example, the Stefan-Boltzmann constant and the self.power could be measured in units, in which the constant is 1. Then you should determine the characteristic time for the object, e.g. the time by which the temperature reaches to a certain degree the equilibrium one. If there are many such objects, you should find the smallest of all characteristic times and use it as unit of measurement for the time. Then the integration timestep should be about an order of magnitude less than the characteristic time, otherwise you completely miss the correct solution to the differential equation and end up with wild oscillations.

Example of what happens now: Let's take an 1 kg iron sphere. With surface area of 3,05.10^(-3) m^2 the radiative heating/cooling power is up to 1,73.10^(-10) W/K^4. With self.power equal to 5 kW, the radiative power equates the internal one when the temperature reaches 2319 K and that's the equilibrium temperature. At low temperatures the radiative heating/cooling is negligible and with the internal heating only you end up with temperature rate of 11,1 K/s. If warp is 1000+, your first integration step results in temperature of 11100 K or more, which overshoots the equilibrium one 5 times. Now the radiative energy is orders of magnitude higher than the internal heating and it leads to huge cool-down rate - multiply by 1000+ and you end up with negative temperature. And then the cycle repeats with higher and higher absolute temperatures until you hit outside the range of the floating-point arithmetic.

Here is a hint for you: if self.power is kept constant, then the equation has an analytical solution. Find it (or use a tool like Maple or Mathematica to find it for you) and then plot the solution. See how your timestep of 1000+ units compares to the timescale of the solution, i.e. the time it takes for the system to reach an almost equilibrium state.

Upvotes: 2

Related Questions