fggascell.cpp

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

      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: Ballonet 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 ballonet fullness value." << endl;
    }
  }  
  if (el->FindElement("valve_coefficient")) {
    ValveCoefficient =
      el->FindElementValueAsNumberConvertTo("valve_coefficient",
                                            "FT4*SEC/SLUG");
    ValveCoefficient = max(ValveCoefficient, 0.0);
  }

  // Initialize state
  if (Temperature == 0.0) {
    Temperature = Parent->GetTemperature();
  }
  if (Pressure == 0.0) {
    Pressure = Parent->GetPressure();
  }
  if (Volume != 0.0) {
    // Calculate initial air 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 air content.
    Contents = Pressure * MaxVolume / (R * Temperature);
  }

  Volume = Contents * R * Temperature / Pressure;

  // Bind relevant properties
  char property_name[80];
  snprintf(property_name, 80,
           "buoyant_forces/gas-cell[%d]/ballonet[%d]/max_volume-ft3",
           Parent->GetIndex(),
           CellNum);
  PropertyManager->Tie( property_name, &MaxVolume );
  PropertyManager->SetWritable( property_name, false );
  snprintf(property_name, 80,
           "buoyant_forces/gas-cell[%d]/ballonet[%d]/temp-R",
           Parent->GetIndex(),
           CellNum);
  PropertyManager->Tie( property_name, &Temperature );
  snprintf(property_name, 80,
           "buoyant_forces/gas-cell[%d]/ballonet[%d]/pressure-psf",
           Parent->GetIndex(),
           CellNum);
  PropertyManager->Tie( property_name, &Pressure );
  snprintf(property_name, 80,
           "buoyant_forces/gas-cell[%d]/ballonet[%d]/volume-ft3",
           Parent->GetIndex(),
           CellNum);
  PropertyManager->Tie( property_name, &Volume );
  snprintf(property_name, 80,
           "buoyant_forces/gas-cell[%d]/ballonet[%d]/contents-mol",
           Parent->GetIndex(),
           CellNum);
  PropertyManager->Tie( property_name, &Contents );
  snprintf(property_name, 80,
           "buoyant_forces/gas-cell[%d]/ballonet[%d]/valve_open",
           Parent->GetIndex(),
           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");
    }
  }
  // Read blower input function
  if (Element* blower = el->FindElement("blower_input")) {
    Element* function_element = blower->FindElement("function");
    BlowerInput = new FGFunction(PropertyManager,
                                 function_element);
  }
}

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

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

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

  delete BlowerInput;
  BlowerInput = NULL;

  Debug(1);
}

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

