fggascell.cpp

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

    double dU = 0.0;
    for (i = 0; i < HeatTransferCoeff.size(); i++) {
      dU += HeatTransferCoeff[i]->GetValue();
    }
    // Don't include dt when accounting for adiabatic expansion/contraction.
    // The rate of adiabatic cooling looks about right: ~5.4 Rankine/1000ft. 
    if (Contents > 0) {
      Temperature +=
        (dU * dt - Pressure * dVolumeIdeal - BallonetsHeatFlow) /
        (Cv_gas() * Contents * R);
    } else {
      Temperature = AirTemperature;
    }
  } else {
    // No simulation of complex temperature changes.
    // Note: Making the gas cell behave adiabatically might be a better
    // option.
    Temperature = AirTemperature;
  }

  //-- Pressure --
  const double IdealPressure =
    Contents * R * Temperature / (MaxVolume - BallonetsVolume);
  if (IdealPressure > AirPressure + MaxOverpressure) {
    Pressure = AirPressure + MaxOverpressure;
  } else {
    Pressure = max(IdealPressure, AirPressure);
  }

  //-- Manual valving --

  // FIXME: Presently the effect of manual valving is computed using
  //        an ad hoc formula which might not be a good representation
  //        of reality.
  if ((ValveCoefficient > 0.0) && (ValveOpen > 0.0)) {
    // First compute the difference in pressure between the gas in the
    // cell and the air above it.
    // FixMe: CellHeight should depend on current volume.
    const double CellHeight = 2 * Zradius + Zwidth;                   // [ft]
    const double GasMass    = Contents * M_gas();                     // [slug]
    const double GasVolume  = Contents * R * Temperature / Pressure;  // [ft砞
    const double GasDensity = GasMass / GasVolume;
    const double DeltaPressure =
      Pressure + CellHeight * g * (AirDensity - GasDensity) - AirPressure;
    const double VolumeValved =
      ValveOpen * ValveCoefficient * DeltaPressure * dt;
    Contents =
      max(0.0, Contents - Pressure * VolumeValved / (R * Temperature));
  }

  //-- Update ballonets. --
  // Doing that here should give them the opportunity to react to the
  // new pressure.
  BallonetsVolume = 0.0;
  for (i = 0; i < no_ballonets; i++) {
    Ballonet[i]->Calculate(dt);
    BallonetsVolume += Ballonet[i]->GetVolume();
  }

  //-- Automatic safety valving. --
  if (Contents * R * Temperature / (MaxVolume - BallonetsVolume) >
      AirPressure + MaxOverpressure) {
    // Gas is automatically valved. Valving capacity is assumed to be infinite.
    // FIXME: This could/should be replaced by damage to the gas cell envelope.
    Contents =
      (AirPressure + MaxOverpressure) *
      (MaxVolume - BallonetsVolume) / (R * Temperature);
  }

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

  //-- Current buoyancy --
  // The buoyancy is computed using the atmospheres local density.
  Buoyancy = Volume * AirDensity * g;
  
  // Note: This is gross buoyancy. The weight of the gas itself and
  // any ballonets is not deducted here as the effects of the gas mass
  // is handled by FGMassBalance.
  vFn.InitMatrix(0.0, 0.0, - Buoyancy);

  // Compute the inertia of the gas cell.
  // Consider the gas cell as a shape of uniform density.
  // FIXME: If the cell isn't ellipsoid or cylindrical the inertia will
  //        be wrong.
  gasCellJ = FGMatrix33();
  const double mass = Contents * M_gas();
  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.
  gasCellJ(1,1) = Ixx;
  gasCellJ(2,2) = Iyy;
  gasCellJ(3,3) = Izz;
  Mass = mass;
  gasCellM.InitMatrix();
  gasCellM(eX) +=
    GetXYZ(eX) * Mass*slugtolb;
  gasCellM(eY) +=
    GetXYZ(eY) * Mass*slugtolb;
  gasCellM(eZ) +=
    GetXYZ(eZ) * Mass*slugtolb;

  if (no_ballonets > 0) {
    // Add the mass, moment and inertia of any ballonets.
    const FGColumnVector3 p = MassBalance->StructuralToBody( GetXYZ() );

    for (i = 0; i < no_ballonets; i++) {
      Mass += Ballonet[i]->GetMass();
       
      // Add ballonet moments.
      gasCellM(eX) +=
        Ballonet[i]->GetXYZ(eX) * Ballonet[i]->GetMass()*slugtolb;
      gasCellM(eY) +=
        Ballonet[i]->GetXYZ(eY) * Ballonet[i]->GetMass()*slugtolb;
      gasCellM(eZ) +=
        Ballonet[i]->GetXYZ(eZ) * Ballonet[i]->GetMass()*slugtolb;
      
      // Moments of inertia must be converted to the gas cell frame here.
      FGColumnVector3 v =
        MassBalance->StructuralToBody( Ballonet[i]->GetXYZ() ) - p;
      // Body basis is in FT. 
      const double mass = Ballonet[i]->GetMass();
      gasCellJ += Ballonet[i]->GetInertia() +
        FGMatrix33( 0,                - mass*v(1)*v(2), - mass*v(1)*v(3),
                    - mass*v(2)*v(1), 0,                - mass*v(2)*v(3),
                    - mass*v(3)*v(1), - mass*v(3)*v(2), 0 );
    }
  }
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//    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 FGGasCell::Debug(int from)
{
  if (debug_lvl <= 0) return;

  if (debug_lvl & 1) { // Standard console startup message output
    if (from == 0) { // Constructor
      cout << "    Gas cell holds " << Contents << " mol " <<
        type << endl;
      cout << "      Cell 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 << "      Manual 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: FGGasCell" << endl;
    if (from == 1) cout << "Destroyed:    FGGasCell" << endl;
  }
  if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
  }
  if (debug_lvl & 8 ) { // Runtime state variables    
      cout << "      " << type << " cell holds " << Contents << " mol " << 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;
    }
  }
}

/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CLASS IMPLEMENTATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
const double FGBallonet::R = 3.4071;              // [lbs ft/(mol Rankine)]
const double FGBallonet::M_air = 0.0019186;       // [slug/mol]
const double FGBallonet::Cv_air = 5.0/2.0;        // [??]

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

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

  ballonetJ = FGMatrix33();

  MaxVolume = MaxOverpressure = Temperature = Pressure =
    Contents = Volume = dVolumeIdeal = dU = 0.0;
  Xradius = Yradius = Zradius = Xwidth = Ywidth = Zwidth = 0.0;
  ValveCoefficient = ValveOpen = 0.0;
  BlowerInput = NULL;
  CellNum = num;
  Parent = parent;

  // NOTE: In the local system X points north, Y points east and Z points down.
  element = el->FindElement("location");
  if (element) {
    vXYZ = element->FindElementTripletConvertTo("IN");
  } else {
    cerr << "Fatal Error: No location found for this ballonet." << 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.
    // FIXME: 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 ballonet shape." << endl;