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 BoundaryType
  24 /*! Boundary condition type  */
  25 enum BoundaryTypeTransport      {InletTransport,  OutletTransport, NeumannTransport, DirichletTransport, NoneBCTransport};//Actually Inlet=Dirichlet and Outlet=Neumann
  26
  27 //! enumeration phaseType
  28 /*! The material phase can be Solid, Gas or liquid  */
  29 enum FluidMaterial
  30 {
  31         Air,/**< Material considered is air */
  32         Water/**< Material considered is water */
  33 };
  34
  35 /** \struct LimitField
  36  * \brief value of some fields on the boundary  */
  37 struct LimitFieldTransport{
  38         LimitFieldTransport(){bcType=NoneBCTransport; T=0; h=0; flux=0; }
  39         LimitFieldTransport(BoundaryTypeTransport _bcType, double _T,   double _h,double _flux  ){
  40                 bcType=_bcType; T=_T; h=_h; flux=_flux;
  41         }
  42
  43         BoundaryTypeTransport bcType;
  44         double T; //for inlet or Dirichlet
  45         double h; //for inlet or Dirichlet
  46         double flux; //for Neumann or outlet
  47 };
  48
  49 class TransportEquation: public ProblemCoreFlows
  50 {
  51
  52 public :
  53         /** \fn TransportEquation
  54                          * \brief Constructor for the enthalpy transport in a fluid
  55                          * \param [in] phase : \ref Liquid or \ref Gas
  56                          * \param [in] pressureEstimate : \ref around1bar or \ref around155bars
  57                          * \param [in] vector<double> : fluid velocity (assumed constant)
  58                          *  */
  59         TransportEquation(FluidMaterial fluid, pressureEstimate pEstimate,vector<double> vitesseTransport, MPI_Comm comm = MPI_COMM_WORLD);
  60
  61         //Gestion du calcul
  62         virtual void initialize();
  63         virtual void terminate();//vide la mémoire et enregistre le résultat final
  64         bool initTimeStep(double dt);
  65         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
  66         void abortTimeStep();//efface les inconnues calculées par solveTimeStep() et reinitialise dt à 0
  67         bool iterateTimeStep(bool &ok);
  68         virtual void save();
  69         virtual void validateTimeStep();
  70
  71         /* Boundary conditions */
  72         /** \fn setIntletBoundaryCondition
  73                          * \brief adds a new boundary condition of type Inlet
  74                          * \details same as setDirichletBoundaryCondition
  75                          * \param [in] string : the name of the boundary
  76                          * \param [in] double : the value of the enthalpy at the boundary
  77                          * \param [out] void
  78                          *  */
  79         void setInletBoundaryCondition(string groupName,double enthalpy){
  80                 _limitField[groupName]=LimitFieldTransport(InletTransport,-1,enthalpy,-1);
  81         };
  82         /** \fn setDirichletBoundaryCondition
  83                          * \brief adds a new boundary condition of type Dirichlet
  84                          * \details same as setInletBoundaryCondition
  85                          * \param [in] string : the name of the boundary
  86                          * \param [in] double : the value of the enthalpy at the boundary
  87                          * \param [out] void
  88                          *  */
  89         void setDirichletBoundaryCondition(string groupName,double enthalpy){
  90                 _limitField[groupName]=LimitFieldTransport(DirichletTransport,-1,enthalpy,-1);
  91         };
  92         /** \fn setDirichletBoundaryCondition
  93                          * \brief adds a new boundary condition of type Inlet
  94                          * \details Reads the boundary field in a med file
  95                          * \param [in] string : the name of the boundary
  96                          * \param [in] string : the file name
  97                          * \param [in] string : the field name
  98                          * \param [in] int : the time step number
  99                          * \param [in] int : int corresponding to the enum CELLS or NODES
 100                          * \param [out] void
 101                          *  */
 102         void setDirichletBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
 103         void setDirichletBoundaryCondition(string groupName, Field bc_field){
 104                 _limitField[groupName]=LimitFieldTransport(DirichletTransport, -1, 0, -1);
 105         };
 106
 107         /** \fn setNeumannBoundaryCondition
 108          * \brief adds a new boundary condition of type Neumann
 109          * \details same as setOutletBoundaryCondition
 110          * \param [in] string the name of the boundary
 111          * \param [in] double : the value of the enthalpy flux at the boundary
 112          * \param [out] void
 113          *  */
 114         void setNeumannBoundaryCondition(string groupName, double flux=0){
 115                 _limitField[groupName]=LimitFieldTransport(NeumannTransport,-1,-1,flux);
 116         };
 117         /** \fn setOutletBoundaryCondition
 118          * \brief adds a new boundary condition of type Outlet
 119          * \details same as setNeumannBoundaryCondition
 120          * \param [in] string the name of the boundary
 121          * \param [in] double : the value of the enthalpy flux at the boundary
 122          * \param [out] void
 123          *  */
 124         void setOutletBoundaryCondition(string groupName, double flux=0){
 125                 _limitField[groupName]=LimitFieldTransport(OutletTransport,-1,-1,flux);
 126         };
 127         /** \fn setNeumannBoundaryCondition
 128                          * \brief adds a new boundary condition of type Neumann
 129                          * \details Reads the boundary field in a med file
 130                          * \param [in] string : the name of the boundary
 131                          * \param [in] string : the file name
 132                          * \param [in] string : the field name
 133                          * \param [in] int : the time step number
 134                          * \param [in] int : int corresponding to the enum CELLS or NODES
 135                          * \param [out] void
 136                          *  */
 137         void setNeumannBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
 138         void setNeumannBoundaryCondition(string groupName, Field bc_field){
 139                 _limitField[groupName]=LimitFieldTransport(NeumannTransport,-1,-1, 0);
 140         };
 141
 142         /** \fn setBoundaryFields
 143          * \brief met à jour  _limitField  ( le type de condition limite )
 144          * \details
 145          * \param [in] string
 146          * \param [out] void
 147          *  */
 148         void setBoundaryFields(map<string, LimitFieldTransport> boundaryFields){
 149                 _limitField = boundaryFields;
 150         };
 151
 152
 153         /*Physical parameters*/
 154         void setLiqSatEnthalpy(double hsatl){
 155                 _hsatl=hsatl;
 156         };
 157         void setVapSatEnthalpy(double hsatv){
 158                 _hsatv=hsatv;
 159         };
 160         void setLiqSatDensity(double rhosatl){
 161                 _rhosatl=rhosatl;
 162         };
 163         void setVapSatDensity(double rhosatv){
 164                 _rhosatv=rhosatv;
 165         };
 166         void setTransportVelocity(Vector v){
 167                 _vitesseTransport=v;
 168         };
 169
 170         /* set input fields to prepare the simulation */
 171         vector<string> getInputFieldsNames();
 172         void setInputField(const string& nameField, Field& inputField );//supply of a required input field
 173
 174         /** \fn setRodTemperatureField
 175          * \brief Set the rod temperature field
 176          * \details
 177          * \param [in] Field
 178          * \param [out] void
 179          *  */
 180         void setRodTemperatureField(Field rodTemperature);
 181
 182         /** \fn setRodTemperature
 183          * \brief Set a constant rod temperature field
 184          * \details
 185          * \param [in] double
 186          * \param [out] void
 187          *  */
 188         void setRodTemperature(double rodTemp){
 189                 _rodTemperature=rodTemp;
 190                 _isStationary=false;//Source term may be changed after previously reaching a stationary state
 191         }
 192
 193         /** \fn getRodTemperatureField
 194          * \brief
 195          * \details
 196          * \param [in] void
 197          * \param [out] Field
 198          *  */
 199         Field& getRodTemperatureField(){ // ?? je ne retrouve pas cet attribut dans le file.cxx
 200                 return _rodTemperatureField;
 201         }
 202
 203         /* get output fields for postprocessing or coupling */
 204         vector<string> getOutputFieldsNames() ;//liste tous les champs que peut fournir le code pour le postraitement
 205         Field&         getOutputField(const string& nameField );//Renvoie un champs pour le postraitement
 206
 207         Field& getFluidTemperatureField(){
 208                 return _TT;
 209         }
 210
 211         Field& getEnthalpyField(){
 212                 return _VV;
 213         }
 214
 215         Field& getVoidFractionField(){
 216                 return _Alpha;
 217         }
 218
 219         Field& getDensityField(){
 220                 return _Rho;
 221         }
 222
 223 protected :
 224         double computeTransportMatrix();
 225         double computeRHS();
 226         void updatePrimitives();
 227
 228         /* Postprocessing fields */
 229         Field   _TT, _Alpha, _Rho;//Fields of temperature, void fraction, density. Unknown field is enthalpy (_VV)
 230         double _rhosatv, _rhosatl;
 231         double _Tref, _href, _cpref;
 232
 233         double temperature(double h){
 234                 return _Tref+(h-_href)/_cpref;
 235         };
 236         double voidFraction(double h){
 237                 double titre=(h-_href)/(_hsatv-_hsatl);
 238                 if (titre<0)
 239                         return 0;
 240                 else if (titre>1)
 241                         return 1;
 242                 else return titre*_rhosatl/(titre*_rhosatl+(1-titre)*_rhosatv);
 243         };
 244         double density(double alpha){
 245                 return alpha*_rhosatv+(1-alpha)*_rhosatl;
 246         };
 247
 248         Vector _vitesseTransport, _normale;
 249         bool _transportMatrixSet;
 250         Vec _Hn, _deltaH, _Hk, _Hkm1, _b0;
 251         Vec _Hn_seq; // Local sequential copy of the parallel vector _Hn, used for saving result files
 252         double _dt_transport, _dt_src;
 253
 254         map<string, LimitFieldTransport> _limitField;
 255
 256         /* source terms */
 257         bool   _rodTemperatureFieldSet;
 258         Field  _rodTemperatureField;
 259         double _rodTemperature;
 260 };
 261
 262 #endif /* TransportEquation_HXX_ */