fggascell.cpp

来自「6 DOF Missle Simulation」· C++ 代码 · 共 873 行 · 第 1/3 页

/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 Header:       FGGasCell.h
 Author:       Anders Gidenstam
 Date started: 01/21/2006

 ----- Copyright (C) 2006 - 2008  Anders Gidenstam (anders(at)gidenstam.org) --

 This program is free software; you can redistribute it and/or modify it under
 the terms of the GNU Lesser 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 Lesser General Public License for more
 details.

 You should have received a copy of the GNU Lesser 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 Lesser General Public License can also be found on
 the world wide web at http://www.gnu.org.

FUNCTIONAL DESCRIPTION
--------------------------------------------------------------------------------
See header file.

HISTORY
--------------------------------------------------------------------------------
01/21/2006  AG   Created

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
INCLUDES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

#include <FGFDMExec.h>
#include <models/FGAuxiliary.h>
#include <models/FGAtmosphere.h>
#include <models/FGInertial.h>
#include <models/FGMassBalance.h>
#include "FGGasCell.h"

using std::cerr;
using std::endl;
using std::cout;

namespace JSBSim {

static const char *IdSrc = "$Id$";
static const char *IdHdr = ID_GASCELL;

/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CLASS IMPLEMENTATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/* Constants. */
const double FGGasCell::R = 3.4071;              // [lbs ft/(mol Rankine)]
const double FGGasCell::M_air = 0.0019186;       // [slug/mol]
const double FGGasCell::M_hydrogen = 0.00013841; // [slug/mol]
const double FGGasCell::M_helium = 0.00027409;   // [slug/mol]

FGGasCell::FGGasCell(FGFDMExec* exec, Element* el, int num) : FGForce(exec)
{
  string token;
  Element* element;

  Auxiliary = exec->GetAuxiliary();
  Atmosphere = exec->GetAtmosphere();
  PropertyManager = exec->GetPropertyManager();
  Inertial = exec->GetInertial();
  MassBalance = exec->GetMassBalance();

  gasCellJ = FGMatrix33();
  gasCellM = FGColumnVector3();

  Buoyancy = MaxVolume = MaxOverpressure = Temperature = Pressure =
    Contents = Volume = dVolumeIdeal = 0.0;
  Xradius = Yradius = Zradius = Xwidth = Ywidth = Zwidth = 0.0;
  ValveCoefficient = ValveOpen = 0.0;
  CellNum = num;

  // NOTE: In the local system X points north, Y points east and Z points down.
  SetTransformType(FGForce::tLocalBody);

  type = el->GetAttributeValue("type");
  if      (type == "HYDROGEN") Type = ttHYDROGEN;
  else if (type == "HELIUM")   Type = ttHELIUM;
  else if (type == "AIR")      Type = ttAIR;
  else                         Type = ttUNKNOWN;

  element = el->FindElement("location");
  if (element) {
    vXYZ = element->FindElementTripletConvertTo("IN");
  } else {
    cerr << "Fatal Error: No location found for this gas cell." << endl;
    exit(-1);
  }
  if ((el->FindElement("x_radius") || el->FindElement("x_width")) &&
      (el->FindElement("y_radius") || el->FindElement("y_width")) &&
      (el->FindElement("z_radius") || el->FindElement("z_width"))) {

    if (el->FindElement("x_radius")) {
      Xradius = el->FindElementValueAsNumberConvertTo("x_radius", "FT");
    }
    if (el->FindElement("y_radius")) {
      Yradius = el->FindElementValueAsNumberConvertTo("y_radius", "FT");
    }
    if (el->FindElement("z_radius")) {
      Zradius = el->FindElementValueAsNumberConvertTo("z_radius", "FT");
    }

    if (el->FindElement("x_width")) {
      Xwidth = el->FindElementValueAsNumberConvertTo("x_width", "FT");
    }
    if (el->FindElement("y_width")) {
      Ywidth = el->FindElementValueAsNumberConvertTo("y_width", "FT");
    }
    if (el->FindElement("z_width")) {
      Zwidth = el->FindElementValueAsNumberConvertTo("z_width", "FT");
    }

    // The volume is a (potentially) extruded ellipsoid.
    // However, currently only a few combinations of radius and width are
    // fully supported.
    if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
        (Xwidth  == 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
      // Ellipsoid volume.
      MaxVolume = 4.0  * M_PI * Xradius * Yradius * Zradius / 3.0;
    } else if  ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
                (Xwidth  != 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
      // Cylindrical volume.
      MaxVolume = M_PI * Yradius * Zradius * Xwidth;
    } else {
      cerr << "Warning: Unsupported gas cell shape." << endl;
      MaxVolume = 
        (4.0  * M_PI * Xradius * Yradius * Zradius / 3.0 +
         M_PI * Yradius * Zradius * Xwidth +
         M_PI * Xradius * Zradius * Ywidth +
         M_PI * Xradius * Yradius * Zwidth +
         2.0  * Xradius * Ywidth * Zwidth +
         2.0  * Yradius * Xwidth * Zwidth +
         2.0  * Zradius * Xwidth * Ywidth +
         Xwidth * Ywidth * Zwidth);
    }
  } else {
    cerr << "Fatal Error: Gas cell shape must be given." << endl;
    exit(-1);
  }
  if (el->FindElement("max_overpressure")) {
    MaxOverpressure = el->FindElementValueAsNumberConvertTo("max_overpressure",
                                                            "LBS/FT2");
  }
  if (el->FindElement("fullness")) {
    const double Fullness = el->FindElementValueAsNumber("fullness");
    if (0 <= Fullness) { 
      Volume = Fullness * MaxVolume; 
    } else {
      cerr << "Warning: Invalid initial gas cell fullness value." << endl;
    }
  }  
  if (el->FindElement("valve_coefficient")) {
    ValveCoefficient =
      el->FindElementValueAsNumberConvertTo("valve_coefficient",
                                            "FT4*SEC/SLUG");
    ValveCoefficient = max(ValveCoefficient, 0.0);
  }

  // Initialize state
  SetLocation(vXYZ);

  if (Temperature == 0.0) {
    Temperature = Atmosphere->GetTemperature();
  }
  if (Pressure == 0.0) {
    Pressure = Atmosphere->GetPressure();
  }
  if (Volume != 0.0) {
    // Calculate initial gas content.
    Contents = Pressure * Volume / (R * Temperature);
    
    // Clip to max allowed value.
    const double IdealPressure = Contents * R * Temperature / MaxVolume;
    if (IdealPressure > Pressure + MaxOverpressure) {
      Contents = (Pressure + MaxOverpressure) * MaxVolume / (R * Temperature);
      Pressure = Pressure + MaxOverpressure;
    } else {
      Pressure = max(IdealPressure, Pressure);
    }
  } else {
    // Calculate initial gas content.
    Contents = Pressure * MaxVolume / (R * Temperature);
  }

  Volume = Contents * R * Temperature / Pressure;
  Mass = Contents * M_gas();

  // Bind relevant properties
  char property_name[80];
  snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/max_volume-ft3",
           CellNum);
  PropertyManager->Tie( property_name, &MaxVolume );
  PropertyManager->SetWritable( property_name, false );
  snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/temp-R",
           CellNum);
  PropertyManager->Tie( property_name, &Temperature );
  snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/pressure-psf",
           CellNum);
  PropertyManager->Tie( property_name, &Pressure );
  snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/volume-ft3",
           CellNum);
  PropertyManager->Tie( property_name, &Volume );
  snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/buoyancy-lbs",
           CellNum);
  PropertyManager->Tie( property_name, &Buoyancy );
  snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/contents-mol",
           CellNum);
  PropertyManager->Tie( property_name, &Contents );
  snprintf(property_name, 80, "buoyant_forces/gas-cell[%d]/valve_open",
           CellNum);
  PropertyManager->Tie( property_name, &ValveOpen );

  Debug(0);

  // Read heat transfer coefficients
  if (Element* heat = el->FindElement("heat")) {
    Element* function_element = heat->FindElement("function");
    while (function_element) {
      HeatTransferCoeff.push_back(new FGFunction(PropertyManager,
                                                 function_element));
      function_element = heat->FindNextElement("function");
    }
  }

  // Load ballonets if there are any
  if (Element* ballonet_element = el->FindElement("ballonet")) {
    while (ballonet_element) {
      Ballonet.push_back(new FGBallonet(exec,
                                        ballonet_element,
                                        Ballonet.size(),
                                        this));
      ballonet_element = el->FindNextElement("ballonet");
    }
  }

}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

