C++ input::StagInletDy方法代码示例

//.........这里部分代码省略.........
      ghostState.data[1] = (*this).U() + normArea.X() * ( (*this).P() - ghostState.P() ) / rhoSoSInt;
      ghostState.data[2] = (*this).V() + normArea.Y() * ( (*this).P() - ghostState.P() ) / rhoSoSInt;
      ghostState.data[3] = (*this).W() + normArea.Z() * ( (*this).P() - ghostState.P() ) / rhoSoSInt;
    }
    else {
      cerr << "ERROR: flow condition for characteristic BC is not recognized!" << endl;
      cerr << "Interior state: " << (*this) << endl;
      cerr << "Ghost state: " << ghostState << endl;
      exit(0);
    }

    if (layer == 2){ //extrapolate to get ghost state at 2nd layer
      ghostState = 2.0 * ghostState - (*this);
    }

  }
  //supersonic inflow boundary condition ------------------------------------------------------------------------------------------------------------
  //this boundary condition enforces the entire state as the specified freestream state
  else if (bcType == "supersonicInflow"){
    //physical boundary conditions - fix everything
    double sos = eqnState.GetSoS(inputVars.PRef(),inputVars.RRef());
    vector3d<double> vel = inputVars.VelRef() / sos;                       //nondimensional velocity
    ghostState.data[0] = 1.0; //nondimensional density
    ghostState.data[1] = vel.X();
    ghostState.data[2] = vel.Y();
    ghostState.data[3] = vel.Z();
    ghostState.data[4] = 1.0 / eqnState.Gamma(); //nondimensional pressure
  }
  //supersonic outflow boundary condition ------------------------------------------------------------------------------------------------------------
  //this boundary condition enforces the entire state as extrapolated from the interior (zeroth order extrapolation)
  else if (bcType == "supersonicOutflow"){
    //do nothing and return boundary state -- numerical BCs for all
    if (layer == 2){ //extrapolate to get ghost state at 2nd layer
      ghostState = 2.0 * ghostState - (*this);
    }
  }
  //stagnation inlet boundary condition --------------------------------------------------------------------------------------------------------------
  //this boundary condition is appropriate for subsonic flow. It is particularly well suited for internal flows. Implementation from Blazek
  else if (bcType == "stagnationInlet"){
    double g = eqnState.Gamma() - 1.0;
    //calculate outgoing riemann invarient
    double rNeg = (*this).Velocity().DotProd(normArea) - 2.0 * (*this).SoS(eqnState) / g;

    //calculate SoS on boundary
    double cosTheta = -1.0 * (*this).Velocity().DotProd(normArea) / (*this).Velocity().Mag();
    double stagSoSsq = pow((*this).SoS(eqnState),2.0) + 0.5 * g * (*this).Velocity().MagSq();

    double sosB = -1.0 * rNeg * g / (g * cosTheta * cosTheta + 2.0) * (1.0 + cosTheta * sqrt( (g * cosTheta * cosTheta + 2.0) * stagSoSsq / 
											      (g * rNeg * rNeg) - 0.5 * g  ) );
    double tb = inputVars.StagInletT0() / inputVars.TRef() * (sosB * sosB / stagSoSsq);
    double aRef = eqnState.GetSoS(inputVars.PRef(),inputVars.RRef());
    double pb = inputVars.StagInletP0() / (inputVars.RRef() * aRef * aRef) * pow(sosB * sosB / stagSoSsq, eqnState.Gamma()/g);
    double vbMag = sqrt(2.0 / g * (inputVars.StagInletT0() / inputVars.TRef() - tb));

    ghostState.data[0] = eqnState.GetDensityTP(tb,pb);
    ghostState.data[1] = vbMag * inputVars.StagInletDx();
    ghostState.data[2] = vbMag * inputVars.StagInletDy();
    ghostState.data[3] = vbMag * inputVars.StagInletDz();
    ghostState.data[4] = pb;

    if (layer == 2){ //extrapolate to get ghost state at 2nd layer
      ghostState = 2.0 * ghostState - (*this);
    }

  }
  //pressure outlet boundary condition ---------------------------------------------------------------------------------------------------------------
  //this boundary condition is appropriate for subsonic flow. Implementation from Blazek
  else if (bcType == "pressureOutlet"){
    double aRef = eqnState.GetSoS(inputVars.PRef(),inputVars.RRef()); //reference speed of sound
    double pb = inputVars.PressureOutletP() / (inputVars.RRef() * aRef * aRef); //nondimensional pressure from input file

    double SoSInt = (*this).SoS(eqnState);
    double rhoSoSInt = (*this).Rho() * SoSInt;
    ghostState.data[4] = pb;
    ghostState.data[0] = (*this).Rho() + ( ghostState.P() - (*this).P() ) / (SoSInt * SoSInt);
    ghostState.data[1] = (*this).U() + normArea.X() * ( (*this).P() - ghostState.P() ) / rhoSoSInt;
    ghostState.data[2] =  (*this).V() + normArea.Y() * ( (*this).P() - ghostState.P() ) / rhoSoSInt;
    ghostState.data[3] =  (*this).W() + normArea.Z() * ( (*this).P() - ghostState.P() ) / rhoSoSInt;

    if (layer == 2){ //extrapolate to get ghost state at 2nd layer
      ghostState = 2.0 * ghostState - (*this);
    }

  }
  //interblock boundary condition ---------------------------------------------------------------------------------------------------------------
  //this boundary condition is appropriate for point matched interfaces between physical blocks or processor blocks
  else if (bcType == "interblock"){
    //do nothing -- assign interior state to ghost state (already done)
    //for second layer of ghost cells interior state should be 2nd interior cell

  }
  else {
    cerr << "ERROR: Error in primVars::GetGhostState ghost state for BC type " << bcType << " is not supported!" << endl;
    cerr << "surface is " << surf << endl;
    exit(0);
  }

  return ghostState;

}