CoreFlows/Models/inc/TransportEquation.hxx

   1 //============================================================================
   2 /**
   3  * \file TransportEquation.hxx
   4  * \author Michael NDJINGA
   5  * \version 1.0
   6  * \date 24 March 2015
   7  * \brief Fluid enthalpy transport equation
   8  * */
   9 //============================================================================
  10
  11 /*! \class TransportEquation TransportEquation.hxx "TransportEquation.hxx"
  12  *  \brief Scalar advection equation for a fluid enthalpy
  13  *  \details see \ref TransportEqPage for more details
  14  */
  15 #ifndef TransportEquation_HXX_
  16 #define TransportEquation_HXX_
  17
  18 #include "ProblemCoreFlows.hxx"
  19
  20 using namespace std;
  21
  22
  23 //! enumeration phase
  24 /*! The fluid type can be LiquidPhase or water  */
  25 enum phase
  26 {
  27         LiquidPhase,/**< Fluid considered is GasPhase */
  28         GasPhase/**< Fluid considered is Gas */
  29 };
  30
  31 //! enumeration pressureEstimate
  32 /*! the pressure estimate needed to fit physical parameters  */
  33 enum pressureMagnitude
  34 {
  35         around1bar300KTransport,/**< pressure is around 1 bar and temperature around 300K (for TransportEquation, SinglePhase and IsothermalTwoFluid) or 373 K (saturation for DriftModel and FiveEqsTwoFluid) */
  36         around155bars600KTransport/**< pressure is around 155 bars  and temperature around 618 K (saturation) */
  37 };
  38
  39 //! enumeration BoundaryType
  40 /*! Boundary condition type  */
  41 enum BoundaryTypeTransport      {InletTransport,  OutletTransport, NeumannTransport, DirichletTransport, NoneBCTransport};//Actually Inlet=Dirichlet and Outlet=Neumann
  42
  43 /** \struct LimitField
  44  * \brief value of some fields on the boundary  */
  45 struct LimitFieldTransport{
  46         LimitFieldTransport(){bcType=NoneBCTransport; T=0; h=0; flux=0; }
  47         LimitFieldTransport(BoundaryTypeTransport _bcType, double _T,   double _h,double _flux  ){
  48                 bcType=_bcType; T=_T; h=_h; flux=_flux;
  49         }
  50
  51         BoundaryTypeTransport bcType;
  52         double T; //for inlet or Dirichlet
  53         double h; //for inlet or Dirichlet
  54         double flux; //for Neumann or outlet
  55 };
  56
  57 class TransportEquation: public ProblemCoreFlows
  58 {
  59
  60 public :
  61         /** \fn TransportEquation
  62                          * \brief Constructor for the enthalpy transport in a fluid
  63                          * \param [in] phase : \ref Liquid or \ref Gas
  64                          * \param [in] pressureMagnitude : \ref around1bar or \ref around155bars
  65                          * \param [in] vector<double> : fluid velocity (assumed constant)
  66                          *  */
  67         TransportEquation(phase fluid, pressureMagnitude pEstimate,vector<double> vitesseTransport, MPI_Comm comm = MPI_COMM_WORLD);
  68
  69         //Gestion du calcul
  70         virtual void initialize();
  71         virtual void terminate();//vide la mémoire et enregistre le résultat final
  72         bool initTimeStep(double dt);
  73         double computeTimeStep(bool & stop);//propose un pas de temps pour le calcul. Celà nécessite de discrétiser les opérateur (convection, diffusion, sources) et pour chacun d'employer la condition cfl. En cas de problème durant ce calcul (exemple t=tmax), renvoie stop=true
  74         void abortTimeStep();//efface les inconnues calculées par solveTimeStep() et reinitialise dt à 0
  75         bool iterateTimeStep(bool &ok);
  76         virtual void save();
  77         virtual void validateTimeStep();
  78
  79         /* Boundary conditions */
  80         /** \fn setIntletBoundaryCondition
  81                          * \brief adds a new boundary condition of type Inlet
  82                          * \details
  83                          * \param [in] string : the name of the boundary
  84                          * \param [in] double : the value of the temperature at the boundary
  85                          * \param [out] void
  86                          *  */
  87         void setInletBoundaryCondition(string groupName,double enthalpy){
  88                 _limitField[groupName]=LimitFieldTransport(InletTransport,-1,enthalpy,-1);
  89         };
  90
  91         /** \fn setNeumannBoundaryCondition
  92          * \brief adds a new boundary condition of type Neumann
  93          * \details
  94          * \param [in] string the name of the boundary
  95          * \param [out] void
  96          *  */
  97         void setNeumannBoundaryCondition(string groupName, double flux=0){
  98                 _limitField[groupName]=LimitFieldTransport(NeumannTransport,-1,flux,-1);
  99         };
 100
 101         /** \fn setBoundaryFields
 102          * \brief met à jour  _limitField  ( le type de condition limite )
 103          * \details
 104          * \param [in] string
 105          * \param [out] void
 106          *  */
 107         void setBoundaryFields(map<string, LimitFieldTransport> boundaryFields){
 108                 _limitField = boundaryFields;
 109         };
 110
 111
 112         /*Physical parameters*/
 113         void setLiqSatEnthalpy(double hsatl){
 114                 _hsatl=hsatl;
 115         };
 116         void setVapSatEnthalpy(double hsatv){
 117                 _hsatv=hsatv;
 118         };
 119         void setLiqSatDensity(double rhosatl){
 120                 _rhosatl=rhosatl;
 121         };
 122         void setVapSatDensity(double rhosatv){
 123                 _rhosatv=rhosatv;
 124         };
 125         void setTransportVelocity(Vector v){
 126                 _vitesseTransport=v;
 127         };
 128
 129         /* set input fields to prepare the simulation */
 130         vector<string> getInputFieldsNames();
 131         void setInputField(const string& nameField, Field& inputField );//supply of a required input field
 132
 133         /** \fn setRodTemperatureField
 134          * \brief Set the rod temperature field
 135          * \details
 136          * \param [in] Field
 137          * \param [out] void
 138          *  */
 139         void setRodTemperatureField(Field rodTemperature){
 140                 _rodTemperatureField=rodTemperature;
 141                 _rodTemperatureFieldSet=true;
 142                 _isStationary=false;//Source term may be changed after previously reaching a stationary state
 143         }
 144
 145         /** \fn setRodTemperature
 146          * \brief Set a constant rod temperature field
 147          * \details
 148          * \param [in] double
 149          * \param [out] void
 150          *  */
 151         void setRodTemperature(double rodTemp){
 152                 _rodTemperature=rodTemp;
 153                 _isStationary=false;//Source term may be changed after previously reaching a stationary state
 154         }
 155
 156         /** \fn getRodTemperatureField
 157          * \brief
 158          * \details
 159          * \param [in] void
 160          * \param [out] Field
 161          *  */
 162         Field& getRodTemperatureField(){ // ?? je ne retrouve pas cet attribut dans le file.cxx
 163                 return _rodTemperatureField;
 164         }
 165
 166         /* get output fields for postprocessing or coupling */
 167         vector<string> getOutputFieldsNames() ;//liste tous les champs que peut fournir le code pour le postraitement
 168         Field&         getOutputField(const string& nameField );//Renvoie un champs pour le postraitement
 169
 170         Field& getFluidTemperatureField(){
 171                 return _TT;
 172         }
 173
 174         Field& getEnthalpyField(){
 175                 return _VV;
 176         }
 177
 178         Field& getVoidFractionField(){
 179                 return _Alpha;
 180         }
 181
 182         Field& getDensityField(){
 183                 return _Rho;
 184         }
 185
 186 protected :
 187         double computeTransportMatrix();
 188         double computeRHS();
 189         void updatePrimitives();
 190
 191         /* Postprocessing fields */
 192         Field   _TT, _Alpha, _Rho;//Fields of temperature, void fraction, density. Unknown field is enthalpy (_VV)
 193         double _rhosatv, _rhosatl;
 194         double _Tref, _href, _cpref;
 195
 196         double temperature(double h){
 197                 return _Tref+(h-_href)/_cpref;
 198         };
 199         double voidFraction(double h){
 200                 double titre=(h-_href)/(_hsatv-_hsatl);
 201                 if (titre<0)
 202                         return 0;
 203                 else if (titre>1)
 204                         return 1;
 205                 else return titre*_rhosatl/(titre*_rhosatl+(1-titre)*_rhosatv);
 206         };
 207         double density(double alpha){
 208                 return alpha*_rhosatv+(1-alpha)*_rhosatl;
 209         };
 210
 211         Vector _vitesseTransport, _normale;
 212         bool _transportMatrixSet;
 213         Vec _Hn, _deltaH, _Hk, _Hkm1, _b0;
 214         Vec _Hn_seq; // Local sequential copy of the parallel vector _Hn, used for saving result files
 215         double _dt_transport, _dt_src;
 216
 217         map<string, LimitFieldTransport> _limitField;
 218
 219         /* source terms */
 220         bool   _rodTemperatureFieldSet;
 221         Field  _rodTemperatureField;
 222         double _rodTemperature;
 223 };
 224
 225 #endif /* TransportEquation_HXX_ */