FGGasCell::~FGGasCell()
{
  unsigned int i;

  for (i = 0; i < HeatTransferCoeff.size(); i++) delete HeatTransferCoeff[i];
  HeatTransferCoeff.clear();

  for (i = 0; i < Ballonet.size(); i++) delete Ballonet[i];
  Ballonet.clear();

  Debug(1);
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

void FGGasCell::Calculate(double dt)
{
  const double AirTemperature = Atmosphere->GetTemperature();  // [Rankine]
  const double AirPressure    = Atmosphere->GetPressure();     // [lbs/ft瞉
  const double AirDensity     = Atmosphere->GetDensity();      // [slug/ft砞
  const double g = Inertial->gravity();                        // [lbs/slug]

  const double OldTemperature = Temperature;
  const double OldPressure    = Pressure;
  unsigned int i;
  const unsigned int no_ballonets = Ballonet.size();

  //-- Read ballonet state --
  // NOTE: This model might need a more proper integration technique. 
  double BallonetsVolume = 0.0;
  double BallonetsHeatFlow = 0.0;
  for (i = 0; i < no_ballonets; i++) {
    BallonetsVolume   += Ballonet[i]->GetVolume();
    BallonetsHeatFlow += Ballonet[i]->GetHeatFlow();
  }

  //-- Gas temperature --

  if (HeatTransferCoeff.size() > 0) {
    // The model is based on the ideal gas law.
    // However, it does look a bit fishy. Please verify.
    //   dT/dt = dU / (Cv n R)