Merge pull request #872 from su2code/feature_radiation

P1 solver for Radiative Heat Transfer

Author: rsanfer

.gitignore

@@ -76,6 +76,11 @@ bin/
 TestData
 *.stl
 
+# Ignore output files if tests are run locally
+*.vtk
+*.vtu
+*.ref
+
 Mercurial
 .hg*
 

Common/include/CConfig.hpp

@@ -500,6 +500,7 @@ class CConfig {
   Kind_TimeStep_Heat,           /*!< \brief Time stepping method for the (fvm) heat equation. */
   Kind_TimeIntScheme_FEA,       /*!< \brief Time integration for the FEA equations. */
   Kind_SpaceIteScheme_FEA,      /*!< \brief Iterative scheme for nonlinear structural analysis. */
+  Kind_TimeIntScheme_Radiation, /*!< \brief Time integration for the Radiation equations. */
   Kind_ConvNumScheme,           /*!< \brief Global definition of the convective term. */
   Kind_ConvNumScheme_Flow,      /*!< \brief Centered or upwind scheme for the flow equations. */
   Kind_ConvNumScheme_FEM_Flow,  /*!< \brief Finite element scheme for the flow equations. */
@@ -1010,27 +1011,6 @@ class CConfig {
   unsigned short Kind_HybridRANSLES;   /*!< \brief Kind of Hybrid RANS/LES. */
   unsigned short Kind_RoeLowDiss;      /*!< \brief Kind of Roe scheme with low dissipation for unsteady flows. */
   bool QCR;                    /*!< \brief Spalart-Allmaras with Quadratic Constitutive Relation, 2000 version (SA-QCR2000) . */
-  su2double *default_vel_inf,  /*!< \brief Default freestream velocity array for the COption class. */
-  *default_eng_cyl,            /*!< \brief Default engine box array for the COption class. */
-  *default_eng_val,            /*!< \brief Default engine box array values for the COption class. */
-  *default_cfl_adapt,          /*!< \brief Default CFL adapt param array for the COption class. */
-  *default_jst_coeff,          /*!< \brief Default artificial dissipation (flow) array for the COption class. */
-  *default_ffd_coeff,          /*!< \brief Default artificial dissipation (flow) array for the COption class. */
-  *default_mixedout_coeff,     /*!< \brief Default default mixedout algorithm coefficients for the COption class. */
-  *default_rampRotFrame_coeff, /*!< \brief Default ramp rotating frame coefficients for the COption class. */
-  *default_rampOutPres_coeff,  /*!< \brief Default ramp outlet pressure coefficients for the COption class. */
-  *default_jst_adj_coeff,      /*!< \brief Default artificial dissipation (adjoint) array for the COption class. */
-  *default_ad_coeff_heat,      /*!< \brief Default artificial dissipation (heat) array for the COption class. */
-  *default_obj_coeff,          /*!< \brief Default objective array for the COption class. */
-  *default_geo_loc,            /*!< \brief Default SU2_GEO section locations array for the COption class. */
-  *default_distortion,         /*!< \brief Default SU2_GEO section locations array for the COption class. */
-  *default_ea_lim,             /*!< \brief Default equivalent area limit array for the COption class. */
-  *default_grid_fix,           /*!< \brief Default fixed grid (non-deforming region) array for the COption class. */
-  *default_htp_axis,           /*!< \brief Default HTP axis for the COption class. */
-  *default_ffd_axis,           /*!< \brief Default FFD axis for the COption class. */
-  *default_inc_crit,           /*!< \brief Default incremental criteria array for the COption class. */
-  *default_extrarelfac,        /*!< \brief Default extra relaxation factor for Giles BC in the COption class. */
-  *default_sineload_coeff;     /*!< \brief Default values for a sine load. */
 
   unsigned short nSpanWiseSections; /*!< \brief number of span-wise sections */
   unsigned short nSpanMaxAllZones;  /*!< \brief number of maximum span-wise sections for all zones */
@@ -1046,8 +1026,6 @@ class CConfig {
   su2double FinalRotation_Rate_Z;       /*!< \brief Final rotation rate Z if Ramp rotating frame is activated. */
   su2double FinalOutletPressure;        /*!< \brief Final outlet pressure if Ramp outlet pressure is activated. */
   su2double MonitorOutletPressure;      /*!< \brief Monitor outlet pressure if Ramp outlet pressure is activated. */
-  su2double *default_body_force;        /*!< \brief Default body force vector for the COption class. */
-  su2double *default_nacelle_location;  /*!< \brief Location of the nacelle. */
   su2double *default_cp_polycoeffs;     /*!< \brief Array for specific heat polynomial coefficients. */
   su2double *default_mu_polycoeffs;     /*!< \brief Array for viscosity polynomial coefficients. */
   su2double *default_kt_polycoeffs;     /*!< \brief Array for thermal conductivity polynomial coefficients. */
@@ -1070,6 +1048,46 @@ class CConfig {
             *top_optim_filter_radius;  /*!< \brief Radius of the filter(s) used on the design density for topology optimization. */
   unsigned short top_optim_proj_type;  /*!< \brief The projection function used in topology optimization. */
   su2double top_optim_proj_param;      /*!< \brief The value of the parameter for the projection function. */
+  bool HeatSource;              /*!< \brief Flag to know if there is a volumetric heat source on the flow. */
+  su2double ValHeatSource;      /*!< \brief Value of the volumetric heat source on the flow (W/m3). */
+  su2double Heat_Source_Rot_Z;    /*!< \brief Rotation of the volumetric heat source on the Z axis. */
+  su2double *Heat_Source_Center,  /*!< \brief Position of the center of the heat source. */
+            *Heat_Source_Axes;  /*!< \brief Principal axes (x, y, z) of the ellipsoid containing the heat source. */
+  unsigned short Kind_Radiation;       /*!< \brief Kind of radiation model used. */
+  unsigned short Kind_P1_Init;         /*!< \brief Kind of initialization used in the P1 model. */
+  su2double Absorption_Coeff,          /*!< \brief Absorption coefficient of the medium (radiation). */
+  Scattering_Coeff;                    /*!< \brief Scattering coefficient of the medium (radiation). */
+  unsigned short nMarker_Emissivity;   /*!< \brief Number of markers for which the emissivity is defined. */
+  string *Marker_Emissivity;           /*!< \brief Wall markers with defined emissivity. */
+  su2double *Wall_Emissivity;          /*!< \brief Emissivity of the wall. */
+  bool Radiation;                      /*!< \brief Determines if a radiation model is incorporated. */
+  su2double CFL_Rad;                   /*!< \brief CFL Number for the radiation solver. */
+
+  su2double default_vel_inf[3],  /*!< \brief Default freestream velocity array for the COption class. */
+  default_eng_cyl[7],            /*!< \brief Default engine box array for the COption class. */
+  default_eng_val[5],            /*!< \brief Default engine box array values for the COption class. */
+  default_cfl_adapt[4],          /*!< \brief Default CFL adapt param array for the COption class. */
+  default_jst_coeff[2],          /*!< \brief Default artificial dissipation (flow) array for the COption class. */
+  default_ffd_coeff[3],          /*!< \brief Default artificial dissipation (flow) array for the COption class. */
+  default_mixedout_coeff[3],     /*!< \brief Default default mixedout algorithm coefficients for the COption class. */
+  default_rampRotFrame_coeff[3], /*!< \brief Default ramp rotating frame coefficients for the COption class. */
+  default_rampOutPres_coeff[3],  /*!< \brief Default ramp outlet pressure coefficients for the COption class. */
+  default_jst_adj_coeff[2],      /*!< \brief Default artificial dissipation (adjoint) array for the COption class. */
+  default_ad_coeff_heat[2],      /*!< \brief Default artificial dissipation (heat) array for the COption class. */
+  default_obj_coeff[5],          /*!< \brief Default objective array for the COption class. */
+  default_geo_loc[2],            /*!< \brief Default SU2_GEO section locations array for the COption class. */
+  default_distortion[2],         /*!< \brief Default SU2_GEO section locations array for the COption class. */
+  default_ea_lim[3],             /*!< \brief Default equivalent area limit array for the COption class. */
+  default_grid_fix[6],           /*!< \brief Default fixed grid (non-deforming region) array for the COption class. */
+  default_htp_axis[2],           /*!< \brief Default HTP axis for the COption class. */
+  default_ffd_axis[3],           /*!< \brief Default FFD axis for the COption class. */
+  default_inc_crit[3],           /*!< \brief Default incremental criteria array for the COption class. */
+  default_extrarelfac[2],        /*!< \brief Default extra relaxation factor for Giles BC in the COption class. */
+  default_sineload_coeff[3],     /*!< \brief Default values for a sine load. */
+  default_body_force[3],         /*!< \brief Default body force vector for the COption class. */
+  default_nacelle_location[5],   /*!< \brief Location of the nacelle. */
+  default_hs_axes[3],            /*!< \brief Default principal axes (x, y, z) of the ellipsoid containing the heat source. */
+  default_hs_center[3];          /*!< \brief Default position of the center of the heat source. */
 
   unsigned short Riemann_Solver_FEM;         /*!< \brief Riemann solver chosen for the DG method. */
   su2double Quadrature_Factor_Straight;      /*!< \brief Factor applied during quadrature of elements with a constant Jacobian. */
@@ -1094,7 +1112,7 @@ class CConfig {
   Restart_Iter;                  /*!< \brief Determines the restart iteration in the multizone problem */
   su2double Time_Step;           /*!< \brief Determines the time step for the multizone problem */
   su2double Max_Time;            /*!< \brief Determines the maximum time for the time-domain problems */
-  su2double *default_wrt_freq;
+
   unsigned long HistoryWrtFreq[3],    /*!< \brief Array containing history writing frequencies for timer iter, outer iter, inner iter */
                 ScreenWrtFreq[3];     /*!< \brief Array containing screen writing frequencies for timer iter, outer iter, inner iter */
   unsigned long VolumeWrtFreq;        /*!< \brief Writing frequency for solution files. */
@@ -4298,6 +4316,15 @@ class CConfig {
    */
   unsigned short GetKind_TimeIntScheme_FEA(void) const { return Kind_TimeIntScheme_FEA; }
 
+  /*!
+   * \brief Get the kind of integration scheme (explicit or implicit)
+   *        for the radiation equations.
+   * \note This value is obtained from the config file, and it is constant
+   *       during the computation.
+   * \return Kind of integration scheme for the radiation equations.
+   */
+  unsigned short GetKind_TimeIntScheme_Radiation(void) const { return Kind_TimeIntScheme_Radiation; }
+
   /*!
    * \brief Get the kind of integration scheme (explicit or implicit)
    *        for the template equations.
@@ -5837,6 +5864,52 @@ class CConfig {
    */
   const su2double* GetBody_Force_Vector(void) const { return Body_Force_Vector; }
 
+  /*!
+   * \brief Get information about the volumetric heat source.
+   * \return <code>TRUE</code> if it uses a volumetric heat source; otherwise <code>FALSE</code>.
+   */
+  inline bool GetHeatSource(void) const { return HeatSource; }
+
+  /*!
+   * \brief Get information about the volumetric heat source.
+   * \return Value of the volumetric heat source
+   */
+  inline su2double GetHeatSource_Val(void) const {return ValHeatSource;}
+
+  /*!
+   * \brief Get the rotation angle of the volumetric heat source in axis Z.
+   * \return Rotation (Z) of the volumetric heat source
+   */
+  inline su2double GetHeatSource_Rot_Z(void) const {return Heat_Source_Rot_Z;}
+
+  /*!
+   * \brief Set the rotation angle of the volumetric heat source in axis Z.
+   * \param[in] val_rot - Rotation (Z) of the volumetric heat source
+   */
+  inline void SetHeatSource_Rot_Z(su2double val_rot) {Heat_Source_Rot_Z = val_rot;}
+
+  /*!
+   * \brief Get the position of the center of the volumetric heat source.
+   * \return Pointer to the center of the ellipsoid that introduces a volumetric heat source.
+   */
+  inline const su2double* GetHeatSource_Center(void) const {return Heat_Source_Center;}
+
+  /*!
+   * \brief Set the position of the center of the volumetric heat source.
+   * \param[in] x_cent = X position of the center of the volumetric heat source.
+   * \param[in] y_cent = Y position of the center of the volumetric heat source.
+   * \param[in] z_cent = Z position of the center of the volumetric heat source.
+   */
+  inline void SetHeatSource_Center(su2double x_cent, su2double y_cent, su2double z_cent) {
+    Heat_Source_Center[0] = x_cent; Heat_Source_Center[1] = y_cent; Heat_Source_Center[2] = z_cent;
+  }
+
+  /*!
+   * \brief Get the radius of the ellipsoid that introduces a volumetric heat source.
+   * \return Pointer to the radii (x, y, z) of the ellipsoid that introduces a volumetric heat source.
+   */
+  inline const su2double* GetHeatSource_Axes(void) const {return Heat_Source_Axes;}
+
   /*!
    * \brief Get information about the rotational frame.
    * \return <code>TRUE</code> if there is a rotational frame; otherwise <code>FALSE</code>.
@@ -9005,6 +9078,49 @@ class CConfig {
    */
   bool GetSinglezone_Driver(void) const { return SinglezoneDriver; }
 
+  /*!
+   * \brief Get the Kind of Radiation model applied.
+   * \return Kind of radiation model used.
+   */
+  unsigned short GetKind_RadiationModel(void) const { return Kind_Radiation; }
+
+  /*!
+   * \brief Get the Kind of P1 initialization method applied.
+   * \return Kind of P1 initialization method used.
+   */
+  unsigned short GetKind_P1_Init(void) const { return Kind_P1_Init; }
+
+  /*!
+   * \brief Get the value of the absorption coefficient of the medium.
+   * \return Value of the absorption coefficient of the medium.
+   */
+  su2double GetAbsorption_Coeff(void) const { return Absorption_Coeff; }
+
+  /*!
+   * \brief Get the value of the scattering coefficient of the medium.
+   * \return Value of the scattering coefficient of the medium.
+   */
+  su2double GetScattering_Coeff(void) const { return Scattering_Coeff; }
+
+  /*!
+   * \brief Get the wall emissivity at a boundary.
+   * \param[in] val_index - Index corresponding to the boundary.
+   * \return The wall emissivity.
+   */
+  su2double GetWall_Emissivity(string val_index) const;
+
+  /*!
+   * \brief Get the value of the CFL condition for radiation solvers.
+   * \return Value of the CFL condition for radiation solvers.
+   */
+  su2double GetCFL_Rad(void) const { return CFL_Rad; }
+
+  /*!
+   * \brief Determines if radiation needs to be incorporated to the analysis.
+   * \return Radiation boolean
+   */
+  bool AddRadiation(void) const { return Radiation; }
+
   /*!
    * \brief Check if the convergence history of each individual zone is written to screen
    * \return YES if the zone convergence history of each individual zone must be written to screen

Common/include/option_structure.hpp

@@ -105,6 +105,8 @@ const su2double FOUR3 = 4.0 / 3.0;  /*!< \brief Four divided by three. */
 
 const su2double PI_NUMBER = 4.0 * atan(1.0);  /*!< \brief Pi number. */
 
+const su2double STEFAN_BOLTZMANN = 5.670367E-08;  /*!< \brief Stefan-Boltzmann constant in W/(m^2*K^4). */
+
 const int MASTER_NODE = 0;			/*!< \brief Master node for MPI parallelization. */
 const int SINGLE_NODE = 1;			/*!< \brief There is only a node in the MPI parallelization. */
 const int SINGLE_ZONE = 1;			/*!< \brief There is only a zone. */
@@ -436,6 +438,8 @@ enum RUNTIME_TYPE {
   RUNTIME_HEAT_SYS = 21,      /*!< \brief One-physics case, the code is solving the heat equation. */
   RUNTIME_ADJHEAT_SYS = 31,   /*!< \brief One-physics case, the code is solving the adjoint heat equation. */
   RUNTIME_TRANS_SYS = 22,     /*!< \brief One-physics case, the code is solving the turbulence model. */
+  RUNTIME_RADIATION_SYS = 23, /*!< \brief One-physics case, the code is solving the radiation model. */
+  RUNTIME_ADJRAD_SYS = 24,    /*!< \brief One-physics case, the code is solving the adjoint radiation model. */
 };
 
 const int FLOW_SOL = 0;     /*!< \brief Position of the mean flow solution in the solver container array. */
@@ -447,6 +451,8 @@ const int ADJTURB_SOL = 3;  /*!< \brief Position of the continuous adjoint turbu
 const int TRANS_SOL = 4;    /*!< \brief Position of the transition model solution in the solver container array. */
 const int HEAT_SOL = 5;     /*!< \brief Position of the heat equation in the solution solver array. */
 const int ADJHEAT_SOL = 6;  /*!< \brief Position of the adjoint heat equation in the solution solver array. */
+const int RAD_SOL = 7;      /*!< \brief Position of the radiation equation in the solution solver array. */
+const int ADJRAD_SOL = 8;   /*!< \brief Position of the continuous adjoint turbulence solution in the solver container array. */
 
 const int FEA_SOL = 0;      /*!< \brief Position of the FEA equation in the solution solver array. */
 const int ADJFEA_SOL = 1;   /*!< \brief Position of the FEA adjoint equation in the solution solver array. */
@@ -1154,6 +1160,30 @@ static const MapType<string, ENUM_DVFEA> DVFEA_Map = {
   MakePair("ELECTRIC_FIELD", ELECTRIC_FIELD)
 };
 
+/*!
+ * \brief Kinds of radiation models
+ */
+enum ENUM_RADIATION {
+  NO_RADIATION = 0,      /*!< \brief No radiation model */
+  P1_MODEL = 1           /*!< \brief P1 Radiation model. */
+};
+static const MapType<string, ENUM_RADIATION> Radiation_Map = {
+  MakePair("NONE", NO_RADIATION)
+  MakePair("P1", P1_MODEL)
+};
+
+/*!
+ * \brief Kinds of P1 initialization
+ */
+enum ENUM_P1_INIT {
+  P1_INIT_ZERO = 0,      /*!< \brief Initialize the P1 model from zero values */
+  P1_INIT_TEMP = 1       /*!< \brief Initialize the P1 model from blackbody energy computed from the initial temperature. */
+};
+static const MapType<string, ENUM_P1_INIT> P1_Init_Map = {
+  MakePair("ZERO", P1_INIT_ZERO)
+  MakePair("TEMPERATURE_INIT", P1_INIT_TEMP)
+};
+
 /*!
  * \brief Kinds of coupling methods at CHT interfaces.
  * The first (temperature) part determines the BC method on the fluid side, the second (heatflux) part determines