void FGBallonet::Calculate(double dt)
{
  const double ParentPressure = Parent->GetPressure();         // [lbs/ft瞉
  const double AirPressure    = Atmosphere->GetPressure();     // [lbs/ft瞉

  const double OldTemperature = Temperature;
  const double OldPressure    = Pressure;
  unsigned int i;

  //-- Gas temperature --

  // The model is based on the ideal gas law.
  // However, it does look a bit fishy. Please verify.
  //   dT/dt = dU / (Cv n R)
  dU = 0.0;
  for (i = 0; i < HeatTransferCoeff.size(); i++) {
    dU += HeatTransferCoeff[i]->GetValue();
  }
  // dt is already accounted for in dVolumeIdeal.
  if (Contents > 0) {
    Temperature +=
      (dU * dt - Pressure * dVolumeIdeal) / (Cv_air * Contents * R);
  } else {
    Temperature = Parent->GetTemperature();
  }

  //-- Pressure --
  const double IdealPressure = Contents * R * Temperature / MaxVolume;
  // The pressure is at least that of the parent gas cell.
  Pressure = max(IdealPressure, ParentPressure);

  //-- Blower input --
  if (BlowerInput) {
    const double AddedVolume = BlowerInput->GetValue() * dt;
    if (AddedVolume > 0.0) {
      Contents += Pressure * AddedVolume / (R * Temperature);
    }
  }

  //-- Pressure relief and manual valving --
  // FIXME: Presently the effect of valving is computed using
  //        an ad hoc formula which might not be a good representation
  //        of reality.
  if ((ValveCoefficient > 0.0) &&
      ((ValveOpen > 0.0) || (Pressure > AirPressure + MaxOverpressure))) {
    const double DeltaPressure = Pressure - AirPressure;
    const double VolumeValved =
      ((Pressure > AirPressure + MaxOverpressure) ? 1.0 : ValveOpen) *
      ValveCoefficient * DeltaPressure * dt;
    // FIXME: Too small values of Contents sometimes leads to NaN.
    //        Currently the minimum is restricted to a safe value.
    Contents =
      max(1.0, Contents - Pressure * VolumeValved / (R * Temperature));
  }

  //-- Volume --
  Volume = Contents * R * Temperature / Pressure;
  dVolumeIdeal =
    Contents * R * (Temperature / Pressure - OldTemperature / OldPressure);

  // Compute the inertia of the ballonet.
  // Consider the ballonet as a shape of uniform density.
  // FIXME: If the ballonet isn't ellipsoid or cylindrical the inertia will
  //        be wrong.
  ballonetJ = FGMatrix33();
  const double mass = Contents * M_air;
  double Ixx, Iyy, Izz;
  if ((Xradius != 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
      (Xwidth  == 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
    // Ellipsoid volume.
    Ixx = (1.0 / 5.0) * mass * (Yradius*Yradius + Zradius*Zradius);
    Iyy = (1.0 / 5.0) * mass * (Xradius*Xradius + Zradius*Zradius);
    Izz = (1.0 / 5.0) * mass * (Xradius*Xradius + Yradius*Yradius);     
  } else if  ((Xradius == 0.0) && (Yradius != 0.0) && (Zradius != 0.0) &&
              (Xwidth  != 0.0) && (Ywidth  == 0.0) && (Zwidth  == 0.0)) {
    // Cylindrical volume (might not be valid with an elliptical cross-section).
    Ixx = (1.0 / 2.0) * mass * Yradius * Zradius;
    Iyy =
      (1.0 / 4.0) * mass * Yradius * Zradius +
      (1.0 / 12.0) * mass * Xwidth * Xwidth;
    Izz =
      (1.0 / 4.0) * mass * Yradius * Zradius +
      (1.0 / 12.0) * mass * Xwidth * Xwidth;
  } else {
    // Not supported. Revert to pointmass model.
    Ixx = Iyy = Izz = 0.0;
  }
  // The volume is symmetric, so Ixy = Ixz = Iyz = 0.
  ballonetJ(1,1) = Ixx;
  ballonetJ(2,2) = Iyy;
  ballonetJ(3,3) = Izz;
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//    The bitmasked value choices are as follows:
//    unset: In this case (the default) JSBSim would only print
//       out the normally expected messages, essentially echoing
//       the config files as they are read. If the environment
//       variable is not set, debug_lvl is set to 1 internally
//    0: This requests JSBSim not to output any messages
//       whatsoever.
//    1: This value explicity requests the normal JSBSim
//       startup messages
//    2: This value asks for a message to be printed out when
//       a class is instantiated
//    4: When this value is set, a message is displayed when a
//       FGModel object executes its Run() method
//    8: When this value is set, various runtime state variables
//       are printed out periodically
//    16: When set various parameters are sanity checked and
//       a message is printed out when they go out of bounds

void FGBallonet::Debug(int from)
{
  if (debug_lvl <= 0) return;

  if (debug_lvl & 1) { // Standard console startup message output
    if (from == 0) { // Constructor
      cout << "      Ballonet holds " << Contents << " mol air" << endl;
      cout << "        Location (X, Y, Z) (in.): " << vXYZ(eX) << ", " <<
        vXYZ(eY) << ", " << vXYZ(eZ) << endl;
      cout << "        Maximum volume: " << MaxVolume << " ft3" << endl;
      cout << "        Relief valve release pressure: " << MaxOverpressure << 
        " lbs/ft2" << endl;
      cout << "        Relief valve coefficient: " << ValveCoefficient << 
        " ft4*sec/slug" << endl;
      cout << "        Initial temperature: " << Temperature << " Rankine" <<
        endl;
      cout << "        Initial pressure: " << Pressure << " lbs/ft2" << endl;
      cout << "        Initial volume: " << Volume << " ft3" << endl;
      cout << "        Initial mass: " << GetMass() << " slug mass" << endl;
      cout << "        Initial weight: " << GetMass()*lbtoslug <<
        " lbs force" << endl;
      cout << "        Heat transfer: " << endl;
    }
  }
  if (debug_lvl & 2 ) { // Instantiation/Destruction notification
    if (from == 0) cout << "Instantiated: FGBallonet" << endl;
    if (from == 1) cout << "Destroyed:    FGBallonet" << endl;
  }
  if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
  }
  if (debug_lvl & 8 ) { // Runtime state variables    
      cout << "        Ballonet holds " << Contents <<
        " mol air" << endl;
      cout << "        Temperature: " << Temperature << " Rankine" << endl;
      cout << "        Pressure: " << Pressure << " lbs/ft2" << endl;
      cout << "        Volume: " << Volume << " ft3" << endl;
      cout << "        Mass: " << GetMass() << " slug mass" << endl;
      cout << "        Weight: " << GetMass()*lbtoslug << " lbs force" << endl;
  }
  if (debug_lvl & 16) { // Sanity checking
  }
  if (debug_lvl & 64) {
    if (from == 0) { // Constructor
      cout << IdSrc << endl;
      cout << IdHdr << endl;
    }
  }
}
}