358 lines
14 KiB
C++
358 lines
14 KiB
C++
|
/*****************************************************************************
|
||
|
|
|
||
|
|
Module: BalloonSim.cpp
|
||
|
|
Author: Christian Mayer
|
||
|
|
Date started: 01.09.99
|
||
|
|
Called by:
|
||
|
|
|
||
|
|
-------- Copyright (C) 1999 Christian Mayer (fgfs@christianmayer.de) --------
|
||
|
|
|
||
|
|
This program is free software; you can redistribute it and/or modify it under
|
||
|
|
the terms of the GNU General Public License as published by the Free Software
|
||
|
|
Foundation; either version 2 of the License, or (at your option) any later
|
||
|
|
version.
|
||
|
|
|
||
|
|
This program is distributed in the hope that it will be useful, but WITHOUT
|
||
|
|
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
||
|
|
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
|
||
|
|
details.
|
||
|
|
|
||
|
|
You should have received a copy of the GNU General Public License along with
|
||
|
|
this program; if not, write to the Free Software Foundation, Inc., 59 Temple
|
||
|
|
Place - Suite 330, Boston, MA 02111-1307, USA.
|
||
|
|
|
||
|
|
Further information about the GNU General Public License can also be found on
|
||
|
|
the world wide web at http://www.gnu.org.
|
||
|
|
|
||
|
|
FUNCTIONAL DESCRIPTION
|
||
|
|
------------------------------------------------------------------------------
|
||
|
|
A hot air balloon simulator
|
||
|
|
|
||
|
|
HISTORY
|
||
|
|
------------------------------------------------------------------------------
|
||
|
|
01.09.1999 Christian Mayer Created
|
||
|
|
03.10.1999 Christian Mayer cleaned the code by moveing WeatherDatabase
|
||
|
|
calls inside the update()
|
||
|
|
*****************************************************************************/
|
||
|
|
|
||
|
|
/****************************************************************************/
|
||
|
|
/* INCLUDES */
|
||
|
|
/****************************************************************************/
|
||
|
|
#include <stdio.h>
|
||
|
|
// #include <conio.h>
|
||
|
|
#include <math.h>
|
||
|
|
|
||
|
|
#include <Aircraft/aircraft.hxx>
|
||
|
|
#include <Include/fg_constants.h>
|
||
|
|
#include <Main/options.hxx>
|
||
|
|
#include <WeatherCM/FGLocalWeatherDatabase.h>
|
||
|
|
|
||
|
|
#include "BalloonSim.h"
|
||
|
|
|
||
|
|
#include "sg.h"
|
||
|
|
|
||
|
|
/****************************************************************************/
|
||
|
|
/********************************** CODE ************************************/
|
||
|
|
/****************************************************************************/
|
||
|
|
|
||
|
|
/****************************************************************************/
|
||
|
|
/* */
|
||
|
|
/* Constructor: */
|
||
|
|
/* Set the balloon model to some values that seem reasonable as I haven't */
|
||
|
|
/* got much original values */
|
||
|
|
/* */
|
||
|
|
/****************************************************************************/
|
||
|
|
balloon::balloon()
|
||
|
|
{
|
||
|
|
dt = 0.1;
|
||
|
|
ground_level = 3400.0;
|
||
|
|
|
||
|
|
sgSetVec3(gravity_vector, 0.0, 0.0, -9.81);
|
||
|
|
sgSetVec3(velocity, 0.0, 0.0, 0.0);
|
||
|
|
sgSetVec3(position, 0.0, 0.0, 0.0);
|
||
|
|
sgSetVec3(hpr, 0.0, 0.0, 0.0);
|
||
|
|
|
||
|
|
/************************************************************************/
|
||
|
|
/* My balloon has a radius of 8.8 metres as that gives a envelope */
|
||
|
|
/* volume of about 2854 m^3 which is about 77000 ft^3, a very common */
|
||
|
|
/* size for hot air balloons */
|
||
|
|
/************************************************************************/
|
||
|
|
|
||
|
|
balloon_envelope_area = 4.0 * (8.8 * 8.8) * FG_PI;
|
||
|
|
balloon_envelope_volume = (4.0/3.0) * (8.8 * 8.8 * 8.8) * FG_PI;
|
||
|
|
|
||
|
|
wind_facing_area_of_balloon = FG_PI * (8.8 * 8.8);
|
||
|
|
wind_facing_area_of_basket = 2.0; //guessed: 2 m^2
|
||
|
|
|
||
|
|
cw_envelope=0.45; //a sphere in this case
|
||
|
|
cw_basket=0.8;
|
||
|
|
|
||
|
|
weight_of_total_fuel = 40.0; //big guess
|
||
|
|
weight_of_envelope = 200.0; //big guess
|
||
|
|
weight_of_basket = 40.0; //big guess
|
||
|
|
weight_of_cargo = 750.0; //big guess
|
||
|
|
|
||
|
|
fuel_left=1.0;
|
||
|
|
max_flow_of_fuel_per_second=10.0*1.0/3600.0; //assuming 10% of one hour of total burn time
|
||
|
|
current_burner_strength = 0.0; //the throttle
|
||
|
|
|
||
|
|
lambda = 0.15; //for plasic
|
||
|
|
l_of_the_envelope = 1.0/1000.0; //the thickness of the envelope (in m): 1mm
|
||
|
|
|
||
|
|
T = 273.16 + 130.6; //Temperature in the envelope => still at ground level
|
||
|
|
}
|
||
|
|
|
||
|
|
void balloon::update()
|
||
|
|
{
|
||
|
|
/************************************************************************/
|
||
|
|
/* I'm simplifying the balloon by reducing the simulation to two */
|
||
|
|
/* points: */
|
||
|
|
/* the center of the basket (CB) and the center of the envelope (CE) */
|
||
|
|
/* */
|
||
|
|
/* ce */
|
||
|
|
/* I */
|
||
|
|
/* I */
|
||
|
|
/* cg (=center of gravity) */
|
||
|
|
/* I */
|
||
|
|
/* cb */
|
||
|
|
/* */
|
||
|
|
/* On each center are forces acting: gravity and wind resitance. CE */
|
||
|
|
/* additionally got the lift (=> I need to calculate the weight of the */
|
||
|
|
/* air inside, too) */
|
||
|
|
/* */
|
||
|
|
/* The weight of the air in the envelope is dependant of the tempera- */
|
||
|
|
/* ture. This temperature is decreasing over the time that is dependant */
|
||
|
|
/* of the insulation of the envelope material (lambda), the gas used */
|
||
|
|
/* (air) and the wind speed. For a plane surface it's for air: */
|
||
|
|
/* */
|
||
|
|
/* alpha = 4.8 + 3.4*v with v < 5.0 m/s */
|
||
|
|
/* */
|
||
|
|
/* The value k that takes all of that into account is defined as: */
|
||
|
|
/* */
|
||
|
|
/* 1 / k = 1 / alpha1 + 1 / alpha2 + l / lambda */
|
||
|
|
/* */
|
||
|
|
/* with 'l' beeing the 'length' i.e. thickness of the insulator, alpha1 */
|
||
|
|
/* the air inside and alpha2 the air outside of the envelope. So our k */
|
||
|
|
/* is: */
|
||
|
|
/* */
|
||
|
|
/* k = 1 / [1/4.8 + 1/(4.8+3.4v) + l/lambda] */
|
||
|
|
/* */
|
||
|
|
/* The energy lost by this process is defined as: */
|
||
|
|
/* */
|
||
|
|
/* dQ = k * A * t * dT */
|
||
|
|
/* */
|
||
|
|
/* with Q being the energy, k that value defined above, A the total */
|
||
|
|
/* area of the envelope, t the time (in hours) and dT the temperature */
|
||
|
|
/* difference between the inside and the outside. */
|
||
|
|
/* To get the temperature of the air in the inside I need the formula: */
|
||
|
|
/* */
|
||
|
|
/* dQ = cAir * m * dT */
|
||
|
|
/* */
|
||
|
|
/* with cAir being the specific heat capacity(?) of air (= 1.00 kcal / */
|
||
|
|
/* kg * degree), m the mass of the air and dT the temperature change. */
|
||
|
|
/* As the envelope is open I'm assuming that the same air pressure is */
|
||
|
|
/* inside and outside of it (practical there should be a slightly */
|
||
|
|
/* higher air pressure in the inside or the envelope would collapse). */
|
||
|
|
/* So it's easy to calculate the density of the air inside: */
|
||
|
|
/* */
|
||
|
|
/* rho = p / R * T */
|
||
|
|
/* */
|
||
|
|
/* with p being the pressure, R the gas constant(?) which is for air */
|
||
|
|
/* 287.14 N * m / kg * K and T the absolute temperature. */
|
||
|
|
/* */
|
||
|
|
/* The value returned by this function is the displacement of the CB */
|
||
|
|
/************************************************************************/
|
||
|
|
|
||
|
|
/************************************************************************/
|
||
|
|
/* NOTE: This is the simplified version: I'm assuming that the whole */
|
||
|
|
/* balloon consists only of the envelope.I will improove the simulation */
|
||
|
|
/* later, but currently was my main concern to get it going... */
|
||
|
|
/************************************************************************/
|
||
|
|
|
||
|
|
sgVec3 v;
|
||
|
|
|
||
|
|
//get the current wind velocity and store it in v
|
||
|
|
Point3D temp = WeatherDatabase->get(position).Wind;
|
||
|
|
sgSetVec3(v, temp.x(), temp.y(), temp.z());
|
||
|
|
|
||
|
|
sgSubVec3(v, velocity);
|
||
|
|
float speed = sgLengthVec3(v);
|
||
|
|
|
||
|
|
// calculate the density of the gas inside
|
||
|
|
double rho = WeatherDatabase->get(position).AirPressure / (287.14 * T);
|
||
|
|
|
||
|
|
// calculate the mass of the air
|
||
|
|
double mAir = rho * balloon_envelope_volume;
|
||
|
|
|
||
|
|
// loss of energy by cooling down:
|
||
|
|
float k = 1.0 / (1.0/4.8 + 1.0/(4.8+3.4*speed) + l_of_the_envelope/lambda);
|
||
|
|
float Q = k * balloon_envelope_area * (dt/3600.0) * (WeatherDatabase->get(position).Temperature - T); //(dt/3600.0) = time since last call in hours
|
||
|
|
|
||
|
|
// gain of energy by heating:
|
||
|
|
if (fuel_left > 0.0) //but only with some fuel left ;-)
|
||
|
|
{
|
||
|
|
float fuel_burning = current_burner_strength * max_flow_of_fuel_per_second * dt * weight_of_total_fuel; //in kg
|
||
|
|
|
||
|
|
//convert to cubemetres (I'm wrongly assuming 'normal' conditions; but that's correct for my special case)
|
||
|
|
float cube_metres_burned = fuel_burning / 2.2; //2.2 is the density for propane
|
||
|
|
|
||
|
|
fuel_left -= fuel_burning / weight_of_total_fuel;
|
||
|
|
|
||
|
|
// get energy through burning.
|
||
|
|
Q += 22250.0 * cube_metres_burned; //22250 for propan, 29500 would be butane and if you dare: 2580 would be hydrogen...
|
||
|
|
}
|
||
|
|
|
||
|
|
// calculate the new temperature in the inside:
|
||
|
|
T += Q / (1.00 * mAir);
|
||
|
|
|
||
|
|
//calculate the masses of the envelope and the basket
|
||
|
|
float mEnvelope = mAir + weight_of_envelope;
|
||
|
|
float mBasket = weight_of_total_fuel*fuel_left + weight_of_basket + weight_of_cargo;
|
||
|
|
|
||
|
|
float mTotal = mEnvelope + mBasket;
|
||
|
|
|
||
|
|
//calulate the forces
|
||
|
|
sgVec3 fTotal, fFriction, fLift;
|
||
|
|
|
||
|
|
sgScaleVec3(fTotal, gravity_vector, mTotal);
|
||
|
|
sgScaleVec3(fFriction, v, cw_envelope * wind_facing_area_of_balloon * WeatherDatabase->getAirDensity(position) * speed / 2.0); //wind resistance
|
||
|
|
sgScaleVec3(fLift, gravity_vector, -balloon_envelope_volume * WeatherDatabase->get(position).AirPressure / (287.14 * WeatherDatabase->get(position).Temperature));
|
||
|
|
|
||
|
|
sgAddVec3(fTotal, fLift);
|
||
|
|
sgAddVec3(fTotal, fFriction);
|
||
|
|
|
||
|
|
//claculate acceleration: a = F / m
|
||
|
|
sgVec3 aTotal, vTotal, dTotal;
|
||
|
|
|
||
|
|
sgScaleVec3(aTotal, fTotal, 1.0 / mTotal);
|
||
|
|
|
||
|
|
//integrate the displacement: d = 0.5 * a * dt**2 + v * dt + d
|
||
|
|
sgScaleVec3(vTotal, velocity, dt);
|
||
|
|
sgScaleVec3(dTotal, aTotal, 0.5*dt*dt); sgAddVec3(dTotal, vTotal);
|
||
|
|
|
||
|
|
//integrate the velocity to 'velocity': v = a * dt + v
|
||
|
|
sgScaleVec3(vTotal, aTotal, dt); sgAddVec3(velocity, vTotal);
|
||
|
|
|
||
|
|
/************************************************************************/
|
||
|
|
/* VERY WRONG STUFF: it's just here to get some results to start with */
|
||
|
|
/************************************************************************/
|
||
|
|
|
||
|
|
// care for the ground
|
||
|
|
if (position[2] < (ground_level+0.001) )
|
||
|
|
velocity[2] = ground_level;
|
||
|
|
|
||
|
|
//return results
|
||
|
|
sgAddVec3(position, dTotal);
|
||
|
|
|
||
|
|
//cout << "BallonSim: T: " << (T-273.16) << " alt: " << position[2] << " ground: " << ground_level << " throttle: " << current_burner_strength << "\n";
|
||
|
|
}
|
||
|
|
|
||
|
|
void balloon::set_burner_strength(const float bs)
|
||
|
|
{
|
||
|
|
if ((bs>=0.0) && (bs<=1.0))
|
||
|
|
current_burner_strength = bs;
|
||
|
|
}
|
||
|
|
|
||
|
|
void balloon::getVelocity(sgVec3 v) const
|
||
|
|
{
|
||
|
|
sgCopyVec3(v, velocity);
|
||
|
|
}
|
||
|
|
|
||
|
|
void balloon::setVelocity(const sgVec3 v)
|
||
|
|
{
|
||
|
|
sgCopyVec3(velocity, v);
|
||
|
|
}
|
||
|
|
|
||
|
|
void balloon::getPosition(sgVec3 v) const
|
||
|
|
{
|
||
|
|
sgCopyVec3(v, position);
|
||
|
|
}
|
||
|
|
|
||
|
|
void balloon::setPosition(const sgVec3 v)
|
||
|
|
{
|
||
|
|
sgCopyVec3(position, v);
|
||
|
|
}
|
||
|
|
|
||
|
|
void balloon::getHPR(sgVec3 angles) const //the balloon isn't allways exactly vertical
|
||
|
|
{
|
||
|
|
sgCopyVec3(angles, hpr);
|
||
|
|
}
|
||
|
|
|
||
|
|
void balloon::setHPR(const sgVec3 angles) //the balloon isn't allways exactly vertical
|
||
|
|
{
|
||
|
|
sgCopyVec3(hpr, angles);
|
||
|
|
}
|
||
|
|
|
||
|
|
void balloon::setGroundLevel(const float altitude)
|
||
|
|
{
|
||
|
|
ground_level = altitude;
|
||
|
|
}
|
||
|
|
|
||
|
|
float balloon::getTemperature(void) const
|
||
|
|
{
|
||
|
|
return T;
|
||
|
|
}
|
||
|
|
|
||
|
|
float balloon::getFuelLeft(void) const
|
||
|
|
{
|
||
|
|
return fuel_left;
|
||
|
|
}
|
||
|
|
|
||
|
|
/*
|
||
|
|
void main(void)
|
||
|
|
{
|
||
|
|
bool burner = true;
|
||
|
|
balloon bal;
|
||
|
|
bool exit = false;
|
||
|
|
sgVec3 pos={0.0, 0.0, 0.0};
|
||
|
|
char c;
|
||
|
|
float acc_dt = 0.0;
|
||
|
|
float alt;
|
||
|
|
|
||
|
|
bool hysteresis = false; // moving up
|
||
|
|
for (;!exit;)
|
||
|
|
{
|
||
|
|
for (int i=0; i<100; i++)
|
||
|
|
{
|
||
|
|
bal.update(0.1); acc_dt += 0.1;
|
||
|
|
bal.getPosition(pos);
|
||
|
|
alt = pos[2];
|
||
|
|
|
||
|
|
if (alt > 3010)
|
||
|
|
{
|
||
|
|
hysteresis = true;
|
||
|
|
burner = false;
|
||
|
|
}
|
||
|
|
if ((alt < 2990) && (hysteresis == true))
|
||
|
|
{
|
||
|
|
hysteresis = false;
|
||
|
|
burner = true;
|
||
|
|
}
|
||
|
|
if ((bal.getTemperature()-273.16)>250.0)
|
||
|
|
burner = false; //emergency
|
||
|
|
}
|
||
|
|
|
||
|
|
// toogle burner
|
||
|
|
c = getch();
|
||
|
|
if (c==' ')
|
||
|
|
burner=!burner;
|
||
|
|
//if (c=='s')
|
||
|
|
// show=!show;
|
||
|
|
|
||
|
|
//printf("Position: (%f/%f/%f), dP: (%f/%f/%f), burner: ", pos[0], pos[1], pos[2], dp[0], dp[1], dp[2]);
|
||
|
|
printf("%f \t%f \t%f \t%f\n", acc_dt/60.0, bal.getTemperature()-273.16, pos[2], bal.getFuelLeft());
|
||
|
|
if (burner==true)
|
||
|
|
{
|
||
|
|
//printf("on\n");
|
||
|
|
bal.set_burner_strength(1.0);
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
//printf("off\n");
|
||
|
|
bal.set_burner_strength(0.0);
|
||
|
|
}
|
||
|
|
|
||
|
|
}
|
||
|
|
}
|
||
|
|
*/
|