src/FDM/Balloon/BalloonSim.cpp

   1 /*****************************************************************************
   2
   3  Module:       BalloonSim.cpp
   4  Author:       Christian Mayer
   5  Date started: 01.09.99
   6  Called by:
   7
   8  -------- Copyright (C) 1999 Christian Mayer (fgfs@christianmayer.de) --------
   9
  10  This program is free software; you can redistribute it and/or modify it under
  11  the terms of the GNU General Public License as published by the Free Software
  12  Foundation; either version 2 of the License, or (at your option) any later
  13  version.
  14
  15  This program is distributed in the hope that it will be useful, but WITHOUT
  16  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  17  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  18  details.
  19
  20  You should have received a copy of the GNU General Public License along with
  21  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  22  Place - Suite 330, Boston, MA  02111-1307, USA.
  23
  24  Further information about the GNU General Public License can also be found on
  25  the world wide web at http://www.gnu.org.
  26
  27 FUNCTIONAL DESCRIPTION
  28 ------------------------------------------------------------------------------
  29 A hot air balloon simulator
  30
  31 HISTORY
  32 ------------------------------------------------------------------------------
  33 01.09.1999 Christian Mayer      Created
  34 03.10.1999 Christian Mayer      cleaned  the code  by moveing  WeatherDatabase
  35                                 calls inside the update()
  36 *****************************************************************************/
  37
  38 /****************************************************************************/
  39 /* INCLUDES                                                                 */
  40 /****************************************************************************/
  41 #include <stdio.h>
  42 // #include <conio.h>
  43 #include <math.h>
  44
  45 #include <simgear/constants.h>
  46
  47 #include <Aircraft/aircraft.hxx>
  48 #ifdef FG_WEATHERCM
  49 # include <WeatherCM/FGLocalWeatherDatabase.h>
  50 #endif
  51
  52 #include "BalloonSim.h"
  53
  54 #include <plib/sg.h>
  55
  56 /****************************************************************************/
  57 /********************************** CODE ************************************/
  58 /****************************************************************************/
  59
  60 /****************************************************************************/
  61 /*                                                                          */
  62 /*  Constructor:                                                            */
  63 /*  Set the balloon model to some values that seem reasonable as I haven't  */
  64 /*  got much original values                                                */
  65 /*                                                                          */
  66 /****************************************************************************/
  67 balloon::balloon()
  68 {
  69     dt = 0.1;
  70     ground_level = 3400.0;
  71
  72     sgSetVec3(gravity_vector, 0.0, 0.0, -9.81);
  73     sgSetVec3(velocity, 0.0, 0.0, 0.0);
  74     sgSetVec3(position, 0.0, 0.0, 0.0);
  75     sgSetVec3(hpr,      0.0, 0.0, 0.0);
  76
  77     /************************************************************************/
  78     /* My balloon  has a  radius of  8.8 metres  as that  gives a  envelope */
  79     /* volume of about 2854 m^3  which is  about 77000 ft^3,  a very common */
  80     /* size for hot air balloons                                            */
  81     /************************************************************************/
  82
  83     balloon_envelope_area = 4.0 * (8.8 * 8.8) * SGD_PI;
  84     balloon_envelope_volume = (4.0/3.0) * (8.8 * 8.8 * 8.8) * SGD_PI;
  85
  86     wind_facing_area_of_balloon = SGD_PI * (8.8 * 8.8);
  87     wind_facing_area_of_basket = 2.0;   //guessed: 2 m^2
  88
  89     cw_envelope=0.45;                   //a sphere in this case
  90     cw_basket=0.8;
  91
  92     weight_of_total_fuel = 40.0;        //big guess
  93     weight_of_envelope = 200.0;         //big guess
  94     weight_of_basket = 40.0;            //big guess
  95     weight_of_cargo = 750.0;            //big guess
  96
  97     fuel_left=1.0;
  98     max_flow_of_fuel_per_second=10.0*1.0/3600.0; //assuming 10% of one hour of total burn time
  99     current_burner_strength = 0.0;      //the throttle
 100
 101     lambda = 0.15;                      //for plasic
 102     l_of_the_envelope = 1.0/1000.0;     //the thickness of the envelope (in m): 1mm
 103
 104     T = 273.16 + 130.6;                 //Temperature in the envelope => still at ground level
 105 }
 106
 107 void balloon::update()
 108 {
 109     /************************************************************************/
 110     /* I'm  simplifying  the  balloon by  reducing the  simulation  to  two */
 111     /* points:                                                              */
 112     /* the center of the basket (CB) and the center of the envelope (CE)    */
 113     /*                                                                      */
 114     /*                                 ce                                   */
 115     /*                                 I                                    */
 116     /*                                 I                                    */
 117     /*                                 cg (=center of gravity)              */
 118     /*                                 I                                    */
 119     /*                                 cb                                   */
 120     /*                                                                      */
 121     /* On each center  are forces acting:  gravity and  wind resitance.  CE */
 122     /* additionally got the lift  (=> I need to calculate the weight of the */
 123     /* air inside, too)                                                     */
 124     /*                                                                      */
 125     /* The weight of the air  in the envelope is  dependant of the tempera- */
 126     /* ture. This temperature is decreasing over the time that is dependant */
 127     /* of the insulation  of the envelope  material (lambda),  the gas used */
 128     /* (air) and the wind speed. For a plane surface it's for air:          */
 129     /*                                                                      */
 130     /*    alpha = 4.8 + 3.4*v   with v < 5.0 m/s                            */
 131     /*                                                                      */
 132     /* The value k that takes all of that into account is defined as:       */
 133     /*                                                                      */
 134     /*    1 / k = 1 / alpha1 + 1 / alpha2 + l / lambda                      */
 135     /*                                                                      */
 136     /* with 'l' beeing the 'length' i.e. thickness of the insulator, alpha1 */
 137     /* the air inside and alpha2 the air outside of the envelope.  So our k */
 138     /* is:                                                                  */
 139     /*                                                                      */
 140     /*    k = 1 / [1/4.8 + 1/(4.8+3.4v) + l/lambda]                         */
 141     /*                                                                      */
 142     /* The energy lost by this process is defined as:                       */
 143     /*                                                                      */
 144     /*    dQ = k * A * t * dT                                               */
 145     /*                                                                      */
 146     /* with Q being  the energy,  k that value  defined above,  A the total */
 147     /* area of the envelope,  t the time (in hours)  and dT the temperature */
 148     /* difference between the inside and the outside.                       */
 149     /* To get the temperature of the air in the inside I need the formula:  */
 150     /*                                                                      */
 151     /*    dQ = cAir * m * dT                                                */
 152     /*                                                                      */
 153     /* with cAir being the specific heat capacity(?)  of air (= 1.00 kcal / */
 154     /* kg * degree),  m the mass of the air and  dT the temperature change. */
 155     /* As the envelope is open I'm  assuming that the  same air pressure is */
 156     /* inside and  outside of  it  (practical there  should  be  a slightly */
 157     /* higher air pressure  in the inside or the envelope  would collapse). */
 158     /* So it's easy to calculate the density of the air inside:             */
 159     /*                                                                      */
 160     /*    rho = p / R * T                                                   */
 161     /*                                                                      */
 162     /* with p being  the pressure,  R the gas constant(?)  which is for air */
 163     /* 287.14 N * m / kg * K and T the absolute temperature.                */
 164     /*                                                                      */
 165     /* The value returned by this function is the displacement of the CB    */
 166     /************************************************************************/
 167
 168     /************************************************************************/
 169     /* NOTE:  This is the simplified version:  I'm assuming that  the whole */
 170     /* balloon consists only of the envelope.I will improove the simulation */
 171     /* later, but currently was my main concern to get it going...          */
 172     /************************************************************************/
 173
 174 #ifdef FG_WEATHERCM
 175     sgVec3 v;
 176
 177     FGPhysicalProperty wdbpos = WeatherDatabase->get(position);
 178
 179     //get the current wind velocity and store it in v
 180     //Point3D temp = wdbpos.Wind;
 181     //sgSetVec3(v, temp.x(), temp.y(), temp.z());
 182     sgCopyVec3(v, wdbpos.Wind );
 183
 184     sgSubVec3(v, velocity);
 185     float speed = sgLengthVec3(v);
 186
 187     // calculate the density of the gas inside
 188     double rho = wdbpos.AirPressure / (287.14 * T);
 189
 190     // calculate the mass of the air
 191     double mAir = rho * balloon_envelope_volume;
 192
 193     // loss of energy by cooling down:
 194     float k = 1.0 / (1.0/4.8 + 1.0/(4.8+3.4*speed) + l_of_the_envelope/lambda);
 195     float Q = k * balloon_envelope_area * (dt/3600.0) * (wdbpos.Temperature - T);   //(dt/3600.0) = time since last call in hours
 196
 197 #else
 198    // I realy don't think there is a solution for this without WeatherCM
 199    // but this is a hack, and it's working -- EMH
 200    double mAir = 0;
 201    float Q = 0;
 202 #endif
 203
 204     // gain of energy by heating:
 205     if (fuel_left > 0.0)        //but only with some fuel left ;-)
 206     {
 207         float fuel_burning = current_burner_strength * max_flow_of_fuel_per_second * dt * weight_of_total_fuel; //in kg
 208
 209         //convert to cubemetres (I'm wrongly assuming 'normal' conditions; but that's correct for my special case)
 210         float cube_metres_burned = fuel_burning / 2.2;  //2.2 is the density for propane
 211
 212         fuel_left -= fuel_burning / weight_of_total_fuel;
 213
 214         // get energy through burning.
 215         Q += 22250.0 * cube_metres_burned;  //22250 for propan, 29500 would be butane and if you dare: 2580 would be hydrogen...
 216     }
 217
 218     // calculate the new temperature in the inside:
 219     T += Q / (1.00 * mAir);
 220
 221     //calculate the masses of the envelope and the basket
 222     float mEnvelope = mAir + weight_of_envelope;
 223     float mBasket   = weight_of_total_fuel*fuel_left + weight_of_basket + weight_of_cargo;
 224
 225     float mTotal = mEnvelope + mBasket;
 226
 227     //calulate the forces
 228     sgVec3 fTotal, fFriction, fLift;
 229
 230     sgScaleVec3(fTotal, gravity_vector, mTotal);
 231 #ifdef FG_WEATHERCM
 232     sgScaleVec3(fFriction, v, cw_envelope * wind_facing_area_of_balloon * WeatherDatabase->getAirDensity(position) * speed / 2.0);  //wind resistance
 233     sgScaleVec3(fLift, gravity_vector, -balloon_envelope_volume * wdbpos.AirPressure / (287.14 * wdbpos.Temperature));
 234 #endif
 235
 236     sgAddVec3(fTotal, fLift);
 237     sgAddVec3(fTotal, fFriction);
 238
 239     //claculate acceleration: a = F / m
 240     sgVec3 aTotal, vTotal, dTotal;
 241
 242     sgScaleVec3(aTotal, fTotal, 1.0 / mTotal);
 243
 244     //integrate the displacement: d = 0.5 * a * dt**2 + v * dt + d
 245     sgScaleVec3(vTotal, velocity, dt);
 246     sgScaleVec3(dTotal, aTotal, 0.5*dt*dt); sgAddVec3(dTotal, vTotal);
 247
 248     //integrate the velocity to 'velocity': v = a * dt + v
 249     sgScaleVec3(vTotal, aTotal, dt); sgAddVec3(velocity, vTotal);
 250
 251     /************************************************************************/
 252     /* VERY WRONG STUFF: it's just here to get some results to start with   */
 253     /************************************************************************/
 254
 255     // care for the ground
 256     if (position[2] < (ground_level+0.001) )
 257         position[2] = ground_level;
 258
 259     //return results
 260     sgAddVec3(position, dTotal);
 261
 262     //cout << "BallonSim: T: " << (T-273.16) << " alt: " << position[2] << " ground: " << ground_level << " throttle: " << current_burner_strength << "\n";
 263 }
 264
 265 void balloon::set_burner_strength(const float bs)
 266 {
 267     if ((bs>=0.0) && (bs<=1.0))
 268     current_burner_strength = bs;
 269 }
 270
 271 void balloon::getVelocity(sgVec3 v) const
 272 {
 273     sgCopyVec3(v, velocity);
 274 }
 275
 276 void balloon::setVelocity(const sgVec3 v)
 277 {
 278     sgCopyVec3(velocity, v);
 279 }
 280
 281 void balloon::getPosition(sgVec3 v) const
 282 {
 283     sgCopyVec3(v, position);
 284 }
 285
 286 void balloon::setPosition(const sgVec3 v)
 287 {
 288     sgCopyVec3(position, v);
 289 }
 290
 291 void balloon::getHPR(sgVec3 angles) const //the balloon isn't allways exactly vertical
 292 {
 293     sgCopyVec3(angles, hpr);
 294 }
 295
 296 void balloon::setHPR(const sgVec3 angles)  //the balloon isn't allways exactly vertical
 297 {
 298     sgCopyVec3(hpr, angles);
 299 }
 300
 301 void balloon::setGroundLevel(const float altitude)
 302 {
 303     ground_level = altitude;
 304 }
 305
 306 float balloon::getTemperature(void) const
 307 {
 308     return T;
 309 }
 310
 311 float balloon::getFuelLeft(void) const
 312 {
 313     return fuel_left;
 314 }
 315
 316 /*
 317 void main(void)
 318 {
 319     bool burner = true;
 320     balloon bal;
 321     bool exit = false;
 322     sgVec3 pos={0.0, 0.0, 0.0};
 323     char c;
 324     float acc_dt = 0.0;
 325     float alt;
 326
 327 bool hysteresis = false; // moving up
 328     for (;!exit;)
 329     {
 330         for (int i=0; i<100; i++)
 331         {
 332             bal.update(0.1); acc_dt += 0.1;
 333             bal.getPosition(pos);
 334             alt = pos[2];
 335
 336             if (alt > 3010)
 337             {
 338                 hysteresis = true;
 339                 burner = false;
 340             }
 341             if ((alt < 2990) && (hysteresis == true))
 342             {
 343                 hysteresis = false;
 344                 burner = true;
 345             }
 346             if ((bal.getTemperature()-273.16)>250.0)
 347                 burner = false; //emergency
 348         }
 349
 350         // toogle burner
 351         c = getch();
 352         if (c==' ')
 353             burner=!burner;
 354         //if (c=='s')
 355         //    show=!show;
 356
 357         //printf("Position: (%f/%f/%f), dP: (%f/%f/%f), burner: ", pos[0], pos[1], pos[2], dp[0], dp[1], dp[2]);
 358         printf("%f         \t%f         \t%f         \t%f\n", acc_dt/60.0, bal.getTemperature()-273.16, pos[2], bal.getFuelLeft());
 359         if (burner==true)
 360         {
 361             //printf("on\n");
 362             bal.set_burner_strength(1.0);
 363         }
 364         else
 365         {
 366             //printf("off\n");
 367             bal.set_burner_strength(0.0);
 368         }
 369
 370     }
 371 }
 372 */